CN111998766A - Surface deformation inversion method based on time sequence InSAR technology - Google Patents

Abstract

The invention relates to a surface deformation inversion method based on time-series InSAR technology, comprising the following steps: geocoding and image registration, generating differential interference pairs, selecting PS points, space-time unwrapping, and error separation. Specifically, in the image registration method, coarse registration and fine registration are used to ensure the registration accuracy required in the TOPS imaging mode; for the atmospheric delay error, different algorithms are used to process the layered atmosphere and the floating atmosphere, respectively. The influence of atmospheric errors is weakened to the greatest extent, and the accuracy of deformation inversion is improved; if GNSS data is provided, the atmospheric delay can be corrected by using the observation data, and a joint adjustment is performed as a constraint condition in the InSAR observation equation to ensure the deformation results. Accuracy and reliability of all-weather, large-area, continuous, high-spatial-resolution surface deformation monitoring.

Description

A Surface Deformation Inversion Method Based on Time Series InSAR Technology

technical field

The invention relates to the field of surface deformation monitoring, in particular to a surface deformation inversion method based on time series InSAR technology.

Background technique

In recent years, geological disasters caused by surface deformation have received extensive attention. The causes include natural factors such as crustal movement, lithology, folds, and fault structures, as well as human factors such as excessive exploitation of groundwater and minerals. Serious surface deformation will cause serious social and economic losses, so it is of great significance to carry out deformation monitoring and inversion.

Traditional monitoring methods of surface deformation, especially seismic coseismic deformation field, urban land subsidence and mine surface subsidence, mostly use first- and second-class precision level or GNSS monitoring. Precision leveling is suitable for engineering surveying with a small area, and has been widely used in early urban subsidence monitoring. With the development of technology, GNSS technology with high accuracy, rapid network deployment and short period has also been successfully applied in deformation monitoring. middle. However, regardless of leveling technology or GNSS technology, traditional monitoring methods still have obvious defects, mainly in (1) observation points are easily damaged, resulting in the lack of subsequent observation results; Scope monitoring requirements. With point-line monitoring, interpolation is usually required to estimate the overall deformation trend. Therefore, the distribution and density of monitoring points greatly affect the accuracy and reliability of monitoring results. If the monitoring points are too sparse, the deformation surface obtained by interpolation cannot be effective. It reflects the overall trend; (3) traditional observation methods are labor-intensive, time-consuming, and inefficient; (4) for complex terrain areas that cannot be reached by manpower, there are fewer monitoring points available. Therefore, how to provide a surface deformation inversion method is of great significance in making up for the insufficiency of traditional monitoring methods and monitoring surface deformation in a large-scale and high-efficiency manner.

Synthetic Aperture Radar Interferometry (InSAR) is an emerging space geodetic technology. Compared with traditional surface deformation monitoring methods, it has made great progress in monitoring range and monitoring efficiency. Combining high-resolution imaging technology, synthetic aperture radar technology and interferometry technology, InSAR has been widely used in civil and military monitoring fields such as earthquakes, volcanoes, landslides, and soil erosion. This technology can solve the problems of low spatial density and long operation cycle of monitoring networks faced by traditional surveys such as leveling and GNSS monitoring.

Due to the traditional differential interferometry (DInSAR) monitoring the deformation information between the two imaging moments, the observation of the two images in a long period of time will be restricted by factors such as different atmospheric environments, spatial baselines, and scatterer decoherence, which limit the detection of the two images. The measurement accuracy of this technique. Spatiotemporal loss of correlation means that the additional noises of the two imaging are different or irrelevant, and the phase difference cannot be canceled, resulting in a low signal-to-noise ratio and inconspicuous interferograms, which affect the final deformation inversion accuracy; atmospheric delay is reflected in the variable atmospheric conditions. will result in different phase delays that may obscure deformation or topographic information. Therefore, how to overcome the above-mentioned temporal and spatial de-correlation, weaken the influence of atmospheric delay on the extraction of deformation information, and maximize the use of SAR image data to establish differential interferograms of long-term series are problems that need to be solved.

Time-series InSAR technologies such as PSInSAR (Permanent Scatter Interferometry) and SBAS (Small BaselineSubset) reduce the impact of spatiotemporal loss of correlation by selecting high coherence points and using short baselines. However, with the development of spaceborne SAR platforms, different imaging methods also put forward higher requirements for image registration, and ordinary polynomial registration methods cannot meet the accuracy requirements. At the same time, for the processing of atmospheric errors, the existing filtering processing methods do not consider the different characteristics of the stratified atmosphere and the floating atmosphere. In addition, considering the complementary characteristics of GNSS and InSAR data in terms of temporal and spatial resolution, how to integrate these two technologies to obtain high-precision deformation sequences is also a problem worth solving.

