CN102163338B - 一种压缩感知系统中的高效重建方法 - Google Patents

Abstract

一种压缩感知系统中的高效重建方法。它涉及一种数据处理方法，它解决了现有重建方法中不能够精度速度同时提高的问题。首先整理测量值Y0为易于重建算法实现的形式，若一维重建，则不整理，若二维重建，则进行矢量化，得到Y；然后，令k＝1，u^k＝0，v^k＝0，得到u^k+1＝δ·shrink(v^k+1，μ)；v^k+1＝v^k+Θ^T(Y-Θu^k)；迭代步骤中出现的无贡献迭代，计算求取无贡献迭代的次数s，则假设，v^k变化s次恰好使得u^k+1有所改变，那么在这些迭代步骤中有如下迭代公式：u^k+s＝u^k+1，进行判定即||u^k+1-u^k||≤ε，再判断是否成立，来确定迭代是否收敛，迭代直至收敛；最后，若一维信号，则直接利用信号稀疏表达重建原始信号，若二维信号，则对稀疏系数u进行逆矢量化，并利用图像的稀疏表达重建原始图像。本发明应用于压缩感知系统中一维或二维信号重建。

Description

一种压缩感知系统中的高效重建方法

技术领域

本发明涉及一种数据处理方法，具体涉及一种压缩感知系统中的信号重建方法。

背景技术

近年来，在国际上出现了一种新的理论——压缩感知(compressive sensing)理论。该理论指出，只要信号是可压缩的或者说信号在某一个变换域是稀疏的，那么就可以用一个与稀疏矩阵不相关的测量矩阵对原始信号进行测量，将高维图像投影到一个低维空间中，然后通过求解一个优化问题就可以从这较少的投影值中完全重建出原始信号。压缩感知无需经过先采样再压缩的过程，这在很大的程度上克服了传统信号获取和处理的缺点。

压缩感知测量和重建过程描述如下：

现有的一维信号X(N×1维)和二维信号X2(N1×N2维)的压缩感知测量和重建过程描述如下：

根据稀疏矩阵Ψ(N×N维)，得到一维信号X的稀疏域表示：

X＝ΨS                        公式一

其中，S是一维信号X的稀疏域表示系数(N×1维)，具有K个非零系数和(N-K)个零；

对于一维信号X的压缩感知测量过程如下：

利用测量矩阵Φ(M×N维)，将一维信号X投影到M(M＜＜N)个测量值Y(M×1维)上；

Y＝ΦX＝ΦΨS＝ΘS                 公式二

其中Θ＝ΦΨ，Θ是个M×N维矩阵；

在测量过程后，先利用公式二从测量值Y中求取出一维信号X的稀疏域表示系数S，然后再利用公式一重建出原始一维信号X，在重建的过程中利用了一维信号X的稀疏性；

对于二维信号X2(N1×N2维)，利用二维信号的行稀疏矩阵Ψ2(N1×N1维)和二维信号的列稀疏矩阵Ψ3(N2×N2维)，得到二维信号X2的稀疏域表示：

X2＝Ψ2S2Ψ3                        公式三

其中，S2是二维信号X2的稀疏域表示系数(N1×N2维)；

利用二维信号的测量矩阵Φ2(M×N1×N2维)，获得M个二维信号的测量值Y2(M×1维)，上述过程用如下公式表示：

Y2＝Φ2*X2＝Φ2*(Ψ2S2Ψ3)                     公式四

上面的公式中主要是表示 $Y2{\left(i\right)}_{i=1:M}=\underset{k=1:N2}{\underset{j=1:N1}{\text{\Σ}}}\mathrm{\Φ}2(i,j,k)X2(j,k),$ 用*符号表示；

由于压缩感知采用非自适应性测量进行欠采样，将信号/图像的获取和压缩融为一体，使得对前台的硬件要求大大降低，获取的主要任务从前台转移到了后台(需要对进行非线性重建，强大的计算机支持)其中，二维信号(本发明中指图像)。因此，压缩感知的重建算法的好坏直接影响到压缩感知理论能否实用。压缩感知的重建实际上就是在满足一定的观测值的条件下，寻求最稀疏解的过程，是个l₀非凸优化的问题。然而，基于l₀范数极小化的方法是个需要寻找最优组合的NP-hard问题，数值计算上无法有效实现。因此一般采用l₁优化逼近l₀优化。重建主要可以分为三类方法：

