Long-Term Satellite Monitoring of the Slumgullion Landslide Using Space-Borne Synthetic Aperture Radar Sub-Pixel Offset Tracking

Kinematic characterization of a landslide at large, small, and detailed scale is today still rare and challenging, especially for long periods, due to the difficulty in implementing demanding ground surveys with adequate spatiotemporal coverage. In this work, the suitability of space-borne synthetic aperture radar sub-pixel offset tracking for the long-term monitoring of the Slumgullion landslide in Colorado (US) is investigated. This landslide is classified as a debris slide and has so far been monitored through ground surveys and, more recently, airborne remote sensing, while satellite images are scarcely exploited. The peculiarity of this landslide is that it is subject to displacements of several meters per year. Therefore, it cannot be monitored with traditional synthetic aperture radar differential interferometry, as this technique has limitations related to the loss of interferometric coherence and to the maximum observable displacement gradient/rate. In order to overcome these limitations, space-borne synthetic aperture radar sub-pixel offset tracking is applied to pairs of images acquired with a time span of one year between August 2011 and August 2013. The obtained results are compared with those available in the literature, both at landslide scale, retrieved through field surveys, and at point scale, using airborne synthetic aperture radar imaging and GPS. The comparison showed full congruence with the past literature. A consistency check covering the full observation period is also implemented to confirm the reliability of the technique, which results in a cheap and effective methodology for the long-term monitoring of large landslide-induced movements.


Introduction
The kinematic estimation of landslide deformation is an important task for investigating the mechanisms affecting their evolution.An improved understanding of changes in landslide behavior as a function of predisposing factors (such as slope, aspect, land use, lithology, etc.) and triggering factors (such as intense or prolonged rainfall, seismicity, human actions, etc.) is crucial for risk mitigation, and thus brings significant social, environmental, and financial benefits [1].
Historically, landslides have been investigated with sparse point measurements acquired through demanding on-field campaigns using traditional instruments (inclinometers, piezometers topographic leveling) and/or GPS.However, recent years have been characterized by a rising interest in and exploitation of satellite technologies for the monitoring of active ground deformation, including landslide motion and hillslope creep, thus allowing for the overcoming of the insufficient spatiotemporal coverage achievable through traditional monitoring methods.In this sense, Differential Synthetic Aperture Radar Interferometry (DInSAR) techniques have largely demonstrated their effectiveness in measuring surface motion and deformation at different scales [2][3][4][5][6].The spread of this methodology is mainly related to Permanent Scatterers (PS) [7] and Small BAseline Subset (SBAS) [3] approaches, and their combination (known as SqueeSAR) [8][9][10].However, these methods are subject to some restrictions, such as the presence of vegetation (causing temporal decorrelation even within very short time intervals) and the 1-D line of sight (LOS) measurement sensitivity, which is also constrained by the adopted wavelength relatively to the maximum measurable displacement rate [11,12].
Today, a new class of methods can provide information complementary to that derived from DInSAR by working with the amplitude channel.These techniques, based on Sub-Pixel Offset Tracking (SPOT) [11,13], allow the measurement of displacements in the South-North and East-West directions without any limitation on the observable rate.This means that a pair of SAR images can be used to detect movements of several meters with a good degree of approximation [13].SPOT methods are generally less precise than conventional DInSAR methods [11], however they are less sensitive to atmospheric effects and are not strictly only applicable to highly coherent targets, i.e., they can measure displacements even in densely vegetated areas [14].
SPOT methods have been successfully applied to estimate movements of glaciers and terrain [15][16][17] due to natural phenomena (landslides) [11,13] and/or human activities (subsidence induced by underground excavation) [18,19].In this work, the suitability of this methodology for the long-term monitoring of the Slumgullion landslide (Colorado, US) is investigated.This landslide has been extensively studied in the past due to the fast/slow related displacements [20][21][22][23].Previous investigations have mainly been implemented with field sensors.Starting from the early 2000s, synthetic aperture radar (SAR) airborne remote sensing was also employed.The first campaign was implemented by the Brigham Young University using an X-band sensor, with the purpose of correlating the measured movements to the soil water content [24].More recently, some studies exploiting airborne L-band remote sensing have been presented [25][26][27].However, so far, the use of satellite technologies in the literature was limited, despite they can provide a cost-effective and continuous update of the landslide state.In fact, kinematic characterization at landslide scale is today still rare, especially for long periods, due to the difficulty in implementing demanding ground surveys with adequate spatiotemporal coverage [22].
To this end, state-of-the-art SAR SPOT was applied to three X-band COSMO-SkyMed spotlight images with about one-meter spatial resolution to monitor the Slumgullion landslide over a time frame of two years, from August 2011 to August 2013.The landslide is investigated at both large and small scale, and, for the first time, the displacements retrieved using very high-resolution space-borne images are validated against ground data provided in the past literature.
The work is organized as follows.The case study and available data are presented in Section 2. Experimental results are presented and validated in Section 3, and discussed in Section 4. Conclusions are drawn at the end of the work, in Section 5.

