CN106251320B - 基于联合稀疏与结构字典的遥感图像融合方法 - Google Patents

Abstract

本发明公开了基于联合稀疏与结构字典的遥感图像融合方法，首先通过在线稀疏字典算法分别从全色图像和多光谱图像学习获得相应的结构字典；然后采用联合稀疏表示获取全色图像不同于多光谱图像的特有成分；最后利用ARSIS融合框架将全色图像的特有成分注入到多光谱图像，得到高分辨率多光谱图像；本发明有效利用字典的稀疏特性和原子间的相关性，进一步降低字典训练的复杂度，提高结构字典的自适应性，使得特有分量的重构更为准确，提高图像融合的质量；同时，本发明充分考虑全色图像的细节信息和低频信息，使得融合更为全面有效。

Description

基于联合稀疏与结构字典的遥感图像融合方法

技术领域

本发明属于图像处理技术领域，具体涉及基于联合稀疏与结构字典的遥感图像融合方法。

背景技术

遥感卫星在拍摄获取多光谱(MS)图像的同时，可获取同场景的全色(PAN)图像，其中，多光谱图像富含光谱信息，但空间分辨率低，清晰度差；全色图像的空间分辨率高，但光谱分辨率低；因此，如何充分利用多光谱图像的光谱信息和全色图像的空间信息获取高分辨率多光谱图像，成为遥感图像融合研究的热点话题。目前，卫星遥感图像的融合方法大致可归为四类：基于替代的方法、基于ARSIS思想的方法、基于多分辨率分析的方法和基于稀疏表示的方法。

分量替代方法最大的特点是计算效率高，能保持空间细节信息，但很难避免出现光谱扭曲现象；基于ARSIS思想的方法通常能够很好的保持多光谱图像的光谱信息，但不足之处是混叠以及局部差异会导致融合质量的下降；而基于多尺度分析方法能更好地保持光谱信息，但是会产生模糊或者信息扭曲，导致空间信息丢失，同时对图像所有几何特征表示不全面。

稀疏性是近些年提出的一种新型的先验假设，它符合人类的视觉特性，且稀疏表示能充分利用图像信号的内部结构信息。最近几年，在字典学习的基础上，越来越多的融合方法被提出。现有方法中有利用稀疏表示和ARSIS思想相结合的融合方法，该类方法首先利用稀疏表示从全色图像获取细节信息，然后通过ARSIS融合框架将细节信息注入到光谱图像各波段中，这类方法的融合效果虽得到极大地提高，但存在的不足是：只融合了全色图像的细节信息而忽略了低频信息。

为充分利用图像内部及图像间存在的结构特性和相关性，以降低图像处理的工作量和提高图像融合的准确度，本文提出基于联合稀疏与结构字典的融合方法(JSRDIP)；一方面，充分考虑全色图像的细节信息和低频信息，使得融合更为全面有效，再者充分利用图像存在的结构特性，以进一步提高图像融合的质量；另一方面，改善遥感图像融合方法普遍存在的效率低的问题。

发明内容

针对上述现有技术中存在的缺陷或不足，本发明的目的在于，提供基于联合稀疏与结构字典的遥感图像融合方法，解决现有技术中遥感图像融合不全面、融合质量和融合效率低的问题。

为了实现上述任务，本发明采用如下技术方案予以实现：

基于联合稀疏与结构字典的遥感图像融合方法，包括以下步骤：

步骤一：获取多光谱图像和全色图像，对多光谱图像进行上采样处理得到上采样图像；

步骤二：对步骤一所得的上采样图像中各波段图像求均值得到强度分量，将全色图像与强度分量进行匹配得到新的全色图像；

步骤三：采用在线稀疏字典学习算法对步骤一得到的上采样图像中各波段图像以及步骤二得到的新的全色图像分别训练，获取对应的结构字典；

步骤四：采用联合稀疏模型，从步骤二得到的新的全色图像中提取不同于步骤一得到的上采样图像中各波段图像的特有分量；

步骤五：基于ARSIS融合框架，将步骤四得到的特有分量融入到步骤一得到的上采样图像的各波段图像中，得到融合多光谱图像。

本发明还有如下区别技术特征：

