CN103761755B - 基于进化多目标优化的非凸压缩感知图像重构方法 - Google Patents

Abstract

本发明涉及基于进化多目标优化的非凸压缩感知图像重构方法，用于对医学图像和自然图像进行重构。本发明提出了分块打散压缩采样方法，克服了现有分块压缩感知重构技术中，有些图像块不稀疏或者弱稀疏的缺点，提高了重构图像的质量。本发明将压缩感知重构中有约束的单目标优化问题转化为将稀疏度也作为优化目标的多目标优化问题，并针对图像利用进化多目标优化思想实现了基于小波域下的压缩感知优化重构方法，克服了现有压缩感知重构技术中，难以确定稀疏度的不足，提高了重构图像的质量。本发明将图像的边缘位置信息作为指导IHT方法求解稀疏系数的位置先验，克服了现有压缩感知重构技术中，较少考虑稀疏系数位置的不足，从而提高了重构图像的质量。

Description

基于进化多目标优化的非凸压缩感知图像重构方法

技术领域

本发明属于图像处理技术领域，确切讲是涉及基于进化多目标优化的非凸压缩感知图像重构方法，用于对医学图像和自然图像进行重构。

背景技术

近几年，在信号处理领域出现了一种新的数据采集理论“压缩感知”CS，该理论在数据采集的同时实现压缩，突破了传统奈奎斯特采样定理的限制，为数据采集技术带来了革命性的变化，使得该理论在压缩成像系统、军事密码学、无线传感等领域有着广阔的应用前景。压缩感知理论主要包括信号的稀疏表示、信号的观测和信号的重构等三个方面。在信号重构方面，通过求解l₀或l₁范数的优化问题来重构图像。

Tropp等人在文献“JoelA.Tropp,AnnaC.Gilbert,SignalRecoveryFromRandomMeasurementsViaOrthogonalMatchingPursuit”中提出基于正交匹配追踪的随机观测的信号恢复方法。该方法对稀疏信号进行低采样的随机观测，从正交的原子库中选择最能匹配信号结构的原子，从而重构出图像。该方法存在的不足是，在重构过程中使用贪婪思想寻找最优解，并不能保证收敛到全局最优解，从而导致重构出的图像不够准确，并且它需要人为指定稀疏度的大小，并不是自适应地寻找适合特定问题的稀疏度，从而导致重构图像质量不好，另外它对压缩感知框架强加了有限等距性RIP约束，从某种意义上讲，限制了压缩感知的应用范围。

HuiLi等人在文献“MOEA/DwithIterativeThresholdingAlgorithmforSparseOptimizationProblems”中提出基于分解的多目标进化算法MOEA/D用于稀疏信号恢复的方法。该方法针对一维随机稀疏信号，将压缩感知重构中有约束的单目标优化问题转化为将稀疏度也作为优化目标的多目标优化问题，且通过结合MOEA/D和迭代硬阈值IHT进行压缩感知优化重构。该方法的不足是，它重构的是一维随机稀疏信号，并没有针对图像来实现该方法，自然就没有引入图像的先验信息来指导，并且IHT也没有使用位置信息来指导求解，其稀疏度范围也是人为定的，从而导致重构的图像质量不好。

另外，对于小波域下的分块压缩感知图像重构方法，其优点是：运算量小，耗时小；其不足之处是：在小波域下，将高频小波系数进行分块后，可能会出现某些块不稀疏或者弱稀疏，从而违背了压缩感知的理论基础，导致重构的图像质量不好。

发明内容

本发明的目的是针对小波域下的分块压缩感知图像重构方法中，可能会出现图像块不稀疏或者弱稀疏的缺点，提供一种基于进化多目标优化的非凸压缩感知图像重构方法，以便优化图像重构算法，提高图像准确重构质量。

实现本发明的技术方案是：基于进化多目标优化的非凸压缩感知图像重构方法，其特征是：至少包括如下步骤：

步骤101：输入测试图像，对测试图像进行小波变换，保留低频小波系数C_l，利用正交随机高斯观测矩阵Φ对高频小波系数每隔8个像素进行分块打散压缩采样，得到每一块的观测向量y；

步骤102：利用保留的低频小波系数C_l，通过小波逆变换、边缘检测和分块打散操作得到对应于高频小波系数的位置矩阵E，位置矩阵E的每一列e对应于每一块的观测向量y；

步骤103：利用位置矩阵E通过统计方法获取每一块的块稀疏度范围[n,n+120]，其中n是位置矩阵E中对应于观测向量y的那一列的非零元素个数；

步骤104：对位置矩阵E执行打散块合并操作，得到分块打散前的位置矩阵E′；

步骤105：对分块打散前的位置矩阵E′进行提取疫苗和注射疫苗操作，得到分块打散前的过渡位置矩阵E″；

步骤106：对分块打散前的过渡位置矩阵E″执行分块打散操作，得到新的位置矩阵新的位置矩阵的每一列对应于每一块的观测向量y；

步骤107：设置种群大小为10，当前进化代数为p，种群中每个个体由编码对（目标1、目标2）组成，目标1对应于要重构的块的小波系数，目标2表示对应的块稀疏度，利用正交随机高斯观测矩阵Φ和每一块的观测向量y通过广义逆，并结合新的位置矩阵中的每一列对种群中每个个体的目标1进行基于块位置信息的初始化操作和根据对应于每一块的观测向量y的块稀疏度范围[n,n+120]对种群中每个个体的目标2进行初始化，得到每一块的种群P＝{(w₁′,s₁),...,(w₁₀′,s₁₀)}，其中w_j′是种群P中的第j个个体的目标1，即对应的块的小波系数，s_j是种群P中的第j个个体的目标2，即对应的块稀疏度，j＝1,…,10；

