4D人体辐射剂量计算技术研究

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

1（中国辐射防护研究院 太原 030006）

2（核药研发转化与精准防护山西省重点实验室 太原 030006）

人体辐射剂量计算是预测性和回溯性人员剂量 评价及健康防护的关键技术，目前主要基于人体数字模型和蒙特卡罗（Monte Carlo，MC）仿真模拟实现［1-4］。早期的人体数字模型为数学模型，简单、粗糙，无法满足当前精细化剂量计算需求；当下应用最广的是人体体素模型［5-9］，然而，限于体素模型的内在构造机制，模型不能进行姿态调整。因此，严格来说，基于体素模型的剂量计算只能用于直立姿态人体的剂量评价，计算得到的是静态三维（3D）结果，无法考虑人体姿态随时间变化的情况，不能计算4D剂量（空间3D及时间，共4D）。这一缺陷使其在实际应用中的适用性和准确性受限。例如，美国的Han等［10］报道，使用单一姿态人体模型进行事故剂量重建时，器官剂量最大低估达78%，有效剂量低估为19%；Yoem等［11］研究显示，非均匀辐射场中，不同姿态下人体红骨髓、肺、胃、结肠、乳腺、性腺等器官的剂量差异显著，其中性腺差异最大，极端情况下相差可超过两个数量级。

图1 基于人体数字模型的剂量计算示意及人体数学、体素模型Fig.1 Human phantom based dose calculation as well as mathematical and voxel phantom

可见，亟须发展4D化的人体辐射剂量计算技术，进一步提高剂量评价精度，为人员精准防护、辐射防护最优化等提供关键基础。

近期一种全新的基于表面约束几何的人体数字面元计算模型被开发出来，为4D剂量计算提供了可能［12-13］。面元模型同时具备可姿态调整和高分辨率两方面优势：模型所有部分均可以进行移动、变形，可人为改变模型姿态，实现与实际人体一致的模型姿态；同时，模型的空间分辨率没有下限，并保证这些结构的轮廓光滑性［14-16］。

图2 人体数字面元模型及Vazquez等[11]实现的动作捕捉引导模型变形效果(a) 人体数字面元模型，(b) 动作捕捉技术与剂量计算融合Fig.2 Mesh phantom and motion capture guided phantom deformation realized by Vazquez(a) Mesh phantom, (b) Merging of motion capture and dose calculation

基于面元模型，美国的Vazquez等［17］利用人体面元模型和动作捕捉系统对一起严重临界事故进行了剂量重建研究。他们使用动作捕捉系统来对人体动作进行记录，然后对面元模型变形重现这些动作，结果显示，新方法得出的剂量值与躯体症状相关性更好。Yoem等［18-19］计算了行走、坐、弯腰、跪、蹲共5种姿态在6种照射条件（即前向（Antero-posterior，AP）、背向（Postero-anterior，PA）、左侧（Left-lateral，LLAT）、右 侧（Right-lateral，RLAT）、旋 转（Rotational，ROT）、各向同性（Isotropic，ISO））下的外照射剂量转换系数，计算结果表明，直立姿态无法考虑手臂和大腿位置对躯干、骨盆附近器官剂量的影响，而4D剂量计算方法则可较好反映人体各处实际所受剂量。目前，国内仍然主要使用基于体素模型的3D静态剂量计算，仅有少数研究者初步开展了面元模型研究，4D剂量计算技术与国外有较大差距。为此，本文在国内率先构建了完整的4D人体辐射剂量计算方法流程。

1 人体数字面元模型姿态调整方法

本文使用了ICRP145号报告给出的参考人数字面元模型MRCP_AF和MRCP_AM，依次建立了模型的骨骼、软组织和内部器官的形变方法。

1.1 骨骼调整

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

旋转时，任意点的坐标按如下三维空间中绕任意轴进行旋转的旋转矩阵进行变换：

式中：Ri为绕轴旋转θ角的变换矩阵；ux、uy、uz分别为向量u⇀在x、y、z坐标轴的分量。

