CN110048723B - 鲁棒稀疏信号恢复方法 - Google Patents

Abstract

本发明提供了一种鲁棒稀疏信号恢复方法，基于相关系数提出的稀疏信号算法通过计算信号之间的相关系数，匹配出相关性强的信号位置以及强度，进而实现对入射信号的完整恢复。本发明的有益效果是：由于相关系数的鲁棒性，当信号受到脉冲噪声影响时，算法也能恢复出完整的信号。

Description

鲁棒稀疏信号恢复方法

技术领域

本发明涉及信号恢复方法，尤其涉及一种鲁棒稀疏信号恢复方法。

背景技术

稀疏信号恢复在压缩感知中是重要的一个环节，它在机器学习、模式识别、图像处理、阵列信号处理以及通信信号处理等领域中有着广泛的应用。因此，对稀疏信号恢复的研究有着重要的理论价值与工程实践前景。对采样后的信号或图像进行恢复，传统的方法要求采样率遵循奈奎斯特采样理论：采样率必须至少是信号频率的2倍。近年来，出现了一类新颖理论―压缩感知(Compressed sensing,CS)理论，在CS理论支撑下的稀疏恢复算法也越来越受到学术界及工程界的重视。在该理论框架下，对一定形式上满足稀疏性的数据，以低于奈奎斯特采样率的频率进行采样时，我们也能够通过特定的算法将原数据进行恢复或重构。所要求的稀疏形式可以表现在时域、频域、空时域等方面，并且可以是近似的稀疏。稀疏恢复算法，目前主要分为三类：凸优化方法，贪婪(greedy)算法，FOCUSS(focalundetermined system solver)算法。

凸优化算法的核心思想是用二阶锥规划方法逼近最稀疏解，算法特点是所需采样少、计算结果稳定，但计算量庞大；贪婪算法的核心思想是每步迭代中进行残差与采样距阵的匹配运算后剔除较小值以逼近稀疏最优解，算法特点是计算量少、收敛速度快，但所需的采样数比较多，典型算法有MP(Matching Pursuit)、OMP(Orthogonal Matching Pursuit)、CoSaMP(Compressive Sampling Matching Pursuit)等；FOCUSS采用加权最小2-范数方法，核心思想是每步迭代中令代价函数2-范数最小，算法特点是计算量小、收敛速度快，但是不能适用于多快拍的数据。

发明内容

为了解决现有技术中的问题，本发明提供了一种鲁棒稀疏信号恢复方法。

本发明提供了一种鲁棒稀疏信号恢复方法，基于相关系数提出的稀疏信号算法通过计算信号之间的/>相关系数，匹配出相关性强的信号位置以及强度，进而实现对入射信号的完整恢复。

作为本发明的进一步改进，一种鲁棒稀疏信号恢复方法，具体包括以下步骤：

S1、残差向量r∈R^M×1分别与感知矩阵A∈R^M×K的每一列计算相关系数，残差向量r∈R^M×1的初始值为观测向量r＝y，并设置误差容限ε，和标签集/>其中，R表示实数集，M表示向量r和矩阵A的行数为M，K表示矩阵A的列数为K；

S2、由求/>相关系数，其中，/>表示对残差向量r求/>范数之后再求p次方；表示寻找最优的(α,θ,z)使増广拉格朗日函数最小，并代入求出这一最小值；α是需要优化的参数；θ、z表示约束向量；/>表示对向量θ求/>范数之后再求p次方；λ₁、λ₂表示惩罚参数；/>表示对θ-z求二范数之后再平方；/>表示对(r-αa_i)-z求二范数之后再平方；a_i表示感知矩阵A的第i列；

S3、由求出/>相关系数，并在依次求出的K个/>相关系数看选出最大的一个/>相关系数，其中，/>表示对残差向量r求/>范数之后再求p次方；a_i表示感知矩阵A的第i列；

S4、取出最大的相关系数在所有K个/>相关系数中所在的位置i，并放入位置索引集合I中；

S5、集合I中包含的元素对应感知矩阵A的相应位置的列向量组成矩阵A_I，用来计算相应位置的信号强度，能使/>最小的u即为x_I，其中，表示对A_Iu-y求/>范数之后再求p次方；/>表示寻找使最小的向量u，并代入求出/>的最小值；

