CN102722896B - 基于自适应压缩感知的自然图像非局部重构方法 - Google Patents

Abstract

本发明公开一种基于自适应压缩感知的自然图像非局部重构方法，主要解决已有技术重构图像信息丢失严重等问题。其实现步骤为：(1)把一幅图像分成N个32×32的子块，根据基本采样率b和感知矩阵Φ得到基本感知矩阵Φ′，利用Φ′对信号采样，得到基本观测向量(2)根据估计出图像的标准差序列{d₁,d₂,....d_N}；(3)根据标准差序列{d₁,d₂,....d_N}，为每个子块自适应的分配采样率a_i并构造自适应感知矩阵利用对信号采样，得到自适应观测向量（4）用基本观测向量与自适应观测向量共同组成子块的观测向量(5)根据观测向量得到图像的初始解x⁰；（6）用x⁰进行迭代，重构原始图像，直至满足终止条件，得到重构图像x′。本发明具有图像重构质量高，原理清晰，操作简单的优点，适用于自然图像的采样和重构。

Description

基于自适应压缩感知的自然图像非局部重构方法

技术领域

本发明属于数字图像处理领域，特别涉及一种自适应的采样策略及非局部的重构方法，用于对自然图像的存储、传输和处理。

背景技术

压缩感知理论是新兴的一种信号采样策略，它将信号的采样过程与压缩过程成功的融为了一体。压缩感知的主要思想是：如果一个信号是稀疏的或者是可压缩的。就可以利用一个观测矩阵将信号投影到一个低维空间上，其中低维空间中的每个观测值都包含了整体信号的信息，它们之间具有很小的冗余性，这样就可以根据少量的观测值，通过求解一个凸规划的优化问题将信号准确的重构回来。显然，在压缩感知理论中，信号的采样与压缩的过程同时以一个较低的速率进行着，使得采样的成本大大的降低。它开创了一个更经济有效的模拟信号数字化的途径。对于任意信号来说，只要能够找到其相应一个稀疏表示的空间，就可以应用压缩感知的理论进行采样与重构。

对于一个二维图像信号来说，由于其自身的高维度，使得压缩感知理论不得不面对这样一个“维数灾难”的问题。感知矩阵的存贮与计算将耗费大量的资源，这将对构造图像信号的实时采样系统带来极大的困难。另一方面，图像的重构过程也同样的面临极高的运算复杂度。针对这一问题，L.Gan提出了二维自然图像的分块压缩感知的算法：将原始图像分为数个尺寸更小的图像块，然后用相同的观测器分别对每一个图像块进行采样。这样的采样方式有三个好处，(1)感知矩阵的规模将大大的缩减，很大程度上减小了采样和重构过程中的计算复杂度；(2)采得的信号更具有实时性，不必等到整幅图像都采样结束后再进行传输；(3)由于每一个分块都进行的独立的采样，很容易得到一个初始解并且信号的重构速度也将大大的加快。

但是在现有的算法中，对于每个子块都使用同样的感知器去观测，即对每一个子块都分配同样的采样率，无论该子块的稀疏度如何。而实际上，由于图像具有不同的结构信息，它的每一个子块的稀疏度都不尽相同。对于一个稀疏度比较高的子块，一个较低的采样率sr，sr＝M/N足以将该子块重构出一个很好的效果，如果使用较高的采样率则会造成采样资源的浪费；而对于一个稀疏度比较低得子块，比如该子块中包含较多的纹理信息等，就需要一个较高的采样率来对该子块进行观测，否则会造成大量的结构信息的丢失。

此外，在信号重构的过程中，迭代的“投影-滤波”方法被广泛的使用。在滤波器的选择上，由于操作便捷，局部滤波的方法被广泛的采用。然而，局部滤波方法其自身有着一定的缺陷，例如会模糊图像的边缘，对噪声去除的不彻底等，大大的影响了图像的重构效果。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于自适应压缩感知的自然图像非局部重构方法，以避免结构信息的丢失和采样资源的浪费，降低图像边缘的模糊和噪声的影响，提高图像的重构效果。

