Investigation on Mining Subsidence Based on Multi-Temporal InSAR and Time-Series Analysis of the Small Baseline Subset — Case Study of Working Faces 22201-1 / 2 in Bu ’ ertai Mine , Shendong Coalfield , China

Chao Ma 1,2, Xiaoqian Cheng 1,2,*, Yali Yang 1, Xiaoke Zhang 3, Zengzhang Guo 2 and Youfeng Zou 2 1 Department of Remote Sensing Science and Technology, Henan Polytechnic University, Jiaozuo 454000, China; mac@hpu.edu.cn (C.M.); d_insar@126.com (Y.Y.) 2 Key Laboratory of Mine Spatial Information Technologies of SBSM, Henan Polytechnic University, Jiaozuo 454000, China; gzc@hpu.edu.cn (Z.G.); zouyf@hpu.edu.cn (Y.Z.) 3 School of Public Administration, Hohai University, Nanjing 210098, China; zhang_xiao_ke@126.com * Correspondence: cxq@hpu.edu.cn; Tel.: +86-391-398-7661


Introduction
Shendong Coalfield is one of the largest coal production bases in the world.The high-density fields and high-intensity mining in Shendong Coalfield will profoundly influence the ecological environment of northwestern China.Therefore, assessing mining-induced damage and understanding mining subsidence in the high-intensity mining environment of Shendong Coalfield are important for the land remediation and ecological reconstruction of mining subsidence areas.Traditional mining subsidence monitoring approaches mainly include geodesy, global navigation satellite system measurements, and electronic distance measurements.These approaches exhibit the following shortcomings.(1) They are labor-intensive, time-consuming, costly, and difficult to keep the monitoring flags in good condition for a long time; (2) The survey crew has to enter the region being monitored, which increases the difficulty and risks of the task; (3) These approaches are inapplicable to fast and accurate large-scale real-time monitoring of mines and cannot monitor subsidence in unknown regions because of the limited spatial extent of monitoring, the low spatial resolution, and the long work cycle; (4) Theoretical analysis is restricted within-data observation of discrete points; hence, identifying the ground deformation characteristics while satisfying the actual requirements of mining subsidence prediction and disaster prevention is difficult [1,2].
Interferometric synthetic aperture radar (InSAR) offers a novel Earth observation approach and can provide all-weather, all-day, and cloud influence free monitoring.Differential InSAR (DInSAR), as an extension of InSAR in terms of monitoring ground deformation, is mainly used to capture centimeter-level or smaller ground deformations along the line of sight (LOS) of a radar satellite.Considerable developments related to DInSAR monitoring of earthquake deformation [3,4], volcanic activities [5,6], glacial shift [7,8], urban water-loss settlement [9,10], mining subsidence [11,12], and landslides [13,14] have been achieved.However, DInSAR technology can easily cause interference decorrelation for mining subsidence with immense deformation and short deformation period.Moreover, for the conventional DInSAR, it is a challenging task to distinguish, and thus, effectively eliminate orbit residue error, residual terrain, atmospheric delay, and other phase errors.Nevertheless, several successful applications of the conventional DInSAR have been reported.In the past 10 years, many scholars have monitored mining subsidence in selected mines using DInSAR and have conducted corresponding experimental research.Perski et al. successfully measured ground settlement caused by underground mining using DInSAR [15][16][17][18][19][20].Other scholars have successfully applied DInSAR to mining subsidence monitoring experiments, thereby confirming the feasibility of using DInSAR in large-scale mining subsidence monitoring and geological disaster assessment [21][22][23][24][25][26][27][28].Therefore, the conventional DInSAR coupled with an appropriate interferometric strategy can recognize ground deformation associated with past-mining activities and atmospheric effects (e.g., water vapor in the troposphere and fluctuations in the ionosphere).It can also reduce data costs and increase observation efficiency [29].
To overcome the shortcomings of DInSAR and to acquire accurate subsidence data, the deformation velocity of highly coherent target points was computed via least squares (LS) estimation [30].In 2002, Berardino et al. proposed the small baseline subset (SBAS) algorithm based on the LS model [31][32][33].SBAS was developed from the traditional DInSAR.It combined multiple small baseline subsets via the singular value decomposition (SVD) method based on InSAR data pairs with small spatial and temporal baselines.Therefore, SBAS does not only inherit the advantages of the traditional DInSAR, but can also effectively solve the time discontinuous problem caused by extremely long spatial baselines among SAR data sets, thereby obtaining subsidence values with high coherence targets during different periods.SBAS becomes increasingly optimized and its monitoring accuracy increases accordingly with the development of three-dimensional (3D) phase unwrapping algorithms, extraction algorithms of coherent target points, and error elimination algorithms, as well as the improvement of the time-series deformation model [34][35][36][37][38][39][40].Nevertheless, it still has low computational efficiency.Moreover, precise observation is limited within a small scale, and thus, a substantial part of mining subsidence (i.e., maximum subsidence) cannot be covered.
Therefore, given the limited data set, macroscopic and microscopic studies on global and local ground displacement features caused by mining subsidence were performed by combining conventional and advanced interferometric techniques (DInSAR and SBAS-InSAR).On the one hand, the large-scale dynamic characteristics of mining subsidence were identified through a comparative analysis of multi-temporal data by combining two DInSAR interferences.On the other hand, a time-series analysis of a typical working face was conducted based on SBAS-InSAR and several ground displacement parameters were determined.These processes not only provide complete control of the maturity and stability of DInSAR to perform a fast assessment of mining-induced damage in the entire coalfield, but also fully utilize the technical advantages of SBAS-InSAR in increasing the time sampling rate as well as in inhibiting terrain and atmospheric delay influences on the algorithm.Consequently, an accurate analysis of the subsidence feature of a typical working face is realized.1).Shendong coalfield is suffering from a large-scale, centralized distribution of mining subsidence and ground damage because it is being large-scale exploitation, the main mineable coal beds is particularly thick, and the ratio of mining depth/mining thickness is smaller.(Figure 2).complete control of the maturity and stability of DInSAR to perform a fast assessment of mining-induced damage in the entire coalfield, but also fully utilize the technical advantages of SBAS-InSAR in increasing the time sampling rate as well as in inhibiting terrain and atmospheric delay influences on the algorithm.Consequently, an accurate analysis of the subsidence feature of a typical working face is realized.1).Shendong coalfield is suffering from a large-scale, centralized distribution of mining subsidence and ground damage because it is being large-scale exploitation, the main mineable coal beds is particularly thick, and the ratio of mining depth/mining thickness is smaller.(Figure 2).

