戴振暉 王學(xué)濤 朱 琳 張 煜 劉小偉

1 (廣東省中醫(yī)院 廣州 510120)

2 (南方醫(yī)科大學(xué) 廣州 510515)

3 (中山大學(xué) 廣州 510275)

蒙特卡羅(Monte Carlo, MC)模擬是放射治療劑量準確計算的常用方法[1,2]。對一個治療射野的模擬，蒙特卡羅程序需追蹤數(shù)十億粒子的輸運和能量沉積過程，非常耗時。因此，在不降低統(tǒng)計誤差的情況下，如何縮短計算時間是一個非常值得研究的課題。OMEGA/BEAM 程序由加拿大國家研究院(National Research Council of Canada, NRC)開發(fā)，是專用于模擬醫(yī)用直線加速器輸出劑量的程序包，本文基于2011年5月18日發(fā)布的版本EGS nrc V42.3.2模擬 Varian 23EX加速器 6 MV X-線，研究BEAMnrc和DOSXYZnrc子程序的參數(shù)設(shè)置對蒙特卡羅計算效率的影響，計算機配置為：英特爾? 酷睿? i7-3770 處理器，主頻3.4 GHz，8 MB緩存。

1 材料與方法

1.1 相空間文件的計算

采用BEAMnrc模擬Varian 23EX醫(yī)用直線加速器6 MV X-線，計算在源皮距為100 cm條件下10 cm×10 cm 射野的相空間文件，加速器治療頭主要部件包括靶、初級準直器、Be窗、均整器、電離室、射野鏡、次級準直器及水體模等[3]。

BEAMnrc中主要參數(shù)設(shè)置為：AE=ECUT=0.70 MeV，光子的截止能量AP=PCUT=0.01 MeV，除鎢靶外全局電子射程截斷能量ESAVE=2.0 MeV(鎢靶0.7 MeV)[4]，邊界穿越算法選擇 PRESTA-I，在EGSnrc中，若粒子能量低于預(yù)先設(shè)置的能量閾值，即沒有足夠的能量使它到最近的區(qū)域邊界，則該粒子的活動就被終止，這樣，將不產(chǎn)生能逃逸該區(qū)域的軔致輻射光子的可能性降到最低。本文主要研究三種方差減少技術(shù)，即均勻軔致輻射分裂(uniform bremsstrahlung splitting, UBS)、選擇軔致輻射分裂(selective bremsstrahlung splitting, SBS)和方向軔致輻射分裂(directional bremsstrahlung splitting, DBS)；在BEAMnrc中，每一次軔致輻射事件產(chǎn)生NBRSPL個光子，每個光子的權(quán)重為 NBRSPL?1，不同方差減小技術(shù)可供選擇的光子分裂數(shù)的范圍不同[5]。在軔致輻射事件中，NBRSPL是基于入射電子的能量和方向，并正比于軔致輻射光子進入射野的幾率。對于UBS，分裂數(shù)NBRSPL設(shè)置為一個固定值，每個軔致輻射事件產(chǎn)生NBRSPL個軔致輻射光子。對于SBS，NBRSPL是變化的，改變分裂數(shù)可最大化進入射野的光子分裂，最小化離開射野的無必要的光子分裂，SBS需設(shè)置背景分裂數(shù)，背景分裂數(shù)一般為分裂數(shù)的十分之一。若打開俄羅斯輪盤賭則UBS和 SBS將分裂更高權(quán)重的軔致輻射光子和來自湮滅事件光子，若關(guān)閉俄羅斯輪盤賭則不分裂更高權(quán)重的光子，進而減小跟蹤消除權(quán)重粒子的CPU時間。DBS需設(shè)置分裂半徑，且分裂半徑要大于射野大小，打開電子分裂可提高相空間文件中電子通量計算的精確性。由于UBS的無方向特性導(dǎo)致很多CPU時間用于跟蹤未進入感興趣射野的分裂光子，另外，SBS也要求額外的CPU時間用于離開射野的軔致輻射光子的背景分裂。DBS則克服了UBS和SBS的缺點，確保射野內(nèi)所有的光子有相同權(quán)重并消除背景分裂。表1為BEAMnrc中三種方差減小技術(shù)的參數(shù)設(shè)置[6,7]，選擇軔致輻射需要輸入背景光子分裂數(shù)，方向軔致輻射的電子分裂有ON和OFF兩種狀態(tài)。

表1 BEAMnrc中三種方差減小技術(shù)參數(shù)設(shè)置表Table 1 Parameters of three variance reduction techniques in BEAMnrc.

1.2 水模體中劑量分布的計算

由 BEAMnrc生成的相空間文件作為DOSXYZnrc模擬的輸入文件，計算水模體中的劑量分布。在DOSXYZnrc中可設(shè)置光子分裂數(shù)Ns和粒子循環(huán)數(shù)Nr，進入水模體的光子初始權(quán)重為W0，分裂 N次后每個光子的權(quán)重為 W0/N，這樣就增加了進入水模體的目標粒子數(shù)。當使用光子分裂時關(guān)閉粒子循環(huán)，使用粒子循環(huán)時關(guān)閉光子分裂，分別設(shè)置分裂數(shù)和循環(huán)數(shù)進行劑量計算。

