Sentinel-1 and Ground-Based Sensors for Continuous Monitoring of the Corvara Landslide ( South Tyrol , Italy )

The Copernicus Sentinel-1 mission provides synthetic aperture radar (SAR) acquisitions over large areas with high temporal and spatial resolution. This new generation of satellites providing open-data products has enhanced the capabilities for continuously studying Earth surface changes. Over the past two decades, several studies have demonstrated the potential of differential synthetic aperture radar interferometry (DInSAR) for detecting and quantifying land surface deformation. DInSAR limitations and challenges are linked to the SAR properties and the field conditions (especially in mountainous environments) leading to spatial and temporal decorrelation of the SAR signal. High temporal decorrelation can be caused by changes in vegetation (particularly in nonurban areas), atmospheric conditions, or high ground surface velocity. In this study, the kinematics of the complex and vegetated Corvara landslide, situated in Val Badia (South Tyrol, Italy), are monitored by a network of three permanent and 13 monthly measured benchmark points measured with the differential global navigation satellite system (DGNSS) technique. The slope displacement rates are found to be highly unsteady and reach several meters a year. This paper focuses firstly on evaluating the performance of DInSAR changing unwrapping and coherence parameters with Sentinel-1 imagery, and secondly, on applying DInSAR with DGNSS measurements to monitor an active and complex landslide. To this end, 41 particular SAR images, coherence thresholds, and 2D and 3D unwrapping processes give various results in terms of reliability and accuracy, supporting the understanding of the landslide velocity field. Evolutions of phase changes are analysed according to the coherence, the changing field conditions, and the monitored ground-based displacements.


Introduction
Natural disasters induced by gravitational mass movements are widespread phenomena of various magnitude caused by geological and climatic conditions or induced by anthropogenic factors [1].Slope displacement can be detected and monitored by several Earth observation techniques [2], while the selection of an adequate monitoring concept depends on the scope of the study, the type and scale of phenomena, the data available, and the skills of the investigators [3,4].The kinematics and spatial and temporal evolution of landslides using InSAR-based techniques have been analysed in a large number of scientific studies [5][6][7][8].
New generations of satellites, such as the Copernicus Sentinel-1 (S1), open up new perspectives for continuous ground surface monitoring, while being characterised by enhancements in terms of revisit time, coverage, timeliness, and reliability of service [9].Indeed, the C-band S1 synthetic aperture radar (SAR) instrument is specifically designed to carry out interferometric analyses over land [10].Recent studies have shown that S1 data allows analysing Earth's surface displacements using differential SAR interferometry (DInSAR) techniques [11][12][13].DInSAR has the capability to precisely monitor surface displacements over time (temporal sampling rate up to 6 days for Sentinel-1A/B) with a wide coverage in a labour-saving, time-, and cost-efficient manner [14].Two main categories of advanced DInSAR processing techniques for displacement time-series generation exist: persistent scatterer interferometry (PSI) and small baseline subsets (SBAS) techniques [15][16][17].Both approaches can be also merged together [18].Whilst these methods have shown great potential for landslide monitoring and detection [19], some limitations remain, especially in the detection and estimation of significant changes in vegetated terrains [20].
According to [21], a common challenge faced in both PSI and SBAS algorithms is the phase unwrapping (PhU) operation that represents the retrieval process of the absolute phase signals from their (measured) modulo-2π restricted components, that is, interferometric fringes.The unwrapping operator estimates an ambiguity that is an integer number of 2π radians in wrapped data.Several approaches have been used to remove this ambiguity and minimise the total length of discontinuities in the unwrapped phase: Branch cut [22] relies on the restriction of integration over images to paths with local phase differences within 2π.Minimum cost flow (MCF) is a statistically based approach to optimize flow in each of the arcs that minimizes the total cost [23].Least-square-based methods [24] are based on a two-dimensional partial differential equation.
Temporal and spatial decorrelation, as well as atmospheric and ionospheric noises, are often present in DInSAR-based results [25,26].However, the area-wide coverage of SAR data can compensate the low spatial sampling of ground-based data, while the high temporal sampling rates of ground-based data overcome the low temporal sampling rates of SAR data.Accordingly, DGNSS and DInSAR could be utilised as complementary methods to fully quantify surface displacement and not only for validation purposes.DGNSS data, for instance, can reach a precision in the order of millimetres at high temporal sampling rates (depending on the configuration mode) on one single point.The combination of DGNSS and DInSAR estimated displacements [27][28][29][30], derived by using well-known interpolation techniques such as minimum curvature [31], Kriging [32], or least-squares collocation [33], can provide an area-wide displacement map and a 3D displacement map [34].DGNSS data could also assist in mitigating the atmospheric artifacts and correcting water vapor effects from the SAR data [35,36].
PSI using COSMO-SkyMed data has only been applied successfully to the Corvara landslide to the X-band corner reflectors (CRs) installed there, while SBAS failed to produce a spatial displacement map due to the high sensitivity of X-band data to vegetation [37].The recent availability of C-band data such as S1 data, also characterised by frequent and constant acquisitions, gives a promising opportunity to highlight the feasibility and accuracy of monitoring a complex and vegetated continuously moving landslide by applying the SBAS algorithm to S1 data.This paper has two main aims: (1) assessing the performance of 2D and 3D phase unwrapping methods in a vegetated landslide and (2) evaluating the application of DInSAR using S1 data and DGNSS to continuously monitor the Corvara landslide.To this aim, the ground-based and remote sensing data are presented as well as the methodological workflow in the following sections.Afterwards, DInSAR results are analysed and compared with DGNSS data in order to discuss the potential and limitations of the new S1 mission for the continuous spatial monitoring of complex and vegetated landslides.

