CN103985145A - 一种基于联合稀疏和先验约束的压缩感知图像重构方法 - Google Patents

Abstract

本发明公开了一种基于联合稀疏和先验约束的压缩感知图像重构方法，过程为：接收低频信息和高频子带分块观测，根据基于边缘信息的先验模型生成各块观测的位置块，并据此对各块观测进行边缘块观测和非边缘块观测的划分；对各非边缘块的观测执行局部聚类操作，对各类的聚类中心块观测，在定义了联合稀疏的适应度函数下通过遗传算法求解最优系数，并将各聚类中心块对应的最优系数作为同类各块观测的最优系数；对各边缘块观测，也使用遗传算法求解对应的最优系数；最后合并所有块的最优系数，并结合低频信息进行小波逆变换获得重构图像。与OMP、BP及IHT方法相比，本发明较好地利用了图像的结构信息，获得了质量较好的重构图像。

Description

一种基于联合稀疏和先验约束的压缩感知图像重构方法

技术领域

本本发明属于图像处理技术领域，更进一步涉及图像重构，具体是一种基于联合稀疏和先验约束的压缩感知图像重构方法。

背景技术

随着信息技术的迅速发展，图像处理技术被越来越多的应用到人们的生产和生活中。如利用卫星所获取的图像进行资源调查、灾害监测、城市规划，利用医学图像进行疾病的检测，利用工业图像对零件进行分类和质量检测等。

奈奎斯特定理指出，只有当采样速率达到信号带宽的两倍以上时，采集到的数字信号才能完整地保留原始信号中的信息，而在现实世界中，由于实际应用中图像数据量巨大，为了降低对信息存储、处理和传输的成本，人们对信号进行高速奈奎斯特采样和压缩编码后再进行存储和传输。但是这种处理方式造成了严重的资源浪费，针对这一问题，D.L.Donoho、CandèsE.J.等人提出了一种新的数据采集技术——压缩感知。压缩感知技术是利用信号的稀疏性，在远小于奈奎斯特采样速率的条件下进行采样，然后通过非线性重构算法准确地重构信号，这样大大降低了设备存储限制和计算的复杂度。目前压缩感知已成为学术界研究的热点，并不断被应用在图像处理领域和无线传感领域中。压缩感知理论主要包括信号的稀疏表示、信号的观测和信号的重构等三个方面。其中，信号重构是压缩感知技术的关键和核心。

在压缩感知技术中，图像信号重构的过程是对数字化信号进行处理的过程，这个过程离不开求解欠定方程组问题。E.Candes等人证明了，如果信号是稀疏或者可压缩的，求解欠定方程组的问题可以转化为最小化l₀范数问题，从而可重构信号。压缩感知重构的本源问题是l₀范数下非凸优化问题，该问题是非确定性多项式问题即NP难问题。目前直接求解l₀范数问题的方法是以正交匹配追踪OMP算法为代表的贪婪算法和以迭代阈值收缩IHT为代表的门限算法两类。

OMP算法是在每次迭代的过程中，基于贪婪的思想并通过局部优化的手段选择最能匹配信号结构的一个原子，并经过一系列逐步递增的方法构建信号的稀疏逼近。OMP算法每一次的迭代主要有两个步骤：原子选择和残差更新。OMP算法通过Gram-Schmidt正交化方法对已选择原子集合进行正交化处理，这样一来每次迭代所选取的最匹配原子均满足一定的条件，残差部分随着迭代次数的增大而迅速减少，因此用少量原子的线性组合来重构原始图像信号，从而有效地减少了迭代次数。但是OMP算法不能对所有图像信号都实现精确重构，重构结果不是很精确，算法也不具有鲁棒性。

IHT算法也是基于l₀范数的重构方法，直接关注稀疏信号的非零元素的个数，寻找最能逼近稀疏信号的K项支撑，迭代过程如下式所示：

xⁿ⁺¹＝H_K(xⁿ+Φ^T(y-Φxⁿ))

其中，xⁿ⁺¹是第n+1次迭代时的重构信号，H_K(θ)是一个非线性算子，它的功能是保留矢量θ中幅度最大的前K个元素，将其他元素都置为零，xⁿ是第n次迭代时的重构信号，Φ是高斯随机观测矩阵，y是观测向量。IHT算法的缺点对测量矩阵的过分依赖，计算复杂度高，运算时间长，阈值的大小对图像信号的重构结果影响较大。

以上两种算法有一个共同的缺点，那就是不能保证收敛到全局最优，造成图像的重构结果不够精确。因此，基于最小化l₀范数的非凸压缩感知重构方法还需要进一步地探索和研究。

西安电子科技大学的专利申请“基于滤波算子的交替优化压缩感知图像重构方法”（公开号：CN102568017A，申请号：201210001645.8，申请日：2012年1月4日）中公开了一种基于滤波算子的交替优化压缩感知图像重构方法,该方法利用稀疏系数位置的先验信息来指导求解稀疏系数的l₀范数，并把滤波和凸集投影作为进化算子引入到进化重构框架中，重构的图像纹理和边缘清晰。但效果的提升很大一部分是由于滤波和凸集操作起作用的，并未充分考虑图像本身的结构信息，因为我们希望通过尽可能地利用图像的结构信息来指导重构过程以获得较好的重构效果。

发明内容

本发明的目的在于针对现有压缩感知重构技术中在重构图像的过程中，并未考虑图像本身的结构特征，只是单纯地从数学角度出发进行图像重构算法的设计，导致重构效果不佳的问题，设计了更加合理的交叉变异算子及基于局部聚类块间相互约束的选择策略，提出了一种基于联合稀疏和先验约束的压缩感知重构方法，提高了重构图像的质量。

