CN105930310A - 一种自适应正则化平滑l0范数方法 - Google Patents

Abstract

本发明公开了一种自适应正则化平滑l₀范数方法，对正则化SL0算法进行了改进，在内循环的最速上升法中以第一次迭代的信号残差项估计值以及该迭代前后的稀疏信号估计的偏差值作为当前正则化参数的选择依据，从而能自适应地调整在每次外循环中的信号稀疏度和误差容许项的权重值，在优化过程中保持两者的平衡性，从而有效降低稀疏信号的重构误差，提高了算法的抗噪声干扰能力；通过引入SVD方法来避免在迭代过程中投影到可行解集的操作中的大规模矩阵求逆运算，有效提高本发明方法对稀疏信号的重构速度。

Description

一种自适应正则化平滑l0范数方法

技术领域

本发明涉及一种自适应正则化平滑l₀范数方法，属于压缩感知恢复技术领域。

背景技术

压缩感知作为信号处理领域的一项新技术，已被广泛地应用于生物医学、图像处理、无线通信和雷达信号处理等领域。压缩感知理论通过求解l₀范数最小化问题，能从少量的非自适应投影测量值中以较高概率重构出稀疏信号。然而，l₀范数最小化问题是NP-hard问题，需要通过组合搜索求解，但信号维度较大时难以实现该问题的求解。平滑l₀范数(Smoothed l₀Norm，SL0)算法是利用一系列高斯函数来近似l₀范数，从而将l₀范数最小化的NP-hard问题转化为易求解的平滑函数最小化问题，然后通过求解该问题能获得稀疏信号的重构值。SL0算法具有计算速度快，而且重构信号所需的测量值较少的优点。SL0算法通过求解以下等式约束最优化问题：

max F_σ(x)s.t.y＝Ax (1)

式中，y为已知的m×1维向量，x为已知的n×1维信号向量，A为已知的m×n维观测矩阵，且σ为高斯函数形状参数，x_i表示向量x的第i个元素，N为向量x中元素个数。在实际环境中，由于噪声的存在导致y和Ax之间存在误差，而SL0算法使用y＝Ax的等式约束条件，因此其稀疏信号重构性能在噪声环境下会严重恶化。为了提高该算法的抗噪能力，文献(Hongxia Bu,Ran Tao,Xia Bai,JuanZhao，Regularized smoothed l0 norm algorithm and its application to CS-basedradar imaging[J],Signal processing,2016,122:115-122.)提出了一种正则化SL0算法，该算法在SL0算法的目标函数中加入一个误差容许项从而将式(1)转换成一个适用于噪声环境下信号重构的目标优化函数，即

$\begin{array}{ccc}\max F_{\mathrm{\σ}}\left(x\right)& s.t.& \left|\right|y=Ax|{|}_2^2\≤\mathrm{\δ}\end{array}---\left(2\right)$

式中，δ为容许的误差值。正则化SL0算法利用拉格朗日乘数法构建无约束目标优化函数，类似于SL0算法，它利用2个嵌套迭代运算来获得稀疏信号x的估计值。该算法在所有迭代运算过程中正则化参数均为一个固定值，然而在实际迭代运算过程中，误差容许项的值会发生较大的改变，因此在不同迭代运算中，固定的正则化参数无法有效保持信号稀疏度和误差容许项之间的平衡，导致该算法的抗噪声能力不强、稳健性较低。针对实际应用场景对压缩感知算法的需求,因此研究一种抗噪声性能强以及正则化参数易于自适应调整的SL0算法是非常有必要的。

发明内容

为了解决上述技术问题，本发明提供了一种自适应正则化平滑l₀范数方法。

为了达到上述目的，本发明所采用的技术方案是：

一种自适应正则化平滑l₀范数方法，包括以下步骤：

步骤1，初始化；