进一步的，步骤一中所述对多光谱图像进行上采样处理得到上采样图像是指：

采用双三次线性插入方法将多光谱图像上采样p倍得到上采样图像X_MS，其中，p为全色图像和多光谱图像空间分辨率的比值。

更进一步的，p＝4。

进一步的，步骤二的具体步骤包括：

步骤2.1：对步骤一所得的上采样图像X_MS中各波段图像求均值，得到强度分量I，即

其中，是X_MS第j(j＝1,2,…,k)个波段的光谱图像，k是波段的个数；

步骤2.2：将步骤一所得的全色图像与步骤2.1所得的强度分量I进行直方图匹配，得到新的全色图像X_PAN。

进一步的，步骤四的具体步骤包括：

步骤4.1：对上采样图像X_MS中各波段图像和X_PAN图像利用联合稀疏建模得到如下式所示的优化问题：

利用OMP算法对上式求解，得到稀疏系数其中，α^Cj为公共稀疏系数，为中第j个波段图像对应的特有稀疏系数，为全色图像X_PAN对应的特有稀疏系数；

步骤4.2：采用类似绝对值取大的方法对步骤4.1获得的全色图像X_PAN对应的特有稀疏系数进行选取，得到融合所需的全色图像特有稀疏系数α_pj，所述类似绝对值取大的方法如下式：

步骤4.3：将步骤4.2获得的全色图像特有稀疏系数α_pj重构，得到所需的特有分量计算公式如下：

其中，D_p为全色图像X_PAN经过训练所获取的结构字典。

进一步的，步骤五的具体步骤包括：

基于ARSIS融合框架，将步骤4.3获得的特有分量融入到步骤一所得的上采样图像X_MS的各波段图像中，得到融合多光谱图像，融入模型为：

其中，为融合后光谱图像的第j个波段图像。

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

(1)本发明采用联合稀疏模型来提取全色图像不同于多光谱图像各波段图像的稀疏特有分量，充分利用图像内部和图像间的相关性，进一步降低对图像处理的工作量，且提高了全色图像特有成分重构的精度；

(2)本发明采用在线稀疏字典学习方法对全色图像和多光谱图像各波段图像进行字典训练，有效利用字典的稀疏特性和原子间的相关性，进一步降低字典训练的复杂度，提高结构字典的自适应性；使得特有分量的重构更为准确，提高图像融合的质量；

(3)本发明在建模求解过程中充分考虑全色图像的细节信息和低频信息，使得融合更为全面有效。

附图说明

图1是本发明的融合框架图。

图2是模拟图像中不同融合方法对QuickBird卫星图像融合结果；其中，图2(a)为源多光谱图像，图2(b)为全色图像，图2(c)为参考图像，图2(d)为CS_BP方法融合图像，图2(e)为CS_OMP方法融合图像，图2(f)为DCS方法融合图像，图2(g)为SRDIP方法融合图像，图2(h)为DIP_OSDL方法融合图像，图2(i)为JSRDIP方法融合图像。

图3是模拟图像中不同融合方法对IKONOS卫星图像融合结果；其中，图3(a)为源多光谱图像，图3(b)为全色图像，图3(c)为参考图像，图3(d)为CS_BP方法融合图像，图3(e)为CS_OMP方法融合图像，图3(f)为DCS方法融合图像，图3(g)为SRDIP方法融合图像，图3(h)为DIP_OSDL方法融合图像，图3(i)为JSRDIP方法融合图像。

图4是实际图像中不同融合方法对QuickBird卫星图像融合结果；其中，图4(a)为源多光谱图像，图4(b)为全色图像，图4(c)为CS_BP方法融合图像，图4(d)为CS_OMP方法融合图像，图4(e)为DCS方法融合图像，图4(f)为SRDIP方法融合图像，图4(g)为DIP_OSDL方法融合图像，图4(h)为JSRDIP方法融合图像。

图5是实际图像中不同融合方法对IKONOS卫星图像融合结果；其中，图5(a)为源多光谱图像，图5(b)为全色图像，图5(c)为CS_BP方法融合图像，图5(d)为CS_OMP方法融合图像，图5(e)为DCS方法融合图像，图5(f)为SRDIP方法融合图像，图5(g)为DIP_OSDL方法融合图像，图5(h)为JSRDIP方法融合图像。

