干涉图像第二类统计Goldstein自适应滤波方法

赵文胜，蒋 弥，何秀凤

河海大学地球科学与工程学院，江苏 南京 211100



干涉图像第二类统计Goldstein自适应滤波方法

赵文胜，蒋 弥，何秀凤

河海大学地球科学与工程学院，江苏 南京 211100

干涉图滤波是InSAR数据处理中的关键步骤之一，滤波结果的优劣会直接影响到相位观测的质量和最终产品精度。本文结合干涉图滤波算法的研究进展，对Goldstein频率域滤波及其经典改进算法进行了系统分析和比较，在此基础上提出了一种基于第二类统计的稳健相干性估计量的Goldstein自适应滤波方法。本文采用模拟数据和Envisat ASAR真实数据与现有方法进行了验证，试验结果表明，新的滤波方法在保持细节和抑制噪声方面优势更加明显。

Goldstein滤波；相位噪声；第二类统计；相干系数

目前，合成孔径雷达干涉测量技术(interferometric synthetic aperture radar,InSAR)在获取数字高程模型(digital elevation model，DEM)、监测大面积地表形变和极地冰川运动等方面极具潜力，但仍有许多因素影响其实际应用的精度和范围。受时间、几何、多普勒中心频率和噪声等失相关因素的制约，InSAR干涉图存在随机噪声，这不同程度地影响了干涉相位质量，在噪声水平严重的情况下，易导致后续相位解缠误差增大，甚至无法进行解缠。基于此，对干涉图进行滤波是保障InSAR产品质量的重要环节[1]。

干涉图滤波一般分为两类：第一类是空域滤波，如圆周期滤波[2]、Lee自适应滤波[3]、堆栈滤波[4]和旋滤波[5]等；第二类为频域滤波，如小波滤波[6-7]、Goldstein滤波及其改进方法[8-9]等。空域滤波继承了传统图像滤波的思想，并结合SAR图像空间特征作了改进。大量试验表明，当SAR图像地表特征比较复杂时，利用该类方法进行干涉图滤波，易造成干涉条纹混叠等失真现象，且常伴随一定的计算复杂度和时间复杂度问题[10]。相比之下，频率域滤波方法中的Goldstein滤波在去噪和保持影像分辨率方面效果更加显著，且计算较为简单，因而广泛用于商业及免费SAR数据处理软件中，已成为当前国际上的主流方法。

随着近年来InSAR数据处理方法的不断发展，在经典Goldstein滤波方法的基础上涌现出许多改进方法[9,11]。为克服Goldstein滤波参数选取的主观性，Baran滤波[12]利用干涉图相干系数以自适应地计算滤波参数，使滤波效果得到了明显改善。文献[13]借鉴Baran滤波方法，利用伪相关图自动计算滤波参数，提出了一种干涉图迭代滤波方法。该方法充分利用了干涉图自身携带的信息来指导滤波参数的计算，取得了较好的滤波结果。然而，这些自适应方法也存在一些不足之处。如Baran滤波[12]在滤波设计时并未考虑相干系数估计量的有偏性和低精度，造成滤波参数的估值偏低。另外，文献[13]的迭代Goldstein滤波算法中估计量的统计性质尚不明确，且伪相关图的质量易受SAR图像中的地形信息的影响[9,14]，这导致很难评估其滤波参数的优劣。从模拟结果来看，伪相干图的值要明显大于相干系数真值，故滤波参数也是被低估的。为解决以上问题，本文引入第二类统计方法对相干系数进行修正，其主要依据是该统计量能够兼顾系统有偏性和方差，从而精确估计滤波参数以便获得最优的滤波干涉图。

1 基于第二类统计相干性估计量的自适应滤波

本节在回顾经典Goldstein滤波方法的基础上，引入第二类统计方法对滤波参数进行稳健估计，并结合待滤波干涉图的干涉视数获得干涉图滤波模型。

1.1 Goldstein滤波

