祝傳廣 范洪冬 鄧喀中 薛繼群

(1)中國(guó)礦業(yè)大學(xué)江蘇省資源環(huán)境信息工程重點(diǎn)實(shí)驗(yàn)室,徐州 221116 2)河南理工大學(xué)礦山空間信息技術(shù)國(guó)家測(cè)繪局重點(diǎn)實(shí)驗(yàn)室,焦作 454000)

Contourlet域 InSAR干涉圖濾波方法研究＊

祝傳廣1,2)范洪冬1,2)鄧喀中1)薛繼群1)

(1)中國(guó)礦業(yè)大學(xué)江蘇省資源環(huán)境信息工程重點(diǎn)實(shí)驗(yàn)室,徐州 221116 2)河南理工大學(xué)礦山空間信息技術(shù)國(guó)家測(cè)繪局重點(diǎn)實(shí)驗(yàn)室,焦作 454000)

研究了基于 Contourlet變換的閾值、循環(huán)平移干涉圖濾波新方法,并與 Goldstein、小波軟閾值濾波方法在殘差點(diǎn)數(shù)目、相干系數(shù)均值等方面進(jìn)行了對(duì)比分析。實(shí)驗(yàn)表明,新算法對(duì)干涉圖噪聲有很強(qiáng)的抑制能力,較好地保持了干涉條紋的細(xì)節(jié),減少了殘差點(diǎn),處理效果好于傳統(tǒng)的小波軟閾值濾波方法,在某些情況下可與 Goldstein濾波效果相比。

干涉圖;時(shí)頻域;Contourlet變換;閾值;循環(huán)平移

1 引言

合成孔徑雷達(dá)干涉測(cè)量技術(shù)(InSAR)是通過(guò)利用同一地區(qū)的多幅 SAR復(fù)圖像獲得的地形相位信息,反演地形高度或地表形變。實(shí)際應(yīng)用中,多幅SAR影像間一般間隔時(shí)間較長(zhǎng),地物的變化、大氣的變化都會(huì)使干涉圖存在不同程度的失相關(guān),在干涉圖中產(chǎn)生噪聲誤差,致使解纏、生成DEM、探測(cè)形變等后續(xù)工作難度增加,結(jié)果可靠性也隨之降低。因此,必須對(duì)干涉圖進(jìn)行濾波。目前的主要濾波方法可分為空間域?yàn)V波和時(shí)頻域?yàn)V波?？臻g域?yàn)V波有均值濾波[1]、Lee濾波[2]等濾波方法;時(shí)頻域?yàn)V波主要是小波域?yàn)V波[3,4]。由于時(shí)頻濾波方法結(jié)合了空間與頻率域?yàn)V波,具有自動(dòng)調(diào)焦功能,對(duì)時(shí)間窗和頻率窗可做調(diào)節(jié)：在低頻段采用高頻率分辨率和低時(shí)間分辨率;在高頻段采用低頻率分辨率和高時(shí)間分辨率[5],濾波效果更為明顯[6,7]。

InSAR干涉圖的小波域?yàn)V波方法雖然對(duì)噪聲有較好的抑制作用,但是由一維小波通過(guò)張量積形成的二維可分離小波變換只能表示一維點(diǎn)奇異信息,而不能有效地描述圖像中的二維或高維奇異信息,如線、輪廓等[8],而干涉條紋圖呈周期變化,輪廓較為分明,采用小波變換對(duì) InSAR干涉圖進(jìn)行降噪處理時(shí),紋理等細(xì)節(jié)信息并不能很好地得到保留。因此,本文將 Contourlet變換引入到 InSAR干涉圖的濾波過(guò)程中,該變換既保留了多尺度分析的一些優(yōu)勢(shì),又具有小波變換所不具備的方向性和各向異性[9]。為獲取更好的濾波效果,得到高精度的解纏結(jié)果,本文研究了基于 Contourlet變換的閾值、循環(huán)平移干涉圖濾波新方法,并與 Goldstein、小波軟閾值濾波進(jìn)行了對(duì)比分析。

2 Contourlet域干涉圖濾波方法

2.1 Contourlet變換

