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

技术领域

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

背景技术

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

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

Jiao Wu等人在文献“Multivariate Compressive Sensing for Image Reconstruction in the  Wavelet 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行 的协方差矩阵

${\Sigma }_{s_i}={({\Sigma }^{-1}+{\Pi }^{-1}{\phi }_{s_i}^T{\phi }_{s_i})}^{-1}$

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

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

${\alpha }_{s_i}=\mathrm{diag}\left({\Sigma }_{s_i}\right)$

β＝diag(Π)

${\mu }_{s_i}={\alpha }_{s_i}(1./\beta ){\phi }_{s_i}^T(Y-\underset{k\ne i}{\overset{M}{\underset{k=1}{\Sigma }}}{\phi }_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\pi \right)}^Q{\left|{\Sigma }_{s_i}\right|}^{1/2}}}\exp \{-\frac{1}{2}{(x_{s_i}-{\mu }_{s_i})}^T{\Sigma }_{s_i}^{-1}(x_{s_i}-{\mu }_{s_i})\}$

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

$x_{s_i}=\mathrm{Gaussian}({\mu }_{s_i},{\Sigma }_{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通过以下的计算得到新的基础协方 差矩阵Σ和残差协方差矩阵Π：

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

$\Pi ={\left(\mathrm{diag}0\left(\mathrm{Gamma}(a_0+\frac{N}{2},b_0+\frac{1}{2}\mathrm{diag}\left({(Y-\mathrm{\Phi X})}^T(Y-\mathrm{\Phi 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算法重构效果好，可用于自然图像的重构。