步骤108：对每一块的种群P中每个个体的目标1根据新的位置矩阵及其目标2进行一致变异操作，得到每一块的子代种群其中是子代种群中第j个个体的目标1，s_j是子代种群中的第j个个体的目标2，j＝1,…,10；

步骤109：对每一块的子代种群中每个个体的目标1根据对应于每一块的新的位置矩阵中的一列及其目标2进行迭代硬阈值局部搜索，得到每一块的新子代种群其中是新子代种群中第j个个体的目标1，s_j是新子代种群中的第j个个体的目标2，j＝1,…,10，并计算每一块的新子代种群中每个个体的适应度；

步骤110：从由每一块的新子代种群中每个个体的目标1组成的集合中进行非支配解选择，并将非支配解对应地保存在每一块的非支配解集合Q中，但不保存弱支配解；

步骤111：根据每一块的非支配解集合Q中的元素个数更新新子代种群中相应个体的目标2，即稀疏度，若新子代种群中相应个体的目标2有更新，则更新每一块的新子代种群中更新了目标2的个体的目标1；若新子代种群中相应个体的目标2没有更新，则保留新子代种群中的个体，在这里为了方便，不管新子代种群中相应个体的目标2有没有更新，都将此操作后的个体记为于是得到每一块的更新子代种群其中是更新子代种群中第j个个体的目标1，s_j′是更新子代种群中第j个个体的目标2，j＝1,…,10，并计算每一块的更新子代种群中每个个体的适应度；

步骤112：如果当前进化代数p大于设置的最大进化代数mp，或者适应度达到要求，则执行步骤113，并保留最终得到的每一块的更新子代种群中适应度最大的个体作为每一块的最优小波系数解；否则p=p+1，并将每一块的更新子代种群覆盖每一块的种群P，返回步骤108；

步骤113：对所有块的最优小波系数解执行打散块合并操作，即恢复成分块打散压缩采样前的形式，结合保留的低频小波系数C_l进行小波逆变换得到测试图像的重构图像。

所述步骤101中的利用正交随机高斯观测矩阵Φ对高频小波系数每隔8个像素进行分块打散压缩采样，包括如下步骤：

步骤201：将高频小波系数分成8×8的块；

步骤202：按从左到右，从上到下的顺序分别取每个8×8的块的一个元素(也是按照从左到右，从上到下的顺序取)组成打散块信号，这样可得到64个打散块信号；

步骤203：利用正交随机高斯观测矩阵Φ对64个打散块信号分别观测。

所述步骤102中的通过小波逆变换、边缘检测和分块打散操作得到对应于高频小波系数的位置矩阵E，包括如下步骤：

步骤301：将高频小波系数设定为零，再利用保留的低频小波系数C_l，通过小波逆变换得到一个边缘模糊的图像I₀；

步骤302：利用Canny边缘检测方法对边缘模糊的图像I₀进行边缘检测，并提取边缘图像得到边缘图I₁；

步骤303：对边缘图I₁进行一层小波变换，得到边缘图I₁的高频小波系数矩阵B；

步骤304：对边缘图I₁的高频小波系数矩阵B每隔8个像素进行分块打散操作得到对应的分块打散操作后的边缘图I₁的高频小波系数矩阵其中的每一列对应于每一块的观测向量y；

步骤305：分别对分块打散操作后的边缘图I₁的高频小波系数矩阵的每一列求均值，并将均值的0.6倍作为阈值对的每一列进行阈值处理，将每一列中大于0.6倍均值的元素置为1；否则，置为0；得到阈值处理后的对应于分块打散操作后的边缘图I₁的高频小波系数矩阵的位置矩阵E，其中位置矩阵E的每一列e对应于每一块的观测向量y。

所述步骤102中的分块打散操作，包括如下步骤：

步骤401：将边缘图I₁的高频小波系数矩阵B分成8×8的块；

步骤402：按从左到右，从上到下的顺序分别取每个8×8的块的一个元素(也是按照从左到右，从上到下的顺序取)组成打散块信号，这样可得到由64个打散块信号组成的矩阵。

所述步骤103中的利用位置矩阵E通过统计方法获取每一块的块稀疏度范围[n,n+120]，包括如下步骤：

步骤501：统计位置矩阵E中每一列的非零元素个数n；

步骤502：将每一列的非零元素个数n加上120，得到对应于每一块的观测向量y的块稀疏度范围[n,n+120]。

所述步骤104中的对位置矩阵E执行打散块合并操作，包括如下步骤：

步骤601：按从左到右的顺序获取位置矩阵E中的每一列；

步骤602：将位置矩阵E中每一列的元素按从左到右、从上到下的顺序每隔8个像素放置在分块打散前的位置上，得到分块打散前的位置矩阵E′。

所述步骤105中的对分块打散前的位置矩阵E′进行提取疫苗和注射疫苗操作，包括如下步骤：

步骤701：在分块打散采样前的位置矩阵E′上用3×3的窗口进行滑动，考虑3×3的窗口的上、下、左、右四个相邻的值，如果这四个值为1的个数大于3，则将1作为该中心位的值；如果四个值为0的个数为4，则将0作为该中心位置的值；若两种情况都不满足，则中心位置上的值不变。

所述步骤107中的结合新的位置矩阵中的每一列对种群中每个个体的目标1进行基于块位置信息的初始化操作和根据对应于每一块的观测向量y的块稀疏度范围[n,n+120]对种群中每个个体的目标2进行初始化，包括如下步骤：

步骤801：设置种群大小为10，种群中每个个体由编码对(目标1，目标2)组成，目标1对应于要重构的块的小波系数，目标2表示对应的块稀疏度，将种群记为P＝{(w₁′,s₁),...,(w₁₀′,s₁₀)}，其中w_j′是种群P中的第j个个体的目标1，即对应的块的小波系数，且每一个w_j′都采用实数编码，s_j是种群P中的第j个个体的目标2，即对应的块稀疏度，且每一个s_j都采用整数编码，j＝1,…,10；

