CN110599562B - 基于多能量系统响应矩阵的放射源定位重建方法 - Google Patents

Abstract

本发明公开了一种基于多能量系统响应矩阵的放射源定位重建方法，将编码孔径成像系统所能探测的能量范围，划分为多个能量窗，在各能量窗下选取特征能量值，通过模拟得到各能量窗下特定能量值下的系统响应矩阵，在实测时，将测得数据按能量窗分别计数，最终形成多组与能量窗对应的投影数据，再采用最大似然期望最大化算法重建各能量窗的重建图像，将所有能量窗下的重建图像相加，得到最终的重建图像。本发明它建立了一种更加准确的系统响应模型，方法适用于单点源、多个相同能量的单点源、多个能量不同的单点源、形状源，并能对上述复杂辐射场景进行准确重建，能提供高质量、高保真度的重建图像，方法可移植性强，适用于不同几何结构、不同材料的设备。

Description

基于多能量系统响应矩阵的放射源定位重建方法

技术领域

本发明涉及一种放射源定位方法，尤其涉及一种基于多能量系统响应矩阵的放射源定位重建方法。

背景技术

在核设施退役、核废物处理、核安全与核安保等领域，相比于传统的辐射剂量仪、伽马谱仪等非成像设备，编码孔径成像探测系统可远距离实现放射源的空间分布成像定位，具备优良的角分辨率及探测灵敏度。其中，对放射源的数量、活度及分布进行高质量、高保真度的图像重建是为相关人员的人身安全提供有力保障的基础。

现有的编码孔径重建算法主要分为两大类：一种是线性重建方法，如维纳滤波、互相关法。其中互相关法利用投影图像与解码矩阵的循环卷积为δ函数这一理想数学特性，并且还具有：位于视野范围内不同方位角的点源，在探测器平面形成的投影具有唯一性的特点。线性方法的优势是重建速度快，但成像条件不理想时的重建图像存在重建伪影，且角分辨率依赖于编码板基本单元的开孔尺寸及编码板与探测器阵列之间的距离。另一种方法为非线性重建方法，如最大似然期望最大化法、压缩感知法等迭代重建算法。在核医学领域，放射源的种类、位置及发射的伽马射线能量事先已知，迭代重建算法可将射线对编码板的透射、晶体内的散射、统计涨落、编码板的自衰减效应、探测器的非线性等非理想效应进行建模及补偿，其中的系统响应矩阵尤为重要。现有的计算系统响应矩阵方法主要有解析计算模型和蒙特卡罗模型，其中蒙特卡罗模型获得的系统响应矩阵已被证实在核医学领域具有更优的空间分辨率及减少伪影的能力。

但不同于核医学领域，放射源的数量、种类及位置分布等在核设施退役、核安全与核安保等领域往往都是未知的，并且不同能量的伽马射线对编码板的穿透能力以及在探测器的响应不同。特别是对于多个不同能量的点源、形状源等复杂的辐射场景，由于不同能量放射源之间的响应函数具有差异，对于同时对多个放射源重建的情况，单个放射源的统计涨落会在其他放射源的投影数据中叠加一个不规则的噪声贡献项；而在对形状源进行重建时，形状源的部分区域也会将统计涨落叠加到该形状源其他区域的投影数据中一个不规则的噪声贡献项。

最大似然期望最大化法（Maximum Likelihood Expectation Maximization，MLEM），是基于泊松统计模型的最大贝叶斯后验概率的图像复原方法，MLEM 算法主要分为两步：1) 计算完全投影数据的似然函数在测定数据及当前参数估计值下的期望；2) 求出使期望最大化的参数估计值。

发明内容

本发明的目的就在于提供一种解决上述问题，能准确重建出高质量、高保真度的放射源空间分布图像的基于多能量系统响应矩阵的放射源定位重建方法。

为了实现上述目的，本发明采用的技术方案是这样的：一种基于多能量系统响应矩阵的放射源定位重建方法，采用编码孔径成像系统进行测量，所述编码孔径成像系统包括编码板和探测器组，其中探测器组为A×A个呈矩阵分布的探测器，探测器组经编码板探测的放射源平面为M×M个像素位置，包括以下步骤：

