核辐射全息定位成像技术

吕振雷， 胡一凡， 王瑞， 朱艺航， 俞彦成， 江年铭， 许天鹏， 乔智红， 马天予

(1.清华大学 工程物理系， 北京 100084； 2.北京永新医疗设备有限公司， 北京 102206； 3.国家核安保技术中心， 北京 102445)

随着核科学与核技术的不断发展和应用推广，放射性物质在医疗、科研、工业、能源等诸多领域发挥着越来越重要的作用。但放射性物质的广泛应用也必然伴随越来越多的核安全风险问题。当放射性物质发生丢失、被窃、泄漏等事故时，有可能会给人民生命健康及社会公共安全带来巨大的影响。根据国际原子能机构(International Atomic Energy Agency，IAEA)的调查报告，自2004年以来，世界范围内每年放射源丢失数目均达100次以上[1]。由于放射性物质所产生的核辐射难以被生物体直接感知到，且严重辐射照射将对生物体产生不可逆的严重损伤，又进一步增加了其风险的隐蔽性和严重性。

为了实现对放射性物质的监测和管理，当前常用的放射性物质监测设备主要包括：计数式设备和伽马相机2类。其中，计数式设备种类众多，但仅提供计数、强度、能量等信息，缺乏定位信息，应用过程中对放射性物质的监控效率低，准确性有限。基于编码板伽马相机[2-3]、康普顿相机[4-6]、全景式伽马相机[7]等成像定位技术，可以直观反映放射源方位，有效提高监测管理准确率，在我国海关口岸、环保、核电、安全等领域广泛应用。但目前的伽马相机设备均为二维方向成像技术，即单独静止的设备只能确定放射源相对于探测器所处的方向信息，而无法确定放射源在三维空间中的绝对位置。在一些复杂应用场景，例如海关出入境通道的密集人流、核电站复杂管道、货物堆放场等，放射源容易被不同物体遮挡，造成二维方向成像难以快速、准确地确定放射源的位置。

当前已有国外学术团队进行了伽马放射源的三维空间成像的研究工作[6-7]。这些工作采用的是康普顿相机，通过移动探测器采集运动路径上多个位置的放射源数据，对三维空间的放射源进行定位成像。这类技术适合于放射源静止不动，且非实时三维成像的应用中(例如货物堆放场、实验室内放射源定位等)，而对于放射源移动，需要实时成像的应用(例如海关出入境通道安检等)则无法适用。

为了解决上述问题，在基于文献[8-9]自主研发的全景式伽马相机技术基础上，本文提出开发核辐射全息定位成像系统，实现了在三维空间中对放射源进行绝对位置成像和三维全息场景匹配定位。

1 核辐射全息定位成像系统设计

核辐射全息定位成像系统主要由4个全景式辐射成像定位仪节点、4个深度相机及图像工作站构成，如图1所示。其中，4个全景式辐射成像定位仪通过数据组网和联合图像重建，实现对放射源分布的三维成像，深度相机实现可见光全息成像，图像工作站实现放射源三维图像重建、三维实时可见光图像生成以及辐射图像与可见光图像的融合显示等功能。

图1 核辐射全息定位成像系统示意Fig.1 Diagram of nuclear radiation holographic imaging system

1.1 系统坐标定义

以海关出入境口岸等人员流动量大的通道作为典型常规应用场景，设定基准系统应用和测试空间分方案。系统成像视野大小为5 m×5 m×2 m，建立三维笛卡尔坐标系，以重力线为Z轴，Z=0平面作为地面，系统坐标原点位于成像视野的底面的中心。根据成像定位精度要求，将成像视野的体素单元定为50 mm×50 mm×50 mm，整个系统像素阵列为100×100×40，如图2所示。

图2 核辐射全息定位成像系统的成像视野与几何坐标系Fig.2 FOV and geometric coordinates of nuclear radiation holographic imaging system

1.2 三维辐射成像定位系统设计

