InSAR Study of Landslides: Early Detection, Three-Dimensional, and Long-Term Surface Displacement Estimation—A Case of Xiaojiang River Basin, China

Landslides, a major natural geohazard, obstruct municipal constructions and may destroy villages and towns, at worst causing significant casualties and economic losses. Interferometric Synthetic Aperture Radar (InSAR) technique offers distinct advantages on landslide detection and monitoring. In this paper, a more systematic workflow is designed for InSAR study of landslides, in terms of three levels: (i) early detection on regional scale, (ii) three-dimensional (3D) surface displacement rates estimation on detailed scale, and (iii) time series analysis on long-term temporal scale. The proposed workflow is applied for landslide research on the Xiaojiang River Basin, China, using ascending and descending Sentinel-1 images acquired from March 2017 to May 2019. First, the landslide inventory has been mapped and updated using InSAR stacking method, supporting geohazard prevention on a regional scale. A total of 22 active landslides are identified, ranging from medium to super large scale. Compared with the existing inventory, three unrecorded landslides are newly detected by our approach, and five recorded landslides are detected significant expansion of their boundaries. Then, specific to a detected landslide, Baobao landslide, a Total Least Squares– Kalman Filter-based approach is presented. Two outcomes are provided for further spatial-temporal pattern analysis: 3D displacement rates, providing an intuitive insight on the spatial characteristics and sliding direction of landslide, which are analyzed to deep the understanding of its kinematic mechanism, and long-term time series, which contribute to deduce the dynamic evolution of landslide, presenting benefits in landslide forecasting.


