一種新的干涉相位圖局部自適應濾波方法

薛海偉 馮大政

?

一種新的干涉相位圖局部自適應濾波方法

薛海偉*馮大政

(西安電子科技大學雷達信號處理國家重點實驗室 西安 710071)

為了有效提高對InSAR干涉相位噪聲的抑制性能并充分保持干涉相位圖細節信息，該文提出一種基于局部地形相位補償和各向異性高斯濾波函數(AGF)的自適應復相位濾波方法。該方法首先利用局部頻率估計方法補償地形相位，以便于消除局部地形相位對濾波窗口內干涉相位的不利影響。然后，構造了尺度和方向自適應的AGF，并對同分布樣本進行局部加權的方向濾波。這里，AGF尺度隨相干系數等級自適應變化：在低相干區域，采用的大尺度AGF能夠充分地抑制相位噪聲；在高相干區域，采用的小尺度AGF能更好地保持相位細節信息。AGF方向根據最大加權相干積累準則確定，以選取同分布的濾波樣本估計中心像素相位值。實驗結果表明，與多種濾波方法相比，該文方法在減少干涉相位圖殘點和保持條紋邊緣等方面均具有更好的性能。

干涉合成孔徑雷達；相位噪聲濾波；相位噪聲模型；頻率估計；各向異性高斯濾波

1 引言

干涉合成孔徑雷達(InSAR)測量技術利用從不同位置不同視角觀測同一場景獲得的兩幅相干SAR圖像來生成復干涉相位圖，并進行高程重建得到觀測場景的數字高程圖(DEM)[1]。由于系統噪聲及各種去相干因素的影響，干涉相位圖中存在明顯的噪聲，嚴重影響了相位展開、高程重建等后續處理的性能。因此，干涉相位濾波是InSAR高程重建過程中的一個關鍵步驟，對于降低相位展開難度、提高高程測量精度都具有重要的意義。

目前，國內外提出了多種干涉相位濾波算法，其中多視濾波器[7]被認為是在最大似然意義下最優的。但是，多視濾波器常采用固定的濾波窗口，在相位細節信息保持和條紋密集區域濾波效果之間很難取得平衡。Goldstein等人[8]針對小塊干涉圖頻譜，提出了一種在保持相位譜情況下將幅度譜進行冪運算，實現條紋連續性保持的自適應濾波器，但因該方法取固定的冪指數，其噪聲抑制效果欠佳。隨后，Baran等人[9]讓該算法中的冪指數隨相干系數自適應變化。文獻[10]根據干涉條紋的方向和密度調整濾波窗口的大小和方向，然后進行中值或均值濾波，但該方法未能實現自適應強度濾波。Trouvé等人[11]在線性相位補償的基礎上采用較大的濾波窗口對殘余相位進行復多視處理，減少了殘點數目。在文獻[11]的工作基礎上，Cai等人[12]根據相干積累準則提出一種窗口大小自適應的局部頻率估計方法，對修正后的相位進行復均值濾波。但這兩種方法[11,12]都沒有進行樣本點的選取，濾波性能欠佳。利用干涉相位條紋的方向信息，Lee等人[13]從16個方向性濾波窗口中選擇具有最大復均值的方向窗口，并將當前窗內的相位平均值與中心像素相位觀測值進行線性組合得到中心像素相位的估計；該方法濾波精度較高，且避免了對高相干區域的過度濾波，但選擇的方向窗口可能不與條紋方向完全平行，從而導致條紋出現斷裂，且在濾波之前需進行局部相位預展開，這會引起較大的誤差。Wu等人[14]根據局部條紋頻率確定局部條紋方向，取方向窗與條紋方向完全平行，先利用插值計算濾波窗內的相位值，然后采用文獻[13]方法實現自適應濾波；該方法能夠更有效地保持條紋邊緣，且不需要局部相位預展開，獲得了良好的濾波性能。文獻[15]通過限定濾波樣本值選取范圍和去除相位外點等措施改進了文獻[13]濾波方法。Suo等人[16]采用自適應尺寸的濾波窗口、相位坡度補償以及濾波參數修正等措施進一步改進了文獻[8]方法，在保持相位邊緣連續性的同時提高了噪聲抑制性能。文獻[17]首先利用分塊后的相位圖頻譜和半功率門限估計主要地形相位，然后采用簡化的文獻[13]濾波形式對殘余相位進行濾波。李錦偉等人[18]首先借助外部數字模型反演干涉條紋并利用方向性窗口估計相干系數，對樣本進行最小噪聲方差意義下的最優加權濾波，取得了較好的濾波效果。

