# Sparsity Optimization Method for Slow-Moving Landslides Detection in Satellite Image Time-Series

Mai Quyen I Pham, Pascal Lacroix, Marie Pierre Doin

► To cite this version:

Mai Quyen I Pham, Pascal Lacroix, Marie Pierre Doin. Sparsity Optimization Method for Slow-Moving Landslides Detection in Satellite Image Time-Series. IEEE Transactions on Geoscience and Remote Sensing, 2018, pp.1-12. 10.1109/TGRS.2018.2871550. hal-02062494

HAL Id: hal-02062494

https://hal.science/hal-02062494v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sparsity Optimization Method for Slow-Moving Landslides Detection in Satellite Image Time-Series

Mai Quyen Pham<sup>1, ∗</sup>, Member, IEEE, Pascal Lacroix, and Marie Pierre Doin

**Abstract**—This paper presents a new method based on recent optimization technique to detect slow-moving landslides (<150 m/year) in time series of displacement field generated by satellite images. Sparse optimization is applied simultaneously on the 3-D data set in space as well as in time. The proposed method takes into account the distinctive signal physical properties in space and time, by enforcing a sparse representation by blocks in space, but a continuing and monotonous representation in time of the landslides. As a result, we show that a mixed $ℓ_{1,2}$-norm is the most suitable norm for this detection problem, compared to pure $ℓ_1$-norm or $ℓ_2$-norm. Moreover, an outlier estimation step is included that sets apart the Gaussian noise from locally sparse processing errors in the data. The performance of this approach is tested by applying it both on synthetic data and on a time series of displacements fields over 16 dates in the Colca Valley, Peru. This detection presents commission and omission errors for landslides of 29% and 14%, respectively, using a medium resolution (10 m) data set of optical satellite images. It detects all important landslides, already known from field investigations. Moreover, it also points out other smaller or unknown landslides, increasing the existing slow-moving landslide inventory by +50%.

***Index Terms***—Detection, optimization, satellite image time-series, segmentation, slow-moving landslides, sparsity.

## I. INTRODUCTION

LANDSLIDES are a major threat in all mountainous environments of the world, claiming thousands of lives each year [1]. Their triggering is controlled by different environmental (rainfalls and earthquakes) or anthropogenic forcings [23], and predicting their time and location is a very complex task. However, in numerous cases, rapid and possibly catastrophic landslides (m/sec) are preceded by a phase of slow creep [18]. One of the best examples of this prefailure motion and activity was found on the Vajont landslide in Northern Italy, which caused 1900 casualties [19]. These observations show the importance of detecting and monitoring gravitational slow motions.

Remote sensing techniques are appropriate tools for that purpose. Both InSAR techniques [2], [20] and optical image

Manuscript received December 12, 2017; revised July 15, 2018; accepted August 7, 2018. This work was supported by European Space Agency through Alcantara Initiative under Grant 15/P26. (Corresponding author: Mai Quyen Pham.)

The authors are with Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTerre, 38000 Grenoble, France (e-mail: mai-quyen.pham@univ-grenoble-alpes.fr).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TGRS.2018.2871550

correlation [10], [16], [22], [37] have been applied for very slow and slow-moving landslide (typically between 1 cm/year and 150 m/year [9]) monitoring. Optical satellites are particularly well adapted to monitor objects of small size such as landslides, moving at centimeter to meter rates per year. Image correlation allows to retrieve displacement fields with one-fifth of a pixel accuracy, thus to detect slow-moving landslides for a large range of image resolution, depending on the landslide velocity and separation time between the acquisitions.

The increasing number of satellite-based earth observation images, with open access to medium resolution (10 m) data at a global scale [Landsat-7/8, Sentinel-2, World Heritage Satellite Pour l'Observation de la Terre (SPOT) Program], enables the creation of multiple stacks of displacement maps that can be the base for gravitational motion detection and monitoring [11], [12], [15], [17]. Series of pairwise correlation maps obtained between any combination of images in a stack can also be inverted to obtain displacement time series, using redundancy to improve measurement uncertainty (e.g., [12]). Whereas many methods exist to automatically detect rapid landslides (m/sec) [9], based for instance on morphology or vegetation changes before/after the event, very few methods have been developed for the automatic detection of active slow-moving landslides. However, it should be noted that some methods based on surface texture and morphology have captured areas slowly deforming under gravity [21]. To our knowledge, only one method has been developed specifically for the semiautomatic detection of slow-moving landslides in optical images [38]. It is a pixel-based analysis of a time series of ground motion measurements that detects slow-moving landslides by thresholding an indicator of the coherency of the displacement vector directions. Its application to very-high resolution (VHR) imagery has shown the method ability to detect several tens of slow-moving landslides in an area of 150 km². As all pixel-based methods of detection (e.g., [24]), it provides a salt and pepper result that is then filtered by postprocessing image analysis. This method, therefore, relies on the tuning of a set of parameters. Here, we propose a new method of active slow-moving landslide detection, based on sparse analysis of image time-series that overcomes the pixel-based analysis and limits the number of parameters.

Sparse representation-based solutions have been successfully used in many fields such as image processing, pattern
recognition, and signal processing. They have also been applied in [26], [28], and [25] and, more specifically, to land- slide studies such as [29] that used the wavelet transform to achieve data fusion between LiDAR and QuickBird imagery. The sparse representation techniques can be classically divided into two steps. The first step is the data transformation into a new domain (Gradient, Fourier transform, Wavelet trans- form, etc.) where the information of importance is obtained. The second step is the evaluation and processing of these information through a sparsity measure such as the ℓ₁-norm, the ℓₚ, 0 < p < 1 quasi-norms or the mixed-norms [26], and the ratio-norms [25], [27].

In this paper, we will first expose the basis of the developed method and then show its applicability to the detection of slow-moving landslides in a time series of medium resolution (10 m) SPOT satellite images, as a proof of concept for its use with Sentinel-2 (10-m resolution) and LandSat8 (15-m resolution) data. Although the method is designed for processing 3-D data, a main development includes the separation of the optimization in space and time based on the physical properties of the searched signal. Especially, this paper uses the mixed-norm ℓ₁,₂ [30]–[32] over the space and time dimensions for measuring sparsity of the signal. In addition, our approach separates the signal related to the landslides from the noise and outliers using the structured pattern of the outliers and the monotony in time of the landslide signals. Taking into account outliers allows avoiding many false alarms.

This paper is organized as follows. Following this Introduc-
tion, Section II describes the assumptions defining the observa-
tion model. Section III shows the optimization method together
with the criterion to minimize and details the proposed algo-
rithm. The validation of the proposed method is assessed in
Section IV. Finally, the conclusions and perspectives are drawn
up in Section V.

II. OBSERVATION MODEL

The proposed method takes as input a time series of ground-
displacement maps. Here, we will use displacement maps
coming from the correlation of diachronic satellite optical
images [13]. Other sources of data, as InSAR, [2] could also
be used. The image correlation leads to a map of displacement
along the rows North–South (NS) and the columns East–West
(EW) of the map. To simplify our data set, we assume that
the movement of slow-moving landslides is oriented along the
topographical slope direction, calculated here over a kernel
size of 30 m. This classical and realistic assumption (e.g., [38])
enables us to project the EW and NS components along the
slope, using a digital elevation model (DEM) with a 10-m
pixel size derived from a stereo pair of Pléiades images [10],
and to process only the projection of the displacement vector
along this local direction.

We consider the time series of displacement fields
z ∈ ℝ^M×N×T as the data observation, where T represents
the number of images and (M, N) are the size of each
image. For simplification reasons, ℝ^M×N×T is denoted 𝒳.
During the time study T, the data set displays, beside the
ground movement that we want to extract (e.g., landslides),
different effects due to important natural or anthropogenic

<figure><img src="image_2.png" /><figcaption>Fig. 1. Cropped displacement time series on a small landslide in the Colca Valley. Left: examples of displacement maps at three time steps. Right: time series of the z (input), x̄ (signal related to the landslide) and, ȳ (outlier) of the green point displayed on the left.</figcaption></figure>

changes (e.g., clouds, new buildings and roads, agricultural
zones, vegetation growth, shadows, etc.), and noise that limits
the correlation accuracy due to radiometric coding, image
structure, and so on. The purpose of the paper is to extract
from the input, here noted **z**, the displacement signals related
to the landslides. We suppose that **z** can be decomposed into
three parts: 1) the first one is the signal of interest **x** ∈ 𝒳,
related to the landslides; 2) the second one, denoted by **y** ∈ 𝒳,
encompasses all the outliers due to bias in the correlation
step or clouds or changes in the ground cover; and 3) the last
unstructured contribution is gathered in the noise term **n** ∈ 𝒳.
More precisely, our problem can be formulated mathematically
as follows (Fig. 1):

