Automatic Mapping of Potential Landslides Using Satellite Multitemporal Interferometry

Abstract

Mapping potential landslides is crucial to mitigating and preventing landslide disasters and understanding mountain landscape evolution. However, the existing methods to map and demonstrate potential landslides in mountainous regions are challenging to use and inefficient. Therefore, herein, we propose a method using hot spot analysis and convolutional neural networks to map potential landslides in mountainous areas at a regional scale based on ground deformation detection using multitemporal interferometry synthetic aperture radar. Ground deformations were detected by processing 76 images acquired from the descending and ascending orbits of the Sentinel-1A satellite. In total, 606 slopes with large ground deformations were automatically detected using hot spot analysis in the study area, and the extraction accuracy rate and the missing rate are 71.02% and 7.89%, respectively. Subsequently, based on the high-deformation areas and potential landslide conditioning factors, we compared the performance of convolutional neural networks with the random forest algorithm and constructed a classification model with the area under the curve (AUC), accuracy, recall, and precision for testing being 0.75, 0.75, 0.82, and 0.75, respectively. Our approach underpins the ability of interferometric synthetic aperture radar (InSAR) to map potential landslides regionally and provide a scientific foundation for landslide risk management. It also enables an accurate and efficient identification of potential landslides within a short period and under extremely hazardous conditions.

1. Introduction

Landslides are widespread natural disasters that occur worldwide, causing extensive losses in terms of economic and natural resources as well as human lives [1]. Long-term records suggest that landslide disasters increase due to climate changes and population growth [2]. Hence, accurate prediction of landslide occurrence is crucial to reduce the damages and losses caused by landslides. According to the entire process of a slope failure, the process can be divided into three stages: stable stage (stable slope), weak-stable stage (potential landslide), and failure stage (landslide failure) [3]. Potential landslides, including unstable slopes and reactivated ancient landslides, are the essential evolutionary stage of landslide failure [4]. Therefore, identifying and mapping potential landslides is the premise of landslide monitoring and early warning at a regional scale.

Field investigation based on expert knowledge is the fundamental method among the various existing methods for identifying potential landslides [5]. The obvious deformation signs (e.g., cracks, scarps, and infrastructure damage) can be obtained on the potential landslides, reflecting the kinematic pattern of slopes [6,7]. The topographic, geomorphic, geological, and hydrological conditions of slopes were manually analyzed in the field to determine the slope state by summarizing the relationship between these condition factors and failed landslides [5]. However, manual investigation for large-scale or steep mountainous regions is unfeasible, and the identification result is subjective, relying on expertise [8]. Landslide susceptibility assessment enables the construction of the relationship between landslide condition factors and landslide occurrence based on spatial data and a geographic information system (GIS) [9] using physics-based models, viewpoint-driven models, statistical models, and machine learning models [10,11,12]. Machine learning is an ideal technique for addressing big-data spatial analyses and can solve non-linear geo-environmental issues [13]. However, the selection of optimal models and preparation of spatial databases regarding landslides and conditioning factors may influence the accuracy, reliability, and applicability of the assessments [14]. In addition, based on the landslide susceptibility result, a reliable spatial location and extent of potential landslide identification is still a challenging task [13]. In order to identify potential landslides quantitatively, the stability analysis of slopes was employed [3,15]. Because the geometry of the slope, physical properties of the material composition, and failure mechanism of the slope are considered [16], the obtained stability number is accurate. However, the method is unsuitable for large-scale monitoring and the calculated process is complex. Therefore, using earth observations to detect and map reliable potential landslides on a large scale has become necessary.

In recent years, the use of multitemporal interferometric synthetic aperture radar (multitemporal InSAR) has become the most efficient way to detect ground deformations and map potential landslides at a regional scale because of the cost effectiveness of the technique and its ability to cover vast areas and handle high-frequency and high-precision surface displacement measurements [17,18]. Furthermore, the multisource synthetic aperture radar (SAR) data and the satellite images obtained from the Sentinel-1 mission of the European Space Agency help perform multitemporal interferometry (MTI), focusing on the mapping of potential landslides [19]. However, mapping potential landslides based on surface displacement needs visual interpretation to delineate the boundaries of the high-deformation areas and field investigations to verify the potential landslides, which is complex and time-consuming [20]. Although some research has proposed methods to identify potential landslides using InSAR deformation results [20], these methods focus on the relationship between deformation and potential landslides, which neglects the effect of geological factors. There are still challenges restricting the efficiency of potential landslide mapping at the regional scale, including (i) mapping slopes with large deformation automatically instead of visual interpretation and (ii) demonstrating that the ground deformation was caused by potential landslides.

Thus, this study aimed to develop an approach combining remote sensing, GIS, and deep learning (DL) methods. We chose Guide-Xunhua Basin in the upper reaches of the Yellow River as our study area, which is prone to severe geohazards [21]; therefore, the number of potential landslides in the area is high. We used MTI and 76 Sentinel-1A data acquisitions to detect ground deformation and conducted a detailed field survey to collect the data required to prepare the potential landslide inventory. We performed hot spot analysis to generate slope polygons with large ground deformations. The potential landslide inventory and a convolutional neural network (CNN) were used to construct a classification model for potential landslide identification.

2. Study Area

The upper reaches of the Yellow River, located at the transitional zone between the Loess Plateau and the Qinghai-Tibet Plateau, have experienced approximately 508 historical landslides [22]. The terrain inclines from the northwest to the southeast, and the landform of the study area is composed of hilly, mountain, plain, and intermountain basin, and is layered with altitude ranging from 4421 m to 1837 m [23]. Because of the uplift of the Qinghai–Tibet Plateau and the erosion of the Yellow River and its tributaries, high-relief and steep slopes have been formed in the area. The tectonic framework of the region is characterized by NW–SE-trending faults dominated by the Lajishan Fault and Xi Qinling-Jishishan Fault [23]. The study area comprised a 15 km buffer zone along both sides of the Yellow River from Guide County to Xunhua County (Figure 1).

The lithology of the study area is mainly Paleogene and Neogene sandstones and mudstone and the Quaternary sediments, characterized by heavy weathering [24]. The climate of the study area is semi-arid with a dry and windy spring, a short and cool summer, a wet and rainy autumn, and a long and dry winter [25]. Rainfall is mainly experienced in summer and autumn with an annual average precipitation of ~342 mm, an evaporation rate of ~1689 mm, and low vegetation development on both sides of the banks [24,25].

