Landslide Mapping and Monitoring Using Persistent Scatterer Interferometry (PSI) Technique in the French Alps

: Continuous geodetic measurements in landslide prone regions are necessary to avoid disasters and better understand the spatiotemporal and kinematic evolution of landslides. The detection and characterization of landslides in high alpine environments remains a challenge associated with di ﬃ cult accessibility, extensive coverage, limitations of available techniques, and the complex nature of landslide process. Recent studies using space-based observations and especially Persistent Scatterer Interferometry (PSI) techniques with the integration of in-situ monitoring instrumentation are providing vital information for an actual landslide monitoring. In the present study, the Stanford Method for Persistent Scatterers InSAR package (StaMPS) is employed to process the series of Sentinel 1-A and 1-B Synthetic Aperture Radar (SAR) images acquired between 2015 and 2019 along ascending and descending orbits for the selected area in the French Alps. We applied the proposed approach, based on extraction of Active Deformation Areas (ADA), to automatically detect and assess the state of activity and the intensity of the suspected slow-moving landslides in the study area. We illustrated the potential of Sentinel-1 data with the aim of detecting regions of relatively low motion rates that be can attributed to activate landslide and updated pre-existing national landslide inventory maps on a regional scale in terms of slow moving landslides. Our results are compared to pre-existing landslide inventories. More than 100 unknown slow-moving landslides, their spatial pattern, deformation rate, state of activity, as well as orientation are successfully identiﬁed over an area of 4000 km 2 located in the French Alps. We also address the current limitations due the nature of PSI and geometric characteristic of InSAR data for measuring slope movements in mountainous environments like Alps. complete cycle of