（1）将编码孔径成像系统所能探测的能量范围，划分为N个能量窗，为每个能量窗开辟一大小为A×A的内存空间，用于存储A×A个探测器的计数数据，且探测器计数初值为零；

在每个能量窗中选取一能量值，作为该能量窗下的特征能量值；

（2）构建一与编码孔径成像系统完全相同的模型，通过模型获取每个特征能量值条件下的系统响应矩阵，作为该能量窗的系统响应矩阵，每个系统响应矩阵的维度为A×A×M×M；

（3）在探测区域内设置编码孔径成像系统，进行一段时间数据采集，获得多个事例数据，每个事例数据包括能量沉积信息和位置信息，判断该事例数据的能量沉积信息位于哪个能量窗，再判断位置信息对应哪个探测器，将该能量窗对应的内存中间中，与该探测器对应的探测器计数加一；

最终形成N组与能量窗对应的投影数据；

（4）结合系统响应矩阵，采用最大似然期望最大化算法重建各能量窗的重建图像；

（5）将所有能量窗下的重建图像相加，得到最终的重建图像。

作为优选：步骤（1）中，在每个能量窗选取一能量值，作为该能量窗下的特征能量值，具体为，在每个能量窗中值处或靠近中值处选取。

作为优选：所述步骤（2）中模型为蒙特卡洛模型，获取每个特征能量值条件下的系统响应矩阵具体为：

（a1）采用蒙特卡罗模拟方法建立一与编码孔径成像系统完全相同的蒙特卡洛模型，所述蒙特卡洛模型包括编码板和探测器，编码板远离探测器的一端为放射源平面，并将放射源平面划分为M×M个像素；

（a2）在放射源平面内设置一发射伽马射线的单点源；

（a3）单点源以第一个特征能量值作为发射能量，历遍M×M个像素位置发射伽马射线，探测器组探测伽马射线，得到第一个特征能量值条件下的系统响应矩阵，将其作为第一个能量窗的系统响应矩阵，所述系统响应矩阵的维度为A×A×M×M；

（a4）按照（a3）的方法，得到其余所有能量窗的系统响应矩阵。

作为优选：所述步骤（4）具体为：

（41）建立以尺寸为M×M的全一矩阵，将其作为第一次的估计放射源图像；

（42）将估计放射源图像与第一个能量窗的系统响应矩阵进行正投影计算，得到正投影估计图像；

（43）将（42）中的正投影估计图像与步骤（3）的投影数据进行对比，得到误差项；

（44）将误差项与系统响应矩阵进行反投影计算得到矫正项；

（45）将矫正项归一化后，与上一次迭代输出的估计放射源图像相乘，得到本次迭代的估计放射源图像；

（46）重复步骤（42）-（45），直到收敛，得到第一个能量窗下的重建图像；

（47）根据步骤（41）-（46），得到所有能量窗下的重建图像。

作为优选：所述编码孔径成像系统的编码方式包括随机编码、非冗余编码、均匀冗余编码、改进均匀冗余编码、自支持编码。

本发明步骤（1）中：不仅将整个能谱划分为多个能量窗，还为每个能量窗开辟内存空间，这样每个能量窗对应的内存空间，都可以用来存储对应探测器计数，且探测器位置与内存空间内元素位置一一对应。

实际操作中，常见的伽马放射源能量一般在50keV~1.5MeV之间，根据实际的需要将能谱划分为N个能量窗，能量窗的数量越多，理论上的重建图像质量越好，但重建消耗的时间也与之成线性正比，需要折衷确定能量窗的数量。在取各个能量窗下的能量值作为该能量窗的特征能量值时，一般选取该窗的中值或与中值接近的放射源特征峰的能量值，当然也不仅限于此。

本发明步骤（2）中，是通过建模，模拟计算获取各个能量窗下的系统响应矩阵；

本发明步骤（3）是实际场景测量，将测量的数据统计到对应能量窗下的内存空间下，最终形成与能量窗对应的投影数据；

本发明步骤（4）（5）再采用最大似然期望最大化算法重建各能量窗的重建图像，并将所有能量窗下的重建图像相加，得到最终的重建图像。

