InSAR礦區形變監測的邊緣保持-Goldstein組合濾波方法

易輝偉 ，朱建軍 ，李 健 ，李 佳 ，嚴新生，陳樹泉

(1.中南大學 地球科學與信息物理學院，長沙 410083；2.湖南省普通高校精密工程測量及形變災害監測重點實驗室，長沙 410083；3.廣東省江門市勘察院，江門 529040)

InSAR技術由于其監測范圍大、精度高、周期短及無接觸式等特點，已廣泛應用于礦區地表形變的監測中[1-2]，然而 InSAR的精度和可靠性在很大程度依賴于由兩幅或多幅SAR圖像所形成干涉圖的質量[3]。受時空失相關、熱噪聲以及大氣等因素影響，干涉圖相位一般存在強噪聲，為解纏順利及精確估計地形形變，必須對SAR干涉圖進行濾波[4]。

目前，SAR干涉圖濾波方法主要分為頻域濾波法和空域濾波法兩類。頻域濾波法有 Goldstein[5]、小波和卡爾曼等基于傅里葉變換對影像譜分量進行處理的濾波方法。Goldstein為頻域經典濾波法，由于去噪能力強且步驟簡單而得到廣泛應用，如Gamma和Doris等許多國際商業或免費 SAR數據處理軟件都采用該方法。Goldstein濾波法將圖像分解到頻率空間，利用相干斑噪聲和信號的頻譜特性不同，即噪聲在頻譜空間屬于寬帶信號且幅值低，而信號則屬于窄帶信號且幅值高，通過在頻率域實施平滑，使得屬于噪聲的寬帶信號幅值均勻而屬于信號的窄帶信號幅值相對增大，從而達到去噪的目的。后來發展的LZW[6]、信噪比[7]及相位偏差[8]等方法都根據噪聲強度不同程度地改善了濾波因子。該方法難以檢測到圖像的邊緣特征，在去除噪聲的同時，也損失了圖像邊緣細節信息。小波去噪有矢量分離式[9]和小波中值濾波[10]等方法，可在一定程度保留細節，但因小波閾值為人工選取，因此，不便于使用，小波-維納組合法[11]進行了較好的改進。卡爾曼[12]和格網值濾波法[13]都是通過對相位值進行估計得到濾波。空域濾波法[14-22]有Lee濾波法、數學形態法、最優化融合法、堆棧法、區域增長法、乘性加性噪聲模型法、中值-自適應濾波、等值線法及邊緣保持法等，是一種在圖像空間借助模板卷積來實現濾波的方法。為避免纏繞相位的影響，通常將InSAR干涉圖的相位圖和灰度圖對應位置的值構成的復數分解為實部和虛部，分別對實部和虛部進行濾波，最后合并得到濾波后的相位圖。Lee濾波法[14]為典型的空域濾波，其利用局部統計特性進行SAR圖像去噪,但目視效果較差，邊緣保持只在同質區比較有效；數學形態法[15]在Lee濾波法中的非同質區域采取膨脹腐蝕增強了邊緣。Lee濾波法基于乘性噪聲模型估計像元的真實值，由于該模型未考慮加性噪聲問題，并且需要求坡度值、檢測條紋邊緣及局部解纏，因而使用受限。針對此缺點，尹宏杰等[16]提出的最優化融合方法根據干涉圖的相干性來選擇融合方向窗口的數量，增強了Lee算法的穩健性。堆棧法[17]將圖像分成若干二進制圖像，利用布爾判斷式去噪再重新組合；區域增長法[18]基于種子點向外增長直至邊緣，取決于邊緣判斷條件；乘性加性噪聲模型法[19]認為對角元素為乘性模型而非對角元素為加性模型，該方法需進一步驗證。以上 3種方法缺乏自適應性。中值-自適應濾波法[20]適用于寬條紋去噪，對于密集的細條紋容易出現混迭現象；等值線濾波法[21]中的等值線窗口的提取受噪聲程度和條紋率的影響。2008年，王興旺等[22]提出了邊緣保持濾波法，該法能夠較好地保留條紋細節。

