基于可变形面元模型的新一代人体辐射剂量计算技术

赵 日,刘兆行,刘 娜,王仙祥,张 静,梁润成,刘 鑫,令狐仁静,戴雨玲

(1.中国辐射防护研究院,太原 030006;2.核药研发转化与精准防护山西省重点实验室,太原 030006)

人体辐射剂量计算是指利用计算机蒙特卡罗(monte carlo,MC)模拟电离辐射在环境和人体中的输运过程从而得到人体各处所受剂量值的技术。该仿真计算无需在人体布设实体剂量测量设备,且可在现场辐射作业开展前或完成后进行预测性或回顾性虚拟计算,同时,能给出全身任意处(包括内部器官、眼晶体等)的剂量结果,剂量分布空间分辨率可达毫米甚至微米级,因此,与传统佩戴剂量计的剂量监测相比更便捷、灵活且精细,近年来越来越受重视,已成为预测性和回溯性职业照射人员剂量评价及健康防护的关键技术。

人体辐射剂量计算中使用的人体数字模型是决定计算效果的关键。目前应用最广的是体素模型[1-3],使用体素(voxel)定义人体结构,类似于人体断层CT图像下使用像素(pixel)描述2维图像。ICRP 116号报告给出的人体外照射剂量转换系数就是基于体素模型计算得到的。然而限于体素模型的内在构造机制,模型不能进行姿态调整,只能保持直立姿态,因此严格来说,当前基于体素模型的剂量计算只适用于直立人体静态剂量评价,无法考虑人体姿态对受照剂量的影响。这一缺陷在均匀且强度不大的辐射场下尚能被接受,但在事故、应急等场景的高度非均匀强场条件下则会带来显著的剂量计算偏差。例如美国的Han等[4]报道,使用单一直立姿态人体模型进行事故剂量重建时,器官剂量最大低估达78%,有效剂量低估为19%;Yoem等[5]研究显示,非均匀辐射场中,不同姿态下人体红骨髓、肺、胃、结肠、乳腺、性腺等器官的剂量差异显著,其中性腺差异最大,极端情况下相差可超过两个数量级。因此,实际应用中亟需一种能够精确计算受照人体在各种姿态下所受剂量的新技术。

近期,ICRP 145号报告给出了一种全新的基于表面约束几何的参考人数字计算模型MRCPs (mesh-type reference computational phantoms),为更高精度的人体辐射剂量计算提供了可能。这种新形式的模型被称为BREP(boundary representative phantom)模型、Mesh模型或面元模型(本文使用“面元模型”代指)[6-7],如图1所示。得益于底层数学形式的灵活性,面元模型同时具备可姿态调整和高分辨率两方面优势,模型所有部分均可以进行自由移动、变形,极大的弥补了体素模型的缺陷,且模型空间分辨率没有下限,可保证组织器官轮廓的光滑性以及对微米级几何结构的精确描述,突破了体素模型毫米级分辨率下限和锯齿形轮廓的限制[8-10]。

图1 人体数字面元模型示意

目前,国外初步开展了基于面元模型的人体辐射剂量计算新技术研究。例如美国的Vazquez等[11]利用人体面元模型和动作捕捉系统对一起严重临界事故进行了剂量重建研究。使用动作捕捉系统来对人体动作进行记录,然后对面元模型变形重现这些动作,结果显示,新方法得出的剂量值与躯体症状相关性更好。Yeom等[12]计算了行走、坐、弯腰、跪、蹲共5种姿态在6种照射条件(即AP、PA、LLAT、RLAT、ROT、ISO)下的外照射剂量转换系数。国内仅有少数相关报道,这些零星研究对面元模型的变形过程比较粗糙,无法保证变形的物理合理性和几何光滑性,同时仍需将变形后的面元模型体素化成体素模型才能开展计算,使得面元模型的优势无法体现[13-14]。为此,本文开展了系统性研究,构建了精细化的面元模型变形方法,实现了面元模型的直接高速MC计算,为国内进一步开展、应用新一代高精度剂量计算技术奠定了基础。

1 人体数字面元模型变形方法

对面元模型进行变形是实现不同姿态人体剂量计算最关键的技术环节。但面元模型的数学结构高度复杂,且软组织(肌肉、皮肤、血管、淋巴组织)、内部器官的真实变形无法被实际观察,而目前已有的研究中的变形都太过粗糙,大多是依据关节的旋转而对软组织、器官进行直接旋转,与真实物理景象差别很大。所以鉴于骨骼、软组织、内部器官三部分的物理、几何特征的差异,分别进行变形,同时考虑到骨骼带动软组织、软组织带动器官变形的物理原理,因此按照骨骼、软组织、器官的顺序依次变形,且为降低变形过程的复杂度,假设变形是时序进行的,对一部分进行变形时,其余部分形状未开始改变。使用ICRP 145号报告给出的男性、女性参考人数字面元模型MRCP_AF和MRCP_AM,依次建立模型的骨骼、软组织和内部器官形变方法。具体算法如下。