SUMMARY OF THE INVENTION

In view of this, in order to overcome the defects existing in the prior art, the present application provides the following technical solutions.

A surface deformation inversion method based on time series InSAR technology, comprising the following steps:

S1. Geocoding and image registration, using an external DEM for geometric registration, and performing related registration based on the amplitude and intensity information of the image to enhance the precise registration of spectral diversity;

S2. Generate a differential interference pair;

S3, select PS points, and use image intensity and phase information to iteratively select PS points;

S4. Space-time unwrapping, adding the time dimension on the basis of the two-dimensional plane, forming a three-dimensional space search path, and determining the residual points for three-dimensional unwrapping;

S5. Error separation. For the layered atmosphere and the floating atmosphere, the methods of parameterization and space-time filtering are used to separate the atmospheric errors.

Further, if GNSS data is provided, the above-mentioned surface deformation inversion method also includes GNSS-InSAR fusion, which is mainly reflected in atmospheric delay correction and InSAR data adjustment with GNSS constraints.

Preferably, in the step S1, an external DEM is used for geometric registration, specifically:

Convert the external DEM to the image coordinate system of the main image, establish a look-up table, and then calculate the position of a pixel through the satellite orbit parameters, set an appropriate window size in the main image, uniformly select pixels, and use the distance equation, Doppler Simultaneously solve the coordinates of the pixel points through the Le equation and the earth model equation, and inversely calculate the coordinates of these pixels in the secondary image coordinate system;

Use Delaunay triangulation interpolation to calculate the offset of all pixels, and use the final offset file to correct the secondary image.

Preferably, in the step S1, correlation registration is performed based on the amplitude and intensity information of the image, specifically:

Select an appropriate window in the main image, and calculate the coherence coefficient of the window;

Use maximum likelihood estimation to perform offset between main and auxiliary images;

The transformation matrix is obtained according to the calculated offset, and the registration of the secondary image is performed, and at the same time, the registration accuracy of 0.05 pixels is ensured.

Preferably, the precise registration of enhanced spectral diversity in the step S1 is specifically:

Extract the overlapping area between the main image and the resampled secondary image burst;

The main and auxiliary images interfere to form the front and back views of the overlapping area, which are interfered and filtered again;

Estimation of ESD phase and registration errors;

The azimuth offset is calculated by fitting a polynomial to all overlapping area deviations.

Preferably, the step S2 generating a differential interference pair includes:

The time baseline, spatial baseline and frequency baseline formed by the main image and the auxiliary image are important factors that affect the interference effect. In order to select the best public main image, it is necessary to maximize the sum of the correlation coefficients of the time series interferogram:

In the formula, f(γ _t ), f(γ _B ), f(γ _DC ) are the time baseline coefficient, vertical baseline coefficient and Doppler centroid frequency difference coefficient, respectively. After finding the best main image, an interference pair is formed;

The external DEM is transferred to the radar coordinate system to form the terrain phase, and the interference phase is subtracted from the terrain phase and the ground phase to form a differential interference pair.

Preferably, the step S3 selecting the PS point includes:

Use the amplitude dispersion coefficient for rough selection, and keep the point targets higher than the threshold;

The coherence coefficient of the phase point in time is calculated by the differential phase, the filtering phase and the topographic residual phase after the differential interference. By selecting an appropriate coherence coefficient threshold, the points higher than the threshold are identified as PS target points. This step The final selection of PS target points needs to be achieved through filtering and iterative operations.

Preferably, the space-time disentanglement in the step S4 includes:

After removing the terrain phase, the single-point phase composition in the interference pair is expressed as:

W{φ}=W{φ _d +φ _τ +φ _a +φ _o +φ _n }

where W{·} is the phase winding, φ is the single-point interference phase, and φ _d , φ _τ , φ _a , φ _o , and φ _n represent the deformation, terrain error, atmosphere, orbit, and noise phase, respectively.

后，将其去除得到：Among them, φ _d , φ _a and φ _o have spatial correlation, which can be separated by spatial low-pass filtering, and the spatial correlation phase can be estimated

Then, remove it to get:

In the formula, the superscript u represents spatial uncorrelation, and δ represents the total amount of unremoved spatially uncorrelated deformation, atmosphere, and orbital errors; because these three have strong spatial correlation, the value of δ is very small;

Construct the coherence coefficient:

Since δ is small, the coherence coefficient γ can reflect the influence of noise on a single point. The larger the noise mean, the larger the coherence coefficient value. The maximum value of γ can be calculated by using the spatial search method, and the corresponding

误差后，进行时空三维解缠。in removing

After the error, the space-time 3D unwrapping is performed.

Preferably, the error separation in step S5 includes:

The unwrapped phase value is:

为残余具有空间相关性的地形残差相位，k为单差模糊度整数；如果三维解缠正确，则大部分点的k应该相等。为分离出形变相位，需要在弧段上进行时空域滤波。where φ _uw is the unwrapped phase,

is the residual terrain residual phase with spatial correlation, k is a single-difference ambiguity integer; if the 3D unwrapping is correct, k should be equal for most points. In order to separate out the deformation phase, it is necessary to perform spatiotemporal filtering on the arc segment.

代表两点间做差，然后构建三角网进行时间域低通滤波，将其分离后得到In the formula,

represents the difference between two points, and then constructs a triangular network for low-pass filtering in the time domain, and separates them to obtain

和

分别表示副影像大气和轨道分量。

在时间域上表现为较强的相关性，利用时间域低通可以获取形变和主影像大气；其余量在时间域上表现为不相关，要获取副影像大气，

在空间表现为相关性，

表现为随机噪声误差，因此可以通过空间的高通滤波分离出副影像的大气延迟相位和轨道误差；最后通过从弧段到点上的积分可以获取每个点的φ_d及

In the formula, [ ] _{LP_time} represents the time domain low-pass filtering,

and

represent the atmospheric and orbital components of the secondary image, respectively.

There is a strong correlation in the time domain, and the deformation and the atmosphere of the main image can be obtained by using the low-pass in the time domain; the remaining quantities are irrelevant in the time domain. To obtain the atmosphere of the auxiliary image,

Shows correlation in space,

It is a random noise error, so the atmospheric delay phase and orbital error of the secondary image can be separated by high-pass filtering in space; finally, the φ _d and

Preferably, if there is still a relatively obvious long-wave trend in the differential phase after the high and low pass filtering in the error separation, that is, the stratified atmospheric delay interference, it is absorbed and removed by using a linear function.

Preferably, the surface deformation inversion method further includes GNSS-InSAR fusion. GNSS has the characteristics of all-weather, high precision, high efficiency, etc. It can continuously observe for a long time and provide high-precision observation data. Due to the long time period of SAR images, it is difficult to provide information with sufficient time resolution. Using the GNSS-InSAR fusion technology, it is possible to correct the errors (such as atmospheric delay, etc.) that are difficult to eliminate in the InSAR data itself, while achieving high GNSS time resolution. , High position accuracy and InSAR high spatial resolution are effectively unified, so as to obtain high reliability and high precision surface deformation results. The technical fusion of GNSS and InSAR is mainly reflected in: GNSS applied to atmospheric delay correction; InSAR data adjustment with GNSS constraints.

Specifically, after obtaining the unwrapped phase, if there is external GNSS data, obtain the ZWD from the GNSS data, invert the atmospheric delay, and remove it from the unwrapped phase;

The DEM error phase is added back and re-parameterized. The least squares method is used to calculate the deformation phase. The GNSS observation equation is used as a constraint, and the joint adjustment is performed with the InSAR observation equation.

Beneficial technical effect obtained by the present invention:

1) The present invention solves the problems that the traditional leveling and GNSS deformation monitoring methods have long monitoring periods, sparse spatial points, manpower and material resources, and the traditional DInSAR monitoring methods are seriously affected by space-time decoherence and atmospheric delay. The invention has strong advantages in the registration method, the atmospheric delay processing algorithm and the GNSS-InSAR fusion technology respectively, and can monitor the surface deformation in an all-weather, large-area, continuous and high-spatial resolution manner, reduce the influence of errors, and improve the monitoring efficiency;

2) In the image registration method, the present invention adopts the method of coarse registration and fine registration. On the basis of rough registration, the registration method of enhanced spectral decomposition is used to ensure the registration accuracy in TOPS imaging mode;

3) For the atmospheric delay error, the present invention uses different algorithms to process the layered atmosphere and the floating atmosphere respectively, so as to weaken the influence of the atmospheric error to the greatest extent and improve the accuracy of deformation inversion;

4) If the GNSS data is provided, the atmospheric delay can be corrected by using the GNSS data, and the joint adjustment can be performed as a constraint condition in the InSAR observation equation to ensure the accuracy and reliability of the deformation results.

The above description is only an overview of the technical solution of the present application, in order to be able to understand the technical means of the present application more clearly, so that it can be implemented according to the content of the description, and in order to make the above and other purposes, features and advantages of the present application more clearly understandable , the preferred embodiments of the present application and the accompanying drawings are described in detail below.

The above and other objects, advantages and features of the present application will be more apparent to those skilled in the art from the following detailed description of the specific embodiments of the present application in conjunction with the accompanying drawings.

Description of drawings

