Characteristics of the Residual Surface Deformation of Multiple Abandoned Mined-Out Areas Based on a Field Investigation and SBAS-InSAR: A Case Study in Jilin, China

: Residual surface deformation resulting from abandoned mined-out areas can lead to severe damage to ground structures (e.g., buildings and infrastructure in mining areas) and the local ecological environment. Long-term monitoring and analyses of surface deformation characteristics of abandoned mined-out areas are signiﬁcant for preventing potential disasters. In this study, a detailed ﬁeld investigation ﬁrst was conducted in Ying’an and Baoshan coal mines located in Jilin Province, China, to survey mining-induced disasters in the mining areas. Based on the 40 Sentinel-1A images acquired from 14 February 2017 to 17 May 2020, small baseline subset interferometry synthetic aperture radar (SBAS-InSAR) technology was employed to obtain the time-series residual surface deformation. Validation of the SBAS-derived results is performed by comparing with the results obtained via leveling measurements. The root mean square error (RMSE) between SBAS-derived and leveling measurements results was found to be 1.144 mm, reﬂecting a fairly good agreement. Furthermore, the ordinary Kriging interpolation approach was adopted to obtain information on the deformation across the entire area. The spatial–temporal evolution characteristics of the derived subsidence bowls in multiple mined-out areas were revealed. The deformation characteristics for the abandoned mined-out areas in di ﬀ erent periods were not completely consistent. Finally, the potential mechanism underlying the inconsistency in the subsidence associated with underground coal exploitation is analyzed. The ﬁndings of this study can provide insights into local construction and ecological improvement as well as guidance for the prediction of deformation in abandoned mined-out areas.