Study Area 2: 22201-1/2 Working Face
The 22201-1/2 long wall belongs to Bu'ertai Mine and is composed of two columns of working faces.22201-1 was exploited in 2012, whereas 22201-2 was exploited in 2013 (Figure 3).The mining width of the working face that corresponded to the subsidence area was about 300 m.The mining length from 20 January 2012 to June 2013 was about 3950 m.The 22# coal bed was exploited.The mining seam thickness was 2.2-2.9 m (2.5 m on average), and the average mining depth was about 260 m [41].In Figure 2, the mining boundaries (the purple vertical line) of the working face, while the RADARSAT-2 were imaging, were calculated according to the monthly stopping lines (the green vertical line marked by the arrow symbol) of the coal excavation plans.

Data used
RADARSAT-2 data from 20 January 2012 to June 2013 were used, and the data revisit period was set to 24 day.The spatial resolution was 5 m.Moreover, this study adopted the C wave band as the working wave band, HH polarization mode, multilook fine beam, and 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).Images that covered the entire Bu'ertai Mine (within 39°24′-39°26′N and 109°58′-110°1′E) were collected by subset 18 SAR scenes via region of interests (ROIs) Data from September 2012, 20 January13, and May 2013 were missing, but the remaining revisit period data were all collected (Table 1).The 22201-1/2 long wall belongs to Bu'ertai Mine and is composed of two columns of working faces.22201-1 was exploited in 2012, whereas 22201-2 was exploited in 2013 (Figure 3).The mining width of the working face that corresponded to the subsidence area was about 300 m.The mining length from 20 January 2012 to June 2013 was about 3950 m.The 22# coal bed was exploited.The mining seam thickness was 2.2-2.9 m (2.5 m on average), and the average mining depth was about 260 m [41].In Figure 2, the mining boundaries (the purple vertical line) of the working face, while the RADARSAT-2 were imaging, were calculated according to the monthly stopping lines (the green vertical line marked by the arrow symbol) of the coal excavation plans.

Data used
RADARSAT-2 data from 20 January 2012 to June 2013 were used, and the data revisit period was set to 24 day.The spatial resolution was 5 m.Moreover, this study adopted the C wave band as the working wave band, HH polarization mode, multilook fine beam, and 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).Images that covered the entire Bu'ertai Mine (within 39°24′-39°26′N and 109°58′-110°1′E) were collected by subset 18 SAR scenes via region of interests (ROIs) Data from September 2012, 20 January13, and May 2013 were missing, but the remaining revisit period data were all collected (Table 1).In Figure 2, the mining boundaries (the purple vertical line) of the working face, while the RADARSAT-2 were imaging, were calculated according to the monthly stopping lines (the green vertical line marked by the arrow symbol) of the coal excavation plans.

Data used
RADARSAT-2 data from 20 January 2012 to June 2013 were used, and the data revisit period was set to 24 day.The spatial resolution was 5 m.Moreover, this study adopted the C wave band as the working wave band, HH polarization mode, multilook fine beam, and 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).Images that covered the entire Bu'ertai Mine (within 39 • 24 -39 • 26 N and 109 • 58 -110 • 1 E) were collected by subset 18 SAR scenes via region of interests (ROIs) Data from September 2012, 20 January13, and May 2013 were missing, but the remaining revisit period data were all collected (Table 1).

Methods
In view of the lack of field data, we will adopt a joint calculation and mutual authentication method, such as consecutive DInSAR, cumulative DInSAR and SBAS-InSAR, in order to achieve a high reliability by the way.

