IPIM-G Dynamic Prediction Model of Mining Subsidence Based on D-InSAR Technical Parameter Inversion

: The probability integral method (PIM) model is more suitable for mining subsidence prediction of underground working-face of coal mine in China and has been widely used. However, the PIM model has the question on too fast convergence in predicting at the edge of subsidence basin. In recent years, many scholars have studied a lot of subsidence monitoring methods in the coal mine area by using the technical advantages of differential interferometry synthetic aperture radar (D-InSAR). But, serious incoherence of interferometry phase occurs because of the large gradient subsidence of mining area, which leads to the inability to accurately obtain large gradient subsidence of surface. Meanwhile, the PIM model is more suitable for static prediction of mining activity in the mining area, yet has certain defects in dynamic subsidence prediction during mining process. In view of the above shortcomings, the improved PIM (IPIM) prediction model was first introduced through improving the PIM model in the paper, and the IPIM-G dynamic prediction model was constructed based on the PIM model and the Gompertz time function for mine-area mining. Then, a method was used to invert the parameters of the IPIM-G dynamic prediction model by using time series superposition results of surface subsidence monitored by the D-InSAR technology, and then the parameters obtained by inversion were applied to the IPIM-G model for mining subsidence prediction. In order to verify the effectiveness of the IPIM-G model, the model was applied to a coal mine in Huaibei mining area, Anhui Province. The experimental results show that the predicted values of the IPIM-G dynamic prediction model are in good agreement with the measured leveling values in field, and the accuracy of the IPIM-G dynamic prediction model is 85.9% higher than the monitoring results of the D-InSAR technology in obtaining the gradient subsidence information of the subsidence basin in the mining area. The accuracy of subsidence prediction at the edge of the subsidence basin is 14.9% higher than the PIM model.


Introduction
Underground mining will lead to large gradient subsidence of the surface and large regional subsidence basin, which will damage the buildings (structures) around the mining area. Therefore, it is particularly important to monitor and predict the large gradient deformation caused by mining and obtain the deformation information at the edge of the subsidence basin. The traditional surface deformation monitoring methods include leveling, satellite navigation and positioning (Li et al., 2015), but these monitoring methods can only monitor local points and lines, and can not realize large-area subsidence observation in the region. Moreover, the observation is easy to be affected by weather, light or terrain, and the observation cost is high. Differential interferometric synthetic aperture radar (D-InSAR) measurement is a new type of earth observation technology, which is widely used for surface subsidence observation with the advantages of all-weather, all-day, wide range and high accuracy (Zhu et al., 2019). Because D-InSAR technology adopts differential interferometry, it has high requirements for the coherence of main and auxiliary images, but the large gradient deformation in the mining area often exceeds the monitoring capacity of D-InSAR (Zhang et al., 2015), resulting in incoherence. It is difficult to obtain the actual deformation of radar line of sight (LOS) after mining subsidence of the working face.
At present, many scholars have conducted some researches for D-InSAR to acquire large gradient deformation of mining area, which can be divided into two main aspects (Chen, 2019): (1)Combining offset-tracking method with D-InSAR technology; (2)Combining mining subsidence prediction model with D-InSAR technology. The first method uses the D-InSAR technique to obtain the cumulative values of time-series deformation in the area with high coherence, while the subsidence in the central area of the subsidence basin is too large to be monitored by the D-InSAR technique, so the offset-tracking method is used to extract the large gradient deformation values Ou et al., 2018;Wang et al., 2020); This method can compensate for the accuracy problems caused by the D-InSAR technique due to phase incoherence in the region of large gradient deformation, but the method is susceptible to the influence of the imagequality itself and the alignment method (Yao, 2016;Zhang et al., 2021), such as changes in surface vegetation cover or artificial disturbance of surface morphology, which will lead to large errors in the final monitoring results. The second method is to dynamically predict the large gradient deformation in the mining area based on the probability integral method (PIM) prediction model or the integration of time function models such as Knothe (Wang et al., 2018;Zhang, 2017), Logistic (Lee et al., 2006;Yang et al., 2017) or Weibull ; In this method, the edge settlement value of large gradient subsidence basin obtained by D-InSAR technology is used to inverse the predicted model parameters, and then the settlement value of large gradient subsidence basin is dynamically predicted, which effectively avoids the loss of coherence in the process of obtaining large gradient deformation information of mining area by D-InSAR technology. A large number of measured studies show that, on the one hand, the above time function model has certain defects. Knothe function is only established when the hypothetical conditions are met , the logistic function model is different from the actual subsidence process, and the physical significance of Weibull function parameters is not clear (Liu, 2010); On the other hand, in the prediction process combined with the time function, PIM has high accuracy in predicting the deformation at the center of surface subsidence, but there is a problem of too fast convergence at the edge of the subsidence basin, resulting in poor accuracy of PIM at the edge of the subsidence basin. Through research, analysis and data access, it is found that this problem is due to the fact that the PIM prediction model is based on certain assumptions, that is, the rock layer above the subsidence basin is regarded as an isotropic, homogeneous and discontinuous medium (Fan et al., 2014), resulting in a deviation between the final calculation process and the actual results.
In this paper, Gompertz function was selected as the time function model and integrated with the improved probability integral method (IPIM) prediction model to construct the IPIM-G dynamic prediction model of large gradient deformation in mining area. Then, D-InSAR technology was used to obtain the surface subsidence in a certain period of time, and the residual function was constructed with the predicted values of IPIM-G, and the optimal parameters were obtained by combining Genetic Algorithm (GA) inversion, which was substituted into the IPIM-G model to predict the forward mining subsidence dynamics. The mining influence range of working face in a coal mining area of Huaibei mining area was selected as the research object. The dynamic prediction results of the IPIM-G model of mining subsidence in the study area were compared with the actual measured values of the mining level in the mining area to verify the applicability and feasibility of the IPIM-G dynamic prediction model.