为实现上述目的，本发明包括如下步骤:

(1)将输入的图像信号x分成N个32×32大小的子块x₁,x₂,...,x_N，给出平均采样率s，基本采样率b和感知矩阵Φ，根据基本采样率b和感知矩阵Φ得到基本感知矩阵Φ′，利用基本感知矩阵Φ′对每个图像子块x_i进行采样，得到每个图像子块的基本观测向量：其中i＝1,2,....N，N为图像子块的个数；

(2)根据基本观测向量估计出图像的标准差序列{d₁,d₂，....d_N}，其中N为图像子块的个数；

(3)根据标准差序列{d₁,d₂,....d_N}，为每一个图像子块自适应的分配一个采样率a_i，i＝1,2,....N，根据自适应采样率和感知矩阵Φ构造自适应感知矩阵利用自适应感知矩阵对每个图像子块进行采样，得到每个图像子块的自适应观测向量：

(4)将基本观测向量与自适应观测向量存放在一个列向量中，组成每个图像子块观测向量 $y_i=\left(\begin{array}{c}{y}_{b_i}\\ y_{a_i}\end{array}\right);$

(5)将每个图像子块的观测向量y_i按列存放在一个向量中，构成整幅图像的观测向量 $y=\left(\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_N\end{array}\right),$ 根据整幅图像观测向量y得到图像信号x的初始解：

$x^0=(x_1^0,x_2^0,.....x_N^0),$

其中 $x_i^0={x^{\′}}_i+{\mathrm{\Φ}}_i^T(y_i-{\mathrm{\Φ}}_i{x^{\′}}_i),$ i＝1,2，....N， ${\mathrm{\Φ}}_i=\left(\begin{array}{c}{\mathrm{\Φ}}^{\′}\\ {\mathrm{\Φ}}_{a_i}\end{array}\right)$ 为每个图像子块对应的感知矩阵，T表示对感知矩阵转置， 为估计出的每个图像子块的均值，i＝1,2,....N，g×z是基本采样矩阵Φ′的维度，var为求解序列的方差；

(6)对初始解x⁰进行迭代，重构出原始信号x′：

6a)设定k为迭代次数，令k＝0；

6b)对上一次迭代图像x^k进行小波变换，得到小波系数：ξ^k＝Ψx^k，Ψ为小波变换基，然后对小波系数ξ^k进行双变量阈值平滑，得到阈值平滑后的小波系数：

${\mathrm{\ξ}}^{k+1}=\frac{{(\sqrt{{\left({\mathrm{\ξ}}^k\right)}^2+{\mathrm{\ξ}}_p^2}-p\frac{\sqrt{{3\mathrm{\σ}}^{\left(i\right)}}}{{\mathrm{\σ}}_{\mathrm{\ξ}}})}_{+}}{\sqrt{{\left({\mathrm{\ξ}}^k\right)}^2+{\mathrm{\ξ}}_p^2}}\·{\mathrm{\ξ}}^k,$

其中ρ为双变量阈值参数，ξ_p为父节点系数，σ⁽ⁱ⁾为子系数的度量方差，σ_ξ为估计出来的图像的方差，+表示值比0大时取本身，比0小时取0，再对系数ξ^k+1进行逆变换得到双变量阈值平滑后的图像：Ψ^-1为小波逆变换基；

6c)对上一次迭代图像x^k进行非局部总变差平滑处理，得到处理后的图像

6d)计算双变量阈值平滑后的图像和非局部总变差平滑处理后的图像的加权均值，得到均值图像：

${\stackrel{\widehat{\‾}}{x}}^{k+1}=c_1{\hat{x}}^{k+1}+c_2{\stackrel{\‾}{x}}^{k+1},$ 其中0＜c₁＜1，0＜c₂＜1，c₁+c₂＝1

