一种优秀的融合算法不仅可以极大地提高影像的空间解析特性，而且还可以最大限度的保持原始影像的光谱信息。然而，现有的高光谱影像融合算法在提高影像空间分辨率的同时，都会或多或少地对影像光谱信息造成一定的影响，造成地物像元光谱曲线失真，限制了影像后续的解译与应用范围[1-2]。为此，笔者从高光谱影像的光谱维入手，深入分析像元光谱曲线各个谐波分量的不同特性，提出了一种基于谐波分析的高光谱影像融合（Harmonic analysis fusion，HAF）算法[3]，该算法可以在提高影像空间分辨率的同时完全地保持光谱曲线的波形形态，但是由于直接采用高空间分辨率影像替换光谱曲线的谐波余项分量，导致原有的地物光谱曲线反射率受到了破坏。

Jakubauskas等人[4]为了对电力系统谐波污染进行有效地监测与控制，提出了谐波分析算法，该算法最初也主要被应用于该领域 [5-7]，在遥感领域的应用还比较少见。Sakamoto和Bradley等人[8-9] 在其研究成果中首次将谐波分析理论应用到AVHRR和MODIS等遥感影像植被的NDVI时间序列分析当中，以便了解农作物的季节变化情况和物候现象；Jakubauskas等人[4,10]在其后续的谐波研究中将该算法应用到了农作物物种识别当中；国内也有学者利用谐波分析理论对山东省典型植被季节变化特征进行了分析[11]；杨可明等人[12]则在谐波分析理论的辅助下对小目标地物进行了目标探测实验；文献[2-3]从高光谱影像的光谱维出发，深入地分析出光谱曲线谐波各个特征分量的物理意义，将谐波分析理论应用到了影像融合当中，并取得了一定的研究成果，但由于初次尝试，很多问题考虑不足，最终造成融合后地物光谱反射率整体出现过多偏移，地物光谱曲线物理解析特性遭到破坏，所以该算法还需进行改进与完善。

Gram-Schmidt（GS）变换融合算法可以较好地改善原始影像的空间细节特征，提高原始影像的空间分辨率，且能最大限度的保持原始影像的光谱物理特性，同时GS变换融合算法对于融合影像的波段数没有限制[13-14]。鉴于此，本文在原有的融合算法基础上，将GS变换引入到HAF算法当中，提出了结合Gram-Schmidt变换的高光谱影像谐波分析融合（Gram-Schmidt transform combined with Harmonic analysis fusion ，GSHAF）改进算法，以便对融合后的光谱曲线反射率进行有效地修正，使其更接近真实地物的光谱反射率，提高谐波分析融合算法的适用性。

1理论与方法

1.1 高光谱影像的谐波分析法

Rayne[15]和Davis[16]两位科学家认为，谐波分析是将目标的时间序列曲线f（t）表示为正（余）弦波相叠加的形式。谐波变换说明了任何连续的周期曲线f（t）都可以表示成为傅里叶级数形式[6]。通过学习高光谱遥感的基础知识可知，高光谱遥感影像具有极高的光谱分辨率，光谱分辨率只有几十甚至几纳米，影像中包含了成千上百个波段，这样多的波段使得像元在每个波段的灰度值可以表示为一条连续的曲线。既然每个像元在各个波段处的光谱值是一条连续的曲线，那么完全可以利用影像连续的光谱曲线v（s）来代替时间序列曲线f（t），这样便可用正（余）弦波相叠加的形式来表示光谱曲线v（s），同时也将像元光谱曲线分解为谐波余项、谐波振幅与谐波相位相叠加的形式，其变换公式如下所示

              （1）

其中，                                            （2）

式中，s代表波段号数；A0/2表示谐波余项；Z代表波段总数；h表示谐波分析次数；Chsin(2hπ/L+φh)表示第h次谐波分量；φh表示第h次谐波的初相位；Ch表示第h次谐波的振幅，Ch2表示功率谱。

