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

Authors: Pu Bu; Chaokui Li; Mengguang Liao; Wentao Yang; Chuanguang Zhu; Jun Fang
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

The goaf left behind after mining has the potential to induce serious geological disasters due to the damaged internal structure of the rock. Estimating the boundary of the underground goaf can effectively control the occurrence of such disasters. However, traditional geophysical methods are inefficient and expensive and are particularly difficult to apply for a wide detection range. This paper proposes a new method for estimating the boundary of underground goaf using the differential interference synthetic aperture radar technique (DInSAR). More specifically, DInSAR is used to obtain the isoline of the subsidence basin above the goaf, and the direction of the two main sections of the goaf is then determined according to the basic law of mining subsidence. Following this, the basic principles of the probability integral and the graphical methods are combined to determine the mining boundary of the strike section and the incline section of the goaf. Finally, six geometric parameters reflecting the boundary of the goaf are obtained. Experiments on simulated and measured data indicate that the proposed method is feasible, with the average relative errors of the simulated and measured data reaching and maintained at 2.2% and 3.7%, respectively.


---

## Body

## 1. Introduction

The presence of underground goaf generally results in the collapse of the surrounding ground. This 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]. Therefore, 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. Thus, 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–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–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.The 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–22]. The 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. This 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. That 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. Thus, 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. The 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. This 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 LOS-deformed observations obtained by DInSAR. This method considers the general law of mining subsidence and makes full use of the surface information of the surface deformation obtained by DInSAR. Thus, 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. This 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]. These 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. Therefore, 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. The 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.The paper is organized as follows. In Section2, 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.


## 2. 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. This 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.The framework of the proposed method for the detection of an underground goaf boundary based on DInSAR is presented in Figure1. 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.Figure 1
Flow chart of the proposed method.


### 2.1. 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 Figure2. 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–36].Figure 2
DInSAR schematic.During the imaging of the two scenes, the surface is deformed. According to the vector relationship, these vector sums are equal to zero. Therefore the following formula can be obtained:(1)R2⟶=R1⟶+ΔRmov⟶−B⟶,where Ri⟶ is the line-of-sight vector, i.e., the vector between the ith subantenna and the ground point, and ΔRmov⟶ is the deformation vector of a ground point during the period from t1 to t2, and B⟶ is the distance vector between the antennas of the two images.The interference phaseϕ of the two images can be expressed as follows:(2)ϕ=4πλρ2−ρ1=4πλR1⟶+D⟶−B⟶,R1⟶+D⟶−B⟶1/2−ρ1,where λ is the radar wavelength, D⟶ is the displacement vector of the ground point during imaging, ρi=Ri⟶ 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. Thus, equation (3) can be expressed as(3)ϕ=4πλR1⟶,D⟶−R1⟶,B⟶=ϕdef+ϕtopo,where ϕdef is the surface deformation phase and ϕtopo is the terrain phase.The phase derived from the differential interference consists of two parts: the terrain phaseϕtopo and the deformation phase ϕdef. The 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:(4)ΔRmov=λ4πϕ−ϕsim,where ϕsim represents the terrain phase simulated by external DEM.


### 2.2. Estimating Underground-Goaf Boundaries Based on the Graphical Approach

A total of four parameters, namely, length (L1), width (L2), depth (H0), 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. This paper aims to estimate these parameters based on the graphical approach.The 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. The 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. The 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 the major and minor semiaxes of the ellipse, respectively. These two lines are the main sections of the subsidence basin. The 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.Figure 3
Subsidence curve of the main section. I-I′ is the strike section, II-II′ is the incline section, the red line denotes the subsidence curve of the strike section, and the green line is the subsidence curve of the incline section.


#### 2.2.1. 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 r0 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:(5)tanβ=Hr,where H is mining depth. Thus, the mining depth of the strike section (i.e., average depth) H0 can be calculated and used to create a horizontal plane below the strike section with depth H0.Following this, two inflection points at both ends of the strike section are found on the subsidence curve. These two inflection points are the computational boundary of the working surface. The deviation of the inflection point (S) can be calculated using the following equation:(6)S=kLH,where kL is the coefficient of the inflection point offset. The mining boundary of the strike section, presented in Figure 4, is obtained by shifting the computational boundary outward by S0, and the distance between these two boundary points is the strike section length L1.Figure 4
Estimate of the strike section boundary. A and D are the actual boundaries of the goaf, B and C are the calculated boundaries of the goaf, and the distance between the two points A and D is the strike section lengthL1.