D-InSAR principle
D-InSAR technology performs differential interference on two SAR images and uses external DEM to remove terrain phase to obtain surface deformation phase information def  . The specific components are as follows: Where ref  is the reference ellipsoid phase, def  is the deformation phase, atm  is the terrain phase, noi  is the atmosphere phase and is the noise phase. The SAR image is processed by differential interferometry, then filtered and coherence calculated, and finally, the deformation phase of the surface is obtained by phase unwrapping and geocoding, and then the Los deformation los W of the surface can be obtained as follows: Since the line of sight of SAR satellite is not perpendicular to the surface, the vertical subsidence of the surface can be obtained according to the trigonometric function relationship: Where  is the incident angle of SAR satellite.

Improved probability integral method (IPIM) prediction model
Probability integral method (PIM) prediction model is based on the random medium theory to predict the surface subsidence. It has been widely used in the subsidence prediction of underground coal mining in Chinese coal mining areas. It regards the rock stratum as a discontinuous ideal loose random medium, and the surface movement and deformation caused by underground mining as a random process. It is concluded that the final surface deformation is the superposition of deformation caused by multiple small unit mining.
The PIM prediction model obtains the final surface subsidence by integrating the unit subsidence influence function within the influence range, and since mining subsidence is a three-dimensional problem, the problem is first transformed to a two-dimensional main section for analysis. The expression of surface unit subsidence basin in the main section is (Deng et al., 2014): Where r is the main influence radius and x is the abscissa of a point on the surface of the main section.
Since the PIM prediction model is based on the assumption that the rock layers are isotropic and homogeneous discontinuous media, the differences between the overlying rock layers are not considered and the influence of loose layers is neglected, thus causing the probability integral method to converge too quickly at the edges of the sinking basin. According to the analysis of Eq. (4), the variability of the overlying rock formations mainly affects the parameters r in the PIM prediction model. Therefore, considering the property differences of overlying bedrock and loose layer, the estimated parameter r in the PIM prediction model is adjusted to control the convergence of the PIM prediction model at the sinking edge. The expression of small unit subsidence value w e (x) in the improved probability integral method (IPIM) model can be (Wang, 2014): Where H is the average mining depth, h is the thickness of loose layer and  is the influence coefficient of loose layer. In Eq. (5), the thickness h of the loose layer can affect the convergence speed at the subsidence boundary, and the influence coefficient of the loose layer  determines that the influence of the loose layer must be considered in the calculation of subsidence, When  infinitely approaches 0, it indicates that the influence of loose layer is very small, and Eq. (5) is basically consistent with Eq. (4). Based on experience and research (Wang, 2014), ( / ) Hh  is often taken to be 0.2r ~ 0.4r. Fig. 1 shows the surface subsidence state of the main section under the semi-infinite mining state, m is the thickness of the mining coal seam, H0 is the mining depth, the abscissa on the main section of unit mining is s, and if the coordinate of any point on the ground on the main section is x, the subsidence value caused by unit mining is w e (x-s), that is, take (x-s) into Eq. (5), then: Assuming that the mining boundary is from 0 to  , the subsidence caused by unit mining is integrated, and the expression of subsidence w (x) on the main section of semi-infinite mining subsidence basin is obtained as follows: Where erf is the error function, 0 W is the maximum surface subsidence value, and the calculation formula is 0 cos W mq   ,  is the coal seam inclination, and q is the subsidence coefficient.

