趙 日 劉兆行 劉 娜 張 靜 梁潤(rùn)成 王仙祥 戴雨玲 令狐仁靜 劉立業(yè)

1（中國(guó)輻射防護(hù)研究院 太原 030006）

2（核藥研發(fā)轉(zhuǎn)化與精準(zhǔn)防護(hù)山西省重點(diǎn)實(shí)驗(yàn)室 太原 030006）

人體輻射劑量計(jì)算是預(yù)測(cè)性和回溯性人員劑量 評(píng)價(jià)及健康防護(hù)的關(guān)鍵技術(shù)，目前主要基于人體數(shù)字模型和蒙特卡羅（Monte Carlo，MC）仿真模擬實(shí)現(xiàn)［1-4］。早期的人體數(shù)字模型為數(shù)學(xué)模型，簡(jiǎn)單、粗糙，無法滿足當(dāng)前精細(xì)化劑量計(jì)算需求；當(dāng)下應(yīng)用最廣的是人體體素模型［5-9］，然而，限于體素模型的內(nèi)在構(gòu)造機(jī)制，模型不能進(jìn)行姿態(tài)調(diào)整。因此，嚴(yán)格來說，基于體素模型的劑量計(jì)算只能用于直立姿態(tài)人體的劑量評(píng)價(jià)，計(jì)算得到的是靜態(tài)三維（3D）結(jié)果，無法考慮人體姿態(tài)隨時(shí)間變化的情況，不能計(jì)算4D劑量（空間3D及時(shí)間，共4D）。這一缺陷使其在實(shí)際應(yīng)用中的適用性和準(zhǔn)確性受限。例如，美國(guó)的Han等［10］報(bào)道，使用單一姿態(tài)人體模型進(jìn)行事故劑量重建時(shí)，器官劑量最大低估達(dá)78%，有效劑量低估為19%；Yoem等［11］研究顯示，非均勻輻射場(chǎng)中，不同姿態(tài)下人體紅骨髓、肺、胃、結(jié)腸、乳腺、性腺等器官的劑量差異顯著，其中性腺差異最大，極端情況下相差可超過兩個(gè)數(shù)量級(jí)。

圖1 基于人體數(shù)字模型的劑量計(jì)算示意及人體數(shù)學(xué)、體素模型Fig.1 Human phantom based dose calculation as well as mathematical and voxel phantom

可見，亟須發(fā)展4D化的人體輻射劑量計(jì)算技術(shù)，進(jìn)一步提高劑量評(píng)價(jià)精度，為人員精準(zhǔn)防護(hù)、輻射防護(hù)最優(yōu)化等提供關(guān)鍵基礎(chǔ)。

近期一種全新的基于表面約束幾何的人體數(shù)字面元計(jì)算模型被開發(fā)出來，為4D劑量計(jì)算提供了可能［12-13］。面元模型同時(shí)具備可姿態(tài)調(diào)整和高分辨率兩方面優(yōu)勢(shì)：模型所有部分均可以進(jìn)行移動(dòng)、變形，可人為改變模型姿態(tài)，實(shí)現(xiàn)與實(shí)際人體一致的模型姿態(tài)；同時(shí)，模型的空間分辨率沒有下限，并保證這些結(jié)構(gòu)的輪廓光滑性［14-16］。

圖2 人體數(shù)字面元模型及Vazquez等[11]實(shí)現(xiàn)的動(dòng)作捕捉引導(dǎo)模型變形效果(a) 人體數(shù)字面元模型，(b) 動(dòng)作捕捉技術(shù)與劑量計(jì)算融合Fig.2 Mesh phantom and motion capture guided phantom deformation realized by Vazquez(a) Mesh phantom, (b) Merging of motion capture and dose calculation