为了在整个成像视野中实现合理的放射源三维定位成像效果，所设计的辐射探测系统由4个辐射探测节点构成。同时兼顾系统在实际应用场景中的安装问题及考虑避免行人、物品、室内设施对射线的遮挡衰减，拟将4个探测节点吊装在距离地面2.5 m的高度上。

为了保证成像视野中成像性能的均匀性，基于空间对称性的基本原则，本研究共设计了4种不同的探测器排布方案，分别命名为design1～design4。同时，另外设计了由8个探测器组成的design5作为高成像性能的对照。这5种设计方案的示意图如图3所示。所有的探测器都位于Z=2.5 m的平面上，图中标明了每种方案下探测器的编号和对应的三维坐标。

1.3 单探测节点设计

辐射探测系统的单个探测节点是一台全景式辐射成像定位仪[5]，其主要核心部件为一个三维位置灵敏的闪烁晶体探测器，如图4所示。

闪烁晶体采用了GAGG(Ce)晶体，晶体阵列为16×16，晶体尺寸为4.05 mm×4.05 mm×20 mm，晶体间采用0.15 mm厚的BaSO4作为反射膜。SiPM板像素阵列亦为16×16，像素间距4.2 mm，从而闪烁晶体和SiPM一对一耦合。

伽马射线在闪烁晶体中沉积能量，转换为可见光子，双端耦合的SiPM阵列将可见光转换为电信号。SiPM阵列输出的多通道信号经过EXYT ASIC芯片[10]的位置加权及信号放大，输出位置信号(X，Y)、能量信号E和时间信号T。ASIC输出信号经过后端电子学数字化和预处理后传输到计算机。通过X，Y信号在预存储的泛场直方图分割表上可以查找到与伽马光子发生作用的晶体编号，如图5。

图3 用于对比三维成像性能的5种系统设计方案示意Fig.3 Illustration of five different system designs for system performance comparison

图4 单个探测节点外观及关键部件照片Fig.4 Appearance of single detection node and photos of key components

图5 探测器三维位置分辨示意Fig.5 Diagram of detector with 3D positioning capability

为了得到在晶体内作用的三维位置，除了晶体编号之外，还需要计算晶体内作用的具体位置。在本研究中，耦合在晶体两端的SiPM的能量信号E1和E2的相对大小来计算伽马事件在晶体内的深度作用位置(depth of interaction，DOI)，其计算公式[13]为：

式中：k和b是预先标定得到的探测器的固有系数，并固化于数据处理电路的FPGA芯片中。经实验测试，探测器模块的DOI分辨率约为3.4 mm，在本研究中将晶体轴向20 mm分为5层，每层4 mm。因此探测器的探测单元阵列为16×16×5。

每个探测节点均具备对放射源二维定向成像的能力，并完成了二维球坐标系下的系统传输矩阵的标定。标定具体步骤如下：

1)4π空间的图像域以二维球坐标系定义，其中θ角的取值范围为0°～ 180°，φ的取值范围为0° ～360°。θ和φ分别以10°为间隔，划分出一个36×19的测量网格。

2)在测量网格的交点处放置一个点源，逐点测量投影，每个点的测量时间在10 s左右，探测到的事件计数约为180万。将每个网格点处的投影以列向量形式组合在一起，获得对应粗测网格的传输矩阵。

3)基于粗测网格传输矩阵，通过样条插值运算得到θ和φ方向上采样间隔均为1°的精细网格传输矩阵(θ为0°～180°，φ为0°～359°)。

系统传输矩阵中的每一列对应于定义在图像域空间的某个探测器单元响应特性，因而称为探测器响应函数。以360×181的数字图像表达。图6展示了在高计数(180万个探测事件计数)情况下，对应于2个不同探测器单元的探测器响应函数，其物理意义是当这个探测器单元记录到一个伽马光子事件时，这个伽马光子来源于4π空间中每一个点的概率值。探测器划分为16×16×5个探测器单元(X-Y平面：16×16闪烁晶体矩阵；Z方向：5个DOI离散化单元)，因此一个完整的系统传输矩阵由16×16×5个同样格式的探测器响应函数组成。

