楊勇,盧航遠,黃淑英,涂偉,李露奕
(1.江西財經大學 信息管理學院,南昌330032; 2.江西財經大學 軟件與物聯網工程學院,南昌330032)
遙感的目的是通過獲取衛星的光譜測量,提取有關地球表面結構和內容的信息[1]。由于傳感器技術的限制,同時在空間域和光譜域上獲取高分辨率難以實現,因此通常利用圖像處理技術融合現有的2類遙感圖像,即高空間分辨率、低光譜分辨率的全色(Panchromatic,PAN)圖像和低空間分辨率、高光譜分辨率的多光譜(Multispectral,MS)圖像。對于在同一區域同時獲得的PAN圖像和MS圖像,通過適當的算法,將這些數據結合起來,生成具有高空間分辨率和高光譜分辨率的圖像,這一過程也稱為PAN銳化[2]。遙感圖像融合廣泛應用于地理、軍事、農業、生態環境等領域。
現有的傳統PAN銳化方法主要分為2類:一類是成分替代法(Component Substitute,CS),即將MS圖像的成分替換成PAN圖像成分,如亮度-色度-飽 和 度 變 換 法 (Intensity-Hue-Saturation,IHS)[3]、主成分分析法(Principle Component Analysis,PCA)[4]、基 于 施 密 特 正 交 化 方 法(Gram-Schmidt,GS)[5]等。基于CS方法的融合結果具有較高的空間分辨率,但其光譜容易失真。另一類是多尺度分析法(Multi-Resolution Analysis,MRA),即對PAN圖像進行多尺度分解后將提取的細節加入到MS圖像中。基于MRA的方法和基于CS的方法的主要區別在于如何提取空間細節[6]。對于MRA方法,通過PAN圖像與其低通濾波之后的圖像之間的差異來獲得空間結構信息。傳統的MRA方法主要包括基于金字塔分解或者小波變換的方法,如離散小波變換[7]、àtrous小波變換[8]、非抽樣的輪廓波變換[9]、非抽樣的剪切波變換[10]等。MRA方法可以較好地保存光譜信息,但可能會損失部分空間信息,因此文獻[11]提出了一種基于引導濾波多尺度分解的方法,效果有所改進。最新的理論研究提出了稀疏表示(Sparse Representation,SR)的方法[12-13],以及非線性的途徑如卷積神經網絡[14]等,取得了較好的效果。但是,這些算法需要大量的計算,運行時間較長,且需要大量的樣本。由于遙感圖像樣本數量的限制,以及算法效率的要求,本文綜合考慮了上述幾種方法,對基于MRA方法的注入模型進行優化,以實現融合時對光譜信息與細節信息的雙保真。
基于MRA方法中關鍵的一步在于如何提取PAN圖像細節。一種較為先進的方法是利用線性時不變的濾波器去匹配MS傳感器的點傳播函數(Point Spread Function,PSF)。Aiazzi等[15]提出使用高斯濾波器去匹配MS傳感器的調制傳輸函數,但需要使用MS傳感器的出廠信息,而該信息不易獲得。另一種較為可靠的方法是從現有圖像中估計濾波器以模擬MS傳感器的PSF。在此框架下,一幅低空間分辨率的MS(Low Resolution MS,LRMS)圖像,可以通過濾波器對一幅高空間分辨率的MS(High Resolution MS,HRMS)圖像進行空間濾波得到。該濾波器的沖擊響應函數可模擬MS傳感器的PSF,且該濾波器通常有與高斯模型相似的形狀[2,15]。使用該濾波器卷積PAN圖像,所得的高頻分量即為HRMS圖像中缺失的細節成分,且與MS圖像高線性相關。此外,綜合考慮光譜的保真和細節的注入,本文提出了一種基于相關性分析的高斯濾波估計和自適應細節注入的遙感圖像融合方法。首先,通過引導濾波多尺度分解的注入模型得到初始融合圖像;其次,用高斯濾波模擬MS傳感器的PSF,并用估計的高斯濾波器卷積PAN圖像得到優化的注入細節;然后,遞歸計算自適應的注入系數,使光譜保真與細節注入達到聯合最佳;最后,通過一種新的邊緣保持策略,得到融合圖像。
本文主要創新點如下:①提出了一種基于高斯濾波估計的細節提取方法,使注入的細節得到優化;②提出了一種綜合考慮光譜信息與細節信息的自適應注入系數模型,實現融合圖像對光譜信息與細節信息的雙保真;③提出了一種減少光譜信息損失的邊緣保持模型,在增強邊緣信息的同時保留光譜信息。
基于MRA方法的注入模型是指通過濾波提取PAN圖像細節,再利用注入系數將其注入到MS圖像中的一種融合方法。該方法可根據實際需要,組合不同的融合技術,以發揮不同融合技術的優點。基于傳統的MRA方法注入模型框架如圖1所示,該注入模型的一般表示形式如下[16]:


圖1 注入模型一般框架Fig.1 Basic framework of injection model

對于注入系數Gk,文獻[16]通過迭代算法,提出基于回歸模型的通道間權重系數;文獻[17]提出基于像素通道間比例的權重系數;文獻[18]通過高通調制(High Pass Modulation,HPM)方法,用MS圖像和PAN圖像的乘法組合作為注入系數。其中,基于通道間比例的注入系數能更好地保持光譜信息,融合效果較好,通道間注入系數定義如下:

由于注入細節的同時不希望改變色度與飽和度,因此圖像融合過程中往往僅對圖像的亮度分量進行處理。MS圖像的亮度分量I可由如下公式計算:

式中:αk為通道的組合系數。
文獻[19-20]提出I分量可自適應地由各通道線性組合表示,即通過求解以下優化問題得到組合系數:

式中:P為PAN圖像。
在得到MS圖像的I分量之后,將P與式(3)的I進行直方圖匹配,得到的PI可由如下公式計算:

式中:mean(P)為P的均值;mean(I)為I分量的均值;std(I)為I分量的標準差;std(P)為P的標準差。
引導濾波最早由文獻[21]提出,其可以保存輸入圖像的主要信息,同時獲得引導圖像的變化趨勢。將PI作為輸入圖像,I作為引導圖像,PI通過I引導濾波產生PL,獲得了與MS圖像線性相關的細節信息。引導濾波可由如下簡化公式而得:

式中:GF(·)表示引導濾波函數。
多尺度引導濾波[22]可表示為

目前,遙感圖像融合注入模型主要存在2類問題:因PAN圖像與MS圖像相關性不高導致的光譜失真;細節注入過多或者注入不足。針對上述問題,本文提出了一種基于高斯濾波估計的細節提取及自適應的注入系數優化方法。本文融合方法框架如圖2所示。
首先,根據式(1),將上采樣后的MS圖像與PAN圖像通過直方圖匹配、引導濾波等傳統注入模型得到初始融合圖像F(1);其次,進行細節優化,即用多尺度高斯濾波器濾波該初始融合圖像以模擬MS傳感器的特性,將得到的高斯濾波估計濾波PAN圖像,得到優化后的細節;然后,優化注入效益,即綜合考慮光譜信息與細節信息計算自適應的注入量系數,具體流程如圖2上虛線框所示,詳細算法見2.2節;最后,通過一種自適應的邊緣保持權重矩陣保持邊緣信息,具體流程如圖2下虛線框所示,詳細算法見2.3節。將優化后的細節與注入系數及邊緣權重系數相乘以后,注入到上采樣后的MS圖像中,得到最終的融合圖像。
文獻[15]提出使用高斯濾波器匹配MS傳感器的調制傳輸函數,即HRMS圖像經過MS傳感器之后成為LRMS圖像,使用高斯濾波模擬該過程,所得的高斯濾波器即為符合MS傳感器PSF的濾波估計。使用估計的濾波器卷積PAN圖像,提取的細節成分即為低分辨率MS圖像所缺失的細節。文獻[15]通過實驗證明了該濾波估計有較強的魯棒性,即使在初始融合圖像中加入白噪聲,也不影響細節的提取。另外,文獻[17]提出了二步驟多尺度分解的方法來精細化PAN銳化。受此啟發,本文提出使用多尺度高斯濾波模擬MS傳感器的PSF,即對初始融合圖像F(1)提取亮度分量,記為I1,并迭代地進行高斯濾波,具體可表示為