Introduction
As a major natural geological hazard, landslides interfere with the transport net-works and infrastructure development or may hold back the expansion of urbanization [1][2][3], and at worst, destroy villages and towns, causing significant casualties and eco-nomic losses [4][5][6][7]. There is a need for increased research on landslides, especially in terms of three levels, i.e., early detection on a regional scale, estimation of three-dimensional displacement rates on a detailed scale, and retrieval of displacement time series on a long-term temporal scale.
The detection of unstable slopes is an essential prerequisite for early warning and geohazard prevention [8], because of their extensive distribution and concealment. Interferometric Synthetic Aperture Radar (InSAR) performs quantitative remote sensing for small displacements over wide areas, with the advantages of centi-or even millimeter precision, high spatial-temporal resolutions, and all-day, all-weather working capabilities [9]. Different InSAR techniques [10][11][12] have demonstrated their distinctive capabilities in regional landslide detection, achieving valuable examples in Italy nationwide [13,14], the Alps [15,16], California [17], and the Jinsha River Catchment, China [18,19]. For instance, Novellino et al. [20] utilized the intermittent Small Baseline Subset (SBAS) to analyze the slope stability in northwestern Sicily, Italy; Bouali et al. [21] combined the Permanent Scatterers Interferometry (PSI) results with the landslide inventory to detect extremely slow (<16 mm/year) landslides across the Palos Verdes Peninsula, California; Liu et al. [19] employed InSAR stacking on Scan Synthetic Aperture Radar (SAR) images, achieving large-scale geohazard identification along the Jinsha River Corridor. Furthermore, the performance of different InSAR techniques for landslide detection has been compared in detail [22], revealing that the stacking method performs better on poor-coherence areas with the advantages of technical requirements and low computation labor. However, conventional InSAR techniques measure only one component of the surface displacement, i.e., in the line-of-sight (LOS) direction of satellite. This poses a challenge to interpretation and detection of landslides, and more than that, narrows the understanding of their displacement patterns and kinematic mechanisms.
Three-dimensional (3D) displacement rates provide intuitive insights on the spatial characteristics and sliding direction of landslides, which can be contributed to kinematic mechanism analysis [23]. To date, various strategies have been developed to estimate 3D displacement rates based on InSAR observations. The most straightforward strategy is to integrate multi-track InSAR observations acquired with different imaging geometries [24,25]. However, since most SAR satellite orbits are near-polar, this strategy may be applied solely for some specific cases in high latitude regions with higher azimuth angle diversity. Another strategy is to combine InSAR measurements with outputs from either offset tracking [26,27], Multiple Aperture Interferometry (MAI) [28,29], or other ancillary sources, e.g., Global Navigation Satellite System (GNSS) [30][31][32]. Moreover, a strategy is to impose a priori information as constraints on deformation process studied, which is typically applied for mining subsidence [33][34][35] and glacier or landslide motions [36][37][38][39][40][41]. Specific to landslide studies, a classical constraint is the surface-parallel flow model, i.e., movement along the steepest downslope, which has been accepted in many representative landslide cases, e.g., the Slumgullion Landslide, USA [38] and the Jiaju Landslide, China [39]. However, this constraint seems too strict for most landslides except translational landslides [40], since it does not allow consideration of the difference between ground surface and slip surface, forcible ignoring potential errors in model design. Overall, it is worth noting that the above strategies are typically developed to measure single deformation episodes with large gradient, e.g., earthquakes, fasting landslides, and volcanic activities. There are few studies on 3D displacement estimation of slow-moving landslides [40,41], particularly on a long-term temporal scale.
In contrast, as the entire landslide process can last from months to years, long-term time series of landslide movement is of great significance to deducing the dynamic evolution, presenting particular benefits in landslide forecasting and risk management [42][43][44]. Some efforts have been devoted for extending aforementioned strategies on temporal scale, e.g., the Pixel-Offset SBAS technique [45], the Multidimensional SBAS (MSBAS) technique [40,46], the Minimum Acceleration (MinA) approach [47], and the approach based on the Kalman Filter (KF) [48]. Among them, MSBAS and MinA are developed for simultaneous post processing of hundreds to thousands of individual interferograms. They regularize the rank-deficient and ill-posed inversion system by imposing further condition (e.g., Tikhonov regularization [40] or the assumption that displacement time series have minimum acceleration [47]), and then solve the 3D displacement time series utilizing the Singular Value Decomposition (SVD) method. In contrast, the KF-based approach is a sequential processing, which combines multi-source InSAR observations at the present epoch in a Kalman filter model with displacement information determined from previous interferograms. This approach allows the displacement time series dynamically update with an improved temporal resolution and limited computational cost. Its performance has been demonstrated in researches on volcanics and faults motions [48,49], but up to now, this approach is rarely applied to the retrieval of the time series of landslide displacements.
In this study, a more systematic workflow is designed for InSAR-based landslide study, including three levels, i.e., (i) early detection on regional scale, (ii) 3D displacement rates estimation on detailed scale, and (iii) time series analysis on long-term temporal scale. First, a standard InSAR stacking procedure is employed to investigate the distribution and activities of landslides over study area. A relatively complete landslide inventory is mapped and updated for geohazard prevention on regional scale. Then, specific to a detected landslide, 3D displacement rates and long-term time series are further estimated for the spatial-temporal pattern analysis in more detail. To this end, a Total Least Squares-Kalman Filter (TLS-KF)-based approach is presented, combining advantages of Total Least Squares (TLS) and advantages of the Kalman Filter (KF). The outcomes can provide intuitive description of sliding direction and displacement characteristics of landslides, deepening the understanding of its kinematic mechanism and dynamic evolution.
Here, the proposed workflow is applied for landslide research on the Xiaojiang River Basin, Eastern Yunnan, China, which is the most serious area of soil erosion and geological disaster in the upper reaches of Yangtze River [50]. This paper is organized as follows. A brief description of the study area, as well as the datasets, is given in Section 2. The proposed workflow is presented in Section 3, with its theory and processing elaborated in detail. In Section 4, a new landslide inventory has been mapped and updated by our approach, with additional information of LOS displacement rates. Moreover, one of the detected super large landslides, the Baobao landslide, is further studied in Section 5. The three-dimensional and long-term displacements are estimated by the TLS-KF-based approach and its spatial-temporal displacement pattern are discussed.

