Renaud Toussaint

Aurore Carrier

Alexandre Mathy

Lionel Lorier

Keywords: Time-lapse photography, optical image processing, stereophotography, deformation analysis, landslide

Cette thèse a été réalisée au sein de l'Université de Strasbourg / IPGS et de la société SAGE grâce au partenariat CIFRE/ANRT. Elle n'aurait pu aboutir sans le soutien et la collaboration de plusieurs personnes. Je tiens tout d'abord à remercier Lionel Lorier et François Blanchet, co-gérants de la société SAGE, pour m'avoir donné l'opportunité de réaliser ma thèse en tant qu'ingénieur-doctorant et pour m'avoir donné les moyens d'y parvenir. Je les remercie de la confiance qu'ils m'ont accordée tout au long de la thèse et de l'accueil au sein de l'entreprise lors de mes passages parmi eux. Je les remercie également pour avoir accepté ma charge d'intervenante extérieure à la Faculté de Géographie de Strasbourg. Ce fut une expérience prenante mais enrichissante ! Chercher sans boussole ne m'aurait menée pas bien loin. . . .Je remercie Jean-Philippe Malet (Directeur de thèse) et Ombeline Brenguier (Responsable scientifique au sein de l'entreprise) pour leurs idées, leurs conseils et leur confiance. Je remercie tout particulièrement Jean-Philippe Malet pour sa patience et les nombreuses heures de relecture et Ombeline Brenguier pour m'avoir initiée à la photogrammétrie. Je

Summary

The term gravitational motion is a generic term that refers to several types of ground destabilization that have different kinematic and geomechanical characteristics. The most common ones are mudflows, blockfalls, landslides and volcanic wall slides. The ground displacement velocities vary according to the type of instability and the state of stability in which the environment is. For example, rock falls are fast moving events and can be sudden, whereas landslides can move over a short time when they are quick (a few m/s) or a long one when very slow (a few mm/year).

One of the movement characteristics that has proven to be most appropriate in a surveillance context is the velocity of the surface movement. There are several geodetic techniques including those that allow the acquisition of positions of points on the ground such as GPS or automatic theodolite surveying.

They provide local but accurate (< mm) measurements. Other techniques have been developed over time to ensure wider spatial coverage. They are called remote sensing techniques and are divided into two categories: active sensors such as LiDAR, radar interferometry, etc., and passive sensors such as satellite, aerial or terrestrial optical imagery, thermal imagery, etc.. In an operational context that requires a certain flexibility in data acquisition and processing, we were interested, in this thesis, in the contribution of terrestrial optical images from a fixed camera for the understanding and analysis of landslides and its complementary contribution to classical in-situ measurements. The use of cameras such as single lens reflex cameras (SLR) is increasingly used in the monitoring of natural surface movements. A significant advantage of this instrument is the low cost of the camera compared to terrestrial LiDARs for example. The resolution of the sensor yet remains high (> 18 MPix). The information extracted from the optical images is both qualitative and quantitative, as it enables changes in surface morphology to be identified and movements to be quantified.

During the temporal evolution of these gravitational phenomena, several issues stand out such as: what are the spatial limits of movement? What are the velocities of deformation? What is the prediction of the rupture date? To answer these questions, we have developed tools and processing methodologies for the analysis of long time series in order to extract a maximum of information.

The work was carried out through the analysis of four time-series acquired on four high-stake sites located in the French Alps. Each time series contributed to the development of new techniques and new approaches to acquisition and processing. Three cameras were installed with the aim of monitoring the long-term dynamics of three unstable slopes in a pre-event phase (a phase of significant acceleration). The fourth camera was installed in a post-event monitoring context. The acquisition of long time-series led to the development of generic tools which are as automated as possible and which can be adapted to different situations. These situations can be geometric (camera distance to a variable object) or geomorphological (size of the movement studied, velocity of surface deformation).

The tools developed have been gathered in a toolbox called TSM (Tracing Surface Motion, (Desrues et al., 2019)). TSM is a modular processing toolbox combining tools for: 1. Camera calibration; 2. Image selection (sorting to eliminate images obstructed by poor weather conditions); 3. Image co-registration in order to correct camera movements induced by wind, temperature change or by the movement of its support; 4. Change detection based on image correlation techniques; 5. Filtering of correlation results in order to detect only significant movements; 6. Metric transformation of the pixel displacements resulting from the correlation; 7. Eventually, a stereoscopic approach can be used instade of (6) for the reconstruction of spatially dense 3D deformation fields; 8. Analysis of displacements from a kinematic and mechanical point of view for a better understanding of the phenomenon.

In order to assess the influence of the parameters of these different steps on the TSM results, a sensitivity analysis was carried out by testing many of them such as the image modality, the correlation parameters, the co-registration accuracy, the correlation strategy applied to the time series as well as the accuracy of the digital terrain model integrated in the calculation and the parameters introduced in the camera calibration. TSM has been applied to two case studies and has demonstrated its robustness and generality in the case of a landslide characterised by velocities of less than 10 cm.day -1 and in the case of a faster landslide with velocities greater than 50 cm.day -1 .

From the displacement fields and velocities, several secondary products can be generated depending on the objectives of the case studies and the phase of the landslide evolution:

• In the case of long-term monitoring, we looked at the general behaviour of the phenomenon: are there relationships between the movement and some triggering factors of its acceleration? Are the observed accelerations spatially homogeneous on the scale of the slope? What type of deformation is involved? What is the estimated thickness of the moving mass?

We have therefore implemented modules in TSM that allow us to link events over the long term with meteorological data such as rainfall, and seismological data. We have also implemented deformation analysis tools, when the approach used is stereoscopic, in order to be able to highlight compression, expansion and shear zones. The thickness of the moving mass results from this last analysis. It is based on the law of mass conservation and integrates two major contributions: a rheology and a geometric contribution. v • In the case of pre-event monitoring, we were interested in the prediction of the time of failure based in particular on the relation linking the inverse of the displacement velocity to time (Fukuzono, 1985) and in the application of physical models to characterise the mechanical behaviour involved in the process (Helmstetter et al., 2004).

• In the case of post-event monitoring, we were once again interested in the characterisation of the slope since other points of interests may arise once the event has passed. Quantifying residual movements is indeed important in order to detect or not areas destabilised by the departure of the masses. Defining the new state of stability of the slope is also crucial before any human intervention on site. vi

Résumé

Le terme mouvement gravitaire est un terme générique qui fait référence à plusieurs types de déstabilisations du sol qui présentent des caractéristiques cinématiques et géomécaniques différentes. Les plus communes sont les coulées de boue, les chutes de blocs, les glissements de terrain ou encore les glissements de parois volcaniques. Les vitesses de déplacements varient suivant le type d'instabilité et suivant l'état de stabilité dans lequel le milieu est. Par exemple, les chutes de blocs sont des événements qui évoluent rapidement et qui peuvent être soudains tandis que les glissements de terrain peuvent évoluer sur des temps courts lorsqu'ils sont rapides (quelques m.s -1 ) comme sur des temps longs lorsqu'ils sont lents (quelques mm.an -1 ).

Une des caractéristiques des mouvements qui s'est révélée être la plus appropriée dans un contexte de surveillance, est la vitesse de déplacement de surface. Il existe plusieurs techniques géodésiques dont celles qui permettent l'acquisition de position de points au sol telles que le relevé par GPS ou par théodolite automatique. Elles fournissent des mesures locales mais précises (< mm). D'autres techniques se sont développées au cours du temps afin d'assurer une couverture spatiale plus large. Elles sont appelés techniques de télédétection et se divisent en deux catégories : les capteurs actifs tels que le LiDAR et l'interférométrie radar, et les capteurs passifs tels que l'imagerie optique satellitaire, aérienne ou terrestre et l'imagerie thermique.

Dans un contexte opérationnel qui demande une certaine flexibilité dans l'acquisition et dans le traitement des données, nous nous sommes intéressés, dans cette thèse, à la contribution des images optiques terrestres issues d'un appareil photographique fixe pour la compréhension et l'analyse des glissements de terrain et de son apport complémentaire à des mesures classiques in-situ. L'utilisation des appareils photographiques tels que les caméras reflex à focale fixe (Single Lens Reflex -SLR) est de plus en plus utilisée dans le monitoring de mouvements de surface naturels. Un des avantages non négligeable de cet outil est le prix bas-coût de l'appareil comparé à des LiDARs terrestres par exemple. La résolution du capteur reste cependant élevé (> 18 MPix). Les informations qui sont extraites des images optiques sont à la fois qualitatives et quantitatives puisqu'elles permettent en effet d'identifier des changements de morphologie de la surface et d'en quantifier les déplacements.

Au cours de l'évolution temporelle de ces phénomènes, plusieurs problématiques se détachent telles que: quelles sont les limites spatiales du mouvement ? Quelles sont les vitesses de déformation ? Quelle est la prédiction de la date de rupture ?

Pour répondre à ces questions, nous avons donc développé des outils et des méthodologies de traitement pour l'analyse de séries temporelles longues afin d'en extraire un maximum d'information.

Le travail s'est réalisé au travers de l'analyse de quatre séries temporelles acquises sur quatre sites à enjeu situés dans les Alpes Françaises. Chaque série temporelle a contribué au développement de vii nouvelles techniques et de nouvelles approches d'acquisition et de traitement. Trois caméras ont été installés dans le but de surveiller les dynamiques à long terme de trois pentes instables dans une phase de pré-événement (phase d'accélération significative). La quatrième caméra a quant à elle, été installée dans un contexte de suivi post événément. L'acquisition de séries temporelles longues a mené au développement d'outils génériques qui soient les plus automatisés possibles et qui soient adaptables à différentes situations. Ces situations peuvent être de l'ordre géométrique (distance entre la caméra et l'objet variable), ou bien géomorphologique (taille du mouvement étudié, vitesse de déformation de surface).

Les outils développés ont été regroupés dans une toolbox appelée TSM pour Tracing Surface Motion (Desrues et al., 2019). TSM est une méthode de traitement modulaire associant des outils pour:

1. La calibration de la caméra; 2. La sélection des images (tri permettant d'éliminer les images obstruées par des conditions météorologiques dégradées); 3. La co-registration des images afin de corriger les déplacements des caméras induits par le vent, le changement de température ou bien par le mouvement de son support; 4. La détection de changement basée sur les techniques de corrélation d'images; 5. Le filtrage des résultats de corrélation afin de ne détecter que les mouvements significatifs; 6. La transformation métrique des déplacements pixels issus de la corrélation. 7. Le module (6) peut être remplacé par le module (7) lorsque nous utilisons une approche stéréoscopique pour la reconstruction de champs de déformation 3D dense spatialement; 8. L"analyse des déplacements d'un point de vue cinématique et mécanique pour une meilleure compréhension du phénomène. Afin d'estimer l'influence des paramètres de ces différentes étapes, sur les résultats de TSM, une analyse de sensibilité a été menée en testant plusieurs de ces paramètres : la modalité de l'image, les paramètres de corrélation, la précision de la co-registration, la stratégie de corrélation appliquée à la série temporelle ainsi que la précision du modèle numérique de terrain intégré dans le calcul et les paramètres introduits dans la calibration de la caméra. TSM a été appliqué sur deux cas d'études et a démontré sa robustesse et sa généricité dans le cas d'un glissement caractérisé par des vitesses de déplacement inférieures à 10 cm.jour -1 et dans le cas d'un glissement plus rapide dont les vitesses sont supérieures à 50 cm.jour -1 .

A partir des champs de déplacements et de vitesses de déplacements, plusieurs produits secondaires peuvent être générés selon les objectifs des cas d'étude et selon la phase d'évolution dans laquelle se trouve le glissement de terrain : viii • Dans le cas d'une surveillance à long terme, nous nous sommes intéressés au comportement général du phénomène : existe-t-il des relations entre le mouvement et des facteurs déclencheurs ? Les accélérations observées sont-elles homogènes spatialement à l'échelle du versant ? Quel type de déformation est impliqué ? Quelle est l'épaisseur estimée de la masse en mouvement ? Nous avons dès lors implémenté dans TSM des modules permettant de mettre en relation les événements sur le long terme avec des données météorologiques telles que la pluie ou encore des données sismologiques. Nous avons également implémenté des outils d'analyse de déformations lorsque l'approche utilisée est stéréoscopique, afin de pouvoir mettre en évidence des zones de compression, de dilatation et de cisaillement. L'épaisseur de la masse en mouvement découle de cette dernière analyse. Elle se base sur la loi de la conservation de la masse et intègre deux contributions majeures : une contribution rhéologique et une géométrique.

• Dans le cas d'une surveillance proche de la rupture du versant, nous nous sommes intéressés à la prédiction de la date de rupture en nous basant notamment sur la relation reliant l'inverse de la vitesse de déplacement au temps (Fukuzono, 1985) et à l'application de modèles physiques pour caractériser le comportement mécanique impliqué dans le processus (Helmstetter et al., 2004).

• Dans le cas d'une surveillance post événement, nous nous sommes intéressés une nouvelle fois à la caractérisation du versant puisque d'autres enjeux peuvent avoir lieu une fois l'événement passé. Quantifier les mouvements résiduels est en effet important afin de détecter ou non des zones déstabilisées par le départ des masses. Définir le nouvel état de stabilité du versant est également crucial avant toute intervention humaine sur site. Processing pipeline of the Tracing Surface Motion (TSM) toolbox (Desrues et al., 2019).

1.1 Schematic surface velocities of a gravitational movement after Mirgon et al. (1993) . . . 1.2 Description of the most common types of instabilities based on Varnes (1978) and Cruden and Varnes (1996) classification. The illustrations are from the United States Geological Survey and are available online (USGS, 2004). DSGSD stands for Deep-Seated Gravitational Slope Deformation and defines sliding movements that affect large slopes over a long period of time. The camera symbol indicates the velocity interval in which time-lapse photography monitoring is most appropriate. . . . . . . . . . . . . . . 1.3 Schematic framework that connect the risk management cycle (orange circle), which is managed by the local authorities, the time-related progress of disasters (white circle), the temporal stages of the monitoring (blue circle) and the observation frequency (gray circle).

Modified after Stumpf (2013) and Alexander (2002) Basento river (Italian Southern Apennines) from (Calcaterra et al., 2012). (b) Photograph of a total station set up in front of the Viella landslide (French Pyrenees) (source: J.-P. Malet). (c) Photograph of an extensometer (fixed site visible) from (Malet et al., 2002).

(d) Principle of the direct levelling form (Schultz, 1995). . . . . . . . . . . . . . . . . . 1.6 Radar principle. a) Radar imaging geometry from Bamler (2000). b) Example of SAR image from Yan et al. (2012). Spatial representation of subsidence over Mexico City (2002 -2007) , 17, 35, 4, 12, 22, 16, 26, 23 (see Figure 3.9 for the localization) according to their directions. Green point shows the snow threshold calculated from the polynomial fit. This is the equivalent of finding the minimum local of the polynomial fit. Some pixels may be considered as snow pixel while they are not. series and the corresponding purpose of the surveillance (defined according to (Mazzanti and Pezzetti, 2013)). The usability is defined as follows: "+++" means that optical terrestrial images may provide relevant results and are useful for landslide understanding.

"++" means that optical terrestrial images are useful but some requirements have to be considered. "+" means that the data and technique could be applicable but only on specific

Introduction

General introduction

Assessing the risk of gravitational processes that may impact slopes, volcano flanks or ice glaciers, is a global challenge since some of these processes can jeopardize either infrastructure or lives. Due to the large numbers of gravitational movements around the world, these events have been recorded to built dense databases to help future assessment (Costa and Schuster, 1991;Pereira et al., 2014). The most frequent gravitational movements encountered are rockfalls (Royán et al., 2014;Kromer et al., 2015;Sala et al., 2019), landslides (Travelletti et al., 2008;Quinn et al., 2011;Valentin, 2018) and volcanic flank slides (Walter, 2011). They may be characterized by several parameters, such as their size and volume, the type of material and the deformation rate (Varnes, 1978;Cruden and Varnes, 1996). The motion can be slow to fast, and continuous or by phases of acceleration/deceleration. The understanding of these events involves several disciplines from geology to geophysics and may call for different monitoring techniques.

Monitoring can be defined as an analysis carried out: (1) by science observatories over long time periods, for the scientific understanding of processes;

(2) to document the pre-event/event dynamics to anticipate hazards and risks. The latter is conducted in operational contexts;

(3) to document the post-event dynamics when residual movements and the long-term stability need to be quantified. The latter case is often associated to mitigation works such as drainage. This means that monitoring requires series of observations acquired over time on a periodic or continuous basis. The analysis of these observations should provide information on temporal changes of the movements.

Motion analysis can be performed at different stages of the temporal evolution of gravitational processes. By representing the cumulative surface displacements in time (Figure 1), we can notice that several states compose and define the motion: a steadily deforming is for example expressed by a linear trend implying a constant velocity. It can be followed either by an occasional reactivation or a failure stage. An occasional reactivation defines a destabilization of the slope essentially caused by external factors such as rainfall, snow melt or seismicity, and reaches a new state called, according to Leroueil et al. (2012), active landslide. Several features are used and combined together to provide information on the temporal evolution of the processes. They include for instance, geomorphological and geological observations (e.g, evolution of crack opening), mudflow detection and/or water flow measurements, if there is for instance a river at the foot of the slide. The latter measurements are usually acquired with infrared cameras giving a picture of what is happening day and night, before an operator actually moves to the site. However, the most relevant parameter for characterizing the slope activity, is the motion (incremental/cumulative displacement, displacement rate) of the ground surface since it provides information on the kinematics and the mechanical behavior (Cruden and Varnes, 1996).

Displacement information can be acquired by various geodetic techniques, by comparing the geodetic positions of given natural or artificial targets between two dates (Wang, 2011;Benoit et al., 2015).

These particular techniques, which are in-situ and point-based, allow the analysis of the geometrical properties of the moving mass and of the mechanical behavior of the slopes by analysing time series of displacements. Similar information is also provided by remote sensing methods which have the advantage of proposing an image of the motion and thus more dense spatial coverage and possibly short revisit times (Mazzanti et al., 2012;Sun et al., 2015;Stumpf et al., 2017). Among the remote sensing methods, there is optical imagery, which has the possibility to provide both qualitative information (detection of geomorphological changes, detection of snow) and quantitative information (such as ground motion calculated by image correlation techniques (Lacroix et al., 2015;Stumpf et al., 2017).

In the case of monitoring, more and more studies are using low-cost remote sensing methods, that are easy to use and to deploy (Pham et al., 2014;Gance et al., 2014;Schwalbe and Maas, 2017;Kromer et al., 2019). Thanks to the large consumer market, digital passive sensors, such as Single Lens Reflex (SLR) cameras, are increasingly available, while the sensor resolution is constantly rising. By combining time series of terrestrial optical imagery and the classical in-situ techniques (such as extensometers, crack-meters), the amount of collected information is diversified, better distributed and denser.

The acquisition of terrestrial optical images can be performed in two principal modes: the first mode is low frequency data acquisition, when acquisition campaigns are spaced out over several weeks or months.

Figure 2: Some approaches that can be adopted for a terrestrial optical image acquisition. (a) Monoscopic approach for the analysis of the displacement of the volcano Mount St. Helens (Walter, 2011). (b) Stereoscopic approach for the analysis of the Argentière glacier displacement (Fallourd et al., 2010). (c) Multi-view approach for the analysis of rock slope deformation within Clear Creek Canyon (Kromer et al., 2019).

The second mode is high-frequency data acquisition with the time between two acquisitions being close to the hour or the day according to the velocity of the slope movement. The latter are therefore more continuous measurements. For each mode, two approaches can be adopted: a monoscopic approach and a stereoscopic or multi-view approach (Figure 2). In the first case, only one camera is used whose position remains fixed for the entire duration of the acquisition. The processing of monoscopic images provides, among others, time series of 2D displacement fields. If more than one camera is used (or more than one line of sight), we are in the stereoscopic or multi-view approach providing access to 3D

surface motion and the possibility to create Digital Elevation Models (DEMs) (Hasegawa et al., 2000;Pierrot-Deseilligny and Paparoditis, 2006;Stumpf et al., 2016). When data acquisition is carried out at high frequency (i.e., continuous measurement) from one or more fixed cameras, we are in the case of time-lapse photography. This acquisition mode provides spatially and temporally dense observations. The time intervals in which time-lapse photography is the most appropriate (see Figure 1), depend on the sensor and acquisition parameters which are linked to the expected ground displacement rates. From the image analysis it is possible to extract surface change information such as the surface velocity evolution over time (Figure 3 (bottom)). A typical trend in surface velocity over time can be broken down into several stages (Mazzanti and Pezzetti, 2013): the long-term monitoring that is composed of the slope characterization and the pre-failure stage, and the early-warning stage. For each stages, there are several technical (metrology) and science questions that can be raised:

• From a geomorphological perspective: what are the spatial boundaries (lateral and top limits) of the motion?

• From a kinematic perspective: what is the 2D/3D velocity field of the movement?

• From a mechanical perspective: what are the deformation regime and rheology of the moving masses?

• From a forecasting/modelling perspective: what is the estimated time of failure?

The work focuses on the analysis of time-series of ground-based photographs acquired on four landslides in the French Alps. Each dataset was used to develop and implement specific processing tools. Three cameras were setup in the field in order to monitor the long-term dynamics of the unstable slopes in their pre-failure stages; one camera was setup on an unstable slope in its post-failure stage and in a monitoring context.

Research context

The research has been supported by the CNRS/Institut de Physique du Globe de Strasbourg (IPGS) and the SAGE company (Société Alpine de GEotechnique) as part of a CIFRE/ANRT1 research contract.

During the thesis, stays at the company (up to 20 percent per year were done) in order to collect data, carry out additional measurements on site and test and optimize the analysis tools.

Some words on the SAGE company...

The SAGE company carries out missions in relation to mountain development, public works and natural risks. This last field, which is essentially based on the assessment of gravitational risks, integrates studies ranging from the geometric slope characterization to the modelling of its evolution, and the design of mitigation measures. For this purpose, the company uses classic in-situ monitoring methods such as extensometers or topographic targets measured by automatic total stations. SAGE also contributes to the analysis of block trajectories, field reconnaissance and modelling. The detection and survey of gravitational instabilities at regional scales are complementary missions. SAGE is continuously looking for new technologies aiming at meeting these expectations and thus improve the understanding of gravitational processes and assist the decision-making of the local authorities.

Objectives

The objectives of this work are to design, test and implement automated image processing pipelines for the analysis of monoscopic/stereoscopic image time series from fixed terrestrial optical cameras.

Low-cost camera system prototypes (e.g. cameras, digitizers, telemetry system, optimized protection boxes) have also been developed for standard deployment on the site. The processing pipeline is built around modules from the selection of the image sequences, to post-processing of change detection and displacement/velocity fields (Figure 4). Providing a versatile, flexible, robust and automated tool for the image time series analysis and interpretation is a real challenge in operational contexts since it is necessary to retrieve a maximum amount of information in short time.

Several issues and questions can be considered:

• What kind of useful information could the analysis of image time series provide?

• Is information retrieved from image time series sufficient to document gravitational motion, especially in the case of monitoring activities?

• How to evaluate the relevance (completeness, resolution, accuracy) of this source of information in the absence in-situ/ground control data?

To answer these questions, we proposed to:

• Develop, implement and test a processing pipeline for the analysis of terrestrial optical image time series acquired from fixed cameras. The processing is based on the open-source photogrammetric library MicMac2 , MPIC3 , COREGIS4 and Python libraries.

• Develop and implement a method to reconstruct the 3D displacement fields from pair-wise cameras considering the correlation results and the camera orientations.

• Develop and implement tools for the thickness estimation of the moving mass.

• Implement tools for the analysis of the landslide mechanical behavior (deformation, rheology).

• Implement post-processing analysis tools for the modelling of slope behaviour (rainfall-velocity relationships, landslide kinematic regime, time to failure forecasts).

Structure of the thesis

The thesis is organized in five chapters.

The first chapter presents the morphology and kinematics of the landslide cases used to develop and test the processing strategy. It further reviews various methods available to study and understand these phenomena, distinguishing between sensors specific to in-situ or remote measurements. An exhaustive state-of-the-art discussion on time-lapse photography is then provided. Finally, the four unstable slopes are presented from geological, kinematic and monitoring points of view.

The second chapter introduces the TSM pipeline (Tracing Surface Motion) (Desrues et al., 2019) developed within the thesis. It is based on image correlation techniques to calculate both qualitative (ground surface changes) and quantitative (displacement/velocity fields) information. The method is modular and associates tools for (1) the selection of the image sequence to process from the initial raw image time series, (2) the co-registration of the image stacks and the correction of the camera movements over time, and (3) the calculation of the terrain motion using change detection approaches. TSM is based on the open-source photogrammetric library MicMac and tailored for the processing of monoscopic images. A sensitivity analysis is conducted by testing four categories of parameters: the image modality, the image matching parameters, the size of the stable area used in the co-registration, and the strategy used to combine the images in the time series. The application of TSM on two case studies -e.g. the Chambon landslide characterized by slow motion (< 10 cm.day -1 ) and the Pas de l'Ours landslide characterized by moderate motion (> 50 cm.day -1 )-demonstrates the generalised applicability of the method. The results provide information on the kinematics and the spatial behavior of both landslides.

The third chapter presents the application of TSM for image time series acquired with pair-wise (stereoscopic) cameras in order to combine the displacement/velocity fields from two different points of view. We propose a technique to reconstruct the 3D displacements from the displacement grids, the correlation coefficients and the camera orientations. This new method allows us to minimize the dispersion of the displacement estimates in the case of the stereoscopic approach. This approach was applied to the slow-moving Montgombert landslide. In addition to this mechanical analysis, we used the constructed 3D displacement time series to estimate the thickness of the moving mass, considering a pure translational slide. The approach is based on the mass conservation law. Applied to the Montgombert landslide, it indicates a thickness estimated to a few meters. This result is consistent with previous geological and geotechnical observations.

Introduction

The fourth chapter presents a multi-sensor analysis combining geological observations, meteorological data, seismicity, topographic measurements and simple physical modelling applied to the Cliets rockslide (Desrues et al., 2020), submitted. We implement post-processing tools to analyse the relationship between displacements and rainfall and to perform kinematics and mechanical analyses. The mechanical analysis is carried out with a frictional model (Helmstetter et al., 2004) which allows us to define the velocity of the sliding regime. We also implement a tool to estimate the time to failure (Voight, 1989).

Additionally, we present the application of TSM to image time series acquired with a monoscopic camera during the pre-failure and the post-event stages.

The fifth chapter proposes a general conclusion and presents some technological and operational perspectives on the use of terrestrial image time series and toolbox and the use of image time series for geohazard analysis.

Chapter 1

Remote sensing methods for the monitoring of gravitational movements

The purpose of this chapter is to introduce time-lapse photography in the context of monitoring gravitational movements. Surface displacement and velocity are characteristic features that allow the monitoring of an instability evolution and the determination of its state of activity. To measure these features, several techniques exist, being either in-situ or remote sensing. In the context of monitoring at the scale of a full slope, time-lapse photography is a complement to existing methods allowing increasing the spatial density and the temporal frequency of data. We first present the different types of gravitational movements that one may encounter, with regard to their kinematic characteristics and, then, review the different monitoring methods, in particular introducing the time-lapse photography method and its challenges. Finally, we present the study sites for which the time-lapse photographs have been acquired during this work.

Notations: The name 'digital elevation model' (or its acronym: DEM) includes, by default in this chapter, digital surface and terrain models. Additionally, MPix will here stand for MegaPixels.

Chapter 1

Remote sensing methods for the monitoring of gravitational instabilities 1 Gravitational movements: risk and challenges

Types of gravitational movements

The term gravitational movements, also commonly referred as landslides, defines a diversity of processes implying a moving mass, caused by the action of gravity. Various types of geological materials can be at play, such as soil, bedrock or a mix of both, in a specific movement type (falls, slides, flows...).

Landslides may be distinguished by their material, geometry, velocity or volume. Therefore, finding a simple classification for these different processes has become necessary. Over time, several classifications were proposed to define the process with more precise vocabulary, thanks to a growing knowledge on the subject. In 1978, the first classification is proposed by Varnes (1978), in which movements are defined according to the type of motion and the type of material. From this starting point, several improvements were especially made by Hutchinson, J.N. (1988) and later by Cruden and Varnes (1996) who append the ground velocity along the slope as a new discriminating criteria. Other classifications where proposed over time to adapt to different purposes such as engineering geological mapping (Nemčok et al., 1972;Hungr et al., 2001), loess mapping (Li and Mo, 2019) or slope instability mapping (Smith and Dixon, 2015).

In the monitoring context, the most appropriate descriptor for the characterization of landslides is the surface velocity over time. Indeed, it provides information on the state of slope stability. Figure 1.1 provides a schematic overview of the classification by ground velocity. It clearly shows that a fast movement cannot be monitored in the same way as a slow movement. In other words, before monitoring, the frequency and the resolution of the measurement acquisition need to be assessed, according to the type of instability.

Based on the Varnes classification, the most commonly encountered types of instabilities according to the velocity descriptor are summarized in Figure 1.2 (Cruden and Varnes, 1996). They are listed hereafter:

• [Extremely rapid velocity: ≥ 5 m/sec] Falls correspond to fast-moving events involving single blocks. These blocks are detached from a steep slope on discontinuities, such as joints. They are characterized by a stage of free fall, which may be accompanied by rolling or bouncing, depending on the slope. The size of typical falling blocks will depend on the geology at play, which includes the bedding thickness, the dip direction and the joint orientations. Falls can be preceded by toppling or sliding motion.

Flows occur when the ground reaches its liquid limit (Hungr et al., 2001;Malet et al., 2003).

There are three types of flows. Debris flows are moving fast (Hungr et al., 2001) and load soil, rock and water as a matrix. Debris avalanches are similar to debris flows except that they more fast and involve very large volumes of material. Earthflows are characterized by saturated material involving fine soils. Varnes (1978) and Cruden and Varnes (1996) classification. The illustrations are from the United States Geological Survey and are available online (USGS, 2004). DSGSD stands for Deep-Seated Gravitational Slope Deformation and defines sliding movements that affect large slopes over a long period of time. The camera symbol indicates the velocity interval in which time-lapse photography monitoring is most appropriate.

Chapter 1

Remote sensing methods for the monitoring of gravitational instabilities

• [Very rapid velocity: ≥ 3 m/min and ≤ 5 m/sec] Topple is the rotation of a rock mass around an axis below its center of gravity. This movement usually takes place towards the front of the slope (Tamrakar et al., 2002;Nichol and Wong, 2007).

This type of movement generally leads to the subsequent fall of the blocks and echoes the movement presented above.

• [Rapid velocity: ≥ 1.8 m/hr and ≤ 3 m/min] Slides are characterized by a downslope movement along a sliding plane called 'failure surface'

or 'surface of rupture' (Petley, 2004). They appear when the stabilizing forces, defined by the shear strength and weight on the sliding surface are exceeded by driving forces including gravity.

We distinguish three types of slides. Rotational landslides have a surface of rupture that draws a concave curve. Thus, the mass revolves around an axis parallel to the slope and remains, more or less, a single coherent mass. It is essentially characterized by the presence of a well-defined scarp at the crown (the top of the landslide) and the presence of a zone of depletion and a zone of accumulation. Some cracks may be observed before the main scarp is formed. A rotational landslide mainly occurs in homogeneous materials. Translational slides have a surface of rupture similar to a roughly planar surface. This surface can be affected by a slight rotational movement and it corresponds usually to a discontinuity such as a fault or a joint. Block slides are translational slides, in which the material is moving as a single mass.

• [Moderate velocity: ≥ 13 m/month and ≤ 1.8 m/hr] Spreads occur on a gentle slope or a flat surface in response to soil liquefaction (Youd, 2018). They are defined by an extension of the soil, with a general subsidence of the fractured mass of cohesive material into softer underlying material.

• [Slow to extremely slow velocity: ≤ 13 m/month] Creep is a deformation supplied by internal shear stress which is not sufficient to reach rupture.

Curved tree trunks are good creep indicators since the base of the tree is moving with the ground soil whereas the top is trying to grow straight up.

All of these processes can be marked by the presence of simultaneous and different types of movements in several slope areas, and/or by a change of type of movement over time. In other words, considering more complex landslides, the environment will not, always and everywhere, evolve at the same ground velocity. Thus, monitoring these gravitational phenomena requires us to be able to detect and track several ranges of velocity both in time and space.

Controlling and triggering factors

The areas prone to landslides are predominantly characterized by controlling factors that are directly part of the environment of study, whether at the scale of the landslide, that of the slope or of the valley. They consist essentially of the geological structures and the hydrological context. The most common natural indicators, to estimate the susceptibility of landslides, are the slope angle and the slope orientation Chapter 1

Remote sensing methods for the monitoring of gravitational instabilities (Donnarumma et al., 2013;Çellek, 2020). The slope angle expresses the surface steepness according to the horizontal and it ranges from 0 • (horizontal surface) to 90 • (vertical surface). It can be computed from Digital Elevation Model (DEM). The angle depends on the local structures and on the lithology (e.g., presence of joints, of discontinuities or of bedding planes). Several studies show that soil landslides mainly occur on gentle slopes, whose angle is comprised between 10 • to 30 • (Zêzere et al., 1999;Iwahashi et al., 2003;Çellek, 2020). Additionally, Beaty (1956) suggests that landslides mainly occur on slopes with the least exposure to sunlight and where temperatures are low, based on the study of the orientations of landslides that have occurred in the Berkeley area (California). 70% of landslides occurred on slopes with both northerly and easterly orientation. This suggests that the areas most prone to landslides are those where there is more moisture and less evaporation, as confirmed by Rotaru et al. (2007). In addition to the slope steepness and orientation, the mechanical properties of the rocks define the stability considering cohesion and angle of friction (Rotaru et al., 2007).

Gravitational movements may be induced or at least disturbed by triggering factors that are related to the climate (i.e., rainfall, snowmelt), or to tectonics and seismicity (Dou et al., 2019). Intense rainfall or snowmelt imply change and/or saturation of the ground-water levels (Zêzere et al., 1999;Polemio and Petrucci, 2000;Jesus et al., 2018), affecting the mechanical process by increasing the rock or soil movement, and leading to a catastrophic scenario. In addition, on a longer time scale, rock alteration and erosion are also natural triggers. Earthquakes and active tectonics also play a role in the destabilization of masses. When the ground shakes, soil materials are subject to dilatation, which allows the infiltration of water into the environment (Rodrıguez et al., 1999). Materials then lose contact with one another involving the loss of frictional forces.

Additionally, the slope environment may be disturbed by the removal of vegetation since it partly regulates the amount of water in the soil and by the human activity since a slope profile is usually modified when there are work-sites nearby (e.g. the construction of roads, houses or transportation corridors) (Baioni, 2011). The human activity may consist of excavation, overloading or fill for instance.

This has the effect of disrupting the equilibrium between the shear force and the shear resistance along a potential failure surface (Rotaru et al., 2007). 2 Landslide monitoring: in-situ measurements and remote sensing

Monitoring a landslide, either in a event context or in a post-event context (Figure 1.3), is distinct from a single landslide survey since the measurement is acquired repeatedly or continuously over time (see Figure 1.4). A survey is a punctual measurement such as drilling (from geotechnical measurements) or a geophysical profile, while monitoring uses techniques that are repeated at more or less regular intervals.

It is often supplemented by field observations. The technique depends on the type of movement, on the size of the object, but also on the possible risks involved. The required precision and the time scale of the measurements are associated with the study, even if the quality and quantity of the measurements depend indirectly on economic constraints.

Selection of the monitoring techniques

We divide the monitoring techniques into two main categories: the in-situ techniques and the remote sensing techniques. These two techniques may involve different platforms (satellite, aerial and terrestrial platforms). The main difference between these two categories is that the first, in contrast to the second, requires a direct contact between the measuring device and the ground surface. Several sensors are presented in Table 1.1 according to their characteristics, such as the type of the measurement that can be obtained (e.g., displacements, DEM comparison) and the precision of the measurement.

Several studies list questions and recap the advantages and disadvantages of the monitoring techniques such as (Malet et al., 2002;Delacourt et al., 2007;Michoud et al., 2010). Stumpf (2013) summarizes in graphical form the most relevant techniques to be used in the case of a study, according to three criteria: the scale of the measurement (point, local or regional), the deformation rate (from extremely slow to extremely rapid) and the type of movement (slide, rockfall, spread ...). For instance, in the case of a rotational landslide, all monitoring techniques are appropriate, whereas for cliff monitoring, the proper techniques are reduced to laser scanning measurements, ground-based radar interferometry (GB-InSAR) and total station measurements. In fact, all aerial platforms are probably the least suitable for cliff monitoring, whatever the technique, because of the line of sight of the measurements which is similar to the main slope aspect (i.e., vertical direction). For a local measurement, the best platforms will be the total station, but also terrestrial photogrammetry and the GB-InSAR, while, for a regional measurement, the satellite platforms will be the best as well as aerial photographic acquisition or LiDAR acquisition.

We can also include, in this list, the frequency and the duration of acquisition with respect to the monitoring stage, which is specific to each instrument: the duration of the acquisition varies from a few seconds (case of a time lapse of images) to several minutes / hours (case of a LiDAR acquisition). In addition, the frequency of the acquisition depends strongly on the type of the platform which is used: satellites have a revisit time of one day minimum while that of terrestrial fixed cameras is of a few minutes or even seconds.

We propose to summarize this list, with these main questions:

• sensor parameters

Chapter 1

Remote sensing methods for the monitoring of gravitational instabilities What measurement accuracy is required?

What is the spatial coverage of the measure?

What is the distance to the object?

What is the temporal and spatial resolution of the measurement?

What is the cost associated to the survey?

• slope parameters

What is the purpose of the measurement?

What types of products are requested?

What is the expected ground deformation rate?

What is the pattern of motion that characterizes the gravitational movement?

• data parameters How reliable is the data?

What is the response time of the measurement?

What is the format of the data storage / data transmission?

Does the data support semi-automatic or even automatic processing?

In-situ measurements

In-situ techniques provide either absolute or relative measurements. Relative measurements between two points are usually distance measurements and may come from an extensometer, crack opening or inclinometer sensors (Corominas et al., 2000;Klimeš et al., 2012). As for absolute measurements, they correspond to angle and distance measurements or directly to 3D coordinate positions. They are also called geodetic measurements.

GPS/GNSS survey

GPS (Global Positioning System) / GNSS (Global Navigation Satellite System) are satellite-based navigation and positioning systems, which provide precise 3D coordinates (with an accuracy that can reach a few millimeters). This measuring system is increasingly used in civil applications. It is suitable both for monitoring slow deformations, such as tectonic movements (Nilforoushan et al., 2003;Devoti et al., 2011) or subsidence Engelkemeir et al. (2010), and for faster deformations, such as landslides (Gili et al., 2000;Malet et al., 2002;Squarzoni et al., 2005;Benoit et al., 2015). GPS/GNSS are often used repeatedly and not continuously over long periods of time mainly because of the cost of such a system and the bias in multi-paths (when the direct path is blocked by trees, buildings...) (Malet et al., 2002). The GPS/GNSS measurement can be carried out in real time (Real Time Kinematic -RTK) as in (Malet et al., 2002). (d) Principle of the direct levelling form (Schultz, 1995).

teletransmission (Leick et al., 2015;Hofmann-Wellenhof et al., 2012). In the latter case, it is not necessary to have a direct line of sight between the GPS base and the measurement point. This flexibility of the measurement is particularly an asset in measurement campaigns conducted in the natural environment.

Another advantage is that it is possible to measure whatever the weather conditions, and whatever the time at which the measurement is made. However, the orbit of the satellites, the ionospheric delays and the environment of the point will influence the quality of the measurement (Gili et al., 2000;Squarzoni et al., 2005).

Tacheometry / Theodolite

A tacheometer is a device used to measure horizontal and vertical angles between two targets. A theodolite, for its part, only measures angles. When the instrument is robotized, the term total station is used. Monitoring landslides from total station measurements requires a dense topographic network, composed of a set of artificial targets, to cover the maximum amount of the area. Additionally, some targets are set up in an assumed stable area to be used for the instrument calibration. The spatial coordinates of each target, expressed in a local reference frame, are deduced from the angle measurements and are, mostly, the deliverable. The angular precision ranges from 1 to 7 . The benefit of such a device is that it can acquire data night and day, in rain and sunshine. Atmospheric conditions may bias the quality of the measurements. As with any in-situ measurement, the total station measurement depends on the ground evolution. Indeed, the displacement of the ground can cause the topographic target to fall or tilt, which then causes misalignment with respect to the total station. As a result, some reliable information may be irreversibly lost.

Precise levelling

Precision levelling is an accurate method of differential levelling that uses precise levels and a procedure of rigorous observation (Ferhat et al., 2017). This method aims to achieve high orders of accuracy (Malet et al., 2002). It is mainly used to monitor areas in which vertical movements are significant (e.g., areas of subsidence) [START_REF] Murase | Volcano-hydrothermal activity detected by precise levelling surveys at the tatun volcano group in northern taiwan during 2006-2013[END_REF]Sabuncu and Ozener, 2014;Ferhat et al., 2015).

