TSM—Tracing Surface Motion: A Generic Toolbox for Analyzing Ground-Based Image Time Series of Slope Deformation

: Passive sensors such as multi-spectral (e.g., Single Lens Reﬂex, SLR) cameras are increasingly being used for geohazards monitoring (landslides, cli ﬀ s a ﬀ ected by rock falls, ice glaciers, and volcano ﬂanks) because of their low cost compared to expensive terrestrial laser scanner (TLS) or radar imaging (GB-InSAR) systems. Indeed, due to the large consumer market, sensor resolution and quality (e.g., gain, dynamic range, and geometry) are increasing rapidly. For gravitational processes, such as landslides, recent research has focused on the development and implementation of image correlation techniques to estimate the spatial shift between at least a pair of images by maximizing a cross-correlation function. A generic and fully automated pipeline is proposed for the processing of long image time series acquired for several site conﬁgurations. 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 respectively 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.


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 [1,2]. Among them, multispectral imagery acquired from satellite [3], airborne [4] or terrestrial [5] platforms can provide distributed information on the terrain motion. These techniques are useful in landslide research for terrain motion detection, mapping and monitoring [6][7][8]. Terrestrial optical imagery is however seldom used for monitoring purposes, though high versatility and high spatial resolution of information can be accessed. With the 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) according to [34]. 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).

Methodology: the TSM Processing Toolbox
The methodology to calculate deformation fields (both in the image plane geometry and in the ground geometry) is presented in Figure 1. The system comprises six modules for (1) sensor orientation, (2) image selection, (3) correction of sensor movement, (4) quantification of ground motion using cross-correlation technique, (5) detection of meaningful motion patterns, and (6) geometrical correction of the displacement fields. The input data consists in image time series acquired from a monoscopic camera (JPG compressed format encoded in 8-bits) with focal lengths in the range 24 to 50 mm. The radial lens distortion is not compensated for this ground-based application as it is assumed that the projection errors are negligible with the used focal lengths (24 to 50 mm) because the errors induced by the radial lens distortion are lower than the image picking precision [5,29]; Section 2.1. Several approaches for image combination (e.g., transformation of the RAW images, sequential vs. redundant sequential correlations, forward vs. backward correlation) are integrated in the processing and are described in detail in Section 3. The registration and matching modules are based on the open source photogrammetric library MicMac [35].

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 orientation is modelled with internal parameters (to simulate the camera optics) and external parameters (to determine the transformation between the relative and absolute coordinate systems). As explained previously, internal parameters are not introduced currently in the processing. External parameters are deduced from the collinearity equations (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 of the three Euler rotation angles [36,37]. Those parameters are computed by applying a Direct Linear Transformation (DLT) on the equations [38]. We assume that the position of the central point (u 0 , v 0 ) is the centre of the image [9] and that the effective focal length and the camera position are known. The error measures are defined with the standard deviation and the Root Mean Square Error (RMSExy) defined in Equations (2)-(3): Remote Sens. 2019, 11, 2189 4 of 23 ∆x i = x proj,i − x picking,i ∆y i = y proj,i − y picking,i , where (x proj,i , y proj,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.

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 orientation is modelled with internal parameters (to simulate the camera optics) and external parameters (to determine the transformation between the relative and absolute coordinate systems). As explained previously, internal parameters are not introduced currently in the processing. External parameters are deduced from the collinearity equations (1): where ( , )are the image-based coordinates, ( , , ) are the world 3-D coordinates ( is the East direction, is the North direction and is the elevation), ( , , ) is the position of the camera in the world coordinate system, ( , ) is the effective focal length and ( , , … ) is

Module 2: Image Selection
Analyzing long image time series (> 30 images) of terrain motion poses the problem of creating sub-sequences 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 selection is, by default, applied for the full image size, but options for masking irrelevant areas (sky, vegetation, water bodies) can be applied.

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 [5], the feature detection method [39] and the SIFT method [40]. 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 [33]. 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). The systematic offsets ( , ) are modelled considering two planes and . As the movement of the camera results mostly from a translation and a rotation, the two planes are modelled with an affine transformation (Equation 4): 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 (4)): 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) [33]. The corrected grids are obtained by subtracting the modelled planes to the disparity grids with Equation (5): where ∆X CORRECTED and ∆Y CORRECTED 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 Image Matching Techniques
The module aims at quantifying changes between two acquisitions through the calculation of disparity vectors using cross-correlation [41]. 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 [42][43][44].
The cost function to minimize is defined by Equation. (6): 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 [43]. 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 [24] in order to highlight meaningful motion patterns. Filtering is applied jointly on the correlation coefficient C (see Section 5.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) is based on the method proposed by [9] (Figure 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) [45], (2) performing a back-projection of the DSM in the image plane according to Equation (1) [46], (3) interpolating the projected DSM coordinates in the image plane, and (4) re-projecting the displacement fields in the geographic coordinate system.

