基于双树复小波变换的多光谱图像重构方法

摘要

本发明公开了一种基于双树复小波变换的多光谱图像重构方法，解决了现有多光谱图像重构技术中重构效果不理想、重构速度慢的问题。本发明的步骤包括：(1)获取混叠光谱图像；(2)数据初始化；(3)降噪处理；(4)判断当前估计值的继续条件是否满足；(5)获取图像重构当前估计值的下一个估计值；(6)判断当前估计值的下一个估计值的继续条件是否满足；(7)更新估计值；(8)判断终止条件是否满足。本发明利用双树复小波变换对图像进行降噪处理，在压缩光谱成像的重构过程中，可获得较好的多光谱图像重构结果和较快的多光谱图像重构速度。

说明书

技术领域

本发明属于图像处理技术领域，特别涉及一种压缩光谱成像技术领域中基于 双树复小波变换的多光谱图像重构方法。本发明可用于多光谱图像的重构，提高 图像的重构质量，提升重构速度。

背景技术

在压缩光谱成像技术领域中，采集到的多光谱图像数据远少于原始的多光谱 图像数据，多光谱图像的重构过程就被转化为一个基于压缩感知理论的逆问题进 行求解。

Jos′e M.Bioucas-Dias，M′ario A.T.Figueiredo，在文献“A New TwIST: Two-Step Iterative Shrinkage/Thresholding Algorithms for Image Restoration”(IEEE  Transactions on Image processing,2007,16(12):2992-3004)中提出通过两个估计值 更新当前值并经过一系列的迭代来实现图像的逼近。该方法首先设置一个迭代初 始值作为第一个估计值，并对初始值进行全变差去噪得到第二个估计值，然后在 全变差去噪对前一估计值进行降噪处理的基础上，利用前两个估计值迭代得到新 的估计值，通过多次迭代最终实现原始图像的逼近。在进行多光谱图像重构时该 方法仍存在一些不足：首先全变差去噪的计算复杂度高、没有快速算法，且多光 谱图像的数据量较大，使得这种方法的重构速度较慢。其次全变差去噪对于图像 平滑区域抑噪不充分，会引入阶梯效应，导致这种方法的重构效果不理想。

西安电子科技大学申请的专利“基于压缩感知的双通道遥感光谱成像系统 和成像方法”(公开号：CN102706449A，申请号：CN201210172731，申请日： 2012年05月24日)中公开了一种利用稀疏变换和非线性优化的光谱图像重构 方法。该方法首先选取一个稀疏变换域，然后获取原始光谱图像在稀疏变换域下 的投影系数，最终利用非线性优化方法重构光谱图像。该专利存在的不足是，所 选用的稀疏变换域为小波变换域，它在平移不变性、方向性等方面仍有一些缺陷， 不具有完全重构性，使得这种方法的重构效果不是很理想。

综上所述，虽然现有方法都可以重构出多光谱图像，但是没有考虑到所用的 去噪方式所存在的缺陷，使得重构速度较慢且难以得到较好的重构图像。

发明内容

本发明针对上述已有技术的不足，提出一种基于双树复小波变换的多光谱图 像重构方法，以加快图像重构速度、提升重构效果。采用这种方法的原因是，双 树复小波变换不但具有小波变换的低熵性、去相关性和选基的灵活性，还具有平 移不变性、较好的方向性和精确的相空间信息；同时对于图像所采用的二维双树 复小波变换，无论其分解层数为多少，数据冗余总是4:1，使得计算复杂度大大 降低。

为实现本发明的目的，其原理是利用Kingsbury于1988年提出的双树复小 波变换对图像进行降噪处理。Kingsbury提出在同一个数据上，用两个平行独立 的小波变换来完成复数小波变换，也就是说，利用一对共轭正交滤波器组对复数 小波的实部进行分解，另一对共轭正交滤波器组对复数小波的虚部进行分解，就 完成双树复小波变换。

实现本发明的目的的具体步骤如下：

(1)获取混叠多光谱图像：

利用压缩光谱成像系统观测静止的图像场景，得到原始的混叠多光谱图像；

(2)数据初始化：

