基于熱傳導模型的像素級遙感圖像融合

許君一,盧秀山,卿熙宏,姚繼鋒

1.山東科技大學測繪科學與工程學院,山東青島266510;2.山東省3S工程技術研究中心,山東青島266510;3.海島(礁)測繪技術國家測繪局重點試驗室,山東青島266510;4.山東科技大學資產處,山東青島266510

基于熱傳導模型的像素級遙感圖像融合

許君一1,2,3,盧秀山1,3,卿熙宏1,4,姚繼鋒1

1.山東科技大學測繪科學與工程學院,山東青島266510;2.山東省3S工程技術研究中心,山東青島266510;3.海島(礁)測繪技術國家測繪局重點試驗室,山東青島266510;4.山東科技大學資產處,山東青島266510

通過熱傳導方程給出一種像素級遙感圖像融合模型和方法:①給出空間域內高分辨率圖像與低分辨率圖像之間的擴散關系,作為特例得到Brovey變換(brovey transform,BT);②給出圖像融合與增強的統一表達式并得到基于亮度平衡的融合方法;③低分辨率多光譜圖像的方差較小情形,指出基于方差的標準圖像融合方法將會丟失高空間分辨率全色圖像信息。試驗表明,除了圖像量化誤差以外,所提議的方法不會丟失已知圖像的空間分辨率和波譜信息。

圖像融合;遙感;HSI;多光譜圖像;熱傳導方程

1 引 言

遙感圖像處理中,高空間分辨率全色圖像與低空間分辨率多光譜圖像的融合技術一直是其研究重點之一。像素級融合方法大體上可分為顏色模型法(HSI、HSB、HSV等)、變換域法(離散小波變換DWT、主分量法 PCA等)、參數估計法[1]。顏色模型法中,HSI變換因其亮度空間與彩色空間的正交性而受到重視。參數估計法中, Bayes法也有比較滿意的結果[2]。文獻[2]通過引入多個觀測模型,增強了待估計圖像各個主成分的空間分辨率,并獲得了線性最小均方誤差意義上的高分辨率多光譜圖像Bayes估計。變換域融合方法中最常見的DWT法以適當損失空間分辨率為代價較好地保留了波譜信息[3-4]。也有作者通過DCT變換進行圖像融合[5]。有一些作者從波譜保持性角度研究多光譜遙感圖像融合問題[6]。從應用角度分析,各種圖像融合方法各具特點[7]。方法研究方面,人們不斷探索新的圖像融合方法使其更具有普遍性和更好的效果。如,偏微分方程理論正逐漸受到重視[8-9]。

現代圖像處理中,非線性擴散方程在圖像的復原、邊緣檢測、分割及去噪等問題中有著重要的應用[7-9]。自從1987年Malik和 Perona把非線性傳播(擴散)方程引入到圖像處理后,人們大量研究了具有對比度不變性、仿射不變性的非線性擴散方程[7-8]。D.A.Socolinsky在他的博士論文中通過偏微分方程研究了圖像融合問題[10]。D.A.Socolinsky的研究為圖像融合問題提供了新的思路和方法[11-12]。

熱傳導方程在研究許多擴散問題時起到借鑒作用。用熱傳導方程來描述圖像的多尺度空間分辨率和圖像的復原問題已經有豐富的成果[6-7]。筆者通過擴散原理和偏微分方程理論研究了低空間分辨率多光譜波段圖像與高空間分辨率全色波段圖像的融合問題。試驗表明,新方法的理論結果是穩定和正確的。

2 基于擴散模型的多光譜和全色波段圖像融合

熱擴散過程與圖像的多尺度分辨率之間有許多相似或相同之處。下面通過熱擴散原理建立多光譜和全色波段圖像融合的熱傳導模型及算法。

2.1 熱傳導擴散模型

設高空間分辨率多光譜圖像在0≤t≤T內發生擴散且滿足熱傳導方程。多光譜圖像的多尺度熱傳導方程與一般熱傳導方程約束條件有所不同的是,模型中增加了終止條件。為了融合后的圖像亮度與多光譜圖像亮度保持一致,給出了均值約束條件。作為可選項,給出了表示圖像細節反差的方差約束作為候選條件。根據以上假設,多尺度圖像熱傳導方程模型為

