The June 2020 Aniangzhai landslide in Sichuan Province, Southwest China: slope instability analysis from radar and optical satellite remote sensing data

A large, deep-seated ancient landslide was partially reactivated on 17 June 2020 close to the Aniangzhai village of Danba County in Sichuan Province of Southwest China. It was initiated by undercutting of the toe of this landslide resulting from increased discharge of the Xiaojinchuan River caused by the failure of a landslide dam, which had been created by the debris flow originating from the Meilong valley. As a result, 12 townships in the downstream area were endangered leading to the evacuation of more than 20000 people. This study investigated the Aniangzhai landslide area by optical and radar satellite remote sensing techniques. A horizontal displacement map produced using cross-correlation of high-resolution optical images from Planet shows a maximum horizontal motion of approximately 15 meters for the slope failure between the two acquisitions. The undercutting effects on the toe of the landslide are clearly revealed by exploiting optical data and field surveys, indicating the direct influence of the overflow from the landslide dam and water release from a nearby hydropower station on the toe erosion. Pre-disaster instability analysis using a stack of SAR data from Sentinel-1 between 2014 and 2020 suggests that the Aniangzhai landslide has long been active before the failure, with the largest annual LOS deformation rate more than 50 mm/yr. The 3-year wet period that followed a relative drought year in 2016 resulted in a 14% higher average velocity in 2018–2020, in comparison to the rate in 2014–2017. A detailed analysis of slope surface kinematics in different parts of the landslide indicates that temporal changes in precipitation are mainly correlated with kinematics of motion at the head part of the failure body, where an accelerated creep is observed since spring 2020 before the large failure. Overall, this study provides an example of how full exploitation of optical and radar satellite remote sensing data can be used for a comprehensive analysis of destabilization and reactivation of an ancient landslide in response to a complex cascading event chain in the transition zone between the Qinghai-Tibetan Plateau and the Sichuan Basin.


Introduction
Landslides are widespread geological hazards in mountainous regions worldwide. Once a landslide mass loses its stability, it could induce strong destructiveness. Landslide processes are complex and often comprise different process types. Some of them move fast (Quecedo et al. 2004), but several other landslides also take place slowly and steadily at the beginning, and then accelerate suddenly terminating in catastrophic avalanche-type or collapselike movement styles (De Blasio 2011). Landslides have occurred more frequently due to increased urbanization, continued deforestation, and increased extreme weather events (Schuster 1996;Biasutti et al. 2016;Lee 2017). To monitor landslide disasters and build effective early warning systems (EWSs), the adopted technical means should meet at least the following requirements: Adequate regional coverage and temporal sampling capacity, sufficient measuring accuracy related to the velocity of the monitored processes, and good cost performance. Ground-based methods, such as continuous GNSS for landslide monitoring, are difficult to set up and implement in mountainous and remote areas (Akbarimehr et al. 2013). Instead, optical and radar satellite remote sensing plays a promising role in driving innovation in large-scale detection, monitoring, and assessment of landslide hazards and can be quite useful to incorporate in the framework of multidisciplinary disaster risk reduction (DRR).
Cross-correlation of optical images can be used to assess the kinematics of large slope failures (Travelletti et al. 2012;Yang et al. 2020). Furthermore, automated and semi-automated approaches using time series of multi-sensor optical images have already been developed to create multi-temporal inventories by identifying landslide areas based on changes in vegetation cover (Behling et al. 2016;Yang et al. 2019). Optical remote sensing data have become more commonly available in recent years and are easily understood and handled by non-experts . However, clear sky images may not be readily available prior to and during a given landslide event. Moreover, the displacement accuracies retrieved from cross-correlation analyses are highly dependent on the resolution of the optical remote sensing acquisitions and the satellite's precise orbit position and orientation posture (Debella-Gilo and Kääb 2011). Hence, optical remote sensing has limited use in reliably supporting near real-time hazard assessments and EWSs.
Synthetic aperture radar (SAR) offers new opportunities to support the systematic detection and monitoring of landslides over extensive regions and for the development of regional-scale landslide warning systems (Colesanti and Wasowski 2006;Bianchini et al. 2013;Herrera et al. 2013;Motagh et al. 2013). With synoptic imaging capabilities, under inclement weather conditions and independent of sunlight conditions, SAR techniques provide invaluable information on landslide locations, boundaries and changes to vegetation within landslide bodies, based on the exploitation of radar amplitude and phase information. Advanced multi-temporal InSAR (MTI) methods, e.g., persistent scatterer interferometry (PSI) and small baseline subsets (SBAS), can be exploited to evaluate subtle changes in landslide creep rates in response to external triggering factors; these changes can indicate impending failures (Teshebaeva et al. 2015;Handwerger et al. 2019;Hu et al. 2020). The new surge in available SAR data via Sentinel-1 (S1) satellites has provided golden opportunities to use SAR sensors as operational instruments for landslide hazard assessments (Solari et al. 2019) and temporal predictions of large failures (Mantovani et al. 2019). In particular, S1 data have higher spatial resolution and global dual-polarization coverage with improved revisit times of 6-12 days over the data from previous C-band SAR missions such as ERS and Envisat. As S1 data are available at no cost, there has also been growing interest from scholars for objective change detections, landslide hazard assessments, and potential techniques for multidisciplinary DRR (Feng et al. 2015; Barra et al. 2016;Dai et al. 2016;Intrieri et al. 2018;Dini et al. 2020;Dai et al. 2020).
On 17 June 2020, close to Aniangzhai village of Danba County in Sichuan Province of Southwest China, a massive landslide of ∼ 6 million m 3 Yan et al. (2021) was partially reactivated. The main triggering factors were the undercutting effects and erosion on the toe of the landslide body from the overflow of a dammed lake (height of nearly 8 ∼ 10 meters), which was created by debris flows coming from the northern Meilong valley comprising a complex cascading event chain. Firstly, the heavy rainfall in summer 2020 induced debris flows in the Meilong valley. With the help of Sentinel-2 (S2) optical images, we observe that the debris flow generated from the valley regions north of the reservoir flowed toward the south. Then, the washed-out stones and soils formed a barrier dam just under the ancient Aniangzhai landslide body and blocked the Xiaojinchuan River, leading to an increase in the water level (seeing supporting material: Fig. S1). Thereafter, the overflow of the barrier dam, influenced by the discharge of the surplus water from the nearby hydropower station to reduce the flood pressure, undercut the toe of the landslide, resulting in partial reactivation of this ancient landslide body. Soon after the lower part of the landslide area collapsed gradually. In this case, this specific cascading event chain of "rainfall-debris flows-dammed lake-outburst floodserosion-landslide" was formed and threatened a dozen villages downstream, resulting in an evacuation of more than 20000 people to abandon and leave their home towns (Yan et al. 2021).
In this study, we report investigations on ground deformation of the Aniangzhai landslide before and during June 2020 failure using optical and radar satellite remote sensing data. Sub-pixel crosscorrelation of high-resolution optical images from Planet is utilized to obtain information on the main landslide failure, e.g., horizontal movement and moving direction. Then, S1 SAR data are analyzed using the MTI techniques to assess the slope instability between 2014 until the time of failure. The results are then analyzed against changes in meteorological conditions to assess the long-term and transient behavior of the Aniangzhai landslide. We also evaluate a method for anticipating the time of failure based on MTI results using a modified inverse-velocity method. Finally, we introduce some findings based on the abnormal behaviors of the Normalized Difference Vegetation Index (NDVI) and interferometric coherence over the landslide mass before the failure.