与现有技术相比，本发明的优点在于：提出了一种新的放射源定位重建方法，本发明不同于现有技术中的整体测量方法只存储一组投影数据用于处理分析的方法，它建立了一种加准确的系统响应模型，方法适用于单点源、多个相同能量的单点源、多个能量不同的单点源、形状源，并能对上述复杂辐射场景进行准确重建，并能提供高质量、高保真度的重建图像，方法可移植性强，适用于不同几何结构、不同材料的设备。

附图说明

图1为本发明中能量窗划分示意图；

图2为编码孔径成像系统成像原理图；

图3为通过蒙特卡洛模型计算的多个系统响应矩阵图；

图4为本发明多组与能量窗对应的投影数据的示意图；

图5a为现有技术中单点源条件下采用互相关法进行重建的图像；

图5b为现有技术中单点源条件下采用单能量系统响应矩阵迭代重建的图像；

图5c为本发明单点源条件下重建的图像；

图6a为现有技术中多个相同能量的点源条件下采用互相关法进行重建的图像；

图6b为现有技术中单多个相同能量的点源条件下采用单能量系统响应矩阵迭代重建的图像；

图6c为本发明多个相同能量的点源条件下重建的图像；

图7a为现有技术中多个不同能量的点源条件下采用互相关法进行重建的图像；

图7b为现有技术中单多个不同能量的点源条件下采用单能量系统响应矩阵迭代重建的图像；

图7c为本发明多个不同能量的点源条件下重建的图像；

图8a为现有技术中两个环状源条件下采用互相关法进行重建的图像；

图8b为现有技术中两个环状源条件下采用单能量系统响应矩阵迭代重建的图像；

图8c为本发明两个环状源条件下重建的图像。

图中：1、放射源平面；2、点源；3、编码板；4、探测器组。

具体实施方式

下面将结合附图对本发明作进一步说明。

实施例1：参见图1到图4，一种基于多能量系统响应矩阵的放射源定位重建方法，采用编码孔径成像系统进行测量，所述编码孔径成像系统包括编码板3和探测器组4，其中探测器组4为A×A个呈矩阵分布的探测器，探测器组4经编码板3探测的放射源平面1为M×M个像素位置，包括以下步骤：

（1）将编码孔径成像系统所能探测的能量范围，划分为N个能量窗，为每个能量窗开辟一大小为A×A的内存空间，用于存储A×A个探测器的计数数据，且探测器计数初值为零；

在每个能量窗中选取一能量值，作为该能量窗下的特征能量值；本实施例中，我们在每个能量窗中值处或靠近中值处选取能量值。

（2）构建一与编码孔径成像系统完全相同的模型，通过模型获取每个特征能量值条件下的系统响应矩阵，作为该能量窗的系统响应矩阵，每个系统响应矩阵的维度为A×A×M×M；

本实施例中，模型为蒙特卡洛模型，获取每个特征能量值条件下的系统响应矩阵具体为：

（a1）采用蒙特卡罗模拟方法建立一与编码孔径成像系统完全相同的蒙特卡洛模型，所述蒙特卡洛模型包括编码板3和探测器，编码板3远离探测器的一端为放射源平面1，并将放射源平面1划分为M×M个像素；

（a2）在放射源平面1内设置一发射伽马射线的单点源2，也就是单放射源；

（a3）单点源2以第一个特征能量值作为发射能量，历遍M×M个像素位置发射伽马射线，探测器组4探测伽马射线，得到第一个特征能量值条件下的系统响应矩阵，将其作为第一个能量窗的系统响应矩阵，所述系统响应矩阵的维度为A×A×M×M；

（a4）按照（a3）的方法，得到其余所有能量窗的系统响应矩阵。

当然，不仅限于蒙特卡洛模型，也可以通过线驱动解析计算方法来获得所有能量窗的系统响应矩阵，此时，需要构造与编码孔径成像系统完全相同的几何构造模型，再通过线驱动解析计算得到。

（3）在探测区域内设置编码孔径成像系统，进行一段时间数据采集，获得多个事例数据，每个事例数据包括能量沉积信息和位置信息，判断该事例数据的能量沉积信息位于哪个能量窗，再判断位置信息对应哪个探测器，将该能量窗对应的内存中间中，与该探测器对应的探测器计数加一；

最终形成N组与能量窗对应的投影数据；

（4）结合系统响应矩阵，采用最大似然期望最大化算法重建各能量窗的重建图像；

