武禎， 路偉， 鄢書暢， 邱睿， 張輝， 李君利

(1.清華大學(xué) 工程物理系， 北京 100084； 2.同方威視技術(shù)股份有限公司， 北京 100084； 3.中國(guó)人民解放軍疾病預(yù)防控制中心核化防護(hù)科， 北京 100071)

輻射劑量計(jì)算是研究輻射效應(yīng)的基礎(chǔ)[1-2]，人體輻射劑量計(jì)算廣泛應(yīng)用于醫(yī)學(xué)放射診斷與治療、事故應(yīng)急照射、職業(yè)照射和環(huán)境所致公眾內(nèi)、外照射等領(lǐng)域。常見人體輻射劑量計(jì)算方法包含4類：確定論方法[2]、解析方法[3]、通用蒙卡方法[4-5]和快速蒙卡方法[6-9]。不同應(yīng)用領(lǐng)域?qū)τ谳椛鋭┝坑?jì)算在計(jì)算精確度和計(jì)算時(shí)間都有一定的要求。解析方法、確定論方法及快速蒙卡方法計(jì)算速度較快，但是在物理模型、幾何、材料方面都做了一定的近似，因此在計(jì)算精確度方面都有一定的妥協(xié)。同時(shí)，作為劑量計(jì)算的“金標(biāo)準(zhǔn)”，蒙卡方法能夠模擬得到問(wèn)題的精確解，但計(jì)算時(shí)間明顯過(guò)長(zhǎng)，這樣又限制了蒙卡程序在對(duì)計(jì)算時(shí)間有一定的要求領(lǐng)域的應(yīng)用，比如臨床放射治療等。

近些年，隨著計(jì)算機(jī)硬件的不斷發(fā)展，蒙卡程序模擬的計(jì)算時(shí)間逐漸減少。依據(jù)摩爾定律，CPU性能每18個(gè)月提高一倍，提高性能的2種主要途徑:1)集成更多的晶體管、提高時(shí)鐘頻率；2)傳統(tǒng)的串行程序加速模式[10]。但每次硬件更新對(duì)程序?qū)嶋H性能的提升依然無(wú)法滿足需求，同時(shí)考慮到能耗、計(jì)算成本等問(wèn)題，使用CPU運(yùn)行蒙卡程序遇到瓶頸。近年來(lái)GPU發(fā)展迅速，與CPU相比，GPU具備了更強(qiáng)的單精度浮點(diǎn)運(yùn)算能力和內(nèi)存帶寬，且GPU硬件系統(tǒng)更易維護(hù)而且成本較低；同時(shí)，雖然2010年以后，國(guó)外依托原快速蒙卡程序或從零開始編寫了一些基于GPU加速的光子輸運(yùn)或光電耦合輸運(yùn)程序，包括ACHEERrt、GPUMCD、goMC和gDPM，并進(jìn)行了大量的正確性驗(yàn)證和加速效率測(cè)試[11-13]，但進(jìn)一步調(diào)研上述文獻(xiàn)可以發(fā)現(xiàn)，光子電子物理過(guò)程做了不同近似處理，因此在一些需要精細(xì)模擬的應(yīng)用中存在一些缺陷。綜合考慮上述情況，本文使用完整的光子電子精細(xì)物理模型開發(fā)了GPU加速光電耦合輸運(yùn)蒙卡程序(GPU-based photon-electron coupled accelerated dose estimation program, Gadep)，并進(jìn)行了多方面的性能優(yōu)化。

1 程序框架設(shè)計(jì)及算法實(shí)現(xiàn)

1.1 程序架構(gòu)設(shè)計(jì)

