Monitoring nonlinear and fast deformation caused by underground mining exploitation using multi-temporal Sentinel-1 radar interferometry and corner reflectors: application, validation and processing obstacles

ABSTRACT In this study, Differential Interferometric Synthetic Aperture Radar Interferometry (DInSAR) of artificial Corner Reflectors (CRs) were validated in the area of fast and nonlinear deformation gradient caused by active coal longwall exploitation. Three Sentinel-1 datasets were processed using conventional DInSAR, Persistent Scatterer Interferometry (PSI), and Small BAseline Subset methods implemented in ENVI SARscape™. For evaluation, leveling and Global Navigation Satellite System (GNSS) measurements were used. Considering the challenge of snow cover, the removal of all winter images was not a successful strategy due to the long temporal baseline and strong movement, which cause phase unwrapping problems and underestimate the real deformation. The results indicate that only conventional DInSAR and SBAS with low network redundancy allow us to capture maximal deformation gradient and the root mean square error calculated between the CRs and the ground truth is on the level of 2–3 cm for the vertical and easting deformation component, respectively. For the small deformation gradient represented by the permanent GNSS station (4 cm/year), all SBAS techniques appeared to be more accurate than DInSAR, which corresponds to higher redundancy and better removal of the atmospheric signal. In contrast, DInSAR results allowed to capture information about two subsidence basins, which was not possible with SBAS and PSI approaches.

In literature, various InSAR stacking-based techniques are used to generate interferometric products including (i) Persistent Scatterer Interferometry (PSI) proposed by e.g.Ferretti, Prati, andRocca 2000, 2001;Hooper et al. 2004;Devanthéry et al. 2014 or by Kampes and Hanssen 2004;(ii) Small BAseline Subset (SBAS) (Berardino et al. 2002) and (iii) mixed methods proposed by Hooper 2008or Ferretti et al. 2011 (SqueeSAR).Each of these methods has certain advantages and disadvantages in the context of accuracy, measurement density, and suitability to monitor vegetated areas.Stacking-based InSAR methods take advantage of large stacks of SAR images (minimum 20) to effectively estimate and eliminate atmospheric or topographic artifacts, and hence to increase the accuracy of the final deformation estimates (Ferretti, Prati, and Rocca 2000).This stack of SAR images is processed using Differential Synthetic Aperture Radar Interferometry (DInSAR) and then various assumptions connected with atmosphere behavior in time and space as well as deformation behavior (model) in time (Ferretti, Prati, and Rocca 2000;Van Leijen and Hanssen 2007) or space (Hooper et al. 2004) are utilized to properly estimate each interferometric components.When the a priori information about the temporal behavior is not known, a linear deformation model in time is usually applied.When the outliers from the linear deformation model are not significant (mm level) then also nonlinear deformation component in time can be estimated as presented in Ferretti, Prati, and Rocca 2000.Moreover, the quality assessment of Persistent Scatterers (PS) is dependent on the deformation model chosen.Persistent Scatterers (PSs) may be falsely rejected due to model imperfections and therefore, not presented in the final results even when many natural scatterers are located in that areas.To accept these PSs, more advanced deformation models should be used (Van Leijen and Hanssen 2007).
Most underground mining activities are carried out in various dynamics depending on the plans of the excavation, mining depth, geological conditions, etc.Therefore, deformation caused by mining activity have also different behavior in time.Generally, the mining-affected deformations around the center of the underground work face are very fast and nonlinear and thus may reach even a couple of meters per year.This causes phase unwrapping problems in DInSAR processing especially for shorter wavelength (X-band) and therefore underestimation ofe the real deformation as presented by Przyłucka et al. 2015 or Wang, Deng, andZheng 2020a.As an example, Figure 1 presents the visualization of InSAR results from EGMS for the deformation basins caused by the active Rydułtowy mine in Poland.Polygons with blue color were marked to point out areas where many natural scatterers exist which are perfect PS candidates.However, interferometric measurements are not available in such regions due to the model imperfection used in PSI processing.Therefore, for that areas, PSI processing with different deformation models\assumptions should be used.However, for EGMS this is understandable that various phenomena represent various deformation behavior in time, and therefore for monitoring deformation caused by various phenomena on a continental scale, one deformation model in time is not enough.For this specific case study with highly nonlinear and fast deformation gradient, the maximal deformation rate detected for one PS in EGMS was on the level of 4 cm/year while authors of Pawluszek-Filipiak andBorkowski (2020a,2020b) presented the deformation rate of that area of approx.0.5 m/year in the center of the subsidence bowls (blue polygons in Figure 1).The most deforming areas with the highest outliers from the linear model are missing PS.Unfortunately, these areas are of the most interest from the mining perspective and associated mining-related hazards.
Therefore, the objective of this study is to evaluate various interferometric techniques and its parameters together with Sentinel-1 data in the context of monitoring fast and nonlinear deformations represented by mining activity to fulfill that missing information in EGSM service.To achieve more measurement points and omit coherence-related problems, five corner reflectors (CRs) have been placed in the area of the study.The application of CRs has been reported in various studies such as landslides (Fu et al. 2010;Schlögel et al. 2017), open pit mines (Gama et al. 2017), and underground mining (Xing et al. 2013).However, so far there is a limited number of studies that demonstrates the application of CR-InSAR (Corner Reflector InSAR) for fast and nonlinear deformation caused by underground mining activity.These conditions significantly limit the successful CRs application and bring many processing obstacles which need to be considered during CR-InSAR processing to get reliable and accurate results.Most of the CR-InSAR studies demonstrate CRs applicability for cm-level deformation (e.g.Xing et al. 2013) usually with linear deformation behavior time and rate of 2-3 cm per year (Gama et al. 2017).Therefore, to overcome this research gap, we present the application of CR-InSAR in the area of fast nonlinear movement reaching approximatelly 50 cm per year with additional challenges represented by snow coverage and vegetation grownth.Each CR was designed and placed in areas with various deformation rate and one of the CR are placed almost in the center of the subsidence basin reaching almost the maximum deformation rate (43 cm/year) represented by one subsidence basin.Estimation of the deformation components for the CRs was monitored using conventional DIn-SAR, PSI (Ferretti, Prati, and Rocca 2000), and SBAS (Berardino et al. 2002) techniques implemented in ENVI SARscape™ software with various connections of interferometric pairs and redundancy.Two leveling and GNSS measurement campaigns of CRs and also high-frequency GNSS allows to evaluate various processing obstacles and associated errors.For the two most accurate CR-InSAR results (DInSAR, low-redundancy SBAS), we performed a 2D decomposition to show the final deformation behavior in the active mining exploitation in Poland.

