CN105354800A - 基于图像结构的粒子群优化非凸压缩感知图像重构方法 - Google Patents

Abstract

本发明公开了一种基于图像结构的粒子群优化非凸压缩感知重构方法，其实现方案是：1.根据观测向量判别图像块结构，将图像块标记为光滑，单方向和多方向；3.对不同类型的图像块所对应的观测向量聚类，并为每一类图像块构造相应的过完备冗余字典；4.为每一类图像块构造初始种群；5.对每一类光滑图像块利用基于分组初始化策略的粒子群算法搜索尺度最优原子组合；6.对每一类单方向和多方向图像块，利用基于交叉和原子方向约束的粒子群算法搜索方向和尺度最优原子组合；7.计算所有图像块的估计值并按顺序拼接成整幅图像输出。本发明重构时间短，重构出的图像具有视觉效果好、峰值信噪比和结构相似度高，可用于图像信号的非凸压缩感知重构。

Description

基于图像结构的粒子群优化非凸压缩感知图像重构方法

技术领域

本发明属于图像处理技术领域，更进一步涉及一种压缩感知重构方法，可用于从图像的压缩观测中获得高质量清晰图像。

背景技术

压缩感知CS理论的出现打破了传统的奈圭斯特采样定理，CS理论指出信号可以进行低速采样和少量采样，并且可以精确重构，采样速率不再决定于带宽，而决定于信息在信号中的结构和内容。压缩感知的研究主要包括三个方面：压缩观测、稀疏表示和重构方法，其中，重构方法研究从信号的压缩观测中精确重构出原信号，是压缩感知中的核心问题。

压缩感知重构的本源问题是l₀范数约束的非凸优化问题。在“F.LiuandL.Lin,etc..Non-convexCompressedSensingbyNature-inspiredOptimizationAlgorithms,IEEETransactionsonCybernetics,DOI:10.1109/TCYB.2014.2343618”一文中提出了一种非凸压缩感知图像重构方法，该方法采用两阶段的重构模型，分别使用遗传优化算法和克隆选择算法获得图像块在字典方向上和尺度位移上的最优原子组合。西安电子科技大学的专利申请“基于冗余字典和结构稀疏的非凸压缩感知图像重构方法”(公开号：CN103295198A，申请号：CN201310192104，申请日：2013年5月13日)中公开了一种基于非凸模型的图像压缩重构方法，该方法采用与文章方法相似的重构模型，提出了另一种非凸压缩感知图像重构方法。上述两种技术都使用了智能优化方法，对l₀范数约束的非凸重构问题采用全局优化策略，获得了比采用局部优化策略的经典贪婪算法，即OMP方法和IHT方法，更优的原子组合和图像重构效果。

但是上述两种方法均存在重构速度慢的不足，两种方法均采用了遗传算法和克隆选择算法两阶段进行优化，两种方法对图像块的结构进行了光滑和非光滑的判别，且第二种方法在种群初始化的时候没有利用图像结构的先验信息对初始化种群进行优化，速度较慢，不利于实时应用。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于图像结构的粒子群优化非凸压缩感知图像重构方法，以提高重构速度，便于实时应用。

本发明的技术方案是挖掘图像块的结构信息，为不同结构特征的图像块采用不同的优化策略，在不损失重构精度的前提下，提高重构速度。其实现步骤包括如下：

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

1)压缩感知接收方接收发送方发送过来的所有图像块的观测向量y₁，y₂,...,y_n，其中，n表示图像块的数目；

2)对每个观测向量所对应的图像块进行结构判别，将图像块标记为光滑块、单方向块和多方向块，并记录单方向块的方向；

3)分别对不同结构类型的图像块所对应的观测向量进行聚类，并根据脊波过完备字典为每一类观测向量所对应的图像块构造一个相应的过完备冗余字典；

4)根据步骤3)得到的过完备字典以及图像块的结构类型分别为每一类光滑类图像块、单方向类图像块和多方向类图像块构造相应的粒子群初始种群；

5)对每一类光滑类图像块，在其对应的过完备冗余字典下，利用基于分组初始化策略的搜索尺度最优原子组合，得到每一类光滑类图像块的最优原子组合；

6)对每一类单方向类和多方向类图像块，分别在其对应的过完备冗余字典下，利用基于交叉和原子方向约束的粒子群算法搜索粒子在方向和尺度上的最优原子组合，得到每一类图像块的最优原子组合；

7)根据步骤5)、6)得到的所有观测向量的最优原子组合，计算所有图像块的估计值，并按顺序拼接成一整幅重构图像输出。

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

1.本发明充分挖掘了图像块的结构特征，针对光滑类图像块和单方向类图像块构造了小规模的过完备字典，减小了压缩感知重构算法的搜索范围，从而有效减少重构时间。

2.本发明充分利用了光滑块对尺度参数更敏感且其方向具有多样性和随机的特点，搜索光滑块的尺度最优原子组合，提高了光滑块的重构精度，缩短了重构时间。