Goldstein频域滤波算法[15]的原理是将干涉图分割为相互重叠的滑动窗口(干涉图小块)，通过二维傅里叶变换计算每个滑动窗口相位的功率谱，并对该功率谱进行平滑

(1)

1.2 基于第二类统计的无偏相干性估计

由于SAR传感器不能在同一时刻连续成像，在估计相干性时，通常用SAR信号空间平均来代替汇集平均，求得的相干系数为

(2)

图1 统计模型相干性统计偏差对比图Fig.1 Coherence bias of different estimators

第二类统计量仅对定义在正半轴上的函数有效，根据定义，概率密度函数p(x)(x∈[0,+∞))的第二类第一特征函数φ(s)和第二类第二特征函数ψ(s)可分别表示为[20]

(3)

ψ(s)=lnφ(s)

(4)

第二类统计的第二类矩定义为

(5)

p(dD,N)=2(N-1)(1-D2)N×

(6)

可得其一阶log矩为

(7)

式中，γ代表真实相干性；2F1表示超几何函数；N表示样本数。通过计算log矩的指数，可以获得第二类统计相干性量级的等效一阶矩

(8)

(9)

(10)

图1展示了蒙特卡罗(Monte Carlo)随机模拟试验的结果。可以看出，在样本N的取值相同时，无论真实相干量级是多少，基于第二类统计的估计量较传统估计量偏差更小。当N=5时，在低相干性区域的偏差大于0.1，暗示了滤波参数差异大于0.1。

1.3 改进的自适应Goldstein滤波

在重新定义估计量之后，可以用估计量式(10)代替Baran滤波模型，获得改进的滤波参数。然而在理论上，相位噪声是相干性和干涉视数N的函数，仅考虑相干性的滤波因子不能完全反映当前干涉图噪声的情况[22]。相位标准方差σφ与干涉图相干性γ和干涉视数N之间的关系为[23]

(11)

式中，φ为干涉图相位；φ0为相位的期望值；PDF为相位概率密度函数

(12)

图2为噪声相位标准方差随干涉图相干性和干涉视数的变化情况。从图中可知，干涉图噪声水平分别随其相干性和干涉视数的增加而减小，故可利用干涉图的相干值和干涉视数重新定义新的滤波参数，从而使其能够更好地适应干涉图噪声的空变特征。

图2 噪声相位标准偏差与相干值及样本数关系Fig.2 Phase noise standard deviation versus the coherence and samples

在实际应用中，由于有偏估计和超几何函数的存在，采用方程式(2)、式(11)和式(12)联合估计相位方差会显著增加估计误差和计算复杂度。根据文献[9]，相位噪声估计的简易形式以及它与最优滤波参数的函数关系可以近似表述为

(13)

(14)

联合方程式(10)、式(13)和式(14)，得到稳健的滤波参数估计量α

(15)

综上，新方法的具体步骤如下：

(1) 利用式(2)计算干涉图的原始相干图，再利用式(10)对原始干涉图进行偏差纠正，得到近似无偏的相干系数图。

(2) 把原始干涉图分成若干小块(如32×32像素)，为保持小块之间条纹的连续性，每块之间有一定的重叠像素。

(3) 利用式(15)和步骤(1)中的相干系数估算有效区域内(小块像素减去重叠像素)的滤波参数α，再利用式(1)对干涉图进行分块滤波，最后通过逆变换得到滤波后的干涉图。

2 滤波试验对比分析

本文采用模拟数据和真实数据对Goldstein滤波、Baran滤波、迭代Goldstein滤波以及本文提出的新滤波方法进行了对比分析。为保证试验公平性，在以下试验中滤波窗口均采用典型滤波窗口，即大小均为32×32，重叠区域为28。传统Goldstein滤波中参数α=0.5，文献[13]的迭代Goldstein滤波为一次迭代滤波的结果。

2.1 模拟数据

