Position Inversion of Goafs in Deep Coal Seams Based on DS-InSAR Data and the Probability Integral Methods

: The goafs caused by coal mining cause great harm to the surface farmland, buildings, and personal safety. The existing monitoring methods cost a lot of workforce and material resources. Therefore, this paper proposes an inversion approach for establishing the locations of underground goafs and the parameters of the probability integral method (PIM), thus integrating distributed scatter interferometric synthetic aperture radar (DS-InSAR) data and the PIM. Firstly, a large amount of surface deformation observation data above the goaf are obtained by DS-InSAR, and the line-of-sight deformation is regarded as the true value. Secondly, according to the obtained surface deformations, the ranges of eight goaf location parameters and three PIM parameters are set. Thirdly, a correlation function between the surface deformation and the underground goaf location is constructed. Finally, a particle swarm optimization algorithm is used to search for the optimal parameters in the range of the set parameters to meet the requirement for minimum error between the surface deformation calculated by PIM and the line-of-sight deformation obtained by DS-InSAR. These optimal parameters are thus regarded as the real values of the position of the underground goaf and the PIM parameters. The simulation results show that the maximum relative error between the position of the goaf and the PIM parameters is 2.11%. Taking the 93,604 working face of the Zhangshuanglou coal mine in the Peibei mining area as the research object and 12 Sentinel-1A images as the data source, the goaf location and PIM parameters of the working face were successfully inverted. The inversion results show that the maximum relative error in the goaf location parameters was 16.61%, and the maximum relative error in the PIM parameters was 26.67%.


Introduction
Coal is currently the main energy source in China, and coal mining has led to the existence of a large number of abandoned goafs. This often results in the reduction of ground elevation, which can cause serious ecological and environmental problems including farmland inundation, soil salinization, and desertification. It can also be detrimental to the safety of people's lives and property [1,2]. Research examining mining subsidence has played a positive role in reducing surface subsidence disasters and economic losses [3][4][5][6]. With the efforts of experts and scholars worldwide, a series of research methods have been established to investigate this problem, including the typical curve method, the profile function method, and the probability integral method (PIM). The latter of these has become a relatively mature method to research the relationship between coal mining and land subsidence [7]. Inversion research on mining subsidence has two main aspects.
Firstly, the PIM parameters are obtained by inversion from the underground mining information and surface subsidence observation data, and the mining subsidence caused by the other working faces can be predicted by the PIM [8]. Secondly, using the known PIM model parameters and surface subsidence observation values, the goaf positions in the mining area can be obtained by inversion [9].
In the process of realizing the inversion of PIM parameters related to goafs or geology, observed values of surface deformation need to be used; the establishment of surface movement observation stations therefore plays an important role in promoting research examining mining subsidence. Leveling, GPS, interferometric synthetic aperture radar (InSAR), and other technologies are widely used in mining subsidence monitoring [10,11]. Due to its large coverage, all-weather capability, and high monitoring accuracy [12], In-SAR is widely used to measure deformations in glaciers [13,14], volcanoes [15,16], earthquakes [17,18], and other landscapes, as well as urban and general land surface deformations caused by resource exploitation. Scholars worldwide have also tried to apply In-SAR technology to the observation of surface deformations caused by mining activities [19,20].
In most existing goaf-location inversion methods, the PIM parameters are given based on experience [21][22][23][24], which leads to large errors. To obtain more accurate model parameters, information should be known about the goaf to invert the model parameters [25][26][27]. However, the existing research still has the following shortcomings: (1) generally, PIM parameters are difficult to obtain accurately; at present, they are mainly given based on experience, which can be arbitrary and may lead to large errors. (2) Even if the PIM parameters are retrieved from surface observation station data, the calculated model parameters are not necessarily reliable due to errors in the monitoring data and the underground working face information. (3) Most previous studies have used differential InSAR (D-InSAR) and leveling and GPS observations from surface movement observation stations to obtain the surface deformation; however, D-InSAR is affected by temporal and spatial incoherence and atmospheric delay, so it is difficult to use this for obtaining accurate deformation information and the point density required for the inversion of underground goaf locations. (4) Most of the existing research cases are related to shallow or horizontal mining, and the influence of the coal seam dip angle is generally ignored; however, when the coal seam is deeply buried, the center coordinate offset is large, and the influence of the dip angle cannot be disregarded.
According to these shortcomings, this paper puts forward the following improvements. Firstly, aiming at shortcomings (1) and (2), both the goaf location parameters and the main PIM parameters are regarded as unknown parameters to participate in the inversion simultaneously; this means that the results are closer to the real values and avoid the large errors caused by the use of empirical values. Secondly, due to the influence of spatial-temporal decoherence, it is difficult to extract high-density surface monitoring points from D-InSAR. This work therefore used the current frontier technology of earth observation: distributed scatter (DS)-InSAR [28]. The core idea of this method is to treat objects (such as bare land, beach, and Gobi) with medium temporal coherence and continuous spatial distribution as homogeneous points. According to the statistical characteristics of these homogeneous points, the phase of distributed objects is optimized and the influence of noise on them is reduced; the density of point objects in non-urban areas is improved, and surface deformation observation values with higher order and accuracy are obtained. Thirdly, in the inversion model, the influence of the coal seam inclination angle and the mining depth on the center coordinate offset are considered. Finally, in view of the fact that there are many inversion parameters, resulting in a high-dimensional nonlinear continuous parameter model, a particle swarm optimization (PSO) algorithm, which has advantages in solving high-dimensional problems, is used for inversion.
The structure of the remainder of this paper is as follows. Section 2 introduces the basic principles and research methods; the results of simulations, and experiments with real data, are presented in Section 3, and a discussion of the inversion errors is given. Conclusions are presented in Section 4.