Therefore, the interplay of neotectonics activities, the erosion of the Yellow River and its tributaries, high relief and steep slopes, and cold and semi-arid climate make the study area susceptible to debris flows and landslides [22]. In addition, many giant landslides and reservoir landslides developed in the study area, like Xijitan landslide [26], Garang landslide [27], and Lijia Gorge Ⅲ landslide [23], which have seriously threatened the lives of residents and the safe operation of hydropower stations.

3. Data and Methods

In this paper, the used approach is illustrated where the MTI technique, hot spot analysis, and CNN algorithm are combined to construct a model for potential landslides classification as shown in Figure 2. Firstly, the surface ground deformation rate was obtained from descending and ascending orbits using MTI techniques. Subsequently, hot spot analysis was employed to detect high ground deformation areas automatically. The high ground deformation areas were determined as modeling databases based on field investigation and remote sensing images. The potential landslide classification model was constructed to analyze the relationship between selected factors and potential landslides based on random forest (RF) and CNN algorithms. Finally, the performance of classification models is evaluated and compared, and a new approach to automatically classify and map potential landslides can be constructed.

3.1. Ground Deformation Detection Using MTI

MTI can be used with multiple SAR images to detect surface ground deformations and identify potential landslides, which can be applied to long time scales with low temporal decorrelation [19]. We used interferometric point target analysis (IPTA) in our study. IPTA can identify the target points with high quality and a particular density [28]. These target points maintain stable scattering characteristics for a long time and are called persistent scatterers (PSs). Spectral diversity and amplitude dispersion can be used to select the PSs. After the initial selection of PS candidates, the differential interferograms are calculated by simulating the unwrapped interferometric phase based on the initial baselines and the available DEM. Then, a two-dimensional regression analysis is taken and height corrections, linear deformation rates, quality measures, residual phases, and the unwrapped interferometric phase are obtained to improve the phase model. Finally, the atmospheric phase delay can be derived from the residual phase using spatial and temporal filtering, and the non-linear deformation is obtained [29,30]. The core idea of IPTA is the stepwise and iterative improvement of different parameters, including baseline refinement, atmospheric phase filter, height correction, and point list extension [28].

Herein, 39 descending Sentinel-1A images of path 135 taken between 16 January 2020 and 28 April 2021 and 37 ascending Sentinel-1A images of path 128 taken between 4 January 2020 and 28 April 2021 were collected from the European Space Agency. The SAR image coverage of the descending and ascending orbit datasets were shown in Figure 1a. Then, 71 (69) interferograms were generated with a perpendicular baseline threshold of 160 m and a temporal baseline of 30 days in the descending (ascending) orbit dataset (Figure 3). A reference point was selected to remove the residual terrain phase, which had a coherence value > 0.8 and was located within a stable region as determined by field investigation (Figure 1). The 1 arc-second (30 m) shuttle radar topography mission data obtained from the United States Geological Survey were used to remove the topographic phase and geocode the InSAR products. Precise orbit data from the European Space Agency were used to correct the orbit error. In the study, we processed the IPTA technology using GAMMA software (v20210701).

3.2. Hot Spot Analysis

Hot spot analysis methods including Getis-Ord Gi* statistics and kernel density estimation were used herein. The Getis-Ord Gi* statistics can be used to evaluate the spatial association within a specified distance of features [31]. If a feature point has a high (low) value and is surrounded by other feature points with high (low) values within a certain distance, the statistically significant hot (cold) spot can be obtained. This process is the same as manually identifying the potential landslides based on InSAR results. For each single CT at a site i and other CT at a site j, the Gi* can be calculated as follows [32]:

$${\mathrm{Gi}}^{*}=\frac{{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i},\mathrm{j}}*{\mathrm{x}}_{\mathrm{j}}-\stackrel{\mathrm{-}}{\mathrm{X}}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i},\mathrm{j}}}{\mathrm{S}*\sqrt{\frac{[\mathrm{n}{\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i},\mathrm{j}}^2-{\left({\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i},\mathrm{j}}\right)}^2]}{\mathrm{n}-1}}}$$

(1)

where ${\mathrm{x}}_{\mathrm{j}}$ is the line-of-sight (LOS) velocity of each interferometric point target (IPT); ${\mathrm{w}}_{\mathrm{i},\mathrm{j}}$ is the spatial weight with 1 for all CTs within a specified distance of CT i and 0 for other CTs including the CT i itself; n is the total number of IPTs within the scale distance d; $\stackrel{\mathrm{-}}{\mathrm{X}}$ is the mean value of the entire IPT velocity; and S is the standard deviation of the whole IPT velocity. After Getis-Ord Gi* statistics, the Z score (standard deviation) and p-value (independence probability) of each IPT can be calculated. The larger the Z score, the more intense the clustering of high values [31].

To better visualize the location of hot spots, the kernel density estimation was performed taking the Z score as the weighting factor. The kernel density estimation calculates the density of the features in their neighborhood by fitting a smoothly tapered surface. The kernel density estimation can be expressed as follows [33,34]:

$$\mathrm{Density}=\frac{1}{{\mathrm{d}}^2}\sum _{\mathrm{i}=1}^{\mathrm{n}}[\frac{3}{\mathsf{\pi }}{\mathrm{p}}_{\mathrm{i}}{(1-{\left(\frac{\mathrm{x}-{\mathrm{x}}_{\mathrm{i}}}{\mathrm{d}}\right)}^2)}^2]$$

(2)

where d is the search radius; $\mathrm{x}-{\mathrm{x}}_{\mathrm{i}}$ is the distance from each calculating pixel to IPT i; n is the total number of IPTs; and ${\mathrm{p}}_{\mathrm{i}}$ is the weighting factor for each IPT. Thereafter, a smooth density map, which converts a large amount of IPTs into several hot spots, was obtained and the mass movement area was highlighted using the Spatial Statistics Tools in ArcGIS v10.6.

3.3. Modeling Algorithms

Conventional machine learning (RF) and deep learning (CNN) algorithms were applied in the current study to map potential landslides. These techniques are briefly described in this section.

3.3.1. CNN Model

