CN102722866B - 基于主成份分析的压缩感知方法 - Google Patents

Abstract

本发明公开一种基于主成份分析的压缩感知方法，主要解决已有技术采样效率较低的问题。其实现步骤为：(1)从灰度自然图像库中取z幅图像，15≤z≤25，对取出的每幅图像延横竖两个方向每隔3个像素取一个32×32大小的子块，组成训练样本集x₁,x₂,...,x_m，m是训练样本的个数；对训练样本集x₁,x₂,...,x_m使用主成份分析的方法训练出满秩的观测矩阵Φ_f；(2)将需要被采样的图像分解成n个32×32大小的子块x₁,x₂,...,x_n，根据采样率s和满秩的观测矩阵Φ_f得到观测矩阵Φ，用观测矩阵Φ对每个图像子块进行采样，得到观测向量y；(3)根据观测向量y得到图像的初始解x⁰；(4)用x⁰进行迭代，直至满足终止条件,得到重构图像x'。本发明具有采样效率高，图像重构质量高，原理清晰，操作简单的优点，适用于自然图像的采样和重构。

Description

基于主成份分析的压缩感知方法

技术领域

本发明属于数字图像处理领域，特别涉及一种使用主成份分析的方法训练采样矩阵，可用于对自然图像进行采样和重构。

背景技术

当今社会是信息社会，人们对信息的追求不只停留在量上，更对信号获取的速度有所要求，即怎样又快又准确地获得需要的信息，而又能尽量地避免那些冗余信息的获取和处理。压缩感知理论是新兴的一种信号采样策略，它将信号的采样过程与压缩过程成功的融为了一体。该理论的前提就是在已知信号具有稀疏性或可压缩性的条件下，对数据进行采集、编解码。其核心主要是降低对信号的测量成本，在这个过程中会利用信号的特点，以便用更小的测量次数恢复这个信号；同时，在这其中包含了许多重要数学理论，涉及到数据采集，信息处理，以及模型优化等领域的知识。压缩感知理论框架主要包括三个方面：稀疏表示，观测矩阵的设计，以及重建优化算法。压缩感知理论开创了一个更经济有效的模拟信号数字化的途径。对于任意信号来说，只要能够找到与其稀疏表示相对应的一个稀疏表示空间，就可以应用压缩感知的理论进行采样与重构。

对于一个二维图像信号来说，由于其自身的高维度，使得压缩感知理论不得不面对这样一个“维数灾难”的问题。感知矩阵的存贮与计算将耗费大量的资源，这将对构造图像信号的实时采样系统带来极大的困难。另一方面，图像的重构过程也同样的面临极高的运算复杂度。因此，图像分块采样的方法在压缩感知领域被广泛的使用。在现有的理论中，由于构造简单，且满足有限等距条件，随机观测矩阵被广泛的使用。事实证明，针对于一个未知的信号，随机观测矩阵具有很好的采样效果，能够尽可能的采得信号的大部分信息。

但是实际上，信号采样所针对的并非是一些未知的信号，而是具有某些共性的一小部分信号，即自然图像块。在这种情况下，使用随机观测的方法，必然会造成一些采样资源的浪费，并且并不能很好的把握信号的特征，会丢失大量的信息。

此外，在信号重构的过程中，信号小波系数的零范数最小先验被广泛的使用。而实际上，小波变换具有弱稀疏性，也就是说信号x在小波基Ψ下的稀疏表示α的l₀范数并不是很小，从某种意义上来讲，α的l₀范数是一个很大的值，因为小波变换具有“高尖峰，长拖尾”的性质，所以α中虽然有很多比较小的系数，但真正为零的系数并不是很多。因此，使用零范数最小先验重构的稀疏系数通常带有很高的误差,许多小系数的丢弃会使图像的一些细节丢失，大大的影响图像的重构质量。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于主成份分析的压缩感知方法，以避免造成采样资源的浪费和图像的细节的丢失，更好的把握信号的特征，提高图像的重构效果。

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

(1)从灰度自然图像库中取z幅常见的灰度自然图像，15≤z≤25，对取出的每幅图像延横竖两个方向每隔3个像素取一个32×32大小的子块，组成训练样本集x₁,x₂,...,x_m，其中m是训练样本的个数，对训练样本集x₁,x₂,...,x_m使用主成份分析的方法训练出满秩的观测矩阵Φ_f；

(2)将需要被采样的图像分成n个32×32大小的图像子块x₁,x₂,...,x_n，n为图像子块的个数，给出信号的采样率s，根据采样率s，取满秩的观测矩阵Φ_f的前M行作为观测矩阵Φ，其中 表示向下取整，N＝1024为图像子块的维数，用观测矩阵Φ对每个图像子块进行采样，得到观测向量 $y=\left[\begin{array}{c}{y}_1\\ y_2\\ \·\\ \·\\ \·\\ y_n\end{array}\right],$ 其中y_i＝Φx_i，i＝1,2,...,n；

