编码成像MLEM互补算法的参数优化

2021-02-10侯英伟宋玉收孙世杰柳若琦黄丽萍胡力元刘辉兰

哈尔滨工程大学学报 2021年12期

侯英伟，宋玉收，孙世杰，柳若琦，黄丽萍，胡力元，刘辉兰

(1.哈尔滨工程大学核科学与技术学院，黑龙江哈尔滨 150001； 2.中国核动力研究设计院，四川成都 610213)

基于编码孔成像的伽马相机可以提供环境中伽马射线源的位置信息与能谱，快速对伽马源项进行定位与分析，因此编码成像技术已成功地应用于核材料监管和防扩散领域。MURA[1]编码具有开孔率大(约50%)、在有限大小内的自相关特性好、可实现中心反对称的方形排列等优点，被应用在各种成像技术中[2]，同时也是现在伽马编码孔成像最常用的编码方式。

当使用编码孔进行近场成像时，编码孔成像的数学模型证明[3]只有在远场条件下，实验设置具有理想相关特性成像关系才完美符合点扩散函数。在近场条件下，点源所放出的伽马射线不能看作是相互平行，因此源项在探测器上的投影具有不均匀性，经典反卷积解码后的重构图像会产生伪影[4]。最大似然期望最大化(maximum likelihood expectation maximization，MLEM)算法[5]可以进一步提高编码孔成像质量，但是当投影图像所含噪声较大时，迭代后噪声的影响也相应放大；另一方面为消除近场成像所产生的伪像与噪声，Accorsi等[6]提出了编码孔成像的双模测量法，采用硬件实现，不需要耗费大量计算时间，操作方式简便不会增加成像难度。为进一步消除伪像与噪声提高信噪比，李汉平等[7]提出基于互补码的改进MLEM算法，该算法在传统MLEM算法公式中加入编码板修正因子[8]，在迭代算法中进行互补编码板相减计算以达到消除噪声的目的。互补编码MLEM虽然能够有效去除伪像提高信噪比，但是互补码方法的修正因子与迭代次数是重要的成像参数，如果修正因子与迭代次数设置不正确，反而会引入噪声[9]，干扰源项位置的判断。并且如果源项在视野范围(field of view，FOV)内的不同位置，最佳的参数设置有所不同。

本文基于MLEM互补算法，利用Geant4对编码孔成像过程进行模拟，通过模拟不同源项位置的成像过程，研究互补编码MLEM成像法的重构图像结果质量与修正参数的关系，选取合适的修正因子大小，对不同源项位置所需迭代次数给出拟合经验公式并加以验证。

1 基于互补编码的MLEM成像算法

为减小由于投影不均匀性产生的重构伪影，MLEM被应用在编码孔成像中。MLEM算法是基于泊松统计Poisson-Likehood概率函数对数最大值的图像重建方法，MLEM算法为：

(1)

式中：fk(x，y)表示第k次迭代后放射源图像的估计值；y(x，y)表示探测器得到的投影数据；h(x，y)表示编码函数。迭代过程中的初始迭代估计值的取值有2种方式：1)f0(x,y)≡1；2)将相关解码得到的重建图像作为初始迭代估计值。

另一方面，成像时由于伪影的存在，将正编码模式和反编码模式对同一放射源进行编码孔成像，并采用相同的相关算法去解码重建图像[10]。将2幅重建图像的图像相加，就能消除编码孔近场成像时由于二阶项导致的伪影，从而增强重建图像的质量。

在使用MLEM迭代算法进行图像重建过程中，如果原始的编码板投影结果包含较强的噪声，随着迭代次数的增加噪声对图像的影响也会相应增大。为了抑制噪声对重建图像的影响，理论证明互补编码去噪方法对MLEM方法同样起作用。对迭代式(1)中的编码函数进行了修正，引入了编码板修正因子β。即：

h(x,y)=ha(x,y)+βhb(x,y)

(2)

式中：ha(x,y)是正编码函数；hb(x,y)是正编码函数旋转90°后得到的反编码函数，β为修正因子，取值范围(-0.8, 0)。式(2)代入式 (1)中可得：

(3)

由式 (2)可知，互补编码的MLEM成像算法是利用反编码修正后的编码函数，对迭代结果fk(x，y)进行投影后与原始投影相比较，并进行反投影运算。在“投影”与“反投影”过程中均利用反编码函数对噪声进行抑制，达到有效减少噪声的目的。

2 编码孔成像过程模拟

