A Goaf-Locating Method Based on the D-InSAR Technique and Stratified Okada Dislocation Model

Abstract

Illegal coal mining is prevalent worldwide, leading to extensive ground subsidence and land collapse. It is crucial to define the location and spatial dimensions of these areas for the efficient prevention of the induced hazards. Conventional methods for goaf locating using the InSAR technique are mostly based on the probability integral model (PIM). However, The PIM requires detailed mining information to preset model parameters and does not account for the layered structure of the coal overburden, making it challenging to detect underground goaves in cases of illegal mining. In response, a novel method based on the InSAR technique and the Stratified Optimal Okada Dislocation Model, named S-ODM, is proposed for locating goaves with basic geological information. Firstly, the S-ODM employs a numerical model to establish a nonlinear function between the goaf parameters and InSAR-derived ground deformation. Then, in order to mitigate the influence of nearby mining activities, the goaf azimuth angle is estimated using the textures and trends of the InSAR-derived deformation time series. Finally, the goaf’s dimensions and location are estimated by the genetic algorithm–particle swarm optimization (GA-PSO). The effectiveness of the proposed method is validated using both simulation and real data, demonstrating average relative errors of 6.29% and 7.37%, respectively. Compared with the PIM and ODM, the proposed S-ODM shows improvements of 19.48% and 52.46% in geometric parameters. Additionally, the errors introduced by GA-PSO and the influence of ground deformation monitoring errors are discussed in this study.

1. Introduction

Total global coal consumption exceeded 8.5 billion tons in 2023, with China accounting for over 50% of the total demand [1]. With the rapid growth in the demand for coal, illegal mining activities have significantly increased, leading to extensive ground subsidence and land collapse [2,3,4]. Besides the induced geological hazards, these activities threaten the interests of legal miners and the lives of local residents, as well as those of illegal miners. According to relevant studies, officially unrecorded mining activities account for 56% of all mining worldwide, highlighting the extreme challenges national governments face in monitoring and managing illegal coal mining [5]. Additionally, with urbanization increasing exponentially, issues related to land resource management have become more prominent, such as housing shortages, land scarcity, and traffic congestion [6]. Therefore, it is crucial to utilize underground space efficiently and scientifically, including managing underground goaves [7]. To this end, detecting and locating underground goaves, especially unrecorded ones, is extremely essential.

Conventional methods for detecting underground goaves are mostly based on geophysical and chemical techniques, including microgravity techniques [8], electromagnetic methods [9], ground-penetrating radar [9], rare gas measurements, and thermal techniques [10,11,12]. These methods can obtain accurate information on both the location and internal structure of underground goaves in small-scale areas. However, they have significant limitations, including inefficient small-scale monitoring, a time-consuming and labor-intensive nature, and dependence on prior geological information. In recent years, numerous unrecorded underground goaves have been found widely distributed across mining industrial areas. Consequently, efficiently defining and locating these spaces over such large-scale areas using geophysical and chemical techniques poses significant challenges in practice.

Compared to conventional geophysical and chemical methods, satellite-based interferometric synthetic aperture radar (InSAR) technology can monitor large areas cost-effectively and efficiently [13]. Based on this, InSAR has emerged as a powerful tool for monitoring various geohazards, including earthquakes [14], volcanoes [15], landslides [16], coal mining [17,18], and coal fires [19,20].

Underground goaf-locating methods can be classified into two main categories. The first category uses key features of ground deformation, such as inflection points and boundaries, to locate underground goaves according to the constructed geometrical relations between coal mining and induced ground subsidence. Hu et al. proposed a D-InSAR-based illegal mining detection system (DIMDS) to locate underground goaves using the empirical spatial characteristics of mining-related deformation [21]. However, this approach can only obtain the central geodetic coordinates of the mine area but cannot obtain the underground geometric parameters. To this end, Du et al. introduced the Feature Point Model (FPM) to estimate the extent of underground goaves. This type of method is computationally efficient but it cannot be universally applied since the geometrical relationships are steady only under supercritical extraction conditions [22]. The second category is based on theoretical models. Yang et al. proposed a model describing the relationship between goaves and ground deformation based on the probability integral model (PIM), employing nonlinear search algorithms for goaf detection (abbreviated as LGM-PIM) [23]. Building upon this approach, Wang et al. proposed the improved probability integral model (IPIM) to extend its applicability to cases involving subcritical underground goaves [24]. Additionally, some studies have applied the Okada Dislocation Model (ODM) to establish nonlinear relationships between underground goaf parameters and ground deformation, achieving notable results [25]. Compared to the first category, these theoretical approaches offer a more robust foundation and are capable of achieving higher accuracy in locating underground goaves.

However, these two theoretical models both have significant limitations. Specifically, the PIM requires detailed mining information to preset model parameters, while the ODM simplifies the coal overburden as a global unit [26,27]. It is extremely challenging to locate underground goaves using these two traditional models in cases of illegal mining without detailed mining information.

To address the issue, a method based on InSAR-derived ground deformation and the Stratified Okada Dislocation Model (hereafter referred to as S-ODM) is proposed for locating underground goaves with basic geological information. The proposed S-ODM can better fit the relationship between an underground goaf’s parameters and ground deformation compared to the conventional ODM and requires less prior information to preset model parameters than the PIM employed in previous studies. In addition, the goaf azimuth is estimated using the ray method based on the textures and trend of the InSAR-derived deformation time series to reduce the influence of the surrounding mining, thereby improving computational efficiency. Finally, the genetic algorithm–particle swarm optimization (GA-PSO) is applied to estimate the optimal values of other characteristic parameters.

This paper is organized as follows. Section 2 explains the study area and the data used for this study. Section 3 gives the details of the proposed method. Simulation and real data experiments and results are presented in Section 4. The discussions and conclusions are presented in Section 5 and Section 6, respectively.

2. Study Area and Datasets

