Slope-Unit Scale Landslide Susceptibility Mapping Based on the Random Forest Model in Deep Valley Areas

Abstract

Landslide susceptibility evaluation is critical for landslide prevention and risk management. Based on the slope unit, this study uses the information value method- random forest (IV-RF) model to evaluate the landslide susceptibility in the deep valley area. First, based on the historical landslide data, a landslide inventory was developed by using remote sensing technology (InSAR and optical remote sensing) and field investigation methods. Twelve factors were then selected as the input data for a landslide susceptibility model. Second, slope units with different scales were obtained by the r.slopeunits method and the information value method- random forest (IV-RF) model is used to evaluate the landslide susceptibility. Finally, the spatial distribution characteristics of landslide susceptibility grade under the optimal scale are analyzed. The results showed that under the slope unit obtained when c = 0.1 and a = 3 × 10⁵ m², the internal homogeneity/external heterogeneity of 8425 slope units extracted by the r.slopeunits method is the best, with an AUC of 0.905 and an F1 of 0.908. In this case, the accuracy of landslide susceptibility evaluation is the highest as well; it is shown that the finer slope units would not always lead to the higher accuracy of landslide susceptibility evaluation results; it is necessary to comprehensively consider the internal homogeneity and external heterogeneity of the slope units. Under the optimal slope unit scale, the number of landslides in the highly and extremely highly susceptible areas in the landslide susceptibility map accounted for 82.60% of the total number of landslides, which was consistent with the actual distribution of landslides; this study shows that the method, combining the slope unit and the information value method- random forest (IV-RF) model, for landslide susceptibility evaluation can obtain high accuracy.

1. Introduction

Landslides are one of the most widely distributed and most serious types of geological disasters [1,2]. In recent years, global warming, extreme weather events, frequent earthquakes, and increasing human construction activities have caused landslides to occur more frequently, which resulted in large numbers of casualties and enormous economic loss [3,4,5,6]. Therefore, the quantitative evaluation of landslide susceptibility is of great significance.

Variations in mapping units may lead to different results in the analysis of landslide susceptibility [7,8,9]. Currently, the most widely applied mapping units include grid and slope units. A grid unit is in a regular rectangular shape with data is stored in rectangles, making it convenient for calculation and processing [10,11,12]. A slope unit, an area formed by a watershed line (or ridge line) and a discharge line (a valley line), is a basic topographical unit used for landslides [13,14,15]. A difference in mapping units can lead to different mapping results, thereby affecting the prediction accuracy of the analysis of regional landslide susceptibility [16,17]. The optimal grid unit scale in landslide susceptibility mapping has been extensively studied in recent years; however, the grid unit is not closely related to the actual geological environment where landslides occurred, which may result in a failure to predict landslide susceptibility accurately [18,19,20]. In contrast, a slope unit is closely related to the geological and geomorphological of landslides; however, the complexity in extracting slope units means that few studies have addressed the optimal slope units in regional landslide susceptibility mapping. Hence, it is necessary to carry out landslide susceptibility mapping based on the optimal slope unit.

Landslide susceptibility mapping models can be divided into three categories: empirical, statistical, and machine learning models. The empirical model includes the analytic hierarchy process (AHP) [21] and fuzzy mathematics methods [22]. The disadvantage of empirical methods is that they are subjective. The statistical model includes the certainty factor (CF) method [23], the information value method [24,25], the weight-of-evidence (WOE) method [26,27] and the frequency ratio method [28]. The disadvantage of statistical methods is that they cannot solve the non-linear relationship in the landslide factors. A machine learning model has a strong ability to learn; common methods include a support vector machine (SVM) [10,12], decision tree (DT) [29], random forest (RF) [30,31,32,33], and an artificial neural network (ANN) [11,34]. In recent years, the combination of statistical and machine learning models has been widely used in landslide susceptibility mapping [35,36,37,38]. The combination of different models can not only optimize the results and accuracy of model prediction, but also ensure the stability of model prediction in alpine valley areas with complicated topographical and geological structures along with frequent seismic activities.

