CN116385296B - 基于自适应软阈值收缩的太赫兹成像方法 - Google Patents

Abstract

本发明提出了一种基于自适应软阈值收缩的太赫兹成像方法，主要解决现有技术数据采集时间长、成像速度慢、成像质量低的问题。其实现步骤为：1)利用现有的太赫兹时域光谱THz‑TDS系统进行数据采集；2)对采集到的数据进行预处理；3)利用预处理后的数据计算得到衰减系数线积分数据；4)利用衰减系数线积分数据，迭代计算求出太赫兹成像结果。本发明在不改动太赫兹时域光谱仪系统硬件的基础上，通过在图像域进行稀疏采样，降低了数据采集时间，并利用太赫兹图像在小波域具有稀疏性的特点，在降低数据采集量的同时，提高了信号的信噪比，保证了成像质量，可用于非金属物质材料和样品的无损检测。

Description

基于自适应软阈值收缩的太赫兹成像方法

技术领域

本发明属于图像技术领域，更进一步涉及一种太赫兹成像方法，可用于非金属物质材料和样品的无损检测。

背景技术

太赫兹波通常指频率在0.1—10THz范围内的电磁波。该频段的电磁波具有穿透性强、辐射能量低以及光谱分辨率高的特点，尤其是许多材料、分子的振动和转动能级落在太赫兹波段范围内，因此太赫兹时域光谱技术作为一种有效的提取材料光学参数的新兴光谱分析手段，具有无损伤和高时间分辨率的特点，非常适合于测量材料的吸收特性，用于进行定性的鉴别分析工作。

太赫兹时域光谱技术作为最成熟和最成功的太赫兹技术，在过去20年里获得了巨大的成功，是最典型的太赫兹技术之一。太赫兹时域光谱技术通过采集每一个像素点的时域波形可以获得丰富的样品信息，选择时域波形的相位或者振幅信息都可以对样品成像，能同时获得样品的结构分布，密度分布，折射率分布等信息。太赫兹时域光谱逐点成像技术由于扫描比较精细，能获得很高的测量精度，同时受噪声的干扰比较小，信噪比能达到10⁴。但是它的缺陷也很明显，由于成像时要进行逐点扫描，导致成像时间极长。因此，通常适合对小尺寸样品进行成像，如何提高成像速度是一个亟待解决的问题。

2002年由B.Ferguson和张希成等人在期刊Optics Letters(27,期15(2002年8月1日):1312-14)上发表的论文《T-ray computed tomography》中提出了一种4f成像系统，该系统的THz光束由传统成像系统的聚焦光斑扩展成平行波束，更大范围的对目标进行照射，经透镜、偏振片、检偏器这些光学元器件处理后，在接收端采用CCD相机直接录取数据。该CCD相机能够一次性录取全部透射的数据，因此可大大提高成像效率，提升太赫兹成像速率。但是此方法缺点也很明显，由于CCD相机不是相干探测，且无法使用锁定放大器来进行降噪处理，导致信号信噪比较低，大大增加了后期图像处理的工作量，降低了图像质量。且CCD相机价格较贵，增加了硬件开销。

为解决上述问题，Mittleman团队在期刊Applied Physics Letters(93,期12(2008年9月22日):121105)上发表的论文《A single-pixel terahertz imaging systembased on compressed sensing》中提出了一种单像素太赫兹成像方法，其用太赫兹连续波系统照射目标，利用压缩感知对目标进行了稀疏恢复。此方法中太赫兹连续波系统同样将光斑扩展为光束以扩大波束的照射面积。不同的是，接收端不采用CCD录取数据，而是采用单像素探测器进行探测透射信号，仅仅记录信号的强度，因此系统较为简单。实验中采用不同的模板对成像区域进行遮挡，透过的像素区域为1，遮挡区域为0，录取的信号则是透过信号的能量和。此方法的关键在于模板的设计，模板的随机性决定了成像的质量，常用的模板包括高斯随机矩阵、哈达玛矩阵、转盘模板。最后通过压缩感知技术恢复图像。该方法虽然系统简单，但由于需要制作大量的模板，会耗费较长的时间，而且无法获取样品的频谱信息，只能获得强度信息，降低了信号的信噪比，无法恢复出较为清晰的太赫兹图像。

