Mining area earth surface deformation monitoring method by fusing space-time Kalman model with InSAR (Interferometric Synthetic Aperture Radar) and GNSS data

Abstract

The invention discloses a mining area earth surface deformation monitoring method based on fusion of a space-time Kalman model and InSAR and GNSS data. For a mining area with non-uniform deformation space change, large gradient and highly nonlinear time change, the advantages of high spatial resolution of InSAR deformation and high temporal resolution of GNSS deformation are fully utilized. The spatial model is adaptively constructed according to the spatial features of the mining area deformation, and the automation and efficiency of the algorithm are improved. Meanwhile, the time feature constraint of the prior deformation of the mining area is added, so that the dependence of a traditional STRE fusion method on GNSS data is reduced, and the deformation fusion precision is improved. According to the method, the InSAR deformation monitoring result, the GNSS deformation monitoring result and the prior deformation characteristics of the mining area can be considered at the same time, the three kinds of information are effectively fused to the maximum degree, and high-temporal-spatial-resolution deformation monitoring of the mining area is achieved. Under the background of exuberant demand on mining area resources in China, the method is of great significance to accurate monitoring, assessment and prevention of mining area geological disasters.

Description

Technical Field

[0001] The present invention relates to Synthetic Aperture Radar Interferometry (InSAR) technology and Global Navigation Satellite System (GNSS) technology in the field of remote sensing imagery, and in particular to a method for monitoring surface deformation in mining areas by fusing InSAR and GNSS data with a spatiotemporal Kalman model. Background Art

[0002] High-intensity coal mining has brought great impacts on the surface, geological structure and surrounding environment of the mining area, and even caused a series of geological disasters such as landslides, mudslides and collapses, which seriously affected the safety of life and property of local people. Deformation is the most intuitive precursor of geological disasters. Therefore, it is extremely important to use scientific methods to monitor the surface deformation of the goaf for accurate assessment and prevention of geological disasters in mining areas. In the past 20 years, the synthetic aperture radar interferometry (InSAR) technology has made great breakthroughs and achievements in various fields such as mining areas, landslides, earthquakes, atmosphere, frozen soil and urban infrastructure with its unique advantages of high spatial resolution, high precision and large-scale rapid monitoring, and has become one of the most effective means of geological disaster monitoring. However, the temporal resolution of InSAR deformation monitoring is limited by the revisit period of SAR satellite sensors, which is generally 6-14 days, resulting in a large lag in the deformation monitored. In the face of sudden deformations with large deformation gradients such as mining areas and landslides, it is impossible to reflect the deformation situation in time, which greatly limits the development and application of InSAR technology. The Global Navigation Satellite System (GNSS) has the advantages of real-time, high-precision and long-term continuous observation, and has been widely used in safety monitoring of various geological disasters and engineering applications. However, the cost of GNSS is relatively high, and it can only monitor at individual important locations, and cannot be deployed over a large area with high density. Therefore, the spatial resolution of its monitoring results is low, and it is difficult to reflect the spatial changes in deformation in the monitored area, especially in mining areas with large deformation gradient changes.

[0003] Spatio-temporal Kalman filtering is a Kalman filtering method that fully considers the spatio-temporal correlation. It has shown great superiority in the field of spatio-temporal deformation analysis. It can perform spatio-temporal dynamic filtering, spatio-temporal interpolation and spatio-temporal prediction on deformation monitoring data. In particular, the spatio-temporal random effects model (STRE) in spatio-temporal Kalman can also reduce the dimension of the data, greatly improving the computational efficiency of spatio-temporal big data. Therefore, some scholars have proposed combining the advantages of high spatial resolution of InSAR deformation and high temporal resolution of GNSS deformation, using the spatio-temporal Kalman filtering method to fuse the monitoring results of InSAR and GNSS to reconstruct the deformation of the surface with high spatio-temporal resolution, and applying it to the restoration of fault deformation in Southern California, etc., and achieved good results. However, there are still two problems in these studies: First, in the process of using the STRE model to reduce the dimension by laying out the spatial basis in the fusion process, the purpose of fast solution is achieved, and the fewer the spatial basis, the better the dimensionality reduction effect. At present, the layout of spatial bases is mainly uniformly laid out at multiple scales, and relies on human experience guidance. Although it can meet the use requirements in uniform deformation such as faults, it is not suitable for mining areas with uneven deformation and large gradients. Second, when InSAR and GNSS are fused to generate surface deformation results with high spatiotemporal resolution, the deformation accuracy of STRE model fusion is closely related to the number of GNSS. The more GNSS, the higher the fusion deformation accuracy, and the fewer GNSS, the fusion results are mainly dominated by InSAR deformation results. However, it is usually impossible to deploy a large number of GNSS stations in mining areas, resulting in the existing spatiotemporal Kalman fusion method being unsuitable for mining areas with complex spatiotemporal deformation changes. Therefore, in the face of the current demand for deformation monitoring in mining areas, it is urgent to consider how to optimize the spatial modeling method of mining area deformation during spatiotemporal Kalman fusion and solve the problem of low fusion deformation accuracy caused by the small number of GNSS stations in mining areas. Summary of the invention

[0004] The purpose of the present invention is to overcome the defect that the current spatiotemporal Kalman method based on the fusion of InSAR and GNSS deformation to reconstruct the deformation process of the surface with high spatiotemporal resolution is difficult to be applied to the large deformation gradient and small number of GNSS in the mining area, and proposes a spatiotemporal Kalman model to fuse InSAR and GNSS data to reconstruct the surface deformation monitoring method of the mining area. The method adaptively constructs a spatial model according to the spatial deformation characteristics of the mining area, thereby improving the automation and efficiency of the algorithm. At the same time, the time characteristic constraints of the prior deformation of the mining area are added, which not only reduces the dependence of the traditional STRE fusion method on GNSS data, but also further improves the accuracy of deformation fusion.

