CN102882530B - 一种压缩感知信号重构方法 - Google Patents

Abstract

本发明公开了一种压缩感知信号重构方法，属于信号处理技术领域。本发明针对 - 正则项化问题求解比较困难，将压缩感知信号的稀疏域的 - 正则化问题通过变量分裂技术转化为与之等价的，比 - 正则项化问题更能体现信号稀疏特性的约束 - 正则项化问题，使得重构信号的精度更高。本发明进一步利用快速交替方向乘子法对约束 - 正则项化问题进行求解，对交替方向乘子法算法的变量进行了二次更新，并更新了乘子，加快了优化求解的收敛速度。相比现有技术，本发明方法具有更高的重构精度及更快的收敛速度。

Description

一种压缩感知信号重构方法

技术领域

本发明涉及一种压缩感知(Compressed Sensing，简称CS)信号重构方法，属于信号处理技术领域。

背景技术

针对具有稀疏特性信号的压缩感知(Compressed Sensing，CS)理论由Donoho等人于2004年提出。它在保证信号不受损失的情况下，用远低于Nyquist采样定理要求的速率采集信号，同时又不损失信息，能够完全恢复信号，是信号处理领域的一大革命性成果。在CS理论框架下，采样速率不再取决于信号的带宽，而取决于信息在信号中的结构和内容。

CS理论中的核心问题是压缩信号的重构问题，如何设计复杂度低、收敛速度快、鲁棒性强的重构算法一直是CS理论重构算法研究的目标。另一方面，实际的信号处理环境中，采集（或者压缩感知得到）的信号，经常是含有噪声的。在使用优化方法重构含噪信号时，不同之处在于重构过程所使用的优化目标函数的形式不同，参数的设置不同。应用不同的优化目标函数，导致信号的重建效果也不尽相同。CS中重构含噪信号的原始模型为l₀-正则化问题，但该问题是一个非凸问题，求解十分困难，迭代硬阈值(Iterative Hard Thresholding,IHT)算法虽然可以通过迭代的方法求解该问题，但IHT算法求解的是该问题局部极小值，重构精度低且算法收敛速度慢。所以，通常情况下，将其转化为l₁-正则化问题进行求解。虽然l₁-正则化问题是一个严格凸问题，可以求得其全局唯一解，但是l₁-正则化问题只有在一定条件下才与l₀-正则化问题等价，所求的解也只有在一定条件下才是原问题的解，此时才能精确的重构原信号。

目前关于含噪信号重构的领域已经有一些研究成果。交替方向乘子法(AlternatingDirection Method of Multipliers,ADMM)以其较好的性能和严格的理论保证，成为众多学者研究的热点。例如用ADMM解决l₁-正则化问题，提出了SALSA(Split augmentedLagrangian Shrinkage Algorithm)和C-SALSA(Constrained-SALSA)算法，并把其应用于图像的去模糊、解卷积等问题，取得了较好的效果，但当观测矩阵为随机矩阵时这两种算法的计算复杂度较高；又如用最速下降法更新迭代变量代替计算量大的求逆运算，但其算法的收敛性没法保证；也有用ADMM来解决TV问题,并把其应用于图像的去模糊、去噪问题；再如用ADMM来解决TVl₁-l₂问题，使用部分Fourier矩阵作为其观测矩阵，但当观测矩阵为随机矩阵时收敛速度较慢；再如用IADMM(Inexact ADMM)来解决TV问题，加快了算法的收敛速度，但其过程中求得的解只是原问题的近似解，导致重构的精度不高。大部分都是将ADMM应用于图像的去噪、去模糊、解卷积等问题，很少有将ADMM应用于投影矩阵为随机矩阵的信号重构问题；另外，目前尚无用ADMM来解决约束l₀-正则化问题。

发明内容

本发明所要解决的技术问题在于克服现有技术的不足，提供一种压缩感知信号重构方法，使得重构的信号的精度更高、算法的收敛速度更快。

本发明具体采用以下技术方案解决上述技术问题：

一种压缩感知信号重构方法，首先通过求解以下优化问题得到原始信号的稀疏系数Θ：

