基于局部學習的遙感圖像融合

高東生（南京航空航天大學，江蘇 南京 210016）王 穎（中國科學院自動化所，北京 100086）

基于局部學習的遙感圖像融合

高東生（南京航空航天大學，江蘇 南京 210016）
王 穎（中國科學院自動化所，北京 100086）

本文提出了一種基于局部學習的遙感圖像融合方法。其基本思想是在局部區域對融合圖像與全色圖像建立對應的局部線性關系。由于圖像數據在局部區域相對簡單，因此局部模型相比全局模型更為合理。在局部學習的基礎上，將全色圖像與融合圖像的全局回歸誤差表示為圖拉普拉斯的形式，其本質是利用局部學習使得融合圖像保持全色圖像的流形結構。同時為了保持多光譜圖像的性質，通過圖像的尺度空間表示，建立融合圖像與多光譜圖像之間的尺度關系。最后通過集成融合圖像的二次拉普拉斯形式和尺度空間表示，構建圖像融合的全局目標函數。為了優化目標函數，本文提出了閉合求解法和快速迭代求解法。實驗結果表明：本文所提出的融合方法比傳統融合方法具有更好的效果。

局部學習；遙感圖像融合

1 簡介

在遙感圖像處理領域，對于光學傳感器,圖像的空間分辨率和光譜分辨率是一對矛盾的因素。要獲得高空間分辨率的圖像就只能以單光譜工作，而要獲得更多光譜就必須降低空間分辨率。為了同時提高圖像的空間分辨率和光譜分辨率，人們提出了利用融合全色圖像（高空間分辨率單光譜圖像）和多光譜圖像（低空間分辨率）來得到高空間分辨率多光譜圖像的方法。總結已有的圖像融合方法,可以將其分為三類:

第一類是基于顏色空間的轉換方法[1,2,3]。這些方法先將多光譜從原顏色空間轉換到另外一種顏色空間，這樣可以將亮度信息分量和顏色信息分量分開，然后用單光譜的全色圖像替換亮度信息分量，由此得到替換后的顏色空間表示，最后再對新的顏色空間進行逆變換,就得到融合后的圖像。如IHS(Intensity, Hue, Saturation)變換法先將原多光譜信息看作是RGB空間的顏色圖像，對它進行IHS變換分別得到亮度通道I、色調通道H和飽和度通道S。之后用全色圖像替代亮度通道I，再進行IHS逆變換從IHS空間變換到RGB空間，從而得到融合后的圖像。基于顏色空間轉換的方法實現原理簡單，但是這種基于顏色空間轉換的方法存在著兩個方面的問題：(1)一般只適應于三通道圖像融合；(2)雖然用全色圖像直接替代了亮度通道，但是色調通道和飽和通道仍是原多光譜的信息的簡單上采樣，這必然會導致空間細節上的損失。

第二類是基于成分分析的方法[3,4,5]這類方法與顏色空間轉換類似，不同的是它可以適用于任何波段的圖像融合，這類方法先將多光譜的每個像素看成是一個多維向量，對其進行統計成分分析(如主成分分析)，得到投影向量,再將多光譜圖像在這些投影向量上進行投影，在主成分投影的值保留了圖像大部分的信息。然后用全色圖像替換主成分圖像，之后進行相應的反變換就可以得到融合圖像。基于統計的主成分方法雖然可以進行任意波段圖像的融合，但融合后的圖像很難保持原多光譜圖像的顏色信息（即發生顏色失真），而且在空間細節上與原全色圖像相比也存在一定的損失。

第三類方法是基于小波分解法[6,7]，這類方法的基本思想是利用小波變換提取全色圖像的高頻信息（細節信息）和多光譜圖像的低頻信息（近似信息）。將通過小波變換提取的全色圖像的高頻信息和多光譜圖像的低頻信息組合成一組新的小波系數，然后對這些小波系數進行小波逆變換，從而得到融合圖像。但是基于小波變換的圖像融合當分解層數太少時，融合結果會存在明顯的造痕跡，從而降低了融合圖像的質量。當分解層數過多則會損失掉多光譜圖像的顏色信息。

