Rainfall Induced Shallow Landslide Temporal Probability Modelling and Early Warning Research in Mountains Areas: A Case Study of Qin-Ba Mountains, Western China

Abstract

The rainfall-induced landslide early warning model (LEWM) is an important means to mitigate property loss and casualties, but the conventional discriminant matrix-based LEWM (DLEWM) leaves room for subjectivity and limits warning accuracy. Additionally, it is important to employ appropriate indicators to evaluate warning model performance. In this study, a new method for calculating the spatiotemporal probability of rainfall-induced landslides based on a Bayesian approach is proposed, and a probabilistic-based LEWM (PLEWM) at the regional scale is developed. The method involves four steps: landslide spatial probability modeling, landslide temporal probability modeling, coupling of spatial and temporal probability models, and the conversion method from the spatiotemporal probability index to warning levels. Each step follows the law of probability and is tested with real data. At the same time, we propose the idea of using economic indicators to evaluate the performance of the multilevel LEWM and reflect its significant and unique aspects. The proposed PLEWM and the conventional DLEWM are used to conduct simulate warnings for the study area day-by-day in the rainy season (July-September) from 2016 to 2020. The results show that the areas of the 2nd-, 3rd-, and 4th-level warning zones issued by the PLEWM account for 60.23%, 45.99%, and 43.98% of those of the DLEWM, respectively. The investment in issuing warning information and the losses caused by landslides account for 54.54% and 59.06% of those of the DLEWM, respectively. Moreover, under extreme rainfall conditions, the correct warning rate of the PLEWM is much higher than that of the DLEWM.

1. Introduction

Landslides have caused severe casualties and economic losses worldwide over the years, and China is among the countries most severely affected by landslides [1,2,3,4]. Rainfall is the primary triggering factor of landslides [5,6]. By analyzing the relationship between rainfall and landslide occurrence, establishing a rainfall-induced landslide early warning system can contribute to issuing useful warning information and mitigating landslide risk [6,7].

Landslide early warning system is a nonstructural passive mitigation risk measure that can generate timely and meaningful information to enable people threatened by landslides to take timely and appropriate action to reduce risk [6]. It is generally considered that the system consists of four interrelated modules, namely, setting, landslide prediction model, warning strategy, and response [6,8,9,10,11]. More generally, the landslide prediction model combined with the warning strategy constitutes the landslide early warning model (LEWM), which is the core part of early warning system [6].

In most cases, the LEWM was achieved by defining rainfall thresholds. The rainfall threshold can be divided into physically based and empirically based thresholds. Physically based thresholds are usually difficult to apply over large areas because detailed knowledge of meteorological, hydrological, and other parameters is not easy to acquire [12,13]. The empirical-based rainfall threshold can take advantage of extensive historical landslide and rainfall data to explore the intrinsic relationships between landslide occurrence and rainfall conditions and is more suitable for large areas [14,15,16,17]. Since Caine [18] used a power-law function to fit 73 rainfall-induced landslides, defining the first global rainfall threshold, the power-law function has become the first choice for that. Several researchers have proposed different methods to determine the parameters of the power-law function and assess the uncertainty of the threshold [19,20,21,22,23,24].

However, power-law-based rainfall threshold may have its practical limitations. First, it is only based on the rainfall data that trigger landslides (triggering rainfall) and ignores the rainfall data that did not trigger landslides (nontriggering rainfall), so part of the information is not included [17,25,26]. Second, rainfall threshold is defined as the level or the value that must be exceeded to produce a given effect or result [17], which means that it can be used to perform categorical forecasts, such as distinguishing whether rainfall can induce landslides or the hazard level of rainfall triggering landslides. However, in general, the distinction between triggering and nontriggering rainfall is not trivial, and it is more appropriate to describe the relationship between rainfall and landslide occurrence by means of probabilities [16,17,27,28]. Third, the power-law function is a relatively simple assumption, which may not fully reflect the complex relationship between rainfall and landslide occurrence in some cases [6,29].

Rainfall threshold can aid in predicting the hazard level of rainfall triggering landslides, but the spatial location of the potential landslides is not explicit, and accurate landslide early warning requires that the model contain geo-environmental information [30,31]. There are two main methods used. The first and most commonly used one is called the “fixed warning zone” [12]. In this method, the entire area is divided into multiple homogeneous zones according to the topography, lithology, elevation, and other factors, in which the spatial probability of landslides is considered to be approximately the same, and for each zone, there is a specific threshold for issuing warning information so that the spatial location for landslide occurrence is roughly explicit, such as the landslide early warning system in Hong Kong, Seattle, Rio de Janeiro, etc. [6,10,11,32,33,34]. The second method utilizes information from the landslide susceptibility map, as it is generally believed that the higher the landslide susceptibility index (LSI) is, the more frequently landslides will occur [35,36,37,38,39]. The LEWM is obtained by combining the landslide susceptibility map and rainfall threshold based on a discriminant matrix (also called a heuristics matrix) or expert knowledge. For example, the LEWM in Zhejiang Province [40,41] uses a 4 × 3 discriminant matrix to combine the susceptibility map (1:100,000 scale, divided into four susceptibility levels) and two rainfall thresholds to issue five levels of warning information. Norway’s LEWM is built based on expert knowledge; a landslide expert on duty (as a member of a rotation team) uses real-time monitoring and forecast data to compare with a predefined threshold and combines it with a susceptibility map to perform a heuristic warning at four levels [42,43,44].

The two mentioned LEWM construction methods, regardless of the method used to select homogeneous zones or design a discriminant matrix, leave room for subjectivity and impair the results comparison, and there are also some other practical limitations, such as the low spatial resolution, the relatively conservative nature of the model, the fact that the area of the issued high-level warning zone accounts for a large proportion, and that the early warning system is difficult to update [6]. In their landmark review paper, Guzzetti et al. [6] clarified that LEWM require an understanding of the areas potentially affected by landslides, and susceptibility maps are an important information source. However, the combination of susceptibility information and rainfall thresholds should follow the laws of probability, e.g., the probabilities should be conditionally independent, or the marginal probabilities should be known and treated appropriately [6].

Inspired by Guzzetti et al. [6], in this research, we attempt to develop an LEWM in a fully probabilistic way. The problem of defining a rainfall threshold is implicitly linked to the problem of estimating the temporal probability of landslide occurrence and choosing an acceptable probability value [17]. Thus, we propose to develop a temporal probability model using frequency analysis combined with machine learning fitting. In this method, we considered both triggering and nontriggering rainfall, which were identified with the same criteria. In addition, to model the spatial probability model of landslides, three most commonly used machine learning algorithms [45,46,47] were used based on 12 geoenvironmental factors to obtain the susceptibility map, the frequency ratio model was used to correct the resulting the LSI to obtain the spatial probability of landslides. Subsequently, the independence of the spatial and temporal probability of landslides was tested, and Bayesian approach was employed to couple the spatial and temporal probability models [48,49], which were used to calculate the spatial-temporal probability index of landslides (P_st). The P_st was used as an indicator to issue various level warning information, and the probabilistic-based LEWM (PLEWM) was realized.