通过观察式（1）和式（2）可知，高光谱遥感影像谐波分析法的处理单元为影像的单个像元，所以HAF算法属于像素级别的融合算法。 

1.2 高光谱影像的谐波分析融合算法

根据式（1）可知，谐波余项（A0/2）分量计算结果为一个常数，该值等于像元在每一个波段光谱值叠加之后的算数平均值。根据数学和物理知识可知，该分量在曲线叠加过程中并不影响光谱曲线波形的形状，只是对谐波振幅与谐波相位叠加之后的曲线进行一个简单的上下平移。谐波余项分量是地物对电磁波反射的综合反映，该分量中包含着影像丰富的空间信息，而光谱曲线的升降陡缓，波峰波谷出现的位置、大小全都由谐波振幅和谐波相位两个谐波分量决定。鉴于此，利用空间信息更为丰富的高空间分辨率影像直接替换谐波余项分量，如式（3）所示，随后对替换后的谐波特征分量进行谐波逆变换，完成高光谱影像与高空间分辨率影像的融合，这便是HAF算法的基本原理。

              （3）

式中，A代表高空间分辨率遥感影像，其他符号意义同式（1）。

1.3 Gram-Schmidt变换融合算法

     Gram-Schmidt（GS）[13-14, 17]变换是统计学中经常用到的一种多维线性正交变换，采用GS变换对高光谱影像多维数据进行正交化处理，可以有效地去除相邻波段间较强的相关性，最大程度地消除影像的信息冗余。Gram-Schmidt变换融合的基本流程如下：

（1）首先利用原始低空间分辨率影像模拟出一幅全色影像；

（2）随后利用该全色影像作为GS变换的第一个分量来对低空间分辨率影像进行GS变换，具体变换公式如下所示：

                 （4）

式中，GST表示经GS变换后的第T个正交分量，表示原始低空间分辨率遥感影像第T波段，为原始低空间分辨率遥感影像第T波段像元灰度值的均值，为原始低空间分辨率影像第K波段与GSl之间的协方差， i和j分别表示原始低空间分辨率影像的行数和列数，M和N表示整幅影像的行数和列数。

（3）用高空间分辨率影像替换GS变换后的第一分量，即GS1分量；

（4）最后再通过式（5）对上述替换后的数据集进行GS逆变换，完成低空间分辨率影像与高空间分辨率影像融合。

                （5）

式中，所有符号意义同式（4）。

1.4 结合Gram-Schmidt变换的谐波分析融合算法

HAF算法在保持融合影像光谱曲线波形形态方面极具优势，但是高光谱影像与高空间分辨率影像光谱范围不同，其所表示的光谱响应均值也不同，直接进行替换会对地物光谱曲线的反射率造成影响，限制了算法的适用性。通过前面的分析可知，原始高光谱影像光谱曲线的波形形态全都受控于谐波振幅与谐波相位，只要保证这两个谐波特征分量不受胁迫，则融合后影像的光谱曲线波形形态将不会受到任何影响，影像的光谱曲线波形可获得完整地保留，所以最后的问题全都归结到如何使谐波相位特征分量的空间分辨率得到提升。至此，我们在保证高光谱影像光谱曲线波形不变的基础上将高光谱影像的融合问题转换成了两幅单波段影像之间的融合，这样便将高光谱影像的融合操作大大简化。

针对之前直接使用高空间分辨率影像替换谐波余项造成HAF算法融合后影像在地物光谱反射率上出现整体失真这个问题，本文提出结合Gram-Schmidt变换的谐波分析融合算法。GS算法具有改善原始影像空间细节特征，提高原始影像空间分辨率，且能最大限度地保持原始影像的光谱物理特性等优点，最重要的一点是GS变换对于融合影像的波段数没有限制，完全满足本文的融合思路。

