CN104459633A - 结合局部频率估计的小波域InSAR干涉相位滤波方法 - Google Patents

Abstract

本发明提供了一种结合局部频率估计的小波域InSAR干涉相位滤波方法。该小波域InSAR干涉相位滤波方法利用局部频率估计实现了对复干涉相位的小波系数中有用信息子带和噪声子带的区分，利用通用阈值收缩和邻域阈值收缩两种方法分别具有去噪效果好和细节保持能力强的特点，对有用信息所在子带的小波系数进行邻域阈值收缩，而对其它子带的小波系数则进行通用阈值收缩，从而尽可能的滤除噪声，同时保持干涉条纹的细节信息不被破坏，实现高精度的干涉相位滤波，为高精度的干涉测量提供了条件。

Description

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

技术领域

本发明涉及电子信息技术雷达技术领域，尤其涉及一种结合局部频率估计的小波域InSAR干涉相位滤波方法。

背景技术

干涉合成孔径雷达(Interferometric Synthetic Aperture Radar，InSAR)是利用合成孔径雷达(SAR)两个通道的干涉相位信息提取地表的高程信息或变化信息，将SAR的测量拓展到三维空间，具有全天时、全天候、高精度的特点，因此在地形测绘、冰川研究、海洋测绘以及地面沉降监测等多个领域都有广泛的应用。

干涉测量的精度和可靠性在很大程度上取决于干涉相位图的质量。然而，在实际系统中，受热噪声去相干、时间去相干、基线去相干、配准误差等多种去相干因素的影响，干涉相位图不可避免的存在相位噪声。低质量的干涉相位将会影响后续的干涉相位解缠及高程反演的准确性。因此，在相位解缠前必须对干涉相位进行滤波，从而获取高质量的干涉相位图。

目前干涉相位的滤波方法可以大致分为空间域滤波和变换域滤波两类。圆周期均值滤波或中值滤波是一种最基本的空间域滤波方法(见参考文献1)，它实现简单，但滤波窗口大小不好确定，在条纹密集时容易破坏相位细节，降低分辨率。变换域滤波方法在实际中应用更为广泛，如Goldstein滤波(见参考文献2)，然而，该方法受分块大小和滤波参数的影响较大，在信噪比很低时，滤波效果较差；小波变换由于其良好的时频分析特性和多分辨率特性，也可以应用于干涉相位滤波中(见参考文献3)，该方法通过增大小波系数中的信号成分来实现滤波，能够很好地保持干涉条纹的细节信息，并在一定程度上提高了图像的信噪比，但由于对噪声信息没有进行有效的抑制，使得去噪效果较差。因此，为满足InSAR应用对干涉相位精度的要求，有必要进一步研究能有效去除噪声并保持相位细节的变换域滤波方法。

参考文献：

[1]R.Lanari.Generation of digital elevation models by usingSIR-C/X-SAR multifrequency two-pass interferometry：The Etna case study.IEEE Transactions on Geoscience and Remote Sensing，1996，34(5)：1097-1114.

[2]R.M.Goldstein，C.L.Werner.Radar Interferogram filtering forGeophysical Application.Geophysical Research Letters.1998，25(21)：4035-4038.

[3]Lopez-Martinez C，Fabregas X.Modeling and reduction of SARinterferometric phase noise in the wavelet domain[J].IEEE Trans.onGeoscience and Remote Sensing，2002，40(12)：2553-2566

发明内容

(一)要解决的技术问题

为解决上述问题，本发明提供了一种结合局部频率估计的小波域InSAR干涉相位滤波的方法，以克服现有的小波域干涉相位滤波方法不能兼顾去噪和细节保持的缺点，从而提高干涉相位的精度。

(二)技术方案

