Goaf Locating Based on InSAR and Probability Integration Method

: Mining goafs can cause many hazards, such as burst water, spontaneous combustion of coal seams, surface collapse, etc. In this paper, a feature-points-based method for the efﬁcient location of mining goafs is proposed. Different interferometric synthetic aperture radar (DInSAR) is used to monitor the subsidence basin caused by mining. Using the principles of the probability integral method (PIM), the inﬂection points and the boundary points of the basin monitored by DInSAR are determined and used as feature points to locate the goaf. In this paper, the necessity of locating goafs and the traditional methods used for this task are discussed ﬁrst. Then, the results of verifying the proposed method by both a simulation experiment and real data experiment are presented. Six RADARSAT-2 images from 13th October 2015 to 5th March 2016 were used to acquire the subsidence basin caused by the 15235 working faces of the Jiulong mining area. The average relative errors of the simulation experiment and real data experiment were about 6.43% and 12.59%, respectively. The average absolute errors of the simulation experiment and real data experiment were about 28 m and 38 m, respectively. In the ﬁnal part of this paper, the error sources are discussed to illustrate the factors that can affect the location result.


Introduction
Coal production and coal consumption are extensive in China, as coal is the country's main energy source. According to the National Energy Administration (http://www.nea.gov.cn), coal production capacity in China is increasing every year. As of the end of June 2017, the national total coal production capacity was 3.41 billion tons/year, and it had increased to 3.49 billion tons/year by June 2018. It is estimated that by 2020, the national energy consumption will still include more than 60% coal. In 2017, the world's coal production was 7.73 billion tons, of which China accounted for 45.6%.
Such large-scale mining results in many goafs. Many of them are recorded, but because of the existence of illegal mining and loss of records, there are many unknown goafs hidden underground with the potential to cause many disasters. Mining under the original goaf activates the overlying rock, which can then lead to secondary disasters, such as ground collapse [1,2]. Mining activities next to goafs can lead to falling roofs, aquifer rupture, and other hazards that can endanger the safety of underground workers [3]. The sudden collapse of a goaf will affect agricultural production, as well as result in damage to the buildings above, landslides, and other geological disasters that endanger

DInSAR Technology
DInSAR uses the phase information of SAR complex images to identify land surface deformation. An interferogram can be obtained by multiplying two SAR images by a conjugate. In the interference fringe pattern of two images, the interferometric phase ϕ is mainly composed of the following parts: where ω(·) is the wrap operator; ϕ top is the phase caused by the terrain, which can be removed by external digital elevation model (DEM) simulation; ϕ f lt is the flat phase, which can be removed by a precise baseline estimate; ϕ de f is the phase caused by the surface deformation in the line-of-sight (LOS) direction; ϕ atm is the phase caused by atmospheric delay. By stacking the interferograms, the signal-to-noise ratio (SNR) of the deformation information can be improved, and the noise of atmospheric interference can be weakened, so the influence of the atmospheric delay can be reduced [23]. The parameter ϕ noise is the noise phase and can be suppressed by low-pass filtering; then, the real phase is restored by the minimum cost flow method to obtain the surface deformation [24].