基于面元模型，美國(guó)的Vazquez等［17］利用人體面元模型和動(dòng)作捕捉系統(tǒng)對(duì)一起嚴(yán)重臨界事故進(jìn)行了劑量重建研究。他們使用動(dòng)作捕捉系統(tǒng)來對(duì)人體動(dòng)作進(jìn)行記錄，然后對(duì)面元模型變形重現(xiàn)這些動(dòng)作，結(jié)果顯示，新方法得出的劑量值與軀體癥狀相關(guān)性更好。Yoem等［18-19］計(jì)算了行走、坐、彎腰、跪、蹲共5種姿態(tài)在6種照射條件（即前向（Antero-posterior，AP）、背向（Postero-anterior，PA）、左側(cè)（Left-lateral，LLAT）、右 側(cè)（Right-lateral，RLAT）、旋 轉(zhuǎn)（Rotational，ROT）、各向同性（Isotropic，ISO））下的外照射劑量轉(zhuǎn)換系數(shù)，計(jì)算結(jié)果表明，直立姿態(tài)無法考慮手臂和大腿位置對(duì)軀干、骨盆附近器官劑量的影響，而4D劑量計(jì)算方法則可較好反映人體各處實(shí)際所受劑量。目前，國(guó)內(nèi)仍然主要使用基于體素模型的3D靜態(tài)劑量計(jì)算，僅有少數(shù)研究者初步開展了面元模型研究，4D劑量計(jì)算技術(shù)與國(guó)外有較大差距。為此，本文在國(guó)內(nèi)率先構(gòu)建了完整的4D人體輻射劑量計(jì)算方法流程。

1 人體數(shù)字面元模型姿態(tài)調(diào)整方法

本文使用了ICRP145號(hào)報(bào)告給出的參考人數(shù)字面元模型MRCP_AF和MRCP_AM，依次建立了模型的骨骼、軟組織和內(nèi)部器官的形變方法。

1.1 骨骼調(diào)整

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

旋轉(zhuǎn)時(shí)，任意點(diǎn)的坐標(biāo)按如下三維空間中繞任意軸進(jìn)行旋轉(zhuǎn)的旋轉(zhuǎn)矩陣進(jìn)行變換：

式中：Ri為繞軸旋轉(zhuǎn)θ角的變換矩陣；ux、uy、uz分別為向量u?在x、y、z坐標(biāo)軸的分量。

待旋轉(zhuǎn)的關(guān)節(jié)部位主要有15個(gè)，即肩關(guān)節(jié)（2個(gè)）、肘關(guān)節(jié)（2個(gè)）、腕掌關(guān)節(jié)（2個(gè)）、髖關(guān)節(jié)（2個(gè)）、膝關(guān)節(jié)（2個(gè)）、踝關(guān)節(jié)（2個(gè)）、脊柱關(guān)節(jié)（簡(jiǎn)化為3個(gè)：頸關(guān)節(jié)、腰關(guān)節(jié)、骶關(guān)節(jié)）。

1.2 軟組織變形

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

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

其中：LM是網(wǎng)格M的離散Laplace算子；εi是形變后網(wǎng)格的Laplace坐標(biāo)；g'是g的子圖；δi（1≤i≤N）是形變后子圖g'中點(diǎn)的Laplace坐標(biāo)。目標(biāo)函數(shù)分為三個(gè)部分，分別刻畫對(duì)網(wǎng)格表面幾何細(xì)節(jié)、用戶指定約束和體圖細(xì)節(jié)的保持程度。

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

最后，根據(jù)骨骼和軟組織的變形來確定器官變形。變形采用的是近似剛體變化（As-Rigid-As-Possible，ARAP）算法［21］。這是因?yàn)槠鞴倜嬖Y(jié)構(gòu)復(fù)雜，變形精細(xì)，與VGL相比，ARAP能更好地應(yīng)用于較為復(fù)雜的網(wǎng)格，且能實(shí)現(xiàn)了更加真實(shí)、自然的變形效果。

設(shè)C~C'為剛體變化，則其變換過程中存在旋轉(zhuǎn)矩陣Ri如下：

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