In this study, it is assumed that on the basis of selecting the optimal slope unit scale, the information random forest (IV-RF) model has high accuracy in landslide sensitivity assessment in deep valley area. In order to verify this hypothesis, Maoxian County, located in the strong tectonic activity zone on the eastern margin of the Qinghai Tibet Plateau, is taken as the study area. Maoxian County, located in the strong tectonic activity zone on the eastern margin of the Qinghai-Tibet Plateau, is selected as the research area. Considering automatic slope unit partition and the homogeneity of the internal slope aspect and slope gradient. used the r.slopeunits method [17] to extract slope units of different scales and carried out landslide susceptibility mapping based on an information value-random forest (IV-RF) model. The results can provide a theoretical basis for landslide prevention and risk management.

2. Materials and Methods

2.1. Study Area

Located in the middle of the upper reaches of the Minjiang River valley and the southeastern rim of the Qinghai-Tibet Plateau (102°56′26″–104°10′32″E and 31°25′06″–32°15′43″N), Maoxian County covers an area of 3903 km². Most of the area lies within the Qionglai and Minshan mountain ranges, and the southeastern part of the area is in the tail section of the Longmen Mountains. Maoxian County is covered by undulating hills, ridges and peaks, steep valley slopes, narrow river valleys and incised rivers with an average elevation of 3021 m (Figure 1).

Maoxian County has a monsoon climate with annual precipitation of 486.3 mm; rainfall from May to October accounts for over 80% of the annual precipitation. The maximum daily precipitation is 75.2 mm. The precipitation increases with an increase in elevation [39]; this county is located in the Minjiang and Tuojiang river systems and includes the main stream of Minjiang River and its tributary, the Heishui River, as well as the Tumen River of the Tuojiang river system.

The strata in the study area were laid down during the Triassic, Permian, Carboniferous, Devonian, Silurian (Maoxian group), Ordovician and Cambrian geological periods; the lithology is mainly phyllite and slate. The geological structure in the study area is mainly composed of the Jiaochang epsilon-type, the Shidaguan arc-shaped and the Longmenshan Cathaysian tectonic zones. The area also has multiple active faults, including Maoxian-Wenchuan, Shidaguan and Jiudingshan faults; these faults cause geological characteristics such as lithologic fracture and inverted attitude of strata in the area [40]. Located in the middle section of the north-south seismic belt, i.e., the Longmenshan and Songpan seismic belts, earthquakes occur frequently. The strongest earthquake in the area was the 7.5-magnitude Diexi earthquake in 1933.

2.2. Data

2.2.1. Landslide Inventory

As an important basis for landslide susceptibility analysis, the accuracy of landslide inventory data directly affects the reliability of the results [41]. Therefore, the higher the accuracy and completeness of a landslide survey, the better the results of landslide susceptibility mapping. A reliable landslide inventory was obtained using multiple methods in the present study, including historical landslide data, Gaofen-2 satellite images of the study area acquired in 2020 (spatial resolution: 1 m), and Google Earth images acquired in 2011. A total of 297 historical landslides were obtained through a detailedinvestigation of geological hazards. Through the visual interpretation of landslides from Gaofen-2 image in 2020 (spatial resolution of 1 m) and Google Earth image in 2011, combined with field verification, 14 new landslides were obtained through optical remote sensing interpretation. In addition, using Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-InSAR) technology [42,43], 68 landslides with continuous deformation were identified based on 61 ascending and 59 descending orbit scenes acquired by the Sentinel-1A satellite in the study area from January 2019 to December 2020, including 31 landslides identified in the ascending orbit, 41 in the descending orbit and four landslides found in both ascending and descending orbits (Figure 2). Through field verification of 68 landslides, 49 of them are historical landslide points, 14 of them coincide with the landslides interpreted by optical remote sensing, and 5 new landslides are obtained through SBAS-InSAR technology. Using the above methods, a total of 316 landslides were identified (Figure 3). In this study, the IVs of landslide factors were input into the RF model, which then output the significance of various factors and a landslide susceptibility index. According to the principles of the model, the same number of positive (landslide point: 1) and negative samples (non-landslide point: 0) should be chosen and divided into a testing set and a training set while using an appropriate ratio. A total of 316 non-landslide points were randomly chosen from the study area about 1 km away from the landslide points. The landslide and non-landslide points were randomly divided into a testing set (70%) and a training set (30%).

