Modeling method for inverting surface high-precision three-dimensional deformation through InSAR and GNSS data fusion

Abstract

The invention discloses a modeling method for inverting surface high-precision three-dimensional deformation through InSAR and GNSS data fusion, and belongs to the technical field of surface deformation monitoring. The method includes the steps of resolving InSAR line-of-sight deformation, unifying InSAR and GNSS data reference systems, establishing a stress model based on an elastic theory, calculating a three-dimensional deformation component for the first time based on a least square theory, obtaining an optimal weight matrix by applying a variance component estimation method, calculating final three-dimensional deformation and the like. According to the method, three-dimensional deformation and strain parameters are estimated by fusing InSAR measurement data and GNSS measurement data based on the elastic theory, the stress relation between adjacent points is considered, and unnecessary spatial interpolation of GNSS measurement values is avoided. The variance component estimation method solves the problem of unreasonable weight ratio of different types of observation values, and according to the adjusted optimal weight matrix, a function model is solved based on the least square method, and high-precision and high-spatial-resolution earth surface three-dimensional deformation information within the research range is obtained.

Description

Technical Field

[0001] The invention belongs to the technical field of surface deformation monitoring, and specifically is a modeling method for inverting high-precision three-dimensional surface deformation by fusing InSAR and GNSS data. Background Art

[0002] In recent years, geological disasters such as landslides and ground subsidence caused by various natural or human factors have occurred frequently, seriously affecting the safety of public life and public property. The construction of large-scale underground projects has a particularly significant impact on surface deformation. As two very mature technologies for monitoring surface deformation, GNSS (Global Navigation Satellite System) and InSAR (Synthetic Aperture Radar Interferometry) technologies have shown different advantages in their respective fields. GNSS technology is a technology that can obtain real-time high-precision three-dimensional positions of monitoring stations, with the highest accuracy reaching millimeter level. Due to the characteristics of GNSS, obtaining high-precision observation results requires long-term deployment and continuous observation of monitoring stations. In addition, the expensive monitoring equipment leads to low spatial resolution of GNSS. InSAR technology has the advantages of all-day, all-weather, large range and high spatial resolution, but the deformation obtained by InSAR technology is line-of-sight (LOS) deformation, which is a one-dimensional observation and has low temporal resolution. It can be seen that the two technologies have different advantages and limitations. If the monitoring results of the two technologies are integrated and fused, and a more accurate three-dimensional inversion of the deformation area is performed through the association of the two data in the spatial model, the deformation characteristics of the study area can be analyzed more deeply.

[0003] At present, there are many methods for inverting the three-dimensional deformation of the surface by fusing InSAR and GNSS data. For example, the combination method uses the Gibbs energy equation as the objective function to obtain the optimal solution of the three-dimensional deformation when the function energy is minimum; another example is the direct decomposition method, which uses the high-precision horizontal deformation displacement provided by GNSS observations as a constraint to decompose the LOS displacement.

[0004] The above methods mainly focus on the deformation calculation of a single point, without considering the spatial correlation between the surface deformations of neighboring points, and cannot meet the relevant requirements in elastic deformation theory. Summary of the invention

[0005] In view of the above problems existing in the prior art, the purpose of the present invention is to provide a modeling method for inverting high-precision three-dimensional deformation of the surface by fusing InSAR and GNSS data. The three-dimensional deformation and strain parameters are estimated by fusing InSAR and GNSS measurement data based on elastic theory, and the stress connection between adjacent points is considered to avoid unnecessary spatial interpolation of GNSS measurement values. The process applies the variance component estimation method to solve the unreasonable weight ratio of different types of observations. According to the adjusted optimal weight matrix and based on the least squares principle, the function model is solved to obtain high-precision and high-spatial-resolution three-dimensional deformation information of the surface within the research range.

[0006] In order to achieve the above object, the technical solution adopted by the present invention is:

[0007] The modeling method of inverting high-precision three-dimensional deformation of the surface by fusing InSAR and GNSS data includes the following steps:

[0008] S1: Use DInSAR technology to process SAR data, set reasonable parameters and perform interference operation to obtain the interference phase;

[0009] S2: Atmospheric delay error correction to obtain the line-of-sight deformation map of the ground surface at the monitoring location;