圖2 本文融合方法框架Fig.2 Framework of proposed fusion method
為了獲取最佳高斯濾波估計,計算式(8)所得的結果Ii1與上采樣后的MS圖像的亮度分量I的相關系數,相關系數最高時,估計的高斯濾波最準確。相關系數計算如下:

式中:cov(·)表示協方差函數。
隨著濾波層次的增加,濾波后的圖像與上采樣的MS圖像越來越接近,達到極大值之后,隨著濾波層次繼續增加,濾波后圖像與上采樣的MS圖像偏離越來越大。通過迭代計算式(8)和式(9)的值,取得式(9)最大時的迭代次數m,即為所估計的高斯濾波器Hm,表示如下:

Hm即表示一幅銳化圖像經MS傳感器模糊化的PSF最佳匹配,將濾波器Hm作用于PAN圖像,所提取的細節成分也應符合MS傳感器PSF,這樣就減少了光譜失真。細節提取D定義如下:

結合式(1)、式(7)和式(11),細節優化后,第k通道融合圖像可更新為

取得優化的注入細節之后,對注入量進行優化。不同的圖像有不同的結構特征和光譜信息,從而對細節注入量有不同的要求,直接等比例注入細節容易導致注入過多從而引起光譜失真,或者注入不足導致圖像模糊,因此需要通過一個注入量系數g來平衡光譜信息與細節信息。加入該系數后,結合式(12),新的融合圖像可以定義為

需要對不同結構特征的遙感圖像確定不同的注入量系數g。文獻[23]認為注入模型中光譜的保真與空間分辨率要求需要保持平衡,可以將兩者的評價指標進行加權求和作為融合圖像指標,權重表示重視程度。設定細節注入量系數g的初值g和步長r,對于g0≤g≤1,得到相應的融合后圖像并計算加權評價指標,取得使該評價指標最高的gmax作為最終的注入量系數。對于評價指標,本文提出了一種基于相關性分析的光譜信息與空間信息綜合評價體系,將光譜信息與空間信息的綜合評價指標Q定義為

式中:ESP和EHF分別為光譜信息評價指標和空間信息評價指標;α為兩者的權重。
對于一幅由式(13)融合后的圖像F(3)
k,其亮度分量記為I3,則ESP用融合后的圖像與上采樣的MS圖像各通道k之間相關系數的平均值來表示,EHF則用I3與PAN圖像的線性相關性表示,定義為

文獻[6]表明,在注入成分與注入量系數g的優化確定上,遵循的前提是MS圖像的亮度成分與PAN圖像相關性越大,融合效果越好。即若兩者相關性高,注入細節的同時能較好地保存光譜信息,此時可以適當增加空間信息的權重,減少光譜信息的權重;反之,若相關性低,則注入的細節容易引起光譜失真,應當增加光譜信息權重而減少空間信息權重。設初值g0產生的融合圖像對應的亮度成分記為I0,則式(14)中的權重α可由下式確定:

平方項的作用為:①避免I0與PI高相關時光譜權重過小;②放大不同圖像空間信息之間的區別。

結合式(13)和式(18),使用優化后的注入量系數,新的融合圖像可更新為

為保持邊緣紋理信息,Leung等[20]提出利用邊緣信息作為各通道融合權重,Yang等[11,24]提出利用改進的邊緣信息作為注入各像素的權重。其中,改進的邊緣信息計算PAN圖像的邊緣檢測矩陣WP,以及上采樣后的MS圖像各通道的邊緣檢測矩陣W,并將兩者進行加權和作為融合圖像的邊緣信息權重。該方法在邊緣信息保護上起到了較好的效果。第k通道的邊緣保持權重矩陣Ek定義如下[20]:

式中:βk為各通道的加權系數,可由式(21)計算得到[11];邊緣檢測矩陣[19]可由式(22)得到。

式中:A為某圖像;ΔA為圖像A的梯度;λ和ε為調制參數。
式(21)獲取各通道的邊緣信息權重,在細節上有更好的表現,但作為注入系數時改變了各通道間的比例,光譜信息略有損失。為了盡量保持原有光譜信息,同時實現邊緣信息的保持,本文提出一種新的邊緣保持權重矩陣。針對式(21),將求取各通道的邊緣檢測矩陣Wk替換為求取MS圖像亮度成分I的邊緣檢測矩陣WI,并計算PAN圖像的邊緣檢測矩陣WP與WI加權和,以綜合考慮PAN圖像和MS圖像的邊緣信息。邊緣保持權重矩陣E可表示為