In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the following briefly introduces the accompanying drawings that need to be used in the description of the embodiments or the prior art. Obviously, the drawings in the following description are For some embodiments of the present application, for those of ordinary skill in the art, other drawings can also be obtained based on these drawings without any creative effort. Similar elements or parts are generally identified by similar reference numerals throughout the drawings. In the drawings, each element or section is not necessarily drawn to actual scale.

1 is a flowchart of a method for inversion of surface deformation based on time-series InSAR technology in an embodiment of the present disclosure;

FIG. 2 is a flowchart of image registration in an embodiment of the present disclosure.

Detailed ways

In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments It is a part of the embodiments of the present application, but not all of the embodiments.

Furthermore, this application may repeat reference numerals and/or letters in different instances. This repetition is for the purpose of simplicity and clarity and does not in itself indicate a relationship between the various embodiments and/or arrangements discussed.

It should also be noted that in this document, relational terms such as first and second are used only to distinguish one entity or operation from another, and do not necessarily require or imply those entities or operations There is no such actual relationship or order between them. Moreover, the terms "comprising", "comprising" or any other variation thereof are intended to encompass non-exclusive inclusion.

Example 1

As shown in accompanying drawing 1, a kind of surface deformation inversion method based on time series InSAR technology, comprises the following steps:

S1. Geocoding and image registration, using an external DEM for geometric registration, performing related registration based on the amplitude and intensity information of the image, and using enhanced spectral diversity for precise registration, as shown in FIG. 2 .

For the TOPS scanning method in Sentinel-1A IW mode, when scanning each burst, the antenna beam rotates and scans from back to front along the azimuth direction, which compresses the azimuth beam pattern and suppresses the scallop effect caused by imaging. Due to the special nature of the TOPS mode, the Doppler centroid shift is required, and the registration accuracy of the data in the azimuth direction is required to be 1/1000. Therefore, the combined method of coarse registration and fine registration is used to process images. For coarse registration, the geometric registration method of SAR imaging geometry, orbit data and external DEM is used to obtain the offset between the main and auxiliary images. In order to avoid the phase winding problem in the process of fine registration using ESD, The coarse registration accuracy should be as close as possible to 0.05 pixels.

Fine registration is mainly to correct the systematic errors caused by the non-parallel orbits and the azimuth particularity of the TOPS mode. Using the enhanced spectral diversity method, this method mainly processes one-dimensional signals in the frequency domain, and uses the response impulse function to perform interference processing on the overlapping areas of the two adjacent bursts in the azimuth of the main and auxiliary images, and then the result of the interference processing. The interference processing is performed again to obtain the interference phase difference in the overlapping area, and finally the registration error in the azimuth direction is obtained.

The geometric registration is performed using an external DEM, specifically:

Convert the external DEM to the image coordinate system of the main image, establish a look-up table, and then calculate the position of a pixel through the satellite orbit parameters, set an appropriate window size in the main image, uniformly select pixels, and use the distance equation, Doppler Simultaneously solve the coordinates of the pixel points through the Le equation and the earth model equation, and inversely calculate the coordinates of these pixels in the secondary image coordinate system;

Use Delaunay triangulation interpolation to calculate the offset of all pixels, and use the final offset file to correct the secondary image.

The correlation registration is performed based on the amplitude and intensity information of the image, specifically:

Select an appropriate window in the main image, and calculate the coherence coefficient of the window;

Use maximum likelihood estimation to perform offset between main and auxiliary images;

The transformation matrix is obtained according to the calculated offset, and the registration of the secondary image is performed, and at the same time, the registration accuracy of 0.05 pixels is ensured.

The fine registration of the enhanced spectral diversity, specifically:

Extract the overlapping area between the main image and the resampled secondary image burst;

The main and auxiliary images interfere to form the front and back views of the overlapping area, which are interfered and filtered again;

Estimation of ESD phase and registration errors;

The azimuth offset is calculated by fitting a polynomial to all overlapping area deviations.

S2. Generate a differential interference pair.

Specifically, the time baseline, spatial baseline and frequency baseline formed by the main image and the auxiliary image are important factors that affect the interference effect. In order to select the best public main image, it is necessary to maximize the sum of the correlation coefficients of the time series interferogram:

In the formula, γ is the coherence coefficient, f(γ _t ), f(γ _B ), f(γ _DC ) are the time baseline coefficient, vertical baseline coefficient and Doppler centroid frequency difference coefficient, respectively. interference pair.

The external DEM is transferred to the radar coordinate system to form the terrain phase, and the interference phase is subtracted from the terrain phase and the ground phase to form a differential interference pair.

S3, select PS points, and use image intensity and phase information to iteratively select PS points.

Rough selection is performed using the amplitude dispersion coefficient, and point targets above the threshold are retained.

