Slope-Scale Evolution Categorization of Deep-Seated Slope Deformation Phenomena with Sentinel-1 Data

Abstract

Deep-seated gravitational slope deformations (DsGSDs) are slope-scale phenomena which are widespread in mountainous regions. Despite interacting with human infrastructures and settlements, only a few cases are monitored with ground-based systems. Remote sensing technologies have recently become a consolidated instrument for monitoring and studying such widespread and slow processes. This paper proposes a three-step novel methodology to analyze the morpho-structural domain of DsGSDs by exploiting the advanced Differential Synthetic Aperture Radar Interferometry (A-DInSAR) technique through (i) the analysis of A-DInSAR measurement point density and distribution defining a coverage threshold; (ii) the assessment of the actual ground deformation with respect to the orientation of phenomena based on slope, aspect, and C-index; and (iii) ground deformation mapping with previously ranked velocity interpolation. The methodology was tested on two differently oriented phenomena: the mainly north–south-oriented Croix de Fana and the mainly east–west-oriented Valtournenche DsGSD, located in the Aosta Valley Region, northern Italy. The results show a variation in the kinematic behavior between the morpho-structural domains, while also considering any other superimposed surficial deformations. This work provides the lines for the implementation of a rapid and low-cost tool based on the use of A-DInSAR measurements which are suitable for assessing the impact of any type of DsGSD on the anthropic facilities and infrastructures in mountainous areas.

1. Introduction

Deep-seated gravitational slope deformations (DsGSDs) widely affect mountainous landscapes worldwide [1,2,3]. Among the natural hazards, these large slow-moving phenomena mainly present a not impulsive evolution, but extremely continuous over a long period, posing a certain risk for human facilities [4,5,6].

These huge phenomena extend for several kilometers in length and hundreds of meters in depth, involving entire valley flanks, dislocating large rock volumes through a continuous slow creeping displacement, and potentially evolving in local collapses [7]. Because of the giant extension and the prolonged continuous deformation over thousands of years [7,8], the evolution of DsGSDs implies a combination of different controlling factors with variable interactions [3]. Distinct geomorphological signs, mainly represented by persistent linear morpho-structural elements, are associated with DsGSD evolution. Generally, double or multiple ridges, scarps and counterscarps, open or infilled trenches, downthrown blocks, graben, and ridge top depressions can be recognized in the upper slope portion, suggesting an extensional deformation regime. Conversely, a compressive deformation style characterizes the middle-lower slope sector, highlighted by toe bulging, buckling folds, intensely fractured rock masses and often secondary nested landslides. These characteristics involve a high spatial and temporal heterogeneity, with strain partitioning in distinct morpho-structural domains, with different modes of deformation and rates of displacement [8,9,10].

The heterogeneity in diverse deformation regimes, combined with the huge volumes involved, make the characterization of DsGSD mechanisms and their activity and behavior a challenge, mainly in terms of impact assessment on anthropic elements and of risk reduction policy definition. Recently, in the literature, diverse approaches and methodologies have been proposed for a proper characterization of these complex phenomena. Advanced remote sensing technologies, spanning from Lidar-based geomorphic analysis [11,12] to Differential Synthetic Aperture Radar Interferometry (DInSAR), have been proposed to spatially and temporally characterize the DsGSDs behavior [13], their impact on structures and infrastructures [5,14], and their state of activity [15,16], which has proven to be extremely effective for DsGSD investigation. However, a comprehensive slope-scale characterization of these massive phenomena to assess the variable displacement pattern and their complexity, distinguishing the distinct kinematic domains, remains an open issue. Advanced DInSAR (A-DInSAR) approaches are consolidated tools to investigate the slow-moving phenomena that affect slopes [17,18], providing the possibility to investigate ground deformation over wide areas with sub-centimeter accuracy [19,20]. Nevertheless, their application in mountainous areas remains challenging due to some limitations arising from the spatial characteristics of the sensor datasets, the residing meteorological condition, and the snow and vegetation cover, which may cause a high reduction in measured points coverage and an uneven distribution, which is inadequate for a comprehensive spatial analysis of displacement patterns [21,22]. High-resolution A-DInSAR analysis [23] has been largely exploited for landslide risk assessment, allowing for the monitoring of Persistent Scatterers (PSs), i.e., targets that are characterized by an electromagnetic response, which remains coherent over observation periods that even span years [24]. Moreover, recent advancements within the high-resolution A-DInSAR processing context have highlighted the possibility of increasing the measurement point density, thus improving the monitoring capabilities. PSs and Distributed Scatterers (DSs), i.e., targets undergoing decorrelation phenomena typically located in rural/sparsely vegetated areas, can be also selected in a more general multi-resolution framework. A possibility for doing so is offered by the SqueeSAR [25] or the Component extrAction and sElection SAR CAESAR [26,27,28] methods, which smartly fit the heterogeneous scattering response of the imaged ground scene. CAESAR [23], used in this work, allows for the preservation of PS information in areas where a high density of PSs is expected, at the same time thickening the measurement points where DSs are typically expected to be by exploiting a spatially variable multilook degree. In addition, it is important to note that the Line of Sight (LOS) of the sensor is, generally, almost perpendicular to the north–south direction, and SAR interferometric data stacks are acquired over ascending and/or descending orbit typologies. This leads to the LOS being directed mainly toward the east or west, respectively. Consequently, the system has an intrinsically low sensitivity to north–south deformation components. A procedure that assumes a prevalent translational movement along the steepest slope direction for each PS [29,30,31] should be considered to partially overcome this problem.

