袁曉冬
(北京航空航天大學 計算機學院,北京 100191)
李 超 盛 浩
(北京航空航天大學 深圳研究院,深圳 518057)
高空間分辨率全色圖像反映空間結構信息,能詳盡地表達地物細節特征;低空間分辨率高光譜圖像的光譜信息豐富,有利于地物識別與解釋.利用全色圖像的空間細節信息和多光譜圖像光譜信息融合得到高空間分辨率的多光譜圖像,如何在提高多光譜圖像空間分辨率的同時,盡可能降低光譜失真度一直是多源遙感圖像融合的難點[1].
針對全色圖像與多光譜圖像的融合,主要有基于主分量分析的方法[2]、基于色彩變換的方法[3]、基于高通濾波的方法[3]等,這些方法由于沒有考慮到多源圖像高頻分量之間的相關性,形成的融合圖像都存在一定的光譜失真.
本文借助多分辨分析思想,利用小波分解具有的良好變換域分頻、方向特性[4],提出一種在圖像小波尺度空間根據頻率分量和局域相關性確定不同融合算子的策略,建立小波分解層融合圖像序列,最后通過逆變換重構得到最終的融合圖像.
多光譜圖像的光譜信息主要體現在低頻分量上,融合增強需要保留低頻分量作為其光譜信息的保持.事實上,低頻分量只是體現了平穩地物的光譜信息,高頻分量同樣也體現了非平穩地物的光譜信息.
Mallat提出的小波分解快速算法,使圖像小波分解的實際應用成為可能[5].圖像在尺度s-1上的小波分解如下:
其中,L和H分別為低、高通一維鏡像濾波算子,下標r和c為對應圖像的行和列;As為圖像s尺度上的近似低頻分量;Vs,Hs,Ds分別為垂直、水平和對角方向的高頻分量.小波逆變換的表示如下:

其中L*和H*分別為L和H的共軛轉置矩陣.
本文將多光譜圖像和全色圖像經小波分解后,對高頻分量的融合定義一定的約束條件,針對局域相關性強弱分別處理,融合的原理如圖1所示.
多光譜圖像的低頻分量代表了平滑地物的光譜特性,與經過配準后的全色圖像低頻分量相似,其融合較為簡單,采用整體圖像能量或空間頻率進行線性加權即可.

圖1 小波多分辨分量相關性圖像融合原理
在小波分解的尺度s上,對分解圖像整體近似分量As基于圖像局域相關系數[6]加權融合計算.假定兩圖像A和B,設局域窗口W的中心像素點坐標為(m,n),則點(m,n)的局域相關系數Co的表達式如下:

其中,A(x,y)和B(x,y)分別表示圖像A和B在(x,y)處的像素值;和分別表示圖像A和B局域分析窗口W的像素平均值.
s尺度上的融合圖像整體近似分量計算如下:

對于各方向上的高頻分量,在s尺度上,設多光譜圖像和全色圖像的小波分解高頻分量分別為Ms和Ps,采用線性融合模型,得到小波域期望融合多光譜圖像的高頻分量為Fs:



其中Σx為x的協方差矩陣:

和分別為多光譜圖像、全色圖像的方差;為多光譜和全色圖像的協方差.采用期望圖像局部方差最大作為融合系數的估計依據,融合系數β的選取應滿足:max(β)=βTΣxβ,并同樣定義約束條件1:

為不引入虛假的高頻信息光譜失真,對圖像高頻分量的融合進行相關性約束,使光譜圖像高頻分量融合前后的相關系數達到門限值T(≥0),即

化簡得約束條件2如下:

同時,由T(≥0)可知,融合系數β的符號應該保證融合分量圖像的對比度變化極性與原始光譜圖像極性相一致,因而可得到約束條件3如下:

其中sign是取符號函數.
由以上分析可知,小波高頻分量融合系數β的求解即轉化為一個約束極值問題,相應的K-T條件(Kuhn-Tucker最優化條件)為

對式(10)等效化簡求解,當極值點不在μ2(β)上,則k2=0,μ2(β)>0,式(10)可表示為

可以看出此時k1是Σx的最大特征值,β為k1對應的特征向量,因此可得分量融合系數的估計值為

其中,v為Σx最大特征值對應的特征向量;abs表示取絕度值;⊙表示矢量乘法.
當極值點在μ2(β)=0上時,實際上約束極值問題轉化為求μ1(β)=0與μ2(β)的交點,代入各自的表達式可得

由于Σ為對稱矩陣,可以進行如下分解:

其中,λ1和λ2為Σ的特征值;U為特征向量形成的單位正交矩陣.則分量融合系數的估計值滿足如下關系式:

實驗平臺選擇為:Intel(R)Core(TM)2Duo CPU E7500 2.93GHz,4.00GB內存;開發環境VS2008,VC++基于OpenCV.實驗中分別對3組衛星遙感多光譜圖像和高空間分辨率全色圖像進行驗證.對圖像進行幾何配準后(大小為1280×960),分別采用色彩空間變換法(HIS,Hue-Intensity-Saturation)[7]、主分量分析法(PCA,Principal Component Analysis)[2]、拉普拉斯塔式分解法(LAP,Laplace)[8]、小波直接替換法(WAV,Wavelet)、小波平均選大方法(WAM,Wavelet Average Max)[9]和本文提出的方法(小波分解層數為3層),選取高頻相關門限值T=0.5,對圖像進行融合.
第1組實驗:
ZY02C衛星所獲取的全色圖像P(分辨率5m,圖2a)和多光譜圖像M(分辨率20m,圖2b),如圖2所示.
第2組實驗:
ZY3衛星所獲取的全色圖像P(分辨率5m,圖3a)和多光譜圖像M(分辨率20m,圖3b),如圖3所示.
第3組實驗:
高分全色圖像P采用UTM zone14N全色灰度圖像(分辨率5m,圖4a);多光譜圖像M采用UTM zone14N多光譜圖像(分辨率20m,圖4b),如圖4所示.
幾種方法的融合圖像在一定程度上保留了原始的光譜信息,空間分辨率也有明顯提高,不同的是3組圖像中的c,d,e,f,g圖與原始多光譜圖像M相比都存在顏色失真.本文方法的融合圖像,在提高光譜圖像空間細節表現能力的同時,圖像顏色偏差不大,較好地保留了原多光譜圖像M的光譜信息.

圖2 第1組各方法融合結果圖像對比

圖3 第2組各方法融合結果圖像對比

圖4 第3組各方法融合結果圖像對比
對于實驗性能,本文采用融合前后多光譜圖像相關系數和光譜扭曲度等指標進行對比.多光譜圖像相關系數[10]定義如下:

其中,M(x,y)和F(x,y)為融合前后多光譜圖像值;μM和μF為融合前后多光譜圖像均值.
多光譜圖像光譜扭曲度[11]定義為(設圖像大小為m×n)

空間頻率FS反映了一幅圖像總體的活躍程度和清晰程度[12],可用于衡量融合圖像的空間細節信息改善情況.
行頻:

列頻:

對角頻:

空間頻率:

為了對融合多光譜圖像整體波段相關性和扭曲程度進行量化分析,經過對單波段評價指標綜合分析,提出了兩個多光譜融合圖像整體評價指標.
多光譜圖像整體維度波普相關度:

其中,ΣF和ΣM分別為融合圖像和原多光譜圖像的單波段方差;ΣFM為融合圖像和多光譜圖像的單波段協方差;n為多光譜圖像波段數.
多光譜圖像整體維度波普扭曲度:

其中,ΣF和ΣM分別為融合圖像和原多光譜圖像的單波段方差;和為圖像的單波段均值;n為多光譜圖像波段數.
以第3組實驗圖像為例,表1中列出了本次實驗中上述幾種方法各自的光譜波段相關系數R(F,M)、光譜扭曲度D(F,M)以及與原高空間分辨率圖像空間頻率的比值RSF等參數,各數值都作(0,1)區間的歸一化處理.本文的方法設置相關系數門限值T=0.5,在使融合前后圖像的空間分辨率得到顯著改善的基礎上,光譜圖像的相關系數和光譜扭曲度都明顯優于其他幾種方法.

表1 本文方法和其他方法的融合性能指標比較
實驗中,針對上述提出的相關性約束條件2,以初值0.1,然后逐次增加0.1,選取不同的相關性門限閾值T,觀察融合多光譜圖像各波段相關系數R(F,M)、扭曲度D(F,M)以及圖像的整體空間頻率FS的變化情況,得到實驗數值見表2.

表2 不同的門限值T下的融合性能指標
由表2可以看出,隨著T選取數值的增大,融合多光譜圖像各波段的R(F,M)和D(F,M)也相應地分別增大和減小,圖像的整體空間頻率FS隨著相關門限值T的增大呈逐漸減小的趨勢,需要針對不同的實際應用場景在二者之間適當平衡.根據本文提出的整體衡量多光譜圖像融合質量的兩個評價指標,驗證了不同的相關系數門限值T對融合圖像的整體維度波普相關度CTC、整體維度波普扭曲度DTS,以及與原高空間分辨率圖像之間的空間頻率比值RSF等的影響,如圖5所示.實驗中選取初始T值為0,逐次增加0.1.