Environmental and geomorphological settings
Danba County is located on the southeastern edge of the Qinghai-Tibet Plateau, and Aniangzhai village is located in the center of Danba County. The geomorphological structure of the region comprises high mountains and narrow valleys with an average elevation of approximately 1800 meters. The June 2020 Aniangzhai landslide is a case of a partial reactivation of an ancient and larger slope failure (Zhao et al. 2021). The original ground surface of the Aniangzhai landslide has an elevation of approximately 2000∼ 2500 meters. The topographic profile and possible thickness of the landslide body were investigated in a recent study exploiting UAVs (Zhao et al. 2021), suggesting a maximum thickness of approximately 60 meters. Based on our field investigation and (Zhao et al. 2021), the vertical component of motion was significant at the head part of the failure, while the landslide slipped down as a whole in the middle-lower part. Therefore, we assume that the Aniangzhai landslide has a rotational-translational mechanism, with rotational component being more significant in the upper part and translational component becoming predominant in the middle-lower part. Moreover, this region is located in the upper reaches of the Yangtze River, which is full of water resources. The foot of the Aniangzhai landslide area reaches the Xiaojinchuan River. There is also a dam nearby, which is very close to the failure region upstream (Fig. 1). We compared the river courses in June over the last 3 years before the 2020 failure (seeing supporting materials: Fig. S2). The river courses demonstrate similar extents and appearances in 2018 and 2019, regulated by the reservoir upstream. In contrast, the river course during the 2020 failure event shows major inundation due to surplus water from the reservoir. The annual mean temperature for this region is approximately 14 • C. However, due to elevation changes, the differences between the top of the mountain and the valley could be greater than 24 • C. Because of the plateau monsoon climate and complicated geomorphology, lots of natural disasters are taken place frequently in this region, especially landslide hazards. Figure 2 illustrates a Skysat optical image and several photos from the fieldwork, in which different zones of the landslide area are highlighted. The high-resolution optical image from Skysat was acquired on 25 November 2020 with an accuracy of half-meter (Fig. 2f). The landslide area lays on the hillside north of the country town, and it was previously described as an ancient rockfall area. The red line in Fig. 2f, indicates the center part of the failure, which had the largest deformation in this event; the orange line indicates the medium motion of approximately 1 ∼ 5 meters, while the yellow line reveals the extent of the whole landslide body. It is obvious to see the main scarp on the head of the failure area, as well as the erosion on the toe of the landslide. As demonstrated in Fig. 2e, the vegetation on the toe was washed away and a steep valley was formed in the front of the landslide body due to enormous mass loss. And this directly triggered the reactivation of the ancient landslide body. The deformation then started in the upper part, and the whole block was moving downwards. Figure 2b, h and i shows the boundaries of the landslide failure. Figure 2b reveals the northern lateral flank of the landslide, where cracks, approximately 1 ∼1.5 meters wide, caused by block motion are clearly visible. Figure 2h displays the southern flank of the central failure, which had the fastest displacement rates during the 2020 failure. Figure 2i shows the main scarp of the failure in the southeast direction. Other pictures reveal examples of damages that occurred in this disaster. Figure 2a and d show the ravaged roads, which were broken and had a drop of a few meters. Figure 2c and g displays a damaged EHV (Extra High Voltage) transmission tower and a cracked house in this event.

