張 濤,劉 軍,楊可明,羅文杉,張育育
1.國(guó)家測(cè)繪地理信息局第一航測(cè)遙感院,陜西 西安710054;2.中國(guó)礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院,北京100083
現(xiàn)有的高光譜影像融合算法在提高影像空間分辨率的同時(shí),都會(huì)或多或少地造成地物光譜曲線失真,影響影像后續(xù)的解譯與應(yīng)用[1-2]。為此,本文提出一種基于諧波分析的高光譜影像融合(harmonic analysis fusion,HAF)算法[3],該算法是在分析并確定了光譜曲線諧波余相分量對(duì)光譜曲線波形形態(tài)無(wú)影響的基礎(chǔ)上,通過(guò)采用空間信息更為豐富的高空間分辨率影像替換諧波余相,然后對(duì)替換后的高空間分辨率影像與諧波振幅與相位進(jìn)行諧波逆變換的一種高光譜融合算法,HAF算法在提高影像空間分辨率的同時(shí)完全保持光譜曲線的波形形態(tài),但是由于直接采用高空間分辨率影像替換光譜曲線的諧波余項(xiàng)分量,導(dǎo)致原有的地物光譜曲線反射率受到了破壞。
文獻(xiàn)[4—7]為了對(duì)電力系統(tǒng)諧波污染進(jìn)行有效地監(jiān)測(cè)與控制,提出了諧波分析算法;文獻(xiàn)[8—9]在其研究成果中首次將諧波分析理論應(yīng)用到AVHRR和MODIS等遙感影像植被的NDVI時(shí)間序列分析當(dāng)中;文獻(xiàn)[10]在其后續(xù)的諧波研究中將該算法應(yīng)用到了農(nóng)作物物種識(shí)別當(dāng)中;國(guó)內(nèi)也有學(xué)者利用諧波分析理論對(duì)山東省典型植被季節(jié)變化特征進(jìn)行了分析[11];文獻(xiàn)[12]則在諧波分析理論的輔助下對(duì)小目標(biāo)地物進(jìn)行了目標(biāo)探測(cè)試驗(yàn);文獻(xiàn)[2—3]通過(guò)分析光譜曲線諧波各個(gè)特征分量的物理意義,將諧波分析理論應(yīng)用到了影像融合當(dāng)中,并取得了一定的研究成果,但由于初次嘗試,很多問(wèn)題考慮不足,最終造成融合后地物光譜反射率整體出現(xiàn)過(guò)多偏移,地物光譜曲線物理解析特性遭到破壞,所以該算法還需進(jìn)行改進(jìn)與完善。
Gram-Schmidt(GS)變換融合算法可以較好地改善原始影像的空間細(xì)節(jié)特征,提高原始影像的空間分辨率,且能最大限度地保持原始影像的光譜物理特性[13-14]。鑒于此,本文將 GS變換引入HAF算法當(dāng)中,提出了結(jié)合Gram-Schmidt變換的高光譜影像諧波分析融合(Gram-Schmidt transform combined with harmonic analysis fusion,GSHAF)改進(jìn)算法。
文獻(xiàn)[15—16]認(rèn)為,諧波分析是將目標(biāo)的時(shí)間序列曲線f(t)表示為正(余)弦波相疊加的形式。諧波變換說(shuō)明了任何連續(xù)的周期曲線f(t)都可以表示成為傅里葉級(jí)數(shù)形式[6]。高光譜遙感影像中,像元在各個(gè)波段的灰度值所組成的光譜曲線v(s)是一條連續(xù)的曲線,那么也可進(jìn)行諧波分解,可將像元光譜曲線v(s)分解為諧波余項(xiàng)、諧波振幅與諧波相位相疊加的形式,其變換公式為
式中,s代表波段號(hào)數(shù);A0/2表示諧波余項(xiàng);Z代表波段總數(shù);h表示諧波分析次數(shù)表示第h次諧波分量;φh表示第h次諧波的初相位;Ch表示第h次諧波的振幅,表示功率譜。
根據(jù)式(1)可知,諧波余項(xiàng)(A0/2)分量計(jì)算結(jié)果為一個(gè)常數(shù),該分量在曲線疊加過(guò)程中并不影響光譜曲線波形的形狀,只是地物對(duì)電磁波反射的綜合反應(yīng),其包含著影像豐富的空間信息。鑒于此,利用空間信息更為豐富的高空間分辨率影像直接替換諧波余項(xiàng)分量,如式(2)所示
式中,A代表高空間分辨率遙感影像,其他符號(hào)意義同式(1)。
隨后對(duì)替換后的諧波特征分量進(jìn)行諧波逆變換,完成高光譜影像與高空間分辨率影像的融合,這便是HAF算法的基本原理。
至此,在保證高光譜影像光譜曲線波形不變的基礎(chǔ)上,將高光譜影像的融合問(wèn)題轉(zhuǎn)換成了各像元光譜曲線的諧波余相特征分量組成的二維影像與高空間分辨率影像之間的融合。
Gram-Schmidt(GS)[13-14,17]變換是統(tǒng)計(jì)學(xué)中經(jīng)常用到的一種多維線性正交變換,采用GS變換對(duì)高光譜影像多維數(shù)據(jù)進(jìn)行正交化處理,可以消除冗余信息。Gram-Schmidt變換融合的基本流程如下。
(1)首先利用原始低空間分辨率影像模擬出一幅全色影像。
(2)隨后利用該全色影像作為GS變換的第一個(gè)分量來(lái)對(duì)低空間分辨率影像進(jìn)行GS變換,具體變換公式為
式中,GST表示經(jīng)GS變換后的第T個(gè)正交分量;BT表示原始低空間分辨率遙感影像第T波段;μT為原始低空間分辨率遙感影像第T波段像元灰度值的均值;φBK,GS()l為原始低空間分辨率影像第K波段與GSl之間的協(xié)方差;i和j分別表示原始低空間分辨率影像的行數(shù)和列數(shù);M和N表示整幅影像的行數(shù)和列數(shù)。
(3)用高空間分辨率影像替換GS變換后的第一分量,即GS1分量。
(4)最后再通過(guò)式(4)對(duì)上述替換后的數(shù)據(jù)集進(jìn)行GS逆變換,完成低空間分辨率影像與高空間分辨率影像融合
式中,所有符號(hào)含義同式(3)。
HAF算法在保持融合影像光譜曲線波形形態(tài)方面極具優(yōu)勢(shì),但是高光譜影像與高空間分辨率影像光譜范圍不同,所表示的光譜響應(yīng)均值也不同,直接進(jìn)行替換會(huì)對(duì)地物光譜曲線的反射率造成影響。針對(duì)該問(wèn)題,本文提出結(jié)合Gram-Schmidt變換的諧波分析融合算法。GSHAF算法具體流程是在原始高光譜影像被諧波分析算法分解為諧波余項(xiàng)、諧波振幅和諧波相位相疊加的形式后,根據(jù)式(3)和式(4),首先對(duì)該光譜曲線的諧波余項(xiàng)分量與高空間分辨率影像進(jìn)行GS變換融合,這樣可以利用GS變換融合算法的優(yōu)點(diǎn)既提高諧波余項(xiàng)的空間分辨率,同時(shí)也有效地修正了高空間分辨率影像與諧波相位像元灰度值之間巨大的差異,使GS變換融合后影像像元灰度值更接近諧波余項(xiàng)分量。隨后對(duì)經(jīng)GS變換融合后的影像與諧波振幅與諧波相位相兩個(gè)特征分量根據(jù)式(5)進(jìn)行諧波逆變換,完成高光譜分辨率與高空間分辨率的影像融合,這樣便使融合后像元光譜曲線反射率更接近真實(shí)地物的光譜反射率,拓寬了融合影像后續(xù)應(yīng)用范圍,增強(qiáng)了GSHAF算法的適用性。公式(5)計(jì)算方法為
式中,G代表原始高光譜影像像元光譜曲線的諧波余項(xiàng)與高空間分辨率遙感影像經(jīng)過(guò)GS變換融合后的影像,其他符號(hào)意義同式(1)。
試驗(yàn)數(shù)據(jù)為美國(guó)EO-1衛(wèi)星在2005年7月底獲取的Hyperion高光譜遙感影像與ALI高空間分辨率遙感影像全色波段影像,如圖1、圖2所示,其中圖1為該地區(qū)經(jīng)過(guò)假彩色合成后的Hyperion高光譜影像立體(RGB合成波段號(hào)為33、63、93),圖2中顯示的是其對(duì)應(yīng)的ALI高空間分辨率全色波段影像。本次試驗(yàn)采用ALI獲取的全色波段充當(dāng)影像融合所需的高空間分辨率影像。
圖1 Hyperion高光譜影像Fig.1 Hyperion hyperspectral image
圖2 ALI全色波段影像Fig.2 ALI panchromatic band image
試驗(yàn)數(shù)據(jù)經(jīng)過(guò)了一系列預(yù)處理操作,詳細(xì)流程見(jiàn)文獻(xiàn)[18—20]。待原始影像完成預(yù)處理之后,首先對(duì)原始Hyperion高光譜影像按照式(1)進(jìn)行最佳分解次數(shù)的諧波分解[3],將其分解為諧波余項(xiàng)、諧波振幅與諧波相位相疊加的表示形式,隨后根據(jù)式(3)、式(4)首先采用GS變換融合算法對(duì)原始高光譜影像的諧波余項(xiàng)分量與ALI高空間分辨率影像進(jìn)行融合,最后再按照式(5)對(duì)融合影像與原始高光譜影像的諧波振幅與諧波相位進(jìn)行諧波逆變換,完成高光譜影像與高空間分辨率影像的融合。融合結(jié)果如圖3所示的高光譜假彩色立體(RGB合成波段號(hào)為33、63、93),為了直觀地說(shuō)明改進(jìn)后的GSHAF算法比HAF算法具有更強(qiáng)的優(yōu)越性,本次試驗(yàn)同樣對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行HAF算法的融合處理,融合結(jié)果如圖4所示(RGB合成波段號(hào)為33、63、93),隨后對(duì) HAF算法與GSHAF算法試驗(yàn)結(jié)果進(jìn)行對(duì)比分析。
圖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
通過(guò)目視評(píng)價(jià)融合后影像可知,GSHAF算法與HAF算法融合后的影像在空間分辨率方面均具有很大的改善,但是在影像色調(diào)方面GSHAF保持的更加貼近原始影像。為了更加客觀地分析經(jīng)過(guò)GSHAF算法與HAF算法融合后影像的光譜信息保持性情況,隨機(jī)地在兩種融合算法融合后影像上選擇植被和建筑物兩種地物的對(duì)應(yīng)像元,并繪制每種地物像元的光譜曲線,如圖5、圖6所示,圖中虛線為HAF算法融合之后影像對(duì)應(yīng)像元的光譜曲線,粗實(shí)線為GSHAF算法融合之后影像對(duì)應(yīng)像元的光譜曲線,細(xì)實(shí)線為原始高光譜影像對(duì)應(yīng)像元的光譜曲線。其中圖5為植被融合前后的光譜曲線,圖6表示的是建筑物融合前后的光譜曲線。從這兩個(gè)圖中可以看出,GSHAF算法與HAF算法在光譜波形形態(tài)方面都完全保留了原始影像的光譜曲線形態(tài),但是GSHAF算法融合后的影像地物像元的光譜反射率與其原始的地物像元反射率更為接近,尤其在建筑物方面,幾乎完全與原始光譜曲線重合。以上分析客觀地證明了GSHAF算法在光譜信息保持方面更為優(yōu)秀,更符合一個(gè)優(yōu)秀融合算法的要求,其融合后的影像具有更為廣泛的應(yīng)用前景。
圖5 融合前后植被光譜曲線Fig.5 The vegetation spectra of original image and fused image
圖6 融合前后建筑物光譜曲線Fig.6 The building spectra of original image and fused image
好的融合算法不僅需要具有優(yōu)良的可行性,同時(shí)也需要具有良好的普適性,為了驗(yàn)證GSHAF算法是否對(duì)其他數(shù)據(jù)源也具有良好的可行性,本節(jié)僅采用國(guó)產(chǎn)HJ-1A衛(wèi)星數(shù)據(jù)為例,利用該影像數(shù)據(jù)對(duì)本文提出的GSHAF算法進(jìn)行普適性驗(yàn)證。本次普適性驗(yàn)證試驗(yàn)數(shù)據(jù)采用2013年4月25日獲取的陜西某城市郊區(qū)及其附近山區(qū)的影像,如圖7所示為經(jīng)過(guò)假彩色合成后的該區(qū)域高光譜影像立體(RGB合成波段號(hào)為60、90、110),本次試驗(yàn)采用多光譜遙感影像的第3波段充當(dāng)高空間分辨率影像,如圖8所示。
圖7 HJ-1A高光譜影像Fig.7 HJ-1Ahyperspectral image
圖8 HJ-1A 衛(wèi) 星 CCD第3波段影像Fig.8The third band image of HJ-1ACCD
同樣對(duì)HJ-1A衛(wèi)星的數(shù)據(jù)進(jìn)行預(yù)處理工作,具體流程參閱文獻(xiàn)[21—22],隨后分別進(jìn)行GSHAF與HAF兩種融合算法的融合,其結(jié)果如圖9、10所示的高光譜假彩色立體(RGB合成波段號(hào)為60、90、110)。通過(guò)目視評(píng)價(jià)可知,經(jīng)過(guò)GSHAF算法與HAF算法融合后的影像空間分辨率都獲得了很大的提高,在其色調(diào)方面兩幅融合后影像整體差異不大,但是仔細(xì)觀察融合前后影像中的河流,還是可以發(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算法在光譜信息保持性方面的情況,同樣繪制出融合前后影像中對(duì)應(yīng)位置隨機(jī)像元A、B的光譜曲線,如圖11、圖12所示,圖中3條曲線所代表的影像像元意義與圖5、6相同。可以明顯看出,該算法所表示出來(lái)的結(jié)果與可行性分析試驗(yàn)完全相同,光譜曲線波形形態(tài)兩種方法都完全地保留下來(lái),反射率方面GSHAF算法融合影像更為接近原始影像,圖12所示GSHAF算法融合后的光譜曲線也幾乎與原始影像光譜曲線重合,以上分析客觀真實(shí)地證明了GSHAF算法不僅具有優(yōu)良的可行性,對(duì)其他數(shù)據(jù)還具有良好的普適性。
圖11 融合前后隨機(jī)像元A光譜曲線Fig.11 The A random pixel spectra of original image and fused image
圖12 融合前后隨機(jī)像元B光譜曲線Fig.12 The B random pixel spectra of original image and fused image
本文在分析基于諧波分析理論的高光譜遙感影像融合算法缺陷的基礎(chǔ)上,提出了結(jié)合Gram-Schmidt變換的諧波分析融合改進(jìn)算法,該算法可以在完全保留融合前后影像光譜曲線波形形態(tài)的基礎(chǔ)上,有效地減少HAF算法融合后像元光譜曲線與原始光譜曲線反射率差異較大、地物像元光譜物理特性失真這一缺陷,使融合后影像獲得更廣泛的用途,并通過(guò)試驗(yàn)來(lái)證明這一假設(shè)。通過(guò)試驗(yàn)與分析可得:①HAF算法和GSHAF算法在影像融合過(guò)程中都可以完全地保留原始影像的光譜曲線波形形態(tài);②GSHAF算法在光譜物理特性保持方面比HAF算法更具優(yōu)勢(shì),融合后地物的光譜反射率值更接近原始影像地物光譜反射率值;③GSHAF算法不僅具有優(yōu)良的可行性,還具有良好的普適性;④雖然GS變換對(duì)HAF算法具有一定的改進(jìn)作用,但是還是無(wú)法完全地保留地物光譜全部的物理特性,然而,本文已經(jīng)在保證融合前后光譜曲線波形不變的情況下,將高光譜遙感影像融合簡(jiǎn)化,只要找到合適的算法,后續(xù)的研究中將會(huì)進(jìn)一步將地物光譜信息保留下來(lái)。
[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].紅外與毫米波學(xué)報(bào),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.(楊可明,張濤,王立博,等.諧波分析法高光譜影像融合及其光譜信息保真度評(píng)價(jià)[J].光譜學(xué)與光譜分析,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].中國(guó)礦業(yè)大學(xué)學(xué)報(bào),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.(任子暉,仇潤(rùn)鶴,張艷.煤礦電網(wǎng)諧波分析模型的建立與濾波器設(shè)計(jì)[J].中國(guó)礦業(yè)大學(xué)學(xué)報(bào),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)諧波——基本原理、分析方法和濾波器設(shè)計(jì)[M].徐政,譯.北京:機(jī)械工業(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.(柴旭崢,文習(xí)山,關(guān)根志,等.一種高精度的電力系統(tǒng)諧波分析算法[J].中國(guó)電機(jī)工程學(xué)報(bào),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.(梁守真,邢前國(guó),施平,等.山東省典型地表覆被NDVI時(shí)間序列諧波分析[J].生態(tài)學(xué)雜志,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.(楊可明,薛朝輝,賈濤濤,等.高光譜影像小目標(biāo)諧波分析探測(cè)模型[J].測(cè)繪學(xué)報(bào),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ì)華,等.兩種高保真遙感影像融合方法比較[J].中國(guó)圖象圖形學(xué)報(bào),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ìn)融合方法[J].地理與地理信息科學(xué),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.(黃登山.像素級(jí)遙感影像融合方法研究[D].長(zhǎng)沙:中南大學(xué),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.(譚炳香,李增元,陳爾學(xué),等.EO-1Hyperion高光譜數(shù)據(jù)的預(yù)處理[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í)陸地成像儀圖像的大氣校正[J].浙江農(nóng)林大學(xué)學(xué)報(bào),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ù)高效大氣校正及應(yīng)用潛力初探[J].遙感技術(shù)與應(yīng)用,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].遙感學(xué)報(bào),2011,15(4):709-721.)