Multi-Temporal InSAR Technology
The interferometric principle of mining subsidence based on DInSAR has been discussed in many studies [11,12,[15][16][17][18][19][20][21][22][23][24][25][26][27][28][29], and will not be duplicated in this work.Instead, multi-temporal InSAR technology was applied and two interferometric processing techniques were used for the entire coalfield: the consecutive (i.e., adjacent acquisitions) DInSAR interferometry and the cumulative DInSAR interferometry.These techniques require various physical values and mathematical methods.They use different temporal and spatial baseline combinations, and thus, the corresponding results provide varying geological significance.The consecutive DInSAR interferometry, which is defined that the two adjacent scenes of time series of SAR data (their temporal baseline is the shortest) carry out an interferometry one another and eventually achieve a complete time series interferemetric results (e.g., φ 1,2 , φ 2,3 , . . ., φ n−1,n ), is convenient for the quantitative analysis of interferometric phase changes and can obtain a subsidence extent within a single period.The cumulative DInSAR interferometry, which is defined that the master image is fixed the first scene of time series of SAR data, meanwhile the next scene is selected (their temporal baseline is gradually lengthened) to implement interferometry with it and eventually achieve a complete time series interferograms(e.g., φ 1,2 , φ 1,3 , . . ., φ 1,n = φ 1,2 + φ 2,3 + . . .+ φ n−1,n ), is convenient for identifying the phase change features in a mining area and can obtain the maximum subsidence magnitude.This technique is conducive to evaluating the variation characteristics for mining subsidence on a deformation field from different spatio-temporal perspectives.

Brief Description of SBAS-InSAR
To improve the monitoring accuracy of the InSAR technique, one of the advanced InSAR techniques termed short baseline subsets (SBAS) InSAR was developed and has been frequently applied since 2002 [31][32][33]42,43].
This algorithm uses interferograms with small baselines that overlap in time in order to reduce spatial decorrelation, and to mitigate atmospheric artifacts and topographic errors in time-sequential interferometric pairs.The differential phase for a generic coherent pixel of the range and azimuth coordinates (x, r) in interferogram j that is generated by combining SAR acquisitions at times t B and t A is where λ is the transmitted radar wavelength; φ(t B ,x,r) and φ(t A ,x,r) are the phases that correspond to times t B and t A ; and d(t B ,x,r) and d(t A ,x,r) are the radar LOS projection of the cumulative deformation referenced to the first scene, which implies φ(t 0 ,x,r) = 0. We also include a phase term related to possible errors ∆z in the applied DEM used to generate differential interferograms.This phase component is proportional to the perpendicular baseline for each interferogram B ⊥j , range distance R, and beam incident angle θ.A possible atmospheric signal is included in the terms φ atm (t B ,x,r) and φ atm (t A ,x,r).Decorrelation effects and other noise sources are included in the last term ∆n j .Therefore, Equation ( 1) can be expressed as follows: where δφ disp is the slant deformation phase between t B and t A ; δφ topo represents the topographic phase error of the external DEM used for differential interferogram generation; δφ atm accounts for temporal atmospheric variation at different SAR acquisition t B and t A ; and δφ noise denotes temporal decorrelation, orbital errors, thermal noise effects, etc.The cumulative deformation between t B and t A can also be expressed as follows: where k indicates the sequence index of SAR acquisition time between t B and t A , and v represents the mean phase velocity in the period from k to k+1.Given M unwrapped interferograms, the cumulative deformation at different SAR acquisition times can be achieved using the LS or SVD methods.

Spatio-Temporal Baselines of Two Interferometric Strategies
The temporal baseline of the consecutive DInSAR interferometry can be controlled efficiently within 24-48 days (Figure 4a), which will produce an ideal interferometric effect.However, the temporal baseline of the cumulative DInSAR interferometry increased continuously from the initial 24 days to the final 504 days (Figure 4b), which indicated that the influence of temporal decorrelation increasingly intensified.

Results of Consecutive and Cumulative DInSAR Interferometries
(1) For the consecutive DInSAR interferometry, no image is fixed to a master image.The advanced 17 images were in turn used as the master image, whereas the subsequent image was used as the slave image.The temporal baseline was the same with the revisit period.The interferogram was characterized by independent differential phase changes in various periods and a forward movement of the footprint of the deformation phase of mining subsidence (Figure 5).

Results of Consecutive and Cumulative DInSAR Interferometries
(1) For the consecutive DInSAR interferometry, no image is fixed to a master image.The advanced 17 images were in turn used as the master image, whereas the subsequent image was used as the slave image.The temporal baseline was the same with the revisit period.The interferogram was characterized by independent differential phase changes in various periods and a forward movement of the footprint of the deformation phase of mining subsidence (Figure 5).

Results of Consecutive and Cumulative DInSAR Interferometries
(1) For the consecutive DInSAR interferometry, no image is fixed to a master image.The advanced 17 images were in turn used as the master image, whereas the subsequent image was used as the slave image.The temporal baseline was the same with the revisit period.The interferogram was characterized by independent differential phase changes in various periods and a forward movement of the footprint of the deformation phase of mining subsidence (Figure 5).(2) For the cumulative DInSAR interferometry, the first image was set as the master image (Master, 20 January 2012) and the subsequent 17 images were alternately used as slave images.The interferogram indicated that the differential phase changes covered all previous phase changes, the deformation phase center of mining subsidence moved forward, and the deformation region expanded monthly (Figure 6).(2) For the cumulative DInSAR interferometry, the first image was set as the master image (Master, 20 January 2012) and the subsequent 17 images were alternately used as slave images.The interferogram indicated that the differential phase changes covered all previous phase changes, the deformation phase center of mining subsidence moved forward, and the deformation region expanded monthly (Figure 6).The cumulative DInSAR interferometry was influenced simultaneously by temporal and spatial decorrelations.The interferogram had many unstable patches and several phase residues.It failed in retaining coherence and phase unwrapping during the last three periods.

SBAS-InSAR Data Processing
To eliminate spatiotemporal decorrelation and atmospheric effects, an experimental plot was selected to test SBAS-InSAR, which was expected to reduce or eliminate the influences of decorrelation factors during DInSAR data processing by selecting appropriate spatiotemporal baseline thresholds.

Spatiotemporal Baseline Optimization of SBAS
All interferometric pairs within the thresholds of the spatial and temporal baselines were screened based on a preset spatial baseline threshold, a temporal baseline threshold, and relevant input parameters, thereby generating the SAR data pair connection diagram (Figure 7).In this experiment, the maximum critical baseline percentage was set at 5%.The maximum perpendicular component of the spatial baseline was 422 m.The maximum temporal baseline was set 200 days.The system automatically screened the super master image (20120823).Interferometric pairs with low coherence and poor unwrapping were eliminated, and finally, 47 interferometric pairs were selected.The cumulative DInSAR interferometry was influenced simultaneously by temporal and spatial decorrelations.The interferogram had many unstable patches and several phase residues.It failed in retaining coherence and phase unwrapping during the last three periods.