Feature Points of PIM
The PIM was first proposed by Polish scholars. It considers the movement of particulates, such as sand or tiny rock masses, as a random event, which makes the relationship between particulates irrelevant. The mining area can be divided into a number of differential units with the assumption that it is mined. The total impact of complete mining on the rock formation and the surface can be regarded as the integral of the small impact generated during the mining of the differential unit [25].
As shown in Figure 1a, when the ball (named a 1 ) on the first floor is taken away, one of the balls on the second floor (a 2 and b 2 ) will fall, so the probability of falling for each of the two balls is 1/2. Then, one of the balls on the third floor will fall, and the falling probabilities for each of the three are 1/4, 2/4, and 1/4, respectively (in Figure 1b). On this basis, the differential equation of the probability distribution function of the sinking event is obtained. For a clear expression, the plane, which is beyond the center of the subsidence basin and the vertical with horizontal, is appointed as the main section. The main section in the strike direction is called the strike main section, and the one in the dip section is called the dip main section. After assuming that the medium is infinitely small, in the main section, the surface subsidence profile equation of semi-infinite mining is where W(x) stands for the subsidence value at x; W max is the max value of the subsidence curve; r is the major influencing radius; er f is the error function, and er f = 2 √ π x 0 e −η 2 dη. The location of x = 0 is the coal mining boundary.
According to differential geometry, the second derivative of subsidence W(x) is the curvature: where K(x) is the curvature at x. In the deformation basin, the value of ( dW(x) d(x) ) 2 is lower by more than three orders of magnitude. When it is neglected, the formula for calculating the curvature can be simplified to According to Equation (4), when x < 0, K(x) > 0, the curve of W(x) is upward convex; when x > 0, K(x) < 0, the curve of W(x) is downward concave. Thus, the point at x = 0 is the turning point between concave and convex geometries, and this is called the inflection point of the subsidence curve. That is to say, when the deformation of the overlying strata fully conforms to the law of a stochastic medium, the projection of the inflection point of the subsidence curve and the mining boundary point on the horizontal plane coincide. So, the inflection point can be used as a feature point to determine the goaf boundary.   The mining depth can be calculated by the inflection points and boundary points. As can be seen in Figure 2, H 0 , δ 0 , S , and S satisfy the following trigonometric function relation: with δ 0 and S obtained from [26]. S can be obtained from the inflection point and the subsidence basin boundary point, so H 0 can be obtained by Equation (5). The four boundaries obtained by DInSAR can be used to calculate the mining depth. Subsidence and curvature curves. Curve (1) represents the distribution law of surface subsidence W(x) caused by mining activities. Curve (2) represents the distribution law of the curvature K(x). It is also the second derivative of the subsidence curve.

Goaf Locating
Using the feature points, the boundaries and depth can be obtained. However, when calculating the boundaries, the actual situation is that the inflection point is located above the mining boundary and slightly deviates to the inside of the goaf. There are two reasons for the phenomenon: first, during mining activities, the roof caves in and forms a voussoir beam structure, which causes the roof near the mining boundary to bend stepwise or not fall, and there is a suspended roof distance; second, the fracture angle significantly affects the inflection point position. The actual surface subsidence and curvature curve caused by underground coal mining are shown in Figure 2.
In Figure 2, O is the basin center, H 0 is the mining depth, E is the inflection point of the subsidence curve, and E is the zero point of the curvature curve. S is the deviation of the inflection point, that is, the distance between the projection of the inflection point on the coal seam plane and the mining boundary. S is the projection of the distance from the boundary point of the basin to the zero point of the curvature curve E on the horizontal plane. The parameter δ 0 is the boundary angle in the main section, that is, the angle between the horizontal line and the line from the basin boundary point to the mining boundary point on the pillar side. As can be seen, the surface subsidence value is at its maximum above the center of the goaf and gradually decreases from the center of the basin to the edge, and it becomes zero at the basin boundary (the position at 10 mm subsidence is considered the boundary point during field measurements). The subsidence curve is symmetrical at the center. The curvature value is zero at the edge of the basin, it first increases and then decreases in the basin center's direction and is zero again at the basin center. S is affected by many factors, such as the overlying rock properties and the existence of an adjacent goaf and geological structures. However, the deviation of the inflection point is usually considered a constant value in the same region in practical applications, because the geological conditions and strata attributes of the same area do not change much [26]. Obviously, this introduces some errors, but they are easily calculated.
After obtaining the subsidence basin by DInSAR, the inflection point E and deviation of the inflection point S can be used as feature points to determine the goaf boundary. For a flat coal seam, the two boundaries in the dip direction are symmetrical with respect to the center of the subsidence basin, so only one inflection point on either side of the dip direction is needed to obtain the two boundaries. The situation in the strike direction is more complicated. When monitoring the surface deformation by InSAR, the mining face is being mined, and because it takes time for the deformation to propagate in the overlying strata, the development of the surface deformation is delayed. This situation is described in Figure 3. In Figure 3, (1/4 ∼ 1/2)H 0 is the excavation distance of the working face when the ground surface begins to move. M1 ∼ M5 are the excavation positions of the working face from different periods, W1 ∼ W5 are the sink curves corresponding to M1 ∼ M5 when the ground surface movement is over, and B1 ∼ B5 are the boundary points of W1 ∼ W5. W5 is the sink curve when the excavation point has just reached M5 while the ground surface is still moving. B5 is the boundary point of W5 . δ 0 is the boundary angle in the main section, and M5 is the excavation point according to δ 0 . M5 is the calculation result for the mining boundary if the basin W5 is used, but it takes months for the basin to expand from B5 to B5. Therefore, in the locating process, according to geological conditions and prior experience, we move the characteristic points (inflection points and boundary points) to a certain distance in the mining direction to counteract the error caused by the surface deformation delay. When determining the depth, the excavating boundary may lead to a large error; the extraction size in the dip direction may not reach critical situation (the reason will be introduced in Section 5.3), so using the back basin boundary in the mining direction to calculate the mining depth will get better results. The sketch of the goaf locating can be represented by Figure 4.