3.本发明采用了一种比遗传算法和克隆选择算法速度更快的粒子群算法重构图像，且不再使用两阶段重构方法，大大缩短了重构时间。

下面结合实施例附图对本发明作进一步说明

附图说明

图1是本发明的图像重构流程示意图；

图2是用本发明及两种现有方法在30％采样率下分别对Lena图的重构结果图；

图3是用本发明及两种现有方法在30％采样率下分别对Barbara图的重构结果图；

图4用本发明及两种现有方法的平均运行时间对比图。

具体实施方式

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

步骤1，压缩感知接收方接收发送方发送过来的所有图像块的观测向量。

发送方通过输入一幅512×512的自然图像，将其分成16×16的不重叠块，得到1024个图像块，再对每个图像块进行高斯随机观测，得到所有图像块的观测向量y₁，y₂,...,y_n，并发送所有图像块的观测向量，其中，n表示图像块的数目。

步骤2，对每个观测向量所对应的图像块进行结构判别，将图像块标记为光滑块、单方向块和多方向块，并记录单方向块的方向。

2a)计算每个观测向量的方差，并设置光滑阈值为0.45σ，其中，σ是所有观测向量的方差的平均值；

2b)对每个观测向量，判断其方差是否小于等于光滑阈值，若是，则将该观测向量对应的图像块标记为光滑块，否则，不做标记；

2c)对每个尚未标记的图像块的观测向量，进行以下判定：

2c1)根据脊波函数生成脊波过完备冗余字典，该字典共有12032个原子，36个方向和16个尺度，每个尺度下大约有20个左右的位移；

2c2)将脊波过完备冗余字典按方向分为36个子字典ψ₁，ψ₂，...，ψ_i,...,ψ₃₆，对每一个图像块，根据其观测向量y和这些子字典，计算得到一个观测残差序列r₁，r₂，...，r_i,...,r₃₆，找到序列中的最小值在序列中的位置i，i＝1,2,...,36；

2c3)按如下公式计算序列中的子字典ψ_i的观测残差r_i：

$r_i=\left|\right|y-{\mathrm{\&Phi;D}}_r\&lsqb;{\left({\mathrm{\&Phi;D}}_r\right)}^{+}y\&rsqb;|{|}_2^2$ <1>

式中，y是待判定图像块的观测向量，Φ是用于观测的高斯矩阵，D_r是子字典ψ_i中与y相关性最大的10个原子组合，(·)⁺表示计算矩阵的伪逆矩阵，是向量二范数的平方；

2c4)利用观测残差序列中位置i-2，i-1，i，i+1和i+2上的五个残差值r_i-2，r_i-1，r_i，r_i+1和r_i+2，对观测向量y对应的图像块进行标记：若r_i-2大于r_i-1，r_i-1大于1.2r_i，r_i+1大于1.2r_i，且r_i+2大于r_i+1，则将观测向量y对应的图像块标记为单方向，并记录该单方向块的方向为第i个方向，否则，不做标记；

其中，在判断条件中，若i为1，r_i-1和r_i-2分别用r₃₆和r₃₅代替，若i为2，r_i-2用r₃₆代替，若i为36，r_i+1和r_i+2分别用r₁和r₂代替，若i为35，r_i+2用r₁代替；

2c5)将其它尚未标记的图像块标记为多方向。

步骤3，分别对不同结构类型的图像块所对应的观测向量进行聚类，并根据脊波过完备字典为每一类观测向量所对应的图像块构造一个相应的过完备冗余字典。

聚类方法有多种，例如C均值聚类方法，模糊聚类方法，仿射聚类方法等，本实施例中，使用的聚类方法为仿射聚类方法，其具体实现如下：

3a)对所有光滑图像块对应的观测向量进行聚类，得到光滑图像块的聚类；

3b)对所有单方向图像块对应的观测向量进行聚类，得到单方向图像块的聚类；

3c)对所有多方向图像块对应的观测向量进行聚类，得到多方向图像块的聚类；

冗余字典有多种，例如主成分分析字典，曲线波过完备冗余字典，脊波过完备冗余字典等，其中，主成分分析字典对图像的主要结构表示较好，曲线波过完备冗余字典对曲线表示较好，而脊波过完备冗余字典对直线表示较好；

在本实施例中，所采用的字典为脊波过完备冗余字典，大小为12032，对这12032个基原子按整数编号为1，2，3，...，12032，字典由36个不同方向的子字典按顺序组成，每个方向子字典有16个尺度，每个尺度有15至25个位移，其步骤如下：

3d)为光滑类图像块构造过完备字典Ψ_s，该过完备字典由脊波冗余字典中包含所有方向的前5个尺度的脊波冗余子字典组成， 表示包含所有方向的尺度为h的脊波冗余子字典，其中，h＝1,2,..,5；依次将中的原子按原子编号从小到大映射到一个连续的正整数空间中，将作为光滑类图像块的过完备冗余字典Ψ_s的映射字典，其中，m是光滑类图像块的过完备冗余字典Ψ_s的大小；

