Coupling Data- and Knowledge-Driven Methods for Landslide Susceptibility Mapping in Human-Modified Environments: A Case Study from Wanzhou County, Three Gorges Reservoir Area, China

Abstract

Landslide susceptibility mapping (LSM) can provide valuable information for local governments in landslide prevention and mitigation. Despite significant improvements in the predictive performance of LSM, it remains a challenge to be carried out in areas with limited availability of data. For example, in the early stage of road construction, landslide inventory data can be particularly scarce, while there is a high need to have a susceptibility map. This study aims to set up a novel procedure for coupling the knowledge-driven and data-driven models for LSM in an area with limited landslide inventory data. In particular, we propose a two-step approach. The first step consists of applying four data-driven models (logistic regression, decision tree, support vector machines, and random forest (RF)) to derive a regional susceptibility map. In the second step, the application of a heuristic model (analytic hierarchy process, AHP) is proposed to calculate a local susceptibility map for the areas with incomplete landslide inventories. The final landslide susceptibility map is obtained by merging the most accurate regional map (RF) with the local map. We apply this novel procedure to a landslide-prone region with developed road construction (National Highway G69) in Wanzhou district, where landslide inventory is difficult to update due to timely recovery from landslide-induced road damage. Results show that the proposed methodology allows identifying new landslide-prone areas, and improving LSM predictive performance, as demonstrated by the fact that two new landslides developed along G69 were perfectly classified in the highly susceptible areas. The results show that implementing the landslide susceptibility assessment with different geographical settings and combining them into best-sensitivity partitions is more accurate than focusing on creating new models or hybrid models.

1. Introduction

Landslides pose great threats to human lives and economic development around the world [1]. The mountainous areas of southwestern China suffer a lot from landslide hazards. For instance, the Shuicheng landslide-debris occurred in Guizhou Province, China, resulting in 21 buildings being buried, 9 missing people, and 42 deaths in 2019 [2]. Another catastrophic landslide in Jiweishan Mountain, Chongqing province, China, caused 79 deaths in 2009 [3].

Landslide susceptibility is a spatial probability of an area to the future occurrence of landslides [4,5]. It predicts “where” landslides are likely to occur, which is the foundation of landslide early warning and risk assessment. A variety of methods have been proposed to ascertain landslide susceptibility. Essentially, most of the approaches and methods are based upon the assumption that slope instability is more likely to occur under the conditions that already led to landslides in the present and past [1].

Characterising the spatial pattern of landslide occurrence and improving the accuracy of landslide susceptibility mapping (LSM) are difficult missions. To date, the landslide susceptibility models can be divided into four categories: inventory-based methods, physically based models, heuristic models, and data-driven models [6]. When the landslide inventory data are limited, we can use heuristic models based on expert knowledge of landslide deformation mechanisms [7,8]. With the fast development of geographic information systems (GISs) and soft computing techniques, many advanced data-driven models have been proposed to conduct LSM [9]. The data-driven models can be mainly divided into statistic models and machine learning (ML) models [10]. The statistical model includes an information value model [11] and logistic regression (LR) [12], while the ML model includes artificial neural networks [13,14], decision tree (DT) [15,16], support vector machine (SVM) [17,18], and random forest (RF) [19], etc. The present scientific literature reports that ML models perform very well in many research areas due to their mathematical flexibility [1,20]. Although nonlinear correlations between landslide and related conditioning factors can be efficiently fitted and predicted by data-driven models, the assessment accuracy is very reliant on the completeness of the landslide inventory.

At present, a reliable landslide inventory can be obtained by many advanced remote sensing techniques. One of these techniques is Interferometric Synthetic Aperture Radar, etc. [21,22]. Nevertheless, due to some inherent limitations, such as dense vegetation and climate and land change, in many cases, landslides are difficult to spot. Missing landslide inventory data affected the accuracy of LSM obtained by data-driven models. In hilly areas, slope-cutting is inevitable during the road-construction period. Slope-cutting often alters the original stress balance of the natural slope [23,24], which in turn leads to the fracture of the rock mass structure in the unloading area, the exposure of weak layers, and so on. Therefore, the land modifications caused by road construction often increase the frequency of landslides occurrence [25,26]. More often than not, slope-cutting-induced landslides do not occur immediately but instead are due to progressive failure progress due to the long-term action of dynamic loads of the road, rainfall, etc.