为实现上述目的，本发明采取以下技术方案：一种基于联合稀疏和先验约束的压缩感知重构方法，包括以下步骤：

(1)输入低频子带Y₀和三个高频子带的分块观测Y_t，其中t∈{v,h,d}，表示子带的方向，其中v表示垂直方向，h表示为水平方向，d表示对角方向；

(2)获得三个高频子带对应的位置矩阵P_t；

(2a)将三个高频子带系数置0，结合低频子带Y₀，做小波逆变换，得到一幅边缘模糊的图像I1；

(2b)用canny算子对边缘模糊的图像I1进行边缘检测，得到一幅只含有边缘信息的图像I2；

(2c)对只含有边缘信息的图像I2执行一层小波变换，得到一个低频子带Y₁和三个含有边缘信息的高频子带为S_t；

(2d)以作为方向为t的子带的阈值,将对应子带S_t中模值大于对应阈值的位置标记为1，意为该位置处为大系数，而将模值小于对应阈值的位置标记为0，意为该位置处为小系数，这样获得三个高频子带对应的位置矩阵P_t，其中μ_t为子带S_t的模值的均值,为伸缩因子，本文中取

(3)对位置矩阵P_t执行提取疫苗和注射疫苗的操作，得到位置矩阵P′_t；

(4)对位置矩阵P′_t执行分块操作，得到三个高频子带下各块观测对应的位置块p′_t,i；

为了记录初始位置块p′_t,i的信息，将p′_t,i保存为变量p″_t,i，即p″_t,i=p′_t,i，后续只对变量p″_t,i进行操作；

其中i表示块号，若图像大小为512×512，则i＝1,2,...256，若图像大小为256×256，则i＝1,2,...128

(5)根据得到的三个子带下各块观测的位置块p″_t,i对各块观测进行边缘块观测和非边缘块观测的划分；

(6)根据第（5）步的边缘块观测和非边缘块观测的划分结果，对三个高频子带的分块观测Y_t中所有非边缘块观测执行局部相似聚类，得到对应于三个子带中非边缘块观测聚类的集合：

$A=\left\{\begin{array}{cc}{a}_{t,1},a_{t,2},L& a_{t,c_t}\end{array}\right\},$

其中，a_t,i表示方向为t的子带的第i类对应的集合，其中i＝1,2,Lc_t，c_t表示方向为t的子带聚类的类别数；

(7)按照种群初始化策略分别初始化三个高频子带下各个聚类中心块观测及各个边缘块观测对应的系数块种群Q＝{q_t,i,j}，其中i表示各聚类中心块观测及各个边缘块观测的块号，j＝1,2,Ln，n为种群规模；

(8)对三个高频子带下的各个系数块种群Q＝{q_t,i,j}执行交叉操作，得到交叉后的子代系数块种群Q′＝{q′_t,i,j}；

(9)对三个高频子带下的各个子代系数块种群Q′执行变异操作，得到变异后的子代系数块种群Q″＝{q″_t,i,j}；

(10)对各聚类中心的系数块和边缘的系数块分别定义两种不同的适应度函数以对子代种群Q″执行相应的种群更新操作；

(11)分别从三个子带下各个子代种群Q″中选择出各个系数块对应的最优系数个体，若该系数块为非边缘块，则将其对应的最优系数个体作为其同类的各个系数块的最优系数块，然后对所有最优系数块执行合并块操作，形成各子带对应的系数B_t；

(12)若进化代数满足停止条件，则转步骤（13），否则，利用各个聚类中心系数块和边缘的系数块选择出的最优系数个体更新该块对应的位置块和系数块种群，转入步骤(8)；

(13)结合保留的低频子带系数Y₀及步骤(11)中得到的三个高频子带系数B_t，进行小波逆变换，得到重构图像。

所述步骤(5)中根据得到的三个子带下各块观测的位置块p″_t,i对各块观测进行边缘块观测和非边缘块观测的划分，若某块观测对应的位置块p″_t,i为全0矩阵，则将该块观测划作非边缘块观测，否则，将其作为边缘块观测。

所述步骤(6)中对三个高频子带的分块观测Y_t中所有非边缘块对应的观测执行局部相似聚类，t∈{v,h,d},以对大小为512×512的图像进行操作为例即块号i＝1,2,...256，具体过程如下：

(6.1)计算Y_t中各个非边缘块观测y_t,i的标准差σ_t,i，i表示块号；

(6.2)对各个非边缘块观测都初始化一个未聚类的标志mark_i，即令mark_i＝0，i表示块号；如果方向子带t下的第i个块观测为某一类的聚类中心块观测的话，则我们用符号a_t,i来表示这个类；初始时令i＝1，表示从第一块开始执行聚类操作；

(6.3)若第i个块观测为非边缘块观测且该块未被聚类，即mark_i＝0，则转(6.4)，否则，转(6.6)

(6.4)将第i个块观测作为类a_t,i的聚类中心块观测；

(6.5)观察第i块的八个邻域块中的所有块观测：假设j表示第i块的八个邻域块中的其中一个块观测，若第j个块观测为非边缘块观测且mark_j＝0，表示第j个块观测未被聚类，则计算第j个块观测的标准差σ_t,j与第i个块观测的标准差σ_t,i的差值，即C_j＝σ_t,i-σ_t,j，若|C_j|≤τ，其中τ为阈值，τ=0.01，则将第j个块观测加入聚类中心为第i个块观测的类a_t,i中，令mark_j＝1，表示该块已被聚类；