為了在有效抑制噪聲的同時避免破壞條紋信息，并選取同分布濾波樣本進行自適應強度的濾波，本文提出了一種基于局部地形相位補償和各向異性高斯濾波函數(AGF)的自適應復相位濾波方法。該方法利用局部地形相位補償獲得同分布樣本，并采用方向和尺度自適應的AGF實現濾波樣本選取和自適應強度濾波，可以充分抑制低相干區域的相位噪聲，并且能保持高相干區域的相位細節信息，即使在條紋密集區域，也能夠兼顧噪聲抑制和細節信息保持。最后采用仿真和實測數據對本文方法的性能進行了驗證。

2 干涉相位模型

InSAR中兩幅SAR圖像的相干性直接影響高程測量的精度，而相干性受到系統噪聲、基線去相干、時間去相干和圖像配準誤差等因素的影響。衡量相干性大小的相干系數定義為

(3)

實際中由于獲得SAR圖像數量的限制，通常對濾波窗口內樣本進行空間平均來實現多視處理。但多視處理要求局部相位滿足同分布，且抑制相位噪聲的同時會降低圖像的分辨率。由于復雜地形的影響，相位圖不同區域的相干性不同，且局部濾波窗口內的濾波樣本不滿足同分布的條件，特別在地形變化劇烈的地區。因此在相位濾波之前需要首先補償地形引起的相位，并選取同分布濾波樣本進行濾波，同時還需要避免對高相干性區域過度濾波導致的分辨率損失。

圖1 干涉相位噪聲標準差與相干系數和視數的關系

3 算法原理

由于InSAR測量中兩個雷達具有不同的觀測位置，獲得的干涉相位會隨著觀測場景的高度起伏和水平距離而變化[1]。本文方法基于局部平穩的假設，首先估計并補償局部地形坡度引起的相位，然后自適應地選取滿足同分布的濾波樣本。為此，將局部地形相位看作由以下兩部分組成：地形坡度引起的線性相位和偏離這一坡度的地形分量引起的相位，同時將式(3)擴展表示為

(6)

3.1 局部地形坡度相位估計

觀測場景內的地形是不斷變化的，但通常可以假定為局部平穩的。即把干涉圖劃分為許多較小的窗口后，每個窗口對應的局部觀測區域內的地形傾斜角是一致的。具有不同視角的雷達觀測傾斜角無變化的地形，所得到的干涉相位在距離方向和方位方向都只包含一個主要頻率，可表示為一個2維正弦模型[19]：

(8)

由于隨機分布噪聲和殘余地形相位的影響，窗口內的濾波樣本不完全滿足同分布。為了獲得更好的濾波性能，還需要采取進一步的措施選取濾波樣本。

3.2 基于AGF的自適應濾波

針對濾波窗口內樣本選取的問題，本文基于2維AGF生成尺度自適應于局部相干系數等級且方向連續變化的支撐區域，然后根據復相位平面內最大加權相干積累準則，采用2維極坐標傅里葉變換快速確定支撐區域的方向信息。方向自適應AGF選取的濾波樣本更加滿足同一分布，提高了相位濾波精度；尺度自適應于局部相干系數等級的AGF能夠有效抑制低相干區域的相位噪聲，同時保護高相干區域的相位條紋細節信息。

水平方向的AGF表示為

(10)

3.2.1 尺度自適應 為了避免對高相干區域過度濾波導致分辨率損失，需要根據相干系數等級來決定濾波的強弱。根據式(4)，高相干系數區域只需要相對較少的濾波樣本就可以達到與低相干區域類似的濾波效果，即高相干系數區域需要較小的支撐區域，因此我們根據局部相干系數等級調整支撐區域大小，實現自適應強度濾波。尺度為,的AGF等效濾波窗口大小為[20]，根據式(4)，和為

圖2 AGF支撐區域示意圖

圖3 自適應尺度AGF濾波效果圖

3.2.2 方向自適應 為了選取更加滿足同分布的濾波樣本，需要估計支撐區域的方向信息。考慮相位受到噪聲的影響及存在以為周期的模糊性，方向信息以最大加權相干積累作為準則自適應地確定。考慮到極坐標變換可以把難以匹配的角度旋轉變換為坐標軸方向的平移，本文利用2維極坐標傅里葉變換快速估計支撐區域的方向信息。首先，將殘余相位窗口中心移到坐標原點，并將直角坐標系中的點變換為極坐標系中的點，變換定義為