Characteristics of the study area
The Rydułtowy mine is the oldest active mining plant in the Upper Silesian Coal Basin (USCB) in Poland, one of the largest hard coal mining areas in Europe.The total area of USCB is about 7400 km 2 (Popiolek et al. 2001).The largest part of this mine region is located in Poland, and the remainder lies in the Czech Republic.It is located in the southwest part of the USCB, and covers around 46 km 2 .The location of the study area is illustrated in Figure 2. The first excavation has been dated to 1792, and the average daily production of the mine ranges between 9000 and 9500 t/day.Coal mines in Poland are very deep mines with an average exploitation depth of 800 m.In a presented mine, a depth of extraction is between 800 and 1200 m (Pawluszek-Filipiak and Borkowski 2020a).The residue from the excavation is taken to the Szarlota slag heap, which is formed from post-mining waste.It is 134 m in height and is the largest in Europe, with an area of 37 hectares and a volume of 13.3 million m 3 (Pilecka and Szermer-Zaucha 2015).According to the Polish Mining Group, coal in the form of deposits makes up around 80.4 million tonnes, and excavation is expected until around 2040.
In addition to significant mining deformation, the intensive extraction of raw materials and the complicated geological structure cause numerous high-energy tremors in the area of interest.Seismic events have typical magnitudes of between 1.0 and 3.5 M and are located at a depth of around 290-1000 m (Pilecka and Szermer-Zaucha 2015).The effects of tremors can be observed at and below the ground surface.Rock bursts lead to damage to the underground infrastructure of the mine, and sometimes to the deaths of miners.The tremors can also be felt by the inhabitants of nearby towns (Dubiński, Stec, and Bukowska 2019).Therefore, deformation monitoring is important in this area to mitigate mining hazards and can be facilitated by InSAR.The area of Rydułtowy mine is mostly covered by urban and rural settlements.In addition to Rydułtowy, there are other larger towns in the region, such as Radlin or Pszów.

Auxiliary data
In the area of interest, various geodetic measurements have been carried out within EPOS-PL and its continuation of EPOS-PL + projects.The EPOS-PL is the Polish realization of the European Plate Observing System (EPOS) initiative, the main goal of which is to integrate existing and newly created research infrastructures to facilitate the use of multidisciplinary data and products in the field of earth sciences in Europe.
Five CRs were placed in the area of the Rydułtowy mine, as shown by the pink symbols in Figure 2(c).Field photographs of some of the CRs are shown in Figure 3.These CRs were located in the area of various deformation gradients.A suitable orientation for the CRs allows the provision of stable backscattering in ascending and descending Sentinel-1 geometries.Four of them (CR01-CR04) were located near the Szarlota slag heap, while CR06 was located around 1 km to the southeast of Szarlota.All CRs were mounted on stable pillars at a height of 1.5 m.
One high-frequency GNSS receiver in the area of interest was also purchased as part of the EPOS-PL project (shown in the pink square in Figure 2), which enabled deformation monitoring in near real-time (NRT) and, in contrary to InSAR, in three dimensions (3D).