Contourlet變換,同小波分解一樣,也是一種多分辨率、局域的、方向的影像表示方法。在分解圖像時(shí),Contourlet基的支撐區(qū)間具有隨尺度而長(zhǎng)寬比變化的長(zhǎng)條形結(jié)構(gòu),因此,對(duì)于細(xì)小的有方向的輪廓和線段的表達(dá)具有獨(dú)特的優(yōu)勢(shì)[8]。Contourlet變換將圖像分解成若干個(gè)尺度,在每個(gè)尺度上的方向數(shù)目(2的整數(shù)次冪)可以根據(jù)需要自適應(yīng)的調(diào)節(jié),因此, Contourlet變換不僅具有小波的優(yōu)良時(shí) (空)頻分析特征,而且表現(xiàn)出了高度的方向性[10],更適合于進(jìn)行多尺度的邊緣增強(qiáng)處理。Contourlet變換首先是對(duì)輸入的影像利用拉普拉斯金字塔分解得到一系列的低通子圖像和帶通子圖像,并“捕獲”不同尺度上的奇異點(diǎn);然后,采用方向?yàn)V波器組對(duì)得到的帶通子圖像進(jìn)行方向性濾波,將頻域分解成 2j個(gè)楔型子帶,并將分布在同一方向的奇異點(diǎn)合成為一個(gè)系數(shù)。通過(guò) Contourlet變換對(duì) SAR圖像進(jìn)行分解時(shí),奇異性特征可以在更多的方向得到自適應(yīng)的多分辨率刻畫,細(xì)節(jié)在方向子帶得到更為顯著的呈現(xiàn)[10]。

2.2 閾值去噪

小波閾值去噪的理論依據(jù)為：對(duì)信號(hào)進(jìn)行變換后,包含有重要信息的系數(shù)其幅值較大,但數(shù)量少;噪聲對(duì)應(yīng)的系數(shù)數(shù)量多,但幅值小[11]。因此,選擇合理的閾值可以獲得較好的去噪效果。

同小波變換類似,Contourlet變換具有去相關(guān)的特性,使得影像變換后的能量集中在有限的變換域系數(shù)上,其余大部分變換域系數(shù)的幅值接近于零,而噪聲經(jīng)變換的能量均勻地分布在所有的變換域系數(shù)上。最常用的閾值去噪函數(shù)有硬閾值函數(shù)和軟閾值函數(shù)公式為[12]：

式中,x是 Contourlet變換之后的系數(shù);T是設(shè)定的閾值,該值由Donoho[13]所給的方法獲取。

2.3 循環(huán)平移濾波方法

Contourlet變換由于缺乏平移不變性,閾值去噪后重構(gòu)的影像在奇異點(diǎn)附近會(huì)出現(xiàn)偽 Gibbs現(xiàn)象,因此,采用閾值去噪會(huì)在重構(gòu)的影像中產(chǎn)生虛假成分。為解決這一問(wèn)題,可以先對(duì)原始影像平移一定的距離,然后對(duì)平移后的影像進(jìn)行Contourlet變換得到變換域系數(shù),并進(jìn)行軟閾值去噪,最后對(duì)去噪后的系數(shù)進(jìn)行 Contourlet逆變換并反向平移得到濾波后的影像。對(duì)于一個(gè)給定的平移范圍[N1,N2],循環(huán)平移去噪算法可表示為[14]:

3 干涉圖濾波效果評(píng)價(jià)

3.1 殘差點(diǎn)

干涉圖中殘差點(diǎn)的數(shù)目反映去噪效果的好壞,殘差點(diǎn)越少,去噪效果越好。殘差點(diǎn)的計(jì)算方法如圖 1[15]：

圖1 殘差點(diǎn)計(jì)算路徑Fig.1 Route of calculation of residuals

沿箭頭所指方向,相鄰兩像素間求相位差值Δ,然后對(duì)Δ求和,即:

式中,若 S=0,則不存在殘差點(diǎn),否則,存在殘差點(diǎn)。對(duì)整幅干涉圖計(jì)算,得到 S≠0的點(diǎn)的總數(shù),即為殘差點(diǎn)數(shù)。

3.2 偽相干系數(shù)均值