其中,k=1,2,3,4分別代表R、G、B三個波段及圖像亮度;φ(x)表示高空間分辨率全色圖像。作為擴散條件,初始圖像φk(x)表示高分辨率圖像(若無特別聲明以下分辨率均指空間分辨率);終止圖像ψk(x)表示低分辨率圖像;uk(x,t)為 t尺度圖像;E(＊)表示圖像平均值;D(＊)表示圖像方差。假定ψk(x)是φk(x)經過 T時間尺度擴散的結果。所以,φk(x)的分辨率高于ψk(x)。對于多光譜k=1,2,3情況,式(1)構成求解初始條件φk(x)和熱源函數 fk(x,t)的反問題。理論上,式(1)的求解是困難的。

圖像亮度一般是各個波段灰度的平均值

所以,多光譜圖像亮度 u4(x,t)在也滿足式(1)的形式。式(1)的Fourier變換為(ξ=(ξ1,ξ2,ξ3)T)

求解式(3)得(k=1,2,3,4)

式(3)中要求滿足條件^uk(ξ,T)=^ψk(ξ),即

令χ(τ)為特征函數

則

做關于t的Fourier變換得

即做逆Fourier變換并求函數^fk(ξ,t)得

其中,δ(t)為Dirac函數,t∈[0,T]。把上式代入下式

得

整理和考慮t=0情況,得

2.2 多光譜與全色波段圖像融合

根據式(4),有

即

復數模|＊|代表了圖像的能量。如果把像素值作為該點的能量,復數模|＊|代表了構成該像素的更小尺度圖像的能量。如圖1所示,數字圖像的獲得可看成如下過程:大圖像是由子圖像拼接得到的。

圖1 大圖像由子圖像拼接得到Fig.1 Sub-image to be large mosaic image

把子圖像縮減成一個光柵點構成新的數字圖像

根據式(6),式(5)表示低分辨率多光譜子圖像、高分辨率全色子圖像及高分辨率多光譜子圖像之間的能量關系。這種能量關系,對大的圖像來說就是像素值之間的關系。通過式(5)、式(6),從子圖像頻率域轉換到大圖像空間域

根據式(7),可通過全色圖像φ4和低分辨率多光譜圖像獲得高分辨率的多光譜圖像(φ1,φ2,φ3)。理想情況下,高分辨率多光譜圖像的亮度φ4和高分辨率全色圖像φ完全一致。但從信息獲取角度看,φ和φ4的區別是φ為已知量,φ4為未知量。用已知的φ代替未知的φ4得到

若假定 ak=a(k=1,2,3),t∈[0,T],根據式(2),自然有 a4=a。此時直接得到 Brovey變換(Brovey transform,BT)[13]。在熱擴散過程中其總能量守恒。所以,擴散過程中可以認定多光譜圖像的各個波段平均亮度不變,E{ψk}=E{φk}(k =1,2,3)。但是,遙感圖像中常見情況是低分辨率的多光譜圖像的亮度與全色圖像的亮度不一致。因此,對式(8)進行如下調整

其中,k=1,2,3。式(11)將使得到的高分辨率多光譜圖像與低分辨率多光譜圖像的亮度保持一致。當ak=a(k=1,2,3),t∈[0,T]時,由式(9)進一步得

經過式(10)給出的高分辨率多光譜圖像的亮度方差為(用D(＊)表示方差)

同理,通過式(8)也能得到常用的標準圖像融合算法

值得注意的是,式(12)的方差 D(IStd)=D(ψ4)而不是D(IStd)=D(φ)。所以,式(12)的分辨率不能很好地代表高分辨率圖像的分辨率。另外,當D(ψ4)≈0時出現 D(IStd)≈0,即得到亮度為E{ψ4(x)}的無圖像細節的同一色圖像。這與事實是不符合的。實際上,多光譜圖像方差為0未必全色圖像沒有細節。在數字圖像中均值為0的圖像方差一定為0,但反過來未必成立。從這個角度,式(7)～式(9)和式(9)～式(10)更具有普遍意義。從實際應用角度,式(12)也不是完全不可用。當低分辨率光譜圖像亮度方差D(ψ4)比較大時,式(12)的效果不僅很好而且表現出對波譜信息的高分辨率特性。所以,更為實際的方案是要依據所獲得的圖像和應用要求決定采用哪一個。

3 試驗與分析

以下試驗原始影像均為256×256的低分辨率多光譜圖像和高分辨率全色圖像。試驗中,除了特別指出參數

