基于异向梯度加权值干涉图滤波

易辉伟^{1, 2}，朱建军^{1, 2}，许兵^{1, 2}，陈建群^{1, 2}，李佳^{1, 2}

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

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

摘要：根据窗口内周围各像素与中心像素之间的异向梯度及像素相干值构成中心像素迭代公式中的权函数，提出一种保持条纹细节的干涉条纹图滤波方法。以模拟数据和意大利Etna火山、土耳其地震的真实干涉图数据为例，将新方法与空间自适应滤波、Lee滤波和最优化方向融合方法进行比较。研究结果表明：新算法能够滤除大量噪声同时保持条纹边缘，模拟干涉图和2幅真实干涉图残余点去除率分别达到99.3%，90.3%和87.9%，更有利于后续的相位解缠。

关键词：InSAR；干涉图滤波；异向梯度；相干值

中图分类号：237.3          文献标志码：A         文章编号：1672-7207(2014)08-2708-08

A filter based on weighting by anisotropic gradient for interferogram

YI Huiwei^{1, 2}, ZHU Jianjun^{1, 2}, XU Bing^{1, 2}, CHEN Jianqun^{1, 2}, LI Jia^{1, 2}

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

2. Key Laboratory of Precise Engineering Surveying & Deformation Hazard Monitoring of Hunan Province,

Changsha 410083, China)

Abstract: The kernel pixel’s iteration formula was weighted by coherence and anisotropic gradients between surrounding pixels of the window and the kernel pixel. By this means a new filter for interferogram was presented, which could keep the fringe detail as well. Taking the simulated data and actual interferograms of Etna volcano in Italy and earthquake in Turkey for example, this new algorithm was compared with complex space domain adaptive filter, Lee filter and optimal integration-based adaptive direction filter. The results show that this new method performs well not only in noise reduction but also in preservation of the sharp sawtooth profile of fringe. The residual points reductions in simulate and actual interferograms are 99.3%，90.3% and 87.9% respectively, in favor of following phase unwrapping.

Key words: InSAR; interferogram filter; anisotropic gradient; coherence

合成孔径干涉雷达(synthetic aperture radar interferometry, InSAR)已成为高分辨率、高精度测定大范围地面形变的强大工具之一。目前，该技术不仅用于地震、火山等自然灾害引起的形变监测，还应用于人类活动，如采矿、油气开采及地下水抽取等造成的地面形变监测^[1]。然而，由于去相关因素及系统热噪声产生的大量相位噪声使得干涉图条纹模糊，严重制约干涉技术的精度，因此，寻求有效滤除干涉图噪声的方法成为利用InSAR技术进行高精度测量的关键问题之一。干涉图像滤波分为频率域和空间域滤波方法。经典频率域滤波方法为Goldstein方法^[2]，该方法将像素复数值进行傅里叶变换，得到功率谱，经指数运算后平滑再反傅里叶变换到图像空间。随后Baran等^[3-6]将其滤波因子改进时考虑了视数及相位标准偏差等影响因素，使其更具自适应性。Goldstein滤波及其改进方法能滤除大部分相位噪声，但会损失边缘细节，给数字高程模型(digital elevation model, DEM)和形变测量带来误差，此外在极强噪区放弃滤波。空间域滤波方法的代表为Lee滤波方法^[7]，其通过对局部干涉图解缠，并根据中心像素邻域内的坡度将中心像素重新赋值，但解缠困难且不准确。在此基础上复数空间自适应^[8]和最优方向融合方法^[9]有所发展，前者根据梯度确定中心像素纵横方向邻域内各像素的权系数，卷积后得到中心像素的新值，但未考虑其他邻域像素对中心像素的影响；后者根据同质区域方差最小的原理确定中心像素所在方向，并选择方向融合平滑，增强了Lee算法的稳健性，计算相对复杂。数学形态法^[10]在Lee滤波方法的非同质区进行膨胀腐蚀得到边缘，有效增强了边缘，但强噪区效果有待提高。近年来，国外对空间域滤波方法研究较多，乘性加性模型法^[11]将对角线相位噪声划为乘性，非对角线噪声划为加性，该假设有待验证；堆栈法^[12]将图像分成若干个二进制文件再通过布尔判断式分别滤波，其布尔判断式还需自适应性；文献[13-14]应用区域增长法并以相位变化最小为准则自适应选择滤波窗口，取得较好滤波效果，但未能考虑条纹率。由于空域滤波边缘保持能力强，在区分水体和陆地^[15]、大气水汽研究乃至目前成为热点的小基线形变监测中得到广泛应用^[16]。组合滤波结合2种以上滤波方法的优点，如小波-维纳组合滤波^[17]、等值线法-Goldstein二级滤波^[18]等方法，前者小波阈值的选择还需智能化决策，后者强噪声区的等值线提取有待进一步研究。复数空间自适应方法在中心像素的水平和垂直方向计算各像素的梯度，采用以e为底的指数函数确定平滑时的权，未考虑周围像素和相位值的可靠性；最优方向融合方法将中心像素列为某方向上一员，然后在该方向上平滑，理论上是可行的，然而由于曲线多样性使得方向的准确界定非常困难。异向梯度加权中心值法毋需确定方向，由周围像素相对于中心像素的异向梯度确定两者之间的接近程度，并根据周围像素相干值确定其可靠程度，综合考虑这2个因素建立权函数式，使得与中心像素相似且相位值可靠的周围像素所赋予的权较大，尽可能参与平滑。本文作者以模拟数据和真实数据作为实验数据，分别与复数空间自适应、Lee滤波、最优化方向融合方法进行比较。