In mountainous areas, in order to meet the functional requirements of the road, most of the slope-cutting induced sliding mass along the road are recovered timely. This makes it challenging to build an accurate and significant landslide inventory that can be effectively used for LSM and opens the need to propose novel landslide susceptibility mapping procedures that allow operating in situations of incomplete landslide inventories.

Given the above issues, in this study, we propose a novel procedure to assess landslide susceptibility in areas with incomplete landslide inventories. The procedure couples data-driven and knowledge-driven approaches. A landslide-prone region with slope-cutting caused by road construction in Wanzhou District was selected as a case study. Four data-driven models (LR, DT, SVM, RF) and heuristic models (analytic hierarchy process) were combined to predict landslide susceptibility of the whole region and the area along National Highway G69, respectively. We produced the final landslide susceptibility map by merging the most accurate regional landslide susceptibility map with the road-affected landslide susceptibility map. Where susceptibility zonation is prepared for land planning, early warning or other practical applications, the authors hoped that the proposed two-step procedure could complete a more accurate susceptibility map to reduce the economic damage and human losses coursed by landslide hazards.

2. Study Area

The study area is located in Wanzhou County, the southwestern of China (108°17′47″ E~108°21′35″ E in longitude and 30°32′02″ N~30°34′36″ N in latitude) (Figure 1). It is composed of hills and mountains (up to 770 m of the elevation) with an area of 12.9 km². Wanzhou belongs to the subtropical monsoon area, with an average annual rainfall of 1191 mm. The roads in the study area are well developed, and the newly built important traffic line, G69 national highway, runs through the whole area (Figure 1c). After field investigations, a 50-m buffer zone of G69 was selected as the target area for the road-affected LSM (this area is herein referred to as “G69”).

In this study area, as the distance from the Yangtze River increases, the exposed bedrock becomes older. Bedrock mostly formed during the Mesozoic Triassic and Jurassic from 137 to 250 Ma. The rock stratum tilts west in general, and the dip angle of the rock stratum gradually decreases to nearly horizontal from the mountain to the Yangtze River (Figure 1d). In addition, the rock strata dip is nearly perpendicular to the flow direction of the Yangtze River, with an average dip angle of 15° [27]. The fluctuation of the reservoir water level, rainfall and hill-cutting are the crucial triggering factors of landslide occurrence. For example, the Jinjinzi landslide was triggered by heavy rainfall and the rapid drawdown of reservoir water level (Figure 2a), while the Tongzilin landslide was induced by slope-cutting during construction activities (Figure 2c). The rock strata of the Tongzilin landslide have the same attitude as the slope surface.

3. Methodology

3.1. Modelling Algorithms

3.1.1. Information Value Model

The information value (IV) model can usefully calculate the importance of landslide causal factors, and the larger IV, the more important for landslide occurrence [11]. Assuming that landslide occurrence (y) is affected by various factors ($x_1,x_2\dots x_n$), the formula for IV can be written as:

$$I\left(y,x_1,x_2\dots x_n\right)={\displaystyle \sum _{i=1}^n{I}_i}={\displaystyle \sum _{i=1}^n{\log }_2\frac{S_0^i/S^i}{A_0/A}},$$

(1)

where $I\left(y,x_1,x_2\dots x_n\right)$ is the IV under the landslide causal factor, I_i is the IV of factor $x_1$, A is the total area, A₀ is the area under the factor $x_1$, S is the total area of landslides, and S₀ is the area of landslides under the factor $x_1$.

3.1.2. Logistic Regression

Logistic regression (LR) is a widely employed regression algorithm [13]. High calculation speed is among its advantages. In the application of LSM, the learning goal of the LR model is to quantitatively describe the relationship between the occurrence probability of a landslide and its causative factors. The calculation formula of LR can be expressed as:

$$\mathrm{y}(\mathrm{x})=\frac{1}{1+e^{-{\theta }^{T}x}},$$

(2)

where the independent variable x refers to the factors that may cause landslides, such as altitude, slope, and aspect. The dependent variable y is landslide susceptibility index returned by LR, with a range of (0,1) [14].

3.1.3. Decision Tree

The decision tree (DT) model is widely used for prediction tasks and consists of tree growth, tree pruning, and boosting steps [15,16]. It works by splitting the sample based on the field that provides the maximum information gain. DT model represents the conditional probability distribution of the landslide occurrence or non-landslides under the conditions of given factors. The main advantage is that its classification is fast and easy to understand.