式中：Ci和C'i分別表示變形前后模型頂點(diǎn)pi和p'i對(duì)應(yīng)得變形單元；N(i)表示pi的1鄰域點(diǎn)的索引，而pi和p'i分別表示pi和p'i的1鄰域頂點(diǎn)；Ri表示Ci到C'i的最優(yōu)旋轉(zhuǎn)矩陣；wij是邊eij=(pi，pj)的權(quán)重。

2 動(dòng)作捕捉引導(dǎo)下的面元模型自動(dòng)調(diào)整技術(shù)

借助動(dòng)作捕捉設(shè)備捕獲人體完整的時(shí)序性姿態(tài)，然后引導(dǎo)面元模型進(jìn)行相應(yīng)姿態(tài)調(diào)整，以高準(zhǔn)確性在劑量模擬計(jì)算中重現(xiàn)人體實(shí)際動(dòng)作情景，達(dá)到4D劑量計(jì)算目的。本文使用了兩套光學(xué)動(dòng)作捕捉系統(tǒng)：一套采用Nokov公司產(chǎn)品（圖3（a）），由8個(gè)高精度光學(xué)攝像頭、可穿戴反光球（Marker）與控制軟件組成，適用于實(shí)驗(yàn)室條件下的精準(zhǔn)捕捉，反光球具有高反光率，可將可見光反射回?cái)z像頭，以此確定每個(gè)攝像頭視角下反光球的二維坐標(biāo)，利用控制軟件結(jié)合攝像頭位置重建后可得到各反光球的三維坐標(biāo)，可穿戴反光球固定在動(dòng)作捕捉模特人體上的位置稱為參考點(diǎn)，軟件可以準(zhǔn)確記錄參考點(diǎn)三維位置的時(shí)間序列，據(jù)此建立動(dòng)作捕捉骨架姿態(tài)的時(shí)間序列；另一套采用微軟公司的Kinect v2.0系統(tǒng)（圖3（b）），這套系統(tǒng)更為小型化和輕便，適用于現(xiàn)場(chǎng)條件。動(dòng)作捕捉系統(tǒng)引導(dǎo)的模型變形效果如圖4所示。

圖3 本研究采用的動(dòng)作捕捉系統(tǒng) (a) Nokov系統(tǒng)，(b) Kinect系統(tǒng)Fig.3 Motion capture systems used in this study (a) Nokov system, (b) Kinect system

圖4 動(dòng)作捕捉引導(dǎo)的模型變形效果Fig.4 Phantom deformation guided by motion capture

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

面元模型雖然能直接輸入Geant4等MC計(jì)算程序，但計(jì)算速度很慢，其根本原因是粒子輸運(yùn)中為確定步長(zhǎng)，每次均需將粒子當(dāng)前位置與所有面元的位置進(jìn)行比較。為此，本文采用了四面體剖分技術(shù)將面元模型幾何空間分解為四面體網(wǎng)格，粒子輸運(yùn)時(shí)每次只需與單個(gè)四面體進(jìn)行位置比較，從而大大縮短計(jì)算時(shí)間。目前空間的四面體剖分算法主要有Delaunay算 法、八 叉 樹（Octree）以 及AFT（Advancing Front Technique）法等，而其中Delaunay三角化方法算法［22］計(jì)算效率較高且剖分單元質(zhì)量好，因此本文使用了Delaunay算法來實(shí)現(xiàn)對(duì)面元模型的四面體剖分。算法流程圖如圖5（a）所示。通過四面體剖分后，人體數(shù)字面元模型被劃分為8.2×106個(gè)四面體，測(cè)試試驗(yàn)表明，在Intel i7個(gè)人電腦平臺(tái)進(jìn)行面元模型單次模擬計(jì)算耗時(shí)約5 min，較四面體切割前提速>50倍，能較好滿足實(shí)際應(yīng)用需求。