1  异向梯度加权值法

根据相干系数与相位值可靠程度的关系，噪声相位标准偏差与相干系数及视数的关系表达式^[19]为：

    (1)

式中：为相位方差；为干涉图相干值；L为视数；为干涉图相位；为相位的期望值；F为相位概率密度函数，其表达式为：

   (2)

式中：；为Gamma函数。为直观起见，图1所示为视数L=1，2，5，10，20，相干值0~1时的相位标准偏差理论值。可见：相干值增大则噪声相位标准偏差减小，相位受噪声污染的程度低，其值更为可靠；反之，其值可靠程度降低。

由上述可知，利用窗口像素滤波应考虑像素的可靠程度。作为空间域滤波方法的一种，异向梯度加权法不失一般性，分别在InSAR干涉复数分解成的实部图和虚部图逐像素进行，算法如下：

 (3)

式中：x和y分别为实部图像素所在的行和列；R^(k+1)(x，y)为实部图(x，y)像素的第k+1次迭代，R⁽⁰⁾为(x，y)像素的初始值；为像素的权系数，m，n为以(x，y)像素为中心的窗口大小，取值为大于1的奇数。式(3)的实质是对以(x，y)像素为中心的窗口像素进行加权平均，将所得新值赋予中心像素(下文称核像素)。加权的策略是使和核像素处于同一条纹的像素得到较大的权值，使其尽量参与平滑；处于边缘点的像素得到小的权，尽量不参与平滑。权系数可定义为：

       (4)

式中：为像素相对于核像素的梯度，因其代表不同方向，在此称为异向梯度；T为条纹梯度临界值，经验取值为梯度最大值的1/2~4/5；为像素对应的相干值；和分别为相干图的最大值和最小值。异向梯度加权值法的策略为：根据异向梯度判断窗口像素与核像素是否位于同一条纹，若某窗口像素的信号梯度大于T，则该像素不与核像素位于同一条纹，其权系数减小，尽量不影响核像素的值；反之，该像素与核像素位于同一条纹，其权系数增大，影响核像素的平滑值。此外，窗口像素值的可靠程度与其对应的相干值有关：相干值较大，则该像素相位值可靠程度较高，应尽量参与平滑。本文采用最大最小值法对相干值进行归一化以表示相对可靠程度，决定其参与平滑的权值。当梯度为0且相干值最大时，其权为最大值1；相干值最小时权为0。新方法的具体步骤如下：

