Potential landslides identification based on temporal and spatial filtering of SBAS-InSAR results

Abstract

Interferometric Synthetic Aperture Radar (InSAR) is an important method for acquiring surface deformation. Considering the difficulty of the identification work, the identification of landslides needs to be combined with the context of the pregnant disaster and the precipitating conditions. To identify potential landslides, we applied spatial and temporal filtering to the InSAR results, which consists of rainfall and landslide susceptibility mapping. In this paper, taking the Badong Ecological Barrier Zone of the Three Gorges reservoir area as the study area, the deformation aggregation areas in the study area were obtained by applying Small Baseline Subset InSAR (SBAS-InSAR) technology and spatial statistical analysis. We screened deformation aggregation areas by combining the susceptibility map and the correlation analysis of rainfall and deformation. Field verification and investigation were conducted on the suspected deformation areas, and 11 landslides were found to have signs of deformation, two of them are newly discovered landslides. In addition, we selected one of the landslides, the Songjiawuchang landslide, and compared the InSAR results with the GPS accumulated displacement to verify the reliability of the results. This research demonstrated the feasibility of combining InSAR results with spatial susceptibility maps and monthly rainfall factors for landslides identification methods.

Key policy highlights

• Applied Spatial filtering and temporal filtering on the InSAR results.

• Evaluated in conjunction with geological hazards themselves.

• Landslide identification accuracy significantly improved.

Keywords:

• Landslides identification
• SBAS-InSAR
• landslide susceptibility mapping
• Pearson correlation coefficient
• hot spot analysis
• Three Gorges Reservoir Area

1. Introduction

Identification of potential landslides is essential for disaster mitigation and prevention. InSAR is a new spatial surface measurement technology, that is unaffected by natural weather effects and can observe the surface at all hours of day and night, allowing it to be used for the identification of potential landslides. InSAR technology can collect high-resolution and high-precision information on small surface deformations over a large area of the ground. Applying InSAR to identify potential hidden hazards in key areas, we can realize the investigation of key areas from the perspective of time and space, thus providing enough reference basis for the prevention and control of regional geological hazards (Akbarimehr et al. 2013; Medhat et al. 2022; Milillo et al. 2022; van Natijne et al. 2022). Many researchers have already applied InSAR technology to identify regional landslides, update landslide inventory maps, determine potential landslide extent, and identify potential landslides precursor features (Intrieri et al. 2018; Devara et al. 2021; Guo et al. 2021; Su et al. 2021).

Some key challenges are still in the use of InSAR for landslides identification, such as the low-density scatterers in mountainous areas, error sources masking the displacement signals, deformation in line-of-sight, define of warning thresholds, and deformation patterns (Schlögl et al. 2022). The SAR image and the range of applications for InSAR technology make it clear that it is not sufficient to identify a geological hazard solely on the basis of InSAR results. It is required to increase the identification accuracy by doing a fundamental analysis of the study region, in addition to the conditions of the study area that foster disasters and the landslide’s contributing causes (Ge et al. 2019).

To investigate a novel idea of universal landslide hazard, it is necessary to create a comprehensive hazard identification system using multi-source data and multi-dimensional monitoring. In this system, we can start from two aspects, one is the accuracy of InSAR for landslide identification, and the other is the development environment of landslide. Most scholars focus on one perspective, and few carry out landslide identification from both perspectives. In terms of the analysis of the pregnant environment, landslide susceptibility is a good choice (Singh et al. 2005; Ciampalini et al. 2015; Rosi et al. 2018; Novellino et al. 2021; Hussain et al. 2022). The combination of InSAR and landslide susceptibility is a regional disaster research method. Landslide susceptibility has been employed as one of the key ways to examine regional landslide differentiation and has important recommendations for the regional control of landslides (Wang et al. 2020). Landslide susceptibility can be used to determine the likelihood of a landslide occurring in the research region and can be used to offer a basic guarantee for landslide detection.

The concept of integrating spatial and temporal filtering was used in this paper to identify potential landslides in the Badong ecological barrier zone. Spatial filtering is the process of implementing landslide susceptibility to filter the research area’s low-susceptibility landslide zones. The use of Pearson correlation coefficients to filter areas where deformation is little correlation with precipitation during the study period is known as temporal filtering. In the identification process, four technical approaches namely, SBAS-InSAR technology, GIS spatial analysis, Pearson correlation coefficient, and landslide susceptibility analysis were mostly used. The identification technique that combines spatial and temporal filtering complements the shortage of single information to identify landslides and improves the accuracy of detecting active landslides. Additionally, based on monthly rainfall and Global Positioning System (GPS) monitoring data, the time-series deformation of the Songjiawuchang landslide was examined, and the findings not only support the accuracy of the identification results but also aid in landslide prevention and warning.

2. Study area and data

2.1. Study area