(13)

(15)

(17)

3.3 算法步驟

本文方法利用各向異性高斯函數生成各向異性的支撐區域，并根據局部噪聲分布自適應地選擇干涉相位濾波樣本實現局域加權的方向濾波。具體步驟如下：

(2)利用式(8)估計局部地形相位，將地形相位從原始干涉相位中去除得到殘余相位；

(3)根據式(11)確定AGF尺度；

(4)基于最大加權相干積累準則，利用2維極坐標傅里葉變換快速估計AGF的方向信息。根據式(16)生成2維高斯濾波函數；

(5)在復數平面上，根據式(17)對殘余相位進行濾波，得到像素的殘余相位估計值；

4 實驗結果

本節分別利用仿真數據和實測數據分析本文濾波方法的性能，并與改進Goldstein濾波器[8]、Lee自適應復濾波器[13]和文獻[11]提出的相位坡度補償方法進行了比較。

圖4所示為對一組仿真數據進行濾波的結果。圖4(a)為不含噪聲的仿真干涉相位圖，由MATLAB中的Peaks函數生成，干涉相位尺寸為像素。根據文獻[13]的加性噪聲模型，加入了標準差為1 rad的相位噪聲，如圖4(b)所示。改進Goldstein濾波器的分塊尺寸為32，重疊尺寸為14，采用窗口平滑頻譜幅度。Lee濾波器采用濾波窗口，16個方向窗口。本文方法中設置為0.2。圖4(c)至圖4(f)為各方法的濾波結果。從圖4所示濾波圖像視覺比較可以發現：改進Goldstein濾波器分塊之間出現不連續性，對噪聲抑制不夠充分，整個圖像內還存在大量的噪聲，性能最差。Lee濾波器在條紋密集處性能下降，多數條紋變得模糊甚至出現斷裂的現象，相位細節信息保持性能較差。文獻[11]方法濾波結果中明顯存在一定數量的孤立點噪聲，性能較好。本文方法不僅有效地抑制了噪聲，而且獲得了完整清晰的干涉相位圖，即使在條紋密集區域相位條紋邊緣也保持了良好的連續性，性能最好。為了更加直觀地比較各方法的濾波性能，我們采用相位偏差和殘點濾除率來衡量濾波算法的性能。圖5所示為上述各方法濾波結果與原始相位圖的相位偏差直方圖，濾波前后的殘點對比結果見表1。可以看出，本文方法濾波結果比其他方法更加接近原始相位圖，并且更有效地降低了殘點數目，具有更好的濾波降噪性能。

采用兩組實測數據以進一步評估本文濾波方法的性能，分別為：X-SAR錄取的Etna火山口的數據和TerraSAR-X衛星獲取的澳大利亞艾爾斯巖石數據，經過圖像配準和去平地相位處理后，生成的原始干涉相位圖如圖6所示，Etna數據和艾爾斯巖石數據的尺寸分別為像素和像素。利用本文方法與改進Goldstein濾波器[8]、Lee濾波器[13]和文獻[11]方法對圖6所示兩組數據進行濾波處理，實驗結果如表1和圖7所示。表1列出了上述各方法濾波后干涉相位圖中殘點的數目，通過對比可以看出，本文方法相對于其他方法更加有效地減少了殘點數目。從圖7(a)至圖7(d)所示對Etna數據濾波結果比較發現：改進Goldstein濾波結果中存在大量的殘點；Lee濾波器和文獻[11]方法能夠比較有效地減少殘點數目，但是在條紋密集區域性能不佳，條紋連續性被破壞，如圖7(b)、圖7(c)中白框內區域所示；本文方法獲得的結果圖像中殘點數目最小，并且更大程度地保持了相位條紋的細節信息，即使是在條紋密集區域，相位邊緣也保持了良好的連續性，如圖7(d)中白色圓框所示。對艾爾斯巖石數據進行濾波處理也得到了同樣的結果，如圖7(e)至圖7(h)所示。綜上，從濾波結果和對剩余殘點的統計結果可以看出，與其他幾種方法相比，本文方法能夠獲得更加清晰的干涉相位圖，在條紋密集區域性能優勢更加明顯，表明本文方法能夠在對噪聲進行有效抑制的同時充分保持干涉相位圖的細節信息。

圖4 仿真干涉相位濾波結果比較

圖5 相位偏差直方圖??????????????圖6 實驗采用的兩組數據

