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

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

Info

Publication number
WO2016086699A1
WO2016086699A1 PCT/CN2015/089141 CN2015089141W WO2016086699A1 WO 2016086699 A1 WO2016086699 A1 WO 2016086699A1 CN 2015089141 W CN2015089141 W CN 2015089141W WO 2016086699 A1 WO2016086699 A1 WO 2016086699A1
Authority
WO
WIPO (PCT)
Prior art keywords
interference phase
wavelet
domain
sub
insar
Prior art date
Application number
PCT/CN2015/089141
Other languages
English (en)
French (fr)
Inventor
丁赤飚
李芳芳
雷斌
林雪
胡东辉
仇晓兰
Original Assignee
中国科学院电子学研究所
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 中国科学院电子学研究所 filed Critical 中国科学院电子学研究所
Priority to AU2015358016A priority Critical patent/AU2015358016B2/en
Priority to EP15865192.7A priority patent/EP3229038B1/en
Publication of WO2016086699A1 publication Critical patent/WO2016086699A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques

Definitions

  • the invention relates to the technical field of electronic information technology radar, in particular to a wavelet domain InSAR interference phase filtering method combined with local frequency estimation.
  • Interferometric Synthetic Aperture Radar extracts the elevation information or change information of the surface using the interference phase information of two channels of Synthetic Aperture Radar (SAR), and extends the measurement of SAR into three-dimensional space with all-day time. All-weather, high-precision features. Therefore, it has been widely used in many fields such as topographic mapping, glacial research, ocean mapping and ground subsidence monitoring.
  • the accuracy and reliability of interferometry depends to a large extent on the quality of the interference phase map.
  • various decoherence factors such as thermal noise decoherence, time decoherence, baseline decoherence, registration error, etc.
  • phase interference is inevitable in the interference phase diagram.
  • the low quality interference phase will affect the accuracy of subsequent interference phase unwrapping and elevation inversion. Therefore, the interference phase must be filtered before phase unwrapping to obtain a high quality interference phase map.
  • the filtering method of the interference phase can be roughly classified into two types: spatial domain filtering and transform domain filtering.
  • Circular periodic mean filtering or median filtering is a basic spatial domain filtering method (see Reference 1). It is simple to implement, but the filtering window size is not well defined. When the stripes are dense, it is easy to destroy the phase details and reduce the resolution.
  • the transform domain filtering method is more widely used in practice, such as Goldstein filtering (see Reference 2). However, this method is greatly affected by the block size and filtering parameters. When the signal-to-noise ratio is low, the filtering effect is poor.
  • Wavelet transform can also be applied to interference phase filtering due to its good time-frequency analysis characteristics and multi-resolution characteristics (see Reference 3).
  • This method can achieve filtering by increasing the signal components in the wavelet coefficients.
  • the ground maintains the detailed information of the interference fringes and improves the signal-to-noise ratio of the image to a certain extent, but the noise reduction effect is poor because the noise information is not effectively suppressed. Therefore, in order to meet the requirements of interferometric phase accuracy of InSAR applications, it is necessary to further study the transform domain filtering method which can effectively remove noise and maintain phase details.
  • the invention provides a wavelet domain InSAR interferometric phase filtering method combined with local frequency estimation, so as to overcome the shortcomings of the traditional wavelet domain interference phase filtering method which can not take into consideration the denoising and detail preservation, thereby improving the accuracy of the interference phase.
  • an interference phase filtering method for wavelet domain interferometric synthetic aperture radar InSAR data combined with local frequency estimation comprises the following steps: Step A: transforming the interference phase ⁇ of the InSAR into the complex domain e j ⁇ , respectively taking the real part and the imaginary part of the complex domain interference phase e j ⁇ ; Step B: performing the localization of the interference phase ⁇ Frequency estimation, obtain the frequency range in which the interference phase is located; Step C: Perform wavelet decomposition with scale s on the real and imaginary parts of the complex domain interference phase to obtain wavelet coefficients of different sub-bands The frequency range, where m, n are the positions of the wavelet coefficients, i is the scale of the decomposition, and the range is 1 s; Step D: the real and imaginary parts of the phase of the complex domain interference, according to the frequency range in which the interference phase is located And the frequency range of the wavelet coefficients of different sub-bands respectively
  • the invention utilizes the local frequency estimation to realize the distinction between the useful information sub-band and the noise sub-band in the wavelet coefficients of the complex interference phase, and the two methods of common threshold shrinkage and neighborhood threshold shrinkage respectively have good denoising effect and strong detail retention capability.
  • the feature is that the wavelet coefficients of the sub-bands with useful information are subjected to neighborhood threshold shrinkage, and the wavelet coefficients of other sub-bands are subjected to general threshold shrinkage, so as to filter out noise as much as possible while keeping the details of the interference fringes intact.
  • High-precision interference phase filtering is provided, which provides conditions for high-precision interferometry.
  • FIG. 1 is a flow chart of a wavelet domain InSAR interference phase filtering method in combination with local frequency estimation according to an embodiment of the present invention
  • Figure 2 is a phase diagram of the onboard InSAR interference measured by the Etna volcano in Italy;
  • 3A to 3D are interference phase diagrams respectively using circular period mean filtering, Goldstein filtering, wavelet filtering proposed by Lopez-Martinez C, and filtering by the method of the embodiment.
  • the sub-band containing the useful information in the wavelet coefficient is determined by using the local frequency estimation, and the wavelet coefficient of the sub-band where the useful information is located is shrunk by using the neighborhood threshold. Processing, while wavelet coefficients for other sub-bands are shrunk using common thresholds Processing, so as to filter out noise as much as possible, without losing the details of the interference fringes, thus providing an effective method for InSAR interference phase filtering.
  • a wavelet domain InSAR interference phase filtering method combining local frequency estimation is provided.
  • 1 is a flow chart of a wavelet domain InSAR interference phase filtering method in combination with local frequency estimation according to an embodiment of the present invention. As shown in FIG. 1, the wavelet domain InSAR interference phase filtering method combined with local frequency estimation in this embodiment includes:
  • the preprocessing here includes: imaging processing of dual-channel InSAR data, registration processing, conjugate multiplication of complex data, and the like. These are common knowledge in the art and will not be described in detail herein.
  • Step B performing local frequency estimation on the complex domain interference phase e j ⁇ to obtain a frequency range in which the interference phase is located;
  • the step B specifically includes:
  • Sub-step B1 For each pixel of the interference phase ⁇ , local frequency estimation is performed to obtain its azimuth interference phase frequency Interference phase frequency details as follows:
  • Sub-step B1a taking an estimation window of (2M+i) ⁇ (2N+1) centered on the pixel, the phase model of the estimation window can be expressed as:
  • the displacement of the pixel in the window relative to the center of the window is k, l, f a , f r are the interference phase frequencies of the window along the azimuth direction and the distance direction, respectively.
  • the values of M and N are 4 and 4, respectively, but the invention is not limited thereto.
  • the values of the M and N satisfy: 2 ⁇ M ⁇ 10, and 2 ⁇ N ⁇ 10.
  • Sub-step B1b Estimating the interferometric phase frequency of the pixel along the azimuth by maximizing the cost function of the phase model corresponding to the estimation window Interferometric phase frequency
  • Sub-step B2 azimuth interference phase frequency estimated for all pixels Take the maximum and minimum values respectively with Distance to interference phase frequency Take the maximum and minimum values respectively with The frequency range in which the entire interference phase map is obtained is
  • Step C Perform wavelet decomposition with scale s on the real and imaginary parts of the complex domain interference phase e j ⁇ to obtain wavelet coefficients of different sub-bands Frequency range, where m, n is the position of the wavelet coefficient, i is the scale of the decomposition, and the range is 1 s;
  • the Symlets function There are many functions for wavelet decomposition in the art, such as: the Symlets function, the Daubechies function, and the Coiflets function, all of which can be applied to the present invention. Further, the scale of s is in the range of 2 to 8, and is not limited to 5 in the embodiment.
  • the subband corresponds to the low frequency portion of the interference phase map, with The high frequency portions corresponding to the vertical, horizontal, and diagonal directions of the respective interference phase maps at respective scales.
  • Step D determining the sub-band where the useful information and the noise are located according to the frequency range of the interference phase and the frequency range of the different sub-band wavelet coefficients respectively for the real part and the imaginary part of the complex domain interference phase;
  • the specific method is: respectively determining the frequency range of wavelet coefficients of different sub-bands versus Whether there is an intersection, if there is an intersection, it is judged that the sub-band contains useful information, the sub-band is a sub-band where the useful information is located; if the intersection is empty, it is judged that the sub-band is a sub-band where the noise is located.
  • Step E performing general threshold shrinkage processing on the wavelet coefficients of the sub-bands in which the noise is located, respectively, on the real and imaginary parts of the complex domain interference phase;
  • T is a general threshold and is obtained by:
  • N is the total number of pixels included in the interference phase diagram
  • ⁇ is the noise standard deviation, which is estimated by the following formula:
  • Median () represents the median.
  • the universal threshold shrinkage has the characteristics of good denoising effect, and can effectively filter the interference phase noise.
  • Step F performing real-time threshold processing on the wavelet coefficients of the sub-bands in which the useful information is located, respectively, on the real and imaginary parts of the complex domain interference phase;
  • Wavelet coefficient after threshold shrinkage for:
  • the subscript + means that the positive value remains unchanged, the negative value is set to zero, and T is the general threshold, which is calculated by the following formula:
  • N is the number of pixels included in the interference phase diagram
  • ⁇ is the noise standard deviation, which is estimated by the following formula:
  • the neighborhood threshold has the characteristics of strong detail retention, which is beneficial to protect the fringe structure in the interference phase map from being destroyed.
  • Step G For the real part and the imaginary part of the complex domain interference phase, respectively, the wavelet coefficient of the subband in which the noise after the general threshold shrink processing in step E is obtained, and the useful information after the neighborhood threshold shrinkage process obtained in step F are respectively The wavelet coefficients of the band are jointly reconstructed by wavelet, and the real part of the filtered complex domain interference phase is obtained. And imaginary
  • Step H The real part of the filtered complex domain interference phase And imaginary Perform the following calculation to obtain the interference phase after InSAR filtering.
  • FIG. 3A is a circular period mean filtering result, and the filtering window is 5 ⁇ 5
  • FIG. 3B is a Goldstein filtering result
  • the block size is 32 ⁇ 32
  • FIG. 3C is a wavelet filtering result
  • FIG. 3D is a filtering result of the method of the present invention.
  • the filtered interference phase of the method of the present invention is the smoothest, and the structure of the stripe densely remains well.
  • the number of residual points in the filtered interference phase map can be used to evaluate the effectiveness of the denoising effect.
  • the number of residual points before and after filtering is calculated as shown in Table 1. It can also be seen from the table that the number of residual points remaining after filtering by the method of the present invention is the least, and the denoising effect is the best.
  • steps B and C in the above embodiment are in the order It can be replaced, and the steps in steps E and F can be performed sequentially or in parallel.
  • the interference phase frequency estimation in step B2 may be replaced by a spectral estimation method such as BF (Beam Forming) or MUSIC (Multi-signal Classification).
  • a spectral estimation method such as BF (Beam Forming) or MUSIC (Multi-signal Classification).
  • the present invention utilizes local frequency estimation to realize the distinction between the useful information subband and the noise subband in the wavelet coefficients of the complex interference phase, and performs neighborhood threshold shrinkage on the wavelet coefficients of the subband in which the useful information is located.
  • the wavelet coefficients of other sub-bands are subjected to general threshold shrinkage to filter out noise as much as possible while maintaining the details of the interference fringes without being destroyed, realizing high-precision interference phase filtering, providing conditions for high-precision interferometry. .

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Image Processing (AREA)
  • Radar Systems Or Details Thereof (AREA)

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 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的干涉相位φ变换到复数域ejφ，分别取复数域干涉相位ejφ的实部和虚部；步骤B：对干涉相位φ进行局部频率估计，得到干涉相位所在的频率范围；步骤C：对复数域干涉相位的实部和虚部，分别进行尺度为s的小波分解，得到不同子带的小波系数
Figure PCTCN2015089141-appb-000001
的频率范围，其中，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预处理后得到的干涉相位φ变换到复数域ejφ＝cosφ+j sinφ中，分别取复数域干涉相位ejφ的实部Re{ejφ}＝cosφ，虚部Im{ejφ}＝sinφ；
此处的预处理包括：双通道InSAR数据的成像处理，配准处理、复数据的共轭相乘等等。这些均为本领域内的公知常识，此处不再详细说明。
步骤B：对复数域干涉相位ejφ进行局部频率估计，得到干涉相位所在的频率范围；
该步骤B具体包括：
子步骤B1：对干涉相位φ的每一像素，对其进行局部频率估计，得到其方位向干涉相位频率
Figure PCTCN2015089141-appb-000002
和距离向干涉相位频率
Figure PCTCN2015089141-appb-000003
具体如下：
子分步骤B1a：以该像素为中心取(2M+i)×(2N+1)的估计窗口，该估计窗口的相位模型可表示为：
Figure PCTCN2015089141-appb-000004
其中，
Figure PCTCN2015089141-appb-000005
为估计窗口中心点的相位，窗口内像素相对窗口中心的位移为k、l，fa，fr分别为该窗口沿方位向和距离向的干涉相位频率。
本实施例中，M和N的取值分别为4和4，但本发明并不以此为限。在本发明的其他实施例中，该M和N的取值满足：2≤M≤10，2≤N≤10即可。
子分步骤B1b：通过最大化对应估计窗口的相位模型的代价函数来估计该像素沿方位向的干涉相位频率
Figure PCTCN2015089141-appb-000006
和距离向的干涉相位频率
Figure PCTCN2015089141-appb-000007
Figure PCTCN2015089141-appb-000008
Figure PCTCN2015089141-appb-000009
子步骤B2：对所有像素估计得到的方位向干涉相位频率
Figure PCTCN2015089141-appb-000010
分别取最大值和最小值为
Figure PCTCN2015089141-appb-000011
和
Figure PCTCN2015089141-appb-000012
距离向干涉相位频率
Figure PCTCN2015089141-appb-000013
分别取最大值和最小值为
Figure PCTCN2015089141-appb-000014
和
Figure PCTCN2015089141-appb-000015
得到整幅干涉相位图所在的频率范围为
步骤C：对复数域干涉相位ejφ的实部和虚部，分别进行尺度为s的小波分解，得到不同子带的小波系数
Figure PCTCN2015089141-appb-000017
的频率范围，其中，m，n为小波系数的位置，i为分解的尺度，其范围为1～s；
在本领域中，已有很多用于小波分解的函数，例如：Symlets函数、Daubechies函数、Coiflets函数，均可以应用到本发明中。此外，尺度为s的范围介于2至8之间，而不局限于本实施例中的5。
不同子带小波系数
Figure PCTCN2015089141-appb-000018
的频率范围为：
Figure PCTCN2015089141-appb-000019
Figure PCTCN2015089141-appb-000020
Figure PCTCN2015089141-appb-000021
Figure PCTCN2015089141-appb-000022
其中，
Figure PCTCN2015089141-appb-000023
为子带对应干涉相位图的低频部分，
Figure PCTCN2015089141-appb-000024
和
Figure PCTCN2015089141-appb-000025
分别为对应干涉相位图在各个尺度的垂直方向、水平方向和对角线方向的高频部分。
步骤D：对复数域干涉相位的实部和虚部，分别根据干涉相位所在的频率范围和不同子带小波系数的频率范围，确定有用信息和噪声所在的子带；
具体方法为：分别判断不同子带小波系数的频率范围
Figure PCTCN2015089141-appb-000026
Figure PCTCN2015089141-appb-000027
与
Figure PCTCN2015089141-appb-000028
是否有交集，如果有交集，则判断该子带即包含有用信息，该子带为有用信息所在的子带；如果交集为空，则判断该子带是噪声所在子带。
步骤E：对复数域干涉相位的实部和虚部，分别对噪声所在子带的小波系数进行通用阈值收缩处理；
其中，通用阈值压缩处理的具体计算如下：
Figure PCTCN2015089141-appb-000029
其中，
Figure PCTCN2015089141-appb-000030
为阈值收缩后的小波系数，sgn(·)为符号函数，下标+表示保持正值不变，将负值置零。T为通用阈值，由下式获得：
Figure PCTCN2015089141-appb-000031
其中，N为干涉相位图包含的全部像素数，σ为噪声标准差，按下式进行估计：
Figure PCTCN2015089141-appb-000032
其中，Median()表示中值。
通用阈值收缩具有去噪效果好的特点，能够有效的滤除干涉相位噪声。
步骤F：对复数域干涉相位的实部和虚部，分别对有用信息所在子带的小波系数进行邻域阈值收缩处理；
其中，邻域阈值收缩处理的具体计算如下：
首先，以当前要处理的小波系数
Figure PCTCN2015089141-appb-000033
为中心，选取大小合适的窗口W，令
Figure PCTCN2015089141-appb-000034
则阈值收缩后的小波系数
Figure PCTCN2015089141-appb-000035
为：
Figure PCTCN2015089141-appb-000036
其中，下标+表示保持正值不变，将负值置零，T为通用阈值，按下式计算：
Figure PCTCN2015089141-appb-000037
其中，N为干涉相位图包含的像素数，σ为噪声标准差，按下式进行估计：
Figure PCTCN2015089141-appb-000038
邻域阈值具有细节保持能力强的特点，有利于保护干涉相位图中的条纹结构不被破坏。
步骤G：对复数域干涉相位的实部和虚部，分别将步骤E得到的通用阈值收缩处理后的噪声所在子带的小波系数和步骤F得到的邻域阈值收缩处理后的有用信息所在子带的小波系数共同进行小波重构，得到滤波后的复数域干涉相位的实部
Figure PCTCN2015089141-appb-000039
和虚部
Figure PCTCN2015089141-appb-000040
关于小波重构的具体计算过程，已经为本领域技术人员所熟知，此处不再赘述。
步骤H：对滤波后的复数域干涉相位的实部
Figure PCTCN2015089141-appb-000041
和虚部
Figure PCTCN2015089141-appb-000042
进行如下计算，得到 InSAR滤波后的干涉相位
Figure PCTCN2015089141-appb-000043
Figure PCTCN2015089141-appb-000044
下面通过实测数据验证了本实施例小波域InSAR干涉相位滤波方法的有效性。图2为意大利Etna火山实测的星载InSAR干涉相位图。图3A～图3D为分别利用圆周期均值滤波、Goldstein滤波、Lopez-Martinez C提出的小波滤波以及本实施例方法滤波后的干涉相位图，其中图3A为圆周期均值滤波结果，滤波窗口为5×5，图3B为Goldstein滤波结果，分块大小为32×32，图3C为小波滤波结果，图3D为本发明方法滤波结果。可以看出，本发明方法滤波后的干涉相位最为平滑，且条纹密集处的结构保持良好。滤波后干涉相位图中的残差点数目可用于评价去噪效果的好坏，计算出滤波前后的残差点数目如表1。从表中同样可以看出本发明方法滤波后剩余的残差点数目最少，去噪效果最好。
表1干涉相位滤波前后残差点数目比较
Figure PCTCN2015089141-appb-000045
此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列且可根据所需设计而变化或重新安排，例如，上述实施例中的步骤B和步骤C为次序可以更换，步骤E和步骤F中的各个步骤可以顺序或并行进行。
至此，已经结合附图对本发明实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明结合局部频率估计的小波域InSAR干涉相位滤波方法有了清楚的认识。
需要说明的是，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换，例如：子分步骤B2中干涉相位频率估计可以采用BF(Beam Forming)或MUSIC(Multi-signal Classification)等谱估计的方法来代替。
综上所述，本发明利用局部频率估计实现了对复干涉相位的小波系数中有用信息子带和噪声子带的区分，对有用信息所在子带的小波系数进行邻域阈值收缩， 而对其它子带的小波系数则进行通用阈值收缩，从而尽可能的滤除噪声，同时保持干涉条纹的细节信息不被破坏，实现高精度的干涉相位滤波，为高精度的干涉测量提供了条件。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

  1. 一种结合局部频率估计的小波域InSAR干涉相位滤波方法，其特征在于，包括：
    步骤A：将InSAR的干涉相位φ变换到复数域ejφ，分别取复数域干涉相位ejφ的实部和虚部；
    步骤B：对干涉相位φ进行局部频率估计，得到干涉相位所在的频率范围；
    步骤C：对复数域干涉相位的实部和虚部，分别进行尺度为s的小波分解，得到不同子带的小波系数
    Figure PCTCN2015089141-appb-100001
    的频率范围，其中，m，n为小波系数的位置，i为分解的尺度，其范围为1～s；
    步骤D：对复数域干涉相位的实部和虚部，分别根据干涉相位所在的频率范围和不同子带小波系数的频率范围，分别确定有用信息和噪声所在的子带；
    步骤E：对复数域干涉相位的实部和虚部，分别对噪声所在子带的小波系数分别进行通用阈值收缩处理；
    步骤F：对复数域干涉相位的实部和虚部，分别对有用信息所在子带的小波系数分别进行邻域阈值收缩处理；
    步骤G：对复数域干涉相位的实部和虚部，分别将通用阈值收缩处理处理后的噪声所在子带的小波系数和邻域阈值收缩处理后的有用信息所在子带的小波系数共同进行小波重构，得到滤波后的复数域干涉相位的实部和虚部；以及
    步骤H：由滤波后的复数域干涉相位的实部和虚部得到InSAR滤波后的干涉相位。
  2. 根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤B包括：
    子步骤B1：对干涉相位φ的每一像素，对其进行局部频率估计，得到其方位向干涉相位频率
    Figure PCTCN2015089141-appb-100002
    和距离向干涉相位频率
    Figure PCTCN2015089141-appb-100003
    子步骤B2：对所有像素估计得到的方位向干涉相位频率
    Figure PCTCN2015089141-appb-100004
    分别取最大值和最小值为
    Figure PCTCN2015089141-appb-100005
    和
    Figure PCTCN2015089141-appb-100006
    距离向干涉相位频率
    Figure PCTCN2015089141-appb-100007
    分别取最大值和最小值为
    Figure PCTCN2015089141-appb-100008
    和
    Figure PCTCN2015089141-appb-100009
    得到整幅干涉相位图所在的频率范围为
    Figure PCTCN2015089141-appb-100010
  3. 根据权利要求2所述的小波域InSAR干涉相位滤波方法，其特征在于，所述子步骤B1中，对干涉相位φ的每一像素得到其方位向干涉相位频率
    Figure PCTCN2015089141-appb-100011
    和距离向干涉相位频率
    Figure PCTCN2015089141-appb-100012
    包括：
    子分步骤B1a：以该像素为中心取(2M+1)×(2N+1)的估计窗口，该估计窗口的相位模型为：
    Figure PCTCN2015089141-appb-100013
    其中，
    Figure PCTCN2015089141-appb-100014
    为估计窗口中心点的相位，k、l为窗口内像素相对窗口中心的位移，fa，fr分别为该窗口沿方位向和距离向的干涉相位频率；
    子分步骤B1b：通过最大化估计窗口的相位模型的代价函数来估计该像素沿方位向的干涉相位频率
    Figure PCTCN2015089141-appb-100015
    和距离向的干涉相位频率
    Figure PCTCN2015089141-appb-100016
    Figure PCTCN2015089141-appb-100017
    Figure PCTCN2015089141-appb-100018
  4. 根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤C中，不同子带小波系数
    Figure PCTCN2015089141-appb-100019
    的频率范围为：
    Figure PCTCN2015089141-appb-100020
    Figure PCTCN2015089141-appb-100021
    Figure PCTCN2015089141-appb-100022
    Figure PCTCN2015089141-appb-100023
    其中，
    Figure PCTCN2015089141-appb-100024
    为子带对应干涉相位图的低频部分，
    Figure PCTCN2015089141-appb-100025
    和
    Figure PCTCN2015089141-appb-100026
    分别为对应干涉相位图在各个尺度的垂直方向、水平方向和对角线方向的高频部分。
  5. 根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤C中：用于小波分解的函数为Symlets函数、Daubechies函数或Coiflets函数，所述尺度s介于2至8之间。
  6. 根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤D根据干涉相位所在的频率范围和不同子带小波系数的频率范围，分别确定有用信息和噪声所在的子带包括：
    分别判断不同子带小波系数的频率范围
    Figure PCTCN2015089141-appb-100027
    与
    Figure PCTCN2015089141-appb-100028
    是否有交集，如果有交集，则判断该子带为有用信息所在的子带；如果交集为空，判断该子带是噪声所在子带。
  7. 根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤E 中，按照下式对噪声所在子带的小波系数进行通用阈值收缩处理：
    Figure PCTCN2015089141-appb-100029
    其中，
    Figure PCTCN2015089141-appb-100030
    为阈值收缩后的小波系数，sgn(·)为符号函数，下标+表示保持正值不变，将负值置零；T为通用阈值。
  8. 根据权利要求1所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤F中对有用信息所在子带的小波系数分别进行邻域阈值收缩处理包括：
    以当前要处理的小波系数
    Figure PCTCN2015089141-appb-100031
    为中心，选取大小合适的窗口W，令
    Figure PCTCN2015089141-appb-100032
    则阈值收缩后的小波系数
    Figure PCTCN2015089141-appb-100033
    为：
    Figure PCTCN2015089141-appb-100034
    其中，下标+表示保持正值不变，将负值置零，T为通用阈值。
  9. 根据权利要求7或8所述的小波域InSAR干涉相位滤波方法，其特征在于，所述通用阈值T按照下式计算：
    Figure PCTCN2015089141-appb-100035
    其中，N为干涉相位图包含的像素数，σ为噪声标准差，按下式进行估计：
    Figure PCTCN2015089141-appb-100036
  10. 根据权利要求1至8中任一项所述的小波域InSAR干涉相位滤波方法，其特征在于，所述步骤H中，按照下式得到InSAR滤波后的干涉相位
    Figure PCTCN2015089141-appb-100037
    Figure PCTCN2015089141-appb-100038
    其中，
    Figure PCTCN2015089141-appb-100039
    和
    Figure PCTCN2015089141-appb-100040
    分别为滤波后的复数域干涉相位的实部和虚部。