以下结合附图和实施例对本发明的具体内容作进一步详细地解释和说明。

具体实施方式

为了使本发明的目的、技术方案及优势更加清楚，结合附图及实施例对本发明做进一步详细说明。

结合图1所示的融合框架，本发明的具体实现包括如下步骤：

步骤一：获取多光谱图像和全色图像，对多光谱图像进行上采样处理得到上采样图像；使多光谱图像与全色图像的尺寸大小相同，以便后续处理；步骤一中所述对多光谱图像进行上采样处理得到上采样图像是指：

采用双三次线性插入方法将多光谱图像上采样p倍得到上采样图像X_MS，其中，p为全色图像和多光谱图像空间分辨率的比值，p＝4；

步骤二：对步骤一所得的上采样图像中各波段图像求均值得到强度分量，将全色图像与强度分量进行匹配得到新的全色图像；进行匹配时是以多光谱图像为标准，调节全色图像的色调与反差，以便做进一步的运算；

步骤2.1：对步骤一所得的X_MS中各波段图像求均值得到强度分量I，以便与全色图像做匹配，其中

其中，是X_MS第j(j＝1,2,…,k)个波段的光谱图像，k是波段的个数，k＝4；

步骤2.2：将步骤一所得的全色图像与步骤2.1所得的强度分量I进行直方图匹配，得到新的全色图像X_PAN；

步骤三：采用在线稀疏字典学习(online sparse dictionary learning，OSDL)算法对步骤一得到的上采样图像中各波段图像以及步骤二得到的全色图像分别训练，获取对应的结构字典；

利用双稀疏模型，采用OSDL算法对步骤一获得的上采样图像X_MS中各波段图像以及步骤二获得的全色图像X_PAN分别训练以获取各自对应的结构字典D_j(j＝1,2,…,k)和D_P。

为获取有效快捷且自适应性强的字典，训练过程中需在字典中添加一些特定约束项，双稀疏模型的约束项即为字典的稀疏特性；在双稀疏模型中，有效结构字典中的原子为基字典Φ中少数任意原子的线性组合所构成，因此，双稀疏模型的表达式为：

其中，A∈R^L×m和X∈R^m×n分别为稀疏字典和系数矩阵，x_i和a_j分别为稀疏系数和稀疏字典原子，p和k分别为稀疏系数和字典的稀疏度，基字典Φ∈R^n×L由剪切小波的Kronecker积构成，因此，双稀疏模型中稀疏字典学习的表达式如下：

采用在线稀疏字典学习算法(OSDL)对图像块进行训练，学习得到稀疏字典，再与基字典线性构成学习型结构字典；因在OSDL算法学习过程中采用的是批量处理替代原来的小图像块处理方式，从而减少了稀疏编码和稀疏字典更新时二者交替变换的次数，进一步降低计算复杂度，训练效率得以提高；再者，OSDL算法采用随机梯度下降法对稀疏字典中每个原子进行更新，可增强字典的自适应性，可使图像稀疏求解及重构更为准确，进一步提高融合图像的质量；

步骤四：采用联合稀疏模型(Joint Sparsity Model，JSM)，从步骤二得到的新的全色图像中提取不同于步骤一得到的上采样图像中各波段图像的特有分量；联合稀疏不仅可充分提取全色图像的特有分量，还可以进一步降低计算的复杂度，且提取出的特有分量较全面地包含多光谱图像所缺失的信息；

不同传感器可获得地物信息之间所存在的冗余信息和互补信息，融合处理就是利用冗余性和互补性来增加图像的信息量；在JSM-1模型中，所有图像都是由一个稀疏共同分量和一个稀疏特有分量组成，因此，对图像做融合处理时可转化为将稀疏特有分量融合到新图像中以尽可能的包含所有待融合图像的信息；

步骤4.1：对上采样图像X_MS中各波段图像和X_PAN图像利用联合稀疏建模得到：

其中，D_cj＝[D,D_p,D_j]，分析字典D由剪切小波(The Cropped Wavelets)构成；