Extensometry

An extensometer is a monitoring device for a relative measurement of the distance between at least two points. A displacement measured by an extensometer does not provide information on the vertical and horizontal components of the movement or information on the thickness of the moving masses. It is a global measure providing information about the ground deformation (Angeli et al., 2000). There are several kind of sensors (Michoud et al., 2010): the probe extensometers which are non-permanent measurements, the fixed borehole extensometers which measure the change of distance of two or more specific points along the axis of a borehole, and the wire extensometers (Corominas et al., 2000) which measure the distance between two points; generally with one point located in the unstable part and one other located on a stable part.

Remote sensing imaging methods

Remote sensing is a set of techniques which detect or monitor the ground features as well as the physical characteristics of an area at a distance. Several platforms may be used: satellite, aerial and ground-based methods. The surface area covered by these techniques varies from a few m 2 to several hundred km 2 . The accuracy of the measurements ranges from few millimeters to meters (Delacourt et al., 2007). The active techniques are based on the emitted-transmitted principle by providing their own energy to illuminate an object. They send energy and receive the response which is reflected by the object. The measure is therefore performed in 1D along the line of sight of the sensor. The operator has control over the characteristics of the emitted signal, such as the wavelength which, depending on its value, may or may not pass through the vegetation. The passive techniques may also be used. They do not interact with the object and detect natural emitted or reflected sunlight radiation. There are part of the group based on optical sensors and provide a 2D measurement, or a 3D one if some additional data are available.

Active imaging sensors

We here introduce two active sensors, which are the most often used in landslide monitoring:

Chapter 1

Remote sensing methods for the monitoring of gravitational instabilities Satellite InSAR, UAVSAR and ground-based SAR interferometry Among passive sensors are radar techniques based on the Synthetic Aperture Radar (SAR) data. Radar is an active and coherent source that emits electromagnetic signals in the [1][2][3][4][5][6][7][8][9][10][11][12] GHz frequency band. Satellite SAR images have regularly been recorded since the 1900s, enabling the monitoring of ground motion. They are used in several applications, such as the measurement of the ground displacements induced by earthquakes, landslides or glaciers (Massonnet et al., 1994;Bawden et al., 2001;Bamber and Rivera, 2007;Nikolaeva et al., 2014;Modeste, 2020). One of the major advantages of this method is the ability to acquire images day and night.

To measure the ground deformation and provide maps of surface deformation, images are acquired repeatedly. This technique, called SAR interferometry, consists in measuring ground motion occurring between two acquisitions by providing the motion as a complex signal expressed by fringes (Rosen et al., 2000;Pathier et al., 2003). The principle of a radar acquisition is shown in Figure 1.6 and explained below.

The radar follows a trajectory, called azimuth, which is perpendicular to the line of sight (LOS) in which the measurement is made. The footprint of the LOS extends from near ranges to far ranges covering areas at the regional scale. Yet, some distortions can be introduced on SAR images, which are caused by some low coherence measure and which may correspond to important vegetation change, snowfall or significant ground surface change. Additionally, in a mountainous context, shadow zones relative to the orientation of the hills in particular can be inserted; the radar method is sensitive to relief.

The Uninhabited Aerial Vehicle Synthetic Aperture Radar (UAVSAR) airborne is a NASA repeat-pass SAR interferometry system. This instrument enables the acquisition of L-band SAR images at various time-scales and different optical imaging geometries (Hensley et al., 2009;Delbridge et al., 2016). The submeter-level resolution (0.6 m in azimuth and 1.9 m in range) allows the imaging of geological structures on small scales (Hensley et al., 2009;Hu et al., 2020a). This technique was successfully used by combining several look directions, for landslide monitoring (Delbridge et al., 2016;Hu et al., 2020a,b) or glacier flow monitoring (Minchew et al., 2015).

The ground-based interferometry (GBInSAR) is often used to monitor the ground in near real-time with a millimetre level of accuracy (Antonello et al., 2004;Bozzano et al., 2011). It covers a local scale, and thus, a smaller area. It is a continuous measurement that does not require the installation of targets. In contrast to the satellite case, this technique is flexible in its application, since the operator controls all the acquisition parameters. However, the instrument needs to be set up on a stable area and to be fixed during the whole acquisition. Thus, the costs related to this method, can be substantial. Chapter 1

Remote sensing methods for the monitoring of gravitational instabilities object : d = (c * δt)/2. 3D coordinates are then recorded as well as the intensity of the pulse and other signal attributes. The instrument can reach a very long range up to 4000 m and can also be embarked in airborne platforms such as airplane or helicopter (Aerial laser scanning -ALS) (Ventura et al., 2011), in mobile platforms such as vehicle or boat (Michoud et al., 2015) or in terrestrial platform (Terrestrial laser scanning -TLS) (Royán et al., 2014;Kromer et al., 2017). The resolution of the measure varies with the distance between the laser and the object, and with the angular spacing between two scanned points. The quality of the measurement can be affected by the atmospheric conditions as well as the ambient humidity.

The technique is used from the laboratory scale to reconstruct rock samples (Wang et al., 2019) to the macroscopic scale to reconstruct digital terrain models (Travelletti et al., 2008;Prokop and Panholzer, 2009) or to monitor rockfalls (Royán et al., 2014;Kromer et al., 2015). Indeed, with a single scan, some slope information such as geometry and orientation are directly available. The 3D model can also be inserted in rockfall simulations software to provide a ground surface on which block trajectories are calculated. (Ondercin et al., 2014;Noël et al., 2018;Sala et al., 2019). This method is used in several fields such as archaeology (Fernandez-Diaz et al., 2014), geosciences (Jaboyedoff et al., 2009;Eitel et al., 2016), construction and rail environment (Assali, 2015;Kedzierski and Fryskowska, 2014) or even in the automobile industry, in order to detect how far away obstacles are from a car (Rasshofer and Gresser, 2005).

To monitor slope changes, multi-temporal laser scanning is often used, either by repeating the measure with an operator or automatically (Jaillet et al., 2011;Kromer et al., 2017;Point et al., 2018). The frequency of acquisition varies from one every few hours to one every days, months, years ( (Eitel et al., 2016)). It depends on the ground velocity of the movement as well as on the objective of the measurement campaign.

As the measure is realized in the light of sight of the instrument, several viewpoints, in the case of a terrestrial acquisition, are necessary in order to avoid occluded areas. The resolution is then increased and spatially heterogeneous. When several point clouds are generated from one or more viewpoints, they are aligned with each other to obtain a single entity. Many algorithms can be used to process point clouds.

The Chapter 1

Remote sensing methods for the monitoring of gravitational instabilities

Optical satellite images Multi-spectral and panchromatic images are, nowadays, acquired with satellite optics with a very high resolution (VHR) compared to the previous decades. The first optical satellites for civilian purposes, SPOT and LANDSAT, were launched in the 1990 s. They are characterised by decametric resolutions and a revisit time of 18 to 30 days respectively. Since then, Quickbird and Ikonos, high-resolution satellites with near-metre resolution (HR), were launched in the 2000 s followed by the Sentinel-2. The resolution of this latest satellite ranges from 10 to 60 m and depends on the spectral band used. Other satellites have been launched since then, enabling even finer resolution to be achieved: 2.5 m for the SPOT satellites and 0.7 m for the Pléiades constellation with a revisit time of 1 day.

Innovations in sensor technology and processing algorithms (Section 3.2.2) have made it possible to use these high and very high resolution images for qualitative and quantitative analysis of surface changes. By improving processing algorithms and developing new ones (e.g., MPIC technique (Stumpf et al., 2017)), change mapping and displacement rates analysis are realized faster with higher confidence. The recent development of new acquisition methods such as stereoscopic or tristereoscopic acquisition has opened up new perspectives. Indeed, thanks to these acquisitions high-resolution DEMs can be generated. Optical satellite images can be found in several applications such as the detection of ground motion as well as that of avalanches (Larsen et al., 2013), earthquakes (Lacroix et al., 2015), landslides (Delacourt et al., 2004;Stumpf et al., 2017Stumpf et al., , 2018;;Desrues et al., 2019a), or the estimation of glacier speed (Heid and Kääb, 2012b). They can also provide a general overview on natural hazards highlighting the consequences of natural disasters by generating DEMs and orthoimages (Shean et al., 2016).

Optical aerial images They are acquired either from a plane, an helicopter or even a paraglider or from a Unmanned Aerial Vehicle (UAV) such as, kites or balloons (Colomina and Molina, 2014). The embedded camera generally points downwards and the images are acquired along the flight path. They are often used in the fields of agriculture, archaeology or in the humanitarian sector (Vigneau et al., 2014;Barge and Régagnon, 2017;Apvrille et al., 2017). UAVs are increasingly being used to map mountainous surface, which, at the same time, allows the human operator to be independent and autonomous in the field. The most common device is the low-cost motorized and piloted UAV, also called drone. This device has the advantage of being very easy to handle and simple in its use. The maximum range can reach 8 km, which is the case for the drone DJI Mavic Pro 2, and the images are acquired in ultra high resolution (>18 MPix). Large areas can be covered (up to several square kilometers) by aerial images with high resolution. It can reach a few centimeters for lower areas (of about 1x1 km 2 ) (Casson et al., 2005).

Optical aerial images are processed by stereo photogrammetry or by Structure from Motion (SfM), which enables the generation of a large-scale digital elevation model. DEMs may achieve a sub-meter accuracy for airborne images (Dewitte et al., 2008), and may reach a resolution of a few centimeters for UAV images (Niethammer et al., 2012;[START_REF] Westoby | structure-frommotion'photogrammetry: A low-cost, effective tool for geoscience applications[END_REF]Lucieer et al., 2014;Gupta and Shukla, 2018).

To provide DEM, ground control points (GCPs), or the position of the cameras, can be introduced in the SfM algorithm in order to convert the model built in a local coordinate system to a geodesic coordinate system. The addition of GCPs improves the quality of the 3D reconstruction by constraining the algorithm (bundle adjustment stage) and by georeferencing the 3D model (Bitelli et al., 2004;Fonstad et al., 2013).

Aerial images also provide extremely important databases that allow the retro-analysis of ground movements over several years, in order to better understand the current dynamics. Indeed, they provide the ability to detect small movements over a large period of time such as a decade (Hapke, 2005;Dewitte et al., 2008). One of the major limitations of the use of airborne images is the geometry of the study object. Indeed, photogrammetry fails when the object represents steep slopes or very marked slope breaks.

As for UAV images, the environment and the weather conditions strongly constrains the acquisition. For instance, most of the drones cannot fly when the wind speed reaches or exceeds 30 or 40 km/h. Beyond this threshold, it becomes more difficult to control the device, the battery loses autonomy and the drone risks crashing. In addition, the required image overlap is compromised. It is also impractical to fly the UAV in a forest environment or in a dense urban environment. In addition, the presence of electrical cables may add some additional constraint.

Optical terrestrial imaging Terrestrial images can be acquired by one or more cameras using two approaches, either the monoscopic approach or the stereoscopic or multi-view approach. Using the terrestrial platform offers the advantage of generating high-resolution DEMs but also of generating surface change maps qualitatively and quantitatively, at low cost with ease of use. It is often used in several domains, such as in archaeology to construct 3D model of fragile objects or inaccessible objects (Drap et al., 2007), in glacier surveillance (Kaufmann, 2012;Piermattei et al., 2015) or in landslide monitoring (Travelletti et al., 2012;Valentin, 2018). This method, used in a repetitive way, is hereafter called time-lapse photography and is presented in Section 3.

Time-lapse photography

In this section, we are interested in time-lapse photography carried out on land-based platforms. The different acquisition approaches are presented as well as the sensors and acquisition systems that can be used. The processing algorithms are also introduced as well as the post-processing algorithms.

State of art

As previously mentioned, images from time-lapse photography can be acquired in different modes: RGB, multi-spectral, etc. depending on the field of investigation and the desired products. Indeed, RGB images give access to quantitative and qualitative surface change detection (Travelletti et al., 2012;Stumpf et al., 2015) whereas infrared images give access to surface temperature maps allowing the detection of thermal anomalies (Frodella et al., 2017;Le Roy, 2020).

In the context of the thesis, we were only interested in images acquired in RGB mode.

Figure 1.9: Schematic table on the possible modes of time-lapse photography.

Principle of time-lapse photography

Time-lapse photography is an acquisition method that allows us to associate a temporal dimension to the photography data. From temporal data, that we will call time series, it is possible to track surface changes and to quantify the different types of deformations. This method can be broken down into two acquisition approaches: a monoscopic approach and a stereoscopic (or multi-view) approach (Figure 1.9).

The first approach consists in acquiring a succession of photographs of a moving object from a single point of view. The acquisition characteristics, such as the focal length, aperture or isometry, and the spatial coverage remain fixed during the entire image acquisition. As for the temporal scale of the measurements, it can vary from milliseconds (video acquisition) to days. This approach generates a set of photographs taken at regular time intervals, that enables changes to be detected and movements to be quantified using image matching methods (see section 3.2.2). Monoscopic time lapse was successfully used in volcanology to monitor dome deformations (Walter, 2011), and in glaciology and landslides (Ahn and Box, 2010;Fallourd et al., 2010;Gance et al., 2014;Schwalbe and Maas, 2017;Altena and Kääb, 2017;Marsy et al., 2020) to monitor surface displacements (Figure 1.10). In addition, Travelletti et al. (2012) have demonstrated the potential of a permanent time-lapse in the study of the Super-Sauze landslide. The displacements obtained by image analysis have a centimetre accuracy.

Like any passive measurement, time-lapse photography is dependent on weather conditions that can render the photographs unusable. Even if the influence of wind remains negligible, rain, snow and fog can obscure the view. The image acquisition requires a proper point of view, which depends directly on the study site. The camera can only capture the parts that are visible from its point of view. Nevertheless, this technique has the advantage of being low cost and easy to use, while offering the possibility of having a quick acquisition that can be applied at different scales in time and space.

Type of camera sensor

The choice of the camera sensor influences two main parameters: the resolution of the ground pixel (i.e., the pixel equivalent size in the real world system), which directly determines the resolution of the displacement maps, and the camera spatial coverage. This leads to the choice of the sensor size and the focal length. There are two types of sensors: CCD (Charge Coupled Device) and CMOS (Complementary metaloxide-semiconductor). The most commonly used, single lens reflex (SLR), cameras employ CMOS sensors. The size of the sensor is related to the size and the resolution of the image. The size of the image pixel is determined as:

δ pixel imagesize = Size sensor Resolution image (1.1)
The focal length of a camera varies from a few millimeters to more than 300 mm depending on the lens used. Fisheye lenses are those with a focal length less than 14 mm. This type of camera can be found in traffic surveillance, for example, and are commonly called 'webcams'. Then comes the ultra wide-angle lenses, whose focal length is between 14 mm and 20 mm, and the wide-range cameras with a focal length between 24 mm and 35 mm. Finally, the more standard lenses are those with a focal length between 45 mm and 70 mm. This range corresponds to the field of vision of the human eye. Above this range, lie the telephoto lenses, with a focal length of more than 300 mm. The focal length controls the view angle:

the longer the focal length, the narrower the angle of view. Thus a substantial advantage of this technique is the ability to use different fixed lenses at different focal lengths.

If we consider f the focal length, d the average camera-object distance, the ground pixel size is then defined as: The use of a very wide angle webcam may also imply strong distortions in the captured images (see Appendix A). To use this type of sensor in the quantification of ground motions, it is thus necessary to apply corrections to the photographs. This requires us to determine the distortion coefficients and, thus, to calibrate the webcam precisely. In the study of the Pas de l'Ours landslide, too few known points on the ground were available to calibrate the webcam. Calibrating a camera provides the opportunity to determine the internal parameters of the camera such as focal length and distortion coefficients, but also the external parameters that enable the camera to be oriented in a chosen geodetic reference frame.

δ pixel groundsize = d f * δ pixel imagesize (1.2)
Appendix A is dedicated to the description of these parameters.

Acquisition system and digitizer

The system that has proven to be the most suitable for a photography time-lapse acquisition is a standalone box in which the camera is placed, as well as a data logger powered by an external battery or by a solar panel. This low cost acquisition system allows us to control the acquisition frequency and its actual moment, according to the sun illumination. It also allows us to automatically upload images onto a server, if a WiFi connection is available in the area. Another advantage of this kind of system is its very low power consumption. Such a system allows to connect a very high resolution camera to the data logger. For instance, cameras with 21 MPix and 50 MPix were, respectively, used by Roncella et al.

(2014) and Kromer et al. (2019) to monitor slopes and rockfalls. This enables the object to be imaged accurately while being at a great distance (i.e., more than 500 m). Some other systems use a Raspberry

Pi as on-board computer, combined with the PiCamera module (Santise et al., 2017). However, the low resolution of the sensor (8 MPix), induces ambient noise in the image and constrains the distance between the object and the camera. In addition to these custom systems, some commercial solutions exist within an all-in-one box (e.g., Devisubox, Boxcam). They are robust to weather conditions but have two major drawbacks: high power consumption and a non-fixed focal length.

Figure 1.12 shows an example of two low cost acquisition systems developed by the SAGE company. In the first months of its installation at the Pas de l'Ours landslide (see Section 4), the captured images were first retrieved directly from the camera's memory card and after modifications of the system, they were Chapter 1

Remote sensing methods for the monitoring of gravitational instabilities sent directly to a server.

Methodology

To extract the information related to a surface change, several image processing algorithms specific to the two acquisition approaches (monoscopic and multi-view approach) are available. The three main steps common to both approaches are: the pre-processing (image selection, camera calibration, ....), the processing to extract the main information such as the image correlation of the DEM generation, and the analysis.

Image pre-processing

Image pre-processing consists of selecting the images to be further processed, transforming them, if necessary, to increase the signal-to-noise ratio and calibrating and determining the internal and external parameters of the camera. Image selection refers to the estimation of the quality of the images which may be altered by the presence of clouds, rain, fog or vegetation. The image quality needs to be as high as possible to avoid the integration of bias during processing. This operation is often carried out manually.

Biases can also be reduced by transforming the images to highlight their texture. The images can be transformed, for example, into grayscale or radiometrically enhanced (transformation on the histogram).

The internal and external parameters of the camera are involved in image processing at several levels: correction of distortion, transformation of image coordinates into geodetic coordinates, orientation of the camera in relation to the object, etc. This step does not affect the image data set as such, but is necessary for pre-processing since it affects the image acquisition as well as the image processing.

Image processing

Image processing depends on the adopted acquisition approach. But in both cases we can have a temporal series of images that can be processed in a similar way (Figure 1.9). For the multi-view approach, the images are processed by stereo-photography or by SfM in order to generate DEMs (Eltner et al., 2017;Parente et al., 2019). The projected DEM is then converted into field data such as shading or roughness, in order to texture the image (Lucieer et al., 2014;Travelletti et al., 2014). Therefore, image matching techniques can be applied to quantify the ground motions.

To track surface changes from two images, taken before and after the ground movement, several algorithms can be used to identify and quantify the displacements. Matching techniques can be differentiated in local matching such as tracking feature methods and in global matching such as image correlation methods.

Other emerging technical matching methods can be used to generate displacement fields which is notably the case with optical flow techniques such as the Lucas-Kanade method and its derivatives (Brox et al., 2004). They assume that the illumination properties of the object are preserved between two dates and that the movements are small.

Chapter 1

Remote sensing methods for the monitoring of gravitational instabilities

Tracking features methods detect time-invariant points (Harris et al., 1988;Lowe, 2004) which are discrete features that can either be natural (blocks) or artificial (topographic targets). There are several detection algorithms based on these specific points, such as the Harris corner detector (Gance et al., 2014) or the Scale Invariant Features Transform (SIFT) (Le et al., 2011). This method has the advantage of being robust and not dependent, for instance, on radiometric changes in objects. Figure 1.13a) shows an example of target tracking whose method is based on target properties (size, geometric luminance) (Gance et al., 2014). This method is limited only to the points which are detected by the algorithm and, therefore, the displacements are only calculated locally. The amount of information is then dependent on the study site and its identifiable landmarks. This is why many researchers use the image correlation technique which allows us to obtain a dense spatial deformation field. This technique has been tested and approved in many cases such as the analysis of glacier movements (Tiwari et al., 2014;Dehecq et al., 2015) or surface deformation (Delacourt et al., 2004;Travelletti et al., 2012;Stumpf et al., 2017).

Image correlation techniques can be used either to detect equivalent points in a pair of images, or to calculate the displacement field between two consecutive images. It is a technique that is sensitive to the displacements along the two spatial directions that are perpendicular to the line of sight of the camera. Let us consider two images of which there is a 'reference' one called here 'master' and there is the one being compared called here 'slave' (see Figure 1.13b)). To compute the displacement of a pixel (u m ,v m ), we center a correlation window, here denoted Cw m , of defined dimensions (nxp) in the master image, and a correlation window Cw s (of the same dimensions as Cw m ) in the slave image. This second window is moved in the slave image and, at each position, the correlation with the first window is calculated. The pixel with the highest correlation peak is considered to be the new position of the considered master location. The correlation results are composed of three components: the row (dv) and column (du) displacements, and the associated correlation coefficient. Figure 1.14 shows an example of correlation results applied to an image time series acquired in the landslide of Super-Sauze.

Correlation calculation can be performed either in the spatial domain or in the frequency domain. To provide a spatial distribution of the displacement field, various algorithms can be used according to the correlation calculation domain (Heid and Kääb, 2012b):

• The Normalized Cross Correlation (NCC) operates in the spatial domain. This algorithm is mostly used because it is not sensitive to high image intensities, and images with strong local anomalies can be better compared (Heid and Kääb, 2012b;Pham et al., 2014). As the correlation coefficient is normalized, the results from several image correlations can also be compared. This coefficient ranges from 0 (anti-correlation) to 1 (perfect match) and it is defined according to Hild and Roux (2012) by Equation ( 1.3): As mentioned, Cw m and Cw s are the spatial average of the respective correlation windows.

NCC = ∑ n,p (Cw m (u m , v m ) -Cw m )(Cw s (u m + du, v m + dv) -Cw m ) ∑ n,p (Cw m (u m , v m ) -Cw m ) 2 ∑ n,p (Cw s (u m + du, v m + dv) -Cw s ) 2 (1.3)
The calculation time of this algorithm can be rather large because the convolution operation is performed in the spatial domain (Bickel et al., 2018). In addition, the algorithm performs badly when the transformation of a point from the master image to the slave image is a rotation (Travelletti et al., 2012;Bickel et al., 2018).

• The Cross Correlation (CC) operates in the frequency domain and it provides the best spatial representation, with a low noise level compared to the NCC algorithm (Bickel et al., 2018). It is defined as Equation (1.4):

CC = IFFT (F m (u m , v m )G * s (u m , v m )) (1.4)
where F m and G s are the Fast Fourier Transform (FFT) of the correlation windows of the master and the slave images, respectively. IFFT denotes the Inverse Fast Fourier Transform.

As the algorithm is not normalized, the highest intensity changes can not be ignored and can cause severe dis-matches. There are several ways to normalize in the Fourier domain and one of the most common methods is the one that considers only the phase information.

The Phase Correlation (PC) allows to normalize the CC algorithm (Brown, 1992). It is defined by Equation (1.5):