Study Area
The study area is situated in the Xiaojiang River Basin, Eastern Yunnan, China ( Figure 1), north from Menggu and south of and nearby Dongchuan, covering an area of approximately 2314 km 2 . It is the most serious area of soil erosion and geological disaster in the upper reaches of the Yangtze River [50], as well known by the Natural Museum of Landslides [51], under combined actions of multiple factors, e.g., complex natural geography, lithology, active tectonic movements, special climate, and human activities. The elevation in most parts of study area is higher than 2500 ma.s.l., reaching more than 4000 ma.s.l. in some regions. Due to violent river downward eroding, valleys feature strong V-shaped topography with height differences up to 3000 m, resulting in slope angles of greater than 25 • in most slopes, thus providing favorable geological conditions for the occurrences of landslides. The tectonic setting is conditioned by a series of nearly north-south trending faults, including the Xiaojiang Faults, the Pudu River Faults, and the Yuanmou Faults (marked as F1, F2, F3 in Figure 1), resulting in frequent seismic activities. There have been approximately 41 earthquakes of Mw ≥ 3.0 in the study area and its surroundings in the past ten years (2009-2019), including three stronger earthquakes greater than Mw = 5.0. Generally, rock masses in active tectonic zones are more prone to develop joints and fissures and to fracture [52], presenting poor stabilities. Therefore, landslide susceptibility tightly correlates with the distribution and activities of faults. Specific to the study area, tectonic movements are mainly affected by the north section of the Xiaojiang Faults, which is the most active sinistral strike-slip fault zone in southeast margin of the Qinghai-Tibet Plateau [53].  Figure 2a shows the geological setting map of the study area. The outcrops are composed of almost all strata from Jurassic to Mesoprotero-zoic, in detail, including middlelower Jurassic, Permian, Triassic, Carboniferous-Devonian, Silurian, Ordovician, Cambrian, Paleozoic, Sinian, and Kunyang Group strata. They mainly include basic volcanic rocks, dolomites, limestones, phyllites, mudstones, and sandstones, and diabase in local regions. The lithology of strata features as soft and hard rocks alternatively distributed, leading to poor compactness and low shear strength [54], and it is an inherent factor in the full development of landslides.
The climate belongs to subtropical plateau monsoon climate [55]. The annual temperature and precipitation are approximately 15 • C and 800 mm, respectively. As a typical hot-dry valley, the climate features relatively high temperature and low humidity and precipitation concentration (88% of precipitation is concentrated during June to October). The former leads to strong physical weathering on the slope surfaces; the latter leads to rain infiltration decreases shear strength of materials, thus driving landslide movements. For instance, Jiangjia Gully debris flow (in Figure 2a, mapped with white curve) is induced by strongly weathered limestones under concentrated heavy rainfalls. In addition, the study area is rich in mineral resources, among which Dongchuan is one of six famous copper productions [56]. Mining activities lasting more than 2000 years have broken the balance of groundwater and destroyed the initial stress states of some slopes, promoting the occurrence of landslides. As noted, multiple factors work together to result in the extensive distribution and strong activity of large-scale landslides in the Xiaojiang River Basin. Potential landslides threaten the safety of cities and towns as well as traffic lines (e.g., Provincial Highway S23). Moreover, the landslides could dump large amounts of rock/soil mass into the Xiaojiang River, largely endangering the operation of Baihetan giant hydropower station on downstream of its mainstream. China Aero Geophysical Survey & Remote Sensing Center for Land and Resources (AGRS) has mapped the distribution of landslides in this area through optical remote sensing (see Figure 2) and verified them by field investigation. However, the detection only based on optical images may omit some active landslides, especially those without significant morphological features (e.g., scarps, sliding masses, and bulging toes). More importantly, complete investigations of landslides in terms of 3D displacements, spatial pattern, and kinematic evolution are absent [41].