Acquiring subsidence basin by DInSAR
Acquiring feature points

Goaf boundary determination
Calculating the distance between feature points Goaf location Goaf depth determination

Simulation Experiment
In order to verify the relationship between the characteristic points of the subsidence basin and the goaf position, the fast lagrangian analysis of continua (FLAC3D) was used to simulate mining activity, and the goaf position was then calculated according to the subsidence basin. FLAC3D adopts an explicit Lagrangian algorithm and a hybrid-discrete partition technique, which can simulate the plastic break and flow simulation parameters very accurately.

Numerical Simulation
In this simulation experiment, a horizontal working face was simulated with a mining depth of 200 m, a length of 1000 m, a width of 300 m, and a mining thickness of 7 m. The larger mining thickness and smaller mining depth ensured that the subsidence basin was fairly typical. The ratio of depth to thickness was close to 30, which ensured the continuity of the surface deformation. The horizontal resolution of the model block was 20 m × 20 m, which ensured the accuracy of the simulation while accounting for operational efficiency. The bottom, left, and right boundaries were constrained by horizontal displacement. The Mohr-Coulomb model was used to simulate the movement of the overlying rocks and surface subsidence. The Mohr-Coulomb model is a general model that is used to simulate mining activities [27]. It is a model applied to plastic materials and the calculation of deformation and stress according to the Mohr-Coulomb failure criterion (Equation (6)) and the maximum stress criterion.
As can be seen in Figure 5, the maximum subsidence point was in the center of the basin, and the maximum subsidence value was 4.76 m. The tilt and curvature were obtained by taking the first-order and two-order derivative of the subsidence basin (shown in Figure 6). To facilitate the display, the acquired tilt and curvature were processed by mirroring. The value represented by the tilting and curvature map was opposite to the true curvature and the tilt value. As can be seen in Figure 6, the wave crest and wave trough were opposite to the main section in the dip direction and the strike direction in graphs (1) and (2), respectively. In graph (1), there is a distance between the edges of the hump and the pit, while there is nearly no such distance in graph (2), which indicates that the mining degree in the strike direction is larger than that in the dip direction. In graphs (3) and (4), the curvature was negative at the edge and positive in the center. For the positive areas, the curvature value was large at the edge and small in the center in graph (3); this is contrary to that in graph (4), which indicates that the mining degree was full in the strike direction and subcritical in the dip direction. In order to more clearly show the variation in the subsidence and curvature curves, we obtained the subsidence curve and the curvature curve on the strike main section and the dip main section. The result is shown in Figure 7. In Figure 7, graph (a) shows that the bottom of the subsidence curve was gentle. The value of the curvature curve decreased after first increasing, and it then continued to increase until it was zero from the edge to the center and zero at the inflection point, indicating that the mining degree in the strike direction was full. In graph (b), the bottom of the subsidence curve is sharp. The value of the curvature curve decreased after increasing; it then fluctuated below the axis and was zero at the inflection point, indicating that the mining degree was subcritical in the dip direction.  (3) and (4).

