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

摘要

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

说明书

技术领域

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

背景技术

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

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

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

参考文献：

[1]R.Lanari.Generation of digital elevation models by using  SIR-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 for  Geophysical Application.Geophysical Research Letters.1998，25(21)： 4035-4038.

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

发明内容

(一)要解决的技术问题

为解决上述问题，本发明提供了一种结合局部频率估计的小波域 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。

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

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

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

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

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

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

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

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

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

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

$\sigma =\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)\in W}{\Sigma }{c_{m,n}^i}^2---\left(7\right)$

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

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

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

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

$\sigma =\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(Beam  Forming)或MUSIC(Multi-signal Classification)等谱估计的方法来代替。

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

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