Data
In this study, we used ascending and descending Sentinel-1A/B TOPSAR images covering the study area, with a six-day revisiting period.The images were acquired between 25th June 2018 and 14th July 2019.Two conventional geodetic measurement campaigns were also carried out to precisely estimate the displacements of the CRs, which were measured using geometric leveling and static GNSS measurements.The first measurement campaign was carried out between the 2nd and 4th of July 2018, and the second between the 28th and 30th of June 2019.GNSS static measurements were performed over three independent sessions.The coordinates of the high-frequency receiver RES1 were also estimated at the postprocessing stage and used for validation; however, RES1 began operation on the 10th of August 2018, and data from this receiver were therefore only available after this date.Additionally, to investigate the quality of the InSAR measurements, weather information was captured from the Institute of Meteorology and Water Management (IMGW), which was acquired by the meteorological station located in Rybnik.Data on snow coverage  were taken from a weather information service provided by www.ogimet.com.Table 1 presents basic information about the data used in this study.

Methodology
An overview of the full methodology applied in this study is presented in Figure 4. We estimated the terrain deformation over the Rydułtowy mine using ENVI SARscape™ interferometric processing software (Sarmap 2009(Sarmap , 2012) ) which offers various stacking-based DInSAR methods including: (i) classical multi-temporal consecutive DInSAR (ii) PSI algorithm (Ferretti, Prati, and Rocca 2000) and (iii) SBAS approach (Berardino et al. 2002).All approaches are described in more details in Section 4.2.PSI and SBAS methods were applied with various interferometric pairs connections.
For the most accurate InSAR solutions, Line of Sight (LOS) deformations was decomposed into vertical and east-west deformation components by using Least Squeres Estimation (LSM) and compared with ground truth data of CRs (GNSS and leveling) and GNSS permanent observations at the RES1 station.Postprocessing, analysis, and visualization were carried out in Python and ArcGIS.

InSAR pre-processing: evaluation of signal quality on CRs
The first step in our research was to evaluate the signal quality of various interferometric pairs on CRs.The main disadvantage of passive CRs is that the stability of the radar backscattering is destroyed if any object or material accumulates on the CR's surface.The main problems are snow coverage during the winter season but also the accumulation of other materials (rubbish, soil, leaves, etc.)For specific periods (January and February), low signal stability and therefore low coherence were observed due to the snow coverage.Moreover, also low coherence was observed for some single interferograms.Since the interferometric pairs in the DInSAR, PSI, and SBAS methods are connected in different ways, coherence analysis was calculated and investigated separately for these methods.We integrated the coherence time series with meteorological data, such as precipitation and snow coverage, and observed a significant loss of coherence within certain specified periods and for particular interferometric pairs.One strategy is to remove all winter images to preserve similar backscattering conditions, however, during the winter season, not all winter pairs are affected by snow coverage in our climate.Therefore, we decided in further processing to test two processing workflows for PSI and SBAS: (i) strategy when only a few images which represent small coherence will be removed (ii) strategy when all winter images from January-February will be removed.In each interferometric processing, new interferograms were then recalculated based on the subsequent image with a larger temporal baseline (e.g.12-18 or 60 days).In DInSAR processing, also low coherence pairs were removed, however, due to the direct integration of subsequent interferograms, the strategy with the removal of all winter images was not tested because of the huge temporal baseline.Additionally, by taking advantage of SBAS method to set various network redundancy represented by maximal temporal baseline, we tested two redundancy strategies to estimate if some phase unwrapping problems associated with the interferograms with high temporal baseline and therefore big deformation gradient are not influenced by the phase unwrapping problem the final displacement estimates.

InSAR processing
As mentioned above, displacement time series were estimated using DInSAR, PSI, and SBAS with various strategies.As an example, some strategies and their interferometric connections are presented in Figure 5.