In this paper, a multi-step methodology based on Sentinel-1 images, processed by using the CAESAR detector (CAESAR-D) [32], for a slope-scale characterization of the morpho-structural domains of DsGSDs has been implemented. The methodology was applied to two phenomena located in the Aosta Valley Region (north-western Italy) characterized by opposite orientation sensitivity, i.e., Croix de Fana DsGSD, mainly north–south-oriented, and Valtournenche DsGSD, mainly east–west-oriented. In order to assess the slope spatial evolution, a PSs coverage analysis of the distribution of measurement points for an appropriate characterization of DsGSD domains behavior was carried out. Subsequently, to investigate the sensitivity of the LOS displacements in slope deformation measurements from the real slope displacement, a SAR data suitability analysis based on the matrix analysis of the measurement points jointly considering slope and aspect, and satellite geometries, was implemented in order to obtain ranked ground deformation maps derived from fileted PSs interpolation.

The proposed methodology offers a novel approach for analyzing and defining the complex behaviors of deep-seated natural phenomena and pinpointing crucial and vulnerable sectors, contributing to assessing geohazards and managing the associated risk in mountainous regions affected by these huge phenomena.

2. Materials and Methods

To characterize the DsGSD evolution and analyze their morpho-structural domains, a methodology based on the exploitation of measurement points retrieved by the multiresolution CAESAR-D processing [32] of Sentinel-1 datasets acquired on ascending and descending orbits over the Aosta Valley Region (north Italy) was developed. The implemented approach is structured in three main steps: (i) SAR data coverage analysis; (ii) SAR data suitability analysis; and (iii) SAR data interpolation. The entire workflow is schematized in Figure 1. To establish the capability of this procedure, the methodology was tested on two slope instabilities characterized by orthogonal orientation to assess the degree of the variability of the slope and aspect of each domain and evaluate the suitability rate of the SAR datasets in the slope surface deformations analysis.

2.1. DsGSD Morpho-Structural Domains

The DsGSD phenomena considered in this work are located in the north-western part of the Italian Alps, in the Aosta Valley Region (Figure 2), i.e., the Croix the Fana DsGSD, located in the middle portion of the main valley, and the Valtournenche DsGSD, located in the namesake lateral valley. These case studies are representative of this mountainous territory, highly affected by DsGSDs (about 13% of the regional territory), and are two of the most extensive active deep-seated phenomena with distinctive easily recognizable morpho-structural domains [9,33].

The Croix the Fana DsGSD covers an area of about 13 km², located on the north side of the east–west-oriented main valley. The phenomenon ranges between 525 m a.s.l. along the Dora Baltea River and 2211 m a.s.l. in correspondence with the Croix-de-Fana peak, involving part of the Quart municipality and affecting a series of strategic infrastructures such as roads, rail networks, and the penstock of the hydroelectric plant of Quart (Figure 2b). From a geological point of view, the bedrock outcropping is mainly represented by lithotypes belonging to the Piedmont Domain, here represented by prasinites and calcschists with the intercalation of mica marbles and metabasites, and to the Austroalpine System represented by orthogneiss, fine-grained gneisses, and subordinate amphibolite lenses. The whole slope has a glacial footprint superimposed with watercourse activities and the gravitative landforms. The DsGSD phenomenon, mainly north–south-oriented, apart from the easternmost sector incised by the Saint Barthélemy torrent, displays a mean slope gradient ranging from 20° up to 30°. An extended scarp of about 3.2 km in the western portion of the DsGSD, variably oriented N40°–N10°, associated with a series of counter scarps forming, forms a graben structure. This setting displays an upper extensional zone (morpho-structural domain A) [35], extended by the Croix-de-Fana peak up to the Avisod and Fonteil hamlets. A series of other typical gravity-induced features, such as trenches, scarps, and elongated depressions, are identifiable along the slope, mainly distributed in the B, D, and E domains, sectors already distinct in the Italian Landslide Inventory IFFI [34,36] based on geomorphological criteria. The lower portion (domain C) is characterized by a compressional kinematic, evidenced by highly fractured rock masses, locally involving and overlapping the glacial deposits, pointing out the complex evolution of this phenomenon [37]. Moreover, evidence of rock dissolution is testified by the presence of local depression and sinkholes [38]. In domain F, mainly west-oriented, it is important to note the complex landslide located in the eastern portion of DsGSDs, i.e., the Vollein landslide, currently monitored in the regional early warning system network [39].

