Mining Area by Integrating Multitemporal D-InSAR and Offset Tracking Technology

Key Laboratory for Evaluation of Active Tectonics and Crustal Stability, Institute of Geomechanics, Chinese Academy of Geological Sciences, Beijing 100081, China State Key Laboratory of Resources and Environmental Information System, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China University of Chinese Academy of Sciences, Beijing 100049, China China Coal Research Institute, Beijing 100013, China


Introduction
The mining and utilization of coal resources are an important support for regional economic development, but the problems of ecological environment and geological disasters follow. Geological disasters such as ground subsidence induced by underground coal mining have become a huge threat to regional water resources, transportation infrastructure, and people's lives and property. It was reported that in China, the total ground subsidence area as a result of coal mining exceeded 700,000 ha during the period from 1949 to 2002, and the consequent economic loss exceeded 50 billion RMB (equivalent to 7 billion USD). Coal mining in northwestern China is an important industry, and the Ordos Basin is known for its rich coal resources. Ordos, located in Inner Mongolia, China, is one of the important locations for super coal mines. In 2018, the production of coal in the Ordos reached 616 million tons. Along with a large number of coal mining regions in Ordos, ground fissures, ground collapse, and land subsidence have undermined the land surface ecosystems [1][2][3]. The land subsidence caused by coal mining in this area has the characteristics of high subsidence rate and large subsidence. Therefore, obtaining the spatial distribution range and deformation amount of the surface subsidence in this area is helpful to study potential environmental catastrophe and is beneficial to the ecological protection and sustainable development of the mining area.
At present, traditional monitoring techniques of coal mining subsidence mainly consist of leveling, GPS, total station system, and 3-D (three-dimensional) laser scanning. These can provide high-precision position and deformation information but have inherent limitations, including point measurements only, large workforce requirements, and cost inefficiency [4][5][6][7]. In recent years, due to the characteristics of high precision, high resolution, continuous, and allweather [8], interferometric synthetic aperture radar (InSAR) technology has developed rapidly and has become one of the most commonly used measurement techniques on large spatial scales. It can overcome the shortcomings of traditional geodetic methods in surface deformation monitoring and has been widely used to monitor surface deformation caused by coseismic motion [9,10], glacier movement [11][12][13], volcanic eruptions [14,15], ground subsidence [16], landslides [17], and underground mining [18,19].
In this paper, we used D-InSAR and offset tracking technology to generate ground subsidence fields in the Ordos Coalfield, China, where a dense GPS network is available. The D-InSAR technology was initially used to extract the subsidence boundary, and offset tracking technology was subsequently used to resolve a substantial displacement near the central subsidence. Finally, we applied the easiest approach which is to utilize SAR pixel offset tracking (OT) with a single pair of SAR images to measure vertical deformation by dividing the LOS displacement by cosine of radar wave incident angle [20][21][22]. Although this method ignored the contribution of horizontal deformation to LOS displacement, considering that the large deformation of the horizontal displacement in this area seemed to be small compared to the vertical deformation, we reconstructed the whole subsidence field in this area by the results of D-InSAR and offset tracking. The deformation calculation results of the mining area are not much different from the GPS data.  (Figure 1(a)), is chosen as the study area. It is a region known for its abundance of coal resources. The topography of the study area is high in the south-west and low in the north-east part, the elevation ranges from 1123 m and 1323 m.

Study Area and SAR Data
The area belongs to the arid-semiarid regions with a mean annual precipitation of 351 mm from 1984 to 2015, which is generally concentrated from July to September in the survey data collected from coal mine factories [3]. Yellow River runs across the study area (Figure 1(b)), and the material of 75% of the deposits is Quaternary aeolian sand and loess. As of the end of 2019, the total mining area near the study area has exceeded 40 km 2 . These mining activities had induced ground collapse and ground fissures. The survey data which was collected from coal mines shows that ground collapses with a depth of tens of meters and ground cracks with a length of tens of meters and a width of several meters can be found there (Figures 1(e)-1(g)).
2.2. Data. In this article, a total of 51 images were chosen to detect the regional land subsidence. The SAR images were collected from TerraSAR and Sentinel-1 on September 03, 2018, and October 26, 2019. The time interval of the SAR data was 418 days. The spatial resolution of the TerraSAR image is 2 m on the ground, and the satellite adopted the X waveband as the working wave band, with a 20 × 35 km 2 coverage area. The spatial resolution of the Sentinel-1 image is 5 × 20 m, and the satellite adopted the C waveband as the working waveband, with a 20 × 35 km 2 coverage area ( Figure 1(b)). Both radar images are in HH polarization mode, and the detailed information of the SAR images is shown in Table 1.
The free digital elevation model (DEM) data required for processing were 30 m elevation data of the Shuttle Radar Topography Mission (SRTM). In addition, we used the CW-15C series unmanned aerial vehicle (UAV) to conduct topography surveys and optical image shooting and finally obtained 0.1 m resolution DSM data. UAV survey data and GPS observation site data are used to verify the accuracy of calculation results.