发明内容

本发明的目的在于针对上述现有技术的不足，提出一种基于自适应软阈值收缩的太赫兹成像方法，以在不改动太赫兹时域光谱仪系统硬件的基础上，降低数据录取时间，提高信号的信噪比，恢复出清晰的太赫兹图像。

为实现上述目的，本发明的技术方案是利用太赫兹图像在小波域稀疏的特点，在图像域进行稀疏采样，再利用稀疏恢复算法恢复太赫兹图像，其实现步骤包括如下：

(1)利用现有的太赫兹时域光谱THz-TDS系统进行数据采集：

将THz-TDS系统固定不动，目标放置于扫描架中，成像平面垂直于波束传播方向，假设目标成像区域大小为Q个像素，样本率为μ，则采样点数为μQ，生成高斯随机采样矩阵H，其中μQ个元素为1，其余数据为0；

按照H对目标区域录取数据，即当元素值为1时录取数据，当为0时，则不录取数据；

设每次测量的太赫兹信号点数为K，则录取数据大小为P＝μQ×K；

(2)对采样数据P进行快速傅里叶变换FFT，将数据从时域变换到频域，再采用900-1000GHz的太赫兹频段信号的积分和作为最终的信号，得到大小为μAB×1的预处理后数据S；

(3)根据比尔定律，利用预处理后的数据S计算得到衰减系数线积分数据Y：

Y＝-log(S/max(S))；

(4)利用衰减系数线积分数据Y，经过迭代计算求出太赫兹成像结果X_final：

(4a)初始化参数：包括总的迭代次数J；当前迭代次数j，其初始值j＝1；停止迭代的误差ε；图像X⁽ⁱ⁾，其初始值X⁽⁰⁾＝0，最终成像结果为X_final；Nesterov加速参数t⁽ⁱ⁾；中间变量Z⁽ⁱ⁾；初始残差r＝HX⁽⁰⁾-Y；初始权重矩阵W为单位阵，其大小与图像真值X的大小相同；

(4b)将7个不同的小波基矩阵Ψ₁,Ψ₂,...,Ψ₇和一个狄拉克函数矩阵Ψ_D拼接成一个变换字典矩阵：

(4c)令Nesterov加速参数初始值t⁽⁰⁾＝1；中间变量初始值Z⁽⁰⁾＝X⁽⁰⁾；

(4d)利用图像X⁽ⁱ⁾，计算稀疏系数其中/>为Ψ的逆变换；

(4e)利用中间变量Z⁽ⁱ⁾、变换字典矩阵Ψ、稀疏系数α⁽ⁱ⁾、转置后的随机采样矩阵H^T、残差r，对图像X⁽ⁱ⁾进行更新，计算得到更新后的图像X⁽ⁱ⁺¹⁾：

其中τ表示迭代步长，表示非负限制；

(4f)更新残差r＝HX⁽ⁱ⁺¹⁾-Y，计算该残差r的中位数绝对偏差量MAD(r)，得到自适应软阈值收缩参数λ＝1.482*MAD(r)*W；

(4g)利用更新前图像X⁽ⁱ⁾、更新后图像X⁽ⁱ⁺¹⁾，计算得到对偶变量u⁽ⁱ⁺¹⁾＝2X⁽ⁱ⁺¹⁾-X⁽ⁱ⁾；

(4h)利用对偶变量u⁽ⁱ⁺¹⁾、变换字典矩阵Ψ的逆变换采用软阈值收缩的方法对稀疏系数α⁽ⁱ⁾进行更新，得到更新后的稀疏系数α⁽ⁱ⁺¹⁾：

其中，η表示迭代步长，表示软阈值函数，λ表示软阈值收缩参数；

(4i)更新Nesterov加速参数t⁽ⁱ⁾，得到更新后的Nesterov加速参数t⁽ⁱ⁺¹⁾；