3e)为单方向类图像块构造过完备冗余字典Ψ_g:从脊波冗余字典中取出第r个方向及与第r个方向左右相邻的4个方向的子字典作为该单方向类图像块的过完备冗余字典Ψ_g，其中，r表示该单方向类图像块的方向，表示第r个方向的脊波冗余子字典，r＝1,2,...,36，若r为1，则r-1和r-2分别用36和35代替，若r为2，则r-2用36代替，若r为36，则r+1和r+2分别用1和2代替，若r为35，则r+2用1代替；

3f)将整个脊波过完备冗余字典作为多方向类图像块的字典。

步骤4，根据步骤3)得到的过完备字典以及图像块的结构类型分别为每一类光滑类图像块、单方向类图像块和多方向类图像块构造相应的粒子群初始种群。

粒子具有速度，对一个规模为n的粒子群A，A＝{b₁,b₂,...,b_i,...,b_n}，对应着一个速度V，V＝{v₁,v₂,...,v_i,...,v_n}，其中b_i表示第i个粒子，v_i表示粒子b_i的速度，b_i与v_i一一对应，i＝1,2,...,n，n为种群大小；

构造粒子群，需要对粒子b_i和粒子的速度v_i进行编码，常用的编码方式有二进制编码，格雷码方法编码，实数编码等；

本实施例中，对粒子b_i采用整型实数编码，并且整型实数范围为脊波过完备冗余字典的原子范围，对粒子的速度v_i采用实数编码，其实现步骤如下：

4a)对每一光滑类图像块，根据光滑类图像块的过完备冗余字典Ψ_s的映射字典按照该映射字典的5个子字典，初始化一个规模为20的粒子群A₁，并将该粒子群A₁均分成5个组,每个组一个尺度；随机初始化每个组的粒子方向，同时保证每个粒子至少包含20个方向，并将每个组中粒子的速度随机初始化为中的任意一个值，其中L^h为映射字典中的子字典的大小，h＝1,2,..,5；

4b)对每一单方向类图像块，根据其过完备冗余字典Ψ_g的5个子字典，初始化一个规模为20的粒子群A₂，并将该粒子群A₂均分成5个组，每个组代表一个方向，随机初始化每个组中的粒子，并将每个粒子上的原子按原子编号从小到大有序排列，再将粒子的速度随机初始化为中的任意一个值，其中L^r为该单方向类过完备冗余字典Ψ_g中的子字典的大小，r＝1,2,...,36；

4c)对多方向类图像块，初始化一个规模为36的粒子群A₃，每个粒子代表两个方向，其中一个方向从脊波过完备冗余字典中选择与该粒子编号相同的方向，另一个方向则从剩下的方向中随机选择一个，并将每个粒子上的原子也按原子编号从小到大有序排列，再将粒子的速度随机初始化为[-Lⁱ,Lⁱ]×0.5中的任意一个值，其中Lⁱ表示所选择的两个代表方向的脊波过完备冗余子字典的大小，i表示粒子编号，i＝1,2,...,36；

4d)对光滑类图像块的粒子群A₁，将其中的粒子映射回光滑类过完备冗余字典Ψ_s中；

4e)按照如下公式，计算光滑类、单方向类和多方向类粒子群中每个粒子的适应度值：

$f\left(b_i\right)=\frac{1}{\underset{j=1}{\overset{q}{\&Sigma;}}\left|\right|y_j-{\mathrm{\&Phi;b}}_i{\mathrm{\&alpha;}}_j|{|}_2^2},$ <2>

其中f(b_i)为该类观测向量对应的图像块的粒子群中第i个粒子的适应度值，b_i是粒子群中第i个粒子所代表的解，对应着脊波过完备冗余字典中的一组基原子组合，即一个子字典，j为观测向量的标号，q为聚类后该类的总测量数，y_j为类中的第j个测量，Φ为高斯观测矩阵，α_j为该类观测向量中第j个测量对应的图像块的稀疏表示系数，它是由传感矩阵的广义逆矩阵和该观测向量y_j相乘得到，传感矩阵由高斯观测矩阵Φ和子字典b_i相乘得到，是向量二范数的平方；

4f)初始化粒子群的历史最优解以及全局历史最优解：将粒子群中每个粒子的历史最优解初始化为粒子自身，找出整个粒子群中适应度值最大的粒子，将其作为全局历史最优解。

步骤5，对每一类光滑类图像块，在其对应的过完备冗余字典下，利用基于分组初始化策略的粒子群算法搜索尺度最优原子组合，得到每一类光滑类图像块的最优原子组合。

本步骤的具体实现如下：

5a)对粒子群算法中的参数进行设置，将粒子群的粒子速度位置更新为：

$v_{ij}^{t+1}={\mathrm{wv}}_{ij}^t+c_1\mathrm{\&xi;}(p_{ij}^t-b_{ij}^t)+c_2\mathrm{\&eta;}(p_{gj}^t-b_{ij}^t),$

<3>

$k_{ij}^{t+1}=b_{ij}^t+v_{ij}^{t+1}.$