近年來，基于梯度場的圖像編輯方法[8,9]在遙感圖像融合方面得到了成功的應用。其基本思想是讓融合圖像保持全色圖像的梯度，同時將低分辨率的多光譜圖像看作是邊界條件，并通過求解Poisson方程來得到融合結果。但是基于梯度場融合的方法由于其邊界條件過于嚴格，使得融合圖像對于多光譜圖像的光譜信息存在一定損失，同時通過梯度算子保持細節信息過于簡單，難以保持全色圖像的所有細節信息。

綜上所述，雖然現有的全色圖像與多光譜圖像融合方法很多，但這些融合方法很難在空間分辨率和光譜分辨率上同時達到最優。為了盡可能使得融合圖像在空間分辨率與光譜分辨率上達到最優。本文在局部學習的基礎上，對融合圖像與全色圖像的局部區域進行建模。具體而言，假設融合圖像與全色圖像在局部區域存在線性關系，對于每個局部區域，將回歸損失函數表示為融合圖像局部像素的二次拉普拉斯形式。通過對局部區域的損失函數進行累加，構建關于融合圖像的全局圖拉普拉斯。同時將多光譜圖像看作是待融合圖像的低分辨率表示。最后，構建總體目標函數。通過一定的轉化，目標函數可以通過求解稀疏線性方程組得到閉合解。為了適應大尺寸圖像融合，本文還提出了一種快速迭代的優化方法來求解融合圖像。實驗結果表明，本文所提出的融合方法相比于已有的算法能得到更高質量的融合圖像。

2 基于局部學習的圖像融合

設p、m1和y分別表示高分辨的全色圖像，低分辨率的多光譜圖像以及待融合圖像。如圖1所示，對于任意像素i所在的局部窗口Ni，假設全色圖像p與待融合圖像y在這個局部窗口滿足：

其中yi中表示融合圖像y在像素點j的像素值，pj表示全色圖像p在像素點j的像素值，wi和bi表示局部線性模型的參數。

圖1 局部線性模型示意圖

在局部區域內，為了使得融合圖像與全色圖像滿足上述線性模型，可以最小化如下局部窗口Ni內的正則化損失函數：

上式的第二項為正則化項，避免過擬合，其中，η表示正則化參數。對于整幅圖像而言，對局部的損失能量進行累加，可得到全局損失函數為

最小化上式能量，使得融合圖像與全色圖像在所有局部區域能盡量滿足局部線性關系。

對于多光譜圖像m,我們將其看作是待融合y的低分辨率表示，即：

其中g_σ表示方差為σ的高斯函數，*表示濾波操作。設全色圖像的分辨率是多光譜圖像的n倍，根據尺度空間理論[10],σ可以根據以下公式來計算來確定：

最小化式(5)可以使得融合圖像盡量滿足等式(4)。

因此，對于圖像融合來說，其目標為最小化(3)和(5)的整體能量。下一節我們將構建圖像融合的總體目標函數，并提出兩種不同的求解方法對目標函數進行優化。

3 目標函數優化

結合式(3)和(5)，圖像融合的目標為最小化如下所示的能量函數：

其中，λ表示權重因子。接下來的兩小節將分別給出兩種不同的方法來最小化上式目標函數(6)。

3.1 閉合求解

為了求解融合圖像y，E(w,b,y)對w,b求導，并令其為0，可求得wi和bi分別為：

式中，Ni表示局部窗口Ni內像素的個數，為全色圖像p在局部窗口Ni內的像素所組成的向量為待融合圖像y在局部窗口Ni內的像素所組成的向量。設(I為單位矩陣)為Ni×Ni的中心化矩陣,(即 pi為pi中心化后的向量)。將式(7)和(8)代入式(6)，那么公式(6)的第一項可以寫為如下形式：