#### 2.2.2. 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. The horizontal distance between these two points on the same side is equal to 0.8r (r1 and r2 are the main influence radii of the rise and dip directions, respectively). The mining depths H1 and H2 of the rise and dip directions, respectively, are then calculated according to r1 and r2 by equation (5). Following this, two horizontal planes below the incline section with depths of H1 and H2, 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. The 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. These 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).Figure 5
Derivation of the dip angle of the goaf. N is the intersection of the vertical line of the rise direction depth, M is the intersection of vertical line of the dip direction depth, andα is the dip angle of the coal seam.According to the definition of the probability integration method, the propagation angle is given as(7)θ=90∘−Kα,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, S1 and S2, of the dip and rise direction, are then calculated using the depths H1 and H2, respectively, based on the geological conditions of the mining area. The mining boundary is obtained by shifting the calculated boundary outwards by S1 and S2, and the distance between these two boundary points is the incline section length L2. The results are presented in Figure 6.Figure 6
Estimate of the incline section boundary. A′ and D′ denote the actual boundaries of the goaf in the dip and rise directions, respectively, while C′ and D′ are the calculated boundaries of the goaf in the dip and rise directions, respectively. The distance between A′ and D′ is the length of the incline section of the faceL2.Based on these results, the boundary of the entire goaf can now be determined. The determined geometric parameters of the goaf are the average depthH0, dip direction depth H1, rise direction depth H2, strike distance L1, slope length L2, and dip angle α.


## 2.1. 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 Figure2. 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–36].Figure 2
DInSAR schematic.During the imaging of the two scenes, the surface is deformed. According to the vector relationship, these vector sums are equal to zero. Therefore the following formula can be obtained:(1)R2⟶=R1⟶+ΔRmov⟶−B⟶,where Ri⟶ is the line-of-sight vector, i.e., the vector between the ith subantenna and the ground point, and ΔRmov⟶ is the deformation vector of a ground point during the period from t1 to t2, and B⟶ is the distance vector between the antennas of the two images.The interference phaseϕ of the two images can be expressed as follows:(2)ϕ=4πλρ2−ρ1=4πλR1⟶+D⟶−B⟶,R1⟶+D⟶−B⟶1/2−ρ1,where λ is the radar wavelength, D⟶ is the displacement vector of the ground point during imaging, ρi=Ri⟶ 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. Thus, equation (3) can be expressed as(3)ϕ=4πλR1⟶,D⟶−R1⟶,B⟶=ϕdef+ϕtopo,where ϕdef is the surface deformation phase and ϕtopo is the terrain phase.The phase derived from the differential interference consists of two parts: the terrain phaseϕtopo and the deformation phase ϕdef. The 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:(4)ΔRmov=λ4πϕ−ϕsim,where ϕsim represents the terrain phase simulated by external DEM.


## 2.2. Estimating Underground-Goaf Boundaries Based on the Graphical Approach

A total of four parameters, namely, length (L1), width (L2), depth (H0), 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. This paper aims to estimate these parameters based on the graphical approach.The 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. The 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. The 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 the major and minor semiaxes of the ellipse, respectively. These two lines are the main sections of the subsidence basin. The 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.Figure 3
Subsidence curve of the main section. I-I′ is the strike section, II-II′ is the incline section, the red line denotes the subsidence curve of the strike section, and the green line is the subsidence curve of the incline section.


### 2.2.1. 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 r0 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:(5)tanβ=Hr,where H is mining depth. Thus, the mining depth of the strike section (i.e., average depth) H0 can be calculated and used to create a horizontal plane below the strike section with depth H0.Following this, two inflection points at both ends of the strike section are found on the subsidence curve. These two inflection points are the computational boundary of the working surface. The deviation of the inflection point (S) can be calculated using the following equation:(6)S=kLH,where kL is the coefficient of the inflection point offset. The mining boundary of the strike section, presented in Figure 4, is obtained by shifting the computational boundary outward by S0, and the distance between these two boundary points is the strike section length L1.Figure 4
Estimate of the strike section boundary. A and D are the actual boundaries of the goaf, B and C are the calculated boundaries of the goaf, and the distance between the two points A and D is the strike section lengthL1.


