# Probability Integral Method Parameter Determination by SBAS-InSAR Technology and GWO Algorithm

Authors: Tieming Liu; Tongkang Zhang; Lichuan Chen; Weiming Liao; Yun Shi; Danchao Xu
Journal: Journal of Sensors
 (2021)

---

## Abstract

This paper proposed a method based on the SBAS-InSAR and gray wolf optimization algorithm aiming at the time-consuming and laborious defects of the traditional method used to obtain the expected parameters of the probability integral method and the shortcomings of the InSAR technology in the field of large gradient deformation detection in the mining area. The fitness function of the algorithm was established based on the geometric relationship between the radar side imaging and the three-dimensional model of the probability integral method. The stable sinking point of the settlement boundary obtained by SBAS-InSAR was used as the input value for the calculation of the predicted parameters of the probability integral method. Firstly, the simulation experiment was employed for the simulation of the direction of the InSAR line of sight combined with the geological mining conditions of the assumed working face, thereby obtaining the probability integral prediction parameters of the working face. Consequently, the maximum relative error ofq, b, tanβ, and θ0 does not exceed 8%, and that of S1, S2, S3, and S4 does not exceed 35.5% (low parameter sensitivity). The error of the LOS-direction deformation fitting is 0.076 m, which meets the tolerance requirements, and the result is trustworthy. At last, the parameter finding method is applied to the engineering example, that is, the 112201 working face of Xiaobaodang Coal Mine in the northern Shaanxi mining area. The settlement value of the stable boundary point is obtained based on the SBAS-InSAR results, which is substituted into the fitness function. And the GWO optimization algorithm is used for optimization and parameter finding; the probability integral expected parameters of the working face are calculated as q=0.63, b=0.37, tanβ=2.76, θ0=83.94, S1=−36.34m, S2=26.69m, S3=−45.64m, and S4=39.62m. Substitute the obtained parameters into the probability integral model for the prediction of the vertical and horizontal displacements of the working face, and verify its accuracy with the GPS measured data. The results showed that the maximum absolute error of vertical displacement reached 116 mm, the median error was 63 mm, and the maximum absolute error of north-south horizontal movement reached 56 mm; meanwhile, the median error was 23 mm, the maximum absolute error of east-west horizontal movement reached 61 mm, and the median error was 29 mm; all the above parameters are within the tolerance range, indicating that the method for the calculation of probability integral parameters proposed in this paper is applicable in actual engineering.


---

## Body

## 1. Introduction

As one of China’s main energy sources, coal has always been a pillar industry of the lifeline of the national economy [1]. China is rich in coal resources, and the proven coal reserves rank the third in the world [2], and coal will continue to be China’s main energy source for a certain period in the future. It is undeniable that the development and utilization of coal resources have played an important role in the development of China’s economy. However, the environmental and geological problems accompanying the mining process are becoming more and more serious. The main environmental problems lie in geological disasters, including ground subsidence, ground fissures, landslides, and mudslides caused by mining [3]. In addition, the movement of rock formations caused by mining often results in the damages of roads, pipelines, villages, and any other infrastructure located in subsidence areas. In addition, it may also cause soil erosion on the ground [4], destroy the stratum structure, and seriously affect the local economic development. Among them, as the most serious problem, the ground subsidence will directly result in the destruction of surface houses, bridges, and other infrastructures, thereby posing a great threat to the safety of lives and property of the residents. In that case, the accurate prediction of these displacements is both necessary and helpful for the assessment of potential damage and the improvement of the design of coal mining. In recent decades [5], the techniques, such as empirical prediction methods [6] and influence function methods, have been developed, which often require a large number of repeated and accurate three-dimensional displacement measurement samples to build a reliable prediction model before performing forward displacement prediction [7]. However, to monitor ground displacement caused by mining. The traditional geodetic methods, such as leveling and global positioning system measurements, are usually based on the measurement at only a few sparse points, which makes it quite difficult to use these sparse point measurements to describe the entire subsidence area [8]. In addition, these methods are costly, time-consuming, and highly labor-intensive which often result in unreliable prediction models and inaccurate forward predictions. Interferometric Synthetic Aperture Radar (InSAR) technology can play a certain role and greatly improve the reliability of predicting mining displacement in an economically feasible way due to its low cost, high accuracy and efficiency, large spatial coverage area, high temporal and spatial resolution [9], and the availability to measure surface deformations related to geophysical and human activities. However, due to the fact that the mining area subsidence is often accompanied by a large number of large gradients, the detection capability of InSAR will be constrained. This is because that the impact of spatial miscorrelation is relatively large [10], resulting in hollow area appeared in InSAR interferograms. This is also a big limitation of the current mining area deformation monitoring by InSAR. However, because the subsidence basin is characterized by small boundary and large center in terms of magnitude, InSAR can detect the settlement information of stable boundary points [11].To better protect the safety of surface buildings and facilities, it is necessary to predict and protect the surface movement caused by mining. At present, the probability integral method is the most widely used mining subsidence prediction method [12], whose prediction quality is mainly determined by the accuracy of the predicted parameters [13]. Because in the probability integral method, the predicted parameters are affected by the complex geological factors and mining conditions; the expression with simple mathematical methods is not available. Moreover, the obtaining of parameters through actual measurement based on the establishment of a mobile surface observation station is also easily restricted by a variety of conditions, and it is difficult to meet the requirements of safe production in mining areas. With the development of computer technology, various artificial intelligence algorithm methods are mastered and applied to the field of probabilistic integration to obtain parameters, which mainly include neural network model [14], support vector machine [15], genetic algorithm [16] (GA), and Particle Swarm Optimization (PSO) [17]. A large number of experiments have proved that the utilization of machine learning to obtain the parameters of the probability integral method shows the advantages of high accuracy, easy operation, and strong practical effect. Gray wolf optimization (GWO) [18] is a population intelligent optimization algorithm, which is mainly characterized by simple operation, few parameters, and easy implementation. Article [19] improved the convergence factor strategy and introduced a dynamic weight strategy to optimize the performance of the gray wolf algorithm; Article [20] combined the advantages of differential evolution and the gray wolf algorithm and proposed a hybrid algorithm for three-dimensional chip testing to further balance the algorithm’s exploration and mining capabilities; Article [21] proposed a hybrid gray wolf optimization algorithm based on chaos principle and elite reverse learning strategy to solve the problem of high-dimensional function optimization; Article [22] embeds an optimal optimization algorithm. The difference reverse learning strategy and a dynamic random difference mutation operator are used to improve the gray wolf optimization algorithm to improve the clustering optimization performance; Article [23] proposed a hybrid algorithm of gray wolf optimization and locust optimization to solve the text mining problem. Article [24] uses the GWO algorithm to invert the parameters of the wind turbine power model. Therefore, this paper constructs the fitness function of the gray wolf optimization algorithm based on the probabilistic integral method model and the principle of InSAR side-view imaging and proposes a probabilistic integral parameter calculation method that combines SBAS-InSAR and GWO algorithms. This method combines the advantages of InSAR technology and machine learning algorithm and reverses the predicted parameters of the probability integral method. The method is novel and provides new method support for the parameter inversion of the probability integral method and the field of mining subsidence prediction.


## 2. Introduction to Gray Wolf Optimization Algorithm