待旋转的关节部位主要有15个，即肩关节（2个）、肘关节（2个）、腕掌关节（2个）、髋关节（2个）、膝关节（2个）、踝关节（2个）、脊柱关节（简化为3个：颈关节、腰关节、骶关节）。

1.2 软组织变形

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

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

其中：LM是网格M的离散Laplace算子；εi是形变后网格的Laplace坐标；g'是g的子图；δi（1≤i≤N）是形变后子图g'中点的Laplace坐标。目标函数分为三个部分，分别刻画对网格表面几何细节、用户指定约束和体图细节的保持程度。

1.3 内部器官变形

最后，根据骨骼和软组织的变形来确定器官变形。变形采用的是近似刚体变化（As-Rigid-As-Possible，ARAP）算法［21］。这是因为器官面元结构复杂，变形精细，与VGL相比，ARAP能更好地应用于较为复杂的网格，且能实现了更加真实、自然的变形效果。

设C~C'为刚体变化，则其变换过程中存在旋转矩阵Ri如下：

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

式中：Ci和C'i分别表示变形前后模型顶点pi和p'i对应得变形单元；N(i)表示pi的1邻域点的索引，而pi和p'i分别表示pi和p'i的1邻域顶点；Ri表示Ci到C'i的最优旋转矩阵；wij是边eij=(pi，pj)的权重。

2 动作捕捉引导下的面元模型自动调整技术

借助动作捕捉设备捕获人体完整的时序性姿态，然后引导面元模型进行相应姿态调整，以高准确性在剂量模拟计算中重现人体实际动作情景，达到4D剂量计算目的。本文使用了两套光学动作捕捉系统：一套采用Nokov公司产品（图3（a）），由8个高精度光学摄像头、可穿戴反光球（Marker）与控制软件组成，适用于实验室条件下的精准捕捉，反光球具有高反光率，可将可见光反射回摄像头，以此确定每个摄像头视角下反光球的二维坐标，利用控制软件结合摄像头位置重建后可得到各反光球的三维坐标，可穿戴反光球固定在动作捕捉模特人体上的位置称为参考点，软件可以准确记录参考点三维位置的时间序列，据此建立动作捕捉骨架姿态的时间序列；另一套采用微软公司的Kinect v2.0系统（图3（b）），这套系统更为小型化和轻便，适用于现场条件。动作捕捉系统引导的模型变形效果如图4所示。

图3 本研究采用的动作捕捉系统 (a) Nokov系统，(b) Kinect系统Fig.3 Motion capture systems used in this study (a) Nokov system, (b) Kinect system

图4 动作捕捉引导的模型变形效果Fig.4 Phantom deformation guided by motion capture

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

面元模型虽然能直接输入Geant4等MC计算程序，但计算速度很慢，其根本原因是粒子输运中为确定步长，每次均需将粒子当前位置与所有面元的位置进行比较。为此，本文采用了四面体剖分技术将面元模型几何空间分解为四面体网格，粒子输运时每次只需与单个四面体进行位置比较，从而大大缩短计算时间。目前空间的四面体剖分算法主要有Delaunay算 法、八 叉 树（Octree）以 及AFT（Advancing Front Technique）法等，而其中Delaunay三角化方法算法［22］计算效率较高且剖分单元质量好，因此本文使用了Delaunay算法来实现对面元模型的四面体剖分。算法流程图如图5（a）所示。通过四面体剖分后，人体数字面元模型被划分为8.2×106个四面体，测试试验表明，在Intel i7个人电脑平台进行面元模型单次模拟计算耗时约5 min，较四面体切割前提速>50倍，能较好满足实际应用需求。

图5 四面体切割原理及效果(a) 基于Delaunay算法的四面体切割流程图，(b) 四面体切割提高计算效率的原理，(c) 四面体切割实际效果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 辐射剂量计算云应用系统开发

基于上述算法，本文建立了4D人体剂量计算云应用系统（图6），服务器端由CPU/GPU计算节点、存储服务器、网络服务器、交换机等构成，动作捕捉系统和VR虚拟操作系统通过网络上传人体时序动作参数至云端，多用户通过Web客户端访问服务器，上传场景建模参数，并获取剂量实时计算结果。该系统的搭建为4D剂量计算的便捷、高效应用提供了可能。