步骤802：根据正交随机高斯观测矩阵Φ和每一块的观测向量y通过广义逆求得初始的小波系数，公式如下：

${\tilde{w}}^0={\mathrm{\Φ}}^{+}y;$

步骤803：我们结合广义逆值和新的位置矩阵中对应于观测向量y的一列来初始化种群P中每个个体的目标1，则种群P中每个个体的目标1(j＝1,…,10)可按下面的公式进行初始化：

${w_j}^{\′}\left(t\right)\left\{\begin{array}{cc}{\tilde{w}}^0\left(t\right),& \mathrm{if}\tilde{e}\left(t\right)=0\\ L\×{\tilde{w}}^0\left(t\right),& \mathrm{if}\tilde{e}\left(t\right)=1\end{array}\right.$

其中，L是伸缩因子，且1≤L≤1.5，t＝1,…,1024；

步骤804：对种群P中每个个体的目标2(块稀疏度)进行初始化，每个个体的目标1w_j′对应的目标2s_j都是从对应于每一块的观测向量y的块稀疏度范围[n,n+120]中随机选取，其中n是位置矩阵E中对应于观测向量y的那一列的非零元素个数，且每个个体的目标2s_j都不同。

所述步骤108中的一致变异操作，包括如下步骤：

步骤901：统计种群P中第j个个体的目标1w_j′的非零元素个数作为第j个个体的目标1w_j′的真实块稀疏度k_j，j＝1,…,10；

步骤902：如果第j个个体的目标1w_j′的真实块稀疏度k_j(j＝1,…,10)小于第j个个体的目标2s_j(j＝1,…,10)，则在第j个个体的目标1w_j′的零元素中随机选取s_j-k_j个位置，并用一次迭代硬阈值(IHT)求得的值来替换这些位置上的零元素；如果第j个个体的目标1w_j′的真实块稀疏度k_j(j＝1,…,10)大于第j个个体的目标2s_j，则在第j个个体的目标1w_j′的非零元素中选取k_j-s_j个值最小的位置，并用0来替换这些位置上的值。

所述步骤109中的迭代硬阈值局部搜索，包括如下步骤：

步骤1001：对于子代种群中第j个个体的目标1(j＝1,…,10)，将其赋值给即 ${\hat{w}}_j^{{\left(1\right)}^{\′}}\←{{\hat{w}}_j}^{\′};$

步骤1002：根据迭代硬阈值的公式 ${\hat{w}}_j^{{(n+1)}^{\′}}=H_{{\tilde{e}}_i,s_j}({\hat{w}}_j^{{\left(n\right)}^{\′}}+{\mathrm{\Φ}}^T(y-\mathrm{\Φ}{\hat{w}}_j^{{\left(n\right)}^{\′}}))$ 计算20次，并将赋值给即其中，是用位置和第j个个体的目标2s_j指导的非线性算子。

所述步骤109中计算每一块的新子代种群中每个个体的适应度，是按照如下公式计算每一个个体的适应度

$f\left({{\hat{w}}_j}^{\′\′}\right)=1/{\left|\right|y-\mathrm{\Φ}{{\hat{w}}_j}^{\′\′}\left|\right|}_2^2$

其中，Φ是正交随机高斯观测矩阵，y是每一个个体对应的块的观测向量，表示向量2范数的平方。

所述步骤111中的根据每一块的非支配解集合Q中的元素个数更新新子代种群中相应个体的目标2，包括如下步骤：

步骤1101：将所有个体的目标2s_j(j＝1,…,10)按升序排序，即其中j₁,…,j₁₀∈{1,…,10}；

步骤1102：根据下面的公式确定待选块稀疏度集

步骤1103：如果从中选择一个值来更新块稀疏度，

如果|Q|＞0.5×10，其中r是从中随机选择的一个值；

如果|Q|≤0.5×10，其中r是从中随机选择的一个值；

步骤1104：如果且

如果|Q|＝10， $s_{j_1}\←\min \{s_{j_{10}}+0.5\×10,n+120\},$ 其中n+120是块稀疏度范围[n,n+120]的最大值；

如果|Q|＝1，其中n是块稀疏度范围[n,n+120]的最小值。

所述步骤111中的更新每一块的新子代种群中更新了目标2的个体的目标1，按如下过程进行：

若新子代种群中相应个体的目标2有更新，则选取每一块的新子代种群中与更新了目标2的个体(j_z∈{1,…,10})的目标2最相近的个体(j_τ∈{1,…,10})，如果选择的个体的目标1的残差比更新了目标2的个体的目标1的残差大，则不更新更新了目标2的个体的目标1；如果选择的个体的目标1的残差比更新了目标2的个体的目标1的残差小，则用替换来更新更新了目标2的个体的目标1。

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

第一，本发明提出了分块打散压缩采样方法，克服了现有分块压缩感知重构技术中，有些图像块不稀疏或者弱稀疏的缺点，提高了重构图像的质量。

第二，本发明将压缩感知重构中有约束的单目标优化问题转化为将稀疏度也作为优化目标的多目标优化问题，并针对图像利用进化多目标优化思想实现了基于小波域下的压缩感知优化重构方法，克服了现有压缩感知重构技术中，难以确定稀疏度的不足，提高了重构图像的质量。

第三，本发明利用图像的先验信息，克服了现有压缩感知重构技术中，需要人为指定稀疏度范围的不足，得到了较好的图像重构效果。

第四，本发明将图像的边缘位置信息作为指导IHT方法求解稀疏系数的位置先验，克服了现有压缩感知重构技术中，较少考虑稀疏系数位置的不足，从而提高了重构图像的质量。

第五，本发明从多目标进化算法和稀疏优化问题的特性出发，将进化多目标优化方法和IHT以及图像先验信息有机结合作为非凸压缩感知图像重构方法，得到了较好的图像重构效果。

下面结合附图对本发明实施例做进一步描述。

附图说明

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

图2是本发明分块打散压缩采样获得总体观测向量的子流程图；

图3是本发明分块打散操作的示意图；

图4是用本发明所采用的打散分块策略与现有压缩感知技术中的分块策略对StomachCT图像的比较图；

图5是用本发明所采用的打散分块策略与现有压缩感知技术中的分块策略对Lena图像的比较图；

图6是用本发明与现有技术在采样率为40%时对StomachCT图像重构结果图及局部放大图；

图7是用本发明与现有技术(OMP,IHT)重构出来的StomachCT图像的峰值信噪比PSNR随采样率变化的趋势图；

图8是用本发明与现有技术在采样率为40%时对Lena图像重构结果图及局部放大图；

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

具体实施方式

如图1所示，基于进化多目标优化的非凸压缩感知图像重构方法，至少包括如下步骤：

步骤101：获得总观测向量。

如图2所示，本步骤的具体实现如下：

1.1)输入测试图像，对测试图像进行小波变换；