将均值图像分成N个32×32大小的子块对每个子块进行投影处理，根据观测值y_i的长度寻找出每个子块对应的感知矩阵Φ_i，将投影到超平面{e|Φ_ie＝y_i,e∈R^32×32}上，e是一个变量，R^32×32表示32×32维的信号，得到投影处理后的图像子块：i＝1,2,...,N，这些投影处理后的图像子块组成迭代后的图像x^k+1；

6e)计算信号在每次迭代后与迭代前的变化情况E(k)＝||x^k+1-x^k||₂，||.||₂表示2范数，如果|E(k)-E(k-1)|≤0.001或者k＞100，迭代终止，重构的原始信号x′＝x^k+1，否则迭代次数k加1，返回步骤6b)。

本发明相比现有技术具有以下优点：

(1)本发明使用了自适应的图像压缩采样方法，针对不同稀疏度的图像子块，采用不同的采样策略，有效的保留了图像信号的大量信息；

(2)本发明将非局部的思想引入到了图像信号的重构过程中来，有效的保护了图像信号的边缘与纹理信息，使重构出的图像信号得到了更高的峰值信噪比，提高了图像的重构效果；

(3)本发明推导出了仅使用观测向量来估计图像信号的标准差的方法；

(4)本发明推导出使用增广拉格朗日法来解决非局部正则化问题的方法；

(5)本发明提出基于标准差最小的图像信号的初始化方法。

附图说明

图1是本发明的总流程图；

图2是本发明中的自适应采样子流程图；

图3是本发明使用基本观测值对信号均值估计的结果图；

图4是本发明使用基本观测值对信号标准差估计的结果图；

图5是本发明中的信号重构子流程图；

图6是用本发明对图像进行初始化的结果图；

图7是用本发明和BCS-SPL-DWT方法得到的重构图像的细节对比图；

图8是用本发明和BCS-SPL-DWT方法得到的重构图像的峰值信噪比及结构相似度的曲线图。

具体实施方式

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

步骤1，将输入的图像信号x分成N个32×32大小的子块x₁,x₂,...,x_N，并对每个子块进行基本采样：

1a）给出平均采样率s，基本采样率b和感知矩阵Φ，根据基本采样率b计算出基本采样行数M＝N_x×s，其中N_x＝1024为信号子块维数，从感知矩阵Φ中取出其前M行构成基本感知矩阵Φ′。

1b)使用基本感知矩阵Φ′对每个图像子块x_i进行采样，得到每个图像子块的基本观测向量：其中i＝1,2,....N，N为图像子块的个数。

步骤2，根据基本观测向量估计出图像的标准差序列{d₁,d₂，....d_N}，其中i＝1,2,...,N，g×z是基本采样矩阵Φ′的维度，var为求解序列的方差，N为图像子块的个数，对标准差的估计效果如图4所示，其中图4（a）为Lena图像估计的标准差效果图，图4（b）为Barbara图像估计的标准差效果图。

步骤3，为每个图像子块进行自适应采样。

参照图2，为了有效的利用采样资源，根据每个图像子块的标准差，为其自适应的分配采样率并进行如下自适应采样：

3a）根据标准差序列{d₁,d₂,....d_N}，为每一个图像子块自适应的分配一个采样率a_i，其中i＝1,2,...,N，s为平均采样率，sum(d)为所有图像子块标准差的和，N为图像子块的个数；

3b）根据自适应采样率a_i计算出自适应采样行数M_i＝N_x×a_i，从采样矩阵Φ中取出其前M_i行构成自适应感知矩阵i＝1,2,....N，N_x＝1024为图像子块的维数。

3c）利用自适应感知矩阵对信号进行采样，得到每个图像子块的自适应观测向量：其中i＝1,2,....N。

