基于EGSnrc程序的錐形束CT模擬仿真研究

付 娟，孟慧鵬，*，艾克文，顧 欣*，王克強(qiáng)，張艷龍，丁紅軍，段敬豪，孫 倩

（1.武警特色醫(yī)學(xué)中心，天津 300162；2.天津大學(xué)精密儀器與光電子工程學(xué)院，天津 300072）

0 引言

機(jī)載式錐形束CT已經(jīng)廣泛用于調(diào)強(qiáng)放射治療（intensity modulated radiotherapy，IMRT）的圖像引導(dǎo)，成為位置驗(yàn)證及跟蹤腫瘤變化等方面的有力工具[1-3]。隨著放療技術(shù)向自適應(yīng)放射治療（adaptive radiotherapy，ART）發(fā)展，越來越多的學(xué)者開始關(guān)注基于錐形束CT進(jìn)行形變配準(zhǔn)、劑量計(jì)算等新的需求[4-6]。然而為方便此類研究的順利開展，首先需要解決的問題是如何對(duì)現(xiàn)有錐形束CT進(jìn)行有效的模擬仿真。因此本研究擬通過Monte Carlo模擬程序EGSnrc對(duì)錐形束CT進(jìn)行模塊化建模，并對(duì)所建立的系統(tǒng)進(jìn)行驗(yàn)證，進(jìn)而探討所建立的仿真系統(tǒng)是否能夠準(zhǔn)確模擬錐形束CT的主要技術(shù)指標(biāo)，及其是否具有進(jìn)一步研究其他錐形束CT相關(guān)問題的潛能。

1 材料和方法

1.1 儀器設(shè)備

實(shí)驗(yàn)設(shè)備為Varian IX直線加速器（美國，瓦里安公司）；三方能譜文件由開源程序Spektr（Version 3.0）生成；測量時(shí)采用型號(hào)為MT-VW的紅固體水，密度為1.03 g/cm3；百分深度劑量測量使用德國PTW公司的MP3-M大水箱，離軸劑量測量使用該公司的OCTAVIUS Detector 729二維電離室矩陣；模擬仿真程序的版本為EGSnrc v2016（2016年1月發(fā)布），在Ubuntu 16.04的64位系統(tǒng)下運(yùn)行[7]。

1.2 方法

1.2.1 主要部件參數(shù)

擬模擬仿真的錐形束CT系統(tǒng)集成于Varian IX直線加速器上，主要組成構(gòu)件為X射線源和非晶硅平板探測器。X射線源球管為Varian公司生產(chǎn)，型號(hào)為G242，峰值電壓和陽極最大熱容量分別為125 kVp和445 kJ，陽極材質(zhì)采用鎢合金，靶光束角為14°，初級(jí)濾波為2 mm厚的鋁合金；平板探測器同樣為Varian公司生產(chǎn)，型號(hào)為4030CB，有效探測面積為39.7 cm×29.8 cm，矩陣尺寸為 1 024×768，像素尺寸為0.388 mm×0.388 mm，探測器可沿探測平面橫向向外偏移 16 cm，最大掃描視野（field of view，F(xiàn)OV）可達(dá)到46.6 cm[8]。

1.2.2 錐形束CT的模擬仿真

本研究通過EGSnrc的BEAMnrc子程序?qū)㈠F形束CT系統(tǒng)分為X射線源、射線源組件、探測器、Monte Carlo模擬用模體共計(jì)4個(gè)模塊建模，所建立的系統(tǒng)示意圖如圖1所示。根據(jù)掃描時(shí)采用的領(lǐng)結(jié)濾波不同，分為全扇束模式和半扇束模式。實(shí)驗(yàn)使用的錐形束CT系統(tǒng)掃描協(xié)議共有3種掃描電壓，而125 kVp對(duì)應(yīng)2種領(lǐng)結(jié)濾波模式，因此共建立4種源模型：100 kVp全領(lǐng)結(jié)濾波、110 kVp半領(lǐng)結(jié)濾波、125 kVp全領(lǐng)結(jié)濾波和半領(lǐng)結(jié)濾波。

圖1 EGSnrc模擬的CBCT系統(tǒng)示意圖

