一种改进的InSAR干涉图复数空间自适应滤波

易辉伟^{1, 2}，朱建军^{1, 2}，陈建群^{1, 2}，齐艳妮¹，李佳^1,2

(1. 中南大学 地球科学与信息物理学院，湖南 长沙，410083；

2. 湖南省普通高校精密工程测量及形变灾害监测重点实验室，湖南 长沙，410083)

摘要：阐述合成孔径干涉雷达(InSAR)干涉条纹图的频率域和空间域滤波的算法及其优点与不足，着重分析复数空间自适应滤波的算法原理，提出2点改进：增加45°方向和135°方向的梯度计算及根据干涉图相干系数确定条纹边缘的尺度函数，从而进一步改善权系数。以模拟数据和我国香港及意大利Etna火山2地的真实干涉图数据为例，用本文方法与空间自适应滤波、Lee滤波、基于信噪比的Goldstein方法、最优化方向融合方法进行比较。实验结果表明：改进算法能够根据噪声强弱相应地控制滤波的强弱，既能更好地抑制相位噪声，又能较好地保持条纹边缘信息，模拟干涉图和香港干涉图残余点改善率分别提高18.6%和17.2%，表现出更强的局部自适应性。

关键词：合成孔径干涉雷达；干涉图；自适应；滤波；相干系数

中图分类号：TP237.3          文献标志码：A         文章编号：1672-7207(2013)02-0632-10

An improved adaptive algorithm for filtering InSAR interferogram in complex plane

YI Huiwei^{1, 2}, ZHU Jianjun^{1, 2}, CHEN Jianqun^{1, 2}, QI Yanni¹, LI Jia^{1, 2}

(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;

2. Key Lab of Precise Engineering Surveying & Deformation Hazard Monitoring of Hunan Province, Changsha 410083, China)

Abstract: Algorithms and their advantages and disadvantages of filters for InSAR fringe images in frequence and space domains were described. The algorithm theory of the complex space domain adaptive filter was discussed. An improvement was made in two points, one was adding the gradient computation in 45° and 135° directions, and the other was the determination of scaling function for fringe edge based on interferogram coherence coefficient, thus modifying weight coefficients. According to the simulation data and actual interferograms of Hongkong and Etna volcano in Italy, the proposed method was compared with the complex space domain adaptive filter, Lee filter, Goldstein based on SNR (signal-noise-ratio), and adaptive direction-dependent based on optimal integration. The experiment results demonstrate that this improved method has capability of tuning the filtering according to noise intensity. It performs well not only in noise reduction but also in preservation of the sharp sawtooth profile of fringe, with an improvement of 18.9% and 17.2% of residual points reduction in simulate and actual experiment respectively, showing stronger local adaptability.

Key words: interferometric synthetic aperture radar; interferogram; adaptive; filter; coherence coefficient