本文開發(fā)的基于GPU加速的光電耦合輸運(yùn)蒙卡程序Gadep主要用于進(jìn)行基于數(shù)字化人體體素模型的輻射劑量蒙卡模擬計(jì)算。程序設(shè)計(jì)框架分為4層，如圖1所示。第1層為硬件環(huán)境層，主要包含CPU和GPU，作為程序的底層硬件支撐；第2層為軟件環(huán)境層，主要包括開發(fā)環(huán)境、開發(fā)語(yǔ)言及相關(guān)運(yùn)行庫(kù)等，其中：開發(fā)環(huán)境為Visual Studio 2015和CUDA 10.0，開發(fā)語(yǔ)言包括C++、C、CUDA C，涉及到運(yùn)行庫(kù)包括CURAND隨機(jī)數(shù)庫(kù)、Thrust庫(kù)以及支撐CUDA程序運(yùn)行所必須的運(yùn)行庫(kù)和顯卡驅(qū)動(dòng)程序；第3層為功能模塊層，包括程序的主要功能模塊，如幾何模塊、源項(xiàng)模塊、輸運(yùn)模塊、物理模塊及統(tǒng)計(jì)模塊等；第4層為用戶層，包括蒙卡模擬參數(shù)的輸入及結(jié)果的輸出及展示等，其中用戶層集成已經(jīng)在本研究團(tuán)隊(duì)開發(fā)三維輻射劑量計(jì)算與防護(hù)設(shè)計(jì)仿真程序THUDose中，并進(jìn)行了大量的正確性驗(yàn)證和加速效率測(cè)試[14]。

圖1 基于GPU加速的光電耦合輸運(yùn)蒙卡程序Gadep設(shè)計(jì)框架Fig.1 Design framework of GPU-based photon-electron coupled accelerated dose estimation program (Gadep)

1.2 程序物理模型選擇

程序主要包含的物理過(guò)程為光子電子耦合輸運(yùn)過(guò)程。光子物理過(guò)程本質(zhì)上是光子與原子核和核外電子發(fā)生反應(yīng)的過(guò)程，前者包含電子對(duì)效應(yīng)和光核反應(yīng)，后者包含彈性散射、光電效應(yīng)和康普頓散射，其中光電效應(yīng)、康普頓散射和電子對(duì)效應(yīng)3種過(guò)程中會(huì)分別產(chǎn)生光電子、散射電子和反沖正負(fù)電子對(duì)；電子物理過(guò)程包含軔致輻射、電離作用、多重散射和正電子湮滅，其中軔致輻射和正電子湮滅可產(chǎn)生次級(jí)光子。全部光電耦合輸運(yùn)過(guò)程如圖2所示。

圖2 光電耦合輸運(yùn)過(guò)程示意Fig.2 Schematic diagram of coupled photon-electron transport process

考慮到本文實(shí)現(xiàn)的光電耦合輸運(yùn)蒙卡程序的適用范圍，為程序選擇了合適的物理模型，軔致輻射和電離輻射包含了連續(xù)和離散過(guò)程，正電子湮滅包含了靜止和離散過(guò)程，多重散射則為連續(xù)過(guò)程。各個(gè)物理過(guò)程對(duì)應(yīng)的物理模型、適用能量范圍如表1所示。其中：光子物理過(guò)程主要參考了Geant4提供的Livermore低能模型[5]，與Geant4標(biāo)準(zhǔn)物理模型在低能部分的截面主要依靠參數(shù)化計(jì)算相比，Livermore低能模型直接利用核素的殼層反應(yīng)截面信息，結(jié)果更為準(zhǔn)確，可模擬光子的能量也更低，對(duì)于Z在1～99的核素范圍能量模擬下限為100 eV；多重散射模型則主要使用了基于Lewis理論的多重散射Urban物理模型，該模型的優(yōu)點(diǎn)是適用于包含電子在內(nèi)的各種帶電粒子，缺點(diǎn)是其在低能區(qū)忽略了碰撞核大小因素，電子模擬的精確度不如Penelope、EGSnrc等采用的Goudsmit-Saunderson模型[15]。

表1 光子電子物理模型Table 1 Photon and electron physical model

1.3 程序?qū)崿F(xiàn)