1.2)保留低频小波系数C_l，对高频小波系数每隔8个像素进行分块打散压缩采样，得到打散块信号：

1.2.1)将高频小波系数分成8×8的块；

1.2.2)按从左到右，从上到下的顺序分别取每个8×8的块的一个元素(也是按照从左到右，从上到下的顺序取)组成打散块信号，这样可得到多个打散块信号；

1.3)利用正交随机高斯观测矩阵Φ对打散块信号分别进行观测得到每一块对应的观测向量y；

1.4)保留每一块的观测向量y和低频小波系数作为总体观测向量。

步骤102：利用保留的低频小波系数C_l，通过小波逆变换、边缘检测和分块打散操作得到对应于高频小波系数的位置矩阵E，位置矩阵E的每一列e对应于每一块的观测向量y。

2.1)将高频小波系数设定为零，再利用保留的低频小波系数C_l，通过小波逆变换得到一个边缘模糊的图像I₀；

2.2)利用Canny边缘检测方法对边缘模糊的图像I₀进行边缘检测，并提取边缘图像得到边缘图I₁；

2.3)对边缘图I₁进行一层小波变换，得到边缘图I₁的高频小波系数矩阵B；

2.4)对边缘图I₁的高频小波系数矩阵B每隔8个像素进行分块打散操作(如图3所示)得到对应的分块打散操作后的边缘图I₁的高频小波系数矩阵其中的每一列对应于每一块的观测向量y；

2.5)分别对分块打散操作后的边缘图I₁的高频小波系数矩阵的每一列求均值，并将均值的0.6倍作为阈值对的每一列进行阈值处理，将每一列中大于0.6倍均值的元素置为1；否则，置为0；得到阈值处理后的对应于分块打散操作后的边缘图I₁的高频小波系数矩阵的位置矩阵E，其中位置矩阵E的每一列e对应于每一块的观测向量y。

步骤103：利用位置矩阵E通过统计方法获取每一块的块稀疏度范围[n,n+120]，其中n是位置矩阵E中对应于观测向量y的那一列的非零元素个数。

3.1)统计位置矩阵E中每一列的非零元素个数n；

3.2)将每一列的非零元素个数n加上120，得到对应于每一块的观测向量y的块稀疏度范围[n,n+120]。

步骤104：对位置矩阵E执行打散块合并操作，得到分块打散前的位置矩阵E′。

如图3所示，本步骤的具体实现如下：

4.1)按从左到右的顺序获取位置矩阵E中的每一列；

4.2)将位置矩阵E中每一列的元素按从左到右、从上到下的顺序每隔8个像素放置在分块打散前的位置上，得到分块打散前的位置矩阵E′。

步骤105：对分块打散前的位置矩阵E′进行提取疫苗和注射疫苗操作，得到分块打散前的过渡位置矩阵E″。

在分块打散采样前的位置矩阵E′上用3×3的窗口进行滑动，考虑3×3的窗口的上、下、左、右四个相邻的值，如果这四个值为1的个数大于3，则将1作为该中心位的值；如果四个值为0的个数为4，则将0作为该中心位置的值；若两种情况都不满足，则中心位置上的值不变。由此，得到分块打散前的过渡位置矩阵E″。

步骤106：对分块打散前的过渡位置矩阵E″执行分块打散操作，得到新的位置矩阵新的位置矩阵的每一列对应于每一块的观测向量y。

6.1)将分块打散前的过渡位置矩阵E″分成8×8的块；

6.2)按从左到右，从上到下的顺序分别取每个8×8的块的一个元素(也是按照从左到右，从上到下的顺序取)组成打散块信号，这样可得到由64个打散块信号组成的矩阵，即是新的位置矩阵

步骤107：设置种群大小为10，当前进化代数为p，种群中每个个体由编码对(目标1，目标2)组成，目标1对应于要重构的块的小波系数，目标2表示对应的块稀疏度，利用正交随机高斯观测矩阵Φ和每一块的观测向量y通过广义逆，并结合新的位置矩阵中的每一列对种群中每个个体的目标1进行基于块位置信息的初始化操作和根据对应于每一块的观测向量y的块稀疏度范围[n,n+120]对种群中每个个体的目标2进行初始化，得到每一块的种群P＝{(w₁′,s₁),...,(w₁₀′,s₁₀)}，其中w_j′是种群P中的第j个个体的目标1，即对应的块的小波系数，s_j是种群P中的第j个个体的目标2，即对应的块稀疏度，j＝1,…,10。