Remote sensing optical images
We use two high-resolution optical satellite images acquired by Planet Lab (Team 2017) satellite constellation to assess the horizontal kinematics of the failure, i.e., extent, direction and magnitude. The satellite data is acquired right before and after the event, i.e., on 15 June 2020 and 24 June 2020, respectively (Fig. 1). The used Planet Lab satellite imagery has a resolution of about 3 meters. Indeed, the horizontal displacements are quite obvious when these two images are superimposed on each other. The Planet Lab images comprise three multi-spectral bands covering the visible part of the spectrum (RGB). The red band is used in this study with the best root-meansquare error in image registration among these three bands. The two Planet images are cropped to the same subset covering the landslide area forming the input to the cross-correlation analysis using the COSI-Corr software package to estimate the horizontal displacements (Leprince et al. 2007). The cross-correlation is estimated with steps of 2 × 2 , which provides the west-east (W-E) and north-south (N-S) horizontal displacements calculated by every two steps. Then a median filter is applied. In the end, the magnitude of displacement in each pixel is calculated as the norm of vectors from the results in two directions.

MTI analysis using Sentinel-1 SAR data
We apply the C-Band SAR images acquired by S1 satellite for MTI analysis in this study (Copernicus 2020), specific information of the acquisition can be found in supporting materials (Table S1). In detail, the InSAR processing is carried out immediately after the failure with 89 descending images of S1 Interferometric Wide (IW) swath mode from October 2014 to July 2020. Among them, a few images in 2014-2015 cover the study area partially and they are stitching together for the exploitation. The spatial resolution is approximately 5m × 20m with a 250 km swath. There are both ascending and descending datasets available. As seen from optical images, the main direction of the slope is toward the north-west, which causes foreshortening effect in ascending data. Thus, the descending data are selected for our analysis in this case. Due to the temporal gap of the original dataset in 2017-2018, the MTI processing is carried out in two temporal frames, i.e., 2014-2017 and 2018-2020. During the processing, the time series of SAR images are wellcoregistered and then cropped to the identified subset of the landslide area. The subset covers an area of approximately 26 square kilometers ( 4.7 × 5.6 km). The 2000 SRTM DEM (30m) is utilized for geocoding and estimating the topographic phase component in InSAR processing (Farr et al. 2007). The processing chain of S1 has already been mentioned in many previous studies (Fattahi et al. 2016;Grandin et al. 2016;Yagüe-Martínez et al. 2016;Haghshenas Haghighi and Motagh 2017). The traditional InSAR has limitations for landslide monitoring. The main limitations are the widespread loss of coherence between consecutive image acquisitions and atmospheric disturbances (Zebker and Villasenor 1992;Wasowski and Bovenga 2014). Thus, we apply the advanced MTI techniques to mitigate the problem and retrieve the information of displacement, i.e., PSI and SBAS methods. The GAMMA and StaMPS software packages are used for the implementation of interferometric and MTI analysis (Hooper et al. 2007(Hooper et al. , 2012Wegnüller et al. 2016), with atmospheric correction obtained using the Generic Atmospheric Correction Online Service (GACOS) product Morishita et al. 2020). In PSI processing, a stack of single-master interferograms is generated and the pixels with the highest signal-to-ratio values are selected (Hooper et al. 2004(Hooper et al. , 2007. Such pixels are regarded as Persistent Scatterers (PS), mostly come from rocks and man-made features. As for the SBAS method, the algorithm exploits a network of small temporal and spatial baselines to minimize the decorrelation between image pairs (Lanari et al. 2007;Anderssohn et al. 2009). The distributed scatterer (DS), which is defined as the pixel that shares similar statistical behavior with its neighbouring pixels, is taken into account. MTI baseline networks and selecting criteria can be found in supporting materials ( Fig. S3-S6). With the help of MTI techniques, we could obtain comparable results between PSI and SBAS methods for analyzing slope instability. In addition, MTI time series are further exploited using inverse-velocity (INV) theory to predict the time of failure.

Auxiliary data
To better understand the dynamics of the Aniangzhai landslide in relation to potential influencing factors, some auxiliaries are exploited (Table 1). The first auxiliary includes precipitation retrieved from the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS). Spanning all longitudes, CHIRPS incorporates 0.05 • resolution ( ∼5km) satellite imagery between 50 • S-50 • N and in-situ station data to create gridded rainfall time series. The precision of the rainfall datasets is sufficient for applications and exploitation at the regional scale. In our study area, available CHIRPS data cover a time span of 20 years between 2000 and 2020, and the precipitation is calculated for Danba County.
The second auxiliary includes multiple optical remote sensing collections to obtain NDVI values. The NDVI value is calculated as follows: where RED is the red portion of the electromagnetic spectrum and NIR is near-infrared light. In this study, NDVI time series from three different satellite datasets are exploited and compared. Details can be found in Table 1. The MODerate-resolution Imaging Spectroradiometer (MODIS) Reflectance product MCD43A4 provides daily reflectance data adjusted using a bidirectional reflectance distribution function (BRDF). Data of both Terra and Aqua satellites are used in the generation of this product, providing the highest probability for quality assurance input data (DAAC 2021). For comparison and validation, two other collections of optical remote sensing satellites are applied, i.e., the Landsat-8 collection (16-day temporal resolution and 30-meter spatial resolution), and S2 data (5-day temporal resolution and 10-meter spatial resolution). To be noticed is that the S2 dataset for this study area is only available since late 2018. The precipitation and NDVI analyses were conducted with the help of the Google Earth Engine (GEE). We developed our own scripts to generate the monthly-mean and yearly precipitation during 2000-2020 for Danba County for further exploration in this study, while NDVI is calculated or obtained for the slope affected by the landslide from the mentioned three satellite collections during different periods (Table 1). The NDVI values are further compared with the interferometric coherence. The purpose of the comparison is to see whether some features could be obtained for early warning without complex MTI processing (Jacquemart and Tiampo 2021).

Inverse-velocity theory for anticipating the time of failure
When landslides, rockfalls and similar hazards are investigated, one of the major interests is to predict a potential time range when a failure might be likely to happen. For this goal, already several approaches have been developed, among them, the INV method which is considered to be a simple but effective method for EWS being used in many studies during recent years (Carlà et al. 2017;Zhou et al. 2020).
In order to apply INV, the first step is to calculate the velocity of LOS displacement from the time series of displacement. The calculation of landslide velocity is always a complicated problem. On the one hand, the strength parameters for different landslide types should be considered in the calculation. On the other hand, the friction coefficient and friction resistance will change with different stages of activities and the volumes of the landslide (De Blasio 2011). The key challenge is that in reality the observations of displacement could be influenced by manmade or systematic noises. Such noises include measurement errors, random instrument noises and noises from periodic changing factors such as rainfall, groundwater, human activities, etc. These noises could lead to outliers and abnormal behaviors for INV, which makes data smoothing necessary. In previous studies, some approaches are exploited to generate the smoothing of the displacement, such as short-term and long-term moving averages and exponential smoothing functions (De Blasio 2011;Carlà et al. 2017). In this study, we have applied these different approaches, but the outcomes have not been satisfying. The reason for this is, if the kernel of the smoothing function is too large, the filtered curve would possibly lose some important features, whilst the noise in displacement could not be improved using a smaller kernel.
In order to obtain ideal fittings which capture the relevant features and to mitigate the influence from noises, we propose a method to smooth the displacement values obtained by MTI processing. The method uses least squares and L1 regression under the assumption that after the main failure has happened, the further displacement occurring within the landslide area remains more or less constant. In this context, we introduce the parameter c to represent the limitations of the measurements, whereas c 1 and c 2 represent the minimum and maximum detectable displacements, respectively. In specific, if the magnitude of displacement is larger than c 1 , then the slope movement is considered to occur in form of sliding. Since the obtained MTI measurements are characterized by cm to mm precision (Osmanoğlu et al. 2011;Wang et al. 2012;Motagh et al. 2017;Haghighi and Motagh 2019), we introduce a relatively generous threshold amounting to c 1 is 0.01 meter. Since we do not want to over-smooth the features caused by the landslide failure in the fitting process, the parameter c 2 is not set in this study. In the result, we calculate the smoothed displacement by the following equation: where y represents the observation, x is variable and A comprises the sparse matrix for the tridiagonal representation of the standard second difference operator, and is the factor balancing the fitting and sparsity. To solve (2), a package for solving convex optimization problems Boyd 2008, 2015) was used, to derive x that minimizes expression (2) being subject to the constraints while using identical parameters. Values of INV will approach zero corresponding to the increasing time as velocities increase asymptotically closer to the failure. Once the smoothed displacements are generated, INV could be derived and thus, a prediction of the failure could be achieved. Figure 3 illustrates the horizontal displacement calculated from Planet optical images for a short period of time comprising the situation right before and after the failure, i.e. 1 day before the failure and 8 days after the failure. The applied two Planet images are shown as background in Fig. 3. The main component of the horizontal displacement occurs in the W-E direction with a maximum displacement of approximately 13.2 meters toward the west. In the N-S direction, the displacement is oriented toward the north with a maximum displacement of approximately 6.9 meters. Overall, the absolute horizontal displacement is estimated as the norm of vectors from displacements in N-S and W-E directions, and the maximum magnitude reaches approximately 14.7 meters in the NW direction. The result also demonstrates that some localized deformations exist out of the main body of the failure, mainly in the northwest and southwest corners of the area shown on the image in Fig. 3. The moving directions of those localized deformations are different compared to the ones obtained for the main failure. This subset area is shown with significant motion comparing to the surrounding areas. From Fig. 3, we can see that the maximum horizontal displacement rate could reach ∼1.6m/day, which is too large to be applied using InSAR monitoring . Figure 4 shows a comparison of the results of MTI processing for the two different periods using both PSI and SBAS methods. The displacement rates have been derived in line-of-sight (LOS) direction, whereas positive values represent motion towards the satellite, whilst the negative values represent motion away from the satellite.