(4j)利用更新前图像X⁽ⁱ⁾、更新后图像X⁽ⁱ⁺¹⁾、更新后Nesterov加速参数t⁽ⁱ⁺¹⁾、计算得到中间变量Z⁽ⁱ⁺¹⁾；

(4k)利用更新前图像X⁽ⁱ⁾、更新后图像X⁽ⁱ⁺¹⁾，计算得到相对误差δ；

(4l)判断误差是否小于设定目标误差ε，即是否满足δ≤ε：

若是，则输出最终图像X_final＝X⁽ⁱ⁺¹⁾；

否则，进行步骤(4m)；

(4m)判断当前迭代次数是否达到最大值N

若i＜N，则令i＝i+1，返回步骤(4c)继续循环；

若i≥N，进行步骤(4n)；

(4n)判断权重矩阵迭代次数是否满足j＜J：

若j＜J，则利用图像X⁽ⁱ⁺¹⁾、软阈值收缩参数λ，计算得到权重矩阵W，令图像初始值X⁽⁰⁾＝X⁽ⁱ⁺¹⁾，令迭代次数j＝j+1，返回步骤(4c)；

否则，得到最终的成像结果X_final＝X⁽ⁱ⁺¹⁾。

本发明与现有技术相比，具有如下优点：

1)本发明由于是在原有THz-TDS系统硬件上实现，因而成本较低，且有利于应用。

2)本发明由于在数据采集过程中进行图像域的随机稀疏采样，可有效降低采样数据量，减小了数据处理过程的计算量，降低了太赫兹成像的时间，提高了太赫兹成像的速度。

3)本发明由于采用能量求和的方法对采样数据P进行预处理，极大的提高了预处理后数据S的信噪比，同时由于采样数据的信噪比对最终的成像质量有重要影响，因此该操作可以使图像质量得到提升。

4)本发明由于通过字典矩阵Ψ将太赫兹图像X⁽ⁱ⁾变换到小波域，利用太赫兹图像在小波域稀疏的特点，在不影响成像质量的情况下，进一步降低了成像时间，提高了成像速度。

5)本发明由于对图像X⁽ⁱ⁾和稀疏系数α⁽ⁱ⁾均进行了更新，且对中间变量Z⁽ⁱ⁺¹⁾进行计算，可避免数据采集过程中随机稀疏采样对成像质量的影响，在采样数据较少的情况下，能得到分辨率较高的太赫兹图像。

6)本发明由于利用残差r和权重矩阵W计算软阈值收缩参数λ，可使软阈值收缩参数λ进行自适应的调整，不仅提高了计算速度，而且提高成像质量。

7)本发明由于利用Nesterov加速参数t⁽ⁱ⁾对迭代成像过程进行加速，减少了成像时间，提高了成像效率。

8)本发明由于利用图像X⁽ⁱ⁺¹⁾和软阈值收缩参数λ计算权重矩阵W，可使权重矩阵W进行自适应的调整，进一步了提高成像质量。

附图说明

图1为本发明的实现流程图；

图2为本发明中使用的太赫兹时域光谱THz-TDS系统图；

图3为实测样品正反两面的照片；

图4为用现有方法对实测样品进行太赫兹全扫描成像的结果图；

图5为分别用本发明和现有方法对实测样品进行太赫兹成像的结果对比图；

图6为用本发明在不同样本率下对实测样品进行太赫兹成像的相对误差图；

图7为用本发明在不同信噪比下对实测样品进行太赫兹成像的相对误差图。

具体实施方式

以下结合附图对本发明的实施例作进一步的详细描述，应当说明的是，具体实施仅仅用以解释本发明，并不用于限定本发明。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供一种基于太赫兹时域光谱技术的成像方法，所述系统包括数据采集设备和数据处理设备。

参照图1，本实例的实现步骤如下：

步骤1，通过太赫兹时域光谱THz-TDS系统完成数据采集。