步骤4，将基本观测向量与自适应观测向量存放在一个列向量中，组成每个图像子块观测向量 $y_i=\left(\begin{array}{c}{y}_{b_i}\\ y_{a_i}\end{array}\right).$

步骤5，对图像进行初始化，得到图像的初始解：

将每个图像子块的观测向量y_i按列存放在一个向量中，构成整幅图像的观测向量 $y=\left(\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_N\end{array}\right),$ 根据整幅图像观测向量y得到图像信号x的初始解：

$x^0=(x_1^0,x_2^0,.....x_N^0),$

图像初始化结果如图6所示，其中i＝1，2,....N， ${\mathrm{\Φ}}_i=\left(\begin{array}{c}{\mathrm{\Φ}}^{\′}\\ {\mathrm{\Φ}}_{a_i}\end{array}\right)$ 为每个图像子块对应的感知矩阵，T表示对矩阵求转置， 为估计出的每个图像子块的均值，估计的均值效果如图3所示，其中图3（a）为Lena图像估计的均值效果图，图3（b）为Barbara图像估计的均值效果图,i＝1,2,....N，g×z是基本采样矩阵Φ′的维度，var为求解序列的方差。

步骤6：根据初始解对图像信号进行重构。

参照图5，为了得到更好的重构效果，引入非局部的思想，联合双变量阈值方法来重构图像，实现步骤如下：

6a)设定k为迭代次数，令k＝0；

6b)对上一次迭代图像x^k进行小波变换，得到小波系数：ξ^k＝Ψx^k，Ψ为小波变换基，然后对小波系数ξ^k进行双变量阈值平滑，得到阈值平滑后的小波系数：

${\mathrm{\ξ}}^{k+1}=\frac{{(\sqrt{{\left({\mathrm{\ξ}}^k\right)}^2+{\mathrm{\ξ}}_p^2}-p\frac{\sqrt{{3\mathrm{\σ}}^{\left(i\right)}}}{{\mathrm{\σ}}_{\mathrm{\ξ}}})}_{+}}{\sqrt{{\left({\mathrm{\ξ}}^k\right)}^2+{\mathrm{\ξ}}_p^2}}\·{\mathrm{\ξ}}^k,$

其中ρ为双变量阈值参数，ξ_p为父节点系数，σ⁽ⁱ⁾为子系数的度量方差，σ_ξ为估计出来的图像的方差，+表示值比0大时取本身，比0小时取0，再对系数ξ^k+1进行逆变换得到双变量阈值平滑后的图像：Ψ^-1为小波逆变换基；

6c)对上一次迭代图像x^k进行非局部总变差平滑处理，得到处理后的图像

6c1）计算x^k的每个像素在搜索框内的非局部权值w_m，n，m，n为图像的坐标值；

6c2）设定优化目标函数 ${\stackrel{\‾}{x}}_{k+1}=\underset{h}{\min }\∫{\left|\right|{\▿}_{\mathrm{NL}}h\left|\right|+\mathrm{\μ}\left|\right|h-x^k\left|\right|}_2^2,$ 其中，h是一个变量，μ为一个调节参数，为非局部梯度，使用增广拉格朗日模型求解上述优化问题，按如下步骤进行：

（1）设定l为迭代次数，令l＝0，l_max＝4，拉格朗日算子λ_m，n，λ_n,m的初始值为1，中间变量f_m，n，f_n,m的初始值为1，m，n为图像的坐标值；

（2）用如下公式对上一次迭代图像x^k的每个像素点进行非局部总变差处理，得到处理后的图像的像素点

${\left(x_m^k\right)}^{l+1}=\frac{{{\mathrm{\μx}}_m^k}^{\′}+2r\underset{n}{\mathrm{\Σ}}{w}_{m,n}{x_n^k}^{\′}-\underset{n}{\mathrm{\Σ}}\sqrt{w_{m,n}}({\mathrm{\λ}}_{m,n}^l-{\mathrm{\λ}}_{n,m}^l)-r\underset{n}{\mathrm{\Σ}}\sqrt{w_{m,n}}(f_{m,n}^l-f_{n,m}^l)}{\mathrm{\μ}+2r\·\underset{n}{\mathrm{\Σ}}{w}_{m,n}}$

