技术领域
[0001] 本
发明属于
图像处理技术领域，涉及压缩感知图像重构方法，可用于
生物传感、遥感图像处理、无线传感网络和
图像采集设备的开发等实际工程领域。
背景技术
[0002] 随着数字化和
信息时代的来临，数字图像处理技术越来越受到人们的关注，特别是在SAR图像处理、医学图像处理以及遥感图像处理等领域。近年来，Donoho等人提出了一种新颖的理论-压缩感知理论CS。在该理论中，
采样速度不再决定于
信号的带宽，而决定于信号中的结构和内容。基于压缩感知理论，信号可以进行低速采样，然后进行编码，大大降低了计算复杂度。该理论一经提出，就受到了广泛关注，不断被应用到数字图像处理中。
[0003] CS理论是一种新的在采样的同时实现压缩目的的理论
框架：假设信号x∈RN在某T个
正交基或紧框架Ψ上是可压缩的，首先求出变换系数Θ＝Ψx，Θ是x的等价或逼近的稀疏表示；然后，设计一个平稳的、与变换基Ψ不相关的M×N维的
观测矩阵Φ，对Θ进T
行观测，即将Θ投影到M维空间，得到观测集合y＝ΦΘ＝ΦΨx，该过程也可以表示为CS CS CS T CS
信号x通过矩阵A 进行非自适应观测：y＝A x，其中A ＝ΦΨ，A 称为CS信息算子；
最后，利用l0范数意义下的优化问题求解x的精确或近似逼近
[0004] min||ΨTx||0 s.t.ACSx＝ΦΨTx＝y 1)
[0005] 求得的向量 就是在基Ψ上的最稀疏表示。
[0006] CS理论表明：利用1)式可以用M＝S+1个独立同分布的高斯测度以大概率精确重构S-
稀疏信号。但对1)式求解时，时间复杂度很高，需要穷举Θ中所有可能的 个非零项的组合。
[0007] 研究表明，1)式的转化问题为：
[0008] min||ΨTx||1 s.t.ACSx＝ΦΨTx＝y 2)
[0009] 即基于l1范数的优化问题。对该问题的求解仍能够精确重构S-稀疏信号，并且仅用M≥CSlog(N/S)个独立同分布的高斯测度就能以大概率逼近可
压缩信号。这是一个3
凸优化问题，可以很方便地简化为基追踪的线性规划问题，它的计算复杂度约为O(N)。
[0010] 针对以上问题，近年来，数学家和工程师们提出了很多压缩感知重构
算法。一类是凸松弛法，通过将非凸问题转化为凸问题来找到原始信号的逼近，如基追踪BP和贪婪基追踪GBP；另一类是贪婪追踪算法，通过每次
迭代时选择一个局部最优解来逼近原始信号，如匹配追踪MP、正交匹配追踪OMP、分段匹配追踪StOMP和正则匹配追踪ROMP。这两类算法都有其固有的缺点，凸松弛法重构信号时所用的观测数较小，重构出来的图像
质量较高，但其运算量大，时间成本高；贪婪追踪算法和凸松弛法相比，克服了计算时间长的问题，适合用来求解大规模的问题，但其需要的观测数较多，重构出来的图像质量不高，并且它对压缩感知框架强加了有限等距性RIP约束，从某种意义上讲，限制了CS的应用范围。
发明内容
[0011] 本发明的目的在于针对上述已有技术的问题，提出一种基于先验模型的压缩感知图像重构方法，通过l0范数，扩展压缩感知CS的应用范围，降低图像重构的计算复杂度，提高图像的重构质量。
[0012] 为实现上述目的，我们首先关注的是如何确定稀疏系数的
位置，然后在该位置信息的指导下去求解稀疏系数值。由此本发明建立了基于先验模型的免疫优化压缩感知重构框架。重构框架的第一部分是利用先验模型用免疫
遗传算法来较快地搜索稀疏系数在小波高频子带中的位置；重构框架的第二部分利用改进的克隆选择算法去求解l0范数意义下该位置上的稀疏系数值。本发明的技术方案包括如下步骤：
[0013] (1)对发送方发送过来的一尺度小波低频子带做逆变换，高频系数全部置零，得到边缘模糊的图像；对该模糊图像进行canny
边缘检测，获取含有边缘信息的图像；对含有边缘信息的图像做
小波变换得到三个含有边缘信息的高频子带，按照系数模值的统计分布确定出每个高频子带中系数的模值为大、小和零的位置，分别标记为1、0.5和0；对于三个含有边缘信息的高频子带，当边缘检测的
阈值参数分别取n个不同的值时，就会产生三个包含n个位置矩阵的位置种群A(k)、B(k)和C(k)：
[0014] A(k)＝{a1(k)，a2(k)，Λ，an(k)}
[0015] B(k)＝{b1(k)，b2(k)，Λ，bn(k)}
[0016] C(k)＝{c1(k)，c2(k)，Λ，cn(k)}，
[0017] 其中，对于256×256的图像，ai(k)、bi(k)和ci(k)分别是128×128的位置编码矩阵，i＝1，2，Λ，n，n为种群中个体数目；
[0018] (2)设定种群进化的终止条件为迭代次数t达到最大迭代次数p，并初始化迭代次数t＝0；根据图像的发送方发送的和观测向量y和观测矩阵Φ，计算初始化的一个高频子+ +带系数矩阵：d0(k)＝Φy，Φ 是观测矩阵Φ的广义逆；
[0019] (3)利用位置种群A(k)，用系数矩阵d0(k)初始化系数种群D(k)，即[0020] D(k)＝{d1(k)，d2(k)，Λ，dn(k)}
[0021]
[0022] 其中，对于256×256的图像，di(k)是一个128×128的高频子带系数矩阵，i＝1，2，Λ，n，L1和L2为1到2之间的随机数，L1＜L2，k＝1，2，...，128；
[0023] (4)对位置种群A(k)中的每个个体ai(k)执行提取
疫苗操作和注射疫苗操作，得到位置种群A′(k)：A′(k)＝{a′1(k)，a′2(k)，Λ，a′n(k)}；
[0024] (5)将位置种群A′(k)和系数种群D(k)作为输入，运行克隆选择，输出新的位置种群A″(k)和新的系数种群D′(k)：
[0025] A″(k)＝{a″1(k)，a″2(k)，Λ，a″n(k)}
[0026] D′(k)＝{d′1(k)，d′2(k)，Λ，d′n(k)}；
[0027] 其中，a″i(k)是第i个位置编码矩阵，d′i(k)是第i个高频子带系数矩阵；
[0028] (6)利用如下公式计算系数种群D′(k)中每个个体d′i(k)的适应度值f(d′i(k))：
[0029]
[0030] 其中，τ取0.00012；
[0031] (7)选择适应度值最大的两个个体d′g(k)和d′e(k)，执行交叉操作，生成两个新的个体d″g(k)和d″e(k)，并计算它们的适应度值，同时在位置种群A″(k)中对d′g(k)和d′e(k)相应的位置矩阵的同一位置进行交换，其中g＝1，2，Λ，n，e＝1，2，Λ，n，且g≠e；
[0032] (8)从d″g(k)、d″e(k)、d′g(k)和d′e(k)中选择适应度值最大的个体，将这个个体与其对应的位置矩阵一同存到记录种群G(k)中；
[0033] (9)如果迭代次数t满足步骤2中的终止条件，则在记录种群G(k)中选出最优的个体，即所求的第一个小波高频子带，转入下一步；否则，返回步骤(1)和步骤(2)，产生n-2个新的系数矩阵和位置矩阵，分别替换原来的系数种群D′(k)中除d′g(k)和d′e(k)之外的n-2个个体和位置种群A″(k)中对应的n-2个个体，转入步骤(3)，t＝t+1；
[0034] (10)分别利用步骤(1)中的另外两个位置种群B(k)和C(k)替换步骤(2)到步骤(9)中的A(k)，重复执行两遍，求出其它两个小波高频子带；
[0035] (11)对得到的三个小波高频子带，以及发送方发送的一个低频子带，做一尺度的小波逆变换得到重构的图像。
[0037] 1.计算复杂度低，运行时间短。
[0038] 压缩感知图像重构的问题本质上是求解l0范数的问题，然而目前的一些算法采用l1范数方法求解，如凸松弛法，通过将非凸问题转化为凸问题来找到原始信号的逼近，虽然重构信号时所用的观测数较小，重构出来的图像质量较高，但其运算量大，时间成本高。本发明的方法在计算复杂度较低的情况下，利用了现有的方法没有用到的稀疏系数的位置分布，确定了稀疏系数在小波高频子带中的准确位置；然后再求解该位置上的稀疏系数值，从而提高了计算速度，缩短了运行时间。
[0039] 2.扩展压缩感知CS的应用范围，提高了图像的重构质量。
[0040] 现有的一些压缩感知图像重构方法，如OMP和ROMP，虽然运行时间比较快，但它本身是用于寻找最优的稀疏字典的方法，若用到重构问题上，则没有考虑图像在小波变换域中的系数具有聚集性的这一特征，也难以得到较好的重构效果。本发明结合现有重构算法中的问题，从l0范数重构问题出发，提出了基于先验模型的免疫优化压缩感知重构框架，考虑了图像在小波域中系数的聚集性，在可接受的时间内能较准确地确定稀疏系数在某尺度内的位置，稳定并精确地求解该位置上的稀疏系数值，从而提高了图像的重构质量。
附图说明
[0042] 图2是本发明中设计先验模型的子流程图；
[0043] 图3是用本发明提取疫苗和注射疫苗的示意图；
[0044] 图4是用本发明进行分
块操作的示意图；
[0045] 图5是用本发明进行互换操作的示意图；
[0046] 图6是用本发明进行交叉操作的示意图；
[0047] 图7是本发明与BP，BCS和MP重构出来的Lena图的仿真对比图；
[0048] 图8是本发明与BP，BCS和MP重构出来的Lena图的三个高频子带的仿真对比图，显示的是系数的模值取反的图像；
[0049] 图9是本发明重构出来的Lena图的三个高频子带的系数分布和位置分布的仿真对比图；
[0050] 图10是本发明与BP，BCS和MP重构出来的Lena图的峰值
信噪比PSNR和重构时间随采样率变化的趋势图；
[0051] 图11是种群大小为10时，本发明重构出来的Lena图的
峰值信噪比PSNR随迭代次数变化的趋势图；
[0052] 图12是最大迭代次数为2时，本发明重构出来的Lena图的峰值信噪比PSNR随种群大小变化的趋势图；
[0053] 图13是本发明重构出来的Lena图的峰值信噪比PSNR随参数L1和L2变化的趋势图。
具体实施方式
[0054] 参照图1，本发明的具体实施步骤如下：
[0055] 步骤1，利用小波低频子带求出三个高频子带对应的三个位置种群。
[0056] 参照图2，本步骤的具体实现如下：
[0057] (1a)对发送方发送过来的一尺度小波低频子带做逆变换，高频系数全部置零，得到边缘模糊的图像；
[0058] (1b)对该模糊图像进行canny边缘检测，获取含有边缘信息的图像；
[0059] (1c)对含有边缘信息的图像做小波变换得到三个含有边缘信息的高频子带，按照系数模值的统计分布确定出每个高频子带中系数的模值为大、小和近似为零的位置，分别标记为1、0.5和0；
[0060] (1d)对于三个含有边缘信息的高频子带，当边缘检测的阈值参数分别取n个不同的值时，就会产生三个包含n个位置矩阵的位置种群A(k)、B(k)和C(k)：
[0061] A(k)＝{a1(k)，a2(k)，Λ，an(k)}
[0062] B(k)＝{b1(k)，b2(k)，Λ，bn(k)}
[0063] C(k)＝{c1(k)，c2(k)，Λ，cn(k)}，
[0064] 其中，对于256×256的图像，ai(k)、bi(k)和ci(k)分别是128×128的位置编码矩阵，i＝1，2，Λ，n，n为种群中个体数目。
[0065] 步骤2，设定种群进化的终止条件为迭代次数t达到最大迭代次数p，并初始化迭代次数t＝0；根据图像的发送方发送的和观测向量y和观测矩阵Φ，计算初始化的一个高+ +频子带系数矩阵：d0(k)＝Φy，Φ 是观测矩阵Φ的广义逆。
[0066] 步骤3，利用位置种群A(k)，用系数矩阵d0(k)初始化系数种群D(k)，即[0067] D(k)＝{d1(k)，d2(k)，Λ，dn(k)}
[0068]
[0069] 其中，对于256×256的图像，di(k)是一个128×128的高频子带系数矩阵，i＝1，2，Λ，n，L1和L2为1到2之间的数，L1＜L2，k＝1，2，...，128；
[0070] 步骤4，对位置种群A(k)中的每个个体ai(k)执行提取疫苗操作和注射疫苗操作，得到位置种群A′(k)。
[0071] 参照图3，本步骤的具体实现如下：
[0072] (4a)用3×3的窗口在位置种群A(k)中的每个个体ai(k)中滑动，根据3×3窗口中中心位置的上、下、左、右四个相邻的值提取疫苗，如果四个值为0.5的个数大于3，则将0.5作为该中心位置的疫苗；如果四个值为1的个数大于3，则将1作为该中心位的疫苗，如图3(a)所示操作；如果四个值为0的个数为4，则将0作为该中心位置的疫苗，如图3(b)所示操作；
[0073] (4b)再次用3×3的窗口在位置种群A(k)中的每个个体ai(k)中滑动，若3×3窗口中心位置上的值和步骤(4a)提取出来的疫苗不相同，则用疫苗的值来替换中心位置上的值，如图3(a)和3(b)所示操作；反之，不进行注射操作；
[0074] 经过上面两步，对位置种群A(k)中的每个个体进行操作后，得到位置种群A′(k)：
[0075] A′(k)＝{a′1(k)，a′2(k)，Λ，a′n(k)}。
[0076] 步骤5，将位置种群A′(k)和系数种群D(k)作为输入，运行克隆选择，输出新的位置种群A″(k)和新的系数种群D′(k)。
[0077] (5a)参照图4，把位置种群A′(k)和系数种群D(k)中的所有个体分成m个16×16大小的图像块，用a′ij(k)表示位置种群A′(k)中第i个个体的第j个块，用dij(k)表示系数种群D(k)中第i个个体的第j个块，其中，i＝1，2，Λ，n，j＝1，2，Λ，m；n是种群中个体数目，m是每个个体分成的块的个数；
[0078] (5b)对系数种群D(k)中第i个个体的第j个块dij(k)进行克隆，克隆次数为2，得到克隆块 和
[0079] (5c)参照图5，对位置种群A′(k)中的块a′ij(k)逐列进行扫描，若有连续两个以上的位置同时出现1，则把克隆块 和 对应位置上的系数交换位置；若有连续两个以上的位置同时出现0.5，则只把克隆块 中对应位置上的系数交换位置；
[0080] (5d)分别计算克隆前的块dij(k)以及克隆后的块 和 的适应度值：
[0081]
[0082]
[0083]
[0084] 其中，y和Φ分别是图像的发送方发送的观测向量和观测矩阵，λ＝0.08；
[0085] (5e)选出上述适应度值f(dij(k))、 和 中最大值对应的块，把该块作为第i个个体的第j个块；合并所有个体的m个块，生成新的系数种群D′(k)：
[0086] D′(k)＝{d′1(k)，d′2(k)，Λ，d′n(k)}，
[0087] 同时合并位置种群中所有个体的m个块，生成新的位置种群A″(k)：
[0088] A″(k)＝{a″1(k)，a″2(k)，Λ，a″n(k)}，
[0089] 其中，a″i(k)是第i个位置编码矩阵，d′i(k)是第i个高频子带系数矩阵，i＝1，2，Λ，n，n是种群中个体数目。
[0090] 步骤6，利用如下公式计算系数种群D′(k)中每个个体d′i(k)的适应度值f(d′i(k))：
[0091]
[0092] 其中，y和Φ分别是图像的发送方发送的观测向量和观测矩阵，τ取0.00012。
[0093] 步骤7，选择步骤6中所有适应度值中最大的两个个体d′g(k)和d′e(k)，执行交叉操作，生成两个新的个体d″g(k)和d″e(k)，其中g＝1，2，Λ，n，e＝1，2，Λ，n，且g≠e。
[0094] 参照图6，本步骤的具体实现如下：
[0095] (7a)将系数种群D′(k)中的两个个体d′g(k)和d′e(k)分成16×16大小的块；
[0096] (7b)从两个个体的所有16×16的块中分别随机选择两个块，以概率0.2交换这两个块，生成新的个体d″g(k)和d″e(k)，图6显示的是两个个体交换第28个块的操作过程。
[0097] 步骤8，在位置种群A″(k)中对d′g(k)和d′e(k)相应的位置矩阵的同一位置进行交换，计算两个新个体d″g(k)和d″e(k)的适应度值，并从d″g(k)、d″e(k)、d′g(k)和d′e(k)中选择适应度值最大的个体，将这个个体与其对应的位置矩阵一同存到记录种群G(k)中。
[0098] 步骤9，如果迭代次数t满足步骤2中的终止条件，则在记录种群G(k)中选出最优的个体，即所求的第一个小波高频子带，转入下一步；否则，返回步骤(1)和步骤(2)，产生n-2个新的系数矩阵和位置矩阵，分别替换原来的系数种群D′(k)中除d′g(k)和d′e(k)之外的n-2个个体和位置种群A″(k)中对应的n-2个个体，转入步骤(3)，t＝t+1。
[0099] 步骤10，分别利用步骤(1)中的另外两个位置种群B(k)和C(k)替换步骤(2)到步骤(9)中的A(k)，重复执行两遍，求出其它两个小波高频子带。
[0100] 步骤11，对得到的三个小波高频子带，以及发送方发送的一个低频子带，做一尺度的小波逆变换得到重构的图像。
[0101] 本发明的优点由以下仿真的数据和图像进一步说明。
[0102] 1.仿真条件
[0103] (1.1)选取256×256的标准测试图像Lena来测试本发明的重构效果，同时将本发明的重构效果与现有的BP，BCS和MP的重构效果进行比较来突出本发明的重构优势。
[0104] (1.2)仿真实验中的种群大小n为10，最大迭代次数p为2。
[0105] (1.3)仿真实验中初始化系数种群时用到的伸缩因子L1为1，L2为1.8。
[0106] (1.4)仿真实验中图像分成的块的大小定为16×16。
[0107] 2.仿真内容
[0108] (2.1)根据上述的具体实施方式，利用本发明的方法来重构Lena图，并用现有的方法BP，BCS和MP作为对比实验，仿真结果如图7所示。其中图7(a)是原图，图7(b)是本方法的重构图，图7(c)、7(d)和7(e)分别是现有BP，BCS和MP方法的重构图。图10(a)给出了本发明与现有BP，BCS和MP方法重构图像的时间随采样率变化的趋势图；图10(b)是本发明与现有BP，BCS和MP方法重构图像的峰值信噪比PSNR随采样率变化的趋势图。
[0109] (2.2)为了说明先验模型设计的重要性和有效性，利用本发明的方法来重构Lena图的
水平方向高频子带，并用现有的BP，BCS和MP重构出来的图像的水平方向高频子带作为对比试验，结果显示的是系数的模值取反的图像，如图8所示。其中图8(a)显示了Lena原图在水平方向高频子带上的小波系数值，图8(b)显示了本发明重构出来的水平方向高频子带上的小波系数值，图8(c)、8(d)和8(e)分别显示了现有的BP，BCS和MP方法重构出来的水平方向高频子带上的小波系数值。为了说明本发明重构出来的小波系数在结构上具有聚集性，用本发明的方法重构出三个高频子带，并显示其对应的位置编码矩阵，如图9所示。其中9(a)、9(c)和9(e)分别是本发明重构出来的三个高频子带中稀疏系数的分布情况；图9(b)、9(d)和9(f)分别是其对应的三个位置编码矩阵，图中白色区域代表“大”系数的位置，灰色区域代表“小”系数的位置，黑色区域代表系数值为“零”的位置。
[0110] (2.3)本发明在实现的过程中用到了很多参数，其中一部分参数是与图像无关的，它们对于所有的输入图像使用的是相同的值，如种群大小n和最大迭代次数p。种群大小n实际上是由边缘检测高阈值的个数决定的，代表着位置种群中位置矩阵的个数，是本发明的一个关键参数。通过实验对种群大小n和最大迭代次数p这两个参数进行分析，其试验是在用这两个参数之一来观察另一个参数对实验结果的影响：
[0111] (2.3.a)当种群大小n固定为10时，观察本发明重构出来的图像的峰值信噪比PSNR随迭代次数p增加的变化趋势，结果如图11所示。
[0112] (2.3.b)当最大迭代次数p固定为2时，观察本发明重构出来的图像的峰值信噪比PSNR随种群大小n增大的变化趋势，结果如图12所示。
[0113] (2.4)另外还有一些参数是与图像相关的，它们是在大量试验下确定的合适的参数值，如伸缩因子L1和L2。为了分析这两个参数，设计实验如下：伸缩因子L2取2.5、1.9、1.8和1.6时，观察峰值信噪比PSNR随L1取值不同的变化趋势，这样可以得到4个不同的曲线，如图13所示。
[0114] 3.仿真结果分析
[0115] (3.1)从图7(c)可以看出，本发明重构出来的图像与原图7(a)非常接近，而现有的BP，BCS和MP方法重构出来的图像在视觉上要差很多。从图10(a)中可看出本方法比现有的BP，BCS和MP方法获得了较高的峰值信噪比PSNR，即重构的图像质量要比其他方法高；从图10(b)看出本发明与现有的BP，BCS和MP方法的运行时间相比，虽然时间不是最短的，但比BCS方法的时间要短很多，随着采样率的增加，本发明与BP和MP方法的运行时间相差不大。
[0116] (3.2)从图8可以看出，本发明重构出来的小波系数值和原图的小波系数值非常吻合，并且能够较准确地
定位原图在三个高频子带中的稀疏系数的位置，而现有的BP，BCS和MP重构出来的结果和原图有很大的差距，甚至连轮廓都没有重构出来。在图9(a)、9(c)和9(e)中，本发明重构出来的三个高频子带上的稀疏系数分布基本上也能够和图9(b)、9(d)和9(f)中位置矩阵的“大”系数和“小”系数的位置相吻合。
[0117] (3.3)由图11和图12可以看出，随着种群大小n和最大迭代次数p的增大，图像的峰值信噪比PSNR也在增加，但是运行时间会变得很长。为了在运行时间和重构效果之间取折中，当种群大小n取10，最大迭代次数p取2时，本发明可以在短时间内取得较好的重构效果。
[0118] (3.4)由图13中的四条曲线可以看出，要获得较大的峰值信噪比PSNR，必须满足1≤L1＜L2≤2，并且当L1＝1，L2＝1.8时，峰值信噪比PSNR相对来说达到最高。
[0119] 综上所述，本发明的基于先验模型和l0范数的压缩感知图像重构方法，能够很好地重构图像在小波域中的高频子带，与现有的其他方法相比，不仅具有较高的重构质量，视觉效果好，而且大大降低了计算量。