EGSnrc的BEAMnrc子程序可提供多種建模組件模塊，通過不同組件的組合可實(shí)現(xiàn)錐形束CT的模擬仿真?？捎糜谀M錐形束CT系統(tǒng)的組件有XTUBE、CONSTAK、SLABS、BLOCK、JAWS和PYRAMIDS等，模擬仿真時(shí)，由SLABS組件構(gòu)建初級(jí)濾波、射束硬化濾波、機(jī)載影像系統(tǒng)窗和探測器4種結(jié)構(gòu)；領(lǐng)結(jié)濾波既可以用PYRAMIDS組件構(gòu)建，也可以用JAWS組件構(gòu)建，本研究中采用PYRAMIDS組件。每個(gè)組件的最終輸入?yún)?shù)如下：

（1）X射線源。

①XTUBE組件構(gòu)建靶：由Varian公司提供的球管參數(shù)文檔可知，陽極靶物質(zhì)為鎢合金，其中主要成分為金屬鎢（占比95%）、金屬銠（占比5%），等效物理密度為18.68 g/cm3。此外，靶光束角為14°，入射電子的焦距為0.4~0.8 mm，本實(shí)驗(yàn)?zāi)M定義的焦斑直徑為0.5mm。②CONSTAK組件構(gòu)建球管出窗：玻璃窗的密度設(shè)置為2.23 g/cm3，厚度設(shè)置為0.9 mm；窗口與靶之間的距離設(shè)置為2.8cm，出窗與靶之間設(shè)置為真空。

（2）射線源組件。

①SLABS組件構(gòu)建初級(jí)濾波：初級(jí)濾波為2 mm厚的鋁合金，與靶的距離為5.72 cm。其中鋁合金的等效物理密度為3.23 g/cm3，包含96.5%的鋁、1.3%的硅、1.2%的鎂和1%的錳。②BLOCK組件構(gòu)建初級(jí)準(zhǔn)直器：初級(jí)準(zhǔn)直器的材質(zhì)為2 cm厚的鉛，其上端開口（靠近源方向）和下端開口（遠(yuǎn)離源方向）分別為3.2cm×3.2cm和4.2cm×4.2cm，與靶的距離為6.52cm。③JAWS組件構(gòu)建擋片：擋片的材質(zhì)為0.2 mm厚的鉛，同準(zhǔn)直器的設(shè)置類似，上下端開口分別設(shè)置為4.72 cm×4.72 cm和4.86 cm×4.86 cm，與靶的距離為9.59 cm。④SLABS組件構(gòu)建射束硬化濾波：射束硬化濾波為0.5 mm厚的鈦，物理密度為4.54 g/cm3，與靶的距離設(shè)置為11.19 cm。⑤SLABS組件構(gòu)建機(jī)載影像系統(tǒng)窗：機(jī)載影像系統(tǒng)窗為1 mm厚的透光有機(jī)物聚碳酸酯，與靶的距離設(shè)置為14.9 cm。⑥PYRAMIDS組件構(gòu)建領(lǐng)結(jié)濾波：構(gòu)建全扇束領(lǐng)結(jié)濾波和半扇束領(lǐng)結(jié)濾波，領(lǐng)結(jié)濾波的材質(zhì)為鋁合金，與初級(jí)濾波的鋁合金成分相同。建模時(shí)由6層鋁層組成（每層中心開孔不同），從上到下的鋁層厚度分別為 1.5、1.5、7、13、3 和 2 mm，與靶的距離為 15 cm。

（3）探測器。

實(shí)驗(yàn)采用的是Varian 4030CB平板探測器，BEAMnrc子程序建模時(shí)主要考慮兩個(gè)部分：一個(gè)是裝配于探測器陣列上方的10∶1的抗散射柵，另一個(gè)是平板探測器銅層和由無定形硅構(gòu)成的平板薄膜晶體管陣列（閃爍體層）。擺位采集到的錐形束CT圖像原始投影信號(hào)實(shí)際上就是探測器陣列對(duì)到達(dá)探測器粒子的響應(yīng)，建模時(shí)X射線的響應(yīng)可通過提取該層的能量沉積實(shí)現(xiàn)。對(duì)于千伏級(jí)別的X射線，抗散射柵自身和平板探測器銅層與粒子的相互作用不可忽略，因此將探測器建模成1 mm厚的銅層，采用SLABS組件構(gòu)建，物理密度設(shè)置為8.9 g/cm3，與靶的距離為150cm，設(shè)置平板探測器的大小為39.7cm×29.8cm，矩陣尺寸為1024×768，像素尺寸為0.388mm×0.388mm。