2.2.2. Landslide Causes

The selection of appropriate factors is the key to building a landslide susceptibility mapping model [44]. Common important factors include earthquakes, topography (lithologic and geological structure), hydrological factors, and human activities. According to the actual conditions of the study area, 12 factors, i.e., seismic intensity, engineering strata, geological structure, slope gradient, slope aspect, curvature, topographic wetness index (TWI), relative slope position (RSP), distance from the nearest river, distance from the nearest road, mean annual precipitation, and normalized difference vegetation index (NDVI) were chosen in this study to build the landslide susceptibility mapping model (Figure 4). According to the characteristics of each factor, slope aspect, engineering strata and seismic intensity were discrete variables, and the other factors were continuous variables.

The slope gradient is an important factor that affects the stability of a slope and directly determines the distribution of stress. The slope aspect determines how the local ground surface receives sunlight and redistributes solar radiation, which is directly related to the differences in the local microclimate of an area, thus affecting the vegetation as well as soil moisture and eventually affecting the stability of the slope. Slope curvature determines whether the surface flow flowing through the grid converges or becomes diffuse, thereby affecting the stability of the slope.

As an important topographical attribute, TWI also has a certain influence on landslides [45]; it is defined as the natural logarithm of the specific catchment area (SCA) to the tangent value of surface slope gradient (tanβ):

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

(1)

where, TWI is the topographic wetness index, and SCA is the upstream catchment area of the surface water in a unit length of the contour.

Relative slope position (RSP) is the relative position of the slope to the valley and ridge, which is calculated using Equation (2):

$$\mathrm{RSP}=\left(\frac{\mathrm{z}\left(\mathrm{s}\right)-\mathrm{z}{(\mathrm{s})}_{\mathrm{v}}}{\mathrm{z}{(\mathrm{s})}_{\mathrm{r}}-\mathrm{z}{(\mathrm{s})}_{\mathrm{v}}}\right),$$

(2)

where RSP is the relative slope position (range: 0–1), z(s) is the elevation of the slope, and $\mathrm{z}{(\mathrm{s})}_{\mathrm{v}}$ and $\mathrm{z}{(\mathrm{s})}_{\mathrm{r}}$ are elevations of the valley and the ridge, respectively. The higher the RSP is, the closer to the ridge.

Slope gradient (°), slope aspect, curvature (m⁻¹), TWI and RSP were extracted using SAGA GIS software with a digital elevation model (DEM) from the Advanced Land Observing Satellite (12.5 m).

Data related to two geological factors, lithology and faults, were obtained by vectorization of a 1:200,000 geological map. Faults and engineering petrofabrics are internal causes of landslides. The physical and mechanical parameters of different lithologies vary widely, resulting in different strengths in various slopes. Generally, the rock and earth mass in a fault zone are fragmentized, which will directly affect the stability of a side slope. The study area was divided into seven types of engineering petrofabrics based on the type of lithology present. The distance from the study area to a fault was obtained based on distance analysis.

Earthquakes are an important driving factor for landslides. Areas of high seismic intensity have a high probability of experiencing large-scale landslides [46]. The study area spans four seismic intensities (VII, VIII, IX and X), with VIII being the most frequent. Scouring of river systems and slope cutting for road repairs can lead to slope instability. Graphs of the distances to rivers and roads to the study area were obtained from ArcGIS vector quantization.

Precipitation is another important factor that affects landslides. Rainfall often infiltrates into a slope and increases the weight of the slope, saturating and softening the soil that is subject to sliding, and reducing its shearing strength and sliding resistance [47,48]. The mean annual precipitation data were obtained using Kriging interpolation of data obtained at seven meteorological stations. The NDVI is an index that reflects the surface vegetation coverage, which is calculated as follows:

$$\mathrm{NDVI}=\left(\mathrm{IR}-\mathrm{R}\right)/\left(\mathrm{IR}+\mathrm{R}\right),$$

(3)