Gray wolf optimization (GWO) [18] is a population intelligent optimization algorithm derived from the wolf pack search and hunting activities. The wolf pack has a strict hierarchical distribution relationship used to determine the optimal parameter value in the optimization process. Similarly, this paper uses the global search ability of the GWO algorithm in the process of searching for prey to perform global optimization, so as to obtain the predicted parameters with an optimal probability integral method. The specific establishment steps are as follows.Firstly, it is assumed that the population size of the wolf pack consists ofn individuals, and then, the position of the ith wolf is determined, where it refers to the d-dimensional position of the ith wolf in a multidimensional space, and the first three optimal solutions of the objective function value are provided. There are three wolves with the highest rank, namely, αwolf, βwolf, and δwolf. And the three wolves are mainly responsible for hunting. Besides, the hunting behavior can be described by the following formula:
(1)D=QXPt−Xt,(2)Q=2r1.Formula (1) represents the distance between the individual and the prey, and formula (2) is the position update formula of the gray wolf. In the above formula, t is the number of iterations, A and Q are coefficient vectors, and XPt and Xt, respectively, represent the positions of the prey and gray wolf after t iterations. The calculation formula of A and Q is
(3)Xt+1=XPt−AD,(4)A=2ar2−a,(5)a=2−2×tT.In the above formula,a is the convergence factor, as the number of iterations decreases linearly from 2 to 0, t is the current iteration number, and T is the maximum iteration. r1 and r2 take random numbers between [0,1].The most important thing for wolves to do in hunting is to be able to determine the location of their prey, so we useαwolves, βwolves, and δwolves to search for the potential location of the prey, and lead other search individuals to constantly update the potential location according to these three optimal solutions, thereby obtaining the position of the prey. The specific steps are as follows:
(6)Dα=Q1XPt−Xt,(7)Dβ=Q2XPt−Xt,(8)Dδ=Q3XPt−Xt,(9)X1t+1=Xαt−A1Dα,(10)X2t+1=Xβt−A2Dβ,(11)X3t+1=Xδt−A3Dδ,(12)Xpt+1=X1+X2+X33,and represent the distance between the position of the candidate wolf and that of the best three wolves, respectively. WhenA>1, gray wolves will spread out to search for prey, while when A<1, gray wolves will gather together to search for prey.


## 3. Probability Integral Model and InSAR Geometric Relationship to Establish Fitness Function

SBAS-InSAR technology is available to monitor and quantify the ground displacement that occurred in adjacent time intervals. Since the radar is based on side-view imaging, it can only measure the path length difference in the line of sight (LOS) direction [25], instead of the true three-dimensional deformation of the surface. To quantify the relationship between the actual three-dimensional ground displacement and the LOS deformation measured by SBAS-InSAR, it is assumed that the east-west, north-south, and up-and-down displacement components of the actual ground are UE, UN, and W, respectively, and then, we project UE, UN, and W to radar LOS direction [26]. The ascending satellite orbit with the azimuth angle αh refers to the angle between the north and the orbit direction (clockwise).As shown in Figure1, the LOS deformation is jointly contributed by the real three-dimensional ground displacement. In the case that the three-dimensional ground displacement in mining engineering can be expressed by the model parameters of the mechanical model, certain functional relationship exists between the model parameters based on equation (15) and the LOS deformation. Developed from the random medium theory, and derived from the probability density function in the prediction equation, PIM presents one of the most commonly used influence function methods in underground mining [27]. According to PIM theory, the vertical settlement Wx,y and horizontal displacement Ux,y,ϕ caused by underground mining can be expressed as
(13)Wx,y=W0xW0yW0,UN=U0xW0ycosφN+U0yW0xsinφN,UE=U0xW0ycosφE+U0yW0xsinφE,among them,
(14)W0=mqcosα,W0x=W′x,r−W′x−l,r,W′x,r=W02erfπrx+1,W0y=W′y,r1−W′y−L,r2,U′x,b,r=bW0exp−πx2r2,U0x=U′x,b,r−U′x−l,b,r,U0y+U′y,b,r+W′y,rcotθ0−U′y−L,b,r−W′y−L,rcotθ0.Figure 1
Schematic diagram of radar side-view imaging.Then, project to the LOS direction through the projection formula:(15)LOS=cosθ−sinθcosαsinθsinαWUEUN,where θ refers to the incident angle of the radar, and LOS represents the deformation derived from SBAS-InSAR, which is defined as the direction of the line of sight of the satellite. W0 represents the maximum possible settlement; q denotes settlement factor; erf represents the probability density.
(16)erfπxr=2π∫0π/rxe−u2du.r=H/tanβ is the main influence radius of the mining area, where tanβ stands for the tangent of the main influence angle; b denotes the horizontal displacement constant.
(17)l=D3−S3−S4,L=D1−S1−S2sinθ0+αsinθ0.D1 and D3 represent the length of the working face’s inclination and strike direction, respectively, where S1, S2, S3, and S4 represent the deviation of the inflection point in the up, down, left, and right directions, respectively, θ0 refers to the mining propagation angle; φ represents the clamp between the given direction and the forward direction of the work panel angle (counterclockwise) [28].The fitness function of the GWO algorithm can be expressed as follows:(18)P=minLOS0x,y−LOSpx,y,where x and y refer to the vectors containing the pixel coordinates of the LOS displacement (indicated by LOS0) measured by the SBAS-InSAR technology, and LOSp represents the predicted ground displacement in the LOS direction of the radar with a predetermined functional relationship. The 8 optimal parameter values of the probability integral method are used as the optimal position to be searched for by the gray wolf.


## 4. Simulation Experiment

### 4.1. Simulation of Geological Mining Conditions

Based on the coal measure strata in the northern Shaanxi mining area as the background, the overlying rock lithology of the simulated working face is set as medium hard, and the horizontal projections ofD1 and D3 of the mining size along the strike and the inclination direction are 200 m and 800 m, respectively. Besides, the mining height H is set as 300 m, the coal seam inclination angle is 2°, and the mining thickness M=6 m; the strike long-arm mining method is adopted for mining, and the all caving method is adopted for roof management. Because D1/H≈0.66<1.2 and D3/H≈2.66>1.2 in both the strike and inclination of the working face, insufficiency mining is used in the inclination, sufficiency mining is used in the strike of the working face, and in general, insufficiency mining is adopted [29]. The characteristics of the subsidence basin in the northern Shaanxi mining area conform to the probability integral method model. Therefore, in this simulation experiment, the main predictive parameters of the probability integral method are set as subsidence coefficient q=0.8, horizontal movement coefficient b=0.25, and tangent value of the main influence angle tanβ=2, maximum sinking angle θ0=85°, and inflection point offset S1=S2=S3=S4=20 m.


### 4.2. InSAR Line-of-Sight Strain Simulation

The main SAR data used in this paper is the Sentinel-1A data of ESA. When the satellite is over northern Shaanxi for geometric imaging, the satellite’s incident angleθ=39.13°, and the orbital inclination angle αh=346.7°. The pixel size is selected according to the spatial resolution of the SAR data during simulation. According to the set probability integral method model parameters, the subsidence value W, the horizontal movement value UN in the north-south direction, and the horizontal movement value UE in the east-west direction of the subsidence basin boundary point are predicted. Then, the orbit parameters of the SAR satellite and the projection relationship are used to calculate the LOS-direction deformation value of each point, with the simulation effect shown in Figure 2.Figure 2
Simulated LOS to deformation value.


### 4.3. Analysis of Experimental Results

