An Approach for Estimating Underground-Goaf Boundaries Based on Combining DInSAR with a Graphical Method

School of Resource & Environment and Safety Engineering, Hunan University of Science and Technology, Xiangtan 411201, China National-Local Joint Engineering Laboratory of Geo-Spatial Information Technology, Hunan University of Science and Technology, Xiangtan 411201, China Hunan Province Key Laboratory of Coal Resources Clean-utilization and Mine Environment Protection, Hunan University of Science and Technology, Xiangtan 411201, China


Introduction
e presence of underground goaf generally results in the collapse of the surrounding ground. is causes damage to houses and various types of ground-level transport routes, as well as polluting groundwater, thus seriously affecting the ecological environment of the mining area [1]. In particular, subsidence may even cause secondary geological disasters, such as mountain cracking, the collapse of constructions, landslides, mudslides, and earthquakes [2]. erefore, underground goaf is currently an important issue that restricts the development of mines and the urbanization of mining areas. Moreover, the distribution range and spatial morphological characteristics of underground goaf are key for the evaluation of the potential hazards and the establishment of countermeasures.
us, the determination of how to quantitatively evaluate the aforementioned range and characteristics of underground goaf is of crucial importance.
Existing goaf detection methods can be divided into geophysical and drilling technology. Geophysical techniques refer to the detection of geological conditions, such as lithology and geological structures, via the investigation of changes in various geophysical fields [3]. Typical geophysical techniques include gravimetric, electromagnetic, and geothermal methods [4][5][6][7]. Drilling technology refers to the use of drilling equipment to drill through rock formations at a predetermined location in order to extract physical samples and other relevant data for experiments [8]. Central drilling technology methods include the flushing fluid, sonic velocity, and ultrasonic imaging methods [9][10][11]. Generally speaking, although these technologies are able to obtain the rough location of the underground goaf along with its spatial distribution, they have several limitations. For example, such technology is only suitable for small-scale detection, which can be costly, time-consuming, and inefficient. At the same time, the approximate geographical location of the goaf is required as prior information before detection. Moreover, over the past few decades, with the booming mining industry (including legal and illegal mining), numerous unknown underground goafs have emerged globally. Traditional geophysical and drilling techniques are unable to meet the current wide-ranging and time-sensitive goaf detection needs.
e interferometric synthetic aperture radar (InSAR) is a popular remote sensing technology that can survey large-scale surface deformation using two-or multi-view synthetic aperture radar (SAR) images. Massonnet et al. combined two-scene ERS-1 images before and after an earthquake over the earthquake-prone area of Landers, California, with the DEM (Digital Elevation Model, DEM) of the region to successfully observe the surface deformation caused by earthquakes. Since then, the differential interference synthetic aperture radar technique (DInSAR) has been widely used for the monitoring of surface deformation caused by geological tectonic movements, such as seismic monitoring [12,13], glacial drift monitoring [14,15], volcanic eruption monitoring [16,17], landslide monitoring [18,19], frozen soil degradation, and expansion monitoring [20][21][22]. e use of DInSAR in mining subsidence is also a hot research topic and includes applications in the detection of subsidence caused by mining [23,24], the prediction of mining-induced geological disasters [25,26], and the assessment of the impact of mining on surface buildings [27,28].
Following the stabilization of the ground subsidence basin caused by underground mining, a close relationship between the size and spatial distribution of the basin and geological mining conditions can be observed.
is is denoted as the subsidence rule. More specifically, the geological mining conditions and the size of the goaf determine the spatial distribution of the ground subsidence basin.
at is, under specific geological mining conditions, the spatial distribution characteristics of the ground subsidence basin can also reflect the size of the goaf. DInSAR can not only survey the movement and deformation of the subsidence area but also determine the distribution state of the whole ground subsidence basin.
us, it can potentially be applied for the detection of underground goaf distributions.
At present, research on the detection of the geometrical distribution of underground goaf using DInSAR is rare. In 2013, Zhe et al. [29] proposed a method known as "DInSAR-Based Illegal Gob Detection System (DIMDS)," which utilizes the DInSAR to obtain the surface deformation value along the radar line-of-sight (LOS) direction. e method then estimates the geodetic coordinates of the center of the ground subsidence basin based on the mining subsidence theory and subsequently determines the extent of underground goaf based on the ground subsidence basin boundary. However, the DIMDS method has a key limitation. It assumes that the center of the underground goaf is on the same vertical line as the center of the ground subsidence basin. Yet for most cases this is not true, as most coal seams are not completely horizontal; as long as the coal seam is tilted, the movement deformation caused by mining will not propagate to the surface in the vertical direction. is discrepancy seriously affects the accuracy of the boundary of the underground goaf estimated using the DIMDS method.
In response to this problem, Yang et al. [30] defined the geometrical distribution characteristics of the underground goaf using eight geometric parameters (length, width, height, dip angle, azimuth angle, mining depth, and two central geodetic coordinates). Based on the probability integral method, the relationship between these geometric parameters and the LOS deformation value of the ground subsidence basin is then established. Following this, a mathematical model based on the simulated annealing algorithm (SAGA) [31] is used to invert the eight geometric parameters from a large number of LOSdeformed observations obtained by DInSAR. is method considers the general law of mining subsidence and makes full use of the surface information of the surface deformation obtained by DInSAR. us, prediction accuracy is greatly improved compared to the DIMDS method. However, due to the large number of inversion parameters, the stability of the mathematical model is low. In particular, the selection of the initial value of the model will have a great impact on the inversion results. Yet in real-world applications, the initial value of a goaf is difficult to determine. Second, due to the coherence losses of the echo signals from a SAR image pair, a large number of null values are observed in the whole basin acquired by DInSAR. is leads to errors in the inversion model solution process and subsequent incorrect results.
At present, most coal mines use long-arm coal mining and full slumping to manage the roof. Under the condition of uniform mining thickness, the main factors affecting the subsidence rule are coal seam dip angle, mining area size (length and width), and mining depth [32]. ese parameters can reflect the boundary range and basic shape of the goaf. How to acquire these unknown parameters quickly and effectively is of great significance, to determine the boundaries of underground goaf. erefore, this paper proposes an underground goaf boundary detection method by combining DInSAR with a graphical method. Our proposed method first uses DInSAR to obtain the subsidence information of the whole basin in the subsidence area and based on this, the subsidence contour map is generated. e two main sections (the strike section and incline section) of the underground goaf are determined according to the contour map, and the subsidence curve of these two sections is drawn. Finally, based on the subsidence rule of coal mining, the length, coal seam dip angle, mining area size, and mining depth of the underground goaf are derived using a graphical method. Moreover, after geocoding the size and location of the underground goaf, the boundary is determined. e paper is organized as follows. In Section 2, we introduce the method used to estimate the boundary range of the underground goaf based on DInSAR. We then verify the proposed method using simulation and measured data in Section 3. Finally, Section 4 concludes the study.