PC = IFFT ( F m (u m , v m )G * s (u m , v m ) |F m (u m , v m )G * s (u m , v m )| ) (1.5)
The maximum of the phase difference enables us to compute the displacement (Heid and Kääb, 2012b;Bickel et al., 2018).

The choice of the size of the correlation window is an important parameter that directly influences the signal-to-noise ratio, the precision as well as the accuracy of the displacements (Delacourt et al., 2007;Travelletti et al., 2012). The larger the correlation window, the better the signal-to-noise ratio and accuracy, but the lower the accuracy. To cope with its influence, sub-pixel techniques can be used. They are considered if the displacements to be detected are smaller than the size of the image pixel. In return, these methods are very sensitive to noise and image quality. These techniques have been developed and applied in different areas such as the analysis of co-seismic deformation (Leprince et al., 2007), of glacier motion (Heid and Kääb, 2012a) and of landslide motion (Travelletti et al., 2012;Delacourt et al., 2004;Stumpf et al., 2017). In case of a sub-pixel hierarchical approach, the physical size of the correlation window varies as well as the size of the window to be explored, commonly named the search window (Figure 1.15). The first steps are performed at low resolution and then at high resolution. The pixel with the highest correlation peak at low resolution becomes the center of the area of interest for the next steps (Travelletti et al., 2012).

Several comparisons were made on the different methods of correlation in order to determine which algorithm was the most robust, the fastest or the most accurate (Heid and Kääb, 2012b;Rosu et al., 2015;Bickel et al., 2018). It was found that, correlation operators in the frequency domain give very successful results when the correlation window is large (bigger than 32x32 pixels) (Leprince et al., 2007). These approaches are embedded in several libraries such as COSI Corr (Leprince et al., 2007) 1 , MEDICIS2 (Berthier et al., 2005) or MicMac3 (Stumpf et al., 2014(Stumpf et al., , 2017) ) or other software like DPIVSoft (Meunier and Leweke, 2003) 4 .

Image post-processing and analysis

Displacements obtained from image correlation provide access to qualitative information that enables the mapping of areas affected by movement, but also quantitative information on displacements, surface velocities and the deformation involved. Several image matching algorithms exist to improve the precision and accuracy of the displacements and thus increase the signal-to-noise ratio, based in particular on the redundancy of information (Stumpf et al., 2017). However, post-processing is often integrated into the overall process in order to eliminate false detection and decorrelation due to large surface changes. The most frequent filters include masks based on the correlation coefficient, filters related to the direction and magnitude of the displacement velocities (Stumpf et al., 2014) and low-pass filters to remove outliers (Heid and Kääb, 2012a). From there, the displacements can be used for mechanical analyses (Helmstetter et al., 2004), to evaluate the dimensions of the object (Booth et al., 2013) or to estimate time of failure (Voight, 1989).

In the case where the images come from the projection of DEMs, other types of information can be extracted from the time series of DEMs. Indeed, the comparison of two point clouds constructed at two different dates, allows us to highlight geomorphological changes, areas of accumulation and loss of materials and in general, stable and unstable areas. The comparison of two point clouds can be carried out according to several algorithms such as the DEM-to-DEM computation (Schürch et al., 2011) and the surface velocity differ from one to another. The dominant deformation process is sliding and the measured average velocity is either slow or moderate (Table 1.2). For each site, optical single-reflex cameras of high-resolution have been permanently installed. The frequency of the acquisitions varies from one site to another since it has been adjusted according to the illumination of the sun, the expected displacement rate and the experience gained (Table 1.3). In addition, two cameras were set up in the field at the Montgombert landslide in order to test and analyse the influence of the line of sight over the measurement and to develop a 3D algorithm to reconstruct the displacements in a world coordinate system directly from the correlation results (see Chapter 3). dataset acquired for the period 8 June to 18 July 2018 (22 images according to the selection criteria), is presented in this work to be in agreement with the ancillary dataset (tacheometer-derived displacement rates on benchmarks) used to assess the results.

Montgombert landslide (Savoie, France)

The Mongombert landslide is located in the Arly Valley (Savoie, France) on the right bank of the Arly River. The bedrock is composed of micaschists and shale screes make up the unconsolidated parts. Since 2016 several mudflows and blockfalls that can reach the departmental road are observed. The observations suggest that the landslide is strongly impacted by weather conditions especially after a event occurred following heavy rains in January 2018 (Mathy and Lorier, 2018). The landslide continues to be monitored by automatic theodolite, extensometer and time lapse photographs. The moving zone is 150 m long and 115 m wide but continues to spread and several geomorphological markers such as cracks and scarps are visible further up the slope (see Chapter 3).

Two cameras were installed as part of a development to reconstruct a 3D displacement field from a stereoscopic system (Figure 1.17c). The conditions required to set up such a system (ratio of the distance between the two cameras to the camera-object distance) could be verified. The two acquisition systems consist of an 18 MPix EOS 100D camera with a focal length of 50 mm, connected to a Paratronics data logger as used previously. Images are acquired every 3 hours. Such a short frequency is due to the fact that the camera is oriented East/Southeast and that being installed in a very narrow valley, the brightness drops rapidly.

Les Cliets rockslide (Savoie, France)

The Cliets rockslide is located in the Arly Valley (Savoie, France) upstream of a departmental road which goes through the 60 meter long tunnel of the Cliets (Figure 1.17d) (downstream of the Montgombert landslide on the right bank of the Arly River). Since a long time, this area undergoes frequent damages (Kasperski, 2008;Pothérat, 2005;Mathy and Lorier, 2013) and in February 2019, 10.000 m 3 of mass fell (Lorier et al., 2019). The rockslide is mainly composed of mica schist which stemmed from an old metamorphic folded detritic sedimentary series. This implies that the rock mass may be generally subject to toppling (see Section Chapter 4). The area extends over a length of about 72 m and a width of 80 m.

A camera Canon EOS 2000 of 24 MPix was used to acquired images at a frequency of 8 images per day.

The camera was set up in the field when significant changes in ground velocity have been measured by in-situ sensors. The focal length used in this case is 24 mm. The camera system was installed on 22 December 2018 and is still operational. Because of mist, rain and snowfall, only a few images were available during the pre-failure stage. From 22 December 2018 to 9 February 2019, date on which the failure occurs, only 17 photographs were able to be used. Following the event of February 2019, several topographic targets tracked by automatic theodolite are no longer available mainly because they have been swept away by unstable masses. The contribution of time-lapse photogrammetry is therefore important in these cases as a support and tool for post-event monitoring.

Conclusions

Time-lapse photography is a tool that is starting to be used operationally as a result of various promising photogrammetric campaigns. It has emerged as a complementary measure to classical total station monitoring thanks to the extent of the spatial coverage and to its low cost price.

Working on the analysis of photographs and the understanding of gravitational movements called for the resolution of several questions during this thesis, on how to acquire the images, process them and extract the maximum amount of information from them. These scientific barriers involve several areas such as the management of data from its acquisition to its conservation and the scientific tools for data processing.

They can be summarized into several major and underlying issues:

• What is the influence of the instrument setting and parameterization on the results?

The acquisition parameters involve the acquisition frequency, the camera line of sight incidence angle, as well as the characteristics of the camera (focal length, aperture...), but also the format in which the image is recorded (raw or compressed format) and the format in which the image is processed (multi-band or single band). These parameters are important since they are set as soon as the camera is installed and subsequently condition the spatial coverage (i.e., the spatial limits of the interpretation) and the resolution of the measurements. How to control and anticipate these parameters? How to calibrate a camera or a couple of cameras once it is installed? In addition, what

is the advantage of one acquisition approach over the other? Is the accuracy improved with the stereoscopic approach compared to the monoscopic approach?

• What is the influence of the processing parameters on the results?

A time series of images can be processed in different correlation combinations, which can be sequential or relative to a unique reference image. These different correlation strategies are mostly derived from the processing of optical satellite images. But are the correlation methods of optical satellite images applicable to terrestrial images? Do we get the same results or do we lose some information with one technique or another? How can the signal-to-noise ratio be increased? Is the method of co-registration influenced by the percentage of stable pixels?

• What are the primary and secondary products that can be generated?

Depending on the acquisition approaches and parameters, several products can be calculated such as maps of metric displacement fields or deformation fields. What types of products can be derived from the different acquisition/processing combinations? What is the accuracy of these products?

• How far can we process data automatically?

Automating the processing as much as possible in order to avoid on-site and numerical interventions by an operator is one of the most complex challenges in our monitoring problem, which requires taking into account the storage conditions of the raw data as well as that of the results, their format and the deliverables/products in an operational context. Many scenarios need to be taken into account in order to provide the user with the best choice of parameters and to allow them to adjust them optimally. So... how far can we automate?

Throughout this thesis, we have tried to answer each of these questions, whether for the monoscopic (Chapter 2) or stereoscopic (Chapter 3) approach. The treatment chains and the influence of the different parameters are discussed and the products resulting from these treatments are presented. Chapter 4 is a foretaste of the limits of the optical method in the case of pre-event monitoring but nevertheless opens up an interesting perspective in the field of visual support when no other type of information is available.

The last question concerning the automation of processing is dealt with implicitly but will be summarised in Chapter 5.

Introduction

Terrain motion monitoring systems are currently integrating technologies to provide data on the space and time distribution of the displacement to complement the classical point-based measurements based on geodetic techniques (Peyret et al., 2008). Among them, multispectral imagery acquired from satellite (Marek et al., 2015), airborne (Nichol and Wong, 2007) or terrestrial (Niethammer et al., 2012) platforms can provide distributed information on the terrain motion. These techniques are useful in landslide research for terrain motion detection, mapping and monitoring (Gance et al., 2014). Terrestrial optical imagery is however seldom used for monitoring purposes, though high versatility and high spatial resolution of information can be accessed. With the development of relatively low-cost acquisition systems (very-high resolution consumer grade camera, 20 MPix digitizers), it becomes easier to generate long-and highfrequency image time series for documenting complex surface motion, and thus discovering specific and transient patterns such as the ones observed on landslides, ice/rock glaciers or volcanoes (Travelletti et al., 2012;Fallourd et al., 2010;Benoit et al., 2015;Hibert et al., 2015;Kromer et al., 2019). Simple permanent installations requiring only solar panels and a stable anchoring can be easily setup (Gance et al., 2014).

Geomorphological processes typically feature complex 3D deformation patterns and the displaced volumes can vary by several orders of magnitudes. Surface displacement rates and volumes of the unstable mass are first-order parameters to estimate dynamic properties (such as the rheology of the media), to quantify relationships with forcing conditions (meteorology, hydrology, earthquake) and identify thresholds for hazard assessment and warning criteria (van Asch et al., 2007). The most commonly used techniques for surface motion measurements are in-situ point-based sensors such as GNSS, inclinometers, extensometers and crack-meters (Malet et al., 2002;Calcaterra et al., 2012), while the most frequently used in-situ remote sensing techniques comprise total station measurements, ground-based interferometric synthetic aperture radar (GB-InSAR) and Terrestrial Laser Scanning (TLS) surveys (Travelletti et al., 2008;Prokop and Panholzer, 2009;Jaboyedoff et al., 2012). Several studies also demonstrated the possibility to estimate horizontal displacement measurements from remote platforms such as optical airborne/satellite images (Leprince et al., 2007;Debella-Gilo and Kääb, 2012) and more recently from SAR amplitude images using offset-tracking techniques (Raucoules et al., 2013). Currently, those measurements are still of limited use for the integration in long-term operational monitoring experiments because of coarse image resolutions (resulting in low accuracies for detecting small displacement) or long-time intervals (mostly related to meteorological conditions, and the absence of acquisitions in the nights) between subsequent observations. For many geomorphological processes, and more specifically for landslides, acquiring ground-based image time series at high temporal frequency (typically 1 h) is a relevant source of information. Terrestrial optical photographs (mono-view, pair-views or multiple views) allow identifying the movement pattern (in either 2D or 3D (Travelletti et al., 2012)), the reconstruction of the surface topography (Kaufmann, 2012;Stumpf et al., 2014) or the identification of discontinuities in the unstable mass (Sturzenegger and Stead, 2009) reveal the spatial and temporal patterns of acceleration and mass accumulation, and thereby help to refine conceptual and numerical models. With only one camera, the approach is not fully 3D and still depends on the availability of frequently updated topographic surface models but could be extended with stereo-and multi-view camera setups. Travelletti et al. (2012) analyzed time series of incremental displacements by cumulating the displacement over the complete epoch at the pixel scale. Such approach is valid for small time series but is not very efficient for detecting changes in longer time series (e.g., 30 time steps) for which fully automated processing and adapted (to different study cases) correlation approaches are needed.

Ground-based camera-derived terrain motion techniques are also applied for river sediment discharge analysis (Bradley et al., 2002;Stumpf et al., 2016) where long-time series of video frames are processed using Particle Image Velocimetry (PIV) approach (Adrian, 1991). In these cases, the fluid flow is free surface and is therefore part of a non-Newtonian approach. This approach cannot be applied to complex movements characterized by several rheologies. The techniques used in the laboratory are the closest to our cases. However, the resulting conditions and assumptions (near fields, controlled illuminations, known magnitude and direction of expected deformation) cannot be applied to a geomorphological use case.

Further, most geomorphological examples currently published were tailored for very specific use case.

Currently no generic pipeline adaptable to any site configuration (slope dimensions, slope deformation rates, distance between the slope and the camera) and any camera type (sensor size, focal length) exists (Pham et al., 2014;Marsy et al., 2018).