The study area is the 15325 working face of Fengfeng Coalfield, situated in the southern region of Hebei Province, China. It has a length of 493 m in the strike direction, a length of 142 m in the dip direction, an average coal mining thickness of 4.5 m, and an average mining depth of 740 m. The inclination of the working face is 13°, and the goaf azimuth is 236°.

Typically, when the inclination of the coal seam is less than 15°, the surface deformation from mining is akin to that from horizontal coal seam mining. However, the S-ODM applied in this study is sensitive to the inclination. Consequently, the coal seam mining in this study is not considered horizontal coal seam mining. As depicted in Figure 1, numerous faults and geological structures surround the study area’s working face. In addition, extensive and prolonged coal mining activities have resulted in several old goaves in the vicinity. These factors significantly increase the challenges associated with locating underground goaves.

As depicted in Table 1, six RADARSAT-2 images from 13 October 2015 to 5 March 2016, covering the coal mining period, were utilized for capturing mining-induced deformation. Additionally, SRTM (Shuttle Radar Topography Mission) data with a 30 m resolution served as the external DEM (digital elevation model) to mitigate topographic phase effects. The atmospheric errors were negligible due to the small scope of the study area. Given the precise measurement and control of the RADARSAT-2 satellite’s attitude, orbit errors were also negligible. The Goldstein filter was then applied to suppress noise phases. Subsequently, the minimum-cost flow approach was employed to unwrap the differential interferograms, and the deformation phases were obtained [28]. Finally, the deformation phases were geocoded into the WGS-84 geographic coordinate system.

3. Methods

Figure 2 depicts the workflow of the proposed method. The proposed algorithm encompasses four primary steps:

• Acquire times-series line-of-sight (LOS) deformation utilizing the D-InSAR technique;

• Determine the S-ODM parameters based on basic geological investigations;

• Apply the ray method to estimate the goaf azimuth with the textures and trends of the InSAR-derived deformation;

• Use GA-PSO to estimate the goaf parameters. The process will stop when either the best fitness value is less than the threshold or the iteration number reaches a preset maximum value.

In the proposed method, the InSAR-derived deformation serves as the reference value. The simulated deformation generated by the S-ODM iteratively approaches this reference value, leading to the derivation of the optimal parameters.

3.1. Differential InSAR

In this study, differential InSAR (D-InSAR) technology is applied to obtain time-series ground deformation. Satellite-based SAR repeatedly acquires single-look complex (SLC) images of the same region at regular intervals [29]. An interferogram is formed by the conjugate multiplication of two registered SAR images. The interferometric phase $\phi$ is mainly composed of the following parts [30]:

(1)$\phi ={\phi }_{def}+{\phi }_{top}+{\phi }_{flt}+{\phi }_{atm}+{\phi }_{noi}$

where ${\phi }_{def}$ is the deformation phase in the LOS direction, ${\phi }_{top}$ is the topographic phase, ${\phi }_{flt}$ is the flat phase, ${\phi }_{atm}$ is the atmospheric phase, and ${\phi }_{noi}$ is the noise phase. In the D-InSAR technique, ${\phi }_{def}$ can be obtained by removing and suppressing the last five phase terms.

3.2. Stratified Okada Dislocation Model

The ODM has been widely utilized in studies of volcanoes and earthquakes [31]. To enhance computational efficiency, the model simplifies the strata as an elastic semi-infinite space confined by the ground surface [32]. Generally, subsidence of the underground goaf can be considered as a rectangular-source dislocation. Thus, the improved model was based on the rectangular-source Okada Dislocation Model. In addition, using the ODM to calculate the goaf-induced ground deformation involves three coordinate systems. Their spatial geometrical relationships are shown in Figure 3.

The coordinate system $(0-xyz)$ is the localized Cartesian coordinate system for the study area, where the x-axis represents the north direction, the y-axis represents the east direction, and the z-axis is upward perpendicular to the xy-plane. The coordinate system $(\widehat{o}-\widehat{x}\widehat{y}\widehat{z})$ is the research object’s ground coordinate system, where the $\widehat{x}$-axis represents the goaf strike direction, the $\widehat{y}$-axis represents the goaf dip direction, and the $\widehat{z}$-axis represents the vertical direction. Specifically, $\phi$ is the angle between the $\widehat{x}$-axis and the x-axis, representing the goaf azimuth angle, and the origin is $(x_0,y_0)$ in $(0-xyz)$. The coordinate system $o_1-\xi \widehat{p}\widehat{q}$ is the research object’s underground coordinate system, where the $\xi$-axis represents the strike direction of the working face, the $\widehat{p}$-axis represents the dip direction of the working face, and the $\widehat{q}$-axis is upward perpendicular to the working face. In addition, $\delta$ is the dip angle of the working face.

Assume that there is a point P on the ground with coordinates $(x,y,0)$ in $(0-xyz)$ and $(\widehat{x},\widehat{y},0)$ in $(\widehat{o}-\widehat{x}\widehat{y}\widehat{z})$, which can be converted according to the following formula:

(2)$\left(\widehat{x},\widehat{y},\widehat{z}\right)=\left(x-x_0,y-y_0,0\right)\left[\begin{array}{ccc}cos\phi & sins\phi & 0\\ sin\phi & -cos\phi & 0\\ 0& 0& 1\end{array}\right]$

In the research object’s underground coordinate system $o_1-\xi \widehat{p}\widehat{q}$, the origin represents the bottom-left corner of the working face with a depth of d. Point P can be converted to the underground coordinate system using the following formula.

(3)$\left(\xi ,p,q\right)=\left(\widehat{x},\widehat{y},d\right)\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\theta & -sin\theta \\ 0& sin\phi \theta & cos\theta \end{array}\right]$

Subsequently, the deformation of point P induced by coal mining can be represented by the following formula, where $d_e$, $d_n$, and $d_u$ represent the deformation in the horizontal and vertical directions [32].

