A Spatio-Temporal Monitoring Method Based on Multi-Source Remote Sensing Data Applied to the Case of the Temi Landslide

Abstract

It is challenging to monitor landslides due to their heavy concealment and the extreme destructiveness during the long development of landslides. Many landslide monitoring tools are somewhat onefold. In this paper, a comprehensive landslide monitoring method involving multiple factors from time-series multi-data sources is proposed. We focus on the changes in three aspects consisting of the vegetation condition, the surface deformation information and the landslide susceptibility. Firstly, the fractional vegetation cover of the landslide is extracted from optical remote sensing Gaofen-1 (GF-1) images using the dimidiate pixel model. Next, the surface deformation information of the landslide is derived from SAR remote sensing Sentinel-1A images applying the SBAS-InSAR method. Then, the landslide susceptibility based on GF-1, Sentinel-1A images and DEM data is computed using the analytic hierarchy process method. Finally, the spatio-temporal correlations of the vegetation condition, the surface deformation information and the landslide susceptibility are compared and interpreted. The Temi landslide is located along the Jinsha River and poses a high risk of blocking the river. Taking the Temi landslide as the study area, it is indicated from the results that the fractional vegetation cover, surface deformation information and landslide susceptibility reveal a consistency in the patterns of changes in spatial and temporal terms. As the surface deformation information improves, the status of the landslide vegetation also deteriorates and the landslide susceptibility becomes high, which indicates an increased probability of the creep and even the occurrence of landslides. In contrast, when the surface deformation information drops, the vegetation condition of the landslide becomes superior and the landslide becomes less susceptible, which means the likelihood of sliding declines. This study provides a new idea for a landslide monitoring method and potential way for natural disaster prevention and mitigation.

1. Introduction

China is prone to geological hazards, with landslides occurring particularly frequently. A landslide is a phenomenon in which part of a slope moves in shear along one or several faces within the slope under the influence of gravity [1]. Because of their suddenness, concealment, destructiveness and wide impact, landslides pose a huge safety hazard to lives and properties of people. A single landslide can cause significant damage, such as the large-scale Xinmo village landslide in 2017 [2] and the Baige landslide in the Jinsha River in Tibet in 2018 [3]. Therefore, the scientific risk assessment and timely monitoring before the occurrence of landslides can effectively reduce losses, which is of great significance for disaster prevention and mitigation.

In recent years, many landslides that have occurred and caused disasters are not catalogued as landslide hazards, but are in inaccessible areas, on dangerous terrain, and highly concealed. Early landslides are currently identified and monitored using long-term monitoring of landslide deformation. The direct measurement monitoring methods include the geodetic method, the close shoot method [4], the global navigation satellite system (GNSS) [5] and the Interferometric Synthetic Aperture Radar (InSAR). The geodetic method uses highly accurate optical and optoelectronic instruments to measure distance and angle. The instruments used in monitoring landslide surface displacements are generally levels, latitude and longitude instruments, total stations, displacement sensors and infrared rangefinders. This method can observe the absolute surface displacement, but it is susceptible to external conditions. The close shoot method uses the close-up surveyor to simultaneously image two fixed observation points at different locations and analyzes the changes in the corresponding locations for the purpose of deformation monitoring, which is easy for storing the measurement information, but vulnerable to weather. The GNSS method monitors the displacement and deformation by placing feature points at fixed locations and obtaining the change value in the feature points based on the difference between 3D coordinates at different periods. This technology is fast for collecting information, but it is difficult to directly reflect the continuous change information of the surface because of the point-by-point basis, and the signal is sensitive to interference.

InSAR technology was used to monitor landslide deformation in 1995. J. Achache et al. used the difference interferometric synthetic aperture radar (D-InSAR) to monitor the Saint-Etienne-de-Tinee landslide in southern France, and the results were consistent with the ground-based measurements [6]. In 2000, the Italian scholar Rizo also analyzed the deformation rate of landslide by the D-InSAR technique and agreed with the GNSS results [7]. In the same year, the persistent scatterer interferometric synthetic aperture radar (PS-InSAR) method was proposed by Ferretti and Rocca et al. in Italy [8]. In 2012, Bianchini S. et al. used the PS-InSAR technique to monitor potential slowly deforming landslides [9]. In 2002, Berardino et al. proposed the small baseline subsets interferometric synthetic aperture radar (SBAS-InSAR) technique [10]. In 2021, Guo et al. used SBAS-InSAR method to obtain the surface deformation characteristics of the Xiongba ancient landslide on the western bank of the Jinsha River [11]. InSAR technology has the advantages of a wide range, being able to be used all day and in all weathers, and high accuracy. It is widely used in surface deformation monitoring. However, there are some problems, such as low density of coherent points and uneven spatial distribution in areas with lush vegetation and steep terrain [5].