Study Site
The Corvara landslide is located in the Eastern Italian Alps, in the Autonomous Province of Bolzano (South Tyrol).The investigated site is located in the Dolomites, a mountain range that consists mainly of sedimentary rocks of Triassic origin, whose peaks can be up to 3000 m high.The nearly vertical dolomitic mountains are surrounded by a gentle, hilly terrain that consists of mainly clay-rich rocks [38].These weak rock masses are the setting of many unstable slopes, among them also the Corvara landslide.Radiocarbon dating has shown that the most ancient signs of activity occurred around 10,000 cal BP (calibrated years before the present) in the present day's source area, where large rotational slides of 50-100 m depth were triggered.A second phase of movement occurred probably between 6500 and 2300 cal BP [38][39][40].The landslide, classified as a deep-seated rotational earth flow, is moving a volume of about 30 million m 3 , affecting a surface area of 2.5 km 2 [41].The movement velocity varies from slow to very slow according to the classification of Cruden and Varnes [42].The slide can be subdivided into source, track, and accumulation zones [41].The source zone covers the eastern, upper part of the landslide and shows an inhomogeneous pattern in terms of movement rates.Mobilised material is transported along the narrow track zone towards the accumulation zone of the landslide in the immediate proximity of the village of Corvara.
Borehole surveys by [41] have shown that the depth of the present-day main active shear surface varies between 10 m and 48 m.Groundwater conditions were found to be very complex, probably due to overlapping confined aquifers in the landslide body.Piezometric surveys found groundwater tables ranging from being in proximity to the surface to 40 m depth [41].The Corvara landslide has also been monitored by radar interferometry techniques such as DInSAR and subpixel offset tracking methods using X-band corner reflectors [37,43,44].

Data
The remote sensing data consists of 41 S1A/B images (i.e., including only single bursts no. 9 (S1A) and 7 (S1B) of the swath no. 3) covering the period from January 2015 to October 2016 (Table 1).S1 imagery were processed by the SBAS algorithm [16], as described in the next section.A DGNSS reference station (Ciampai), operated by the South Tyrolean Positioning Service (STPOS), delivers the data to correct the measured raw DGNSS data to an accuracy in the millimetre range.Raw data measured during monthly field campaigns is processed with the Leica © Geo Office software.The processing of permanent DGNSS data has recently been operationalised with the Leica © Spider software.The current DGNSS monitoring network was established in 2013 and consists of three permanent DGNSS stations and 13 periodically measured benchmarks (landslide movement destroyed the stations 51 and 55), which are equipped with a Leica GM10 receiver and an AS10 antenna.The 13 benchmarks are surveyed during field visits with DGNSS devices (a Leica © GS10 receiver and AS10 antenna) once a month with a measurement time of at least 1.5-2 hours, allowing a precision of 1 mm (1 standard deviation).The DGNSS monitoring network (Figure 1A) collects pointwise information about the landslide displacement and can also be used for validation purposes of the results obtained by DInSAR processing.At each DGNSS-surveyed point, both permanently and monthly, a corner reflector with a support system for the DGNSS antenna is installed (Figure 1B).The CRs are of different sizes in order to be captured by SAR sensors of different wavelengths.X-band CRs were installed in 2013, and ground surface displacements were measured by multitemporal imagery coming from the Cosmo Sky-Med satellite [37].

Methods
The selected parameters and processing steps carried out with the software SARscape [45] to estimate the quantitative landslide movement are summarized in Figure 2. First, the diagram of connections to selected S1 image pairs was created before the generation of the interferograms.It defines the combination of SAR pairs to be processed choosing the appropriate master and slave images, after calculating their normal and temporal baselines.