设定图像重构的初始估计值为0，图像重构的初始迭代步长为1，图像重构 的迭代终止系数为10^-6；

(3)降噪处理：

对图像重构的初始估计值进行降噪处理，将得到的图像重构的初始估计值的 降噪处理结果，作为图像重构的当前估计值；

所述的降噪处理的具体步骤如下：

第一步，按照下式，对图像重构的估计值进行优化：

$x_{\mathrm{opt}}=x+\frac{A^{-1}(y-\mathrm{Ax})}{s}$

其中，x表示图像重构的估计值，x_opt表示优化后的图像重构的估计值，y表 示原始的混叠多光谱图像，A表示压缩光谱成像系统所采用的观测算子，A^-1表 示对观测算子进行逆操作，s表示图像重构的迭代步长；

第二步，对优化后的图像重构的估计值进行双树复小波变换，得到变换系数；

第三步，对变换系数进行阈值寻优，得到最佳的阈值；将低于该阈值的变换 系数置为0，幅值高于该阈值的变换系数做相应的“收缩”处理，得到降噪后的 变换系数；

第四步，对降噪后的变换系数进行双树复小波逆变换，得到图像重构估计值 的降噪处理结果；

(4)判断图像重构当前估计值的继续条件是否满足，若是，则执行步骤(5)， 否则，图像重构的迭代步长加倍，执行步骤(3)；

(5)获取图像重构当前估计值的下一个估计值：

(5a)采用步骤(3)所述的降噪处理方法，对图像重构的当前估计值进行 降噪处理，得到图像重构当前估计值的降噪处理结果：

(5b)按照下式，获取图像重构当前估计值的下一个估计值：

x_t＝(1-α)·x_t-2+(α-β)·x_t-1+β·Γ(x_t-1)

其中，x_t表示图像重构当前估计值的下一个估计值，x_t-2表示图像重构当前 估计值的上一个估计值、x_t-1表示图像重构的当前估计值，Γ(x_t-1)表示图像重构 当前估计值的降噪处理结果，α表示当前估计值的上一个估计值的迭代系数，β 表示当前估计值降噪处理结果的迭代系数；

(6)判断图像重构当前估计值的下一个估计值的继续条件是否满足，若是， 则执行步骤(7)，否则，用图像重构的当前估计值更新图像重构的初始估计值， 执行步骤(3)；

(7)更新估计值：

(7a)用图像重构的当前估计值更新图像重构当前估计值的上一个估计值；

(7b)用图像重构当前估计值的下一个估计值更新图像重构的当前估计值；

(8)判断图像重构的终止条件是否满足，若是，则结束，此时图像重构的 结果为图像重构的当前估计值，否则，执行步骤(5)。

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

第一，本发明利用双树复小波变换对图像进行降噪处理，由于双树复小波变 换具有低熵性、去相关性、选基的灵活性、平移不变性、较好的方向性和精确的 相空间信息的特点，克服了现有小波变换域在平移不变性、方向性上的不足，使 得本发明提高了重构图像的质量，可获得峰值信噪比较高的重构结果。

第二，本发明利用双树复小波变换对图像进行降噪处理，充分发挥双树复小 波数据冗余有限和计算效率高的特点，克服了全变差去噪技术计算复杂度高、没 有快速算法的不足，使得本发明极大的减少了重构过程的迭代次数，提高了多光 谱图像的重构速度。

附图说明

图1为本发明的流程图；

图2为采用本发明与全变差去噪方法分别对balloons图重构结果的示意图；

图3为采用本发明与全变差去噪方法分别对fake_and_real_food图重构结果 的示意图。

具体实施方式

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

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

步骤1，获取混叠多光谱图像。

利用压缩光谱成像系统观测静止的图像场景，得到原始的混叠多光谱图像。

压缩光谱成像系统的基本成像过程如下：

首先多光谱图像的光束通过编码模板，编码模板对光束进行随机编码，然后 使编码后的光束经过色散棱镜，色散后的光束照射到阵列传感器上，便得到原始 的混叠多光谱图像。

步骤2，数据初始化。

