CN104992422A - 基于离散剪切波正则化的低剂量ct图像统计重建方法 - Google Patents

Abstract

本发明公开了一种基于离散剪切波正则化的低剂量CT图像统计重建方法。本发明在数据保真项加入符合数据统计特性的系数加权，以降低噪声对重建结果的影响，并利用增广拉格朗日方法将图像在剪切波域可以稀疏表示的先验信息作为正则化项加入目标函数。本发明的正则化系数是根据误差自适应变化的，进一步提高了算法的鲁棒性。重建图像在辐射剂量降低到FBP算法的10％甚至更低时仍然能够得到清晰保留结构细节的重建图像。本发明适用于少视角和低管电流扫描两种情况下的低剂量CT图像统计迭代重建方法。

Description

基于离散剪切波正则化的低剂量CT图像统计重建方法

技术领域

本发明涉及低剂量CT的图像统计重建技术，更具体的说，是面向少视角投影或低X-射线管电流投影的CT图像统计重建方法。

背景技术

由于CT技术具有快速、准确、无创伤、无痛苦等特点，被越来越多地应用于临床诊断。但是X射线在透射过程中，会将部分能量转移到人体，引起身体损伤甚至致癌。近年来，由于CT的广泛使用，辐射风险也越来越受到关注，低剂量CT问题逐渐成为研究热点。

常见的降低CT剂量的方法主要有两种，一是降低X射线管电流强度来降低每个视角下的曝光剂量。当前临床上为获得高质量的CT重建图像，所用CT设备X射线管电流强度都较高，一般临床诊断X射线管电流强度一般超过200mA。降低管电流强度能降低单次扫描的辐射剂量，但会使投影数据信噪比下降，噪声强度呈指数倍增长。二是在不改变单个视角曝光剂量的前提下，减少扫描视角的数量。临床上一次CT扫描通常需要上千视角，扫描视角的减少可以降低一次扫描的总辐射剂量，但会使投影数据的数目远远小于待建CT图像的像素数目，使图像重建问题变成一个欠定问题。

对于第一种降低剂量的策略，通常对投影数据进行降噪预处理或使用统计迭代方法(Statistical Iterative Reconstruction，SIR)，在迭代重建方法中加入投影数据的统计特性。如部分学者采用平稳高斯模型对投影数据中的噪声进行建模，并采用最大后验概率估计(maximum a posteriori，MAP)方法对代价函数寻优。还有学者则采用具有边界保持特性的先验模型对重建图像进行建模，并同样采用MAP估计方法进行迭代图像重建。同时，还有学者指出低剂量CT的投影数据在滤除少量孤立噪声点后服从高斯分布，因此提出基于分割的自适应滤波方式，提高图像重建质量。也有学者在充分考虑了二维投影图像中相邻像素的差异性后，提出了采用各向异性二次惩罚加权最小二乘算法实现低剂量CT图像去噪，较好保持了重建图像的边缘细节信息。

线性代数中，将不适定问题转化为适定问题称为正则化。因此可以通过引入待建图像的先验知识来作为正则化项，使高度病态的不完备投影数据获得稳定而准确的重建。其中最常见的是全变分(Tatal Variation，TV)正则化迭代方法。2006年已经证明一个信号如果是可稀疏表示的，则可以利用全变分最小化作为正则化项，从少量测量数据中精确重建该信号。全变分最小化的稀疏角度CT图像重建方法每次迭代都由凸集投影和梯度下降两项组成。在此基础上，很多文献提出了一些改进的全变分正则化迭代CT图像重建方法。但是当投影数据含噪声时，利用上述方法得到的重建图像中出现了较为严重的噪声干扰和条形伪影，且部分细节被噪声掩盖，不能得到满意的重建图像。有学者提出构造离散梯度变换和离散差分变换的伪逆形式，将软阈值滤波的思想应用到基于TV最小化的CT重建中，加快了方法的收敛速度。但在重建临床诊断CT图像重建时使用TV最小化这一先验仍然存在一些问题。首先，TV最小化只是对数据保真项所重建出的图像的离散梯度变换之和的大小进行约束，属于一种具有较好的边缘保存能力的光滑性约束，但是表达图像特征的能力有限。因此，对其进行稀疏性约束时容易导致在消除噪声的情况下损失信号。其次，TV最小化约束建立在图像分片光滑的基础上，实际中的CT图像并不能精确满足这一条件，基于TV最小化约束的不完备数据CT重建中，重建结果常常存在图像边缘不清晰、表达细小结构能力差、高噪声下产生块状伪影等问题。因此，寻找更合适的图像稀疏变换和稀疏表达方式，是进一步提高不完备投影数据重建质量的关键。