S6、更新残差值r＝y-A_Ix_I，然后继续重复步骤S1至S6，当时迭代结束，输出标签集I和相应的x_I，标签集I即表示恢复出信号的位置，x_I表示相应位置信号的强度，即完成了对稀疏信号的完整恢复。

作为本发明的进一步改进，在步骤S1中，计算相关系数时，主要按照以下步骤依次求出/>相关系数内的未知参数(θ、z、α)：

1.1、首先给予参数α和z赋初始值，再求参数θ，当p＝1时，通过对θ求导，使导数值为0时的θ即为满足条件的θ；当p≠1时，通过对θ求导，使导数值为0时的θ即为满足条件的θ，/>表示对向量θ求/>范数之后再求p次方；/>表示对向量θ-z求二范数之后再平方；R、θ表示约束向量；/>表示对向量Rθ求二范数之后再平方；λ₁为惩罚参数；/>表示对向量θ-z求二范数之后再平方；

1.2、参数z通过ADMM理论求出，具体求解过程为z＝λ₁(r-αa_i)+λ₂θ，而a_i是感知矩阵A的第i列向量，i＝1,2,…,K，λ₁、λ₂为惩罚参数；α是需要优化的参数；r表示残差向量；a_i表示感知矩阵A的第i列；θ表示约束向量；

1.3、继续由最小二乘法可以求解出参数α，具体表达式为a_i表示感知矩阵A的第i列；r表示残差向量；z表示约束向量；/>表示对向量a_i求转置；/>表示对/>求逆；

1.4、重复1.1至1.3的步骤，直至跳出循环，r表示残差向量；a_i表示感知矩阵A的第i列；/>表示对r-αa_i求二范数后再平方；

α是需要优化的参数；ε表示误差容限。

本发明的有益效果是：基于相关系数提出的稀疏信号算法通过计算信号之间的/>相关系数，匹配出相关性强的信号位置以及强度，进而实现对入射信号的完整恢复。由于/>相关系数的鲁棒性，当信号受到脉冲噪声影响时，算法也能恢复出完整的信号。

附图说明

图1是本发明一种鲁棒稀疏信号恢复方法的理想情况下的稀疏信号恢复示意图。

图2是本发明一种鲁棒稀疏信号恢复方法的混合高斯噪声下的稀疏信号恢复示意图。

具体实施方式

下面结合附图说明及具体实施方式对本发明作进一步说明。

本发明提出了一种基于相关系数的稀疏信号恢复方法。相关系数是一种常见的统计指标，是研究变量之间线性相关程度的量。相关系数为1时，说明变量之间绝对正相关；相关系数为-1时，说明变量之间绝对负相关。总之，相关系数的绝对值越靠近1，线性相关程度越高，若相关系数绝对值越接近于0，那么就说明线性相关程度越低。/>矩阵相关系数的提出主要是用于计算二维数据之间的相关性。当二维数据的某一行受到脉冲噪声的影响时，/>相关系数也能正确的计算出数据之间原有的相关程度。基于/>相关系数提出的稀疏信号算法通过计算信号之间的/>相关系数，匹配出相关性强的信号位置以及强度，进而实现对入射信号的完整恢复。由于/>相关系数的鲁棒性，当信号受到脉冲噪声影响时，算法也能恢复出完整的信号。

对于一个矩阵B＝[b₁,b₂,…,b_M]^T∈C^M×K，b₁，b₂，…，b_M表示M×1阶向量，“[·]^T”表示转置。首先可以定义计算出它的范数。计算时首先定义一个向量/>此向量/>的元素是矩阵B的每一行元素所组成的行向量做/>范数计算之后的结果。向量/>的表达式为符号“||·||₂”表示对向量计算/>范数，/>中共包含M个元素，第m个元素即为/>接着对于向量/>可计算出它的l_p范数。至此，矩阵B的/>范数就已通过计算得出，总结起来矩阵B的l_2,p范数的计算过程可以表示如下

“||·||_p”表示l_p范数，“||·||_2,p”表示范数。已知式(3-1)所示的矩阵/>范数的定义之后就不难提出/>相关系数的概念。若有两个矩阵A＝[a₁,a₂,…,a_M]^T,B＝[b₁,b₂,…,b_M]^T∈C^M×K，a₁，a₂，…，a_M也表示M×1阶向量，为了衡量这两个矩阵之间的相关程度，这里定义了两个矩阵之间的/>相关系数，表示如下

