張 濤，劉 軍，楊可明，羅文杉，張育育

1.國家測繪地理信息局第一航測遙感院，陜西 西安710054；2.中國礦業(yè)大學（北京）地球科學與測繪工程學院，北京100083

現(xiàn)有的高光譜影像融合算法在提高影像空間分辨率的同時，都會或多或少地造成地物光譜曲線失真，影響影像后續(xù)的解譯與應用［1－2］。為此，本文提出一種基于諧波分析的高光譜影像融合（harmonic analysis fusion，HAF）算法［3］，該算法是在分析并確定了光譜曲線諧波余相分量對光譜曲線波形形態(tài)無影響的基礎(chǔ)上，通過采用空間信息更為豐富的高空間分辨率影像替換諧波余相，然后對替換后的高空間分辨率影像與諧波振幅與相位進行諧波逆變換的一種高光譜融合算法，HAF算法在提高影像空間分辨率的同時完全保持光譜曲線的波形形態(tài)，但是由于直接采用高空間分辨率影像替換光譜曲線的諧波余項分量，導致原有的地物光譜曲線反射率受到了破壞。

文獻［4—7］為了對電力系統(tǒng)諧波污染進行有效地監(jiān)測與控制，提出了諧波分析算法；文獻［8—9］在其研究成果中首次將諧波分析理論應用到AVHRR和MODIS等遙感影像植被的NDVI時間序列分析當中；文獻［10］在其后續(xù)的諧波研究中將該算法應用到了農(nóng)作物物種識別當中；國內(nèi)也有學者利用諧波分析理論對山東省典型植被季節(jié)變化特征進行了分析［11］；文獻［12］則在諧波分析理論的輔助下對小目標地物進行了目標探測試驗；文獻［2—3］通過分析光譜曲線諧波各個特征分量的物理意義，將諧波分析理論應用到了影像融合當中，并取得了一定的研究成果，但由于初次嘗試，很多問題考慮不足，最終造成融合后地物光譜反射率整體出現(xiàn)過多偏移，地物光譜曲線物理解析特性遭到破壞，所以該算法還需進行改進與完善。

Gram－Schmidt（GS）變換融合算法可以較好地改善原始影像的空間細節(jié)特征，提高原始影像的空間分辨率，且能最大限度地保持原始影像的光譜物理特性［13－14］。鑒于此，本文將 GS變換引入HAF算法當中，提出了結(jié)合Gram－Schmidt變換的高光譜影像諧波分析融合（Gram－Schmidt transform combined with harmonic analysis fusion，GSHAF）改進算法。

1 理論與方法

1.1 高光譜影像的諧波分析法

文獻［15—16］認為，諧波分析是將目標的時間序列曲線f（t）表示為正（余）弦波相疊加的形式。諧波變換說明了任何連續(xù)的周期曲線f（t）都可以表示成為傅里葉級數(shù)形式［6］。高光譜遙感影像中，像元在各個波段的灰度值所組成的光譜曲線v（s）是一條連續(xù)的曲線，那么也可進行諧波分解，可將像元光譜曲線v（s）分解為諧波余項、諧波振幅與諧波相位相疊加的形式，其變換公式為

式中，s代表波段號數(shù)；A0/2表示諧波余項；Z代表波段總數(shù)；h表示諧波分析次數(shù)表示第h次諧波分量；φh表示第h次諧波的初相位；Ch表示第h次諧波的振幅，表示功率譜。

1.2 高光譜影像的諧波分析融合算法

根據(jù)式（1）可知，諧波余項（A0/2）分量計算結(jié)果為一個常數(shù)，該分量在曲線疊加過程中并不影響光譜曲線波形的形狀，只是地物對電磁波反射的綜合反應，其包含著影像豐富的空間信息。鑒于此，利用空間信息更為豐富的高空間分辨率影像直接替換諧波余項分量，如式（2）所示

式中，A代表高空間分辨率遙感影像，其他符號意義同式（1）。

隨后對替換后的諧波特征分量進行諧波逆變換，完成高光譜影像與高空間分辨率影像的融合，這便是HAF算法的基本原理。

至此，在保證高光譜影像光譜曲線波形不變的基礎(chǔ)上，將高光譜影像的融合問題轉(zhuǎn)換成了各像元光譜曲線的諧波余相特征分量組成的二維影像與高空間分辨率影像之間的融合。