图像的稀疏度对图像重建质量有着重要的影响，CT图像中常常包含了大量的曲线和局部信息，所以需要寻找更易于有效地捕捉局部图像特征的变换域。剪切波(Shearlet)是采用具有合成膨胀的仿射系统构造的一种多维函数，继承了曲波变换和轮廓波变换平移不变性和方向选择性的优点，通过对基函数缩放、剪切和平移等仿射变换，生成具有不同特征的剪切波基函数，是多尺度几何分析的最新发展，可以对图像的细节信息给出最优的表示性能。目前在图像去噪、压缩、融合等领域取得了很好的效果。对于二维信号，剪切波变换不仅可以检测到所有的奇异点，而且还可以自适应跟踪奇异曲线的方向，且随着尺度参数变化，可精确描述函数的奇异性特征，从而获得更好的图像稀疏表示能力。因此，基于离散剪切波变换的稀疏表示约束比基于离散梯度变换的TV最小化约束具有更好的欠完备投影CT图像重建性能。

发明内容

本发明所要解决的技术问题是，克服上述现有技术中存在的不足，提供一种基于离散剪切波正则化的低剂量CT图像统计重建方法。

将待建断层的X射线线性衰减系数图像离散化为N个像素，用列向量f＝[f₀,f₁，…,f_N-1]^T表示，探测器测得的投影数据为p(l,θ)，其中l表示像素距离原点的长度，θ表示扫描视角。假设共采集到M个投影数据，用列向量p＝[p₀,p₁,…,p_M-1]^T表示，数据采集系统的几何结构用矩阵A＝{a_ij}∈R^M×N表示，称为系统矩阵，其中a_ij表示第j个像素对第i条投影路径的几何贡献。传统迭代重建方法就是在已知A和p的基础上通过迭代方法由Af＝p求解出f，但是当辐射剂量降低时，由于不满足完备条件或者信噪比过低，重建效果差。

剪切波变换具有很强的二维图像稀疏表示及去噪性能，基于剪切波变换的低剂量CT图像重建方法就是在求解Af＝p时加入待建图像的剪切波稀疏表示作为正则化项，以缩小解空间，提高重建质量。该方法可以表示为求解如下问题

$\begin{array}{ccccc}\underset{\mathrm{\α}}{min}\left|\right|\mathrm{\α}|{|}_0& s.t.& f=\mathrm{\ψ}\mathrm{\α}& and& Af=p\end{array}---\left(1\right)$

其中，α为稀疏系数向量，其中只有少数非零元素，||||₀表示l₀范数。ψ表示离散剪切波变换。

式(1)是一个约束最优化问题，它所对应的增广拉格朗日方程为

$\mathrm{\Φ}(f,\mathrm{\α})=\underset{\mathrm{\α}}{min}\left|\right|\mathrm{\α}|{|}_0-\mathrm{\λ}(f-\mathrm{\ψ}\mathrm{\α})+\frac{\mathrm{\μ}}{2}||f-\mathrm{\ψ}\mathrm{\α}|{|}_2^2-\mathrm{\β}(Af-p)+\frac{\mathrm{\η}}{2}\left|\right|Af-p|{|}_2^2---\left(2\right)$

