结合局部频率估计的小波域InSAR相位滤波新方法

2015-06-05李芳芳胡东辉丁赤飚

系统工程与电子技术 2015年12期

关键词：子带条纹小波

李芳芳，林雪，胡东辉，丁赤飚

（1.中国科学院电子学研究所，北京100190；2.中国科学院空间信息处理与应用系统技术重点实验室，北京100190）

结合局部频率估计的小波域InSAR相位滤波新方法

李芳芳1，2，林雪1，2，胡东辉1，2，丁赤飚1，2

（1.中国科学院电子学研究所，北京100190；2.中国科学院空间信息处理与应用系统技术重点实验室，北京100190）

干涉相位噪声的存在直接影响相位解缠的效果及最终干涉测量的精度。针对这一问题，本文提出了一种结合局部频率估计的小波域干涉合成孔径雷达（interferometric synthetic aperture radar，InSAR）相位滤波方法。根据局部频率估计可以判断出干涉相位小波系数中包含有用信息的子带。结合VisuShrink和NeighShrink两种小波阈值收缩方法分别具有去噪效果好和细节保持能力强的特点，对有用信息所在子带的小波系数利用NeighShrink方法进行阈值收缩，而对其他子带的小波系数则利用VisuShrink方法进行阈值收缩，从而尽可能地滤除噪声，同时保持干涉条纹的细节信息不被破坏。仿真和实测数据实验验证了本文提出的滤波方法的有效性。

干涉合成孔径雷达；相位滤波；小波变换；局部频率估计；阈值收缩

0 引言

合成孔径雷达（synthetic aperture radar，SAR）干涉测量技术通过从同一场景多次观测的SAR复图像中获取相位信息，从而反演地表的高程信息和地形变化信息［1］。干涉测量的精度和可靠性在很大程度上取决于干涉相位图的质量。然而，在实际系统中，受热噪声去相干、时间去相干、基线去相干、信号处理去相干等多种去相干因素的影响，干涉相位噪声难以避免［2］。相位噪声不仅给相位解缠带来困难，而且影响最终的干涉测量精度。因此，在相位解缠前必须对干涉相位进行滤波，从而获得较为准确的干涉相位估计值。

目前，干涉相位的滤波方法大体可以分为空间域滤波和频率域滤波两类。圆周期均值或中值滤波［3］是最常用的空间域滤波方法，它实现简单，运算速度快，但是滤波窗口难以确定，在条纹密集时容易破坏相位细节，降低分辨率。在此基础上，发展出了基于局部频率补偿的自适应滤波方法［45］和自适应选择窗口大小和方向的滤波方法［67］。近年来，发展出了基于非局部的滤波方法［89］，并有学者将其应用到了干涉相位处理中［1011］，与传统基于局部窗口的滤波方法不同，该方法基于整幅图像，利用图像的纹理相似性来滤除噪声，可以更好地保留图像细节，但该方法的计算复杂度过高，运算效率较低。频率域滤波方法中，Goldstein滤波［12］应用最为广泛，但该方法受分块大小和滤波参数的影响较大，在信噪比很低时，滤波效果较差。另外，文献［13］提出了基于信号子空间投影的干涉相位估计方法，可以在配准存在误差的情况下获得较好的估计结果，然而在相干性较低时，噪声子空间的维数难以判断，影响估计的准确性。

小波变换由于其良好的时频分析特性和多分辨率特性，在信号分析和图像处理领域具有广泛的应用。Carlos等人提出利用离散小波变换（discrete wavelet transform，DWT）进行干涉相位滤波［14］，通过增大小波系数中的信号成分来实现，该方法能够很好地保持干涉条纹的细节信息，并在一定程度上提高了图像的信噪比。但由于对噪声信息没有进行有效的抑制，使得去噪效果不够好。因此，本文提出一种结合局部频率估计和小波阈值收缩的干涉相位滤波方法，通过局部频率估计结果区分干涉相位小波变换后条纹信息及噪声所在的子带，分别利用NeighShrink和VisuShrink方法进行阈值收缩，从而在滤除噪声的同时保持干涉条纹的细节信息不被破坏。仿真和实测数据实验验证了本文算法的有效性和适应性。

本文的内容安排如下：第1节对小波域中的干涉相位模型进行简单介绍；第2节阐述了两种小波阈值收缩方法的阈值选取原则及各自的优缺点；第3节在介绍局部频率估计方法的基础上，给出了本文提出的新的滤波方法流程；第4节利用仿真和实际数据对本文提出的干涉相位滤波方法进行了验证；最后对文章进行了总结。

1 小波域干涉相位模型

