Prediction Method for Dynamic Subsidence Basin in Mining Area Based on SBAS-InSAR and Time Function

Abstract

Dynamic predictions of surface subsidence are crucial for assessing ground damage and protecting surface buildings. Based on Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-InSAR) technology, a method for making dynamic predictions of large-scale surface subsidence in mining areas can be established; however, the problem of phase coherence loss in InSAR data makes it impossible to predict the complete dynamic subsidence basin. In this study, a method combining the WeiBull time function and the improved probabilistic integral method (IPIM) model was established based on the PIM model, and a method for predicting the dynamic subsidence basin in the mining area was proposed by integrating the IPIM and the combined WeiBull time function. Time-series subsidence data, obtained using SBAS-InSAR, were used as fitting data, and the parameters of the combined WeiBull function were inverted, pixel by pixel, to predict the dynamic subsidence of the working face in the study area. Based on the predicted surface subsidence results of a certain moment in the working face, the parameters of the IPIM model were inverted to predict the subsidence value in the incoherent region. The subsidence predictions of the combined WeiBull time function and the IPIM model were fused using inverse distance weighting (IDW) interpolation to restore the complete subsidence basin in the mining area. This method was tested at the Wannian Mine in Hebei, and the obtained complete subsidence basin was compared with the measured data, with an absolute error range of 0 to 10 mm. The results show that the dynamic subsidence basin prediction method for the SBAS-InSAR mining area, involving the combination of the IPIM model and the combined WeiBull model, can not only accurately fit the time series of surface observation points affected by mining but also accurately restore the subsidence data in the incoherent region to obtain complete subsidence basin information in the mining area.

1. Introduction

Coal resources still account for over 50% of the energy mix in China. However, the high-intensity mining of coal resources can cause the redistribution of overlying stratum stress in goaf areas, leading to the movement and deformation of rock layers and resulting in geological hazards such as cracks and collapses on the surface [1,2]. The subsidence of a land surface caused by the extraction of coal resources is the process of the redistribution of overlying rock stress and long-term dynamic changes on the surface. It is inadequate to predict only the damage due to ground stability following mining operations. The extent of the damage to surface buildings or structures caused by mining subsidence not only depends on the final results of various indicators (static values), but also on the changes in the process of various indicators (dynamic values) [3]. Currently, the damage assessment of surface buildings primarily utilizes indicators such as inclinations, curvature, and horizontal deformation, and assesses damage levels based on the measured data of subsidence and horizontal movement. Therefore, restoring the complete dynamic subsidence basin on the surface is of significant theoretical and practical importance for studying the dynamic changes caused by mining subsidence [4,5].

The surface dynamic subsidence caused by coal mining is a complex surface deformation process. Under conditions of uniform mining, the dynamic subsidence curve of the surface presents a sigmoid shape, and can be divided into three stages: initial subsidence, main subsidence, and residual subsidence [6]. In response to this characteristic, Knothe proposed a time function to describe the dynamic subsidence process of mining activities. However, the shape of its subsidence–time curve is not an “S” shape; therefore, it cannot fully reflect the three stages of surface subsidence. Based on this, some scholars have proposed improved models, such as the power Knothe and segmented Knothe models [7,8,9], as well as growth curve models, such as the WeiBull, Usher, MMF, and logistic models [10,11,12]. The inversion of time-function parameters based on limited Global Navigation Satellite System (GNSS) and leveling measurement data and the prediction of future surface subsidence trends have been widely applied in the field of mining subsidence. However, this method requires a significant amount of manpower and resources due to the limited measured data for predicting the large-scale subsidence caused by mining activities [13].

