带有松弛因子和惩罚项的最小二乘PET图像重建

刘 祎 桂志国

（中北大学电子测试技术国家重点实验室，太原 030051）

引言

正电子放射断层成像（positron emission tomography，PET）技术是一种放射性核素显像技术，是核医学成像技术和断层成像技术的最新发展。PET的最大优势是能定量评价人体组织的生理、生化功能，被誉为活体的分子断层图像［1］，也使得PET成为核医学诊断中的重要技术手段。

PET重建算法分为解析法和迭代法，经典解析法是滤波反投影算法（filtered back projection，FBP）［2］。FBP法简单，重建速度快，但FBP算法过分简化了PET的物理模型，要求有完整的投影数据，无法有效地抑制噪声［3］。而迭代重建算法适用于信息量相对不足的情况。经典的迭代算法是基于泊松模型的最大似然算法（maximum likelihood，ML）和基于高斯模型的加权最小二乘算法（weighted least squares，WLS）。但这两种算法都存在着由于噪声而导致的图像退化和收敛速度慢的问题，而且随着迭代次数的增加，图像质量下降，噪声非常严重［5］。为了加速迭代收敛，Hudson和Larkin提出了有序子集（ordered subset，OS）方法，应用于MLEM算法，取得了较好的图像效果。自此，OS算法不仅用于EM算法，更与其他的算法结合派生出新的“OS”型算法［6］，最典型的是OSLS算法。其他的加速收敛速度的算法还有Fessler提出的广义的空间更新期望值最大化（space alternating generalized expectation maximization algorithm，SAGE）算法［7-8］。虽然这些加速方法可以加快算法的收敛速度，但依然具有较差的噪声特性。为了克服噪声的影响，许多学者针对算法引入了约束条件和惩罚函数，以改善图像质量。

为解决噪声和收敛问题，本研究带有松弛因子和惩罚项的最小二乘算法。首先从最小二乘目标函数出发，为平滑图像和抑制噪声加入二次惩罚项，其次，为减少运算量而采用修改的逐次超松弛（successive over-relaxation，SOR）迭代法求解最优解［9］，并利用可变有序子集（MOS）对其进行加速收敛，获得基于可变有序子集的惩罚最小二乘重建算法，文中记为MOS-PLS+MSOR。仿真实验中，将MOS-PLS+MSOR算法与最小二乘（least squares，LS）、惩罚最小二乘（penalized least squares，PLS）等算法的重建结果作比较，实验证明，MOS-PLS+MSOR有更好的重建效果。

1 算法理论

在PET重建中，常采用的离散模型形式为

其中P=［p1，p2，…pN］T表示观测的投影数据，G为N×M维的系统矩阵，λ=［λ1，λ2，…λM］T表示图像向量。gij表示第j个像素被第i对探测器检测到的概率。

基于最小二乘估计的目标函数为

写成矩阵的形式为

因此，最小二乘估计为：

显然，对式（1）的求解相当于解下面的方程组：

由实际的系统模型可知G是大型的稀疏矩阵，直接求解相当困难，可以通过共轭梯度法等方法求解，这里使用SOR迭代算法［10］，因为同时SOR算法高频分量收敛速度比低频分量快，本质上也加快了整体的收敛速度［11］。而且OS方法虽然能加快重建速度，但OS重建的图像往往不能收敛，随着迭代次数的增加，噪声大大增加。解决OS发散的一个方法就是使用松弛参数［12］。

记A=GTG，B=GTP，则式（5）可以写成

很显然，A=（aij）∈RM×M是对称的非奇异矩阵，且满足aii≠0，i=1，2，…，N，由SOR迭代公式得：

式中，ω为松弛因子，取值范围为0＜ω＜2。将A=GTG、B=GTP带入式（7）得：