c_2,p(A,B)表示的是范数空间的/>相关系数。式(3-2)中，min(·)表示取到最小值，L_2,p将在式(3-4)给出详细表达式。若式min(L_2,p)值为0，那么/>相关系数θ_2,p(A,B)为1，说明这两个矩阵之间具有高的相关性；若min(L_2,p)值为||B||_2,p，则/>相关系数θ_2,p(A,B)为0，说明两个矩阵之间的相关性最弱。不难看出，解决式(3-2)中L_2,p的最小化问题是/>相关系数正确与否的关键，所以L_2,p的最小化过程我们重点考虑的问题。

式(3-1)中计算矩阵的范数时，首先要考虑的是矩阵每一行的/>范数，经过这一步之后的矩阵变成一维的向量/>进而可以对降维后的向量/>计算其/>范数。/>相关系数在计算两个矩阵间的相关性之外，也要保证对出现在两个矩阵间的异常值的鲁棒性。此时，我们认为两个矩阵之间的异常值是存在于矩阵的整行数据中，与/>相关系数中处理向量某一点出现异常值的情况不一样，/>相关系数考虑的异常值是某一行全部被其影响的情况。因而，此时不能简单的认为/>表示寻找最优的α值使/>取到最小值)，由于求解/>范数时首先计算了矩阵每一行的2范数，2范数对异常值是非常敏感的，所以这一步忽略了异常值带来的不利影响，直接使用此式缺少应用价值。对此，我们可先对矩阵A和B转置后再矢量化，即定义两个向量a＝vec(A^T)以及b＝vec(B^T)，vec(·)表示矢量化。然后，给出/>相关系数定义如下：

式中，λ₁+λ₂＝1。若要求解l_2,p相关系数，由定义的l_2,p相关系数的表达式可以看出关键在于求解最优的(α,θ,z)组合，使得L_2,p(α,θ,z)的值最小。

接下来我们要讨论相关系数求解方法，由式(3-3)可知求解/>相关系数的关键在于对下式

的最小化。由于在p＝2时，相关系数即为/>相关系数这种特殊情况，因此在求解/>相关系数时仅考虑0<p<2的情况。

而在0<p<2的范围内我们把求解两个矩阵间相关系数的算法分为p＝1与p≠1两种情况来讨论。L_2,p(α,θ,z)由三个参数共同决定，表达式包含三个部分，算法将按次顺求出三个参数，求第一个参数时，需要对另外两个参数设置初始值，然后通过多次迭代使得最终的三个参数是同时满足使L_2,p(α,θ,z)最小的最优解。

首先讨论p＝1的情况，算法先对θ进行求解，求θ时，我们认为L_2,p(α,θ,z)式中不包含参数θ的部分值为0。即此时要求解使

最小的θ，问题转化为式(3-5)对θ求导，并使L_2,p′(θ)＝0。

满足L_2,p′(θ)＝0的θ可由如下方法计算

max{·}表示取最大的值参与计算，对于式中出现的参数z，由L_2,p(α,θ,z)式中的第三部分的值为0可知z＝(b-αa)，而α可对其设置一个初始值。所以由式(3-6)，θ的值得以解出。

然后，可继续求解余下的参数z，并且已求出的参数θ也将用于以下对z的求解过程之中。由于θ已为定值，此时相当于求能使式最小的z。在机器学习中，ADMM(Alternating Direction Method of Multipliers)是一种广泛使用的约束问题最优方法，这种方法将无约束优化的部分用块坐标下降法来分别优化，是目前比较成熟，比较受欢迎的约束问题最优化通用框架。而ADMM算法在此处刚好可以应用于求解z的模型之中，由ADMM可得

z＝λ₁(b-αa)+λ₂θ (3-7)

最后再更新α的值，虽然在之前的过程中α作为设定的初始值参与了θ和z求解，但是初始设置的α显然不是最优解。同样在θ和z已求出的前提下，对α的求解即为寻找使最小的α值。相比于前两个参数的求解，此时α的计算相对简单，只需满足(b-αa)-z这一个表达式取到最小值即可，也就是令(b-αa)-z＝0后，从中求解α。而面临这种形式的方程时，不难得出运用最小二乘法(Least Square Method)就可以求出满足条件的α。因此，由最小二乘法我们可以由下式计算出α

α＝(a^Ta)^-1a^T(b-z) (3-8)