其中，表示第t+1次迭代更新后粒子b_i第j个分量对应的速度，表示第t+1次迭代更新后粒子b_i第j个分量的取值，b_i表示第i个粒子，v_i表示第i个粒子b_i对应的速度，b_i与v_i具有相同的维数，两者一一对应；表示粒子群中第i个粒子b_i的第j个分量，表示粒子群中第i个粒子对应的速度v_i的第j个分量，表示第i个粒子的历史最优位置的第j个分量，表示全局历史最优位置的第j个分量，t表示当前是第几次迭代；w表示惯性权重，对w采取了线性递减策略，w＝w_max-(w_max-w_min)×t/iter_max，其中iter_max表示总迭代次数，w_max表示w的最大值，w_min表示w的最小值，w从0.9线性递减至0.4；c₁表示粒子跟踪自己本身历史最优值的权重系数，它表示粒子自身的认识，叫做“认知”，c₁＝c_{1_max}-(c_{1_max}-c_{1_min})×t/iter_max，其中c_{1_max}和c_{1_min}分别表示c₁的最大值和最小值，c₁从2.5线性递减至0.5；c₂表示粒子跟踪群体最优值的权重系数，即粒子对整个群体知识的认识，叫做“社会知识”，也叫做“社会”，c₂＝c_{2_min}+(c_{2_max}-c_{2_min})×t/iter_max，其中c_{2_max}和c_{2_min}分别表示c₂的最大值和最小值，c₂从0.5线性递增至2.5；ξ、η是[0,1]区间内均匀分布的两个随机数；

5b)根据粒子群的粒子速度位置更新公式计算粒子的新速度，对得到的粒子的新速度进行控制：若 $v_{ij}^{t+1}<(1-b_{ij}^t)\&times;0.1,$ 则令 $v_{ij}^{t+1}=(1-b_{ij}^t)\&times;0.1;$ 若 $v_{ij}^{t+1}>(m-b_{ij}^t)\&times;0.1,$ 则令然后再计算粒子的新位置，对得到的粒子的新位置进行取整操作，使其映射在光滑类图像块的过完备冗余字典的映射字典中，其中，m是映射字典的大小；

5c)计算更新后的粒子的适应度值，将更新后粒子新的适应度值与粒子自身的历史最优值进行比较，若大于历史最优值，则用更新后的粒子代替历史最优解，并从粒子群中所有粒子的历史最优解中找出适应度值最大的历史最优解，若其适应度值大于全局历史最优值，则用该粒子的历史最优解代替全局历史最优解；

5d)判断是否满足粒子群算法的迭代停止条件，若满足，则执行步骤5e)，若不满足，返回步骤5b)，继续对粒子进行速度位置更新操作；

5e)将粒子群的全局历史最优解映射回光滑类图像块的过完备冗余字典Ψ_s中，得到真实的全局历史最优解，将该全局历史最优解原子组合作为该光滑类图像块的重构基。

步骤6，对每一类单方向类和多方向类图像块，分别在其对应的过完备冗余字典下，利用基于交叉和原子方向约束的粒子群算法搜索粒子在方向和尺度上的最优原子组合，得到每一类图像块的最优原子组合。

本步骤的具体实现如下：

6a)对粒子群算法中的参数进行设置：将惯性权重w从0.9线性递减至0.4，将粒子跟踪自己本身历史最优值的权重系数c₁从2.5线性递减至0.5，将粒子跟踪群体最优值的权重系数c₂从0.5线性递增至2.5，设ξ、η分别为[0,1]区间内均匀分布的两个随机数，增设一个控制粒子速度的参数P_r，P_r随着进化代数从1.0线性递减至0.1；

6b)用粒子速度位置更新公式计算粒子的速度和位置，计算第t次迭代每个粒子每维分量上的原子所在的方向r，得到第r个方向的原子范围[g^r,h^r]；

6c)对更新后粒子的新速度进行控制，约束粒子的方向：若则控制 $v_{ij}^{t+1}=(g^r-b_{ij}^t)\&times;0.5,$ 若 $v_{ij}^{t+1}>(h^r-b_{ij}^t)\&times;0.5,$ 则控制 $v_{ij}^{t+1}=(h^r-b_{ij}^t)\&times;0.5,$ 即粒子的新位置被约束在第r个方向的原子范围[g^r,h^r]内；

6d)生成一个大小为n×k的随机矩阵R＝{r₁,r₂,..,r_i,...,r_n}，r_i＝(r_i1,r_i2,...,r_ij,...,r_ik)，按照下式，计算粒子的新位置：

$b_{ij}^{t+1}=\left\{\begin{array}{cc}{b}_{ij}^t,& r_{ij}\&GreaterEqual;P_r\\ \&lsqb;b_{ij}^t+v_{ij}^{t+1}\&rsqb;,& r_{ij}<P_r\end{array}\right.,$ <4>

其中，表示第t+1次迭代更新后的粒子b_i第j个分量，表示第t次迭代粒子群中第i个粒子b_i的第j个分量，表示对应的速度，i＝1,2,...,n，j＝1,2,...,k，n为粒子群的大小，k为粒子位置的编码长度，即稀疏度，[·]表示取整；

6e)按照公式<2>，对更新后的粒子计算粒子的新适应度值，若粒子的新适应度值大于粒子上一次迭代的适应度值，则用更新后的粒子代替当前粒子，并将每个粒子上的原子按原子编号从小到大排序，否则，粒子保持不变，同时，对更新粒子的速度按粒子上的原子排序后的索引重新排列，使得每个粒子的每维位置与其速度一一对应；