3.1.4. Random Forest

Random forest (RF) is a representative Bagging ensemble algorithm based on the decision tree as the evaluator, following the principle of average or minority to the majority to determine the ensemble result [28,29]. The higher the accuracy of a single decision tree, the higher the accuracy of the RF. Many studies show that the RF model provides a more effective ability to avoid over-fitting problems and generally achieves better generalisation performance than other ML techniques [1]. The performance of the RF model is mainly controlled by the number of trees in the forest N, the number of random factors M employed in each random tree, and the parameters of the maximum depth and a minimum number of samples. This study applies the grid optimisation function to select the optimal parameter combination.

3.1.5. Support Vector Machine

Support vector machine (SVM) is one of the algorithms that attracts most people’s attention in the ML [30]. The learning goal of SVM model is to find the separation hyperplane in the feature space, which can classify the training data set in the feature space into different classes. Radial Basis Function was selected to establish the SVM model in LSM.

3.1.6. Analytic Hierarchy Process

The analytic hierarchy process (AHP) is a semi-quantitative decision-making method of subjective weight assignment [31]. It adopts the simplification strategy of complex problems layer by layer and determines the weight of complex index systems through pairwise comparison. AHP model includes four steps: (a) decomposing the complex problems to be solved into components at different levels; (b) arranging causal factors hierarchically and establishing a hierarchical relationship model among indicators; (c) establishing a judgment matrix based on the important relationship between each index and verifying its rationality; (d) calculating the maximum eigenvalue of the judgment matrix, and giving the index weight in the corresponding level.

3.2. Proposed Framework for LSM

In this study, a procedure coupling heuristic and data-driven models was used to obtain high-quality landslide susceptibility map in complex geographic setting. This procedure includes the following six stages (Figure 3): (a) collect data required for LSM, including historical landslides data and a suite of landslide causal factors; (b) select the landslide causal factors for regional modelling, and construct a regional geospatial database with collected historical landslides; (c) construct a geospatial database of road-affected area, which including the causal factors of landslides developed under the slope-cutting conditions; (d) producing the preliminary landslide susceptibility map based on data-driven model; (e) generating the landslide susceptibility map of the road affected area (G69) by heuristic model based on expert knowledge; (f) producing the final susceptibility map by coupling the two landslide susceptibility maps of steps d and e.

4. Data Preparation

4.1. Regional Geospatial Database

4.1.1. Landslide Inventory

An accurate landslide inventory map can boost the understanding of the relationship between historical landslides distribution and landslide causal factors. Through field surveys and historical data collection, we built an inventory of 23 recorded landslide events with the attribution of location, type, boundary, material, and deformation characteristics. The landslide sizes range from 9987 m² to over 482,976 m², which covers 9.7% of the total area.

4.1.2. Landslide Conditioning Factors

Based on the field survey and previous studies [32,33,34], 12 factors were initially selected for LSM in this paper, namely altitude, slope, aspect, curvature, plan curvature, profile curvature, topographic wetness index (TWI) [35], terrain roughness index (TRI) [36], lithology, bedding structure [37], distance to rivers, and distance to roads (Figure 4). Altitude, slope, aspect, curvature, TWI, and TRI were obtained by digital elevation model (DEM). Distance to rivers and distance to roads were calculated from the Euclidean distance analysis tool. Lithology and bedding structures were provided by the modified geological map. All the processing of factors in this research were conducted in ArcGIS 10.6 software.

