New Perspectives in Landslide Displacement Detection Using Sentinel-1 Datasets

: Space-borne radar interferometry is a fundamental tool to detect and measure a variety of ground surface deformations, either human induced or originated by natural processes. Latest development of radar remote sensing imaging techniques and the increasing number of space missions, speciﬁcally designed for interferometry analyses, led to the development of new and more e ﬀ ective approaches, commonly referred to as Advanced DInSAR (A-DInSAR) or Time Series Radar Interferometry (TS-InSAR). Nevertheless, even if these methods were proved to be suitable for the study of a large majority of ground surface dynamic phenomena, their application to landslides detection is still problematic. One of the main limiting factors is related to the rate of displacement of the unstable slopes: landslides evolving too fast decorrelate the radar signal making the interferometric phase useless. This is the reason why A-DInSAR techniques have been successfully applied exclusively to measure very slow landslides (few centimetres per year). This study demonstrates how the C-band data collected since 2014 by the Sentinel-1 (S1) mission and properly designed interferometric approaches can pull down this restriction allowing to measure rate of displacements ten times higher than previously done, thus providing new perspectives in landslides detection. The analysis was carried out on a test site located in the Cortina d’Ampezzo valley (Eastern Italian Alps), which is a ﬀ ected by several earth ﬂows characterized by di ﬀ erent size and kinematics.


Introduction
The application of imaging radar remote sensing to geosciences is nowadays a fast moving field. Since the launch of Seasat in 1978, the first Earth-orbiting satellite equipped with a Synthetic Aperture Radar (SAR), numerous space missions have followed. The developments of radar remote sensing led to the study of ground surface dynamic phenomena and in 1993, Massonet [1] reported the first successful application of space-borne differential interferometic SAR (DInSAR) to map the co-seismic deformation of the Landers earthquake. In recent times, DInSAR has experienced a continuous growth mainly related to the increased number of missions specifically designed for interferometric analyses and the large amount of acquisitions collected through the years that have favored the generation of a wide number of data processing methods and analysis tools e.g., [2,3]. These techniques, commonly referred to as Advanced DInSAR (A-DInSAR) or Time Series Radar Interferometry (TS-InSAR) are similar implementations of two main approaches: the Persistent Scatterer Interferometry (PSI) [4], based on the identification of natural point-like targets in stacks of interferograms, and the Small Baseline Subsets (SBAS) [5], centered on the privileged spatial and temporal relations of the acquisitions. So far, A-DInSAR analyses have proven to be effective in the detection and measurement of landslide classified as very slow (few centimeters per year) [11][12][13][14][15][16]. This study demonstrates that processing data collected by the S1 mission, with properly designed interferometric approaches, can extend the field of application to the upper class of slow landslides (few metres per year), thus providing new perspectives in landslides detection. The analysis was carried out on the test site of Cortina d'Ampezzo, a town located in a wide valley in the hearth of the Dolomites (Italian Eastern Alps), on whose surrounding slopes more than 30 landslides, differing in type, size, state and stile of activity and rate of displacement have been identified and mapped [17]. The analysis focused on the landslides affecting large portions of the western slope of the Cortina d'Ampezzo valley. Mass movements different in size and rate of displacement by an order of magnitude have been discriminated, mapped and measured, establishing new boundaries for the application of radar remote sensing to landslides hazard assessment.

