王慧賢，江萬壽，雷呈強，張一鳴

1.武漢大學 測繪遙感信息工程國家重點實驗室，湖北 武漢 430079；2.軍械工程學院 電子與光學工程系，河北 石家莊 050003；3.北京市遙感信息研究所，北京 100192

1 引 言

隨著對地觀測技術的發(fā)展，大多數對地觀測衛(wèi)星如Landsat－7、Spot 5、IKONOS、QuickBird、GeoEye－1和WorldView－2可以同時提供高光譜分辨率的多光譜影像和高空間分辨率的全色影像。為了充分利用光譜和空間特性，影像融合技術得到了迅速的發(fā)展。由文獻［1—2］可知，目前的遙感影像融合方法可以分為投影替換類、相對光譜分布類、結構注入空間分辨率增強類ARSIS（amélioration de la résolution spatiale par injection de structures）及混合模型法4類。以IHS（intensity－h(huán)ue－saturation）［3］、PCA （principle component analysis）［4］、gram－schmidt spectral［5］、HCS（hyperspherical color sharpening）［6］和偏微分替換［7］為代表的投影替換類算法，利用全色波段或其變體替換多光譜影像投影變換后得到的亮度分量，然后進行投影逆變換獲得融合結果；以Brovey［8］、合成變量比（synthetic variable ratio，SVR）［9］為代表的相對光譜分布類算法，基于全色影像與多光譜影像之間的線性組合假設，把全色影像信息按比例分布到多光譜影像的各個波段；以 Wavelets［10－11］、Curvelets［12］、HPF（high－pass filtering）［13］、multiscale toggle contrast operator［14］、Markov random field（MRF）models［15］和多尺度光譜增益調制［16］為代表的ARSIS類算法在不同尺度上提取影像的空間信息和光譜信息進行影像融合。投影替換和相對光譜分布這兩類算法能夠較好地保持空間分辨率，但都依賴于全色影像與多光譜影像的相關性，相關性越好融合結果越好，且光譜信息的保持依賴于傳感器的配置和影像信息本身。ARSIS類算法能夠較好地保持光譜特性，但是全色影像和多光譜影像空間信息的不一致容易造成空間細節(jié)信息的混疊，導致空間分辨率的降低。混合模型算法［17］則試圖吸收前3類算法在空間信息保持或光譜信息保持的優(yōu)勢，并引入 Map估計［18］、稀疏表達［19－20］、壓縮感知［21－22］等新理論，該類算法研究處于起步階段，算法也比較復雜，離實用化還有一定距離。

總之，相對光譜分布類的代表性算法SVR算法具有明確的物理意義和理論基礎，在全色波段和多光譜波段（或部分波段）能夠相互覆蓋的條件下（如IKONOS、QuickBird），不僅能保持全色影像的空間細節(jié)，也能很好地保持多光譜影像（multispectral image，MS）的光譜信息。然而，對于不滿足光譜覆蓋條件的影像，SVR會造成融合影像的光譜扭曲［23］。另外，對于滿足光譜覆蓋條件的影像，也可能會由于影像上地物類型分布不均衡，導致組合系數求解的精度不夠，使合成的全色波段（panchromatic image，PAN）影像與實際不符，造成融合影像的光譜失真。針對以上兩個問題，本文提出一種基于光譜和空間局部相關的改進算法SSCSVR（spectral and spatial correlation－based synthetic variable ratio）。

2 基于光譜與空間局部相關融合方法原理

2.1 SVR算法原理與問題分析

文獻［24］提出的SVR影像融合方法，出發(fā)點是低分辨率全色影像能夠利用多光譜影像合成

式中，PanSyn是合成的全色影像灰度值；MSi是要融合的多光譜影像中第i波段的灰度值；φi是第i波段的權系數值。

另一方面，SVR方法認為高分辨率全色影像和合成的低分辨率全色影像的比值反映了全色影像和多光譜影像的細節(jié)差異，如果把這些差異按比例分布到多光譜影像的各個波段，就可以生成和全色影像分辨率一致的高分辨率多光譜影像

式中，F(xiàn)usedi是融合后影像中第i波段的灰度值；Panori是原始的全色影像灰度值。