式中:WP和WI計算方法同式(22);自適應的權重因子βnew由式(24)計算,并由式(25)歸一化得到。

根據式(23)得到的邊緣保持權重矩陣E,再結合式(19),可得到最終融合圖像,表示如下:

與式(20)的邊緣保持權重矩陣相比,本文得到的邊緣保持權重矩陣E僅對細節成分D產生權重變化,以盡可能增強邊緣信息,并不改變通道間的比例,從而減少光譜損失。
為評估本文方法的有效性,使用來自不同遙感衛星的四大數據集,包括Quick Bird、IKONOS、WorldView-2及pleiades,這些遙感圖像在空間分辨率、光譜波長等方面都有不同的特征,且都包括紅、綠、藍、近紅外4個通道。
參數設置中,設置注入量系數初值g0=0.1,步長r=0.05,以綜合考慮計算效率與計算精度;高斯濾波器窗口設置為5×5,參考文獻[25]中的默認設置;式(22)中設置參數λ=10-9,ε=10-10,同樣參考文獻[19]中的默認設置。
本文進行了2種類型的實驗。一種是對真實圖像進行下采樣之后的仿真圖像,也稱為有參考遙感圖像融合實驗,即遵循Wald協議[26],原圖像作為參考圖像,仿真圖像包含四大數據集共計80組圖片。另一種是真實圖像的實驗,也稱為無參考遙感圖像融合實驗,真實圖像包括四大數據集共計50組圖片。
用于對比的方法包括最新的幾種方法,其中基于MRA模型的方法有基于調制傳遞函數模型的多尺度分析方法(CBD)[15]、加性小波亮度比的方 法(Additive Wavelet Luminance Proportion,AWLP)[27]、基于雙邊濾波亮度比的方法(Bilateral Filter Luminance Proportional,BFLP)[28];將MRA和CS兩種模型相結合的方法有基于多尺度引導濾波和改進的成分替代法(Improved IHS and Multi-Scale Guided Filter,IMG)[11];將注入模型與其他算法相結合的方法有基于摳圖模型和多尺度變換的方法(Matting Model and Multi-scale Transform,MMMT)[24]、基于自適應光譜與亮度調制的方法(Adaptive Spectral Intensity Modulation Pansharpening,ASIMP)[25]。對這些方法分別進行主觀評價與客觀評價,也即定性和定量評價。
評價包括肉眼視覺的主觀評價和客觀評價指標的定量評價。常用的有參考圖像的客觀評價指標如下:
1)相關系數(Correlation Coefficient,CC)[12]。主要用于評價融合圖像與參考圖像的空間相似度,其值區間為[0,1],越接近1表示2幅圖像越相近,融合效果越好。
2)通用圖像質量指數(Universal Image Quality Indexes,UIQI)[29]。UIQI聯合相關度損耗、亮度失真和對比度失真3個因子來評價融合圖像的效果,值區間為[0,1]。UIQI的值越接近1,說明結構失真越小。
3)均方根誤差(Root Mean Square Error,RMSE)[30]。RMSE計算融合圖像與參考圖像各像素之間均方根誤差,用于表示兩者的差異,其理想值為0。
4)相對平均光譜誤差(Relative Average Spectral Error,RASE)[31]。RASE反 映 融 合 圖 像在光譜方面的誤差,其理想值為0。
5) 光 譜 角 映 射(Spectral Angle Mapper,SAM)[32]。SAM用于度量2幅圖像之間的光譜信息接近程度,其值越小,光譜失真越小,理想值為0。
6)全 局 相 對 光 譜 損 失(ERGAS)[33]。ERGAS用于評估融合圖像的空間和光譜質量,其值越小說明總體性能越好,理想值為0。
常用的無參考圖像的客觀評價指標包括Ds、Dλ和QNR[34]。其中,QNR值由Ds和Dλ共同計算而得,并基于UIQI測度評價融合圖像和參考圖像間的局部相關性、亮度和對比度。Ds和Dλ越小越好,理想值為0,而QNR值越大越好,理想值為1。
為了驗證2.1節~2.3節提出的優化方法的有效性,定量說明本文方法所帶來的性能提升,分別計算各優化步驟的平均客觀評價指標。4個優化步驟對應的融合圖像表示為:F1表示初始融合圖像;F2表示經高斯濾波優化細節后圖像;F3表示經自適應優化注入量系數后的圖像;F4表示加入自適應邊緣保持后的圖像。并分別計算80組圖像融合結果的平均客觀評價指標。這4個步驟是依次遞進的過程。為了更直觀地展示各指標的相對變化,各客觀評價指標分別進行最大值歸一化處理,結果如圖3所示。可知,對于初始融合圖像F1,注入細節是通過PAN圖像與MS圖像的I分量進行直方圖匹配再由引導濾波得到,細節存在與MS圖像線性非相關的問題,且未考慮注入量,從而導致融合結果與參考圖像的光譜相關度和空間相關度較低;通過估計的高斯濾波器優化細節后,所得的細節與MS圖像有更高的相關性,各項指標明顯提升;而優化注入量系數后,有效防止了注入過多或注入不足,所有的指標進一步提升;進行邊緣保護后,改進的邊緣保持方法在保護邊緣信息的同時,減少了光譜失真。隨著F2、F3、F4優化步驟的依次進行,總體指標也越來越優。其中SAM指標基本保持平穩,隨著優化的進行,SAM指標先小幅上升再小幅降低,說明本文方法對光譜信息損失的影響很小。