The coherence coefficient of the phase point in time is calculated by the differential phase, the filtering phase and the topographic residual phase after the differential interference. By selecting an appropriate coherence coefficient threshold, the points higher than the threshold are identified as PS target points. This step The final selection of PS target points needs to be achieved through filtering and iterative operations. Specifically, for the selected point of amplitude dispersion, when the pixel has a main scatterer, the phase standard deviation and amplitude have the following relationship:

是相位标准差；σ_A是幅度A的标准差；m_A为时序SAR影像的幅度均值；D_A为离差指数，当D_A小于阈值时，选取为PS待定点。In the formula,

is the phase standard deviation; σ _A is the standard deviation of the amplitude A; m _A is the amplitude mean of the time series SAR image; D _A is the dispersion index. When D _A is less than the threshold, it is selected as the PS to be determined.

Amplitude dispersion has a high success rate in selecting PS points with strong amplitude, but for scatterers with low signal-to-noise ratio, the relationship between amplitude dispersion and phase stability is low, so further selection is required.

The phase analysis method mainly uses the temporal and spatial correlation characteristics of different phase components, uses low-pass filtering to separate the noise phase, and selects the pixels less affected by noise as PS points.

为滤波得到的相位；

为由基线引起的误差；N为SAR影像数量。γ反映了目标点在时间序列上的稳定性，γ越大表明待定PS点的噪声越小，因此选取高于阈值的点PS点。该阈值的确定通常需要迭代及滤波运算实现。In the formula,

is the phase obtained by filtering;

is the error caused by the baseline; N is the number of SAR images. γ reflects the stability of the target point in the time series. The larger the γ, the smaller the noise of the PS point to be determined. Therefore, the point PS point higher than the threshold is selected. Determination of the threshold usually requires iteration and filtering operations.

S4, space-time unwrapping, adding the time dimension on the basis of the two-dimensional plane, forming a three-dimensional space search path, and determining residual points for three-dimensional unwrapping. Specifically include:

After removing the terrain phase, the single-point phase composition in the interference pair is expressed as:

W{φ}=W{φ _d +φ _τ +φ _a +φ _o +φ _n } (4)

后，将其去除得到：where W{·} is the phase winding, φ is the single-point interference phase, and φ _d , φ _τ , φ _a , φ _o , and φ _n represent the deformation, terrain error, atmosphere, orbit, and noise phase, respectively. Among them, φ _d , φ _a and φ _o have spatial correlation, which can be separated by spatial low-pass filtering, and the spatial correlation phase can be estimated

Then, remove it to get:

In the formula, the superscript u represents spatial uncorrelation, and δ represents the total amount of unremoved spatially uncorrelated deformation, atmosphere, and orbital errors; because these three have strong spatial correlation, the value of δ is very small;

Construct the coherence coefficient:

Since δ is small, the coherence coefficient γ can reflect the influence of noise on a single point. The larger the noise mean, the larger the coherence coefficient value. The maximum value of γ can be calculated by using the spatial search method, and the corresponding

误差后，进行时空三维解缠。in removing

After the error, the space-time 3D unwrapping is performed.

Traditional time-series InSAR methods mostly use time-dimension (1D) and space-dimension (2D) unwrapping algorithms alone, or 1D+2D pseudo-three-dimensional unwrapping, which cannot consider the spatiotemporal characteristics of the winding phase at the same time. This embodiment adopts a true three-dimensional (3D) unwrapping algorithm, that is, the residual point search path is extended from the plane to the space. This method searches for the residual path in the three-dimensional space and generates residual points. Compared with the plane branch-cut method, the three-dimensional integral path generated by this method can better ensure the accuracy of unwrapping.

The three-dimensional unwrapping is reflected in the space and time dimensions respectively. Based on the minimized L-P norm framework, the phase unwrapping is regarded as an optimization problem, and the minimum value of the objective function is found to obtain the unwrapped phase value. The objective function is:

In the formula, Δφ is the unwrapping phase gradient, Δψ is the winding phase gradient, x and y represent the two two-dimensional directions, z represents the third-dimensional direction, i and j represent the point coordinates, and w represents the weight.

The difference from the traditional L-P norm distribution is that the phase gradient information of the third dimension is added here, usually time series information. Here, the parameter P determines how to deal with the difference between Δφ and Δψ. When P=0, if only two-dimensional information is considered, the above formula is the L-0 norm objective function similar to the branch-cut unwrapping algorithm; when P=2 is the least squares unwrapping algorithm.

For nonlinear deformation, in the case of insufficient prior information, the deformation is usually assumed to be a linear model in the deformation calculation. However, it is difficult to ensure that the assumption is established in the actual deformation process, which affects the calculation of the elevation and deformation parameters. In the case of periodic deformation and complex deformation, the function model has a greater influence on the DEM error estimation. Aiming at the problem of model error, iterative processing is adopted, and the spatial and temporal characteristics of terrain, atmosphere and deformation are used to separate the deformation signal by means of high and low pass filtering in the temporal and spatial domain, which reduces the influence of the model on parameter estimation, thereby ensuring the reliability of deformation extraction. . In addition, the deformation signal has the characteristics of spatial correlation, and the deformation of adjacent points is similar. In order to avoid the errors introduced by the prior deformation model, space deformation constraints can also be added to directly invert the deformation variables.

