Integration of Digital Image Correlation of Sentinel-2 Data and Continuous GNSS for Long-Term Slope Movements Monitoring in Moderately Rapid Landslides

: This work explores the advantages and drawbacks of the application of Digital Image Correlation (DIC) to Sentinel-2 Multi Spectral Instrument (MSI) data in conjunction with continuous Global Navigation Satellite System (GNSS) monitoring. The goal is to retrieve a spatially distributed and long-term time-series of slope movements in large-scale moderately rapid landslides. The short revisit time of Sentinel-2 satellites (5 days since March 2017 and 10 days before) increases the availability of cloud and snow free satellite acquisitions of the area of interest, which is a prerequisite for the extrapolation of slope movement time-series using DIC techniques. Despite the Sentinel-2 limited spatial resolution, the derived long time-series can be integrated with—and validated by—continuous GNSS monitoring data. This allows to e ﬀ ectively monitor landslide movements that are too fast for the application of interferometric approaches. In this study, we used the Normalized Cross Correlation (NCC) digital image correlation technique by 51 Sentinel-2 MSI scenes (band 4 with 10 m spatial resolution), acquired between 19 February 2016 and 16 July 2019, to derive the slope movement time-series of the Ca’ Lita earthslide-earthﬂow in the northern Apennines (Italy). During the period considered, the landslide experienced two to three months-long phases of moderately rapid velocity (around 10 m / month) and, in between, prolonged periods of slow movements (approx. 10 cm / month). NCC results have been integrated with, and are compared to, time series from three continuous GNSS devices located in di ﬀ erent geomorphic zones of the landslide. On this basis, the errors and limitations associated to NCC time series are analysed and discussed together with their advantages and potentialities for assessing the spatial distribution and monitoring slope movements during moderately rapid reactivation events.


Case Study
The Ca' Lita landslide is located in the northern Apennines of Italy, in the Emilia-Romagna Region (Figure 1). The landslide is a roto-translational rock slide in the head zone, where flysch rock masses are involved, and evolves a downslope as a translational earthslide-earthflow affecting clayey complexes and debris material from the degradation of flysch rock masses [34]. The total length is approximately 3 km and the maximum width is about 1.4 km, with a maximum depth of the sliding surface of about 42 m measured by inclinometers in the head zone [35].
In April 2004, the Ca' Lita landslide underwent a first paroxysmal reactivation that caused significant retrogression and advancement of the landslide [36,37]. After that, the landslide was subject to the design and implementation of structural mitigation works [34][35][36][37], consisting of draining trenches, draining wells, check dams and pile-funded anchored retaining walls, which were mostly completed by 2011. The implementation of a monitoring network associated to deep drainage systems evidenced the influence of deep-groundwater inflow in the evolution of the landside [38]. In March 2016, a new paroxysmal reactivation of the landslide took place, resulting in displacements up to tens of meters in the upper landslide and up to more than hundred meters in the earthflow zone. The structural mitigation measures were severely damaged; retaining walls were overturned or compromised and part of the deep-drainage wells lost functionality. From a few days after reactivation, i.e., from 24 March 2016 until June 2016, a GNSS monitoring system, composed of three rover units and one reference master unit [1], was installed in order to support civil protection activities ( Figure 2) (rovers of this first monitoring period are indicated with the subscript "a"). The deceleration of the landslide during the following months was recorded by GNSS rovers. By the end of April 2019, cumulated planimetric displacements in the order of 44 m in the upper portion of the earthslide-earthflow (ROV2a) and 2.7 m in the lower portion of the earthslide-earthflow (ROV1a) were reached. In this same period, the roto-translational slide in the head zone did not show relevant movements (ROV3a). In December 2017 to April 2018, more than 40 m of movement took place in the earthslide-earthflow zone. This is known on the basis of the GNSS systems reinstalled in December 2016 that is still operative today (rovers in this second monitoring period are referred to with the subscript "b"). In May to June 2019, another large-scale paroxysmal reactivation took place. In the head zone, 9.6 m of planar and −2 m of vertical displacements were reached in the course of 12 h (ROV3b). Soon after, movements propagated downslope, cumulating up to 70 m planimetric displacement in the upper track zone (ROV2b) and up to 40 m in the lower track zone (ROV1b). The specific characteristics of the GNSS dataset will be presented in details in the "Methods and Data" paragraph.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 18 earthflow zone. This is known on the basis of the GNSS systems reinstalled in December 2016 that is still operative today (rovers in this second monitoring period are referred to with the subscript "b"). In May to June 2019, another large-scale paroxysmal reactivation took place. In the head zone, 9.6 m of planar and −2 m of vertical displacements were reached in the course of 12 h (ROV3b). Soon after, movements propagated downslope, cumulating up to 70 m planimetric displacement in the upper track zone (ROV2b) and up to 40 m in the lower track zone (ROV1b). The specific characteristics of the GNSS dataset will be presented in details in the "Methods and Data" paragraph.