其中，参数μ、η、λ和β是增广拉格朗日系数，μ和η根据经验取固定值，而向量λ和β分别随稀疏表示误差以及重建图像的投影误差的减小而增大。这样可以在保证重建结果的同时提高收敛速度。

当X射线管电流强度下降，投影数据信噪比降低时，考虑投影数据的统计特性可以提高方法的抗噪性能。根据低管电流条件下真实投影数据的实验分析结果可知，经系统校准及对数变换后的低剂量CT投影数据可近似认为是理想投影数据经加性噪声n污染的结果，其中n近似服从空间非平稳高斯分布，其各个投影数据上的噪声均值为0，方差与各投影数据自身的统计均值呈现非线性解析关系，其解析式可用如下公式描述

${\mathrm{\σ}}_{p_i}^2=f_i\×\exp ({\stackrel{\‾}{p}}_i/\mathrm{\η})---\left(3\right)$

其中，为第i个探测器(detector)上获得的数据均值，为对应的噪声方差；η和f_i为与所扫描物体无关的参数，其值完全由所用的CT设备的配置决定，f_i为对应于第i个探测器的参数，其值可以事先从重复扫描获得的投影数据中计算得到，参数η的作用是作为一个描述系统校准过程的尺度系数，而f_i则主要用于反映弓形滤波器(bowtie filter)的作用。

由式(3)可知，噪声强度随着投影数值的增长呈指数增长，为了降低噪声对重建效果的影响，在式(2)的目标函数中加入统计加权

$w_i=\frac{1}{{\mathrm{\σ}}_{p_i}}---\left(4\right)$

其中，p_i为第i个探测器上获得的投影数据。那么式(2)的目标函数变为

$\mathrm{\Φ}(f,\mathrm{\α})=\underset{\mathrm{\α},f}{\min }\[\begin{array}{c}\left|\right|\mathrm{\α}|{|}_0-{\mathrm{\λ}}^T(f-\mathrm{\ψ}\mathrm{\α})+\frac{\mathrm{\μ}}{2}||f-\mathrm{\ψ}\mathrm{\α}|{|}_2^2\\ -\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\β}}_i{w}_i({\[Af\]}_i-p_i)+\frac{\mathrm{\η}}{2}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i\left|\right|{\[Af\]}_i-p_i|{|}_2^2\end{array}\]---\left(5\right)$

其中，μ和η取固定值，而向量λ和β分别随稀疏表示误差以及重建图像的投影误差的减小而增大，β_i为β中的第i个元素，1≤i≤M。λ和β分别如式(6)和式(7)所示

λ^t＝λ^t-1-μ(f^t-1-ψα^t-1)   (6)

${\mathrm{\β}}_i^t={\mathrm{\β}}_i^{t-1}-{\mathrm{\ηw}}_i({\[{\mathrm{Af}}^{t-1}\]}_i-p_i)---\left(7\right)$

其中，t≥1表示迭代次数。

该方法的实现过程是使用交替最小化的方法来求解这两个变量；首先对f^t-1(t≥1)进行稀疏表示，目标函数为

${\mathrm{\α}}^t=\underset{{\mathrm{\α}}^t}{min}\[\left|\right|{\mathrm{\α}}^t|{|}_0-\mathrm{\λ}(f^{t-1}-{\mathrm{\ψ\α}}^t)+\frac{\mathrm{\μ}}{2}||f^{t-1}-{\mathrm{\ψ\α}}^t|{|}_2^2\]---\left(8\right)$

然后固定式(8)求得的剪切波系数α^t，更新重建图像f^t，这时的优化目标为

$f^t=\underset{f^t}{\min }\[\left|\right|{\mathrm{\α}}^t|{|}_1-\mathrm{\λ}(f^t-{\mathrm{\ψ\α}}^t)+\frac{\mathrm{\μ}}{2}||f^t-{\mathrm{\ψ\α}}^t|{|}_2^2-\mathrm{\β}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i({\[{\mathrm{Af}}^t\]}_i-p_i)+\frac{\mathrm{\η}}{2}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}w{i}_{}\left|\right|{\[{\mathrm{Af}}^t\]}_i-p_i|{|}_2^2\]---\left(9\right)$