where IR and R represent the measured reflectivity of the red and near-infrared spectra, respectively. An NDVI graph of the study area was obtained from Sentinel-2 satellite imagery with a cloud amount of less than 2% in 2020 using Google Earth Engine software.

2.3. Methods

A flow chart of the work in this study was divided into six stages (Figure 5): (1) A landslide inventory was established by integrating data collected by multiple types of remote sensing technology; (2) The K-Means clustering method was used to discretize the causes for continuous landslides; (3) Different scales of slope units were extracted using the r.slopeunits method; (4) Landslide susceptibility mapping was carried out under different slope unit scales based on an information value-random forest (IV-RF) model; (5) The optimal slope unit scale was determined based on the global aspect variance (V), global weighted Moran’s index, AUC and F1 values, and; (6) A landslide susceptibility map using the optimal slope unit was obtained.

2.3.1. Mapping Unit

The selection of mapping units was used as the basis for evaluating regional landslide stability. Considering automatic slope unit partitioning and the homogeneity of the internal slope aspect, the r.slopeunits method was used to extract the slope units [17].

In most landslide stability analyses, researchers assume that the units (i.e., internal slope aspect, slope gradient, mechanical parameters) are homogeneous [15,49,50]. Therefore, the maximum internal homogeneity and the maximum heterogeneity between units should be guaranteed during regional landslide susceptibility mapping based on slope units. The circular variance (c) controls the slope uniformity of slope units in the r.slopeunits method; thus, the global aspect variance (V) and global slope weighted Moran’s index (wMI) were used to represent the global heterogeneity of slope units [51]. The definitions are as follows:

$$V={\sum }_{i=1}^k\frac{S_i}{S}·{\sigma }_i,$$

(4)

where k is the total number of slope units, S is the total area, $S_i$ is the area of the ith slope unit, and ${\sigma }_i$ is the variance of slope aspect of the ith slope unit. The smaller the V value is, the better the slope uniformity.

Moran’s index is often used to measure the heterogeneity between objects. Over- and under-division of slope units have a great impact on Moran’s Index (MI); thus, the global weighted Moran’s index (wMI) was used in this study:

$$wMI=\frac{k^2\sum _{i=1}^k\frac{{\sigma }_i}{S}\sum _{j=1}^k{w}_{ij}\left(y_i-\overline{y}\right)\left(y_j-\overline{y}\right)}{\sum _{i=1}^k{\left(y_i-\overline{y}\right)}^2\left(\sum _{i\ne j}\sum w_{ij}\right)}$$

(5)

where $y_i$ is the mean slope aspect of the ith slope unit, $\overline{y}$ is the mean slope aspect of the entire study area, $w_{ij}$ represents the spatial adjacency relationship of slope units. In this equation, $w_{ij}$ = 1 or = 0 when the units are or not adjacent, respectively. The smaller the value of $w_{ij}$ is, the lower the spatial correlation of the adjacent SU space is, and the stronger the heterogeneity is.

A good partition should lead to low intra-object variance (i.e., low intra-object heterogeneity) and low inter-object Moran’s index (i.e., high inter-object heterogeneity). The global heterogeneity measurement needs to take into account these two indices. In order to balance the intra- and inter-object heterogeneity measurements, the measurements were normalized:

$$\mathrm{GH}=\frac{V-V_{min}}{V_{max}-V_{min}}+\frac{wMI-wM{I}_{min}}{wM{I}_{max}-wM{I}_{min}},$$

(6)

where, $V_{min}$, $V_{max}$, $wM{I}_{max}$ and $wM{I}_{min}$ are the maximum and minimum values of the two indices at different scales. The smaller the GH value is, the better the partitioning result is.

2.3.2. K-Means Clustering

K-Means clustering is an unsupervised learning method that finds the optimal clustering center through multiple iterations such that the sum of the distance between the eigenvalue of each sample and the mean value is minimal [52]:

$$E=min\left\{{\sum }_{i=1}^c{\sum }_{x\in Q_j^{(j)}}{∥g(x)-{\mu }_j^{(i+1)}∥}^2\right\},$$

(7)

