CN114594479B - Full scatterer FS-InSAR method and system - Google Patents

Abstract

The invention discloses a method and a system for a full scatterer FS-InSAR, which comprises the following steps: selecting a monitoring area, and combining time sequence SAR images of the monitoring area to generate a differential interference image based on the first interference image; phase unwrapping is carried out on the differential interference pattern to obtain an unwrapped phase; obtaining a time sequence phase corresponding to each SAR imaging moment according to the unwrapping phase; according to the double-scale time domain low-pass filtering, suppressing an atmospheric phase and a noise phase in a time sequence phase to obtain a time sequence low-frequency phase; recalculating the interference image pair combination to obtain a second interference image pair combination, and differentiating the time sequence low-frequency phase according to the second interference image pair combination to obtain a differential low-frequency phase; obtaining an average amplitude according to the time sequence SAR image, and extracting a Full Scatterer (FS) of a non-water body in a monitoring area according to the average amplitude; the earth surface deformation information is obtained based on the full scatterer and the differential low-frequency phase, so that the phase quality is greatly improved, and the earth surface deformation monitoring density is improved.

Description

A full-scatterer FS-InSAR method and system

technical field

The invention relates to the technical field of spaceborne synthetic aperture radar time series interference, in particular to a full scatterer FS-InSAR method and system.

Background technique

Time-series InSAR technology includes PS-type InSAR technology and DS-type InSAR technology, and has been widely used in deformation monitoring fields such as land subsidence, gob collapse, earthquakes, volcanoes, landslides, and the safety of important infrastructure. However, time-series InSAR high-precision surface deformation monitoring still faces two challenges: one is the interference of tropospheric atmospheric effects that are highly uncorrelated in time, especially small-scale atmospheric turbulence effects, whose delayed phase even exceeds the deformation phase. Second, the surface in high vegetation coverage areas is seriously decoherent, resulting in sparse monitoring density and insufficient deformation information.

At present, atmospheric correction methods can be divided into two categories:

The first is to use external atmospheric data modeling to estimate the atmospheric delay phase, such as meteorological numerical models, GPS data, imaging spectrometer data, etc. These methods can correct atmospheric turbulence effects on large and medium scales to a certain extent, but it is difficult to correct small-scale atmospheric turbulence effects due to the interference of weather conditions or spatial and temporal resolution.

The second is atmospheric correction based on the InSAR data itself, such as Stacking InSAR, interferogram overlay atmospheric estimation methods, and PS or DS methods. Both the Stacking InSAR method and the interferogram stacking atmospheric estimation method assume that the surface deformation is linear, so they are not suitable for nonlinear deformation areas. The PS or DS method assumes that the atmospheric phases are spatially correlated within a certain distance (usually less than 2 km), and suppresses the influence of atmospheric effects by the phase difference between adjacent points.

However, when a small-scale atmospheric effect occurs, the atmospheres separated by hundreds of meters or even tens of meters may have large differences, and it is difficult to remove atmospheric effects through phase difference.

The surface decoherence problem is currently mainly solved by distributed target (DS) InSAR technology. There are two kinds of DS extraction methods, one is non-parametric statistical methods, such as KS (Kolmogorov-Smirnov) test, AD (Anderson-Darling) test, BWS (Baumgartner-Weiβ-Schindler) test, etc., which have strong scene adaptability, but computational efficiency Low; the second is parameter statistical methods, such as likelihood ratio test, FaSHPS (Fast Statistically Homogeneous Pixel Selection), HTCI (Hypothesis Test of Confidence Interval) and other methods, although the calculation efficiency is high, but the surface is required to obey a Gaussian distribution. It is difficult for the region to satisfy this assumption. Therefore, the existing DS methods are mainly applied to moderately coherent targets such as bare land, grassland, and farmland, and are less applied to low-coherence areas covered by high-density vegetation.

Therefore, the prior art has the following problems:

1. When using external atmospheric data modeling to estimate the atmospheric delay phase, it is easy to be interfered by factors such as weather conditions or spatial and temporal resolution, and it is difficult to correct the small-scale atmospheric turbulence effect.

2. Atmospheric correction based on InSAR data itself, when using Stacking InSAR and interferogram overlay atmospheric estimation methods, is not suitable for nonlinear deformation areas; PS or DS methods are difficult to suppress the atmosphere by means of phase difference when small-scale atmospheric effects occur. influences.

