*Article* **Displacement Characterization and Spatial-Temporal Evolution of the 2020 Aniangzhai Landslide in Danba County Using Time-Series InSAR and Multi-Temporal Optical Dataset**

**Jianming Kuang 1,2, Alex Hay-Man Ng 1,2,3,\* and Linlin Ge <sup>2</sup>**


**Abstract:** On 17 June 2020, a large ancient landslide over the Aniangzhai (ANZ) slope, Danba County, Sichuan Province, China, was reactivated by a series of multiple phenomena, including debris flow triggered by heavy rainfall and flooding. In this study, Synthetic Aperture Radar (SAR) images acquired by the Sentinel-1A/B satellite and optical images captured by the PlanetScope satellites were jointly used to analyze and explore the deformation characteristics and the Spatial-Temporal evolution of the ANZ landslide before and after the multi-hazard chain. Several areas of pre-failure movements were found from the multi-temporal optical images analysis before the reactivation of the ANZ landslide. The large post-failure surface deformation over the ANZ slope was also retrieved by the optical pixel offset tracking (POT) technique. A major northwest movement with the maximum horizontal deformation of up to 14.4 m was found. A time-series InSAR technique was applied to analyze the descending and ascending Sentinel-1A/B datasets spanning from March 2018 to July 2020, showing that the maximum magnitudes of the Line of Sight (LoS) displacement velocities were −70 mm/year and 45 mm/year, respectively. The Spatial-Temporal evolution over the ANZ landslide was analyzed based on the time-series results. No obvious change in acceleration (precursory deformation) was detected before the multi-hazard chain, while clear accelerated deformation can be observed over the slope after the event. This suggested that heavy rainfall was the most significant triggering factor for the generation and reactivation of the ANZ landslide. Other preparatory factors, including the deformation behavior, the undercutting and erosion of the river and the outburst flood, the local terrain conditions, and earthquakes, might also have played an important role in the generation and reactivation of the landslide.

**Keywords:** Aniangzhai landslide; Danba County; multi-hazard chain; time-series InSAR; optical image analysis

#### **1. Introduction**

The movement of a wide range of ground elements, such as rock masses, soil, debris, or garbage, is referred to as a landslide [1]. Landslides are natural phenomena distributed all around the world; they pose a severe threat to lives and infrastructures and have caused enormous loss to our society in terms of safety and economy. A series of geological, hydrometeorological, or even human-related factors, such as seismic and volcanic activity [2,3], climate change [4], heavy rainfall [5–7], subsurface and surface engineering work [8,9], deforestation [10], agricultural activities, and urbanization [11], have significantly contributed to the occurrence of landslides. The southwest part of China is prone to the occurrence of landslides because of the abundant precipitation and mountainous topography. Sichuan

**Citation:** Kuang, J.; Ng, A.H.-M.; Ge, L. Displacement Characterization and Spatial-Temporal Evolution of the 2020 Aniangzhai Landslide in Danba County Using Time-Series InSAR and Multi-Temporal Optical Dataset. *Remote Sens.* **2022**, *14*, 68. https://doi.org/10.3390/rs14010068

Academic Editor: Francesca Cigna

Received: 19 October 2021 Accepted: 21 December 2021 Published: 24 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Province is the one of the regions in the southwest part of China which frequently suffers from medium- and large-scale landslides. On 10 July 2013, days of heavy rainfall and floods triggered a landslide in Sanxi village, Dujiangyan city, Sichuan Province, which caused 44 fatalities, 117 missing people, and the destruction of 12 houses [12,13]. The catastrophic 2017 Maoxian landslide in Xinmo village, Sichuan Province, buried 64 houses, killed 10 people, and left 73 people missing [6,14]. Landslide mapping and monitoring over these landslide-prone regions is essential for understanding their failure mechanisms and their long- and short-term evolutions. The measurement of the surface displacement caused by landslides is the most intuitive indication for understanding the failure mechanisms and evolution patterns over a long and short period. Conventional field investigation or in situ measurements can indeed provide a high accuracy of measurements at centimetre to millimetre levels [15–17]. However, the several limitations of these techniques are that: (1) only measurements on a point-by-point basis can be provided; (2) there is a high demand for the labour force, the techniques are time-consuming, and there are expensive costs for the installation and maintenance of the survey network; and (3) plenty of landslides in remote and rural areas are inaccessible.

The Spaceborne Interferometric Synthetic Aperture Radar (InSAR) technique, as a radar remote-sensing technique, can measure surface deformation along the radar Line of Sight (LoS) direction at high accuracy with all-weather and all-day operational capability, wide spatial coverage, and high spatial resolution (up to 1 m) [18–22]. Multi-temporal InSAR, also called time-series InSAR, is an improved InSAR technique. This technique can significantly minimize the temporal and spatial decorrelation issues in the conventional InSAR technique. A variety of time-series InSAR approaches were developed by different researchers, for example the Permanent Scatterer InSAR (PSInSARTM) [23], the Small Baseline Subset Approach (SBAS) [24,25], and the SqueeSARTM [26], which provide a chance for researchers to detect the subtle movement of slow-moving landslides and gain a retrospective view of the historic displacements in the temporal and spatial domain. Time-series InSAR has been widely and successfully applied in landslide investigation, including landslide detection [27–29], precursory deformation detection, and pre- and post-failure analysis [6,30,31], as well as landslide inventory, susceptibility, and hazard assessment [32].