A CNN is a biologically inspired DL technique that demonstrates outstanding performance for landslide prediction modeling and feature extraction, unlike traditional machine learning [35]. The architecture of the used CNN architecture is shown in Figure 4, which refers to the LeNet-5 [36]. The novel CNN model comprises multiple layers of neurons, including input, convolution, pooling, fully connected, and output layers [37,38]. The input layer comprises several neurons, each representing a modeling factor [39]. The convolution layer can extract various feature maps related to the target using convolutional kernels with different weights and biases [35,40]. The activation function following the convolution layer is used to present nonlinear relationships, and the rectified linear unit (ReLU) function is the most common and effective activation function [41]. The pooling layer can reduce the dimensionality of feature maps to decrease the number of parameters and computational complexity. The widely used pooling operation is max pooling, which can compute the peak values of the local unit patches in the feature maps while the pooling operation retains the essential information [37]. The fully connected layer is connected to the output layer, which comprises two neurons and functions as a classifier that categorizes landslides as potential and non-potential landslides [35]. All the parameters in the CNN algorithm are optimized based on the ideas of back propagation and stochastic gradient descent, which are iteratively updated to minimize the loss value [37].

This study implements the CNN model based on Kears with TensorFlow as a backend. In the CNN structure, the input layer can be made into a 14 × 1 array format, which refers to the number of input features. The number of channels of the convolutional layer was set to 14, which means that 14 feature maps can be extracted. The kernel size of the convolutional layer was set to 2 × 1. A max-pooling layer with the pooling size of 2 × 1 is added after the convolutional layer. The fully connected layer with 20 neural units follows the previous pooling layer. In addition, a dropout operation at a rate of 0.2 is employed in the fully connected layer to reduce any substantial overfitting [42]. The Adam optimizer is employed with standard b values of 0.01, and the learning rate was set as 0.001. The categorical cross-entropy loss function is applied with the number of epochs set as 500.

3.3.2. RF Model

RF is an ensemble learning method of training multiple independent decision trees based on different data subsets and combining several decision trees for classification [43]. The operating process of the RF algorithm is shown in Figure 5. First, a large number of unrelated training subsets are constructed using the bagging sampling technique to select samples in playback randomly [44]. Based on the training subsets, n decision trees can be established, which combine into the final random forest. The output of the random forest is determined by decision tree voting, which can avoid the weak classification of minority trees. RF model has shown strong robust and accurate performance in processing complex data [45].

In this study, the Scikit-Learn package in Python was used to develop the RF model. The RF model is composed of 50 binary decision trees, which are constructed using the classification and regression tree (CART) model [46]. The Gini coefficient minimization is used as a criterion to split nodes, and the depth of the CART is set as 3 to avoid overfitting. RF has a high tolerance for outliers and noise, and high prediction accuracy and stability [44]. Therefore, we selected RF as one of our models to compare the performance of potential landslide identification with CNN.

3.4. Modeling Factors

The selection of factors for modeling depends on the assumption that the development of potential landslides is subjected to the same environments as the previous landslides [47]. According to field investigation, data analysis, and previous researches, multiple modeling factors including geomorphologic, topographic, hydrological, geological, and human activity are determined [48,49,50]. Hence, potential landslide conditioning factors were selected as initial modeling factors including slope, aspect, plan curvature, profile curvature, elevation, topographic wetness index (TWI), stream power index (SPI), topographic position index (TPI), terrain ruggedness index (TRI), lithology, land use, InSAR time series trend, and distance to faults, rivers and roads. It is well known that rainfall is an important factor that triggers landslides. However, this factor is not considered in this study because the variability of rainfall data along the study area is not clear as the study area covers only a small area.

3.4.1. Discretizing of Continuous Data

The discretizing continuous data can effectively and substantially increase the performance of the model [51]. ChiMerge is a supervised method for continuous data discretization [52]. In the algorithm, each distinct value of a continuous variable is assumed as an independent interval, and then χ² statistic is tested for whether the adjacent intervals to be merged or not. If the χ² statistic for adjacent intervals is smaller than the predefined threshold, adjacent intervals are merged because they are assumed statistically similar [53]. In the study, this algorithm was applied to automatically determine the best threshold value which can reduce the impact of statistical noise and increase the stability of the classification model [54]. So, the continuous causal factors (elevation, slope, and so on) were discretized into small intervals based on the ChiMerge algorithm.

3.4.2. Topographic and Geomorphic Factors

In this study, the Advanced Land Observing Satellite (ALOS) Digital elevation model (DEM) with 12.5 m resolution, downloaded from the Alaska Satellite Facility (ASF) Distributed Active Archive Center (DAAC) was used to create the topographic factors. The slope, aspect, elevation, plan curvature, profile curvature, and TPI were extracted through spatial analysis with GIS software.

Slope is often considered as one of the most important factors affecting slope stability [55], which can affect the stress distribution of the slope [48]. As the slope increases, the driving forces, and the probability of landslide hazard also increase. The slope map of the study area is shown in Figure 6a. Based on the ChiMerge algorithm, the slope was divided into eight classes: 0–2°, 2–9°, 9–15°, 15–19.5°, 19–21.6°, 21.6–21.8°, 21.8–39°, and >39°.

Aspect indicates the direction of the slope [56], reflecting the amount of exposure to sunlight, soil moisture, and precipitation, which leads to different land cover, water evaporation, and weathering [57]. The aspect map produced from DEM is shown in Figure 6b. The aspect was divided into eight classes: north (0–22.5°, 337.5–360°), northeast (22.5–67.5°), east (67.5–112.5°), southeast (112.5–157.5°), south (157.5–202.5°), southwest (202.5–247.5°), west (247.5–292.5°), and northwest (292.5–337.5°).

Elevation is another factor for understanding the landslide behavior, which is the general altitude of land above the mean sea level [58]. The degree and type of erosion and the concentration of human activities always vary with elevation, which causes landslides to spread in a certain range of elevations [59]. The elevation ranges from 1837 m to 4421 m as shown in Figure 6c. A total of eight categories were obtained: <1990 m, 1990–2354 m, 2354–2362 m, 2362–2447 m, 2447–2473 m, 2473–2493 m, 2493–4277 m and >4277 m.