3. For the problem of surface decoherence, it is easy to solve the problem of low computational efficiency and difficult to meet the assumption conditions for vegetation-covered areas when solved by distributed target InSAR technology.

SUMMARY OF THE INVENTION

In view of the above deficiencies of the prior art, the present application provides a full scatterer FS-InSAR method and system.

In the first aspect, the present application proposes a full scatterer FS-InSAR method, which includes the following steps:

Selecting a monitoring area, and generating a differential interferogram based on the first interferometric image pair combination and combining the time series SAR images of the monitoring area;

performing phase unwrapping on the differential interferogram to obtain an unwrapped phase;

Obtain the time sequence phase corresponding to each SAR imaging moment according to the unwrapping phase;

According to the dual-scale time-domain low-pass filtering, the atmospheric phase and the noise phase in the time-series phase are suppressed to obtain the time-series low-frequency phase;

Recalculate the combination of interference image pairs to obtain a second combination of interference image pairs, and differentiate the time series low-frequency phase according to the second combination of interference image pairs to obtain a differential low-frequency phase;

Obtain an average amplitude according to the time series SAR image, and extract the total scatterers of non-water bodies in the monitoring area according to the average amplitude;

Surface deformation information is obtained based on the total scatterer and the differential low frequency phase.

In some embodiments, in the above-mentioned full-scatterer FS-InSAR method, in the selected monitoring area, a differential interferogram is generated based on the first interferometric image pair combination and in combination with the time-series SAR images of the monitoring area, including:

和第一垂直基线阈值

得到所述第一干涉像对组合，所述第一干涉像对组合为连通网络，所述第一时间基线阈值小于等于两个卫星重访周期。Set the first time baseline threshold

and the first vertical baseline threshold

The first interference image pair combination is obtained, the first interference image pair combination is a connected network, and the first time baseline threshold is less than or equal to two satellite revisit periods.

In some embodiments, in the above-mentioned full scatterer FS-InSAR method, the phase unwrapping of the differential interferogram to obtain the unwrapped phase includes:

Phase unwrapping is performed on the differential interferogram by using a phase unwrapping method to obtain an unwrapped phase;

The phase unwrapping methods include least cost flow method, branch cutting method and region growing method.

In some embodiments, in the above-mentioned full-scatterer FS-InSAR method, obtaining the time-series phase corresponding to each SAR imaging moment according to the unwrapped phase includes:

The least squares method is used to calculate the unwrapped phase at each SAR imaging moment to obtain the corresponding time series phase, which can be expressed as:

（1）

(1)

为形变相位，

为DEM高程误差相位，

为轨道误差相位，

为大气延迟相位，

为噪声相位。in,

is the deformation phase,

is the DEM elevation error phase,

is the orbital error phase,

is the atmospheric delay phase,

is the noise phase.

In some embodiments, in the above-mentioned full scatterer FS-InSAR method, before the low-pass filtering in the dual-scale time domain, the atmospheric phase and the noise phase in the time-series phase are suppressed to obtain the time-series low-frequency phase, including: :

Determine the dual-scale time window;

The dual-scale time window includes a first-scale time window and a second-scale time window, the first-scale time window being smaller than the second-scale time window.

In some embodiments, in the above-mentioned full-scatterer FS-InSAR method, the time-series low-frequency phase is obtained by suppressing the atmospheric phase and the noise phase in the time-series phase according to dual-scale time-domain low-pass filtering, including:

The first scale time window is adopted and the time sequence phase is low-pass filtered by a low-pass filter to obtain the first low-frequency phase:

（2）

(2)

（3）

(3)

为时序相位，

为高斯函数，

为第一尺度时间窗口下输出的

时刻的低通滤波相位，

为N幅SAR影像的成像日期，

为第一尺度时间窗口；in,

is the timing phase,

is a Gaussian function,

is the output under the first scale time window

the low-pass filtered phase at time,

is the imaging date of N SAR images,

is the first scale time window;

The first low-frequency phase is subtracted from the timing phase to obtain a residual phase, where the residual phase is:

（4）

(4)

Using the second scale time window and low-pass filtering the residual phase again through a low-pass filter, the second low-frequency phase is obtained:

（5）

(5)

（6）

(6)

为第二尺度时间窗口，

为第二尺度时间窗口下输出的

时刻的低通滤波相位；in,

is the second scale time window,

is output under the second scale time window

Low-pass filter phase at time;

The first low frequency phase and the second low frequency phase are superimposed to obtain the time series low frequency phase, and the time series low frequency phase is:

（7）

(7)