For example, Intrieri et al. [30] have successfully applied the SqueeSARTM to detect the precursors of failure before the 24 June 2017 Maoxian landslide in China; Zhou et al. [33] used the Stanford Method for Persistent Scatterer (StaMPS) small baselines subset method with Sentinel-1 images to detect the obvious spatiotemporal deformation of the Muyubao landslide in the Three Gorges Reservoir area in China; Ciampalini et al. [34] have demonstrated that the PSInSAR and SqueeSAR™ techniques, with the COSMO-SkyMed (CSK) images, can be used to further refine the Landslide susceptibility maps (LSM) of the Messina Province (Sicily, Italy); Dong et al. [35] developed a time-series InSAR method, known as coherent scatterer InSAR (CSI), which can jointly analyse the persistent scatterers (PS) and distributed scatterers (DS), and applied it to an investigation of the deformation characteristics of the Jiaju landslide in Danba County, China, using the ALOS PALSAR and ENVISAT ASAR datasets.

On the other hand, the optical remote-sensing technique can be used for visual interpretation of the Spatial-Temporal evolution of the landslide boundary and the detection of precursory features such as scarps and cracks over the slope surface. Besides, the optical pixel offset tracking (POT) technique can be applied to retrieve the horizontal displacement caused by slope movement with an accuracy of up to about 0.2 pixels of the optical images used [36–38]. The relatively large displacement (up to tens of metres) of landslide is not likely to be seen by the time-series InSAR technique due to high decorrelations. More importantly, the time-series InSAR technique is not sensitive to the movement in the north– south direction due to the near-polar orbits of SAR satellites [39]. Horizontal displacements measured by the optical POT technique are therefore an excellent supplement for the time-series InSAR technique in this regard. Hence, the time-series InSAR technique and the

optical remote-sensing technique are complementary in their detection of Spatial-Temporal evolution and their displacement measurements of landslides at different stages.

On 17 June 2020, an intense, short-duration rainfall struck the region of the Aniangzhai and Guanzhou villages in Danba County, Sichuan Province, China, which triggered a serious multi-hazard chain along the Xiaojinchuan River and, most importantly, reactivated the ancient Aniangzhai (ANZ) landslide, posing a serious threat to 12 towns and 112 villages along the downstream of the Xiaojinchuan River, including over 20,000 people in 5000 families. Emergency field investigations were conducted shortly after the occurrence of the event, which provided valuable and reliable information for our study [40,41]. The multi-hazard chain was also interpreted and characterized with the use of seismic signal data [42]. However, these field survey results could not provide a deformation map with dense measurement points and only post-event information could be obtained. By using satellite remote-sensing data for monitoring the landslide, the characteristics of the historic deformation with dense measurement points could be reviewed, and unstable slopes could be identified in the early stage, which provided significant information for the early warning of landslides and the emergency response.

In this study, the time-series InSAR technique was applied to process two stacks of Sentinel-1A/B data (descending and ascending tracks) to map the surface displacement and explore the Spatial-Temporal evolution of the ANZ landslide. Besides that, a timeseries optical images analysis, based on multi-temporal, high-resolution optical images acquired from the PlanetScope satellite, was also conducted to investigate the evolution of the landslide features over the ANZ slope. The large post-failure displacement over the slope was retrieved using the optical POT technique. The displacement characteristics, failure mechanism, and triggering factors of the ANZ landslide were further explored and discussed, based on the radar and optical remote-sensing results.

#### **2. Study Area**

As shown in Figure 1a, the ANZ landslide is located in the central-west section of the Sichuan Province, which is in the southwest part of China, with the location at 30.979◦ N, 102.024◦ E. The region lies at the transition zone between the Qinghai-Tibetan Plateau and the Sichuan Basin where the elevation difference is up to 7000 m (Figure 1b). At the same time, this region is located at the interaction zone of two large, active fault belts: the Xianshuihe-Xiaojiang fault in the west and the Longmenshan thrust fault belt in the east; this is a seismically active region with plenty of active faults around and numerous geohazards occurring. According to the United States Geological Survey (USGS), there have been 13 large historical earthquakes with a magnitude of over 6.0 on the Richter Scale recorded in this region, including the 1933 MW 7.5 Diexi Earthquake on the Songpinggou fault [43], the 1973 MS 7.6 Luhuo earthquake on the Xianshuihe fault [44], and the 2008 MW 7.9 Wenchuan Earthquake on the Longmenshan fault [45], as marked by the yellow and blue stars in Figure 1b. The most recent event was recorded on 20 April 2013, with a distance of 110 km in the southeast and a magnitude of 6.6 (MW).

**Figure 1.** Geological setting of the study area: (**a**) overall location of the study area; (**b**) regional tectonic setting and topography map derived from the 1-arc-second Digital Elevation Model (DEM) of the Shuttle Radar Topography Mission (SRTM) [46]; (**c**) overview of the multi-hazard chain in the study area. Black lines indicate the contour lines.

Figure 1c depicts the detailed overview of the multi-hazard chain on 17 June 2020. The ANZ landslide is located on the left bank of the Xiaojinchuan River. The ANZ landslide is an ancient landslide that once blocked the Xiaojinchuan River, and it has experienced slow-creep deformation in recent years, according to local residents [40]. According to local records, at about 3:00 am on 17 June 2020, the intense rainfall at midnight had triggered the Meilong debris flow on the right bank of the Xiaojinchuan River. The hillslope failure with a high elevation brought the deposits of the Meilong debris flow rapidly rushing down to the Xiaojinchuan River. According to the field investigation [40], the riverbed was uplifted by approximately 8–12 m, which blocked the Xiaojinchuan River in the upstream and formed a large dammed lake just behind the Guanzhou power station. A new river channel was then formed over the deposits, and it flowed down to the toe of the ANZ landslide, which seriously washed and eroded the toe section. At about 11:00 p.m. local time on the same day, the dammed lake began to breach and caused a flood [41]. Under the continuous downcutting and erosion of the new river channel and outburst flood, the whole toe section of the ANZ slope collapsed on 18 June 2020. Then, the lower and middle sections of the ancient ANZ landslide were reactivated, and the affected area was named the reactivated ANZ landslide. That formed a steep face with a slope of about 60–70◦ and a height of 60 m at the toe of the ANZ landslide, according to the field investigation. In contrast, the slope in the middle and upper section was relatively gentle, ranging from 35◦ to 45◦. The elevation difference of the ancient ANZ landslide is up to 1000 m, with about 3100 m at the top. According to the rainfall data recorded from the Danba Meteorological Bureau, the annual rainfall in the area of the ANZ landslide is about 532.7–823.3 mm, 80% of which is mainly concentrated between May and September [40]. The average annual temperature in this area is 14.2 ◦C. According to the filed investigation and published work, the lithology of the ANZ is mainly composed of the Lower Carboniferous Yaoji Formation (D1q) and the Quaternary sedimentary layers [41]. The bedrock is primarily composed of limestone, with a generally dipping attitude of N50 − 60◦W/NE∠70 − 80◦ and a visible thickness of 5–20 m. The avalanche deposits are widely spread over the surface of the slope, with the thickness ranging from 10 to 50 m. Fifty to eighty percent of this unit is