根據全色影像和多光譜影像的線性組合假設，全色影像合成的權系數φi可通過式（3）得到

文獻［24］利用城區(qū)、土壤、水、草地和樹木5類地物來求取權系數φi。

文獻［9］對SVR的方法進行進一步研究，明確了合成全色影像的多光譜波段應該是全色波段覆蓋的波段，而不一定是多光譜影像的全部波段。同時，文獻［9］認為求合成系數φi不能僅僅由指定的5類地物來求解，而應該利用影像的全部整體信息來求解。為此，文獻［9］的改進算法采用了直方圖歸一化、融合后影像直方圖重新拉伸等處理，使SVR成為一類具有獨特意義的實用算法。然而對于不滿足光譜覆蓋條件的影像，SVR［9］會造成融合影像的光譜扭曲［23］。另外，對于滿足光譜覆蓋條件的影像，也可能會由于影像上地物類型分布的不均衡，導致組合系數求解的精度不夠，使合成的PAN影像與實際不符，造成融合影像的光譜失真。

綜合分析文獻［9］的SVR算法，導致偏色的原因有：① 合成的PAN只考慮光譜相關而沒有考慮空間相關，導致合成PAN不準確，引起偏色，尤其當多光譜波段與PAN波段不能完全相互覆蓋時，偏色情況更為嚴重；② 對多光譜和全色波段進行直方圖的歸一化，這樣可能會引起光譜信息的損失，引起偏色問題；③ 算法中要求合成的PAN與原始的PAN灰度值盡量接近，但是利用影像的整體信息來計算權值，會引起局部模擬的全色和原始全色影像光譜有較大不同，從而引起偏色問題。

本文基于光譜與空間局部相關的算法就是針對上述問題提出的。2.2.1節(jié)提出的光譜與空間相關模型能夠更加嚴密地合成PAN值，克服第1個偏色原因；2.2.2節(jié)提出的自適應的局部處理和2.2.3節(jié)提出的非負約束最小二乘求解及異常處理方法可以有效地解決第2、3個偏色問題。

2.2 改進的基于光譜與空間局部相關的SVR方法融合原理

2.2.1 改進的光譜與空間相關模型

融合中光譜的保持與多光譜波段結構有關，而空間細節(jié)的保持與波段內的空間結構有關。計算PanSyn時，傳統(tǒng)的SVR只考慮了PAN與多光譜波段的關系，但影像中像元之間是空間相關的，如果一個像元對融合影像有貢獻，這個像元周圍的像元也會對融合影像有貢獻。PAN波段覆蓋范圍一般較大，不僅含有多光譜影像的光譜信息而且含有多于多光譜影像的空間細節(jié)信息?；谝陨戏治?，本文提出了改進的空間與光譜相關的回歸模型，模型見公式（4）。該模型能夠更好地描述高空間分辨率的PAN與低分辨率的多光譜之間的關系，得到更準確的權值，從而得到更加準確的合成PAN值

式中，β是空間相關成分的權值；Gs是空間相關成分信息即高頻成分。式（4）的矩陣形式如式（5）所示

式中，X＝［φ1φ2…φNβ］T為待求的權系數矩陣；N為多光譜影像的波段數，B＝［MS1MS2…MSNGs］為已知的光譜和空間成分。

為了更好地獲得高頻成分Gs，筆者選用了高斯函數卷積的方法。高斯函數具有低通性質，并且可以模擬人眼視覺機理，通過不同的尺度參數σ設定，能夠模擬人眼遠近處所看到的圖像。高斯函數表達式為

式中，系數K為歸一化值。將高斯函數與影像進行卷積，可以得到低頻信息，從原始影像中減去低頻信息，即可得到高頻信息Gs，?表示卷積操作，Gs的計算見式（7）

高斯核窗口r和尺度參數σ會直接影響到提取空間細節(jié)成分，進而影響光譜與空間相關權系數的求解。尺度參數σ越小提取的空間細節(jié)信息越少，反之尺度越大，提取的空間細節(jié)信息越多，但有可能出現(xiàn)過提取的情況。綜合考慮空間細節(jié)提取的程度，本文采用標準高斯分布的σ值即σ＝1。高斯核窗口r越大，Gs高頻細節(jié)信息越多，但是計算量大，反之亦然，根據誤差3σ規(guī)則，本文高斯核窗口r＝3σ。

