朱衛(wèi)東,李全海,徐克科,李天子

(1.同濟(jì)大學(xué) 測量與國土信息工程系,上海200092;2.河南理工大學(xué) 測繪與國土信息工程學(xué)院,河南 焦作454003)

遙感圖像融合是將同一地區(qū)的不同類型的遙感圖像進(jìn)行有機(jī)變換結(jié)合的技術(shù).融合后得到的新圖像具有更高的可信度,模糊較少,可理解性更好,更適合于人的視覺或計(jì)算機(jī)視覺的檢測、分類、識別等后處理.遙感圖像融合技術(shù)在土地資源調(diào)查、洪水監(jiān)測、地形測繪以及地物分類等方面有著廣泛的應(yīng)用.

遙感圖像融合方法主要有4類:基于彩色空間的方法,如IHS(intensity hue saturation)方法;基于數(shù)理統(tǒng)計(jì)的方法,如PCA方法、非負(fù)矩陣分解方法;基于計(jì)算智能的方法,如基于神經(jīng)網(wǎng)絡(luò)、模糊理論的方法;基于圖像多尺度分析的方法,如金字塔變換方法、小波變換法、Bandelet變換、Curvelet變換和Contourlet變換方法[1].單一的IHS方法容易產(chǎn)生光譜失真,所以很多學(xué)者將IHS方法和其他類型的方法結(jié)合起來,加以改進(jìn)[2].在基于數(shù)理統(tǒng)計(jì)的方法中,PCA和非負(fù)矩陣方法[3]要求圖像具有強(qiáng)相關(guān)性,較適合應(yīng)用于多光譜和高光譜影像.當(dāng)前,小波變化、金字塔變換、Bandelet、Contourlet、Curvelet等是融合領(lǐng)域研究的熱點(diǎn),這些多尺度分析方法對圖像的邊緣、紋理等方面的信息捕捉非常好[4].由Peyré和Mallat提出的二代Bandelet變換是一種基于邊緣的圖像表示方法,能自適應(yīng)地跟蹤圖像的幾何正則方向,具有最大限度保留圖像細(xì)節(jié)的優(yōu)良特性[5].文獻(xiàn)[4-7]利用Bandelet變換研究光學(xué)圖像的融合,幾何流和Bandelet系數(shù)采用最大值、絕對值最大規(guī)則,或平均與選擇相結(jié)合.這樣的融合規(guī)則不適合高光譜圖像的融合.高光譜具有多個(gè)波段,波段之間相關(guān)性很強(qiáng),數(shù)據(jù)冗余度較大,光譜特征差異很小,存在噪聲干擾等特點(diǎn)[8].傳統(tǒng)的光譜融合方法處理效果不理想,不能同時(shí)克服去相關(guān)性和去噪等特點(diǎn)[9].因此,探討利用二代Bandelet變換方法獲得圖像的局部最佳幾何流和Bandelet系數(shù),再利用PCA方法融合最佳幾何流和 Bandelet系數(shù),經(jīng)過Bandelet逆變換重構(gòu)圖像,達(dá)到多尺度幾何分析融合圖像,同時(shí)減小了高光譜波段間冗余.

1 二代Bandelet變換算法

2000年,法國學(xué)者Pennec和M allat提出一種多尺度幾何分析方法——Bandelet變換.為了改進(jìn)第一代 Bandelet變換計(jì)算復(fù)雜的不足,2005年, Peyré和Mallat提出了第二代的Bandelet變換的思想,即把圖像中的幾何特征定義為矢量場,而不是普通邊緣的集.二代Bandelet變換能夠自適應(yīng)跟蹤圖像幾何結(jié)構(gòu),捕捉幾何正則性,給出圖像的最優(yōu)表達(dá).基于二代Bandelet變換的圖像融合是將變換得到的幾何流和Bandelet系數(shù)按照一定的規(guī)則進(jìn)行結(jié)合,再進(jìn)行逆變換得到融合圖像[4].二代Bandelet變化的步驟如下.

(1)正交小波對圖像進(jìn)行二維分解.

(2)對各子帶分別用自底向上的全局優(yōu)化算法建立最佳四叉樹分割,同時(shí)計(jì)算各分割區(qū)域內(nèi)的幾何流方向,記為:Gj(i)(j=1,2,…,N,N是圖像個(gè)數(shù);i=1,2,…,n,n代表Bandelet塊數(shù)).