SAR干涉圖濾波要達到的濾波效果需在盡量去除相位噪聲的同時，保留圖像的條紋信息。一般情況下，頻域濾波法的去噪效果好而空域濾波法的細節信息保持能力強，因此，本文作者以偽相干值為橋梁，并用具有局部統計特性的相位標準偏差改善Goldstein方法的濾波因子，將邊緣保持和Goldstein這兩種濾波方法組合起來，形成邊緣保持-Goldstein組合濾波方法，該方法充分結合頻域濾波法和空域濾波法的優點，能有效控制噪聲并保持條紋細節，且適用于條紋密集的干涉圖。

1 基于偽相干值的邊緣保持-Goldstein噪聲控制原理

1.1 邊緣保持濾波法

邊緣保持濾波法以方差作為各個鄰域灰度均勻性的量度。若像素鄰域含有尖銳的邊緣，則其灰度方差較大；而不含邊緣或灰度均勻的鄰域，方差則較小。通過找到最小方差的鄰域，對鄰域內像素取均值，可達到既消除部分噪聲、又不破壞鄰域邊界細節的效果。考慮到條紋邊界形狀的多樣性，對圖像上任意像素(x,y)采用如圖1所示的模板，其中包括一個 3×3的正方形，4個五邊形和4個六邊形共9個鄰域。由于五邊形和六邊形在(x , y)處都有銳角，因此，即使像素(x,y)位于一個復雜形狀的邊緣處，也能找到均勻的鄰域。

1.2 Goldstein濾波法

Goldstein 濾波法主要通過對圖像頻譜進行平滑處理而達到去噪目的，其算法如下：

式中：Z(u,v)為對圖像進行傅里葉變換得到的頻譜值；S{ }是平滑函數；Z’(u,v)是平滑處理后的頻譜值；α的取值范圍為0～1，濾波的強度隨α的增大而增大，α為0時不濾波。BARAN等[23]將α改進為為有效窗口(通常取32×32像素窗口為滑動窗口，其中間的4行4列4×4為有效窗口)內相干值平均值，這種自適應的Goldstein方法濾波效果更好。

經邊緣保持法濾波后的干涉圖，其圖像質量可用偽相干值[24]來衡量，其計算公式如下：

圖1 邊緣保持法對應模板Fig.1 Templates for selective edge-preservation filter: (a)One square; (b)Four pentagons; (c)Four hexagons

式中：N是計算窗口內像素的個數；φ(i,j)是干涉圖相位值。偽相干值由干涉圖統計得到，隨著相位噪聲強度的降低而增大，最大值不超過 1。偽相干值是邊緣保持濾波結果的反映，濾波后圖像相干性得到增強。此外，相位標準偏差具有較強的局部統計特性，表征圖像噪聲的強度，且與圖像相干性存在如下關系[6]：

式中：σφ為相位標準偏差；γ為干涉圖相干值；φ為干涉圖相位值；φ0為相位的期望值；Li2(·)為一種 Euler算法，表示為

1.3 邊緣保持-Goldstein組合濾波方法

綜合考慮具體兩個因素，提出用式(5)所示改進Goldstein法的濾波因子，使其對強噪聲區進行強濾波，對弱噪聲區進行弱濾波，從而使 Goldstein 方法具有更強的局部自適應能力。

式中：為有效窗口內偽相干性平均值；為有效窗口內相位標準偏差；σmax為的最大值。邊緣保持-Goldstein組合濾波方法首先采用邊緣保持法對干涉圖進行初步濾波，在去除一定噪聲的同時，保持每個像素與其鄰域的關系；再采用式(5)進行Goldstein濾波，由于和均能有效反映噪聲分布，本文作者采用控制濾波強度，不僅濾波效果好，而且條紋信息能得到保持；最后按照式(6)將 Goldstein濾波結果轉換到空間域。