Study area and available data
The study area concerns the Slumgullion landslide, US.It is depicted in Figure 1, with a LiDAR Digital Elevation Model (DEM) with a 0.5-meter spatial resolution superimposed and developed by the US National Center for Airborne Laser Mapping (NCALM).The landslide is located in the San Juan Mountains in southwestern Colorado, near the town of Lake City.It has been moving for about 350 years, with a maximum measured velocity of six meters per year [21].According to [28], it has been classified as a debris flow which involves deeply weathered tertiary volcanic rocks.It extends for about seven kilometers from the Cannibal Plateau to Lake San Cristobal, with a mean width and depth of 300 meters and 14 meters, respectively [20].
Inside the landslide, an area of about one square kilometer (see black curve in Figure 1) is still active, and is characterized by a low ground-surface inclination (about eight degrees).Based on its average displacement rate (between 0.1-2.0cm/day as reported in [2]), it can be considered a very slow landslide according to [28].The study area concerns the Slumgullion landslide, US.It is depicted in Figure 1, with a LiDAR Digital Elevation Model (DEM) with a 0.5-meter spatial resolution superimposed and developed by the US National Center for Airborne Laser Mapping (NCALM).The landslide is located in the San Juan Mountains in southwestern Colorado, near the town of Lake City.It has been moving for about 350 years, with a maximum measured velocity of six meters per year [21].According to [28], it has been classified as a debris flow which involves deeply weathered tertiary volcanic rocks.It extends for about seven kilometers from the Cannibal Plateau to Lake San Cristobal, with a mean width and depth of 300 meters and 14 meters, respectively [20].
Inside the landslide, an area of about one square kilometer (see black curve in Figure 1) is still active, and is characterized by a low ground-surface inclination (about eight degrees).Based on its average displacement rate (between 0.1-2.0cm/day as reported in [2]), it can be considered a very slow landslide according to [28].As observed in [29] and [30], the Slumgullion landslide is composed of multiple kinematic elements, each of them moving like a rigid block sliding along faults.Accordingly, the landslide has been segregated into 11 primary kinematic elements (see labels in Figure 1) that can be assumed as temporally consistent [22].As observed in [29,30], the Slumgullion landslide is composed of multiple kinematic elements, each of them moving like a rigid block sliding along faults.Accordingly, the landslide has been segregated into 11 primary kinematic elements (see labels in Figure 1) that can be assumed as temporally consistent [22].
Between 1985 and 1990, velocity of elements from 1 to 4, constituting the flattest part of the landslide, increased nearly linearly in the downslope direction, as reported in [18,20,21].Displacement rate abruptly increased from element 4 to 5 and increased almost linearly up to element 7. The latter was the fastest and steepest element of the landslide.Then, it widens and flattens, slowing downslope to elements 8 and 9. Finally, the landslide moved along the oldest ground surface forming elements 10 and 11, whose speed was much slower than that of elements 8 and 9.As indicated in [31], these regions can be grouped into head (Region 1 to 4), upper neck (Region 5 to 7), lower neck (Region 8 and 9), and toe (Region 10 and 11) (see box (b) in Figure 1).
Satellite monitoring has been implemented by exploiting three COSMO-SkyMed spotlight images with about one-meter spatial resolution.The dataset covers a time frame of two years, from August 2011 to August 2013.Images were acquired with one-year frequency.They were combined in two co-registered pairs, the first covering the period August 2011-August 2012 and the second covering the period August 2012-August 2013.A third pair, covering the time frame August 2011-August 2013, was also considered to implement a consistency check of the obtained results (see Section 3).

