Use of Multiplatform SAR Imagery in Mining Deformation Monitoring with Dense Vegetation Coverage: A Case Study in the Fengfeng Mining Area, China

: Ground deformation related to mining activities may occur immediately or many years later, leading to a series of mine geological disasters, such as ground ﬁssures, collapses, and even mining earthquakes. Deformation monitoring has been carried out with techniques, such as multitemporal interferometric synthetic aperture radar (MTInSAR). Over the past decade, MTInSAR has been widely used in monitoring mining deformation, and it is still difﬁcult to retrieve mining deformation over dense vegetation areas. In this study, we use multiple-platform SAR images to retrieve mining deformation over dense vegetation areas. The high-quality interferograms are selected by the coherence map, and the mining deformation is retrieved by the MSBAS-InSAR technique. SAR images from TerraSAR-X, Sentinel-1A, Radarsat-2, and PALSAR-2 over the Fengfeng mining area, Heibei, China, are used to retrieve the deformation of mining activities covered with dense vegetation. The results show that the subsidence in the Fengfeng mining area reaches up to 90 cm over the period from July 2015 to April 2016. The root-mean-square error (RMSE) between the results from InSAR and leveling is 83.5 mm/yr at two mining sites, i.e., Wannian and Jiulong Mines. by mining activities during the observed period. The largest subsidence reached 88.8 cm. Joint analysis of the InSAR-derived subsidence and the ground true leveling data proved the high accuracy of the MSBAS results. The study suggests that using multiplatform SAR datasets helps to retrieve the ground deformation in mining areas with dense vegetation.


Introduction
Mining activities can break the original stress balance of the overlying rock layer, and lead to a series of geological and environmental problems, e.g., aquifer destruction, landslides, and ground structural damage [1][2][3], which inevitably lead to ground deformation. Thus, detecting and predicting such mining deformation is essential for analyzing the mechanism of mining activities and potential geological hazards. In general, mining deformation can be assessed quantitatively, by empirical and theoretical prediction models [3]. The accuracy of the predicted mining deformation of the models is limited by the input parameters related to the geological condition, geotechnical information, and deformation observed at sparse points. For example, the accuracy of the predicted deformation could range from 48% to 773%, without actual ground deformation information [4]. Therefore, monitoring ground deformation plays an essential role in exploring the deformation mechanism of the mining area.
Ground deformation can be measured by in situ observations from traditional geodetic surveying techniques, such as leveling and global navigation satellite system (GNSS).
These techniques can monitor ground deformation with high accuracy; however, the spatial resolutions are normally poor, and they are often labor intensive. Mining prediction models highly rely on the resolution, sampling rate, and coverage of measurement points. Spaceborne interferometric synthetic aperture radar (InSAR) is an active remote sensing technique, and has become one of the most useful geodetic techniques for ground deformation monitoring [5][6][7][8][9]. Compared with traditional techniques, InSAR has proven its unique advantages through broad coverage and high spatial resolution [10,11]. To further improve the accuracy of traditional InSAR, multitemporal InSAR (MTInSAR) has been developed for analyzing the time series phase of point targets by a series of interferograms. For example, persistent scatter InSAR (PSInSAR) and small baseline subset (SBAS) methods were proposed to reduce the effect of unwanted phase contributions, such as atmosphere, decorrelation noise, and digital elevation model (DEM) phase residual, and to retrieve more accurate ground deformation [12,13]. Currently, MTInSAR has been widely applied in monitoring the time series of deformation in different scenarios, such as seismic motion, landslides, and ground subsidence [14][15][16][17][18][19][20][21][22]. In addition to mining activities, MTInSAR offers great potential for promoting studies of mining deformation, such as mechanism interpretation, parameter inversion, and deformation prediction [3].
In 1996, Carnec et al. [23] generated the first deformation map caused by underground mining, using InSAR. Then, many applications have been carried out to monitor mining deformation with InSAR in recent decades. For example, in 1999, Wright et al. [4] obtained elevation changes induced by mining activity with two ERS SAR datasets. In 2010, Milan et al. [24] retrieved a series of mining deformations in the area of the Ostrava-Karviná Revier in the Czech Republic using the Stanford method for persistent scatterers (StaMPS). Since mining deformation often occurs in three-dimensional (3D) space, studies with individual SAR datasets using the SAR techniques, such as multi-aperture InSAR [25], and pixel offset tracking [26], can only be dedicated to measuring one-/two-dimensional (1D/2D) deformation along the radar line of sight (LOS) and/or azimuth directions, by which the real pattern of mining deformation is difficult to elucidate. To solve this problem, several methods have been proposed to reconstruct 3D mining displacements, by combining InSAR observations from multitrack or multiplatform systems. For example, in 2011, Hay-Man Ng et al. [27] retrieved the 3D deformation of the southern highland coalfield in New South Wales, Australia, with multiplatform SAR datasets. In 2013, Samsonov et al. [2] proposed the MSBAS method to retrieve the 2D time series of ground deformation of postmining activities along the French-German border.
MTInSAR has been widely used to monitor mining deformation, but it is still limited by environmental issues, such as dense vegetation and extreme weather. Mining areas are often located in rural regions, which are usually covered with dense vegetation, which raises the significant decorrelation of InSAR interferograms and leads to spatial and temporal gaps for mining deformation monitoring [28]. In this article, we carried out a case study of deformation monitoring in a mining area with dense vegetation coverage. We use multiplatform SAR images and MSBAS-InSAR technology to retrieve the mining deformation. The dataset from different SAR platforms is integrated into the MTSBAS-InSAR to estimate the mining deformation. High-quality interferograms are selected according to their coherence maps and weather conditions. The generalized deformation maps of 36 mining sites in the Fengfeng mining area, Hebei, China are obtained by combining the SAR images of TerraSAR-X, Sentinel-1A, Radarsat-2, and PALSAR-2. The deformation results are verified with the measurements from leveling and indicate that the high accuracy of the deformation can be achieved with multiplatform SAR images.
The remainder of this paper is structured as follows. In Section 2, the background of the Fengfeng mining area and SAR datasets used are introduced. Section 3 gives the strategy of data processing and the MSBAS-InSAR method. The result and discussion are shown in Sections 4 and 5 respectively. Section 6 gives the conclusions of this study.

