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

基于局部學(xué)習(xí)的遙感圖像融合

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

本文提出了一種基于局部學(xué)習(xí)的遙感圖像融合方法。其基本思想是在局部區(qū)域?qū)θ诤蠄D像與全色圖像建立對(duì)應(yīng)的局部線性關(guān)系。由于圖像數(shù)據(jù)在局部區(qū)域相對(duì)簡單，因此局部模型相比全局模型更為合理。在局部學(xué)習(xí)的基礎(chǔ)上，將全色圖像與融合圖像的全局回歸誤差表示為圖拉普拉斯的形式，其本質(zhì)是利用局部學(xué)習(xí)使得融合圖像保持全色圖像的流形結(jié)構(gòu)。同時(shí)為了保持多光譜圖像的性質(zhì)，通過圖像的尺度空間表示，建立融合圖像與多光譜圖像之間的尺度關(guān)系。最后通過集成融合圖像的二次拉普拉斯形式和尺度空間表示，構(gòu)建圖像融合的全局目標(biāo)函數(shù)。為了優(yōu)化目標(biāo)函數(shù)，本文提出了閉合求解法和快速迭代求解法。實(shí)驗(yàn)結(jié)果表明：本文所提出的融合方法比傳統(tǒng)融合方法具有更好的效果。

局部學(xué)習(xí)；遙感圖像融合

1 簡介

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

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

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

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

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

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

2 基于局部學(xué)習(xí)的圖像融合

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

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

圖1 局部線性模型示意圖

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

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

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

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

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

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

因此，對(duì)于圖像融合來說，其目標(biāo)為最小化(3)和(5)的整體能量。下一節(jié)我們將構(gòu)建圖像融合的總體目標(biāo)函數(shù)，并提出兩種不同的求解方法對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化。

3 目標(biāo)函數(shù)優(yōu)化

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

其中，λ表示權(quán)重因子。接下來的兩小節(jié)將分別給出兩種不同的方法來最小化上式目標(biāo)函數(shù)(6)。

3.1 閉合求解

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

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

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

式(10)的優(yōu)化問題是關(guān)于y的二次型，那么其全局最優(yōu)解可以對(duì)y求導(dǎo)并令其為0得到，即等價(jià)于求解以下稀疏線性方程：

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

3.2 快速迭代求解

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

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

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

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

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

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

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

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

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

式(11)對(duì)y求導(dǎo)，并令其為0，可求得：

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

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

對(duì)Y進(jìn)行傅立葉逆變換就可求得y。

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

4 實(shí)驗(yàn)結(jié)果及分析

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

圖2 與Poisson融合的比較

圖3 與Poisson融合的比較

為了定量評(píng)價(jià)本文所提出融合方法的性能，本文通過對(duì)融合圖像與原始全色圖像和多光譜圖像進(jìn)行比較來評(píng)價(jià)融合方法的有效性。本文從兩個(gè)方面來評(píng)價(jià)融合圖像：圖像空間細(xì)節(jié)的保持度和光譜信息保持度?？臻g細(xì)節(jié)的保持度通過計(jì)算融合圖像與全色圖像的高頻信息的差異來度量，即：

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

圖4 本文融合方法與Poisson

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

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

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

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

5 總結(jié)

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

圖5 各種融合方法在QuickBird圖像上的結(jié)果

[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]. 電子學(xué)報(bào), 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.

高東生，南京航空航天大學(xué)信息工程學(xué)院電力與通信工程專業(yè)碩士研究生，主要研究方向?yàn)檫b感圖像處理。

王穎，中國科學(xué)院自動(dòng)化研究所助理研究員，主要研究方向?yàn)檫b感圖像處理、模式識(shí)別和計(jì)算機(jī)視覺。

更正聲明

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

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