CN103337087B - 一种基于伪逆自适应算法的压缩感知图像重构方法 - Google Patents

Abstract

本发明涉及一种基于伪逆自适应算法的压缩感知图像重构方法，技术特征在于：首先选择一个高斯随机矩阵对图像进行投影以得到一个测量信号矩阵，然后构造图像稀疏基，并针对测量信号矩阵的每列采用本发明提出的伪逆自适应算法(Pseudo-inverse？Adaptive？Matching？Pursuit，PIAMP)进行重构以得到图像每列的稀疏系数，最后将图像每列稀疏系数进行排列并进行余弦变换得到重构图像。本发明方法在重构稀疏系数过程中具体考虑了字典相关度和迭代阶段划分等问题，故使稀疏系数的选择更精确，图像重构的精度更高。

Description

一种基于伪逆自适应算法的压缩感知图像重构方法

技术领域

本发明涉及一种压缩感知图像重构的新方法，即基于伪逆自适应算法(Pseudo-inverseAdaptiveMatchingPursuit，PIAMP)的压缩感知图像重构。其可以广泛的应用于各种以压缩感知理论作为支撑的图像处理系统中。

背景技术

在传统的图像编码传输过程中，首先要根据奈奎斯特采样定理（采样频率的设定要大于等于模拟信号中最高频率的两倍）对图像进行采样，之后再在一组稀疏基下对采样后的数据进行变换，得到图像的稀疏表达系数，然后再采用阈值处理的方法舍弃零或接近零的小数据，只对图像变换得到的大系数进行传输和处理。由上可以看出，其中先获取大量数据而后又舍弃的处理过程既浪费时间又浪费了硬件存储等物理资源。然而，2006年Candes等人提出的压缩感知理论（CompressiveSensing,CS）却解决了这一难题，CS理论的核心思想就是于在图像的前期编码中将采样和压缩过程合二为一，使其同时进行，如此便打破了传统的奈奎斯特采样定律。目前，基于CS理论的图像重构研究主要集中在三个方面，分别是测量矩阵的构造、图像的稀疏表示和重构算法的设计，其中重构算法的设计是图像重构研究的核心。

至今为止，很多种重构算法已被提出，主要有最小l₁范数类算法、贪婪类算法、迭代阈值法及基于最小全变差模型的算法等，其中贪婪类算法是这四种算法中重构速度最快，应用最广泛的一类算法。在贪婪类算法中，最基本的有正交匹配追踪算法(OrthogonalMatchingPursuit，OMP)，正则化正交匹配追踪算法(RegularizedOrthogonalMatchingPursuit，ROMP)、压缩采样匹配追踪算法(Compressivesamplingmatchingpursuit,CoSaMP)等。虽然基于这些算法的重构方法都可以精确的重构原图像，但它们有一个共同的缺点，即都要已知信号的稀疏度，而这在实际应用中这往往是不可达到的。于是在2008年Do等人提出了针对实际信号重构的自适应匹配追踪算法(Sparsityadaptivematchingpursuitalgorithm,SAMP)，之后赵瑞珍等人又在此基础上提出了正则化自适应匹配跟踪算法（RegularizationAdaptiveMatchingPursuitAlgorithm,RAMP）。而本发明则是在此算法的基础上，从最优原子的选择方式和支撑集的更新过程两方面进行改进，提出了PIAMP算法。由于PIAMP算法在迭代过程进行了更细致的阶段划分，故基于PIAMP算法的压缩感知图像重构结果在理论上可以达到更精确。

发明内容

要解决的技术问题

为提高压缩感知框架下图像重构的精度，减少重构图像的误差，本发明提出了一种压缩感知图像重构的新方法，即基于PIAMP算法的压缩感知图像重构。

技术方案

一种基于伪逆自适应算法的压缩感知图像重构方法，其特征在于步骤如下：