式中：r為距離向向量；a為方位向向量；F-1為傅立葉反變換。

邊緣保持-Goldstein組合濾波的具體步驟如下：

1)分別提取干涉圖的實部和虛部，對于每個像素，分別計算上述9個對應模板的均值和方差；

2)選取方差最小的模板，將該模板的均值作為該像素的實部和虛部輸出值；

3)將兩個新矩陣重新組合，生成干涉圖；

4)取32×32像素的滑動窗口，每次滑動4個像素，滑動窗口中間的4行4列為有效窗口，計算有效窗口的和相位標準偏差均值；

5)在有效窗口中對復數進行離散傅立葉變換；

6)對取出的部分按改進后的濾波因子進行平滑處理；

7)將平滑后的圖像反變換到空間域。

2 模擬數據驗證

圖2 模擬干涉圖濾波結果Fig.2 Results of simulated interferograms by different filters:(a)Simulated phase diagram; (b)Corrupted phase diagram;(c)Once Goldstein filter; (d)Twice Goldstein filter; (e)Lee filter; (f )Median- adaptive filter; (g)This method

利用模擬數據，去噪結果可與無噪圖直接比較。采用多分形技術模擬數字地面模型(DEM)，再由模擬的DEM根據ERS1/2的成像幾何參數和給定垂直基線長度(本實驗中為200 m)模擬出纏繞的相位，最后根據相位噪聲與視數及相干性的關系模擬出給定視數的相位噪聲圖，該圖既含乘性噪聲又含加性噪聲。圖2(a)所示為模擬的512×512像素原始相位圖；圖2(b)所示為視數為2的模擬含噪相位圖，對其分別用一次、二次Goldstein濾波法、Lee濾波法、中值-自適應濾波法及本研究所提出的新方法進行濾波，結果分別為圖2(c)～(g)所示。從目視效果看，一次 Goldstein濾波去噪效果不好，很多區域仍殘存較多噪聲，而有些區域條紋細節未能有效保持；二次Goldstein過度濾波，條紋高度平滑，丟失了更多邊緣信息，如圖2(d)中右上、左下及右下矩形區域所示；Lee濾波噪聲分布不均勻，有些區域仍滯留了噪聲，如圖2(e)中左下矩形所示；中值-自適應濾波去噪均勻，邊緣保持良好，但仍有少量噪聲殘留并滲入相位，條紋比原始圖更粗糙。本文提出的方法根據邊緣保持濾波結果的偽相干值和相位標準偏差估計噪聲強度，并用調節Goldstein濾波強弱，在去除大量噪聲的同時較好地保護了條紋的曲折度，濾波結果更接近原始相位圖。

圖3所示為圖2對應的第420行(如圖2(g)黑線所示)剖面線圖。圖3(a)和(b)所示分別為不含噪聲和含噪聲的相位剖面線；圖3(c)～(g)所示分別為圖2中對應方法的相位剖面線。可見，一次Goldstein濾波剖面線與原始相位圖相比，仍有較多噪聲影響其相位值，局部有毛刺；二次 Goldstein剖面線損失了許多相位細節，且隨著濾波次數的增加，此現象更加明顯；Lee濾波在局部區域仍保留噪聲；中值-自適應濾波條紋保持較好，但仍有少量毛刺，相位受殘余噪聲影響略顯粗糙；本方法的剖面圖幾乎沒有毛刺，噪聲也大量減少，整條曲線比較平滑，不僅去噪效果良好，而且有效保持了相位細節，其剖面圖與原始圖更好地保持了一致。

圖3 模擬干涉圖濾波結果第420行剖面圖Fig.3 Cross sections of line 420 of simulated interferograms filtered by different filters: (a)Simulated phase; (b)Corrupted phase; (c)Once Goldstein filter; (d)Twice Goldstein filter;(e)Lee filter; (f)Median-adaptive filter; (g)This method

