Semi-automated regional classification of the style of activity of slow rock-slope deformations using PS InSAR and SqueeSAR velocity data

Large slow rock-slope deformations, including deep-seated gravitational slope deformations and large landslides, are widespread in alpine environments. They develop over thousands of years by progressive failure, resulting in slow movements that impact infrastructures and can eventually evolve into catastrophic rockslides. A robust characterization of their style of activity is thus required in a risk management perspective. We combine an original inventory of slow rock-slope deformations with different PS-InSAR and SqueeSAR datasets to develop a novel, semi-automated approach to characterize and classify 208 slow rock-slope deformations in Lombardia (Italian Central Alps) based on their displacement rate, kinematics, heterogeneity and morphometric expression. Through a peak analysis of displacement rate distributions, we characterize the segmentation of mapped landslides and highlight the occurrence of nested sectors with differential activity and displacement rates. Combining 2D decomposition of InSAR velocity vectors and machine learning classification, we develop an automatic approach to characterize the kinematics of each landslide. Then, we sequentially combine principal component and K-medoids cluster analyses to identify groups of slow rock-slope deformations with consistent styles of activity. Our methodology is readily applicable to different landslide datasets and provides an objective and cost-effective support to land planning and the prioritization of local-scale studies aimed at granting safety and infrastructure integrity.