根据本发明的一个方面，提供了一种结合局部频率估计的小波域InSAR干涉相位滤波方法。该小波域InSAR干涉相位滤波方法包括：步骤A：将InSAR的干涉相位φ变换到复数域e^jφ，分别取复数域干涉相位e^jφ的实部和虚部；步骤B：对干涉相位φ进行局部频率估计，得到干涉相位所在的频率范围；步骤C：对复数域干涉相位的实部和虚部，分别进行尺度为s的小波分解，得到不同子带的小波系数的频率范围，其中，m，n为小波系数的位置，i为分解的尺度，其范围为1～s；步骤D：对复数域干涉相位的实部和虚部，分别根据干涉相位所在的频率范围和不同子带小波系数的频率范围，分别确定有用信息和噪声所在的子带；步骤E：对复数域干涉相位的实部和虚部，分别对噪声所在子带的小波系数分别进行通用阈值收缩处理；步骤F：对复数域干涉相位的实部和虚部，分别对有用信息所在子带的小波系数分别进行邻域阈值收缩处理；步骤G：对复数域干涉相位的实部和虚部，分别将通用阈值收缩处理处理后的噪声所在子带的小波系数和邻域阈值收缩处理后的有用信息所在子带的小波系数共同进行小波重构，得到滤波后的复数域干涉相位的实部和虚部；以及步骤H：由滤波后的复数域干涉相位的实部和虚部得到InSAR滤波后的干涉相位。

(三)有益效果

本发明利用局部频率估计实现了对复干涉相位的小波系数中有用信息子带和噪声子带的区分，利用通用阈值收缩和邻域阈值收缩两种方法分别具有去噪效果好和细节保持能力强的特点，对有用信息所在子带的小波系数进行邻域阈值收缩，而对其它子带的小波系数则进行通用阈值收缩，从而尽可能的滤除噪声，同时保持干涉条纹的细节信息不被破坏，实现高精度的干涉相位滤波，为高精度的干涉测量提供了条件。

附图说明

图1为根据本发明实施例结合局部频率估计的小波域InSAR干涉相位滤波方法的流程图；

图2为意大利Etna火山实测的星载InSAR干涉相位图；

图3A～图3D为分别利用圆周期均值滤波，Goldstein滤波，Lopez-Martinez C提出的小波滤波以及本实施例方法滤波后的干涉相位图

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。需要说明的是，在附图或说明书描述中，相似或相同的部分都使用相同的图号。附图中未绘示或描述的实现方式，为所属技术领域中普通技术人员所知的形式。另外，虽然本文可提供包含特定值的参数的示范，但应了解，参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应的值。

本发明将将干涉相位变换到小波域后，利用局部频率估计判断出小波系数中包含有用信息的子带，对有用信息所在子带的小波系数利用邻域阈值进行收缩处理，而对其他子带的小波系数利用通用阈值进行收缩处理，从而尽可能的滤除噪声，同时不损失干涉条纹的细节信息，从而为InSAR干涉相位滤波提供了一种有效的方法。

根据本发明的一个方面，提供了一种结合局部频率估计的小波域InSAR干涉相位滤波方法。图1为根据本发明实施例结合局部频率估计的小波域InSAR干涉相位滤波方法的流程图。如图1所示，本实施例结合局部频率估计的小波域InSAR干涉相位滤波方法包括：

步骤A：将InSAR预处理后得到的干涉相位φ变换到复数域e^jφ＝cosφ+jsinφ中，分别取复数域干涉相位e^jφ的实部Re{e^jφ}＝cosφ，虚部Im{e^jφ}＝sinφ；

此处的预处理包括：双通道InSAR数据的成像处理，配准处理、复数据的共轭相乘等等。这些均为本领域内的公知常识，此处不再详细说明。

步骤B：对复数域干涉相位e^jφ进行局部频率估计，得到干涉相位所在的频率范围；

该步骤B具体包括：

子步骤B1：对干涉相位φ的每一像素，对其进行局部频率估计，得到其方位向干涉相位频率和距离向干涉相位频率具体如下：

子分步骤B1a：以该像素为中心取(2M+1)×(2N+1)的估计窗口，该估计窗口的相位模型可表示为：

其中，为估计窗口中心点的相位，窗口内像素相对窗口中心的位移为k、l，f_a，f_r分别为该窗口沿方位向和距离向的干涉相位频率。