参照图2，本实例的数据采集设备为太赫兹时域光谱系统THz-TDS，其由发射器、接收器和转台组成。工作时发射器首先将飞秒激光通过分光镜产生两路激光，一路照射到砷化镓晶体上用于产生带宽为0-4THz、动态范围为70dB的太赫兹波束，作为探测波束；另一路作为泵浦波束用于检测。探测波束经过透镜聚焦到转台上的样品上，反射或透射的太赫兹波束再通过透镜聚焦到接收器的探测晶体上，产生电场，接收器对电场信号进行采样，再将采样后的数据发送到信号处理器中进行进一步的处理。

本步骤利用所述THz-TDS系统进行数据采集的实现如下：

.1)将THz-TDS系统固定不动，目标放置于扫描架中，成像平面垂直于波束传播方向，假设目标成像区域大小为Q个像素。设样本率为μ，利用软件根据样本率μ生成大小为Q的高斯随机采样矩阵H，其中H中μQ个元素数值为1，其余数据为0；

1.2)通过扫描架控制目标在成像平面内根据矩阵H进行移动，当H中的元素值为1时，录取数据，当为0时，则不录取数据，移动到下一个位置，直到完成整个目标区域的数据录取，设每次测量的太赫兹信号点数为K，则最终得到大小为μQ×K的采样数据P。

步骤2，对录取数据进行预处理。

太赫兹时域光谱THz-TDS系统录取的数据包含多种信息，因此需要在成像之前进行预处理，现有的预处理方法包括：最大值法、最小值法、峰峰值法、时间延迟法、时域信号能量法、频率选取法、频域能量求和法。在本实施例中，采用频域能量求和法，从而通过提高信号信噪比来增强成像质量，具体实现如下：

首先，对数据进行快速傅里叶变换FFT，将数据从时域变换到频域；

再在频域中选取900-1000GHz的太赫兹频段信号进行积分，最终得到大小为P×1的预处理后数据S。

步骤3，计算衰减系数线积分数据。

在太赫兹成像中，衰减系数线积分是一种用于测量材料吸收和散射的方法。它提供了样品的化学成分、密度、厚度、缺陷和结构这些信息，可用于对材料、医学成像和安全检查领域的成像和检测。本实例成像过程需衰减系数线积分数据，其根据比尔定律的关系式：S＝S₀exp(-Y)，得到衰减系数线积分数据：

Y＝-log(S/max(S))

其中,max(S)为S的最大值，Y为衰减系数线积分，exp(·)表示以自然常数为底的指数函数。

步骤4，利用衰减系数线积分数据Y，经过迭代计算求出太赫兹成像结果X_final。

4a)初始化参数：包括总的迭代次数N；当前迭代次数i，其初始值i＝1；总的权重迭代次数J；当前权重迭代次数j，其初始值j＝1；停止迭代的误差ε；图像X⁽ⁱ⁾，其初始值X⁽⁰⁾＝0；设最终成像结果为X_final；Nesterov加速参数初始值t⁽⁰⁾＝1；中间变量初始值Z⁽⁰⁾＝X⁽⁰⁾；初始权重矩阵W为单位阵，其大小与图像真值X的大小相同；计算初始残差r＝HX⁽⁰⁾-Y；

(4b)将7个不同的小波基矩阵Ψ₁,Ψ₂,...,Ψ₇与一个狄拉克函数矩阵Ψ_D进行拼接，得到变换字典矩阵Ψ及其逆变换矩阵

其中，为Ψ_i的逆变换，/>为Ψ_D的逆变换；

(4c)令Nesterov加速参数初始值t⁽⁰⁾＝1；中间变量初始值Z⁽⁰⁾＝X⁽⁰⁾；

(4d)利用图像X⁽ⁱ⁾，计算稀疏系数其中/>为Ψ的逆变换；

(4e)利用中间变量Z⁽ⁱ⁾、变换字典矩阵Ψ、稀疏系数α⁽ⁱ⁾、转置后的随机采样矩阵H^T、残差r，对图像X⁽ⁱ⁾进行更新，计算得到更新后的图像X⁽ⁱ⁺¹⁾：

其中τ表示迭代步长，表示非负限制；