[0005] The present invention first provides the following technical solution:

[0006] A method for monitoring surface deformation in mining areas by fusing InSAR and GNSS data with a spatiotemporal Kalman model includes the following steps:

[0007] S1. InSAR and GNSS data processing: The acquired InSAR interferograms are filtered, unwrapped, subjected to large-scale trend removal, atmospheric correction, and space-time filtering before being geocoded. The average deformation rate and time series deformation results of the study area are then calculated based on the SBAS-InSAR technology. The GNSS three-dimensional deformation monitoring results are converted to the InSAR line of sight. Finally, the GNSS deformation and InSAR deformation are unified in terms of time and space.

[0008] S2. Establish a spatiotemporal Kalman model that integrates InSAR and GNSS deformations. First, establish the STRE model, and then establish a fused spatial model and temporal model based on the acquired InSAR and GNSS data to form a spatiotemporal Kalman fusion model.

[0009] S3. Construct an adaptive spatial model for large-scale major deformations based on the spatial deformation characteristics of the mining area: Based on the average deformation rate of InSAR, first lay out a uniform spatial basis in the entire map to capture global smooth deformations; then calculate the residual after the first-layer spatial basis modeling, identify the part with a larger median value of the residual, and perform a morphological dilation operation on it, and lay out the second-layer spatial basis for the obtained irregular area to capture local deformations. Then calculate the residual after the first and second-layer total spatial basis modeling, identify the part with a larger median value of the residual, and perform a morphological dilation operation on it, and lay out the third-layer spatial basis for the obtained irregular area to capture local small deformations. Finally, merge the three-layer InSAR spatial basis and the GNSS spatial basis to complete the adaptive spatial modeling of large-scale major deformations.

[0010] S4. According to the temporal deformation characteristics of the mining area, the prior information of the mining area is introduced, and the Logistic model constraints are added to the state quantity to improve the estimation accuracy of the state quantity without InSAR observation time, while reducing the dependence on GNSS data. Finally, a spatiotemporal Kalman model of fusion InSAR and GNSS deformation with spatiotemporal characteristics of deformation is constructed.

[0011] S5. Solve the model parameters in the fusion process: First, determine the variance of InSAR small-scale subtle deformation, the variance of InSAR observation noise, and the variance of GNSS observation noise through semivariogram fitting and variance component estimation; then determine the model parameters in the fusion process through least squares fitting, forward filtering, and backward smoothing; finally, use EM iterative estimation for the entire fusion process to improve the estimation accuracy of model parameters.

[0012] S6. Reconstruct the deformation results of the surface with high temporal and spatial resolution: Reconstruct and restore the deformation results of the surface with high temporal and spatial resolution based on the spatial basis, state quantities and small-scale subtle deformations solved in the above steps, and evaluate the accuracy of the deformation results.

[0013] According to some preferred embodiments of the present invention, S1 comprises:

[0014] S11. The acquired radar images are geocoded to form an interferogram after registration, small baseline meshing, interference, flat land and terrain removal, phase filtering, phase unwrapping, large-scale trend removal, atmospheric correction and spatiotemporal filtering. It is worth noting that since the trend in deformation is not modeled later, the trend item in the interferogram must be removed in this step.

[0015] S12. Assuming that there are N SAR images in the study area, M interferograms are generated according to the temporal and spatial baseline threshold, and there is only one small baseline set, the following equation can be constructed by the short baseline set (Small BAseline Subsets, SBAS) method:

[0016]

[0017] By least squares, we can get the above formula:

[0018]

[0019] in is the final deformation sequence result of the adjustment; A is the coefficient matrix in the above formula; P is the weight matrix, which can be determined according to the coherence of the interference pattern; L is the phase of the interference pattern.

[0020] After the deformation sequence is obtained, the least squares fitting can be used to obtain the average deformation rate of the entire time period.

[0021] S13. Since GNSS has deformation in three directions, and InSAR can generally only obtain deformation in the radar line of sight (LOS), it is necessary to convert the 3D direction of GNSS into the radar line of sight:

[0022]

[0023] in Denotes the transformation of GNSS 3D deformation to LOS backward deformation, D E , D N and D U They represent the east-west, north-south and vertical deformations of the GNSS point respectively, θ represents the local incident angle of the satellite radar, and α represents the radar azimuth.

[0024] S14. Convert the InSAR points and GNSS points to the same coordinate frame.

[0025] In space, find the InSAR point closest to the GNSS point as the GNSS point of the same name.

[0026] Then, one of the stable GNSS points is selected to correct the time and space dimensions of the remaining GNSS points and InSAR, so as to achieve the unification of the time and space benchmarks of the two.

[0027] According to some preferred embodiments of the present invention, S2 comprises:

[0028] S21. For a series of spatiotemporal deformation data Z(s, t), it can be expressed as:

[0029] Z(s, t) = Y(s, t) + ω(s, t)

[0030] Where s = 1, 2, ..., n and t = 1, 2, ..., T represent the spatial position and observation time of the midpoint in region D, respectively, Y(s, t) represents the true value of the spatiotemporal observation data, ω(s, t) ~ N(0, σ n 2 ) represents the noise of spatiotemporal observation data.

[0031] The true value Y(s, t) of the further spatiotemporal observation data can be expressed as follows:

[0032] Y(s, t) = v t (s)+ξ t (s)

[0033] where v t (s) represents the main large-scale deformation in space, ξ t (s) represents a small-scale deformation in space. Note that if the deformation contains a trend term, the trend term can also be spatially modeled. Since the trend is generally removed in advance in mining area deformation monitoring, there is no trend term deformation.