MTI analysis
The reference point, representing a stable area during the whole time period is selected outside of the landslide region in the northern hill slope. From the MTI analysis, it is deducted that the area of the June 2020 failure had already experienced movements prior to the actual failure, especially in the center part of the landslide body. As seen in Fig. 4, the creeping movement could already be revealed within the landslide body for the time period of 2014-2017. In this period, the maximum displacement rate in LOS direction amounts to approximately −24 and −38 mm/yr, for the PSI and SBAS methods, respectively. For the 2018-2020 period up to the failure, the displacement rates within the area of the June 2020 failure have significantly increased compared to the 2014-2017 period, reaching the maximum of approximately −40 and −55 mm/yr for the PSI and SBAS results, respectively, in the center of the landslide body. Moreover, the areas outside of the landslide failure turn out to be stable in general. The relevant statistics for the PSI and SBAS results for the ancient slope failure reactivated in 2020 can be found in Table 2. Figure 5 shows a comparison of displacements time series along a selected topographic profile (marked by blue in Fig. 4d) for the SBAS results in 2018-2020. Point T1 is situated on the channel floor and points T2-T4 are located on the partial failure part; while T5 and T6 are situated on the upper part of the ancient landslide and outside of the failure body. The displacements rates of selected points in the top (T5 and T6) and bottom (T1 and T2) zones are smaller compared to the displacements rates in the central failure zones (T3 and T4). Meanwhile, the central failure zones have a relatively larger slope inclination compared to the top and bottom zones. The above results are further elaborated in the Discussion. Figure 6a demonstrates the annual precipitation amounting to 855.5 mm in 2014, then increasing by 5.7% in 2015, but decreasing by 9.0% in 2016. In contrast, for the period of 2017-2019, we observe a constant increase by 12.3%, 1.1% and 5.2% respectively compared to the previous year. Figure 6b shows that from April to June 2020, rainfall is 30.5%, 4.0% and 26.4% higher than the long-term average for the last 20 years, respectively.