In addition to the above direct surface displacement monitoring, landslides can also be monitored indirectly using optical remote sensing. Optical remote sensing technology has characteristics such as non-contact, direct visibility and a large range. Guo et al. [12] identify the vegetation change information in optical remote sensing images before the occurrence of landslides and have found that the vegetation change has a significant spatio-temporal correlation with the landslide creep. This method makes up for the shortcomings of traditional methods in complex environments and saves manpower and resources effectively. However, landslide monitoring only through vegetation change information is relatively onefold. According to the above analysis, different monitoring methods monitor from different perspectives. Each has its own strengths and weaknesses, and each monitoring method is not comprehensive enough. Moreover, the landslide is a complex process of interdisciplinary interaction of multiple factors. It is necessary to find a new method to consider multiple factors from multiple perspectives with multi-source data for monitoring landslides under complex conditions effectively and economically.

The Temi landslide is located along the Jinsha River, which has a great risk of blocking the river. Once the giant landslide body slides, a dammed lake will be formed, which will seriously threaten the safety of lives and properties of people. This study focuses on the Temi landslide using multi-source data including time-series optical remote sensing data, time-series synthetic aperture radar (SAR) data, and digital elevation model (DEM). Considering multiple factors affecting landslides, we identify the vegetation change information on the landslide body, extract the surface deformation information and construct the landslide susceptibility model to comprehensively monitor landslides from spatial and temporal perspectives.

2. Study Area and Data

2.1. Study Area

The Temi landslide along the Jinsha River is under 99°3′6.91″ E and 29°58′31.35″ N, respectively. It is at the junction of Batang County (Ganzi Prefecture, Sichuan Province) and Mangkang County (Changdu City, Tibet Autonomous Region). The altitude is between 2530 m and 3150 m. Due to the constraints of terrain factors, the large highway nearby is G318 (see Figure 1). The front edge of the landslide is in the air and the two sides are cracked. The front edge is the steepest, followed by the gentle trailing edge and the middle part. The Temi landslide is about 1300 m in length and 450–800 m in width, with a total area of 750 thousand m² and a total volume of 60 million m³. The main landslide body is on the left bank with about 33 million m³, and about 14 million m³ of accumulation remains on the right bank [13]. This study focuses on the main landslide on the left bank. The steep terrain and vegetation cover in the Temi landslide make it difficult to reach, thus making the direct measurement monitoring tools impossible.

2.2. GF-1 Optical Image

Optical remote sensing technology has the advantages of wide range and non-contact. High spatial resolution images are needed for the single landslide. Gaofen-1 (GF-1) images from 2018 to 2021 are selected from the Land Observation Satellite Data Service Platform of the China Center for Resource Satellite Data and Application. The GF-1 image carries two cameras, a pomatic (PAN) camera with 2 m spatial resolution and a multispectral (MS) camera with 8 m spatial resolution. In order to avoid the influence of too large a time difference for the vegetation, images with similar time are selected in the vegetation growing season as far as possible. However, due to factors such as clouds and rain in the study area, the final images available are shown in

Table 1. GF-1 image information of the Temi landslide.

ID	Satellite	Sensor	Date	Spatial Resolution/m
1	GF-1	PMS1	2018-06-08	2, 8
2	GF-1	PMS1	2019-06-12	2, 8
3	GF-1	PMS1	2020-07-26	2, 8
4	GF-1	PMS1	2021-10-20	2, 8

.

2.3. Sentinel-1A SAR Data

SAR images facilitate all-weather observations. Sentinel-1A satellite launched in 2014 by the Global Monitoring for Environment and Security (GMES), carries the C-band synthetic aperture radar with a distance resolution of 5 m and a direction resolution of 20 m. In this study, the single look complex (SLC) images from 27 October 2017 to 11 November 2021 with interferometric wide swath (IW) mode, descending orbit direction and VV polarization are used, for a total of 62 images, as shown in

Table 2. Sentinel-1A image information of the Temi landslide.

