Three-Dimensional Deformation Prediction Based on the Improved Segmented Knothe–Dynamic Probabilistic Integral–Interferometric Synthetic Aperture Radar Model

Abstract

Coal is the main mineral resource, but over-exploitation will cause a series of geological disasters. Interferometric synthetic aperture radar (InSAR) technology provides a superior monitoring method to compensate for the inadequacy of traditional measurements for mine surface deformation monitoring. In this study, the whole process of mining a working face in Huaibei Mining District, Anhui Province, is taken as the object of study. The ALOS PALSAR satellite radar image data and ground measurements were acquired, and the ISK-DPIM-InSAR deformation monitoring model with the dynamic probabilistic integral model (DPIM) was proposed by combining the probabilistic integral method (PIM) and the improved segmented Knothe time function (ISK). The ISK-DPIM-InSAR model constructs the inversion equations of InSAR line-of-sight deformation, north–south and east–west horizontal movement deformation, vertical deformation, inverts the optimal values of the predicted parameters of the workforce through the particle swarm algorithm, and substitutes it into the ISK-DPIM-InSAR model for predicting the three-dimensional dynamic deformation of a mining face. Simulated workface experiments determined the feasibility of the model, and by comparing the level observation results of the working face, it is confirmed that the ISK-DPIM-InSAR model can accurately monitor the three-dimensional deformation of the surface in the mining area.

1. Introduction

After the mineral resources underneath the ground have been extracted, the initial stress equalization of the bunker around the working face of the mine is disturbed [1,2,3]. At the same time, the surface of the working face will reach a new equilibrium that is different from the previous one [4], which causes warping, sustained motion, and discontinuity damage. The rock layers and the ground surface in the process [5] give rise to many geological problems in the mining area, such as surface subsidence and the destruction of buildings. Traditional mine surface monitoring establishes surface rock movement stations of observation along the underground working face’s incline and strike. It uses leveling and GPS survey techniques to measure the elevation and coordinates of the surface rock movement observation stations at certain intervals to obtain the outcomes of the horizontal movement and subsidence of the mine surface [6,7]. However, most of the mine surface rock movement observation stations show line distribution, which makes it challenging to depict the entire surface movement accurately [8,9,10]. Therefore, relying on the line data obtained from the surface rock movement observation stations, the accuracy of the inversion of the predicted parameters will be reduced, and the predicted results will be somewhat affected [11,12].

Because of its extensive spatial coverage, excellent precision, and low cost, interferometric synthetic aperture radar has emerged as a useful tool for monitoring surface deformation data [13,14]. Among them, the differential synthetic aperture measurement technique technology (D-InSAR) has been widely used to monitor surface subsidence in mining areas [15,16,17,18,19], but due to its limited monitoring gradient, it is not able to monitor the large-scale subsidence of the working face during the active period of mining [20,21], which severely restricts the implementation of InSAR technology for ground subsidence monitoring in mining regions [22,23]. Li and Yang et al. constructed the relationship equations between mining-induced vertical and horizontal deformation based on single InSAR interferometric pairs and based on geometric InSAR datasets, using the one-dimensional deformation to solve the time-series three-dimensional deformation [24,25], but this method is too reliant on high-quality and high-resolution synthetic aperture radar (SAR) imagery, and the spatial and temporal baselines, which will affect the accuracy of the final monitoring findings, and time-series InSAR techniques are similarly limited in their ability to monitor large-scale settlement in the working face [26,27]. Yang et al. constructed line-of-sight (LOS) deformation formulae according to the offset tracking technique of single interferometric pairs to acquire the north–south, east–west, and vertical deformation of the mining workforce [28,29], but the offset tracking technique requires high-resolution radar imagery and has a limited range of application [1]. In addition, there are limitations to the methods described above, as surface deformation caused by working face mining can be monitored [30] but the extent of the impact of underground mining of mineral resources on surface deformation at subsequent workings cannot be predicted, and it is impossible to realize the assessment of the potential risk of the mine [31,32].

Probabilistic integration methods are commonly applied to mine subsidence due to the ability to forecast the ultimate horizontal movement and surface settling brought on by subsurface mining [2]. The PIM-InSAR technique combines the InSAR technique and the probability integration method. Firstly, the monitoring of established LOS deformation in the working face area using differential synthetic aperture radar interferometry (D-InSAR) is conducted and then the deformation equations are constructed by combining the vertical, east–west, and north–south deformations predicted by the probability integration method, and then the predicted parameters of the probability integration method are inverted to obtain the three-dimensional deformations induced by the underground mining. Yang et al. constructed deformation equations utilizing PIM-InSAR technology and inverted the expected parameters of the probability integral method through a genetic algorithm [33,34], but PIM predicts the final deformation state of the ground surface and cannot predict the deformation results during the mining process. Yang and Jiang et al. combined the probabilistic integration method and the improved probabilistic integration method with the Knothe time function to acknowledge the dynamic prediction of three-dimensional deformation due to subterranean mining [35,36,37]. However, the Knothe time function itself has theoretical deficiencies, and its surface mobility which is derived from the real process of surface movement and deformation is inconsistent with the subsidence velocity and acceleration [38,39].

In this paper, a DPIM-InSAR method combining the probabilistic integration method and improved segmented Knothe function of time is proposed. D-InSAR technique has lower monitoring accuracy for large-scale deformation in the process of working face mining; so, after a period of time at the end of the mining of the working surface, the sinking rate slows down and the surface deformation stabilizes, and then the D-InSAR approach is used to obtain the surface’s LOS deformation, and then the vertical deformation, north–south deformation, and east–west deformation are predicted by DPIM. The deformation function equation is constructed. The particle swarm algorithm is utilized to reverse the parameters of the PIM and the improved segmented Knothe time function, which are predicted to acquire the dynamic three-dimensional deformation of the ground surface.

2. Theoretical Method

2.1. Dynamic Probability Integral Method

The random medium theory serves as the foundation for the probabilistic integral method’s mathematical theory, which considers the creation of rocks as an arbitrary medium and the surface movement and subsidence during rock mining as a random process [40]. The probabilistic integral method considers that the deformation of the mining face can be considered the deformation’s superposition brought on by the random movement of the unit random medium. The main idea is to divide the working surface of length L into n small rectangles, assigning each rectangle a corresponding time function value, thus realizing the anticipated deformation at any given time of any spot on the working surface. As shown in Figure 1, along the operating face’s mining direction, the length of the strike of the working face, L, minus the offset of the right boundary’s inflection point, S₃, and the offset of the left boundary’s inflection point, S₄, the length (L-S₃-S₄), is divided into small rectangles of the same shape according to the equal distances and numbered as 1, 2, 3……, n. This gives n rectangles with length $(L-S_3-S_4)/n$ and width D-S₁-S₂. The deformation value of each small rectangle is expected using the probabilistic integration method, and then the deformation value of n small rectangles can be superimposed to obtain the deformation value of the whole working surface.

In accordance with the probability integral method’s principle, the equation for calculating the surface subsidence profile of a semi-mining working face is as follows:

(1)$W(x)={\displaystyle \frac{W_0}{2}}[\mathrm{erf}({\displaystyle \frac{\sqrt{\pi }}{r}}x)+1]$

In the equation above, W₀ is the maximum value that can be reached by the sinking of the working surface under full mining. $W_0=qm\mathit{cos}\alpha$, q is the subsidence factor, α is the dip angle of the coal seam, and erf() is the error function. The sinking value at point A is as follows when the working face advances a distance l:

(2)$W(x,y,l)={\displaystyle \frac{1}{W_0}}(W(x)-W(x-l))(W(y)-W(y-D_0))$

In the above equation, D₀ is the calculated length of dip, which is calculated by projecting the actual length of the dip seam onto the horizontal plane:

(3)$D_0={\displaystyle \frac{(D-S_1-S_2)\mathit{sin}(\pi -\alpha -\theta )}{\mathit{sin}\theta }}$

The dynamic probabilistic integral model DPIM is obtained by combining Equation (2) with the time function. A single rectangle mining will have a dynamic effect on the working surface at point A. Point A’s subsidence value is shown below:

(4)$W(x,y,t)={\displaystyle \frac{1}{W_0}}(W(x)-W(x-l_t))(W(y)-W(y-D_0))\phi (t)$

When the working face mining advances at any given time, this provides the value of subsidence at any location on the surface:

(5)$W(x,y,t_i)={\displaystyle \frac{1}{W_0}}{\sum }_{i=1}^n\phi (t_i)(W(x)-W(x-l_{t_i}))(W(y)-W(y-D_0))$

Similarly, the surface point (x, y)’s dynamic horizontal movement values along the east–west and north–south directions at the moment t during the mining of the working face are as follows:

(6)$U_{\mathrm{ew}}(x,y,t_i,{\mu }_{\mathrm{ew}})={\displaystyle \frac{1}{W_0}}{\displaystyle \sum _{i=1}^n(U(x,t_i)W(y,t_i)\cos }{\mu }_{\mathrm{ew}}+U(y,t_i)W(x,t_i)\sin {\mu }_{\mathrm{ew}}$

(7)$U_{\mathrm{ew}}(x,y,t_i,{\mu }_{\mathrm{ew}})={\displaystyle \frac{1}{W_0}}{\displaystyle \sum _{i=1}^n(U(x,t_i)W(y,t_i)\cos }{\mu }_{\mathrm{ew}}+U(y,t_i)W(x,t_i)\sin {\mu }_{\mathrm{ew}}$

(8)$U(x)=b{W}_0{e^{-\pi ({\displaystyle \frac{x}{r}})}}^2$

In the formula above, U_ns is the working surface’s horizontal movement deformation in the north–south direction; U_ew is the working surface’s horizontal movement deformation in the east–west direction; b is the horizontal movement factor; μ_ns is the east–west angle to the projected coordinate system’s x-axis of the working face; and μ_ew is the north–south angle to the projected coordinate system’s x-axis of the working face.

2.2. Improved Time Function Model

The mine surface dynamic subsidence prediction model was fitted using the Knothe time function model, which was first proposed by Polish scholar Knothe in 1952. It is predicated mostly on the idea that the difference between the point’s eventual subsidence value and the subsidence value at this time in the subsidence process is related to the instantaneous sinking rate of a surface point at a given moment:

(9)$\psi \left(t\right)=1-\mathit{exp}(-ct)$

In the above equation, $\psi \left(t\right)$ is the time function. The time function’s fixed coefficient is denoted by c and t is the difference in time between the working face’s mining time and the anticipated time. An expression for the surface point’s dynamic subsidence value is as follows:

(10)$W\left(t\right)=W_0\times y\left(t\right)=W_0(1-exp(-ct\left)\right)$

W₀ is the surface point’s ultimate subsidence value in the equation above.

However, it is discovered that the theoretical surface subsidence process, as defined by the Knothe time function, is not entirely compatible with the actual subsidence process in the practice of dynamic subsidence prediction inside the mining region, which is mainly manifested in the fact that the rate of surface point sinking described by the time function is not consistent with the actual deformation process. The main thing is that the Knothe time function curve climbs rapidly to large values in a very short period of time, then rises slowly, gradually stabilizes, and eventually converges to a certain value. But the rule of surface deformation is that it does not rise very fast for a short time but reaches a maximum rate of subsidence in the intermediate process. Thus, Chang and other scholars raise the Knothe segmented time function, which divides the described surface deformation process into two phases: With the mining of the working face, the surface subsidence rate increases from 0 to the highest value in the first phase. Then, it gradually decreases from the maximum value to 0, which is expressed as:

(11)$\psi (t)=\left\{\begin{array}{l}{\psi }_1\left(t\right)=0.5\left[\mathit{exp}\right(-c(\tau -t))-\mathit{exp}(-c\tau \left)\right], 0<t\le \tau \\ {\psi }_2\left(t\right)=1-0.5\left[\mathit{exp}\right(-c(t-\tau ))+\mathit{exp}(-ct\left)\right], \tau <t\le T\end{array}\right.$

In the above equation, $\psi \left(t\right)$ is Knothe segmented time function, ${\psi }_1\left(t\right) \mathrm{a}\mathrm{n}\mathrm{d} {\psi }_2\left(t\right)$ are the values of the functions corresponding to the first and second stages at the moment t, c is a fixed coefficient of the time function, the time interval between the beginning of surface subsidence and its maximum rate is denoted by τ, and the surface subsidence time is denoted by T. The expressions for dynamic surface subsidence based on segmented time functions are given as follows:

(12)$W\left(t\right)=W_0\times \psi (t)=\left\{\begin{array}{l}{W}_0\times {\psi }_1(t), 0<t\le \tau \\ W_0\times {\psi }_2\left(t\right), \tau <t\le T\end{array}\right.$

When predicting surface deformation, the segmented time function increases the Knothe time function’s accuracy, but there is a discrepancy between the predicted subsidence boundary and the actual subsidence. Theoretically, when t is 0, the surface point’s subsidence value should be 0, that is, the time function value should be 0. However, when t = 0, the time function value of this segmented time function is not equal to 0. Meanwhile, theoretically, when τ moment is 0, the surface point’s rate of subsidence attains the maximum value, and the value of subsidence falls to half of the maximum value. At this time, the value of the time function at τ moment is 0.5 but in this segmented time function, when the τ value is determined, the time function’s value at the τ moment deviates from the theoretical value with the change in the c value. Specifically, as c gradually decreases, $\psi \left(\tau \right)$ gradually increases from the theoretical value of 0.5. Meanwhile, theoretically, with the extension of time, when the surface deformation is gradually stabilized, the time function value is theoretically 1 when the surface sinking ultimately reaches the theoretical maximum value. However, when varying the time function fixed coefficient c and the maximum sinking velocity moment τ, the final time function value does not converge to 1 but to a function value less than 1. Specifically, with the gradual decrease in the product of the maximum sinking velocity moment τ and the fixed coefficient c of the time function, the converged value’s deviance from the theoretical value of 1 gradually increases, and only if the value of the product of the maximum sinking velocity moment τ and the fixed coefficient c of the time function is greater than a specific amount, the final converged value of the time function will be close to the theoretical value of 1. These issues have a significant impact on the model’s applicability during mine surface deformation as well as the accuracy of the Knothe segmented time function in estimating the maximum subsidence value.

To address the above problem, solving the issue where the function value is not 0 while t is at 0 moments is the first step. After comparative calculations, it is found that at the moment t = 0, the functional value of $\psi \left(0\right)$ is $0.5\mathit{exp}(-c\tau )$, and the functional expression for the 0–τ time period is corrected to obtain the functional expression for the time of the first stage as:

(13)${\psi }_1(t)=0.5[\mathit{exp}(-c(\tau -t)-{\displaystyle \frac{(\tau -t)\mathit{exp}(-ct)}{\tau }})],0\le t\le \tau$

When correcting for the function value of the first stage, the function value of the second stage should theoretically be corrected equally, so the expression for the time function of the second stage can be obtained as:

(14)${\psi }_2(t)=1-0.5[\mathit{exp}(-c(\tau -t)-{\displaystyle \frac{(t-\tau )\mathit{exp}(-ct)}{\tau }})],\tau <t\le T$

The corrected time function expression satisfies that the time function value is 0 at the moment t = 0, and also satisfies that the time function value is 0.5 at the moment τ. However, with the extension of the time, the final value of the segmented time function is still not guaranteed to converge to 1. In order to deal with this problem, the second stage of the time function expression is corrected:

(15)${\psi }_2(t)=1-0.5[\mathit{exp}(-c(\tau -t)-{\displaystyle \frac{(t-\tau )\mathit{exp}(-ct)}{\tau }})]-0.5{\displaystyle \frac{(t-\tau )\mathit{exp}(-ct)}{\tau }},\tau <t\le T$

Simplifying Equations (13) and (15) obtains the improved Knothe segmented time function’s final expression as follows:

(16)$\psi (t)=\left\{\begin{array}{l}{\psi }_1(t)=0.5[{\displaystyle \frac{t}{\tau \mathit{exp}(c\tau )}}-{\displaystyle \frac{1-\mathit{exp}(ct)}{\mathit{exp}(c\tau )}}], 0\le t\le \tau \\ {\psi }_2\left(t\right)=1-0.5\mathit{exp}(c\tau -ct), \tau <t\le T\end{array}\right.$

The improved Knothe segmented time function corrects the problems of the Knothe time function and segmented Knothe time function, which is more in line with the law of the surface deformation process and closer to the ideal time function curve. As a result, the dynamic subsidence expression of surface deformation according to the improved Knothe segmented time function is:

(17)$W\left(t\right)=W_0\times \psi (t)=\left\{\begin{array}{l}{\psi }_1(t)=0.5W_0[{\displaystyle \frac{t}{\tau \mathit{exp}(c\tau )}}-{\displaystyle \frac{1-\mathit{exp}(ct)}{\mathit{exp}(c\tau )}}], 0\le t\le \tau \\ {\psi }_2\left(t\right)=W_0[1-0.5\mathit{exp}(c\tau -ct)] \tau <t\le T\end{array}\right.$

2.3. Constructing the DPIM-InSAR Equation

As shown in Figure 2, the LOS-oriented deformation of the surface monitored by the InSAR is the sum of the three-dimensional deformation projected into LOS from the vertical deformation, the north–south horizontally moving deformation, and the east–west horizontally moving deformation, as shown below:

(18)$D_{\mathrm{los}}=W\mathit{cos}\theta +U_{\mathrm{ns}}\mathit{sin}\phi \mathit{sin}\theta -U_{\mathrm{ew}}\mathit{cos}\phi \mathit{sin}\theta$

In the above equation, D_los is LOS deformation value, W is the vertical deformation value, U_ns is the north–south horizontal movement deformation value, φ is the satellite heading angle, U_ew is the east–west horizontal movement deformation value, and the satellite sensor signal’s angle of incidence with respect to the vertical direction is denoted by θ.

The mining subsidence prediction model cannot directly predict the LOS deformation of surface points but can predict the vertical deformation value, north–south horizontal movement deformation, and east–west horizontal movement deformation. Section 2.1 provides a detailed introduction to the prediction method and formula. Here, for the convenience of expression and simplification of the formula, a mapping function is used to express it. The expression is as follows:

(19)$\left\{\begin{array}{l}W(x,y,t)=f_W(x,y,t)\\ U_{\mathrm{ns}}(x,y,t,\mu )={f_U}_{\mathrm{ns}}(x,y,t,\mu )\\ U_{\mathrm{ew}}(x,y,t,\mu )={f_U}_{\mathrm{ew}}(x,y,t,\mu )\end{array}\right.$

In the above equations, $W(x,y,t)$ is the expression of subsidence function, $f_W(x,y,t)$ is the mapping relation of subsidence function, $U_{\mathrm{ns}}(x,y,t,\mu )$ is the expression of the north–south horizontal movement deformation function, ${f_U}_{\mathrm{ns}}(x,y,t,\mu )$ is the mapping relation of the north–south horizontal movement deformation function, $U_{\mathrm{ew}}(x,y,t,\mu )$ is the expression of the east–west horizontal movement deformation function, and ${f_U}_{\mathrm{ew}}(x,y,t,\mu )$ is the mapping relation of the east–west horizontal movement deformation function.

The geometric projection relationship between the three-dimensional deformations and the LOS-oriented deformation that InSAR monitors serves as the foundation for the DPIM-InSAR equations:

(20)$D_{\mathrm{los}}(x,y,t)={\left[\begin{array}{l}\mathit{cos}\theta \\ \mathit{sin}\theta \mathit{sin}\phi \\ -\mathit{sin}\theta \mathit{cos}\phi \end{array}\right]}^T\left[\begin{array}{l}{f}_W(x,y,t)\\ f_{U_{\mathrm{ew}}}(x,y,t,{\mu }_{\mathrm{ew}})\\ f_{U_{\mathrm{ns}}}(x,y,t,{\mu }_{\mathrm{ns}})\end{array}\right]$

This leads to surface point LOS deformation values at a given time period:

(21)$D_{\mathrm{los}}(x,y,t_i)=\left[\begin{array}{l}\mathit{cos}\theta \\ \mathit{sin}\theta \mathit{sin}\phi \\ -\mathit{sin}\theta \mathit{cos}\phi \end{array}\right]\left[\begin{array}{l}{f}_W(x,y,t_i)\\ f_{U_{\mathrm{ew}}}(x,y,t,{\mu }_{\mathrm{ew}})\\ f_{U_{\mathrm{ns}}}(x,y,t,{\mu }_{\mathrm{ns}})\end{array}\right]$

(22)$D_{\mathrm{los}}(x,y,t_{i-1})=\left[\begin{array}{l}\mathit{cos}\theta \\ \mathit{sin}\theta \mathit{sin}\phi \\ -\mathit{sin}\theta \mathit{cos}\phi \end{array}\right]\left[\begin{array}{l}{f}_W(x,y,t_{i-1})\\ f_{U_{\mathrm{ew}}}(x,y,t,{\mu }_{\mathrm{ew}})\\ f_{U_{\mathrm{ns}}}(x,y,t,{\mu }_{\mathrm{ns}})\end{array}\right]$

(23)$\Delta D_{\mathrm{los}}(x,y,t_i-t_{i-1})=D_{\mathrm{los}}(x,y,t_i)-D_{\mathrm{los}}(x,y,t_{i-1})$

In the above equation, D_los(x, y, t_i) is the surface point (x, y) LOS deformation value at the moment of t_i, D_los(x, y, t_i−1) is the surface point (x, y) LOS deformation value at the moment of t_i−1, and ΔD_los(x, y, t_i − t_i−1) is the surface point (x, y) LOS deformation value for the period of time from the moment of t_i to the moment of t_i−1.

2.4. Parameter Inversion for Particle Swarm Algorithm

When the improved function is utilized to predict the dynamic deformation, the required parameters are the probability integration parameters [q, b, tanβ, θ, S₁, S₂, S₃, S₄] and the improved segmented time function [c, τ].

The PSO algorithm is based on the phenomenon of birds feeding in flocks and its basic idea is to provide each person’s and the group’s roles their due. While searching for the optimal solution, the individual shares the acquired information with the group, compares the values of each individual among the group to determine the optimal solution, and repeats the cycle until the best optimal solution is found. The specific process is shown in Figure 3.

According to available LOS deformation results from D-InSAR monitoring, the eigenvalues are extracted to the inversion equation and the required parameters are inverted using the particle swarm algorithm in the following steps:

1. Initialize particles. Set a certain number of particles to form a particle swarm and randomly assign values to this swarm of particles to obtain the current global best optimal solution.

2. Iterating particle values. Each particle updates itself by tracking two poles (pbest, gbest) and the particle uses Equations (24) and (25) to update its position and velocity after determining the two poles.

(24)$v_i=v_i+c_1\times rand\left(\right)\times (pbes{t}_i-x_i)+c_2\times rand\left(\right)\times (gbes{t}_i-x_i)$

(25)$x_i=x_i+v_i$

In the above equations, v_i is the velocity of the particle, rand() is a random number between 0 and 1, c₁ and c₂ are the learning factors, and x_i is the current position of the particle.

• 3.. Output optimal solution. Compare the particle extremes and output the global optimal solution.

3. Simulation Experiment

3.1. Simulated Experimental Data

To verify that the ISK-DPIM-InSAR model is accurate, the three-dimensional deformation field of the mine’s working surface at a certain time period was simulated using the dynamic probability integration method. With a strike length of 500 m, a dip length of 200 m, a coal seam thickness of 4 m, a mining depth of 350 m, a coal seam dip angle of 7 degrees, and an azimuth angle of 180 degrees along the strike, the simulated working face is depicted in Figure 4. According to the mining geological conditions, the probability integral parameters were simulated with a horizontal motion coefficient b of 0.3, a subsidence coefficient q of 0.7, a mining influence propagation angle θ of 89°, a main influence angle tangent tanβ of 1.5, and the uphill, downhill, left boundary, and bounded inflection point offsets are all 20 percent from S1, S2, S3, and S4. In the improved Knothe segmented time function, the fixed parameter c of the time function is 0.02, and the maximum sinking velocity moment τ is 200. All the above parameters are the true values of the parameters for this simulated workface and the values inverted by the DPIM-InSAR model are the observed values of the parameters. To facilitate validation of the accuracy of the proposed ISK-DPIM-InSAR model, two observation lines are set up on the simulated working surface, as shown in Figure 4, and monitoring points are set up every 10 m along the direction AB, which is parallel to the strike, and the direction CD, which is parallel to the tendency, and the true values and observed values of the monitoring points are compared, so as to analyze the accuracy of the DPIM-InSAR model.

Based on the DPIM-InSAR model with improved Knothe segmentation time function, the surface deformation results when the working face mining time is 300 days and 350 days are expected and the surface subsidence results, the surface east–west horizontal movement deformation results, and the surface north–south horizontal movement deformation results can be obtained in the time period between 300 and 350 days, as shown in Figure 5. The results obtained from the DPIM model are assumed to be sedimentation results from D-InSAR monitoring in this time period. It is also assumed that the satellite radar has an incidence angle of 38.7° and an azimuth of 350° and that the revisit period of the satellite time is 12 days.

According to Figure 5, the maximum vertical deformation of the simulated working surface is −405 mm during the period from 300 days to 350 days from the mining moment. There is symmetry in the horizontal movement deformation of the surface in both the east–west and north–south directions. The maximum value of horizontal movement deformation in the east–west direction is −162 mm. The maximum value of horizontal movement deformation in the north–south direction is −72 mm. The maximum value of the 3D deformation projection of the simulated working surface into the LOS direction is −331 mm.

The values of deformation W in the vertical direction, horizontal moving deformation U_ns along the north–south direction, horizontal moving deformation U_ew along the east–west direction, and LOS deformation of the simulated working face predicted by the ISK-DPIM model were combined with the DPIM-InSAR equations of Equations (20) and (23), and a particle swarm algorithm is utilized to invert the parameters in the DPIM model [q, b, tan β, θ, S₁, S₂, S₃, S₄] and [c, τ].

3.2. Simulation Experiment Results and Accuracy Analysis

From the simulated LOS-oriented deformation of the working face, 80 observation point deformation values were randomly selected uniformly at the subsidence basin’s edge as the deformation data for the ISK-DPIM-InSAR model inversion parameters, and the point locations are shown in Figure 6. Utilizing the 80 distinct observation stations’ LOS-oriented deformation values, combined with the DPIM-InSAR equations, the particle swarm algorithm is used to invert the model parameters [q, b, tanβ, θ, S₁, S₂, S₃, S₄] and [c, τ], and the inversion results and the parameter truth values are as displayed in Table 1.

The inverted parameter values were brought into the DPIM model in order to anticipate the 3D deformation of the simulated working face at the moment t = 350, i.e., surface deformation in the vertical direction and horizontal movement deformations in the east–west direction and north–south direction are as seen in Figure 7. The predicted results are also compared and analyzed with the true values, as shown in Figure 8.

For a comparison of the predicted values from the ISK-DPIM-InSAR model and the true values from the simulated workings, observation points were set up along the simulated working surface along the observation line of the AB section and the observation line of the CD section at an interval of 10 m. The AB section observation had 101 observation points, while the CD section observation line had 66 observation points, for a total of 167 spots.

The subsidence of the observation line in section AB, the subsidence of the observation line in section CD, and the horizontal movement deformations along the east–west direction and the north–south direction predicted by the inversion parameters of the model and the real value parameters of the simulated working surface are compared, as shown in Figure 8. It goes to show that the deformation trends of the predicted model inversion parameters and the predicted results of the true value are basically the same, although the maximum subsidence value predicted by the model inversion parameters is slightly larger than the true value of the simulated working surface; the maximum absolute error is 87 mm, which is only 5.3% of the maximum subsidence.

The highest possible absolute error in subsidence of the observation line in section AB is 78 mm, which is only 4.8% of the maximum subsidence, and the highest possible absolute error in subsidence of the observation line in section CD is 66 mm, which is only 4.0% of the maximum subsidence. The predicted horizontal deformation values of the model inversion parameters are generally consistent with the true values: the observation line’s horizontal movement in section AB has a maximum absolute inaccuracy of 36 mm and the observation line’s horizontal movement in section BC has a maximum absolute inaccuracy of 31 mm. Overall the error between the deformation values predicted by the proposed model and the true values is small, so the ISK-DPIM-InSAR model can effectively predict the 3D deformation of mining subsidence in mining areas. The predicted accuracy of the model was further verified by calculating the root mean square error (RMSE) of the 3D deformation for the predicted and true values, as shown in Table 2. The RMSE of the subsidence of the observation line in the AB section is 40.10 mm, which is only 2.8% of the maximum subsidence, and that of the observation line in the CD section is 30.87 mm, which is only 1.8% of the maximum subsidence. The root mean square error is only 10.00 mm for horizontally moving deformations in the north–south direction and 18.59 mm for horizontally moving deformations in the east–west direction. The results show that the proposed ISK-DPIM-InSAR model can accurately predict the three-dimensional surface deformation caused by mining subsidence in mining areas.

4. ISK-DPIM-InSAR Modeling for Mine Monitoring

4.1. Materials and Methods

The study area is a working face situated in the Huaibei mining area, Anhui Province. There is a strike length of 890 m, a dip length of 120 m, an average coal seam thickness of 7.27 m, a dip angle of 17°, and an average mining depth of 592 m on this working face. The mining date of the face is 12 December 2006, the stopping date is 25 February 2008, and the mining time lasts 431 days. The working face is surrounded by farmland, ditches, and villages, which will be affected by mine mining, and the geographical location is shown in Figure 9.

The SAR images required by the D-InSAR technology are the two-scene images acquired by the ALOS PALSAR satellite sensor, which transmits L-band radar with multiple imaging modes, which is more effective for monitoring the mining area. Considering the working period for face mining, the primary image data used are from 10 December 2007 to 12 January 2008, which correspond to the 354th to 399th days of the working face mining time, respectively. Table 3 details parameters such as the wavelength, imaging mode, and polarization mode of the SAR image data used.

4.2. D-InSAR Results and Parameter Inversion

After processing by D-InSAR technology, the surface subsidence results of the working face from 354 days to 399 days can be obtained, as shown in Figure 10. Excluding the area with serious decoherence, the maximum deformation value of D-InSAR in LOS direction is −180 mm.

In this paper, the D-InSAR technique with two images + external DEM method is used to monitor the surface deformation, and the processing flow is divided into eight steps as follows:

1. Data preprocessing

The original SAR image data are focused into the single-view complex format (SLC), and the complex data include the modulus of the complex number and the phase of the complex number, and the modulus of the complex number represents the SAR signal strength. The complex phase represents the time delay of the SAR echo signal.

• 2.. Main and secondary image alignment and resampling

The image interference requires that the main and auxiliary images should have the same radar coordinate image space, and it is generally considered that the image after the date is the auxiliary image, and the auxiliary image is resampled to the radar coordinate image space of the main image. Finally, the main and secondary images are coarsely aligned and finely aligned.

• 3.. Interference between primary and secondary images

After the second step, the main and auxiliary images already have the same name image element, the interference phase can be obtained by conjugating the complex data with the same name phase, and the interfered data are still complex data. The interference phase includes the surface deformation phase and other influence phases, which need to be eliminated in turn.

• 4.. Removing the reference ellipsoid phase

The geometric relationship between the satellite’s spatial attitude and the Earth’s reference ellipsoid at the moment of shooting the primary and secondary images forms the reference ellipsoid phase, and it is necessary to obtain the relevant parameters of the satellite at the moment of shooting and construct a model to remove the reference ellipsoid phase.

• 5.. Removal of terrain phase

To obtain the surface deformation phase of the study area through external DEM data, firstly, the geographic coordinates of the DEM data should be converted to under the SAR image coordinate system, and the converted data should be aligned to obtain the elevation value of the same name pixel of the primary and secondary images, so that the surface deformation phase can be removed.

• 6.. Interferogram filtering and masking

The atmospheric delay phase and noise phase can be filtered by means of the interferogram to reduce their influence on the D-InSAR results, and the Goldstein adaptive filtering method is used here. At the same time, the coherence of the interferogram is calculated, and the regions with serious incoherence are masked.

• 7.. Phase unwrapping

The phase of the interferogram is entangled in a certain interval, which is not the true phase value, and the phase is de-entangled by the minimum cost flow method MCF to recover the true phase of the interferogram.

• 8.. Coordinate conversion

The interference results after processing in the above steps are still in the SAR image coordinate system and need to be converted to the geographic coordinate system.

Based on the ISK-DPIM-InSAR equation constructed in Section 2.3, the LOS deformation value monitored by D-InSAR needs to be extracted as the left input value of the equation, and the right side of the equation is composed of DPIM model functions with parameters. In theory, the left and right sides of the equation are equal, so that the parameters can be obtained by particle swarm optimization. Based on the geological conditions around the working face and the results of D-InSAR monitoring, the LOS deformation values of 50 feature points were extracted around the subsidence basin of the working face.

The ISK-DPIM-InSAR equations are constructed from the extracted feature point coordinates and deformation values, the initial particle swarm size is set to 100, the iteration is 100 times, the learning factors c₁ and c₂ are both 2.0, and the initial value of the particle swarm is randomly assigned. The optimal solution of the inversion parameters is output after iteration, as shown in Table 4. The model parameters based on the inversion can realize the three-dimensional deformation prediction of vertical deformation W, north–south horizontal movement deformation U_ns, and east–west horizontal movement deformation U_ew of the working face.

4.3. Analysis of Experimental Results and Accuracy

Using the dynamic prediction parameters of mining subsidence obtained by the inversion described in Section 4.2, the three-dimensional deformation dynamic prediction of vertical deformation W, north–south horizontal movement deformation U_ns, and east–west horizontal movement deformation U_ew in the mining process of the working face is carried out. Simultaneously, in order to verify the accuracy of the ISK-DPIM-InSAR model, 7 leveling data of the working face at different stages of the mining process were selected to compare with the predicted deformation results of the model.

The vertical deformation W results predicted by the DPIM-InSAR model (as shown in Figure 11) show that with the continuous advancement of working face mining, the vertical deformation of the surface gradually increases, forming a trend of progressively declining from the basin’s center to its two edges, the subsidence basin bias workings downhill show symmetry, and the maximum surface subsidence deformation reaches −1558 mm.

The U_ns deformation results of the north–south horizontal movement (as shown in Figure 12) and the U_ew deformation results of the east–west horizontal movement (as shown in Figure 13) predicted by the ISK-DPIM-InSAR model show that the U_ns deformation value of the working face is between −624 mm and −624 mm, and the U_ew deformation value is between −197 mm and −157 mm. With the advancement of the working face mining process, the surface horizontal movement deformation gradually increases and presents symmetry.

A total of 71 leveling points were laid on the strike observation line of the working face, and 55 leveling points were laid on the dip observation line of the working face. The corresponding point level data are extracted and compared with the predicted value of the ISK-DPIM-InSAR model (as shown in Figure 14). The results show that the predicted value of the ISK-DPIM-InSAR model is basically consistent with the trend of strive subsidence and dip subsidence of leveling measured data in the process of mining. The root mean square error of strive subsidence is 45 mm, and the root mean square error of dip subsidence is 32 mm. The ISK-DPIM-InSAR model can effectively predict the 3D surface deformation of the mining area.

5. Discussion and Conclusions

5.1. Discussion

During the mining process, the Probability Integral Model (PIM) has been widely applied. This model allows operators to predict the final subsidence during the excavation of the working face, providing an important basis for ensuring safety in mining operations. However, the model requires numerous parameters, with seven basic ones alone. These parameters include subsidence factors, mining area, and seam thickness, among others. Changes in these factors directly influence the model’s accuracy and its practical applicability. While the traditional Probability Integral Model can reflect the general trend of subsidence to a certain extent, its calculation process relies heavily on large amounts of field data and historical experience, making the parameter acquisition process relatively cumbersome and limited in precision. To address this issue, this paper proposes a method that combines InSAR (Interferometric Synthetic Aperture Radar) with the Probability Integral Model, significantly improving the efficiency of parameter acquisition. Furthermore, by applying machine learning techniques to invert the parameters, the accuracy of the parameter estimation is improved, making the model’s predictions more aligned with the actual conditions. The introduction of the dynamic Probability Integral Model further enhances the safety of mining operations, as it allows for an understanding of subsidence at different time intervals. This enables timely adjustments to safety measures and problem-solving after a collapse of the working face based on the observed subsidence.

Overall, the main factors influencing the accuracy can be summarized into two aspects. The first factor is the precision of the InSAR monitoring results. InSAR monitoring accuracy depends on the quality of satellite data and the frequency band used. According to the model’s predicted results, the vertical deformation at the mining face can reach a maximum of −1558 mm, which exceeds the range of InSAR monitoring. Although the images used in this study are from the ALOS satellite, which can monitor larger deformation values than C-band satellite imagery, it still cannot accurately monitor deformation in the central region of the working face. Therefore, we can only use the data from the closest points to the central area that fall within the monitorable range to invert the parameters. Ideally, if the deformation data in the central area could be monitored more precisely, the accuracy of the InSAR inversion would be greatly improved. At the same time, parameter inversion does not require many deformation data periods. Therefore, compared to the time-series InSAR monitoring method, the D-InSAR method can meet the model’s requirements. Additionally, since the D-InSAR method has a lower computational load, it saves a significant amount of image data cost, making it more cost-effective in practical applications.

The second factor depends on the accuracy of the DPIM. First of all, although the PIM model is currently the most extensively utilized mining subsidence prediction model in the mining area, Liu and the scholars who proposed the model have published an article indicating that there are certain limitations in adapting the PIM itself to the mining area [40]. The mathematical theory of the PIM is a random medium, which itself does not conform to the complex rock medium of the mine. The mining subsidence model suitable for the mine can be considered from the perspective of physical mechanics. Secondly, the time function model combined with the PIM is also an important factor affecting the accuracy. The advantages of the ISK time function model over the other previous time function models are explained in depth in this paper’s Section 2.2, and the improved model enhances the accuracy of mining subsidence prediction.

5.2. Conclusions

This study focuses on the three-dimensional deformation of subsidence caused by mining activities and proposes an ISK-DPIM-InSAR model for dynamic subsidence prediction, which integrates InSAR technology and an improved piecewise Knothe time function-based Probability Integral Method (PIM). Firstly, this study develops geometric relationship equations based on D-InSAR technology to monitor different deformation directions, including Line of Sight (LOS) deformation, vertical deformation, and horizontal displacements along the north–south and east–west directions. Secondly, a simulated mining face is designed to verify the feasibility of the model. The true values of the inverted parameters are used to predict the LOS deformation of the simulated mining face. To optimize the model parameters, Particle Swarm Optimization (PSO) is employed to iteratively seek the optimal solution for the inversion equations. These optimized parameters are then input into the Dynamic Probability Integral Model (DPIM) for dynamic three-dimensional subsidence prediction of the mining face. The model’s predictions are compared with the true three-dimensional deformation of the simulated mining face to assess the accuracy and reliability of the results. The study’s findings demonstrate that the ISK-DPIM-InSAR model effectively predicts the dynamic three-dimensional deformation during mining subsidence. Not only does it enhance the accuracy of deformation prediction, it also provides real-time updates on subsidence changes during the mining process, offering valuable insights for mine safety management.

In this paper, the ISK-DPIM-InSAR model, which has been verified to be feasible, is applied to an example in a mining area. Firstly, the LOS deformation was obtained by processing two views of L-band SAR pictures in the Huaibei mining area by D-InSAR technology. Together with the corresponding vertically orientated deformations and U_ns, U_es horizontal moving deformations, it forms the geometrical equations according to the projected functions. The predicted parameters are obtained by particle swarm inversion and substituted for the ISK-DPIM-InSAR model to predict the dynamic 3D deformation, and the deformation values are also compared with the level values. The proposed model can accurately monitor the three-dimensional deformation in the mining process, which effectively makes up for the limitation that only the LOS deformation can be captured using InSAR technology, and solves the problem that finding the model’s parameters to forecast the dynamics of mining subsidence is challenging.

Author Contributions

Conceptualization, S.W. and G.L.; methodology, S.W.; software, Z.S.; validation, M.L., Y.C. and M.W.; formal analysis, S.W.; investigation, Y.C.; resources, K.Y.; data curation, Z.S.; writing—original draft preparation, S.W.; writing—review and editing, S.W.; visualization, S.W.; supervision, K.Y.; project administration, G.L.; funding acquisition, G.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

SRTM DEM data were downloaded from the EARTHDATA (https://search.earthdata.nasa.gov/search?q=SRTM, accessed on 2 March 2022).

Acknowledgments

The authors thank Keming Yang, China University of Mining and Technology, Beijing, for his guidance. The authors thank the leadership of the Center for Geophysical Survey for their support. The authors thank the remote sensing team for their contributions.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Dynamic projection model.

Figure 2. InSAR side-view imaging geometric relationships and angular parameters.

Figure 3. Particle swarm algorithm flow.

Figure 4. Simulated working surface plan.

Figure 5. Simulated work face of surface deformation between 300 and 350 days. (a) Simulated of work face of surface subsidence between 300 days and 350 days; (b) simulated work face of horizontal ground movement deformation along the east–west direction between 300 days and 350 days; (c) simulated work face of horizontal ground movement deformation along the north–south direction between 300 days and 350 days; and (d) simulated work face surface deformation values along the LOS between 300 days and 350 days.

Figure 6. Simulation of LOS deformation results of extracted feature points.

Figure 7. DPIM-InSAR model predicts 3D surface deformation values for simulated working surfaces. (a) Simulated working face of vertical surface deformation on day 350. (b) Simulated working face of horizontal ground movement deformation along the north–south direction on day 350. (c) Simulated working face of horizontal ground movement deformation along the east–west direction on day 350.

Figure 8. Comparison of predicted and true values for simulated workface. (a) Comparison of subsidence along the AB observation line. (b) Comparison of subsidence along the CD observation line. (c) Comparison of horizontal movement deformation along east–west direction. (d) Comparison of horizontal movement deformation along the north–south direction.

Figure 9. Location of the study area.

Figure 10. LOS-oriented deformation results and characteristic points monitored by D-InSAR.

Figure 11. Predicted vertical deformation results for DPIM-InSAR model inversion parameters.

Figure 12. DPIM-InSAR model inversion parameters predicted north–south horizontal movement deformation results.

Figure 13. DPIM-InSAR model inversion parameters predicted east–west horizontal movement deformation results.

Figure 14. Comparison of DPIM-InSAR model projections and level data.

Table 1

Comparison of true and inverse values of parameters.

Parameters	q	b	tan β	θ (°)	S1 (m)	S2 (m)	S3 (m)	S4 (m)	c	τ (d)
Truth value	0.7	0.3	1.5	89	20	20	20	20	0.02	200
Inversion value	0.73	0.28	1.53	89.4	15.7	15.9	19.5	19.3	0.0195	195

Table 2

RMSE of the predicted 3D deformation of the simulated workface.

	Segment AB Subsidence	Segment CD Subsidence	North–South Horizontal Deformation	East–West Horizontal Deformation
RMSE (mm)	40.10	30.87	10.00	18.59

Table 3

SAR image specific parameters.

Platform	Date	Wavelength(cm)	Mode	PolarizationMode	Data Type	Average Incidence Angle (°)	Spatial Resolution (m)		Satellite Posture
							Direction	Azimuthal
ALOS PALSAR	10 December 200725 January 2008	23.6	FBS	HH	SLC	38.7	7.5	3.2	elevated orbit right vision

Table 4

Projected parameters for the DPIM-InSAR model inversion of the workface.

Parameter	q	b	tan β	θ (°)	S1 (m)	S2 (m)	S3 (m)	S4 (m)	c	τ
Inversion value	0.83	0.38	1.55	89.8	−21.7	20.1	57.8	57.9	0.0145	230

References

1. Dong, L.; Wang, C.; Tang, Y.; Tang, F.; Zhang, H.; Wang, J.; Duan, W. Time Series InSAR Three-Dimensional Displacement Inversion Model of Coal Mining Areas Based on Symmetrical Features of Mining Subsidence. Remote Sens.; 2021; 13, 2143. [DOI: https://dx.doi.org/10.3390/rs13112143]

2. Jiang, K.; Yang, K.; Zhang, Y.; Li, Y.; Li, T.; Zhao, X. An Extraction Method for Large Gradient Three-Dimensional Displacements of Mining Areas Using Single-Track InSAR, Boltzmann Function, and Subsidence Characteristics. Remote Sens.; 2023; 15, 2946. [DOI: https://dx.doi.org/10.3390/rs15112946]

3. Xing, X.; Zhang, T.; Chen, L.; Yang, Z.; Liu, X.; Peng, W.; Yuan, Z. InSAR Modeling and Deformation Prediction for Salt Solution Mining Using a Novel CT-PIM Function. Remote Sens.; 2022; 14, 842. [DOI: https://dx.doi.org/10.3390/rs14040842]

4. Li, Y.-X.; Yang, K.-M.; Zhang, J.-H.; Hou, Z.-X.; Wang, S.; Ding, X.-M. Research on time series InSAR monitoring method for multiple types of surface deformation in mining area. Nat. Hazards; 2022; 114, pp. 2479-2508. [DOI: https://dx.doi.org/10.1007/s11069-022-05476-8]

5. Zhu, J.; Yang, Z.; Li, Z. Recent progress in retrieving and predicting mining-induced 3D displace-ments using InSAR. Acta Geod. Cartogr. Sin.; 2019; 48, pp. 135-144. [DOI: https://dx.doi.org/10.1111/1755-6724.14223]

6. Chen, H.; Zhao, C.; Tomas, R.; Chen, L.; Yang, C.; Zhang, Y. Retrieving the Kinematic Process of Repeated-Mining-Induced Landslides by Fusing SAR/InSAR Displacement, Logistic Model, and Probability Integral Method. Remote Sens.; 2023; 15, 3145. [DOI: https://dx.doi.org/10.3390/rs15123145]

7. Xie, Y.; Bagan, H.; Tan, L.; Te, T.; Damdinsuren, A.; Wang, Q. Time-Series Analysis of Mining-Induced Subsidence in the Arid Region of Mongolia Based on SBAS-InSAR. Remote Sens.; 2024; 16, 2166. [DOI: https://dx.doi.org/10.3390/rs16122166]

8. Xu, Y.; Li, T.; Tang, X.; Zhang, X.; Fan, H.; Wang, Y. Research on the Applicability of DInSAR, Stacking-InSAR and SBAS-InSAR for Mining Region Subsidence Detection in the Datong Coalfield. Remote Sens.; 2022; 14, 3314. [DOI: https://dx.doi.org/10.3390/rs14143314]

9. Shi, J.; Yang, Z.; Wu, L.; Qiao, S. Improving the Robustness of the MTI-Estimated Mining-Induced 3D Time-Series Displacements with a Logistic Model. Remote Sens.; 2021; 13, 3782. [DOI: https://dx.doi.org/10.3390/rs13183782]

10. Shi, M.; Yang, H.; Wang, B.; Peng, J.; Gao, Z.; Zhang, B. Improving Boundary Constraint of Probability Integral Method in SBAS-InSAR for Deformation Monitoring in Mining Areas. Remote Sens.; 2021; 13, 1497. [DOI: https://dx.doi.org/10.3390/rs13081497]

11. Yang, Z.; Li, Z.; Zhu, J.; Yi, H.; Hu, J.; Feng, G. Deriving Dynamic Subsidence of Coal Mining Areas Using InSAR and Logistic Model. Remote Sens.; 2017; 9, 125. [DOI: https://dx.doi.org/10.3390/rs9020125]

12. Chen, Y.; Dong, X.; Qi, Y.; Huang, P.; Sun, W.; Xu, W.; Tan, W.; Li, X.; Liu, X. Integration of DInSAR-PS-Stacking and SBAS-PS-InSAR Methods to Monitor Mining-Related Surface Subsidence. Remote Sens.; 2023; 15, 2691. [DOI: https://dx.doi.org/10.3390/rs15102691]

13. Yang, Z.; Zhu, J.; Xie, J.; Li, Z.; Wu, L.; Ma, Z. Resolving 3-D Mining Displacements from Multi-Track InSAR by Incorporating with a Prior Model: The Dynamic Changes and Adaptive Estimation of the Model Parameters. IEEE Trans. Geosci. Remote Sens.; 2022; 60, 4504610. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3093058]

14. Hu, J.; Yan, Y.; Dai, H.; He, X.; Lv, B.; Han, M.; Zhu, Y.; Zhang, Y. Prediction Method for Dynamic Subsidence Basin in Mining Area Based on SBAS-InSAR and Time Function. Remote Sens.; 2024; 16, 1938. [DOI: https://dx.doi.org/10.3390/rs16111938]

15. Yang, Z.; Li, Z.; Zhu, J.; Wang, Y.; Wu, L. Use of SAR/InSAR in Mining Deformation Monitoring, Parameter Inversion, and Forward Predictions: A Review. IEEE Geosci. Remote Sens. Mag.; 2020; 8, pp. 71-90. [DOI: https://dx.doi.org/10.1109/MGRS.2019.2954824]

16. Huang, C.J.; Xia, H.M.; Hu, J.Y. Surface Deformation Monitoring in Coal Mine Area Based on PSI. IEEE Access; 2019; 7, pp. 29672-29678. [DOI: https://dx.doi.org/10.1109/ACCESS.2019.2900258]

17. Blachowski, J.; Jirankova, E.; Lazecky, M.; Kadlecik, P.; Milczarek, W. Application of satellite radar interferometry (PSINSAR) in analysis of secondary surface deformations in mining areas case studies from czech republic and poland. Acta Geodyn. Geomater.; 2018; 15, pp. 173-185. [DOI: https://dx.doi.org/10.13168/AGG.2018.0013]

18. Oliveira, S.C.; Zezere, J.L.; Catalao, J.; Nico, G. The contribution of PSInSAR interferometry to landslide hazard in weak rock-dominated areas. Landslides; 2015; 12, pp. 703-719. [DOI: https://dx.doi.org/10.1007/s10346-014-0522-9]

19. Monika,; Govil, H.; Chatterjee, R.S.; Bhaumik, P.; Vishwakarma, N. Deformation monitoring of Surakachhar underground coal mines of Korba, India using SAR interferometry. Adv. Space Res.; 2022; 70, pp. 3905-3916. [DOI: https://dx.doi.org/10.1016/j.asr.2022.05.018]

20. Li, T.; Zhang, H.; Fan, H.; Zheng, C.; Liu, J. Position Inversion of Goafs in Deep Coal Seams Based on DS-InSAR Data and the Probability Integral Methods. Remote Sens.; 2021; 13, 2898. [DOI: https://dx.doi.org/10.3390/rs13152898]

21. Yang, Z.; Xu, B.; Li, Z.; Wu, L.; Zhu, J. Prediction of Mining-Induced Kinematic 3-D Displacements from InSAR Using a Weibull Model and a Kalman Filter. IEEE Trans. Geosci. Remote Sens.; 2022; 60, 4500912. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3055854]

22. Liu, H.; Yuan, M.; Li, M.; Li, B.; Chen, N.; Wang, J.; Li, X.; Wu, X. TDFPI: A Three-Dimensional and Full Parameter Inversion Model and Its Application for Building Damage Assessment in Guotun Coal Mining Areas, Shandong, China. Remote Sens.; 2024; 16, 698. [DOI: https://dx.doi.org/10.3390/rs16040698]

23. Yang, Z.; Li, Z.; Zhu, J.; Feng, G.; Wang, Q.; Hu, J.; Wang, C. Deriving time-series three-dimensional displacements of mining areas from a single-geometry InSAR dataset. J. Geod.; 2018; 92, pp. 529-544. [DOI: https://dx.doi.org/10.1007/s00190-017-1079-x]

24. Li, Z.W.; Yang, Z.F.; Zhu, J.J.; Hu, J.; Wang, Y.J.; Li, P.X.; Chen, G.L. Retrieving three-dimensional displacement fields of mining areas from a single InSAR pair. J. Geod.; 2015; 89, pp. 17-32. [DOI: https://dx.doi.org/10.1007/s00190-014-0757-1]

25. Ding, X.; Yang, K.; Wang, S.; Hou, Z.; Li, Y.; Li, Y.; Zhao, H. Time series monitoring and prediction of coal mining subsidence based on multitemporal InSAR technology and GSM-HW model. J. Appl. Remote Sens.; 2022; 16, 038505. [DOI: https://dx.doi.org/10.1117/1.JRS.16.038505]

26. Jiang, C.; Wang, L.; Yu, X. Retrieving 3D Large Gradient Deformation Induced to Mining Subsidence Based on Fusion of Boltzmann Prediction Model and Single-Track InSAR Earth Observation Technology. IEEE Access; 2021; 9, pp. 87156-87172. [DOI: https://dx.doi.org/10.1109/ACCESS.2021.3089160]

27. Yang, Z.; Li, Z.; Zhu, J.; Preusse, A.; Yi, H.; Hu, J.; Feng, G.; Papst, M. Retrieving 3-D Large Displacements of Mining Areas from a Single Amplitude Pair of SAR Using Offset Tracking. Remote Sens.; 2017; 9, 338. [DOI: https://dx.doi.org/10.3390/rs9040338]

28. Wang, R.; Wu, K.; He, Q.; He, Y.; Gu, Y.; Wu, S. A Novel Method of Monitoring Surface Subsidence Law Based on Probability Integral Model Combined with Active and Passive Remote Sensing Data. Remote Sens.; 2022; 14, 299. [DOI: https://dx.doi.org/10.3390/rs14020299]

29. Chen, Y.; Yu, S.; Tao, Q.; Liu, G.; Wang, L.; Wang, F. Accuracy Verification and Correction of D-InSAR and SBAS-InSAR in Monitoring Mining Surface Subsidence. Remote Sens.; 2021; 13, 4365. [DOI: https://dx.doi.org/10.3390/rs13214365]

30. Yao, S.Y.; Yang, K.M.; Dong, X.L.; Shi, X.Y.; Wang, J. Large-gradient deformation monitoring and parameter inversion in a mining area using a method combining a dynamic prediction model and InSAR. Remote Sens. Lett.; 2021; 12, pp. 838-847. [DOI: https://dx.doi.org/10.1080/2150704X.2021.1942583]

31. Wang, J.; Yan, L.; Yang, K.M.; Tang, W.; Xie, H.; Yao, S.Y.; Xu, Z.H.; Yang, J.B. Deriving Mining-Induced 3-D Deformations at Any Moment and Assessing Building Damage by Integrating Single InSAR Interferogram and Gompertz Probability Integral Model (SII-GPIM). IEEE Trans. Geosci. Remote Sens.; 2022; 60, 4709817. [DOI: https://dx.doi.org/10.1109/TGRS.2022.3226216]

32. Yang, Z.F.; Li, Z.W.; Zhu, J.J.; Hu, J.; Wang, Y.J.; Chen, G.L. InSAR-Based Model Parameter Estimation of Probability Integral Method and Its Application for Predicting Mining-Induced Horizontal and Vertical Displacements. IEEE Trans. Geosci. Remote Sens.; 2016; 54, pp. 4818-4832. [DOI: https://dx.doi.org/10.1109/TGRS.2016.2551779]

33. Chen, B.; Mei, H.; Li, Z.; Wang, Z.; Yu, Y.; Yu, H. Retrieving Three-Dimensional Large Surface Displacements in Coal Mining Areas by Combining SAR Pixel Offset Measurements with an Improved Mining Subsidence Model. Remote Sens.; 2021; 13, 2541. [DOI: https://dx.doi.org/10.3390/rs13132541]

34. Jiang, C.; Wang, L.; Yu, X.X.; Wei, T.; Chi, S.S.; Guo, Q.B. Prediction of 3D deformation due to large gradient mining subsidence based on InSAR and constraints of IDPIM model. Int. J. Remote Sens.; 2021; 42, pp. 188-219. [DOI: https://dx.doi.org/10.1080/01431161.2020.1804088]

35. Chuang, J.; Lei, W.; Yu, X.Y.; Chi, S.S.; Tao, W.; Guo, Z.C. A DPIM-InSAR method for monitoring mining subsidence based on deformation information of the working face after mining has ended. Int. J. Remote Sens.; 2021; 42, pp. 6333-6361. [DOI: https://dx.doi.org/10.1080/01431161.2021.1931540]

36. Yang, Z.; Li, Z.; Zhu, J.; Preusse, A.; Hu, J.; Feng, G.; Wang, Y.; Papst, M. An InSAR-Based Temporal Probability Integral Method and its Application for Predicting Mining-Induced Dynamic Deformations and Assessing Progressive Damage to Surface Buildings. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2018; 11, pp. 472-484. [DOI: https://dx.doi.org/10.1109/JSTARS.2018.2789341]

37. Cui, X.; Deng, K. Research review of predicting theory and method for coal mining subsidence. Coal Sci. Technol.; 2017; 45, pp. 160-169.

38. Cui, X.; Liao, X.; Zhao, L.; Jin, R. Discussion on the time function of time dependent surface movement. J. China Coal Soc.; 1999; 24, pp. 453-456.

39. Li, P.; Wan, H.; Xu, Y.; Yuan, X.; Zhao, Y. Parameter inversion of probability integration method using surface movement vector. Chin. J. Geotech. Eng.; 2018; 40, pp. 767-776.

40. Liu, B.; Dai, H. Research Development and Origin of Probability Integral Method. Coal Min. Technol.; 2016; 21, pp. 1-3.