The Valtournenche DsGSD extends for about 4 km² in the middle portion of the Valtournenche valley, affecting the namesake village up to the Paquier locality. The phenomenon ranges from 300 m a.s.l. along the valley bottom up to about 2800 m a.s.l., with a mean slope gradient from about 15° to 30°. From a geo-lithological point of view, the area affected by the DsGSD presents a series of mafic rocks as gabbroid rocks, and lithotypes in the greenschist facies as serpentinites, micaschscists, calcschists, and prasinites pertaining to the Piedmont Domain. The slope has been reshaped by glacial and gravitative processes that have hidden or revealed the basement rocks. The upstream sector (domain A) corresponds to a steep slope sector superimposed by two east–west-extended landslides. Along the watershed, extended debris coverages are also observable. The middle portion corresponds to a gentle slope sector (domain B), partially affected by landslide deposit accumulation. A noticeable secondary scarp, north–south-oriented, separates this sector from the lower one (domain C), which affects the Paquier hamlet. The downstream sector (domain D) displays a gentle slope corresponding to the most urbanized area, with a local inactive landslide (L3), catalogued in the IFFI, which is recognizable in the northernmost sector.

2.2. Sentinel-1 Datasets and A-DInSAR Processing Methodology

The SAR datasets consist of 336 scenes along the ascending orbit (track 88) and 278 scenes for the descending orbit (track 139), acquired using Sentinel-1 A and B satellite sensors in the C-band, over an area of the Western Italian Alps. The acquired images cover the period from late 2014 (i.e., October for the descending dataset and November for the ascending one) to June 2021. Both datasets were processed with a two-step A-DInSAR processing procedure. The first step, implemented with a Small BAseline Subset (SBAS) approach [40,41], allows for the compensation of the atmospheric phase screen and the estimation of small-scale (low-resolution) displacement. The second step implements a model-based detection of PSs/DSs via the multi-resolution component extraction and selection SAR (CAESAR) detector (CAESAR-D) [32]. This detector is able to increase the density of measurement points via a spatially variable multilook. It allows preserving Permanent Scatter (PS) information in areas where a high density of PSs is expected, at the same time improving the density of measurement points where Distributed Scatters (DSs) are typically expected to be. In this way, for each measurement point in the study areas, the Line of Sight (LOS) mean deformation velocity (V_LOS) relative to the period of observation (i.e., 2014–2021) is obtained. Figure 3 shows the distribution of the measurement points with the LOS velocities for the ascending (V_LOSa) and descending (V_LOSd) geometries.

The ascending and descending data were also exploited to project the LOS velocities along the steepest slope (V_SLOPE) based on the C-index [29,31] application using Equation (1):

V_SLOPE = V_LOS/C

(1)

where C is given by:

C = {N × [COS(s) × SIN (a − 90)]} + {E × [−1 × COS(s) × COS(a − 90)]} + [H × SIN(s)]

(2)

where N, E, and H are the directional cosines of the radar LOS, and a and s are aspect and slope gradients, respectively. For V_SLOPE computation, to avoid any exaggeration, all the coherent pixels characterized by a C value between 0.2 and −0.2 were removed, as thoroughly explained in Section 2.3. Additionally, only measurement points with a slope gradient above 5° have been projected, considering that low slope degree may be affected by other processes unrelated to slope movement, such as subsidence or local deformation; consequently, their projection will be biased. For C-index computation, slope and aspect values were derived from the STRM 30 m digital terrain model (DTM) to reduce and smooth out the slope morphology variations.

Moreover, ascending and descending datasets allowed us to compute the E-W and vertical components of LOS velocity using Equations (3) and (4), respectively:

V_ew = (V_LOSd/h_d − V_LOSa/h_a)/(e_d/h_d − e_a/h_a)

(3)

V_v = (V_LOSd/e_d − V_LOSa/e_a)/(h_d/e_d − h_a/e_a)

(4)

where h_a, h_d, e_a, and e_d correspond to the directional cosines of ascending (a) and descending (d) LOS vectors derived from the incident angle and LOS azimuth of the satellite used.

2.3. SAR Data Coverage Analysis

