Abstract

High-intensity underground mining generates considerable surface subsidence in mining areas, including ground cracks and collapse pits on roads and farmland, threatening the safety of buildings. Large-amplitude subsidence (e.g., >2 m) is usually characterized by a large phase gradient in interferograms, leading to severe phase decorrelation and unwrapping errors. Therefore, the subsidence on the surface cannot be well derived simply using conventional differential interferometric synthetic aperture radar (DInSAR) or other geodetic measurements. We propose a new method that combines both DInSAR and subpixel offset-tracking technology to improve mine subsidence monitoring over large areas. We utilize their respective advantages to extract both the spatial boundaries and the amplitude of displacements. Using high-resolution RADARSAT-2 SAR images (5 m) acquired on February 13, 2012, and November 27, 2012, in the Shendong Coalfield located at the border between Shaanxi Province and Inner Mongolia Province, China, we obtain the subcentimetre-level subsidence of the mine boundary by DInSAR and resolve the metre-level mine subsidence centre based on subpixel offset tracking. The whole subsidence field is obtained by combining and analyzing the subcentimetre-level and the metre-level subsidence. We use the probability integral method (PIM) function model to fit the boundary and central mine subsidence to reconstruct the spatial distribution of the mine subsidence. Our results show that the maximum central subsidence reaches ~4.0 m (beyond the monitoring capabilities of DInSAR), which is generally in agreement with the maximum subsidence of ~4.0-5.0 m from field investigation. We also model the boundary and the central subsidence (the final fitting coefficient is 0.978). Our findings indicate that the offset-tracking method can compensate for the deficiency of DInSAR in large-amplitude subsidence extraction, and the inclusion of the PIM technique helps reconstruct the whole subsidence field in mining areas.

1. Introduction

Large-scale exploitation of coal resources meets the increasing demands of regional economic development. However, it also results in serious ecological, environmental, and geological problems. Ground subsidence in mining areas has caused a series of geological disasters, which brings potential dangers related to the construction of mining areas, including threats to water resources, transportation infrastructure, and farmland. The major technical characteristics of high-intensity mining exploitation are a large working face, the rapidly advancing speed of the working face, a small ratio of mining depth to thickness, serious damage to the overlying strata, and substantial surface movement and deformation [1, 2]. The Shendong mining area, located at the border between Shaanxi Province and Inner Mongolia Province, has sufficient coal resources with shallow buried depth (average mining depth of ~260 m) and small coal seam thickness (2.2–2.9 m, 2.5 m on average). Ground subsidence caused by coal mining in this area has the same characteristics as that caused by high-intensity mining, such as high subsidence rates (maximum subsidence reaching ~4 m/yr). Therefore, obtaining the spatial distribution and amplitude of surface subsidence in these areas is helpful for understanding the potential geological hazards and providing data for comprehensive management support. This information is also of great importance for ecological protection and sustainable development in mining regions.

Interferometric synthetic aperture radar (InSAR) is a useful geodetic tool for deformation monitoring that has been rapidly developed in recent years. Differential interferometric synthetic aperture radar (DInSAR) can overcome the shortcomings of traditional geodetic measurements for surface deformation and can thus be used to monitor ground subsidence in large areas [3, 4]. A great number of applications related to DInSAR deformation monitoring of earthquake cycles [5, 6], volcanic activities [7, 8], glacial shift [9, 10], urban water-loss settlement [11, 12], and landslides [13, 14] have been documented. The application of DInSAR to monitoring mining subsidence has also been studied by many authors [1518]. However, DInSAR is susceptible to space-time phase decorrelation in areas with large displacement gradients, which may induce unwrapping errors. There are still technical bottlenecks in monitoring rapid and dynamic deformation processes caused by mining subsidence [19, 20]. Small baseline subset InSAR (SBAS InSAR) [2125] and permanent scatter InSAR (PS InSAR) [2630] have recently been used to study subsidence in mining areas. To a certain extent, the accuracy has been greatly improved and the observation period of mining subsidence can reach months to years; however, there are still some deficiencies in measuring the substantial deformation due to the rapid pace of subsidence (i.e., maximum subsidence reaching metre-level deformation) [31, 32]. Offset-tracking technology based on SAR image intensity information is a powerful complement to DInSAR technology and can be used to resolve the substantial displacement caused by rapid subsidence in mining areas [33].

