Inversion and Analysis of Mining Subsidence by Integrating DInSAR, Offset Tracking, and PIM Technology

School of Environmental Science and Engineering, Suzhou University of Science and Technology, Suzhou 215009, China Department of Remote Sensing Science and Technology, Henan Polytechnic University, Jiaozuo 454000, China Key Laboratory of Spatio-temporal Information and Ecological Restoration of Mines (MNR), Henan Polytechnic University, Jiaozuo 454000, China State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China


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 highintensity 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 [15][16][17][18]. However, DInSAR is susceptible to spacetime 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) [21][22][23][24][25] and permanent scatter InSAR (PS InSAR) [26][27][28][29][30] 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.

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 Coal-field 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 northsouth 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 km 2 . 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.

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 50 * 50 km 2 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°06′E, 39°23′-39°31′N) ( Figure 1). The basic information of the SAR images is shown in Table 1.

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 Δr.
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 ϕ top represents the topographic phase obtained from the external DEM; ϕ def represents the slant deformation phase between different transit times of the radar satellite; ϕ flat represents the flattened phase obtained from the external DEM; ϕ atm represents atmospheric phase variation; and ϕ noi 2 Journal of Sensors represents the noise phase, which contains InSAR processing errors, thermal noise effects, etc. The differential interferogram phase is obtained by subtracting the topographic phase ϕ top and the flattened phase ϕ flat from the interferogram phase ϕ int .
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 Δr along the line of sight of radar satellite [19].

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, offsettracking 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 crosscorrelation (NCC) [40]. This algorithm finds the peak value of the intensity correlation coefficient of feature pixels in  3 Journal of Sensors 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].

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 [47][48][49]. The concrete theoretical formulas are shown as follows: The subsidence of any point ðx, yÞ on the surface caused by mining of underground random discrete mining units ðs, tÞ is set as follows: The PIM for calculating the ground subsidence at any point caused by mining of all underground discrete elements ðs i , t j Þ on the whole mining face (strike length D1, dip width D2) is as follows: where W 0 is the maximum mining subsidence value; r is the main influencing radius of subsidence, and r = H 0 /tan β, where H 0 is the average mining depth and tan β is the main influencing angle β tangent; l = H * c tan θ, H is the calculation unit  Journal of Sensors depth, and θ is the maximum subsidence angle; ðs, tÞ is the plane coordinate of the working face at the centre of the mining unit; and ðx, yÞ is the terrain coordinate at any point. Based on the above probability integral formula, the relationship between underground movement and ground subsidence is established. 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  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 D.

Journal of Sensors
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. (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 D u = D S /cos θ (D u represents vertical deformation, D S 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.

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 offsettracking 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 7 Journal of Sensors 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 The red parts in the black ellipses in Figure 9(a) represent the DInSAR results, and the blue parts represent the offsettracking results. As shown, the red subsidence boundary acquired by DInSAR surrounds the blue central subsidence acquired by offset tracking. The characteristics of the subsi-dence 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.

Journal of Sensors
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 largescale 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.

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~2 103.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 q, 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 offsettracking 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 offsettracking 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 distri-bution along the mining face, and the inversion of the subsidence location coincides well with the observed subsidence location. The inversion maximum value is approximatelỹ 1.75 m, which is slightly different from the maximum subsidence obtained by the offset tracking. Because we use the

12
Journal of Sensors 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 R 2 = 0:978 is also calculated, and the high coefficient R 2 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.

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 max-imum 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 R 2 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.