6f)用当前正在运动的粒子群作为待交叉种群，执行交叉操作：

若该类图像块为单方向类图像块，则设置交叉概率Pc＝0.5，若为多方向类图像块，则设置交叉概率Pc＝0.6，为每个粒子产生一个(0,1)之间的随机数，若该随机数小于等于Pc，则将该粒子选作待交叉粒子，否则，不选作待交叉粒子；

若该类图像块为单方向类图像块，则从其它组中随机挑选另一个粒子作为另一待交叉粒子，若该类图像块为多方向类图像块，则从待交叉种群中随机挑选另一个粒子作为另一待交叉粒子；

为两个待交叉粒子随机产生一个(1,k-1)的交叉点，将两个待交叉粒子及两个粒子的速度进行单点交叉，产生两个新的子代粒子；

从两个新的子代粒子中选择一个适应度值较大的粒子，将其适应度值与产生它的两个父代粒子的适应度值进行比较：若比父代粒子的适应度值大，则代替父代粒子，并将交叉后的每个粒子上的原子按原子编号从小到大进行排序，将粒子的速度按粒子上的原子排序后的索引重新排列；否则，不替代；

6g)更新粒子群中粒子的历史最优位置以及全局历史最优位置；

6h)判断是否满足粒子群算法的迭代停止条件，若满足，则执行步骤6i)，若不满足，返回步骤6b)，继续对粒子进行速度位置更新操作；

6i)将粒子群的全局历史最优解作为该类图像块的重构基。

步骤7，根据步骤5)、6)得到的所有观测向量的最优原子组合，计算所有图像块的估计值，并按顺序拼接成一整幅重构图像输出。

7a)通过如下公式计算计算所有图像块的估计值：

x_n＝D_n[(ΦD_n)⁺y_n]，

其中，D_n表示该图像块的最优原子组合，y_n表示该图像块的观测向量，Φ为高斯观测矩阵，(·)⁺表示计算矩阵的伪逆矩阵；

7b)将所有图像块的估计值按顺序拼接起来，得到重构图像并输出；

本发明的效果由以下仿真的图像和数据进一步说明：

1.仿真条件

1)仿真实验使用的是512×512的标准测试图像库中的Lena图、Barbara图，图像块大小定为16×16；

2)本实验的观测矩阵为随机高斯观测矩阵，采样率分别为40％、30％、20％；

3)本实验采用的Ridgelet冗余字典规模为12032，一共36个方向；

4)本实验光滑块的稀疏度值固定为20，单方向块和多方向块的稀疏度值定为32；

5)本实验光滑类重构时粒子群大小为20，单方向类重构时粒子群大小为20，多方向类重构时粒子群大小为36；

6)本实验光滑类重构时迭代20次，单方向类重构时迭代次数100次，多方向类重构时迭代次数100次；

7)本实验选择CPU为Interi5-3470，主频为3.2GHZ，内存为4G，操作系统为32位Win7，仿真平台为Matlab2012a。

2.仿真内容与结果

仿真1，在采样率为30％的条件下，用本发明方法和现有方法分别对Lena图进行重构，仿真结果如图2所示，其中：

图2(a)为Lena原图，图2(b)为图2(a)的局部放大图；

图2(c)为用结构稀疏方法得到的重构图，图2(d)为图2(c)的局部放大图；

图2(e)为用两阶段重构方法得到的重构图，图2(f)为图2(e)的局部放大图；

图2(g)为用本发明得到的重构图，图2(h)为图2(g)的局部放大图；

图2的实验结果说明，使用本发明方法得到的重构图像比使用结构稀疏方法和两阶段重构方法得到的重构图像，在视觉效果上更佳。从这些图像的局部放大图的对比可以看出，本发明在Lena肩膀的光滑部分显得更平滑，且边缘较清晰，本发明方法能够获得更清晰的图像的边缘，更平滑的光滑区域，说明本发明对自然图像具有较好的重构性能。

仿真2，在采样率为30％的条件下，用本发明方法和现有方法分别对Barbara图进行重构，仿真结果如图3所示，其中：

图3(a)为Barbara原图，图3(b)为图3(a)的局部放大图；

图3(c)为用结构稀疏方法得到的重构图，图3(d)为图3(c)的局部放大图；

图3(e)为用两阶段重构方法得到的重构图，图3(f)为图3(e)的局部放大图；

图3(g)为用本发明得到的重构图，图3(h)为图3(g)的局部放大图；