基于上述程序設(shè)計(jì)架構(gòu)及所選擇的光電耦合物理過(guò)程，程序?qū)崿F(xiàn)流程如圖3所示，具體描述如下：首先，在GPU上預(yù)分配內(nèi)存空間，通過(guò)CUDA系統(tǒng)數(shù)據(jù)傳輸函數(shù)，將在CPU程序初始化生成的源項(xiàng)、截面、常量、隨機(jī)數(shù)種子、幾何、統(tǒng)計(jì)和材料等數(shù)據(jù)傳輸?shù)筋A(yù)分配的內(nèi)存上；然后，在顯卡設(shè)備端(device)實(shí)現(xiàn)光子、正電子和負(fù)電子3個(gè)粒子輸運(yùn)內(nèi)核函數(shù)，可根據(jù)入射粒子或次級(jí)粒子數(shù)組粒子類型調(diào)用相應(yīng)的輸運(yùn)函數(shù)；上述內(nèi)核函數(shù)實(shí)現(xiàn)后，在CPU主機(jī)端(host)選取合適的設(shè)備內(nèi)核函數(shù)配置(即執(zhí)行的線程塊、線程網(wǎng)格的個(gè)數(shù))，并預(yù)先給每個(gè)線程分配若干(一般取50)次級(jí)粒子數(shù)組，然后根據(jù)初始入射粒子及其次級(jí)粒子類型調(diào)用輸運(yùn)內(nèi)核函數(shù)，按照先進(jìn)后出的原則進(jìn)行模擬。

圖3 Gadep具體實(shí)現(xiàn)流程Fig.3 Specific implementation process of Gadep

1.4 步長(zhǎng)計(jì)算方法選擇

在通用蒙卡程序中，粒子的輸運(yùn)是以步(step)循環(huán)進(jìn)行模擬，因此步長(zhǎng)計(jì)算是粒子輸運(yùn)的關(guān)鍵環(huán)節(jié)之一。通常步長(zhǎng)大小計(jì)算公式為：

s=nλλ

(1)

式中：λ為粒子的平均自由程；nλ為粒子在介質(zhì)中輸運(yùn)的平均自由程的個(gè)數(shù)，其值計(jì)算為：

nλ=-logη

(2)

式中η為(0，1)上的隨機(jī)數(shù)。

(3)

為了比較2種步長(zhǎng)計(jì)算算法的優(yōu)劣，假設(shè)2種算法下都經(jīng)過(guò)了n步，穿越了m次幾何邊界，同時(shí)假設(shè)每次獲取分步長(zhǎng)平均自由程時(shí)間都是Tλ、生成單個(gè)隨機(jī)數(shù)時(shí)間是Tδ，則總截面法和分截面法時(shí)間總計(jì)分別為n[3Tλ+2Tδ]+m[3Tλ+2Tδ]和[n[3Tλ+Tδ]+2Tδ]+m×3Tλ?？梢钥吹?，2種算法在平均自由程數(shù)據(jù)訪問(wèn)次數(shù)上一致，但是在生成隨機(jī)數(shù)次數(shù)上，同一幾何體內(nèi)總截面法和分截面法之比為2n/(n+2)，穿越幾何邊界時(shí)總截面法和分截面法之比為1+2Tδ/3Tλ。因此在使用數(shù)字體素模型進(jìn)行劑量計(jì)算的實(shí)際應(yīng)用中，人體模型中的體素個(gè)數(shù)通常在百萬(wàn)量級(jí)，模擬的粒子數(shù)在千萬(wàn)量級(jí)，采用分截面法需要的隨機(jī)數(shù)將大大減少，從而可以進(jìn)一步減少時(shí)間開銷，提高計(jì)算速度。

1.5 基于顯卡內(nèi)存的數(shù)據(jù)存儲(chǔ)方案設(shè)計(jì)

程序中用到的數(shù)據(jù)類型按用途分為基本物理常量、程序輸入輸出量和初始化數(shù)據(jù)?；疚锢砹堪锢砘締挝?、物理常量、物理模型參數(shù)、原子結(jié)合能數(shù)據(jù)表、原子殼層表等。程序輸入輸出量包含入射源項(xiàng)、幾何數(shù)據(jù)(體素器官ID、器官材料號(hào)、體素密度)。初始化數(shù)據(jù)主要是指原子反應(yīng)截面數(shù)據(jù)，根據(jù)輸入材料信息生成的康普頓散射、電子對(duì)效應(yīng)和光電效應(yīng)產(chǎn)生的平均自由程數(shù)據(jù)庫(kù)，以及碰撞原子核抽樣截面庫(kù)等。

