基于分量替換高分辨率遙感圖像融合方法的對比研究

李春華

（福建師范大學 地理科學學院，福州350007）

Aiazzi等［1］在對25a來遙感圖像融合的回顧中最新提出了遙感圖像融合可以分為兩類方法：分量替換法和多分辨率分析法。分量替換融合方法是其中一類重要的融合方法，該方法基于新的向量空間對原始的多光譜圖像進行光譜變化。分量替換融合思想最早出現在Carper等［2］提出的IHS融合方法中。該方法首先對多光譜數據（僅限3個波段）進行IHS變換，再用Pan波段數據替換I分量，最后將新的數據集進行IHS逆變換即可獲得融合圖像。通常在進行波段替換之前，將Pan波段進行直方圖匹配，使得Pan波段總體的均值與方差和I分量一致。隨后出現的主成分融合法（PCA）和IHS融合法相似，區別在于PCA融合法是利用Pan波段數據替換第一主成分（PC1），并可以對多個波段進行融合處理［3］。

由 Laben 和 Brover［4］提 出 的 Gram－Schmidt（GS）融合算法是另一個基于分量替換思想的融合算法。該方法首先模擬低分辨率全色波段，將模擬的低分辨率全色波段作為第一分量參與GS正交變換，再將Pan波段替換第一分量構成新的數據集進行GS逆變換得到融合圖像。在進行波段替換之前，將Pan波段進行直方圖匹配，使得Pan波段總體的均值與方差和模擬的低分辨率全色波段一致。

以上提到的IHS，PCA和GS分量替換法的基本方法都是在進行逆變換之前，將Pan波段替換經變換后的第一分量（稱為可替換分量），在替換前首先對Pan波段進行直方圖匹配，使得Pan波段總體的均值、方差和所替換的分量一致，但是由于局部均值和方差的不一致往往導致融合圖像出現光譜失真。為了改善光譜失真，Aiazzi等［5］提出了改進的分量替換法，該方法首先通過對全色波段和多光譜波段進行線性回歸獲取擬合系數來構造可替換分量，再將Pan波段替代可替換分量進行逆變換獲取融合圖像，改進的分量替換法可以應用在IHS和GS融合算法中。分量替換法由于其簡單快速的特點，在高分辨率遙感圖像融合中一直得以運用，但涉及該方法關鍵性原理的闡述在國內還鮮有報道，在眾多的應用中往往由于缺乏對該融合算法機理的理解而未能得以正確使用。

為此，本文首先從線性代數的角度闡述分量替換融合算法的實質，并對基于光譜響應函數模擬低分辨率全色波段和基于多元線性回歸模擬低分辨率全色波段的這兩種GS融合算法進行原理分析。其次，以QuickBird數據為例，進行三種基于分量替換的融合方法的對比分析，分別是PCA融合法、基于光譜響應函數模擬低分辨率全色波段的GS融合法［6］、基于多元一次線性回歸擬合低分辨率全色波段的GS融合法［5］。最后，對上述三種融合算法的光譜保真性進行評價，并分析融合影像光譜失真的原因，以此來探明高分辨率遙感圖像融合的研究方向。

1 原理分析

1.1 分量替換融合算法的實質

Dou等［7］分析了分量替換融合方法的理論框架，本研究從線性代數的角度結合Dou的理論框架更為具體地闡述分量替換融合算法的技術細節和實質。

首先，假設原多光譜數據有n個波段，將多光譜數據用矩陣形式表示為：

式中：MSi（i＝1，2，…，n）——用行向量的形式表示的原始多光譜波段數據。

將多光譜波段數據進行線性變換（只須可逆變換即可）得到新的數據集，表示為：

式中：MSi′（i＝1，2，…，n）——經過線性變換后的各多光譜波段數據；Q——線性變換對應的可逆矩陣。

用全色波段Pan替換MS1′得到新的數據集，對新數據集進行逆變換，得到融合影像，見式（3）：

式中：MSfi（i＝1，2，…，n）——融合后的各多光譜波段數據；qi1（i＝1，2，…，n）——Q－1的第一列向量的各分量元素。以上的推理過程說明，任意的非退化矩陣（可逆矩陣）都可以定義相應的線性變換，式（3）表明融合的結果實質上是由兩部分構成，一部分是原始的多光譜圖像，另一部分Pan－MS1′代表Pan波段數據的高頻部分。融合的效果取決于線性變換的第一分量MS1′，MS1′稱為被替換波段，在IHS融合算法中MS1′即為I分量，PCA融合算法中MS1′表示第一主成份PC1，在GS融合算法中 MS1′稱為模擬的低分辨率全色波段。