偽相干值的高低表明了圖像的相干性的好壞,值越大表明相干性越好。計(jì)算公式為[16]:

式中,k為窗口大小。計(jì)算出每一個(gè)像元的偽相干值,即可求得一副干涉圖像的偽相干系數(shù)均值。

3.3 相位導(dǎo)數(shù)均值

相位導(dǎo)數(shù)變化計(jì)算公式為[16]

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

實(shí)驗(yàn)所使用的干涉數(shù)據(jù)為 ERS1/2影像。第一幅干涉圖位于江蘇平原地區(qū) (因平原地區(qū)高差小,干涉條紋不夠明顯,實(shí)驗(yàn)是針對(duì)未去平地效應(yīng)的干涉圖操作的),第二幅干涉圖位于青藏高原地區(qū)。4景影像均為 C波段,VV極化方式,影像參數(shù)見(jiàn)表 1。

表1 影像參數(shù)Tab.1 Parameters of the i mage

Goldstein濾波效果與選取的冪有關(guān),范圍一般在 0.4～0.7,冪越大,濾波程度大,但也有可能過(guò)度濾波。本文采用Matlab語(yǔ)言實(shí)現(xiàn)了所研究的 Contourlet域軟閾值及循環(huán)平移干涉圖濾波新算法,并與Goldstein濾波方法及小波軟閾值濾波方法進(jìn)行了對(duì)比,其中,Goldstein濾波方法的冪取 0.5與 0.7。原始干涉圖、Goldstein(0.5)、Goldstein(0.7)、DWT軟閾值、Contourlet軟閾值及 Contourlet循環(huán)平移濾波方法獲取的干涉圖分別見(jiàn)圖 2、圖 3中的 (a)～(f)。圖 3影像大小為 290×290個(gè)像元,圖 4影像大小為256×256個(gè)像元。對(duì)干涉圖進(jìn)行了 3層 LP分解,每層的細(xì)節(jié)方向數(shù)分別為 4、8和 16。

由干涉圖 2、3可以看出：在平原與山區(qū),本文研究的 Contourlet循環(huán)平移濾波方法所得到的結(jié)果在清晰度上有明顯提高,干涉條紋邊緣層次感強(qiáng),比較連續(xù),均要優(yōu)于 Contourlet、小波軟閾值濾波算法得到的結(jié)果;而同 Goldstein濾波方法相比,在不同地形區(qū)域時(shí),Contourlet軟閾值、小波軟閾值及 Goldstein濾波算法則表現(xiàn)出不同的效果。

1)由圖 2可以看出,在平原地區(qū),Contourlet軟閾值與小波軟閾值濾波方法的結(jié)果各有優(yōu)劣：在干涉較差的地區(qū)(圖框內(nèi)部),小波軟閾值濾波效果要優(yōu)于 Contourlet軟閾值濾波;而在干涉較好的地區(qū)(圖框外圍),Contourlet軟閾值濾波則要優(yōu)于小波軟閾值濾波,紋理更為清楚。Contourlet循環(huán)平移濾波方法略優(yōu)于 Goldstein(0.5),而與 Goldstein(0.7)相比則較差。

2)由圖 3可以看出,在山區(qū),在地形陡峭 (即條紋密集處,圖框內(nèi)部所示)及平緩區(qū)域 (圖框外部所示),Contourlet軟閾值濾波方法均要優(yōu)于小波軟閾值濾波方法;并且 Contourlet循環(huán)平移濾波方法明顯優(yōu)于 Goldstein(0.5),與 Goldstein(0.7)相近,在地形陡峭的區(qū)域,甚至優(yōu)于 Goldstein(0.7)。

濾波后干涉圖,從目視角度看,Contourlet循環(huán)平移濾波方法在平原地區(qū)不如選取適當(dāng)冪值的 Goldstein濾波方法;但是在山區(qū),Contourlet循環(huán)平移濾波方法則要優(yōu)于 Goldstein濾波方法。