顯卡內(nèi)存分全局內(nèi)存(global)、共享內(nèi)存(shared)、常量?jī)?nèi)存(constant)和紋理內(nèi)存(texture)幾種，其大小、讀寫速度、位置上具有不同的特性，因此也有不同的適用性。通過(guò)分析程序現(xiàn)有數(shù)據(jù)，并結(jié)合顯卡內(nèi)存各種類特點(diǎn)，本文將程序中數(shù)據(jù)類型設(shè)計(jì)為常量、一維數(shù)組、枚舉和結(jié)構(gòu)體，并放置于不同的顯卡內(nèi)存中。

此外，由于設(shè)備端的自定義函數(shù)不能直接訪問(wèn)主機(jī)端內(nèi)存，因此，需要利用CUDA內(nèi)置函數(shù)將主機(jī)端的數(shù)據(jù)拷貝到設(shè)備端才能訪問(wèn)。程序中使用的CUDA內(nèi)置內(nèi)存拷貝函數(shù)拷貝全局內(nèi)存變量數(shù)據(jù)和常量?jī)?nèi)存變量數(shù)據(jù)，一般初始化后在主機(jī)端調(diào)用這些函數(shù)就可以將變量值、數(shù)組等簡(jiǎn)單數(shù)據(jù)結(jié)構(gòu)拷貝到預(yù)先分配的設(shè)備端變量中。但是，對(duì)于自定義包含指針的結(jié)構(gòu)體類型，CUDA并沒有內(nèi)置結(jié)構(gòu)體內(nèi)存拷貝函數(shù)，因此還需要在程序中編寫相應(yīng)的結(jié)構(gòu)體拷貝函數(shù)實(shí)現(xiàn)結(jié)構(gòu)體的內(nèi)存數(shù)據(jù)轉(zhuǎn)移到GPU內(nèi)存功能。

2 正確性檢驗(yàn)

為了驗(yàn)證本文開發(fā)的GPU加速蒙卡程序Gadep的正確性，選取ICRP 116號(hào)報(bào)告[16]中人體器官和組織受照劑量計(jì)算6種照射情景之中的前后照射(antero-posterior，AP)和后前照射(posterior-antero，PA)照射作為模擬對(duì)象，如圖4所示，統(tǒng)計(jì)人體內(nèi)光電耦合輸運(yùn)下的能量沉積，計(jì)算人體器官或組織劑量轉(zhuǎn)換系數(shù)，并和ICRP 116號(hào)報(bào)告結(jié)果進(jìn)行對(duì)比。ICRP 116號(hào)報(bào)告結(jié)果是程序FLUKA、MCNPX、PHITS、Geant4、EGSnrc等計(jì)算結(jié)果經(jīng)過(guò)平均、平滑和擬合處理后得到的。

圖4 AP和PA照射示意Fig.4 AP and PA irradiation

使用本程序分別模擬了AP、PA照射下入射粒子分別為光子、電子4種情形下人體器官或組織的外照射器官劑量轉(zhuǎn)換系數(shù)，采用的體素模型為ICRP成年男性參考人體模(AM)，模擬粒子數(shù)為1×108，所使用顯卡為NVIDIA GeForce 1080Ti，除深部、小器官外其他器官或組織的外照射劑量轉(zhuǎn)換系數(shù)不確定度小于1.0%。圖5～8分別為一些代表性器官光子、電子外照射劑量轉(zhuǎn)換系數(shù)計(jì)算結(jié)果和ICRP 116號(hào)報(bào)告中的外照射劑量轉(zhuǎn)換系數(shù)對(duì)比圖。