(3)根据观测向量y得到图像信号x的初始解：

$x^0={[x_1^0,x_2^0,....x_n^0]}^T$

其中i＝1,2，...,n，T表示矩阵的转置；

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

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

4b)对上一次迭代图像x^k进行三维块匹配BM3D滤波平滑处理，得到滤波平滑后的图像

4c)将滤波平滑后的图像分成n个32×32大小的子块对每个子块进行投影处理，将投影到超平面{e|Φ_ie＝y_i,e∈R^32×32}上，e是一个变量，R^32×32表示32×32维的信号，得到投影处理后的图像子块：i＝1,2,...,n，这些投影处理后的图像子块组成迭代后的图像

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

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

(1)本发明使用训练的感知矩阵来对信号进行采样，能够获取被感知信号的绝大部分能量，有效的利用采样资源，保留图像信号的大量信息；

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

(3)本发明成功的将能量最小初始解与标准差最小初始解统一起来，使用伪逆法就可以得到一个很好的初始解，用这个初始解可以很好的重构出原始图像。

附图说明

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

图2是用本发明训练出的变换空间示意图；

图3是用本发明训练出的满秩观测矩阵示意图；

图4是用本发明训练出的观测矩阵示意图；

图5是用本发明和BCS-SPL-DWT方法得到的重构图像的对比图。

具体实施方式

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

步骤1，训练出满秩的观测矩阵Φ_f。

为了寻求图像子块的共性，使用主成份分析的方法来训练出满秩的观测矩阵Φ_f，其步骤如下：

1a）从灰度自然图像库中取z幅常见的灰度自然图像，15≤z≤25，对取出的每幅图像延横竖两个方向每隔3个像素取一个32×32大小的子块，组成训练样本集x₁,x₂,...,x_m，其中m是训练样本的个数，本次实验中取z＝19，m＝4935；

1b)求解出训练样本集x₁,x₂,...,x_m的均值向量μ以及协方差矩阵m是训练样本的个数，T表示矩阵的转置；

1c)求解出协方差矩阵E的特征值λ_j，j＝0,1,...,r-1，r是特征值的个数，将它的特征值从大到小排列：λ₀≥λ₁≥...≥λ_r-1，其对应的特征向量为u_j，则满秩的观测矩阵Φ_f＝[u₀,u₁，...u_r-1]^T，T表示矩阵的转置，训练出的满秩观测矩阵Φ_f示意图如图3所示，满秩的观测矩阵即为变换空间矩阵，变换空间示意图如图2所示，其中图2（a）是信号空间，图2（b）是变换空间。

步骤2，给出信号的采样率s，对需要被采样的图像进行观测，得到观测向量y。

2a）将需要被采样的图像分解成n个32×32大小的图像子块x₁,x₂,...,x_n，n为图像子块的个数；

2b）根据采样率s，取Φ_f的前M行作为观测矩阵Φ，其中 表示向下取整，N＝1024为子块的维数，当M＝64，观测矩阵Φ的示意图如图4所示；

2c)用观测矩阵Φ对每个图像子块进行采样，得到观测向量 $y=\left[\begin{array}{c}{y}_1\\ y_2\\ \·\\ \·\\ \·\\ y_n\end{array}\right],$ 其中y_i＝Φx_i，i＝1,2,...,n。

步骤3，根据观测向量y得到图像信号x的初始解：

$x^0={[x_1^0,x_2^0,....x_n^0]}^T$

其中i＝1,2，...,n，T表示矩阵的转置。

步骤4，使用迭代方法对图像信号进行重构。

为了得到更好的重构效果，引入三维块匹配BM3D的思想来重构图像，其实现步骤如下：

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

4b)对上一次迭代图像x^k进行滤波平滑处理，得到滤波平滑后的图像滤波平滑的方法有多种方法，例如三维块匹配BM3D和非局部总变差NLTV，本实例中采用三维块匹配BM3D的方法，具体平滑过程参考文献“Image and video denoising bysparse 3D transform-domain collaborative filtering”，简单步骤如下：

4b1）把上一次迭代图像x^k分成一定大小的块，根据图像块之间的相似性，把具有相似结构的二维图像块组合在一起形成三维信号；

4b2）用联合滤波的方法对这些三维信号进行滤波处理，通过逆变换得到三维信号中所有图像块的估计值，把这些图像块的估计值返回到它们的在上一次迭代图像x^k中的原始位置，对得到的有重叠的局部块进行加权平均，得到滤波平滑后的图像

4c)将滤波平滑后的图像分成n个32×32大小的子块对每个子块进行投影处理，将投影到超平面{e|Φ_ie＝y_i,e∈R^32×32}上，e是一个变量，R^32×32表示32×32维的信号，得到投影处理后的图像子块：i＝1,2,...,n，这些投影处理后的图像子块组成迭代后的图像x^k+1；

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

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