Simulated Goaf Locating
Equation (3) shows that the inflection point was close to the goaf boundary after projection on the horizontal plane. If the red lines (a, b, c, d in Figure 6) are moved by the distance of the deviation of the inflection point toward the side of the coal pillar, they become the goaf boundary. Therefore, according to Equation (4) and the position of the inflection point and boundary point in the main section, the mining depth can be determined. Because the simulation data were ideal (as shown in Figure 5), there was no need to fit the subsidence curve.
According to [26], the boundary angle δ 0 is 51 • , and the deviation of the inflection point S was 20 m. Because the basin was stable, there was no need to consider the delay distance. The result is shown in Figure 8. As can be seen in Figure 8, the length of the calculated goaf was smaller than the real one, but the calculated width and depth are larger. According to Equation (5), the depth is directly proportional to the distance from the mining boundary to the basin boundary: thus, the larger the length, the smaller the depth. This conclusion agrees with the experimental result.

Accuracy Evaluation
In order to accurately evaluate the error of this method, the calculated goaf location was compared with the position of the goaf in the simulation model. The error of the simulation experiment is shown in Table 2. Table 2 shows that the length of the calculated goaf was smaller than the real length, and the relative error was 5.00%. Because there was no need to consider the delay distance and the mining degree in the strike direction was full, the error of the goaf boundary in the strike direction can be assumed to be mainly caused by the deviation of the inflection point and simulated block resolution. The crushing and deformation during mining caused the plasticity of the overlying strata to decrease, which led to a decrease in the length of the voussoir beam structure and a slight change in the deformation propagating in the overlying strata. The empirical value of the deviation inflection point was obtained with the measured data from monitoring the mining subsidence deformation, so the change in the overlying strata properties were reflected in the empirical value. However, once the mechanical parameters of the overlying strata were set in the simulation experiment, it was difficult to modify them during the simulation calculation process, and this resulted in a larger inflection point deviation. Since the inflection point deviation was not reduced properly in the goaf boundary calculation process, the calculated strike boundary was shifted to the basin center because of the parameter change. The resolution of the ground observation points was 20 m, and the inflection point can only be located at one point and cannot be between two points. So, the maximum error caused by the resolution could have reached 10 m, and the direction was uncertain. The mining degree in the dip direction was subcritical, which causes the inflection point to shift away from the basin center, and this error was in the opposite direction of the parameter change. Compared with the real goaf boundary, the calculated boundary shifted away from the basin center, indicating that the error caused by subcritical mining was greater than that caused by the parameter change. The mining depth was calculated by the boundary behind the mining direction and inflection point, so it was not affected by subcritical mining and delay distance. The boundary angle, inflection point deviation, strike boundary position, and observation point resolution were the main error sources for calculating the goaf depth. Similar to the error source of the inflection point deviation, the property of the overlying strata changed with deformation, but the mechanical parameters in the numerical simulation model were always constant; this led to a smaller boundary angle in the simulation experiment than the empirical value used in the real experiment, resulting in a smaller calculated value of the mining depth. In general, the combined influence of the boundary angle, inflection point deviation, and basin observation point resolution results in a calculated mining depth that is larger than the actual value. Compared with the real data experiment, the simulation experiment was simple and not affected by complex geological conditions. The error sources were all model errors. Therefore, the average model errors of the proposed method were considered to be 4.15% and 11.00% when calculating the boundary and depth of the goaf, respectively.

Study Area
The 15235 working face of the Jiulong mining area is located in Fengfeng coalfield, Hebei Province, China. It has a length of 935 m in the strike direction, a length of 142 m in the dip direction, an average coal mining thickness of 4.5 m, and an average mining depth of 740 m. The inclination of the working face is 13 • . When the inclination of the coal seam is less than 15 • , the surface deformation caused by mining can be regarded as being the same as that caused by horizontal coal seam mining [20]. As can be seen in Figure 9, there is a large fault and some old goafs located in the northwest and west, respectively, which may affect the location accuracy.