(3)根據(jù)幾何流方向,再對各Bandelet塊實(shí)施Bandelet化,得到Bandelet系數(shù),記為Cj(i)(j=1, 2,…,N,N是圖像個(gè)數(shù);i=1,2,…,n;n代表Bandelet塊數(shù)).

1.1 四叉樹分割

首先對原圖像進(jìn)行二維正交小波變換.假定原圖像為f(x,y),用二尺度差分方程φ和ψ,沿著x和y方向先后進(jìn)行分解,得到平滑逼近和細(xì)節(jié)的部分,然后分解對角高頻信息.二維正交小波分解的公式如下:

采用文獻(xiàn)[10]中使用的二進(jìn)四叉樹分割方法,先將圖像等分成4個(gè)子帶,每一子帶在下一層的分割中又被分成4個(gè)子帶.依次分割下去,直到最底層的子帶的分解尺度達(dá)到預(yù)先設(shè)定的最小尺度Jmin.用四叉樹表示分解過程,四叉樹的葉節(jié)點(diǎn)對應(yīng)的就是Bandelet塊[11],如圖1所示.

圖1 二進(jìn)四叉樹分割示意圖Fig.1 Binary quadtree segmentation

1.2 最佳幾何方向的計(jì)算

各子帶的幾何流是從最小尺度子帶自下而上進(jìn)行計(jì)算的,計(jì)算步驟如下[12].

(1)角度的預(yù)采樣.為計(jì)算幾何流,先在子帶內(nèi)進(jìn)行角度的采樣,假設(shè)子帶寸大小為L×L,那么將圓周角[0,π)等角度離散為L2-1個(gè),采用角度可能的取值為

對于無幾何流的情形標(biāo)記為Inf,那么在L×L的子帶內(nèi)獲取了L2個(gè)采樣角度即(θ,I nf).

(2)曲波變換.曲波變換可自適應(yīng)地跟蹤具有曲線奇異特征的目標(biāo),能夠稀疏地表達(dá)目標(biāo).首先構(gòu)造一個(gè)與子帶同樣大小網(wǎng)格點(diǎn)(設(shè)子帶大小為L× L),然后計(jì)算每個(gè)網(wǎng)格點(diǎn)在采樣角度上的正交投影誤差為

式中:θ是采樣角度;x(i),y(j)為網(wǎng)格點(diǎn)坐標(biāo).網(wǎng)格點(diǎn)按誤差值從小到大排序成為一維數(shù)組,然后將與網(wǎng)格點(diǎn)位置相對應(yīng)的子帶系數(shù)進(jìn)行一維小波變換后即實(shí)現(xiàn)了曲波變換.

(3)選擇最優(yōu)幾何流方向.計(jì)算幾何流是通過最小化Lagrange函數(shù)獲取的,Lagrange函數(shù)的表達(dá)式如下:

計(jì)算各采樣角度的 Lagrange系數(shù),獲得最小Lagrange系數(shù)的采樣角度即為最優(yōu)幾何正則方向,也就是要進(jìn)行融合變換的最佳幾何流Gj(i).

1.3 Bandelet系數(shù)計(jì)算

Bandelet系數(shù)計(jì)算的過程又被稱為 Bandelet化,在第二代Bandelet變換中,Bandelet化的對象是小波變換的各高頻子帶的Bandelet塊,也就是經(jīng)過二進(jìn)制分割和四叉樹優(yōu)化之后的小區(qū)域;因?yàn)槌叽缱銐蛐?在每個(gè)分割區(qū)域內(nèi),圖像的幾何正則可以用直線逼近;然后在每個(gè)分割區(qū)域內(nèi)分別進(jìn)行Bandelet化.Bandelet化的實(shí)施過程包含了兩步:一是沿幾何流方向?qū)π〔ㄏ禂?shù)進(jìn)行重排,得到一維小波系數(shù);二是對一維小波系數(shù)做一維小波變換,得到Bandelet系數(shù).

假設(shè)某一高頻子帶Ωi里,二維小波系數(shù)為

式中:f[n1,n2]為子帶離散圖像;Φj,k1為尺度因子; Ψj,k2為彎曲小波函數(shù).對二維小波系數(shù)進(jìn)行重新排列,記為V′i,對 V′i進(jìn)行一維信號小波變換(使用haar小波)得到Bandelet系數(shù)Cj(i).

