0 引言

随着遥感技术不断发展，对同一地区拍摄的遥感图像数量越来越多。由于成像原理和成像条件存在差异，不同来源的遥感图像各具优势也各有局限。图像融合技术则可以最大限度地提取各个来源图像的有利信息，进而综合形成高质量图像^[1]。采用遥感图像融合技术，将遥感图像的全色波段图像和多光谱图像融合，可以充分利用图像互补性，有效提高多光谱图像的空间分辨率，为进一步实际应用提供高质量的彩色遥感影像。许多经典的遥感图像融合算法，如IHS(Intensity，Hue，Saturation)变换法^[2]、主成分分析(Principal Component Analysis,PCA)法^[3]、小波变换(Wavelet Transform，WT)法^[4-5]和Brovey算法^[6-9]等，在遥感图像融合领域都有广泛应用。但是，通过对上述算法原理分析和实验研究发现，算法在综合考虑空间信息增强和光谱信息保持方面，通常只能满足某一特定方面的应用需求^[10]。

为此，研究学者提出了诸多改进算法。其中，针对IHS融合算法的改进已有许多成果发表。伍娟等^[11]利用多光谱影像I分量(亮度分量)的直方图为参照匹配全色影像，通过建立I分量和全色图像(PANchromatic image，PAN)的灰度分布概率密度函数之间的映射关系，完成PAN对I分量的替代，获得了光谱失真较少、清晰度较高的融合影像。肖化超等^[12]对多光谱图像进行IHS变换，然后将代表空间分辨率的I分量与全色图像进行Curvelet变换，并进行标准差融合，获得了优于传统方法的融合图像； Huang等^[13]对IHS空间中I分量的计算进行重新定义，提出了一种快速的IHS融合方法，实现了多光谱与全色影像的有效融合； 徐佳等^[14]通过线性回归的方法重新构造了亮度分量，并利用空间细节注入方法，调整与原始多光谱图像的空间细节信息比例，实现了融合图像的质量提高； 高永刚等^[15]将全色图像与多光谱图像的三个波段分别作直方图匹配，并对匹配图像求均值后替代I分量参与反变换，获得了光谱失真度较小的高质量融合图像。以上改进算法能够较好地保持融合图像的空间信息，也同时在一定程度上减少了融合图像光谱失真和扭曲的发生。但综合两者在融合图像中的表现，仍无法同时达到最优。

因此，如何兼顾融合图像空间信息保持度和光谱信息保真度，使两者在融合图像中表现最优，是遥感图像融合技术中有待解决的难点问题。针对这一问题，本文提出了一种可以保持融合图像空间分辨率信息，又能兼顾抑制融合图像光谱失真的算法，并将其集成在经典的IHS算法框架中。实验结果表明，采用本文算法获得的融合图像光谱失真度小、空间信息保持度高。

1 IHS融合图像光谱失真分析

根据图像的色彩变换理论可知，某种颜色既可以用RGB色彩空间内的R、G、B(红、绿、蓝)来描述，也可以用IHS色度空间内的I、H、S(亮度、色调、饱和度)描述。IHS变换就是在RGB空间与IHS空间之间进行变换^[16]。IHS融合在IHS变换的基础上实现，算法原理如图 1所示。

首先，将多光谱图像的三个波段由彩色空间变换到IHS空间；然后，用高分辨率的全色图像替代IHS空间的亮度分量；最后，经IHS反变换得到融合图像。具体计算步骤如下：

1) 根据式(1) 将需要融合的多光谱图像由彩色空间变换到IHS空间。

$\left[ \begin{array}{l} I\\ {v_1}\\ {v_2} \end{array} \right] = \left[ {\begin{array}{*{20}{c}} {1/\sqrt 3 }&{1/\sqrt 3 }&{1/\sqrt 3 }\\ {1/\sqrt 6 }&{1/\sqrt 6 }&{ - 2 \times \sqrt 6 }\\ {1/\sqrt 2 }&{ - 1/\sqrt 2 }&0 \end{array}} \right]\left[ \begin{array}{l} \mathit{\boldsymbol{R}}\\ \mathit{\boldsymbol{G}}\\ \mathit{\boldsymbol{B}} \end{array} \right]$

