一種改進的InSAR干涉圖復數(shù)空間自適應濾波

易輝偉 ，朱建軍 ，陳建群 ，齊艷妮，李佳

(1. 中南大學 地球科學與信息物理學院，湖南 長沙，410083；2. 湖南省普通高校精密工程測量及形變?yōu)暮ΡO(jiān)測重點實驗室，湖南 長沙，410083)

合成孔徑干涉雷達(InSAR)能夠全天候、全天時、無接觸式地獲取地球表面形變的精密信息，已廣泛用于地震、滑坡、冰移、火山及礦山沉降等現(xiàn)象監(jiān)測中。然而，受熱噪聲、時間去相關(guān)、空間去相關(guān)等因素影響[1]，INSAR干涉圖像不可避免地存在大量噪聲，難以解纏以至于無法得到數(shù)字地面高程模型及形變結(jié)果。因此，如何更好地對干涉圖進行去噪，成為INSAR實際應用中亟待解決的問題[2]。干涉圖濾波方法分為頻率域和空間域2類，頻率域濾波是將圖像的滑動小窗口(通常為 32像素×32像素)內(nèi)像素通過傅里葉變換為功率譜，并對該功率譜進行平滑再通過反傅里葉變換到圖像空間。頻率域濾波方法的代表為Goldstein[3]濾波。Baran等[4]應用相干值改善Goldstein的濾波因子，使其更具有自適應性。Li等[5]引入干涉圖的視數(shù)及相位標準偏差來調(diào)整濾波因子，使干涉圖相位信息保持和算法穩(wěn)健性得到有效改善；孫倩等[6]應用信噪比控制濾波強度，對低信噪比區(qū)域強濾波，高信噪比區(qū)域弱濾波。于曉歆等[7]應用最大最小值法對相位標準偏差進行歸一化處理后采用其指數(shù)作為濾波因子，得到較好濾波效果。Goldstein方法去噪能力較強，但會損失大量邊緣細節(jié)信息，給 DEM和形變測量帶來誤差。近年發(fā)展的矢量分離式小波濾波[8]和小波-維納組合濾波[9]能在一定程度上保留細節(jié)，其小波基的構(gòu)造及小波軟硬閾值的選取還需要良好的智能決策。卡爾曼濾波[10]和基于格網(wǎng)值法[11]都是通過估計相位來濾波圖像，不同的估計方法將得到不同的濾波結(jié)果?？臻g域是逐像素進行的，通過周圍像素來影響中心像素值達到去噪目的，將模板(通常為3×3像素)與中心像素周圍對應位置進行卷積，得到新的中心像素。均值濾波和中值濾波簡單易行但未考慮干涉圖噪聲強弱的變化及圖像的條紋特征，易引起條紋混迭。作為空間域濾波方法代表的Lee濾波[12]通過局部解纏估計中心像素鄰域內(nèi)的坡度，從而根據(jù)乘性噪聲模型改善中心像素。由于帶噪聲的局部解纏計算復雜、費時且不準確，Lee提出了修正Lee濾波，即將包含邊緣的局部區(qū)域分為8個方向，利用局域梯度邊緣檢測方法確定邊緣的方向，并判定中心像素所在方向，并用該方向的像素取值改變中心像素，此方法能更好地保護條紋邊緣；尹宏杰等[13]根據(jù)同質(zhì)區(qū)域方差最小的原理確定中心像素所在方向，并根據(jù)干涉圖的相干性來選擇融合方向窗口的數(shù)量，增強Lee算法的穩(wěn)健性，但算法極為復雜。數(shù)學形態(tài)法[14]在Lee濾波的非同質(zhì)區(qū)采用膨脹腐蝕法找到邊緣，能有效增強邊緣。王興旺等[15]提出邊緣保持法，其原理也是找到中心像素所在方向，該方法受窗口大小和噪聲強弱影響較大。文獻[16-17]應用區(qū)域增長法并以相位變化最小為準則自適應選擇濾波窗口，取得較好濾波效果，但未能考慮條紋率。李佳等[18]結(jié)合頻率域和空間域濾波的優(yōu)點設計二級濾波，并取得一定成效。近年來發(fā)展的堆棧法[19]將圖像分成若干個二進制文件再通過布爾判斷式分別濾波。乘性加性模型法[20]將對角線相位噪聲劃為乘性，非對角線噪聲劃為加性，該方法有待進一步驗證。廖明生等[21]提出的復數(shù)空間自適應濾波是根據(jù)梯度確定中心像素縱橫方向鄰域內(nèi)各像素的權(quán)系數(shù)，卷積后得到中心像素的新值，其算法簡單，效果明顯，迭代3～5次即可。本文作者首先介紹復數(shù)空間自適應濾波算法，通過增加 45°和 135°梯度計算進行斜方向條紋的探索，并根據(jù)相干系數(shù)確定條紋邊緣的尺度函數(shù) k，進而改善權(quán)系數(shù)的表達式；然后，以模擬數(shù)據(jù)和真實數(shù)據(jù)作為實驗數(shù)據(jù)，分別與復數(shù)空間自適應、Lee濾波、基于信噪比的Goldstein及最優(yōu)化方向融合方法進行比較，結(jié)果表明本文方法在去噪和邊緣信息保持方面較前4種方法都有明顯改善，局部自適應性更強。