(1) 选择一个合适大小(如5×5)的窗口，求取核像素所在行像素的平均值作为核像素初始值R₀。

(2) 计算窗口各像素与核像素之间距离，

               (5)

           (6)

式中：m为表示窗口大小的奇数；l为核像素位置值；i，j为窗口像素位置值；D(i，j) 为窗口像素(i，j)与核像素之间的距离。此处核像素(l，l)与自身的距离设为1。

(3) 计算各像素相对于核像素的梯度，

             (7)

式中：为像素(i，j)相对于核像素的梯度；R(i，j)为像素(i，j)的实部值。

(4) 按式(4)计算各像素的权；按式(3)计算各像素加权平均值，赋予核像素。

(5) 虚部图滤波步骤同上，将滤波后的实部和虚部重新组合成复数，计算幅角即可得到滤波后的相位条纹图。

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

Fig. 1  Relationship between phase standard deviation and coherence and look number

2  实验结果及分析

2.1  模拟干涉图验证

根据分形理论模拟多分形数字高程模型(大小为512像素×512像素，分辨率30 m)，然后以ERS1/2的成像几何参数，并假定垂直基线为200 m，模拟“真实”干涉图，如图2(a)所示。并将根据式(2)模拟得到视数为2的相位噪声添加到真实干涉图，得到含噪声干涉图，如图2(b)所示。在模拟相干性的过程中，考虑了时间去相关、空间去相关等因素的影响。

将图2(b)分解为实部图和虚部图滤波，参数选取因考虑到条纹稀疏而取7×7窗口，T取最大梯度的3/5，按前述步骤滤波，迭代2次，将滤波结果组合成复数并求出相位。图2(c)~2(f)所示分别为复数空间自适应滤波、Lee滤波、最优方向融合滤波以及本文方法滤波的结果。比较可知：自适应滤波结果(图2(c))也经过2次迭代，但仍残存较多噪声，滤波强度不如本文方法，主要因为只考虑纵横向像素梯度，且只取3×3的小窗口无法顾及到条纹宽度，导致较大的独立斑点被当作信号保留，如右下部有明显斑点。Lee滤波结果图大部分区域光滑，但左下角仍残留噪声，右中部条纹内部细节反被滤掉，造成偏差。由于相位值含噪声，使得局部解缠困难甚至错误。方向融合滤波结果图和原图基本接近，但右中下部仍潴留若干点状区域，且条纹边缘像素较为弥散，此为边缘曲率较大时用较大模板判断边缘点所在方向存在不确定性所致；本方法能够在邻域中找到处于同一边缘的点进行降噪，同时排斥了其他条纹的点，既能达到条纹内部光滑，又能使条纹外边缘鲜明，结果保持了原来条纹边缘的收敛、连贯，无论从去噪程度还是条纹细节保留方面都与原图最接近。

图2  模拟干涉图滤波结果

Fig. 2  Results of different filters for simulated interferogram

图3所示为与图2对应的各滤波结果图中第54行剖面图，原始无噪相位剖面线光滑，在±p之间缠绕，加噪后曲线相位十分杂乱，将原有曲线完全遮蔽，由此可见噪声模型较真实地模拟了随机噪声。由图3(c)~3(f)可以看出，复数空间自适应滤波、Lee滤波、最优方向融合滤波相位曲线比较粗糙，甚至出现毛刺，表明其滤波结果仍不同程度地残留了噪声，而本文方法的剖面线光滑程度高，沿线的相位起伏与原始无噪图更吻合，其去噪能力和保真度优于前3种方法。

2.2  真实干涉图验证

分别用意大利Etna火山和土耳其地震差分干涉图为实验数据，Etna火山SAR数据获取时间分别为2000-09-06和2009-10-11，间隔35 d，垂直基线为305 m；土耳其地震SAR数据获取时间分别为2011-10-10和2011-10-26，间隔16 d，垂直基线为673 m。

图4(a)所示为Etna火山的原始干涉局部图，图4(b)所示为其相干系数图，大小为300像素×300像素，由相干图可见该区域干涉图相干性较好，相干值均值达到0.57。