[0010] S3: After removing invalid values ​​from the sight deformation map, output it as a raster and encode it so that it is in the same reference system as the GNSS data used;

[0011] S4: Calculate the conversion parameters between the line-of-sight deformation and the three-dimensional surface deformation components, and establish the observation equation;

[0012] S5: Select the GNSS data and InSAR ascending and descending orbit data points of the target point, construct the three-dimensional stress decomposition equation, and solve the three-dimensional deformation components and stress parameters of the target point based on the least squares method;

[0013] S6: Using the initially calculated observation corrections, estimate the pre-test variances of each observation after the test until the unit weight variances of each type of observation are equal or approximate to obtain the optimal weight matrix;

[0014] S7: Calculate the final three-dimensional deformation components of the target point based on the least squares method;

[0015] S8: Traverse the target area points and repeat steps S5 to S7 to invert the three-dimensional deformation.

[0016] As a further improvement of the above technical solution:

[0017] In step S1, the DInSAR technology is used to perform interference processing on the SAR image pair, the minimum interference flow method is used to unwrap the interference phase, a polynomial model is used for the unwrapped phase, and the polynomial coefficients are solved in combination with the ground stable point phase, and the optimized unwrapped phase is obtained after re-flattening.

[0018] In step S1, a multi-view ratio of 3:1 is adopted when performing interference processing.

[0019] In step S2, GACOS uses an iterative tropospheric decomposition model to decompose the tropospheric delay to obtain high-resolution tropospheric zenith delay data, converts it to the line-of-sight direction and removes it from the unwrapped phase, and finally converts the unwrapped phase into a line-of-sight deformation map that is deformed to the surface of the monitoring location.

[0020] In step S5, a model with stress relationship is established between local observations based on elastic theory, and the InSAR orbit raising data, InSAR orbit lowering data, and GNSS data are fused, and the selection ratio of various data points can be adjusted.

[0021] The variance component estimation method in step S6 estimates the unit weight variance of each type of observation value by using the correction numbers after the initial adjustment of each type of observation value and makes them equal or approximate, and the observation value weight matrix obtained in this way is the optimal weight matrix.

[0022] In step S7, the updated optimal weight matrix is ​​used to recalculate the vector according to the least square method to obtain the three-dimensional deformation components and various stress parameters of the target point.

[0023] In step S8, the points to be solved in the entire target area are traversed, and the calculated deformation results are combined to obtain high-precision three-dimensional deformation information of the entire study area.

[0024] The beneficial effects of the present invention are:

[0025] (1) The modeling method is highly scientific and applicable. It obtains the overall three-dimensional deformation information by inverting the spatial stress-related deformation and combining the advantages of the two types of data.

[0026] (2) The modeling method not only combines the advantages of multiple data, but also takes into account the spatial relationship between adjacent data points in the surface space. It can obtain high-precision three-dimensional deformation in the study area based on the inversion of elastic deformation theory.

[0027] (3) Based on elastic theory, the three-dimensional deformation and strain parameters are estimated by fusing InSAR and GNSS measurement data, taking into account the stress relationship between adjacent points and avoiding unnecessary spatial interpolation of GNSS measurement values. The variance component estimation method is used to solve the unreasonable weight ratio of different types of observations. The function model is solved based on the least squares principle according to the adjusted optimal weight matrix to obtain high-precision and high-spatial-resolution three-dimensional surface deformation information within the research scope. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 It is a schematic diagram of the process of the present invention. DETAILED DESCRIPTION

[0029] The specific implementation of the present invention is described in detail below in conjunction with the accompanying drawings. It should be understood that the specific implementation described here is only used to illustrate and explain the present invention, and is not used to limit the present invention.

[0030] For ease of description, spatially relative terms such as "above", "above", "on the upper surface of", "above", etc. may be used here to describe the spatial positional relationship between a device or feature and other devices or features as shown in the figure. It should be understood that spatially relative terms are intended to include different orientations of the device in use or operation in addition to the orientation described in the figure. For example, if the device in the accompanying drawings is inverted, the device described as "above other devices or structures" or "above other devices or structures" will be positioned as "below other devices or structures" or "below other devices or structures". Thus, the exemplary term "above" can include both "above" and "below". The device can also be positioned in other different ways (rotated 90 degrees or in other orientations), and the spatially relative descriptions used here are interpreted accordingly.

