兩種改進小波算法的衛星多波段數據融合①

許 晨,康 雪,張春同,徐 洋,呂達仁

1(成都市氣象局,成都 611130)

2(四川省氣象局,成都 610072)

3(中國科學院 大氣物理研究所,北京 100029)

基于衛星多波段數據融合算法,將多傳感器和多信息源的數據、信息加以聯合、相關和組合,以獲得精確的位置估計和一致性估計.基本策略是由低層到高層對多源信息進行整合、抽象的信息處理過程[1].通過多源信息的互補,消除觀測對象的不確定性、減小數據冗余,可改善待觀測對象的精確性和定量可靠性,從而有效提高數據利用率.

現有融合算法主要有以下幾類[1]:

(1)彩色空間變換融合法.將低分辨率的RGB 影像數據經過變換映射至HIS 空間,然后采用特定的融合算法使其與高分辨率數據進行融合處理,進而置換相應的部分,最后經過逆變換重構融合數據.根據映射的顏色空間不同,彩色空間融合法可分為HIS 變換、YIQ 變換、HLS 變換融合法等.

(2)加權融合或基于信息量的融合.可根據經驗值或相關系數設置權重函數,以減少冗余.此算法的優點在于簡單易行,可結合基于特征的融合方法,針對不同要求,靈活改變信息特征提取方法.

(3)高通濾波融合(High-Pass Filter,HPF).通過將高分辨率數據中的幾何信息逐像素疊加到低分辨率數據中而進行.先對高分辨率的全色數據和低分辨率的多光譜數據各波段進行直方圖匹配,對匹配后的全色數據進行高通濾波,再將其加入多光譜的各個波段,最后將多光譜各個波段的數據進行彩色合成.

(4)主成分變換融合(Principal Component Analysis,PCA).通過該變換,使多光譜數據在各波段具有統計獨立性,便于在各波段采用相應的融合策略.

(5)小波融合(Wavelet Transforms,WT).小波變換作為新興的數學分析方法日益受到廣泛重視,是分析和處理非平穩信號的一種有效工具.小波變換以局部化函數所形成的小波基作為基底而展開,是一種窗口大小固定不變但其形狀可改變的時頻局部化分析方法.它已在圖像編碼領域、計算機視覺、語音信號處理等領域取得了突破性進展.

在小波融合算法中,最常用到二進制小波變換,常用來處理來自同一傳感器的遙感數據,通過對數據進行分解、重構,在高、低頻采用相應的融合策略來實現數據融合[2].但在實際應用中,常需要融合來自不同傳感器的數據,需要擴展到多進制小波算法,但為保證融合處理速度,小波分解的階數取得并不高,這時融合結果空間細節的表現受到影響.考慮到HSV 融合法能突出空間信息表達能力、小波包融合法能克服小波變換對高頻信息處理的缺陷,所以構建兩種改進算法:基于HSV的小波融合法和基于區域特征的自適應小波包融合算法分別實現不同傳感器的數據融合.

1 基于HSV的多進制小波融合(HSV-WT)

多進制小波可理解為頻率域的分解問題,二進制小波把頻率域分解成兩個通道,多進制小波把頻率域分解成多個通道,分別對每個分解層次的低、高頻部分按各自的融合策略進行處理,綜合各組數據的特征信息,最后根據多進制小波重構公式進行小波逆變換,重構融合數據[3].因此它能突破待融合數據的分辨率比值限制,實現分辨率之比非2n的數據融合.

1.1 HSV 融合法

HSV 分別代表 Hue (色調)、Saturation (飽和度)、Value (亮度),融合模型為一圓錐體,如圖1所示,圓錐底面對應色調H,表示所處的顏色位置,以繞V 軸的旋轉角度來表示不同顏色,0°對應紅色,120°對應綠色,240°對應藍色,互補色之間相差180°;飽和度S從低到高表示為由軸心向錐體圓周過渡,表示所選顏色的純度和該色最大純度之間的比率,范圍為0-1.當S=0 時,只有灰度;明度V表示色彩明亮程度,范圍為0-1[4,5].

圖1 HSV 變換模型

從RGB 空間轉換到HSV 空間[6]:首先,歸一化RGB 色彩空間的值,求其最大值和最小值.

令m=MAX(R,G,B),n=MIN(R,G,B),其轉換公式為式(1).

1.2 數據配準

由于各傳感器通過的光路不同或成像機制不同等原因,多源遙感數據間可能出現相對平移、旋轉或比例縮放等現象,必須先對圖像進行配準.配準的關鍵在于從多個圖像中找到具有共性特征的控制點,包括相對配準和絕對配準[1].

(1)相對配準:先在兩幅待配準圖像上選擇同名控制點,用二次多項式模型建立兩個同名像素的關系,最后重采樣成相同分辨率的圖像.

