ISTRE-based deformation monitoring fusion algorithm

Abstract

The invention belongs to the technical field of deformation monitoring, and relates to an ISTRE-based deformation monitoring fusion method, which comprises the following steps: performing vertical correction on SBAS-InSAR data, and estimating ISTR E model parameters by using SBAS-InSAR deformation data; the corresponding SBAS-InSAR data is covered by using the known GNSS data; re-weighting the covariance matrix of the observation equation; obtaining vertical correction data through forward filtering, and calculating the optimal estimation of the state variable at the current moment through a recursive mode by using the state estimation value at the previous moment and the observation value at the current moment; the correction data and the GNSS data are fused, and ISTRE model parameters are estimated by using the correction data and the GNSS data; and carrying out fusion calculation through forward filtering. According to the method, the ISTRE model is provided, fusion calculation of GNSS vertical direction deformation and SBAS-InSAR line-of-sight direction deformation is achieved, and high precision is achieved.

Description

Technical Field

[0001] The invention belongs to the field of deformation monitoring, and in particular relates to a GNSS and SBAS-InSAR deformation monitoring fusion method based on ISTRE. Background Art

[0002] At present, GNSS RTK technology can identify sub-centimeter surface deformation with extremely high temporal resolution (seconds) and can obtain sub-centimeter deformation features. However, it is difficult to deploy GNSS receivers at a high density and cannot obtain deformation fields with high spatial resolution. Satellite-based InSAR deformation monitoring has high accuracy and spatial resolution, but the satellite revisit period determines that the temporal resolution of this technology is low, which limits its application in deformation monitoring with high timeliness requirements. Based on this, GNSS and InSAR can be fused to achieve information complementarity and generate data with high temporal and spatial resolution, which can effectively increase the application scope of deformation monitoring data.

[0003] In the field of GNSS and InSAR fusion research, traditional methods such as analytical optimization method improve the accuracy of InSAR based on GNSS inversion of InSAR orbit error and atmospheric delay, and variance component estimation (VCE) improves the accuracy of GNSS and InSAR fusion by calculating the prior variance of GNSS and InSAR to determine the model weight. However, both types of fusion algorithms focus on the accuracy improvement of InSAR by GNSS technology, and the fusion result is the same as the time resolution of InSAR, abandoning the high time resolution advantage of GNSS. Compared with traditional methods, Kalman filtering algorithm has significant advantages. It establishes state equations and observation equations based on physical mechanisms and the relationship between state variables and observation variables. Kalman filtering can achieve the minimum variance estimation of state variables based on the least squares principle. This fusion method also improves the deformation monitoring accuracy of InSAR points, and the fusion result is the same as the time resolution of GNSS, but it cannot maintain the spatial resolution of InSAR.

[0004] To this end, based on Klalman filtering, researchers proposed a fusion algorithm combined with a spatial autoregressive model (AR). Although the algorithm retains the high temporal resolution of the Kalman filter and the high spatial resolution of the autoregressive model, it ignores the temporal correlation of the data. The spatio-temporal Kalman filter model combined with the kriging technology takes into account both the temporal correlation and spatial correlation of the deformation, and can quickly simulate the spatio-temporal process of the three-dimensional surface. It is an effective tool for constructing the spatio-temporal model of the GNSS network and InSAR deformation field. Based on this, the researchers proposed the Spatio-Temporal Random Effect Model (STRE). Based on the idea of ​​spatio-temporal random effects, the algorithm introduces a fixed-rank spatial basis in the observation equation to reduce the complexity of the model, and can finely analyze both large-scale local variations and small-scale subtle variations. The construction process includes: establishing a spatial field based on InSAR data, correcting the spatial field in combination with GNSS data, and making the model dynamic.

[0005] However, due to the fact that it is usually impossible to obtain vertical deformation through dual-track InSAR data in areas with severe deformation or complex terrain, the STRE model needs to project the three-dimensional deformation monitored by GNSS to the line of sight (LOS) direction of InSAR before fusion to unify the measurement geometry of the two. Deformation in the LOS direction cannot directly reflect the real deformation characteristics of the surface. The weakening of radar signal scattering and interference effects in areas with complex terrain will reduce the measurement accuracy of InSAR, making it difficult to apply InSAR to fields such as geological disaster risk assessment and settlement mechanism interpretation. Using the incident angle of the satellite sensor to project the LOS measurement value to the vertical direction to unify the measurement geometry to the GNSS direction may result in large errors, thus affecting the fusion accuracy.