(4)$\left\{\begin{array}{c}{d}_e=f_x\left(p,q\right)-f_x\left(p,q-W\right)-f_x\left(p-L,q\right)-f_x\left(p-L,q-W\right)\\ d_n=f_y\left(p,q\right)-f_y\left(p,q-W\right)-f_y\left(p-L,q\right)-f_y\left(p-L,q-W\right)\\ d_u=f_z\left(p,q\right)-f_z\left(p,q-W\right)-f_z\left(p-L,q\right)-f_z\left(p-L,q-W\right)\end{array}\right.$

with

(5)$f_x(\widehat{x},\widehat{y})={\displaystyle \frac{U}{2\pi }}\left[{\displaystyle \frac{q^2}{R(R+\widehat{y})}}-I_{2}sin{\delta }^2\right]$

(6)$f_y(\widehat{x},\widehat{y})={\displaystyle \frac{U}{2\pi }}\left[{\displaystyle \frac{-\left(\widehat{y}sin\delta -qcos\delta \right)q}{R(R+\widehat{x})}}-sin\delta \left({\displaystyle \frac{\widehat{x}q}{R(R+\widehat{y})}}-{tan}^{-1}{\displaystyle \frac{\widehat{x}\widehat{y}}{qR}}\right)\right]-I_{1}sin{\delta }^2$

(7)$f_z(\widehat{x},\widehat{y})={\displaystyle \frac{U}{2\pi }}\left[{\displaystyle \frac{-\left(\widehat{y}sin\delta -qcos\delta \right)q}{R(R+\widehat{x})}}+cos\delta \left({\displaystyle \frac{\widehat{x}q}{R(R+\widehat{y})}}-{tan}^{-1}{\displaystyle \frac{\widehat{x}\widehat{y}}{qR}}\right)\right]-I_{4}sin{\delta }^2$

(8)$I_1={\displaystyle \frac{\mu }{\lambda +\mu }}\left[{\displaystyle \frac{-1}{cos\delta }}{\displaystyle \frac{\xi }{R+(\eta sin\delta -qcos\delta )}}\right]-{\displaystyle \frac{sin\delta }{cos\delta }}{I}_4$

(9)$I_2={\displaystyle \frac{\mu }{\lambda +\mu }}\left[{\displaystyle \frac{-1}{cos\delta }}{\displaystyle \frac{\widehat{y}cos\delta +qsin\delta }{R+(\widehat{y}sin\delta -qcos\delta )}}-ln(R+\widehat{y})\right]+{\displaystyle \frac{sin\delta }{cos\delta }}{I}_3$

(10)$I_3={\displaystyle \frac{\mu }{\lambda +\mu }}{\displaystyle \frac{1}{cos\delta }}\left[ln(R+(\widehat{y}sin\delta -qcos\delta ))-sin\delta ln(R+\widehat{y})\right]$

(11)$I_4={\displaystyle \frac{\mu }{\lambda +\mu }}{\displaystyle \frac{2}{cos\delta }}{tan}^{-1}{\displaystyle \frac{\widehat{y}(X-qcos\delta )+X(R+X)sin\delta }{\widehat{x}(R+X)cos\delta }}$

(12)$p=ycos\delta +dsin\delta$

(13)$q=-ysin\delta +dcos\delta$

(14)$R=\sqrt{{\xi }^2+p^2+q^2}$

(15)$X=\sqrt{{\xi }^2+q^2}$

(16)$\mu ={\displaystyle \frac{E}{2(1+\nu )}}$

(17)$\lambda ={\displaystyle \frac{E\nu }{(1-2\nu )(1+\nu )}}$

where E and $\nu$ are the elastic moduli and Poisson’s ratio of the overlying stratum and $\delta$ is the dip angle of the surface source. Specifically, U is the dislocation in the vertical direction, and in the case of coal mining, it represents the subsidence of the coal seam roof. Generally, this value differs from the thickness h of the coal seam, and the relationship between them can be expressed as follows:

(18)$U=h\ast S$

where S is the subsidence factor, which is determined by the geological conditions and mining techniques. In order to improve the performance of the model in the layered structure, the Stratified Okada Dislocation Model (S-ODM) is proposed. Assuming that there are i layers of geological structure from the mining face to the ground, the deformation of the $i\mathrm{t}\mathrm{h}$ layer can be mathematically expressed as follows:

(19)$D_i=\left\{\begin{array}{c}f\left(\xi ,\eta ,h,S_1,{\nu }_1\right),i=1\hfill \\ f\left(\xi ,\eta ,D_{i-1},S_i,{\nu }_i\right),i>1\hfill \end{array}\right.$

As depicted in Figure 4, the transfer of the deformation considers the coal overburden stratum as a layered structure rather than a simple global unit. It is evident that the S-ODM is more appropriate for real conditions.

3.3. Underground Goaf Locating

3.3.1. Mining Azimuth Estimation

Most coal mining activities are surrounded by other mining working faces. Consequently, the cumulative deformation tends to deviate toward old goaves. Under these conditions, using the cumulative ground deformation to estimate the azimuth angle can result in significant errors. To this end, an azimuth angle detection method using the ray method based on the deformation binary conversion of time-series InSAR-derived deformation, named RM-DBC, is proposed. The main process is shown in Figure 5.

First, the LOS deformation images are converted into binary images by setting the deformation threshold, which is determined by the average value of a stable region without subsidence. Then, the binary image is scanned in a round using a ray at a step of 1°. The sum of the pixel values intersected by each ray is recorded. Finally, the direction with the highest sum of pixel values is regarded as the azimuth angle. Additionally, it should be pointed out that the origin point of the ray method is referred to as the maximum subsidence point.

3.3.2. Estimation of Goaf Parameters

The $LOS$ direction deformation, observed using the InSAR technique, is composed of three-dimensional surface deformations $(W,U_E,U_N)$. The model can be formulated as follows:

(20)$D_{LOS}=Wcos{\theta }_s-sin{\theta }_s\left(U_{N}cos{\alpha }_s+U_{E}sin{\alpha }_s\right)$

where ${\theta }_s$ and ${\alpha }_s$ represent the incidence angle and the heading angle of the SAR satellite. The S-ODM establishes a relationship between underground deformation and ground subsidence. The theoretical equations are derived by combining the InSAR-derived deformation and the S-ODM:

(21)$D_{LOS}=Wcos{\theta }_s-sin{\theta }_s\left(f_N(\xi ,\eta )cos{\alpha }_s+f_E(\xi ,\eta )sin{\alpha }_s\right)$

The equation contains nine parameters: central geodetic coordinates $(X,Y)$, strike length L, dip length W, depth H, underground goaf height h, dip angle $\delta$, and azimuthal angle $\alpha$. Due to the highly nonlinear relationship between these parameters and ground deformation, the genetic algorithm–particle swarm optimization (GA-PSO) is employed for rapid computation and global optimal parameter search. This algorithm inherits the strong global search performance of the GA and the rapid convergence efficiency of PSO. Compared to non-fused algorithms, GA-PSO exhibits a more efficient optimization process and higher accuracy [33].

During processing, GA-PSO searches for the optimal parameters, minimizing the difference between the ODM estimates and the InSAR-derived deformation. The difference is defined as the fitness value:

(22)$fit=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^n{\left(D_{{LOS}_i}-D_{{los}_i}\right)}^2}{n}}}$