图6 2个探测器单元对应的探测器响应函数结果Fig.6 Detector response function figures of two detector bins

1.4 深度相机系统设计

为了实现三维可见光全息图像，系统使用了4个微软的Kinect深度相机。为了保证三维系统5 m×5 m×2 m成像视野可被完整覆盖，并考虑减少行人、物品、室内设施的遮挡，故需将深度相机进行吊装，将4个深度相机对准视野中心，再考虑深度相机的成像视角等。最终确定了深度相机的安装坐标及外观如图7所示。每台深度相机可以独立采集和实时生成点云图像。

图7 深度相机成像系统设计及深度相机外观Fig.7 Design of depth camera imaging system and appearance of camera

采用基于实时多点云融合的三维重建技术，完成场景点云匹配。首先需要进行预标定，对不同视角捕捉到的点云使用Colored ICP方法[14]进行匹配，得到相对临时坐标系的变换矩阵，并将变换后的点云叠加为场景完整点云。同时在实际场景中放置多个标记物，通过标记物的临时坐标和物理坐标，计算得到坐标系间的变换矩阵。

由于目前的深度相机探测范围有限，且场景中部分物体如墙面和地面固定不变，因此可以对这部分进行离线预建模。离线建模有2种选择:1)使用单台深度相机在场景不同位置拍摄点云，然后使用Colored ICP方法进行融合;2)采用Open3D库提供的工具流程，手持单台深度相机连续移动拍摄场景，然后以TSDF方法进行离线建模。

点云实时渲染使用OpenGL，将每台相机捕捉到的点云、预建模点云和预标定得到的变换矩阵在GPU内变换至物理坐标，再通过渲染即可得到三维可见光全息图像。

2 三维辐射图像重建算法

本研究采用了基于统计估计原理和迭代修正策略的最大似然期望最大化(maximum likelihood expectation maximization, MLEM)算法[11-12]完成图像重建。其迭代重建公式为：

在MLEM算法中，精确的系统传输矩阵的获取是保证重建图像质量的最关键因素。一般而言，成像系统的系统传输矩阵的获取方法包括直接测量、几何计算、蒙特卡罗模拟等。直接测量法能够获得较为准确系统传输矩阵，但往往需要大量的实验测量工作；几何计算可以实现系统的快速计算，但因对系统物理过程、设备各种误差等问题的简化而造成计算的不精确；蒙特卡罗模拟可以较为精确地描述物理过程，但计算量较大且仍旧无法避免设备各种误差带来的问题。

对于本系统，由于成像视野大、图像体素数目多，难以通过实验直接测量整个系统传输矩阵。为了获得较为精确的系统传输矩阵，本文中采用在二维球坐标中测量单探测节点的定向系统传输矩阵，然后再计算整个系统的三维笛卡尔坐标下的系统传输矩阵的方法，具体过程描述如下。

首先通过实验的方式测量出4个探测节点在二维球坐标系下的二维系统传输矩阵S(θ,φ)。二维传输矩阵的标定实验[5]是将探测节点安装在一个二维旋转平台上，探测模块中心与旋转中心重合，放射源放置在距离探测器模块80 cm的位置上固定不动，利用平台的二维旋转，即可任意改变放射源与探测器模块之间的相对方向。以10°的采样间隔对整个二维球坐标系进行实验数据采集。再通过利用样条曲线插值的方法，将系统传输矩阵S(θ,φ)的间隔插值到1°间隔，提升系统传输矩阵的分辨精度。

然后针对每个探测节点，计算笛卡尔坐标系中每个体素单元相对于该探测节点在球坐标系中的坐标(θ,φ)和距离d，然后在该探测节点的二维系统传输矩阵中通过二维插值的方式计算出(θ,φ)的投影，再根据利用投影值与距离d之间平方反比关系，即可计算出笛卡尔坐标系下该体素单元对应的该探测节点的投影。按顺序遍历所有探测节点和所有的体素，即可以计算出完整的系统传输矩阵。