The main purpose of our work was to propose a novel PLEWM and to verify whether it performs better than the conventional method. Therefore, we also employed the method proposed by Brunetti et al. [20] to define four percentage rainfall thresholds, and a discriminant matrix was adopted to combine the thresholds and susceptibility map to develop a conventional discriminant matrix-based LEWM (DLEWM) for comparison.

Performance evaluation plays an important role in the process of constructing a LEWM, at the same time, it is a complex process because many aspects need to be considered [50]. The most typical method is to construct a 2 × 2 contingency table, which considers the warning event and landslide event as dichotomous variables, and statistical indicators, such as the missed warning rate, correct warning rate and false warning rate, are then used to evaluate the model performance [51,52,53,54]. However, the PLEWM and DLEWM applied in this research are multilevel warning models, and the warning zone is not fixed. The existed performance evaluation method cannot reflect the important and unique aspects of the models, such as: (1) the difference between each warning level; (2) the number of landslides located in each warning zone; (3) the area of each warning zone; and (4) the interaction between the landslide number, area of the warning zone and warning level [55].

Considering that many researchers and organizations regard an early warning system as a cost-effective risk mitigation measure [6,12,56,57,58,59,60,61], the risk caused by geohazards can also be measured by economics [56,62,63]. Therefore, we think the most ideal is this kind of warning model, they issue a small area of high-level warning zones, but can predict most landslides because the small investment results in a large gain. Simultaneously, economic indicators can quantify the four mentioned aspects of the multilevel warning model. Therefore, in this research, we propose economic indicators, such as investment and loss, to evaluate early warning model performance. At the same time, we simulated warnings for the study area day-by-day in the rainy season (July–September) from 2016 to 2020 through PLEWM and DLEWM, and compared these two models.

The structure of this manuscript is outlined as follows: Section 2 describes the methodologies used in this research; in Section 3, we give the experimental results, including landslide spatial probability modeling results, temporal probability modeling results, warning strategies and simulated warning results; Section 4 and Section 5 provide the discussion and conclusions.

2. Methodology

The construction process of the PLEWM and DLEWM can be separated into five steps: (1) construction of the landslide and geoenvironmental factor database; (2) landslide spatial probability modeling; (3) landslide temporal probability modeling and definition of the rainfall threshold; (4) construction of the PLEWM and DLEWM; and (5) verification and comparison of the two warning models. The detailed modeling flowchart is shown in Figure 1.

2.1. Introduction to the Qin-Ba Mountain Area

The Qin-Ba Mountain area covers 83,200 km², and is located in the southern part of Shaanxi Province, China (Figure 2). The area is characterized by high altitude in the south and north and low altitude in the middle region. Depression basins formed by stratigraphic faults are widely distributed in the central area, which is the main human settlement area. At the same time, the study area is located at the junction of the North China Plate and the Yangtze Plate, after hundreds of millions of years of plate extrusion and multistage tectonic activity, resulting in abundant fissures, joints, and folds in the rock. Thus, the rock in the study area always has a low mechanical strength and is prone to weathering. These poor structures can form unstable slopes that are prone to collapse or supply thick and loose gravelly soil with a few meters that is prone to sliding.

Qin-Ba Mountain is the dividing line between the northern and southern climates of China, and the area is influenced by a continental monsoon climate in the warm temperature zone. The towering mountains in the north trap the warm and humid monsoon blowing from the south, resulting in heavy and continuous rainfall from July to September each year. Due to the unique rainfall patterns and complex geological conditions, the study area has become one of the most active areas for landslide occurrence in China.

2.2. Construction of the Landslides and Geoenvironmental Factor Database