圖3 優化過程平均客觀評價指標對比Fig.3 Average objective evaluation indicator comparison of optimization steps
如2.3節所述,本文提出的一種新的自適應邊緣保持方法是對文獻[11]的改進。為了驗證本文方法的有效性,在保持前3個步驟不變的前提下,對改進前[11]與改進后的方法在80組遙感圖像上進行平均客觀評價指標對比,指標同樣進行最大值歸一化處理,結果如圖4所示。可知,改進后的邊緣保持方法在所有指標上都有所提升,其中,在RASE、RMSE及ERGAS等指標上提升效果明顯,可知本文方法能有效地提升圖像融合的質量,減少光譜損失。

圖4 邊緣保持方法改進前后平均客觀評價指標對比Fig.4 Comparison of average objective evaluation indicators before and after improvement of edge preserving method
本文采用3組來自不同衛星的遙感圖像作為展示,分別進行主觀效果和客觀評價指標的評價,然后對80組圖像測量平均客觀評價指標。
第1組圖像來自QuickBird數據集,參考圖像大小為256×256,經不同方法融合后主觀效果如圖5所示。小紅色矩形框區域經放大后在大矩形框中顯示。可知,AWLP方法、CBD方法、IMG方法森林區域存在明顯光譜失真,與參考圖像相比,顏色偏亮;BFLP方法、MMMT方法有稍微的光譜失真,且存在過度注入問題,引入了不必要的噪聲。ASIMP方法與本文方法在海灘細節與森林的光譜信息上與參考圖像較為接近,難以用肉眼區別,用表1所示的客觀評價指標加以區分。表中:黑體數據表示最優值,下劃線數據表示次優值。可知,本文方法除了在SAM 指標上比BFLP方法略低以外,在所有其他指標上都達到最優。

圖5 QuickBird數據集遙感圖像融合結果Fig.5 Fusion results of remote sensing images from QuickBird dataset
第2組圖像來自World View-2數據集,參考圖像大小為256×256,經不同方法融合后主觀效果如圖6所示,客觀評價指標如表2所示。由圖中紅色放大矩形框可知,IMG方法有較為明顯的光譜真失,亮度偏亮;AWLP方法、BFLP方法和CBD方法有輕微的光譜失真,表現為屋頂區域顏色更淺,且在屋頂、中間白色地面有較為明顯的噪聲;MMMT方法、ASIMP方法與本文方法能較好地保持光譜信息,但MMMT方法在邊緣上更加模糊,而本文方法與參考圖像不論在空間信息上還是光譜信息上都更為接近。為更好地區分融合效果,可通過表2查看融合結果客觀評價指標,可知,本文方法在所有指標上都達到最優。
第3組圖像來自pleiades數據集,經不同方法融合后主觀效果如圖7所示,客觀評價指標如表3所示。從放大的紅色矩形區域可知,AWLP方法、BFLP方法和CBD方法存在較為明顯的光譜失真,且引入了較多的噪聲;IMG方法、MMMT方法和ASIMP方法有較輕微的光譜失真,引入少量的噪聲;而本文方法能在保持空間細節的同時,較好地保留光譜信息,色彩更加真實。對于難以用肉眼明顯區分的方法,可通過表3查看融合結果客觀評價指標,可知,本文方法除了在UIQI指標上比ASIMP略低以外,在所有其他指標上都達到最優。