S5. Error separation. For the layered atmosphere and the floating atmosphere, the methods of parameterization and space-time filtering are used to separate the atmospheric errors.

Among the many factors affecting the accuracy of InSAR land subsidence monitoring, atmospheric delay is one of the main errors. Considering the complex and changeable meteorological conditions in the survey area, there may be serious atmospheric delay errors in the InSAR interferogram. The atmospheric delay is mainly affected by the troposphere. For the tropospheric layered atmosphere, it can be absorbed by an elevation-related parameter because of its strong regularity. The floating atmosphere is affected by the change of water vapor, and the model estimation accuracy is low. This embodiment provides processing solutions for the two parts of the delay respectively.

Stratified atmospheres: Terrain-dependent atmospheric delay phases caused by vertical stratification effects can be effectively separated using elevation phase correlation analysis. According to the principle of the vertical stratification effect, it is generally believed that the atmospheric delay phase has a linear relationship with the elevation. Therefore, a linear regression model is used to establish the relationship between the elevation and the phase to invert the atmospheric delay phase, or a quadratic or high-order polynomial can be used for fitting. The traditional method estimates the stratified atmosphere separately, but due to the existence of deformation and terrain error phase, there will be errors in the estimation of the delay parameter K; in addition, the spatial variation of the stratified atmosphere is not consistent, so a single parameter cannot be applied to Estimated delay for the entire interferometric pair. Aiming at the above two problems, the method of joint estimation of layered atmospheric delay and deformation and terrain parameters is adopted, and atmospheric parameters are added to the model for absorption.

Δφ=Δφ _d+τ +K×Δh

In the formula, Δ represents double-difference observation, Δφ _d+τ is the sum of deformation and terrain phase, and Δh is the elevation difference. Interference solves the partition, which not only reduces the computational burden but also ensures the accuracy of parameter estimation.

Floating atmosphere: This part of the atmosphere is random in space and time, and it is difficult to describe it with a deterministic model. In particular, the low latitude areas are seriously affected by atmospheric images in summer, and are mostly affected by the floating atmosphere. Therefore, on the basis of ensuring a sufficient number of images, individual severely de-correlated interference pairs can be appropriately deleted to reduce the influence of atmospheric delay. The residual atmosphere is separated by its statistical characteristics, that is, according to the correlation difference between the atmospheric delay and the deformation phase in the space-time dimension, mainly through the method of space-time filtering.

The error separation specifically includes:

The unwrapped phase value is:

为残余具有空间相关性的地形残差相位，k为单差模糊度整数。如果三维解缠正确，则大部分点的k应该相等。为分离出形变相位，需要在弧段上进行时空域滤波。where φ _uw is the unwrapped phase,

is the residual terrain residual phase with spatial correlation, and k is a single-difference ambiguity integer. If the 3D unwrapping is correct, the k should be equal for most of the points. In order to separate out the deformation phase, it is necessary to perform spatiotemporal filtering on the arc segment.

代表两点间做差，然后构建三角网进行时间域低通滤波，将其分离后得到In the formula,

represents the difference between two points, and then constructs a triangular network for low-pass filtering in the time domain, and separates them to obtain

和

分别表示副影像大气和轨道分量；

在时间域上表现为较强的相关性，利用时间域低通可以获取形变和主影像大气；其余量在时间域上表现为不相关，要获取副影像大气，

在空间表现为相关性，

表现为随机噪声误差，因此可以通过空间的高通滤波分离出副影像的大气延迟相位和轨道误差；最后通过从弧段到点上的积分可以获取每个点的φ_d及

In the formula, [ ] _{LP_time} represents the time domain low-pass filtering,

and

represent the atmospheric and orbital components of the secondary image, respectively;

There is a strong correlation in the time domain, and the deformation and the atmosphere of the main image can be obtained by using the low-pass in the time domain; the remaining quantities are irrelevant in the time domain. To obtain the atmosphere of the auxiliary image,

Shows correlation in space,

It is a random noise error, so the atmospheric delay phase and orbital error of the secondary image can be separated by high-pass filtering in space; finally, the φ _d and

Further, the surface deformation inversion method also includes GNSS-InSAR fusion. After obtaining the unwrapped phase, if there is external GNSS data, obtain the ZWD from the GNSS data, invert the atmospheric delay, and remove it from the unwrapped phase.