本方法设置的终止目标为 $\left|\right|{\mathrm{Af}}^t-p|{|}_2^2\≤{10}^{-4}$ 或迭代次数达到100次，当达到上述目标之一时得到最终重建CT图像。

本发明主要是针对低剂量CT投影数据欠完备或低信噪比的特点，根据投影数据所特有的空间非平稳噪声特性，设计一种基于离散剪切波正则化的低剂量CT统计迭代重建方法。本方法在辐射剂量降低到不足原来10％时仍能准确重建CT图像。其有益效果主要体现在：

(1)新颖性：首次在数据保真项加入符合数据统计特性的系数加权，以降低噪声对重建结果的影响，并将待建图像在剪切波域可以稀疏表示作为先验信息，利用增广拉格朗日方法将此先验信息作为正则化项加入目标函数，缩小了解空间，提出一种离散剪切波变换正则化的低剂量CT图像统计迭代重建方法。

(2)有效性：通过实验证明了重建图像在投影数据远远不满足完备性条件，或投影数据信噪比急剧下降的情况下本方法能够重建出高质量图像。在辐射剂量降低到滤波反投影(Filtered-back Projection，FBP)方法的10％甚至更低时仍然能够得到清晰保留结构细节的重建图像。

(3)实用性：简单可行，可以用于少视角扫描和低管电流扫描两种情况下的低剂量CT重建中。

本发明该方法通过在数据保真项加入符合投影数据统计特性的加权系数，将待建图像在剪切波域的稀疏表示作为先验信息加入到目标函数中，并设置自适应的正则化参数以增加方法的鲁棒性。该方法利用增广拉格朗日方法将此先验信息作为正则化项加入目标函数缩小解空间，实现了基于离散剪切波变换的CT图像统计迭代重建方法的设计，适用于少视角和低管电流扫描两种情况下的低剂量CT图像统计迭代重建方法。

附图说明

图1是本方法验证过程使用的仿真模体；

图2a是FBP方法对100个视角的投影数据进行重建得到的CT图像；

图2b是SIRT方法对100个视角的投影数据进行重建得到的CT图像；

图2c是TV正则化方法对100个视角的投影数据进行重建得到的CT图像；

图2d是本发明方法对100个视角的投影数据进行重建得到的CT图像；

图3a是图2a得到的CT图像左肺部细节；

图3b是图2b得到的CT图像左肺部细节；

图3c是图2c得到的CT图像左肺部细节；

图3d是图2d得到的CT图像左肺部细节；

图4a是FBP方法对低X-射线管电流投影重建图像；

图4b是SIRT方法对低X-射线管电流投影重建图像；

图4c是TV正则化方法对低X-射线管电流投影重建图像；

图4d是本发明方法对低X-射线管电流投影重建图像；

图5a是图4a得到的CT图像左肺部细节；

图5b是图4b得到的CT图像左肺部细节；

图5c是图4c得到的CT图像左肺部细节；

图5d是图4d得到的CT图像左肺部细节。

具体实施方式

该方法的具体过程如下：

步骤1初始化CT图像为投影数据经FBP方法重建的图像f⁰，^μ和^η分别取固定值2⁶和2¹²。

步骤2对f^t-1(t≥1)进行稀疏表示。假设，这时求剪切波系数的目标函数为

${\mathrm{\α}}^t=\underset{{\mathrm{\α}}^t}{min}\[\left|\right|{\mathrm{\α}}^t|{|}_0-\mathrm{\λ}(f^{t-1}-{\mathrm{\ψ\α}}^t)+\frac{\mathrm{\μ}}{2}||f^{t-1}-{\mathrm{\ψ\α}}^t|{|}_2^2\]---\left(8\right)$

由于求解l₀范数问题是NP-hard问题，计算量太大，计算过程中可由l₁范数代替，目标函数变为

