遥感图像融合是图像融合的一个重要分支，针对全色光图像与多光谱图像的融合始终是遥感图像处理领域的研究热点。随着传感器技术的飞速发展，许多在轨卫星，如SPOT、IKONOS、QuickBird等能够同时提供高分辨率全色光图像与多光谱图像。通常，多光谱图像包含多个波段，具有较高的光谱分辨率，但其空间分辨率相比于单波段全色光图像较低。图像融合技术就是利用高分辨率全色光图像去提高多光谱图像的空间分辨率。通过融合得到的具有高空间分辨率的多光谱图像可用于诸如谷歌地球、城市规划、高精度地物分类等众多领域。

新型WorldView-2(WV-2)星载图像代表了下一代超高分辨率遥感图像的发展趋势，同时也对现有的图像融合方法提出了挑战^[1-2]。为了突破输入波段数量的限制，文献[3]提出了一种基于超球面彩色变换(hyperspherical color sharpening, HCS)的融合方法，将WV-2所有8波段多光谱数据一次性地从原始空间转换到超球面彩色空间进行处理。尽管该方法已成为ERDAS IMAGINE 2011商业软件的模块，但是其融合结果仍会出现明显的光谱失真。文献[4]提出一种带有可调参数(adjustable pan-sharpening, APS)的图像融合方法，以处理包含WV-2在内的多种高分辨率遥感图像。APS方法通过分析不同星载传感器强度分量匹配关系，利用两个可调参数整合了传统的色度-强度-保和度(intensity-hue-saturation, HIS)、Bravery变换(Bravery transform, BT)以及光滑调制滤波(smoothing-filter-based intensity modulation, SFIM)融合方法，实现了在5种融合模式中的切换。文献[5]利用WV-2多光谱数据地物分类的先验信息设计了一种基于广义HIS(generalized HIS, GIHS)变换的融合方法，然而地物分类精度对融合质量的影响并未展开讨论。

本文针对WV-2传感器各波段的光谱响应特点，提出了一种基于对应分析的融合新方法。利用多光谱和全色光波段间的相对光谱响应关系构造出合适的低分辨率全色光图像，进而提取出原始全色光图像的空间细节信息，并注入至对应分析变换的成分空间中，最终得到高质量的融合结果。

1 WorldView-2遥感图像

2009年10月8日发射升空的WV-2卫星能够同时提供8波段多光谱图像(Band₁=Coastal、Band₂= Blue、Band₃=Green、Band₄=Yellow、Band₅=Red、Band₆=RedEdge、Band₇=NIR₁、Band₈=NIR₂)和单波段全色光图像(panchromatic，PAN)。多光谱传感器的光谱覆盖范围为400~1 050 nm，空间分辨率为1.84 m；全色光传感器光谱范围为450~800 nm，空间分辨率为0.46 m。其相对光谱响应如图 1所示。

与QuickBird、IKONOS、GeoEye-1等高分辨率遥感图像相比较：1) WV-2数据包含的多光谱波段划分更细、通道数更多；2)全色光波段的光谱覆盖范围变窄；3)多光谱与全色光波段光谱响应之间的匹配发生了明显的变化，特别是Band₁和Band₈基本上在PAN波段光谱响应之外。正是这些新特点给遥感图像融合研究带来了新的挑战，现有的许多融合方法并不适合于WV-2遥感图像，若直接套用，则会出现严重的光谱和空间信息失真现象。

2 对应分析变换

对应分析(correspondence analysis, CA)是近年来发展起来的一种多元相依变量统计分析技术，通过分析由定性变量构成的数据列联表来揭示同一变量的各个类别之间的差异，以及不同变量各个类别之间的对应关系。若数据列联表中的元素是同源的且非负，那么也可以用对应分析方法对其分析^[6]。针对多光谱遥感图像，可以利用${\chi ^2}$距离计算波段之间的关联性，即变换之后的数据列联表具有Pearson(${\chi ^2}$)统计特性。

首先，将列联表中每个元素${x_{ij}}$除以数据集中所有元素之和${x_{ + + }}$得到对应的比例系数${p_{ij}}$，由${p_{ij}}$组成新的数据列联表$P$。行权系数${p_{i + }}$等于${{{x_{i + }}} \mathord{\left/ {\vphantom {{{x_{i + }}} {{x_{ + + }}}}} \right. } {{x_{ + + }}}}$，其中${x_{i + }}$表示第$i$行所有元素之和。列权系数${p_{ + j}}$等于${{{x_{ + j}}} \mathord{\left/ {\vphantom {{{x_{ + j}}} {{x_{ + + }}}}} \right. } {{x_{ + + }}}}$，其中${x_{ + j}}$表示第$j$列所有元素之和，列联表中每个元素的值为：