where fit denotes the particle fitness value used to evaluate the quality of particles. $D_{{LOS}_i}$ is the LOS deformation of point j observed using the InSAR technique, and $D_{{los}_i}$ is the simulated LOS deformation of point j calculated by the ODM. Finally, the particle with the minimum fitness value is selected as the final parameter. The main process of GA-PSO is shown in Figure 6.

4. Results

4.1. Simulation Results

To examine the correlation between ground subsidence and the dimensions of the coal mining source, the Fast Lagrangian Analysis of Continua (FLAC3D) was employed. FLAC3D utilizes an explicit Lagrangian algorithm and a hybrid-discrete partition technique for simulating the deformation transfer. Additionally, to assess the accuracy of the S-ODM under simulated conditions, the simulation results were applied to estimate the geometrical parameters of the underground space based on the S-ODM.

In this simulation experiment, an inclined working face was simulated with a mining depth (H) of 500 m, a strike length (L) of 600 m, a dip width (W) of 100 m, an azimuth ($\alpha$) of 45°, a dip angle $\phi$ of 30°, and a mining thickness (m) of 5 m. The simulation results are shown in Figure 7. Additionally, in order to synthesize the deformation in the LOS direction, it is assumed that the incident angle and heading angle of the satellite are 35.5° and 349.6°, respectively.

Then, the LOS deformation was considered as the true value, and the initial inversion ranges of the goaf parameters were established as L ∈ [0 m,1000 m], W ∈ [0 m,1000 m], H ∈ [0 m,1000 m], $X\in [0,2000]$, $Y\in [0,2000]$, h ∈ [0 m,20 m], $\psi \in [0,{90}^{\circ }]$, and $\alpha \in [0,{90}^{\circ }]$. Next, 50 calculations of 1000 iterations each were performed, and the results are presented in Table 2.

In Table 2, it can be observed that the length parameters (L, W, H, h) exhibit favorable alignment with the simulated values, boasting an average relative error of 6.29%. Notably, the estimation of the dip length has the highest relative error of 9.64%, while the strike length records the maximum absolute error of 27.72 m. Furthermore, the absolute errors of the center coordinates amount to 8.61 m. In conclusion, these results demonstrate the feasibility of underground space localization using the S-ODM and GA-PSO.

4.2. Real Data Results

4.2.1. Deformation Results

The deformation results obtained by the processing described in Section 2 are presented in Figure 8.

It can be seen that the deformation expands toward the southwest direction, reaching a maximum subsidence value of −267 mm. In order to verify the accuracy of the InSAR-derived deformation, 50 leveling points were compared with the InSAR-derived deformation for the same period. As illustrated in Figure 9, it is evident that the InSAR monitoring results agree well with the leveling data. Specifically, the maximum deviation is only 39.1 mm (located at point 41), and the root-mean-square error (RMSE) between the InSAR and leveling measurements is 14.167 mm. This confirms the reliability of the InSAR-derived deformation in the study area.

4.2.2. Goaf Azimuth Estimate

A binary filter with a threshold of $2.6$ mm was constructed. This threshold was determined based on a stable region identified within the time-series cumulative InSAR-derived deformation, where the average deformation was $2.6$ mm. Subsequently, the deformation results were filtered using this threshold, as depicted in Figure 10.

Subsequently, the ray method was applied to estimate the goaf method. Specifically, the maximum subsidence point was set as the origin point. The results are shown in Figure 11.

In Figure 11, the maximum value of the ray method is 422, occurring at a direction of 234°. Consequently, 234° was identified as the azimuth indicating the fastest development trend of the goaf.

4.2.3. Underground Goaf Locating

According to the results in Section 4.2.2, the goaf azimuth was determined to be 234°. Additionally, some geological parameters required by the S-ODM should be defined before the inversion process. Based on prior geological investigations, the coal overburden of the study area can be classified into three main layered structures. According to related geological research [34], the first layer is soft topsoil with a thickness of 20 m and a Poisson’s ratio of 0.35, the second layer is a mixed stratum of sandy mudstone and siltstone with a thickness of 400 m and a Poisson’s ratio of 0.22, and the third layer is a mixed stratum of sandstone and sandy mudstone with a thickness of 320 m and a Poisson’s ratio of 0.18. After applying the method described in Section 4.2.3, the inversion results and errors were calculated and are shown in Table 3 and Figure 12.