至此，L_2,p(α,θ,z)中的三个参数全部求出一轮，但由于在求解的过程中，赋给α的初始值较最优α值有一定的误差，因此需要用式(3-8)所求出的α继续迭代求解。迭代求解的目的是找出使式(3-4)最小的(α,θ,z)组合解，而当式(3-4)中b-αa这一部分值达到最小时，L_2,p(α,θ,z)值达到最小。在迭代的过程中，可设置一个误差容限值ε，当满足如下式所示的条件时，迭代结束。

在p＝1时，由上述迭代求解方法求出最小的L_2,p(α,θ,z)之后，再将其依次代入式(3-3)与式(3-2)中，即可求得矩阵A与B之间的相关系数，算法3-1给出了整个迭代求解的详细过程。

算法3-1：p＝1时的L_2,p相关系数求解过程：

初始化 α₀＝1.5，t＝0

计算 z₀＝(b-αa)

迭代

迭代次数 t←t+1

第一步，计算θ_t

第二步，通过ADMM法求z_t

z_t＝λ₁(b-αa)+λ₂θ

第三步，通过最小二乘法法求α

α＝(a^Ta)^-1a^T(b-z)

直到

接下来继续分析p≠1时的相关系数计算方法，求解的思路也是按照p＝1时的迭代顺序来寻找最优解，但具体的操作与p＝1时不同，且由式(3-4)的形式可知，算法的不同主要体现在对θ求解上。在(3-9)式中，向量a与b是由矩阵Α与B按行矢量化得来的，因此a与b是MK×1阶向量，即θ也是MK×1阶的向量，求解θ时我们首先随机生成一个初始的MK×1阶向量θ。考虑定义/>相关系数的目的是计算矩阵之间相关性，现将随机生成的MK×1阶向量θ重塑成M×K阶新矩阵

其中ζ_mK是1×K阶的列向量。向量θ的元素按照一行接一行的顺序依次填满这个M×K矩阵δ，借鉴MATLAB的程序的表达方式表示θ的第(m-1)K+1至mK个元素，则这一过程可以表示成ζ_mK＝θ(1+(m-1)K:1:mK),其中m＝1,2,…M。然后由范数的计算方法，对此矩阵δ的每一行组成的M个1×K阶行向量ζ_1K、ζ_2K、…,ζ_MK，要依次求出其/>范数，这样我们将得到一个M×1阶列向量/>其中第m个元素计算如下

接着对向量求/>范数，就完成了对式(3-4)中/>的求解。所以，在p≠1时，求解θ的式(3-5)可以写成如下形式

式中，m∈[1,2…,M]。

对于式(3-12)中的这一部分，因为此时p≠1，计算难度较大，若能将此处的p范数变形成为与第二部分/>一样的2范数时，问题将大为简化。因此，我们对做如下巧妙的变化：

变化之后的值并没有发生改变，只有形式上的变化。若定义一个M×1阶的向量d，且满足

显然，式(3-14)中有将其代入式(3-13)就可使式(3-13)形式更加简化，原来的指数为p，代入后指数变为常数2，即

又由2范数的定义可知

所以，结合式(3-15)至式(3-16)的推理，最初的式(3-12)通过这一系列的代换，可得如下的表达式

式中，符号“.”表示d与之间的点乘。

观察式(3-17)可知，向量外部的p范数通过公式变换已经成为更加常规的2范数，但是正因为寻求更加简洁的公式表达式进行了一系列变换，使得此时式中第一部分的表达式中没有要求解的向量θ，而我们的算法思想正是通过对θ求导进而来寻找使式(3-4)最小的θ值，所以还需要对(3-17)式做进一步推导。

对于式中M×1阶向量如式(3-11)所述有/>所以式(3-17)中出现的/>可以表示如下

式中k∈[1,2,…,K]。

式(3-18)中d_m是是向量d中的第m个元素，且式(3-18)中最后一步d_m与向量ζ_mK相乘之后再求和。这一步也可以这样理解：已知ζ_mK是1×K阶向量，将d_m自身重复扩展为1×K阶的向量，再与向量ζ_mK点乘，对点乘之后的1×K阶向量求和。而这样理解的好处是保证了尺寸上的一致。

因此可定义一个1×K阶向量c_mK，使

式中1_1×K表示1×K阶单位向量，符号表示克罗内克积，m∈[1,2…,M]。由克罗内克积的定义可知，c_mK包含的K个元素相同且都是d_m，即c_mk＝d_m。因此，式(3-18)又可以写成