[0034] Further major deformations on large scales in space v t (s) is generally modeled as a spatial field S t (s) and the time-varying state quantity x t It can be described by a linear combination of:

[0035]

[0036] Where S t (s)=[S 1，t (s), S 2，t (s),...,S r,t (s)] T, represents the established r×n dimensional space basis function, x t =[x 1，t (·), x 2，t (·),...,x r,t (·)] T Represents a zero-mean Gaussian random state of r×m. Note that the basis function can change over time, but for a mining area, the deformation area can be easily captured by the average deformation rate. In order to reduce the complexity of the model, it is modeled as a function that does not change over time, that is, S(s) = S t (s)=[S1(s),S2(s),...,S r (s)] T .

[0037] The basis function is the bisquare function commonly used in space-time analysis, which is expressed as follows:

[0038]

[0039] Among them, S l，j (s) represents the spatial basis function value of the jth point in the lth spatial basis on the space s, j = 1, 2, ..., r; p l,j represents the spatial basis position of the jth point in the lth spatial basis; d l It is the range of spatial variation in the lth level of space (generally 1.5 times the shortest spatial base distance at the lth level of scale).

[0040] S22, now the state quantity x t Evolution in time constitutes the state transfer equation:

[0041] x t =Hx t-1 +e t

[0042] Where H is the r×r state transition matrix, is the noise in the state transition process. The final form of the STRE model can be established through the above:

[0043]

[0044] From the above analysis, we can see that the model parameters of STRE are small-scale subtle deformations ξ t (s), state quantity x t , the state transfer matrix H and the variance in the state transfer process

[0045] S23. Assume that in the study area D, a total of n1 InSAR points are obtained, with m1 observation dates T1; and a total of n2 GNSS points are obtained, with m2 observation dates Assuming that all the time to be reconstructed is T=T2, and the time correspondence between InSAR and GNSS is Id, then T1=T2(Id)=T(Id). Let n=n1+n2, m=m2.

[0046] According to the spatiotemporal random effect model of InSAR and GNSS fusion established by STRE in S21 above:

[0047]

[0048] The meanings of each variable are as follows in Table 1:

[0049] Table 1 Meaning of variables

[0050]

[0051] Note that the variances of small-scale subtle deformations are the same for InSAR and GNSS during the fusion process.

[0052] Finally, the same as step S22, the state quantity x t The evolution in time constitutes the state transfer equation, forming the final fused space-time Kalman model.

[0053] According to some preferred embodiments of the present invention, S3 includes:

[0054] S31, in step S23, in order to fully capture the spatial variation of deformation, multiple layers of spatial bases are generally arranged at multiple scales, which is equivalent to the fitting of multiple variograms. However, the deformation does not change uniformly in space, so it is necessary to make some adaptive adjustments to the position of the spatial base.

[0055] First, a uniform spatial basis is laid out in the entire deformation field according to the InSAR average deformation rate, and then the residual deformation remaining in this uniform spatial basis model is calculated:

[0056]

[0057] Among them, V1 represents the residual after the first layer of spatial basis modeling is deployed, and Z a represents the average deformation rate of InSAR, and S1 represents the r1 spatial basis arranged in the first layer.

[0058] S32, identify the part with larger value in V1. In this process, some isolated islands and small areas will appear. Therefore, it is necessary to delete the small area first, and perform morphological expansion operation at the same time to fill the holes and expand the boundaries to obtain the area where the first layer of spatial basis cannot fully capture the deformation.

[0059] S33. For the obtained irregular area, a second layer of uniform spatial basis (assuming r2) is arranged with a higher density, and combined with the first spatial basis to calculate a new residual:

[0060]

[0061] Use S 1,2 represents the r1+r2 space basis of the first and second layouts, then S 1,2 =[S1,S2] T .

[0062] S34. Since the deformation gradient of the mining area is relatively large, the two-layer space basis still cannot fully capture all deformations, so it is necessary to lay out a third-layer space basis with a smaller spacing on this basis. Similar to step S32, identify the part with a larger value in V2, and some isolated islands and smaller areas will appear in the process. Therefore, it is necessary to delete the small area first, and perform a morphological expansion operation at the same time, fill the holes and expand the boundaries to a certain extent, and obtain the areas that the first and second layers of space basis cannot fully capture.

[0063] S35. Arrange the third layer of uniform spatial basis for the obtained irregular area, and combine it with the previous two spatial basis to calculate the new residual:

[0064]

[0065] Let S3 represent the r3 space basis of the third arrangement, then S ins =[S1,S2,S3] T , there are a total of r=r1+r2+r3 space bases.

[0066] In the above process, you only need to specify the density of each layer of spatial basis to obtain the position and number of spatial basis for the entire area. This not only reduces human intervention and improves the automation of the algorithm, but also reduces the number of spatial basis and improves the efficiency of the algorithm. Especially when performing large-scale deformation calculations, it can greatly save time and space costs.

[0067] It is worth noting that in general, after three layouts, the residual will be very small, which can meet the needs of subsequent calculations. However, when the deformation gradient is very large, the three-scale spatial basis may still not be able to fully capture all deformations. In this case, more scale spatial basis can be laid out on the basis of these three layers to further improve the accuracy of spatial modeling. However, too many spatial basis will consume a lot of computing time and may cause overfitting, so a compromise is needed.

[0068] S36. Use all GNSS points to calculate the GNSS spatial basis S gnss So far, the spatial bases of InSAR and GNSS have been obtained respectively.

[0069] According to some preferred embodiments of the present invention, S4 includes:

[0070] S41. Since the deformation of the mining area has obvious time characteristics, it generally conforms to the Logistic model. Therefore, the nonlinear least squares fitting of the following three-parameter Logistic model can be used for each InSAR time series:

[0071]

[0072] Where F(t ins ) indicates t ins is the InSAR shape variable at time , and a, b, c are the model parameters to be estimated.

[0073] After obtaining the parameters, the predicted deformation of the model can be obtained for any time t:

[0074]

[0075] in are the parameters estimated by nonlinear least squares.

[0076] S42. At this time, the relationship between the model prediction value and the state quantity can be established by the observation equation in step S22:

[0077]

[0078] In general, the magnitude of subtle deformation is relatively small, so when adding constraint equations, small-scale subtle deformations can be ignored, and the constraint conditions of the state quantity can be obtained:

[0079] Gx t =F s (t)

[0080] in Represents the spatial basis for evaluation.

[0081] Ordinary STRE does not contain additional information about the system, but for deformations such as mining areas, the state quantities should satisfy these constraints, so this additional information is used to obtain better filtering performance than the Kalman filter.

[0082] According to the above, we can get the STRE model with state constraints, referred to as cSTRE:

[0083]

[0084] where F s (t,θ) represents the established spatially uncorrelated and temporally correlated function constraint model, and θ represents the model parameters.

[0085] S43. The original STRE model was built based on a linear system. However, the deformation of the mining area has highly nonlinear characteristics, and its state quantity also exhibits highly nonlinear characteristics. Therefore, there are errors beyond the system in the original STRE model, that is, the complete system description is different from the description established by the standard space-time Kalman filter. At this time, the additional state constraints can improve the performance of the original STRE.

[0086] Kalman filtering is a state optimal estimation method. After adding model constraints, this optimal estimation characteristic should be maintained, so it is necessary to solve the optimal constraint state.

[0087] A lot of research has been done on the Kalman filtering problem with additional prior constraints. Generally, the solution of standard Kalman filtering is converted into a combination of standard Kalman filtering and quadratic programming problem, and good performance has been achieved.

[0088] use Represents the constraint state quantity at time t, and represents the optimal state quantity of the unconstrained Kalman estimate, and is used Unconstrained Kalman estimates the posterior variance. According to existing research, one solution is to project the unconstrained state estimate onto the constraint surface. Therefore, the following objective function is established based on the minimum variance criterion:

[0089]

[0090] where W is a positive definite weight matrix, which can be set as Or the identity matrix I.

[0091] Then the optimal solution of the constrained Kalman can be expressed as:

[0092]

[0093] Since an uncertain model is introduced, it is also necessary to re-estimate the variance of the state quantity with model constraints. According to the error propagation rate, it can be obtained:

[0094]

[0095] Where J = WG'(GWG') -1 .

[0096] The above constraints indicate that when the variance of the state quantity evaluated by Kalman is small, it can be considered more reliable, and a smaller weight is given to the model constraint at this time; when the variance of the state quantity evaluated by Kalman is large, it is considered that it needs to be corrected, and a larger weight is given to the model constraint at this time. In this way, it can be ensured that the Kalman filter result does not deviate from the model fitting result. When designing, it should be noted that the state quantity close to the GNSS point should be kept unchanged as much as possible to avoid damaging the accuracy of GNSS.

[0097] At this point, we have obtained an estimate of the state quantity and its variance of the additional constraints. When performing the next round of Kalman filtering, we only need to use replace use replace Then perform a new iteration.

[0098] Through the above process, the spatiotemporal Kalman model of InSAR and GNSS fusion constrained by deformation spatiotemporal characteristics can be obtained, that is, the construction process of the cSTRE model.

[0099] According to some preferred embodiments of the present invention, S5 includes:

[0100] S51. The residual information after InSAR spatial modeling contains small-scale subtle deformation and noise, and the variance of these two types of information needs to be separated. The residual after InSAR deformation spatial modeling of each scene is fitted with a semivariogram to obtain the variance value at different distances:

[0101]

[0102] Where h is the distance between deformation points; N(h) is the number of pairs of all observation points with h as the distance; z(x i ) and z(x i +h) represent the InSAR deformation observation values ​​at the relative distance h; is the variance with a spacing of h, i.e., the variation value, which increases with the increase of h within a certain range. When the distance between the measured points is greater than the maximum correlation distance, the value tends to be stable. Through this step, the variance of the InSAR residual deformation can be obtained. with sample values ​​at spacing h.

[0103] S52. Use the classic spherical function model to fit the InSAR data samples at each time to obtain model parameters such as nugget, range, sill and partial sill values. The nugget reflects the variance of the InSAR spatially uncorrelated noise ω(s,t) The base reflects the small-scale subtle deformation related to space. t Variance of (s) Through the above operations, the variance of InSAR spatially correlated small-scale subtle deformation and the variance of observation noise are separated.

[0104] S53. Since the observation accuracy of InSAR and GNSS is different, the accuracy ratio k of the two observation methods can be determined by using the variance component estimation based on the variance of the InSAR small-scale subtle deformation and the variance of the observation noise obtained in step S52, and then the observation noise of GNSS is calculated according to the following formula:

[0105]

[0106] The GNSS observation noise can be obtained as:

[0107]

[0108] S54. From the above analysis, we can know that the unknown parameters of the cSTRE model are: small-scale subtle deformation ξ t (s), the state quantity x of additional constraints t , the state transfer matrix H and the variance in the state transfer process

[0109] Before solving, the remaining parameters can be set to initial empirical values, and the initial value of the state transfer matrix H is generally set to:

[0110] H=ρE

[0111] where 0<ρ<1 is the multiplication factor and E is the identity matrix.