为所述时序低频相位。in,

is the low frequency phase of the timing sequence.

In some embodiments, in the above-mentioned full scatterer FS-InSAR method, the interferometric image pair combination is recalculated to obtain a second interferometric image pair combination, and the time series low frequency phase is determined according to the second interference image pair combination. Differentiate to get the differential low-frequency phase, including:

和第二垂直基线阈值

得到所述第二干涉像对组合，所述第二干涉像对组合的数目大于所述第一干涉像对组合的数目；Set Second Time Baseline Threshold

and the second vertical baseline threshold

obtaining the second interference image pair combination, where the number of the second interference image pair combination is greater than the number of the first interference image pair combination;

According to the combination of the second interference image pair, the time series low-frequency phase is differentiated to obtain the differential low-frequency phase:

（8）

(8)

为主影像对应的时序低频相位，

为辅影像对应的时序低频相位，

,

。in,

is the time series low frequency phase corresponding to the main image,

is the time-series low-frequency phase corresponding to the auxiliary image,

,

.

In some embodiments, in the above-mentioned full scatterer FS-InSAR method, the obtaining an average amplitude according to the time series SAR image, and extracting the total scatterer of the non-water body in the monitoring area according to the average amplitude, comprising:

Performing relative normalization processing on the time series SAR image to obtain an average amplitude;

A screening threshold is set. If the average amplitude of the time series SAR image is greater than the screening threshold, the corresponding pixel is extracted as a full scatterer, and the full scatterer is expressed as:

（9）

(9)

为时序SAR影像的平均幅度，

为筛选阈值。in,

is the average amplitude of time-series SAR images,

is the filter threshold.

In some embodiments, in the above-mentioned full-scatterer FS-InSAR method, the obtaining surface deformation information based on the full-scatterer and the differential low-frequency phase includes:

extracting the differential low frequency phase corresponding to the full scatterer according to the full scatterer and the differential low frequency phase;

Using the conventional time-series InSAR method, the differential low-frequency phase of the full scatterer is processed to obtain surface deformation information.

In the second aspect, this application proposes a full-scatterer FS-InSAR system, including a differential interferogram generation module, a phase unwrapping module, a time-series phase generation module, a dual-scale time-domain low-pass filtering module, a differential low-frequency phase generation module, a full Scatter extraction module and surface deformation acquisition module;

The differential interferogram generation module is used to select a monitoring area, and generate a differential interferogram based on the first interferometric image pair combination and in combination with the time-series SAR images of the monitoring area;

the phase unwrapping module, configured to perform phase unwrapping on the differential interferogram to obtain an unwrapped phase;

The sequential phase generation module is used to obtain the sequential phase corresponding to each SAR imaging moment according to the unwrapped phase;

The dual-scale time-domain low-pass filtering module is used for suppressing the atmospheric phase and the noise phase in the time-series phase according to the dual-scale time-domain low-pass filtering to obtain the time-series low-frequency phase;

The differential low-frequency phase generation module is configured to recalculate the combination of interference image pairs to obtain a second combination of interference image pairs, and to differentiate the time-series low-frequency phases according to the second combination of interference image pairs to obtain a differential low-frequency phase;

The total scatterer extraction module is configured to obtain the average amplitude according to the time series SAR image, and extract the total scatterer of the non-water body in the monitoring area according to the average amplitude;

The surface deformation acquisition module is configured to obtain surface deformation information based on the total scatterer and the differential low-frequency phase.

Beneficial effects of the present invention:

For those skilled in the art, it is generally recognized that using a low-pass filter to process the residual deformation of a point target (that is, the phase after removing the main deformation and the elevation error phase from the original differential interference phase) requires not only time Low-pass filtering, and spatial low-pass filtering processing is also required, so as to achieve the effect of suppressing the atmospheric phase and suppressing the noise-removing phase;

This scheme can directly suppress the elevation error phase, orbit error phase, atmospheric phase and noise phase from the differential interference phase only by using temporal low-pass filtering processing, and does not require spatial low-pass filtering processing, and the filtered phase quality is greatly improved. In addition, the processing object is all pixels of non-water bodies, no need to select point targets or distributed targets, and it is difficult to correct the small-scale atmospheric turbulence effect in the existing technology; PS or DS method When small-scale atmospheric effects occur, the phase difference method is used to solve the problem. It is difficult to remove the atmospheric influence by means of the distributed target InSAR technology; it is prone to problems of low computational efficiency and difficulty in meeting the assumptions for vegetation-covered areas.

