Worldview 3全色与短波红外影像两步式融合框架

保 盈，王雷光，2，3，郭梦晓，代沁伶，郑 晨

1.西南林业大学 林学院，昆明650224 2.西南林业大学 大数据与人工智能研究院，昆明650224 3.西南林业大学 国家林业与草原局森林生态大数据重点实验室，昆明650224 4.西南林业大学 艺术与设计学院，昆明650224 5.河南大学 数学与统计学院，河南 开封475004

遥感影像像素级融合是利用算法将同一场景的两幅或多幅影像的互补性信息组合起来，形成一幅新影像的过程[1]，其主要目的是为了解决传感器设备获取遥感数据时，时间-空间-光谱分辨率之间的固有矛盾，获得同时具有更高时-空-谱分辨率的合成数据。像素级融合为影像地图制作、变化检测、目视解译及分类等应用提供了重要的数据资源[2]。

相关统计表明[3-4]，超过70%的光学对地观察卫星，同时负载高空间分辨率和多/高光谱分辨率传感器。因此，现有研究多围绕相同卫星平台同步获取的全色与光谱数据开展。典型数据源有QuickBird影像、Worldview-2影像以及国产高分一号、二号影像数据等。这类数据的特点是全色与多光谱波段的波谱范围重叠度高，波段间信息冗余度高，且空间分辨率相差较小，一般全色与多光谱波段的分辨率比率多为4倍。

围绕全色与多光谱影像的融合问题，现有融合方法主要有成分替代融合法（Component Substitution，CS）、多分辨率分析融合法（Multiresolution Analysis，MRA）、基于稀疏表达的方法和基于深度学习的融合方法等[5-6]。CS类方法和MRA类方法都是通过计算低分辨率全色分量与全色分量之间的差异提取出结构信息，然后通过合适的注入方案，将提取出的高空间结构信息注入到多光谱图像中。两类的差异在于低分辨率全色分量的计算方式[5]。前者一般采用多光谱波段的线性加权计算，后者多借助全色波段的滤波实现。

典型的CS方法包括PCA[7]（Principal Component Analysis，主成分变换法）、GS[8]（Gram-Schmidt Transform，GS变换法）和Aiazzi等人提出的GSA（Adaptive GS，自适应GS法）[9]等。典型的MRA方法包括HPF[10]（High Pass Filter，高通滤波融合法）、ATWT[11]（àtrous Wavelet Transform，小波变换法）等。这些方法在空间信息增强和光谱信息保持等方面各具优势，其相似之处在于事先给定融合规则[12]。基于稀疏表达的方法将融合问题变成稀疏表达下最优化复原问题，融合的影像可以从训练样本中重建，减少光谱畸变[13]。基于深度学习的方法通过深度神经网络的训练建立输入影像和融合结果的映射模型[14]，模型的建立需要大量数据的标记数据集。

短波红外（Shortwave Infrared，SWIR）属于红外数据，一般指波长范围在1 100～2 400 nm[15]的影像波段。红外遥感仍属被动、光学遥感范畴，因其成像波长较长，不仅能够反映这一波段范围内地物的独特反射特性，还具有大气散射作用小，穿透雾、烟和云层的能力强的特点。

现有的光学遥感卫星如Worldview 3、Landsat-7、Landsat-8、Sentinel-2，以及高分四号等都带有短波红外传感器。Worldview 3是其中一颗先进的、商业化应用的、超光谱、超高空间分辨率卫星。除了提供0.31 m空间分辨率、覆盖波谱范围为400～800 nm的全色影像（Panchromatic，Pan）和8波段1.2 m多光谱（Multispectral，MS）影像外，还提供了8个波谱覆盖总范围为1 195～2 365 nm的SWIR波段。但受入射电磁波能量强度、传感器辐射分辨率以及政策因素的影响，可用的短波红外波段的空间分辨率仅为7.5 m。全色与多光谱波段的空间分辨率比率为4，与SWIR波段空间分辨率比率为25。