1 復數(shù)空間自適應濾波方法

InSAR干涉圖由復數(shù)構(gòu)成，圖像顯示的數(shù)據(jù)為[-π，+π]之間的纏繞相位存在跳躍性，因此不能對其直接平滑，而應將復數(shù)的實部和虛部分別濾波，再進行復數(shù)重構(gòu)，復數(shù)空間自適應濾波方法的算法可用下式表示。

其中： f(t)(x + i , y + i )為 (x + i, y + j )處像素的第t次迭代， ω(t)(x+ i , y +j)為(x + i,y + j )處像素的權(quán)系數(shù)。其思想是：(x，y)處像素值的濾波值是以其周圍 3×3像素為窗口，對窗口內(nèi)每個像素進行加權(quán)平均而得。權(quán)系數(shù)可定義為：

其中：s′(t)(x)為信號 s(t)(x)的一階導數(shù)，在此為梯度計算，見式(3)和(4)；k為尺度函數(shù)，確定在平滑過程中可保留的邊緣幅度，通常取梯度最大值的1/3～1/2。如果某點的信號梯度大于k，則認為該點處于邊緣上，邊緣點在濾波窗口對應的權(quán)系數(shù)減小，即盡量不參與平滑，使邊緣得以保留；反之，若信號梯度小于 k，則對應非邊緣的點，其權(quán)系數(shù)就會變大，在非邊緣區(qū)域梯度低于 k，窗口內(nèi)各權(quán)系數(shù)接近，在較大程度上參與平滑，從而達到去除噪聲的效果。

其中：Gx(x，y)為(x，y)像素的水平方向梯度；Gy(x，y)為(x，y)像素的垂直方向梯度。

2 改進的復數(shù)空間自適應濾波

復數(shù)空間自適應算法存在不足，如梯度計算只考慮水平和垂直方向，不夠全面；尺度函數(shù)k的取值為1/3～1/2的最大梯度這一經(jīng)驗值。事實上，噪聲強弱會直接影響到梯度的變化，條紋尺度函數(shù)的固定取值不能適應噪聲的強弱變化，會降低對條紋判斷的準確性，其局部適應性不強；此外，人工選取k費時且主觀性強。針對以上不足，作如下改進：

(1) 梯度計算中加入斜方向梯度。3×3像素的濾波窗口如圖1所示。從圖1可見：像素5周圍有8個像素，梯度方向不僅包括4，5和6水平方向以及2，5和8垂直方向，還應包括3，5和7所形成的45°以及1，5和9 所形成的135°這2個斜方向，其梯度計算見式(5)和(6)。因而，濾波窗口權(quán)系數(shù)相應地改為式(7)；