（4）Monte Carlo模擬用模體。

為了驗(yàn)證BEAMnrc子程序建立系統(tǒng)的準(zhǔn)確性及后續(xù)實(shí)驗(yàn)，建立的固體水模尺寸為30 cm×30 cm×30 cm，同時(shí)建立相同尺寸的標(biāo)準(zhǔn)水模體。

1.2.3 仿真系統(tǒng)的驗(yàn)證

為了驗(yàn)證建立系統(tǒng)的幾何和光束特性是否正確，通過X射線能譜、百分深度劑量、離軸劑量等進(jìn)行分析比較。Monte Carlo模擬上述參數(shù)的結(jié)果分兩步生成：第一，在BEAMnrc子程序中相應(yīng)平面生成相空間文件；第二，由相應(yīng)子程序分析相空間文件得到對(duì)應(yīng)結(jié)果，其中BEAMDP子程序可獲得X射線能譜，DOSXYZnrc子程序可獲得百分深度劑量和離軸劑量。比對(duì)數(shù)據(jù)的生成如下：由開源程序Spektr生成X射線能譜，由實(shí)測數(shù)據(jù)獲取百分深度劑量和離軸劑量。

（1）X射線能譜。

用于比對(duì)的能譜文件由開源程序Spektr（Version 3.0）生成，分別生成100、110和125 kVp的能譜文件。具體設(shè)置時(shí)，在Spektr程序的UI界面中選擇Tube_09，并在tubeSettings函數(shù)中將Tube_09的值自定義為[13 2；74 0]，在“ADDED FILTERATION”模塊中，增加Ti（鈦，Z=22），“Thickness”設(shè)置為0.5 mm，以確保生成的能譜與本實(shí)驗(yàn)設(shè)備的參數(shù)盡可能一致。

Monte Carlo模擬X射線能譜時(shí)，對(duì)應(yīng)錐形束CT系統(tǒng)掃描協(xié)議的3種X射線管電壓，實(shí)驗(yàn)驗(yàn)證進(jìn)行了3次主要的Monte Carlo模擬。每次運(yùn)行均得到對(duì)應(yīng)的相空間文件，該文件生成于模擬仿真中的一個(gè)虛擬平面，可存儲(chǔ)所有通過該平面的粒子的能量、位置、方向等相關(guān)數(shù)據(jù)，由BEAMDP子程序分析后得到X射線能譜。由于Spektr程序只能生成純鎢靶X射線能譜，且只能生成距離源100 cm處的能譜文件，為了便于比較，在生成分析X射線能譜的相空間文件模擬中，本實(shí)驗(yàn)將所建立系統(tǒng)中的靶材料替換為純鎢（其余模擬均為鎢合金），相空間文件記錄平面設(shè)置為距離靶100 cm處。X射線能譜的Monte Carlo模擬采用的粒子數(shù)為1×109，在距離靶100 cm處模擬的射野大小為50 cm×50 cm，模擬時(shí)未加載領(lǐng)結(jié)濾波組件。

（2）百分深度劑量和離軸劑量。

實(shí)測數(shù)據(jù)均需錐形束CT射線源在0°處，因此選擇設(shè)備的服務(wù)模式，該模式可自定義相關(guān)掃描參數(shù)，如管電壓、管電流、掃描時(shí)間、有無領(lǐng)結(jié)濾波、擋片開口大小等，所有數(shù)據(jù)均在透視模式下源皮距（source skin distance，SSD）100 cm處測量得到。為驗(yàn)證所建立的仿真系統(tǒng)的準(zhǔn)確性，分別實(shí)測了100kVp全領(lǐng)結(jié)濾波、110 kVp半領(lǐng)結(jié)濾波、125 kVp全領(lǐng)結(jié)濾波和半領(lǐng)結(jié)濾波的百分深度劑量和離軸劑量。百分深度劑量測量時(shí)設(shè)置擋片野大小為26 cm×26 cm，在射束穩(wěn)定狀態(tài)下，采集射束中心軸0~20 cm范圍內(nèi)的劑量；離軸劑量測量時(shí)將0.5 cm厚的固體水放置于電離室矩陣之上，采集固體水下0.75 cm平面內(nèi)等中心位置橫向和縱向的劑量，射野大小為20 cm×20 cm。