(1)凸优化方法，包括基追踪BP、梯度追踪GPSR、Bregman迭代、最小角回归(LARs)等方法。此类方法重建精度高，需要的压缩测量个数少，但是计算复杂度相对较高。

(2)贪婪方法，主要包括：匹配追踪(MP)、正交匹配追踪(OMP)、的梯度追踪(GP)、分段匹配追踪(StOMP)、子空间追踪(SP)等。此类方法计算复杂度相对较低，运算速度快、但是总的来说与凸优化方法相比，需要更多的压缩测量，重建精度相对较低。

(3)非凸方法，主要包括FOCUSS算法、迭代重新加权算法、多层Bayesian压缩感知方法等。该类方法所需的压缩测量个数、计算复杂度、重构精度总体上介于上述两类方法之间。

总结上述的重建方法，要么重建速度快但是需要的测量数据多且精度不高，如贪婪算法中的OMP等；要么需要的测量数据少，精度高但是重建速度慢，如GPSR等。而在实际应用中，重建精度高是最重要的，因为其决定了压缩感知系统应用的效果。而重建速度快也是必须的，因为获取的往往是具有数据量大的特点，重建的计算量很大，速度决定了应用效率。因此，我们需要寻找一种重建精度高且速度快的方法。

发明内容

本发明为了解决现有重建方法中不能够精度速度同时提高的问题，而提出了一种压缩感知系统中的高效重建方法。

本发明一种压缩感知系统中的高效重建方法的步骤如下：

步骤一：输入测量值Y0，以及一维信号重建的信息和二维信号重建的信息；

步骤二：根据步骤一输入的一维信号重建的信息和二维信号重建的信息判断是一维信号重建，还是二维信号重建，若是一维信号重建，则执行步骤三，若是二维信号重建，则执行步骤四；

步骤三：记Y＝Y0，并输入测量矩阵Φ，稀疏矩阵Ψ，执行步骤五；

步骤四：对二维信号的测量值Y0执行矢量化操作vec：

$Y=\mathrm{vec}\left(Y0\right)=\mathrm{vec}(\mathrm{\Φ}2*X2)=\stackrel{\‾}{\mathrm{\Φ}2}\mathrm{vec}\left(\mathrm{\Ψ}2S2\mathrm{\Ψ}{3}^T\right)$                       公式五

公式五中，Ψ2为二维信号的行稀疏矩阵，Ψ3为二维信号的列稀疏矩阵，是二维信号的测量矩阵Φ2的变换形式：

$\stackrel{\‾}{\mathrm{\Φ}2}\left(r\right)=\mathrm{vec}\left(\mathrm{\Φ}\left(r\right)\right)$                        公式七

r是二维信号的测量矩阵Φ2(r，j，k)中的一个元素，并且1≤r≤M，

利用Kronecker积将公式五变换为如下形式：

$\mathrm{vec}\left(Y0\right)=\stackrel{\‾}{\mathrm{\Φ}2}(\mathrm{\Ψ}2\⊗\mathrm{\Ψ}{3}^T)\mathrm{vec}\left(S2\right)$                     公式八

令Y＝vec(Y0)，再令 Θ＝ΦΨ，u＝vec(S2)，则公式八等价于公式十一：

Y＝ΦΨu＝Θu                      公式十一

执行步骤五；

步骤五：令k＝1，u^k＝0，v^k＝0，得到v^k+1和u^k+1的公式为：

u^k+1＝δ·shrink(v^k+1，μ)                    公式十二

v^k+1＝v^k+Θ^T(Y-Θu^k)                    公式十三

其中，μ是拉格朗日常数，

δ是固定步长因子，取值范围为：

shrink是个软阈值；则假设，v^k变化s次恰好使得u^k+1有所改变，那么在这些迭代步骤中有如下迭代公式：

u^k+s＝u^k+1                       公式十五

步骤六：判断u^k+1-u^k||≤ε是否成立，其中ε是设定的阈值；若不成立则执行步骤七，若成立则执行步骤八；