$$
\mathbf{z} = \bar{\mathbf{x}} + \bar{\mathbf{y}} + \mathbf{n}. \tag{1}
$$

Detecting the landslides **x** by using the ill-posed problem
of only one equation (1) with three unknowns requires the
incorporation of all plausible knowledge on those three parts.
We can list the *a priori* regularization conditions according
to the physical properties of each contribution as follows.

1) Noise **n**: The additive noise **n** is a realization of a zero-mean white Gaussian noise.

2) Landslide **x**:

a) Limited number of landslides whose pixels are
connected (sparse in space).

b) Landslide volume is preserved during the evolu-
tion.

c) Landslide moves monotonically in time, i.e., the derivative of the displacement with time is always either positive or negative.

3) Structured noise **y**:

a) Sparse in time.

b) Structured in space due to the large size of the cor-
relation operator used to compute z, often between
16 and 64 pixels [12].

In Section III, we will explain in detail the proposed method to
solve the problem by using the above-mentioned assumptions.

III. PROPOSED METHOD

A. Criterion

The purpose of this paper is to separate **x** and **y** within **z**
through (1), which leads to an ill-posed problem. To solve
this, we propose an optimization approach that minimizes the following criterion:

$$ \text{Find } (\hat{\mathbf{x}}, \hat{\mathbf{y}}) \in \underset{\mathbf{x} \in C, \mathbf{y} \in \mathcal{X}}{\operatorname{argmin}} (\psi(\mathbf{x}, \mathbf{y}) + \varphi(W\mathbf{x}) + \rho(\mathbf{y})) \quad (2) $$

where the data fidelity term is defined by function $\psi$. Functions $\varphi$ and $\rho$ design the regularizations on the landslide $\mathbf{x}$ and the structured noise $\mathbf{y}$, respectively. Finally, the other a priori on the landslides are modeled by the nonempty closed convex set $C$. We will now turn our attention to the choice of each term.

1) Data fidelity term: Function, $\psi: \mathcal{X} \to \mathbb{R}$, accounts for the noise $n$. As mentioned earlier, this noise is assumed to be additive, zero-mean, white, and Gaussian. This leads to the following quadratic form:

$$ \psi(\mathbf{x}, \mathbf{y}) = \frac{1}{2} \|\mathbf{z} - (\mathbf{x} + \mathbf{y})\|^2. \quad (3) $$

2) A regularization on $\mathbf{x}$: Function, $\varphi: \mathcal{X} \to \mathbb{R}$, takes into account the sparse spatial representation of the landslides $\mathbf{x}$ via a linear operator $W$. More precisely, because of the structured aspect of the landslide signal in space and time, we choose to represent the sparsity on the detected landslide contour by defining an operator $W$ that operates locally on patches of each image. The simplest choice for this weighted operator is expressed as follows:

$$ (W\mathbf{x})_{m,n}^t = x_{m,n}^t - \frac{1}{P^2} \sum_{i=m-\left\lfloor \frac{P}{2} \right\rfloor}^{m+\left\lfloor \frac{P}{2} \right\rfloor-1} \sum_{j=n-\left\lfloor \frac{P}{2} \right\rfloor}^{n+\left\lfloor \frac{P}{2} \right\rfloor-1} x_{i,j}^t \quad (4) $$

where $\lfloor \cdot \rfloor$ rounds the element to the nearest integer less than or equal to that element, $(m, n) \in \{1, 2, \dots, M\} \times \{1, 2, \dots, N\}$ and $t \in \{1, 2, \dots, T\}$ express the spatial and time increments, respectively, and $P$ is the size of the patches. In the following, we denote $\mathcal{M}$, $\mathcal{N}$, and $\mathcal{T}$ by $\{1, 2, \dots, M\}$, $\{1, 2, \dots, N\}$, and $\{1, 2, \dots, T\}$, respectively. Next, we propose to use the two-level mixed-norm $\ell_{1,2}^2$ for the regularization $\varphi$. This choice is conducted because such norm has the ability to structure the prior, in order to incorporate the different sparse representation level of the landslide along space and time. More precisely, for every $\mathbf{s} = W\mathbf{x} \in \mathcal{X}$, we define the two-level mixed-norm value by:

$$ \varphi(\mathbf{s}) = \lambda \sum_{t \in \mathcal{T}} \left( \sum_{(m,n) \in \mathcal{M} \times \mathcal{N}} |s_{m,n}^t| \right)^2 \quad (5) $$

with $\lambda \in ]0, +\infty[$, the regularization parameter related to the regularization term of $\mathbf{x}$. If the pixels represent the same connected landslide, then the corresponding weight will be small. Sparse spatial regularization is thus applied to landslide contours.

3) A regularization on $\mathbf{y}$: Function, $\rho: \mathcal{X} \to \mathbb{R}$, is a proper, lower semicontinuous, convex function that introduces the prior knowledge on the structured noise $\mathbf{y}$. As mentioned in Section II, $\mathbf{y}$ appears only a few times during the period of study, in contrast to the landslide signal.

Therefore, the sparsity of $\mathbf{y}$ is stronger in time than in space. Thus, we choose a two-level mixed-norm for this regularization $\rho$, with $\ell_1$ applied on the time and $\ell_2$ on the space dimension

$$ \rho(\mathbf{y}) = \tau \sum_{(m,n) \in \mathcal{M} \times \mathcal{N}} \left( \sum_{t \in \mathcal{T}} |y_{m,n}^t| \right)^2 \quad (6) $$

where $\tau \in ]0, +\infty[$ is the regularization parameter related to the regularization term of $\mathbf{y}$.

4) Apriori on the landslide: $C$ is a nonempty closed convex constraint set. This set takes into account the information concerning the physical properties of the landslides. Here, we consider the following two attributes.

a) First, the landslide volume is preserved during the evolution and the landslide moves monotonically in time. We thus can introduce

$$ C_1 = \{ \mathbf{x} \in \mathbb{R}^{M \times N \times T} \mid x_{m,n} \text{ is monotone in time} \\ \text{and, } |x_{m,n}^t - x_{m,n}^{t-1}| \le \kappa, \\ \forall (m, n, t) \in \mathcal{M} \times \mathcal{N} \times \mathcal{T} \} \quad (7) $$

where $\kappa \in [0, +\infty[$. For every $(m, n) \in \mathcal{M} \times \mathcal{N}$, the sign of the difference between $z_{m,n}$ at the date 1 and at the date $T$ is used to specify the monotony in time of $x_{m,n}$. Thus, $x_{m,n}$ can be an increasing or decreasing sequence in time. We denote $H$ the linear operator from $\mathcal{X}$ to $\mathbb{R}^{M \times N \times (T-1)}$ such that: for every $(m, n, t) \in \mathcal{M} \times \mathcal{N} \times \mathcal{T}$, $(Hx)_{m,n}^t = h_{m,n}(x_{m,n}^t - x_{m,n}^{t+1})$ with $h_{m,n} = \text{sign}(z_{m,n}^1 - z_{m,n}^T)$.

b) Second, we limit our analysis to the landslides that move smoothly for the first and the last dates. This attribute serves us to remove the outliers at the limits of the time series that can degrade the monotony analysis. This attribute is expressed in a mathematical way: for every $(m, n, t) \in \mathcal{M} \times \mathcal{N} \times \mathcal{T}$ and $\mathbf{x} \in \mathbb{R}^{M \times N \times T}$,

i) $x_{m,n}^1 - x_{m,n}^2 = a_{m,n}$

ii) $x_{m,n}^{T-1} - x_{m,n}^T = a_{m,n}$

where

$$ a_{m,n} = \frac{1}{T-3} (x_{m,n}^2 - x_{m,n}^{T-1}) $$

These a priori are equivalent to introducing the two following constraints:

$$ C_2 = \{ \mathbf{x} \in \mathcal{X} \mid x_{m,n}^1 - x_{m,n}^2 = a_{m,n}, \\ \forall (m, n, t) \in \mathcal{M} \times \mathcal{N} \times \mathcal{T} \} \quad (8) $$

and

$$ C_3 = \{ \mathbf{x} \in \mathcal{X} \mid x_{m,n}^{T-1} - x_{m,n}^T = a_{m,n}, \\ \forall (m, n, t) \in \mathcal{M} \times \mathcal{N} \times \mathcal{T} \}. \quad (9) $$

Consequently, we define $C = C_1 \cap C_2 \cap C_3$.
## B. Proximity Operators

In this section, we recall some mathematical tools in optimization used in this paper. Let $\Gamma_0(\mathbb{R}^Q)$ denote the class of proper lower semicontinuous convex functions from $\mathbb{R}^Q$ to $]-\infty, +\infty]$. The proximity operator [3], [4] of a function $f \in \Gamma_0(\mathbb{R}^Q)$ at $\bar{x} \in \mathbb{R}^Q$, denoted $\text{prox}_f : \mathbb{R}^Q \to \mathbb{R}^Q$, is defined as

$$ \forall \bar{x} \in \mathbb{R}^Q, \quad \text{prox}_f(\bar{x}) := \underset{x \in \mathbb{R}^N}{\text{minimize}} \frac{1}{2} \|x - \bar{x}\|^2 + f(x). \quad (10) $$

**Definition 1 (Indicator function):** Let $C$ be a nonempty subset of $\mathbb{R}^Q$. Its indicator function $\iota_C$ is defined as

$$ (\forall x \in \mathbb{R}^Q) \quad \iota_C(x) = \begin{cases} 0 & \text{if } x \in C \\ +\infty & \text{otherwise.} \end{cases} $$

In consequence, if $C$ is a nonempty closed convex set of $\mathbb{R}^Q$, then $\text{prox}_{\iota_C}$ is reduced to the projection $\Pi_C$ onto $C$. The proximity operator has numerous mathematical properties that can be found in [8] and [4].

Inserting the choice functions into the optimization criteria (Section III-A) and taking into account the above-mentioned definitions, problem (2) can be expressed again as

$$ \text{Find } (\hat{\mathbf{x}}, \hat{\mathbf{y}}) \in \underset{(\mathbf{x}, \mathbf{y}) \in \mathcal{X}^2}{\operatorname{argmin}} \theta(\mathbf{x}, \mathbf{y}) \quad (11) $$

where

$$ \theta(\mathbf{x}, \mathbf{y}) = \psi(\mathbf{x}, \mathbf{y}) + \varphi(W\mathbf{x}) + \rho(\mathbf{y}) + \iota_{\tilde{C}_1}(H\mathbf{x}) + \iota_{C_2}(\mathbf{x}) + \iota_{C_3}(\mathbf{x}) \quad (12) $$

$$ \text{and } \tilde{C}_1 = \{\tilde{\mathbf{x}} \in \mathbb{R}^{M \times N \times (T-1)} : \forall (m, n, t), 0 \le \tilde{x}_{m,n}^t \le \kappa\}. $$

## C. Algorithm

To solve our convex optimization problem, which involves the sum of four convex functions and various linear operators, we propose to use the Monotone+Lipschitz Forward-Backward-Forward (M + L FBF) algorithm (presented in Algorithm 1, and recently proposed in [5]). This algorithm belongs to the class of primal-dual algorithms and avoids some matrix inversions, necessary in many other approaches [6], [7]. It requires to compute the proximity operators of each of the terms in (11) (note that when one term has a Lipschitz continuous gradient, the gradient can be computed instead of the proximity operator).

We will denote by $W^*$ and $H^*$ the adjoint operators of $W$ and $H$, respectively. At each iteration $k$, the step size $\gamma_k$ must be chosen to satisfy the following rule: let $\beta = 1 + (\|W\|^2 + \|H\|^2 + 3)^{1/2}$ and $\epsilon \in ]0, (1/\beta + 1)[$, then $\gamma_k \in [\epsilon, (1 - \epsilon/\beta)]$. It must be emphasized that the choice of the step size is crucial for the convergence speed of the algorithm [33]–[36]. The functional (12) to be minimized being composed of convex functions as well as indicator functions of convex set, the algorithm typically requires to compute, in parallel, proximal operators and projection onto the closed convex set. Concerning proximity operators, closed-form expressions for a wide class of convex functions can be found in [8]. The projection onto $\tilde{C}_1$ is explicit and reduces to projections onto hypercube, that is,

$$ \Pi_{\tilde{C}}(\tilde{\mathbf{x}}) = (\Pi_{[0,\kappa]}(\tilde{x}_{m,n}^t))_{(m,n,t) \in \mathcal{M} \times \mathcal{N} \times \mathcal{T}} \quad (13) $$

## Algorithm 1 Primal-Dual Algorithm M + L FBF to Solve (11)

$$ \text{Let } \gamma_k \in [\epsilon, \frac{1-\epsilon}{\beta}], (\mathbf{x}_0, \mathbf{y}_0, \mathbf{v}_0^1, \mathbf{v}_0^2, \mathbf{v}_0^3, \mathbf{v}_0^4, \mathbf{u}_0) \in \mathcal{X}^7. $$

**for** $k = 0, 1, \dots, \text{do}$

**Gradient computation**

$$ \mathbf{r}_k^1 = \mathbf{x}_k - \gamma_k (\nabla \psi(\mathbf{x}_k) + W^* \mathbf{v}_k^1 + H^* \mathbf{v}_k^2 + \mathbf{v}_k^3 + \mathbf{v}_k^4) $$

$$ \mathbf{r}_k^2 = \mathbf{y}_k - \gamma_k (\nabla \psi(\mathbf{y}_k) + \mathbf{u}_k) $$

**Proximal compute**

$$ \mathbf{s}_k^1 = \mathbf{v}_k^1 + \gamma_k W \mathbf{x}_k $$

$$ \mathbf{p}_k^1 = \mathbf{s}_k^1 - \gamma_k \text{prox}_{\varphi} \left( \frac{\mathbf{s}_k^1}{\gamma_k} \right) $$

$$ \mathbf{e}_k = \mathbf{u}_k + \gamma_k \mathbf{y}_k $$

$$ \mathbf{q}_k = \mathbf{e}_k - \gamma_k \text{prox}_{\rho} \left( \frac{\mathbf{e}_k}{\gamma_k} \right) $$

**Projection compute**

$$ \mathbf{s}_k^2 = \mathbf{v}_k^2 + \gamma_k H \mathbf{x}_k $$

$$ \mathbf{p}_k^2 = \mathbf{s}_k^2 - \gamma_k \Pi_{\tilde{C}_1} \left( \frac{\mathbf{s}_k^2}{\gamma_k} \right) $$

$$ \mathbf{s}_k^3 = \mathbf{v}_k^3 + \gamma_k \mathbf{x}_k $$

$$ \mathbf{p}_k^3 = \mathbf{s}_k^3 - \gamma_k \Pi_{C_2} \left( \frac{\mathbf{s}_k^3}{\gamma_k} \right) $$

$$ \mathbf{s}_k^4 = \mathbf{v}_k^4 + \gamma_k \mathbf{x}_k $$

$$ \mathbf{p}_k^4 = \mathbf{s}_k^4 - \gamma_k \Pi_{C_3} \left( \frac{\mathbf{s}_k^4}{\gamma_k} \right) $$

**Averaging**

$$ \bar{\mathbf{p}}_k^1 = \mathbf{p}_k^1 + \gamma_k W \mathbf{r}_k^1 $$

$$ \mathbf{v}_{k+1}^1 = \mathbf{v}_k^1 - \mathbf{s}_k^1 + \bar{\mathbf{p}}_k^1 $$

$$ \bar{\mathbf{q}}_k = \mathbf{q}_k + \gamma_k \mathbf{r}_k^2 $$

$$ \mathbf{u}_{k+1} = \mathbf{u}_k - \mathbf{e}_k + \bar{\mathbf{q}}_k $$

$$ \bar{\mathbf{p}}_k^2 = \mathbf{p}_k^2 + \gamma_k H \mathbf{r}_k^1 $$

$$ \mathbf{v}_{k+1}^2 = \mathbf{v}_k^2 - \mathbf{s}_k^2 + \bar{\mathbf{p}}_k^2 $$

$$ \bar{\mathbf{p}}_k^3 = \mathbf{p}_k^3 + \gamma_k \mathbf{r}_k^1 $$

$$ \mathbf{v}_{k+1}^3 = \mathbf{v}_k^3 - \mathbf{s}_k^3 + \bar{\mathbf{p}}_k^3 $$

$$ \bar{\mathbf{p}}_k^4 = \mathbf{p}_k^4 + \gamma_k \mathbf{r}_k^1 $$

$$ \mathbf{v}_{k+1}^4 = \mathbf{v}_k^4 - \mathbf{s}_k^4 + \bar{\mathbf{p}}_k^4 $$

**Update**

$$ \mathbf{x}_{k+1} = \mathbf{x}_k - \gamma_k (\nabla \psi(\mathbf{r}_k^1) + W^* \mathbf{p}_k^1 + H^* \mathbf{p}_k^2 + \mathbf{p}_k^3 + \mathbf{p}_k^4) $$

$$ \mathbf{y}_{k+1} = \mathbf{y}_k - \gamma_k (\nabla \psi(\mathbf{r}_k^2) + \mathbf{q}_k) $$

**end**

where $\Pi_{[0,\kappa]}(\tilde{x}_{m,n}^t) \text{ is the projection of } \tilde{x}_{m,n}^t \in \mathbb{R} \text{ onto the compact interval } [0, \kappa]$, defined by

$$ (\forall x \in \mathbb{R}) \quad \Pi_{[0,\kappa]}(\tilde{x}) = \begin{cases} 0 & \text{if } x \le 0 \\ \kappa & \text{if } x \ge \kappa \\ x & \text{otherwise.} \end{cases} \quad (14) $$

Finally, the projections onto hyperplanes $C_2$ and $C_3$ are given explicitly as follows. For every $\mathbf{x} \in \mathbb{R}^{M \times N \times T}$ and $(m, n)$, we denote $(\Pi_{C_2}(\mathbf{x}))_{m,n}$ and $\Pi_{C_2}(\mathbf{x})_{m,n}$ the vectors that contain the projection coefficients of $\mathbf{x}$ onto $C_2$ at point $(m, n)$ in space at every time step. Then,

$$ 1) \quad (\Pi_{C_2}(\mathbf{x}))_{m,n} = \mathbf{x}_{m,n} - (1/\|L_1\|^2)L_1\mathbf{x}_{m,n}L_1^T $$