Monte Carlo模擬百分深度劑量和離軸劑量時(shí)，對(duì)應(yīng)錐形束CT系統(tǒng)實(shí)測數(shù)據(jù)的4種電壓-濾波模式進(jìn)行4次主要的Monte Carlo模擬，共記錄4組相空間文件，記錄位置為源到等中心距離（source isocentre distance，SID）100 cm 處，Monte Carlo 模擬使用的粒子數(shù)為1×109。記錄的相空間文件使用DOSXYZnrc子程序分析計(jì)算，使用的模體大小為30 cm×30 cm×30 cm，計(jì)算體素大小為 0.5 cm×0.5 cm×0.5 cm，粒子數(shù)為5×108，光子分裂數(shù)設(shè)置為100，對(duì)所有感興趣區(qū)域的體素計(jì)算中統(tǒng)計(jì)不確定性控制在2%以內(nèi)。DOSXYZnrc子程序的輸出文件格式為*.3ddose，該文件可由STATDOSE子程序分析。提取出射束中心軸0~20 cm范圍內(nèi)的劑量和模體深度為1.25 cm平面內(nèi)等中心位置橫向和縱向的劑量，繪制生成百分深度劑量和橫向、縱向離軸劑量曲線與實(shí)測結(jié)果進(jìn)行比較。

2 結(jié)果

在SID 100 cm平面處，3種管電壓100、110和125 kVp的Monte Carlo模擬結(jié)果與Spektr程序生成的對(duì)應(yīng)比較結(jié)果如圖2所示，所有能譜信息均與峰值信息做了歸一化處理。結(jié)果表明，2種能譜具有良好的一致性，差異為0.6%~3.0%，平均差異小于2.0%。由圖2還可知，3種管電壓的最大差異區(qū)域均出現(xiàn)在30~55 keV范圍內(nèi)，其中125 kVp管電壓的36 keV能量處差異達(dá)到了3.0%。

圖2 仿真系統(tǒng)生成的能譜信息與Spektr程序生成的能譜信息的比較

Monte Carlo模擬的百分深度劑量和離軸劑量結(jié)果與實(shí)測的結(jié)果一致性良好，如圖3~5所示，所有的數(shù)據(jù)均對(duì)最大值做了歸一化處理。圖3為百分深度劑量分布的比較，4種電壓-濾波模式的射束中心軸劑量差異均較小，劑量差異平均值在±2.0%以內(nèi)。其中100 kVp全領(lǐng)結(jié)濾波差異最大，且表現(xiàn)為全領(lǐng)結(jié)濾波模式[如 3 圖（a）和（c）所示）]的差異大于半領(lǐng)結(jié)濾波模式[如圖 3（b）和（d）所示]，且水下 5 cm厚度以內(nèi)的Monte Carlo模擬值與實(shí)測值的差異大于5 cm以外的。差異最大處在100 kVp全領(lǐng)結(jié)濾波的水下5 cm范圍內(nèi)，達(dá)到了2.7%；4種模式水下5 cm以外差異值均逐漸減小，均小于1.0%。圖4為橫向離軸劑量分布，用于驗(yàn)證系統(tǒng)加載全領(lǐng)結(jié)濾波和半領(lǐng)結(jié)濾波過濾器的幾何設(shè)計(jì)，4種電壓-濾波模式橫向離軸劑量的差異平均值在±1.1%以內(nèi)。由圖4可知，110 kVp半領(lǐng)結(jié)濾波模式Monte Carlo模擬值與實(shí)測值在4種模式中差異最大，最大處達(dá)到了4.0%，該模式中的差異平均值為1.3%。圖5為縱向離軸劑量分布比較，4種電壓-濾波模式縱向離軸劑量的差異平均值在±0.8%以內(nèi)，最大差異為2.7%。由圖5可知，4種模式的曲線有類似的形狀，其中125 kVp全領(lǐng)結(jié)濾波和半領(lǐng)結(jié)濾波模式的曲線肉眼無法區(qū)分其差異。