rock fragments, which are mainly angular and sub-angular blocks from the detachment of the rock avalanche. The ANZ landslide is located on a northwest-facing slope, with an inclination of 50◦ at the top and 20◦ at the lower part.

#### **3. Materials and Methods**

*3.1. Datasets*

#### 3.1.1. Optical Images

The PlanetScope satellite can acquire optical images with high resolution (3 m in pixel size) in near daily global coverage, which can be used as a useful tool for landslide detection and monitoring [47]. In this study, Planetscope images with four bands (blue, green, red, and near infrared) were collected from the Planet Labs Company under a research and education license. The level 3B surface reflectance product of PlanetScope was chosen as it has been orthorectified and pre-processed with geometric, radiometric, and atmospheric corrections [48]. A total of nine scenes of Planetscope images from 26 March 2020 to 27 July 2020 were selected with cloud-free (i.e., 0% cloud percentage) and minimal shadow condition over the ANZ slope. The temporal coverage spans both the pre-failure and post-failure stages of the multi-hazard chain. The detailed information of the PlanetScope images is shown in Table 1.

**Table 1.** Optical images acquired by the PlanetScope satellite for this study.


#### 3.1.2. Satellite SAR Datasets

In this study, both the descending and the ascending tracks of the C-band Sentinel-1A/B SAR images in the Terrain Observation by Progressive Scans (TOPS) mode were collected for time-series InSAR analysis. The TOPS SAR data can cover a swath width of 250 km at about 5 m by 20 m resolution in the range and azimuth directions, respectively [49]. The basic parameters of the SAR images are shown in Table 2. Fifty-eight images were included in the descending stack spanning from 27 March 2018 to 2 July 2020 and seventy-one images in the ascending stack for the period from 8 March 2018 to 7 July 2020. The spatial coverages of both tracks are shown in Figure 1b, with black and purple rectangles for the descending and ascending tracks, respectively.

**Table 2.** Basic parameters of SAR data.


#### 3.1.3. Rainfall Data

The monthly rainfall spanning from 2018 to 2020 was calculated from the daily rainfall based on the CHIRPS (Climate Hazards Group InfraRed Precipitation with Station) data from the Google Earth Engine [50,51]. In addition, the daily rainfall data between 1 May 2020 and 17 June 2020, according to the Anianggouer station, were collected from the published paper [40], as shown in Figure 2a. The daily rainfall data recorded after the multi-hazard chain, from 22 June 2020 to 19 August 2020, were also collected from the locally installed rainfall gauge (Figure 2b) [41]. The locations of the Anianggouer station and locally installed rainfall gauge are shown in Figure 1c.

**Figure 2.** (**a**) Daily rainfall data between 1 May 2020 and 17 June 2020, before the multi-hazard chain, according to the Anianggouer weather station from the Danba Meteorological Bureau. The daily rainfall data were collected from Zhao et al., 2021 [40]. (**b**) Daily rainfall data between 22 June 2020 and 19 August 2020, after the multi-hazard chain, according to the installed rainfall gauge over the ANZ slope. These data were collected from Zhu et al., 2021 [41].

#### *3.2. Optical Analysis Method*

Three main procedures were conducted based on the collected optical images. Firstly, true color PlanetScope images were generated by compositing three bands (red, green, and blue) for the initial visual interpretation of the evolution of the surface features over the ANZ slope. Secondly, by using the vegetation index, the pre- and post-failure landslide features can be effectively mapped and detected, and the effect of the mountain shadows can be minimized. The Atmospherically Resistant Vegetation Index (*ARVI*) time series was derived and was then classified using unsupervised classification to investigate the Spatial-Temporal evolution of the landslide in the pre- and post-failure stages [52]. Finally, the optical pixel offset tracking (POT) technique was applied to obtain the post-failure displacement of the ANZ slope.

The interpretation of the landslide features based on true color images could often be influenced by the shadows in the images as the landslide generally occurred in the mountain region. Using the Normalized Difference Vegetation Index (NDVI) for the landslide interpretation could minimize the influence of the mountain shadows [53]. The *ARVI*, an improved index compared to the NDVI, uses blue light reflectance measurements to correct for the atmospheric scattering effects through a self-correction process. The *ARVI* calculation can be expressed as [54]:

$$\begin{array}{l}\text{ARVI} = (\rho\_{nir} - \rho\_{rb}) / (\rho\_{nir} + \rho\_{rb})\\\rho\_{rb} = \rho\_{red} - \gamma(\rho\_{blw} - \rho\_{red})\end{array} \tag{1}$$

where *ρnir*, *ρred*, and *ρnir* are the near-infrared, red, and blue bands, respectively. The assumption that *γ* = 1 is often used unless the aerosol model is known a priori.