合成孔径干涉雷达(InSAR)能够全天候、全天时、无接触式地获取地球表面形变的精密信息，已广泛用于地震、滑坡、冰移、火山及矿山沉降等现象监测中。然而，受热噪声、时间去相关、空间去相关等因素影响^[1]，INSAR干涉图像不可避免地存在大量噪声，难以解缠以至于无法得到数字地面高程模型及形变结果。因此，如何更好地对干涉图进行去噪，成为INSAR实际应用中亟待解决的问题^[2]。干涉图滤波方法分为频率域和空间域2类，频率域滤波是将图像的滑动小窗口(通常为32像素×32像素)内像素通过傅里叶变换为功率谱，并对该功率谱进行平滑再通过反傅里叶变换到图像空间。频率域滤波方法的代表为Goldstein^[3]滤波。Baran等^[4]应用相干值改善Goldstein的滤波因子，使其更具有自适应性。Li等^[5]引入干涉图的视数及相位标准偏差来调整滤波因子，使干涉图相位信息保持和算法稳健性得到有效改善；孙倩等^[6]应用信噪比控制滤波强度，对低信噪比区域强滤波，高信噪比区域弱滤波。于晓歆等^[7]应用最大最小值法对相位标准偏差进行归一化处理后采用其指数作为滤波因子，得到较好滤波效果。Goldstein方法去噪能力较强，但会损失大量边缘细节信息，给DEM和形变测量带来误差。近年发展的矢量分离式小波滤波^[8]和小波-维纳组合滤波^[9]能在一定程度上保留细节，其小波基的构造及小波软硬阈值的选取还需要良好的智能决策。卡尔曼滤波^[10]和基于格网值法^[11]都是通过估计相位来滤波图像，不同的估计方法将得到不同的滤波结果。空间域是逐像素进行的，通过周围像素来影响中心像素值达到去噪目的，将模板(通常为3×3像素)与中心像素周围对应位置进行卷积，得到新的中心像素。均值滤波和中值滤波简单易行但未考虑干涉图噪声强弱的变化及图像的条纹特征，易引起条纹混迭。作为空间域滤波方法代表的Lee滤波^[12]通过局部解缠估计中心像素邻域内的坡度，从而根据乘性噪声模型改善中心像素。由于带噪声的局部解缠计算复杂、费时且不准确，Lee提出了修正Lee滤波，即将包含边缘的局部区域分为8个方向，利用局域梯度边缘检测方法确定边缘的方向，并判定中心像素所在方向，并用该方向的像素取值改变中心像素，此方法能更好地保护条纹边缘；尹宏杰等^[13]根据同质区域方差最小的原理确定中心像素所在方向，并根据干涉图的相干性来选择融合方向窗口的数量，增强Lee算法的稳健性，但算法极为复杂。数学形态法^[14]在Lee滤波的非同质区采用膨胀腐蚀法找到边缘，能有效增强边缘。王兴旺等^[15]提出边缘保持法，其原理也是找到中心像素所在方向，该方法受窗口大小和噪声强弱影响较大。文献[16-17]应用区域增长法并以相位变化最小为准则自适应选择滤波窗口，取得较好滤波效果，但未能考虑条纹率。李佳等^[18]结合频率域和空间域滤波的优点设计二级滤波，并取得一定成效。近年来发展的堆栈法^[19]将图像分成若干个二进制文件再通过布尔判断式分别滤波。乘性加性模型法^[20]将对角线相位噪声划为乘性，非对角线噪声划为加性，该方法有待进一步验证。廖明生等^[21]提出的复数空间自适应滤波是根据梯度确定中心像素纵横方向邻域内各像素的权系数，卷积后得到中心像素的新值，其算法简单，效果明显，迭代3~5次即可。本文作者首先介绍复数空间自适应滤波算法，通过增加45°和135°梯度计算进行斜方向条纹的探索，并根据相干系数确定条纹边缘的尺度函数k，进而改善权系数的表达式；然后，以模拟数据和真实数据作为实验数据，分别与复数空间自适应、Lee滤波、基于信噪比的Goldstein及最优化方向融合方法进行比较，结果表明本文方法在去噪和边缘信息保持方面较前4种方法都有明显改善，局部自适应性更强。

1  复数空间自适应滤波方法

InSAR干涉图由复数构成，图像显示的数据为[-p，+p]之间的缠绕相位存在跳跃性，因此不能对其直接平滑，而应将复数的实部和虚部分别滤波，再进行复数重构，复数空间自适应滤波方法的算法可用下式表示。

(1)

其中：为处像素的第t次迭代，为处像素的权系数。其思想是：(x，y)处像素值的滤波值是以其周围3×3像素为窗口，对窗口内每个像素进行加权平均而得。权系数可定义为：

         (2)

其中：为信号s^(t)(x)的一阶导数，在此为梯度计算，见式(3)和(4)；k为尺度函数，确定在平滑过程中可保留的边缘幅度，通常取梯度最大值的1/3~1/2。如果某点的信号梯度大于k，则认为该点处于边缘上，边缘点在滤波窗口对应的权系数减小，即尽量不参与平滑，使边缘得以保留；反之，若信号梯度小于k，则对应非边缘的点，其权系数就会变大，在非边缘区域梯度低于k，窗口内各权系数接近，在较大程度上参与平滑，从而达到去除噪声的效果。

     (3)

     (4)

其中：G_x(x，y)为(x，y)像素的水平方向梯度；G_y(x，y)为(x，y)像素的垂直方向梯度。

2  改进的复数空间自适应滤波

复数空间自适应算法存在不足，如梯度计算只考虑水平和垂直方向，不够全面；尺度函数k的取值为1/3~1/2的最大梯度这一经验值。事实上，噪声强弱会直接影响到梯度的变化，条纹尺度函数的固定取值不能适应噪声的强弱变化，会降低对条纹判断的准确性，其局部适应性不强；此外，人工选取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)。因而，滤波窗口权系数相应地改为式(7)；

   (5)

图1  3像素×3像素的滤波窗口

Fig.1  3×3 pixels filter window

   (6)

     (7)

其中：G₄₅(x，y)为45°方向梯度；G₁₃₅(x，y)为135°方向梯度。