Based on the literature review of landslide susceptibility and the available data in the study area, four categories (hydrological, geological, morphological and land cover) with a total of 12 geoenvironmental factors were selected, namely: distance to river, topographic wetness index (TWI), lithology, fault density, elevation, slope angle, aspect, profile curvature, soil type, distance to road, land cover type, and NDVI, which were shown in Figure 3. Among them, the DEM used in this study was obtained from the ALOS PALSAR dataset, with a spatial resolution of 12.5 m, and downloaded using the Google Earth engine platform (https://developers.google.cn/earth-engine, accessed on 1 July 2022). Four other factors, such as slope angle, aspect, plan curvature, and TWI, were calculated from the DEM. It should be noted that the aspect factor was calculated by the “aspect analysis” tool of ArcGIS pro, and the value is distributed from 0 to 360°. It was divided into 2 categories, namely the southern slope (90° to 270°) and the northern slope (0° to 90°, and 270° to 360°). The river network and road network data were extracted from the OpenStreetMap project (https://download.geofabrik.de/, accessed on 1 July 2022). Lithology and faults were extracted from the geological map at a scale of 1:500,000, and the geological map was obtained by manually digitizing a paper map. The lithology was grouped into 9 classes according to the characteristics related to landslide occurrence, such as bedding thickness, strength, weathering degree, and composition. Land cover data were obtained from the open data source supported by ESRI [64], with a spatial resolution of 10 m; NDVI factors were calculated in the Google earth engine platform using ESA sentinel-2 images, with a spatial resolution of 10 m. The soil type factor was obtained from 1: 1,000,000 soil maps of the People’s Republic of China, which is provided by the Resource and Environment Science and Data Center (http://www.resdc.cn, accessed on 1 July 2022). Finally, all of the factors were converted into rasters and resampled to a spatial resolution of 30 m.

The successful development of a PLEWM is dependent on preparing a complete landslide dataset with a long-time span. The Shaanxi Institute of Geo-environmental Monitoring maintains a complete inventory of geological hazards in the study area, which contains two datasets, one for historical geological hazards and the other for high-risk slopes (HRS), which have been certified individually by experts. In the study area, rainfall-induced shallow landslides are the most common type, accounting for more than 80% of the total number of geohazards. A total of 3127 landslides were recorded over the past 20 years (2001–2020), part of them reported by local authorities and partly from remote sensing interpretation after each extreme rainfall and confirmed in the field (about 700 landslides in this part), 2843 landslides with exact occurrence times and spatial locations were selected. Simultaneously, a total of 7436 HRSs were also contained in the database. HRSs were used to train and test the machine learning models for the landslide susceptibility map. Landslides were used to validate the corrected susceptibility model, develop a landslide temporal probability model, and define the empirical rainfall threshold.

There are 41 count-level weather stations that are homogeneously distributed over the study area, which can provide daily rainfall data from 2001 to 2020, and we downloaded these data from the National Meteorological Center (https://data.cma.cn/, accessed on 1 January 2021).

2.3. Landslide Spatial Probability Model

In previous related studies, researchers used the LSI to express landslide spatial probability. Machine learning models have advantages in landslide susceptibility modeling over large regions because they can take advantage of extensive historical landslide data to explore the complex nonlinear relationship between landslide occurrence and geoenvironmental factors. Reichenbach et al. [45] and Merghadi et al. [46] reviewed the application of machine learning models in the landslide susceptibility assessment in detail. According to their review, we choose three most commonly used machine learning models, namely, the logistic model (LR), support vector machine (SVM), and artificial neural network (ANN), to conduct landslide susceptibility assessment in this research.

Considering that the value ranges, meanings and data types of different geoenvironmental factors differ, preprocessing is necessary to ensure that all the factors have the same meanings, value ranges, and data types. The frequency ratio model is commonly used to address this issue and can be used to quantitatively express the correlation of each subclass of factors and landslides. The frequency ratio value can be calculated by Equation (1):

$$FR=\frac{N_i/N}{A_i/A}$$

(1)

where N_i is the landslide number or landslide area within the i-th subclass of a certain geo-environment factor, N is the landslide number or landslide area in the whole study area, A_i is the area of the i-th subclass of the certain geo-environment factor, and A is the area of the whole study area, and FR is the frequency ratio value and is greater than zero, FR >1 indicates that the subclass of factors is conducive to landslide occurrence. The larger the FR is, the greater the contribution to landslide occurrence, and vice versa.

In different machine learning models, the model architectures, adopted assumptions, optimization algorithms, and mechanisms for generating probability values are different from each other; as a result, the LSIs generated from machine learning models are typically very different from each other [46,65,66,67,68]. However, few studies have verified whether the LSI is consistent with the spatial probability of landslides. Considering that N_i/A_i is the conditional probability of landslide occurrence in the i-th subclass of a certain factor, and N/A is the prior probability, FR actually refers to the average spatial probability of landslide occurrence in the i-th subclass of a certain factor. Therefore, in this research, the frequency ratio model was used to verify the LSI generated from the machine learning model. If there was a linear relationship between the LSI and FR, the LSI was selected to express the landslide spatial probability; otherwise, the landslide susceptibility map needs to be corrected according to the fitted function between the LSI and FR.

2.4. Landslide Temporal Probability Model

Defining the rainfall threshold and temporal probability model includes the following steps: (1) rainfall-induced landslide inventory; (2) identification of the triggering and nontriggering rainfall events; (3) selection of the optimal combination of rainfall variables; and (4) definition of the rainfall threshold and temporal probability model. The detailed process is shown in Figure 1b.

2.4.1. Identification of Rainfall Events

To identify triggering and nontriggering rainfall events, researchers generally use cumulative rainfall (E_T) over a period of time (D_T) as a threshold to truncate long-term rainfall sequences into multiple independent rainfall events [22,24,69,70,71]. We used Python language to compile a complete processing flow, which can accept input variables D_T, E_T, rainfall record, and landslide data, and output triggering and nontriggering rainfall events, and the detailed processing can be seen in the Supplementary Materials (Figure S1). We refer to the study of Berti et al. [17], and set a threshold of D_T = 3 d, and E_T = 5 mm. Ultimately, we identified 941 triggering rainfall events and 114,908 nontriggering rainfall events.

2.4.2. Effective Rainfall Model

Rainfall-induced landslides refer to the fact that due to the infiltration of rainfall, the water content of the landslide body increases and the suction force of the matrix decreases, which leads to a decrease in the shear strength and instability of the landslide. However, due to the effect of runoff, discharge, and evapotranspiration processes, only a portion of the rainfall is involved in the process of triggering the landslide. This part of the rainfall that infiltrates the hillslope is called effective rainfall [67,72,73,74].

The effective rainfall is related to landslide occurrence more than the cumulative rainfall. There are different models to simulate effective rainfall [29,73,75], and the most commonly used model is the power-law-based model, which can be written as (Equation (2)):

$$EE={\displaystyle \sum _{i=0}^D{K}^i\cdot R_i}$$

(2)

where EE is the cumulative effective rainfall, R₀ is the rainfall on the day of landslide occurrence, and R_n is the daily rainfall on the nth day before the landslide occurrence. K is the decay factor representing the outflow of the regolith [29]. It is commonly determined empirically, varies in different areas, and is generally considered to be in the range of 0.7–0.9 [29,67,73]; when K = 1, EE is the same as the cumulative rainfall.

2.4.3. Optimal Rainfall Variable Combination Selection Based on Sensitivity Analysis

The rainfall event is characterized by various rainfall variables, such as the rainfall duration (D), cumulative rainfall (E), average rainfall intensity, daily rainfall, normalized cumulative rainfall, antecedent rainfall (AE), and EE (more details can be found in Guzzetti et al. [14]), whereas due to the lack of a reasonable rainfall variable selection method, the proposed rainfall thresholds are mainly defined by I-D or E-D. Simultaneously, for rainfall variables, such as EE, there is a parameter K that varies regionally; as a result, it is necessary to use a reasonable method to determine the optimal rainfall variable combination and parameter K to define the optimal rainfall threshold.

Sensitivity analysis studies the relationships between the variation in the input variable values and the response (output) variation in a mathematical model [76]. Sensitivity analysis methods can quantify the influence of parameter variation on the model response to achieve the purpose of variable selection. The basic motivation of sensitivity analysis -based rainfall variable selection is to compare the difference in the probability of a specific rainfall variable combination in triggering and nontriggering rainfall events; a large difference indicates the high significance of the considered variable combinations. The equation is expressed as follows (Equation (3)):

$$SAI=d\left(P\left(R_1,R_2\left|L\right.\right), P\left(R_1,R_2\left|NL\right.\right)\right)$$

(3)

where SAI is the sensitivity index for rainfall variable combination R₁–R₂, L, and NL are landslide and nonlandslide events, respectively, and d is the means by which to measure the difference between the two probability distributions. In this research, we use Jensen–Shannon divergence because it is symmetric and suitable for probability distributions, which is expressed as follows:

$$d\left(p\left(x\right),q\left(x\right)\right)=\frac{1}{2}{\displaystyle \sum _{i=1}^n\left(p\left(x_i\right)\cdot \log \left(\frac{2p\left(x_i\right)}{p\left(x_i\right)+q\left(x_i\right)}\right)\right)}+\frac{1}{2}{\displaystyle \sum _{i=1}^n\left(q\left(x_i\right)\cdot \log \left(\frac{2q\left(x_i\right)}{p\left(x_i\right)+q\left(x_i\right)}\right)\right)}$$

(4)

where $p\left(x\right)$ and $q\left(x\right)$ are the two probability distributions to be calculated for divergence. For continuous probability distributions, discretization is needed, and n refers to the amount of discretization.

2.4.4. Temporal Probability Model of Landslide Occurrence

Bayesian analysis can be used to evaluate the conditional probability of event A under the condition of event B. The formula used to calculate the conditional probability of landslide occurrence under the condition of R₁-R₂ is given as follows:

$$P\left(L\left|R_1,R_2\right.\right)=\frac{P\left(R_1,R_2\left|L\right.\right)\cdot P\left(L\right)}{P\left(R_1,R_2\right)}$$

(5)

where $P\left(R_1,R_2\left|L\right.\right)$ is the likelihood of a specific R₁-R₂ combination under the condition of landslide occurrence; $P\left(R_1,R_2\right)$ is the marginal probability; $P\left(L\right)$ is the prior probability of landslide occurrence, which is a constant; and $P\left(L\left|R_1,R_2\right.\right)$ is the posterior probability of landslide occurrence. Finally, we use ANN to fit the posterior probability to obtain the landslide temporal probability model.

2.4.5. Power-Law-Based Rainfall Threshold

As mentioned in the Introduction, power-law functions are often used to define rainfall threshold, and the formula is shown in Equation (6), where $\alpha$ and $\beta$ are fitting parameters, which can be determined by fitting the triggering rainfall; details can be found in Brunetti et al. [20].

$$R_1=\alpha \cdot R_2{}^{\beta }$$

(6)

In this study, we defined a temporal probability model and rainfall threshold using landslide and rainfall data from 2001 to 2015, which contain 2347 landslides, 762 triggering rainfall events, and 85,890 nontriggering rainfall events. Landslide and rainfall data from 2016 to 2020 were used for verification, including 496 landslides, 179 triggering rainfall events, and 29,018 nontriggering rainfall events.

2.5. Landslide Early Warning Model

2.5.1. Probabilistic-Based Landslide Early Warning Model

The PLEWM refers to the spatial-temporal probability model of landslide occurrence, where the temporal probability of landslides can be expressed by Equation (5), and the spatial probability of landslide occurrence can be expressed by the corrected landslide susceptibility index. Therefore, there are three variables in the PLEWM, namely, R₁, R_2, and the landslide spatial probability. In this part, a three-dimensional Bayesian formula was employed to construct the PLEWM, which is expressed as follows in Equation (7):

$$P\left(L{\left|S,R_1,R\right.}_2\right)=\frac{P\left(S,R_1,R_2\left|L\right.\right)\cdot P\left(L\right)}{P\left(S,R_1,R_2\right)}$$

(7)

where $P\left(S,R_1,R_2\right)$ is the marginal probability of the combination of rainfall variables R₁-R₂ and the landslide spatial probability S, irrespective of the landslide occurrence; considering that in natural situations, rainfall events and the spatial probability of landslide occurrence are independent, $P\left(S,R_1,R_2\right)$ can be written as $P\left(S\right)\cdot P\left(R_1,R_2\right)$. If given the landslide events, rainfall variables were also independent of the landslide spatial probability, $P\left(S,R_1,R_2\left|L\right.\right)$ can be written as Equation (8):

$$P\left(S,R_1,R_2\left|L\right.\right)=P\left(S\left|L\right.\right)\cdot P\left(R_1,R_2\left|L\right.\right)$$

(8)

and Equation (7) can be written as follows:

$$P\left(L{\left|S,R_1,R\right.}_2\right)\propto \frac{P\left(S\left|L\right.\right)\cdot P\left(R_1,R_2\left|L\right.\right)}{P\left(S\right)\cdot P\left(R_1,R_2\right)}\propto P\left(L\left|S\right.\right)\cdot P\left(L{\left|R_1,R\right.}_2\right)$$

(9)

where $P\left(L\left|S\right.\right)$ is the spatial probability. In the existing research, Equation (8) holds by default, and Equation (9) is used to forecast landslides [13,77,78]. However, the following situation may exist: in areas with high spatial probability of landslide occurrence, slopes are less stable and light rainfall may trigger landslides, while in areas with low spatial probability, slopes are very stable and heavy rainfall is necessary to trigger landslides. In this situation, the landslide spatial probability is not independent of the rainfall variables, and Equation (7) cannot be decomposed into Equation (9). Therefore, before using Equation (9) to define the LEWM, the independence of the spatial probability and temporal probability must be tested.

Furthermore, we still need a conversion criterion to accept the spatial-temporal probability and issue appropriate levels of warning information [6]. The process is typically expert-driven and may affect the results comparison. Therefore, we take the frequency of historical landslides as the criterion and uniformly determine the conversion method for the DLEWM and PLEWM to warning levels.

Table 1. The criteria of the response measures that need to be taken, the investment per unit area against the response measure, and the resulting loss per landslide.

Warning Level	Corresponding Measures	Investment (10,000 CNY/A)	Loss (10,000 CNY/Landslide)
4th-level	There is a great probability of landslide occurrence in the warning zone, and the HRS should be continuously monitored, emergency evacuation routes should be prepared, and residents should be evacuated near the HRS. The dangerous road potentially affected by landslides should be closed.	8	0.1
3rd-level	Landslides occur with a large probability, and HRS monitoring should be strengthened, residents should be notified to prepare for emergency evacuation, and warning signs should be placed on dangerous roads potentially affected by landslides.	4	2
2nd-level	Regular monitoring of the HRS, and residents should be notified to report to staff if there is any abnormality in the HRS.	2	4
1st-level	No preventive measures are needed	0	8

shows the recommended measures at each warning level, and the detailed conversion process can be found in Section 3.3.2.

2.5.2. Discriminant Matrix-Based Landslide Early Warning Model

A discriminant matrix was used to combine the rainfall thresholds, categorical landslide spatial probability, and the recommended warning levels. The construction of a discriminant matrix is generally based on expert experience, so it is also called a heuristics matrix. Similarly, to objectively compare the performances of the PLEWM and DLEWM, we used the same method mentioned in Section 2.5.1 to determine the elements in the discriminant matrix.

2.6. Early Warning Model Performance Evaluation

During the operation process of the LEWM, it is necessary to take reasonable action or corresponding measures (such as continuous monitoring, avoidance, personnel transfer, etc.) for the issued warning information. The larger the proportion of high-level early warning zones is, the greater the economic investment that is needed, and the corresponding measures will mitigate the loss caused by landslides. Therefore, the investment and the loss within one warning event represent the cost and effectiveness of the LEWM, respectively. The lower investment and loss of the LEWM suggests a higher performance. Table 1 shows the corresponding response measures at each warning level and the investment per unit area (referring to the area of the study region, abbreviated as A in the table) against the measures and the estimated loss per landslide. Notably, the investment and loss criteria in this research are defined by the authors through an exponential function. The basic idea or assumption is that as the warning level increases, the investment will increase exponentially by a factor of 2, while the loss caused by each landslide decreases exponentially by a factor of 2. As the warning level increases, the expected investment in Table 1 is set to 0, 2, 4, 8, respectively, and considering that landslides certainly will cause losses when they occur, we think that the losses should be a smaller value instead of 0 when landslides occur in the 4th-level warning area, so the expected losses for each landslide are set to 8, 4, 2, 0.1 in Table 1, respectively. In particular, it should be emphasized that this criterion is not based on any official standard, but is set by the authors through intuition (it should be following an exponential function) in order to instantiate the idea of using investments and losses to evaluate the performance of multilevel warning models, which in practice should vary from place to place.

3. Results

3.1. Landslide Spatial Probability in the Qin-Ba Mountain Area

3.1.1. Geoenvironmental Factors and Their Relationship with Landslides

As mentioned in Section 2.2, we considered 12 landslide-related geoenvironmental factors (Figure 3), and the frequency ratio model was employed to reassign the values to subclasses of the factors. In addition, we also analyzed the relationship between raw values of continuous factors and the frequency ratio values (Figure 4a–g). There is a specific pattern between the value of each factor and the spatial probability of landslides. For example, the spatial probability of landslides is exponentially related to elevation (red fitted curve in Figure 4a), indicating that landslides are more likely to occur at lower elevations. The spatial probability shows a segmented linear relationship with the slope (Figure 4b), with the highest probability at a slope of 10°. Other factors, such as the distance to rivers (Figure 4c), fault density (Figure 4f), and distance to roads (Figure 4g), are also segmentally linearly related to the spatial probability of landslides. Although there is a specific relationship between each factor and the spatial probability of landslides, the presence of uncertainty in it cannot be ignored, as it indicates the correlation between the factor and landslide occurrence. Factors with a smaller degree of uncertainty indicate a clearer relationship with landslide occurrence and a stronger correlation. The frequency ratio value of the discrete factors can be seen in Figure 4h–l. All reassigned factors were taken as input variables for each machine learning model.

3.1.2. The Spatial Probability of Landslide Occurrence

The training process was carried out in Scikit-learn 1.0.1 [79]. A grid search method was employed to find the optimal hyperparameters in the model. Subsequently, to avoid overfitting, a 10-fold cross-validation method was used to train the machine learning model, and the results are shown in Figure 5.

A test set was employed to verify the accuracy of the three machine learning models, and the ROC curves are presented in Figure S2. The area under the ROC curve indicates that the ANN model achieved an optimal performance with a value of 0.809; however, it was not much higher than that of the other two models because the area under the ROC curve of the SVM was 0.802 and that of the LR was 0.798. Simultaneously, the ROC curves of the three models are staggered, and none of the curves completely cover the other one. Therefore, these three models have a similar performance.

The next step is to verify whether these three obtained landslide susceptibility maps can express landslide spatial probability. First, the frequencies of LSI in the whole study area and landslides in the test set were counted separately (beige and blue boxes in Figure 6). Then, Equation (1) was adopted to calculate the frequency ratio values for each interval (black dots in Figure 6). The Logistic Regression-based LSI has a roughly linear relationship with the frequency ration value (Figure 6a), but exhibits a large degree of uncertainty. The R-square is only 0.6856, indicating that although the LSI can be used to express the spatial probability of landslides, it deviates significantly from the actual spatial probability in most areas. Figure 6b shows that the SVM-based LSI has a segmented linear relationship with FR and the R-square value is higher than that of the LR model, which indicates that using the SVM-based LSI to express the spatial probability of landslides will have less uncertainty. The ANN-based LSI has an exponential relationship with the frequency ratio value (Figure 6c) and has the highest R-squared value among these three models. Meanwhile, after frequency analysis, we confirm that there is indeed a different nonlinear relationship between the LSI generated by different machine learning models and the spatial probability of landslide, and in practice applications, the LSIs need to be corrected so that they show a monotonically increasing linear relationship with the spatial probability of landslides.

Since the LSI generated by ANN show a least uncertainty with the spatial probability of landslide, so we choose the ANN-based landslide susceptibility map for the correction. The procedure is simple, and we substitute all LSIs into the fitted exponential function (red fitted line in Figure 6c) and perform normalization so that we obtain corrected susceptibility map. The corrected susceptibility map and the statistical information of the LSI are shown in Figure 7. The frequency ratio value calculated using the test data showed a linear relationship with the LSI, with an R-square of 0.997. Simultaneously, the frequency ratio value calculated using validation data (landslide data from 2001 to 2015) also showed a linear relationship with the LSI (Figure 7b), which indicates that the corrected LSI can express the landslide spatial probability in a right way.

The susceptibility map in Figure 7a was used as the landslide spatial probability model in the PLEWM. During the process of developing the DLEWM, the susceptibility map was divided into very low, low, moderate, high, and very high susceptibility levels (S1–S5) with LSIs of 0–0.2, 0.2–0.4, 0.4–0.6, 0.6–0.8, and 0.8–1, respectively.

3.2. Rainfall-Induced Landslide Temporal Probability Model

3.2.1. Selecting the Optimal Combination of Rainfall Variables

Five variables were used to characterize the rainfall events, namely, cumulative effective rainfall (EE), rainfall duration (D), effective rainfall intensity (I = EE/D), antecedent rainfall in 5 days (AE5) and 10 days (AE10). Considering the correlations between rainfall variables, we list 9 possible combinations in

Table 2. Possible combinations of rainfall variables.

Variable Combination	D	I	EE	AE5	AE10
D	-	I-D	D-EE	D-AE5	D-AE10
I	-	-	I-EE	I-AE5	I-AE10
EE	-	-	-	EE-AE5	EE-AE10

. Except for D-AE5 and D-AE10, there was a parameter K in the other combinations; thus, we picked up 50 parameters evenly from 0.1 to 1 and determined the optimal parameter based on sensitivity analysis.

The results are shown in Figure 8. It can be seen that EE-D has a higher sensitivity than the other combinations, and significantly higher than I-D combination, indicating that E-D combination is more related to landslide occurrence. This conclusion is identical to that of Gariano et al. [80], they find EE-D combination have the best performance compare with other rainfall combinations like I-D, I-EE, and they suggest that it may be due to the fact that the two variables EE-D are measured independently, while there is an interdependence between other combinations resulting in. However, further, we noticed that the sensitivity of all variable combinations with EE participation are very high, like EE-I, EE-AE5, and EE-AE10 (Figure 8). This may imply that the reason for the high sensitivity of EE-D is due to the highest correlation of EE with landslide occurrence not the cause of dependence between variables as suggested by Gariano et al. [80]. When K is in the interval from 0.78 to 0.84, the variable combinations involving EE will have the highest sensitivity; thus, in this research, we chose EE-D as the optimal variable combination and K = 0.816 as the optimal parameter.

3.2.2. Empirical Rainfall Threshold

To define the empirical rainfall threshold, we use the triggering rainfall events from 2001 to 2015 and adopt the method proposed by Brunetti et al. [20] to determine the parameters in Equation (6). Here, we set four equally spaced exceedance probabilities: 20% (T_{20, S}), 40% (T_{40, S}), 60% (T_{60, S}), and 80% (T_{80, S}). The rainfall condition can be divided into five hazard levels from T1 to T5, and the equation of the threshold is shown in Figure 9.

3.2.3. Temporal Probability Model of Landslide Occurrence

To develop a temporal probability model of landslide occurrence, we performed frequency statistics on triggering and nontriggering rainfall events in log (D)-log (EE) space. The Bayesian formula (Equation (5) was employed to calculate the posterior probability P(L|D, EE) of landslide occurrence. The results are displayed as a bar chart in 3D space (Figure 10). Under different values of D, the posterior probability of landslide occurrence always increases exponentially with increasing EE (red fit line in Figure 10), while under different values of EE, the posterior probability does not show a significant increase or decrease with increasing D. This finding indicates that EE plays a leading role in the process of triggering landslides. Therefore, it also explains why the sensitivity index of the rainfall variable combination with EE participation is much higher than that of other combinations (Figure 8).

Based on the histogram in Figure 10, the ANN model was employed to fit the posterior probability of landslide occurrence to obtain the temporal probability model of landslides. Since the number of rainfall event samples in each grid shown in Figure 10 is very different, and generally, the number of rainfall event samples will decrease with the increase in EE and D, the posterior probabilities calculated in each grid will have a different degree of reliability, which is determined by the number of samples. Therefore, in the process of fitting, we used the number of samples in the grids to weight the corresponding posterior probability, and the fit result is shown in Figure 11.

The contours in Figure 11 can be used to represent the thresholds derived from the temporal probability model, and the values indicate the probability of landslides triggered by rainfall conditions (P (L|D, EE)). Meanwhile, the slope of the threshold gradually decreases as P (L|D, EE) increases, and when it reaches 0.05, the slope is approximately zero. This phenomenon means that under different rainfall conditions, the contributions of EE and D in triggering landslides gradually change, and when P (L|D, EE) > 0.05, the triggering process is nearly dominated by EE. This point is very different from the thresholds defined by power-law function (Figure 9) because it assumes that the contributions of EE and D do not change, and each threshold curve is parallel to each other (Figure 9). Such assumptions may not be suitable for some complex situations.

3.2.4. Verification of the Temporal Probability Model

After we developed the landslide temporal probability model, the accuracy and the uncertainty of the model were verified using rainfall data from 2001 to 2015 (training data) and 2016 to 2020 (test data), respectively. First, we inputted the EE and D of all rainfall events for each of the two time periods into the temporal probability model and calculated the temporal probability of each rainfall event. Subsequently, we summarized the number of recorded triggering rainfalls and model-predicted triggering rainfalls by month (Figure 12a). It can be seen that the trend between the number of predicted triggering rainfall and the recorded triggering rainfall is consistent. During part of the training data, the correlation coefficient reached 0.886, and in part of the test data, the correlation coefficient was 0.845, which indicates that the model prediction of the temporal probability of rainfall-induced landslides is approximately the same as the actual situation, and the model does not show overfitting. However, the model also has some uncertainties due to various factors, such as insufficient temporal resolution and precision of the recorded rainfall data and incomplete rainfall-induced landslide inventory. The bias between the number of predicted triggering rainfall events and actually recorded events were plot in Figure 12b. It can be seen that the largest bias reaches −17.8, but the mean value is approximately zero.

To investigate the uncertainty in the temporal model, we performed statistical analysis of the bias in the training data (Figure 13a) and test data (Figure 13b). For the training set, the mean value of the bias is 0.049, which means that the model produces an average bias of 0.049 over a month, with a daily average bias of approximately 0.0016. Considering that there are 42 weather stations in the study area, the bias in prediction for each weather station is 0.0016/42, so we can assume that the model does not produce a significant bias. At the same time, $\sigma =1.975$, which means that the predicted bias of 95% of the data is between $\pm 2\sigma$; in other words, most of the errors generated by the forecasting model each month are concentrated between 3.94 and −3.94. If we assume that the bias is cumulative, averaged over each weather station, the model produces a bias between $\pm 0.003$ per day. The same situation exists for the test data set. Considering that the model, although it does have some random bias, does not show systematic bias, we still use it as a temporal probability model to predict the temporal probability of landslide occurrence.

3.3. Warning Strategy

3.3.1. Independence Test

According to Equations (8) and (9), we know that the independence test for the spatial and temporal probability model was equivalent to testing whether Equation (8) holds. Thus, we took the spatial probability as the x-axis and the temporal probability as the y-axis for the 2D histogram statistics. Among them, spatial probability was equally divided into 10 intervals, and temporal probability was divided into 11 intervals according to Figure 11. A total of 110 grids were obtained. Subsequently, we counted the frequency of landslide occurrence in each grid and calculated the product of the marginal probabilities ($P\left(S\left|L\right.\right)\cdot P\left(D,EE\left|L\right.\right)$) and the joint probability ($P\left(S,D,EE\left|L\right.\right)$) within each grid (Figure 10). The result can be seen in Figure S3. The joint probability and the product of the marginal probability satisfied a linear relationship, the slope of the fitted line was close to one, and the correlation coefficient reached 0.9827. Therefore, we can argue that $P\left(S\left|L\right.\right)\cdot P\left(D,EE\left|L\right.\right)$ and $P\left(S,D,EE\left|L\right.\right)$ were statistically consistent, and the spatial and temporal probabilities of landslides in the study area were independent.

According to Equation (9), we defined a spatiotemporal probability index (P_st) to forecast landslides. P_st can be calculated by the following equation and is proportional to the spatiotemporal probability of the landslide.

$$P_{st}=P\left(L\left|S\right.\right)\cdot P\left(L\left|D,EE\right.\right)$$

(10)

3.3.2. Conversion of the Landslide Forecast Model to Warning Levels

For the PLEWM, we first calculated the P_st of all landslides from 2001 to 2015 by Equation (10) and counted the frequency of landslides against different intervals of P_st (Figure 14b). For the DLEWM, we quantified the landslide susceptibility levels S1 to S5 and rainfall hazard levels T1 to T5 as 1, 2, 3, 4, and 5, respectively. The five numerically susceptible levels multiplied by five rainfall hazard levels provided 14 discrete spatiotemporal warning indexes. Subsequently, we counted the frequency of historical landslides in each spatiotemporal warning index, and the results are shown in Figure 14a.

We used the criterion that the historical landslide frequencies for each warning level are as equal as possible to define the discriminant matrix and link the P_st with the warning levels. As a result, for the DLEWM, spatiotemporal warning indexes 1, 2, and 3 are classified as 1st-level warning, including 21.7% of the historical landslides; 4 and 5 are classified as 2nd-level warning, including 28.83% of the historical landslides; 6, 8, 9, 10, and 12 are classified as 3rd-level warning, including 29.63% of landslides; and 15, 16, 20, 25 classified as 4th-level warning, including 19.80% of landslides, and the corresponding discriminant matrix is shown in

Table 3. Discriminant matrix determined by the frequency method.

Warning Levels	S1 (0–0.2)	S2 (0.2–0.4)	S3 (0.4–0.6)	S4 (0.6–0.8)	S5 (0.8–1.0)
T1 (<20%)	1st-level warning	1st-level warning	1st-level warning	2nd-level warning	2nd-level warning
T2 (20–40%)	1st-level warning	2nd-level warning	3rd-level warning	3rd-level warning	3rd-level warning
T3 (40–60%)	1st-level warning	3rd-level warning	3rd-level warning	3rd-level warning	4th-level warning
T4 (60–80%)	2nd-level warning	3rd-level warning	3rd-level warning	4th-level warning	4th-level warning
T5 (>80%)	2nd-level warning	3rd-level warning	4th-level warning	4th-level warning	4th-level warning

. For the PLEWM, P_st ≤ 0.004 corresponds to the 1st-level warning, 0.004 < P_st ≤ 0.02 corresponds to the 2nd-level warning, 0.02 < P_st ≤ 0.064 corresponds to the 3rd-level warning and 0.064 < P_st corresponds to the 4th-level warning. The black dashed lines in Figure 14 are the dividing lines of the warning levels, which ensures the same number of landslides in the same-level warning zone in the DLEWM and PLEWM; in other words, the landslide forecast models in the DLEWM and PLEWM were converted to warning levels according to the same criterion and equal correct warning rates are set in advance for both models.

3.4. Simulated Warning Using the PLEWM and DLEWM

We use the PLEWM and DLEWM to conduct simulated warnings in the rainy season (July to September) each day from 2016 to 2020, and we counted the area of different warning zones issued by the two warning models and the number of landslides occurring within each warning zone. The results are summarized monthly (Table S1) and shown in Figure 15 and Figure 16.

From Figure 15, we can see that in each month, the area of warning zones that need to take corresponding response measures (2nd-, 3rd-, and 4th- level warnings) issued by the PLEWM is significantly smaller than that issued by the DLEWM. The total area of warning zones of the 2nd, 3rd, and 4th levels issued from 2016 to 2020 only accounts for 60.23%, 45.99%, and 43.98% of the same-level warning zone issues by the DLEWM, respectively. During the process of converting the landslide forecast model to warning levels, we used the criterion that the same number of historical landslides will occur within the same-level warning zone for both the DLEWM and PLEWM, which means that we set the same correct warning rate for both the DLEWM and PLEWM in advance. This means that PLEWM will achieve the same warning effect as DLEWM by issuing a smaller warning area and is more spatially accurate.

As mentioned in Section 3.3.2, we preset the correct warning rate to be the same for both warning models, so the simulated warning results show that the number of landslides occurring within the same level of warning zone is approximately the same in the PLEWM and DLEWM (Figure 16). However, in July 2018, there was an anomaly, the number of landslides occurring in 1st-, 2nd-, 3rd-, and 4th-level warning zones are 88, 57, 84, and 6, respectively (37.45%, 24.26%, 35.74%, and 2.55% of the total), a large number of landslides occurred within the low-level warning zone, and only 6 landslides were located in the 4th-level warning zone (Figure 16d). This indicates that the correct warning rate of DLEWM during this period is below the expected level, because when we defined the discriminant matrix in Section 3.3.2, the expected number occurred in each level of warning zone should be approximately. In contrast, for the PLEWM, the number of landslides occurring in 1st-, 2nd-, 3rd-, and 4th-level warning zones is 13, 42, 64, 116, respectively (5.54%, 17.87%, 27.23%, 49.36% of the total), almost half of the landslide occurred within the 4th-level warning zone, only 13 landslides occurred within the 1st-level warning zone, indicating that the model’s correct warning rate far exceeds the expected level.

To find out the cause of this anomaly, we examined information on landslides and the rainfall data during this period, and we found that the landslides that occurred in July 2018 were mainly concentrated in LveYang County on the 11th and 14th, with a total 182 landslides occurring on these two days. The rainfall features of LveYang County in these two days were D = 16 d and EE = 143.02 mm, and D = 19 d and EE = 145.22 mm, respectively. These were extreme rainfall events of long duration and high intensity. The rainfall threshold defined based on power-law function showed that the rainfall events in these two days was high hazard (T4) and medium hazard (T3), respectively, not reaching the very high hazard (T5) level. However, according to the temporal probability model (Figure 11), the temporal probability of rainfall on these two days was 0.5855 and 0.6422, respectively, ranking in the first 1% of all triggering rainfall events in 2001–2015. This finding indicates that there is a significant difference between the power-law-based rainfall threshold and temporal probability model when judging the hazard degree of prolonged and intense rainfall. The reason is that the power-law-based rainfall threshold curves are parallel to each other, so it is believed that EE and D always have the same contribution in the process of triggering landslides; therefore, the hazard degree will be underestimated when the rainfall duration is long, while the temporal probability model shows that EE has a greater contribution in the process of triggering landslides when the rainfall duration is long, and it makes a correct estimate in this warning process.

Figure 17 shows the location of the warning zones at each level issued by the PLEWM and DLEWM on July 11 and 14 in the study area and the location of the landslide occurrence. On 11 July, the 3rd-level and 4th-level warning zones issued by the PLEWM were mainly concentrated in Lveyang and NingQiang Counties (Figure 17a), while 83 landslides occurred in LveYang on that day, of which 75 (90.36%) were located in the 3rd-level and 4th-level warning zones, and the rest were located in the 2nd-level warning zone. On the same day, the DLEWM mainly issued a 2nd-level warning zone in LveYang, where the areas of the 3rd- and 4th-level warning zones accounted for only 2.49% (Figure 17c), resulting in 45 landslides (54.22%) located in the 2nd-level warning zones and just 6 landslides (7.23%) located in the 4th-level warning zone, which did not achieve the best warning effectiveness. Similarly, on July 14, because the rainfall threshold underestimated the hazard degree of the rainfall in LveYang, the DLEWM only issued 10.3% of the 3rd-level warning zones (Figure 17d), and the rest were 1st-level warning zones. On that day, 99 landslides occurred in LveYang County, of which 64 landslides (64.65%) were located in the 1st-level warning zone, and just 35 landslides were located in the 3rd-level warning zone, caused a significant missed warning. However, for the PLEWM, 8, 29, and 62 landslides occurred in the 2nd-, 3rd- and 4th-level warning zones, respectively, and no landslides occurred in the 1st-level warning zone (Figure 17b, Table S2), indicating a very high correct warning rate.

As seen, since the power-law based rainfall threshold model is a relatively simple universal model, and it may not be able to make a reasonable judgment about the hazard of some complex meteorological conditions, such as extreme rainfall with long duration and high intensity. It is here that extreme meteorological events can often trigger a large number of geological disasters, thus causing serious casualties and property damage. The temporal probability model, on the other hand, is relatively more flexible and can fully exploit the relationship between rainfall and landslide occurrence in historical data, so a more reasonable estimate of the hazard of some extreme rainfall can be made.

Finally, since the early warning model we study can issue multiple levels of warning information, it is difficult and not intuitive to evaluate the performance of the model using the area of the warning zones issued at each level or the number of landslides occurring within the warning zones at each level. Therefore, we propose to use Investment and Loss (Table 1) to quantify the difference of each warning level, the area of warning zones at each level, and the number of landslides occurring within the warning zones at each level, and these two indicators are more intuitive. We calculated the expected investment required by the two warning models for issuing early warning information and the losses caused by landslides, based on the criteria in Table 1. The result is shown in Figure 18.

It can be seen that the expected investment required is almost half of the DLEWM in all months (Figure 18a) because of the high spatial accuracy of the PLEWM, and with a total expected investment of 54.54% of the DLEWM. In this sense, we can say that the spatial accuracy of the PLEWM is 1.834 (1/0.5454) times higher than that of the DLEWM. We can also see that the losses caused by landslides are roughly the same with the help of the two warning models (Figure 18b), except for July 2018. In July 2018, because the DLEWM underestimated the hazard of two extreme rainfall events, the expected loss caused by landslides was 11.06 million CNY, while with the help of early warning information issued by the PLEWM, the expected loss caused by landslides was 4.12 million CNY, a reduction of62.75% in comparison (Figure 18b).

Based on the above analysis, the PLEWM has the following characteristics compared with the conventional DLEWM. First, the model is highly spatially accurate. On the premise of achieving the same correct warning rate, the areas of the 2nd-, 3rd-, and 4th-level warning zones are generally smaller, accounting for 60.23%, 45.99%, and 43.98% of the same-level warning zones in the DLEWM, respectively, and the investment required to issue warning information, according to the criteria of Table 1, accounts for only 54.54% of the DLEWM. Second, for extreme rainfall events with long durations and high intensities, the correct warning rate of the PLEWM is far higher than that of the conventional DLEWM, which can greatly reduce property losses.

4. Discussion

4.1. Spatial Probability Model of Landslide Occurrence

In most previous studies, the spatial probability of landslide occurrence was directly represented by the LSI and without verification or correction. In this research, we first used training data to test the relationship between the LSI and frequency ratio values, and the functional relationship between the LSI and frequency ratio values was determined by curve fitting. Then, we corrected the susceptibility map through the fitting functions and ensured that the corrected LSI had a linear increasing relationship with the spatial probability of landslide. Whether the two processing methods will have a significant impact on the performance of the LEWM has not been studied in detail, but from the perspective of probability, what is required in the LEWM is the spatial probability of landslides rather than the landslide susceptibility, and this correction process is necessary.

4.2. The PLEWM or DLEWM in Practical Applications?

The analysis in Section 3.4 shows that the DLEWM performance is relatively poor. Under the condition of setting the same correct warning rate, the issued area of the 2nd-, 3rd-, and 4th-level warning zones is approximately twice that of the PLEWM. Simultaneously, due to the defects of the power-law-based rainfall threshold itself, the DLEWM may not be able to issue the appropriate warning information for some extreme rainfall weather events, resulting in serious consequences. However, as a knowledge-driven model, the DLEWM is widely used in practice due to its simple form and clear physical meaning. With the support of expert experience, the performance will be guaranteed.

As a data-driven model, the PLEWM has a reasonable probability theoretical basis and data to support a series of steps, such as spatial and temporal probability modeling and the coupling of spatial-temporal probability. Therefore, the model has high accuracy and credibility, and can also provide reasonable estimates of temporal probability for some extreme rainfall events. More importantly, probabilistic approaches are commonly used in quantitative risk assessment, which provides further application of PLEWM in regional landslide risk assessment. However, its limitations are also evident in that it requires a good understanding of the geoenvironmental conditions in the study area, as well as landslide and monitoring data over a long-time span. Therefore, it is suitable for areas with adequate research.

According to comprehensive analysis, for the regions in the early stage of early warning, due to the lack of systematic and comprehensive disaster, meteorological, and hydrological monitoring data, DLEWM can provide a basic service of geo-hazard early warning. For regions with long time of geo-hazard early warning foundation, and need to further achieve more accurate warning effective, PLEWM is a scientific and reasonable alternative.

4.3. Performance Evaluation of the LEWM

The performance evaluation of multilevel early warning models is an important and complex task, and many important and unique aspects need to be considered, such as how to quantify the differences between warning levels, the relationship between the number of landslides, the performance of warning models, and the relationship between the area of each warning zone and the performance. LEWSs are often regarded as cost-effective risk mitigation measures. The economy or capital can perfectly quantify all the above aspects. Therefore, it is a new way to evaluate LEWM performance from an economic perspective.

Since the high-level warning zone issued by the PLEWM is smaller than that of the DLEWM, the number of landslides occurring in the high-level warning zone is also more than that of the DLEWM. Therefore, the performance comparison between the PLEWM and DLEWM can be completed by only two indicators of Investment and Loss. However, in practical applications, there may be more complex situations. For example, warning model A has more investment and less loss caused by landslides, while warning model B has less investment and greater loss caused by landslides. In this situation, the performance of the two models is not absolutely superior or inferior, but rather their applicability or inapplicability to the specific situation. Therefore, there is a need for more detail indicators to describe each aspect of the warning model. We have initially defined a series of indicators in economic perspective, and applicated them in the simulate warning result of PLEWM and DPLEWM (

Table 4. Economic indicators system for evaluating the performance of multi-level warning models.

Symbol	Performance Indicator	DLEWM	PLEWM	Description or Formula
LP	Potential losses	3088	3088	losses caused by landslides if no warning information is given, which is calculated based on Table 1.
L	Loss	1887.9	1115.1	Losses caused by landslides with the help of warning information
Inv	Investment	84.16	45.90	Cost inputs required for the response measures
Eff	Effectiveness	1200.1	1972.9	Mitigated losses, LP − L
ER	Effectiveness rate	0.3886	0.6389	E/LP
ECW	Effectiveness of correct warning	844.8	1558	Losses mitigated in case of correct warning
ERCA	Effectiveness rate of correct warning	0.2736	0.5045	ECW/LP
LE	Losses caused by error warning	1388.7	780.5	Losses caused by missed and false warning
LER	Loss rate of error waning	0.4497	0.2528	LE/LP
CE	Total cost-effectiveness	35.85	151.39	λ·E − I, here λ is taken as 0.1
CER	Cost-effective conversion rates	1.4257	4.2983	λ·E/I, here λ is taken as 0.1
CL	Total costs and losses	272.97	157.41	λ·Loss + I, here λ is taken as 0.1

).

In Table 4, the indicator of “Potential losses” means the losses caused by landslides without the help of LEWM. Inthis situation, all landslides will occur in 1st-level warning zone. The indicator of “effectiveness” refers to the loss mitigated with the help of the early warning model, which should be equal to the “Potential loss” minus the actual loss caused. The “effectiveness rate” equals “effectiveness” divided by “potential losses”. The indicator of “effectiveness of correct warning” refers to the loss of mitigation under the condition of correct warning, and here we refer to the research of Calvello and Piciullo [7] to define the concept of correct warning, missed warning and false warning for multi-level early warning mode, with detailed criteria shown in Tables S3 and S4. The “Losses caused by error warning” refer to the loss caused by landslide under the condition of false warning and missed warning. Similarly, indicators of “Effectiveness rate of correct warning” and “Loss rate of error waning” have the same means. “Total cost-effectiveness”, “Cost-effective conversion rates”, and “Total costs and losses” are three comprehensive indicators, which can be used when considering the investment and loss together, where the parameter λ reflects the weighting of the cost and effectiveness in the model performance.

With the help of these indicators, decision makers can get a more comprehensive understanding of an early warning model’s performance and select an early warning model from a variety of perspectives.

5. Conclusions

The primary objective of this research was to develop a probabilistic-based LEWM at the regional scale and compare it with the conventional DLEWM. Three machine learning models, namely, the Logistic Regression, Support Vector Machine, and Artificial Neural Network, combined with the frequency ratio model, were first used to obtain the landslide spatial probability model. Then, based on a detailed inventory of rainfall-induced landslides that occurred from 2001 to 2020, we implemented a Bayesian approach considering the variable EE-D to determine the temporal probability of rainfall-induced landslides, and the Bayesian formula was used to couple spatial and temporal models to develop a PLEWM. Finally, we simulate warnings day-by-day in the rainy season (July–September) from 2016 to 2020 using the PLEWM and DLEWM, separately, and economic indicators were used to evaluate the model performance. The results are described as follows:

The Logistic Regression, Support Vector Machine, and Artificial Neural Network have similar prediction results for landslide susceptibility in the study area, but the Artificial Neural Network has the least uncertainty. In addition, the LSIs calculated by the Logistic Regression, Support Vector Machine, and Artificial Neural Network have linear, piecewise linear, and exponential relationships with the spatial probability, respectively. Thus, the landslide susceptibility generated by the machine learning models may sometimes be different from the landslide spatial probability.

The sensitivity analysis shows that the optimal rainfall variable combination is EE-D, and the optimal coefficient K in EE is 0.816. In addition, the contributions of EE and D to inducing landslides gradually change; therefore, the threshold curve determined by the temporal probability model is not a parallel straight line in the logD-logEE space. When P (L|D, EE) > 0.05, the triggering process is nearly dominated by EE, and this point is very different from the power-law-based rainfall threshold.

The simulated warning results show that the areas of the 2nd-, 3rd-, and 4th-level warning zones issued by the PLEWM account for 60.23%, 45.99%, and 43.98% of that of the DLEWM, respectively. The investment in issue warning information and the losses caused by landslides account for 54.54% and 59.06% of the DLEWM, respectively. Moreover, under extreme rainfall conditions, the correct warning rate of the PLEWM is much higher than that of the DLEWM.