Introduction
Landslides are among the most widespread geological hazards that pose a major threat to human life, infrastructure and natural environment in most mountainous regions with steep slopes. The slope failure refers to a variety of processes that results in mass wasting down a slope that includes a broad spectrum of ground movements such as rock falls, mudflows, debris flow, artificial fill or a combination of these materials. It may occur in different modes of movements: falling, toppling, sliding, spreading, and flowing, or a combination of these [1]. Gravity is the primary driving factor but there are also other factors controlling the landslide triggering mechanism that favours slope to failure such as heavy rainfall, river erosion, earthquake, snowmelt, volcanic eruption, slope cut etc. The terms "landslide", "mass movement", "slope movement", "slope instability", and "slope failure" are used as synonyms in this study.
A high quality surface displacement field of landslides is important to provide insights into their geometries and mechanical properties [2]. Therefore, monitoring its surface motion over a wide area of space is a key challenge in order to quantify any potential landslide risk [3]. Landslide inventory databases and maps are documenting the identification of landslide processes including the spatial distribution and characteristics such as geometry, volume, total length of area, causal factors, temporal frequencies, and sliding rates of different types of past landslide activities [4,5]. This information might also provide a clue to the locations of future mass movements that are very helpful in predicting the landslide susceptibility and vulnerability in order to support local authorities for hazard mitigation [6][7][8].
In addition, it can also provide input for other applications including land-use planning; developing and refining emergency response plans; as well as prioritizing mitigation measures [9]. However, landslide maps may have different resolutions and details depending on the purpose of the inventory, scales, data collection method, and characteristics of the available imagery [10,11]. In France, for example, national landslide databases provided by French Geological Survey (BRGM) are mostly based on historical records, field inventory and detailed geomorphological formation surveys and aerial orthophotos.
Landslide inventories can be prepared using a variety of conventional and innovative techniques depending on the extent of the study area. Visual interpretation of landslides, geomorphological field mapping as well as stereoscopic imagery techniques are considered as traditional methods based on landform analysis and are only useful over small areas. In addition, these old methods are inevitably subjective [12], cumbersome, time-consuming, prone to error, and most importantly difficult to perform over densely vegetated areas [13,14]. Especially in hilly and mountainous terrain, in-situ methods for analyzing the stability of landslides are not always practical for systematic investigation of slope movements at the regional scale [5]. For these very reasons, the innovative methods and recent technologies help to overcome some of the drawbacks of conventional mapping practices. The modern application of multiple remote sensing technologies including Synthetic Aperture Radar (SAR), optical as well as Light Detection and Ranging (LiDAR) measurements, which represents very valuable complementary data sources relative to conventional mapping and monitoring methods. Especially, the emerging space industry and new technologies will help in reducing the time and resources for systematic updating of landslide inventories over very large areas.
The application of InSAR techniques for detecting and monitoring long-term mass movements has been well documented for different landslide typologies in the last decade owing to its broad spatial coverage, high spatio-temporal resolution and its operational ability during all-weather conditions. For example, the traditional differential InSAR method has been employed for monitoring slow-moving landslides on the order of dm/yr [15][16][17][18][19][20], whereas multi-temporal InSAR methodologies based on the analysis of SAR data stack such as small baseline subsets (SBAS) InSAR [21][22][23] and Persistent Scatterer Interferometry (PSI) [4,[24][25][26][27][28][29] techniques enable detecting extremely slow slope movements on the order of mm/yr. On the other hand, despite its success in landslide kinematics investigations, these multi-temporal InSAR techniques still have certain limitations including severe temporal decorrelation due to the presence of vegetation and soil surface change, unwrapping errors, DEM errors, and atmospheric artifacts. Furthermore, several studies have demonstrated the effect of interferometric pairs selection on the InSAR displacement estimates [30,31].
Different methodologies based on InSAR-derived displacements combined with in-situ observations have been proposed by the scientific community for better assessment of landslide activities on local and regional scales (e.g., [24,25,27]). These methodologies combine thematic data such as topographic indications, geological units, land use, and optical images with ground deformation rates obtained by in-situ measurements and PSI analyses [32,33]. Jointly using InSAR-derived deformation maps and DEM gradient maps is also an effective and systematic way of detecting potential active Remote Sens. 2020, 12, 1305 3 of 22 slow-moving landslide candidates by setting a threshold of topographic slopes based on in-situ geological investigations (e.g., [34][35][36]).
Other proposed methodology aim at fast and systematic identification of actives clusters of the deformation areas derived from InSAR products. These active clusters are corresponding to the spatially aggregated moving points so-called "anomalous areas" [37]. This approach has been proposed by Meisina et al. [37] for slope movements and subsidence detection and its geological interpretation for civil protection purposes. These areas have been called Active Deformation Areas (ADA) since Barra et al. [38]. ADA maps represent critical areas characterized by active ground motions [37,39]. This approach has been followed by Solari et al. 2018 [40] for geohazard mapping by deriving moving areas in Canary Islands (Spain) using a simple and reproducible approach. Tomas et al. [41] presented a methodology to extract and classify ADA [38] and relate them to different potential types of geological procedures in a semi-automatic way based on several thematic layers.
The selection procedure of active deformation areas proposed by [37,[42][43][44][45] is based on a standard deviation (σ vel ) threshold of the deformation map obtained by PS-InSAR technique. Those areas affected by active deformation velocities higher than a critical threshold are considered as potential geohazard risk zones.
Numerous slope instabilities have occurred at different scales in the French portion of the European Alpes over more than a century. Shallow and deep-seated slope failures of different types and rapid/slow mass movements have been documented over large areas in this region. Especially, slow-moving landslide deformations at the level of millimeter or centimeter are giving a good indication for the future slope instabilities. Within the frame of landslide reduction strategies and risk assessment, it is very essential to work for a cost-effective way of zoning susceptible areas for systematic and proper territorial land-use planning. The systematic information on type, pattern, its relation with the morphological and geological characteristics, spatio-temporal evolution, its level of risk, orientation, rate of sliding, and boundaries is lacking in the study area.
For this reason, in the present study we first introduce mean velocity field for the study area ( Figure 1) and the method for automatic detection of slow-moving landslides using these InSAR-PSI products. We then investigate spatial pattern and boundaries of slow-moving mass movements in the high alpine environments in France using InSAR-PSI time series algorithm by combining ADA extraction methodology and classifying them by PSI-derived ground deformation velocity and DEM gradient based on terrain slope and aspect. We also discuss the limitations due to InSAR viewing geometry and local topography for landslide detection and monitoring in mountainous regions. With this goal, we investigated the effect of the principal morphological indices such as local terrain slope and terrain aspect angles on the detection sensitivity of down-slope motions on the direction of line-of-sight by simulating downslope unit deformations evenly distributed along the different slope and aspect angles.

Previous Measurements of Slope Movements in the French Alps
Over the last few decades, numerous studies on detection and monitoring of landslide activities based on different techniques have been conducted by various specialists in the French Alps. In this section, we will highlight only the studies including remote sensing-based approaches carried out in the study area. Most of these studies are focused on either an individual landslide activity on a local scale or multiple slope movements in specific locations. Among these studies, the space-borne techniques using optical and radar images provide valuable insights to investigate landslide kinematics. For instance, the monitoring of deformation pattern of La Clapière landslide has been documented by Schlöger et al., [46] by using L-Band ALOS-PALSAR images with a traditional D-InSAR approach. Another individual La Valette landslide (South French Alps) activity was first documented by differential SAR-interferometry using ERS-1/2 [15] then using the same technique; the landslide superficial displacement rate over nine years has been demonstrated by Squarzoni et al., [47]. The same landslide also studied by Leprince et al., [48] using the subpixel correlation of SPOT-5 images and by Raucoules et al., [49], using the image correlation of X-band SAR images. In all case studies, the observed displacements are in agreement with the ground-truth measurements. Besides this, there are also high-resolution optical images analyzed for the landslide of the Tunnel du Chambon in the French Alps by Desrues et al., [50] who have shown the potential landslide failure. Another large-scale landslide of Harmaliere and new movements and slope

Previous Measurements of Slope Movements in the French Alps
Over the last few decades, numerous studies on detection and monitoring of landslide activities based on different techniques have been conducted by various specialists in the French Alps. In this section, we will highlight only the studies including remote sensing-based approaches carried out in the study area. Most of these studies are focused on either an individual landslide activity on a local scale or multiple slope movements in specific locations. Among these studies, the space-borne techniques using optical and radar images provide valuable insights to investigate landslide kinematics. For instance, the monitoring of deformation pattern of La Clapière landslide has been documented by Schlöger et al., [46] by using L-Band ALOS-PALSAR images with a traditional D-InSAR approach. Another individual La Valette landslide (South French Alps) activity was first documented by differential SAR-interferometry using ERS-1/2 [15] then using the same technique; the landslide superficial displacement rate over nine years has been demonstrated by Squarzoni et al., [47]. The same landslide also studied by Leprince et al., [48] using the subpixel correlation of SPOT-5 images and by Raucoules et al., [49], using the image correlation of X-band SAR images. In all case studies, the observed displacements are in agreement with the ground-truth measurements. Besides this, Remote Sens. 2020, 12, 1305 5 of 22 there are also high-resolution optical images analyzed for the landslide of the Tunnel du Chambon in the French Alps by Desrues et al., [50] who have shown the potential landslide failure. Another large-scale landslide of Harmaliere and new movements and slope deformations around it quantified using aerial photography and image correlation [51,52] reported creeping behavior and transient motion as well as a precursory pattern before failure for the Harmalière case study using freely available Sentinel-2 optical images. The review paper by Delacourt et al. [2] provides the techniques of image analysis to measure the displacement on the slopes with the case studies in French South Alps. UAV-based (unmanned aerial vehicle) high-resolution remote-sensing technique was also carried out for monitoring tempo-spatial dynamics of the mudslide in Super-Sauze in the French Alps [53].

InSAR Datasets
The Copernicus Sentinel-1 C-band SAR mission is a joint initiative of the European Commission (EC) and the European Space Agency (ESA), who is operating the Sentinel-1A and -1B twin spaceborne segments. The Sentinel-1 Terrain Observation with Progressive Scan (TOPS) mode represents an important advantage when compared to other sensors' modes as it provides wide area coverage and short revisit time (6 days over Europe and 12 days globally). Sentinel-1 data are free and openly accessible via various sources, in our case the Copernicus Open Access Hub. In the present study, we processed two independent sets of Sentinel-1A/B imagery, acquired along descending and ascending orbits (Track 88 [Asc] and 66 [Dsc]) to map the surface displacement field of slow-moving landslides in the French Alpine slopes. We discarded data obtained during the snowing season from November to April that reduces the detection capacity due to temporal decorrelation in high altitudes [40]. Hence, this limitation severely reduced the amount of usable interferograms approximately by half. On the ascending track 88, we computed 58 interferograms and on the descending track 66, we formed 50 pairs (see Table 1 and Figure 2 for detail). We employed an individual set of interferograms for each track with a sufficiently high phase coherence pattern over the study area even in the mountainous area and discarded the ones that were affected by temporal decorrelation and/or atmospheric noise correlated with topography, by visual inspection. Even so, due to the nature of the adapted data processing technique, PSI is based on single master interferograms that are using amplitude and phase criterion for persistent scatterer (PS) selection; it does not work optimally when coherence is low (temporal and spatial decorrelation effects), leading to low PS density [54]. For this reason, we have different PS density for the ascending and descending tracks as seen in Table 1 and Figure 4.
The results obtained from the InSAR time series are compared to the official landslide inventory catalogue of the French Geological Services [BRGM] for the selected areas.

Methodology of Landslide Mapping
For investigating the landslide distribution and geohazard activity mapping as well as displacement monitoring using Sentinel-1 data over the selected study area, we followed the procedure based on the estimation of the mean velocity and the time-series of deformation using the PSI approach followed by the extraction of active deformation areas (ADAs) implementing certain basic statistical criteria and a density-based clustering algorithm. Here, we integrated morphological and satellite PSI data with the pre-existing landslide inventory. After estimating the PSI velocities, within a PSI post-processing stage, density and the rate of the local ground deformations are assessed and the PS displacements in LOS direction are projected along the direction of the steepest terrain slope. In the third stage, the result of the velocity projection is used as an input for detecting the active deformation areas by cluster analysis. Figure 3 shows the flowchart of the adopted methodology and outputs used at each processing step. Finally, the results are integrated with the landslide inventory to generate landslide activity map and its distribution with respect to morphological indices.

Methodology of Landslide Mapping
For investigating the landslide distribution and geohazard activity mapping as well as displacement monitoring using Sentinel-1 data over the selected study area, we followed the procedure based on the estimation of the mean velocity and the time-series of deformation using the PSI approach followed by the extraction of active deformation areas (ADAs) implementing certain basic statistical criteria and a density-based clustering algorithm. Here, we integrated morphological and satellite PSI data with the pre-existing landslide inventory. After estimating the PSI velocities, within a PSI post-processing stage, density and the rate of the local ground deformations are assessed and the PS displacements in LOS direction are projected along the direction of the steepest terrain slope. In the third stage, the result of the velocity projection is used as an input for detecting the active deformation areas by cluster analysis. Figure 3 shows the flowchart of the adopted methodology and outputs used at each processing step. Finally, the results are integrated with the landslide inventory to generate landslide activity map and its distribution with respect to morphological indices.

PSI Processing Methodology
All interferograms were generated using the open source (GNU General Public License) InSAR processing system called "Generic Mapping Tools Synthetic Aperture Radar (GMTSAR)" [55]. To correct the topographic contribution to the radar phase, we used the Shuttle Radar Topography Mission (SRTM) 3-arcsecond digital elevation model [56]. All interferograms were computed based on a single master network for PSI analysis. The choice of the master images minimizes the spatial and temporal baselines. The single master stacks of interferograms were processed using the StaMPS software package [57,58]. StaMPS allows the identification of PS points, using both amplitude and phase information. In the first step, the initial selection of PS points is performed based on their noise characteristics, using amplitude analysis. The amplitude dispersion criterion is defined by D Amp = (σ Amp m Amp ⁄ ), where σ Amp and m Amp are the standard deviation and mean of the amplitude in time, respectively [59]. We selected a threshold value of 0.27 for D Amp , which minimizes the random amplitude variability and eliminates highly decorrelated pixels in some areas covered with vegetation, agricultural fields, or snow. Once stable PS targets have been selected based on amplitude analysis, the PS probability is refined by phase analysis in a series of iterations. This process allows the detection of stable pixels even with low amplitude. Once the final selection of PSs has been done, the residual topographic component can be removed. Then, phase unwrapping is performed both spatially and temporally. This analysis enables retrieval of the average Line-Of-Sight (LOS) surface deformation rate maps.
To remove atmospheric effects from interferograms, we used the freely available Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) [60]. This toolbox uses ERA-Interim (ERA-I, Europiean Center for Medium-Range Weather Forecast) numerical weather model datasets [61].

InSAR-Derived Mean Velocity Maps
This study focused on a mountainous area of approximately 4000 km 2 that covers a high density of slow-moving landslides in the French Alps. Figure 4 shows the InSAR-derived mean line-of-sight velocity fields calculated from PSI time-series analysis for the selected areas

PSI Processing Methodology
All interferograms were generated using the open source (GNU General Public License) InSAR processing system called "Generic Mapping Tools Synthetic Aperture Radar (GMTSAR)" [55]. To correct the topographic contribution to the radar phase, we used the Shuttle Radar Topography Mission (SRTM) 3-arcsecond digital elevation model [56]. All interferograms were computed based on a single master network for PSI analysis. The choice of the master images minimizes the spatial and temporal baselines. The single master stacks of interferograms were processed using the StaMPS software package [57,58]. StaMPS allows the identification of PS points, using both amplitude and phase information. In the first step, the initial selection of PS points is performed based on their noise characteristics, using amplitude analysis. The amplitude dispersion criterion is defined by D Amp = (σ Amp /m Amp ), where σ Amp and m Amp are the standard deviation and mean of the amplitude in time, respectively [59]. We selected a threshold value of 0.27 for D Amp , which minimizes the random amplitude variability and eliminates highly decorrelated pixels in some areas covered with vegetation, agricultural fields, or snow. Once stable PS targets have been selected based on amplitude analysis, the PS probability is refined by phase analysis in a series of iterations. This process allows the detection of stable pixels even with low amplitude. Once the final selection of PSs has been done, the residual topographic component can be removed. Then, phase unwrapping is performed both spatially and temporally. This analysis enables retrieval of the average Line-Of-Sight (LOS) surface deformation rate maps.
To remove atmospheric effects from interferograms, we used the freely available Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) [60]. This toolbox uses ERA-Interim (ERA-I, Europiean Center for Medium-Range Weather Forecast) numerical weather model datasets [61].