步骤七：根据如下公式计算获得v^k+1：

k＝k+1                     公式十六

v^k+1＝v^k+Θ^T(Y-Θu^k)                      公式十七

获得v^k+1后进入步骤十；

步骤八：计算获得v^k的变化次数s：

                          公式十八

s＝min{s_i}i∈I₀                      公式十九

其中I₀表示u^k+1中所有分量为0的下标，即I₁表示u^k+1中所有分量不为0的下标，即

步骤九：将步骤八获得的v^k的变化次数s带入如下公式计算获得v^k+1：

k＝k+1                       公式十六

$v_i^{k+1}=v_i^k+s\·{\left({\mathrm{\Θ}}^T(Y-\mathrm{\Θ}{u}^k)\right)}_i,\∀i\∈I_0$                       公式二十

$v_i^{k+1}=v_i^k,\∀i\∈I_1$                      公式二十一

获得v^k+1后进入步骤十；

步骤十：将获得的v^k+1带入公式十二计算获得u^k+1：

u^k+1＝δshrink(v^k+1，μ)                 公式十二

获得u^k+1后进入步骤十一；

步骤十一：判断是否成立，其中γ是判断收敛的阈值；若成立则执行步骤十二，若不成立则返回步骤六；

步骤十二：稀疏系数u＝u^k+1；

步骤十三：根据步骤一输入的一维信号重建的信息和二维信号重建的信息判断是一维信号重建，还是二维信号重建，若是一维信号重建，则执行步骤十四，若是二维信号重建，则执行步骤十五；

步骤十四：根据稀疏矩阵Ψ，然后再利用公式二十二重建出原始一维信号X：

X＝Ψu                      公式二十二

步骤十五：对系数向量u执行逆矢量化操作ivec，获得二维信号X2的稀疏域表示系数S2：

S2＝ivec(u)                    公式二十三

再利用X2＝Ψ2S2Ψ3重建出原始二维信号X2。

在上述的重建方法中，凸优化方法的重建精度高，其中的Bregman方法相对于其他凸优化算法重建速度也较快。我们的研究主要针对Bregman方法，在它重建精度高的优势基础上，进一步提高其重建速度，实现一种高效的重建方法。由于Bregman方法中最关键的迭代部分具有可以进行优化的潜力，为实现上述目的提供了保障。为了验证本发明算法的性能，我们进行了计算机仿真实验，这里我们使用的测量基是部分傅里叶测量基，使用的稀疏基是db2小波基，我们在两种获取率30％和40％下对同一图像进行仿真。在进行加速算法之前先预设一个判定条件：u^k+1≈u^k是否成立。因此，我们必须设定一个阈值来对条件进行判定即||u^k+1-u^k||≤ε。ε的大小直接影响到迭代所需的时间和迭代精度。此外，还要预设一个判断收敛的阈值γ用以判断来确定迭代是否收敛，此处设γ＝10^-1。实验中，采用的测试图像为一幅标准的1024×1024大小的遥感图像，如图2所示。仿真实验结果如图3至图6所示。由图3至图6可知，在不同的获取率下，本发明所涉及的方法(记为加速后)都比未改进的Bregman方法(记为加速前)迭代时间有着明显地降低，尽管代价是恢复图像的质量有所下降，即信噪比下降。输出时间曲线呈下凸曲线，即随着阈值的变化，时间一开始就下降的很快，只是在阈值变得很大时，时间的下降才缓慢下来。相反，输出信噪比曲线则呈上凸曲线，即随着阈值的变化，输出信噪比一开始下降的很慢，只是在阈值变得很大时，输出信噪比的的下降才会变快。因此，我们可以在满足足够高的信噪比的条件下，选择一个最合适的阈值，在这一阈值下，我们只需最短的迭代时间。也就是说，我们总可以找到这样的一个阈值，使得由于信噪比的降低带来的成本损失通过时间的缩短来挽回。比如，在这里，我们在30％的采样率下，采用阈值ε＝5；而在40％的采样率下，采用阈值ε＝3或者ε＝5。

附图说明