[0112] The parameters obtained in the above steps are brought into the fusion model for solution. The fusion model is generally solved by forward filtering and backward smoothing. That is, the following are calculated in sequence: state quantity prediction, state quantity prior variance estimation, calculation of spatiotemporal Kalman gain, posterior optimal estimation of state quantity and its variance, and constraint state quantity and its variance estimation. Then repeat the above process for all times, and finally use EM iteration for the entire fusion process to further improve the accuracy of parameter estimation.

[0113] According to some preferred embodiments of the present invention, S6 includes:

[0114] S61. Evaluate and obtain model parameters according to the above steps to reconstruct the high temporal and spatial resolution deformation of the surface:

[0115]

[0116] Among them, s represents the position of the spatial point after reconstruction, which is generally consistent with the distribution of InSAR points, and t represents the time of the reconstruction point, which is generally consistent with the GNSS observation time.

[0117] S62. The accuracy of the reconstructed deformation monitoring results can be verified in time and space respectively.

[0118] In terms of time, a part of the real surface GNSS data that is not involved in the fusion is used to calculate the root mean square error between it and the fusion result.

[0119] In order to facilitate accuracy verification in space, a part of the InSAR data can be deleted before reconstruction so that it does not participate in the fusion. The fusion result is then used to calibrate this part of the InSAR data and calculate its root mean square error.

[0120] Through the above two processes, the reconstruction accuracy of the deformation in time and space of the method proposed in the present invention can be verified respectively.

[0121] According to the above-mentioned method for monitoring surface deformation of mining areas by fusing InSAR and GNSS data with spatiotemporal Kalman model, a deformation monitoring system for reconstructing the surface of mining areas with high spatiotemporal resolution by fusing InSAR and GNSS data with spatiotemporal Kalman model based on prior deformation feature constraints can be obtained.

[0122] The above fusion method or system can be used for deformation monitoring in mining areas and other scenarios.

[0123] The present invention aims at mining areas with uneven spatial deformation changes, large gradients, and highly nonlinear temporal changes. The proposed method of reconstructing mining area surface deformation monitoring by fusing InSAR and GNSS data with a spatiotemporal Kalman model makes full use of the advantages of high spatial resolution of InSAR technology and high temporal resolution of GNSS technology. The spatial model is adaptively constructed according to the spatial characteristics of the mining area deformation, which improves the automation and efficiency of the algorithm. At the same time, the temporal characteristic constraints of the prior deformation of the mining area are added, which not only reduces the dependence of the traditional STRE fusion method on GNSS data, but also further improves the accuracy of deformation fusion. This method can take into account the InSAR deformation monitoring results, GNSS deformation monitoring results and prior deformation characteristics of the mining area at the same time, effectively fuse these three types of information to the greatest extent, and realize deformation monitoring of the mining area with high spatiotemporal resolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0124] Figure 1 A schematic diagram of a specific implementation process of the method of the present invention;

[0125] Figure 2 The InSAR deformation sequence results simulated in the embodiment;

[0126] Figure 3 The distribution result of GNSS in the embodiment;

[0127] Figure 4 yes Figure 2 The accuracy comparison results of 11 GNSS and InSAR;

[0128] Figure 5 It is the InSAR deformation rate in the embodiment and the 419 spatial basis of the traditional three-layer uniform layout;

[0129] Figure 6 In the embodiment Figure 5 residual distribution;

[0130] Figure 7 359 spatial bases of the InSAR deformation rate and the three-layer arrangement of the method of the present invention in the embodiment;

[0131] Figure 8 yes Figure 7 residual distribution;

[0132] Fig. 9 It is the InSAR time series deformation results reconstructed by the traditional method and the method of the present invention and the residual compared with the observed value;

[0133] Fig.10 The comparison results of the deformation monitored by the two GNSS stations DT01 and DT07 of the embodiment and the deformation reconstructed with high temporal and spatial resolution are shown;

[0134] Fig.11 The reconstruction results of the deleted three-scene InSAR deformation of the embodiment are compared with the SBAS-InSAR deformation results. DETAILED DESCRIPTION

[0135] The present invention is described in detail below in conjunction with the embodiments and drawings, but it should be understood that the embodiments and drawings are only used to exemplify the present invention and do not constitute any limitation on the protection scope of the present invention. All reasonable changes and combinations within the scope of the inventive concept of the present invention fall within the protection scope of the present invention.

[0136] This example is carried out using 21 scenes of Sentinel-1SAR data from July 20, 2020 to March 17, 2021 in a certain area. Figure 1 A specific spatiotemporal Kalman model is shown to fuse InSAR and GNSS data to reconstruct surface deformation monitoring in mining areas.

[0137] Figure 2 is the deformation rate obtained by the SBAS-InSAR method (in cm / year), where the red five-pointed star represents the common reference point of InSAR and GNSS, and the temporal and spatial benchmarks of InSAR and GNSS are unified based on this point, and the 11 triangles represent the remaining GNSS stations. The black box is the selected area of ​​interest.

[0138] Figure 3 The 24 SBAS-InSAR images of the area of ​​interest were used to monitor the deformation in real time. Then, the DT01 and DT07 images in GNSS and the three images with InSAR dates of 20201024, 20201211 and 20210128 were selected and not involved in the fusion, and used to verify the accuracy of reconstructing the surface deformation with high temporal and spatial resolution.

[0139] Figure 4 yes Figure 2 From the accuracy comparison results of 11 GNSS and InSAR, we can see that the RMSE of each point is higher than 0.65cm, which fully demonstrates the reliability of SBAS-InSAR deformation monitoring results.

[0140] Figures 5 to 8 It is the comparison result of the spatial base layout between the traditional method and the method of the present invention. Figure 5 It is the InSAR deformation rate and the 419 spatial bases of the traditional three-layer uniform layout. Figure 6 yes Figure 5 The residual distribution has RMS = 1.18 cm. Figure 7 It is the InSAR deformation rate and the 359 spatial bases of the three-layer layout of the method of the present invention. Figure 8 yes Figure 7 The residual distribution has RMS = 1.17 cm. Therefore, the spatial basis arranged by the new method is not only less in number but also more accurate, so it is more suitable for deformation monitoring in mining areas.