2.2.2 自適應的局部處理方法

在地物類型分布嚴重不均衡的情況下，利用影像的整體信息來計算權值，會引起合成的全色和原始全色影像光譜存在較大不同，從而引起偏色問題。考慮到影像像元的空間差異性，局部處理可能會比全局處理具有更大的優(yōu)勢。因此，本文將整個影像分為很多影像塊，利用每個影像塊里面的像素分別建立2.2.1中的光譜與空間相關回歸模型，并采用非負的最小二乘方法進行求解，得到每個影像塊里像素的光譜與空間相關權值，從而得到合成的全色PAN值。設h和l分別為PAN和MS的空間分辨率，則r＝l／h為兩者空間分辨率的比值，綜合考慮融合效果與處理效率，分塊大小w本文試驗中取w＝5r＋1。

局部分塊優(yōu)化處理往往會產生分塊效應，綜合考慮融合效果和時間復雜度的影響，采用雙線性內插來獲得每個像素的權值，從而得到每個像素的合成PAN值。

2.2.3 非負約束最小二乘求解及異常情況處理

本文采用非負約束最小二乘方法進行求解，可以使求解的系數非負，使系數具有物理意義。非負約束最小二乘是一種較為成熟的優(yōu)化理論，原理描述如下

式中，LSE為非負約束的優(yōu)化誤差；X≥0表示對于所有的1≤k≤N＋1都有xk≥0。為了應用拉格朗日方法解求優(yōu)化問題，引入N＋1維未知的被動約束常量c＝［c1c2…cN＋1］T＞0，ck＞0，k的取值為1，2，…，N＋1。通過c，可以形成拉格朗日式J

從式（11）可得出式（12）和式（13）

式（12）和（13）可以用來解求非負約束優(yōu)化解＾XNCLS和拉格朗日系數λ＝［λ1λ2…λN＋1］T。

本文求解系數采用了快速非負約束最小二乘（fast non－negativity－constrained least squares，F(xiàn)NNLS）算法，詳細介紹見文獻［25］。該算法的主要思想是將估計量分離為活動集R和被動集P，前者包括中所有負值和零值，后者包括所有剩余的正值。FNNLS處理開始時被動集P為空集，假設活動集R包括中所有元素，通過式（12）和（13）迭代處理調節(jié)活動集和被動集。具體實施時，迭代結束條件可設為活動集R為空集或者活動集R不為空集且所有拉格朗日系數λk，k∈R為負值。

利用FNNLS算法求解權系數矩陣過程中可以看出需對實對稱矩陣求逆，此時，要考慮該矩陣是否可逆的問題，如果該矩陣可逆，則可以利用FNNLS方法較容易地求得準確的權系數，反之則需進行異常情況的處理。通過試驗發(fā)現(xiàn)實對稱矩陣不可逆是由于用來求解系數的像素灰度值變化不大，非常接近，如大片水的區(qū)域容易發(fā)生這種情況。當系數求解不出時，采用式（14），即系數等于塊內所有全色像素的和的比值與所有多光譜像素的和

式中，np為用于求解系數的像素個數；j的取值范圍為1，2，…，np；Panj是全色影像第j個像素的灰度值；MSi，j是多光譜影像第i波段第j個像素的灰度值。一般用于求解系數的像素采用影像塊內所有像素，np＝w×w。

3 試驗結果與分析

選取意大利羅馬地區(qū)的 WorldView－2影像作為試驗數據，影像數據包括1個空間分辨率0.5m的全色波段和8個空間分辨率2m的多光譜波段，全色影像大小為4600像素×4604像素，多光譜影像大小為1150像素×1151像素。試驗區(qū)域為城市區(qū)域，包括植被、建筑物、道路、水體等地物類型。