Introduction
Coal resources are still one of the main energy sources in many countries [1], which are widely employed in industrial production, power generation and so on. Although coal resources have been central to promoting economic development, extensive underground coal exploitation inevitably results in some negative effects [2], such as surface subsidence [3,4]. Mining-induced surface subsidence gives rise to cracks and other types of damage to infrastructure and resident buildings within mining areas [5,6]. Such surface subsidence also may destroy cultivated land and even give rise to a threat to the local ecological environment [7][8][9]. Therefore, long-term monitoring of surface deformation in mining areas and analysis of surface deformation characteristics are of great significance for preventing and mitigating regional disasters.
Traditional approaches for monitoring and analyzing the surface deformation resulting from mining activities include leveling measurements, global navigation satellite systems, global positioning systems, total station surveys and three-dimensional laser scanning, etc. [10][11][12]. Over the past few decades, the aforementioned geodetic measurement approaches have proven to be reliable geodetic techniques with centimeter-to-millimeter accuracy [10,13], and have been extensively utilized to investigate surface deformation [14,15]. Since these techniques are implemented on a point-by-point basis, the collection of monitoring data is labor-intensive, time-consuming and costly, especially for large study regions [6,16]. Therefore, obtaining effectively the deformation characteristics of the entire mining areas becomes a major challenge [4,17].
Over the past few decades, the interferometric synthetic aperture radar (InSAR) technique has received increased attention. Differential interferometry synthetic aperture radar (D-InSAR) was proposed based on InSAR and has been developed as an available method for investigating surface deformation. As a novel earth observation technology, D-InSAR has the advantages of high precision, high spatial-temporal resolution, continuous coverage, and all-weather [18]. D-InSAR has thus become one of the most widely used measurement techniques for large spatial regions. D-InSAR has been employed to monitor subtle surface deformation resulting from underground mining activity [10,19], earthquakes [20,21], landslide activity [22][23][24], groundwater withdrawal [25,26] and so on. Moreover, the recently popular machine learning methods have been employed to identify related deformation by combining with InSAR data [27][28][29]. However, since the influence of atmospheric delay and temporal-spatial decorrelation on D-InSAR technique is significant, it has the deficiency of not being capable of extracting the long time surface deformation [10,30].
Afterwards, substantial progress in the development of advanced InSAR techniques has been made. Time-series InSAR techniques, such as persistent scatterer InSAR (PS-InSAR) [31] and small baseline subset InSAR (SBAS-InSAR) [32], have been successively proposed. At present, the SBAS-InSAR technique has been an extensively applied approach in monitoring time-series surface deformation. The SBAS-InSAR technique can reduce the influence of temporal-spatial decorrelation by recognizing highly coherent pixels [33,34]. Furthermore, it employs the least squares (LS) and singular value decomposition (SVD) algorithms to acquire surface deformation information at the high coherent pixels based on the multi-scene SAR images and thus characterize the time-series surface deformation [11,35].
For a region with abundant coal reserves, the exploitation time in different coal seams is not synchronous. Consequently, such regions may consist of multiple abandoned mined-out areas. A large number of monitoring cases related to mined-out areas have been effectively conducted using SBAS-InSAR technology. However, to our knowledge, the characteristics of surface deformation in this type of multiple abandoned mined-out areas varying in exploitation time have been rarely examined in the previous works. Nevertheless, the surface deformation problems caused by coal exploitation are complicated, given that the intensity of deformation is simultaneously affected by several factors, such as the mechanical properties of overlying rock strata, mining depth, coal seam thickness and mining method [17,36]. Furthermore, previous studies have indicated that surface deformation caused by coal exploitation tends to occur over longer timescales, lasting for several years or even decades [16]. Therefore, reasonable analyses of the spatial-temporal deformation characteristics of multiple abandoned mined-out areas based on the results of long-time monitoring are important for preventing and mitigating local disasters.
In this study, the deformation characteristics of multiple abandoned mined-out areas in two adjacent coal mines, located in Jilin Province, China, are analyzed. First, the study area and details of field investigation were described (Section 2). Subsequently, a brief description of the datasets (Sentinel-1A) and methodologies, i.e., SBAS-InSAR and ordinary Kriging method, was carried out (Section 3). The spatial-temporal evolution characteristics of the typical surface deformation were then analyzed (Section 4). Finally, the reliability of the SBAS-InSAR results was assessed via a comparison The terrain in the mining areas is relatively flat, and only the northeastern side is gradually transiting to low mountainous areas with undulating topography. The elevation ranges from 30 m to 113 m, and the landform in the mining areas is mainly alluvial plain ( Figure 1b).
As is shown in Figure 2, the sedimentary sequences in the study area include the Lower Permian Kedao Formation (P 1 k), Paleogene Hunchun Formation (E 2-3 h) and Quaternary (i.e., Q 3 and Q 4 ). The Kedao Formation is mainly composed of black-gray slate and metamorphic sandstone. The Hunchun Formation is dominated by siltstone, conglomerate, mudstone and sandstone. The Quaternary primarily consists of loose deposits exposed in the form of silty clay and gravel. The Hunchun Formation is the major coal-bearing sequence in the mining areas. It is divided into six subsections from bottom to top, namely E 2-3 h 1 , E 2-3 h 2 , E 2-3 h 3 , E 2-3 h 4 , E 2-3 h 5 and E 2-3 h 6 . The E 2-3 h 2 , with an average thickness of 145 m, is a major coal-bearing unit in the Ying'an and Baoshan coal mines ( Figure 2). It is composed of argillaceous siltstone, dark gray mudstone, medium/coarse-grained sandstone and coal seams.

Surface Deformation Features
The field engineering geological survey was carried out in the mining areas from 20 April to 8 May 2019. The investigation found that the field surface presented a series of non-uniform deformation phenomena due to the frequent mining activities.
Mining-induced collapses make the force exerted on the surface exceed its tensile strength, resulting in tensile cracks in farmlands and roads. As shown in Figure 3a,b, many tensile cracks with a nearly parallel distribution had developed in farmlands and roads. The field survey found that the cracks in the farmlands are by step-shaped with lengths ranging from 100 to 300 m and widths ranging from 0.6 to 3.0 m. The cracks are shallow in depth ranging from 0.5 to 1.0 m. These cracks had resulted in crop yield to be reduced, and in some case even destroyed cultivated land. The produced cracks directly cut off roads and pose a threat to the road surface landscape and traffic safety.
Mining-induced surface subsidence also increased the inclination of the foundation of houses, leading to cracks in the walls of the houses. As shown in Figure 3c,d, these cracks exhibited obvious radial distribution, which is a potential danger to the property and lives of local residents. (Note that the yellow circles represent the location of observation points for leveling measurements, and will be introduced in Section 5.1).
On the southwest side of the mining areas, the mining-induced subsidence produced several pits with depths between 1.0 and 2.0 m. Rainfall over long periods has gathered into the pits to form local lakes. Here, a lake is selected as representative example as shown in Figure 3f. The lake is in the shape of an approximate ellipse with a major axis of 418 m and a minor axis of 175 m. These lakes directly damage farmlands. In addition, the accumulation of rainwater in the pits over long time has an impact on the hydrogeological conditions of the area.
The surface collapse caused by the mining activities also occurred in the forest. Affected by the surface collapse, the surrounding trees tend to be bent (Figure 3e,g), which is similar to the bent trees arising from the movement of the landslide [38]. Over time, these trees may become dried-up. Therefore, mining activities can bring about some damage to the ecological environment to some extent.
In summary, mining-induced surface subsidence in the mining areas is obvious. The subsidence is a potential threat to the life and property of local residents, farmlands and ecological environment. Furthermore, filed investigations indicate that the subsidence in the mining areas is still going on. As a result, there is a need to detect the subsidence zones and quantitatively analyze the deformation characteristics of the mining area so that a valuable reference can be provided for local ecological construction and subsidence treatment.