Worldview 3全色影像和短波红外影像不仅波长范围完全不重叠，并且两者的空间分辨率差异较大。因此，相比传统的全色与多光谱影像的融合，短波红外影像的空间分辨率增强更具挑战性。有鉴于此，本文提出了一种两步式融合框架：通过渐进式的空间信息注入，旨在提高短波红外影像的空间分辨率，同时减少短波红外影像的光谱失真。

1 提出的融合框架与方法

按照融合的一般流程，首先需要将短波红外与全色影像进行配准，将短波红外影像重采样到全色影像相同大小，然后再将抽取的空间细节以一定方法到重采样后的低空间分辨率影像中。但是，由于本文中短波红外影像与全色影像两者分辨率差别较大，如果直接将短波红外影像与全色影像直接进行融合，受到SWIR高频细节的缺失，与全色波段亮度值范围、空间结构的不匹配等因素的影响，可能会导致融合时块状效应的产生和空间信息的损失，其具体融合的有效性仍有待评价。有鉴于此，本研究提出了如下两步式的融合框架。

首先对相关公式标识做如下约定：将Worldview 3卫星原始单波段全色影像和8个波段的短波红外影像分别记为Pan和，SWk为对应的第k个SWIR波段；采用和分别代表波段bd进行下采样、上采样获得的数据；和代表融合结果以及第i步的融合结果。

1.1 两步式融合框架

本框架定义的融合过程为：采用高空间分辨率的Pan和空间分辨率低的，通过两步融合，获得与全色影像空间分辨率一致、且保留与原始短波红外光谱特性的新影像。具体而言，如图1所示，该框架主要包括五个步骤：

（1）先对Pan进行下采样，使其空间分辨率由原始0.31 m降到3.6 m，结果记为。

上述融合框架中，第一步融合将SWIR波段的空间分辨率提升至3.6 m，第二步融合提升至0.31 m。通过两次空间细节的注入，在一定程度上缓解了SWIR波段低频信息与全色波段高频信息不匹配导致的融合块状效应，并有效增强了SWIR的空间分辨率。

提出的融合框架涉及两步融合处理，现有的像素级融合方法均可以作为处理的候选方法。下一节将介绍几种常用的成分替代和多分辨分析融合方法的原理及应用到本框架中的过程。

1.2 采用的融合方法

基于成分替代（CS）和多分辨分析（MRA）的融合方法多具有计算复杂度低的特点。考虑到这些方法在空间分辨率增强和光谱保真等方面各具优劣[16]，因此选择了6种典型的CS和MRA融合方法进行组合实验，以分析评价两步式框架的有效性。6种融合方法的选择综合考虑了算法的计算复杂度、算法对参数依赖程度、融合质量等因素。另外，由于本节只涉及单步融合方法，故公式中省略标识融合步数的上标。

1.2.1 成分替代融合方法

CS融合方法的思路是将多光谱影像投影变换到另一个向量空间，然后采用全色波段替代变换空间内的某个多光谱分量，最后再逆变换到原始图像空间。上述过程也可以表示为空间细节抽取和细节分波段注入的通用融合模型[17]，表示为：

其中，Fk为融合后的影像；为重采样后的SWIR影像；Pan为校正后的全色波段；gk为注入增益；IL为待估计的低分辨率全色波段，由式（2）计算：

ωk是待估计的权重系数。Pan-IL代表从全色波段中抽取的有效细节。

实现上，CS方法包括四个主要步骤：

（1）对SWIR影像进行上采样，使其与Pan具有相同的影像大小，记为。

（2）估计或给定ωk，然后根据式（2）对SWIR各波段进行加权求和，计算出强度分量IL。

（3）将Pan波段与计算出的强度分量IL进行直方图匹配，更新Pan波段，以增强二者的相关性，减少光谱失真，即：

（4）根据式（1），由Pan-IL提取细节信息，并按照一定的权重系数gk注入到SWIR各波段中，得到融合结果Fk。

