Fengjie County of Chongqing, with typical mountainous landforms, is located in the east of the Sichuan Basin, the upper reaches of the Yangtze River (Figure 1). It is the junction of the Dabashan arc fold fault zone and east Sichuan arc concave fold zone, with the sophisticated structural stress field. The range of elevation is 87–2,125 m, which is higher in the terrain in the north than that in the south in general. It is located in the subtropical monsoon climate zone, with an annual average precipitation of 1,132 mm in many years.

The data of landslides in Fengjie County were obtained from the Chongqing Geological Environment Monitoring Station. Its attributes contained the landslide name, coordinates, elevation, and time of occurrence. The 1:200,000/1:50,000 geological map of the China Geological Information Library (available online: http://ngac.org.cn/) provides a reference basis for the site investigation of landslides in terms of regional geological environment, geological structure, and stratigraphic lithology. The data of Chongqing points of interest (POI) were obtained using a python program or foot written according to certain rules that can automatically capture World Wide Web information. These activity points are standard for various types of commercial and educational activities that can represent human engineering activities, such as hospitals, primary and secondary schools, business centers, parks and squares, and so on. Also, the time range was consistent with the historical landslides. Other primary data information is shown in Table 1.

In this study, Fengjie County, a typical mountainous county in the Three Gorges Reservoir area of China, was chosen as a test site to carry out the research. The methodological flowchart is shown in Figure 2. Firstly, we selected 22 factors as the initial conditioning factors for LSM, which were derived from a satellite image, a DEM, geological data including lithology and faults, and other multisource data. Secondly, a geospatial database was established based on all data (including landslides and nonlandslides). By including geographic factor detection into a stepwise logistic regression (LR) procedure, screening out the dominant factors was then performed. Then, 10-fold cross validation was used to select the optimized training sample to generate the LSM in the study area. Then, a stable LR model with dominant factors was applied for higher accuracy and stability. Finally, the accuracy of the two models before and after optimization was evaluated and used for LSM.

Logistic regression (LR) is a generalized linear regression analysis method suitable for multivariable control. Different from the general linear regression model, the logical regression model restricts the output value to the interval [0,1] through the sigmoid function. Therefore, f(z), representing the probability of landslide susceptibility, can be expressed by the following equation: f ( z ) = 1 1 + e − z , (1) where z = w1x1+ w2x2+…+wMxM + b is a weighted linear combination model. b is a constant number, which is the intercept of a function. wM (M = 1,2,3,…,22) is the correlation coefficient of the function. The independent variable xM (M = 1,2,3,…,22) represents 22 landslide conditioning factors. When the f(z) value is 1, it means that a landslide will surely occur, and when it is 0, it means no landslide will occur (Kalantar et al., 2018).

The LR model includes the stepwise regression method to screen and eliminate the variables, thus eliminating multicollinearity and obtained significant factors. The stepwise regression is of three steps. The first step is introducing variables into the model one by one. Second, an F test is performed after introducing each explanatory variable. The final step is to perform a t-test on the explanatory variables already in the model. Once the newly added explanatory variable causes the previous explanatory variable to be no longer significant, it shall be eliminated to ensure that only essential variables are included in the regression equation before introducing each new variable.

The GeoDetector is a statistical method proposed by Wang in 2010 (Wang et al., 2010), which is used to detect spatial differentiation, revealing its driving force based on four detectors: differentiation and factor detection, interaction detection, risk-area detection, and ecological detection. This method can work both on numerical and qualitative data. In this study, to improving geospatial agreement between landslide and its conditioning factors, the factor detection function was used to screen the landslide dominant factor.

The general assumption of landslide factors screening using the GeoDetector can be expressed as follows: if a conditioning factor controls or contributes to a landslide, the geospatial distribution should be similar to a certain extent between the conditioning factor and the landslide. The similarity between the two is determined by the local and global variance (Luo and Liu, 2018). Let X be a layer of data representing the impact factors (e.g., lithology or slope) of our interests and the binary variable Y quantitatively represent the spatial distribution of the landslide. X needs to be either a category layer (for example, lithology) or can be transformed into a category zone (for example, a continuous slope can be divided into three areas: gentle, moderate, and steep). Overlaying Y and X layers subdivide Y according to the area of X.

If the factor X is related to Y, then Y will exhibit a spatial distribution similar to that of X. In the perfect case, factor X thoroughly explains the pattern of Y, the value of Y is uniform in each region of X, and the spatial variance of Y in all regions is 0. In reality, the degree of spatial correspondence between X and Y is measured by the degree of interpretation of factor X, which is defined as q = 1 − ∑ m = 1 S N m σ m 2 N σ 2 = 1 − W S S T S S (2) WSS = ∑ m = 1 S N m σ m 2 ,   TSS=   N σ 2 , (3) σ m 2 = ∑ i = 1 N m ( Y m , i − Y m ¯ ) N m − 1 , (4) σ 2 = ∑ j = 1 N ( Y j − Y ¯ ) N − 1 , (5) where m = 1,…, S is the layering of variable Y or factor X, that is, classification or partition; N_m and N are the layer m and the number of units in the whole area; and σ m 2 and σ 2   are the layer m and the variance of the Y value for the entire region, respectively. Within Sum of Squares (WSS) is the sum of the variances within the stratum, and the Total Sum of Squares (TSS) is the total variance of the entire region. Y_{m, i} is the value of Y in the ith unit of layer m, and   Y m ¯ is the average value of variable Y in layer m. Y j is the Y value of the jth unit of the whole study area, and Y ¯ is the total average of the variable Y in the whole study area. The range of q is [0, 1].

10-fold cross validation, used to test the accuracy of the algorithm, is a common test method. The dataset was divided into ten parts, nine of which were used as training data and one as test data in turn.

The correct rate (or error rate) will be obtained for each test. The average value of the accuracy of the results 10 times is used as the estimation of the accuracy of the algorithm. Generally, multiple 10-fold cross validations (for example, 10-fold cross validation) are needed, and then, the average value is calculated as the estimation of the accuracy of the algorithm.

The landslide’s occurrence is affected by both the internal geographical environmental factors and external disturbance environmental factors (Tsangaratos et al., 2017). The internal factors include topographic and geological conditions, while the external factors include human engineering activities, rainfall, and reservoir water level changes. Based on comprehensive consideration of various factors and the actual conditions of the study area, 22 factors were selected as the initial conditioning factors including topography (elevation, slope, degree of relief, aspect, slope position, landforms, curvature, profile curvature, plan curvature, terrain roughness index (TRI), topographic wetness index (TWI), sediment transport index (STI), and stream power index (SPI)), geology (lithology, distance from fault, and combination reclassification of the stratum dip direction and slope aspect (CRDS)), environmental conditions (normalized vegetation index (NDVI), distance from hydrographic net, annual average rainfall, and land cover), and human activities (distance from roads and POI kernel density). All the 22 conditioning factors mentioned above were used to create the geospatial database with GIS software (Sun et al., 2021c).

Here, slope, degree of relief, aspect, slope position, landforms, curvature, profile curvature, plan curvature, terrain roughness index (TRI) (Althuwaynee et al., 2014), topographic wetness index (TWI) (Yilmaz 2009; Hong et al., 2016), sediment transport index (STI) (Pourghasemi et al., 2012), and stream power index (SPI) (Moore and Wilson, 1992) were based on the processing of the DEM with 30 m resolution. Lithology and faults were extracted by geological maps with scales of 1:50,000–1:200,000. The distances from faults, hydrographic nets, and roads were generated after buffering the faults, river networks, and roads, respectively. The selection of buffer distance was based on field surveys, imagery resolution, and previous research (Xie et al., 2018; Sun et al., 2021a). CRDS (Sun et al., 2021b) was generated by subtraction and reclassification of aspect and tendency. NDVI was calculated using landsat8 OLI data with a resolution of 30 m. Annual average rainfall was generated by the spatial interpolation method from the original data. POI kernel density was generated by kernel density calculation of POI point data, which refers to any nongeographical meaningful points on the map (Sun et al., 2021b).

Moreover, to reduce the disadvantageous effects on the model caused by data discretization, reclassification was performed for the 13 continuous variables, including elevation, slope, degree of relief, curvature, profile curvature, plan curvature, TRI, TWI, STI, SPI, NDVI, annual average rainfall, and POI kernel density. Here, based on research experiences, the reclassification threshold value of each factor obtained by the natural breakpoint method was followed to slightly adjust by counting the number of historical landslide points under each category to agree with the actual situation as much as possible. The other nine factors, which were originally classified, could directly be processed. A classification scheme used was established for each continuous factor. To sum up, a spatial database of landslide conditioning factors after reclassification was constructed in a 30 m resolution grid cell. The thematic map of landslide influencing factors is shown in Figure 3.

To reduce the affection of the logistic regression model by factor data discreteness, all the 22 factors after reclassification were normalized to [0,1] by the following equation: X ∗ = ( X − X m i n / ( X m a x − X m i n ) , (6) where X* is the normalized data; X is the original data; X _min is the minimum value of data; and X _max is the maximum value of data.

In the study, landslide cells (positive cells) and no-landslide cells (negative cells) made up all datasets. Landslide cells consisted of 1,522 historical landslide events, each of which was regarded as a single cell. Because of the geometric effect, all landslide cells were excluded by setting a 500 m buffer zone for all the 1,522 landslide points, and the remaining areas were regarded as the no-landslide areas (Xie et al., 2018). Researchers have different opinions on the ratio between landslide and nonlandslide cells, but the most widely used ratios are 1: 1, 1: 5, and 1:10 between landslide and no-landslide cells, and the last ratio of 1:10 was selected by this study after many experiments. 15,220 nonlandslide cells were randomly extracted from the no-landslide area.

Any landslide susceptibility assessment will have no scientific significance without validation, so it is necessary to evaluate the validity of the models used. The landslide is prone to a typical binary classification problem, and the confusion matrix can be used to analyze the accuracy further. Instances are divided into positive and negative categories. That is to say, if the sample point is a landslide, it is positive, and if the sample point is a nonlandslide, it is negative. We predicted that four situations will occur in results: 1) an instance status is “landslide” and also predicted as “landslide,” recorded as True Positive (TP); 2) “nonlandslide” but predicted as “landslide,” recorded as False Positive (FP); 3) an instance status “nonlandslide” and predicted as “nonlandslide,” and it is recorded as True Negative (TN); and 4) “landslide” but predicted as “nonlandslide,” recorded as False Negative class (FN). Furthermore, we used the “OptimalCutoff” function of the “InformationValue” package of R language to find the optimal threshold for partitioning the prediction results of the model. If the predicted value is greater than the threshold, a landslide will occur. Otherwise, the landslide will not occur.