图1本发明的流程图；图2是实验所采用的测试图像；图3是采样率为30％时不同阈值下算法的输出信噪比曲线；图4是采样率为30％时不同阈值下算法的输出信噪比曲线时间曲线；图5是采样率为40％时不同阈值下算法的输出信噪比曲线；图6是采样率为40％时不同阈值下算法的输出信噪比曲线时间曲线。

具体实施方式

具体实施方式一：结合图说明本实施方式，本实施方式的步骤如下：

步骤一：输入测量值Y0，以及一维信号重建的信息和二维信号重建的信息；

步骤二：根据步骤一输入的一维信号重建的信息和二维信号重建的信息判断是一维信号重建，还是二维信号重建，若是一维信号重建，则执行步骤三，若是二维信号重建，则执行步骤四；

步骤三：记Y＝Y0，并输入测量矩阵Φ，稀疏矩阵Ψ，执行步骤五；

步骤四：对测量值Y0执行矢量化操作vec，记Y＝vec(Y0)：

如果考虑到压缩感知过程中测量矩阵中系数的随机性，只要二维信号的测量矩阵Φ2不是非结构性，那么对测量值Y0执行矢量化操作vec：

$Y=\mathrm{vec}\left(Y0\right)=\mathrm{vec}(\mathrm{\Φ}2*X2)=\stackrel{\‾}{\mathrm{\Φ}2}\mathrm{vec}\left(\mathrm{\Ψ}2S2\mathrm{\Ψ}{3}^T\right)$                 公式五

公式五中，矢量化操作vec的含义如下：

设A＝(a_i，j)_M×N，称M×N维的列向量vec(A)为矩阵A的按行展开，或矩阵A按行拉直的列向量，如下式：

vec(A)＝(a₁₁，…，a₂₁，…，a_m1，…，a_mn)^T            公式六

公式五中，是二维信号的测量矩阵Φ2的变换形式(维的矩阵，)，如果令中的每一个元素为1≤r≤M，，令二维信号的测量矩阵Φ2中的每一个元素为Φ2(r，j，k)，1≤r≤M，，1≤j≤N1，1≤k≤N2，那么是由二维信号的测量矩阵Φ2确定的，关系如下：

$\stackrel{\‾}{\mathrm{\Φ}2}\left(r\right)=\mathrm{vec}\left(\mathrm{\Φ}\left(r\right)\right)$                公式七

利用Kronecker积将公式五变换为如下形式：

$\mathrm{vec}\left(Y0\right)=\stackrel{\‾}{\mathrm{\Φ}2}(\mathrm{\Ψ}2\⊗\mathrm{\Ψ}{3}^T)\mathrm{vec}\left(S2\right)$                 公式八

Kronecker积的具体解释如下：

根据矩阵分析中的理论，设A∈C^m×n，X∈C^n×p，B∈C^p×q，那么则有下式成立

$\mathrm{vec}\left(\mathrm{AXB}\right)=(A\⊗B^T)\mathrm{vecX}$              公式九

其中，表示矩阵A和矩阵B的Kronecker积，

其定义如公式十，其中设A＝(a_i，j)∈C^m×n，B＝(b_i，j)∈C^p×q。

$A\⊗B=\left[\begin{array}{cccc}{a}_{11}B& a_{12}B& \·\·\·& a_{1n}B\\ a_{21}B& a_{22}B& \·\·\·& a_{2n}B\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ a_{m1}B& a_{m2}B& \·\·\·& a_{\mathrm{mn}}B\end{array}\right]=\left(a_{\mathrm{ij}}B\right)\∈C^{\mathrm{mp}\×\mathrm{nq}}$            公式十

令Y＝vec(Y0)，再令 Θ＝ΦΨ，u＝vec(S2)，则公式八等价于公式十一：

Y＝ΦΨu＝Θu               公式十一

执行步骤五；

步骤五：令k＝1，u^k＝0，v^k＝0，得到v^k+1和u^k+1的公式为：

u^k+1＝δ·shrink(v^k+1，μ)                 公式十二

v^k+1＝v^k+Θ^T(Y-Θu^k)                公式十三

其中，μ是拉格朗日常数，应使μ在足够小的范围内越大，收敛越快，

δ是固定步长因子，取值范围为：