Description of drawings

FIG. 1 is an overall flow chart of the present invention.

FIG. 2 is a spatiotemporal baseline distribution diagram of the first interferometric image pair combination.

FIG. 3 is a spatiotemporal baseline distribution diagram of the second interferometric image pair combination.

Figure 4 is a differential interference phase diagram.

Figure 5 is a differential high-frequency phase diagram after dual-scale time-domain low-pass filtering.

Figure 6 is a differential low-frequency phase diagram after dual-scale time-domain low-pass filtering.

Figure 7 is a graph of the average amplitude.

Figure 8 is a plot of the point target sedimentation rate.

Figure 9 is a plot of the total scatterer sedimentation rate.

FIG. 10 is a schematic block diagram of the system of the present invention.

Detailed ways

The present application proposes a full scatterer FS-InSAR method and system, which solves the problem that it is difficult to correct small-scale atmospheric turbulence effects in the prior art; when small-scale atmospheric effects occur in PS or DS methods, it is difficult to suppress atmospheric effects by means of phase difference ; The problem of low computational efficiency and difficult to meet assumptions for vegetation coverage areas is easy to solve through distributed target InSAR technology. Only through dual-scale time-domain low-pass filtering can greatly improve the phase quality and improve the surface deformation monitoring technology. Effect.

The present application will be further described below with reference to the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solutions of the present invention, and cannot be used to limit the protection scope of the present application.

In the first aspect, the present application proposes a full-scatterer FS-InSAR method, as shown in Figure 1, including the following steps:

S100 selects a monitoring area, and generates a differential interferogram based on the first interferometric image pair combination and in combination with the time-series SAR images of the monitoring area;

Taking the Sentinel-1 images obtained from 20201002 to 20211009 in Tianjin as an example, there are 30 images in total, and the time baselines are shown in the following table:

和第一垂直基线阈值

得到所述第一干涉像对组合，其中，设置的第一时间基线阈值小于等于两个卫星重访周期，本实施例中设置

≤24天，

≤200米，从而得到53个第一干涉像对组合，即小基线干涉像对组合，如图2所示；In order to facilitate the subsequent phase unwrapping, on the premise of ensuring the connectivity of the differential interferogram space-time baseline, for the first interferometric image pair combination, the time baseline and vertical baseline should be reduced as much as possible, so the first time baseline threshold is set.

and the first vertical baseline threshold

The first interference image pair combination is obtained, wherein the set first time baseline threshold is less than or equal to two satellite revisit periods, and in this embodiment, set

≤24 days,

≤200 meters, so as to obtain 53 first interference image pair combinations, that is, small baseline interference image pair combinations, as shown in Figure 2;

就需要大于两个重访周期；值得一提的是，对于不同类型的卫星，其重访周期不同，所以需要根据卫星类型具体确定对应的重访周期。In special cases, when the time interval between two adjacent SAR images is relatively large, the first time baseline threshold may also be greater than the two revisit periods. In this embodiment, when the time interval between two adjacent SAR images is 36 days, then

It needs more than two revisit periods; it is worth mentioning that for different types of satellites, the revisit periods are different, so the corresponding revisit periods need to be specifically determined according to the satellite type.

S200: Perform phase unwrapping on the differential interferogram to obtain an unwrapped phase;

Phase unwrapping is performed on the differential interferogram by using a phase unwrapping method to obtain an unwrapped phase;

It is worth noting that the phase unwrapping method includes but is not limited to the least-cost flow method, the branch cutting method and the region growing method, and the skilled person can choose the phase unwrapping method adaptively according to the actual application.

Wherein, in this embodiment, the phase unwrapping method of minimum cost flow is used to unwrapped the differential interferogram to obtain the unwrapped phase.

S300: Obtain a time sequence phase corresponding to each SAR imaging moment according to the unwrapped phase;

The least squares method is used for the unwrapped phase to calculate the phase at each SAR imaging moment, and the corresponding time sequence phase is obtained, and the time sequence phase can be expressed as:

（1）

(1)

为形变相位，

为DEM高程误差相位，

为轨道误差相位，

为大气延迟相位，

为噪声相位。in,

is the deformation phase,

is the DEM elevation error phase,

is the orbital error phase,

is the atmospheric delay phase,

is the noise phase.

Before obtaining the low-frequency phase of the time series according to the low-pass filtering in the dual-scale time domain, a dual-scale time window needs to be determined. The method for determining the dual-scale time window is as follows:

The dual-scale time window is determined through the dual-scale time-domain low-pass filtering simulation experiment;

Among them, the dual-scale time-domain low-pass filtering simulation experiment is:

Simulation data is used to verify the effectiveness of dual-scale time-domain low-pass filtering to separate the deformation phase:

The time sampling interval is set to 12 days, the deformation occurrence period is 730 days, and the sinusoidal deformation phase is simulated with an amplitude of 10rad and a period of 365 days. Referring to the real SAR time series differential interference phase of multiple non-deformation regions, the Gaussian noise with mean value of 0rad and standard deviation of 8rad is used as a high-frequency signal, which is superimposed with the deformation phase to obtain the simulated time series differential interference phase.

First, a first-scale time window is used for low-pass filtering. Set the time window size from 1 to 365 days, filter the time series phase respectively, and calculate the standard deviation and mean of the difference between the low-frequency phase after each filtering and the simulated deformation phase;

The results show that when the time window is set to 34 days, the standard deviation of the error between the low-frequency phase and the simulated deformation phase is at least 1.85rad, and the mean is -1.34rad.

Based on the filtering results of 34 days, subtract the first low-frequency phase from the simulated deformation phase to obtain the residual deformation phase; subtract the first low-frequency phase from the simulated time series differential phase to obtain the residual differential phase, and still set the time window size to 1~ For 365 days, filter the residual differential phase respectively, and calculate the standard deviation and mean of the difference between the secondary low-frequency phase and the residual deformation phase;

The results show that when the time window is 77 days, the error standard deviation of the secondary filtering is at least 3.85rad, and the mean is -1.35rad.

The time series low frequency phase of the two filters is superimposed to obtain the total time series low frequency phase. The standard deviation of the difference between the total time series low frequency phase and the simulated deformation phase is 1.74rad, and the mean value is -1.32rad. Compared with the first filtering, the standard deviation of the difference is reduced by 5.43%, indicating that the dual-scale temporal low-pass filtering can better fit the real deformation than the single filtering.

The above simulation experiments show that in the dual-scale time-domain low-pass filtering, the time window of the first filtering should be small, which can be set to 30-40 days; the time window of the second filtering should be large, which can be set to 60-40 days. 80 days.

Therefore, the range of the time window of the first scale is 30-40 days, and the range of the time window of the second scale is 60-80 days.

S400: According to dual-scale time-domain low-pass filtering, suppress the atmospheric phase and the noise phase in the time-series phase to obtain the time-series low-frequency phase;

表现为低频信号；轨道误差相位、大气延迟相位、噪声相位均为高频信号；高程误差相位与垂直基线有关，表现为高频信号；因此，根据不同相位成分的时域频谱特征，利用低通滤波器，可将形变相位与其他相位分离，从而有效抑制大气延迟相位和噪声相位的影响；In the time series, the deformation phase

It appears as a low-frequency signal; the orbital error phase, atmospheric delay phase, and noise phase are all high-frequency signals; the elevation error phase is related to the vertical baseline and appears as a high-frequency signal; therefore, according to the time-domain spectral characteristics of different phase components, the low-pass The filter can separate the deformation phase from other phases, thereby effectively suppressing the influence of atmospheric delay phase and noise phase;

The low-pass filter may be a Gaussian filter or other low-pass filters. This embodiment takes the Gaussian filter as an example for detailed description:

The first scale time window is adopted and the time sequence phase is low-pass filtered by a low-pass filter to obtain the first low-frequency phase:

（2）

(2)

（3）

(3)

为时序相位，

为高斯函数，

为第一尺度时间窗口下输出的

时刻的低通滤波相位，

为N幅SAR影像的成像日期，

为第一尺度时间窗口；in,

is the timing phase,

is a Gaussian function,

is the output under the first scale time window

the low-pass filtered phase at time,

is the imaging date of N SAR images,

is the first scale time window;

Since the Gaussian filter is implemented in a weighted average manner, the closer the center point is to the greater the phase weight, and the farther away the center point is, the smaller the phase weight is. Therefore, although the primary low-frequency phase after filtering is dominated by deformation phase, it also contains a small amount of elevation error phase, orbit error phase, atmospheric phase and noise phase;

Therefore, the first low-frequency phase is subtracted from the timing phase to obtain a residual phase, and the residual phase is:

（4）

(4)

