Ground Stability Monitoring of Undermined and Landslide Prone Areas by Means of Sentinel-1 Multi-Temporal InSAR, Case Study from Slovakia

Multi-temporal synthetic aperture radar interferometry techniques (MT-InSAR) are nowadays a well-developed remote sensing tool for ground stability monitoring of areas afflicted by natural hazards. Its application capability has recently been emphasized by the Sentinel-1 satellite mission, providing extensive spatial coverage, regular temporal sampling and free data availability. We perform MT-InSAR analysis over the wider Upper Nitra region in Slovakia, utilizing all Sentinel-1 images acquired since November 2014 until March 2017. This region is notable for its extensive landslide susceptibility as well as intensive brown coal mining. We focus on two case studies, being impaired by recent activation of these geohazards, which caused serious damage to local structures. We incorporate a processing chain based on open-source tools, combining the current Sentinel Application Platform (SNAP) and Stanford Method for Persistent Scatterers (StaMPS) implementation. MT-InSAR results reveal substantial activity at both case studies, exceeding the annual displacement velocities of 30 mm/year. Moreover, our observations are validated and their accuracy is confirmed via comparison with ground truth data from borehole inclinometers and terrestrial levelling. Detected displacement time series provide valuable insight into the spatio-temporal evolution of corresponding deformation phenomena and are thus complementary to conventional terrestrial monitoring techniques. At the same time, they not only demonstrate the feasibility of MT-InSAR for the assessment of remediation works, but also constitute the possibility of operational monitoring and routine landslide inventory updates, regarding the free Sentinel-1 data.


Introduction
Slope deformations are the most significant geohazards in Slovakia which annually cause an extensive economic damage, seriously limit the rational use of land and in rare cases also threaten human lives. This is conditioned by the complex geological setting of Slovakia, especially in Flysch regions and periphery of Neovolcanic regions, where landslides threaten up to 60% of the territory. Slopes within Neogene and Paleogene depressions are also substantially affected. Moreover, the activation or re-activation of relatively large amounts of landslides has been witnessed in recent years. Since 2010, more than 700 new slope failures have been registered and their activation was driven mainly by climatic anomalies, such as extraordinary rainfalls and melting snow cover. Many of these landslides currently represent a direct threat to the lives, health and property of the residents in the affected areas [1]. This fact has initiated their inclusion in the monitoring system built within the scope of the project "Partial-monitoring system-geological factors" [2,3]. The stability condition of these sites is being observed by systematic geological, geotechnical (groundwater table level, borehole inclinometers), and geodetic (GNSS, levelling) measurements, as well as regular evaluation of climatic factors. All those conventional in-situ terrestrial monitoring methods, however, provide only point-wise information about temporal evolution of deformation phenomena, while requiring laborious and regular measurement campaigns. Another geohazard posing a threat for the residents of Slovakia is a subsidence due to undermining, related to the present or historical exploitation of brown coal, various ores or salts. Brown coal mining is still in operation in the Upper Nitra region in Central Slovakia, where a collapse due to undermining has led to the abandonment of two municipalities, or their parts-Laskar and Kos. Coal mining in the Upper Nitra region has been operated by the Hornonitrianske bane (HB) Prievidza Mining Company, which has also been monitoring the undermining subsidence using the classical geodetic survey at the surface and extensive underground monitoring by the methods of mining geodesy.
Monitoring of mass-wasting geohazards like landslides or land subsidence due to undermining has been marked by new perspective from remote sensing, especially thanks to satellite-borne synthetic aperture radar (SAR) interferometry (InSAR), which successfully attains geodetic precision at extended spatial coverage [4]. Multi-temporal SAR interferometry techniques (MT-InSAR) are nowadays a well-developed tool for ground stability monitoring. They aim to overcome the major limitations of conventional differential InSAR, employing the time series of SAR images to identify the radar targets, which do not change their backscattering characteristic over time and remain phase coherent. Those are referred to as persistent scatterers (PS) and generally correspond to man-made structures as well as stable natural reflectors (e.g., bare rocks or outcrops). For such points, temporal evolution of deformation can be retrieved, thus forming a natural, opportunistic geodetic network to monitor seismic activity, landslides, cities, or even individual structures [5][6][7]. This application capability has recently been emphasized by European Space Agency's Copernicus operational programme of the Sentinel-1 satellite mission, providing guaranteed and freely available SAR images of the entire Earth's landmass within a 6 day repeat cycle.