Methods
After the underground ore body is mined out, a cavity (i.e., a goaf ) will be left inside the rock strata, and the original stress equilibrium will be destroyed, which will result in the stress redistribution in country-rock taking place to reach a new equilibrium. is is a very complex variation of physical and mechanical process, causing movement and damage to the overlying strata. With the continuous progress of mining activities, the goaf has expanded to a certain extent, and then the movement of rock strata will develop to the surface. For most coal mines, the ratio of mining depth and thickness is usually large. At this time, the surface deformation is continuous in space and gradual in time, and it has obvious regularity. Because the geometry of the subsidence basin is closely related to the spatial distribution of the goaf, it can be obtained by measurement to infer the spatial distribution characteristics of the goaf. e framework of the proposed method for the detection of an underground goaf boundary based on DInSAR is presented in Figure 1. First, DInSAR is used to obtain the subsidence information of the whole basin in the subsidence area. Based on this, the subsidence contour map is created. Second, the strike and incline sections of the underground goaf are determined from the contour map, and the subsidence curves of these two sections are drawn. Finally, based on the coal mining subsidence rule, the length, coal seam dip angle, mining area size, and mining depth of the underground goaf are calculated using a graphical method.

Extracting Subsidence Data Using
DInSAR. Synthetic aperture radar differential interference (DInSAR) acquires surface deformation information via the differential interference processing of two SAR images of the same area in different phases. If the spatial baseline between the SAR image pair is small enough, the repeated deformation observation can be used to survey surface deformation, as shown in Figure 2. More specifically, the surveyed deformation is the direction of the line-of-sight of the sensor. In theory, it is possible to survey changes within the millimeter range using DInSAR [33][34][35][36].
During the imaging of the two scenes, the surface is deformed. According to the vector relationship, these vector sums are equal to zero. erefore the following formula can be obtained: where R i → is the line-of-sight vector, i.e., the vector between the ith subantenna and the ground point, and ΔR mov ������→ is the deformation vector of a ground point during the period from t 1 to t 2 , and B → is the distance vector between the antennas of the two images. e interference phase ϕ of the two images can be expressed as follows: where λ is the radar wavelength, D → is the displacement vector of the ground point during imaging, ρ i � |R i → | is the distance between the antenna and the ground point, and < > indicates point multiplication.
For a space-based SAR system, the deformation of the ground point and the spatial baseline of the two scene images is much smaller than the distance between the ground point and the satellite. us, equation (3) can be expressed as where ϕ def is the surface deformation phase and ϕ topo is the terrain phase. e phase derived from the differential interference consists of two parts: the terrain phase ϕ topo and the deformation phase ϕ def . e terrain phase determined via DEM simulation can be applied to eliminate the influence of terrain factors (known as the two-track method), and thus the ground point deformation of the radar line of sight can be obtained as follows: where ϕ sim represents the terrain phase simulated by external DEM.