Methods
The selected parameters and processing steps carried out with the software SARscape [45] to estimate the quantitative landslide movement are summarized in Figure 2. First, the diagram of connections to selected S1 image pairs was created before the generation of the interferograms.It defines the combination of SAR pairs to be processed choosing the appropriate master and slave images, after calculating their normal and temporal baselines.In the interferometric processing section, coregistration, interferogram generation, multilooking, filtering, and phase unwrapping were performed successively.The slave images were coregistered with the master using a 2.5-m resolution digital elevation model (DEM) provided by the Autonomous Province of Bolzano.In the coregistration step, (1) a local nonparametric shift estimate is calculated by using DEM and orbital information, (2) a set of windows (64 × 64 in range and 128 × 128 in azimuth) was defined on the master image to compute the cross-correlation function, and (3) the residual parametric shift was estimated and the proper shifts were applied.To increase the signalto-noise ratio (SNR) of the interferograms (more reliable coherence estimation) and obtain squared pixels, a multi-looking factor of 4 × 1, leading to the pixel size of 13.27 m × 13.8 m, was used.To smoothen the interferograms, the Goldstein filter was applied before the final interferograms' unwrapping.Two coherence thresholds (CC) of 0.2 and 0.35 were used for phase unwrapping (PhU).The 2D PhU is based on the MCF algorithm with a regular grid covering the data extent [46], while the Delaunay MCF (hereafter named 3D PhU) utilizes a Delaunay triangular network, only considering pixels with coherence values above the unwrapping (CC) set in the unwrapping step.To this end, 3D PhU is performed through two steps: (1) a "temporal" unwrapping operation based on the MCF approach for each arc, connecting neighboring pixels on the perpendicular baseline and time grid ( ⊥ × T), and (2) utilizing the results of the first step to operate a "spatial" unwrapping on each single interferogram on the range and azimuth grid (R×A) [47].In the 3D PhU, each pixel is unwrapped if its coherence value is above the CC in at least 75% of the pairs in the Delaunay connections, leading to an increase in the robustness and reliability of the unwrapping process.In order to reduce unwrapping error due to the low coherence area, a decomposition level (DL) In the interferometric processing section, coregistration, interferogram generation, multi-looking, filtering, and phase unwrapping were performed successively.The slave images were coregistered with the master using a 2.5-m resolution digital elevation model (DEM) provided by the Autonomous Province of Bolzano.In the coregistration step, (1) a local nonparametric shift estimate is calculated by using DEM and orbital information, (2) a set of windows (64 × 64 in range and 128 × 128 in azimuth) was defined on the master image to compute the cross-correlation function, and (3) the residual parametric shift was estimated and the proper shifts were applied.To increase the signal-to-noise ratio (SNR) of the interferograms (more reliable coherence estimation) and obtain squared pixels, a multi-looking factor of 4 × 1, leading to the pixel size of 13.27 m × 13.8 m, was used.To smoothen the interferograms, the Goldstein filter was applied before the final interferograms' unwrapping.Two coherence thresholds (CC) of 0.2 and 0.35 were used for phase unwrapping (PhU).The 2D PhU is based on the MCF algorithm with a regular grid covering the data extent [46], while the Delaunay MCF (hereafter named 3D PhU) utilizes a Delaunay triangular network, only considering pixels with coherence values above the unwrapping (CC) set in the unwrapping step.To this end, 3D PhU is performed through two steps: (1) a "temporal" unwrapping operation based on the MCF approach for each arc, connecting neighboring pixels on the perpendicular baseline and time grid (B ⊥ × T), and (2) utilizing the results of the first step to operate a "spatial" unwrapping on each single interferogram on the range and azimuth grid (R×A) [47].In the 3D PhU, each pixel is unwrapped if its coherence value is above the CC in at least 75% of the pairs in the Delaunay connections, leading to an increase in the robustness and reliability of the unwrapping process.In order to reduce unwrapping error due to the low coherence area, a decomposition level (DL) processing was exploited.The DL operates data multi-looking and undersampling in an iterative way to unwrap the interferograms at a lower resolution and then reconstruct them back at the original resolution [45].A third-degree polynomial was applied to the selected 30 reference points to estimate and remove the phase constant and ramp during the refinement and reflattening steps.The reference points were selected outside the actively moving landslide area on highly coherent pixels (i.e., more than 0.8), where unwrapped values showed no residuals and jumps.The selected points were mainly located in the permanent settlement area of Corvara and on the surrounding mountains of the landslide, where no displacement is expected.Then, 15 out of the 30 points (including 12 points over the urban area) were chosen as stable points.Finally, the pixel values were interpolated using a 4th cubic convolution interpolation, considering 16 surrounding pixels with interpolation and mean window sizes of 7 and 3, respectively, to derive an area-wide spatially comprehensive displacement map.
Two steps of inversion were performed in order to estimate the landslide velocity field.In the first inversion step, the SBAS inversion kernel was implemented to retrieve the first estimate of the displacement rate and residual topography (∆φ topo ) using a linear model [16].In the second inversion step, after the initial estimation of the displacement time series, a custom atmospheric filtering was defined and applied to the results of the previous step in order to estimate and remove the atmospheric contribution (∆φ atm ). with The filter window sizes of 1200 m and 365 (days) were set to take into account spatial and temporal atmospheric variations, respectively.Finally, the displacement time series and mean displacement map obtained were geocoded according to the UTM reference system.
Finally, the spatial displacement patterns were analysed and compared to DGNSS data from field observations.Based on the DGNSS time series, 3D displacements were computed and projected in the 1D line of sight (LOS) of the satellite, enabling direct comparison with DInSAR results.The decomposition of the LOS displacement vector described in the East-North-Up reference system is [0.64, 0.17, 0.74], according to a descending S1 incidence angle θ inc = 42.2 • and a satellite path azimuth α = −165 • (counter-clockwise to the North) [48].As the decomposition LOS vector shows, the InSAR system is more sensitive to the Up component (i.e., 0.74) rather than East one (i.e., 0.64) and presents the least sensitivity to the North component, whereas DGNSS system indicates more sensitivity to the East-North plane rather than the Up component.