The boundary point of the LOS deformation value of the upper simulated working face was extracted as the input value of the algorithm fitness function, and the probability integral prediction parameter calculation model of the InSAR fusion gray wolf optimization algorithm was established for the calculation of the probability integral model parameters under the conditions of the mining area, thereby eliminating the influence of gross errors in the calculation process. Five experiments were carried out on the algorithm, and the average value was calculated to obtain a reliable result.The specific steps for obtaining the parameters are as follows:Step 1. Input the initial parameters of PIM. According to the geological mining conditions of the Loess mining area, the initial values of each parameter are set asq=0.8, b=0.25, tanβ=2, θ0=85, S1=20, S2=20, S3=20, and S4=20.Step 2. Determine the initial parameters of the differential gray wolf optimization algorithm. The parameter values set in this paper are 30 populations, 500 maximum iterations, and 8 independent variable dimensions; besides, the boundaries are set.Step 3. Determine the initial position of the population, calculate the value of the objective function in the parent population, and determine the three best individuals in the population asαwolf, βwolf, and δwolf which are denoted as Xα, Xβ, and Xδ.Step 4. Use equations (6)–(8) to find the distance between other individuals and the best three individuals, and then, update the position of each individual in the population according to equations (9)–(11).Step 5. Iterate, calculate, and update the values of parametersQ and A according to formulas (2), (4), and (5), and update the position of Xα, Xβ, and Xδ to prevent border crossing.Step 6. Update the positions of the optimalαwolf, βwolf, and δwolf in the parent population.Step 7. Determine the number of iterations; in the case that the number of iterations reaches the maximum value, select the optimal value of PIM parameter; otherwise, return to Step 4 to continue.The results of parameter calculation are shown in Table1.Table 1
GWO expected parameter inversion results.
Serial numberqbtanβθ0S1 (m)S2 (m)S3 (m)S4 (m)10.820.311.8882.9519.4513.65-14.318.4820.860.252.1386.1217.2916.22-14.9119.4230.860.271.6184.0222.1914.74-12.8116.2640.830.262.0683.2223.812.81-9.3413.7550.850.282.0885.2624.5812.65-13.1517.47Average0.840.271.9584.3121.4614.01-12.9015.08Parameter true value0.80.2528520202020Absolute error0.040.02-0.05-0.691.46-5.997.1-4.92Relative error%58-2.5-0.817.3-29.9535.5-24.6The effect of parameter fitting is shown in Figure3.Figure 3
Comparison of simulated LOS deformation and fitted LOS deformation.As shown in Table1, the results of the probability integral predicted parameter obtained by the LOS deformation based on InSAR combined with the gray wolf optimization algorithm are more accurate. Among them, the maximum relative error of q, b, tanβ, and θ0 fails to exceed 8%, and that of the inflection point offset S1-S4 does not exceed 35.5% (the parameter sensitivity is low and has little effect on the overall effect). Figure 3 illustrates the LOS fitting diagram in the process of parameter obtaining. According to the figure, the algorithm shows a better overall fitting effect on the simulated values of the InSAR boundary points, and the trends are more consistent. The residual error range fluctuates from -0.13 m to 0.109, the median error is 0.076 m, and the result is trustworthy. It suggests that the method to obtain the parameters of the probability integral combining the SBAS-InSAR and GWO algorithms is applicable in practice.


## 4.1. Simulation of Geological Mining Conditions

Based on the coal measure strata in the northern Shaanxi mining area as the background, the overlying rock lithology of the simulated working face is set as medium hard, and the horizontal projections ofD1 and D3 of the mining size along the strike and the inclination direction are 200 m and 800 m, respectively. Besides, the mining height H is set as 300 m, the coal seam inclination angle is 2°, and the mining thickness M=6 m; the strike long-arm mining method is adopted for mining, and the all caving method is adopted for roof management. Because D1/H≈0.66<1.2 and D3/H≈2.66>1.2 in both the strike and inclination of the working face, insufficiency mining is used in the inclination, sufficiency mining is used in the strike of the working face, and in general, insufficiency mining is adopted [29]. The characteristics of the subsidence basin in the northern Shaanxi mining area conform to the probability integral method model. Therefore, in this simulation experiment, the main predictive parameters of the probability integral method are set as subsidence coefficient q=0.8, horizontal movement coefficient b=0.25, and tangent value of the main influence angle tanβ=2, maximum sinking angle θ0=85°, and inflection point offset S1=S2=S3=S4=20 m.


## 4.2. InSAR Line-of-Sight Strain Simulation

The main SAR data used in this paper is the Sentinel-1A data of ESA. When the satellite is over northern Shaanxi for geometric imaging, the satellite’s incident angleθ=39.13°, and the orbital inclination angle αh=346.7°. The pixel size is selected according to the spatial resolution of the SAR data during simulation. According to the set probability integral method model parameters, the subsidence value W, the horizontal movement value UN in the north-south direction, and the horizontal movement value UE in the east-west direction of the subsidence basin boundary point are predicted. Then, the orbit parameters of the SAR satellite and the projection relationship are used to calculate the LOS-direction deformation value of each point, with the simulation effect shown in Figure 2.Figure 2
Simulated LOS to deformation value.


## 4.3. Analysis of Experimental Results

The boundary point of the LOS deformation value of the upper simulated working face was extracted as the input value of the algorithm fitness function, and the probability integral prediction parameter calculation model of the InSAR fusion gray wolf optimization algorithm was established for the calculation of the probability integral model parameters under the conditions of the mining area, thereby eliminating the influence of gross errors in the calculation process. Five experiments were carried out on the algorithm, and the average value was calculated to obtain a reliable result.The specific steps for obtaining the parameters are as follows:Step 1. Input the initial parameters of PIM. According to the geological mining conditions of the Loess mining area, the initial values of each parameter are set asq=0.8, b=0.25, tanβ=2, θ0=85, S1=20, S2=20, S3=20, and S4=20.Step 2. Determine the initial parameters of the differential gray wolf optimization algorithm. The parameter values set in this paper are 30 populations, 500 maximum iterations, and 8 independent variable dimensions; besides, the boundaries are set.Step 3. Determine the initial position of the population, calculate the value of the objective function in the parent population, and determine the three best individuals in the population asαwolf, βwolf, and δwolf which are denoted as Xα, Xβ, and Xδ.Step 4. Use equations (6)–(8) to find the distance between other individuals and the best three individuals, and then, update the position of each individual in the population according to equations (9)–(11).Step 5. Iterate, calculate, and update the values of parametersQ and A according to formulas (2), (4), and (5), and update the position of Xα, Xβ, and Xδ to prevent border crossing.Step 6. Update the positions of the optimalαwolf, βwolf, and δwolf in the parent population.Step 7. Determine the number of iterations; in the case that the number of iterations reaches the maximum value, select the optimal value of PIM parameter; otherwise, return to Step 4 to continue.The results of parameter calculation are shown in Table1.Table 1
GWO expected parameter inversion results.
Serial numberqbtanβθ0S1 (m)S2 (m)S3 (m)S4 (m)10.820.311.8882.9519.4513.65-14.318.4820.860.252.1386.1217.2916.22-14.9119.4230.860.271.6184.0222.1914.74-12.8116.2640.830.262.0683.2223.812.81-9.3413.7550.850.282.0885.2624.5812.65-13.1517.47Average0.840.271.9584.3121.4614.01-12.9015.08Parameter true value0.80.2528520202020Absolute error0.040.02-0.05-0.691.46-5.997.1-4.92Relative error%58-2.5-0.817.3-29.9535.5-24.6The effect of parameter fitting is shown in Figure3.Figure 3
Comparison of simulated LOS deformation and fitted LOS deformation.As shown in Table1, the results of the probability integral predicted parameter obtained by the LOS deformation based on InSAR combined with the gray wolf optimization algorithm are more accurate. Among them, the maximum relative error of q, b, tanβ, and θ0 fails to exceed 8%, and that of the inflection point offset S1-S4 does not exceed 35.5% (the parameter sensitivity is low and has little effect on the overall effect). Figure 3 illustrates the LOS fitting diagram in the process of parameter obtaining. According to the figure, the algorithm shows a better overall fitting effect on the simulated values of the InSAR boundary points, and the trends are more consistent. The residual error range fluctuates from -0.13 m to 0.109, the median error is 0.076 m, and the result is trustworthy. It suggests that the method to obtain the parameters of the probability integral combining the SBAS-InSAR and GWO algorithms is applicable in practice.