Study Area
The aim of this study is to test the feasibility of MT-InSAR technology within the Upper Nitra region, Central Slovakia, exploiting all to-date available Sentinel-1 data. The investigated territory spans approximately 440 km 2 -22 km in the longitudinal direction and 20 km in the latitudinal direction ( Figure 1). This site is notable for its landslide susceptibility, as well as extensive coal mining activity. The mass movement processes are generated by the specific geological setting of the Vtacnik mountain range. Rigid volcanic rocks made of epiclastic volcanic conglomerates and breccias, pyroclastics and lava flows of Older to Middle Sarmatian age, form the summit parts and flanks of the mountain range [8]. They overlay softer complexes of Neogene sediments, prevailingly claystones, sandstones and siltstones, intercalated by brown coal seams. Slope deformations have a strong representation here, with almost all basic types of slope failures. In the summit parts of the mountain range, block ridges are present. Down the slope, detached blocks form block fields, and towards the Upper Nitra depression, landslides have evolved. Since Badenian complexes are known for rich deposits of brown coal, the unfavourable slope stability conditions have been impaired by underground mining, which at several sites has led to the acceleration of mass movements or to the subsidence of the territory above the mined-out underground spaces [9,10]. We closely investigate the above-mentioned geohazards on two case studies. The first study focuses on the suburbs of Prievidza city: Hradec, Velka and Mala Lehotka (Case study 1: Figure 2). The 2012-2013 reactivation of slope deformations in Hradec and Velka Lehotka caused serious damage to local infrastructure and buildings (Figure 3), leading to a declaration of an emergency situation. Consequently, landslides survey [11,12] and remediation [13] were implemented as inevitable steps, accomplished in 2014 and followed by the establishment of a monitoring network [3]. The remedial works mostly involved underground drainage by horizontal boreholes, combined with sand piles, drainage trenches and surface drains.   Figure 4) and spatially devastates the countryside. Since 1950, a building closure was declared in the north-western part of the Kos village, and mining began in 1988. Notable negative influences of mining are related not only to the ecosystems of the original landscape, but also to demographic changes due to the evacuation and withdrawal of former inhabitants. Population of the village has since decreased from the original 3500 to today's approximately 1200 inhabitants. The undermined areas of ≈158 ha are affected by strong subsidence movements, and the formation of several depressions (over area of ≈58 ha) is evidenced. Number of slope failures are induced over the undermined fields, which are continually subsiding during the mining activities. Therefore, the location, total area and depth of the depressions are changing dynamically over time. The predicted maximum decrease is in the range of 4-7 m with 3-5 m depth of the groundwater table level, what has resulted for example in formation of wetlands [15]. The 3D model of the respective brown coal deposit was created on the basis of mining data of the HB Company, a.s. within the project [16]. The subsidence of the area, which was reflected in the creation of ponds, was monitored via satellite and aerial images acquired in years 1987, 2001, 2004, 2006 and 2007. In 2007, these areas were surveyed by geodetic methods, revealing that the subsidence had begun within 1-2 years after the coal seams were extracted. The extent of the worked out spaces and related sinkholes, often indicated by lakes and wetlands, is depicted in the Figure 4b.

Sentinel-1 Dataset
We performed MT-InSAR ground stability analysis over the wider Upper Nitra region, utilizing 150 C-band Sentinel-1 SAR images per both of its ascending (78) and descending (72) orbits, which cover the time period from November 2014 till March 2017. Figure 1 depicts subsets of Sentinel-1 IWS scenes for ascending track No. 175 and descending track No. 51 over the area of interest (AOI). Due to very small perpendicular baseline distribution secured by stable Sentinel-1 orbital tube, as well as regular temporal sampling of acquisitions, forming single-master image combinations provides stable and rigorous geodetic time series [17]. Interferometric image combinations employed in processing are represented by Figure