GSHAF算法具体流程是在原始高光谱影像被谐波分析算法分解为谐波余项、谐波振幅和谐波相位相叠加的形式后，根据式（4）和（5），首先对该光谱曲线的谐波余项分量与高空间分辨率影像进行GS变换融合，这样可以利用GS变换融合算法的优点既提高谐波余项的空间分辨率，同时也有效地修正了高空间分辨率影像与谐波相位像元灰度值之间巨大的差异，使GS变换融合后影像像元灰度值更接近谐波余项分量。随后对经GS变换融合后的影像与谐波振幅与谐波相位相两个特征分量根据式（6）进行谐波逆变换，完成高光谱分辨率与高空间分辨率的影像融合，这样便使融合后像元光谱曲线反射率更接近真实地物的光谱反射率，拓宽了融合影像后续应用范围，增强了GSHAF算法的适用性，GSHAF算法的具体流程如图1所示。

                 （6）

式中，G代表原始高光谱影像像元光谱曲线的谐波余项与高空间分辨率遥感影像经过GS变换融合后的影像，其他符号意义同式（1）。

图1 结合Gram-Schmidt的谐波分析融合算法流程

   Fig.1 GSHAF algorithm process

2 实验与分析

2.1 实验数据

实验数据为美国EO-1卫星在2005年7月底获取的Hyperion高光谱遥感影像与ALI高空间分辨率遥感影像全色波段影像。该影像区域为美国艾奥瓦州的一个小镇，如图2、3所显示，其中图2为

该地区经过假彩色合成后的Hyperion高光谱影像立体，图3中显示的是其对应的ALI高空间分辨率全色波段影像。实验区域植被覆盖率较高，主要植被为玉米（影像中黄色区域）。

                       图2 Hyperion高光谱影像

 图3 ALI 全色波段影像

 Hyperion高光谱影像共有242个波段，光谱波长范围为400～2500nm，地面空间分辨率为30m，光谱分辨率达到10nm。高级陆地成像仪ALI影像共有10个波段，其多光谱波段为30m的空间分辨率，ALI影像具有1个10m空间分辨率的全色波段，光谱范围为440～740nm，本次实验便采用ALI获取的全色波段充当影像融合所需的高空间分辨率影像。

2.2 谐波分析融合处理

实验数据经过了大气校正、垂直条纹去除、坏线修复、未定标波段剔除、配准、裁剪等一系列操作，详细预处理流程可参阅文献[18-20]。待原始影像完成预处理之后，首先对原始Hyperion高光谱影像按照式（1）进行最佳分解次数的谐波分解[3]，将其分解为谐波余项、谐波振幅与谐波相位相叠加的表示形式，随后根据式（4）、（5）首先采用GS变换融合算法对原始高光谱影像的谐波余项分量与ALI高空间分辨率影像进行融合，最后再按照式（6）对融合影像与原始高光谱影像的谐波振幅与谐波相位进行谐波逆变换，完成高光谱影像与高空间分辨率影像的融合。本实验所采用的算法均使用IDL（interactive data language）编程实现。融合结果如图4所示，为了直观地说明改进后的GSHAF算法比HAF算法具有更强的优越性，本次实验同样对实验数据进行HAF算法的融合处理，融合结果如图5所示，随后对HAF算法与GSHAF算法实验结果进行对比分析。 

图4  GSHAF算法融合后Hyperion影像

图5  HAF算法融合后Hyperion影像

           
2.3 GSHAF算法融合影像光谱信息保持性分析

    通过目视评价融合后影像可知， GSHAF算法与HAF算法融合后的影像在空间分辨率方面均具有很大的改善，影像分辨率明显得到提高，但是在影像色调方面GSHAF保持的更加贴近原始影像，HAF算法融合后的影像出现了色差，影像整体有点变色。这就直观地说明了GSHAF算法确实具有改进作用。为了更加客观地分析经过GSHAF算法与HAF算法融合后影像的光谱信息保持性情况，随机地在两种融合算法融合后影像上选择植被和建筑物两种地物的对应像元，并绘制每种地物像元的光谱曲线，如图6、7所示，图中虚线为HAF算法融合之后影像对应像元的光谱曲线，粗实线为GSHAF算法融合之后影像对应像元的光谱曲线，细实线为原始高光谱影像对应像元的光谱曲线。其中图6为植被融合前后的光谱曲线，图7表示的是建筑物融合前后的光谱曲线。从图中可以清晰的看出，GSHAF算法与HAF算法在光谱波形形态方面都完全地保留了原始影像的光谱曲线形态，但是GSHAF算法融合后的影像地物像元的光谱反射率与其原始的地物像元反射率更为接近，尤其在建筑物方面，几乎完全与原始光谱曲线重合。以上分析客观的证明了GSHAF算法在光谱信息保持方面更为优秀，更符合一个优秀融合算法的要求，其融合后的影像具有更为广泛的应用前景。