本实施例中，M和N的取值分别为4和4，但本发明并不以此为限。在本发明的其他实施例中，该M和N的取值满足：2≤M≤10，2≤N≤10即可。

子分步骤B1b：通过最大化对应估计窗口的相位模型的代价函数来估计该像素沿方位向的干涉相位频率和距离向的干涉相位频率

子步骤B2：对所有像素估计得到的方位向干涉相位频率分别取最大值和最小值为和距离向干涉相位频率分别取最大值和最小值为和得到整幅干涉相位图所在的频率范围为

步骤C：对复数域干涉相位e^jφ的实部和虚部，分别进行尺度为s的小波分解，得到不同子带的小波系数的频率范围，其中，m，n为小波系数的位置，i为分解的尺度，其范围为1～s；

在本领域中，已有很多用于小波分解的函数，例如：Symlets函数、Daubechies函数、Coiflets函数，均可以应用到本发明中。此外，尺度为s的范围介于2至8之间，而不局限于本实施例中的5。

不同子带小波系数的频率范围为：

$\begin{array}{c}{a}_s^{\mathrm{LL}}\∈[0,2^{-s}\mathrm{\π})\×[0,2^{-s}\mathrm{\π})\\ \begin{array}{cc}{d}_i^{\mathrm{LH}}\∈[0,2^{-i}\mathrm{\π})\×{[2}^{-i}\mathrm{\π},2^{1-i}\mathrm{\π})& i=1,...,s\end{array}\\ \begin{array}{cc}{d}_i^{\mathrm{HL}}\∈[2^{-i}\mathrm{\π},2^{1-i}\mathrm{\π})\×[0,2^{-i}\mathrm{\π})& i=1,...,s\end{array}\\ \begin{array}{cc}{d}_i^{\mathrm{HH}}\∈[2^{-i}\mathrm{\π},2^{1-i}\mathrm{\π})\×[2^{-i}\mathrm{\π},2^{1-i}\mathrm{\π})& i=1,...,s\end{array}\end{array}---\left(3\right)$

其中，为子带对应干涉相位图的低频部分，和分别为对应干涉相位图在各个尺度的垂直方向、水平方向和对角线方向的高频部分。

步骤D：对复数域干涉相位的实部和虚部，分别根据干涉相位所在的频率范围和不同子带小波系数的频率范围，确定有用信息和噪声所在的子带；

具体方法为：分别判断不同子带小波系数的频率范围 与是否有交集，如果有交集，则判断该子带即包含有用信息，该子带为有用信息所在的子带；如果交集为空，则判断该子带是噪声所在子带。

步骤E：对复数域干涉相位的实部和虚部，分别对噪声所在子带的小波系数进行通用阈值收缩处理；

其中，通用阈值压缩处理的具体计算如下：

其中，为阈值收缩后的小波系数，sgn(·)为符号函数，下标+表示保持正值不变，将负值置零。T为通用阈值，由下式获得：

$T=\mathrm{\σ}\sqrt{2\log N}---\left(5\right)$

其中，N为干涉相位图包含的全部像素数，σ为噪声标准差，按下式进行估计：

$\mathrm{\σ}=\frac{\mathrm{Median}\left(\right|c_{m,n}^i\left|\right)}{0.6745}---\left(6\right)$

其中，Median()表示中值。

通用阈值收缩具有去噪效果好的特点，能够有效的滤除干涉相位噪声。

步骤F：对复数域干涉相位的实部和虚部，分别对有用信息所在子带的小波系数进行邻域阈值收缩处理；

其中，邻域阈值收缩处理的具体计算如下：

首先，以当前要处理的小波系数为中心，选取大小合适的窗口W，令

${S_{m,n}^i}^2=\underset{(m,n)\∈W}{\mathrm{\Σ}}{c_{m,n}^i}^2---\left(7\right)$

则阈值收缩后的小波系数为：

其中，下标+表示保持正值不变，将负值置零，T为通用阈值，按下式计算：

$T=\mathrm{\σ}\sqrt{2\log N}---\left(9\right)$