${x_{ij}} = \sqrt {{x_{ + + }}} \frac{{{p_{ij}} - {p_{i + }}{p_{ + j}}}}{{\sqrt {{p_{i + }}{p_{ + j}}} }}$

(1)

若使用${q_{ij}} = {{{x_{ij}}} \mathord{\left/ {\vphantom {{{x_{ij}}} {\sqrt {{x_{ + + }}} }}} \right. } {\sqrt {{x_{ + + }}} }}$替代${x_{ij}}$，则由${q_{ij}}$组成的矩阵可以表示为：

$\mathit{\boldsymbol{Q}} = [{q_{ij}}] = \left[ {\frac{{{p_{ij}} - {p_{i + }}{p_{ + j}}}}{{\sqrt {{p_{i + }}{p_{ + j}}} }}} \right]$

(2)

其特征向量矩阵为：

$\mathit{\boldsymbol{U}} = {\mathit{\boldsymbol{Q}}^{\bf{T}}}\mathit{\boldsymbol{Q}}$

(3)

这样，利用特征向量矩阵${\mathit{\boldsymbol{U}}}$可以将多光谱数据${\mathit{\boldsymbol{X}}}$转换到成分空间${\mathit{\boldsymbol{Y}}}$，即：

${\mathit{\boldsymbol{Y}}} = {{\mathit{\boldsymbol{U}}}^\bf{T}}{\mathit{\boldsymbol{X}}}$

(4)

其反变换为：

${\mathit{\boldsymbol{X}}} = {({{\mathit{\boldsymbol{U}}}^\bf{T}})^{ - 1}}{\mathit{\boldsymbol{Y}}}$

(5)

3 融合方法

一个好的融合方法在提高多光谱图像空间分辨率的同时，还应该尽可能地保持多光谱图像的原始光谱信息。本文针对WV-2遥感图像提出了一种基于对应分析的融合方法，将输入波段数扩展至任意多个，其融合框架如图 2所示。

由图 2可知，该融合方法分以下7步实施。

1) 将8波段多光谱图像(Band₁, Band₂, $ \cdots $, Band₈)重采样放大至与全色光图像尺寸一致。采用最近邻插值法得到放大的低分辨率多光谱图像，记为LRM₁, LRM₂, $ \cdots $, LRM₈。最近邻插值法使用原始像素点复制，不会引入新的像素值，因此采样后的多光谱图像在光谱保持能力方面优于双线性插值或双三次插值。

2) 利用多元线性回归方法预测多光谱图像与全色光图像间的关系。理论上可以假设全色光图像PAN与高分辨率多光谱图像HRM满足线性关系，即${\rm PAN} = \sum\limits_{i = 1}^8 {{\omega _i}{\rm HR}{{\rm M}_i}} + b$。因此，低分辨率的全色光图像LRP与低分辨率的多光谱图像也应满足线性关系，即${\rm LRP} = \sum\limits_{i = 1}^8 {{\omega _i}{\rm LR}{{\rm M}_i}} + b$。

考虑到WV-2数据中全色光图像与多光谱图像的Band 1和Band 8的光谱响应几乎没有重叠(如图 1所示)，因此本文只选择LRM₂, LRM₃, $ \cdots $, LRM₇共6幅多光谱图像参与拟合，即：

${\rm PAN} = \sum\limits_{i = 2}^7 {{\omega _i}{\rm LR}{{\rm M}_i}} + b$

(6)

3) 全色光图像的空间细节提取。利用式(6)计算得到的权系数ω_i及偏置b可进一步构造出低分辨率的全色光图像，即：

${\rm LRP} = \sum\limits_{i = 2}^7 {{\omega _i}{\rm LR}{{\rm M}_i}} + b$

(7)

构造出的LRP与LRM在空间分辨率上是相同的，可以通过原始全色光图像与构造出的低分辨率全色光图像的比值来获取空间细节信息，即：

$\delta = {\rm PAN}/{\rm LRP}$

(8)

4) 光谱图像的对应分析变换。利用重采样的多光谱数据构造输入矩阵：

${\mathit{\boldsymbol{X}}_{{{\rm IN}}}} = [{\rm LR}{{\rm M}_{\rm 1}},{\rm LR}{{\rm M}_{\rm 2}}, \cdots ,{\rm LR}{{\rm M}_{\rm 8}}]$