Combination Strategies for Processing Large Image Datasets
Several combinations of images are possible to process large time series [32,47]. TSM implements three approaches ( Figure 3) to be flexible to any type of use cases. The results are compared and discussed in Section 5.4: (1) the Common Master Correlation method (CMC) consists in forming all pairs of images from a common master image [31] (Figure 3a); (2) the Variable Sequential Correlation method (VSC) consists in correlating the images sequentially along the time series (Figure 3b). The correlations are thus performed for each interval [I, I + n] 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 3c). The average displacements are then calculated. This approach allows obtaining redundant information [48] 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 [33] for processing time series of satellite imagery. Correlations can be performed for each interval [I, I + n] in forward and backward mode ( Figure 3d).

Application to Use Cases: the Chambon and the Pas de l'Ours Landslides
TSM is designed as an end-to-end toolbox and its performance in terms of detection of terrain motion is evaluated on two image time series acquired on two unstable slopes characterized by slow Remote Sens. 2019, 11, 2189 8 of 23 motion (Chambon use case) and moderate motion (Pas de l'Ours use case) and by different viewing geometries. For the two use cases, the cameras are placed horizontally in front of the slope at distances (in the line of sight) comprised between 440 m and 640 m for the Chambon landslide and between 80 m and 560 m for the Pas de l'Ours landslide. The ground pixel sizes vary from 12.10 −2 m 2 to 31.10 −2 m 2 for the Chambon landslide and from 2.10 −2 m 2 to 18.10 −2 m 2 for the Pas de l'Ours landslide. The image acquisition frequency is variable and designed according to the landslide velocity. The choice of the focal length depends on the requested ground pixel size taking into account that focal lengths below 24 mm induces high distortion if the order of magnitude of the surface displacements is greater than all the errors related to the camera. To estimate the sensor orientation, six GCPs are available for the Chambon use case, and nine GCPs are available for the Pas de l'Ours use case. The GCPs are picked manually in an image from the fixed camera. The picking accuracy is about 2 pixels for the Chambon and 5 pixels for Pas de l'Ours. The standard deviations are 0.4 pixels for the Chambon and 5.6 pixels for the Pas de l'Ours use cases, while RMSE xy are respectively 1.4 pixels and 10.0 pixels. This difference is explained by the picking accuracy and the distribution of GCPs in the images of the two use cases.
The ISO sensitivity and the aperture opening are fixed at 100 and 5.6 respectively, and are kept constant for all acquisitions.  [49]. An acceleration of the landslide occurred in July 2015 with displacement in the range of 60 centimeter per day [50]; on the 26th of June 2015, mitigation measures with the purge of the moving mass was realized.
The images are acquired with a digital single-lens reflex camera (Canon 100D) of 18 MPix connected to a Paratronic LNS datalogger ( Figure 4c). For the period 15 February to 23 May 2017, eighth images were acquired per day in JPG file format, corresponding to a total of 830 images. To be in agreement with the terrestrial laser scanning point clouds used to validate the displacement rates, the image dataset used in this work corresponds to 26 images acquired at 14h UTC from 9 March to 17 May 2017 (Section 6.1). During this period, the landslide displacement rate is in the range of a centimeter per day. Figure 4a indicates the footprint of the camera view on the ground.

Application to a Moderate Motion Landslide: the Pas de l'Ours Use Case
The Pas de l'Ours landslide is located in the Guil valley (Hautes-Alpes, France). The landslide is about 1 km wide and 2 km long and is overhanging the Guil River. It involves lustrated clay-shales of lower Cretaceous and unconsolidated moraine formations. The main scarp visible on Figure 4b corresponds to the boundary of a paleo-landslide. In 2014, the landslide started to reactivate with some rockfalls triggered from the main scarp. In March 2017, acceleration of the central and lower part of the slope started with velocities up to several tens of centimeter per day. Currently, the landslide motion extends uphill. The depth of the rupture surface is on average 30 m, and the total volume is estimated at almost 15M m 3 .
The images are acquired with the same acquisition setup as for the Chambon use case. In order to track fast motion of the slope from 8 June to 4 October 2018, 10 images were acquired per day in the Nikon native CR2 file format further converted in JPG file format. The dataset consists of a total of 757 images for this period. The dataset acquired for the period 8 June to 18 July 2018 (22 images according to the selection criteria, Section 2.2) is presented in this work to be in agreement with the ancillary dataset (tacheometer-derived displacement rates on benchmarks) used to assess the results (Section 6.2). Figure 4b indicates the footprint of the camera view on the ground.

Sensitivity Analysis
The sensitivity of TSM to several processing parameters is evaluated. The influences of the image modalities (Section 5.1), of the image matching parameters (Section 5.2), of the size of the stable part used for the co-registration (Section 5.3) and of the influence of the correlation technique (Section 5.4) are analyzed on both datasets.

Sensitivity to Image Modalities
Sensitivity 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 [51]: where I is the brightness image and R, G and B are the three bands [10,31]; 3. the "average of p images A*" modality corresponds to a stack of p = 3 images used to reduce radiometric noise; 4. the "texture-derived image T*" corresponds to a Sobel filter applied to the grayscale converted images using a 3 × 3 convolution matrix [9,52]. The Sobel filter is used to highlight natural edges and to improve image quality which can be degraded by weather conditions or illumination [53].
These four modalities can be applied either on the initial RGB image or on the enhanced image. The 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.
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)). As outliers due to matching errors need to be filtered to avoid bias in the co-registration accuracy [33], the RMSExy is calculated for a Gaussian distribution filtered from the outliers at a 99% confidence threshold. The percentage of pixels with a C > C min is also calculated. C min corresponds to the correlation threshold: pixels with C < C min have no influence on the matching cost function. We estimate the impacts of the image modality on non-ideal images affected by some illumination changes ( Figure 5). The analysis indicates that compared to the first modality G, the modality A decreases the RMSExy parameter in both cases and increases the percentage of pixels with C > C min . The used of enhanced images have little influence on the two parameters with a difference smaller than 1/10th of a pixel. Same results are obtained for the modality B*. The modality T* decreases the percentage of pixels with C > C min ; less than 10% of the pixels has impacts on the correlation calculation.