where $g\left(x\right)$ is the eigenvalue of sample point x, ${\mu }_j$ is the mean value of the jth characteristic, $Q_j^{\left(i\right)}$ is the collection of the jth type sample points after i iterations, and c is the number of categories.

The natural breakpoint method is often used in the grade classification of factors in traditional landslide susceptibility mapping models. The method is subjective and easily ignores the internal characteristics of units for interlaced topographical factors. Therefore, the K-Means clustering method was used to discretize the continuous factors in slope units to handle factors with interlaced geographical space. As for the discrete factors, the type with the largest proportion of slope units was taken as the slope unit factor. The steps of K-Means clustering are as follows:

(1) Calculate the mean value (MEAN) and standard deviation (STD) for the continuous factors in each slope unit.

(2) The MEAN and STD are input into the K-Means clustering method. The slope aspect, curvature, TWI, RSP, NDVI, mean annual precipitation, distance to the fault, distance to the river, distance to the road, initial number of clustering centers were set to 5, 3, 5, 5, 5, 5, 5, 6 and 5, respectively. Areas with similar topography were classified into the same category.

2.3.3. Information Value (IV) Model

As one of the Bayesian probability models, the information value model [53] presents a quantitative description in the form of a probability, which reflects the contribution of different disaster element intervals to the formation of geological disasters. The expression of geological disaster IV I_Aj→B is as follows:

$${\mathrm{I}}_{\mathrm{Aj}\to \mathrm{B}}=\ln \frac{\mathrm{P}\left(\mathrm{B}/{\mathrm{A}}_{\mathrm{j}}\right)}{\mathrm{P}\left(\mathrm{B}\right)}=\ln \frac{{\mathrm{N}}_{\mathrm{j}}/{\mathrm{S}}_{\mathrm{j}}}{\mathrm{N}/\mathrm{S}}\left(\mathrm{j}=1,2,\cdots \mathrm{n}\right),$$

(8)

where I_Aj→B represents the IV of geological disaster B in interval j of the disaster element A, n is the number of secondary intervals of disaster element A, N_j is the area of a geological disaster or the number of disaster points in interval j of disaster element A, S_j is the distribution area of interval j of disaster element A, N is the total distribution area or total points of regional geological disasters, and S is the total area.

2.3.4. Random Forest (RF)

An RF model is a machine learning algorithm composed of multiple decision trees; it was proposed by Breiman based on the Bagging ensemble learning theory [54]. An RF model includes multiple decision trees trained by the Bagging ensemble learning algorithm. For a sample to be classified, the final result is determined via a voting mechanism of many decision trees; this method has solved the performance problem of decision trees, and it has good parallelism and expansibility for the classification of high-dimensional data and excellent tolerance for noise and abnormalities [55,56,57,58]. Unlike decision trees, the decision-making boundaries of an RF are several segmented functions (Figure 6). The Gini index was used in this study to analyze the importance of input factors in the RF algorithm.

2.3.5. Accuracy Validation

In this study, the area under the ROC curve (AUC) and F1 were jointly used to measure the accuracy of the landslide prediction model [11,23,33]. The F1 is the harmonic mean of precision and recall, which is calculated as follows:

$$\mathrm{F}1=2·\frac{precision·recall}{precision+recall}$$

(9)

$$\mathit{recall} =\frac{TP}{TP+FN}$$

(10)

$$\mathit{precision} =\frac{TP}{TP+FP}$$

(11)

where TP (true positive) means that a positive sample was predicted to be positive, and FP (false positive) means that a negative sample was to be predicted positive.

An FN (false negative) result means that a positive sample was predicted to be negative, and a TN (true negative) result means that a negative sample was predicted to be negative. The greater the values of AUC and F1 are, the higher the accuracy of the model is.

3. Results

3.1. Factor Multicollinearity Analysis

If a correlation exists between the factors in landslide susceptibility mapping, the accuracy and efficiency of the model will be reduced. Therefore, it was necessary to analyze the correlation between the factors and remove any factors that were strongly correlated. Multicollinearity analysis was performed in SPSS software to obtain the values of tolerance (TOL) and variance inflation factors (VIFs) for correlation analysis. When TOL > 0.2, i.e., VIF < 5, the factors are independent of each other; thus, both factors can be used to build the regression model [59,60].

