GPU加速光電耦合輸運蒙卡程序研發(fā)及應用

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

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

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

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

1 程序框架設計及算法實現(xiàn)

1.1 程序架構設計

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

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

1.2 程序物理模型選擇

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

圖2 光電耦合輸運過程示意Fig.2 Schematic diagram of coupled photon-electron transport process

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

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

1.3 程序實現(xiàn)

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

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

1.4 步長計算方法選擇

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

s=nλλ

(1)

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

nλ=-logη

(2)

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

(3)

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

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

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

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

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

2 正確性檢驗

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

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

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

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

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

圖5 光子AP外照射器官劑量轉換系數(shù)計算結果和ICRP 116號報告結果對比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照射器官劑量轉換系數(shù)計算結果和ICRP 116號報告結果對比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外照射器官劑量轉換系數(shù)計算結果和ICRP 116號報告結果對比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外照射器官劑量轉換系數(shù)計算結果和ICRP 116號報告結果對比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)化及加速效果比對

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

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

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

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

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

3.2 GPU硬件結構優(yōu)化

GPU計算性能一般用浮點數(shù)運算峰值(FLOPS)來表示，一般GPU的單精度FLOPS要比雙精度FLOPS要高。因此在程序計算精度允許的情況下，使用單精度浮點類型可以提高GPU程序的計算效率。同樣以AP照射情形下全身體素徑跡長度分布計算為例，測試了將計算精確度要求較高的長度相關物理量分別設為單、雙精度浮點類型，對計算結果和計算時間的影響，見表3?？梢钥闯觯?0%左右的體素徑跡長度計算結果相對差異小于1%，相應的計算時間減少了44.4%～58.1%。

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

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

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

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

3.4 程序加速效果比對

為對程序Gadep的加速效果進行驗證，對該程序計算AP、PA標準照射條件下光子、電子外照射劑量轉換系數(shù)的計算時間進行統(tǒng)計，并與同等計算條件下使用MCNP進行計算的計算時間進行了比對。其中MCNP中取粒子模擬類型mode P E,為光子電子耦合輸運，光子、電子截斷能量為程序默認值；測試CPU為Intel i7-4710U，GPU顯卡為NVIDIA GeForce 1080Ti。圖11給出了4種計算情景下Gadep程序與MCNP程序的加速效率對比結果。

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

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

4 程序應用

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

表4 不同程序器官吸收劑量結果對比

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

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

表5 不同程序器官吸收劑量計算時間對比

5 結論

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

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

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

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