In this paper, we use DInSAR and offset-tracking technology to generate ground subsidence fields in the Shendong Coalfield. To this end, DInSAR technology is initially used to extract the subsidence boundary, and offset-tracking technology is subsequently used to resolve a substantial displacement near the central subsidence. While DInSAR and offset tracking are able to obtain the discrete displacement measurements of ground subsidence in boundary and central mining area, they cannot derive the whole shape of the subsidence field in a straightforward way. In order to further analyse the overall pattern, spatial distributions, and deformation characteristics of the subsidence field, we reconstruct the whole subsidence field in this area by the probability integral method (PIM) using InSAR and offset-tracking results. Finally, the results from these two methods are then combined to analyse and reconstruct the ground subsidence form of the goaf according to the probability integral method.

2. Study Area and SAR Data

2.1. Study Area

Shendong Coalfield is the general name for the Shenfu Coalfield in Shaanxi Province and the Dongsheng Coalfield in Inner Mongolia Province. The Shendong Coalfield is located in the northern region of Shenmu County, the western region of Fugu County, the southern region of Yijin Holuo Banner in Inner Mongolia, and the southwest region of Zhungeer Banner [34]. The coalfield includes Daliuta, Bulianta, and other mines with 10 million tons of coal. The Shendong Coalfield is one of the largest coal production mines in China and one of the eight largest coal mines in the world [35, 36]. The mining area is mainly located on both banks of the Ulanmulun River. The north–south length of the mining area is approximately 38–90 km, and its east–west width is approximately 35–55 km, covering an area of approximately 3481 km2. The proven reserves are substantial. The study area (shown as the yellow triangle position in Figure 1(c)) is located in Buertai mine and Cuncaota no. 1 and no. 2 mines in the Shendong Coalfield and is a typical near-horizontal shallowly buried coal seam. The average thickness of the coal seam is 2.5 m, and the mining depth is 265 m in the 22201 working face of Buertai mine, which is currently the largest designed output of a single shaft in the world [37]. The mining subsidence of this mine has a large influence range and causes serious ground subsidence, including ground cracking.

2.2. Data

The SAR images were collected from RADARSAT-2 on February 13, 2012, and November 27, 2012. The time interval of the SAR data was 289 days. The spatial resolution of the image is 5 m on the ground. The satellite adopted the C wave band as the working wave band, HH polarization mode, multilook fine beam, and a  km2 coverage area. The digital elevation model (DEM) data required for processing were 90 m elevation data of the Shuttle Radar Topography Mission (SRTM DEM). The SAR image was cut via a region of interest to obtain an image covering the Buertai mining, Cuncaota no. 1 and no. 2 mine areas (109°55-110°06E, 39°23-39°31N) (Figure 1). The basic information of the SAR images is shown in Table 1.

3. Methods

3.1. DInSAR

DInSAR technology mainly uses phase information of SAR images to extract ground deformation [5, 7]. The basic principle of the algorithm is to multiply two SAR images covering the same area at different times. The differential interferometry phase is processed, and the complex interferogram is converted to the ground deformation (containing subsidence, etc.) information .

After registration of the master and slave images, the interferometric phase of the same pixels in the image pair can be calculated. The interferometric phase usually consists of five parts, as shown in where represents the topographic phase obtained from the external DEM; represents the slant deformation phase between different transit times of the radar satellite; represents the flattened phase obtained from the external DEM; represents atmospheric phase variation; and represents the noise phase, which contains InSAR processing errors, thermal noise effects, etc. The differential interferogram phase is obtained by subtracting the topographic phase and the flattened phase from the interferogram phase .