圖1 3像素×3像素的濾波窗口Fig.1 3×3 pixels filter window

其中：G45(x，y)為 45°方向梯度；G135(x，y)為 135°方向梯度。

(2) 用相干系數(shù)改善尺度函數(shù)k。如前所述，k取值大則平滑程度大，即強去噪；取值小則平滑程度小，即弱去噪。而INSAR干涉圖中噪聲分布并不均勻，有的區(qū)域噪聲強，有的區(qū)域噪聲弱，對其濾波較理想的結(jié)果是強噪?yún)^(qū)得到強濾波，弱噪?yún)^(qū)得到弱濾波并能保持條紋形狀，顯然k取固定值不能滿足這一要求。噪聲的強度由噪聲相位標準偏差決定，噪聲相位標準偏差大，則噪聲強；噪聲相位標準偏差小，則噪聲弱，噪聲相位標準偏差與相干系數(shù)和視數(shù)存在以下關(guān)系[22]：

其中： σφ2為相位方差；γ為干涉圖相干系數(shù)；L為視數(shù)；φ為干涉圖相位；φ0為相位的期望；fPDF為相位概率密度函數(shù)。

圖2 噪聲相位標準偏差與相干系數(shù)及視數(shù)關(guān)系Fig.2 Relationship between phase standard deviation and coherence&look number

α與視數(shù)L有關(guān)，通過大量數(shù)據(jù)模擬試驗得到：

3 實驗結(jié)果及分析

3.1 模擬干涉圖實驗結(jié)果

利用多分形技術(shù)模擬數(shù)字高程模型，再根據(jù)ERS1/2的成像幾何參數(shù)和200 m垂直基線模擬出纏繞的相位，最后根據(jù)式(9)模擬出給定視數(shù)的相位噪聲。圖3(a)所示為512像素×512像素的模擬無噪相位圖，圖3(b)所示為加噪相位圖，圖3(c)～(g)所示為分別用復數(shù)空間自適應濾波(取尺度函數(shù)為1/3kmax)、Lee濾波、基于信噪比的Goldstein、最優(yōu)方向融合濾波以及本文方法對其進行濾波的結(jié)果圖。從圖3可見：5種方法濾波效果都不錯，基本復原無噪相位圖。本文方法結(jié)果(圖3(g))與圖3(c)相比有明顯改善，噪聲大量減少，條紋更清晰。整體效果與圖 3(d)及圖 3(f)中的結(jié)果相似，但條紋邊緣細微處信息保持得更好，與圖 3(a)更為接近，且運行時間要比Lee濾波及最優(yōu)方向融合濾波快很多。比較圖 3(e)與圖 3(a)可見：基于信噪比的Goldstein濾波去噪能力較強，特別是同質(zhì)區(qū)內(nèi)的顯著斑點噪聲，但同時將一些細節(jié)真實信息也當成噪聲被濾掉；此外，由于該算法本身的原因還需處理圖邊像素。無論從噪聲的消除程度還是邊緣信息的保護程度都可看出本文方法具有明顯的優(yōu)勢。圖4所示為與圖3中對應的第280行剖面圖。從該剖面圖也可見：本文方法的剖面線光滑程度更高，噪聲的抑制效果非常明顯，更好地恢復原始無噪相位，與無噪相位圖的剖面圖4(a)更為吻合。

3.2 真實干涉圖實驗結(jié)果