GNSS has the characteristics of high temporal sampling rate, which can provide higher temporal resolution and higher precision ZWD estimation, and reduce or even avoid the uncertainty caused by the inconsistency between water vapor products and radar imaging time. However, GNSS reflects the single-point ZTD estimation, and applying it to InSAR atmospheric phase correction requires spatial interpolation, so there are higher requirements on the density and distribution of GNSS stations.

The DEM error phase is added back and re-parameterized. The least squares method is used to calculate the deformation phase. The GNSS observation equation is used as a constraint, and the joint adjustment is performed with the InSAR observation equation. Specifically, the GNSS prior monitoring results are used as strong constraints to construct the constraints of deformation calculation, which can improve the calculation accuracy of least squares. There are three schemes for the fusion of this constraint and InSAR data: (1) use GNSS observations as conditional constraints to constrain the function model; (2) use the covariance matrix of GNSS observations to solve additional random model constraints; (3) ) performs the fusion solution of the dual constraints of the function and the stochastic model at the same time. The constrained error equation can be written as:

In the formula, L _GNSS and L _InSAR are the GNSS and InSAR observation values, respectively, ε _GNSS and ε _InSAR are the observation errors, X is the parameter to be estimated, including the deformation and terrain parameters, and the final deformation information is obtained by the least squares solution.

The above-mentioned surface deformation inversion method based on time series InSAR technology can monitor the surface deformation in all weather, large area, continuous and high spatial resolution. Specifically, in the image registration method, coarse registration and fine registration are used to ensure the registration accuracy required in different imaging modes; for the atmospheric delay error, different algorithms are used to process the layered atmosphere and the floating atmosphere, respectively. The influence of atmospheric errors is weakened to the greatest extent, and the accuracy of deformation inversion is improved; if GNSS data is provided, the atmospheric delay can be corrected by using GNSS data, and the joint adjustment is performed as a constraint in the InSAR observation equation to ensure the accuracy of the deformation results. Precision and reliability. Compared with the traditional time sequence technology, the present invention has strong advantages in the registration method, the atmospheric delay processing algorithm and the GNSS-InSAR fusion technology respectively.

The above descriptions are only preferred embodiments of the present invention, which are not intended to limit the protection scope of the present invention. For those skilled in the art, the present invention may have various modifications and changes. Any changes, modifications, substitutions, integrations and parameter changes to these embodiments without departing from the principles and spirit of the present invention, through conventional substitutions or capable of achieving the same function within the spirit and principles of the present invention, all fall within the scope of the present invention. into the protection scope of the present invention.

Claims (10)

1. A surface deformation inversion method based on a time sequence InSAR technology is characterized by comprising the following steps:

s1, geocoding and image registration, wherein an external DEM is used for geometric registration, correlation registration is carried out based on the amplitude and intensity information of an image, and a method for enhancing spectral diversity is used for fine registration;

s2, generating a differential interference pair;

s3, selecting a PS point, and performing iterative selection of the PS point by respectively using the image intensity and the phase information;

s4, performing space-time unwrapping, increasing time dimension on the basis of a two-dimensional plane to form a three-dimensional space search path, and determining residual points to perform three-dimensional unwrapping;

and S5, error separation, namely, carrying out atmospheric error separation by adopting parameterization and space-time filtering modes respectively for layered atmosphere and floating atmosphere.

2. The surface deformation inversion method based on the time series InSAR technology as claimed in claim 1, characterized in that in the step S1, an external DEM is used for geometric registration, specifically:

converting the external DEM into an image coordinate system of a main image, establishing a lookup table, then calculating the position of a certain pixel through satellite orbit parameters, setting a proper window size in the main image, uniformly selecting pixel points, performing simultaneous solution on the coordinates of the pixel points through a distance equation, a Doppler equation and an earth model equation, and inversely calculating the coordinates of the pixel points in an auxiliary image coordinate system;

and calculating the offsets of all the pixel points by utilizing Delaunay triangulation network interpolation, and correcting the secondary image by utilizing a final offset file.

3. The earth surface deformation inversion method based on the time series InSAR technology as claimed in claim 1, wherein the step S1 is performed with correlation registration based on the amplitude and intensity information of the image, specifically:

selecting a proper window in the main image, and calculating the coherence coefficient of the window;

utilizing maximum likelihood estimation to carry out offset between the main image and the auxiliary image;

and obtaining a conversion matrix according to the calculated offset, carrying out registration on the secondary image, and simultaneously ensuring that the registration precision reaches 0.05 pixel.

4. The surface deformation inversion method based on the time series InSAR technology as claimed in claim 1, wherein the fine registration for enhancing the spectrum diversity in step S1 specifically comprises:

extracting an overlapping area between the main image and the resampled auxiliary image burst;