Implemented Sub-Pixel Offset Tracking (SPOT) technique
The flow diagram of the implemented methodology is depicted in Figure 2. In particular, in the upper part, the general processing chain, from the input to the final displacement vector field, is reported.In the lower part, an exploded view of the SPOT processing block is displayed.Between 1985 and 1990, velocity of elements from 1 to 4, constituting the flattest part of the landslide, increased nearly linearly in the downslope direction, as reported in [18], [20] and [21].Displacement rate abruptly increased from element 4 to 5 and increased almost linearly up to element 7. The latter was the fastest and steepest element of the landslide.Then, it widens and flattens, slowing downslope to elements 8 and 9. Finally, the landslide moved along the oldest ground surface forming elements 10 and 11, whose speed was much slower than that of elements 8 and 9.As indicated in [31], these regions can be grouped into head (Region 1 to 4), upper neck (Region 5 to 7), lower neck (Region 8 and 9), and toe (Region 10 and 11) (see box (b) in Figure 1).
Satellite monitoring has been implemented by exploiting three COSMO-SkyMed spotlight images with about one-meter spatial resolution.The dataset covers a time frame of two years, from August 2011 to August 2013.Images were acquired with one-year frequency.They were combined in two co-registered pairs, the first covering the period August 2011-August 2012 and the second covering the period August 2012-August 2013.A third pair, covering the time frame August 2011-August 2013, was also considered to implement a consistency check of the obtained results (see Section 3).

Implemented Sub-Pixel Offset Tracking (SPOT) technique
The flow diagram of the implemented methodology is depicted in Figure 2. In particular, in the upper part, the general processing chain, from the input to the final displacement vector field, is reported.In the lower part, an exploded view of the SPOT processing block is displayed.General processing starts from standard co-registration [32].If precise orbit data are available, an only orbit-based co-registration can be implemented.Otherwise, full co-registration (i.e., including cross-correlation and coherence refinement steps) is suggested, provided that the scene is much larger than the area of interest.In this case, Ground Control Points (GCPs) located on the landslide will be automatically discarded since they typically exhibit very low cross-correlation and interferometric coherence values.The first acquired image is assumed to be the master image, i.e., the reference for the estimation of displacements.The image in which displacements are evaluated is the slave image.
The SPOT algorithm exploits cross-correlation measures on several windows extracted from the co-registered image pair to estimate the shift between the master patch and the slave patch.Windows are extracted around grid points usually regularly distributed across the images (see first two blocks of the lower diagram in Figure 2).General processing starts from standard co-registration [32].If precise orbit data are available, an only orbit-based co-registration can be implemented.Otherwise, full co-registration (i.e., including cross-correlation and coherence refinement steps) is suggested, provided that the scene is much larger than the area of interest.In this case, Ground Control Points (GCPs) located on the landslide will be automatically discarded since they typically exhibit very low cross-correlation and interferometric coherence values.The first acquired image is assumed to be the master image, i.e., the reference for the estimation of displacements.The image in which displacements are evaluated is the slave image.
The SPOT algorithm exploits cross-correlation measures on several windows extracted from the co-registered image pair to estimate the shift between the master patch and the slave patch.Windows are extracted around grid points usually regularly distributed across the images (see first two blocks of the lower diagram in Figure 2).
The cross-correlation between two null-mean patches M and S, taken, respectively, on the master and slave image, is computed as follows: where FFT and IFFT indicate the fast Fourier transformation and the inverse fast Fourier transformation, respectively, the apex * the complex conjugation operation, and the symbol • the mean operator.
In this equation M and S are oversampled by a factor f (which must be a power of two in order to optimize the FFT calculation) to take into account sub-pixel movements, being the minimum detectable displacement (in pixel units) equal to 1/f.In Equation ( 1), C is a matrix.Its maximum identifies the amount of shift to be applied to the slave patch to have it superimposed onto the master patch.The higher (and sharper) the peak identifying the matrix maximum, the more reliable is the estimated shift.Note that C is a circular matrix, and therefore the maximum detectable shift is equal to ±d/2, where d is the patch dimension (usually square).
In order to identify reliable shifts, two quality parameters are considered, i.e., the peak value c max of the cross-correlation matrix and the ratio q = c max / C .For both of them, a pre-determined user-defined threshold is used to exclude invalid GCPs.
Accepted GCPs are then subjected to a filtering procedure to minimize noisy displacement patterns.To this end, shifts are classified based on their sign; specifically, a three-class classification map (positive shift, negative shift, and null shift) is generated.A connected component labeling algorithm [33] is used to segment this classification map.Small regions surrounded by a homogeneous background of shifts with the same sign are assumed to be noisy patterns and are discarded.Finally, filtered GCPs are interpolated into the final displacement map.