将y＝Ax最小二乘解作为算法运行的初始值，定义j为外循环的迭代计数值，初始时j＝1，j＝1时的高斯函数形状参数σ_j为第j次外循环迭代中高斯函数形状参数，σ_J外循环终止时的高斯函数形状参数，ρ为收缩因子，0＜ρ＜1，L为内循环的最大迭代次数；

步骤2，在可行解集χ＝{x|||y-Ax||₂≤δ}上利用最速上升法求使F_σ(x)最大值的解；

具体如下：

A1)令σ＝σ_j，l＝1，转至步骤A2；

A2)判断l是否小于等于L，如果是，则转至步骤A3，如果不是，则转至步骤A8；

A3)x←x-μxexp[-x²/(2σ²)]，转至步骤A4；

其中，μ为大于0的常数；

A4)判断l是否等于1，如果是，则转至步骤A5，如果不是，则转至步骤A6；

A5)自适应调整正则化参数转至步骤A6；

A6)将x投影到可行解集χ＝{x|||y-Ax||₂≤δ}上，即x←x-A^HU[(b1_1×m)⊙U^H](Ax-y)，转至步骤A7；

其中，U∈C^m×m称为左奇异矩阵，1_1×m为1×m维全1向量，b为m×1维向量，[·]^T表示转置，α₁,α₂,…,α_m为对角矩阵Σ'＝Σ·∑^H的对角元素值，∑∈C^m×n为由特征值组成的对角矩阵，U和∑是通过计算A的SVD变换获得，

A7)l＝l+1，转至步骤A2；

A8)令

A9)判断σ_j是否小于σ_J，如果是，则为使F_σ(x)最大值的解，即为信号最稀疏表示解，如果不是，则j＝j+1，σ_j＝ρσ_j-1，转至步骤A1。

自适应调整正则化参数的推导过程为，

在给定σ＝σ_j的迭代循环过程中的优化问题可采用拉格朗日形式表示为，

${\hat{x}}^{\left(j\right)}=\arg \underset{x^{\left(j\right)}}{min}\left(\right||x^{\left(j\right)}-x|{|}_2^2+\mathrm{\λ}\left|\right|y-{\mathrm{Ax}}^{\left(j\right)}\left|{|}_2^2\right)---\left(3\right)$

利用加权最小二乘法求解(3)得，

${\hat{x}}^{\left(j\right)}=x-A^H{({\mathrm{AA}}^H+{\mathrm{\λ}}^{-1}{I}_m)}^{-1}(A\hat{x}-y)---\left(4\right)$

其中，I_m为m×m维单位矩阵，

式(4)，将投影到可行解集χ＝{x|||y-Ax||₂≤δ}上，

令函数

由式(3)可知，当时，函数H(x^(j),λ)达到最小值，

函数H(x(^j),λ)可以改写为，

H(x^(j),λ)＝(x^(j)-x)^H(x^(j)-x)+λ(y-Ax^(j))^H(y-Ax^(j)) (6)

＝x^(j)Hx^(j)-x^(j)Hx-x^Hx^(j)+x^Hx+λ(y^Hy-y^HAx^(j)-x^(j)HA^Hy+x^(j)HA^HAx^(j))

将式(6)关于x^(j)进行复数求导，

$\frac{\∂H(x^{\left(j\right)},\mathrm{\λ})}{\∂x^{\left(j\right)}}=2(x^{\left(j\right)}-x)+2\mathrm{\λ}(A^H{\mathrm{Ax}}^{\left(j\right)}-A^{H}y)---\left(7\right)$

令即可得，

(x^(j)-x)+λ(A^HAx^(j)-A^Hy)＝0 (8)

将带入式(8)，得到，

$\mathrm{\λ}(A^H{\mathrm{Ax}}^{\left(j\right)}-A^{H}y)=x-{\hat{x}}^{\left(j\right)}---\left(9\right)$

将式(9)等号两边取2范数操作，即，

$\mathrm{\λ}\left|\right|A^H{\mathrm{Ax}}^{\left(j\right)}-A^{H}y|{|}_2=||x-{\hat{x}}^{\left(j\right)}|{|}_2---\left(10\right)$