The PS spatial distribution plays a key role in ground deformation analysis, particularly for the spatial characterization of extensive deep-seated phenomena such as DsGSDs. This poses a challenge in DInSAR technique application in mountainous areas, mainly related to complex topography with variable high slope gradients, land use coverage characterized by widespread forested areas, and snow cover for a long period of the year. Hence, in the first phase of the implemented methodology, a spatial analysis of the measurement points has been executed at the morpho-structural domain scale. The SAR data coverage analysis aims to assess the suitable distribution and density of measurement points, excluding areas with an uneven distribution of points, to ensure a proper characterization of ground deformation in a reliable way. Firstly, spatial data coverage was computed by dividing the areas of interest into a 100 × 100 m grid. Afterwards, the density of the measurement points related to all available datasets internally to each cell was calculated. In accordance with [42,43], all the cells with a density value lower than 30 PSs/km² were discarded. The imposed density threshold allows us to better estimate the slope deformation, computing the PS Coverage (C_PS) for each domain with the following equation:

C_PS = V_c/T_c

(5)

where V_c are the valid cells with density value ≥30 PS/km², and T_c is the total number of the cells within each area of interest. An empirical threshold corresponding to a coverage value of 25%, under which the result may be considered unreliable, was fixed. First, we exploited the V_SLOPE dataset for the data coverage analysis, being a derivable datum from all cases of slope orientation, as well as being the dataset with the highest measurement point abundance. The same analysis was also carried out for the other datasets available to verify whether each sector reached a C_PS > 25% (Supplementary Material, Table S1).

2.4. SAR Data Suitability Analysis

The adoption of DInSAR techniques allows us to measure the average velocity along the LOS of the satellite only. Thus, with the LOS mainly belonging to the plane containing the vertical and east–west directions, an intrinsic low sensitivity to north–south deformation components occurs. Hence, the capability to detect ground surface deformation is notoriously affected by some limitations, mainly related to the slope topography setting (i.e., slope and aspect) and the acquisition geometries of the considered satellite (i.e., incident angle θ and azimuth φ) [21,29,44].

To evaluate the LOS displacement sensitivity compared to actual slope displacement, a Dataset Suitability Ranking (DSR) was computed. Firstly, the matrixes of C-index values for both ascending (Ca) (Figure 4a) and descending (Cd) (Figure 4b) geometries were computed, considering every possible combination of aspect and slope. The value represented in each cell of the matrixes Ca and Cd is the C index that was calculated with Equation (2) using the slope gradient (y-axis), the aspect (x-axis), and the LOS cosine for both ascending (Ca, Figure 4a) and descending (Cd, Figure 4b) geometries. We grouped the C index values into 5 classes. In particular, the classes with C-index values over 0.2 (green, cyan, and blue cells) represent the V_LOS progressively near to the ground truth velocity, the class with the range −0.2–0.2 (red cells) represents the slopes with the geometry combination settings that cause over-estimated V_LOS projection, and finally, the C-index values below −0.2 (grey cells) represent the slopes with reverse values of V_LOS (e.g., positive descending V_LOS on slope toward east). Then, combining the Ca and Cd matrixes, a third matrix representing the classified DSR values was obtained (Figure 4c) using the parameters defined in

Table 1. Combination of the possible cases identifying the relative qualitative velocity ranking associated with the most suitable type of velocity component.

Aspect Class	Slope Threshold	C-Index Combination	Case	Qualitative Ranking	Velocity Type Priority
45–135°	<30°	Ca > 0.2 and Cd < −0.2	Case 1a	High	Vew–Vv; VLOSa; VSLOPE, VLOSd
225–315°	<30°	Ca < −0.2 and Cd > 0.2	Case 1d	High	Vew–Vv; VLOSd; VSLOPE, VLOSa
45–135°	>30°	Ca > 0.2	Case 2a	Medium–High	VLOSa; VSLOPE
225–315°	>30°	Cd > 0.2	Case 2d	Medium–High	VLOSd; VSLOPE
(135–225°) or (0–45°) or (315–360°)	-	Ca > 0.2 or Cd > 0.2	Case 3	Low	VSLOPE; Vv; Vew
-	-	−0.2 ≤ Ca ≥ 0.2 and −0.2 ≤ Cd ≥ 0.2	Case 4	Worst case	None

. Through a step-by-step conditional procedure, the DSR matrix is compiled. Each class of DSR is defined with (i) a combination of Ca and Cd thresholded values; (ii) the slope orientation, predefined aspect classes, N, S, E and W; and (iii) a slope gradient threshold (where necessary) that depends on the LOS incidence angle, which is based on the maximum slope values for the class C < −0.2.

The classified DSR suggests the most suitable velocity type to use for the ground deformation analysis. Therefore, once the average orientation of the phenomenon of interest is known, it is possible to check the suitability of the dataset.

According to Figure 4c and Table 1, four main cases based on the combinations of slope and aspect were defined, with their relative quality rankings.