Sensitivity to the Image Matching Parameters
The size of the research window (Inc), the size of the correlation window (SzW) (related to the spatial resolution of the image; [43,54]), the value of the regularization parameter (Reg) of the cost function [35] and the value of the correlation threshold (C min ) are evaluated. The parameter sub-pixel precision (spp) and the parameter γ are considered non-sensitive [48] and thus are set at spp = 0.1 pixel (by default the value equals 0.05 pixel but to reduce the computational costs, it is increased) and at γ = 2. changes ( Figure 5). The analysis indicates that compared to the first modality G, the modality A decreases the RMSExy parameter in both cases and increases the percentage of pixels with C > Cmin. The used of enhanced images have little influence on the two parameters with a difference smaller than 1/10 th of a pixel. Same results are obtained for the modality B*. The modality T* decreases the percentage of pixels with C > Cmin; less than 10% of the pixels has impacts on the correlation calculation.

Sensitivity to the Image Matching Parameters
The size of the research window (Inc), the size of the correlation window (SzW) (related to the spatial resolution of the image; [43,54]), the value of the regularization parameter (Reg) of the cost function [35] and the value of the correlation threshold (Cmin) are evaluated. The parameter sub-pixel precision (spp) and the parameter γ are considered non-sensitive [48] and thus are set at spp = 0.1 pixel (by default the value equals 0.05 pixel but to reduce the computational costs, it is increased) and at γ = 2.
To analyze the influence of those parameters on the correlation results and to choose the best combination to extract the most relevant and accurate information as possible, several correlations were carried out on one image pair of both landslides ( Figure 6). Value fixed arbitrarily at = Figure 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.
To analyze the influence of those parameters on the correlation results and to choose the best combination to extract the most relevant and accurate information as possible, several correlations were carried out on one image pair of both landslides ( Figure 6). Value fixed arbitrarily at Inc = 20, SzW = 4, Reg = 0.2, and Cmin = 0.3 are used as reference for the Chambon case and at Inc = 20, SzW = 8, Reg = 0.6, and Cmin = 0.5 for the Pas de l'Ours case.
Three tests were carried out to evaluate the sensitivity of the window size of the research area at the Chambon landslide. The sizes 4 × 4, 8 × 8 and 20 × 20 were tested. Seven tests were realized at the Pas de l'Ours landslide where the sizes 8 × 8, 12 × 12, 16 × 16, 20 × 20, 24 × 24, 32 × 32 and 40 × 40 were tested. This parameter is directly linked to the detection accuracy [9]. As shown in Figure 6b (left), higher values of Inc increase the percentage of pixels with C > C min . Above Inc = 20, the influence of this parameter diminishes. In Figure 6a (left), the parameter has no influence. Four tests (Chambon) and six tests (Pas de l'Ours) were carried out to define the influence of the correlation window size (Figure 6a,b middle left). The sizes 4 × 4, 6 × 6, 8 × 8, 10 × 10 were tested for both cases and the sizes 12 × 12 and 14 × 14 were also tested for the Pas de l'Ours case. As the influence of the correlation window is well-known [6], the parameter is represented against the percentage of pixels with C > C min . Compared to the Inc parameter, the SzW parameter has a small influence on the quality of the results. Four (Chambon) and five tests (Pas de l'Ours) were carried out to quantify the influence of the correlation threshold C min (Figure 6a, 6b middle right), by varying the parameter value from 0.3 to Remote Sens. 2019, 11, 2189 12 of 23 0.8 (Chambon) and from 0.1 to 0.9 with a step of 0.2 (Pas de l'Ours). 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 6d). The amount of noise and outliers is decreasing when the parameter increases because of the spatial smoothing of the displacement field [43] 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, Section 3). 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 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's shown in Figure 7a, the simultaneous use of stable areas in the near-and far-fields allows to increase the standard deviation and the average shift.  Then, the stability of the RMSE xy parameter (Equation (1)) along time, is used as criteria to compare the three combination strategies along the time series. Figure 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 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 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. 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.

Sensitivity to the Image Combination Strategy
The influence of forward and backward calculation (Figure 3) 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, 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 9. The RMSExy parameter (Equation (2)-(3)) 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 9a). As the time series stability of the Pas de l'Ours use case (Figure 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. 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 9. The RMSExy parameter (Equation (2)-(3)) 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 9a). As the time series stability of the Pas de l'Ours use case (Figure 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 1) which seemed the most appropriate to the use cases.  Afterward, we used a combination of the parameters (Table 1) which seemed the most appropriate to the use cases.

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 image spanning from 8 to 27 March 2017 is processed to illustrate the TSM application and characterize the motion of the Chambon landslide. Figure 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 Figure 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 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 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) increase. This pendulum motion is in agreement with a rotational movement which is confirmed by topographic surveys [49,55]. 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 by Point et al. (submitted) consisting in (1) registering the point clouds using an Iterative-Closest Point method (ICP), (2) generating 3D displacements fields on the point clouds. The point-clouds derived displacements are compared to the TSM-derived displacement on images acquired between the 9 March and the 8 April 2017. As there is no detectable movement between the 8th of April and the 10th of April according to the Figure 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. 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 by Point et al. (submitted) consisting in (1) registering the point clouds using an Iterative-Closest Point method (ICP), (2) generating 3D displacements fields on the point clouds. The point-clouds derived displacements are compared to the TSM-derived displacement on images acquired between the 9 March and the 8 April 2017. As there is no detectable movement between the 8 th of April and the 10 th of April according to the Figure 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 and 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 12 shows an example of displacement