Influence of precipitation on the kinematics of the landslide
To better quantify the role of the 2020 heavy precipitation in influencing the kinematics of the landslide, we focused on the 2018-2020 period and analyzed the time series of LOS displacements at different parts of the landslide. Figure 7 illustrates the locations of ten points selected arbitrarily over the whole landslide body: points P1-P3 from the head of the failure part (Zone I); points P4-P6 from the central failure body (Zone II); points P7-P9 from the foot of the landslide (Zone III); and point P10 from the landslide body, but outside of the 2020 failure.
The results of MTI time series for the selected points are displayed in Fig. 8; the dot lines show the time series retrieved from the SBAS processing, while the blue lines represent the fitting curves obtained based on the smoothing methodology described in Sect. 3.4. As seen in Fig. 8, points P1-P3, which are located in Zone I, show accelerations and decelerations throughout the entire time series with a clear accelerating trend as of spring season in 2020. In contrast to points P1-P3, the displacement rates at points P7-P9 are relatively lower; they exhibit periods of acceleration in the step-wise pattern and a constant velocity at the end of the time series. Points P4-P6 in the center region show the highest displacement rates among all the selected points with fewer variations in overall velocity compared to other parts of the landslide, although some periods of slowing-downs and accelerations occur in the same period as in case of points P1-P3 (e.g., in June 2019). Point P10, located outside of the 2020 failure, reveals a steady movement at the beginning, punctuated by a short episode of acceleration in March-April 2019 and another obvious acceleration since spring 2020. Figure 9 reveals the comparison of the LOS velocity of points P1-P3 on the image of Fig. 9a and the corresponding precipitation in the period of 2018-2020. These three points are selected for INV processing due to the strong correlation of the simultaneous speeding up of their displacement rates in response to the increasing cumulative precipitation in the rainfall season. Obviously, there are three rainfall seasons in Fig. 9b, corresponding to the rainy months of May to September. We calculate the increments of the precipitation in different years and then compare this with the amount of changes in the corresponding LOS displacements of the selected points in the same duration. The results are displayed in Table 3. Here we observe that for the first time period of mid-April to mid-June in 2018, where rainfall reaches 248 mm, approximately 6.5, 6.5 and 2.3 mm increases in displacement values are observed at points P1-P3, respectively. Interestingly, compared to the period t 1 , the increments of displacement in t 2 of these points decline corresponding to the reduced rainfall, i.e., when the rainfall in t 2 drops by approximately 12%, the variations of displacement of points P1-P3 also drop by approximately 36%, 53% and 79%, respectively. However, for the third time span t 3 , we can see that since the rainfall raises by 22% compared to t 2 , the displacement increments of points P1-P3 also increase by 59%, 160% and 456% than the values in t 2 .
INV results for anticipating the time of failure Figure 10 demonstrates the results of INV analyses. In Section 3.4, we described anticipating the time of failure using the modified INV theory. The period considered in the INV analyses is from April 2020 to mid-June 2020, which is affected by the heavy precipitation before the failure, and is the same for all the selected points. It is worth noting that, the areas with accelerating LOS displacement before failure are needed to be evaluated for INV analysis (Manconi and Giordan 2016). Otherwise, the method leads to underestimation or overestimation of the failure time. This has been shown in Fig. 10b for points P5 and P7, which are located in an area, where no acceleration was found before the failure. Applying the INV method to these two points resulted in overestimated prediction times, i.e., approximately 20 and 66 days after failure. The results presented in Section 4.3 shows that among the whole landslide body, only the top of the failure area indicated accelerated creep since spring 2020 in response to the heavy rainfall (Zone I in Fig. 7). By performing INV analysis for the points in this area, we observe that the INV can predict the time of failure properly (Fig. 10a). Figure 11a shows the negative correlation between coherence and NDVI (MODIS) in 2014-2020. Regardless of the temporal gap in SAR data availability, the changes in coherence and NDVI show quite an obvious opposite trend. NDVI indicates whether or not the target region being observed is covered by vegetation, while coherence is used to describe changes in backscattering properties Table 2 Relevant statistics of MTI results for the ancient slope failure reactivated in 2020 (Fig. 1). Parameters for the PSI and SBAS results (Fig. 4)   and similarities between radar echoes (Wang et al. 2009). The high NDVI values occur in summer, when the area is covered by more vegetation and thus the coherence becomes low. In contrast, less vegetation and thus, less volume scattering in winter results in higher coherence values for that season. An interesting result from Fig. 11a is that coherence drops to its lowest values ( ∼0.22) over the past six years before the landslide failure. There are another two minimum values for June 2015 and July 2016 in the coherence plot. However, these minima are caused by long temporal baselines between SAR image pairs during these periods. Figure 11b shows the comparison of NDVI values from three different satellite collections pronounced in Sect. 3.3. Due to the limited temporal resolution, NDVI values generated from Landsat-8 could not reveal promising results for detailed monitoring of vegetation dynamics. In contrast, NDVI time series generated from MODIS and S2 show good temporal correlation and agreement.  Fig. 7 Moreover, the two NDVI time series show some declines before the failure, i.e., MODIS-derived NDVI shows two declines of approximately 50 and 10 days before the failure; whilst the one derived from S2 indicates a drop of approximately 20 days before the failure. Reasons for these drops are elaborated in more detail in Discussion.

