非下采样剪切波变换域下遥感图像融合

熊璋玺， 刘明亮， 郭擎， 李安

(1.黑龙江大学 电子工程学院,黑龙江 哈尔滨 150080； 2.中国科学院空天信息研究院， 北京 100101； 3.黑龙江省信息融合估计与检测重点实验室，黑龙江 哈尔滨 150080)

遥感图像融合是将全色图像的高分辨率空间细节信息和多光谱或高光谱图像中丰富的光谱信息进行融合，以获得空间和光谱质量都高的图像。其算法可以分成成分替换法、多分辨率分析法、基于模型的方法3类。成分替换法是对多光谱图像进行光谱变换，目的是分离出相对于光谱对应物的空间信息。然后，用全色图像代替该空间分量，得到空间增强的多光谱图像。成分替换法主要包括IHS变换[1]、主成分分析方法(principal component analysis,PCA)[2]、GS方法[3]、自适应GS法[4]、基于物理约束的波段相关空间细节(band-dependent spatial-detail with physical constrains, BDSD_PC)方法[5]、基于部分置换的自适应分量替换(partial replacement-based adaptive component substitution, PRACS)法[6]等。成分替换法具有空间细节清晰、运行效率高、对误配准和混叠错误具有一定鲁棒性而被广泛使用，不过往往会带来光谱失真。

多分辨率分析方法是将全色图像多分辨率分解得到的空间细节注入到上采样的多光谱图像中，该方法通常通过线性分解方法获得，如小波变换[7]、超小波变换[8]、加性小波亮度比例(additive wavelet luminance proportional, AWLP)方法[9]、基于调制传递函数的广义拉普拉斯金字塔(modulation transfer function based generalized laplacian pyramid, MTF_GLP)方法[10]等。与成分替换方法相比，多分辨率分析方法能更好地保持光谱特性，但在空间增强方面有时令人不太满意。

基于模型的方法一般分为基于稀疏表示的方法和基于深度学习的方法。基于稀疏表示的方法[11]首先从低空间分辨率数据中学习光谱字典，然后结合已知的高空间分辨率数据预测高空间分辨率和高光谱分辨率数据。例如，Li等[12]提出了一种基于稀疏诱导先验信息的压缩感知方法，通过构建从MS图像中随机抽样的图像块字典来实现稀疏性。为了避免字典构建的成本，Zhu等[13]提出了一种稀疏图像融合算法。Cheng等[14]提出了一种融合小波变换和稀疏表示的融合框架。与多分辨分析法相比，这些方法具有超分辨能力和鲁棒性，能够获得更高的空间分辨率和光谱分辨率，且光谱失真较小。

近年来，人们对深度学习融合算法越来越感兴趣，例如Giuseppe等[15]提出一种名为基本卷积神经网络的全色图像锐化(CNN-based Pan-sharpening，PNN)的网络，它直接学习低分辨率图像和高分辨率图像之间的关系，冗余度高，学习会比较困难。Rao等[16]提出了一种基于残差卷积神经网络的图像融合方法，直接学习输入与输出之间的残差，但浅层网络很难学习到深层次丰富的地物特征，数据量增多时还会出现欠拟合。Li等[17]使用超分辨率网络来做融合，然而缺乏真实高空间分辨率多光谱图像做标签，只能通过降尺度方法制作模拟数据集。Xiong等[18]通过将质量评价函数应用于损失函数解决了制作模拟数据集问题，并为深度学习做遥感图像融合提供了新思路[19-20]。

虽然现有的融合算法在许多方面表现良好，但仍有一些方面需要改进。例如基于深度学习的方法往往需要非常大的训练集，而这些专业的遥感数据训练集通常很少。而且，不同的卫星有不同的数据类型，目前还很难做到不同卫星数据的同时训练。另外，网络训练需要大量的时间，使得网络调整的实时性较差[21-22]。在基于稀疏表示的融合算法中，如何找到一个最优的变换基并使在其上的表示最稀疏是较困难的问题。此外，稀疏表示有时会忽略图像的内在几何结构[23-24]。