如果區(qū)域內(nèi)沒有明顯的幾何流方向,就沒有必要進(jìn)行Bandelet化,表明塊內(nèi)圖像灰度值的變化是一致正則的,直接保留二維小波變換系數(shù)即可,但仍統(tǒng)一稱作Bandelet系數(shù)G j(i)[13].得到Bandelet系數(shù)和幾何流,就可以按照一定融合規(guī)則進(jìn)行圖像融合.

2 PCA變換規(guī)則

PCA即主成分分析法,是一種常用降維方法. PCA變換的基本原理是將原始存在較高相關(guān)性的N維特征的高維數(shù)據(jù),通過線性變換投影到M維正交子空間(M?N),用新子空間中少數(shù)幾個(gè)不相關(guān)的、方差大的、信息量大的特征來表示原始數(shù)據(jù)信息[14],減小數(shù)據(jù)冗余[15].

經(jīng)過二代 Bandelet變換后得到Bandelet系數(shù)Cj(i)和最佳幾何流Gj(i).按照PCA變換,選取它們的主成分進(jìn)行融合.這里以Bandelet系數(shù)為例,論述PCA變換融合的規(guī)則;幾何流的融合規(guī)則與此相同.

N幅圖像二代Bandelet系數(shù)按照波段順序排列矩陣為C j(i),用向量的形式表示為C=(C1,C2,…,CN),對C進(jìn)行線性變換,形成新的綜合的變量Y,即

(1)Yi和Yj不相關(guān)(i≠j,i,j=1,2,…,p).

(2)Y1為C1,C2,…,CN的一切滿足公式(6)的線性組合中方差最大的,Y2是與Y1不相關(guān)的C1,C2,…,CN所有線性組合中方差最大的,YN是與Y1,Y2,…,YN-1不相關(guān)的C1,C2,…,CN所有線性組合中方差最大的.

Y1是第一主成分,其方差最大,包含原始數(shù)據(jù)的信息最多.公式(6)中的系數(shù)(u 1 j,u 2 j,…,u Nj)是C的特征向量,其相對應(yīng)的特征值為(λ1,λ2,…, λN).

在求解主成分時(shí),通常從分析原始變量C=(C1, C2,…,CN)的協(xié)方差矩陣和相關(guān)矩陣出發(fā),協(xié)方差陣和基于相關(guān)陣求出的主成分存在著較大的差異,但是在數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理后,結(jié)果是一致的[16].標(biāo)準(zhǔn)化之后,計(jì)算的速度也可提高.標(biāo)準(zhǔn)化公式如下:

式中:Cij是代表第i幅圖像的第j個(gè)Bandelet系數(shù),是所有圖像第j個(gè) Bandelet系數(shù)的平均值,是相應(yīng)的標(biāo)準(zhǔn)差.

3 基于二代Bandelet和PCA的高光譜遙感圖像融合方法

針對高光譜遙感圖像的特點(diǎn),采用二代Bandelet變換對圖像進(jìn)行多尺度幾何變換,利用PCA分析各個(gè)波段圖像的高頻和低頻信息,根據(jù)PCA變換分析結(jié)果,再重構(gòu)融合圖像.二代Bandelet多尺度幾何分析可彌補(bǔ)PCA在空域分析方面的不足,PCA可處理高光譜多波段數(shù)據(jù),減小冗余,兩者有機(jī)結(jié)合,能夠很好地融合高光譜圖像,取得良好的圖像處理效果.該方法的步驟如下.

(1)N個(gè)波段的高光譜圖像預(yù)處理、配準(zhǔn).

(2)分別對N個(gè)波段的圖像進(jìn)行Bandelet變換,其過程中二維àtrous[2]小波變換為三層分解;二進(jìn)制分割子帶的大小為4×4;Bandelet化時(shí)采用的是一維haar小波三層分解.得到每一個(gè)波段圖像的Bandelet系數(shù)Cj(i)和幾何流Gj(i).

(3)用 PCA變換分析 N個(gè)波段圖像的Bandelet系數(shù)和幾何流,獲得各自的主成分,其主成分分別記為C F(i),G F(i).

(4)利用Bandelet系數(shù)CF(i)和新的幾何流GF(i),Bandelet逆變換重構(gòu),得到融合新圖像.

4 實(shí)驗(yàn)分析