(1)

2) I分量被PAN替代；

3) 根据式(2) 完成IHS反变换。

$\begin{array}{l} F(\mathit{\boldsymbol{X}}) = \left[ \begin{array}{l} F(\mathit{\boldsymbol{R}})\\ F(\mathit{\boldsymbol{G}})\\ F(\mathit{\boldsymbol{B}}) \end{array} \right] = \left[ {\begin{array}{*{20}{c}} {1/\sqrt 3 }&{1/\sqrt 6 }&{1/\sqrt 2 }\\ {1/\sqrt 3 }&{1/\sqrt 6 }&{ - 1/\sqrt 2 }\\ {1/\sqrt 3 }&{ - 2 \times \sqrt 6 }&0 \end{array}} \right]\left[ \begin{array}{l} \mathit{\boldsymbol{PAN}}\\ \;\;{\mathit{\boldsymbol{v}}_1}\\ \;\;{\mathit{\boldsymbol{v}}_2} \end{array} \right] = \\ \left[ {\begin{array}{*{20}{c}} {1/\sqrt 3 }&{1/\sqrt 6 }&{1/\sqrt 2 }\\ {1/\sqrt 3 }&{1/\sqrt 6 }&{ - 1/\sqrt 2 }\\ {1/\sqrt 3 }&{ - 2 \times \sqrt 6 }&0 \end{array}} \right]\left[ \begin{array}{l} \mathit{\boldsymbol{I}} + (\mathit{\boldsymbol{PAN}} - \mathit{\boldsymbol{I}})\\ \;\quad {\kern 1pt} \quad \;{\mathit{\boldsymbol{v}}_1}\\ \;\;\quad \quad {\mathit{\boldsymbol{v}}_2} \end{array} \right] = \left[ \begin{array}{l} \mathit{\boldsymbol{R}} + (\mathit{\boldsymbol{PAN}} - \mathit{\boldsymbol{I}})\\ \mathit{\boldsymbol{G}} + (\mathit{\boldsymbol{PAN}} - \mathit{\boldsymbol{I}})\\ \mathit{\boldsymbol{B}} + (\mathit{\boldsymbol{PAN}} - \mathit{\boldsymbol{I}}) \end{array} \right] \end{array}$

(2)

其中:F(X)表示融合图像，X为R、G、B三个分量的组合。从式(2) 可以看出，当全色图像PAN与多光谱图像I之间信息差异较大时，将导致融合图像产生光谱信息损失^[17]。实际上，全色图像的空间分辨率高，光谱分辨率较低；而多光谱图像的光谱分辨率高，空间分辨率较低。全色图像PAN与多光谱图像I携带的信息存在固有差别，光谱损失在基于IHS方法的融合图像中不可避免。因此，遥感图像融合过程要综合考虑提高空间分辨率与降低光谱失真程度，在满足实际应用的前提下，两者尽可能同时达到最优。

2 本文算法

本文提出了一种直方图中轴化融合算法，通过折中策略调整全色图像和多光谱图像亮度分量两者的灰度直方图，降低了两者的信息差异，最终达到提高融合图像空间分辨率和减少光谱失真之间的有效平衡。

2.1 算法原理

在IHS融合算法的框架下，假设用I_new表示融合后图像的亮度I分量，I_old表示融合前图像的亮度分量，那么融合前后的亮度差异图像即可表示为|I_new-I_old|，该值的大小表征了融合前后图像的变化程度。有以下两种情况：①当用全色图像完全替代多光谱亮度分量，即I_new=PAN时，|I_new-I_old|值最大，融合结果的空间分辨率最高，光谱扭曲最严重；②当多光谱图像亮度分量等于其本身，即I_new=I_old时，|I_new-I_old|值最小为0，融合结果光谱扭曲程度最小，但空间分辨率最低。因此，同时兼顾融合图像分辨率提高程度最高和光谱失真最小条件时，IHS算法可以表示为式(3) 的优化问题：