Results
An example of the results obtained with the SPOT technique is reported in Figure 3.It refers to the time frame August 2011-August 2012.As expected from the past literature, the highest velocity values have been recorded in the central part of the landslide (upper and lower necks).The arrows in Figure 3 represent the direction of the estimated vector field retrieved from the North-South and East-West component.All the results presented here have been obtained after re-projection of SAR images into a cartographic system using the NCALM LIDAR DEM.
The result for the time frame August 2012-August 2013 is very similar to that for the time frame August 2011-August 2012, and therefore is omitted for brevity.In both experiments, the size of the correlation window was set to 64 pixels, the cross-correlation threshold to 0.1, the threshold on the ratio q between the peak and the mean of the correlation matrix to 4, and the oversampling factor to 4. This means that the minimum retrievable displacement is on the order of 17 centimeters in the azimuth direction and 10 centimeters in the range direction.In Table 1, the parameters set to run the experiments are summarized.The third one (Run 3) concerns a consistency check of the retrieved displacements fields.
In Figure 4, the maps of quality parameters relative to the pair August 2011-August 2012 are shown.Figure 4a represents the map of the maximum correlation coefficient and Figure 4b is the q map, e.g., the ratio between the peak and the mean of the cross-correlation matrix.The noisy displacement patterns displayed in Figure 3 correspond to areas characterized by a low value of the considered quality parameters for the GCPs selection.
In Table 2, quantitative data about the estimated displacements rate are reported.These values were obtained by averaging the estimated velocities in the kinematic regions indicated in [22].Data concerning the standard deviation of the measurements are also provided.
Two different window dimensions (64 and 128 pixels) were experimented, and produced similar results.The last columns of the table are reserved for the consistency check, which was conceived as follows.A is the result for the period 2011-2012, B that for the period 2012-2013, and C that for the period 2011-2013; reliable results should return A + B − C = 0.This is not strictly achieved, however the deviation from zero of the aforementioned equation is generally small.The distribution of the consistency check is Gaussian with a mean of about 1 cm/year and standard deviation of 42 cm/year in the 64-pixel window case.In the 128-pixel window case, the mean and standard deviation of the distribution are on the order of 2 cm/year and 42 cm/year, respectively.deviation from zero of the aforementioned equation is generally small.The distribution of the consistency check is Gaussian with a mean of about 1 cm/year and standard deviation of 42 cm/year in the 64-pixel window case.In the 128-pixel window case, the mean and standard deviation of the distribution are on the order of 2 cm/year and 42 cm/year, respectively.

Validation
The Slumgullion landslide has been widely investigated in the past literature.Most of the studies rely on field surveys [1,[20][21][22][23]29,30,34].Recently, some papers exploiting remote sensing data have been presented [2,[25][26][27]31].In this context, only Reference [2] made use of space-borne SAR data.However, this work was focused on spotlight DInSAR methods covering one year of observations, with small insights in long-term displacement monitoring and limited validation of the presented results against literature data.All other works reviewed rely on airborne images acquired by the NASA/JPL L-band UAVSAR with 0.6-and 1.9-meter spatial resolution in the azimuth and slant range directions, respectively [35].
The validation of the obtained results is implemented both at landslide scale and at point scale.A perfectly consistent validation set (i.e., ground measurements acquired over the same time span covered by the SAR images used in this study) is not available, especially concerning the landslide scale.In this case, the most referenced data are relevant to the period 1985-1990 and to the year 2010.They are reported in Table 3 for the ease of the reader.
Data concerning the period 1985-1990 were produced in [29] and [30] using photogrammetry and field surveys.Data relevant to the year 2010 were collected in [22] through Ground-Based SAR Table 2. Summary of the obtained results.Data are reported in meters per year and have been averaged using the regions indicated in [22].V-displacement velocity; CC-consistency check.