Datasets
In order to comprehensively characterize the surface deformation of the study area, a total of 40 Sentinel-1A SAR images, collected by the European Space Agency (ESA), were utilized in this study. The time span was 14 February 2017 to 17 May 2020, with a sampling interval of approximately 30 days. The digital elevation model, i.e., Shuttle Radar Topography Mission (SRTM) is offered by the National Aeronautics and Space Administration (NASA) for removing the topographic phase component. Precise orbital ephemeris data, obtained by the ESA, were applied for orbital refinement and phase re-flattening. Table 1 shows the specific information of the datasets.

SBAS-InSAR
The SBAS is an advanced InSAR post-processing algorithm with millimeter accuracy [39], so it was employed to process the obtained Sentinel-1A SLC images to generate deformation velocity and accumulated surface deformation displacement.

Basic Theory of SBAS-InSAR
N + 1 SAR images covering the fixed area were obtained at the ordered times (t 0 , t 1 , . . . , t N ). Accordingly, the number (M) of possible differential interferometric pairs were produced based on baseline threshold values. The M should meet the following relation [40,41]: Assuming that j-interferogram is produced by combining two SAR at times t B and t A (t B > t A ). Accordingly, the interferometric phase in the pixel of azimuth and range coordinates (x, r) can be obtained as follows: where φ(t B , x, r) and φ(t A , x, r) stand for phase values of two SAR images at times t B and t A , respectively. φ disp,j (x, r) corresponds to the deformation phase, which is generated due to the change in distance between the object and the radar along the direction of the line-of-sight (LOS). In addition, phases φ topo,j (x, r), φ orb,j (x, r), φ atm,j (x, r) and φ noise,j (x, r) are caused by the terrain, satellite orbit, atmosphere and other noise (e.g., thermal noise and spatial decorrelation). To obtain deformation Remote Sens. 2020, 12, 3752 7 of 17 information, the deformation phase φ disp,j (x, r) can be obtained by removing the aforementioned phases in SBAS-InSAR. Afterwards, a system of M equations containing N unknown parameters can be expressed as ] T stands for the N × 1 vector of unknown deformation phases on the basis of measurement points, and thus δφ = [δφ(t 1 ), δφ(t 2 ), · · · δφ(t N )] T represents the vector of the unwrapped phase values related to differential interferograms. In order to acquire the deformation velocity, the Equation (3) can be rewritten as follows: Bv =δφ (4) where B stands for an M × N coefficient matrix, and v is the mean phase velocity expressed as Next, the deformation velocity can be obtained using the LS and singular value SVD. Finally, the accumulated displacement of each image relative to the first image can be calculated by multiplying the corresponding velocity and the time span [32].

Data Processing
In this study, the Environment for Visualizing Images (ENVI) and SARscape software packages were employed to obtain the surface subsidence deformation velocity and displacement time series in the study area. The steps of data processing are summarized below.

•
Step 1 Generation of the connection graph. Spatial and temporal baseline thresholds (121 m and 120 days, respectively) were first specified following the small baseline principles. The image obtained on 21 June 2018 was chosen as the super master image for interferometric combinations. The 40 SAR images were then separated into different subsets, and a total of 190 interferometric pairs were generated as shown in Figure 4.