针对成分替换和多分辨率分析方法的优点，本文提出了一种混合的遥感图像融合新方法。在NSST域，首先利用IHS变换提取多光谱图像的亮度分量，然后将全色图像和提取的亮度分量在NSST域分解成高频子带和低频子带，采用导向滤波器对高频子带进行融合，对低频子带进行区域清晰度加权融合。最后将融合后的高频子带和低频子带在进行NSST逆变换和IHS逆变换得到融合的高分辨率多光谱图像。为了验证方法的有效性，本文选取Quick Bird, SPOT&TM 以及 WorldView-2卫星图像进行实验验证。融合结果与6种典型先进方法比较，本文提出的融合方法在有参考和无参考质量评价指标上都表现优良。

1 非下采样剪切波变换与融合规则

1.1 非下采样剪切波变换

由Easley等[25]提出的NSST变换是小波变换在多维和多向情况下的扩展，它将非下采样金字塔(non-subsampled pyramid，NSP)变换与几种不同的剪切滤波器(shearlet filter, SF)相结合，以提供多尺度和多方向的分解。

NSST变换主要分为以下两步：多尺度分解和方向局部化。多尺度分解采用非下采样金字塔滤波器组(non-subsampled pyramid filter banks, NSPFs)实现。用剪切滤波器组(shearlet filter banks, SFBs)实现方向局部化。对源图像进行NSPFs分解，得到由细到粗的低频和高频子带图像。图1示出了两级NSST分解的示意图。

图1 两级NSST分解示意Fig.1 The schematic of a two level NSST decomposition

1.2 导向滤波

导向滤波器[26](guided filter, GF)是一种图像滤波技术，该技术通过利用一幅引导图x，对目标图像p(输入图像)进行滤波操作，从而使得输出图像与p基本相似，但纹理特征与x一致。本文将导向滤波器应用于高频子带的融合。

在本文中x为全色图像经NSST变换后得到的高频子带，p为多光谱图像的亮度分量经NSST变换后得到的高频子带，融合的高频子带为导向滤波的输出q。导向滤波的一个重要假设是输出图像q和引导图像x在滤波窗口ωk上存在局部线性关系：

qi=akxi+bk,∀i∈ωk

(1)

qi=pi-ni

(2)

最终的目标就是最小化这个噪声。对于每一个滤波窗口，该算法在最小二乘意义上的最优化可表示为:

(3)

式中ε是给定的正则化参数。系数ak和bk可通过线性回归直接求解：

(4)

(5)

(6)

导向滤波器具有许多优良的特性。首先，导频滤波器是局部线性滤波器，如式(1)所示，输出图像q可以理解为在局部范围内具有权重系数ak和偏移量bk的向导图像x之和，这使得向导图像的细节信息传输到输出图像成为可能。导向滤波器的这一结构传递特性可以最大限度地消除上采样带来的边缘遮挡伪影，从而提高融合图像的清晰度。其次，导向滤波器具有梯度保持特性，可以有效避免梯度翻转现象，这种特性可以应用于图像增强和高动态范围压缩。最后，导向滤波器模型具有快速、高效、易于实现的特点。基于这些特性，本文选择导向滤波作为高频融合规则。

1.3 融合规则

图2示出了所提出融合方法的示意图。在融合之前，源图像已经被配准。

图2 本文方法流程Fig.2 Flow char of proposed method

影响多光谱与全色图像融合质量的主要因素一方面是全色图像空间细节的融入度，另一方面是多光谱图像的保持度。源图像经NSST变换被分解为不同的频率分量，低频分量部分是近似源图的子图，其表示图像强度或灰度值变化缓慢的区域，即图像中较大的平坦区域，它描述了图像的主要部分并决定其轮廓。低频融合规则的选取对于提升融合图像的视觉效果，很好的保持源图像光谱信息具有十分重要的意义。高频分量对应于图像中变化剧烈的部分，即图像的纹理、边缘(轮廓)、噪声或细节，是图像突变性的反映。高频融合规则的选取对于空间细节的保留具有十分重要的意义。高频与低频的融合规则是本文融合方法中核心，所以下文着重对融合规则进行描述。

1.3.1 低频分量融合规则

MS和Pan图像低频子带融合的目的是为了获得更清晰的MS图像。清晰图像和模糊图像之间的区别是我们需要研究的。从逻辑上考虑，图像模糊是因为图像中物体的轮廓不明显，轮廓边缘灰度变化不强烈，层次感不强造成的。反之，如果轮廓边缘灰度变化明显，层次性强，则图像清晰。对于一幅图像，可以把图像看作一个二维函数f(x,y)，由微积分可知，图像的清晰度可用图像的梯度来表示。