Plan curvature is defined as the curvature of a contour line formed by the intersection of the terrain surface with a horizontal plane which is used to measure the shape of the curve [55], and profile curvature is defined as the curvature of a kinked line formed by intersecting the earth’s surface with a vertical plane, which is the direct reflection of geometric characteristics of slope profile [49]. The former usually has a direct impact on the convergence and dispersion of surface runoffs [60] while the latter material has an impact on deposition by managing the acceleration or deceleration of these materials on slope [50], which were selected as landslide conditioning factors as shown in Figure 6d,e. Plan curvature was arranged into eight categories: <−5.39 × 10⁻³, −5.39 × 10⁻³–−2.38×10⁻⁴, −2.38 × 10⁻⁴–6.17 × 10⁻⁵, 6.17 × 10⁻⁵–1.01 × 10⁻⁴, 1.01 × 10⁻⁴–9.6 × 10⁻⁴, 9.6 × 10⁻⁴–1.72 × 10⁻³, 1.72 × 10⁻³–4.3 × 10⁻³, and >4.3 × 10⁻³. In the same way, eight categories of profile curvature values were generated: <−8.6 × 10⁻³, −8.6 × 10⁻³–−2.7 × 10⁻³, −2.7 × 10⁻³–−1.8 × 10⁻³, −1.8 × 10⁻³–−1.7 × 10⁻³, −1.7 × 10⁻³−3.7 × 10⁻⁵, 3.7 × 10⁻⁵−6.9 × 10⁻⁴, 6.9 × 10⁻⁴−3.5 × 10⁻³, and >3.5 × 10⁻³.

TPI was used to describe the slope position of a topographic feature, which can indicate the difference between the focal cell elevation and the mean elevation of all cells in the neighborhood [61]. Positive values typically represent the location higher than surrounding areas, values near zero are flat areas and negative values represent lower areas [62]. The TPI has been calculated as follows [63]:

$$\mathrm{T}\mathrm{PI}={\mathrm{M}}_0-\frac{{\sum }_{\mathrm{n}-1}{\mathrm{M}}_{\mathrm{n}}}{\mathrm{n}}$$

(3)

where ${\mathrm{M}}_0$ is the elevation of the middle point, ${\mathrm{M}}_{\mathrm{n}}$ is the elevation of the grid, and n is the total number of points employed in the evaluation. The TPI map is shown in Figure 6f and the TPI values were reclassified as being <−11.4, −11.4–−4.5, −4.5–−4.2, −4.2–−3.9, −3.9–−3.5, −3.5–−0.85, −0.85–6.2, and >6.2.

TRI was used to illustrate the amount of elevation difference between the surface and surrounding area [64], which can describe the effect of erosion in the hydrological processes. The TRI was considered as a morphometric measure which describes the heterogeneous condition of a land surface [65]. The TRI was calculated using the following equation [66]:

$$\mathrm{TRI}=\sqrt{\left|\mathrm{x}\right|({\max }^2-{\min }^2)}$$

(4)

where max and min represent the highest and lowest values within a 3 $\times$ 3 cell window, respectively. The TRI map is shown in Figure 6g. The TRI values were reclassified as being <0.4, 0.4–1.6, 1.6–2.2, 2.2–3.3, 3.3–4.1, 4.1–5.2, 5.2–6.5, and >6.5.

3.4.3. Hydrological Factors

Hydrological factors also play an important role in the occurrence of landslides. The change in hydrological environment will change the pore water pressure in the slope and reduce the shear strength of the slope. At the same time, variations of hydrological conditions will continuously erode the slope, and make the slope unstable.

SPI can be employed to measure the erosive power of flowing water on the surface [67]. A stream with a high flow power will erode the surrounding slopes much more along the direction of flow. As a result of the abrasion of the slope, the stability of the slope will deteriorate and the slopes will become more sensitive to landslides [68]. SPI was computed using the following formula [69]:

$$\mathrm{SPI}={\mathrm{A}}_{\mathrm{S}}\times \tan \mathsf{\beta }$$

(5)

where ${\mathrm{A}}_{\mathrm{S}}$ and $\mathsf{\beta }$ are the specific basin area and the local slope gradient (in degrees), respectively. The SPI map is shown in Figure 6h. The SPI values were divided into eight classes, such as <1.4, 1.4–20.6, 20.6–45.0, 45.0–1397.3, 1397.3–2523.5, 2523.5–2730.0, 2730.0–31,158.0, and >31,158.0.

TWI is another vital landslide conditioning factor used for the assessment of the wetness at any given point of the catchment area [70]. TWI was computed using the following formula [71]:

$$\mathrm{TWI}=\ln \left(\frac{\mathsf{\alpha }}{\tan \mathsf{\beta }}\right)$$

(6)

where $\mathsf{\alpha }$ is the cumulative upslope area draining through a point (per unit contour length) and $\tan \mathsf{\beta }$ is the angle of slope at the point. The TWI is a physically based index of the effect of local topography on rainfall runoff patterns, areas of potentially increased soil moisture, and ponding areas [72]. TWI can be used to reflect the spatial distribution of soil saturation in response to a rainfall of a specified duration [73]. The TWI map is shown in Figure 6i. And the TWI values were reclassified as being: <4.2, 4.2–5.8, 5.8–5.9, 5.9–6.4, 6.4–6.8, 6.8–6.9, 6.9–10.6, and >10.6.

The distance to rivers is regarded as one of the most crucial factors in landslide occurrence. Rivers contribute to the instability of the slopes by scouring the slope toes and saturating the materials on the slopes [48]. It can be perceived that the relationship between distance to rivers and slope stability has a negative correlation [74]. The river distribution map of the study area was extracted from Google Earth high-resolution images and then was digitized in ArcGIS (Figure 6j). The distance to rivers in the study area was grouped into eight classes, including 0–41.4 m, 41.4–85.0 m, 85.0–1657 m, 1657–1688 m, 1688–2041 m, 2041–5449 m, 5449–10,154 m, and >10,154 m.

3.4.4. Geological Factors

The geological factors including lithology and distance to faults were considered as landslide conditioning factors to reflect the effect of the geological process on landslide development.