圖5 四面體切割原理及效果(a) 基于Delaunay算法的四面體切割流程圖，(b) 四面體切割提高計(jì)算效率的原理，(c) 四面體切割實(shí)際效果Fig.5 Tetrahedralization principle and effect(a) Flow chart of Delaunay algorithm based tetrahedralization, (b) Principle of calculation efficiency improvement by tetrahedralization, (c) Real effect of tetrahedralization

4 輻射劑量計(jì)算云應(yīng)用系統(tǒng)開發(fā)

基于上述算法，本文建立了4D人體劑量計(jì)算云應(yīng)用系統(tǒng)（圖6），服務(wù)器端由CPU/GPU計(jì)算節(jié)點(diǎn)、存儲(chǔ)服務(wù)器、網(wǎng)絡(luò)服務(wù)器、交換機(jī)等構(gòu)成，動(dòng)作捕捉系統(tǒng)和VR虛擬操作系統(tǒng)通過網(wǎng)絡(luò)上傳人體時(shí)序動(dòng)作參數(shù)至云端，多用戶通過Web客戶端訪問服務(wù)器，上傳場(chǎng)景建模參數(shù)，并獲取劑量實(shí)時(shí)計(jì)算結(jié)果。該系統(tǒng)的搭建為4D劑量計(jì)算的便捷、高效應(yīng)用提供了可能。

圖6 本研究搭建的4D人體劑量計(jì)算云應(yīng)用系統(tǒng)架構(gòu)圖Fig.6 The architecture diagram of 4D dose calculation cloud application system built in this study

5 應(yīng)用實(shí)驗(yàn)

本文在秦山三期核電廠開展了現(xiàn)場(chǎng)應(yīng)用實(shí)驗(yàn)。實(shí)驗(yàn)地點(diǎn)選擇了其R/B四樓402廠房，場(chǎng)景選擇了對(duì)該廠房?jī)?nèi)開展臨時(shí)輻射屏蔽架設(shè)人員的劑量測(cè)量與計(jì)算，如圖7（a）所示。實(shí)驗(yàn)開展前為作業(yè)人員佩戴了Hp（10）（胸口、性腺處）、Hp（3）（眼部附近）、Hp（0.07）（指端）劑量計(jì)，如圖7（b）所示，劑量計(jì)均為中國(guó)輻射防護(hù)研究院自制的熱釋光劑量計(jì)（Thermo Luminescent Dosemeter，TLD），材質(zhì)為L(zhǎng)iF（Mg，Cu，P），使用前均進(jìn)行了性能篩選和劑量刻度；同時(shí)，通過測(cè)量得到場(chǎng)景幾何信息，包括現(xiàn)場(chǎng)環(huán)境中主要管道、設(shè)備的幾何參數(shù)以及人體與物體相對(duì)位置信息等；另外，場(chǎng)景中的主要輻射源項(xiàng)為作業(yè)人員正前方的管道內(nèi)的放射性核素，現(xiàn)場(chǎng)通過γ能譜儀測(cè)量能譜來確定放射性核素種類及表面活度，結(jié)果如表1所示。實(shí)驗(yàn)開始后，現(xiàn)場(chǎng)通過架設(shè)Kinect動(dòng)作捕捉設(shè)備捕捉人體15個(gè)關(guān)節(jié)點(diǎn)的時(shí)序化坐標(biāo)?，F(xiàn)場(chǎng)數(shù)據(jù)采集結(jié)束后，將源核素信息、場(chǎng)景信息、動(dòng)作捕捉數(shù)據(jù)等上傳本文搭建的4D人體劑量計(jì)算云應(yīng)用系統(tǒng)，系統(tǒng)自動(dòng)對(duì)人體數(shù)字面元模型進(jìn)行相應(yīng)的時(shí)序化變形，然后利用四面體分割技術(shù)對(duì)變形后的模型進(jìn)行分割，最后根據(jù)核素和場(chǎng)景信息進(jìn)行蒙特卡羅模擬計(jì)算，得到最終的個(gè)人劑量結(jié)果。

表1 現(xiàn)場(chǎng)的放射性核素種類與表面活度值Table 1 Radionuclides and their surface activities in the workplace