图3的实验结果说明，与已有的两种方法相比，使用本发明方法得到的重构图像，在视觉效果上更好。从各局部放大图的对比可以看出，本发明对Barbara裤腿上的纹理重构得更清晰，说明本发明对图像的边缘、单方向纹理有更准确的估计。

仿真3，在不同采样率下，用本发明方法和现有方法分别对Lena，Barbara和Boat图进行重构，并将得到的数值结果进行比较，结果如表1，表2和表3以及图4所示，其中：

图4是三种方法的平均运行时间对比图。

表1三种方法的图像峰值信噪比PSNR指标的结果

表2三种方法的图像结构相似性SSIM指标的结果

表3三种方法的平均运行时间Time指标的结果

表1和表2的实验结果说明，在三种实验方法中，本发明在绝大部分的测试采样率和测试图像上，能取得最好的高峰值信噪比PSNR和结构相似SSIM指标值，说明本发明对自然图像有较好的重构性能。

表3和图4的实验结果说明，本发明方法的运行时间远远少于已有的两种方法，更适合实时应用。

综上所述，本发明实现了通过充分利用图像的方向结构特性，采用速度较快的不同重构策略的粒子群算法进行重构，能够以较少的时间代价，获得对自然图像良好的压缩感知重构效果。

Claims (7)

1.基于图像结构的粒子群优化非凸压缩感知图像重构方法，包括如下步骤：

1)压缩感知接收方接收发送方发送过来的所有图像块的观测向量y₁，y₂,...,y_n，其中，n表示图像块的数目；

2)对每个观测向量所对应的图像块进行结构判别，将图像块标记为光滑块、单方向块和多方向块，并记录单方向块的方向；

3)分别对不同结构类型的图像块所对应的观测向量进行聚类，并根据脊波过完备字典为每一类观测向量所对应的图像块构造一个相应的过完备冗余字典；

4)根据步骤3)得到的过完备字典以及图像块的结构类型分别为每一类光滑类图像块、单方向类图像块和多方向类图像块构造相应的粒子群初始种群；

5)对每一类光滑类图像块，在其对应的过完备冗余字典下，利用基于分组初始化策略的搜索尺度最优原子组合，得到每一类光滑类图像块的最优原子组合；

6)对每一类单方向类和多方向类图像块，分别在其对应的过完备冗余字典下，利用基于交叉和原子方向约束的粒子群算法搜索粒子在方向和尺度上的最优原子组合，得到每一类图像块的最优原子组合；

7)根据步骤5)、6)得到的所有观测向量的最优原子组合，计算所有图像块的估计值，并按顺序拼接成一整幅重构图像输出。

2.根据权利要求1所述基于图像结构的粒子群优化非凸压缩感知图像重构方法，其中所述步骤2)中对每个观测向量所对应的图像块进行结构判别，按如下步骤进行：

2a)计算每个观测向量的方差，并设置光滑阈值为0.45σ，其中，σ是所有观测向量的方差的平均值；

2b)对每个观测向量，判断其方差是否小于等于光滑阈值，若是，则将该观测向量对应的图像块标记为光滑块，否则，不做标记；

2c)对每个尚未标记的图像块的观测向量，进行以下判定：

2c1)根据脊波函数生成脊波过完备冗余字典，该字典共有12032个原子，36个方向和16个尺度，每个尺度下大约有20个左右的位移；

2c2)将脊波过完备冗余字典按方向分为36个子字典ψ₁，ψ₂，...，ψ_i,...,ψ₃₆，对每一个图像块，根据其观测向量y和这些子字典，计算得到一个观测残差序列r₁，r₂，...，r_i,...,r₃₆，找到序列中的最小值在序列中的位置i，i＝1,2,...,36；

2c3)按如下公式计算序列中的子字典ψ_i的观测残差r_i：

式中，y是待判定图像块的观测向量，Φ是用于观测的高斯矩阵，D_r是子字典ψ_i中与y相关性最大的10个原子组合，(·)⁺表示计算矩阵的伪逆矩阵，是向量二范数的平方；

2c4)利用观测残差序列中位置i-2，i-1，i，i+1和i+2上的五个残差值r_i-2，r_i-1，r_i，r_i+1和r_i+2，对观测向量y对应的图像块进行标记：若r_i-2大于r_i-1，r_i-1大于1.2r_i，r_i+1大于1.2r_i，且r_i+2大于r_i+1，则将观测向量y对应的图像块标记为单方向，并记录该单方向块的方向为第i个方向，否则，不做标记；

其中，在判断条件中，若i为1，r_i-1和r_i-2分别用r₃₆和r₃₅代替，若i为2，r_i-2用r₃₆代替，若i为36，r_i+1和r_i+2分别用r₁和r₂代替，若i为35，r_i+2用r₁代替；

2c5)将其它尚未标记的图像块标记为多方向。

3.根据权利要求1所述基于图像结构的粒子群优化非凸压缩感知图像重构方法，其中所述步骤3)中为每一类观测向量所对应的图像块构造一个相应的过完备冗余字典，按如下步骤进行：