圖3 仿真系統(tǒng)中心軸百分深度劑量曲線的Monte Carlo模擬數(shù)據(jù)與實(shí)測數(shù)據(jù)比較

圖4 仿真系統(tǒng)等中心位置橫向離軸劑量曲線的Monte Carlo模擬數(shù)據(jù)與二維電離室矩陣實(shí)測數(shù)據(jù)比較

3 討論

ART技術(shù)最終向劑量引導(dǎo)放射治療（dose-guided radiotherapy，DGRT）[9]發(fā)展，實(shí)現(xiàn)實(shí)際患者體內(nèi)劑量沉積的實(shí)時(shí)追蹤[4]。通過治療時(shí)采集的患者的擺位錐形束CT圖像可間接實(shí)現(xiàn)患者體內(nèi)累積劑量的追蹤，是ART技術(shù)當(dāng)前較易實(shí)現(xiàn)的一種方法。然而，基于錐形束CT進(jìn)行形變配準(zhǔn)、劑量計(jì)算的研究[5-6]必須基于一套準(zhǔn)確高效的模擬仿真系統(tǒng)。本研究在充分了解Varian OBI系統(tǒng)（錐形束CT系統(tǒng)）的構(gòu)造及特點(diǎn)的基礎(chǔ)上，基于Monte Carlo模擬軟件的EGSnrc程序?qū)崿F(xiàn)了錐形束CT系統(tǒng)的模擬仿真。所建立系統(tǒng)計(jì)算得到的百分深度劑量、離軸劑量等物理參數(shù)與實(shí)測值誤差平均值均在±2.0%以內(nèi)，與文獻(xiàn)[8，10-11]中其他學(xué)者的建模方法相比，各種幾何結(jié)構(gòu)及物質(zhì)組成更接近真實(shí)情況，如采用鎢合金和鋁合金代替其他學(xué)者采用純鎢和純鋁的模擬、考慮平板探測器的1 mm銅層的模擬等。此外，通過改變各組件有限的參數(shù)，所建立的系統(tǒng)可對(duì)其他主流錐形束CT系統(tǒng)進(jìn)行模擬，具有良好的通用性，可用于模擬研究其他錐形束CT系統(tǒng)及錐形束CT系統(tǒng)的其他問題。

值的一提的是，實(shí)驗(yàn)中模擬生成的能譜未能與實(shí)際測量的能譜數(shù)據(jù)進(jìn)行比對(duì)，主要原因是能譜的測量復(fù)雜，實(shí)際測量較困難，且測量精度受多種因素制約，因此本實(shí)驗(yàn)采用很多研究使用的成熟的第三方能譜生成軟件來生成比對(duì)能譜信息。能譜生成程序采用的是Spektr最新版本，其基于三次樣條插值算法的鎢陽極光譜模型（tungsten anode spectral model using interpolating cubic splines，TASMICS）的算法，相對(duì)于之前采用多項(xiàng)式插值的鎢陽極光譜模型（tungsten anode spectral model using interpolating polynomials，TASMIP）算法的老版本進(jìn)行了改進(jìn)，3.0版本允許用戶自定義球管參數(shù)、不同厚度材料的過濾器及對(duì)生成的能譜文件優(yōu)化，可生成20~150 keV能量范圍內(nèi)的X射線能譜，其模擬出的能譜信息與實(shí)際數(shù)據(jù)的差異在4.0%以內(nèi)，平均差異在2.7%以內(nèi)[10]。本中心Ming等[11-12]的前期工作表明，10%的能譜差異可導(dǎo)致約4%的射束中心軸百分深度劑量差異和小于1%的器官劑量差異，因此仿真系統(tǒng)生成的能譜文件中射束中心軸百分深度劑量差異均較小，可滿足錐形束CT系統(tǒng)的相關(guān)研究。對(duì)于百分深度劑量曲線，在5 cm以內(nèi)深度處差異較大的原因與模擬時(shí)未考慮實(shí)際照射的多向散射及電離室的精度有關(guān)，其均會(huì)帶來淺表區(qū)域的差異增加；對(duì)于能譜信息30~55 keV能量范圍內(nèi)的差異大于其他區(qū)域的可能原因是初級(jí)濾波和射束硬化濾波對(duì)低能射線的過濾放大了該區(qū)域的差異，此外本實(shí)驗(yàn)建模的1 mm銅層探測器也對(duì)低能部分影響較大。