$\underset{\mathrm{\Θ},\mathrm{\Ξ}}{\min }\frac{1}{2}\left|\right|\mathrm{A\Θ}-y|{|}_2^2+\mathrm{\τ}{\left|\right|\mathrm{\Ξ}\left|\right|}_0$

s.t.Θ=Ξ

其中，A为压缩感知采样的观测矩阵，Ξ为辅助变量，y为对原始信号进行压缩感知采样获得的观测向量，τ∈[0,∞)为正则化参数；

然后根据得到的稀疏系数Θ重构出原始信号。

针对l₀-正则项化问题求解比较困难，本发明将压缩感知信号的稀疏域的l₀-正则化问题通过变量分裂(Variable Splitting,VS)技术转化为与其等价的约束l₀-正则项化问题。由于l₀-正则化问题比l₁-正则化问题更能体现信号的稀疏特性，因此使得重构的信号的精度更高。对于转化后的约束l₀-正则项问题，可采用交替方向乘子法、对偶交替方向乘子法、乘子法、非精确交替方向乘子法、交替线性化法等现有的方法进行求解。本发明为了加快优化求解算法的迭代速度，提高信号重构的效率和实时性，进一步对交替方向乘子法进行改进，具体包括以下步骤：

步骤1、初始化：设置初始惩罚参数μ>0，初始辅助变量Ξ₀，初始乘子d₀，初始变量t₀=1，初始迭代次数k=0；

步骤2、更新稀疏系数变量Θ：

${\stackrel{\‾}{\mathrm{\Θ}}}_{k+1}=\frac{1}{\mathrm{\μ}}[I-A^T{(\mathrm{\μI}+{\mathrm{AA}}^T)}^{-1}A](A^{T}y+\mathrm{\μ}({\mathrm{\Ξ}}_k+d_k)),$ 其中I为单位矩阵；

步骤3、更新变量t： $t_{k+1}=\frac{1+\sqrt{1+{4t}_k^2}}{2};$

步骤4、二次更新稀疏系数变量Θ： ${\mathrm{\Θ}}_{k+1}={\stackrel{\‾}{\mathrm{\Θ}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\‾}{\mathrm{\Θ}}}_{k+1}-{\stackrel{\‾}{\mathrm{\Θ}}}_k);$

步骤5、更新辅助变量Ξ：其中hard(,τ)为阈值为τ的硬阈值函数；

步骤6、二次更新辅助变量Ξ： ${\mathrm{\Ξ}}_{k+1}={\stackrel{\‾}{\mathrm{\Ξ}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\‾}{\mathrm{\Ξ}}}_{k+1}-{\stackrel{\‾}{\mathrm{\Ξ}}}_k);$

步骤7、更新乘子d：d_k+1=d_k-(Θ_k+1-Ξ_k+1)；

步骤8、判断是否满足终止条件，满足则停止；否则，令k=k+1，转步骤4。

相比现有技术，本发明具有以下有益效果：

(1)通过变量分裂技术将压缩感知信号重构的l₀-正则化问题转化为约束l₀-正则化问题，该正则化问题比l₁-正则化问题更体现信号的稀疏特性，使得重构的信号的精度更高；

(2)在对约束l₀-正则化问题求解时，用加速策略对ADMM的变量进行了二次更新，加快了算法的收敛速度。

附图说明

图1为本发明的压缩感知信号重构系统框图；

图2为采用各种重构算法得到的含噪图像压缩感知重构结果。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

为了便于公众理解本发明的技术方案，下面先对压缩感知的基本知识做一简要介绍。