The lithology of an area will have a significant effect on landslide occurrence [75]. The obvious distinctions of lithological units were reflected in physical and mechanical characteristics including type, strength, degree of weathering, durability, density, and permeability [76]. In terms of the lithology, a total of thirteen groups were obtained in the study area according to the 1:500,000 geological map from the China Geological Survey. According to the rock weakness, the lithology was divided into four classes: hard (granite, gneiss of Proterozoic Era), sub-hard (sandstone, siltstone, and slate of Triassic, and sandstone, siltstone, and mudstone of Cretaceous), sub-soft (glutenite, sandstone, siltstone, and silty mudstone of Paleogene and Neocene), soft (sediments of Quaternary) [77] (Figure 6k).

The presence of structural discontinuities including faults, fractures, folds, joints, and shear zones plays an important role in decreasing the rock strength and causing landslides [78]. In the study area, the distribution of faults includes Lajishan Fault, Jianzhadong Fault, Dehenglong Fault, Songba Fault, Zhamashan Fault, Wendudasi Fault, and Xi Qinling-Jishishan Fault [79].

Therefore, the distance to faults was regarded as a potential indicator of the landslides. For this research, the faults distribution map (Figure 6j) was also obtained from a 1:500,000 geological map and the distance to faults in the study area was grouped into eight classes, including 0–6.6 m, 6.6–1208 m, 1208–2001 m, 2001–2073 m, 2073–7768 m, 7768–9542 m, 9542–13,071 m, and >13,071 m.

3.4.5. Human Activity Factors

In some mountainous regions, the establishment of roads is accompanied with intense engineering activities such as cutting slopes or removing vegetation which lead to the change in original geological conditions and weakening the natural support of rock slopes [76]. In this study, the distance to road was used to reflect the influence of constructing roads. For this research, the road lines were digitized from Google Earth high-resolution images into GIS (Figure 6j). According to the ChiMerge algorithm, the distance to roads in the study area was grouped into eight classes, including 0–766 m, 766–795 m, 795–1227 m, 1227–1304 m, 1304–2352 m, 2352–4604 m, 4604–10,983 m, and >10,983 m.

Land use has been a significant predisposing factor for landslide prediction because different vegetation has different root cohesions and hydrological conditions [80], which can prevent soil erosion and enhance slope stability. The land use map was extracted from the Environmental Systems Research Institute (Esri) 2020 Land Cover in 10 m resolution [81], and the land use map of the study area comprises eleven major groups including water, trees, grass, flooded vegetation, crops, scrub/shrub, built area, bare ground, and snow/ice (Figure 6l).

3.4.6. IPTA Time Series Deformation Trend

The time series deformation trend can represent different styles of ground deformation which are related to different slope movement processes and can also provide useful information on slope dynamics [82,83]. In this study, PS-Time, a free statistical software for time series analysis [84], was used to classify the deformation trend into three types: uncorrelated, linear, and non-linear. Uncorrelated time series trends usually denote random fluctuations of displacements around zero and typically indicate no significant movement [84]. A linear trend presents ground displacements at a constant rate, which can indicate the processes of creep, slope erosion, and natural subsidence [84]. A non-linear trend usually represents a change in the displacement rate which can be regarded as a consequence of landslide processes [82] and should be paid more attention. For the research, the proportions of each trend in a high-deformation area were used as modeling factors.

4. Results

4.1. Ground Deformation Detection Using InSAR Technology

Figure 7 shows the mean displacement rate map obtained from the datasets collected from the descending and ascending orbits of Sentinel-1A. We optimized the coherence threshold (0.7) to balance the coverage density and ground deformation quality of the IPTs. The IPTs (930,193) obtained from the descending orbit covered the study area (approximately 3516 km²) with an average density of 264 IPTs km⁻², and the IPTs (811,419) obtained from the ascending orbit covered the study area with an average density of 230 IPTs km⁻². The displacement was along the direction of radar line-of-sight (LOS) with the negative values (marked in red) indicating that the surface moved away from the satellite and the positive values (marked in blue) indicating that the surface moved toward the satellite.

Previous studies have demonstrated that the statistical characteristics of all IPT displacement rates indicate the stability threshold of potential landslide mapping. The velocity thresholds used to distinguish potential landslides are case-specific and depend on the mechanical properties of the failed material, failure mechanism, sensor measurement precision, and investigation objectives [19]. Herein, approximately 95% of the measured displacement rates were between −10 and +10 mm/y with the largest displacement rate exceeding ±300 mm/y. Most of the slopes with large deformations, which indicate the reactivation of many historical landslides, were distributed in the mountain areas of the Qunke–Xunhua section [25,26]. We demonstrated the results of MTI measurements using a ground movement survey conducted on slopes with large displacement rates. The results show that ground movements (e.g., cracks, fissures, scarps, and fragmented surfaces) are rarely discovered on slopes with displacement rates below ±10 mm/y and that most of the ground movements investigated are on slopes with displacement rates exceeding ±10 mm/y. Therefore, the stability threshold was taken as ±10 mm/y when mapping potential landslides.

4.2. Extraction of Potential Landslide Candidates via Hot Spot Analysis

We performed hot spot analysis to process the ground deformations obtained using MTI and automatically map the areas with clustered IPTs whose displacement rates exceeded ±10 mm/y. First, we calculated the Gi* index taking the displacement rate of the IPTs as the weighting factor, which could reflect the clustered and random patterns of IPTs with large displacement rates [33,85]. In order to obtain an easier and more straightforward visualization of the clustering of high positive and low negative velocity, we estimated the kernel density of the Gi* index to produce hot spot maps indicating slopes with ground deformations exceeding the stability threshold, which all represent the location of high ground deformation (Figure 8).

After that, we generated the potential landslide candidates (PLCs) using the classification tools of ArcGIS to convert each cell of a raster into a polygon. Finally, 606 PLCs were mapped based on ground deformations obtained from both the ascending and descending orbits of Sentinel-1A (Figure 9). In hot spot analysis, the search radius setting required to calculate the Gi* index, which determines the boundary between stable and unstable areas, is crucial. Unsuitable values could produce large or fragmented polygons. The search radius was 200 m, considering that the study area had experienced primarily medium and large historical landslides and that the medium spatial resolution of Sentinel 1A images was approximately 14 m × 4 m. The method is effective for extracting PLCs when the number of IPTs in an area is large.