### 2.2.2. 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. The horizontal distance between these two points on the same side is equal to 0.8r (r1 and r2 are the main influence radii of the rise and dip directions, respectively). The mining depths H1 and H2 of the rise and dip directions, respectively, are then calculated according to r1 and r2 by equation (5). Following this, two horizontal planes below the incline section with depths of H1 and H2, 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. The 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. These 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).Figure 5
Derivation of the dip angle of the goaf. N is the intersection of the vertical line of the rise direction depth, M is the intersection of vertical line of the dip direction depth, andα is the dip angle of the coal seam.According to the definition of the probability integration method, the propagation angle is given as(7)θ=90∘−Kα,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, S1 and S2, of the dip and rise direction, are then calculated using the depths H1 and H2, respectively, based on the geological conditions of the mining area. The mining boundary is obtained by shifting the calculated boundary outwards by S1 and S2, and the distance between these two boundary points is the incline section length L2. The results are presented in Figure 6.Figure 6
Estimate of the incline section boundary. A′ and D′ denote the actual boundaries of the goaf in the dip and rise directions, respectively, while C′ and D′ are the calculated boundaries of the goaf in the dip and rise directions, respectively. The distance between A′ and D′ is the length of the incline section of the faceL2.Based on these results, the boundary of the entire goaf can now be determined. The determined geometric parameters of the goaf are the average depthH0, dip direction depth H1, rise direction depth H2, strike distance L1, slope length L2, and dip angle α.


## 2.2.1. 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 r0 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:(5)tanβ=Hr,where H is mining depth. Thus, the mining depth of the strike section (i.e., average depth) H0 can be calculated and used to create a horizontal plane below the strike section with depth H0.Following this, two inflection points at both ends of the strike section are found on the subsidence curve. These two inflection points are the computational boundary of the working surface. The deviation of the inflection point (S) can be calculated using the following equation:(6)S=kLH,where kL is the coefficient of the inflection point offset. The mining boundary of the strike section, presented in Figure 4, is obtained by shifting the computational boundary outward by S0, and the distance between these two boundary points is the strike section length L1.Figure 4
Estimate of the strike section boundary. A and D are the actual boundaries of the goaf, B and C are the calculated boundaries of the goaf, and the distance between the two points A and D is the strike section lengthL1.


## 2.2.2. 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. The horizontal distance between these two points on the same side is equal to 0.8r (r1 and r2 are the main influence radii of the rise and dip directions, respectively). The mining depths H1 and H2 of the rise and dip directions, respectively, are then calculated according to r1 and r2 by equation (5). Following this, two horizontal planes below the incline section with depths of H1 and H2, 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. The 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. These 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).Figure 5
Derivation of the dip angle of the goaf. N is the intersection of the vertical line of the rise direction depth, M is the intersection of vertical line of the dip direction depth, andα is the dip angle of the coal seam.According to the definition of the probability integration method, the propagation angle is given as(7)θ=90∘−Kα,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, S1 and S2, of the dip and rise direction, are then calculated using the depths H1 and H2, respectively, based on the geological conditions of the mining area. The mining boundary is obtained by shifting the calculated boundary outwards by S1 and S2, and the distance between these two boundary points is the incline section length L2. The results are presented in Figure 6.Figure 6
Estimate of the incline section boundary. A′ and D′ denote the actual boundaries of the goaf in the dip and rise directions, respectively, while C′ and D′ are the calculated boundaries of the goaf in the dip and rise directions, respectively. The distance between A′ and D′ is the length of the incline section of the faceL2.Based on these results, the boundary of the entire goaf can now be determined. The determined geometric parameters of the goaf are the average depthH0, dip direction depth H1, rise direction depth H2, strike distance L1, slope length L2, and dip angle α.


## 3. Experimental Analysis

### 3.1. Simulation Experiment

Assume that a coal face has the following geometric parameters: coal thicknessm = 2.2 m, dip angle α = 18°, strike length L1 = 240 m, inclined length L2 = 138.5 m, dip section boundary depth H1 = 220 m, and rise direction boundary depth H2 = 187 m. Moreover, the parameters of the probability integral method are as follows: subsidence coefficient q = 0.65, tan β = 2.0, inflection point offset S = 0.15 H, and mining reflection coefficient kL = 0.8. Matlab (2015b, MathWorks) is used to apply the probability integral method to simulate the vertical displacement of the surface as a result of mining. After adding a random error to the simulated value, the subsidence contour of the area is obtained.In order to estimate the boundary range of the goaf, the direction of the two main sections is determined according to the subsidence contour. The 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.


#### 3.1.1. 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. The subsidence values, 0.84Wm and 0.16Wm, are extracted from the subsidence curve. The distance between these two points is 0.8r0, which is subsequently used to calculate r0. According to the probability integral method parameter tan β, the strike depth is determined as H0 = r0·tan β = 197 m.The inflection offset point of the simulated working face is determined asS0 = 0.15H0 = 29.55 m. Thus, the actual boundary of the simulated working face can be determined following an outward shift of the calculated boundary by S0.


#### 3.1.2. Incline Section Boundary Estimate