ID	Date	ID	Date	ID	Date	ID	Date
1	27 October 2017	17	15 November 2018	33	16 December 2019	49	3 January 2021
2	20 November 2017	18	21 December 2018	34	9 January 2020	50	27 January 2021
3	14 December 2017	19	14 January 2019	35	2 February 2020	51	20 February 2021
4	7 January 2018	20	2019-02-07 7 February 2019	36	26 February 2020	52	16 March 2021
5	31 January 2018	21	3 March 2019	37	21 March 2020	53	9 April 2021
6	24 February 2018	22	27 March 2019	38	14 April 2020	54	3 May 2021
7	20 March 2018	23	20 April 2019	39	8 May 2020	55	27 May 2021
8	13 April 2018	24	14 May 2019	40	1 June 2020	56	20 June 2021
9	7 May 2018	25	7 June 2019	41	25 June 2020	57	14 July 2021
10	31 May 2018	26	1 July 2019	42	19 July 2020	58	7 August 2021
11	24 June 2018	27	25 July 2019	43	12 August 2020	59	31 August 2021
12	18 July 2018	28	18 August 2019	44	5 September 2020	60	24 September 2021
13	11 August 2018	29	11 September 2019	45	29 September 2020	61	18 October 2021
14	4 September 2018	30	5 October 2019	46	23 October 2020	62	11 November 2021
15	28 September 2018	31	29 October 2019	47	16 November 2020
16	22 October 2018	32	22 November 2019	48	10 December 2020

.

2.4. DEM Data

The Advanced Land Observation Satellite, Phased Array Type L-band Synthetic Aperture Radar (ALOS-PALSAR) DEM is selected as our DEM data. ALOS is a land observation satellite launched by Japan in 2006, which carries an L-band synthetic aperture radar. It has a high spatial resolution of 12.5 m, which facilitates the extraction of terrain information.

3. Methodology

The following Figure 2 is the flowchart of our methodology in this study:

3.1. Vegetation Information Extraction

3.1.1. GF-1 Image Pre-Processing

Pre-processing of GF-1 images including radiometric calibration, atmospheric correction through the fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) method [14] and geometric rectification is helpful to obtain true ground reflectance and eliminate geometric distortions. Then, in order to strength the vegetation information extraction, the nearest neighbor diffusion method with high spectral fidelity is used to fuse the registered PAN and MS images to increase the spatial resolution to 2 m while keeping the original spectral information of MS unchanged as much as possible. Then, GF-1 images are aligned to eliminate the image mismatching between different phases.

3.1.2. Study Area Zoning

In order to ensure the relative consistency of natural conditions within the same geographic unit, we have zoned the Temi landslide. Integrating GF-1 images, DEM data, slope map and Google Earth, and referring to the developmental characteristics and geographical conditions of Temi landslide, we roughly divide the Temi landslide into three zones, including the upper sliding source zone, the middle sliding accumulation zone and the lower main accumulation zone (see Figure 3) [13]. The upper sliding source zone has a clear back edge boundary, where the source of the gully develops, with an elevation range of 2990–3190 m and a slope degree range of 40–50°. Moreover, the weathering is strong here, and there is a winding road at 200 m-distance from the back edge. The middle sliding accumulation zone has an elevation range of 2850–2980 m and a slope of 20–30°, within which the vegetation is well developed. The lower main accumulation zone ranges in elevation from 2550 m to 2850 m, with a steep front edge from 50° to 60°, and the rest of the slope is mostly 10–20°. Thus, the obvious platform area is formed there, the gully is well developed, and the catchment area is also formed.

3.1.3. Calculation of FVC

As the cover of the landslide surface, vegetation reflects the internal activities of the landslide and the influence of various factors. Fractional vegetation cover (FVC) is an indicator that directly characterizes the growth status of vegetation. The normalized difference vegetation index (NDVI) is often used to estimate FVC with the dimidiate pixel model, which can weaken the interference of atmospheric, soil and vegetation type, and further greatly improves the calculation accuracy of FVC [15]. The formula for calculating NDVI is as follows:

$$\mathrm{NDVI}=\frac{{\mathsf{\rho }}_{\mathrm{NIR}}-{\mathsf{\rho }}_{\mathrm{R}}}{{\mathsf{\rho }}_{\mathrm{NIR}}{+\mathsf{\rho }}_{\mathrm{R}}} ,$$

(1)

where ${\mathsf{\rho }}_{\mathrm{NIR}}$ represents the reflectance of the near infrared band, and ${\mathsf{\rho }}_{\mathrm{R}}$ represents the reflectance of the red band.

The dimidiate pixel model divides the information acquired by the sensor into vegetation and soil components. The sum of the vegetation area ${\mathrm{A}}_{\mathrm{veg}}$ and the soil area ${\mathrm{A}}_{\mathrm{soil}}$ is 1 over the cell area composition, which is shown in Equation (2):

$${\mathrm{A}}_{\mathrm{veg}}{+\mathrm{A}}_{\mathrm{soil}}=1,$$

(2)