In this study, the *ARVI* was calculated based on Equation (1) using the collected PlanetScope images. These images were divided into two groups: the pre- and post-failure groups. The *ARVI* images in the pre-failure and post-failure groups were stacked to form the time-series *ARVI* with 6 and 2 layers, respectively, in a chronological order using the Layer Stack module in the ENVI software. Then, the K-means unsupervised classification method was applied to classify the pre-failure time-series *ARVI* and the post-failure timeseries *ARVI* image. The number of classes was set to 5 with other the default parameters in the classification process.

COSI-Corr, a useful POT software for measuring sub-pixel ground surface deformation based on the image correlation method [55], was used to map the post-failure displacement over the ANZ slope. The PlanetScope image pair between 16 June 2020 (pre-failure) and 24 June 2020 (post-failure) was chosen as the input data, and the red band was used for the POT calculation. The combinations of the initial and final window sizes used in this study are 128 and 64, respectively, which can reduce the background noise and uncertainties [56]. Three layers of the POT results can be derived, including displacements in the east–west and north–south directions, and the signal-to-noise ratio (SNR). Positive values in the east– west and north–south displacement layers indicate eastward and northward movements, respectively. Displacements to the west and south are represented by negative values. The range of SNR values is from 0 to 1, and the higher value suggests the measurement with the higher quality. Measurements with an SNR value of less than 0.95 were masked in the east–west and north–south displacement datasets. To reduce the mismatch error for the image pair, a stable area was selected near the ANZ slope. The displacements derived from the stable area can be used to depress the background noise through removing the co-registration errors within the image pair used. In this case, the mean east–west and north–south displacement derived from stable area was used to correct the image mismatch and generate the final corrected east–west and north–south displacement. The standard deviations of the east–west and north–south displacements in the stable area were used justify the uncertainties of the results [56]. The overall horizontal displacement was then calculated based on the corrected east–west and north–south displacements.

#### *3.3. Time-Series InSAR Analysis Method*

The Persistent Scatterer (PS) and Small Baseline (SB) approaches were both exploited and jointly processed for the time-series analysis. The combination of PS and Coherent Scatterer (CS) points was applied to increase the density of the measurement points (MPs) in the rural and mountainous region. Three major processing steps are included: PS processing, SB processing, and combination processing of the PS and CS targets. The Sentinel-1 datasets were pre-processed using the TOPSAR stack processor module in the Interferometric synthetic aperture radar Scientific Computing Environment (ISCE) software to generate single-master and multiple-master differential interferograms stacks [57]. Timeseries analysis was then conducted using the Stanford Method for Persistent Scatterers (StaMPS) software based on the differential interferograms generated [58].

As for the PS processing, the PS candidates were initially selected based on the Amplitude Dispersion Index (ADI) [23]. The pixels with ADI values over 0.4 were initially selected as the PS candidates. Then, the final PS targets are further selected based on the phase stability criterion [59].

For the case of the SB processing, the Spatial-Temporal distribution of the interferometric combination of the Small Baseline network is shown in Figure 3. The perpendicular baseline for over 90% of pairs was less than 100 m due to the precise orbit control of the Sentinel-1A/B satellites. To ensure low temporal decorrelation, the maximum temporal baseline was set to 120 days. The CS points were selected based on the amplitude difference dispersion (*D*Δ*A*), which is similar to *DA* [58].

For both the PS and the SB processing, the topographic effects were removed using the 1-arc-second SRTM DEM [46]. The overall wrapped differential phase *ϕDif f* of a pixel cell after flattening and removing the topographic phase can be expressed as follows:

$$\varphi\_{Diff} = W \left\{ \varphi\_{Defo} + \varphi\_{Attmos} + \varphi\_{Orb} + \varphi\_{DEM} + \varphi\_N \right\} \tag{2}$$

where *W*{} denotes the wrapping operator. *ϕDef o* represents the Line of Sight (LoS) deformation phase term. *ϕAtmos* indicates the atmospheric delay phase and *ϕOrb* is the phase caused by orbit error. *ϕDEM* represents the residual topographic phase caused by the DEM error and *ϕ<sup>N</sup>* is the random noise phase, including thermal noise and co-registration error. The deformation phase can be estimated and separated from the other phase terms by iteratively filtering in the temporal–spatial domain. The detailed processing chain can be referred to in [58,60,61].

The displacement velocities measured from the times-series InSAR are along the radar Line-of-Sight (*LoS*) direction. Therefore, the displacement measured from the time-series InSAR is a composite of vertical, easting, and northing displacement components, which can be written as [62–66]:

$$
\begin{bmatrix}
\cos(\theta) & -\sin(\theta)\cos(\mathfrak{a}) & \sin(\theta)\sin(\mathfrak{a})
\end{bmatrix}
\begin{bmatrix}
D\_V\\D\_E\\D\_N
\end{bmatrix} = D\_{LoS} \tag{3}
$$

where *θ* is the incidence angle and *α* is the heading angle of the satellite (positive clockwise from the north). *DV*, *DE*, and *DN* are the vertical, eastward, and northward displacement, respectively. *DLoS* is the displacement measured in the radar *LoS* direction. Ideally, measurements from at least three viewing geometries are needed to resolve the 3D deformation.

#### **4. Results**

#### *4.1. Landslide Evolution and Deformation Mapped by Optical Datasets*

#### 4.1.1. Landslide Evolution Revealed by the Time-Series Optical Analysis

Visual interpretation based on multi-temporal optical images can provide useful information from identifying surface moving features, which could directly reflect the formation and evolution of landslides. As shown in Figure 4, true color images were synthesized with three bands: red (B3 at 590–670 nm), green (B2 at 500–590 nm), and blue (B1 at 455–515 nm). Before 17 June 2020, there were three flows on the pre-failure surfaces observed, with two small masses over the main slope and a long flow at the upper neck. After the multi-hazard chain occurred on 17 June 2020, the collapse area could be clearly observed at the toe of the landslide, and the Xiaojinchuan River extended much wider because of the flooding caused by the debris flow at the upstream.