[0031] A modeling method for inverting three-dimensional surface deformation by fusing InSAR and GNSS data, such as Figure 1 As shown, the method includes the following steps: solving the InSAR line-of-sight deformation, unifying the InSAR and GNSS data reference systems, establishing a stress model and initially calculating the three-dimensional deformation components, obtaining the optimal weight matrix and calculating the final three-dimensional deformation results.

[0032] The following Step 1 and Step 2 are the steps to solve the InSAR line of sight deformation:

[0033] Step 1: Use DInSAR (differential interferometry) technology to perform interference processing on the SAR (synthetic aperture radar) image pairs before and after the deformation of the target area, and set a reasonable multi-view ratio to meet the mapping resolution. In this embodiment, the multi-view ratio is 3:1, and the terrain phase is removed with the help of DEM (digital elevation model) data. After phase filtering, the minimum flow method is used for unwrapping, the unwrapping decomposition level is set to level 0, and the coherence threshold is set to 0.3. A polynomial model is used for the unwrapped phase, and the polynomial coefficients are solved in combination with the ground stable point phase. After re-leveling, the optimized unwrapped phase is obtained. Among them, the terrain phase and the unwrapped phase belong to the interference phase, and the terrain phase is the error that needs to be removed in the interference phase; the unwrapped phase is calculated from the interference phase after the error is removed.

[0034] Step 2: Use GACOS data for atmospheric correction. GACOS uses an iterative tropospheric decomposition model to decompose the tropospheric delay, which can obtain high-resolution tropospheric zenith delay data, and then convert it to the line of sight to remove it from the unwrapped phase. Finally, the unwrapped phase is converted into a line of sight deformation LOS (Line of sight) map of the surface at the monitoring location.

[0035] The following Step 3 and Step 4 are unified steps for the InSAR and GNSS data reference systems:

[0036] Step 3: Remove the null or invalid values ​​from the obtained line-of-sight deformation map and perform geocoding to make it in the same reference system as the GNSS data used. After removing the invalid values, output it as a raster and encode it. The raster output format assigns the deformation data to the pixels within the graphics range in a planar manner, which is opposite to the vector format.

[0037] Assume that there is a target point in the target area Its shape Neighbors Its shape in, The target point x 0 The coordinates in the E, N, and U directions respectively, The target point x 0 The deformation components in the E, N, and U directions, is the neighboring point x i The coordinates in the E, N, and U directions respectively, is the neighboring point x i The deformation components in the E, N, and U directions respectively. E, N, and U are three mutually perpendicular directions. Specifically, E is the east-west direction, N is the north-south direction, and U is the perpendicular direction perpendicular to E and N.

[0038] Step 4: Assume that the InSAR line-of-sight deformation observation is L 1 , the target point observation value is L0 GNSS observation L 2 , InSAR observation value and its corresponding position three-dimensional deformation u i The relationship is

[0039]

[0040] In the above formula, T represents the transposed matrix of the corresponding matrix, They are x i The deformation components of the point in the directions of E, N, and U. The above formula is the established observation equation.

[0041] Among them, S 1 ,S 2 ,S 3 are conversion parameters determined by the satellite azimuth and incident angle when acquiring InSAR data, and their respective calculation formulas are:

[0042]

[0043] In the formula, α and θ represent the azimuth and the angle of incidence, respectively. To ensure the accuracy of the conversion parameters, the geocoded overall line of sight incidence angle graphic data ILOS and line of sight azimuth data ALOS are used in the calculation. This process calculates parameters for each target point separately, which can reduce the systematic error caused by coordinate conversion.

[0044] The following Step 5 is the steps for establishing the stress model and initial calculation:

[0045] In Step 5, a stress model is established based on elasticity theory, and the three-dimensional deformation components are initially calculated based on the least squares theory.

[0046] In Step 5, the selected neighbor point x i The data include InSAR ascending data, InSAR descending data, and GNSS data.

[0047] Step 5: Establish a model with stress relationship between local observations based on elastic theory i =H·Δ i +u 0 , where Δ i =x i -x 0 =[Δx E(i) Δx N(i) Δx U(i) ] represents the coordinate increment matrix, H = Ε + Ω represents the strain parameter matrix, where Ε and Ω represent the strain tensor and rigid body rotation tensor respectively, u iis the shape variable of the neighboring points involved in the solution. The formula is rewritten as L = B x, where B represents the design matrix in the stress-strain model, and x represents the parameter vector at the target point. Their respective expressions are:

[0048] B=[B 1 T ,B 2 T ,B 3 T ] T

[0049] B 1 is the design matrix corresponding to InSAR data. E(i) , Δx N(i) , Δx U(i) is the midpoint x of the coordinate increment matrix i The components of the increment in the E, N, and U directions respectively.

[0050]

[0051] B 2 Design matrix corresponding to GNSS data:

[0052]

[0053] B 3 Design matrix for target point correspondence:

[0054] B 3 =[S 1 S 2 S 3 000 000 000]

[0055] Solution vector x:

[0056] x=[u E(0) u N(0) u U(0) ξ 11 ξ 12 ξ 13 ξ 22 ξ 23 ξ 33 ω 1 ω 2 ω 3 ] T

[0057] Vector x represents the target point x 0 The parameter vector at which u E(0) 、u N(0) 、u U(0) are the displacements of the target point in the east-west, north-south and vertical directions, ξ 11 ,12 , 13 , 22 , 23 , 33 represents the strain parameter, ω 1 ,ω 2 ,ω 3 is the rigid body rotation parameter.

[0058] Observation L i They represent different types of observations:

[0059]

[0060] L 3 =[L 0 ] T

[0061] Establish the observation error equation:

[0062] V i =B i xL i

[0063] Where i = 1, 2, 3, V i Table error.

[0064] The equation established for each data point in the region is solved based on the least squares theory.

[0065] Establish the normal equation:

[0066] N i =B i T P i B i

[0067] W i =B i T P i L i

[0068] Among them, N i , Wi is the formula of the normal equation, Pi is the weight matrix, i=1,2,3.

[0069] Where P is the weight matrix of each observation. Considering the difference between the horizontal observation accuracy and the vertical observation accuracy of GNSS, the weight matrix of GNSS observation value can be multiplied by the proportional coefficient.

[0070]

[0071] Where d i is the i-th ground neighbor point or GNSS station x i With the target point x0 The distance between.

[0072]

[0073] Where d0 is the distance decay constant, which defines the estimated locality level.

[0074]

[0075] d ij Represents the distance between the i-th GNSS station and the j-th neighboring data point. M represents the number of GNSS stations, and K represents the number of neighboring data points of the target point.

[0076] Solve the algorithm equation to get the stress parameter (ξ 11 , 12 , 13 , 22 , 23 , 33 ,ω 1 ,ω 2 ,ω 3 ) and the three-dimensional deformation component (u E(0) 、u N(0) 、u U(0) ):

[0077] x=N -1 W

[0078] in

[0079] N=N 1 +N 2 +N 3

[0080] W=W 1 +W 2 +W 3

[0081] Where N is N 1 、N 2 、N 3 The sum of the three, W is W 1 , W 2 , W 3 The sum of the three.

[0082] In Step 5, u i =H·Δ i +u 0 and B 1 , B 2 , B 3 is the three-dimensional stress decomposition equation. By designing the matrix B 1 , B 2 , B 3Various types of data, such as InSAR ascending data, descending data, and GNSS data, are integrated, and various data are linked together to uniformly solve the target point parameter vector x; the selection ratio of various data points is adjusted by adjusting M and K.

[0083] The following Step 6 and Step 7 are steps for applying the variance component estimation method to obtain the optimal weight matrix and calculate the final three-dimensional deformation:

[0084] Step 6: Apply the variance component estimation algorithm, that is, use the sum of squares of the corrections after adjustment of various observations To estimate is the unit weight variance of the observation vector of type i. To do this, we first need to establish the relationship between the residual sum of squares and the unit weight variance:

[0085] According to the expectation theorem of quadratic forms:

[0086] E(V i T P i V i )=tr(P i D(V i ))

[0087] Among them, tr represents the trace number, E is the expectation, D(V i ) is the variance-covariance matrix.

[0088] When the number of observed values ​​is n:

[0089]

[0090] From the covariance propagation rate we get:

[0091]

[0092] Get the relationship:

[0093]

[0094] In the above formula, ni is equal to D(V i )Matrix dimensions.