现代信号处理的一个关键基础是Shannon采样理论：一个信号可以无失真重建所要求的离散样本数由其带宽决定。但是Shannon采样定理是一个信号重建的充分非必要条件。而压缩感知作为一个新的采样理论，它可以在远小于Nyquist采样率的条件下获取信号的离散样本，保证信号的无失真重建。压缩感知理论的核心思想主要包括两点。第一个是信号的稀疏结构。传统的Shannon信号表示方法只开发利用了很少的被采样信号的先验信息，即信号的带宽。但是，现实生活中很多广受关注的信号本身具有一些结构特点。相对于带宽信息的自由度，这些结构特点是由信号的更小的一部分自由度所决定。换句话说，在很少的信息损失情况下，这种信号可以用很少的数字编码表示。所以，在这种意义上，这种信号是稀疏信号（或者可压缩信号）。另外一点是不相关特性。稀疏信号的有用信息的获取可以通过一个非自适应的采样方法将信号压缩成较小的样本数据来完成。理论证明压缩感知的采样方法只是一个简单的将信号与一组确定的波形进行相关的操作。这些波形要求是与信号所在的稀疏空间不相关的。压缩感知可广泛应用于信息论、图像处理、地球科学、光学/微波成像、模式识别、无线通信、大气、地质等领域。

具体地讲，对于任意信号x∈R^N，可以用N×1维正交基向量Ψ的线性组合来表示

其中，为投影系数，Θ=Ψ^Tx为投影系数矢量。若Θ是K-稀疏的，可以用一个与Ψ不相关的矩阵A∈R^M×N(其中M＜＜N)对Θ进行线性测量，得到观测向量y∈R^M为

y=AΘ+e                (2)

其中，e为噪声(当不含有噪声时，令e=0即可)，A∈R^M×N为与Ψ不相关的观测矩阵。重构的方法是通过求解如下最优化问题

$\underset{\mathrm{\Θ}}{\min }\frac{1}{2}{\left|\right|\mathrm{A\Θ}-y\left|\right|}_2^2+\mathrm{\τ}{\left|\right|\mathrm{\Θ}\left|\right|}_0---\left(3\right)$

其中τ∈[0,∞)为正则化参数(当τ=0即为不含噪声时的重构模型)。通过求解该优化问题即可得到稀疏系数Θ，进而可得到原始信号x。

正如背景技术部分所述，式（3）的l₀-正则化问题是一个非凸问题，求解十分困难，因此通常将其转化为l₁-正则化问题进行求解。虽然l₁-正则化问题是一个严格凸问题，可以求得其全局唯一解，但其具有较大局限性。

本发明的思路是利用变量分裂技术，引入辅助变量Ξ，将稀疏域含噪信号的重构问题转化为如下的约束l₀-正则化问题

$\underset{\mathrm{\Θ},\mathrm{\Ξ}}{\min }\frac{1}{2}{\left|\right|\mathrm{A\Θ}-y\left|\right|}_2^2+\mathrm{\τ}{\left|\right|\mathrm{\Ξ}\left|\right|}_0---\left(3\right)$

s.t.Θ=Ξ

此时即可用具有较好的性能和严格的理论保证的交替方向乘子法（ADMM）对式(4)的约束l₀-正则化问题进行求解，得到信号的更加精确的稀疏系数，进而得到原始信号。

约束l₀-正则化问题(4)的乘子函数为

$\mathrm{\Γ}(\mathrm{\Θ},\mathrm{\Ξ},\mathrm{\λ})=\frac{1}{2}{\left|\right|\mathrm{A\Θ}-y\left|\right|}_2^2+\mathrm{\τ}{\left|\right|\mathrm{\Ξ}\left|\right|}_0-{\mathrm{\λ}}^T(\mathrm{\Θ}-\mathrm{\Ξ})+\frac{\mathrm{\μ}}{2}{\left|\right|\mathrm{\Θ}-\mathrm{\Ξ}\left|\right|}_2^2---\left(5\right)$

其中，λ∈R^M为乘子，μ∈[0,∞)为惩罚参数。

对于变量Θ得到：

${\stackrel{\‾}{\mathrm{\Θ}}}_{k+1}=\underset{\mathrm{\Θ}}{\min }\frac{1}{2}{\left|\right|\mathrm{A\Θ}-y\left|\right|}_2^2-{\mathrm{\λ}}_k^T(\mathrm{\Θ}-{\mathrm{\Ξ}}_k)+\frac{\mathrm{\μ}}{2}{\left|\right|\mathrm{\Θ}-{\mathrm{\Ξ}}_k\left|\right|}_2^2=\underset{\mathrm{\Θ}}{\min }\frac{1}{2}{\left|\right|\mathrm{A\Θ}-y\left|\right|}_2^2+\frac{\mathrm{\μ}}{2}{\left|\right|\mathrm{\Θ}-({\mathrm{\Ξ}}_k+d_k)\left|\right|}_2^2$