(1)骨骼调整

读取待旋转骨骼的整个网格点数据,确定骨骼关节中心,以该处为旋转中心,确定旋转方向和旋转角度,然后计算旋转矩阵,基于该矩阵计算骨骼中各网格顶点的坐标旋转变换。

旋转时,任意点的坐标按式(1)三维空间中绕任意轴进行旋转的旋转矩阵进行变换:

(1)

待旋转的关节部位主要有15个,即肩关节(2个)、肘关节(2个)、腕掌关节(2个)、髋关节(2个)、膝关节(2个)和踝关节(2个)、脊柱关节(简化为3个,颈关节、腰关节、骶关节)。变形效果如图2所示。

图2 本文实现的骨骼旋转、变形效果

(2)软组织变形

接下来对各关节周围的软组织(肌肉、皮肤、血管)进行变形。采用体积图拉普拉斯算子法[15](volumetric graph Laplacian,VGL)对每个关节周围区域的网格进行变形。VGL算法的最大优点是可以保持曲面内体积不变并避免曲面局部自交,另外其计算速度也较快。

VGL算法的主要原理是:对于待变形网格M,构造一个填充网格内部的体图和一个覆盖网格外侧的体图,用来防止曲面内部体积的收缩和曲面的自交,然后通过类似泊松变形的传播方法将控制曲线的变换显式地传播到感兴趣区域, 最后通过线性变分求解变形后网络坐标。算法的实施过程描述为求解式(2)目标函数的极值:

(2)

式中,LM为网格M的离散Laplace算子;vi为形变后网格的顶点坐标;εi为形变后网格的Laplace坐标;ui为控制点坐标;g为M围成的体积图;δi(1≤i≤T)为形变后体积图的Laplace坐标;n为网格顶点总数;m是控制点数目;T为体积图顶点总数;α、β均为约束强度调节参数。可见,目标函数分为三个部分,分别刻画对网格表面几何细节、用户指定约束和体图细节的保持程度。变形效果如图3所示。

图3 实现的软组织变形效果

(3)内部器官变形

最后,根据骨骼和软组织的变形来确定器官变形。变形采用的是近似刚体变化(as-rigid-as-possible,ARAP)算法,这是因为器官面元结构复杂、变形精细,与VGL相比ARAP[16]能更好地应用于较为复杂的网格,且能实现更加真实、自然的变形效果。

设C至C′为刚体变化,则其变换过程中存在旋转矩阵Ri如下:

(3)

ARAP变形算法的核心能量函数如式(4)所示,通过最小化该能量函数实现模型的尽可能刚性变形,此为形状匹配问题的加权实例。

(4)

图4 本文实现的器官变形效果

2 人体数字面元模型高速MC计算技术

面元模型虽然能直接输入Geant4等MC计算程序,但计算速度极慢,一般单次模拟耗时为数十小时,无法满足实际应用需求,其根本原因是每一步的粒子输运均需将粒子位置与所有面元位置比较以确定输运步长。本文建立了一种基于四面体剖分的计算技巧,大大提高了计算速度。具体来说,使用特定的四面体剖分算法将面元模型中由三角形面构成的几何空间分解为由无数单个四面体构成的网格,这样每步粒子输运只需将粒子位置与单个四面体进行位置比较,从而大大缩短计算时间。空间四面体剖分算法主要有Delaunay算法[17]、八叉树(octree)以及AFT算法等,而其中Delaunay三角化方法算法计算效率较高且剖分单元质量好,因此本文使用Delaunay算法来实现面元模型的四面体剖分。算法流程如图5(a)所示,四面体分割加速计算的原理如图5(b),面元模型四面体切割后的效果如图5(c)。

(a)基于Delaunay算法的四面体切割流程图;(b)四面体切割提高计算效率的原理;(c)四面体切割实际效果。

3 与不可变形体素模型剂量计算的比较

为分析和验证可变形人体模型在剂量计算上的优势,本文基于上述建立的面元模型剂量计算方法并利用开源MC程序Geant4计算了两种姿态人体在两种照射条件下的剂量,并与体素模型的计算结果进行了比较。

计算条件如下:

(1)人体模型数据来源:面元模型为MRCP(来源ICRP 145),体素模型为ARCP(adult reference computational phantom,来源ICRP 110),均选择其中的男性模型。