在实数域中，干涉相位噪声符合加性噪声模型的特征［7］，可以表示为

式中，φz为实际的干涉相位；φx为理想的干涉相位；v为相位噪声，均值为零。受三角测量的限制，干涉相位缠绕在（－π，π］之间，为了避免相位缠绕可能导致的误差，大多数干涉相位滤波算法在复数域进行，即有

根据式（1）可知，式（2）的实部和虚部［14］可以分别表示为

由此可见，Nc值与干涉相干性有关。当相干性高时，小波系数中的有用相位信息应占主导，此时由于Nc值较大，故小波系数的幅度值也较大；当相干性差时，小波系数中的噪声占主要成分，由于Nc值较小，故小波系数的幅度值也较小。因此，通过设置阈值可以区分相位信息和噪声对应的小波系数，分别进行处理，从而实现干涉相位滤波［15］。

小波变换的多分辨率特性使得不同的子带对应不同的频率范围。图1给出了尺度为2时的二维DWT示意图，在进行第i尺度的变换后，各个子带对应的频率范围如下：

图1 尺度为2的二维DWT示意图

2 小波阈值收缩方法

小波变换具有将信号或图像进行稀疏表达的能力。也就是说，含噪声的信号或图像经过小波变换后，由有用信息产生的小波系数仅集中在较少的位置和尺度上，且通常具有较大的幅值；而由噪声产生的小波系数在时频空间分布较广，但幅值很小。小波阈值收缩算法就是基于这一特征，通过抑制低于一定阈值的小波系数，而保留高于该阈值的小波系数，从而达到去除噪声的目的。选择合适的阈值对于小波阈值收缩算法至关重要，本文将结合VisuShrink和NeighShrink两种阈值收缩算法进行干涉相位滤波，下面分别介绍其阈值选取原则。

2.1 VisuShrink方法

由Donoho和Johnstone提出的VisuShrink方法是目前最为经典的阈值收缩方法［16］，也称为通用阈值收缩方法。对于二维图像而言，该方法的实现过程为：首先通过小波变换获得噪声图像的小波系数，其中i为小波尺度，m，n为小波系数的位置。然后根据式（8）计算通用阈值T：

式中，N为信号长度或图像尺寸；σ为噪声标准差，按式（9）进行估计

由于高频子带上的小波系数主要由噪声引起，因此这里仅利用第一层小波分解后的对角高频子带上的小波系数估计噪声标准差。得到通用阈值后，利用软阈值函数对小波系数进行阈值收缩，即有

式中，sgn（·）为符号函数，下标“＋”表示保持正值不变，将负值置零。最后，将收缩后的小波系数进行小波重构得到去噪后的图像。

VisuShrink方法实现简单，能够有效地去除噪声。然而，由于通用阈值较大，该方法容易导致图像的过度平滑，使细节保持效果不好。而且VisuShrink方法仅对高频子带的小波系数进行阈值收缩，而保留低频子带的小波系数，因而难以选择合适的小波变换尺度。当小波变换尺度较小时，低频子带中的噪声无法滤除，小波变换尺度过大时，又容易破坏细节。

2.2 NeighShrink方法

小波变换具有能力集中和系数聚簇的特点。也就是说，一个小波系数若包含有用信息，具有较大幅值时，则其邻域的小波系数包含有用信息的可能性也很大。根据小波变换的这一性质，文献［17］提出了利用邻域信息的小波系数收缩方法。文献［18］将其扩展到二维，用于图像去噪，提出了NeighShrink方法。该方法以当前要处理的小波系数为中心，选取大小合适的窗口，令

小波系数的收缩方法为

式中，T为通用阈值，计算方法见式（8）。

由此可见，与VisuShrink方法在阈值处理时仅考虑当前待处理的小波系数不同，NeighShrink方法充分利用了其周围小波系数的分布特点并将其引入收缩策略，具有比VisuShrink方法较低的阈值，能够更好地保持图像的细节。选择的邻域窗口越大，收缩阈值越低，细节保持效果更好，但去噪能力相应会下降。为保证噪声滤除的效果，通常选取的窗口大小为3×3或者5×5。

3 结合局部频率估计的小波域干涉相位滤波方法