The NDVI values of pure vegetation (covered only by the vegetation) and pure soil (covered only by the soil) are represented by ${\mathrm{NDVI}}_{\mathrm{veg}}$ and${ \mathrm{NDVI}}_{\mathrm{soil}}$, respectively, then the NDVI value of the mixed pixel is shown in Equation (3).

$${\mathrm{NDVI}=\mathrm{A}}_{\mathrm{veg}}·{\mathrm{NDVI}}_{\mathrm{veg}}{+\mathrm{A}}_{\mathrm{soil}}·{\mathrm{NDVI}}_{\mathrm{soil}} ,$$

(3)

Combining the above, we can know:

$${\mathrm{A}}_{\mathrm{veg}}=\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\mathrm{soil}}}{{\mathrm{NDVI}}_{\mathrm{veg}}-{\mathrm{NDVI}}_{\mathrm{soil}}} ,$$

(4)

i.e.,

$$\mathrm{FVC}=\{\begin{array}{c}0\\ \frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\mathrm{soil}}}{{\mathrm{NDVI}}_{\mathrm{veg}}-{\mathrm{NDVI}}_{\mathrm{soil}}}\hfill \\ 1\end{array}\begin{array}{c},\mathrm{NDVI}\le {\mathrm{NDVI}}_{\mathrm{soil}}\hfill \\ { ,\mathrm{NDVI}}_{\mathrm{soil}}\le \mathrm{NDVI}\le {\mathrm{NDVI}}_{\mathrm{veg}}\hfill \\ ,\mathrm{NDVI}\ge {\mathrm{NDVI}}_{\mathrm{veg}}\hfill \end{array}$$

(5)

Therefore, determining ${\mathrm{NDVI}}_{\mathrm{veg}}$ and ${\mathrm{NDVI}}_{\mathrm{soil}}$ is the key to calculating FVC. Generally, ${\mathrm{NDVI}}_{\mathrm{veg}}$ and ${\mathrm{NDVI}}_{\mathrm{soil}}$ are determined through on-the-spot investigation, but it is difficult to reach the Temi landslide due to the special geographical environment. Many studies have shown that it is reliable to determine ${\mathrm{NDVI}}_{\mathrm{veg}}$ and ${\mathrm{NDVI}}_{\mathrm{soil}}$ by setting confidence intervals through the NDVI statistical table [16,17,18]. We try different confidence intervals to calculate FVC. The confidence intervals that best matches with the features in GF-1 images are selected.

We try different confidence intervals to calculate FVC of different phases and different sub-areas. According to the matching degree between the maps of FVC and the actual features, the best matching confidence intervals are selected. Due to the differences in the distribution of ground objects in different sub-areas, the selected confidence intervals are different. There are also some deviations due to the subjectivity of visual interpretation in images of different phases in the same sub-areas. In order to make the analysis convenient and accurate, the confidence intervals in this case should be consistent as much as possible. Ultimately, 1–99% confidence intervals are chosen for the upper sliding source zone and the lower main accumulation zone; 0.5–99.5% confidence intervals are chosen for the middle sliding accumulation zone. The FVC values in three sub-areas are calculated by Equation (5).

3.2. Surface Deformation Information Extraction

The commonly used time-series InSAR methods are mainly PS-InSAR and SBAS-InSAR. The PS-InSAR method obtains surface deformation information by identifying the deformation rate and elevation difference for the persistent scatterer (PS). However, PS points are usually man-made feature targets that are suitable in urban environments [8]. The SBAS-InSAR method uses the distributed scatterers (DSs), which are highly coherent over a short time period and are commonly found in nature including sand, grassland, bare soil, etc. [19]. Hence, the SBAS-InSAR method is generally used in natural scenes. Moreover, the SBAS-InSAR method is still effective in areas covered by vegetation [20]. Therefore, the SBAS-InSAR method is selected to obtain the time-series deformation information of the Temi landslide.

The SBAS-InSAR method forms multiple small baseline subsets of SAR images based on temporal and spatial baselines to improve the interferometric coherence. Each image is interferometrically measured to generate the interferogram, and then the phase unwrapping and singular value decomposition (SVD) are performed. Finally, the surface deformation information and terrain elevation errors are obtained using the temporal filtering method. Specifically, The SBAS-InSAR processing performed in SARSCAPE software includes the generation of connection graph, the interferometric process, the orbit refinement and re-flattening, the twice-step inversions, and the geocoding. The initial step is the generation of a connection graph, which is achieved by setting up a temporal baseline and a spatial baseline to obtain interferometric pairs. Next, the interferometric process is performed on all paired interferometric image pairs. We use the Goldstein method for filtering and the minimum cost flow (MCF) method for phase unwrapping. For the orbit refinement and re-flattening, the ground control points (GCPs) are used to remove noisy phases, such as track errors. Afterwards, the deformation rate and the residual topographic phase are estimated in the first-step inversion. In the second-step inversion, we remove the atmospheric phase to obtain accurate time-series surface displacement results. Finally, the SAR coordinate system is converted to the geographic coordinate system by geocoding, which is consistent with the GF images we use.