Discussion
This study has shown the great potential of applying highresolution optical and radar satellite remote sensing data and related techniques for the quantitative multi-temporal assessment of surface kinematics related to the 2020 Aniangzhai landslide failure in the mountainous region of Danba County. Generally, optical and radar remote sensing are two methods with their own advantages and weaknesses that can be used to monitor different stages or types of landslides, e.g., from initial slow creep motion to accelerated stage. These methods are complementary to each other and are exploited together in this study.
Based on the optical remote sensing data, the dynamics of the presented cascading events leading to Aniangzhai landslide failure could be clearly observed in their consecutive steps, allowing a comprehensive understanding of the resulting disaster chain. The failure with large displacement as shown in Fig. 3 is beyond the detective capability of the InSAR technique ). Our optical results are similar to the deformation vector distribution of reactivated deposits obtained by (Zhao et al. 2021). However, (Zhao et al. 2021) presented the results from tens of monitoring points from 22 June to 12 July, while in our study we derived the horizontal deformation of the whole failure area using optical remote sensing data between 15 June and 24 June. They also concluded that the Aniangzhai landslide slipped down as a whole; as the movement did not change the microtopography, rather the relative positions of points. Our fieldwork supports this observation, e.g., the houses on the failure body were only cracked but not collapsed (Fig. 2). The flooded areas, as well as the sediments of debris flows, are well demonstrated in the Planet and S2 images, indicating the direct cause and the sources of this event. The main direction of horizontal displacement from Planet optical images also provides a guidance for choosing the 6-year descending dataset in MTI analysis with a better observation geometry. Moreover, the results obtained from optical images of the June 2020 failure indicate larger deformation on the lower part of the slide compared to its middle and head parts in the large failure zone (Fig. 3). The explanation for this result could be that the foot of the landslide is closer to the river and is influenced by both water release from the upstream dam and barrier lake to relieve flood pressure. All these factors make the slope more vulnerable to undercutting and erosion, resulting in the large failure on the lower part of the slide.  Table 3. (b) The daily and cumulative precipitation from 2018 to 2020 Table 3 Increments of precipitation and the corresponding increments of LOS displacement of points P1-P3 from mid-April to mid-June in 2018-2020 As for the undercutting effects and erosion on the toe of the landslide body, debris flows from the Meilong valley caused by higher rainfall in summer 2020 played a vital role. This is consistent with the conclusions as revealed in (Yan et al. 2021). They suggest that the continuous rainfall in 2020 increased groundwater content and reduced the stability of loose sediments in Meilong valley. Then the debris flows formed in Meilong valley washed away loose sediments and bedrocks in the Xiaojinchuan River, a barrier dam was formed just under the Aniangzhai landslide toes and eventually the overflowed current from the barrier dam washed and eroded the foot of the ancient landslide body. By applying radar remote sensing and MTI techniques, temporal and spatial variability in the kinematics of the Aniangzhai landslide from 2014 until the 2020 failure could be comprehensively investigated ( Fig. 8 and Table 2). Mean LOS deformation rates during pre-disaster stage clearly identifies instability of the landslide, with the largest deformation rates higher than 50 mm/yr in Zone II. The deformation rates, however, change spatially in the entire slope with points located on larger slope angles, i.e., T3 and T4 in Fig. 5, showing higher LOS displacement rates compared to the other points. Zhao et al. (2021) identified the Aniangzhai slope as showing characteristics of a landslide with a constant deformation state that requires certain prevention measures before entering into next phase of catastrophic failure. Our multi-temporal interferometric results confirm this observation by Zhao et al. (2021) and further show that the long-term displacement rates before the June 2020 failure were not constant; rather, they changed over time. Influenced by above-average precipitation in summer and the 3-year wet period that followed a relative drought year in 2016, we observe that the landslide moved in 2018-2020 approximately 14% faster than in 2014-2017.

