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

Abstract

Passive sensors such as multi-spectral (e.g., Single Lens Reflex, SLR) cameras are increasingly being used for geohazards monitoring (landslides, cliffs affected by rock falls, ice glaciers, and volcano flanks) 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 configurations. 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⁻¹), and the Pas de l’Ours landslide (Hautes-Alpes, France) characterized by moderate motion (> 50 cm·day⁻¹). 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.

1. 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 development of relatively low-cost acquisition systems (very-high resolution consumer grade camera, > 20 MPix digitizers), it allows generating long- and high-frequency 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 [9,10,11,12,13]. Simple permanent installations requiring only solar panels and a stable anchoring can be easily setup [5].

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 [14]. The most commonly used techniques for surface motion measurements are in-situ point-based sensors such as GNSS, inclinometers, extensometers and crack-meters [15,16], 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 [17,18,19]. Several studies also demonstrated the possibility to estimate horizontal displacement measurements from remote platforms such as optical airborne/satellite images [20,21] and more recently from SAR amplitude images using offset-tracking techniques [22]. 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; [9]), the reconstruction of the surface topography [23,24] or the identification of discontinuities in the unstable mass [25]. Recently, Gance et al. [5], Travelletti et al. [9], Roncella et al. [26] and Mishra et al. [27] demonstrated that a fixed installation of a camera allows obtaining time series of the displacement fields from multiple pairwise image correlation. The obtained displacement fields are able to 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. [9] 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 [28,29] where long-time series of video frames are processed using Particle Image Velocimetry (PIV) approach [30]. 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 possibly adapted to any site configuration (slope dimensions, slope deformation rates, distance between the slope and the camera) and any camera type (sensor size, focal length) do exist [31,32].

Our works target the presentation of the “TSM—Tracing Surface Motion” processing toolbox. TSM is a generic and fully-automated toolbox, built on previous work [5,9,29,33], 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) 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).

2. 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].

2.1. 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):

(1)$\left\{\begin{array}{c}u=u_0-f_u\frac{m_{11}\left(X-X_{cam}\right)+m_{12}\left(Y-Y_{cam}\right)+m_{13}\left(Z-Z_{cam}\right)}{m_{31}\left(X-X_{cam}\right)+m_{32}\left(Y-Y_{cam}\right)+m_{33}\left(Z-Z_{cam}\right)}\\ v=v_0-f_v\frac{m_{21}\left(X-X_{cam}\right)+m_{22}\left(Y-Y_{cam}\right)+m_{23}\left(Z-Z_{cam}\right)}{m_{31}\left(X-X_{cam}\right)+m_{32}\left(Y-Y_{cam}\right)+m_{33}\left(Z-Z_{cam}\right)}\end{array}\right.,$

where $\left(u,v\right)$ are the image-based coordinates, $\left(X,Y,Z\right)$ are the world 3-D coordinates ($X$ is the East direction, $Y$ is the North direction and $Z$ is the elevation), $\left(X_{cam},Y_{cam},Z_{cam}\right)$ is the position of the camera in the world coordinate system, $\left(f_u,f_v\right)$ is the effective focal length and $\left(m_{11},m_{12},\dots m_{33}\right)$ 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 $\left(u_0,v_0\right)$ 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):

(2)$RMS{E}_{xy}=\sqrt{\frac{1}{n}\times \sum \mathrm{∆}{x}^2+\mathrm{∆}{y}^2},$

(3)$\left\{\begin{array}{c}\mathrm{∆}{x}_i=x_{proj,i}-x_{picking,i}\\ \mathrm{∆}{y}_i=y_{proj,i}-y_{picking,i}\end{array}\right.,$

where $\left(x_{proj,i},y_{proj,i}\right)$ are the coordinates of the $i$th re-projected GCP in the image-based coordinate system, $\left(x_{picking,i},y_{picking,i}\right)$ are the coordinates of the $i$th picked GCP in the image-based coordinate, and $n$ is the number of GCPs.

2.2. 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.

2.3. 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 $\Delta X_{DISP}$ and $\Delta Y_{DISP}$ and a Correlation Coefficient (C) grid (Figure 2).

The systematic offsets ($\Delta X_{DISP}$,$\Delta Y_{DISP}$) are modelled considering two planes $\Delta X_{model }$ and $\Delta 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)):

(4)$\left\{\begin{array}{c}\Delta x_{model,i}=a+b\times x_{DISP,i}+c\times y_{DISP,i}\\ \Delta y_{model,i}=a^{\prime }+b^{\prime }\times x_{DISP,i}+c^{\prime }\times y_{DISP,i}\end{array}\right.,$

where $\Delta x_{model,i}$ and $\Delta y_{model,i}$ represent the modelled offsets of the $i$th pixel and $x_{DISP,i},y_{DISP,i}$ are the raw disparities of the $i$th pixel. The parameters $a$, $b$, $c$, $a^{\prime }$, $b^{\prime }$, $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):