The Badong Ecological Barrier Zone is in the northern part of Badong County, Hubei Province, China, and covers the Yangtze River’s main channel and tributaries, with longitudes ranging from 110°5′E to 110°29′E and latitudes ranging from 31°16′N to 30°57′N, covering an area of around 309 km². The barrier zone has a well-developed water infrastructure and is a highly populated area as well as an important ecological zone. The research area is part of the Chuan-E folded mountainous area, with high topography from the north and south and low terrain in the centre. The geological conditions in the area are complex, with abundant precipitation and strong human engineering activities, resulting in frequent geological disasters (Fourniadis et al. 2007; Juang 2021). The lithology of the exposed strata is mudstone, sandstone interspersed with marl and tuff of the Middle Triassic Badong Formation (T₂b). The slope structure type is soft and hard rocky with a compliant slope. Due to the deep cutting of the Yangtze River and the strong human economic and engineering activities, the landslide is extremely developed, like most large and medium-sized landslides. In total, there are 242 landslides in the area (Figure 1). Statistically, 174 landslides’ deformation was mainly caused by rainfall. Considering that the InSAR method can be useful for extremely slow and very slow landslide deformation monitoring during the primary and secondary stages (Xiong et al. 2020), we aimed to identify landslides that produce slow deformation under the influence of rainfall.

Figure 1. (a) Location of the Three Gorges reservoir area and Badong County. The star indicates the location indication of the study area in the Three Gorges reservoir area. (b) Badong County and ALOS-2 PALSAR-2 images coverage. (c) Shape location of the study area, and the distribution of landslides and water systems in the study area.

2.2. SAR data

This study used ALOS-2 PALSAR-2 images in the L-band (1.2 GHz band) with Ultra-Fine mode, HH + HV polarization, greater than 3 m resolution, and a 55 km × 70 km width. PALSAR-2 images have the advantage of high resolution and strong penetration, which is useful for detecting small deformations in the disaster body and has better application benefits in monitoring areas with dense foliage cover (Suzuki et al. 2011). The SAR data is available from 2020/4/11 to 2021/8/14 and a total of 14 frames of SAR images were utilized (Table 1), with a centre incident angle of 42.9°.

Table 1. ALOS-2 PALSAR-2 images schedule.

Data number	Date	Data number	Date
1	2020/4/11	8	2020/11/7
2	2020/5/23	9	2021/1/16
3	2020/6/20	10	2021/3/13
4	2020/7/4	11	2021/4/10
5	2020/8/1	12	2021/6/5
6	2020/9/12	13	2021/7/3
7	2020/10/10	14	2021/8/14

2.3. Landslide causative factors

The Geological Environmental Center of Hubei Province provided the research area’s landslide inventory map, which included 182 sites, and this landslide inventory is official data that has been manually fielded and verified. The 182 landslides were then digitized and rasterized in the Environmental Systems Research Institute (ESRI)’s ArcGIS software (version 10.5.0) at a spatial resolution of 12.5 m, which was the same as the Digital Elevation Model (DEM) data used in this study.