Our works target the presentation of the "TSM-Tracing Surface Motion" processing toolbox. TSM is a generic and fully-automated toolbox, built on previouswork (Gance et al., 2014;Travelletti et al., 2012;Stumpf et al., 2016Stumpf et al., , 2018)), for the processing of ground-based optical image time series integrating several modules from the pre-processing of the raw images to the calculation of displacement fields, and the identification of slope units of similar kinematics. The method is applied and evaluated on two landslide use cases: a very slow movement (Chambon landslide; Romanche valley, France) and a moderate movement (Pas de l'Ours landslide; Guil valley, France) based on rates published by Cruden and Varnes (1996). The sensitivity of the method to the processing parameters and to the strategy used to combine the images ('leap frog', Common Master approach) are assessed. The potential of the system to detect the motion pattern of the two landslides compared to reference measurements (laser scanning point clouds and on-site GNSS acquisitions) is discussed.

First, we present the methodology applied to process large image time series (30 images) (Section 2) from the estimation of the camera parameters to the extraction of relevant information in time. Second, the possible strategies used for processing large datasets are discussed (Section 3). Third, the datasets acquired on the use cases are presented (Section 4). Fourth, the sensitivity of the method is assessed (Section 5). Finally, the results obtained on large time series are discussed (Section 6 and Section 7).

Methodology

The methodology to calculate deformation fields (both in the image plane geometry and in the ground geometry) is presented in

Module 1: Orientation estimation

This module aims at calculating the external orientation of the camera installed in the field from a series of Ground Control Points (GCP). The camera calibration is essential to get metric information for threedimensional applications. To obtain precise measures, it's necessary to take into account the camera model used. The most common are the spherical model (Feng et al., 2012;Aghayari et al., 2017) and the pinhole model (Figure 2.2). As the camera lens of a spherical model introduced high distortion, it's more easy to deal with the pinhole model which is commonly used. Here, we propose two methods for the orientation estimation in the framework of a pinhole model.

Field calibration without a priori

As explained previously, internal parameters are not introduced currently in the processing. External parameters are deduced from the collinearity equation's, Equation (2.1):

               u = u 0 -f u m 11 (X -X cam ) + m 12 (Y -Y cam ) + m 13 (Z -Z cam ) m 31 (X -X cam ) + m 32 (Y -Y cam ) + m 33 (Z -Z cam ) v = v 0 -f v m 21 (X -X cam ) + m 22 (Y -Y cam ) + m 23 (Z -Z cam ) m 31 (X -X cam ) + m 32 (Y -Y cam ) + m 33 (Z -Z cam ) , ( 2
RMSE xy = 1 n × ∑ ∆x 2 + ∆y 2 ,
(2.2)

∆x i = x pro j,i -x picking,i ∆y i = y pro j,i -y picking,i , (2.3)
where (x pro j,i , y pro j,i ) are the coordinates of the ith re-projected GCP in the image-based coordinate system, (x picking,i , y picking,i ) are the coordinates of the ith picked GCP in the image-based coordinate, and n is the number of GCPs.

Because the estimated parameters are dependent on each other, there is no information about the external parameters and the measure of the camera position is biased.

Field calibration with a priori

The calibration with a priori is used to estimate each parameter independently of each other. As introduced previously, the camera orientation is modelled with internal parameters (to simulate the camera optics)

and external parameters (to determine the transformation between the image and the world coordinate systems). The relation between both coordinate systems can be expressed by Equation (2.4). Let's consider P in the image plane and Q the associated point in 3 D as drawn in Figure 2.2:

sP = A R|t Q (2.4)
where A is the intrinsic matrix commonly called the matrix calibration, R|t is the extrinsic matrix composed of a rotational matrix R and a translation matrix t, and s is an arbitrary scaling factor which is independent of the camera model. The camera parameters are described with more details in the annex A. We use the Python toolbox OpenCV© (Zhang, 2000;Bradski and Kaehler, 2008;Wang et al., 2010) to run the calibration by integrating a priori on the intrinsic and extrinsic parameters. OpenCV©is an Open source Computer Vision library written in C-language which is provides several image processing algorithms and computer vision. Calibration algorithm is based on Zhang's calibration method which assumed a two-dimensional template as study object. Camera parameters are estimated by the DLT method as previously done. Using these parameters, the distortion coefficients are then estimated by minimizing the distances between the image points and the reprojection points. To constrain the algorithm to determine the parameters that best fit the position of the camera relative to its true position, the camera matrix, the distortion coefficients and the rotation and translation matrices are inverted by applying a differential evolution algorithm (Storn and Price, 1997). selection is, by default, applied for the full image size, but options for masking irrelevant areas (sky, vegetation, water bodies) can be applied. Masks are based on the principle that each area gets is own radiometric signature and can be extracted from the three components.

Module 2: Image selection

Module 3: Correction of sensor movement

The module aims at correcting the image time series from camera movements. For many on-site applications, the cameras are affected by significant movements induced by temperature, wind or terrain motion. This induces non-negligible systematic offsets between two images of the time series. To correct sensor movement, techniques based on characteristic points are used, such as the Harris detection method (Gance et al., 2014), the feature detection method (Zitova and Flusser, 2003) and the SIFT method (Lowe, 2004). However, the number of points taken into account in the transformation calculation and the distribution of these points in the image vary from one image pair to another. This has a direct impact on the accuracy of co-registration of the time series. To overcome these limitations, the module implements a statistical co-registration initially tailored for satellite imagery (Stumpf et al., 2018). The approach consists in correlating two images in order to construct disparity grids ∆X DISP and ∆Y DISP and a Correlation Coefficient (C) grid (Figure 2.3). The systematic offsets (∆X DISP , ∆Y DISP ) are modelled considering two planes ∆X model and ∆Y model . As the movement of the camera results mostly from a translation and a rotation, the two planes are modelled with an affine transformation, Equation (2.5):

∆x model,i = a + b * x DISP,i + c * y DISP,i ∆y model,i = a + b * x DISP,i + c * y DISP,i , (2.5)
where ∆x model,i and ∆y model,i represent the modelled offsets of the ith pixel and x DISP,i , y DISP,i are the raw disparities of the ith pixel. The parameters a, b, c, a , b , c, are estimated to reduce residuals between both offsets using an Iteratively Re-weighted Least Square (IRLS) (Stumpf et al., 2018). The corrected grids are obtained by subtracting the modelled planes to the disparity grids with Equation (2.6):

∆X CORRECT ED = ∆X DISP -∆X MODEL ∆Y CORRECT ED = ∆Y DISP -∆Y MODEL , (2.6)
where ∆X CORRECT ED and ∆Y CORRECT ED are the corrected disparity maps. To disregard the areas in the images where real terrain motion needs to be quantified, the systematic offsets are calculated on sub-portions of the images considered as stable. A mask, defined by the user, is applied for selecting only the relevant sub portion of the images.

Module 4: Quantification of terrain motion using matching techniques

The module aims at quantifying changes between two acquisitions through the calculation of disparity vectors using cross-correlation (Hild and Roux, 2012). The disparity vectors are first expressed in pixels as they are determined in the image plane. Three grids are computed with the open source photogrammetric library MicMac corresponding to the relative displacement in line, the relative displacement in row and the correlation coefficient. The image correlation algorithm uses a multi-scale and hierarchical approach based on Normalized Cross-Correlation (NCC) and a non-linear cost function. A regularization parameter is introduced in order to avoid matching ambiguities and to have the ability to process small matching windows (Pierrot-Deseilligny and Paparoditis, 2006;Rosu et al., 2015;Rupnik et al., 2017). The cost function to minimize is defined by Equation (2.7):

E α = ∑ (1 -Corr) + α × F, (2.7)
where (1 -Corr) is the similarity measure with Corr being the NCC score; α being a weighting parameter and F being a positive function representing the regularization term. The matching is computed only in the spatial domain (Rosu et al., 2015). The parameters used in the module are presented in section 5.

Module 5: Filtering and detection of meaningful motion patterns

The module aims at filtering ambient noise induced by changes in illumination (affecting the reflectance of the pixel) or changes in soil surface state (humidity, large displacement between two images) creating loss of coherence. The objective is to filter the displacement fields of Module 4 using the criteria proposed by Stumpf et al. (2014) in order to highlight meaningful motion patterns. Filtering is applied jointly on the correlation coefficient C (see section 5.1.2) and on the direction and norm of the disparity vectors.

Indeed, disparity vectors which are believed inconsistent with a gravitational movement (e.g., for instance disparity vector pointing towards the top of the terrain) are removed. Agreements between the direction of the disparity vectors and the aspect of the terrain slope are determined. The values of terrain direction angles are adapted for each use case.

Module 6: Geometrical correction of the displacement fields

This module aims at converting the relative displacement fields (e.g., in pixels) in absolute displacement fields (e.g., in a projected geographic coordinate system). The projection from the relative coordinate system in the image geometry (∆u, ∆v) to the absolute coordinate system in terrain geometry (∆X, ∆Y, ∆Z) 3 Combination strategies for processing large image datasets Several combinations of images are possible to process large time series (Marsy et al., 2018;Pham, 2015).

is
TSM implements three approaches (Figure 2. with the incremental parameter n being variable. In the case of n being fixed (n = 1), this method is known in the literature as "leap frog";

• (3) the Redundant Variable Sequential Correlation method (RVSC), which consists in correlating p images at I with p other images at I + n, with the incremental parameter n not being fixed. In total, p 2 correlations are performed for the forward matching (Figure 2.4c). The average displacements are then calculated. This approach allows obtaining redundant information (Stumpf et al., 2017) to increase the Signal-to-Noise ratio and improving the detection of meaningful surface motion.

This technique builds on the Multi-Pairwise Correlation method (MPIC) proposed by (Stumpf et al., 2018) for processing time series of satellite imagery. Correlations can be performed for each interval [I, I + n] in forward and backward mode (Figure 2.4d). The ISO sensitivity and the aperture opening are fixed at 100 and 5.6 respectively, and are kept constant for all acquisitions.

Sensitivity analysis

The sensitivity of TSM to several processing parameters is evaluated. Among the internal parameters, the influences of the image modalities (section 5. 4. the "texture-derived image T*" corresponds to a Sobel filter applied to the grayscale converted images using a 3 × 3 convolution matrix (Chambon, 2003;Travelletti et al., 2012). The Sobel filter is used to highlight natural edges and to improve image quality which can be degraded by weather conditions or illumination (Berthier et al., 2005).

These four modalities can be applied either on the initial RGB image or on the enhanced image. The For each pair of image modality, block matching correlation is carried out with fixed parameters. The co-registration residuals expressed by the RMSExy for all pixels are calculated (Equation (2.2)). As outliers Large variations are observed: the higher C min is, the smaller the percentage of pixels is. This explains why most published studies used a value close to 0.5 as a compromise on the quality and the quantity of the information. For landslide analysis, it is also necessary to take into account the change in intensity of pixel radiometry that can decrease the correlation coefficient. Finally, five tests were carried out in both cases to document the sensitivity of the regularization parameter, in the range from 0.2 to 1.0 with a step of 0.2 (Figure 2.6d). The amount of noise and outliers is decreasing when the parameter increases because of the spatial smoothing of the displacement field (Rosu et al., 2015) that translates in an increase of the amount of pixels with information. The parameter has small influence in our case because the C values are already high.

Sensitivity to the Size of the Stable Area Used for the Co-registration

The size of the stable area used for evaluating the quality of the co-registration is assessed (Figure 2.2).

The size of the area is expressed by the number of pixels in the stable area according to the total number of pixels in the image. As described in section 3, a mask is generated before the co-registration step in order to calculate the statistical fitting plane. Figure 2.7 presents the distribution of residuals in the (u, v) directions for the stable area for a given image pair. The size of the stable area introduced in the statistical modelling varies in the range 36 % (full stable area) to 6.5 % for the Chambon landslide and from 9 % (full stable area) to 1 % for the Pas de l'Ours landslide. In both cases, the larger the stable area is, the more accurate is the quality of the co-registration. However, above 4.5 % of stable area for the Pas de l'Ours landslide and 9.7 % for the Chambon landslide, the distributions of residuals are similar. As it is shown in Figure 2.7a, the simultaneous use of stable areas in the near-and far-fields allows us to increase the standard deviation and the average shift.

Sensitivity to the Image Combination Strategy

Several image combination strategies (CMC, VSC, and RVSC) are evaluated (Figure 2.4) through three comparisons: (1) CMC vs. VSC;

(2) VSC vs. RVSC; and (3) forward vs. backward correlation. Then, the stability of the RMSE xy parameter (Equation (2.2)) along time, is used as criteria to compare the three combination strategies along the time series.

Figure 2.8 shows the cumulative displacements of four selected points, located either in the stable or in the unstable parts, calculated with the VSC method and with the CMC method using two masters (master 1, master 2) and one slave image. Because the CMC method can be strongly dependent of the choice of the master image, two masters are selected (no meteorological alterations such as rain, mist or shadow). For the Chambon use case (Figure 2.8a) characterized by slow motion (displacement rates less than 1 pix.day -1 ), both methods give similar results since the accuracy is about a half-pixel. For the Pas de l'Ours use case (Figure 2.8b) characterized by moderate motion (displacement rates higher than 10 pix.day -1 ), both methods give the same results except for Point 2. The displacement calculated between the master 1 and the slave image differs from those calculated by the VSC method and those calculated with the master 2 (> 10 pixels). This can be explained by the quality of the measure given by the correlation coefficient. Actually, C < 0.3 for the displacement measured between the master 1 and the slave image and C > 0.9 from those measured between the master 2 and the slave. For both use cases, the VSC method performs better for displacements > 1 pixel among two dates; applied to use cases with small displacements (< 1 pixel), the VSC method gives noisy results as small displacements are cumulated. The CMC method gives satisfactory results in both cases (low or high displacements). However, it can lead to artefacts when the master and slave images are heterogeneous in terms of radiometry.

The VSC and the RVSC methods are compared for the full time series. The standard deviations of all residuals located in the stable area are calculated and averaged. The standard deviations (σ u , σ v ) in both direction (u, v) for the Chambon and the Pas de l'Ours landslides are respectively (0.22, 0.21) pixels and (0.97, 0.65) pixels for the VSC method, and respectively (0.16, 0.16) pixels and (0.96, 0.67) pixels for the RVSC method. Even if some residuals are still observed due to some defectiveness of the co-registration, the RVSC method tends to improve the detection accuracy and reduces false detections.

The influence of forward and backward calculation (Figure 2.4) is evaluated for both use cases on a pair of images. The standard deviations of all residuals are calculated. The standard deviations (σ u , σ v ) in both direction (u, v) are respectively, for the Chambon and the Pas de l'Ours landslides, (0.14, 0.19) and (0.22, 0.17) pixels in the forward mode and (0.14, 0.19) and (0.21, 0.17) pixels in the backward mode in the forward-backward mode, they are respectively (0.13, 0.18) pixels and (0.20, 0.15) pixels.

Compared to the RVSC method, the forward-backward approach is less powerful to improve the detection accuracy significantly.

The results of the RMSExy stability are shown in Figure 2.9. The RMSExy parameter (Equation (2.2)) was calculated for the three methods and two master images were selected for having two time series in a CMC mode. The co-registration errors vary according to the correlation method used; in particular the RMSExy relative to the CMC method seems to depend on the defined master. Indeed the temporal stability can be disrupted by the surface-changes which can lead to significant variations on the stability (Figure 2.9a). As the time series stability of the Pas de l'Ours use case (Figure 2.9b) is not affected by the master choice, we assume that the stability of the calculation depends on both the method used and the study sites.

Afterward, we used a combination of the parameters (Table 2.1) which seemed the most appropriate to the use cases.

Role of the external parameters

To estimate the influence of external parameters on TSM results, we conducted a statistical analysis based on the OAT (One-at-a-time) method. This method consists of varying only one parameter at a time, the others being fixed. The parameter variation is a Gaussian noise applied randomly to the data and calculated according to the measurement errors associated with the data. The analysis was carried out on a data set other than that of the Chambon and the Pas de l'Ours landslides for timing reasons. For more clarity, the results are in the Appendix. Nevertheless, it results that the errors related to the noise relative to the in-situ data are of the order of several centimeters, those relative to the DEM are of the order of centimeters and those relative to the picking are of the order of meters for the noise applied in the direction perpendicular to the movement and less than one millimeter for the direction of the movement.

Results and Discussion

The quality of the camera-derived displacement fields is evaluated in comparison to reference geodetic observations and the pros and cons of the TSM toolbox are discussed.

Displacement Fields at the Chambon Landslide

A time series of images spanning from 8 to 27 March 2017 is processed to illustrate the TSM application and characterize the motion of the Chambon landslide. Figure 2.10 shows an example of a displacement fields (in cm.day -1 ) for the image pair of 9-12 March 2017. The image pair was processed with the VSC method in the forward direction taking the first image as master. The most active part of the landslide is clearly depicted with displacements rates > to 3 cm.day -1 ; the stable parts are also highlighted with very small displacement rates lower than 0.6 cm.day -1 indicating the sensibility and robustness of the method.

Two surface velocity profiles are presented in increase. This pendulum motion is in agreement with a rotational movement which is confirmed by topographic surveys (Mathy and Lorier, 2013;Desrues et al., 2019). The profile P2 identifies two distinct compartments, with a slide-type behavior (from 2 to 50 m) and a subsidence type behavior (from 60 to 92 m).

To validate the displacements derived from the image time series, we compared the TSM outputs to displacements calculated from two high-density point clouds acquired on 9 March and 10 April 2017 with a terrestrial laser scanner RIEGL VZ-2000. Point clouds were processed using the approach developed As there is no detectable movement between the 8 th of April and the 10 th of April according to the Figure 2.11, we can assume that the two results are similar. To highlight the movement and reduce the ambient noise, images were correlated according to the RVSC method in the forward direction. 79 points were chosen randomly (with 43 points inside the active area of the landslide, and 36 points outside the landslide area) and displacements were calculated over a window of 16 × 16 pixels (which represent a window size comprised between 56 × 56 cm and 96 × 96 cm). The mean and standard deviations between the two displacement grids are respectively -0.013 m and 0.043 m.

Displacement Fields at the Pas de l'Ours Landslide

A time series of images from 6 to 22 June 2018 is processed to illustrate the TSM application and characterize the motion of the Pas de l'Ours landslide. Figure 2.12 shows an example of displacement fields (in m.day -1 ) for the image pair 16-18 June 2018. Images were processed with the VSC method in the forward direction. Two sections in the landslide are clearly depicted: the lower part which presents low displacement rates < 0.10 m.day -1 and the upper part which presents high displacements > 0.10 m.day -1 to > 2.0 m.day -1 . The accuracy of the measurements is evaluated on a stable slope located in the first quarter of the image; the standard deviations of the displacement rates are 0.07 m.day -1 .

Advantages and Limitations of TSM Pipeline

The TSM methodology for the massive processing of time series has proved to be robust and applicable to two sites with different field setup and terrain motion. The coupling of several correlation techniques has proven to be an asset to detect small change in velocity but several factors limit the capacity of the technique to detect terrain motion above a cm.day -1 . For instance, in most cases the choice of the camera position depends on the topography of the site (e.g., viewing geometry) and on the access path. This constrains the camera-to-object distance, the viewing geometry (e.g., camera line-of-sight vs.

direction of the terrain motion) and the timing of image acquisition (depending on the incidence of the sun).

The hypothesis of negligible radial-distortion is not true if the camera is equipped with very wide view lenses (focal length < 24 mm) for which lens distortion should be systematically corrected. However, the analysis of very wide view images like the ones acquired with webcam should provide more stable information over time due to the system's robustness to climate, its lower energy consumption compared to SLR cameras, and the possibility of integrating larger stable areas in the image view.

For the operational use of TSM in monitoring conditions, several factors need to be taken into account to leverage the detectability of small movements and increase the accuracy:

• The implementation of GNSS field targets for the estimation of the orientation parameters of the camera. Most of the time, some areas of the images are not covered by GCPs because of the field access. This has direct consequences on the calibration accuracy since it is necessary to have the Chapter 2

Image time series analysis for a monoscopic model most homogeneous distribution of GCPs in the image.

• The optimization of the image matching calculation; using a Debian 64 bits computer with a 64 Go RAM and a 3.50 GHz processor, and for images acquired with 18Mpix cameras (corresponding to 5184 × 3456 pixels), one correlation between an image pair takes about 20 min (without any pre-processing such as greyscale reduction) for an 8 × 8 pixels search window size and can exceed 60 min if the pixel search window increases (> 12 × 12 pix). High-Performance Calculation (HPC) is compulsory for processing massive datasets and calculating all the possible combination of image pairs.

• The optimization of the change detection calculation. The filtering of the stacks of correlation grids (e.g., RVSC approach, combination of forward and backward correlation) is also time-consuming and HPC is also compulsory for operational uses.

The process presented here requires some a-priori information (dimensions, surface velocity estimation)

to adapt the camera field setup and the parameterization of TSM.

Conclusions

This work targeted the development, testing and operational use of the TSM toolbox as a generic pipeline for the analysis of large time series of terrestrial optical images for the detection and quantification of terrain motion. The TSM toolbox is composed of several modules to process images acquired with a fixed monoscopic camera. It includes six modules that can be used independently of each other or sequentially:

the sensor orientation module, the image selection module, the sensor movement correction module, the quantification of ground motion module, the detection module and the geometrical correction module.

The sensitivity analysis carried out on image dataset acquired on the landslide use cases revealed the low influence of image modalities (G*, B*, A*, T*) on the matching results. The TSM method is sensitive to the correlation parameters such as the search window and the correlation threshold which need to be tuned for each case studies. For landslide applications, the parameterization of TSM is constrained by the amplitude of the displacements and the type of deformation which will affect the quality of the results.

Several calculation strategies are possible to combine the image time series; the implemented VSC and RVSC combination methods performed well for the two use cases. The stability and accuracy of the CMC combination method is highly dependent of the choice of the master image.

Comparison with displacements calculated from high density point clouds for the Chambon landslide and in-situ measurements for the Pas de l'Ours landslide show a standard deviation of respectively 4.3 cm and 16 cm. The displacement field is further consistent with previous studies and technical reports (Mathy and Lorier, 2013;Laurent et al., 2016), (Point et al. (submitted)).

The two use cases demonstrated the robustness and large applicability of TSM since it was applied on both slow and moderate motion landslides with several configurations (slope dimensions, slope deformation rates, distance camera-slope).

Chapter 3

A stereoscopic model for landslide analysis: Application to the Montgombert landslide (Savoie, French Alps)

This chapter is devoted to the development of a method to infer landslide 3D displacement fields and further landslide geometry from time series of optical images acquired from two remote locations. The method is developed and tested for the Montgombert landslide (Savoie, France) use case and the quality of the information obtained by image analysis is assessed with in-situ measurements. The sliding properties are approached thanks to the 3D surface displacements which give both access to tangential and normal displacements.

The chapter is divided into two main sections: (1) the first section presents the methodology developed to construct dense 3D displacements from a stereoscopic camera system and (2)

the second one presents a technique to infer a field of depth of mass movement by using the 3D displacement field reconstructed, together with the digital elevation model (DEM).

In this chapter a stereoscopic model refers to a system of two fixed cameras for the acquisition of the time series.

This chapter is based on an article in preparation: Desrues M, Toussaint R, Malet JP, Brenguier O (2021) Landslide thickness estimation inferred from terrestrial stereoscopic optical images (Savoie, French Alps).

1 Foreword

Context setting

The analysis of displacements and velocity fields is a useful tool for monitoring gravitational instabilities.

It allows to define the state of stability of the slope, its evolution over time, its spatial coverage and to identify if the movement is progressing laterally or retrogressively. Other types of information can also be extracted from this analysis such as the tangential and normal deformation calculated locally, together with a DEM, or the strain field that highlights areas of compression and extension (Travelletti et al., 2014). The kinematic analysis, that can be qualified as classical approach, becomes therefore a deformation analysis that enables to assess the mechanical properties of the moving mass such as its thickness, providing a better understanding of the movement (Baum et al., 1998). To do so, it is consequently important to ensure that the computation of the 3D displacements provides both information on normal and tangential displacements.

Several methods give access to complete 3D movement information including the monoscopic approach thanks to the integration of a DEM in the image plane. However, this method suggests that the motion is mostly tangential to the surface. Complete information can also be obtained by integrating two DEMs in the camera image plane at two different dates, Displacements are then calculated by projecting the pixel points onto the respective DEM. This technique is limited to specific cases in which normal displacements are significant to be detected and pixel displacements trackable to avoid de-correlation in the image correlation process. This technique requires both a time series of images and a time series of DEMs. Thus, to take into account the full information with the few possible data, and to detect tangential and normal displacements, we developed a methodology to construct the 3D displacements directly from the correlation results from a pairwise of fixed camera and tools to perform deformation analysis.

Most of the pre-existing methods used for estimating the sliding surface depth are simple techniques, based on in-situ methods that provide sparse data without any surface information. This is the case of geotechnical techniques (e.g., boreholes) or some geophysical prospecting that may be time consuming since the fieldwork is realized on hazardeous areas. Other methods combine in-situ data with surface information (DEM and field observations) and use, for instance, some polynomial fits to reconstruct the slip surface in between the localised information (Tang et al., 2020). In recent years, the thickness of the moving mass has been calculated from the surface displacements derived either from surface topography (Farinotti et al., 2009;Aryal et al., 2015) or from satellite synthetic-aperture interferometry (InSAR) (Morlighem et al., 2011;Nikolaeva et al., 2014) of from orthorectified optical images (Delacourt et al., 2007;Booth et al., 2013). The thickness of the moving mass can be estimated by several methods. Bishop (1999) introduced the application of the 2D balance cross section (Woodward et al., 1989) for the estimation of the depth of translational landslide. This geometrical method suggests that deformation is purely in-compressible. Booth et al. (2013) introduced also a new method to estimate the thickness of a moving mass by inverting 3D displacements calculated from the correlation of aerial optical images.

This approach suggests in a similar way than the previous method, that deformation is in-compressible.

Contrary to previous methods, the estimated thickness is calibrated with in-situ measurements.

Thickness can also be estimated assuming that deformation evolves in a homogeneous elastic halfspace. This method refers to the elastic dislocation models introduced by Okada (1985) that suggests that slip is more important at the sliding surface than at the ground surface. It infers depth of sliding surface but also fault orientation. Nikolaeva et al. (2014), for instance, applied the elastic dislocation model to find a possible dislocation plane by inverting displacements derived from InSAR measurements. Aryal et al.

(2015) tested and compared the 2D balance cross section approach with the elastic dislocation method on displacements derived from terrestrial laser scanning. They highlighted some limitations such as the elastic dislocation plane is only applied at the initial phase of the landslide where there is no significant inelastic displacement and the balance cross section is restrained to translational landslides. Landslides can also be modelled by considering various velocity profile such as, in particular, a viscoplastic model [START_REF] Savage | A model for creeping flow in landslides[END_REF].

To estimate the thickness of the moving mass from 3D displacements derived from stereoscopic optical images, we propose a methodology based on the law of conservation (i.e., displacement in-compressible) by invoking the rheology of the material involved (Booth et al., 2013). In order to take into account, in this model, a more complex slip geometry, we introduced a disbonding parameter that marks the presence of a dislocation area at the top limit of the moving mass which implies a non-zero velocity at the sliding surface.

In this chapter, we propose to test the different methodologies (i.e., reconstruction of the 3D displacements with a stereoscopic approach and estimation of the depth of the sliding surface) by applying them to the Montgombert use case presented in the next section.

The Montgombert use case

The Mongombert landslide is located in the Arly valley in Savoie (France) and is mainly composed of mica schist scree. The landslide has been monitored since 2016 since it undergoes frequent accelerations.

The most active parts are located at the base of the landslide above a stable 50 m high cliff, crossed by three propagation corridors. Field observations indicate that the movement retrogresses uphill, that the lateral boundaries are susceptible to change, and that the movement is composed of several overlapping compartments (Mathy and Lorier, 2018). As shown in Figure 3.1, these observations are consistent in time. Because of the dense vegetation and the difficult access to the lower part of the movement, the observations were concentrated on the upper part of the landslide. In March-April 2020, heavy rains destabilized the lower part of the movement (Figure 3.2). Some displacements in the order of a meter were measured by total station. However, as shown in Figure 3.2, highest displacements could not be measured, because of the loss of some topographic targets due to their tilting or tearing off during the movement. In total, 32 topographic targets were measured in a hourly basis before the event and only 13 could be measured after. Chapter 3

A stereoscopic model for landslide analysis: Application to the Montgombert landslide

In addition to topographic measurements, two inclinometers were set up in the field as well as two piezometers (Figure 3.1). These two geotechnical equipments were installed uphill of the visible area covered by the time lapse photograph. As mentioned by Mathy and Lorier (2018), the inclinometer at the 700 m elevation (Figure 3.21) is sheared at three different depths: at 8 m, 15 m and 30 m from the natural terrain. As for the second inclinometer at the 794 m elevation, it is sheared at 13.5 m. This suggests well the complexity of the movement which is declined in several lobe slides which evolve and retrogress upstream (see section 4.2.2).

We first present the stereoscopic model used to generate the 3D displacement fields. A comparison with the in-situ measurements and the measurements from a monoscopic model is then led to assess the relevance of such a reconstruction. Finally, some thickness and volume estimates are discussed.

2 Stereoscopic model

Two-camera model

Stereophotogrammetry is the ability to build a 3D object from its observation from several points of view.

In the case of a single camera, the reconstruction is dependent on the used camera model (pinhole, thin lens, etc). A pinhole model relates the coordinates in the world system to those in the image system by a simple relationship as shown in Appendix A. In the case of a two-camera model, this relation is based on an epipolar geometry (Hartley and Zisserman, 2004), which is dependent on the intrinsic and extrinsic parameters of the cameras, but does not depend on the structure or content of the object. The epipolar geometry is constructed from an epipolar plane defined by the baseline (i.e., the straight line connecting the projections of the camera centers) and a given point M (see Figure 3.3). The geometry assumes that the point Q, homologous to P, is located on the epipolar line in the image on the right. All epipolar lines intersect at the epipoles.

Let P be a point belonging to the left image and Q to the right image. For each homologous pair (P, Q) there is a point M of coordinates (x,y,z) in the world coordinate system. The relation that governs this geometry is the one linking each pair of homologous points (P,Q) by the fundamental matrix F:

P T FQ = 0 (3.1)
where

P =    x p y p 1    Q =    x q y q 1    F =    F11 F12 F13 F21 F22 F23 F31 F32 F33   
P and Q are defined in the camera coordinate system. If this equation is not satisfied for a pair of points, then they are not homologous. The fundamental matrix can be determined from a set of at least eight known homologous points by using the eight-points algorithm proposed by Hartley and Zisserman (2004).The fundamental matrix is essentially used to reconstruct the 3D model in the case of a pair of uncalibrated cameras. Indeed, by calculating F, we refrain from calculating the intrinsic parameters of the cameras.

Camera positioning

The position of the cameras is, first, defined by the difference in the angle of view of the cameras, which must not be too large, otherwise there will be very little overlap in their capture images and, therefore, very few similar points. This difference in the angle of view is constrained by the distance between the two cameras and the distance between the cameras and the object. These two parameters are combined together to form a ratio which, when close to 1, means little overlap but distinct angles of view and, when low, means high overlap but similar angles of view. This ratio is called the b/H ratio and is shown below.

The position of the cameras is also defined by the angles that characterize the external orientation of the cameras in the world coordinate system.

b/H ratio

The b/H ratio is the ratio of the distance between the two cameras, commonly called baseline (b), to the distance between the baseline and the object called height (H) (see Figure 3.4 a)). This ratio is important to consider in the case of aerial or satellite stereoscopy. Several studies recommend a base-height ratio in the range of 0.5 to 0.9 in order to obtain the most accurate DEM possible (Hasegawa et al., 2000). Indeed, this ratio directly influences the image matching, since it constrains the difference in the angle of view of the cameras (Hasegawa et al., 2000;Aguilar et al., 2013). As far as terrestrial stereoscopy is concerned, few studies recommend that cameras be installed in a precise ratio. Gance et al. (2014) experimented with a b/H ratio between 1.6 and 2.1 and it would seem that these values are too high for a terrestrial stereoscopic system. In addition, a baseline that is too large may bias the matching (Moser, 2012; Gance The cameras are fixed to the bedrock with the help of poles (see Chapter 1) in order to minimize their displacement over time. The difference of respective distance between the cameras and the object implies that the ground pixel size differs from one camera to the other. Its size ranges from 5 to 6.5 cm for the left camera, and from 4 to 5.5 cm for the right one. Therefore the resolution of the images, as well as that resulting from the correlation measurements, are different.

External camera orientations

Each camera is oriented in space using three angles commonly named Euler angles (ω, φ , κ) (see Appendix A); in the airborne photogrammetry domain, one might rather use pitch, roll and yaw. These angles are zero when the camera is horizontal and pointing North. In the world coordinate system, defined by the geographic coordinates (X, Y, Z), ω is the rotation around the X-axis; φ is the rotation around the Y-axis and κ is the rotation around the Z-axis. The sign of these angles are defined according to the right-hand rule. Euler angles are involved in the rotation matrix used to link the image coordinate system to the coordinates of the world system.

The cameras were installed without constraining Euler's angles. Indeed the installation of the cameras was a complex, since the installation was mainly dependent on the environment (restricted access requiring a rope, road, presence of protective grid). It was therefore necessary to determine the orientation of the cameras, a posteriori, by using the camera calibration introduced in Chapter 2 and showed in section 2.3.2.

Image pairing

The matching process consists in finding all pairs of homologous points (P,Q) in stereoscopic images.

These points are specific features that are unique and easy to track such as a corner, a line intersection or an edge (Harris et al., 1988;Zitova and Flusser, 2003). They are characterized by distinctive descriptors, which correspond to regions located on the features. The neighborhood allows the algorithms to situate the feature in the right place on both images. As mentioned in Chapiter 1, several algorithms can be used to detect these features. Among the most well-known are the classical feature-based algorithms such as the Harris Corner Detector or the SIFT detector (Le et al., 2011;Hassaballah et al., 2016) and the image-based matching such as the image correlation technique.

Classical methods

The image-based matching is a robust method that provides results with high spatial overlap. However, the results of the image matching proved to be unsatisfactory as very few pixels could be matched. The difference in viewing angle is probably too high in our case to ensure a dense matching image. Therefore we turned to the application of a feature-based algorithm.

We used the ORB (Oriented FAST and Rotated BRIEF) algorithm to detect the features (Rublee et al., 2011). Faster than SIFT detector, this algorithm is also rotation invariant and resistant to noise. The number of homologous points found on the series of stereo pairs varies from 150 to 185 points. As shown in Figure 3.5, the detection of homologous points does not lead to a unique result (several points matched in one image for a given point in the other one). This can be explained by the presence of similar objects (vegetation, tree trunks) and heterogeneous objects such as scree or soil areas distributed throughout the image. In fact, this method is often used in the case of object related to streets infrastructure since the features are very well defined in these cases. Here, the characteristics of the features are too volatile. The few homologous points compared to the number of pixels in the image, as well as their spatial distribution in the image with respect to the area of interest, makes this technique unsatisfactory for our study. Indeed, although only few points are necessary, when one wishes to register an image compared to another, it is not the case when these points become features to be tracked in time. In our case, this technique does not allow us to obtain the densest possible network to be used afterwards.

Back-Projection method

In order to obtain dense pairs of homologous points spatially distributed homogeneously in the image, we decided to use a Digital Elevation Model (DEM), which is involved in the calculation of pixel to meter transformation (see Chapter 2). The workflow used is organised as follows:

1. Camera calibration: Nine artificial targets, with absolute positions measured by GNSS!, were used for the camera calibration. The orientation of the cameras was first roughly estimated (as an initial parameter for a later inversion process). We do similar estimations for the calibration matrix which is composed of the effective focal length ( f x , f y ), and the coordinates of the principal point, which are commonly estimated to the first order as those of the center of the image (the skew parameter, which is usually equal to zero, is not expressed in OpenCV© (Zhang, 2000) -(see Chapter 2). The distortion coefficients have been set to zero. The calibration algorithm was developed in the case of planar objects. In our case, regarding the camera-object distances as well as the geometrical shape of the object, we can assume that we are in a similar case. Indeed, if we consider a plane formed by the three targets A, B and C, the distances of the other targets to this plane deviate by a maximum of 2 % (Figure 3.6).

The inversion is performed by applying a differential evolution algorithm which is a stochastic algorithm (Storn and Price, 1997). This implies that the results from the calibration are not unique.

In order to determine the most appropriate combination of parameters, several inversions have been performed and the one that minimizes the error on the camera position and gives Euler angles such that the camera's reference system is orthonormalized, is selected. The advantage of this technique is to be able to match the maximum number of pixels over a large area.

The image pairing by searching for the pixel according to the coordinates X, Y, Z is carried out by the nearest neighbor algorithm. The calculation time is therefore directly related to the number of pixels considered, since the calculation is performed pixel by pixel. To get rid of the long computing time (> 4h) per image pair, we processed the data at low resolution considering a step of 25 pixels. For an average physical pixel size of about 0.05 m, this means calculating the 3D displacement every 1.25 m if the displacements are distributed on a regular grid.

Inversion of the 3D displacement fields from the 2D displacement fields in pixels

Let us consider that (u 1 ,v 1 ) are the pixel displacements of a point P in the left image and (u 2 ,v 2 ), those of the homologous point, Q, in the right image. Thus, the pixel displacements are related to the metric displacements (dx, dy, dz) by the matrix A, Equation (3.2):

      u 1 v 1 u 2 v 2       = A •    dx dy dz    (3.2) Let us call u = (u 1 , v 1 , u 2 , v 2 )
T and dr = (dx, dy, dz) T . Equation (3.2) is thus written:

u = A • dr (3.3)
which can be solved by inverting the matrix A:

dr = (A T A) -1 A T * •u (3.4)
The metric displacement dr can be expressed as a function of the geodetic unit vectors (t u , t v ) that orientates the sensor and the geodetic unit vector perpendicular to the image plane n. They can therefore be given by the Euler angles. dr = ( dr. t u ) t u + ( dr. tv ) tv + ( dr. n) n (3.5)

By applying Thales' theorem to the system, we can express the projection of the metric displacements on one of the unit axes using the physical size of the pixel (l i , l j ):

( dr. t u ) = l j * u ( dr. tv ) = l i * v (3.6) Equation (3.6
) is therefore written as follows: (3.7)

         u = dr x *
Chapter 3

A stereoscopic model for landslide analysis: Application to the Montgombert landslide Applying Equation (3.7) to the two cameras, the result is that the matrix A can be written as follows:

A =               t u1x l j1 t u1y l j1 t u1z l j1 t v1x l i1 t v1y l i1 t v1z l i1 t u2x l j2 t u2y l j2 t u2z l j2 t v2x l i2 t v2y l i2 t v2z l i2               (3.8)
This system is over-constrained since we have 4 pixel-displacement inputs to determine only 3 metricdisplacement outputs. It was attempted to take only 3 out of 4 values to estimate the contribution of each input. The 4x4 matrix thus becomes a 3x3 matrix which is written according to 4 possible combinations.

It turns out that, in our specific case, the matrix A does not always invert since it was too similar to a singular matrix. We have therefore decided to keep the 4x4 matrix A for the rest of the computations to ensure some stability in the calculation along the time series.

The 3D displacement grids are generated from the correlation results filtered by the correlation coefficient, by the standard deviation of the stable parts and by the morphological filter (see Chapter 2).

To complete the filtered results and thus obtain a dense deformation field, we used the Inverse Distance Weighted (IDW) method which allows us to use an adaptive step interpolation. The IDW method is based on the fact that the points furthest away from the missing value have little influence on it while the closest points have the greatest influence (Chen and Liu, 2012). It is a deterministic method that is local (Adhikary and Dash, 2017). The advantage of using this method is that a missing value can be interpolated either by setting the number of points to be used (variable radius of the search circle) or by setting the maximum radius of the search circle (variable number of input points to be used). In our study, we opted to use the IDW method so that the number of points used in the calculation is fixed and is chosen by the user in order to adapt the interpolation for different study cases. In order to avoid using completely unrelated information in areas with a poor input density, the maximum radius is also chosen by the user.

Let us consider Q the empty cell, Q i the neighboring cell and i the variable from 1 to the number of cells taken into account, n. This last parameter n, ranges from a minimum number n min to a maximum required number n max . The interpolated value v(Q) is calculated if within a maximum radius there is the minimum number of points required. It is defined as Equation (3.9) :

v(Q) =        ∑ n i=1 w p i v i ∑ n i=1 w p i w i = 1 d(Q, Q i ) (3.9)
where w i is the weight and p the IDW power value. The larger the power value is, the more influence the close points have on the result of the interpolation. In other words, the larger the power values is, the more localized the interpolation is. The selected interpolation method implies three parameters that are defined by the user: the minimum number of points to consider, the maximum radius of the search circle and the power parameter. As often happens when it comes to interpolating, there is no straightforward guide to choose the power setting. It must be adapted for each situation since it influences the smoothing of the results. The minimum number of points is by default (in our study) equal to 3 and the power parameter is set to one. The maximum radius of the search circle is by default set to 30 pixels.

Landslide displacement estimation

In order to ensure the validity of the results and to estimate the associated quality and uncertainty, the time series of the two cameras were first processed independently with the TSM toolbox (see Figure 3.7).

The correlation results, corrected for the camera motion, were then used for the 3D reconstruction. To remain consistent, the same prior calibration was applied in all three cases. Finally, the obtained time series were compared to the in-situ data acquired using the total station.

Monoscopic displacement fields

The images were processed with the CMC approach (Common Master Correlation) (Desrues et al., 2019) and the first image, dating from February 11 2020, was used as the master image. 12 images were processed between February 11 and April 20 2020. This correlation strategy was adopted to minimize the influence of the device on the calculations by integrating a DEM acquired at a date prior to the time series. The stable area that is taken into account corresponds to the apparent bedrock. Figure 3.8 shows an example of 3D displacements calculated for each camera between February 11 and March 13 2020.

Displacements are interpolated according to the interpolation method presented previously in section 2.4.

The displacements calculated by the two monoscopic approaches show a good agreement between them that could indicate that they may be independent of the two positions of the cameras (Figure 3.8). To compare both monoscopic results, we fit the distribution of the difference between the results obtained by the right camera and those obtained by the left camera by a normal distribution. The mean and the standard deviation associated to this distribution are respectively -1.5 cm and 9 cm for the EW component, -2 cm and 10 cm for the NS component and 0.8 cm and 6 cm for the Up component. These results show that the displacements calculated by the two monoscopic approaches are of the same order of magnitude. The main differences can be explained by the geographical position of the camera in relation to the object under study, as suggested by Gance et al. (2014). It affects the stable part localization where radiometric properties are stable along time. Stable areas have to be distributed in order to be as spatially homogeneous as possible in the image plane. It also influences the line of sight of the camera that must be optimized in order to limit errors due to incorrectly estimated distortion parameters, and the back-lighting influence. The green color refers to the stereoscopic approach. The term Module refers to the TSM processing (Desrues et al., 2019) introduced in Chapter 2. The term Motion time series stands for the 3D displacements expressed in a geodetic reference system. Finally, the process DEM back projection is situated after Module 1 because, in our case, images are processed according to a Common Master approach (CMC). If not, the process should be performed after Module 3. .

Stereoscopic displacement fields

The evolution of the pattern of the 3D cumulative displacements is shown in Figure 3.9. From February 11 to March 3, the movement is relatively stable with velocities of less than one centimeter per day. Only a small part, at the base of the slope, is in motion with an average movement rate of 5 cm.day -1 . A deceleration is then observed in this zone on February 25. From March 6, a reactivation of the localized zone marked by velocities greater than 8 cm.day -1 and a destabilization of the lower part of the movement (2-3 cm.day -1 ) are observed. This unstable zone is highlighted by a clear boundary: the lateral limits coincide with the observed limits and the upper limit is set back from the observed scarp. The movement then evolves homogeneously from March 6. A clear acceleration is observed until March 14, followed by a deceleration of the movement which suggests that the slope reaches a new state of stability. The upper part of the slide is slightly impacted since it has average speed of about 2 cm.day -1 .

We can also notice that there is no information on the lower part of the slope, which, in particular, includes the apparent bedrock. This is due to the filtering of the data based on the residual displacements of the stable parts. At the base of the localized part which is in motion, the lack of information is, in this case, due to the loss of coherence from one image to another. This suggests a rapid movement or a flow type deformation that induces a loss of soil texture. Furthermore, the footprint of the information was constrained by the interpolation that was performed on a reduced part of the image. This is why the lateral edges are cut off. In addition, the direction of the 3D displacements reconstructed from a stereoscopic approach is consistent with that of the topographic targets measured by total station as we can see in Figure 3.10 for a set of nine topographic targets. The azimuth represented is that of the displacements calculated over the period February 11, 2020 -April 20, 2020 (period of the time-lapse photography). A difference of about 3 • separates the average directions of the two measurements. The directions of the displacements are quite similar for each of the topographic targets which suggests a spatially homogeneous deformation.

Let us now consider the local normal and the tangential displacements that are shown in Figure 3.11 with thumbnails from C to H. We can see that a normal deformation particularly affects the central part from the beginning of March. This deformation increases over time but not homogeneously over the whole surface. This deformation draws, around targets 23 and 26, several bands perpendicular to the preferential direction of the tangential displacements. On the other hand, it is homogeneous at the level of target 22.

Comparison of the stereoscopic approach with in-situ measurements and the monoscopic approach

The results from the 3D reconstruction are compared to some in-situ measurements, realized by a total station. Figure 3.12 shows the cumulative displacements for the 9 topographic targets located in Figure 3.9 in thumbnail A. There are three types of trends shown in Figure 3.12:

• A trend towards a stable state in which the measured displacements do not change over time. This is the case of target 12 which is located in the part which was already considered stable according to field observations, and which is thus confirmed by the photographic and in-situ measurements.

On the five known-to-be stable targets, the standard deviation of the displacements, calculated from the photography time lapse, is less than 1 cm for the EW component (i.e., approximately along the width direction of the slide), is less than 5 cm for the NS component (along to the length direction)

and is less than 4 cm for the Up component. • A trend of moderate homogeneous movement. Targets 17, 4, and 8 were affected by displacements of several centimeters during the event but are still in function. The error of the photography measurements relative to the total station measurements, expressed by the RMSE (Eq. (2.2)), is calculated over 13 targets which show a trend of moderate homogeneous movement. The error varies from 0.7 to 3.2 cm for the EW direction, from 2.7 to 6.5 cm for the NS direction and from 1.4 to 4.4 cm for the up direction. There is thus an average difference of about 3-4 + / -2.7 cm between the total station and photographic measurements. This is rather satisfactory for a low-cost method providing 3D displacement fields only on pairs of images.

• A trend of significant movement that results in the loss of signal from some topographic targets and/or their physical destruction. This is the case for the targets located in the lower part of the movement. The contribution of time-lapse photography, in these cases, is extremely important, since we no longer have in-situ information on the behavior of the environment in these areas: Was there a failure? Has there been a change in the predominant type of deformation (from sliding to toppling, from sliding to mudflow, for example)? The results of the 3D reconstruction shows that the environment has been destabilized without loss of coherence to reach a new state of stability.

Figure 3.13 shows three cumulative displacements (for targets 8, 17 and 6) according to the three components EW, NS and Up. These targets were here chosen in order to illustrate the accuracy of the measurements through different behaviours: target 8 was displaced in all three directions in the order of several centimeters; target 17 underwent greater displacements but kept a low velocity less than 0.5 m.day -1 ; and, finally, target 26 moved by about one meter. We can highlight that:

• First, the results from the monoscopic system shows an average dispersion in position, calculated over the five known-to-be stable targets, equaling, respectively for the right and the left cameras, 2.6 cm and 1.6 cm for the EW component, 4.0 cm and 2.2 for the NS component and 2.5 cm and 0.9 cm for the Up component. The RMSE is computed by taking into account photography displacements and those measured by the total station along all of the time series for each topographic target (i.e., 13 available targets). The results show that the RMSE associated with the EW component ranges from 1 cm to 13 cm, the one associated with the NS component varies from 2 cm to 24 cm and the one associated with the Up component ranges from 4 cm to 17 cm.

• Secondly, the results from the stereoscopic system shows an average dispersion in position, calculated over the same five known-to-be stable targets, equals 0.5 cm for the EW component, 0.7 cm for the NS component and 0.2 cm for the Up component, which is an improvement compared to the monoscopic approach. In addition, they are in agreement with the in-situ data. The calculated RMSE computed on the EW component equals 1.7 cm, the one computed on the NS component equals 3.9 cm and the one computed on the Up component equals 2.6 cm.

• Finally, the plots that are relative to target 26, the large displacements resulting from the 3D reconstruction are overestimated compared to those calculated by the monoscopic approach for the NS and Up components. This deviation is induced by several factors such as the downsampling of the DEM in the image matching process or the non-optimal parameters of the cameras, but it also can be influenced by the type of deformations which is monitored.

Discussion

Source of potential errors

The reconstruction of the 3D displacements involves a DEM during the image matching process but not during the inversion of the displacements. This implies that the 3D displacement computation depends on the geometry in which the scene is contained. In other words, we eliminate the DEM in the computation but the results remain dependent on the DEM. The errors that potentially influence the calculation are therefore: the DEM resolution in the image matching process and the estimation of the camera parameters.

Image matching and DEM resolution

The image matching is performed by using a DEM backprojected in the image plane. To evaluate the influence of the resolution of the DEM (Figure 3.14), we have generated five DEMs at resolutions (valid for both directions): 100 cm, 50 cm, 25 cm, 15 cm and 5 cm. Each corresponding coordinate grid was interpolated linearly and the image matching technique was applied.

Two RMSE were computed for each matching image, either from the left to the right image or vice versa. They have been calculated on the 9 ground control points (GCPs) visible on Figure 3.6 on both the geographical coordinates (X,Y ,Z) and the pixel coordinates (u,v). The impact of the DEM resolution on the matching image is relatively small since the RMSE calculated on the pixel coordinates varies at most 2 pixels. The RMSE calculated on the geographic coordinates varies at most by 2 cm. A lower variability of the results is however noted in the right-to-left direction. If other GCPs had been available, it would have been interesting to estimate the influence of the number of GCPs in addition to the resolution of the DEM. However, a more detailed study would have more impact on an object whose surface presents counter-slopes or different preferential orientations.

Calibration parameters and camera position The computation of the 3D displacements is based on the external orientation of the cameras that defined the matrix A (Equation 3.8). The estimation of external and internal properties is directly dependent on the accuracy of the calibration. However, these parameters may be subject to biases induced by variations of the environmental conditions (e.g., temperature, humidity) or by camera motion. In our study, the calibration was carried out at the beginning of the installation of the cameras but was not renewed. In the case of long time series, it would be interesting to re-estimate these parameters after periods of rain and/or cold.

Moreover, we make the strong hypothesis that the physical pixel size (l i , l j ), which depends on the camera focal length and the camera-object distance, remains constant for the computation. Indeed this statement is valid for small displacements where the physical pixel size does not change significantly from the initial position to the moved pixel position. This may be not true for larger displacements.

For displacements larger than 500 mm, the physical pixel size varies from a fraction of a centimeter.

Consequently, displacements could be overestimated or underestimated by a few centimeters.

Type of displacements and computation precision

As introduced previously, the technique of 3D displacement construction by stereoscopic method has the major advantage of being able to reconstruct tangential surface and normal displacements, as opposed to the monoscopic method which loses the information related to the normal displacements by assuming that all displacements are tangential. Figure 3.11 shows that the normal displacements mainly affect the lower part of the landslide in which the calculated displacements are largest. But if we look at the normal to tangential ratios along time, we can see, in our case, that the tangential displacements predominate the motion. However, even if the normal displacements are small, they are far from negligible since they represent about 23 % of the total displacement over the full period. This therefore suggests that the monoscopic approach may truncate the results in the case where the normal components are significant.

In order to be able to conclude on the observed bias in the large displacements calculated by stereoscopy compared to the monoscopic approach, further investigations could lead to a better understanding.

Landslide deformation estimation

Landslide deformation analysis is complementary to the displacements analysis since it uses the surface displacements to generate mechanical information. Deformation analysis was performed using the 3D displacement reconstructed from the stereophotography.

Deformation analysis

A deformation analysis is carried out in order to highlight different mechanical behaviors within the landslide such as shear zones or compression/extension zones. In the case of the Montgombert landslide few significant changes are visible by photo-interpretation due to the dense vegetation, even though trees in the upper parts of the landslide continue to indicate the activity of the movement due to their tilting.

Strain analysis

Let us consider, u and v as the displacements along the x, y directions in the local reference system such that the resulting displacement is planimetric. The Cauchy 2D deformation tensor, ε, is then written:

ε =      ∂u x ∂x 1 2 ( ∂u x ∂x + ∂u y ∂y ) 1 2 ( ∂u x ∂x + ∂u y ∂y ) ∂u y ∂y      (3.10)
where the diagonal components, A stereoscopic model for landslide analysis: Application to the Montgombert landslide tensor ε, involves two important points: (1) the displacements u and v are considered to be small and (2) the displacements are in the form of a regular grid resolved so that the results are filtered out of outliers to avoid error propagation and to increase the accuracy of the results. To do so, the least-square approach that consists in fitting a linear plane locally and inverting its parameters allows to build a more robust gradient and it is generally used (Pan et al., 2009;Travelletti et al., 2014;Shen et al., 2015).

To take into account the local variations of the slope geomporhology (aspect and slope), we have chosen to consider the displacements/velocities expressed in the local tangent-normal reference system. Locally, displacements are projected on a plane centered on a point. The local divergence is then calculated for the specific point. The divergence corresponds to the trace of the deformation tensor; it results that the slide is extended when the divergence is positive and is in compression when the divergence is negative.

A zero value indicates a homogeneous motion. As the displacement grids have been interpolated and are therefore regular with a certain predefined step size in the previous steps, we have directly applied a central discrete operator for the calculation of the deformation. It is defined as follows:

div( v) = v 2 (i, j + n) -v 2 (i, j -n) 2n∆ j + v 1 (i + n, j) -v 1 (i -n, j) 2n∆i (3.11)
where (v 1 , v 2 ) are the projected displacements in the local tangential plane centered in (i, j), n is the step and ∆i, ∆ j are the lattice steps of the grid.

Results

The deformation calculation module implemented in TSM proposes two approaches, both realized in the tangential plane: a first approach which consists in calculating the deformation tensor and the associated eigenvalues. This assumes that the deformations are small and represent the majority of 3D displacements.

The surface deformation, ε, is then represented and is equal to the sum of the eigenvalues: ε = ε 1 + ε 2 (Travelletti et al., 2014) and the shear strain γ is defined as γ = |ε 1ε 2 |. In the second approach, the deformation is directly calculated as a divergence. In our case the deformations calculated in the tangential plane represent about 77% of the total displacements (after the event by considering the full surface).

These displacements remain large and therefore will not be considered as small displacements. This is why we preferred to turn to an analysis of the divergence of the displacement field which still allows us to highlight the compression and extension zones. • the movement can be considered spatially homogeneous in the slide except in the most active zone;

• a top limit is observed below the scarp observed from 2018. It is mainly formed by extension zones but also by an alternation of extension-compression features. These zones have the particularity of being perpendicular to the local tangential movement;

Thickness estimation

In this section we propose to determine the sliding depth from a low-cost and affordable system based on essentially kinematic considerations. We used a similar approach to the one proposed by Booth et al.

(2013) by considering the 3D displacements computed with optical image correlation and by invoking the material rheology. However, our methodology does not compute the thickness by calibrating it with in-situ measurements. To consider a more complex geometry of the sliding surface and to get closer to a more realistic geometry close to the area where the sliding surface rises to the ground surface, we introduce a geometrical parameter (i.e., disbonding parameter) presented below.

Landslide and conservation of mass

Considering that the landslide is similar to an in-compressible fluid, sliding along a slope parallel to the surface (Figure 3.16a), the conservation of the mass is written:

-h(x) ū(x) + h(x + dx) ū(x + dx) + u top n = 0 (3.12)
where u top n is the normal component of the surface deformation field, h is the thickness of the unstable mass and ū is the average of the velocity inside the unstable mass (Figure 3.16a). Equation (3.12) can be rewritten as follows:

u top n + ∇.(h ū) = 0 (3.13)
Assuming that the thickness is constant over the length of the moving layer, Eq. (3.13) becomes: Conservation of mass applied in the case of a presence of a zone of disbonding, in which h is variable and the normal velocity at the bottom is non-zero (Eq. (3.17)). θ is the angle between the surface of direction vector u tang bottom the disbonding surface of directing vector u slip . δ is the slope angle.

h = - u top n ∇.( ū) (3.14)
ū can also be expressed as being proportional to the surface velocity that corresponds to the tangential component of the deformation field measured at the surface, u top tang , by integrating a factor α such as:

ū = αu top tang (3.15)
The parameter α is a constant related to the landslide rheology. As explained by Booth et al. (2013), this parameter is very poorly quantified and the assumption that it is a constant implies that the rheology of the movement is spatially uniform. A value equal to 1 indicates that the moving mass is similar to a rigid block. A value between 2/3 and 1 indicates motion is consistent with that of a plug flow, a value equal to 2/3 to a viscous Newtonian fluid and a value equal to 1/2 is consistent with a linear vertical velocity profile.

If the sliding surface is not plane and shows a disbonding zone which results in the presence of a non-zero velocity at the base of the slide, called from now interface (Figure 3.16b), Equation (3.12) then becomes:

-h(x) ū(x) + h(x + dx) ū(x + dx) + u top n (x)dx -u bottom n (x)dx = 0 (3.16)
where u bottom n corresponds to the displacement at the interface along the direction of u top n . It means that u bottom n is not related to the direction of the normal to the interface. By integrating Eq. (3.15), Equation (3.16) is written as follows:

αh∇.u top tang + u top tang ∇.(αh) = u bottom n -u top n (3.17)
By considering that the variations at the interface are negligible with time, the relationship between displacement u bottom n and the one along the direction of the interface, u slip is:

u bottom n = -u slip ∇.(h) (3.18)
Let us introduce the parameter β related to the zone of disbonding being the ratio of the slide velocity at the interface to the tangential velocity at the surface:

β = u slip u top tang (3.19)
Thus, the relation between the normal velocity at the bottom and the tangential velocity at the top can be written as follows:

u bottom n = -β u top tang ∇.(h) (3.20)
Assuming that α is constant and h varies along the layer, equation 3.17 then becomes:

αh∇.u top tang + (α + β )u top tang ∇.(h) = -u top n (3.21)
To solve this equation, we decided to apply an iterative system that can only be valid if u bottom n is small enough. To do this, we rewrote equation 3.21 as a set of n equations using a perturbation approach so that

h = h 0 + h 1 + ... + h n with h n-1 << h n :                              h 0 = - u top n α∇.(u top tang ) h 1 = - u top n + (α + β )u top tang ∇.(h 0 ) α∇.(u top tang ) ... h n = - u top n + (α + β )u top tang ∇.(h n-1 ) α∇.(u top tang ) (3.22)
This set of equations is implemented in TSM so that the user sets the α and β parameters and the number of iterations n from which the solution converges. A more complex solution that varies the pair of parameters in a stochastic way would also be interesting to implement. This system is solved by minimizing the RMS calculated from the values h n and h n-1 :

RMS = 1 n × ∑ (h n -h n-1 ) 2 , (3.23)
This relation is valid only if the divergence of tangential velocity at the surface is significant. Indeed, low values (which reflect a homogeneous deformation) would lead to an overestimation of the values of h. It is therefore necessary to have at least significant divergence values in order to get consistent results of the landslide thickness. For the next calculation, the divergence of the tangential field was filtered from its low values by considering any value less than two standard deviations as homogeneous deformation.

Figure 3.17: Profiles (located in Figure 3.18) of the ground surface (black line) and profiles of the sliding surface (red dots). The thickness was calculated without taking into account a disbonding zone. α is equal to 1 which refers to a single sliding block. It was computed from the correlation results of the pair-wise images dating from February 11, 2020 and April 13, 2020. N, S, W, E stand for North, South, West, East respectively.

Results and discussion

Thickness calculated between February 11, 2020 and April 13, 2020 Three profiles, outlined in the preferred tangential direction and located in Figure 3.18, are shown in Figure 3.17. The thickness values have been averaged in an ellipse set by the user in order to be able to quickly average in a preferred direction. Here, the values were averaged in a radius of 1 m around the central point and extracted so that there is one point every meter along the profile (user-specified parameter). We used α = 1 with no disbonding zone as a first approximation. The maximum thickness calculated along the profiles is 6.3 m. Remember that the thickness is measured perpendicular to the topographic surface. A sliding lobe profile is drawn along the P3-P4 profile between points Q4 and Q5. However, it is not visible on profile P5-P6 because too few points are transferred there. From point Q6 to the end of the profile in P5, a small thickness of h (< 1 m) is calculated and echoes the highly deformed zone in which compression and extension alternates. We can therefore argue that from point Q6 to P5, the zone has undergone a very superficial movement which may be similar to a translational slide.

The thickness values were also estimated by using several time intervals in which the event occurred. For time intervals larger than 31 days (i.e., thumbnails from H to L), the estimated thicknesses show a similar distribution. For time intervals less than 31 days, the estimated thicknesses shows a finer distribution.

Further analysis will lead to the analysis of the stability of the moving mass.

Thickness sensitivity Two a priori unknown parameters are introduced in the equation of the calculation of the thickness h (Equation 3.21): the α parameter which is a parameter related to the rheology of the material and thus to its mechanical behavior, and the β parameter which is a geometrical parameter related to the disbonding zone. In order to evaluate the influence of these parameters on the thickness estimation, we were interested in the P4-P5 profile in which the sliding lobe is the most visible. We then modified the parameters one after the other. Figure 3.19 shows the influence of the rheological parameter α on the maximum and average thickness calculated along the P3-P4 profile without taking into account any zone of disbonding (Equation 3.14). As we can see the parameter α has an important role on the calculation of the thickness and thus on the estimation of the volume involved (section 4.3). A difference of 6 m is calculated between a mass with a linear velocity profile and a mass assumed to be a rigid sliding block.

However, whatever the α parameter chosen, the spatial pattern is the same as also observed by Booth et al. (2013). This means that even if the measurement cannot be constrained by an a priori on α, it will still be possible to spatially highlight the most active, shallowest and deepest areas, thus allowing the landslide to be mapped. The influence of the parameter β has been tested on some points of the P3-P4

profile which seem to belong to the upper part of the sliding surface. Indeed, this parameter only matters when there is a large depth gradient. The depth of this zone was calculated at different values of parameter β and at α arbitrarily equal to 1. β ranges from 0 to 1 to be in accordance with the fact that the tangential displacement at the surface is higher than the displacement along the interface. A value close to zero therefore means that there is no or little discontinuity (i.e., θ ≈ 0, Figure 3.16b). As shown in Figure 3.20, few modifications are introduced by β but they confirm the physical meaning of this extra term. The introduction of this refinement is only for a more complete method. As explained, this parameter only matters in particular cases, as close to regions where the scarp is cutting the surface.

Interpolation The calculation proposed in TSM is a discrete calculation which happens to be very sensitive to noise. Several interpolation and smoothing methods are proposed in the literature such as polynomial fitting (Tang et al., 2020) or statistical operator such as Krigging method. However, they are not integrated by default in the calculation of thickness in TSM. Indeed, the calculation of the sliding surface depth requires high divergence values in order not to integrate extreme values that would have no real meaning. This implies that the calculation is performed in areas where the deformations are not homogeneous. Therefore, interpolating values could lead to misinterpretation.

Used in-situ information The results presented above were realized with a set of α and β parameters that designate a specific behavior and geometry. Finding out which parameter is best suited to the data and the study site requires external data such as geophysical or geotechnical data even if it is often very difficult to obtain those point data. (3) a shallow thickness is consistent with an interface rising to the surface. In addition, the piezometric survey revealed the presence of a groundwater level at a depth of 18 m and a resurgence was noticed at Q1. As observed by Mathy and Lorier (2018), the material is mainly composed of mica shale scree. When the latter is saturated with water it gives rise to regular falls of material. The movements detected at the bottom of point Q4 and thus of point Q1 could then be detected as a result of surface changes due to water infiltration. It is therefore consistent with a shallow sliding.

Volume estimation

Estimating the volume of a landslide or rock mass is important for anticipating disaster scenarios by quantifying the potential damage and implications. Several techniques have been developed and approved.

Geometrical relationships

The geometrical relations are based on the geometrical shape of the sliding surface in the case of purely rotational movement. The most common way to calculate a volume is to consider the landslide as having a half ellipsoid shape (Cruden and Varnes, 1996). Let us consider V e the volume of the moving mass, h max the maximum depth of the landslide surface perpendicular to the surface topography, L and W the length and width of a half ellipsoid along the slope. The relationship between these variables is expressed in the following manner:

V e = 1 6 π.h max .W.L (3.24)
By converting these variables into horizontal values, equation 3.24 becomes:

V e ≈ 2 3 .A h .z max (3.25)
where z max is the maximum depth along the vertical assuming that z max ≈ h max /cos(β ) with β the average slope angle. A h is the surface area defined as A h = π.L h .W h which are respectively the horizontal length and width as shown in Figure 3.22. Another geometrical shape can also be taken into account, that of an elliptical paraboloid which has no quasi vertical limits. This is an approach that is often considered to be the most suitable and realistic (Jaboyedoff et al., 2019;Meier et al., 2020). Let V p be the volume calculated from an elliptical paraboloid shape. The volume is defined as:

V p = 1 8 .π.W h .L h .z max (3.26)
which can be approximated as:

V p ≈ 1 2 .π.A h .z max (
V t = π.A.B.(h - h 3 3C 2 ) (3.28)
Chapter 3

A stereoscopic model for landslide analysis: Application to the Montgombert landslide where h is the sliding surface depth, A and B the two major axes of the ellipsoid in the XY plane, and C is the major axis in the Z direction.

Empirical relationship between the volume and the surface area

The volume can also be calculated by empirical laws that relate it to the surface area. This method is based on landslide erosion rates and it is governed by the power law below (Larsen et al., 2010;Nikolaeva et al., 2014):

V r = η.A γ (3.29)
where η is the intercept and γ is the scaling component. We can already see that this relationship, although geometric, does not depend on the topography or the heterogeneity of the material under consideration (Guzzetti et al., 2008;Nikolaeva et al., 2014). Measurements carried out on a number of data sets (Larsen et al., 2010) have highlighted two ranges of values for γ, the first of which is [1.1-1.3] and applies to shallow soil-based landslides. The second value range is [1.3-1.6] and applies to landslides with bedrock failure. It is noted that a small variation in γ can drastically change the volume estimate. Considering the empirical relationships, we used several γ and η constants derived from power laws established by different authors (Table 3.1) so that the area of the slip surface is integrated into the interval that the authors used. The surface area calculated is equal to 3619 m 2 . Figure 3.23 displays the volumes estimated from the different methods. Measurements with a γ constant less than 1.3 underestimate the volume compared to geometric methods. Only the power laws whose gamma is greater than 1.3 estimate with the same order of magnitude as the geometrical measurements. According to Larsen et al. (2010), a surface area of the order of 10 3 m 2 with a moving mass thickness of several meters would correspond to a bedrock landslide scar. From observations and knowledge of the terrain, the landslide that affected the Montgombert landslide in 2020 is more similar to a bedrock landslide deposits.

Conclusions

This chapter was devoted to the development of two new methods allowing: (1) to estimate the 3D displacement fields using a steresoscopic system and the image correlation technique applied to time series of photographs, and (2) to estimate the deformation and the depth of the sliding mass only from surface imaging.

The displacements calculated either from a monoscopic approach or from a stereoscopic system, are consistent with in-situ data. The two monosocopic approaches, with different angle of views and distance to the slope, gave results to within ten centimeters. The position of the camera and the estimation of the intrinsic camera parameters may explain the difference. In addition, we have noticed that the metric displacements reconstructed by stereoscopy are overestimated compared to those calculated by monoscopic approach when the displacements are large (greater than 500 mm). A study conducted purely on synthetic data would be relevant in order to control all the input parameters of the processing and understand the overestimation. To better understand the landslide behavior, we carried out an analysis of the deformation field performed through the tangential vector fields and the components normal to the local surfaces. It showed that the major advantage of the stereoscopic technique is to be able to take into account the tangential and normal displacements as opposed to the monoscopic approach which only considers displacements that are mostly tangential to the surface. The normal displacements being not negligible in our case, the monoscopic method proves to be less adapted than the stereoscopic method.

Determining the depth field by a technique based on terrestrial optical tools is, by contrast with boreholes, much lighter, and relies on easily affordable material compared to geophysical monitoring. We proposed a new methodology to estimate the thickness of a moving mass by considering an in-compressible fluid. Two parameters come into play: the rheological factor α and the geometrical parameter β . The α parameter determines the velocity profile along the slope and controls depth values. The β parameter seems to have little influence in our case, assuming that a disbonding zone may not be appropriate here. A maximum depth of about 6 m is calculated for α equal to 1 (i.e., sliding rigid block) whereas it reach 12 m for α equal to 1/2 (i.e., linear vertical velocity profile). These values remains of the order of those estimated by geotechnical profiles. The volume involved in the March 2020 acceleration was calculated using existing geometrical and empirical methods. Considering a γ greater than 1.3, an averaged volume was calculated and is 12000 m 3 .

This method showed the potential of using terrestrial optical imaging for the understanding of landslide behavior. Nerveless, to demonstrate better the influence of the β parameter, further analysis should be performed either by using synthetic data or by considering a slide for which there are information, and a priori, on slope deformation.

Chapter 4

Pre-and post-event monitoring analysis: application to the Cliets rockslide (Savoie, French Alps)

This chapter presents a multi-technique analysis applied to the Cliets rockslide. The analysis is conducted by combining geological observations, meteorological data analysis, active seismics, topographic measurements and simple physical modelling. We implement post-processing tools to perform analysis of the relationship between the displacement and rainfall and to perform kinematics and mechanical analysis. Additionally, we present the application of TSM to image time series acquired with a monsocopic camera during the pre-failure stage and the post-event monitoring.

This chapter is divided in two sections: a first section deals with the pre-event study of the movement of Cliets using the multi-technique approach. The in-situ data, meteorological data, seismological and geophysical data were acquired in an operational context and were generated by the SAGE company in collaboration with the Savoie Departmental Council.

The second section deals with the photography time lapse monitoring that was carried out in parallel with the pre-event study. We present the data acquired at the different stages of the evolution of the rockslide (event-monitoring, failure monitoring and post-event monitoring) and the type of information we can extract from it.

This first section is based on a submitted article entitled: Desrues M, Malet JP, Brenguier O, Carrier A, Mathy A, Lorier L (2020) Landslide dynamics inferred from in-situ measurements: the Cliets rockslide (Savoie, French Alps).

In this section, the term purge refers to the removal of unstable blocks and masses without causing destabilisation of the surrounding environment.

1 Case study in the context of monitoring and early-warning

Introduction

Gravitational instabilities affecting hard rocks have various forms such as topples, rockslides and rockfalls (Varnes, 1978;Hoek and Bray, 1981;Hutchinson, J.N., 1988). Toppling is mainly observed in large slopes in response to Quaternary glacial retreat and slope stress redistribution, river erosion and thus slope debuttressing (Grämiger et al., 2017). Toppling initiated by glacial retreat was observed at the La Clapière landslide (France) which involved about 50 millions m 3 of material (Follacci et al., 1988;Vengeon, 1998), at the Siwalik Hills (Nepal) (Tamrakar et al., 2002) or at the Mount Breakenridge area (Canada) (Nichol et al., 2002).

This gravitational phenomena occurs on rock with regularly spaced discontinuities whose dip is opposed to the slope. This geological structure has the particularity of forming a series of columns. Goodman and Bray (1976) (3) a block flexure toppling which is defined by a pseudo-continuous flexure along columns divided by numerous joints. Column stability is then controlled by the slope geometry, the rock properties and the extent of the discontinuities (Rotaru et al., 2007). Consequently toppling can occur if the force applied by the center of gravity overhangs the base of the columns while sliding can occur if the thrust becomes higher than the ratio of normal strength applied to the columns and if it exceeds the angle of friction (Goodman and Bray, 1976;Wyllie, 1980).

Generally, toppling are analyzed using numerical modelling to analyse the column stability based on Goodman and Bray's limit equilibrium technique (Goodman and Bray, 1976). It suggests a discretization of the columns which leads us to consider them as free to rotate around their base (Nichol, 2000;Nichol et al., 2002;[START_REF] Merrien-Soukatchoff | Modélisation par éléments distincts du phénomène de fauchage gravitaire. application au glissement de la Clapière (Saint-Etienne-de-Tinée, Alpes-Maritimes)[END_REF]Pereira et al., 2013). An analytical solution is proposed by Guo et al. (2017) for a toppling failure triggered by an earthquake. Here we are only interested in the contribution of a kinematic study in the understanding of toppling movements which is the common way to monitor an instability in an operational context.

Slope motion can reveal the predominant type of movement. For sliding-type processes, surface velocity indicates an overturning of the block topples according to the toppled beds (Cruden and Varnes, 1996;Nichol et al., 2002;Glueer et al., 2019) whereas for toppling processes, surface velocity depends on the geometry of the surface of rupture (Glastonbury and Fell, 2010;Glueer et al., 2019). Several forcing factors control the slope destabilization possibly leading to another kinematic regime of even catastrophic failures (Cruden and Varnes, 1996;Rose and Hungr, 2007;Oppikofer et al., 2008;Jesus et al., 2018).

Among them are the geological effects driven by earthquakes (material dilatation leading to fast water Chapter 4

Pre-and post-event monitoring analysis ingress), water-level changes (rainfall, snowmelt) and by the presence of discontinuities or faults; the human effects such as deforestation or drainage. These events are reflected in changes in trend when looking at the evolution of displacement velocity over time. These events can be periodical (e.g annual or seasonal) or sudden (extreme events) (Rose and Hungr, 2007;Guo et al., 2017). When the velocity trend is clear and is only slightly affected by external factors, we can hypothesize that the motion is governed by a mechanical action in which external destabilizing factors play only a small role.

Forecasting the time to failure is crucial to take decisions such as closing a road or evacuating citizens. A robust approach is proposed by Voight (1989) considering that surface displacements present a vertical asymptote close to the failure time. It suggests that deformation evolves linearly with time and becomes independent of external factors (Voight, 1989;Favre et al., 1992;Rose and Hungr, 2007).

Generally, if slope materials undergo deformation, changes of material properties or changes in slope geometry, the measured displacements take the form of a standard creep curve (Petley, 2004;Federico et al., 2012;Crosta and Agliardi, 2003). This curve is divided into three stages, based on the strain-time relationships for the creep of materials: the first stage corresponds to strain hardening, the second stage corresponds with a steady state where velocity is constant and the third stage corresponds to accelerating creep. These stages are summarized by introducing the notions of primary, secondary and tertiary creep (Emery, 1978). Recently, Scoppettuolo et al. (2020) propose another approach based on the partitioning of the displacement curve in time according to triggering factors. 4 trends are highlighted: the first one corresponds to the linear behavior, the second one to a transition phase between two stable stages, the third trend is a response of an occasional event and the last trend is the failure stage.

To understand the mechanisms involved which depend on the geometry of the instability and the material properties, physical approaches are considered (Mirgon et al., 1993). A physical basis is applied to simulate the kinematics such as viscoplastic models for clayey materials (Angeli et al., 1996;Van Asch and Malet, 2009;Picarelli et al., 2008) and rate-weakening friction models (Helmstetter et al., 2004;Handwerger et al., 2016).

In this work, we analyse the mechanical behavior of the Cliets unstable slope which is representative of gravitational instabilities observed in the micaschists of the Arly Valley, Central French Alps (Kasperski, 2008). This unstable slope is very interesting because several acceleration phases were observed before the final acceleration toward failure in a toppling-sliding mode of motion (Mathy andLorier, 2013, 2018). The surface displacement time series is analyzed for the last two failures which occurred on 23 January 2014 and 9 February 2019, and is compared to the long-term topographic measurement data sets acquired since 2004. A simple mechanical analysis is proposed by using a combination of friction models and forecasting methods.

We first present the geological setting of the unstable slope (section 1.2) and the observations collected in-situ (section 1.3). The role of, respectively, the geological structures and the triggering factors is analyzed (section 1.4). Time to failure is then analyzed (section 1.5) and is followed by a more general discussion (section 1.6).

1.2 Geological setting of the study area (Pothérat, 2005;Lorier et al., 2019)). This narrow valley is the object of many landslides, regular and sometimes very important which can affect the road and the Arly riverbed. The first documented event (the 'Moulin Ravier' landslide) occurred in 1955 involving about 17 millions m 3 (Goguel, 1957). Since then, several events generally affecting the RD1212 road have been recorded (Baudin, 2015) According to the 1/50000 • geological map (Figure 4.1), the left bank of the Arly valley is located in the western branch of the Belledonne crystalline massif. The "satin series" is composed of varied micaschists with carbonaceous schists coming from an old metamorphosed detrital series. Because of its mineralogical composition marked by the presence of phyllite minerals, the rock is very brittle and very altered in the presence of water (Pothérat, 2005;Baudin, 2015). The right bank of the Arly valley is composed of sedimentary coal sedimentary rocks that are above these crystalline rocks (Jeannin, 2001;Dussauge-Peisser et al., 2002).

The Cliets unstable slope in the Arly Valley

The Cliets unstable slope (45.77 • , 6.486 • ) from 650 to 775 m in elevation is crossed by a 60 m long tunnel (Figure 4.2) in order to protect the road. Many rockfall events have been documented in December 1996 (1500 m 3 ), October 2003 (3000 m 3 ), January 2014 (9000 m 3 ) and February 2019 (10000 m 3 ). The source areas of these rockfalls are delineated in Figure 4.2 (Kasperski, 2008;Pothérat, 2005;Mathy and Lorier, 2013). The area extends over a length of about 72 m and a width of 80 m. To protect the road, geotechnical screws were drilled in 2003 and rock purges were carried out between 2003 and 2004.

Local geology and discontinuity of the Cliets unstable slope

The Cliets unstable slope is composed of micaschist which formed from an old metamorphic folded detritic sedimentary series (cf 1.2.1). The alternation of sandstone with black micaschist provides a mix of hardly consolidated and soft rocks sensitive to internal deformation and weathering. The rock masses are subject to toppling as shown in Figure 4.2. The sub-vertical schistosity tilts by gravity towards West and the dip is towards East with an average of 45 • . Consequently, the deep deformation induces a failure surface constituting a "rock topple -rock slide" complex (Cruden et al., 1993). The area concerned by toppling is estimated to be around 30 m width, 40 m long along the line of largest slope and on average 10 m depth (Mathy and Lorier, 2018). Furthermore, no surface flow networks are observed for the periods 2004 -2014 and 2018 -2019 thus suggesting that water infiltrates directly into the sub-vertical discontinuities.

1.3 Data (Pothérat, 2005;Kasperski, 2008).

Long-term monitoring of the Cliets provides a very precise vision of the activity of the rock movement and the various hazards show that the road remains subject and vulnerable to the threat of landslides of several thousand cubic meters. In our analysis of the two last events of the Cliets rockslide (periods 3 and 6 according to Table 4.1), we relied on geological and geophysical measurements, on topographical data and data that would indicate triggering conditions.

Chapter 4

Preand post-event monitoring analysis

Rock discontinuity detection and mapping

Geological observations The geological model of the slope is inferred from field observations of the discontinuities, stereoscopic image analysis and geophysical surveys (Goguel, 1957;Pothérat, 2005;Kasperski, 2008;Mathy and Lorier, 2013;Baudin, 2015;Mathy and Lorier, 2018). Rock discontinuities were mapped at different periods in 2004, 2013 -2014 and 2018 and are presented in section 1.4.1. Three major zones are highlighted:

• an active zone uphill where the two rockfalls of 1996 and 2003 occurred. In 2013, the presence of voids are noticed and observations concluded to a constant toppling motion. Benchmark 4 is situated in this zone (topographic network of the period 3 -see Figure 1.3.3).

• a rocky spur where the last rockfalls occurred. Considered as a stable area in 2013, some discontinuities are noticed in 2018 and deformation of the geotechnical screws (benchmark 12 -topographic network of the period 6) are observed.

• an active zone downhill where visible cracks are observed (benchmarks 6 and 7 -topographic network of the period 3 ). In 2018, an unstable area is observed close to the rocky spur.

Seismic survey In 2018, a geophysical survey has been conducted by the SAGE company to determine the depth of the soft rock prone to rock falls (Mathy and Lorier, 2018). Three seismic profiles (P1, P2 and P3) of 60 ml were carried out using hammer blow as seismic source. The data were processed using the classical seismic refraction method (plus-minus method) (Figure 4.3). Analysis of the P-wave velocity highlights a seismic interface between 10 to 15 m among 730 and 745 m of altitude and between 5 to 8 m among 700 and 730 m of altitude. This layer is characterized by velocities close to 800 m.s -1 .

Topographic measurements

Displacements of 10 benchmarks have been monitored at a monthly frequency from the 7 September 2004 to the 31 January 2014 (period 3). The array of benchmark is shown in and rock slide events. The 24-h meteorological stations Col des Saisies and Albertville, available for the period 2004 -2018 are used for the analysis. For the period 2018 -2019, the 1-h meteorological station

Les Cliets is added. To integrate all data, hourly data are converted to daily data and a weighted average is calculated by the Inverse Distance Weighting method. Table 4.2 summarizes the characteristics to each station. The Mean Annual Precipitation (MAP) associated with the total weighted rainfall (R w ) for the period 2004 -2019 is 1227 mm. The effective rainfall is calculated according to Oudin et al. (2005) in order to estimate the contribution of meteorological water. The water balance is written as:

R E = R w -ET P -Ru max (4.1)
where R E is the effective rainfall (mm), R w is the total weighted rainfall (mm), ET P is the evapotranspiration (mm) and Ru max is the ground water storage (mm). The MAP of the weighted effective rainfall is 646 mm. We can point out that the rockslides at the beginning of the year (Jan. 2014 and Feb. 2019) precede a year with a calculated annual effective rainfall of more than 800 mm.

Water infiltration can also be related to snowmelt. The degree-day method (Van Mullem et al., 2004) is used to approximate snowmelt fluxes at daily time steps (Moussav et al., 1989;Kustas et al., 1994;Hock, 2003):

M = C(t * -T b ) (4.2)
where M is the total melt rate (cm.day -1 ), C is the melt factor expressed (cm.day -1 . • C), t * is the average daily temperature ( • C) and T b is the base temperature (usually equals to 0 • C). The melt factor varies (Kustas et al., 1994) and has a high spatial and temporal variability. According to Moussav et al. (1989), the melt factor is taken to 0.47 cm.day -1 . • C -1 . The snowmelt effective rainfall is calculated adding the melt contribution to the rainfall excess.

Seismic forcing is estimated by analysing the temporal occurrence of earthquakes for the period 2004 -2019 in a 80 km 2 area centered on the unstable slope (Engels and Grunberg, 2013) (Figure 4.4b)).

The maximum earthquake magnitude recorded is 3.6 close to Albertville (Figure 4.4b)) in October 2015.

During the first period (Dec. 2014 to Jan. 2014) 48 earthquakes of any magnitude are recorded among 28 with ML > 2. During the second period (Aug. 2018 to Feb. 2019) 12 seismic events of any magnitude are recorded among 3 with ML > 2. The mean of the magnitudes over the all periods reaches 1.8 and the distance to the Cliets landslide is in average 16 km. Based on these results, no direct relationship can be made between seismic and gravitational events and we will continue with the rainfall-motion relationship.

1.4 Role of geological structures and of the triggers 1.4.1 Predisposing factors: geological structures From a regional point of view, three major families of structures were highlighted by stereo photo interpretation (Kasperski, 2008): ( 1) faults (dip direction) N50 -60 • on the right bank of the Arly river,

(2) incidents N110 -120 • of greater length and (3) faults N80 • which are transverse to the Arly river. Kasperski (2008) also noticed a fault N20 -30 • which runs the Arly river in the east and passes upstream of the unstable zone.

From a local point of view, according to the field observations three main joint sets are observed (Kasperski, 2008;Mathy andLorier, 2013, 2018) and are summarized in Figure 4.5:

• S0 dips sub-vertically towards the SE (strike -dip angle: N30-90 • ). This discontinuity set matches the schistosity and correlates with the direction of the main scarp. Its direction is consistent with those of the regional faults.

• F1 dips towards the W with mostly directions of N140-50 • at the base of the slope and N140-90 • at the top.

• F2 dips towards the S (N70-70 • ). This set was only observed during the period 3. Its direction gets closer to regional faults.

• D1 dips towards W (N142-46 • ). This discontinuity is observed as of 2018.

Further several discontinuities were observed in 2013 prior to the rockslide of 2014 in Figure 4.5a (Mathy and Lorier, 2018). There are summarized with the notation Dn and are documented in the stereonet (Figure 4.5b). The green plan is related to the toppling foliation. The stereonet representation highlights two major families of discontinuities: F2 and the majority of Dn from the first family of discontinuities observed in 2013. It indicates that the preferred sliding direction is SSE. The second family of discontinuities is the one formed by F2 and D1 (observed only in 2018). Their direction is specific to the slip of the Cliets in view of the regional structures. The sliding direction is therefore preferably WSW and follows the direction of the greatest slope.

Triggering factors: rainfall and gravitational events

Searching for correlation among rainfall and the occurrence of gravitational events is mostly carried out on catchment to regional scales (Keefer, 1994;Zêzere et al., 1999;Polemio and Petrucci, 2000;Marc et al., 2018). This leads to the definition of rainfall thresholds in terms of high-intensity rainfalls or long-lasting rainfall and can be expressed by several rain properties (intensity, duration, antecedent rainfall over several time scales; (Zêzere et al., 1999;Guzzetti et al., 2008;Giannecchini et al., 2012;Teja et al., 2019). The same approach is used to establish some links between earthquake magnitude and landslide occurrences (Keefer, 1994;Rodrıguez et al., 1999;Marc et al., 2017;Martha et al., 2017). At the slope scale, linking forcing factors to slope instability requires us to introduce a new definition for a gravitational event. We thus define a gravitational event as a positive surface velocity change between two dates. The following method and criteria are used to define what we called here a change motion: (1) in-situ velocity observations are smoothed using a moving average; (2) knee points are detected in the smoothed time series using Satopaa's algorithm (Satopaa et al., 2011). The knee points are calculated over a moving window and are always part of a convex portion of the velocity curve with an increasing slope;

(3) criteria are applied to discriminate local knee points using ∑ k i=k-n/2 a(i) < ∑ k+n/2 i=k a(i) where n is the length of the window in which the acceleration a is calculated. n depends on the data sampling. With this method, there is no discrimination of the results on the value of the calculated accelerations around the knee point.

The length of the moving window depends on the data resolution. For period 6, the length of the window is 30 samples which corresponds to a 30-day window. The window is shifted by 12 days in order to have overlap and the possibility to detect other knee points. Figure 4.6 shows the events detected by the knee point method for all targets for the period 6. An event is detected within 2 days, which corresponds to the uncertainty of the method in our case related to the size of the moving window and the shift. We can see that there are a total of 6 events detected around: October 24 2018, November 20 and 29, December 16, January 9 2019 and January 28. A spatial representation of the date of the event shows us that: 47% of the targets have a concomitant change in velocity (targets surrounding the rocky spur), 12% of the targets are within the uncertainty of the nearest two days, 41% of the targets have no change in velocity detected in the 5 days before and after November 29. For period 3, the frequency of the data is not regular. The most appropriate combination of parameters proved to be: a sample size of 13 samples for the moving window and a step size of 5 samples. As shown in Figure 4.6 several events were not detected and this is mainly due to the lack of data. Nevertheless a dozen events are highlighted from 2005 to 2014 including at least three at the end of 2013 (visible on all targets).

Cross-correlation techniques are used to analyze the relation between effective rainfall and surface velocity. Figure 4.7 shows the results of the cross-correlation (expressed with the Pearson coefficient) for benchmarks whose events are the most marked for period 3 (benchmarks 6 and 7) and for benchmarks with successful event detection for period 6 (benchmarks 12 and 36). The cross-correlation is realized over a rolling window which size corresponds to 10 samples. Dashed lines indicate gravitational events determined by the knee-point method. No relationship emerges for period 3; at the opposite for period 6, each gravitational event is associated with a progressive rise in the Pearson coefficient up to nearly 1.0. This suggests that motion is controlled by daily rainfall rather than by antecedent rainfall. These observations are similar if we consider the daily gross rainfall or the snowmelt equivalent.

Analysis of slope behaviour

As already identified, the gravitational events of period 6 are mostly related to daily effective rainfall.

However, a general trend in the surface velocity profile is identified for both previous periods independently of the velocity changes generated by the triggering factor (Rose and Hungr, 2007). This is why we are interested in analysing the slope behavior from a kinematic and mechanical points of view.

Forecasting landslide failure

Considering a ductile material in conditions where stress and temperature are constant (e.g. third creep stage), the time to slope failure can be determined, at first stage, by the relation linking the inverse velocity to time (Fukuzono, 1985):

1 V = [A(α -1)] 1 1-α (t f -t) 1 α-1 (4.3)
where A and α are constants, t f is the predicted time to slope failure and t is time. In the case where α = 2, the fit matches with a linear trend. This simple approach often provides reliable forecasts when used shortly before failure. (Rose and Hungr, 2007;Voight, 1989).

Here we used the surface velocity time series to predict time to slope failure (Crosta and Agliardi, 2003;Carlà et al., 2018). As mentioned by Carlà et al. (2017), the predicted time to failure does not represent an exact date but a possible interval time. The forecasts are impacted by the parameterization of the time series smoothing. Here, the frequency acquisition is different from one campaign to another which requires that data are filtered for the period 3 by applying a centered moving average windows of width n = 5 and raw data are used for the period 6 thanks to the data quality.

For clarity, Figure 4.8 shows the inverse velocity for the benchmarks 6 and 7 for the period 3 and benchmarks 12, 29 and 36 for the period 6 (Figure 4.8). The prediction realized from data from 2011-03-06 to 2013-07-07, corresponding to a period of two and a half years, underestimates the time to failure (2013-12-04 for benchmark 6 and 2014-01-11 for benchmark 7). For the period 6, only a forecast realized a week before the failure gives a relevant estimate of the time to failure within three days and this, for the 19 benchmarks for which an acceleration is observed. A prediction on the same number of days for the two periods could unfortunately not be made since the number of data is different from one period to the other (e.g., five values for 48 days for period 3 and one value per day for period 6).

The period predicted for the failure ranges over 3 days from 2019-01-31 to 2019-02-09. This estimated period of 3 days is narrow, and is similar for the 19 benchmarks suggesting that all portions of the slope behave similarly from a kinematic point of view.

Physical model of slope behaviour

The kinematic analysis indicates the possibility to predict a realistic time interval for the slope failure. This suggests that motion is governed primarily by gravity and that the influence of external forcing factors is similar over time. To explain the acceleration of the displacements, we apply a simple physical model based on sliding instability. This model permits us to differentiate a stable slow slide from an unstable fast slide, knowing that a slope rupture occurs when the sliding surface is characterized by a "rate-weakening" friction (Handwerger et al., 2016).

Clustering In order to model slope behaviour with a physical model, first, clusters of benchmarks with similar motion and kinematic trend (e.g. surface displacement magnitude, azimutal displacement direction, ratio of horizontal vs. vertical displacement, range of velocity and acceleration) and corresponding to kinematic sub-units have to be identified. To define robust clusters, a dense network of benchmarks (at least 10 samples) is compulsory. For the period 6, 5 clusters and thus 5 kinematic units are identified from the statistical analysis of the motion of 37 benchmarks (Figure 4.9). All clusters show a RMSE < 10 -3 . It can be concluded that, for the period between the 1 st December 2018 and the day before failure: (1) the first cluster (red color in Figure 4.9) corresponds to the upper part of the rocky spur defined in Section 1.2.3 characterized with average velocity of 7 mm.day -1 ; (2) the second cluster (yellow) corresponds to the area surrounding the rocky spur and is characterized by an averaged velocity of 2 mm.day -1 ; (3) the third cluster (lime) corresponds to the lower part of the spur characterized with an average velocity of 0.6 mm.day -1 ; (4) the fourth cluster (dark green) corresponds to conditionally stable units with an averaged velocity of 0.6 mm.day -1 . The main difference between this cluster and the third cluster results in the shape of the velocity curves before the rupture. The third cluster shows an acceleration of displacements one week before the rupture date while the fourth cluster does not, confirming its stability over time. Finally, (5) the fifth cluster (blue) corresponds to stable units with an averaged velocity of the order of 10 -5 mm.day -1 . This simple spatial and temporal classification reveals that the top limit of the first cluster corresponds to the fracture sets D1 observed in 2018 (Section 4.5). The lack of observations surrounding this area is consistent with stable zones which suggests that the motion is highly localized. Additionally, having three clusters in the collapsed area also suggests that the motion has been initialized at the top.

Plunge analysis As mentioned in section 1.2.3, the Cliets unstable slope presents two mechanical behaviors of toppling and sliding along a plane. In order to estimate if one behavior predominates over time, we use the displacement vectors to calculate the plunge angle. As detailed in Figure 4.10, the plunge angle corresponds to the angle between the horizon and the displacement vector. This is a simplified method compared to the plunge computation of Glueer et al. (2019) who used structural information such as the slope dip direction, the friction angle and the toppling depths. For a toppling mode, the plunge angle should increase with time, while for a sliding mode, the plunge angle should be constant over time.

The plunge angle is calculated for 4 benchmarks (period 3) and for the 3 kinematic clusters (period 6) of high velocities (Figure 4.10). At first order, two patterns are emerging:

• For the period 3, the plunge angle of benchmarks 6 and 7 increases until 2010 and then becomes constant. This indicates a change of toppling to sliding over time. This trend is not sufficiently pronounced for benchmarks 4 and 8 to assume that sliding is predominant over toppling. This is consistent with observations realized in July 2013 (see section 1.3.2, (Mathy and Lorier, 2013)).

• For the period 6, the plunge angle decreases until December 2018 and then becomes relatively constant. This decrease may correspond to a transition between toppling and sliding along an initial foliation acting as sliding discontinuity.

Frictional rate and state model of the unstable slope Considering the period where the plunge angle is constant, we assume that the unstable slope is sliding on a plane over time and that the motion can be assumed to a single block sliding on an inclined plane. Consequently, we apply a frictional model based on the Dieterich-Ruina frictional law ( (Dieterich, 2007) which relates the solid friction coefficient Chapter 4

Pre-and post-event monitoring analysis

Discussion

To determine the contribution of modelling, two assessments were conducted: ( 1 The question to answer is to identify if the stability of the system is modified by each gravitational event. Segmentation is realized according to the main gravitational events presented in section 1.4.2. Two models are applied for the period 3 and three models for the period UAV photogrammetry of 20 cm of resolution. Images were converted into brightness images. The sub vertical cliffs were used to correct the camera movements. The camera was calibrated according to the none apriori method using 30 topographic targets as ground control points whose geodetic positions were known (section 2.1). The associated RMSE xy is of the order of 3.1 pixels. The distance between the camera and the object ranges from 160 to 225 m implying that the ground pixel size varies from 2.4 to 3.6 cm.

Images were processed according to the CMC and V SC methods (the correlation parameters were adapted to the method used). Because of low motion in the time lapse interval, small movements are considered as noise by the V SC method and are therefore filtered out, while they are detected by the CMC method (Figure 4.16). This is highlighted by the amount of information detected for the targets 36, 29, 18, 8 and 1 over time. A pixel recurrence analysis corresponding to the number of times the pixel has significant information (detected and unfiltered) was performed (Figure 4.17). Even if the amount of information differs depending on the method used, it can be seen that all the information is clustered in the central part of the landslide. Using the V SC method, if we assume that it detects the most significant displacements, we can highlight two sub-areas: one located uphill (pixel recurrence > 14) and one located downhill (pixel recurrence < 10). This difference can be interpreted as an indication of the dynamics of the motion: a high pixel recurrence in the VSC method is similar to a fast motion whose displacements measured by correlation are significant along the time series; a lower recurrence is similar to a state in which the displacements were significant at a given date in the time series.

Visual image interpretation

The fact that an image is not exploitable for the calculation of displacements does not necessarily imply that it is useless. Indeed, optical images are also often used for image interpretation in order to identify weather conditions (cloud detection in the case of satellite images, snow and rain detection for terrestrial, aerial and satellite images (Bossu et al., 2011)). Snow cover is also a key factor in understanding the dynamics of climate change and water resources. In the case of the study of the Cliets rockslide, several images show significant snow cover that varies over time. Assessing the presence of snow cover can be an important information. Indeed, during thawing, water runoff can also penetrate the ground and add additional water to the slope. If the slope is saturated, the environment is unstable and the conditions are checked for destabilization of the environment.

Determining the percentage of pixels in an image that correspond to snow is constrained by the illumination of the object over time since the radiometric properties of the pixel can vary from one image to another. There are several calculation methods to estimate whether a pixel represents snow or not using threshold methods. The snow pixel has a specific radiometric signature since the color white has a high digital number (DN). Salvatori et al. (2011) assumes, by default, that pixels with a DN greater than 127 correspond to snow. We can distinguish: the snow-to-snow method, which analyzes the blue component of the image since it is in this channel that the reflectance of the snow is the highest (Salvatori et al., 2011); the RGB Normalized Difference Snow Index (RGBNDSI) method developed by Hinkler et al. (2002) which is based on the calculation of a band derived from the R, G, B bands and the empirical approximation of the MIR band (Mid-Infrared); and finally the Gaussian Mixture Model (GMM) method (Rüfenacht et al., 2014). These three methods have been applied to digital camera images in the context of snow monitoring in mountainous areas. To approximate the percentage of snow in an image, we applied the snow-to-snow method by normalizing the images ( We can plot the percentage of snow pixels identified during the time series. The associated curve What can be seen is that the slope failure is preceded by a significant snow melt episode in which the snow cover decreases by about 45 %.

Image processing during the event of February 9, 2019

The event of February 9, 2019 of the Cliets rockslide led to the rupture of the slope and the fall of around 10000 m 3 of rocks. Fortunately, this event occurred after the road at the base of the slope was closed. The correlation of the images associated with this event highlights the spatial and kinematic evolution of the destabilization of the slope (Figure 4.20). As we have seen previously, the unstable zone is very well delineated and is detached from the rest of the slope by moderate displacements. During the February event, we can see that only the central rocky spur is in movement and that the stable part is not affected, at least in the short term. On date February 9, from 8 am to 9 am, displacements of more than 4 pixels are detected in the upper zone of the unstable area. A velocity gradient from top to bottom is also observed.

From 9 am to 12 am the movements are amplified and we notice that some areas are already de-correlated (gray areas). This de-correlation indicates a substantial change of texture from one image to another which suggests that the movements were too important (moved or collapsed material). Between noon and 1 pm, the movement accelerates but becomes homogeneous within the unstable area. In section 1.5.2

we have seen that the predominant first-order behavior was slippage from December 2019. Looking at a finer scale, we see that the velocity gradient detected at the beginning of the day (between 8 and 11 am)

suggests that the movement is much more complex than what happens at first order. Between 1 and 2 pm

Conclusions

In this chapter we have covered the different stages of a multi-technique and modelling monitoring with the study of the Cliets rockslide.

We carry out a pre-event analysis by conducting a kinematic and mechanistic analysis of unstable slope in order to characterize the influence of gravitational events on the kinematic regime and define how to forecast time to failure for complex slope movement associating toppling and sliding. A set of topographic measurements over a period of several years (1996 -2019) is used to understand the long-term kinematics and define possible relationships with triggering factors. The analysis was realized for the two last gravitational events which occurred the 23 rd January 2014 in the downhill section of the slope and the 9 th February 2019 impacting the central section (rocky spur) of the slope. The main specificity of this slope is the mechanism involved for two rock failures linked to the structure of the slope composed of schists and nearly vertical discontinuities subject to both toppling and sliding (sections 1.2.3 and 1.5).

The kinematic analysis indicates that there are no direct relationships between slope activity and environmental forcing factors (effective rainfall, snowmelt and crustal seismicity). However, a clear velocity trend is observed and the plunge analysis indicate that the slope has been affected by a transition from toppling to sliding (section 1.5.2). A rate and state frictional model suggests that velocity trends are dominant and consistent with the behavior of the material (section 1.5.2) and that rain and seismicity plays a second role in destabilizing the slope.

This study demonstrates that the velocity analysis may be applied to complex rockslides by providing information on both the dynamics of the movement and the prediction of time to failure. As shown in sections 1.5.1 and 1.6, the prediction analysis is constrained by the density of available data close to the rupture date. We suggest performing a prediction analysis by considering only the last values of the time Chapter 4

Pre-and post-event monitoring analysis series, thus avoiding biasing the results by adding data. The application of a physical model to explain the displacements requires a good quality of data close to a potential slope failure. The advantage of applying this model to a dense set of benchmarks is to be able to spatially compartmentalise the whole area in terms of movement dynamics. This method may be applied to all types of time series of displacements from different platforms such as GPS or ground-based interferometry radar (GB-InSAR). It therefore opens many perspectives in the use of new measurement techniques for the study and understanding of the dynamics of rockslide movements.

The time-lapse photography unfortunately did not allow us to complete this pre-event study from a kinematic point of view. However, the results of the processing of the time series of images are promising and suggest that a mechanical and kinematic study could have been carried out on these data. Time-lapse photography proved to be a complementary tool to classical monitoring techniques since it provided a visual image interpretation during the pre-event and event. It has become indispensable for post-event studies in order to understand in-situ techniques that could not be used after the event. As a result of this study, it is shown that the Cliets rockslide has reached a new state of stability in which it evolves in a linear way. The post-event movements were only detected in the central collapsed zone, which supposes

(1) that there are no zones destabilized by the event, (2) the movements detected in the central zone are linked to residual surface movements induced strongly by surface erosion.

Chapter 5

Conclusions and perspectives

consider a simple mechanical model with an incompressible media moving on an inclined plane which sliding surface can present a disbonding zone. The method, based on the law of mass conservation (Bishop, 1999;Booth et al., 2013), has shown, in the case of the Montgombert landslide, the strong influence of the rheology parameter on the estimated depth values and the weak influence of the geometrical parameter on the geometry of the sliding surface. The inverted depths are in the same order of magnitude as those measured by geotechnical methods. The method is however dependent on the possible evolution of the sliding geometry (for instance with the retrogression of the slide) and to the fact that the movements only affect the subsurface.

In order to complete the post-processing analysis tools, we have proposed (Chapter 4) several tools for the kinematic and mechanical landslide analysis. For this purpose, we used rainfall-displacement relationships (Guzzetti et al., 2008;Giannecchini et al., 2012) for estimating the time to failure (Voight, 1989;Rose and Hungr, 2007) and a friction model to determine the kinematic regime of the landslide (Helmstetter et al., 2004;Sornette et al., 2004). The tools are applied to the dataset of the Les Cliets landslide, (Savoie, France). The tools are adapted both for the analysis of displacement time series acquired both from classical geodetic techniques (total stations, GNSS) or imaging techniques (terrestrial photographs, ground-based InSAR). The methods allow to detect and cluster kinematic areas of homogeneous behavior (displacement amplitude and direction) and the transition between mechanical behavior (in this case, from toppling to sliding). In addition, the analysis of the time-lapse photographs acquired at several development stages of the landslide allowed us to estimate the snow cover and melting influence on the slope evolution ten days before the failure, to delineate the extent of the sliding before the failure, and to visualize and thus better understand the progressive development of the failure, and finally to detect and map possible residual movements after the main failure. In this latter application, the analysis of the displacement fields allow us to estimate that the slope reached a stability plateau two months after the failure.

Several scientific and operational questions were addressed in this work (see Chapter 1) and a summary is proposed below.

On the influence of the instrument setting and parameterization

The acquisition parameters may vary since they relate to the technical specifications of the camera, to the orientation and position of the instrument in relation to the object, and to the acquisition mode such as the measurement frequency or the image modality. All these parameters are interrelated and will constrain the type of motion that can be detected (Chapter 2, Chapter 3): the spatial coverage, the quality of the image co-registration (linked to the stability of the pole/concrete pillar, and to the percentage of stable areas in the image plane), the size of the pixel on the ground (and therefore the accuracy of the measurement), the internal calibration of the camera. In addition, we showed that defining the acquisition approach (monoscopic or stereoscopic) is of importance (Chapter 3) if the normal and tangential components need to be detected, or if many zones are masked from one unique viewpoint. The analysis of the image time series resulting from these two approaches Chapter 5

Pre-and post-event monitoring analysis processing workflow developed for time series of optical images GDM-OPT from ForM@Ter and MPIC-OPT from ESA/GEP (Provost et al., 2020)), and 2) as an on-demand service in which the operator can control most of the processing options, thus allowing us to obtain the most accurate information for a site (Chapter 2). Implementing TSM as fully automated service (e.g. operational surveillance mode) or as on-demand service (e.g. science mode) is an objective of both the research lab and the SAGE Society in the coming months.

Perspectives

2.1 Technical and science development of time-lapse optical imagery for landslide research

The TSM-based analysis of optical image time series acquired on several use cases calls for further works on instrument, processing and post-processing code developments:

In terms of instrument developments:

• The growing development of stand-alone acquisition systems extends the number of sensors and acquisition devices available. A more in-depth analysis on the contribution of webcams (and in general low-cost IoT based camera sensors) for the generation of displacement fields could offer an interesting alternative in operational contexts in which webcams are currently only used for visual interpretation.

• The use of infrared cameras, recently successfully applied in various domains such as volcano surveillance (Spampinato et al., 2011;Guerin et al., 2019), rockfall monitoring (Le Roy, 2020) or landslide monitoring (Frodella et al., 2017) should be considered. This approach paves the way to new application such as the indirect estimation of surface temperature (and thus heat transfer and friction) and soil/rock humidity owing to increased slope characterization and understanding of the pre-failure/failure landslide stages.

• The use of triggered-cameras should be considered in order to either store/analyse sequences of image time series useful, or to change automatically the frequency of measurements in case of changes in the slope regime. For instance, Hibert et al. (2015) used two seismologicallytriggered video cameras to extract images only when a seismic event associated to rockfalls are recorded. Similarly, cameras may be triggered with meteorological information in order to adapt the acquisition frequency to the purpose of the monitoring.

• The cameras could be linked to on-board computers such as RaspberryPi (Santise et al., 2017), with the possibility of a real-time (or near-real-time) onboard processing with TSM; this would allow us to transfer the model outputs only (e.g., displacement fields) and not the images.

• Hierarchical approaches could be considered for the image time series analysis by processing, first, images at a low resolutions for a first evaluation of the landslide kinematics. Then, when the (Mazzanti and Pezzetti, 2013)). The usability is defined as follows: "+++" means that optical terrestrial images may provide relevant results and are useful for landslide understanding. "++" means that optical terrestrial images are useful but some requirements have to be considered. "+" means that the data and technique could be applicable but only on specific cases. Possible monitoring alternatives should be considered.

Use of terrestrial optical time series for operational landslide applications

For uptake of the technology (data and processing) by users (e.g. private engineering and consultancy companies, stakeholders in charge of disaster risk management), operational monitoring products derived from terrestrial optical time series should consider the customer operational needs and the budget available.

The monitoring system and data acquisitions should be adapted to the main purpose of the study based on criteria related to the current knowledge (kinematics) of the unstable slope, the available geological and geotechnical information, and the environmental forcing conditions. Choosing the most suitable sensor and acquisition parameters requires consideration of to several criteria (Mazzanti and Pezzetti, 2013) that can be summarized into three items:

• the definition of the main purpose of the surveillance, which can be for instance the slope characterization, its long-term evolution or the prediction of time to failure (i.e. early-warning stage);

• the definition of the monitored parameters which justifies the use of a unique or several instruments according to the number of other (in-situ) parameters being monitored. In some cases, fully automated systems are not necessarily the most suitable;

• the definition of the acquisition period of the monitored parameters which define the acquisition frequency (i.e., continuous real-time or long-term campaign monitoring). As underlined by Mazzanti and Pezzetti (2013), monitoring is not a static procedure. It has to consider field feedback and experts' points of view.

In the case of a monitoring carried out with terrestrial optical sensors, the above-presented three criteria are applicable. Time-lapse photography provides a large sequence of optical images and it can be used for landslide studies where the main purposes are slope characterization by giving access to change detection maps and long-term displacement time series useful to document the pre-failure stage (Figure 5.1).

However optical techniques remain dependent of the meteorological conditions limiting their usability as monitoring sensor for the early-warning stage.

Products that can be generated from fixed cameras are listed in Table 5.1 in relation to operational monitoring purposes, and specific requirements to take into account are addressed. In-situ forcing factor datasets. Quant à la distorsion tangentielle, celle-ci est définie par ses deux premiers coefficients (P 1 , P 2 ) de la manière suivante : δ x = P 1 .(r 2 + 2.x 2 ) + 2P 2 .x.y δ y = P 2 .(r 2 + 2.x 2 ) + 2P 1 .x.y (A.6)

Paramètres externes

Les paramètres externes definissent l'orientation et la localisation de la caméra dans un système de coordonnées géodésiques. En d'autres termes, ils décrivent le mouvement de la caméra autour de la scène.

La relation qui relie le système de coordonnées géodésique (X w ,Y w , Z w ) et le système de coordonnées caméra s'exprime de la manière suivante :

      X c Y c Z c 1       =   R t 0 1   *       X w Y w Z w 1       (A.7)
Où R est la matrice de rotation et t la matrice de translation. La matrice R s'exprime à partir des angles d'Euler (ω, Φ, κ) tels que, dans le système de coordonnées géodésiques, ω tourne autour de l'axe the X; Φ tourne selon l'axe Y axis and κ selon l'axe Z. Le sens de ces angles respecte la règle de la main droite.

La matrice de rotation est donc une combinaison de ces trois angles et s'exprime de la manière suivante: La matrice des erreurs ∆Err x s'écrit dès lors : X

R = R z R y R x (A.
∆Err x =                ∂Resultat x ∂DEM x

  . . . . . . . . . . . . . . . . . . . . 1.4 Monitoring techniques for landslide understanding. . . . . . . . . . . . . . . . . . . . . xiv 1.5 In situ techniques. (a) Photograph of a fixed GPS station located at the valley of the

  Photograph of the image acquisition system set up at the Montgombert landslide (Savoie, France). In this configuration, the camera sensor and the data logger are in separate cases. 1.13 Image matching techniques. (a) Tracking features flowchart from Gance et al. (2014). The automated feature detection uses the properties of the target such as the object size, its geometry and its luminance. (b) Image correlation principle. . . . . . . . . . . . . . . 1.14 Example of correlation results applied to the image time series acquired in the landslide of Super-Sauze. Displacement rates in pixel/day. From Travelletti et al. (2012). . . . . . 1.15 Principle of normalized hierarchical image correlation from Travelletti (2011). . . . . . . 1.16 Location of the four studies sites presented in this work. Photographs are taken from the fixed camera. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Figure 3.2) and red lines show the lateral and top limits of the landslide (see Figure 3.1). The black line shows the forestry track while the brown polygons highlight the apparent bedrock. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

  cases. Possible monitoring alternatives should be considered. . . . . . . . . . . . . . . . 166 A.1 Effet de la distorsion sur une grille régulière[START_REF] Noury | Calibrage de caméra plénoptique à partir des images brutes[END_REF]. a) Distorsion radiale posiive (à gauche) et négative (à droite). b) Distorsion tangentielle. . . . . . . . . . . . IV C.1 Vitesse surfacique exprimée m.day -1 et calculée le long de la série temporelle d'images acquises après la crise de Février (Glissement des Cliets, Savoie, Fr.). Les images onté été acquises entre le 15 Février 2019 et le 4 Octobre 2019. Elles ont été traitées par TSM et corrélées selon l'approche V SC. Les flèches noires indiquent la direction du mouvement. X xxiii List of Tables 1.1 Summary of landslide monitoring methods modified after Malet et al. (2002); Delacourt et al. (2007); Michoud et al. (2010). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Summary of the study sites characteristics presented in these works. Location corresponds to the French departments The average velocity refers to the classification proposed by Cruden and Varnes (1996). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Acquisitions characteristics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Combination of parameters used to process the two time series for respectively the Chambon and the Pas de l'Ours landslides. . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Constants from the surface area -volume power law established by several authors after Guzzetti et al. (2009) and Michel et al. (2020). . . . . . . . . . . . . . . . . . . . . . . . 4.1 Chronology of the site's instrumentation . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Characteristics of the meteorological stations used in the Cliets rockslide analysis: time sampling, distance from the slip and weighting coefficients. . . . . . . . . . . . . . . . . 4.3 Results of the friction model realized on two simulations for the period 3 and three simulations for the period 6. They are expressed by the m and Xi estimated parameters and the NRMSE associated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Summary of the products that can be generated from fixed cameras respect to the purposes of the study and some specific requirements. . . . . . . . . . . . . . . . . . . . . . . . . B.1 Ecart type utilisé dans le calcul du bruit gaussien. Ce bruit est appliqué de manière alétaoire à la variable d'entré. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VII

Figure 1 :

 1 Figure 1: Typical sequences of the cumulative displacements of the ground surface over time after Scoppettuolo et al. (2020). (a) Landslide reactivation due to triggering factors. (b) Landslide evolution to failure from a pre-failure stage.

Figure 3 :

 3 Figure 3: At the top, photographs acquired at the Pas de l'Ours landslide in Aiguilles (Hautes-Alpes, France) from a fixed camera on June 8, 12, 16 and 18, 2018. The camera is a Canon EOS 100D of 18 Mega Pixels (MPix) with a focal length equal to 24 mm.At the bottom, schematic trend of the surface velocity over time that can be extracted from the image analysis. Several questions that can be raised at each specific stage of the landslide evolution are presented.

Figure 3 (

 3 Figure 3 (top) shows an example of a temporal series of terrestrial optical images composed of four photographic acquisitions made from a fixed camera at the Pas de l'Ours landslide (Hautes-Alpes, France).

Figure 4 :

 4 Figure 4: Processing pipeline of the Tracing Surface Motion (TSM) toolbox (Desrues et al., 2019).

Figure 1 . 1 :Figure 1 . 2 :

 1112 Figure 1.1: Schematic surface velocities of a gravitational movement afterMirgon et al. (1993)

Figure 1 .

 1 Figure 1.3: Schematic framework that connect the risk management cycle (orange circle), which is managed by the local authorities, the time-related progress of disasters (white circle), the temporal stages of the monitoring (blue circle) and the observation frequency (gray circle). Modified after Stumpf (2013) and Alexander (2002).

Figure 1 . 4 :

 14 Figure 1.4: Monitoring techniques for landslide understanding.

Figure 1 . 5 :

 15 Figure 1.5: In situ techniques. (a) Photograph of a fixed GPS station located at the valley of the Basento river (Italian Southern Apennines) from (Calcaterra et al., 2012). (b) Photograph of a total station set up in front of the Viella landslide (French Pyrenees) (source: J.-P. Malet). (c) Photograph of an extensometer (fixed site visible) from(Malet et al., 2002). (d) Principle of the direct levelling form(Schultz, 1995).

  Terrestrial and Aerial laser scanning Laser scanning or light detection and ranging (LiDAR), is a measure providing 3D dense point clouds. Measures are realized thanks to an infrared laser pulse, and are based most often on the time-of-flight principle (see Figure1.7). Knowing the speed of light, c ,and the response time δt, it is possible to convert the signal into the distance d between the instrument and the

Figure 1 .

 1 Figure 1.6: Radar principle. a) Radar imaging geometry from Bamler (2000). b) Example of SAR image from Yan et al. (2012). Spatial representation of subsidence over Mexico City (2002 -2007).

Figure 1 . 7 :

 17 Figure 1.7: Laser scanning principle. a) Laser scanning acquisition principle from Jaboyedoff et al. (2012). b) 3D point cloud acquired by TLS at the Chambon landslide in 2016 (France) (Desrues et al., 2018).

  most common algorithm is the Iterative Closest Point algorithm (ICP) successfully used byHebel and Stilla (2007);Travelletti et al. (2014);Bellekens et al. (2015) andMoon et al. (2019).2.3.2 Passive imaging sensorsPassive sensors include optical images from several platforms: satellite, airborne and terrestrial platforms 1.8. Depending on the type of the embedded sensor, images may be acquired in several modes: in mono band (grayscale image), in spectral bands (colored image with three bands in the visible spectre: Red, Green Blue (i.e., RGB image)) and in multi-or hyper-spectral bands (several bands in the visible spectrum with the infrared band).