Datasets
The detection of active landslides is an essential task for disaster prevention. More than that, the estimation of 3D landslide displacement rates and long-term time series is of great significant to deepen the understanding of its kinematic mechanism and dynamic evolution. To these ends, a total of 117 SAR images are collected from Sentinel-1 ascending Track 26 and descending Track 62, covering the period from March 2017 to May 2019. The spatial coverage of the SAR images provided by each track is illustrated in Figure 1. Other basic parameters of the SAR images are summarized in Table 1, in detail. Moreover, a 1-arc-second SRTM digital elevation model (DEM) is employed in the standard differential InSAR processing, and in the calculation of the LOS visibility maps of ascending and descending tracks [57]. The holes (non-value pixels) in the DEM are filled by bicubic spline interpolation.

Methodology
The InSAR stacking method is employed to detect active landslides in the study area. A new approach combining Total Least Squares (TLS) and Kalman Filter (KF) is presented to estimate 3D displacement rates and long-term time series of detected landslides. The workflow is shown in Figure 3, which consists of three steps, as follows. Step 1. Generation of LOS displacement rate maps and detection of active landslides. First, each SAR dataset is handled with a standard differential InSAR procedure, mainly including co-registration, interferograms generation, phase unwrapping, errors correction, and quality screening. Then, high-quality unwrapped interferograms are outputted for the following processing, stacking interferograms [12]. Relative InSAR processing is elaborated in Section 3.1. As outputs, the LOS displacement rate maps are derived from ascending and descending datasets, with a spatial resolution of 15 m and in a Universal Transverse Mercator (UTM) geodetic system. These rate maps are superimposed onto optical images, combined with morphological features (e.g., scarps, sliding masses, and bulging toes), to interpret and identify active landslides. Finally, suspected landslides detected from each dataset are mosaiced [21] to map a relatively complete landslide inventory.
Step 2. Estimation of 2D/3D displacement rates and analysis of spatial pattern. Specific to a detected landslide, if the north-south displacement component can be neglected, the 2D displacement rates are derived by fusing ascending and descending unwrapping interferograms [40]. Otherwise, the 3D displacement rates are estimated in total least squares sense by imposing a relaxed constraint, i.e., movement along slip surface, which is more general than the ground-surface parallel flow constraint [40,58]. The methodology is presented in Section 3.2. Based on the outputs, the optimal sliding direction of landslide is determined for the analysis of the spatial displacement pattern.
Step 3. Retrieval of long-term displacement times series and temporal evolution. Ascending and descending unwrapped interferograms are fused to retrieve the long-term displacement time series in the optimal sliding direction based on the Kalman Filter [48]. The methodology is described in Section 3.3. Then, the temporal evolution of landslide is investigated based on the OSD displacement times series and the possible driving factor, precipitation.

InSAR Processing
A standard differential InSAR procedure is employed to process all SAR images, as follows. First, single look complex (SLC) images from ascending and descending Sentinel-1 are co-registered, separately, with the aid of auxiliary data (e.g., precise orbit files and SRTM DEM). Second, the small baseline subset (SBAS) combination strategy [10] is adopted to generate high-quality interferograms with less temporal and spatial decoherence. Specific to our datasets, the temporal baseline and spatial baseline are restricted to 40 days and 200 m, respectively. In processing, the multi-look factor is set to 4 (4 pixels in range and 1 pixel in azimuth), resulting in interferograms characterized by a ground resolution of approximately 15 m × 15 m. Similar resolution data have proven their performance in regional landslide mapping in many previous works [59,60]. For a given interferogram, the minimum cost flow algorithm [61] is applied to unwrap phases. A semi-adaptive procedure is performed to correct various items of phase errors. The unwrapping error is corrected by either a phase compensation algorithm or re-unwrapping at a high coherence point [19]. The error related to residual baseline is estimated and reduced based on a quadratic polynomial model [62]. The error caused by DEM is corrected using the method presented in [63]. The atmospheric phase is corrected by the spatially variable power law tropospheric correction technique [64], preliminarily. It is worth noting that the climate in the study area features complex and ever-changing local micro-climate, so the atmospheric phase needs further correction in the follow-up. Finally, a total number of 186 and 131 highquality unwrapped interferograms are screened from ascending and descending Sentinel-1 datasets, respectively. Their spatiotemporal distributions are illustrated in Figure 4. Then, the stacking interferograms method [12] is employed to retrieve the displacement rate of each dataset in LOS direction. Equation (1) shows its kernel, a linear combination of interferograms (ϕ i ) with the temporal baseline (∆T i ) as weights: While the atmospheric phase is random for independent interferograms, the related error was effectively corrected by stacking, and the signal-to-noise ratio was improved by a factor of √ N [65]. Specific to our datasets, compared with the descending dataset, the ascending dataset generates a larger number of interferograms (approximately 1.5 times) for stacking, thus inferring that atmospheric errors will be corrected more effectively and a higher accuracy will be achieved. In addition, the LOS displacement rate maps are geocoded into UTM grid system (e.g., World Geodetic System 1984), preparing for the integration of multi-source InSAR results to detect active landslides. Moreover, invalid measurement scatterers in flat areas and in areas presenting low degree of LOS visibility (e.g., layover and shadow) are masked to avoid possible misjudgments.