根据上节对两种阈值收缩方法的分析可以看出，VisuShrink方法的阈值较高，去噪能力较好，而NeighShrink方法考虑了当前小波系数的邻域信息，阈值低于VisuShrink方法，细节保持能力更强。根据以上特点，本文希望结合VisuShrink和NeighShrink方法进行干涉相位滤波，即对于干涉相位有用信息所在频率子带的小波系数利用NeighShrink方法进行阈值收缩，而对于其他主要包含噪声的频率子带的小波系数则利用VisuShrink方法进行阈值收缩，这样可以在尽可能地滤除噪声的同时保持干涉条纹的细节信息不被破坏。因此，在进行小波系数阈值收缩之前，需要首先确定有用信息所在的子带。局部频率估计可以获得干涉条纹的频率范围，根据式（7）给出的小波变换后不同子带的频率范围即可确定干涉条纹信息所在的子带。因此本节首先介绍局部频率估计的方法，之后给出本文提出的滤波处理算法流程。

3.1 局部频率估计

干涉相位局部频率估计的方法有多种。其中，最大似然频率估计方法是应用最为广泛的方法之一，它可以通过二维快速傅里叶变换（fast Fourier transform，FFT）实现。在复干涉相位图的一个小窗口内，地形的坡度可以认为是一个常数，即复干涉相位图中仅包含一个主要频率，可以用一个二维复正弦信号表示。在一个矩形窗内，干涉相位可写为

φ（m＋k，n＋l）＝2πk·fa＋2πl·fr＋φ（m，n）（13）式中，φ（m，n）为窗口中心像素的相位；k和l为窗口中的像素相对于中心像素沿方位向和距离向上的偏移；fa和fr分别为方位向和距离向的干涉相位频率。最大似然频率估计方法按下式进行

为保证窗口内信号的平稳性，窗口的选择不能过大，这样在具体实现过程中，由于二维FFT的量化效应，会使估计结果的方差较大。因此，在不增加窗口大小的前提下，为了提高估计精度，可以通过对窗口内的信号补零后，再进行二维FFT，从而达到减小频率采样间隔的目的。根据二维局部频率估计方差的克拉美罗界［19］，可以利用式（15）近似计算出补零后的窗口大小

式中，La，Lr为补零后方位向和距离向的窗口长度，rs／n为干涉相位图的信噪比，不等式左侧为离散傅里叶变换后输出频率分辨率的平方，右侧为局部频率估计方差的克拉美罗界。当信噪比为0 dB，截取干涉相位图的窗口大小为9×9时，根据式（15）可计算出需要进行256×256点的补零FFT。然而，对干涉相位图中的每一像素都进行如此多点的补零FFT是非常耗时的。实际上，只有在频谱的主瓣内需要做内插，这样可以大大提高计算效率。而Chirp-Z变换（Chirp-Z transform，CZT）能够在任意指定的频率范围内调整频率采样间隔。已知x（n），0≤n≤N－1，则它的CZT为

式中，M表示欲分析的复频谱点数，不一定等于N；A和W均为任意复数，A为复频谱的起始取样点，W为取样点之间的间隔。因此，可以首先利用16×16点的补零FFT找出

这样就达到了256×256点补零FFT的估计精度，但相比直接补零FFT而言，计算效率提高了约30倍。

3.2 滤波算法流程

根据上述分析，本文提出了结合小波阈值收缩和局部频率估计的干涉相位滤波方法，算法流程如图2所示。

图2 本文提出相位滤波方法流程图

具体实现步骤如下：

步骤1利用最大似然方法对干涉相位图进行局部频率估计，得到干涉条纹所在的频率范围。

步骤2对式（2）中的实部和虚部分别利用DWT进行分解，得到各自的小波系数。根据式（7）和步骤1的局部频率估计结果，可以判断出干涉相位中有用信息所在的子带。

步骤3对于有用信息所在子带的小波系数，利用NeighShrink方法进行阈值收缩；而对于其他主要为噪声的子带中的小波系数，则利用VisuShrink方法进行阈值收缩。

步骤4将阈值收缩后的小波系数，利用IDWT进行重构，分别得到滤波后的复干涉相位的实部和虚部。

步骤5计算出滤波后的干涉相位为＝arctan（I︵m／）。

4 实验结果分析

4.1 仿真实验