从式(10)中可以得出正则化参数的估计值，

$\widehat{\mathrm{\λ}}=\left|\right|x-{\hat{x}}^{\left(j\right)}|{|}_2/|\left|A^H(A{\hat{x}}^{\left(j\right)}-y)\right|{|}_2---\left(11\right)$

由式(11)和式(4)可知，的解依赖于的值，在利用式(11)计算参数值时是未知的，因此，的值可由在σ＝σ_j-1时的稀疏信号重构值来代替，

即，正则化参数

x←x-A^HU[(b1_1×m)⊙U^H](Ax-y)的推导过程为，

计算A的SVD变换可得，

A＝UΣV^H (13)

其中，V∈C^m×m称为右奇异矩阵，

利用U、V的特性可得，

AA^H＝UΣΣ^H U^H (14)

UU^H＝I_m (15)

将式(12)中新的正则化参数、式(13)和(14)带入式(4)得，

$\begin{array}{c}{\hat{x}}^{\left(j\right)}=x-A^H{({\mathrm{U\Σ\Σ}}^H{U}^H+{\widehat{\mathrm{\λ}}}^{-1}{\mathrm{UU}}^H)}^{-1}(Ax-y)\\ =x-A^H{\[U({\mathrm{\Σ\Σ}}^H+{\widehat{\mathrm{\λ}}}^{-1}{I}_m)U^H\]}^{-1}(Ax-y)\\ =x-A^{H}U{({\mathrm{\Σ\Σ}}^H+{\widehat{\mathrm{\λ}}}^{-1}{I}_m)}^{-1}{U}^H(Ax-y)\end{array}---\left(16\right)$

其中，ΣΣ^H＝diag(α₁,α₂,…,α_m)，diag(·)为对角化操作，

令m×1维向量

则，

(ΣΣ^H+λ^-1I_m)U^H＝(b1_1×m)⊙U^H (17)

式中，⊙表示Hadamard积，

那么式(16)可改写为，

A的SVD变换在信息重构处理前进行离线计算。

本发明所达到的有益效果：1、本发明对正则化SL0算法进行了改进，在内循环的最速上升法中以第一次迭代的信号残差项估计值以及该迭代前后的稀疏信号估计的偏差值作为当前正则化参数的选择依据，从而能自适应地调整在外循环迭代中的信号稀疏度和误差容许项的权重值，在优化过程中保持两者的平衡性，从而有效降低稀疏信号的重构误差，提高了算法的抗噪声干扰能力；2、本发明提高了对稀疏信号的重构速度，通过引入SVD方法来避免在迭代过程中投影到可行解集的操作中的大规模矩阵求逆运算，并对A的SVD分解以及与奇异矩阵相乘操作进行离线计算并存储，有效提高本发明方法对稀疏信号的重构速度。

附图说明

图1为本发明的流程图。

图2为不同算法的的重构信噪比SER与噪声功率σ′_n的变化关系。

图3为不同算法的信号重构均方误差MSE与噪声功率σ′_n的变化关系。

图4为不同算法的信号重构速度与噪声功率σ′_n的变化关系。

图5为不同算法的重构信噪比SER与信号稀疏度K的变化关系。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

本发明类似于SL0算法和正则化SL0算法，也采用2个嵌套循环运算从式(2)中来获得信号x的最稀疏表示解。当σ较小时，函数F_σ(x)出现高度非光滑现象，从而导致许多局部极小值出现，不易进行优化；而当σ较大时，虽然函数F_σ(x)较为光滑，则有利于进行优化，但是稀疏信号x的重构误差较大。因此采取逐步减小σ的策略来避免优化F_σ(x)过程中陷入局部最大值，针对每个σ值，在可行解集χ＝{x|||y-Ax||₂≤δ}上寻找使得F_σ(x)最大值的x值，并将该x值作为下一次迭代的初值。与正则化SL0算法类似，本发明方法的内循环过程包含最速上升方法步骤x←x+μσ²▽F_σ(x)(其中μ为大于0的常数、▽F_σ(x)为F_σ(x)的梯度值)以及优化问题的求解。