The residual phase includes not only high frequency phases such as elevation error phase, orbit error phase, atmospheric phase and noise phase, but also a small amount of deformation phase. In order to extract this part of the deformed phase, a second scale time window is used and the residual phase is again low-pass filtered by a low-pass filter to obtain the second low-frequency phase:

（5）

(5)

（6）

(6)

为第二尺度时间窗口，

为第二尺度时间窗口下输出的

时刻的低通滤波相位；in,

is the second scale time window,

is the output under the second scale time window

Low-pass filter phase at time;

The first low frequency phase and the second low frequency phase are superimposed to obtain the time series low frequency phase, and the time series low frequency phase is:

（7）

(7)

为所述时序低频相位。in,

is the low frequency phase of the timing sequence.

S500: Recalculate the interference image pair combination to obtain a second interference image pair combination, and differentiate the time series low frequency phase according to the second interference image pair combination to obtain a differential low frequency phase;

和第二垂直基线阈值

得到所述第二干涉像对组合，所述第二干涉像对组合的数目大于所述第一干涉像对组合的数目；其中，

≤100天，

≤200米，从而得到194个第二干涉像对组合，如图3所示；Take the Sentinel-1 images obtained from 20201002 to 20211009 in Tianjin as an example, there are 30 images in total, and the time baseline table will not be repeated, and the second time baseline threshold is set.

and the second vertical baseline threshold

The second interference image pair combination is obtained, and the number of the second interference image pair combination is greater than the number of the first interference image pair combination; wherein,

≤100 days,

≤200 meters, thus obtaining 194 second interference image pair combinations, as shown in Figure 3;

According to the combination of the second interference image pair, the time series low-frequency phase is differentiated to obtain the differential low-frequency phase:

（8）

(8)

为主影像对应的时序低频相位，

为辅影像对应的时序低频相位，

,

。in,

is the time series low frequency phase corresponding to the main image,

is the time-series low-frequency phase corresponding to the auxiliary image,

,

.

S600: Obtain an average amplitude according to the time series SAR image, and extract the total scatterers of non-water bodies in the monitoring area according to the average amplitude;

Performing relative normalization processing on the time series SAR images (ie, dividing the amplitude of each SAR image by its mean amplitude) to obtain an average amplitude;

Set a screening threshold. If the average amplitude of the time series SAR image is greater than the screening threshold, the corresponding pixel is extracted as a full scatterer. The screening threshold is set so that the full scatterer does not contain water pixels, and the full scatterer Expressed as:

（9）

(9)

为时序SAR影像的平均幅度，

为筛选阈值，筛选阈值可设为0.5，由于不同计算方法会得到不同的平均幅度，因此筛选阈值可根据具体情况而定。in,

is the average amplitude of time-series SAR images,

For the screening threshold, the screening threshold can be set to 0.5. Since different calculation methods will obtain different average amplitudes, the screening threshold can be determined according to the specific situation.

It should be noted that since the low-frequency phase of the time series is dominated by the deformation phase, and other high-frequency phases have been basically removed, we can choose non-water pixels as the full scatterers (FS) for surface deformation monitoring, which is similar to PS Compared with the -InSAR or DS-InSAR method, this will greatly increase the monitoring density and obtain global surface deformation information.

Optionally, due to the huge number of total scattering volume pixels, when the selected monitoring area is large, the calculation time is long, and conventional stable point target selection methods can be used, such as amplitude dispersion method, amplitude difference separation method, The average coherence coefficient method is used to extract the required monitoring objects.

S700: Obtain surface deformation information based on the total scatterer and the differential low-frequency phase.

extracting the differential low frequency phase corresponding to the full scatterer according to the full scatterer and the differential low frequency phase;

Using the conventional time-series InSAR method, the differential low-frequency phase of the full scatterer is processed to obtain surface deformation information.

Among them, the specific processing process includes full scatterer target triangulation network connection, solving the relative linear deformation rate and relative height difference error of adjacent targets, solving the linear deformation rate and elevation error of each target, residual phase calculation, nonlinear deformation estimation, Accumulated deformation information acquisition, etc.

The following describes the core effects of this scheme with reference to Figures 4-9:

Among them, Figure 4 is the differential interference phase of the Sentinel-1 image in Tianjin, Figure 5 is the differential high-frequency phase after dual-scale temporal low-pass filtering of the Sentinel-1 image in Tianjin, and Figure 6 is the Sentinel-1 image in Tianjin. The differential low-frequency phase after dual-scale time-domain low-pass filtering; it can be seen that the atmospheric phase and the noise phase are effectively suppressed only by the time low-pass filtering, and the phase quality is greatly improved.