Case 1a–d: In this case, both Ca and Cd have an absolute value > 0.2 and the slope is oriented toward the east (±45°) for sub-case 1a or the west (±45°) for sub-case 1d. The slope gradient is <30° (for the used Sentinel-1 incidence angle). This case has a high-quality ranking because it has the best orientation and allows for the use of all velocity datasets. When both geometries are available, the resolved E-W and vertical velocities (V_ew, V_v) are the most suitable dataset, followed by ascending V_LOSa for Case 1a and V_LOSd for Case 1d. The V_SLOPE also can be used as a reliable dataset.

Case 2a–d: In this case, the slope orientation is toward the east (±45°) and Ca > 0.2 (Case 1a), or west (±45°) and Cd > 0.2 (Case 2d). For the used Sentinel-1 incidence angle, the slope gradient is >30°. This case has a medium-high quality ranking and allows for the use of the ascending dataset (V_LOSa) for the eastern slope and descending dataset (V_LOSd) for the western slope as the most reliable results, followed by the V_SLOPE that still has reliable values. The resolved velocities (V_ew, V_v) are usually difficult to compute because of the high slope gradient. On such steep slopes, one of the two geometries is affected by well-known topographic effects (e.g., foreshortening, layover) that decrease the density of measurement points.

Case 3: In this case, at least one of the geometries should have C > (|0.2|) and the slope orientation toward the south or north (±45°). The slope gradient is not a constraint. This case has a low-quality ranking, hence the V_SLOPE or V_v (when both ascending and descending velocities are available) allow us to obtain a reliable velocity. The east–west component can also be computed, considering that it strongly underestimates displacement.

Case 4: In this case, both datasets have a C < |0.2|. This condition is the worst case because all datasets are affected by low reliability. However, this is limited to a small combination of slope/aspect: it corresponds to slope-oriented N-S with a low slope gradient.

Following the schemas shown in Figure 4 and Table 1, the DSR of each polygon representing a sector of DSGSD or a landside was computed using the field calculator in the GIS environment.

It is important to note that, since the LOS velocity projection along the slope is resolved, assuming that the displacement occurs along the maximum slope direction, it could be affected by some biases. Instead, the V_ew and V_v components are resolved without any additional hypothesis on the basis of a vectorial composition of ascending and descending V_LOS. Hence, for Case 3, the V_SLOPE is an estimated velocity that probably corresponds to ground truth. At the same time, the vertical component, even if it underestimates the total displacement, is not affected by a priori assumptions or bias.

2.5. SAR Data Interpolation

A spatial analysis was performed, interpolating the filtered measurement points by exploiting the results obtained in the previous stages. Following the priority obtained with DSR, each velocity type was interpolated. Using an open-source GIS environment, i.e., QGIS, we interpolated the diverse velocity typologies with the Inverse Distance Weighting (IDW) tool. The IDW technique is particularly well suited for constructing ground deformation maps through the exploitation of SAR data because the tool operates on the premise that neighboring PSs exhibit stronger correlations compared to distant measuring points [45,46]. The interpolation radius was set to 100 m, while the minimum number of points was set to 2. These values were chosen to be coherent with the previous density threshold of 30 PS/km². To obtain smoother results, we also empirically set the ‘weighting power’ parameter to 2, and the ‘smoothing’ one to 10. Moreover, using the SAGA (System for Automated Geoscientific Analyses) library available on QGIS, we applied the ‘Simple filter’ tool by setting the radius to 10 to further smooth the results obtained via the IDW interpolation process. Smoothing operations serve to obtain a more readable spatial distribution of deformations and the parameters controlling the intensity of smoothing must be set according to the degree of detail desired. It should be noted that over-smoothing means a loss of data. Depending on the cases previously identified (Section 2.4), the interpolation was performed following the priority given by the DSR qualitative ranking and considering only the inherent velocity types.

3. Results

The proposed methodology was applied and tested on the two orthogonal-oriented phenomena, i.e., Croix de Fana, north–south-oriented, and Valtournenche, east–west-oriented, selected as representatives of the best and weakest cases in terms of the orientation of the investigated phenomena.

3.1. PS Density Map

The SAR data coverage analysis allowed us to assess the finest distribution of the measurement points, focusing on each domain of the two considered phenomena. Figure 5 shows the percentage of the coverage (C_PS) defined for each domain, considering as valid the cells with PS density equal to or greater than 30 PS/km², obtained in the first stage of the implemented approach.