shrink是个软阈值：

$\mathrm{shink}(y,\mathrm{\α})=\mathrm{sign}\left(y\right)\·\max \left\{\right|y|-\mathrm{\α},0\}=\left\{\begin{array}{cc}y-\mathrm{\α},& y\∈(\mathrm{\α},\∞)\\ 0& y\∈[-\mathrm{\α},\mathrm{\α}]\\ y+\mathrm{\α},& y\∈(-\∞,-\mathrm{\α})\end{array}\right.$             公式十四

在上述迭代中，由于shrink软阈值操作的影响，当出现v∈[-μ，μ]时，u不发生变化，只有v出了[-μ，μ]范围内时，u才发生变化。也即是说u将要保持不变直到v迭代使u出现一个非零值之后。而在此期间v的迭代会占用一定的时间，从而影响整体重建方法的速度，这个期间成为无贡献迭代。这里设I₀表示u^k+1中所有分量为0的下标，即同理，我们设I₁表示u^k+1中所有分量不为0的下标，即对于I₀中的所有值，由于shrink操作的存在，可能导致多个v^k的不断变化，均使得u^k+1没有任何改变，这样在计算v^k时，无需每次都计算一次Θ^T(Y-Θu^k)，只需计算一次即可，这样就可以对线性算法进行加速。假设，v^k变化s次恰好使得u^k+1有所改变，那么在这些迭代步骤中有如下迭代公式：

u^k+s＝u^k+1                 公式十五

步骤六：判断||u^k+1-u^k||≤ε是否成立，其中ε是设定的阈值，大小直接影响到迭代所需的时间和迭代精度；若不成立则执行步骤七，若成立则执行步骤八；

步骤七：根据如下公式计算获得v^k+1：

k＝k+1                  公式十六

v^k+1＝v^k+Θ^T(Y-Θu^k)                    公式十七

获得v^k+1后进入步骤十；

步骤八：计算获得s：

           公式十八

s＝min{s_i}i∈I₀              公式十九

其中I₀表示u^k+1中所有分量为0的下标，即I₁表示u^k+1中所有分量不为0的下标，即

步骤九：将步骤八获得的s带入如下公式计算获得v^k+1：

k＝k+1                    公式十六

$v_i^{k+1}=v_i^k+s\·{\left({\mathrm{\Θ}}^T(Y-\mathrm{\Θ}{u}^k)\right)}_i,\∀i\∈I_0$            公式二十

$v_i^{k+1}=v_i^k,\∀i\∈I_1$                公式二十一

获得v^k+1后进入步骤十；

步骤十：将获得的v^k+1带入公式十二计算获得u^k+1：

u^k+1＝δshrink(v^k+1，μ)                   公式十二

获得u^k+1后进入步骤十一；

步骤十一：判断是否成立，其中γ是判断收敛的阈值；若成立则执行步骤十二，若不成立则返回步骤六；

步骤十二：稀疏系数u＝u^k+1；

步骤十三：根据步骤一输入的一维信号重建的信息和二维信号重建的信息判断是一维信号重建，还是二维信号重建，若是一维信号重建，则执行步骤十四，若是二维信号重建，则执行步骤十五；

步骤十四：根据稀疏矩阵Ψ，然后再利用公式一重建出原始一维信号X：

X＝Ψu                         公式二十二

步骤十五：对系数向量u执行逆矢量化操作ivec，获得二维信号X2的稀疏域表示系数S2：

S2＝ivec(u)                    公式二十三

再利用X2＝Ψ2S2Ψ3重建出原始二维信号X2。

Claims (1)

1.一种压缩感知系统中的高效重建方法，其特征在于它的步骤如下：

步骤一：输入测量值Y0，以及一维信号重建的信息和二维信号重建的信息；

步骤二：根据步骤一输入的一维信号重建的信息和二维信号重建的信息判断是一维信号重建，还是二维信号重建，若是一维信号重建，则执行步骤三，若是二维信号重建，则执行步骤四；

步骤三：记Y=Y0，并输入测量矩阵Φ，稀疏矩阵Ψ，执行步骤五；

步骤四：对二维信号的测量值Y0执行矢量化操作vec：

$Y=\mathrm{vec}\left(Y0\right)=\mathrm{vec}(\mathrm{\Φ}2*X2)=\stackrel{\‾}{\mathrm{\Φ}2}\mathrm{vec}\left(\mathrm{\Ψ}2S2\mathrm{\Ψ}{3}^T\right)$    公式五

公式五中，Ψ2为二维信号的行稀疏矩阵，Ψ3为二维信号的列稀疏矩阵，是二维信号的测量矩阵Φ2的变换形式：

$\stackrel{\‾}{\mathrm{\Φ}2}\left(r\right)=\mathrm{vec}\left(\mathrm{\Φ}\left(r\right)\right)$    公式七

r是二维信号的测量矩阵Φ2中的一个元素，并且1≤r≤M，

利用Kronecker积将公式五变换为如下形式：

$\mathrm{vec}\left(U0\right)=\stackrel{\‾}{\mathrm{\Φ}2}(\mathrm{\Ψ}2\⊗{\mathrm{\Ψ}3}^T)\mathrm{vec}\left(S2\right)$    公式八

令Y=vec(Y0)，再令 $\mathrm{\Φ}=\stackrel{\‾}{\mathrm{\Φ}2},\mathrm{\Ψ}=\mathrm{\Ψ}2\⊗{\mathrm{\Ψ}3}^T,\mathrm{\Θ}=\mathrm{\Φ\Ψ},u=\mathrm{vec}\left(S2\right),$ 则公式八等价于公式十一：

Y=ΦΨu=Θu   公式十一

执行步骤五；

步骤五：令k=1,u^k=0,v^k=0，得到v^k+1和u^k+1的公式为：

u^k+1＝δ·shrink(v^k+1，μ)   公式十二

v^k+1＝v^k+Θ^T(Y-Θu^k)   公式十三

其中，μ是拉格朗日常数，

δ是固定步长因子，取值范围为：

shrink是个软阈值；则假设，v^k变化s次恰好使得u^k+1有所改变，那么在这些迭代步骤中有如下迭代公式：

u^k+s=u^k+1   公式十五

步骤六：判断||u^k+1-u^k||≤ε是否成立，其中ε是设定的阈值；若不成立则执行步骤七，若成立则执行步骤八；

步骤七：根据如下公式计算获得v^k+1：

k=k+1   公式十六

v^k+1=v^k+Θ^T(Y-Θu^k)   公式十七

获得v^k+1后进入步骤十；

步骤八：计算获得v^k的变化次数s：

   公式十八

s=min{s_i}i∈I₀   公式十九

其中I₀表示u^k+1中所有分量为0的下标，即I₁表示u^k+1中所有分量不为0的下标，即 $I_1=\left\{i\right|u_i^{k+1}\≠0\};$

步骤九：将步骤八获得的v^k的变化次数s带入如下公式计算获得v^k+1：

k=k+1   公式十六

$v_i^{k+1}=v_i^k+s\·{\left({\mathrm{\Θ}}^T(Y-{\mathrm{\Θu}}^k)\right)}_i,\∀i\∈I_0$    公式二十

$v_i^{k+1}=v_i^k,\∀i\∈I_1$    公式二十一

获得v^k+1后进入步骤十；

步骤十：将获得的v^k+1带入公式十二计算获得u^k+1：

u^k+1=δshrink(v^k+1,μ)   公式十二

获得u^k+1后进入步骤十一；

步骤十一：判断是否成立，其中γ是判断收敛的阈值；若成立则执行步骤十二，若不成立则返回步骤六；

步骤十二：稀疏系数u=u^k+1；

步骤十三：根据步骤一输入的一维信号重建的信息和二维信号重建的信息判断是一维信号重建，还是二维信号重建，若是一维信号重建，则执行步骤十四，若是二维信号重建，则执行步骤十五；

步骤十四：根据稀疏矩阵Ψ，然后再利用公式二十二重建出原始一维信号X：

X=Ψu   公式二十二

步骤十五：对系数向量u执行逆矢量化操作ivec，获得二维信号X2的稀疏域表示系数S2：

S2=ivec(u)   公式二十三

再利用X2=Ψ2S2Ψ3重建出原始二维信号X2。