2.4 GSHAF算法的普适性验证

一个好的融合算法不仅需要具有优良的可行性，同时也需要具有极强的普适性，为了验证GSHAF算法是否具有极强的普适性，本节采用国产HJ-1A卫星获取的影像数据对本文提出的GSHAF算法进行普适性验证。HJ-1A卫星搭载了两种传感器，一种是可获取多光谱影像的电荷藕合器件（Charge Coupled Device，CCD）相机，另外一种是我国自主研制的空间调制型干涉高光谱成像光谱仪（HSI）。所以，HJ-1A星也具有同时获取高光谱遥感影像和多光谱遥感影像的能力，且其HSI光谱仪获取的高光谱遥感影像共有115个波段，光谱波长范围为459～956nm，光谱分辨率达到5nm，地面空间分辨率为100m；CCD传感器获取的多光谱遥感影像有4个波段，地面空间分辨率为30m。本次普适性验证实验数据采用2013年4月25日获取的陕西某城市郊区及其附近山区的影像，如图8所示为经过假彩色合成后的该区域高光谱影像立体，图中红色带状地物为水体，蓝色区域为山上植被；经观察发现，CCD获取的多光谱遥感影像第三波段的空间解析特性相对较为突出，该波段影像具有较为丰富的空间信息，所以采用多光谱遥感影像的第三波段充当高空间分辨率影像，如图9所示。

同样对HJ-1A卫星的数据分别进行GSHAF与HAF两种融合算法的融合，其结果如图10、11所示。通过目视评价可知，经过GSHAF算法与HAF算法融合后的影像空间分辨率都获得了很大的提高，在其色调方面两幅融合后影像整体差异不大，但是仔细观察融合前后影像中的河流的话，还是可以发现GSHAF算法融合后的影像色调更接近原始影像。

为了客观了解其在光谱信息保持性方面的情况，同样绘制出融合前后影像中对应位置随机像元A、B的光谱曲线，如图12、13所示，图中三条曲线所代表的影像像元意义与图6、7相同。可以明显的看出，其所表示出来的结果与可行性分析实验完全相同，光谱曲线波形形态两种方法都完全地保留下来，反射率方面GSHAF算法融合影像更为接近原始影像，图13所示GSHAF算法融合后的光谱曲线也几乎与原始影像光谱曲线重合，以上分析客观真实地证明了GSHAF算法不仅具有优良的可行性，还具有极强的普适性。    

3 结论

本文在分析了基于谐波分析理论的高光谱遥感影像融合算法缺陷的基础上，提出了结合Gram-Schmidt变换的谐波分析融合改进算法，该算法可以在完全保留融合前后影像光谱曲线波形形态的基础上，将高光谱影像融合简化为两幅单波段影像之间的融合，相比HAF算法，GSHAF算法可以有效地减少HAF算法融合后像元光谱曲线与原始光谱曲线反射率差异较大、地物像元光谱物理特性失真这一缺陷，使融合后的高空间分辨率与高光谱分辨率影像获得更广泛的用途，并通过实验来证明这一假设。通过实验与分析可得：1）HAF算法和GSHAF算法在影像融合过程中都可以完全地保留原始影像的光谱曲线波形形态；2）GSHAF算法在光谱物理特性保持方面比HAF算法更具优势，融合后地物的光谱反射率值更接近原始影像地物光谱反射率值；3）GSHAF算法不仅具有优良的可行性，还具有极强的普适性；4）虽然GS变换对HAF算法具有一定的改进作用，但是其还是无法完全地保留地物光谱全部的物理特性，然而，本文已经在保证融合前后光谱曲线波形不变的情况下，将高光谱遥感影像融合简化为两幅单波段影像之间的融合，只要找到合适的算法，相信完全可以最大限度地将地物光谱信息保留下来。

