CN104392412B - 基于进化正交匹配追踪的压缩感知信号恢复方法 - Google Patents

Abstract

本发明公开了一种基于进化正交匹配追踪的信号恢复方法，主要解决压缩感知中传统追踪算法过于贪婪，回溯能力差和恢复准确率低的问题。其技术方案是：将进化计算的框架引入到压缩感知信号恢复当中；将原子选择的问题转化为基于启发式搜索的种群寻优的过程；结合传统贪婪追踪算法中观测误差与原子的相关性，定义了一种活性函数来度量每个原子被选择的可能性；通过活性函数，设计出了弱贪婪的交叉和变异算子，从而使得更多的原子有可能被选择，增加了信号恢复中原子搜索的可达空间。实验表明，对于信号的压缩感知恢复，本发明比传统的贪婪追踪算法有更高的恢复概率和更小的恢复误差，可用于一维信号和二维图像信号在低采样率随机观测下的恢复问题。

Description

基于进化正交匹配追踪的压缩感知信号恢复方法

技术领域

本发明属于信号处理技术领域，特别涉及一种压缩感知信号恢复方法。可用于一维信号和二维图像信号在低采样率随机观测下的恢复问题。

背景技术

压缩感知是一种高效的信号获取手段，它仅通过少量的观测即可恢复出原始信号。压缩感知的基本假设是：信号可以被某组基原子稀疏地线性表示，即大多数原子系数为零，只有极少数原子系数为非零。在这种假设下，通过满足一定条件的观测，即可准确地恢复出信号。由于压缩感知是非凸欠定问题，很难给出一个闭合解形式。解决这类问题最常用的技术是贪婪追踪，经典的算法有匹配追踪、正交匹配追踪、子空间追踪、分段正交匹配追踪和正则化正交匹配追踪。

正交匹配追踪算法的基本思想是：首先初始化观测误差为观测向量，在每次迭代中，计算观测误差与压缩感知矩阵中每个原子的内积；然后，选择对应内积绝对值最大的原子到支撑集当中去，并用最小二乘法计算支撑集中原子的系数；最后，更新观测误差，进入下一次迭代。通过这个迭代过程，使观测误差逐步减少，实现对原始信号的逼近。

子空间追踪算法由正交匹配追踪算法改进而来，它引入了回溯技术来淘汰掉不可信的原子，其基本思想为：假设原子非零系数的个数为K，首先，初始化观测误差为观测向量，并计算观测误差与压缩感知矩阵中每个原子的内积，将内积绝对值最大的K个原子添加到支撑集中去，然后利用最小二乘法求这K个原子的系数，并更新观测误差；在每次迭代中，再次计算观测误差与压缩感知矩阵中每个原子的内积，并且将内积绝对值最大的K个原子添加到支撑集中去，这样得到了2K个原子的支撑集；然后，利用最小二乘法求解这2K个原子的系数；之后，将支撑集中2K个原子系数绝对值最大的K个保留，将剩下的原子删除，这样得到了一个K个原子的支撑集；再利用最小二乘法求解这K个原子的系数；最后，更新观测误差，进入下一次迭代。

虽然许多基于正交匹配追踪的改进算法对于信号的恢复概率和误差获得了很大的提升，但是由于匹配追踪算法自身过于贪婪和搜索范围有限等局限性，经常会使搜索陷入局部最优。

发明内容

本发明的目的在于改善传统追踪算法过于贪婪，回溯能力差的问题，来提高信号的恢复概率和恢复精度。是一种基于进化正交匹配追踪的压缩感知信号恢复方法，用于对一维信号和二维图像信号在低采样率随机观测下的恢复问题。本发明在进化计算的框架下来求解压缩感知信号恢复问题，根据贪婪算法中观测误差与原子的相关性大小引入了基因活性的概念，利用基因活性作为整个种群的启发式知识，来操作不同个体之间的交叉和变异，从而使得整个种群能够不断地迭代寻优。

