Space-time kalman model fusion insar and gnss data reconstruction mine area surface deformation monitoring method

Abstract

The application discloses a kind of space-time Kalman model fusion InSAR and GNSS data reconstruction mining area surface deformation monitoring method.Aiming at the uneven spatial variation of deformation, large gradient, time variation of highly nonlinear mining area, the application makes full use of the advantages of InSAR deformation high spatial resolution and GNSS deformation high time resolution.According to the spatial characteristics of mining area deformation, the spatial model is adaptively constructed, and the automation and efficiency of the algorithm are improved.Meanwhile, the time characteristic constraint of prior deformation of mining area is added, which not only reduces the dependence of traditional STRE fusion method on GNSS data, but also improves the accuracy of deformation fusion.The method of the application can consider InSAR deformation monitoring result, GNSS deformation monitoring result and prior deformation characteristics of mining area at the same time, effectively fuse the three kinds of information to the greatest extent, and realize high space-time resolution deformation monitoring of mining area.Under the background of strong demand for mining resources in China, the application has great significance for accurate monitoring, evaluation and prevention of mining 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 geodetic surveying based on remote sensing images, and in particular to a method for monitoring surface deformation in mining areas by fusing InSAR and GNSS data using a spatiotemporal Kalman model. Background Art

[0002] Intensive coal mining has a significant impact on the surface, geological structure, and surrounding environment of mining areas, even triggering a series of geological disasters such as landslides, debris flows, and collapses, seriously affecting the lives and property of local residents. Deformation is the most intuitive precursor to geological disasters. Therefore, using scientific methods to monitor surface deformation in goaf areas is extremely important for accurately assessing and preventing geological disasters in mining areas. Over the past 20 years, Synthetic Aperture Radar (InSAR) technology, with its unique advantages of high spatial resolution, high precision, and rapid monitoring over large areas, has rapidly achieved significant breakthroughs and achievements in various fields such as mining areas, landslides, earthquakes, the atmosphere, permafrost, and urban infrastructure, becoming 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. This results in a significant lag in the deformation monitoring. This makes it difficult to timely reflect the deformation status of sudden and severe deformation gradients in mining areas and landslides, significantly limiting the development and application of InSAR technology. The Global Navigation Satellite System (GNSS) boasts advantages such as real-time, high-precision, and long-term continuous observation, making it widely used for safety monitoring of various geological hazards and engineering applications. However, GNSS is relatively expensive and can only monitor a few key locations, preventing large-scale, high-density deployment. Consequently, its monitoring results have low spatial resolution, making it difficult to reflect spatial variations in deformation within the monitored area, especially in mining areas with large deformation gradients.

[0003] The spatiotemporal Kalman filter (STF) is a Kalman filtering method that fully accounts for spatiotemporal correlations. It has demonstrated significant advantages in the field of spatiotemporal deformation analysis, enabling dynamic spatiotemporal filtering, spatiotemporal interpolation, and spatiotemporal prediction of deformation monitoring data. In particular, the spatiotemporal random effects (STRE) model within STF can also reduce data dimensionality, significantly improving the computational efficiency of spatiotemporal big data. Therefore, some researchers have proposed combining the advantages of the high spatial resolution of InSAR deformation with the high temporal resolution of GNSS deformation, using a STF method to fuse InSAR and GNSS monitoring data to reconstruct surface deformation with high spatiotemporal resolution. This method has been applied to fault deformation reconstruction in Southern California and other areas, achieving promising results. However, these studies still face two challenges: First, the STRE model employs spatial bases during the fusion process to reduce dimensionality in order to achieve a fast solution, and the fewer spatial bases used, the better the dimensionality reduction. Currently, spatial bases are primarily deployed uniformly at multiple scales and rely on empirical guidance. While this approach can meet the needs of uniform deformations such as faults, it is not suitable for mining areas with non-uniform deformations and large gradients. Secondly, when InSAR and GNSS are fused to generate surface deformation results with high spatiotemporal resolution, the deformation accuracy of the STRE model fusion is closely related to the number of GNSSs. The more GNSSs there are, the higher the fusion deformation accuracy, while the fewer GNSSs there are, the more the fusion results are dominated by the InSAR deformation results. However, mining areas are generally unable to deploy a large number of GNSS stations, making existing spatiotemporal Kalman fusion methods unsuitable for mining areas with complex spatiotemporal deformation variations. Therefore, in response to the current demand for mining area deformation monitoring, it is urgent to consider how to optimize the spatial modeling method of mining area deformation during the spatiotemporal Kalman fusion process and address the problem of low fusion deformation accuracy caused by the small number of GNSS stations in mining areas. Summary of the Invention