7.1)设置种群大小为10，种群中每个个体由编码对(目标1，目标2)组成，目标1对应于要重构的块的小波系数，目标2表示对应的块稀疏度，将种群记为P＝{(w₁′,s₁),...,(w₁₀′,s₁₀)}，其中w_j′是种群P中的第j个个体的目标1，即对应的块的小波系数，且每一个w_j′都采用实数编码，s_j是种群P中的第j个个体的目标2，即对应的块稀疏度，且每一个s_j都采用整数编码，j＝1,…,10；

7.2)根据正交随机高斯观测矩阵Φ和每一块的观测向量y通过广义逆求得初始的小波系数，公式如下：

${\tilde{w}}^0={\mathrm{\Φ}}^{+}y;$

7.3)我们结合广义逆值和新的位置矩阵中对应于观测向量y的一列来初始化种群P中每个个体的目标1，则种群P中每个个体的目标1w_j′(j＝1,…,10)可按下面的公式进行初始化：

${w_j}^{\′}\left(t\right)\left\{\begin{array}{cc}{\tilde{w}}^0\left(t\right),& \mathrm{if}\tilde{e}\left(t\right)=0\\ L\×{\tilde{w}}^0\left(t\right),& \mathrm{if}\tilde{e}\left(t\right)=1\end{array}\right.$

其中，L是伸缩因子，且1≤L≤1.5，t＝1,…,1024；

7.4)对种群P中每个个体的目标2(块稀疏度)进行初始化，每个个体的目标1w_j′对应的目标2s_j都是从对应于每一块的观测向量y的块稀疏度范围[n,n+120]中随机选取，其中n是位置矩阵E中对应于观测向量y的那一列的非零元素个数，且每个个体的目标2s_j都不同。

步骤108：对每一块的种群P中每个个体的目标1根据新的位置矩阵及其目标2进行一致变异操作，得到每一块的子代种群其中是子代种群中第j个个体的目标1，s_j是子代种群中的第j个个体的目标2，j＝1,…,10。

8.1)统计种群P中第j个个体的目标1w_j′的非零元素个数作为第j个个体的目标1w_j′的真实块稀疏度k_j，j＝1,…,10；

8.2)如果第j个个体的目标1w_j′的真实块稀疏度k_j(j＝1,…,10)小于第j个个体的目标2s_j(j＝1,…,10)，则在第j个个体的目标1w_j′的零元素中随机选取s_j-k_j个位置，并用一次迭代硬阈值(IHT)求得的值来替换这些位置上的零元素；如果第j个个体的目标1w_j′的真实块稀疏度k_j(j＝1,…,10)大于第j个个体的目标2s_j，则在第j个个体的目标1w_j′的非零元素中选取k_j-s_j个值最小的位置，并用0来替换这些位置上的值；由此，得到每一块的子代种群其中是子代种群中第j个个体的目标1，s_j是子代种群中的第j个个体的目标2，j＝1,…,10。

步骤109：对每一块的子代种群中每个个体的目标1根据对应于每一块的新的位置矩阵中的一列及其目标2进行迭代硬阈值局部搜索，得到每一块的新子代种群其中是新子代种群中第j个个体的目标1，sj是新子代种群中的第j个个体的目标2，j＝1,…,10，并计算每一块的新子代种群中每个个体的适应度。

9.1)对于子代种群中第j个个体的目标1(j＝1,…,10)，将其赋值给即 ${\hat{w}}_j^{{\left(1\right)}^{\′}}\←{{\hat{w}}_j}^{\′};$

9.2)根据迭代硬阈值的公式 ${\hat{w}}_j^{{(n+1)}^{\′}}=H_{{\tilde{e}}_i,s_j}({\hat{w}}_j^{{\left(n\right)}^{\′}}+{\mathrm{\Φ}}^T(y-\mathrm{\Φ}{\hat{w}}_j^{{\left(n\right)}^{\′}}))$ 计算20次，并将赋值给即其中，是用位置和第j个个体的目标2s_j指导的非线性算子；

9.3)按照如下公式计算每一块的新子代种群中的每一个个体的适应度

$f\left({{\hat{w}}_j}^{\′\′}\right)=1/{\left|\right|y-\mathrm{\Φ}{{\hat{w}}_j}^{\′\′}\left|\right|}_2^2$

其中，Φ是正交随机高斯观测矩阵，y是每一个个体对应的块的观测向量，表示向量2范数的平方。

步骤110：从由每一块的新子代种群中每个个体的目标1组成的集合中进行非支配解选择，并将非支配解对应地保存在每一块的非支配解集合Q中，但不保存弱支配解。

步骤111：根据每一块的非支配解集合Q中的元素个数更新新子代种群中相应个体的目标2，即稀疏度，若新子代种群中相应个体的目标2有更新，则更新每一块的新子代种群中更新了目标2的个体的目标1；若新子代种群中相应个体的目标2没有更新，则保留新子代种群中的个体，在这里为了方便，不管新子代种群中相应个体的目标2有没有更新，都将此操作后的个体记为于是得到每一块的更新子代种群其中是更新子代种群中第j个个体的目标1，s_j′是更新子代种群中第j个个体的目标2，j＝1,…,10，并计算每一块的更新子代种群中每个个体的适应度。

11.1)将所有个体的目标2s_j(j＝1,…,10)按升序排序，即其中j₁,…,j₁₀∈{1,…,10}；

11.2)根据下面的公式确定待选块稀疏度集

11.3)如果从中选择一个值来更新块稀疏度，

11.3-1)如果|Q|＞0.5×10，其中r是从中随机选择的一个值；

11.3-2)如果|Q|≤0.5×10，其中r是从中随机选择的一个值；

11.4)如果

11.4-1)如果|Q|＝10， $s_{j_1}\←\min \{s_{j_{10}}+0.5\×10,n+120\},$ 其中n+120是块稀疏度范围[n,n+120]的最大值；