本发明的技术方案是：首先初始化父代种群，计算其中所有个体的适应度和基因活性，并记录下最优个体。然后对父代种群中的个体随机两两配对，根据基因活性对每一对个体实施交叉操作，将得到的新个体存入子代种群，并计算它们的适应度和基因活性。之后，根据整体基因活性对子代种群中的每一个个体进行变异，每变异一个个体都要及时更新其适应度、基因活性和整体基因活性。最后，根据适应度好坏，在父代和子代中选择出新的父代种群到下一次迭代中，直到迭代停止。其具体步骤包括如下：

(1)输入压缩感知矩阵D_cs和观测向量y，初始化一个含有S个个体的父代种群P_f＝{p_i}_1≤i≤S，置计数器t＝0；

(2)计算父代种群的适应度和基因活性，并且记录适应度最大的个体为最优个体p_best；

(3)对父代种群进行交叉操作，得到子代种群

(4)计算子代种群的适应度和基因活性；

(5)对子代种群进行变异操作，变异过程中更新每个个体的适应度、基因活性和整体基因活性；

(6)从父代种群和子代种群中选择出新的父代种群并且更新最优个体p_best，置计数器t＝t+1；

(7)设最大迭代次数为T_max，若t<T_max，返回步骤(3)；否则，输出最优结果。

本发明结合了进化计算和贪婪追踪技术，将压缩感知中稀疏系数的求解转化为弱贪婪的启发式种群搜索过程，并且引入贪婪算法中观测误差与原子的相关性作为种群搜索的启发式知识，它具有如下优点：

(A)、基于群体的启发式搜索改善了传统贪婪追踪算法陷入局部最优的问题；

(B)、对于一维信号的恢复概率和恢复精度有明显的提高；

(C)、对于二维图像的恢复，具有更小的恢复误差和更少的块效应。

实验证明，对于一维模拟信号的恢复，本发明比传统的贪婪追踪算法有更高的恢复概率和更小的恢复误差。对于二维图像的恢复，本发明比传统的贪婪追踪算法具有更高峰值信噪比PSNR以及更小的恢复误差。

附图说明

图1是本发明的整体实现流程图；

图2是本发明与其它方法对于系数服从标准正太分布的一维信号恢复效果对比；

图3是本发明与其它方法对于系数服从-1到1均匀分布的一维信号恢复效果对比；

图4是本发明与其它方法对于二维图像Lena的恢复效果对比；

图5是本发明与其它方法对于二维图像Lena的恢复误差对比；

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤1，初始化一个含有16个个体的父代种群

1.1)输入压缩感知矩阵D_cs与观测向量y，置父代种群为空集

1.2)计算压缩感知矩阵D_cs中原子与观测向量y的相关性c：

其中，|·|表示求向量元素的绝对值；

1.3)将向量c中的元素由大到小排序：找到一组索引序列{λ₁,λ₂,…,λ_N}，使得c[λ₁]≥c[λ₂]≥…≥c[λ_N]，其中c[λ_n]表示向量c的第λ_n个元素，N是信号维数；

1.4)生成一个全0个体p_i＝[0,0,…,0]，置p_i[λ_i]＝1，其中p_i[λ_i]表示个体p_i的第λ_i个基因；

1.5)将个体p_i添加到父代种群中置i＝i+1；

1.6)若i≤16，返回步骤1.4；否则，输出父代种群

步骤2,计算父代种群的适应度和基因活性，并且记录适应度最大的个体为最优个体p_best。

2.1)对于父代种群中的每个个体p，求解其稀疏系数向量α，按如下公式计算：

其中α_p表示以个体p中基因为1的位置为索引从向量α中抽取的子向量，表示以个体p中基因为0的位置为索引从向量α中抽取的子向量，表示以个体p中基因为1的位置为索引从矩阵D_cs中抽取的列向量构成的子矩阵，上标表示求矩阵的伪逆，0表示一个元素全为0的向量；

2.2)求个体p的适应度值，按如下公式计算：

其中，||·||₂为计算向量的2范数；

2.3)计算个体p的基因活性，按如下公式计算：

其中U是归一化常数，等于中元素的最大值；