The subsidence values of 0.84Wm and 0.16Wm are extracted from the subsidence curve. The 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 r1 and r2, respectively. Following this, tan β is used to calculate the dip direction depth of H1 = 219.6 m and the rise direction depth of H2 = 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. The vertical line of the rise direction intersects the rise direction depth, while the other line intersects with the dip direction depth. These 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 kL, 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 directionS1 and the rise direction S2 allows for the derivation of the actual boundary of the incline section working face after the calculated boundary is shifted outward by S1 and S2, respectively.


#### 3.1.3. Experimental Results

Using the estimated geometric parameters of the goaf, the mining boundary can be determined. The final estimate results are shown in Figure7.Figure 7
Estimate results. A and D are the actual boundaries of the strike section, C and B are the calculated boundaries of the strike section, and the distance between A and D is the length of the strike sectionL1. A′ is the actual boundary of the dip direction, B′ is the calculated boundary of the dip direction, C′ is the actual boundary of the rise direction, D' is the calculated boundary of the rise direction, and the distance between A′ and D′ is the length of the incline section L2. The final result is the rectangle AA′-DA′-DD′-AD′ in the vertical projection of the surface.


#### 3.1.4. Accuracy Assessment

Next, we can compare the predicted boundary to the simulated boundary. Figure8 shows the difference between the simulated and predicted boundary. The predicted boundary (red rectangle) and the simulated boundary (cyan rectangle) almost overlap, and the observed degree of coincidence is high. This indicates that the boundary predicted by the proposed method is reliable.Figure 8
Comparison of prediction results. The red rectangle represents the predicted boundary and the cyan rectangle represents the simulated boundary.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 Table1), as well as the relative error between the predicted and simulated parameters. The relative error K is calculated as follows:(8)K=fGsim,where f is the deviation between the predicted and simulated parameters, and Gsim is the simulated geometrical parameter of the goaf.Table 1
Comparison between the estimated and simulated parameters of the underground goaf.
ParametersL1 (m)L2 (m)H0 (m)H1 (m)H2 (m)α (°)Simulation value240138.5203.522018718Predictive value242.5138.3197219.6176.218.3Difference−2.5+0.2+6.5+0.4+10.8−0.3Relative error (%)1.01.43.20.25.81.7The results of the comparison are reported in Table1. The predicted geometric parameters are generally consistent with the simulated geometric parameters, with an average relative error of 2.2%. This again demonstrates the ability of the proposed method to accurately predict the boundary range of the underground goaf.


### 3.2. Real Data Experiment

#### 3.2.1. Study Area

The Pangzhuang Coal mine is located in the Jiuli District of Xuzhou City, 13 km away from the center of Xuzhou City. The surface is mainly covered by vegetation and buildings (Figure9). 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. The working face adopts the fully mechanized long arm mining process and completely manages the roof. The mining area began to operate in the 1980s and has a sufficient amount of geological and measured data of the underground goaf (Table 2). Thus, this working face possesses the necessary conditions required for the accuracy verification of the proposed method.Figure 9
Geographic location of the Pangzhuang Coal Mine.Table 2
Geological and geometric parameters of No. 7503 working face of the Pangzhuang Coal Mine.
ParametersValueParametersValueWorking face7503Strike length (m)716Subsidence coefficient (q)0.85Incline length (m)188Main influence angle tangent (tanβ)2.71Thickness (m)6.87Inflection point offset (kL)0.13Dip angle (°)4Propagation coefficient (K)0.8Average depth (m)315.5


#### 3.2.2. 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. The wavelength of the C-band is 5.6 cm, the incident angle is 22.78°, and the ground resolution is 20.12 m. The image coverage area is shown in Figure10. Due to the time-space correlation, only five images generated interference patterns, taken on the following dates: 2009-12-01, 2010-01-5, 2010-02-09, 2010-03-16, and 2010-04-20. The parameters of the interference pairs are reported in Table 3.Figure 10
Footprint of the ASAR imagery used in this experiment.Table 3
Parameters of the chosen interferometric pairs.
No.Master imageSlave imageTemporal baseline (day)Normal baseline (m)Incidence angle of master image12009-12-012010-01-535−31.8622.816 422010-01-52010-02-0935360.8622.802 032010-02-092010-03-1635−333.9522.795 142010-03-162010-04-2035230.6222.817 2ENVI 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. The 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. The GCPs (ground control points) are selected based on their suitability for most unwrapped phases. These 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. The 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 Figure11.Figure 11
DInSAR-derived LOS deformation maps over four different time periods.From Figure11, it can be seen that several deformation abnormalities can be observed. But only the area of 7503 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. Therefore, we can attempt to estimate the boundary extent of an underground goaf in the area using the proposed method.


#### 3.2.3. Experimental Results