4.1 原始圖像分析

本文采用了8個(gè)波段、大小為256像素×256像素的高光譜數(shù)據(jù)(圖2),嚴(yán)格配準(zhǔn)后,再進(jìn)行實(shí)驗(yàn)分析.首先分析8個(gè)波段的原始圖像,直接觀察可以發(fā)現(xiàn):各幅圖像沒有明顯的噪聲,不同波段的圖像所包含的信息不同,前幾幅圖像波段對道路河流較為敏感;而后幾幅圖像信息量更大,容易判讀森林等植被信息.圖3是各個(gè)圖像的灰度直方圖,可看出多個(gè)波段的圖像之間存在著強(qiáng)的相關(guān)性,圖像之間的數(shù)據(jù)冗余量較大.

4.2 融合結(jié)果主觀評價(jià)

圖4a為本文方法融合的圖像;圖4b是利用Bandelet變換融合得到的圖像(8個(gè)圖像的Bandelet變換和本文第3節(jié)所闡述的過程相同,融合規(guī)則不同,幾何流和Bandelet系數(shù)分別采用最大值和絕對值最大規(guī)則);圖4c是PCA變換融合得到的圖像.3幅融合圖像比較,圖4a的目視效果最好,邊緣清晰,富含的信息更多,地物地貌的邊緣也較為清晰,驗(yàn)證了Bandelet變換善于捕捉圖像邊緣的優(yōu)點(diǎn);另外圖中道路、河流水體、樹木植被等都能較為清晰地識別,說明該法融合的圖像利于判讀和特征提取.而圖4b中河流水體和道路能夠清晰識別,森林植被較清晰.圖4c中道路、河流水體邊緣清晰,其他地物并不清晰.

圖3 8個(gè)波段圖像的直方圖Fig.3 Histogram s of the eight or iginal bands images

圖4 融合結(jié)果Fig.4 The fused images

4.3 融合結(jié)果客觀評價(jià)

為了從客觀的角度評價(jià)3種融合方法,計(jì)算了融合圖像的灰度平均值、標(biāo)準(zhǔn)差、熵、平均梯度,計(jì)算結(jié)果見表1.表中灰度均值比較接近,其他的指標(biāo)都有很大的提高.這樣的結(jié)果表明了本文方法的視覺效果最好,圖像包含的信息最豐富,圖像局部紋理特征分明,這從主觀視覺上也可驗(yàn)證,也說明了本文方法在融合性能提高方面具有很大的優(yōu)勢.另外,結(jié)合本文方法和二代Bandelet方法的原理,比較圖4a和4b的實(shí)驗(yàn)結(jié)果,說明PCA規(guī)則好于傳統(tǒng)的規(guī)則;但是比較3幅圖像的紋理信息,會發(fā)現(xiàn)PCA方法不如其他兩種方法,說明PCA融合在圖像的空域上分析效果有不足之處,所以應(yīng)該借助二代Bandelet這一多尺度幾何分析工具進(jìn)行彌補(bǔ).綜合來說,利用本文方法融合的圖像具有較好的效果和質(zhì)量.

表1 不同方法的統(tǒng)計(jì)評價(jià)Tab.1 Statistics evaluation of different algorithms

5 結(jié)語

針對高光譜圖像的特點(diǎn),利用具有自適應(yīng)跟蹤目標(biāo)幾何特征的二代Bandelet變換,以及具有減小數(shù)據(jù)冗余特點(diǎn)的PCA變換,進(jìn)行高光譜圖像融合.二代Bandelet變換靈活捕捉圖像的邊緣信息,自適應(yīng)地跟蹤圖像的幾何正則方向,彌補(bǔ)了PCA在空域變換上的不足;而PCA變換減小了高光譜波段數(shù)據(jù)的冗余性.經(jīng)試驗(yàn)分析,該方法取得了比其他兩種方法質(zhì)量更高的融合圖像,是處理高光譜圖像融合的一種理想方法,并且對紋理邊緣處理效果顯著.客觀評價(jià)指標(biāo)有較大的提高,尤其是熵的提高,說明該方法豐富了融合圖像的光譜信息.綜合主客觀評價(jià)可知,本文方法融合結(jié)果圖像層次明顯,地物地貌較清晰,邊緣輪廓分明,有利于圖像分類和提取的實(shí)現(xiàn).

