Applicability of Susceptibility Model for Rock and Loess Earthquake Landslides in the Eastern Tibetan Plateau

: It is crucial to explore a suitable landslide susceptibility model with an excellent prediction capability for rapid evaluation and disaster relief in seismic regions with different lithological features. In this study, we selected two typical seismic events, the Jiuzhaigou and Minxian earthquakes, which occurred in the Alpine karst and loess regions, respectively. Eight influencing factors and five models were chosen to calculate the susceptibility of landslide, including the information (I) model, certainty factor (CF) model, logistic regression (LR) model, I + LR coupling model, and CF + LR coupling model. Then, the accuracy and the landslide susceptibility distribution of these models were assessed by the area under curve (AUC) and distribution criteria. Finally, the model with high accuracy and good applicability for the rock landslide or loess landslide regions was optimized. Our results showed that the accuracy of the coupling model is higher than that of the single models. Except for the LR model, the landslide susceptibility distribution for the above-mentioned models is consistent with universal cognition. The coupling models are generally better than their single models. Among them, the I + LR model can obtain the best comprehensive results for assessing the distribution and accuracy of both rock and loess landslide susceptibility, which is helpful for disaster relief and policy-making, and it can also provide useful scientific data for post-seismic reconstruction and restoration.


Introduction
Landslides are the main secondary disaster induced by seismic activities. According to the lithology of landslide material, landslides can be divided into rock landslides and soil landslides [1]. Many large seismic events occurred in the eastern and northeastern margin of Tibet Plateau, such as the Haiyuan Ms 8.5 earthquake on 16 December 1920 [2], the Diexi Ms 7.5 earthquake on 25 August 1933 [3], the Wenchuan Ms 8.0 earthquake on 12 May 2008 [4], the Lushan Ms 7.0 earthquake on 20 April 2013 [5], the Minxian Ms 6.6 earthquake on 22 July 2013 [6], the Ludian Ms 6.5 earthquake on 3 August 2014 [7], and Jiuzhaigou Ms 7.0 earthquake on 8 August 2017 [8] (Figure 1). Usually, the topography with deep-incised valleys and high mountains developed in these high seismic risk regions ( Figure 1c) can easily induce a series of coseismic landslides. The rock and loess landslides have developed in the eastern and northeastern parts of the Tibet Plateau, respectively. The earthquake-induced landslides not only threaten the safety of humans and property but also pose great safety risks to the reconstruction of the disaster area. Thus, it is important to explore and choose a model with excellent prediction capability and suitable for quickly assessing earthquake-induced landslide susceptibility in different lithology regions.
At present, landslide susceptibility methods can be divided into qualitative and quantitative analyses [9,10]. The qualitative methods mainly rely on expert experiences to estimate susceptibility [11]. In contrast, the quantitative methods are more objective because they mainly rely on data rather than expert experiences [12]. The Newmark method [13], which is one of the quantitative methods, needs to obtain clear geotechnical physical mechanics and ground motion parameters. However, it is difficult to obtain these parameters in a large area, and sometimes there are errors between the obtained parameters and the true values. So, it is hard to use the Newmark model to calculate regional landslide susceptibility in many cases [14,15]. The statistical models are indirect quantitative methods with high accuracy that have been used in a large number of studies. Previous studies show that these methods are effective in landslide susceptibility calculation [16]. The common statistical models include the frequency ratio (FR) model [17], information (I) model [18], certainty factor (CF) model [19], support vector machine (SVM) [20], weight of evidence method [21], and logistic regression (LR) model [22]. The LR model has been widely used in the study of earthquake landslide susceptibility because of its simple calculation and clear physical meaning, such as in the Wenchuan earthquake [23], the Yushu earthquake [24], and so on.
However, these statistical models still have some limitations. The SVM model needs to choose the appropriate kernel function. Some methods, such as FR, I, CF, and weight of evidence models, do not consider the differences of influencing factors and add each influence factor layer with the same weight for landslide susceptibility. Other methods, such as the LR model, can not solve the quantitative problem of each influencing factor. Hence, it is difficult to objectively and accurately calculate susceptibility to landslides on a consistent basis [25]. Coupling these single models can solve the problems of quantification and weight of influencing factors, complement their respective advantages, optimize the susceptibility model, improve accuracy, and enhance rationality; model coupling is one of the hotspots of current landslide susceptibility model research [26][27][28][29].
Moreover, due to differences in regional geological environments, different landslide susceptibility calculation models still have problems in the application process. Therefore, before calculating landslide susceptibility, it is necessary to couple these single models and select a model with high accuracy and strong applicability [30]. To select a rapid susceptibility model suitable for different lithology landslides, the Jiuzhaigou earthquake area and the Minxian earthquake area were selected as the study areas. In this study, the I, CF, LR, I + LR, and CF + LR models were used to calculate landslide susceptibility, respectively. Then, these models were assessed by the AUC value and distribution criteria. The best method with a more reasonable distribution and higher accuracy suitable for rock landslides and loess landslides could be selected. That method could quickly evaluate landslide susceptibility, helping disaster relief and reconstruction in earthquakestricken areas. This study will provide a reference for landslide susceptibility calculations in areas with similar geological environments in the future.