MT-InSAR Processing
We incorporated the processing chain (please refer to flowchart in Figure 6) based solely on open-source tools. First processing segment was carried out under current Sentinel Application Platform (SNAP), being developed for ESA by Array Systems Computing under GPL (free, libre and open-source) license [19]. Thanks to its graph processing tool utility, we wrote series of scripts to automate the processing of whole stacks of Sentinel-1 SLC images in batch fashion.
Regarding the large size of original SLC data and associated high computational demands, the prime step involved the extraction of bursts covering the AOI. Then, the orbit state vectors of individual scenes were refined using Sentinel-1 precise orbit ephemerides, published by ESA on [20]. Sentinel-1 operates in a unique acquisition mode TOPS (Terrain Observation by Progressive Scans) [21], which captures images consisting of series of bursts with mutual overlaps. This provides wide swath width with a reasonable spatial resolution of 5 m × 20 m. Sentinel-1 TOPS SLC images, however, require a specific coregistration approach regarding the correct burst stitching (according to [22], an azimuth coregistration accuracy of 0.003 samples is required to achieve a phase error of less than 10 degrees). The main advantage of SNAP inclusion is its current state-of-the-art performance in TOPS coregistration. Note that TOPS coregistration operator in second step (consisting of back-geocoding and enhanced-spectral diversity correction) was carried out burst-wise, in order to make use of the entire burst overlap as well as retaining a sufficiently coherent area. Coregistration was assisted by SRTM 1 arc second digital elevation model (DEM) [23]. The AOI was subsequently cropped from de-bursted stacks, which was followed by differential interferogram formation, again employing SRTM 1 DEM. The output data are SLC stacks and differential interferograms coregistered onto a common master image frame, plus a geocoding mask. Geocoding based on orbital data and DEM is limited to positional accuracy within several metres, depending on the resolution of SLC data. Multi-temporal PS processing was carried out in StaMPS (Stanford Method for Persistent Scatterers) software implementation [24]. Standard PSInSAR technique (originally presented by [5]) selects PS candidates based on their phase variation in time, presuming a temporally linear deformation model; whilst StaMPS methodology exploits the spatial correlation of their phase without prior assumptions about its temporal nature [25]. Therefore StaMPS succeeds in the identification of less bright scatterers with lower SCR (signal-to-clutter), and thus provides denser PS coverage even in sparsely urbanized and natural terrain. The processing procedure followed steps, closely described in [24,26,27]: • Apriori PS candidates selection based on the 0.4 amplitude dispersion threshold.

•
Iterative phase stability estimation of PS candidates (analyzing their phase noise term). • PS selection based on ensemble temporal coherence in a probabilistic fashion.

•
Elimination of phase contribution due to residual topography (DEM error).
Estimation of spatially-correlated nuisance terms (atmosphere and orbital phase contributions).
Careful selection of processing parameters was required to avoid aliasing of the deformation signal, especially the size of resampling grid for unwrapping and window sizes of temporal and Goldstein filter [28] applied prior to unwrapping. The unwrapping procedure was followed by thorough detection and correction of possible phase unwrapping errors with regard to phase time series, as illustrated by e.g., [29]. Such errors occurred where the deformation pattern was not smooth.
After reasonable mitigation of phase unwrapping errors and spatially-correlated error terms, some atmospheric signal still persisted in the estimate. To date, several strategies for the elimination of tropospheric phase delays were developed, e.g., employing weather models, GNSS measurements or spectrometer data [30,31]. However, these are often limited by the spatio-temporal resolution or a non-simultaneous acquisition of sensors. With respect to the smaller extent of our AOI, we estimated the tropospheric phase screens and their stratification for each acquisition date empirically from the relationship between the interferometric phase and the topography, as proposed by [32] and implemented in Toolbox for Reducing Atmospheric InSAR Noise (TRAIN). It is important to note that strong tropospheric turbulence, albeit not expected in our AOI, cannot be accounted for using only phase-based correction methods and would leak into the residual phase as an incorrect signal [32].
Finally, unwrapped residual phase time series were converted to displacement via sensors wavelength. Mean annual velocities of displacement were then estimated in the least squares sense by inversion of the linear movement model.