In order to evaluate the extraction accuracy of high-deformation areas quantificationally, the visual interpretation results of high-deformation areas based on field surveys and Google Earth images are compared with the hot spot analysis extraction results. The extraction accuracy and missing rate of overall PLCs are 71.02% and 7.89%, respectively. Figure 10a shows the extraction result of the Xijitan landslide, which is partly reactivated. In the visual evaluation, most reactivated areas have been extracted, with the extraction accuracy rate and missing rate being 86% and 2%, respectively. Figure 10b shows the extraction result of the unstable slope where the scarp and secondary landslides have been developed, where the extraction accuracy rate and missing rate are 90% and 9.8%, respectively. The extraction error is induced by the boundary discontinuity, but the majority of the high-deformation areas has been extracted, which can satisfy the extraction requirements. Therefore, this approach can extract the location of high-deformation areas with high accuracy, which is an effective way of extracting high ground deformation areas from a large number of target points.

4.3. Modeling Process

4.3.1. Preparation of Potential Landslide Inventory

Ground deformation generated using InSAR is not the only factor that can be used to assess the state of slopes and map potential landslides. The slow millimeter to centimeter (per year) ground surface deformations may arise from different causes and may not reflect shear movements or landslide occurrences [19,86]. Thus, the polygons mapped using hot spot analysis should be accurately interpreted in the investigation of potential landslides. The detailed field survey of each polygon has been effective in many previous studies [18]. However, its efficiency is low, and it will not be appropriately useful in an emergency.

The successful application of the classification model depends on the availability of a detailed and reliable potential landslide inventory map. Thus, the state of the PLCs mapped using hot spot analysis should be accurately determined to train and test the classification model. Combining field investigation and optical image interpretation, obvious deformation signatures were identified as evidence to judge the state of PLCs. As shown in Figure 11a1 and a2, obvious cracks and scarps were observed on the road at the trailing edge of the slope, which indicates the active state of the slope. Some ancient landslides also show reactivated features of secondary landslides and fresh cracks (Figure 11b1 and b2). Figure 11c1 and c2 show PLCs with no suitable terrain conditions for the development of landslides [87]. The surface deformation may be caused by human activities of land cultivation and engineering construction, so these areas cannot be identified as a potential landslide. Although some PLCs satisfy the terrain conditions for landslide development, obvious deformation signatures were not observed in the field investigation, and the high deformation was induced by natural slope erosion (Figure 11d1 and d2).

We investigated 299 potential landslides, of which 35 were reactivated large ancient landslides and others were on unstable slopes in the study area (Figure 12). Amounts of 224 potential landslides and 213 non-potential landslides were selected for use in the modeling dataset, and 70% of the completed inventory was set as the training sample, while the remaining 30% was set as the testing sample. Other potential and non-potential landslides were included in the validation dataset used to optimize the model parameters.

4.3.2. Selection of Conditioning Factors

Before constructing the model, selecting proper positive landslide conditioning factors is essential because factors with noise or poor quality may negatively impact potential landslide spatial prediction results [88]. In the present case study, the multicollinearity and prediction abilities of the selected conditioning factors were quantified by applying the Pearson correlation coefficient [89] and Information Gain Ratio (IGR) [90], respectively.

The Pearson correlation coefficient can be used to measure the linear correlation between two variables. The Pearson correlation coefficient can be calculated as follows [89]:

$$\mathsf{\gamma }=\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{x}}_{\mathrm{i}}-\stackrel{\mathrm{-}}{\mathrm{x}}\right)\sum _{\mathrm{j}=1}^{\mathrm{n}}\left({\mathrm{y}}_{\mathrm{j}}-\stackrel{\mathrm{-}}{\mathrm{y}}\right)}{\sqrt{\sum _{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{x}}_{\mathrm{i}}-\stackrel{\mathrm{-}}{\mathrm{x}})}^2\sum _{\mathrm{j}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{j}}-\stackrel{\mathrm{-}}{\mathrm{y}}\right)}^2}}$$

(7)

where n is the sample size, ${\mathrm{x}}_{\mathrm{i}}$ and ${\mathrm{y}}_{\mathrm{j}}$ represent the variable values for ${\mathrm{X}}_{\mathrm{i}}$ and ${\mathrm{Y}}_{\mathrm{j}}$, $\stackrel{\mathrm{-}}{\mathrm{x}}$ and $\stackrel{\mathrm{-}}{\mathrm{y}}$ are the mean value of X and Y. The value of $\mathsf{\gamma }$ is between −1 and +1, the positive of which means the positive linear correlation, and the negative of which means the negative linear correlation between two factors [91]. In this study, a value of $\mathsf{\gamma }$ greater than 0.7 indicates a high level of collinearity between two factors, which should be removed from model building [92]. Figure 13 shows the correlation coefficients among the 17 initial conditioning factors. According to the correlation analysis, there appears a strong positive correlation between slope and TRI ($\mathsf{\gamma }$ = 0.90), a positive correlation between profile curvature and TPI ($\mathsf{\gamma }$ = 0.72), and a negative correlation between non-linear and linear ($\mathsf{\gamma }$ = −0.94). Therefore, the TRI, TPI, and linear trend were excluded from the conditioning factors.

The IGR was used to determine the importance of selected conditioning factors. The conditioning factors with a high information gain rate mean that they form a good prediction. Additionally, the low and null predictive ability of factors must be eliminated to obtain a result with higher accuracy [93]. Given that a training datum D consists of n input samples, n(L_i, D) is the number of samples in the training datum D belonging to the class L_i (potential landslide, non-potential landslide). The entropy needed to classify D is calculated as

$$\mathrm{Info}\left(\mathrm{D}\right)=-{\sum }_{\mathrm{i}=1}^2\frac{\mathrm{n}({\mathrm{L}}_{\mathrm{i}},\mathrm{D})}{\left|\mathrm{S}\right|}{\log }_2\frac{\mathrm{n}({\mathrm{L}}_{\mathrm{i}},\mathrm{D})}{\left|\mathrm{S}\right|}$$

(8)

The amount of information that needs to split D into (D₁, D₂, $\cdots$, D_m) regarding the conditioning factor A is estimated as

$$\mathrm{Info}(\mathrm{D},\mathrm{A})={\sum }_{\mathrm{j}=1}^{\mathrm{m}}\frac{{\mathrm{S}}_{\mathrm{j}}}{\left|\mathrm{S}\right|}\mathrm{Info}\left(\mathrm{D}\right)$$

(9)

