趙 日,劉兆行,劉 娜,王仙祥,張 靜,梁潤成,劉 鑫,令狐仁靜,戴雨玲

(1.中國輻射防護(hù)研究院,太原 030006;2.核藥研發(fā)轉(zhuǎn)化與精準(zhǔn)防護(hù)山西省重點實驗室,太原 030006)

人體輻射劑量計算是指利用計算機(jī)蒙特卡羅(monte carlo,MC)模擬電離輻射在環(huán)境和人體中的輸運(yùn)過程從而得到人體各處所受劑量值的技術(shù)。該仿真計算無需在人體布設(shè)實體劑量測量設(shè)備,且可在現(xiàn)場輻射作業(yè)開展前或完成后進(jìn)行預(yù)測性或回顧性虛擬計算,同時,能給出全身任意處(包括內(nèi)部器官、眼晶體等)的劑量結(jié)果,劑量分布空間分辨率可達(dá)毫米甚至微米級,因此,與傳統(tǒng)佩戴劑量計的劑量監(jiān)測相比更便捷、靈活且精細(xì),近年來越來越受重視,已成為預(yù)測性和回溯性職業(yè)照射人員劑量評價及健康防護(hù)的關(guān)鍵技術(shù)。

人體輻射劑量計算中使用的人體數(shù)字模型是決定計算效果的關(guān)鍵。目前應(yīng)用最廣的是體素模型[1-3],使用體素(voxel)定義人體結(jié)構(gòu),類似于人體斷層CT圖像下使用像素(pixel)描述2維圖像。ICRP 116號報告給出的人體外照射劑量轉(zhuǎn)換系數(shù)就是基于體素模型計算得到的。然而限于體素模型的內(nèi)在構(gòu)造機(jī)制,模型不能進(jìn)行姿態(tài)調(diào)整,只能保持直立姿態(tài),因此嚴(yán)格來說,當(dāng)前基于體素模型的劑量計算只適用于直立人體靜態(tài)劑量評價,無法考慮人體姿態(tài)對受照劑量的影響。這一缺陷在均勻且強(qiáng)度不大的輻射場下尚能被接受,但在事故、應(yīng)急等場景的高度非均勻強(qiáng)場條件下則會帶來顯著的劑量計算偏差。例如美國的Han等[4]報道,使用單一直立姿態(tài)人體模型進(jìn)行事故劑量重建時,器官劑量最大低估達(dá)78%,有效劑量低估為19%;Yoem等[5]研究顯示,非均勻輻射場中,不同姿態(tài)下人體紅骨髓、肺、胃、結(jié)腸、乳腺、性腺等器官的劑量差異顯著,其中性腺差異最大,極端情況下相差可超過兩個數(shù)量級。因此,實際應(yīng)用中亟需一種能夠精確計算受照人體在各種姿態(tài)下所受劑量的新技術(shù)。

近期,ICRP 145號報告給出了一種全新的基于表面約束幾何的參考人數(shù)字計算模型MRCPs (mesh-type reference computational phantoms),為更高精度的人體輻射劑量計算提供了可能。這種新形式的模型被稱為BREP(boundary representative phantom)模型、Mesh模型或面元模型(本文使用“面元模型”代指)[6-7],如圖1所示。得益于底層數(shù)學(xué)形式的靈活性,面元模型同時具備可姿態(tài)調(diào)整和高分辨率兩方面優(yōu)勢,模型所有部分均可以進(jìn)行自由移動、變形,極大的彌補(bǔ)了體素模型的缺陷,且模型空間分辨率沒有下限,可保證組織器官輪廓的光滑性以及對微米級幾何結(jié)構(gòu)的精確描述,突破了體素模型毫米級分辨率下限和鋸齒形輪廓的限制[8-10]。

圖1 人體數(shù)字面元模型示意

目前,國外初步開展了基于面元模型的人體輻射劑量計算新技術(shù)研究。例如美國的Vazquez等[11]利用人體面元模型和動作捕捉系統(tǒng)對一起嚴(yán)重臨界事故進(jìn)行了劑量重建研究。使用動作捕捉系統(tǒng)來對人體動作進(jìn)行記錄,然后對面元模型變形重現(xiàn)這些動作,結(jié)果顯示,新方法得出的劑量值與軀體癥狀相關(guān)性更好。Yeom等[12]計算了行走、坐、彎腰、跪、蹲共5種姿態(tài)在6種照射條件(即AP、PA、LLAT、RLAT、ROT、ISO)下的外照射劑量轉(zhuǎn)換系數(shù)。國內(nèi)僅有少數(shù)相關(guān)報道,這些零星研究對面元模型的變形過程比較粗糙,無法保證變形的物理合理性和幾何光滑性,同時仍需將變形后的面元模型體素化成體素模型才能開展計算,使得面元模型的優(yōu)勢無法體現(xiàn)[13-14]。為此,本文開展了系統(tǒng)性研究,構(gòu)建了精細(xì)化的面元模型變形方法,實現(xiàn)了面元模型的直接高速M(fèi)C計算,為國內(nèi)進(jìn)一步開展、應(yīng)用新一代高精度劑量計算技術(shù)奠定了基礎(chǔ)。