InSAR-Derived Mean Velocity Maps
This study focused on a mountainous area of approximately 4000 km 2 that covers a high density of slow-moving landslides in the French Alps. Figure 4 shows the InSAR-derived mean line-of-sight Remote Sens. 2020, 12, 1305 8 of 22 velocity fields calculated from PSI time-series analysis for the selected areas characterized by a high density of slow-moving landslides using ascending (track 88) and descending (track 66) orbits. As the deformation is illustrated in LOS direction, positive velocities (cold colors) represent the motion of the ground toward the satellite, and negative velocities (warm colors) represent the motion away from the satellite.
Remote Sens. 2020, 03, x FOR PEER REVIEW 8 of 23 characterized by a high density of slow-moving landslides using ascending (track 88) and descending (track 66) orbits. As the deformation is illustrated in LOS direction, positive velocities (cold colors) represent the motion of the ground toward the satellite, and negative velocities (warm colors) represent the motion away from the satellite.

Limitation of PSI technique in Landslide Detection
The PSI technique for mapping and monitoring of slope mass movements is capable of estimating slow-moving movements with millimeter precision only over the available PSs. According to the International Union of Geological Sciences Working Group on Landslides, landslide classification by velocities consists of seven velocity classes ranging from an extremely slow moving class less than 16 mm/year rate of movement, to extremely rapid slope movements over 5 meter/second [62]. However, in contrast to this landslide classification, the PSI technique can only provide the mean velocity over several years long data in the LOS direction and is not capable of detecting the maximum velocity during the moment of active sliding. Thus, the distribution of the detected active deformations areas in this study mostly focused on slow-moving phenomena such as extremely slow (velocity < 16 mm/yr) and very slow landslides (16mm/yr < velocity < 1.6 m/year).