SBAS-InSAR Data Processing
To eliminate spatiotemporal decorrelation and atmospheric effects, an experimental plot was selected to test SBAS-InSAR, which was expected to reduce or eliminate the influences of decorrelation factors during DInSAR data processing by selecting appropriate spatiotemporal baseline thresholds.

Spatiotemporal Baseline Optimization of SBAS
All interferometric pairs within the thresholds of the spatial and temporal baselines were screened based on a preset spatial baseline threshold, a temporal baseline threshold, and relevant input parameters, thereby generating the SAR data pair connection diagram (Figure 7).In this experiment, the maximum critical baseline percentage was set at 5%.The maximum perpendicular component of the spatial baseline was 422 m.The maximum temporal baseline was set 200 days.The system automatically screened the super master image (20120823).Interferometric pairs with low coherence and poor unwrapping were eliminated, and finally, 47 interferometric pairs were selected.

Interferometry Results of SBAS
Considering significant computational efforts and necessary disk space, SBAS was applied finely in a single working face of the 22201-1/2 long wall.Interferometry was performed on each image pair according to their connection relationship based on the aforementioned optimized baselines of SBAS, and 47 interferograms were obtained.A total of 17 time-series cumulative phase deformation diagrams were collected through orbit refining, reflatting, phase unwrapping, and geocoding.Subsequently, a −10 mm subsidence contour map was acquired based on the further processing of these diagrams (Figure 8a-q).

Interferometry Results of SBAS
Considering significant computational efforts and necessary disk space, SBAS was applied finely in a single working face of the 22201-1/2 long wall.Interferometry was performed on each image pair according to their connection relationship based on the aforementioned optimized baselines of SBAS, and 47 interferograms were obtained.A total of 17 time-series cumulative phase deformation diagrams were collected through orbit refining, reflatting, phase unwrapping, and geocoding.Subsequently, a −10 mm subsidence contour map was acquired based on the further processing of these diagrams (Figure 8a-q).

Interferometry Results of SBAS
Considering significant computational efforts and necessary disk space, SBAS was applied finely in a single working face of the 22201-1/2 long wall.Interferometry was performed on each image pair according to their connection relationship based on the aforementioned optimized baselines of SBAS, and 47 interferograms were obtained.A total of 17 time-series cumulative phase deformation diagrams were collected through orbit refining, reflatting, phase unwrapping, and geocoding.Subsequently, a −10 mm subsidence contour map was acquired based on the further processing of these diagrams (Figure 8a-q).

The Reliability Validation of the Interferometry
Because of the past mining subsidence and lack of field data, the reliability of interferometric measurements needs to be verified and validated.We randomly selected three subsidence basins (basin A,B,C) and some understanding are gotten by comparing the results of the two interferometry measurements(i.e., SABS-InSAR and the consecutive DInSAR).
See Figures 9-11, using the same profile line, we extract the subsidence values in the two subsidence basins by SABS-InSAR and the consecutive DInSAR, and then perform the overlay analysis.
The unwrapping of the central area of the subsidence basin is erroneous due to the influence of decoherence.So we only compare the information at the edge of the subsidence basins.We found

The Reliability Validation of the Interferometry
Because of the past mining subsidence and lack of field data, the reliability of interferometric measurements needs to be verified and validated.We randomly selected three subsidence basins (basin A,B,C) and some understanding are gotten by comparing the results of the two interferometry measurements(i.e., SABS-InSAR and the consecutive DInSAR).
See Figures 9-11, using the same profile line, we extract the subsidence values in the two subsidence basins by SABS-InSAR and the consecutive DInSAR, and then perform the overlay analysis.
The unwrapping of the central area of the subsidence basin is erroneous due to the influence of decoherence.So we only compare the information at the edge of the subsidence basins.We found that the results of two interferometric measurements coincide each other fairly well.For example, the Root Mean Square Error (RMSE) of the subsidence values on both sides of basin A are ±0.6 mm (left) and ±0.3 mm (right) respectively (see Figure 9), the RMSE of the subsidence values on both sides of basin B are ±0.6 mm (left) and ±0.2 mm (right) respectively (see Figure 10), and the RMSE of the subsidence values on both sides of basin C are ±0.4 mm (left) and ±0.1 mm (right) respectively (see Figure 11).Therefore, through the above cross-validation, we believe that our data processing scheme is highly reliable, and it can be used for large-scale investigation of past mining subsidence.
that the results of two interferometric measurements coincide each other fairly well.For example, the Root Mean Square Error (RMSE) of the subsidence values on both sides of basin A are ±0.6 mm (left) and ±0.3 mm (right) respectively (see Figure 9), the RMSE of the subsidence values on both sides of basin B are ±0.6 mm (left) and ±0.2 mm (right) respectively (see Figure 10), and the RMSE of the subsidence values on both sides of basin C are ±0.4 mm (left) and ±0.1 mm (right) respectively (see Figure 11).Therefore, through the above cross-validation, we believe that our data processing scheme is highly reliable, and it can be used for large-scale investigation of past mining subsidence.