又考虑到向量d中有M个元素，M个元素依次与1_1×K进行克罗内克积运算后，将有M个1×K阶向量c_1K、c_2K、…c_MK，若将这M个向量组合起来，则有M×K阶矩阵C＝[c_1K,c_2K,…c_MK]^T。矩阵C实质上就是由下式计算得出

结合式(3-20)，将推出的的关系代入/>有其中/>相当于矩阵δ＝[ζ_K1,ζ_K2,…,ζ_KM]^T的元素平方后先按行求和，然后按列求和，那么，/>也就表示对向量θ的元素平方后再求和，也即求/>同理，/>即为对矩阵C的元素平方后先按行求和，然后按列求和。若设MK×1阶的向量r＝vec(C^T)，那么/>也表示对向量r的元素平方后再求和。所以，此时式(3-17)可表示为

式中，符号“.”表示θ与r之间的点乘。式(4-27)可进一步写成

式中，矩阵R表示以向量r为对角线元素的对角矩阵，即R＝diag(r)。

对于式(3-23)，目前形式上与式(3-5)相当，因此接下来即可按照表3-1中p＝1时的相关系数方法来继续求解。

对于稀疏信号的恢复可以用-OMP算法实现，/>-OMP算法是OMP算法在/>范数空间的应用。

下面将介绍-OMP算法过程。考虑在信号处理过程中，经常会遇到这样一个模型：

y＝Ax+v (3-24)

而-OMP算法正是要从这一模型中恢复未知信号x。

式(3-24)中，向量y∈R^M是观测向量(observation vector)。x＝[x₁,x₂,…,x_K]^T∈R^K是真实信号，包含要求解的信息。向量v∈R^M是信号观测或传输时产生的噪声误差。

矩阵A＝[a₁,a₂,…,a_K]∈R^M×K(M<K)是系数矩阵，也称作过完备字典矩阵，且矩阵A的每一列均是正则化的，即||a_k||₂＝1,k＝1,2,…K。该模型的目标就是通过观测值y和矩阵A，把未知向量x求出来。

-OMP算法的主要过程可以理解成输入观测数据向量y和字典矩阵A后，通过内部迭代，再输出稀疏系数向量x。/>-OMP算法首先要初始化一个标签集/>同时定义初始残差向量r⁰＝y，并令t＝1，t是对迭代次数的统计。详细迭代过程如下：

第一步，求矩阵A中与残差向量r^t-1的相关性最高的列，找出该列的位置索引i_t，i的取值满足i＝1,2,…M；

i_t＝argmaxC_2,p(a_i,r^t-1) (3-25)

并取出该列在矩阵A中的位置索引标签i_t，并放入标签集I中。

I^t＝I^t-1∪i_t (3-26)

第二步，求解最小化的问题

其中表示标签集I^t内索引位置对应的矩阵A中的列向量组成的新矩阵。

第三步，由式(3-27)求出的结果，可更新残差值；

残差值在不断的迭代过程中会越来越小，在理想情况下，残差值最终结果为0。

第四步，令t＝t+1，再重复第一步至第三步的整个过程，直到满足停止循环迭代的条件。最后输出系数向量和标签集I^t，即完成了对稀疏信号x的完整恢复。

在-OMP算法时，通常有几种判断迭代终止的方法。可以运行到某个固定的迭代步数后停止，而这个步数一般要大于稀疏系数向量x中的非零值的总个数。但是更为严谨的是预先给定一个极小的残差容限值ε，当迭代过程中求出的残差值小于ε时则可停止迭代，即

也可停止迭代。

本发明提供的一种鲁棒稀疏信号恢复方法，具体包括以下步骤：

1、残差向量r∈R^M×1分别与感知矩阵A∈R^M×K的每一列的计算相关系数，残差向量r∈R^M×1的初始值为观测向量r＝y，并设置误差容限ε，和标签集/>其中，R表示实数集，M表示向量r和矩阵A的行数为M，K表示矩阵A的列数为K。计算/>相关系数时主要按照以下步骤依次求出/>相关系数内的未知参数(θ、z、α)。

1.1、首先给予参数α和z赋初始值，再求参数θ，当p＝1时，通过对θ求导，使导数值为0时的θ即为满足条件的θ；当p≠1时，通过对θ求导，使导数值为0时的θ即为满足条件的θ，/>表示对向量θ求/>范数之后再求p次方；/>表示对向量θ-z求二范数之后再平方；R、θ表示约束向量；表示对向量Rθ求二范数之后再平方；λ₁为惩罚参数；/>表示对向量θ-z求二范数之后再平方。