图5所示为复数空间自适应滤波、Lee滤波、最优方向融合滤波方法结果。由于真实图噪声的复杂性，滤波效果相比模拟图而言稍弱，但去噪明显，条纹增强了。复数空间自适应滤波方法的边缘细节较丰富，但边缘处仍残存少量噪声，条纹完整性和连贯性不够突出；Lee滤波方法的噪声分布不均匀，上部较为“干净”，右下部却有较多颗粒状噪声，曲度与原图有较大差异，说明其解缠结果受噪声影响较大；最优化融合滤波方法结果的条纹明显增强，但局部区域出现条纹粘连，高噪区条纹断裂现象严重；本文方法结果在低噪区条纹清晰，曲度适中，细节保持，在高噪区(左下部)能恢复部分被噪声淹没的条纹，保持了条纹的连续性和完整性，其局部自适应性更强。

图3  模拟干涉图滤波结果第54行剖面图

Fig. 3  Profile of line 54 of filtering results of simulated interferogram

图4  Etna火山原始干涉图及相干图

Fig. 4  Original interferogram and coherence map over Etna volcano

图5  Etna地区干涉图滤波结果

Fig. 5  Results of different filters over Etna region

图6所示为土耳其地震的原始干涉局部图和相干系数图，均为300像素×300像素，相干性较低，相干值均值仅为0.36。土耳其地震干涉图滤波结果如图7所示。从图7可知：复数空间自适应滤波和Lee滤波结果残余噪声明显较多；最优方向融合滤波结果条纹受强噪声的影响而不够连贯；本文方法几乎完整地复原所有条纹，条纹完整连续，再次验证了该方法的优越性。

图6  土耳其地震原始干涉图及相干图

Fig. 6  Original earthquake interferogram and coherence map over Turkey region

图7  土耳其地震干涉图滤波结果

Fig. 7  Results of different filters over Turkey region

本方法的窗口大小应根据条纹密集程度而定。实验证明，当条纹稀疏时，如土耳其地震干涉条纹图，取较大窗口(如9×9)不仅提高了滤波效率，还可以过滤较大的独立斑点噪声，运行该程序时一般迭代2次即可；当条纹较密集时，如意大利Etna火山干涉条纹图,较大窗口容易引起条纹混叠，因此需取小窗口(3×3)多次迭代；条纹尺度T的取值与噪声程度有关，噪声弱则T取值小，窗口像素权值小，参与平滑程度低，因而得到弱去噪；噪声强则T取值大，窗口像素权值大，参与平滑程度高，因而得到强去噪。窗口大小和系数梯度临界值的选取有待进一步研究。

2.3  定量指标评价

模拟和真实干涉图的滤波结果采用常用定量评价指标有残差点数、相位标准偏差(PSD)及相位梯度和值(SPD)等，结果如表1和表2所示。从表1可知：模拟图加噪后3项指标均为最高，经滤波后均大幅度降低，程度由高到低依次为本文方法、最优化方向融合方法、Lee滤波及复数空间滤波方法，其中本文方法的残差点去除率达到99.3%，PSD和SPD略高于最优化方向融合方法，改善率分别达到73.4%和59.0%。

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

Table 1  Quantative comparison among filtering results of simulated interferograms

从表2可知：意大利Etna火山和土耳其地震真实干涉图的残差点去除率分别为90.3%和87.9%，稍低于模拟干涉图的结果。原因是真实干涉图的噪声比模拟图的更为复杂，此外，模拟图根据几何参数模拟的噪声分布更均匀，相干值均值较高，条纹更稀疏，因而模拟图滤波效果更好。

采用本文方法滤波，意大利Etna火山图的PSD和SPD改善率分别为25.5%和53.7%，土耳其地震干涉图PSD和SPD改善率分别为17.4%和44.3%，略低于前者。主要是因为两者的噪声程度和条纹密集程度不同而导致滤波结果的差异。经本文方法滤波后的干涉图的3项指标明显优于前3种方法。

表2  真实干涉图滤波结果定量比较