The method proposed in this paper aims to estimate the boundary of the underground goaf in the study area. The time series surface subsidence can be obtained by stacking the images (as in Figure12). 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. Then allow for the generation of the subsidence contour map after the tool finishing running. The 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).Figure 12
DInSAR-derived LOS deformation maps.Figure 13
Contour generation, the extraction of two main sections, and the plotting of the subsidence curves.The 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 ofH0 = 305.6 m, the dip direction depth of H1 = 318.9 m, and the rise direction depth of H2 = 300.2 m. Then 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., S0=0.13×H0=39.7m, S1=0.13×H1=41.5m, S2=0.13×H2=39.0m). The inflection point of the strike section is used to determine the calculated boundary of the strike section, and S0 is used to determine the mining boundary of the strike section. The distance between these two boundary points is the strike section length L1 = 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°. Then 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. The mining boundary of the incline section is obtained by shifting the calculated boundary outwards by S1 and S2. The distance between these two boundary points is the strike section length L2 = 176.9. Finally, the estimated results of the boundaries of the underground goaf can be obtained. The boundaries of the probe are marked with a red rectangle in Figure 14.Figure 14
Estimate result by DInSAR. The red cuboids represent the estimated boundary of the underground goaf. Points A, D, A′ and D′ are the actual boundaries. Points B, C, B′ and C′ are the calculated boundaries. The distance between A and D is the length of the strike sectionL1. And the distance between A′ and D′ is the width of the incline section L2.


#### 3.2.4. 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 Figure15. The 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), as well as the relative error between the predicted and actual parameters (see Table 4).Figure 15
Differences between the boundary of the goaf predicted by DInSAR and standard geophysical techniques. The red rectangle represents the boundary obtained by the proposed method and the blue rectangle represents the boundary obtained by the geophysical techniques.Table 4
Comparison between the estimated and actual geometric parameters of the underground goaf.
ParametersL1 (m)L2 (m)H0 (m)α (°)Actual value716188315.54Predictive value703.4176.9305.64.3Deviation+12.6+11.1+9.9−0.3Relative error (%)1.75.93.17.5Relative error in Yang (2018) (%)12.46.58.820As reported in Table4, the relative errors of the four predicted geometric parameters are below 10%, with an average value of 4.55%. The 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 L1, L2, H0, and α are reduced by 7.375%, suggesting that the method proposed in this paper is highly accurate for estimating two-dimensional boundaries.


## 3.1. Simulation Experiment

Assume that a coal face has the following geometric parameters: coal thicknessm = 2.2 m, dip angle α = 18°, strike length L1 = 240 m, inclined length L2 = 138.5 m, dip section boundary depth H1 = 220 m, and rise direction boundary depth H2 = 187 m. Moreover, the parameters of the probability integral method are as follows: subsidence coefficient q = 0.65, tan β = 2.0, inflection point offset S = 0.15 H, and mining reflection coefficient kL = 0.8. Matlab (2015b, MathWorks) is used to apply the probability integral method to simulate the vertical displacement of the surface as a result of mining. After adding a random error to the simulated value, the subsidence contour of the area is obtained.In order to estimate the boundary range of the goaf, the direction of the two main sections is determined according to the subsidence contour. The 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.


### 3.1.1. 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. The subsidence values, 0.84Wm and 0.16Wm, are extracted from the subsidence curve. The distance between these two points is 0.8r0, which is subsequently used to calculate r0. According to the probability integral method parameter tan β, the strike depth is determined as H0 = r0·tan β = 197 m.The inflection offset point of the simulated working face is determined asS0 = 0.15H0 = 29.55 m. Thus, the actual boundary of the simulated working face can be determined following an outward shift of the calculated boundary by S0.


### 3.1.2. Incline Section Boundary Estimate

The subsidence values of 0.84Wm and 0.16Wm are extracted from the subsidence curve. The 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 r1 and r2, respectively. Following this, tan β is used to calculate the dip direction depth of H1 = 219.6 m and the rise direction depth of H2 = 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. The vertical line of the rise direction intersects the rise direction depth, while the other line intersects with the dip direction depth. These 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 kL, 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 directionS1 and the rise direction S2 allows for the derivation of the actual boundary of the incline section working face after the calculated boundary is shifted outward by S1 and S2, respectively.


### 3.1.3. Experimental Results

Using the estimated geometric parameters of the goaf, the mining boundary can be determined. The final estimate results are shown in Figure7.Figure 7
Estimate results. A and D are the actual boundaries of the strike section, C and B are the calculated boundaries of the strike section, and the distance between A and D is the length of the strike sectionL1. A′ is the actual boundary of the dip direction, B′ is the calculated boundary of the dip direction, C′ is the actual boundary of the rise direction, D' is the calculated boundary of the rise direction, and the distance between A′ and D′ is the length of the incline section L2. The final result is the rectangle AA′-DA′-DD′-AD′ in the vertical projection of the surface.