Fig. 1 Surface subsidence of the main section of semi-infinite mining
The subsidence of semi-infinite mining in one direction of the main section is analyzed above, but the mining of the working face is often finite. During finite mining, the subsidence profile function of the main section can be regarded as the superposition of two semi-infinite mining profile functions, so the subsidence of any point A(x, y) along the strike main section and the dip main section during finite mining are as follows: Where l is the calculated length in the x direction of the mining area, and L is the calculated length in the y direction of the mining area.

Gompertz function model
Gompertz function model has high accuracy in fitting "S" type data (Tjorve and Tjorve, 2017). It was first used to study the evolution of an animal population. Through continuous research by scholars, it was found that this function model can well express the change process of surface subsidence with time (Wu and Hu, 2006;Xu et al., 2019), so it can be used to study the law of surface subsidence. The Gompertz function line has an asymmetric "S" shape, and its functional expression is: Where W f is the upper asymptote of growth, k is the growth rate coefficient, and c is the inflection point of the function. It can be seen from Fig. 2 that in the Gompertz time function model, the function curve is "S" type. The subsidence speed of the function before the inflection point gradually increases with time, and the subsidence speed after the inflection point gradually decreases to 0 with time, indicating that the final surface subsidence value tends to be a fixed value. Therefore, the change of the function conforms to the law of surface mining subsidence. When the working face is pushed to a certain surface point, it is regarded as the starting time of subsidence at this point. In this way, the time, subsidence and the advancing position of the working face can be changed synchronously, which is convenient for application in practical engineering and subsequent research.

IPIM-G dynamic prediction model integrating IPIM and Gompertz
Combined with the analysis idea of mining subsidence law and differential theory, as shown in Fig. 3, the working face is divided into n small units, the coal seam inclination is  , the W f in Gompertz time function is regarded as the final subsidence caused by each small unit, and then the IPIM-G dynamic prediction model for subsidence at any point on the surface is established by associating the relationship between IPIM prediction model and Gompertz function.
, then the working face is mined along the x-axis direction, and the subsidence of each small unit at time t is: Where t i is the mining time of the small unit i, x is the abscissa of the mining position, Assuming that the working face advances at a uniform speed, since the length of each small unit divided is equal, the mining time of each small unit is equal, which is recorded as t  , then Eq. (11) becomes: n small units whose subsidence is: At the same time, the calculation formula of y-axis subsidence of working face mining adopts: Where y is the ordinate of the mining position in the dip direction, and L is the calculated length in the y-axis direction of the working face mining.
Then the calculation model of subsidence for any point on the ground surface is: Eq. (15) combines the IPIM prediction model with the Gompertz function to form the IPIM-G model for predicting dynamic subsidence of mining area. Unknown parameters in the equation include parameters of the IPIM prediction model and the Gompertz time function. After the parameters are determined, the subsidence at any position and at any time on the surface can be predicted according to the above equation.

Predicted parameters of IPIM-G model 3.2.1 IPIM prediction model parameters and determination
IPIM prediction model parameters include settlement coefficient q, main influence angle tangent tan  , inflection point offset distance of downhill, uphill, left boundary and right boundary S 1 , S 2 , S 3 and S 4 , mining propagation influence angle  , horizontal displacement constant b (mining influence coefficient), average mining depth H, loose layer thickness h and loose layer influence coefficient  . The IPIM prediction model only studies the surface subsidence of the mining area, so the parameters to be determined are q, tanβ, S 1 , S 2 , S 3 , S 4 , H, h and  . In the application of the above parameters in different mining areas, there will be different parameter values due to different actual conditions. Therefore, the parameter value provided by the mine is selected for the IPIM prediction model parameters q, tanβ, S 1 , S 2 , S 3 , S 4 , H, and the value of ( / ) Hh  is 0.2r (r is the main influence radius) according to experience.