Estimating Underground-Goaf Boundaries Based on the Graphical Approach.
A total of four parameters, namely, length (L 1 ), width (L 2 ), depth (H 0 ), and dip angle (α), are needed to characterize the goaf boundary, where length and width are used to determine the size of the goaf, and depth and dip angle are used to determine the position of the goaf. is paper aims to estimate these parameters based on the graphical approach. e subsidence data of the mining area used in this study was obtained using DInSAR. Preprocessing of the data includes the geocoding and conversion to ArcGIS (10.0, ERSI) or CASS (9.0, SOUTHIS) readable formats. e subsidence contour map of the mining area is then obtained by ArcGIS or CASS. According to the distribution law of full-area mining subsidence [32], once the subsidence basin is stable, the subsidence contour is approximately elliptical. e subsidence value peaks at the center of the ellipse, while the farther away from the center, the smaller the subsidence value. When the working surface is approximately rectangular, the long semiaxis direction of the ellipse points to the strike section of the working surface, and the direction perpendicular to the strike section is taken as the incline section of this working surface. So we can calibrate the position of the maximum subsidence point on the subsidence contour map, and then make two straight lines along Advances in Civil Engineering 3 the major and minor semiaxes of the ellipse, respectively. ese two lines are the main sections of the subsidence basin. e subsidence values at each point on the two main sections are extracted and fitted based on the basic principle of the probability integration method. Following this, the subsidence curves of the main section are drawn, as shown in Figure 3.

Strike Section Boundary
Estimate. Two points of subsidence values of 0.84Wm and 0.16Wm (where Wm is the maximum subsidence value) are extracted from the strike section, and the horizontal distance between these two points is equal to 0.8r (r is the main influence radius, and r 0 is the main influence radius of the strike section). If the r values obtained on both sides of the maximum subside point are different, the average of the two values is taken. And then we can calculate the depth (H) according to the definition of the parameter tan β (tangent of main effect angle) from the probability integration method, as shown in the following equation: where H is mining depth. us, the mining depth of the strike section (i.e., average depth) H 0 can be calculated and used to create a horizontal plane below the strike section with depth H 0 . Following this, two inflection points at both ends of the strike section are found on the subsidence curve. ese two inflection points are the computational boundary of the working surface. e deviation of the inflection point (S) can be calculated using the following equation: where k L is the coefficient of the inflection point offset. e mining boundary of the strike section, presented in Figure 4, is obtained by shifting the computational boundary outward by S 0 , and the distance between these two boundary points is the strike section length L 1 .