Before modelling, landslide conditioning factors should be classified and normalised [19]. In this study, the method proposed by Yu et al. [37] was selected to discretise the continuous factors; after that, all factors were normalised between 0.01 and 0.99 based on the IV model. Finally, all the factors were divided into 5~8 classes (

Table 1. Spatial relationships between landslide causal factors and landslides using IV model.

Factors	Class	PL	PTD	PLTL	PDTD	IV	Normalised Class
Altitude (m)	<250	7427	24,137	0.64	0.20	1.15	0.99
	250~425	3907	44,856	0.33	0.37	−0.11	0.65
	425~500	194	22,932	0.02	0.19	−2.44	0.01
	500~550	164	9773	0.01	0.08	−1.76	0.19
	>550	0	18,344	0.00	0.15	∞	0.01
Slope (°)	<6	1992	7874	0.17	0.07	0.95	0.99
	6~15	5706	38,648	0.49	0.32	0.42	0.83
	15~24	2961	38,401	0.25	0.32	−0.23	0.63
	24~36	938	24,642	0.08	0.21	−0.94	0.41
	36~48	95	9145	0.01	0.08	−2.24	0.01
	48~69	0	1332	0.00	0.01	∞	0.01
Aspect (°)	0~45	951	17,278	0.08	0.14	−0.57	0.67
	45~90	212	9860	0.02	0.08	−1.51	0.39
	90~135	22	3713	0.00	0.03	−2.80	0.01
	135~180	331	8321	0.03	0.07	−0.89	0.57
	180~225	1495	17,495	0.13	0.15	−0.13	0.8
	225~270	2145	18,064	0.18	0.15	0.19	0.9
	270~315	4793	28,871	0.41	0.24	0.53	0.99
	315~360	1743	16,440	0.15	0.14	0.08	0.86
Curvature	−15 to −125	107	2973	0.01	0.02	−1.00	0.01
	0 to −15	2029	24,461	0.17	0.20	−0.16	0.71
	0~5	8299	71,077	0.71	0.59	0.18	0.99
	5~66	1257	21,531	0.11	0.18	−0.51	0.41
Plan curvature	−58 to −35	0	35	0.00	0.00	∞	0.01
	−35 to −5	262	5881	0.02	0.05	−0.78	0.01
	−5 to 0	6159	55,862	0.53	0.47	0.12	0.99
	0~25	5271	58,249	0.45	0.49	−0.07	0.78
	25~30	0	15	0.00	0.00	∞	0.01
Profile curvature	−53 to −40	0	16	0.00	0.00	∞	0.01
	−2 to −40	2261	32,064	0.19	0.27	−0.32	0.07
	−2 to 5	8588	75,401	0.73	0.63	0.15	0.99
	5~30	843	12,423	0.07	0.10	−0.36	0.01
	30~76	0	138	0.00	0.00	∞	0.01
Topographic roughness index	0~1.1	10,748	86,559	0.92	0.72	0.24	0.99
	1.1~1.2	751	19,303	0.06	0.16	−0.92	0.67
	1.2~1.3	170	8041	0.01	0.07	−1.53	0.49
	>1.3	23	6139	0.00	0.05	−3.26	0.01
Topographic wetness index	−3.4 to −2	714	13,159	0.06	0.11	−0.58	0.01
	−2 to 1	1308	7989	0.11	0.07	0.52	0.99
	1~5	3760	53,213	0.32	0.44	−0.32	0.24
	5~15	5910	45,661	0.51	0.38	0.28	0.78
	15~17	0	20	0.00	0.00	∞	0.01
Distance to rivers (m)	0~100	1119	16,968	0.10	0.14	−0.39	0.29
	100~250	5542	44,602	0.47	0.37	0.24	0.99
	250~500	4500	48,041	0.38	0.40	−0.04	0.68
	>500	531	10,431	0.05	0.09	−0.65	0.01
Distance to roads (m)	0~100	3919	38,865	0.34	0.32	0.03	0.95
	100~250	5690	54,083	0.49	0.45	0.08	0.99
	250~450	1859	21,821	0.16	0.18	−0.13	0.76
	>450	224	5273	0.02	0.04	−0.83	0.01
Lithology	L1	0	6666	0.00	0.06	∞	0.01
	L2	2722	27,915	0.23	0.23	0.01	0.46
	L3	6020	50,556	0.51	0.42	0.20	0.60
	L4	750	15,173	0.06	0.13	−0.67	0.01
	L5	545	10,365	0.05	0.09	−0.61	0.04
	L6	413	1908	0.04	0.02	0.80	0.99
	L7	1242	7500	0.11	0.06	0.53	0.82
Bedding structure	BS1	705	6852	0.06	0.06	0.06	0.72
	BS2	5111	19,535	0.44	0.16	0.99	0.99
	BS3	1883	24,231	0.16	0.20	−0.23	0.64
	BS4	1747	24,366	0.15	0.20	−0.31	0.61
	BS5	1930	29,408	0.17	0.24	−0.39	0.58
	BS6	268	10,457	0.02	0.09	−1.33	0.30
	BS7	48	5193	0.00	0.04	−2.35	0.01

Note: PL = pixels in landslide, PTD = pixels in total domain, PLTL = proportion of landslide in total landslide, PDTD = proportion of domain in total domain.

). The lithology was divided into seven classes: limestone and dolomite rock (L1), mudstone (L2), sandstone (L3), sandstone with mudstone (L4), shale with dolostone (L5), mudstone with silty mudstone (L6), and sand–mud interbed (L7). According to the China Geological Survey, bedding structure can be classified into eight categories: horizontal strata slope (BS1), over-dip slope (BS2), under-dip slope (BS3), dip-oblique slope (BS4), transverse slope (BS5), anaclinal oblique slope (BS6), and anaclinal slope (BS7) [37]. All indicator layers were converted into a 10 m spatial resolution raster format as independent variables in preparation for spatial modelling.