Jiuzhaigou Earthquake
Jiuzhaigou is located in the eastern Tibetan Plateau. It is a topographic transition zone from the Tibetan Plateau to the Sichuan Basin ( Figure 1). The main faults include the Minjiang fault, the Tazang fault, and the northwestern segment of the Huya fault in the Jiuzhaigou area. This region features deep-incised valleys and high mountains [31,32]. The maximum peak is over 4700 m a.s.l. and the minimum elevation is about 1160 m a.s.l. in the river valley. The slope gradient is higher than 30° [8]. Carbonate rocks are widely distributed in this area. There are a large number of folds and faults. Neotectonic movements are strong, and a large number of joints exist in the rocks. The rocks are relatively broken [33].  Figure 6f). The most representative sedimentary rocks are the thick sequence of deep marine deposits, including limestone, flysch complex, and sandstone. The lithologic characteristics of Devonian strata are stratiform limestone and massive dolomite; the lithologic characteristic of Carboniferous stratum is layered and dense massive limestone ( Figure 3e); the lithologic characteristics of Permian strata are composed of siliceous limestone and sandy limestone; Triassic lithologic characteristics are mainly composed of limestone ( Figure 3d) and sandstone, and Quaternary loose deposits are distributed along the river valleys (Table 1) [32,34].  20°N, 103.82°E). The area affected by the Jiuzhaigou earthquake was more than 4000 km 2 . The seismogenic fault of the earthquake was a previously unknown blind fault, the Xiongmaohai-Shangsizhai fault, which is likely the northwestern extension of the Huya fault [35]. This earthquake triggered thousands of rock landslides [8,36,37].

Minxian Earthquake
Minxian county is located on the transition zone between the Tibetan Plateau and Loess Plateau (Figure 1). The main fault is the Hetuo fault, which is a secondary fault of the Lintan-Tanchang fault zone in the Minxian area [38]. The study area is of plateau mountainous landforms. The mountain area accounts for nearly 90% of the study area. Its elevation is mostly in the range of 2190-3360 m. The slope gradient ranges from 0° to 72.3563°, 19.37° on average [39].
In the study area, the exposed strata from the surface to bottom are mainly covered by late Pleistocene Malan loess (Q3), middle Pleistocene Lishi loess (Q2), loess-like silt, weathered sandstone, and other loose sediments [40]. Field investigation showed that the loess is thicker than 10 m [41]. The structure of the Malan loess is weak and loose, and the vertical joints are well developed.
Landslides are mainly distributed in the area where the underlying lithology is tertiary red sandstones and siltstones (Guyuan Formation, Gansu Group) (Table 2, Figure 7f) [42]. The boundary between the overlying loess and the underlying tertiary red beds is the zone, which is easy to slip [6,41]. Moreover, the structure of loess under the influence of the earthquake can easily cause loess landslides, collapses, and subsidence [43].