[1] 郭雷,李暉暉,鮑永生.圖像融合[M].北京:電子工業(yè)出版社,2008.

GUO Lei,LI Huihui,BAO Yongsheng.Image fusion[M]. Beijing:Publishing House Electronics Industry,2008.

[2] 后斌,喬偉峰,孫在宏.基于IHS變換與àtrous小波分解的遙感影像融合[J].南京師大學(xué)報(bào):自然科學(xué)版,2006,29 (1):116.

HOU Bin,QIAO W eifeng,SUN Zaihong.Rem ote-sensing im age fusion based on IHS transform andàtrous w avelet decomposition[J].Journal of Nanjing Normal University: NaturalScience,2006,29(1):116.

[3] 陳鷹,郭睿.非負(fù)矩陣分解在遙感圖像融合中的應(yīng)用[J].計(jì)算機(jī)工程與應(yīng)用,2007,43(20):68.

CHEN Ying,GUO Rui.Non-negative matrix factorizations applied to im age fusion in remote sensing[J].Computer Engineering and Applications,2007,43(20):68.

[4] 閆敬文,屈小波.超小波分析及應(yīng)用[M].北京:國防工業(yè)出版社,2008.

YAN Jingw en,QU X iaobo.Analysis and application beyond w avelet[M].Beijing:NationalDefen ce Industry Press,2008.

[5] PeyréG,Mallat S.Discretebandelets with geometri orthogonal filters[C] ∥ IEEE International Conference on Im age Processing.Vancouver:[s.n.]2005:65-68.

[6] QU X iaobo,YAN Jinw en,X IE Guofu,et al.A novel im age fusion algo rithm based on bandelet transform[J].Chinese Optics Letters,2007,5(10):569.

[7] 焦李成,侯彪,王爽,等.圖像多尺度幾何分析理論與應(yīng)用[M].西安:西安電子科學(xué)技術(shù)大學(xué)出版社,2008.

JIAO Licheng,HOU Biao,WANG Shuang,et al.Im age mu ltiscale geometric analysis:theory and application[M].Xi' an:X i'an Electronic Science and Technology University Press, 2008.

[8] 張良培,張立福.高光譜遙感[M].武漢:武漢大學(xué)出版社,2005.

ZHANG Liangpei,ZHANG Lifu.H yperspectral remote sensing[M].Wuhan:W uhan University Press,2005.

[9] Shaw G,M ankolakis D.Signal p rocessing for hyperspectral spectral im age exploitation[J].IEEE Signal Processing Magazine,2002,19(1):12.

[10] Donoho D L.Wedgelets:nearly-m inimax estimation of edges [J].A nnals of Stat,1999,27:859.

[11] 梁銳華,成禮智,朱志誠,等.基于第二代Bandelet變換的圖像編碼[J].現(xiàn)代電子技術(shù),2007,4:61.

LIANG Ruihua,CHENG Lizhi,ZHU Zhicheng.Image coding based on second generation Bandelet transform[J].Modern Electronics Technique,2007,4:61.

[12] Pennec E L,Mallat S.Sparse geometrie im age representation w ith Bandelets[J].IEEE T rans on Image Processing,2005,14 (4):423.

[13] 龍奕.Bandelet變換在圖像去噪與增強(qiáng)中的應(yīng)用[D].成都:西南交通大學(xué)信息科學(xué)與技術(shù)學(xué)院,2008.

LONG Y i.The application of Bandelet transform in im age denoise and enhancem ent processing[D].Chengdu:Southw est Jiaotong University.College of Information Science and Technolgy,2008.

[14] Vogt F,Tacke M.Fast principal component analysis of large data sets[J].Chemometrics and Intelligent Laboratory Sy stems,2001,59(1):1.

[15] 戴昌達(dá),姜小光,唐伶俐.遙感圖像應(yīng)用處理與分析[M].北京:清華大學(xué)出版社,2004.

DAI Changda,JIANG Xiaoguang,TANG Lingli.Remote sensing image app lication processing and analysis[M].Beijing: Tsinghua University Press,2004.

[16] 李靜萍,謝邦昌.多元統(tǒng)計(jì)分析方法與應(yīng)用[M].北京:中國人民大學(xué)出版社,2008.

LI Jingping,XIE Bangchang.Method and application of mu ltivariate statistical analysis[M].Beijing:China Renm in University Press,2008.