（41）建立以尺寸为M×M的全一矩阵，将其作为第一次的估计放射源图像；

（42）将估计放射源图像与第一个能量窗的系统响应矩阵进行正投影计算，得到正投影估计图像；

（43）将（42）中的正投影估计图像与步骤（3）的投影数据进行对比，得到误差项；

（44）将误差项与系统响应矩阵进行反投影计算得到矫正项；

（45）将矫正项归一化后，与上一次迭代输出的估计放射源图像相乘，得到本次迭代的估计放射源图像；

（46）重复步骤（42）-（45），直到收敛，得到第一个能量窗下的重建图像；

（47）根据步骤（41）-（46），得到所有能量窗下的重建图像。

（5）将所有能量窗下的重建图像相加，得到最终的重建图像。

本实施例中：所述编码孔径成像系统的编码方式包括随机编码、非冗余编码、均匀冗余编码、改进均匀冗余编码、自支持编码。

实施例2：参见图1到图8，为了更好的说明本实施例的内存空间，我们具体举例如下：

（1）常见的伽马能量范围为50keV-1.5MeV,我们将编码孔径成像系统所能探测的能量范围也设置为50keV-1.5MeV，并将该范围内的整个能谱划分为4个能量窗，分别为能量窗1：50-200keV, 能量窗2 ：200-500keV, 能量窗3： 500-800keV, 能量窗4：800-1500keV；

我们为每个能量窗开辟一内存空间，假设探测器组4为3×3的阵列，也就是9个探测器，那么内存空间也为3×3的矩阵，且9个探测器与内存空间的九个元素一一对应；

在每个能量窗中，我们都要选择一个能量值，作为该能量窗下的特征能量值，在这里，我们分别取122keV、356keV、662keV、1173keV作为四个能量窗的特征能量值。

（2）通过模型获取每个特征能量值条件下的系统响应矩阵，我们选用蒙特卡罗模拟方法，并以122keV这个特征能量值举例说明：

假设，我们将放射源平面1划分为30×30个像素位置，再让单点源2以122keV的能量，历遍30×30个像素位置发射伽马射线，探测器组4探测伽马射线，得到122keV这个特征能量值条件下的系统响应矩阵，也就是第一个能量窗的系统响应矩阵，系统响应矩阵的维度为3×3×30×30； 同理，得到其余3个能量窗的系统响应矩阵待用；

（3）开始实地测量，进行一段时间数据采集，我们获得很多个事例数据，每个事例数据包括能量沉积信息和位置信息，判断该事例数据的能量沉积信息位于哪个能量窗，例如，若该事例数据能量沉积信息为170keV，位置信息为7，则对应能量窗1，对应第三行第一个探测器，那么就在能量窗1对应的内存空间中，第三排第一个位置的探测器计数值加一。若该事例数据能量沉积信息为320keV，位置信息为5，那么就在能量窗2对应的内存空间中，第二排第二个位置的探测器计数值加一。依次类推。这样，就能得到四组3*3的数据，与四个能量窗一一对应。这和传统方法中，只是整体测量，并未通过能量窗分类的方法，是完全不同的，具体可参见图4。

（4）结合系统响应矩阵，采用最大似然期望最大化算法重建各能量窗的重建图像；

（5）将所有能量窗下的重建图像相加，得到最终的重建图像。

为了验证方法的有效性，对单点源2、多个相同能量的单点源2、多个不同能量的单点源2、双环状源等复杂辐射场景进行了测试，并采用了互相关法、单能量系统响应矩阵以及本发明的方法进行了对比。

参见图5，对于单点源2，采用迭代法获取的重建点源2图像的全峰半高宽明显优于互相关法的重建结果，由于系统响应矩阵与成像探测系统对射线的响应不匹配的原因，采用单能系统响应矩阵的结果发生了弥散，而采用本发明的方法可以获得收敛于单个像素的重建热点。

参见图6，对于5个相同能量的点源2，互相关法的重建结果已无法分清有几个放射源，而是重建出了一个棒源；采用单能系统响应矩阵的方法重建出了6个热点，并且一些热点的位置存在明显错误；采用本发明的方法，5个相同能量的放射源可被准确重建。