Material and Method
The steps of landslide susceptibility calculation are listed as follows (

Acquisition of Landslide Database and Influencing Factors
For the Jiuzhaigou earthquake landslide database, we used the database established by Guo et al. (Figure 3) [32]. They integrated space-air-ground monitoring technology to extract landslides ( Figure 4). Firstly, support vector machine (SVM) classification was used to preliminarily identify the bare land (exposed ground without plants, snow, or water) within 2055 km 2 of the Jiuzhaigou area by Sentinel-2 (10 m) and Unmanned Aerial Vehicle (UAV) remote sensing images (1 m and 0.16 m). Then, visual interpretation was used to verify that the bare land represented the landslide results or not, and eliminate false detection areas. For key areas or uncertain areas, field investigations were carried out to verify the results. A total of 4456 landslides in September 2017 were extracted. The total area of the landslides was 13.7 km 2 [32].   For the influencing factors, seismic factors, topographic factors, geological factors, and human activity data were selected for landslide susceptibility calculations in previous studies (Figures 6 and 7) [14,35,[45][46][47]. They included peak ground acceleration (PGA), slope, aspect, elevation, strata, distance from fault, river, and road. Among them, terrain data such as slope, aspect, and elevation were obtained by ALOS DEM (12.5 m). The existing tools of ArcGIS were used to extract topographic parameters. Spatial Analyst Tools-Surface-Slope tool was used for slope extraction. Spatial Analyst Tools-Surface-Aspect tool was used for aspect extraction. Strata and fault data were from the 1:200,000 geological map (http://dcc.ngac.org.cn/geologicalData, accessed on 8 May 2021). Seismic data such as PGA were provided by the USGS website (http://earthquake. usgshakemap.gov, accessed on 8 May 2021). The data on human activities such as roads and rivers came from the national data of roads and rivers (http://www.gscloud.cn/sources, accessed on 8 May 2021).

Preprocessing of Model Input Data
For Jiuzhaigou landslide susceptibility calculation, 3000 landslide training samples and 1000 landslide testing samples were randomly selected from the 4456 landslide points as positive samples of the model. In the non-landslide area, 4456 non-sliding points were randomly selected from outside the landslide area. A total of 3000 training samples and 1000 testing samples were randomly selected from the 4456 non-sliding points as negative samples of the model. Finally, we obtained the training dataset containing 6000 training samples and the testing dataset containing 2000 testing samples.
For Minxian landslide susceptibility calculation, 1650 landslide training samples and 550 landslide testing samples were randomly selected from the Minxian earthquake landslides database. A total of 2330 non-sliding samples were randomly selected in non-landslide areas, and 1650 training samples and 550 testing samples were randomly selected from non-sliding samples. Finally, we obtained the training dataset containing 3300 training samples and the testing dataset containing 1100 testing samples.
For influencing factors, all samples were converted into raster form with a grid cell size of 12.5 m × 12.5 m, which was used to calculate the area and extract the attribute values of sample points by ArcGIS tools in the next step. Then, this study classified the selected factors (Figures 6 and 7). For example, the aspect was divided into nine categories: flat, east, southeast, north, northeast, west, northwest, south, and southeast ( Figures  6b and 7b). Finally, we calculated the area and the landslide area density of each category. The specific formula is as follows [34]: where is the landslide area density; is the factor (slope, aspect, elevation, PGA, fault, strata, river, and road); is a subset of ; is the landslide area of the category " " in the factor " "; is the area of the category " " in the factor " ".

Landslide Susceptibility Model
The influencing factors of the training dataset and the testing dataset, which were classified and calculated in Section 3.2, will be used in the following susceptibility model.

Information Model (I)
In the information model, the landslide probability is represented by the decrease of entropy in the landslide prediction process. The average reduction of uncertainty caused by the combination of the factors is equal to the entropy change of the landslide system. The information model considers that the landslide probability is related to the quantity and quality of information obtained in the process of prediction, which is measured by the information value. The greater the information value, the greater the landslide possibility [48]. The specific formula is as follows [49]: where ( , ) is the information value of the category " " in the factor " " contributed to the landslide event ; is the landslide area in the category " "; is the total area of landslides in the whole study area; is the area of the category " "; is the total area of the whole study area. The larger the information value, the more likely a landslide is to occur. Conversely, the smaller the information value, the less likely a landslide is to occur.

Certainty Factor Model (CF)
The CF model is a probability function. It was first proposed by Shortliffe and Buchanan [50] to analyze the susceptibility of various factors affecting the occurrence of an event. Lan et al. applied the CF method to evaluate the landslide susceptibility [19]. The formula is listed as follows: where is the conditional probability of landslide occurrence in category " "; is the prior probability of landslide occurrence in the whole study area; is the landslide area density of the category " "; is the landslide area density of the total landslide. The range of is [−1, 1]. When > 0, the larger is, the more likely the slope will slide; when < 0, the smaller is, the less likely the slope will slide; when = 0, it indicates that the probability of the slope sliding may not be determined.

Logistic Regression Model (LR)
The LR model predicts the probability of events by establishing a regression relationship between a dependent variable and multiple independent variables. It does not need to meet the assumption that the error distribution tends to be the normal distribution, nor does it need to meet the condition that the independent variables conform to the normal distribution. It can be used to predict the probability of dependent variables with binomial characteristics. So, it is more suitable for the landslide susceptibility assessment [34,51]. In analysis, the dependent variable is a binary variable that represents whether there is a landslide. The value "1" and value "0" of indicate the occurrence and non-occurrence of a landslide, respectively. Z is the intermediate parameter, varying from −∞ to +∞. There are independent variables ( , , ,…, ) in an LR model. Influenced by factors, the conditional probability of a landslide is = ( = 1 | , , ,…, ). Then the model can be expressed as [52,53]: where Z is the linear combination; ( = 1, 2, 3,…, n) is an independent factor; is the regression constant; is the coefficient that measures the contribution of independent factor ( = 1, 2, 3,…, n); is the prediction value of landslide probability.

Coupling Model I + LR
In the I + LR coupling model, the I values calculated by the formula (2) substitute for the of formula (4). The data of the training dataset are imported into SPSS 25 software to calculate the constant and regression coefficients of the logistic regression. The constant and regression coefficients are substituted into formula (5), and then P is obtained.

Coupling Model CF + LR
In the CF + LR coupling model, the CF values calculated by the formula (3) substitute for of the formula (4). The data of the training dataset are imported into SPSS 25 software to calculate the constant and regression coefficient of the logistic regression. The constant and regression coefficients are substituted into formula (5), and then P is obtained.

Rationality and Accuracy Verification of Susceptibility Results
In this study, the Raster Calculator tool of ArcGIS was used to obtain the landslide susceptibility distribution map. Then the landslide susceptibility of the Jiuzhaigou earthquake and Minxian earthquake were both divided into five susceptibility classes using Jenks natural breaks optimization: very low, low, moderate, high, and very high. The Jenks natural breaks optimization follows the law of data, which can make the classification more objective. Many studies have shown that this classification method is more reasonable [17,34].
In each susceptibility class, the percentage of landslides number, the percentage of the susceptibility classification area, and the landslide number density (LND) were counted and calculated. We analyzed the rationality of the landslide susceptibility distribution by two criteria in each model. For the susceptibility distribution map, we evaluated the result of the susceptibility distribution map by: (1) Landslides should locate in the areas with high susceptibility as much as possible; (2) The areas with very high susceptibility in the distribution map should be as small as possible [54].
Additionally, this study quantitatively evaluated the accuracy of the landslide susceptibility model by the ROC curve. The ROC curve is regarded as an effective method to evaluate the prediction ability of models [55]. The abscissa of the curve is the specificity or the false positive rate (FPR), which is the probability that the non-landslide is judged as a landslide by the model; the ordinate is the sensitivity or the true positive rate (TPR), which is the probability that the landslide is judged as a landslide. The area under curve (AUC) is a quantitative index of the ROC curve to evaluate the model. The AUC value ranges from 0.5 to 1. The closer the AUC value is to 1, the better the prediction effect of the model is. If AUC = 0.5, the prediction result of the model is random; if AUC is at 0.5-0.7, the accuracy of the model is low; if AUC is at 0.7-0.9, the accuracy of the model is high; if AUC > 0.9, the accuracy of the model is very high [56]. In this study, the data of the testing dataset and model prediction results were input into SPSS 25 software, and the ROC curves and AUC values of each model were output.

Distribution Maps of Landslide Susceptibility
The landslide susceptibility maps of Jiuzhaigou (I, CF, LR, I + LR, and CF + LR) show the following (Figure 8): (1) Except for the LR model, the percentages of landslide points in the very high susceptibility area were the largest, which were 74.02%, 56.01%, 73.27%, and 57.28%, respectively. Except for the LR model, the percentage of landslide points decrease with the decrease of susceptibility, which was consistent with the evaluation criteria in each model. (2) The area proportions of the very high susceptibility class in each model were small, but the area proportions of moderate, low, and very low susceptibility classes were large, which also met the evaluation criteria. (3) The LNDs of the very high susceptibility area were the largest, which were 21.94/km 2 , 26.47/km 2 , 26.64/km 2 , 22.73/km 2 , and 27.26/km 2 . The LNDs of the coupling model were increased compared to their single models and the LR model. The LND decreased with the decrease of susceptibility in each model, which was consistent with the universal cognition ( Table 3).
The landslide susceptibility maps of Minxian (I, CF, LR, I + LR, and CF + LR) show the following (Figure 9): (1) The percentages of landslide points in the very high susceptibility area were the largest, which were 64.72%, 45.41%, 53.52%, 60.73%, and 58.97%, respectively. The percentage of landslide points decreased with the decrease of susceptibility, which was consistent with the evaluation criteria in each model. (2) The area proportions of the very high susceptibility class in each model were small, but the area proportions of moderate, low, and very low susceptibility classes were large, which also met the evaluation criteria. (3) The LND in the very high susceptibility area were the largest, which were 32.96/km 2 , 36.48/km 2 , 36.31/km 2 , 37.56/km 2 , and 39.35/km 2 . The LND of the coupling model was increased, compared to their single models and the LR model. The LND decreased with the decrease of susceptibility in each model, which was consistent with the universal cognition (Table 4).

Accuracy of Landslide Susceptibility Models
For the ROC curves of the Jiuzhaigou earthquake landslide (Figure 10), the AUC values of I, CF, LR, I + LR, and CF + LR models were 0.911, 0.895, 0.899, 0.915, and 0.909, respectively. The prediction accuracy of the I, I + LR, and CF + LR models were higher. For the ROC curves of the Minxian earthquake landslide, the AUC values of I, CF, LR, I + LR, and CF + LR models were 0.841, 0.837, 0.836, 0.844, and 0.845, respectively. The prediction accuracy of the I, I + LR, and CF + LR models were also higher. In single models, the prediction accuracy of the I model was the highest. In the accuracy improvement, the prediction accuracy of the coupling model was generally higher than that of the single models. Overall, the comprehensive effect of the I + LR model was the best.

Discussion
The fault, slope, and lithologic characteristics play important roles in the spatial distribution and susceptibility of landslides. The Jiuzhaigou landslides were mainly distributed in areas of limestone and dolomite. Most landslides were concentrated within 2 km from the fault (NW) (Figure 6e) [8,14,32]. The landslides were mainly distributed on the slopes of 30°−55°. The Minxian landslides were concentrated in an elongated region, with a 5 km width, parallel with the fault strike (NWW) (Figure 7e). The underlying sedimentary rocks were mainly conglomerate and sandstone. Most of the landslides were distributed on the slopes of 5°−25°.
For landslide susceptibility calculation, the previous studies suggest that I and CF models do not consider the differences of influencing factors on landslide susceptibility. It is unreasonable to add each influence factor layer with the same weight. The LR model can not solve the quantitative problem of each influencing factor, which is sometimes called the multi-source data combination problem. The quality of data quantitative results is directly related to the reliability of the results. So, it is difficult to objectively and quantitatively calculate landslide susceptibility with the I, CF, and LR models [28,57].
For evaluation of the model, the AUC value is used to evaluate the accuracy of the landslide susceptibility model. However, it is unable to evaluate whether the susceptibility distribution is reasonable. The model must be evaluated not only from the accuracy but also from the rationality of the susceptibility distribution. The evaluation criteria proposed by Tolga are consistent with universal cognition. So, the model should be evaluated from the AUC value and the distribution criteria.
For the I + LR and CF + LR models, the I model and CF model quantify the influence of each factor grade on landslides into a unified I value or CF value. They not only can contribute to the combination of subsequent influencing factors but also can specifically represent the susceptibility of landslides relative to the whole study area in the classification of each factor. Then, the I or CF values of influencing factors in the training dataset are put into the LR model to obtain the regression coefficients, which are the weights of the influencing factors. They consider the differences of influencing factors on landslide susceptibility.
The coupling models are affected by the single models that make them up. However, the coupling models are generally better than their single models. In the calculation process of the coupling models, the I and CF values of each factor grade are imported into the LR model rather than their original data. For example, the values of the slope or the distance to the fault are replaced by their I or CF values. The calculation formula for I and CF is simple and easy to implement, which makes the implementation of the coupling model very simple.
The coupling model in this study reasonably integrates the advantages of different models, solves the problems of quantification and weight (combination) of multi-source influencing factors for landslide susceptibility, and reasonably describes the complex nonlinear relationship between natural phenomena. Overall, the coupling models not only have higher accuracy than their single models, but fit the rules proposed by Tolga for landslide susceptibility evaluation. They can objectively and accurately evaluate the susceptibility of earthquakes and landslides in rock and loess areas.

Conclusions
In this study, the rock landslide and loess landslide related to the 2017 Jiuzhaigou earthquake and the 2013 Minxian earthquake, which have different lithologies, respectively, were chosen as the study targets. Eight influencing factors, including PGA, slope, aspect, elevation, stratum, distance from the fault, river, and road, were employed to calculate landslide susceptibility by the I, CF, LR, I + LR, and CF + LR models. The models were evaluated by the AUC value and the distribution criteria to choose the best one. The main results and conclusions obtained from this research were as follows: (1) The Jiuzhaigou landslides were mainly distributed in areas of limestone and dolomite. The Minxian landslides were mainly distributed in the area where the underlying bedrocks are conglomerate and sandstone. (2) The influencing factors adopted to calculate the susceptibility of rock landslides or loess landslides were reasonable, and the high AUC value suggested that they are suitable for the universal model.
(3) For the distribution of susceptibility, most models fit the rules proposed by Tolga [56], except for the LR model. The landslide susceptibility distribution map calculated from the coupling models was more reasonable than that derived from their single models. (4) For prediction accuracy, the coupling models were generally more accurate than their single models. The prediction accuracy of the I + LR model was high in both rock and loess areas, which have high or moderate ground motion parameters.
In general, the comprehensive effect of the I + LR model was the best one. The coupling models had excellent prediction capability and strong robustness, and not only showed higher accuracy but also fit the rules proposed by Tolga for landslide susceptibility. They could objectively, accurately, and quickly calculate susceptibility to earthquake landslides in both rock and loess areas.
Author Contributions: Conceptualization, X.G. and B.F.; methodology, X.G., Q.C., and W.Z.; writing-original draft preparation, X.G.; writing review and editing, B.F., J.D., and P.S. All authors have read and agreed to the published version of the manuscript.