The deformation phase can be obtained by subtracting the atmospheric and noise phases from the differential interferometric phase. Finally, the unwrapped deformation phase can be obtained by unwrapping the deformation phase. The unwrapped deformation phase can be converted to the ground deformation value along the line of sight of radar satellite [19].

3.2. Offset Tracking

Subpixel offset-tracking technology mainly uses the intensity information of two SAR images covering the same area for subpixel registration to obtain a large number of coordinate offsets for the same pixels [38, 39]; then, the subpixel offsets of the same pixels are decomposed into the slant range (along the line of sight of the satellite) and the azimuth (along the satellite orbit direction) offsets, as shown in Figure 2. There are two ways to implement the offset-tracking algorithm, i.e., the intensity tracking algorithm and coherence tracking algorithm, and these two methods have different application conditions. Because of the multiple farms and the vegetation coverage on the surface of the mining area, the coherence level of SAR images in the mining area is often low. Therefore, offset-tracking technology based on the intensity tracking algorithm is more suitable for monitoring mining subsidence. The intensity tracking algorithm calculates the intensity information of two SAR images by normalized cross-correlation (NCC) [40]. This algorithm finds the peak value of the intensity correlation coefficient of feature pixels in two SAR images. When the feature elements of the feature pixels in the image are highly correlated, a high monitoring accuracy can be obtained using a small search window. If the pixel feature elements are not coherent, it is necessary to resample the multilook images to reduce the image resolution at the cost of increasing the coherence of the surface features for deformation monitoring [41, 42].

3.3. PIM

Mining subsidence prediction methods mainly include empirical methods, theoretical simulation methods, and influence function methods [43]. After decades of research on mining subsidence prediction methods, the PIM is widely used and relatively mature in China. The PIM considers an underground rock mass as a random medium, and the motion process of the rock mass is considered a random motion process obeying statistical law [44]. The theoretical basis of this method is random medium theory. The prediction of surface subsidence in the mining area is determined by the PIM. It is considered that the movement law of underground rock mass is similar to that of particles in a random medium, and a discontinuous medium movement model is used to study the movement law of rock mass [45, 46]. The movement model of random medium particles is simplified to a simple sphere model in Figure 3(a). It is assumed that the rock media are composed of spheres with similar shape and equal size and weight, which are evenly arranged between the surface and the mining coal seams. When the sphere in the first floor moves, the probability that the two spheres in the second layer will fall into the first layer due to gravity is 1/2. Based on the probability statistics method, the probability that the third layer balls will fall into the first layer is 1/4, 2/4, and 1/4, respectively, and the probability that the fourth layer balls will fall into the first layer is 1/8, 3/8, 3/8, and 1/8, respectively. Accordingly, the probability distribution function of the balls falling into the first layer is deduced as shown by the dotted line in Figure 3(a). If the particle size is equivalent to an infinitesimal size, the probability distribution histogram will evolve into a smooth approximate normal distribution curve, as shown by the blue solid line in Figure 3(a). The PIM is used to predict ground deformation by integrating the probability density function to obtain subsidence basin curves [4749]. The concrete theoretical formulas are shown as follows:

The subsidence of any point on the surface caused by mining of underground random discrete mining units is set as follows:

The PIM for calculating the ground subsidence at any point caused by mining of all underground discrete elements on the whole mining face (strike length , dip width ) is as follows: where is the maximum mining subsidence value; is the main influencing radius of subsidence, and , where is the average mining depth and is the main influencing angle tangent;, is the calculation unit depth, and is the maximum subsidence angle; is the plane coordinate of the working face at the centre of the mining unit; and is the terrain coordinate at any point. Based on the above probability integral formula, the relationship between underground movement and ground subsidence is established.

4. Results and Analysis

4.1. DInSAR Results and Analysis