$$ 2) \quad (\Pi_{C_3}(\mathbf{x}))_{m,n} = \mathbf{x}_{m,n} - (1/\|L_2\|^2)L_2\mathbf{x}_{m,n}L_2^T $$

where

$$ L_1 = [(T - 3), -(T - 2), 0, \dots, 0, 1, 0]^\top \in \mathbb{R}^T $$
<table><caption>TABLE I
SATELLITE OPTICAL IMAGES USED IN THIS PAPER</caption><thead><tr><th>Satellites</th><th>Dates</th><th>Resolution (m)</th></tr></thead><tbody><tr><td>SPOT 1</td><td>15/06/1986 - 10/06/2001</td><td>10</td></tr><tr><td>SPOT 2</td><td>08/07/1990 - 14/09/1990 - 03/07/1998 - 07/10/2001</td><td>10</td></tr><tr><td>SPOT 3</td><td>20/10/1996</td><td>10</td></tr><tr><td>SPOT 4</td><td>16/11/2000 - 16/04/2001 - 16/05/2002</td><td>10</td></tr><tr><td>SPOT 5</td><td>21/07/2003 - 05/05/2010</td><td>2.5 - 5</td></tr><tr><td>Pléiades</td><td>21/03/2013 - 15/04/2013 - 29/07/2013 - 21/04/2014</td><td>0.7</td></tr></tbody></table>

and

$$L_2 = [0, 1, 0, \dots, 0, -(T-2), (T-3)]^\top \in \mathbb{R}^T.$$

## IV. METHOD VALIDATION

### A. Study Area

We apply the proposed method to detect the slow-moving landslides on the Colca Valley, a volcanic area in southern Peru. This area is affected by several slow-moving landslides of different sizes detected on the field and by satellite imagery [10], some of them being active for at least 30 years [12].

We select 16 satellite images from SPOT1-5 and Pléiades satellites covering the time period 1986–2013. Description of the images is given in Table I. The data set was first orthorectified using a DEM at 2-m resolution obtained from Pléiades imagery [10] and resampled to the same resolution of 10-m pixel size, before being processed using a time-series analysis to get a time series of displacement fields. The main steps of this preprocessing step are described in [12]. It includes: 1) the correlation of all possible pairs of images using the COSI-CORR correlator [13]; 2) the removal of possible bias and stripes in the displacement maps [14]; and 3) the inversion of the redundant system per pixel to provide a robust and unbiased time series of displacement [12]. This process leads to a series of ground-displacement maps for the 16 time steps with errors on displacement between 0.7 and 2.6 m [12], which is about one-seventh of a pixel size. To simplify our problem, and as explained earlier (see Section II), we project the displacement vector along the local topographical slope direction.

This data set of size 800 × 1200 × 16 is denoted $D_R$ in the following. A synthetic data set is created from $D_R$ that will be used to choose the regularization parameters $\lambda$ and $\tau$, and estimate the performance of the method.

### B. Synthetic Data Set

We build a simulated data set on a small area of size 100 × 100 pixels (1 × 1 km²). We first extract from $D_R$ an area where active landslide boundaries are known from the field observations [10], [39]. An image of the landslide activity only is created from the last image of the time series with displacement kept inside the boundary region of the active landslide, and set to 0 elsewhere.

We use this image as the last image $\mathbf{x}^T$ with $T = 16$ of the synthetic landslides $\mathbf{x}$. We then create a monotonic landslide

evolution, defining the remaining $\mathbf{x}^t$ with $t \in \{1, \dots, T-1\}$ of $\mathbf{x}$ as follows:

$$\forall (m, n) \in \mathcal{M} \times \mathcal{N}, \quad \mathbf{x}_{m,n}^t = \left(\frac{t}{T}\right)^2 \mathbf{x}_{m,n}^T. \quad (15)$$

The input data $\mathbf{z}$ is then created by adding to this synthetic landslide signal a cropped portion of the real data in an area without landslide signals. To improve the representativity of the proposed method, we crop three different areas that provide three different noise levels. Fig. 2 presents the last image ($\mathbf{x}^T$) of the input simulated data for the high-noise case and the results for different choices of the regularization functions ($\varphi$ and $\rho$). The zones indicated with red ellipses in Fig. 2(A) correspond to the landslide input signal.

### C. Results on the Synthetic Data Set

The algorithm is launched on the synthetic data set, using 10000 iterations at most. It can stop earlier at iteration $k$ if $\|\hat{\mathbf{x}}_k - \hat{\mathbf{x}}_{k-1}\| \le \sqrt{MNT} \times 10^{-6}$. The regularization parameters $\lambda$ and $\tau$ are adjusted to minimize the detection error between the original ($\bar{\mathbf{x}}$) and the reconstructed signals of landslide ($\hat{\mathbf{x}}$). The size $P$ of the patches of the $W$ operator is equal to 5.

We compare qualitatively the results for the different choices of norm for the regularization terms, the $\ell_2$-norm, $\ell_1$-norm, and the $\ell_{1,2}$ mixed-norm. We observe that the results obtained by using the $\ell_2$-norm do not detect the landslides [Fig. 2(B)], because the data set is too noisy. The use of the $\ell_1$-norm gives satisfactory results with low level of remaining noise [Fig. 2(C)]. However, the $\ell_{1,2}$ mixed-norm shows even better results [Fig. 2(D)]. This first analysis justifies the choice of the $\ell_{1,2}$ mixed-norm in a qualitative way.

We now compare quantitatively the results for these choices of the regularization terms. Since the landslide signals are stronger in the last image of the displacement time series, we analyze our results only on this last time step. This last image is transformed on a binary image where value 1 represents landslides if the absolute value of the pixel is greater than or equal to 0.5 m, and 0 otherwise (no-landslides). The images produced using the three different norms were evaluated with the reference landslide data set. The efficiency of the detection is estimated through different classical indicators that are described as follows.