Background of Fengfeng Mining Area
The Fengfeng mining area is located in the southwestern part of Handan city, China, as shown in Figure 1. Approximately 36 sites and significant subsidence (see Figure 1b) have been reported in and around this area [29,30]. According to the field work photo of Fengfeng mining area in Figure 1b, the ground surface is covered with dense vegetation, e.g., crops and weeds, leading the InSAR interferogram to be easily decorrelated. Especially in the summer period, vegetation grows quickly, and some vegetation is over the height of a man. In addition, the ongoing mining activities has generated large gradient subsidence that cause difficulties to InSAR for monitoring of deformation.

SAR Datasets
Multisensor spaceborne SAR satellite datasets from both ascending and descending tracks over the Fengfeng area were collected to measure the mining deformation, as outlined in Table 1. The coverage of these data and the distribution of the main mining regions are shown in Figure 2. A digital elevation model (DEM) from shuttle radar topography missions (SRTM), provided by the United States Geological Survey (USGS) with a spatial resolution of 30 m, is adopted for topography compensation and geocoding. We also collected measurements of 152 leveling stations in the 6 phases at two mining sites, i.e., Jiulong and Wannian Mines, to validate the InSAR results.

Interferometric Processing
The data preprocessing of each SAR dataset used in this study was carried out using GAMMA software [31], including single look complex (SLC) image extraction and coregistration, interferogram generation, flat-earth phase removal, phase noise filtering, phase unwrapping, and geocoding.

Retrieval of Time-Series Deformation
In this study, the MSBAS technique in [2] was used to retrieve the time-series deformation of the Fengfeng mining area. MSBAS is an advanced extension of the traditional SBAS [13], that extends the deformation retrieval from the line-of-sight (LOS) direction to the two-dimensional direction, i.e., the horizontal (east-west direction) and vertical (updown direction) deformation. In addition to the advantage of 2D deformation, MSBAS can significantly improve the temporal sampling of deformation time series with multiplatform SAR images. The difference from the SBAS method is that before time-series analyzing, all unwrapped interferograms must be geocoded and interpolated into a common grid with the same spatial resolution, as shown in Figure 3. The lower resolution of the SAR dataset is usually chosen as the resolution of the common grid, which is the final spatial resolution of the deformation result. Both SBAS and MSBAS are briefly introduced in the following sections.