再考慮式(6)的第二項，由于濾波操作為線性操作，可以表示為稀疏矩陣-向量乘法的形式[13]，這里假設高斯濾波核gσ對應的稀疏矩陣表示形式為G。通過消除變量w,b后，圖像融合的目標函數為：

式(10)的優化問題是關于y的二次型，那么其全局最優解可以對y求導并令其為0得到，即等價于求解以下稀疏線性方程：

(L+λG^TG)y=λG^Tm.

3.2 快速迭代求解

采用上述方法消除中間變量w,b，使得目標函數的優化等價于求解線性方程組。這種求解方式可以得到精確的閉合解。但是這種求解方法的計算代價過大。假設融合圖像的大小為M×N，局部窗口的大小為ω×ω，那么構造拉普拉斯矩陣L的計算復雜度為：O(MNω6)，當圖像尺寸較大時，構造矩陣L的復雜度過高。同時，雖然矩陣L為稀疏矩陣，但其大小為MN×MN，且每行的非零元素個數為(2ω-1)2。因此，當圖像尺寸較大時，存儲矩陣L這樣一個大規模的稀疏矩陣所占用的內存仍然很大。

這說明第3.1小節的閉合求解方法存在兩個方面的問題。一方面由于構建拉普拉斯矩陣的復雜度較高，使得求解速度仍然較慢，另一方面，當圖像比較大時，稀疏矩陣L需要占用大量的內存資源。因此，閉合求解的方法不適應于大幅面圖像的融合。為此，我們提出一種快速的迭代求解方法。

首先，我們可以將方程(6)寫為如下積分形式：

式中，k(i,j)表示一個窗口函數，假設窗口大小為3×3，那么k(i,j)的取值是：當像素j在3×3窗口內部時k(i,j)=1，否者為0。為了后續計算方便，將k(i,j)的值進行歸一化，使其求和等于1，即k(i,j)=1/9。

可以采用類似EM的交替迭代方法，來最小化目標函數(11)。

（1）固定y，求解w,b。

利用變分法，式(11)對w,b求導并令其為0，可得：

式中，由于濾波核k的元素都相等，因此，上述的濾波操作可以采用積分圖像的技巧[14]進行加速。

（2）固定w,b，求解y。

式(11)對y求導，并令其為0，可求得：

直接求解方程(13)比較困難，這里我們將其轉換到頻率上進行求解。方程(13)兩邊同時進行快速傅立葉變換可得：

式中，Y,P,M,W,B分別表示y,p,m,w,b的傅立葉變換，G,K分別表示濾波核gσ,k的傅立葉變換。通過式(14)求得Y為：

對Y進行傅立葉逆變換就可求得y。

重復上述兩個步驟直到算法收斂。由于交替迭代的優化方法只能保證局部最優，因此其最終解依賴于初始值y的給定，對于全色圖像與多光譜圖像融合來說，多光譜圖像m是y的一個非常好的初始解，因此y的初值設定為多光譜圖像m。

4 實驗結果及分析

為了驗證本文融合算法的有效性，本節將進行高分辨率全色圖像和低分辨率多光譜的融合實驗。在本節實驗中，本文所提出的融合算法的參數設置為：局部窗口大小為33，局部線性回歸的正則化參數η=10-4。首先，我們與Poisson融合方法做比較，圖2和3表示在QuickBird與LandsatETM+衛星圖像上的融合結果。圖2為2.4米多光譜與0.6米全色圖像融合的結果（為了體現融合的對比效果,對多光譜圖像進行了一次下采樣），圖3為30米多光譜與15米全色圖像的融合結果。從圖2（注意紅色矩形框所標注的圖像）和圖3可以看出本文所提出的圖像融合方法在光譜保持和空間細節上要明顯優于Poisson融合方法。

圖2 與Poisson融合的比較

圖3 與Poisson融合的比較