[0006] Based on this, in order to solve the problem of inconsistent fusion data directions at the current stage and based on the advantages of STRE, it is very necessary to construct a fusion algorithm based on STRE for the fusion calculation of GNSS vertical direction and SBAS-InSAR LOS direction deformation, which has wide application value. Summary of the invention

[0007] The purpose of the present invention is to provide a GNSS and SBAS-InSAR deformation monitoring fusion method based on ISTRE to solve the problem of fusion algorithms with different structural deformation directions and maintain good accuracy.

[0008] The objective of the present invention is achieved through the following technical solutions:

[0009] A GNSS and SBAS-InSAR deformation monitoring fusion method based on ISTRE includes the following steps:

[0010] S1. Perform vertical correction on SBAS-InSAR data. First, use SBAS-InSAR deformation data to estimate ISTRE model parameters. Then, overlay the corresponding SBAS-InSAR data with known GNSS data. Then, re-weight the covariance matrix of the observation equation.

[0011] S2, obtaining vertical correction data through forward filtering, and recursively calculating the optimal estimate of the state variable at the current moment using the state estimate at the previous moment and the observed value at the current moment;

[0012] S3: The correction data is fused with the GNSS data. The ISTRE model parameters are estimated using the correction data and GNSS data. Then, the fusion calculation is performed through forward filtering.

[0013] Furthermore, in step S1, the model parameters are estimated based on the ISTRE model using the SBAS-InSAR deformation data as observation data. After the model parameters are determined, the corresponding SBAS-InSAR data are overlaid with the known GNSS data:

[0014] Z(t)=[Z G，1|t … Z G，j|t Z S,k|t … Z S，N|t ] (4)

[0015] Among them, Z G is the vertical velocity data of GNSS, Z S is the velocity data in the SBAS-InSAR LOS direction, j is the number of GNSS points, k=j+1; the spacing between GNSS deformation monitoring points is small, and the distance between them and the base station is less than 1km, so it is assumed that the monitoring accuracy of each point is the same; SBAS-InSAR monitoring also assumes that the accuracy of each pixel is the same, and the deformation results obtained by the two technologies are highly correlated; therefore, when the GNSS data of the same site is used to replace the SBAS-InSAR data, the spatial field H remains unchanged.

[0016] Furthermore, when the observation data changes, the covariance matrix of the observation equation needs to be re-weighted to ensure the accuracy of the model; the difference between the SBAS-InSAR LOS directional deformation and the vertical correction deformation can be divided into systematic error and noise error. The noise error is corrected by formula (5), and the systematic error is weighted by the covariance matrix of the observation noise based on the concept of robustness to reduce its impact on fusion:

[0017]

[0018] Among them, γ G and γ SThe variance inflation factors for GNSS and SBAS-InSAR are respectively, which can be calculated based on the robustness principle.

[0019] Furthermore, the γ G and γ S GNSS and SBAS-InSAR are the variance inflation factors, respectively. The calculation formula of γ is:

[0020]

[0021] Among them, k 0 and k 1 are the quantile parameter and elimination point respectively. ii is the robustness factor:

[0022]

[0023] in, is the i-th element of the residual vector of the observation vector at time k, is the covariance matrix of the residual vector. 0 and k 1 The k value selected for GNSS and SBAS-InSAR data with different accuracy is determined according to the size of the standardized residual. 0 and k 1 different.

[0024] Furthermore, in step S2, the vertical correction result is obtained by forward filtering solution of spatiotemporal Kalman filtering; the optimal estimate of the state variable at the current moment is calculated recursively using the state estimate value at the previous moment and the observation value at the current moment.

[0025] Furthermore, the process of calculating the optimal estimate is divided into two stages: prediction and correction:

[0026] Prediction stage:

[0027]

[0028] Correction phase:

[0029]