在给定σ＝σ_j的迭代循环过程中的优化问题可采用拉格朗日形式表示为，

${\hat{x}}^{\left(j\right)}=\arg \underset{x^{\left(j\right)}}{min}\left(\right||x^{\left(j\right)}-x|{|}_2^2+\mathrm{\λ}\left|\right|y-{\mathrm{Ax}}^{\left(j\right)}\left|{|}_2^2\right)---\left(3\right)$

式中，λ为正则化参数，

由于的解即为使F_σ(x)达到最大的稀疏信号重构值，因此，λ其实是调节信号稀疏度和信号残差项在目标函数值所占的比重，

利用加权最小二乘法求解(3)得，

${\hat{x}}^{\left(j\right)}=x-A^H{({\mathrm{AA}}^H+{\mathrm{\λ}}^{-1}{I}_m)}^{-1}(Ax-y)---\left(4\right)$

其中，I_m为m×m维单位矩阵，

式(4)，将投影到可行解集χ＝{x|||y-Ax||₂≤δ}上，

由于在实际迭代运算过程中，信号残差项的值会发生较大改变，为了保持信号稀疏度和信号残差项在目标函数值中的平衡性，以增强算法的抗噪声性能，可提出了一种正则化参数的自适应调整方法，即

令函数

由式(3)可知，当时，函数H(x^(j),λ)达到最小值，

函数H(x^(j),λ)可以改写为，

H(x^(j),λ)＝(x^(j)-x)^H(x^(j)-x)+λ(y-Ax^(j))^H(y-Ax^(j)) (6)

＝x^(j)Hx^(j)-x^(j)Hx-x^Hx^(j)+x^Hx+λ(y^Hy-y^HAx^(j)-x^(j)HA^Hy+x^(j)HA^HAx^(j))

为了能使H(x^(j),λ)达到最小值，将式(6)关于x^(j)进行复数求导，

$\frac{\∂H(x^{\left(j\right)},\mathrm{\λ})}{\∂x^{\left(j\right)}}=2(x^{\left(j\right)}-x)+2\mathrm{\λ}(A^H{\mathrm{Ax}}^{\left(j\right)}-A^{H}y)---\left(7\right)$

令即可得，

(x^(j)-x)+λ(A^HAx^(j)-A^Hy)＝0 (8)

将带入式(8)，得到，

$\mathrm{\λ}(A^H{\mathrm{Ax}}^{\left(j\right)}-A^{H}y)=x-{\hat{x}}^{\left(j\right)}---\left(9\right)$

将式(9)等号两边取2范数操作，即，

$\mathrm{\λ}\left|\right|A^H{\mathrm{Ax}}^{\left(j\right)}-A^{H}y|{|}_2=||x-{\hat{x}}^{\left(j\right)}|{|}_2---\left(10\right)$

从式(10)中可以得出正则化参数的估计值，

$\widehat{\mathrm{\λ}}=\left|\right|x-{\hat{x}}^{\left(j\right)}|{|}_2/|\left|A^H(A{\hat{x}}^{\left(j\right)}-y)\right|{|}_2---\left(11\right)$

由式(11)和式(4)可知，的解依赖于的值，在利用式(11)计算参数值时是未知的，因此，的值可由在σ＝σ_j-1时的稀疏信号重构值来代替，

即，正则化参数

由式(12)可知，在迭代过程中正则化参数能根据信号残差项估计值以及最速上升法的迭代前后的稀疏信号估计的偏差值进行自适应调整，平衡信号稀疏度和信号残差项在目标函数值中的比重。