Displacement Fields and 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 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 .
Two surface velocity profiles are presented in Figure 13: P1 is a longitudinal profile of 198 m in the N-S direction and P2 is a transverse profile of 132 m in the W-E direction (Figure 13a). The displacements are presented for the period 6-22 June 2018. The image time series was processed with the VSC method in the forward direction. One image per day is used. The profile P1 reveals the presence of two zones with > 0.5 m·day −1 ): the first one situated along the first 25 m and the second one from 50 to 175 m (Figure 13b). The higher displacements are measured during the period 6-12 June 2018. Significant changes (loss of spectral coherence) of the ground surface is responsible of the lack of information between 80 and 135 m for the period 12-16 June 2018. In Profile 2, the main active zone is clearly identified from 20 to 100 m. The stable part cannot be determined directly in Profile 2 because low displacement rates are measured close to the Eastern limit. 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 .  One target was measured by a total station during the same time span than the image time series. It corresponds here to the point Pt0 in Figure 8. Its cumulative displacement is represented in Figure  14 according to other specific points (located in Figure 8) whose displacements were calculated with TSM. As we can see, displacements measured by the total station at the Pt0 are in the error bar of the TSM processing. The standard deviation linked to the difference between the displacements One target was measured by a total station during the same time span than the image time series. It corresponds here to the point Pt0 in Figure 8. Its cumulative displacement is represented in Figure 14 according to other specific points (located in Figure 8) whose displacements were calculated with TSM. As we can see, displacements measured by the total station at the Pt0 are in the error bar of the TSM processing. The standard deviation linked to the difference between the displacements measured at Pt0 and calculated by the TSM is equal to 16 cm.

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 only 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 constrain 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 < 24mm) 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 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  Figure 3) measured by TSM and by a total station at the Pas de l'Ours landslide.

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 only 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 constrain 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 < 24mm) 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 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 ( [49,50] 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).