本文分別選取香港地區(qū)和意大利 Etna火山局部區(qū)域的干涉圖作為實驗數(shù)據(jù)，香港地區(qū)數(shù)據(jù)為 1幅1996-03-18—19 ERS—1/2衛(wèi)星香港地區(qū)的InSAR干涉圖，其垂直基線長為100 m。該圖條紋密度適中，在此取700像素×1 000像素矩形區(qū)域作為實驗區(qū)域，其濾波結(jié)果如圖5所示。由于真實圖噪聲的復雜性，上述幾種濾波方法的效果與模擬圖相比稍弱，但都明顯地去除噪聲，增強條紋。為便于觀察，選定橫向800～925像素、縱向340～425像素區(qū)域局部放大，如圖6所示。圖6(a)所示為分布大量噪聲的原始干涉圖，圖6(b)～(f)所示分別為仍是復數(shù)空間自適應濾波、Lee濾波、基于信噪比的Goldstein、最優(yōu)方向融合濾波以及本文方法的濾波結(jié)果。從圖6可見：本文方法和復數(shù)空間自適應濾波基本去除噪聲，而其他3種方法還留有殘余噪聲。而圖 6(f)與圖 6(b)相比，既保持條紋的曲度又適當平滑條紋邊緣，效果更好。運行速度方面，Lee濾波所需時間8 min，最優(yōu)方向融合需要4 min，其他方法只需2 min左右。從第240行的相位剖面圖7(a)～(f)可進一步看出：圖 7(c)～(e)中還有明顯較多的噪聲，黑框中的圖7(b)還殘存毛刺；圖7(c)～(e)中的相位受噪聲影響存在明顯的迂回停滯，本文方法的剖面圖 7(f)幾乎沒有毛刺，噪聲也大量減少，整條曲線比較平滑，不僅去噪效果良好，而且有效保持相位細節(jié)。

意大利Etna火山100×200像素區(qū)域的密集條紋干涉圖及濾波結(jié)果如圖8所示?？臻g自適應濾波的條紋稍顯粗糙，而Lee濾波部分條紋明顯斷開，條紋“僵直”；最優(yōu)方向融合法條紋很不明顯，且光滑度不好；基于信噪比的 Goldstein方法的條紋光滑度及曲度保持也不如本文方法明顯。

圖3 模擬干涉圖濾波結(jié)果Fig.3 Simulated interferograms and results of different filters

圖4 模擬干涉圖濾波結(jié)果第280行剖面圖Fig.4 Profiles of line 280 of filtering results of simulated interferogram

圖5 香港地區(qū)干涉圖濾波結(jié)果Fig.5 Filtering results of interferogram in Hong Kong region

圖6 香港濾波結(jié)果局部放大圖Fig.6 Enlarged views of interferogram in local Hong Kong region

圖7 香港干涉圖濾波結(jié)果第240行剖面圖Fig.7 Profiles of line 240 of filtering results in Hong Kong region

圖8 意大利Etna火山截取干涉圖濾波結(jié)果Fig.8 Filtering results of interferogram over part of Etna volcano in Italy

表1 干涉圖濾波結(jié)果定量比較Table 1 Quantative comparison among filtering results of simulated interferogram

仍以模擬圖和香港干涉圖為例，用殘余點數(shù)、相位梯度和(簡稱為SPD)對上述5種濾波方法進行定量評價，結(jié)果見表 1。模擬圖的各種濾波方法的指標均較理想，而復數(shù)空間自適應濾波由于模擬圖的噪聲較多而受影響，殘余點改善率稍低，說明kmax/3的尺度函數(shù)未能適應模擬圖的噪聲水平。從表1可知：雖然香港干涉圖噪聲不強，但由于真實噪聲分布不均勻，隨機性大，因而各項指標均有所下降。此處復數(shù)空間自適應濾波下降不大，說明該尺度函數(shù)更適應于香港干涉圖，從側(cè)面說明尺度函數(shù)應該隨噪聲水平而變。無論模擬圖還是真實圖，本文方法與復數(shù)空間自適應方法相比都有較大改善，且優(yōu)于其他3種方法，模擬圖殘余點改善率由 79.6%提高到 98.2%，香港干涉圖改善率由 74.6%提高到 91.8%，分別提高 18.6%和17.2%，說明本文方法更具有優(yōu)勢。