1 人體數(shù)字面元模型變形方法

對面元模型進(jìn)行變形是實現(xiàn)不同姿態(tài)人體劑量計算最關(guān)鍵的技術(shù)環(huán)節(jié)。但面元模型的數(shù)學(xué)結(jié)構(gòu)高度復(fù)雜,且軟組織(肌肉、皮膚、血管、淋巴組織)、內(nèi)部器官的真實變形無法被實際觀察,而目前已有的研究中的變形都太過粗糙,大多是依據(jù)關(guān)節(jié)的旋轉(zhuǎn)而對軟組織、器官進(jìn)行直接旋轉(zhuǎn),與真實物理景象差別很大。所以鑒于骨骼、軟組織、內(nèi)部器官三部分的物理、幾何特征的差異,分別進(jìn)行變形,同時考慮到骨骼帶動軟組織、軟組織帶動器官變形的物理原理,因此按照骨骼、軟組織、器官的順序依次變形,且為降低變形過程的復(fù)雜度,假設(shè)變形是時序進(jìn)行的,對一部分進(jìn)行變形時,其余部分形狀未開始改變。使用ICRP 145號報告給出的男性、女性參考人數(shù)字面元模型MRCP_AF和MRCP_AM,依次建立模型的骨骼、軟組織和內(nèi)部器官形變方法。具體算法如下。

(1)骨骼調(diào)整

讀取待旋轉(zhuǎn)骨骼的整個網(wǎng)格點數(shù)據(jù),確定骨骼關(guān)節(jié)中心,以該處為旋轉(zhuǎn)中心,確定旋轉(zhuǎn)方向和旋轉(zhuǎn)角度,然后計算旋轉(zhuǎn)矩陣,基于該矩陣計算骨骼中各網(wǎng)格頂點的坐標(biāo)旋轉(zhuǎn)變換。

旋轉(zhuǎn)時,任意點的坐標(biāo)按式(1)三維空間中繞任意軸進(jìn)行旋轉(zhuǎn)的旋轉(zhuǎn)矩陣進(jìn)行變換:

(1)

待旋轉(zhuǎn)的關(guān)節(jié)部位主要有15個,即肩關(guān)節(jié)(2個)、肘關(guān)節(jié)(2個)、腕掌關(guān)節(jié)(2個)、髖關(guān)節(jié)(2個)、膝關(guān)節(jié)(2個)和踝關(guān)節(jié)(2個)、脊柱關(guān)節(jié)(簡化為3個,頸關(guān)節(jié)、腰關(guān)節(jié)、骶關(guān)節(jié))。變形效果如圖2所示。

圖2 本文實現(xiàn)的骨骼旋轉(zhuǎn)、變形效果

(2)軟組織變形

接下來對各關(guān)節(jié)周圍的軟組織(肌肉、皮膚、血管)進(jìn)行變形。采用體積圖拉普拉斯算子法[15](volumetric graph Laplacian,VGL)對每個關(guān)節(jié)周圍區(qū)域的網(wǎng)格進(jìn)行變形。VGL算法的最大優(yōu)點是可以保持曲面內(nèi)體積不變并避免曲面局部自交,另外其計算速度也較快。

VGL算法的主要原理是:對于待變形網(wǎng)格M,構(gòu)造一個填充網(wǎng)格內(nèi)部的體圖和一個覆蓋網(wǎng)格外側(cè)的體圖,用來防止曲面內(nèi)部體積的收縮和曲面的自交,然后通過類似泊松變形的傳播方法將控制曲線的變換顯式地傳播到感興趣區(qū)域, 最后通過線性變分求解變形后網(wǎng)絡(luò)坐標(biāo)。算法的實施過程描述為求解式(2)目標(biāo)函數(shù)的極值:

(2)