2.4)将适应度值最大的个体记录为最佳个体p_best。

步骤3，对父代种群进行交叉操作，得到子代种群

3.1)设子代种群为空集：将父代种群中的个体随机两两配对，得到8对个体(p_i,p_j)；

3.2)对每一对个体(p_i,p_j)，找到它们同一位置下基因不同的索引集合Λ：

其中，p_i[λ]表示个体p_i的第λ位基因；

3.3)找到个体p_i和p_j在索引Λ中基因活性最大值的位置和

其中，l_i[λ]表示个体p_i第λ位基因的活性值；

3.4)交换个体p_i和p_j的基因位:如果置得到新个体p_i′，置得到新个体p′_j；否则，置得到新个体p_i′，置得到新个体p′_j；

3.5)将两个新个体p_i′和p′_j添加到子代种群重复步骤3.2到3.5操作，直到所有个体完成交叉。

步骤4，计算子代种群的适应度和基因活性。

对子代种群中的所有个体，按照步骤2.1到2.3计算它们的适应度和基因活性。

步骤5，对子代种群进行变异操作，变异过程中更新每个个体的适应度、基因活性和整体基因活性。

5.1)置i＝1；

5.2)计算子代种群的平均稀疏度K_mean，按如下公式计算：

其中p[j]代表个体p的第j个基因，N是信号长度；

5.3)计算子代种群的整体基因活性，按如下公式计算：

其中l_j表示第j个个体的基因活性；

5.4)计算个体p_i的稀疏度K_i：

5.5)判断个体p_i的稀疏度K_i与种群的平均稀疏度K_mean的大小，并根据大小关系进行变异操作；

5.5.1)如果K_i≤K_mean，记录个体p_i中基因值为0的索引集合：

Λ＝{λ|p_i[λ]＝0,1≤λ≤N}，

找到整体基因活性L在索引集合Λ中最大值的位置λ_max：

置个体p_i的第λ_max个基因为1，即p_i[λ_max]＝1；

5.5.2)如果K_i>K_mean，记录个体p_i中基因值为1的索引集合：

Λ＝{λ|p_i[λ]＝1,1≤λ≤N}，

找到整体基因活性L在索引集合Λ中最小值的位置λ_min：

置个体p_i的第λ_min个基因为0，即p_i[λ_min]＝0；

5.6)按照步骤2.1)到2.3)更新个体p_i的适应度值和基因活性，置i＝i+1，若i≤16，返回步骤5.2)；否则，输出变异后的子代种群

步骤6，从父代种群和子代种群中选择出新的父代种群并且更新最优个体p_best，置计数器t＝t+1。

在父代种群和子代种群的集合中，选择16个适应度值最大的个体组成新的父代种群将新的父代种群中的适应度值最大的个体更新为当前最优个体p_best，并且置计数器t＝t+1。

步骤7，设最大迭代次数为T_max＝100，若t<T_max，返回步骤3；否则，输出最优结果。

本发明的效果可以通过仿真实验具体说明：

1.实验条件

实验所用微机CPU为Intel Core22.66GHz内存4GB，编程平台是Matlab2013。实验中采用的图像数据为标准图像测试库中的Lena图像，大小为512×512。

2.实验内容

实验1对非零系数服从标准正态分布的一维模拟信号进行恢复：

设信号长度N＝256，观测维数M＝100，稀疏系数向量α的稀疏度从10增加到50，每个稀疏度做1000次独立的实验，每次实验中，系数向量α的非零系数用标准正太分布随机产生，压缩感知矩阵D_cs的元素用均值为0，标准差为1/N的高斯分布产生。然后观测不同方法在1000次实验中正确恢复信号的概率和平均恢复误差。用本发明和4种经典贪婪追踪算法对上述仿真信号进行恢复，4种经典贪婪追踪算法分别为：匹配追踪MP、正交匹配追踪OMP、分段正交匹配追踪StOMP和子空间追踪SP。实验结果如图2所示，其中：