$=\underset{\mathrm{\Θ}}{\min }{\left|\right|\mathrm{A\Θ}-y\left|\right|}_2^2+\mathrm{\μ}{\left|\right|\mathrm{\Θ-}({\mathrm{\Ξ}}_k+d_k)\left|\right|}_2^2={(A^{T}A+\mathrm{\μI})}^{-1}(A^{T}y+\mathrm{\μ}({\mathrm{\Ξ}}_k+d_k))$

(6)

其中，也称为乘子，I为单位矩阵。当矩阵A∈R^M×N为随机矩阵时，(A^TA+μI)^-1的计算量为O(N³)，为了减少其计算量，根据Sherman-Morrison-Woodburg定理得到

${(A^{T}A+\mathrm{\μI})}^{-1}=\frac{1}{\mathrm{\μ}}[I-A^T{(\mathrm{\μI}+{\mathrm{AA}}^T)}^{-1}A]---\left(7\right)$

问题(6)是一个严格凸问题，其解唯一。从而得到

${\stackrel{\‾}{\mathrm{\Θ}}}_{k+1}=\frac{1}{\mathrm{\μ}}[I-A^T{(\mathrm{\μI}+{\mathrm{AA}}^T)}^{-1}A](A^{T}y+\mathrm{\μ}({\mathrm{\Ξ}}_k+d_k))---\left(8\right)$

对于变量Ξ得到：

${\stackrel{\‾}{\mathrm{\Ξ}}}_{k+1}=\underset{\mathrm{\Ξ}}{\min }\mathrm{\τ}{\left|\right|\mathrm{\Ξ}\left|\right|}_0-{\mathrm{\λ}}_k^T({\mathrm{\Θ}}_{k+1}-\mathrm{\Ξ})+\frac{\mathrm{\μ}}{2}{\left|\right|{\mathrm{\Θ}}_{k+1}-\mathrm{\Ξ}\left|\right|}_2^2=\underset{\mathrm{\Ξ}}{\min }\mathrm{\τ}{\left|\right|\mathrm{\Ξ}\left|\right|}_0+\frac{\mathrm{\μ}}{2}{\left|\right|\mathrm{\Ξ}-({\mathrm{\Θ}}_{k+1}-d_k)\left|\right|}_2^2---\left(9\right)$

$=\mathrm{hard}({\mathrm{\Θ}}_{k+1}-d_k,\mathrm{\τ}/\mathrm{\μ})$

其中hard(,τ)为阈值为τ的硬阈值函数。为了加速ADMM的收敛速度，本发明对Θ和Ξ进行了二次更新，其二次更新的形式分别为

${\mathrm{\Θ}}_{k+1}={\stackrel{\‾}{\mathrm{\Θ}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\‾}{\mathrm{\Θ}}}_{k+1}-{\stackrel{\‾}{\mathrm{\Θ}}}_k)---\left(10\right)$

${\mathrm{\Ξ}}_{k+1}={\stackrel{\‾}{\mathrm{\Ξ}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\‾}{\mathrm{\Ξ}}}_{k+1}-{\stackrel{\‾}{\mathrm{\Ξ}}}_k)---\left(11\right)$

其中t₀=1，为算法迭代更新的步长。且对乘子进行了更新

d_k+1=d_k-(Θ_k+1-Ξ_k+1)                  (12)

本发明称上述的改进交替方向乘子法算法为快速ADMM，即FADMM。之所以FADMM会加快算法的收敛速度，是因为在迭代过程中不仅仅用到当前点的信息(即算法ADMM)，而且用到的信息，使得每步的迭代结果更准确，在终止准则不变时，使得算法的迭代次数减少；而且FADMM的主要计算量仍在式(6)和(9)，对Θ和Ξ进行的二次更新并无实质性的增加计算量，即FADMM每次迭代的运算量并无实质性的增加。从而导致FADMM的运行时间因其收敛速度快而减少。之所以取t₀=1是为了保证对于任意的k≥0，均有从而保证FADMM的超线性收敛速度。