Gompertz time function parameters and determination
The parameters to be determined for the Gompertz time function include growth rate coefficient k and inflection point c of the function. In the IPIM-G dynamic prediction model, k is the shape parameter, which affects the slope of subsidence curve of each small unit; c is the position parameter, which is the time when each small unit reaches the maximum sinking speed. The parameters k and c need to be inversed by D-InSAR technology. The surface subsidence values are extracted from D-InSAR according to the single line of sight, and the objective function is constructed by combining the surface subsidence values extracted by the D-InSAR technology with the predicted values of the IPIM-G dynamic prediction model, and then the parameter values of Gompertz time function are finally inversed according to genetic algorithm (GA). The specific steps are as follows: 1) Initialize the population. Select the range of parameters to be optimized to generate the initial population, that is, initialize the population. 2) Establish fitness function. The monitoring results of D-InSAR technology are converted into vertical deformation according to Eq. (3), which is recorded as W⊥. According to the parameter values of the IPIM prediction model provided by the actual condition of the mining area and the Gompertz function parameters generated by the initialization population, the IPIM-G dynamic prediction model is substituted to predict the surface subsidence of the mining area, and the predicted results are recorded as ' W . Then the estimated residual (v) is: Sum the squares of these residuals: Finding the minimum value of expected residual can be transformed into finding the minimum value of function f (v). For this, the fitness function is constructed as follows: max( ) fitness vv  (18) 3) Binary encoding. Genetic algorithm is an algorithm constructed according to the principle of genetics. The algorithm needs to encode the population into binary numbers. 4) Fitness calculation. The fitness of the population is calculated according to the fitness function, and then the population is selected according to the calculated fitness. The selected population is crossed and mutated according to the set crossover rate and mutation rate to generate a new generation of population. 5) Iterative computation and decoding. Iterative calculation is to repeat the above 3) and 4) processes. After reaching the set number of iterations, the population is decoded to obtain the final parameter result.

Overview of the study area
The mining impact of Ⅱ513 working face of a coal mine in Huaibei mining area of Anhui Province was selected for the study. There are important buildings such as the drying water tower and machine house of a power plant on the south side of the working face, and Ⅱ311 working face on the right side. The strike length of working face Ⅱ513 is 280m, the dip length is 220m, the average mining depth is 263m, the elevation of working face is -177.7～-280m, the average dip angle of the coal seam is 13°, and the average mining thickness is 3.2m.
Above the working face, there are county-level roads passing through and covered with vegetation and farmland. The hydrogeological conditions are single, containing more sandstone and mudstone, and the texture is soft, which brings some challenges to the monitoring and prediction of mining subsidence of the working face. Meanwhile, working face Ⅱ513 is close to the old goaf, which increases the difficulty of research.

Research data
Sentinel-1A satellite was launched by ESA. The working band of its sensor is C-band, the wavelength is 5.6cm, the polarization mode is VV, the imaging mode is IW, the time baseline of the image is 12d, the azimuth resolution is 20m and the range resolution is 5m. A total of 17 scenes of Sentinel-1A satellite in the study area from 2021-03-12 to 2021-09-20 were selected data images, in which the interference results of 4 images on 2021-03-12, 2021-03-24, 2021-04-05 and 2021-04-17 are used as the input data of model parameter inversion, and the topographic phase in the interference process is removed with the help of SRTM DEM data with 90m resolution. The time-series superposition results of 17 images are used as the D-InSAR monitoring results of mining subsidence.
At the same time, real-time leveling was also carried out in the study area. The leveling observation data is mainly used to compare and verify the prediction results of the dynamic prediction model and the time-series superposition results of D-InSAR monitoring.