Results and Discussion
The basic and most commonly used representation of MT-InSAR analysis are ground deformation maps, depicting mean line-of-sight (LOS) annual velocities of individual PSs. Concerning the reliability of the estimate, its indicator is ensemble temporal coherence (as a maximum of periodogram). Hence the most common, albeit harsh practice for outlier elimination is thresholding on ensemble temporal coherence only. Though in order to preserve as much information as possible, we applied outlier elimination approach proposed by [33], which allows for the minimisation of outliers in final results while preserving spatial and statistical dependency among the observations. This approach preserved a valuable amount of 28,101 points for the ascending track and 25,957 points for the descending track (representing approximately 25 percent of whole point-cloud), that would have been otherwise discarded by standard coherence thresholding. Outlier-free mean LOS velocity deformation maps for both ascending and descending datasets are in Figure 7. Those are relative to a common reference area (depicted by a black rectangle). Reference area is situated in the old central part of Prievidza city. Its stability was discussed with geologists and verified in-situ. Although the small values of standard deviations may seem outstanding, they actually correspond to precision, not accuracy. Precision is strongly influenced by the number of images, i.e., the degree of freedom of the estimate. Higher standard deviation does not necessarily correspond to noisier scatterers, as non-linear deformation could also push this value up [4]. To assess the accuracy of the results, calibration methods employing artificial corner reflectors or compact active transponders are required [34].
We achieved solid coverage over urban and peri-urban areas, yet several clusters of points were successfully retrieved even in rural parts of the AOI (amidst arid fields or outcrops). Deformation maps provide a global outlook on the relative movements in the whole area of interest, and thus show that no large spatial wavelength deformation signal is present. At the same time, they demarcate distinct small-areas of detected spatially clustered motion. We focus on the two already introduced and most significant geohazards within the Upper Nitra region.