由于系数矩阵的特殊性，对于较小的图像，上式中的计算量已经非常大，不仅耗费的时间长，而且还要占用很大的内存空间。为了计算简便，本文考虑雅可比迭代的思想，将式（8）中的换成，即把式（8）中的后两项合并起来，这样就不用计算完后再计算，从而大大减少了迭代时间。从而得到简便的并行迭代公式，写成矩阵的形式：

很显然，式（9）属于同时迭代重建法，有利于克服图像对噪声的敏感性。PET图像重建实质是一个病态解的求解问题，式（9）的方法虽然有利于克服图像对噪声的敏感性，但仍然不能克服求解病态解得困难。正则化是求解病态问题的常用方法。因此，本文中引入正则项（惩罚项），许多学者已经提出很多惩罚函数［13］，一些惩罚函数目的在于保持边缘锐化的同时平滑灰度一致的区域。本文引入简单的二次平滑惩罚项［14］，如下所示：

其中Sj表示第j个像素的8邻近像素构成的像素集合，ωjk表示权重量，水平方向和垂直方向上ωjk=1，对角线方向上ωjk=1/。即ωjk的值同像素b和j之间的距离成反比。不考虑边缘效果，优化矩阵R可以表示为：

此时目标函数为

其中β是平滑参数，考虑二次惩罚项，式（9）的迭代公式变为

文中将这种方法记为PLS+MSOR方法。

2 基于可变有序子集的PLS算法

MSOR算法虽然减少了迭代过程中的计算量，使算法较之于传统的SOR使用了较少的时间，但并不能改变收敛速度。OS方法作为加速收敛速度的重要方法，已经得到广泛的应用。在理论上，随着子集划分个数增加，收敛速度增加，高频分量得到恢复。但并不是子集个数越多越好，因为一般要遵循“子集平衡”的条件，即每一个子集都含有相等的图像放射性参数信息，子集个数增加使噪声随迭代次数的增加逐步放大［15-17］。然而，若采用较小的子集个数，图像的高频信息就难以恢复，并且收敛速度也会变慢。

为了解决这个问题，采用可变有序子集的方法（MOS），先选择大的子集数，以恢复高频信息，再适当减小子集数，减轻噪声对重建图像的影响，然后再选择稍大的子集数，以恢复没有恢复的高频信息，如此循环往复，即加快了收敛速度，同时也减少了噪声的影响。MOS算法在迭代过程中改变子集的大小，克服了OS算法只采用固定子集造成的不能同时顾全噪声和高频信息的缺点。

3 实验方法

为验证MOS-PLS+MSOR算法的有效性，采用大小为128像素×128像素的Sheep-Logan头部模型（λ*）进行仿真，如图1所示。图像像素点取值范围为0到8，模拟的投影参数为N=128×128，即在0°～180°内均匀分布128个投影方向，每个方向有128个探测器对。系统矩阵G不考虑探测器的效率和衰减校正等因素，利用公式P*=G·λ*产生无噪声数据，并以P*为泊松变量噪声均值来生成所需要的观测数据P。为了公平比较，设置所有的算法均使用相同的均匀正值起始图像，初始值设为0.5，实验步骤如下：

步骤1：置k=0，给未知量λk赋初值；

步骤2：在t（t=1，...，T）个子集内顺序执行下面的迭代更新公式

步骤3：T个子集都更新，即完成一次迭代；

步骤4：设置迭代停止条件，满足条件则停止迭代，不满足条件则将更新后的像素值λk再重复执行步骤2。

图1 Sheep-Logan模型Fig.1 Sheep-Logan model

4 实验仿真

4.1 图像质量评价参数

为评价重建图像质量，需对重建结果进行定量描述。采用一些常见的质量评价参数定量分析图像的重建质量，分别是相关系数、信噪比、归一化均方误差，定义如下：

1）相关系数（correlation coefficient，CORR）

2）信噪比（signal to noise ratio，SNR）

信噪比越高表示重建图像质量越好。

3）归一化均方误差（root mean squared error，RMSE）