InSAR technology can be used for deformation monitoring over large areas with short revisit periods, which cannot be achieved by traditional leveling measurements or GNSS measurements [14]. Differential InSAR (D-InSAR)/Time-Series InSAR (TS-InSAR) has been widely used for monitoring surface deformations and has played an important role in monitoring seismic activities, volcano eruptions, landslides, groundwater-induced subsidence, mining-induced deformations, and land deformation nearby airport areas [15,16,17,18,19,20]. In 2002, Berardino proposed the SBAS-InSAR technology; compared to traditional measurement methods, it has the characteristics of high accuracy, all-weather monitoring, and wide coverage in a short period of time [21]. This technology has provided an effective monitoring means for predicting large-scale subsidence in coal-mining processes, and has further developed pixel-offset SBAS and multidimensional SBAS [22]. Currently, there are two approaches for predicting large-scale subsidence caused by mining activities based on SBAS-InSAR technology. The first approach is to invert the parameters of the time-function model based on the monitoring results of SBAS-InSAR technology, and forwardly predict the subsidence or line-of-sight (LOS) deformation caused by mining-induced goaf areas. Zhang et al. identified the conditions necessary for making coal-mining subsidence predictions based on SBAS-InSAR technology and proposed a dynamic prediction method for mining subsidence by combining SBAS-InSAR with a logistic time function [23]. Yang et al. combined the WeiBull model with Kalman filtering to predict the line-of-sight displacement caused by mining activities, and based on prior knowledge, they performed three-dimensional deformation decomposition to obtain information on the three-dimensional surface deformation caused by mining activities [24]. The second approach involves using deep learning or machine learning methods, such as Long Short-Term Memory (LSTM) models or Support Vector Machine (SVM) models, for predicting the dynamic subsidence in mining areas [25,26,27]. However, due to issues such as overfitting and an unclear physical significance, it is difficult to widely apply these methods in engineering. Currently, the focus of research remains on accurately predicting complete surface subsidence by integrating SBAS-InSAR technology with time-function models.

However, the aforementioned research still faces some challenges. Changes in the surface deformation or movement, as well as variations in the backscattering characteristics of soil and vegetation caused by seasonal climate changes, are the main reasons for the temporal and spatial decorrelation in InSAR [28]. In China, mining areas are often located in remote rural areas where the surface is covered with extensive farmland and dense vegetation. The seasonal changes in vegetation reduce the coherence of SAR image interferograms. The decorrelation phenomenon leads to inaccurate monitoring data in D-InSAR/TS-InSAR monitoring results, resulting in the inability to monitor the complete dynamic subsidence basin. On the other hand, current research often uses a single time function to make surface dynamic predictions. However, studies have shown that the subsidence caused by coal mining in thick, loose-stratum mining areas mainly consists of mining-induced subsidence and additional settlements of the soil, resulting from the combined effects of overlying bedrock and loose strata. A single time function may not accurately describe the superimposed effects of the bedrock and soil [29,30].

In order to address the aforementioned issues and obtain accurate and complete dynamic subsidence basins, this study established a combined WeiBull time function and the IPIM model and proposed a method for extracting dynamic subsidence basins in mining areas by integrating the IPIM and the combined WeiBull time function. Using SBAS-InSAR technology, the surface subsidence caused by coal mining was monitored to obtain time-series subsidence data in the mining area. The Particle Swarm Optimization (PSO) algorithm was employed to calculate the parameters of the combined WeiBull model, pixel by pixel, and the combined WeiBull model was used to dynamically predict surface subsidence. The dynamic prediction results were used to invert the parameters of the IPIM model and forwardly predict the subsidence values in the incoherent region. Based on the IDW algorithm, the prediction results from the IPIM and the combined WeiBull model were fused to restore the complete dynamic subsidence basin in the mining area, which was then compared with actual leveling data for validation.

2. Study Area and Data Description

The study area was located in the northern part of the Wannian Mine, with the Nanming River crossing the entire area, which is usually dry. The terrain is mostly flat, primarily consisting of farmland and some villages, as shown in Figure 1. The space for coal mining, where surface subsidence resulted from the excavation of the coal seam during underground coal mining, is referred to as the coal mining working face. The 132,157 working face is a coal mining face located in the Wannian Mine. Taking the 132,157 working face as an engineering example, the area above the working face is flat, with the villages and vegetation distributed. The mining conditions of the working face are described in

Table 1. Mining conditions of the working face in this research.

Length in Strike Direction (m)	Length in Dip Direction (m)	Coal Seam Dip (°)	Average Mining Depth (m)	Average Coal-Mining Thickness (m)
512	137	27	900	3.3

. The extraction of coal from the working face began in January 2018, with a mining direction depth ratio (n₃ = D₃/H) of 0.7 < (1.2~1.4), and a mining inclination depth ratio (n₁ = D₁/H) of 0.57 < (1.2~1.4), both falling under the category of insufficient mining. The coal seam being extracted is the Permian 2# coal seam, and the full-caving roof management method is being employed.