(4f)更新残差r＝HX⁽ⁱ⁺¹⁾-Y，计算该残差r的中位数绝对偏差量MAD(r)：

MAD(r)＝median(|r-median(r)|)

其中，median(r)表示残差r的中值，根据残差r的中位数绝对偏差量MAD(r)，得到自适应软阈值收缩参数λ＝1.482*MAD(r)*W；

(4g)利用更新前图像X⁽ⁱ⁾、更新后图像X⁽ⁱ⁺¹⁾，计算得到对偶变量u⁽ⁱ⁺¹⁾＝2X⁽ⁱ⁺¹⁾-X⁽ⁱ⁾；

(4h)利用对偶变量u⁽ⁱ⁺¹⁾、逆变换字典采用软阈值收缩的方法对稀疏系数α⁽ⁱ⁾进行更新，得到更新后的稀疏系数：/>

其中为软阈值函数，其表示如下：

式中，α⁽ⁱ⁾表示稀疏系数，η表示迭代步长，u⁽ⁱ⁺¹⁾表示对偶变量，表示变换字典矩阵Ψ的逆变换，λ表示软阈值收缩参数；

(4i)更新Nesterov加速参数t⁽ⁱ⁾，得到更新后的Nesterov加速参数t⁽ⁱ⁺¹⁾：

(4j)利用更新前图像X⁽ⁱ⁾、更新后图像X⁽ⁱ⁺¹⁾、更新后Nesterov加速参数t⁽ⁱ⁺¹⁾、计算得到中间变量Z⁽ⁱ⁺¹⁾：

(4k)利用更新前图像X⁽ⁱ⁾、更新后图像X⁽ⁱ⁺¹⁾，计算得到相对误差δ：

其中，X⁽ⁱ⁾表示更新前图像、X⁽ⁱ⁺¹⁾表示更新后图像，‖ ‖₂表示二范数；

(4l)判断误差是否小于设定目标误差ε，即是否满足δ≤ε：

若是，则输出最终图像X_final＝X⁽ⁱ⁺¹⁾；

否则，进行步骤(4m)；

(4m)判断当前迭代次数是否达到最大值N：

若i＜N，则令i＝i+1，返回步骤(4c)继续循环；

若i≥N，则执行步骤(4n)；

(4n)判断权重矩阵迭代次数是否满足j＜J：

若j＜J，则根据更新后的图像X⁽ⁱ⁺¹⁾，软阈值收缩参数λ，计算得到权重矩阵W：令图像初始值X⁽⁰⁾＝X⁽ⁱ⁺¹⁾，令迭代次数j＝j+1，返回步骤(4c)继续循环；

否则，得到最终的成像结果X_final＝X⁽ⁱ⁺¹⁾。

下面结合仿真实验结果和实测结果对本发明的效果做进一步说明：

1、实验条件：

本发明的实验是在主频2.1GHz、内存64GB的Dell计算机和MATLAB R2021a的环境中进行编程实现的。

实测样品如图3所示，其中图3a为实测样品的正面，图3b为实测样品的反面。

2、实验评价指标：

为了评估本发明的实验效果，对实测样品进行太赫兹全扫描成像，成像结果如图4所示，将该结果作为参考图像，评估准则采用相对误差，计算公式为：

其中X₀为参考图像，X_final为恢复图像。

3、实验内容和结果：

实验一：用本发明和现有方法分别在样本率为0.1、0.3、0.5下，对实测样品进行太赫兹成像，对比其成像效果，结果如图5所示。

从图5可以看出，在不同的样本率下，本发明的成像结果均优于现有方法；即使当样本率降低到0.1时，本发明仍能够很好地恢复图像的整体结构，尽管细节会有所丢失。

实验二：用本发明和现有线性插值法、OMP算法这两种方法分别在样本率为0.1、0.3、0.5下，对实测样品进行太赫兹成像，对比其成像相对误差，结果如表1所示：

表1本发明和现有方法对实测样品进行太赫兹成像的相对误差对比图

样本率	线性插值法	OMP算法	本发明
				0.1	0.1476	0.8617	0.0943