## 5. Case Analysis

### 5.1. Overview of the Study Area

Xiaobaodang Minefield is located in the southwest of Shenmu County, Yulin City, Shaanxi Province, bordering the Inner Mongolia Plain in the northwest, facing the Yellow River in the east and Xing County in Shanxi Province across the river, and adjacent to Jiaxian County, Yulin City, Shaanxi Province in the south, which is under the jurisdiction of Dabaodang Town, Shenmu County, Yulin City. Its geographical coordinates are as follows: east longitude 109°49′09.6553″ ~ 109°59′41.9675″; north latitude 38°39′15.8090″ ~ 38°47′32.0113″. The Xiaobaodang Minefield is delineated by the connecting lines of four inflection points, with a width of 7.7 km from north to south and a length of 15.2 km from east to west, covering an area of 97.18 km2. The traffic around Xiaobaodang Minefield is well developed, specifically, the Yushen Class II Highway runs in the southeast of the minefield, and the Baomao Expressway and the 210 National Highway pass through the west. Xiaobaodang Mine Field is located at the northern end of the Loess Plateau in northern Shaanxi, on the southeastern edge of the Mu Us Desert. The area is eroded by the Quaternary aeolian sand all year round, with most of the area covered by fixed dunes and semifixed dunes. The landforms are mainly characterized by wind-eroded aeolian desert hilly landforms. From the topographical point of view, the overall trend is higher in the southwest and lower in the east. The highest point is 1346.60 m. Shibuliang is located in the south. The lowest point in the east is 1198 m east. The elevation range is between 1250 and 1300 m, and the relative maximum height difference is 148.60 m.The 112201 working face of Xiaobaodang Coal Mine is shown as the black rectangle mark in Figure4. The 112201 working face was in operation from September 2018 to November 2019 adopting long-arm mining, and the all-caving method for roof management. Besides, sufficient mining was used in the strike direction of the working face, and insufficient mining was utilized in the inclination direction; in general, insufficient mining was used. The geometric parameters are summarized as follows: average mining depth H≈340 m, average mining thickness m≈6 m, mining length D3≈4560 m, mining width D1≈350 m, and dip angle α≈1°. Two monitoring lines, i.e., Z and Q, are arranged along the direction of the working face and the main section of the inclination; besides, 37 monitoring points are arranged in the strike direction, 40 monitoring points are arranged in the inclined direction, and 77 monitoring points are arranged in total. The monitoring line is shown in Figure 2. Based on the constructed model and the simulated mining subsidence prediction parameters (Table 1), the subsidence and horizontal movement after the mining subsidence and stabilization are predicted.Figure 4
Location map of the study area.The location of the study area is shown in Figure4.


### 5.2. SBAS-InSAR Subsidence Basin Boundary Information Extraction

#### 5.2.1. Data Processing

The experimental data are mainly based on the data of Sentinel-1A of ESA. The Sentinel satellite is an earth observation satellite launched by ESA in 2014, which is equipped with a C-band synthetic aperture radar, and available for all-weather earth observation throughout the day. The revisit period of a single constellation is 12 days, in which the interference width can reach 250 km, and the ground resolution is (5∗20) m. In this paper, 41 scenes of the ascending orbit were extracted from the archived data between September 3, 2018, and January 8, 2020.The data processing of SBAS-InSAR is mainly divided into two parts, and the first part is data preprocessing. The preprocessing part mainly includes the extraction of single-view complex image (SLC) data, and the selection of super master images for subpixel precision registration of the data. The second step is to remove the error phase. Firstly, set the space baseline and time baseline as 200 m and 33 d, respectively, according to the geological conditions of the study area. After phase unfolding, the maximum spatial baseline is 113.7 meters, and the minimum is 2.8 meters. As shown in Figure5, select 20190113 as the main image. The generation of the interferogram is the most critical step; in this process, the GoIdstein filter is used to improve the phase expansion accuracy and measurement accuracy, thereby maximizing the signal-to-noise ratio of the interferogram. The minimum cost flow method is employed to unwrap the phase, and the coherence threshold is set at 0.2. With the aid of precision orbits and DEM, flat ground, phase, and noise are removed, and a high-coherence region unwrapping map is obtained. The coherence coefficient method is used for the extraction of high-coherence points, and 62 high-quality interference pairs are selected. The linear model equations of linear deformation rate and elevation error are established in the region of high coherence point, and the SVD algorithm is used to solve the model parameters, thereby obtaining the linear deformation phase. After the removal of the linear deformation phase, the residual error contains atmospheric phase, nonlinear deformation, and noise. The residual can be subjected to high-pass filtering in the time domain and low-pass filtering in the spatial domain according to the characteristics of each component, thereby removing the atmospheric phase components and noise to obtain a nonlinear deformation phase. Finally, the least square method is utilized to obtain the deformation rate and cumulative deformation value of each coherent point.Figure 5
SBAS-InSAR time-space baseline diagram.


#### 5.2.2. Result Analysis