4.2. Geospatial Database of G69

At the G69 region, there are no data on historical landslides induced by slope cutting. Based on the knowledge of Geotechnics and Rock Mechanics, we selected six assessment factors for modelling, including slope, bedding structure, rock association, rock mass structure, cutting height, and overburden thickness. Slope and bedding structure are the same factor of the regional database, while slope refers to the natural slope angle before cutting. Cutting height refers to the height of the slope surface formed after the excavation of the slope foot when the road is constructed, which is a decisive factor in controlling the stability of the slope. Rock mass structure refers to the spatial spread of structural surfaces in a rock mass and their combination characteristic [23]. In addition, rock association is one of the factors which reflect the structural characteristics of slopes [38].

The cutting height data along the road were obtained by linear interpolation of the survey points. They were divided into 0 m, 0~5 m, 5~10 m, 10~20 m, and >20 m, respectively (Figure 5). Similarly, the collected data were interpolated by kriging to obtain the overburden thickness, which was divided into four categories: <0.5 m, 0.5~1 m, 1~1.5 m, and >1.5 m, respectively. The rock mass was classified into a blocky structure, layer structure, thin layer structure, interlayer structure, fragmented structure, and scattered structure [23]. The rock association is divided into nine categories by combining the distribution characteristics of soft layer location (Figure 6).

5. Results

5.1. Regional Landslide Susceptibility Modelling

5.1.1. Factor Selection

Multicollinearity Analysis

Numerous studies have contended that eliminating collinearity factors can improve the accuracy of LSM [39,40]. In this paper, the variance inflation factors (VIF) and the tolerances were used to quantify the correlation between landslide causal factors. The result is shown in

Table 2. Multicollinearity of the causal factors.

Factor	Original Factor		New Factor
	Tolerances	VIF	Tolerances	VIF
Altitude	0.797	1.254	0.629	1.590
Slope	0.234	4.271	0.632	1.583
Aspect	0.909	1.100	0.824	1.214
Plan curvature	0.873	1.145	0.731	1.369
Profile curvature	0.643	1.555	0.647	1.544
TRI	0.252	4.969	/	/
TWI	0.807	1.240	0.899	1.112
Lithology	0.918	1.090	0.951	1.051
Bedding structure	0.748	1.337	0.639	1.565
Distance to rivers	0.975	1.026	0.844	1.185
Distance to road	0.944	1.059	0.879	1.138

. The statistical theory states that when VIF is ≥5 or tolerance is ≤0.2, the factors have collinearity. The VIF and tolerance of TRI are 0.252 and 4.969, and the VIF and tolerance of slope are 0.234 and 4.271. That means there is collinearity between TRI and slope. After removing TRI, the minimum tolerance is 0.629, and the maximum VIF is 1.590 (Table 2). There is no collinearity among the new landslide causal factors.

Importance Analysis

Importance analysis is usually used to eliminate redundant landslide causal factors during the landslide susceptibility modelling process [41,42]. Compared to two other methods (filter and wrapper methods), the embedded feature selection (EFS) method can improve computational efficiency and prediction performance of data-driven models by allowing the learner to obtain its own optimal feature subset [43,44].

In this study, EFS was used to estimate the relative importance of factors and select the best subset of the factors for LSM. The RF model was selected as the base learner, and the area under the curve (AUC) calculated by receiver operating characteristic curve (ROC) was used as the performance metric [1]. The feature importance obtained by RF reflects the prediction ability of conditioning factors on landslide occurrence (Figure 7). The closer it is to 1, the greater the role of this factor. As shown in