SBAS Technique
For each SAR data processing, it is assumed that N SAR images from a single platform are obtained at the time order of t 1 , t 2 , t 3 , · · · , t N . M differential interferograms (where ) with a short multimaster baseline are generated and well phase unwrapped. The unwrapped phase (∆ϕ) of one pixel on an interferogram can be written as: where ∆ϕ de f o is the deformation phase in the LOS direction between the period of two SAR acquisitions; ∆ϕ DEMerror is the topographic residual phase, caused by the imperfection of the external DEM; ∆ϕ orb is the phase related to the orbital error that can be compensated by the orbital model and external information [32,33]; ∆ϕ atm and ∆ϕ dec are the atmospheric phase screen (APS) and decorrelation noise phase, respectively; and λ, D LOS , θ, B ⊥ , and R are the radar wavelength, deformation, incident angle, perpendicular baseline, and slant range distance between the radar sensor and ground target, respectively.
After compensating for the orbital error (∆ϕ orb ) and part of the atmospheric phase (i.e., altitude-related), Equation (1) can be expressed in the form of matrices when all interferograms are taken into consideration: where ∆ϕ is the unwrapped phase vector of a pixel in all interferograms with the dimension of M × 1, X is the P × 1 unknown parameter vector that is used to model the deformation phase (D LOS ) and DEM error (H) [13], e is the noise phase component mainly including the atmospheric phase and decorrelation noise, A is a design matrix linking the unknows and unwrapped phase, which has M rows and P columns. Generally, the solution can be estimated using least-square inversion, i.e., Regularization methods, such as Tikhonov, can also be used to solve the problem when design matrix A is rank deficient [34].
The LOS deformation from an individual SAR geometry or platform can be expressed as [35], where β is the azimuth angle of the flight direction. D U , D N , and D E are the deformation vectors along the up-down, north-south, and east-west directions, respectively. Equation (3) can be rewritten in the form of a matrix as: where S is the LOS unit vector with north, east, and up components, i.e., S N = sin(β)sin(θ), S E = −cos(β)sin(θ), and S U = cos(θ). For individual SAR geometry, the result in Equations (3) or (4) cannot be decomposed into 3D displacements. The SBAS methodology can be modified to generate an approximate solution of time series displacement in two or three dimensions, when SAR images are obtained from more than one acquisition geometry [36].

MSBAS Technique
In MSBAS method, each SAR dataset should have an overlap period in time and with different observation geometries. After the data processing with SBAS method, the DEM errors, orbital errors, and atmospheric errors related to altitude are compensated from the unwrapped interferograms for each SAR dataset to achieve the deformation interferograms. Then, the unwrapped deformation interferograms of each SAR dataset are interpolated into the common grid, as shown in Figure 3. The 3D displacement time series on the common grid point can then be resolved as: where ∆ϕ is the unwrapped phase vector generated by the multiplatform SAR images, A is the coefficient matrix, and depending on the interferograms generated, S is the LOS unit vector containing geometric information. Currently, most spaceborne SAR satellites have near-polar sun-synchronous orbits, the SAR sensor is relatively insensitive to displacement in the north-south direction, compared to the vertical and east-west directions. The deformation in the north-south direction can be often ignored [36], and Equation (5) can be then rewritten as: where the coefficient matrix is updated toÂ = {S E A S U A}, ζ is the regularization factor, and I is the unit matrix. Similar to the SBAS method, equation system (6) may be a rank-deficient problem. Transactive singular value decomposition (TSVD) and Tikhonov regularization can be used to find the least square solution in the sense of the minimum norm. In MSBAS, Tikhonov regularization is applied. However, due to the effect of atmospheric phase and decorrelation noise, the solution of Equation (6) can be biased when the time period between each SAR dataset does not overlap or match very well [36]. Therefore, the displacement can be modelled by the deformation model with unknown parameters. For example, an interval deformation velocity model is used to model the 3D displacement, Equation (5) can be rewritten as where V E and V U are the interval deformation velocities in consecutive SAR acquisitions along the east-west and vertical directions, respectively.Â is the coefficient matrix, and contains the time intervals between consecutive SAR acquisition. The solution of Equation (7) is solved by applying Tikhonov regularization of a second order. In the regularization parameter, ζ = 0.1 is used to produce a smooth solution. The cumulative displacement of the 2D time series can be recovered by the deformation model according to time. For the interval linear deformation velocity model, the 2D displacement time series can be recovered by Equation (8).
The accuracy of MSBAS is highly dependent on the quality of the observation, i.e., ∆ϕ. The phase quality of an interferogram is easily affected by decorrelation and atmospheric artifacts. The selection of interferograms is therefore essential in MSBAS data processing. For this reason, all interferograms in which the spatiotemporal baselines were within the given thresholds (see the detailed threshold value in Table A1 in Appendix A) were generated from the available SAR images in this case study. A multilook operation was carried out to suppress the phase noise. Due to the dense vegetation cover during the observation period, we removed interferograms with low coherence in the mining area (mean coherence smaller than 0.28). Phase unwrapping was performed using the MCF method in GAMMA software. In addition, since the decorrelation effect caused by dense vegetation heavily affected the usage of SAR images with different wavelengths, in this case study, we attempted to gradually integrate the multiplatform SAR interferograms with high phase quality into the MSBAS estimator, to reduce the effect of decorrelation for retrieval of mining deformation. Detailed results are presented in Section 4.