以外,都默認為λk=1(k=1,2,3)。稱參數λk為反擴散系數。

試驗1 圖2為原始影像。圖3(a)、(b)通過式(8)、式(9)得到。目視效果上,圖3(a)亮度與圖2(b)一致,圖3(b)亮度與圖2(a)一致。這一點同理論分析一致。圖4通過式(12)得到。細節反差上,圖4不如圖3(b)。圖4整體亮度相對較低,在高亮度區可分辨部分細節。

圖2 低分辨率多光譜圖像(a)和高分辨率全色圖像(b)Fig.2 The low-resolution multispectral images(a)and the high-resolution panchromatic image(b)

圖3 (a)、(b)分別為通過式(8)和式(9)得的結果Fig.3 (a)and(b)are obtained by method(8)and (9),respectively

圖4 根據式(12)得的結果Fig.4 The result obtained by method(12)

圖5 原圖像Fig.5 The original image

試驗2 圖5、圖6為原始影像。圖5為高分辨率多光譜圖像。圖7(a)直接通過式(9)得到,圖7(b)通過亮度與彩色空間正交的HSI變換后,對亮度分量用式(12)處理,然后進行HSI反變換得到。表1為均方誤差

其中,(si,j)為原始高分辨率多光譜圖像,(f usioni,j)為融合后的圖像。圖8(a)、(b)通過式(9)得到。其中反擴散系數λk=0.2,1.8(k=1,2,3)。

圖6 低分辨率多光譜圖像和高分辨率全色圖像Fig.6 The low-resolution multispectral images(a)and the high-resolution panchromatic image(b)

圖7 (a)、(b)分別為通過式(9)和式(12)得的結果Fig.7 (a)and(b)are obtained by method(9)and (12),respectively

圖8 (a)、(b)經式(9)得到,λk=0.2,1.8(k=1,2,3)Fig.8 (a)and(b)are obtained by method(9)using parametersλk=0.2,1.8(k=1,2,3),respectively

試驗3 原始影像如圖9所示。圖10(a), (b)通過式(10),式(12)得到。圖10(b)沒有能夠體現出高分辨率全色圖像的信息。這與理論分析一致。

圖9 低分辨率多光譜圖像和高分辨率全色圖像Fig.9 The low-resolution multispectral images(a)and the high-resolution panchromatic image(b)

圖10 (a)、(b)分別為通過式(9)和式(12)得到的結果Fig.10 (a)and(b)are obtained by method(9)and (12),respectively

表1 試驗2的RGB波段均方誤差Tab.1 Themeansquared errorofRGB bandin second experiment

試驗4 采用試驗1的圖像數據進行基于4層小波分解的兩種融合算法試驗。

融合方法1 ①對圖2(a)進行HSI變換,獲得亮度分量I;②對I和圖2(b)進行4層小波分解,小波采用db4(daubechies 4);③在第4層,低頻取I的低頻分量,高頻取I和圖2(b)相互對應的高頻分量中系數絕對值大的作為新的高頻分量,并把這些高低頻分量作為圖2(b)的第4層小波分解分量;④對圖2(b)進行逐層小波重構,其結果作為新多光譜圖像的HSI亮度分量I;⑤進行HSI反變換,得到圖11(a)。

融合方法2 ①進行融合方法1中的過程①和②;②在第4層,通過I和圖2(b)的第4層低通分量進行主成分分析(PCA)得到權值;③對I和圖2(b)的第4層分量進行加權求和,并把它作為圖2(b)的第4層分量;④對圖2(b)進行逐層小波重構和HSI反變換,得到圖11(b)。試驗中,對應全色和多光譜圖像的權值分別為0.376 72和0.623 28。

圖11 圖像(a)是試驗4融合方法1的結果,而圖像(b)是融合方法2的結果Fig.11 (a)and(b)are the results of fusion method 1 and 2 in experiment 4,respectively