基于此，对于低频分量，设计基于区域清晰度加权的融合策略。梯度能量(energy of gradient, EOG)作为清晰度的评价指标，在空间域具有明确的物理意义，具有良好的时间性能和较高的灵敏度。EOG使用相邻点之间的差来计算一个点的梯度：

[f(x,y+1)-f(x,y)]2}

(7)

用CPL、CIL和CFL分别表示Pan、MS的亮度分量I和融合图像F在同一子带、相同方向和相同位置的低频系数，RPL和RIL分别表示Pan和MS的亮度分量I在相同子带、相同方向和相同位置的NSST域低频子带系数的EOG，则融合策略可以表示为：

CFL=wPCPL+wICIL

(8)

1.3.2 高频分量融合规则

对于高频子带融合，其目的是尽可能多地保留空间细节、边缘和轮廓。在第1.2节中讨论了导向滤波器的特性，这使得它作为一种高频子带融合策略成为可能。由于Pan图像比MS图像更清晰，具有丰富的空间细节和锐利的边缘特征，因此，将MS的亮度分量经NSST分解后的高频子带作为输入图像，Pan经过NSST分解后的高频子带作为向导图像，输出图像即为融合后的高频子带图像。此外，输入图像和向导图像必须保证在同一层和同一方向上。融合策略可以表示为：

(9)

2 融合实验与对比分析

为了验证本文方法的有效性与可靠性，用真实遥感图像进行实验验证。对同源、异源、高分辨率、中分辨率遥感图像进行实验，并用主观视觉和客观定量评价进行结果分析，设置算法参数后，对实验结果进行展示和讨论。

2.1 实验设置

2.1.1 实验图像

为了验证该方法的有效性和通用性，选取了3组不同数据类型的图像。第1组是QuickBird卫星拍摄的空间分辨率分别为0.7 m和2.8 m的Pan图像和MS图像。第2组是空间分辨率为10 m的SPOT的 Pan图像和空间分辨率为30 m的TM图像。第3组是worldwiew-2卫星拍摄的空间分辨率分别为0.5 m和2 m的Pan图像和MS图像。其中，第1组和第3组是同源卫星图像，第2组是异源卫星图像。图像地物场景包括城市、农田和草地。

2.1.2 所选对比方法

将所提出的融合方法与GSA、BDSD_PC、PRACS、AWLP、MTF_GLP以及基于卷积神经网络的方法PNN等6种有代表性的融合方法进行比较。其中GSA、BDSD_PC、PRACS是分量替换类，AWLP与MTF_GLP是多分辨率分析类，PNN是深度学习类。

2.1.3 客观评价指标

客观评价通过定量评价方法和指标，科学、客观地评价图像融合算法的性能，它不受人的视觉特征和心理状态的影响。在本文中，信息熵(information entropy, EN)、标准差(standard deviation, STD)、相关系数(correlation coefficient, CC)、结构相似度[27](structural similarity index, SSIM)、相对整体维数合成误差[28](erreur relative globale adimensionnelle de synthese, ERGAS)、光谱角匹配指数[29](spectral angle mapper, SAM)、无参考质量评价指标[30](quality with no reference index, QNR,Ds,Dλ)这9个被广泛认可的融合度量指标被应用到本文的实验中。其中EN与STD越大越好，CC、SSIM、ERGAS、SAM、QNR、Ds和Dλ的理想值分别为1, 1, 0, 0, 1, 0, 0。其中CC、ERGAS、SAM和Dλ度量融合图像的光谱质量，EN、STD、SSIM和Ds反映融合图像的空间细节，QNR是综合评价指标。

2.1.4 算法参数设置

本文提出的方法有3个参数需要设置：NSST分解层数L、滤波器半径r和导向滤波器的模糊度ε。实验按照经验设置L=4，r=1，ε=0.0012。

2.2 实验结果与分析