首先，采用多分形技术和ERS卫星参数(垂直基线30 m)模拟500×500像素的无噪声相位图(见图3(b))。接着，在时空去相关、体散射去相关的影响下，模拟无偏相干系数图(见图4(a))。然后根据相位标准偏差公式和模拟的相干系数矩阵获得各像素的相位噪声(N=9)，最后把模拟的噪声逐点加到无噪声的相位图上进行缠绕，得到含噪干涉图(见图3(a))。本文将相干性的理论偏差分别加入无偏相干图后得到有偏的相干图(见图4(b)、图4(d)),伪相关系数(见图4(c))采用样本估计公式(参考文献[13]中的式(7))直接获得。

将图4(b)的相干值代入Baran滤波公式，获得Baran的滤波结果；将图4(c)的相干图代入文献[13]的迭代Goldstein滤波公式，得到迭代Goldstein滤波结果；最后，将图4(d)的相干值代入式(15)获得新方法的滤波结果。从图5可知，所有滤波方法都能得到条纹相对清晰的干涉图，表明能够去除强噪声。但定性地看，新提出的滤波方法滤波效果明显优于其他3种滤波方法，最大限度地还原了原始相位信息。

图3 含噪声和无噪声干涉图Fig.3 Noise added and free interferograms

为定量评价滤波结果，本文选取了相位差值和(thesumofphasedifference,SPD)[24]及相位均方误差(RMS)作为干涉图质量评价标准，滤波后干涉图的相位梯度值越小，表征干涉图越平滑，滤波效果越好；均方误差越小，代表滤波器的保真性越好。从表1中可以看出新方法较其他3种滤波具有更好的去噪能力,使模拟噪声干涉图的SPD减少了87.5%。本文还针对每幅干涉图的同一位置(图5中标记为A黑线所示)作了剖面分析。如图6所示，从滤波细节上看，相比其他滤波方法，新方法的整体滤波结果更接近原始无噪声相位值，从而说明了新方法在去噪的同时具有更好的信号保真性。

图4 真实相干图和有偏相干图Fig.4 True coherence and bias added coherence maps

滤波方法评价指标SPDRMS数量(×105)减少率/(%)原干涉图3.77931.2186Goldstein1.589058.00.5538Baran1.035672.60.3685迭代Goldstein1.645956.40.5698新方法0.491787.50.1950

2.2 真实数据

2.2.1 夏威夷火山区域

为验证上述各种算法在真实数据中干涉条纹稠密区域的滤波效果，本文首先选取了美国夏威夷火山区域的两幅ASAR数据并截取了其中200×200的子区域，其中，图7(a)为原始干涉图；图7(b)、图7(d)分别为基于两类统计模型的相干系数图；图7(c)为伪相关系数图；图8为原始干涉图黑色方框区域经过4种滤波方法滤波后的干涉图。

图8中4种滤波方法都能滤除一定的干涉图噪声，但是经过迭代Goldstein滤波的干涉图中含有较多噪声，这是因为估计量伪相干值也是有偏的，并产生高估偏差。因此，这一高估现象导致迭代Goldstein滤波的滤波参数被低估，在强噪声处的滤波效果降低。同理，基于传统统计的Baran滤波在低相干区域的滤波效果也较差。而在精确估计相干性后，本文提出的方法较好地抑制了噪声水平，得到了更加清晰的干涉条纹。本文采用国际上通用的滤波评价指标对结果进行评定，即相位奇异点(residue point number)[15]以及SPD。干涉图中的相位奇异点(残差点)越少，表明噪声程度越小，更利于相位解缠处理。如表2所示，新滤波可以获得最优的结果，使原来噪声干涉图的SPD值降低了36.7%，相位奇异点值则降低了71.7%。这说明在强噪区域，更加精确的相干性估计有益于提高相位质量。

图5 模拟数据的滤波结果对比图Fig.5 Filtered interferograms

图6 模拟数据滤波结果局部横断面Fig.6 Cross-sections of interferograms filtered with different methods