[0141] Fig. 9 The InSAR time series deformation results reconstructed by the traditional method and the new method and their residuals compared with the observed values. It can be seen that the traditional method has a larger residual in the deformation area, while the new method has a smaller error due to the introduction of prior deformation features as constraints. Overall, the root mean square of the residual of the new method is reduced by 10.97% compared with the traditional method, which fully demonstrates the superiority of the new method.

[0142] Fig.10The results of the comparison of deformation monitoring and reconstruction with high spatiotemporal resolution at two GNSS stations, DT01 and DT07, are shown. At DT07, the accuracy of STRE and cSTRE reconstructions is comparable due to the high accuracy of InSAR itself. At DT01, although the accuracy of the new method is slightly higher than that of the traditional method, the reconstruction accuracy is not very high due to the influence of InSAR observation errors.

[0143] Fig.11 The comparison between the reconstruction results of the deleted three-scene InSAR deformation and the SBAS-InSAR deformation results is further shown. It can be seen that the deformation error reconstructed by the traditional method is larger, and the deformation will be underestimated in local areas. The RMSE of the residuals between it and the deformation monitored by SBAS-InSAR are 0.47cm, 0.61cm and 0.63cm respectively. The new method has a very small error due to the introduction of the logistic model. The RMSE of the reconstructed deformation residuals are 0.42cm, 0.51cm and 0.55cm respectively. Compared with the traditional method, the accuracy of the new method has increased by 10.64%, 16.39% and 12.70% respectively, with an average increase of up to 13.24%.

[0144] The above analysis fully demonstrates the effectiveness of the method proposed in the present invention, which can be used to fuse GNSS and InSAR to reconstruct surface deformation with high temporal and spatial resolution in complex deformation situations such as mining areas.

[0145] The above embodiments are only preferred embodiments of the present invention, and the protection scope of the present invention is not limited to the above embodiments. All technical solutions under the concept of the present invention belong to the protection scope of the present invention. It should be pointed out that for ordinary technicians in this technical field, improvements and modifications without departing from the principle of the present invention should also be regarded as the protection scope of the present invention.

Claims

1. A method for monitoring surface deformation in mining areas by fusing InSAR and GNSS data with a spatiotemporal Kalman model, characterized in that: The following steps are involved: S1. InSAR and GNSS data processing: The acquired InSAR interferograms are filtered, unwrapped, trend removed over a large area, atmospheric corrected, and spatially and temporally filtered before being geocoded. The average deformation rate and time series deformation results of the study area are then calculated based on the SBAS-InSAR technology. Convert the GNSS 3D deformation monitoring results to the InSAR line of sight; finally, unify the GNSS deformation and InSAR deformation in time and space; S2. Establish a spatiotemporal Kalman model that integrates InSAR and GNSS deformations. First, establish a STRE model, and then establish a fused spatial model and temporal model based on the acquired InSAR and GNSS data to form a spatiotemporal Kalman fusion model. S3. Construct an adaptive spatial model for large-scale major deformations based on the spatial deformation characteristics of the mining area: Taking the average deformation rate of InSAR as the benchmark, first lay out a uniform spatial basis in the whole map to capture the global smooth deformation; then calculate the residual after the first-layer spatial basis modeling, identify the part with a larger median value of the residual, and perform a morphological dilation operation on it, and lay out the second-layer spatial basis for the obtained irregular area to capture local deformation; then calculate the residual after the total spatial basis modeling of the first and second layers, identify the part with a larger median value of the residual, and perform a morphological dilation operation on it, and lay out the third-layer spatial basis for the obtained irregular area to capture local small deformation; finally, merge the three-layer InSAR spatial basis and the GNSS spatial basis to complete the adaptive spatial modeling of large-scale major deformations; S4. According to the temporal deformation characteristics of the mining area, the prior information of the mining area is introduced, and the Logistic model constraints are added to the state quantity to improve the estimation accuracy of the state quantity without InSAR observation time, while reducing the dependence on GNSS data; finally, a spatiotemporal Kalman model of InSAR and GNSS deformation is constructed with the constraints of the spatiotemporal characteristics of deformation; S5. Solve the model parameters in the fusion process: first, determine the variance of the InSAR small-scale subtle deformation, the variance of the InSAR observation noise, and the variance of the GNSS observation noise through semivariogram fitting and variance component estimation; then determine the model parameters in the fusion process through least squares fitting, forward filtering, and backward smoothing; finally, use EM iterative estimation to improve the estimation accuracy of the model parameters in the entire fusion process; S6. Reconstruct the deformation results of the surface with high temporal and spatial resolution: Reconstruct and restore the deformation results of the surface with high temporal and spatial resolution based on the spatial basis, state quantities and small-scale subtle deformations solved in the above steps, and evaluate the accuracy of the deformation results.