0.3	0.0633	0.5356	0.0264
				0.5	0.0218	0.0768	0.0178

从表1可以看出，在不同的样本率下本发明的相对误差均低于现有方法，,即使当样本率降低到0.1时，本发明的相对误差仍保持在10％以下。

实验三：用本发明在样本率为0.1-0.8下对实测样品进行太赫兹成像，统计不同样本率下本发明的太赫兹成像相对误差，评估本发明的数据压缩能力，结果如图6所示。

从图6可以看出，随着样本率的提高，本发明的太赫兹成像效果越来越好，当样本率达到0.3时，相对误差降至2.5％以下，此时样本率的进一步提高对成像质量的提升不明显。因此，在利用本发明进行实际成像过程中，建议选择0.3的样本率进行太赫兹成像。

实验三：用本发明在信噪比为10dB-60dB下对实测样品进行太赫兹成像，统计不同信噪比下本发明的太赫兹成像相对误差，评估本发明的抗噪能力，结果如图7所示。

从图7可以看出，随着信噪比的提高，本发明的太赫兹成像效果越来越好。在信噪比比较低的情况下，成像效果提升明显，但当信噪比达到40dB后，成像效果提升缓慢。由于在THz-TDS系统中，容易实现40dB的信噪比，因此，在利用本发明进行实际成像过程中，建议选择不低于40dB的信噪比进行太赫兹成像。

上述仿真实验结果表明，本发明提出的基于自适应软阈值收缩的太赫兹稀疏成像算法，充分利用了太赫兹图像在小波域具有稀疏性的特点，通过在图像域进行稀疏采样，并采用稀疏恢复算法，可以在低信噪比和低样本率的条件下稳健地恢复图像，并且在降低数据采集量的同时保证了成像质量。

Claims (6)

1.一种基于自适应软阈值收缩的太赫兹成像方法，其特征在于，包括如下步骤：

(1)利用现有的太赫兹时域光谱THz-TDS系统进行数据采集：

将THz-TDS系统固定不动，目标放置于扫描架中，成像平面垂直于波束传播方向，假设目标成像区域大小为Q个像素，样本率为μ，则采样点数为μQ，生成高斯随机采样矩阵H，其中μQ个元素为1，其余数据为0；

按照H对目标区域录取数据，即当元素值为1时录取数据，当为0时，则不录取数据；

设每次测量的太赫兹信号点数为K，则录取数据大小为P＝μQ×K；

(2)对采样数据P进行快速傅里叶变换FFT，将数据从时域变换到频域，再采用900-1000GHz的太赫兹频段信号的积分和作为最终的信号，得到大小为μAB×1的预处理后数据S；

(3)根据比尔定律，利用预处理后的数据S计算得到衰减系数线积分数据Y：

Y＝-log(S/max(S))；

(4)利用衰减系数线积分数据Y，经过迭代计算求出太赫兹成像结果X_final：

(4a)初始化参数：包括总的迭代次数J；当前迭代次数j，其初始值j＝1；停止迭代的误差ε；图像X⁽ⁱ⁾，其初始值X⁽⁰⁾＝0，最终成像结果为X_final；Nesterov加速参数t⁽ⁱ⁾；中间变量Z⁽ⁱ⁾；初始残差r＝HX⁽⁰⁾-Y；初始权重矩阵W为单位阵，其大小与图像真值X的大小相同；

(4b)将7个不同的小波基矩阵Ψ₁,Ψ₂,...,Ψ₇和一个狄拉克函数矩阵Ψ_D拼接成一个变换字典矩阵：

(4c)令Nesterov加速参数初始值t⁽⁰⁾＝1；中间变量初始值Z⁽⁰⁾＝X⁽⁰⁾；

(4d)利用图像X⁽ⁱ⁾，计算稀疏系数其中/>为Ψ的逆变换；

(4e)利用中间变量Z⁽ⁱ⁾、变换字典矩阵Ψ、稀疏系数α⁽ⁱ⁾、转置后的随机采样矩阵H^T、残差r，对图像X⁽ⁱ⁾进行更新，计算得到更新后的图像X⁽ⁱ⁺¹⁾：