(6.6)令i＝i+1,若i≤256，则转(6.3),否则，表示所有块都已聚类完毕，统计该方向上聚类的类别数，记作c_t；

所述步骤(10)中对各聚类中心的系数块和边缘的系数块分别定义两种不同的适应度函数以对子代种群Q″执行相应的种群更新操作按如下过程进行：

(10.1)若当前块为非边缘块，也即当前块为一个聚类中心块，则按以下适应度函数计算各个个体的适应度f(q″_t,i,j)，表示如下：

$f\left(q_{t,i,j}^{\′\′}\right)=\frac{1}{\left|a_{t,i}\right|}\underset{k\∈a_{t,i}}{\mathrm{\Σ}}\frac{1}{{\left|\right|y_{t,k}-\mathrm{\Φ}{q}_{t,i,j}^{\′\′}\left|\right|}_2^2},\mathrm{t\∈}\{v,h,d\}$

其中，t∈{v,h,d}，表示的是该块所属子带的方向，其中v表示垂直方向，h表示为水平方向，d表示对角方向；由步骤（6）可得，a_t,i表示方向为t的子带中以第i块作为聚类中心块的类；y_t,k为方向为t的子带中类a_t,i中的第k块对应的观测，q″_t,i,j为该聚类中心块所对应的种群Q″中的第j个系数个体，f(q″_t,i,j)即为个体q″_t,i,j的适应度；

若当前块为边缘块，则按以下适应度函数计算各个个体的适应度f(q″_t,i,j)，表示如下：

$f\left(q_{t,i,j}^{\′\′}\right)=\frac{1}{{\left|\right|y_{t,i}-\mathrm{\Φ}{q}_{t,i,j}^{\′\′}\left|\right|}_2^2},t\∈\{v,h,d\}$

其中，t∈{v,h,d}，表示的是该块所属子带的方向，其中v表示垂直方向，h表示为水平方向，d表示对角方向；y_t,i为方向为t的子带中第i块对应的观测，q″_t,i,j为该块所对应的种群Q″中的第j个系数个体，f(q″_t,i,j)即为个体q″_t,i,j的适应度；

(10.2)将Q″中所有个体的适应度与上一代种群Q中所有个体的适应度相比较，从中选择前n个适应度较大的个体更新种群Q，并将适应度最大的个体作为本代进化的最优个体。

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

第一，本发明的重构方法是求解小波域下三个高频子带系数的l₀范数，克服了现在压缩感知框架中强加的有限等距性的约束条件，从而扩展了压缩感知的应用范围；

第二，本发明在重构图像时，构造了基于分块的先验模型，并根据所得的位置块对图像的各块进行了边缘块和非边缘块的划分；

第三,本发明在重构图像时对非边缘块进行了局部聚类操作，在遗传框架下只对各类的聚类中心块进行学习得到其最优系数，其中选择操作需同类所有块的联合约束，利用了局部图像块间相互约束的结构信息，并将各聚类中心块对应的最优系数作为同类各块的最优系数，对各边缘块也在遗传框架下通过进化学习得到对应的最优系数，这种方法克服了现有压缩感知技术中没有关注结构信息的缺点，获得了较好的效果。

附图说明

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

图2(a)是Barbara原始图像；

图2(b)是图2(a)的局部放大图；

图2(c)是本发明得到的Barbara重构图；

图2(d)是图2(c)的局部放大图；

图2(e)是OMP得到的Barbara的重构图；

图2(f)是图2(e)的局部放大图；

图2(g)是BP得到的Barbara重构图；

图2(h)是图2(g)的局部放大图；

图2(i)是IHT得到的Barbara重构图；

图2(j)是图2(i)的局部放大图；

图3是用本发明与现有技术重构出来的Barbara图像的峰值信噪比PSNR随采样率变化的趋势图；

图4(a)是Lena原始图像；

图4(b)是图4(a)的局部放大图；

图4(c)是本发明得到的Lena重构图；

图4(d)是图4(c)的局部放大图；

图4(e)是OMP得到的Lena的重构图；

图4(f)是图4(e)的局部放大图；

图4(g)是BP得到的Lena重构图；

图4(h)是图4(g)的局部放大图；

图4(i)是IHT得到的Lena重构图；

图4(j)是图4(i)的局部放大图；

图5是用本发明与现有技术重构出来的Lena图像的峰值信噪比PSNR随采样率变化的趋势图。

具体实施方式

下面结合附图和实施例对本发明的进行详细的描述。

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

步骤一，输入低频子带Y₀和三个高频子带的分块观测Y_t，其中t∈{v,h,d}，表示子带的方向，其中v表示垂直方向，h表示为水平方向，d表示对角方向；

步骤二，获得三个高频子带对应的位置矩阵P_t；

(2a)将三个高频子带系数置0，结合低频子带Y₀，做小波逆变换，得到一幅边缘模糊的图像I1；

(2b)用canny算子对边缘模糊的图像I1进行边缘检测，得到一幅只含有边缘信息的图像I2；

(2c)对只含有边缘信息的图像I2执行一层小波变换，得到一个低频子带Y₁和三个含有边缘信息的高频子带为S_t；

(2d)以作为方向为t的子带的阈值,将对应子带S_t中模值大于对应阈值的位置标记为1，意为该位置处为大系数，而将模值小于对应阈值的位置标记为0，意为该位置处为小系数，这样获得三个高频子带对应的位置矩阵P_t，其中μ_t为子带S_t的模值的均值,为伸缩因子，本文中取

步骤三，对位置矩阵P_t执行提取疫苗和注射疫苗的操作，得到位置矩阵P′_t；；

本步骤的具体实现如下（以垂直方向为例）：