Table 3. The prediction accuracy after eliminating factors.

Model	Eliminating Less Important Factors	Accuracy
Model 1	Without eliminating any factor	0.824
Model 2	Profile curvature	0.914
Model 3	Plan curvature	0.917
Model 4	Profile curvature, plan curvature	0.907
Model 5	Curvature	0.904
Model 6	TWI, profile curvature, plan curvature, curvature	0.901
Model 7	TWI, slope, profile curvature, plan curvature, curvature	0.901

, when eliminating plan curvature, the AUC of susceptibility modelling is the highest at 0.917. Hence, ten important factors were utilised for LSM at last.

5.1.2. Modelling

Before modelling, the spatial database should be divided into a training dataset and a test dataset. In this study, 50% of landslide grid cells were randomly selected for model training, while the remaining 50% were applied to verify the model performance. The same number of non-landslide pixels was randomly selected as training samples. Furthermore, to illustrate the effectiveness of four data-driven models (LR, DT, RF, SVM) in the study area, the AHP model was treated as a control. Finally, five models of AHP, LR, DT, RF, and SVM were used for LSM. The modelling process was completed in Python 3.7.

Parameters have a significant influence on the model performance. In this study, the parameters for the AHP model was shown in

Table 4. Description of 10 causal factors and grade scores of variables in AHP model.

Factors	Weight	Categories (Normalised Values)
Altitude (m)	0.108	<250(0.35), 250~425(0.3), 425~500(0.17), 500~550(0.15), >550(0.03)
Slope (°)	0.161	<6(0.24), 6~15(0.32), 15~24(0.21), 24~36(0.1), 36~48(0.05), 48~69(0.05)
Aspect (°)	0.096	0~45(0.07), 45~90(0.03), 90~135(0.07), 135~180(0.11), 180~225(0.11), 225~270(0.17), 270~315(0.27), 315~360(0.17)
Curvature	0.021	−125 to −15(0.21), −15 to 0(0.27), 0~5(0.31), 5~66(0.21)
Profile curvature	0.014	−53 to −40(0.11), −40 to −2(0.17), −2 to 5(0.4), 5~30(0.21), 30~76(0.11)
TWI	0.021	−3.4 to −2(0.07), −2 to 1(0.37), 1~5(0.25), 5~15(0.2), 15~17(0.11)
Distance to rivers (m)	0.156	0~100(0.45), 100~250(0.31), 250~500(0.15), >500(0.09)
Distance to roads (m)	0.041	0~100(0.35), 100~250(0.25), 250~450(0.22), >450(0.18)
Lithology	0.161	L1(0.05), L2(0.11), L3(0.11), L4(0.05), L5(0.11), L6(0.37), L7(0.3)
Bedding structure	0.221	BS1(0.2), BS2(0.32), BS3(0.15), BS4(0.11), BS5(0.07), BS6(0.11), BS7(0.05)

, and the parameters for LR, DT, RF, and SVM model were determined using the grid search module provided by the ski-learn library. The optimised parameters of the four models are shown in

Table 5. The parameters of four models.

Model	Parameters	Notes
LR	C = 0.3	C is the reciprocal of regularisation strength.
DT	max_depth = 4	Max_depth is the maximum depth of the tree.
RF	n_estimators = 73, max_depth = 30	n_estimators is the number of trees in random forests, max_depth is the maximum depth of the tree.
SVM	C = 5, gamma = 0.45	C is the penalty term coefficient of the relaxation coefficient, and gamma is the coefficient of the kernel function.

.

The landslide susceptibility index obtained from used models was divided into five levels according to the proportion: very high (10%), high (20%), moderate (20%), low (20%), and very low (30%). The susceptibility maps are shown in Figure 8.

5.1.3. Accuracy Comparison

After modelling, we used the ROC curve to compare the accuracy of these five models. As shown in Figure 9, the AUC of the RF model was the highest at 0.922; the others were 0.916 (SVM), 0.873 (DT), 0. 805 (LR) and 0. 722 (AHP). Therefore, based on the ROC curve, the RF model presents the best generalisation and accuracy in the regional LSM. In the following, “the original landslide susceptibility map” refers to the map produced by the RF model.

5.2. Landslide Susceptibility Modelling of G69