分析 試驗1中,整體上亮度約束效果比方差好。在高亮度區域,由于低分辨率多光譜圖像高亮度區的過高亮度和融合圖像的量化截斷,導致圖3(b)在高亮度區域難于分辨空間細節。基于方差的融合圖像圖4在高亮度區體現了一定的細節。試驗2給出了融合圖像和原始圖像的差別程度。目視效果圖7(a)比圖7(b)好一些,但差別不明顯。表1數據表明,各波段的平均離散程度都很小。由式(9)和式(10)得到圖像的均方誤差完全一樣。這也從另一個方面說明了亮度空間與彩色空間正交的 HSI變換的合理性。圖8(a)、(b)表明,融合中擴散條件的不同而表現出圖像的逐漸增強。不同的參數代表了不同的擴散程度,所以當加大逆擴散程度時表現出了使圖像增強的特點。試驗1～試驗3表明了理論與試驗的一致性,同時指出了用亮度作為約束條件比用方差作為約束條件更有效和合理。試驗3直接指出了方差約束在極端情形下的不合理性。試驗4給出了兩種典型的小波融合試驗。目視效果上,圖11(a)的波譜信息與圖2(a)比較接近。在細節上和清晰程度上,圖3(b)的空間分辨率明顯好于圖11(a)。雖然圖11(b)在細節保持性上有一定的改進,但較多地丟失了波譜信息。所以,無論是波譜信息保持性上還是空間分辨率上,圖11(b)都不如圖3(b)。目前,小波融合算法的多數融合策略是取多光譜的低頻和全色波段的高頻作為融合策略,其目的是保持多光譜的波譜信息和全色波段的高空間分辨率信息。研究表明,小波融合算法中的分解層數不是越多或越少就越好。因此,許多小波圖像融合算法對空間分辨率和波譜信息損失是不可避免的。如何使空間分辨率和波譜信息的損失減少到最少是小波圖像融合技術中的一個核心問題。試驗4的結果也說明了這個問題。在空間分辨率和波譜信息損失中,不僅有目視效果上可感覺得到的損失,也有感覺不到的損失。這些損失是小波分解下的圖像融合策略導致的。這一點上,新方法具有空間分辨率和波譜信息保持性上的優勢。

4 結 論

筆者根據熱傳導方程得到了圖像融合算法。根據熱傳導方程,理論上得到了BT算法和融合與增強統一描述的更一般的算法。試驗表明,該算法在波譜保持性和細節保持性上表現出比較好的特性。分析和試驗表明,標準圖像融合算法在方差比較小時丟失高分辨率全色圖像信息。從多光譜波譜信息上,除了調整波譜亮度以外,新方法在不同波譜之間的相對比值關系沒有發生變化。所以,從波譜之間的比值角度來說,融合結果沒有發生波譜畸變。新方法適合應用于遙感圖像目視解譯與分析。新方法不需要做 HSI變換,這將有效減少計算量。相對于Bayes法,新方法無需進行迭代計算,也沒有像Bayes法中出現的光譜信息和空間細節信息等方面的預測或估計問題。新方法是基于熱擴散理論的方法,理論上是嚴格而精確的算法,無須進行最優性估計。融合原理試驗結果表明,新方法在細節保持性上優于多層小波分解方法。

新算法精度在于,不同分辨率的遙感圖像融合是否符合擴散原理。作為新方法特例的BT方法文獻和本文試驗表明,多尺度意義上的擴散模型滿足熱傳導方程在內的擴散方程。所以,新方法不僅具有一定的普適性,更具有進一步進行理論發掘的潛力。多光譜融合圖像的波譜真偽與低分辨率的多光譜圖像的擴散程度有直接關系。如果擴散程度過大,其光譜不能代表高分辨率空間位置的波譜信息。因此,在適當的波譜空間分辨率下的融合才有實際意義。

根據以上分析和結論,所提議的方法具有較強的理論意義和較好的實際意義。根據本文研究結果和思路,將進一步研究多焦距圖像融合的擴散模型和非線性融合模型。

[1] ZHANG Yujin.Image Engineering[M].2nd ed.Beijing: Tsinghua University Press,2007:315-343.(章毓晉.圖像工程[M].2版.北京:清華大學出版社,2007:315-343.)

[2] GE Zhirong,WANG Bin,ZHANG Liming.Remote Sensing Image Fusion Based on Bayesian Linear Estimation [J].Science in China:Series F,Information Sciences, 2007,50(2):227-240.(葛志榮,王斌,張立明.基于Bayes線性估計的遙感圖像融合[J].中國科學:F輯,信息科學,2007,37(4):501-513.)