This study used the European Space Agency’s Sentinel-1A data with VV polarization and the IW working mode in single look complex (SLC) format. The data had a spatial resolution range of 5 m, an azimuth spatial resolution of 20 m, and a revisit period of 12 days. The selected data covered the period from 8 January 2018 to 17 September 2018, including 21 scenes of Sentinel-1A data covering the 132,157 working face of the Wannian Mine, Digital Elevation Model (DEM) data with a spatial resolution of 30 m, and precise orbit data corresponding to the Sentinel-1A data. Detailed information on the SAR data is shown in

Table 2. Acquisition dates of Sentinel-1A SAR data.

No.	Time	No.	Time	No.	Time
1	8 January 2018	8	2 April 2018	15	7 July 2018
2	20 January 2018	9	14 April 2018	16	19 July 2018
3	1 February 2018	10	26 April 2018	17	31 July 2018
4	13 February 2018	11	8 May 2018	18	12 August 2018
5	25 February 2018	12	1 June 2018	19	24 August 2018
6	9 March 2018	13	13 June 2018	20	5 September 2018
7	21 March 2018	14	25 June 2018	21	17 September 2018

. The distribution of measured leveling points is shown in Figure 1c. The leveling measurements were conducted on 11 January 2018, 7 April 2018, 19 May 2018, 26 June 2018, and 28 September 2018; thus, subsidence data for four consecutive time periods was obtained.

Based on the SBAS-InSAR monitoring results from 8 January 2018 to 5 August 2018, the parameters of the combined WeiBull model were calculated, pixel by pixel, and used to predict the future subsidence amounts. The SBAS-InSAR monitoring results from 5 September 2018 and 17 September 2018, covering two time periods, were selected to validate the accuracy of the combined WeiBull model predictions. Based on the dynamic prediction results from 17 September 2018, the IPIM parameters were inverted and used to forwardly predict the surface subsidence basin. Finally, by combining the dynamic prediction results with the IPIM prediction results, the complete dynamic subsidence basin in the study area was restored.

3. Methods

The proposed dynamic subsidence basin extraction model is mainly composed of three parts. Section 3.1 describes the construction of the combined WeiBull time-function model. Section 3.2 describes the construction of the IPIM model. Section 3.3 introduces the workflow of the dynamic subsidence extraction method based on SBAS-InSAR.

3.1. The Combined WeiBull Time-Function Model

The WeiBull model is a continuous probability distribution model that was proposed by Waloddi Weibull in 1951. This model exhibits central symmetry and includes three characteristic inflection points (the point of maximum acceleration, the point of maximum velocity, and the point of minimum acceleration), which describe the dynamic subsidence process induced by mining activities [31]. It can be expressed as follows:

$$W(t)=W_m(1-e^{-a{t}^b})$$

(1)

where a and b are the model parameters related to the properties of the overlying strata above the mining face and W_m represents the maximum subsidence at that point.

In deep-buried, thick, loose-stratum mining areas, the subsidence caused by coal mining is mainly induced by mining activities and the additional settlement of the soil, which is the result of the combined effects of the overlying bedrock and alluvial layers. Wang represented the subsidence basin caused by mining activities by combining units of subsidence basins with different main influence radius in a certain proportion. Following the superposition principle, they established a model for predicting the subsidence and deformation at any point within the surface movement basin. Additionally, based on the data from multiple observation stations, they provided empirical formulas for the proportion coefficients in relation to the thick, loose strata and the average depth [32], as shown in Equation (2):