The IGR for a certain potential landslide conditioning factor A is computed as

$$\mathrm{IGR}\left(\mathrm{D},\mathrm{A}\right)=\frac{\mathrm{Info}\left(\mathrm{D}\right)-\mathrm{Info}\left(\mathrm{D},\mathrm{A}\right)}{\mathrm{SplitInfo}\left(\mathrm{D},\mathrm{A}\right)}$$

(10)

where SplitInfo represents the potential information generated by dividing the training data S into m subsets. SplitInfo is calculated as

$$\mathrm{SplitInfo}\left(\mathrm{D},\mathrm{A}\right)=-{\sum }_{\mathrm{j}=1}^{\mathrm{m}}\frac{\left|{\mathrm{D}}_{\mathrm{j}}\right|}{\left|\mathrm{D}\right|}{\log }_2\frac{\left|{\mathrm{D}}_{\mathrm{j}}\right|}{\left|\mathrm{D}\right|}$$

(11)

After removing the TRI, TPI, and linear trend, the importance of each factor in the modeling was quantitatively calculated using IGR; the results are shown in Figure 14. The results show that the non-linear tread has the highest predictive capability with a value of 0.0689, which demonstrates that the non-linear deformation tread is associated with landslide movement. Other factors, including slope (AM = 0.0678), land use (AM = 0.0668), uncorrelated tread (AM = 0.057), and SPI (AM = 0.0412) show high predictive capability. Although the distance to faults (AM = 0.0152), elevation (AM = 0.0158), distance to roads (AM = 0.0163), and TWI (AM = 0.0165) have a low effect on potential landslide classification, the contribution of these factors is positive. Finally, the 14 potential landslide conditioning factors were selected to establish the model.

4.3.3. Model Performance Evaluation Index

Validating the performance of models is an essential step to attest to the effectiveness and scientific value of the classification models [94]. The validation samples were not used to construct the potential landslides classification model, so they were employed to verify the performance of the classification model. The accuracy, recall, precision, and the area under the curve (AUC) of the receiver operating characteristic (ROC) curve were applied to assess the performance of the potential landslide classification model.

The accuracy can measure the accuracy of the prediction result of the potential and non-potential landslide. The recall shows the percentage of potential landslides correctly classified in the potential landslide dataset. The precision implies the percentage of potential landslide of correctly predicted in the total number correctly classified datasets [95]. The accuracy, recall, and precision were calculated using the following equation [95,96]:

$$\mathrm{Accuracy}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}$$

(12)

$$\mathrm{Recall}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$$

(13)

$$\mathrm{Precision}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}$$

(14)

$$\mathrm{FPR}=\frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{FN}}$$

(15)

where True Positive (TP) means the areas classified potential landslides observed potential landslides, True Negative means the areas classified non-potential landslides observed non-potential landslides, False Positive means the areas classified potential landslides observed non-potential landslides, and False Negative means the areas classified non-potential landslides observed potential landslides. The ROC curve can reflect the corrections between false positive rate (FPR) (Equation (15)) and true positive rate (TPR) (also recall). If the ROC curve is nearest to the upper left corner, the model will have higher goodness-of-fit and excellent accuracy [97]. The value of AUC can be used to quantitatively reflect the accuracy of model performance with higher values indicating a better predictive capability of the model [44].

4.4. The Potential Landslides Classification Model

Although hot spot analysis enables to extract PLCs automatically, classifying potential landslides effectively is still difficult. Landslide susceptibility mapping (LSM) in a GIS-integrated environment can predict the locations of landslide occurrences by analyzing landslide conditioning factors [16]. Based on the assumptions that potential landslides in a typical area can develop on slopes with common geotechnical behaviors and identical conditioning factors, we applied the advantages of LSM to construct a classification model for potential landslide mapping based on PLCs produced by hot spot analysis.

We constructed a classification model for classifying potential landslides using the CNN and RF algorithms with optimized parameters to link potential landslides to the 14 conditioning factors. Testing the classification model was essential to verify its effectiveness and scientific value [94]. We assessed the model performance using the test samples from the potential landslide inventory.

The performance of the established models is shown in

Table 1. Performance of the models with the test dataset.

Models	Accuracy	Precision	Recall	AUC
RF	0.67	0.64	0.74	0.73
CNN	0.75	0.75	0.82	0.75

. Compared with the RF model, the accuracy, precision, and recall of CNN are 0.75, 0.75, and 0.82, respectively, which indicates the CNN model has stronger predictive capability on the test dataset. The ROC curve of the models is shown in Figure 15. The AUC of RF and CNN are 0.73 and 0.75, respectively. A high AUC value (0.75) indicates that the performance of the classification model is satisfactory [95].

The main aim of the study is to identify as many potential landslides as possible, because potential landslides are very likely to evolve into landslide geohazards, which would cause great losses of lives and properties, and the cost of identifying errors is very low. Recall denotes the ability of a model to find out all the positive samples, so a high recall value means that all potential landslides have been detected, regardless of how many objects have been mistakenly classified as potential landslides [98,99]. Therefore, the recall of the potential landslide prediction model being 0.82 means the performance of the model can meet our potential landslide identification purpose.

The evaluation result demonstrates that the CNN model is more suitable for establishing a potential landslide classification model. The convolution layer of the CNN model can not only consider the features of the current position, but also take into account the features within a certain range of the surrounding neighbors [100], which is consistent with the first law of geography [101], and is also the reason why the CNN performs well in many geographic problems [102].

5. Discussion

5.1. Advantages of the Approach

Reviewing existing researches about InSAR monitoring and potential landslide mapping at a regional scale, investigators require a lot of time and effort to delineate high-deformation areas and identify the potential landslides [2,18,103]. Compared with the traditional methods used in potential landslide investigation and mapping, the automatic approach we proposed combines the GIS spatial analysis, DL method, and MTI technique based on the strategy of LSM. The approach can improve the efficiency of potential landslide mapping and expand the capacity of InSAR for geohazard mitigation significantly.