•
Step 2 The interferometric process for generated pairs. During this process, the external SRTM was employed to eliminate the topographical phase. Furthermore, the noise phase was removed by means of the Goldstein filter, and the coherence coefficient map for each pair was obtained. Finally, the minimum cost flow (MCF) approach was utilized for phase unwrapping, and the unwrapping coherence threshold was set to 0.3 for calculation. Subsequently, 25 low-quality interferometric pairs (i.e., low coherence and poor unwrapped phase) were removed based on the above process. Finally, 165 interferometric pairs were utilized for the further time-series analysis.

•
Step 3 Generation of surface deformation. In the current study, 90 ground control points (GCPs) were selected to remove the change of phase resulting from satellite orbit and the residual topographic phase. Furthermore, the atmosphere delay was alleviated by spatial-temporal filter with the time filter window of 365 days and the space filter window of 1200 m to acquire the accumulated surface deformation value and deformation velocity along the LOS direction.

•
Step 4 Generation of surface deformation in the vertical direction. Vertical deformation is generally thought to be a dominant component in the mining-induced surface deformation. Therefore, the vertical surface subsidence (d v ) was obtained according to the deformation (d LOS ) along the LOS direction and incidence angle (θ) as presented in Equation (6) [42,43]. Finally, the deformation data were geocoded to gain the result in the form of the WGS84 coordinate system.
Remote Sens. 2020, 12, 3752 8 of 17  Step 1 Generation of the connection graph. Spatial and temporal baseline thresholds (121m and 120 days, respectively) were first specified following the small baseline principles. The image obtained on 21 June 2018 was chosen as the super master image for interferometric combinations. The 40 SAR images were then separated into different subsets, and a total of 190 interferometric pairs were generated as shown in Figure 4.  Step 2 The interferometric process for generated pairs. During this process, the external SRTM was employed to eliminate the topographical phase. Furthermore, the noise phase was removed by means of the Goldstein filter, and the coherence coefficient map for each pair was obtained. Finally, the minimum cost flow (MCF) approach was utilized for phase unwrapping,

Ordinary Kriging Model
Because of the inevitable technical problem, such as the incoherence and unwrapping problems [4], the SBAS-derived deformation results usually contain some pixels without data. However, researchers have demonstrated that the surface deformation caused by mining activities is usually spatially continuous [44]. In order to generate a continuous deformation map of the entire study area for comprehensively presenting the global deformation characteristics, the ordinary Kriging interpolation method was adopted in this study.
The ordinary Kriging model is one of the most popular approaches for spatial local prediction or spatial local interpolation in the field of geostatistics [45,46]. The fundamental assumption of the Kriging model is that closer sampled points tend to be more similar [47]. The model is based on variogram theory and structural analysis [48,49]. In order to obtain the weight coefficients, the fit of a variogram is indispensable to the spatial covariance structure of the known sampled points [50]. Therefore, the weight coefficients are not only associated with the distance between the known points and the prediction points but also with the spatial distribution among the known points [4].
In the model, the original known sampled data of the regionalized variables and variogram characteristics are utilized to handle the linear unbiased optimal estimation for the unmeasured points in the limited region. Specifically, the predicted value Z * (x 0 ) at the unmeasured point x 0 can be obtained according to the following equation [47,50]: where λ i is the Kriging weight coefficient of point x i , and n is the number of measured points within the defined region around the unmeasured point x 0 . Furthermore, the weight coefficients meet the following equation: For more detailed introduction to the ordinary Kriging model, the readers can refer to the literatures [44,51].