DS-InSAR Method
Recognition of homogeneous points and phase optimization are the prerequisites for extracting the deformation phases of distributed targets. The common homomorphic point-recognition methods mainly include two types: nonparametric and parametric hypothesis testing methods. Among these, the fast statistically homogeneous pixel selection (FaSHPS) method [29], which is a parametric test, has the advantage of small heterogeneity and a fast recognition speed, and it has been widely used.
According to the characteristics that the single-look SAR amplitude image obeys Rayleigh distribution in homogeneous regions, and their coefficients of variation are constants, the confidence interval based on the amplitude mean can be constructed as described in the FaSHPS method: where P{·} is the probability, μ(L) is the expectation value of pixel L, Z1−α/2 is the (1 − α/2) quantile in the standard normal distribution, N is the number of SAR images in the time series, and

 
A L is the mean value of SAR image amplitude. According to the above formula, we only need to calculate the average amplitude of a SAR image pixel in a time series and then establish whether the value falls into the above confidence interval to judge whether it belongs to the set of homogeneous points.
After the identification of homogeneous points, our approach uses the coherent matrix decomposition method [30] to optimize the phase. In this method, the eigenvalues of the pixel covariance matrix of time-series SAR data are used to separate different scattering signals. The eigenvector corresponding to the largest eigenvalue is regarded as the main scatterer, and the eigenvectors corresponding to the other eigenvalues are regarded as noise. The phase component of the dominant scattering mechanism is separated to optimize the phase. The model is expressed as where signal T  and noise T  are the main scatterer signal and the noise signal, respectively; the eigenvector μ1, corresponding to the maximum eigenvalue λ1, is the optimized phase of the distributed target; λi are the n eigenvalues of the coherence matrix; and μi is the eigenvector corresponding to the different eigenvalues λi and thus corresponding to different scattering mechanisms. The processing flow of DS-InSAR is as follows: (1) the interferometric pairs are set up according to the small spatio-temporal baseline, and the differential interferograms are generated by removing the flat and terrain phase based on satellite orbit parameters and external DEM. (2) The homogeneity points are identified by the amplitude images, and the phases of the differential interferogram are optimized. Then, the spatio-temporal coherences of the homogeneous points are calculated by the optimized phase, and the DS points are obtained by setting the threshold.