1) Producer's accuracy: the number of correctly classified pixels of each category ("landslide," "no landslide") divided by the total number of pixels detected in each category. It measures the error of commission.
<figure><img src="image_3.png" /><figcaption>Fig. 2. (A) Input data (z) with high noise level. Output (x̂) by using (B) ℓ₂-norm, (C) ℓ₁-norm, and (D) ℓ₁,₂ mixed-norm for regularization.</figcaption></figure>

<table><caption>TABLE II
COMPARISON OF THE RESULTS OBTAINED BY USING THE ℓ₂-NORM, ℓ₁-NORM, AND ℓ₁,₂ MIXED-NORM FOR THE REGULARIZATION TERMS BASED ON φ AND ρ FUNCTIONS WITH (‖x̂T - zT‖₁/MN) = 0.2407</caption><thead><tr><th rowspan="2">Reference data</th><th rowspan="2"></th><th colspan="2">ℓ₂-Regularization</th><th colspan="2">ℓ₁-Regularization</th><th colspan="2">ℓ₁,₂-Regularization</th></tr><tr><th>landslides</th><th>no-landslides</th><th>landslides</th><th>no-landslides</th><th>landslides</th><th>no-landslides</th></tr></thead><tbody><tr><td rowspan="2"></td><td>landslides (222)</td><td>222</td><td>0</td><td>199</td><td>23</td><td>213</td><td>9</td></tr><tr><td>no-landslides (9778)</td><td>1080</td><td>8698</td><td>6</td><td>9772</td><td>45</td><td>9733</td></tr><tr><td rowspan="2">Producer's accuracy</td><td></td><td>17.05</td><td>100*</td><td>97.07*</td><td>99.77</td><td>82.56</td><td>99.91</td></tr><tr><td>User's accuracy</td><td>100*</td><td>88.95</td><td>89.64</td><td>99.94*</td><td>95.95</td><td>99.54</td></tr><tr><td colspan="2">Average accuracy (producer's)</td><td></td><td>58.53</td><td></td><td><b>98.42</b></td><td></td><td>91.24</td></tr><tr><td colspan="2">Average accuracy (user's)</td><td></td><td>94.48</td><td></td><td>94.79</td><td></td><td><b>97.75</b></td></tr><tr><td colspan="2">Overall accuracy</td><td></td><td>89.2</td><td></td><td><b>99.71</b></td><td></td><td>99.46</td></tr><tr><td colspan="2">Combined accuracy (producer's)</td><td></td><td>73.87</td><td></td><td><b>99.07</b></td><td></td><td>95.35</td></tr><tr><td colspan="2">Combined accuracy (user's)</td><td></td><td>91.84</td><td></td><td>97.25</td><td></td><td>98.60</td></tr><tr><td colspan="2">Detection error: <sup>||x̂T - x̂T||<sub>1</sub></sup>/<sub>MN</sub></td><td></td><td>0.1014</td><td></td><td>0.0466</td><td></td><td><b>0.0248</b></td></tr></tbody></table>

2) User's accuracy: the number of correctly classified pixels of each category divided by the total number of real pixels of each category. It measures the error of omission.

3) Average accuracy: average of the accuracy of the two categories.

4) Overall accuracy: the total number of pixels correctly detected divided by the total number of pixels.

5) Combined accuracy: the average of the overall accuracy and average accuracy.

6) Detection error (‖x̂T - ˆx̂T‖₁/MN): this is the only accuracy indicator of these tables that takes into account the original pixel values. It, therefore, evaluates the displacement amplitude of the final image.

Those indicators are computed for three levels of noise: low, mild, and high noise levels. Results are indicated in Tables II–IV for the three levels of noise. To improve the reading, numbers with * or * indicate the best indexes (in the same line) for the landslides or nonlandslides categories, respectively. The numbers in bold correspond to the best indexes for the average/overall/combined accuracy and the detection error.

These tables emphasize that using the ℓ₂-norm enables us to find almost all the pixels related to the landslides but also adds many false alarms. Despite large producer's accuracy for “no landslide” category and large user's accuracy for “landslides” category (greater than 99% for both indexes), the indexes of
<table><caption>TABLE III
COMPARISON OF THE RESULTS OBTAINED BY USING THE $ℓ_2$-NORM, $ℓ_1$-NORM, AND $ℓ_{1,2}$ MIXED-NORM FOR THE REGULARIZATION TERMS BASED ON $φ$ AND $ρ$ FUNCTIONS WITH $(||\bar{x}^T - z^T||_1/MN) = 0.6246$</caption><thead><tr><th rowspan="2">Results obtained</th><th rowspan="2"></th><th colspan="2">ℓ<sub>2</sub>-Regularization</th><th colspan="2">ℓ<sub>1</sub>-Regularization</th><th colspan="2">ℓ<sub>1,2</sub>-Regularization</th></tr><tr><th>landslides</th><th>no-landslides</th><th>landslides</th><th>no-landslides</th><th>landslides</th><th>no-landslides</th></tr></thead><tbody><tr><td rowspan="2">Reference data</td><td>landslides (222)</td><td>221</td><td>1</td><td>183</td><td>39</td><td>206</td><td>16</td></tr><tr><td>no-landslides (9778)</td><td>2447</td><td>7331</td><td>0</td><td>9778</td><td>95</td><td>9683</td></tr><tr><td colspan="2">Producer's accuracy</td><td>8.28</td><td>99.99*</td><td>100.00*</td><td>99.60</td><td>68.44</td><td>99.84</td></tr><tr><td colspan="2">User's accuracy</td><td>99.55*</td><td>74.97</td><td>82.43</td><td>100*</td><td>92.79</td><td>99.03</td></tr><tr><td colspan="2">Average accuracy (producer's)</td><td colspan="2">54.14</td><td colspan="2"><b>99.80</b></td><td colspan="2">84.14</td></tr><tr><td colspan="2">Average accuracy (user's)</td><td colspan="2">87.26</td><td colspan="2">91.22</td><td colspan="2"><b>95.91</b></td></tr><tr><td colspan="2">Overall accuracy</td><td colspan="2">75.52</td><td colspan="2"><b>99.61</b></td><td colspan="2">98.89</td></tr><tr><td colspan="2">Combined accuracy (producer's)</td><td colspan="2">64.83</td><td colspan="2"><b>99.71</b></td><td colspan="2">91.52</td></tr><tr><td colspan="2">Combined accuracy (user's)</td><td colspan="2">81.39</td><td colspan="2">95.41</td><td colspan="2"><b>97.40</b></td></tr><tr><td colspan="2">Detection error: ||x̄<sup>T</sup> - ẋ̄<sup>T</sup>||<sub>1</sub>/MN</td><td colspan="2">0.3320</td><td colspan="2">0.1239</td><td colspan="2"><b>0.1225</b></td></tr></tbody></table>

<table><caption>TABLE IV
COMPARISON OF THE RESULTS OBTAINED BY USING THE $ℓ_2$-NORM, $ℓ_1$-NORM, AND $ℓ_{1,2}$ MIXED-NORM FOR THE REGULARIZATION TERMS BASED ON $φ$ AND $ρ$ FUNCTIONS WITH NOISE LEVEL $(||\bar{x}^T - z^T||_1/MN) = 1.2454$</caption><thead><tr><th rowspan="2">Results obtained</th><th rowspan="2"></th><th colspan="2">ℓ<sub>2</sub>-Regularization</th><th colspan="2">ℓ<sub>1</sub>-Regularization</th><th colspan="2">ℓ<sub>1,2</sub>-Regularization</th></tr><tr><th>landslides</th><th>no-landslides</th><th>landslides</th><th>no-landslides</th><th>landslides</th><th>no-landslides</th></tr></thead><tbody><tr><td rowspan="2">Reference data</td><td>landslides (222)</td><td>220</td><td>2</td><td>212</td><td>10</td><td>182</td><td>40</td></tr><tr><td>no-landslides (9778)</td><td>3725</td><td>6053</td><td>516</td><td>9262</td><td>94</td><td>9684</td></tr><tr><td colspan="2">Producer's accuracy</td><td>5.58</td><td>99.97*</td><td>29.12</td><td>99.89</td><td>65.94*</td><td>99.59</td></tr><tr><td colspan="2">User's accuracy</td><td>99.10</td><td>61.90</td><td>95.5</td><td>94.72</td><td>81.98</td><td>99.04*</td></tr><tr><td colspan="2">Average accuracy (producer's)</td><td colspan="2">52.78</td><td colspan="2">64.51</td><td colspan="2"><b>82.77</b></td></tr><tr><td colspan="2">Average accuracy (user's)</td><td colspan="2">80.50</td><td colspan="2"><b>95.11</b></td><td colspan="2">90.51</td></tr><tr><td colspan="2">Overall accuracy</td><td colspan="2">62.73</td><td colspan="2">94.74</td><td colspan="2"><b>98.66</b></td></tr><tr><td colspan="2">Combined accuracy (producer's)</td><td colspan="2">57.75</td><td colspan="2">79.63</td><td colspan="2"><b>90.72</b></td></tr><tr><td colspan="2">Combined accuracy (user's)</td><td colspan="2">71.62</td><td colspan="2"><b>94.93</b></td><td colspan="2">94.59</td></tr><tr><td colspan="2">Detection error: ||x̄<sup>T</sup> - ẋ̄<sup>T</sup>||<sub>1</sub>/MN</td><td colspan="2">0.8766</td><td colspan="2">0.7554</td><td colspan="2"><b>0.3895</b></td></tr></tbody></table>