(2) 用相干系数改善尺度函数k。如前所述，k取值大则平滑程度大，即强去噪；取值小则平滑程度小，即弱去噪。而INSAR干涉图中噪声分布并不均匀，有的区域噪声强，有的区域噪声弱，对其滤波较理想的结果是强噪区得到强滤波，弱噪区得到弱滤波并能保持条纹形状，显然k取固定值不能满足这一要求。噪声的强度由噪声相位标准偏差决定，噪声相位标准偏差大，则噪声强；噪声相位标准偏差小，则噪声弱，噪声相位标准偏差与相干系数和视数存在以下关系^[22]：

    (8)

其中：为相位方差；γ为干涉图相干系数；L为视数；f为干涉图相位；f₀为相位的期望；f_PDF为相位概率密度函数。

 (9)

其中：；Γ为Gamma函数。图2所示为根据数值分析的方法计算得到的视数L=1，2，5，10，20以及相干性从0~1变化时的相位标准偏差理论值。从图2可见：当视数一定时，相位标准偏差与相干值成反比，即相干系数大时则噪声相位标准偏差小，此时噪声弱；相干系数小时则噪声相位标准偏差大，此时噪声强。因此，根据相干系数可反映出噪声的强弱。考虑到噪声相位标准偏差理论值难以直接求取，本文提出用下式来改善k的取值。

             (10)

其中：k_max为k在滑动窗口中的最大梯度；为计算梯度的3×3像素窗口中对应相干系数的均值。

图2  噪声相位标准偏差与相干系数及视数关系

Fig.2  Relationship between phase standard deviation and coherence&look number

α与视数L有关，通过大量数据模拟试验得到：

         (11)

由式(10)可知，当大即噪声弱时，k变小，滤波强度小；当小即噪声强时，k变大，滤波强度大；与取固定值k相比，具有更强的局部自适应性，能取得更好的滤波效果。

3  实验结果及分析

3.1  模拟干涉图实验结果

利用多分形技术模拟数字高程模型，再根据ERS1/2的成像几何参数和200 m垂直基线模拟出缠绕的相位，最后根据式(9)模拟出给定视数的相位噪声。图3(a)所示为512像素×512像素的模拟无噪相位图，图3(b)所示为加噪相位图，图3(c)~(g)所示为分别用复数空间自适应滤波(取尺度函数为1/3k_max)、Lee滤波、基于信噪比的Goldstein、最优方向融合滤波以及本文方法对其进行滤波的结果图。从图3可见：5种方法滤波效果都不错，基本复原无噪相位图。本文方法结果(图3(g))与图3(c)相比有明显改善，噪声大量减少，条纹更清晰。整体效果与图3(d)及图3(f)中的结果相似，但条纹边缘细微处信息保持得更好，与图3(a)更为接近，且运行时间要比Lee滤波及最优方向融合滤波快很多。比较图3(e)与图3(a)可见：基于信噪比的Goldstein滤波去噪能力较强，特别是同质区内的显著斑点噪声，但同时将一些细节真实信息也当成噪声被滤掉；此外，由于该算法本身的原因还需处理图边像素。无论从噪声的消除程度还是边缘信息的保护程度都可看出本文方法具有明显的优势。图4所示为与图3中对应的第280行剖面图。从该剖面图也可见：本文方法的剖面线光滑程度更高，噪声的抑制效果非常明显，更好地恢复原始无噪相位，与无噪相位图的剖面图4(a)更为吻合。

3.2  真实干涉图实验结果

本文分别选取香港地区和意大利Etna火山局部区域的干涉图作为实验数据，香港地区数据为1幅1996-03-18—19 ERS—1/2卫星香港地区的InSAR干涉图，其垂直基线长为100 m。该图条纹密度适中，在此取700像素×1 000像素矩形区域作为实验区域，其滤波结果如图5所示。由于真实图噪声的复杂性，上述几种滤波方法的效果与模拟图相比稍弱，但都明显地去除噪声，增强条纹。为便于观察，选定横向800~925像素、纵向340~425像素区域局部放大，如图6所示。图6(a)所示为分布大量噪声的原始干涉图，图6(b)~(f)所示分别为仍是复数空间自适应滤波、Lee滤波、基于信噪比的Goldstein、最优方向融合滤波以及本文方法的滤波结果。从图6可见：本文方法和复数空间自适应滤波基本去除噪声，而其他3种方法还留有残余噪声。而图6(f)与图6(b)相比，既保持条纹的曲度又适当平滑条纹边缘，效果更好。运行速度方面，Lee滤波所需时间8 min，最优方向融合需要4 min，其他方法只需2 min左右。从第240行的相位剖面图7(a)~(f)可进一步看出：图7(c)~(e)中还有明显较多的噪声，黑框中的图7(b)还残存毛刺；图7(c)~(e)中的相位受噪声影响存在明显的迂回停滞，本文方法的剖面图7(f)几乎没有毛刺，噪声也大量减少，整条曲线比较平滑，不仅去噪效果良好，而且有效保持相位细节。

