CN104463927A - 基于多变量观测和相关性约束的压缩感知图像重构方法 - Google Patents

Abstract

本发明提供一种基于多变量高斯分布的压缩感知图像重构方法，重构过程：根据观测矩阵和测量矩阵的相关性，确定求解系数顺序的非零索引集；初始化基础协方差矩阵和残差协方差矩阵及系数矩阵；根据观测矩阵、观测向量、基础协方差和残差协方差矩阵得到索引集中对应的系数矩阵中行的均值向量和协方差矩阵，建立多变量高斯模型，求解该行系数，最终得到系数矩阵；根据保留的低频系数和重构的高频系数进行小波逆变换，得到重构图。本发明利用小波系数的统计特性和聚集性，建立多变量高斯模型，对模型参数自适应修正，有效提高了图像的重构质量。本发明重构效果好，可用于自然图像重构。

Description

基于多变量观测和相关性约束的压缩感知图像重构方法

技术领域

本发明属于图像处理技术领域，主要涉及统计压缩感知图像重构方法，具体是一种基于多变量高斯模型的压缩感知图像重构方法，可用于对自然图像进行重构。

背景技术

在图像处理技术领域，Nyquist曾提出奈奎斯特采样定理，该定理指出要从观测数据中精确重构信号，采样速率至少是信号带宽的两倍。近几年，出现了一种新的信号处理理论压缩感知(Compressed Sensing CS)，该理论在采样的同时实现压缩，大大降低了信号的采样频率，同时可以精确重构信号。压缩感知理论主要包括信号的稀疏表示、信号的观测和信号的重构等三个方面。在信号重构方面，通过求解零范数问题来重构图像。

Lihan He等人在文献“Exploiting Structure in Wavelet-Based Bayesian CompressiveSensing”中提出基于小波树结构Bayesian压缩感知图像重构方法。该方法对多尺度小波系数构造分层Bayesian模型，即单变量高斯分布模型，通过MCMC采样重构图像。该方法存在的不足是，将图像展开成列向量，由于单变量之间相互独立，破坏了图像的结构性和小波系数的聚集性，并且对计算机内存要求很高，限制了处理图像的大小。

Jiao Wu等人在文献“Multivariate Compressive Sensing for Image Reconstruction in theWavelet Domain:Using Scale Mixture Models”中提出基于混合尺度模型的多变量压缩感知图像重构(MPA)。该方法对小波系数构造多变量先验分布模型，使用多种不同的先验模型，对小波系数的统计相关性进行建模。该方法的不足是，虽然构造的是多变量模型，但在求解处理时认为变量之间是相互独立的，没有充分利用小波系数的聚集性。

综上，对于小波域下的单高斯压缩感知图像重构，其优点是：运算量小，耗时小，操作简单；其不足之处是：在小波域下，系数具有聚集性，展成列向量破坏了系数的聚集性，导致重构的图像质量不好。对于基于混合尺度模型的压缩感知图像重构方法，其有优点是节省内存，计算简单，缺点同样是破坏了小波系数的聚集性。

发明内容

本发明的目的是针对小波域下的单高斯分布压缩感知图像重构方法中，没有充分利用小波系数聚集性的缺点，提出一种多变量高斯分布的压缩感知图像重构方法，以便优化图像重构算法，提高图像重构质量。

实现本发明的技术方案是：基于多变量观测和相关性约束的压缩感知图像重构方法，包括如下具体步骤：

步骤1：接收方接收图像发送方发送的正交随机高斯观测矩阵Φ、低频小波系数矩阵C_N×N、三个高频子带的测量矩阵Y₁,Y₂,Y₃，由于三个高频子带的重构方法相同，统一用Y_M×N'表示，其中测量矩阵Y＝Φ*B，矩阵B_N'×Q有由矩阵A_N×N变换得到，A为小波分解得到的高频子带系数矩阵；