(5)$\left\{\begin{array}{c}\Delta X_{CORRECTED}=\Delta X_{DISP}-\Delta X_{MODEL}\\ \Delta Y_{CORRECTED}=\Delta Y_{DISP}-\Delta Y_{MODEL}\end{array}\right.,$

where $\Delta X_{CORRECTED}$ and $\Delta 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.

2.4. 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):

(6)$E_{\alpha }=\sum \left(1-Corr\right)+\alpha \times F,$

where $\left(1-Corr\right)$ is the similarity measure with $Corr$ being the NCC score; $\alpha$ 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.

2.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.

2.6. 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 ($\Delta u$, $\Delta v$) to the absolute coordinate system in terrain geometry ($\Delta X,\Delta Y,\Delta 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.

3. 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 $\left[I,I+n\right]$ 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 $\left[I,I+n\right]$ in forward and backward mode (Figure 3d).

4. 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 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⁻² m² to 31.10⁻² m² for the Chambon landslide and from 2.10⁻² m² to 18.10⁻² m² 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.

4.1. Application to a Slow Motion Landslide: The Chambon Use Case

The Chambon landslide (Figure 4a) is located in the Romanche valley (Isère, France), downhill of a tunnel allowing the Road RD1091 to cross the valley above the Chambon lake. The landslide involves sedimentary rocks from the Lias and extends over a length of about 170 m and a width of 92 m. The depth of rupture surface was estimated by geophysical surveys at 25 m and the volume at about 6000.000 m³ [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.

4.2. 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³.

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.

5. 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.

5.1. 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]: $I=0.30\times R+0.59\times G+0.11\times B$ 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.

5.2. 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.

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,\mathrm{and} Cmin=0.3$ are used as reference for the Chambon case and at $Inc=20,SzW=8,Reg=0.6,\mathrm{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 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.

5.3. 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 $\left(u,v\right)$ 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.

5.4. Sensitivity to the Image Combination Strategy

Several image combination strategies (CMC, VSC, and RVSC) are evaluated (Figure 3) 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 (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⁻¹), 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⁻¹), 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 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, 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 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.

6. 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.

6.1. 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⁻¹) 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⁻¹; the stable parts are also highlighted with very small displacement rates lower than 0.6 cm·day⁻¹ 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.

6.2. 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⁻¹) 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⁻¹ and the upper part which presents high displacements > 0.10 m·day⁻¹ to > 2.0 m·day⁻¹. 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⁻¹.

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⁻¹): 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.

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.

7. 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⁻¹. 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.

8. 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).

Author Contributions

M.D. developed the processing pipeline and performed the experiments. J.-P.M., O.B. and L.L. defined the functions of TSM and designed the experiments. A.S. supports the development of the code. M.D., O.B. and J.-P.M. collected the data in the field. J.P. processed the terrestrial laser scanning point clouds of the Chambon landslide. M.D. and J.-P.M. prepared the manuscript with contributions from all co-authors. L.L. and J.-P.M. provided financial support and resources. Access to the code is possible for science applications by contacting J.-P.M.

Funding

This research received no external funding.

Acknowledgments

This work has been supported by CNRS/Institut de Physique du Globe de Strasbourg (IPGS) and SAGE Society as part of a CIFRE/ANRT research contract. The authors acknowledge P. Bornemann (IPGS/LIVE) for contributing to the terrestrial laser scanning acquisition at the Chambon landslide, P. Diot (RTM-05) for support in the field at the Pas de l’Ours landslide, and S. Masse (CCGQ) and V. Ségel (RTM-05) for access to the total station dataset at the Pas de l’Ours landslide.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Table

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

Figure 2. Details of the TSM processing toolbox for the Modules 3 to 6.

Figure 3. 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. h1, h2 and h3 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.

Figure 4. Landslide use cases: (a) (left) The Chambon landslide (Isère, France) with a terrestrial laser scanner derived Digital Surface Model (DSM) of the slope and (right) an example of image acquired by a fixed monoscopic camera. (b) (left) The Pas de l’Ours landslide (Hautes-Alpes, France) with an airborne LiDAR-derived Digital Surface Model (DSM) of the slope and (right) an example of image acquired by a fixed monoscopic camera. (c) Example of the image acquisition system developed by SAGE.

Figure 5. Sensitivity of co-registration quality to the image modality expressed by the criteria RMSExy and the criteria percentage of pixels with C > Cmin 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 6. Sensitivity of co-registration quality to the correlation parameters expressed by the criteria percentage of pixels with C > Cmin 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 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 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 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 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 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 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 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 14. Cumulated displacements [m] of specific points (cf. Figure 3) measured by TSM and by a total station at the Pas de l’Ours landslide.

Table 1

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

Use Case	Image Modality	Image Matching Parameters [pix]				Size of the Stable Area [%]	Image Combination Strategy
		Inc	SzW	Reg	Cmin
Chambon landslide	*G	20 × 20	4 × 4	0.2	0.5	36	VSC
Pas de l’Ours landslide	*G	24 × 24	4 × 4	0.2	0.5	9	VSC