Structural detection study of mine blanket area based on dynamic 3D deformation prediction

3.1. Design of inversion method for structural parameters of mining airspace area

The inversion method of structural parameters of mining airspace is a method to infer the structural parameters of the underground mining airspace in the mining area through the observed phenomena such as surface settlement and changes in the water table in the mining area. Interferometric synthetic aperture radar (InSAR) technique is a method that uses SAR interferometry to detect and monitor ground deformation. The Post-Inversion Method (PIM) model is a commonly used method for processing InSAR data in the inversion of geometric parameters (GP) of the subsurface voided area.The PIM model obtains the GPs of the subsurface voided area, such as the size, shape, and location of the voided area, by inverting the InSAR data. However, the PIM model also has some disadvantages. For example, PIM models need to pre-set some parameters, such as the smoothness, coherence and noise level of the surface and subsurface media, etc. In addition to the key voiding zone GPs, PIM also requires a set of modelling parameters. These model parameters usually need to be determined from actual measurements and will vary with the mining geological conditions. In the PIM model, the angle of the main influence angle tangent (IAT) depends on the dominant nature of the overlying rock formation, which makes the difference between the main IAT in the strike and tendency downhill and uphill directions smaller, and it can be approximated by assuming that the main IAT is equal in all three directions [19]. In addition, the inflection deviation in the model is only used to calculate the calculated boundary lengths of strike and tendency, which are positively correlated with the strike boundary length and tendency boundary length. In the PIM model, the tangent of the main influencing angles of strike and dip can be approximately assumed to be equal, while the inflection point offset is only used to calculate the boundary length of strike and dip, which is positively correlated with the length of strike and dip boundaries. To simplify the calculation, it can be assumed that these parameters are equal. In the process of acquiring the surface deformation field by InSAR technology, the geometric 3D projection relationship of image imaging is shown in Fig. 1.

Fig. 1Image imaging geometry 3D projection relationship

In three-dimensional projection, the vector projection sum in the direction of oblique moment is the line-of-sight direction deformation. Based on the three-dimensional geometric relationship, a linear equation is established as shown in Eq. (1) [20]:

where, the line-of-sight deformation serves as $D_{LOS}$ and the surface deformation components in the north-south, east-west and vertical directions are $D_N$, $D_E$ and $D_W$, respectively. The radar incidence angle serves as $\theta$ and the satellite heading azimuth is $\sigma$. The relationship between line-of-sight deformation and the GPs of the mine site is shown in Eq. (2) [21]:

where, the horizontal deformations of any point on the surface of the mining area from north to south and from east to west are $U(x,y,{\phi }_N,GP,MP)$ and $U(x,y,{\phi }_E,GP,MP)$, and the vertical deformation is $W(x,y,GP,MP)$, respectively. However, InSAR data are usually high-resolution and high-noise, so they need to be downsampled to improve the computational efficiency and reduce the noise effect. Therefore, the study chooses the adaptive quadtree downsampling method to process the InSAR results, which can obtain the surface deformation information more efficiently and provide more accurate data support for the construction of the inverse model of the GPs of the mining area. Quadtree is a data structure that can dynamically adjust its structure according to changes in data to achieve optimal storage and query performance. The quadtree downsampling algorithm is an image sampling method based on the quadtree data structure. Firstly, divide the original image into equally sized quadtree structures, and then sample the nodes in the quadtree according to the set sampling strategy. During the sampling process, different sampling methods can be used, including random sampling, minimum coding sampling, clustering sampling, etc. In the quadtree downsampling algorithm, the selection of sampling strategy has a significant impact on the sampling results, and needs to be selected according to specific application scenarios and requirements. In InSAR processing, the adaptive quadtree downsampling method can effectively reduce data volume, improve data processing efficiency, and ensure data accuracy and integrity. Considering the coupling characteristics between the parameters, the study will improve the GP inversion method of mining airspace area by cross iteration, and reduce the dependence of model parameters by dividing the parameters. The crossover genetic fitness function is shown in Eq. (3):