[0030] Among them, a(t-1|t-1), a(t|t-1) and a(t|t) are the optimal estimate at time t-1, the predicted value and the optimal estimate at time t, respectively; and are the corresponding covariance matrices respectively; K is the Kalman gain matrix, which is obtained based on the covariance minimum principle:

[0031]

[0032] Among them, the covariance matrix

[0033] Furthermore, in step S3, the correction data and the model parameters of the GNSS data are combined, and the fusion result is obtained by forward filtering.

[0034] Furthermore, the frequency of SBAS-InSAR data is 12 days, while the frequency of the selected GNSS data is 3 days. Therefore, there are data null values ​​in the process of model parameter estimation and solution. In the process of parameter estimation, if the value of any group of data is null, the data of this period will not be used in the parameter estimation. In the process of forward filtering calculation using the estimated parameters, if the observation value Z(t) at time t is i is a null value, then Z(t) i Does not participate in the state estimation at time t. Before the state estimation:

[0035] Z(t) i =Ha(t|t-1) i (14)

[0036] Compared with the prior art, the present invention has the following beneficial effects:

[0037] The present invention provides a GNSS and SBAS-InSAR deformation monitoring fusion method based on ISTRE, which is used for fusing GNSS vertical deformation with SBAS-InSAR line of sight deformation, including vertical correction of SBAS-InSAR data, firstly estimating ISTRE model parameters using SBAS-InSAR deformation data; then covering corresponding SBAS-InSAR data with known GNSS data; and then re-weighting the covariance matrix of the observation equation. The vertical correction data is obtained by forward filtering. The correction data is fused with GNSS data, firstly estimating ISTRE model parameters using the correction data and GNSS data; then, fusion calculation is performed by forward filtering. The present invention proposes an ISTRE model, which realizes the fusion calculation of GNSS vertical deformation and SBAS-InSAR line of sight deformation, and has high accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments are briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without creative work.

[0039] Figure 1 It is a flow chart of the GNSS and SBAS-InSAR deformation monitoring fusion method based on ISTRE of the present invention;

[0040] Figure 2 A comparison diagram of vertical correction between the present invention and different spatiotemporal Kalman fusion algorithms;

[0041] Figure 3 This is a comparison chart of the fusion results of the present invention and different spatiotemporal Kalman fusion algorithms at GNSS points.

[0042] Figure 4 This is a comparison chart of the fusion results at the vertical correction point between the present invention and different spatiotemporal Kalman fusion algorithms. DETAILED DESCRIPTION

[0043] The present invention will be further described below in conjunction with embodiments:

[0044] The present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. It is to be understood that the specific embodiments described herein are only used to explain the present invention, rather than to limit the present invention. It should also be noted that, for ease of description, only parts related to the present invention, rather than all structures, are shown in the accompanying drawings.

[0045] It should be noted that similar reference numerals and letters represent similar items in the following drawings, so once an item is defined in one drawing, it does not need to be further defined and explained in the subsequent drawings. At the same time, in the description of the present invention, the terms "first", "second", etc. are only used to distinguish the description and cannot be understood as indicating or implying relative importance.

[0046] The present invention provides a GNSS and SBAS-InSAR deformation monitoring fusion method based on ISTRE, comprising the following steps:

[0047] S1. Perform vertical correction on SBAS-InSAR data. First, convert the SBAS-InSAR cumulative deformation data into deformation velocity data, and then use the SBAS-InSAR deformation velocity data to estimate the ISTRE model parameters; extract GNSS data according to the time interval of SBAS-InSAR data, unify the time scale, and then use the known GNSS data to cover the corresponding SBAS-InSAR data; finally, re-weight the covariance matrix of the observation equation;

[0048] S2, obtaining vertical correction data through forward filtering, and recursively calculating the optimal estimate of the state variable at the current moment using the state estimate at the previous moment and the observed value at the current moment;

[0049] S3: The correction data is fused with the GNSS data. The ISTRE model parameters are estimated using the correction data and GNSS data. Then, the fusion calculation is performed through forward filtering.