PCT/CN2015/089141 2014-12-01 2015-09-08 一种结合局部频率估计的小波域InSAR干涉相位滤波方法 WO2016086699A1 (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
AU2015358016A AU2015358016B2 (en) 2014-12-01 2015-09-08 Method for insar interferometric phase filtering in wavelet domain in conjunction with local frequency estimation
EP15865192.7A EP3229038B1 (en) 2014-12-01 2015-09-08 Wavelet domain insar interferometric phase filtering method in combination with local frequency estimation

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201410717610.3A CN104459633B (zh) 2014-12-01 2014-12-01 结合局部频率估计的小波域InSAR干涉相位滤波方法
CN201410717610.3 2014-12-01

Publications (1)

Publication Number Publication Date
WO2016086699A1 true WO2016086699A1 (zh) 2016-06-09

Family

ID=52905986

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2015/089141 WO2016086699A1 (zh) 2014-12-01 2015-09-08 一种结合局部频率估计的小波域InSAR干涉相位滤波方法

Country Status (4)

Country Link
EP (1) EP3229038B1 (zh)
CN (1) CN104459633B (zh)
AU (1) AU2015358016B2 (zh)
WO (1) WO2016086699A1 (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109697314A (zh) * 2018-12-18 2019-04-30 国网西藏电力有限公司 一种冰崩的最大运动距离的计算方法及其应用
CN111025294A (zh) * 2019-12-12 2020-04-17 南昌大学 基于均方容积卡尔曼滤波器的InSAR相位解缠绕方法
CN113640797A (zh) * 2021-08-09 2021-11-12 北京航空航天大学 一种用于参考条带模式InSAR的前斜视测高方法
CN114355348A (zh) * 2022-01-10 2022-04-15 交通运输部路网监测与应急处置中心 Sar干涉图小波降噪处理方法及其处理装置
CN114859346A (zh) * 2022-04-21 2022-08-05 桂林电子科技大学 一种基于Insarbm3d滤波算法的迭代反演滤波方法
CN115128548A (zh) * 2022-05-27 2022-09-30 西安电子科技大学杭州研究院 一种sar射频干扰检测方法
CN116224332A (zh) * 2023-01-17 2023-06-06 中国矿业大学 一种协调多元指标的雷达干涉相位质量估计方法
CN117213443A (zh) * 2023-11-07 2023-12-12 江苏省地质调查研究院 一种天地深一体化地面沉降监测网建设与更新方法

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104459633B (zh) * 2014-12-01 2016-08-17 中国科学院电子学研究所 结合局部频率估计的小波域InSAR干涉相位滤波方法
CN105976340A (zh) * 2016-05-20 2016-09-28 山东师范大学 改进的基于小波分解的旋滤波算法
CN106097404B (zh) * 2016-05-27 2019-03-12 山东科技大学 利用非线性矢量曲面构建InSAR相位图像模型的方法
CN106680785B (zh) * 2017-03-06 2019-03-12 浙江工业大学 基于小波变换空间变迹的sar图像旁瓣抑制方法
CN108303735A (zh) * 2018-01-30 2018-07-20 单新建 基于最优参数设置的子带干涉测量的地震形变获取方法
CN108416165B (zh) * 2018-03-26 2021-08-27 合肥工业大学 一种适用于弹性成像的高抗噪局部频率估计方法
CN108898155A (zh) * 2018-05-18 2018-11-27 浙江工业大学 一种小波阈值降噪结合卷积神经网络的sar图像目标识别方法
CN109472834B (zh) * 2018-10-23 2023-04-14 桂林电子科技大学 一种基于小波变换的卡尔曼滤波相位展开方法
CN109669182B (zh) * 2018-12-05 2022-05-17 南京邮电大学 无源双基地sar动/静目标联合稀疏成像方法
EP3866105A1 (en) * 2020-02-17 2021-08-18 Paris Sciences et Lettres - Quartier Latin Method for processing insar images to extract ground deformation signals
CN111239736B (zh) * 2020-03-19 2022-02-11 中南大学 基于单基线的地表高程校正方法、装置、设备及存储介质
CN112419205B (zh) * 2020-11-27 2024-04-19 中国矿业大学 一种用于SAR干涉图处理的Goldstein金字塔构建方法
CN113177887A (zh) * 2021-04-16 2021-07-27 中国科学院精密测量科学与技术创新研究院 一种基于多尺度时频分析的快速非局部均值InSAR相位滤波方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6097328A (en) * 1998-07-02 2000-08-01 Raytheon Company Averaging-area-constrained adaptive interferometric filter that optimizes combined coherent and noncoherent averaging
CN103208101A (zh) * 2013-03-28 2013-07-17 中国科学院对地观测与数字地球科学中心 一种基于局部信噪比的干涉图滤波方法
CN103226194A (zh) * 2013-03-26 2013-07-31 中国科学院电子学研究所 一种基于经验模式分解的InSAR干涉相位滤波方法
CN103823219A (zh) * 2014-03-14 2014-05-28 中国科学院电子学研究所 自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法
CN104459633A (zh) * 2014-12-01 2015-03-25 中国科学院电子学研究所 结合局部频率估计的小波域InSAR干涉相位滤波方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5952957A (en) * 1998-05-01 1999-09-14 The United States Of America As Represented By The Secretary Of The Navy Wavelet transform of super-resolutions based on radar and infrared sensor fusion

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6097328A (en) * 1998-07-02 2000-08-01 Raytheon Company Averaging-area-constrained adaptive interferometric filter that optimizes combined coherent and noncoherent averaging
CN103226194A (zh) * 2013-03-26 2013-07-31 中国科学院电子学研究所 一种基于经验模式分解的InSAR干涉相位滤波方法
CN103208101A (zh) * 2013-03-28 2013-07-17 中国科学院对地观测与数字地球科学中心 一种基于局部信噪比的干涉图滤波方法
CN103823219A (zh) * 2014-03-14 2014-05-28 中国科学院电子学研究所 自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法
CN104459633A (zh) * 2014-12-01 2015-03-25 中国科学院电子学研究所 结合局部频率估计的小波域InSAR干涉相位滤波方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
YUE, YINHUAN ET AL.: "SAR Interferogran Filtering Based on Stationary Wavelet Decomposition", HIGH TECHNOLOGY LETTERS, vol. 1, no. 5, 31 May 2002 (2002-05-31), pages 5 - 9, ISSN: 1002-0470 *

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109697314A (zh) * 2018-12-18 2019-04-30 国网西藏电力有限公司 一种冰崩的最大运动距离的计算方法及其应用
CN109697314B (zh) * 2018-12-18 2022-11-01 国网西藏电力有限公司 一种冰崩的最大运动距离的计算方法及其应用
CN111025294A (zh) * 2019-12-12 2020-04-17 南昌大学 基于均方容积卡尔曼滤波器的InSAR相位解缠绕方法
CN113640797A (zh) * 2021-08-09 2021-11-12 北京航空航天大学 一种用于参考条带模式InSAR的前斜视测高方法
CN113640797B (zh) * 2021-08-09 2022-04-12 北京航空航天大学 一种用于参考条带模式InSAR的前斜视测高方法
CN114355348A (zh) * 2022-01-10 2022-04-15 交通运输部路网监测与应急处置中心 Sar干涉图小波降噪处理方法及其处理装置
CN114859346B (zh) * 2022-04-21 2024-05-03 桂林电子科技大学 一种基于Insarbm3d滤波算法的迭代反演滤波方法
CN114859346A (zh) * 2022-04-21 2022-08-05 桂林电子科技大学 一种基于Insarbm3d滤波算法的迭代反演滤波方法
CN115128548A (zh) * 2022-05-27 2022-09-30 西安电子科技大学杭州研究院 一种sar射频干扰检测方法
CN116224332B (zh) * 2023-01-17 2023-10-20 中国矿业大学 一种协调多元指标的雷达干涉相位质量估计方法
CN116224332A (zh) * 2023-01-17 2023-06-06 中国矿业大学 一种协调多元指标的雷达干涉相位质量估计方法
CN117213443A (zh) * 2023-11-07 2023-12-12 江苏省地质调查研究院 一种天地深一体化地面沉降监测网建设与更新方法
CN117213443B (zh) * 2023-11-07 2024-03-19 江苏省地质调查研究院 一种天地深一体化地面沉降监测网建设与更新方法

Also Published As

Publication number Publication date
CN104459633B (zh) 2016-08-17
CN104459633A (zh) 2015-03-25
EP3229038A4 (en) 2018-08-15
EP3229038A1 (en) 2017-10-11
EP3229038B1 (en) 2020-05-27
AU2015358016A1 (en) 2017-07-20
AU2015358016B2 (en) 2018-06-28

Similar Documents

Publication Publication Date Title
WO2016086699A1 (zh) 一种结合局部频率估计的小波域InSAR干涉相位滤波方法
CN109633648B (zh) 一种基于似然估计的多基线相位估计装置及方法
CN104007439B (zh) 一种干涉圆迹sar高程估计处理方法
Natsuaki et al. InSAR local co-registration method assisted by shape-from-shading
CN111856459B (zh) 一种改进的DEM最大似然约束多基线InSAR相位解缠方法
CN103823219B (zh) 自适应迭代的非局部干涉合成孔径雷达干涉相位滤波方法
CN110109100B (zh) 一种基于质量图加权的多基线最小二乘相位解缠方法
CN115060208A (zh) 基于多源卫星融合的输变电线路地质灾害监测方法及系统
CN113340191A (zh) 时间序列干涉sar的形变量测量方法及sar系统
CN104730519A (zh) 一种采用误差迭代补偿的高精度相位解缠方法
CN103824287A (zh) 一种基于稳健平面拟合的相位相关亚像素匹配方法
Sun et al. Improved Goldstein filter for InSAR noise reduction based on local SNR
CN103809180A (zh) 用于InSAR地形测量的方位向预滤波处理方法
CN103226194A (zh) 一种基于经验模式分解的InSAR干涉相位滤波方法
CN103942775B (zh) 基于最大核密度估计的相位相关亚像素匹配方法
Yan et al. A fast non-local means filtering method for interferometric phase based on wavelet packet transform
CN105116410B (zh) 基于线性模型匹配的干涉相位图自适应滤波算法
CN105678716A (zh) 一种地基sar大气干扰相位校正方法及装置
CN104361596A (zh) 一种基于Contourlet变换和Frobenius范数半参考图像质量评价方法
Schmitt et al. Adaptive multilooking of airborne Ka-band multi-baseline InSAR data of urban areas
Natsuaki et al. Performance improvement of InSAR local co-registration method with multiresolution interferogram
Belhadj-Aissa et al. Contextual filtering methods based on the subbands and subspaces decomposition of complex SAR interferograms
CN111696207B (zh) 一种基于引导滤波的多基线dem融合方法
Natsuaki et al. Phase property in complex-correlation and real-imaginary-correlation filtered SAR interferograms and its influence on DEM quality
Wan et al. Phase correlation based local illumination-invariant method for multi-tempro remote sensing image matching

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 15865192

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

REEP Request for entry into the european phase

Ref document number: 2015865192

Country of ref document: EP

ENP Entry into the national phase

Ref document number: 2015358016

Country of ref document: AU

Date of ref document: 20150908

Kind code of ref document: A