Deformation with Single SAR Platform
The mining area suffering ground deformation was identified first, before InSAR time-series data processing. We used two TerraSAR-X images with a 22-day revisit period to generate the deformation interferogram. The results in Figure 4 show that many de-formation areas probably induced by mining activities were identified. Five segmentation areas (i.e., A, B, C, D, and E, as shown in Figure 4 with white rectangles) were divided to clarify the mining regions, which was used to identify the phase quality of the interferograms in the mining area. After the identification of the mining site, we used 14 TerraSAR-X images from June 2015 to April 2016 to estimate the ground deformation over the Fengfeng mining area. Thirteen interferograms with short baseline information were generated, as shown in Figure 5. The coherence maps of each segmentation area were used to select the highquality interferograms for the MSBAS method, and some examples of coherence maps are given in Figure A1 of Appendix A. According to the coherence maps, the interferogram with 11-day and 22-day time intervals was more effective in monitoring mining deformation with dense vegetation. In addition, due to the heavy weather conditions shown in Figure A2, the surface crops may have changed rapidly during the observation period, leading the backscattering coefficient of the object to vary greatly, and some interferograms with short spatiotemporal baselines may still be affected. If interferograms with longer temporal baselines were included, we could not obtain many coherent pints in the central area of mines, due to the decorrelation caused by dense vegetation. Therefore, we have made corresponding adjustments to only use the data from December 2015 to April 2016 with mainly 11-day baseline interference pairs, to eliminate severe decorrelation in large areas. Seven interferograms with good phase quality were generated according to the spatiotemporal baseline in Figure 5 with a red network. Since only the TerraSAR-X SAR dataset was used in this section, the SBAS method was used to retrieve the mining deformation. Figure 6 displays the LOS deformation velocity map of the investigated area. The whole area remained stable over the observed period, except for the significant mining deformation over the mine area. The mining deformation of five segmentation regions was well retrieved. Denser mining areas were identified separately, and the deformation boundary of two closer mining areas were also identified, indicating the applicability of the TerraSAR-X in mining deformation areas during the winter period. The maximum deformation among the mining sites reached approximately 100 cm/yr. Two mining sites, Wannian and Jiulong, are displayed in the subfigure of Figure 6. The maximum LOS deformations of these two sites were approximately 23 cm and 16 cm, respectively. The time-series deformation of the Fengfeng mine is shown in Figure 7. During the observed period, all the mining sites suffered significant deformation, and the deformation of each mining site started at different moments, which indicates that mining activities were carried out in all mining areas during the observation period, and the scale of mining activities and working hours varied. The maximum LOS deformation related to mining activities was approximately 30 cm.