Inversion of model parameters
The dynamic prediction parameters of different mining areas will have different values due to different goaf sizes, mining depth, coal seam mining method or hydrogeology. Therefore, the IPIM parameters q, tanβ, S 1 , S 2 , S 3 , S 4 and H were selected as the values provided by the mining area, and the values of Gompertz function parameters k and c suitable for the study area were inversed according to the actual situation of the study area. In the inversion process of this paper, firstly, the obtained three-phase D-InSAR interferometric results were superimposed according to the time sequence to obtain the LOS deformation results from 2021-03-12 to 2021-04-17, which are projected to the vertical direction to obtain the vertical surface deformation map, and the vertical deformation results are recorded as W⊥. As shown in Fig. 5, it can be seen that the D-InSAR monitoring results of the post mining subsidence range of the working face have strong coherence in the area with less vegetation coverage at the edge of the mining subsidence basin of the working face. Select the subsidence feature points at the places with strong coherence, and extract the subsidence values of these selected feature points as the input data of the inversion model parameters. The coordinate system used in the inversion process is different from the working face coordinate system. Before inversion, the selected feature points shall be transformed from the working face coordinate system to the expected coordinate system.
In the process of establishing the coordinate system, the dynamic prediction shall follow the principle that the downhill direction of the coal seam must point to the negative direction of the longitudinal axis y of the coordinate system (as shown in Fig. 2). Therefore, the upper left corner of working face II513 was selected as the coordinate origin of the predicted coordinate system. After the system was established, the coordinates of the feature points were transformed according to Eqs. (18) Where x and y are the coordinates under the projected coordinate system, ' x and ' y are the coordinates under the working surface coordinate system, 0 ' x and 0 ' y are the coordinates of the projected coordinate system origin O under the working surface coordinate system, and  are the angles of counterclockwise rotation of the working surface coordinate system to the working surface coordinate system.
After the coordinate conversion was completed, the number of initial populations of the GA algorithm was set to 100, that is, the parameters k and c have 100 different values respectively within the set range, and then the IPIM-G model was used to predict the subsidence values of the corresponding feature points recorded as ' W . The predicted residuals were calculated according to Eq. (16), the fitness function was constructed according to Eq. (17), the number of iterations is set to 10000, the crossover rate is 0.95 and the variation rate is 0.05, the initial population was encoded with the fitness calculation, and the decoding calculation was performed after the completion of the iteration. The final results of the inversion parameters were obtained as k = 0.0521 and c = 40.0474.

Dynamic prediction of mining subsidence based on inversion parameters
For the model parameters of the inversion using single InSAR interference pairs, in order to verify the applicability of the inversion parameters that can be used for model prediction, the time nodes of 51d, 118d, 146d and 192d after mining of working face II513 were selected to dynamically predict the subsidence caused by mining in this working face.
The dynamic prediction results of vertical subsidence of working face II513 are shown in Fig. 5. It can be seen from the figure that with the advancement of the working face, the center of the subsidence basin continues to move forward, and the working face has large gradient subsidence. Through analysis, it is known that the upper part of the working face is mostly covered with sandstone and mudstone, and the soil is soft, resulting in rapid changes in surface subsidence. At the same time, the scope of the subsidence basin produced by working face mining expands rapidly, and the maximum subsidence value of the subsidence basin increases continuously. The dynamic prediction of working face mining results is in line with the general law of surface subsidence of mine mining. The maximum subsidence in this area reached -3153mm after 192d of working face mining.