Based on the confusion matrix, Receiver Operating Characteristic (ROC) curves and Area Under the Curve (AUC) values are commonly used to comprehensively test and evaluate model accuracy (Naghibi and Moradi Dashtpagerdi, 2016). When the AUC value is more significant than 0.5 and the closer it is to 1, the higher the accuracy of the model prediction.

The optimal thresholds of LR_ao and LR_bo models were very similar, which were 0.534 and 0.592 (Table 6). The accuracy, precision, and recall of the LR_ao and the LR_bo were almost the same. The total accuracy of the two models had little difference, which was 0.911 and 0.910, respectively. Nevertheless, the accuracy of landslide and nonlandslide of the LR_ao model (nonlandslide: 0.914; landslide 0.611) was slightly higher than that of the LR_bo model (nonlandslide: 0.911; landslide 0.595). The two models had extremely high recall rates for nonlandslide, and the LR_ao (0.996) is slightly lower than the LR_bo (0.998). However, for landslides, the recall rate was lower, and the LR_ao (0.065) was slightly higher than the LR_bo (0.031). The difference between the two models was insignificant, and both had exceptional reliability.

The AUC values of the ROC curves of the training dataset of the two models were 0.843 and 0.835, the test dataset were 0.834 and 0.840, and all datasets were 0.842 and 0.835 (Figure 5), indicating that the LR model before and after optimization in this study both had high reliability and good prediction ability.