Deformation with Two SAR Platforms
To retrieve the 2D displacement of the Fengfeng mine, the Sentinel-1A dataset from the ascending orbital direction between December 2015 and April 2016 was used for estimation, together with the TerraSAR-X dataset. The spatiotemporal baseline configuration of the two SAR datasets is shown in Figure 8a, in which the temporal baseline of both datasets was well matched. Thus, the MSBAS technique in Section 3.2 was used to retrieve the 2D mining displacement. Since the coverage of Sentinel-1A is larger than the coverage of TerraSAR-X, the common area of the displacement has the same coverage as the TerraSAR-X dataset. The retrieved 2D displacement of the Fengfeng mine is shown in Figure 9. The maximum accumulated subsidence in Figure 9a reached 33 cm during the observed period. Due to the surface subsidence of the mining area, the edges of the subsidence area collapsed toward the center area, resulting in a horizontal east-west displacement that moved to the center area. Therefore, each mining site suffered from not only subsidence but also horizontal displacement. Figure 9b shows that the relative horizontal displacement of the Fengfeng mining area reached a maximum of 20 cm. Since InSAR is insensitive to horizontal deformation in the north-south direction, the horizontal displacement corresponding to the north-south direction was difficult to retrieve.

Deformation with Three SAR Platforms
Since Sentinel-1A is susceptible to the influence of the atmosphere, we added a set of Radarsat-2 images with the C-band to detect the corresponding 2D displacement results. The baseline configuration is displayed in Figure 8b. The coverage of Radarsat-2 could not fully cover TerraSAR-X data, as shown in Figure 2, and the 2D displacement of some research areas could not be retrieved. The final 2D displacement is shown in Figure 10. The maximum accumulated subsidence in Figure 10a reached 33 cm during the observed period. The vertical subsidence of the Wannian mining site was enlarged, and the horizontal displacement of the B mining region was more significant.

Deformation with All SAR Platforms
PALSAR-2 data use a long radar wavelength in L-band, and can penetrate part of the vegetation. We therefore used three available PALSAR-2 data to retrieve the mining deformation in the summer period, and extended the deformation time series for MSBAS. The ground deformation in mining areas often produced regional subsidence that funnel along its dig working face, as sketched in Figure 11. Largest vertical subsidence suffered in the regions near the center axis of the dig working face, while horizontal displacements of these regions was close to zero, as the results of the Wannian mining area show in Figures 9 and 10. For ascending and descending orbital geometry, the horizontal deformation of the area near the mining center caused completely opposite effects, as shown by points A and C in the Figure 11. With multiplatform SAR images from opposite observed geometries in a common period, the influence on the combined deformation results will be weakened. Therefore, we ignored the horizontal displacement and converted all the LOS direction deformations into vertical, to solve the longer time series of mining subsidence. The subsidence of the Fengfeng mine from July 2015 to April 2016 is shown in Figure 12a. Due to the difference between the coverage of each SAR dataset, the deformation map was reshaped again. The subsidence of most of the Fengfeng mining regions was retrieved, and the maximum subsidence was approximately 90 cm. The subsidence map larger than 5 cm is displayed in Figure 12b. In total, 36 deformed mining sites were revealed, and their time-series subsidence is shown in Figures A3 and A4 of the Appendix A. The maximum subsidence of each mining site is listed in Table A2 of the Appendix A. The maximum subsidence occurred in the D2 mining region and reached 88.8 cm. The ground coverage of the subsidence with different intervals is listed in Table A2. Compared with the official report investigation in 2016 [37], the deformation areas from the InSAR results were larger than the report provided by the government, indicating that the subsidence area has increased, due to underground mining activities. It also demonstrated the usability of the InSAR method in the mining area with denser vegetation coverage.