步骤2：计算观测矩阵Φ和测量矩阵Y的相关性矩阵U＝Φ'*Y＝(u_ij)_N'×Q，按行求和得到相关性向量u＝(u₁,u₂,…,u_N')^T，其中对u的各个分量排序，设置阈值c，c为非零行的行数，得到索引集S＝{s₁,s₂,…,s_i,…,s_c}，使得

步骤3：初始化多变量高斯模型的基础协方差矩阵Σ和残差协方差矩阵Π，初始化待重构的系数矩阵X＝(x_ij)_N'×Q＝(x₁,x₂,…,x_N')^T，x_i＝(x_i1,x_i2,…,x_iQ)(i＝1,2,…,N')为系数矩阵的第i行，N'为系数矩阵的行数，设置初始迭代次数n＝1；

步骤4：根据观测矩阵Φ、测量矩阵Y、整体协方差矩阵Σ、残差协方差矩阵Π和索引集得到系数矩阵X的第s_i行的行系数的均值向量和协方差矩阵建立对应的多变量高斯模型，生成行号不在索引集内的行系数为零，得到本次迭代的系数矩阵X＝(x₁,x₂,…,x_N')^T；

步骤5：根据观测矩阵Φ、测量矩阵Y以及迭代生成的系数矩阵X更新基础协方差矩阵Σ和残差协方差矩阵Π；

步骤6：根据迭代次数n判断迭代状态，如果满足终止条件，得到最终的系数矩阵；

步骤7：根据保留的低频子带系数C和迭代产生的小波系数矩阵X，进行小波逆变换，得到原图的重构图。

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

本发明利用图像小波系数的统计特性和聚集性，建立多变量高斯模型，并对模型的参数进行自适应修正，根据相关性确定非零索引集，提高了图像重构质量。本发明对图像数据进行一层小波变换后，得到一个低频子带系数矩阵和三个高频子带系数矩阵，低频子带包含图像的轮廓信息，三个高频子带包含图像的更多的细节信息，如光滑、纹理、边缘等信息，高斯分布能够很好地描述高频子带的统计特性，多变量分布可以描述高频子带的聚集性，因此，有效提高了图像的重构质量。

附图说明

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

图2是系数矩阵变换的示意图，A为原系数矩阵，B为多变量系数矩阵；

图3是多变量观测的示意图，Φ为观测矩阵，B为多变量系数矩阵，Y为观测矩阵；

图4(a)为测试图像，4(b)为测试图像的局部放大图；

图5(a)为本发明方法在采样率为30％时的重构图，图5(b)为对应的局部放大图，图5(c)为MPA算法在采样率为30％时的重构图，图5(d)为对应的局部放大图；

图6(a)为本发明方法在采样率为40％时的重构图，图6(b)为对应的局部放大图，图6(c)为MPA算法在采样率为40％时的重构图，图6(d)为对应的局部放大图；

具体实施方式

下面结合附图和实施例对本发明详细说明：

实施例1

本发明基于多变量观测和相关性约束的压缩感知图像重构方法，如图1所示，具体实施方式包括有如下步骤：

步骤1：接收方接收图像发送方发送的正交随机高斯观测矩阵Φ、低频小波系数矩阵C_N×N、三个高频子带的测量矩阵Y₁,Y₂,Y₃，由于三个高频子带的重构方法相同，统一用Y_M×N'表示，其中测量矩阵Y＝Φ*B，参见图3，矩阵B_N'×Q由矩阵A_N×N变换得到，参见图2，A为小波分解得到的高频子带系数矩阵；

步骤2：计算观测矩阵Φ和测量矩阵Y的相关性矩阵U＝Φ'*Y＝(u_ij)_N'×Q，按行求和得到相关性向量u＝(u₁,u₂,…,u_N')^T，其中对u的各个分量排序，设置阈值c，c根据实验测试得到，与图像的非光滑区域有关，本例中c为向量u中为零元素的个数，c＝1800，得到索引集S＝{s₁,s₂,…,s_i,…,s_c}，使得

步骤3：初始化多变量高斯模型的基础协方差矩阵Σ和残差协方差矩阵Π，初始化待重构的系数矩阵X＝(x_ij)_N'×Q＝(x₁,x₂,…,x_N')^T，x_i＝(x_i1,x_i2,…,x_iQ)(i＝1,2,…,N')为系数矩阵的第i行，N'为系数矩阵的行数，设置初始迭代次数n＝1。

步骤4：根据观测矩阵Φ、测量矩阵Y、整体协方差矩阵Σ、残差协方差矩阵Π和索引集得到系数矩阵X的第s_i行的行系数的均值向量和协方差矩阵建立对应的多变量高斯模型，生成行号不在索引集内的行系数为零，得到本次迭代的系数矩阵X＝(x₁,x₂,…,x_N')^T。建立多变量高斯模型并产生系数矩阵X的具体步骤包括：

4.1初始化行系数求解顺序i＝1，待求解行系数的行号为s_i＝s₁；

4.2根据基础协方差矩阵Σ、残差协方差矩阵Π，观测矩阵Φ得到系数矩阵X的第s_i行的协方差矩阵

${\mathrm{\Σ}}_{s_i}={({\mathrm{\Σ}}^{-1}+{\mathrm{\Π}}^{-1}{\mathrm{\φ}}_{s_i}^T{\mathrm{\φ}}_{s_i})}^{-1}$

其中(·)^-1表示矩阵的逆，为观测矩阵Φ的第s_i列，为的转置；

4.3根据残差协方差矩阵Π和系数矩阵X的第s_i行的协方差矩阵得到该行的均值向量

${\mathrm{\α}}_{s_i}=\mathrm{diag}\left({\mathrm{\Σ}}_{s_i}\right)$

β＝diag(Π)

${\mathrm{\μ}}_{s_i}={\mathrm{\α}}_{s_i}(1./\mathrm{\β}){\mathrm{\φ}}_{s_i}^T(Y-\underset{k\≠i}{\overset{M}{\underset{k=1}{\mathrm{\Σ}}}}{\mathrm{\φ}}_k{x}_k)$

其中diag(·)表示矩阵对角线元素组成的向量，1./β为向量β的每个元素分别取倒数组成的向量，φ_i为观测矩阵Φ的第i列，x_k为系数矩阵X的第k行；

4.4根据系数矩阵X的第s_i行的均值向量和协方差矩阵建立对应的多变量高斯模型：

$p\left(x_{s_i}\right)=\frac{1}{\sqrt{{\left(2\mathrm{\π}\right)}^Q{\left|{\mathrm{\Σ}}_{s_i}\right|}^{1/2}}}\exp \{-\frac{1}{2}{(x_{s_i}-{\mathrm{\μ}}_{s_i})}^T{\mathrm{\Σ}}_{s_i}^{-1}(x_{s_i}-{\mathrm{\μ}}_{s_i})\}$

4.5根据高斯模型，生成系数矩阵X第s_i行系数

$x_{s_i}=\mathrm{Gaussian}({\mathrm{\μ}}_{s_i},{\mathrm{\Σ}}_{s_i})$

其中，表示生成一个服从均值向量为协方差矩阵的多变量高斯分布的向量；

4.6如果行系数求解顺序i＜c，则i＝i+1，行系数行号s_i＝s_i+1，返回4.2，参见图1，否则，系数矩阵X中行号不在索引集内的行系数为零，得到系数矩阵X＝(x₁,x₂,…,x_N')^T。

本发明将小波系数矩阵变换为多变量矩阵，建立多变量高斯模型，充分体现了图像小波系数的聚集性，准确刻画了小波系数，为图像的精确重构打下了基础。

步骤5：根据观测矩阵Φ、测量矩阵Y以及迭代生成的系数矩阵X更新基础协方差矩阵Σ和残差协方差矩阵Π。

步骤6：根据迭代次数n判断迭代状态，如果满足终止条件，得到最终的系数矩阵，判断迭代状态的具体步骤包括：

6.1若迭代次数n≤N₁，则n＝n+1，若没有达到最大迭代次数，转到步骤4，参见图1；

6.2若N₁＜n≤N₁+N₂，累计迭代结果S＝S+X，迭代次n＝n+1，若没有达到最大迭代次数，转到步骤4，参见图1；

6.3若n＝N₁+N₂，累计迭代结果S＝S+X，当达到最大迭代次数时，停止迭代，系数矩阵X＝S/N₂，输出系数矩阵X。

步骤7：根据保留的低频子带系数C和迭代产生的小波系数矩阵X，进行小波逆变换，得到原图的重构图，参见图4(g)。

本发明利用图像小波系数的统计特性和聚集性，建立多变量高斯模型，并对模型的参数进行自适应修正，根据相关性确定索引集，提高了图像重构质量。本发明对图像数据进行一层小波变换后，得到一个低频子带系数矩阵和三个高频子带系数矩阵，高斯分布能够很好地描述高频子带的统计特性，多变量分布可以描述高频子带的聚集性，因此，采用本发明对图像重构，有效提高了图像的重构质量。

实施例2

基于多变量观测和相关性约束的压缩感知图像重构方法同实施例1，其中步骤5中根据迭代产生的系数矩阵更新基础协方差矩阵Σ和残差协方差矩阵Π，具体步骤如下：

5.1已知系数矩阵X以及观测矩阵Φ和测量矩阵Y，基础协方差矩阵Σ和残差协方差矩阵Π分别均服从伽马分布。

5.2初始化伽马分布的参数a₀，b₀，c₀，d₀为给定1×Q常向量，其中Q为系数矩阵X行向量的维数，伽马分布参数的每个向量的每个元素的值均为0.000001。

5.3根据系数矩阵X以及观测矩阵Φ和测量矩阵Y通过以下的计算得到新的基础协方差矩阵Σ和残差协方差矩阵Π：

$\mathrm{\Σ}={\left(\mathrm{diag}0\left(\mathrm{Gamma}(c_0+N^{\′},d_0+\frac{1}{2}\mathrm{diag}\left(X^{T}X\right))\right)\right)}^{-1}$

$\mathrm{\Π}={\left(\mathrm{diag}0\left(\mathrm{Gamma}(a_0+\frac{N}{2},b_0+\frac{1}{2}\mathrm{diag}\left({(Y-\mathrm{\ΦX})}^T(Y-\mathrm{\ΦX})\right))\right)\right)}^{-1}$

其中diag0(c)为方阵，方阵的对角线元素为c，非对角线元素为0，Gamma(a,b)表示生成一个服从形状参数向量为a，尺度参数向量为b的伽马分布向量，其中向量a和b的维数相同，伽马分布产生的向量维数与a相同。本发明利用系数矩阵自适应的调整多变量高斯模型的参数，使得模型更准确的描述图像的小波系数，进而提高了图像的重构精度。

实施例3

基于多变量观测和相关性约束的压缩感知图像重构方法同实施例1-2。

本发明对图像的重构效果可以通过用本发明方法进行重构的图像和MPA算法获得重构图像对比说明。

图5和图6是本发明方法与MPA算法的实验结果对比图。在采样率均为30％时，由图5(a)、图5(c)和图4(a)对比可以看出，本发明方法的重构图即图5(a)的光滑区域与原图即图4(a)的光滑区域十分相似，而MPA算法重构图即图5(c)有明显的斑点，并且本发明的重构图边缘更为清晰。

图6(a)和图6(c)均为采样率为40％时的重构图，图5(a)为本发明的重构图，与MPA算法的重构图对比，其结果仍是本发明的重构图像更加清晰，细节更丰富。从局部放大图的对比，也可以明显看出来本发明的重构图即图6(b)比同样采样率下的MPA算法的的重构图即图6(d)边缘清晰，对比度和表现力都处于优势。

实施例4

基于多变量观测和相关性约束的压缩感知图像重构方法同实施例1-2。

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

1.仿真条件：

本发明的仿真在windows7，，SPI，CPU Intel(R)Core(TM)2，基本频率3.00GHz，软件平台为Matlab R2012a上运行，仿真选用的是512×512的Lena图像。

2.仿真内容与结果：

本仿真中，使用MPA算法和本发明采用的多变量高低模型对图像进行重构，采样率分别为为30％，35％和40％，对比图如图4。重构结果的峰值信噪比PSNR的平均值如表1所示。

表1 Lena图像在不同采样率下用MPA算法和本发明方法的重构结果

简而言之，本发明的基于多变量高斯分布的压缩感知图像重构方法，解决了充分利用小波系数的相关特性和聚集性的技术问题。本发明的重构过程为：根据观测矩阵和测量矩阵的相关性，确定求解系数顺序的索引集；初始化基础协方差矩阵和残差协方差矩阵；根据观测矩阵、观测矩阵、基础协方差和残差协方差矩阵得到索引集中对应的系数矩阵中行的均值向量和协方差矩阵；根据该行的均值向量和协方差矩阵，建立多变量高斯模型，求解该行系数，得到系数矩阵；根据保留的低频系数和重构的高频系数进行小波逆变换，得到重构图。本发明利用图像小波系数的统计特性和聚集性，建立多变量高斯模型，并对参数自适应修正，根据相关性确定非零索引集，有效提高了图像的重构质量。试验表明，本发明方法比MPA算法重构效果好，可用于自然图像的重构。

Claims (3)

1.基于多变量观测和相关性约束的压缩感知图像重构方法，其特征在于，包括有如下具体步骤：

步骤1：接收方接收图像发送方发送的正交随机高斯观测矩阵Φ、低频小波系数矩阵C_N×N、三个高频子带的测量矩阵Y₁,Y₂,Y₃，由于三个高频子带的重构方法相同，统一用Y_M×N'表示，其中测量矩阵Y＝Φ*B，矩阵B_N'×Q有由矩阵A_N×N变换得到，A为小波分解得到的高频子带系数矩阵；

步骤2：计算观测矩阵Φ和测量矩阵Y的相关性矩阵U＝Φ'*Y＝(u_ij)_N'×Q，按行求和得到相关性向量u＝(u₁,u₂,…,u_N')^T，其中对u的各个分量排序，设置阈值c，c为非零行的行数，得到索引集S＝{s₁,s₂,…,s_i,…,s_c}，使得 $\left|u_{s_1}\right|>\left|u_{s_2}\right|>...>\left|u_{s_c}\right|;$

步骤3：初始化多变量高斯模型的基础协方差矩阵Σ和残差协方差矩阵Π，初始化待重构的系数矩阵X＝(x_ij)_N'×Q＝(x₁,x₂,…,x_N')^T，x_i＝(x_i1,x_i2,…，x_iQ)(i＝1,2,…,N')为系数矩阵的第i行，N'为系数矩阵的行数，设置初始迭代次数n＝1；

步骤4：根据观测矩阵Φ、测量矩阵Y、整体协方差矩阵Σ、残差协方差矩阵Π和索引集得到系数矩阵X的第s_i行的行系数的均值向量和协方差矩阵建立对应的多变量高斯模型，生成行号不在索引集内的行系数为零，得到本次迭代的系数矩阵X＝(x₁,x₂,…,x_N')^T；

步骤5：根据观测矩阵Φ、测量矩阵Y以及迭代生成的系数矩阵X更新基础协方差矩阵Σ和残差协方差矩阵Π；

步骤6：根据迭代次数n判断迭代状态，如果满足终止条件，得到最终的系数矩阵；

步骤7：根据保留的低频子带系数C和迭代产生的小波系数矩阵X，进行小波逆变换，得到原图的重构图。

2.根据权利要求1所述的基于多变量高斯分布的压缩感知图像重构方法，其特征在于，所述步骤4中计算索引集对应的系数矩阵X第si行的行系数的均值向量和协方差矩阵具体步骤如下：

4.1初始化行系数求解顺序i＝1，待求解行系数的行号为s_i＝s₁；

4.2根据基础协方差矩阵Σ、残差协方差矩阵Π，观测矩阵Φ得到系数矩阵X的第s_i行的协方差矩阵

${\mathrm{\Σ}}_{s_i}={({\mathrm{\Σ}}^{-1}+{\mathrm{\Π}}^{-1}{\mathrm{\φ}}_{s_i}^T{\mathrm{\φ}}_{s_i})}^{-1}$

其中(·)^-1表示矩阵的逆，为观测矩阵Φ的第s_i列，为的转置；

4.3根据残差协方差矩阵Π和系数矩阵X的第s_i行的协方差矩阵得到该行的均值向量

${\mathrm{\α}}_{s_i}=\mathrm{diag}\left({\mathrm{\Σ}}_{s_i}\right)$

β＝diag(Π)

${\mathrm{\μ}}_{s_i}={\mathrm{\α}}_{s_i}(1./\mathrm{\β}){\mathrm{\φ}}_{s_i}^T(Y-\underset{k\≠i}{\overset{M}{\underset{k=1}{\mathrm{\Σ}}}}{\mathrm{\φ}}_k{x}_k)$

其中diag(·)表示矩阵对角线元素组成的向量，1./β为向量β的每个元素分别取倒数组成的向量，φ_i为观测矩阵Φ的第i列，x_k为系数矩阵X的第k行；

4.4根据系数矩阵X的第s_i行的均值向量和协方差矩阵建立对应的多变量高斯模型：

$p\left(x_{s_i}\right)=\frac{1}{\sqrt{{\left(2\mathrm{\π}\right)}^Q}{\left|{\mathrm{\Σ}}_{s_i}\right|}^{1/2}}\exp \{-\frac{1}{2}{(x_{s_i}-{\mathrm{\μ}}_{s_i})}^T{\mathrm{\Σ}}_{s_i}^{-1}(x_{s_i}-{\mathrm{\μ}}_{s_i})\}$

4.5根据高斯模型，生成系数矩阵X第s_i行系数

$x_{s_i}=\mathrm{Gaussian}({\mathrm{\μ}}_{s_i},{\mathrm{\Σ}}_{s_i})$

其中，表示生成一个服从均值向量为协方差矩阵的多变量高斯分布的向量；

4.6如果行系数求解顺序i＜c，则i＝i+1，行系数行号s_i＝s_i+1，返回4.2，否则，行号不在索引集内的行系数为零，得到系数矩阵X＝(x₁,x₂,…,x_N')^T。

3.根据权利要求1所述的基于多变量高斯分布的压缩感知图像重构方法，所述步骤5中的根据迭代产生的系数矩阵更新基础协方差矩阵Σ和残差协方差矩阵Π，具体步骤如下：

5.1已知系数矩阵X以及观测矩阵Φ和测量矩阵Y，基础协方差矩阵Σ和残差协方差矩阵Π分别服从伽马分布；

5.2初始化伽马分布的参数a₀，b₀，c₀，d₀为给定1×Q常向量，每个向量的每个元素的值均为0.000001；

5.3系数矩阵X以及观测矩阵Φ和测量矩阵Y得到新的基础协方差矩阵Σ和残差协方差矩阵Π：

$\mathrm{\Σ}={\left(\mathrm{diag}0\left(\mathrm{Gamma}(c_0+N^{\′},d_0+\frac{1}{2}\mathrm{diag}\left(X^{T}X\right))\right)\right)}^{-1}$

$\mathrm{\Π}={\left(\mathrm{diag}0\left(\mathrm{Gamma}(a_0+\frac{N}{2},b_0+\frac{1}{2}\mathrm{diag}\left({(Y-\mathrm{\ΦX})}^T(Y-\mathrm{\ΦX})\right))\right)\right)}^{-1}$

其中diag0(c)为方阵，方阵的对角线元素为c，非对角线元素为0，Gamma(a,b)表示生成一个服从形状参数向量为a，尺度参数向量为b的伽马分布向量，其中向量a和b的维数相同，伽马分布产生的向量维数与a相同。