Figure6 shows the annual average settlement rate of the study area. It can be seen from the figure that the center of the subsidence basin has been seriously out of coherence, resulting in the lack of most of the settlement information in the center of the working face, and only the marginal settlement with a smaller magnitude is still preserved better. The settlement magnitude ranges from -0.105 m to 0.035 m, which is relatively small. The blue line is the two cross-sections of the direction and the inclination. The ground GPS observation points are arranged along the direction and the inclination, respectively, which are the green points in the figure. The direction point number is Z1-Z37, and the inclination point number is Q1-Q40. There are labels in the figure. In order to further study the time series settlement trend of the boundary points, the time series diagram obtained by SBAS-InSAR is derived as shown in Figure 7.Figure 6
The annual average settlement rate of the study area.Figure 7
SBAS-InSAR monitoring time series deformation diagram.As shown in Figure7, the black area refers to the working surface, the blue line represents the two section lines along the strike and the inclination, and the black mark points around the working surface are the extracted boundary points. The boundary points of B12 and B14 coincide with the points of Z2 and Q3 on the oblique observation line in the strike and inclination direction, which can be used to detect the observation accuracy of boundary points. During the study period, the 112201 working face was being mined, so the settlement area was mainly concentrated in the 112201 working face. It can be seen from the time sequence diagram that due to the unique geological conditions of the study area, the surface deformation due to coal mining in the Loess Plateau area is characterized by large gradients and large scale, which causes the spatial decorrelation of the SAR interferogram. In that case, high-coherence points are rarely selected above the working surface, which results in the absence of goaf center in the InSAR topographic map. However, based on the basin effect of the goaf surface subsidence and the law of surface deformation of mining subsidence, it can be found that the mining center of the mining area has the characteristics of large deformation speed. However, the basin boundary is a slow linear deformation point, which is in line with the deformation point that can be detected by InSAR. Therefore, only the time series settlement of the boundary point is realized in the timing diagram, and the mining center is not reflected in the InSAR results. Therefore, in this paper, we only discuss the deformation law of the boundary point of the subsidence basin. It can be seen that the settlement funnel extends from southeast to northwest, which is consistent with the mining direction. To verify the monitoring accuracy of SBAS-InSAR, B12 and B14 are utilized to verify their accuracy. Due to the fact that the deformation in the LOS direction is obtained by InSAR, while the GPS observation obtains the vertical and horizontal displacement, firstly, the displacements in three directions, i.e., the vertical direction, east-west, and north-south, obtained by the GPS observation are projected to the LOS upward, and the projection formula is shown in equation (15).As shown in Figure8, the GPS measurement results are projected to the LOS backward and compared with SBAS-InSAR. The two GPS observations at points B12 and B14 in the early stage of mining are two. Comparing the results under the conditions that the time and space references are unified, it can be seen that the GPS measured values at points B12 and B14 are slightly larger than the InSAR monitoring value. And the maximum residual error is less than 5 mm, while the maximum residual error at point Q3 is 7 mm. It can be seen that the overall fit between the InSAR results and the GPS results is good, the correlation is high, and the residual line fluctuates in a small range near the 0 value; besides, the fluctuations are stable. It shows that the results of SBAS-InSAR at boundary points are trustworthy, which can be used to invert model parameters in mining subsidence.Figure 8
GPS/InSAR comparison results of boundary points.


#### 5.2.3. Analysis of GWO Algorithm Parameter Inversion Results

Since the mining in 112201 working face was started in September 2018 and ended in November 2019, the subsidence basin was basically stable till the last scene selected in this paper during the imaging in January 2021. Therefore, the cumulative settlement in the last scene of the above SBAS-InSAR is used as the LOS input value of the fitness function, and the error fitting diagram of the fitness function is shown in Figure9.Figure 9
LOS fitting residual plot of fitness function.Due to the instability of the algorithm, the parameter value will be affected by the gross error. To mitigate the influence of these errors and the instability of the algorithm, we especially solve the parameters for 5 times and find the average value. The result of the solution is shown in Table2.Table 2
GWO expected parameter inversion results.
Serial numberqbtanβθ0S1 (m)S2 (m)S3 (m)S4 (m)10.620.352.3685.12-39.0624.78-46.6241.2520.580.392.9383.19-35.1526.16-44.8937.5430.650.373.0583.24-34.9428.57-49.1843.1640.670.372.7483.22-37.9926.96-38.5233.7650.640.352.7184.92-34.5827.00-48.9942.39Average0.630.372.7683.94-36.3426.69-45.6439.62As shown in the table above, the prediction results of the probability integral method areq=0.63, b=0.37, tanβ=2.76, θ0=83.94, S1=−36.34 m, S2=26.69 m, S3=−45.64 m, and S4=39.62 m. To further verify the accuracy of the calculated parameters, we substituted the result parameters into the probability integral model, combined with the geological and mining conditions of the working face, to determine the vertical and horizontal displacements of the observation points on the observation line in the strike direction and the inclination direction of the working face. In addition, predictions and analyses of the accuracy of each observation point based on the GPS measured data were made, with the results shown as follows.It is found in Figure10 that the surface movement obtained by the inversion is in good agreement with the measured results. The point with the largest absolute error of sinking is Point Q24 on the oblique observation line, the absolute error of this point reaches 116 mm, and the median error is 63 mm. The point with the largest absolute error of north-south horizontal movement is Point Z15, with the error of 56 mm, and the median error of 23 mm. The point with the largest absolute error of east-west horizontal movement is Point Z34, with the maximum absolute error of 61 mm, and the median error of 29 mm, both of them are within the tolerance range. It illustrates that the inversion of the probabilistic integration method parameters based on the SBAS-InSAR and gray wolf optimization algorithm used in this paper is feasible, which is available to obtain more reliable surface movement probabilistic integration method parameters.Figure 10
Comparison between the predicted result of probability integral and the measured result of GPS.
(a)
Accuracy comparison chart in the vertical direction of the observation point(b)
Accuracy comparison chart for the north and south of the observation point(c)
Accuracy comparison chart in the east-west direction of the observation pointSince the prediction results of a small number of points on the edge of the subsidence basin with smaller settlement magnitudes are not representative, we only consider the points with settlement magnitudes greater than 100 mm in our accuracy evaluation and analyze the accuracy of the results, mainly from the absolute error (AE), relative error (RE), and average absolute percentage error (MAPE) evaluated. And a one-samplet test with a significance of 0.05 and a test value of 60 was done for the sample point error. The specific test results are shown in Table 3.Table 3
Statistical test table of forecast results.
DirectionMAX.AE (m)MAX.RE (%)Mape (%)sigConfidence intervalVertical0.1163590.015[-27,-3]East-west0.0562570[-33,-27]North-south0.0612180[-27,-21]The error is tested and found to be in accordance with the normal distribution. A one-samplet test with a significance of 0.05 is performed. The test value is set to 50. The significance levels in the three directions of the test are 0.015, 0, and 0, indicating rejection null hypothesis. There is a 95% probability that the error value is less than 50, and the confidence interval is estimated. Further, verify the accuracy of the predicted value.


## 5.1. Overview of the Study Area

Xiaobaodang Minefield is located in the southwest of Shenmu County, Yulin City, Shaanxi Province, bordering the Inner Mongolia Plain in the northwest, facing the Yellow River in the east and Xing County in Shanxi Province across the river, and adjacent to Jiaxian County, Yulin City, Shaanxi Province in the south, which is under the jurisdiction of Dabaodang Town, Shenmu County, Yulin City. Its geographical coordinates are as follows: east longitude 109°49′09.6553″ ~ 109°59′41.9675″; north latitude 38°39′15.8090″ ~ 38°47′32.0113″. The Xiaobaodang Minefield is delineated by the connecting lines of four inflection points, with a width of 7.7 km from north to south and a length of 15.2 km from east to west, covering an area of 97.18 km2. The traffic around Xiaobaodang Minefield is well developed, specifically, the Yushen Class II Highway runs in the southeast of the minefield, and the Baomao Expressway and the 210 National Highway pass through the west. Xiaobaodang Mine Field is located at the northern end of the Loess Plateau in northern Shaanxi, on the southeastern edge of the Mu Us Desert. The area is eroded by the Quaternary aeolian sand all year round, with most of the area covered by fixed dunes and semifixed dunes. The landforms are mainly characterized by wind-eroded aeolian desert hilly landforms. From the topographical point of view, the overall trend is higher in the southwest and lower in the east. The highest point is 1346.60 m. Shibuliang is located in the south. The lowest point in the east is 1198 m east. The elevation range is between 1250 and 1300 m, and the relative maximum height difference is 148.60 m.The 112201 working face of Xiaobaodang Coal Mine is shown as the black rectangle mark in Figure4. The 112201 working face was in operation from September 2018 to November 2019 adopting long-arm mining, and the all-caving method for roof management. Besides, sufficient mining was used in the strike direction of the working face, and insufficient mining was utilized in the inclination direction; in general, insufficient mining was used. The geometric parameters are summarized as follows: average mining depth H≈340 m, average mining thickness m≈6 m, mining length D3≈4560 m, mining width D1≈350 m, and dip angle α≈1°. Two monitoring lines, i.e., Z and Q, are arranged along the direction of the working face and the main section of the inclination; besides, 37 monitoring points are arranged in the strike direction, 40 monitoring points are arranged in the inclined direction, and 77 monitoring points are arranged in total. The monitoring line is shown in Figure 2. Based on the constructed model and the simulated mining subsidence prediction parameters (Table 1), the subsidence and horizontal movement after the mining subsidence and stabilization are predicted.Figure 4
Location map of the study area.The location of the study area is shown in Figure4.