3 真實數據驗證

以意大利Marini礦區的一幅1 200×1 200像素的干涉雷達降軌圖像為實驗對象，如圖4所示，該圖由歐洲航天局ERS-2衛星于2000年9月6日和10月11日分別獲取(幅框號：2583, 軌道號：222)，其干涉對垂直基線為305 m，方位向經過5個像素的多視處理，分辨率約為25 m×25 m，其相干值最大為0.78，大多在0.2～0.5之間，屬于強噪聲高密度條紋干涉圖。

圖4所示為Marini礦區原始相位干涉圖。為便于分析，分別選取 300×300像素的平緩區A區和200×200像素的條紋密集區B區的濾波結果放大，結果分別如圖5和6所示。A和B區的相干值均值分別為0.42和0.31，屬強噪區。結果顯示，二次Goldstein濾波條紋的曲折度明顯減小，變得強直；Lee濾波和中值-自適應濾波均有殘留噪聲，但條紋保護較好；而本方法的去噪和細節保持均良好。特別是B區條紋密集且多不連續，幾乎被噪聲遮蔽。二次Goldstein濾波在強噪聲處仍殘留若干斑點，且原圖中可觀察到的條紋曲度消失，圖6(c)黑框所示處的條紋連接甚至與原圖不一致，表現為“過濾波”；Lee濾波和中值-自適應濾波在強噪聲處效果不夠明顯；本方法濾波使干涉條紋更加清晰，并檢索出部分被遮蔽的條紋，使斷裂條紋呈現“愈合”趨勢，原先不連續的條紋在一定程度上連續了，條紋曲度保持得很好，其保真性和自適應性表現更加突出。

圖4 Marini礦區原始相位干涉圖Fig.4 Original phase interferogram of Marini mining region

2.3 定量評價

對模擬圖和真實圖采用正殘差點數(Positive residual points, PRPs)、負殘差點數(Negative residual points, NRPs)、相位標準偏差(PSD)、相位梯度和值(SPD)及邊緣保持指數(EPI)等指標進行比較。

圖5 圖4中A區濾波結果放大圖Fig.5 Enlarged filtering results of zone A in Fig.4 : (a)Original phase diagram; (b)Once Goldstein filter; (c)Twice Goldstein filter;(d)Lee filter; (e)Median-adaptive filter; (f)This method

就單個指標而言，殘差點數、PSD及SPD都是越小越好，而EPI越接近1越好，如表1和2所列。表1所列為模擬相位圖濾波結果，由殘差點數可以看出，對于上述5種濾波方法，后4種方法均在較大程度上去除了噪聲，二次 Goldstein和本方法的去噪能力最強，殘差點改善率分別為 97.7%和 97.2%。本方法的PSD、SPD與二次Goldstein的接近，但二次Goldstein的 EPI僅為 0.48，表現為“過濾波”，本方法的 EPI為0.85，更好地保留了邊緣細節，具有較強條紋保持能力，表現為“適濾波”。

表2所列為真實相位圖A區和B區的濾波結果，在去噪方面，一次Goldstein、Lee濾波和中值-自適應濾波的殘差點改善率均有所降低，二次Goldstein和本方法的殘差點改善率仍保持在90%以上，表明這2種方法能保持很強的去噪能力。本方法在B區的殘差點改善率比一次Goldstein的提高了23.7%，且略高于二次Goldstein的，表明其在強噪聲條紋密集區自適應性更強。PSD和SPD均以二次Goldstein和本方法的為最小，而本方法在A區和B區仍具有較高的邊緣保持指數，分別為 0.77和 0.74，相比較而言，本方法對于條紋平滑區和條紋密集區都有很好的濾波效果。

圖6 圖4中B區濾波結果放大圖Fig.6 Enlarged filtering results of zone B in Fig.4 : (a)Original phase diagram; (b)Once Goldstein filter; (c)Twice Goldstein filter;(d)Lee filter; (e)Median-adaptive filter; (f)This method