1.2、参数z可通过ADMM理论求出，具体求解过程为z＝λ₁(r-αa_i)+λ₂θ，而a_i是感知矩阵A的第i列向量，i＝1,2,…,K，λ₁、λ₂为惩罚参数；α是需要优化的参数；r表示残差向量；a_i表示感知矩阵A的第i列；θ表示约束向量。

1.3、继续由最小二乘法可以求解出参数α，具体表达式为a_i表示感知矩阵A的第i列；r表示残差向量；z表示约束向量；/>表示对向量a_i求转置；/>表示对/>求逆。

1.4、重复1.1至1.3的步骤，直至跳出循环，r表示残差向量；a_i表示感知矩阵A的第i列；/>表示对r-αa_i求二范数后再平方；α是需要优化的参数；ε表示误差容限。

2、在三个参数全部求出之后，即可由求/>相关系数，其中，/>表示对残差向量r求/>范数之后再求p次方；/>表示寻找最优的(α,θ,z)使増广拉格朗日函数/>最小，并代入求出这一最小值；α是需要优化的参数；θ、z表示约束向量；/>表示对向量θ求/>范数之后再求p次方；λ₁、λ₂表示惩罚参数；/>表示对θ-z求二范数之后再平方；/>表示对(r-αa_i)-z求二范数之后再平方；a_i表示感知矩阵A的第i列。

3、相关系数求出之后，即可由/>求出/>相关系数，并在依次求出的K个/>相关系数选出最大的一个/>相关系数，其中，/>表示对残差向量r求范数之后再求p次方；a_i表示感知矩阵A的第i列。

4、取出最大的相关系数在所有K个/>相关系数中所在的位置i，并放入位置索引集合I中。

5、集合I中包含的元素对应感知矩阵A的相应位置的列向量组成矩阵A_I，用来计算相应位置的信号强度，能使/>最小的u即为x_I，其中，表示对A_Iu-y求/>范数之后再求p次方；/>表示寻找使最小的向量u，并代入求出/>的最小值。

6、更新残差值r＝y-A_Ix_I。然后继续重复第1步至第6步。当时迭代结束。输出标签集I和相应的x_I，标签集I即表示恢复出信号的位置，x_I表示相应位置信号的强度。

实验中将用混合高斯模型来模拟脉冲噪声。我们可以定义q来表示整个混合高斯噪声过程中，冲激出现的频率。并设脉冲噪声的冲激强度为P_o，显然，q越大，噪声中冲激的出现越频密，而P_o的值越大，则噪声中冲激的强度越强。表示一般零均值高斯白噪声的方差，那么对于冲激部分则服从零均值方差为/>的高斯分布。对M×K阶的矩阵A加混合高斯分布的噪声，我们首先要随机生成一个元素取值在[0，1]内的M×1阶的向量γ，若γ中第m个元素的值大于等于q，则矩阵A的第m行将叠加方差为/>高斯白噪声，若第m个元素的值小于q，则矩阵A的第m行将叠加方差为/>的冲激噪声。若矩阵A的能量由/>表示，那么信噪比定义如下

在实验中我们设置需要重建的稀疏信号为一100×1阶向量x，且x中有T＝10个非零量。x中10个非零量出现的位置是随机的，非零量的数值大小[-2,-1]∪[1,2]的范围内随机生成。

而对x的观测由一个随机生成的服从高斯分布的矩阵A来感知，称A为感知矩阵，同时设置A是150×100阶的矩阵。在设置矩阵A与向量x都设置好之后，可以得到所观测的信息为y＝Ax+v，y称为观测向量，其中向量v代表混合高斯噪声。在本实验中，同样设置混合高斯噪声模型中的q＝0.2，P_o＝100，设信噪比为2dB。实验中，Ax是150×1阶的向量，在叠加脉冲噪声时，将Ax看成是30×5阶的矩阵，即每连续5个元素组成一行，再对其叠加噪声，叠加混合高斯噪声之后的矩阵转置后矢量化成向量之后就可以开始应用-OMP算法对x进行的恢复。实验将/>-OMP算法做对比，分无噪的理想情况和叠加混合高斯噪声两种情况讨论。分别验证了p＝1、p＝1.2和p＝0.7这三种取值下的稀疏信号恢复情况。