Study Area
The Dolomites are a mountainous group of the Italian Eastern Alps, numbering 18 peaks above 3000 metres a.s.l. and covering an area of 231,000 hectares. Known for their stunning beauty, featuring spectacular landscapes with one of the best examples of the preservation of Mesozoic carbonate platform systems, the Dolomites were listed as a UNESCO World Heritage Site in 2009. The town of Cortina is located in the Ampezzo Valley, crossed by the Boite River and surrounded by the grandeur of the dolomitic peaks ( Figure 1). The town is certainly one of the most known and fashionable touristic resort of the Dolomites. Cortina hosted the Winter Olympics game in 1956 and, after that a number of international winter sports events. Recently, it has obtained the assignment of the Alpine Ski World Championships, scheduled in 2021 and the 2026 Winter Olympics. The impressive location of Cortina d'Ampezzo, its glamorous legacy as a tourist destination at any time of the year and the natural predisposition at hosting relevant sports events, led to a constant growth of the anthropic pressure and hence the vulnerability in a valley, known to be intensely affected by geo-hydrological hazards.   Since the retreat of the Würmian glaciers, some 11,000 years ago, the Ampezzo valley was largely interested by mass movements [17]. The renewals of these extensive and ancient landslides continuously endanger and in some cases damage houses, roads and skiing facilities. The proneness of the Ampezzo valley to instability phenomena can be attributed to several factors [18][19][20]. First of all, the structural conditions of the area have to be considered. The stratigraphical sequence consists in an alternation of dolomitic rocks and clayey geological formations. The presence of lithotypes with different geomechanical properties represent a predisposing factor wherever the rigid and resistant rocks with a brittle behavior overlie the plastic rocks characterized by a ductile deformation. Another important factor derives from the effects of tectonics. Folds and faults stressed and dismembered the rock masses that became much weaker and prone to erosion. Finally, the glacio-pressure effect [19] has to be considered. It is likely that the pressure of ice on the valley sides during the Würmian glaciation determined rock deformations in correspondence with surfaces of structural discontinuity, favoring the formation of potential sliding surfaces. The glacial debuttressing, permafrost melting and high relief energy have then favored the occurrence of slope instability processes [21][22][23].
The characteristics of the landslides generally depend on the lithological and structural conditions: where the dolomitic rock outcrops, in the highest part of the valley, rock falls occur, while slides and flows are frequent in the medium and lower sectors of the slopes in correspondence of the marly and clayey formations [20,24]. On the western slope of the valley, the earth flow (or mudslide according to the British usage), known as Lacedel landslide, was presumably triggered in the Late Glacial and is characterized by a recurrent activity that continues to date. The mass movement extends for an approximate length of 3.5 km, from the crown located at 1750 m a.s.l. to the toe, in correspondence to the Boite river bed at 1150 m a.s.l. Subsurface investigations estimated the thickness of the landslide between 40 and 60 m, most of which was accumulated in the early postglacial [20]. The landslide body is composed by an alternation of silty clays, originated by the weathering of the claystone bedrock and interdigitated dolomitic gravel lenses deriving from Late Glacial debris flow events [25]. A partial reactivation of the earth flow, named Rio Roncatto, is at present the most active sector of the western valley slope. Damaged and abandoned buildings and disrupted roads witness the power release of this landslide ( Figure 2).

135
(Terrain Observation with Progressive Scans in azimuth) modes to provide large swath width with 136 enhanced image performance [27]. The main characteristics of S1 mission, instrumentation and 137 products are resumed in Table 2.  Recent GNSS (Global Navigation Satellite System) surveys assessed the rate of displacement of the Rio Roncatto landslide body to about 1 m per year, while for the Lacedel landslide the measured deformations are much lower, ranging from in few millimetres per year to several centimetres [26]. However, these measurements are point-based and do not allow to extend their information without incorporating some uncertainty in a bi-dimensional space, specifically a map. Therefore, the definition of the contours of the landslide characterized by different levels of activity, at the time being, must necessary rely also on the qualitative information represented by the geomorphological analysis.