Limitation of PSI Technique in Landslide Detection
The PSI technique for mapping and monitoring of slope mass movements is capable of estimating slow-moving movements with millimeter precision only over the available PSs. According to the International Union of Geological Sciences Working Group on Landslides, landslide classification by velocities consists of seven velocity classes ranging from an extremely slow moving class less than 16 mm/year rate of movement, to extremely rapid slope movements over 5 m/s [62]. However, in contrast to this landslide classification, the PSI technique can only provide the mean velocity over several years long data in the LOS direction and is not capable of detecting the maximum velocity during the moment of active sliding. Thus, the distribution of the detected active deformations areas in this study mostly focused on slow-moving phenomena such as extremely slow (velocity < 16 mm/yr) and very slow landslides (16mm/yr < velocity < 1.6 m/year). Figure 5 shows the geometry of satellite's heading and azimuth look direction (ALD) for the ascending satellite pass in right looking, as well as the projection of a 3D displacement vector D with components D N , D E and D U in north, east and vertical up directions, respectively, onto the line-of-sight direction [63]. LOS velocities represent the deformation in the one-dimensional projection onto the satellite's Line of Sight (LOS) of a target's motion that actually occurs in all three dimensions. Thus, LOS velocities can be computed as:

Projection from V LOS to V SLOPE
(1) Figure 5 shows the geometry of satellite's heading and azimuth look direction (ALD) for the ascending satellite pass in right looking, as well as the projection of a 3D displacement vector D with components DN, DE and DU in north, east and vertical up directions, respectively, onto the line-of-sight direction [63]. LOS velocities represent the deformation in the one-dimensional projection onto the satellite's Line of Sight (LOS) of a target's motion that actually occurs in all three dimensions. Thus, LOS velocities can be computed as: As seen in Figure 6a, the azimuth look direction (ALD) comprises information about multiple horizontal components (east-west & north-south motions) and is always perpendicular to the satellite heading [64]. Assuming an incidence angle and satellite's heading azimuth of relative to north, the component of the azimuth look direction (D ALD ) can be computed (Figure 5a) as: And the component in LOS direction can be computed as (Figure-5b) : After calculating the direction cosines of LOS, we convert the projection of the LOS deformation rate onto the maximum slope direction (VSLOPE). As PS mean velocity measurements represent the deformation in the one-dimensional projection onto the satellite's LOS of a target's motion that actually occurs in all three dimensions, we partially overcome this limitation given by ground geometry and satellite's viewing geometry by converting LOS velocity vector (VLOS) into slope direction (VSLOPE) [43,[65][66][67]. Thus, this downslope projection facilitates the data interpretation and enriches VLOS information. In translational type landslides, most of the ground-surface deformation occurs mainly along the direction of the steepest slope and, therefore, the most representative component of the mass movements is considered to be parallel to this direction [37,68]. It is, therefore, the projection that is based on the assumption that the sliding velocity is purely parallel to the steepest slope within every single ADA cluster and is performed by using the following equations [45]: As seen in Figure 5a, the azimuth look direction (ALD) comprises information about multiple horizontal components (east-west & north-south motions) and is always perpendicular to the satellite heading [64]. Assuming an incidence angle γ and satellite's heading azimuth of α relative to north, the component of the azimuth look direction (D ALD ) can be computed (Figure 5a) as: And the component in LOS direction can be computed as (Figure 5b): After calculating the direction cosines of LOS, we convert the projection of the LOS deformation rate onto the maximum slope direction (V SLOPE ) ( Figure 6). As PS mean velocity measurements represent the deformation in the one-dimensional projection onto the satellite's LOS of a target's motion that actually occurs in all three dimensions, we partially overcome this limitation given by ground geometry and satellite's viewing geometry by converting LOS velocity vector (V LOS ) into slope direction (V SLOPE ) [43,[65][66][67]. Thus, this downslope projection facilitates the data interpretation and enriches V LOS information. In translational type landslides, most of the ground-surface deformation occurs mainly along the direction of the steepest slope and, therefore, the most representative component of the mass movements is considered to be parallel to this direction [37,68]. It is, therefore, the projection that is based on the assumption that the sliding velocity is purely parallel to the steepest slope within every single ADA cluster and is performed by using the following equations [45]: The coefficient C gives information about the portion (percentage) of "real" displacement (V SLOPE ) [70>69] and represents the sensitivity of the LOS vector against its projection on the slope direction, which is a function of the satellite's viewing geometry and most principal morphological indices: local terrain aspect and local slope angles. The conversion factor C is calculated by: where A is the local terrain aspect with respect to North, and S is the local slope angle [45,69]. These parameters require a Digital Elevation Model (DEM) of the area of interest. For the selected area in this paper, we used 30 m resolution ALOS (AW3D30) DEM. N, E, and H are the directional cosines of the LOS vector, and are calculated from the incidence angle (γ) and line-of-sight azimuth (α) in degrees by using the equations below as calculated in Equations (2)- (4): For certain areas with very low sensitivity for detecting any deformation in the line-of-sight direction, C values approach zero and in this case the velocities in the slope direction tend to infinity. In order to prevent such artificial exaggeration of the V SLOPE , C values fixed to −0.3 if −0.3 ≤ C < 0 and to +0.3 if 0 ≤ C < +0.3 as proposed by Kalia [67].