圖1（a）為試驗區(qū)的待融合全色影像，圖1（b）為經過嚴格空間配準和重采樣后的待融合多光譜影像。為了評價本文方法，采用Gram－Schmidt［5］、HCS［6］和 SVR［9］方法進行對比。圖1（c）為本文SSCSVR方法融合后的影像。鑒于篇幅限制，在此只列出本文方法融合后的結果。圖1（d）～（i）分別為待融合全色影像（圖中綠框區(qū)域）、多光譜影像以及本文SSCSVR、HCS、Gram－Schmidt和SVR［9］融合后影像的局部放大圖。影像以波段（5，3，2）顯示的（R，G，B）。與原始多光譜影像相比，4種融合方法的空間分辨率都有明顯的提高，細節(jié)得到了增強，達到了很好的空間分辨效果。但是從光譜色彩來看，Gram－Schmidt和SVR［9］方法整體光譜扭曲比較嚴重，整個影像偏紅，灰色的道路顏色偏白，暗黑色區(qū)域偏藍。分析SVR［9］方法偏色的原因主要是：影像中房屋紅色偏多，Zhang SVR利用整體信息求解出的權系數紅色成分偏大，導致道路和暗黑色小區(qū)域嚴重偏色，并且使影像整體偏色嚴重。WorldView－2的全色波段與多光譜（或部分波段）不能互相完全覆蓋，導致SVR［9］求解出的系數不能準確描述全色與多光譜影像之間的關系。本文采用光譜與空間相關局部優(yōu)化策略，克服了以上缺點。HCS是處理 WorldView－2衛(wèi)星影像的官方處理方法，本文方法不僅光譜保持程度可以和HCS方法相媲美，能保持地物的原來色彩，達到了很好的視覺效果，而且在空間細節(jié)方面比HCS要好，如圖紅框區(qū)域，通過仔細查看可以看出HCS方法處理結果丟失了部分細節(jié)，而本文SSCSVR方法細節(jié)明顯。

ERGAS（relative global－dimensional synthesis error）［26］可以客觀地評價融合影像的整體光譜質量，偏差指數反映影像的偏色水平，SSIM（structure similarity index Measurement）［27］可以較好地評價光譜結構特征，信息熵反映影像空間細節(jié)情況。這4個客觀評價指標可以較好地評價融合后影像的光譜和空間特征，且較為常用。因此，為了客觀定量評價文中各方法，本文選取了這4個評價指標。ERGAS和偏差指數值越小說明光譜保持越好；SSIM值越大說明光譜結構越相似，質量越好；信息熵反映影像空間細節(jié)情況，熵值越大說明影像所含的信息量越多。原始多光譜影像平均信息熵5.382，全色影像信息熵5.567。表1列出了各種融合方法定量評價的結果，從表中可以看出各種方法的融合結果信息熵都有所提高，從而到達了融合的目的。Gram－Schmidt和SVR［9］的各項指標比較相近，ERGAS和偏差指數較高，而SSIM值偏低，說明這兩種方法光譜扭曲程度較大，而信息熵值較低，說明影像包含信息較少。從ERGAS、SSIM和偏差指數可以看出HCS和SSCSVR兩個方法都能較好地保持光譜特征且能夠保持原來光譜的內部結構，從信息熵可以看出SSCSVR比HCS信息含量稍微高一點，充分驗證了目視評價中HCS空間細節(jié)稍遜于本文方法SSCSVR的說法。

圖1 WorldView－2影像融合結果Fig.1 Fusion results of WorldView－2images

表1 各方法融合結果的定量評價Tab.1 Quantitative assessments of fusion results

從目視和客觀評價可以看出本文SSCSVR方法較SVR［9］光譜和空間細節(jié)方面都有一定提高，充分克服了SVR方法本身光譜扭曲的缺陷。

為了充分說明本文方法的正確性，筆者又選用了另外一組QuickBird的印度地區(qū)數據，影像包括5個波段，4個多光譜波段，空間分辨率2.44m，1個全色波段，空間分辨率0.61m。全色影像大小為6000像素×6000像素，多光譜影像大小為1500像素×1500像素。試驗區(qū)域為城市區(qū)域，包括水體、植被、建筑物和道路等地物類型。融合結果如圖2所示，以波段（3，2，1）顯示的（R，G，B）。

圖2 QuickBird影像融合結果Fig.2 Fusion results of QuickBird images