本文采用了GPU并行加速技术，基于本设计中系统传输矩阵尺寸4×16×16×5×100×100×40，迭代次数5 000次的计算量条件下，单次三维图像重建的时间约为5～10 s。

3 系统成像结果与精度分析

为了比较design1～design4这几种排布方案的优劣性，本研究选取了一个位于(-1.275 m，-1.275 m，0.975 m)的测试点进行了图像重建结果比较。图8给出5种设计方案对位于该位置处点源的成像结果，其中FWHM代表该重建图像在X，Y和Z方向上的半高宽平均值。

图8 5种设计方案对(-1.275 m，-1.275 m，0.975 m)位置处点源的图像重建结果Fig.8 Reconstructed images of the point source at (-1.275 m，-1.275 m，0.975 m) for five system designs

FWHM的值越小，代表图像分辨率越高。相应成像质量的排序为：design5>design2>design1>design4>design3。由此可以判断在由4个探测器组成的系统中，design2的表现是相对最好且最稳定的，且与探测器数目加倍的design5系统相比，图像分辨率差距不大。因此选用design2作为最终的设计方案。

为了实际评估本系统对放射源的三维定位精度，本文在测试场地完成了系统的实际搭建，如图9所示，并使用活度为6 mCi的57Co点源进行了实验测试。

图9 核辐射全息定位成像系统及测试场地Fig.9 The nuclear radiation holographic imaging system and testing space

将放射源放置在系统的成像视野中特定空间位置，评估成像结果；通过几何位置测量，确定放射源在空间中真实位置坐标；然后在系统重建三维图像中通过加权求重心的方法确定成像三维位置坐标，并计算2个坐标之间距离作为定位偏差。

实验测试了成像视野中的10个位置，对每个位置放射源采集了10 s的数据，每个探测节点的计数约为70 000个计数。计算每个位置的平均定位偏差结果，如表1所示。从表中数据可以计算出10个位置平均定位偏差为18.37 mm。图10中展示了表1中序号1的放射源重建图像的3个方向的剖面曲线，可以看到各个方向曲线均成高斯分布。

表1 系统定位偏差评估结果Table 1 Evaluation of system location bias

图10 三维重建点源(表1中序号1点源)各方向剖面曲线Fig.10 Profile curve of 3D reconstruction of point source(No.1 in Table 1)

将系统三维辐射图像与可见光图像融合定位成像结果如图11所示。融合图像清晰还原了实验场地的周边场景及放射源热点所在位置。

图11 系统三维全息成像示例Fig.11 Example of 3D holographic imaging

4 结论

1)本文完成了核辐射全息定位成像系统的设计、研发和测试。

2)系统实现了放射源的三维绝对位置定位和可见光三维全息成像，系统的平均定位偏差为18.37 mm。

3)本系统在海关、核电、公共安全等领域具有应用前景，对保障核安全、反核恐等具有重要意义。

目前，本文只对单点源定位成像进行了测试，其主要原因该系统的设计应用目标为海关出入境口岸等人员流动量大的通道，在这一场景下，一般监控范围内出现超过单个放射源的概率可以忽略；且旅客、行包内可能携带的放射性物质尺寸较小，可以合理地近似为点源。在其他放射源监控场景下，可能出现多点源甚至分布复杂的体源分布。

原则上只要多点源的间距大于系统的空间分辨率，本系统可以有效完成多点源定位成像。而对于分布复杂的体源成像，目前4节点系统容易出现采样完备性不足的情况，成像性能会受到限制，可以通过增加探测器节点数目或是增加探测器运动等方式来增加采样完备性，从而提高对体源的成像能力。未来将进一步对系统进行优化，并通过实验测试系统对多点源甚至是复杂体源的成像效果，从而明确系统未来的优化和升级方向。