Sentinel-1 DInSAR Analysis
S1 acquisitions were analysed according to meteorological data and perpendicular baseline values for reducing the temporal decorrelation of the DInSAR results.Area-wide snow cover was observed before 10 January 2015 and between 15 January and 15 May 2016 in the meteorological time series.
Figure 3A-D presents the S1 images at their particular acquisition date and their connections between each other to create interferometric pairs.The selected acquisitions are displayed as slave datasets (green points) and the master (yellow point).Due to the difficulty in retrieving phase information from pixels with coherence values lower than 0.2, the pairs with a mean coherence lower than 0.2 (calculated within the landslide boundary) were removed from the graph presented in Figure 3A.Six acquisitions were identified as snowy data, including 28 January 2015, 4-28 February 2016, 11-23 March 2016, and 4 April 2016, and the final connection graph was obtained, as shown in Figure 3C.After applying the mean coherence threshold of 0.2 for discarding the incoherent pairs and removing the snowy acquisitions (end of 2015 and spring 2016), we observed a few connections at this part of the network (red points in Figure 4C), which could potentially lead to phase ambiguity in the phase unwrapping step.To overcome this problem as best as possible with the data available, we needed to have more connections with respect to the other part of the network due to the longer temporal baseline caused by removing the snowy data.To this end, we decided in this time interval to lose the constraints (i.e., 0.2 threshold) in order to preserve more connections.Therefore, the connections with a mean coherence value lower than 0.2 (i.e., 0.18 up to 0.2 within the landslide boundary) were used (Figure 3C).Delaunay connections for 3D phase unwrapping are shown in Figure 3B. the phase unwrapping step.To overcome this problem as best as possible with the data available, we needed to have more connections with respect to the other part of the network due to the longer temporal baseline caused by removing the snowy data.To this end, we decided in this time interval to lose the constraints (i.e., 0.2 threshold) in order to preserve more connections.Therefore, the connections with a mean coherence value lower than 0.2 (i.e., 0.18 up to 0.2 within the landslide boundary) were used (Figure 3C).Delaunay connections for 3D phase unwrapping are shown in Figure 3B.   Figure 5A,B presents the coherence (with a threshold of 0.2) and displacement maps along with the locations of 8 DGNSS stations and 30 reference points.The following convention for presenting displacement values is used in Figure 5B: (i) the negative values refer to an increase of sensor-toground distance and (ii) the positive values indicate a decrease of sensor-to-ground distance.Station 13, 28, 53, 54, and 56 were excluded from SBAS processing because their displacement rate was too high to be quantified using the presented data and SBAS technique.The results show that the spatial dispersion of phase and coherence values are mainly located on the toe of the landslide and urban areas, where the absence of data is represented by the black color (Figure 5A,B).The positive effect of using a shorter temporal baseline (6 days obtained by adding the Sentinel-1B data) in increasing the interferogram quality and coherence is observed in Figure 4A,B.The low-coherence values observed in the Corvara area shows the negative effects of decorrelation caused by vegetation (Figure 4C,D) and snow (Figure 4E,F).
Figure 5A,B presents the coherence (with a threshold of 0.2) and displacement maps along with the locations of 8 DGNSS stations and 30 reference points.The following convention for presenting displacement values is used in Figure 5B: (i) the negative values refer to an increase of sensor-to-ground distance and (ii) the positive values indicate a decrease of sensor-to-ground distance.Station 13, 28, 53, 54, and 56 were excluded from SBAS processing because their displacement rate was too high to be quantified using the presented data and SBAS technique.The results show that the spatial dispersion of phase and coherence values are mainly located on the toe of the landslide and urban areas, where the absence of data is represented by the black color (Figure 5A,B).Figure 6 shows the differences of mean LOS velocities obtained with four combinations of CC and PhU techniques.Spatial displacement patterns are similar at the middle part of the landslide (blue colour) and differ at the depletion areas (eastern part of the landslide in Figure 6), depending on the chosen processing parameters.According to the negative velocity difference values, the propagation of the displacement flow from east to west is clearly visible in Figure 6A-D.Whilst a higher velocity peak (i.e., around 120 mm/year) is obtained with the lower CC and 3D PhU processing, warm colours are visible on the northern depletion area and very little displacement is observed in the eastern depletion area (Figure 6D).However, the frequency of negative values (mainly considering western direction displacement) is more pronounced considering a higher CC and 2D PhU (Figure 6A,C).It seems that "more successful" is due to "frequency of negative values", while in reality, the motivation is the agreement with DGNSS measurements: "frequency of negative values" is a consequence of that point.The a-b cross-section comparison of interferograms confirms that using 2D PhU (followed tightly by 3D PhU) induced a higher frequency of negative values, especially in the eastern depletion area, while considering a higher CC (green line in Figure 7).Figure 6 shows the differences of mean LOS velocities obtained with four combinations of CC and PhU techniques.Spatial displacement patterns are similar at the middle part of the landslide (blue colour) and differ at the depletion areas (eastern part of the landslide in Figure 6), depending on the chosen processing parameters.According to the negative velocity difference values, the propagation of the displacement flow from east to west is clearly visible in Figure 6A-D.Whilst a higher velocity peak (i.e., around 120 mm/year) is obtained with the lower CC and 3D PhU processing, warm colours are visible on the northern depletion area and very little displacement is observed in the eastern depletion area (Figure 6D).However, the frequency of negative values (mainly considering western direction displacement) is more pronounced considering a higher CC and 2D PhU (Figure 6A,C).It seems that "more successful" is due to "frequency of negative values", while in reality, the motivation is the agreement with DGNSS measurements: "frequency of negative values" is a consequence of that point.The a-b cross-section comparison of interferograms confirms that using 2D PhU (followed tightly by 3D PhU) induced a higher frequency of negative values, especially in the eastern depletion area, while considering a higher CC (green line in Figure 7).the DGNSS benchmarks on the Corvara landslide.Station number 8 had been manipulated several times between October 2015 and November 2016 and therefore it had to be excluded.Time series of measured 3D displacements at 13 stations from the source, track, and accumulation zones are plotted in two different charts, according to the velocity classes of Cruden and Varnes [42] (Figure 8).The 3D displacements measured in 2015 and 2016 range from 176 mm (at station 6 in the accumulation zone of the landslide) to 4956 mm (at station 56 in the source zone of the landslide).For the DGNSS stations in the source zone (23,28,49,56, 57, and 58), a deceleration during winter and an acceleration during spring and summer months can be observed (Figure 9).The displacement follows the rainfall rates, which could be explained by the fact that at the upper zone of the landslide, the authors of [41] found the groundwater table to be close to the surface.It can be assumed that this favors a quick response of the terrain to rainfall in terms of movement.Regarding the track zone, the movement pattern of the stations 13, 53, and 54 does not follow a seasonal trend, but shows a more steady movement behavior.