(2)人体模型姿态:蹲姿和跪姿,如图6所示,其中面元模型采用了第2节所述方法进行姿态变形。

图6 人体姿态与照射条件

(3)照射条件:661 keV、1 173 keV、1 332 keV三种能量的伽马射线从正面和底面垂直均匀照射人体,三种射线数量比为1∶1∶1,以模拟作业环境中存在137Cs和60Co核素污染情况(两种核素活度比为1∶1),模拟时伽马射线面通量取3×109/cm2。

(4)统计方法:计算了ICRP 116号报告给出器官权重因子的所有器官的当量剂量,以及眼晶体剂量和有效剂量。其他计算参数列于表1。计算结果如图7所示,具体结果列于表2。

表1 面元、体素模型的若干参数比较

表2 面元模型计算结果与体素模型结果的差异

图7 面元和体素模型在两种姿态、两种照射条件下计算结果的比较

从表2可见,前向照射条件下面元模型的有效剂量均低于体素模型结果,且器官剂量结果普遍显著低于后者。具体来说:蹲姿前向照射条件下,面元模型的有效剂量比体素模型低11.0%,器官剂量较体素模型结果最大差异达-33.2%,其中性腺低17.7%,红骨髓低23.0%;跪姿前向照射条件下,面元模型的有效剂量比体素模型低4.4%,器官剂量较体素模型结果最大差异达-20.1%,其中性腺低15.1%,红骨髓低9.3%。这是由于相比于直立姿态,蹲和跪时手臂、腿对前向入射的射线均有不同程度的遮挡,进而减少了人体胸腹部及性腺等的受照剂量,可见在这些条件下使用传统基于体素的剂量估算方法会高估实际剂量值。

底向照射时则正好相反。此时面元模型的有效剂量显著高于体素模型结果,且器官剂量普遍大幅高于后者。蹲姿底向照射条件下,面元模型的有效剂量比体素模型高51.2%,器官剂量较体素模型结果最大差异达98.6%,其中性腺高63.1%,红骨髓高78.3%;跪姿底向照射条件下,面元模型的有效剂量比体素模型高58.7%,器官剂量较体素模型结果最大差异达98.0%,其中性腺高59.9%,红骨髓高75.9%。这是由于相比于直立姿态,蹲和跪时双腿叉开,使得有更多底向入射的射线直接照射人体胸腹部及性腺等,可见在这些条件下使用传统基于体素的剂量估算方法会大大低估实际剂量值。

另外,面元模型的眼晶体剂量计算结果均高于体素模型,造成这样差异的主要原因是体素模型使用体素表示人体器官组织,在眼晶体等较精细器官的表示上较为粗糙,空间分辨率较差,而面元模型则能更好地表征精细结构,这一差异体现在剂量结果上。

通过比较可见,目前基于体素模型的人体剂量计算技术难以表征人体姿态对受照剂量的影响,给出的有效剂量可低估超过50%,器官剂量低估可能接近100%,计算精度无法满足日益发展的精准防护概念的需求,而基于面元模型的剂量计算技术则可较好地解决不同人体姿态下的剂量计算问题,具备在事故剂量重建、应急监测等场景中付诸应用的潜力。

4 结论

当前基于体素模型的人体辐射剂量计算只能评价直立固定姿态下的人体剂量,无法评估人体姿态变化对受照剂量的影响,计算结果可能会有显著偏差,限制了其在精准防护实践中的应用。为解决可变人体姿态下剂量估算难题,发展新一代的人体辐射剂量计算技术,本文在国内率先建立了完整的人体数字面元模型变形算法和面元模型直接高速MC计算方法,其中面元模型变形采用刚体旋转矩阵、体积图拉普拉斯算子和近似刚体变换三种算法,分别实现了骨骼、软组织和内部器官的变形;面元模型直接高速MC计算则基于Delaunay四面体切割技术。上述算法、技术的构建为剂量精准计算的应用奠定了基础。基于此,本文比较了面元模型与体素模型的剂量计算结果,人体姿态为蹲、跪两种,射线从前向和底向照射。计算结果显示,前向照射条件下面元模型得到的有效剂量低于体素模型结果,器官剂量也普遍低于后者;而底向照射时正好相反,面元模型的结果均显著高于后者。该对比试验表明,目前基于体素模型的人体剂量计算技术在应用中可能带来较大偏差,有效剂量可能低估超过50%,器官剂量可能低估接近100%。综上所述,本文建立了完整的基于可变形面元模型的新一代人体辐射剂量计算方法,未来有望在核事故剂量重建、核电厂检修、核设施退役治理作业剂量预测及医学介入治疗医护人员精准防护等场景下实现人员剂量的精确评价。