本小节利用Matlab中的peaks函数生成理想的干涉相位图如图3（a）所示，数据大小为800×800。根据单视时的干涉相位概率密度函数［2］，仿真相干系数为0.85时的干涉相位如图3（b）所示。干涉相位方位向和距离向的局部频率估计结果分别如图3（c）和图3（d）所示，从中得到干涉相位图的频率范围为［0，0.133π］×［0，0.094π］。实验中，小波变换的尺度为5，由于第i个尺度分解后的低频子带范围为［0，2－iπ］，因此本文方法对前两个尺度分解后的高频子带中的小波系数均利用VisuShrink方法进行阈值收缩，而其他子带的小波系数利用NeighShrink方法进行阈值收缩，滤波结果如图3（i）所示。图3（e）～图3（h）分别给出了WInPF滤波［14］、Lee滤波［7］、单独用VisuShrink方法及NeighShrink方法阈值收缩的滤波结果。为了更清楚地显示滤波后干涉条纹细节保持效果，将图3（e）～图3（i）中白色方框的区域进行放大显示，如图3（j）～图3（n）所示。可以看出，WInPF和Lee滤波后，干涉相位图中还残余一定的噪声。VisuShrink方法在条纹密集处一定程度上会破坏条纹的结构，而NeighShrink方法的细节保持能力较好但残余了部分噪声。相比而言，本文方法获得了最好的滤波效果。为了进一步量化评价不同方法的滤波效果，表1给出了不同方法滤波后的残差点数目以及滤波后干涉相位与理想干涉相位之间的均方根误差（root-meam-square error，RMSE），可以看出本文方法较其他方法均有较大的改善。

另外，表1还给出了不同滤波方法处理该仿真数据的运行时间。可以看出，本文方法由于需进行局部频率估计比WInPF滤波、VisuShrink滤波及NeighShrink滤波运行时间长，但是相比Lee滤波而言运算效率有较大提高，能够处理较大维数的数据。

表1 不同滤波方法滤波效果及运行时间对比

图3 仿真干涉相位滤波结果

4.2 实测数据实验

本小节利用SIR-C／X-SAR获取的意大利Etna火山口干涉相位进行相位滤波实验。实际的干涉相位如图4（a）所示，数据大小为512×512。可以看出，该干涉相位图条纹非常密集，而且信噪比较低。分别利用WInPF滤波、Lee滤波、VisuShrink滤波、NeighShrink滤波以及本文方法进行处理，滤波结果如图4（b）～图4（f）所示。可以看出，WInPF滤波和Lee滤波的去噪效果较差，Neigh-Shrink方法去噪效果稍好，而VisuShrink方法对条纹的细节信息破坏严重。本文滤波方法的小波变换尺度仍为5，根据图4（g）和图4（h）的频率估计结果得到干涉相位图的频率范围为［0，0.273π］×［0，0.438π］，因而对第一个尺度分解后的高频子带中的小波系数用VisuShrink方法处理，而其他子带的小波系数用NeighShrink方法处理。从图4（f）可以看出，条纹密集处的干涉相位仍得到了很好的保持。表1显示了不同方法滤波后干涉相位的残差点数，可以看出本文方法滤波后残差点数目最少。

图4 实际干涉相位滤波结果

5 结论

本文提出了一种局部频率估计与小波阈值收缩相结合的干涉相位滤波新方法。该方法针对VisuShrink和NeighShrink两种小波阈值收缩算法各自的优点和缺陷，将二者结合使用进行干涉相位滤波。该方法的关键在于通过最大似然估计得到干涉相位的局部条纹频率，根据估计结果，对估计出的频率范围所在子带的小波系数利用NeighShrink方法进行阈值收缩，而对其他子带的小波系数则利用VisuShrink方法进行阈值收缩，从而可以有效地滤除干涉相位噪声，并保证干涉条纹结构尽量不被破坏。本文方法继承了基于小波变换的滤波方法效率高的优点，可以用于较大维数数据的处理。仿真和实测数据实验验证了本文方法的有效性。

［1］Rosen P A，Hensley S，Joughin I R，et al.Synthetic aperture radar interferometry［J］.Proceedings of the IEEE，2000，88（3）：333 -382.

［2］Hanssen R F.Radar interferometry，data interpretation and error analysis［M］.The Netherlands：Kluwer Academic Publishers，2001.

［3］Lanari R.Generation of digital elevation models by using SIR-C／X-SAR multifrequency two-pass interferometry：the Etna case study［J］.IEEE Trans.on Geoscience and Remote Sensing，1996，34（5）：1097- 1114.

［4］Suo Z，Li Z，Bao Z.A new strategy to estimate local fringe frequencies for InSAR phase noise reduction［J］.IEEE Geoscience and Remote Sensing Letters，2010，7（4）：771- 775.

［5］Cai B，Liang D，Dong Z.A new adaptive multiresolution noise-filtering approach for SAR interferometric phase images［J］.IEEE Geoscience and Remote Sensing Letters，2008，5（2）：266- 270.

［6］Wu N，Feng D，Li J.A locally adaptive filter of interferometric phase images［J］.IEEE Geoscience and Remote Sensing Letters，2006，3（1）：73- 77.

［7］Lee J S，Papathanassiou K P，Ainsworth T L，et al.A new technique for noise filtering of SAR interferometric phase images［J］.IEEE Trans.on Geoscience and Remote Sensing，1998，36（5）：1456- 1465.