不同CS融合方法主要体现在IL合成方式的差异和空间细节注入权重gk的计算差异。下面将介绍本文选取的典型的CS方法。

（1）PCA[7]法借助正交变换将具有一定相关性的原始变量转化成一组新的彼此之间互不相关的综合性变量，同时尽可能多地保留原有信息。将其应用于影像融合，理想结果是多光谱影像的空间信息（所有通道共享）集中在第一个成分（PC1），作为IL，故其权重ωk和注入增益gk都由变换后的第一个主成分求得。一般ωk为正变换矩阵的第一行，记为Il，gk是后向变换矩阵的第一列，记为Ih。

（2）GS[8]变换融合法是多元统计分析中常用的方法，它可以对矩阵或多维影像进行成交变换，消除冗余信息，同时保持良好的数值稳定性。GS变换法是PCA变换法的推广，两者区别在于PCA变换后第一主成分的信息量最多，其他成分的信息依次减少，而GS变换后，各分量只是正交，但信息量均匀分布，这种方法避免了PCA中信息过分集中的问题，对多光谱的各个波段使用相同的权重系数ωk，更好地保持了光谱信息。

（3）GSA[9]是对传统的GS方法的改进。与GS融合方法的融合过程相同，区别在于IL的计算方式不同。GS变换基于全局模型，对多光谱的各个波段使用相同的权重系数ωk，而GSA方法是根据均方根误差最小化原理，对影像的各个波段进行自适应权重系数的选择，从而更好地保持光谱连续性，即：

上述三种方法的光谱权值ωk和注入增益gk的计算方法具体如表1所示。

表1 6种融合方法的权重和增益确定方法

1.2.2 多分辨率分析融合方法

MRA方法也可以采用空间细节抽取和细节分波段注入的通用融合模型式（1）表达[5]，实现过程也与上文CS方法类似。两者的主要区别是：相比CS类方法，MRA类方法的低分辨率版本影像IL，往往借助小波变换，高、低通滤波器等多尺度信号分析工具，在不同分解尺度、不同频域子带上实现，往往具有更佳的光谱保真效果[16]。本文采用的MRA方法包括HPF[10]、SFIM[18]、ATWT[11]三种。

（1）HPF[10]是高通滤波器的简称，一幅影像中的高频分量对应空间细节信息，对影像进行高通滤波器操作可以滤除影像的低频信息，同时将高频细节信息保留下来。HPF融合方法的原理就是利用高通滤波器来提取Pan影像中的高频分量，再通过特定的权重模型将提取的空间高频信息注入到另一幅影像中，引入空间信息。

（2）SFIM[18]（Smoothing Filter-Based Intensity modulation，基于平滑滤波的亮度调制方法），基于太阳辐射和地表反射模型提出的，并用SPOT和TM数据进行了融合实验。与经典Brovey融合相比，SFIM能提高多光谱图像的空间细节且光谱保真度更好。

（3）ATWT[11]是一种基于多孔小波变换的影像融合方法。首先，对严格空间配准后的各原始影像分别进行多尺小波分解，生成对应不同频率子带的子图像；然后，建立适当的融合规则，分别对全色和光谱图像的高、低频子带进行融合；最后，进行小波逆变换，得到最终的融合影像。小波变换的优点是具有时域和频域的局部性、冗余信息较少。该方法将图像在空间频域内分解成几乎不相交的带通信道，特别适用于卫星图像的融合。

1.3 两步式融合过程

表1列出了上述两类、共6种融合算法的权重和注入增益的计算方法。这些方法的差异主要在于如何从Pan波段中提取空间细节，以及如何分配空间细节到对应的SWIR波段中。许多研究表明[16]，CS类的算法可以优化融合图像的空间细节质量，而MRA类的算法融合图像的光谱失真较低。而具体选用的这些方法都具有计算复杂度低、算法预设参数少的特点。

图1 全色与短波红外影像融合的两步框架