The AHP model was used to map the landslide susceptibility of G69 without landslide inventory. When using the AHP model for LSM, its decisions are taken using weights through pairwise relative comparisons. The final weights and quantitative values of all indicators are shown in

Table 6. Causal factors and grade scores of variables in AHP model.

Factors	Weight	Categories	Normalised Values	Factors	Weight	Categories	Normalised Values
Slope/(°)	0.1194	<10	0.130	Overburden thickness/(m)	0.081	<0.5	0.123
		10~20	0.193			0.5~1	0.198
		20~30	0.252			1~1.5	0.316
		>30	0.425			>1.5	0.363
Rock association	0.2287	RA1	0.155	Rock mass structure	0.148	Massive structure	0.068
		RA2	0.072			Interlayer structure	0.149
		RA3	0.145			Layer structure	0.115
		RA4	0.07			Fragmented structure	0.163
		RA5	0.101			Thin layer structure	0.145
		RA6	0.103			Blocky structure	0.079
		RA7	0.075			Scattered structure	0.281
		RA8	0.087	Bedding structure	0.2829	BS1	0.086
		RA9	0.053			BS2	0.110
Cutting height/(m)	0.14	0	0.061			BS3	0.119
		0~5	0.104			BS4	0.131
		5~10	0.165			BS5	0.139
		10~20	0.251			BS6	0.179
		>20	0.418			BS7	0.234

.

Then, the weighted summation formula was used to calculate the landslide susceptibility index. It is given in Equation (3).

$$I={\displaystyle {\sum }_{i,j=1}^{m,n}{\omega }_i{x}_j},$$

(3)

where n is the number of landslide causal factors, m is the categories number of each factor, ${\omega }_i$ is the weight of the ith factor, and x_j is the normalised values of the jth category of the factor.

Subsequently, the landslide susceptibility index was divided into five grades: very high (10%), high (20%), moderate (20%), low (20%), and very low (30%). In the construction of G69, many treatment projects such as retaining structures or soil-improvement techniques have been applied, which can greatly reduce landslide hazards (Figure 10).

Finally, the method of downgrading to a very low susceptibility level if there is a treatment project was used to revise the landslide susceptibility of G69. The susceptibility maps of G69 are shown in Figure 11.

5.3. Fusion of Landslide Susceptibility Maps

The final landslide susceptibility map was produced by combining the original landslide susceptibility map and the landslide susceptibility map of G69. The original landslide susceptibility map produced by the RF model shows the best accuracy, which was applied for the combination. During the fusion, we followed the principle of taking the maximum susceptibility value. The formula is shown in Equation (4).