[3] LIU Bin,PENGJiaxiong.Fusion Method of Multi-spectral Image and Panchromatic Image Based on Four Channels Non-separable Additive Wavelets[J].Chinese Journal of Computers,2009,32(2):350-356.(劉斌,彭嘉雄.基于四通道不可分加性小波的多光譜圖像融合[J].計算機學報,2009,32(2):350-356.)

[4] LIANG Dong,LI Yao,SHEN Min,et al.An Algorithm for Multi-focus Image Fusion Using WaveletBased Contourlet Transform[J].Acta Electronica Sinica,2007, 35(2):320-322.(梁棟,李瑤,沈敏,等.一種基于小波-Contourlet變換的多聚焦圖像融合算法[J].電子學報, 2007,35(2):320-322.)

[5] 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.)

[6] WU Lianxi,LIANG Bo,LIU Xiaomei,et al.A Spectral Preservation Fusion Technique for Remote Sensing Images [J].Acta Geodaetica et Cartographica Sinica,2005,34 (2):118-128.(吳連喜,梁波,劉曉梅,等.保持光譜信息的遙感圖像融合方法研究[J].測繪學報,2005,34(2): 118-128.)

[7] CHENG Yinglei,ZHAO Rongchun,LI Weihua.et al. Overview of Methods of Data Fusion for Images Based on Pixel-level[J].Application Research of Computers,2004, 4(5):169-172.(程英蕾,趙榮椿,李衛華,等.基于像素級的圖像融合方法研究[J].計算機應用研究,2004,4 (5):169-172.)

[8] SAPIRO G.Geometric Partial Differential Equations and Image Analysis[M].Cambridge:Cambridge University Press,2001:260-340.

[9] MAITRE H.Le Traitement des Images[M].Translated by SUN Hong.Beijing:Publishing House of Electronics Industry,2006:156-173.(MAITRE H.現代數字圖像處理[M].孫紅,譯.北京:電子工業出版社,2006: 156-173.)

[10] SOCOLINSKY D A.A Variational Aproach to Image Fusion[D].Baltimore:Johns Hopkins University,2000.

[11] SOCOLINSKY D A.Dynamic Range Constraints in Image Fusion and Visualization[C]∥Proceedings of Signal and Image Processing 2000.Las Vegas:[s.n.],2000: 349-354.

[12] SOCOLINSKY D A,WOLFF L B.A New Visualization Paradigm for Multispectral Imagery and Data Fusion[C]∥Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition:Vol 1.Fort Collins:IEEE,1999:319-324.

[13] POHL C.Tools and Methods for Fusion of Images of Different Spatial Resolution[C]∥International Archives ofPhotogrammetry and Remote Sensing:Vol32. Valladolid:[s.n.],1999:76-82.

(責任編輯:叢樹平)

Pixel-level Remote Sensing Images Fusion Based on Heat Conduction Model

XU J unyi1,2,3,LU Xiushan1,3,QING Xihong3,YAO Jifeng1
1.College of Geomatics,Shandong University of Science and Technology,Qingdao 266510,China;2.“3S”Research Center of Engineering and Technology,Qingdao 266510,China;3.KeyLaboratory of Surveying and Mapping Technology on Island and Reef, State Bureau of Surveying and Mapping,Qingdao 266510,China;4.Assets Division,Shandong University of Science and Technology, Qingdao 266510,China

A novel pixel-level remote sensing image fusion based on heat conduction equation is proposed.The main contributions are as follows:①the diffusion relationship between the high-resolution image and the low-resolution image are obtained in the space domain,where Brovey transform(BT)is one of its special cases;②a unified expression of pixel-level image fusion and image enhancement is obtained,and a brightness balanced-based image fusion is obtained by using the expression;③the disadvantage that the standard variance-based image fusion method will result in the loss of high spatial resolution panchromatic image information,for smaller variance of lowresolution multispectral images is pointed out.Experimental results show that this approach does not lose the image spatial resolution and spectral information and outperforms the existing method.

image fusion;remote sensing;HSI;multispectral image;heat conduction equation

XU J unyi(1960—),male,PhD,professor,majors in the image processing,pattern recognition,digital watermarking,digital photogrametry and marine charting.

E-mail:xjyxxl@yahoo.com.cn

1001-1595(2010)06-0566-06

P237

A

國家自然科學基金(30700107)

2009-08-01

2009-09-28

許君一(1960—),男,博士,教授,主要研究方向為圖像處理、模式識別、數字水印、數字攝影測量和海洋測繪。