從比對(duì)結(jié)果可以看出，使用Gadep計(jì)算得到的外照射器官劑量轉(zhuǎn)換系數(shù)與ICRP 116號(hào)報(bào)告結(jié)果基本一致，驗(yàn)證了程序的正確性。對(duì)于入射粒子為光子的AP和PA照射(圖5和圖6)，與有效劑量相關(guān)的器官劑量轉(zhuǎn)換系數(shù)計(jì)算結(jié)果和ICRP 116號(hào)報(bào)告[16]之間的差異多數(shù)可以達(dá)到1%以內(nèi)，部分小器官或深部器官，由于粒子進(jìn)入的概率低，不確定度相對(duì)較大，但與ICRP116號(hào)報(bào)告的差異也基本在3%以內(nèi)。

對(duì)于入射粒子為電子的AP和PA照射(圖7和圖8)，在低能部分Gadep和ICRP 116號(hào)報(bào)告計(jì)算結(jié)果之間差異較大。這是由于對(duì)于同樣的粒子數(shù)，入射粒子為電子時(shí)深部器官計(jì)算結(jié)果不確定度會(huì)比光子要大很多。ICRP 116號(hào)報(bào)告中也指出，600 keV的入射電子沿AP照射情況下，大器官和小器官的相對(duì)統(tǒng)計(jì)不確定度分別為4%和10%。

圖5 光子AP外照射器官劑量轉(zhuǎn)換系數(shù)計(jì)算結(jié)果和ICRP 116號(hào)報(bào)告結(jié)果對(duì)比Fig.5 Comparison of calculation results of dose conversion coefficient of external photon AP irradiation by Gadep program with the results reported by ICRP 116

圖6 光子PA照射器官劑量轉(zhuǎn)換系數(shù)計(jì)算結(jié)果和ICRP 116號(hào)報(bào)告結(jié)果對(duì)比Fig.6 Comparison of calculation results of dose conversion coefficient of external photon PA irradiation by Gadep program with the results reported by ICRP 116

圖7 電子AP外照射器官劑量轉(zhuǎn)換系數(shù)計(jì)算結(jié)果和ICRP 116號(hào)報(bào)告結(jié)果對(duì)比Fig.7 Comparison of calculation results of dose conversion coefficient of external electron AP irradiation by Gadep program with the results reported by ICRP 116

圖8 電子PA外照射器官劑量轉(zhuǎn)換系數(shù)計(jì)算結(jié)果和ICRP 116號(hào)報(bào)告結(jié)果對(duì)比Fig.8 Comparison of calculation results of dose conversion coefficient of external electron PA irradiation by Gadep program with the results reported by ICRP 116

3 程序優(yōu)化及加速效果比對(duì)

3.1 幾何輸運(yùn)算法優(yōu)化

跨邊界算法是一種根據(jù)體素模型所做的幾何輸運(yùn)算法的優(yōu)化，如圖9所示。其原理是相鄰體素如果材料相同，則光子可以穿越幾何邊界而不停頓，減少輸運(yùn)過(guò)程中對(duì)同種材料進(jìn)入相鄰體素后反應(yīng)截面重新抽樣次數(shù)，達(dá)到減少計(jì)算時(shí)間、提高加速效率的目的。

圖9 非跨幾何和跨幾何邊界輸運(yùn)算法示意Fig.9 Schematic diagram of the cross and non-cross geometric boundary transport algorithms

下面以AP照射情形下全身體素徑跡長(zhǎng)度分布計(jì)算為例，對(duì)使用跨邊界算法的正確性以及優(yōu)化效果進(jìn)行了驗(yàn)證。表2給出了不同能量下是否使用跨邊界算法的情況下每個(gè)體素的計(jì)算差異及優(yōu)化效果，其中σ是采用跨幾何邊界和非跨幾何邊界計(jì)算結(jié)果相對(duì)差異，所對(duì)應(yīng)的數(shù)值為相對(duì)差異小于σ的體素百分比，可以看出，在不同的能量點(diǎn)上，約80%的體素計(jì)算結(jié)果差異在1%以內(nèi)；表格最后一行給出了是否使用跨邊界算法的優(yōu)化效果，可以看出計(jì)算時(shí)間減少了29.4%～41.8%。

表2 跨幾何邊界算法下每個(gè)體素計(jì)算結(jié)果差異及優(yōu)化效果

3.2 GPU硬件結(jié)構(gòu)優(yōu)化