在利用互补编码的MLEM成像算法进行迭代成像过程中，需要选取合适的重构参数——迭代次数N与修正因子β，使得利用算法进行重构后的图像具有更好的图像质量。为了寻找合适的迭代参数，使用4.10.4版本Geant4进行蒙卡模拟，Geant4的PhysicsList采用QGSP_BERT，其包含所有的典型物理过程。

模拟采用137Cs伽马点状放射源，放出的伽马能量为662 keV。MURA编码板设置为19×19维扩展编码板，单质钨作为编码孔材料，钨的原子序数是74，密度为19.30 g/cm3，其对能量为662 keV的伽马射线的线性衰减系数约为1.7 cm-1，编码板几何厚度设置为1.5 cm。如图1所示为本次模拟的编码板具体设计形状，单位小孔大小为2 mm×2 mm。

图1 嵌套MURA (19×19)编码板开孔情况Fig.1 MURA (19×19) coding board

CsI晶体作为位置灵敏探测器，单位编码孔的采样率为4尺寸为47.5 mm×47.5 mm×50 mm，即探测器在XOY平面被平均分为76×76个像素，通过Geant4的灵敏探测器功能获取伽马射线入射至CsI晶体的位置信息，并划归至对应的CsI像素前端面中心位置，以获得离散化的投影图像。

将伽马源所在位置设置为坐标原点，伽马点源、编码板、探测器沿Z轴依次放置，源项所在位置即探测系统的FOV为坐标系的XOY平面，原点为探测器与编码板的中心投影位置。伽马源距编码板中心80 cm，编码板中心与探测器前端面的距离为20 cm。投影成像的放大倍数为1.25倍，此编码板探测系统的FOV为19 cm×19 cm。图2为编码孔成像模拟的几何设置示意图。

图2 编码孔成像几何设置示意Fig.2 Diagrammatic of coding aperture imaging

3 成像结果分析与优化方法

以极坐标观察伽马源所在XOY平面，改变伽马源在成像平面的位置，将源项分别设置于极角θ为0°、15°、30°、45°处，距原点距离R即源项所在极径分别设置为50、80、160、200、240和290 mm处分别进行成像模拟。使用基于互补编码板的MLEM方法进行源项位置重构，测量不同源项位置合适迭代次数N与修正因子β的关系。

3.1 成像结果质量评价参数

为更好地判断重构参数的选取是否正确，引入峰噪比(peak contrast to noise ratio，PCNR)与源项重构强度的半高全宽(full width half maximum，FWHM)作为成像质量的评价参数。峰噪比PCNR为：

(4)

式中：S为重构图像中计数强度的最大值；B为图像背景均值；σ为图像背景标准差。

PCNR越大证明图像的峰值相对于背景噪声的对比度越大，在成像结果源项定位位置正确的前提下PCNR间接表示了图像的成像质量。在进行互补编码板的MLEM迭代计算过程中，为选取合适的迭代次数N与修正因子β，选取峰噪比PCNR与FWHM作为判断标准。如图3为源项位置位于(141.4, 141.4)时，不同修正因子β设置条件下FWHM与PCNR随迭代次数增加的变化情况，不同形状点代表不同的β取值，每个图像的右侧端点为首次迭代开始点。

由图3可知无论选取多大的修正因子，随着迭代次数的增加，PCNR逐渐变大与FWHM逐渐变小也就是重构图像的质量逐渐提高。

图3 不同修正因子β设置条件下FWHM与PCNR随迭代次数增加的变化情况Fig.3 With different correction factors β, the change of FWHM and PCNR with the increase of iteration number

随着迭代次数的增加，图像质量的提升幅度尤其是FWHM变小的幅度越来越小，过多的迭代次数会极大降低计算效率消耗计算性能，因此计算过程中不推荐使用过多的迭代次数。在某些修正因子下当重构图像质量提升至一定大小后，如图3中β>0.4情况，图4所示继续进行迭代计算将导致图像出现紊乱，使得重构的源项位置出现偏差，并且图像质量降低，甚至图像变形。图4为源项位置位于(141.4, 141.4)，β取值为-0.5时不同迭代次数的成像结果。

图4 不同迭代次数的成像结果Fig.4 Imaging results of different iterations

为避免成像结果的质量出现震荡及提高计算效率，在模拟伽马源成像的过程中，设置当成像结果第一次达到PCNR>30时，便达到了合适的迭代次数。由于本文算研究内容的均是单点源，PCNR是对整个图像的评价参数，而FWHM是对重构的源项位置分辨的评价参数，所以没有将FWHM参数作为迭代次数选取的限制条件。

3.2 迭代修正参数的计算方法