Total Least Squares Estimation of Three-Dimensional Landslide Displacement Rates
If only two independent SAR datasets are available and the north-south displacement component is non-negligible, a prior model should be imposed as the constraint on landslide movement to reduce the freedom degree for inversion of 3D displacement rates. A classical constraint is the surface-parallel flow model [58], which has proved its applicability on translational landslides. However, this constraint is too restrictive for most landslides, as their mechanisms and types are more complex in reality, and their buried slip surfaces are usually undeterminable to be parallel to topography. To generalize further, the constraint is relaxed for non-translational landslides, as follows: where V U , V E , and V N are the displacement rate parameters in vertical, east-west, and north-south direction, respectively, ∂(H+∆h) ∂x and ∂(H+∆h) ∂y represent the first derivatives of slip surface in the east and north directions, respectively, H is the elevation of the ground surface, which was derived from a smoothed external DEM [66], and ∆h is the elevation difference between the ground surface and slip surface.
Under the relaxed constraint, i.e., movement along slip surface, the observation model is constructed as follows: where L is the observation vector in LOS direction derived from ascending and descending tracks, S E , S N , and S U are the east-west, north-south, and vertical components of the LOS unit vector S [40], and G represents the time intervals between consecutive SAR acquisitions. It is worth noting that the slip surface is roughly deduced using the strain model [67] and/or field exploration [68], and is sometimes even inaccessible; thus, ∆h is unknown with considerable error. Furthermore, whereas ∆h is generally small enough relative to H, the terms composed of ∆h are regarded as Errors-In-Variables (EIV) and the observation model is converted into the EIV frame, as follows: where X = [V E , V N ] T is the unknown parameter vector,Ĝ = G S E + ∂H ∂x S U , S N + ∂H ∂y S U is the updated design matrix, and E G = S U G ∂∆h ∂x , ∂∆h ∂y and V represent errors of model design and observation vector, respectively, with stochastic properties characterized by: where vec indicates the operator that stacks one column of a matrix under the previous column and Σ V and Σ G are the variance-covariance matrixes of observation error and design error, respectively. It is well known that errors of model design are ignored in the classical Least Squares (LS) estimation; thus, the solutions in LS sense are no longer optimal. To adjust to the EIV observation model, the Total Least Squares (TLS) estimation takes into account both errors in observations and model design, resolving based on the criterion of minimizing the objective function with Equation (6). Their performances have been proved in some similar problems, e.g., 3D coordinated transformation [69], in GNSS: Furthermore, Σ G can be approximately expressed in form of the Kronecker-Zehfuss product (⊗), as Σ G := I 2 ⊗ Σ 0 , where I 2 is the identity matrix and Σ 0 is the variance matrix of each column of E G , composed of variance of elevation difference (σ 2 ∆h ), as Σ 0 := σ 2 ∆h * S 2 U * diag G(1, 1) 2 . . . G(N, 1) 2 . Thus, a fast convergent iterative algorithm presented in [70] is adopted to obtain the optimal solution in the TLS sense (see Table 2). Table 2. Equations of iterative total least squares estimation presented in [70].