Table 1. Multicollinearity of Factors.

Factor	Tolerance	VIF
Slope	0.40	2.48
Aspect	0.95	1.06
Curvature	0.47	2.14
TWI	0.59	1.69
RSP	0.81	1.24
NDVI	0.84	1.19
Rainfall	0.57	1.76
Distance to faults	0.71	1.40
Distance to rivers	0.90	1.11
Distance to Roads	0.73	1.38
petrofabric	0.48	2.10
Seismic intensit	0.40	2.48

shows that the VIFs of variables of all 12 disaster-inducing factors were less than 5, and the TOLs were greater than 0.2, indicating the independence between these factors. Hence, all 12 factors could be used in the model calculation.

3.2. Selection of the Optimal Slope Unit

In the r.slopeunits method, the threshold (t) of the flow area (FA) and iterative coefficient (r) only control the numerical convergence of the data and have no geomorphological significance, whereas the minimum area (a) and circular variance (c) determine the size and direction of slope units. Because most landslides in the study area were medium-sized, the FA threshold was set to t = 1 × 10⁶ m², the iterative coefficient r = 10, and cleansize = 2 × 10⁴ m². The minimum area a was changed between a = 0.5 × 10⁵ m², 1 × 10⁵ m², 1.5 × 10⁵ m², 2 × 10⁵ m², 2.5 × 10⁵ m² and 3 × 10⁵ m², the circular variance c had a uniform distribution with an increment of 0.05 in the range of 0.1 ≤ c ≤ 0.4. A total of 42 groups of slope units at varying scales were obtained through different combinations of a and c. Figure 7 shows three slope units.

The number of slope units was used in this study to reflect different partition methods in the study area.

Table 2. The number of Slope Units at Different Scales.

Number of Slope Units		a (105 m2)
		0.5	1	1.5	2	2.5	3
c	0.1	23,498	16,906	13,444	11,132	9537	8425
	0.15	11,466	9499	8208	7217	6439	5882
	0.2	6181	5624	5132	4770	4476	4237
	0.25	3840	3635	3481	3362	3428	3154
	0.3	2734	2810	2734	2684	2635	2547
	0.35	2450	2400	2358	2356	2308	2242
	0.4	2268	2226	2206	2190	2182	2182

shows that the maximum partition number of slope units is 23,498 and the minimum number is 2242. Figure 8 shows the influence of different a and c combinations on the number of slope units. According to Figure 8a, when c ≥ 0.3, the number of slope units converges almost to the same value, suggesting that when c ≥ 0.3, the number of slope units does not change with c and a. From Figure 8b, (1) when c = 0.1, 0.15 or 0.2, the number of slope units changes as the value of a changes; the change is the smallest when c = 0.2, and the number of slope units starts to converge when a ≥ 1.5 × 10⁵ m². (2) When a = 3 × 10⁵ m², the number of slope units converges with different c values.

Thus, the following 23 groups are of value in landslide susceptibility mapping in the study area, (1) c = 0.1, 0.15, and 0.2, combined with different a values, (2) c = 0.25 and a = 0.5 × 10⁵ m², 1 × 10⁵ m², 1.5 × 10⁵ m² and (3) c = 0.3 and a = 0.5 × 10⁵ m², 1 × 10⁵ m².

Susceptibility mapping was carried out for the 23 groups of slope units based on the IV-RF model, and the AUC and F1 values of the model were calculated (Figure 9a,b). Equation (6) was used to calculate the GH value of internal homogeneity and external heterogeneity of slope units (Figure 9c). In this paper, the maximum value of AUC was 0.905, and the minimum value was 0.857, with an average value of 0.890. Based on the conclusion of Hosmer and Lemeshow that AUC > 0.7 indicates high prediction accuracy, the proposed model in this study shows excellent spatial prediction ability.