### 3.1.4. Accuracy Assessment

Next, we can compare the predicted boundary to the simulated boundary. Figure8 shows the difference between the simulated and predicted boundary. The predicted boundary (red rectangle) and the simulated boundary (cyan rectangle) almost overlap, and the observed degree of coincidence is high. This indicates that the boundary predicted by the proposed method is reliable.Figure 8
Comparison of prediction results. The red rectangle represents the predicted boundary and the cyan rectangle represents the simulated boundary.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 Table1), as well as the relative error between the predicted and simulated parameters. The relative error K is calculated as follows:(8)K=fGsim,where f is the deviation between the predicted and simulated parameters, and Gsim is the simulated geometrical parameter of the goaf.Table 1
Comparison between the estimated and simulated parameters of the underground goaf.
ParametersL1 (m)L2 (m)H0 (m)H1 (m)H2 (m)α (°)Simulation value240138.5203.522018718Predictive value242.5138.3197219.6176.218.3Difference−2.5+0.2+6.5+0.4+10.8−0.3Relative error (%)1.01.43.20.25.81.7The results of the comparison are reported in Table1. The predicted geometric parameters are generally consistent with the simulated geometric parameters, with an average relative error of 2.2%. This again demonstrates the ability of the proposed method to accurately predict the boundary range of the underground goaf.


## 3.1.1. 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. The subsidence values, 0.84Wm and 0.16Wm, are extracted from the subsidence curve. The distance between these two points is 0.8r0, which is subsequently used to calculate r0. According to the probability integral method parameter tan β, the strike depth is determined as H0 = r0·tan β = 197 m.The inflection offset point of the simulated working face is determined asS0 = 0.15H0 = 29.55 m. Thus, the actual boundary of the simulated working face can be determined following an outward shift of the calculated boundary by S0.


## 3.1.2. Incline Section Boundary Estimate

The subsidence values of 0.84Wm and 0.16Wm are extracted from the subsidence curve. The 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 r1 and r2, respectively. Following this, tan β is used to calculate the dip direction depth of H1 = 219.6 m and the rise direction depth of H2 = 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. The vertical line of the rise direction intersects the rise direction depth, while the other line intersects with the dip direction depth. These 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 kL, 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 directionS1 and the rise direction S2 allows for the derivation of the actual boundary of the incline section working face after the calculated boundary is shifted outward by S1 and S2, respectively.


## 3.1.3. Experimental Results

Using the estimated geometric parameters of the goaf, the mining boundary can be determined. The final estimate results are shown in Figure7.Figure 7
Estimate results. A and D are the actual boundaries of the strike section, C and B are the calculated boundaries of the strike section, and the distance between A and D is the length of the strike sectionL1. A′ is the actual boundary of the dip direction, B′ is the calculated boundary of the dip direction, C′ is the actual boundary of the rise direction, D' is the calculated boundary of the rise direction, and the distance between A′ and D′ is the length of the incline section L2. The final result is the rectangle AA′-DA′-DD′-AD′ in the vertical projection of the surface.


## 3.1.4. Accuracy Assessment

Next, we can compare the predicted boundary to the simulated boundary. Figure8 shows the difference between the simulated and predicted boundary. The predicted boundary (red rectangle) and the simulated boundary (cyan rectangle) almost overlap, and the observed degree of coincidence is high. This indicates that the boundary predicted by the proposed method is reliable.Figure 8
Comparison of prediction results. The red rectangle represents the predicted boundary and the cyan rectangle represents the simulated boundary.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 Table1), as well as the relative error between the predicted and simulated parameters. The relative error K is calculated as follows:(8)K=fGsim,where f is the deviation between the predicted and simulated parameters, and Gsim is the simulated geometrical parameter of the goaf.Table 1
Comparison between the estimated and simulated parameters of the underground goaf.
ParametersL1 (m)L2 (m)H0 (m)H1 (m)H2 (m)α (°)Simulation value240138.5203.522018718Predictive value242.5138.3197219.6176.218.3Difference−2.5+0.2+6.5+0.4+10.8−0.3Relative error (%)1.01.43.20.25.81.7The results of the comparison are reported in Table1. The predicted geometric parameters are generally consistent with the simulated geometric parameters, with an average relative error of 2.2%. This again demonstrates the ability of the proposed method to accurately predict the boundary range of the underground goaf.


## 3.2. Real Data Experiment

### 3.2.1. Study Area