${\mathrm{\α}}^t=\underset{{\mathrm{\α}}^t}{min}\[\left|\right|{\mathrm{\α}}^t|{|}_1-\mathrm{\λ}(f^{t-1}-{\mathrm{\ψ\α}}^t)+\frac{\mathrm{\μ}}{2}||f^{t-1}-{\mathrm{\ψ\α}}^t|{|}_2^2\]---\left(9\right)$

步骤3固定剪切波系数α^t，更新重建图像f^t。这时的优化目标为

$f^t=\underset{f^t}{\min }\[\left|\right|{\mathrm{\α}}^t|{|}_1-\mathrm{\λ}(f^t-{\mathrm{\ψ\α}}^t)+\frac{\mathrm{\μ}}{2}||f^t-{\mathrm{\ψ\α}}^t|{|}_2^2-\mathrm{\β}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i({\[{\mathrm{Af}}^t\]}_i-p_i)+\frac{\mathrm{\η}}{2}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}w{i}_{}\left|\right|{\[{\mathrm{Af}}^t\]}_i-p_i|{|}_2^2\]---\left(10\right)$

步骤4验证是否满足终止条件。

本文设置终止目标为或迭代次数达到100次。若不满足条件，再回到步骤2。

图1是一幅充足剂量CT图像的数字模体图像。模体图像共512×512像素，显示窗为[-1000HU，1000HU]，该图主要为了验证少视角投影情况下该方法的重建效果。

仿真使用扇形投影，735个等角探测器单元，X线源扫描轨迹半径为400mm，探测器阵列与X-射线源呈同心圆配置。临床上一次CT扫描通常需要上千视角，为了验证本方法在少视角扫描时的有效性，本实验分别在0～360度内均匀选取了100个视角，并对每组投影数据分别用FBP、SIRT、TV正则化以及本方法进行了图像重建，重建图像均为512×512像素，显示窗为[-1000HU，1000HU]。

图2a图2b图2c图2d给出了用FBP、SIRT、TV正则化以及本方法分别对100个视角的投影数据进行重建得到的CT图像。从重建效果图可以看出，当扫描视角为100时，虽然由于扫描视角过少而在边缘出现了一些伪影，但左右肺部细节结构仍然比其它方法清晰。TV正则化方法由于是一种平滑算法，因此没有出现明显条纹状伪影，但是出现严重的组织结构细节丢失问题。为了更清晰地看到细节保留效果。图3a图3b图3c图3d是对图2a图2b图2c图2d中各子图分别取左肺的一部分进行放大对比结果。

由图3a图3b图3c图3d可以看出，本发明对于细小的血管和气管的重建最为清晰，细节最为明显，重建效果最好。FBP和SIRT方法伪影严重，TV方法虽然没有明显伪影，但是过于平滑，细节丢失过多。

为了验证低X-射线管电流强度投影时的重建效果，实验中是在10mAs的管电流强度扫描得到的噪声模型。仿真使用扇形投影，共888个等角探测器单元，旋转中心与探测器的距离是408.075mm，探测器阵列与X-射线源呈同心圆配置，X-射线源与旋转中心距离541.00mm。本实验＝均匀选取了1000、800和500个视角的投影数据，加入均值为0，方差由式(3)计算得到的高斯噪声，并对每组投影数据分别用FBP、SIRT、软阈值TV正则化以及本文提出方法进行了图像重建，重建图像均为像素，显示窗为[-1000HU，1000HU]。图4a图4b图4c图4d给出了用FBP、SIRT、TV正则化以及本方法分别对低X-射线管电流强度投影数据进行重建得到的CT图像。

由于噪声方差与投影强度呈指数关系，图4a图4b图4c的FBP、SIRT以及TV正则化方法得到的重建结果均受噪声影响较大，重建图像模糊。而本方法由于考虑了投影数据的统计特性，在数据保真项中加入了统计加权，并且正则化系数和在重建过程中是自适应变化的，降低了噪声对重建结果的影响，并且剪切波变换本身具有很强的去噪性能，因此本方法重建出的图像质量远远高于其他重建方法。