Methods and Data
In this work, the Normalized Cross Correlation procedure by Kääb, Vollmer and Heid [40,41], implemented in the Correlation Image Analysis (CIAS) software has been used. The CIAS software works with pairs of precisely co-registered images in geotiff or tiff-world format acquired from any sensors. Images must have equal vertical and horizontal resolution and need to be single channel 8-bit. The dataset used in this work includes 51 Sentinel-2 Multi Spectral Instrument (MSI) Level-2A acquisitions with a median temporal granularity of 21 days from February 2016 to July 2019 (see Table 1 for details and Figure 3a for processing flowchart). Such data, freely distributed through the SciHub portal of the European Space Agency (ESA) [42], are ortho-images in UTM/WGS84 projection provided as Bottom of Atmosphere (BOA) reflectance images derived from the associated Level-1C products [43]. In this work, data belonging to band 4 of the MSI sensor, characterised by a 10 m spatial resolution, have been used. By iteratively applying the NCC analysis to subsequent pairs of images, a total of 50 displacement maps have been obtained, from which displacement time-series could be derived for every pixel within the AOI. By taking into consideration the dimensions of the landslide (minimum width of 150 m in the track zone) and the expected maximum displacement between two scenes, the settings adopted in CIAS have been a reference block size of 100 m and search area size of 150 m, which allow for maximum theorical planimetric displacements between the two analysed scenes of 70 m. Prior to NCC analysis, the co-registration pre-analysis routine inside CIAS has been run using 14 points outside the landslide as stable benchmarks ( Figure 2). More in detail, these benchmarks are represented by elements at ground height (in order to avoid any parallax bias) such as road intersections and corners of cultivation fields located outside the landslide that show no evidence of movements. This allows us to bypass limitations related to the absolute geolocation uncertainty of 20 m without Ground Control Points (GCP) and 12.5 m with GCPs and a 0.3 pixel multi-temporal registration [43]. scenes of 70 m. Prior to NCC analysis, the co-registration pre-analysis routine inside CIAS has been run using 14 points outside the landslide as stable benchmarks ( Figure 2). More in detail, these benchmarks are represented by elements at ground height (in order to avoid any parallax bias) such as road intersections and corners of cultivation fields located outside the landslide that show no evidence of movements. This allows us to bypass limitations related to the absolute geolocation uncertainty of 20 m without Ground Control Points (GCP) and 12.5 m with GCPs and a 0.3 pixel multi-temporal registration [43].  Continuous GNSS receivers have been operated in the landslide during most of the period analyses with the NCC of Sentinel 2 datasets. GNSS receivers were located at the head zone (ROV3), the upper track zone (ROV2) and the lower track zone (ROV1). The GNSS monitoring dataset covers two time intervals: period "a" from 24 March 2016 to 1 June 2016; period "b" from 28 November 2016 to 15 July 2019 ( Figure 3b). During period "a", the GNSS monitoring array was composed by: an external reference station (double frequency Leica L1-L2 receiver mod. GMX902), and three monitoring rovers (Leica L1 receivers mod. GMX901) indicated as ROV1a, ROV2a and ROV3a. During period "b" the GNSS monitoring array was composed by: two reference stations (a Leica L1-L2 receiver mod. GMX902 in exactly the same position as during period "a", plus an Emlid L1 receiver model Reach RS installed right next to it, Figure 2a), and three monitoring rovers (two Leica L1 receivers mod. GMX901, indicated as ROV1b and ROV2b, plus an Emlid L1 receiver model "Reach RS" indicated as ROV3b) (Figure 4). With the Leica GNSS array, positioning solutions are computed on the basis of rinex products of different duration (1 h, 30 min and 10 min), and are post-processed by means of the "Leica Spider Software" installed in an industrial PC in the field. Position data are then sent to a remote data server. With the Emlid GNSS array, the ROV3b computes 5 Hz Real Time Kinematic (RTK) solutions by using Radio Technical Commission for Maritime Services (RTCM) stream corrections sent by the Emlid reference station. RTK positioning solutions are calculated by the receiver, and hourly positioning solutions are sent to a remote data server. This type of RTK array, is a cost-effective alternative to the use of high-end Geodetic receivers [44]. A detailed description of the GNSS positioning dataset and of the technical characteristics of the arrays is given in Mulas et al. [39]. The monthly stability check of the reference station is presented in Appendix A. In the blue square is the Emlid single frequency low cost receiver (model Reach RS) acting as a reference station. (b) Reach RS receiver is installed as ROV3b rover unit and its photovoltaic power unit. (c) Leica GMX901 receiver is installed as ROV2b rover unit and its photovoltaic power unit. During period "a", ROV1a, ROV2a and ROV3a were equipped with Leica GMX901 receivers. During period "b", ROV1b and ROV2b were equipped with Leica GMX901 receivers, while ROV3b was equipped with the low-cost Emlid Reach RS.