令E＝blkdiag(D_j,D_p)Ψ＝[B,E]，则上式转化为如下式所示的优化问题：

利用OMP算法对上式求解，得到稀疏系数其中，α^Cj为公共稀疏系数，为中第j个波段图像对应的特有稀疏系数，为全色图像X_PAN对应的特有稀疏系数；联合稀疏求解得到的特有稀疏系数用来重构融合所需的全色图像的特有分量；

步骤4.2：由于全色图像的特有分量无需全部融入到多光谱图像各波段中去，因此选取最大活跃度对应的特有稀疏系数α_pj作为融合系数；为了更好地提取所需的特有稀疏系数，采用类似绝对值取大的方法对步骤4.1获得的全色图像X_PAN对应的特有稀疏系数进行选取，得到融合所需的全色图像特有稀疏系数α_pj，选取方式如下式所示：

步骤4.3：将步骤4.2获得的全色图像特有稀疏系数α_pj重构，得到所需的特有分量重构出的特有分量是为了融入到多光谱图像中去，以提高多光谱图像的分辨率，特有分量计算公式如下：

其中，D_p为全色图像X_PAN经过训练所获取的结构字典；

步骤五：基于ARSIS融合框架，将步骤四得到的特有分量融入到步骤一得到的上采样图像的各波段图像中，得到融合多光谱图像；其中，采用ARSIS融合框架是为了使融合图像很好地保留源图像的光谱信息；

基于ARSIS融合框架，将步骤4.3获得的特有分量融入到步骤一所得的上采样图像X_MS的各波段图像中，得到高分辨率多光谱图像，融入模型为：

其中，为融合后光谱图像的第j个波段图像。

实施例

本发明采用四组卫星遥感图像来验证所提融合算法的有效性；QuickBird卫星可提供空间分辨率为0.7米的全色图像和空间分辨率为2.8米的多光谱图像；IKONOS卫星提供空间分辨率为1米的全色图像和空间分辨率为4米的多光谱图像；其中，多光谱图像包括红、绿、蓝和近红外四个波段；为更好评估融合方法的实用性，本发明给出模拟图像实验和实际图像实验，其中，模拟图像实验中采用的模拟图像是通过对实际图像依次进行MTF滤波和下采样4倍获得，实际图像实验是直接将真实图像进行融合。

实验中采用已配准好的全色图像和多光谱图像，其大小分别为256×256和64×64，本发明主要与基于压缩感知(CS)、基于分布压缩感知(DCS)以及同类采用注入模型(SRDIP)的融合方法进行比较；一方面为使实验对比分析更准确，实验部分将给出原CS方法(CS_BP)以及CS改进方法(CS_OMP)作为对比实验；另一方面为突出结构字典的性能，增加DIP_OSDL融合方法作为对比试验。CS方法和DCS实验中，全色图像和多光谱图像对应的图像块大小分别是8×8和2×2，重叠像素值个数分别为1和4，实验重构误差ε设为1，其中，CS方法采用的字典为随机抽取的图像块组成，且字典大小为1000，而DCS所采用的字典包含K-SVD算法训练得到学习型字典和分析型字典；SRDIP融合方法实验中，全色图像的图像块大小为8×8，重叠像素值个数为7，字典大小为64，trous小波分解层数为2，稀疏表示求解误差为0.01，K-SVD训练字典中，稀疏编码的稀疏度为6，迭代次数为10；DIP_OSDL方法选用的全色图像的图像块大小为3×3，像素重叠个数为2，字典大小为64，trous小波分解层数为2，稀疏求解误差为0.01，采用OSDL算法训练字典时，稀疏编码的稀疏度为6，稀疏字典的稀疏度为40，迭代次数为3，其中，高频结构字典D_HF为64×128，低频结构字典D_IF大小为64×64；本文方法选取的全色图像块大小为4×4，像素重叠个数为3，稀疏求解允许误差为0.01，字典训练中稀疏字典的稀疏度为40，稀疏编码的稀疏度为6。

模拟图像实验效果分析：