首先，执行提取疫苗的操作，用大小为3×3的窗口在位置矩阵P_v上滑动，窗口的中心位置需保证滑过位置矩阵P_v的每一个位置，考虑该中心位的上、下、左、右四个邻域的值，如果其中有不少于3个邻域的值为1，则将该中心位置的疫苗取为1；如果4个邻域的值均为0，则将该中心位置的疫苗取为0；如果以上两种情况都不满足，则该中心位置的疫苗取其自身的值为疫苗，所有的疫苗值组成一个新的矩阵C_v，称作对应的疫苗矩阵。

然后，执行注射疫苗的操作，比较位置矩阵P_v和疫苗矩阵C_v的每一个元素，如对应位置处的元素相同，则不进行操作，否则，用疫苗矩阵C_v中的值代替位置矩阵中对应位置处的值，这样就得到对应于垂直方向的整个子带的位置矩阵P′_v；

对应于水平方向和对角方向的位置矩阵P_h，P_d也用上述步骤求得其对应的注射疫苗后的位置矩阵P′_h，P′_d。

步骤四，对位置矩阵P′_t执行分块操作，得到三个高频子带下各块观测对应的位置块p′_t,i；

为了记录初始位置块p′_t,i的信息，将p′_t,i保存为变量p″_t,i，即p″_t,i=p′_t,i，后续只对变量p″_t,i进行操作；

其中i表示块号；

步骤五，根据得到的三个子带下各块观测的位置块p″_t,i对各块观测进行边缘块观测和非边缘块观测的划分，若某块观测对应的位置块p″_t,i为全0矩阵，则将该块观测划作非边缘块观测，否则，将其作为边缘块观测；

步骤六，根据第（5）步的边缘块观测和非边缘块观测的划分结果，对三个高频子带的分块观测Y_t中所有非边缘块观测执行局部相似聚类，得到对应于三个子带中非边缘块观测聚类的集合：

$A=\left\{\begin{array}{cc}{a}_{t,1},a_{t,2},L& a_{t,c_t}\end{array}\right\},$

其中，a_t,i表示方向为t的子带的第i类对应的集合，其中i＝1,2,Lc_t，c_t表示方向为t的子带聚类的类别数；

以对大小为512×512的图像进行操作为例，即块号i＝1,2,...256，具体过程如下：

(6.1)计算Y_t中各个非边缘块观测y_t,i的标准差σ_t,i，i表示块号；

(6.2)对各个非边缘块观测都初始化一个未聚类的标志mark_i，即令mark_i＝0，i表示块号；如果方向子带t下的第i个块观测为某一类的聚类中心块观测的话，则我们用符号a_t,i来表示这个类；初始时令i＝1，表示从第一块开始执行聚类操作；

(6.3)若第i个块观测为非边缘块观测且该块未被聚类，即mark_i＝0，则转(6.4)，否则，转(6.6)

(6.4)将第i个块观测作为类a_t,i的聚类中心块观测；

(6.5)观察第i块的八个邻域块中的所有块观测：假设j表示第i块的八个邻域块中的其中一个块观测，若第j个块观测为非边缘块观测且mark_j＝0，表示第j个块观测未被聚类，则计算第j个块观测的标准差σ_t,j与第i个块观测的标准差σ_t,i的差值，即C_j＝σ_t,i-σ_t,j，若|C_j|≤τ(j∈E)，其中τ为阈值，τ=0.01，则将第j个块观测加入聚类中心为第i个块观测的类a_t,i中，令mark_j＝1，表示该块已被聚类；

(6.6)令i＝i+1,若i≤256，则转(6.3),否则，表示所有块都已聚类完毕，统计该方向上聚类的类别数，记作c_t；

步骤七，按照种群初始化策略分别初始化三个高频子带下各个聚类中心块观测及各个边缘块观测对应的系数块种群Q＝{q_t,i,j}，其中i表示各聚类中心块观测及各个边缘块观测的块号，j＝1,2,Ln，n为种群规模；；

种群初始化策略表述如下：

(7.1)首先需要说明的是，本文中的观测矩阵采用正交高斯随机观测矩阵。

为了更清楚地表明种群初始化的过程，我们以对大小为512×512的图像A为例来说明。首先，对图像A执行一层小波变换操作，得到一个大小为256×256的低频子带系数和三个大小为256×256的高频子带系数，对这三个高频子带系数执行分块操作，块大小为16×16，并将每一块变换为256×1的列向量，按块号从小到大排列表示为S_v＝[s_v,1 s_v,2L s_v,256]，S_h＝[s_h,1 s_h,2L s_h,256]，S_d＝[s_d,1 s_d,2L s_d,256]，对应的观测为Y_v＝[y_v,1 y_v,2L y_v,256]，Y_h＝[y_h,1 y_h,2L y_h,256]，Y_d＝[y_d,1 y_d,2L y_d,256]，表示如下：

$\begin{array}{c}{Y}_v=\mathrm{\Φ}{S}_v\\ =\left[\begin{array}{ccc}\mathrm{\Φ}{s}_{v,1}& \mathrm{\Φ}{s}_{v,2}L& \mathrm{\Φ}{s}_{v,256}\end{array}\right]\\ =\left[\begin{array}{ccc}{y}_{v,1}& y_{v,2}L& y_{v,256}\end{array}\right]\end{array}$