GPU計(jì)算性能一般用浮點(diǎn)數(shù)運(yùn)算峰值(FLOPS)來(lái)表示，一般GPU的單精度FLOPS要比雙精度FLOPS要高。因此在程序計(jì)算精度允許的情況下，使用單精度浮點(diǎn)類型可以提高GPU程序的計(jì)算效率。同樣以AP照射情形下全身體素徑跡長(zhǎng)度分布計(jì)算為例，測(cè)試了將計(jì)算精確度要求較高的長(zhǎng)度相關(guān)物理量分別設(shè)為單、雙精度浮點(diǎn)類型，對(duì)計(jì)算結(jié)果和計(jì)算時(shí)間的影響，見表3。可以看出，60%左右的體素徑跡長(zhǎng)度計(jì)算結(jié)果相對(duì)差異小于1%，相應(yīng)的計(jì)算時(shí)間減少了44.4%～58.1%。

3.3 GPU內(nèi)核函數(shù)配置優(yōu)化

CUDA中線程執(zhí)行的最小單元是線程束(warp)，每個(gè)線程束中包含32個(gè)線程，每個(gè)線程塊包含多個(gè)線程束。GPU程序在執(zhí)行內(nèi)核時(shí)，需要分配線程塊的大小，線程塊的大小和GPU硬件架構(gòu)相關(guān)，通常取32的整數(shù)倍；同時(shí)，CUDA內(nèi)核所占資源(寄存器、共享內(nèi)存、紋理內(nèi)存等)也會(huì)限制GPU線程塊池(pool)可執(zhí)行線程塊的個(gè)數(shù)。上述2個(gè)方面的限制都對(duì)GPU程序加速效率有影響，圖10給出了采用不同的線程執(zhí)行模型進(jìn)行AP照射情形下全身體素徑跡長(zhǎng)度分布計(jì)算時(shí)，不同能量的光子輸運(yùn)運(yùn)行時(shí)間比較。由圖可知，當(dāng)配置參數(shù)取32時(shí)，各能量點(diǎn)計(jì)算時(shí)間最少，32為程序Gadep在顯卡GeForce 1080Ti的最優(yōu)內(nèi)核函數(shù)配置參數(shù)。

表3 單精度和雙精度浮點(diǎn)類型Gadep的計(jì)算結(jié)果差異及優(yōu)化效果

圖10 不同內(nèi)核執(zhí)行模型配置參數(shù)下計(jì)算時(shí)間對(duì)比Fig.10 Comparison of calculation time under different kernel configuration parameters

3.4 程序加速效果比對(duì)

為對(duì)程序Gadep的加速效果進(jìn)行驗(yàn)證，對(duì)該程序計(jì)算AP、PA標(biāo)準(zhǔn)照射條件下光子、電子外照射劑量轉(zhuǎn)換系數(shù)的計(jì)算時(shí)間進(jìn)行統(tǒng)計(jì)，并與同等計(jì)算條件下使用MCNP進(jìn)行計(jì)算的計(jì)算時(shí)間進(jìn)行了比對(duì)。其中MCNP中取粒子模擬類型mode P E,為光子電子耦合輸運(yùn)，光子、電子截?cái)嗄芰繛槌绦蚰J(rèn)值；測(cè)試CPU為Intel i7-4710U，GPU顯卡為NVIDIA GeForce 1080Ti。圖11給出了4種計(jì)算情景下Gadep程序與MCNP程序的加速效率對(duì)比結(jié)果。

從圖中可以看出，對(duì)于光子照射情景，Gadep在低能區(qū)加速效率可以達(dá)到224倍，中能區(qū)加速效率在50～100倍；當(dāng)入射光子能量增加，次級(jí)粒子數(shù)多，線程歧離變大，加速效率值總體逐漸降低，當(dāng)光子能量為100 MeV時(shí)，加速效率值約為25倍。對(duì)于電子照射情景，Gadep在電子能量為0.6 MeV時(shí)加速效率最大值達(dá)到300倍，增加或降低電子能量，加速效率值變?。籄P和PA照射下最小加速效率值分別為48倍和53倍，對(duì)應(yīng)電子能量分別為0.01 MeV和100 MeV。