在给定σ＝σ_j所对应的内循环过程中，类似于SL0方法，本发明方法也是利用最速上升方法求解式(3)所表示的最小化优化问题。假设最速上升方法的迭代次数为L，即内循环的最大迭代次数为L，由于在L次迭代内x根据公式x←x+μσ²▽F_σ(x)进行更新，则正则化参数和值也应分别更新L次，计算量较大。考虑到内循环中并不要求最速上升法达到收敛即迭代次数L可以很小，那么正则化参数在L次迭代内变化不会太大，即可认为在L次迭代中正则化参数固定不变，因此在内循环的L次迭代中的大规模矩阵求逆只需要运算一次即可，提高了计算效率。

在逐步减小σ的外循环中，正则化参数需要重新计算更新，则在外循环的每次迭代中的大规模矩阵求逆也需要重新计算。为了进一步提高算法的实时性，本发明在引入SVD方法来降低在迭代过程中投影到可行解集时的大规模矩阵求逆运算。对A进行SVD变换可得，

A＝U∑V^H (13)

其中，V∈C^m×m称为右奇异矩阵，

利用U、V的特性可得，

AA^H＝UΣΣ^H U^H (14)

UU^H＝I_m (15)

将式(12)中新的正则化参数、式(14)和(15)带入式(4)得，

$\begin{array}{c}{\hat{x}}^{\left(j\right)}=x-A^H{({\mathrm{U\Σ\Σ}}^H{U}^H+{\widehat{\mathrm{\λ}}}^{-1}{\mathrm{UU}}^H)}^{-1}(Ax-y)\\ =x-A^H{\[U({\mathrm{\Σ\Σ}}^H+{\widehat{\mathrm{\λ}}}^{-1}{I}_m)U^H\]}^{-1}(Ax-y)\\ =x-A^{H}U{({\mathrm{\Σ\Σ}}^H+{\widehat{\mathrm{\λ}}}^{-1}{I}_m)}^{-1}{U}^H(Ax-y)\end{array}---\left(16\right)$

其中，ΣΣ^H＝diag(α₁,α₂,…,α_m)，diag(·)为对角化操作，

令m×1维向量

则，

(∑∑^H+λ^-1I_m)U^H＝(b1_1×m)⊙U^H (17)

式中，⊙表示Hadamard积，

那么式(16)可改写为，

式(18)已不存在大规模矩阵求逆运算，因此其计算效率要，而且由于在稀疏信号重构前观测矩阵A是已知的，因此A的SVD变换可在信息重构处理前进行离线计算，这样能有效提高本发明方法对稀疏信号的重构速度。

综上所述，本发明的具体流程如图1所示，一种自适应正则化平滑l₀范数方法，包括以下步骤：

步骤1，离线计算；

离线计算A的SVD变换，获得U、V和∑，并计算A^HU和U^H；

其中，U∈C^m×m和V∈C^n×n为正交的酉矩阵，分别称为左奇异矩阵和右奇异矩阵，Σ∈C^m×n为对角矩阵；

步骤2，初始化；

将y＝Ax最小二乘解作为算法运行的初始值，定义j为外循环的迭代计数值，初始时j＝1，j＝1时的高斯函数形状参数σ_j为第j次外循环迭代中高斯函数形状参数，σ_J外循环终止时的高斯函数形状参数，ρ为收缩因子，0＜ρ＜1，L为内循环的最大迭代次数；

步骤3，在可行解集χ＝{x|||y-Ax||₂≤δ}上利用最速上升法求使F_σ(x)最大值的解；

具体如下：

A1)令σ＝σ_j，l＝1，转至步骤A2；

A2)判断l是否小于等于L，如果是，则转至步骤A3，如果不是，则转至步骤A8；

A3)x←x-μxexp[-x²/(2σ²)]，转至步骤A4；

其中，μ为大于0的常数；

A4)判断l是否等于1，如果是，则转至步骤A5，如果不是，则转至步骤A6；

A5)自适应调整正则化参数转至步骤A6；

A6)将x投影到可行解集χ＝{x|||y-Ax||₂≤δ}上，即x←x-A^HU[(b1_1×m)⊙U^H](Ax-y)，转至步骤A7；