Deformation Velocity
Based on the SBAS-InSAR technique, the mean deformation velocity map of the mining areas was obtained during the period from 14 February 2017 to 17 May 2020 as shown in Figure 5. A total of 16,078 coherence points were extracted by SBAS-InSAR in the mining areas, which are represented by the dots in Figure 5. The velocity results were superimposed on the Google Earth image. The map was represented using a color scale from red to green. In the deformation velocity map, positive values indicate that the points are uplifting in the vertical direction during the study period, whereas negative values denote subsidence in the vertical direction.
Based on the SBAS-InSAR technique, the mean deformation velocity map of the mining areas was obtained during the period from 14 February 2017 to 17 May 2020 as shown in Figure 5. A total of 16,078 coherence points were extracted by SBAS-InSAR in the mining areas, which are represented by the dots in Figure 5. The velocity results were superimposed on the Google Earth image. The map was represented using a color scale from red to green. In the deformation velocity map, positive values indicate that the points are uplifting in the vertical direction during the study period, whereas negative values denote subsidence in the vertical direction.
The SBAS-derived results indicate that the mean deformation velocity in the study area ranges from -36 mm / y to 5 mm / y. Three major subsidence zones were observed in the study areas, which are defined as Zone A, Zone B and Zone C and are marked by the black dashed ellipses in Figure 5. These three subsidence zones are located above the abandoned mined-out areas. It can be clearly seen that all the subsidence zones are characterized by uneven deformation. The maximum subsidence deformation occurred in Zone C, with a mean deformation velocity of -35.8 mm / y. More information on the relationship between uneven subsidence and coal exploitation will be discussed in detail in Section 5. In addition, Figure 5 also illustrates that the other regions approximatively remain stable from February 2017 to May 2020 and the mean deformation velocity ranged from −6 mm / y to 5 mm / y. The mean deformation velocity of the entire study area was predicted based on the ordinary Kriging interpolation method. The exponential function is specified as the variogram in this study. Finally, the obtained mean velocity map is shown in Figure 6a. Based on the data of 16,078 coherence points, the linear regression analysis of mean deformation velocity between SBAS-derived and predicted results was carried out. The results indicate that there was an excellent correlation between them, with a correlation coefficient value of up to 0.9566 as shown in Figure 6b. Furthermore, the absolute errors of the mean deformation velocity between SBAS and the predicted results were also analyzed. The distribution map of the absolute errors is shown in Figure 6c. The absolute error of the 95% coherence points was less than 3 mm / y. The above analysis demonstrates that the interpolation results of the mean deformation velocity using the ordinary Kriging model are appropriate in the mining areas. The SBAS-derived results indicate that the mean deformation velocity in the study area ranges from −36 mm/y to 5 mm/y. Three major subsidence zones were observed in the study areas, which are defined as Zone A, Zone B and Zone C and are marked by the black dashed ellipses in Figure 5. These three subsidence zones are located above the abandoned mined-out areas. It can be clearly seen that all the subsidence zones are characterized by uneven deformation. The maximum subsidence deformation occurred in Zone C, with a mean deformation velocity of −35.8 mm/y. More information on the relationship between uneven subsidence and coal exploitation will be discussed in detail in Section 5. In addition, Figure 5 also illustrates that the other regions approximatively remain stable from February 2017 to May 2020 and the mean deformation velocity ranged from −6 mm/y to 5 mm/y.
The mean deformation velocity of the entire study area was predicted based on the ordinary Kriging interpolation method. The exponential function is specified as the variogram in this study. Finally, the obtained mean velocity map is shown in Figure 6a. Based on the data of 16,078 coherence points, the linear regression analysis of mean deformation velocity between SBAS-derived and predicted results was carried out. The results indicate that there was an excellent correlation between them, with a correlation coefficient value of up to 0.9566 as shown in Figure 6b. Furthermore, the absolute errors of the mean deformation velocity between SBAS and the predicted results were also analyzed. The distribution map of the absolute errors is shown in Figure 6c. The absolute error of the 95% coherence points was less than 3 mm/y. The above analysis demonstrates that the interpolation results of the mean deformation velocity using the ordinary Kriging model are appropriate in the mining areas.    It can be seen that the above three zones have been undergoing subsidence over time during the period from February 2017 to May 2020. The maximum subsidence values of Zones A, B and C were 70.1, 94.6 and 109.3 mm, respectively. Except for the above three zones, the accumulated deformations of the other regions in the mining areas were relatively small, with accumulated deformation values ranging from -15 to 10 mm. Three coherent points were selected along the profiles AA', BB' and CC' as shown in Figure 7, respectively, and these points are marked as A1, A2, A3, B1, B2, B3, C1, C2, C3. Then, their cumulative displacement curves over time were plotted as shown in Figure 8, showing that the cumulative subsidence deformation increases gradually.