Interpretation based on the consecutive DInSAR interferometry
The results of the 17-temporal the consecutive DInSAR interferometry were interpreted using the subsidence value of −10 mm recommended by a literature as the outer boundary of mining subsidence [44].Polygons with low speckle in the interpreted subsidence were further screened.Some polygons, which did not change during three revisit periods (According [44], the zone, where mining subsidence is less than −50 mm during 3 months, is considered to be the residual deformation zone.), will be removed, thereby obtaining the time-series subsidence map of the consecutive DInSAR interferometry (Figure 12).
The interpretation results indicated that Huhewusu Mine, Erlintu Mine, Zhugaita Mine, and Dahaize Mine in Shendong Coalfield had been hardly any exploited, whereas the remaining mines produced an average of 45 ± 3 dynamic subsidence regions during each observation period from 20 January 2012 to 7 June 2013.The average mining extent during each observation period reached as high as 57.00 ± 0.24 km 2 , which indicated the active production activities in the coalfield.Approximately 45 mining areas were exploited simultaneously, which caused serious mining subsidence.Interpretation was a difficult task accomplished through team cooperation.The concerned subsidence quantity and area would differ to a certain extent given that the team members had varying work experiences.

Interpretation based on the consecutive DInSAR interferometry
The results of the 17-temporal the consecutive DInSAR interferometry were interpreted using the subsidence value of −10 mm recommended by a literature as the outer boundary of mining subsidence [44].Polygons with low speckle in the interpreted subsidence were further screened.Some polygons, which did not change during three revisit periods (According [44], the zone, where mining subsidence is less than −50 mm during 3 months, is considered to be the residual deformation zone.), will be removed, thereby obtaining the time-series subsidence map of the consecutive DInSAR interferometry (Figure 12).
The interpretation results indicated that Huhewusu Mine, Erlintu Mine, Zhugaita Mine, and Dahaize Mine in Shendong Coalfield had been hardly any exploited, whereas the remaining mines produced an average of 45 ± 3 dynamic subsidence regions during each observation period from 20 January 2012 to 7 June 2013.The average mining extent during each observation period reached as high as 57.00 ± 0.24 km 2 , which indicated the active production activities in the coalfield.Approximately 45 mining areas were exploited simultaneously, which caused serious mining subsidence.Interpretation was a difficult task accomplished through team cooperation.The concerned subsidence quantity and area would differ to a certain extent given that the team members had varying work experiences.

Interpretation Based on the cumulative DInSAR Interferometry
Similarly, the results of the 14-temporal the cumulative DInSAR interferometry were interpreted using the subsidence value of −10 mm recommended by [44] as the outer boundary of mining subsidence.The polygons with low speckle in the interpreted subsidence were further screened, and those that did not change successively will be removed.Finally, the time-series subsidence map was obtained (Figure 13).Interpretation became more difficult given the influence of temporal decorrelation, and the interpretable degree was lower than that of the consecutive DInSAR interferometry.Evidently, the

Interpretation Based on the cumulative DInSAR Interferometry
Similarly, the results of the 14-temporal the cumulative DInSAR interferometry were interpreted using the subsidence value of −10 mm recommended by [44] as the outer boundary of mining subsidence.The polygons with low speckle in the interpreted subsidence were further screened, and those that did not change successively will be removed.Finally, the time-series subsidence map was obtained (Figure 13).

Interpretation Based on the cumulative DInSAR Interferometry
Similarly, the results of the 14-temporal the cumulative DInSAR interferometry were interpreted using the subsidence value of −10 mm recommended by [44] as the outer boundary of mining subsidence.The polygons with low speckle in the interpreted subsidence were further screened, and those that did not change successively will be removed.Finally, the time-series subsidence map was obtained (Figure 13).Interpretation became more difficult given the influence of temporal decorrelation, and the interpretable degree was lower than that of the consecutive DInSAR interferometry.Evidently, the Interpretation became more difficult given the influence of temporal decorrelation, and the interpretable degree was lower than that of the consecutive DInSAR interferometry.Evidently, the characteristic of the time series of the cumulative DInSAR interferometry was that polygons with low speckle in the subsidence range expanded to a fixed direction.An average of 46 ± 5 subsidence areas expanded stably during each study period from 20 January, 2012 to 21 December 2012.The subsidence area increased from an initial value of 54.98 km 2 /24 day to 225.20 km 2 /504 day throughout the observation period (20 January 2012 to June 2013), which indicated an average rate of 13.09 km 2 .
In general, it is still difficult to obtain a continuous series of subsidence areas.On the one hand, mine production regarded mine as a unit, which was not uniform.Several subsidence areas were caused by ground deformation associated with post-mining activities, others resulted from the preparation of alternate long walls, whereas some areas were produced by trial mining.On the other hand, mining developments in different mines varied because of the different geological conditions of the mines.Moreover, surface construction might also cause misjudgment of interpretations.

Mining Subsidence Characteristics in Shendong Coalfield
To identify the mining subsidence characteristics in high-intensity mines, 19 continuous deformation areas with mining subsidence in the west bank of Wulanmulun River were investigated.These study areas were divided into A (A1, A2, A3, A4, and A5), B (B1, B2, and B3), C (C1), D (D1, D2, D3, and D4), E (E1, E2, and E3), F (F1 and F2), and G (G1) (Figure 14 A-G).characteristic of the time series of the cumulative DInSAR interferometry was that polygons with low speckle in the subsidence range expanded to a fixed direction.An average of 46 ± 5 subsidence areas expanded stably during each study period from 20 January, 2012 to 21 December 2012.The subsidence area increased from an initial value of 54.98 km 2 /24 day to 225.20 km 2 /504 day throughout the observation period (20 January 2012 to June 2013), which indicated an average rate of 13.09 km 2 .
In general, it is still difficult to obtain a continuous series of subsidence areas.On the one hand, mine production regarded mine as a unit, which was not uniform.Several subsidence areas were caused by ground deformation associated with post-mining activities, others resulted from the preparation of alternate long walls, whereas some areas were produced by trial mining.On the other hand, mining developments in different mines varied because of the different geological conditions of the mines.Moreover, surface construction might also cause misjudgment of interpretations.