[0050] Furthermore, in step S1, the ISTRE model achieves model dimensionality reduction and models large-scale local variations and small-scale subtle variations, taking into account more sophisticated spatiotemporal deformations. Different from traditional spatiotemporal Kalman filter modeling, the ISTRE model needs to first use SBAS-InSAR deformation data to estimate the ISTRE model parameters. The mathematical model of the ISTRE model is:

[0051]

[0052] Among them, Z t (s) is the observation value at time t at position s, a(t) is the state vector at time t, ε is the observation noise, ε~N(0, R), η is the state noise, η~N(0, Q), is the state transfer matrix, ξ(s) is a small-scale variation with spatial irrelevance, μ t (s) is the trend term, and S(s) is the spatial basis describing the spatial correlation.

[0053] The trend term needs to be calculated and eliminated by the least squares method before the ISTRE model is built. The spatial basis function S(s) of the model has no orthogonality restriction. The appropriate spatial basis function can be selected based on experience, such as the bisquare spatial basis, and the spatial basis function can be adjusted until the fitting residual V reaches a stable state.

[0054] Specifically, the observation data input to the model is the deformation velocity data. Since the settlement process in the detection area is stable, the trend field can be identified as μ t (s) = 1, the spatial basis function S(s) of the model has no orthogonality restriction, and a suitable spatial basis function can be selected based on experience, such as the bisquare spatial basis, and the spatial basis function is adjusted until the fitting residual V reaches a stable state:

[0055] V=Z t (s)-μ t (s)-S(s)(S(s)′S(s)) -1 S(s)′(Z t (s)-μ t (s)) (2)

[0056] The ISTRE model establishes small-scale subtle variations ξ(s) separately in the observation equation. ξ(s) is a spatially uncorrelated Gaussian white noise, and its covariance matrix can be calculated using the spherical model and other variance functions. After obtaining the spatial field H, the unknown parameters in the model are obtained based on the maximum likelihood estimation method. Best estimate:

[0057]

[0058] The state vector α(t) can be iteratively calculated using the EM (expectation-maximum) algorithm. After the model parameters are determined, the corresponding SBAS-InSAR data is overlaid with known GNSS data:

[0059] Z(t)=[Z G，1|t … Z G，j|t Z S,k|t … Z S,N|t ] (4)

[0060] Among them, Z G is the vertical velocity data of GNSS, Z S is the velocity data in the SBAS-InSAR LOS direction, j is the number of GNSS points, and k=j+1. The spacing between GNSS deformation monitoring points is small, and the distance between them and the base station is less than 1km. It is assumed that the monitoring accuracy of each point is the same. SBAS-InSAR monitoring also assumes that the accuracy of each pixel is the same, and the deformation results obtained by the two technologies are highly correlated. Therefore, when the GNSS data of the same site is used to replace the SBAS-InSAR data, the spatial field S(s) remains unchanged.

[0061] When the observation data changes, the covariance matrix of the observation equation needs to be re-weighted to ensure the accuracy of the model. The difference between the SBAS-InSAR LOS directional deformation and the vertical correction deformation can be divided into systematic error and noise error. The noise error is corrected by formula (5). The systematic error is based on the concept of robustness to weight the covariance matrix of the observation noise to reduce its impact on the fusion:

[0062]

[0063] Among them, γ G and γ S GNSS and SBAS-InSAR are the variance inflation factors, respectively. The calculation formula of γ is:

[0064]

[0065] Among them, k 0 and k 1 are the quantile parameter and elimination point respectively. ii is the robustness factor:

[0066]

[0067] in, is the i-th element of the residual vector of the observation vector at time k, is the covariance matrix of the residual vector. 0 and k 1The k value selected for GNSS and SBAS-InSAR data with different accuracy is determined according to the size of the standardized residual. 0 and k 1 different.

[0068] Furthermore, in step S2, the vertical correction result is obtained by forward filtering of the spatiotemporal Kalman filter. The optimal estimate of the state variable at the current moment is calculated recursively using the state estimate at the previous moment and the observed value at the current moment. The process is divided into two stages: prediction and correction:

[0069] Prediction stage:

[0070]

[0071] Correction phase:

[0072]