In Table 3, it can be seen that the average relative errors for geometric parameters are 7.37%, generally remaining within 10%, except for the dip length, where the error reaches 13.98%. This indicates the effectiveness of the proposed method in locating underground goaves. It is noteworthy that the relative error for the dip length (13.98%) exceeds that for the strike length (5.92%). This is due to the shorter mining length in the dip direction and the impact of dip angle errors primarily affecting the dip length rather than the strike length. Furthermore, the absolute error of the goaf azimuth angle $\alpha$ is only 2°. These results demonstrate that the proposed RM-DBC can obtain reliable outcomes even with the influence of surrounding old goaves.

5. Discussion

5.1. Comparison with Existing Methods

The LGM-PIM and ODM are methods previously employed for locating underground goaves using D-InSAR. Therefore, a comparison was made between the proposed method and these two methods. The results are shown in Table 4.

It can be observed that compared with the PIM and ODM, the S-ODM improved accuracy by 19.482% and 52.46% for geometric parameters. This indicates that among these three methods, the proposed S-ODM performed the best. In order to further analyze the superiority of the S-ODM, the deformation profile lines in the strike and dip directions generated by the three methods (S-ODM, PIM, ODM) were produced to compare with the InSAR-derived data.

As shown in Figure 13, the average absolute errors in the strike direction for the three methods (S-ODM, PIM, and ODM) using D-InSAR were 14.1 mm, 24.1 mm, and 23.2 mm. Compared with the PIM and ODM, the S-ODM improved accuracy by 41.49% and 39.22%, respectively. In the dip direction, these errors were 18.6 mm, 35.1 mm, and 28.9 mm. Compared with the PIM and ODM, the S-ODM improved accuracy by 47.01% and 35.64%, respectively. This indicates that the proposed S-ODM can more accurately describe the relationship between ground deformation and underground goaf.

5.2. Errors Introduced by GA-PSO

To evaluate the performance of GA-PSO, a simulation experiment based on the S-ODM was conducted. Initially, the ground deformation was generated using the error-free S-ODM. Subsequently, GA-PSO was executed 50 times with 1000 iterations each, and the outcomes are presented in Table 5.

As shown in Table 5, the length parameters (L, W, H, h) exhibit a favorable alignment with the simulated values, demonstrating an average relative error of 5.77%. The dip length estimation exhibits the highest relative error of 8.64%, while the depth records the maximum absolute error of 19.51 m. Furthermore, the absolute errors of the center coordinates amount to 8.61 m. In conclusion, these results affirm the feasibility of localizing underground spaces using the S-ODM and GA-PSO.

5.3. Influence of Ground Deformation Monitoring Errors

To evaluate the robustness of the S-ODM against observational errors in the InSAR-derived deformation, 40 groups of normally distributed errors, ranging from −20 cm to 20 cm in 10 mm intervals with a standard deviation of 5 mm, were simulated. These errors were then incorporated into the LOS deformation obtained in Section 4.1. The subsequent estimation of the geometric parameters with varying errors was conducted using the S-ODM, and the results are presented in Figure 14.

As depicted in Figure 14, as the deformation error (−20 cm–20 cm) increases, the geometric parameters (i.e., L, W, H, and m) exhibit a decreasing trend since the range and maximum subsidence value of the deformation can be directly affected by the monitoring error, consequently influencing the accuracy of parameter estimation. Conversely, the azimuth shows random fluctuations around the true value with increasing errors, indicating its superior resistance to errors. However, the estimated value of the dip angle shows a maximum deviation of 14.89${(}^{\circ })$ with increasing errors, highlighting its sensitivity to monitoring errors in the S-ODM. The maximum localization error of the center coordinates reaches 60.11 m, which is acceptable within the kilometer-level range of the underground space. The trends in the X and Y coordinates are related to the value of the azimuth angle, and the observed trend is influenced by the 45${(}^{\circ })$ azimuth angle. In conclusion, low monitoring errors minimally impact the shape and position of the underground space. However, the accuracy of parameter estimation, excluding the azimuth, significantly decreases with increasing monitoring errors, which underscores the importance of high-accuracy deformation data as a prerequisite for effectively locating underground space.

6. Conclusions

In this study, a method based on D-InSAR and the Stratified Okada Dislocation Model, named S-ODM, has been proposed for locating underground goaves. Through a comparison of both simulation and real data results, several conclusions were drawn.

Specifically, the proposed S-ODM demonstrated good accuracy in both simulation and real conditions, achieving average relative errors of 6.29% and 7.37%, respectively, in underground goaf locating. Furthermore, for the proposed S-ODM, a new azimuth angle estimation method called RM-DBC was introduced, which can significantly reduce the influence of surrounding mining and improve the accuracy of goaf inversion. In addition, the deformation results generated by the S-ODM, PIM, and ODM were different from those obtained by D-InSAR. The results showed that when fitting to real InSAR-derived deformation, the proposed S-ODM outperformed both the PIM and ODM, reducing the absolute errors by 44.2% and 37.4%, respectively. It was demonstrated that the dimensions and locations of underground goaves could be effectively estimated using the proposed S-ODM. Furthermore, the S-ODM exhibited robustness under conditions of minimal ground deformation monitoring errors. However, the estimation accuracy of the dimensions and locations decreased significantly with increasing monitoring errors. Therefore, it is critical to obtain highly accurate deformation monitoring results.

Author Contributions

Conceptualization, K.Z., S.D., F.Z. and Y.W.; methodology, K.Z., S.D. and Y.W; software, K.Z.; validation, K.Z., S.D., F.Z., T.W. and Y.W.; formal analysis, K.Z., S.D. and N.Z.; Investigation, D.Z. and X.D.; writing—original draft preparation, K.Z.; writing—review and editing, K.Z., S.D., F.Z. and Y.W.; visualization, K.Z.; supervision, K.Z.; project administration, K.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Acknowledgments