圖5 仿真系統(tǒng)等中心位置縱向離軸劑量曲線的Monte Carlo模擬數(shù)據(jù)與二維電離室矩陣實(shí)測數(shù)據(jù)比較

橫向離軸劑量曲線全領(lǐng)結(jié)濾波模式在離軸距離[-10 cm，10 cm]區(qū)間內(nèi)差異略大于[-13 cm，-10 cm）和（10 cm，13 cm]區(qū)間，可能的原因是模塊化建模時(shí)忽略支撐領(lǐng)結(jié)濾波的鋼支架（厚度小于1 mm）導(dǎo)致的。全領(lǐng)結(jié)濾波結(jié)構(gòu)為中間開孔，兩邊厚，鋼支架對(duì)中心位置的劑量影響被放大，從而導(dǎo)致領(lǐng)結(jié)濾波開孔區(qū)域內(nèi)的劑量差異變大。半領(lǐng)結(jié)濾波模式在離軸距離[-13 cm，6 cm]區(qū)間內(nèi)的差異略大于（6 cm，13 cm]區(qū)間，原因與全領(lǐng)結(jié)濾波類似。盡管存在差異，但是橫向離軸劑量的差異整體在較小的范圍內(nèi)，對(duì)使用該仿真系統(tǒng)進(jìn)行錐形束CT的其他研究影響甚微。

縱向離軸劑量曲線中，125 kVp全領(lǐng)結(jié)濾波和125 kVp半領(lǐng)結(jié)濾波模式的曲線肉眼無法區(qū)分其差異，究其原因與領(lǐng)結(jié)濾波的結(jié)構(gòu)有關(guān)。領(lǐng)結(jié)濾波的中心開孔方向與縱向平行，導(dǎo)致125 kVp全領(lǐng)結(jié)濾波采集數(shù)據(jù)時(shí)領(lǐng)結(jié)濾波對(duì)其造成的影響與半領(lǐng)結(jié)濾波相同（僅與鋼支架發(fā)生作用）。領(lǐng)結(jié)濾波在[-10 cm，10 cm]區(qū)間內(nèi)縱向離軸劑量模擬值與實(shí)測值吻合性更好，其余區(qū)域差異略大，究其原因是實(shí)際測量使用的OCTAVIUS Detector 729二維電離室矩陣相鄰2個(gè)電離室的中心間距為1 cm，而模擬時(shí)的體素計(jì)算間距為0.5 cm，分辨力的差異造成劑量變化較快的區(qū)域出現(xiàn)更大誤差，此外差異可能還與二維電離室矩陣的測量精度有關(guān)。分析該組數(shù)據(jù)還發(fā)現(xiàn)，Monte Carlo模擬值與實(shí)測值差異最大的模式為125 kVp全領(lǐng)結(jié)濾波，出現(xiàn)在離軸距離約-12 cm處，差異值為2.7%。同樣，縱向離軸劑量的差異整體上也在較小的范圍內(nèi)，4種電壓-濾波模式縱向離軸劑量的差異平均值在±0.8%以內(nèi)，對(duì)使用該仿真系統(tǒng)進(jìn)行錐形束CT的其他研究時(shí)的影響比橫向的影響更小。

綜上所述，本研究建立的模擬仿真系統(tǒng)可實(shí)現(xiàn)錐形束CT的準(zhǔn)確模擬，可為研究錐形束CT的相關(guān)問題奠定基礎(chǔ)。下一步將基于該系統(tǒng)對(duì)如何提高錐形束CT圖像質(zhì)量及基于該系統(tǒng)的形變配準(zhǔn)算法等展開研究。