[0073] Among them, a(t-1|t-1), a(t|t-1) and a(t|t) are the optimal estimate at time t-1, the predicted value and the optimal estimate at time t, respectively; and are the corresponding covariance matrices respectively; K is the Kalman gain matrix, which is obtained based on the covariance minimum principle:

[0074]

[0075] Among them, the covariance matrix

[0076] Furthermore, in step S3, the correction data and GNSS data are used as observation data, and the model parameters are obtained using formula (1) to formula (3), and then the fusion result is obtained using forward filtering.

[0077] The frequency of SBAS-InSAR data is 12 days, while the frequency of the selected GNSS data is 3 days. Therefore, there are data null values ​​in the process of model parameter estimation and solution. During the parameter estimation process, if the value of any group of data is null, the data of that period will not be used in the parameter estimation. In the process of forward filtering calculation using the estimated parameters, if the observation value Z(t) at time t is i is a null value, then Z(t) i Does not participate in the state estimation at time t. Before the state estimation:

[0078] Z(t) i =Ha(t|t-1) i (14)

[0079] Example 1

[0080] like Figure 1As shown, a GNSS and SBAS-InSAR deformation monitoring fusion method based on ISTRE includes the following steps:

[0081] S1. Perform vertical correction on SBAS-InSAR data. First, convert the SBAS-InSAR cumulative deformation data into deformation velocity data, and then use the SBAS-InSAR deformation velocity data to estimate the ISTRE model parameters; extract GNSS data according to the time interval of SBAS-InSAR data, unify the time scale, and then use the known GNSS data to cover the corresponding SBAS-InSAR data; finally, re-weight the covariance matrix of the observation equation;

[0082] S2, obtaining vertical correction data through forward filtering, and recursively calculating the optimal estimate of the state variable at the current moment using the state estimate at the previous moment and the observed value at the current moment;

[0083] S3: The correction data is fused with the GNSS data. The ISTRE model parameters are estimated using the correction data and GNSS data. Then, the fusion calculation is performed through forward filtering.

[0084] In order to demonstrate the effect of the present invention, a fusion experiment was conducted using the descending image data provided by the Sentinel-1A satellite and the GNSS data provided by five permanent GNSS monitoring stations, and a comparison was made using the new fusion method at this stage. The results are shown in Figure 2 , Figure 3 and Figure 4 . Figure 2 The first row is the vertical correction image of the accumulated deformation, and the second row represents the vertical correction map of the deformation speed. Figure 3 The first row is the fusion image of the cumulative deformation of GNSS points, and the second row represents the fusion image of the deformation velocity of GNSS points. Figure 4 The first row is the fusion image of the accumulated deformation of the vertical correction points, and the second row represents the fusion image of the deformation speed of the vertical correction points. The quantitative evaluation of vertical correction and fusion is shown in Tables 1 and 2. RMSE is the root mean square error between the GNSS data and AIR is the accuracy improvement rate of the model. The results once again show that the present invention has certain advantages.

[0085] Table 1 Quantitative evaluation of vertical correction experiments of the present invention and the comparison algorithm

[0086]

[0087]

[0088] Table 2 Quantitative evaluation of vertical correction experiments of the present invention and the comparison algorithm

[0089]

[0090] Note that the above are only preferred embodiments of the present invention and the technical principles used. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and that various obvious changes, readjustments and substitutions can be made by those skilled in the art without departing from the scope of protection of the present invention. Therefore, although the present invention has been described in more detail through the above embodiments, the present invention is not limited to the above embodiments, and may include more other equivalent embodiments without departing from the concept of the present invention, and the scope of the present invention is determined by the scope of the appended claims.

Claims

1. A GNSS and SBAS-InSAR deformation monitoring fusion method based on ISTRE, characterized in that: The following steps are involved: S1. Perform vertical correction on SBAS-InSAR data. First, use SBAS-InSAR deformation data to estimate ISTRE model parameters. Then, overlay the corresponding SBAS-InSAR data with known GNSS data. Then, re-weight the covariance matrix of the observation equation. S2, obtaining vertical correction data through forward filtering, and recursively calculating the optimal estimate of the state variable at the current moment using the state estimate at the previous moment and the observed value at the current moment; S3, the correction data is integrated with the GNSS data, and the ISTRE model parameters are estimated using the correction data and GNSS data; Then, the fusion calculation is performed through forward filtering.