We would like to thank the ESA (European Space Agency) for supplying the Sentinel-1 datasets. Some figures were prepared using the public domain GMT v6.4 software (Wessel and Smith, 1998). Optical images were provided by Google Earth. The writing and proofreading of this manuscript were greatly assisted by Wanting Bao, who graduated from Nanjing University of Aeronautics and Astronautics.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Geographic location of the study area.

Figure 2. Scheme of the proposed S-ODM. The boxes outlined with dotted lines represent parameters estimated by the proposed algorithm, which form the goaf parameters.

Figure 3. Spatial geometrical relationships of three coordinate systems. (a) Spatial relationships among the three coordinate systems. (b) Top view of the two ground-based coordinate systems.

Figure 4. The layered structure of the coal overburden. The deformation is transmitted to the ground surface through multiple strata.

Figure 5. Scheme of RM-DBC. (a) Simulation deformation example; the red rectangle is the stable area without subsidence. (b) The deformation binary conversion result.

Figure 6. The optimization process of GA-PSO. [Forumla omitted. See PDF.] represents the initial particle swarm and is constructed using seven parameters (i.e., [Forumla omitted. See PDF.], L, W, H, m, [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]). [Forumla omitted. See PDF.] (i is the iteration number, [Forumla omitted. See PDF.]) represents the [Forumla omitted. See PDF.] particles sorted by fitness value, and [Forumla omitted. See PDF.] is the particle swarm optimized by PSO and the GA for the next iteration.

Figure 7. Ground deformation of simulated mining. (a) Ground subsidence. (b,c) Horizontal deformation in the north–south and east–west directions, respectively. (d) Simulated LOS deformation. The blue and pink rectangles represent the simulated and estimated work faces.

Figure 8. Time-series cumulative subsidence map of the study area. The black line is the 15325 working face.

Figure 9. (a) Comparison between the InSAR measurements from 13 October 2015 to 5 March 2016 and the leveling measurements from 10 October 2015 to 4 March 2016. (b) The InSAR-derived deformation for the same period.

Figure 10. (a–e) Five time-series results filtered using a binary filter with a threshold of 2.6 mm, encompassing the coal mining period.

Figure 11. The RM-DBC results for the study area. (a) The subsidence profile along the 234° direction. (b) The performance for all directions, where the length of the polar axis represents the sum of the pixel values in each direction and the angle represents the angle between the ray and the north direction.

Figure 12. The goaf-locating results using the S-ODM. The upper and bottom layers are the LOS deformation based on the S-ODM-inversed and InSAR-derived information, respectively. The pink and blue cuboids represent the in situ measured goaf and the S-ODM-estimated goaf, respectively.

Figure 13. Comparison of the deformation results generated by the S-ODM, PIM, ODM, and D-InSAR for the strike and dip directions, respectively.

Figure 14. Influence of monitoring errors on the estimated geometric parameters of underground goaves.

Table 1

Parameters of interferogram pairs.

Interferogram Pair	Spatial Baseline	Time Baseline
20150919–20151013	101.8260 m	24 d
20151013–20151106	47.2675 m	24 d
20151106–20151130	3.1222 m	24 d
20151130–20151224	−30.4078 m	24 d
20151224–20160117	−43.6049 m	24 d
20160117–20160305	−37.7825 m	48 d

Table 2

Comparison between the simulated results generated using FLAC3D and the estimated parameters using GA-PSO.

Parameter	Simulated Value	Estimated Value	Absolute Error	Relative Error
X(m)	1000	1007.26	7.26	-
Y(m)	1000	995.38	4.62	-
L(m)	600	627.72	27.72	4.62%
W(m)	100	90.36	9.64	9.64%
H(m)	500	519.51	19.51	3.9%
$\phi {(}^{\circ })$	30	26.87	3.13	-
$\alpha {(}^{\circ })$	45	45.77	0.77	-
m(m)	5	4.65	0.35	7%

Table 3

Estimation results and errors of the underground goaf parameters for the study area.

Parameter	Real Value	Estimated Value	Absolute Error	Relative Error
X(m)	1057	1078.66	21.66	-
Y(m)	565	550.38	14.62	-
L(m)	493	522.21	19.21	5.92%
W(m)	142	161.86	19.86	13.98%
H(m)	740	701.6	38.4	5.19%
$\phi {(}^{\circ })$	13	10.98	2.02	-
$\alpha {(}^{\circ })$	236	234	2	-
m(m)	4.5	4.70	0.20	4.4%

Table 4

Comparison between the true values and the S-ODM, PIM, and ODM results.

Parameter	True Value	S-ODM		PIM		ODM
		Estimated Value	Absolute Error	Estimated Value	Absolute Error	Estimated Value	Absolute Error
X(m)	1057	1078.66	21.66	1089.31	32.31	1134.14	74.14
Y(m)	565	550.38	14.62	545.15	19.85	512.64	52.36
L(m)	493	522.21	19.21	510.62	17.62	465.14	27.86
W(m)	142	161.86	19.86	167.14	25.14	183.87	41.87
H(m)	740	701.6	38.4	687.83	52.17	680.15	59.85
$\phi {(}^{\circ })$	13	10.98	2.02	10.08	2.92	10.12	2.88
$\alpha {(}^{\circ })$	236	234	2.0	241	5.0	221	15
m(m)	4.5	4.7	0.2	4.64	0.14	4.86	0.36

Table 5

Comparison between the simulated results generated using the S-ODM and the estimated parameters using GA-PSO.

Parameter	Simulated Value	Estimated Value	Absolute Error	Relative Error
X(m)	2000	2007.26	7.26	-
Y(m)	2000	1995.38	4.62	-
L(m)	500	517.72	17.72	3.54%
W(m)	100	91.36	8.64	8.64%
H(m)	500	519.51	19.51	3.9%
$\phi {(}^{\circ })$	15	15.87	3.87	-
$\alpha {(}^{\circ })$	45	45.77	0.77	-
m(m)	5	4.65	0.35	7%