其中，N为干涉相位图包含的像素数，σ为噪声标准差，按下式进行估计：

$\mathrm{\σ}=\frac{\mathrm{Median}\left(\right|c_{m,n}^i\left|\right)}{0.6745}---\left(10\right)$

邻域阈值具有细节保持能力强的特点，有利于保护干涉相位图中的条纹结构不被破坏。

步骤G：对复数域干涉相位的实部和虚部，分别将步骤E得到的通用阈值收缩处理后的噪声所在子带的小波系数和步骤F得到的邻域阈值收缩处理后的有用信息所在子带的小波系数共同进行小波重构，得到滤波后的复数域干涉相位的实部和虚部

关于小波重构的具体计算过程，已经为本领域技术人员所熟知，此处不再赘述。

步骤H：对滤波后的复数域干涉相位的实部和虚部进行如下计算，得到InSAR滤波后的干涉相位

下面通过实测数据验证了本实施例小波域InSAR干涉相位滤波方法的有效性。图2为意大利Etna火山实测的星载InSAR干涉相位图。图3A～图3D为分别利用圆周期均值滤波，Goldstein滤波，Lopez-Martinez C提出的小波滤波以及本实施例方法滤波后的干涉相位图，其中图3A为圆周期均值滤波结果，滤波窗口为5×5，图3B为Goldstein滤波结果，分块大小为32×32，图3C为小波滤波结果，图3D为本发明方法滤波结果。可以看出，本发明方法滤波后的干涉相位最为平滑，且条纹密集处的结构保持良好。滤波后干涉相位图中的残差点数目可用于评价去噪效果的好坏，计算出滤波前后的残差点数目如表1。从表中同样可以看出本发明方法滤波后剩余的残差点数目最少，去噪效果最好。

表1 干涉相位滤波前后残差点数目比较

此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列且可根据所需设计而变化或重新安排，例如，上述实施例中的步骤B和步骤C为次序可以更换，步骤E和步骤F中的各个步骤可以顺序或并行进行。

至此，已经结合附图对本发明实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明结合局部频率估计的小波域InSAR干涉相位滤波方法有了清楚的认识。

需要说明的是，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换，例如：子分步骤B2中干涉相位频率估计可以采用BF(BeamForming)或MUSIC(Multi-signal Classification)等谱估计的方法来代替。

综上所述，本发明利用局部频率估计实现了对复干涉相位的小波系数中有用信息子带和噪声子带的区分，对有用信息所在子带的小波系数进行邻域阈值收缩，而对其它子带的小波系数则进行通用阈值收缩，从而尽可能的滤除噪声，同时保持干涉条纹的细节信息不被破坏，实现高精度的干涉相位滤波，为高精度的干涉测量提供了条件。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种结合局部频率估计的小波域InSAR干涉相位滤波方法，其特征在于，包括：

步骤A：将InSAR的干涉相位φ变换到复数域e^jφ，分别取复数域干涉相位e^jφ的实部和虚部；

步骤B：对干涉相位φ进行局部频率估计，得到干涉相位所在的频率范围；

步骤C：对复数域干涉相位的实部和虚部，分别进行尺度为s的小波分解，得到不同子带的小波系数的频率范围，其中，m，n为小波系数的位置，i为分解的尺度，其范围为1～s；

步骤D：对复数域干涉相位的实部和虚部，分别根据干涉相位所在的频率范围和不同子带小波系数的频率范围，分别确定有用信息和噪声所在的子带；

步骤E：对复数域干涉相位的实部和虚部，分别对噪声所在子带的小波系数分别进行通用阈值收缩处理；

步骤F：对复数域干涉相位的实部和虚部，分别对有用信息所在子带的小波系数分别进行邻域阈值收缩处理；

步骤G：对复数域干涉相位的实部和虚部，分别将通用阈值收缩处理处理后的噪声所在子带的小波系数和邻域阈值收缩处理后的有用信息所在子带的小波系数共同进行小波重构，得到滤波后的复数域干涉相位的实部和虚部；以及