Mining Subsidence Characteristics in Shendong Coalfield
To identify the mining subsidence characteristics in high-intensity mines, 19 continuous deformation areas with mining subsidence in the west bank of Wulanmulun River were investigated.These study areas were divided into A (A1, A2, A3, A4, and A5), B (B1, B2, and B3), C (C1), D (D1, D2, D3, and D4), E (E1, E2, and E3), F (F1 and F2), and G (G1) (Figure 14 A-G).In order to ensure that the section lines of each period would pass through the maximum subsidence point, the coherence images of the consecutive DInSAR interferometry were used, in which mining subsidence areas were determined as low-coherence areas.The centroids (i.e., geometric centers) of low-coherence areas were connected, thereby forming the principal sectional line of the maximum subsidence.Subsequently, this sectional line was used to extract the subsidence values under the consecutive DInSAR interferometry and the cumulative DInSAR interferometry.Then, the subsidence characteristic curve of 19 continuous deformation areas was drawn (Figure 15).
Only several representative working faces were illustrated and analyzed in this paper given the limit in article length.
In Figure 14, the subfigures('A'-'G') show 7 sub-regions of the focused investigation in the west bank of Wulanmulun River, and they contain a total of 19 continuous deformation areas.The Subfigure 'H' is a regional aerial view that represents the relative spatial relationships of the 7 sub-regions.
In Figure 15, in a sectional graph way, interferometric results from 6 working faces (A3, B2, C1, D2, E1, and G1) are shown concentrately, which comprise that of the consecutive DInSAR interferometry located in left and that of the cumulative DInSAR interferometry located in right.The 6 working faces are located in Cuncaota No. 2 Mine, Cuncaota No.1 Mine, Bu'ertai Mine, Buliantai Mine, Shangwan Mine, and Huojitu Mine respectively (see also Figures 1 and 14).
In order to ensure that the section lines of each period would pass through the maximum subsidence point, the coherence images of the consecutive DInSAR interferometry were used, in which mining subsidence areas were determined as low-coherence areas.The centroids (i.e., geometric centers) of low-coherence areas were connected, thereby forming the principal sectional line of the maximum subsidence.Subsequently, this sectional line was used to extract the subsidence values under the consecutive DInSAR interferometry and the cumulative DInSAR interferometry.Then, the subsidence characteristic curve of 19 continuous deformation areas was drawn (Figure 15).
Only several representative working faces were illustrated and analyzed in this paper given the limit in article length.
In Figure 14, the subfigures('A'-'G') show 7 sub-regions of the focused investigation in the west bank of Wulanmulun River, and they contain a total of 19 continuous deformation areas.The Subfigure 'H' is a regional aerial view that represents the relative spatial relationships of the 7 sub-regions.
In Figure 15, in a sectional graph way, interferometric results from 6 working faces (A3, B2, C1, D2, E1, and G1) are shown concentrately, which comprise that of the consecutive DInSAR interferometry located in left and that of the cumulative DInSAR interferometry located in right.The 6 working faces are located in Cuncaota No. 2 Mine, Cuncaota No.1 Mine, Bu'ertai Mine, Buliantai Mine, Shangwan Mine, and Huojitu Mine respectively (see also Figure 1 and Figure 14).Along the centroids of subsidence basin of the working faces (see also Figure 14), 6 observation lines (length difference from 2200 m (B2) to 4840 m (G1)) were set, and re-sampled at an interval of 20 m along the line.So a series of subsidence curves were drawn along the time axis (interval difference from 336 day (A3, B2, and E1) to 504 day (C1 and G1)).
In Figure 15, a common attribute is in that the results of the consecutive DInSAR interferometry reappear some phenomena of a forward movement of the footprint of the deformation area with a series of perfect subsidence basins of subcritical mining such as the Figure 15a,c,e,g.Only in Figure 15 i,k it represents a bit of fluctuate due the influence of near mining.Mean over, another common attribute is in that the results of the cumulative DInSAR interferometry represent a pattern of the subsidence basin area continuously expanded and the center of mining subsidence continuously moved forward (such as the Figure 15b,d,f,h,j,l).
Nevertheless, the Figure 15 also represents the characteristic of high-intensity mining subsidence.Firstly, it is in that maximum subsidence cannot be obtained from the results of the consecutive DInSAR interferometry due to strong phase discontinuities caused by the large deformation gradient.In other words, the subsidence curve is incorrect in the bottom of the bowl.Secondly, the results of the consecutive DInSAR interferometry express three different patterns of subsidence that the subcritical mining (Figure 15d,j), the critical mining (Figure 15b,l), and the supercritical mining respectively (Figure 15f,h).Finally, the results of the average subsidence rate during the observation period were range from 240 mm/year to 720 mm/year, which is a far less than that of the Poland [43] and of the Australia (a private communications, Ge, L.L., 2016.7 [45]).Along the centroids of subsidence basin of the working faces (see also Figure 14), 6 observation lines (length difference from 2200 m (B2) to 4840 m (G1)) were set, and re-sampled at an interval of 20 m along the line.So a series of subsidence curves were drawn along the time axis (interval difference from 336 day (A3, B2, and E1) to 504 day (C1 and G1)).
In Figure 15, a common attribute is in that the results of the consecutive DInSAR interferometry reappear some phenomena of a forward movement of the footprint of the deformation area with a series of perfect subsidence basins of subcritical mining such as the Figure 15a,c,e,g.Only in Figure 15i,k it represents a bit of fluctuate due the influence of near mining.Mean over, another common attribute is in that the results of the cumulative DInSAR interferometry represent a pattern of the subsidence basin area continuously expanded and the center of mining subsidence continuously moved forward (such as the Figure 15b,d,f,h,j,l).
Nevertheless, the Figure 15 also represents the characteristic of high-intensity mining subsidence.Firstly, it is in that maximum subsidence cannot be obtained from the results of the consecutive DInSAR interferometry due to strong phase discontinuities caused by the large deformation gradient.In other words, the subsidence curve is incorrect in the bottom of the bowl.Secondly, the results of the consecutive DInSAR interferometry express three different patterns of subsidence that the subcritical mining (Figure 15d,j), the critical mining (Figure 15b,l), and the supercritical mining respectively (Figure 15f,h).Finally, the results of the average subsidence rate during the observation period were range from 240 mm/year to 720 mm/year, which is a far less than that of the Poland [43] and of the Australia (a private communications, Ge, L.L., 2016.7 [45]).