4 結(jié)論

(1) 復數(shù)空間自適應濾波的目的是找到與原始像素值相近的鄰域像素盡可能參與平滑，因此，梯度小則認為該點與原始像素相近，相應地將其權(quán)重加大，與原始像素一起平滑；而梯度大則認為該點與原始像素不在同一條邊緣上，將其權(quán)重減小，盡可能不參與平滑。實質(zhì)上是給原始像素找到同質(zhì)或相近點，并對這些點根據(jù)接近程度進行加權(quán)平均。

(2) 由于條紋走向有可能是斜方向，因此，考慮45°和 135°方向的梯度，可以更準確地判斷原始中心像素與周圍像素的關(guān)系，有助于找到遺漏的同質(zhì)點，從而更合理地確定權(quán)系數(shù)。

(3) 噪聲強弱決定濾波強弱，噪聲可由干涉圖相干系數(shù)表征，相干系數(shù)小時則噪聲強，此時較大的尺度函數(shù)才能檢測出條紋；相反，相干系數(shù)大時則噪聲強，此時較小的尺度函數(shù)就能檢測出條紋。因而加入相干系數(shù)的干預，使得復數(shù)空間自適應濾波的穩(wěn)健性更強，局部自適應能力更強。

[1] Zebker H A, Villasenor J. Decorrelation in interferometric radar echoes[J]. IEEE Transactions on Geoscience and Remote Sensing, 1992, 30(5): 950-959.

[2] LI Zhiwei, DING Xiaoli, LIU Guoxiang. Modeling atmospheric effects on InSAR with meteorological and continuous GPS observations: Algorithms and some test results[J]. Journal of Atmospheric and Solar-Terrestrial Physics, 2004, 66: 907-917.

[3] Goldstein R M, Werner C L. Radar interferogram filtering for geophysical application[J]. Geophysical Research Letters, 1998,25(21): 4035-4038.

[4] Baran I, Stewart M P, Kampes B M, et al. A modification to the Goldstein radar interferogram filter[J]. IEEE Transaction Geoscience and Remote Sensing, 2003, 41(9): 2114-2118.

[5] LI Zhiwei, DING Xiaoli. Improved filtering parameter determination for the goldstein radar interferogram filter[J].ISPRS Journal of Photogrammetry and Remote Sensing, 2008,63: 621-634.