為了定量評價本文所提出融合方法的性能，本文通過對融合圖像與原始全色圖像和多光譜圖像進行比較來評價融合方法的有效性。本文從兩個方面來評價融合圖像：圖像空間細節的保持度和光譜信息保持度。空間細節的保持度通過計算融合圖像與全色圖像的高頻信息的差異來度量，即：

其中，c表示第c個波段，C為波段數，M,N為圖像尺寸，y'c,p'分別表示融合圖像和全色圖像的拉普拉斯響應。光譜信息的保持度通過計算融合圖像的低分辨率圖像與多光譜圖像之間的差異來度量，即：

圖4 本文融合方法與Poisson

圖4表示本文所提出的融合方法和Poisson融合在圖2和3所計算得到的平均E_spatial和E_spectral。從圖4可以看出本文所提出的方法在空間細節和光譜保持度上都要明顯由于Poisson融合，而快速迭代求解算法與閉合求解方法的差異非常小，這說明迭代求解算法在快速求解的同時，融合效果與閉合求解方法基本一致。

除了Poisson融合相比較外，本文還與一些傳統的融合方法進行了比較，包括：IHS、PCA、Brovey變換、小波變換（這里采用db4小波進行4層分解，并將多光譜圖像小波變換的近似信息與全色圖像小波變換的細節信息進行組合，再進行小波逆變換）。

為了客觀評價融合方法的有效性，我們采用文獻[15]中所提出的方法，先對原始多光譜圖像和全色圖像都進行尺度退化，并將原多光譜圖像作為參考圖像。通過計算融合圖像與參考圖像之間的均方根誤差(root-mean-squareerror,RMSE),相關性系數(correlationcoefficient,CC),相對均值偏差(AbsoluteMeanofBias,AMB)來評價融合方法的性能。圖5為各種融合方法在2米多光譜與0.6米全色圖像上的融合結果。表1為各種融合方法在圖5所得結果的平均指標。從表1的各種指標可以看出，本文所提出的融合方法具有最小的RMSE，最大的CC以及最小的AMB，這說明本文方法相比于其他方法具有更好的融合效果。同時注意到迭代求解的方法與閉合求解方法的性能指標非常近似，這說明迭代求解的方法不僅與閉合求解法相當的融合效果。

表1 圖5中各種融合方法的RMSE, CC和Bias

5 總結

本文提出了一種基于局部學習的遙感圖像融合方法。該方法一方面利用局部學習方法使融合圖像保持全色圖像的流形結構。另一方面將多光譜圖像看作是融合圖像的低分辨率表示。結合這兩方面構建圖像融合的目標函數，目標函數的優化可以通過求解稀疏線性方程組來得到閉合解，同時為了提高求解效率，本文還提出了一種快速迭代的求解方法。實驗結果表明，本文所提出的融合方法相比傳統融合方法具有更好的效果。同時，本文所提出的快速迭代算法將多光譜圖像作為初始解，只需要少量的迭代次數（實驗中一般小于10次）就達到穩定解。特別地，對于融合一幅1024×1024的圖像而言，快速迭代求解方法一般只需要1秒左右的時間就可以得到與閉合求解（10秒）相當的融合結果。

圖5 各種融合方法在QuickBird圖像上的結果

[1] C. PohlandJ. V. Genderen. Multisensorimagefusioninremotesensing: Concepts, methodsandapplications[J]. InternationalJournalofRemoteSensing, 1998 (19) 5, 823-854

[2] Z. Wang, D. Ziou, C. Armenakis, D. Li, Q. Li.Acomparativeanalysisofimagef usionmethods[J]. IEEETransactionsonGeoscienceandRemoteSensing, 2005 (43) 6, 1991-1402

[3] F. A. Al-Wassai, N. Kalyankar, andA. A. Al-Zuky. TheIHStransformationsba sedimagefusion[J]. JournalofGlobalResearchinComputerScience, 2011 (2) 5.