Mining Subsidence Characteristics in Working Face 22201-1/2 in Bu'ertai Mine
Time-series analysis was conducted to discriminate ground deformation associated with past-mining activity and atmospheric effects, and consequently, acquire the spatial distribution characteristics, deformation pattern, development direction, evolution property, and temporal correlation of the multi-temporal mining subsidence phase-changing areas.For the determined boundaries of mining-induced damages, the classification of ground displacement and the calculated ground displacement parameters have practical significance.The experimental research only focused on working face 22201-1/2 in Bu'ertai Mine considering the heavy computing and analysis workloads.Multi-temporal DInSAR interferometry and SBAS-InSAR interferometry were adopted to calculate the available parameters.The subsidence patterns obtained using these two interferometric strategies were compared.The results are shown in Figure 16.

Mining Subsidence Characteristics in Working Face 22201-1/2 in Bu'ertai Mine
Time-series analysis was conducted to discriminate ground deformation associated with past-mining activity and atmospheric effects, and consequently, acquire the spatial distribution characteristics, deformation pattern, development direction, evolution property, and temporal correlation of the multi-temporal mining subsidence phase-changing areas.For the determined boundaries of mining-induced damages, the classification of ground displacement and the calculated ground displacement parameters have practical significance.The experimental research only focused on working face 22201-1/2 in Bu'ertai Mine considering the heavy computing and analysis workloads.Multi-temporal DInSAR interferometry and SBAS-InSAR interferometry were adopted to calculate the available parameters.The subsidence patterns obtained using these two interferometric strategies were compared.The results are shown in Figure 16.The subsidence curves of multi-temporal DInSAR interferometry with SBAS-InSAR interferometry were compared, and the following conclusions were drawn.
(1) In Figures 16a,b, the 17-temporal subsidence curves were obtained from the superposition of multi-temporal DInSAR interferometry, which could reflect general mining subsidence characteristics.That is, underground mining caused subsidence.With the expansion of the mining The subsidence curves of multi-temporal DInSAR interferometry with SBAS-InSAR interferometry were compared, and the following conclusions were drawn.
(1) In Figure 16a,b, the 17-temporal subsidence curves were obtained from the superposition of multi-temporal DInSAR interferometry, which could reflect general mining subsidence characteristics.That is, underground mining caused subsidence.With the expansion of the mining area, the subsidence area gradually expanded toward the advancing direction and subsidence continuously increased.The maximum subsidence became stable when it reached a certain extent, thereby forming a critical mining-induced subsidence bowl.However, rapid considerable subsidence caused by high-intensity mining (large mining height, thin bedrock, and fast advancing) brought serious decorrelation, and the principal value of phase change could not be effectively extracted.Accordingly, only the first temporal subsidence curve was connected smoothly because subsidence was within the unwrapping capacity of DInSAR (λ/4 = 14 mm) [46].
By contrast, in the subsequent observation the large deformation gradient resulted unwrapping failure and had a noise phase that was beyond the calculation capacity of the DInSAR.Such subsidence curves had no practical significance and could not obtained maximum subsidence.Furthermore, the effect of temporal decorrelation increased with temporal baseline, and the subsidence curve became insensitive.This condition was accompanied by a certain delay in the middle-temporal and late-temporal subsidence curves.The accumulation of residual phase changes also caused unexpected subsidence in non-mining areas.
(2) As shown in Figure 16c,d, SBAS-InSAR interferometry optimized the spatio-temporal baselines and improved the spatio-temporal sampling rates.The obtained results were more reliable than those of DInSAR.The 17-temporal subsidence curves of SBAS-InSAR did not only reflect the mining subsidence characteristics more accurately, but was also closer to the subsidence curve of probability integral estimation and maintained good continuity.However, it still failed to extract the principal value of phase changes because of the large deformation caused by mining subsidence.The acquired maximum subsidence was smaller than 200 mm, which was significantly different from that in practical situations.
(3) The mining area (1.23 km 2 ) and subsidence area (5.58 km 2 ) of the 22201-1/2 working face were calculated by further interpreting the results of SBAS-InSAR.The ratio of the mining area to the subsidence area was 1:4.5.The daily average advance rate of the working face was 8.1 m/day, and the daily average advance rate of subsidence area was 7.9 m/day.The average maximum subsidence angle was 67.3 • .The advance angles of influence of DInSAR and SABS-InSAR were 53.5 • and 50.5 • , respectively.