$\begin{array}{c}{Y}_h=\mathrm{\Φ}{S}_h\\ =\left[\begin{array}{ccc}\mathrm{\Φ}{s}_{h,1}& \mathrm{\Φ}{s}_{h,2}L& \mathrm{\Φ}{s}_{h,256}\end{array}\right]\\ =\left[\begin{array}{ccc}{y}_{h,1}& y_{h,2}L& y_{h,256}\end{array}\right]\end{array}$ $\begin{array}{c}{Y}_d=\mathrm{\Φ}{S}_d\\ =\left[\begin{array}{ccc}\mathrm{\Φ}{s}_{d,1}& \mathrm{\Φ}{s}_{d,2}L& \mathrm{\Φ}{s}_{d,256}\end{array}\right]\\ =\left[\begin{array}{ccc}{y}_{d,1}& y_{d,2}L& y_{d.256}\end{array}\right]\end{array}$

那么，我们就可以根据所得的观测Y_v＝[y_v,1 y_v,2L y_v,256]，Y_h＝[y_h,1 y_h,2L y_h,256]，Y_d＝[y_d,1 y_d,2L y_d,256]和正交高斯随机观测矩阵Φ，得到三个高频子带的系数S′_v，S′_h，S′_d，表示如下：

S′_v＝Φ⁺Y_v＝Φ⁺[y_v,1 y_v,2L y_v,256]＝[s′_v,1 s′_v,2L s′_v,256]

S′_h＝Φ⁺Y_h＝Φ⁺[y_h,1 y_h,2L y_h,256]＝[s′_h,1 s′_h,2L s′_h,256]

S′_d＝Φ⁺Y_d＝Φ⁺[y_d,1 y_d,2L y_d,256]＝[s′_d,1 s′_d,2L s′_d,256]

其中，Φ⁺是正交高斯随机观测矩阵Φ的广义逆。

将三个高频子带的系数S′_v，S′_h，S′_d的各列转化为16×16的块，各块系数的集合表示为S″_v＝{s″_v,i}，S″_h＝{s″_h,i}，S″_d＝{s″_d,i}，i表示块号。

(7.2)具体地，以对垂直方向第m类（对垂直方向而言m＝1,2,L c₁，对水平方向而言m＝1,2,L c₂，对对角方向而言m＝1,2,L c₃）的聚类中心块（块号为i）对应的系数块种群Q＝{q_v,i,j}（j＝1,2,L n）的初始化为例，该块对应的位置块为p″_v,i，从(8.1)可以提取到该块对应的系数为s″_v,i，那么系数块种群Q_v,i中第j个个体是按下式初始化的：