where, the genetic fitness function is $f$, and the line-of-sight direction deformation of the predicted point is $D_{LOS}^{\mathrm{\text{'}}}(x,y)$. The parameters to be solved by the genetic algorithm are shown in Eq. (4):

where, the solution of the parameter to be solved is $P_{req}$, the estimated value of the parameter to be solved is $P_i^1$, and the number of parameter valuation groups is $k$. The specific cross iteration is as follows: firstly, the first iteration is carried out, assuming that the initial PIM model parameters are known and their empirical values are determined, and then the genetic algorithm is used to invert another set of parameters. Then the second iteration is carried out, the parameters obtained from the inversion are regarded as known, and the remaining parameters are inverted. In the iteration process, iteration termination judgement is required, and the iteration is judged to be terminated when the set threshold is met, at which time the GPs of the extraction zone are known. If it does not meet the judgement conditions, the initial PIM model parameters obtained are used to repeat the iterative update until the threshold conditions are met. The whole process of updating parameters through continuous iteration, and finally get the GPs of the mining zone that meet the set threshold. The overall flow of the structural parameter inversion method for the mining air-mining zone is shown in Fig. 2.

Fig. 2Overall flow of the inversion method for structural parameters of the mining area's air-mining zone

3.2. Design of dynamic three-dimensional deformation prediction method for mining airspace area

In recent years, ground deformation and movement caused by mining has become a serious geological disaster. In order to prevent and control this disaster, the study proposes a method to transform the three-dimensional deformation prediction of the mining area into a one-dimensional line-of-sight orientated deformation prediction, which is based on the constraint relationship between horizontal movement and vertical deformation gradient, and uses spatio-temporal Kalman filtering for spatio-temporal data processing. At the same time, the adoption will also use the structural parameter inversion method for the mine's air-mining area to obtain a priori constraint parameters and solve the dynamic 3D deformation of the mine based on the predicted dynamic line-of-sight deformation. The time-varying linear combination of spatial field and temporal components in the dynamic 3D model $L_t$ is shown in Eq. (5) [22]:

where, the spatial field $H=[h_1,h_2,\cdots ,h_p]$ dimension of the observation data is $n\times p$; the dimension of the dynamic time component $X\left(t\right)=\left[x_1\right(t),x_2(t),\cdots ,x_p(t\left)\right]$ is $p\times 1$, and the system state of the mean component of the observation data serves as $X\left(t\right)$. The representation of the observational data is showcased in Eq. (6) [23]:

where, the observation data is $Y_t$ and the observation error component is ${\omega }_t$. After Kalman filter modelling, the state equation is calculated as showcased in Eq. (7):

where, the system state vector (SV) is $X\mathrm{\text{'}}\left(t\right)$; the state transfer matrix is $K$, and the state update at the time of $t$ is $\epsilon \left(t\right)$. The initial parameter estimation flow is shown in Fig. 3.

Firstly, the initial SV and state transfer matrix are set to meet the requirements of the spatio-temporal Kalman filtering model. Define the observation vector and establish the relationship with the SV. According to the nature of multivariate Gaussian distribution, set the mean vector and covariance matrix of the SV and observation vector. Use the maximum expectation algorithm for iterative parameter estimation to obtain the estimated value and covariance matrix of the SV. In this process, first initialize the SV, initialize the observation sequence, and initialize the model parameters, and then perform multiple iterations. In each iteration process, the probability of the state sequence is calculated based on the observation sequence and model parameters, and then the SV and covariance matrix are updated using the maximum expectation principle. This process will continue until convergence occurs. Finally, the estimated SV values and covariance matrix are obtained. The SV is filtered using the Kalman filter formula and the observation vector, and the filtered SV is interpolated to obtain the estimates of the missing data. Finally, the interpolated SVs and the state transfer matrix are used to predict the future SVs [24]. The kriging prediction for a point in the spatio-temporal field construction ${\widehat{L}}_t\left(s\right)$ is shown in Eq. (8):

where, given the spatial trend as $f\left(s\right)$ and the row vector of the spatial trend $F$ as ${\sigma }_{\xi }\left(s\right)$, the trend matrix $A$ and the bending energy matrix $B$ are calculated as shown in Eq. (9):

Fig. 3Initial parameter estimation flow

In Eq. (9), the spatial covariance matrix constructed by ${\sigma }_{\xi }$ is ${\sum }_{\xi }$. The spatial trend component is removed by least squares as shown in Eq. (10):

where, the spatial trend component is $Y_d\left(t\right)$, and the experimental semi-variance is calculated as shown in Eq. (11):

where, the experimental semivariance fitting function is ${\widehat{\gamma }}_d$, which is usually chosen to be expressed in the form of an exponential or spherical semivariance function. The spatial covariance is calculated as shown in Eq. (12):

where, the abutment value is $C$ and the distance operator is $‖\cdot ‖$. The spectral decomposition of the bending energy matrix $B$ is calculated as shown in Eq. (13):

where, the eigenvectors of the matrix $B$ are the column vectors $u_i$ of $U$, and its eigenvalues are the diagonal elements ${\lambda }_i$ of the diagonal matrix $\mathrm{\Lambda }$. After the time-varying linear combination of the spatial field and the time component is represented by the linear combination of the eigenvectors, the spatial field can be shown in Eq. (14):

where, the spatial trend $F=\left[f^T\right(s_1),f^T(s_2),\cdots ,f^T(s_n\left)\right]$. In order to solve for the observation valuation, it is first necessary to compute an estimate of the SV based on the known observations and the SV prediction. Then, based on the estimated value of the SV and the corresponding gain, the predicted value of the state covariance is calculated. Next, for the missing observations, the filtering results of the observations at moment $t$ can be derived by setting the corresponding matrix value to 0 and then based on the SV estimate. Finally, if it is necessary to predict the observation value at a certain moment, the state prediction at that moment is calculated, and then the prediction of the observation value is calculated using the state prediction and the corresponding gain. When the basic model is constructed, the cumulative line-of-sight to deformation sequence is shown in Eq. (15):

where, the cumulative line-of-sight deformation sequence with respect to the $t_1$ time is $d_{los}$, and the cumulative line-of-sight deformation sequence after spatio-temporal interpolation filtering is shown in Eq. (16):

where, the cumulative line-of-sight deformation sequence after spatio-temporal filtering difference is $d_{los}^1$. After spatio-temporal interpolation filtering, as well as state transfer matrix and spatial field processing, the dynamic prediction of the cumulative line-of-sight deformation of the mine surface is shown in Eq. (17):

The cumulative line-of-sight to deformation dynamics in Eq. (17) is predicted as $d_{los}^{predict}$. Conventional methods require InSAR monitoring data from at least three tracks, this study performs dynamic 3D deformation prediction of the mine site by basing on single-track InSAR monitoring data. The method reduces the number of unknowns to $m$ by constructing the constraint relationship between horizontal deformation and vertical deformation, which enables the decomposition of 3D deformation. The study establishes a basic model from InSAR time series processing to dynamic 3D deformation prediction through the functional relationship between horizontal deformation and vertical deformation gradient. In summary, the dynamic three-dimensional deformation prediction method for mining airspace is shown in Fig. 4.

4.1. Analysis of the application of structural parameter inversion methods for mining airspace areas

Aiming at studying the influence of the proposed method of inversion of structural parameters in the mining area, the experiment takes the San Yuan coal mine as the research object and analyses the structural parameters of the empty mining area. In GeoStudio, a geophysical exploration software, a three-dimensional geological model is constructed according to the actual mining area, reflecting the geological structure and lithological distribution of the mining area. The structural parameters of the goaf in the mining area were experimentally analyzed, and the accuracy of cross iterative estimation and non iterative estimation was compared. At the same time, the Kalman filtering model was used to predict the dynamic three-dimensional deformation of the goaf in the mining area, and compared with the Weibull model. The experimental factors mainly include the specific geological conditions of the mining area, the accuracy and performance of geophysical exploration equipment, etc. During the experiment, the measurement equipment mainly includes geophysical exploration equipment, including GeoStudio software, Global Positioning System (GPS), etc. The seismic instrument model is Trimble T8, the electromagnetic instrument model is GEM Systems GSM-19T, and the magnetic instrument model is Bruker GSM-29, The GPS receiver model is Trimble R10. GeoStudio software version is 2021.1, including Seismic Electrical, Magnetic, and GPS modules. The accuracy and reliability of GeoStudio software are directly related to the construction of geological models, while the accuracy and performance of GPS directly affect the accuracy of deformation prediction. The statistical analysis of measurement data is shown in Table 1.

In Table 1, the statistical analysis results of the measurement data indicate that the errors in length, width, and thickness are ±0.01 m, ±0.01 m, and ±0.07 m, respectively, with corresponding sensitivities of 83 %, 92 %, and 91 %, and confidence intervals of 95 %. Measurement data has high accuracy and reliability, and can be used as a basis for practical applications. The experiment evaluates the accuracy and precision of the inversion method for structural parameters of goaf in mining areas by comparing the three-dimensional surface maps of simulated goaf and actual estimated goaf, as well as parameter prediction results. Then, cross iterative estimation and non-iterative estimation methods were used to predict geometric parameters, and their accuracy was compared. The three-dimensional surface map and parameter prediction of the mine's empty mining area are shown in Fig. 5.

Fig. 53D surface map and parameter prediction of the mine's air extraction zone

a) 3D surface diagram

b) Parameter prediction

Fig. 5(a) shows the 3D surface map of the mining airspace area, and the aggregate shapes of the simulated airspace area and the actual estimated airspace area are similar, and the vertical projections basically coincide with each other. Fig. 5(b) shows the combined parameter prediction results of the mining airspace area, with MV denoting the measured data and PV denoting the predicted data. The absolute error of parameter prediction is in the interval of 0-63.8 m, and the absolute error is in the range of 0-8 %. The results show that the structural parameter inversion method for mining airspace zones proposed by the study has high accuracy and precision in both the shape of the airspace zones and parameter prediction. For verifying the reliability of the inversion of the GPs of the mining airspace. The study analyses the accuracy of cross iterative inversion based on the predicted and measured values of GPs. The comparison of predicted and measured values of GPs is shown in Fig. 6.

Fig. 6Comparison of predicted and measured values of geometrical parameters

a) Length

b) Width

c) Thickness