[0004] The present invention aims to overcome the drawbacks of the current spatiotemporal Kalman method, which fuses InSAR and GNSS deformation data to reconstruct high spatiotemporal resolution of the surface deformation, making it difficult to apply to mining areas with large deformation gradients and a small number of GNSS data. This method proposes a spatiotemporal Kalman model to fuse InSAR and GNSS data to reconstruct surface deformation monitoring in mining areas. This method adaptively constructs a spatial model based on the spatial deformation characteristics of the mining area, improving the algorithm's automation and efficiency. Furthermore, the addition of a priori temporal characteristic constraints on the mining area's deformation not only reduces the traditional STRE fusion method's reliance on GNSS data but also further improves the accuracy of deformation fusion.

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

[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 spatiotemporal filtering before being geocoded. The average deformation rate and time series deformation results of the study area are then calculated using 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 spatiotemporal benchmarking.

[0008] S2. Establish a spatiotemporal Kalman model that fuses InSAR and GNSS deformations. First, establish the STRE model. Then, based on the acquired InSAR and GNSS data, establish the fused spatial and temporal models 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: Using the InSAR average deformation rate as a benchmark, first lay out a uniform spatial basis across the entire map to capture global smooth deformations. Then, calculate the residuals after modeling the first spatial basis, identify the parts with larger median values, and perform a morphological dilation operation on them. A second spatial basis is then laid out for the resulting irregular areas to capture local deformations. Next, calculate the residuals after modeling the total of the first and second spatial basis models, identify the parts with larger median values, and perform a morphological dilation operation on them. A third spatial basis is then laid out for the resulting irregular areas to capture small local deformations. Finally, merge the three InSAR spatial basis layers and the GNSS spatial basis to complete the adaptive spatial modeling of large-scale major deformations.

[0010] S4. Based on the temporal deformation characteristics of the mining area, we introduce prior information about the mining area and apply logistic model constraints to the state variables. This improves the estimation accuracy of the state variables during periods without InSAR observations while reducing reliance on GNSS data. Finally, we construct a spatiotemporal Kalman model that fuses InSAR and GNSS deformation, constrained by the spatiotemporal characteristics of the deformation.

[0011] S5. Determine the model parameters during 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 during 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 during the entire fusion process.

[0012] S6. Reconstruct the deformation results of the surface with high temporal and spatial resolution: Reconstruct 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 includes:

[0014] S11. The acquired radar images are geocoded to form an interferogram after registration, small baseline meshing, interferometry, flatland and terrain removal, phase filtering, phase unwrapping, large-scale trend removal, atmospheric correction, and spatiotemporal filtering. Note that since no deformation trends are subsequently modeled, this step requires removing the trend term from the interferogram.

[0015] S12. Assume that there are N SAR images in the study area, and M interferograms are generated according to the spatiotemporal baseline threshold. There is only one small baseline set. Then the following equation can be constructed using the Small Baseline Subsets (SBAS) method:

[0016]

[0017] The least square method of the above formula can be obtained:

[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 obtaining the deformation sequence, 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 only obtain deformation in the radar line of sight (LOS), it is necessary to convert the GNSS three-dimensional direction to the radar line of sight:

[0022]

[0023] in Indicates 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 deformation of the GNSS point respectively, θ represents the local incidence 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, the InSAR point closest to the GNSS point is found as the GNSS synonym point.

[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 includes:

[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 point 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 small-scale, subtle deformations in space. Note that if the deformation contains a trend term, this term can also be spatially modeled. Since trends are generally removed in advance during 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 x t It can be described by a linear combination of:

[0035]

[0036] Among them S t (s)=[S 1，t (s), S 2，t (s),..., S r,t (s)] T, represents the established r x n dimensional space basis function, x t = [x 1，t (·), x 2，t (·),..., x r,t (·)] T represents the r x m zero mean Gaussian random state quantity. Note that the basis function can change over time, but for a mine area it can be very easy to capture the deformation area by averaging the deformation rate, and also to reduce the complexity of the model, so it is modeled as a function that does not change over time, i.e. S(s) = S t (s) = [S1(s), S2(s),..., S r (s)] T .

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

[0038]

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

[0040] S22, the state quantity x t will evolve over time to form the state transition equation:

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

[0042] where H is an r x r state transition matrix, is the noise in the state transition process. From the above, the final form of the STRE model can be established:

[0043]

[0044] From the above analysis, the model parameters of STRE are the small-scale fine deformation ξ t (s), the state quantity x t , the state transition matrix H and the variance

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

[0046] Then the spatio-temporal random effect model of the fusion of InSAR and GNSS established according to the above S21 STRE is:

[0047]

[0048] The meanings of the respective variables are as follows in Table 1:

[0049] Table 1 meanings of variables

[0050]

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

[0052] Finally, in step S22, the state quantity x t Evolution in time constitutes a state transition equation, forming a final fused spatio-temporal Kalman model.

[0053] According to some preferred embodiments of the present application, the S3 comprises:

[0054] S31, in order to fully capture the spatial variation of deformation in step S23, a plurality of layers of spatial bases are arranged at multiple scales, which is equivalent to the nesting of a plurality of variogram functions. However, deformation is not uniformly varied in space, so some adaptive adjustment needs to be made to the position arrangement of the spatial base.

[0055] First, a layer of uniform spatial base is arranged according to the average deformation rate of InSAR in the entire deformation field, and then the residual deformation remaining after modeling by the layer of uniform spatial base is calculated:

[0056]

[0057] Wherein, V1 represents the residual error after modeling by the first layer of spatial base, Z a represents the average deformation rate of InSAR, and S1 represents the r1 spatial bases arranged in the first layer.

[0058] S32: Identify the areas with larger values ​​in V1. This process will result in some isolated islands and smaller areas. Therefore, it is necessary to first remove these smaller areas and then perform a morphological dilation operation to fill the holes and expand the boundaries to obtain areas where the first spatial basis fails to fully capture the deformation.

[0059] S33. For the obtained irregular area, a second layer of uniform spatial basis (assuming r2) is laid out with a higher density, and combined with the first spatial basis to calculate the 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 spatial basis still cannot fully capture all deformations. Therefore, it is necessary to lay out a third-layer spatial basis with a smaller spacing on this basis. Similar to step S32, the part with a larger value in V2 is identified. In this process, some isolated islands and small areas will appear. Therefore, it is necessary to first delete the small areas and perform a morphological dilation operation at the same time to fill the holes and expand the boundaries to obtain the areas that the first and second-layer spatial 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 , a total of r=r1+r2+r3 space bases.

[0066] In the above process, simply specifying the density of each spatial basis layer will yield the position and number of spatial bases for adaptation across the entire region. This not only reduces human intervention and improves the algorithm's automation, but also reduces the number of spatial bases, improving algorithm efficiency. This can significantly save time and space when calculating large-scale deformations.

[0067] It's worth noting that, in general, after three placements, the residuals are very small, sufficient for subsequent calculations. However, when the deformation gradient is very large, a three-scale spatial basis may still not be able to fully capture all deformations. In this case, additional spatial basis scales can be added to the three-scale basis to further improve the accuracy of spatial modeling. However, too many spatial basis scales consume a lot of computation time and may cause overfitting, so a compromise is necessary.

[0068] S36. Use all GNSS points to calculate the GNSS spatial basis S gnss At this point, the spatial bases of InSAR and GNSS are obtained respectively.

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

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

[0071]

[0072] Where F(t ins ) represents 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 expected 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 point, the relationship between the model prediction value and the state quantity can be established by the observation equation in step S22:

[0077]

[0078] Generally speaking, 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 variables 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] Among them F s (t,θ) represents the established spatially independent and time-dependent function constraint model, and θ represents the model parameters.

[0085] S43. The original STRE model was constructed based on a linear system. However, mining area deformation is highly nonlinear, and its state variables also exhibit highly nonlinear characteristics. Therefore, the original STRE model contains errors that are beyond the system. In other words, the complete system description differs from the description established by the standard space-time Kalman filter. In this case, 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 Denotes the optimal state quantity of the unconstrained Kalman estimate, and is expressed as Unconstrained Kalman estimation posterior variance, according to existing research, one solution is to project the unconstrained state estimate onto the constraint surface, so the following objective function is established according to 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, we can get:

[0094]

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

[0096] The above constraints indicate that when the variance of the Kalman-evaluated state variables is small, they are considered more reliable, and the model constraints are assigned a smaller weight. When the variance of the Kalman-evaluated state variables is large, they are considered to require correction, and the model constraints are assigned a larger weight. This ensures that the Kalman filter results do not deviate from the model fitting results. During design, care should be taken to keep state variables close to the GNSS points as constant as possible to avoid compromising GNSS accuracy.

[0097] At this point, we have obtained the 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 proceed with 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 deformations and noise, and the variance of these two types of information needs to be separated. The residual information after InSAR deformation spatial modeling of each scene is fitted with a semivariogram to obtain the variance values ​​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 the InSAR spatially correlated small-scale subtle deformation and the variance of the 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. Then, the GNSS observation noise can be 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 , state transfer matrix H and variance in the state transfer process

[0109] Before solving, the rest of the 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 fed into the fusion model for solution. This fusion model is generally solved using forward filtering and backward smoothing. Specifically, the following steps are performed: state prediction, prior variance estimation of the state, calculation of the spatiotemporal Kalman gain, a posteriori optimal estimate of the state and its variance, and estimation of the constrained state and its variance. This process is then repeated for all time steps. Finally, the entire fusion process is iterated using EM iteration 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 surface deformation with high temporal and spatial resolution:

[0115]

[0116] Where s represents the position of the reconstructed spatial point, 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 terms of space, in order to facilitate accuracy verification, 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 with this part of the InSAR data to calculate its root mean square error.

[0120] The above two processes can respectively verify the reconstruction accuracy of the method proposed in the present invention in terms of temporal and spatial deformation.

[0121] According to the above-mentioned method for reconstructing mining area surface deformation monitoring by fusing InSAR and GNSS data with spatiotemporal Kalman model, a deformation monitoring system for reconstructing mining area surface with high spatiotemporal resolution by fusing InSAR and GNSS data 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 targets mining areas with uneven spatial deformation variations, large gradients, and highly nonlinear temporal variations. The proposed method for reconstructing mining area surface deformation monitoring by fusing InSAR and GNSS data with a spatiotemporal Kalman model fully utilizes the advantages of the high spatial resolution of InSAR technology and the high temporal resolution of GNSS technology. The spatial model is adaptively constructed according to the spatial characteristics of the mining area deformation, thereby improving 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 simultaneously take into account the InSAR deformation monitoring results, GNSS deformation monitoring results, and prior deformation characteristics of the mining area, effectively fusing these three types of information to the greatest extent, and realizing 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 results of the InSAR deformation sequence 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 The InSAR deformation rate and the 419 spatial basis of the traditional three-layer uniform layout in the embodiment;

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

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

[0131] Figure 8 yes Figure 7 residual distribution;

[0132] Figure 9 is the InSAR time series deformation results reconstructed by the traditional method and the method of the present invention and their residuals compared with the observed values;

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

[0134] Figure 11 The comparison between the reconstruction results of the deleted three-scene InSAR deformation and the SBAS-InSAR deformation results of the embodiment is shown. DETAILED DESCRIPTION

[0135] The present invention is described in detail below with reference to the embodiments and accompanying drawings. However, it should be understood that the embodiments and accompanying drawings are merely exemplary descriptions of the present invention and do not constitute any limitation on the scope of protection of the present invention. All reasonable variations and combinations within the scope of the inventive concept of the present invention fall within the scope of protection of the present invention.

[0136] This example uses 21 scenes of Sentinel-1SAR data from July 20, 2020 to March 17, 2021 in a certain area. Figure 1 The figure shows a specific spatiotemporal Kalman model fusing InSAR and GNSS data to reconstruct mining area surface deformation monitoring method.

[0137] Figure 2 The deformation rate (in cm / year) obtained using the SBAS-InSAR method is shown in the figure. The red star indicates the common reference point for InSAR and GNSS, which serves as the basis for the temporal and spatial unification of the InSAR and GNSS data. The 11 triangles represent the remaining GNSS stations. The black boxes indicate the selected area of ​​interest.

[0138] Figure 3 The following are 24 SBAS-InSAR time-series deformation monitoring results for the area of ​​interest. GNSS images DT01 and DT07, and InSAR images dated 20201024, 20201211, and 20210128, were subsequently excluded from the fusion process to verify the accuracy of the high-resolution spatial and temporal reconstruction of surface deformation.

[0139] Figure 4 yes Figure 2 The accuracy comparison results of 11 GNSS and InSAR show 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 This 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 an RMS of 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 an RMS of 1.17 cm. Therefore, the new method not only has fewer spatial bases but also has higher accuracy, making it more suitable for mining area deformation monitoring.

[0141] Figure 9 Figure 2 shows the InSAR time series deformation reconstructed using the traditional and new methods, along with their residuals compared to the observed values. The traditional method exhibits larger residuals in deformed areas, while the new method, due to the introduction of prior deformation features as constraints, exhibits smaller errors. Overall, the root mean square (RMS) residual of the new method is 10.97% lower than that of the traditional method, demonstrating its superiority.

[0142] Figure 10Comparisons of deformation measurements and reconstructions using high temporal and spatial resolution at two GNSS stations, DT01 and DT07, are presented. At DT07, due to the inherently high accuracy of InSAR, the STRE and cSTRE reconstructions are comparable. However, at DT01, although the new method achieves slightly higher accuracy than traditional methods, the accuracy is not very high due to InSAR observation errors.

[0143] Figure 11 The results of the reconstruction of the deleted three-view InSAR deformation are compared with those of the SBAS-InSAR deformation. The traditional method reconstructs larger deformation errors and underestimates deformation in local areas. The RMSE of the residuals between the traditional method and the SBAS-InSAR deformation are 0.47cm, 0.61cm, and 0.63cm, respectively. However, the new method, due to the introduction of the logistic model, has a very small error. The RMSE of the reconstructed deformation residuals are 0.42cm, 0.51cm, and 0.55cm, respectively. Compared with the traditional method, the new method improves accuracy by 10.64%, 16.39%, and 12.70%, respectively, with an average improvement of 13.24%.

[0144] The above analysis fully demonstrates the effectiveness of the method proposed in this paper, 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 merely preferred embodiments of the present invention, and the scope of protection of the present invention is not limited thereto. All technical solutions within the scope of protection of the present invention are within the scope of protection of the present invention. It should be noted that improvements and modifications that can be made by persons of ordinary skill in the art without departing from the principles of the present invention are also considered within the scope of protection of the present invention.

Claims

1. A method for monitoring surface deformation in mining areas using a spatiotemporal Kalman model that fuses InSAR and GNSS data, characterized by: The following steps are involved: S1. InSAR and GNSS data processing: The acquired InSAR interferograms are filtered, unwrapped, trended, atmospherically corrected, and spatially filtered before being geocoded. The average deformation rate and time series deformation results of the study area are then calculated using 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 fuses InSAR and GNSS deformations. First, establish a STRE model whose model parameters are small-scale subtle deformations, state quantities, state transition matrices, and variances during state transitions. Then, establish a fused spatial 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 InSAR average deformation rate as the benchmark, first lay out a uniform spatial basis on the entire map to capture global smooth deformations; then calculate the residuals after modeling the first-layer spatial basis, identify the parts with larger median values ​​of the residuals, and perform a morphological dilation operation on them. A second-layer spatial basis is laid out on the obtained irregular areas to capture local deformations; then calculate the residuals after modeling the total of the first and second-layer spatial basis, identify the parts with larger median values ​​of the residuals, and perform a morphological dilation operation on them. A third-layer spatial basis is laid out on the obtained irregular areas 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; S4. Based on the temporal deformation characteristics of the mining area, the prior information of the mining area is introduced to add Logistic model constraints to the state quantity, while reducing the dependence on GNSS data. Finally, a spatiotemporal Kalman model of InSAR and GNSS deformation is constructed with the spatiotemporal characteristics of deformation constraints. S5. Determine the model parameters during 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 during the fusion process through least squares fitting, forward filtering, and backward smoothing. Finally, use EM iterative estimation for the entire fusion process. S6. Reconstruct the deformation results of the surface with high temporal and spatial resolution: Reconstruct 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 mining area surface deformation by fusing InSAR and GNSS data with a spatiotemporal Kalman model according to claim 1 is characterized in that: The S1 includes the following sub-steps: S11. The acquired radar images are geocoded to form an interferogram after registration, small baseline meshing, interferometry, flatland and terrain removal, phase filtering, phase unwrapping, large-scale trend removal, atmospheric correction, and spatiotemporal filtering. This step removes trend terms from the interferogram. S12. Assume that there are N SAR images in the study area, and generate M interferograms based on the temporal and spatial baseline thresholds. If there is only one small baseline set, the following equation is constructed using the short baseline set SBAS method: The least square method of the above formula is: 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 by the coherence of the interference pattern; L is the phase of the interference pattern; T is the transpose of the matrix; After obtaining the deformation sequence, the least squares fitting is used to obtain the average deformation rate of the entire time period; S13. Convert the GNSS three-dimensional direction to the radar line of sight: in Indicates 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 radar incidence angle of the satellite, 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 synonym point; 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 mining area surface deformation by fusing InSAR and GNSS data with a spatiotemporal Kalman model according to claim 2 is characterized in that: The S2 includes 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 point 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 deformation in space; If the deformation contains a trend term, perform spatial modeling on the trend term; The main large-scale deformation in space v t (s) is modeled as a spatial field S t (s) and time-varying state quantity x t It can be described by a linear combination of: Among them 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 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 j-th point in the l-th spatial basis on the space s, j = 1, 2, ..., r; p l,j Indicates the spatial basis position of the jth point in the lth layer of 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 , variance during state transition S23. In the study area D, a total of n1 InSAR points are obtained, which have m1 observation dates T1; at the same time, a total of n2 GNSS points are obtained, which have m2 observation dates 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 InSAR small-scale subtle deformation and observation noise; 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 state transfer equation is formed by evolving in time, and the final fused space-time Kalman model is formed.

4. The method for monitoring mining area surface deformation by fusing InSAR and GNSS data with a spatiotemporal Kalman model according to claim 3 is characterized in that: The S3 includes the following sub-steps: S31. Adaptively adjust 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 laid out, Z a represents the InSAR average deformation rate, S1 represents the r1 spatial basis laid out in the first layer; S32: Identify the parts with large median values ​​in V1, which will result in isolated islands and small areas. First, delete the small areas and perform a morphological dilation operation to fill the holes and expand the boundaries, obtaining areas where the first spatial basis cannot fully capture the deformation. S33. For the obtained irregular area, a second layer of uniform spatial basis is laid out with a higher density, the number of which is r2, and the new residual is calculated by combining it with the first spatial basis: Use S 1,2 represents the r1+r2 space basis of the first and second layouts, then S 1,2 =[S1,S2] T ; S34. Lay out a third layer of spatial bases with small spacing; first delete the small areas, and then perform a morphological dilation operation to fill the holes and expand the boundaries, thereby obtaining areas that were not fully captured by the first and second layers of spatial bases; 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 4 is characterized in that: The S4 includes 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 ) represents t ins The InSAR shape variable at time , a, b, c are the model parameters to be estimated; After obtaining the parameters, the expected deformation of the model is obtained at any time t: in are the parameters estimated by nonlinear least squares; 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 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: Among them F s (t, θ) represents the established spatially independent and temporally dependent function constraint model, and θ represents the model parameters; S43. The original STRE model was constructed based on a linear system. However, the deformation in the mining area is highly nonlinear, and its state variables are also highly nonlinear. The original STRE model contains errors beyond the system, that is, the complete system description is different from the description established by the standard space-time Kalman filter. In this case, the additional state constraints improve the performance of the original STRE. use Represents the constraint state quantity at time t, and uses Denotes the optimal state quantity of the unconstrained Kalman estimate, and is expressed as Unconstrained Kalman estimation of the posterior variance, according to the minimum variance criterion, establishes the following objective function: 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, we get: Where J = WG'(GWG') -1 ; So far, we have obtained the estimate of the state quantity and its variance of the additional constraint. 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 made; 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 mining area surface deformation by fusing InSAR and GNSS data with a spatiotemporal Kalman model according to claim 5, 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. The variance of these two types of information is separated. The residual information 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, that is, 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 precision ratio k of the two observation methods, and then the GNSS observation noise is calculated according to the following formula: The GNSS observation noise is obtained as: The unknown parameters of the S54 and cSTRE models are: small-scale subtle deformation ξ t (s), the state quantity x of additional constraints t , state transfer matrix H and 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, and constraint state quantity and its variance estimation; then the above process is repeated for all time, 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, characterized in that: The S6 comprises the following sub-steps: S61. Obtain model parameters and reconstruct surface deformation with high temporal and spatial resolution: Where s represents the position of the reconstructed spatial point, 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. Verify the accuracy of the reconstructed deformation monitoring results in time and space respectively; In terms of time, a portion of the ground truth 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 the 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