其中，为中间变量，

${x_m^k}^{\′}=\frac{\mathrm{\μ}{\left(x_m^k\right)}^l+2r\underset{n}{\mathrm{\Σ}}{w}_{m,n}{\left(x_n^k\right)}^l-\underset{n}{\mathrm{\Σ}}\sqrt{w_{m,n}}({\mathrm{\λ}}_{m,n}^l-{\mathrm{\λ}}_{n,m}^l)-r\underset{n}{\mathrm{\Σ}}\sqrt{w_{m,n}}(f_{m,n}^l-f_{n,m}^l)}{\mathrm{\μ}+2r\·\underset{n}{\mathrm{\Σ}}{w}_{m,n}},$ μ为一个调节参数，r是一个正常数，为拉格朗日算子，w_m,n为非局部权值， 为中间变量，表示x^k的像素点，m，n为图像的坐标值；

（3）对中间变量以及拉格朗日算子的值进行更新，得到更新后的中间变量和拉格朗日算子

$f_{m,n}^{l+1}=\left\{\begin{array}{cc}\frac{1}{r}(1-\frac{1}{\left|{\mathrm{\ω}}_{m,n}\right|}){\mathrm{\ω}}_{m,n},& \mathrm{if}\left|{\mathrm{\ω}}_{m,n}\right|>1\\ 0,& \mathrm{otherwise}\end{array}\right.,$

${\mathrm{\λ}}_{m,n}^{l+1}={\mathrm{\λ}}_{m,n}^l+r(f_{m,n}^l-{\▿}_{\mathrm{NL}}{\left(x_m^k\right)}^{l+1})$

其中 ${\mathrm{\ω}}_{m,n}=r\sqrt{w_{m,n}}(x_m^{l+1}-x_n^{l+1})-{\mathrm{\λ}}_{m,n}^l,$ w_m，n为非局部权值，r是一个正常数，为非局部梯度；

（4）如果l＞l_max，迭代终止， 为平滑处理后图像的像素点，m为图像的坐标值，否则迭代次数l加1，返回步骤（2）；

6d)计算双变量阈值平滑后的图像和非局部总变差平滑处理后的图像的加权均值，得到均值图像：

${\stackrel{\widehat{\‾}}{x}}^{k+1}=c_1{\hat{x}}^{k+1}+c_2{\stackrel{\‾}{x}}^{k+1},$ 其中0＜c₁＜1，0＜c₂＜1，c₁+c₂＝1

将均值图像分成N个32×32大小的子块对每个子块进行投影处理，根据观测值y_i的长度寻找出每个子块对应的感知矩阵Φ_i，将投影到超平面{e|Φ_ie＝y_i,e∈R^32×32}上，e是一个变量，R^32×32表示32×32维的信号，得到投影处理后的图像子块：i＝1,2,...,N，这些投影处理后的图像子块组成迭代后的图像x^k+1；

6e)计算信号在每次迭代后与迭代前的变化情况E(k)＝||x^k+1-x^k||₂，||.||₂表示2范数，如果|E(k)-E(k-1)|≤0.001或者k＞100，迭代终止，重构的原始信号x′＝x^k+1，否则迭代次数k加1，返回步骤6b)。

本发明的效果可以通过下面的仿真结果进一步说明。

1．仿真条件

用本发明对2幅自然图像进行仿真实验。

在两幅512×512大小的图像上进行测试，它们分别为Lena图像，Barbara图像。图像采样时所用的基本采样率b＝0.4×s，s为采样率。