因此，減少光譜失真的關鍵性技術就在于構造被替換的波段 MS1′。只有當 MS1′能真實地模擬出Pan波段的灰度信息，Pan－MS1′才能確切地反映出Pan波段數據的高頻信息，否則其中必然包含不一致的灰度信息，從而造成融合影像光譜失真。

1.2 兩種典型的構造可替換波段的分量替換法

基于光譜響應函數模擬低分辨率全色波段的GS融合算法［6］和基于多元線性回歸方法模擬低分辨率全色波段的GS融合算法［5］是兩種十分有代表性的通過構造可替換波段改進光譜失真的分量替換方法，以下對這兩種方法的原理進行說明。

1.2.1 基于光譜響應函數模擬低分辨率全色波段的GS融合算法 理解光譜響應函數對傳感器接收到的輻射亮度值的影響是說明基于光譜響應函數模擬低分辨率全色波段的GS融合算法的前提，因此，本節首先說明什么是光譜響應函數以及光譜響應函數對影像融合的影響，再解釋該融合算法的具體步驟。光譜響應函數反映的是傳感器的一種靜態性能，遙感器的光譜響應特征是波長的函數，是遙感器在每個波長處接收輻亮度與入射輻亮度的比值。不同傳感器光譜響應函數曲線的形狀、波段中心波長位置、帶寬、波段之間交迭程度等的不同都將導致它們觀測得到的波段反射率不一致［8］。遙感器每個波段都有一定的波長響應寬度，遙感器的每個波段實際接收到的能量是該波段的波長范圍內各個波長處接收的能量的和。通常采用對遙感器接受的輻亮度按照光譜響應函數進行歸一化來降低光譜響應特征不同的影響［8］，如（4）式所示。

式中：f（λ）——遙感器某波段的光譜響應函數；L（λ）——入射到遙感傳感器所接收到并記錄的輻亮度值；λmax，λmin——波段范圍的上下界。

為了說明光譜響應函數對遙感圖像融合的影響，本文以QuickBird數據為例，表1表示QuickBird各波段波長范圍和中心波長，圖1表示QuickBird的光譜響應特征［9］，從表1和圖1可以看出Blue、Green、Red和NIR各波段之間都有交迭，而且Pan波段上界超出了NIR的上界，全色波段與多光譜波段光譜響應函數不完全一致。如果按照傳統的GS融合方法，僅僅通過簡單地選取多光譜圖像的均值來模擬低分辨率全色波段，該方法雖然快速簡單，但不能準確獲取Pan波段數據在低分辨率情況下的灰度信息。

表1 QuickBird波段范圍和中心波長［9］

為了避免上述問題，Gonzales［6］通過對 MS與Pan傳感器的光譜響應函數進行線性回歸，獲取使用MS模擬低分辨率Pan的線性組合系數，再將模擬的低分辨率全色波段應用到GS融合模型中。具體步驟如下：

（1）首先對MS和Pan的光譜響應函數進行線性回歸，獲取各波段的擬合系數 ，通過對MS各波段的光譜響應函數的線性組合模擬Pan波段的光譜響應函數，見式（5）：

其中：RPAN′（λ）為模型的Pan波段的光譜響應函數，RB，RG，RR，RNIR分別為紅、綠、藍、近紅外波段的光譜響應函數。

圖1 QuickBird多光譜波段和全色波段的光譜響應函數［9］

（2）由于模擬Pan波段光譜響應函數和Pan波段很接近，因此由模擬的Pan波段光譜響應函數所記錄的光譜輻射值與Pan波段數據具有較高的可比性，這一模擬的低分辨率全色波段輻射亮度值可以表示為

式中：LPAN′（λ）——模擬的Pan波段的輻射亮度值；LB，LG，LR，LNIR——紅、綠、藍、近紅外波段的輻射亮度值；ci（i＝1，2，3，4）——擬合系數。

（3）最后，將模擬的低分辨率Pan波段數據應用到GS融合模型中。