3a)为光滑类图像块构造过完备字典Ψ_s，该过完备字典由脊波冗余字典中包含所有方向的前5个尺度的脊波冗余子字典组成， 表示包含所有方向的尺度为h的脊波冗余子字典，其中，h＝1,2,..,5；依次将中的原子按原子编号从小到大映射到一个连续的正整数空间中，将作为光滑类图像块的过完备冗余字典Ψ_s的映射字典，其中，m是光滑类图像块的过完备冗余字典Ψ_s的大小；

3b)为单方向类图像块构造过完备冗余字典Ψ_g:从脊波冗余字典中取出第r个方向及与第r个方向左右相邻的4个方向的子字典作为该单方向类图像块的过完备冗余字典Ψ_g，其中，r表示该单方向类图像块的方向，表示第r个方向的脊波冗余子字典，r＝1,2,...,36，若r为1，则r-1和r-2分别用36和35代替，若r为2，则r-2用36代替，若r为36，则r+1和r+2分别用1和2代替，若r为35，则r+2用1代替；

3c)将整个脊波过完备冗余字典作为多方向类图像块的字典。

4.根据权利要求1所述的方法，其特征在于，步骤4)所述的为每一类光滑类图像块、单方向类图像块和多方向类图像块构造相应的粒子群初始种群，按如下步骤进行：

4a)对每一光滑类图像块，根据光滑类图像块的过完备冗余字典Ψ_s的映射字典按照该映射字典的5个子字典，初始化一个规模为20的粒子群A₁，并将该粒子群A₁均分成5个组,每个组一个尺度；随机初始化每个组的粒子方向，并将每个组中粒子的速度随机初始化为[-L^h,L^h]×0.1中的任意一个值，其中L^h为映射字典中的子字典的大小，h＝1,2,..,5；

4b)对每一单方向类图像块，根据其过完备冗余字典Ψ_g的5个子字典，初始化一个规模为20的粒子群A₂，并将该粒子群A₂均分成5个组，每个组代表一个方向，随机初始化每个组中的粒子，并将每个粒子上的原子按原子编号从小到大有序排列，再将粒子的速度随机初始化为[-L^r,L^r]×0.5中的任意一个值，其中L^r为该单方向类过完备冗余字典Ψ_g中的子字典的大小，r＝1,2,...,36；

4c)对多方向类图像块，初始化一个规模为36的粒子群A₃，每个粒子代表两个方向，其中一个方向从脊波过完备冗余字典中选择与该粒子编号相同的方向，另一个方向则从剩下的方向中随机选择一个，并将每个粒子上的原子也按原子编号从小到大有序排列，再将粒子的速度随机初始化为[-Lⁱ,Lⁱ]×0.5中的任意一个值，其中Lⁱ表示所选择的两个代表方向的脊波过完备冗余子字典的大小，i表示粒子编号，i＝1,2,...,36；

4d)对光滑类图像块的粒子群A₁，将其中的粒子映射回光滑类过完备冗余字典Ψ_s中；

4e)按照如下公式，计算光滑类、单方向类和多方向类粒子群中每个粒子的适应度值：

其中f(b_i)为该类观测向量对应的图像块的粒子群中第i个粒子的适应度值，b_i是粒子群中第i个粒子所代表的解，对应着脊波过完备冗余字典中的一组基原子组合，即一个子字典，j为观测向量的标号，q为聚类后该类的总测量数，y_j为类中的第j个测量，Φ为高斯观测矩阵，α_j为该类观测向量中第j个测量对应的图像块的稀疏表示系数，它是由传感矩阵的广义逆矩阵和该观测向量y_j相乘得到，传感矩阵由高斯观测矩阵Φ和子字典b_i相乘得到，是向量二范数的平方；

4f)初始化粒子群的历史最优解以及全局历史最优解：将粒子群中每个粒子的历史最优解初始化为粒子自身，找出整个粒子群中适应度值最大的粒子，将其作为全局历史最优解。

5.根据权利要求1所述的方法，其特征在于，所述步骤5)中对每一类光滑类图像块，在其对应的过完备冗余字典下，利用基于分组初始化策略的搜索尺度最优原子组合，按如下步骤进行：

5a)对粒子群算法中的参数进行设置，将粒子群的粒子速度位置更新为：

表示第t+1次迭代更新后粒子b_i第j个分量对应的速度，表示第t+1次迭代更新后粒子b_i第j个分量的取值，b_i表示第i个粒子，v_i表示第i个粒子b_i对应的速度，b_i与v_i具有相同的维数，两者一一对应；表示粒子群中第i个粒子b_i的第j个分量，表示粒子群中第i个粒子对应的速度v_i的第j个分量，表示第i个粒子的历史最优位置的第j个分量，表示全局历史最优位置的第j个分量，t表示当前是第几次迭代；w表示惯性权重，对w采取了线性递减策略，w从0.9线性递减至0.4；c₁表示粒子跟踪自己本身历史最优值的权重系数，c₁从2.5线性递减至0.5；c₂表示粒子跟踪群体最优值的权重系数，c₂从0.5线性递增至2.5；ξ、η是[0,1]区间内均匀分布的两个随机数；