ctualTable 2  Quantative comparison among filtering results of actual interferograms

3  结论

(1) 基于异向梯度加权值滤波方法避开寻找核像素所在方向这一复杂的问题，直接以异向梯度判断它与周围各像素的亲疏关系并定权，且根据各像素对应相干值调节权大小以增强权函数式的稳健性，效果良好。

(2) 模拟和真实干涉图滤波实验验证结果表明，无论从目视效果还是定量指标来看，本文方法在去噪和条纹保持方面都有较大优越性，甚至能检索出部分隐蔽条纹，局部适应性强。

(3) 像素窗口大小及梯度临界值的自适应选取有待进一步研究。

参考文献：

[1] Samsonov S. Topographic correction for ALOS PALSAR interferometry[J]. IEEE Geoscience and Remote Sensing, 2010, 48(7): 3020-3027.

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

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

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

[5] 孙倩, 朱建军, 李志伟, 等. 基于信噪比的干涉图自适应滤波[J]. 测绘学报, 2009, 38(5): 437-442.

SUN Qian, ZHU Jianjun, LI Zhiwei, et al. A new adaptive interferogram filter based on SNR[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(5): 437-442.

[6] 于晓歆, 杨红磊, 彭军还. 一种改进的 Goldstein 干涉图滤波算法[J]. 武汉大学学报(信息科学版), 2011, 36(9): 1051-1054.

YU Xiaoxin, YANG Honglei, PENG Junhuan. A modified Goldstein algorithm for interferogram filtering[J]. Geometrics and Information Science of Wuhan University, 2011, 36(9): 1051-1054.

[7] LeeJS,Papathanassiou KP, Ainsworth TL, et al. Anew technique for noise filtering of SAR interferogram phase images[J]. IEEE Transaction Geoscience and RemoteSensing, 1998, 36(5): 1456-1465.

[8] 廖明生, 林珲, 张祖勋, 等. 干涉条纹图的复数空间自适应滤波[J]. 遥感学报, 2003, 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.

[9] 尹宏杰, 李志伟, 丁晓利, 等. 干涉图最优化方向融合滤波[J]. 遥感学报, 2009, 13(6): 1099-1105.

YIN Hongjie, LI Zhiwei, DING Xiaoli, et al. Optimal integration-based adaptive direction filter for interferogram[J]. Journal of Remote Sensing, 2009, 13(6): 1099-1105.

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

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

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

[13] 郭交, 李真芳, 刘艳阳, 等.一种干涉相位图的自适应滤波算法[J]. 西安电子科技大学学报, 2011, 38(4): 77-88.

GUO Jiao, LI Zhenfang, LIU Yanyang, et al. New adaptive noise suppressing method for interferometric phase images[J]. Journal of Xidian University, 2011, 38(4): 77-88.

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

[15] Fugura A A, Billa L, Estuarine B P. Semi-automated procedures for shoreline extraction using single RADARSAT-1 SAR image[J]. Estuarine, Coastal and Shelf Science, 2011, 95(4): 395-400.

[16] Goel K, Adam N. An advanced algorithm for deformation estimation in non-urban areas[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2012, 73: 100-110.

[17] 蔡国林, 李永树, 刘国祥. 小波-维纳组合滤波算法及其在干涉图去噪中的应用[J]. 遥感学报, 2009, 13(1): 129-136.

CAI Guolin, LI Yongshu, LIU Guoxiang. Wavelet-wiener combined filter and its application on interferogram[J]. Journal of Remote Sensing, 2009, 13(1): 129-136.

[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] Franceschetti G, Lanari R. Synthetic aperture radar processing[M]. Florida: CRC Press, 1999: 185-190.

(编辑  赵俊)

收稿日期：2013-08-07；修回日期：2013-11-05

基金项目：国家高技术研究发展计划(“863”计划)项目(2012AA121301)；国家自然科学基金资助项目(41222027)；教育部博士点基金资助项目(20130162110015)；湖南省杰出青年基金资助项目(13JJ1006)；湖南省教育厅平台开放基金资助项目(12K010)

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