[0095] Written in matrix form is the variance estimation equation:

[0096]

[0097] Where S represents the transformation matrix:

[0098]

[0099] W θ Represents a quadratic vector of observation corrections:

[0100]

[0101] is an estimate of the unit weight variance of each type of observation vector.

[0102]

[0103] The solution is:

[0104]

[0105] After obtaining the variance estimate, make the first estimate of the observation weights:

[0106]

[0107] Where c is any constant, usually chosen After the weight is updated, the adjustment is repeated. Then perform component estimation; iterate repeatedly until Or when the ratio of each variance is approximately 1, the optimal weight matrix is ​​determined. This method of estimating the variance of each type of observation and determining the weight by the a posteriori method has a certain reliability.

[0108] Step 7: Use the updated optimal weight matrix P to recalculate the vector x according to the least square method to obtain the three-dimensional deformation components and various stress parameters of the target point. Traverse the points to be solved in the entire target area, and combine the deformation results calculated according to the above steps to obtain high-precision three-dimensional deformation information of the entire study area.

[0109] The above disclosed embodiments are described to illustrate the present invention and enable professionals in the field to implement or use the present invention. Some embodiments do not describe the details in detail, and various modifications of these embodiments and implementations in other embodiments are foreseeable to those skilled in the art. Therefore, any modifications made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A modeling method for inverting high-precision three-dimensional deformation of the surface by fusing InSAR and GNSS data, characterized in that: The following steps are involved: S1: Use DInSAR technology to process SAR data, set reasonable parameters and perform interference operation to obtain the interference phase; S2: Atmospheric delay error correction to obtain the line-of-sight deformation map of the ground surface at the monitoring location; S3: After removing invalid values ​​from the sight deformation map, output it as a raster and encode it so that it is in the same reference system as the GNSS data used; S4: Calculate the conversion parameters between the line-of-sight deformation and the three-dimensional surface deformation components, and establish the observation equation; S5: Select the GNSS data and InSAR ascending and descending orbit data points of the target point, construct the three-dimensional stress decomposition equation, and solve the three-dimensional deformation components and stress parameters of the target point based on the least squares method; S6: Using the initially calculated observation corrections, estimate the pre-test variances of each observation after the test until the unit weight variances of each type of observation are equal or approximate to obtain the optimal weight matrix; S7: Calculate the final three-dimensional deformation components of the target point based on the least squares method; S8: Traverse the target area points and repeat steps S5 to S7 to invert the three-dimensional deformation.

2. The modeling method according to claim 1, characterized in that: In step S1, the DInSAR technology is used to perform interference processing on the SAR image pair, the minimum interference flow method is used to unwrap the interference phase, a polynomial model is used for the unwrapped phase, and the polynomial coefficients are solved in combination with the ground stable point phase, and the optimized unwrapped phase is obtained after re-flattening.

3. The modeling method according to claim 2, characterized in that: In step S1, a multi-view ratio of 3:1 is adopted when performing interference processing.

4. The modeling method according to claim 1, characterized in that: In step S2, GACOS uses an iterative tropospheric decomposition model to decompose the tropospheric delay to obtain high-resolution tropospheric zenith delay data, converts it to the line-of-sight direction and removes it from the unwrapped phase, and finally converts the unwrapped phase into a line-of-sight deformation map that is deformed to the surface of the monitoring location.

5. The modeling method according to claim 1, characterized in that: In step S5, a model with stress relationship is established between local observations based on elastic theory, and the InSAR orbit raising data, InSAR orbit lowering data, and GNSS data are fused, and the selection ratio of various data points can be adjusted.

6. The modeling method according to claim 1, characterized in that: The variance component estimation method in step S6 estimates the unit weight variance of each type of observation value by using the correction numbers after the initial adjustment of each type of observation value and makes them equal or approximate, and the observation value weight matrix obtained in this way is the optimal weight matrix.

7. The modeling method according to claim 1, characterized in that: In step S7, the updated optimal weight matrix is ​​used to recalculate the vector according to the least square method to obtain the three-dimensional deformation components and various stress parameters of the target point.

8. The modeling method according to claim 1, characterized in that: In step S8, the points to be solved in the entire target area are traversed, and the calculated deformation results are combined to obtain high-precision three-dimensional deformation information of the entire study area.

Images