11.4-2)如果|Q|＝1，其中n是块稀疏度范围[n,n+120]的最小值；

11.5)若新子代种群中相应个体的目标2有更新，则选取每一块的新子代种群中与更新了目标2的个体(j_z∈{1,…,10})的目标2最相近的个体(j_τ∈{1,…,10})，如果选择的个体的目标1的残差比更新了目标2的个体的目标1的残差大，则不更新更新了目标2的个体的目标1；如果选择的个体的目标1的残差比更新了目标2的个体的目标1的残差小，则用替换来更新更新了目标2的个体的目标1，在这里为了方便，不管新子代种群中相应个体的目标2有没有更新，都将此操作后的个体记为于是得到每一块的更新子代种群其中是更新子代种群中第j个个体的目标1，s_j′是更新子代种群中第j个个体的目标2，j＝1,…,10；

11.6)计算每一块的更新子代种群中每个个体的适应度。

步骤112：如果当前进化代数p大于设置的最大进化代数mp，或者适应度达到要求，则执行步骤十三，并保留最终得到的每一块的更新子代种群中适应度最大的个体作为每一块的最优小波系数解；否则p=p+1，并将每一块的更新子代种群覆盖每一块的种群P，返回步骤八。

步骤113：对所有块的最优小波系数解执行打散块合并操作，即恢复成分块打散压缩采样前的形式，结合保留的低频小波系数Cl进行小波逆变换得到测试图像的重构图像。

本发明的具体实施从多目标进化算法和稀疏优化问题的特性出发，通过设计基于图像先验信息的系数种群和稀疏度集合的种群初始化方案和基于图像先验信息的一致变异操作和基于图像先验信息的迭代硬阈值(IHT)局部搜索策略，将基于分解的多目标进化算法(MOEA/D)和IHT以及图像先验信息有机结合作为非凸压缩感知优化重构方法，实现在小波域下从全局出发求解最优的小波系数解。

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

1.仿真条件：

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

2.仿真内容与结果：

(1)仿真1：

本仿真中，使用现有压缩感知技术中的分块策略和本发明所采用的分块打散策略对大小为512×512的StomachCT图像和Lena图像，在小波域下进行分块策略比较，其比较结果如图4和图5所述。其中：

图4(a)是StomachCT图像在分块策略下每块维数是256的情况下的大系数统计图，

图4(b)是StomachCT图像在分块打散策略下每块维数是256的情况下的大系数统计图，

图4(c)是StomachCT图像在分块策略下每块维数是1024的情况下的大系数统计图，

图4(d)是StomachCT图像在分块打散策略下每块维数是1024的情况下的大系数统计图；

图5(a)是Lena图像在分块策略下每块维数是256的情况下的大系数统计图，

图5(b)是Lena图像在分块打散策略下每块维数是256的情况下的大系数统计图，

图5(c)是Lena图像在分块策略下每块维数是1024的情况下的大系数统计图，

图5(d)是Lena图像在分块打散策略下每块维数是1024的情况下的大系数统计图。

从图4和图5中可以看出，现有压缩感知技术中的分块策略在小波域下会出现不稀疏或者弱稀疏的情况，违背了压缩感知的理论基础，而本发明所采用的分块打散策略并不会出现不稀疏或者弱稀疏的情况，所以我们采用本发明的分块打散策略能够得到更好的重构图像。

(2)仿真2：

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

图6(a)是StomachCT原始图像，

图6(b)是图6(a)的局部放大图，

图6(c)是OMP得到的StomachCT的重构图，

图6(d)是图6(c)的局部放大图，

图6(e)是IHT得到的StomachCT重构图，

图6(f)是图6(e)的局部放大图，

图6(g)是本发明得到的StomachCT重构图，

图6(h)是图6(g)的局部放大图。

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

用现有的OMP，IHT和本发明方法分别在采样率为30%、35%、40%、45%、50%的情况下，对大小为512×512的医学图像StomachCT的10次重构结果的峰值信噪比PSNR的平均值如表1所示。

表1各采样率下PSNR值

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

根据表1数据得到OMP，IHT和本发明方法重构出的StomachCT图像的PSNR随采样率变化的趋势图，其结果如图7所示，图7中的横坐标表示采样率，纵坐标表示峰值信噪比PSNR(dB)值。

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

(3)仿真3：

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

图8(a)是Lena原始图像，

图8(b)是图8(a)的局部放大图，

图8(c)是OMP得到的Lena的重构图，

图8(d)是图8(c)的局部放大图，

图8(e)是IHT得到的Lena重构图，

图8(f)是图8(e)的局部放大图，

图8(g)是本发明得到的Lena重构图，

图8(h)是图8(g)的局部放大图。

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

用现有的OMP，IHT和本发明方法分别在采样率为30%、35%、40%、45%、50%的情况下，对大小为512×512的自然图像Lena的10次重构结果的峰值信噪比PSNR的平均值如表2所示。

表2各采样率下PSNR值

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

根据表2数据得到OMP，IHT和本发明方法重构出的Lena图像的PSNR随采样率变化的趋势图，其结果如图9所示，图9中的横坐标表示采样率，纵坐标表示峰值信噪比PSNR(dB)值。

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

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

Claims (8)

1.基于进化多目标优化的非凸压缩感知图像重构方法，其特征是：至少包括如下步骤：

步骤101：输入测试图像，对测试图像进行小波变换，保留低频小波系数C_l，利用正交随机高斯观测矩阵Φ对高频小波系数每隔8个像素进行分块打散压缩采样，得到每一块的观测向量y；

利用正交随机高斯观测矩阵Φ对高频小波系数每隔8个像素进行分块打散压缩采样,包括如下步骤：