The spatial location of landslide occurrence is predicted via landslide susceptibility mapping. For the susceptibility analysis of landslides, six factors that are more important for landslides were selected based on the information from previous studies, including elevation, slope, Topographic Wetness Index (TWI), lithology groups, land use/land cover, and Normalized Difference Vegetation Index (NDVI) (Chen et al. 2018; Mallick et al. 2018; Rong et al. 2020). The study’s DEM data sources were ALOS PALSAR 12.5 m data from Alaska Satellite Facility-Distributed Active Archive Centers (ASF DAAC, https://search.asf.alaska.edu/#/) for extracting elevation, slope, and TWI; a 1:200,000 geological map for extracting lithology groups factor; 2017 landuse change survey data for extracting land use/land cover data; and NDVI was produced from two Landsat-8 Operational Land Imager (OLI) images (path/row 125/38, 125/39) obtained on September 18, 2020. To facilitate machine learning model training, the six variables were categorized according to the landslide frequencies (Niu et al. 2014; Chen et al. 2016), displayed hierarchically (Figure 2), and assigned to the training datasets (Table 2).

Figure 2. Landslide causative factors maps: (a) Elevation; (b) slope; (c) NDVI; (d) lithology groups; (e) TWI; (f) land use/land cover. The lithology groups are described as follows: I: Hard medium to thick-bedded massive dolomite metamorphic conglomerate rock group; II-1: Hard medium to thickly laminated quartz sandstone quartzite fine sandstone rock group; II-3: Hard harder medium to thick-bedded sandstone mud siltstone interbedded with shale coal seams and mudstone shale interbedded rock group; II-4: Harder to softer thinly to moderately thickly bedded shale sandstone mudstone rock group; II-5: Soft thin to medium-thick bedded muddy siltstone shale rock group; III-1: Hard medium to thick-bedded strongly karsted carbonate rock group; III-2: Harder harder medium to thickly bedded strongly to moderately karsted carbonate rock group; III-3: Harder thin to medium-thick bedded weakly karsted carbonate rock group.

Table 2. Values for each category of causative factors.

Factors	Class	Landslides	Frequency of landslide	Classification value
Elevation	0–200 m	81	37%	1
	200–400 m	105	48%	2
	400–600 m	19	9%	3
	>600 m	16	7%	4
Slope	0–10°	6	3%	1
	10–20°	56	25%	2
	20–30°	99	45%	3
	30–40°	50	23%	4
	>40°	10	5%	5
NDVI	<0.1	30	14%	1
	0.1–0.3	34	15%	2
	0.3–0.5	40	18%	3
	0.5–0.7	74	33%	4
	>0.7	43	19%	5
Lithology groups	I	0	0%	1
	II-1	0	0%	1
	II-3	4	2%	2
	II-4	3	1%	3
	II-5	113	51%	5
	III-1	33	15%	3
	III-2	6	3%	2
	III-3	62	28%	4
TWI	<4.5	8	4%	1
	4.5–5	44	20%	2
	5–5.5	95	43%	3
	5.5–6	50	23%	4
	6–6.5	19	9%	5
	>6.5	5	2%	6
Land use/land cover	Arable land	48	22%	1
	Urban	20	9%	2
	Rural	13	6%	3
	Waters	27	12%	4
	Woodland	113	51%	5
	Grassland	0	0%	6
	Others	0	0%	7

As a precipitating factor, rainfall plays a crucial influence in the deformation of landslides (Bayer et al. 2018). The Hubei Provincial Hydrological and Water Resources Center supplied monthly precipitation data from 23 hydrological stations. The monthly precipitation distribution maps from April 2020 to July 2021 were created by interpolating the hydrological stations using the Spline tool in ArcGIS 10.5.0 software to assist the spatial calculation with the monthly cumulative deformation maps of InSAR. Interpolation predicts values for cells in a raster from a limited number of sample data points. The Spline tool uses an interpolation method that estimates values using a mathematical function that minimizes overall surface curvature, resulting in a smooth surface that passes exactly through the input points (Hadi and Tombul 2018).

3. Methodology

3.1. General method and technical route

In this research, we highlight SBAS-InSAR, rainfall and deformation correlation computation, obtaining landslide susceptibility, and identification as the four crucial steps in the process of identifying landslides. In general, the Pearson correlation coefficient can be used to characterize the temporal correlation between the deformation location of rainfall-type landslides and the rainfall (Liu et al. 2022). When we performed correlation analysis for each point of the study, the Pearson correlation coefficient result of rainfall and deformation can reflect the areas that vary over time. The spatial distribution probability of landslides in the studied area can be determined using landslide susceptibility zoning. The article filters the InSAR results in time and space using the Pearson correlation coefficient and the landslide susceptibility distribution map to obtain the final identification results.

To acquire the deformation aggregation areas in the research area, the spatial statistical analysis was first done utilizing the deformation rates obtained by SBAS-InSAR. Then, Pearson correlation coefficient analysis was done using monthly rainfall data and monthly cumulative deformation variables to create the distribution map of deformation correlation coefficients in the area. The GBDT model was used to calculate the influencing factors and landslide data set at the spatial distribution level to generate the landslide susceptibility distribution map of the barrier zone. Two steps make up the filtering procedure for the deformation aggregate region, one for the temporal level and the other for the spatial level. The first screening of deformation aggregation areas was performed at the time change level to reserve deformation areas with rainfall correlation. The second screening of deformation aggregation areas was performed at the spatial change level to reserve deformation areas with medium or high susceptibility. Finally, an outdoor geological survey of the final reserved deformation aggregation areas was conducted to obtain the landslides identification results of the study area (Figure 3).

Figure 3. Flow chart for identifying potential landslides in Badong based on spatial and temporal filtering of SBAS-InSAR results.

3.2. SBAS-InSAR

The SBAS-InSAR technique is used to extract surface deformations. It was formally proposed by Berardino et al. (2002) to solve the decoherence problem caused by too-long spatial baselines in the D-InSAR technique. SBAS has become one of the most representative techniques for extracting time series surface deformation (Tizzani et al. 2007; Dong et al. 2014; Zhao et al. 2016; Novellino et al. 2017).

From the principle standpoint, the adoption of the least-squares principle and singular value decomposition is the most fundamental feature of the small baseline set technique (Berardino et al. 2002). The small baseline set technique screens the SAR images with short time and spatial baselines according to the threshold value solves the time-series surface deformation of each set using the least-squares principle, and then solves the least-squares solution in the sense of the smallest parametric value for all small baseline sets jointly, to the time-series surface deformation covering the entire monitoring period using the method of singular value decomposition. This technique interferometrically processes SAR images with a short time and space, which effectively avoids temporal and spatial decoherence of radar interferometric processing in the Badong Mountain area, increases the temporal sampling rate, improves the quality of interferograms, and also improves the comprehensive utilization rate of data in the current project with fewer time-series radar images, and ensures the timeliness and accuracy of surface deformation monitoring.

The following are the main processes in the SBAS-InSAR approach to extracting ground subsidence (Berardino et al. 2002). The experimental data were processed by Swiss GAMMA commercial software, and we did the atmospheric correction with GACOS (Jolivet et al. 2014). Figure 4 depicts the spatial and temporal baselines for SAR data, with a spatial baseline of 100 m and an unset temporal baseline. After ensuring that data alignment is within 0.1 pixels, differential interference is performed on all image pairs (60 pairs) with an adaptive filtering window of 128, and a stable point is chosen as the starting point for the deconvolution with a deconvolution threshold of 0.3, followed by phase transformation processing to obtain the results.

Figure 4. Spatial and temporal baselines of SAR datasets.

3.3. Spatial statistical analysis of InSAR rate result

We employed hot spot analysis ($\text{Getis}-\text{Ord}\mathrm{ }{\mathit{G}}_{\mathit{i}}^{\mathrm{*}}$) and kernel density analysis techniques to perform spatial statistical analysis of InSAR deformation rate results to identify deformation aggregation locations and find strong deformation sites (Morelli et al. 2011; Lu et al. 2019; Zhu et al. 2021).

The $\text{Getis}-\text{Ord}\mathrm{ }{G}_i^{\mathrm{*}}$ statistic is a local spatial statistic that measures the degree of spatial correlation within a defined range and can be used to evaluate the degree of spatial aggregation of deformation rates. $G_i^{\mathrm{*}}$ is calculated as follows (1) $$G_i\left(d\right)=\frac{\sum _{j=1}^n w_{ij}\left(d\right)x_j}{\sum _{j=1}^n x_j}(j\ne i)$$(1) $x_j$ is the deformation rate value of the point target; $w_{ij}$ is the spatial of the weight of point targets i and j; n is the number of point targets in the total dataset. When $G_i^{\mathrm{*}}\mathrm{ }$ is positive (hot spot) and highly significant, it suggests a high value spatial aggregation; when $G_i^{\mathrm{*}}$ is negative (cold spot) and very significant, it indicates a low value spatial aggregation. The $G_i^{\mathrm{*}}$ of each point can be determined using hot spot analysis (Getis-Ord $G_i^{\mathrm{*}}$). The greater the absolute value of $G_i^{\mathrm{*}},$ the stonger the aggregation relationship of the deformation points.

Kernel density analysis is frequently used to compute the density of elements in their immediate surroundings, and the results can reflect the degree of aggregation of the input data in terms of spatial arrangement. Based on $G_i^{\mathrm{*}},$ a smooth surface can be fitted by applying kernel density analysis to generate the final deformation aggregation map, with the kernel density expression as (2) $${\mathrm{f}}_n=\frac{1}{nh}\sum _{i=1}^n k\left(\frac{X-X_i}{h}\right)$$(2) ${\mathrm{f}}_n$ is the kernel density estimate; k () is the kernel density function; n is the number of existing points; $X-X_i$ is the distance between the deformation rate point $X$ and $X_i;$ h represents the search radius. The kernel density analysis of $G_i^{\mathrm{*}}$ score of deformation points can turn the area of significant movement of the aggregation of deformation points into an area of hot and cold spots with good visibility, highlighting the deformation area and more closely matching the actual landslide movement. High absolute density readings indicate unusual surface activity in this region, which could be caused by geological risks such as landslides.

3.4. Pearson correlation coefficient

According to the statistics, rainfall is the primary cause of the Badong ecological barrier. Locating the deformation area affected by rainfall in the high landslide area can improve the accuracy of landslides identification. Pearson’s correlation coefficient method is a statistical approach to measure the proximity of the relationship between two variables, which can be used for the spatial correlation between two variables in a given spatial environment, and then to provide an objective interpretation of the relationship between the variables (Buytaert et al. 2006; Di et al. 2013; Singh et al. 2020). In this work, the correlation between rainfall and deformation is determined using Pearson’s coefficient. Several sets of data can be obtained throughout time for two variables, rainfall $\mathit{x}$ and cumulative form variable $\mathbf{y},$ denoted as $\left({\mathit{X}}_{\mathit{i}},{\mathit{Y}}_{\mathit{i}}\right)\mathrm{ }\left(\mathbf{i}=1,2,\dots ,\mathrm{ }\mathbf{n}\right),$ then the mathematical expression of the correlation coefficient is: (3) $$r=\frac{\sum _{i=1}^n \left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}{\sqrt{\sum _{i=1}^n {\left(X_i-\overline{X}\right)}^2}\sqrt{\sum _{i=1}^n {\left(Y_i-\overline{Y}\right)}^2}}$$(3)

In the formula, $\mathrm{ }x$ and $\mathrm{y}$ are the mean values of the $\mathrm{n}$ data, respectively. The correlation coefficient $\mathrm{r}$ ranges between $-1$ and $+1.$ The closer $\left|\mathrm{r}\right|$ is to 1, the higher the degree of linear correlation between $x$ and $\mathrm{y}.$ If $r=1,$ it shows that $x$ and $\mathrm{y}$ are fully linearly correlated; if $\mathrm{r}=-1,$ it indicates a completely negative linear correlation between $x$ and $\mathrm{y}.$ The degree of correlation can be categorized into the following cases based on the study scale: when $0.6\le r<$1, it is regarded as high correlation; when $0.4\le r<0.6,$ it is regarded as medium correlation; when $0.2\le r<0.4,$ it is regarded as low correlation; when $r<0.2,$ it means that the correlation between the two variables is extremely weak and can be regarded as non-linear correlation.

3.5. Landslide susceptibility mapping

Landslide susceptibility analysis is to analyze under what conditions landslides are most prone to occur, where the conditions include various characteristics, such as geological structure, topography, land use, slope structure, and meteorology and hydrology (Marjanović et al. 2011; Ghorbanzadeh et al. 2018).

Machine learning can be used to predict the probability of landslide occurrence and can achieve good results, which can provide a basis for landslide monitoring and forecasting. This experiment used the gradient boosting decision tree (GBDT) model (Yao et al. 2008; Tehrany et al. 2015; Tien Bui et al. 2016). GBDT is one of the best fitting algorithms for real-world distributions among traditional machine learning methods, with strong generalization capabilities. It can be used in both classification and regression problems, as well as in regular functions to improve training results and reduce model overfitting (Chen et al. 2020; Rong et al. 2020). Unlike the AdaBoost method, GBDT generates a weak classifier through multiple iterations, with each classifier being trained based on the residuals of the previous iteration. The weak classifiers are weighted and summed in each iteration to get the result. The GBDT model can be described as: (4) $$F_m=\sum _{m=1}^{M}T(x,{\theta }_m)$$(4) where $M$ is the number of iterations, $T(x,{\theta }_m)$ is the weak classifier generated at each iteration, and ${\theta }_m$ is the loss function, which can be described as: (5) $${\theta }_m=\mathit{argmin}\sum _{i=1}^{N}L(y_i,F_{m-1}(x_i)+T(x_i,{\theta }_m\left)\right)$$(5) $F_{m-1}\left(x_i\right)$ is the current iteration, and the GBDT algorithm builds the parameters of the next classifier by minimizing ${\theta }_m.$ Each training round attempts to lower the loss function as much as possible to reach a locally optimal solution or a globally optimal solution. GBDT is highly generalizable and well-suited for classification and prediction. The GBDT algorithm was developed in Python by the scikit-learn library (Song et al. 2018).

After the landslide susceptibility values of each pixel in the study area are calculated by machine learning models, the landslide susceptibility is generally classified into four or five categories according to the natural breaks classification method, such as (a) no susceptibility, (b) low susceptibility, (c) moderate susceptibility, (d) high susceptibility, and (e) very high susceptibility (Zhao et al. 2020; Youssef and Pourghasemi 2021).

We implemented the GBDT model in python and evaluated the prediction accuracy of the landslide susceptibility model using the receiver operating characteristic (ROC) curve, also known as the perceptibility curve. The ROC curve is a composite indicator, showing the relationship between two continuous variables. Vertical coordinate means susceptibility and horizontal coordinate is specificity (Cantarino et al. 2019; Goyes-Peñafiel and Hernandez-Rojas 2021; Sun et al. 2021; Fustos-Toribio et al. 2022). The area under the curve is referred to as AUC, and the larger the AUC, the more accurate the model is.

4. Results

4.1. InSAR deformation results

The SBAS-InSAR approach was used to monitor the deformation of the barrier zone. Figure 5(a) shows the annual average deformation rate along the radar line of direction of sight in the research region from April 2020 to August 2021, with the maximum settlement rate of −148 mm/a. The locations with larger deformation are mostly found on the south bank of the Yangtze River and the eastern section of the north bank of the Yangtze River. The average coherence value of the 60 interference pairs was 0.580, with a maximum of 0.621 and a minimum of 0.541. The overall coherence of the images was good, demonstrating that L-band SAR images were more applicable in the lush vegetation of Badong. The deformation results are distributed throughout the study region. We performed hot spot analysis and kernel density analysis of deformation rate and obtained the kernel density of deformation rate. We circled the aggregation areas and obtained 27 deformation aggregation areas. We then performed landslides identification in these 27 deformation areas (Figure 5(b)).

Figure 5. (a) Annual average deformation rate along the radar line of sight direction. The negative value represents the displacement of the target feature away from the satellite direction, which is shown in red in the figure, and the positive value represents the displacement of the target feature close to the satellite direction, which is shown in blue in the figure. (b) Nuclear density map.

4.2. Pearson correlation analysis of rainfall and InSAR deformation results

Previously, the correlation calculation between rainfall and deformation was for a single point. In this experiment, we used GIS spatial analysis to calculate the correlation between rainfall and deformation for the whole region, which improves the visualization effect and facilitates the identification of potential landslides. The rainfall station data were interpolated to obtain the rainfall distribution maps for each month, together with the monthly accumulated deformation for the corresponding months as shown in Figure 6. A raster calculation was performed for each pixel of the study area using the formula of Pearson’s correlation coefficient to obtain the correlation distribution of monthly accumulated deformation with rainfall for the period from April 2020 to July 2021 and as shown in Figure 7 the study area was classified into four categories according to the definition of Pearson coefficient. Eight of the 27 deformation aggregates are uncorrelated, seven have a low correlation, ten have a medium correlation and two have a high correlation. We performed the first filtering at the temporal level for the 27 deformation clusters. 19 clusters with a correlation greater than or equal to Low correlation were retained (No. 1, 2, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19, 23, 26).

Figure 6. (a–c) The monthly rainfall distributions for April 2020, May 2020 and July 2021, respectively. (d–f) The monthly accumulate deformation for April 2020, May 2020 and July 2021, respectively. n = 10 represents the dates of 10 SAR images.

Figure 7. Correlation distribution map. The darker the colour, the higher the correlation between the deformation here and the monthly rainfall.

4.3. Landslide susceptibility mapping results and analysis

This experiment generates a total of 1,741,207 points, of which 84,689 are landslide points and 1,656,518 are non-landslide points, based on a raster size of 12.5 m. We randomly selected 150,000 points for model training and testing (70% for training the GBDT model and 30% for verifying the model performance). A 2:3 ratio of slippery and non-slippery points (60,000 slippery points and 90,000 non-slippery points) came with the best results.

The obtained model was then used to calculate landslide susceptibility for the entire area points, yielding the largest area of the ROC curve and AUC = 0.89 (Figure 8).

Figure 8. ROC curve for GBDT model.

According to Figure 9, highly susceptibility areas for landslides are mainly concentrated along the Yangtze River in the eastern part of the study area and the western bank of the Yangtze River tributaries in the north part of the study.

Figure 9. Landslide susceptibility map.

Of the 27 deformation aggregation areas, four were in the no-susceptibility area, six in the low susceptibility area, three in the moderate susceptibility area, two in the high-susceptibility area, and twelve in the very high susceptibility area. Of the 19 temporally filtered deformation aggregates, 13 deformation aggregates that are greater than or equal to moderate susceptibility were retained. According to Figure 10, the secondary filtering was completed (No. 2, 3, 5, 6, 7, 8, 10, 11, 12, 15, 16, 19, 23).

Figure 10. The map of the results of deformation landslides identification.

4.4. Potential landslides identification results

27 deformation aggregation zones were circled using hot spot analysis of InSAR rate results. We extracted the Pearson correlation coefficients and sensitivities of the deformation aggregation zones, and the deformation zones with Pearson correlation coefficients greater than or equal to 0.2 and higher than medium susceptibility were retained. Through Pearson correlation coefficient results and susceptibility results, we screened 13 deformation aggregation areas caused by rainfall. Then those 13 areas were verified in the field, with two areas not being reached (No. 3, 5). The remaining 11 (No. 2, 6, 7, 8, 10, 11, 12, 15, 16, 19, 23) aggregation areas all had larger or smaller deformations after field verification and comparison with the survey data in 2020. Finally, after determining the locations of these 11 landslides with deformation, we refined the boundaries of the 11 landslides based on the survey data and high-resolution images and other information, and obtained the final identification results, as shown in Figure 10.

In Table 3, we give the basis for our identification, with the minimum, maximum and mean values of the annual average deformation rate. The identification results include two new landslides, No. 12 and No. 19, along with nine historical landslides. In the direction away from the satellite, the absolute values of deformation rates of Dugongci landslide (Shijiapo) (No. 8), Songjiazhuang landslide (No. 23) and Songjiawuchang landslide (No. 11) are larger, exceeding 75 mm/a.

Table 3. Identifying landslides information table.

Deformation no.	Landslide name	Susceptibility	Correlation	Latitude, longitude	Maximum LOS rate (mm/year)	Field validation description
2	Wujiayuanzi landslide	Very high	0.2–0.4	31°12′30.281″ N, 110°17′52.852″ E	−68	Debris had rolled down on the highway side of the site. After investigation, the landslide was in a less stable state in the middle of the road cut slope zone at the back edge. Debris had rolled down on the highway side of the site. GPS data showed that the monitoring point increased the accumulated displacement deformation by 0–3 mm in 2019. The horizontal accumulated deformation increased 0–35.9 mm in 2020. The landslide was currently in the peristaltic deformation stage.
6	Hulubao landslide	High	0.4–0.6	31°4′18.771″ N 110°21′8.551″ E	−48	The site was under construction. Many anti-slip piles could be seen, and the surface soil showed obvious signs of excavation, with a soil layer of about 6–8 m. The front of the house on the north side of the survey site slipped in April or May 2021. The original landslide (landslide body No. 5) was deformed in 2012. It threatened the main body of the station building under construction and the road below.
7	Pengjiawuchang landslide	Very high	0.6–0.8	31°5′3.030″ N 110°24′0.178″ E	−66	The front edge of the road, orange tree development under-seat cracks, direction 188°, about 8 m long, the maximum under-seat about 0.5 m, in a straight line, seam width 1–2 cm. There is a local collapse in the multi-stage dry masonry retaining wall. Every year, the soil layer has a downward deformation, with the largest place being about 0.5 m, according to the households. It posed a threat to the slope’s residents and the road.
8	Dugongci landslide (Shijiapo)	Very high	0.4–0.6	31°3′22.433″ N 110°18′52.963″ E	−113	The point was in the south slope of the topography of the near east-west wash trough. Under the influence of tectonics, the rock fissures were more developed and the rock is more broken.
10	Tangu landslide	Very high	0.4–0.6	31°2′26.941″ N 110°26′17.352″ E	−64	The site soil was soft and the trees were skewed.
11	Songjiawuchang landslide	Very high	0.6–0.8	31°1′54.426″ N 110°27′7.280″ E	−79	There were cracks on the road with a direction of 220°and a depth of about 40 cm. The road was repaired in 2018 and 2019 according to the occupants, but there was a new cracking. Soil would slip during the rainy season. It posed a threat to the nearby residents and the road.
12	Chapengzi landslide (new)	High	0.2–0.4	31°1′23.924″ N 110°25′21.891″ E	−39	There were landslide management grids at the site, but they have collapsed several times and caused damage to bridges. The landslide and the bridge were both repaired at the same time. Construction began in October 2020 and was completed around September 15, 2021.
15	Huangtupo landslide	Very high	0.4–0.6	31°2′20.728″ N 110°23′9.516″ E	−38	There were many settlement deformation cracks on the highway. The walls of abandoned buildings cracked. The deformation length of the road was about 10 m. Those road cracks occurred about 2–3 years ago.
16	Zhaoshuling landslide	Very high	0.2–0.4	31°2′42.249″ N 110°20′44.047″ E	−31	Landslide control grids were visible on-site and under construction.
19	Ziduanmai landslide (new)	Very high	0.4–0.6	31°0′4.550″ N 110°19′16.393″ E	−37	Collapse deposited on a large scale. The rock was soft and had a thickness of about 15–25 m and a length of about 150 m. Localized sliding and collapse could be seen at the site along the road.
23	Songjiazhuang landslide	Very high	0.4–0.6	31°0′42.626″ N 110°15′15.877″ E	−83	Roadbed collapsed because rocks in the slope body above crumbled and hit the roadbed and highway guardrail. We could see the airside of the slope. The rock was in an unstable state and may crumble down at any time.

Notes on the minimum, maximum and mean values of the annual average deformation rate: the negative value represents the displacement of the target feature away from the satellite direction, and the positive value represents the displacement of the target feature close to the satellite direction.

5. Discussion

5.1. Analysis of Songjiawuchang landslide deformation

The Songjiawuchang landslide in deformation zone No. 11, which has the larger deformation, was selected for deformation analysis (Figure 11(a)). The landslide of Songjiawuchang is located along the Yangtze River, with a slope direction of about 175°. Scholars used ground truth for results verification, e.g. GPS, levelling data, and Global Navigation Satellite System (GNSS) (Yin et al. 2010; Wasowski and Bovenga 2014; Dong et al. 2019). We selected the GPS points P1, P2, and P3 in the landslide body for the analysis, and the profiles are shown in Figure 11(b). The landslide is rocky, and the landslide body is primarily made up of landslide debris interspersed with soil, forming a multi-phase sliding accumulation. As demonstrated in Figure 11(c–e), the cumulative deformation curves of InSAR and GPS exhibit good agreement, and the trends are consistent, demonstrating the validity of the InSAR results of this experiment.

Figure 11. Songjiawuchang landslide deformation analysis figures. (a) The map of the whole landslide and the indication of the location of the A–A′ profile line; (b) the map of the Songjiawuchang landslide profile; (c–e) combined maps of GPS, InSAR, and rainfall; (f–h) geological investigation site maps.

The precipitation in the region is mostly concentrated between May and October. The highest monthly rainfall occurs in July 2020, reaching 232 mm, while the lowest monthly precipitation is less than 10 mm in December 2020. We compared the data of InSAR with GPS and found that the subsidence in the region is consistent with the trend of precipitation. The slope of the deformation curves observed by InSAR and GPS increases after the surge of precipitation in May each year, when the flood season begins, indicating that the deformation trend accelerates with the increase of precipitation, reflecting that the surge of precipitation accelerates the deformation trend of the monitored sites in the region.

From the deformation curves, the GPS and InSAR results showed that the landslide had large deformation in April to September 2020 and April to August 2021, indicating that the surface deformation in the region was significantly affected by the increase or decrease of precipitation. We should strengthen the monitoring density for the slope during the rainy season. According to the cumulative deformation of the three points, P1 has a larger displacement than P2 and P3, indicating that the higher part of the slope body has a larger displacement. The higher part is more sensitive to rainfall. It is a pushed landslide. We need to strengthen the monitoring of the trailing edge of the landslide to ensure the safety of the place.

We conducted two field surveys of the Songjiawuchang landslide during the study period and found tensile deformation in the soil of orange trees (Figure 11(f)), a 2 cm-wide road crack (Figure 11(g)), and house cracks in villagers’ homes (Figure 11(h)). We photographed the tensile cracks under orange trees and the road cracks on December 12, 2020, and photographed house cracks in the villager’s house on April 20, 2021. According to the findings of the visit, earthquakes struck the site in 2019 and 2020. There are 1 cm-wide house cracks at No. 23 of Group 6 and No. 25 of Group 7 of Green Bamboo Raft Village. June 2020 to the present, the road is being built within this landslide, which has had an impact on the road surface in front of the houses, causing cracks in the road and making them unstable. As a result, it appears that, in addition to rainfall, earthquakes and human engineering activities play a role in the Songjiawuchang landslide.

5.2. The reliability and applicability conditions of the identification method

InSAR provides direct evidence of deformation for identifying landslides. However, InSAR results have redundant information on deformation. Considering that the deformation of landslides is aggregated and non-fragmented, we introduced two spatial analysis tools, hot spot analysis, and density analysis. The results of the spatial analysis showed that the method successfully filtered some impurity information compared with the original deformation rate map, which helped us to select 27 deformation aggregation areas in the study area. The work avoided some subjectivity and helped us to retain objectivity.

We need to analyze the probability of spatial occurrence and temporal triggering factors of landslides based on the deformation zones. Therefore, we introduced landslide susceptibility analysis and Pearson correlation coefficients based on rainfall and deformation. We used those to filter and screen out the hidden dangerous areas at spatial and temporal levels. Finally, 11 landslide hazards were found at the 11 different locations. Our investigation found that all these 11 landslide hazards had either large or small deformation. The study also carefully analyzed the Songjiawuchang landslide, which has large deformation. The comparison of the curves from the landslide development and InSAR and GPS monitoring data all proved the effectiveness of the identification method.

The spatial filter, landslide susceptibility is relatively easy to obtain, based on the corresponding method. And it is generally more closely integrated with geology and does not change in a short time. On the other hand, the temporal filter must be considered carefully along with the study area. Since the study area is mainly along the Yangtze River, where the rock layer is significantly affected by rainfall and reservoir water level, a time filter method based on rainfall is more effective. Corresponding adjustments can be made for different study areas based on the main triggering factors of the corresponding area, which can include human engineering activities, earthquakes, groundwater, etc.

6. Conclusions

Identifying landslides is essential for sustainable human development, and InSAR is a good method to monitor surface deformation, especially for large-scale studies. How to explore the relationship between surface deformation and landslide is in building a bridge between InSAR and landslide identification. Therefore, this paper proposed an identification process method that integrated five means of InSAR, spatial analysis, landslide susceptibility, rainfall correlation, and field survey for identifying potential rainfall landslides. We proposed a rainfall-based landslide identification method that included temporal and spatial filtering. We used the SBAS-InSAR technique to process fourteen PALSAR-2 images acquired between April 14, 2020, and August 7, 2021, to generate ground displacement rates and identify potential landslides in the Badong barrier area.

Using InSAR for landslides identification necessitates taking into account both risk factors and hazards-pregnant background conditions. In this study, the temporal and spatial filtering of the results of the deformation aggregation area of InSAR was actually a method that takes into account the pregnant conditions and predisposing factors. Landslide susceptibility is a study to evaluate the spatial probability of landslide occurrence, which incorporates geological conditions, geomorphological factors, etc., and can realize the assessment and filtering of hazards-pregnant background conditions in the study area. In addition to background conditions, predisposing factors need to be considered. There were more rainfall induced landslides developed in the Badong ecological barrier area, so the rainfall factor was paid attention to in our research. Correlation analysis of monthly rainfall and InSAR results can further target potential impact areas. The addition of background conditions and predisposing factors to the InSAR results greatly improves the success rate of landslides identification.

Data availability statement

The data that support the findings of this study are available from the Hubei Geological Environment Station but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the authors upon reasonable request and with permission of the Hubei Geological Environment Station.

Additional information

Funding

This work was supported by the 2019 Central Natural Disaster Prevention and Control System Construction Grant (Regional Integrated Remote Sensing Monitoring and Hidden Danger Identification in Enshi Prefecture) (HBYM-GTJ-2020-GK-1002) and the Geohazard Prevention and Control in Three Gorges Follow-up Work – Guidance on Professional Monitoring Construction and Early Warning Analysis of Geohazards in Three Gorges Reservoir Area (Year 2022) (000121-2012A-C50-021).