$\mathop {\min }\limits_{{\mathit{\boldsymbol{I}}_{new}}} \left\{ {{{\left| {\mathit{\boldsymbol{PAN}} - {\mathit{\boldsymbol{I}}_{new}}} \right|}^2} + {{\left| {{\mathit{\boldsymbol{I}}_{new}} - {\mathit{\boldsymbol{I}}_{old}}} \right|}^2}} \right\}$

(3)

其中：当I_new与PAN图像接近时，可以达到空间信息保持度最大；当I_new与I_old图像相近时，可以保证光谱失真程度最小。那么，当I_new是PAN和I_old图像的折中时，式(3) 的值达到最小。

基于以上分析，计算PAN和I_old图像的折中I_new，需要选择一个合理的映射规则H。本文发现，如果H是位于PAN和I_old图像累积直方图的中轴线，那么将两幅图像的累积直方图都向这一中间曲线作变换，可以得到PAN和I_old图像最为理想的折中图像。

下面用两幅人工图像说明这一过程。首先，假设PAN和I_old图像为两幅已知图像，表示为I_i(i=1,2) ，它们的累积直方图为H_i(i=1,2) 。求取两幅图像的折中图像，需要选取一个最优映射规则，涉及的累积直方图及其变换的定义和性质^[18]如下。

定义1 令I_i:Ω→[0,1]¹ (Ω⊆R²,i=1,2) 是一对图像，h_i是图像I_i的直方图。那么，图像I_i的累积直方图函数H_i:[0,1]→[0,1]。可以定义如式(4) ：

${H_i}(x) = \int_0^x {{h_i}(t){\rm{d}}t} $

(4)

性质1 设H_i:[0,1]→[0,1]是图像的累积直方图，是一个连续单调递增函数，则有变换H_i′=H^-1[H_i]，变换后的累积直方图为H_i′。

如图 2(a)~(f)是两幅人工图像以及变换图像的比较。其中，图 2(a)、(b)是两幅原始图像，图 2(c)和(d)是原始图像的灰度直方图，图 2(e)是累积直方图，其中H₁和H₂分别是原始图像的累积直方图，H_Aver是H₁和H₂的算数平均轴，由式(5) 求得；H_Midway是H₁和H₂的几何中间轴，由式(6) 求得。

${\mathit{\boldsymbol{H}}_{Aver}} = \frac{{{\mathit{\boldsymbol{H}}_1} + {\mathit{\boldsymbol{H}}_2}}}{2}$

(5)

${\mathit{\boldsymbol{H}}_{Midway}} = {\left( {\frac{{{\mathit{\boldsymbol{H}}_1}^{ - 1} + {\mathit{\boldsymbol{H}}_2}^{ - 1}}}{2}} \right)^{ - 1}}$

(6)

图 2(f)是根据H_Aver(算数平均轴)对原始图像的变换，变换后图像产生了原来没有的结构；图 2(g)是根据H_Midway(几何平均轴)对原始图像的变换，变换后图像信息介于两幅原始图像之间，且没有产生额外的结构信息。可以看出，H_Midway是两幅原始图像折中图像最佳变换规则。

本文把根据H_Midway进行的图像直方图变换的方法称为“直方图中轴化”。在IHS算法框架下，对原始多光谱图像亮度分量与全色图像替换之前，先将两者作直方图中轴化变换，将变换后的图像作为新的亮度分量，再参与IHS反变换得到融合图像。由于新的亮度分量图像同时实现了将全色图像携带的空间信息和多光谱图像亮度分量携带的光谱信息的均衡化，所以本文算法可以综合提高融合图像的空间信息和光谱信息的质量。

2.2 算法流程

本文算法流程如图 3所示，执行过程如下：

① 首先对多光谱图像进行重采样，使之与全色图像PAN尺寸一致，之后对多光谱图像进行IHS变换，得到各分量表示为I_old、H、S；

② 根据式(6) 计算PAN和I_old的几何中间轴，再根据几何中间轴变换PAN和I_old，获得理想折中图像I_new；