1．仿真条件

从灰度自然图像库中取19幅常见的灰度自然图像，用来训练满秩的观测矩阵Φ_f；用本发明对4幅512×512大小的标准自然图像进行测试，它们分别为Lena图像，Barbara图像，Peppers图像和Goldhill图像。

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

仿真1、用本发明与现有BCS-SPL-DWT方法在采样率s=0.2的情况下，对两幅大小为512×512的标准自然图像Lena，Barbara进行重构，其重构结果如图5所示，其中图5（a）是用BCS-SPL-DWT方法对Lena图像的重构结果图，图5（b）是用本发明对Lena图像的重构结果图，图5（c）是用BCS-SPL-DWT方法对Barbara图像的重构结果图，图5（d）是用本发明对Barbara图像的重构结果图，从图5可以看出，本发明可以大幅度的提高重构图像的质量，无论在边缘区域还是在平滑区域重构效果都得到了大幅度的提升。

仿真2、用本发明与现有BCS-SPL-DWT方法分别在采样率s=0.2,0.3,0.4和0.5的情况下，对4幅大小为512×512的标准自然图像Lena，Barbara，Peppers和Goldhill进行重构，重构结果的峰值信噪比的对比结果如表1所示：

表1：峰值信噪比的比较结果（单位：db）

从表1可以看出，本发明可以大幅度的提高重构图像的峰值信噪比，表明重构图像的质量好。

综上，本发明提出了学习压缩感知的思想，打破了稀疏先验对压缩感知理论的统治地位，使用图像的一些空域先验与结构信息，如：标准差最小与三维块匹配BM3D最小，达到良好的重构效果。同时，本发明利用训练的方式来构造感知矩阵，这样的感知矩阵对于图像子块有着很好的采样效果，它能够采得图像信号的大多数的能量，达到很好的重构效果。

Claims (3)

1.一种基于主成份分析的压缩感知方法，包括如下步骤：

(1)从灰度自然图像库中取z幅常见的灰度自然图像，15≤z≤25，对取出的每幅图像沿横竖两个方向每隔3个像素取一个32×32大小的子块，组成训练样本集x₁,x₂,...,x_m，其中m是训练样本的个数，对训练样本集x₁,x₂,...,x_m使用主成份分析的方法训练出满秩的观测矩阵Φ_f；

(2)将需要被采样的图像分成n个32×32大小的图像子块x₁,x₂,...,x_n，n为图像子块的个数，给出信号的采样率s，根据采样率s，取满秩的观测矩阵Φ_f的前M行作为观测矩阵Φ，其中  表示向下取整，N＝1024为图像子块的维数，用观测矩阵Φ对每个图像子块进行采样，得到观测向量其中y_i＝Φx_i，i＝1,2,...,n；

(3)根据观测向量y得到图像信号x的初始解：

其中 i＝1,2，...,n，T表示矩阵的转置；

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

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

4b)对上一次迭代图像x^k进行三维块匹配BM3D滤波平滑处理，得到滤波平滑后的图像 

4c)将滤波平滑后的图像 分成n个32×32大小的子块 对每个子块进行投影处理，将 投影到超平面{e|Φ_ie＝y_i,e∈R^32×32}上，e是一个变量，R^32×32表 示32×32维的信号，得到投影处理后的图像子块： i＝1,2,...,n，这些投影处理后的图像子块组成迭代后的图像x^k+1；

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

2.根据权利要求1所述的基于主成份分析的压缩感知方法，其中步骤(1)所述的对训练样本集x₁,x₂,...,x_m使用主成份分析的方法训练出满秩的观测矩阵Φ_f，按如下步骤进行：

2a)求解出训练样本集x₁,x₂,...,x_m的均值向量μ以及协方差矩阵 m是训练样本的个数，T表示矩阵的转置；

2b)求解出协方差矩阵E的特征值λ_j，j＝0,1,...,r-1，r是特征值的个数，将它的特征值从大到小排列：λ₀≥λ₁≥...≥λ_r-1，其对应的特征向量为u_j，则满秩的观测矩阵Φ_f＝[u₀,u₁，...u_r-1]^T，T表示矩阵的转置。

3.根据权利要求1所述的基于主成份分析的压缩感知方法，其中步骤4b)所述的对上一次迭代图像x^k进行三维块匹配BM3D滤波平滑处理，按如下步骤进行：

4b1）将上一次迭代图像x^k分成一定大小的块，根据图像块之间的相似性，把具有相似结构的二维图像块组合在一起形成三维信号；

4b2）用联合滤波的方法对这些三维信号进行滤波处理，再通过逆变换得到三维信号中所有图像块的估计值，把这些图像块的估计值返回到它们的在上一次迭代图像x^k中的原始位置，对得到的有重叠的局部块进行加权平均，得到滤波平滑后的图像