图2(a)、(b)和图3(a)、(b)分别为QuickBird和IKONOS源多光谱图和的全色图，图2(c)和图3(c)为参考图像，图2(e)～(i)和图3(e)～(i)分别为QuickBird和IKONOS卫星融合结果图(只显示RGB三个波段)，每幅图像中的左下角为局部放大图；通过对原始图像和融合图像进行视觉比较，可以看出所有方法都能提高源多光谱图像的分辨率，但是CS方法融合图像的颜色发生较大改变；从图2(f)～(h)可以看出，其他方法均能很好地保留源多光谱图像的光谱信息，但这些方法的局部放大图像出现有块状效应；图3(g)和图3(h)可以看出，SRDIP和DIP_OSDL的融合图像虽然清晰度很高，但是部分地方融合过度，部分地方融合不到位；而从图2(i)和图3(i)可以看到，本发明在提高多光谱图像空间分辨率的同时，很好地保留了源图像的光谱信息，获得更好、更为自然的融合图像。

通过对图像的分析可以对融合结果有比较直观的认识，但单纯依靠主观评价很难给出正确判断，因此，需结合客观指标共同对融合结果进行评估；本发明采用SSIM、MI、CC、SAM、Q4这五种客观评价指标对图像进行全面地评估；其中，SSIM通过从亮度、对比度和结构这三个方面的比较，实现对融合结果中每个波段图像与参考图像中每个波段图像结构相似的度量；MI用以测量参考图像与融合图像的互相关信息；CC表示相关系数，用于测量两个图像对应的波段图像的光谱和空间信息相关程度；SAM表示全局光谱扭曲测量，最优值为0；Q4是一种综合性更强的融合图像质量标准，最优值为1。

表1和表3所示为不同融合方法下模拟图像融合结果的客观指标；从表1和表3可以看到，本发明的各方面客观指标较优于其他方法，其中，本发明的SAM值远远低于CS方法和DCS方法，这是因为本发明采用基于ARSIS融合框架，所以可以很好地保留源多光谱图像的光谱信息；同时，本发明的CC、MI以及SSIM值均远远高于CS方法，即本文方法融合图像与参考图像的相关性以及所含的互信息量都要远远高于CS方法。

相对于SRDIP和DIP_OSDL方法，本文方法也呈现出相对的优越性；表2和表4给出各种融合方法的计算复杂度，从表中可以看到，本发明的耗时远远低于其他方法，即本发明的计算效率有了明显的提高；综上，本发明在保证融合质量最佳的情况下，计算效率也得到极大提高。

表1 不同融合方法对QuickBird卫星获得的图像融合结果的客观指标

表2 对QuickBird卫星获得的图像进行融合所耗时间

表3 不同融合方法对IKONOS卫星获得的图像融合结果的客观指标

表4 对IKONOS卫星获得的图像进行融合所耗时间

实际图像实验效果分析：

图4(a)、(b)以及图5(a)、(b)分别表示QuickBird和IKONOS卫星图像的源多光谱图和全色图，图4(c)～(h)和图5(c)～(h)分别为QuickBird和IKONOS卫星融合结果图(只显示RGB三个波段)，每幅图左下角为融合图像的局部放大图；其中，图4(c)、(d)和图5(c)、(d)是CS方法的融合结果，可以看到，相比于其他方法，图4(c)、(d)和图5(c)、(d)的对比度较差，有少许光谱失真，且图4(d)从整体上看去有些模糊；从图4(e)和图5(e)整体来看，相比于DSC方法，本发明的融合图像的空间分辨率较高；对比图5(f)～(h)的局部放大区域可以看到，本发明得到的融合图像的轮廓较SRDIP方法和DIP_OSDL方法更稍清晰。

实际图像实验中，由于没有参考图像，为有效客观评估各融合方法的性能，采用无参考图像的客观评价指标QNR来评价图像融合质量；QNR包含光谱信息损失指标D_λ和空间信息损失指标D_s，其中，QNR的最优值为1，而D_λ和D_s的最优值为0。

表5和表7所示为不同融合方法下实际图像融合结果的客观指标；从表5和表7可以看到，本发明融合过程中产生的光谱信息损失最少，虽然在表5中DIP_OSDL方法的空间信息损失略低于本发明，但本发明的无参评价客观指标QNR较其他所有方法是最优的；而对表6和表8的中融合图像耗时进行分析，本发明的计算效率明显高于其他方法，且本发明的计算效率相比无字典训练的CS方法提高将近一倍。