Detection of Active Deformation Areas (ADA)
The first step in this procedure is to define a threshold value in order to identify active PS points to be used for ADA selection. Standard deviation (σ) of the mean PS velocities is a measure of how far the mean PS velocities fluctuate from the mean and provides information about the noise level of the deformation velocity map. Therefore, it is considered as an indicator of general sensitivity of the deformation map to identify the active deformation areas [37]. Through visual inspection, we experimented with different multiplication factors between 1σ and 2σ of mean velocities. For the selected areas, 2σ is defined as a stability threshold. As a second step, we created two velocity classes by using the above-mentioned threshold; 'Active' as Class 1 and 'Stable' as Class 2. A PS point is considered moving and hereby active if |v| > 2σ vel , where |v| is the absolute value of the mean velocity of each PS point. Thus, Active Deformation Areas belongs to Class 1. If the absolute LOS velocity is smaller than the stability threshold, PS points are classified as 'stable' and belong to Class 2. Following the selection of the moving active points, we used density-based clustering algorithm [70,71] to identify distinctive deformation areas in the active (Class1) data. This algorithm is based on the idea that a cluster (ADA) in the active PS data space (|v| > 2σ vel ) with a high point density is separated from other similar clusters of low point density. By this way, we can cluster the data points into groups and separate them from noise. It requires two parameters, eps (ε) and MinPts. Eps (ε) defines the neighborhood around a data point. If the Cartesian distance between two points is lower or equal to 'eps' then they are considered neighbors and they will belong to the same cluster. In our case, eps parameter defines the radius of influence to consider the active PSs as members of an individual ADA cluster. Through visual inspection, we experimented with different 'eps' values in the selected areas and a value of 30 m has been adopted. MinPtS defines the minimum number of active contiguous data points (PS) within the eps radius to be an ADA. In this work, a value of 5 is adopted. Thus, if the grouped PSs are less than five, they are considered as noise or outliers and cannot represent an ADA cluster.
Finally, for each active deformation area (ADA), we estimated the number of aggregated active PSs; their minimum, maximum and average velocities; as well as the local steepest slope and morphological aspect angle (i.e., local directional mean).
Following the detection of the active deformation areas, we performed local terrain slope and local terrain aspect analysis using the 30 m resolution ALOS Global Digital Surface Model data for the study area. Local terrain slope, which can be regarded as the most basic and important topographic parameter, not only describes the relief and morphology of the earth's surface but is also widely used for many types of environmental applications including landslide monitoring and analysis of mass movements [72]. Local terrain aspect identifies the compass direction of the place with respect to north. It can also be thought of as the slope direction. In this study, we used local terrain slope to classify landslides and distinguish them from other subsidence related signals. Here, we set a threshold of 10 • for the slope gradient based on terrain slopes of known inventoried landslides. The detected ADA with slope angles less than 10 degrees are not considered as landslides in this case. We used the parameter of local terrain aspect to compare landslide velocities with different slope orientations by considering the geometric limitations of InSAR in the mountainous areas and geometric distortions due to viewing geometry. These geometric limitations of monitoring landslides in the areas characterized by steep terrain can be partially overcome by integrating multi-sensor or multi-track SAR data. For this reason, in this study we used ascending and descending tracks over the study area to complement detection capabilities on different slope orientations. This is because the same landslide signal will be detected at different levels of activities (velocities) due to the orientation of the slope with respect to the satellite acquisition geometry. Finally, by integrating these morphological parameters, we provided the statistical distribution of the landslides in the study area.
In order to better illustrate the conversion from V LOS into V SLOPE and the various steps involved for the extraction of the active deformation areas, we focused on a small area characterized by relatively dense slow moving landslides (see Figure 7a).  ) extracted from the projection of the LOS velocity field into slope direction as seen in Figure 6d. (b-c) Identification of distinctive deformation areas using a density-based clustering algorithm using the clustering radius parameter ( ) as 30 meter and with a minimum of five active contiguous data points within this radius. Each colored point represents a core cluster member. Small black points represent outliers. (d) Histogram of landslide number calculated for each ADA cluster with respect to their velocities in the steepest slope direction. Calculated maximum velocities of down-slope deformation for every Figure 7. Extraction of active deformation areas (ADAs) for the selected test area and their classification with respect to primary topographic parameters: slope and aspect. (a) Selected PS velocities that are considered active based on the velocity threshold (|v| > 2σ vel ) extracted from the projection of the LOS velocity field into slope direction as seen in Figure 6d. (b,c) Identification of distinctive deformation areas using a density-based clustering algorithm using the clustering radius parameter (ε) as 30 m and with a minimum of five active contiguous data points within this radius. Each colored point represents a core cluster member. Small black points represent outliers. (d) Histogram of landslide number calculated for each ADA cluster with respect to their velocities in the steepest slope direction. Calculated maximum velocities of down-slope deformation for every single ADA cluster are plotted against mean local terrain aspect (e) and mean local terrain slope (f) within these clusters.