$$\rho =\left\{\begin{array}{c}0.0053\exp (\frac{25.344h}{H_0-h}),\frac{h}{H_0-h}<0.2\\ 0.0796\frac{h}{H_0-h}+0.201,\frac{h}{H_0-h}\ge 0.2\end{array}\right.$$

(2)

The proportion coefficient ρ is mainly related to the thickness of the loose layer h and the average mining depth H₀. For a certain mining depth, the larger the thickness of the loose layer (that is, the greater the ratio of the thickness of the loose layer to the bedrock H₀ − h), the larger the proportion coefficient p. When the ratio of the thickness of the loose layer to the bedrock is less than 0.2, the proportion coefficient ρ exponentially increases with an increase in the ratio. However, when the ratio is greater than 0.2, ρ increases approximately linearly with the ratio [32].

Therefore, based on the proportion relationship given by Equation (2) above, this study constructed a combined WeiBull model, which involved combining two WeiBull time functions with different parameters in a certain proportion to describe the dynamic surface subsidence process. This model, known as the combined WeiBull model, is theoretically applicable not only for dynamic model predictions under the conditions of thick loose-stratum mining, but also for general mining conditions to establish a dynamic prediction model for surface subsidence, as shown in Formula (3):

$$w(t)=\left\{\begin{array}{c}{w}_1(t)(1-0.0053\exp (\frac{25.344h}{H_0-h}))+w_2(t)0.0053\exp (\frac{25.344h}{H_0-h}),\frac{h}{H_0-h}<0.2\\ w_1(t)(1-0.0796\frac{h}{H_0-h}+0.201)+w_2(t)0.0796\frac{h}{H_0-h}+0.201,\frac{h}{H_0-h}\ge 0.2\end{array}\right.$$

(3)

3.2. Building the IPIM Model

The probability integral method is a widely used subsidence prediction method for coal mining in China. The calculation formula for subsidence in the prediction model of the probability integral method is as follows [33]:

$$w(x)=\frac{m_{c}q\cos \theta }{2}(erf(\frac{\sqrt{\pi }}{r}x)+1)$$

(4)

In the above formula, m_c represents the coal seam thickness; x is the mining coordinate system’s strike coordinate; w(x) is the predicted subsidence value in the main section; θ denotes the coal seam dip angle; r is the influencing radius of mining; and q represents the subsidence coefficient.

The prediction formula for surface subsidence within the strike and dip of the main section is as follows, the geometric relationship is shown in Figure 2

$$\left\{\begin{array}{c}{w}^0(x)=w(x)-w(x-l_3)\\ w^0(y)=w(y)-w(y-l_1)\\ l_1=(D_1-S_1-S_2)\frac{\sin ({\theta }^0+\theta )}{\sin {\theta }^0}\\ l_3=D_3-2S\end{array}\right.$$

(5)

In the above formula, D₁ represents the dip mining length; D₃ is the strike mining length; θ⁰ denotes the mining influence propagation angle; y is the dip coordinate in the mining coordinate system; w⁰(x) is the predicted subsidence value in the strike of the main section; w⁰(y) is the predicted subsidence value in the dip of the main section; l₁ stands for the dip face’s computational length; l₃ is the strike face’s computational length; S represents the deflection distance of the strike inflection point; S₁ is the deflection distance of the inflection point in the downhill direction; and S₂ is the deflection distance of the inflection point in the uphill direction. The prediction formula for the subsidence and horizontal movement of any point on the surface is as follows:

$$D_v=\frac{1}{w_0}{w}^0(x)w^0(y)$$

(6)

In the above formula, w₀ represents the maximum subsidence in the main section and D_v is the subsidence value at any point on the surface.

Currently, the probability integral method remains a widely used theoretical method for subsidence predictions in China. However, in recent years, most coal mines have transitioned to deep mining, making inadequate mining conditions more common. The probability integral method is suitable for adequate mining conditions, but not appropriate for predicting subsidence under inadequate mining conditions [34]. According to the parameter modification method for predicting subsidence under inadequate mining conditions established by Wu [35], which involves introducing a mining intensity factor to modify the prediction parameters, this study established a formula for predicting the surface subsidence under inadequate mining conditions.

$$\frac{q_d}{q}=\left\{\begin{array}{cc}0.97n^2-0.07n+0.39& 0.1<n\le 0.83\\ 1.0& n>0.83\end{array}\right.$$

(7)

In the above formula, n is the mining intensity factor; q_d is the subsidence coefficient under inadequate mining conditions; and q is the subsidence coefficient under adequate mining conditions. According to Equation (7), the probability integral method model in Equation (4) can be modified as follows:

$$w(x)=\frac{m(0.97n^2-0.07n+0.39)\mathrm{q}\cos \theta }{2}(erf(\frac{\sqrt{\pi }}{r}x)+1)$$

(8)

After modifying the probability integral method, the mining intensity factor becomes one of the prediction parameters in the probability integral method. The expressions for the corrected geological parameter GP and the prediction parameter P are as follows:

$$\left\{\begin{array}{c}GP{=[\mathrm{H},\mathrm{H}}_1{,\mathrm{H}}_2,D_1,D_2,\theta ,m_c]\\ P=[q,\tan \beta ,{\theta }^0,S,S_1,S_{2,},n]\end{array}\right.$$

(9)

In the above formula, H is the mining depth in the strike direction; H₁ is the downslope mining depth; H₂ is the upslope mining depth; D₁ is the dip mining length; D₃ is the strike mining length; θ is the dip angle of the coal seam; m_c is the thickness of the coal seam; q is the subsidence coefficient; tan β is the tangent of the main influence angle; θ⁰ is the mining influence propagation angle; n is the mining intensity factor; S is the deflection distance of the strike inflection point; S₁ is the deflection distance of the downslope inflection point; and S₂ is the deflection distance of the upslope inflection point.

Based on Equations (4)–(8), the probabilistic integral method provides a modified relationship between the subsidence amount D_v and the prediction parameters GP and P.

$$D_v=f(GP,P)$$

(10)

This paper refers to the modified probabilistic integral method as the IPIM (improved probabilistic integral method). The IPIM model is applicable for subsidence prediction under both inadequate and adequate mining conditions.

3.3. The SBAS-InSAR Mining Area Dynamic Subsidence Extraction Method, Combining the IPIM and the Combined WeiBull Model

Based on the established IPIM and the combined WeiBull time-function model, a method for extracting the dynamic subsidence in mining areas was combined with SBAS-InSAR. The dynamic subsidence W₁(t_n) (n = 1, 2, 3…), predicted by the combined WeiBull model at time t_n, was used as the fitting data for the IPIM. The PSO algorithm was used to invert the IPIM parameters, and a forward prediction of the subsidence basin was carried out.

$$\left\{\begin{array}{c}{\mathrm{W}}_1(t_n)\stackrel{\mathrm{PSO}}{\to }{t}_n(GP,P)(n=1,2,3\dots )\\ t_n(GP,P)\stackrel{\mathrm{IPIM}}{\to }{\mathrm{W}}_2(t_n)\end{array}\right.$$

(11)

In the above equation, t_n (GP, P) represents the inverted IPIM model parameters and W₂(t_n) represents the dynamic subsidence basin predicted by the IPIM.

The inverse distance weighting interpolation method assumes that each input point has a local influence that weakens with increasing distance, which is consistent with the distribution of contour lines for surface subsidence [36]. The inverse distance weighting interpolation method can be used to integrate the predicted results of the IPIM and the combined WeiBull model, as shown in Equation (12).

$$\mathrm{P}(x_i,y_i)={\displaystyle \sum _{i=1}^n{\lambda }_i{P}_i(x_i,y_j)}$$

(12)

In the above equation, P(x_i, y_i) represents the area of pixels for which the solution needs to be determined.

The formula for calculating the weight λ_i of each point is as follows:

$${\lambda }_{\mathrm{i}}=\frac{\frac{1}{d_i}}{{\displaystyle \sum _{i=1}^n\frac{1}{d_i}}}$$

(13)

where d_i represents the absolute distance from each effective pixel to the pixel being solved.

By integrating the IPIM’s predicted result W₂(t_n) and the combined WeiBull model’s predicted result W₁(t_n) through IDW interpolation fusion, the complete dynamic subsidence basin W(t_n) of the mining area at time t_n can be restored. This method’s process is shown in Figure 3.

$$W({\mathrm{t}}_n)=IDW({\mathrm{W}}_1(t_n),{\mathrm{W}}_2(t_n))$$

(14)

Based on the models constructed in Section 3.1, Section 3.2 and Section 3.3, a flowchart depicting the dynamic subsidence basin extraction process within the mining area can be established as shown in Figure 3. The specific procedures are as follows:

• SBAS-InSAR processing: The SBAS-InSAR technology was used to acquire the cumulative LOS deformation of the mined-out area, and the LOS deformation was converted to vertical deformation based on the radar incidence angle.
• Forward subsidence prediction: Parameters P were calculated based on the geological conditions of the working face, following which a combined WeiBull model was constructed. The Singular Spectrum Analysis (SSA) algorithm was utilized to eliminate error data, and subsequently, the parameters of the combined WeiBull model were calculated on a pixel-by-pixel basis to predict forward subsidence.
• Integrating InSAR and IPIM data: Based on the forward prediction results of the combined WeiBull model in step 2, the subsidence values w(x,y) of the predicted basin strike and dip were extracted as global fitting data. The objective function f of the PSO algorithm was constructed as shown in Equation (15). The predicted parameters P were inverted based on the PSO algorithm, and the subsidence values were subsequently predicted using the IPIM model. Based on the IDW algorithm, the prediction results from the IPIM and the combined WeiBull model were fused to restore the complete dynamic subsidence basin in the mining area.

$$f=\min \left|\right|w(x,y)-w(x,y,GP,P)\left|\right|$$

(15)

4. Results

The Sentinel-1A data were processed using the interferometric overlap data-processing module in the ENVI SARscape 5.2.1 software. The software was configured with a spatial–temporal baseline threshold to generate connection pairs for maintaining the coherence of interferometric pairs and creating as many as possible. The date selected for the master image was 2 April 2018, with the spatial–temporal baseline thresholds set at 45 days and 100 m. The spatial–temporal baseline connectivity map is shown in Figure 4. The filtering method was set to Goldstein’s filtering method, and the phase-unwrapping method was the minimum cost flow method.

Firstly, after interferogram formation, phase unwrapping, filtering, a coherence calculation, and unwrapping, a total of 66 interferometric pairs were generated. After inspection, five pairs with a low coherence were removed. Secondly, in regions without residual topographic fringes, deformation fringes, or phase discontinuities, 20 ground control points were appropriately selected for orbit refinement and re-interferogram formation. Next, SBAS-InSAR was performed with two inversion steps: in the first inversion, the displacement rates and residual topography were estimated, and a second phase unwrapping was conducted to optimize the results. The second inversion involved computing the displacements in the time series. Based on the deformation rates obtained, a customized atmospheric filter was applied to estimate and remove atmospheric phase contributions, yielding a purer time-series displacement [21].

InSAR is most sensitive to vertical displacement, followed by east–west displacement, and it is least sensitive to north–south displacement [37]. Therefore, the line-of-sight deformation obtained from SBAS-InSAR monitoring was converted to vertical deformation. The calculation formula for the subsidence value w(x, y) at any point on the ground is as follows:

$$w(x,y)=\frac{w_{los}(x,y)}{\cos \theta }$$

(16)

In the above equation, w_los(x,y) represents the deformation value in the line-of-sight direction at any point and θ represents the radar incidence angle.

The SBAS-InSAR monitoring period was from 8 January 2018 to 17 September 2018. The results for six time points were selected, as shown in Figure 5. From 8 January 2018 to 17 September 2018, the SBAS-InSAR monitoring results indicated a maximum surface subsidence of 164 mm. Using empirical parameters for the mining area, the forward prediction of the static subsidence basin caused by mining at this working face indicated a maximum subsidence of approximately 169 mm. Due to a coal seam dip angle of 27° at this working face, the subsidence basin caused by mining was inclined in the downslope direction. The maximum subsidence position observed through InSAR monitoring was in good agreement with the predicted maximum subsidence position according to the IPIM, as shown in Figure 5f.

In order to quantitatively verify the InSAR monitoring results, they were compared with single-point, multi-time-period, measured level data. Since the InSAR monitoring date was not exactly the same as the GNSS monitoring date, the SBAS-InSAR monitoring values were linearly interpolated to obtain the settlement values corresponding to the InSAR monitoring on the measured date. The linear interpolation calculation method used is shown in Equation (17):

$$\mathrm{InSAR}(t_m)=\frac{\mathrm{InSAR}(t_2)-\mathrm{InSAR}(t_1)}{t_2-t_1}(t_m-t_1)+\mathrm{InSAR}(t_1)$$

(17)

where InSAR(t_m) represents the measured settlement value corresponding to the measured date; t₁ and t₂ represent the nearest SAR satellite observation times before and after the measured date; and InSAR(t₁) and InSAR(t₂) are the InSAR monitoring values corresponding to the SAR satellite observation time before and after the measured date.

The comparison of level data with the InSAR monitoring results is shown in Figure 6, from which it can be seen that the trends in the InSAR results and the measured data were consistent. Due to the fact that the times used to collect the InSAR monitoring data and the measured level data did not exactly correspond to each other, linear interpolation was used, and the method used to calculate the InSAR monitoring values that corresponded to the time of the actual measurements is shown in Equation (17). The results of the comparison are shown in

Table 3. Comparison of measured data with InSAR monitoring values.

Point Name	Field Measurement Time	Field Measurement Value (mm)	InSAR Measurement Value (mm)	Difference (mm)
q7	7 April 2018	−53.0	−51.8	−1.2
	19 May 2018	−98.9	−98.6	−0.3
	26 June 2018	−102.9	−106.8	3.9
	28 September 2018	−112.7	−114.9	2.2
q8	7 April 2018	−59.2	−57.9	−1.3
	19 May 2018	−107.7	−103.2	−4.5
	26 June 2018	−112.3	−107.0	−5.3
	28 September 2018	−123.8	−118.0	−5.8
q10	7 April 2018	−68.2	−75.9	7.7
	19 May 2018	−118.4	−125.0	6.6
	26 June 2018	−123.6	−133.9	10.3
	28 September 2018	−137.2	−141.5	4.3
q12	7 April 2018	−75.6	−85.4	9.8
	19 May 2018	−123.8	−130.1	6.3
	26 June 2018	−129.2	−138.5	9.3
	28 September 2018	−143.9	−145.8	1.9

. The absolute error distribution between the InSAR and the actual measurement results ranged from 0.3 to 10.3 mm, with an average absolute error of 5 mm. An analysis of the difference between the InSAR-monitored settlements and the measured level settlements showed that these two methods had good consistency.

4.1. Combining the WeiBull Projected Results with the IPIM Projected Results

To obtain the thickness of the loose layer and the thickness of the bedrock, Formula (2) can be used to calculate a scale factor of 0.22 for the combined WeiBull model, according to the geological data of the mining area; then, the p-value of the combined WeiBull method can be calculated. In this study, the time-series settlement data from 8 January to 24 August 2018 were used to calculate the parameters of the combined WeiBull model pixel by pixel, with a final settlement greater than 50 mm, using the PSO global optimization algorithm. Based on the combined WeiBull time function obtained by fitting the time-series settlement data from 8 January to 24 August 2018, the cumulative surface settlement of the mine site on 5 September 2018 and 17 September 2018 was predicted, and the expected results are shown in Figure 7.

Based on the combined WeiBull model, the subsidence values of the strike and inclination of the working face predicted on 17 September 2018 were extracted as the fitted values of the IPIM model. The inversion of the IPIM parameters yielded the following results: subsidence coefficient q = 0.78; mining coefficient n = 0.45; mining influence angle tangent tanb = 1.6; and a mining influence propagation angle of 80°. The positively predicted settlement basin of the IPIM is shown in Figure 7b.

4.2. Integration of the InSAR and IPIM Projected Results

Based on the IDW interpolation fusion of the joint projected IPIM results calculated using Equation (14) and the projected results of the combined WeiBull model on 17 September, the complete settlement basin of the mine surface was recovered. Most of the buildings in the village were affected by the mining of the 132,157 working face, as shown in Figure 8a. In order to quantify the damage to the buildings in the village, surface movement deformation, curvature, and tilt values can be determined based on the above IPIM parameter inversion results and the empirical parameters of the horizontal movement coefficient in this mining area, as shown in Figure 8b–g. According to the assessment criteria for the class I of buildings, the maximum values for inclinations, curvature, and horizontal deformation are 3 mm/m, 0.2 mm/m², and 0.6 mm/m, respectively. This working face was subjected to deep and very insufficient mining, the surface movement deformation was small, the settlement was slow, the maximum tilt of the surface was 0.48 mm/m, the maximum curvature was 0.003 mm/m², the maximum horizontal deformation was 0.6 mm/m, and the degree of damage to the surface village was categorized as class I.

In order to test the accuracy of the dynamic subsidence basin extraction method described in this study, the expected results were compared with the measured results of the corresponding real points. The results showed that the range of absolute error was 0~10 mm and the root–mean–square error (RMSE) was 4 mm, as shown in Figure 9.

5. Discussion

Compositional WeiBull Functions vs. Other Time Functions

Comparison experiments were conducted using logistic, WeiBull, Usher, MMF, and power Knothe models. The predicted results of these six time functions were validated separately based on the SBAS-InSAR monitoring results on 5 September. The 5 September InSAR monitoring results and the surface settlement values predicted based on the above six time functions were plotted on the x-axis and y-axis, respectively, to create a consistency validation plot, as shown in Figure 10. The root–mean–square errors calculated using each model are shown in

Table 4. Expected accuracy of model results according to RMSE (mm).

	This Study’s Method	WeiBull	Usher	Power Knothe	MMF	Logistic
5 September 2018	2.7	3.4	29.0	3.1	3.5	7.7
17 September 2018	3.1	3.8	29.0	3.6	4.2	8.8

. Figure 10 shows that most of the points from the combined WeiBull method were distributed close to a straight line along y = x, indicating that the predicted settlement values were consistent with the measured settlement values. Compared with other time-function models such as the logistic and WeiBull models, the combined WeiBull time-function model proposed in this study can accurately describe the surface settlement process caused by quarrying, and the predicted accuracies of the combined WeiBull model were all better than those of other models.

As can be seen from Table 4, with an increase in the forward prediction time interval, the predicted accuracies all slightly decreased, but the combined WeiBull time-function accuracy was still optimal. The time-series monitoring results of two SBAS-InSAR monitoring sites were extracted and the fitting effects of multiple time-function models are shown in Figure 11; these fitting effects indicate that the combined WeiBull model can better describe the dynamic process of surface settlement than other models.

Although the research method in this paper has achieved excellent experimental results compared to other time functions, more information needs to be considered under certain geological conditions due to the complexity of the mining area. In the engineering case discussed in this paper, the mining working face involved deep and insufficient mining. The InSAR/TS-InSAR technology obtained relatively complete monitoring information. However, under certain shallow-buried and high-intensity mining conditions, X-band SAR satellites may not have sufficient monitoring capabilities to obtain complete deformation information. In this scenario, the use of multi-source data (drones, LiDAR) is necessary to enhance the research method presented in this paper.

The PIM/IPIM model can be used for subsidence prediction in most mining working faces, but complex geological conditions such as faults within the strata may exist [38]. In such cases, the PIM/IPIM model no longer conforms to the subsidence laws. Other models will be utilized to replace the PIM/IPIM model and restore the complete subsidence in the mining area.

6. Conclusions

With the aim of addressing the shortcomings of SBAS-InSAR technology in the dynamic prediction of large-scale subsidence in mining areas, a combined WeiBull model with the IPIM was established based on PIM and WeiBull time functions. As a result, an SBAS-InSAR method for making dynamic subsidence basin predictions in mining areas by integrating the IPIM and the combined WeiBull model was proposed. The study area was located in the Wannian mining area, and the method was able to recover the dynamic subsidence of the lost coherence area in the study area and determine the complete dynamic subsidence basin of this working face.

The time-function model based on the combined WeiBull model was more suitable for describing the dynamic subsidence of the surface caused by mining. Comparison experiments with logistic, WeiBull, Usher, MMF, and power Knothe models were conducted, and the predicted accuracies of the combined WeiBull model were all better than those of the other models. This indicates that the combined WeiBull model can fit the settlement time series of the surface observation points in the mining subsidence area with a high accuracy, and it is fully consistent with the settlement process law of the observation points in the surface subsidence area caused by underground mining.

Due to the complexity of the surface conditions in the mining area, the subsidence basins obtained by InSAR were usually “out of coherence”, which made it difficult to obtain complete subsidence basins. The subsidence extraction method that included integrating the IPIM and the combined WeiBull model was able to make accurate dynamic prediction results and recover the subsidence data in the out-of-coherence areas. The subsidence basins obtained using the method proposed in this study had an absolute error distribution of 0~10 mm when compared with the measured level data, which made up for the shortcomings of the SBAS-InSAR technique alone in predicting the dynamic subsidence.

In an assessment of the damage to surface village buildings using the model in this study, the maximum tilt of the surface was 0.48 mm/m, the maximum curvature was 0.003 mm/m², and the maximum horizontal deformation was 0.6 mm/m. According to the coal-mining regulations, the degree of damage to the surface village in this working face was categorized as class I. With the increasing demand for higher accuracy in predicting the surface dynamic deformation of large-scale infrastructure in current coal-based cities, the study method proposed will have better practical application.

Although the complete dynamic subsidence basin of the mining area can be extracted using the methods in this study, mining-induced surface movement involves a three-dimensional spatial change, including subsidence and horizontal movement in the east–west and north–south directions. The insensitivity of InSAR technology to deformation in the north–south direction makes it difficult to accurately acquire mining-induced 3D deformation with single-view, line-of-sight InSAR. Further in-depth research is needed to combine multi-source data to assess three-dimensional dynamic changes in the ground surface. On the other hand, the surface deformation caused by mining is the movement of spatial vectors, and the LOS deformation monitored by InSAR is essentially the movement of spatial vectors as well. The next step will be to further investigate the relationship between them.