Kalman Filter for Long-Term Landslide Displacement Time Series
The retrieval of long-term time series of landslide displacement is particularly meaningful to the deduction of landslide dynamic evolution. Toward this end, the Kalman filter, an effective approach to fuse multi-source SAR observations obtained at different times [48], is employed to estimate the landslide displacement time series for each pixel in the optimal sliding direction (OSD). The OSD is determined in Step 2 using the TLS estimations of 3D displacement rates, V TLS,U :V TLS,E :V TLS,N = a:b:c. The observation model and the state model of the Kalman filter are constructed for a sliding pixel at acquisition date t k , as follows: where L k is the LOS observation, X k = d k v k T is the unknown state vector, d k and v k are the OSD displacement and velocity, respectively, H = aS U +bS E +cS N √ a 2 +b 2 +C 2 0 is the projection matrix, A = 1 ∆T 0 1 is the state transition matrix, ∆T is the time interval between t k−1 and t k , and V k and W k−1 are the vectors of measurement noise and process noise, respectively, with stochastic properties characterized by: where R k is the observation variance of the kth interferogram and Q k is the variancecovariance matrix of process noise, which can be recursed using Sage-Husa adaptive filter [71], as Equation (10). The OSD displacement time series can be dynamically estimated by the Kalman filter in near-real time. All relevant equations are listed in Table 3. Table 3. Equations of the Kalman filter presented in [72].
Time update ("predict"): Sage-Husa adaptive filter of Q k : where x − k and x k are the prior and posteriori state estimations, with P − k and P k being the corresponding variance-covariance matrixes, K k is the Kalman gain, f k = 1−b 1−b k is the forgetting factor, and b is a preset parameter, recommended to be set at 0.95-0.99.