2. The method for monitoring surface deformation in mining areas by fusing InSAR and GNSS data with a spatiotemporal Kalman model according to claim 1 is characterized in that: The S1 comprises the following sub-steps: S11. The acquired radar images are geocoded to form an interferogram after registration, small baseline meshing, interference, flat land and terrain removal, phase filtering, phase unwrapping, large-scale trend removal, atmospheric correction and spatiotemporal filtering. This step is to remove the trend item in the interferogram. S12, let there be N SAR images in the study area, generate M interferograms according to the temporal and spatial baseline threshold, and there is only one small baseline set, then the short baseline set SBAS method is used to construct the following equation: By least squares, we can get the above formula: in is the final deformation sequence result of the adjustment; A is the coefficient matrix in the above formula; P is the weight matrix, which is determined according to the coherence of the interference pattern; L is the phase of the interference pattern; T is the transpose of the matrix; After the deformation sequence is obtained, the average deformation rate of the entire time period is obtained by least square fitting; S13. Convert the 3D direction of GNSS to the radar line of sight: in Denotes the transformation of GNSS 3D deformation to LOS backward deformation, D E , D N and D U They represent the east-west, north-south and vertical deformations of the GNSS point, θ represents the local incident angle of the satellite radar, and α represents the radar azimuth; S14, converting the InSAR points and the GNSS points to the same coordinate frame; In space, find the InSAR point closest to the GNSS point as the GNSS point of the same name; Then, one of the stable GNSS points is selected to correct the time and space dimensions of the remaining GNSS points and InSAR, so as to achieve the unification of the time and space benchmarks of the two.

3. The method for monitoring surface deformation in mining areas by fusing InSAR and GNSS data with a spatiotemporal Kalman model according to claim 1 is characterized in that: The S2 comprises the following sub-steps: S21. For a series of spatiotemporal deformation data Z(s,t), it can be expressed as: Z(s,t)=Y(s,t)+ω(s,t) Where s = 1, 2, ..., n and t = 1, 2, ..., T represent the spatial position and observation time of the midpoint in region D respectively, Y(s, t) represents the true value of the spatiotemporal observation data, ω(s, t) ~ N(0,σ n 2 ) represents the noise of spatiotemporal observation data; The true value Y(s,t) of spatiotemporal observation data is expressed as follows: Y(s,t)=υ t (s)+ξ t (s) where υ t (s) represents the main large-scale deformation in space, ξ t (s) represents a small-scale subtle deformation in space; If the deformation contains a trend term, perform spatial modeling on the trend term; The main large-scale deformation in space t (s) is modeled as a spatial field S t (s) and the time-varying state quantity x t It can be described by a linear combination of: Where S t (s)=[S 1,t (s),S 2,t (s),…,S r,t (s)] T , represents the r×n dimensional space basis function, x t =[x 1,t (·),x 2,t (·),…,x r,t (·)] T Represents r×m zero-mean Gaussian random state quantity; Model it as a function that does not change with time, that is, S(s) = S t (s)=[S1(s),S2(s),…,S r (s)] T . ; The basis function is the bisquare function commonly used in space-time analysis, which is expressed as follows: Among them, S l,j (s) represents the spatial basis function value of the jth point in the lth spatial basis on the space s, j = 1, 2, ..., r; p l,j represents the spatial basis position of the jth point in the lth spatial basis; d l is the range of spatial variation in the l-th layer of space; S22, now the state quantity x t Evolution in time constitutes the state transfer equation: x t =Hx t-1 +e t Where H is the r×r state transition matrix, is the noise in the state transition process; the final form of the STRE model is established: The model parameters of STRE are small-scale subtle deformations ξ t (s), state quantity x t , the variance during the state transition S23, let n1 InSAR points be obtained in the study area D, with m1 observation date T1; and n2 GNSS points be obtained, with m2 observation date ; Let all the moments to be reconstructed be T=T2, and the time correspondence between InSAR and GNSS be Id, then T1=T2(Id)=T(Id). ; Let n=n1+n2, m=m2. ; According to the spatiotemporal random effect model of InSAR and GNSS fusion established by STRE in S21 above: The meanings of the variables are: Z ins (s, t) is the deformation observation value of InSAR; Z gnss (s, t) is the deformation observation value of GNSS; is the spatial basis of InSAR; is the spatial basis of GNSS; ξ ins (s, t) is the small-scale subtle deformation of InSAR; ξ gnss (s, t) is the small-scale subtle deformation of GNSS; ω ins (s, t) is the InSAR observation noise; ω gnss (s, t) is the GNSS observation noise; is the total variance of the small-scale subtle deformation and observation noise of InSAR; is the total variance of the small-scale subtle deformation and observation noise of GNSS; is the variance of small-scale subtle deformation; the variance of small-scale subtle deformation of InSAR and GNSS is the same during the fusion process; is the variance of InSAR observation noise; is the variance of GNSS observation noise; Finally, the same as step S22, the state quantity x t The evolution in time constitutes the state transfer equation, forming the final fused space-time Kalman model.

4. The method for monitoring surface deformation in mining areas by fusing InSAR and GNSS data with a spatiotemporal Kalman model according to claim 1 is characterized in that: The S3 comprises the following sub-steps: S31. Make some adaptive adjustments to the position layout of the space base; First, a uniform spatial basis is laid out in the entire deformation field according to the InSAR average deformation rate, and then the residual deformation remaining in this uniform spatial basis model is calculated: Among them, V1 represents the residual after the first layer of spatial basis modeling is deployed, and Z a represents the average deformation rate of InSAR, S1 represents the r1 spatial basis arranged in the first layer; S32, identify the part with large median value in V1, some isolated islands and small areas will appear; first delete the small area, and perform morphological expansion operation at the same time, fill the holes and expand the boundaries, and obtain the area where the first layer of spatial basis cannot fully capture the deformation; S33. For the obtained irregular area, a second layer of uniform spatial basis is arranged with a higher density, let it be r2, and combined with the first spatial basis to calculate the new residual: Use S 1,2 represents the r1+r2 space basis of the first and second layouts, then S 1,2 =[S1,S2] T . ; S34, laying out the third layer of spatial bases with small spacing; firstly deleting the small area, and performing morphological dilation operation at the same time, filling the holes therein and expanding the boundaries to a certain extent, so as to obtain the areas that the first and second layers of spatial bases cannot fully capture; S35. Arrange the third layer of uniform spatial basis for the obtained irregular area, and combine it with the previous two spatial basis to calculate the new residual: Let S3 represent the r3 space basis of the third arrangement, then S ins =[S1,S2,S3] T , a total of r = r1 + r2 + r3 space bases; S36. Use all GNSS points to calculate the GNSS spatial basis S gnss . ; So far, the spatial bases of InSAR and GNSS have been obtained respectively.