设定图像重构的初始估计值为0，图像重构的初始迭代步长为1，图像重构 的迭代终止系数为10^-6。本发明实施例中，设定图像重构的初始迭代次数为1， 图像重构的迭代终止次数为1000，图像重构的迭代终止步长为2³⁰。

步骤3，降噪处理。

对图像重构的初始估计值进行降噪处理，将得到的图像重构的初始估计值的 降噪处理结果，作为图像重构的当前估计值。本发明实施例中，若执行完步骤3 时，图像重构的迭代步长大于图像重构的迭代终止步长，则结束，此时图像重构 的结果为图像重构的当前估计值，否则，继续执行步骤4。

所述的降噪处理的具体步骤如下：

第一步，按照下式，对图像重构的初始估计值进行优化：

$x_{\mathrm{opt}}=x+\frac{A^{-1}(y-\mathrm{Ax})}{s}$

其中，x表示图像重构的初始估计值，x_opt表示优化后的图像重构的估计值， y表示原始的混叠多光谱图像，A表示压缩光谱成像系统所采用的观测算子，A^-1表示对观测算子进行逆操作，s表示图像重构的迭代步长。

第二步，对优化后的图像重构的估计值进行双树复小波变换，得到变换系数。

第三步，对变换系数进行阈值寻优，得到最佳的阈值；将低于该阈值的变换 系数置为0，幅值高于该阈值的变换系数做相应的“收缩”处理，得到降噪后的 变换系数。

第四步，对降噪后的变换系数进行双树复小波逆变换，得到图像重构估计值 的降噪处理结果。

双树复小波变换的方法如下：

第一步，选取near_sym波作为复小波基函数的实部，qshift波作为复小波基 函数的虚部。

第二步，根据复小波基函数固定尺度函数。

第三步，利用复小波基函数的实部和虚部分别对图像重构的估计值进行尺度 分解，得到变换系数。

步骤4，判断图像重构当前估计值的继续条件是否满足，若是，执行步骤5， 否则，图像重构的迭代步长加倍，执行步骤3。本发明实施例中，图像重构当前 估计值的继续条件满足时，将迭代次数加1。

所述的图像重构当前估计值的继续条件是指，图像重构当前估计值的目标函 数值小于图像重构初始估计值的目标函数值。

按照下式，获取图像重构估计值的目标函数值：

$f\left(x\right)=\frac{1}{2}{\left|\right|\mathrm{Ax}-y\left|\right|}_1^2+{\left|\right|{\Psi }^{T}x\left|\right|}_1$

其中，f(x)表示图像重构估计值的目标函数值，x表示图像重构的估计值， y表示原始的混叠多光谱图像，A表示压缩光谱成像系统所采用的观测算子，λ 表示目标函数的平衡参数，Ψ^Tx表示图像重构估计值的全变差，表示二范数 的平方，||·||₁表示l₁范数。

步骤5，获取图像重构当前估计值的下一个估计值。

第一步，采用步骤3所述的降噪处理方法，对图像重构的当前估计值进行降 噪处理，得到图像重构当前估计值的降噪处理结果。

第二步，按照下式，获取图像重构当前估计值的下一个估计值：

x_t＝(1-α)·x_t-2+(α-β)·x_t-1+β·Γ(x_t-1)

其中，x_t表示图像重构当前估计值的下一个估计值，x_t-2表示图像重构当前 估计值的上一个估计值、x_t-1表示图像重构的当前估计值，Γ(x_t-1)表示图像重构 当前估计值的降噪处理结果，α表示当前估计值的上一个估计值的迭代系数，β 表示当前估计值降噪处理结果的迭代系数。

按照下式，计算图像重构的本征系数：

$\rho =\frac{1-{\lambda }_N}{1+{\lambda }_N}$

其中，ρ表示图像重构的本征系数，λ_N表示图像重构的本征值，λ_N＝10^-4。

按照下式，计算当前估计值的上一个估计值的迭代系数：

$\alpha =\frac{2}{1+\sqrt{1-{\rho }^2}}$

其中，α表示当前估计值的上一个估计值的迭代系数，ρ表示图像重构的本 征系数。

按照下式，计算当前估计值降噪处理结果的迭代系数：

$\beta =\frac{2\alpha }{1+{\lambda }_N}$

