劉穎,孟婷,路鶴晴,章浩偉△
(1.上海理工大學 健康科學與工程學院, 上海 200093;2. 同濟大學附屬第一婦嬰保健院 設備科,上海 200092)
醫(yī)學輻射是電離輻射的主要來源之一[1]。目前,X射線被用于醫(yī)療照射所產(chǎn)生的輻射劑量在人體的整體劑量中占比最大[2],其中,CT檢查是輻射劑量的主要來源[3]。國際放射防護委員會ICRP102號報告指出,在過去20年內(nèi),CT檢查的使用頻率在世界范圍內(nèi)增加了逾8倍[4]?!笆晃濉逼陂g,關于上海市醫(yī)學輻射水平的調研表明,2009年比1996年的X-CT檢查,大幅度凈增3.17倍[5]。因此,如何準確計算CT受檢者所受的輻射劑量顯得尤為重要。
在各類輻射劑量估算方法中,蒙特卡洛方法是一類以概率論為中心思想、具備隨機性和不確定性的統(tǒng)計方法[6],被廣泛應用于生物醫(yī)學、粒子輸運計算、核工程等領域。常見的蒙特卡洛程序包括MCNP (Monte Carlo Neutron and Photo Transport Code)程序、ESG (Electron-Gamma Shower)程序、FLUKA (FLUktuierende KAskade) 程序等。利用蒙特卡洛程序準確計算人體在CT檢查時所受的輻射劑量,需要建立準確的CT輻射系統(tǒng)模型,簡稱CT模型。CT模型的具體構建順序為獲取X射線能譜、設計蝶形過濾器形狀、構建準直器模型和模擬球管的旋轉運動。隨著CT技術的不斷提升,國內(nèi)外有許多研究者提出了CT輻射系統(tǒng)模型在MCNP模擬中的構建方法,并且建立了不同類型的CT模型。本文按照CT輻射系統(tǒng)的建模順序分別對這些方法進行綜述,為相關領域的研究者以及后續(xù)的研究方向提供參考。
首個完整的CT輻射系統(tǒng)模型由澳大利亞的Caon等[7]于1999年建立,他們利用EGS程序建立了一臺GE Hi-Speed Advantage CT模型,并對MIRD (Medical Internal Radiation Dose)人體模型進行了劑量估算。隨后,英國的Khursheed等[8]利用MCNP分別建立了西門子DRH、GE 9800和飛利浦LX模型,并對從新生兒到成年人的各階段人群進行了劑量估算與對比。2003年,美國的Jarry等[9]提出了單排螺旋CT的建模方法,并利用MCNP構建了GE HiSpeed CT/i CT模型,且對MIRD 體模進行了劑量驗證。在此基礎上,同研究團隊的DeMarco等[10]于2005年首次構建了一臺GE LightSpeed 16排CT模型,并對標準體模進行了劑量驗證。同年,希臘的Tzedakis等[11]利用MCNP建立了西門子 Sensation多排螺旋CT的輻射系統(tǒng)模型,并根據(jù)實際掃描條件計算了體模所受的劑量。西班牙的Salvado[12]利用EGS建立了GE HiSpeed LX/i CT模型,并對仿真體模和人體體素模型的輻射劑量進行估算。2009年,美國的Gu等[13]首次詳細描述了GE LightSpeed 16排CT的建模技術路線,代表著完整CT建模體系的形成。隨后的研究都是在此基礎上進行建模方法的改進。2011-2013年間,美國的Lee及其團隊[14-15]不斷提出了CT模型的改進方法,并增強了西門子Sensation 16排CT的準確度。直到2014年,清華大學的潘羽晞等[16]在Gu等[13]的技術路線基礎上,利用MCNP建立了GE LightSpeed 16排CT模型,并首次估算了一個符合中國人體征的兒童模型所受的輻射劑量。2014年至今,研究者們致力于研究CT建模的每個環(huán)節(jié),分別從精準的測量儀器、物理測量方法和理論算法角度提升模型的準確度。
X射線譜包含軔致輻射和特征輻射。當入射高速電子與原子核碰撞,其部分或全部動能變?yōu)楫a(chǎn)生的光子的能量,由此產(chǎn)生的X射線稱為軔致輻射。X射線管的管電壓超過某一臨界電壓時,在某些特定能量值處會出現(xiàn)強度很高、非常狹窄的譜線,它們疊加在連續(xù)譜上,即為特征輻射[17]。在蒙特卡洛模擬中,為簡化放射源模型,通常將初始X射線穿過固有濾過后的能譜視作放射源能譜,用以描述光子能量的抽樣分布。獲取X射線能譜的方法主要分為以下兩種:
(1)由廠家提供球管幾何參數(shù),輸入并由軟件生成X射線能譜。
從上個世紀80年代以來,陸續(xù)有研究者開發(fā)了X射線能譜軟件(見表1)。這些軟件的理論原理分別來自半經(jīng)驗模型、蒙特卡洛模擬和兩種方法的結合??赏ㄟ^在軟件中輸入廠家提供的陽極靶角、固有濾過厚度等幾何參數(shù)來獲得X射線能譜。Caon等[7]利用Tucker等[18]提出的半經(jīng)驗模型來獲得能譜。Khursheed等[8]、Salvado等[12]和Lee等[14-15]利用IPEM78[19]和廠家提供的數(shù)據(jù)生成CT模型的X射線能譜。Tzedakis等[11]利用TASMIP[20]來生成X射線能譜;Gu等[13]和潘羽晞等[16]利用了XCOMP5R[21]生成X射線能譜。
表1 現(xiàn)有的X射線能譜軟件及其原理
Rogers等[22]通過對比實驗值和蒙特卡洛模擬值,驗證了由蒙特卡洛模擬生成的能譜精度要高于半經(jīng)驗模型生成的能譜。Poludniowski等[17]通過在不同能譜軟件中輸入相同的球管參數(shù),發(fā)現(xiàn)基于蒙特卡洛模擬的能譜軟件之間的輸出結果最大差異小于3%,而基于半經(jīng)驗模型的能譜軟件的最大差異為6%。由此可知,在模擬X射線源時,采用半經(jīng)驗模型生成的能譜會給結果帶來較大誤差。
(2)測量X射線球管在一定管電壓下的半值層,在能譜軟件中匹配測量的半值層來獲得X射線等效能譜。
半值層(half-value layer,HVL)為使輻射束的空氣比釋動能減少到其初始值一半時,指定材料的厚度,通常用mmAL表示[23]。半值層是X射線能量表達的重要參數(shù),因此可以通過匹配實驗室測量與軟件輸出的半值層來獲得X射線等效能譜。等效能譜的概念由Turner等[24]于2009年提出。他們對比了基于廠家給出的數(shù)據(jù)所構建出來的源模型和基于半值層獲得的源模型,發(fā)現(xiàn)前者與實際測量所得的CT劑量指數(shù)(CTDI)的均方根誤差最高可達20%,而后者最高僅為7%,這代表等效能譜可在很大程度上提升模型的準確性,而僅靠廠家提供的數(shù)據(jù)會給模型帶來較大誤差。
人體的橫截面類似于一個橢圓,當X射線穿過人體時,中心處吸收的能量要高于邊緣處,這導致邊緣皮膚劑量增大,因此,需把過濾器設計成蝶形,才能使能量分布均勻。由于蝶形過濾器的形狀復雜,且不同CT之間差異巨大,因此該結構為CT建模中最難的部分。為了在蒙特卡洛模擬中設計出蝶形過濾器的具體形狀,需要獲得其厚度角分布,即過濾器厚度關于X射線角度的函數(shù)。獲得蝶形過濾器厚度角分布依賴于實驗室測量,方法大致有3種:
(1)測量標準體模中心與邊緣的CTDI比值,簡化蝶形濾波器模型,再通過調整簡化模型來使模擬值貼近測量值。
Gu等[13]和潘羽晞等[16]在構建GE LightSpeed 16排CT模型時,利用長方體減去一個橢圓柱的方式來構建蝶形過濾器的簡化模型,見圖1。通過不斷調整橢圓的長短軸來使模擬與測量的標準體模中心與邊緣的CTDI比值偏差小于5%,這時認為該簡化模型可代替蝶形過濾器模型。
圖1 蝶形過濾器簡化模型(藍色部分)
(2)測量靜態(tài)劑量曲線,計算獲得蝶形濾波器的厚度角分布。
Turner等[24]將X射線球管固定在90°的位置來測量多個偏離等中心點處的空氣比釋動能。各個測量點的間距為5~10 mm,見圖2。根據(jù)測量的結果可以得到橫坐標為X射線角,縱坐標為空氣比釋動能的靜態(tài)劑量曲線。根據(jù)曲線可以得到每個測量點與等中心點的空氣比釋動能比值。該蝶形過濾器的材料為鋁,可以通過查詢NIST報告來獲得鋁對于不同能量光子的衰減系數(shù),再根據(jù)上一步所獲得的等效能譜,利用能量衰減公式來獲得蝶形過濾器厚度角分布。通過模擬實驗測量發(fā)現(xiàn),實驗值與模擬值的均方根誤差小于5%?;谏鲜隼碚?,Lee等[14-15]獲得了西門子Sensation 16排螺旋CT的蝶形濾波器模型,發(fā)現(xiàn)實驗值與模擬值的均方根誤差小于3%;Hassan等[25]也利用了該方法構建了西門子Definition 64排CT的蝶形過濾器模型,且由于劑量計的精度增高,實驗值與模擬值的均方根誤差僅小于1.16%。該CT模型具有兩種不同形狀的蝶形過濾器,通過計算發(fā)現(xiàn),這兩種過濾器的形狀都與橢圓相差甚遠。由此可知,Gu等[13]提出的模型簡化法并不適用于所有CT的蝶形過濾器模型,反而會給模型帶來較大誤差,而通過靜態(tài)劑量曲線的方式來獲得的模型會更加精確。
圖2 靜態(tài)劑量曲線的測量示意圖
(3)測量旋轉劑量積分曲線,計算獲得蝶形濾波器的厚度角分布。
有別于靜態(tài)劑量曲線,旋轉積分劑量曲線是不固定球管的位置,并利用探測器測量球管旋轉過程中的劑量變化,獲得蝶形過濾器的厚度角分布。該測量方法由Boone提出,命名為COBRA(characterization of bowtie relative attenuation,COBRA)法[26]。測量方式見圖3,劑量計位于視野的邊緣,通過X射線的平方反比定律可以獲得測量點與等中心點的空氣比釋動能關系,見式(1)。
圖3 旋轉積分劑量測量圖Fig.3 Schematic diagram of rotating integral dose measurement
(1)
其中,M(α)為測量點的空氣比釋動能,F(xiàn)(θ)為蝶形過濾器導致X射線在θ角的衰減系數(shù),s為X射線源點到等中心點的距離,g為X射線源點到測量點的距離,I0為等中心點的空氣比釋動能。
隨著球管旋轉一圈,對式(1)兩邊求取角度為0到+2π的積分,根據(jù)劑量計的測量值,可以得到F(θ)的函數(shù)關系式,從而獲得蝶形過濾器的厚度角分布。
在Boone的基礎上, McKenney等[27]將該方法應用在一臺西門子乳腺CT的模型構建中,將構建的蝶形過濾器與廠家提供的模型進行對比,發(fā)現(xiàn)二者的厚度角分布差異均小于2%,該建模方法的優(yōu)越性在于可以分析計算出多種材料組成的蝶形過濾器,并且可以達到較好的準確度;而不足在于對劑量計的時間靈敏度要求很高,無法大范圍應用于其他的研究工作,且無法在球管旋轉的情況下獲得半值層。
隨后,Whiting等[28]、Randazzo等[29]和 Hassan等[25]分別從劑量計、測量方式和理論算法來彌補COBRA法的不足,其模擬值和實際測量值均有高度一致性。
在蒙特卡洛模擬中主要是針對前準直器進行建模,CT的前準直器能夠減少受檢者的輻射劑量,并能限定X射線的掃描范圍??梢酝ㄟ^以下方式來建模:
(1)利用蒙特卡洛程序內(nèi)的功能函數(shù),直接實現(xiàn)粒子截斷。
MCNP程序中內(nèi)置有cookie-cutter函數(shù),可以記錄一定范圍內(nèi)的粒子,舍棄超出范圍內(nèi)的粒子。利用這一點,Gu等[13]和潘羽晞等[16]在放射源附近建立了一個立方體來限定X射線的照射范圍,即圖1中的黃色部分。根據(jù)源到長方體右側底面和到等中心點的距離之比可以調整前準直器寬度,因此該方法比較便利。
(2)建立實際的幾何體來限定光子飛行的范圍。
Lee等[14-15]對此準直器的建模方法進行了詳細闡述,并提供了模型示意圖。首先建立幾個圓柱并進行布爾運算,其次在放射源附近形成僅容一定光子通過的縫隙,且不允許光子穿過縫隙外的空間,最后形成準直器模型,該方法的原理與第一種方法相似。
為了生成CT圖像,球管需要通過旋轉采集各角度的信號?,F(xiàn)有的蒙特卡洛程序無法描述一個不斷運動著的放射源,因此,研究者們通過對不同位置源進行抽樣來模擬球管旋轉一周的情況。
Khursheed等[8]分別對比了一周布置18和72個抽樣點的CTDI模擬結果,發(fā)現(xiàn)二者差異較小,認為18個角度抽樣點足以等效模擬球管旋轉一圈的情況。Gu等[13]同樣對比了24和32個抽樣點的CTDI模擬結果,發(fā)現(xiàn)二者差異小于5%,綜合考慮下,他們采用了16個角度抽樣點來模擬球管旋轉一圈的情況,最后獲得了與實驗值相差較小的劑量模擬結果。
本文針對CT輻射系統(tǒng)模型在蒙特卡洛模擬中的構建方法展開了綜述。按照CT建模順序,即獲取X射線能譜、設計蝶形過濾器形狀、構建準直器模型和模擬球管的旋轉運動,分別對國內(nèi)外研究者構建的CT模型進行了討論和對比。為了獲得準確的CT模型,對于X射線能譜,研究者應盡量使用基于蒙特卡洛模擬建立的能譜軟件,并通過匹配實測半值層的方式來獲得X射線等效能譜;對于蝶形過濾器,研究者應結合自身實驗條件,選擇靜態(tài)劑量法或旋轉積分劑量法來獲得厚度角分布;對于準直器模型和球管的掃描運動,目前的模擬方法比較統(tǒng)一,因此研究者可自行選擇。
另外,通過分析發(fā)現(xiàn),目前研究的CT大多數(shù)來自于GE和西門子,不同CT種類內(nèi)部的結構,尤其是蝶形過濾器差異較大,急需提出一種通用的、準確的模型構建方法,應用于各類CT中。該建模方法應融合并改進國內(nèi)外研究者提出的方法,形成一套完整的CT建模系統(tǒng)。考慮到蒙特卡洛程序的輸入文件結構復雜,且CT模型幾何參數(shù),尤其是蝶形過濾器幾何參數(shù)的計算過程繁瑣,因此,該建模系統(tǒng)不僅需要支持CT模型實時三維可視化,還需要根據(jù)輸入的實測物理量來計算幾何參數(shù),自動生成蒙特卡洛程序的輸入文件,降低CT建模工作難度。本文雖然對比討論了目前CT建模較為準確的方法,但是隨著劑量儀測量精度的不斷提高以及CT掃描技術的不斷發(fā)展,CT建模方法也應不斷改進。