Accuracy verification 5.3.1 Prediction accuracy analysis of large gradient subsidence in subsidence basin
The prediction results of the IPIM-G dynamic prediction model and the superposition results of D-InSAR monitoring time-series are extracted and compared with the measured leveling subsidence values of surface deformation observation points arranged on the strike and dip main section of working face II513, which can be used to verify the prediction accuracy of the constructed IPIM-G dynamic prediction model. According to the layout regulations of surface rock movement observation stations in the mining area, the strike and dip observation lines of the surface above the working face can basically monitor the approximate shape of the subsidence basin, and the subsidence accuracy of the two main sections can also reflect the subsidence accuracy of the whole basin. However, in the process of continuous surface subsidence caused by the advance of the working face, some measured observation points on the surface will be damaged due to ponding in the subsidence basin, surface operation or man-made damage, and the actual subsidence can not be obtained. Fig. 6 shows the comparison of the prediction results of the IPIM-G dynamic prediction model and the subsidence between the D-InSAR time-series superposition results and the actual measured values. According to Fig. 6, it can be seen that the prediction results of the IPIM-G dynamic prediction model are basically consistent with the actual subsidence, with high accuracy and good coincidence at the edge of the subsidence basin. The monitoring results of D-InSAR technology at the edge of the subsidence basin are good, but the monitoring results at the center of the subsidence basin are completely incorrect. This is because the severe surface deformation in the center of the subsidence basin leads to serious incoherence in D-InSAR monitoring, so it is impossible to monitor the subsidence in the area with large gradient subsidence in the mining area.
In order to further quantify the error of the prediction results of the IPIM-G dynamic prediction model, the root mean square error (RMSE) of the predicted value of the IPIM-G dynamic prediction model with respect to the measured subsidence value can be calculated. The RMSE calculation formula is: Where W y is the predicted value of IPIM model at point j and W s is the measured leveling value at point j. On the strike and dip observation line, as the subsidence value predicted by the IPIM-G dynamic prediction model increases with the mining time, the increase of subsidence value leads to the gradual increase of RMSE. Maximum error of 239mm and minimum error of 88mm for the strike line and maximum error of 237mm and minimum error of 47mm for the dip line for the prediction results of the IPIM-G dynamic prediction model; The maximum error of the strike line of D-InSAR monitoring results is 1765mm and the minimum error is 421mm, and the maximum error of the dip line is 2077mm and the minimum error is 58mm. Therefore, it can be calculated that the accuracy of the IPIM-G prediction result is 85.9% higher than that of D-InSAR time-series superposition results, which is in good agreement with the measured leveling subsidence value. It has better applicability than D-InSAR technology in extracting large gradient deformation information of mining area, and makes up for the deficiency of D-InSAR Technology to a great extent. The average RMSE predicted by the IPIM-G model is 167mm, which can meet the engineering application requirements of most large gradient settlement mining areas.

Prediction accuracy analysis at the edge of subsidence basin
When the PIM prediction model is applied to predict the subsidence at the edge of the subsidence basin, the prediction effect is not ideal because the PIM prediction model converges rapidly at the edge. In contrast, the IPIM prediction model introduces new parameters and modifies the influence radius of each small unit. The IPIM-G dynamic prediction model combined with the Gompertz function improves the accuracy of prediction results at the edge. The PIM prediction model and IPIM-G dynamic prediction model were used to extract the predicted deformation information of the subsidence basin edge points and compared with the measured subsidence data. As shown in Table 2, the subsidence value of the edge points predicted by IPIM-G is closer to the measured subsidence value of the leveling, and the prediction accuracy at the edge of the subsidence basin is improved.  In order to more clearly compare the accuracy of the IPIM-G dynamic prediction model and the PIM prediction model, the RMSE magnitude of them was calculated. According to Table 3, the prediction accuracy of the IPIM-G dynamic prediction model in each phase is higher than that of the PIM prediction model. The average RMSE of the PIM prediction model for four phases is 7.81mm, and the average RMSE of the IPIM-G model is 6.64mm, which improves the accuracy by 14.9%. The IPIM-G dynamic prediction model achieves higher accuracy in the deformation prediction at the edge of the subsidence basin, which provides a more accurate method for determining the boundary of the subsidence basin in the mining area, and has certain guiding significance for the mining of the working face in the mining area.

Conclusion
In view of the deficiency of conventional D-InSAR in obtaining the information of large gradient subsidence in the mining area and the problem that the PIM prediction model converges too fast at the edge of the subsidence basin, this paper is formed the dynamic prediction model IPIM-G for large gradient subsidence in the mine area by introducing the IPIM prediction model and establishing its relationship with the Gompertz function. The model can not only dynamically predict the large gradient subsidence in the mining area, but also effectively predict the subsidence at the edge of the subsidence basin. The model was applied to the prediction of mining subsidence at working face II513 of a coal mine in Huaibei mining area, and the accuracy was improved by 85.9% compared with the D-InSAR technology in acquiring the information of large gradient subsidence in the mining area, and the accuracy was improved by 14.9% compared with the PIM prediction model in the prediction of subsidence at the edge of the subsidence basin. Because this paper uses C-band low-resolution data, there will be some errors in inversion parameters, and the IPIM-G dynamic prediction model established does not predict the inclination, curvature, horizontal movement and horizontal deformation of the subsidence basin. Subsequently, based on the predicted subsidence value, the inclination, curvature, horizontal movement and horizontal deformation of the subsidence basin can be predicted and analyzed.