Case Study 1: Landslides
InSAR results reveal that Hradec and Velka Lehotka (suburbs of Prievidza city) are affected by severe slope deformation processes. While those are undergoing movement towards the satellite on ascending track (Figure 8a), respective PSs exhibit movement away from the satellite for a descending track results (Figure 8b). This feature is a consequence of prevailing horizontal direction of actual displacement vector, as illustrated by Figure 9. Interpretation of surface movements is generally challenging when dealing with LOS displacement velocities. Assuming major sensitivity to vertical displacement and ignoring the presence of a potential horizontal one can introduce large errors into the final estimates if using mere projection via incidence angle [35]. Moreover, Sentinel-1 TOPS acquisition look angles are way larger than those of legacy ERS or Envisat satellites, what further expands this vulnerability. Yet, thanks to the availability of both acquisition geometries, LOS velocities can be further decomposed into vertical component and horizontal component in descending azimuth look direction (Figure 9).
Considering that ascending and descending stacks are related to different reference frames and individual PSs might actually correspond to different structures due to the strong directional dependence of the scattering mechanism, such decomposition encounters problem of PS pixel-matching. Hence we apply a simplified, three step approach: 1. Mean LOS velocities are resampled into regular 50 m × 50 m grid. 2. The grid is subsequently filled by the weighted average of all PS's velocities allocated within individual cells. Temporal coherence is assumed a weighting factor. 3. Finally, decomposition is evaluated using the system of equations in accordance to [36].
Its outcome is visualised by Figure 10. The detected movement activity is in a good agreement with the slope failure delineation and its direction mapped by [14]. Minor discrepancies occur only in the northern part of Velka Lehotka village (the southern one). In 2013, several buildings were damaged, which, according to the above-mentioned map, are situated beyond the boundary of the potential landslides. InSAR observations, however, show ongoing activity, especially in the eastern sector. Most unstable parts lie on the Handlova Formation (tuffs and sandstones) and correspond mainly to the landslide deluvia. Several PSs in central part of Hradec village (north-eastern one) even exceed the values of 3 cm/year. Stable northern section of the village, which rests upon a block formed by lava flows of pyroxene andesite, is clearly separated by landslide boundary. Structural failure of buildings is most common in these areas. On the contrary, in the central part of the village, where the greatest velocities are observed, the buildings are "carried" by landslide body without the impact on their construction. Similar, but smoother division is also observable in Velka Lehotka village. Results visualization in terms of displacement velocity maps is for many instances insufficient, as they do not reveal the valuable content in the displacement time series of each individual PS. In order to validate InSAR observations, PS time series were compared with independent ground measurements in borehole inclinometers. Since the measurements are of a different nature, they are displayed in individual graphs. Moreover, considering the geocoding uncertainty, collocation of points is only approximate. Therefore the closest, well coherent PSs are chosen for comparison. Figure 11 shows two examples of such comparison with precise borehole inclinometers codenamed IGH-5i and IGH-1i. The two example boreholes characterize distinct landslide bodies (see Figure 8). Borehole IGH-5i is located in the most active part of the slope deformation with a clearly developed shear zone (Figure 11a). Installation of a stationary inclinometer in the nearby IGP-10i borehole has confirmed that the territory is permanently in motion. Average displacement velocities range from 7 to 68 mm/year, dominantly depending upon climatic factors. Although the shear zone is less evident at borehole IGH-1i (Figure 11c), the presence of landsliding activity is recognizable at nearby buildings ( Figure 11d) and fences, broken by fractures.
While the method of precise inclinometry provides information on the deformation evolution at particular depths of borehole and directly determines the shear zone, the surface associated InSAR observations enable to assess the deformation phenomena across the whole sliding areas. Therefore, those information are complementary in terms of their spatial scope. Comparison with borehole displacement values confirmed the presence, extent and magnitude of slope deformations. Periods of larger movement activity and sudden activation coincide between epochs of inclinometry measurements (vertical colour lines in PS displacement time series). Variations of slope movement activity could be presumably caused by seasonal changes of climatic factors and the related groundwater table level, being one of its major activators [3]. Besides, the magnitude of sub-surface cumulative displacement at shallow depths in IGH-5i borehole adheres to surface ones, as apparent from the time series. Displacement time series tend to gradually stabilize towards the end of an observation period. Finally, note the PS time series variance decrease, stiffened by higher sample rate since October 2016 as a benefit of Sentinel-1B integration. Altogether 30 Sentinel-1B images were employed within processing for ascending (15) and for descending (15) track.
Implementation of InSAR observations within the conventional landslide monitoring techniques is of a major significance, as it helps to better understand the surface deformation pattern throughout the sliding body, its spatial and temporal evolution, as well as to assess the success of remediation works (which have been, in this particular case, deemed insufficient). This underscores the operational capability of Sentinel-1 data for routine updates of landslide inventory maps on a weekly basis, which could even contribute to development of early warning systems.