Impact of the Satellite Sensor Viewing Geometry on Monitoring Down-Slope Movement
In order to show the limitations of the satellite sensor looking geometry and the impact of changing morphological indices such as slope and aspect angles on the observed differences on LOS direction, we simulated a translational slide with changing slope (fixed aspect) ( Figure 8) and aspect (fixed slope) (Figure 8) angles specified in terms of single unit vectors from ascending and descending geometries. Spaceborn SAR systems are capable of measuring only the component of the deformation vector projected on the radar line-of-sight; a favorable orientation of the slope is required [73], which makes data interpretation challenging while adding ambiguity to the modeling process [74]. As depicted in Figure 8, the component of the unit deformation vector is projected into the LOS for different slope values higher than 10 • (gray arrows); the true entity of the unit deformation vector can be fully obtained when the slope moves exactly parallel to the sensor-target direction, whereas it is nearly insensitive when the slope (thus the true displacement) is perpendicular to the LOS (azimuth direction). Moreover, acquisitions in both ascending and descending orbits captured a larger portion of the actual deformation of the slope moving either way from the sensor (hot colors), whereas only a small portion of the real deformation for lower slopes have been captured when the slope faces to the sensor (cold colors). That is also the reason why most of the downslope movements detected by both orbits are more pronounced in terms of motion away from the sensors (negative velocities), as only a few targets show motion towards the satellite. In addition, slope angle also affects the spatial resolution, not only the sensitivity in terms of the percentage of the real movements. The satellite-facing site of the slope will have a lower spatial resolution than the slope facing away from the satellite, which is known as the foreshortening and layover effect.
Remote Sens. 2020, 03, x FOR PEER REVIEW 14 of 23 single ADA cluster are plotted against mean local terrain aspect (e) and mean local terrain slope (f) within these clusters.