Materials and Methods
S1 satellites are equipped with an active phased array antenna that supports four imaging modes characterized by different resolution and coverage: Interferometric Wide Swath (IW), Extra Wide Swath (EW), StripMap (SM), and Wave (WV). Both the IW and EW mode are implemented as TOPS (Terrain Observation with Progressive Scans in azimuth) modes to provide large swath width with enhanced image performance [27]. The main characteristics of S1 mission, instrumentation and products are resumed in Table 2.
The interferometric analyses were carried out processing 104 S1 images in IW mode along the ascending track n. 44, from the end of October 2014 until late November 2018. Wintertime acquisitions (from December to March) were discarded due to the likely presence of snow that would have increased signal noise. The images were multi-looked by a factor 10 in range and 2 in azimuth direction in order to reduce the speckle noise. This procedure increased the pixel spacing to 23.2 m and to 27.8 m respectively. The datasets were co-registered to a single image acquired on the 17 th of July 2017. Even if the S1 orbital parameters are very reliable and despite the fact that the deviations are maintained within a tube of +/-50 m radius (RMS) [27], precision state vectors, downloaded from the European Space Agency (ESA) Quality Control website (https://qc.sentinel1.eo.esa.int/aux_poeorb/), were used to guarantee highly accurate results. IW images are composed by 3 swaths made of bursts which are affected by large Doppler centroid variations as a consequence of the implementation of the TOPS acquisition mode. The resulting steep azimuth spectrum ramp within each burst requires misregistration errors smaller than 5/1000 of sample in order to avoid residual interferometric phase jumps [28]. As a consequence, a stringent co-registration procedure is obliged. To achieve this task, an iterative procedure was used. A lookup table that permits to resample the data among the master and the slave geometries was derived by exploiting the Shuttle Radar Topography Mission (SRTM) DEM sampled at 3 arc-seconds [29]. After the resampling, the offsets among the images were determined through a matching procedure that uses iteratively cross correlation optimization to reduce mis-registration errors to about 1/100 of sample. Afterwards, a spectral diversity method, which considers the double difference phase in the bursts overlap regions, was applied. The double differences interferograms were multi-looked and unwrapped. The algorithm contemplates only overlap regions with a minimal coverage of high coherence pixels (e.g., 0.8) and the statistics of the unwrapped phases to avoid errors. The accepted values were then combined using a weighted averaging. The offset parameter file is updated, the slave images are resampled and the spectral diversity method applied once again. The procedure was iterated until the mis-registration errors result smaller than 5/1000 of sample. Afterwards, for each burst, the azimuth phase ramp of the master acquisition was calculated, subtracted and then used to deramp the co-registered slave images. The stack of the 104 co-registered and deramped images was then resized to a square of 600 x 600 pixels, covering area of interest for an extension of more than 230 km 2 . At this stage, two different interferometric procedures were implemented in order to detect both the Lacedel and Rio Roncatto landslides, whose displacement rates differ considerably. The investigation of the former was carried on using an A-DInSAR approach known as Interferometric Point Target Analysis (IPTA) [30] while, for the latter, standard differential interferometric approaches were applied.

Interferometric Point Target Analysis
IPTA technique exploits the temporal and spatial characteristics of interferometric signatures back-scattered by natural point-like targets to assess their displacements histories [30]. Point target candidates were selected using two approaches. The first one is based on the spectral properties of prevalent targets within each individual image. Spectral characteristics were averaged over the set of images and targets that show low spectral phase diversity were selected as point target candidates (PTC). The second one is based on the very likely assumption that dominant targets have a lower temporal variability throughout the acquisitions stack compared to distributed scatterers. PTC were hence selected by fixing a threshold to the ratio between the temporal average of backscattering and its standard deviation. Every image in the stack was processed along with the master to generate 103 interferograms with an average perpendicular baseline of 23 m and a maximum temporal baseline of 858 days (Figure 3).   Afterwards, the interferometric phase at each PTC position was flattened by subtraction with the unwrapped synthetic phase derived from the SRTM DEM. The resulting differential interferometric phases between pairs of PTC were analyzed by means of a two-dimensional regression that accounts for the perpendicular baseline and the time elapsed between the interferometric pairs. The parameters resulting from the regression such as the height corrections, the deformation rate, the baseline refinement, the atmospheric phase term and the extension of the PT (Point Target) list, were iteratively improved until the phase model resulted satisfactory.

DInSAR Analysis
To detect and isolate fringes compatible with the surface displacements of the Rio Roncatto earth flow, a multi-baseline approach was applied. By setting a maximum temporal baseline equal or minor to 24 days, 508 differential interferograms were generated. In several highly coherent phase images a deformation pattern attributable to the mass movement was clearly visible even if the noise and the atmospheric phase screen contribution made it often difficult to discriminate the signal (Figure 4). In order to map the Rio Roncatto landslide and measure its rate of movement, an interferometric approach based on two very likely assumptions was designed. Firstly, we assumed the stability of the area surrounding the Rio Roncatto landslide and specifically the Cortina d'Ampezzo urban area and the remaining portion of the western slope (by inspection of the IPTA results, considering a maximum temporal baseline of 24 days, the rate of deformation of the latter can be considered negligible). Under this assumption, we masked the Rio Roncatto landslide and applied a spatial filtering and an interpolation of the interferometric phases over the "stable" areas. The result was interpreted as atmospheric disturbances and subtracted from the differential interferogram. The filtered and flattened complex valued interferograms were unwrapped using a minimum cost flow algorithm [31] and visually inspected to assess the presence of phase errors. The second assumption was that the rate of deformation of the earth flow can be considered constant as demonstrated by the GNSS (Global Navigation Satellite System) surveys operated by Bossi et al. (2016) [26]. Between July 2008 and April 2012, eight fast static relative measurements were carried out. The inferred rate of displacements for the four benchmarks installed over the Rio Roncatto landslide were assessed between 1 and 2 m per year ( Figure 5).
Due to temporal decorrelation, interferograms generated from a short time intervals image pairs are generally more coherent. In other words, it is likely that the phase unwrapping would results more reliable and correct for 6-day interval interferograms rather than 12, 18 or 24. Based on this assumption, we added the unwrapped phases of the 6-day interferograms to create phase models for 12 or more days. The complex valued interferograms generated from images with a temporal baseline longer than 6 days were unwrapped using an algorithm forcing the phase of each pixel to the value in the (-π, π) interval of the model. After completion of the unwrapping, a phase constant was subtracted from all unwrapped phase values such that the deformation phase at the indicated stable reference location (the same used in the IPTA processing) would equal 0. The linear rate of differential phase backscattered was estimated from the Rio Roncatto landslide from a stack of six unwrapped differential interferograms. The reference image, acquired on September 22, 2017, was processed along with the 6 subsequent images taken 6 days apart from each other. Phase rate estimation was obtained by weighting the individual interferogram phases for the time interval under the assumption that atmospheric statistics were stationary for the stack of interferograms, according to the following equations: Where: phrate is the phase linear rate, var(phrate) is the variance of phase linear rate, N is the number of interferograms (in this case 6), t is the time interval, φ is the differential interferometric phase, λ is the signal wavelength.
In order to map the Rio Roncatto landslide and measure its rate of movement, an interferometric 207 approach based on two very likely assumptions was designed. Firstly, we assumed the stability of

Results
The outcomes of the IPTA processing are reported in Figure 6. Patterns of ground displacements driven by the Lacedel landslide on the western slope of the Cortina d'Ampezzo valley are clearly visible in correspondence of the PT located in the hamlets of Colfiere, Lacedel, Mortisa and Grignes. The town of Cortina d'Ampezzo, at the bottom of the valley, results to be stable, as expected. The majority of the PT are detected in correspondence of built-up areas, ski lifts pylons and a few rock blocks in the upper part of the slope. No displacement information can be inferred over the body of the Rio Roncatto landslide. The maximum displacement rate calculated at PT location varies from 10 (Colfiere) to 6 mm/yr (Lacedel).    A quantitative assessment of the IPTA analysis outcomes can be inferred from the histogram in Figure 7. PT have been grouped in displacement rate classes with a bin-width of 1 mm/yr. Negative values of the velocity witness an increasing distance between the sensor and the ground target throughout time. It can be noticed that the data arrange as a normal distribution slightly skewed to the left, towards negative values. The large majority of the PT are grouped around the null value of the rate of displacement and reflects the kinematics behavior of the radar target located in the town of Cortina d'Ampezzo. Negative values are related to the mass movements affecting the western slope of the valley (i.e. Lacedel landslide). On the contrary, the majority of the PT spread on the right side of the histogram can be attributable to the eastern slope instabilities. The positive values of the rate of deformations are recoded by targets moving toward the sensor and are hence suitable considering the ascending geometry of acquisition of the processed dataset and the prevalent direction of the displacements.
The assessment of the accuracy of the IPTA analysis outcomes can be inferred from the histogram in Figure 8. PT have been grouped in estimated displacement rate uncertainty classes with a bin-width of 0.05 mm/yr. Except very few outliers with values larger than 0.5 mm/yr, most of the targets are grouped in the classes characterized by a displacement rate accuracy ranging from 0.25 to 0.45 mm/yr.
The results of the DInSAR analysis are shown in Figure 9. The pattern of deformation delineates quite well the shape of the Rio Roncatto landslide body mapped from geomorphological evidences. The rate of displacement of the earth flow, assessed along the satellite line of sight (LOS), reaches a peak of 0.5 m/yr, more than an order of magnitude with respect to the values measured with the IPTA technique. The most active sectors of the earth flow are located at the source area of the landslide and to the southern sector of the toe. The sector developing north of Grignes recorded slower rate of displacement rate. Deformation patterns are also clearly evident outside the mapped landslide, specifically north of the source area and uphill with respect to the hamlet of Lacedel (white dotted contours in Figure 9). These findings might suggest an updating of the landslide inventory maps.       Figure 10 sums up the results of IPTA processing performed with S1 and ENVISAT datasets, restricted to the unstable area as mapped by geomorphological evidence.

Discussion
The most obvious remark is that the number of PT in the S1 dataset is much larger than the one in the ENVISAT. The PT located inside the landslide contours are 512 and 101 respectively, meaning that the number of natural reflectors suitable for interferometric analysis are five times larger in the S1 dataset. S1 acquisition mode and the short revisiting time favors the identification of PT leading to an improved spatial sampling and the possibility to infer a greater amount of displacement data in the area of interest. The cluster of PT located over the hamlet of Colfiere in the S1 dataset and absent in the ENVISAT analysis clearly witness the upgrade of the effectiveness of the latest ESA space mission for interferometric purposes. A second consideration can be drawn on the consistency of the S1 dataset. Patterns of deformation are detected in the same sectors of the unstable slope as defined by the ENVISAT data, although with rates slightly different that can be attributable to a change in the style of activity of the landslide and marginally to the different acquisition geometries between S1 and ENVISAT. A comparison between rate of displacement of quasi-homologous PT is shown the graphs reported in Figure 11.

309
The most obvious remark is that the number of PT in the S1 dataset is much larger than the one  Beside the interesting results of the IPTA processing the outcomes of the DInSAR analysis draw particular attention. If the capability to detect and accurately map the earth flow is clearly highlighted in Figure 8, the reliability of the estimation of the displacement rate needs to be further investigated. In order to achieve this purpose, the cumulated displacement recorded at GNSS benchmarks locations, between the first reading in 2008 and the last one in 2012, has been decomposed and projected along the S1 sensor LOS according to the following formula: of the S1 dataset. Patterns of deformation are detected in the same sectors of the unstable slope as 318 defined by the ENVISAT data, although with rates slightly different that can be attributable to a 319 change in the style of activity of the landslide and marginally to the different acquisition geometries 320 between S1 and ENVISAT. A comparison between rate of displacement of quasi-homologous PT is 321 shown the graphs reported in Figure 11.

325
Beside the interesting results of the IPTA processing the outcomes of the DInSAR analysis draw 326 particular attention. If the capability to detect and accurately map the earth flow is clearly highlighted 327 in Figure 8, the reliability of the estimation of the displacement rate needs to be further investigated.

328
In order to achieve this purpose, the cumulated displacement recorded at GNSS benchmarks 329 locations, between the first reading in 2008 and the last one in 2012, has been decomposed and 330 projected along the S1 sensor LOS according to the following formula:

339
Broadly assuming a linear deformation ( Figure 5), the rate of displacement along the LOS can Broadly assuming a linear deformation ( Figure 5), the rate of displacement along the LOS can simply be calculated by the following equation: Where LOS rate is the rate of deformation along the LOS, t is the time elapsed between the first and last GNSS measurement expressed in days. Table 3 summarizes the input data used in the equation and the results calculated for the 4 GNSS benchmarks localized in the lower sector of the Rio Roncatto landslide. Each benchmark location has then been transformed in range-Doppler coordinates (i.e. range and azimuth). The rate of displacement calculated through the interferometric phase of the corresponding image pixel is reported in Table 4.
The results presented in Table 3 are in good agreement with the values inferred by the DInSAR analysis (last column in Table 4), showing that the approach, beside an excellent capability and accuracy of mapping and detection, provides also a reliable estimation of the rate of displacement. By comparing the LOS velocities calculated from the GNSS measurements (GNSS-LOS rate) with the values inferred from the interferometric analysis (LOS rate), a good agreement can be noticed. The maximum difference of 0.155 m/yr is found at ID 20 GNSS benchmark location. It has however to be noticed that, beside the assumption of a constant rate of displacement of the Rio Roncatto earth flow, GNSS-LOS rate were inferred from a 4-year measurement (2008-2012) while LOS rate is calculated on a monthly basis (36 days in 2017), hence discrepancies between the values were expected.

Conclusions
The impact of climate change, which results in increasing frequency of extreme meteorological events, makes landslides a worldwide problem of great concern, especially considering the growth of population that leads to the construction of new settlements, the enlargement of cities and the increasing need for infrastructures [32]. Constructions in hazardous areas might lead to unsafe situations and project uncertainties, causing large economic consequences and sometimes casualties. Aiming at an effective risk management and mitigation, the crucial tasks to be achieved are [7]: • Detection: where and when deformation occurred (or is occurring); • Assessment of the deformation; • Evolution (past and future) of the deformation.
Monitoring techniques have always been considered one of the most common tools for risk mitigation and, among these, space-borne interferometry plays a key role. A-DInSAR techniques are currently the only cost-effective approaches for frequent, large-scale monitoring of ground deformation, able to provide the information needed to accomplish the aforementioned objectives. Still these approaches might fail when the inherent characteristics of the mass movements interfere with physical constraints of the technique. Typically, the size, the rate and direction of the mass movements represent the limiting factors to investigate landslides by means of A-DInSAR. Recent space missions, specifically designed for interferometric purposes, such as the ESA S1 satellites, are providing the scientific community with new effective datasets. S1 mission has certainly increased the potential radar interferometry techniques to detect landslides, mainly for the short repeat cycle and the small perpendicular baselines among the acquisitions, besides the timeliness and the reliability of a free of charge service.
Monitoring very slow landslides by means of A-DInSAR analysis is nowadays a well-consolidated practice, but this study has demonstrated how properly designed approaches can extend the field of application of space-borne radar interferometry to the upper class of slow landslides (few metres per year). The landslides affecting the western slope of the Ampezzo valley have been correctly and accurately characterized beside their wide differences and terms of spatial extension and kinematics.
Knowing that the assessment of landslide risk includes both the area affected and the velocity and the product of these two parameters is approximately proportional to the power release of the landslide, the results obtained in the Rio Roncatto earth flow investigations are relevant. The capability to accurately map landslides contours and to carry out reliable displacement measurements allow to produce a data-driven distributed activity map, providing high quality information for geotechnical modelling, which are fundamental to predict risks scenarios and to design mitigation works. It is worth noticing that, prior to the launch of S1, it would have been unfeasible to characterize an earth flow displacing at a velocity in the order of decimetres to metres per year in such a continuous and extensive way. Ultimately, S1 mission has proved to provide a fundamental contribution in the challenging field of landslide risk mitigation and management Author Contributions: M.M. conceived the paper, performed the interferometric analyses and interpreted the results. G.B., L.S. and G.T. organized and analyzed the data; G.M. and G.Te. supervised the geomorphological field investigations; A.P. led the research. M.M., G.B. and G.Ti. wrote the paper. All the authors contributed to editing and revisions of the paper.
Funding: This research received no external funding.