The Pangzhuang Coal mine is located in the Jiuli District of Xuzhou City, 13 km away from the center of Xuzhou City. The surface is mainly covered by vegetation and buildings (Figure9). 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. The working face adopts the fully mechanized long arm mining process and completely manages the roof. The mining area began to operate in the 1980s and has a sufficient amount of geological and measured data of the underground goaf (Table 2). Thus, this working face possesses the necessary conditions required for the accuracy verification of the proposed method.Figure 9
Geographic location of the Pangzhuang Coal Mine.Table 2
Geological and geometric parameters of No. 7503 working face of the Pangzhuang Coal Mine.
ParametersValueParametersValueWorking face7503Strike length (m)716Subsidence coefficient (q)0.85Incline length (m)188Main influence angle tangent (tanβ)2.71Thickness (m)6.87Inflection point offset (kL)0.13Dip angle (°)4Propagation coefficient (K)0.8Average depth (m)315.5


### 3.2.2. 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. The wavelength of the C-band is 5.6 cm, the incident angle is 22.78°, and the ground resolution is 20.12 m. The image coverage area is shown in Figure10. Due to the time-space correlation, only five images generated interference patterns, taken on the following dates: 2009-12-01, 2010-01-5, 2010-02-09, 2010-03-16, and 2010-04-20. The parameters of the interference pairs are reported in Table 3.Figure 10
Footprint of the ASAR imagery used in this experiment.Table 3
Parameters of the chosen interferometric pairs.
No.Master imageSlave imageTemporal baseline (day)Normal baseline (m)Incidence angle of master image12009-12-012010-01-535−31.8622.816 422010-01-52010-02-0935360.8622.802 032010-02-092010-03-1635−333.9522.795 142010-03-162010-04-2035230.6222.817 2ENVI 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. The 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. The GCPs (ground control points) are selected based on their suitability for most unwrapped phases. These 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. The 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 Figure11.Figure 11
DInSAR-derived LOS deformation maps over four different time periods.From Figure11, it can be seen that several deformation abnormalities can be observed. But only the area of 7503 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. Therefore, we can attempt to estimate the boundary extent of an underground goaf in the area using the proposed method.


### 3.2.3. Experimental Results

The method proposed in this paper aims to estimate the boundary of the underground goaf in the study area. The time series surface subsidence can be obtained by stacking the images (as in Figure12). 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. Then allow for the generation of the subsidence contour map after the tool finishing running. The 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).Figure 12
DInSAR-derived LOS deformation maps.Figure 13
Contour generation, the extraction of two main sections, and the plotting of the subsidence curves.The 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 ofH0 = 305.6 m, the dip direction depth of H1 = 318.9 m, and the rise direction depth of H2 = 300.2 m. Then 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., S0=0.13×H0=39.7m, S1=0.13×H1=41.5m, S2=0.13×H2=39.0m). The inflection point of the strike section is used to determine the calculated boundary of the strike section, and S0 is used to determine the mining boundary of the strike section. The distance between these two boundary points is the strike section length L1 = 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°. Then 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. The mining boundary of the incline section is obtained by shifting the calculated boundary outwards by S1 and S2. The distance between these two boundary points is the strike section length L2 = 176.9. Finally, the estimated results of the boundaries of the underground goaf can be obtained. The boundaries of the probe are marked with a red rectangle in Figure 14.Figure 14
Estimate result by DInSAR. The red cuboids represent the estimated boundary of the underground goaf. Points A, D, A′ and D′ are the actual boundaries. Points B, C, B′ and C′ are the calculated boundaries. The distance between A and D is the length of the strike sectionL1. And the distance between A′ and D′ is the width of the incline section L2.


### 3.2.4. 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 Figure15. The 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), as well as the relative error between the predicted and actual parameters (see Table 4).Figure 15
Differences between the boundary of the goaf predicted by DInSAR and standard geophysical techniques. The red rectangle represents the boundary obtained by the proposed method and the blue rectangle represents the boundary obtained by the geophysical techniques.Table 4
Comparison between the estimated and actual geometric parameters of the underground goaf.
ParametersL1 (m)L2 (m)H0 (m)α (°)Actual value716188315.54Predictive value703.4176.9305.64.3Deviation+12.6+11.1+9.9−0.3Relative error (%)1.75.93.17.5Relative error in Yang (2018) (%)12.46.58.820As reported in Table4, the relative errors of the four predicted geometric parameters are below 10%, with an average value of 4.55%. The 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 L1, L2, H0, and α are reduced by 7.375%, suggesting that the method proposed in this paper is highly accurate for estimating two-dimensional boundaries.


## 3.2.1. Study Area