为了更清晰地看到细节保留效果，本文对图4中各子图分别取左肺的一部分进行放大对比，如图5a图5b图5c图5d所示。

从图5a图5b图5c图5d可以看出，本方法对于细小气管的重建最为清晰，细节最为明显，重建效果最好。FBP和SIRT方法伪影严重、图像模糊，几乎无法显示任何细节。TV方法虽然伪影较轻，但是过于平滑，细节丢失严重。

为了对本发明进行定量的评价，对图3a图3b图3c图3d中的重建图像与理想模体图像从均方根误差(Root Mean Square Error，RMSE)和图像质量评价指标SSIM(Structural Similarity，SSIM)两个指标来判断重建图像质量，分别见表1和表2。SSIM反映了两幅图像之间的结构相似性，SSIM值越靠近1，图像相似度越高。RMSE定义为

$RMSE=\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}{(f_i-{\tilde{f}}_i)}^2}---\left(11\right)$

其中，f_i表示理想模体图像的第i个像素的HU值，表示重建图像的第i个像素的HU值，N为总的像素数目。

表1 FBP、SIRT、TV正则化和本发明在少视角投影人体肺部仿真实验中重建结果的RMSE对比(单位：HU)

视角数目	FBP	SIRT	TV正则化	本发明方法
					580	132.0781	113.3887	43.3520	29.2340
290	182.0981	154.4146	45.8935	32.3491
					180	272.1809	209.2168	55.7561	36.8280
100	382.5267	327.0988	108.0976	49.3347

RMSE越小，说明重建图像与理想模体的差距越小，由表1可知，本发明方法的RMSE要远远低于另外三种方法。虽然随着投影视角的减少，RMSE有所增加，但性能仍然好于其他三种方法。

表2 FBP、SIRT、TV正则化和本发明在少视角投影人体肺部仿真实验中的重建结果的SSIM对比

视角数目

FBP

SIRT

TV正则化

本发明方法

580	0.6436	0.7045	0.9498	0.9909
					290	0.6416	0.6740	0.9229	0.9901
180	0.5631	0.6348	0.9114	0.9823
					100	0.3717	0.5989	0.8248	0.9622

由表2可以看到，本发明得到的重建图像与理想模体图像的结构相似度最高。

为了对本发明在低X-射线管电流投影进行定量的评价，本实验所得的重建图像和理想模体图像的均方根误差和结构相似度分别见表3和表4。

表3 FBP、SIRT、TV正则化和本发明在10mAs电流强度投影人体肺部仿真实验中的重建结果的RMSE对比(单位：HU)

视角数目	FBP	SIRT	TV正则化	本发明方法
					1000	450.1357	276.3132	108.2189	39.7394
800	452.1827	278.1633	113.7742	41.2592
					500	456.4367	286.0472	149.7718	84.5909

表4 FBP、SIRT、TV正则化和本发明在低X-射线管电流强度投影人体肺部仿真实验中的重建结果的SSIM对比

视角数目	FBP	SIRT	TV正则化	本发明方法
					1000	0.1366	0.3695	0.8279	0.9797
800	0.1271	0.3639	0.8152	0.9634
					500	0.1158	0.3531	0.7507	0.8492

由表3和表4可知，本发明在X-射线管电流强度降低、投影数据信噪比下降时重建的图像的RMSE要远远低于另外三种方法，而结构相似度要远远高于其它三种方法。

Claims (1)

1.一种基于离散剪切波正则化的低剂量CT图像统计重建方法，其特征在于，按照下述步骤进行：

(1)将待建图像在剪切波域可以稀疏表示作为先验信息，在求解Af＝p时加入待建图像的剪切波稀疏表示作为正则化项，该方法可以表示为求解如下问题

$\begin{array}{ccccc}\underset{\mathrm{\α}}{min}\left|\right|\mathrm{\α}|{|}_0& s.t.& f=\mathrm{\ψ}\mathrm{\α}& and& Af=p\end{array}---\left(1\right)$