表征重建图像与原始图像的近似程度，RMSE越小则表明重建结果越接近原始图像。

4.2 几种重建算法的比较

将本研究的PLS+MSOR算法与LS，PLS算法做比较，图2分别给出了计量数为1×106时3种算法重建的结果。LS迭代21次，PLS迭代27次，PLS+MSOR迭代35次（ω=0.25β=0.002 12）。从图2的结果可以看出，在视觉效果上，3种算法中PLS+MSOR算法重建的图像较之于LS算法和PLS有较少的伪影，而且中间区域比较平滑，重建质量最好。同时，图3给出了3种算法的CORR曲线变化和SNR曲线变化。

表1中列出了各种算法的质量参数和迭代时间。从表中可以看出，由于PLS中含有惩罚项，因此迭代时间比LS的迭代时间长。但是PLS+SOR的迭代时间却比PLS的迭代时间短，甚至比LS的迭代时间还要短，足以看出PLS+MSOR算法在时间上的绝对优势。

图2 LS，PLS，PLS+MSOR算法重建图像。（a）LS算法；（b）PLS算法；（c）PLS+MSOR算法（计量数1×106）Fig.2 Reconstruction images of LS，PLS，PLS+MSOR.（a）LS；（b）PLS；（c）PLS+MSOR

图3 LS，PLS，PLS+MSOR算法的质量参数曲线。（a）CORR曲线；（b）SNR曲线Fig.3 Quality parameter curve of LS，PLS，PLS+MSOR.（a）CORR；（b）SNR

然而在实际的PET检测中，光子数量与摄入体内的放射核素的浓度有关，光子计数值不可能很大，如果要得到较大的计数值，就需要注射过量的放射性核素，这将给患者造成伤害。另外可获得的最大的计数值受探测器晶体的计数能力和后项处理电路的限制，因此在图4中给出了低计数情况下的重建结果。从图2和图3对比可以看出，低计数量情况下，LS和的PET图像重建质量远远不如高计数量情况下的质量，但是PLS+MSOR算法在低剂量情况下仍然有较好的重建效果。在表2中列出了计数量在1.26×105情况下3种算法重建图像的质量参数，可以看出，PLS+MSOR算法在低级数量时的重建效果要远好于LS算法和PLS算法。

图4 LS，PLS，PLS+MSOR算法重建图像。（a）LS算法；（b）PLS算法；（c）PLS+MSOR算法（计量数为1.26×105）Fig.4 Reconstruction images of LS，PLS，PLS+MSOR.（a）LS；（b）PLS；（c）PLS+MSOR

MOS-PLS+MSOR的子集个数序列为16，4，4，8，4，4，OS-LS和OS-PLS的子集数为4，并使相邻子集满足正交关系，重建结果如图5所示。表3中列出了3种算法的质量参数。相对于OS-LS和OS-PLS算法，本算法在有限的迭代中表现出较好的噪声特性，说明了OS-PLS不仅能够重建出质量较好的图像，而且具有收敛速度上的优势。为了更能说明本算法的可靠性，图6给出了3种算法的RMSE曲线和SNR曲线。

图5 OS-LS，OS-PLS和MOS-PLS+MSOR算法的重建图像。（a）OSLS算法；（b）OS-PLS算法；（c）MOS-PLS+MSOR算法；Fig.5 Reconstruction images of OS-LS，OS-PLS and MOS-PLS+MSOR.（a）OS-LS；（b）OS-PLS；（c）MOS-PLS+MSOR

表1 LS，PLS和PLS+SOR的质量参数（计量数为1×106）Tab.1 The quality parameter of LS PLS and PLS+SOR algorithm（The number of measurement：1×106）

表2 LS，PLS和PLS+SOR的质量参数（计量数为1.26×105）Tab.2 The quality parameter of LS PLS and PLS+SOR algorithm（The number of measurement：1.26×105）