2.2.2 徐州矿区

本文进一步选取了中国徐州矿区的两幅PALSAR数据验证新方法在干涉条纹稀疏区域的滤波效果。为展现滤波细节，本次试验截取了干涉图中512×512像素的子区域，图9(a)—(d)分别给出了其中的原始干涉图、基于两类统计模型的相干系数图和伪相关系数图，该图进一步验证了基于第二类统计的相干性估计可获得较传统统计模型偏差更小的估计量。图10为原始干涉图经过4种滤波方法处理后的滤波结果图。

表2 夏威夷干涉图SPD和Residues滤波结果统计

Tab.2 Statistics of SPD and Residues of the Interferogram of Hawaii from different filters

滤波方法评价指标SPDresiduesnumber数量(×104)减少率/(%)数量(×103)减少率/(%)原干涉图7.30377.711Goldstein6.169215.54.38943.1Baran5.561423.93.11459.6迭代Goldstein6.219614.84.20945.4新方法4.514936.72.16471.7

图10同样给出了笔者期望的滤波结果，与其他3种滤波方法所获滤波结果相比，经过新方法处理的滤波结果图10(d)再一次展现了较好的平滑效果。此处，本小节仍然采用了相位奇异点和SPD两种滤波评价指标对来自上述4种滤波方法的滤波结果进行定量评定，所得结果已列于表3，

图7 夏威夷干涉图和干涉质量图Fig.7 Interferogram of Hawaii and estimated coherence maps using different estimators

图9 徐州地区干涉图和干涉质量图Fig.9 Interferogram of Xuzhou and estimated coherence maps using different estimators

图10 徐州干涉图滤波结果对比图Fig.10 Interferogram of Xuzhou filtered with different methods

本次试验新的滤波方法使原来噪声干涉图的SPD值和相位奇异点值分别降低了89.7%和99.5%。图11给出了图10中4幅干涉图黑色横线位置的剖面图，由于所选区域位于平原区，且为城市区域，地势较平坦，故滤波后的干涉图应该较平滑，而新滤波方法得到的结果比其他3种方法得到的都理想，无明显相位波动。因此，结合图10和图11的定量分析结果，不难发现新方法在条纹稀疏的噪声区域，也能得到较好的滤波结果。

表3 徐州干涉图SPD和Residues滤波结果统计

Tab.3 Statistics of SPD and Residues of the Interferogram of Xuzhou from different filters

滤波方法评价指标SPDResiduesnumber数量(×105)减少率/(%)数量(×104)减少率/(%)原干涉图4.8515-3.804-Goldstein1.501669.00.10497.3Baran0.982779.90.05198.7迭代Goldstein2.305752.50.14496.2新方法0.499589.70.02099.5

图11 徐州干涉图滤波结果局部横断面Fig.11 Cross-sections of interferograms of Xuzhou filtered with different methods

3 结 论

本文提出了一种新的Goldstein自适应滤波方法，该方法可以在第二类统计的基础之上获得更加精确的相干性估计量，进而获得更准确的滤波参数。从模拟和真实数据的测试情况来看，相比于现有的自适应方法，该方法能在维护分辨率的同时兼顾滤波效果，降低噪声能力较其他3种滤波方法提高了近一倍。在低相干区域，新方法因纠正了滤波参数的低估偏差而抑制了强噪声，进而解决了现有方法在低相干区域欠滤波的问题。

[1] 黄长军,郭际明,喻小东,等.干涉图EMD-自适应滤波去噪法[J].测绘学报,2013,42(5): 707-714.

HUANG Changjun,GUO Jiming,YU Xiaodong,et al.The Study of Interferogram Denoising Method Based on EMD and Adaptive Filter[J].Acta Geodaetica et Cartographica Sinica,2013,42(5): 707-714.