Figure 1 . 8 :

 18 Figure 1.8: Example of optical images acquired by the three platforms and result from the image correlation technique. (a) Left Panchromatic Pléiades image acquired in October 12, 2012. Right Horizontal displacements (in meter) measured at the La Valette landslide (France) between August 7 and October 5, 2012. From Stumpf et al. (2014). (b) Left Aerial image of the La Clapière landslide (France). Right Horizontal velocity (in cm.day -1 ) calculated between 1995 and 1999. From Delacourt et al. (2007). (c) Left Terrestrial image of the Sanières landslide (France) acquired in July 10, 2015 by a PENTAX K200D with a 21 mm of focal length. Right Surface velocity (in pixel.day -1 ) between July 10 and July 30, 2015.

Figure 1 .Figure 1 .

 11 Figure 1.10: Example of a photography time-lapse at the Argentière glacier. (Left) Photograph of the camera view point. (Right) Displacement rates in pixel/day. From Fallourd et al. (2010)

Figure 1 .

 1 Figure 1.11 relates the ground pixel size to the camera-object distance, according to the focal length, for two particular cameras, a Dinion IP ultra 8000MP 1/2.3" CMOS Progressive Scan webcam (12 Mpix of resolution) and a Canon EOS 100D camera (24 MPix of resolution). These two cameras were installed on the Pas de l'Ours site (see Section 4). The advantage of using a webcam in an operational context is the adjustment of the day/night modes, which allows us to cover a whole time interval. On the other hand, many webcams do not have a fixed focal length. In the case of this webcam the focal length varies from 3.2 mm (120 • ) to 5 mm (71 • ) and a simple screw fixes it. It is therefore subject to potential unwanted variations over time.