Incline Section Boundary
Estimate. For this step, four points of the subsidence values of 0.84Wm and 0.16Wm are extracted in the rise and dip direction of the incline section, respectively. e horizontal distance between these two points on the same side is equal to 0.8r (r 1 and r 2 are the main influence radii of the rise and dip directions, respectively). e mining depths H 1 and H 2 of the rise and dip directions, respectively, are then calculated according to r 1 and r 2 by equation (5). Following this, two horizontal planes below the incline section with depths of H 1 and H 2 , respectively, are determined.
Next, two inflection points at both sides of the incline section on the subsidence curve are found. Two vertical lines are drawn through these two points. e vertical line of the rise direction intersects the rise direction depth at point N, and the vertical line of the dip direction intersects the dip direction depth at point M. ese two intersection points are then connected, whereby the angle between the intersection line and the horizontal direction is the dip angle of the coal seam (α) (see Figure 5).
According to the definition of the probability integration method, the propagation angle is given as where K is the propagation coefficient of extraction. An auxiliary line is then drawn in the depth direction through the inflection point, such that the angle between the auxiliary line and the horizontal direction is θ. Moreover, the intersection of the auxiliary and horizontal lines of the depth is the calculated boundary of the goaf. Following this, the offset distances, S 1 and S 2 , of the dip and rise direction, are then calculated using the depths H 1 and H 2 , respectively, based on the geological conditions of the mining area. e mining boundary is obtained by shifting the calculated boundary outwards by S 1 and S 2 , and the distance between these two boundary points is the incline section length L 2 . e results are presented in Figure 6. Based on these results, the boundary of the entire goaf can now be determined. e determined geometric parameters of the goaf are the average depth H 0 , dip direction depth H 1 , rise direction depth H 2 , strike distance L 1 , slope length L 2 , and dip angle α. In order to estimate the boundary range of the goaf, the direction of the two main sections is determined according to the subsidence contour. e subsidence value of each point on the main section is then extracted. Following this, the subsidence value is fitted using the probability integral method, and the subsidence and slope curves of the strike and inclined sections are created.

Strike Section Boundary Estimate.
Two inflection points at both ends of the incline section are determined on the subsidence curve, i.e., the calculation boundary of the working face. e subsidence values, 0.84Wm and 0.16Wm, are extracted from the subsidence curve. e distance between these two points is 0.8r 0 , which is subsequently used to calculate r 0 . According to the probability integral method parameter tan β, the strike depth is determined as H 0 � r 0 ·tan β � 197 m. e inflection offset point of the simulated working face is determined as S 0 � 0.15H 0 � 29.55 m.
us, the actual boundary of the simulated working face can be  Advances in Civil Engineering determined following an outward shift of the calculated boundary by S 0 .

Incline Section Boundary Estimate.
e subsidence values of 0.84Wm and 0.16Wm are extracted from the subsidence curve. e distance between these two points is 0.8r, which is then used to calculate the dip and rise directions of the main influence radii of r 1 and r 2 , respectively. Following this, tan β is used to calculate the dip direction depth of H 1 � 219.6 m and the rise direction depth of H 2 � 176.2 m. Next, two inflection points at both ends of the incline section on the subsidence curve are determined. Two vertical lines are drawn to the horizontal line through these two points, respectively. e vertical line of the rise direction intersects the rise direction depth, while the other line intersects with the dip direction depth. ese two intersection points are connected to obtain the dip angle α, i.e., the angle between the intersection line and the horizontal direction, at α � 18.3°. Based on the dip angle and mining reflection coefficients k L , the propagation angle θ is determined, and the incline section boundary of the working face is thus identified.
Calculating the inflection offset points of the dip direction S 1 and the rise direction S 2 allows for the derivation of the actual boundary of the incline section working face after the calculated boundary is shifted outward by S 1 and S 2 , respectively.

Experimental Results.
Using the estimated geometric parameters of the goaf, the mining boundary can be determined. e final estimate results are shown in Figure 7.