1.2.2 基于多元線性回歸方法模擬低分辨率全色波段的GS融合算法 Aiazzi等［5］在研究改進的分量替換方法中利用最小二乘法，直接通過全色波段和多光譜波段進行線性回歸獲取擬合系數來模擬出低分辨率全色波段。算法的具體步驟如下：

（1）將原始Pan波段數據通過低通濾波重采樣為分辨率和原始多光譜波段分辨率一致的圖像，記為P。假設

式中：I——所要模擬的低分辨率全色波段I的估計量，通過線性回歸方法最小化P和I之間的PMSE獲取估計量w1，w2，w3，w4和b。

（2）通過以上的擬合系數模擬出低分辨率全色波段

（3）最后，將模擬的低分辨率全色波段數據應用到GS融合模型中。

2 融合實驗結果與分析

為了進一步理解分量替換融合思想，并檢驗改進的分量替換融合方法在光譜保真性上的改善。本文進行三種典型的分量替換融合方法的對比分析，分別是PCA融合法、基于光譜響應函數模擬低分辨率全色波段的GS融合法（以下簡稱GS1）以及基于多元一次線性回歸擬合低分辨率全色波段GS融合法（以下簡稱GS2）。

2.1 數據說明

遙感數據融合前的輻射水準歸一化以及參與融合的各波段數據的精確配準是融合前的重要基礎工作。論文所選取的數據為DigitalGlobe公司的Demo數據（qb＿boulder＿msi和qb＿boulder＿Pan），多光譜圖像大小為1 024×1 024像素，多光譜圖像空間分辨率為2.8m，全色圖像空間分辨率為0.7m。圖像的輻射分辨率為11bit，以16bit的格式發放，數據的動態范圍為0～2 047。影像級別為標準產品（Standard Imagery Products）：具有地理參考和地圖投影，經過了輻射校正、傳感器和衛星平臺引起的系統誤差校正，利用粗DEM做了地形誤差校正，經過了全色和多光譜波段之間的配準。由于所選取的數據已經經過了輻射校正和各波段間的配準，可直接用于本文的融合實驗。

本文選取ENVI 5.0軟件實現PCA和GS融合實驗。目前集成在ENVI 5.0軟件中的GS算法有4種可選擇的模式用來模擬低分辨率全色波段，模式一：選取多光譜波段的均值來模擬低分辨率全色波段；模式二：輸入現有的單波段數據來模擬低分辨率全色波段；模式三：基于不同傳感器光譜響應函數模擬低分辨率全色波段；模式四：自定義低通濾波函數來模擬低分辨率全色波段。本實驗選擇模式三進行GS1融合實驗，選擇模式四進行GS2融合實驗。

2.2 融合效果對比

融合效果評價是一個復雜的問題，均值差（Mean Bias，MB）、方差差異（Variance Difference，VD）、標準偏差差異（Standard Deviation Difference，SDD）、相 關 系 數 （Correlation Coefficient，CC）、光 譜 角（Spectral Angle Mapper，SAM）、相對量綱全局誤差（Relative Dimensionless Global Error，ERGAS）、Q4質量評價指標（Q4Quality Index，Q4）是目前遙感領域常用的7個融合質量評價指［10－13］。Zhang［14］通過比較實驗發現這7個指標的定量評價結果和目視比較分析結果以及數字分類結果存在嚴重的不一致，他認為目前國際上常用上述7個指標并不能作為評價遙感圖像融合質量的可靠標準。

因此本研究沒有使用7個評價指標，而是從融合前后圖像的總體和細節（典型地物）的光譜特征展開比較：（1）從總體上進行比較，包括融合影像與原始影像的相關系數、RMSE和灰度直方圖的比較；（2）融合影像與原始影像典型地物的光譜特性比較，所選取的典型地物既有高反射率地物（如沙地、水泥地），也有低反射率地物（如水體），還包括植被。依此來考察經過融合處理后的影像在整個譜段范圍的光譜保真性。