從圖2可以看出，SVR［9］方法存在較為嚴重的光譜偏差，整體影像偏綠，偏白，HCS綠色有點偏藍，Gram－Schmidt方法影像偏綠，圍欄區(qū)域偏白而SSCSVR光譜保持較好，4種方法空間細節(jié)上看不出什么明顯的差別。定量比較結果見表2，SSCSVR各項指標都最優(yōu)。從綜合上來看，SSCSVR的融合結果是最好的，這一結論與WorldView－2融合試驗的結論一致。通過試驗可以證明，即使當全色波段與多光譜波段完全覆蓋時，由于影像質量本身的原因（QuickBird影像中地物類型分布不均衡，全色影像細節(jié)較多），利用整體信息和只考慮光譜相關求解出的權系數不準確，導致SVR［9］的融合結果不是很理想，而本文方法則較好地克服了這一點。

表2 各方法融合結果的定量評價Tab.2 Quantitative assessments of fusion results

筆者也對其他傳感器影像如Spot－5、IKONOS、Landsat－7、GeoEye－1作了類似試驗，可以得出與文中一樣的結論，驗證了本文方法的優(yōu)越性。

4 結 論

本文針對傳統(tǒng)SVR方法存在光譜扭曲問題，提出了一種基于光譜與空間局部相關的光譜和空間保持型影像融合方法。與傳統(tǒng)SVR方法不同的是，該方法不需要將全色和多光譜影像進行直方圖歸一化，因此能更好地保護空間和光譜信息，也無須嚴格滿足全色與多光譜相互覆蓋的條件。該方法用于 WorldView－2等衛(wèi)星的全色和多光譜影像融合試驗中，能夠得到空間細節(jié)與光譜信息保持度高的融合影像，能適應不同傳感器的特點和地物分布情況。

［1］ THOMAS C，RANCHIN T，WALD L，et al.Synthesis of Multispectral Images to High Spatial Resolution：A Critical Review of Fusion Methods Based on Remote Sensing Physics［J］.IEEE Transactions on Geoscience and Remote Sensing，2008，46（5）：1301－1312.

［2］ BARONTI S，AIAZZI B，SELVA M，et al.A Theoretical Analysis of the Effects of Aliasing and Misregistration on Pansharpened Imagery［J］.Journal of Selected Topics in Signal Processing，2011，5（3）：446－453.

［3］ CHIKR EL－MEZOUAR M，TALEB N，KPALMA K，et al.An IHS－Based Fusion for Color Distortion Reduction and Vegetation Enhancement in IKONOS Imagery［J］.IEEE Transactions on Geoscience and Remote Sensing，2011，49（5）：1590－1602.

［4］ PAT S C J，STUART S C，JEFFREY A A.Comparison of Three Different Methods to Merge Multiresolution and Multispectral Data：Landsat TM and SPOT Panchromatic［J］.Photogrammetric Engineering and Remote Sensing，1991，57（3）：295－303.

［5］ LABEN C A，BROWER B V.Process for Enhancing the Spatial Resolution of Multispectral Imagery Using Pan－Sharpening：United States，6011875［P］.2000.

［6］ PADWICK C，PACIFICI F，SMALLWOOD S.WorldView－2 Pan－Sharpening［C］∥ASPRS 2010Annual Conference.San Diego：［s.n.］，2010.

［7］ JAEWAN C，KIYUN Y，YONGIL K.A New Adaptive Component－Substitution－Based Satellite Image Fusion by Using Partial Replacement［J］.IEEE Transactions on Geoscience and Remote Sensing，2011，49（1）：295－309.

［8］ GILLESPIE A R，KAHLE A B，WALKER R E.Color Enhancement of Highly Correlated Images.II.Channel Ratio and"chromaticity"Transformation Techniques［J］.Remote Sensing of Environment，1987，22（3）：343－365.

［9］ ZHANG Yun.System and Method for Image Fusion：United States，7340099B2［P］.2008.

［10］ OTAZU X，GONZALEZ－AUDICANA M，F(xiàn)ORS O，et al.Introduction of Sensor Spectral Response into Image Fusion Methods.Application to Wavelet－based Methods［J］.IEEE Transactions on Geoscience and Remote Sensing，2005，43（10）：2376－2385.