Validation
The Slumgullion landslide has been widely investigated in the past literature.Most of the studies rely on field surveys [1,[20][21][22][23]29,30,34].Recently, some papers exploiting remote sensing data have been presented [2,[25][26][27]31].In this context, only Reference [2] made use of space-borne SAR data.However, this work was focused on spotlight DInSAR methods covering one year of observations, with small insights in long-term displacement monitoring and limited validation of the presented results against literature data.All other works reviewed rely on airborne images acquired by the NASA/JPL L-band UAVSAR with 0.6-and 1.9-m spatial resolution in the azimuth and slant range directions, respectively [35].
The validation of the obtained results is implemented both at landslide scale and at point scale.A perfectly consistent validation set (i.e., ground measurements acquired over the same time span covered by the SAR images used in this study) is not available, especially concerning the landslide scale.In this case, the most referenced data are relevant to the period 1985-1990 and to the year 2010.They are reported in Table 3 for the ease of the reader.
Data concerning the period 1985-1990 were produced in [29,30] using photogrammetry and field surveys.Data relevant to the year 2010 were collected in [22] through Ground-Based SAR Interferometry (GBInSAR) measurements.The latter study highlighted that, by 2010, the landslide's velocity halved compared to its values in the 1985-1990 period, and that the landslide head was affected by the largest decrease in velocity.The authors ascribed this behavior to both geomorphological and climatic factors.At the climatic level, they suggested that the average increase in temperature and decrease in precipitation could have induced an overall slowing of the landslide.Reference [1] showed that the movement of the Slumgullion landslide is strongly correlated to the soil moisture, and that its decreasing trend reflected the general slowing of the landslide.From a geomorphological point of view, the observed thinning of the landslide head [36] caused its stability to increase.The slowing of the head should have favored the overall slowing of the landslide by decreasing downslope-directed transfer of shear stresses [22].
The displacement rates retrieved in this study are congruent with those reported in Table 3 in the column relevant to the GBInSAR survey implemented in 2010 [22].The direction of the estimated vector field (see arrows in Figure 3) mainly follows the landslide slope profile and is qualitatively quite similar to that presented in the past literature (see as an example [31]).
The obtained results showed that the effect of the variation of the correlation window is negligible from the viewpoint of the estimated displacements, being for all the regions below the theoretical sensitivity of the method, which, as previously stated, is given by 1/f, where f is the applied oversampling factor.On the other hand, defining a smaller correlation window (e.g., 32 pixels) makes the frequency-domain cross-correlation less reliable, and this increases the standard deviation of the estimated displacement field (not reported here for brevity), which results in noisiness and physical inconsistency.Therefore, it is suggested to operate with the 64-pixel window.This allows a lower computational time compared to the 128-pixel window (for these experiments, the computational time was about 2.1 h per run, compared to about 45 min for each 64-pixel window run), as well as a higher level of detail and a better preservation of the landslide edges.
For the 64-pixel window, the registered values of the standard deviation range from 0.09 m/year in Region 1 (pair 2012-2013) to 0.67 m/year in Region 7 (pair 2011-2012).They are similar to those indicated in [29,30].
It is remarkable that the noisier displacement patterns, see as an example that north of landslide Region 8, are characterized by very low values of the quality parameters considered.This means that the peak of the correlation matrix is not sharp (i.e., it is not well-defined compared with the background), thus leading to an unreliable estimate of the displacements.In the landslide area, even though the peak of the correlation matrix c max is not very pronounced (as expected based on the characteristics of the phenomenon under investigation), its ratio q with respect to the mean is quite high.The average values registered for the pair 2011-2012 for these quantities within the defined kinematic regions are, respectively: Region 1-0.33, 12.34; Region 2-0.39, 12.65; Region 3-0.32, 12.14; Region 4-0.30, 12.33; Region 5-0.30, 12.38; Region 6-0.28, 10.21; Region 7-0.24, 9.70; Region 8-0.23, 9.92; Region 9-0.25, 11.17; Region 10-0.28,11.42; and Region 11-0.24,10.27.
The consistency check involves the three pairs: August 2011-August 2012; August 2012-August 2013; and August 2011-August 2013.A consistent result should pose that the sum of the displacements is zero.As stated this is not strictly achieved, however the resulting distribution has a mean very close to zero either at landslide scale or within the 11 kinematic regions.The registered deviations from zero (using the 64-pixel window) range (in absolute value) from less than 1 cm/year (region 2) to 12 cm/year (region 5).Similar results were obtained using the 128-pixel window (see Table 2).
In the following, the obtained results will be discussed at a finer scale exploiting data concerning 19 measurement points (MPs) installed on the landslide by the US Geological Survey (USGS) [34].Reference data for comparison were extracted from Reference [26], in which the kinematics of the landslide was analyzed in the time frame August 2011-April 2012 using airborne L-band remote sensing and was compared with GPS data collected at the USGS MPs.Note that these data were reproduced through graph digitalization, and it is therefore possible that they exhibit (negligible) differences with respect to their original version.
In Figure 5, the comparison between the results presented here (time frame August 2011-August 2012) and those from [26] is shown.The position of the MPs with respect to one of the available SAR acquisitions is also reported (note that MP1 and MP2 are not shown in this graphic since their position falls outside the image cut considered).Airborne SAR-derived data are rendered with blue bars and GPS data with gray bars.Outcomes from this study (orange bars) were produced by taking the maximum displacement in a window of 100 × 100 meters around each MP position (only for MP1, which is expected to be immobile, the mean displacement was considered).This choice was made in order to compensate geocoding errors (some of the MPs are at the landslide borders and/or at the transition between different kinematic zones) and SPOT maps inhomogeneity/noise.A consistent result should pose that the sum of the displacements is zero.As stated above, this is not strictly achieved, however the resulting distribution has a mean very close to zero either at landslide scale or within the 11 kinematic regions.The registered deviations from zero (using the 64-pixel window) range (in absolute value) from less than 1 cm/year (region 2) to 12 cm/year (region 5).Similar results were obtained using the 128-pixel window (see Table 2).
In the following, the obtained results will be discussed at a finer scale exploiting data concerning 19 measurement points (MPs) installed on the landslide by the US Geological Survey (USGS) [34].Reference data for comparison were extracted from Reference [26], in which the kinematics of the landslide was analyzed in the time frame August 2011-April 2012 using airborne L-band remote sensing and was compared with GPS data collected at the USGS MPs.Note that these data were reproduced through graph digitalization, and it is therefore possible that they exhibit (negligible) differences with respect to their original version.
In Figure 5, the comparison between the results presented here (time frame August 2011-August 2012) and those from [26] is shown.The position of the MPs with respect to one of the available SAR acquisitions is also reported (note that MP1 and MP2 are not shown in this graphic since their position falls outside the image cut considered).Airborne SAR-derived data are rendered with blue bars and GPS data with gray bars.Outcomes from this study (orange bars) were produced by taking the maximum displacement in a window of 100 × 100 meters around each MP position (only for MP1, which is expected to be immobile, the mean displacement was considered).This choice was made in order to compensate geocoding errors (some of the MPs are at the landslide borders and/or at the transition between different kinematic zones) and SPOT maps inhomogeneity/noise.Comparison between space-borne COSMO-SkyMed displacement rates (orange), airborne UAVSAR displacement rates (blue) and GPS displacement rates (gray) for the 19 USGS measurement points.UAVSAR and GPS displacement rates have been extracted from [26].
The three datasets qualitatively show a good agreement.Disagreements are concentrated, as expected, in the neck sector, where the landslide is faster.In this area, estimates from this work exhibit the highest discrepancy when compared to airborne SAR and GPS data, in particular in MP8, MP9, and MP14.For those points, the landslide displacement rate is underestimated when compared to Comparison between space-borne COSMO-SkyMed displacement rates (orange), airborne UAVSAR displacement rates (blue) and GPS displacement rates (gray) for the 19 USGS measurement points.UAVSAR and GPS displacement rates have been extracted from [26].
The three datasets qualitatively show a good agreement.Disagreements are concentrated, as expected, in the neck sector, where the landslide is faster.In this area, estimates from this work exhibit the highest discrepancy when compared to airborne SAR and GPS data, in particular in MP8, MP9, and MP14.For those points, the landslide displacement rate is underestimated when compared to GPS data of about 0.17 cm/day, about 0.3 cm/day, and about 0.16 cm/day, respectively.If the measurements extracted from airborne SAR images are considered, the underestimation in correspondence of the aforementioned points is on the order of 0.13 cm/day, 0.15 cm/day, and 0.17 cm/day.
Assuming that points MP1 to MP7 belong to the head, that MP8 to MP14 belong to the neck, and that MP15 to MP19 belong to the toe, the registered Root-Mean-Square Error (RMSE) of the displacement rates here estimated with respect to GPS measurements is about 0.05 cm/day for the head, 0.15 cm/day for the neck, and 0.09 cm/day for the toe.