In Fig. 6, MV represents the measured value and IV represents the inversion value. The absolute error is the difference between the measured value and the true value, and the relative error is the ratio of the absolute error to the true value. Fig. 6(a) showcases the comparison results of the length estimation and inversion value of the open area, which shows that the absolute error of the length inversion value is less than 0.2 m, and the relative error is less than 1 %. Fig. 6(b) showcases the comparison results of the estimated width of the open area and the inversion value, the absolute error of the inversion value of the visible width is less than 6m, and the relative error is less than 5 %. Fig. 6(c) showcases the comparison results of the estimated thickness of the hollow zone and the inversion value, which shows that the absolute error of the thickness inversion value is less than 0.3 m, and the relative error is less than 8 %. The results show that the estimates of the goaf geometric parameters (length, width and thickness) obtained by the inversion method have high accuracy and reliability compared with the estimated value. The comparison of the accuracy of cross iterative estimation and non-iterative estimation is shown in Fig. 7.

Fig. 7Comparison of accuracy of cross iterative estimation and non-iterative estimation

a) Length indicators

b) Angle indicators

In Fig. 7, With denotes cross iterative estimation and Non denotes non-cross iterative estimation. Fig. 7(a) showcases the comparison of the inversion errors for the length class parameters, and the non-crossed iterative estimation increases the relative error of the length class parameter estimation by about 18 %. Fig. 7(b) showcases the comparison of the inversion errors of the angle class parameters, and the non-cross iteration estimation increases the relative error of the angle class parameter estimation by about 24 %. Comprehensively, it showcases that the cross iteration estimation has a better effect compared with the non-cross iteration estimation. The effect of the PIM initial model error on the inversion results of the air-sea zone is shown in Fig. 8.