1.3 模擬效率

蒙特卡羅模擬的效率e定義為：e = 1/s2T，式中，s表示目標數(shù)據(jù)在數(shù)值上的不確定度，T表示CPU的計算時間。

Tp是 DOSXYZnrc模擬過程消耗 CPU的時間，Tj是BEAMnrc模擬過程消耗CPU的時間，平均不確定度為[5]：

式中，Di表示第i個體素的劑量，DDi表示相關(guān)的統(tǒng)計不確定度， Di3 Dmax/2表示只有劑量大于最大劑量50%的體素才用于統(tǒng)計計算。

1.4 X-線射野劑量分布的測量和蒙特卡羅模擬

在BEAMnrc中設(shè)置電子的截止能量AE=ECUT=0.70 MeV，光子的截止能量AP=PCUT=0.01 MeV，電子射程截斷能量ESAVE=2.0 MeV(鎢靶0.7 MeV)，計算采用 DBS方差減少技術(shù)，軔致輻射粒子分裂數(shù) NBRSPL=1000，分裂野半徑 FS=30 cm，在DOSXYZnrc中設(shè)置Ns=40。

為驗證優(yōu)化計算模型的精確度，我們使用電離室進行了相同條件下的測量，并將計算結(jié)果與測量結(jié)果進行比較。在三維水模體中計算 Varian 23EX加速器6 MV X-線在SSD=100 cm、10 cm×10 cm射野的百分深度劑量和不同深度的離軸比曲線，水模體尺寸48 cm×48 cm×41 cm，體素尺寸0.5 cm×0.5 cm×0.5 cm。模擬的粒子數(shù)范圍為5×108，以保證計算結(jié)果的統(tǒng)計誤差均≤3%。在相同條件下，用Blue Phantom三維水箱(IBA公司，德國)測量百分深度劑量和離軸比曲線。圖1為BEAMnrc和DOSXYZnrc之間的數(shù)據(jù)傳輸流程圖。

圖1 BEAMnrc和DOSXYZnrc之間數(shù)據(jù)傳輸流程圖Fig.1 Flow chart of data transmission between BEAMnrc and DOSXYZnrc.

2 結(jié)果

2.1 BEAMnrc的計算效率

BEAMnrc的計算效率如圖2所示，以無軔致輻射分裂的計算效率為基準，由圖可見，無電子分裂的DBS計算效率高于電子分裂的 DBS，當分裂數(shù)為 2500時，計算效率最高；使用俄羅斯輪盤賭的SBS計算效率高于不使用其的 SBS，當分裂數(shù)為2500時，計算效率最高；使用俄羅斯輪盤賭的UBS計算效率高于不使用其的UBS，當分裂數(shù)為750時，計算效率最高。三種軔致輻射的最高計算效率相比，DBS約為SBS的8倍、UBS 的17倍。軔致輻射分裂可大大降低不確定度，延長CPU計算時間，提高了計算效率。

圖2 不同分裂方式對應(yīng)的相對計算效率Fig.2 Relative computing efficiency of different manner.

2.2 DOSXYZnrc的計算效率

DOSXYZnrc的計算效率如圖3所示，以光子分裂數(shù) Ns和粒子循環(huán)數(shù) Nr為零時的計算效率為基準，當光子分裂數(shù)為40時計算效率最高。

圖3 光子分裂數(shù)和粒子循環(huán)數(shù)與相對劑量計算效率的關(guān)系圖Fig.3 Relative efficiency of different splitting numberand cycling number.

2.3 使用與不使用最優(yōu)參數(shù)的效率比較

表2為使用優(yōu)化和未使用優(yōu)化技巧的模擬效率對比，粒子數(shù)為 5×108。由表 2可見，當不使用優(yōu)化技巧時計算時間是76 h，平均不確定度是3%，使用最優(yōu)參數(shù)設(shè)置技巧的計算時間是13 h，平均不確定度是0.25%，計算效率是不使用優(yōu)化技巧的490倍。

2.4 劑量分布測量和計算結(jié)果的比較

圖 4(a)是 Varian 23EX加速器 6 MV X-線10 cm×10 cm射野的百分深度劑量分布圖。由圖可見，在劑量建成區(qū)域，蒙特卡羅模擬和測量值差異明顯，估計與測量電離室本身的特性和側(cè)向電子不平衡有關(guān)；在深度大于劑量最大點深度的區(qū)域，蒙特卡羅模擬和測量值無明顯差異，誤差小于 1%。圖4(b)是6 MV X-線10 cm×10 cm射野不同深度的離軸比，按中心點劑量歸一(為顯示方便，深度為1.5、5、10、20 cm 的離軸比曲線分別乘以 1.00、0.85、0.65、0.35)。由圖可見，在半影區(qū)域(10%–90%的劑量區(qū)域)，蒙特卡羅模擬和測量值差異明顯，這是由于測量電離室本身的特性、測量條件等造成測量結(jié)果存在誤差，在射野半影區(qū)計算劑量低是因為忽略了源電子束的寬度，未考慮源射束在射野半影區(qū)造成的多焦點影響；在劑量平坦區(qū)域，蒙特卡羅模擬和測量值未明顯差異，大部分點的誤差小于1.5%。