(2)絕對配準:在統一地理坐標系下對待配準圖像進行幾何糾正,常用多項式糾正法,再重采樣為相同分辨率圖像.

多項式糾正法的關鍵是采用一定數量的具有空間坐標的地面控制點,利用這些點對應的坐標,通過平差原理計算多項式系數,數學模型如下[7]:

式中,(x,y)是圖像中的坐標,(X,Y)是對應的地面坐標,多項式控制點個數N與多項式次數n之間的關系為:

根據上式計算多項式系數:a0,a1,a2,a3,···,通過變換關系式計算原影像上的坐標,將該點亮度值替換為輸出影像的坐標[8].

1.3 融合算法實現

將兩種算法相結合,具體處理流程如圖2所示[9-11].

圖2 基于HSV的小波融合法實現框圖

(1)將待融合的高分辨率數據I1和多光譜數據I2配準,用三次卷積內插法采樣.在配準過程中,先通過控制點位法進行相對配準,再通過多項式糾正法進行絕對配準.

(2)將低分辨率多光譜數據I2由RGB 空間轉換到HSV 空間,得到數據I2的H,S,V三分量信息.

(3)將V分量和高分辨率數據I1進行直方圖匹配,使V分量數據和I1的幅度值保持一致.

(4)利用多進制小波融合算法對V分量信息和高分辨率數據I1進行小波融合從而得到新的V'分量.

(5)將得到的V'分量和前面的H、S分量數據進行HSV 逆變換,輸出高分辨率的融合后數據.

2 基于區域特征的自適應小波包融合算法(AWP)

小波包融合算法通過對不同分辨率的高頻部分進一步劃分,進行遞歸分解,突出高分辨率數據的細節層次信息.該算法能突出細節區域特征,且采用自適應算法,故稱之為基于區域特征的自適應小波包融合算法[12](Adaptive Wavelet Packet based on region features,AWP).

2.1 算法介紹

多分辨率分析可以對信號進行有效的時效分解,但由于其尺度函數時按二進制變化的,因此在高頻段其頻率分辨率較差,只能對信號的頻段進行指數等間隔劃分.小波包分解為信號提供一種更加精細的分解方法,通過把頻帶進行多層次劃分,對多分辨分析中沒有細分的高頻部分進一步分解,并根據被分析信號特征,自適應地選擇相應頻段,使之與信號頻譜相匹配,從而提高時頻分辨率,因此具有更廣泛的應用價值[13].

這種算法可同時對數據的高、低頻部分進行任意尺度的遞歸分解和融合處理,克服了傳統小波融合算法對高分辨率數據高頻細節信息處理的不足,突出高分辨率數據細節層次信息,獲取高分辨率的多光譜融合數據,從而提高數據可用度和研究對象解譯的可靠性[1].

2.2 算法實現

首先對配準的遙感數據進行小波包分解,分解時,對上一層的各高、低頻部分進行全方位遞歸分解,采用基于自適應能力和區域特征的融合策略對子數據融合,最后逆變換重構數據.一般情況下,小波包分解層次越多,融合結果中包含的細節信息就越豐富.但隨著分解層次的增多,計算量增加,而且易造成頂層融合數據損失增大[14].因此折中考慮兩點,一般分解層次取3-5之間.本文的分解層次取3,Haar 小波基.

如圖3,基于自適應能力和區域特征的融合策略為:按照不同融合策略分別對每一級小波系數上的高、低頻信息進行處理,并合并系數.加權平均法的基本原理是[1]:

其中,CA,CB,CF分別表示源數據A、B及融合后數據F的在最后一個分解層的子數據,m,n為最后一個分解層子數據的像素位置,α為權值,值在0-1 之間.

圖3 融合策略

平均梯度法的融合規則為[1]:利用該像素的局部平均梯度確定融合后高頻子數據的像素值.設待融合數據為A(x,y),B(x,y),由不同分辨率的高頻子數據得到的梯度數據分別為(x,y),(x,y),則不同分辨率上的高頻子數據(x,y)為:

最后使用重構濾波器逆向逐層二插值重構數據,獲取融合后數據.

算法步驟[15,16]如圖4.

圖4 基于區域特征的自適應小波包融合法

算法步驟具體說明如下:

(1)輸入高分辨率的全色數據和多光譜數據.

(2)初始化分解系數、分解層次數目等參數.

(3)設計尺度函數? (x)對應的低通分解濾波器和小波函數 ψ (x)對應的高通分解濾波器,以及相應的重構濾波器.

(4)根據小波包分解算法,使用分解濾波器逐層抽樣分解全色和多光譜數據,獲取高、低頻分解數據.

(5)對小波包分解后的多光譜數據和全色數據進行融合.采用基于區域特征的自適應小波包融合算法,通過直方圖均衡化,求分解窗口的方差、能量和信息熵,計算像素權值,按照公式獲取融合子數據.