表1 模擬干涉圖濾波結果比較Table 1 Comparison among filtering results of simulated interferogram

表2 Marini礦區干涉圖濾波結果比較Table 2 Comparison among filtering results of InSAR interferogram over Marini mining region

4 結論

1)單次Goldstein濾波去噪能力不強，重復濾波難以保留條紋細節信息。

2)本研究提出的邊緣保持-Goldstein組合濾波方法結合空域和頻域濾波的優點，運用邊緣保持法既進行了初步濾波，又保護了邊緣，將其結果作為Goldstein的輸入，并根據反映其初步結果噪聲強度的偽相干值和相位標準偏差來改善 Goldstein的濾波因子，使得弱噪區弱濾波，強噪區強濾波，濾波強度得到有效調整。

3)用模擬數據和真實數據對經典的 Goldstein、Lee、中值-自適應方法和本文作者提出的新方法進行對比研究表明，無論從目視解譯還是定量指標評價，邊緣保持-Goldstein組合濾波法優于其他方法。模擬圖和真實圖的殘余點數去除率均達到90%以上，在真實圖的強噪區甚至高于二次Goldstein，而表示條紋細節的邊緣保持指數達到0.74，高于其他方法的。此法在去噪和條紋細節保持之間能得到更好的平衡，具有更強的局部自適應能力。

[1]陳國滸, 劉云華, 單新建.PSInSAR技術在北京采空塌陷區地表形變測量中的應用探析[J].中國地質災害與防治學報,2010, 21(2): 59-63.CHEN Guo-hu, LIU Yun-hua, SHAN Xin-jian.Application of PSInSAR technique in the deformation monitoring in mining collapse areas in Beijing[J].The Chinese Journal of Geological Hazard and Control, 2010, 21(2): 59-63.

[2]朱建軍, 邢學敏, 胡 俊, 李志偉.利用InSAR技術監測礦區地表形變[J].有色金屬學報, 2011, 21(10): 2564-2575.ZHU Jian-jun, XING Xue-min, HU Jun, LI Zhi-wei.Monitoring of ground surface deformation in mining area with InSAR technique[J].The Chinese Journal of Nonferrous Metals, 2011,21(10): 2564-2575.

[3]李 陶.重復軌道星載SAR差分干涉監測地表形變研究[D].武漢: 武漢大學測繪學院, 2004: 37-45.LI Tao.Study on surface deformation monitoring by repeat pass spaceborne SAR differential interferograms [D].Wuhan: Wuhan University, 2004: 37-45.

[4]廖明生, 林 琿.雷達干涉測量—原理與信號處理基礎[M].北京: 測繪出版社, 2003: 84-103.LIAO Ming-sheng, LIN Hui.Synthetic aperture radar interferometry—Principle and signal processing[M].Beijing:Surveying and Mapping Press, 2003: 84-103.

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

[6]LI Zhi-wei, DING Xiao-li, HUANG Chen, ZHU Jian-jun, CHEN Yan-li.Improved filtering parameter determination for the Goldstein radar interferogram filter[J].ISPRS Journal of Photogrammetry and Remote Sensing, 2008, 63: 621-634.

[7]孫 倩, 朱建軍, 李志偉, 尹宏杰.基于信噪比的InSAR干涉圖自適應濾波[J].測繪學報, 2009, 38(5): 437-442.SUN Qian, ZHU Jian-jun, LI Zhi-wei, YIN Hong-jie.A new adaptive InSAR interferogram filter based on SNR[J].Acta Geodaetica et Cartographica Sinica, 2009, 38(5): 437-442.

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