Introduction
Slow rock-slope deformations are common in mountain ranges worldwide. They affect entire hillslopes and displace volumes up to hundreds of millions of cubic meters (Bovis 1990;Chigira 1992;Saroli et al. 2005;Audemard et al. 2010;Agliardi et al. 2013;Crosta et al. 2013;Lin et al. 2013). Although characterized by very low displacement rates (mm/yr to cm/yr; Rott et al. 1999;Wasowski and Bovenga 2014;Frattini et al. 2018;Crippa et al. 2020), they damage infrastructures and may promote secondary failures. Heterogeneity and segmentation of these phenomena result in differential deformations of both surface (e.g. buildings, roads, pipelines) and underground structures (tunnels), strongly affecting the vulnerability of elements at risk with different structural characteristics, topologies and relative positions to different landslide sectors .
The mechanism and internal segmentation of individual slow rock-slope deformations are constrained by inherited tectonic features (Agliardi et al. 2001;Ambrosi and Crosta 2006;Agliardi et al. 2009;Stead and Wolter 2015) both at regional (e.g. major tectonic lineaments) and local scales (e.g. faults, folds, master fractures). In paraglacial environments, slow rock-slope deformations are promoted by stress and hydrological perturbations associated to deglaciation. These trigger progressive slope failure until the development of differentiated rockslides, sensitive to hydrological forcing and mirrored by complex creep behaviour Riva et al. 2018;Agliardi et al. 2020).
Slow rock-slope deformations are recognized by peculiar morpho-structural evidence (Zischinsky 1966;Radbruch-Hall et al. 1976;Massart 1983;Agliardi et al. 2001;Bovis 2013), including both extensional (double-crested ridges, trenches, scarps and counterscarps also associated in half-graben systems) and compressional features (toe bulging, thrusting and folding) that mirror deep deformation patterns and localization. These mass movements include a full spectrum of phenomena, ranging from deep-seated gravitational slope deformations (DSGSD) to large landslides (here LL). DSGSDs affect entire high-relief mountain slopes (>1000 m) involving the ridge crest and exhibit sharp morpho-structural features often controlled by inherited tectonic features. Large landslides are generally smaller than DSGSD (few square kilometres) and can evolve independently or nucleate within DSGSDs. In this case, large landslides show evidence of larger internal deformation and morphological maturity than thehost DSGSD (Agliardi et al. 2012).
Because of the variability of their mechanisms, slow rock-slope deformations exhibit diverse styles of activity, defined by the interplay of displacement rate, kinematics and internal segmentation (Agliardi et al. 2012;Frattini et al. 2018). For example, landslides with translational global kinematics (Cruden and Varnes 1996;Hungr et al., 2014;Rainer et al., 2005;Rotaru et al., 2007) are intrinsically more unstable and prone to higher mobility than landslides with rotational kinematics. On the other hand, slope deformations characterized by similar kinematics and rates can evolve as coherent masses or break up in smaller sectors, depending on their structural segmentation and internal damage.
This complexity, associated with the low typical displacement rates, makes these phenomena difficult to tackle in a risk management perspective. Thus, a regional-scale classification is needed to gather a preliminary characterization of recognized phenomena.
Because of the depth and lack of characterization of basal shear zones, rarely reached by geotechnical and geophysical investigations (Brückl et al. 2006(Brückl et al. , 2013, the assessment of slow rock-slope deformation kinematics is not obvious and strongly relies on surface displacement measurements. In the past, these were limited to sparse traditional ground-based measurements (Bovis, 1990), unable to capture spatial displacements patterns. In the last decades, spaceborne SAR radar interferometry (InSAR) became a powerful tool to characterize ground deformation rates from few millimetres to centimetres per year, acquiring information along the satellite line-of-sight (LOS; Gabriel et al. 1989). The technique maximizes the spatial and temporal coverage of wide areas and proved useful to map and assess the activity of landslides (Colesanti and Wasowski 2006;Del Soldato et al., 2019). In particular, multi-temporal persistent scatterer interferometry (PSI) techniques (e.g. PS-InSAR TM , Ferretti et al. 2001;SqueeSAR TM , Ferretti et al. 2011) allow measuring ground deformations with millimetric precision, making them suitable for regional-scale landslide mapping and inventory studies Colesanti and Wasowski 2006;Rosi et al. 2014;Wasowski and Bovenga 2014;Frattini et al. 2018). However, as an intrinsic limitation of spaceborne SAR platforms, when true landslide deformation vectors deviate from the satellite line-of-sight (LOS), radar sensitivity decreases and slow movements become difficult to resolve (Schlögel et al. 2015;Eriksen et al. 2017;Crippa et al. 2020).
InSAR application to slow rock-slope deformations is further complicated by low signal-to-noise ratio and atmospheric and scatter (e.g. leafs) disturbances, limiting the number of PS, and by the ambiguities related to landslide mechanisms (Notti et al. 2012). Evaluating landslide activity in terms of mean along-slope projected velocity (Meisina et al. 2008;Notti et al. 2014;Tofani et al. 2014) is unsuitable for complex phenomena with kinematics unknown a priori. This can be better tackled using multi-geometry SAR processing combining ascending and descending data (2DInSAR; Eriksen et al. 2017).
In this paper, we use PSI products products (TRE Altamira TM ), acquired over the alpine sector of the Lombardia region (Northern Italy; Fig. 1a), to propose a novel approach to characterize and classify slow rock-slope deformations based on their style of activity. We developed algorithms and Matlab/GIS tools that capture the kinematics, internal segmentation and activity of slow rockslope deformations for a rapid regional-scale analysis. We use multivariate statistical analyses to achieve a robust style of activity classification of mapped slow rock-slope deformations and implement the procedures in a semi-automated workflow that can be readily applied to other landslide datasets in a fast and costeffective way.

Study area
The study area (Fig. 1a) includes the alpine and prealpine sectors of Lombardia region. Here, N-verging Austroalpine and Penninic units are separated from the S-verging Southalpine units by the Insubric Line, a steep E-W trending fault zone active during Oligo-Miocene (Schmid 2017). Penninic units include oceanic lithosphere remnants (e.g. Malenco-Forno units) and portions of European margin (e.g. Adula, Tambò and Suretta nappe) and outcrop in western Lombardia (Valchiavenna and Valmalenco; Fig 1a). Overlying Austroalpine units mainly consist of continental lithosphere of the Apulian margin, including metamorphic basement, intrusives and Mesozoic sedimentary successions (Florineth and Froitzheim 1994;Schmid 2017). They extend from Valmalenco to the eastern part of Valtellina (Fig 1a) and are intruded by two major Tertiary tonalitic and granodioritic bodies outcropping in Valchiavenna-Val Masino (Masino-Bregaglia) and Val Camonica (Adamello). Southalpine units, south of the Insubric Line, form a fold-and-thrust belt made, from N to S, of a Variscan metamorphic basement, Permian volcaniclastic and sedimentary successions and Mesozoic sedimentary succession including carbonate and terrigenous rocks (Fig. 1a).
The alpine sector is characterized by high peaks (up to 4000 m a.s.l.), high local relief and relatively steep topography, all progressively decreasing moving south towards the Southalpine sector and the foothills . This bulk topography is the result of complex interplay between rock type, Cenozoic tectonic forcing and Quaternary exhumation in fluvial and glacial environments (Sternai et al. 2012). During the Last Glacial Maximum (LGM; Ivy-Ochs et al. 2008;Grämiger et al. 2017), a thick ice cap covered the axial alpine sector with ice stages reaching 2800 m a.s.l. (Bini et al. 2009), feeding major valley glaciers (e.g. Adda and Oglio) that carved the topography reaching the Po Plain to the south. Post-LGM deglaciation caused major stress and hydrological perturbations that triggered paraglacial slope response and progressive failure (Grämiger et al. 2017;Riva et al. 2018, Spreafico et al. 2020. Fluvial incision reshaped the valleys especially in the prealpine and foothill areas. These perturbations promoted the onset and spatial distribution of large, slow rockslope deformation Agliardi et al. 2013;Pánek et al. 2019), modulated by lithology and preconditioned by inherited tectonic structure on different scales, from clustering around regional tectonic boundaries  to slope scale controls by master fractures and folds (Agliardi et al., 2001(Agliardi et al., , 2009(Agliardi et al., , 2019aPedrazzini et al., 2016). Slow rock-slope deformations are widespread in axial alpine areas characterized by medium-strength anisotropic (e.g. foliated) rocks, high local relief and valley incision and tend to cluster along major tectonic structures. In areas with sedimentary rocks and relatively low relief (e.g. Southern Alps), giant slow rock-slope deformations are less frequent and often controlled by second-order stratigraphic and structural features (Crosta and Zanchi 2000;Agliardi et al. 2013) that can only be accounted for in site-specific studies.

Materials and methods
A robust characterization of the style of activity of slow rock-slope deformations must consider the interplay between their displacement rates, kinematics and complexity (e.g. segmentation, heterogeneity, internal damage, structural controls). To this aim, we combine an original geomorphological mapping dataset with different PS-InSAR TM and SqueeSAR TM datasets to develop a novel approach to characterize and classify slow rock-slope deformations on a regional scale using deterministic and multivariate statistical analyses.

Geomorphological mapping
We consider 208 slow rock-slope deformations selected from different inventories (GeoIFFI landslide inventory, Ceriani and Fossati 2005;DSGSD inventory, Crosta et al. 2013, including 134 deep-seated gravitational slope deformations (DSGSDs) and 74 large landslides. The latter are extracted from the GeoIFFI database based on areal extent (area >1.5 km 2 , >1 km 2 when interacting with elements at risk) and practical relevance (e.g. actively monitored or undergoing remediation works).

InSAR data
We use 3 different PS-InSAR TM and SqueeSAR TM datasets (TRE Altamira TM ) derived from ERS 1-2, Radarsat and Sentinel 1A/B Cband SAR imagery. Images were acquired between 1992 and 2017 Velocity of the smallest peak (v_Pm) Less frequent peak velocity Velocity corresponding to the smallest velocity frequency peak identified Dispersion of the velocity frequency distribution in each landslide, proxy of velocity heterogeneity Δ Skewness( Δ_SK) Skewness of Δ distribution Asymmetry of the frequency distribution curve of Δ angles (Fig. 4) Δ Median( Δ_M) Median of Δ distribution Median value of Δ angle frequency distribution ( Fig. 4) along ascending and descending tracks ( Table 2) and processed over the entire study area. For each mapped rock-slope deformation (Fig. 1a), we select a reference best dataset, characterized by the highest spatial density of PS (permanent scatterers) and DS (distributed scatterers) and the most suitable acquisition geometry. The latter is selected depending on landslide slope and aspect parameters (Figs. 2a,b) and satellite orbit parameters (LOS and azimuth angle), according to the C index (Notti et al., 2012(Notti et al., , 2014Fig. 2c). This index quantifies the amount of along-slope movement that can be recorded by a SAR sensor with specific orbital and LOS directions. Datasets derived from Sentinel 1A/B images always provide the best spatial coverage and highest number of PS and DS and, when available, are preferred in the analysis. Since velocity values associated to PS and DS data can be locally related to shallow movements of slope deposits (e.g. scree, periglacial features) and not directly related to deep-seated landslides (Meisina et al. 2008;Frattini et al. 2018), we remove PS inside slope deposits. To this aim, we use geomorphological maps accompanying the GeoIFFI inventory, refined by interpretation of Google Earth TM imagery and recent orthophotos.
The robustness of InSAR data analyses is affected by the density and distribution of active PS or DS within individual landslides (Frattini et al. 2018; online Resource Fig.S1), which depend on slope orientation (aspect and inclination) and vegetation cover. Fewer PS and DS limit the statistical robustness of velocity value samples, resulting in a more uncertain definition of landslide segmentation, representative velocities and kinematics. As a limiting condition, less than three PS or pseudo-PS hamper the spatial interpolation of LOS velocity, required to retrieve information on landslide segmentation, degree of activity and kinematics.

Deterministic analysis: segmentation, heterogeneity and kinematics
The simplest InSAR-based descriptor of landslide activity is the mean LOS velocity of active scatterers inside each landslide area. We consider "active" all the PS and DS with coherence > 0.7 and associated mean velocity outside the non-specific uncertainty (noise) range of ±2 mm/yr, commonly used in regional-scale studies (Bianchini et al., 2015;Del Ventisette et al., 2015;Raspini et al., 2018;. However, the mean LOS velocity computed from point-like PS and DS over a landslide area may not represent its actual activity (Frattini et al. 2018). In fact, large slope deformations are often characterized by nested sectors faster than the main body, as well as heterogeneities related to rock mass damage and active morpho-structures (Crippa et al. 2020), resulting in large LOS velocity variability.
We thus develop an objective method to characterize the internal segmentation of each landslide, based on the identification of peaks on the frequency distribution curve of LOS velocity (Fig. 3). In general, landslide bodies with little internal damage or segmentation ( (Notti et al., 2012(Notti et al., , 2014 indicate the amount of along-slope movement that can be recorded by a SAR sensor on each slope according to its morphometry and orbital parameters. High values of C correspond to favorably oriented slopes (green in the map) that are better sensed by the satellites, low values point out landslides with unfavorable orientation (red values) In our analysis, peaks are defined as local maxima exceeding a certain probability density threshold that can be set arbitrarily to detect peaks with different relevance. In this study, we set its value to 1/10 of the maximum probability density. Each detected peak represents the distribution of a specific cluster of velocity values, without providing clues on its spatial distribution. Distinct peaks emerging from the threshold indicate the occurrence of clearly segmented domains, each characterized by a given mean velocity and dispersion. Minor peaks (frequency < threshold) are usually related to outliers or noise. We quantify the activity heterogeneity (i.e. dispersion of velocity values) of each landslide in terms of quartile deviation Q_Dev=(Q3-Q1)/(Q3+Q1), where Q1 and Q3 are the 25th and 75th percentiles of the velocity distribution, respectively. Larger heterogeneity is mirrored by a larger value of Q_Dev.
Slope deformation kinematics can be estimated using a multigeometry analysis of InSAR data, combining measurements provided by ascending and descending tracks to decompose the 2D velocity vector in the E-W trending vertical plane (2D-InSAR; Dalla Via et al. 2012;Eriksen et al. 2017;Crippa et al. 2020). We developed a Matlab TM script to automatically discretize each landslide area into regular square cells (size: 25 m), in which the average LOS velocity of PS and DS from the same acquisition geometry is assigned to cell centroids. For each "pseudo-PS" (i.e. cell centroids for which information from both ascending and descending data are available), the script derives the vertical (V z ) and horizontal (V e ) displacement rate components, the 2D displacement rate vector V T in the E-W vertical plane, as well as its inclination τ (Manzo et al. 2006;Dalla Via et al. 2012;Eriksen et al. 2017;Crippa et al. 2020). The cell size (25 m) was selected to optimize the number of pseudo-PS and the accuracy of computed mean velocity values, avoiding mixing distant scatterers. In slopes with spatially heterogeneous deformations, large cells (100 m) can extend over different deformation domains, resulting in unrealistic estimates of local kinematics. On the contrary, small grid cells (10 m, below the spatial resolution of the sensor, about 20 m for Sentinel-1) result in a reduced probability of finding pseudo-PS.
The local (cell-scale) slope kinematics can be readily inferred (Agliardi et al. 2019b) by observing the difference between the 2D velocity vector inclination (τ) and the local slope dip (α) in each square cell, namely Δ=τ-α (Fig. 4a). Where Δ is positive, we have a dominant downward movement dipping in the slope (Figs. 4b,c). Values of Δ close to zero indicate slope-parallel sliding, while negative values indicate daylighting or bulging movements. The global (i.e. slope-scale) landslide kinematic is inferred by analysing the statistical distribution of the parameter Δ within each entire landslide area (Figs. 4d,e). In particular, the frequency distribution of Δ values of translational landslides tends to be symmetrical and centred around zero or slightly shifted towards positive ("dip-in") values, accounting for the local kinematics of scarp areas (Fig. 4d). Instead, rotational landslides have either bimodal distributions or distributions skewed towards positive values, as the displacement vector becomes steeper within the entire mass and daylights (negative Δ) at the toe (Figs. 4e).
In order to exploit this concept in a regional-scale analysis, we use Δ descriptive statistics (i.e. mean, mode, median, skewness, kurtosis) as inputs to a supervised machine learning analysis through the Matlab TM "classification learner" tool. This trains different classification models (e.g. decision trees, discriminant analysis, support vector machines, logistic regression, nearest neighbours, naïve Bayes and ensemble classification) to select the best-performing ones. We train the model on a subset of 16 landslides (Online Resource Fig. S2) with known kinematics (Agliardi et al. 2001;Allievi et al. 2003;Ambrosi and Crosta 2006;Frattini et al. 2018) and generate predictions for the remaining data. Linear discriminant (Fisher, 1936) proves to be the best predictive model, with a mean accuracy higher than 80% over several (>20) iterations. We also validate the result using a subset of 6 selected case studies (Online Resource Fig. S2) and by local-scale detailed observation of morpho-structures, providing surface expression of the deep mechanisms. Our results show that the median (Δ_M) and in minor portion skewness (Δ_SK) of the frequency distribution of Δ values are effective predictors of global landslide kinematics (Fig. 4f).

Multivariate statistical analysis: PCA and cluster analysis
We run a series of multivariate statistical analyses, including principal component analysis (PCA, Pearson 1901;Hotelling 1933;Fig. 4 Results of the kinematic analysis. (a) quantities considered in the analysis. α: local slope angle; τ: inclination of the 2D velocity vector derived by 2D SAR decomposition; Δ: difference between α and τ. Δ values distinguish (b) along slope movement and (c) dip-in or daylight displacements. Δ value distribution over a landslide area provides an indicator of its main deformation style: (d) translational; (e) rotational. (f) machine learning classification of mapped landslide kinematics according to Δ distribution skewness and median. The background coloured fields correspond to the kinematic groups (R, RT, T) and are defined using linear discriminant decision surfaces Ballabio 2015;Cohen-Waeber et al. 2018) and cluster analysis, to identify the lithological, morphometric, morpho-structural and InSAR-derived activity variables (Table 1 and Table 2) that influence the style of activity of slow rock-slope deformations at the regional scale for classification. After a careful evaluation of the activity of landslides with different dominant lithologies (Online Resource Fig. S1), we excluded this variable from the subsequent analyses. In fact, although rock type is a well-recognized control on the occurrence of slow rock-slope deformations Crosta et al. 2013;Pedrazzini et al. 2016), on regional scale, it appears uncorrelated with activity and kinematics (Online Resource Fig. S3).
In the different analysis steps, we keep principal components with eigenvalues > 1. Cases missing values of some variables, e.g. due to lack of PS, DS or pseudo-PS, are automatically excluded from the PCA according to a listwise deletion procedure.
We use PCA results as inputs for an unsupervised cluster analysis, aimed at detecting similarities and classifying cases into few groups (Massart 1983;Ballabio and Consonni 2013;Ballabio 2015), representative of different styles of activity. For the cluster analysis, we use the K-medoids method instead of the traditional K-means. K-medoids uses actual sample points as cluster centres (Kaufman and Rousseeuw 1990;Jin and Han 2010). Medoids are selected randomly from data objects to form n clusters, and remaining data objects are grouped to minimize the sum of dissimilarities of the surrounding points, making the approach more robust to noises and outliers (Arora and Varshney 2016).
We perform the following multivariate analyses (MV), i.e. PCA and cluster analysis, considering the following subsets of landslides and variables:

Results
Slow rock-slope deformation segmentation, activity and kinematics Our results, validated using field data, show that 57 slow rockslope deformations move as coherent blocks ("homogenous" class, e.g. Mt. Letè in Fig. 5a), while most mapped cases undergo a variable degree of internal segmentation. Segmented landslides are usually characterized by one (55 cases) or two (81 cases) distinct nested sectors, with different LOS velocity than the background, e.g. Mt. Padrio Varadega DSGSD  and Mt. Mater (Crippa et al. 2020). Fifteen cases are very segmented, with more than three nested sectors ( Fig. 5a; e.g. Corna Rossa DSGSD, Agliardi et al. 2018). Modal values of LOS velocity representing the state of activity of main landslide bodies are generally in the range of 3-5 mm/yr (absolute values). When landslide movements are homogeneous, mean and modal velocities tend to converge. Mean landslide velocity (Fig. 5b) usually exceeds the modal one (Fig. 5c) when the main landslide body hosts smaller and faster nested sectors. On opposite, modal velocity exceeds the mean one (Fig. 5d) when faster nested sectors, with abundant PS and DS, involve a high percentage of the total landslide extent (e.g. Piz Groppera large landslide, Valchiavenna n°4 in Fig. 11). Differences between modal and mean velocity are usually in the range of few mm/yr, but can reach several mm/yr (e.g. Belviso n°12 in Fig. 11: mean velocity~13 mm/yr, modal velocity~15 mm/yr ) and up to 1 cm (e.g. Piz Groppera: mean velocity~33 mm/yr, modal velocity~42 mm/yr).

Style of activity classification
MV1: DSGSD and large landslides covered by InSAR data, all variables PCA on the bulk inventory (including both DSGSDs and large landslides) allows characterizing only 166 landslides out of 208. In fact, 42 cases with less than 3 active PS or lacking pseudo-PS are automatically discarded from PCA by listwise deletion.
In this analysis, we select the first 6 principal components (PCs), accounting for about 70% of the entire multivariate space variance, and focus on the relationships between the first three principal components (PC1, PC2 and PC3, Fig. 6), which explain about 49% of the variance. PC1 and PC2 loadings (Fig. 6a) are mainly related to morphometric and morpho-structural variables. The strongest control on PC1 is exerted by the density of immature nested landslides (NB) and deformed nested landslides (DB), which distribution is complementary as suggested by their opposite direction along PC1. PC2 is mainly related to the hypsometric integral (Hi), the relief energy (Δh) and the density of morphostructures (DM). PC3 (Fig.  6b) accounts for InSAR-derived parameters such as modal LOS velocity (v_PM), quartile deviation (Q_dev) and kinematics statistics (Δ_SK, Δ_M).
The PC1/PC2 biplot (Fig. 6a) shows a clear trend along PC1, corresponding to different distributions of cases pre-classified as "DSGSD" and "large landslide" in the inventory. The distinction between the two classes is dominated by their morphological features, such as a different L/W (i.e. more elongated shape for large landslides), and the abundance of deformed nested bodies (DB) in large landslides with respect to DSGSDs, the latter hosting abundant immature nested bodies (NB). The same distinction between DSGSD and large landslide groups is evident in the PC1/ PC3 plot (Fig. 6b). We checked the statistical consistency of the inventory pre-classification into these 2 classes by performing a Kmedoids cluster analysis on the first 3 PCs (Fig. 6c,d), repeating the analysis 20 times to improve the robustness of the classification.

Landslides
The resulting two clusters displayed in PC1/PC2 (Fig. 6c) and PC1/ PC3 (Fig. 6d) plots are in very good agreement with inventorybased classification into large landslides and DSGSDs.

MV2: all DSGSD and large landslides, morphometric and morpho-structural variables
The second bulk-dataset PCA, considering only morphometric and morpho-structural variables, allows characterizing all the 208 landslides with respect to the first 3 PCs, which account for the 67% of the multivariate space variance.
Cluster analysis with PC1, PC2 and PC3 as input variables and a 3-cluster partition of the dataset define three main groups, namely: gm1, gm2 and gm3. The first cluster (gm1, Fig. 7c,d) includes landslides with high density of deformed nested bodies (DB) and well-developed scarp areas (LS), testifying significant accumulated deformation. The cluster gm2 (Fig. 7c,d) includes both large landslides and DSGSD with high values of elongation ratio L/W and affecting a relatively immature topography (high Hi). The third cluster (gm3) mainly includes DSGSDs affecting entire high-relief slopes and characterized by high density of immature nested bodies (NB). These clusters are mainly classified according to PC1, as their boundaries are almost parallel to PC3 axis in the PC1-PC3 plot (Fig. 7c,d).

MV3: DSGSD covered by InSAR data, all variables
This PCA allows establishing links between morphometric, morpho-structural and activity characteristics in DSGSD areas covered by Sentinel InSAR data. These are characterized by high areal PS and pseudo-PS (mean density 43 PS/km 2 ), with a maximum of 259/km 2 in a single DSGSD. We analyse PCA results considering the first three PCs, which account for 50.2% of variance.
Again, PC1 and PC2 are related to morpho-structural and morphometric information (Fig. 8a). When cases are classified according to their kinematics, the PC1/PC2 plot (Fig. 8b) shows a clear trend from the left to the right hand, corresponding to a shift from translational to roto-translational and rotational mechanisms, in agreement with the orientation of Δ skewness (Δ_SK) and Δ median (Δ_M) eigenvectors. A trend in the PC space can be also found by classifying cases by their modal velocity, respectful of segmentation effects and ranging between 2 and 21 mm/yr towards  (Fig. 8c). In general, translational DSGSDs move faster than rotational and roto-translational ones.
These relationships are outlined in a five-class K-medoids classification, based on PCs as input variables in which the resulting groups (gc1 to gc5) have consistent morphometric and morphostructural characteristics and represent different styles of activity, illustrated by the arrows in Fig. 8d, which are oriented as the corresponding variables eigenvectors. Cluster gc1 includes the fastest DSGSDs, characterized by translational kinematics and hosting deformed nested bodies undergoing differential evolution. Clusters gc2 and gc3 include slower DSGSDs, significantly segmented but with different kinematics (gc2 mainly translational, gc3 mainly rotational) and internal deformation (gc2 characterized by distributed damage, gc3 hosting deformed nested sectors). Finally, clusters gc4 and gc5 include slow-moving DSGSDs with mainly roto-translational to rotational kinematics, limited internal damage and immature nested bodies.

MV4: large landslides covered by InSAR data, all variables
The same procedure used in MV3 allows classifying large landslides. These usually affect densely vegetated middle-lower slope sectors; therefore, the areal density of PS and pseudo-PS in these areas is quite low (35PS/km 2 on average), with a maximum of 163 PS/km 2 in a single landslide. PCA interpretation is based on the first three PCs accounting for the 47.2% of variance (Fig. 9a).
The PC1-PC2 plot shows a weak trend with respect to landslide kinematics, with a dominantly translational behaviour on the left and a mainly rotational/roto-translational one to the right (Fig.  9b). Similar distribution is found for other variables such as the densities of immature nested bodies (NB) and deformed nested bodies (DB), mainly related to PC1. Activity-related variables, like velocity quartile deviation (Q_dev) and modal LOS velocity (v_PM), are oriented almost perpendicular to PC1 (Fig. 9a) without any clear velocity trend (Fig. 9c).
Because of the small number of input variables, K-medoids cluster classification only considers two groups, to avoid over splitting the sparse dataset. The groups (Fig. 9d) are mainly discriminated by PC1 values and defined by the density of immature nested bodies (NB) and deformed nested bodies (DB).

MV5: DSGSD and large landslides lacking InSAR data, proximity analysis
Since 42 out of 208 rock-slope deformation lack significant InSAR data coverage, they are not classified in the MV3 and MV4 analyses but only in MV2, based on morphometric and morpho-structural variables and lacking InSAR activity measurements. A two-step proximity analysis is then performed to obtain a complete classification of all the mapped landslides. First, all the mapped landslides were classified in the three main groups extracted from MV2 (gm1, gm2, gm3; Fig. 10) and displayed in a PC1-PC2 plot. Then, all landslides with InSAR data belonging to the different gm groups were further reclassified in terms of the style of activity classes resulting from MV3 and MV4 (Figs. 8 and 9). Eventually, the 42 SAR-blind cases were reclassified through a proximity analysis considering the surrounding points (minimum Euclidean distance) and their class signature (Fig. 10).

Discussion
Redefining the activity of slow rock-slope deformations Despite the apparent simplicity of the concept (Cruden and Varnes 1996), a comprehensive definition of "activity" for complex, slow rock-slope deformations remains elusive. A definition only based on a representative displacement rate can be insufficient to capture the behaviour of slope instabilities evolving over thousands of years (Agliardi et al. , 2019aPánek and Klimeš 2016), for which present-day displacement rates are just a snapshot of longer or variable trends (Riva et al. 2018). Moreover, different phenomena moving at similar rates can have different impacts on elements at risk, depending on their size, mechanisms, maturity and predisposition to collapse Peduto et al. 2017;Nappo et al. 2019).
We propose that a complete definition of the style of activity of large slow rock-slope deformations must include different aspects, namely: displacement rate, segmentation/heterogeneity, kinematics, internal damage and accumulated strain.
Since slow rock-slope deformations are widespread in mountain ranges and especially in alpine areas, we need to quantify such style of activity on the regional scale, to perform a cost-effective screening of the different landslides' evolutionary characters and orient further hazard assessment analyses or local-scale engineering geological studies. To do this, we develop a novel approach that combines information derived from persistent scatterer interferometry and from ad hoc systematic geomorphological mapping, to obtain an objective description of the mapped phenomena.
Our analysis benefits from both deterministic and statistical components. The deterministic component consists of original routines in Matlab TM (i.e. peak analysis for the segmentation analysis and 2DInSAR decomposition and machine learning approach for a kinematic assessment) and GIS (semi-detailed mapping and extraction of morphometric and morphological variables) that allow quantifying the internal segmentation, heterogeneity, kinematics and representative displacement rates of each mapped landslide, through a refined post-processing of persistent-scatterer datasets validated using geomorphological mapping. These routines can be used as stand-alone tools and provide variables that, together with morphometric and morphostructural variables, feed a multivariate statistical analysis aimed at Objective characterization of segmentation and kinematics from PSI data Despite the common application of persistent scatterer data to the assessment of landslide activity Frattini et al. 2018), a full exploitation of their potential for slow rock-slope deformations remains challenging, because of the low signal-to-noise ratio, heterogeneity and complex displacement patterns that can be hardly caught by sparse or unfavourably oriented PS data (Fig. 2). Some authors (Meisina et al. 2008;Notti et al. 2014;Frattini et al. 2018;Aslan et al. 2020) deal with these problems under the assumption of simplified mechanisms (i.e. dominant slope-parallel displacements) not respectful of the complexity of slow rock-slope deformations. The common "unbiased" representation of landslide activity by means of LOS velocity projection along slope (V slope ; Notti et al. 2012Notti et al. , 2014Aslan et al. 2020) facilitates the interpretation of V LOS data, but hampers any interpretation of global landslide kinematics (Meisina et al. 2008;Frattini et al. 2018).
Also, when exploiting multi-geometry 2DInSAR techniques (Manzo et al. 2006;Dalla Via et al. 2012;Eriksen et al. 2017;Crippa et al. 2020) that allow deriving the vertical and horizontal components of the 2D velocity vector lying in the E-W vertical plane, the assessment of landslide kinematics is not straightforward. In fact, (a) objective criteria are needed to integrate point-like scatterer information within a global interpretation; (b) internal segmentation affects landslide displacement patterns. We tackle these problems sequentially in a semiautomated workflow.
First, we objectively assess landslide segmentation, i.e. partitioning of slope deformation into discrete nested sectors with differential evolution (Fig. 3). Since a mean landslide LOS velocity (Fig. 5b), computed on all the active PS, potentially mixes signals from discrete landslide sector as well as scattered noise, we use a peak analysis of the velocity frequency distribution to identify the occurrence of segmentation within the main landslide mass (Fig.  5a). The highest peak (Fig. 3), corresponding to the modal LOS velocity, is a representative descriptor of the displacement rate of the main landslide mass (Fig. 5c). Lower (i.e. less frequent) peaks (Fig. 3) outline the signals of discrete nested sectors emerging from the background displacement signals and from scattered noise.
Our method's ability to outline segmentation is influenced by the abundance and homogeneity of PS and DS in nested sectors (Online Resource Fig. S1, Fig. S4) and by the frequency threshold values used in the peak analysis. Thus, the number of sectors in Fig. 5a must be regarded as a lower-bound estimate. At the same time, our analysis is unable to detect fast-moving landslides within slow rock-slope deformation, due to the temporal baselines of PS-InSAR TM and SqueeSAR TM analyses processed on the regional scale. To achieve an objective evaluation of landslide kinematics, we combine measures from ascending and descending datasets to reconstruct the 2D total displacement vector in the vertical E-W plane. The difference Δ between the inclination (τ) of this displacement vector and the local slope angle allows identifying cells with movements dipping in the slope, parallel to the slope, or daylighting. The statistical distribution of Δ (Fig. 4) over each landslide area provides a signature of its global kinematics. This evaluation of landslide kinematics is biased at slopes where real displacements diverge from the LOS directions, due to slope aspect (i.e. north-or south-facing slopes; Fig. 2b) or to oblique slope movements, and cannot be completely read by the radar sensor (Fig. 2c) due to the intrinsic geometrical limitations of spaceborne SAR platforms. A precise quantification of such bias requires information from additional LOS direction (e.g. other satellite tracks; UAV-SAR ad hoc acquisitions, Hu et al., 2020;GPS or pixel offset tracks), usually unavailable in regional-scale studies, or the a priori knowledge of landslide kinematics, which is our unknown. As we depart from the LOS plane, we can generally expect to underestimate the 2D displacement vector inclination, resulting in a less effective detection of rotational movements. Nonetheless, our proposed approach is suitable for a regional-scale objective screening of landslide kinematics, providing guidance to further local-scale studies.
Our analysis approach can be applied to other large landslide datasets, provided that InSAR data with suitable quality is available and, at the same time, provides information useful for the site-specific characterization of individual landslides.
Style of activity classification of slow rock-slope deformations A practical evaluation of the style of activity of slow rock-slope deformation for land planning, hazard assessment and engineering purposes must account for a suite of kinematic, morphometric and morpho-structural descriptors. Our multivariate statistical analysis allows classifying the studied dataset into 7 groups with fairly homogeneous styles of activity, even including cases for which InSAR data are unavailable to describe the present-day patterns of movement.
First, our results (MV1) outline clear differences between DSGSD and large landslide phenomena according to the considered variables. This is not obvious, because the distinction between these two classes of slope instability is debated and usually only relies on morphometric and morpho-structural criteria, including the affected relief and slope sectors, morpho-structural expression and total accumulated strain ). On the other hand, the very similar magnitude-frequency scaling of hundreds of DSGSD and large landslides suggests that the two classes of phenomena belong to a continuous spectrum of slow rock-slope deformations, where the diagnostic features used for their distinction are constrained by size (Agliardi et al. 2012). Interestingly, multivariate statistical analysis provides a sharp distinction in two groups, well-fitting the geomorphological pre-classification into DSGSDs and large landslides (Fig. 6). More elongated shapes and high density of deformed nested bodies, testifying higher accumulated internal deformation, characterize large landslides with respect to DSGSDs (Fig. 7). This suggests that DSGSDs and large landslides undergo different mechanisms and/or evolutionary stages and should be treated separately.
Multivariate statistical analyses, performed on DSGSD and large landslides separately (MV3 and MV4) and expanded through a similarity analysis to include SAR-blind landslides (MV5), led to the classification of the landslide inventory into seven style-ofactivity classes, five for DSGSDs and two for large landslides. The distributions of the different variables in each group (boxplots in Online Resource Fig. S5) demonstrate the impossibility to use single variables as reliable indicators of the style of activity of slow rock-slope deformations.
Furthermore, the concept of "inactive" landslides implies a capability to correctly measure displacements at all sites, not guaranteed by InSAR techniques due to (a) the inability of radar sensors to illuminate unfavourably oriented slopes and (b) the lack of coherence of interferograms in densely vegetated or debris covered slopes. Since we cannot be sure that SAR-blind landslides are not moving, to avoid underestimation of the related risks, we dropped the concept of "inactive landslides" and classified all the cases in the seven style-of-activity classes.
Our classification provides an accurate mapping (Fig. 11) of slow rock-slope deformations with consistent style of activity. This is useful to identify critical phenomena to prioritize site-specific analyses and to analyse the damage potential of slow rock-slope deformations on specific classes of elements at risk, depending on the most important interaction factors (e.g. rate, volume, heterogeneity).
The class gc1 includes the fastest DSGSDs (5, all covered by InSAR data; Fig. 11), typically affecting high relief slopes (>1200 m) and characterized by dominant translational kinematics ( Fig. 8 and Online Resource Fig. S5). These DSGSDs evolve at modal rates up to 15-20 mm/yr and host deformed nested bodies that undergo faster differential evolution at displacement rates exceeding 25 mm/yr (Online Resource Fig. S5). This class includes the most active DSGSDs in Lombardia (i.e. Corna Rossa, Agliardi et al. 2018Agliardi et al. , 2019aSaline, Agliardi et al. 2001;Mt. Mater, Crippa et al. 2020;11, 10 and 3, respectively, in Fig. 11), that pose major risk to infrastructures (roadways, hydroelectric facilities and lifelines). These are the most important candidates for site-specific studies to quantify their potential for slow-to-fast evolution and catastrophic failure (Crosta et al. 2017;Agliardi et al. 2020).
The other classes include far slower landslides (i.e. modal and mean LOS velocities lower than 5-10 mm/yr) that cannot be distinguished in terms of displacement rates (Online Resource Fig. S5), but are characterized by different mechanisms and degrees of segmentation, internal damage and accumulated strain (Online Resource Figs. S5 and 11). Clusters gc2 (32 DSGSD, including 7 SAR-blind cases) and gc3 (21 DSGSD; Fig. 11) include phenomena affecting middle-high relief slopes and characterized by significant segmentation, but with different kinematics and internal deformation style. gc2 are dominantly translational and characterized by significant distributed damage (e.g. Bosco del Conte DSGSD; (Agliardi et al. 2009); Fuipiano DSGSD; (Forcella, 1984); 8 and 14, respectively, in Fig. 11), while gc3 are dominantly rotational and segmented into deformed nested sectors (e.g. Masuccio DSGSD, (De Finis et al., 2015); 6 in Fig. 11). Although these DSGSDs are slow, their strong heterogeneity must be considered in practical engineering geological problems including (a) the definition of potential collapse scenarios, i.e. actually evolving nested volumes constrained by segmentation; (b) the quantification of rock mass properties for stability modelling, strongly depending on internal damage; and (c) the prediction of damage to infrastructures that is significantly complicated by differential displacements in heterogeneous landslide masses. Finally, classes gc4 (38 DSGSD, including 10 SAR-blind cases) and gc5 (40 DSGSD, including 2 SAR-blind cases) include slow-moving DSGSDs characterized by either dominant roto-translational (gc4; e.g. Lake Palù DSGSD ; 5 in Fig. 11) or rotational kinematics (gc5; e.g. Albenza DSGSD; (Forcella, 1987) 13 in Fig. 11). Both classes include phenomena with limited internal damage and segmentation into poorly deformed nested bodies (Online Resource Fig. S5).
Style of activity classes defined for large landslides, i.e. LL1 (32 large landslides, including 12 SAR-blind cases) and LL2 (40 large landslides, including 11 SAR-blind cases), differs according to morpho-structural characteristics that witness different geomorphological maturities and accumulated deformations (Online Resource Fig. S5). In particular, class LL1 (e.g. Mt. Solena; 9 in Fig. 11) includes translational landslides affecting slopes with relatively low relief (Online Resource Fig.   S5), characterized by limited scarp areas and segmented into immature nested masses, suggesting a limited accumulated strain. Instead, class LL2 (e.g. Pizzo Groppera; 4 in Fig. 11) mainly includes roto-translational and rotational landslides, affecting high-relief slopes and characterized by large accumulated deformations.

Conclusions
Slow rock-slope deformations are widespread in alpine valleys. Despite their limited displacement rates (mm to few cm per year), these phenomena are long-lived and deform slowly, possibly until massive sector collapse, posing major risks to infrastructures and human lives. Our regional-scale analysis allows performing a complete screening of ongoing slow rock-slope deformations, with the following results: -New semi-automated, objective methods and tools to characterize the internal segmentation and kinematics of slow rockslope deformations. These methods, readily applicable to other datasets, highlight phenomena subjected to a certain degree of segmentation, with the presence of internal nested bodies that  (Forcella 1987); 14: Fuipiano DSGSD (Forcella and Rossi 1987) can potentially evolve into faster landslides, and their global kinematics. -A statistically based evidence that DSGSD and large landslides have different mechanisms and/or evolutionary stages, with different associated morphological features that testify higher accumulated internal deformation for large landslides with respect to DSGSD. -A statistically based classification of slow rock-slope deformations based on their style of activity, including the contributions of displacement rate, segmentation, kinematics and internal damage. This classification supports regional-scale land planning and the prioritization of detailed site-specific mapping, monitoring and modelling studies, required to understand and deal with issues posed by slow rock-slope deformations. to geomorphological mapping. We thank Dr. Andrea Manconi for fruitful discussions on the potential and limitations of regionalscale InSAR applications and Marco Previtali for its tips and tricks with Matlab programming. We are also indebted to Massimo Ceriani (Regione Lombardia) for providing some PS datasets.
Author contribution C. Crippa and F. Agliardi conceived the research and wrote the paper. E. Valbuzzi, F. Agliardi and M. Spreafico edited the geomorphological mapping dataset. C. Crippa performed InSAR data analysis and multivariate statistical analysis. P. Frattini contributed to statistical analysis. P. Frattini and G.B. Crosta revised the paper and contributed to the discussion.

Funding
Open access funding provided by Università degli Studi di Milano -Bicocca within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/ 4.0/.