DInSAR Data Processing
Six images of RADARSAT-2, collected from 13th October 2015 to 5th March 2016, were used to obtain the subsidence basin using the multi-look fine mode. Five interferograms were formed with the parameters shown in Table 3. The coverage of the images is shown in Figure 9 with the red rectangle. SRTM data with a resolution of 3 seconds (90 m) were used as an external DEM to remove the terrain phase from the interferograms. The study area was small, so the atmospheric errors can be ignored [28]. The orbit error can also be ignored because of the highly precise measurement and control of the attitude of the RADARSAT-2 satellite [29]. Adaptive filtering and low-pass filtering were used to reduce noise effects. The interference results were used to calculate the surface LOS deformation. The time series surface subsidence can be obtained by projecting the LOS deformation in the vertical direction and stacking the images. If we assume that there is no horizontal movement, then the vertical deformation equals the LOS deformation divided by the cosine of the incident angle. The influence of this assumption is analyzed in Section 5.1.
As can be seen from Figure 10, the subsidence basin continued to expand, and the maximum subsidence value increased continuously and developed toward the southwest. Perhaps as a result of the influence of the surrounding geological conditions and the old goaf, the edge of the subsidence basin developed toward the northwest and southwest. This was inconsistent with the development direction of the subsidence basin center and will affect the final location results.

Goaf Location Process and Result
According to the development trend of the subsidence basin center in Figure 10, the mining strike azimuth was about 236 • . Two observing lines were made that were perpendicular to each other along the strike and dip directions and passed through the subsidence center. Sixty sampling points (black dots in Figure 10) were set at equal distances on each observing line. The subsidence value of the sampling point position was extracted from Figure 10 to obtain the subsidence curves in the strike and dip directions.
In Figure 11, it can be seen that the subsidence value at each sampling point increased over time, and the closer the point was to the subsidence basin center, the greater the subsidence value. The maximum subsidence value was 148 mm. Some points on the subsidence curves in the strike and dip directions showed positive values in the first three observation periods, indicating that there was uplift on the surface of the corresponding sampling points. The subsidence curves fluctuate in the strike and dip directions. The five subsidence curves in each graph fluctuate with the same trend, indicating that the cause of the fluctuation was the uneven surface subsidence, not the DInSAR observation error. In order to display the law of the subsidence curve more clearly and determine the inflection point more accurately, the subsidence curve was fitted. The end of the fitting curve in the strike direction did not tend to zero, which coincides with the subsidence graph in Figure 10. It further confirms that the edge of the subsidence basin was affected by geological structures or the goaf, which led to changes in its distribution. Therefore, when locating the boundary in the dip direction, the inflection point in the southeast was calculated first to determine the southeast boundary in the dip direction, and then the northwest boundary in the dip direction could be determined by the symmetry of the subsidence basin to avoid the error caused by the edge shifting of the subsidence basin. In this experiment, the proposed method was used to locate the goaf. First, the planimetric position of the working face 15,235 was determined according to the inflection points. With reference to the geological conditions of the Fengfeng mining area and other working faces, the deviation of the inflection point was set to 20 m, the boundary angle was set to 57 • , and the delay distance was set to 100 m. If the basin boundaries in the dip direction were used to calculate the mining depth, the result can be affected by the subcritical mining degree. If the southeast boundary was used to calculate the mining depth, the result can be affected by the subsidence delay. Therefore, the mining depth was calculated by the northeast boundary in the strike direction, and the final location result was obtained (Figure 12). The sampling points were arranged from northeast to southwest and from small to large in the strike direction, from southeast to northwest and from small to large in the dip direction.
It can be seen in Figure 12 that the calculated northeast mining boundary in the strike direction and the mining boundary in the dip direction were close to the real mining boundaries. The calculated southwest mining boundary was far from the real mining boundary, and this was because the geological structure and old goaf affect the subsidence basin. The calculated northeast mining boundary in the strike direction was located inside the real mining boundary, while the calculated mining boundary in the dip direction is located outside the real mining boundary; these findings were consistent with the simulation experiment, indicating that insufficient mining causes the calculated mining boundary to shift outward.