1.3 Gram－Schmidt變換融合算法

Gram－Schmidt（GS）［13－14，17］變換是統(tǒng)計學中經(jīng)常用到的一種多維線性正交變換，采用GS變換對高光譜影像多維數(shù)據(jù)進行正交化處理，可以消除冗余信息。Gram－Schmidt變換融合的基本流程如下。

（1）首先利用原始低空間分辨率影像模擬出一幅全色影像。

（2）隨后利用該全色影像作為GS變換的第一個分量來對低空間分辨率影像進行GS變換，具體變換公式為

式中，GST表示經(jīng)GS變換后的第T個正交分量；BT表示原始低空間分辨率遙感影像第T波段；μT為原始低空間分辨率遙感影像第T波段像元灰度值的均值；φBK，GS（）l為原始低空間分辨率影像第K波段與GSl之間的協(xié)方差；i和j分別表示原始低空間分辨率影像的行數(shù)和列數(shù)；M和N表示整幅影像的行數(shù)和列數(shù)。

（3）用高空間分辨率影像替換GS變換后的第一分量，即GS1分量。

（4）最后再通過式（4）對上述替換后的數(shù)據(jù)集進行GS逆變換，完成低空間分辨率影像與高空間分辨率影像融合

式中，所有符號含義同式（3）。

1.4 結(jié)合Gram－Schmidt變換的諧波分析融合算法

HAF算法在保持融合影像光譜曲線波形形態(tài)方面極具優(yōu)勢，但是高光譜影像與高空間分辨率影像光譜范圍不同，所表示的光譜響應均值也不同，直接進行替換會對地物光譜曲線的反射率造成影響。針對該問題，本文提出結(jié)合Gram－Schmidt變換的諧波分析融合算法。GSHAF算法具體流程是在原始高光譜影像被諧波分析算法分解為諧波余項、諧波振幅和諧波相位相疊加的形式后，根據(jù)式（3）和式（4），首先對該光譜曲線的諧波余項分量與高空間分辨率影像進行GS變換融合，這樣可以利用GS變換融合算法的優(yōu)點既提高諧波余項的空間分辨率，同時也有效地修正了高空間分辨率影像與諧波相位像元灰度值之間巨大的差異，使GS變換融合后影像像元灰度值更接近諧波余項分量。隨后對經(jīng)GS變換融合后的影像與諧波振幅與諧波相位相兩個特征分量根據(jù)式（5）進行諧波逆變換，完成高光譜分辨率與高空間分辨率的影像融合，這樣便使融合后像元光譜曲線反射率更接近真實地物的光譜反射率，拓寬了融合影像后續(xù)應用范圍，增強了GSHAF算法的適用性。公式（5）計算方法為

式中，G代表原始高光譜影像像元光譜曲線的諧波余項與高空間分辨率遙感影像經(jīng)過GS變換融合后的影像，其他符號意義同式（1）。

2 試驗與分析

2.1 試驗數(shù)據(jù)

試驗數(shù)據(jù)為美國EO－1衛(wèi)星在2005年7月底獲取的Hyperion高光譜遙感影像與ALI高空間分辨率遙感影像全色波段影像，如圖1、圖2所示，其中圖1為該地區(qū)經(jīng)過假彩色合成后的Hyperion高光譜影像立體（RGB合成波段號為33、63、93），圖2中顯示的是其對應的ALI高空間分辨率全色波段影像。本次試驗采用ALI獲取的全色波段充當影像融合所需的高空間分辨率影像。

圖1 Hyperion高光譜影像Fig.1 Hyperion hyperspectral image

圖2 ALI全色波段影像Fig.2 ALI panchromatic band image

2.2 諧波分析融合處理

試驗數(shù)據(jù)經(jīng)過了一系列預處理操作，詳細流程見文獻［18—20］。待原始影像完成預處理之后，首先對原始Hyperion高光譜影像按照式（1）進行最佳分解次數(shù)的諧波分解［3］，將其分解為諧波余項、諧波振幅與諧波相位相疊加的表示形式，隨后根據(jù)式（3）、式（4）首先采用GS變換融合算法對原始高光譜影像的諧波余項分量與ALI高空間分辨率影像進行融合，最后再按照式（5）對融合影像與原始高光譜影像的諧波振幅與諧波相位進行諧波逆變換，完成高光譜影像與高空間分辨率影像的融合。融合結(jié)果如圖3所示的高光譜假彩色立體（RGB合成波段號為33、63、93），為了直觀地說明改進后的GSHAF算法比HAF算法具有更強的優(yōu)越性，本次試驗同樣對試驗數(shù)據(jù)進行HAF算法的融合處理，融合結(jié)果如圖4所示（RGB合成波段號為33、63、93），隨后對 HAF算法與GSHAF算法試驗結(jié)果進行對比分析。