步骤1、投影测量：对于一幅维数为N×K的原始图像A，选择一个维数为M×N，且M≤N的高斯随机矩阵Φ对图像A进行投影测量，得到一个测量信号矩阵Y＝Φ·A;所述高斯随机矩阵Φ的每一个元素都服从均值为0、方差为的高斯分布；并且Φ的维数满足M≤N，得到的测量信号矩阵Y的维数为M×K,低于原始图像A的维数；

步骤2、构造图像的稀疏基：以图像的行数构造一个维数为N×N的离散余弦基Ψ作为图像的稀疏基；

离散余弦基Ψ的构造公式为：

$\mathrm{\Ψ}=\sqrt{\frac{2}{N}}{\left[C\left(k\right)\cos \left(\frac{(2m+1)\mathrm{k\π}}{2N}\right)\right]}_{N\×N},k,m=\mathrm{0,1},...,N-1$

$C\left(k\right)=\left\{\begin{array}{cc}\frac{\sqrt{2}}{2}& k=0\\ 1& k=\mathrm{1,2},....,N-1\end{array}\right.$

其中：m和k表示从0到N-1、包括0和N-1的整数值，C(k)表示构造每一个矩阵元素时的系数值，cos(·）表示对括号内的变量求余弦，[·]_N×N表示一个N×N维的矩阵；

步骤3、利用PIAMP重构图像每一列的稀疏系数：

步骤3：提取步骤1中得到的测量信号矩阵Y的每一列作为一维信号单独进行重构处理；以y_i表示矩阵Y的第i列，1≤i≤K，θ_i表示重构出图像第i列的稀疏系数，寻找每一列稀疏系数的具体过程如下：

(a)初始化参数设置：首先利用高斯随机矩阵Φ和离散余弦基Ψ的乘积求得图像的重构字典D＝Φ·Ψ;然后设置重构θ_i过程中控制阶段转换的阈值为ε₁，每一次迭代得到的余量为r，且初始余量为y_i，设置θ_i初始的稀疏度估计为p，第一次筛选中从重构字典D中所选择列向量的个数S＝p，并令迭代过程中存放从重构字典D中所选列向量列号的集合L、J₀、J₁都为空集；

(b)伪逆处理：将步骤(a)得到的重构字典D进行伪逆处理得到的矩阵Ω，

Ω＝D⁺＝(D·D^T)^-l·D

其中D⁺表示对矩阵D求伪逆，符号D^T表示对矩阵D求转置，(·)^-1表示对括号内的矩阵求逆运算；

(c)第一次筛选：首先通过求内积的方式来获得Ω中每一个列向量w_j(1≤j≤N)和余量r的相关系数u_j＝|r^T·W_j|，以表示每一个列向量与余量的相关度，其中|·|表示求绝对值符号；然后将求得的u_j由大至小排序，从中找出前S个值的索引值，并将其存入集合J₀；

(d)第二次筛选：首先将集合J₀中索引值对应的相关系数进行正则化处理，然后将处理得到的相关系数所对应的索引值放入集合J₁，并令L＝L∪J₁；

(e)阶段转换：利用下述两个公式求得本次迭代中所得到的稀疏系数θ_new和余量r₁，

θ_new＝(D_L)⁺·y

r₁＝y-D_Lθ_new

其中D_L是由矩阵D所构成的矩阵，其维数与D相同，且集合L中索引值所对应矩阵D_L的列向量与D相同；

然后检查||θ_new-θ_pre||₂≤ε₁是否满足，若满足则停止迭代，令θ_i＝θ_new；否则继续判断||θ_new-θ_pre||₂≤5ε₁是否满足，若满足则说明迭代过程已趋于停止，此时令θ_pre＝θ_new，r＝r₁，S＝S+p/3，返回步骤(c)继续进行迭代；否则进一步判断||r₁||₂≤||r||₂是否满足，若满足则令θ_pre＝θ_new，r＝r₁，且返回步骤(c)，否则令θ_pre＝θ_new，r＝r₁，S＝S+p，再返回步骤(c)，直到搜寻出θ_i；其中θ_pre表示上次迭代过程得到的稀疏系数；

步骤4：稀疏系数变换：将步骤3中搜索出的θ_i按照i的顺序进行排列，以得到一个维数为N×K的图像稀疏矩阵Θ；然后利用离散余弦基Ψ进行变换求得重构图像 $\stackrel{\‾}{A}=\mathrm{\Ψ}\·\mathrm{\Θ}.$

所述正则化处理的过程是：首先按照|u_t|≤2|u_h|(t,h∈J₀)的规则将J₀中索引值对应的相关系数进行分组，其中u_t和u_h分别表示J₀中任意两个不同索引值所对应的相关系数；然后选出能量最大的一组相关系数作为正则化的结果进行输出，其中能量的定义为e表示由任意一组相关系数所组成的一个一维向量，e₁,e₂,…e_z表示向量e的z个分量。||·||₂表示向量的l₂范数，表示向量l₂范数的平方。

所述阶段转换的阈值ε₁范围是0～0.01。

所述初始稀疏度估计p的范围为0.01M～0.2M，M为矩阵Φ的行数。

有益效果

本发明提出的一种的基于伪逆自适应算法的压缩感知图像重构方法，首先选择一个高斯随机矩阵对图像整体进行投影以得到一个测量信号矩阵，然后构造稀疏基并提取测量信号矩阵的每一列作为一个初始的一维信号进行处理。之后在具体的重构迭代中，首先利用离散余弦基和随机高斯矩阵的乘积来构造重构字典，并将字典进行伪逆处理以降低其相关度。然后通过求余量和伪逆处理后字典中原子内积过程和正则化过程分别对原子进行筛选。而且针对图像稀疏度未知的问题，本发明采用双阈值方法来控制原子搜索过程，以重构出图像每一列的稀疏系数。最后将每列的稀疏系数进行组合，并采用离散余弦变换得到重构图像。

本发明中基于伪逆自适应算法图像重构的优点在于：第一：对图像进行重构时，由于图像每一列的稀疏度不同，如果设置统一的稀疏度值进行重构，势必会造成信号的过匹配或者匹配不足，从而降低重构的质量，而本发明提出的算法保留了自适应匹配追踪算法的优势，采用一个预先设定的步长来估计不同一维信号的稀疏度，以提升重构的质量。第二：在图像重构的过程中加入了伪逆处理过程以降低重构字典的相关度，相较于直接采用重构字典进行原子的第一次筛选而言，更易搜索出较优的原子以备第二次筛选。第三：在原子的搜索过程中，双阈值的设置可以使搜索过程更平稳，搜索出的原子趋于最优，并且重构时所用的初始步长越小，图像可以精确重构的概率越高。

附图说明

图1：本发明方法重构图像的流程图

图2：PIAMP算法的流程图

图3：采样率为0.6时基于不同重构算法的Lena图像重构结果

(a)源图像；(b)基于OMP算法的重构结果；(c)基于RAMP算法的重构结果；(d)基于PIAMP算法的重构结果；

具体实施方式

现结合实施例、附图对本发明作进一步描述:

用于实施的硬件环境是：Pentium-43G计算机、2GB内存、128M显卡，运行的软件环境是：Mat1ab7.0和windowsXP。我们用Matlab程序设计语言实现了本发明提出的新算法。图像数据采用512×512的Lena、Pepper图像和Barbara图像。

本发明具体实施如下：

步骤1投影测量：对于一幅维数为N×K的原始图像A，选择一个维数为M×N(M≤N)的高斯随机矩阵Φ对图像A进行投影测量，以得到一个测量信号矩阵Y。在此所选矩阵Φ的每一个元素都服从均值为0，方差为的高斯分布。投影测量的具体过程如下所示：

Y＝Φ·A

其中由于Φ的维数满足M≤N，因此得到的测量信号矩阵Y的维数为M×K,低于原始图像A的维数。

步骤2构造图像的稀疏基：根据图像的行数构造一个维数为N×N的离散余弦基Ψ作为图像的稀疏基。离散余弦基Ψ的构造公式如下：

$\mathrm{\Ψ}=\sqrt{\frac{2}{N}}{\left[C\left(k\right)\cos \left(\frac{(2m+1)\mathrm{k\π}}{2N}\right)\right]}_{N\×N},k,m=\mathrm{0,1},...,N-1$

$C\left(k\right)=\left\{\begin{array}{cc}\frac{\sqrt{2}}{2}& k=0\\ 1& k=\mathrm{1,2},....,N-1\end{array}\right.$

其中m和k表示从0到N-1（包括0和N-1）的整数值，C(k)表示构造每一个矩阵元素时的系数值，cos(·)表示对括号内的变量求余弦，[·]_NXN表示一个N×N维的矩阵。

步骤3利用PIAMP重构图像每一列的稀疏系数：提取步骤1中得到的测量信号矩阵Y的每一列作为一维信号单独进行重构处理，假设y_i表示矩阵Y的第i列（1≤i≤K），θ_i表示重构出图像第i列的稀疏系数。则寻找每一列稀疏系数的具体过程如下：

(a)初始化参数设置：首先利用高斯随机矩阵Φ和离散余弦基Ψ的乘积求得图像的重构字典D，即D＝Φ·Ψ。然后设置重构θ_i过程中控制阶段转换的阈值为ε₁＝0.01，每一次迭代得到的余量为r，且初始余量为y_i。另外，设置θ_i初始的稀疏度估计为p＝0.05M，第一次筛选中从重构字典D中所选列向量的个数S＝p，并令迭代过程中存放从重构字典D中所选列向量列号（索引值）的集合L、J₀、J₁都为空集。

(b)伪逆处理：将步骤(a)得到的重构字典D进行伪逆处理以得到一个相关度较低的矩阵Ω，伪逆处理过程如下所示：

Ω＝D⁺＝(D·D^T)^-l·D

其中D⁺表示对矩阵D求伪逆，符号D^T表示对矩阵D求转置，(·)^-1表示对括号内的矩阵求逆运算。

(c)第一次筛选：首先通过求内积的方式来获得Ω中每一个列向量w_j(1≤j≤N)和余量r的相关系数u_j，以表示每一个列向量与余量的相关度，即u_j＝|r^T·w_j｜其中|·｜表示求绝对值符号。然后将求得的u_j由大至小排序，从中找出前S个值的索引值，并将其存入集合J₀。

(d)第二次筛选：首先将集合J₀中索引值对应的相关系数进行正则化处理，然后将处理得到的相关系数所对应的索引值放入集合J₁，并令L＝L∪J₁，其中正则化处理的过程是：首先按照|u_t|≤2|u_h|(t,h∈J₀)的规则将J₀中索引值对应的相关系数进行分组，其中u_t和u_h分别表示J₀中任意两个不同索引值所对应的相关系数。然后选出能量最大的一组相关系数作为正则化的结果进行输出。其中能量的定义为e表示由任意一组相关系数所组成的一个一维向量，e₁,e₂,...，e_z表示向量e的z个分量。||·||₂表示向量的l₂范数，表示向量l₂范数的平方。

(e)阶段转换：利用如下两个式子求得本次迭代中所得到的稀疏系数θ_new和余量r₁。

θ_new=(D_L)⁺·y

r₁＝y-D_Lθ_new

其中D_L是由矩阵D所构成的矩阵，其维数与D相同，且集合L中索引值所对应矩阵D_L的列向量与D相同。然后检查||θ_new-θ_pre||₂≤ε₁是否满足，θ_pre表示上次迭代过程得到的稀疏系数，若满足则停止迭代，令θ_i＝θ_new。否则继续判断||θ_new-θ_pre||₂≤5ε₁是否满足，若满足则说明迭代过程已趋于停止，此时令θ_pre＝θ_new，r＝r₁，S＝S+p/3，返回步骤(c)继续进行迭代。否则进一步判断||r₁||₂≤||r||₂是否满足，若满足则令θ_pre＝θ_new，r＝r₁，且返回步骤(c)，否则令θ_pre＝θ_new，r＝r₁，S＝S+p，再返回步骤(c)，直到搜寻出θ_i。

步骤4稀疏系数变换：将步骤3中搜索出的θ_i按照i的顺序进行排列，以得到一个维数为N×K的图像稀疏矩阵Θ。然后利用离散余弦基Ψ进行变换求得重构图像。稀疏系数变换的过程如下式所示，其中表示重构出的图像。

$\stackrel{\‾}{A}=\mathrm{\Ψ}\·\mathrm{\Θ}$

基于此，上述四步就是采用伪逆自适应算法重构图像的详细过程。

图2是采样率为0.6时基于不同重构算法的Lena图像重构结果，将本发明所得的重构结果与其他重构方法所得的结果进行对比，可以直观得看出，在测量矩阵和稀疏基选择相同的情况下，基于RAMP算法和PIAMP算法重构出的图像质量明显优于基于OMP算法重构的图像，尤其在Lena图像的帽子边缘及纹理部分，这种差距尤为明显，采用OMP算法重构得出的图像条纹现象明显。另外进一步比较RAMP算法和PIAMP算法重构得出的图像，视觉上差距较小，但可以根据其各自的客观指标值辨识出差距。

选用图像的峰值信噪比(PeakSignaltoNoiseRatio,PSNR)和相对误差(error)为客观评价指标，定量的评价不同图像采用不同算法重构效果。客观评价结果分别如表1、表2所示。其中峰值信噪比和相对误差的定义式如下所示。

PSNR＝10·log(255²/MSE)

$\mathrm{MSE}={\left|\right|I_{\mathrm{RE}}\left|\right|}_F/\sqrt{m\×n}$

$\mathrm{error}={\left|\right|I_{\mathrm{RE}}-I\left|\right|}_F^2/{\left|\right|I\left|\right|}_F^2$

其中I始图像，I_RE构出的图像，m和n是图像的尺寸大小，MSE构图像的均方误差。将OMP、RAMP和PIAMP算法应用到多幅图像重构时，由在不同采样率下得到的客观指标值可以看出，RAMP算法与PIAMP算法重构图像的PSNR比OMP算法重构图像的PSNR高约1.9～3dB，相对误差约降低了0.003～0.02，而进一步比较RAMP算法和PIAMP算法的重构结果，PIAMP算法重构图像的PSNR最高可比RAMP高0.6dB，相对误差最大可降低0.01。故综上所述，无论从主观视觉还是从客观指标上分析，本发明所提的PIAMP重构算法可以进一步提高图像的重构精度。

表1不同图像重构的峰值信噪比

表2不同图像重构的相对误差

Claims (4)

1.一种基于伪逆自适应算法的压缩感知图像重构方法，其特征在于步骤如下：

步骤1、投影测量：对于一幅维数为N×K的原始图像A，选择一个维数为M×N，且M≤N的高斯随机矩阵Φ对图像A进行投影测量，得到一个测量信号矩阵Y＝Φ·A；所述高斯随机矩阵Φ的每一个元素都服从均值为0、方差为的高斯分布；并且Φ的维数满足M≤N，得到的测量信号矩阵Y的维数为M×K,低于原始图像A的维数；

步骤2、构造图像的稀疏基：以图像的行数构造一个维数为N×N的离散余弦基Ψ作为图像的稀疏基；

离散余弦基Ψ的构造公式为：

$\mathrm{\Ψ}=\sqrt{\frac{2}{N}}{\[C\left(k\right)cos\left(\frac{(2m+1)k\mathrm{\π}}{2N}\right)\]}_{N\×N},k,m=0,1,...,N-1$

$C\left(k\right)=\left\{\begin{array}{cc}\frac{\sqrt{2}}{2}& k=0\\ 1& k=1,2,\mathrm{....},N-1\end{array}\right.$

其中：m和k表示从0到N-1、包括0和N-1的整数值，C(k)表示构造每一个矩阵元素时的系数值，cos(·)表示对括号内的变量求余弦，[·]_N×N表示一个N×N维的矩阵；

步骤3、利用伪逆自适应算法PIAMP重构图像每一列的稀疏系数：

提取步骤1中得到的测量信号矩阵Y的每一列作为一维信号单独进行重构处理；以y_i表示矩阵Y的第i列，1≤i≤K，θ_i表示重构出图像第i列的稀疏系数，寻找每一列稀疏系数的具体过程如下：

(a)初始化参数设置：首先利用高斯随机矩阵Φ和离散余弦基Ψ的乘积求得图像的重构字典D＝Φ·Ψ；然后设置重构θ_i过程中控制阶段转换的阈值为ε₁，每一次迭代得到的余量为r，且初始余量r₀为y_i，设置θ_i初始的稀疏度估计为p，第一次筛选中从重构字典D中所选择列向量的个数S＝p，并令迭代过程中存放从重构字典D中所选列向量列号的集合L、J₀、J₁都为空集；

(b)伪逆处理：将步骤(a)得到的重构字典D进行伪逆处理得到矩阵Ω，

Ω＝D⁺＝(D·D^T)^-1·D

其中D⁺表示对矩阵D求伪逆，符号D^T表示对矩阵D求转置，(·)^-1表示对括号内的矩阵求逆运算；

(c)第一次筛选：首先通过求内积的方式来获得Ω中每一个列向量w_j，1≤j≤N和余量r的相关系数u_j＝|r^T·w_j|，以表示每一个列向量与余量的相关度，其中|·|表示求绝对值符号；然后将求得的u_j由大至小排序，从中找出前S个值的索引值，并将其存入集合J₀；

(d)第二次筛选：首先将集合J₀中索引值对应的相关系数进行正则化处理，然后将处理得到的相关系数所对应的索引值放入集合J₁，并令L＝L∪J₁；

(e)阶段转换：利用下述两个公式求得本次迭代中所得到的稀疏系数θ_new和余量r₁，

θ_new＝(D_L)⁺·y_i

r₁＝y_i-D_Lθ_new

其中D_L是由矩阵D所构成的矩阵，其维数与D相同，且集合L中索引值所对应矩阵D_L的列向量与D相同；

然后检查||θ_new-θ_pre||₂≤ε₁是否满足，若满足则停止迭代，令θ_i＝θ_new；否则继续判断||θ_new-θ_pre||₂≤5ε₁是否满足，若满足则说明迭代过程已趋于停止，此时令θ_pre＝θ_new，r＝r₁，S＝S+p/3，返回步骤(c)继续进行迭代；否则进一步判断||r₁₂||≤||r₂||是否满足，若满足则令θ_pre＝θ_new，r＝r₁，且返回步骤(c)，否则令θ_pre＝θ_new，r＝r₁，S＝S+p，再返回步骤(c)，直到搜寻出θ_i；其中θ_pre表示上次迭代过程得到的稀疏系数，||·||₂表示向量的l₂范数；

步骤4：稀疏系数变换：将步骤3中搜索出的θ_i按照i的顺序进行排列，以得到一个维数为N×K的图像稀疏矩阵Θ；然后利用离散余弦基Ψ进行变换求得重构图像 $\stackrel{\‾}{A}=\mathrm{\Ψ}\·\mathrm{\Θ}.$

2.根据权利要求1所述基于伪逆自适应算法的压缩感知图像重构方法，其特征在于：所述正则化处理的过程是：首先按照|u_t|≤2|u_h|，t、h∈J₀的规则将J₀中索引值对应的相关系数进行分组，其中u_t和u_h分别表示J₀中任意两个不同索引值所对应的相关系数；然后选出能量最大的一组相关系数作为正则化的结果进行输出，其中能量的定义为e表示由任意一组相关系数所组成的一个一维向量，e₁,e₂,…,e_z表示向量e的z个分量，表示向量l₂范数的平方。

3.根据权利要求1所述的基于伪逆自适应算法的压缩感知图像重构方法，其特征在于：所述阶段转换的阈值ε₁范围是0～0.01。

4.根据权利要求1所述的基于伪逆自适应算法的压缩感知图像重构方法，其特征在于：所述初始稀疏度估计p的范围为0.01M～0.2M，M为矩阵Φ的行数。