Methods
We calculated the deformation over the Buliangou mine by using two SAR techniques, namely multitemporal consecutive D-InSAR and the offset tracking technique. Figure 2 presents the methodology that we applied in this paper and the research process. Deformation calculating from two diverse detection techniques can make up for the shortcoming, in which a single processing method cannot fully obtain the deformation of the entire area, and allowed us to estimate the surface deformation of the entire area. LOS and range displacements from the OT and D-InSAR technique were integrated by using an adaptive spatial filtering method. Afterward, we carried out deformation decomposition from both results above. Finally, we evaluated the accuracy of obtained results, using UAV measurement and GPS leveling monitoring. Interferometric processing was carried out in the GAMMA software of the Linux system, while modeling and postprocessing were carried out in MATLAB and ArcGIS. A more detailed description of each intermediate step was provided in the following subsections.   3 Journal of Sensors images to extract ground deformation [23,24]; the interferometric phase usually consists of five parts: 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 represents the noise phase, which contains InSAR processing errors and thermal noise effects. 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 [25].
The multitemporal consecutive D-InSAR processing has been used in many studies [26]. In contrast, in the consecutive D-InSAR technique and differential interferograms of adjacent SAR acquisitions are calculated and accumulated to provide complete time series interferometric results (e.g.,    Journal of Sensors temporal baselines (Sentinel-1 24 days, TerraSAR 50 days) and perpendicular baseline (Sentinel-1 220 m, TerraSAR 460 m), which keep the interference pairs continuous in time and minimize temporal decorrelation [27]. In addition, due to the absence of TSX data in the first half of 2019, we added two longer baseline interference pairs to maintain time continuity; the distribution of the SAR images in the timeperpendicular baseline plane is shown in Figure 3. First, we extracted SLCs (single look complexes) of the Ter-raSAR and Sentinel-1 images and performed the registration. The accuracy of the registration is below 0.05 pixel error. After registration of the master and slave images, the interferometric phase of the same pixels in the image pair can be calculated. Due to the wide coverage of the image, in order to improve the calculation efficiency, we have carried out the corresponding clipping of the image. After coregistration and clipping, the Ter-raSAR and Sentinel-1 data were calculated using the method of multitemporal D-InSAR process; the simulated phase corresponding to the topography was subtracted from Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) with 30 m pixel resolution, and the minimum-cost flow function is used for phase unwrapping. At last, the unwrapped phase was calculated by the following equation: Due to the small area of the clipped image, ionospheric errors and atmospheric turbulence errors can often be ignored, and only the elevation atmospheric errors are removed.
3.1.2. Offset Tracking Processing. The offset tracking method based on SAR intensity uses normalized cross-correlation (NCC) of the input SAR intensity images to estimate the extent of deformation in both the slant range (line-of-sight of the satellite) and the azimuth (along the orbit of the satellite) directions [28]. Compared with interferometric phase measurements, OT is based on the amplitude information of SAR data and has relatively lower accuracy, typically on the order of 1/30th of a pixel for both range and azimuth directions [29], but OT has a unique advantage in calculating large deformations at the meter level produced by coal mines and is less vulnerability to temporal and spatial decorrelation. This paper selected two SAR images (TerraSAR, Sentinel-1) corresponding to key time nodes in coal mining and uses the offset tracking method to calculate the deformation of the mining area. The deformation results obtained by the offset tracking method are largely affected by the size of the cross-correlation window. Considering the selected data and the spatial resolution, to maximize the correlation values while taking into account the computation efficiency, this paper selects 256 × 256 pixels, and the oversampling factor was set to 2. Although the azimuth deviation is calculated, it is not further analyzed in this study due to the serious influence of ionospheric disturbance.