3.3. Landslide Susceptibility

3.3.1. Landslide Susceptibility Factors

Landslide susceptibility is the magnitude of probability or vulnerability of landslides occurring in a given area. It is generally considered as a spatial concept rather than a temporal concept. However, landslides are, in fact, constantly changing with time. It is difficult to accurately predict the landslide occurrence by considering the landslide susceptibility as a static concept. Hence, in this paper, we extend the landslide susceptibility to the temporal domain through the spatial susceptibility in time-series.

The landslide susceptibility assessment directly depends on different factors [21]. It is widely recognized that the numerous influence factors of the landslide susceptibility include the terrain and geomorphological conditions, the hydrological factors, the surface cover factors, the surface deformation information, and the human activities. Terrain and geomorphological factors include the elevation, the slope, the relative relief and the slope direction. Since the slope direction is unique in the single landslide, it is not used for modelling. The elevation, slope, and relative relief are selected for modelling in terrain and geomorphological factors. Hydrological factors cover the distance from rivers, and the variability of groundwater. The topographic wetness index (TWI) is generally used to reflect the hydrological situation of landslide. We also use TWI for the hydrological factor. NDVI is adopted for the surface cover factor in modelling. Surface displacement of a landslide accurately describes the deformation characteristics of landslides and have a direct impact on landslide susceptibility, which is used for the surface deformation information in modelling. In addition, human factors include the distribution of residential and agricultural land, and engineering projects such as the construction of roads and bridges. Human influences are not included in the modelling as they are not easily quantified and evaluated.

Elevation is an essential component affecting the landslide [22]. The climatic conditions and rock conditions of the landslide differ at different elevation values. In general, as the elevation increases, the environment of the landslide becomes more vulnerable and the susceptibility to the landslide becomes greater. The elevations of the Temi landslide are distributed between 2486 m and 3154 m, as shown in Figure 4a.

Slope is identified to act a vital role for triggering landslide phenomena [23]. The slope is the angle between the tangent plane and the horizontal plane at any point on the slope. The landslide slides under the influence of gravity, and it is more likely to occur as the slope increases. The susceptibility of landslide is considered to be positively correlated with the slope. The slope of the Temi landslide ranges from 0.810231° to 62.6101°, as shown in Figure 4b.

Relative relief is the difference between the highest elevation and the lowest elevation within an area [24]. It reflects the differences in terrain. The greater the relative relief, the greater the susceptibility to the landslide. The relative relief of the Temi landslide is between 1 m and 60 m, as shown in Figure 4c.

TWI is related to the hydrological flow accumulation at the given terrain [25]. It quantifies the role of topography on hydrology [26], which is efficiently used to find out the soil moisture condition and other related phenomena [27], as shown in Equation (6) [28]. The stability of a landslide is influenced by the erosive effect of hydrology on the slope, which acts as a lubricant for the landslide. TWI exerts a positive influence on the landslide susceptibility. The TWI for the Temi landslide ranges from 1.78798 to 15.0784, as shown in Figure 4d.

$$\mathrm{TWI}=\ln \frac{\mathrm{CA}}{\tan \mathrm{slope}} ,$$

(6)

in which $\mathrm{CA}$ stands for the catchment area and $\mathrm{slope}$ represents the slope gradient.

Several studies have concluded that the vegetation change can indicate landslide creep. Vegetation is considered to be a protective layer for landslides, which means that phenomenon of landslides decreases with increasing vegetation cover [24,29]. NDVI indicates the vegetation condition. The higher the NDVI value, the lower the probability of landslides occurring. Therefore, it is believed that NDVI plays a negative role in landslide susceptibility, as shown in Figure 4e for NDVI (2018), Figure 4f for NDVI (2019), Figure 4g for NDVI (2020) and Figure 4h for NDVI (2021).

The surface displacement of landslide directly affects the landslide susceptibility [30]. Several studies have shown that taking the surface displacement into consideration can effectively evaluate the landslide susceptibility. Furthermore, the surface displacement plays a positive role in the landslide susceptibility.