2.2.1 融合影像與原始影像的總體比較 從融合目視效果看，三種方法融合處理后的影像都紋理清晰，地物清晰可辨，PCA和GS2融合影像顏色上和原多光譜圖像非常一致，但GS1融合后影像上植被的顏色明顯比原多光譜圖像亮。下文從融合影像與原始影像的灰度直方圖、相關系數和RMSE三方面來做比較分析。（1）圖2顯示，從總體上，三種方法融合后的圖像和原多光譜圖像各波段的灰度直方圖較為一致，并且各融合圖像之間各波段的灰度直方圖極其一致，幾乎重疊在一起。（2）從所獲取的相關系數（表2）來看，相關系數都接近1，其中GS2融合后Band1，Band2，Band3和原多光譜圖像的相關系數最高（0.917 1，0.919 8，0.921 4），PCA融合后Band4和原多光譜圖像的相關系數最高（0.931 0）。（3）表3中所有RMSE≤2.5%灰度范圍，其中PCA融合后的Band1、Band2和原多光譜圖像的RMSE最小（10.906 0，22.321 3），GS2融合后的 Band3、Band4和原多光譜圖像的RMSE最小（28.572 7，46.876 8）。

從總體的比較分析發現，三種融合方法都具有很高的光譜保真性，其中GS2具有最優的光譜保真性，PCA和GS1融合算法次之。

表2 融合影像與原始影像的相關系數

表3 融合影像與原始影像的RMSE

圖2 原始圖像和融合后圖像的灰度直方圖

2.2.2 融合影像與原始影像典型地物的光譜特性比較 為了進一步更全面考察經過融合處理后影像的光譜保持性，本研究選取純凈的水體（代表低反射率地物）、植被和建筑（代表高反射率地物）三種典型地物來探究融合前后影像光譜特性的變化。盡量選取同質均勻的典型地物區域，以利于融合前后典型地物光譜特性的對比。

（1）PCA、GS1、GS2融合前后水體的光譜特性比較。根據水體的光譜反射率隨著波長的增大逐漸降低的特性，本研究選取水體的綠光波段（Green）和近紅外波段（NIR）進行光譜分析，圖3表示融合前后水體在Green和NIR波段水平方向的光譜剖面曲線，分析圖3可以得出以下結論：① 經過PCA融合處理后圖像水體的Green和NIR波段灰度值略高于原圖像；② 經過GS1、GS2融合后的圖像和原始圖像水體的Green波段灰度值幾乎完全一致；③ 經過GS1融合后的圖像水體的NIR波段灰度值和融合前最為接近。因此，GS1融合方法具有最優的水體光譜保真性，GS2次之，PCA融合方法在水體的綠光和近紅外波段的光譜反射率略高于原圖像。

（2）PCA、GS1、GS2融合前后植被的光譜特性比較。根據植被在紅光（Red）吸收而在近紅外（NIR）波段反射劇增的光譜特性，論文中選取Red和NIR波段進行植被的光譜分析，圖4表示融合前后植被在Red和NIR波段的空間光譜剖面曲線（利用ENVI軟件繪制圖像的線向量的空間光譜剖面曲線），分析融合前后Red和NIR的變化可以得出以下結論：①經過GS1融合處理后的圖像植被的Red和NIR波段明顯高于原圖像和其他融合結果，這一結論和Konstantinos［15］的結論相一致，Konstantinos認為融合后植被覆蓋地區能夠更容易從耕地和建筑物中區分出來；②GS2融合處理后的圖像植被的Red和NIR波段和原圖像最接近，PCA融合方法次之。

（3）PCA，GS1，GS2融合前后建筑物的光譜特性比較。選取Green和NIR波段進行建筑物的光譜分析，圖5為融合前后建筑在Green和NIR波段的空間光譜剖面曲線（利用ENVI軟件繪制圖像的線向量的空間光譜剖面曲線），經過分析有如下結論：① 經過GS1融合處理后的圖像建筑物在Green和NIR波段的灰度值都明顯低于原圖像；② 經過PCA和GS2融合處理后的圖像建筑物在Green波段和NIR波段灰度值略低原圖像，其中PCA融合結果和原圖像最為接近。

圖3 融合前后水體在Green和NIR波段水平方向的光譜剖面曲線

圖4 融合前后植被在Red和NIR波段的空間光譜剖面曲線

圖5 融合前后建筑在Green和NIR波段的空間光譜剖面曲線

3 結論

（1）從總體的比較（融合影像與原始影像的相關系數、RMSE和灰度直方圖）和典型地物在融合前后的光譜保持性的比較一致發現，三種融合方法都具有很高的光譜保真性，其中基于多元一次線性回歸擬合低分辨率全色波段的GS2融合方法具有最優的光譜保真性，PCA和GS1融合方法次之。尤其值得注意的是GS1融合處理后的圖像存在部分光譜失真的現象。經過GS1融合處理后的圖像，植被覆蓋地區Red和NIR波段光譜反射率明顯高于原始的多光譜圖像，而高反射率區域（如沙地、水泥建筑）Green和NIR波段光譜反射率則明顯低于原始的多光譜圖像。