参考文献

[1] 董广军，张永生，范永弘. PHI高光谱数据和高空间分辨 率遥感图像融合技术研究[J]. 红外与毫米波学报，2006， 25(2)：123-126.

[2] 杨可明，张涛，王立博，等. 谐波分析法高光谱影像融合及其光谱信息保真度评价[J]. 光谱学与光谱分析，2013，33(9)：2498-2501.

[3]杨可明，张涛，王立博，等. 高光谱影像的谐波分析融合算法研究[J]. 中国矿业大学学报，2014，43(3)：547-553.

[4]JAKUBAUSKAS M E，LEGATES D R，KASTENS J H. Harmonic Analysis of Time-Series AVHRR NDVI Data [J]. Photogrammetric Engineering & Remote Sensing，2001，67(4)：461-470.

[5]任子晖，仇润鹤，张艳. 煤矿电网谐波分析模型的建立与滤波器设计[J]. 中国矿业大学学报，2004，33(1)：45-49.

[6]GEORGE J，WAKILEH G J. 电力系统谐波——基本原理、分析方法和滤波器设计[M]. 徐 政，译. 北京: 机械工业出版社，2003：4-13.

[7] 柴旭铮，文习山，关根志，等. 一种高精度的电力系统谐波分析算法[J]. 中国电机工程学报，2003，23(9)：68-71.

[8] SAKAMOTO T，YOKOZAWA M，TORITANI H，et al. A crop phenology detection method using time-series MODIS data [J]. Remote Sensing of Environment，2005，96：366-374.

[9]BRADLEY B A，JACOB R W，HERMANCE J F，et al. A curve fitting procedure to derive inter-annual phenologies from time series of noisy satellite NDVI data [J]. Remote Sensing of Environment，2007，106：137-145.

[10] JAKUBAUSKAS M E，LEGATES D R，KASTENS J H. Crop identification using harmonic analysis of time-series AVHRR NDVI data [J]. Computers and electronics in agriculture，2002，37：127-139.

[11] 梁守真，邢前国，施平，等. 山东省典型地表覆被NDVI时间序列谐波分析[J]. 生态学杂志，2011，30(1)：59-65.

[12] 杨可明，薛朝辉，贾涛涛，等. 高光谱影像小目标谐波分析探测模型[J]. 测绘学报，2013，42（1）：34-43. 

[13]李存军，刘良云，王纪华，等. 两种高保真遥感影像融合方法比较[J]. 中国图象图形学报，2004，9(11)：1376~1385.

[14] 于海洋，闫柏琨，甘甫平，等. 基于Gram-Schmidt变换的高光谱遥感图像改进融合算法[J]. 地理与地理信息科学，2007，23(5)：30-42

[15] RAYNER J N. An Introduction to Spectral Analysis [M]. London：Pion Ltd，1971：15-21

[16] DAVIS J C. Statistics and Data Analysis in Geology (2nd edition) [M]. New York：John Wiley & Sons，1986：23-34.

[17] 黄登山. 像素级遥感影像融合方法研究[D]. 长沙：中南大学，2011：18.

[18] 谭炳香，李增元，陈尔学，等. EO-1 Hyperion高光谱数据的预处理[J]. 遥感信息，2005，(6)：36-41.

[19] 陈建珍，何超，岳彩荣. 基于FLAASH模块的高级陆地成像仪图像的大气校正[J]. 浙江农林大学学报，2011，28(4)：590-596.

[20] CHANDER C，MARKHAM B L，HELDER D L. Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors [J]. Remote Sensing of Environment，2009，113：893-903.