**Figure 4.** Time-series true color images derived from multi-temporal optical images from the PlanetScope satellite over the ANZ Landslide, with acquisition dates annotated in the upper-right corners.

A long stripping flow was detected close to the left boundary of the ancient landslide area. At the same time, multiple dislocations could be clearly observed over the roads at the middle of the slope, which was also consistent with the published field investigation [40]. Two obvious head scarps at the middle were detected on the post-failure image acquired on 24 June 2020, and it possibly penetrated down to the toe through the multiple road dislocations at the upstream section. Those two head scarps were connected from the image acquired on 27 July 2020. No obvious changes were seen at the pre-failure flowing mass on the upper neck.

Time-series *ARVI* images were generated based on high-resolution PlanetScope images acquired from 26 March 2020 to 27 July 2020, as shown in Figure 5. The pre- and post-failure flows over the slope can be obviously detected over the slope due to the bright background and reduced influence of the mountain shadows. It is noted that the head scarps were clearly observed to run through and down to the toe of the slope in the post-failure *ARVI* image. For further investigation of the evolution of the landslide features over the slope, the K-mean unsupervised classification based on the time-series *ARVI* images was conducted. Two groups were divided, corresponding to two different stages: the pre-failure stage (10 May 2020–16 June 2020) and the post-failure (24 June 2020–27 July 2020) stage, as shown in Figure 6. The *ARVI* image on 30 May 2020 was discarded in the K-mean unsupervised classification due to limited ground coverage. Five major classes were identified over the slope. Class 1, with a red colour, consists of ground features with the lowest *ARVI* value, including river and bare soil land. Class 2, with a pink colour, represents features with a medium-low *ARVI* value, including road, landslide scarps, and buildings. The landslide area includes flows, bare soil land, and landslide scarps. Therefore, the flows could be jointly covered by class 1 and 2 in this case. Class 4 and 5 represent ground features with a high *ARVI* value, including evergreen, tall, and dense vegetation over the slope. Class 3 is the transition zone among classes 1–2 and classes 4–5, possibly representing unstable vegetation cover over this zone.

Landslide scarps and flows can be distinguished from rivers, roads, and buildings as they are irregular features and located on the same slope. As shown in Figure 6a, the Xiaojinchuan River and roads, as linear features, are detected in class 1 and class 2, respectively. At the same time, three pre-failure flows were clearly detected during the pre-failure stage, with two at the centre slope and one at the upper slope. As for the post-failure stage on 24 June 2020 and 27 July 2020 in Figure 6b, the Xiaojinchuan River was obviously much wider than before due to the occurrence of flooding caused by the debris flow on the upstream. The post-failure flow next to the left boundary at the downstream section was observed, with an area up to 22,000 m2. The pre-failure flow at the centre upstream was extended further when compared with that from the pre-failure stage, while no obvious change was observed for the pre-failure flow at the upper slope.

#### 4.1.2. Post-Failure Displacement Detected from Optical Pixel Offset Tracking

Figure 7 shows the optical POT results between one pre-failure image on 16 June 2020 and one post-failure image on 24 June 2020. Positive values in the east–west and north–south measurements indicated that the slope moved towards the west and north directions, respectively. Westward and southward movements are represented by negative values. The standard deviation of the east–west and north–south displacements in the stable are 0.43 m and 0.49 m, respectively. These minor standard deviations indicated that the selected stable area is reliable and can be used to derive the corrected east–west and north–south displacements [56]. As shown in Figure 7a,b, post-failure movement over the reactivated landslide area can be clearly observed from both the west–east and north–south directions. The horizontal displacement is generated by combining these two measurements in different directions (Figure 7c). The ANZ slope mainly moved towards the northwest direction, with the maximum west and north displacement up to −13.7 m and 10.0 m, respectively. The horizontal displacement peaks at 14.4 m, and most of the large displacement ranging from 12.0 to 15.0 m, were found at the middle and upstream section of the reactivated landslide area. However, less than 4 m of horizontal displacements were observed in the upper section of the ancient landslide area during this period. The overall reactivated boundary of the ANZ landslide can be clearly observed from the ancient landslide area, showing an hourglass-shaped boundary from the top to the bottom.

**Figure 5.** Time-series *ARVI* images derived from high-resolution optical images from the PlanetScope satellite over the ANZ landslide, with acquisition dates annotated in the upper-right corners.

**Figure 6.** K-mean unsupervised classification based on the time-series *ARVI* images: (**a**) pre-failure stage (10 May 2020–16 June 2020); (**b**) post-failure stage (24 June 2020–27 July 2020).

**Figure 7.** Post-failure deformation measured from image pair between 16 June 2020 and 24 June 24, acquired from the PlanetScope satellite: (**a**) west–east displacement; (**b**) north–south displacement; and (**c**) magnitude of horizontal displacement covered by the proportional black arrows indicating moving direction.

#### *4.2. Deformation Velocity Measured by Time-Series InSAR Analysis*