Multi-temporal consecutive DInSAR processing
In this processing workflow, interferograms were generated from consecutive pairs of SAR images (i.e.φ 1-2 , φ 2-3 , φ 3-4 , … φ n-1,n ) as presented in Figure 5(a).The results were combined to create a complete time series.For the Sentinel-1 A/B data, most of the pairs had a time span of six days.The use of small temporal baselines minimizes the temporal decorrelation, and this is a significant advantage of this approach.However, this processing method has also certain limitations.When errors arise in any of the calculated deformations (e.g.due to residual terrain, atmospheric delay, or other phase errors) within the interferograms, these will affect the subsequent time-series deformation results (Pawluszek-Filipiak and Borkowski 2020b; Yang et al. 2017).
Figure 5. Examples of various interferometric pairs used for CR-InSAR processing with winter period marked with light blue color (a) DInSAR-direct integration of SAR images with some images removed, (b) PSI-some images were removed, (c) PSI-all images from Jan-Feb were removed, (d) SBAS-small redundancy with some images removed, (e) SBAS-high redundancy (f) SBAS-high redundancy with Jan-Feb images removed.Corresponding results for each strategy can be found in Figure 8.
In this study, we calculated 61, 63 and 62 differential interferograms for all three datasets captured by various Sentinel-1 geometries.To improve results in terms of pixel density and interferogram quality, we also took into account SAR signal quality.Therefore, as presented in section 4.1 any interferograms with low coherence were removed from the DInSAR processing chain (red squares in Figure 5(a)).For InSAR processing, ALOS-World-3D DEM was used to remove the topographic interferometric component.Goldstein Fast Fourier transform with a window size of 64 was used for phase smoothing and the Minimum Cost Flow (MCF) function with a regular grid covering the data extent was applied for phase unwrapping (PhU) (Costantini 1998).To enhance the signal-to-noise ratio (SNR) of the interferograms and to generate squared pixels, a multi-looking factor of 4 × 1 was applied.For the PhU, the threshold of 0.3 was set for coherence.Additionally, five various Ground Control Points (GCP) points located in stable areas and demonstrating highly coherent pixels were selected to estimate APS by using 3rd-degree polynomial estimation within refining and reflattening step implemented in ENVI SARscape™ software.