References

1. Mao, J.; Tong, H. Coal resources, production, and use in China. Coal Handb.; 2023; 2, pp. 431-454. [DOI: https://dx.doi.org/10.1016/B978-0-12-824327-5.00004-1]

2. Fan, H.D.; Cheng, D.; Deng, K.Z.; Chen, B.Q.; Zhu, C.G. Subsidence monitoring using D-InSAR and probability integral prediction modelling in deep mining areas. Emp. Surv. Rev.; 2015; 47, pp. 438-445. [DOI: https://dx.doi.org/10.1179/1752270614Y.0000000153]

3. Yang, Z.; Li, Z.; Zhu, J.; Wang, Y.; Wu, L. Use of SAR/InSAR in mining deformation monitoring, parameter inversion, and forward predictions: A review. IEEE Geosci. Remote Sens. Mag.; 2020; 8, pp. 71-90. [DOI: https://dx.doi.org/10.1109/MGRS.2019.2954824]

4. Li, T.; Zhang, H.; Fan, H.; Zheng, C.; Liu, J. Position inversion of goafs in deep coal seams based on DS-InSAR data and the probability integral methods. Remote Sens.; 2021; 13, 2898. [DOI: https://dx.doi.org/10.3390/rs13152898]

5. Maus, V.; Werner, T.T. Impacts for half of the world’s mining areas are undocumented. Nature; 2024; 625, pp. 26-29. [DOI: https://dx.doi.org/10.1038/d41586-023-04090-3]

6. Wang, X.; Zhen, F.; Huang, X.; Zhang, M.; Liu, Z. Factors influencing the development potential of urban underground space: Structural equation model approach. Tunn. Undergr. Space Technol.; 2013; 38, pp. 235-243. [DOI: https://dx.doi.org/10.1016/j.tust.2013.06.005]

7. Wang, J.; Wang, C.; Huang, X.; Zou, J. Research on multi-frequency ultrasonic scanning detecting technology of cavity in the test borehole. Bull. Eng. Geol. Environ.; 2021; 80, pp. 1249-1264. [DOI: https://dx.doi.org/10.1007/s10064-020-01979-5]

8. Zahorec, P.; Pašteka, R.; Papčo, J.; Putiška, R.; Mojzeš, A.; Kušnirák, D.; Plakinger, M. Mapping hazardous cavities over collapsed coal mines: Case study experiences using the microgravity method. Near Surf. Geophys.; 2021; 19, pp. 353-364. [DOI: https://dx.doi.org/10.1002/nsg.12139]

9. Chang, J.; Su, B.; Malekian, R.; Xing, X. Detection of water-filled mining goaf using mining transient electromagnetic method. IEEE Trans. Ind. Inform.; 2019; 16, pp. 2977-2984. [DOI: https://dx.doi.org/10.1109/TII.2019.2901856]

10. Yuan, G.; Wang, Y.; Zhao, F.; Wang, T.; Zhang, L.; Hao, M.; Yan, S.; Dang, L.; Peng, B. Accuracy assessment and scale effect investigation of UAV thermography for underground coal fire surface temperature monitoring. Int. J. Appl. Earth Obs. Geoinf.; 2021; 102, 102426. [DOI: https://dx.doi.org/10.1016/j.jag.2021.102426]

11. Yuan, G.; Wang, Y.; Zhao, F.; Yan, S.; Zhang, H.; Lang, F.; Hao, M.; Cao, F.; Peng, B.; Dang, L. et al. Spatiotemporal Correlation Characteristics between Thermal Infrared Remote Sensing Obtained Surface Thermal Anomalies and Reconstructed 4D Temperature Fields of Underground Coal Fires. IEEE Trans. Geosci. Remote Sens.; 2023; 61, pp. 1-18. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3313196]

12. Wang, T.; Wang, Y.; Zhao, F.; Feng, H.; Liu, J.; Zhang, L.; Zhang, N.; Yuan, G.; Wang, D. A spatio-temporal temperature-based thresholding algorithm for underground coal fire detection with satellite thermal infrared and radar remote sensing. Int. J. Appl. Earth Obs. Geoinf.; 2022; 110, 102805. [DOI: https://dx.doi.org/10.1016/j.jag.2022.102805]

13. Zhu, J.; Li, Z.; Hu, J. Research Progress and Methods of InSAR for Deformation Monitoring. Acta Geod. Et Cartogr. Sin.; 2017; 46, pp. 1717-1733. [DOI: https://dx.doi.org/10.11947/j.AGCS.2017.20170350]

14. Gao, H.; Liao, M.; Liu, X.; Xu, W.; Fang, N. Source geometry and causes of the 2019 Ms6. 0 Changning earthquake in Sichuan, China based on InSAR. Remote Sens.; 2022; 14, 2082. [DOI: https://dx.doi.org/10.3390/rs14092082]

15. Anantrasirichai, N.; Biggs, J.; Albino, F.; Bull, D. A deep learning approach to detecting volcano deformation from satellite imagery using synthetic datasets. Remote Sens. Environ.; 2019; 230, 111179. [DOI: https://dx.doi.org/10.1016/j.rse.2019.04.032]

16. Cai, J.; Liu, G.; Jia, H.; Zhang, B.; Wu, R.; Fu, Y.; Xiang, W.; Mao, W.; Wang, X.; Zhang, R. A new algorithm for landslide dynamic monitoring with high temporal resolution by Kalman filter integration of multiplatform time-series InSAR processing. Int. J. Appl. Earth Obs. Geoinf.; 2022; 110, 102812. [DOI: https://dx.doi.org/10.1016/j.jag.2022.102812]

17. Zhang, N.; Wang, Y.; Zhao, F.; Wang, T.; Zhang, K.; Fan, H.; Zhou, D.; Zhang, L.; Yan, S.; Diao, X. et al. Monitoring and Analysis of the Collapse at Xinjing Open-Pit Mine, Inner Mongolia, China, Using Multi-Source Remote Sensing. Remote Sens.; 2024; 16, 993. [DOI: https://dx.doi.org/10.3390/rs16060993]

18. Liu, X.; Wang, Y.; Yan, S.; Shao, Y.; Zhou, H.; Li, Y. Ground subsidence characteristics associated with urbanization in East China analyzed with a Sentinel-1A-based InSAR time series approach. Bull. Eng. Geol. Environ.; 2019; 78, pp. 4003-4015. [DOI: https://dx.doi.org/10.1007/s10064-018-1383-6]

19. Zhang, Y.; Wang, Y.; Huo, W.; Zhao, F.; Hu, Z.; Wang, T.; Song, R.; Liu, J.; Zhang, L.; Fernández, J. et al. Ground Deformation Monitoring over Xinjiang Coal Fire Area by an Adaptive ERA5-Corrected Stacking-InSAR Method. Remote Sens.; 2023; 15, 1444. [DOI: https://dx.doi.org/10.3390/rs15051444]

20. Wang, Y.; Yuan, G.; Wang, T.; Liu, J.; Zhao, F.; Feng, H.; Dang, L.; Peng, K.; Zhang, L. Research on Multi-source Remote Sensing Detection of Concealed Fire Sources in Coalfields. Geomat. Inf. Sci. Wuhan Univ.; 2022; 47, pp. 1651-1661. [DOI: https://dx.doi.org/10.13203/j.whugis20220184]

21. Hu, Z.; Ge, L.; Li, X.; Zhang, K.; Zhang, L. An underground-mining detection system based on DInSAR. IEEE Trans. Geosci. Remote Sens.; 2012; 51, pp. 615-625. [DOI: https://dx.doi.org/10.1109/TGRS.2012.2202243]

22. Du, S.; Wang, Y.; Zheng, M.; Zhou, D.; Xia, Y. Goaf locating based on InSAR and probability integration method. Remote Sens.; 2019; 11, 812. [DOI: https://dx.doi.org/10.3390/rs11070812]

23. Yang, Z.; Li, Z.; Zhu, J.; Yi, H.; Feng, G.; Hu, J.; Wu, L.; Preusse, A.; Wang, Y.; Papst, M. Locating and defining underground goaf caused by coal mining from space-borne SAR interferometry. ISPRS J. Photogramm. Remote Sens.; 2018; 135, pp. 112-126. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2017.11.020]

24. Wang, T.; Zhao, F.; Wang, Y.; Zhang, N.; Zhou, D.; Diao, X.; Zhao, X. An Algorithm for Locating Subcritical Underground Goaf Based on InSAR Technique and Improved Probability Integral Model. IEEE Trans. Geosci. Remote Sens.; 2023; 61, pp. 1-14. [DOI: https://dx.doi.org/10.1109/TGRS.2023.3308565]

25. He, L.; Wang, X.; Liu, C.; Tang, Y.; Zhang, X.; Kang, J.; Guo, C.; He, R. Time series interferometric synthetic aperture radar-based modeling and analysis of complex land subsidence caused by multi-seam coal mining on the Liaohe Plain, China. J. Appl. Remote Sens.; 2022; 16, 024512. [DOI: https://dx.doi.org/10.1117/1.JRS.16.024512]

26. Wang, R.; Wu, K.; He, Q.; He, Y.; Gu, Y.; Wu, S. A novel method of monitoring surface subsidence law based on probability integral model combined with active and passive remote sensing data. Remote Sens.; 2022; 14, 299. [DOI: https://dx.doi.org/10.3390/rs14020299]

27. Kinoshita, Y.; Furuta, R. Slow slip event displacement on 2018 offshore Boso Peninsula detected by Sentinel-1 InSAR time-series analysis with numerical weather model assistance. Geophys. J. Int.; 2024; 237, pp. 75-89. [DOI: https://dx.doi.org/10.1093/gji/ggae028]

28. Pepe, A.; Lanari, R. On the extension of the minimum cost flow algorithm for phase unwrapping of multitemporal differential SAR interferograms. IEEE Trans. Geosci. Remote Sens.; 2006; 44, pp. 2374-2383. [DOI: https://dx.doi.org/10.1109/TGRS.2006.873207]

29. Crosetto, M.; Monserrat, O.; Cuevas-González, M.; Devanthéry, N.; Crippa, B. Persistent scatterer interferometry: A review. ISPRS J. Photogramm. Remote Sens.; 2016; 115, pp. 78-89. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2015.10.011]

30. Zheng, M.; Deng, K.; Fan, H.; Du, S. Monitoring and analysis of surface deformation in mining area based on InSAR and GRACE. Remote Sens.; 2018; 10, 1392. [DOI: https://dx.doi.org/10.3390/rs10091392]

31. Aspri, A.; Beretta, E.; Mazzucato, A.L.; De Hoop, M.V. Analysis of a model of elastic dislocations in geophysics. Arch. Ration. Mech. Anal.; 2020; 236, pp. 71-111. [DOI: https://dx.doi.org/10.1007/s00205-019-01462-w]

32. Okada, Y. Surface deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am.; 1985; 75, pp. 1135-1154. [DOI: https://dx.doi.org/10.1785/BSSA0750041135]

33. Sheikhalishahi, M.; Ebrahimipour, V.; Shiri, H.; Zaman, H.; Jeihoonian, M. A hybrid GA–PSO approach for reliability optimization in redundancy allocation problem. Int. J. Adv. Manuf. Technol.; 2013; 68, pp. 317-338. [DOI: https://dx.doi.org/10.1007/s00170-013-4730-6]

34. Gercek, H. Poisson’s ratio values for rocks. Int. J. Rock Mech. Min. Sci.; 2007; 44, pp. 1-13. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2006.04.011]