图2(a)为本发明和所述4种经典贪婪追踪算法在不同稀疏度下，对非零系数服从标准正态分布的信号准确恢复的概率对比曲线，

图2(b)为本发明和所述4种经典贪婪追踪算法在不同稀疏度下，对非零系数服从标准正态分布的信号平均恢复误差的对比曲线。

实验2对非零系数服从-1到1均匀分布的一维模拟信号进行恢复：

设信号长度N＝256，观测维数M＝100，稀疏系数向量α的稀疏度从10增加到50，每个稀疏度做1000次独立的实验，每次实验中，系数向量α的非零系数用-1到1的均匀分布随机产生，压缩感知矩阵D_cs的元素用均值为0，标准差为1/N的高斯分布产生。然后观测不同方法在1000次实验中正确恢复信号的概率和平均恢复误差。用本发明和4种经典贪婪追踪算法对上述仿真信号进行恢复，4种经典贪婪追踪算法分别为：匹配追踪MP、正交匹配追踪OMP、分段正交匹配追踪StOMP和子空间追踪SP。实验结果如图3所示，其中：

图3(a)为本发明和所述4种经典贪婪追踪算法在不同稀疏度下，对非零系数服从-1到1均匀分布的信号准确恢复的概率对比曲线，

图3(b)为本发明和所述4种经典贪婪追踪算法在不同稀疏度下，对非零系数服从-1到1均匀分布的信号平均恢复误差的对比曲线。

从图2和图3中可以看出，对于一维仿真信号，相比于其它4种经典的贪婪追踪算法，本发明具有更高的恢复概率和更小的恢复误差。

实验3对二维图像Lena进行分块压缩感知恢复,比较不同算法的峰值信噪比PSNR：

将Lena图像分成8×8的小块，设置观测维数M＝32，最大迭代次数为16。观测矩阵Φ中的元素用均值为0，标准差为1/N的高斯分布产生，基字典Ψ采用DCT字典，压缩感知矩阵D_cs＝ΦΨ。观察不同方法的恢复结果，并且比较它们的峰值信噪比PSNR。用本发明和现有匹配追踪MP与正交匹配追踪OMP算法对图像Lena进行恢复。实验结果如图4所示，其中：

图4(a)为原始的Lena局部图像，

图4(b)为本发明对Lena局部图像的恢复结果，峰值信噪比PSNR＝31.3dB，

图4(c)为MP算法对Lena局部图像的恢复结果，峰值信噪比PSNR＝29.2dB，

图4(d)为OMP算法对Lena局部图像的恢复结果，峰值信噪比PSNR＝29.6dB。

从图4可以看出，对于二维图像信号的恢复，本发明恢复结果的边界最接近原图，峰值信噪比PSNR最高，并且人工块效应最小。

实验4对二维图像Lena进行分块压缩感知恢复,比较不同算法的恢复误差：

将Lena图像分成8×8的小块，设置观测维数M＝32，最大迭代次数为16。观测矩阵Φ中的元素用均值为0，标准差为1/N的高斯分布产生，基字典Ψ采用DCT字典，压缩感知矩阵D_cs＝ΦΨ。观察不同方法的恢复误差图。比较本发明和子空间追踪SP、匹配追踪MP与正交匹配追踪OMP算法的恢复误差。实验结果如图5所示，其中：

图5(a)为SP算法对Lena局部图像的恢复误差图，

图5(b)为本发明对Lena局部图像的恢复误差图，

图5(c)为MP算法对Lena局部图像的恢复误差图，

图5(d)为OMP算法对Lena局部图像的恢复误差图。

从图5中的误差图来看，本发明恢复结果的误差幅度是最小的。

综上，无论在一维模拟信号还是二维图像信号的恢复中，本发明的恢复效果均优于其它4中经典的贪婪追踪算法。

Claims (5)

1.一种基于进化正交匹配追踪的压缩感知信号恢复方法，其流程如下：

(1)输入压缩感知矩阵D_cs和观测向量y，初始化一个含有S个个体的父代种群置计数器t＝0；