本实验的目的是展示在相同采样率的情况下，本发明与现有BCS-SPL-DWT方法对自然图像的重构结果；在不同采样率情况下，本发明与现有BCS-SPL-DWT方法对自然图像重构的峰值信噪比及结构相似度的变化趋势。

仿真1，用本发明与现有BCS-SPL-DWT方法分别在采样率s=0.1,0.2,0.3,0.4和0.5的情况下，对大小为512×512的标准自然图像Barbara进行重构，其结果如图7所示，其中图7（a）是BCS-SPL-DWT方法在各个采样率下的细节重构结果图，图7（b）是本发明在各个采样率下的细节重构结果图。从图7可以看出，在各个采样率下，本发明能够很好的将图像的细节和纹理重构出来。

仿真2，用本发明与现有BCS-SPL-DWT方法分别在不同采样率下，对两幅大小为512×512的标准自然图像Lena，Barbara进行重构，其重构结果的峰值信噪比和结构相似度的曲线图如图8所示，其中图8（a）是Lena图像重构结果的峰值信噪比的曲线图，8（b）是Lena图像重构结果的结构相似度的曲线图，8（c）是Barbara图像重构结果的峰值信噪比的曲线图，8（d）是Barbara图像重构结果的结构相似度的曲线图。从图8可以看出，本发明在不同采样率下重构结果的峰值信噪比和结构相似度均远高于BCS-SPL-DWT方法，其中峰值信噪比最多提升了4.13db。

综上所述，本发明不仅有效的保留了图像信号的大量信息，而且能够很好的将图像的细节与纹理重构出来，特别适合于对自然图像的采样和重构。

Claims (6)

1.一种基于自适应压缩感知的自然图像非局部重构方法，包括如下步骤：

(1)将输入的图像信号x分成N个32×32大小的子块x₁,x₂,...,x_N，给出平均采样率s，基本采样率b和感知矩阵Φ，根据基本采样率b和感知矩阵Φ得到基本感知矩阵Φ'，利用基本感知矩阵Φ'对每个图像子块x_i进行采样，得到每个图像子块的基本观测向量：其中i＝1,2,....N，N为图像子块的个数；

(2)根据基本观测向量估计出图像的标准差序列{d₁,d₂,....d_N}，其中N为图像子块的个数；

(3)根据标准差序列{d₁,d₂,....d_N}，为每一个图像子块自适应的分配一个采样率a_i，i＝1,2,....N，根据自适应采样率和感知矩阵Φ构造自适应感知矩阵，利用自适应感知矩阵对每个图像子块进行采样，得到每个图像子块的自适应观测向量：

(4)将基本观测向量与自适应观测向量存放在一个列向量中，组成每个图像子块观测向量 $y_i=\left(\begin{array}{c}{y}_{b_i}\\ y_{a_i}\end{array}\right);$

(5)将每个图像子块的观测向量y_i按列存放在一个向量中，构成整幅图像的观测向量 $y=\left(\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_N\end{array}\right),$ 根据整幅图像观测向量y得到图像信号x的初始解：

$x^0=(x_1^0,x_2^0,...x_N^0),$

其中 $x_i^0={x^{\′}}_i+{\mathrm{\Φ}}_i^T(y_i-{\mathrm{\Φ}}_i{x^{\′}}_i),i=\mathrm{1,2},....,N,{\mathrm{\Φ}}_i=\left(\begin{array}{c}{\mathrm{\Φ}}^{\′}\\ {\mathrm{\Φ}}_{a_i}\end{array}\right)$ 为每个图像子块对应的感知矩阵，T表示对感知矩阵转置， ${x^{\′}}_i=\{{\stackrel{\‾}{x}}_1,{\stackrel{\‾}{x}}_2,...,{\stackrel{\‾}{x}}_N\},{\stackrel{\‾}{x}}_i=\sqrt{\frac{g-1}{z}\·\mathrm{var}\left(y_{a_i}\right)}$ 为估计出的每个图像子块的均值，i＝1,2,....N，基本采样矩阵Φ'的维数为g×z，g为矩阵行数，z为矩阵列数，var为求解序列的方差；