(9)

按照对应分析变换式(1)~式(4)将输入矩阵X_IN变换到成分空间，得到成分空间[C₁, C₂, $ \cdots $, C₈]。

5) 空间细节的注入。由于空间细节图像$\delta $的均值接近1，因此采用乘性注入方式所产生的光谱畸变较小。将细节信息注入至第八成分，即：

${C_8}_{{{\rm NEW}}} = \delta {C_{\rm 8}}$

(10)

6) CA反变换。用C_8NEW替换C₈，再按照式(5)对新成分空间[C₁, C₂, $ \cdots $, C_8NEW]进行对应分析反变换，即可得到高分辨率的多光谱图像(HRM₁, HRM₂, $ \cdots $, HRM₈)。

7) 小波变换。重采样多光谱图像由于采用最近邻插值会产生方块或锯齿效应，特别是在轮廓和边缘部分较为明显，因此可以引入冗余小波变换方法去除步骤6)得到的高分辨率多光谱图像(HRM₁, HRM₂, $ \cdots $, HRM₈)中的锯齿效应。本文采用à trous小波算法实现^[7]：首先对PAN和HRM_i分别进行小波分解，得到各自的近似分量和高频细节；然后用PAN的高频细节去替换HRM_i的高频细节，利用HRM_i的近似分量和替换的高频细节进行重构，得到最终的融合结果SB i(i = 1, 2, $ \cdots $, 8)。

4 实验及结果分析 4.1 实验数据及评价指标

WV-2实验数据集包含空间分辨率为2.0 m的8波段多光谱图像和0.5 m的全色光图像，原始尺寸1600×1600像素，场景是2011年4月拍摄的澳大利亚悉尼港区域，其RGB彩色合成图像(Band 5-3-2)如图 3所示。多光谱图像与全色光图像已经过配准。为了验证本文方法的有效性和先进性，先后与文献[4]的APS方法、文献[3]的HCS方法以及文献[8]的IAWP方法的实验结果进行对比，并采用相关系数(correlation coefficient, CC)、空间相关系数(spatial CC, SCC)、相对全局维数综合误差(relative global dimensional synthesis error, ERGAS)以及光谱角映射(spectral angle mapper, SAM)4种常用客观评价指标对各方法的融合结果进行客观评价^[9-12]。为方便显示，图 4给出了400×400像素的子场景。融合实验的硬件环境为Intel CPU 3.10 GHz，4 GB内存，64 bit Windows 7操作系统；编程环境为ENVI 4.1/IDL 6.1。

4.2 结果分析

由图 4可知，该子场景包含地物内容较丰富，有植被、建筑物、车辆、道路等。4种方法的融合结果与原始多光谱图像(图 4a)相比，其空间分辨率均有不同程度地提高。本文方法与APS方法在空间分辨率提升方面效果显著，图像的清晰程度基本能与全色光图像(图 4b)保持一致。但APS方法在光谱保持方面却略显不足，特别是在绿色植被区域将原有的深绿色扭曲为浅黄绿色，出现了色彩失真。

HCS方法的空间分辨率提高不够，并且在高亮度的白色屋顶部分出现了明显的光谱失真。IAWP方法尽管其光谱保持能力较强，但细节信息注入不够，导致图像的空间信息失真严重，视觉效果较差。

4种融合方法的客观评价结果如表 1所示。从表中可以看出，IAWP方法的CC数值说明其光谱保持能力强，但是SCC评价最差说明其空间分辨率增强效果最差。APS方法的SCC评价最高，但是其ERGAS、SAM和CC数值说明该方法的光谱保持能力以及总体质量不是最优。本文方法在ERGAS和SAM数值优势明显，SCC评价非常接近APS方法，说明本文方法在空间分辨率提高和光谱信息保持方面达到了更好的平衡，融合结果的整体质量最好，主观评价与客观分析结果能够达到一致。

5 结论

针对新型高分辨率WorldView-2遥感图像，本文从传感器各波段光谱响应特点入手，有选择性地利用多光谱图像构造低分辨率全色光图像，进一步提取出原始全色光图像的空间信息并将其注入至多光谱图像的成分空间中，经对应分析反变换，再利用冗余小波变换去除各波段的块效应得到最终融合结果。融合实验对几种方法进行验证，实验结果的主观评价和客观分析表明，本文提出的新方法相比于其他几种融合方法在提高空间分辨率和保持光谱信息方面能够达到更好的平衡。