[2] EICHEL P,GHIGLIA D,JAKOWATZ C,et al.Spotlight SAR Interferometry for Terrain Elevation Mapping and Interferometric Change Detection[R].Sand: Sandia National Labs Technology,1993: 2539-2546.

[3] LEE J S,PAPATHANASSIOU K P,AINSWORTH T L,et al.A New Technique for Noise Filtering of SAR Interferometric Phase Images[J].IEEE Transactions on Geoscience and Remote Sensing,1998,36(5): 1456-1465.

[4] BUEMI M E,JACOBO J,MEJAIL M.SAR Image Processing Using Adaptive Stack Filter[J].Pattern Recognition Letters,2010,31(4): 307-314.

[5] 于晶涛,陈鹰.一种新的INSAR干涉条纹图滤波方法[J].测绘学报,2004,33(2): 121-126.DOI: 10.3321/j.issn:1001-1595.2004.02.006.

YU Jingtao,CHEN Ying.A New Filter on InSAR Interferogram Fringe[J].Acta Geodaetica et Cartographica Sinica,2004,33(2): 121-126.DOI: 10.3321/j.issn:1001-1595.2004.02.006.

[6] LOPEZ-MARTINEZ C,FABREGAS X.Modeling and Reduction of SAR Interferometric Phase Noise in the Wavelet Domain[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(12): 2553-2566.

[7] 何儒云,王耀南.一种基于小波变换的InSAR干涉图滤波方法[J].测绘学报,2006,35(2): 128-132.DOI: 10.3321/j.issn:1001-1595.2006.02.007.HE Ruyun,WANG Yaonan.InSAR Interferogram Filtering Based on Wavelet Transform[J].Acta Geodaetica et Cartographica Sinica,2006,35(2): 128-132.DOI: 10.3321/j.issn:1001-1595.2006.02.007.

[8] JIANG Mi,DING Xiaoli,LI Zhiwei,et al.The Improvement for Baran Phase Filter Derived from Unbiased InSAR Coherence[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2014,7(7): 3002-3010.

[9] JIANG Mi,DING Xiaoli,TIAN Xin,et al.A Hybrid Method for Optimization of the Adaptive Goldstein Filter[J].ISPRS Journal of Photogrammetry and Remote Sensing,2014,98: 29-43.

[10] 尹宏杰.基于InSAR 的矿区地表形变监测研究[D].长沙: 中南大学,2009.YIN Hongjie.Ground Subsidence Monitoring in Mining Area Using InSAR[D].Changsha: Central South University,2009.

[11] LI Z W,DING X L,HUANG C,et al.Improved Filtering Parameter Determination for the Goldstein Radar Interferogram Filter[J].ISPRS Journal of Photogrammetry and Remote Sensing,2008,63(6): 621-634.

[12] BARAN I,STEWART M P,KAMPES B M,et al.A Modification to the Goldstein Radar Interferogram Filter[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(9): 2114-2118.

[13] ZHAO Chaoying,ZHANG Qin,DING Xiaoli,et al.An Iterative Goldstein SAR Interferogram Filter[J].International Journal of Remote Sensing,2012,33(11): 3443-3455.

[14] GHIGLIA D C,PRITT M D.Two-dimensional Phase Unwrapping: Theory,Algorithms,and Software[M].New York: Wiley-Interscience,1998.

[15] GOLDSTEIN R M,WERNER C L.Radar Interferogram Filtering for Geophysical Applications[J].Geophysical Research Letters,1998,25(21): 4035-4038.

[16] BARAN I.Advanced Satellite Radar Interferometry for Small-scale Surface Deformation Detection[D].Perth: Curtin University of Technology,2004.

[17] 蒋弥,丁晓利,李志伟,等.基于时间序列的InSAR相干性量级估计[J].地球物理学报,2013,56(3): 799-811.JIANG Mi,DING Xiaoli,LI Zhiwei,et al.InSAR Coherence Magnitude Estimation Based on Data Stack[J].Chinese Journal of Geophysics,2013,56(3): 799-811.

[18] TOUZI R,LOPES A,BRUNIQUEL J,et al.Coherence Estimation for SAR Imagery[J].IEEE Transactions on Geoscience and Remote Sensing,1999,37(1): 135-149.

[19] ABDELFATTAH R,NICOLAS J M.Interferometric SAR Coherence Magnitude Estimation Using Second Kind Statistics[J].IEEE Transactions on Geoscience and Remote Sensing,2006,44(7): 1942-1953.

[20] NICOLAS J M,ANFINSEN S N.Introduction to Second Kind Statistics: Application of Log-Moments and Log-Cumulants to the Analysis of Radar Image Distributions[R].Technical Note,Translation from French,2012.

[21] TOUZI R,LOPES A.Statistics of the Stokes Parameters and of the Complex Coherence Parameters in One-Look and Multilook Speckle Fields[J].IEEE Transactions on Geoscience and Remote Sensing,1996,34(2): 519-531.

[22] SUO Zhiyong,ZHANG Jinqiang,LI Ming,et al.Improved InSAR Phase Noise Filter in Frequency Domain[J].IEEE Transactions on Geoscience and Remote Sensing,2016,54(2): 1185-1195.

[23] HANSSEN R F.Radar Interferometry: Data Interpretation and Error Analysis[M].Netherlands: Springer,2001.

[24] LI Zhilin,ZOU Weibao,DING Xiaoli,et al.A Quantitative Measure for the Quality of InSAR Interferograms Based on Phase Differences[J].Photogrammetric Engineering & Remote Sensing,2004,70(10): 1131-1137.

(责任编辑：陈品馨)

Improved Adaptive Goldstein Interferogram Filter Based on Second Kind Statistics

ZHAO Wensheng,JIANG Mi,HE Xiufeng

School of Earth Sciences and Engineering,Hohai University,Nanjing 211100,China

Interferometric filtering is one of the most important procedures in InSAR data processing as it refers to the quality of phase observations and the final products.A variant under the framework of Goldstein filter is presented in this paper,which is based on the robust coherence estimator from the second kind statistics.Compared with the state-of-the-art,the significant advantage of the new method is that more accurate filtering parameter alpha can be deduced and therefore better performance of the Goldstein filtering can be expected.Experimental results from both simulation and real Envisat ASAR data demonstrate the value of the method.

Goldstein filter; phase noise; second kind statistics; coherence

The National Natural Science Foundation of China (Nos.41404009; 41274017); The National Key Technology Research and Development Program(No.2015BAB07B10); The Surveying and Mapping Basic Research Program of National Administration of Surveying,Mapping and Geoinformation (No.15-01-04); The Transportation Technology and Achievements Transformation Project of Jiangsu Province(No.16Y08)

ZHAO Wensheng(1985—)，male,PhD candidate,majors in processing and application of InSAR dataset.

JIANG Mi

赵文胜，蒋弥，何秀凤.干涉图像第二类统计Goldstein自适应滤波方法[J].测绘学报，2016,45(10)：1200-1209.

10.11947/j.AGCS.2016.20150457.

ZHAO Wensheng,JIANG Mi,HE Xiufeng.Improved Adaptive Goldstein Interferogram Filter Based on Second Kind Statistics[J].Acta Geodaetica et Cartographica Sinica,2016,45(10):1200-1209.DOI:10.11947/j.AGCS.2016.20150457.

P237

A

1001-1595(2016)10-1200-10

国家自然科学基金(41404009; 41274017)；国家科技支撑计划(2015BAB07B10)；国家测绘地理信息局测绘基础研究基金(15-01-04)；江苏省交通运输科技与成果转化项目(16Y08)

2015-09-10

修回日期：2016-07-28

赵文胜(1985—)，男，博士生，研究方向为InSAR数据处理与应用。

E-mail：wensheng_zh@hhu.edu.cn

蒋弥

E-mail：mijiang@hhu.edu.cn