Results
By applying NCC to Sentinel 2 data, a total of 50 deformation maps of the Ca' Lita landslide have been obtained for the period February 2016 to July 2019 with a median periodicity of 21 days (Table 1). The resulting displacements maps allow us to identify changes occurring in the AOI, such as landslide activity or the Secchia river channel migration ( Figure 5). We recall that Sentinel-2 scenes were selected by discarding images with cloud cover above the active landslide area. Despite that, cloud cover effect in the remaining frame is still sometimes an issue (Figure 5d). Taking this into consideration, and in order to focus the attention on the Ca' Lita landslide dynamics, in Figure 6 and  In the blue square is the Emlid single frequency low cost receiver (model Reach RS) acting as a reference station. (b) Reach RS receiver is installed as ROV3b rover unit and its photovoltaic power unit. (c) Leica GMX901 receiver is installed as ROV2b rover unit and its photovoltaic power unit. During period "a", ROV1a, ROV2a and ROV3a were equipped with Leica GMX901 receivers. During period "b", ROV1b and ROV2b were equipped with Leica GMX901 receivers, while ROV3b was equipped with the low-cost Emlid Reach RS.
The analysis of the DIC displacement time-series (DIC displ ) has been carried out by comparing them to the planimetric continuous GNSS displacement time-series (GNSS displ ), derived by computing for each rover the median value of GNSS positioning solutions on 1-day long time windows. The DIC displ time series have been retrieved by computing, for each Sentinel-2 acquisition time-step, the mean value of DIC displacement (direction and magnitude) in 60 × 60 m computation matrices centred on the median coordinate values of each GNSS rover. Consequently, the DIC displ time series have been created for the head zone (ROV3a and ROV3b), upper track zone (ROV2a and ROV2b) and lower track zone (ROV1a and ROV1b). The discrepancy between the DIC displ and GNSS displ time series at variable minimum displacement thresholds has been quantified using the relative % error analysis and the linear fit between time series.