Probability Integral Model
As shown in Figure 1, in the coal seam coordinate system, sO1t corresponds to the surface coordinate system xOy; in the surface coordinate system, xOy is the projection of sO1t in the coal seam coordinate system. If a unit B(s, t) is exploited at (s, t), the land subsidence value of A(x, y) at any point on the surface caused by this micro unit is: If the strike length of the mining working face is D1 and the dip length is D2, the land subsidence value of point A(x, y) caused by coal mining of the whole working face is: where m is the average mining thickness of the working face; q is the subsidence coefficient; α is the inclination angle; tanβ is the tangent of the main influence angle; H, H1, and H2 are the mining depths of the strike main section, the lower boundary of the working face, and the upper boundary of the working face, respectively; l is the calculated length of the working face strike, l = D1 − S1 − S3; and L is the calculated length of the working face dip, L = (D2 − S2 − S4) sin(θ + α)/sinθ, in which S1, S2, S3, and S4 are the offset distances in the left, down, right, and upper directions and θ is the propagation angle of the mining influence.
The horizontal movement U(x, y, φ) of any point A(x, y) on the surface in the φ direction (an angle value that rotates counterclockwise in the positive direction of the x axis to the specified direction), can be expressed as: where r is the main influence radius and b is the coefficient of horizontal movement. According to Equation (8), when φ is an angle value in the north-south direction, the calculated horizontal movement is UN. When φ is an angle value in the east-west direction, the calculated horizontal movement is UE. It can be seen from Figure 2 that the surface deformation caused by the working face is biased to the inclination direction of the coal seam, and the offset value is where k1 and k2 are proportional constants.

Establishment of Inversion Model
The LOS surface deformation obtained by DS-InSAR is recorded as DLOS; the LOS deformation is calculated by PIM according to the relationship between the three-dimensional surface deformation W, UN, and UE, and this is recorded as Dlos. The model for this is where θ is the incident angle of the SAR satellite radar beam and αh is the satellite heading angle.
The DLOS obtained by DS-InSAR is taken as the true value, and the Dlos calculated by PSO and PIM each time is taken as the observed value. When the error is less than a given threshold, the calculated goaf location and PIM parameters are taken as the final results:

Numerical Simulations
In this work, the simulated mining depth was 900 m, the length was 1000 m, the width was 500 m, the mining thickness was 6 m, and the dip angle was 25°. In the simulations, the scale constants were k1 = 0.1 and k2 = 0.6, and the horizontal resolution of the model block was 20 m × 20 m. Figure 3 shows the surface subsidence and LOS deformations caused by this simulated coal mining. It can be seen that the maximum settlement point is in the center of the basin, and the maximum settlement value is 2.8 m. Due to the influence of the inclination angle and the coal seam burial depth, the center coordinate of the working face deviates from the settlement center.   Figure 3. It can be seen from Figure 4a that the horizontal movement between the boundary of the subsidence basin and its inflection point increases gradually; it is the largest at the inflection point of the subsidence curve, and it then decreases gradually from the inflection point to the point of greatest subsidence, where the horizontal movement is zero; the second half of the curve is a symmetrical reflection of the first. The difference between Figure 4a,b is that the horizontal movement of the coal seam in the downhill direction is greater than that in the uphill direction, and the curve is asymmetric. It can be seen from Figure 4 that the surface burial depth is large, and it is not easy to realize full mining after mining.  Based on the simulated underground goaf data in the previous section and the model parameters determined by geological conditions, the subsidence W of each point on the surface was generated by the PIM, providing the horizontal movement UN in the northsouth direction and the horizontal movement UE in the east-west direction. Simulated InSAR surface deformation observation values were generated using Equation (12) and recorded as DLOS, and taking the underground goaf data and PIM parameters as unknown values, an inversion experiment was carried out. The inversion results were compared with the simulated values, and the inversion values of each parameter were obtained. The absolute and relative errors are shown in Table 1. The following can be seen from Table 1. (1) Among the length parameters (D1, D2, Xc, Yc, H, m), the maximum relative error is 1.21% in the coal seam thickness m, the absolute error is 0.07 m, and the maximum absolute error is 6.5 m in the coal seam depth H. (2) The maximum absolute error and the maximum relative error of the two angle parameters appear in the inclination angle α, and the maximum absolute error and the maximum relative error are 0.53° and 2.11%, showing that the two angle parameters have a good inversion effect and can better approximate the true value. (3) For the PIM parameters (q, b, tanβ), the maximum relative error and absolute error appear at the tangent of the major influence angle tanβ; in this example, the absolute error is 0.01 and the relative error is 0.45%, which is similar to the angle parameters and can better approximate the real parameters. (4) The parameters (Xc and Yc, D1 and D2) have the same levels of error under the conditions of using simulated data for inversion. At the same time, the errors of Xc and D2 are slightly larger than those of Yc and D1, which is caused by the complexity of tendency calculation model being greater than that of the trend, which is more consistent with the actual situation.
These simulations show that the proposed method is feasible. For areas with unknown PIM parameters and underground goafs, the time-series surface deformation field is first obtained by InSAR imaging, and the inversion of the PIM parameters and the locations of goafs can then be obtained simultaneously.