The Croix de Fana DsGSD domains (Figure 5a,c) show good coverage, with PS percentage coverage values computed for all the datasets above 37% for all the morpho-structural domains. The highest coverage values match the most urbanized areas, variably distributed along the slope with widespread hamlets and locality and extended natural areas with little or no vegetation. Quite good coverage has been observed for the Valtournenche DsGSD (Figure 5b,d). Most of the domains present a C_PS greater than 40%, apart from domain B, which has poor coverage on all datasets, and domain A, which shows poor PS coverage in the V_LOSa, V_ew, and V_v datasets. The lack of data in domain B is mainly due to a large, forested area that limits the number of coherent pixels. Instead, the high density of the lower portion of the DsGSD (domains C and D) was observed, corresponding to a densely urbanized area. At the same time, a cluster with high PS density is localized between sectors A and B, related to unvegetated landslide accumulation (L1) and to vast debris coverage.

The C_PS values computed for the other velocity datasets (V_ew, V_v, V_LOSa, V_LOSd) are reported in the supplementary material (Table S1, Figure S1).

3.2. Data Suitability Ranking

Through the joint analysis of slope/aspect and C-index parameters, combined with the matrix analysis of the ascending and descending orbit, the DSR was determined for each domain (Figure 6).

The Croix de Fana phenomenon shows a variable degree in slope values, with an average slope ranging from 20 to 25° in the lower-central portions, i.e., B, C, and F, 25–30° in the upper domain A, up to 30–35° in the minor central domains D and E. Analyzing the aspect variability locally, the entire DsGSD is generally oriented toward the south, apart from the domain F, which is mainly oriented toward the south-east. Based on the sectors’ orientation toward the south (±45°) and the moderate slope gradient, all the domains of the Croix de Fana DsGSD were classified with low DSR (Case 3), corresponding to the V_SLOPE and V_v component as the most suitable velocity type to be considered for reliable ground deformation analysis.

For the Valtournenche phenomenon, a slightly variable slope gradient was observed, with a steeper sector corresponding to the upstream portion of the DsGSD, domain A, with slope values of 20–25°, against the slope values of 15–20° of the other domains. The orientation of the DsGSD is generally toward the west, with a slight variation toward the north-west in the middle slope sector. Considering the general orientation toward the west and the low–moderate gradient slope, it was possible to classify all domains in the high DSR class (Case 1d). Therefore, excluding domain B that does not meet the minimum density threshold imposed, for the other domains of this DsGSD, the V_LOS and the resolved V_ew, V_v, and V_SLOPE can be jointly considered and analyzed to assess the actual ground slope deformation and investigate the evolution of the phenomena.

3.3. Ranked Ground Deformation Maps

Following the criteria and the thresholds imposed in the SAR data coverage and SAR data suitability analysis, a series of ranked ground deformation maps were obtained for each considered phenomenon, using the IDW interpolation tool.

Based on the DSR classes, the highest sensitivity of ground deformation measurement for the main north–south-oriented Croix de Fana DsGSD can be interred from V_v and then V_SLOPE datasets (Figure 7). The interpolation of the V_v dataset (Figure 7a) displays a prevalent downward component, with a general decrease from the upper to the lower portion of the DsGSD. Specifically, the highest deformation occurs in the upper part (domain A), mainly in correspondence of the extended graben structure, with a downward movement of −4.44 mm/yr. This trend gradually decreases from the medium portion (domains B-D-E) of about −3 mm/yr, to the lower portions (domains C-F) with movements close to the stability range. This pattern was confirmed by the ground deformation map derived from the interpolation of the V_SLOPE (Figure 7b), with higher mean velocity in domain A (−11.7 mm/yr), in the range of −7–−9 mm/yr in the B-F-D domains to a lower deformation (about −3 mm/yr) in the C and F domains. The V_SLOPE shows a high value compared to V_v, while the E-W component (not reported in the figure) shows values close to 0, which are coherent with the slope geometry. Thus, the high values of V_SLOPE are likely related to a high north–south displacement component. The effect of the secondary landslide of Vollein (L_v), currently monitored as early warning in the regional monitoring network [39], displays a velocity range of V_SLOPE of −5–−10 mm/yr. Note a relevant component toward the east around 6 mm/yr in correspondence to the L_v.

For the main east–west-oriented Valtournenche DsGSD, the highest sensitivity of ground deformation measurement can be inferred primarily from V_LOSd jointly to the V_v, V_ew, and V_SLOPE datasets (Figure 8). Referring to the SAR data coverage analysis, the results obtained for domain B were unreliable due to poor coverage (lower than 25%). Consequently, these points were discarded in the interpolation procedure.