Results: Line-of-Sight Displacement Rates between March 2007 to May 2019
Following Step 1, the LOS displacement rate maps are derived from ascending and descending Sentinel-1 SAR datasets, respectively, as shown in Figure 5. It is worth noting that the negative values (red color) indicate the movements away from the satellite, and the positive values (blue color) indicate the movements toward the satellite.  Measurement scatterers (MSs) with a total number of 9,584,080 and 9,179,644 are identified from the ascending dataset ( Figure 5a) and descending dataset (Figure 5b), respectively, thus producing an overall spatial density of approximately 4000 MSs/km 2 . These scatterers are selected by a fixed standard (a coherence threshold of 0.2 and an intensity threshold of −1 dB) and generally distributed on the roads, bridges, artificialities, and rocks/soils with sparse vegetations. Compared with the ascending track, the heading angle of descending track is more likely to cause severe geometric distortions on eastward slopes, resulting in a relatively sparse MSs for landslide detection in the study area. Moreover, the distribution histograms of MSs in stable regions are calculated for accuracy assessment, as shown in Figure 6. The mean values (red line) are lower than 0.1 cm/year, and the standard deviations (blue line) are approximately 1 cm/year. Refocusing on Figure 5a,b, the locations of active displacement regions are roughly consistent between ascending and descending results, while the displacement magnitude and their detailed extents are locally different. This is mainly attributed to their different LOS sensitivities to mass movements in each slope orientation. To complement each other, the ascending and descending LOS displacement rate maps are combined to map a relatively complete landslide inventory. Finally, a total of 22 active landslides, ranging from medium scale to super large scale, are identified, mapped, and labeled in Figure 5, and their detailed information are summarized in Table 4. These landslides are concentrated near the Xiaojiang River and some of them are re-active parts of historical debris flows.    Super large S1A. S1D Low Debris flow source 1 Notes: S1A and S1D represent ascending and descending Sentinel-1 SAR images, respectively. Newly detected landslides (No. [20][21][22] and recorded landslides with an obvious expansion (No. 14, [16][17][18][19] are highlighted in red. At this point, the landslide inventory of the study area has been mapped and updated by our approach with addition information of the LOS displacement rate, supporting the regional disaster prevention and risk control. However, for the detected landslides, their 3D displacement rates and long-term time series remain to be further studied, in order to deepen the understanding of their kinematic mechanism and dynamic evolution. In Section 5, the Baobao landslide (labeled No.14 in Table 4) is taken as an example case to retrieve its 3D displacement rates and long-term time series. Furthermore, the spatialtemporal displacement pattern of Baobao landslide will be discussed in the following sections.

Three-Dimensional Displacement Rates for Spatial Displacement Pattern Analysis
The Baobao landslide is one of detected landslides on super large scale, the extent of which is~0.76 km in width and~1.57 km in length, located on the east bank of Pucha River and in the Xiaojiang Fault zone. Its optical remote sensing image acquired in April 2019 is shown in Figure 7a, where polygons with different colors represent four Blocks B1-B4 of the landslide, which are divided according to the geomorphological analysis. Following Step 2, three-dimensional displacement rates are estimated in the vertical, north-south, and east-west directions during March 2017 to May 2019, as shown in Figure 7b-d, respectively.
It is worth noting that negative values (red color) indicate downward movements in Figure 7b, southward movements in Figure 7c, and westward movements in Figure 7d, and positive values (blue color) indicate northward movements in Figure 7c. The maximum displacement rates in the vertical, south-north, and east-west directions are more than −6.6 cm/year, ±11.0 cm/year, and −16.2 cm/year, respectively. The two most active movements are observed at the middle of block B1, i.e., Region R1, and lower-left edge of block B2, i.e., Region R2, with averaging displacement rates of 13.5 cm/year and 14.9 cm/year in the optimal sliding direction (OSD), respectively. Overall, the entire landslide presents a trend of westward and downward movement. As evidence in Figure 7c, northward displacements mainly occurred in block B1, and are sliding toward the northwest direction, while southward displacements mainly occurred in Region R2 and the lower half of block B3, and are sliding toward the southwest direction. A distinct sliding boundary can be clearly observed between block B1 and blocks B2-B3. For spatial displacement pattern analysis, displacement rates and elevation are extracted along two representative Profiles AA' and BB', as shown in Figure 8. The Profile AA' is along the terrain gradient descent, approximately parallel to the main sliding direction of Blocks B1-B2, and the Profile BB' transversely passes through Blocks B2-B3 at similar elevation (see Figure 7d). In Figure 8a, the sliding movement in OSD of Blocks B1-B2 (along AA') presents a negative correlation with the elevation, that is, displacement in the front part (Block B2) is larger than the posterior part (Block B1). This evidence indicates that Blocks B1-B2 are developed under a pull-type mechanism [73], which can be preliminarily verified by geomorphological features (e.g., tensile cracks and scratches at back edge of B1) presented in Figure 7a. In Figure 8b, evidence illustrates that Block B3 is more stable than Block B2 (along BB'). Moreover, as evidenced in Figure 7b-d, 3D displacements in Blocks B3-B4 are significantly smaller than its above Block B1 and present positive correlations with elevation, which suggest that Blocks B3-B4 is very likely to be deformed with a thrust-type mechanism [73], under the continuous extrusion of Block B1. Furthermore, Figure 9 shows the optimal sliding direction (OSD) for each pixel of landslide, determined using the estimated 3D displacement rates. It worth noting that OSD vectors indicate that Block B3 is sliding towards the southwest, imposing a direct threaten on the Baobao Village. Under the pressure of Block B1, the sliding movement of Block B3 will be more active. It is strongly recommended to enhance the monitoring, prevention, and control of the Baobao landslide, in order to avoid potential losses of life and property in the surroundings.
Overall, through Figures 7-9, an intuitive description of the sliding direction and displacement characteristics of each block is provided for the study of Baobao landslide from a three-dimensional perspective. The spatial pattern and deformation mechanism is analyzed on a detailed scale. In particular in Figure 9, the OSD vectors superimposed on the optical image most intuitively exhibit the potential threat of landslide to surrounding villages, traffics, rivers, etc., playing a positive role in disaster risk management.

Long-Term Displacement Time Series for Temporal Displacement Pattern Analysis
In order to reconstruct the dynamic evolution of Baobao landslide, four points, P1-P4 (marked in Figure 7b), located in different blocks of the landslide are selected to retrieve the long-term time series of displacement in optimal sliding direction from ascending and descending Sentinel-1 images, according to Step 3.   Except for gravity, the primary driving factor, the acceleration of landslide movement may also be caused by some external environment factors, including earthquakes, heavy precipitation, and river and groundwater level fluctuations. Specific to the study area, two recorded earthquakes marked in Figure 1 occurred in 2009 and 2013, respectively, outside the temporal coverage of SAR acquisition. With this in mind, focusing on the correlation between displacement time series and monthly precipitations, the evolution process of Baobao landslide is analyzed in combination with the standard three-stage creep rupture curve model (refer to Figure 14a in [41]). From March 2017 to May 2018, the displacements of Points P1-P4 all present linear evolution trends. Landslide displacement during this period exhibits creep behavior of the secondary creep stage, i.e., constant state with a steady displacement rate. After June 2018, a significant signal of displacement acceleration is observed at all points P1-P4, which may be correlated with heavy rainfalls in flood season. Large amounts of rain infiltration can increase the pore-water pressures of slopes and, thus, reduce the shear strength of materials, resulting in an increase in sliding movement. The Baobao landslide has entered the tertiary creep stage. Then, with the arrival of dry season (from October 2018 to May 2019), precipitation gradually decreased and, thus, the Baobao landslide reached a new limit equilibrium.

Conclusions
In this paper, a full-scale workflow is designed for the current increasingly systematic and multi-level InSAR study of landslides, in terms of three levels: (i) early detection on a regional scale, in favor of landslide inventory and risk assessment, (ii) 3D displacement rates estimation on a detailed scale, enhancing understanding on kinematic mechanism and failure mode, and (iii) time series analysis on long-term temporal scale, presenting benefits in landslide forecasting. This workflow is applied for landslide research on the Xiaojiang River Basin, Eastern Yunnan, China, using ascending and descending Sentinel-1 images acquired from March 2017 to May 2019. A relatively complete landslide inventory is mapped and updated for geohazard prevention. New landslides are discovered, and known ones are re-evaluated on their activities and boundaries. Furthermore, three-dimensional and long-term displacement information of Baobao Landslide are estimated through the proposed Total Least Squares-Kalman Filter-based approach. Based on these outcomes, the spatial-temporal pattern is further analyzed to deduce the failure mechanism and current evolution stage.
Through this study, the potential of Sentinel-1 images in InSAR study of landslides has been fully explored and demonstrated. As a rare open resource, Sentinel-1 satellites provide a large amount of free valuable observations since their launch, greatly promoting the development of InSAR technology and its applications in landslide studies. However, Sentinel-1 data with C band are relatively sensitive to seasonal variation in vegetation coverage, often resulting in unsatisfactory performance on landslide monitoring in natural environments with lush flora. Moreover, medium ground resolution of Sentinel-1 data (~15 m) can hardly detect small to medium landslides, let alone to characterize motion mechanism in detail. These two points limit the capability of Sentinel-1 in the InSAR study of landslides to some extent. In contrast, ALOS-2 data with L band and high ground resolution (~3 m under ultra-fine mode) could have played a greater role in this filed. However, its unstable data acquisition results in that the actual observation interval is several times longer than the revisit period (14 days) in most regions. Previous experience indicates that stable data acquisition with high observation frequency is of great significance to closely monitor the landslide development and capture the omen of failure. A hot prospect of new and forthcoming SAR satellites (e.g., LT-1, Tandem-L), with longer wavelengths, higher resolution, and shorter revisit period, is predictable in InSAR study of landslides, focusing on two aspects: (i) regional early detection based on multi-source and multi-mode data integration and (ii) 3D displacement time series estimation and prediction combing high-resolution high-frequency observations and landslide models.