Study Area
The surface of the Zhangshuanglou coal mine in the Peibei mining area ( Figure 5) belongs to the Yellow River alluvial plain, which has ground elevations of 37-39 m. The terrain is high in the west and low in the east, and the surface water system is dense. Longterm and high-intensity underground coal mining has resulted in a large area of surface subsidence and even collapse, which not only affects the ecological environment of the area but also endangers the lives and property of residents. The coal seam in this area is stable and of medium thickness, with an average total thickness of 5.82 m, and the surface of the region is mostly covered by vegetation. It is difficult to obtain high-precision and long time-series surface deformations using the D-InSAR method. The shape of the 93,604 working face from April 2017 to July 2018 is shown in Figure 5. The coal seam in the north of 93,604 has not been mined, and the geology there is stable. The southern coal seam was mined in July 2016, and the surface deformation has been stable; this will therefore not affect the deformation monitoring of the 93,604 working face.

SAR Data Processing
Twelve repeatable-track sentinel images were selected for the experiment, and the pixel sizes in the distance and azimuth directions were 2.33 m and 13.91 m, respectively.
The parameters of the observation time, vertical baseline, and time baseline of the images are shown in Table 2. To reduce the computation time, the images were first cropped to 2400 × 500 pixels. The surface deformations obtained by the accumulated phase based on DS-InSAR are shown in Figure 5. In this experiment, the mask threshold was 0.4, and the multi-view ratio was 4:1 in the range and azimuth directions. A total of 78,605 high-coherence points were selected. The maximum surface deformation value generated by DS-InSAR was about 323 mm.
The locations of the mining subsidence basins can be judged according to the monitoring results in Figure 6. It can be seen that there are three areas with obvious surface deformation. In this experiment, it is necessary to verify the accuracy of three PIM parameters and eight goaf location parameters. Therefore, two groups of experiments were designed. (1) The cyan area was used to obtain accurate PIM parameters [8]; the aim of this group of experiments was to provide the PIM parameters of this area for the second experiment to verify the experimental results. (2) The red was used to obtain the positioning parameters and PIM parameters of the underground goaf to verify the feasibility of the proposed method. To reduce the effect of decoherence, SAR data obtained in spring and winter was used to ensure that the study area would not be covered by large-scale vegetation.

Inversion of Goaf Parameters
The real values of the goaf parameters (D1, D2, Xc, Yc, H, φn, α, m) can be accurately obtained from the mining map, and the errors in these data can be ignored. The parameters of the PIM model (q, b, tanβ) related to the geological mining conditions need to be obtained by geological survey, and they are often measured as the mean value of parameters in a large range of the region. Although this value is appropriate for general practical applications, it is not good enough for small-scale geological inversions due to the complexity of the geological environment. Therefore, in this work, we used the area indicated by the cyan frame in Figure 6 to obtain the PIM parameters of the area using accurate goaf parameters and surface observations of deformation obtained by InSAR. The inversion method is shown in [8], and the inversion result was q = 0.83, b = 0.3, tanβ = 1.8.
To evaluate the accuracy of this experiment, the absolute and relative errors between the inversion parameters and the real values were calculated, and the calculation results are shown in Table 3. The following can be observed from Table 3. (1) In the inversion results of the goaf parameters (D1, D2, Xc, Yc, H, φn, α, m), the maximum relative error is 16.61%, and the absolute error is 26.57 m, which appears on the dip length D2. (2) Despite being the same type of parameters, the absolute error and relative error of the strike and dip length D1 and D2 are quite different. The main reason for this is that in PIM, the dip length needs to be projected onto the horizontal plane to calculate the horizontal distance, and in this process the parameter is also affected by the parameters α, H, while the strike length is only affected by the parameter H. (3) The absolute and relative errors of the angle parameter φn and α are small. This is because these two parameters are highly correlated with the surface deformation morphology, and the surface deformation morphology obtained by DS-InSAR technology can easily determine the approximate value range. (4) The inverted position of the center of the inclined coal seam and the deformation center observed by DS-InSAR deviates by ∆Y in the dip direction, as shown in Equation (9). Due to the offset considered in PIM, the inversion parameters Xc and Yc can approach the real values with better accuracy. Figure 7 shows a projection of the real LOS deformation (upper layer) and the actual position of the mining working face on the horizontal plane (lower layer). The limited differences between the two rectangles indicates that the inverted result is similar to the real working face in terms of shape; the position deviation is small, and this demonstrates a good inversion effect.