其中，U∈C^m×m称为左奇异矩阵，1_1×m为1×m维全1向量，b为m×1维向量，[·]^T表示转置，α₁,α₂,…,α_m为对角矩阵∑'＝∑·∑^H的对角元素值，Σ∈C^m×n为由特征值组成的对角矩阵，U和Σ是通过计算A的SVD变换获得，

A7)l＝l+1，转至步骤A2；

A8)令

A9)判断σ_j是否小于σ_J，如果是，则为使F_σ(x)最大值的解，即为信号最稀疏表示解，如果不是，则j＝j+1，σ_j＝ρσ_j-1，转至步骤A1。

为了进一步说明上述方法，做以下仿真实验。

在该仿真实验中，为了检验SL0算法(SL0)、表面欠定方程组求解(FocalUndetermined System Solver，FOCUSS)方法、正则化SL0算法(ReSL0)、自适应正则化SL0算法(AReSL0)以及引入SVD的自适应正则化SL0算法(AReSL0_SVD)性能，分别它们将应用于在噪声环境下的稀疏信号重构问题中。仿真参数设置：L＝3，μ＝2.5，所有算法处理的原始信号都为长度n＝1500的高斯稀疏信号，其中非零元素个数为K＝120。稀疏信号的非零元素值服从高斯分布，且在每次独立实验中非零元素的位置随机选取。在噪声环境下，观测信号可以表示为y＝Ax+n。每次实验的观测信号y长度为m＝500。在每次实验中，A的列元素服从标准的独立同分布的高斯分布。重构信号比SER(dB)定义为信号重构均方误差定义为其中为原始信号x的估计值。以下仿真实验分别独立进行100次，重构信号比SER、重构均方误差MSE以及运行时间都为200次独立实验的平均值。

仿真内容1：噪声对算法重构性能的影响

图2和图3分别为不同算法的重构信噪比SER和信号重构均方误差MSE与噪声功率σ′_n的变化关系图。由图2和图3可知，当噪声功率σ′_n逐渐增加即信噪比降低时，所有算法的信号重构性能都会随之下降；SL0算法和FOCUSS算法的重构信噪比SER和信号重构均方误差MSE基本类似；由于在目标函数里加入了误差容许项，ReSL0算法在噪声环境下的重构性能要优于SL0算法和FOCUSS算法；在ReSL0算法基础上，提出了一种自适应正则化SL0算法(AReSL0算法)，在每次迭代过程中自适应地调整信号稀疏度和误差容许项的权重值，保持两者的平衡性，从而有效降低稀疏信号的重构误差，其稀疏信号重构性能优于ReSL0算法。为了降低AReSL0算法的运算量，通过引入SVD方法来避免了在迭代过程中的大规模矩阵求逆运算，AReSL0_SVD算法的性能与AReSL0算法几乎一致，因此SVD变换的引入并没有使得AReSL0算法的重构性能降低。

仿真内容2：算法的运算时间对比

以CPU运算时间作为不同算法复杂度的判断依据，虽然CPU运算时间不能对算法的复杂度进行准确测量，但是可以粗略评价算法的复杂度。本实验在MATLAB R2013a中完成，计算机配置为：Intel(R)Core(TM)i5-4570处理器、主频为3.2GHz、内存为4GB。图4描述了不同算法的信号重构速度。由图4可知，SL0算法的运算速度最快，ReSL0算法的重构速度其次，AReSL0算法虽然在重构性能上要明显优于SL0算法和ReSL0算法，但是其运算复杂度要高于两者，通过引入SVD方法来提高重构速度，相比AReSL0算法，AReSL0_SVD算法的稀疏信号重构时间约降低了27％左右。但是AReSL0算法和AReSL0_SVD算法的重构时间都比FOCUSS算法要低得多。

仿真内容3：算法的算法重构性能与信号稀疏度K的关系