综上所述，利用FADMM求解约束l₀-正则化含噪信号重构问题的算法步骤为

1)初始化：设置初始惩罚参数μ>0，初始辅助变量Ξ₀，初始乘子d₀，初始变量t₀=1，迭代次数k=0；

2)更新稀疏系数变量Θ： ${\stackrel{\‾}{\mathrm{\Θ}}}_{k+1}=\frac{1}{\mathrm{\μ}}[I-A^T{(\mathrm{\μI}+{\mathrm{AA}}^T)}^{-1}A](A^{T}y+\mathrm{\μ}({\mathrm{\Ξ}}_k+d_k));$

3)更新变量t： $t_{k+1}=\frac{1+\sqrt{1+{4t}_k^2}}{2};$

4)二次更新稀疏系数变量Θ： ${\mathrm{\Θ}}_{k+1}={\stackrel{\‾}{\mathrm{\Θ}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\‾}{\mathrm{\Θ}}}_{k+1}-{\stackrel{\‾}{\mathrm{\Θ}}}_k);$

5)更新辅助变量Ξ： ${\stackrel{\‾}{\mathrm{\Ξ}}}_{k+1}=\mathrm{hard}({\mathrm{\Θ}}_{k+1}-d_k,\mathrm{\τ}/\mathrm{\μ});$

6)二次更新辅助变量Ξ： ${\mathrm{\Ξ}}_{k+1}={\stackrel{\‾}{\mathrm{\Ξ}}}_{k+1}+\left(\frac{t_k-1}{t_{k+1}}\right)({\stackrel{\‾}{\mathrm{\Ξ}}}_{k+1}-{\stackrel{\‾}{\mathrm{\Ξ}}}_k);$

7)更新乘子d：d_k+1=d_k-(Θ_k+1-Ξ_k+1)；

8)判断是否满足终止条件，满足则停止；否则，令k=k+1，转步骤4)。

定理1：令设为FADMM生成的序列，则对任意的k≥0，有

$\frac{F\left(\begin{array}{c}{\mathrm{\Θ}}_k\\ {\mathrm{\Ξ}}_k\end{array}\right)-F\left(\begin{array}{c}{\mathrm{\Θ}}^{*}\\ {\mathrm{\Ξ}}^{*}\end{array}\right)}{\left|\right|{\left(\begin{array}{c}{\mathrm{\Θ}}_0\\ {\mathrm{\Ξ}}_0\end{array}\right)-\left(\begin{array}{c}{\mathrm{\Θ}}^{*}\\ {\mathrm{\Ξ}}^{*}\end{array}\right)\left|\right|}_2^2}=O\left(\frac{1}{{(k+1)}^2}\right)---\left(13\right)$

其中， $\left(\begin{array}{c}{\mathrm{\Θ}}_k\\ {\mathrm{\Ξ}}_k\end{array}\right)\→\left(\begin{array}{c}{\mathrm{\Θ}}^{*}\\ {\mathrm{\Ξ}}^{*}\end{array}\right).$

定理1说明FADMM具有超线性收敛性。

图1显示了一种采用本发明的方法的压缩感知信号重构系统，如图所示，该系统包括编码端和解码端，解码端对原始的含噪信号先进行小波变换，提取出原始含噪信号的小波系数，然后对小波系数进行压缩感知采样，得到观测量；解码端采用本发明的FADMM算法对编码端发送的观测量进行优化求解，得到重构信号的小波系数，然后进行反小波变换，即得到重构后的信号。