$q_{v,i,j}(m_1,m_2)=\left\{\begin{array}{ccc}{s}_{v,i}^{\′\′}(m_1,m_2)& \mathrm{if}& p_{v,i}^{\′\′}(m_1,m_2)=0\\ L\×s_{v,i}^{\′\′}(m_1,m_2)& \mathrm{if}& p_{v,i}^{\′\′}(m_1,m_2)=1\end{array}\right.$

其中，q_v,i,j(m₁,m₂)为第j个系数个体q_v,i,j第m₁行第m₂列的值，m₁＝1,2,L 16,m₂＝1,2,L 16,L是取自区间[1,1.5]的随机数，正是由于L的随机性，当L取不同的值时，便可获得不同的系数个体，组成该块对应的系数块种群Q_v,i。

对水平方向和垂直方向，我们也用相应的位置块p″_h,i和系数s″_h,i及位置块p″_d,i和系数s″_d,i，通过上述方法得到各块对应的系数块种群Q＝{q_h,i,j}和Q＝{q_d,i,j}，j＝1,2,L n。

(7.3)各个边缘块对应的种群的初始化策略同上，只需将各个聚类中心块替换为各个边缘块，位置块也调用各边缘块对应的位置块；

步骤八，对三个高频子带下的各个系数块种群Q＝{q_t,i,j}执行交叉操作，得到交叉后的子代系数块种群Q′＝{q′_t,i,j}；

以对垂直方向第i块对应的系数块种群Q＝{q_v,i,j}的交叉操作为例来说明，假设块大小为16×16，具体过程为：将该种群中的n个个体，两两配对作为待交叉的配对个体，假设两个系数个体和凑为一对，其中，j₁,j₂∈{1,2,L n}，并从或中随机选择一个基因位，交换这两个个体和的以该基因位为中心的5×5的区域，得到两个新的个体和，组成对应于垂直方向第i块的子代系数块种群Q′＝{q′_v,i,j}，j＝1,2,L n+2。

对水平方向和对角方向的各个块对应的系数块种群Q＝{q_h,i,j}和Q＝{q_d,i,j}，也按照上述方法，生成对应的交叉操作后的子代系数块种群Q′＝{q′_h,i,j}和Q′＝{q′_d,i,j}，j＝1,2,L n+2。

步骤九，对三个高频子带下的各个子代系数块种群Q′＝{q′_t,i,j}执行变异操作；

需要说明的是，变异操作是对系数块种群中的每个个体的某个基因位按概率进行的。以对垂直方向第i块对应的子代系数块种群Q′＝{q′_v,i,j}的变异操作为例来说明，假设块大小为16×16，对子代系数块种群Q′＝{q′_v,i,j}中第j个系数个体q′_v,i,j的变异操作过程为：从系数个体q′_v,i,j中随机选择一个基因位，再从区间(0,1)中随机生成一个数，若随机数小于给定的概率p，则观察该块对应的位置块P′_v,i中所选基因位处的值，若位置块中该处的值为1，则从区间[min_v,i,-λ_v,i]或者[λ_v,i,max_v,i]中随机选择一个数作为系数个体q′_v,i,j当前位置处的值，而若位置块中该处的值为0，则从区间[-λ_v,i,λ_v,i]中随机选择一个数作为系数个体q′_v,i,j当前位置处的值;反之，如果生成的随机数大于给定的概率p，则不进行任何操作。其中，λ_v,i为步骤二中所得的垂直方向第i块的分界参数，min_v,i为种群初始化时，第i块对应的系数块种群所有系数个体所含元素的最小值，max_v,i种群初始化时，第i块对应的系数块种群所有系数个体所含元素的最大值。对种群中各个个体执行变异操作，得到变异操作后的新子代系数块种群Q″＝{q″_v,i,j}。

对水平方向和对角方向的各个块对应的子代系数块种群Q′＝{q′_h,i,j}和Q′＝{q′_d,i,j}，也按照上述方法，生成对应的变异操作后的新子代系数块种群Q″＝{q″_h,i,j}和Q″＝{q″_d,i,j}。

步骤十，对各聚类中心的系数块和边缘的系数块分别定义两种不同的适应度函数以对子代种群Q″执行相应的种群更新操作；

(10.1)若当前块为非边缘块，也即当前块为一个聚类中心块，则按以下适应度函数计算各个个体的适应度f(q″_t,i,j)，表示如下：

$f\left(q_{t,i,j}^{\′\′}\right)=\frac{1}{\left|a_{t,i}\right|}\underset{k\∈a_{t,i}}{\mathrm{\Σ}}\frac{1}{{\left|\right|y_{t,k}-\mathrm{\Φ}{q}_{t,i,j}^{\′\′}\left|\right|}_2^2},\mathrm{t\∈}\{v,h,d\}$

其中，t∈{v,h,d}，表示的是该块所属子带的方向，其中v表示垂直方向，h表示为水平方向，d表示对角方向；由步骤（6）可得，a_t,i表示方向为t的子带中以第i块作为聚类中心块的类；y_t,k为方向为t的子带中类a_t,i中的第k块对应的观测，q″_t,i,j为该聚类中心块所对应的种群Q″中的第j个系数个体，f(q″_t,i,j)即为个体q″_t,i,j的适应度；

若当前块为边缘块，则按以下适应度函数计算各个个体的适应度f(q″_t,i,j)，表示如下：

$f\left(q_{t,i,j}^{\′\′}\right)=\frac{1}{{\left|\right|y_{t,i}-\mathrm{\Φ}{q}_{t,i,j}^{\′\′}\left|\right|}_2^2},t\∈\{v,h,d\}$

其中，t∈{v,h,d}，表示的是该块所属子带的方向，其中v表示垂直方向，h表示为水平方向，d表示对角方向；y_t,i为方向为t的子带中第i块对应的观测，q″_t,i,j为该块所对应的种群Q″中的第j个系数个体，f(q″_t,i,j)即为个体q″_t,i,j的适应度；

(10.2)将Q″中所有个体的适应度与上一代种群Q中所有个体的适应度相比较，从中选择前n个适应度较大的个体更新种群Q，并将适应度最大的个体作为本代进化的最优个体。

步骤十一，分别对三个子带下的各个块选择出的最优系数，若该块为非边缘块，则将其对应的最优系数作为其同类的各个块的最优系数，然后对所有块对应的最优系数执行合并块操作，形成各个子带对应的系数B′_v（垂直方向），B′_h（水平方向），B′_d（对角方向）；

步骤十二，若进化代数满足停止条件，则转步骤十三，否则，利用各聚类中心块和边缘块选择出的最优个体更新该块对应的位置块和系数块种群，转入步骤九；

以对垂直方向第i块的位置块p″_v,i和系数块种群Q＝{q_t,i,j}的更新为例。整个更新过程分三步进行：

12.1）首先，由步骤十中，完成了对系数块种群Q的更新。

12.2）假设q″_v,i,j为步骤十中选择出的该块对应的最优系数个体，我们希望该个体的某些优秀的特性能够扩展影响整个系数块种群，所以首先利用该个体对其对应的位置矩阵更新，而系数块种群中的所有个体都需要保证始终在对应位置矩阵的约束之下，所以利用更新后的位置矩阵对系数块种群中的其他个体进行更新，从而，使得整个种群的特性向好的方向发展。

具体来说，首先初始化一个全零的大小与位置块大小相等的矩阵p″_v,i，逐个检查最优个体q″_v,i,j中的每个值，若某位置处的模值大于步骤三中得到的该块对应的分阈值λ_v,i，则将矩阵p″′_v，i中该位置处的值置为1，否则，将矩阵p″′_v，i中该位置处的值置为0，将所有值都比较完成后，得到由0和1组成的矩阵p″′_v,i；然后再逐个位置比较矩阵p″′_v,i与步骤六中该块初始化时的位置块p′_v,i，若初始化位置块p′_v,i中某位置处的值为1，而矩阵p″′_v,i中对应位置处的值为0，则将矩阵p″′_v,i中该位置处的值置1。如此，得到的矩阵p″′_v,i作为该块对应的位置块更新p″_v,i，即p″_v,i=p″′_v,i，这样做即保证了保证边缘提取得到位置信息得以保留，而且加入了新的最优个体带来的新的位置信息。

12.3）接着需要利用新的位置块p″_v,i对系数块种群Q＝{q_t,i,j}中的其他个体进行更新，具体分为以下几种情况：

（1）若位置块p″_v,i中某位置处值为1，意指该位置处应为一个大系数，而系数个体q_v,i,j中对应位置处为一个小系数，则需要对该位置处的系数进行调整：

（1a）若该位置处系数数个体q_v,i,j中对应位置处为一个正的小系数，则将对应位置处的系数加上该块对应的分界参数λ_v,i形成一个正的大系数，代替原有的小系数；

（1b）若该位置处系数数个体q_v,i,j中对应位置处为一个负的小系数，则将对应位置处的系数减去该块对应的分界参数λ_v,i形成一个负的大系数，代替原有的小系数；

（2）若位置块p″_v,i中某位置处值为0，意指该位置处应为一个小系数，而系数数个体q_v,i,j中对应位置处为一个大系数，则需要对该位置处的系数进行调整：

（2a）若该位置处系数数个体q_v,i,j中对应位置处为一个负的大系数，则从区间[-λ_v,i,0]中随机选择一个数形成一个负的小系数，代替原有的大系数；

（2b）若该位置处系数数个体q_v,i,j中对应位置处为一个正的小系数，则从区间[0,λ_v,i]中随机选择一个数形成一个正的小系数，代替原有的大系数；

如此，就完成了对位置块和对应系数块种群的更新，对水平方向和对角方向的各位置块和系数块种群的更新也按上述方法逐步进行即可。

步骤十三，结合保留的低频子带系数Y₀及步骤十一中得到的三个高频子带系数B_t进行小波逆变换，得到重构图像；

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

1.仿真条件：

本发明的仿真在windows7，SPI，CPUIntel(R)Core(TM)2，基本频率3.00GHz，软件平台为MatlabR2011b上运行，仿真选用的是512×512的标准Barbara图像和Lena图像。

2.仿真内容与结果：

本仿真中，使用正交匹配追踪算法OMP，基追踪算法BP，迭代硬阈值算法IHT和本发明方法对大小为512×512的Barbara图像和Lena图像在采样率分别为40%和35%的条件下进行图像重构，其中OMP和IHT方法都是在小波域下进行重构，重构结果如图所述：

图2(a)是Barbara原始图像；

图2(b)是图2(a)的局部放大图；

图2(c)是本发明得到的Barbara重构图；

图2(d)是图2(c)的局部放大图；

图2(e)是OMP得到的Barbara的重构图；

图2(f)是图2(e)的局部放大图；

图2(g)是BP得到的Barbara重构图；

图2(h)是图2(g)的局部放大图；

图2(i)是IHT得到的Barbara重构图；

图2(j)是图2(i)的局部放大图；

图4(a)是Lena原始图像；

图4(b)是图4(a)的局部放大图；

图4(c)是本发明得到的Lena重构图；

图4(d)是图4(c)的局部放大图；

图4(e)是OMP得到的Lena的重构图；

图4(f)是图4(e)的局部放大图；

图4(g)是BP得到的Lena重构图；

图4(h)是图4(g)的局部放大图；

图4(i)是IHT得到的Lena重构图；

图4(j)是图4(i)的局部放大图；

从重构图和局部放大图可以看出，本发明的重构图像的边缘部分保持的较好，平滑部分也优于OMP，BP和IHT的重构图像。

用现有的OMP，BP，IHT和本发明方法分别在采样率为30%、35%、40%、45%的情况下，对大小为512×512的Barbara图像、Lena像和Boat图像做仿真实验，各算法对应各图10次重构结果的峰值信噪比PSNR的平均值如表1所示。

表1各采样率下PSNR值

从表1数据可以看出，本发明方法在采样率为30%、35%、40%、45%下得到的结果图的峰值信噪比PSNR都要高于OMP，BP和IHT方法得到的PSNR，即本发明的方法的重构图像质量要比OMP，BP和IHT方法高。

从图2中本发明重构结果图与OMP，BP和IHT方法重构结果图以及各自的局部放大图可以看出，本发明有围巾等纹理处的重构细节更好，重构质量更高。

由图3和图5可以看出，本发明方法得到的重构结果图的PSNR值高于其他方法。

综上，本发明能够很好地得到清晰的图像，与现有的其他重构方法相比，本发明提高了图像的重构质量。

Claims (4)

1.一种基于联合稀疏和先验约束的压缩感知图像重构方法，其特征在于：包括如下步骤： 

(1)输入低频子带Y₀和三个高频子带的分块观测Y_t，其中t∈{v,h,d}，表示子带的方向，其中v表示垂直方向，h表示为水平方向，d表示对角方向； 

(2)获得三个高频子带对应的位置矩阵P_t； 

(2a)将三个高频子带系数置0，结合低频子带Y₀，做小波逆变换，得到一幅边缘模糊的图像I1； 

(2b)用canny算子对边缘模糊的图像I1进行边缘检测，得到一幅只含有边缘信息的图像I2； 

(2c)对只含有边缘信息的图像I2执行一层小波变换，得到一个低频子带Y₁和三个含有边缘信息的高频子带为S_t； 

(2d)以作为方向为t的子带的阈值,将对应子带S_t中模值大于对应阈值的位置标记为1，意为该位置处为大系数，而将模值小于对应阈值的位置标记为0，意为该位置处为小系数，这样获得三个高频子带对应的位置矩阵P_t，其中μ_t为子带S_t的模值的均值,为伸缩因子，本文中取

(3)对位置矩阵P_t执行提取疫苗和注射疫苗的操作，得到位置矩阵P′_t； 

(4)对位置矩阵P′_t执行分块操作，得到三个高频子带下各块观测对应的位置块p′_t,i； 

为了记录初始位置块p′_t,i的信息，将p′_t,i保存为变量p″_t，i，即p″_t,i=p′_t,i，后续只对变量p″_t,i进行操作； 

其中i表示块号，若图像大小为512×512，则i＝1,2,...256，若图像大小为256×256，则i＝1,2,...128 

(5)根据得到的三个子带下各块观测的位置块p″_t,i对各块观测进行边缘块观测和非边缘块观测的划分； 

(6)根据第（5）步的边缘块观测和非边缘块观测的划分结果，对三个高频子带的分块观测Y_t中所有非边缘块观测执行局部相似聚类，得到对应于三个子带中非边缘块观测聚类的集合： 

其中，a_t,i表示方向为t的子带的第i类对应的集合，其中i＝1,2,Lc_t，c_t表示方向为t的子带聚类的类别数； 

(7)按照种群初始化策略分别初始化三个高频子带下各个聚类中心块观测及各个边缘块观测对应的系数块种群Q＝{q_t,i,j}，其中i表示各聚类中心块观测及各个边缘块观测的块号，j＝1,2,Ln，n为种群规模； 

(8)对三个高频子带下的各个系数块种群Q＝{q_t,i,j}执行交叉操作，得到交叉后的子代系数块种群Q′＝{q′_t,i,j}； 

(9)对三个高频子带下的各个子代系数块种群Q′执行变异操作，得到变异后的子代系数块种群Q″＝{q″_t,i,j}； 

(10)对各聚类中心的系数块和边缘的系数块分别定义两种不同的适应度函数以对子代种群Q″执行相应的种群更新操作； 

(11)分别从三个子带下各个子代种群Q″中选择出各个系数块对应的最优系数个体，若该系数块为非边缘块，则将其对应的最优系数个体作为其同类的各个系数块的最优系数块，然后对所有最优系数块执行合并块操作，形成各子带对应的系数B_t； 

(12)若进化代数满足停止条件，则转步骤（13），否则，利用各个聚类中心系数块和边缘的系数块选择出的最优系数个体更新该块对应的位置块和系数块种群，转入步骤(8)； 

(13)结合保留的低频子带系数Y₀及步骤(11)中得到的三个高频子带系数B_t，进行小波逆变换，得到重构图像。 

2.根据权利要求1所述的一种基于联合稀疏和先验约束的压缩感知图像重构方法，其特征在于：所述步骤(5)中根据得到的三个子带下各块观测的位置块p″_t,i对各块观测进行边缘块观测和非边缘块观测的划分，若某块观测对应的位置块p″_t,i为全0矩阵，则将该块观测划作非边缘块观测，否则，将其作为边缘块观测。 

3.根据权利要求1所述的一种基于联合稀疏和先验约束的压缩感知图像重构方法，其特征在于：所述步骤(6)中对三个高频子带的分块观测Y_t中所有非边缘块对应的观测执行局部相似聚类，t∈{v,h,d},以对大小为512×512的图像进行操作为例即块号i＝1,2,...256，具体过程如 下： 

(6.1)计算Y_t中各个非边缘块观测y_t,i的标准差σ_t,i，i表示块号； 

(6.2)对各个非边缘块观测都初始化一个未聚类的标志mark_i，即令mark_i＝0，i表示块号；如果方向子带t下的第i个块观测为某一类的聚类中心块观测的话，则我们用符号a_t,i来表示这个类；初始时令i＝1，表示从第一块开始执行聚类操作； 

(6.3)若第i个块观测为非边缘块观测且该块未被聚类，即mark_i＝0，则转(6.4)，否则，转(6.6) 

(6.4)将第i个块观测作为类a_t,i的聚类中心块观测； 

(6.5)观察第i块的八个邻域块中的所有块观测：假设j表示第i块的八个邻域块中的其中一个块观测，若第j个块观测为非边缘块观测且mark_j＝0，表示第j个块观测未被聚类，则计算第j个块观测的标准差σ_t,j与第i个块观测的标准差σ_t,i的差值，即C_j＝σ_t,i-σ_t,j，若|C_j|≤τ，其中τ为阈值，τ=0.01，则将第j个块观测加入聚类中心为第i个块观测的类a_t,i中，令mark_j＝1，表示该块已被聚类； 

(6.6)令i＝i+1,若i≤256，则转(6.3),否则，表示所有块都已聚类完毕，统计该方向上聚类的类别数，记作c_t。

4.根据权利要求1所述的一种基于联合稀疏和先验约束的压缩感知图像重构方法，其特征在于：所述步骤(10)中对各聚类中心的系数块和边缘的系数块分别定义两种不同的适应度函数以对子代种群Q″执行相应的种群更新操作按如下过程进行： 

(10.1)若当前块为非边缘块，也即当前块为一个聚类中心块，则按以下适应度函数计算各个个体的适应度f(q″_t,i,j)，表示如下： 

其中，t∈{v,h,d}，表示的是该块所属子带的方向，其中v表示垂直方向，h表示为水平方向，d表示对角方向；由步骤（6）可得，a_t,i表示方向为t的子带中以第i块作为聚类中心 块的类；y_t,k为方向为t的子带中类a_t,i中的第k块对应的观测，q″_t,i,j为该聚类中心块所对应的种群Q″中的第j个系数个体，f(q″_t,i,j)即为个体q″_t,i,j的适应度； 

若当前块为边缘块，则按以下适应度函数计算各个个体的适应度f(q″_t,i,j)，表示如下： 

其中，t∈{v,h,d}，表示的是该块所属子带的方向，其中v表示垂直方向，h表示为水平方向，d表示对角方向；y_t,i为方向为t的子带中第i块对应的观测，q″_t,i,j为该块所对应的种群Q″中的第j个系数个体，f(q″_t,i,j)即为个体q″_t,i,j的适应度； 

(10.2)将Q″中所有个体的适应度与上一代种群Q中所有个体的适应度相比较，从中选择前n个适应度较大的个体更新种群Q，并将适应度最大的个体作为本代进化的最优个体。