为研究修正参数与源项所在极径的关系，将同一极角的数据进行对比分析。图5为源项所在位置的极径不同时，最佳迭代次数N随着修正因子β取值的变化情况。图5(a)与(b)分别为源项所在位置极角0°与30°，不同形状点为不同源项所在位置的极径大小。由图5可知，当采用的修正因子绝对值太小时，如图中当β>-0.3时，噪声的抑制效果不明显，因此需要更多的迭代次数才能获得较好的成像结果。需要注意的是，过多的迭代次数不仅有可能放大噪声对成像的影响，而且增加了成像过程的不稳定性。当源项所在位置极径较大位于FOV边缘时，源项在探测器上的投影更不均匀，相应重构后的伪影也较强。因此当修正参数绝对值较低时较弱的去噪效果不能完全消除伪影，无论经过多少次迭代计算也无法获得合适的重构图像结果。如图5中当R=290 mm，β>-0.3时，无法达到PCNR>30质量要求，没有有效的数据点。

图5 源项所在位置的极径不同时，最佳迭代次数N随着修正因子β的变化情况Fig.5 When the polar diameter of the source is different, the optimal number of iterations N changes with the correction factor β

当修正因子的绝对值很大时，伪影消除效果明显，所需的迭代次数也相应减少。但如图3所示过强的伪影消除效果，导致迭代次数增加至一定数量反而使得成像质量出现震荡，甚至无法成像定位。且每次迭代结果相差较大不易控制迭代次数，不易寻找到普适性的迭代次数。通过对比成像结果可知，由于伪像主要由源项在探测器上投影不均匀造成的，因此相比于极径，角度对参数选取的影响较小。

总的来说，相比于源项所在极角，极径对于迭代参数的选取影响更大；不应选取绝对值过小的修正系数，会降低迭代计算效率以及无法重构出较大的极径位置的源项；也不应选取绝对值过大的修正系数，使得所需迭代次数过小而不易控制成像质量并且导致图像质量变差。

综上所述，过高的的或过低的修正因子选取，都会相应的产生图像重构问题，因此综合考虑，在当前模拟设置下，选取-0.5作为固定的修正因子，此数值与文献[8]中提供的修正因子计算公式的计算结果接近。图6为修正因子为-0.5时，迭代次数随源项所在极径的变化情况，不同形状的数据点为不同的源项所在极角大小。趋势再次证明，在确定的修正因子条件下，迭代次数与源项极角关系很小，且随着极径的增加所需迭代次数也相应增加。

利用3次多项式对源项所在极径R与最佳迭代次数N进行拟合，拟合结果如图6所示，关系式为：

图6 最佳迭代次数N与源项所在极径的关系Fig.6 The relationship between the optimal number of iterations N and the polar diameter of the source

N=[-1.40+0.15R-1.1×10-3R2+

2.77×10-6R3]

(5)

式中R为源项位置的极径大小，需要注意的是由于迭代次数不能是小数，因此需要将迭代次数的计算结果四舍五入取整。并且由于编码板以及探测器是都是正方形，并且编码板形状是以正方形对角线为对称轴的轴对称图形，因此通过模拟0～45°的源项位置可以代表大多数的源项位置情况。利用此关系式，即可完成针对此设置下的不同成像位置的迭代参数选取。

3.3 优化方法验证与对比

为验证这一方法的准确性，在不改变其他实验模拟设置的条件下，模拟(40 mm, 40 mm)与(200 mm, 200 mm)位置的伽马点源的成像结果，其位置所在极径分别为56.57、282.84 mm，若β取值为-0.5，通过式 (5)计算得到的理想迭代次数为4次与16次，计算如图7所示，图7(a)与(b)为源项位置(40 mm, 40 mm)时，未进行迭代计算与迭代4次之后的成像结果，图7(c)为对应利用传统δ解码[1]成像结果，图7(d)与(e)为源项位置(200 mm, 200 mm)时，未进行迭代计算与迭代16次之后的成像结果，图7(f)为对应利用传统δ解码[11-12]成像结果。评价参数如表1所示。

表1 成像结果评价参数PCNR对比表Table 1 Comparison table of evaluation parameter PCNR of imaging results

图7 利用计算得到的迭代次数所重构的图像结果Fig.7 The result of image reconstruction using the number of iterations calculated

由图7可知，经过计算公式所得次数进行迭代后的位置分辨率达到1.2 mm，并且如表 1所示，虽然未经迭代的成像结果PCNR参数均小于传统δ解码方法，但经过特定次数的迭代之后，成像结果的PCNR参数明显提高，成像质量相应有效增加。结果说明利用此方法计算得到的迭代次数与修正因子的大小，针对本次模拟条件，进行基于MLEM互补算法成像可以有效提高成像质量。