综上所述，本发明在很好保留待融合图像光谱信息的同时，大大提高了融合图像的空间分辨率，且进一步降低融合过程的计算复杂度。

表5 不同融合方法对QuickBird卫星获得的图像融合结果的客观指标

表6 对QuickBird卫星获得的图像进行融合所耗时间

表7 不同融合方法对IKONOS卫星获得的图像融合结果的客观指标

表8 对IKONOS卫星获得的图像进行融合所耗时间

Claims (5)

1.基于联合稀疏与结构字典的遥感图像融合方法，其特征在于，包括以下步骤：

步骤一：获取多光谱图像和全色图像，对多光谱图像进行上采样处理得到上采样图像；

步骤二：对步骤一所得的上采样图像X_MS中各波段图像求均值得到强度分量，将全色图像与强度分量进行匹配得到新的全色图像X_PAN；

步骤三：采用在线稀疏字典学习算法对步骤一得到的上采样图像中各波段图像以及步骤二得到的新的全色图像分别训练，获取对应的结构字典；

步骤四：采用联合稀疏模型，从步骤二得到的新的全色图像中提取不同于步骤一得到的上采样图像中各波段图像的特有分量，具体包括：

步骤4.1：对上采样图像X_MS中各波段图像和X_PAN图像利用联合稀疏建模得到如下式所示的优化问题：

利用OMP算法对上式求解，得到稀疏系数其中，α^Cj为公共稀疏系数，为中第j个波段图像对应的特有稀疏系数，为全色图像X_PAN对应的特有稀疏系数；Ψ＝[B,E]，B＝[D_cj ^T D_cj ^T]^T，D_cj＝[D,D_p,D_j]；E＝blkdiag(D_j,D_p)；是X_MS第j个波段的光谱图像；D_P为全色图像X_PAN经过训练所获取的结构字典，D_j为上采样图像X_MS中各波段图像经过训练所获得的结构字典，j＝1,2,…,k；

步骤4.2：采用类似绝对值取大的方法对步骤4.1获得的全色图像X_PAN对应的特有稀疏系数进行选取，得到融合所需的全色图像特有稀疏系数α_pj，所述类似绝对值取大的方法如下式：

步骤4.3：将步骤4.2获得的全色图像特有稀疏系数α_pj重构，得到所需的特有分量计算公式如下：

其中，D_p为全色图像X_PAN经过训练所获取的结构字典；

步骤五：基于ARSIS融合框架，将步骤四得到的特有分量融入到步骤一得到的上采样图像的各波段图像中，得到融合多光谱图像。

2.如权利要求1所述的基于联合稀疏与结构字典的遥感图像融合方法，其特征在于，步骤一中所述对多光谱图像进行上采样处理得到上采样图像是指：

采用双三次线性插入方法将多光谱图像上采样p倍得到上采样图像X_MS，其中，p为全色图像和多光谱图像空间分辨率的比值。

3.如权利要求2所述的基于联合稀疏与结构字典的遥感图像融合方法，其特征在于，p＝4。

4.如权利要求1所述的基于联合稀疏与结构字典的遥感图像融合方法，其特征在于，步骤二的具体步骤包括：

步骤2.1：对步骤一所得的上采样图像X_MS中各波段图像求均值，得到强度分量I，即

其中，是X_MS第j个波段的光谱图像，j＝1,2,…,k，k是波段的个数；

步骤2.2：将步骤一所得的全色图像与步骤2.1所得的强度分量I进行直方图匹配，得到新的全色图像X_PAN。

5.如权利要求1所述的基于联合稀疏与结构字典的遥感图像融合方法，其特征在于，步骤五的具体步骤包括：

基于ARSIS融合框架，将步骤4.3获得的特有分量融入到步骤一所得的上采样图像X_MS的各波段图像中，得到融合多光谱图像，融入模型为：

其中，为融合后光谱图像的第j个波段图像；是X_MS第j个波段的光谱图像，j＝1,2,…,k，k是波段的个数。