圖3 GSHAF算法融合后Hyperion影像Fig.3 The fused Hyperion image based on the GSHAF

圖4 HAF算法融合后Hyperion影像Fig.4 The fused Hyperion image based on the HAF

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

通過目視評價融合后影像可知，GSHAF算法與HAF算法融合后的影像在空間分辨率方面均具有很大的改善，但是在影像色調(diào)方面GSHAF保持的更加貼近原始影像。為了更加客觀地分析經(jīng)過GSHAF算法與HAF算法融合后影像的光譜信息保持性情況，隨機地在兩種融合算法融合后影像上選擇植被和建筑物兩種地物的對應像元，并繪制每種地物像元的光譜曲線，如圖5、圖6所示，圖中虛線為HAF算法融合之后影像對應像元的光譜曲線，粗實線為GSHAF算法融合之后影像對應像元的光譜曲線，細實線為原始高光譜影像對應像元的光譜曲線。其中圖5為植被融合前后的光譜曲線，圖6表示的是建筑物融合前后的光譜曲線。從這兩個圖中可以看出，GSHAF算法與HAF算法在光譜波形形態(tài)方面都完全保留了原始影像的光譜曲線形態(tài)，但是GSHAF算法融合后的影像地物像元的光譜反射率與其原始的地物像元反射率更為接近，尤其在建筑物方面，幾乎完全與原始光譜曲線重合。以上分析客觀地證明了GSHAF算法在光譜信息保持方面更為優(yōu)秀，更符合一個優(yōu)秀融合算法的要求，其融合后的影像具有更為廣泛的應用前景。

圖5 融合前后植被光譜曲線Fig.5 The vegetation spectra of original image and fused image

圖6 融合前后建筑物光譜曲線Fig.6 The building spectra of original image and fused image

2.4 GSHAF算法的普適性驗證

好的融合算法不僅需要具有優(yōu)良的可行性，同時也需要具有良好的普適性，為了驗證GSHAF算法是否對其他數(shù)據(jù)源也具有良好的可行性，本節(jié)僅采用國產(chǎn)HJ－1A衛(wèi)星數(shù)據(jù)為例，利用該影像數(shù)據(jù)對本文提出的GSHAF算法進行普適性驗證。本次普適性驗證試驗數(shù)據(jù)采用2013年4月25日獲取的陜西某城市郊區(qū)及其附近山區(qū)的影像，如圖7所示為經(jīng)過假彩色合成后的該區(qū)域高光譜影像立體（RGB合成波段號為60、90、110），本次試驗采用多光譜遙感影像的第3波段充當高空間分辨率影像，如圖8所示。

圖7 HJ－1A高光譜影像Fig.7 HJ－1Ahyperspectral image

圖8 HJ－1A 衛(wèi) 星 CCD第3波段影像Fig.8The third band image of HJ－1ACCD

同樣對HJ－1A衛(wèi)星的數(shù)據(jù)進行預處理工作，具體流程參閱文獻［21—22］，隨后分別進行GSHAF與HAF兩種融合算法的融合，其結(jié)果如圖9、10所示的高光譜假彩色立體（RGB合成波段號為60、90、110）。通過目視評價可知，經(jīng)過GSHAF算法與HAF算法融合后的影像空間分辨率都獲得了很大的提高，在其色調(diào)方面兩幅融合后影像整體差異不大，但是仔細觀察融合前后影像中的河流，還是可以發(fā)現(xiàn)GSHAF算法融合后的影像色調(diào)更接近原始影像。

圖9 GSHAF算法融合HJ－1A 影像Fig.9 The fused HJ－1A image based on the GSHAF

圖10 HAF算法融合HJ－1A 影像Fig.10 The fused HJ－1A image based on the HAF