Persistent scatterer interferometry processing
In this study, we applied the classical PSI variant introduced by Ferretti, Prati, andRocca (2000, 2001).The PSI technique relies on phase information from single isolated objects called persistent scatterers (PSs), which are characterized by a high temporal phase stability (building, infrastructure objects) evaluated based on e.g. the amplitude dispersion index (<0.25).The main goal of this method is to compute differential interferograms for all images (secondary) with the same reference (primary) image (e.g.(Crosetto et al. 2016;Wasowski and Bovenga 2014).In PSI, the primary (φ m ) image is selected usually in the middle of the investigated time series to limit temporal decorrelation (Figure 5(b,c)) The use of multiple SAR images means that both the topographic interferometric components and the APS can be effectively estimated and removed.For the removal of the APS, high pass filtering in time with a size of 365 days in time and low pass filtering with a size of 1200 m in space was used.All secondary images (φ 0,1,2,..,n-2,n-1,n ) were coregistered to the one primary images (φ m ).For topographic phase removal, similarly as in DInSAR, ALOS World 3D DEM was used.PhU was performed in a sparse grid, which is feasible only in a multi-interferogram framework, where all phase data relative to the same target are jointly analyzed (Ferretti, Prati, and Rocca 2000).
ENVI SARscape™ -based implementation of PSI utilized a priory linear deformation model in time to properly unwrapped interferometric phase.This fitness between real defamation and modeled can be described by temporal coherence (TC) which can be presented as: where m is the number of SAR images, f i def is the phase component related to the displacement including modeled and un-modeled deformation and f i model is the model phase.The TC ranges between 0 and 1 and low coherence number indicate large unmodeled deformation and/or large phase noise.Having considered mining-related deformation behavior in time, in our processing we set TC on a threshold equal to 0.6 to better capture nonlinear deformation pattern.In the first try, we selected TC equal to 0.7 but it provide no information for CRs in the AoI, therefore we decreased that value.

Small BAseline subset
SBAS methods take the advantage of interferometric pairs connected based on the small temporal and spatial baseline.This method allows to capture distributed scatterers (DS) and therefore more information in the agricultural and vegetated areas.Interferometric processing chain of SBAS includes coregistration, interferogram generation, multi-looking filtering and PhU.The secondary images were coregistered according to the primary images with assistance of the DEM.However within the SBAS processing, many subsets represented by temporal baseline threshold are created and many primary images can be connected to secondary images and therefore provides redundancy of the observation as presented in Figure 5(d-f).In our study, 1st SBAS strategy created interferometric pairs between images with temporal baseline of maximal 24 days while in the 2nd and 3rd strategy interferograms were created between images with maximal temporal baseline of 66 days (in 3rd strategy images from Jan-Feb were additionally removed).In the coregistration step, (1) a local nonparametric shift estimate is computed by using orbital information and DEM (2) a set of windows (64 × 64 in range and 128 × 128 in azimuth) were selected on the primary image to compute the cross-correlation function and (3) the residual parametric shift was approximated and the appropriate shifts were used.Another processing steps, such as multi-looking, Goldstein filtering, PhU was carried with the same parameters as in conventional DInSAR processing.
To limit PhU errors caused by the low coherence, a decomposition level (DL) was exploited.The DL performs data multi-looking and under-sampling in an iterative way to unwrap the interferograms at a lower resolution and then reconstruct them back to the original resolution (ENVI SARs-cape™, technical description 2009, 2012).Refinement and reflattening step was performed based on the same GCP which were used in DInSAR.After PhU and refinement and reflattening step, in the first inversion step, the SBAS inversion kernel with periodical linear deformation model in time was implemented to provide the first displacement rate estimation and residual topography (Berardino et al. 2002).In the second inversion step, APS filtering (the same parameter as in PSI processing) was carried out based on previous results and was removed from the first-estimated results.Finally, the displacement time series for the TC were set as 0.2 and were geocoded according to the WGS84 reference system.

Decomposition of LOS deformation into vertical and easting deformation components
Displacements estimated with SAR are represented as 1D deformations along the LOS direction, made up of a combination of vertical (d u ) easting (d e ) and northing (d n ) deformation components as shown in Equation (2) (Hanssen 2001): where d LOS is the deformation along the LOS direction, a h is the heading (azimuth) angle, and u inc is the incidence angle.This means that a deformation analysis of a stack of LOS directions using InSAR is not able to fully capture the magnitude and direction of the displacements.Decomposition is therefore necessary.Mathematically speaking, it is possible to calculate d u , d e , d n from the LOS deformation based on at least three various geometries.However, the sensitivity of SAR sensors in terms of measuring d u , d e and d n can be represented as (Ng et al. 2011): Based on Equation (3), in the case of an ascending pair with an incidence angle of 38.35°and an azimuth angle of 81.77°, the coefficients are 0.78, 0.61 and 0.09 for the vertical, easting and northing deformations, respectively.The contribution of the northing displacement to the total LOS displacement is insignificant (0.09), and since there is a lack of variety in the geometric configurations of different SAR missions, it is very difficult to measure the northing deformation accurately (Ng et al. 2011).In view of the low d n sensitivity, only the vertical and east-west deformation components were therefore estimated by assuming north component as zero, even if we utilized three various viewing geometries.An interpolation of the deformation in time was also needed, since the image acquisition dates were different in each viewing geometry (with offsets of two and three days).When interpolating the deformation in time, we assumed a linear relationship between two neighboring deformation maps (e.g. 6 or 12 days).Finally, interpolate time series displacement where decomposed by Least Square Estimation (LSM) by using the following functional model ( 4): Using Equation (3) for data from each orbit, it was possible to transform to the functional model of the least squares method (4), we obtained (5): The final values of the unknowns were determined based on Equation ( 6):

Acquisition of auxiliary data using classical geodetic techniques
Five CRs were measured using GNSS static observation and leveling, and three independent GNSS measurement sessions were carried out every three days.More specifically, in the first measurement campaign, the first, second and third measurement sessions were carried out on the 2nd, 3rd and 4th of July 2018, and each session involved least 10 h of GNSS observations.Over the same three days, two measurement teams carried out leveling on the CRs.The same processes were repeated for the second measurement campaign carried out between 28th and 30th June 2019: three GNSS sessions were carried out, and leveling was conducted.The leveling measurements were performed in a geometric manner on 10 km of leveling lines with an accuracy of 0.2 mm/km.GNSS postprocessing calculations were performed using Leica Infinity software, based on precise satellite orbit and clock corrections: CODE final (Dach et al. 2018).The CR coordinates were determined by estimating independent vectors assigned to four reference points: KATO, TAR1, WOD1 and ZYWI.The reference GNSS stations belonged to the Polish European Position Determination System (EUPOS) network, known as ASG-EUPOS (www.asgeupos.pl).The results of the postprocessing calculations were obtained using an ETRF2000 reference frame.Permanent GNSS measurements were also carried out via RES1 permanent station, which was provided by a postprocessing service developed as part of the EPOS-PL project, based on Bernese GNSS Software v. 5.2 (Dach et al. 2018).This was used in conjunction with observations from the multi-GNSS satellite systems of GPS, GLONASS and Galileo.To ensure the highest quality of the solution, the final products (orbits, satellite clocks and Earth rotation parameters) were acquired from the International GNSS Service (IGS).The postprocessing service used a double-difference method to perform calculations on observations from 39 stations located in Poland and neighboring countries.Of these, 25 sites formed part of the IGS and EPN (EUREF Permanent GNSS Network) networks, and represented reference data sources.The remaining 14 stations (including RES1) were located in regions affected by mining deformation, i.e. the USCB.

Accuracy measures used to evaluate the results
Since height coordination estimated by using leveling and GNSS measurements are the same, we evaluated the InSAR results using GNSS and RES1 measurements.We evaluated our results in two stages.In the first stage, we evaluated the results internally in LOS geometry for each Sentinel-1 dataset.In this case, the reference data (GNSS and RES1) were reprojected onto the LOS direction and compared with InSAR measurements, based on DInSAR, PSI and SBAS results.The absolute errors were then used to calculate the RMSE according to the following equation: where REF may be either the conventional geodetic GNSS measurements on CRs or the RES1 permanent station measurements, respectively.InSAR may be either the DInSAR, PSI or SBAS processing results.RL is the relative orbit number.P is the number of validated observations.In the case of InSAR methods, this will be five CRs, respectively; in the case of RES1, it will be N−1, where N is the number of SAR scenes.
The second stage of evaluation was based on the decomposed results from the most reliable CR-InSAR results.Decomposition of the InSAR-derived LOS deformation map was carried out as described in Section 4.3, providing vertical (Up), and horizontal (easting) components.The RMSEs were calculated using the following equation: where DC is the specific deformation component (easting, vertical), represented as either as Up or east.L is the number of validated observations-five in case of CRs, in the case of RES1, it will be N −1, where N is the total number of SAR scenes in the three different geometries.

Evaluation of the InSAR signal quality using multi-temporal coherence
As mentioned in Section 4.1, we observed a significant loss of coherence in some specific periods, which were identified mostly at the end of January and February.More specifically, coherence decreased between the 10th and 28th of January 2018 for all geometries.The meteorological data (Figure 6(a-b)) indicates snow coverage of around 15 cm during this period; this makes it likely that snow accumulated on the CRs leading to a problem with the coherence between snow and snow-free acquisition.Figure 6(c,d) show an example of the results of the original coherence for a descending orbit (124) for the DInSAR and PSI method, respectively.Coherence loss during winter season but also in some individual cases can be observed.This non-winter coherence loss can be associated with another environmental disturbances (e.g.frozen rainfall, water, falling leaves).Figure 6(e,f) show the coherence values for this dataset after the removal of these bad images, using the DInSAR and PSI methods, respectively.This step allows to restore the coherence values for all InSAR processing techniques (DInSAR, PSI and SBAS).
The same analysis of signal quality was carried out for additional two Sentinel-1 datasets, acquired using ascending 175 and descending 51 geometries, and the final multi-temporal coherence is shown in Figure 7.It can be seen that a low coherence value for CR04 is represented for the whole time series for ascending geometry-175 (Figure 7(c,d)) and fluctuates around 0.6 and 0.4 for the DInSAR and 0.1 and 0.4 for PSI results, respectively.Since this situation appears only one geometry, it can be foreseen that some obscures occurred between the satellite LOS and CR04.This is likely to have been caused by the proximity of this CR to tall trees (see Figure 2(c)).

Ability to capture maximal deformation gradient by using various techniques and strategies of InSAR processing
The zoomed views of the LOS results captured for all three geometries and for all processing techniques and strategies are visualized in Figure 8 and superimposed on CR01 and CR02 location  (CR01-02 represent the biggest deformation gradient).Despite the removal of low-quality images, the PSI failed in the estimation of this fast deformation gradient in all strategies.In the first strategy when no images were removed (Figure 8(a-c)), there were no information for the center of the subsidence basins in the location of CR01 and CR02.In the second strategy, when some images were removed (Figure 8(a ′ -c ′ )), underestimation of the real deformation is observed (CR01 represents subsidence of approx.43 cm).The same situation appears for the third PSI strategy, when all winter images were removed (Figure 8(a ′′ -c ′′ )).In case of SBAS, the same underestimation of the real deformation appeared for two SBAS strategies with high interferograms redundancy without winter images removal (Figure 8(d ′ -f ′ )) and with high redundancy and winter images removal (Figure 8 (d ′′ -f ′′ )).Estimation of the reasonable deformation in the CR01 was only possible for SBAS with small interferograms redundancy determined by temporal threshold of max.24 days (Figure 8  (d-f)).In contrary, DInSAR succeed in the estimation of the maximum deformation gradient but also allowed to detect two separate subsidence basins (Figure 8(g-h)) which were not possible or not visible in SBAS results even with low interferometric redundancy.For a better representation of the problem, RMSE and the InSAR-estimated and ground truth deformation for the most deforming area-CR01 are presented in Table 2.
For a better representation of the active mining deformation phenomena, results for the whole AoI and for the strategies which allow to capture reasonable results on CRs (SBAS with low redundancy and DInSAR strategy) are presented in Figures 9-10.As can be seen, SBAS and DInSAR results in all geometries correspond with each other quite fairly, however by using SBAS method, clear detection of the border between the two subsidence basins is not possible.In DInSAR, the borders of the deforming areas can be easily observed and therefore differentiation of the two subsidence bowls can be observed.These NaN pixels in the area of the maximal deformation are mostly due to the nonlinear temporal deformation model in time.For a better representation of the problem, CR01 and CR02 time series displacement are represented in Figure 9(a ′′ -c ′′ ).It can be seen that the deformations for CR01 and CR02 were highly nonlinear (particularly CR01), with a large gradient (reaching 43 cm in descending geometry), and the deformation estimation using the PSI with linear model in time was therefore unsuccessful for these two CRs.

Vertical and horizontal deformation maps derived using DInSAR and SBAS
Deformations from all three geometries were decomposed into 2D deformation components (vertical and easting) for both techniques (DInSAR and SBAS).Figure 11 shows the vertical (a-a ′ , b-b ′ ) and horizontal (c-c ′ , d-d ′ ) components obtained from DInSAR (a-a ′ , c-c ′ ) and SBAS (b-b ′ , d-d ′ ) processing, respectively.The maximum deformation gradient, which was detected in one subsidence basins for the vertical component was -54 cm and was only detected with DInSAR method.The maximal subsidence that could be detected with SBAS was almost -41.8 cm.A similar situation arose for the easting deformation component: with the DInSAR method, the maximal detected horizontal movement was approximately 22 cm, while in SBAS, it was 17 cm.RMSE of decomposed results for the location of CRs and RES1 permanent station are presented in Table 3.

Processing obstacles of CR-InSAR
The high level of coherence between the SAR images directly affects the number of detected PS or DS and therefore InSAR product density.Placement of CRs in the field does not guarantee the acquisition of interferometric information.Firstly, CRs should be appropriately positioned and oriented considering the geometry of the satellites and some obscures which can result in backscattering stability (e.g.trees see Figure 3 and Figure 7(c-d)).Additional factors that are important include atmospheric conditions, and especially snow coverage, as this limits the possibilities of the practical use of CRs.Variations in multi-temporal coherence over time indicate that snow coverage has a strong influence in terms coherence loss.Temporary snow cover on CRs causes the coherence between images taken before/after and during the snow cover to decrease significantly.In this study, this effect was seen in some specific period of January and February 2019.The approach which removes all January-February images is not the best approach in Polish conditions, where snow coverage is not continually preserved in the whole period.Moreover, removal of so  many images can lead to significant loss of coherence but also to the PhU problems caused by phase aliasing due to the large deformation gradient.Although, the quality of the solution can be improved by excluding only these low-coherence images.Similar problem connected with long temporal baseline which corresponds to the significant deformation gradient and PhU problems can be associated with the SBAS approach with high interferograms redundancy.Deformation estimates delivered from the interferograms calculated with longer temporal baseline (e.g. 60 days) are leading to the PhU problems and underestimation of the real deformation.In contrary, based on Table 2, we can observe that SBAS method with high redundancy provides the most accurate SBAS results for areas with small deformation gradient (RES1).Based on the achieved results it is advisable to utilize high network but for the phenomena with small and moderate deformation gradient.The main limitation of the PSI implementation in ENVI SARscape™ software is only one deformation model in time (linear) which for mining phenomena is not appropriate.This led to the situation that in the case of the PSI technique, it was impossible to determine the deformation values for CR01 and CR02.Periodically linear model in time used in SBAS processing partially addresses that problem, but for better representation of the mining deformation another deformation model in time needs to be tested in the future.
Considering the InSAR final product density and accuracy for mining application, only two approaches were considered reasonable -DInSAR and SBAS with low redundancy network.The great advantage of DInSAR is that its results allowed to observe two separate subsidence basins which were not visible in the results of any other technique due to the missing PS/DS.In contrary, the accuracy of DInSAR estimations was lower than SBAS with RMSEs of 27 and 34 mm for the East and Up directions, respectively.Additionally, for a high deformation gradient, the vertical deformations are usually underestimated in any InSAR results, and this was especially true in the case of CR01 when compared with classical geodetic measurements.

Time series variations in the PSI and DInSAR products observed in relation to the GNSS RES1 permanent station
Figure 11 shows that both InSAR results generally correspond to each other fairly well in the area with a small deformation gradient.However, the deformations estimated from the SBAS results are smoother due to the temporal filtering of APS, while DInSAR is biased by atmospheric influence, which was only eliminated using polynomial trend removal.
There is good agreement between InSAR estimated deformation and GNSS in RES01 (Figure 11).The strongest correspondence can be observed between the PSI (with the exception of PSI strategy with the removal of winter images) all SBAS strategies and GNSS results, for which the difference between time series was less than 5 mm over the whole period of analysis (Table 3).The DInSAR results for the descending orbit also correspond well with the GNSS measurement, although less strongly than SBAS.The results for the ascending orbit showed a moderate consistency between the DInSAR and GNSS results, while those from SBAS showed fairly high  consistency.It can therefore be concluded that some of the errors that accumulated in the ascending DInSAR results influenced the time series.

Conclusions
Deformations caused by underground coal exploitation are characterized by high nonlinearity, large deformation gradient, and, very often, loss of coherence.The application of CRs in such areas may therefore allow for the restoration of stability to the radar signal.However, monitoring such ground surfaces is challenging, and it is necessary to use a method that can determine the deformation scale over large areas with high spatial and temporal resolution and with high accuracy.Various SAR techniques have been developed, although they have certain advantages and limitations in terms of APS removal, PS/DS density or accuracy.
Another obstacle was related to weather conditions, in particular snow coverage and unsuitable locations of CRs close to tall trees, which destroy double-bounce radar scattering and decrease the coherence of the CRs.While the impact of snow coverage can be eliminated by the removal of snow-covered images, limitations related to incorrect locations of CRs, such as close proximity to tall trees, could not be mitigated and need to be considered beforehand.Additionally, in the case of monitoring large motion, it is advisable to remove only images which destroy signal stability, not all winter images.This allows to restore the coherence and minimize temporal baseline and also facilitates the PhU step, which for significant deformation gradients and short wavelength can be challenging.
The analyses performed in this study, show that the installation of artificial CRs in undeveloped and heavily vegetated areas allows radar signal stability to be restored and a high level of coherence between SAR images to be achieved.However, we observed significant obstacles in matching linear (PSI) and partially linear (SBAS) deformation model with real mining displacement behavior.Linear model in PSI fails in the estimation of fast and nonlinear deformation while SBAS method with periodically linear model partially addresses that problem.The advantageous approach proved to be conventional consecutive DInSAR without any a priori deformation model, however the accuracy is lower when compared with other InSAR strategies.However, APS removal in PSI and SBAS method provides more accurate results and smoother time series variations that show a stronger correspondence with the RES1 GNSS permanent measurements than those obtained using DIn-SAR.Therefore, it was concluded that DInSAR and SBAS with small redundancy should be jointly applied in monitoring such deformation patterns and there is still a need for more dedicated approaches and deformation models to better representation of real mining-related deformation.This will fulfill the EGMS gap and therefore will be beneficial for the mining sector.

Figure 1 .
Figure 1.Deformation over Rydułtowy Mine presented by EGMS (accessed 11.10.2022):(a) general location of the area in EGMS; (b) visualization of the deformation estimated by using PSI in that region (blue polygons marked indicate areas where many PS candidates are located but are not taken into the final processing due to the model imperfection (linear vs. nonlinear) (c) timeseries representation of the persistent scatterer with the maximal detected deformation rate (white circle on fig.1b).

Figure 2 .
Figure 2. Area of study: (a) general location of area of interest (AoI); (b) location with the surrounding cities; (c) AoI and the CRs and the GNSS permanent station.

Figure 4 .
Figure 4. Methodology applied in the presented study ('imgs' is abbreviation of images).

Figure 6 .
Figure 6.Example of coherence time series for descending 124 orbit in case of (c,e) DInSAR processing and (d,f) PSI processing.Subfigures (c,d) show the original processing results, and (e,f) results after removal of specific SAR images.Subfigures (a) and (b) represent weather information.

Figure 8 .
Figure 8. Deformation estimated for all InSAR processing techniques with various strategies.Red color indicated the strategies which allows to capture reliable (not underestimated) information for all five CRs.

Figure 11 .
Figure 11.Deformation estimates in the LOS direction for the RES1 station for (a) 51 descending, (b) 124 descending, and (c) 175 ascending orbits.Figure 11(d-f) presents northing, easting and vertical deformation component in comparison to the RES1 station, respectively.
GNSS campaigns on CRs Time spans of the first and second measurements 2nd to 4th July 2018, 28th to 30th June 2019 Leveling data Time spans of the first and second measurements 2nd to 4th July 2018, 28th to 30th June 2019 Meteorological data Station Rybnik Source https://danepubliczne.imgw.pl/Type of data precipitation Station Katowice Source https://www.ogimet.com/Type of data snow cover 51 ) − sin (u inc 51 )cos(a h 51 ) cos(u inc 124 ) − sin (u inc 124 )cos(a h 124 ) cos(u inc 175 ) − sin (u inc 175 )cos(a h 175 )

Table 2 .
Real and InSAR estimated deformation for CR01 together with RMSE calculated for all CRs (NaN indicates that not for all CRs it was possible to capture information by this particular InSAR strategy) and RMSE for RES1 permanent station.Bold font indicate the processing strategy when not significant underestimation was observed (phase unwrapping errors).

Table 3 .
RMSE calculated for CRs and RMSE for RES1 permanent station for the DInSAR and SBAS results after decomposition.