## 5.2. SBAS-InSAR Subsidence Basin Boundary Information Extraction

### 5.2.1. Data Processing

The experimental data are mainly based on the data of Sentinel-1A of ESA. The Sentinel satellite is an earth observation satellite launched by ESA in 2014, which is equipped with a C-band synthetic aperture radar, and available for all-weather earth observation throughout the day. The revisit period of a single constellation is 12 days, in which the interference width can reach 250 km, and the ground resolution is (5∗20) m. In this paper, 41 scenes of the ascending orbit were extracted from the archived data between September 3, 2018, and January 8, 2020.The data processing of SBAS-InSAR is mainly divided into two parts, and the first part is data preprocessing. The preprocessing part mainly includes the extraction of single-view complex image (SLC) data, and the selection of super master images for subpixel precision registration of the data. The second step is to remove the error phase. Firstly, set the space baseline and time baseline as 200 m and 33 d, respectively, according to the geological conditions of the study area. After phase unfolding, the maximum spatial baseline is 113.7 meters, and the minimum is 2.8 meters. As shown in Figure5, select 20190113 as the main image. The generation of the interferogram is the most critical step; in this process, the GoIdstein filter is used to improve the phase expansion accuracy and measurement accuracy, thereby maximizing the signal-to-noise ratio of the interferogram. The minimum cost flow method is employed to unwrap the phase, and the coherence threshold is set at 0.2. With the aid of precision orbits and DEM, flat ground, phase, and noise are removed, and a high-coherence region unwrapping map is obtained. The coherence coefficient method is used for the extraction of high-coherence points, and 62 high-quality interference pairs are selected. The linear model equations of linear deformation rate and elevation error are established in the region of high coherence point, and the SVD algorithm is used to solve the model parameters, thereby obtaining the linear deformation phase. After the removal of the linear deformation phase, the residual error contains atmospheric phase, nonlinear deformation, and noise. The residual can be subjected to high-pass filtering in the time domain and low-pass filtering in the spatial domain according to the characteristics of each component, thereby removing the atmospheric phase components and noise to obtain a nonlinear deformation phase. Finally, the least square method is utilized to obtain the deformation rate and cumulative deformation value of each coherent point.Figure 5
SBAS-InSAR time-space baseline diagram.


### 5.2.2. Result Analysis

Figure6 shows the annual average settlement rate of the study area. It can be seen from the figure that the center of the subsidence basin has been seriously out of coherence, resulting in the lack of most of the settlement information in the center of the working face, and only the marginal settlement with a smaller magnitude is still preserved better. The settlement magnitude ranges from -0.105 m to 0.035 m, which is relatively small. The blue line is the two cross-sections of the direction and the inclination. The ground GPS observation points are arranged along the direction and the inclination, respectively, which are the green points in the figure. The direction point number is Z1-Z37, and the inclination point number is Q1-Q40. There are labels in the figure. In order to further study the time series settlement trend of the boundary points, the time series diagram obtained by SBAS-InSAR is derived as shown in Figure 7.Figure 6
The annual average settlement rate of the study area.Figure 7
SBAS-InSAR monitoring time series deformation diagram.As shown in Figure7, the black area refers to the working surface, the blue line represents the two section lines along the strike and the inclination, and the black mark points around the working surface are the extracted boundary points. The boundary points of B12 and B14 coincide with the points of Z2 and Q3 on the oblique observation line in the strike and inclination direction, which can be used to detect the observation accuracy of boundary points. During the study period, the 112201 working face was being mined, so the settlement area was mainly concentrated in the 112201 working face. It can be seen from the time sequence diagram that due to the unique geological conditions of the study area, the surface deformation due to coal mining in the Loess Plateau area is characterized by large gradients and large scale, which causes the spatial decorrelation of the SAR interferogram. In that case, high-coherence points are rarely selected above the working surface, which results in the absence of goaf center in the InSAR topographic map. However, based on the basin effect of the goaf surface subsidence and the law of surface deformation of mining subsidence, it can be found that the mining center of the mining area has the characteristics of large deformation speed. However, the basin boundary is a slow linear deformation point, which is in line with the deformation point that can be detected by InSAR. Therefore, only the time series settlement of the boundary point is realized in the timing diagram, and the mining center is not reflected in the InSAR results. Therefore, in this paper, we only discuss the deformation law of the boundary point of the subsidence basin. It can be seen that the settlement funnel extends from southeast to northwest, which is consistent with the mining direction. To verify the monitoring accuracy of SBAS-InSAR, B12 and B14 are utilized to verify their accuracy. Due to the fact that the deformation in the LOS direction is obtained by InSAR, while the GPS observation obtains the vertical and horizontal displacement, firstly, the displacements in three directions, i.e., the vertical direction, east-west, and north-south, obtained by the GPS observation are projected to the LOS upward, and the projection formula is shown in equation (15).As shown in Figure8, the GPS measurement results are projected to the LOS backward and compared with SBAS-InSAR. The two GPS observations at points B12 and B14 in the early stage of mining are two. Comparing the results under the conditions that the time and space references are unified, it can be seen that the GPS measured values at points B12 and B14 are slightly larger than the InSAR monitoring value. And the maximum residual error is less than 5 mm, while the maximum residual error at point Q3 is 7 mm. It can be seen that the overall fit between the InSAR results and the GPS results is good, the correlation is high, and the residual line fluctuates in a small range near the 0 value; besides, the fluctuations are stable. It shows that the results of SBAS-InSAR at boundary points are trustworthy, which can be used to invert model parameters in mining subsidence.Figure 8
GPS/InSAR comparison results of boundary points.


### 5.2.3. Analysis of GWO Algorithm Parameter Inversion Results