為了客觀了解GSHAF算法在光譜信息保持性方面的情況，同樣繪制出融合前后影像中對應位置隨機像元A、B的光譜曲線，如圖11、圖12所示，圖中3條曲線所代表的影像像元意義與圖5、6相同?？梢悦黠@看出，該算法所表示出來的結(jié)果與可行性分析試驗完全相同，光譜曲線波形形態(tài)兩種方法都完全地保留下來，反射率方面GSHAF算法融合影像更為接近原始影像，圖12所示GSHAF算法融合后的光譜曲線也幾乎與原始影像光譜曲線重合，以上分析客觀真實地證明了GSHAF算法不僅具有優(yōu)良的可行性，對其他數(shù)據(jù)還具有良好的普適性。

圖11 融合前后隨機像元A光譜曲線Fig.11 The A random pixel spectra of original image and fused image

圖12 融合前后隨機像元B光譜曲線Fig.12 The B random pixel spectra of original image and fused image

3 結(jié) 論

本文在分析基于諧波分析理論的高光譜遙感影像融合算法缺陷的基礎(chǔ)上，提出了結(jié)合Gram－Schmidt變換的諧波分析融合改進算法，該算法可以在完全保留融合前后影像光譜曲線波形形態(tài)的基礎(chǔ)上，有效地減少HAF算法融合后像元光譜曲線與原始光譜曲線反射率差異較大、地物像元光譜物理特性失真這一缺陷，使融合后影像獲得更廣泛的用途，并通過試驗來證明這一假設。通過試驗與分析可得：①HAF算法和GSHAF算法在影像融合過程中都可以完全地保留原始影像的光譜曲線波形形態(tài)；②GSHAF算法在光譜物理特性保持方面比HAF算法更具優(yōu)勢，融合后地物的光譜反射率值更接近原始影像地物光譜反射率值；③GSHAF算法不僅具有優(yōu)良的可行性，還具有良好的普適性；④雖然GS變換對HAF算法具有一定的改進作用，但是還是無法完全地保留地物光譜全部的物理特性，然而，本文已經(jīng)在保證融合前后光譜曲線波形不變的情況下，將高光譜遙感影像融合簡化，只要找到合適的算法，后續(xù)的研究中將會進一步將地物光譜信息保留下來。

［1］DONG Guangjun，ZHANG Yongsheng，F(xiàn)AN Yonghong.Image Fusion for Hyperspectral Data of PHI and High Resolution Aerial Image［J］.Journal of Infrared and Millimeter Waves，2006，25（2）：123－126.（董廣軍，張永生，范永弘.PHI高光譜數(shù)據(jù)和高空間分辨率遙感圖像融合技術(shù)研究［J］.紅外與毫米波學報，2006，25（2）：123－126.）

［2］YANG Keming，ZHANG Tao，WANG Libo，et al.Harmonic Analysis Fusion of Hyperspectral Image and Its Spectral Information Fidelity Evaluation［J］.Spectroscopy and Spectral Analysis，2013，38，33（9）：2498－2501.（楊可明，張濤，王立博，等.諧波分析法高光譜影像融合及其光譜信息保真度評價［J］.光譜學與光譜分析，2013，33（9）：2498－2501.）

［3］YANG Keming，ZHANG Tao，WANG Libo，et al.A New Algorithm on Hyperspectral Image Fusion Based on the Harmonic Analysis［J］.Journal of China University of Mining ＆ Technology，2014，43（3）：547－553.（楊可明，張濤，王立博，等.高光譜影像的諧波分析融合算法研究［J］.中國礦業(yè)大學學報，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］REN Zihui，QIU Runhe，ZHANG Yan.Formation of Power System Harmonic Model of Mine and Design of Harmonic Filter for Mine Hoister［J］.Journal of China University of Mining ＆ Technology，2004，33（1）：45－49.（任子暉，仇潤鶴，張艷.煤礦電網(wǎng)諧波分析模型的建立與濾波器設計［J］.中國礦業(yè)大學學報，2004，33（1）：45－49.）

［6］GEORGE J W.Power Systems Harmonics Fundamentals，Analysis and Filter Design［M］.Translated by XU Zheng.Beijing：China Machine Press，2003：4－13.（GEORGE J W.電力系統(tǒng)諧波——基本原理、分析方法和濾波器設計［M］.徐政，譯.北京：機械工業(yè)出版社，2003：4－13.）