Accuracy Assessment.
Next, we can compare the predicted boundary to the simulated boundary. Figure 8 shows the difference between the simulated and predicted boundary. e predicted boundary (red rectangle) and the simulated boundary (cyan rectangle) almost overlap, and the observed degree of coincidence is high. is indicates that the boundary predicted by the proposed method is reliable.
In order to quantitatively evaluate the accuracy of the predicted boundary, we calculate the deviation between the predicted and simulated goaf geometric parameters (as reported in Table 1), as well as the relative error between the predicted and simulated parameters. e relative error K is calculated as follows: where f is the deviation between the predicted and simulated parameters, and G sim is the simulated geometrical parameter of the goaf. e results of the comparison are reported in Table 1. e predicted geometric parameters are generally consistent with the simulated geometric parameters, with an average relative error of 2.2%. is again demonstrates the ability of the proposed method to accurately predict the boundary range of the underground goaf.

Study Area.
e Pangzhuang Coal mine is located in the Jiuli District of Xuzhou City, 13 km away from the center of Xuzhou City. e surface is mainly covered by vegetation and buildings (Figure 9). Over the years, the coal mine has produced a large number of underground goafs, causing serious problems in the area, including damage to buildings and roads, and water inrush from mines. In order to minimize the adverse effects of the underground goaf, its boundary must be identified. We selected No. 7503 working face (represented in Figure 9 by a blue rectangular) of the Pangzhuang 7 coal seam to test the proposed method. e working face adopts the fully mechanized long arm mining process and completely manages the roof. e mining area began to operate in the 1980s and has a sufficient amount of geological and measured data of the underground goaf (Table 2). us, this working face possesses the necessary conditions required for the accuracy verification of the proposed method.

Data Processing.
A total of 13 C-band ENVISAT ASAR images covering the Pangzhuang coal mining area were selected, with acquisition times from 2009-01-20 to 2010-10-12.
e wavelength of the C-band is 5.6 cm, the incident angle is 22.78°, and the ground resolution is 20.12 m. e image coverage area is shown in Figure 10 Table 3. ENVI SARscape (2014, sarmap) is used to implement the DInSAR processing. First, the SAR images are cut under the boundary of the study area, and the size of the clipping data is 884 × 4717 pixels. After that, based on the characteristics of the ASAR imagery, the looks of azimuth and distance are set to 6 and 1, respectively, in multi-looking processing, to get a cartographic resolution of 25 m. SARscape produces the interference image pair and uses the 90 m resolution SRTM-3 digital elevation data (V.4) to remove the horizon effect of the interferogram. Following this, the Goldstein method is used to filter the deleveling interferogram to reduce coherent noise. e phase unwrapping process is then applied using the minimum cost flow method, these pixels with low coherence will be masked, and the coherence threshold is set to 0.2. e GCPs (ground control points) are selected based on their suitability for most unwrapped phases. ese GCP points are used to redefine the baseline parameters for orbital refinement. A polynomial model is then used to calculate the phase offset, completing the releveling process. Finally, the residual phase is transformed into a shape variable and geocoded to obtain the deformation result of the LOS direction. e subsidence value can be obtained by projecting the LOS deformation in the vertical direction, which is equal to the LOS deformation divided by the cosine of the incident angle, as shown in Figure 11.
From Figure 11, it can be seen that several deformation abnormalities can be observed. But only the area of 7503 Advances in Civil Engineering 7 working face (a black rectangular in Figure 11) presents a complete subsidence process during the imaging period (i.e., subsidence begins, subsidence development, and subsidence stability). It indicates that underground mining operations are present in the area during this period and the subsidence has been stable. erefore, we can attempt to estimate the boundary extent of an underground goaf in the area using the proposed method.