Precision Analysis
In order to evaluate the accuracy of this real data experiment, we calculated the absolute and relative error between the calculated mining boundary and the real mining boundary shown in Table 4 (the relative error is the ratio of the absolute error to the mining length in the strike or dip direction). The absolute errors of the southwest boundary and the boundaries in the dip direction were about 75 and 23 m, respectively, which were larger than the absolute errors of 50 m and 10 m in the simulation experiment. It was shown that the simulation experiment's neglected factors, such as DInSAR observation error and geological structure, had a great influence on the calculated results. Although the northwest boundary was obtained by symmetry, its error was different from that of the southeast boundary because the center of the subsidence basin was not located in the center of the working face. The error of the calculated mining depth was 9.07%, which was close to the error of 11.00% in the simulated experiment. Overall, the average relative error was 12.59%, and the average absolute error was about 38 m. The average error of the planar positioning was about 31 m. This accuracy was sufficient for engineering applications that require the positioning of the goaf (illegal mining location, goaf records, etc.). With an error of 67 m in depth determination and the geological data of the study area, we can determine the coal seam where the mining activity was located. Therefore, the accuracy of the method in this paper is acceptable.

Influence of DInSAR on the Location Result
Observed errors inevitably occur when using DInSAR to obtain surface deformation. In addition, the LOS deformation monitored by DInSAR is directly transformed into vertical deformation without considering the horizontal movement in the real data experiment. In order to verify the accuracy of the DInSAR monitoring results, the leveling data above the working face were used. The blue triangles represent the positions of the 51 leveling points arranged from northeast to southwest in Figure 10. The dates on which the first and last leveling data were measured were 10th October 2015 and 4th March 2016, respectively, which are similar to the acquisition times of the first and last SAR images. If the subsidence rate is regarded as a constant value, a time difference of two days will lead to an error of only 2.3%, so the influence of time inconsistency is very small. The comparison is shown in Figure 13. It can be seen from Figure 13 that the DInSAR monitoring result coincides with the leveling data. The maximum deviation is 31.5 mm, located at point 38; the minimum deviation is less than 0.1 mm, located at point 4; the average error is 8.5 mm; the root-mean-square error is 15.7 mm. With the horizontal movement ignored, the LOS deformation detected by DInSAR is considered to consist only of the component of the vertical surface deformation in the LOS direction. In fact, the components of the horizontal surface deformation in the LOS direction are also the components of the LOS deformation. Taking a point in the dip main section as an example (Figure 14), the LOS deformation obtained by DInSAR is composed of the component of vertical deformation in the LOS direction and the component of the horizontal deformation in the LOS direction (the blue arrow in Figure 14). In this example, the two components have the same direction. The LOS deformation is a fixed value that is obtained by DInSAR. When horizontal deformation is neglected, the vertical deformation will increase and become larger than the true value so that its component will be the same as the LOS deformation. The converse is also true: that is to say, if the component of the horizontal movement in the LOS direction is the same as the direction of the satellite incident angle, ignoring the horizontal movement will lead to the vertical deformation being larger than the true value; if the component of the horizontal movement in the LOS direction is opposite to the direction of the satellite incident angle, ignoring the horizontal movement will lead to the vertical deformation being smaller than the true value. The SAR data used in this experiment were orbit-lifting. The heading angle of the satellite is about −10 • , and the azimuth angle of the working face is about 236 • . The horizontal movement caused by mining the working face is toward the center of the subsidence basin in both the strike and dip directions. Therefore, the subsidence value is larger in the western part of the subsidence basin, while, in the east part of subsidence basin, the subsidence value is smaller. In Figure 13, points 1-40 are located in the eastern part of the subsidence basin, and most of the DInSAR data are smaller than the leveling data. Points 40-51 are located in the western part of the subsidence basin, and most of the DInSAR data are larger than the leveling data, which is consistent with the above conclusions. On the whole, the subsidence curve migrates to the southwest after neglecting the horizontal movement, resulting in the calculated mining boundary shifting toward the southwest. With respect to the real data experiment, the calculated northeast boundary is located inside the real boundary, which also confirms this conclusion. Therefore, with multi-platform InSAR data to calculate the three-dimension deformation, the precision of InSAR measurements will increase.