图5为不同算法的重构信噪比SER与信号稀疏度K的变化关系，其中信号稀疏度K在30～380之间取值，其它仿真参数同仿真实验1。由图5可知，原始信号具有不同稀疏度K的情况下，AReSL0算法和AReSL0_SVD算法重构性能均要好于SL0算法、ReSL0算法和FOCUSS算法，尤其是在信号稀疏度K较小时它们要明显优于上述三种对比算法。

综上所述，本发明在内循环的最速上升法中以第一次迭代的信号残差项估计值以及该迭代前后的稀疏信号估计的偏差值作为当前正则化参数的选择依据，从而能自适应地调整在外循环迭代中的信号稀疏度和误差容许项的权重值，在优化过程中保持两者的平衡性，从而有效降低稀疏信号的重构误差，提高了算法的抗噪声干扰能力；通过引入SVD方法来避免在迭代过程中投影到可行解集的操作中的大规模矩阵求逆运算，并对A的SVD分解以及与奇异矩阵相乘操作进行离线计算并存储，有效提高本发明方法对稀疏信号的重构速度。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (4)

1.一种自适应正则化平滑l₀范数方法，其特征在于：包括以下步骤，

步骤1，初始化；

将y＝Ax最小二乘解作为算法运行的初始值，定义j为外循环的迭代计数值，初始时j＝1，j＝1时的高斯函数形状参数σ_j为第j次外循环迭代中高斯函数形状参数，σ_J外循环终止时的高斯函数形状参数，ρ为收缩因子，0＜ρ＜1，L为内循环的最大迭代次数；

步骤2，在可行解集χ＝{x|||y-Ax||₂≤δ}上利用最速上升法求使F_σ(x)最大值的解；

具体如下：

A1)令σ＝σ_j，l＝1，转至步骤A2；

A2)判断l是否小于等于L，如果是，则转至步骤A3，如果不是，则转至步骤A8；

A3)x←x-μxexp[-x²/(2σ²)]，转至步骤A4；

其中，μ为大于0的常数；

A4)判断l是否等于1，如果是，则转至步骤A5，如果不是，则转至步骤A6；

A5)自适应调整正则化参数转至步骤A6；

A6)将x投影到可行解集χ＝{x|||y-Ax||₂≤δ}上，即转至步骤A7；

其中，U∈C^m×m称为左奇异矩阵，1_1×m为1×m维全1向量，b为m×1维向量，[·]^T表示转置，α₁,α₂,…,α_m为对角矩阵Σ'＝Σ·Σ^H的对角元素值，Σ∈C^m×n为由特征值组成的对角矩阵，U和Σ是通过计算A的SVD变换获得，

A7)l＝l+1，转至步骤A2；

A8)令

A9)判断σ_j是否小于σ_J，如果是，则为使F_σ(x)最大值的解，即为信号最稀疏表示解，如果不是，则j＝j+1，σ_j＝ρσ_j-1，转至步骤A1。

2.根据权利要求1所述的一种自适应正则化平滑l₀范数方法，其特征在于：自适应调整正则化参数的推导过程为，

在给定σ＝σ_j的迭代循环过程中的优化问题可采用拉格朗日形式表示为，

${\hat{x}}^{\left(j\right)}=\arg \underset{x^{\left(j\right)}}{min}\left(\right||x^{\left(j\right)}-x|{|}_2^2+\mathrm{\λ}\left|\right|y-{\mathrm{Ax}}^{\left(j\right)}\left|{|}_2^2\right)---\left(3\right)$

利用加权最小二乘法求解(3)得，

${\hat{x}}^{\left(j\right)}=x-A^H{({\mathrm{AA}}^H+{\mathrm{\λ}}^{-1}{I}_m)}^{-1}(A\hat{x}-y)---\left(4\right)$

其中，I_m为m×m维单位矩阵，

式(4)，将投影到可行解集χ＝{x|||y-Ax||₂≤δ}上，

令函数

由式(3)可知，当时，函数H(x^(j),λ)达到最小值，