Case Study 2: Land Subsidence Due to Undermining
Major part of Kos municipality is afflicted by significant loss of coherence due to fast deformation rates and formation of wetlands, as a consequence of land subsidence down below average groundwater level. According to [15], depth of some depressions reaches up to several meters. Theoretical maximum deformation rate observable by InSAR is half a wavelength per revisit time, i.e., 42.6 cm/year for C-band Sentinel-1 [4]. This actually also depends on the noise level and the phase unwrapping technique [37]. Therefore, identifications of PSs in northwestern part of Kos village fails.
Nonetheless, we managed to retrieve several reliable PSs in the southeastern part of village from the descending track, which exhibit coherent deformation as depicted by Figure 12. While the smaller-magnitude subsidence is noticeable even at PSs remote from mining affected areas, displacement values slowly increase towards the directly undermined territory, reaching rates higher than 3 cm/year. Since results from the ascending track are generally noisier and without coherent points over the affected area even after multi-variate outlier elimination (see Figure 7), LOS decomposition was not performed.  Independent displacements time series strongly correlate, especially at benchmark VKB-20, and thus provide the notion about the relative accuracy of InSAR observations. Sole levelling time series of benchmark VBK-14 are also presented in order to illustrate the magnitude of fast deformation rates (exceeding 3 m/year) occurring within the western part of the village, lying directly above the mining corridors. Displacement evolution of some points tends to stabilize (VBK-3, VBK-22), while for example VBK-20 exhibits a non-linear trend with activation since summer 2015, which could have been caused by the restoration of mining activities and the increase of their extent.
Proof of an extreme subsidence occurring in the old westernmost part of Kos village, devoid of PSs, is provided by differences between UAV and LiDAR-derived digital elevation models in Figure 14. The former DEM was derived from Trimble R UX5 (Geotronics Slovakia and Slovak University of Technology) measurement campaign carried out in March 2016 and processed by Agisoft R PhotoScan (Agisoft LLC, St. Petersburg, Russia). Average density of ground filtered point cloud reached 20 points/m 2 . The latter DEM was generated on the basis of Trimble R Harrier 68i LiDAR system (Zilinska University and Slovak University of Technology) measurements which took place in November 2016. Through processing by LAStools (rapidlasso GmbH, Germany), we achieved average density of ground filtered point cloud of 15 points/m 2 . For both cases, DEM pixel size is 50 cm, and vertical accuracy better than 10 cm was confirmed via independent ground control points. DEM comparison reveals a distinct north-western body of Kos village subsiding down 7 m over 8 months period. This gives a clear notion why the identification of coherent PSs using C-band data was impossible for this particular area. The consequences of undermining are evident not only geomorphologically, but also on structures and infrastructure ( Figure 15). Several houses in northwestern and central parts of the village have already been demolished; the remaining, still inhabited houses in central part are continually exposed to risk. Sentinel-1 MT-InSAR proves to be able to secure remote and cost-effective stability monitoring of those structures, comparable to conventional monitoring techniques in terms of accuracy, yet maintaining the vast spatial coverage.
Its limitation concerning the maximum observable deformation rate could be partially overcome by different MT-InSAR techniques, e.g., Small Baseline Subset (SBAS) [38], employing data from longer wavelength sensors [39], or even by exploitation of the amplitude information of SAR images [40].

Conclusions
We have demonstrated the capabilities of open-source MT-InSAR tools SNAP & StaMPS in conjunction with the open-access Sentinel-1 data to monitor the most significant geohazards in the Upper Nitra region, Central Slovakia. Compared to point-wise terrestrial monitoring techniques, the advantages of MT-InSAR combine reasonable spatial coverage, assured temporal sampling and the capability to remotely assess the historic as well as ongoing deformation phenomena at low-cost merit, regarding the Sentinel data. We have investigated currently active slope deformations in Hradec and Velka Lehotka suburbs of the Prievidza city, where major part of villages exhibits coherent, dominantly westward displacements, while annual velocities of some PSs even exceed the values of 3 cm/year even after remediation works. Thanks to the application of multi-variate outlier post-processing analysis, we retained valuable PSs in affected areas, which would have been otherwise (imposing the standard coherence threshold) undersampled. Moreover, we have managed to identify several rapidly subsiding PSs in the dynamically changing environment of the undermined area of Kos municipality. Our observations have been validated via comparison with ground truth data from borehole inclinometers and terrestrial levelling. Outcomes of this analysis have enhanced our understanding of deformation patterns and their spatio-temporal progression. Assuming the detailed analysis of PS kinematic time series along with routine updates of deformation maps, those could contribute to development of early warning systems. Thus the classification of deformation trends based on conditional sequence of statistical tests (as proposed by e.g., [41] or [42]) is among the planned steps for respective case studies. Since both case study deformation phenomena can still be considered substantially active and a number of buildings are subject to continual damage, assessment of remediation works and operational routine monitoring is of utmost importance. Exploiting the free Sentinel-1 data by means of multi-temporal InSAR techniques constitutes to be capable of such task. Therefore, the main future objective involves pilot operational monitoring of at least hereby presented areas, supported by the communication with State Geological Institute of Dionyz Stur, HB mining company and local authorities. Finally, efficient installation of artificial corner reflectors or active radio transponders [43] shall further increase the density of PS at critical sites lacking any strong natural scatterers.