Impact of the Satellite Sensor Viewing Geometry on Monitoring Down-slope Movement
In order to show the limitations of the satellite sensor looking geometry and the impact of changing morphological indices such as slope and aspect angles on the observed differences on LOS direction, we simulated a translational slide with changing slope (fixed aspect) ( Figure 8) and aspect (fixed slope) (Figure 8) angles specified in terms of single unit vectors from ascending and descending geometries. Spaceborn SAR systems are capable of measuring only the component of the deformation vector projected on the radar line-of-sight; a favorable orientation of the slope is required [73], which makes data interpretation challenging while adding ambiguity to the modeling process [74]. As depicted in Figure 8, the component of the unit deformation vector is projected into the LOS for different slope values higher than 10° (gray arrows); the true entity of the unit deformation vector can be fully obtained when the slope moves exactly parallel to the sensor-target direction, whereas it is nearly insensitive when the slope (thus the true displacement) is perpendicular to the LOS (azimuth direction). Moreover, acquisitions in both ascending and descending orbits captured a larger portion of the actual deformation of the slope moving either way from the sensor (hot colors), whereas only a small portion of the real deformation for lower slopes have been captured when the slope faces to the sensor (cold colors). That is also the reason why most of the downslope movements detected by both orbits are more pronounced in terms of motion away from the sensors (negative velocities), as only a few targets show motion towards the satellite. In addition, slope angle also affects the spatial resolution, not only the sensitivity in terms of the percentage of the real movements. The satellite-facing site of the slope will have a lower spatial resolution than the slope facing away from the satellite, which is known as the foreshortening and layover effect.  Using a similar approach, this time we simulated translational slides with respect to a complete cycle of aspect angle from ascending and descending viewing geometries. The simulated translational slide vectors oriented evenly to all directions from the summit of an idealized cone-shaped synthetic mountain. We examined the sensitivity of the projection of unit deformation vector onto the LOS direction with changing aspect angle by choosing a uniform slope value in all directions (10 • ) that can be considered as a control parameter. We selected a low slope angle in order to show the sensitivity of the motion with respect to satellite: towards the satellite and away from the satellite. As seen in Figure 9, the component of the unit deformation vector is projected onto LOS for different aspect angles ranging from 0 to 360 • ; the true entity of the unit deformation vector can be fully detected when the motion is parallel to the LOS direction, whereas it is nearly insensitive when the displacement is perpendicular to the LOS direction. This limitation arises from the polar orbiting geometry of the satellite. That is why InSAR method is not very sensitive to the motion in north-south direction, and Figure 8c,d depicts this limitation very clearly.
Remote Sens. 2020, 03, x FOR PEER REVIEW 15 of 23 slopes that can be detected on the direction of LOS for ascending (c) and descending (d) viewing geometries. The gray arrows define the deformation of the unit vector. Negative values (hot colors) indicate an increase in the distance from target to satellite, and positive values (cold colors) show the motion towards satellite. The gray zone represents the unfavorable slope angles where the satellite sensor is insensible to measure any deformation due to shadow and layover effect.
Using a similar approach, this time we simulated translational slides with respect to a complete cycle of aspect angle from ascending and descending viewing geometries. The simulated translational slide vectors oriented evenly to all directions from the summit of an idealized cone-shaped synthetic mountain. We examined the sensitivity of the projection of unit deformation vector onto the LOS direction with changing aspect angle by choosing a uniform slope value in all directions (10°) that can be considered as a control parameter. We selected a low slope angle in order to show the sensitivity of the motion with respect to satellite: towards the satellite and away from the satellite. As seen in Figure 8, the component of the unit deformation vector is projected onto LOS for different aspect angles ranging from 0 to 360°; the true entity of the unit deformation vector can be fully detected when the motion is parallel to the LOS direction, whereas it is nearly insensitive when the displacement is perpendicular to the LOS direction. This limitation arises from the polar orbiting geometry of the satellite. That is why InSAR method is not very sensitive to the motion in north-south direction, and Figure 8c and 8d depicts this limitation very clearly.
In such an idealized topographical condition, integration of ascending and descending InSAR measurements over the same region is necessary to resolve the ambiguity in data interpretation.  In such an idealized topographical condition, integration of ascending and descending InSAR measurements over the same region is necessary to resolve the ambiguity in data interpretation.