［8］Buades A，Coll B，Morel J M.Nonlocal image and movie denoising［J］.International Journal of Computer Vision，2008，76（2）：123- 140.

［9］Parrilli S，Poderico M，Angelino C V，et al.A non-local SAR image denosing algorithm based on LLMMSE wavelet shrinkage［J］.IEEE Trans.on Geoscience and Remote Sensing，2012，50（2）：606- 616.

［10］Deledalle C A，Denis L，Tupin F.NL-InSAR：non-local interferogram estimation［J］.IEEE Trans.on Geoscience and Remote Sensing，2011，49（4）：1441- 1452.

［11］Lin X，Li F F，Meng D D，et al.Nonlocal SAR interferometric phase filtering through higher order singular value decomposition［J］.IEEE Geoscience and Remote Sensing Letters，2015，12（4）：806- 810.

［12］Goldstein R M，Werner C L.Radar interferogram filtering for geophysical application［J］.Geophysical Research Letters，1998，25（21）：4035- 4038.

［13］Li Z，Bao Z，Li H，et al.Image autocoregistration and InSAR interferogram estimation using joint subspace projection［J］.IEEE Trans.on Geoscience and Remote Sensing，2006，44（2）：288- 297.

［14］Lopez-Martinez C，Fabregas X.Modeling and reduction of SAR interferometric phase noise in the wavelet domain［J］.IEEE Trans. on Geoscience and Remote Sensing，2002，40（12）：2553- 2566.

［15］Wang P，Wang Y F，Zhang B C，et al.InSAR interferogram filtering based on Bayesian threshold in stationary wavelet domain［J］.Journal of Electronics&Information Technology，2007，29（11）：2706- 2710.（汪沛，王岩飞，张冰尘，等.基于贝叶斯门限的静态小波域干涉相位图滤波［J］.电子与信息学报，2007，29（11）：2706- 2710.）

［16］Donoho D L，Johnstone I M.Ideal spatial adaptation by wavelet shrinkage［J］.Biometrika，1994，81（3）：425- 455.

［17］Cai T T，Silverman B W.Incorporating information on neighboring coefficients into wavelet estimation［J］.The Indian Journal of Statistics，Series B，2001，63（2）：127- 148.

［18］Chen G Y，Bui T D，Krzyak A.Image denoising with neighbor dependency and customized wavelet and threshold［J］.Pattern Recognition，2005，38（1）：115- 124.

［19］Apagnolini U.2-D phase unwrapping and instantaneous frequency estimation［J］.IEEE Trans.on Geoscience and Remote Sensing，1995，33（3）：579- 589.

Novel InSAR phase filtering method in wavelet domain integrated with local frequency estimation

LI Fang-fang1，2，LIN Xue1，2，HU Dong-hui1，2，DING Chi-biao1，2
（1.Institute of Electronics，Chinese Academy of Sciences，Beijing 100190，China；2.Key Laboratory of Technology in Geo-spatial Information Processing and Application System，Chinese Academy of Sciences，Beijing 100190，China）

Phase noise in the interferogram directly affects the result of phase unwrapping and the precision of interferometirc measurement.A novel approach combining the local frequency estimation with wavelet transform is presented to reduce interferometric phase noise for interferometric synthetic aperture radar（InSAR）.According to the local frequency estimation result，the subband of wavelet coefficients containing useful information can be found out.Based on the fact that the VisuShrink method has effective noise reduction ability and the NeighShrink method performs well in detail preservation，the wavelet coefficients containing useful information are shrunk by the NeighShrink method while the others are shrunk by the VisuShrink method.As a result，the noisy interferogram can be efficiently filtered without destroying the fringes.The experiments using the simulated data and the real data demonstrate the validity of the proposed method.

interferometric synthetic aperture radar（InSAR）；phase filtering；wavelet transform；local frequency estimation；threshold shrinkage

TP 391

10.3969／j.issn.1001-506X.2015.12.09

李芳芳（1986- ），女，助理研究员，博士，主要研究方向为干涉SAR信号处理。

E-mail：iecas_lab8＠126.com

林雪（1986-），女，博士研究生，主要研究方向为干涉SAR信号处理。

E-mail：linxue862002＠163.com

胡东辉（1970-），男，副研究员，硕士，主要研究方向为SAR／ISAR成像。

E-mail：dhhu＠mails.ie.ac.cn

丁赤飚（196-9- ），男，研究员，博士，主要研究方向为遥感信息处理和应用系统。

E-mail：cbding＠mails.ie.ac.cn

1001-506X（2015）12-2719-06