The LoS displacement annual-velocity map can be obtained from the time-series InSAR processing of the Sentinel-1A/B datasets from the descending and ascending orbits, as shown in Figure 8. The measurement points with a red colour (negative values) represent targets moving away from the satellite along the LoS direction, while those with a blue colour represent motion toward the satellite. Due to dense vegetation cover and steep topography over the region, only limited measurement points were obtained. According to the deformation characteristics and local topography, the ANZ landslide can be divided into four zones: the central deformation zone (I), the landslide toe collapse zone (II), the upstream deformation zone (III), and the landslide shoulder deformation zone (IV). Deformation caused by the slope movement can be clearly observed from the descending measurement at the central deformation zone (I), the landslide toe collapse zone (II), and the upstream deformation zone (III), mostly ranging from −20 to −50 mm/year, as shown in Figure 8a. The maximum magnitude of deformation velocity was up to −70 mm/year, which was found at landslide toe collapse zone (II) near the Xiaojinchuan River. On the other hand, the land surface was evolving less over the landslide shoulder deformation zone (IV), and the deformation velocity ranging from 15 to −30 mm/year was observed. As for the ascending measurement shown in Figure 8b, the overall detected deformation velocity is smaller than the descending measurement because of the lower sensitivity of the ascending SAR data for the northwest-facing slope [67]. The density of measurement points is also lower, relatively, than that from the descending orbit. The deformation rate of most measurement points was between 5 to 15 mm/year along the LoS direction. Similarly, the largest deformation was found at zone (II), with a rate of 45 mm/year along the LoS direction. In the upstream section of zone (II), where there was no measurement point detected from the descending SAR data, it was detected to be unstable, with most of the points ranging from 20–45 mm/year from the ascending data, which could probably be uplifted accumulation at the left bank of the Xiaojinchuan River. In addition, much smaller deformations were observed at zone IV, with the deformation rate ranging from −10 to 7 mm/year.

**Figure 8.** LoS deformation velocity detected from the Sentinel-1: (**a**) descending dataset and (**b**) ascending datasets. Black star represents the reference point, and the black circles represent the selected measurement points.

Beside the ANZ landslide, an unstable slope was detected 200 m away from the right side of the ANZ landslide, with its maximum deformation velocity up to 19 and 30 mm/year from the descending and ascending measurements, respectively. As both measurements showed positive movement over the detected unstable slope, such an area is likely to be an uplifted area. The location of this unstable slope is just above the dammed lake and next to the Guanzhou power station, which could form a secondary landslide and causes a serious threat to the station.

In many InSAR studies, the LoS displacement measured from time-series InSAR can be projected into the vertical direction by assuming that there was no displacement in horizontal direction [63,68], whereas in this case the horizontal displacement was obviously detected over the ANZ slope from the optical POT result (Figure 7) and could not be neglected. In order to conduct a cross-validation between the descending and ascending measurements, a stable area over the Guanzhou power station was selected. This power station is expected to be relatively stable (i.e., no deformation) during the InSAR acquisition period. Both the ascending and the descending measurements were first projected into the vertical direction by dividing the cosine of the local incidence angle, and the differences between the projected vertical displacements from both tracks were calculated. The standard deviation of the calculated differences is 5.5 mm/year. As shown in Figure 9, the differences were mainly distributed in the range of −5 mm/year to 5 mm/year, accounting for over 60% of the measurement points, indicating that the two measurements are highly correlated with each other and that the time-series InSAR results are reliable.

**Figure 9.** Cross validation based on the differences between the projected vertical displacements from the descending and ascending InSAR measurements.

#### **5. Discussion**

#### *5.1. Temporal–Spatial Evolution of the ANZ Landslide*

To further investigate the temporal evolution of the pre-failure slope stability and whether there were any precursory signals detected before the multi-hazard chain, a time-series analysis of selected measurement points was conducted (Figure 10). Several measurement points at four deformation zones of the ANZ landslide were selected, as marked by the black circles in Figure 8. In order to investigate the effect of local rainfall on the ANZ landslide, the monthly rainfall was calculated from the daily rainfall based on the CHIRPS (Climate Hazards Group InfraRed Precipitation with Station) data from the Google Earth Engine [50,51]. The mean of the measurement points in a specific spatial distance (with a radius of 30 m) was calculated to eliminate possible gross errors. Figure 10a shows the time-series deformation of selected points (P1 and P2) at zone (I). Due to different observation geometries of the descending and ascending data for the direction of the slope movement, the LoS time-series deformations at point P1 show opposite trends. That is consistent with the northwest movement detected from the optical POT result. No obvious acceleration in surface movement was detected from the descending and ascending measurements before the multi-hazard chain on 17 June 2020. Before the occurrence of the multi-hazard chain, only a 6 mm change in cumulative deformation was observed at point P1 between 1 June 2020 and 13 June 2020 from the ascending track. A sudden acceleration in surface movement can be observed at point P1 from the ascending track between 25 June 2020 and 7 July 2020, moving at approximately 14.5 mm along the LoS direction within the final acquisition (12 days). According to the local rainfall data collected from 22 June 2020 to 19 August 2020, a maximum daily rainfall of 10.2 mm was recorded on 5 July 2020 (Figure 2b). This sudden change in surface movement was responding to the intense daily rainfall recorded on 5 July 2020. An accelerated displacement was also observed from a field survey on 5 July [41], showing that our InSAR result is in good agreement with the field data. Conversely, a linear growth trend can be observed from the descending measurement at point P1 until 2 July 2020. A gentle acceleration can be observed at point P2 from the descending measurement between 8 June 2020 and 2 July 2020, and the cumulative deformation over the observation period is up to 110 mm.

At zone (II), the cumulative deformations at point P3 and P4 are 83 mm and −130 mm for the observation period, respectively, which are significantly larger than those at point P5 from both datasets, as shown in Figure 10b, and are probably caused by serious washout and erosion of the river and the outburst flood caused by the Meilong debris flow at the upstream and middle section of zone (II), where point P3 and P4 were located (Figure 8). The rainfall season started in June with the intense rainfall of 38.1 mm recorded before the multi-hazard chain on 17 June 2020. During this period, a clear acceleration in deformation can also be observed at point P4 from the descending track, falling about 17 mm between 8 June 2020 and 2 July 2020. A sharp drop in cumulative deformation was obtained from the ascending measurement at point P3 between 13 June 2020 and 25 June 2020 when the multi-hazard chain occurred, which could have been caused by the collapse of the landslide toe at the upstream section after the multi-hazard chain.