Validation with Leveling Measurements
To confirm the reliability of the results obtained by the SBAS-InSAR technique, precise in-situ monitoring based on widely used leveling measurement was conducted [52,53]. Four-order leveling was employed in this study, whose error of height difference per kilometer is 10 mm [16]. A total of twenty observation points (marked P1～P20) were set up around the study areas. The points P1～P4 were distributed in Baoshan coal mine and P5～P20 were in Ying'an coal mine. The specific locations of observation points are presented in Figure 3a and are indicated by yellow circles. The leveling data were acquired on 5 January 2019, 10 April 019, 20 July 2019, 20 October 2019 and 15 December 2019, respectively.
In the overlapping period from January 2019 to December 2019, the accumulative deformations obtained by leveling measurements and the SBAS-InSAR technique were extracted and plotted together in Figure 9. The comparison figure intuitively indicates that the SBAS-InSAR results agree well with the results obtained by the leveling measurements. Furthermore, the root mean square error (RMSE) between the results of the SBAS-InSAR and the leveling measurements was calculated by the following equation.
In Equation (9), dlev, i and dSBAS, i stand for the deformation values of the ith observation point obtained by leveling measurements and SBAS-InSAR, respectively, and n is the number of observation points. Finally, the RMSE derived from Equation (9) was 1.144 mm, which quantitatively proves that the results obtained by SBAS-InSAR technique were reliable in this study.

Validation with Leveling Measurements
To confirm the reliability of the results obtained by the SBAS-InSAR technique, precise in-situ monitoring based on widely used leveling measurement was conducted [52,53]. Four-order leveling was employed in this study, whose error of height difference per kilometer is 10 mm [16]. A total of twenty observation points (marked P1~P20) were set up around the study areas. The points P1~P4 were distributed in Baoshan coal mine and P5~P20 were in Ying'an coal mine. The specific locations of observation points are presented in Figure 3a

The Residual Subsidence Associated with Coal Exploitation
To further explore the spatial deformation evolution characteristics associated with coal exploitation, the profiles AA', BB' and CC' shown in Figure 7 were arranged in to Zone A, Zone B and Zone C, respectively. Over the period between February 2017 and May 2020, the accumulative surface deformation values of 35 points in each profile were extracted and were presented in Figure  10. The uneven subsidence phenomenon along the three profiles could be clearly observed. The In Equation (9), d lev,i and d SBAS,i stand for the deformation values of the ith observation point obtained by leveling measurements and SBAS-InSAR, respectively, and n is the number of observation points. Finally, the RMSE derived from Equation (9) was 1.144 mm, which quantitatively proves that the results obtained by SBAS-InSAR technique were reliable in this study.

The Residual Subsidence Associated with Coal Exploitation
To further explore the spatial deformation evolution characteristics associated with coal exploitation, the profiles AA', BB' and CC' shown in Figure 7 were arranged in to Zone A, Zone B and Zone C, respectively. Over the period between February 2017 and May 2020, the accumulative surface deformation values of 35 points in each profile were extracted and were presented in Figure 10. The uneven subsidence phenomenon along the three profiles could be clearly observed. The maximal subsidence magnitude of −109.3 mm occurred on profile CC'.