为了验证本发明方法的效果，以含噪图像的压缩感知重构为例，对本发明的重构方法与其他现有方法进行了比较。实验对象为Lena图像(大小为256×256)，正交基矩阵Ψ为小波基db3，观测矩阵A为服从高斯分布的随机矩阵，压缩比r=0.5，使用imnoise函数对实验对象加入均值为0，方差为σ²的高斯噪声。分别采用峰值信噪比PSNR(单位为dB)和运行时间t(单位为S)作为算法的重构性能和收敛速度的评价指标。重复进行10次实验，各种算法的峰值信噪比PSNR和运行时间t均为10次实验的平均结果。并将本发明的FADMM算法与SALSA、AIHT(Accelerated Iterative Hard Thresholding)算法、IHT算法、TwIST(Two-step Iterative Shrinkage/Thresholding)算法、GPSR(Gradient Projection forSparse Reconstruction)算法的性能和收敛速度进行了比较与分析。

对实验对象的原图像加入均值为0，方差为σ²=0.001的高斯噪声，图2显示了实验对象的原图像、含噪图像及各种算法重构的图像，可以看出，当压缩比r=0.5时，对原图像加入均值为0，方差为σ²=0.001的高斯噪声，各种算法重构的图像效果均较好。

表1显示了各种重构算法在不同噪声程度下，采用PSNR(单位为dB)作为衡量它们重构性能指标的结果。

表1不同噪声程度各种重构算法下的PSNR

σ2	FADMM	SALSA	AIHT	IHT	TwIST	GPSR
							0.001	32.055	32.042	29.879	29.736	31.969	32.038
0.005	28.039	27.771	25.286	25.038	27.953	27.685
							0.01	25.668	24.531	22.804	22.641	25.601	25.646
0.05	19.838	18.937	16.932	16.602	19.082	18.883
							0.1	17.608	15.763	14.647	14.354	17.491	17.537

随着噪声方差σ²的增大，各种算法重构图像的PSNR减小。此外在相同噪声程度下，FADMM算法重构图像的PSNR最大，IHT算法重构的图像的PSNR最小。而PSNR在一定程度上反映算法的重建质量，同等条件下PSNR越大，算法重建质量越高。所以在相同噪声程度下，FADMM算法重构图像的质量高于其他算法(即SALSA、AIHT算法、IHT算法、TwIST算法、GPSR算法)重构图像的质量。

表2显示了各种重构算法在不同噪声程度下，采用运行时间t(单位为S)作为衡量它们重构性能指标的结果。

表2不同噪声程度各种重构算法下的t

σ2	FADMM	SALSA	AIHT	IHT	TwIST	GPSR
							0.001	3.860	4.695	4.561	11.767	24.031	56.375
0.005	4.249	5.454	5.079	11.924	30.352	58.239
							0.01	4.203	7.218	4.960	11.999	38.281	59.996
0.05	4.405	11.852	5.454	12.010	42.399	60.925
							0.1	4.842	12.406	5.639	12.313	48.420	62.817

随着噪声方差σ²的增大，各种重构算法重构图像的运行时间t增大。此外相同噪声程度下，FADMM算法的运行时间最短，GPSR算法的运行时间最长，各种算法的运行时间快慢可以表示为t_FADMM<t_AIHT<t_SALSA<t_IHT<t_TwIST<t_GPSR。

综合考虑峰值信噪比PSNR和运行时间t，FADMM算法的性能最好。

Claims (1)

1.一种压缩感知信号重构方法，其特征在于，首先通过求解以下优化问题得到原始信号的稀疏系数                                               ：

其中，为压缩感知采样的观测矩阵，为辅助变量，为对原始信号进行压缩感知采样获得的观测向量，为正则化参数；

然后根据得到的稀疏系数重构出原始信号；

采用快速交替方向乘子法进行所述优化问题的求解，具体包括以下步骤：

步骤1、初始化：设置初始惩罚参数，初始辅助变量，初始乘子，初始变量，初始迭代次数；

步骤2、更新稀疏系数变量：

，其中为单位矩阵；

步骤3、更新变量：；

步骤4、二次更新稀疏系数变量：；

步骤5、更新辅助变量：，其中为阈值为的硬阈值函数；

步骤6、二次更新辅助变量：；

步骤7、更新乘子：；

步骤8、判断是否满足终止条件，满足则停止；否则，令，转步骤 4。