③ 用折中图像I_new替换I_old分量；

④ 进行IHS反变换，得到融合图像。

由于I_new与全色图像PAN和多光谱图像的亮度分量I_old接近程度均匀一致，因而可以最大化地保持PAN图像的结构信息，并使图像光谱信息损失最小。

3 实验结果与分析

为验证本文算法的正确性和有效性，将本文算法与四种传统图像融合算法(IHS、PCA、WT、Brovey)进行主观和客观两方面的比较。实验结果如图 4所示，相应的客观评价指标如表 1所示。本文实验数据来自中国资源卫星应用中心(http://www.cresda.com/CN/Satellite)提供的卫星图像以及标准测试数据。为了方便进行处理，从全色图像中选取500像素×500像素的区域，在多光谱图像中选取同样区域进行融合。所有实验在具有8 GB内存的PC上进行，使用了Matlab R2012a编程语言。

3.1 主观评价

从图 4可以看出，IHS、PCA、WT、Brovey以及本文算法对图像的融合能力各不相同。以图 4(a)为例分析说明：IHS算法融合图像的视觉效果虽然较亮，但是地物色彩有较大畸变，这说明图像的光谱信息在融合过程中产生了损失；PCA方法融合图像的光谱色调较为符合实际，但是当待融合图像本身比较暗时，图像的光谱保真度并不理想；WT方法融合图像在光谱失真方面有一定程度的减少，并且色彩也更加接近融合前低分辨图像，但是融合图像的清晰度略有下降，清晰度不高的融合图像并不能满足后续使用的要求；Brovey算法融合图像的光谱失真程度比较低，但是其视觉亮度较暗，仍不利于后续使用；本文算法融合图像的色彩与原始低分辨图像最为相似，这说明本文算法对光谱损失实现了有效的抑制，同时图像亮度较亮、纹理清晰。可以看出本文算法可以有效抑制光谱失真，同时对图像的亮度、纹理等方面的影响较小。

3.2 客观评价

从表 1的客观评价指标数据也能够得出与上述视觉效果相同的结论。从各个实验结果客观评价指标以及均值来看，本文算法在“光谱扭曲度”指标中的数值最低，说明通过调整替换前后图像有效减少了融合图像的光谱损失，提高了融合图像光谱信息量，符合本文理论分析结果；“峰值信噪比”指标中本文算法的均值最高，说明本文算法融合的图像信息损失最低； “标准偏差”“平均梯度”“相关系数”“信息熵”四个指标中本文算法的均值同样高于其他算法，说明融合图像的清晰度较好，融合图像与多光谱图像相似程度最高、信息量最大，有利于视觉分析；另外，本文算法在“灰度均值”这项指标中，虽没有达到最高，但是与其他高灰度均值的指标值相当，从视觉效果也可以看到本文算法的融合图像亮度完全可以满足视觉对图像明暗程度的要求。综合各项指标来看，本文算法生成的融合图像光谱失真最低，融合图像包含的信息最丰富，同时图像清晰度和图像亮度也能够满足应用要求。 从表 1的算法计算耗时比较可以看出：在融合同一组图像时，本文算法与其他算法相比用时略少。但是，本文算法融合真实卫星图像的时间过长，是融合标准图像所用时间的近15倍，原因是算法在实施过程可能存在较大的数据冗余，算法还可以进一步优化提速。

综合实验结果可知，本文方法生成的融合图像在图像光谱保留、空间细节表现能力等方面表现优良。本文基于直方图中轴化的图像融合算法可以有效降低融合图像的光谱扭曲度，并保证融合图像具有较高的清晰度。

4 结语

在全色图像和多光谱图像融合过程中，直接利用IHS变换进行融合时，由于存在光谱特征扭曲，使融合后图像的部分光谱信息产生失真，对图像的识别和解译造成了不利影响。本文采取直方图中轴化策略，改进了以往的IHS变换中直接舍弃多光谱图像I分量的处理，最大程度地保留了原多光谱影像的光谱信息，同时也提高了融合图像空间分辨率。但本文算法的计算速度还有待优化和提高。