意大利Etna火山100×200像素区域的密集条纹干涉图及滤波结果如图8所示。空间自适应滤波的条纹稍显粗糙，而Lee滤波部分条纹明显断开，条纹“僵直”；最优方向融合法条纹很不明显，且光滑度不好；基于信噪比的Goldstein方法的条纹光滑度及曲度保持也不如本文方法明显。

仍以模拟图和香港干涉图为例，用残余点数、相位梯度和(简称为SPD)对上述5种滤波方法进行定量评价，结果见表1。模拟图的各种滤波方法的指标均较理想，而复数空间自适应滤波由于模拟图的噪声较多而受影响，残余点改善率稍低，说明k_max/3的尺度函数未能适应模拟图的噪声水平。从表1可知：虽然香港干涉图噪声不强，但由于真实噪声分布不均匀，随机性大，因而各项指标均有所下降。此处复数空间自适应滤波下降不大，说明该尺度函数更适应于香港干涉图，从侧面说明尺度函数应该随噪声水平而变。无论模拟图还是真实图，本文方法与复数空间自适应方法相比都有较大改善，且优于其他3种方法，模拟图残余点改善率由79.6%提高到98.2%，香港干涉图改善率由74.6%提高到91.8%，分别提高18.6%和17.2%，说明本文方法更具有优势。

图3  模拟干涉图滤波结果

Fig.3  Simulated interferograms and results of different filters

图4  模拟干涉图滤波结果第280行剖面图

Fig.4  Profiles of line 280 of filtering results of simulated interferogram

图5  香港地区干涉图滤波结果

Fig.5  Filtering results of interferogram in Hong Kong region

图6  香港滤波结果局部放大图

Fig.6  Enlarged views of interferogram in local Hong Kong region

图7  香港干涉图滤波结果第240行剖面图

Fig.7  Profiles of line 240 of filtering results in Hong Kong region

图8  意大利Etna火山截取干涉图滤波结果

Fig.8  Filtering results of interferogram over part of Etna volcano in Italy

表1  干涉图滤波结果定量比较

Table 1  Quantative comparison among filtering results of simulated interferogram

4  结论

(1) 复数空间自适应滤波的目的是找到与原始像素值相近的邻域像素尽可能参与平滑，因此，梯度小则认为该点与原始像素相近，相应地将其权重加大，与原始像素一起平滑；而梯度大则认为该点与原始像素不在同一条边缘上，将其权重减小，尽可能不参与平滑。实质上是给原始像素找到同质或相近点，并对这些点根据接近程度进行加权平均。

(2) 由于条纹走向有可能是斜方向，因此，考虑45°和135°方向的梯度，可以更准确地判断原始中心像素与周围像素的关系，有助于找到遗漏的同质点，从而更合理地确定权系数。

(3) 噪声强弱决定滤波强弱，噪声可由干涉图相干系数表征，相干系数小时则噪声强，此时较大的尺度函数才能检测出条纹；相反，相干系数大时则噪声强，此时较小的尺度函数就能检测出条纹。因而加入相干系数的干预，使得复数空间自适应滤波的稳健性更强，局部自适应能力更强。

参考文献：

[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 RemoteSensing, 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. Anew technique for noise filtering of SAR interferogram phase images[J]. IEEE Transaction Geoscience and Remote Sensing, 1998, 36(5): 1456-1465.

[13] 尹宏杰, 李志伟, 丁晓利, 等. InSAR干涉图最优化方向融合滤波[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干涉条纹图的复数空间自适应滤波[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.

(编辑  邓履翔)

收稿日期：2012-08-01；修回日期：2012-09-29

基金项目：国家自然科学基金资助项目(40974006，40901172)；教育部新世纪优秀人才支持计划(NCET-08-0570)；中南大学学位论文创新基金资助项目(2010ssxt228)

通信作者：易辉伟(1974-)，女，江西萍乡人，博士研究生，从事INSAR数据处理及变形监测研究；电话：18684757287；E-mail：yhw74@163.com