Goaf Inversion with Known PIM Parameters
Most previous studies have used known PIM parameters to invert the goaf position. In this section, we make a comparison of this approach with the method of this paper. To ensure that the experimental data are not affected by other errors, the LOS deformation observed by InSAR was generated based on the simulated data in Table 1 and Equation (12). Since the LOS deformation has no errors, when the PIM parameters are assumed to be known, the inverted goaf position parameters should be consistent with the simulated values. Based on this assumption, the inversion results were compared with the simulation results in this paper to explore the interaction between the PIM parameters and the goaf position. The inversion results are shown in Table 4.  Table 1, but this improvement is not obvious. (2) For the parameter φn, the two inversion methods have the same effect and can obtain an accurate final result, which shows that this parameter has strong robustness. (3) Comparing the results of the two groups of simulation data, we can see that whether the PIM parameters (q, b, tanβ) are known or not has a little effect on the inversion accuracy of the goaf position parameters; the relative and absolute errors decrease slightly, which shows that the research method in this paper is feasible.

Influence of the LOS Deformation Obtained by DS-InSAR
There are many factors influencing DS-InSAR observations, including DEM errors, orbit errors, atmospheric delay, phase unwrapping errors, and noise. The final deformations generated by DS-InSAR may thus contain observation errors. To theoretically analyze the influence of LOS observation errors on the results of parameter inversion, deep coal seam data was used in simulations to verify the method in this paper. To observe changes in the parameters, random errors from −0.4 m to 0.4 m with an average interval of 10 mm were added to LOS observations without errors. Scatter plots of the variations in the 11 parameters with different error levels are shown in Figure 8. The following can be seen from Figure 8. (1) The parameters D1, D2, Xc, Yc, φn, α, q, b, and tanβ change linearly with increasing LOS error; the inversion errors of parameters H and m do not change synchronously with the LOS error, but rather they fluctuate around the zero axis. (2) The LOS deformation errors directly affect the inversion parameters, so the LOS observations used should be highly accurate. This is why DS-InSAR was used here to provide high-quality observations.

Conclusions
This paper presents a method to simultaneously retrieve the parameters of the position and PIM parameters of underground goafs using DS-InSAR data. The method improves the target density of points in non-urban areas by reducing the influence of noise on distributed targets using DS-InSAR. It then establishes a correlation model between underground mining and surface deformation. Finally, a PSO algorithm is used to quickly optimize the parameters. The following conclusions were reached.
(1) The amount of surface deformation data used is very important to constrain the inversion model. Therefore, DS-InSAR is applied to reduce the influence of spatio-temporal incoherence, and this can effectively increase the density and accuracy of the fieldobserved deformations.
(2) For mining subsidence caused by exploitation of deep coal seams, the center of the surface deformation deviates from the position of the underground goaf due to the coal seam dip angle. Thus, the center deviation caused by inclined coal seams in the model must be considered, or the inversion error in the center coordinate will be larger.
(3) Simulation results show that it is feasible to use the PIM parameters as unknown parameters in goaf inversion, and the inversion errors are relatively small. Experiments with real data verified the results of the simulations. In the case of the PIM parameters being involved in the inversion, the goaf position parameters can still be obtained with a high accuracy. Because the PIM parameters are difficult to obtain accurately, the method in this paper avoids the need for their selection according to experience in goaf location inversion.
(4) The maximum relative error of the simulations was 2.11%, the maximum relative error in the experiments with real data was 26.67%, and the errors in other inversion parameters were relatively small. The experimental results show that the method has a good effect on the inversion of underground goaf locations and has the advantage of a large range of non-contact measurements.