(2)计算父代种群的适应度和基因活性，并且记录适应度最大的个体为最优个体p_best，计算基因活性的公式如下：

$l=\frac{\left|D_{cs}^T(y-D_{cs}\mathrm{\&alpha;})\right|}{U}$

其中U是归一化常数，等于中元素的最大值，y表示观测向量，D_cs表示压缩感知矩阵，α表示稀疏系数向量；

(3)对父代种群进行交叉操作，得到子代种群

(4)计算子代种群的适应度和基因活性；

(5)对子代种群进行变异操作，变异过程中更新每个个体的适应度、基因活性和整体基因活性，计算整体基因活性的公式如下：

$L=\underset{j=1}{\overset{S}{\&Sigma;}}{l}_j$

其中l_j表示第j个个体的基因活性，S表示种群中含有个体的数量；

(6)从父代种群和子代种群中选择出新的父代种群并且更新最优个体p_best，置计数器t＝t+1；

(7)设最大迭代次数为T_max，若t＜T_max，返回步骤(3)；否则，输出最优结果。

2.根据权利要求1所述的基于进化正交匹配追踪的压缩感知信号恢复方法，其步骤(1)所述的输入压缩感知矩阵D_cs和观测向量y，初始化一个含有S个个体的父代种群置计数器t＝0具体步骤如下：

(1.1)输入压缩感知矩阵D_cs和观测向量y；

(1.2)置父代种群为空集置i＝1；

(1.3)计算压缩感知矩阵D_cs中原子与观测向量y的相关性c:

$c=\left|D_{cs}^{T}y\right|$

其中，|·|表示求向量元素的绝对值；

(1.4)将向量c中的元素由大到小排序：即找到一组索引{λ₁,λ₂,...λ_n...,λ_N}，使得c[λ₁]≥c[λ₂]≥...c[λ_n]...≥c[λ_N]，其中c[λ_n]表示向量c的第λ_n个元素，N是信号维数；

(1.5)生成一个全0个体p_i＝[0,0,...,0]，置p_i[λ_i]＝1，其中p_i[λ_i]表示个体p_i的第λ_i个元素；

(1.6)将个体p_i添加到父代种群中置i＝i+1；

(1.7)设S为父代种群中的个体数量，若i≤S，返回(1.5)；否则，输出父代种群置计数器t＝0。

3.根据权利要求1所述的基于进化正交匹配追踪的压缩感知信号恢复方法，其步骤(2)所述的计算父代种群的适应度和基因活性，并且记录适应度最大的个体为最优个体p_best具体步骤为：

(2.1)对于父代种群中的每个个体p，求解其稀疏系数向量α：

其中α_p表示以个体p中元素为1的位置为索引从向量α中抽取的子向量，表示以个体p中元素为0的位置为索引从向量α中抽取的子向量，表示以个体p中元素为1的位置为索引从矩阵D_cs中抽取的列向量构成的子矩阵，上标表示求矩阵的伪逆，0表示一个元素全为0的向量；

(2.2)求个体p的适应度值：

$f=\frac{1}{\left|\right|y-D_{cs}\mathrm{\&alpha;}|{|}_2}$

其中，||·||₂为计算向量的2范数；

(2.3)计算个体p的基因活性：

$l=\frac{\left|D_{cs}^T(y-D_{cs}\mathrm{\&alpha;})\right|}{U}$

其中U是归一化常数，等于中元素的最大值。

4.根据权利要求1所述的基于进化正交匹配追踪的压缩感知信号恢复方法，其步骤(3)所述的对父代种群进行交叉操作的方法如下：

(3.1)置子代种群为空集将父代种群中的个体两两配对，得到S/2对个体；

(3.2)对每一对个体(p_i,p_j)，找到它们同一位置下基因不同的索引集合Λ，即：

$\mathrm{\&Lambda;}=\left\{\mathrm{\&lambda;}\right|p_i\&lsqb;\mathrm{\&lambda;}\&rsqb;-p_j\&lsqb;\mathrm{\&lambda;}\&rsqb;\&NotEqual;0,\&ForAll;1\&le;\mathrm{\&lambda;}\&le;N\}$