Figure 10c shows the time-series deformation at zone (III), the upstream section of the ANZ landslide. The cumulative deformations at point P6 are 46 mm and 71 mm for the descending and ascending measurements, respectively, which are both smaller than those at zones (I-II). Similarly, there is a sudden acceleration in displacement observed at point P6 from the ascending datasets from 1 June 2020, rising about 30 mm in about one month. This is consistent with the rainfall season in June and the intense daily rainfall recorded on 5 July 2020 (Figure 2). However, no obvious acceleration was detected from the descending datasets.

The landslide shoulder zone (IV) is located at the upper section of the ancient landslide area and close to the downstream section. Before June 2019, the cumulative displacements for all of the two selected points in both tracks were fluctuating within a range of −20 mm to 20 mm, which indicated that this section was in a relatively stable stage.

At point P7 (Figure 10d), both the descending and the ascending measurements captured a clear acceleration at the last three acquisitions, mainly after the multi-hazard chain and with the increasing monthly rainfall in June and July, moving 22 mm and 30 mm within one month, respectively. Point P8 is located at the pre-failure flow in the upper section detected from the optical image analysis (Figures 4 and 5). The cumulative deformation is smaller at point P8 in both datasets, mostly ranging from −30 to 30 mm. It is noted that a similar acceleration of deformation was also captured by the ascending measurements at point P8, rising about 25 mm in the last two acquisitions, which is consistent with the ascending measurement at point P7.

**Figure 10.** Pre-failure time-series plotting for selected measurement points from the different zones: (**a**) central zone (I); (**b**) landslide toe zone (II); (**c**) upstream zone (III); (**d**) landslide shoulder zone (IV); (**e**) unstable slope (V) zone in the pre-failure mean velocity map (Figure 8). Red dashed line indicates the reactivation of the ANZ landslide.

The time-series deformation over the detected unstable slope, located on the right side of the ANZ landslide, is shown in Figure 10e. The maximum cumulative deformations are 63 mm and 126 mm for the descending and ascending measurements, respectively. It was found that acceleration in deformation was detected from the descending measurement since the end of May in 2020, moving 22 mm within the final four acquisitions, which agreed with the increasing monthly rainfall in June. A dramatic change of cumulative deformation was also found at the final two acquisitions (25 June 2020 to 7 July 2020) from the ascending measurement, surging 30 mm within just 12 days, which is a response to the intense rainfall recorded on 5 July 2020, suggesting that this unstable slope may develop to failure after this event and that continuous monitoring should be further conducted over this site.

To further explore the deformation evolution in the spatial domain, the cumulative LoS deformation maps of the descending and ascending measurements are depicted in Figures 11 and 12, respectively. Sixteen acquisitions were selected and referred to the initial acquisition dates on 27 March 2018 and 8 March 2018 for the descending and ascending maps, respectively. From the joint analysis of the descending and ascending cumulative deformation maps, the central and toe sections of the ANZ slope were found to move initially in January 2019. It is noted from the descending map that the deformation detected at the central section gradually expanded to the surrounding areas with increasing cumulative deformation within the reactivated ANZ landslide area. Most of the points in the upper section were in a relatively stable state with much smaller cumulative deformations. Most importantly, no obvious accelerations were found in the cumulative deformation before the multi-hazard chain on 17 June 2020 from either of the measurements. However, clear changes in deformation can be observed over several sections of the ANZ slope from both measurements after the multi-hazard chain.

#### *5.2. Deformation Mechanism and Triggering/Preparatory Factor Analysis*

When comparing the Spatial-Temporal evolution between the different zones above, it is clearly seen that the movement of zones (I-II) at the central and toe sections of the ANZ landslide has occurred since September 2018, long before the obvious movement detected at zone (IV) in the upper section since June 2019. After the expansion of deformation at the central section and following the collapse of the toe section, the upper section gradually lost support at the bottom and presented movement. It is clearly seen that the deformation at the upper section was much smaller than that from the central and toe section, which indicated that the ANZ landslide was characteristic of the retrogressive failure mechanism [69]. After the occurrence of the multi-hazard chain on 17 June 2020, obvious accelerations of deformation from several sections of the ANZ landslide were clearly captured by the Spatial-Temporal analysis of the time-series InSAR result. It is evident that the reactivated ANZ landslide presented an acceleration between 17 June 2020 and the beginning of July 2020, which is highly consistent with the field survey [40,41]. At the same time, an unstable slope next to the ancient ANZ landslide was detected from the time-series InSAR analysis, and it also presented an accelerated deformation similar to that of the reactivated ANZ landslide after the multi-hazard chain. Therefore, continuous monitoring over this area is essential for further assessment of the post-failure stability and the detection of secondary landslides.

**Figure 11.** Descending cumulative LOS deformation of the measurement points located at the Aniangzhai slope for the 16 selected acquisition dates of Sentinel-1 descending datasets, referring to the first SAR image acquired on 27 March 2018.

It is noted that a recently published work [70] found a clear acceleration over the upper parts of the ANZ landslide starting in spring 2020 from the multi-temporal InSAR analysis with the descending Sentinel-1A dataset [70]. However, the InSAR results obtained in this study suggested that no obvious accelerations were detected before 17 June 2020, with both the descending and the ascending Sentinel-1A datasets, which is consistent with another published InSAR result with shorter temporal coverage (study period from 20 April 2018 to 20 May 2020) [41]. Such inconsistencies between the InSAR results could be caused by different datasets and the pre-processing and post-processing strategies. Most importantly, it is consistent that all three InSAR results confirmed the clear displacements over the ANZ slope before the event, ranging from −50 to −80 millimetres/year.

**Figure 12.** Ascending cumulative LOS deformation of the measurement points located at the Aniangzhai slope for the 16 selected acquisition dates of Sentinel-1 descending datasets, referring to the first SAR image acquired on 8 March 2018.