图3(a)、4(a)和5(a)分别显示QuickBird、SPOT&TM和WorldView-2卫星的上采样低分辨率MS图像。图3(b)～(h)、4(b)～(h)和5(b)～(h)分别给出了AWLP，BDSD_PC，GSA，MTF_GLP，PRACS，PNN和本文提出方法的融合结果。表1～3分别列出了图3～5融合结果的客观质量评价。每个客观质量评价指标的最佳结果都以粗体标记。为了更好地比较各种方法的融合结果，还对融合后的图像进行局部放大。

在图3中，通过将融合结果与原始图像进行视觉比较发现：BDSD_PC方法和GSA方法都发生了一定程度的光谱失真，特别表现在局部放大图像中，其树看起来比MS亮，除此之外，在融合图像左下角的一片灌木区光谱失真也严重。PRACS方法的融合图像右下角的屋顶出现严重光谱失真，其颜色变成了暗棕色。对于生长在河岸上的树木，可以发现AWLP方法和MTF_GLP方法的光谱失真很小，但在右下角的屋顶出现了较严重的光谱失真，与MS的屋顶对比，这2种方法屋顶看起来偏亮。虽然PNN方法保持了较好的光谱特性，但其整体空间特性不如本文的方法。本文提出的方法改善空间细节信息的同时保持了更好的光谱特性。

图3 第1组QuickBird数据不同方法的融合结果(RGB真彩色合成)Fig.3 The R, G, B bands for natural color composition by different methods of the first group QuickBird data

从表1可以看出，本文提出的方法在指标SSIM、CC、ERGAS、SAM、Ds和QNR上表现最好，EN和Dλ排名第3，STD排名第2，但只比最优的结果差一点。

表1 第1组QuickBird数据融合结果的客观评价Table 1 Objective evaluation of fusion results using the first group QuickBird data

在图4的局部放大图像中，BDSD_PC、PRACS和PNN方法看起来很模糊，尤其是PNN方法，空间细节没有太大的增强。GSA方法保持了良好的空间细节，但是光谱失真严重，从图像的左中部可以看到，绿色看起来更亮。与本文提出的方法相比，AWLP方法和MTF_GLP方法的差别不明显。

图4 第2组TM&SPOT数据不同方法的融合结果(SWIR, B, NIR假彩色合成)Fig.4 The SWIR, B, NIR bands for false color composition by different methods of the second group TM and SPOT data

从表2可以看出，本文提出的方法在指标EN、STD、Ds和QNR方面表现最好。SSIM、CC、ERGAS、SAM和Dλ排名第2，这些指标表现最优的方法是PNN方法，然而PNN方法的空间细节不如本文提出的方法好。

表2 第2组TM与SPOT数据融合结果的客观评价Table 2 Evaluation of fusion results using the second group TM and SPOT data

从图5的局部放大图像可以看出，BDSD_-PC方法、PRACS方法和PNN方法丢失了大量的空间细节。与本文提出的方法相比，GSA方法在图像的右上角看起来更明亮，而AWLP和MTF_GLP方法的纹理特征没有本文方法丰富。

图5 第3组WorldView-2数据不同方法的融合结果(RGB真彩色合成)Fig.5 The R, G, B bands for natural color composition by different methods of third group WorldView-2 data

表3中，本文提出的方法在EN、STD、SAM、Ds和QNR方面表现最好。SSIM和CC接近最优值1，ERGAS和Dλ值都较小。

表3 第3组WorldView-2数据融合结果的客观评价Table 3 Evaluation of fusion results using the third group WorldView-2 data

结合视觉融合结果和质量评价指标，可以看出本文提出的方法优于目前常用的对比方法。此外，在3组实验中，本文提出的方法在大多数质量评价指标上都优于其他方法，在保持原始多光谱图像光谱信息的同时提高了空间分辨率。

3 结论

1)本文方法在对源图像进行NSST变换的基础上，针对高低频子带信息的不同特点设计了基于区域清晰度加权的低频子带融合规则，以及基于导向滤波的高频子带融合规则。

2)利用3组高、中空间分辨率的遥感卫星图像进行同源、异源融合实验，融合结果对空间细节信息增强的同时更好地保留了光谱信息，验证了该方法的有效性。与其他6种有代表性的融合方法通过九种评价指标进行比较和评价，结果表明本文方法对于融合多光谱图像的光谱信息和全色图像的空间信息的有效性和优越性，是一种行之有效的多光谱图像与全色图像融合方法。