ANALYSIS OF MINING-INDUCED TERRAIN DEFORMATION USING MULTITEMPORAL DISTRIBUTED SCATTERER SAR INTERFEROMETRY

. This work addresses a methodology based on the Interferometric Synthetic Aperture Radar (InSAR) to analyse and monitor ground motion phenomena induced by underground mining activities, in the Legnica-Glogow Copper District, south-western Poland. Two stacks of ascending and descending Sentinel-1 Synthetic Aperture Radar (SAR) images are processed with a small baseline multitemporal approach. A simple method to select interferograms with high coherence and eliminated images with low redundancy is implemented to optimize the interferogram netwrork. The estimated displacement maps and time series show the effect of both linear and impulsive ground motion and are validated against Global Navigation Satellite System (GNSS) measurements.


INTRODUCTION
Underground mining activity often triggers substantial ground surface displacements both vertical, i.e. subsidence, and horizontal directions. This work is devoted to analyse such ground motion phenomena using satellite radar interferometry. Interferometric Synthetic Aperture Radar (InSAR) has been employed to study terrain deformation for decades (Massonnet and Feigl 1998) and has seen relevant developments in terms of accuracy and coverage with the introduction of techniques based on persistent scatterers (Persistent Scatterer InSAR, PSInSAR) and distributed scatterers (DSInSAR) (Ferretti et al., 2000(Ferretti et al., , 2001Van Leijen, 2016;Crosetto et al., 2016). The availability of SAR datasets with increasing spatial and temporal coverage, and decreasing temporal intervals between two subsequent acquisitions, such as the ones collected by the Copernicus Sentinel-1 mission, gives the opportunity to analyse and monitor ground-surface deformation phenomena, of natural origin (such as landslides) or man-induced. In this work, we employ a multitemporal DSInSAR technique on a satellite SAR dataset collected by the Sentinel-1 A&B sensors, using the interferometric wide swath (IW) acquisition mode and both descending and ascending passes. relevant challenges have been faced in developing robust PS solutions due to the decorrelation between SAR complex images, resulting in poor interferometric phase stability. Several approaches have been proposed to decrease the decorrelation among interferometric phases, such as interferogram filtering (Goldstein and Werner 1998), or the use of multitemporal techniques, which form a network of multi-baseline interferograms aimed at decreasing the effect of noisy interferograms (Pepe et al., 2015). The area of interest (AOI) is located around the town of Polkowice, south-western Poland, which is located within the Legnica-Glogow Copper District (LCGD). The presence of copper mines within the AOI causes ground-surface deformations, with both a vertical component (subsidence) and horizontal component, detectable in the East-West direction. (Hejmanowski et al., 2020;Witkowski et al., 2021;Malinowska et al., 2018;Antonielli et al., 2021;Mazzanti et al., 2021). The geographical extension of the area affected by terrain deformation varies depending mainly on the depth of the excavation room, with deformation cone footprints ranging from a few hundreds of meters to kilometers. On the other hand, the * corresponding author: rpalama@cttc.es temporal deformation behaviour changes the deformation rate caused by small collapsing events or heavy mining-induced seismic tremors, which are considered the most serious threat to the built environment due to their impulsive nature The results show the presence of seven main areas that are affected by ground surface deformation. The UD component of the linear velocity reaches maximum values of about 40-50 mm/year, common to all the affected areas, whereas the EW horizontal component is maximum for the area at the south of Polkowice (about 30 mm/year eastward). This paper is organized as follows: section 2 illustrates the dataset and study area, section 3 describes the employed methodology, section 4 illustrates the main results, main remarks and future developments are discussed in section 5.

DATASET AND STUDY AREA
The proposed method was tested on SAR SLC images collected by the Sentinel-1 A/B on-board sensors, using the Interferometric Wide Swath (IW) acquisition mode. Vertically polarized transmit-receive (VV) data were analysed, as they are expected to provide a more extended DS coverage with respect to (w.r.t.) the other available polarization (vertical-transmit horizontalreceive, VH). The area of interest (AOI) is located around Polkowice, south-western Poland. The area consists of urbanized environments, rural areas, bare soil, forested and vegetated land. This mining activity causes deformation that extends within the south-east area of Polkowice, also affecting the nearby villages of Biedrzychowa and Pieszkowice. In order to study the deformation phenomena affecting the AOI, we processed the SAR images collected at both ascending and descending passes, which makes feasible the estimation of the horizontal and vertical deformation components. The analysed datasets were collected by the Sentinel-1 on-board radar system from September 2019 to September 2020, with an interval of 6 days between two subsequent acquisitions. Sentinel-1 A/B SAR SLC images are available at the European Space Agency (ESA) Open Access Data Hub (https://scihub.copernicus.eu) and NASA's Alaska Satellite Facility platform (https://search.asf.alaska.edu). A network of GNSS stations is also located within the AOI, which provides robust measurements of the three-dimensional displacement vector.

METHODOLOGY
The main processing steps leading to the formation of the DS grid consist of (i) image co-registration, where all the SAR images are coregistered against the selected super-master (SM) image; ; (ii) interferogram generation, achieved using a multitemporal network, with temporal baselines ranging from 6 to 42 days; (iii) multilooking, using a window whose azimuth and range lengths are 2 and 10, respectively, yielding a pixel size of about 28 m in azimuth and 23 m in range; (iv) coherence computation, i.e. for each interferogram, the interferometric coherence (Hanssen, 2001) is computed using a 2-by-2 window; (v) interferogram selection: the interferograms generated in step 2 form a multitemporal network where each master image is associated with each slave image through an edge which is rep-resented by the interferogram itself. Typical interferometric techniques make use of the full interferogram network, formed by all the interferograms whose temporal baseline falls between predefined minimum and maximum values. In this work we adopt an interferogram selection procedure to eliminate low-coherence interferograms, based on constraints on the minimum total coherence and minimum image redundancy. All the interferograms whose average spatial coherence is lower than a defined threshold are eliminated. Then, all the images whose redundancy is lower than a minimum value are also eliminated. In fact, removing some images causes the elimination of those interferograms they are involved into, which makes necessary an iterative procedure to converge to a reduced interferogram network, for which the two conditions, i.e. minimum average coherence and minimum image redundancy, are met. (v) DS selection based on a temporal coherence metric (Pepe et al. 2015) which is calculated as where the index P identifies a pixel, M is the number of interferograms selected in the previous step, k is the interferogram index, φ k is the wrapped interferometric phase of the k-th interferogram and φ k,LP is the low-pass estimate of the wrapped interferometric phase of the k-th interferogram. The calculation is performed on the wrapped interferograms. First, the estimation of the spatial low-pass interferometric phase is computed using a boxcar averaging window (in this work we used a 2-by-2 window), then the spatial low-pass phase is subtracted modulo-2π from the original interferometric phase, giving a high-pass estimate of the interferometric phase. For each pixel, the high-pass phasors are averaged coherently, in order to mitigate the effect of the phase noise, then the module of the resultant vector is computed. The values of the equivalent temporal coherence are included between 0 and 1. The higher the value of the equivalent temporal coherence, ( ) , the more coherent a pixel. This multitemporal point selection criterion highlights the overall behaviour of the interferometric phases of a pixel, ignoring isolated noisy interferograms, that can otherwise hinder the selection of a highly coherent pixel.
(vi) phase unwrapping; (vii) time series generation: this step is described in (Devanthéry et al., 2014a) and produces a 2D (spatial) plus 1D (temporal) grid of DS points.; (viii) filtering of the Atmospheric Phase Screen, which is previously estimated using a combination of spatial low-pass and temporal high-pass filters (Palamà et al., 2021); (viii) geocoding. This procedure was applied to both ascending and descending data, thus enabling the extraction of time series of the East-West (EW) and Up-Down (UD) displacements. The software used for image coregistration, interferogram generation, coherence computation, phase unwrapping and time series generation is contained within the CTTC Persistent Scatterer Interferometry processing chain (PSIG) (Devanthéry et al., 2014a(Devanthéry et al., , 2014b.

RESULTS
This section is devoted to the description of the main results obtained using the processing chain explained in Section 3 on the SAR SLC datasets listed in Section 2. It should be noted that the analysed area represents a sub-area of one burst for both ascending and descending data. For the ascending data, we selected an area of 970 range (rg) bins by 720 azimuth (az) bins, whereas for the descending data the area was 810 (rg) by 730 (az).

Interferogram and point selection
The initial interferogram network was generated using temporal baselines from 6 to 42 days, separated by an interval of 6 days, i.e. the revisit time of Sentinel-1 constellation. The interferogram network is represented in this paper using a graph, where each image is associated with a point in the twodimensional plane given by the temporal and perpendicular baselines with respect to the SM image (i.e. the image taken as reference for coregistration). The graph of the full network for the Sentinel-1 ascending datasets is shown in Figure 2a. The original network consists of 399 interferograms involving 62 images, with a maximum redundancy of 14 achieved in the middle of the monitored period. The interferogram selection method described in section 3 yields a reduced network of 298 interferograms for 61 images, shown in Figure 2b. This result was obtained using a minimum average coherence of 0.2 (typical values are included between 0.15 and 0.55) and a minimum image redundancy of 4. These values were tuned by a trial-error procedure to eliminate the maximum number of interferograms while keeping as many images as possible, in order not to decrease the temporal resolution for the time series analysis. Figure 3 shows the histogram of the temporal coherence values for the ascending and descending data obtained using the full and reduced interferogram networks. We observe that using the reduced interferogram network increases the values of the temporal coherence, which yields a larger number of selected points (the threshold used for point selection is equal to 0.8).

Vertical and horizontal deformation components
The UD and EW deformation components extracted from the ascending and descending LOS displacement datasets are shown in Figure 4, displaying strong deformation rates. The UD component reaches maximum deformation values of about -130 mm in the Rudna mine (Polkowice south, marked with 1). Regarding the EW component, it is characterized by an area with positive values moving eastward (maximum deformation around 60 mm) which is close to an area displaying negative values moving westward (maximum deformation about 40 mm). A similar behaviour, with significantly stronger deformation values, is observed in the area de-limited by the towns of Przesieczna and Borów (about 7-8 km north-west of Polkowice, area marked with indexes 3 and 4). Two deformation spots moving in opposite directions in the EW deformation map are also noticed.

Comparison with GNSS data
In this section the LOS displacement time series are compared with the displacement time series measured by four selected GNSS stations projected onto the LOS directions. The GNSS device measures the 3D displacement vector along the vertical (height), latitude (East-West) and longitude (North-South) directions, which are projected onto the ascending and descending LOSs. This procedure allows a basic validation of the obtained LOS displacement time series, as the 3D GNSS measurements provide for a robust benchmark when projected onto the LOS direction. The locations of the GNSS devices in the AOI are summarized in Table 2 and indicated in Figure 1.
In this work we analyse four different GNSS measurements, collected by the PIES (Pieszkowice), WYZY (south of Polkowice), ZHPR and ZHSR (located within Polkowice urban area) stations. Comparisons between DSInSAR results and GNSS LOS-projected time series are shown in Figure 5 for both ascending and descending trajectories. It should be noted that GNSS data are available from October 2019, which means that this validation is carried out during an 11-month period, as Sentinel-1 data were analysed until September 2020. We observe that the best match is achieved for the WYZY station ( Figure 5 g-h) which shows a quite significant deformation of the order of 120 mm. The WYZY time series highlights the presence of a slow deformation trend of about 5 mm/month for both ascending and descending data. A large discontinuity interrupts both the GNSS and DS time series in July 2020. Noticeably, the ascending DS time series shows a small discontinuity in May 2020, which causes a slight increase in the deformation rate, corresponding to a minor seismic event recorded at the beginning May 2020. A quite good agreement between DS and GNSS time series is achieved for the ZHPR station ( Figure 5 e-f), especially for ascending data, that shows similar deformation rates, whereas the descending DS data and GNSS data differ for about 9 mm at the end of the monitored period. A slight increase in the descending deformation rate is noticed in May 2020, corresponding probably to the seismic event recorded in the same period. The Pieszkowice (PIES) station shows good agreement for ascending data (deformation rate of about 3.6 mm/month). The descending data show a linear deformation rate faster than the ascending data (about 6 mm/month) until a large discontinuity occurs in July 2020. The ZHSR station shows a good match for descending data (Figure 5b), whose deformation rate is about 3.8 mm/month. A small discontinuity (smaller than 10 mm) is recorded by the GNSS data in July 2020, but it results smoothed in DS data. The ascending data show similar shapes of the deformation curves between DS and GNSS data, but significantly different deformation rates (0.9 mm/month and 3.6 mm/month for DS and GNSS data, respectively). The observed deviations can be caused by a not-perfect overlapping of the detected DS points and GNSS locations and also by the smoothing effect of the APS filter.

CONCLUSIONS
This work has addressed a DSInSAR method to analyse mining induced ground motion phenomena, implementing a simple interferogram and point selection technique. The interferogram selection method applies a dual condition on the minimum average coherence of each interferogram forming the network and the minimum redundancy of each image in the analysed stack. The proposed method was applied to two stacks of Sentinel-1 SAR SLC images in descending and ascending orbits, to study the deformations caused by underground mining activity in the Legnica-Glogow Copper District, Lower Silesian region, Poland. The adopted method has shown good performance in detecting both slow and fast deformations. The results were validated against in-situ GNSS measurements. Future developments will include implementing an interferogram filtering stage (Pepe et al., 2015) before the interferogam selection stage, as well as analysing the quality of the produced results using proper indicators  and implementing time series classifiers .