Conclusions
In this study, InSAR and SBAS-InSAR techniques allow us understand the annual growing characteristics of subsidence disasters under high-intensity mining in Shendong Coalfield, and recognize the ground deformation features above a single working face.
(1) Conventional DInSAR is an effective approach for investigating ground damages caused by mining subsidence in coal mining areas.It demonstrates our tremendous destructive capability of the high-intensity coal exploitation.Over 85 phase change regions are identified within 504 days (20 January 2012 to June 2013) in Shendong Coalfield, including at that time and in the past, and including open-pit mining and underground mining areas.Annual subsidence area is about 150km 2 , approximately, which accounted for 4% of the total area of Shendong Coalfield (3481 km 2 ).
(2) The time-series analysis of SBAS-InSAR achieves more reliable results and more accurately reflects the mining subsidence characteristics of a single working face.It makes some important parameters (such as the advance distance of influence and the advance angle of influence, which parameters can be used to predict the extent of mining subsidence so that action can be taken to prevent buildings from being damaged or other measures taken to control potential geohazards, etc.) is possible to be obtained.
(3) We find that the ground damage rate(7.9m/day) caused by the high-intensity mining, is almost equal to mining rate of the working face(8.1 m/day), which should be a unique phenomenon occurred only at high-intensity mining conditions.In addition, we find that the advance distance of influence is much greater than the distance of the exploitation, but the advance angle of influence(50.5 • -53.5 • ) is much smaller than the angle of conventional mining(67.0• -70.0 • ), this may be due to the rate of water-loss settlement is much than that of mining subsidence.
(4) We find that some uplift are occurred at the front of working face in the direction of advancing, which may be caused systematically by orbit error, or by an inappropriate selection of the reference points.It is more likely a unique phenomenon associated with high-intensity mining.Although we have not yet to find a well-received explanation, it should be the focus of our future work, and we will try to understand and explain this phenomenon by means of the in-situ field measurements.

2. 1 . 1 .
Study Area 1: Shendong Coalfield Shendong Coalfield, an immense coal production base in China, is located in the northern region of Shenmu County, the western region of Fugu County (Yunlin City, Shaanxi Province) and Ejin Horo County, and the southern region of Dongsheng District and Junggar County (Ordos City, Inner Mongolia).Its geographical coordinates are 38 • 52 -39 • 41 N and 109 • 51 -110 • 46 E. Its north-south length is approximately 38-90 km and its east-west width is approximately 35-55 km, thereby covering an area of approximately 3481 km 2 .The proven reserves are 727 Gt.Shendong Coalfield covers a cluster of the mines of 10-million-tons annual production that comprise 14 underground coal mines, including Bu'ertai (20.0 Mt/year), Daliuta (25.0 Mt/year), Bulianta (21.0 Mt/year), Halagou (12.0 Mt/year), Shangwan (10.0 Mt/year), and Shigetai (10.0 Mt/year) mine, as well as 5 open-pit mines, including Ha'erwusu (30.0 Mt/year) and Heidaigou (31.0 Mt/year) mine.Shendong Coalfield is one of the largest combined open-pit-underground mining enterprises with the highest degree of modernization in the world (Figure

Figure 2 .
Figure 2. The photographs of subsidence damage types in Shendong Coalfield (2009-2015).2.1.2.Study Area 2: 22201-1/2 Working FaceThe 22201-1/2 long wall belongs to Bu'ertai Mine and is composed of two columns of working faces.22201-1 was exploited in 2012, whereas 22201-2 was exploited in 2013 (Figure3).The mining width of the working face that corresponded to the subsidence area was about 300 m.The mining length from 20 January 2012 to June 2013 was about 3950 m.The 22# coal bed was exploited.The mining seam thickness was 2.2-2.9 m (2.5 m on average), and the average mining depth was about 260 m[41].

Figure 4 .
Figure 4. Distribution of the spatiotemporal baseline.(a) Spatiotemporal baseline distribution of consecutive DInSAR; (b) Spatiotemporal baseline distribution of cumulative DInSAR.

Figure 4 .
Figure 4. Distribution of the spatiotemporal baseline.(a) Spatiotemporal baseline distribution of consecutive DInSAR; (b) Spatiotemporal baseline distribution of cumulative DInSAR.

Figure 4 .
Figure 4. Distribution of the spatiotemporal baseline.(a) Spatiotemporal baseline distribution of consecutive DInSAR; (b) Spatiotemporal baseline distribution of cumulative DInSAR.

Figure 7 .
Figure 7. Spatiotemporal baseline distribution of the connection diagram and unwrapping results: (a) Spatiotemporal baseline distribution of the connection diagram; (b) 3D unwrapping sub-pairs.

Figure 7 .
Figure 7. Spatiotemporal baseline distribution of the connection diagram and unwrapping results: (a) Spatiotemporal baseline distribution of the connection diagram; (b) 3D unwrapping sub-pairs.

Figure 7 .
Figure 7. Spatiotemporal baseline distribution of the connection diagram and unwrapping results: (a) Spatiotemporal baseline distribution of the connection diagram; (b) 3D unwrapping sub-pairs.

Figure 9 .
Figure 9. Example I: verification and validation of Interferometric Measurement.Figure 9. Example I: verification and validation of Interferometric Measurement.

Figure 9 .Figure 10 .
Figure 9. Example I: verification and validation of Interferometric Measurement.Figure 9. Example I: verification and validation of Interferometric Measurement.

Figure 11 .
Figure 11.Example III: verification and validation of Interferometric Measurement.

Figure 11 .
Figure 11.Example III: verification and validation of Interferometric Measurement.

Figure 14 .
Figure 14.(A-H) Subsidence areas obtained from the InSAR time-series analysis (Western region to the Wulanmulun River).Figure 14. (A-H) Subsidence areas obtained from the InSAR time-series analysis (Western region to the Wulanmulun River).

Figure 14 .
Figure 14.(A-H) Subsidence areas obtained from the InSAR time-series analysis (Western region to the Wulanmulun River).Figure 14. (A-H) Subsidence areas obtained from the InSAR time-series analysis (Western region to the Wulanmulun River).

Table 1 .
Imaging parameters and geometry of RADARSAT-2 MF6 data used in this work.