Influence of Geological Structure on the Location Result
It can be seen from Figures 10 and 12 that although the position of the maximum subsidence point of the surface subsidence basin is still inside the working face, the edge of the basin shifts toward the northwest and southwest. As can be seen in Figure 9, the position of the faults and old goafs are consistent with the basin's shifting direction, so the reason for the basin shifting may be that the mining activity of 15235 activates the overlying strata of faults and old goafs, which causes the subsidence beyond the surface. The shifting of the basin edge results in the migration of inflection points and boundary points. If the inflection point on the northwest side is used to locate the mining boundary, the absolute error is 123.95 m, and the relative error is 87.29%. This indicates that it is reasonable to directly use symmetry to locate the mining boundary.

Influence of Subcritical Extraction on the Location Result
According to experience from mining subsidence studies, the relationship between the mining degree and the ratio of mining size to depth is as follows [20]: when D 1 /H 0 , D 3 /H 0 < 1.2 ∼ 1.4, the mining degree is subcritical; when D 1 /H 0 , D 3 /H 0 = 1.2 ∼ 1.4, the mining degree is full; when D 1 /H 0 , D 3 /H 0 > 1.2 ∼ 1.4, the mining degree is super full, where D 1 and D 3 are the lengths in the dip and strike direction, respectively, and H 0 is the mining depth. The range 1.2 ∼ 1.4 is the experience value related to the properties of the overlying rock. The ratio of the length in the dip and strike directions to the mining depth of 15235 is 0.2 and 1.09, respectively, which means that the mining degree of 15235 is close to full in the strike direction and is subcritical in the dip direction. With the increase in mining degree, the inflection point moves from the outside to the inside of the mining boundary [20]. Therefore, in the simulation experiment and the real data experiment, the mining degrees in the dip direction are subcritical, so the calculated mining boundary is outside of the real mining boundary in the dip direction.

Conclusions
The proposed method uses the feature points of the basin monitored by DInSAR to locate the goaf. The method has the advantages of having a low economic cost and a high degree of automation compared with the traditional geophysical approaches. Although the PIM is taken as the theoretical basis, this method does not depend on a complex model algorithm or the model parameters of the probability integral method (such as subsidence coefficient, tangent of the main influence angle, etc.), and fewer parameters mean less uncertainty and better prospects for engineering applications. A simulation experiment and real data experiment were carried out to verify the reliability and accuracy of the proposed approach. The average relative errors of the simulation experiment and real data experiment are 6.43% and 12.59%. Because the basin is directly calculated in the simulated experiment, the errors come from the model. The relative error of the real data experiment is twice that of the simulated experiment; thus, the errors of the model and of the rest of the factors have the same order of magnitude. The experiments prove that the results of DInSAR monitoring can be used to locate the goaf because the DInSAR precision is good enough to accurately obtain the shape of the basin.
In this method, the accuracy of the shape and the feature points of the basin is more important than the deformation accuracy obtained by DInSAR monitoring. The mining degree affects the relative position of the inflection points and the goaf boundaries. The lower the mining degree, the farther the inflection points from the boundaries, so a subcritical extraction leads to the calculated goaf width being smaller than the actual value. Ignoring the horizontal movement causes the migration of the basin, and then the whole calculated goaf has an offset. Therefore, a high-precision observation method helps to improve the accuracy of this method. In our future work, we will focus on the influence of coal seam inclination, geological structure, and subcritical mining on the location results.