［7］CHAI Xuzheng，WEN Xishan，GUAN Genzhi，et al.An Algorithm with High Accuracy for Analysis of Power System Harmonics［J］.Proceedings of the CSEE，2003，23（9）：68－71.（柴旭崢，文習山，關(guān)根志，等.一種高精度的電力系統(tǒng)諧波分析算法［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（3－4）：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（2）：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（1－3）：127－139.

［11］LIANG Shouzhen，XING Qianguo，SHI Ping，et al.Harmonic Analysis on NDVI Time Series of Typical Land Covers in Shandong Province［J］.Chinese Journal of Ecology，2011，30（1）：59－65.（梁守真，邢前國，施平，等.山東省典型地表覆被NDVI時間序列諧波分析［J］.生態(tài)學雜志，2011，30（1）：59－65.）

［12］YANG Keming，XUE Zhaohui，JIA Taotao，et al.A Harmonic Analysis Model of Small Target Detection of Hyperspectral Imagery［J］.Acta Geodaetica et Cartographica Sinica，2013，42（1）：34－43.（楊可明，薛朝輝，賈濤濤，等.高光譜影像小目標諧波分析探測模型［J］.測繪學報，2013，42（1）：34－43.）

［13］LI Cunjun，LIU Liangyun，WANG Jihua，et al.Comparison of Two Methods of Fusing Remote Sensing Images with Fidelity of Spectral Information［J］.Journal of Image and Graphics，2004，9（11）：1376－1385.（李存軍，劉良云，王紀華，等.兩種高保真遙感影像融合方法比較［J］.中國圖象圖形學報，2004，9（11）：1376－1385.）

［14］YU H aiyang，YA N Bokun，GA N Fuping，et al.Hyperspectral Image Fusion by an Enhanced Gram Schmidt Spectral Transformation［J］.Geography and Geoinformation Science，2007，23（5）：30－42.（于海洋，閆柏琨，甘甫平，等.基于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［M］.2nd ed.New York：John Wiley ＆Sons，1986：23－34.

［17］HUANG Dengshan.Research of Pixel－level Remote sensing Image Fusion Method［D］.Changsha：Central South University，2011：18.（黃登山.像素級遙感影像融合方法研究［D］.長沙：中南大學，2011：18.）

［18］TAN Bingxiang，LI Zengyuan，CHEN Erxue，et al.Preprocessing of EO－1Hyperion Hyperspectral Data［J］.Remote Sensing Information，2005（6）：36－41.（譚炳香，李增元，陳爾學，等.EO－1Hyperion高光譜數(shù)據(jù)的預處理［J］.遙感信息，2005（6）：36－41.）

［19］CHEN Jianzhen，HE Chao，YUE Cairong.Atmospheric Correction of an Advance Land Imager（ALI）Image Based on the FLAASH Module［J］.Journal of Zhejiang A＆F University，2011，28（4）：590－596.（陳建珍，何超，岳彩榮.基于FLAASH模塊的高級陸地成像儀圖像的大氣校正［J］.浙江農(nóng)林大學學報，2011，28（4）：590－596.）

［20］CHANDER G，MARKHAM B L，HELDER D L.Summary of Current Radiometric Calibration Coefficients for Landsat MSS，TM，ETM＋，and EO－1ALI Sensors［J］.Remote Sensing of Environment，2009，113（5）：893－903.

［21］MA Lingling，WANG Xinhong，TANG Lingli.A Highly Efficient Atmospheric Correction Method for HJ－1A/HSI and the Exploration on Its Application Capability［J］.Remote Sensing Technology and Application，2010，25（4）：525－531.（馬靈玲，王新鴻，唐伶俐.HJ－1A高光譜數(shù)據(jù)高效大氣校正及應用潛力初探［J］.遙感技術(shù)與應用，2010，25（4）：525－531.）

［22］ZHENG Sheng，ZHAO Xiang，ZHANG Hao，et al.Atmospheric Correction on CCD Data of HJ－1Satellite and Analysis of Its Effect［J］.Journal of Remote Sensing，2011，15（4）：709－721.（鄭盛，趙祥，張顥，等.HJ－1衛(wèi)星CCD數(shù)據(jù)的大氣校正及其效果分析［J］.遙感學報，2011，15（4）：709－721.）