其中τ表示迭代步长，表示非负限制；

(4f)更新残差r＝HX⁽ⁱ⁺¹⁾-Y，计算该残差r的中位数绝对偏差量MAD(r)，得到自适应软阈值收缩参数λ＝1.482*MAD(r)*W；

(4g)利用更新前图像X⁽ⁱ⁾、更新后图像X⁽ⁱ⁺¹⁾，计算得到对偶变量u⁽ⁱ⁺¹⁾＝2X⁽ⁱ⁺¹⁾-X⁽ⁱ⁾；

(4h)利用对偶变量u⁽ⁱ⁺¹⁾、变换字典矩阵Ψ的逆变换采用软阈值收缩的方法对稀疏系数α⁽ⁱ⁾进行更新，得到更新后的稀疏系数α⁽ⁱ⁺¹⁾：

其中，η表示迭代步长，表示软阈值函数，λ表示软阈值收缩参数；

(4i)更新Nesterov加速参数t⁽ⁱ⁾，得到更新后的Nesterov加速参数t⁽ⁱ⁺¹⁾；

(4j)利用更新前图像X⁽ⁱ⁾、更新后图像X⁽ⁱ⁺¹⁾、更新后Nesterov加速参数t⁽ⁱ⁺¹⁾、计算得到中间变量Z⁽ⁱ⁺¹⁾；

(4k)利用更新前图像X⁽ⁱ⁾、更新后图像X⁽ⁱ⁺¹⁾，计算得到相对误差δ；

(4l)判断误差是否小于设定目标误差ε，即是否满足δ≤ε：

若是，则输出最终图像X_final＝X⁽ⁱ⁺¹⁾；

否则，进行步骤(4m)；

(4m)判断当前迭代次数是否达到最大值N

若i＜N，则令i＝i+1，返回步骤(4c)继续循环；

若i≥N，进行步骤(4n)；

(4n)判断权重矩阵迭代次数是否满足j＜J：

若j＜J，则利用图像X⁽ⁱ⁺¹⁾、软阈值收缩参数λ，计算得到权重矩阵W，令图像初始值X⁽⁰⁾＝X⁽ⁱ⁺¹⁾，令迭代次数j＝j+1，返回步骤(4c)；

否则，得到最终的成像结果X_final＝X⁽ⁱ⁺¹⁾；

所述权重矩阵W，表示如下：

其中，X⁽ⁱ⁺¹⁾表示更新后的图像，λ表示软阈值收缩参数。

2.根据权利要求1所述的方法，其特征在于，步骤(4f)中计算残差r的中位数绝对偏差量MAD(r)，公式如下：

MAD(r)＝median(|r-median(r)|)

其中，median(r)表示残差r的中值。

3.根据权利要求1所述的方法，其特征在于，步骤(4h)中的软阈值函数表示如下：

其中，α⁽ⁱ⁾表示稀疏系数，η表示迭代步长，u⁽ⁱ⁺¹⁾表示对偶变量，表示变换字典矩阵Ψ的逆变换，λ表示软阈值收缩参数。

4.根据权利要求1所述的方法，其特征在于，步骤(4i)中得到更新后的Nesterov加速参数t⁽ⁱ⁺¹⁾，表示如下：

其中，t⁽ⁱ⁾表示更新前的Nesterov加速参数。

5.根据权利要求1所述的方法，其特征在于，步骤(4j)中得到的中间变量Z⁽ⁱ⁺¹⁾,表示如下：

其中，X⁽ⁱ⁾表示更新前图像，X⁽ⁱ⁺¹⁾表示更新后图像，t⁽ⁱ⁺¹⁾表示更新后Nesterov加速参数。

6.根据权利要求1所述的方法，其特征在于，步骤(4k)中计算相对误差δ，公式如下：

其中，X⁽ⁱ⁾表示更新前图像、X⁽ⁱ⁺¹⁾表示更新后图像，‖‖₂表示二范数。