Experimental
Results. e method proposed in this paper aims to estimate the boundary of the underground goaf in the study area. e time series surface subsidence can be obtained by stacking the images (as in Figure 12). ArcGIS is used to process the subsidence results obtained by InSAR. At first, load the subsidence data in ArcGIS. After that, select the "Raster Surface-Contour" tool from the "3D Analyst" menu in the toolbox. en allow for the generation of the subsidence contour map after the tool finishing running. e maximum subsidence point in the subsidence basin is then identified. Second, according to the geometry of the contour, the direction of the main section of the underground goaf is determined, and the two main sections are drawn through the maximum subsidence point. Finally, the subsidence values of various ground points are extracted from the main section, and the subsidence curves of the two main sections are drawn (as in Figures 13).   e subsidence values of 0.84Wm and 0.16Wm are extracted from the strike and inline sections, respectively. Next, the distance between the two points is measured to calculate the average depth of H 0 � 305.6 m, the dip direction depth of H 1 � 318.9 m, and the rise direction depth of H 2 � 300.2 m. en determine the position of an inflection point on the subsidence curve. According to the rock motion report (unpublished) of the mining area, the coefficient of the inflection point offset is 0.13 (see from Table 2); thus we can calculate the inflection point offset by formula (6) (i.e., S 0 � 0.13 × H 0 � 39.7 m, S 1 � 0.13 × H 1 � 41.5 m, S 2 � 0.13 × H 2 � 39.0 m). e inflection point of the strike section is used to determine the calculated boundary of the strike section, and S 0 is used to determine the mining boundary of the strike section. e distance between these two boundary points is the strike section length L 1 � 703.4. Following this, the inflection point of the incline section and the depth of the dip and rise directions are used to determine the coal seam inclination as α � 4.3°. en the propagation angle was calculated by equation (7). And the calculated boundary of the incline section could be determined by the propagation angel and the inflection point. e mining boundary of the incline section is obtained by shifting the calculated boundary outwards by S 1 and S 2 . e distance between these two boundary points is the strike section length L 2 � 176.9. Finally, the estimated results of the boundaries of the underground goaf can be obtained. e boundaries of the probe are marked with a red rectangle in Figure 14.

Precision Evaluation.
In order to evaluate the accuracy of the estimate results, we compared the boundary obtained by standard geophysical techniques (blue) with that obtained by the proposed method (red) as shown in Figure 15. e boundary of the goaf predicted by DInSAR is consistent with that predicted by the geophysical techniques, indicating that the method proposed in this paper is reliable. To quantitatively evaluate the reliability of the method, we calculate the difference between the geometric parameters of the underground goaf estimated by our method and the actual geometric parameters of the goaf (see Table 2    Advances in Civil Engineering relative error between the predicted and actual parameters (see Table 4). As reported in Table 4, the relative errors of the four predicted geometric parameters are below 10%, with an average value of 4.55%. e average relative error is approximately 3.7%, respectively. Compared with the method proposed by Yang et al. [30], the relative errors of the four geometric parameters L 1 , L 2 , H 0 , and α are reduced by 7.375%, suggesting that the method proposed in this paper is highly accurate for estimating two-dimensional boundaries.

Conclusion
is paper reports an approach for estimating the boundary range of underground goaf using the DInSAR by combining a graphical method. e proposed method fully considers the spatial distribution relationship between surface deformation obtained by DInSAR and the geometry of    underground goaf and combines the relevant principles of the probability integral method model. Our results show that the predicted boundary of the underground goaf strongly agrees with the measured boundary. e average relative errors of the estimated and measured data are approximately 2.2% and 3.7%, respectively.
However, it is worth noting that DInSAR acquires the surface deformation of the LOS direction; thus, the geometry of the generated contour will be deflected toward the radar line of sight, leading to a deviation in the direction of the main section of the predicted goaf. A higher deviation between the predicted and measured data is observed compared to the simulated data. In addition, the strike section length L 1 , the incline section length L 2 , and the depths H 0 , H 1 , and H 2 (predicted by the measured data) are all underestimated. is is because the deformation caused by underground mining has not fully propagated to the surface when SAR images are acquired, resulting in the underestimation of the boundary size. Future research should focus on how to solve these problems in order to further improve the prediction accuracy. Considering the low cost and large spatial coverage of DInSAR, it can be an effective tool for the estimate of the boundary range of underground goaf.

Data Availability
e data used to support the findings of this study are included in the article, which is based on the SAR images and the external DEM of the study area.

Conflicts of Interest
e authors declare no conflicts of interest regarding the publication of this paper.