the main and auxiliary images are interfered to form a front view and a rear view of an overlapped area, and are interfered again and filtered;

estimating ESD phase and registration errors;

and fitting all the deviation of the overlapping areas by using a polynomial to calculate the azimuth offset.

5. The surface deformation inversion method based on the time series InSAR technology as claimed in claim 1, wherein the step S2 of generating differential interference pairs comprises:

the time base line, the space base line and the frequency base line formed by the main image and the secondary image are important factors influencing the interference effect, and in order to select the optimal public main image, the sum of the correlation coefficients of the time sequence interference pattern needs to be maximized:

wherein γ is a coherence coefficient, f (γ)_t)，f(γ_B)，f(γ_DC) Respectively obtaining a time baseline coefficient, a vertical baseline coefficient and a Doppler mass center frequency difference coefficient, and forming an interference pair after finding the optimal main image;

and the external DEM is transferred to a radar coordinate system to form a terrain phase, and the terrain phase and the flat ground phase are subtracted from the interference phase to form a differential interference pair.

6. The surface deformation inversion method based on the time sequence InSAR technology as claimed in claim 1, wherein the step S3 of selecting the PS point comprises:

roughly selecting by using the amplitude dispersion coefficient, and reserving a point target higher than a threshold value;

calculating a coherence coefficient of a phase point in time through a differential phase, a filtering phase and a terrain residual phase after differential interference, selecting a proper coherence coefficient threshold, and identifying a point higher than the threshold as a PS target point.

7. The surface deformation inversion method based on the time series InSAR technology as claimed in claim 1, wherein the space-time unwrapping in the step S4 includes:

after removing the topographic phase, the single point phase in the interferometric pair is represented as:

W{φ}＝W{φ_d+φ_τ+φ_a+φ_o+φ_n}

wherein W {. is phase winding, phi is single-point interference phase, phi_d、φ_τ、φ_a、φ_o、φ_nRespectively representing deformation, terrain error, atmosphere, orbit, and noise phase; wherein phi_d、φ_aAnd phi_oHaving correlation in space, the correlation phase can be estimated by separating it with space low-pass filtering

After that, it was removed to give:

in the formula, the superscript u represents spatial irrelevance and represents the total amount of deformation, atmosphere and orbit errors which are not removed and are spatially uncorrelated; the three have strong spatial correlation, so the value is small;

constructing a coherent coefficient:

because the interference coefficient gamma is very small, the interference coefficient gamma can reflect the influence of noise on a single point, the larger the noise mean value is, the larger the interference coefficient value is, the maximum value of gamma can be calculated by using a space search method, and the corresponding maximum value can be calculated

In removing

And after the error, performing space-time three-dimensional unwrapping.

8. The surface deformation inversion method based on the time series InSAR technology as claimed in claim 1, wherein the error separation in the step S5 includes:

the phase values after unwrapping were:

in the formula, phi_uwIn order to be able to unwind the phase,

the residual terrain residual error phase with spatial correlation is obtained, and k is a single-difference ambiguity integer; if the three-dimensional unwrapping is correct, k should be equal for most points; in order to separate out the deformation phase, time-space domain filtering is required to be carried out on the arc section;

in the formula (I), the compound is shown in the specification,

representing the difference between two points, then constructing a triangular network to carry out time domain low-pass filtering, and separating the triangular network to obtain

In the formula [ ·]_{LP_time}Which represents a time domain low-pass filtering,

and

respectively representing the atmosphere and the orbit components of the secondary image;

the method has strong correlation in a time domain, and deformation and main image atmosphere can be obtained by utilizing time domain low pass; the remaining amount of the image is irrelevant in the time domain, the sub-image atmosphere is acquired,

the correlation is shown in the space, and,

the error is expressed as a random noise error, so that the atmospheric delay phase and the orbit error of the secondary image can be separated through spatial high-pass filtering; finally, the phi of each point can be obtained by integrating from the arc segment to the point_dAnd

9. the earth surface deformation inversion method based on the time sequence InSAR technology as claimed in claim 1, characterized in that if the differential phase still has a relatively obvious long wave trend after the high-low pass filtering in the error separation, i.e. the layered atmosphere delay interference, the linear function is adopted for absorption and removal.

10. The earth surface deformation inversion method based on the time sequence InSAR technology as claimed in claim 1, characterized in that the earth surface deformation inversion method further comprises GNSS-InSAR fusion, after the unwrapping phase is obtained, if there is GNSS external data, ZWD can be obtained preferably through GNSS data, atmospheric delay is inverted and removed from the unwrapping phase;

and adding back the DEM error phase, carrying out parameterization again, carrying out deformation phase solution by using least square, and carrying out combined adjustment on a GNSS observation equation and an InSAR observation equation to obtain final deformation information.

Images