其中，α表示当前估计值的上一个估计值的迭代系数，β表示当前估计值降 噪处理结果的迭代系数，λ_N表示图像重构的本征值，λ_N＝10^-4。

在本发明实施例中，执行步骤5时，首先判断图像重构的迭代次数是否大于 图像重构的迭代终止次数，若是，则结束，此时图像重构的结果为图像重构的当 前估计值，否则，继续执行步骤5。

步骤6，判断图像重构当前估计值的下一个估计值的继续条件是否满足，若 是，则执行步骤7，否则，用图像重构的当前估计值更新图像重构的初始估计值， 执行步骤3。

所述的图像重构当前估计值的下一个估计值的继续条件是指，图像重构当前 估计值的下一个估计值的目标函数值小于图像重构当前估计值的目标函数值。

步骤7，更新估计值。

第一步，用图像重构的当前估计值更新图像重构当前估计值的上一个估计 值。

第二步，用图像重构当前估计值的下一个估计值更新图像重构的当前估计 值。

步骤8，判断图像重构的终止条件是否满足，若是，则结束，此时图像重构 的结果为图像重构的当前估计值，否则，执行步骤5。本发明实施例中，若图像 重构的终止条件不满足，则图像重构的迭代次数加1。

所述的图像重构的终止条件是指，终止函数值小于图像重构的迭代终止系 数。

按照下式获取终止函数值：

$C(x_t,x_{t-1})=\frac{|f\left(x_t\right)-f\left(x_{t-1}\right)|}{f\left(x_{t-1}\right)}$

其中，x_t表示图像重构的当前估计值，x_t-1表示图像重构当前估计值的上一 个估计值，t表示图像重构的当前迭代次数，C(x_t,x_t-1)表示终止函数值，f(x_t)表 示当前估计值的目标函数值，f(x_t-1)表示当前估计值的上一个估计值的目标函数 值，|·|表示取绝对值操作。

下面结合仿真图对本发明的效果做进一步的说明。

1.仿真条件：

本发明仿真实验的硬件测试平台是：Intel Core i73770K CPU，主频3.50GHz， 内存16GB；软件仿真平台为：windows 7 64位操作系统和Matlab 2013a；测试 图像为：哥伦比亚大学公开的光谱图像，空间分辨率为(512,512)，谱间分辨率 为31。

2.仿真内容与结果分析：

为了验证本发明的有效性，采用本发明做了两个仿真实验，根据仿真实验的 重构结果计算出多光谱图像重构结果各个谱段的峰值信噪比，并与全变差去噪方 法的重构结果进行比较。

仿真实验1，以哥伦比亚大学公开的光谱图像balloons图作为原始多光谱图 像，分别利用现有技术的全变差去噪方法和本发明进行多光谱图像重构。本仿真 实验对所有谱段的图像进行了重构，由于多光谱图像的谱段数较多，故本仿真结 果只列举了谱段7和谱段23的图像。多光谱图像重构结果的部分结果示意图如 图2所示。图2(a)为balloons图的原始多光谱图像谱段7的图像。图2(b) 为balloons图的原始多光谱图像谱段23的图像。图2(c)为全变差去噪方法的 多光谱图像重构结果谱段7的图像。图2(d)为全变差去噪方法的多光谱图像 重构结果谱段23的图像。图2(e)为本发明的多光谱图像重构结果谱段7的图 像。图2(f)为本发明的多光谱图像重构结果谱段23的图像。多光谱图像重构 结果各个谱段的峰值信噪比如表2所示。

仿真实验2，以哥伦比亚大学公开的光谱图像fake_and_real_food图作为原 始多光谱图像，分别利用现有技术的全变差去噪方法和本发明进行多光谱图像重 构。本仿真实验对所有谱段的图像进行了重构，由于多光谱图像的谱段数较多， 故本仿真结果只列举了谱段23和谱段31的图像。多光谱图像重构结果的部分结 果示意图如图3所示。图3(a)为fake_and_real_food图的原始多光谱图像谱段 23的图像。图3(b)为fake_and_real_food图的原始多光谱图像谱段31的图像。 图3(c)为全变差去噪方法的多光谱图像重构结果谱段23的图像。图3(d)为 全变差去噪方法的多光谱图像重构结果谱段31的图像。图3(e)为本发明的多 光谱图像重构结果谱段23的图像。图3(f)为本发明的多光谱图像重构结果谱 段31的图像。多光谱图像重构结果各个谱段的峰值信噪比如表3所示。