图6 本研究搭建的4D人体剂量计算云应用系统架构图Fig.6 The architecture diagram of 4D dose calculation cloud application system built in this study

5 应用实验

本文在秦山三期核电厂开展了现场应用实验。实验地点选择了其R/B四楼402厂房，场景选择了对该厂房内开展临时辐射屏蔽架设人员的剂量测量与计算，如图7（a）所示。实验开展前为作业人员佩戴了Hp（10）（胸口、性腺处）、Hp（3）（眼部附近）、Hp（0.07）（指端）剂量计，如图7（b）所示，剂量计均为中国辐射防护研究院自制的热释光剂量计（Thermo Luminescent Dosemeter，TLD），材质为LiF（Mg，Cu，P），使用前均进行了性能筛选和剂量刻度；同时，通过测量得到场景几何信息，包括现场环境中主要管道、设备的几何参数以及人体与物体相对位置信息等；另外，场景中的主要辐射源项为作业人员正前方的管道内的放射性核素，现场通过γ能谱仪测量能谱来确定放射性核素种类及表面活度，结果如表1所示。实验开始后，现场通过架设Kinect动作捕捉设备捕捉人体15个关节点的时序化坐标。现场数据采集结束后，将源核素信息、场景信息、动作捕捉数据等上传本文搭建的4D人体剂量计算云应用系统，系统自动对人体数字面元模型进行相应的时序化变形，然后利用四面体分割技术对变形后的模型进行分割，最后根据核素和场景信息进行蒙特卡罗模拟计算，得到最终的个人剂量结果。

表1 现场的放射性核素种类与表面活度值Table 1 Radionuclides and their surface activities in the workplace

图7 4D剂量现场应用 (a) 实验场景及动作捕捉示意，(b) 作业人员佩戴剂量计部位Fig.7 4D dose field application(a) Experiment scenario and motion capture demonstration, (b) Positions where dosimeters worn by workers

个人剂量计算值与实测值的比较如表2所示。结果显示：二者Hp（10）偏差小于10%，Hp（3）与Hp（0.07）的偏差小于15%。实验验证了在现场进行人体4D剂量计算的可行性与优势，检验了通过人体动作捕捉、面元模型变形、四面体切割、蒙特卡罗模拟计算得到个人剂量的可靠性。

表2 实测与计算剂量值的比较Table 2 Comparation of measured and calculated dose values

6 结语

针对当前人体辐射剂量计算只能进行直立固定姿态3D计算的缺陷，本文建立了完整的4D剂量计算方法，包括人体数字面元模型姿态调整方法、动作捕捉引导下的面元模型自动变形、面元模型四面体切割算法、高速MC计算等，并在现场进行了方法验证，结果表明，计算剂量值与实测值的偏差小于15%，能够满足实际应用需求。

未来，本文4D剂量计算技术可能的应用场景包括：1）核电厂运行或核设施退役中特定辐射作业的剂量预测或剂量重建。具体来说，现场使用动作捕捉设备记录人员完整时序姿态，或剂量重建时捕捉人员凭记忆还原的动作变化，然后进行人体数字面元模型的自动变形，并根据现场几何、物理条件及辐射源项信息进行高速模拟计算，给出人员全身剂量精准分布。2）医学介入治疗过程的医护人员的精准防护。医学介入治疗中，医护人员由于长时间近距离暴露于射线下，会受到较大辐射剂量，严重时可能会超过国家规定的年剂量限值。利用本文所建4D剂量计算技术，则可实时计算和显示医护人员全身各处剂量，从而实现精准防护。

作者贡献声明赵日：负责研究的提出及设计、数据的收集和整理、文章的起草和最终版本的修订；刘兆行、刘娜：负责算法开发、数据处理；张静、梁润成、王仙祥、戴雨玲、令狐仁静：负责实验测试及结果验证；刘立业：负责研究的提出及设计、文章审核、项目监督和管理。