average/overall/combined accuracy are always smaller to the ones obtained with the $ℓ_1$ and $ℓ_{1,2}$ regularizations norms, for all noise level cases. Notably, the average accuracies (producer's) of $ℓ_2$-norm are 58.53%, 54.14%, and 52.78% for three noise level cases (from low to high noise), while $ℓ_{1,2}$ always gives the indexes greater than 82%. Moreover, the detection errors of $ℓ_2$-norm are at least two times larger than the ones of $ℓ_{1,2}$.

Both $ℓ_1$-norm and $ℓ_{1,2}$ mixed-norm give satisfactory results for the cases with low and mild noise level. $ℓ_{1,2}$ detects slightly more landslide pixels than $ℓ_1$, the average accuracies (user's) of $ℓ_{1,2}$ are 97.75% for low noise and 95.91% for mild noise level, while the ones of $ℓ_1$ are 94.79% for low noise and 91.22% for mild noise level. However, $ℓ_{1,2}$ also adds slightly more false alarms than $ℓ_1$. Thus, the average accuracies (producer's) of $ℓ_1$ (98.42% and 99.80%) are greater than the ones of $ℓ_{1,2}$ (91.24% and 84.14%) for these two noise level cases. When the noise is high, the $ℓ_1$-norm creates many false alarms (the average accuracy (producer's) is 64.51%), while the $ℓ_{1,2}$ still provides good results with 82.77%. Regarding the energy estimation, the $ℓ_{1,2}$-norm always provides better results than other norms. For instance, the energy detection error is twice larger using the $ℓ_1$-norm than $ℓ_{1,2}$ for the high-noise case study. To summarize, the $ℓ_{1,2}$ has the smallest error for the landslides detection in terms of sparsity and displacement energy. Furthermore, it is robust against noise. Also, note that the undetected landslide pixels lie along the landslide contour, a result that is expected from the use of the W operator in (4).

### D. Results on the Real Data Set

In this section, we aim at validating the performance of the proposed method for the real data set. The parameters $λ$ and $τ$ are chosen equal to the optimized parameters from the synthetic data for the mild noise level case. The size $P$ of the patches of the W operator in this case is equal to 10.

First, to qualitatively evaluate the performance of the $ℓ_{1,2}$ mixed-norm on the real data set, we display the outlier signal ($\hat{y}$) as estimated by the different norms [Fig. 3(B)–(D)]. The input data [Fig. 3(A)], despite the high noise content, clearly shows the important landslides (indicated by the black arrows). We note that the outliers obtained by using the $ℓ_2$- or $ℓ_1$-norm contain many signals related to the landslides. On the contrary, the $ℓ_{1,2}$ mixed-norm allows to clearly separate the outliers from the landslide signals. Thus, we consider, in the following, only the results obtained with the $ℓ_{1,2}$ mixed-norm.

The proposed method detects several patches of displacement, but some of them are very small in size and/or in energy. For these reasons, we select only the contiguous patches greater than 20 pixels (that is 2000 m²). The proposed method allows us to discover the motion increasing with time (that is a displacement oriented in the direction of the local slope), but also the motion decreasing with time (opposite direction to the local slope). In landslides, this opposite direction motion can correspond to upper parts of tilted blocks along listric faults. We assume that within a landslide, pixels presenting this opposite-slope motion are much less numerous
<figure><img src="image_4.png" /><figcaption>Fig. 3. (A) Real input data (z) at the end of the time series (April 2014). Outlier output (ŷ) (April 2014) by using (B) $ℓ_2$-norm, (C) $ℓ_1$-norm, and (D) $ℓ_{1,2}$ mixed-norm as regularization terms.</figcaption></figure>

than the downward motion. Therefore, among these patches, we removed from our analysis the ones whose number of pixels that follow the slope direction is less than 20% of the total pixels of the entire patch. This simple postprocessing leads to 22 displacement patches greater than 2000 m² that may be slow-moving landslides.

To validate this detection, we investigate, one by one, each patch by comparing it with previous knowledge on the slow-moving landslides in the area [10], [12] (Fig. 4, red circles), or visual inspection on field conducted in May 2017 (blue circles), or a check using images from Google Earth when field investigation was missing. Based on this investigation, some patches were either confirmed or infirmed as landslides, or their status could be left unknown. Indeed, signs of slow motion are not always visible even on high-resolution images from Google Earth. In Fig. 4, we summarize our results. Almost all the active landslides detected by [10] coincide with our detected landslides, except for the landslides No. 23 and No. 24. Landslide No. 23 corresponds to a landslide small in size (50 × 50 m²), easily detected on VHR images of 70-cm pixel size [10], but not on the 10-m resolution data set. Indeed, the image correlation steps use windows of 32 pixel size [12] which is much larger than the landslide size. Landslide No. 24 corresponds to a landslide of lower velocity than the other ones detected [10]. Uncertainties on the

displacement fields of about one-seventh of a pixel size [37] make the detection fairly easy with 70-cm pixel size images, but much more tricky with 10-m resolution images.

The field investigations conducted in May 2017 confirmed four of the newly detected landslides and infirmed two of them. The use of Google Earth images allows to confirm two other landslides and infirmed three of them. Interestingly, these false alarms correspond to real surface cover changes, as rapid vegetation growth, new roads, or new reservoir lakes. We denote these false alarms by the red circles in Fig. 4. However, for five cases (black circles), the detected landslide status is left unknown and further field investigations should be conducted. In total, over the 22 landslides detected, 12 were confirmed, 5 infirmed, and 5 have an unknown status. This detection presents commission and omission errors for landslides of 29% and 14%, respectively. Compared to the previous study in this area [10], six new slow-moving landslides were thus detected (+50%).

In Figs. 5 and 6, we show two examples of landslides detected by our method and confirmed by field investigations. First, the Maca landslide (No. 9) is considered as the biggest of all landslides in the Colca Valley [39]. As displayed in the bottom of Fig. 5, this landslide extends approximately 1 km North to South, and about 1.5 km wide East to West. Its motion is more rapid in the early 1990's due to a local
<figure><img src="image_5.png" /><figcaption>Fig. 4. Map of the studied area. Color scales with the elevation. The red rectangle contours the footprint of the processed satellite images. The blue curve corresponds to the Colca river. Symbols represent active slow-moving landslides detected in [10] with VHR images (green circles), new landslides verified by field investigations conducted in May 2017 and by using images from Google Earth (blue circles), landslides detected by the proposed method (yellow circles).</figcaption></figure>

earthquake that accelerated it [12]. During 2003–2010 period, it is fairly stable, which is confirmed by the GPS data [39]. Since 2011, it accelerated again due to intense precipitations (Fig. 5). The spatial pattern of motion reflects a complex slide, with a general motion toward the north, faster closer to the river than in the upper part, and imbricated blocks where the displacement direction is opposite to the local slope topography (blue patches). Those patterns reflect a typical mechanism of tilted blocks along listric faults [39]. These observations show that one should neither threshold the topographical slope nor select the motion direction to detect slow-moving landslides.

Second, the landslide No. 16 shown in Fig. 6 is a landslide of medium size (150 × 100 m²), situated in a meander of the Colca river. Clear scarps delimit this landslide, showing its high degree of maturity. Its time evolution shows an activation in the years 1998–2000, and high activity between 2010 and 2013, when intense precipitations occurred [39].

All these observations highlight the role of river erosion in its development.

This analysis and comparison of the results with VHR displacement maps show that the ability of the method to detect landslides is controlled by the input data set and, particularly, its resolution. Obviously, the image resolution controls the landslide size and velocity can be detected. The minimum landslide velocity can be detected depends on both the image resolution (that controls the displacement field uncertainty) and the time-separation between images. Fig. 6 shows that landslides with velocity as low as 0.2 m/year can be detected by our method using 10-m resolution images, taking the advantage of a 30-year long data set. Finally, the performance of the method is also dependent on the noise of the input time series of displacement with respect to the searched landslide signals. Uncertainties of one-fifth of a pixel size have been evaluated for many sensors of different resolutions [10]. This shows the importance to preprocess
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5. Landslide No. 9 (Maca) detected in the Colca Valley. (A) Detected landslide. Time series of displacement averaged over all the pixels of the landslide that (B) follow the slope direction and (C) reverse to the slope direction. (D) Picture of the corresponding area.</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 6. Landslide No. 16 detected in the Colca Valley. (A) Detected landslide. Time series of displacement averaged over all the pixels of the landslide that (B) follow the slope direction and (C) reverse to the slope direction. (D) Picture of the corresponding area.</figcaption></figure>

the data to improve the measurement uncertainty, as done in [12]. This also suggests that in theory, the results of the developed method are quite sensor dependent for sensors of the same resolution. Unfortunately, this assumption could not be

evaluated with our data set, as the Pléiades and SPOT images do not span the same observation period.

V. CONCLUSION

In this paper, we propose a new approach in optimization for the slow-moving landslide detection problem. We formulate
with mathematical terms the different properties of the signal including the monotony of the landslide motion. Furthermore, in this approach, an outlier estimation step is included to remove them from the solution. The physical properties of the signals are processed all together, combined in space as well as in time by using the $ℓ_{1,2}$ mixed-norm. The proposed approach, formulated in 3-D, overcomes the pixel-based analysis and limits the number of parameters. It depends only on two regularization parameters, $\lambda$ and $\beta$, and we propose a method to estimate them based on synthetic data generated directly from the real data.

The proposed approach demonstrates good detection ability, despite the high-noise-level conditions, both for simulated and real field data of the Colca Valley. In particular, it allowed detecting +50% of new slow-moving landslides, mostly small ones, compared to previous knowledge in the area, despite the medium resolution data set (10 m) and its long duration (28 years), where lots of noise, therefore, arise from ground cover changes. This data set is, therefore, a perfect proof of concept for its use on most recent open-source and global medium-resolution satellite data, as Sentinel-2 or Landsat-7/8. This method can also be tuned to process displacement fields from InSAR by taking into account a different structuring noise of radar measurements. This opens interesting perspectives, as time series of ground motion from InSAR from short revisited time satellites are able to measure small accelerating patterns of landslides before failure [40], [41].

## REFERENCES

[1] D. Petley, “Global patterns of loss of life from landslides,” *Geology*, vol. 40, no. 10, pp. 927–930, 2012.

[2] A. L. Handwerger, J. J. Roering, and D. A. Schmidt, “Controls on the seasonal deformation of slow-moving landslides,” *Earth Planet. Sci. Lett.*, vols. 377–378, pp. 239–247, Sep. 2013.

[3] R.-T. Rockafellar, “Monotone operators and the proximal point algorithm,” *SIAM J. Control Optim.*, vol. 14, no. 5, pp. 877–898, Aug. 1976.

[4] H.-H. Bauschke and P.-L. Combettes, “Convex Analysis and Monotone Operator Theory in Hilbert Spaces,” (CMS Books in Mathematics), Springer, 2011.

[5] P. L. Combettes and J.-C. Pesquet, “Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallel-sum type monotone operators,” *Set-Valued Variational Anal.*, vol. 20, no. 2, pp. 307–330, Jun. 2012.

[6] M. V. Afonso, J. M. Biucas-Dias, and M. A. T. Figueiredo, “An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems,” *IEEE Trans. Image Process.*, vol. 20, no. 3, pp. 681–695, Mar. 2011.

[7] J.-C. Pesquet and N. Pustelnik, “A parallel inertial proximal optimization method,” *Pacific J. Opt.*, vol. 8, no. 2, pp. 273–305, Apr. 2012.

[8] P. L. Combettes and J.-C. Pesquet, “Proximal splitting methods in signal processing,” in *Fixed-Point Algorithms for Inverse Problems in Science and Engineering*, H. H. Bauschke, R. Burachik, P. L. Combettes, V. Elser, D. Luke, and H. Wolkowicz, Eds. New York, NY, USA: Springer-Verlag, 2011, pp. 185–212.

[9] D. Cruden and D.-L. Varnes, “Landslide types and processes,” in *Landslides: Investigation and Mitigation*. Washington, DC, USA: National Academy Press, 1996, pp. 36–75.

[10] P. Lacroix, E. Berthier, and E. T. Maquerhua, “Earthquake-driven acceleration of slow-moving landslides in the Colca valley, Peru, detected from Pléiades images,” *Remote Sens. Environ.*, vol. 165, pp. 148–158, Aug. 2015, doi: 10.1016/j.rse.2015.05.010.

[11] P. Lacroix, G. Bièvre, E. Pathier, U. Kniess, and D. Jongmans, “Use of Sentinel-2 images for the detection of precursory motions before landslide failures,” *Remote Sens. Environ.*, vol. 215, pp. 507–516, Sep. 2018.

[12] N. Bontemps, P. Lacroix, and M.-P. Doin, “Inversion of deformation fields time-series from optical images, and application to the long term kinematics of slow-moving landslides in Peru,” *Remote Sens. Environ.*, vol. 210, pp. 144–158, Jun. 2018.

[13] S. Leprince, S. Barbot, F. Ayoub, and J.-P. Avouac, “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements,” *IEEE Trans. Geosci. Remote Sens.*, vol. 45, no. 6, pp. 1529–1558, Jun. 2007.

[14] S. Leprince, P. Mus, and J.-P. Avouac, “In-flight CCD distortion calibration for pushbroom satellites based on subpixel correlation,” *IEEE Trans. Geosci. Remote Sens.*, vol. 46, no. 9, pp. 2675–2683, Sep. 2008.

[15] A. Dehecq, N. Gourmelen, and E. Trouve, “Deriving large-scale glacier velocities from a complete satellite archive: Application to the Pamir-Karakoram-Himalaya,” *Remote Sens. Environ.*, vol. 162, pp. 55–66, Jun. 2015.

[16] C. Delacourt, P. Allemand, B. Casson, and H. Vadon, “Velocity field of the ‘La Clapière’ landslide measured by the correlation of aerial and QuickBird satellite images,” *Geophys. Res. Lett.*, vol. 31, no. 15, 2004.

[17] M. Fahnestock, T. Scambos, T. Moon, A. Gardner, T. Haran, and M. Klinger, “Rapid large-area mapping of ice flow using Landsat 8,” *Remote Sens. Environ.*, vol. 185, pp. 84–94, Nov. 2016, doi: 10.1016/j.rse.2015.11.023.

[18] A. Federico, M. Popescu, G. Elia, C. Fidelibus, G. Internò, and A. Murianni, “Prediction of time to slope failure: A general framework,” *Environ. Earth Sci.*, vol. 66, no. 1, pp. 245–256, doi: 10.1007/s12665-011-1231-5.

[19] A.-J. Hendron and F.-D. Patton, “The Vaiont slide, a geotechnical analysis based on how geologic observations of the failure surface,” U.S. Army Corps Eng., Tech. Rep. GL-85-5, Jun. 1985, vol. 1, p. 104.

[20] G. E. Hilley, R. Bürgmann, A. Ferretti, F. Novali, and F. Rocca, “Dynamics of slow-moving landslides from permanent scatterer analysis,” *Science*, vol. 304, no. 5679, pp. 1952–1955, 2004, doi: 10.1126/science.1098821.

[21] D. Hölbling, P. Füreder, F. Antolini, F. Cigna, N. Casagli, and S. Lang, “A semi-automated object-based approach for landslide detection validated by persistent scatterer interferometry measures and landslide inventories,” *Remote Sens.*, vol. 4, no. 5, pp. 1310–1336, 2012, doi: 10.3390/rs4051310.

[22] A. Kääb, “Monitoring high-mountain terrain deformation from repeated air- and spaceborne optical data: Examples using digital aerial imagery and ASTER data,” *ISPRS J. Photogramm. Remote Sens.*, vol. 57, pp. 39–52, Nov. 2002.

[23] Q. Sun, J. Hu, L. Zhang, and X. L. Ding, “Towards slow-moving landslide monitoring by integrating multi-sensor InSAR time series datasets: The Zhouqu case study, China,” *Remote Sens.*, vol. 8, no. 11, p. 908, 2016, doi: 10.3390/rs8110908.2016.

[24] P. L. Rosin and J. Hervás, “Remote sensing image thresholding methods for determining landslide activity,” *Int. J. Remote Sens.*, vol. 26, no. 6, pp. 1075–1092, 2005.

[25] A. Repetti, M. Q. Pham, L. Duval, E. Chouzenoux, and J.-C. Pesquet, “Euclid in a taxicab: Sparse blind deconvolution with smoothed $ℓ_1/ℓ_2$ regularization,” *IEEE Signal Process. Lett.*, vol. 22, no. 5, pp. 539–543, May 2015.

[26] M. Q. Pham, L. Duval, C. Chaux, and J.-C. Pesquet, “A primal-dual proximal algorithm for sparse template-based adaptive filtering: Application to seismic multiple removal,” *IEEE Trans. Signal Process.*, vol. 62, no. 16, pp. 4256–4269, Aug. 2014.

[27] M. Q. Pham, B. Oudompheng, J.-I. Mars, and B. Nicolas, “A noise-robust method with smoothed $ℓ_1/ℓ_2$ regularization for sparse moving-source mapping,” *Signal Process.*, vol. 135, pp. 96–106, Jun. 2016.

[28] E. A. Robinson and S. Treitel, *Geophysical Signal Analysis*. Englewood Cliffs, NJ, USA: Prentice-Hall, 1980.

[29] B. Pradhan, M. N. Jebur, H. Z. M. Shafri, and M. S. Tehrany, “Data fusion technique using wavelet transform and Taguchi methods for automatic landslide detection from airborne laser scanning data and quickbird satellite imagery,” *IEEE Trans. Geosci. Remote Sens.*, vol. 54, no. 3, pp. 1610–1622, Mar. 2016, doi: 10.1109/TGRS.2015.2484325.

[30] M. Kowalski, K. Siedenburg, and M. Dorfler, “Social sparsity! Neighborhood systems enrich structured shrinkage operators,” *IEEE Trans. Signal Process.*, vol. 61, no. 10, pp. 2498–2511, May 2013.

[31] M. Kowalski, “Sparse regression using mixed norms,” *Appl. Comput. Harmon. Anal.*, vol. 27, no. 3, pp. 303–324, 2009.

[32] M. Kowalski and B. Torrésani, “Sparsity and persistence: Mixed norms provide simple signal models with dependent coefficients,” *Signal, Image Video Process.*, vol. 3, no. 3, pp. 251–264, Sep. 2009.
[33] K. Chen, *Matrix Preconditioning Techniques and Applications* (Cambridge Monographs on Applied and Computational Mathematics), vol. 19. Cambridge, U.K.: Cambridge Univ. Press, 2005.

[34] J. Nocedal and S. J. Wright, *Numerical Optimization*. New York, NY, USA: Springer-Verlag, 1999.

[35] T. Pock and A. Chambolle, "Diagonal preconditioning for first order primal-dual algorithms in convex optimization," in *Proc. IEEE Int. Conf. Comput. Vis. (ICCV)*, Barcelona, Spain, Nov. 2011, pp. 1762–1769.

[36] A. Repetti, E. Chouzenoux, and J.-C. Pesquet, “A penalized weighted least squares approach for restoring data corrupted with signal-dependent noise,” in *Proc. 20th Eur. Signal Process. Conf. (EUSIPCO)*, Bucharest, Romania, Aug. 2012, pp. 1553–1557.

[37] A. Stumpf, J.-P. Malet, P. Allemand, and P. Ulrich, “Surface reconstruction and landslide displacement measurements with Pléiades satellite images,” *ISPRS J. Photogramm. Remote Sens.*, vol. 95, pp. 1–12, Sep. 2014.

[38] A. Stumpf, J.-P. Malet, and C. Delacourt, “Correlation of satellite image time-series for the detection and monitoring of slow-moving landslides,” *Remote Sens. Environ.*, vol. 189, pp. 40–55, Feb. 2017.

[39] S. Zerathe et al., "Morphology, structure and kinematics of a rainfall controlled slow-moving Andean landslide, Peru," *Earth Surf. Processes Landforms*, vol. 41, no. 11, pp. 1477–1493, 2016.

[40] E. Intrieri et al., "The Maoxian landslide as seen from space: Detecting precursors of failure with Sentinel-1 data," *Landslides*, vol. 15, no. 1, pp. 123–133, 2017.

[41] J. Dong et al., "Measuring precursory movements of the recent Xinmo landslide in Mao County, China with Sentinel-1 and ALOS-2 PALSAR-2 datasets," *Landslides*, vol. 15, no. 1, pp. 135–144, 2017.

<figure><img src="image_14.png" /></figure>

**Pascal Lacroix** received the M.S. degree in remote sensing from the University of Paris VI, Paris, France, in 2004, the Engineer degree in geophysics from the University of Strasbourg, Strasbourg, France, in 2004, and the Ph.D. degree in earth surface observation from LEGOS, University of Toulouse, Toulouse, France, in 2007.

Lima, Peru, where he was involved in the use of satellite images for ice-caps and landslide studies. Since 2013, he has been a Researcher with the Institut de Recherche pour le Développement, Marseille, France. He is involved in the development of geodetic and seismologic activities to monitor landslides in Peru, France, and Nepal. His research interests include the use of optical remote sensing to detect and monitor landslides, and understand the mechanics of landslides under seismic forcings.

He was a Post-Doctoral Researcher with ISTerre Laboratory, Gières, France, with the University of Tasmania, Hobart, TAS, Australia, with the University of Oslo, Oslo, Norway, and with INGEMMET,

<figure><img src="image_15.png" /></figure>

**Mai Quyen Pham** (S'14–M'15) received the B.S. degree in mathematics from the Hanoi National University of Education, Hanoi, Vietnam, in 2007, the M.S. degree in software engineering and statistics from Aix Marseille University, Marseille, France, in 2011, and the Ph.D. degree in signal and image processing from the University of Paris-Est, Paris, France, and the IFP Energies Nouvelles, Paris, in 2015.

From 2015 to 2017, she was a Post-Doctoral Fellow with the GIPSA-Lab, Saint-Martin-d'Hères, France, and also with the ISTerre Laboratory, Gières, France. In 2018, she joined the SIROCCO Research Team, INRIA Rennes-Bretagne Atlantique, Rennes, France, as a Post-Doctoral Fellow. Her research interests include inverse problems, particularly image restoration and reconstruction, seismic data processing, acoustic data, satellite imagery, and aspects of sparse representations, penalized criteria, and optimization methods.

<figure><img src="image_16.png" /></figure>

**Marie Pierre Doin** received the M.S. degree in earth sciences from École Normale Supérieure, Paris, France, in 1992, and the Ph.D. degree in geophysics from the University Paris VI, Paris, in 1995.

From 1996 to 2011, she was a Researcher with the Laboratoire de Géologie, Centre National de la Recherche Scientifique, École Normale Supérieure, Paris. She is currently Research Director with the ISTerre Laboratory, Gières, France. Her research interests include the tectonic deformation and seismic cycle, the thermomechanical interaction between lithosphere and asthenosphere, and the deformation mechanisms of the lithosphere. She is currently involved in interferometric synthetic aperture radar time series to monitor transient ground deformations of small amplitude, particularly around active faults and anthropogenic surface loads, in order to bring constraints on the crust and uppermost mantle rheology.