$$I=\left\{\begin{array}{c}{I}_{overallregion}{I}_{roadregion}<I_{overallregion}\\ I_{roadregion}{I}_{overallregion}<I_{roadregion}\end{array}\right.,$$

(4)

where I is the final susceptibility index, I_{overall region} is the original susceptibility index of the overall region, and I_{road area} is the susceptibility index of the road area. After calculating by raster calculator in ArcGIS 10.6, the final landslide susceptibility map was produced (Figure 12).

6. Discussion

6.1. Landslide Distribution Partten

Landslide patterns usually show a great distinction in different areas due to the variation of the environment and materials involved; hence, the susceptibility models often perform variously [37]. In the study area of Wanzhou country, we explored the performance of four data-driven models (LR, DT, SVM, RF) and compared them with the AHP model. The results show that the RF model performed the best in LSM than SVM, DT, LR, and AHP (Figure 9). Similar to many studies [19,28,29], the RF model can be considered an effective tool in LSM because of its fast-modelling speed, adequate robustness, and accurate performance. In the case study, four regional landslide susceptibility maps show that most of the areas along the Yangtze River are in the high-susceptibility class (Figure 8). The reason for this is that the impoundment of reservoir water raised the susceptibility of the bank areas. Furthermore, the low-susceptibility class pixels are mainly located in the middle- and high-altitude areas. These areas are covered by limestone and quartz sandstone lithology, which do not easily cause landslide hazards. Hence, the evaluation results are more in line with reality. However, the G69 areas are classified as low susceptibility from the obtained maps, which is not consistent with the results of the field investigation.

The factor importance analysis by the RF model indicates that bedding structure, aspect, and lithology contribute the most to the LSM (Figure 7). Due to frequent disturbance at the foot, including lateral erosion and excavation, dip slopes are more prone to instability [45]. It considerably improves the factor importance of aspect. To understand laws governing landslides and underlying geomorphological process, the information values of each factor’s classes were calculated by Equation (2). In different bedding structures, the IV of dip slope (BS2) was the largest 0.988 (Table 1), having a promoting effect on the development of a landslide. Alternatively, the anaclinal slope (BS7) has the least promoting effect with the IV of −2.355. The above results are in good accordance with the investigation results, as landslide activity is clearly related to the shear strength parameters of the hillslope materials and slope structures. The failure of the dip slope is generally along the structural plane or weak intercalation, whereas the failure of the anaclinal slope is generally controlled by rock joints [46,47]. In a similar geological environment, the shear strength of the material in the dip slope, in general, is weaker than the anaclinal slope, which is more prone to landslide deformation.

6.2. Comparison of the Final Map and the Original Map

This paper explores a novel procedure coupling data and a heuristic model to LSM with limited landslide inventory. To determine whether there are statistically differences between the final map and the original map, the number of grids at each susceptibility level and their proportions are listed in

Table 7. Statistical analysis of the proportion of each susceptibility class.

Susceptibility Level	Original Map (RF)		Final Map (AHP + RF)		Variation
	Pixels	Proportion	Pixels	Proportion
Very low	82, 421	68.66%	82,738	68.92%	+0.23%
Low	13, 589	11.32%	10,950	9.12%	−2.24%
Moderate	8, 415	7.01%	8251	6.87%	−0.14%
High	5, 210	4.34%	4815	4.01%	−0.25%
Very high	10, 408	8.67%	13,289	11.07%	+2.4%

. The results show that the proportion of the very-high-susceptibility class has increased by 2.4%, while the proportions of the low-, moderate-, and high-susceptibility class have decreased by 2.24%, 0.25%, and 0.14%, respectively. Therefore, this study effectively corrects the susceptibility class of the G69 area, which was underestimated due to the lack of landslide inventory in the map only obtained by the data-driven model.

A field investigation was also conducted for comparing and verifying the original and final maps. During the period from July to August 2020, two new landslides (D1 and D2) occurred along the road G69 (Figure 13). As can be seen in the final landslide susceptibility map (Figure 13), these two landslides are located in the high- and extremely-high-susceptibility levels, respectively, while in the original map produced by only the RF model, these two landslides are located in the low- and very-low-susceptibility levels (Figure 13). In view of the comparison results, we argue that a more accurate landslide susceptibility map can be obtained by modifying the AHP model. Hence, we recommend using the region-road-fusion process for LSM to overcome the limited landslides and obtain a more detailed susceptibility map, rather than focusing on creating new models or hybrid models.

In field work, we found that the impact of reservoirs water fluctuation on landslides is gradually weakening, which is in line with research by Zhang et al. [48]. At the same time, landslides induced by anthropogenic actions frequently occur in places where residents are concentrated and significantly impact people’s life, health, and economic development. Hence, in order to reduce the losses caused by landslides, we should pay more attention to the probability of landslide occurrence in places where residents are concentrated, especially in accurately calculating the landslide susceptibility after road construction, landslide engineering mitigation, and land-use changes.

7. Conclusions

In this research, a novel method has been proposed for landslide susceptibility mapping, suitable for areas with heterogeneous completeness of landslides inventories. In particular, the method combines data-driven and knowledge-driven approaches. The proposed approach has been applied to an area receiving road construction in Wanzhou, China. The data-driven approach has been applied to derive a first regional landslide susceptibility map valid for the whole area, except in the nearby of the road. The knowledge-driven method has then been used to derive susceptibility in the vicinity of the road (lower completeness of the landslide inventory). The following conclusions can be drawn from the present analysis: (a) For the regional map, among the four tested data-driven models (RF, SVM, DT, LR), RF has the greatest accuracy. (b) Due to the lack of landslide samples, the landslide susceptibility class in G69 areas cannot be accurately produced by the data-driven models. A knowledge-driven approach (AHP) is more suitable to obtain landslide susceptibility in the vicinity of road G69. (c) In the final map, combining the original regional map (RF) with the AHP map, the proportion of the very-high-susceptibility class has increased by 2.4%. Two new landslides developed along G69 are the robust evidence to prove it. Hence, the proposed fusion process for LSM has been successful in overcoming the issue of incomplete landslide inventories along human-modified areas within a broader region.