表1 各方法對仿真數據、Etna數據和艾爾斯巖石數據濾波后的殘點數目比較

濾波方法仿真數據Etna數據艾爾斯巖石數據 殘點數殘點減少(%)殘點數殘點減少(%)殘點數殘點減少(%) 原始干涉相位圖22598-97657-46494- 改進Goldstein濾波器[8]1272243.702670572.65 949779.57 Lee濾波器[13] 159892.931310586.58 288793.79 文獻[11]方法 89096.06 732592.50 397391.45 本文方法 9499.58 142598.54 81998.24

計算復雜度方面，改進Goldstein濾波器由于逐塊進行濾波，運行時間最短。Lee濾波器、文獻[11]方法和本文方法均是逐像素進行濾波，Lee濾波器需要定位方向窗和計算加權系數，文獻[11]方法僅估計局部相位坡度，本文方法需要估計局部地形相位和自適應AGF的參數，因此在局部窗口大小相同的條件下，本文方法計算復雜度高于改進Goldstein濾波器和文獻[11]方法，比Lee濾波器略高。表2所示為各方法的運行時間(CPU: 3.6 GHz，內存：32.0 G)。

表2 各方法運行時間

濾波方法運行時間(s) 仿真數據Etna數據艾爾斯巖石數據 改進Goldstein濾波器 6.859 26.699 23.206 Lee濾波器119.067465.891372.797 文獻[11]方法 56.254213.550187.089 本文方法132.114502.934435.183

5 結束語

本文提出了一種有效抑制干涉相位噪聲并同時保持干涉相位圖細節信息的自適應復相位濾波方法。與傳統空域濾波相比，本文方法具有以下優勢：其一，方法首先對局部地形斜面引起的相位進行估計并補償，降低了局部地形對干涉相位的影響，更有利于選取同分布樣本進行濾波；其二，AGF尺度自適應于相干系數提供不同程度的平滑效果[20]，具體來講：低相干區域采用大尺度AGF生成較大的有效支撐區域，提供了較大的平滑度，以確保相位噪聲得到足夠的抑制；而高相干區域采用小尺度AGF生成較小有效支撐區域，能夠更好地保持相位細節信息。其三，AGF方向根據最大加權相干積累準則確定，使得AGF支撐區域內選取的濾波樣本更加滿足同分布，提高了濾波精度，同時AGF方向可利用2維極坐標傅里葉變換快速確定。仿真和實測數據處理結果驗證了本文方法的性能。

參考文獻

[1] ROSEN P A, HENSLEY S, JOUGHIN I R,. Synthetic aperture radar interferometry[J]., 2000, 88(3): 333-382. doi: 10.1109/5.838084.

[2] LIN X, LI F F, MENG D D,. Nonlocal SAR interferometric phase filtering through higher order singular value decomposition[J]., 2015, 12(4): 806-810. doi: 10.1109/LGRS. 2014.2362952.

[3] CAO M Y, LI S Q, WANG R,. Interferometric phase denoising by median patch-based locally optimal wiener filter[J]., 2015, 12(8): 1730-1734. doi: 10.1109/LGRS.2015.2422788.

[4] LI H Y, SONG H J, WANG R,. A modification to the complex-valued MRF modeling filter of interferometric SAR phase[J]., 2015, 12(3): 681-685. doi: 10.1109/LGRS.2014.2357449.

[5] LI J W, LI Z F, BAO Z,. Noise filtering of high- resolution interferograms over vegetation and urban areas with a refined nonlocal filter[J]., 2015, 12(1): 77-81. doi: 10.1109/LGRS.2014. 2326462.

[6] SONG R, GUO H D, LIU G,. Improved Goldstein SAR interferogram filter based on adaptive-neighborhood technique[J]., 2015, 12(1): 140-144. doi: 10.1109/LGRS.2014.2329498.

[7] SEYMOUR M S and CUMMING I G. Maximum likelihood estimation for SAR interferometry[C]. International Geoscience and Remote Sensing Symposium, 1994, 4: 2272-2275. doi: 10.1109/IGARSS.1994.399711.

[8] GOLDSTEIN R M and WERNER C L. Radar interferogram filtering for geophysical applications[J]., 1998, 25(21): 4035–4038. doi: 10.1029/1998GL900033.

[9] BARAN I, STEWART M P, KAMPES B M,. A modification to the Goldstein radar interferogram filter[J]., 2003, 41(9): 2114-2118. doi: 10.1109/TGRS.2003.817212.