具 体 选 用 的PCA[7]、GS[8]、GSA[9]、HPF[10]、SFIM[18]、ATWT[11]，均需要以严格配准的Pan和SWIR波段为输入，输出结果为空间分辨率与输入Pan一致、波段数与输入SWIR一致的新SWIR波段。理论上，6种方法可以任意选择，带入到图1所示的框架内，作为（3）和（5）融合步骤的算法候选，从而形成该框架下共计36种融合设置。本研究将通过实验选择最佳的融合方法组合，进而验证提出框架的可行性。

2 实验与分析

2.1 实验数据与质量评价方法

为了验证该融合框架的有效性，本研究采用的实验数据为巴西里约热内卢的Worldview 3影像，成像时间为2016年2月5日。其中，全色影像的大小为13 632×11 244像素，空间分辨率为0.31 m，短波红外影像的分辨率为7.5 m。影像包含丰富的地物类型，有水体、建筑物、植被、道路等。为便于计算和结果展示，选取了两组1 024×1 024像素的子影像进行融合实验，实验区涵盖了上述主要地物类型。

为更客观地分析融合结果，采用反应影像亮度信息变换的标准差、反应影像清晰度的平均梯度和空间频率、反应影像信息量的信息熵以及广义无参考质量评价（Generalized Quality With No Reference，GQNR）[19]，进行定量评价。

（1）标准差（Standard Deviation，SD）可衡量图像信息的丰富程度，融合图像F的标准差定义为：

其中，Fˉ是图像F的像素均值，一般来讲，图像F的标准差越大，则图像内所含信息越丰富。若SD趋近于0，则图像各像素点的像素值在其均值上下几乎无波动，即图像被同一灰度级填充。

（2）空间频率（Spatial Frequency，SF）反映图像灰度的变化率，其计算公式为：

其中，RF和CF分别为图像F的行频率和列频率，即：

空间频率越高，图像越清晰，故空间频率可用于度量图像的清晰度。

（3）平均梯度（Average Gradient，AG）和空间频率一样，平均梯度也能反应影像的清晰度，其计算公式为：

其中，∆x F(i,j)=F(i,j)-F(i+1,j)为水平梯度，∆y F(i,j)=F(i,j)-F(i,j+1)为垂直梯度。

（4）信息熵（Information Entropy，IE）反应了影像的信息量，用于衡量融合图像中信息丰富程度，能量分布得越均匀，熵越大。在影像中，熵越大，表示影像包含的信息量越多，影像的融合效果越好。

F()

a是某一灰度等级a在图像中出现的频率。

（5）广义无参考质量评价指标（GQNR）[19]用于评价Worldview 3影像融合质量。该指标由光谱失真度(Dλ)和空间偏差度(DS)的乘积构成：

空间偏差DS较小，空间偏差越小，光谱失真度Dλ越小，光谱失真越小，GQNR越小越好。

光谱失真度(Dλ)定义为：

其中，SW是上采样后待融合SWIR影像，是对应融合结果；N为波段数；p为给定正整数，一般取1[20]；di,j度量融合前后i波段和j波段相关性的变化情况，即SW影像中i与j波段的相关性与影像中i与j波段的相关性的一致性程度：

且

其中，σxy表示样本x和y的协方差，和表示样本x和y的方差，和分别表示样本的均值[20]。Q函数的范围为[0，1]，Q越大代表个波段之间相关性越强。因此，Dλ表示融合结果和原始SWIR波段特征之间的光谱相似性，Dλ越小，光谱失真越小。

空间偏差度(DS)定义为：

其中，IL是与原始SWIR大小相同的低分辨率全色波段，q通常也设置为1[20]。当DS达到最小值（接近于0），说明融合前后的全色与短波红外波段的相关性接近相同。故DS值越小，影像空间质量越好。

2.2 实验设置

为了验证提出框架的有效性，将PCA[7]、GS[8]、GSA[9]、HPF[10]、SFIM[18]、ATWT[11]共6种方法代入两步式融合框架，形成36种组合配置（详见表2）。实验结果评价采用目视对比和定量评价结合的方式，首先分析两步式融合的必要性，然后对不同的两步融合方法进行对比和评价。