Using RADARSAT-2 SAR image pairs acquired on February 13, 2012, and November 27, 2012 (Table 1), and SRTM DEM v3.0 (90 m) data, the deformation field due to surface subsidence in the mining area that occurred during the observation period was derived using the radar interferometry software GAMMA (GAMMA Remote Sensing Research and Consulting AG, Gümligen, Switzerland; http://www.gamma-rs.ch/). The ground subsidence of the goaf can be resolved by transforming the LOS (line of sight) deformation phase into vertical deformation assuming that the observed displacement field is generated by pure subsidence movement in the vertical direction.

As shown in Figure 4, there are seven obvious subsidence fields (marked by A, B, C, D, E, F, and G) in the study area, each of which is approximately elliptical. Due to the characteristics of large-scale land subsidence, the high gradient of subsidence, and the sudden onset of high-intensity mining subsidence, serious spatial and temporal decoherence in the interferogram occurs in the central area of the subsidence area (the coloured clutter area enclosed by red ellipses in Figure 4(a)). The long-time interval of the SAR images further aggravates the decoherence phenomenon. DInSAR technology requires high coherence to ensure the accuracy of deformation monitoring. Therefore, the coherence threshold is set to 0.4, and the part with a coherence of less than 0.4 is masked in phase unwrapping (the black area enclosed by red ellipses in Figure 4(b)), which is not involved in the final deformation calculation. Therefore, DInSAR technology cannot obtain subsidence information in the central area of the subsidence field; only the deformation information at the edge of the subsidence field is monitored and is between -0.01 and -0.09 m (Figure 5). The maximum subsidence of the subsidence edge is -0.091 m at field .

Although the subsidence from -0.01 to -0.09 m at the edge of the subsidence field was monitored by DInSAR, because of the high surface subsidence gradient and sudden onset of high subsidence rates, unwrapping failure or even errors occurred in some subsidence values in the subsidence range from -0.01 to -0.09 m at the edge of the subsidence field. Therefore, we consider that the subsidence between -0.01 and -0.02 m at the edge of the mining area extracted by DInSAR is more accurate in this study. According to mining subsidence [43], the subsidence field boundary is delineated by taking the -0.01 m subsidence position as the boundary of the subsidence area, such as the light blue boundary in the red ellipses in Figure 6. Using DInSAR technology, the boundary of the mine goaf can be calculated, and the subsidence boundary of each goaf can be delineated by the subsidence position of -0.01 m.

4.2. Offset-Tracking Results and Analysis

The DInSAR displacement field (Figures 5 and 6) shows that there are still monitoring blind areas in the central area of the mine goaf due to phase decorrelation, and large-magnitude subsidence in the central mining area cannot be well obtained. Therefore, we use the offset-tracking method to better estimate subsidence with large displacements to address the shortcoming of the DInSAR method. Additionally, offset tracking resolves displacements along the slant range and azimuth range. Considering that the deformation produced by mining activities is mainly characterized by vertical subsidence, we transformed slant-range displacements into vertical deformation by ( represents vertical deformation, represents slant range deformation, and is the beam incidence angle of the satellite). During the offset-tracking processing, a 1 : 2 multilook factor was used to reduce the image resolution as well as to increase the correlation of characteristic pixels. The results show long strips of deformation in the subsidence field (blue strips of deformation are indicated by the red ellipses shown in Figure 7). The surface subsidence in Figure 8 was almost between -1.0 m and -4.0 m, which was generally in agreement with the field investigations of mining subsidence mostly from -1.0 to -3.0 m; the peak subsidence reaches -4.0 ~ -5.0 m in this mining area [32]. The results are generally consistent with the maximum mining subsidence of ~4.478 m in the same mining area presented by Fan et al. [50] using 11 TerraSAR-X SAR images acquired from 13 December 2012 to 2 April 2013 based on a subpixel offset-tracking method. It is also in agreement with the result by Wang et al. [51] which obtained the maximum subsidence value of ~4.5 m based on the offset-tracking technology using TerraSAR-X images acquired in the same time period. From the offset-tracking result, most of the surface subsidence is approximately -2.0 m.

4.3. Joint Analysis of DInSAR and Offset-Tracking Results

Comparison and analysis of the results of DInSAR and offset tracking show that DInSAR can obtain the boundary of subsided mining areas more effectively, and offset tracking can measure the large-amplitude subsidence (metre-level) in the central area. To better utilize the advantages of DInSAR and offset tracking, the results of DInSAR and offset-tracking are merged into a subsidence map in the following analysis based on the band math method. The inferred spatial boundary in the final subsidence map from InSAR is used as the first band image (Figure 6), and the subsidence centre results obtained from offset tracking are regarded as the second band image (Figure 7). Actually, these two results from InSAR and offset tracking share the same geographical reference frame; thus, the geographical locations in the two measurements are also exactly the same. Using the band math method, the above two-band images are directly superposed to generate our final subsidence map shown in Figure 9. The boundary deformation of the subsidence field obtained by DInSAR within -0.01 to -0.02 m (Figure 6) is combined with the metre-level central subsidence obtained by offset tracking (Figure 7). The fusion results shown in Figure 9 indicate that the geographic location and spatial distribution of mining subsidence monitored by the DInSAR and offset-tracking methods are in relatively good agreement. The red parts in the black ellipses in Figure 9(a) represent the DInSAR results, and the blue parts represent the offset-tracking results. As shown, the red subsidence boundary acquired by DInSAR surrounds the blue central subsidence acquired by offset tracking. The characteristics of the subsidence centre surrounded by the subsidence boundary of the goaf are obvious and are shown in almost every subsidence field. To study the spatial relationship of the subsidence field acquired by the two monitoring methods in detail, local enlargements of subsidence fields B, D, and G in Figure 9(a) are used to display the layout details of the subsidence field. The light green-dotted lines in Figure 9(b), 9(d) and 9(g) delineate the boundary of the subsidence field monitored by DInSAR, while the black-dotted box delineates the large-scale deformation of the centre of the subsidence field obtained by offset tracking. The black-dotted line frames in Figure 9(b), 9(d) and 9(g) are located in the centre of the subsidence field. The location and scope of the two technology results are also generally consistent.

4.4. PIM Prediction and Analysis

To reconstruct the whole spatial pattern of the mining subsidence based on the results from DInSAR and offset tracking, a 22201-1/2 working face (black-dashed box in Figure 9(b), 9(d)) was selected to conduct a preliminary study of mining subsidence. The abovementioned PIM was adopted for the mining subsidence analysis. According to the 22201-1/2 mining plan and SAR data acquisition time, the position and progress of the working face are obtained. The working face has a strike length of ~2103.6 m, a dip length of ~311.4 m, a coal seam dip of ~1°, an average mining depth of ~260 m, and a coal seam mining thickness of ~2.5 m. The predicted parameters refer to the parameters obtained by the adjacent surface movement observations [35]. The subsidence coefficient , tangent value of the main influence angle tanβ, and propagation angle of the mining influence are set to 0.62, 1.63, and 89.0°, respectively. The surface observations in this model are from 2556 surface observation points acquired by the DInSAR and offset-tracking methods. These points are distributed at the boundary and centre of the subsidence field (as shown by the blue points indicated by arrows in Figure 10(a)). We transformed the coordinate system of the subsidence measurement into the coordinate system of the mining face through coordinate translation and rotation. Then, the result of the mining subsidence analysis was obtained through multiple least square regression operations [52, 53]. Figure 10(a) shows that the centimetre-level subsidence boundary obtained by InSAR well constrains the inversion of the subsidence field boundary, while the subsidence centre value obtained by the offset-tracking method is used to fit the maximum subsidence value of the subsidence field centre. The inversion of the subsidence field (Figure 10(b)) shows an elongated funnel-shaped distribution along the mining face, and the inversion of the subsidence location coincides well with the observed subsidence location. The inversion maximum value is approximately ~1.75 m, which is slightly different from the maximum subsidence obtained by the offset tracking. Because we use the probabilistic integral method to invert the average value of whole subsidence field, while the offset-tracking method is used to obtain the local subsidence value of a single discrete point, there will be differences in the maximum subsidence value. The inversion value of the subsidence field can be seen in the contour map of the subsidence in Figure 10(c).

Two sections (labelled longitudinal profile and transverse profile in Figure 11(b)) along the vertical and horizontal directions are made in the subsidence field. The sections include the centimetre-level deformation of the goaf boundary monitored by DInSAR, the central metre-level deformation monitored by the offset tracking, and the inversion profile of the subsidence field. Additionally, the surface subsidence profile of the goaf is obtained, as shown in Figures 11(a) and 11(c). Figures 11(a) and 11(c) show that the observed values (red dots) have a high fitting degree with the inversion subsidence curves. To further analyse the fitting accuracy of the subsidence values and inversion values of the observation subsidence points, the monitoring values of all observation points are substituted into a probability integral model to calculate the corresponding inversion value, and the residual of the inversion value and observation value is calculated. The fitting determination coefficient is also calculated, and the high coefficient indicates that the observation values and inversion values have a high fitting degree. It is concluded that inversion of the goaf subsidence based on a probabilistic integral model has high reliability. Thus, a probabilistic integral function can be used to evaluate goaf subsidence and recover the subsidence value of any point in the subsidence field.

5. Conclusions

To better monitor the large-scale subsidence in mining areas, which is featured by high displacement gradients in interferograms, we take technical advantages of the DInSAR and offset-tracking methods in this study. This method is able to extract microscale subsidence boundary information by employing DInSAR, and it can map large-amplitude deformation in the central mine by employing offset-tracking technology. We combine both phase and intensity information of SAR images to observe patterns and features of the subsidence. We use the PIM model to restore the complete spatial subsidence field in the mining area. The main conclusions are derived as follows: (1)We use DInSAR to extract the subsidence boundary of the mining basin in the goaf, and -0.01 m subsidence is used in the calculation. The long-term (289 days) surface subsidence pattern and amplitude of the goaf are obtained(2)To compensate for the disadvantage of DInSAR in measuring high phase gradients and large-amplitude deformation in the centre of the goaf, we use the offset-tracking method to generate the central subsidence. The surface subsidence was almost between -1.0 m and -4.0 m. The maximum subsidence reaches ~4.0 m, which is generally in agreement with the maximum subsidence of -4.0 to -5.0 m from the field investigations. Our subsidence observations from offset tracking are also consistent with the peak mining subsidence value of ~4.478 m in the almost same mining area presented by Fan et al. [50] using 11 TerraSAR-X SAR images from 13 December 2012 to 2 April 2013 based on a subpixel offset-tracking method. It is also in agreement with the offset-tracking measurements by Wang et al. [51] which obtained the maximum subsidence value of ~4.5 m using TerraSAR-X images acquired in the same time period. The boundary subsidence from DInSAR and the central subsidence from offset tracking are fitted by the probability integral function model (the fitting coefficient is 0.978)(3)Our results indicate that the combined DInSAR and offset-tracking technology can be used to map the surface subsidence subjected to severe decorrelation due to high subsidence magnitudes and large deformation gradients in these areas. The PIM model can restore the whole subsidence field. Combining DInSAR and offset-tracking technology is a practical tool for mining subsidence measurements. Effective monitoring of mining subsidence has potential to serve as an important reference for geological hazard assessments in mining areas, prevention of geological hazards in mines, and ecological restoration of mining areas.

Data Availability

The data used to support the study is available upon request to the corresponding author.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors would like to thank AJE® (American Journal Experts, Co.) for her English editing in enhancing the quality of the manuscript. This research was jointly funded by the National Natural Science Foundation of China (Grant No. 41701515), the National Natural Science Foundation of China and the Shenhua Coal Industry Group Co., Ltd. (Grant No. U1261106 and Grant No. U1261206), and the State Key Laboratory of Earthquake Dynamics (Grant No. LED2018b04).