DGNSS Monitoring Results
Figure 8 shows the mean annual horizontal and vertical velocity as well as the azimuth angle of the DGNSS benchmarks on the Corvara landslide.Station number 8 had been manipulated several times between October 2015 and November 2016 and therefore it had to be excluded.Time series of measured 3D displacements at 13 stations from the source, track, and accumulation zones are plotted in two different charts, according to the velocity classes of Cruden and Varnes [42] (Figure 8).The 3D displacements measured in 2015 and 2016 range from 176 mm (at station 6 in the accumulation zone of the landslide) to 4956 mm (at station 56 in the source zone of the landslide).
For the DGNSS stations in the source zone (23, 28, 49, 56, 57, and 58), a deceleration during winter and an acceleration during spring and summer months can be observed (Figure 9).The displacement follows the rainfall rates, which could be explained by the fact that at the upper zone of the landslide, the authors of [41] found the groundwater table to be close to the surface.It can be assumed that this favors a quick response of the terrain to rainfall in terms of movement.Regarding the track zone, the movement pattern of the stations 13, 53, and 54 does not follow a seasonal trend, but shows a more steady movement behavior.
Figure 10 shows a decreasing trend of velocity measured at some benchmarks (e.g., station 54).This might be attributed to the excavation of a drainage channel network, which was initiated in 2014 by the Corvara municipality and repeated on a yearly basis.However, more detailed surveys would be necessary to provide empirical evidence on the effectiveness of this drainage system.For many stations, with the exception of 56, a trend of deceleration can be observed over the shown period.Already, the authors of [37] have noticed that the landslide's kinematic intensity has decreased in general.Figure 10 shows a decreasing trend of velocity measured at some benchmarks (e.g., station 54).This might be attributed to the excavation of a drainage channel network, which was initiated in 2014 by the Corvara municipality and repeated on a yearly basis.However, more detailed surveys would be necessary to provide empirical evidence on the effectiveness of this drainage system.For many stations, with the exception of 56, a trend of deceleration can be observed over the shown period.Already, the authors of [37] have noticed that the landslide's kinematic intensity has decreased in general.

DInSAR and DGNSS Results Comparison
As an indication of measurement differences by ground-based instruments and DInSAR, the cumulative displacement of the stations 4 (located at the toe of the landslide) and 57 (located in the eastern depletion area) are represented in Figure 11.The CC of 0.35 with 2D and 3D PhU methods was used to create the plots.The frequency and acquisition dates of monthly DGNSS and SAR data were not identical; therefore, the DGNSS measurements were approximated by a linear model to make them comparable with the SBAS results.As observed previously, the 2D PhU gives

DInSAR and DGNSS Results Comparison
As an indication of measurement differences by ground-based instruments and DInSAR, the cumulative displacement of the stations 4 (located at the toe of the landslide) and 57 (located in the eastern depletion area) are represented in Figure 11.The CC of 0.35 with 2D and 3D PhU methods was used to create the plots.The frequency and acquisition dates of monthly DGNSS and SAR data were not identical; therefore, the DGNSS measurements were approximated by a linear model to make them comparable with the SBAS results.As observed previously, the 2D PhU gives measurements closer to the ground truth (red and green triangles in Figure 11).
As an indication of measurement differences by ground-based instruments and DInSAR, the cumulative displacement of the stations 4 (located at the toe of the landslide) and 57 (located in the eastern depletion area) are represented in Figure 11.The CC of 0.35 with 2D and 3D PhU methods was used to create the plots.The frequency and acquisition dates of monthly DGNSS and SAR data were not identical; therefore, the DGNSS measurements were approximated by a linear model to make them comparable with the SBAS results.As observed previously, the 2D PhU gives measurements closer to the ground truth (red and green triangles in Figure 11).
The cumulative LOS displacement and the validation plot for the estimated velocity using different coherence thresholds and PhU parameters are presented in Figures 12A and B, respectively.SBAS results underestimated the displacements of the stations 4 and 58; estimated the displacements of the stations 6, 49, and 57 with a close agreement; and estimated the displacements of the stations 11, 23, and 25 with a relatively low accuracy (Figure 12A).Generally, 2D and 3D PhU with the CC of 0.35 outperformed the other settings used.Compared to other parametric configurations, a 2D unwrapping in higher coherence areas (i.e., blue points in Figure 12B) better represents the real displacement of 8 selected DGNSS benchmarks according to measurement accuracy, while being located closer to the 1:1 line.The cumulative LOS displacement and the validation plot for the estimated velocity using different coherence thresholds and PhU parameters are presented in Figure 12A,B, respectively.SBAS results underestimated the displacements of the stations 4 and 58; estimated the displacements of the stations 6, 49, and 57 with a close agreement; and estimated the displacements of the stations 11, 23, and 25 with a relatively low accuracy (Figure 12A).Generally, 2D and 3D PhU with the CC of 0.35 outperformed the other settings used.Compared to other parametric configurations, a 2D unwrapping in higher coherence areas (i.e., blue points in Figure 12B) better represents the real displacement of 8 selected DGNSS benchmarks according to measurement accuracy, while being located closer to the 1:1 line.
The main potential causes of the discrepancy observed between the SBAS and DGNSS results (Figure 12) are (1) the high decorrelation due to the 12-day revisit time of S1, (2) the phase delay caused by tropospheric artifacts, (3) the nonlinear and non-LOS-oriented landslide movement, and (4) the DGNSS measurement errors due to possible antenna tilting.The comparison between DGNSS and SBAS displacement values projected in the LOS direction is presented in Table 2.The RMSE (Root Mean Square Error) value is maximal (i.e., equal to 9 mm) when considering a CC of 0.2 and 2D PhU and minimal (i.e., 6.1 mm) when considering a CC of 0.35 and 2D PhU.The RMSE value of 3D PhU with CC = 0.2 (i.e., 7.8 mm) is lower than 3D PhU with CC = 0.35 (i.e., 7.9 mm).From the PhU method point of view, 2D PhU generally outperformed 3D PhU where a higher coherence threshold was used, while 3D PhU provided a nearly identical result for both CC threshold thresholds with a relatively superior value of CC = 0.2.This means that 2D PhU is a coherence-sensitive operator (referring to the quality of pixel phase information) and 3D PhU is a pixel-dependent operator (referring to the quantity of pixel numbers), where the number of pixels is masked by setting the lower coherence threshold.