Figure 1 .

 1 Figure 1.12: Example of the image acquisition system developed by the SAGE Society. (a) Photograph of the image acquisition system set up at the Pas de l'Ours landslide (Hautes-Alpes, France), after Desrues et al. (2019b). (b) View of the inside of the box shown in (a). (c) Photograph of the image acquisition system set up at the Montgombert landslide (Savoie, France). In this configuration, the camera sensor and the data logger are in separate cases.

Figure 1 .

 1 Figure 1.13: Image matching techniques. (a) Tracking features flowchart from Gance et al. (2014). The automated feature detection uses the properties of the target such as the object size, its geometry and its luminance. (b) Image correlation principle.

Figure 1 .

 1 Figure 1.14: Example of correlation results applied to the image time series acquired in the landslide of Super-Sauze. Displacement rates in pixel/day. From Travelletti et al. (2012).

Figure 1 .

 1 Figure 1.15: Principle of normalized hierarchical image correlation from Travelletti (2011).

  Figure 1.16: Location of the four studies sites presented in this work. Photographs are taken from the fixed camera.

4. 1 1

 11 Figure 1.17: Landslide use cases: (a) The Chambon landslide (Isère, France) with a terrestrial laser scanner derived Digital Surface Model (DSM) of the slope. (b) The Pas de l'Ours landslide (Hautes-Alpes, France) with an airborne LiDAR-derived DSM of the slope. (c) The Montgombert landslide with an airborne LiDAR-derived DSM of the slope. (d) The Cliets landslide with an airborne UAV-derived DSM of the slope. Colored areas corresponds to the footprint of the camera. Two colours are used for c) the Montgombert landslide since two cameras were set up to provide a stereoscopic model. Red lines indicate limits of the movement. Blue circle indicates the position of the camera. Black dots are spot elevations.