Among them, Figure 7 is the average amplitude in the range of 400*400, Figure 8 is the sedimentation rate of the point target, and Figure 9 is the sedimentation rate of the full scatterer. Comparison of results: the number of point targets is 23782, and the monitoring density is: 14.86%; The number is 157,685, and the monitoring density is: 98.55%, which is 6.63 times that of the point target. The comparison shows that the monitoring object is a non-water body Full Scatterers (FS). Compared with the existing point target or distributed target method, monitoring Density has been breakthroughly improved.

In the second aspect, the present application proposes a full-scatterer FS-InSAR system, as shown in Figure 10, including a differential interferogram generation module, a phase unwrapping module, a sequential phase generation module, a dual-scale time-domain low-pass filtering module, a differential interferogram generation module, and a differential interferogram generation module. Low frequency phase generation module, full scatterer extraction module and surface deformation acquisition module;

The differential interferogram generation module is used to select a monitoring area, and generate a differential interferogram based on the first interferometric image pair combination and in combination with the time-series SAR images of the monitoring area;

the phase unwrapping module, configured to perform phase unwrapping on the differential interferogram to obtain an unwrapped phase;

The sequential phase generation module is used to obtain the sequential phase corresponding to each SAR imaging moment according to the unwrapped phase;

The dual-scale time-domain low-pass filtering module is used for suppressing the atmospheric phase and the noise phase in the time-series phase according to the dual-scale time-domain low-pass filtering to obtain the time-series low-frequency phase;

The differential low-frequency phase generation module is configured to recalculate the combination of interference image pairs to obtain a second combination of interference image pairs, and differentiate the time-series low-frequency phases according to the second combination of interference image pairs to obtain a differential low-frequency phase;

The total scatterer extraction module is configured to obtain the average amplitude according to the time series SAR image, and extract the total scatterer of the non-water body in the monitoring area according to the average amplitude;

The surface deformation acquisition module is configured to obtain surface deformation information based on the total scatterer and the differential low-frequency phase.

The above are only the preferred embodiments of the present invention. It should be pointed out that the technical solutions made by those skilled in the art without departing from the technical solutions should be regarded as falling within the requirements of the claims. scope of protection.

Claims (10)

1. a full scatterer FS-InSAR method, is characterized in that: comprise the following steps: Selecting a monitoring area, and generating a differential interferogram based on the first interferometric image pair combination and combining the time series SAR images of the monitoring area; performing phase unwrapping on the differential interferogram to obtain an unwrapped phase; Obtain the time sequence phase corresponding to each SAR imaging moment according to the unwrapping phase; According to the dual-scale time-domain low-pass filtering, the atmospheric phase and the noise phase in the time-series phase are suppressed to obtain the time-series low-frequency phase; Recalculate the combination of interference image pairs to obtain a second combination of interference image pairs, and differentiate the time series low-frequency phase according to the second combination of interference image pairs to obtain a differential low-frequency phase; Obtain an average amplitude according to the time series SAR image, and extract the total scatterers of non-water bodies in the monitoring area according to the average amplitude; Surface deformation information is obtained based on the total scatterer and the differential low frequency phase. 2 . The full scatterer FS-InSAR method according to claim 1 , wherein, in the selected monitoring area, differential interferometry is generated based on the combination of the first interferometric image pair and combined with the time-series SAR images of the monitoring area. 3 . Figures, including: Set the first time baseline threshold

and the first vertical baseline threshold

The first interference image pair combination is obtained, the first interference image pair combination is a connected network, and the first time baseline threshold is less than or equal to two satellite revisit periods. 3. A full scatterer FS-InSAR method according to claim 2, wherein: performing phase unwrapping on the differential interferogram to obtain an unwrapped phase, comprising: Phase unwrapping is performed on the differential interferogram by using a phase unwrapping method to obtain an unwrapped phase; The phase unwrapping methods include least cost flow method, branch cutting method and region growing method. 4. A full scatterer FS-InSAR method according to claim 3, wherein: obtaining the time sequence phase corresponding to each SAR imaging moment according to the unwrapping phase, comprising: The least squares method is used to calculate the unwrapped phase at each SAR imaging moment to obtain the corresponding time series phase, which can be expressed as:

(1) in,

is the deformation phase,

is the DEM elevation error phase,

is the orbital error phase,

is the atmospheric delay phase,