Results
By applying NCC to Sentinel 2 data, a total of 50 deformation maps of the Ca' Lita landslide have been obtained for the period February 2016 to July 2019 with a median periodicity of 21 days (Table 1). The resulting displacements maps allow us to identify changes occurring in the AOI, such as landslide activity or the Secchia river channel migration ( Figure 5). We recall that Sentinel-2 scenes were selected by discarding images with cloud cover above the active landslide area. Despite that, cloud cover effect in the remaining frame is still sometimes an issue (Figure 5d). Taking this into consideration, and in order to focus the attention on the Ca' Lita landslide dynamics, in Figures 6 and 7  The displacements map referring to February-March 2016, bracketing the onset of the March 2016 reactivation event, is presented in Figure 6 together with "rose-wind like" plots of displacement vectors, to highlight preferential directions of movement within morphological zones. The map (masked outside the active landslide) evidences higher displacement rates in the head zone and in the higher track zone. Conversely, the lower track zone did not experience relevant displacements. The rose diagrams show that the recorded directions of DIC displacement vectors are coherent with the expected kinematic of the landslide in the different zones. A selection of displacement maps referring to different activity phases in the entire analysis period is presented in Figure 7  The displacements map referring to February-March 2016, bracketing the onset of the March 2016 reactivation event, is presented in Figure 6 together with "rose-wind like" plots of displacement vectors, to highlight preferential directions of movement within morphological zones. The map (masked outside the active landslide) evidences higher displacement rates in the head zone and in the higher track zone. Conversely, the lower track zone did not experience relevant displacements. The rose diagrams show that the recorded directions of DIC displacement vectors are coherent with the expected kinematic of the landslide in the different zones. A selection of displacement maps referring to different activity phases in the entire analysis period is presented in Figure 7   The results of the DICdispl time-series extraction with the 60 × 60 m computation matrices centred on the median coordinates values of each GNSS rover are presented in Figure 8. An issue with the construction of the DICdispl time-series is that errors can cumulate at each time step resulting, in the long-term, in unreliable cumulated displacement values. In order to assess this issue, prior to DICdispl time-series extraction, we have filtered the DIC movement maps by using different minimum displacement thresholds ranging from 1 to 10 m. This has allowed us to generate different DICdispl time-series, at variable filtering thresholds, for slope movements around the GNSS (Figure 8). Results show, qualitatively, that the higher the displacement threshold is, the better is the fit between the DICdispl and the GNSSdispl time-series. For instance, by adopting a displacement threshold of ≥10 m (i.e., the pixel size of Sentinel data), the DICdispl over the entire period of analysis reach totals of around 290 m in the upper track zone, 90 m in the lower track zone and 20 m in the head zone, which are quantitatively comparable with GNSSdispl results. On the other hand, by using minimum displacement thresholds ˂10 m (i.e., sub-pixel accuracies), the DICdispl time series tend to overestimate movements, as errors are cumulated over time. In the high mobility GNSS rovers ROV1 and ROV2, the offset between the DICdispl and the GNSSdispl time-series becomes quantitatively significant for minimum displacement thresholds lower than 5 to 7 m. In the low mobility ROV3, the discrepancy is already evident for minimum displacement thresholds lower than 8 m. The results of the DIC displ time-series extraction with the 60 × 60 m computation matrices centred on the median coordinates values of each GNSS rover are presented in Figure 8. An issue with the construction of the DIC displ time-series is that errors can cumulate at each time step resulting, in the long-term, in unreliable cumulated displacement values. In order to assess this issue, prior to DIC displ time-series extraction, we have filtered the DIC movement maps by using different minimum displacement thresholds ranging from 1 to 10 m. This has allowed us to generate different DIC displ time-series, at variable filtering thresholds, for slope movements around the GNSS (Figure 8). Results show, qualitatively, that the higher the displacement threshold is, the better is the fit between the DIC displ and the GNSS displ time-series. For instance, by adopting a displacement threshold of ≥10 m (i.e., the pixel size of Sentinel data), the DIC displ over the entire period of analysis reach totals of around 290 m in the upper track zone, 90 m in the lower track zone and 20 m in the head zone, which are quantitatively comparable with GNSS displ results. On the other hand, by using minimum displacement thresholds <10 m (i.e., sub-pixel accuracies), the DICd ispl time series tend to overestimate movements, as errors are cumulated over time. In the high mobility GNSS rovers ROV1 and ROV2, the offset between the DIC displ and the GNSS displ time-series becomes quantitatively significant for minimum displacement thresholds lower than 5 to 7 m. In the low mobility ROV3, the discrepancy is already evident for minimum displacement thresholds lower than 8 m.  A quantitative assessment of the errors associated with the DIC displ time-series retrieval at sub-pixel resolution has been performed by computing the relative % error between DIC displ and GNSS displ for each minimum displacement threshold considered. The plot in Figure 9a evidences that relative errors are way above 100% for minimum displacement thresholds lower than 2 m (i.e., 0.2 pixels), while they are reduced to less than 15% for minimum displacement thresholds higher than 6 m. In Figure 9b, the coefficients of determination (R 2 ) for the linear fits between DIC displ and GNSS displ at different minimum displacement thresholds are presented together with the resulting equations (y 1 to y 10 in relation to thresholds from 1 to 10 m). The R 2 indexes are all relatively good, between 0.90 and 0.93, evidencing that even if the relative errors might be very high, the discrepancy between the absolute value of the measured movements is limited even at sub-pixel resolutions. The frequency distribution of relative errors associated with DIC displ at different minimum displacement thresholds is presented in Figure 9c. The results are that the median values of relative errors vary in three groups inside a quite small range. One group is that of results associated to thresholds <3 m, with the higher median relative errors and the larger interquartile ranges (IQR). Another group is associated to thresholds from 3 to 7 m, with slightly lower median error values and a smaller IQR. Finally, a third group is that of results for thresholds higher than 7 m, with the lowest median relative error values and, also, the lowest IQR.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 18 A quantitative assessment of the errors associated with the DICdispl time-series retrieval at subpixel resolution has been performed by computing the relative % error between DICdispl and GNSSdispl for each minimum displacement threshold considered. The plot in Figure 9a evidences that relative errors are way above 100% for minimum displacement thresholds lower than 2 m (i.e., 0.2 pixels), while they are reduced to less than 15% for minimum displacement thresholds higher than 6 m. In Figure 9b, the coefficients of determination (R 2 ) for the linear fits between DICdispl and GNSSdispl at different minimum displacement thresholds are presented together with the resulting equations (y1 to y10 in relation to thresholds from 1 to 10 m). The R 2 indexes are all relatively good, between 0.90 and 0.93, evidencing that even if the relative errors might be very high, the discrepancy between the absolute value of the measured movements is limited even at sub-pixel resolutions. The frequency distribution of relative errors associated with DICdispl at different minimum displacement thresholds is presented in Figure 9c. The results are that the median values of relative errors vary in three groups inside a quite small range. One group is that of results associated to thresholds <3 m, with the higher median relative errors and the larger interquartile ranges (IQR). Another group is associated to thresholds from 3 to 7 m, with slightly lower median error values and a smaller IQR. Finally, a third group is that of results for thresholds higher than 7 m, with the lowest median relative error values and, also, the lowest IQR.  . Error analysis: (a) relative error of DIC displ as a function of GNSS displ , where the grey area is the 95% confidence band for the red dashed regression curve; (b) DIC displ vs. GNSS displ : y 1 to y 10 are the equations resulting from linear fit of points belonging to thresholds from 1 to 10 m and the respective coefficient of determination. Displacements smaller than 1 m are excluded; (c) Percentage error according to displacement threshold classes. Boxplots ranges: extremes of the boxes are the 1st and 3rd quartile (Q 1 and Q 2 ). The black lines depict median values, the lower whisker is "Q 1 -1.5 × IQR" and the upper whisker is "Q 2 -1.5 × IQR" (where IQR stands for "interquartile range" and is the difference between Q 2 and Q 1 ).