All the deformation maps obtained for the DsGSD of Valtournenche show a similar deformation pattern. Specifically, the maximum values of the LOS velocity related to the descending orbit are concentrated in domain C, near Paquier village, represented by a mean velocity value of −8.03 mm/yr, and in domain A, where the mean value of deformation is −5.19 mm/yr. Generally, it is possible to appreciate that the lower sectors of deformation (domains C and D) show widespread movement, with a slight gradient, while the upper part (domain A) shows a more irregular pattern. The V_SLOPE dataset shows values slightly higher than V_LOSd, which agrees with the high values of Cd on this west-oriented slope. It is worth observing, in the map of V_SLOPE, a large cluster of high movement along the slope in correspondence with the most northerly landslide body (L1), with a mean of V_SLOPE of about −20 mm/yr. The map of resolved velocities (Figure 8c,d) was also created for sectors C and D, excluding A and B domains that do not reach the C_PS threshold. In the downstream sectors, moderate values of the downward component of the LOS velocity were detected, with mean velocity values ranging between −2 mm/yr (domain D) and −5 mm/yr (domain C). It is possible to note (Figure 8d) the presence of a predominant V_ew component toward the west, compared to the vertical one with mean values ranging from −5 mm/yr (domain C) to −11 mm/yr (domain D). Moreover, the vectorial sum of the V_ew and the V_v shows values near the V_SLOPE, suggesting that the north–south component for this case is not significant and the datasets agree with slope geometries.

4. Discussion

An in-depth analysis of giant deep-seated slow-moving phenomena able to characterize their complex evolution and heterogeneous displacement pattern points out the need for a slope-scale investigation to analyze, in detail, the morpho-structural domains of each DsGSD and identify the most potential critical sector impacting the functionality of anthropic structures and infrastructures.

Due to the compound behavior related to multiple causal factors and the long-term evolution of these phenomena, the characterization of DsGSD deformation represents a challenging issue. Considering the extremely slow or very slow deformation trend, on-site monitoring is rarely used, apart from specific cases of point-wise ground-based measurements along sensitive infrastructures [47,48]. In such a context, the application of spaceborne DInSAR techniques proves to be a sustainable solution to remotely monitor and characterize these slow-moving phenomena, allowing us to measure ground surface displacements over large areas with millimeter precision [20,44,49]. However, the location of these phenomena in high mountainous areas entails well-known difficulties in effective DInSAR application. Local terrain complexities and high slope gradients, together with the acquisition direction of the considered satellite, could cause the geometric distortions (e.g., foreshortening, layover, shadowing) of SAR images. Moreover, extended forested areas and snow coverage for long periods of the year represent a constraint in SAR measurement point detection and its uniform distribution in the alpine territories. In the literature [5,42,50,51], it is stated that the number, distribution, and density of SAR data impact the correct assessment of the local-scale kinematic behavior. A proper SAR data coverage analysis led to identifying areas or domains with a small number of measurement points and uneven distribution, poorly representative of a proper interpretation of local deformation.

Moreover, it is known that DInSAR techniques allow us to measure only the component of the actual displacement along the satellite LOS, with a generally low sensitivity along the north–south direction [21]. Additionally, the complex topography, typical of mountainous areas, poses an additional challenge. To partially overcome this limit, previous studies assumed that the main movement occurs along the steepest slope and derived the V_SLOPE values based on the C-index [31,43], thus providing a reasonable measure of slope deformation. In recent years, some studies have tried to specifically investigate the relationship between LOS displacement and the real displacement on the slope through a SAR visibility analysis [52,53], highlighting the importance of acquiring and defining the actual slope displacement from DInSAR-derived LOS displacement to guarantee a proper analysis of slope instability interpretation and characterization.