Discussion
SPOT methods allow for the estimation of two-dimensional displacements by using only the intensity channel.They are less sensitive to atmospheric effects and land cover types (i.e. they can be applied also in presence of vegetation) if compared with DInSAR.This is due to the lack of the phase information, which also allows for avoiding unwrapping procedures, often leading to the failure of DInSAR because of the low density of valid pixels [37].
With a lack of any facilities (i.e., corner reflectors) installed in situ, the accuracy of the estimated displacements depends on the resolution of input images, the magnitude of the real displacements, the acquisition geometry, and the correlation coefficient between the scatterers [38].These parameters are clearly interconnected, since higher-resolution images can provide a better description of the scene features, thus maximizing their correlation at fine scale.On the other hand, when displacements are significant, the hypothesis of rigid shift is weakened, i.e., the correlation tends to be lower, as does the accuracy of the estimated offsets, as testified by the comparison with the GPS data provided in Section 3.1, which gave the worst results in correspondence with areas characterized by faster movements.
Spatial decorrelation can also be due to orbital issues.Differences in the acquisition geometry can significantly change speckle patterns, thus affecting the amplitude and the sharpness of the correlation peak [39].Therefore, it is crucial that images ingested in the SPOT algorithm are acquired with similar orbit, i.e., with contained normal baselines.Their values for the examined image pairs were on the order of 30 m for the couple 2011-2012 and 460 m for the couples 2012-2013 and 2011-2013, which values are significantly below the X-band critical baseline (i.e., the baseline value causing complete image decorrelation), which is more than 4 km [40].
The implemented SPOT method, validated against the available literature data both at landslide and at point scale, was able to successfully reconstruct the kinematics of the Slumgullion landslide using very high-resolution space-borne SAR remote sensing.At landslide scale, the estimated average displacement rates in relevant kinematic were consistent with the closest available data, which were acquired a couple of years before the first considered observation period.The obtained displacement maps resulted quite homogeneous in those regions, as testified by the low values of the standard deviation reported in Table 2.As expected, its peaks were registered where the landslide is faster, i.e., where both the correlation value and the q parameter exhibited the lowest value (see Figure 4).
At the point scale, the performance of the implemented SPOT was fully comparable to that achieved by dedicated airborne campaigns, whose results showed a slightly better agreement with reference GPS data.However, this was expected, since UAVSAR L-band images exploited in [26], having higher penetration depth, are less sensitive to land cover (e.g., vegetation) compared to the X-band COSMO-SkyMed data used for this study.This means that the correlation between different images is expected to be higher in the first case [41] and, in turn, the displacement estimation is more precise.
However, the acquisition of space-borne images is surely cheaper than that of airborne ones and, since differences in performance are small, it can be argued that the first solution is more beneficial, being also independent from weather conditions, especially when monitoring is implemented with high temporal frequency.