(6)对初始解x⁰进行迭代，重构出原始信号x'：

6a)设定k为迭代次数，令k＝0；

6b)对上一次迭代图像x^k进行小波变换，得到小波系数：ξ^k＝Ψx^k，Ψ为小波变换基，然后对小波系数ξ^k进行双变量阈值平滑，得到阈值平滑后的小波系数：

${\mathrm{\ξ}}^{k+1}=\frac{{(\sqrt{{\left({\mathrm{\ξ}}^k\right)}^2+{\mathrm{\ξ}}_p^2}-\mathrm{\ρ}\frac{\sqrt{3{\mathrm{\σ}}^{\left(i\right)}}}{{\mathrm{\σ}}_{\mathrm{\ξ}}})}_{+}}{\sqrt{{\left({\mathrm{\ξ}}^k\right)}^2+{\mathrm{\ξ}}_p^2}}\·{\mathrm{\ξ}}^k,$

其中ρ为双变量阈值参数，ξ_p为父节点系数，σ⁽ⁱ⁾为子系数的度量方差，σ_ξ为估计出来的图像的方差，+表示值比0大时取本身，比0小时取0，再对系数ξ^k+1进行逆变换得到双变量阈值平滑后的图像：Ψ^-1为小波逆变换基；

6c)对上一次迭代图像x^k进行非局部总变差平滑处理，得到处理后的图像

6d)计算双变量阈值平滑后的图像和非局部总变差平滑处理后的图像的加权均值，得到均值图像：

其中0＜c₁＜1，0＜c₂＜1，c₁+c₂＝1

将均值图像分成N个32×32大小的子块对每个子块进行投影处理，根据观测值y_i的长度寻找出每个子块对应的感知矩阵Φ_i，将投影到超平面{e|Φ_ie＝y_i,e∈R^32×32}上，e是一个变量，R^32×32表示32×32维的信号，得到投影处理后的图像子块： $x_i^{k+1}={\widehat{\stackrel{\‾}{x}}}_i^{k+1}+{\mathrm{\Φ}}_i^T(y_i-{\mathrm{\Φ}}_i{\widehat{\stackrel{\‾}{x}}}_i^{k+1}),i=\mathrm{1,2},...,N,$ 这些投影处理后的图像子块组成迭代后的图像x^k+1；

6e)计算信号在每次迭代后与迭代前的变化情况E(k)＝||x^k+1-x^k||₂，||.||₂表示2范数，如果E(k)-E(k-1)≤0.001或者k＞100，迭代终止，重构的原始信号x'＝x^k+1，否则迭代次数k加1，返回步骤6b)。

2.根据权利要求1所述的基于自适应压缩感知的自然图像非局部重构方法，其中步骤(1)所述的根据基本采样率b和感知矩阵Φ得到基本感知矩阵Φ'，是根据基本采样率b计算出基本采样行数M＝N_x×b，从感知矩阵Φ中取出前M行构成基本感知矩阵Φ'，其中N_x＝1024为图像子块的维数。

3.根据权利要求1所述的基于自适应压缩感知的自然图像非局部重构方法，其中步骤(2)所述的根据基本观测向量估计出图像标准差序列{d₁,d₂,....d_N}，是通过下式进行的：

其中i＝1,2,...,N，基本采样矩阵Φ'的维数为g×z，g为矩阵行数，z为矩阵列数，var为求解序列的方差，N为图像子块的个数。

4.根据权利要求1所述的基于自适应压缩感知的自然图像非局部重构方法，其中步骤(3)所述的根据标准差序列{d₁,d₂,....d_N}，为每一个图像子块自适应的分配一个采样率a_i，是通过下式进行的：