表1 圖5對應的遙感圖像融合結果客觀評價指標Table 1 Objective evaluation indicators of remote sensing image fusion results in Fig.5

圖6 WorldView-2數據集遙感圖像融合結果Fig.6 Fusion results of remote sensing images from WorldView-2 dataset

表2 圖6對應的遙感圖像融合結果客觀評價指標Table 2 Objective evaluation indicators of remote sensing image fusion results in Fig.6
為了更客觀地評價各方法的性能,對來自QuickBird、IKONOS、WorldView-2及pleiades四大數據集的共計80組有參考遙感圖像融合后計算平均客觀評價指標,結果如表4所示。可知,本文方法在所有的平均客觀評價指標上達到了最優。

圖7 pleiades數據集遙感圖像融合結果Fig.7 Fusion results of remote sensing images from pleiades dataset

表3 圖7對應的遙感圖像融合結果客觀評價指標Table 3 Objective evaluation indicators of remote sensing image fusion results in Fig.7

表4 80組遙感圖像仿真實驗的平均客觀評價指標Table 4 Average objective evaluation indicators of simulation experiments from 80 groups of remote sensing images
真實圖像的實驗采用2組圖像進行主觀和客觀評價,并對所有來自四大數據集的50組真實圖像計算平均客觀評價指標。
第1組圖像來自IKONOS數據集,圖像大小為780×608。不同方法融合后的主觀效果如圖8所示,相應的客觀評價指標如表5所示。由圖8可知,BFLP方法和MMMT方法在橙色屋頂區域有明顯光譜失真,其他方法與本文方法難以用肉眼直接判別,可通過表5中圖8的客觀評價指標來評價,可知,本文方法除了在Ds指標上略低于ASIMP方法以外,其他指標都達到最優。
第2組圖像來自QuickBird數據集,圖像大小為512×512,融合后視覺效果如圖9所示。可知,AWLP方法、CBD方法和IMG方法亮度偏亮,存在過度注入的問題,其他方法難以用肉眼直接觀察,可由表6中圖9的客觀評價指標測量,可知,本文方法在所有客觀評價指標上都達到最優。
為更加客觀評價不同方法的有效性,本文對來自四大數據集的50組真實圖像進行實驗,并測算平均客觀評價指標,如表7所示。可知,本文方法所有指標都達到最優。

表5 圖8對應的真實圖像定量評價結果Table 5 Quantitative evaluation results of real image in Fig.8

圖9 QuickBird數據集真實圖像融合結果Fig.9 Fusion results of real images from Quick Bird dataset

表6 圖9對應的真實圖像定量評價結果Table 6 Quantitative evaluation r esults of real image in Fig.9

表7 50組真實圖像的平均定量評價結果Table 7 Average quantitative evaluation results of 50 groups of real images
1)本文方法在優化注入細節的過程中使用高斯濾波估計去匹配MS傳感器的特性,并用該估計的高斯濾波器去卷積PAN圖像,得到優化后的細節。
2)對于注入效益的優化,本文方法綜合考慮光譜信息和細節信息,并建立綜合評價指標,以確定自適應的注入量系數。
3)為加強邊緣細節,本文方法提出一種光譜信息保真的自適應邊緣提取方法,以對注入效益進一步優化。
4)通過對來自QuickBird、IKONOS、World-View-2及pleiades四大數據集的80組仿真圖像與50組真實圖像進行融合實驗,并將本文方法與一些先進的融合方法進行對比,結果表明本文提出的方法在主觀和平均客觀評價指標上能達到最優。
在未來的工作中,可能需要精細化高斯濾波估計,以得到更準確的濾波器。這些工作將在以后的研究中進一步嘗試。