圖5 相關系數、光譜扭曲度與門限值T的關系曲線
由圖5可以看出,融合多光譜圖像的CTC和DTS隨著相關門限值T的增大分別呈逐漸增大和減小的趨勢,其與原全色圖像的空間頻率比值RSF也相應減小.對于實驗圖像,同時兼顧光譜信息和空間細節的情況下,最佳取值為T=0.75.因此針對不同的應用場景,可以通過參數T的調整,在融合圖像的光譜信息和空間細節兩方面進行有效調節.
針對多光譜圖像與高空間分辨率全色圖像的融合,提出了基于小波分解的多分辨分量相關性圖像融合方法.在原始圖像有效幾何配準的前提下,先對原始多光譜圖像和全色圖像進行小波分解,再對分解得到的低頻分量按局域相關系數、高頻分量按擬定的相關性約束條件分別進行融合,最后通過小波逆變換重構得到最終的融合多光譜圖像.實驗中,將本文的方法與幾種典型方法就相關系數和光譜扭曲度等進行對比,表明本方法在光譜信息保持和空間分辨率提高方面的優勢.為了對融合多光譜圖像整體的波段相關性和扭曲程度進行量化分析,提出了多光譜圖像整體維度波普相關度、多光譜圖像整體維度波普扭曲度兩個多光譜融合圖像整體評價指標.最后,實驗分析了不同相關系數門限值T的選取對最終融合效果的影響,可以針對不同的應用場景,適當地平衡融合圖像在光譜信息和空間分辨率兩方面要求.
(References)
[1]Tu Teming,Cheng Wenchun,Chang Chienping,et al.Best tradeoff for high-resolution image fusion to preserve spatial details and minimize color distortion[J].IEEE Geosciences and Remote Sensing Letters,2007,4(2):302-306
[2]ShahV P,Younan N H,King R L.An efficient pan-sharpening method via a combined adaptive PCA approach and con-tourlets[J].Geosciences and Remote Sensing,IEEE,2008,46(5):1323-1335
[3]Metwalli M R,Nasr A H,Farag Allah O S,et al.Image fusion based on principal component analysis and high-pass filter[C]//Computer Engineering & System,ICCES.Piscataway,NJ:IEEE Computer Society,2009:63-70
[4]Demirel H,Ozcinar C,Anbarjafari G.Satellite image contrast enhancement using discrete wavelet transform and singular value decomposition[J].IEEE Geosciences and Remote Sensing Letters,2010,7(2):333-337
[5]Kim Yonghyun,Lee Changno,Han Dongyeob.Improved additive-wavelet image fusion[J].IEEE Geosciences and Remote Sensing Letters,2011,8(2):263-267
[6]Zhu Junjie,Fan Xiangtao,Ding Chibiao,et al.Afusion method of high-resolution images considering spectral distortions[J].Journal of Wuhan University,2006,31(10):858-861
[7]Yang Jun,Zhao Zhongming.Aremote sensing images fusion method based on HIS and brightness control[J].Research of Computer Application,2007,24(4):195-197
[8]Cao Guangzhen,Hou Peng,Jin Yaqiu,et al.Alaplacetower remote sensing images fusion method based ontexture characteristic[J].Application of Remote Sensing Technology,2007,22(5):628-632
[9]Yang Ya,Wang Zheng,Zhang Sulan,et al.Amulti-focusremote sensing images fusion method based on wavelet transform[J].Computer Technology and Development,2010,20(3):56-58
[10]Yang Senlin,Gao Jinghuai.Seismicquality factors estimation from spectral correlation[J].IEEE Geosciences and Remote Sensing Letters,2008,5(4):740-744
[11]Al-Wassai F A,Kalyankar N V,Al-Zaky A A.Spatial and spectral quality evaluation based on edges regions of satellite:imagefusion[C]//Advanced Computing & Communication Technologies(ACCT).Piscataway,NJ:IEEE Computer Society,2012:265-276
[12]Beak Du Hyun,Yoon Jin Woo,Shin Jae Sung,et al.Restoration of high spatial frequency at the image formed by stimulated Brillouin scattering with a prepulse[J].Applied Physics Letters,2008,93(23):231113