Discussion
DGNSS measurements show different velocities in (i) the area of depletion (e.g., 174-2633 mm/yr), with acceleration phases observed in late summer/fall 2015 and spring 2016; (ii) the track zone (e.g., 356-814 mm/yr), with a relatively constant displacement behavior; and (iii) the accumulation zone (e.g., 231-261 mm/yr), with an acceleration phase during the first half of 2016.This variation of the landslide displacement in relation to snow melt, precipitation, and groundwater conditions should be investigated more in depth.
The application of 2D PhU provided results closer to the ground-truth measurements, while 3D PhU has been proven to avoid phase jumps [21].Hooper and Zebker [49] also used a stepwise 3D PhU algorithm finding reasonable results, whereas for cases without multiple cycles of phase discontinuities (or jumps), no improvement was measured with 3D PhU.In this research, many pairs in the connection graph were discarded (Figure 3) due to the low coherence found in the landslide area (e.g., Figure 4F).Therefore, the 3D PhU, including temporal information to help unwrapping in low coherence interferogram areas, was not beneficial due to the low redundancy of the connection graph.Hence, the standard 2D unwrapping with a coherence threshold of 0.35 leads to fewer pixels, but to a more reliable pixel selection.A resulting lower effect of the eventual phase jumps leads to measurements closer to the ground truth.
The RMSE values indicate a low accuracy in several points and SBAS results are affected by the considerable errors, leading to an overall underestimation and a failure to detect fast-moving areas.Several factors led to propagation of the error into the results, which could considerably affect the precision of the DGNSS measurements and accuracy of the SBAS results.These error sources can be expressed as falling into the three following categories: (i) DGNSS-related error, (ii) non-LOS landslide displacement, and (iii) atmospheric artifacts.As the DGNSS antenna was mounted on the metal bar of the CR structure, slight changes resulting from the landslide movement led to the antenna tilting within the timespan of our study.Since the phase center variation (PCV) of the DGNSS antenna is calculated in the vertical position; any deviations from the vertical position of the antenna induce a bias in the measurement, leading to a low precision.This bias could manifest itself as a higher discrepancy (low accuracy) between the SBAS and DGNSS results (i.e., lower accuracy or higher RMSE values).This could be explained by the low accuracy obtained from SBAS results (e.g., stations 4, 11, 23, and 25).
The second error source is related to the SAR system geometry in relation to the real displacement field.The SBAS displacement has been estimated along the LOS of the satellite, while the direction of the landslide movement of depletion areas (in the Eastern and Northern parts) is not aligned to the LOS.Therefore, only a fraction of the occurred displacement is retrieved by the SAR system (e.g., stations 23 and 25).In particular, north-south-oriented movements of the northeastern part of the landslide are barely detectable by the DInSAR technique (e.g., station 58).Due to the complexity of the motion pattern, multi-aperture interferometry (MAI) techniques could be an alternative method to extract nearly north-south movements to be integrated to the LOS-based results [36].
Another important error source propagated in the SBAS results is the phase delay caused by the atmosphere that often affects the displacement signal.Indeed, the spatial and temporal variations of atmospheric properties over the case study cannot be precisely established by the temporal and spatial window sizes of the high-and low-pass filters in the SBAS processing.Therefore, the delay cannot be fully estimated and removed from interferograms by a filtering approach, especially in high regions with highly variable topography such as alpine areas, whereas the turbulent and stratified components of the atmosphere are constantly changing in time and space.In such circumstances, some sophisticated methods for atmospheric correction can be used, such as phase-based methods (e.g., linear model) and using auxiliary weather data such as from the Moderate Resolution Imaging Spectroradiometer (MODIS) and global atmospheric reanalysis data (i.e., ERA) provided by the ECMWF (European Center for Medium-Range Weather Forecasts) to mitigate the atmospheric artifacts.

Conclusions
Multi-temporal SBAS using Sentinel-1 imagery works relatively well for recognizing recently moving areas even if decorrelation might be high, mostly due to vegetation and the changing meteorological conditions in mountainous environments.The combination of a coherence threshold of 0.35 and 2D PhU gives results closer to the ground-truth measurements with maximal velocities of around 115 mm/yr.Therefore, using 2D PhU seems to minimize the total length of discontinuities in the unwrapped phase for this landslide site.However, no reliable InSAR velocities were detectable on south-facing slopes, because the displacement orientation is nearly perpendicular to the LOS direction.As the western part of the landslide is aligned to the LOS of the descending mode, including the national roads and urban infrastructures, this mode was intentionally selected.To overcome the limitations of the only 1D LOS measurements of the InSAR system, a combination of ascending and descending S1 data could provide a better displacement overview of the landslide kinematics.This could also reduce the uncertainty induced in the displacement map, especially for the eastern part of the landslide, which is mainly aligned to the LOS of the ascending mode.
During the period of time investigated, it was not possible to completely exploit both S1A and B, because S1B was not operational yet (at the time of data processing).Further analyses are ongoing combining both SBAS and PSI with complete Sentinel-1A/B time series in order to reduce the decorrelation occurring over time.New artificial 1-m edge CRs were installed on the landslide in late summer 2017.Despite an initial backscattering signal of 6.9 dB and monthly field campaigns to verify their orientations, their unavoidable tilting seems to spoil their systematic recognition over time.Therefore, further PSI processing using these CRs is expected to be challenging, mainly due to the landslide activity.The integration of proximal data covering parts of the landslide with remote sensing data is also tested for representing the complex changing slope kinematics assessment in (near) real-time.
Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 19X-band CRs were installed in 2013, and ground surface displacements were measured by multitemporal imagery coming from the Cosmo Sky-Med satellite[37].