其中i＝1,2,...,N，s为平均采样率，sum(d)为所有图像子块标准差的和，N为图像子块的个数。

5.根据权利要求1所述的基于自适应压缩感知的自然图像非局部重构方法，其中步骤(3)所述的根据自适应采样率和感知矩阵Φ构造自适应感知矩阵是根据自适应采样率a_i计算出自适应采样行数M_i＝N_x×a_i，从采样矩阵Φ中取出其前M_i行构成自适应感知矩阵i＝1,2,....N，N_x＝1024为图像子块的维数。

6.根据权利要求1所述的基于自适应压缩感知的自然图像非局部重构方法，其中步骤6c)中所述的对上一次迭代图像x^k进行非局部总变差平滑处理，是按如下步骤进行：

6c1）计算x^k的每个像素在搜索窗口的非局部权值w_m,n，m，n为图像的坐标值；

6c2）设定优化目标函数其中，h是一个变量，μ为一个调节参数，▽_NL为非局部梯度，使用增广拉格朗日模型求解上述优化问题，按如下步骤进行：

（1）设定l为迭代次数，令l＝0，l_max＝4，拉格朗日算子λ_m,n，λ_n,m的初始值为1， 中间变量f_m,n，f_n,m的初始值为1，m，n为图像的坐标值；

（2）用如下公式对上一次迭代图像x^k的每个像素点进行非局部总变差处理，得到处理后的图像的像素点

${\left(x_m^k\right)}^{l+1}=\frac{{\mathrm{\μx}}_m^{k\′}+2r\underset{n}{\mathrm{\Σ}}{w}_{m,n}{x}_n^{k\′}-\underset{n}{\mathrm{\Σ}}\sqrt{w_{m,n}}({\mathrm{\λ}}_{m,n}^l-{\mathrm{\λ}}_{n,m}^l)-r\underset{n}{\mathrm{\Σ}}\sqrt{w_{m,n}}(f_{m,n}^l-f_{n,m}^l)}{\mathrm{\μ}+2r\·\underset{n}{\mathrm{\Σ}{w}_{m,n}}}$

其中，为中间变量，

$x_m^{k\′}=\frac{\mathrm{\μ}{\left(x_m^k\right)}^l+2r\underset{n}{\mathrm{\Σ}}{w}_{m,n}{\left(x_n^k\right)}^l-\underset{n}{\mathrm{\Σ}}\sqrt{w_{m,n}}({\mathrm{\λ}}_{m,n}^l-{\mathrm{\λ}}_{n,m}^l)-r\underset{n}{\mathrm{\Σ}}\sqrt{w_{m,n}}(f_{m,n}^l-f_{n,m}^l)}{\mathrm{\μ}+2r\·\underset{n}{\mathrm{\Σ}}{w}_{m,n}},$

μ为一个调节参数，r是一个正常数，为拉格朗日算子，w_m,n为非局部权值，为中间变量，表示x^k的像素点，m，n为图像的坐标值；

（3）对中间变量以及拉格朗日算子的值进行更新，得到更新后的中间变量和拉格朗日算子

$f_{m,n}^{l+1}=\left(\begin{array}{c}\frac{1}{r}(1-\frac{1}{\left|{\mathrm{\ω}}_{m,n}\right|}){\mathrm{\ω}}_{m,n},\mathrm{if}\left|{\mathrm{\ω}}_{m,n}\right|>1\\ 0,\mathrm{otherwise}\end{array}\right.,$

${\mathrm{\λ}}_{m,n}^{l+1}={\mathrm{\λ}}_{m,n}^l+r(f_{m,n}^l-{\▿}_{\mathrm{NL}}{\left(x_m^k\right)}^{l+1})$

其中w_m,n为非局部权值，r是一个正常数，▽_NL为非局部梯度；

（4）如果l＞l_max，迭代终止， 为平滑处理后图像的像素点，m为图像的坐标值，否则迭代次数l加1，返回步骤（2）。