[10] FU S H, LONG X J, YANG X,. Directionally adaptive filter for synthetic aperture radar interferometric phase images[J]., 2013, 51(1): 552-559. doi: 10.1109/TGRS.2012.22. 2911.

[11] TROUVE E, NICOLAS J M, and MAITRE H. Improving phase unwrapping techniques by the use of local frequency estimates[J]., 1998, 36(6): 1963-1972. doi: 10.1109/36.729368.

[12] CAI B, LIANG D N, and DONG Z. A new adaptive multiresolution noise-filtering approach for SAR interferometric phase images[J]., 2008, 5(2): 266-270. doi: 10.1109/ LGRS.2008.915942.

[13] LEE J S, PAPATHANASSIOU K P, AINSWORTH T L,. A new technique for noise filtering of SAR interferometric phase images[J]., 1998, 36(5): 1456-1465. doi: 10.1109/36. 718849.

[14] WU N, FENG D Z, and LI J X. A locally adaptive filter of interferometric phase images[J]., 2006, 3(1): 73-77. doi: 10.1109/ LGRS.2005.856703.

[15] CHAO C F, CHEN K S, and LEE J S. Refined filtering of interferometric phase from InSAR data[J]., 2013, 51(12): 5315-5323. doi: 10.1109/TGRS.2012.2234467.

[16] SUO Z Y, ZHANG J Q, LI M,. Improved InSAR phase noise filter in frequency domain[J]., 2016, 54(2): 1185-1195. doi: 10.1109/TGRS.2015.2476355.

[17] WANG Q S, HUANG H F, YU A X,. An efficient and adaptive approach for noise filtering of SAR interferometric phase images[J]., 2011, 8(6): 1140-1144. doi: 10.1109/LGRS.2011. 2158289.

[18] 李錦偉, 李真芳, 劉艷陽, 等. 一種相干系數加權的最優干涉相位濾波[J]. 西安電子科技大學學報, 2014, 41(2): 25-31. doi: 10.3969/j.issn.1001-2400.2014.02.005.

LI J W, LI Z F, LIU Y Y,. Coherence-weighted optimum interferometric phase filtering method[J]., 2014, 41(2): 25-31. doi: 10.3969/ j.issn.1001-2400.2014.02.005.

[19] SUO Z Y, LI Z F, and BAO Z. A new strategy to estimate local fringe frequencies for InSAR phase noise reduction[J]., 2010, 7(4): 771-775. doi: 10.1109/LGRS.2010.2047935.

[20] SHUI P L and ZHANG W C. Noise-robust edge detector combing isotropic and anisotropic Gaussian kernels[J]., 2012, 45(2): 806-820. doi: 10.1016/ j.patcog.2011.07.020.

薛海偉： 男，1985年生，博士生，研究方向為干涉SAR數據處理、SAR成像.

馮大政： 男，1959年生，教授，博士生導師，研究方向為雷達成像、陣列信號處理、盲信號處理、神經網絡等.

New Locally Adaptive Method for InSAR Phase Noise Filtering

XUE Haiwei FENG Dazheng

(,,710071,)

In order to effectively suppress the noise of InSAR phase images and preserve the detailed fringe information, an adaptive phase filtering method based on local slope compensation and the Anisotropic Gaussian Filter (AGF) is proposed. Firstly, the topography-induced phase is approximately measured by local frequency estimation and removed from the original phase to eliminate the effect of the terrain topography. Secondly, the AGF with adaptive scale and orientation is developed to directionally filter out the noisy phase for the pixels with more homogeneous phase values. The scale of the AGF varies adaptively with the local coherence: a large-scaled AGF can better smooth the noise of low coherence areas, whereas a small-scaled AGF can better preserve the phase details of high coherence areas. Moreover, the orientation angle of the AGF is fast determined to select the identically distributed samples according to the maximum weighted coherent summation principle. The experimental results obtained via simulated and real data show that compared with commonly used filters, the proposed method achieves better performance in terms of residue reduction and fringe preservation.

Interferometric SAR (InSAR); Phase noise filtering; Phase noise modeling; Frequency estimation; Anisotropic Gaussian Filter (AGF)

TN957.51

A

1009-5896(2016)12-3085-08

10.11999/JEIT161013

2016-09-30；改回日期：2016-11-25；

2016-12-14

薛海偉 xuehw001@163.com

國家自然科學基金(61271293)

The National Natural Science Foundation of China (61271293)