5b)根据粒子群的粒子速度位置更新公式计算粒子的新速度，对得到的粒子的新速度进行控制，若则令若则令然后再计算粒子的新位置，对得到的粒子的新位置进行取整操作，使其映射在光滑类图像块的过完备冗余字典的映射字典中，其中，m是映射字典的大小；

5c)计算更新后的粒子的适应度值，将更新后粒子新的适应度值与粒子自身的历史最优值进行比较，若大于历史最优值，则用更新后的粒子代替历史最优解，并从粒子群中所有粒子的历史最优解中找出适应度值最大的历史最优解，若其适应度值大于全局历史最优值，则用该粒子的历史最优解代替全局历史最优解；

5d)判断是否满足粒子群算法的迭代停止条件，若满足，则执行步骤5e)，若不满足，返回步骤5b)，继续对粒子进行速度位置更新操作；

5e)将粒子群的全局历史最优解映射回光滑类图像块的过完备冗余字典Ψ_s中，得到真实的全局历史最优解，将该全局历史最优解原子组合作为该光滑类图像块的重构基。

6.根据权利要求1所述的方法，其特征在于，所述步骤6)中对每一类单方向类和多方向类图像块，分别在其对应的过完备冗余字典下，利用基于交叉和原子方向约束的粒子群算法搜索粒子在方向和尺度上的最优原子组合，按如下步骤进行：

6a)对粒子群算法中的参数进行设置：惯性权重w从0.9线性递减至0.4，粒子跟踪自己本身历史最优值的权重系数c₁从2.5线性递减至0.5，粒子跟踪群体最优值的权重系数c₂从0.5线性递增至2.5，ξ、η分别为[0,1]区间内均匀分布的两个随机数，增设一个控制粒子速度的参数P_r，P_r随着进化代数从1.0线性递减至0.1；

6b)用粒子速度位置更新公式计算粒子的速度和位置，计算第t次迭代每个粒子每维分量上的原子所在的方向r，得到第r个方向的原子范围[g^r,h^r]；

6c)对更新后粒子的新速度进行控制，约束粒子的方向：若则控制若则控制即粒子的新位置被约束在第r个方向的原子范围[g^r,h^r]内；

6d)生成一个大小为n×k的随机矩阵R＝{r₁,r₂,..,r_i,...,r_n}，r_i＝(r_i1,r_i2,...,r_ij,...,r_ik)，按照下式，计算粒子的新位置：

其中，表示第t+1次迭代更新后的粒子b_i第j个分量，表示第t次迭代粒子群中第i个粒子b_i的第j个分量，表示对应的速度，i＝1,2,...,n，j＝1,2,...,k，n为粒子群的大小，k为粒子位置的编码长度，即稀疏度，[·]表示取整；

6e)按照公式<2>，对更新后的粒子计算粒子的新适应度值，若粒子的新适应度值大于粒子上一次迭代的适应度值，则用更新后的粒子代替当前粒子，并将每个粒子上的原子按原子编号从小到大排序，否则，粒子保持不变，同时，对更新粒子的速度按粒子上的原子排序后的索引重新排列，使得每个粒子的每维位置与其速度一一对应；

6f)用当前正在运动的粒子群作为待交叉种群，执行交叉操作：

若该类图像块为单方向类图像块，则设置交叉概率Pc＝0.5，若为多方向类图像块，则设置交叉概率Pc＝0.6，为每个粒子产生一个(0,1)之间的随机数，若该随机数小于等于Pc，则将该粒子选作待交叉粒子，否则，不选作待交叉粒子；

若该类图像块为单方向类图像块，则从其它组中随机挑选另一个粒子作为另一待交叉粒子，若该类图像块为多方向类图像块，则从待交叉种群中随机挑选另一个粒子作为另一待交叉粒子；

为两个待交叉粒子随机产生一个(1,k-1)的交叉点，将两个待交叉粒子及两个粒子的速度进行单点交叉，产生两个新的子代粒子；

从两个新的子代粒子中选择一个适应度值较大的粒子，将其适应度值与产生它的两个父代粒子的适应度值进行比较：若比父代粒子的适应度值大，则代替父代粒子，并将交叉后的每个粒子上的原子按原子编号从小到大进行排序，将粒子的速度按粒子上的原子排序后的索引重新排列；否则，不替代；

6g)更新粒子群中粒子的历史最优位置以及全局历史最优位置；

6h)判断是否满足粒子群算法的迭代停止条件，若满足，则执行步骤6i)，若不满足，返回步骤6b)，继续对粒子进行速度位置更新操作；

6i)将粒子群的全局历史最优解作为该类图像块的重构基。

7.根据权利要求1所述的方法，其特征在于，所述步骤7)中计算所有图像块的估计值，通过如下公式计算：

x_n＝D_n[(ΦD_n)⁺y_n]，<5>其中，D_n表示该图像块的最优原子组合，y_n表示该图像块的观测向量，Φ为高斯观测矩阵，(·)⁺表示计算矩阵的伪逆矩阵。