As the factors above are derived from different data sources, there are differences in coordinate systems, numbers of rows and columns, and spatial resolution, which affects the landslide modelling. Therefore, it is necessary to unify them from multiple sources. In order to make the modelling independent of the magnitude and scale of different factors and ensure that each pixel has a value, the factors need to be normalized separately. The factors that have positive values and contribute to the landslide susceptibility are normalized directly, including elevation, slope, relative relief and TWI. However, the surface displacement values are negative. The larger the absolute value of the surface displacement, the greater the impact on the landslide. Moreover, the time series is considered. So, the absolute value of the surface displacement needs to be normalized relatively. Firstly, we compute the absolute values of all the surface displacements. Secondly, we select the maximum and minimum values of all the absolute values. Then, we normalize all the absolute surface displacement values using the maximum and minimum values. NDVI has the opposite effect on the landslide susceptibility compared to other factors. With the increase in NDVI, the landslide susceptibility decreases. The NDVI in all the time series is relatively inverse normalized, i.e., the larger the NDVI value, the smaller the effect on the landslide susceptibility.

3.3.2. Landslide Susceptibility Modeling

The factors for the landslide susceptibility model influence the landslide and also influence each other. Through this relationship, the analytic hierarchy process (AHP) method is used to determine the weights of different factors. AHP is a decision-making method proposed by Saaty in 1990 [31]. It is still widely used at present, and it simplifies the complex decision-making process into the weights assessing process of factors. According to Kaliraj [32] and Basu [24], the relationship between various factors is first considered, where the 0.5–1 weighting scoring method is used to determine the value of the vertical factors on the horizontal ones. If an indicator has a direct impact on another, the score of 1 is assigned. For an indirect impact, the assigned score of the indicator is 0.5. For example, the slope directly affects the vegetation (NDVI), which is assigned a score of 1, while the vegetation (NDVI) affects the slope indirectly, which is assigned a score of 0.5. The sum of the scores constitutes the scale weight. Thus, the AHP pairwise matrix table is shown in

Table 3. Assignment of score for AHP pairwise matrix table.

Factors	Elevation	Slope	Relative Relief	TWI	NDVI	Surface Displacement	Scale Weight
Elevation	1	0.5	0.5	0.5	1	0.5	4
Slope	0.5	1	1	0.5	1	1	5
Relative Relief	0.5	1	1	1	0.5	0.5	4.5
TWI	0.5	0.5	0.5	1	0.5	0.5	3.5
NDVI	0.5	0.5	0.5	0.5	1	1	4
Surface Displacement	0.5	0.5	0.5	0.5	1	1	4

.

Another aspect of the AHP method is the calculation of the consistency index (CI) and consistency ratio (CR). CI is calculated as follows:

$${\mathrm{CI}=(\mathsf{\lambda }}_{\max }-\mathrm{n})/(\mathrm{n}-1) ,$$

(7)

in which ${\mathsf{\lambda }}_{\max }$ stands for the maximum eigenvalue, $\mathrm{n}$ represents the number of factors in the AHP pairwise matrix.

The larger the CI value, the greater the deviation of the pairwise matrix from the consistency. The smaller the CI value, the better the consistency of the pairwise matrix. However, the evaluation of the pairwise matrix by CI alone may not meet the actual situation. CR is another widely used index. When CR is less than 0.1, the matrix is considered to have an excellent consistency. CR is calculated as follows [33]:

$$\mathrm{CR}=\mathrm{CI}/\mathrm{RI} ,$$

(8)

in which RI is the random index representing the consistency of a randomly generated pairwise matrix [34]. RI is usually obtained by the lookup table [35], as shown in

Table 4. Values of the random index (RI).

N	1	2	3	4	5	6	7	8	9	10
RI	0	0	0.58	0.90	1.12	1.24	1.32	1.41	1.45	1.49

.

According to the pairwise matrix in Table 3, the maximum eigenvalue is 9.3747 and n is 6, so CI equals 0.0468. Based on Table 4, RI is 1.24 when n is 6. As a result, CR is 0.04258 which is less than 0.1. Therefore, the matrix has satisfactory consistency and is feasible. The normalized weights are obtained by normalizing the scale weights. The normalized weights of all factors are shown in

Table 5. Normalized weights of factors.

Factors	Elevation	Slope	Relative Relief	TWI	NDVI	Surface Displacement
Weights	0.16	0.2	0.18	0.14	0.16	0.16

. From Table 5, the slope has the greatest contribution to the landslide among all the factors, followed by the relative relief.

The landslide susceptibility index (LSI) is used to represent the landslide susceptibility, which is calculated by Equation (9).

$$\mathrm{LSI}=\sum { \mathrm{w}}_{\mathrm{i}}\times \mathrm{factor},$$

(9)

in which${ \mathrm{w}}_{\mathrm{I}}$ is the weight of each factor, and $\mathrm{factor}$ represents the factors of LSI.