Figure 1 .

 1 Figure 1.17: Landslide use cases: (a) The Chambon landslide (Isère, France) with a terrestrial laser scanner derived Digital Surface Model (DSM) of the slope. (b) The Pas de l'Ours landslide (Hautes-Alpes, France) with an airborne LiDAR-derived DSM of the slope. (c) The Montgombert landslide with an airborne LiDAR-derived DSM of the slope. (d) The Cliets landslide with an airborne UAV-derived DSM of the slope. Colored areas corresponds to the footprint of the camera. Two colours are used for c) the Montgombert landslide since two cameras were set up to provide a stereoscopic model. Red lines indicate limits of the movement. Blue circle indicates the position of the camera. Black dots are spot elevations. (cont.)

Figure 2 . 1 .

 21 photogrammetric library MicMac (Pierrot-Deseilligny and Paparoditis, 2006).

. 1 )

 1 where (u, v) are the image-based coordinates, (X,Y, Z) are the world 3-D coordinates (X is the East direction, Y is the North direction and Z is the elevation), (X cam ,Y cam , Z cam ) is the position of the camera in the world coordinate system, ( f u , f v ) is the effective focal length and (m 11 , m 12 , . . . m 33 ) is the orientation matrix which depends on the three Euler rotation angles[START_REF] Heikkila | A four-step camera calibration procedure with implicit image correction[END_REF] Kraus and Waldhäusl,

Figure 2 . 1 :

 21 Figure 2.1: Methodology of the Tracing Surface Motion (TSM) processing toolbox from data collection to analysis of motion time series.

Figure 2 . 2 :

 22 Figure 2.2: Definition of the three coordinate systems implied in the pinhole model.

  Analyzing long image time series (< 30 images) of terrain motion poses the problem of creating subsequences of images of equal-level quality. The images may be altered by the presence of snow, rain, mist or shadow. As the selection of images can be very time-consuming, this module allows identifying automatically the best images according to thresholds in radiometry (red band R, green band, G, and blue band B). The sequential steps of the module are(1) the selection of a reference image referred to as the master image, not altered by meteorological or illumination conditions, (2) the calculation of the statistical properties (standard deviation, median and skew of each band of the master and slave images), and (3) the creation of a new image sequence after comparison of the statistical properties of the images. Image

Figure 2 . 3 :

 23 Figure 2.3: Details of the TSM processing toolbox for the Modules 3 to 6.

  based on the method proposed byTravelletti et al. (2012) (Figure2.2). It consists in(1) calculating the line-of-sight of the portion of the terrain visible from the camera position using a high-resolution Digital Surface Model (DSM)( Čučković, 2016), (2) performing a back-projection of the DSM in the image plane according to Equation (2.1)(Corripio, 2004), (3) interpolating the projected DSM coordinates in the image plane, and (4) re-projecting the displacement fields in the geographic coordinate system.

  4) to be flexible to any type of use cases. The results are compared and discussed in section 5.1.4: • (1) the Common Master Correlation method (CMC) consists in forming all pairs of images from a common master image (Pham et al., 2014) (Figure 2.4a); • (2) the Variable Sequential Correlation method (VSC) consists in correlating the images sequentially along the time series (Figure 2.4b). The correlations are thus performed for each interval [I, I + n]

Figure 2 . 4 :

 24 Figure 2.4: Combination strategies for processing large image dataset: (a) Common Master Correlation method (CMC), (b) Variable Sequential Correlation Method (VSC), (c) Redundant Variable Sequential Correlation Method (RVSC) Each gray square represents an image acquired at day t. h 1 , h 2 and h 3 corresponds to images acquired at hour h. The combination Forward-Backward, (d), is illustrated for the VSC method but can be applied to the other two combinations.

  1.1), of the image matching parameters (section 5.1.2), of the size of the stable part used for the co-registration (section 5.1.3) and of the influence of the correlation technique (section 5.1.4) are analyzed on both datasets. For external parameters, the influences of the DEM accuracy as well as the point picking accuracy are analyzed on a third dataset.5.1 Role of the internal parameters5.1.1 Sensitivity to Image ModalitiesSensitivity to the image modality (grayscale image, brightness image, average of p images, texture-derived image) used as input to the image matching is evaluated:1. the "grayscale image G*" modality corresponds to a neo-band derived using an averaging of the three initial RGB bands with a weighting equal to 1; 2. the "brightness image B*" modality corresponds to the transformation of the RGB images using Craig's formula(Pratt, 2013): I = 0.30 × R + 0.59 × G + 0.11 × B where I is the brightness image and R, G and B are the three bands(Fallourd et al., 2010; Pham et al., 2014); 3. the "average of p images A*" modality corresponds to a stack of p = 3 images used to reduce radiometric noise;

  radiometric enhancement is expressed by the histogram transformations Ha, Hb and Hc which are applied in the initial RGB bands: the Ha transformation corresponds to a linear contrast enhancement by stretching the values from min = 0 to max = 255; the Hb transformation corresponds to a histogram normalization to adjust the contrast; and the Hc transformation corresponds to the application of a Contrast Limited Adaptive Histogram Equalization (CLAHE) to improve the local contrast and highlight natural edges.

Figure 2 . 5 :

 25 Figure 2.5: Sensitivity of co-registration quality to the image modality expressed by the criteria RMSE xy and the criteria percentage of pixels with C > C min at Chambon landslide (a) and Pas de l'Ours landslide (b). G* is the grayscale image modality, B* is the brightness image modality, A* is the average of p images modality and T* is the textured-image modality. Ha, Hb and Hc are the histogram transformations applied on the RGB image.

Figure 2

 2 Figure 2.6: Sensitivity of co-registration quality to the correlation parameters expressed by the criteria percentage of pixels with C > C min at (a) the Chambon landslide and (b) the Pas de l'Ours landslide: (left) Inc (size of the research window), (middle left) SzW (size of the correlation window), (middle right) Cmin (minimum threshold in the correlation coefficients) and (right) Reg (image matching regularization parameter). The range of tested parameter values is indicated in brackets.

Figure 2 . 7 :

 27 Figure 2.7: Sensitivity of the co-registration quality to the size of the stable area. Only pixels located on stable areas are taken into account. (a) Offset distribution along the U axis (left) and along the V axis (right) for the Chambon landslide. (b) Offset distribution along the U axis (left) and along the V axis (right) for the Pas de l'Ours landslide.

Figure 2

 2 Figure 2.8: TSM-derived displacement time series [pix] for specific points at the two landslides, at Chambon landslide (a) and Pas de l'Ours landslide (b). The image combination strategies VSC (colored lines) and CMC (dates of the master and slave images indicated on the graphs) are presented.

Figure 2 . 9 :

 29 Figure 2.9: Comparison of the combination strategies Variable Sequential Correlation method (VSC), Redundant Variable Sequential Correlation method (RVSC) and Common Master Correlation method (CMC) processed in a forward mode. Stability of the criteria RMSExy in time for the Chambon landslide (a) and the Pas de l'Ours landslide (b).

  Figure 2.11: P1 is a longitudinal profile of 137 m in the N-S direction and P2 is a transverse profile of 92 m in the W-E direction (Figure 2.11a). The displacements are presented for the period 8-27 March 2017. The image time series was processed with the VSC method in forward direction. One image per day is used. The profile P1 reveals two active parts from 0 to 20 m and from 80 to 100 m (Figure 2.11b). They have contrasting behavior: the displacement rates of the upper part (from 0 to 20 m) decrease over time while those of the lower part (from 80 to 100 m)

Figure 2 .

 2 Figure 2.10: Displacement rates [cm.day -1 ] measured by TSM from 9 to 12 March 2017 at the Chambon landslide. The continuous white line corresponds to the main scarp and the discontinuous white lines corresponds to the lateral flanks of the landslide.

Figure 2 .

 2 Figure 2.11: Displacement rates [cm.day -1 ] measured by TSM from 8 to 27 March 2017 at the Chambon landslide. (a) Profiles location on a terrestrial laser-scanner derived DSM and an ortho-image from 12 March 2017. (b) Displacement rates [m.day -1 ] along Profile P1. (c) Displacement rates [m.day -1 ] along Profile P2. The displacement rates are calculated over a band of 20 pixels in width.

Figure 2 .

 2 Figure 2.12: Displacement rates [m.day -1 ] measured by TSM from 16 to 18 June 2018 at the Pas de l'Ours landslide. The discontinuous white lines corresponds to the Eastern flank of the landslide.

Figure 2 .

 2 Figure 2.13: Displacement rates [m.day -1 ] measured by TSM from 6 to 22 June 2018 at the Pas de l'Ours landslide. (a) Profiles location on a terrestrial laser-scanner derived DSM and an ortho-image from 12 June 2018. (b) Displacement rates [m.day -1 ] along Profile P1. (c) Displacement rates [m.day -1 ] along Profile P2. The displacement rates are calculated over a band of 20 pixels in width.

Figure 2 .

 2 Figure 2.14: Cumulated displacements [m] of specific points (cf. Figure 2.8 measured by TSM and by a total station at the Pas de l'Ours landslide.

Figure 3 . 1 :

 31 Figure 3.1: The Mongombert landslide: (Left) Topographic map after Mathy and Lorier (2018) and Mathy and Lorier (2019). (Right) Topographic profile corresponding to the black line of the topographic map.

Figure 3

 3 Figure 3.2: Cumulative 3D displacements measured by total station. Targets whose displacements time series are represented on the right, are located on the georeferenced aerial image on the left. The red dots without numbers correspond to a position of a topographic target. The accumulated rainfall comes from a weather station in the Arly valley about 3 km from the Montgombert landslide upstream.

Figure 3

 3 Figure 3.3: A schematic illustration of the epipolar geometry that governs the stereoscopic model.

Figure 3 . 4 :

 34 Figure 3.4: Field configuration for the two-camera model with the position of the two cameras at the Mongtbombert landslide. The b/H ratio is equal to 0.66. Both cameras are imaging the lower part of the landslide which is also the deforested part. The camera footprint corresponding to the Left camera is blue and the one corresponding to the right camera is red. The white frame refers to Figure 3.1 for more details about the boundaries of the landslide (red lines).

Figure 3

 3 Figure 3.5: Image pairing by using the ORB detector. 171 features are detected for the pair-image acquired the February 11 2020. Drawn lines connect the pairs of homologous points. As we can noticed with the red circles and the red rectangular, several features from one image get the same homologous point in the second image and conversely.

Figure 3 . 6 :

 36 Figure 3.6: Reprojection of the points used to determine the intrinsic and extrinsic parameters of the cameras by applying the a priori calibration. Points A, B and C indicate the targets used to construct a mean plane of the object and thus show that the Zhang calibration method can be applied.

Figure 3 . 7 :

 37 Figure 3.7: General workflow applied in the case of the Mongtombert study case. It allows to assess the contribution of the 3D reconstruction relative to a monoscopic approach. The blue and red colors refer to Figure 3.4 and to the monoscopic approach.The green color refers to the stereoscopic approach. The term Module refers to the TSM processing(Desrues et al., 2019) introduced in Chapter 2. The term Motion time series stands for the 3D displacements expressed in a geodetic reference system. Finally, the process DEM back projection is situated after Module 1 because, in our case, images are processed according to a Common Master approach (CMC). If not, the process should be performed after Module 3.

Figure 3

 3 Figure 3.8: Monoscopic approach. Displacements are calculated between February 11 and March 13, 2020. (a) 3D displacements [meter], calculated from the image time series from the two monoscopic cameras. Red lines correspond to scarps and limits of the movement observed in 2019 and 2020 (see Figure 3.4). The black line shows the forestry track while the brown polygons highlight the apparent bedrock. (b) Distribution of the difference between the results obtained by the right camera and those obtained by the left camera.

Figure 3

 3 Figure3.9: 3D displacements, in meters, calculated from the image time series taken from a pair of stereoscopic cameras. Red lines correspond to scarps and limits of the movement observed in 2019 and 2020 (see Figure3.4). The black line shows the forestry track while the brown polygons highlight the apparent bedrock. Red points show the localization of some reference topographic targets.

Figure 3 .

 3 Figure 3.10: Cumulative surface displacements of targets 8, 17, 35, 4, 12, 22, 16, 26, 23 (see Figure 3.9 for the localization) according to their directions. (a) Azimuth calculated vs the displacements measured by total station [mm] over the period from February 11, 2020 to April 20, 2020. The direction is calculated every 7 days. (b) Azimuth calculated vs the displacements derived from the steresocopic model [mm] (over the same period as (a)).

Figure 3 .

 3 Figure 3.11: Normal and tangential displacements from the stereoscopic system corresponding to thumbnails from C to H (Figure 3.9). (a) Displacements of the normal component in meters. (b) Displacements of the tangential component by vector fields. Red lines correspond to scarps and limits of the movement observed in 2019 and 2020 (see Figure 3.4). The black line shows the forestry track while the brown polygons highlight the apparent bedrock. Red points show the localization of some reference topographic targets.

Figure 3 .

 3 Figure 3.12: Cumulative 3D displacements of some topographic targets. Lines correspond to the measurements acquired by the total station. Data were collected on an hourly basis. Squares correspond to the displacements calculated from the stereoscopic image time series. Images were processed with a Common Master Correlation approach. (a) Cumulative displacements corresponding to the all image time series. (b) Zoom of the black square located in (a).

Figure 3 .

 3 Figure 3.13: Cumulative 3D displacements of topographic targets 8, 17 and 26 according to the three components WE, NS and Up. Doted lines correspond to the measurements acquired by the total station. Squares correspond to the displacements calculated from the stereoscopic image time series. Images were processed with a Common Master Correlation approach either in monoscopic or in stereoscopic way.

Figure 3 .

 3 Figure 3.14: Part of the image acquired by the left camera with the projected points of the DEM (red dots). a) DEM resolution 15 cm and b) 50 cm resolution.

  change in length per unit length in the x and y directions respectively. A positive value means an extension. The deformation tensor ε is symmetrical and represents the internal deformation of the moving mass. Therefore, eigenvalues ε 1 and ε 2 can be calculated as well as the eigenvectors whose axes indicate the main strain directions. Applying the deformation Chapter 3

Figure 3 .

 3 Figure 3.15 shows an example of the Montgombert landslide deformation analysis. The results come from the displacement fields between the image of February 11, 2020 and April 13, 2020. They thus include the acceleration of March 2020 in which several meters of displacement were measured by terrestrial photography. The normal component, represented in Figure 3.15a, highlights: the active part of the landslide in which the most important displacements were measured, a zone further uphill which is detached from the landslide part by zero values of the normal, a part to the West of the active part in which

Figure 3 .

 3 Figure 3.15: Deformation analysis performed by TSM applied to the pair of images 2020-02-11 -2020-04-13. a) Normal component of the surface displacements expressed in meter, b) Tangential vector field of the surface displacements according to the slope orientation expressed in degrees, c) Ratio of the tangential component to the normal component taken in absolute value, d) Tangential vector field of the surface displacements according to its divergence. Red dots indicate specific topographic target locations (see Figure 3.2) and red lines show the lateral and top limits of the landslide (see Figure 3.1). The black line shows the forestry track while the brown polygons highlight the apparent bedrock.

Figure 3 .

 3 Figure 3.16: Schematic deformation model where the law of the conservation of mass is applied. (a) Conservation of mass applied to a sliding parallel to the surface (Eq. (3.14)). The thickness h is constant along the layer from x to x + dx. (b)Conservation of mass applied in the case of a presence of a zone of disbonding, in which h is variable and the normal velocity at the bottom is non-zero (Eq. (3.17)). θ is the angle between the surface of direction vector u tang bottom the disbonding surface of directing vector u slip . δ is the slope angle.

Figure 3 .

 3 Figure 3.18: Thickness computation with the parameters α = 1 without including a zone of disbonding and localisation of the profiles P1-P2, P3-P4, P5-P6 and the specific points Q i added as reference points. The thickness was calculated from the displacement fields of the pair-wise images dating from February 11, 2020 and April 13, 2020. Green lines indicate the measured dispersion. Red dots indicate specific topographic target locations (see Figure 3.2) and red lines show the lateral and top limits of the landslide (see Figure 3.1). The black line shows the forestry track while the brown polygons highlight the apparent bedrock.

Figure 3 .

 3 Figure 3.19: Maximum and mean thickness calculated over the profile P3-P4. Each pair of values are related to a rheological parameter. A value of 1/2 refers to a linear vertical velocity profile, a 2/3 value stands for a Newtonian viscous flow, a value between 2/3 and 1 means a plug flow and value equal to 1 means a rigid sliding block.

Figure 3 .

 3 Figure 3.20: Zoom on the P3-P4 profile at the top of the sliding surface. The thickness is calculated according to several beta values. Higher β values (superior to 1) are mathematically possible but they do not have physical meaning.

  Figure3.21 shows the geotechnical profiles drawn in 2018 corresponding to the moving zone of 2020 detected by time lapse photography. As a result, at least five sliding surfaces are involved in the Montgombert slide. The thickness of the sliding surfaces at points Q1, Q2 is measured, from the geotechnical profiles, at 3 m and 10 m respectively. The associated thickness calculated by photography time lapse, if we consider a movement that is similar to a rigid block (α = 1), are respectively, 1.3 m and 2.4 m. The values for 2018 and 2020 are more or less of the same order of magnitude, but there remains one notable difference. This can be explained by several interpretations. The most likely of which is to assume that the first slide lobe has been retrogressing since 2018. Indeed several

Figure 3 .

 3 Figure 3.21: Geotechnical profile carried out by Mathy and Lorier (2018) as part of the monitoring of the Montgbombert landslide. The continuous line represents the topographic surface of the natural terrain in 2016. The dotted lines represent slide lobes. Points Q1, Q2 and Q3 are markers whose location is given in Figure 3.18. Incl. means inclinometer and Pz. is for piezometer. The blue dotted line indicates a potential water level revealed by the piezometric measurements.

Figure 3 .

 3 Figure 3.22: Schematic representation of the various parameters included in the geometrical relationships. (a) Schematic representation from Meier et al. (2020) for the computation of the volume by the half ellispoid and the elliptical paraboloid methods. W is the width of a half ellipsoid along the slope. L rh and W rh are respectively the horizontal length and width of the shape along the horizontal. A rh is the associated surface area. z max is the maximum depth along the vertical. (b) Schematic representation from Nikolaeva et al. (2014). A, B and C are the ellipsoid major axis, respectively, in the X, Y and Z directions. a is the major axis of the sliding plane in the X direction. h is the distance between the plane of the landslide surface.

  3.27)Nikolaeva et al. (2014) propose a new, more complex configuration by not considering a purely rotational movement. They add a translational component to the geometry of an ellipsoid. This means that movement consists of a rotational and a translational component. This new component is an ellispoidal segment constrained by two parallel planes which are the surface of the slide and the plane of detachment. Let us consider V t the volume associated to this new configuration. It is defined as:

4. 3 . 3

 33 Volume estimation of the Montgombert landslideTo estimate the volume of the Montgombert slip of 2020 we applied both the geometrical relationship and the empirical relationship (Figure3.23).To calculate the volume we considered the dimensions of the unstable part highlighted by image analysis: L h = 64 m, W h = 72 m. The maximum vertical surface depth was averaged to z max = 7.2 m. Hence, the volume estimated by the half ellipsoid method equal to 17261 m 3 and the one estimated by the elliptic paraboloid equal to 12946 m 3 .

Figure 3 .

 3 Figure 3.23: Estimated volume of the moving mass during the acceleration 2020 at the Montgombert landslide based on previous published empirical relationships.

  define three modes of toppling: (1) a flexural toppling which is characterized by a continuous break in flexure of the columns. It is a retrogressive movement which involves predominantly schists and foliated (phyllite) rocks. (2) a block toppling which involves individual columns divided by spaced joints.

  such as the landslides of the Panissats Dessous in 2014 and 2016, regular landslides at the Nant Cortay localization (≤ 10000 m 3 ) and the Montgombert landslide in 2016 and 2018 (Figure 4.1). Jeannin (2001) counted in 2001 a total of 111 landslides from 1 to more than 10000 m 3 (overhanging micaceous shales) in the high gorges; in the middle gorges 19 landslides of more than 1 m 3 and in the lower gorges 23 landslides. Today, the road is still closed because of the last landslide event which occurred in February 2019 involving 10000 m 3 . This landslide called the Cliets landslide is situated on the left bank of the Arly river. Generally the presence of steep rock slopes and of dense forest cover make the valley difficult to monitor and gravitational instabilities still represent the major issue.

Figure 4 . 1 :

 41 Figure 4.1: Geology context of the Arly valley with (a) the localization of the Cliets rock slide and its geological context simplified from (Baudin, 2015) and (b) a geological cross-section of the upper Arly gorges (which localization is indicated in (a)) modified from (Dussauge-Peisser et al., 2002; Jeannin, 2001) (source http://infoterre.brgm.fr/).

Figure 4 . 2 :

 42 Figure 4.2: Mapping of the antecedent rockfalls at the Cliets unstable slope. (a) Rockfall sources are shown on the topographic map of July 2018. (b) Photographs of the July 2013 and the December 2018 rockfalls. ml stands for linear meter.

1. 3 . 1

 31 Chronology of the site's instrumentationAs shown in Figure1.2.2, several areas of the Cliets site have been destabilized over time. It was therefore necessary to adjust the techniques according to the spatial extent of the movement and its dynamics. This is why the Cliets site has gradually been instrumented. Table4.1 summarizes the instruments put in place from 1996 to 2019. The instruments generally used in an operational context are the topographic targets whose position is surveyed on a regular basis (annually, monthly or hourly). This survey can be done manually if the frequency of measurement is low enough (annual survey). If the frequency increases (hourly survey), an automatic theodolite monitoring is preferable. Theodolite instrumentation is often accompanied by extensometric measurements to monitor the opening or closing of rock fracture. They are often connected to road closure lights to protect lives. At the same time, a weather station close to the rockslide was deployed and geological observations as well as stereo-interpretation were realized in the course of 2003

  Figure 4.3. Because of an acceleration of benchmark displacements in 2013, a new topographical network composed of 10 benchmarks was setup on 3 December 2013 and displacements were measured hourly by an automatic total station till the 14 January 2014 (period 4). The total station is located in the right bank of the Arly river in front of the Cliets tunnel. On 28 July 2018, another new topographic network of 40 benchmarks was installed (Figure4.3), and displacements are measured hourly by an automatic total station till June 2019 (period 6). To be consistent with the rest of data (1.3.4), the displacement and velocity are expressed in a daily time step. The associated measurement error is calculated on reference targets that are stable over time. The error is of the order of few millimeters for both periods.1.3.4 Forcing factorsThe daily gross rainfall observed by three Meteo-France meteorological stations (Col des Saisies, Albertville, Cliets, Figure4.4a)) over the period 2004 -2019 is used to define the relation between rainfall

Figure 4 . 3 :

 43 Figure 4.3: Topographic network set up to monitor the landslide activity. (Left) Network associated to the periods 3 and 4 on the photography from December 2013. (Right) Network associated to the periods 5 (approximated positions) and 6 on the photography from July 2018. P3, P2 and P1 represent seismic profiles realized in 2018 (cf. section 1.3.2).

Figure 4

 4 Figure 4.4: Forcing factors related to the Cliets rockslide. (a) Weather station locations. (b) Location of recorded earthquakes around the Cliets site. The size of the seismic event indicators (circles) corresponds to the magnitude and their colors show the year in which they were recorded.

Figure 4

 4 Figure 4.5: Structural analysis of the Cliets rockslide. (a) 3D view of the unstable slope and of the discontinuities for the two last failure events. (b) Stereonet of the major discontinuities: S0 and F1 are observed in the main scarp area; F2 and Dn are only observed during the period 3 (2004-2014) and D1 is observed only in 2018. Colored points are the associated poles. (c) and (d) Topographic profiles along the slope.

Figure 4 . 6 :

 46 Figure 4.6: Effective daily rainfall (mm) reconstructed according to section 1.3.4, compared to the cumulative 3D displacements (m) for both periods 3 (a) and 6 (b). The numbers refer to the targets located in Figure 4.3. The circles correspond to gravitational events detected thanks to the knee point method. An example is given for period 6 in the inset of (b), over the photography of the landslide in 2018. The points indicate topographic target positions. The colors correspond to the time of detection.

Figure 4

 4 Figure 4.7: Dependency analysis between the daily effective rainfall, reconstructed according to section 1.3.4, and the 3D surface velocity. Correlograms are calculated over a rolling window for the two periods 3 (a) and 6 (b). A vertical dashed line indicates a gravitational event, as defined in section 1.4.2.

Figure 4

 4 Figure 4.8: Inverse velocity and time to failure forecast. Forecasts for (a) benchmarks 6 and 7 for the period 3; (b) benchmarks 12, 29, 36 for the period 6. The red arrows indicate the estimated time to failure computed within the time span represented by a grey rectangle. The black arrow shows the real time of failure.

Figure 4

 4 Figure 4.9: Clustering of the slope into kinematic units based on the behavior of the benchmark velocity properties.

Figure 4 .

 4 Figure 4.10: Plunge angle versus time calculated for a series of benchmarks for the period 3 (2004-2014) and for the kinematic units for the period 6 (2018-2019). The bottom figure corresponds to a zoom of the bold rectangle highlighted after 2018.

Figure 4 .

 4 Figure 4.11: Rate and state models applied to benchmarks 6 and 7 of the period 3.

Figure 4 .

 4 Figure4.12: Rate and state models applied to clusters 1, 2 and 3 for the period 6.

  ) an estimate of the time to failure by including more measurements (Figure4.13 a), (2) an estimate of the time to failure by taking into account a rolling window of a fixed number of measurements (Figure4.13 b). In Figure4.13, tc stands for the time to failure and t * to the time of the last date used for the 1/v model. For both cases, the modelled time to failure converges towards a single value while the estimated time to failure of the inverse velocity applied directly to the data deviates by about twenty days for the first case and by about ten days for the second case. Including more data in the estimation of the time to failure leads thus to more dispersion in the calculation than considering a fixed number of data. But, by considering the second case, the prediction of time to failure tc, converges to a plateau ten days before the failure.We apply the rate and state model on different sections of the velocity curve to estimate the evolution phases of the slope(Figure 4.14).

Figure 4 .

 4 Figure 4.13: Estimated time to failure tc versus the time t * of the last point used for the fit. The days are indexed relative to the first data point in the series. The green dashed line indicates the real time to failure. (a) Estimation of the time to failure by including more measurements and (b) Estimation of the time to failure by taking into account a rolling window of a fixed number of measurements.

Figure 4 .

 4 Figure4.14: Application of the rate and state models before the main rock slide events on benchmark 7 of the period 3 and benchmark 1 of the period 6.

Figure 4 .

 4 Figure 4.15: Images from the image time series of the Cliets rockslide. The images are acquired with a Canon 2000 D with a focal length of 24 mm. Snow cover can be thick on the slope but also in front of the camera in the foreground.

  Figure 4.18). Since the illumination of the images changes over time, we chose to apply a variable threshold, established for each image. The presence of snow is reflected in the histogram of the normalized blue band by the presence of two local maximums (bi-modal distribution). The area of the curve closest to the value 255 corresponds to white pixels and therefore to snow (Figure 4.18a)). A mask has been applied to hide the lower part of the images since snow has accumulated over time in the foreground. The calculation of the number of pixels was therefore estimated only on the upper part of the image. 81 images were manipulated in the period from January 17, 2019 to February 14, 2019.