為精確評(píng)估幾種濾波方法的效果,表 2列出了3種評(píng)價(jià)濾波效果的指標(biāo)值。由表 2可以看出, Contourlet濾波及其循環(huán)平移濾波方法和DWT方法都能夠較好地提高濾波后干涉圖的相干性,相干系數(shù)得到提高,殘差點(diǎn)大量減少,相位導(dǎo)數(shù)變化程度得到有效抑制。本文研究的 Contourlet循環(huán)平移濾波效果要明顯優(yōu)于 Contourlet軟閾值濾波與小波軟閾值濾波。并且,目視效果與數(shù)據(jù)顯示一致,即：

1)在平原地區(qū),Contourlet軟閾值與小波軟閾值濾波方法效果等價(jià),但均不如Contourlet循環(huán)平移濾波方法;Contourlet循環(huán)平移同 Goldstein(0.5)濾波效果等價(jià),但不如 Goldstein(0.7)。

2)在山區(qū),小波軟閾值濾波、Contourlet軟閾值、Contourlet循環(huán)平移濾波效果依次變好;并且,Contourlet循環(huán)平移濾波優(yōu)于 Goldstein(0.5),同 Goldstein(0.7)相近。

圖 2 平原地區(qū)各種方法的濾波結(jié)果Fig.2 Filtering resultswith differentmethods in the plain area

圖3 山區(qū)各種方法的濾波結(jié)果Fig.3 Filtering resultswith differentmethods in the mountainous area

Contourlet循環(huán)平移濾波之所以在第二幅干涉圖中濾波效果較 Goldstein好,分析原因有：

1)由表 1可見(jiàn),生成第二幅干涉圖的兩景影像時(shí)間基線僅為 1天,空間基線為 103 m,而第一幅干涉圖的兩景影像時(shí)間基線為 71天,空間基線為 55 m。因此,第一幅干涉圖受到的時(shí)間失相干較第二幅干涉圖嚴(yán)重,而兩幅干涉圖受到的空間失相干則相差不多。

2)Contourlet的優(yōu)勢(shì)在于可以較好描述線、輪廓。在平原地區(qū),輪廓不明顯;在山區(qū)則不同,地形起伏較大,輪廓非常明顯,Contourlet的優(yōu)勢(shì)可以得到較好的發(fā)揮。

表 2 濾波前后評(píng)價(jià)指標(biāo)的比較Tab.2 Comparison between esti mate i ndexes before and after filteri ng

5 結(jié)語(yǔ)

將 Contourlet變換引入到 InSAR干涉圖的濾波過(guò)程中,研究了基于 Contourlet變換的軟閾值、循環(huán)平移濾波新方法,并與 Goldstein、小波軟閾值濾波方法進(jìn)行了對(duì)比。實(shí)驗(yàn)表明,Contourlet循環(huán)平移濾波在有效去除噪聲的同時(shí),較大程度地保留了細(xì)節(jié)和邊緣信息,在高、低相干區(qū)域的處理都要明顯好于小波軟閾值方法,在某些情況下可與 Goldstein相當(dāng)。為了得到更好的濾波效果,在后續(xù)研究中還要對(duì)閾值的選取方法、邊緣保持、尺度相關(guān)去噪等方面做進(jìn)一步的研究。

1 Eichel P H and Guiglia D C.Spotlight SAR interferometry for terrain elevation mapping and interometric change detection[R].Sand：Sandia Nationl Labs Technology,1993, 2 593-2 546.

2 Lee J S.Digital image enhancement and noise filtering by use of local statistics[J]. IEEE Trans. PAM I,1980, 2(2)：165-168.

3 W itkin A.Scale space filtering：A new approach to multiscale description[A].Acoustics,speech,and signal processing,IEEE InternationalConference on I CASSP’84[C]. IEEE Conferences,1984,150-153.

4 徐晨,等.小波分析應(yīng)用算法[M].北京：科學(xué)出版社, 2004.(Xu Chen,et al.Wavelet analysis and its application algorithm[M].Beijing：Science Press,2004)

5 Mallat S.A theory for multiresolution signal decomposition wavelet representation[J]. IEEE Trans.Pattern Anal. Mach.Intell.,1989,11(7)：674-691.

6 靳國(guó)旺,等.InSAR干涉圖的矢量分離式小波濾波[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),2008,33(2)：132-135.(Jin Guowang,et al.Filtering for InSAR interferograms by vector decomposing and wavelet transfor mation[J].Geomatics and Information Science of Wuhan University,2008,33(2)：132-135)