Discussion
The set of results obtained indicates that this technique can be a useful tool to monitor the dynamics of large-scale slope movements characterised by moderate velocity ranges. With this respect, the following issues can be discussed.

•
Quantitative errors associated to displacement values retrieved by NCC. Our results prove that NCC movement estimates between two subsequent scenes are: (i) affected by large relative errors for movements lower than 0.3 pixels; (ii) affected by significant errors-yet indicative of ongoing slope dynamics trends-for movements between 0.3 to 0.7 pixels; (iii) quite reliable for movements larger than 0.7 pixels, providing time-series that can be used to monitor slope movements evolution over long periods of time, with median relative errors in the range of 11.66%. The error analysis shows that the CIAS algorithm applied to Sentinel-2 band 4 scenes reaches good performances (approx. 10% error) when displacements are of the same order of the scene geometric resolution (10 m). These results are somehow at odds with the conclusion of other authors, which have reported higher sub-pixel accuracies. Bickel et al. [16], by presenting displacement maps obtained by DIC coherent with ones derived by using GNSS, argue that the velocity of the observed displacement did not seem to have an influence on the correlation accuracy. Mazzanti at al. [25] report sub-pixel accuracies using a DIC methodology based on redundant correlations with multiple master images, and declare precisions up to 1% of the pixel size. Stumpf et al. [30] also reports accuracies in between 1/9th and 1/20th of a pixel by using a two-step method that calibrates 1st grade DIC displ residuals against ground measurements using High Performance Computing (HPC) infrastructures.

•
Applicability of NCC to large displacements and moderate velocity ranges. Results prove that NCC goes beyond the range of landslide velocities that are normally retrieved with the application of other satellite-based techniques, such as Synthetic Aperture Radar (SAR) interferometry. This study indicates that it is possible to detect displacements during activation phases of large landslides. This is true even with displacements of several tens of meters and velocities in ranges from slow to moderate according to landslide movement rate scale by Cruden and Varnes [45]. Similar findings were achieved by other authors in monitoring ice flow measurements in glaciers [28].

•
Identification of spatially distributed displacement patterns. Another advantage is that the output of DIC is spatial displacements maps that can reveal interesting information on the reactivation patterns of large landslides that otherwise are generally monitored in a few specific points only. Results shown in Figure 5 demonstrate that outside the landslide area, the movements computed by the DIC algorithm are substantially nil, strengthening the interpretation of movement patterns retrieved inside the landslide area.