圖7 4D劑量現(xiàn)場(chǎng)應(yīng)用 (a) 實(shí)驗(yàn)場(chǎng)景及動(dòng)作捕捉示意，(b) 作業(yè)人員佩戴劑量計(jì)部位Fig.7 4D dose field application(a) Experiment scenario and motion capture demonstration, (b) Positions where dosimeters worn by workers

個(gè)人劑量計(jì)算值與實(shí)測(cè)值的比較如表2所示。結(jié)果顯示：二者Hp（10）偏差小于10%，Hp（3）與Hp（0.07）的偏差小于15%。實(shí)驗(yàn)驗(yàn)證了在現(xiàn)場(chǎng)進(jìn)行人體4D劑量計(jì)算的可行性與優(yōu)勢(shì)，檢驗(yàn)了通過人體動(dòng)作捕捉、面元模型變形、四面體切割、蒙特卡羅模擬計(jì)算得到個(gè)人劑量的可靠性。

表2 實(shí)測(cè)與計(jì)算劑量值的比較Table 2 Comparation of measured and calculated dose values

6 結(jié)語

針對(duì)當(dāng)前人體輻射劑量計(jì)算只能進(jìn)行直立固定姿態(tài)3D計(jì)算的缺陷，本文建立了完整的4D劑量計(jì)算方法，包括人體數(shù)字面元模型姿態(tài)調(diào)整方法、動(dòng)作捕捉引導(dǎo)下的面元模型自動(dòng)變形、面元模型四面體切割算法、高速M(fèi)C計(jì)算等，并在現(xiàn)場(chǎng)進(jìn)行了方法驗(yàn)證，結(jié)果表明，計(jì)算劑量值與實(shí)測(cè)值的偏差小于15%，能夠滿足實(shí)際應(yīng)用需求。

未來，本文4D劑量計(jì)算技術(shù)可能的應(yīng)用場(chǎng)景包括：1）核電廠運(yùn)行或核設(shè)施退役中特定輻射作業(yè)的劑量預(yù)測(cè)或劑量重建。具體來說，現(xiàn)場(chǎng)使用動(dòng)作捕捉設(shè)備記錄人員完整時(shí)序姿態(tài)，或劑量重建時(shí)捕捉人員憑記憶還原的動(dòng)作變化，然后進(jìn)行人體數(shù)字面元模型的自動(dòng)變形，并根據(jù)現(xiàn)場(chǎng)幾何、物理?xiàng)l件及輻射源項(xiàng)信息進(jìn)行高速模擬計(jì)算，給出人員全身劑量精準(zhǔn)分布。2）醫(yī)學(xué)介入治療過程的醫(yī)護(hù)人員的精準(zhǔn)防護(hù)。醫(yī)學(xué)介入治療中，醫(yī)護(hù)人員由于長(zhǎng)時(shí)間近距離暴露于射線下，會(huì)受到較大輻射劑量，嚴(yán)重時(shí)可能會(huì)超過國(guó)家規(guī)定的年劑量限值。利用本文所建4D劑量計(jì)算技術(shù)，則可實(shí)時(shí)計(jì)算和顯示醫(yī)護(hù)人員全身各處劑量，從而實(shí)現(xiàn)精準(zhǔn)防護(hù)。

作者貢獻(xiàn)聲明趙日：負(fù)責(zé)研究的提出及設(shè)計(jì)、數(shù)據(jù)的收集和整理、文章的起草和最終版本的修訂；劉兆行、劉娜：負(fù)責(zé)算法開發(fā)、數(shù)據(jù)處理；張靜、梁潤(rùn)成、王仙祥、戴雨玲、令狐仁靜：負(fù)責(zé)實(shí)驗(yàn)測(cè)試及結(jié)果驗(yàn)證；劉立業(yè)：負(fù)責(zé)研究的提出及設(shè)計(jì)、文章審核、項(xiàng)目監(jiān)督和管理。