The proposed three-phase methodology tries to partially examine and overcome these issues by combining the analysis of the coverage of the available Sentinel-1 datasets with the evaluation of SAR suitability to measure ground deformation, allowing us to define a ranked investigation of the morpho-structural domain behavior of DsGSDs. Phase one of our methodology ensured the identification of such domains with a suitable PS distribution based on an empirical threshold (C_PS) of measurement points density considered suitable to delineate potential variations and heterogeneities in the DsGSDs movement. The empirical C_PS threshold can vary depending on the compromise between accurate analysis and the possibility of having information about many phenomena; in our cases, it was set as C_PS > 25%. The proposed C_PS ensures that the DInSAR data is spatially distributed and is above a specific density for each domain. The DSR led to the classified ranking of the morpho-structural domains based on their orientation suitability for DInSAR analysis. We defined four qualitative ranking classes, assigning the suitable velocities dataset for each one. The application of the implemented methodology on the two case studies shows that, for low-quality ranking Case 3 (i.e., Croix de Fana DsGSD, north–south-oriented), only the vertical velocity and the V_SLOPE give fairly reliable displacement; instead, the high-quality ranking case (i.e., Valtournenche DsGSD, east–west-oriented) gives the possibility to use more velocity type with ground truth reliability (vertical and E-W velocities, V_SLOPE and LOS velocities). In general, the proposed DSR matrix (Figure 4c; Supplementary Material Figure S2) shows that about half of the possible slope–aspect combinations have high or medium-high ranking, and only a tiny percentage (<5%) is the worst scenario. The C-index matrixes also show that the LOS incidence angle plays a role in defining the slope thresholds of high-quality class; in the case of Sentinel-1, it is relatively high (more than 30°), while for sensors with a low incidence angle (e.g., ERS, 23°), this decreases. In particular, for the Sentinel-1 datasets used in this work, the incidence angle equals 34° for the descending dataset and 39.2° for the ascending one. This characteristic and its global free coverage qualify Sentinel-1 as the best sensor to perform this analysis. Finally, the ground deformation maps are obtained from the interpolated velocity of the selected dataset (from DSR) only on the sectors with adequate C_PS. This constraint aims to improve kinematic analysis only where data is sufficient and reliable with ground truth to reasonably reflect the heterogeneity and complexity of the DsGSD, with a strong segmentation of displacement. The obtained interpolated velocity maps revealed that the higher density provided by the CAESAR-D approach and Sentinel-1 sensor allowed us to locally characterize the kinematic of the whole phenomenon and the variability of deformation between morpho-structural domains, jointly considering measured point distribution and slope orientation with respect to the sensor LOS. The ranked maps may represent a valuable tool for identifying areas of high spatial variability in terms of deformation, specifically concerning their interaction with anthropic structures and infrastructures. Both DsGSDs affect partially urbanized areas; therefore, the slow and widespread deformations represent an element of risk for infrastructure and settlements. Specifically, the settlement of Paquier (Figure 9a), located in the lower part of the DsGSD of Valtournenche, revealed significant horizontal (to the west) movements in correspondence with domain C. On the other hand, the facilities of the hydroelectric plant of Quart cross multiple morpho-structural domains characterized by different deformation rates, subjecting the infrastructure to non-negligible deformation stress (Figure 9b) [54].

The ranked ground deformation maps also pointed out some surficial displacements not directly related to DsGSD movements but, for instance, attributable to minor secondary landslides or debris coverage. In contrast to some authors [43] that impose specific thresholds to entirely remove these surficial deformations, we operated through a geomorphological interpretation of each specific outlier observed. In the Valtournenche DsGSD, the local well-recognizable cluster within domain A, with a peak velocity (V_SLOPE) of −73 mm/yr, was distinctly recognized. The local kinematic analysis revealed that this sector may be separated from the rest of the domain and geomorphologically referred to a nested landslide (i.e., L1), which the actual movement of the DsGSD can distinguish. Conversely, the landslides L2 and L3 appear revegetated and stabilized, without contributing to surface deformation. Instead, the movement registered over the debris coverage in the upper part of sector A should not be excluded from the velocity averaging because it is not possible to discriminate if it is related to shallow creep or deep-seated displacement due to its homogeneous rate compared to the rest of the domain.

5. Conclusions

This paper presents an innovative methodology for the study of ground deformation related to DsGSD phenomena through the effective exploitation of Sentinel-1 data. The three-step approach we propose allows us to locally characterize the morpho-structural domains of these massive phenomena, which typically exhibit heterogeneous spatial deformation, by assessing a suitable PS density and analyzing the sensitivity of LOS displacement measurement to slope orientation. Assigning quality classes to the datasets from the slope, aspect, and C-index parameters provided a suitable approach for the targeted analysis, specifically considering the main orientation of the investigated phenomena. This approach can help partially overcome some limitations related to DInSAR techniques (i.e., the relationship between satellite geometries and topography). In particular, the DSR provides a suitable procedure for ranking the velocity types to interpolate that better suits the considered geometry of the phenomenon. In general, while LOS and resolved east–west and vertical velocities are suitable for high DSR ranking (e.g., DsGSD domains on the western or eastern slope and low gradient) with a low DSR (e.g., a DsGSD that is north–south-oriented), the velocity projected along the slope or the vertical velocity better describe the ground truth. In addition, the ground deformation maps allowed for the detection and discrimination of the surface movements compared to the movement of the DsGSD, leading to a correct interpretation of the results. Notable results are the strong horizontal components (11 mm/yr to the west) of the lower sectors of Valtournenche DsGSD, or the moderate movement (about 12 mm/yr) of the upper part of the Croix de Fana DsGSD revealed by the V_SLOPE. In the case of the Valtournenche DsGSD, the analysis of velocity maps allowed us to ascribe to landslides L1 the high deformation rates (up to 70 mm/yr) observed in domain A. Overall, the methodology has proven to be a powerful and useful tool in the risk assessment and management of deep slope deformation phenomena. In the case of the Croix de Fana site, slow deformations impact several strategic infrastructures, including the penstock. In a risk management scenario, deformation maps can help identify the critical spots where the infrastructures are subject to significant stresses. The proposed methodology has proven to be applicable in any mountainous context and to variably oriented phenomena with good SAR data coverage, providing a functional tool for characterizing ground deformation, and for preparatory element-to-land use planning purposes.