步骤H：由滤波后的复数域干涉相位的实部和虚部得到InSAR滤波后的干涉相位。

2.根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤B包括：

子步骤B1：对干涉相位φ的每一像素，对其进行局部频率估计，得到其方位向干涉相位频率和距离向干涉相位频率

子步骤B2：对所有像素估计得到的方位向干涉相位频率分别取最大值和最小值为和距离向干涉相位频率分别取最大值和最小值为和得到整幅干涉相位图所在的频率范围为

3.根据权利要求2所述的小波域InSAR干涉相位滤波方法，其特征在于，所述子步骤B1中，对干涉相位φ的每一像素得到其方位向干涉相位频率和距离向干涉相位频率包括：

子分步骤B1a：以该像素为中心取(2M+1)×(2N+1)的估计窗口，该估计窗口的相位模型为：

其中，为估计窗口中心点的相位，k、l为窗口内像素相对窗口中心的位移，f_a，f_r分别为该窗口沿方位向和距离向的干涉相位频率；

子分步骤B1b：通过最大化估计窗口的相位模型的代价函数来估计该像素沿方位向的干涉相位频率和距离向的干涉相位频率

4.根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤C中，不同子带小波系数的频率范围为：

$a_s^{\mathrm{LL}}\∈[0,2^{-s}\mathrm{\π})\×[0,2^{-s}\mathrm{\π})$

$d_i^{\mathrm{LH}}\∈[0,2^{-i}\mathrm{\π})\×[2^{-i}\mathrm{\π},2^{1-i}\mathrm{\π})i=1,...,s$

$d_i^{\mathrm{HL}}\∈[2^{-i}\mathrm{\π},2^{1-i}\mathrm{\π})\×[0,2^{-i}\mathrm{\π})i=1,...,s$

$d_i^{\mathrm{HH}}\∈[2^{-i}\mathrm{\π},2^{1-i}\mathrm{\π})\×[2^{-i}\mathrm{\π},2^{1-i}\mathrm{\π})i=1,...,s$

其中，为子带对应干涉相位图的低频部分，和分别为对应干涉相位图在各个尺度的垂直方向、水平方向和对角线方向的高频部分。

5.根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤C中：用于小波分解的函数为Symlets函数、Daubechies函数或Coiflets函数，所述尺度s介于2至8之间。

6.根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤D根据干涉相位所在的频率范围和不同子带小波系数的频率范围，分别确定有用信息和噪声所在的子带包括：

分别判断不同子带小波系数的频率范围与是否有交集，如果有交集，则判断该子带为有用信息所在的子带；如果交集为空，判断该子带是噪声所在子带。

7.根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤E中，按照下式对噪声所在子带的小波系数进行通用阈值收缩处理：

其中，为阈值收缩后的小波系数，sgn(·)为符号函数，下标+表示保持正值不变，将负值置零；T为通用阈值。

8.根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤F中对有用信息所在子带的小波系数分别进行邻域阈值收缩处理包括：

以当前要处理的小波系数为中心，选取大小合适的窗口W，令 ${S_{m,n}^i}^2=\underset{(m,n)\∈W}{\mathrm{\Σ}}{c_{m,n}^i}^2;$

则阈值收缩后的小波系数为：其中，下标+表示保持正值不变，将负值置零，T为通用阈值。

9.根据权利要求7或8所述的小波域InSAR干涉相位滤波方法，其特征在于，所述通用阈值T按照下式计算：

$T=\mathrm{\σ}\sqrt{2\log N}$

其中，N为干涉相位图包含的像素数，σ为噪声标准差，按下式进行估计： $\mathrm{\σ}=\frac{\mathrm{Median}\left(\right|c_{m,n}^i\left|\right)}{0.6745}.$

10.根据权利要求1至8中任一项所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤H中，按照下式得到InSAR滤波后的干涉相位

$\widehat{\mathrm{\φ}}=\arctan (\widehat{\mathrm{Im}}/\widehat{\mathrm{Re}})$

其中，和分别为滤波后的复数域干涉相位的实部和虚部。