According to the rainfall data collected from the Anianggouer weather station, the cumulative rainfall from 1 May to 17 June 2020 was over 250 mm and the daily rainfall on 17 June 2020 was up to 38.1 mm, as shown in Figure 2a [40]. The Spatial-Temporal evolution of the ANZ landslide from the time-series InSAR analysis clearly revealed that there was no obvious acceleration (precursory signal) detected before the failure of the multi-hazard chain on 17 June 2020, which is consistent with the published results [41]. This also suggested that the intense rainfall from May to June 2020 was one of the most significant triggering factors for the ANZ landslide. In addition, several preparatory factors also contributed to the occurrence of the ANZ landslide. Firstly, the short, intense rainfall event on 17 June 2020 indeed triggered the failure of the Meilong debris flow, and then, the undercutting river and outburst flood from the dammed lake seriously damaged the stability of the toe section of the ANZ slope. The following collapse of the toe section further reactivated the ANZ landslide and formed the accelerations of deformation. Secondly, the ANZ landslide was located on a region with numerous active faults where numbers of strong earthquakes occurred. The frequent occurrence of earthquakes would significantly shake the ANZ slope and develop fractures and cracks over the slope, which could make the slope looser and slide more easily. At the same time, the dip angle of the landslide toe is up to 70◦. It is prone to fail under the effect of gravity and long-term washout and erosion of the underlying river. The ANZ landslide body mainly consists of decimetre-scale angular and sub-angular blocks, which make it easier for the rainwater and groundwater

to infiltrate into the surface layers of the slope. This would increase the weight of the soil in the surface layer and result in the reduction in the shear strength of the soil and gradually contribute to the generation of a sliding surface over the slope. Besides that, based on the high-resolution true colour images shown in Figure 4, multiple man-built roads laterally crossed the middle section of the ANZ slope. When intense rainfall occurs in a short duration, a large amount of rainwater will infiltrate the surface layer of the slope. These roads in the middle section could temporally block the flow of rainwater over the surface, which would probably have an influence on the movement of the superficial layers over the slope. It is also consistent with the larger magnitude of deformation in the middle section detected from the optical POT result in Figure 7. At the same time, from the true colour images, there are no roads built in the upper section. Therefore, the magnitude of deformation is much smaller from the optical POT result. The time-series InSAR result also confirmed that the cumulative deformation in the upper section is smaller than the lower section. Thus, it is believed that human activities might also have some influence on the more superficial movement of the slope. Above all, the ANZ landslide was likely triggered by the heavy rainfall, and several preparatory factors, including wash-out and erosion of the river and outburst flood, seismic activities, and local terrain conditions over the slope, also make significant contributions to this event.

#### **6. Conclusions**

In this study, the deformation characteristics and Spatial-Temporal evolution of the ANZ landslide were investigated in detail by the joint use of multi-temporal optical images and time-series InSAR analysis. The analysis of multi-temporal optical images clearly showed the spatial evolution of pre- and post-failure landslide features over the slope. Several pre-failure flows over the slope were observed. The post-failure displacement was retrieved with the optical POT technique. These optical POT results revealed that the maximum horizontal displacement was up to 14.4 m, and the ANZ slope mainly moved towards the northwest direction. The deformation velocity and time series of the reactivated ANZ landslide was also retrieved using time-series InSAR analysis with Sentinel-1A/B datasets from the descending and ascending tracks. The maximum magnitude of the LoS deformation velocity was up to −70 mm/year and 45 mm/year from the descending and ascending measurements, respectively. The large deforming areas were found at the landslide toe collapse zone. A large deforming area with a velocity of −20 to −50 mm/year was clearly observed from the descending measurement at the central zone over the ANZ slope. At the same time, an unstable slope was detected on the right side of the ANZ landslide. Most importantly, the time-series analysis of multiple sections over the reactivated ANZ landslide showed that no obvious accelerations of deformation (precursory deformation) were detected before the multi-hazard chain on 17 June 2020. However, after the occurrence of the multi-hazard chain, obvious accelerations can be observed over the ANZ slope, which suggested that the ANZ landslide entered the acceleration status after the reactivation. Through comparing the Spatial-Temporal evolution between different zones of the slope, it is evident that the ANZ landslide presented retrogressive failure mode. The Spatial-Temporal evolution and deformation behaviour derived from the time-series InSAR result suggested the occurrence of the ANZ landslide is not the result of a single triggering factor due to heavy rainfall, but a joint effect of several preparatory factors, including frequent seismic activities, serious wash-out, and erosion of the river and the outburst flood and local terrain conditions.

Satellite SAR interferometry demonstrated its effectiveness and usefulness for mapping and monitoring the surface deformation caused by landslides over a long period. At the same time, the optical remote-sensing technique can provide valuable supplementary information about landslide features and extremely large deformation-of-slope movements. On-going monitoring of the reactivated ANZ landslide is recommended as the ANZ slope presented a clear acceleration after the multi-hazard chain. Therefore, continuous monitoring of the slope is essential for the detection of further failure and for emergency

response. Thus, further investigation of post-failure monitoring with the time-series InSAR technique and optical time-series analysis with the POT technique will be studied in our future research.

**Author Contributions:** Conceptualization, J.K. and A.H.-M.N.; methodology, J.K. and A.H.-M.N.; validation, J.K. and A.H.-M.N.; formal analysis, J.K. and A.H.-M.N.; investigation, J.K.; resources, J.K., A.H.-M.N. and L.G.; data curation, J.K.; writing—original draft preparation, J.K.; writing—review and editing, J.K., A.H.-M.N. and L.G.; supervision, A.H.-M.N. and L.G.; project administration, A.H.-M.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Program for Guangdong Introducing Innovative and Entrepreneurial Teams (2019ZT08L213), Natural Science Foundation of Guangdong Province (grant number 2018A030310538 and 2021A1515011483).

**Acknowledgments:** The authors thank ESA for providing the Sentinel-1A/B data and the Planet Labs Company for providing the PlanetScope data under a research and education license.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