其中，α为稀疏系数向量，其中只有少数非零元素，||||₀表示l₀范数，ψ表示离散剪切波变换；

(2)利用增广拉格朗日方法将此先验信息作为正则化项加入目标函数，式(1)所对应的增广拉格朗日方程为

$\mathrm{\Φ}(f,\mathrm{\α})=\underset{\mathrm{\α}}{min}\left|\right|\mathrm{\α}|{|}_0-\mathrm{\λ}(f-\mathrm{\ψ}\mathrm{\α})+\frac{\mathrm{\μ}}{2}||f-\mathrm{\ψ}\mathrm{\α}|{|}_2^2-\mathrm{\β}(Af-p)+\frac{\mathrm{\η}}{2}\left|\right|Af-p|{|}_2^2---\left(2\right)$

其中，参数μ、η、λ和β是增广拉格朗日系数，μ和η根据经验取固定值，而向量λ和β分别随稀疏表示误差以及重建图像的投影误差的减小而增大；

(3)噪声强度随着投影数值的增长呈指数增长，为了降低噪声对重建效果的影响，在式(2)的目标函数中加入统计加权其中那么式(2)的目标函数变为

$\mathrm{\Φ}(f,\mathrm{\α})=\underset{\mathrm{\α},f}{\min }\left[\begin{array}{c}\left|\right|\mathrm{\α}|{|}_0-{\mathrm{\λ}}^T(f-\mathrm{\ψ}\mathrm{\α})+\frac{\mathrm{\μ}}{2}||f-\mathrm{\ψ}\mathrm{\α}|{|}_2^2\\ -\underset{i=1}{\overset{M}{\Σ}}{\mathrm{\β}}_i{w}_i({\[Af\]}_i-p_i)+\frac{\mathrm{\η}}{2}\underset{i=1}{\overset{M}{\Σ}}{w}_i\left|\right|{\[Af\]}_i-p_i|{|}_2^2\end{array}\right]---\left(3\right)$ λ和β分别如式(4)和式(5)所示

λ^t＝λ^t-1-μ(f^t-1-ψα^t-1)                      (4)

${\mathrm{\β}}_i^t={\mathrm{\β}}_i^{t-1}-{\mathrm{\ηw}}_i({\[{\mathrm{Af}}^{t-1}\]}_i-p_i)---\left(5\right)$

其中，t≥1表示迭代次数；

(4)式(3)存在两个变量f和α，该方法的实现过程是使用交替最小化的方法来求解这两个变量；首先对f^t-1(t≥1)进行稀疏表示，目标函数为

${\mathrm{\α}}^t=\underset{{\mathrm{\α}}^t}{min}\[\left|\right|{\mathrm{\α}}^t|{|}_0-\mathrm{\λ}(f^{t-1}-{\mathrm{\ψ\α}}^t)+\frac{\mathrm{\μ}}{2}||f^{t-1}-{\mathrm{\ψ\α}}^t|{|}_2^2\]---\left(6\right)$

然后固定式(6)求得的剪切波系数α^t，更新重建图像f^t，这时的优化目标为

$f^t=\underset{f^t}{min}\[\left|\right|{\mathrm{\α}}^t|{|}_1-\mathrm{\λ}(f^t-{\mathrm{\ψ\α}}^t)+\frac{\mathrm{\μ}}{2}||f^t-{\mathrm{\ψ\α}}^t|{|}_2^2-\mathrm{\β}\underset{i=1}{\overset{M}{\Σ}}{w}_i({\[{\mathrm{Af}}^t\]}_i-p_i)+\frac{\mathrm{\η}}{2}\underset{i=1}{\overset{M}{\Σ}}{w}_i\left|\right|{\[{\mathrm{Af}}^t\]}_i-p_i|{|}_2^2\]---\left(7\right)$

(5)本方法设置的终止目标为或迭代次数达到100次，当达到上述目标之一时得到最终重建图像。