步骤201：将高频小波系数分成8×8的块；

步骤202：按从左到右，从上到下的顺序分别取每个8×8的块的一个元素也是按照从左到右，从上到下的顺序取组成打散块信号，这样可得到64个打散块信号；

步骤203：利用正交随机高斯观测矩阵Φ对64个打散块信号分别观测；

步骤102：利用保留的低频小波系数C_l，通过小波逆变换、边缘检测和分块打散操作得到对应于高频小波系数的位置矩阵E，位置矩阵E的每一列e对应于每一块的观测向量y；

所述的分块打散操作，包括如下步骤：

步骤401：将边缘图I₁的高频小波系数矩阵B分成8×8的块；

步骤402：按从左到右，从上到下的顺序分别取每个8×8的块的一个元素组成打散块信号，这样可得到由64个打散块信号组成的矩阵；

步骤103：利用位置矩阵E通过统计方法获取每一块的块稀疏度范围[n,n+120]，其中n是位置矩阵E中对应于观测向量y的那一列的非零元素个数；

步骤104：对位置矩阵E执行打散块合并操作，得到分块打散前的位置矩阵E′；

所述对位置矩阵E执行打散块合并操作，包括如下步骤：

步骤601：按从左到右的顺序获取位置矩阵E中的每一列；

步骤602：将位置矩阵E中每一列的元素按从左到右、从上到下的顺序每隔8个像素放置在分块打散前的位置上，得到分块打散前的位置矩阵E′；

步骤105：对分块打散前的位置矩阵E′进行提取疫苗和注射疫苗操作，得到分块打散前的过渡位置矩阵E″；

步骤106：对分块打散前的过渡位置矩阵E″执行分块打散操作，得到新的位置矩阵新的位置矩阵的每一列对应于每一块的观测向量y；

步骤107：设置种群大小为10，当前进化代数为p，种群中每个个体由编码对由目标1、目标2组成，目标1对应于要重构的块的小波系数，目标2表示对应的块稀疏度，利用正交随机高斯观测矩阵Φ和每一块观测向量y的广义逆，并结合新的位置矩阵中的每一列对种群中每个个体的目标1进行基于块位置信息的初始化操作和根据对应于每一块的观测向量y的块稀疏度范围[n,n+120]对种群中每个个体的目标2进行初始化，得到每一块的种群P＝{(w₁′,s₁),...,(w₁₀′,s₁₀)}，其中w_j′是种群P中的第j个个体的目标1，即对应的块的小波系数，s_j是种群P中的第j个个体的目标2，即对应的块稀疏度，j＝1,…,10；

步骤108：对每一块的种群P中每个个体的目标1根据新的位置矩阵及其目标2进行一致变异操作，得到每一块的子代种群其中是子代种群中第j个个体的目标1，s_j是子代种群中的第j个个体的目标2，j＝1,…,10；

步骤109：对每一块的子代种群中每个个体的目标1根据对应于每一块的新的位置矩阵中的一列及其目标2进行迭代硬阈值局部搜索，得到每一块的新子代种群其中是新子代种群中第j个个体的目标1，s_j是新子代种群中的第j个个体的目标2，j＝1,…,10，并计算每一块的新子代种群中每个个体的适应度；

步骤110：从由每一块的新子代种群中每个个体的目标1组成的集合中进行非支配解选择，并将非支配解对应地保存在每一块的非支配解集合Q中，但不保存弱支配解；

步骤111：根据每一块的非支配解集合Q中的元素个数更新新子代种群中相应个体的目标2，即稀疏度，若新子代种群中相应个体的目标2有更新，则更新每一块的新子代种群中更新了目标2的个体的目标1；若新子代种群中相应个体的目标2没有更新，则保留新子代种群中的个体，在这里为了方便，不管新子代种群中相应个体的目标2有没有更新，都将此操作后的个体记为于是得到每一块的更新子代种群其中是更新子代种群中第j个个体的目标1，s_j′是更新子代种群中第j个个体的目标2，j＝1,…,10，并计算每一块的更新子代种群中每个个体的适应度；

步骤112：如果当前进化代数p大于设置的最大进化代数mp，或者适应度达到要求，则执行步骤113，并保留最终得到的每一块的更新子代种群中适应度最大的个体作为每一块的最优小波系数解；否则p＝p+1，并将每一块的更新子代种群覆盖每一块的种群P，返回步骤108；

步骤113：对所有块的最优小波系数解执行打散块合并操作，即恢复成分块打散压缩采样前的形式，结合保留的低频小波系数C_l进行小波逆变换得到测试图像的重构图像。

2.根据权利要求1所述的基于进化多目标优化的非凸压缩感知图像重构方法，其特征是：所述步骤102中的通过小波逆变换、边缘检测和分块打散操作得到对应于高频小波系数的位置矩阵E，包括如下步骤：

步骤301：将高频小波系数设定为零，再利用保留的低频小波系数C_l，通过小波逆变换得到一个边缘模糊的图像I₀；

步骤302：利用Canny边缘检测方法对边缘模糊的图像I₀进行边缘检测，并提取边缘图像得到边缘图I₁；

步骤303：对边缘图I₁进行一层小波变换，得到边缘图I₁的高频小波系数矩阵B；

步骤304：对边缘图I₁的高频小波系数矩阵B每隔8个像素进行分块打散操作得到对应的分块打散操作后的边缘图I₁的高频小波系数矩阵其中的每一列对应于每一块的观测向量y；

步骤305：分别对分块打散操作后的边缘图I₁的高频小波系数矩阵的每一列求均值，并将均值的0.6倍作为阈值对的每一列进行阈值处理，将每一列中大于0.6倍均值的元素置为1；否则，置为0；得到阈值处理后的对应于分块打散操作后的边缘图I₁的高频小波系数矩阵的位置矩阵E，其中位置矩阵E的每一列e对应于每一块的观测向量y。