Increment
To better investigate the role of precipitation in influencing the kinematics of landslides, the statistics on the 2018-2020 period from mid-April to mid-June are analyzed (Table 3). This clearly shows how temporal changes in precipitation are correlated with the kinematics of motion of points in Zone I (Fig. 7). Several InSAR studies have shown that ancient landslides reveal instabilities or even precursors in the form of accelerated creep before the failure (Teshebaeva et al. 2015;Handwerger et al. 2019;Ao et al. 2020), but the source of acceleration could be different. In some cases, e.g., Teshebaeva et al. (2015), small creeps and accelerations have been correlated well with the increasing seismicity. We have checked the Chinese Earthquake Catalog and searched for earthquakes around Aniangzhai with a radius of 20 km over the last 1 year before the failure. Results, however, show no big earthquake ( > 2.0 ) occurred in the region; the nearest seismic activity during this period was around 25 km away from Aniangzhai in the southwest direction on 9 January 2020 (M 1.8 and a depth of 10 km). Therefore, we exclude tectonic forces as the source of the accelerated creep in Zone I.
It is worth noting that deep-seated landslides such as Aniangzhai cannot directly react to rainfall, since the changing of groundwater conditions towards activation of the sliding plane requires some time until the surface runoff has been infiltrated to a certain depth (Iverson 2000;Vallet et al. 2016). In the normal non-flooding seasons, toe erosion of the landslide should be a constantly ongoing process as well, leading to a backward propagation of deformation in the upslope direction until the time of failure (Teshebaeva et al. 2015;Leroueil and Locat 2020). In this process, rainfall does not play a direct role, but as an indirect one usually occurring with a lag in time (Iverson 2000;Teshebaeva et al. 2015;Haghshenas Haghighi and Motagh 2016;Vallet et al. 2016). In our case, however, results show that the lag time in the head is very short as the acceleration there occurs almost at approximately the same time when the precipitation increases. More, in-depth geophysical analysis will be needed to investigate the reason behind this short time lag between rainfall and motions in Zone I of this deep-seated landslide. Besides, our results of precipitation analysis also reveal that the formation of debris flows in the Meilong valley could trace back to April 2020 ( Fig. 6 and Fig. 9), when rainfall was approximately 30% higher than the mean values of the last 20 years. This is consistent with (Chen et al. 2005), that the formation of debris flow in this area, i.e., Danba County requires a longer preparatory phase of increased precipitation before a larger rainstorm eventually triggers the onset of the debris flow. All these observations suggest that the ancient Aniangzhai landslide was already active, which eventually failed partially following the undercutting effect in 2020.
The MTI analysis also provides a good basis for investigating evolution and kinematics of motion in different parts of the Aniangzhai landslide. The smoothing method developed in Section 3.4 helped us to properly smooth the data in order to detect the long-term and transient deformation without losing significant information from the data. These results have important implications for developing an early warning system for the Aniangzhai landslide and highlight that InSAR techniques can be used as an operational monitoring system in Aniangzhai to track progressive deformation and potential release areas in near real time in order to mitigate hazards associated with landslide failure (Hu et al. 2018;Carlà et al. 2018Carlà et al. , 2019Ao et al. 2020). Future works should focus on comparing the performance of moderate resolution S1 images with higher-resolution SAR images from missions like TerrasAR-X or CosmoSky-Med to better investigate the potential and existing limitations in Sentinel-1 data for landslide analysis Bovenga et al. NDVI analysis using MODIS and Sentinel-2 data reveals some interesting patterns, which may have the potential for early landslide warnings in the Aniangzhai landslide. As expected, the NDVI values increase while the coherence decreases due to the corresponding increases in the volume scattering. This behavior has been illustrated well in Fig. 11a, where a negative seasonal correlation is clearly observed between the coherence and NDVI values. Despite having different spatial resolutions, time series of NDVI from MODIS and S2 show good consistency with each other (Fig. 11b). Two interesting drops are observed in NDVI retrieved from MODIS data before the 2020 failure; one in May, and the other in June 2020, i.e., 50 and 10 days before the large failure. For the S2 data, a single drop is observed, i.e., 20 days before the failure. To investigate whether these signals are related to the behavior of the Aniangzhai landslide or whether they occur seasonally and independent of the slope motion, we plotted the NDVI values (MODIS) in 2020 against its historical values of the 2017-2019 period in Fig. 12a. Examples of LOS displacements in 2020 from points selected arbitrarily from Zone I are also shown in Fig. 12a   moisture, surface erosion, plant degradation and slope deformation Farrar et al. 1994;Liu et al. 2015;Jacquemart and Tiampo 2021). One explanation could be that the drops are related to errors in data production. As mentioned in Section 3.3, the MODIS dataset is generated with adjustment through BDRF, which might contribute to interpolation errors. Alternatively, the drops could be related to changes in vegetation structure before the failure as it occurs at approximately the same time, when a distinct acceleration in landslide active deformation is seen; no similar behaviours were observed in the historical values. Unfortunately, the S2 data do not have a dense temporal coverage for the first drop to be used as the validation. However, based on the above discussion and our results, we do not exclude the interpretation that the increase in rate of active deformation or the occurrence of small landslides before the main failure, as observed in seismic data (Yan et al. 2021), could have altered the scattering properties and vegetation structure in the landslide region, in turn leading to lower NDVI values. Similar observations of shallow soil erosion and vegetation damage near fault zones and river networks have been reported elsewhere (Gan et al. 2019;Jacquemart and Tiampo 2021;Geitner et al. 2021). As for the post-failure behaviors of NDVI time series, signification changes are expected compared to the pre-failure situation ( Fig. 11 and 12a). Similar observations of post-event behaviors can be found in (Behling et al. 2016). The structure of vegetation over the slope area could be greatly influenced by the large failure, and would last for a certain period in the future.
In a similar manner as above, we also investigated the changes in the time series of coherence values for the Aniangzhai landslide. A recent study indicated the linear relationship between NDVI and coherence (Bai et al. 2020). As longer temporal baselines could lead to lower coherence values, only SAR images with a 12-day temporal baseline were used here. However, as shown in Fig. 12b, the coherence-based results are not very promising. A possible reason is that coherence is usually already low during summer in this region, since this time of the year is characterized by maximal vegetation growth. Therefore, it seems difficult to accurately discern a further drop in the time series due to low coherence and high uncertainty. Nevertheless, we can still observe a relatively stronger declining trend for coherence in 2020 compared to the previous years, which might be related to active deformation and changes in volume scattering of the landslide. Unfortunately, there are no available Sentinel-1B data for this study area. Otherwise, we could study the coherence from 6-day image pairs to determine whether the drop in coherence caused by slope destabilization would be more accessible. It can be assumed that, if performed in winter, coherence analysis might have resulted in a better performance for ongoing slope activation with less vegetation compared to the analysis in this study.
The stretch to an EWS is hypothetical at this stage, and additional case studies are needed to further analyze the key factors in changing NDVI and coherence values. We suggest that the consideration of both parameters might lead to possible observations of signs of slope activation. With more experiments in the future, these results might contribute to a potential EWS for landslide hazards.