D-InSAR and OT Fusion.
In coal mining areas, due to problems such as the surface coherence, large deformation gradient on the central surface, and unwrapping problems, a single InSAR method can hardly obtain a continuous deformation map of the mining area. Therefore, in this article, we combined the advantages of both OT and D-InSAR and retained the measured value of D-InSAR in the small deformation area and the measured value of OT large deformation in the center of the mining area.
Through the calculation of the D-InSAR and OT method ("InSAR Processing"), we obtained four parts of data processing results (S1-D-InSAR, TSX-D-InSAR, S1-OT, and TSX-OT). Subsequently, we adopted a strategy to retain only their respective advantageous measurement areas for the four parts of data; the conditions for the application of them are OT measurements have accuracy on the order of several centimeters (on the order of 1/30th of a pixel for both range and azimuth directions), and it is worth remarking that the deformation retrieval capability of the OT measurements obviously improves when the azimuth and range spatial resolution increases [29]. So in order to retain deformation data with an adequate accuracy from both the datasets and obtain a complete deformation field at the same time, we adopted a strategy in which, when deformation values of OT are larger than 0.47 m and 0.33 m (1/30th of a pixel, respectively), the deformation values are retained, while the deformation values smaller than 0.47 m and 0.33 m are discarded; considering that consecutive D-InSAR of short periods were used and its results are more reliable, so the remaining relatively small deformation area was supplemented by D-InSAR results. Although this may introduce some unwrapping errors that lose the deformation phase period, the error of D-InSAR results located in the small deformation area is far small compared with the deformation of meter-grade coal mine subsidence area. Finally, the Kriging interpolation algorithm was applied to link the OT and D-InSAR data above; the two deformation fields were superimposed into a complete field. Kriging is a geostatistical forecasting method that minimizes the error variance by using a weighted linear combination of data. The weight is not only based on the distance between the measuring point and the predicted position but also based on the overall spatial arrangement between the measuring points. This fusion procedure was independently performed for the ascending and descending orbit obtained from the D-InSAR and OT results, respectively. For the geostatistical study, we used the ArcGIS software.

Vertical Deformation Decomposition.
Based on geometric principles of SAR imaging, the relationship between LOS displacement and real 3-D ground surface deformation [30] can be expressed as where θ and α are the incidence angle of radar wave and heading angle of satellite, respectively; the detailed information is shown in Table 2. W, U N , and U E represent real vertical, north-south direction, and west-east direction deformations, respectively. The LOS displacement was composed of the above three parts. The sensitivity of LOS displacement in the vertical component is higher than that in the horizontal component (north-south and west-east), but SAR data are not sensitive to north-south deformation due to a smaller heading angle between the north direction and satellite flight direction; therefore, we ignored the influence of the north-south direction on LOS displacement and assumed that U N = 0; this allowed us to calculate the values of W based on the LOS values or range values from two SAR data.

3.4.
Validation. This paper used two methods to verify the accuracy of the InSAR deformation results. One is GPS observation station data, which is mainly distributed above the surface of the coal mining area and arranged along the mining strike direction, which can monitor the vertical deformation values of the surface; the other is UAV terrain measurement. Through UAV photography and aerial triangulation, it provided the latest surface elevation data. By comparing SRTM DEM data, the elevation change of the mining area can be calculated for verification of InSAR results (Figure 4).  Figure 5 shows the surface deformation values in the LOS direction in various periods of coal mining calculated by TSX data using the D-InSAR method, the maximum measurable deformation value in each period is around 10 cm, and the settlement of the deformation center area has not been measured. Figures 5(a)   For the processing of Sentinel-1 data, we used a method similar to the TSX processing. The difference is that Sentinel-1 images of 39 scenes ensure the continuity of the deformation acquisition period, but due to the lower image resolution, the processing effect is worse than that of TSX (Figures 7 and 8). Figure 7 is a superimposed LOS deformation map of continuous images of multiple scenes. The time baseline is 12 days, but there will still be a large area of deco-herence, the maximum measurable value is about 34 cm, and the maximum settlement calculated in Figure 8 is 11.5 m.

Results and Analysis
Through the fusion of the two methods, we calculated the total surface deformation of the two mining faces during the mining time in the directions of ascending and descending orbits, respectively. The offset tracking method has a great calculation error in the azimuth direction, so this article does not use the deformation value in this direction but regards the range value as the projection value of the settlement in the incident direction: The part that the D-InSAR method cannot obtain the value is replaced by the range value calculated by the offset tracking. And we used a spatial filter method to smooth the deformation distribution. From this, we have obtained the entire deformation field of the coal mining area (Figure 9). 8