Figure 9a,b show that there is little difference between AUC and F1, which had a similar trend. When c = 0.1, AUC and F1 reached the maximum values of 0.905 and 0.908 simultaneously at a = 3 × 10⁵ m². When c = 0.15, AUC and F1 reach the maximum values of 0.905 and 0.908 simultaneously at a = 2 × 10⁵ m². When c = 0.2, AUC and F1 reach the maximum values of 0.905 and 0.906 simultaneously at a = 1.5 × 10⁵ m² and a = 2 × 10⁵ m², respectively. When c = 0.25, AUC and F1 reach the maximum values of 0.894 and 0.899 simultaneously at a = 1 × 10⁵ m². When c = 0.3, AUC and F1 reach the maximum values of 0.884 and 0.892 simultaneously at a = 1 × 10⁵ m². Figure 9c shows that when c is a constant, the GH value decreases with the increase of a, indicating good partition results. The minimum value of GH was 0.872 when c = 0.1 and a = 3 × 10⁵ m². In summary, the combination of c = 0.1 and a = 3 × 10⁵ m² leads to the minimum GH and the maximum values of AUC and F1, i.e., optimal internal homogeneity and external heterogeneity and the best accuracy of landslide susceptibility mapping.

3.3. Landslide Susceptibility Map

The landslide susceptibility indices of slope units extracted by the combination of c = 0.1 and a = 3 × 10⁵ m² were divided into five grades: extremely high, high, medium, low and extremely low landslide susceptibility (Figure 10). Most areas of Maoxian County are located in the extremely low and low landslide susceptibility areas. Extremely high and high landslide susceptibility areas were mainly distributed in river valleys and both sides of major roads (e.g., Minjiang, Heishui and Tumen rivers, and road G213). The topography in these areas is fragmented and cracked, and the rocks are mostly phyllite and metamorphic rocks with an interbedded soft and hard rock structure, which can easily form a weak structural surface that is typically affected by rainfall and groundwater activities. The slope vegetation coverage in these areas was less than 10%; moreover, human economic activities and neotectonics are commonly present here. All these adverse factors contribute to the occurrence of landslides.

Table 3. Statistics of Different Susceptibility Areas.

Susceptibility Class	Number of Slope Units	Slope Units (%)	Area (km2)	Area (%)	Landslide (Points)	Landslides (%)
Very Low	1076	12.77%	481.78	12.10%	1	0.32%
Low	4815	57.15%	2338.20	58.74%	15	4.75%
Medium	834	9.90%	387.83	9.74%	39	12.34%
High	1108	13.15%	516.56	12.98%	108	34.18%
Very High	592	7.03%	256.07	6.43%	153	48.42%

and Figure 11 shows the number of slope units in each area with a different level of susceptibility and the number of previous landslides that occurred historically in each area. The proportion of previous landslides increased with an increase in the susceptibility grade. Areas of extremely low and low susceptibility accounted for 70.84% of the total area, and the number of previous landslides in these areas only accounted for 5.07% of the total number. In addition, areas of extremely high and high susceptibility only accounted for 19.41% of the total area, yet previous landslides here accounted for 82.60% of the total number. The results are in line with the actual situation.

4. Discussion

4.1. Reliability of the IV-RF Model

Multi-phase high-resolution remote sensing images, interferometric synthetic aperture radar (InSAR), historical landslide and field surveys were integrated to obtain 316 landslide samples with high accuracy; this study combines optical remote sensing with radar remote sensing to identify the potential landslide sites. Twelve factors were included, namely, seismic intensity, engineering strata, geological structure, slope gradient, slope aspect, curvature, TWI, RSP, distance from the nearest river, distance from the nearest road, mean annual precipitation and NDVI to carry out landslide susceptibility mapping based on the IV-RF. The IV model has been frequently used in previous studies [24,25]; however, it ignores the complex interactive relationship between the occurrence of landslides and factors influencing landslides. The RF model can capture the complex non-linear relationship between the factors and occurrence of landslides and has high prediction accuracy. In a study by Dou et al. [33], two groups of training samples were used to compare the landslide susceptibility mapping results of DT and RF models in the volcanic island area of Izu-Oshima in Japan [33]. The results showed that the RF model was superior to the DT model in different training sets; moreover, using logistic regression (LR), RF and support vector machine (SVM), Chang et al. [61] carried out landslide susceptibility mapping and found that RF has the best performance and is effective for DEMs with different resolutions. Numerous studies have shown that the mixed model has some advantages such as high accuracy and efficiency. Hence, the IV-RF model was adopted in this study to generate landslide susceptibility maps. In this study, the maximum value of AUC was 0.905, the minimum was 0.857, and the mean value was 0.890. According to Hosmer and Lemeshow, when AUC is greater than 0.7, the model has good prediction accuracy; thus, the proposed model in this study has excellent spatial prediction ability. The method based on slope unit scale and IV-RF model adopted in this study is suitable for landslide Susceptibility evaluation in areas with large topographic relief. In the future, we will strengthen the research on error conduction, numerical simulation of landslide and the correlation between landslide sensitivity evaluation results and landslide characteristics [62,63].