Conclusion
This paper focused on exploiting multi-sensor remote sensing technology to investigate the June 2020 Aniangzhai slope failure and the active deformation prior to the event since late 2014. Crosscorrelation analysis of high-resolution optical data from Planet provides detailed information about the spatial pattern of slope kinematics. Moreover, the undercutting effects on the toe of the landslide body, which played a vital role in the toe erosion and reactivation of this ancient landslide body, are also clearly visible in the optical data. The toe erosion was triggered by overflow of a dammed lake, created due to heavy rainfall and the resulting debris flows coming from the Meilong valley to the Xiaojinchuan River, and was influenced also by the discharge of the surplus water from a nearby hydropower station to reduce the flood pressure. Complementary analyses using multi-temporal SAR satellite remote sensing shows that the Aniangzhai landslide was not dormant. Rather, it was already active before the failure, with a maximum LOS displacement rate of around 38 mm/yr in 2014-2017, reaching approximately 55 mm/yr in 2018-2020. Our findings indicate that not the whole landslide body was subject to accelerating creep before the June 2020 failure; rather, only the points situated on the upper parts of the landslide failure sustained pronounced acceleration of the creep starting in spring 2020. As a result, the time series of displacements derived from these points could be utilized to forecast the potential window of failure. Moreover, we observed the sign in which an acceleration of creep on the head part of the failure region and a decrease in NDVI values took place almost at the same time, opposite to the prevailing trends in this area. We discussed the likely causes to interpret this phenomenon and suggested that this sign may be regarded as a parameter to be integrated into an EWS. With more case studies in the future, the methods proposed in this paper can be utilized under the framework of multidisciplinary DRR for a comprehensive analysis of the cascading event chain influencing the instability of the ancient landslides.