In Fig. 8, S, T and O denote the subsidence coefficient, tangent angle and offset in the initial parameters of the PIM model, respectively. Fig. 8(a) showcases the results of the effect of the error content of the initial parameters of the PIM model on the inversion of coal seam thickness. Fig. 8(b) showcases the results of the influence of the error content of the initial parameters of the PIM model on the inversion of the length of the mining zone. Fig. 8(c) showcases the results of the influence of the error content of the initial parameters of the PIM model on the inversion of the dip angle of the air-mined zone. When the error between the sinking coefficient and the tangent of the main influencing angle increases, the differences between the boundary length, boundary width, mining direction azimuth and the true value can be basically ignored. When the initial model parameter error of any PIM increases, there is no significant change in the accuracy of coal seam thickness and average mining depth. It demonstrates that the error content of the initial parameters of the PIM model has no obvious influence on the inversion accuracy of the geometrical parameters of the extraction zone.

Fig. 8Effect of PIM initial model error on the inversion results of the extraction zone

a) Thickness

b) Length

c) Angle

4.2. Analysis of the application of dynamic 3D deformation prediction methods for mining airspace areas

Aiming at verifying the prediction effect of the Kalman filter model (KFM) in the dynamic three-dimensional deformation prediction method of the mining area, the KFM was utilized for forecasting the line-of-sight direction deformation after the spatio-temporal field matrix and the unknown parameters were established. The difference-like results predicted by the KFM are shown in Fig. 9.

Fig. 9(a) showcases the prediction effect of the KFM for date 1, with a peak prediction difference of less than 40 percent and an average prediction error of less than 7 per cent. Fig. 9(b) shows the prediction effect of the KFM for date 2, with a peak prediction difference of less than 40 % and an average prediction error of less than 8 %. Fig. 9(c) shows the prediction effect of the KFM for date 3, with the peak prediction difference less than 40 % and the average prediction error less than 10 %. In summary, it can be seen that the KFM prediction error will gradually increase with the increase between, which is a common problem of model prediction, but the overall error is less than 10 %. Overall, the KFM shows a strong prediction ability. The results of the deformation prediction of the Global Positioning System (GPS) position in different directions are shown in Fig. 10.