Figure 1 .
Figure 1.Landslide monitoring network location and field impressions: (A) Landslide monitoring network; (B) monthly measurement station consisting of an X-band corner reflector and a support for the DGNSS antenna; (C) permanent DGNSS station, solar panel for power supply, and X-band corner reflector.DGNSS: differential global navigation satellite system.

Figure 1 .
Figure 1.Landslide monitoring network location and field impressions: (A) Landslide monitoring network; (B) monthly measurement station consisting of an X-band corner reflector and a support for the DGNSS antenna; (C) permanent DGNSS station, solar panel for power supply, and X-band corner reflector.DGNSS: differential global navigation satellite system.

Figure 2 .
Figure 2. Workflow of the Small BAseline Subsets (SBAS) processing for the slope displacement analysis.Atm. and IFG refer to atmospheric and interferograms, respectively.UTM: Universal Transverse Mercator

Figure 2 .
Figure 2. Workflow of the Small BAseline Subsets (SBAS) processing for the slope displacement analysis.Atm. and IFG refer to atmospheric and interferograms, respectively.UTM: Universal Transverse Mercator.

Figure 3 .
Figure 3. S1 SAR acquisition connections.(A) Perpendicular baselines (indicating a maximal spatial baseline of 174 m) between different image acquisitions according to their relative position (570 interferograms; yellow and green points indicate the master and slaves, respectively); (B) Delaunay connections for 3D phase unwrapping; (C) final connection graph after discarding pairs with a mean coherence less than 0.2 (the gray dashed lines present the snow time span, the black lines between acquisitions show the remaining connections, and the red dots indicate the discarded data); and (D) min/max (orange/red colors), mean-Std/mean+Std (green colors), and mean (blue color) coherence (within the landslide boundary) of the remaining pairs (i.e., 60 interferograms) after the graph editing.

Figure 3 .
Figure 3. S1 SAR acquisition connections.(A) Perpendicular baselines (indicating a maximal spatial baseline of 174 m) between different image acquisitions according to their relative position (570 interferograms; yellow and green points indicate the master and slaves, respectively); (B) Delaunay connections for 3D phase unwrapping; (C) final connection graph after discarding pairs with a mean coherence less than 0.2 (the gray dashed lines present the snow time span, the black lines between acquisitions show the remaining connections, and the red dots indicate the discarded data); and (D) min/max (orange/red colors), mean-Std/mean+Std (green colors), and mean (blue color) coherence (within the landslide boundary) of the remaining pairs (i.e., 60 interferograms) after the graph editing.

Figure 4 .
Figure 4.The effect of seasonal, temporal baseline, and surface scatterers' decorrelation on phase (left) and coherence (right).(A) Coherent interferogram showing the effect of the short baseline (i.e., 6 days) on the phase and (B) coherence related to the pair of 19-25 Oct 2016.(C) Phase values affected by decorrelation caused by vegetation and (D) coherence related to the pair of 6 June-3 July 2015.(E) Phase values affected by decorrelation caused by snow and (F) coherence related to the pair of 2-28 February 2016.The white polygon shows the boundary of the Corvara landslide.All images are presented in the SAR geometry.

Figure 4 .
Figure 4.The effect of seasonal, temporal baseline, and surface scatterers' decorrelation on phase (left) and coherence (right).(A) Coherent interferogram showing the effect of the short baseline (i.e., 6 days) on the phase and (B) coherence related to the pair of 19-25 Oct 2016.(C) Phase values affected by decorrelation caused by vegetation and (D) coherence related to the pair of 6 June-3 July 2015.(E) Phase values affected by decorrelation caused by snow and (F) coherence related to the pair of 2-28 February 2016.The white polygon shows the boundary of the Corvara landslide.All images are presented in the SAR geometry.

Figure 5 .
Figure 5. Coherence and displacement maps before the geocoding step.(A) Coherence map (threshold of 0.2).(B) Coherence map (threshold of 0.35).The positions of 8 DGNSS stations and 30 reference points selected for the refinement step (ramp and phase constant removal) are presented on the coherence map.(C) Cumulative line of sight (LOS) displacement map and (D) mean LOS velocity map (both created with the coherence threshold of 0.2 and 2D unwrapping method).The locations of 8 DGNSS stations and 30 reference points are presented on the displacement map.The images are presented in the SAR geometry.

Figure 5 .
Figure 5. Coherence and displacement maps before the geocoding step.(A) Coherence map (threshold of 0.2).(B) Coherence map (threshold of 0.35).The positions of 8 DGNSS stations and 30 reference points selected for the refinement step (ramp and phase constant removal) are presented on the coherence map.(C) Cumulative line of sight (LOS) displacement map and (D) mean LOS velocity map (both created with the coherence threshold of 0.2 and 2D unwrapping method).The locations of 8 DGNSS stations and 30 reference points are presented on the displacement map.The images are presented in the SAR geometry.