式中,LM為網(wǎng)格M的離散Laplace算子;vi為形變后網(wǎng)格的頂點坐標(biāo);εi為形變后網(wǎng)格的Laplace坐標(biāo);ui為控制點坐標(biāo);g為M圍成的體積圖;δi(1≤i≤T)為形變后體積圖的Laplace坐標(biāo);n為網(wǎng)格頂點總數(shù);m是控制點數(shù)目;T為體積圖頂點總數(shù);α、β均為約束強(qiáng)度調(diào)節(jié)參數(shù)。可見,目標(biāo)函數(shù)分為三個部分,分別刻畫對網(wǎng)格表面幾何細(xì)節(jié)、用戶指定約束和體圖細(xì)節(jié)的保持程度。變形效果如圖3所示。

圖3 實現(xiàn)的軟組織變形效果

(3)內(nèi)部器官變形

最后,根據(jù)骨骼和軟組織的變形來確定器官變形。變形采用的是近似剛體變化(as-rigid-as-possible,ARAP)算法,這是因為器官面元結(jié)構(gòu)復(fù)雜、變形精細(xì),與VGL相比ARAP[16]能更好地應(yīng)用于較為復(fù)雜的網(wǎng)格,且能實現(xiàn)更加真實、自然的變形效果。

設(shè)C至C′為剛體變化,則其變換過程中存在旋轉(zhuǎn)矩陣Ri如下:

(3)

ARAP變形算法的核心能量函數(shù)如式(4)所示,通過最小化該能量函數(shù)實現(xiàn)模型的盡可能剛性變形,此為形狀匹配問題的加權(quán)實例。

(4)

圖4 本文實現(xiàn)的器官變形效果

2 人體數(shù)字面元模型高速M(fèi)C計算技術(shù)

面元模型雖然能直接輸入Geant4等MC計算程序,但計算速度極慢,一般單次模擬耗時為數(shù)十小時,無法滿足實際應(yīng)用需求,其根本原因是每一步的粒子輸運(yùn)均需將粒子位置與所有面元位置比較以確定輸運(yùn)步長。本文建立了一種基于四面體剖分的計算技巧,大大提高了計算速度。具體來說,使用特定的四面體剖分算法將面元模型中由三角形面構(gòu)成的幾何空間分解為由無數(shù)單個四面體構(gòu)成的網(wǎng)格,這樣每步粒子輸運(yùn)只需將粒子位置與單個四面體進(jìn)行位置比較,從而大大縮短計算時間?？臻g四面體剖分算法主要有Delaunay算法[17]、八叉樹(octree)以及AFT算法等,而其中Delaunay三角化方法算法計算效率較高且剖分單元質(zhì)量好,因此本文使用Delaunay算法來實現(xiàn)面元模型的四面體剖分。算法流程如圖5(a)所示,四面體分割加速計算的原理如圖5(b),面元模型四面體切割后的效果如圖5(c)。

(a)基于Delaunay算法的四面體切割流程圖;(b)四面體切割提高計算效率的原理;(c)四面體切割實際效果。

3 與不可變形體素模型劑量計算的比較

為分析和驗證可變形人體模型在劑量計算上的優(yōu)勢,本文基于上述建立的面元模型劑量計算方法并利用開源MC程序Geant4計算了兩種姿態(tài)人體在兩種照射條件下的劑量,并與體素模型的計算結(jié)果進(jìn)行了比較。

計算條件如下:

(1)人體模型數(shù)據(jù)來源:面元模型為MRCP(來源ICRP 145),體素模型為ARCP(adult reference computational phantom,來源ICRP 110),均選擇其中的男性模型。

(2)人體模型姿態(tài):蹲姿和跪姿,如圖6所示,其中面元模型采用了第2節(jié)所述方法進(jìn)行姿態(tài)變形。

圖6 人體姿態(tài)與照射條件

(3)照射條件:661 keV、1 173 keV、1 332 keV三種能量的伽馬射線從正面和底面垂直均勻照射人體,三種射線數(shù)量比為1∶1∶1,以模擬作業(yè)環(huán)境中存在137Cs和60Co核素污染情況(兩種核素活度比為1∶1),模擬時伽馬射線面通量取3×109/cm2。

(4)統(tǒng)計方法:計算了ICRP 116號報告給出器官權(quán)重因子的所有器官的當(dāng)量劑量,以及眼晶體劑量和有效劑量。其他計算參數(shù)列于表1。計算結(jié)果如圖7所示,具體結(jié)果列于表2。

表1 面元、體素模型的若干參數(shù)比較

表2 面元模型計算結(jié)果與體素模型結(jié)果的差異

圖7 面元和體素模型在兩種姿態(tài)、兩種照射條件下計算結(jié)果的比較