[4] G. Simone, A. Farina, F. C. Morabito, S. B. Serpico, andL. Bruzzone. Imagefusi ontechniquesforremotesensingapplications[J]. InformationFusion, 2002 (3) 1, 3-15. [5] J.Sun, Y. Jiang, S. Zeng. AstudyofPCAimagefusiontechniquesonremotesensi ng. inProceedingsoftheSPIEInternationalConferenceonSpaceInformationTechno logy. 2005.

[6] H. Li, B. Manjunath, S. Mitra. Multi-sensorimagefusionusingthewavelettrans form. inIEEEInternationalConferenceImageProcessing.

[7] 晁銳, 張科, 李言俊. 一種基于小波變換的圖像融合算法[J]. 電子學報, 2002 (32) 5, 750-753.

[8] J. Wen, Y. Li, H. Gong.Remotesensingimagefusionongradientfield.inProceedi ngsofthe18thInternationalConferenceonPatternRecognition (ICPR04) .

[9] Z. Zhou, S. Peng, B. Wang, Z. Hao, andS. Chen.Anoptimizedapproachforpan sharpeningveryhighresolutionmultispectralimages. IEEEGeoscienceandRemote SensingLetters.

[10] T. Lindeberg. Scale-SpaceTheoryinComputerVision. Norwell. MA, USA: KluwerAcademicPublishers, 1994.

[11] S. Xiang, F. Nie, C. Pan, C. Zhang.RegressionreformulationsofLLEan dITSAwithlocallylineartransformation. IEEETransactionsonSystems, Man, andCybernetics, PartB: Cybernetics, vol. 41, no. 5, pp. 1250–1262, 2011.

[12] A. Levin, D. Lischinski, Y. Weiss. Aclosed-formsolutiontonaturalimagematt ing. IEEETransactionsonPatternAnalysisandMachineIntelligence.

[13] Y. Wang, H. Yan, C. Pan, S. Xiang. Imageeditingbasedonsparsematrixvectormultiplication. in2011IEEEInternationalConferenceonAcoustics, SpeechandSignalProcessing (ICASSP) .

[14] F. C. Crow. Summed-areatablesfortexturemapping. inProceedingsofthe11th AnnualConferenceonComputerGraphicsandInteractiveTechniques.

[15] L. Wald, T. Ranchin, M. Mangolini. Fusionofsatelliteimagesofdifferentspati alresolutions: Assessingthequalityofresultingimages. Photogrammetricengineeri ngandremotesensing, 1997.

高東生，南京航空航天大學信息工程學院電力與通信工程專業碩士研究生，主要研究方向為遙感圖像處理。

王穎，中國科學院自動化研究所助理研究員，主要研究方向為遙感圖像處理、模式識別和計算機視覺。

更正聲明

刊登在《自動化博覽》雜志2013年12月刊第77頁的張萍的作者簡介中“甘肅寧夏人”有誤，應為“寧夏銀川人”。特此聲明。

Remote Sensing Image Fusion Based on Local Learning

This paper presents a local learning based method for remote sensing image fusion. The key idea of our method is to construct the local linear model between the panchromatic image and the fused image in local region. Compared with the complexity in global region, image data in local region are much simpler. This means that local model is more reasonable than global model. Based on local learning, the regression error between the panchromatic image and the fused image can be formulated with Laplacian quadratic form, which makes the fused image preserve the manifold structure of the panchromatic image. In order to preserve the properties of the multispectral image, we build the scale space relationship between the fused image and the multispectral image via scale space representation. Consequently, the objective function of our model is proposed by combining the Laplacian quadratic form and the scale space representation. Meanwhile, a close form solution method and a fast iterative method are proposed to optimize the proposed model. Experimental results demonstrate the effectiveness of our model in comparison to the state-of-the-art methods.

Local Learning; Remote Sensing Image Fusion

B

1003-0492(2014)01-0086-05

TP212