The Residual Subsidence Associated with Coal Exploitation
To further explore the spatial deformation evolution characteristics associated with coal exploitation, the profiles AA', BB' and CC' shown in Figure 7 were arranged in to Zone A, Zone B and Zone C, respectively. Over the period between February 2017 and May 2020, the accumulative surface deformation values of 35 points in each profile were extracted and were presented in Figure  10. The uneven subsidence phenomenon along the three profiles could be clearly observed. The maximal subsidence magnitude of -109.3 mm occurred on profile CC'.
As shown in Figure 10, the subsidence phenomena of profiles BB' and CC' were characterized by a single subsidence bowl, whereas profile AA' was characterized by double subsidence bowls. The time interval of the three mined-out areas from formation to the present state is recorded as TA, TB and TC, respectively. The data of the mining area indicate that the order of the three time intervals was TA > TB > TC, i.e., the mined-out area A was the earliest formed. These three profiles all correspond to the locations of the abandoned mined-out area, but their deformation characteristics were not completely consistent. The above characteristics can be explained through the following mechanisms based on the SBAS-InSAR results and subsidence process.  As shown in Figure 10, the subsidence phenomena of profiles BB' and CC' were characterized by a single subsidence bowl, whereas profile AA' was characterized by double subsidence bowls. The time interval of the three mined-out areas from formation to the present state is recorded as T A , T B and T C , respectively. The data of the mining area indicate that the order of the three time intervals was T A > T B > T C , i.e., the mined-out area A was the earliest formed. These three profiles all correspond to the locations of the abandoned mined-out area, but their deformation characteristics were not completely consistent. The above characteristics can be explained through the following mechanisms based on the SBAS-InSAR results and subsidence process.
Generally, the original stress states of the overlying strata are disturbed and redistributed after the exploitation of underground coal seams, leading to the corresponding movement and deformation of the overlying strata. Meanwhile, numerous voids, including spaces between broken rocks, fractures in the overlying strata and fissures in the upper loose sediments, are produced [54]. Originally, the movement and deformation of overlying strata are mainly concentrated in the middle of the mined-out area since the middle provides the weakest support force for the overlying strata. As the overlying strata sink, the continuous deformation gradually extends upwards and eventually gives rise to the surface settlement, which was characterized by a single subsidence bowl as shown in Figure 11a. As time goes on, the regions close to the middle of the mined-out area are the first to reach a stable state. However, there is always a certain space that needs to be compacted in the boundary regions of the mined-out area as presented in Figure 11a. Consequently, later subsidence of the mined-out area mainly focuses on the boundary regions, causing the subsidence values of the boundary to be greater than those of the center. Once all types of voids are compacted, the subsidence is completed as shown in Figure 11b. The aforementioned processes can account for the surface deformation phenomenon with a characteristic of double subsidence bowls along profile AA' as marked by the black curve in Figure 10. Overall, a desirable spatial-temporal correlation was detected between surface subsidence and coal exploitation. rocks, fractures in the overlying strata and fissures in the upper loose sediments, are produced [54]. Originally, the movement and deformation of overlying strata are mainly concentrated in the middle of the mined-out area since the middle provides the weakest support force for the overlying strata. As the overlying strata sink, the continuous deformation gradually extends upwards and eventually gives rise to the surface settlement, which was characterized by a single subsidence bowl as shown in Figure 11a. As time goes on, the regions close to the middle of the mined-out area are the first to reach a stable state. However, there is always a certain space that needs to be compacted in the boundary regions of the mined-out area as presented in Figure 11a. Consequently, later subsidence of the mined-out area mainly focuses on the boundary regions, causing the subsidence values of the boundary to be greater than those of the center. Once all types of voids are compacted, the subsidence is completed as shown in Figure 11b. The aforementioned processes can account for the surface deformation phenomenon with a characteristic of double subsidence bowls along profile AA' as marked by the black curve in Figure 10. Overall, a desirable spatial-temporal correlation was detected between surface subsidence and coal exploitation.

Conclusions
In this study, a field investigation and SBAS-InSAR technology were adopted to obtain the characteristics of residual surface deformation caused by the underground mining activities in Ying'an and Baoshan coal mines. Time-series analysis of 40 Sentinel-1A SAR images between 14 February 2017 and 17 May 2020 indicated that the study area has been influenced by coal exploitation

Conclusions
In this study, a field investigation and SBAS-InSAR technology were adopted to obtain the characteristics of residual surface deformation caused by the underground mining activities in Ying'an and Baoshan coal mines. Time-series analysis of 40 Sentinel-1A SAR images between 14 February 2017 and 17 May 2020 indicated that the study area has been influenced by coal exploitation for a long time. Three dominant subsidence zones were detected in the study area, and they were all the location of abandoned mined-out areas. The maximum subsidence values of Zones A, B and C were 70.1 mm, 94.6 and 109.3 mm, respectively, during the study period. Comparison with in-situ leveling measurements demonstrated that SBAS-InSAR technology was a reliable monitoring method with high precision for characterizing mining-induced surface deformation.
According to the time-series deformation results obtained by SBAS-InSAR, mining-induced residual deformation was analyzed in the study area. Overall, the post-mining subsidence processes can be separated into two stages. In the first stage, the subsidence occurs mainly in the middle of the mined-out areas. In the second stage, the subsidence is concentrated on the boundary of the mined-out areas. Therefore, in the establishment of the subsidence prediction models of mined-out areas, the surface deformation characteristics of different stages should be fully taken into consideration in addition to the results of field investigations.