(3.3)找到个体p_i和p_j在索引Λ中基因活性最大值的位置和即:

$\begin{array}{cc}{l}_i\&lsqb;{\mathrm{\&lambda;}}_i^{max}\&rsqb;\&GreaterEqual;l_i\&lsqb;\mathrm{\&lambda;}\&rsqb;& \&ForAll;{\mathrm{\&lambda;}}_i^{max},\mathrm{\&lambda;}\&Element;\mathrm{\&Lambda;}\end{array}$

$\begin{array}{cc}{l}_j\&lsqb;{\mathrm{\&lambda;}}_j^{max}\&rsqb;\&GreaterEqual;l_j\&lsqb;\mathrm{\&lambda;}\&rsqb;& \&ForAll;{\mathrm{\&lambda;}}_j^{\max },\mathrm{\&lambda;}\&Element;\mathrm{\&Lambda;}\end{array}$

其中l_i[λ]表示个体p_i第λ位的基因活性值；

(3.4)交换个体p_i和p_j的基因位：如果置得到新个体p′_i，置得到新个体p′_j；否则，置得到新个体p′_i，置得到新个体p′_j；

(3.5)将两个新个体p′_i和p′_j添加到子代种群重复步骤(3.2)到(3.5)操作，直到所有个体完成交叉。

5.根据权利要求1所述的基于进化正交匹配追踪的压缩感知信号恢复方法，其步骤(5)所述的对子代种群进行变异操作，过程如下：

(5.1)置i＝1；

(5.2)计算子代种群的平均稀疏度K_mean：

$K_{mean}=\frac{1}{S\&times;N}\underset{p\&Element;P_s}{\&Sigma;}\underset{j=1}{\overset{N}{\&Sigma;}}p\&lsqb;j\&rsqb;$

其中p[j]代表个体p的第j个基因，N是信号长度；

(5.3)计算子代种群的整体基因活性，公式如下：

$L=\underset{j=1}{\overset{S}{\&Sigma;}}{l}_j$

其中l_j表示第j个个体的基因活性；

(5.4)计算个体p_i的稀疏度K_i，公式如下：

$K_i=\underset{j=1}{\overset{N}{\&Sigma;}}{p}_i\&lsqb;j\&rsqb;$

(5.5)判断个体p_i的稀疏度K_i与种群的平均稀疏度K_mean的大小，并根据大小关系进行变异操作；

(5.5.1)如果K_i≤K_mean，记录个体p_i中基因值为0的索引集合：

Λ＝{λ|p_i[λ]＝0,1≤λ≤N}，

找到整体基因活性L在索引集合Λ中最大值的位置λ_max：

$\begin{array}{cc}L\&lsqb;{\mathrm{\&lambda;}}_{max}\&rsqb;\&GreaterEqual;L\&lsqb;\mathrm{\&lambda;}\&rsqb;& \&ForAll;{\mathrm{\&lambda;}}_{max},\mathrm{\&lambda;}\&Element;\mathrm{\&Lambda;},\end{array}$

置个体p_i的第λ_max个基因为1，即p_i[λ_max]＝1；

(5.5.2)如果K_i＞K_mean，记录个体p_i中基因值为1的索引集合：

Λ＝{λ|p_i[λ]＝1,1≤λ≤N}，

找到整体基因活性L在索引集合Λ中最小值的位置λ_min：

$\begin{array}{cc}L\&lsqb;{\mathrm{\&lambda;}}_{\min }\&rsqb;\&le;L\&lsqb;\mathrm{\&lambda;}\&rsqb;& \&ForAll;{\mathrm{\&lambda;}}_{min},\mathrm{\&lambda;}\&Element;\mathrm{\&Lambda;},\end{array}$

置个体p_i的第λ_min个基因为0，即p_i[λ_min]＝0；

(5.6)更新个体p_i的适应度值，基因活性，置i＝i+1，若i≤S，返回(5.2)；否则，输出变异后的子代种群