圖11 程序Gadep加速效率隨入射粒子能量、類型和照射情景變化趨勢(shì)Fig.11 Trend chart of Gadep program acceleration effect over incident particle energy, type and irradiation situation

4 程序應(yīng)用

對(duì)事故進(jìn)行劑量重建是評(píng)價(jià)放射源造成輻照危害大小的重要手段。作者曾使用THUDose和中國(guó)參考人體素模型庫(kù)對(duì)2014年5月南京工業(yè)探傷機(jī)放射源192Ir丟失所致人員受照事故進(jìn)行了精細(xì)的物理劑量重建，模擬了病人全身光子通量、全身器官吸收劑量分布情況，并統(tǒng)計(jì)了計(jì)算時(shí)間[17-18]。本文使用了Gadep程序，采用相同的模擬條件，對(duì)南京放射源丟失事故中的病人器官吸收劑量進(jìn)行計(jì)算，并與MCNP程序的計(jì)算結(jié)果進(jìn)行了比對(duì)，部分比對(duì)結(jié)果見表4，其中為達(dá)到器官吸收劑量統(tǒng)計(jì)量不確定度在1%以內(nèi)，兩者模擬粒子數(shù)均取1×108。

表4 不同程序器官吸收劑量結(jié)果對(duì)比

對(duì)比程序Gadep和MCNP計(jì)算結(jié)果，多數(shù)器官劑量值相對(duì)差異小于1%，少數(shù)器官(小器官、較遠(yuǎn)器官)相對(duì)差異在3%左右，從而驗(yàn)證了2種程序計(jì)算結(jié)果的一致性。事故中病人受到最大照射的器官是性腺，MCNP5和Gadep程序性腺吸收劑量計(jì)算值在10.46 Gy左右，超過(guò)了輻射致永久性不育的吸收劑量限值，而郭凱琳等評(píng)估病人睪丸受照劑量約為9.16 Gy[19]，在受照后92 d精子存活率降為0，表明臨床診斷和物理劑量重建計(jì)算結(jié)果一致。

在加速效率方面，如表5所示，與MCNP單核CPU計(jì)算模擬相比，Gadep程序計(jì)算加速效率達(dá)到50倍以上。

表5 不同程序器官吸收劑量計(jì)算時(shí)間對(duì)比

5 結(jié)論

1)本文開發(fā)了一套GPU加速光電耦合輸運(yùn)蒙卡程序Gadep，采用完整的光子、電子精細(xì)物理模型，并實(shí)現(xiàn)了跨幾何邊界輸運(yùn)、結(jié)構(gòu)體表訪問(wèn)等算法優(yōu)化。

2)驗(yàn)證了Gadep程序的正確性。以ICRP 116號(hào)報(bào)告成年男性參考人外照射器官劑量轉(zhuǎn)換系數(shù)計(jì)算為例，程序Gadep計(jì)算結(jié)果和參考數(shù)據(jù)一致。

3)相比單核CPU程序MCNP5計(jì)算時(shí)間，Gadep程序在外照射器官劑量轉(zhuǎn)換系數(shù)計(jì)算中可以實(shí)現(xiàn)最高300倍的加速效果，但加速效率是和入射粒子類型、能量相關(guān)。

程序還存在較大的優(yōu)化空間，需要探索進(jìn)一步提高程序加速效率的方法，例如：通過(guò)光子、電子分開模擬進(jìn)一步減少線程歧離；在異構(gòu)快速蒙特卡羅程序中實(shí)現(xiàn)串行程序成熟的減方差技巧等加速算法。此外，目前程序中幾何模塊主要實(shí)現(xiàn)了基于長(zhǎng)方體體素化建模，基本能滿足放射治療、人員劑量評(píng)估等模擬需求，但對(duì)于復(fù)雜幾何建模，程序還存在一定的局限，后續(xù)將開展復(fù)雜幾何模型的實(shí)現(xiàn)工作，從而進(jìn)一步拓展程序的適用范圍。