表2 两步式融合框架的36种融合配置

图2 三种融合方法直接用于全色与短波红外影像的结果

2.3 实验分析

2.3.1 两步式融合的必要性分析

首先，考虑直接采用0.31 m分辨率Pan和7.5 m的SWIR直接融合方法。图2和图3展示了两个典型区域的融合结果。其中，图2采用的融合方法为PCA[7]、GS[8]和GSA[9]，图3采用的融合方法为PCA[7]、GS[8]、GSA[9]、HPF[10]、SFIM[18]、ATWT[11]。

图3 6种融合方法直接用于全色与短波红外图像融合的结果

由图2可以看出，短波红外直接采用立方插值上采样到全色波段空间大小后（图2（b）），其空间细节未见明显增加，边界存在严重的马赛克现象，建筑物与植被区域存在块状效应和空间细节丢失。采用三种融合方法后（图2（c）～（e）），空间细节均有一定程度的增强，但光谱出现了扭曲，特别是图中与建筑物邻接的植被区域。

由图3可以看出：GS[8]、GSA[9]和PCA[7]方法的建筑物与植被区域存在块状效应和空间细节丢失，GSA[9]法存在严重的光谱失真与空间细节的丢失。HPF[10]、SFIM[18]、ATWT[11]三种融合方法空间信息引入较少，且块状效应明显（图3（f）～（h））。6种融合方法中，空间信息增强的目视效果最好的是GS[8]法和PCA[7]法，图3（i）还展示了两种方法的局部细节图对比。

由图2和图3可知：直接用6种单步式融合方法进行全色与短波红外影像融合，会造成严重的块状效应和空间信息的损失，不能达到提高短波红外影像空间分辨率的目的。这种问题的出现主要原因是Worldview 3的全色和SWIR波段两者空间分辨率的差异过大。直接将短波红外影像的分辨率重采样到与全色影像同样大小后，由于SWIR高频细节的缺失，以及与全色波段亮度值范围、空间结构的不匹配，导致融合时块状效应的产生和空间信息的损失。因此，直接融合的方法不适用于Worldview 3影像的全色与SWIR影像融合。

2.3.2 两步式融合方法的对比分析

考虑的36种融合配置中，部分融合结果明显视觉效果较差，故首先将其筛选出去，不参与后续的定量评价。

如图4（a）～（e）所示，两步融合均采用MRA类方法的9种配置（对应表2中紫色区域），在线状或面状地物上出现锯齿状或波浪状的边缘。这是因为MRA方法对于数据源配准误差非常敏感，两步融合均采用MRA类方法容易放大配准误差，进而造成严重的光谱失真或空间信息严重缺失。因此，这9组MRA-MRA配置不再参与后续讨论。

如图4（f）～（j）所示，第一步使用PCA法，第二步选用其他方法时，由于第一步的PCA法原理为保留第一主分量，即空间信息，舍去较多光谱信息，再进行第二步融合时，会造成光谱传递误差，也会导致出现严重的光谱失真，此类方法（表2中绿色区域）后续也不再评价和讨论。由图4（f）还可以看出：PCA-SFIM方法融合结果较差，不仅整体空间细节增强不够，光谱也出现了一定扭曲。原因在于：SFIM融合中的低分辨率影像合成涉及对全色影像较大范围的分块均值滤波操作，导致分块间的融合效果差异明显，进而导致整体的光谱失真；另一方面，PCA具有较好的空间增强能力也容易出现光谱失真，进一步将块与块之前的差异放大，产生了明显的块状效应和光谱扭曲，因此，PCA-SFIM方法也不参与后续讨论。

图4（k）展示了GSA组合方法结果较好的一组（ATWTGSA）。GSA方法要求待融合的全色影像和多光谱影像的波谱范围有所重叠，且满足线性关系（如式（2）所示）。而由于全色与短波红外影像波段范围不重叠，并不满足此线性假设，故容易造成严重的光谱失真，故将GSA方法用于两步式融合容易造成光谱失真。因此，与GSA相关的组合方法（表2中绿色区域）也不再参与后续讨论。