設小波分解后的多光譜數據為A(x,y),全色數據為B(x,y),將多光譜數據進行色彩分離,得到三波段的子數據分別為(x,y)(k=1,2,3),j是尺度函數,將全色數據分別對多光譜數據的3 個波段進行直方圖均衡化.

其中,P(i)表示一像元i在數據中的出現概率,其范圍是[0,1,···,L].

其中,G(i,j)為像元灰度值,為均值.

其中,P[i][j]是特征提取算子,本文取經驗值.

其中,a,b,c為各特征的權值,默認為1.

按照式(10)獲取融合后數據(x,y):

(6)根據上述的小波包重構濾波器逆向逐層插值重構數據,獲取融合后數據.

3 數據介紹

3.1 多光譜LandSat TM 數據

LandSat是美國陸地探測衛星系統,TM (Thematic Mapper)是LandSat-4、LandSat-5 攜帶的專題繪圖儀[17],LandSat TM 波段信息如表1所示.本文采用Level 1T標準地形校正產品,選用Band 4、Band 3、Band 2 進行融合,有利于分辨植被、土壤和道路等信息[18].

表1 LandSat TM 波段信息

3.2 全色SPOT-5 數據

SPOT 衛星是法國空間研究中心(CNES)研制的地球觀測衛星系統,SPOT-5 星上載有2 臺高分辨率幾何成像裝置(HRG)、1 臺高分辨率立體成像裝置(HRS)、1 臺寬視域植被探測儀(VEG)等,前后模式實時獲得立體像對[19],其波段信息如表2所示.本文選用SPOT-5 星的全色波段PAN的10 m 分辨率HRS 數據.

表2 SPOT-5 衛星的波段信息

3.3 合成孔徑雷達SAR 數據

合成孔徑雷達SAR (Synthetic Aperture Radar)采用搭載在飛機或衛星上的移動雷達,達到和大型天線同樣精度的雷達系統.特點是分辨率高,能全天候工作,能有效地識別偽裝和穿透掩蓋物[20].本文選用歐洲空間局的ERS-2 (the second European Remote sensing Satellite)雷達數據,平均軌道高度為780 km,其參數如表3所示[20].

4 結果分析

本文采用HSV-WT、AWP 算法對兩組數據進行處理,選取的濾波器窗口為3×3,小波變換法采用了3 層小波分解和重構算法,Haar 小波基.

4.1 多光譜LandSat TM 數據和全色SPOT 數據融合

采用HSV-WT和AWP 對兩組數據進行融合.

4.1.1 讀取源數據

先選擇TM的Band 3 與SPOT 進行配準,可發現兩圖左邊部分的山脈區域清晰可見,圖5(b)中有一塊白色云區.

表3 合成孔徑雷達成像模式特性參數

圖5 源數據

4.1.2 數據配準校正

將圖5(a)作為基準數據,對圖5(b)數據進行校正.通過控制點位法進行相對配準,再通過二次多項式糾正法進行絕對配準,如圖6所示.同理校正另外兩個波段數據.

采用多項式糾正法時,誤差如表4.由表4所知,采用一次多項式糾正時誤差較大,因為一次項糾正僅對旋轉、平移、縮放帶來的誤差進行了糾正,并沒有對非線性變化引起的誤差進行改正.采用二次多項式時,中誤差明顯變小,不僅糾正了線性變形,還進一步改正了非線性變形.

4.1.3 融合TM Band 4、Band 3、Band 2

采用二進制小波變換法融合校正后的TM 三通道數據,如圖7所示.

圖6 校正后的TM Band 3

表4 多項式糾正誤差(單位:m)

圖7 TM Band 4、Band 3、Band 2 融合結果

分析融合結果,左邊區域為山脈,左邊有一塊白色區域為云,旁邊與其形狀相似的黑色區域為云影.衛星觀測的角度以及山脈的凸凹不平等因素導致云和云影之間的位置差異;深色區域為水體.

4.1.4 融合TM和SPOT 數據

分別采用HSV-WT、AWP 對圖5(a)和圖7進行融合,如圖8所示.

4.1.5 融合評價

觀察圖8發現,使用兩種方法融合后的結果既融入了全色SPOT 數據的細節信息,又保留了TM 數據的光譜信息.融合后數據在保持光譜信息和增強空間信息兩方面得到提高.對比圖8(a)、圖8(b)發現,圖8(a)的光譜失真較大,視覺效果不如圖8(b).在圖8(b)中,山脈的紋理特性和右邊部分的細節特征更加豐富,更有利于地物判別.比較兩種融合結果的性能參數如表5.

圖8 融合結果

表5 TM和SPOT 融合后數據的參數結果