5. The method for monitoring mining area surface deformation by fusing InSAR and GNSS data with a spatiotemporal Kalman model according to claim 1 is characterized in that: The S4 comprises the following sub-steps: S41. For each InSAR time series, a nonlinear least squares fit is used to fit the following three-parameter Logistic model: Where F(t ins ) indicates t ins The InSAR shape variable at time , a, b, c are the model parameters to be estimated; After obtaining the parameters, the predicted deformation of the model is obtained for any time t: in are the parameters obtained by nonlinear least squares estimation; S42. At this time, the relationship between the model prediction value and the state quantity is established by the observation equation in step S22: When adding constraint equations, we ignore small-scale subtle deformations and obtain the constraint conditions of the state quantity: Gx t =F s (t) in represents the spatial basis of the evaluation; The STRE model with state constraints is obtained, referred to as cSTRE: where F s (t,θ) represents the established spatially uncorrelated and temporally correlated function constraint model, and θ represents the model parameters; S43. The original STRE model was built based on a linear system. However, the deformation of the mining area has highly nonlinear characteristics, and its state quantity also exhibits highly nonlinear characteristics. There are errors beyond the system in the original STRE model, that is, the complete system description is different from the description established by the standard space-time Kalman filter. At this time, the additional state constraints improve the performance of the original STRE. use Represents the constraint state quantity at time t, and represents the optimal state quantity of the unconstrained Kalman estimate, and is used Unconstrained Kalman estimates the posterior variance and establishes the following objective function based on the minimum variance criterion: Where W is a positive definite weight matrix, set it to Or the identity matrix I. ; The optimal solution of the constrained Kalman is expressed as: Re-estimate the variance of the state quantity with model constraints, and according to the error propagation rate: Where J = WG'(GWG') -1 . ; At this point, the state quantity and its variance estimate of the additional constraint are obtained. When the next round of Kalman filtering is performed, it is used in the state transition prediction. replace use replace Then a new iteration is performed; Through the above process, the spatiotemporal Kalman model of InSAR and GNSS fusion constrained by deformation spatiotemporal characteristics is obtained, that is, the construction process of the cSTRE model.

6. The method for monitoring surface deformation in mining areas by fusing InSAR and GNSS data with a spatiotemporal Kalman model according to claim 1 is characterized in that: The S5 comprises the following sub-steps: S51. The residual information after InSAR spatial modeling contains small-scale subtle deformation and noise, and the variance of these two types of information is separated; the residual after InSAR deformation spatial modeling of each scene is fitted with a semivariogram to obtain the variance values ​​at different distances: Where h is the distance between deformation points; N(h) is the number of pairs of all observation points with h as the distance; z(x i ) and z(x i +h) represent the InSAR deformation observation values ​​at the relative distance h; is the variance with a spacing of h, i.e., the variation value; when the distance between the measured points is greater than the maximum correlation distance, the value tends to be stable; the variance of the InSAR residual deformation is obtained through this step Sample values ​​with spacing h; S52. Use the classic spherical function model to fit the InSAR data samples at each time to obtain the nugget, range, sill and partial sill value model parameters; the nugget reflects the variance of the InSAR spatially uncorrelated noise ω(s,t) The base reflects the small-scale subtle deformation related to space. t Variance of (s) Through the above operations, the variance of the InSAR spatially correlated small-scale subtle deformation and the variance of the observation noise are separated; S53, the variance of the InSAR small-scale subtle deformation and the variance of the observation noise obtained in step S52 are used to determine the accuracy ratio k of the two observation methods by using variance component estimation, and then the GNSS observation noise is calculated according to the following formula: The GNSS observation noise is obtained as: The unknown parameters of S54 and cSTRE models are: small-scale subtle deformation ξ t (s), the state quantity x of additional constraints t , the state transfer matrix H and the variance in the state transfer process Before solving, the rest of the parameters are set to initial empirical values, and the initial value of the state transfer matrix H is set to: H=ρE Where 0<ρ<1 is the multiplication factor and E is the identity matrix; The parameters obtained in the above steps are brought into the fusion model for solution; the fusion model is solved by forward filtering and backward smoothing; that is, the following are calculated in sequence: state quantity prediction, state quantity prior variance estimation, calculation of spatiotemporal Kalman gain, posterior optimal estimation of state quantity and its variance, constraint state quantity and its variance estimation; then the above process is repeated for all times, and finally the entire fusion process is iterated using EM.

7. The method for monitoring mining area surface deformation by fusing InSAR and GNSS data with a spatiotemporal Kalman model according to claim 1 is characterized in that: The S6 comprises the following sub-steps: S61. Obtain model parameters and reconstruct surface deformation with high temporal and spatial resolution: Among them, s represents the position of the spatial point after reconstruction, which is consistent with the distribution of InSAR points, and t represents the time of the reconstruction point, which is consistent with the GNSS observation time; S62, verifying the accuracy of the reconstructed deformation monitoring results in time and space respectively; In terms of time, a part of the real surface GNSS data that is not involved in the fusion is used to calculate the root mean square error between it and the fusion result; In terms of space, a part of InSAR data is deleted before reconstruction so that it does not participate in fusion. Then the fusion result is used to calibrate with this part of InSAR data and calculate its root mean square error. The above two processes are used to verify the reconstruction accuracy of deformation in time and space respectively.

Images