Conclusions
Satellite technologies have widely demonstrated their potential in environmental monitoring and forecasting, response, and recovery of natural hazards.In this paper, the suitability of space-borne SAR observations for fast/slow landslide monitoring was explored.This application cannot be addressed using classic differential SAR interferometry due to the magnitude of the displacement rate and scarce preservation of interferometric coherence when movements are on the order of meters.Therefore, sub-pixel offset tracking methodology was exploited.This technique, using only the amplitude channel, is less sensitive to atmospheric and scene decorrelation factors and does not have particular restrictions on the observable displacement rates.Thus, it can be applied to fast/slow landslides such as the Slumgullion landslide (Colorado, US) which represented the objective of this study.
This landslide has been extensively investigated in the past through field surveys.Recently, airborne remote sensing has been exploited for monitoring displacements, however, as far as we know, there have been no results obtained using space-borne images validated against ground data.Therefore, state-of-the-art sub-pixel offset tracking was applied to pairs of SAR images acquired with about one-year temporal baseline in spotlight modality with sub-meter resolution by the COSMO-SkyMed constellation.These were combined in two pairs covering two years, from August 2011 to August 2013, and the obtained results were validated with an inter-sensor consistency check and against the temporally closest available ground data.
First, the cross-consistency of the estimated displacement fields was checked by comparing the movements of the period 2011-2013 with the sum of those estimated for the time frames 2011-2012 and 2012-2013.The comparison was satisfactory, since the resulting displacement rates showed discrepancies on the order of a few centimeters per year.
The obtained displacement fields were then compared with available literature data, both at landslide scale and at point scale.
At landslide scale, the average displacement rates computed via ground-based SAR interferometry in the aforementioned kinematic regions identified in [18] were used as benchmark.As a result, the displacement fields estimated with satellite imagery were fully comparable with those provided by field surveys.
At point scale, literature data regarding 19 measurement points installed by the US Geological Survey were exploited.For that location, both GPS data and airborne remote sensing derived data were available.The comparison was satisfactory even in this case, since maximum deviations from literature data were on the order of fractions of a centimeter localized in the central part of the landslide.
The results discussed in this study demonstrated the reliability of space-borne radar imagery for large landslide-induced movement monitoring and its consistency with the results given by more expensive and time-demanding methodologies such as field surveys and airborne remote sensing.Future research will address the validation of the displacements estimated using this methodology over an extended time frame.