7 蔡國(guó)林,等.小波—維納組合濾波算法及其在 InSAR干涉圖去噪中的應(yīng)用 [J].遙感學(xué)報(bào),2009,13(1)：129-136.(Cai Guolin,et al.Wavelet-W iener combined filter and its application on InSAR interferogram[J].Journal of Remote Sensing,2009,13(1)：129-136)

8 林立宇,等.Contourlet變換-影像處理應(yīng)用[M].北京：科學(xué)出版社,2008.(Lin Liyu,et al.Contourlet transfor mimage precessing and applying[M].Beijing：Science Press, 2008)

9 Do M N and VetterliM.The contourlet transform：An efficient directional multiresolution image representation[J]. IEEE Transactions on I mage Processing,2003,14(12)：2 091-2 106.

10 焦李成,等.智能 SAR圖像處理與解譯[M].北京：科學(xué)出版社,2008.(Jiao Licheng,et al.Intelligent SAR image processing and interpretation[M].Beijing：Science Press, 2008)

11 Donoho D L and Johnstone IM.Ideal spatial adaptation by wavelet shrinkage[J].Biometrika,1994,81(3)：425-455.

12 Chang S C,et al.Adaptive wavelet thresholding for image denoising and compression[J].IEEE Transaction on I mage Processing,2000,9(9)：1 532-1 546.

13 Donoho D L.De-noising by soft-thresholding[J].IEEE Trans.Information Theory,1995,41(5)：613-627.

14 王貞儉.一種 Contourlet循環(huán)平移自適應(yīng)閾值 SAR圖像去噪算法[J].電光與控制,2009,16(4)：51-54.(Wang Zhenjian.A contourlet cycle-spinning and adaptive threshold SAR speckle reduction algorithm[J].ElectronicsOptics amp;Control,2009,16(4)：51-54)

15 靳國(guó)旺.InSAR獲取高精度 DEM關(guān)鍵處理技術(shù)研究[D].解放軍信息工程大學(xué),2007.(Jin Guowang.Research on key processing techniques for deriving accurate DEM from InSAR[D].Infor mation EngineeringUniversity, 2007)

16 王超,等.星載合成孔徑雷達(dá)干涉測(cè)量[M].北京：科學(xué)出版社,2002.(Wang Chao,et al.Spaceborne synthetic aperture radar interferometry[M].Beijing：Science Press, 2002)

STUDY ON FILTERING M ETHODS OF InSAR INTERFEROGRAM IN CONTOURLET DOMA IN

Zhu Chuanguang1,2),Fan Hongdong1,2),Deng Kazhong1)and Xue Jiqun1)

(1)China University of M ining and Technology,Jiangsu Key Laboratory of Resources and Environm ental Infor m ation Engineering,Xuzhou221116 2)Henan Polytechnic University,Key Laboratory of M ine Spatial Infor m ation Technologies,Jiaozuo 454000）

Some new ti me-frequency filtering methods such as threshold,cycle spinning filtering algorithms based on Contourlet transform were studied.After compared and analyzed with the filteringmethods based on DWT (DiscreteWavelet Transform)and the Golstein in the number of residues and the coherence value it can be seen that the new algorithms keep the details of interferogram verywell,reduce more residues,and the filtering effect is much better than the classical soft threshold method based onDWT,and even as good as Goldstein in a certain condition.

interferogram;ti me-frequency domain;Contourlet Transfor m;threshold;cycle spinning

1671-5942(2010)05-0116-05

2010-03-31

礦山空間信息技術(shù)國(guó)家測(cè)繪局重點(diǎn)實(shí)驗(yàn)室(河南理工大學(xué),河南省測(cè)繪局)開(kāi)放基金(KLM200909);國(guó)家教育部博士點(diǎn)基金(20090095110002);國(guó)家自然科學(xué)基金(41071273)

祝傳廣,男,1984年生,碩士研究生,主要研究方向?yàn)檫b感數(shù)據(jù)處理及變形監(jiān)測(cè).E-mail：zhuchuanguang@163.com

P225.1

A