函数H(x^(j),λ)可以改写为，

$\begin{array}{c}H\left(x^{\left(j\right)},\mathrm{\λ}\right)={\left(x^{\left(j\right)}-x\right)}^H\left(x^{\left(j\right)}-x\right)+\mathrm{\λ}{\left(y-{\mathrm{Ax}}^{\left(j\right)}\right)}^H\left(y-{\mathrm{Ax}}^{\left(j\right)}\right)\\ =x^{\left(j\right)H}{x}^{\left(j\right)}-x^{\left(j\right)H}x-x^H{x}^{\left(j\right)}+x^{H}x+\mathrm{\λ}\left(y^{H}y-y^H{\mathrm{Ax}}^{\left(j\right)}-x^{\left(j\right)H}{A}^{H}y+x^{\left(j\right)H}{A}^H{\mathrm{Ax}}^{\left(j\right)}\right)\end{array}---\left(6\right)$

将式(6)关于x^(j)进行复数求导，

$\frac{\∂H(x^{\left(j\right)},\mathrm{\λ})}{\∂x^{\left(j\right)}}=2(x^{\left(j\right)}-x)+2\mathrm{\λ}(A^H{\mathrm{Ax}}^{\left(j\right)}-A^{H}y)---\left(7\right)$

令即可得，

(x^(j)-x)+λ(A^HAx^(j)-A^Hy)＝0 (8)

将带入式(8)，得到，

$\mathrm{\λ}(A^H{\mathrm{Ax}}^{\left(j\right)}-A^{H}y)=x-{\hat{x}}^{\left(j\right)}---\left(9\right)$

将式(9)等号两边取2范数操作，即，

$\mathrm{\λ}\left|\right|A^H{\mathrm{Ax}}^{\left(j\right)}-A^{H}y|{|}_2=||x-{\hat{x}}^{\left(j\right)}|{|}_2---\left(10\right)$

从式(10)中可以得出正则化参数的估计值，

$\widehat{\mathrm{\λ}}=\left|\right|x-{\hat{x}}^{\left(j\right)}|{|}_2/|\left|A^H(A{\hat{x}}^{\left(j\right)}-y)\right|{|}_2---\left(11\right)$

由式(11)和式(4)可知，的解依赖于的值，在利用式(11)计算参数值时是未知的，因此，的值可由在σ＝σ_j-1时的稀疏信号重构值来代替，

即，正则化参数

3.根据权利要求2所述的一种自适应正则化平滑l₀范数方法，其特征在于：的推导过程为，

计算A的SVD变换可得，

A＝UΣV^H (13)

其中，V∈C^m×m称为右奇异矩阵，

利用U、V的特性可得，

AA^H＝UΣΣ^H U^H (14)

UU^H＝I_m (15)

将式(12)中新的正则化参数、式(14)和(15)代入式(4)得，

$\begin{array}{c}{\hat{x}}^{\left(j\right)}=x-A^H{\left({\mathrm{U\Σ\Σ}}^H{U}^H+{\widehat{\mathrm{\λ}}}^{-1}{\mathrm{UU}}^H\right)}^{-1}\left(Ax-y\right)\\ =x-A^H{\[U\left({\mathrm{\Σ\Σ}}^H+{\widehat{\mathrm{\λ}}}^{-1}{I}_m\right)U^H\]}^{-1}\left(Ax-y\right)\\ =x-A^{H}U{\left({\mathrm{\Σ\Σ}}^H+{\widehat{\mathrm{\λ}}}^{-1}{I}_m\right)}^{-1}{U}^H\left(Ax-y\right)\end{array}---\left(16\right)$

其中，ΣΣ^H＝diag(α₁,α₂,…,α_m)，diag(·)为对角化操作，

令m×1维向量

则，

式中，表示Hadamard积，

那么式(16)可改写为，

4.根据权利要求3所述的一种自适应正则化平滑l₀范数方法，其特征在于：A的SVD变换在信息重构处理前进行离线计算。