The Pangzhuang Coal mine is located in the Jiuli District of Xuzhou City, 13 km away from the center of Xuzhou City. The surface is mainly covered by vegetation and buildings (Figure9). 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. The working face adopts the fully mechanized long arm mining process and completely manages the roof. The mining area began to operate in the 1980s and has a sufficient amount of geological and measured data of the underground goaf (Table 2). Thus, this working face possesses the necessary conditions required for the accuracy verification of the proposed method.Figure 9
Geographic location of the Pangzhuang Coal Mine.Table 2
Geological and geometric parameters of No. 7503 working face of the Pangzhuang Coal Mine.
ParametersValueParametersValueWorking face7503Strike length (m)716Subsidence coefficient (q)0.85Incline length (m)188Main influence angle tangent (tanβ)2.71Thickness (m)6.87Inflection point offset (kL)0.13Dip angle (°)4Propagation coefficient (K)0.8Average depth (m)315.5


## 3.2.2. 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. The wavelength of the C-band is 5.6 cm, the incident angle is 22.78°, and the ground resolution is 20.12 m. The image coverage area is shown in Figure10. Due to the time-space correlation, only five images generated interference patterns, taken on the following dates: 2009-12-01, 2010-01-5, 2010-02-09, 2010-03-16, and 2010-04-20. The parameters of the interference pairs are reported in Table 3.Figure 10
Footprint of the ASAR imagery used in this experiment.Table 3
Parameters of the chosen interferometric pairs.
No.Master imageSlave imageTemporal baseline (day)Normal baseline (m)Incidence angle of master image12009-12-012010-01-535−31.8622.816 422010-01-52010-02-0935360.8622.802 032010-02-092010-03-1635−333.9522.795 142010-03-162010-04-2035230.6222.817 2ENVI 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. The 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. The GCPs (ground control points) are selected based on their suitability for most unwrapped phases. These 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. The 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 Figure11.Figure 11
DInSAR-derived LOS deformation maps over four different time periods.From Figure11, it can be seen that several deformation abnormalities can be observed. But only the area of 7503 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. Therefore, we can attempt to estimate the boundary extent of an underground goaf in the area using the proposed method.


## 3.2.3. Experimental Results

The method proposed in this paper aims to estimate the boundary of the underground goaf in the study area. The time series surface subsidence can be obtained by stacking the images (as in Figure12). 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. Then allow for the generation of the subsidence contour map after the tool finishing running. The 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).Figure 12
DInSAR-derived LOS deformation maps.Figure 13
Contour generation, the extraction of two main sections, and the plotting of the subsidence curves.The 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 ofH0 = 305.6 m, the dip direction depth of H1 = 318.9 m, and the rise direction depth of H2 = 300.2 m. Then 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., S0=0.13×H0=39.7m, S1=0.13×H1=41.5m, S2=0.13×H2=39.0m). The inflection point of the strike section is used to determine the calculated boundary of the strike section, and S0 is used to determine the mining boundary of the strike section. The distance between these two boundary points is the strike section length L1 = 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°. Then 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. The mining boundary of the incline section is obtained by shifting the calculated boundary outwards by S1 and S2. The distance between these two boundary points is the strike section length L2 = 176.9. Finally, the estimated results of the boundaries of the underground goaf can be obtained. The boundaries of the probe are marked with a red rectangle in Figure 14.Figure 14
Estimate result by DInSAR. The red cuboids represent the estimated boundary of the underground goaf. Points A, D, A′ and D′ are the actual boundaries. Points B, C, B′ and C′ are the calculated boundaries. The distance between A and D is the length of the strike sectionL1. And the distance between A′ and D′ is the width of the incline section L2.


## 3.2.4. 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 Figure15. The 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), as well as the relative error between the predicted and actual parameters (see Table 4).Figure 15
Differences between the boundary of the goaf predicted by DInSAR and standard geophysical techniques. The red rectangle represents the boundary obtained by the proposed method and the blue rectangle represents the boundary obtained by the geophysical techniques.Table 4
Comparison between the estimated and actual geometric parameters of the underground goaf.
ParametersL1 (m)L2 (m)H0 (m)α (°)Actual value716188315.54Predictive value703.4176.9305.64.3Deviation+12.6+11.1+9.9−0.3Relative error (%)1.75.93.17.5Relative error in Yang (2018) (%)12.46.58.820As reported in Table4, the relative errors of the four predicted geometric parameters are below 10%, with an average value of 4.55%. The 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 L1, L2, H0, and α are reduced by 7.375%, suggesting that the method proposed in this paper is highly accurate for estimating two-dimensional boundaries.


## 4. Conclusion

This paper reports an approach for estimating the boundary range of underground goaf using the DInSAR by combining a graphical method. The 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. The 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 lengthL1, the incline section length L2, and the depths H0, H1, and H2 (predicted by the measured data) are all underestimated. This 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.



---

*Source: 9375056-2020-05-29.xml*