Landslide Identification Results
The landslide monitoring is performed through the comparision of two acquisition geometries by applying the procedure discussed in Section 3. One-hundred-and-thirty-four ADA have been extracted using ascending track 88, and 111 ADA have been detected using the descending track 66. Figure 10a,b show the impact of local terrain slope on the sensitivity of the line-of-sight measurement along the hill-slope movements; the majority of the landslides in both acquisition geometries are detected in the direction away from the satellite and are mostly concentrated along the slope angles parallel to the line-of-sight direction. These results are in agreement with the LOS-sensitivity analysis obtained from the simulation of the unit deformation vector for different slope angles as depicted in Figure 8c,d. When ADA slope direction velocities estimated from both orbits are plotted against local slope orientation (Figure 10c,d), the results are in a good agreement with the LOS-sensitivity analysis with respect to the aspect angle as seen in Figure 9c,d. The majority of the down-slope movements are in the direction away from the sensor and are concentrated greatly in the azimuth look direction as the simulation results showed. The estimated few landslides in the direction of the flight geometries are arising from the imperfect DEM data and can be considerted as defects. We selected every landslides that shows both directions of motion. Hot colors are those moving away from the satellite and cold colors are those approaching the sensor. The selected landslides in Figure 3a are enlarged and superimposed on perspective Google Earth images shown in Figure 11. From these figures, it can be found that the obvious features of landslides, their boundaries and orientations of slidings can be seen clearly from 3D terrain data of Google Earth. Figure 11 a,b,c and d show unstable slopes moving in two opposite senses with respect to the sensor-target direction. We selected every landslides that shows both directions of motion. Hot colors are those moving away from the satellite and cold colors are those approaching the sensor. The selected landslides in Figure 4a are enlarged and superimposed on perspective Google Earth images shown in Figure 11. From these figures, it can be found that the obvious features of landslides, their boundaries and orientations of slidings can be seen clearly from 3D terrain data of Google Earth. Figure 11a-d show unstable slopes moving in two opposite senses with respect to the sensor-target direction. For example, among them, considereing the topography of the region of interest, the deformation pattern in Figure 11b looks very much like the results obtained from the simulation of the unit deformation vector around an idealized synthetic cone-shape mountain. The unstable slope on the left side of the mountain moves toward a western direction and shows positive deformation values (cold colors); on the other hand, another unstable down-slope motion on the right side of the mountain moves towards an eastern direction and shows negative deformation values.
Remote Sens. 2020, 03, x FOR PEER REVIEW 18 of 23 includes only the cases in high velocities, none of the suspected landslides have been found in the relevant landslide catalogues. In this sense, InSAR PSI technique provides an important contribution for updating the inventory maps mostly in term of slow-moving slope instabilities in the French Alps.

Conclusions
In this work, mapping and monitoring of slow-moving landslides over an area of 4000 km 2 located in the French Alps has been carried out using the PSI approach applied to Sentinel 1 TOPSAR data in ascending and descending orbits. We presented a PSI post-processing procedure for the identification of the suspected active landslides on a regional scale. The second aim of this In order to confirm the detected slope instabilities within the study area, we compared them with the published literature and inventory maps, but since the landslide archive of this region includes only the cases in high velocities, none of the suspected landslides have been found in the relevant landslide catalogues. In this sense, InSAR PSI technique provides an important contribution for updating the inventory maps mostly in term of slow-moving slope instabilities in the French Alps.

Conclusions
In this work, mapping and monitoring of slow-moving landslides over an area of 4000 km 2 located in the French Alps has been carried out using the PSI approach applied to Sentinel 1 TOPSAR data in ascending and descending orbits. We presented a PSI post-processing procedure for the identification of the suspected active landslides on a regional scale. The second aim of this study was to investigate the effects of satellite observation geometry and local topography: slope angle and slope orientation on the applicability of the SAR technique and down-slope motion detection sensibility on the LOS direction.
A combination of the detected landslides from both orbits and their relationship with the local morphology was investigated in order to better picture the capabilities and the abovementioned limitaons of the InSAR for systematic detection and monitoring and failure forecast of alpine slope movements. A comparison of the detected landslides velocitis with respect to the local morphological parametres shows that the results are in accordiance with the sensitivity analysis performed herein, showing the distribution of the slope and aspect angles for the detected landslides activities.
The various PSI post-processing methodologies including the projection of LOS velocity onto the steepest slope and extraction of the active deformation areas (ADA) are explained. A total of 134 ADA from the ascending and 111 ADA in the descending orbit have been detected over an area of 4000 km 2 in a high Alpine environment. The significant majority of the ADA that have been detected from both orbits are predominantly characterized by a slow down-slope motion (<1.6cm/yr) This study confirms that regional scale InSAR is key in improving the active landslides catalogue. In complement with in situ measurements and modeling, this tool could be further used for landslides mechanical characterization. As a next step, in order to develop a full picture of the geohazard activity in terms of landslides ranging from very slow to rapid in velocity and its potential for future movements, additional studies will be needed. With this purpose, future plans include the processing free and open Sentinel-1 SAR data for the whole of the Alpine region in France using both PSI and SBAS approaches for systematically mapping the slope instabilities ranging from very slow to moderate in velocities and "sub-pixel offset tracking" method for detecting fast landslides by using the amplitude channel of SAR images. Besides, in order to better quantify and characterise the landslide and reveal the underlying causes, the integration of InSAR-derived deformation time series with the continuous in-situ monitorings such as seasonal variations in the water table level, water pressure conditions and precipications need to be adressed. Finally, such results will be a key element for the identification and mapping of vulnerable populated areas following the land use data integration.