Fig. 10(a) shows the deformation prediction of the GIS location in vertical direction with the root mean square error of 1.5 cm, Fig. 10(b) shows the deformation prediction of the GIS location in east-west direction with the root mean square error of 1.7 cm, and Fig. 10(c) shows the deformation prediction of the GIS location in north-south direction with the root mean square error of 0.8 cm.

Fig. 9Disparity-like results predicted by the KFM

a) Data 1

b) Data 2

c) Data 3

Fig. 10Deformation prediction results for different orientations of GPS positions

a) Vertical direction

b) East-west direction

c) North-south direction

The results show that the proposed dynamic 3D deformation prediction method for mining airspace area shows high accuracy in predicting the deformation of mining area. deformation prediction method shows high accuracy in predicting the deformation of the mining area, which helps to better understand and predict the deformation of the mining airspace. In order to verify the application of the proposed dynamic 3D deformation prediction method based on the KFM, the Weibull model is used as a comparison in the experiment, and the fitting of the Weibull model is shown in Fig. 11. The Weibull model is a continuous probability distribution model used to describe the probability of a random variable reaching or exceeding a certain value at a certain moment. This model is based on a Weibull function, which is a combination of a power function and a logarithmic function.

Fig. 11The fit of the Weibull model

a) Point 1

b) Point 2

Fig. 11 showcases the fit of the Weibull model before the time-zero correction. Fig. 11(a) shows the fit of the Weibull model to measurement point 1. Fig. 11(b) showcases the fit of the Weibull model to measurement point 2. It demonstrates that the fitting accuracy of the Weibull model is poor and easily limited by the time zeros. Although the Weibull model can characterise the distribution of random variables, it may not be able to adequately capture the dynamics of the data when dealing with time series data. A comparison of the model's line-of-sight deformation prediction results is shown in Fig. 12.

Fig. 12Comparison of line-of-sight deformation prediction results of models

a) GPS

b) InSAR

In Fig. 12, KF denotes the KFM. Fig. 12(a) shows the fitting of GPS data for the KFM and the Weibull model, and the KFM has a better fit. Compared to the Weibull model, the prediction error of the KFM is reduced by 8 %. Fig. 12(b) shows the comparison between the fitting of the KFM and the InSAR measurement data, which shows that both of them have the same trend and the spatial characterisation is basically the same. It can be seen that the dynamic 3D deformation prediction method based on the KFM has better application and practicability when dealing with problems such as dynamic 3D deformation prediction in mining areas. During the period of 20210915, the deformation prediction comparison of the goaf in the mining area is shown in Fig. 13.

Fig. 13Prediction and comparison of goaf deformation in mining area

a) InSAR-horizontal

b) KFM-horizontal

c) Weibull-horizontal

d) Logisitc-horizontal

e) InSAR-vertical

f) KFM-vertical

g) Weibull-vertical

h) Logisitc-vertical

Fig. 13(a-d) represent the horizontal gradients predicted by each model, while Fig. 13(e-h) represent the vertical gradients predicted by each model. The black dashed line surrounding the area indicates a difference between the predicted value and the actual value. It can be seen that there is no significant difference between the predicted deformation of the KFM proposed in the study and the InSAR measurement data, indicating that the spatial continuity between the two is basically the same. However, there is a difference between the predicted deformation of Weilbull and Logisitc models and the InSAR measurement data. The results indicate that the proposed KFM has the highest prediction accuracy and can accurately predict the dynamic three-dimensional deformation of mining areas. In order to verify the effectiveness of the dynamic 3D deformation prediction method based on the KFM proposed in the study, the experiment will compare the mine subsidence prediction method based on InSAR technology and LSTM algorithm proposed in 2023 [25]. The experiment will use maximum absolute error (MAE), maximum relative error (MRE), mean absolute deviation (MAD), and Wilmot consistency index as evaluation indicators. The comparison results of the two methods are shown in Table 2.

From Table 2, it can be seen that the two methods have similar performance in prediction accuracy and both have high prediction accuracy. The methods proposed in the study on MAE, MRE, and MAD are slightly lower than those proposed in 2023. The Wilmot consistency index of the two methods is equal. The results indicate that the proposed method performs slightly better in prediction and evaluation than the method proposed in 2023, but the two methods are comparable in reliability.