本发明与全变差去噪方法重构多光谱图像所消耗的时间如表1所示，从表1 中可以看出，本发明相比于全变差去噪方法，多光谱图像重构所消耗的时间缩短 了近一半，表明本发明的多光谱图像的重构速度较快。

表1重构耗时

耗时/s balloons图 ake_and_real_food图 全变差去噪方法 2310.5098 2430.0513 本发明 1479.1343 1206.4893

从表2和表3中可以看出，本发明相比于全变差去噪方法，多光谱图像重构 结果单个谱段的峰值信噪比有2-5dB的提高，多光谱图像重构结果的平均峰值信 噪比有2-3dB的提高，表明本发明的多光谱图像的重构效果较好这两方面充分说 明，本发明在进行多光谱图像重构时有更好的性能。

表2 balloons图重构结果各个谱段的峰值信噪比

谱段数 全变差去噪 本发明 差值 1 32.4802 38.1640 5.6839 2 31.8843 36.2071 4.3228 3 30.6707 34.9792 4.3084 4 30.4354 34.5914 4.1560 5 30.3348 34.3710 4.0361 6 30.3945 34.5187 4.1241 7 30.7402 34.8957 4.1555 8 31.1165 35.0360 3.9195 9 31.5732 35.2971 3.7239 10 32.1491 35.4428 3.2937 11 32.1298 34.8602 2.7304 12 32.2287 34.5689 2.3402 13 32.2382 34.4358 2.1976 14 32.2104 34.3210 2.1105 15 32.3278 34.3828 2.0550 16 32.4797 34.4873 2.0076 17 32.5129 34.4550 1.9421 18 32.4059 34.4290 2.0231 19 32.3837 34.3540 1.9703 20 32.5323 34.4901 1.9578 21 32.1909 34.4643 2.2735 22 31.7861 34.3561 2.5700 23 31.3612 34.1240 2.7627 24 31.1480 34.0563 2.9082 25 31.1071 34.1972 3.0901 26 30.8021 34.0592 3.2572 27 30.4019 33.8273 3.4254 28 30.1901 33.7249 3.5347 29 29.7807 33.4629 3.6822 30 29.6683 33.5733 3.9050 31 29.9386 33.8465 3.9080 平均 31.4066 34.5800 3.1734

表3 fake_and_real_food图重构结果各个谱段的峰值信噪比

谱段数 全变差去噪 本发明 差值 1 32.2512 35.2615 3.0103 2 30.6533 34.0775 3.4242 3 30.2628 33.4740 3.2112 4 29.9311 33.0701 3.1390 5 30.0304 32.9341 2.9036 6 30.1855 33.0100 2.8244 7 30.0668 32.8590 2.7922 8 29.9132 32.6567 2.7435 9 29.8053 32.3959 2.5906 10 29.9311 32.4361 2.5050 11 30.1413 32.4042 2.2630 12 30.2514 32.1753 1.9239 13 30.4493 32.0065 1.5572 14 30.6509 31.9729 1.3220 15 30.7311 31.9976 1.2665 16 30.8080 32.1297 1.3217 17 30.6843 32.0176 1.3334 18 30.5011 31.9104 1.4093 19 30.3611 31.8537 1.4926 20 30.2193 31.7772 1.5580 21 29.9408 31.5620 1.6212 22 29.5325 31.3607 1.8282 23 29.0162 31.2076 2.1914 24 28.7171 31.1013 2.3842 25 28.4108 31.0069 2.5961 26 28.1436 30.9256 2.7820 27 27.9722 30.8966 2.9244 28 27.6832 30.8091 3.1259 29 27.3315 30.6942 3.3627 30 27.0453 30.6409 3.5957 31 27.2524 30.6193 3.3669 平均 29.6411 32.0401 2.3990