is the noise phase. 5. A full scatterer FS-InSAR method according to claim 4, characterized in that: in the two-scale time domain low-pass filtering, the atmospheric phase and the noise phase in the time sequence phase are suppressed to obtain the time sequence Before the low frequency phase, including: Determine the dual-scale time window; The dual-scale time window includes a first-scale time window and a second-scale time window, the first-scale time window being smaller than the second-scale time window. 6 . The full scatterer FS-InSAR method according to claim 5 , wherein: according to the dual-scale time-domain low-pass filtering, the atmospheric phase and the noise phase in the time-series phase are suppressed to obtain the time-series low-frequency phase. 7 . Phase, including: Using the first scale time window and performing low-pass filtering on the timing phase through a low-pass filter, the first low-frequency phase is obtained:

(2)

(3) in,

is the timing phase,

is a Gaussian function,

is the output under the first scale time window

the low-pass filtered phase at time,

is the imaging date of N SAR images,

is the first scale time window; The first low-frequency phase is subtracted from the timing phase to obtain a residual phase, where the residual phase is:

(4) Using the second scale time window and low-pass filtering the residual phase again through a low-pass filter, the second low-frequency phase is obtained:

(5)

(6) in,

is the second scale time window,

is the output under the second scale time window

Low-pass filter phase at time; The first low frequency phase and the second low frequency phase are superimposed to obtain the time series low frequency phase, and the time series low frequency phase is:

(7) in,

is the low frequency phase of the timing sequence. 7 . The full scatterer FS-InSAR method according to claim 6 , wherein the recalculation of the interference image pair combination is performed to obtain a second interference image pair combination, and the second interference image pair combination is obtained according to the second interference image pair combination. 8 . The time series low-frequency phase is differentiated to obtain a differential low-frequency phase, including: Set Second Time Baseline Threshold

and the second vertical baseline threshold

obtaining the second interference image pair combination, where the number of the second interference image pair combination is greater than the number of the first interference image pair combination; According to the combination of the second interference image pair, the time series low-frequency phase is differentiated to obtain the differential low-frequency phase:

(8) in,

is the time series low frequency phase corresponding to the main image,

is the time-series low-frequency phase corresponding to the auxiliary image,

,

. 8 . The full scatterer FS-InSAR method according to claim 7 , wherein the average amplitude is obtained according to the time series SAR image, and the total scatterer of the non-water body in the monitoring area is extracted according to the average amplitude. 9 . ,include Performing relative normalization processing on the time series SAR image to obtain an average amplitude; A screening threshold is set. If the average amplitude of the time series SAR image is greater than the screening threshold, the corresponding pixel is extracted as a full scatterer, and the full scatterer is expressed as:

(9) in,

is the average amplitude of time-series SAR images,

is the filter threshold. 9 . The full scatterer FS-InSAR method according to claim 8 , wherein the obtaining surface deformation information based on the full scatterer and the differential low-frequency phase comprises: 10 . extracting the differential low frequency phase corresponding to the full scatterer according to the full scatterer and the differential low frequency phase; Using the conventional time-series InSAR method, the differential low-frequency phase of the full scatterer is processed to obtain surface deformation information. 10. A full-scatterer FS-InSAR system, characterized in that it includes a differential interferogram generation module, a phase unwrapping module, a time-series phase generation module, a dual-scale time-domain low-pass filtering module, a differential low-frequency phase generation module, and a full-scattering module. Volume extraction module and surface deformation acquisition module; The differential interferogram generation module is used for selecting a monitoring area, and generating a differential interferogram based on the first interferometric image pair combination and in combination with the time-series SAR images of the monitoring area; the phase unwrapping module, configured to perform phase unwrapping on the differential interferogram to obtain an unwrapped phase; The sequential phase generation module is used to obtain the sequential phase corresponding to each SAR imaging moment according to the unwrapped phase; The dual-scale time-domain low-pass filtering module is used for suppressing the atmospheric phase and the noise phase in the time-series phase according to the dual-scale time-domain low-pass filtering to obtain the time-series low-frequency phase; The differential low-frequency phase generation module is configured to recalculate the combination of interference image pairs to obtain a second combination of interference image pairs, and differentiate the time-series low-frequency phases according to the second combination of interference image pairs to obtain a differential low-frequency phase; the total scatterer extraction module, configured to obtain the average amplitude according to the time series SAR image, and extract the total scatterer of the non-water body in the monitoring area according to the average amplitude; The surface deformation acquisition module is configured to obtain surface deformation information based on the total scatterer and the differential low-frequency phase.

Images