Since the mining in 112201 working face was started in September 2018 and ended in November 2019, the subsidence basin was basically stable till the last scene selected in this paper during the imaging in January 2021. Therefore, the cumulative settlement in the last scene of the above SBAS-InSAR is used as the LOS input value of the fitness function, and the error fitting diagram of the fitness function is shown in Figure9.Figure 9
LOS fitting residual plot of fitness function.Due to the instability of the algorithm, the parameter value will be affected by the gross error. To mitigate the influence of these errors and the instability of the algorithm, we especially solve the parameters for 5 times and find the average value. The result of the solution is shown in Table2.Table 2
GWO expected parameter inversion results.
Serial numberqbtanβθ0S1 (m)S2 (m)S3 (m)S4 (m)10.620.352.3685.12-39.0624.78-46.6241.2520.580.392.9383.19-35.1526.16-44.8937.5430.650.373.0583.24-34.9428.57-49.1843.1640.670.372.7483.22-37.9926.96-38.5233.7650.640.352.7184.92-34.5827.00-48.9942.39Average0.630.372.7683.94-36.3426.69-45.6439.62As shown in the table above, the prediction results of the probability integral method areq=0.63, b=0.37, tanβ=2.76, θ0=83.94, S1=−36.34 m, S2=26.69 m, S3=−45.64 m, and S4=39.62 m. To further verify the accuracy of the calculated parameters, we substituted the result parameters into the probability integral model, combined with the geological and mining conditions of the working face, to determine the vertical and horizontal displacements of the observation points on the observation line in the strike direction and the inclination direction of the working face. In addition, predictions and analyses of the accuracy of each observation point based on the GPS measured data were made, with the results shown as follows.It is found in Figure10 that the surface movement obtained by the inversion is in good agreement with the measured results. The point with the largest absolute error of sinking is Point Q24 on the oblique observation line, the absolute error of this point reaches 116 mm, and the median error is 63 mm. The point with the largest absolute error of north-south horizontal movement is Point Z15, with the error of 56 mm, and the median error of 23 mm. The point with the largest absolute error of east-west horizontal movement is Point Z34, with the maximum absolute error of 61 mm, and the median error of 29 mm, both of them are within the tolerance range. It illustrates that the inversion of the probabilistic integration method parameters based on the SBAS-InSAR and gray wolf optimization algorithm used in this paper is feasible, which is available to obtain more reliable surface movement probabilistic integration method parameters.Figure 10
Comparison between the predicted result of probability integral and the measured result of GPS.
(a)
Accuracy comparison chart in the vertical direction of the observation point(b)
Accuracy comparison chart for the north and south of the observation point(c)
Accuracy comparison chart in the east-west direction of the observation pointSince the prediction results of a small number of points on the edge of the subsidence basin with smaller settlement magnitudes are not representative, we only consider the points with settlement magnitudes greater than 100 mm in our accuracy evaluation and analyze the accuracy of the results, mainly from the absolute error (AE), relative error (RE), and average absolute percentage error (MAPE) evaluated. And a one-samplet test with a significance of 0.05 and a test value of 60 was done for the sample point error. The specific test results are shown in Table 3.Table 3
Statistical test table of forecast results.
DirectionMAX.AE (m)MAX.RE (%)Mape (%)sigConfidence intervalVertical0.1163590.015[-27,-3]East-west0.0562570[-33,-27]North-south0.0612180[-27,-21]The error is tested and found to be in accordance with the normal distribution. A one-samplet test with a significance of 0.05 is performed. The test value is set to 50. The significance levels in the three directions of the test are 0.015, 0, and 0, indicating rejection null hypothesis. There is a 95% probability that the error value is less than 50, and the confidence interval is estimated. Further, verify the accuracy of the predicted value.


## 5.2.1. Data Processing

The experimental data are mainly based on the data of Sentinel-1A of ESA. The Sentinel satellite is an earth observation satellite launched by ESA in 2014, which is equipped with a C-band synthetic aperture radar, and available for all-weather earth observation throughout the day. The revisit period of a single constellation is 12 days, in which the interference width can reach 250 km, and the ground resolution is (5∗20) m. In this paper, 41 scenes of the ascending orbit were extracted from the archived data between September 3, 2018, and January 8, 2020.The data processing of SBAS-InSAR is mainly divided into two parts, and the first part is data preprocessing. The preprocessing part mainly includes the extraction of single-view complex image (SLC) data, and the selection of super master images for subpixel precision registration of the data. The second step is to remove the error phase. Firstly, set the space baseline and time baseline as 200 m and 33 d, respectively, according to the geological conditions of the study area. After phase unfolding, the maximum spatial baseline is 113.7 meters, and the minimum is 2.8 meters. As shown in Figure5, select 20190113 as the main image. The generation of the interferogram is the most critical step; in this process, the GoIdstein filter is used to improve the phase expansion accuracy and measurement accuracy, thereby maximizing the signal-to-noise ratio of the interferogram. The minimum cost flow method is employed to unwrap the phase, and the coherence threshold is set at 0.2. With the aid of precision orbits and DEM, flat ground, phase, and noise are removed, and a high-coherence region unwrapping map is obtained. The coherence coefficient method is used for the extraction of high-coherence points, and 62 high-quality interference pairs are selected. The linear model equations of linear deformation rate and elevation error are established in the region of high coherence point, and the SVD algorithm is used to solve the model parameters, thereby obtaining the linear deformation phase. After the removal of the linear deformation phase, the residual error contains atmospheric phase, nonlinear deformation, and noise. The residual can be subjected to high-pass filtering in the time domain and low-pass filtering in the spatial domain according to the characteristics of each component, thereby removing the atmospheric phase components and noise to obtain a nonlinear deformation phase. Finally, the least square method is utilized to obtain the deformation rate and cumulative deformation value of each coherent point.Figure 5
SBAS-InSAR time-space baseline diagram.


## 5.2.2. Result Analysis

Figure6 shows the annual average settlement rate of the study area. It can be seen from the figure that the center of the subsidence basin has been seriously out of coherence, resulting in the lack of most of the settlement information in the center of the working face, and only the marginal settlement with a smaller magnitude is still preserved better. The settlement magnitude ranges from -0.105 m to 0.035 m, which is relatively small. The blue line is the two cross-sections of the direction and the inclination. The ground GPS observation points are arranged along the direction and the inclination, respectively, which are the green points in the figure. The direction point number is Z1-Z37, and the inclination point number is Q1-Q40. There are labels in the figure. In order to further study the time series settlement trend of the boundary points, the time series diagram obtained by SBAS-InSAR is derived as shown in Figure 7.Figure 6
The annual average settlement rate of the study area.Figure 7
SBAS-InSAR monitoring time series deformation diagram.As shown in Figure7, the black area refers to the working surface, the blue line represents the two section lines along the strike and the inclination, and the black mark points around the working surface are the extracted boundary points. The boundary points of B12 and B14 coincide with the points of Z2 and Q3 on the oblique observation line in the strike and inclination direction, which can be used to detect the observation accuracy of boundary points. During the study period, the 112201 working face was being mined, so the settlement area was mainly concentrated in the 112201 working face. It can be seen from the time sequence diagram that due to the unique geological conditions of the study area, the surface deformation due to coal mining in the Loess Plateau area is characterized by large gradients and large scale, which causes the spatial decorrelation of the SAR interferogram. In that case, high-coherence points are rarely selected above the working surface, which results in the absence of goaf center in the InSAR topographic map. However, based on the basin effect of the goaf surface subsidence and the law of surface deformation of mining subsidence, it can be found that the mining center of the mining area has the characteristics of large deformation speed. However, the basin boundary is a slow linear deformation point, which is in line with the deformation point that can be detected by InSAR. Therefore, only the time series settlement of the boundary point is realized in the timing diagram, and the mining center is not reflected in the InSAR results. Therefore, in this paper, we only discuss the deformation law of the boundary point of the subsidence basin. It can be seen that the settlement funnel extends from southeast to northwest, which is consistent with the mining direction. To verify the monitoring accuracy of SBAS-InSAR, B12 and B14 are utilized to verify their accuracy. Due to the fact that the deformation in the LOS direction is obtained by InSAR, while the GPS observation obtains the vertical and horizontal displacement, firstly, the displacements in three directions, i.e., the vertical direction, east-west, and north-south, obtained by the GPS observation are projected to the LOS upward, and the projection formula is shown in equation (15).As shown in Figure8, the GPS measurement results are projected to the LOS backward and compared with SBAS-InSAR. The two GPS observations at points B12 and B14 in the early stage of mining are two. Comparing the results under the conditions that the time and space references are unified, it can be seen that the GPS measured values at points B12 and B14 are slightly larger than the InSAR monitoring value. And the maximum residual error is less than 5 mm, while the maximum residual error at point Q3 is 7 mm. It can be seen that the overall fit between the InSAR results and the GPS results is good, the correlation is high, and the residual line fluctuates in a small range near the 0 value; besides, the fluctuations are stable. It shows that the results of SBAS-InSAR at boundary points are trustworthy, which can be used to invert model parameters in mining subsidence.Figure 8
GPS/InSAR comparison results of boundary points.