[6] 孫倩, 朱建軍, 李志偉, 等. 基于信噪比的InSAR干涉圖自適應濾波[J]. 測繪學報, 2009, 38(5): 437-442.SUN Qian, ZHU Jianjun, LI Zhiwei, et al. A new adaptive InSAR interferogram filter based on SNR[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(5): 437-442.

[7] 于曉歆, 楊紅磊, 彭軍還. 一種改進的 Goldstein InSAR干涉圖濾波算法[J]. 武漢大學學報: 信息科學版, 2011, 36(9):1051-1054.YU Xiaoxin, YANG Honglei, PENG Junhuan. A modified Goldstein algorithm for InSAR interferogram filtering[J].Geomatics and Information Science of Wuhan University, 2011,36(9): 1051-1054.

[8] 靳國旺, 韓曉丁, 賈博, 等. InSAR干涉圖的矢量分離式小波濾波[J]. 武漢大學學報: 信息科學版, 2008, 33(2): 132-135.JIN Guowang, HAN Xiaoding, JIA Bo, et al. Filtering for InSAR interferograms by vector decomposing and wavelet transformation[J]. Geomatics and Information Science of Wuhan University, 2008, 33(2): 132-135.

[9] 蔡國林, 李永樹, 劉國祥. 小波-維納組合濾波算法及其在InSAR干涉圖去噪中的應用[J]. 遙感學報, 2009, 13(1):129-136.CAI Guolin, LI Yongshu, LIU Guoxiang. Wavelet-wiener combined filter and its application on InSAR interferogram[J].Journal of Remote Sensing, 2009, 13(1): 129-136.

[10] Loffeld O, Nies H, Knedlik S, et al. Phase unwrapping for SAR interferometry: A data fusion approach by Kalman filtering[J].IEEE Transactions on Geoscience and Remote Sensing, 2009,47(4): 1197-1211.

[11] Martinez-Espla J J, Martinez-Marin T, Lopez-Sanchez J M.Using a grid-based filter approach to solve InSAR phase unwrapping[J]. IEEE Transactions on Geoscience and Remote Sensing letter, 2008, 5(2): 147-151.

[12] Lee J S, Papathanassiou K P, Ainsworth T L, el al. A new technique for noise filtering of SAR interferogram phase images[J]. IEEE Transaction Geoscience and Remote Sensing,1998, 36(5): 1456-1465.

[13] 尹宏杰, 李志偉, 丁曉利, 等. InSAR干涉圖最優(yōu)化方向融合濾波[J]. 遙感學報, 2009, 13(6): 1099-1105.YIN Hongjie, LI Zhiwei, DING Xiaoli, et al. Optimal integration-based adaptive direction filter for InSAR interferogram[J]. Journal of Remote Sensing, 2009, 13(6):1099-1105.

[14] Mashaly A S, AbdElkawy E E F, Mahmoud T A. Speckle noise reduction in SAR images using adaptive morphological filter[C]//Intelligent Systems Design and Applications ISDA 2010 10th International Conference, Publisher: IEEE, 2010,41(3): 260-265.

[15] 王興旺, 易輝偉, 朱建軍, 等. 基于有選擇保邊緣平滑法的SAR干涉圖像濾波算法[J]. 工程勘察, 2008(11): 50-53.WANG Xingwang, YI Huiwei, ZHU Jianjun, et al. A new phase noise reduction for INSAR interferogram based on edge-preservation smoother[J]. Geotechnical Investigation &Surveying, 2008(11): 50-53.

[16] Martinez-Espla J J, Martinez-Marin T, Lopez-Sanchez J M. An optimized algorithm for insar phase unwrapping based on particle filtering, matrix pencil, and region-growing techniques[J]. IEEE Geoscience and Remote Sensing Letters,2009, 6(4): 835-839.

[17] 郭交, 李真芳, 劉艷陽, 等. 一種InSAR干涉相位圖的自適應濾波算法[J]. 西安電子科技大學學報, 2011, 38(4): 77-88.GUO Jiao, LI Zhengfang, LIU Yanyang, et al. New adaptive noise suppressing method for interferometric phase images[J].Journal of Xidian University, 2011, 38(4): 77-88.

[18] 李佳, 李志偉, 丁曉利, 等. 強噪聲 SAR干涉圖的等值線中值: Goldstein二級濾波[J]. 遙感學報, 2011, 15(4): 758-764.LI Jia, LI Zhiwei, DING Xiaoli, et al. Filtering strong noisy synthetic aperture radar interferogram with integrated contoured Median and Goldstein two-step filter[J]. Journal of Remote Sensing, 2011, 15(4): 758-764.

[19] Buemime, M E, Jacobo J, Mejail M. SAR image processing using adaptive stack filter[J]. Pattern Recognition Letters, 2010,31(4): 307-314.

[20] Lopez-Martinez C, Fabregas X, Pipia L. Forest parameter estimation in the Pol-InSAR context employing the multiplicative-additive speckle noise model[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66(5): 597-607.

[21] 廖明生, 林琿, 張祖勛, 等. InSAR干涉條紋圖的復數(shù)空間自適應濾波[J]. 遙感學報, 7(2): 98-105.LIAO Mingsheng, LIN Hui, ZHANG Zuxun, et al. Adaptive algorithm for filtering interferometric phase noise[J].Journal of Remote Sensing, 2003, 7(2): 98-105.

[22] Franceschetti G, Lanari R. Synthetic aperture radar processing[M]. Florida: CRC Press, 1999: 185-190.