•
Integration with in situ monitoring systems. As a matter of fact, when displacement rates are significant, the maintenance in operation of in situ-systems is quite demanding and potentially very costly (as GNSS or other systems can possibly be lost during the landslide activity). In practice, with NCC, we have been able to obtain long-term time-series that can be used to extend existing-often discontinuous-monitoring time-series from in situ instruments.

•
Main limitations in the application of DIC techniques to Sentinel-2 scenes and the individuation of alternative sensing solutions. One limitation that it is worthwhile to recall is that the MSI onboard the Sentinel-2 constellation is susceptible to cloud cover in the AOI. Even if the 5-day acquisition frequency can help to reduce this issue, this limits the number of MSI acquisitions that can effectively be used in order to retrieve displacements rates by the application of DIC techniques. Even MSI acquisition that is cloud-free over the landslide area can, however, produce false displacement information in areas outside the landside where cloud coverage exists (see for example Figure 5d). Moreover, the presence of snow cover also inhibits the application of DIC algorithms, making the identification of movements during snow-cover or snow melt phases quite problematic. Such issue could be tackled by using satellite constellation with sub metric resolution Remote Sens. 2020, 12, 2605 13 of 17 and up to 12 h revisit time (i.e., "Planet Scope" from "planet" company) [25]. An alternative solution to reduce the influence of cloud cover could be represented by the use of the SAR amplitude scenes. The use of high-resolution COSMO-SkyMed or TerraSAR-X X-band SAR constellations has been already proven by several authors by applying DIC or radargrammetry algorithms [17,26]. Despite the technical advantages of X-band SAR and of high spatial and temporal resolution earth observing satellite constellations, we should recall that they are commercial products while Sentinel-2 scenes are freely distributed by the ESA.

•
Possibility to identify other geomorphic processes using DIC algorithms. The application to the entire AOI has highlighted that the adopted technique can also be used to reveal river dynamics. Results in Figure 5 evidence changes along the Secchia river, which can be ascribed to bedload transport and channel migration. The analysis of such fluvial dynamics goes beyond the scope of this paper.

Conclusions
This paper presents a successful application of DIC applied to Sentinel-2 MSI scenes to retrieve slope movement maps and displacement time-series for a large-scale moderate-velocity landslide during a more than three years period. During this time interval, three major reactivations have been back analysed, which resulted in cumulated displacements up to 300 m in the upper track zone. The relative discrepancy between displacements assessed using DIC compared to continuous GNSS measurements has an 11.66% median value using minimum displacement resolutions of 0.7 pixels. This allows us to retrieve time-series with a good fit with displacements obtained by in situ GNSS sensors. In perspective, thanks to the availability of data acquired by Sentinel-2 and Sentinel-1 constellations with high temporal frequency and global coverage, the integrated use of interferometry based (differential interferometry and multi-temporal interferometry) and DIC processing techniques can allow us to bypass flaws of both techniques by fully exploiting the specific ranges of movement measurement capabilities. In the case of the Ca' Lita landslide, DIC can be applied to monitor rapid accelerations that cannot be resolved by interferometry based techniques, while SAR interferometry can be used to monitor slow movements phases. For future developments, the application of redundant NCC algorithms [25,30] to the Sentinel-2 stack used in this work would be of great interest in order to measure the accuracy improvement that could be achieved by using these techniques. Future research will also focus on the error estimation of DIC applied to Landsat-7 scenes (15 m, panchromatic band) in order to verify the feasibility of analysing historical datasets and following the evolution of large landslides along decades.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In order to verify the stability of the area where both the GNSS reference stations are located, coordinates of the Leica station have been periodically computed. Differential solutions with respect to a permanent GNSS station belonging to the "Rete Integrata Nazionale GNSS" (RING) national network [46] are obtained. More in detail, daily solutions have been processed, in the period between December 2016 and July 2019, on a monthly basis using as reference the MODE permanent station (30 km NE from Ca' Lita landslide). Solutions ( Figure A1) show a small standard deviation (around 0.01 m for the planimetric components and 0.02 m for the elevation) with respect to their median value (Table A1).
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 18 December 2016 and July 2019, on a monthly basis using as reference the MODE permanent station (30 km NE from Ca' Lita landslide). Solutions ( Figure A1) show a small standard deviation (around 0.01 m for the planimetric components and 0.02 m for the elevation) with respect to their median value (Table A1).