参见图7，对于4个不同能量的点源2，采用互相关法仍然仅能重建出一个棒源；而采用单能量系统响应矩阵仅能重建出3个热点，并且有的热点位置发生了偏移；而采用本方法可准确重建出4个热点及位置。

参见图8，对于两个不同能量的环状源，互相关法的重建结果已经混叠，无法将两个环形分开，并且背景具有较强的伪影；采用单能系统响应矩阵重建的结果尽管背景伪影较少，但由于系统响应矩阵与实际成像探测系统对伽马射线的响应不匹配的原因，重建的两个环形仍然是模糊且无法分辨；而采用本发明提出的方法，两个环形可被准确重建出，并且内测环状源的空心是个标准的圆形，可证明本方法的重建结果具有更高的角分辨率。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于多能量系统响应矩阵的放射源定位重建方法，采用编码孔径成像系统进行测量，所述编码孔径成像系统包括编码板和探测器组，其中探测器组为A×A个呈矩阵分布的探测器，探测器组经编码板探测的放射源平面为M×M个像素位置，其特征在于：包括以下步骤：

（1）将编码孔径成像系统所能探测的能量范围，划分为N个能量窗，为每个能量窗开辟一大小为A×A的内存空间，用于存储A×A个探测器的计数数据，且探测器计数初值为零；

在每个能量窗中选取一能量值，作为该能量窗下的特征能量值；

（2）构建一与编码孔径成像系统完全相同的模型，通过模型获取每个特征能量值条件下的系统响应矩阵，作为该能量窗的系统响应矩阵，每个系统响应矩阵的维度为A×A×M×M；

（3）在探测区域内设置编码孔径成像系统，进行一段时间数据采集，获得多个事例数据，每个事例数据包括能量沉积信息和位置信息，判断该事例数据的能量沉积信息位于哪个能量窗，再判断位置信息对应哪个探测器，将该能量窗对应的内存空间中，与该探测器对应的探测器计数加一；

最终形成N组与能量窗对应的投影数据；

（4）结合系统响应矩阵，采用最大似然期望最大化算法重建各能量窗的重建图像；

（5）将所有能量窗下的重建图像相加，得到最终的重建图像；

所述步骤（2）中模型为蒙特卡洛模型，获取每个特征能量值条件下的系统响应矩阵具体为：

（a1）采用蒙特卡罗模拟方法建立一与编码孔径成像系统完全相同的蒙特卡洛模型，所述蒙特卡洛模型包括编码板和探测器，编码板远离探测器的一端为放射源平面，并将放射源平面划分为M×M个像素；

（a2）在放射源平面内设置一发射伽马射线的单点源；

（a3）单点源以第一个特征能量值作为发射能量，历遍M×M个像素位置发射伽马射线，探测器组探测伽马射线，得到第一个特征能量值条件下的系统响应矩阵，将其作为第一个能量窗的系统响应矩阵，所述系统响应矩阵的维度为A×A×M×M；

（a4）按照（a3）的方法，得到其余所有能量窗的系统响应矩阵。

2.根据权利要求1所述的基于多能量系统响应矩阵的放射源定位重建方法，其特征在于：步骤（1）中，在每个能量窗选取一能量值，作为该能量窗下的特征能量值，具体为，在每个能量窗中值处或靠近中值处选取。

3.根据权利要求1所述的基于多能量系统响应矩阵的放射源定位重建方法，其特征在于：所述步骤（4）具体为：

（41）建立以尺寸为M×M的全一矩阵，将其作为第一次的估计放射源图像；

（42）将估计放射源图像与第一个能量窗的系统响应矩阵进行正投影计算，得到正投影估计图像；

（43）将（42）中的正投影估计图像与步骤（3）的投影数据进行对比，得到误差项；

（44）将误差项与系统响应矩阵进行反投影计算得到矫正项；

（45）将矫正项归一化后，与上一次迭代输出的估计放射源图像相乘，得到本次迭代的估计放射源图像；

（46）重复步骤（42）-（45），直到收敛，得到第一个能量窗下的重建图像；

（47）根据步骤（41）-（46），得到所有能量窗下的重建图像。

4.根据权利要求1所述的基于多能量系统响应矩阵的放射源定位重建方法，其特征在于：所述编码孔径成像系统的编码方式包括随机编码、非冗余编码、均匀冗余编码、改进均匀冗余编码、自支持编码。

Images