The final model is as follows:

$$\mathrm{LSI}=0.16\times \mathrm{Elevation}+0.2\times \mathrm{Slope}+0.18\times \mathrm{Relative} \mathrm{Relief}+0.14\times \mathrm{TWI}+0.16\times \mathrm{NDVI}+0.16\times \mathrm{Surface} \mathrm{Displacement}$$

(10)

4. Results and Discussion

Maps of FVC from 2018 to 2021 are shown in Figure 5. Based on the results, it can be seen that FVC in the upper sliding source zone is the lowest, followed by the middle sliding accumulation zone, and the lower main accumulation zone has the highest FVC, except for the circular arc front edge area, which has almost no vegetation cover.

We extract surface displacement images from 2018 to 2021. Then, the images in which time phases are closest to the GF-1 images are shown in Figure 6. The maps of surface deformation rate are dropped onto Google Earth due to the lack of obvious feature information in the SAR images, as shown in Figure 7. The surface deformation rate of the landslide subsidence is negative, so only negative values are displayed. Most of the coherent points are distributed in the upper and middle part of the landslide. According to the classification criteria of the landslide surface deformation rate [36], the deformation rate of the Temi landslide is between 0 and 10 mm/a, which is less than 16 mm/a. Hence, the Temi landslide is in an extremely slow creep stage.

The LSI maps for 2018 to 2021 are shown in Figure 8. Based on the maps of LSI over these years, it is clear that LSI is the greatest in the upper sliding source zone, followed by the middle sliding accumulation zone, and the lower main accumulation zone is lowest, except for the circular arc front edge area. The upper sliding source zone is at a higher elevation and slope, with poorer vegetation cover than the other two areas, and has a higher probability of sliding. The circular arc front edge area in the lower main accumulation area created a clear demarcation due to the previous slides. For the front edge area, the slope angle is high with no vegetation and low elevation. The LSI is high, which shows the dominant role of slope in the assessment of landslide susceptibility.

4.1. Analysis of Spatial Variability of Landslide

The variation patterns of FVC, surface displacement and LSI over the entire landslide body are first considered. As the SBAS-InSAR method acquires the surface deformation information along the line of sight, not along the direction of landslide movement, the longitudinal profile of the landslide is made to obtain the surface displacement information from the top to the bottom of the landslide. Since there is no vegetation cover on the circular arc front edge of the landslide, this area is not considered. Therefore, using the latest images, a profile line is made along the longitudinal profile from the back edge to the circular arc front edge of the landslide, as shown in Figure 9a. Then, the variation in the three aspects along the profile is analyzed. A, B and C are the typical positions along the longitudinal profile line to illustrate the change regularity. When FVC decreases, the vegetation condition deteriorates. As the surface displacement decreases, its value increases but its absolute value decreases. A decline in LSI means that the landslide becomes less susceptible and more stable. It is evident that FVC shows the same trend as the surface displacement, while being the opposite of LSI. As shown in Figure 9b–d, the FVC situation is at a small trough at position A, while the surface displacement is also at a low point and LSI is at a small peak. The pattern at position C is the same as the position A. At position B, FVC is at a small height, and the surface displacement is at a relative high place, while LSI is at a small low level.

4.2. Analysis of Temporal Variability of Landslide in Sub-Areas

Three sub-areas are selected in the upper sliding source zone, the middle sliding accumulation zone and the lower main accumulation zone, named Area1, Area2 and Area3, respectively from top to bottom, as shown in Figure 10.

There are five classes of FVC, level 1 (FVC < 10%), level 2 (10% ≤ FVC < 30%), level 3 (30% ≤ FVC < 50%), level 4 (50% ≤ FVC < 70%), level 5 (FVC ≥ 70%) [37]. Among the three sub-areas, Area1 uses the threshold value of 50% to calculate FVC due to its poor vegetation cover. Then, the percentage of the number of pixels with FVC greater than 50% is calculated, i.e., the proportion of the number with FVC of level 4 and level 5. Area2 and Area3 use the threshold value of 70%, and the percentage of the number of pixels with FVC of 70% is calculated, i.e., the proportion of the pixel number with FVC of level 5. The surface displacement and LSI are represented by their mean values to represent their average levels. In all three sub-areas, FVC and surface displacement show the trend of “increasing—increasing—decreasing” year by year, while LSI shows a “decreasing—decreasing—increasing” curve. When FVC rises, the pixel numbers with high FVC increase. At the same time, the surface displacement also shows an upward trend, which means the surface displacement becomes smaller because its value is negative. On the other hand, LSI shows a downward trend when the vegetation condition becomes better and the surface displacement becomes smaller. All these three aspects show that the probability of sliding becomes less. Accordingly, there is a relationship between landslide vegetation, surface displacement and landslide susceptibility, which indicates that the smaller the surface displacement, the better the vegetation and the less the landslide susceptibility. Meanwhile, the larger the surface displacement, the worse the vegetation and the greater the landslide susceptibility.