Figure 1 .
Figure 1.Study area.The Slumgullion landslide is represented by the area within the purple line.In box (a), the LIDAR Digital Elevation Model (DEM) by the US National Center for Airborne Laser Mapping (NCALM) has been superimposed onto the Google Earth view of the landslide.

Figure 1 .
Figure 1.Study area.The Slumgullion landslide is represented by the area within the purple line.In box (a), the LIDAR Digital Elevation Model (DEM) by the US National Center for Airborne Laser Mapping (NCALM) has been superimposed onto the Google Earth view of the landslide.

Figure 2 .
Figure 2. Flow diagram of the implemented methodology.The upper part of the picture shows the general processing chain, while in the lower part the flux represents an exploded view of the Sub-Pixel Offset Tracking (SPOT) processing block.

Figure 2 .
Figure 2. Flow diagram of the implemented methodology.The upper part of the picture shows the general processing chain, while in the lower part the flux represents an exploded view of the Sub-Pixel Offset Tracking (SPOT) processing block.

Figure 3 .
Figure 3. Slumgullion landslide displacement rate, in meters per year, for the time frame August 2011-August 2012.The superimposed arrows represent the direction of the estimated vector field.

Figure 3 .
Figure 3. Slumgullion landslide displacement rate, in meters per year, for the time frame August 2011-August 2012.The superimposed arrows represent the direction of the estimated vector field.

Figure 4 .
Figure 4. Ground Control Point (GCP) quality parameter maps for the satellite pair August 2011-August 2012.(a) Maximum of the cross-correlation matrix.(b) Ratio between the peak and the mean of the cross-correlation matrix.In both pictures, black contours indicate the kinematic regions as defined in [22].

Figure 4 .
Figure 4. Ground Control Point (GCP) quality parameter maps for the satellite pair August 2011-August 2012.(a) Maximum of the cross-correlation matrix.(b) Ratio between the peak and the mean of the cross-correlation matrix.In both pictures, black contours indicate the kinematic regions as defined [22].

Figure 5 .
Figure 5.Comparison between space-borne COSMO-SkyMed displacement rates (orange), airborne UAVSAR displacement rates (blue) and GPS displacement rates (gray) for the 19 USGS measurement points.UAVSAR and GPS displacement rates have been extracted from[26].

Figure 5 .
Figure 5.Comparison between space-borne COSMO-SkyMed displacement rates (orange), airborne UAVSAR displacement rates (blue) and GPS displacement rates (gray) for the 19 USGS measurement points.UAVSAR and GPS displacement rates have been extracted from[26].

Table 1 .
Summary of the adopted Sub-Pixel Offset Tracking (SPOT) parameters.

Table 1 .
Summary of the adopted Sub-Pixel Offset Tracking (SPOT) parameters.

Table 2 .
[22]ary of the obtained results.Data are reported in meters per year and have been averaged using the regions indicated in[22].V-displacement velocity; CC-consistency check.