## 5.2.3. Analysis of GWO Algorithm Parameter Inversion Results

Since the mining in 112201 working face was started in September 2018 and ended in November 2019, the subsidence basin was basically stable till the last scene selected in this paper during the imaging in January 2021. Therefore, the cumulative settlement in the last scene of the above SBAS-InSAR is used as the LOS input value of the fitness function, and the error fitting diagram of the fitness function is shown in Figure9.Figure 9
LOS fitting residual plot of fitness function.Due to the instability of the algorithm, the parameter value will be affected by the gross error. To mitigate the influence of these errors and the instability of the algorithm, we especially solve the parameters for 5 times and find the average value. The result of the solution is shown in Table2.Table 2
GWO expected parameter inversion results.
Serial numberqbtanβθ0S1 (m)S2 (m)S3 (m)S4 (m)10.620.352.3685.12-39.0624.78-46.6241.2520.580.392.9383.19-35.1526.16-44.8937.5430.650.373.0583.24-34.9428.57-49.1843.1640.670.372.7483.22-37.9926.96-38.5233.7650.640.352.7184.92-34.5827.00-48.9942.39Average0.630.372.7683.94-36.3426.69-45.6439.62As shown in the table above, the prediction results of the probability integral method areq=0.63, b=0.37, tanβ=2.76, θ0=83.94, S1=−36.34 m, S2=26.69 m, S3=−45.64 m, and S4=39.62 m. To further verify the accuracy of the calculated parameters, we substituted the result parameters into the probability integral model, combined with the geological and mining conditions of the working face, to determine the vertical and horizontal displacements of the observation points on the observation line in the strike direction and the inclination direction of the working face. In addition, predictions and analyses of the accuracy of each observation point based on the GPS measured data were made, with the results shown as follows.It is found in Figure10 that the surface movement obtained by the inversion is in good agreement with the measured results. The point with the largest absolute error of sinking is Point Q24 on the oblique observation line, the absolute error of this point reaches 116 mm, and the median error is 63 mm. The point with the largest absolute error of north-south horizontal movement is Point Z15, with the error of 56 mm, and the median error of 23 mm. The point with the largest absolute error of east-west horizontal movement is Point Z34, with the maximum absolute error of 61 mm, and the median error of 29 mm, both of them are within the tolerance range. It illustrates that the inversion of the probabilistic integration method parameters based on the SBAS-InSAR and gray wolf optimization algorithm used in this paper is feasible, which is available to obtain more reliable surface movement probabilistic integration method parameters.Figure 10
Comparison between the predicted result of probability integral and the measured result of GPS.
(a)
Accuracy comparison chart in the vertical direction of the observation point(b)
Accuracy comparison chart for the north and south of the observation point(c)
Accuracy comparison chart in the east-west direction of the observation pointSince the prediction results of a small number of points on the edge of the subsidence basin with smaller settlement magnitudes are not representative, we only consider the points with settlement magnitudes greater than 100 mm in our accuracy evaluation and analyze the accuracy of the results, mainly from the absolute error (AE), relative error (RE), and average absolute percentage error (MAPE) evaluated. And a one-samplet test with a significance of 0.05 and a test value of 60 was done for the sample point error. The specific test results are shown in Table 3.Table 3
Statistical test table of forecast results.
DirectionMAX.AE (m)MAX.RE (%)Mape (%)sigConfidence intervalVertical0.1163590.015[-27,-3]East-west0.0562570[-33,-27]North-south0.0612180[-27,-21]The error is tested and found to be in accordance with the normal distribution. A one-samplet test with a significance of 0.05 is performed. The test value is set to 50. The significance levels in the three directions of the test are 0.015, 0, and 0, indicating rejection null hypothesis. There is a 95% probability that the error value is less than 50, and the confidence interval is estimated. Further, verify the accuracy of the predicted value.


## 6. Discussion

The prediction of mining subsidence has always been a problem worthy of discussion. The InSAR-PIM prediction algorithm proposed in this paper can be widely used in the field of mining subsidence prediction and has the following advantages. First, this method is based on the most accurate prediction method in the field of mining subsidence, and it has certain reliability. The second method is an improvement on the traditional probability integration method, which is obtained through machine learning methods. The predicted parameters of the probability integral method avoid the cumbersomeness of traditional field observations. Third, this method combines InSAR technology, machine learning algorithms, and probability integration theory. It is a combination of multiple methods and advanced mining subsidence prediction algorithms, which can provide a basis for future research in the field of mining subsidence prediction.This method is based on the SBAS-InSAR technology and the GWO algorithm to obtain the predicted parameters of the probability integral method, which avoids the cumbersomeness of obtaining the predicted parameters of the probability integral method based on setting the rock movement observation station of the working face, but it also has certain limitations. First, this method is restricted by the resolution and wavelength of the SAR data, as well as the unwrapping error and atmospheric delay error in the data processing process. These factors may cause the accuracy of the edge points of the subsidence basin extracted by SBAS-InSAR to be insufficient. This directly affects the accuracy of the later parameters. Secondly, the gray wolf optimization algorithm is not fast enough to converge and the global optimization ability is not strong enough, which may cause the parameter inversion effect to be not good enough. The next step is to improve the optimization ability of the algorithm to make it more accurate. The expected accuracy.


## 7. Conclusion

This paper proposed a method to obtain the predicted parameters of the probability integral method based on the SBAS-InSAR and GWO algorithm, which makes up for the defect of obtaining the predicted parameters of the probability integral method by the traditional method, and the deficiency of InSAR technology in the large gradient deformation of the mining area. The following conclusions are drawn:(1)
The SBAS-InSAR technology is available for the extraction of the stable boundary point sinking information of the mining area, and the verification of the accuracy of the point that coincides with the trend observation line(2)
The maximum error of the strike coincidence point B12 is 4 mm, and that of the tendency coincidence point B14 is 7 mm, indicating that the boundary information of the subsidence basin extracted by this technology is trustworthy(3)
Establish the fitness function based on the geometric relationship between the InSAR lateral imaging and the probability integral model, and obtain the predicted parameters of the probability integral method for the 112201 working face of Xiaobaodang Coal Mine in the northern Shaanxi mining area through the GWO algorithm. Besides, the obtained results include thatq=0.63, b=0.37, tanβ=2.76, θ0=83.94, S1=−36.34 m, S2=26.69 m, S3=−45.64 m, and S4=39.62 m. To verify its accuracy, the calculated parameters were substituted into the parameters of the probability integral model for the prediction of the sinking and horizontal movement of the working face. In addition, the predicted results were compared with the GPS measured data. It showed that the sinking error was 63 mm. The horizontal movement errors of mm, north-south, and east-west are 61 mm and 29 mm, respectively, and the accuracy met the actual production requirements. Therefore, this paper proposed a method for obtaining probability integral prediction parameters that can be applied in actual engineering



---

*Source: 9376711-2021-08-31.xml*