Figure 11, Figure 12 and Figure 13 show the change analysis in Area1, Area2 and Area3, respectively. The changes in the three sub-areas are similar, so Area2 can be taken as an example for concrete analysis. As shown in Figure 12, the blue area in the FVC map becomes smaller from 2018 to 2019, which indicates a better vegetation condition, while the purple area in the surface displacement map becomes larger, which suggests a smaller surface displacement, and the red area in the LSI map decreases, which indicates lower landslide susceptibility. From 2019 to 2020, the red area in the FVC map expands and vegetation conditions improve, while the purple area in the surface displacement map enlarges and surface displacement declines. Meanwhile, the purple area in the LSI map expands and the landslide becomes less susceptible. In the period of 2020–2021, the red region in FVC shrinks, indicating the deteriorated vegetation status; the purple region in surface displacement decreases, indicating the improved surface displacement; and the red region in LSI increases indicating the increased landslide susceptibility.

4.3. Analysis of Temporal Variability of Landslide at Center Points

In each of the three sub-areas, three center points are selected to illustrate the change mode of vegetation, surface deformation rate and landslide susceptibility from 2018 to 2021 (see Figure 14). The surface deformation rate is negative in almost all the 62 time-phased sentinel-1A images, suggesting that it is the sinking movement of the landslide. The Sentinel-1A images closest to the GF-1 images are chosen for the change analysis of the surface deformation rate from 2018 to 2021, as shown in Figure 12. In the yearly variation in the three aspects from 2018 to 2021, it is obvious that FVC and surface deformation rate both follow the “increase—increase—decrease” mode, while LSI shows the contrary “decrease—decrease—increase” mode. The change modes at center points are consistent with the three sub-areas.

5. Conclusions

Considering the limitations of current landslide monitoring tools, this paper proposes a comprehensive spatio-temporal monitoring method for landslides by using the combination of landslide vegetation, landslide surface deformation and landslide susceptibility based on the time-series multi-source remote sensing data. The extracted vegetation information from multitemporal optical images indirectly reflects the landslide creep. The calculated surface deformation information from multitemporal interferometry SAR images directly reflects the landslide creep. The designed susceptibility model considering multiple factors overall reflects the landslide status. The landslide susceptibility is extended to the spatio-temporal domain from only the spatial domain. The Temi landslide with high hazardous and highly concealed characteristics is used as our study area. The results of this study show that:

(1) The vegetation covering the landslide surface is changing as the landslide develops. The FVC of the Temi landslide shows an increasing trend from top to bottom except the circular arc front edge in the lower part of the lower main accumulation zone without vegetation.

(2) The surface deformation rate owing to movement of the Temi landslide is between −10 and 0 mm/a except for the uplift caused by rivers and other reasons, which shows that the landslide is currently in an extremely slow creep period.

(3) Landslide susceptibility is not purely a spatial concept. As the various factors affecting landslide are dynamic, landslide susceptibility can be extended to the temporal dimension to recognize the variability of the landslide with time. Among the factors for modelling the landslide susceptibility, the slope has the most considerable contribution on the landslide susceptibility. The susceptibility of the Temi landslide derived by the AHP method suggests that the landslide susceptibility of the upper part is the highest, followed by the middle zone and the lower area is the lowest.

(4) The variation in vegetation condition, surface deformation and landslide susceptibility of the Temi landslide shows the consistency in space and time. Spatially, there is relatively small surface deformation with high vegetation cover along the direction of the landslide subsidence, where there is less susceptibility to landslide, meaning that sliding is less likely to occur. In temporal terms, the vegetation cover and the surface displacement from 2018 to 2021 follow the same trend with the landslide susceptibility. It is quite apparent that when the surface deformation increases, the vegetation situation is worse, while at the same time the land becomes more vulnerable to sliding.

Based on the consequences of the Temi landslide, we can conclude that considering multiple factors affecting landslides with multi-source data to analyze the spatial and temporal correlation of landslide susceptibility, surface deformation and vegetation cover, allows for fine monitoring of the landslide variability. This idea provides a feasible perspective for landslide monitoring and aids natural disaster prevention and mitigation. Nevertheless, as only detailed analyses of one single landslide have been carried out, it may not be possible to draw general conclusions. The extended experiments for other landslides should perform in the future.