［11］ PRADHAN P S，KING R L，YOUNAN N H，et al.Estimation of the Number of Decomposition Levels for a Wavelet－Based Multiresolution Multisensor Image Fusion［J］.IEEE Transactions on Geoscience and Remote Sensing，2006，44（12）：3674－3686.

［12］ CUNHA A L，ZHOU Jianping，DO M N.The Nonsubsampled Contourlet Transform：Theory，Design，and Applications［J］.IEEE Transactions on Image Processing，2006，15（10）：3089－3101.

［13］ GANGKOFNER U G，PRADHAN P S，HOLCOMB D W.Optimizing the High－Pass Filter Addition Technique for Image Fusion［J］.Photogrammetric Engineering and Remote Sensing，2008，74（9）：1107－1118.

［14］ BAI Xiangzhi，ZHOU Fugen，XUE Bindang.Edge Preserved Image Fusion based on Multiscale Toggle Contrast Operator［J］.Image and Vision Computing，2011，29（12）：829－839.

［15］ XU Min，CHEN Hao，VARSHNEY P K.An Image Fusion Approach Based on Markov Random Fields［J］.IEEE Transactions on Geoscience and Remote Sensing，2011，49（12）：5116－5127.

［16］ LIU Jun，SHAO Zhenfeng.Remote Sensing Image Fusion Using Multi－scale Spectrum Gain Modulation［J］.Acta Geodaetica et Cartographica Sinica，2011，40（4）：470－476.（劉軍，邵振峰.基于多尺度光譜增益調制的遙感影像融合方法［J］.測繪學報，2011，40（04）：470－476.）

［17］ CHU Heng，WANG Ruyan，ZHU Weile.Remote Sensing Image Fusion Algorithm in the DCT Domain［J］.Acta Geodaetica et Cartographica Sinica，2008，37（1）：70－76.（楚恒，王汝言，朱維樂.DCT域遙感影像融合算法［J］.測繪學報，2008，37（1）：70－76.）

［18］ JOSHI M，JALOBEANU A.MAP Estimation for Multiresolution Fusion in Remotely Sensed Images Using an IGMRF Prior Model［J］.IEEE Transactions on Geoscience and Remote Sensing，2010，48（3）：1245－1255.

［19］ YU Nannan，QIU Tianshuang，BI Feng，et al.Image Features Extraction and Fusion Based on Joint Sparse Representation［J］.IEEE Journal of Selected Topics in Signal Processing，2011，5（5）：1074－1082.

［20］ YANG Bin，LI Shutao.Pixel－Level Image Fusion with Simultaneous Orthogonal Matching Pursuit［J］.Information Fusion，2012，13（1）：10－19.

［21］ LI Shutao，YANG Bin.A New Pan－Sharpening Method Using a Compressed Sensing Technique［J］.IEEE Transactions on Geoscience and Remote Sensing，2011，49（2）：738－746.

［22］ JIANG Cheng，ZHANG Hongyan，SHEN Huanfeng，et al.A Practical Compressed Sensing－Based Pan－Sharpening Method［J］.Geoscience and Remote Sensing Letters，2012，9（4）：629－633.

［23］ ALPARONE L，WALD L，CHANUSSOT J，et al.Comparison of Pansharpening Algorithms：Outcome of the 2006GRS－S Data－Fusion Contest［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45（10）：3012－3021.

［24］ MUNECHIKA C K，WARNICK J S，SALVAGGIO C，et al.Resolution Enhancement of Multispectral Image Data to Improve Classification Accuracy［J］.Photogrammetric Engineering and Remote Sensing，1993，59（Compendex）：67－72.

［25］ BRO R，JONG S.A Fast Non－Negativity－Constrained Least Squares Algorithm［J］.Journal of Chemometrics，1997，11：393－401.

［26］ WALD L.Data Fusion：Definitions and Architectures Fusion of Images of Different Spatial Resolutions［M］.Paris：Presses de l′Ecole，Ecole des Mines de Paris，2002.

［27］ WANG Zhou，BOVIK A C，SHEIKH H R，et al.Image Quality Assessment：from Error Visibility to Structural Similarity［J］.IEEE Transactions on Image Processing，2004，13（4）：600－612.