[9]靳國旺, 韓曉丁, 賈 博, 李 豪.InSAR干涉圖的矢量分離式小波濾波[J].武漢大學學報: 信息科學版, 2008, 33(2):132-135.JIN Guo-wang, HAN Xiao-ding, JIA Bo, LI Hao.Filtering for InSAR interferograms by vector decomposing and wavelet transformation[J].Journal of Wuhan University: Geomatics and Information Science, 2008, 33(2): 132-135.

[10]汪魯才, 王耀南, 毛六平.基于小波變換和中值濾波的InSAR干涉圖濾波方法[J].測繪學報, 2005, 34(2): 108-112.WANG Lu-cai, WANG Yao-nan, MAO Liu-ping.An algorithm of interferometric phase filter of InSAR based on wavelet analysis and median filter algorithm[J].Acta Geodaetica et Cartographica Sinica, 2005, 34(2): 108-112.

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

[12]LOFFELD O, NIES H, KNEDLIK S, WANG Y.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.

[13]MARTINEZ-ESPLA J J, MARTINEZ-MARIN T, LOPEZSANCHEZ 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.

[14]LEE J S, PAPATHANASSIOU K P, AINSWORTH T L, GRUNE S M R, REIGBER A.A new techniques for noise filtering of SAR interferogram phase images[J].IEEE Transaction Geoscien ce and Remote Sensing, 1998, 36(5): 1456-1465.

[15]MASHALY A S, ABDELKAWY E E F, MAHMOUD T A.Speckle noise reduction in SAR images using adaptive morphological filter[C]//Proceedings of the 10th International Conference on Intelligent Systems Design and Applications ISDA 2010.New York: IEEE, 2010: 260-265.

[16]尹宏杰, 李志偉, 丁曉利, 蔣 彌, 孫 倩, 王 平.InSAR干涉圖最優化方向融合濾波[J].遙感學報, 2009, 13(6):1099-1105.YIN Hong-jie, LI Zhi-wei, DING Xiao-li, JIANG Mi, SUN Qian,WANG Ping.Optimal integration-based adaptive direction filter for InSAR interferogram[J].Journal of Remote Sensing, 2009,13(6): 1099-1105.

[17]BUEMI M E, JACOBO J, MEJAIL M.SAR image processing using adaptive stack filter[J].Pattern Recognition Letters, 2010,31(4): 307-314.

[18]MARTINEZ-ESPLA J J, MARTINEZ-MARIN T, LOPEZSANCHEZ 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.

[19]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.

[20]廖明生, 林 琿, 張祖勛, 楊 文. InSAR 干涉條紋圖的復數空間自適應濾波[J].遙感學報, 2003, 7(2): 98-105.LIAO Ming-sheng, LIN Hui, ZHANG Zu-xun, YANG Wen.Adaptive algorithm for filtering interferometric phase noise[J].Journal of Remote Sensing, 2003, 7(2): 98-105.

[21]YU Qi-feng, YANG Xia, FU Si-hua, LIU Xiao-lin, SUN Xiang-yi.An adaptive contoured window filter for interferometric synthetic aperture radar[J].IEEE Geoscience and Remote Sensing Letters, 2007, 4(1): 23-26.

[22]王興旺, 易輝偉, 朱建軍, 張長書, 胡 俊.基于有選擇保邊緣平滑法的InSAR干涉圖像濾波算法[J].工程勘察, 2008(11):50-53.WANG Xing-wang, YI Hui-wei, ZHU Jian-jun, ZHANG Chang-shu, HU Jun.A new phase noise reduction for InSAR interferogram based on edge-preservation smoother[J].Geotechnical Investigation and Surveying, 2008(11): 50-53.

[23]BARAN L, STEWART M P, KAMPES B M, PERSKI Z, LILI P.A modification to the Goldstein radar interferogram filter[J].IEEE Transaction on Geoscience and Remote Sensing, 2003,41(9): 2114-2118.

[24]GHIGLIA D C, PRITT M D.Two-dimensional phase unwrapping: Theory, algorithms, and software[M].New York:John Wiley & Sons, 1998: 53-102.