图1与2中圆圈表示真实信号的，星号表示算法恢复出的信号；图1中在理想情况下，-OMP算法和/>-OMP算法都能正确恢复出真实信号x，图中圆圈的位置和高度与星号代表的恢复出的信号完全一致。而在混合高斯噪声情况下，如图2所示，/>-OMP算法在p的三种取值下，恢复出的信号与真实信号不管是位置还是强度都发生了严重的偏离，而/>-OMP此时依然能保证恢复出的信号与原真实信号的一致。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (2)

1.一种鲁棒稀疏信号恢复方法，其特征在于，基于l_2,p相关系数提出的稀疏信号算法通过计算信号之间的l_2,p相关系数，匹配出相关性强的信号位置以及强度，进而实现对入射信号的完整恢复；

所述方法包括以下步骤：

S1、残差向量r∈R^M×1分别与感知矩阵A∈R^M×K的每一列计算l_2,p相关系数，残差向量r∈R^M×1的初始值为观测向量y＝残差向量r，并设置误差容限ε，和标签集其中，R表示实数集，M表示向量r和矩阵A的行数为M，K表示矩阵A的列数为K；

S2、由求l_2,p相关系数，其中，c_2,p(a_i,r)表示的是l_2,p范数空间的l_2,p相关系数，/>表示对残差向量r求l_2,p范数之后再求p次方；/>表示寻找最优的(α,θ,z)使増广拉格朗日函数/>最小，并代入求出这一最小值；α是需要优化的参数；θ、z表示约束向量；/>表示对向量θ求l_2,p范数之后再求p次方；λ₁、λ₂表示惩罚参数；表示对θ-z求二范数之后再平方；/>表示对(r-αa_i)-z求二范数之后再平方；a_i表示感知矩阵A的第i列；

S3、由求出l_2,p相关系数，并在依次求出的K个l_2,p相关系数看选出最大的一个l_2,p相关系数，其中，θ_2,p(a_i,r)表示l_2,p相关系数，/>表示对残差向量r求l_2,p范数之后再求p次方；a_i表示感知矩阵A的第i列；

S4、取出最大的l_2,p相关系数在所有K个l_2,p相关系数中所在的位置i，并放入位置索引集合I中；

S5、集合I中包含的元素对应感知矩阵A的相应位置的列向量组成矩阵A_I，用来计算相应位置的信号强度，能使/>最小的u即为x_I；其中，表示对A_Iu-y求l_2,p范数之后再求p次方；/>表示寻找使最小的向量u，并代入求出/>的最小值；

S6、更新残差值r＝y-A_Ix_I，然后继续重复步骤S1至S6，当时迭代结束，输出标签集I和相应的x_I，标签集I即表示恢复出信号的位置，x_I表示相应位置信号的强度，即完成了对稀疏信号的完整恢复。

2.根据权利要求1所述的鲁棒稀疏信号恢复方法，其特征在于，在步骤S1中，计算l_2,p相关系数时，主要按照以下步骤依次求出l_2,p相关系数内的未知参数(θ、z、α)：

1.1、首先给予参数α和z赋初始值，再求参数θ，当p＝1时，通过对θ求导，使导数值为0时的θ即为满足条件的θ；当p≠1时，通过/>对θ求导，使导数值为0时的θ即为满足条件的θ；/>表示对向量θ求l_2,p范数之后再求p次方；表示对向量θ-z求二范数之后再平方；R、θ表示约束向量；/>表示对向量Rθ求二范数之后再平方；λ₁为惩罚参数；/>表示对向量θ-z求二范数之后再平方；

1.2、参数z通过ADMM理论求出，具体求解过程为z＝λ₁(r-αa_i)+λ₂θ，而a_i是感知矩阵A的第i列向量，i＝1,2,…,K；λ₁、λ₂为惩罚参数；α是需要优化的参数；r表示残差向量；a_i表示感知矩阵A的第i列；θ表示约束向量；

1.3、继续由最小二乘法可以求解出参数α，具体表达式为a_i表示感知矩阵A的第i列；r表示残差向量；z表示约束向量；/>表示对向量a_i求转置；/>表示对/>求逆；

1.4、重复1.1至1.3的步骤，直至跳出循环，r表示残差向量；a_i表示感知矩阵A的第i列；/>表示对r-αa_i求二范数后再平方；α是需要优化的参数；ε表示误差容限。