Compared with the LR_ao, the LR_bo (the optimized model) had absolute stability and high accuracy. So, we applied the LR_bo to the entire case study area for LSM. Expert empirical methods and dynamic setting threshold methods were used to select the appropriate division-level threshold. The selection of the optimal threshold should minimize the density of historical landslides in the low-susceptibility region and maximize the density of historical landslides in the high-susceptibility region. After repeated trials and errors, five susceptibility levels were finally divided (very-low–susceptibility region with p < 0.05 vs. low-susceptibility region with 0.05 ≤ P < 0.10 vs. moderate-susceptibility region with 0.10 ≤ P< 0.17 vs. high-susceptibility region with 0.17 ≤ P < 0.28 vs. very-high–susceptibility region with p>=0.28) (Figure 6). It shows the very-low– and low-susceptibility regions were concentrated in the south and southeast of Fengjie County. Meanwhile, the high- and very-high–susceptibility regions, which agreed with the geospatial of the historical landslides prone to occurring, were concentrated on both sides of the Yangtze River and its tributaries, mainly in the northern and central parts of Fengjie County.

Table 7 is the statistical result of landslide susceptibility at different levels. The proportion of historical landslides increased gradually with the increase of the susceptibility level, and the density of landslides was positively correlated with the susceptibility level. The area of very-low– and low-susceptibility regions accounted for 65.43% of the total area of the study area, while the number of historical landslides only accounted for 23.85% of the total landslides; the area of very-high– and high-susceptibility regions accounted for 20.51% of the total area, while the number of historical landslides accounted for 57.42% of the total landslides. This also shows the agreement of the geospatial between the high-susceptibility regions and the historical landslides occurred.

The historical landslide data used in this study were from 2001 to 2016. To verify the LSM results further, we also collected 61 new landslide data in the study area in 2017. All the new landslides were triggered by rainfall, which means that the main factor inducing landslides is rainfall in this study area. The geospatial distribution of the new landslides was overlaid to the LSM (Figure 7), and Table 7 shows the details including the location susceptibility level of all the new landslide events. 65% of the new landslide events were located in the high- and very-high–susceptibility regions, and 10% were in the moderate-susceptibility region in general. For case study purposes, we analyzed three typical landslides: Qulongxiaoxue landslide, Sujiawan landslide, and Tianjiadawu landslide. Qulongxiaoxue landslide was a medium-sized landslide, which was the largest of all new landslides, covering an affected area of 28.179 m². Although the area and volume of the Sujiawan landslide were not significant, the number of people affected and the number of threats were second. Tianjiadawu landslide was a medium-sized landslide that occurred in northwestern Fengjie County in October 2017. These three landslides were located in the high-susceptibility region, indicating that the LSM has a good geospatial agreement with the actual landslide events and the model had a good prediction ability.