從表5可看出,融合結果圖8(a)、圖8(b)在多光譜TM 數據的基礎上,提高了方差,其信息熵和清晰度都比源數據大;對比融合結果圖8(a)、圖8(b),在方差、信息熵、清晰度和相關指數4 項指標上,圖8(b)較大,且圖8(b)的扭曲程度、偏差指數較小.說明采用AWP 方法較HSV-WT 法更優,采用AWP 法處理過的圖像更加清晰,紋理信息更豐富,對象的幾何特征更完整.

4.1.6 融合另一組倫敦地區數據

對另一組倫敦地區的數據采用相同方法進行處理,選用TM的R、G、B 波段合成彩色影像.

圖9(a)采用TM的R、G、B 波段融合表示彩色圖像,圖像色彩更接近自然色彩,更符合人的視覺特性.圖中,黑色區域為水體,綠色區域為地表.放大觀察圖9(c)、圖9(d),圖像的細節信息得到明顯改善,光譜信息的增加提高地物的紋理特性.目視觀察比較圖9(c)、圖9(d),圖9(d)的色彩更為自然,紋理細節更豐富,目視效果更加清楚.表6列出圖9融合結果的性能參數,

綜合比較,圖9(d)數據優于圖9(c),說明采用AWT效果較好.

圖9 倫敦地區的數據融合結果

表6 圖9融合后數據的參數結果

4.2 多光譜LandSat TM 數據和SAR 數據融合

對意大利羅馬地區的TM和ERS-2 SAR 數據進行融合,使用AWP 實現.

4.2.1 讀取源數據

讀取羅馬地區的源數據,將TM Band 4、Band 3、Band 2 進行二進制小波融合,如圖10所示,SAR 數據有一定的立體視覺性,圖10(b)中黑色區域為水體.

圖10 意大利羅馬地區的源數據

4.2.2 數據配準

將圖10(a)的SAR 數據作為基準數據,對圖10(b)進行校正.校正法同上,圖11所示為校正后的TM 數據.

4.2.3 數據融合

由4.1.5 節和4.1.6 節知,基于自適應小波包融合算法較優,故采用此算法融合圖10的SAR 數據以及圖11的校正后TM 數據,融合結果如圖12.

圖11 校正后的TM 數據

4.2.4 融合評價

由圖12,融合后影像具有SAR 數據的立體視覺性,在融入TM的光譜信息后,增加了紋理的清晰度,有助于進一步區分地表的覆蓋類型.表7列出圖12融合結果的性能參數,

由表7知,融合后的方差、信息熵均高于源數據,說明融合豐富了信息量.融合后的清晰度提高較為明顯,說明融合后數據增強了細節特征.

圖12 SAR和TM 融合結果

表7 圖12融合后數據的參數結果

4.3 與其他融合算法的對比分析

基于4.1.1 節的多光譜LandSat TM 數據和全色SPOT數據,分別采用常用融合算法HSV 變換、主成分變換(PCA)、多進制小波變換(WT)對同一組數據進行融合,與采用本文算法HSV-WT、AWP的融合結果進行比較,如表8所示.

從表8可看出,融合結果均在TM 數據的基礎上,提高了方差,其信息熵和清晰度都比源數據大;對比各項參數,本文算法AWP、HSV-WT 在各項指標上,均優于其他3 種傳統算法,驗證了本文改進算法的有效性.

表8 不同融合算法的參數結果

5 結論

本文將低分辨率數據的多光譜信息、全色數據的高分辨率信息或SAR 數據的高空間信息有機結合起來,關注在盡可能保持原光譜信息的同時,提高融合結果的空間分辨率.主要結論如下:

(1)兩種改進算法均優于傳統融合算法,融合后的數據在保持光譜信息和提高空間細節信息兩方面均得到提高.這兩種算法相對傳統小波算法,能克服對高頻信息處理的缺陷,突破待融合數據的分辨率比值限制,實現分辨率之比非2n的數據融合.

(2)HSV-WT 將HSV的空間信息表達能力與多進制小波變換的良好局部化性質有機結合,從而使融合后的數據在保持光譜信息和提高空間分辨率兩方面的綜合性能達到平衡.實驗證明此方法簡單有效,融合結果既融入了SPOT 數據的細節信息,又保留了TM 數據的光譜信息.

(3)AWP 克服小波變換容易受分解階數影響的缺陷,對高頻部分進一步劃分,通過遞歸分解突出高分辨率數據的細節層次信息.用此方法處理兩組數據,效果優于HSV-WT.

(4)下階段應進一步探討針對要突出的不同特征,怎樣選取波段組合;SAR 數據具有穿透云層、植被的能力,它所獲得的信息取決于物體的幾何特性和介電特性,可思考怎樣發揮其優點,最大限度地提取地物特征信息.