Figure 6 .
Figure 6.Comparison of interpolated SBAS LOS velocity maps (A) with an unwrapping coherence threshold (CC) of 0.35 and 2D PhU, (B) with CC of 0.35 and 3D PhU, (C) with CC of 0.2 and 2D PhU, and (D) with CC of 0.2 and 3D PhU.The spatial interpolation applied to the results of the DInSAR processing chain is meant to spatially preserve the results' continuity, avoiding discontinuities between pixels of low coherence.The figures show the terrain-corrected results in the UTM reference system.

Figure 7 .
Figure 7. Cross-section comparison of DInSAR results and landslide movement rates.CC represents the coherence thresholds and PhU the phase unwrapping process.The topography profile a-b is divided into (1) an accumulation area, (2) a first transition, (3) a transit area, (4) a second transition, and (5) a source area.Since the movement direction of the right side of the landslide, corresponding to the "b" profile, is not aligned to LOS (based on the DGNSS measurements of Figure8), the SBAS velocity for this part of the landslide was projected to the "b" profile alignment to avoid underestimating the velocity.The green and purple lines overlapped each other at the left part and then separated at the distance of 1500 m.

Figure 6 . 19 Figure 6 .
Figure 6.Comparison of interpolated SBAS LOS velocity maps (A) with an unwrapping coherence threshold (CC) of 0.35 and 2D PhU, (B) with CC of 0.35 and 3D PhU, (C) with CC of 0.2 and 2D PhU, and (D) with CC of 0.2 and 3D PhU.The spatial interpolation applied to the results of the DInSAR processing chain is meant to spatially preserve the results' continuity, avoiding discontinuities between pixels of low coherence.The figures show the terrain-corrected results in the UTM reference system.

Figure 7 .
Figure 7. Cross-section comparison of DInSAR results and landslide movement rates.CC represents the coherence thresholds and PhU the phase unwrapping process.The topography profile a-b is divided into (1) an accumulation area, (2) a first transition, (3) a transit area, (4) a second transition, and (5) a source area.Since the movement direction of the right side of the landslide, corresponding to the "b" profile, is not aligned to LOS (based on the DGNSS measurements of Figure8), the SBAS velocity for this part of the landslide was projected to the "b" profile alignment to avoid underestimating the velocity.The green and purple lines overlapped each other at the left part and then separated at the distance of 1500 m.

Figure 7 .
Figure 7. Cross-section comparison of DInSAR results and landslide movement rates.CC represents the coherence thresholds and PhU the phase unwrapping process.The topography profile a-b is divided into (1) an accumulation area, (2) a first transition, (3) a transit area, (4) a second transition, and(5) a source area.Since the movement direction of the right side of the landslide, corresponding to the "b" profile, is not aligned to LOS (based on the DGNSS measurements of Figure8), the SBAS velocity for this part of the landslide was projected to the "b" profile alignment to avoid underestimating the velocity.The green and purple lines overlapped each other at the left part and then separated at the distance of 1500 m.

Figure 8 .
Figure 8. Vector direction and velocity rate of DGNSS benchmarks for 2015 and 2016.

Figure 8 .
Figure 8. Vector direction and velocity rate of DGNSS benchmarks for 2015 and 2016.

19 Figure 10 .
Figure 10.3D displacement velocities between September 2013 and December 2016 at selected DGNSS stations.

Figure 10 .
Figure 10.3D displacement velocities between September 2013 and December 2016 at selected DGNSS stations.

Figure 11 .
Figure 11.Comparison of SBAS and DGNSS time series results.LOS cumulative displacement of the stations 4 and 57 is with a coherence threshold of 0.35.The DGNSS-fitted line (DGNSSL) is indicated in the figure, and the data-free area (i.e., from December 2015 to April 2016) indicates the snow period.SBAS-2D and -3D present the phase unwrapping method applied to the SBAS processing.

Figure 11 .
Figure 11.Comparison of SBAS and DGNSS time series results.LOS cumulative displacement of the stations 4 and 57 is with a coherence threshold of 0.35.The DGNSS-fitted line (DGNSSL) is indicated in the figure, and the data-free area (i.e., from December 2015 to April 2016) indicates the snow period.SBAS-2D and -3D present the phase unwrapping method applied to the SBAS processing.

Figure 11 .
Figure 11.Comparison of SBAS and DGNSS time series results.LOS cumulative displacement of the stations 4 and 57 is with a coherence threshold of 0.35.The DGNSS-fitted line (DGNSSL) is indicated in the figure, and the data-free area (i.e., from December 2015 to April 2016) indicates the snow period.SBAS-2D and -3D present the phase unwrapping method applied to the SBAS processing.

Figure 12 .
Figure 12.LOS DGNSS and SBAS results comparison.(A) Cumulative LOS displacement of the SBAS results versus the DGNSS measurement presented for each DGNSS benchmark.(B) Validation of the

Figure 12 .
Figure 12.LOS DGNSS and SBAS results comparison.(A) Cumulative LOS displacement of the SBAS results versus the DGNSS measurement presented for each DGNSS benchmark.(B) Validation of the SBAS velocity as a function of the chosen coherence threshold and PhU techniques.The error bars indicate the standard error for each displacement measurement with their related parameters.

Table 2 .
Comparison of DGNSS and SBAS LOS displacement (D) in mm.RMSE values show the accuracy of each unwrapping method with the different coherence threshold.