Journal of Sensors
The fusion results of TSX data were used to verify the accuracy of the measured values ( Figure 10). Four deformation profiles were drawn around the 0.33 m contours; as shown in Figure 10, the red line represents Kriging interpolated data, the green line represents D-InSAR data, and the blue line represents OT data. The dotted line represents data with possible calculation errors, such as D-InSAR unwrapping problems. After Kriging interpolation, the D-InSAR and OT data were effectively connected, and the lost D-InSAR unwrapping phase was restored. The curves before  Journal of Sensors and after interpolation show that the loss of deformation value due to the unwrapping errors is less than 0.1 m.
To validate the results calculated, we compared the deformation profiles showed in Figure 9, and the deformation curve is shown in the Figures 11 and 12. In Figure 11, the blue line in the figure represents the deformation value of the GPS observation station, the yellow and red represent the Sentinel-1 and TerraSAR measurements, and the green line represents the difference between the UAV and SRTM elevation after mining. At 400 m in the strike direction of the F6211 mining face, the maximum settlement reached about 12 m. GPS monitoring values and two deformation calculation methods have captured this point. In the direction of 200 m and 950 m, there are another two larger settlement points with a peak value of about 9 m and 7 m. These peaks are completely calculated by Sentinel-1 and TerraSAR with higher resolution. In terms of the overall deformation trend, TerraSAR data matches the GPS curve well. Although the calculation results of Sentinel-1 and UAV are generally consistent, the overall error is large.
In Figure 12, at the F6207 mining surface, the two SAR calculation results are basically the same as the GPS monitoring results. In the range of 800 m to 1200 m, the deformation curve of 4-5 m is calculated, and the overall trend is not much different from GPS. Considering that the SRTM measurement time has been too long and the surface changes are large, the UAV measurement results are quite different from GPS.

Comparison of OT and D-InSAR
Methods with GPS and UAV. In this study, the complex stratigraphic conditions of the coal mine goaf, the limitation of the SAR data, and methods of data processing are something worthy of improvement in this article.

Journal of Sensors
The analysis of both the offset tracking and D-InSAR results identified two subsidence basins in the Ordos mine. However, these two methods have some advantages and disadvantages, D-InSAR cannot estimate deformation with strong displacement in the center of the mining area in Figures 5 and 7, and this is because of the large deformation gradient. In contrast, multitemporal consecutive D-InSAR was able to detect fast deformation using multitemporal deformation results, so the D-InSAR results at the cm level are completely reliable, which can be used to calculate the range of ground deformation images in coal mine areas.
Offset tracking technology can measure meter-level surface deformation, but it has a large calculation error in Figures 6 and 8. The multiple aperture SAR interferometry (MAI), with the deformation resolution between the above two methods, can obtain deformation values with higher accuracy on the order of several centimeters to decimeters, which can get an estimate of the surface displacement along the azimuth (sensor line-of-sight (LOS)) direction [31] and has been used to measure the along-track components of the Earth's surface deformation [32][33][34][35][36]; the introduction of the MAI method may improve the accuracy of deformation field calculation. Moreover, in terms of SAR data, the resolution of Sentinel-1 data is too low, so that the OT results of Sentinel-1 are not reliable and can only be used for quantitative analysis. TerraSAR high-resolution data has high credibility, and it has better fitting results with GPS data, but the error in some areas is still large.
The height difference calculated by UAV DSM and STRM DEM is not reliable. The possible reason is that the time of SRTM DEM is too long, and the surface has undergone major changes. If UAV was used for terrain measurement before and after mining, the result will be more accurate. In summary, it is feasible to use offset tracking and D-InSAR results for deformation measurement in mining areas. The deformation trend of Sentinel-1 and TerraSAR data results are basically correct, and the high-resolution data may produce more accurate results that can be used for quantitative analysis.

Vertical Deformation
Component. SAR orientation affects the radar sensor's sensitivity to deformation, and vertical displacement is usually associated with two LOS deformation measurements (ascending and descending), which requires that we need the same SAR data in both directions and in the same time period. The data in this paper cannot meet this requirement, so we only assumed that the deformation in the LOS or range direction is completely caused by the vertical settlement of the ground surface, starting from the principle of geometric imaging to calculate the vertical deformation, and this will undoubtedly introduce a certain error. In the subsequent research, we will fully consider the purchase of data and accurately calculate the three-dimensional deformation of the mining area through the TSX image of both ascending and descending directions.

Conclusions
In this paper, TerraSAR, Sentinel-1 SAR data, and two methods were used to calculate the deformation basin in the Ordos mining area. Offset tracking and D-InSAR methods were fused by a spatial filter method, where large deformation value in the basin selected offset tracking results and the smaller chose multitemporal consecutive D-InSAR results. During the observation period, the maximum deformation value in the total field was nearly 12 m and 4 m in the F6211 and F6207 mining areas, respectively; they are all located at the center of the mining face. By comparing the experimental results with GPS and UAV data, we found that both the S1 data and the TSX data have remained roughly the same in terms of deformation trend and maximum deformation; however, there is a certain gap in the local deformation curve, and the elevation difference calculated by UAV has a large error, which is only used for trend analysis. Considering the accuracy of the offset tracking method itself and the resolution of the S1 data, the error of the calculation result is within an acceptable range.

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

Conflicts of Interest
The authors declare no conflict of interest.