4.2. Optimal Mapping Units

The mapping unit is the basis of landslide susceptibility mapping. Different scales of mapping units will lead to different mapping results. Cama et al. [19] used forward stepwise binary logic regression (BLR) to model the susceptibility of an area to landslides with 2, 4, 16 and 32 m grid cells, the results of which showed that the 8 and 16 m groups have higher accuracy. The slope unit is the basic topographic unit of landslides, and currently few studies have analyzed the optimal slope unit scale in regional landslide susceptibility mapping. Positive and negative DEM hydrological analysis is the most common method used to extract slope units, yet this method is inefficient and is not conducive to the extraction of large-scale slope units, while the internal homogeneity of slope units is not considered. In comparison, the r.slopeunits method proposed by Alvioli et al. [17] can ensure the homogeneity of the slope aspect of the extracted slope units and controls the size and direction of slope units by the minimum area threshold a and circular variance c. The r.slopeunits method was used in this study to obtain 23 groups of slope units with different scales using different parameter combinations (a, c). The internal homogeneity of slope units and heterogeneity between different units were analyzed based on global aspect variance (V) and global weighted Moran index (wMI). The AUC and F1 were calculated to analyze the accuracy of the 23 groups of slope units in the IV-RF landslide susceptibility mapping model. Last, the optimal slope unit was decided based on V, wMI, ROC and F1. Specifically, when c = 0.1 and a = 3 × 10⁵ m², the internal homogeneity and external heterogeneity of slope units were optimal, resulting in the best accuracy of landslide susceptibility mapping. The predicted spatial distribution of landslides in extremely high and high susceptibility areas was in line with the actual conditions.

The current study has some limitations. First, K-Means clustering was used to discretize continuous factors, so the clustering results may be ineffective if the initial clustering center in the K-Mean algorithm was not good. In the next step, we will compare the results of different clustering methods, such as hierarchical clustering, self-organizing-map (SOM) clustering, and Gaussian mixture model [64]. Second, non-landslide points in this study were random samples in areas about 1 km away from landslide points (buffer zone); however, there might be potential landslide points in the selected non-landslide points. Zhu et al. [65] has explored different methods such as a one-class support vector machine (one-class SVM), kernel density estimation, an artificial neural network (ANN) and a two-class support vector machine (two-class SVM) to generate non-landslide units. In a study by Huang et al. [66], a self-organizing-map (SOM) network was applied to generate the initial landslide susceptibility map and non-landslide grid cells were chosen from the extremely low susceptibility areas. A target spatial externalized sampling method was proposed by Xiao et al. [67] to generate non-landslide data directly using the eigenspace of existing landslides. We believe there will be more and more effective methods created that are designed to improve the quality of non-landslide points.

5. Conclusions

The objective of the present study was to carry out landslide susceptibility mapping in a deep valley area based on the optimal slope unit. The following conclusions are drawn.

(1) The mixed IV-RF model has excellent spatial prediction ability and accuracy in landslide susceptibility mapping.

(2) When c = 0.1 and a = 3 × 10⁵ m², 8425 slope units were extracted by the r.slopeunits method, the internal homogeneity and external heterogeneity were optimal, indicating the model had the best level of accuracy. In other words, the optimal slope unit of a landslide susceptibility model is related to the category and size of landslides, rather than the size of slope units.

(3) The spatial distribution of high and extremely high susceptibility areas is consistent with that of existing landslides under the optimal slope units.