2. The ISTRE-based GNSS and SBAS-InSAR deformation monitoring fusion method according to claim 1, characterized in that: In step S1, the model parameters are estimated based on the ISTRE model using the SBAS-InSAR deformation data as observation data. After the model parameters are determined, the corresponding SBAS-InSAR data are overlaid with the known GNSS data: Z(t)=[Z G,1|t … Z G,j|t Z S,k|t … Z S,N|t ] (4) Among them, Z G is the vertical velocity data of GNSS, Z S is the velocity data in the SBAS-InSAR LOS direction, j is the number of GNSS points, k=j+1; the spacing between GNSS deformation monitoring points is small, and the distance between them and the base station is less than 1km, so it is assumed that the monitoring accuracy of each point is the same; SBAS-InSAR monitoring also assumes that the accuracy of each pixel is the same, and the deformation results obtained by the two technologies are highly correlated; therefore, when the GNSS data of the same site is used to replace the SBAS-InSAR data, the spatial field H remains unchanged.

3. The GNSS and SBAS-InSAR deformation monitoring fusion method based on ISTRE according to claim 2 is characterized in that: When the observation data changes, the covariance matrix of the observation equation needs to be re-weighted to ensure the accuracy of the model; the difference between the SBAS-InSAR LOS directional deformation and the vertical correction deformation can be divided into systematic error and noise error. The noise error is corrected by formula (5). The systematic error is based on the concept of robustness to weight the covariance matrix of the observation noise to reduce its impact on the fusion: Among them, γ G and γ S The variance inflation factors for GNSS and SBAS-InSAR are respectively, which can be calculated based on the robustness principle.

4. The ISTRE-based GNSS and SBAS-InSAR deformation monitoring fusion method according to claim 3, characterized in that: γ G and γ S GNSS and SBAS-InSAR are the variance inflation factors, respectively. The calculation formula of γ is: Among them, k0 and k1 are the quantile parameter and elimination point respectively; s ii is the robustness factor: in, is the i-th element of the residual vector of the observation vector at time k, is the covariance matrix of the residual vector; k0 and k1 are determined according to the size of the standardized residual. Different GNSS and SBAS-InSAR data with different accuracies have different k0 and k1.

5. The ISTRE-based GNSS and SBAS-InSAR deformation monitoring fusion method according to claim 1, characterized in that: In step S2, the vertical correction result is obtained by forward filtering solution of spatiotemporal Kalman filtering; the optimal estimate of the state variable at the current moment is calculated recursively using the state estimate value at the previous moment and the observation value at the current moment.

6. The ISTRE-based GNSS and SBAS-InSAR deformation monitoring fusion method according to claim 5, characterized in that: The process of calculating the optimal estimate is divided into two stages: prediction and correction: Prediction stage: Correction phase: a(t|t)=a(t|t-1)+K(Z(t)-Ha(t|t-1)) (11) Among them, a(t-1|t-1), a(t|t-1) and a(t|t) are the optimal estimate at time t-1, the predicted value and the optimal estimate at time t, respectively; and are the corresponding covariance matrices respectively; K is the Kalman gain matrix, which is obtained based on the covariance minimum principle: Among them, the covariance matrix 7. The GNSS and SBAS-InSAR deformation monitoring fusion method based on ISTRE according to claim 1, characterized in that: In step S3, the correction data and GNSS data are combined to estimate the parameters of the ISTRE model, and then the forward filtering solution is used to obtain the fusion result.

8. The ISTRE-based GNSS and SBAS-InSAR deformation monitoring fusion method according to claim 7, characterized in that: The frequency of SBAS-InSAR data is 12 days, while the frequency of the selected GNSS data is 3 days. Therefore, there are data null values ​​in the process of model parameter estimation and solution. In the process of parameter estimation, if the value of any group of data is null, the data of this period will not be used in the parameter estimation. In the process of forward filtering calculation using the estimated parameters, if the observation value Z(t) at time t is i is a null value, then Z(t) i Does not participate in the state estimation at time t. Before the state estimation: Z(t) i =Ha(t|t-1) i (14)。

Images