（2）造成GS2融合算法明顯優于GS1的原因是GS1融合算法使用的只是實驗室理想環境下所獲取的名義上的光譜響應特征。傳感器的實際成像過程受到諸多因素的影響，如傳感器在軌工作環境、大氣的影響、觀測角度不同等的影響。GS1算法單純通過不同波段的光譜響應函數的線性組合來模擬低分辨率全色波段并不十分準確，GS2直接利用MS和Pan波段像元灰度值進行線性回歸，克服了上述不確定性問題。

（3）通過對PCA、GS1和GS2三種融合算法的光譜保真性的比較分析說明，如何利用多光譜數據準確地模擬低分辨率全色波段，直接影響到融合后影像的光譜保真性，是目前高分辨率遙感圖像融合的關鍵技術。由多光譜數據模擬低分辨率全色波段的方法引起了越來越多研究人員的重視，也代表了目前高分辨率遙感圖像融合的研究方向。

［1］ Aiazzi B，Alparone L，Baronti S，et al.25years of pansharpening：a Critical Review and New Developments［M］∥Signal and Image Processing for Remote Sensing.2nd E－dition.Boca Raton，FL：CRC Press，2011：533－548.

［2］ Carper W，Lillesand T，Kiefer R.The use of intensityhue－saturation transformations for merging SPOT panchromatic and multispectral image data［J］.Photogrammetic Engineering and Remote Sensing，1990，56（4）：459－467.

［3］ Chavez P S，Sides J S C，Anderson J A.Comparison of three different methods to merge multiresolution and mutispectral data：Landsat TM and SPOT panchromatic［J］.Photogrammetic Engineering and Remote Sensing，1991，57（3）：295－303.

［4］ Lanben C A，Brower B V.Process for enhancing the spatial resolution of multispectral imagery using Pansharpening［P］.Tech.Rep.，Eastman Kodak Company：US，6011875.2000.

［5］ Aiazzi B，Baronti S，Selva M.Improving component substitution pansharpening through multivariate regression of MS＋Pan data［J］.Geoscience and Remote Sensing，IEEE Transactions on，2007，45（10）：3230－3239.

［6］ González－Audícana M，Otazu X，Fors O，et al.A low computational－cost method to fuse IKONOS images using the spectral response function of its sensors［J］.Geoscience and Remote Sensing，IEEE Transactions on，2006，44（6）：1683－1691.

［7］ Dou W，Chen Y，Li X，et al.A general framework for component substitution image fusion：An implementation using the fast image fusion method［J］.Computers＆ Geosciences，2007，33（2）：219－228.

［8］ 梁順林.定量遙感［M］.北京：科學出版社，2009：157－176.

［9］ DigitalGlobe，Inc.Spectral Response for DigitalGlobe Earth Imaging Instruments［EB／OL］.2011.http：∥www.digitalglobe.com／downloads／DigitalGlobe＿Spectral＿Response.pdf.

［10］ Wald L，Ranchin T，Mangolini M.Fusion of satellite images of different spatial resolutions：assessing the quality of resulting images［J］.Photogrammetric Engineering ＆ Remote Sensing，1997，63（6）：691－699.

［11］ Wang Z，Ziou D，Armenakis C，et al.A comparative analysis of image fusion methods［J］.Geoscience and Remote Sensing，IEEE Transactions on，2005，43（6）：1391－1402.

［12］ Wang Z，Bovik A C，Sheikh H R，et al.Image quality assessment：From error visibility to structural similarity［J］.Image Processing，IEEE Transactions on，2004，13（4）：600－612.

［13］ Wang Z，Bovik A C.A universal image quality index［J］.Signal Processing Letters，IEEE，2002，9（3）：81－84.

［14］ Zhang Y.Methods for image fusion quality assessment－a review，comparison and analysis［J］.The International Archives of the Photogrammetry，Remote Sensing and Spatial Information Sciences，2008，37：1101－1109.

［15］ Nikolakopoulos K G.Spatial resolution enhancement of Hyperion hyperspectral data［C］∥Hyperspectral Image and Signal Processing：Evolution in Remote Sensing，2009.WHISPERS′09.First Workshop on.IEEE，2009：1－4.