表2 使用優(yōu)化和未使用優(yōu)化技巧的模擬效率對比Table 2 Efficiency of the optimized method in beam and dose simulation relative to simulation without optimized method.

圖4 10 cm×10 cm射野百分深度劑量和不同深度的離軸比蒙特卡羅模擬和測量值的比較Fig.4 Comparison of the MC and measured percent depth dose and off-axis ratio in different depths for 10 cm×10 cm field size.

3 討論

對于 Varian 23EX加速器 6 MV X-線，在BEAMnrc中使用無電子分裂的DBS方式，分裂數(shù)為2500時計算效率最高；在DOSXYZnrc中，當光子分裂數(shù)為40，粒子循環(huán)數(shù)為15時，計算效率最高。其他因素(如射線能量、射野尺寸、體素尺寸、準直器鉛門和多葉光柵的設(shè)置、EGSnrc的版本等)對計算效率的影響較小[8,9]。

EGSnrc是功能強大的軟件程序包，是目前蒙特卡羅模擬醫(yī)用加速器的主要程序[10]。VMC/XVMC/VMC++、MCDOSE/MCSIM、DPM 等[11–14]蒙特卡羅程序包使用各種工具提高運算速度，但帶來較大誤差；因此，用戶可按照計算速度和統(tǒng)計誤差的不同要求，選擇不同的蒙特卡羅程序。

本文所述方法適用于一般的蒙特卡羅模擬程序，針對不同程序，要選擇合適的改進效率的方法和參數(shù)，以提高計算效率。

1 Chetty I J, Curran B, Cygler J E, et al. The AAPM task group report No.105: issues associated with clinical implementation of Monte Carlo-based external beam treatment planning[J]. Medical Physics, 2007, 34(12):4818–4853

2 Reynaert N, Van der Marck S C, Schaart D R, et al.Monte Carlo treatment planning for photon and electron beams[J]. Radiation Physics and Chemistry, 2007, 76(4):643–686

3 王學(xué)濤, 陳少文, 戴振暉, 等. MatriXX 二維電離室陣列劑量分布的角度響應(yīng)[J]. 核技術(shù), 2012, 35(2):126–129 WANG Xuetao, CHEN Shaowen, DAI Zhenhui, et al.Angular dependence study on dose distribution of MatriXX 2-D chamber array[J]. Nuclear Techniques,2012, 35(2): 126–129

4 Ding G X, Duggan D M, Coffey C W, et al. First macro Monte Carlo based commercial dose calculation module for electron beam treatment planning—new issues for clinical consideration[J]. Physics in Medicine and Biology,2006, 51(11): 2781–2799

5 Popple R A, Weinberg R, Antolak J A, et al.Comprehensive evaluation of a commercial macro Monte Carlo electron dose calculation implementation using a standard verification data set[J]. Medical Physics, 2006,33(6): 1540–1551

6 Kawrakow I, Rogers D W O, Walters B R B. Large efficiency improvements in BEAMnrc using directional splitting[J]. Medical Physics, 2004, 31(10): 2883–2898

7 Sheikh-Bagheri D, Rogers D W O, Ross C K, et al.Comparison of measured and Monte Carlo calculated dose distributions from the NRC linac[J]. Medical Physics,2000, 27(10): 2256–2266

8 Kawrakow I, Walters B R B. Efficient photon beam dose calculations using DOSXYZnrc with BEAMnrc[J].Medical Physics, 2006, 33(8): 3046–3056

9 Kawrakow I, Mainegra-Hing E, Rogers D W O.EGSnrcMP, the new multi platform version of EGSnrc[J].Medical Physics, 2004, 31: 1731–1735

10 Rogers D W O, Faddegon B A, Ding G X, et al. BEAM:A Monte Carlo code to simulate radiotherapy treatment units[J]. Medical Physics, 1995, 22(5): 503–524

11 Kawrakow I, Fippel M. Investigation of variance reducetion techniques for Monte Carlo photon dose calculation using XVMC[J]. Physics in Medicine and Biology, 2000,45(8): 2163–2183

12 Neuenschwander H, Born J E. A macro Monte Carlo method for electron beam dose calculations[J]. Physics in Medicine and Biology, 1992, 37(1): 107–125

13 Ma C M, Li J, Pawlicki T, et al. MCSIM: A Monte Carlo dose verification tool for radiation therapy treatment planning and beam delivery[J]. Medical Physics, 2002,29(6): 1316–1320

14 張貴英, 曹磊, 鄧君, 等. 15 MV醫(yī)用電子直線加速器光核中子劑量分布的MC模擬及測量[J]. 核技術(shù), 2010,33(1): 35–38 ZHANG Guiying, CAO Lei, DENG Jun, et al. Simulation and measurement of photon neutron dose distribution from a 15 MV X-ray medical electronic linear accelerator[J]. Nuclear Techniques, 2010, 33(1): 35–38