Figure 4 .

 4 Figure 4.16: Cumulative 3 D displacements in meter computed from the time series processed in a) and b) CMC mode and c) and d) in V SC mode. Displacements were averaged in a 10x10 window size centered on the target. The RMSE calulated on the magnitude of the 3 D displacement is represented on a photography taken from the fixed camera. Colors correspond to the value of the RMSE and the circle size indicates visually the pixel recurrence. The larger the diameter of the circle, the greater the pixel recurrence.

Figure 4 .

 4 Figure4.17: Pixel recurrence in time series computed on the filtered correlation results as presented in Chapter 2. a) Recurrence calculated on the time series processed according to the CMC approach. The first image of the time series was used as master image. b) Recurrence calculated on the one processed according to the V SC method. Number of information is more important for the first case a) in the central area as well as on the edges. The information acquired by the V SC method is all the more filtered as the residual movement (on the edges) is similar to noise.

Figure 4 .

 4 Figure 4.18: Adapted snow-to-snow method. a) Histogram of the normalized blue band of a photograph.Green point shows the snow threshold calculated from the polynomial fit. This is the equivalent of finding the minimum local of the polynomial fit. Some pixels may be considered as snow pixel while they are not. Indeed, high values can also indicate pixels with high sun illumination. Therefore, this method may overestimate the number of pixels.

Figure 4 .

 4 Figure 4.19: Snow covery calculated over the image time series and expresses in %.

  Figure 4.19 shows the evolution of the snow cover over time. It is represented with the effective rainfall expressed in millimeter and the cumulative 3D displacements of 4 targets located in Figure 4.16e) and d). Several snowfall and snow melt events are highlighted by image analysis. Each snow episode is accompanied by a peak of rain.

Figure 4 .

 4 Figure 4.21: Mean displacement rates in m.day -1 calculated over the post event image time series. Images were processed in the V SC mode. The velocity is averaged over the rocky spur area in order to ensure the new stability of this zone.

Figure 5 . 1 :

 51 Figure 5.1: Typical phases of landslide activity, periods of usability of terrestrial optical image time series and the corresponding purpose of the surveillance (defined according to(Mazzanti and Pezzetti, 2013)). The usability is defined as follows: "+++" means that optical terrestrial images may provide relevant results and are useful for landslide understanding. "++" means that optical terrestrial images are useful but some requirements have to be considered. "+" means that the data and technique could be applicable but only on specific cases. Possible monitoring alternatives should be considered.

Figure A. 1 :

 1 Figure A.1: Effet de la distorsion sur une grille régulière (Noury et al., 2017). a) Distorsion radiale posiive (à gauche) et négative (à droite). b) Distorsion tangentielle.

  cos(φ )sin(κ) cos(ω) + sin(φ ) cos(κ) sin(ω) sin(κ) sin(ω) + sin(κ) cos(φ ) cos(ω) sin(κ) cos(φ ) cos(κ)cos(ω) + sin(κ) sin(φ ) sin(ω)cos(κ) sin(ω) + sin(κ) sin(φ ) cos(ω)sin(φ ) cos(φ ) sin(ω) cos(φ ) cos(ω) post-event monitoring analysis partir des erreurs de mesure propres à chaque paramètre pour être ensuite appliqué de manière aléatoire sur la donnée. L'écart type associé au bruit est renseigné dans la table B.1. L'erreur est calculée de la manière suivante pour les trois composantes : EW notée X, NS notée Y et la vertical notée Z : ∆Err = Resultat bruite -Resultat re f erence σ gauss (B.1)

  également valable pour les erreurs en Y et en Z notées respectivement ∆Err y et ∆Err z .Il en résulte que les matrices d'erreur par composante égalent (valeurs métriques) :Nous pouvons constater que les erreurs liées au MNT bruité sont de l'ordre de plusieurs centimètres tandis que celles liées aux positions sont légèrement plus faibles. Les erreurs les plus élevées sont celles qui correspondent au bruit introduit dans le picking et donc dans la calibration de la caméra et il est intéressant de voir qu'une direction est plus impactée que l'autre. Dans le cas des images des Cliets (acquises en mode portrait), la direction u correspond à la direction perpendiculaire au mouvement et la direction v à la direction du mouvement. Un erreur de quelques pixels dans la direction perpendiculaire au mouvement a donc plus d'impact qu'une erreur introduit dans la direction du mouvement.Pour quantifier de manière plus précise ces erreurs, deux approches seraient à envisager: (1) une approche statistique dans laquelle une grande quantité de résultats bruités sont générés. On obtiendrait alors une distribution dense d'erreurs. (2) une approche plus direct dans laquelle le bruit introduit est connu.

Figure C. 1 :

 1 Figure C.1: Vitesse surfacique exprimée m.day -1 et calculée le long de la série temporelle d'images acquises après la crise de Février (Glissement des Cliets, Savoie, Fr.). Les images onté été acquises entre le 15 Février 2019 et le 4 Octobre 2019. Elles ont été traitées par TSM et corrélées selon l'approche V SC. Les flèches noires indiquent la direction du mouvement.

  Indeed, high values can also indicate pixels with high sun illumination. Therefore, this method may overestimate the number of pixels.150 4.19 Snow covery calculated over the image time series and expresses in %. . . . . . . . . . 151 4.20 Displacement rates in pixel.h -1 calculated from the correlation technique. Images were taken on February 9, 2019 on the day of the slope failure. . . . . . . . . . . . . . . . . . 153 4.21 Mean displacement rates in m.day -1 calculated over the post event image time series. Images were processed in the V SC mode. The velocity is averaged over the rocky spur area in order to ensure the new stability of this zone. . . . . . . . . . . . . . . . . . . . . 154 5.1 Typical phases of landslide activity, periods of usability of terrestrial optical image time

Table 1

 1

		.1: Summary of landslide monitoring methods modified after Malet et al. (2002); Delacourt et al. (2007); Michoud et al.
	(2010).
		Monitoring methods	Measurement type	Precision	Some References
		Global Positioning System	Target displacements	1-2 mm	Malet et al. (2002); Squarzoni et al. (2005); Gili et al. (2000)
	In-situ	Tacheometry, theodolite	Displacements along ver-tical	20 mm	Angeli et al. (2000)
		Precise levelling	Target displacement	0.5 mm	Malet et al. (2002)
		Extensometry	Local displacement	0.5 mm	Angeli et al. (2000); Corominas et al. (2000)
		Passive sensors -Optical Imaging
		Satellite imaging	Displacement	0.6 -80 m	Kääb (2002); Delacourt et al. (2004); Stumpf et al. (2018)
	Remote sensing	Aerial imaging Terrestrial imaging Passive sensors -Thermal imaging	Displacement Displacement Temperature map	0.5 -2 m cm -dm 1 -2 • C	Niethammer et al. (2012); De-witte et al. (2008) Travelletti et al. (2012); Stumpf et al. (2015) Guerin et al. (2019)
		Active sensors -Laser scanning
		Aerial laser scanning	DEM comparison	dm -m	Jaboyedoff et al. (2012)
		Terrestrial laser scanning	DEM comparison	mm -cm	Kromer et al. (2017); Noël et al. (2018)
		Active sensors -Radar interferometry
		Satellite	LOS displacement	2 mm	Strozzi et al. (2004); Nikolaeva et al. (2014)
				1.7x10 -2
		Aerial	LOS displacement	displ.	Hu et al. (2020a)
				gradient
		Ground-based	LOS displacement	∼ mm	Antonello et al. (2004); Bozzano et al. (2011)

Table 1 .

 1 2: Summary of the study sites characteristics presented in these works. Location corresponds to the French departments The average velocity refers to the classification proposed byCruden and Varnes (1996).

	Landslide name	Chambon	Pas de l'Ours	Les Cliets	Montgombert
	Location	Isère	Hautes-Alpes	Savoie	Savoie
	Process	Rotational landslide	Complex landslide	Toppling rockslide	Complex landslide
	Lithology	Sedimentary rocks from the Lias	Lustrous clay-shales of lower dated moraine formations Cretaceous and unconsoli-	Micachists and shale screes	Micachists
	Average velocity dur-ing the acquisition	slow motion	moderate motion	moderate motion	slow motion
	Dimension of the area
	of interest (length x	170 m x 92 m	500 m x 500 m	72 m x 80 m	150 m x 115 m
	width)
	Complementary instru-mentation	Terrestrial LiDAR (non-permanent)	Total station, Terrestrial Li-DAR (permanent), GBInSAR	Total station, Extensometer	Total station, Terrestrial Li-DAR (non-permanent), Ex-tensometer
	Context of study	Post event monitoring	Pre/Post-failure monitoring	Pre/Post-failure monitoring	Pre/Post-failure monitoring

  . (2014) demonstrated that a fixed installation of a camera allows us to obtain time series of the displacement fields from multiple pairwise image correlation. The obtained displacement fields are able to

. Recently, Gance et al. (2014), Travelletti et al. (2012), Roncella et al. (2014) and Mishra et al

Table 2 .

 2

	Use case	Image modality	Image matching parameters [pix]	Size of the stable area [%]	Correlation strategy
			Inc	SzW	Reg	C min
	Chambon landslide *G	20x20 4x4	0.2	0.5	36	VSC
	Pas de l'Ours land-slide	*G	24x24 4x4	0.2	0.5	9	VSC

1: Combination of parameters used to process the two time series for respectively the Chambon and the Pas de l'Ours landslides.

Table 3 . 1 :

 31 Constants from the surface area -volume power law established by several authors afterGuzzetti et al. (2009) andMichel et al. (2020).

	η	γ	Range of values for the surface area [m 2 ]	Number of landslide	Source
	0.1479 1.368 2.3x10 0 -1.9x10 5	207	Simonett (1967)
	1.826	0.898 5x10 1 -1.6x10 4	1019	Larsen and Torres-Sanchez (1998)
	0.155	1.091 7.0x10 2 -1.2x10 5	124	Guthrie and Evans (2004)
	0.0844 1.4324 2.0x10 0 -1.0x10 9	539	Guzzetti et al. (2008)
	0.074	1.145 2.0x10 0 -1.0x10 9	677	Guzzetti et al. (2009)
	0.05	1.5	4.0x10 2 -1.0x10 5	4984	Hovius et al. (1997)

  1.2.1 Context of the Arly ValleyThe Arly Valley is located in the Savoie department in the Central French Alps (Figure4.1). The valley is an important infrastructure connecting to many ski resorts. About 4650 vehicles per day on average are using this road (data set in 2001

Table 4 .

 4 1: Chronology of the site's instrumentation

	Period	Label Instrumentation	Averaged velocity (from topographic measurements) Source
				Vertical	Planimetric
	1996 -2003	1	5 geodetic references + 9 extensome-ters for monitoring the transverse opening of rock corridors	200 mm.year -1 (upper part)	100 mm.year -1 (North part)	Kasperski (2008)
				20 mm.year -1 (middle part)	30-40 mm.year -1 (South part)
			10 geodesic references to complete
	2000 -2003	2	the previous network (annually mea-
			sured)
	Dec. 2003 and Feb. 2004	Purges		Loss of instruments
	2004 -2014	3	9 topographic targets (monthly mea-sured)	35 mm.year -1	47 mm.year -1	Mathy and Lorier (2013)
			10 additional topographic targets
	2013 -2014	4	measured by total station + exten-
			someters
	June 2014		Purges		Loss of instruments
	2014 -2018	5	10 topographic targets measured by total station (measured every 3hours) + extensometers	1.8 mm.month -1 (during 2014)	2.6 mm.month -1 (during 2014)	Mathy and Lorier (2014)
	2018 -2019	6	30 topographic targets measured by total station + extensometers + fixed permanent camera	2.8 mm.day -1 (rocky spur)	6.2 mmday -1 (rocky spur)	Mathy and Lorier (2018)

Table 4 .

 4

	Characteristics	Col des Saisies Albertville Cliets
	Observation time sampling	daily	daily	hourly
	Horizontal distance to unstable slope (km) 3.9	13	0.8
	Elevation difference to unstable slope (m) 904	377	0
	Weight (2004-2018)	0.77	0.23	/
	Weight (2018-2019)	0.15	0.23	0.8
	typically between 0.35 and 0.60 cm.day -1 .

2: Characteristics of the meteorological stations used in the Cliets rockslide analysis: time sampling, distance from the slip and weighting coefficients. • C -1

Table 4 .

 4 3: Results of the friction model realized on two simulations for the period 3 and three simulations for the period 6. They are expressed by the m and Xi estimated parameters and the NRMSE associated.

	Time span Data set	Xi	m	NRMSE
	Period 3	benchmark 6 0.99	4.99	1.2
	Period 3	benchmark 7 0.95	4.99	1.4
	Period 6	cluster 1	0.0102 1.0001 3.1
	Period 6	cluster 2	0.0017 1.001	4.5
	Period 6	cluster 3	0.0076 1.004	3.5

Table 5 .

 5 1: Summary of the products that can be generated from fixed cameras respect to the purposes of the study and some specific requirements.

	Landslide infor-mation type	Purposes	Products	Requirements
		Long-term	Geomorphological change de-	Monoscopic or multi-view im-
		monitoring	tection maps	ages.
	Geomorphological	Slope character-	Environmental condition maps	Textured images with homoge-
	information	ization	(e.g. for instance snow cover)	neous radiometric properties.
				Presence of object with rigid type
				deformation.
				Above requirements.
	Kinematic	Long-term monitoring	Surface displacement / velocity fields and time series	Calibrated cameras.
	information			Precise	time	series	co-
				registration.
		Long-term monitoring	Deformation / Strain fields maps	Above requirements.
		Slope character-		Dense spatial coverage of infor-
		ization		mation.
		Long-term monitoring	Identification of landslide me-chanical regime	Above requirements.
	Mechanical information	Early-warning		Dense image time series close to the main acceleration / failure event.
		Long-term	Landslide thickness / depth	Stereo or multi-view datasets+ all
		monitoring	maps	above requirements.
		Slope character-
		ization
		Early-warning	Estimation of time to failure	Mono or multi-view datasets + all above requirements.
	Forecasting /
	modelling information	Long-term monitoring	Relationships velocity-rainfall	Mono or multi-view datasets + above requirements.
		Slope character-
		ization

The CIFRE/ANRT research contract -Conventions Industrielles de Formation par la REcherche -is an agreement which enables a French company to develop a research collaboration with a public research laboratory through a Private-Public Partnership (PPP).

MicMac: Multi Images Correspondances par Méthodes Automatiques de Corrélation. It is a free open-source software, implemented at IGN (Institut National de l'Information Géographique et Forestière, France).

MPIC: Multiple Pairwise Image Correlation. It is a processing strategy for time series of stereoscopic very-high resolution (Pléiades) images(Stumpf et al., 2017), and time series of Sentinel-2 and Landsat-8 images(Stumpf et al., 2018).

COREGIS: Python implementation to improve the co-registration of Sentinel-2 and Landsat-8 images(Stumpf et al., 2018).

COSI Corr: non-commercial software developed for measuring ground deformation using optical satellite and aerial images. Correlations are performed in the frequency domain.

MEDICIS: Moyen d'Evaluation de Décalages entre Images, Commun à l'Imagerie Spatiale. This correlator is developed by CNES (Centre National d'Etudes Spatiales, France). Correlations are process either in the spatial or in the frequency domain.

MicMac: Multi Images Correspondances par Méthodes Automatiques de Corrélation. It is a free open-source software, implemented at IGN (Institut National de l'Information Géographique et Forestière, France). Correlations are processed in the spatial domain.

DPIVSoft: software for image correlation. It based on the Particle Image Velocity (PIV) algorithms(Keane and Adrian, 1992; Aryal et al., 2012)

Remerciements

This work is presented around the analysis of four landslides located in the French Alps

Chapter 2

Image time series analysis for a monoscopic model

This chapter introduces a generic pipeline proposed for the processing of long image time series acquired by Single Lens Reflex (SLR) cameras. The system associates modules for (1) the selection of the image sequences, (2) the registration of the image stacks and the correction of the camera movements, and (3) the calculation of the terrain motion using change detection approaches.

The system is based on the open-source photogrammetric library MicMac and tailored for the processing of monoscopic images. A sensitivity analysis is conducted to design and test the image processing for two use cases; the Chambon landslide (Isère, France) characterized by slow motion (< 10 cm.day -1 ), and the Pas de l'Ours landslide (Hautes-Alpes, France) characterized by moderate motion (> 50 cm.day -1 ). Four categories of parameters are tested: the image modality, the image matching parameters, the size of the stable area used in the co-registration stage, and the strategy used to combine the images in the time series. The application of the pipeline on the two use cases provides information about the kinematics and the spatial behavior of the landslides. This chapter is based on the article: Desrues M, Malet JP, Brenguier O, Point J, Stumpf A, Lorier L (2019b) TSM-Tracing Surface Motion: A Generic Toolbox for Analyzing Ground-Based Image Time Series of Slope Deformation. Remote Sensing 11(19):218 9 (https://doi.org/10.3390/rs11192189).

Comments: Some additional analyses are inserted in this chapter to keep some continuity in the reading of the manuscript.

Chapter 4

Pre-and post-event monitoring analysis according to the velocity and the cumulative displacement:

+ Bln θ θ 0 (4.4)

where A and B are constants and θ is the frictional state variable which evolves with time according to Equation (4.5):

where Dc is a representative slip distance. The resolution of Equation (4.4) is proposed by Helmstetter et al. (2004) and summarized in Sornette et al. (2004) by introducing reduced variables. Equation (4.4) is then written as:

where X, D and m are the reduced variables expressed with:

where τ and σ are, respectively, the average shear and the normal stresses. Equation (4.6) is transformed to Equation (4.8):

where t = t/T is a dimensionless time which means that time is usually expressed in the units of T = Dc D (Helmstetter et al., 2004;Sornette et al., 2004). Helmstetter et al. (2004) indicate that four possible regimes can be identified and are defined by the parameters Xi and m: (1) m > 1 and X i < 1 indicate an acceleration which leads to a finite-time singularity of the velocity expressed by Equation (3) with α = 2; (2) m < 1 and X i > 1 indicate an acceleration evolving towards a stable regime; the two other combinations (3) m < 1 and X i > 1 and (4) m > 1 and X i > 1 indicate a deceleration of the motion.

The velocity model is then determined by inverting the parameters m, D, T and the initial condition of the parameter X i . The cost function used is the average Normalized Root Mean Square Error defined as:

where δ is the average velocity. As the inversion is very poorly constrained, several iterations (more than 50 iterations) are used and the best fit corresponds to the lower value of the NRMSE.

The frictional model is applied to benchmarks 6 and 7 for the period 3 (Figure 4.11) and to clusters 1, 2 and 3 for the period 6 (4.12). Values of the m and Xi parameters are shown in Table 4.3. A high value of m is obtained for the two benchmarks of the period 3 while according to Helmstetter et al. (2004) Chapter 4

Pre-and post-event monitoring analysis the slope collapses and the de-correlation of several zones within the stable part is potentially linked to the movement of the wire mesh placed on top of it that was carried away by the collapse. This hypothesis can be verified by the absence of movement detected between 2 and 3 pm. On the other hand, a strong signal is measured on the fracture plane, which indicates a surface movement due to erosion.

Landslide dynamics in a post-event context

After the main event, a post-event analysis is performed to determine residual movements. The slope may present stable areas as well as moving areas that may be destabilized by the departure of masses. In the case of landslides that affect roads or river beds for example, human intervention is required to clear the paths and/or carry out entropic purges. But before any human intervention on site, it is important to assess the stability of the slope.

In the case of the Cliets, purges and reconstruction work on the tunnel were carried out a few months after the event. A series of 15 images (acquired between February 15 and October 4, 2019) were processed to ensure that there were no possible departures of new unstable masses. The time interval between each image varies over time: one image every two days was used first, then every two weeks. A DEM from a UAV acquisition realized after the event was used in the process. Since little signal is detected outside the collapsed area (see Annex) we opted to represent here only the average speed of the area which correspond to the rocky spur. Figure 4.21 represents the averaged velocity in m.day -1 of this area.

From February to March 2019 the average velocity decreases from a few centimeters to a few millimeters per day. In April/May it reaches a plateau where the average velocity ranges between 2 and 4 mm.day -1 .

The behavior is similar to that of a linear state in which conditions are stable, speeds are constant and acceleration is zero (Scoppettuolo et al., 2020).

Chapter 4

Pre-and post-event monitoring analysis

General conclusions

Time-lapse photography is a low-cost remote-sensing method allowing the acquisition of high resolution and high value information but which is still under-exploited in the field of natural hazards. In this work, we developed several processing tools for the analysis of time series of terrestrial optical images acquired from fixed cameras. The tools allow the analysis of RGB images acquired from mono (or stereo-) single lens cameras in order to extract primary (geomorphology, kinematics) and secondary (geometry, mechanics) information relevant for landslide monitoring and hazard assessment. To this end, we have developed and implemented the TSM modular toolbox, built on previous works of CNRS/EOST (Travelletti et al., 2012;Stumpf et al., 2017). The tool was developed to be as versatile as possible by taking into account data from several acquisition systems and for imaging several types of motion, with also the possibility to be used both for science and operational applications. The toolbox has been designed to be modular and upgradeable.

This work is a contribution to the development of new remote sensing methods for monitoring and understanding gravitational movements from terrestrial photographs (Chapter 1). The accessibility of the method and the low cost of the instrumentation make this method attractive compared to other types of similar instruments such as Laser Scanning, ground-based InSAR and satellite imagery. The processing of image time series acquired from fixed cameras irequires planning since instrumentation, site configuration and processing parameters must be taken into account to measure displacement fields accurately both in amplitude and direction.

Chapter 2 proposes a processing chain which is intended to be as exhaustive as possible for the processing of images acquired in monoscopic mode. It is broken into several modules such as the camera calibration, the image selection, the correction of camera movement and finally the calculation of displacement fields.

The modules were intensively tested on two use cases, namely the Pas de l'Ours landslide (Hautes-Alpes, France) characterized by moderate displacement rates (ca. 50 cm.day -1 ) and the Chambon landslide (Isère, France) characterized by low displacement rates (< 10 cm.day -1 ). The tests indicate that relevant displacement above 1 cm.day -1 can be detected and that the errors relative to reference in-situ measurements is of the order of ca. 4 cm for the Chambon landslide and ca. 15 cm for the Pas de l'Ours landslide.

The influences of internal and external parameters on the displacement results were also discussed (see paragraph below: On the influence of the parameterization of the processing).

Chapter 3 implements a module able to process stereoscopic acquisition in order to allow the dense reconstruction of 3D displacement fields directly from the results of image correlation. Such mode of acquisition is suitable for the analysis of complex movements for which displacements along the normal are not negligible. The errors calculated in reference to in-situ data is less than 4 cm for the three EW, NS and Up components. We further calculate the tangential and normal displacements in order to estimate the strain fields and invert the 3D displacement fields for estimating the thickness of the moving mass. We

Chapter 5

Pre-and post-event monitoring analysis highlights: (1) the weak influence of the camera position in relation to the object for monoscopic acquisitions, and (2) the increased accuracy of the measurements for stereoscopic acquisitions.

2. On the influence of the processing parameters

In our approach, the displacement fields were calculated using NCC-based and hierarchical image correlation methods. The choice of the correlation parameters has a direct influence on the quality (accuracy, precision, redundancy of information, spatial coverage) of the results. TSM is based on the MicMac open source image correlation library for which the parameters that show the highest sensitivity are the size of the correlation window and the size of the search window (Chapter 2).

In order to increase the signal-to-noise ratio to extract a maximum of reliable information, several correlation strategies have been implemented and tested (Chapter 2-4). These strategies are the (3) the quantity of information varies according to the coherence between two correlated images and its quality according to the detected noise. For long-term monitoring, these strategies can be applied seasonally in order to highlight transitory phenomena characterized by a modification of the equilibrium of the slope by triggering factors. This would allow a better characterization and understanding of the long-term behavior of the slope.

3. On the primary and secondary products that can be generated TSM allows computing several outputs in order to analyze landslide motion both qualitatively and quantitatively. Primary and secondary products can be calculated, both from the monoscopic and the stereoscopic approach:

(1) primary products are either qualitative (such as maps of affected areas) and quantitative (such as surface displacement fields; Chapter 2, Chapter 4), (2) secondary products are results of postprocessing analyses such as spatio-temporal clustering (amplitude, direction) of the displacement fields in order to delineate zones of homogeneous behaviour and mechanical/predictive assessments such as time of failure models or dynamic evolution models (Chapter 4). With a stereoscopic approach, additional products are obtained since the tangential and normal surface fields can be calculated in addition to the surface displacement fields (i.e., a monoscopic approach assumes that the movements are mostly in the plane tangential to the surface). It is therefore possible to generate strain maps and to infer the geometry of the moving mass from the calculation of depth maps (Chapter 3).

4. On the automation of the processing Automating the processing workflow means an easier use of the algorithm but also to have standard solutions. This standardization can be applied to the input/output file formats, to the selection of images, to the correlation parameters (correlator and correlation strategy). In practice, the TSM toolbox has been built to be used both 1) as a moderately tunable automated processing service (such as for instance the automated image computed displacements exceed a threshold, processing is performed at higher resolution. The calculation time could thus be shortened and the data transfer/storage optimized.

• Finally, as we have demonstrated, previously, that acquiring and processing terrestrial optical images requires a priori information about the object being studied (e.g., dimensions, surface velocity), combining a priori information with various external and internal camera parameters could help in designing a decision support tool for the "best" choice of camera position, focal length and acquisition frequency (Gance et al., 2014).

In terms of image processing:

• Novel image correlation methods are being developed and in particular optical flow methods. Such approaches could be applied to the analysis of fast gravitational processes and/or video type data (Hadhri et al., 2019;Marsy et al., 2020;Qiao et al., 2020) which could improve the detection of behaviour changes as well as the measurement accuracy or the completeness of information for complex movements. Methods allowing smart fusion of the displacement fields calculated using several correlation strategies and/or several correlators would be interesting in order to obtain redundant, dense and reliable information.

• The analysis of image time series can also go further than the generation of displacement fields.

More and more studies are applying data mining/machine learning methods to detect transitory information. These approaches categorize spatial sectors and periods of similar kinematic behaviours either from satellite images (Julea et al., 2010;Méger et al., 2018) or terrestrial photogrammetry (Pericault et al., 2015). Information reduction and clustering methods are applied to both the displacement fields and the correlation coefficients in order to identify possible deformation sources.

Several dimensional reduction methods can be applied (PCA -Principal Component Analysis; ICA -Independent Component Analysis , LLE -Linear Local Embedding).

• Time series of displacements can show gaps induced by missing data as we saw especially in Chapter 4. To obtain a better temporal resolution, new methods could be implemented in TSM.

Such approaches were successfully applied to time series of displacements from the analysis of optical satellite images (Bontemps et al., 2018) by inversion method, and from the analysis of SAR images (Samsonov and d'Oreye, 2017) by simultaneous processing method (Multidimensional Small Baseline Subset (MSBAS)).

• Finally, in order to improve the usability of the toolbox, both a Graphical User Interface (GUI) and the implementation of the code on High-Performance Computers as webservice would be beneficial for the near-real time analysis of the data, and the processing of many instrumental datasets at the same time.

Appendix

This appendix is written in French in order to fulfill the requirements of the University of Strasbourg on the writing of English thesis manuscripts.

L'annexe comporte trois principales sections qui ont été écrites dans le but d'apporter des précisions et des informations complémentaires aux différents chapitres de la thèse : (1) une section portant sur les caractéristiques des caméras (modèle caméra, paramètres internes et externes);

(2) une deuxième section portant sur l'influence des paramètres externes au processus TSM sur les résultats qui sont les déplacements de surface dans un repère géodésique; (3) une dernière section portant sur l'analyse post-événement du glissement des Cliets à partir de la série temporelle d'images optiques terrestres.

I

Appendix A

Caractéristiques des caméras

Cette première annexe présente les caractéristiques des appareils photographiques, appelés par la suite caméras, utilisées dans cette thèse. Elle introduit la notion de modèle de caméra et définit les paramètres externes et internes du capteur.

1 Modèle de caméra Soient (x p , y p , z p ) les coordonnées d'un point P dans le système de coordonnées caméra et (X p , Y p , Z p ) ses coordonnées dans un système géodésique métrique. Les deux systèmes de coordonnées sont reliées par la relation suivante :

où f est la distance focale exprimée en millimètre. Celle-ci définie la distance entre le centre de projection et le plan image.

Paramètres externes et internes des caméras

Les paramètres qui interviennent dans le modèle sténopé sont divisés en deux catégories : les paramètres externes qui définissent la position de la caméra dans l'espace 3D et les paramètres internes qui décrivent les caractéristiques internes de la caméra et qui sont indépendants de sa position.

II Chapter A

Pre-and post-event monitoring analysis

Paramètres intrinsèques

Les paramètres intrinsèques permettent de relier le système de coordonnées dans le plan image à celui de la caméra. C'est pourquoi, ils sont indépendants de la scène mais dépendent du capteur optique utilisé.

Parmis ces paramètres, on retrouve la distance focale f , les densités de pixels en direction des axes u et v, k u , k v en pix.mm -1 ; et le point principal (x 0 , y 0 ) usuellement défini comme le centre de l'image.

Ensemble, ils définissent la matrice caméra appelée K :

Tel que :

Un terme supplémentaire peut être ajouté à la matrice K sous le nom de skew (terme anglais généralement utilisé). Il permet d'indiquer si les axes u et v sont perpendiculires entre eux. Il est introduit de la manière suivante :

De manière quasi systèmatique, ce terme est considéré comme étant nul.

Paramètres de distorsion

Selon l'optique utilisée, l'image peut être sujette à des déformations géométriques appelées distorsions. Il existe deux types de déformations (Figure A.1):

• Les distorsions radiales qui traduisent un défaut de courbure radiale dand les lentilles de la caméra;

• Les distorsions tangentielles qui traduisent un défaut d'alignement entre le capteur et la lentille.

Ce phénomène se traduit mathématiquement par l'ajout d'un terme correctif δ aux coordonées images non corrigées (x, y) :

Dans le cas d'une distrosion radiale, le terme correctif s'écrit tel un polynôme dépendant de la variable r défini comme étant r 2 = x 2 + y 2 :

Détermination des angles d'Euler

Les angles d'Euler interviennent dans l'étape de calibration de la caméra en étant intégrés dans les équations de collinéarité (Chapitre 2) mais également dans le calcul des déplacements 3D en considérant un modèle stéréoscopique (Chapitre 3).

En généralisant la matrice de rotation, l'équation (A.9) s'écrit :

Pre-and post-event monitoring analysis Ainsi les deux combinaisons possibles sont :

Analyse de sensibilité des paramétres externes sur les résultats de TSM Cette deuxième annexe présente l'analyse de sensibilité des paramètres externes sur les résultats de TSM réalisée sur le jeu de donnés du glissement des Cliets. Cette section se réfère au Chapitre 2.

L'analyse de sensibité a été menée en appliquant la méthode OAT (One-A-Time) qui consiste à modifier un paramètre, à la fois, en lui appliquant une certaine perturbation. Dans le cas de notre analyse, nous avons utilisée une pertubration aléatoire calculée à partir du bruit déterminé pour chaque type mesure.

Nous rappelons que les paramètres externes mis en jeu sont:

• le Modèle Numérique de Terrain (MNT) qui intervient dans le calcul des valeurs métriques de déplacements (transformation du repère plan image (déplacements pixels) en repère géodésiques (déplacements métriques)).

• le pointage des cibles de références utilisées dans le calcul de la calibration.

• La position des cibles utilisées dans le calcul de la calibration. Celle-ci est donnée par la position mesurée par la station totale, qui est ensuie convertie en position géodésique. L'erreur associée correspond donc à la variabilité de la mesure calculée sur des cibles de référence supposées fixes.

Pour cela, un couple d'images a été traités par TSM avec les paramètres externes et internes les plus optimums possible afin de générer les valeurs de références. Un bruit gaussien a ensuite été calculé à