On the one hand, delineating high ground deformation areas derived from InSAR results is vital in mapping potential landslides over a large scale. However, the dominant method in this process is still manual interpretation through human–computer interaction, which is time-consuming and labor-consuming [104]. Other researchers used object-based image analysis (OBIA) to extract high ground deformation areas [105], but the modeling process has high requirements regarding the segmentation scale and the used features, implying that multiple experiments and practical experience are needed to obtain optimal results [106], limiting the applicability of the OBIA method. In addition, directly using InSAR deformation may ignore the effect of sparse IPTs with large deformation velocities induced by data-processing errors, which can be excluded during visual interpretation. In this study, applying a hot spot map is similar to the characteristic of manually delineating high ground deformation areas, which only identifies clustering high ground deformation and helps decrease the effect of data-processing error.

On the other hand, the state of slopes was determined manually based on field investigation and has to consider the relationship between topographic, geomorphic, geological, and hydrological conditions and landslide occurrence [5]. LSM relies on models to establish this relationship and is used to predict the landslide-prone areas [12,13,14]. However, using a landslide inventory map before assessing landslide susceptibility in practice means that the landslide database was composed of historical landslides [107]. The landslide conditioning factors extracted from historical landslides cannot reflect the original topographic and geomorphic conditions, which would result in inaccurate prediction results of potential landslides. The proposed method absorbs the advantage of LSM, but the method focuses more on PLCs as objects, making the identification more accurate.

Therefore, we attempted to combine the advantages of field investigation, InSAR technology, and landslide susceptibility to build the potential landslide classification model based on InSAR technology. The high ground deformation areas were automatically extracted using the hot spot analysis which absorbed the advantages of InSAR technology and can save a lot of time and labor. The potential landslide classification model can express the expert experience as a deep learning algorithm taking potential landslides as a training database and can accurately and efficiently classify the state of the slope accurately and efficiently.

5.2. Limitations and Further Directions of the Approach

Some limitations regarding InSAR deformation monitoring, high ground deformation extraction, and landslide condition factors selection may affect the final identification results. To generate detailed PLCs using hot spot analysis, an accurate velocity map with sufficiently dense IPTs covering all the considered slopes must first be prepared. Due to the influence of SAR imaging geometry and microwave electromagnetic characteristics, the application of InSAR technology in landslide monitoring is still limited by factors such as geometric distortion caused by terrain fluctuations, uncorrelated noise generated by dense vegetation cover, and atmospheric delay errors caused by water vapor changes [108]. The visibility of the target area to satellite sensors depends on the combination of local slope direction and the acquisition geometry and mode of satellite sensor LOS [109]. The caused geometric distortions, including layover, foreshortening, and shadow, can be improved using ascending and descending data monitoring [110]. Many advanced InSAR processing methods and other improved methods continue to emerge, such as distributed scatterer (DS-InSAR) technology [111], stacking InSAR [112], and phase-decomposition-based persistent scatterer (PS) InSAR (PD-PSInSAR) technology [113], which enabled a higher density of measurement points in areas with complex topographic environment conditions [114]. In case of very fast displacements, the temporal coherence decreases leading to the loss of measurement points [19]. D-InSAR and offset tracking can be used to obtain deformation areas [115]. The Generic Atmospheric Correction Online Service (GACOS) atmospheric delay products provide a new approach for the atmospheric correction of repeat-pass InSAR, which can improve the accuracy of InSAR deformation monitoring [116]. Moreover, the L band has more penetration than the C band in landslide deformation monitoring under high vegetation cover, which can better resist the effects of decoherence and obtain higher quality deformation results [117].

Although hot spot analysis is an ideal method to highlight the high-deformation area, the extracted boundaries are incomplete with most high-deformation areas extracted, while the boundaries are inconsistent with the manual delineation. Relying on the surface geomorphological features (e.g., scarps, sliding masses, and bulging toes) reflected on optical images [118], the boundaries of high-deformation areas can be determined. Semantic segmentation deep learning model is a typical application in computer vision that aims to annotate every pixel within imagery with a specific semantic label [119]. Combining hot spot map with high-resolution optical images enables the more accurate extraction of boundaries of high-deformation areas. This process is similar to providing geomorphic features for researchers to delineate boundaries of high-deformation areas. The modeling factors should represent all the influences related to landslide development in the modeling process. Therefore, the modeling factors have to be selected after considering historical landslide mechanisms, study area characteristics, and the availability of geospatial databases. The precipitation, earthquake, or other triggering factors should be considered as critical factors for modeling [14]. Moreover, the spatial resolution of the DEM is important to derive topographic variables used in landslide classification. More precise DEM can be obtained based on unmanned aerial vehicle (UAV) or other high-precision sensors. The spatial resolution impacts the sensitivity and availability of the morphological features of slope objects and the memory storage and computational load associated with modeling. All the variables generated from the DEM are used to assess the local geomorphological characteristics of the hill slopes at a mesoscale. The balance between satisfactory model performance and practical computing time could be determined based on the magnitudes of recent landslides [16].

Additional methods for monitoring ground deformation, applying semantic segmentation deep learning models, experimenting DEM scale effects, and conditioning factors of the training model should be developed in the future. The proposed method has proven applicability in automatically mapping potential landslides and provided significant practical outcomes.

6. Conclusions

This paper presents a new approach to map regional-scale potential landslides automatically based on MTI, hot spot analysis, and DL in a 15 km buffer zone on both sides of the Yellow River from Xunhua County to Guide County of Qinghai Province, China. Our principal conclusions are as follows.

The 39 descending orbit and 37 ascending orbit Sentinel-1A images from 4 January 2020 to 28 April 2021 were processed by IPTA technique to generate the ground displacement rates to support the extraction of high deformation. The maximum displacement rate exceeded 300 mm/y was detected in the study area. The extraction accuracy rate and the missing rate of extracting the boundaries of high-deformation areas via hot spot analysis are 71.02% and 7.89%, respectively, indicating the effectiveness of hot spot analysis in this process. Combining the investigated potential landslide inventory and the selected potential landslide conditioning factors, the potential landslide classification model was constructed based on CNN and RF algorithms. The compared result demonstrated that the CNN algorithm has a better classification performance, where the AUC, accuracy, recall, and precision for testing are 0.75, 0.75, 0.82, and 0.75, respectively.

The proposed method expands the application of InSAR technology, reduces the time, economic and labor costs associated with the potential landslide mapping process, and provides a strategic basis to automatically map potential landslides in large regions. Though some problems associated with the method need to be resolved in the future, our result can provide insights for delineating potential landslides automatically and provide valuable information for landslide disasters mitigation and risk management.