图4 部分两步式融合方法用于全色与短波红外影像融合的（较差）结果

综上所述，还剩余10种组合的融合结果如图5所示（对应表2中粉色区域）。

由图5可看出：对应的10种组合方法融合没有明显的光谱失真，较好地实现了细节信息的增强。左下角为对应的10种融合方法的局部细节，可以直观看出：HPF-PCA、SFIM-GS、HPF-GS三种组合方法，边缘细节信息丰富，光谱信息也保持较好。整体上看，两步式融合框架能有效地将全色影像的空间细节信息引入短波红外影像，并保留原始影像的光谱信息，有效地减少融合过程中存在的光谱失真和块状效应。

图5 两步式融合方法用于全色与短波红外影像融合的（较好）结果

表3为目视解译较优的2种单步融合方法和10种两步融合方法的5种评价指标的评价结果。由表可得：标准差最大的是SFIM-GS组合方法；平均梯度和空间频率最大的是SFIM-GS组合方法；GQNR指标最好的是HPF-GS。单一的CS和PCA方法，GQNR指标最差，为了能更好地进行组合方法选取，表4进行了各个信息的统计筛选，找到各个指标的前四名，记为较好、良好、一般。经过对比，最后综合选取的组合方法是HPF-GS法，即第一步融合使用HPF方法，第二步融合选用GS方法。

整体而言，本文设计的两步融合组合方式可分为两类。第一类为两步均使用同一类方法，即采用CS-CS或MRA-MRA类型的设置，具体如ATWT-HPF、SFIM-HPF、GS-GS等。这类组合方法融合效果普遍较差，这是因为CS或MRA方法，都有特定的优劣，如光谱保真的优劣或空间信息增强能力的强弱，两步均使用同一类方法，容易造成光谱或空间信息误差的累积传递，造成融合结果效果较差。第二类为两步分别使用不同CS-MRA或MRA-CS配置的融合方法，如HPF-GS、GS-HPF、GS-SFIM等，这类组合方法融合效果整体较好。这是因为两种方法中在光谱保真和空间效果增强具有较好互补性。实验结果还表明：HPF-GS法与GS-HPF两种组合方法中，效果更好的是HPF-GS，即先采用MRA方法，融合结果更优。

表3 （较优）单步融合和两步式融合方法融合结果定量评价

表4 综合评价结果

综上所述，针对Worldview 3影像的全色与短波红外波段融合，提出的两步式融合框架，采用合适的融合方法组合，能获得全色影像的空间细节信息，并保留原短波红外影像的光谱信息，达到增强空间分辨率并保持其光谱信息的目的。

3 结论与展望

本研究针对Worldview 3影像数据，选取了CS法和MRA法中的6种融合方法对其全色影像和短波红外影像进行两步式融合处理，针对两幅影像波谱范围不重叠且分辨率差别较大的问题，提出了两步式融合框架，并采用目视评价与定量评价方法对融合结果进行了评价。与单一的融合方法相比，本文提出的两步式框架目视效果与定量评价指标均取得良好结果，缓解了单一融合方法存在的块状效应、空间信息丢失和光谱失真的问题。在此框架下，两步分别采用属于MAR和CS类的方法，有助于获得更好的融合效果，综合最优的融合结果由HPS-GS方法获得。

在后期的研究中，将采用一些新的融合方法进行研究，如深度学习方法和基于局部自适应参数选择的方法等，来探寻更适用于全色与短波红外波段间的光谱关系的模型，从而提高短波红外影像的空间分辨率，并拓展到其他卫星传感器的全色与SWIR波段融合问题中。

致谢感谢DigitalGlobe公司提供本研究全部Worldview 3实验数据。感谢审稿人和编辑老师提出的宝贵意见，极大改善了论文的质量。