從表2可見,前向照射條件下面元模型的有效劑量均低于體素模型結(jié)果,且器官劑量結(jié)果普遍顯著低于后者。具體來說:蹲姿前向照射條件下,面元模型的有效劑量比體素模型低11.0%,器官劑量較體素模型結(jié)果最大差異達(dá)-33.2%,其中性腺低17.7%,紅骨髓低23.0%;跪姿前向照射條件下,面元模型的有效劑量比體素模型低4.4%,器官劑量較體素模型結(jié)果最大差異達(dá)-20.1%,其中性腺低15.1%,紅骨髓低9.3%。這是由于相比于直立姿態(tài),蹲和跪時手臂、腿對前向入射的射線均有不同程度的遮擋,進(jìn)而減少了人體胸腹部及性腺等的受照劑量,可見在這些條件下使用傳統(tǒng)基于體素的劑量估算方法會高估實際劑量值。

底向照射時則正好相反。此時面元模型的有效劑量顯著高于體素模型結(jié)果,且器官劑量普遍大幅高于后者。蹲姿底向照射條件下,面元模型的有效劑量比體素模型高51.2%,器官劑量較體素模型結(jié)果最大差異達(dá)98.6%,其中性腺高63.1%,紅骨髓高78.3%;跪姿底向照射條件下,面元模型的有效劑量比體素模型高58.7%,器官劑量較體素模型結(jié)果最大差異達(dá)98.0%,其中性腺高59.9%,紅骨髓高75.9%。這是由于相比于直立姿態(tài),蹲和跪時雙腿叉開,使得有更多底向入射的射線直接照射人體胸腹部及性腺等,可見在這些條件下使用傳統(tǒng)基于體素的劑量估算方法會大大低估實際劑量值。

另外,面元模型的眼晶體劑量計算結(jié)果均高于體素模型,造成這樣差異的主要原因是體素模型使用體素表示人體器官組織,在眼晶體等較精細(xì)器官的表示上較為粗糙,空間分辨率較差,而面元模型則能更好地表征精細(xì)結(jié)構(gòu),這一差異體現(xiàn)在劑量結(jié)果上。

通過比較可見,目前基于體素模型的人體劑量計算技術(shù)難以表征人體姿態(tài)對受照劑量的影響,給出的有效劑量可低估超過50%,器官劑量低估可能接近100%,計算精度無法滿足日益發(fā)展的精準(zhǔn)防護(hù)概念的需求,而基于面元模型的劑量計算技術(shù)則可較好地解決不同人體姿態(tài)下的劑量計算問題,具備在事故劑量重建、應(yīng)急監(jiān)測等場景中付諸應(yīng)用的潛力。

4 結(jié)論

當(dāng)前基于體素模型的人體輻射劑量計算只能評價直立固定姿態(tài)下的人體劑量,無法評估人體姿態(tài)變化對受照劑量的影響,計算結(jié)果可能會有顯著偏差,限制了其在精準(zhǔn)防護(hù)實踐中的應(yīng)用。為解決可變?nèi)梭w姿態(tài)下劑量估算難題,發(fā)展新一代的人體輻射劑量計算技術(shù),本文在國內(nèi)率先建立了完整的人體數(shù)字面元模型變形算法和面元模型直接高速M(fèi)C計算方法,其中面元模型變形采用剛體旋轉(zhuǎn)矩陣、體積圖拉普拉斯算子和近似剛體變換三種算法,分別實現(xiàn)了骨骼、軟組織和內(nèi)部器官的變形;面元模型直接高速M(fèi)C計算則基于Delaunay四面體切割技術(shù)。上述算法、技術(shù)的構(gòu)建為劑量精準(zhǔn)計算的應(yīng)用奠定了基礎(chǔ)?；诖?本文比較了面元模型與體素模型的劑量計算結(jié)果,人體姿態(tài)為蹲、跪兩種,射線從前向和底向照射。計算結(jié)果顯示,前向照射條件下面元模型得到的有效劑量低于體素模型結(jié)果,器官劑量也普遍低于后者;而底向照射時正好相反,面元模型的結(jié)果均顯著高于后者。該對比試驗表明,目前基于體素模型的人體劑量計算技術(shù)在應(yīng)用中可能帶來較大偏差,有效劑量可能低估超過50%,器官劑量可能低估接近100%。綜上所述,本文建立了完整的基于可變形面元模型的新一代人體輻射劑量計算方法,未來有望在核事故劑量重建、核電廠檢修、核設(shè)施退役治理作業(yè)劑量預(yù)測及醫(yī)學(xué)介入治療醫(yī)護(hù)人員精準(zhǔn)防護(hù)等場景下實現(xiàn)人員劑量的精確評價。