图6 OS-LS，OS-PLS，MOS-PLS+MSOR算法的质量参数曲线。（a）RMSE曲线；（b）SNR曲线Fig.6 Quality parameter curve of OS-LS，OSPLS，MOS-PLS+MSOR.（a）curveof RMSE；（b）curve of SNR

表3 OS-LS，OS-PLS和MOS-PLS+MSOR的质量参数（计量数为1×106）Tab.3 The quality parameter of OS-LS，OS-PLS and MOS-PLS+SOR algorithm（The number of measurement：1×106）

4 结论与讨论

MOS-PLS+MSOR算法具有加速收敛和抑制噪声的作用，较之于传统的最小二乘算法有较好的重建效果。然而本研究中使用的二次先验惩罚项虽然比较简单且在理论上更加易于分析，但是对图像边缘区域造成过度平滑。同时因此使用更为合理的先验惩罚项，这也是进行下一步研究的内容。另外需要研究关于迭代过程中如何选取松弛因子的问题，但至今仍未能有效地解决，这也是下一步要研究的内容。

［1］张泽宝，医学影像物理学［M］.北京：人民卫生出版社，2000：189-191.

［2］Nam-Yong Lee，Yong Choi.A modified OSEM algorithm for PET reconstruction using wavelet processing［J］.Computer Methods and Programs in Biomedicine，2005，80：236-245.

［3］赵凌云.三维正电子发射断层成像理论及实现［D］.浙江：浙江大学，2001.

［4］Shepp SA，Vardi Y.Maximum likelihood restoration for emission tomography［J］.IEEE Trans Med Imageing，1982，1：113-122.

［5］Fessler J.A.，Hakan E.A paraboloidal surrogates algorithm for convergent penalized-likelihood emission image reconstruction［C］//Fessler J.A..IEEE Nuclear Science Symposium and Medical Imaging Conference Record.Toronto：IEEE，1998：1132-1135.

［6］周健.正电子发射断层的图像重建方法研究［D］.南京：东南大学，2006.

［7］Fessler JA，Hero AQ.Space-alternating generalized expectation maximization algorithm［J］.IEEE Trans on Signal Processing，1994，42（10）：2664-2676.

［8］Fessler JA，Hero AQ.Penalized maximum-likelihoodimage reconstruction using space-alternating generalized EM algorithms［J］.IEEE Trans Imag Proc，1995，4（10）：1417-1429.

［9］孔慧华，潘晋孝.序列子集联合代数重建技术［J］.CT理论与应用研究，2008，17（2）：40-45.

［10］封建胡，车刚明.数值分析原理［M］.北京：科学出版社，2001.

［11］Fessler JA.Penalized weighted least-squares image reconstruction for positron emission tomography［J］.IEEE Medical Imaging，1994，13（2）：290-300.

［12］颜建华.正电子发射断层图像重建算法研究［D］.武汉：华中科技大学，2007.

［13］Geman S，McClure D.E.，Statistical methods for tomographic image reconstruction［C］//Geman S.Processings 46th Session.ISI，Bulletin of the ISI.Tokyo：Legal Information Institute，1987：5-21.

［14］Sangtae A，Fessler JA.Globally convergent ordered subsets algorithms：Application to tomography［J］.IEEE Transactions on Medical Imaging，2003，22（5）：613-626.

［15］张书文，陈英茂，田嘉禾.OSEM子集及迭代次数对PET图像质量的影响［J］.中国医学影像学杂志，2003，11（3）：199-201.

［16］Liang Z，Jaszczak R，Creer K.On Bayesian image reconstruction from projections：uniform and nonuniform a priorisource information［J］.IEEE Transactions on Medical Imaging，1989，8（3）：227-235.

［17］Hudson HM，Larkin RS.Accelerated image reconstruction using ordered subsets of projection data［J］.IEEE Trans Medical Imaging，1994，13（4）：601-609.