3.根据权利要求1所述的基于进化多目标优化的非凸压缩感知图像重构方法，其特征是：所述步骤103中的利用位置矩阵E通过统计方法获取每一块的块稀疏度范围[n,n＋120]，包括如下步骤：

步骤501：统计位置矩阵E中每一列的非零元素个数n；

步骤502：将每一列的非零元素个数n加上120，得到对应于每一块的观测向量y的块稀疏度范围[n,n＋120]。

4.根据权利要求1所述的基于进化多目标优化的非凸压缩感知图像重构方法，其特征是：所述步骤105中的对分块打散前的位置矩阵E′进行提取疫苗和注射疫苗操作，包括如下步骤：

步骤701：在分块打散采样前的位置矩阵E′上用3×3的窗口进行滑动，考虑3×3的窗口的上、下、左、右四个相邻的值，如果这四个值为1的个数大于等于3，则将1作为该窗口的中心位置的值；如果四个值为0的个数为4，则将0作为中心位置的值；若两种情况都不满足，则该中心位置上的值不变。

5.根据权利要求1所述的基于进化多目标优化的非凸压缩感知图像重构方法，其特征是：所述步骤107中的结合新的位置矩阵中的每一列对种群中每个个体的目标1进行基于块位置信息的初始化操作和根据对应于每一块的观测向量y的块稀疏度范围[n,n＋120]对种群中每个个体的目标2进行初始化，包括如下步骤：

步骤801：设置种群大小为10，种群中每个个体由编码对(目标1，目标2)组成，目标1对应于要重构的块的小波系数，目标2表示对应的块稀疏度，将种群记为P＝{(w₁′,s₁),...,(w₁₀′,s₁₀)}，其中w_j′是种群P中的第j个个体的目标1，即对应的块的小波系数，且每一个w_j′都采用实数编码，s_j是种群P中的第j个个体的目标2，即对应的块稀疏度，且每一个s_j都采用整数编码，j＝1,…,10；

步骤802：根据正交随机高斯观测矩阵Φ和每一块的观测向量y通过广义逆求得初始的小波系数，公式如下：

步骤803：我们结合广义逆值和新的位置矩阵中对应于观测向量y的一列来初始化种群P中每个个体的目标1，则种群P中每个个体的目标1w_j′(j＝1,…,10)可按下面的公式进行初始化：

其中，L是伸缩因子，且1≤L≤1.5，t＝1,…,1024；

步骤804：对种群P中每个个体的目标2进行初始化，每个个体的目标1w_j′对应的目标2s_j都是从对应于每一块的观测向量y的块稀疏度范围[n,n+120]中随机选取，其中n是位置矩阵E中对应于观测向量y的那一列的非零元素个数，且每个个体的目标2s_j都不同。

6.根据权利要求1所述的基于进化多目标优化的非凸压缩感知图像重构方法，其特征是：所述步骤108中的一致变异操作，包括如下步骤：

步骤901：统计种群P中第j个个体的目标1w_j′的非零元素个数作为第j个个体的目标1w_j′的真实块稀疏度k_j，j＝1,…,10；

步骤902：如果第j个个体的目标1w_j′的真实块稀疏度k_j(j＝1,…,10)小于第j个个体的目标2s_j(j＝1,…,10)，则在第j个个体的目标1w_j′的零元素中随机选取s_j-k_j个位置，并用一次迭代硬阈值(IHT)求得的值来替换这些位置上的零元素；如果第j个个体的目标1w_j′的真实块稀疏度k_j(j＝1,…,10)大于第j个个体的目标2s_j，则在第j个个体的目标1w_j′的非零元素中选取k_j-s_j个值最小的位置，并用0来替换这些位置上的值。

7.根据权利要求1所述的基于进化多目标优化的非凸压缩感知图像重构方法，其特征是：所述步骤109中的迭代硬阈值局部搜索，包括如下步骤：

步骤1001：对于子代种群中第j个个体的目标1将其赋值给即

步骤1002：根据迭代硬阈值的公式计算20次，并将赋值给即其中，是用位置和第j个个体的目标2s_j指导的非线性算子；

所述步骤109中计算每一块的新子代种群中每个个体的适应度，是按照如下公式计算每一个个体的适应度

其中，φ是正交随机高斯观测矩阵，y是每一个个体对应的块的观测向量，表示向量2范数的平方。

8.根据权利要求1所述的基于进化多目标优化的非凸压缩感知图像重构方法，其特征是：所述步骤111中的根据每一块的非支配解集合Q中的元素个数更新新子代种群中相应个体的目标2，包括如下步骤：

步骤1101：将所有个体的目标2s_j(j＝1,…,10)按升序排序，即其中j₁,…,j10∈{1,…,10}；

步骤1102：根据下面的公式确定待选块稀疏度集

步骤1103：如果从中选择一个值来更新块稀疏度，

如果|Q|＞0.5×10，其中r是从中随机选择的一个值；

如果|Q|≤0.5×10，其中r是从中随机选择的一个值；

步骤1104：如果且

如果|Q|＝10，其中n+120是块稀疏度范围[n,n+120]的最大值；

如果|Q|＝1，其中n是块稀疏度范围[n,n+120]的最小值；

所述步骤111中的更新每一块的新子代种群中更新了目标2的个体的目标1，按如下过程进行：

若新子代种群中相应个体的目标2有更新，则选取每一块的新子代种群中与更新了目标2的个体的目标2最相近的个体如果选择的个体的目标1的残差比更新了目标2的个体的目标1的残差大，则不更新更新了目标2的个体的目标1；如果选择的个体的目标1的残差比更新了目标2的个体的目标1的残差小，则用替换来更新更新了目标2的个体的目标1。