Accuracy Validation
We used the subsidence results of ground leveling measurements at two mining sites, Wanjian and Jiulong Mines, to validate the results derived from the multiplatform SAR datasets. The location of the leveling points is shown with the white dots on in Figure 13. The path of the leveling stations was along the side part of the mining area, to obtain ground subsidence. We compared the time-series subsidence between the MSBAS-derived results and the leveling on twelve selected points, which were distributed uniformly over the mining area, as labeled in Figure 13. Because the resolution of the MSBAS-derived results depends on the common reference grid of the multiplatform SAR dataset, we selected MSBAS measurements with distances to leveling stations smaller than 100 m, to compare with the leveling data of the two sites. All qualified MSBAS measurements were used to calculate the uncertainty and are displayed in Figures 14 and 15. The results of MSBAS agreed well with the leveling measurements. The differences between them might relate to the uncertainty of the atmospheric delays, decorrelation effect, and differences in location and time between the InSAR and leveling data, such as the unmatched data in the last period of the two leveling stations at the Jiulong site, i.e., JL03 and JLQ39.   Since the date of the acquisition of leveling and the SAR images were quite different, it was difficult to compare the exact value between them, thus, the deformation velocity was used to make the cross validation. We therefore calculated the subsidence velocity with the leveling data. The points of MSBAS measurements with a distance smaller than 100 m were compared with the leveling data, to reduce the effect of the shift in geolocation. Figure 16 gives the scatterplots between the leveling and the results of InSAR at the two mining sites. The corresponding RMSE is 83.51 mm/year. Due to the limited number of interferograms used in retrieving the deformation, the RMSE between InSAR and leveling results is a little bit large. This value could be further reduced if more high-quality interferograms with short temporal baseline can be obtained from multiplatform SAR, however, it is difficult for the monitoring scenario like the region in this study with the state-of-the-art SAR missions. The significant difference of leveling stations that is framed with a dashed black rectangle, may be caused by the large gradient deformation in the center of the mining site that is leading to the unwrapping error [3]. According to the presented studies in [4,38], it has been suggested that the accuracy of mining subsidence measured by InSAR are from millimeters to few centimeters. Therefore, the accuracy of mining deformation detected by the MSBAS and multiplatform SAR datasets is acceptable for this study.

Advantages of Multiplatform SAR Data in Mining Deformation Monitoring
Because mining areas are often located in rural areas covered with dense vegetation, ground coverage, extreme weather, and large gradients of deformation will raise the decorrelation of InSAR interferograms and lead to incomplete spatial coverage with information gaps for mining deformation monitoring. It is worth noting from this study that the L-band PALSAR-2 data were effective for mining deformation monitoring, even in summer, and pairs with long baselines can be included in the interferometric analysis; therefore, it would be good to use the L-band data for mining deformation monitoring. However, the revisit period (i.e., from 40 to 90 days) of PALSAR-2 is too long to achieve enough temporal samples of the deformation for further analysis of the mechanism of mining activities. The objective of our contribution was to show that important improvements to these limitations are still possible. In denser vegetation areas, the mining deformation can be retrieved by properly selecting the high-quality interferograms from multiplatform SAR datasets. The fusion processing of multiplatform data can not only obtain the vertical subsidence of the mining area but also obtain the horizontal displacement. Meanwhile, a more satisfactory coverage of spatial monitoring and the increased temporal sampling of the deformation can be obtained, compared with the single SAR platform dataset.

Conclusions
Deformation associated with mining activities is usually difficult to precisely monitor in regard to both time and space because of the diverse topography, large-gradient ground subsidence, and dense vegetation coverage. These problems possibly result in decorrelation in InSAR, and therefore limit the application of InSAR in mining areas. In this study, we used multiplatform SAR data and MSBAS technique to monitor mining deformation over the Fengfeng mining area, China, from July 2015 to April 2016. The derived results indicated that the Fengfeng mining area experienced serious ground subsidence caused by mining activities during the observed period. The largest subsidence reached 88.8 cm. Joint analysis of the InSAR-derived subsidence and the ground true leveling data proved the high accuracy of the MSBAS results. The study suggests that using multiplatform SAR datasets helps to retrieve the ground deformation in mining areas with dense vegetation.

Data Availability Statement: Not applicable.
Acknowledgments: The authors would like to thank the DLR for providing the TerraSAR-X dataset under the science proposal LAN2979 and GEO3603, and the ESA/Copernicus for providing the Sentinel-1A SAR images, respectively. The Radarsat-2 was from CSA. The PALSAR-2 was provided by the JAXA under the science proposal PI3380. The authors would like to thank the China University of Mining and Technology for providing leveling data. The authors would also like to thank the anonymous reviewers for their constructive comments and suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
This appendix has shown the details and data supplemental to the main text.