Framework for risk assessment of economic loss from structures damaged by rainfall-induced landslides using machine learning

ABSTRACT

Given the increased frequency of extreme rainfall events, pre-disaster countermeasures against landslides triggered by heavy rainfall are important to enhance disaster resilience. This study presents a methodology for economic risk assessment of structures affected by rainfall-induced landslides using machine learning (ML). Random Forest and LightGBM algorithms were applied to develop ML-based landslide prediction models considering the spatial distributions of landslide conditioning and triggering factors. The rainfall index was calculated considering the temporal variation in rainfall and was used as a feature associated with rainfall intensity. The rainfall hazard curve, representing the relationship between the rainfall index and its annual exceedance probability, was statistically estimated using a generalised extreme value distribution. Rainfall-induced landslide susceptibility was assessed using an ML-based landslide prediction model and rainfall hazard curve. Finally, the risk curve associated with the economic loss from structures damaged by rainfall-induced landslides was estimated based on landslide susceptibility and structure distribution maps. In this study, LightGBM showed better prediction performance for evaluating rainfall-induced landslide susceptibility than Random Forest. An illustrative example is presented to demonstrate that the proposed methodology can be used to develop an appropriate risk-based disaster mitigation strategy.

KEYWORDS:

• Machine learning
• landslide susceptibility
• rainfall hazard
• economic risk
• disaster mitigation strategy

1. Introduction

Rainfall-induced landslides are one of the most dominant hazards to communities in the mountainous areas of Japan because of high precipitation and severe topography (Yoshimatsu and Abe 2006). Unprecedented heavy rainfall occurred over a wide area of Japan between 28 June and 8 July 2018. This rainfall event, which is officially referred to as the Heavy Rain Event of July 2018 (Japan Meteorological Agency (JMA) 2018), not only caused severe damage to structures and infrastructure, but also resulted in significant loss of life (Tsuguti et al. 2019). Typhoon Hagibis struck Japan in October 2019, resulting in serious property damage, deterioration of the living environment, and economic loss due to floods and landslides (Wang et al. 2023). Landslides triggered by rainfall can cause significant delays in the post-disaster recovery processes. Considering the long-term trend of the increase in the intensity and frequency of rainfall events (Sun et al. 2023), countermeasures against heavy rainfall and associated landslides are critical for minimising the negative social impacts and enhancing the disaster resilience of communities in potentially affected areas.

The time and location of landslides due to rainfall can be predicted using process-based models based on slope stability analysis, considering various physical and mechanical properties (Kim et al. 2021). The JMA has provided a real-time rainfall-induced landslide hazard map estimated based on the soil water index calculated using a three-layer tank model with fixed parameters throughout Japan (Matsuyama, Saito, and Zemtsov 2021). However, it is difficult to apply these models over a wide area since a large dataset with detailed information on hydrological, lithological and soil characteristics is required to consider the local variation in landslide conditioning factors, such as topography and vegetation (Guzzetti et al. 2007). Data-driven approaches based on machine learning (ML) have been studied to develop landslide prediction models applicable to large and complex areas (Huang and Zhao 2018; Kainthura and Sharma 2022; Mondal and Mandal 2018; Ng et al. 2021; Taalab, Cheng, and Zhang 2018). The spatial distributions of various landslide conditioning and triggering factors can be correlated with landslide susceptibility using ML algorithms, considering the complex nonlinear problems of landslide occurrence. These ML-based landslide prediction models are expected to support decision-making in landslide disaster management over wide areas.

A hazard curve representing the relationship between hazard intensity and occurrence probability should be considered when evaluating the safety of structures against future disasters (Ishibashi et al. 2021a). In the field of earthquake engineering, probabilistic seismic hazard analysis has been performed with the development of ground motion prediction models for future earthquakes considering aleatory and epistemic uncertainties, which represent the natural randomness in the physical system and lack of knowledge in the model due to limited data, respectively (Bommer 2022). Statistical approaches based on historical data are generally used to estimate rainfall hazards since the occurrence time, location, and intensity of rainfall events depend on complex relationships among various meteorological factors. Moreover, solid technologies for evaluating the uncertainties associated with meteorological factors during extreme rainfall events have not been developed. Several methodologies have been studied on landslide susceptibility mapping based on the occurrence probability of extreme rainfall events estimated using statistical approaches (Kawagoe, Kazama, and Sarukkalige 2010; Lee et al. 2021; Mallick et al. 2021). In these previous studies, rainfall intensity measures obtained without considering the temporal variation in rainfall (i.e. both long- and short-term rainfall) before landslide occurrence were used in the prediction models. In general, deep-seated landslides occur because of long-term rainfall, whereas shallow landslides and debris flows are triggered by short-term heavy rainfall (Abraham et al. 2020). Liu and Wang (2022; 2023) investigated a critical rainfall pattern curve representing the threshold for slope stability and quantified the effects of rainfall intensity and duration on the annual occurrence probability of rainfall-induced landslides at a specific slope. It is important to develop methodologies for landslide susceptibility assessment that incorporate rainfall hazard analysis, considering the effects of temporal variation in rainfall on landslide occurrence.

Risk is a performance indicator of the social impacts of future disasters. Fu et al. (2020) proposed a methodology for the landslide risk assessment of property damage in areas with landslide data and damage records. Alqadhi et al. (2022) used ML to assess a distribution map of ecosystem services affected by rainfall-induced landslides with various return periods. Although methodologies for quantifying the impacts of landslides have been investigated, very few studies have been conducted on landslide risk assessment considering rainfall hazards, particularly for estimating the risk curve representing the relationship between the consequence of landslides and its exceedance probability. The risk curve supports both long- and short-term decision-making regarding disaster countermeasures (Ishibashi et al. 2021b). Frameworks for estimating the risk of landslides triggered by rainfall based on hazard analysis should be established to improve disaster management systems in landslide-prone areas.

This study presents a framework for the risk assessment of structures damaged by rainfall-induced landslides using ML. The rainfall index, calculated considering the effects of both long- and short-term rainfall on landslide occurrence, was used as a rainfall intensity measure. The performances of landslide prediction models developed using the Random Forest and LightGBM algorithms were compared to determine the optimal ML algorithm for landslide susceptibility assessment. Rainfall-induced landslide susceptibility was assessed using the ML-based prediction model and the rainfall hazard curve estimated based on a statistical approach. The risk curve associated with the economic loss from structures affected by rainfall-induced landslides was then estimated by combining the landslide susceptibility map and the spatial distribution of structures. In an illustrative example, the economic risk of structures damaged by rainfall-induced landslides in an area of Japan was analysed to discuss efficient disaster management.

2. Procedure for risk assessment of economic loss from structures affected by rainfall-induced landslides

Figure 1 shows the procedure for estimating the economic risk of structures damaged by rainfall-induced landslides. The proposed framework consists of three phases: (i) development of a rainfall-induced landslide prediction model using ML, (ii) landslide susceptibility assessment, and (iii) risk assessment of structures affected by rainfall-induced landslides. In phase (i), the Random Forest and LightGBM algorithms are applied to develop the landslide prediction models. The optimal prediction model is determined by comparing the prediction performance. In phase (ii), landslide susceptibility in the analysed area is assessed using the ML-based landslide prediction model and the hazard curve obtained using a statistical method. Finally, in phase (iii), the risk curve associated with the economic loss from structures damaged by rainfall-induced landslides is estimated by combining landslide susceptibility and structure distribution maps.

Figure 1. Procedure for estimating economic risk of structures damaged by rainfall-induced landslides.

2.1. Development of landslide prediction model using machine learning

2.1.1. Data collection for machine learning

The areas affected by the Heavy Rain Event of July 2018 and Typhoon Hagibis in Japan were selected as data collection areas. The ML dataset was created by combining the data obtained from these areas to ensure the applicability of the landslide prediction models to various areas with different landslide conditioning factors. The locations of the data collection areas are shown in Figure 2, which also presents the location of the area used to assess the economic risk of structures affected by rainfall-induced landslides as the illustrative example described in Section 3. The distribution maps of the ML dataset in the areas affected by the Heavy Rain Event of July 2018 and Typhoon Hagibis are shown in Figures 3 and 4, respectively. The spatial distribution data of landslide location, rainfall index, elevation, slope angle, terrain ruggedness index (TRI) (Riley, DeGloria, and Elliot 1999), general curvature, vegetation, topography, soil and lithology were collected and used as features (i.e. explanatory variables) in ML. The Geospatial Information Authority of Japan (GSI) (2018; 2019) provided data on landslide distributions identified from aerial photographs taken after the Heavy Rain Event of July 2018 and Typhoon Hagibis. The data collection areas were divided into 10 × 10 m grid cells using a geographic information system (GIS). A grid resolution of 10 m is sufficient to capture local topographic characteristics and can reduce the uncertainties associated with geological conditions in each grid cell to insignificant levels (Kim et al. 2021; Tofani et al. 2017).

Figure 2. Locations of data collection areas.

Figure 3. Distribution maps of machine learning dataset in the area affected by the Heavy Rain Event of July 2018.

Figure 4. Distribution maps of machine learning dataset in the area affected by Typhoon Hagibis.

Landslide occurrence depends on temporal variation in rainfall. In this study, considering the effects of both long- and short-term rainfall, the rainfall index proposed by Nakai et al. (2006) was used as a feature to represent the rainfall intensity. The rainfall index has a good correlation with the occurrence of rainfall-induced landslides (Moriwaki et al. 2019). The rainfall index R^’ is calculated as follows: (1) $R_{\mathrm{fw}}=\sqrt{{(R_1-R_w)}^2+a^2{(r_1-r_w)}^2}$(1) (2) $R^{\mathrm{\prime }}=R_{\mathrm{fw}0}-R_{\mathrm{fw}}$(2) where R_w is the long-term effective rainfall with a half-time value period of 72 h, r_w is the short-term effective rainfall with a half-time value period of 1.5 h, R₁ and r₁ are the reference values for long- and short-term effective rainfall, respectively, a is the weighting factor, and R_fw0 is the R_fw calculated when R_w = r_w = 0. The half-time value periods of the long- and short-term effective rainfall represent the regional characteristics of the surface and subsurface rainfall-runoff processes, respectively. Takahashi et al. (2000) studied on the relationship among landslide occurrences, regional differences and half-time value period of effective rainfall, and revealed that the long-term effective rainfall with a half-time value period of 72 h is appropriate for the granite area, which is widely distributed in the data collection areas for ML, as described below. In addition, Terada and Nakaya (2001) reported that short-term effective rainfall with a half-time value period of 1.5 h can be applied regardless of regional characteristics. R₁, r₁ and a were set to 600, 200 mm and 3, respectively. The rainfall index calculated with these values of R₁, r₁, a and half-time value periods is correlated with landslide occurrences triggered by rainfall in Japan (Moriwaki et al. 2019).

The rainfall index distributions in the data collection areas, shown in Figures 3(b) and 4(b), were computed based on the rainfall data observed by the eXtended RAdar Information Network (XRAIN) using X- and C-band multi-parameter radars provided by the Ministry of Land, Infrastructure and Transport (MLIT) (2010). Distribution maps of elevation, slope angle, TRI, and general curvature were obtained based on the digital elevation model with a 10 m resolution provided by GSI (2016). A 1/25,000 vegetation map based on field investigations and aerial photographs provided by the Biodiversity Center of Japan was processed using GIS. Distribution maps of the topography, soil and lithology were developed based on the 1/200,000 national land survey provided by the MLIT. Decomposed granite soil (i.e. Masado), which is formed by weathered granite belonging to plutonic rock and can significantly lose strength when submerged, is widely distributed in the two data collection areas (Irasawa et al. 2020; Miura and O-hara 1979). Finally, 109,703 landslide and 19,696,111 non-landslide grid data points were obtained to develop ML-based landslide prediction models. As shown in Region A of Figure 3, although the rainfall index was relatively low, many landslides occurred because of a combination of landslide-prone conditions, such as a large slope angle and plutonic rock, which are recognised as major factors in rainfall-induced landslides (Kaihara et al. 2022). The occurrence of rainfall-induced landslides depends on complex interactions between various conditioning and triggering factors.

2.1.2. Machine learning algorithm

The ML approach is suitable for handling large datasets and identifying complex landslide phenomena hidden in the data (Tehrani et al. 2022). In this study, models for solving the binary classification problem of landslide/non-landslide classes were developed using the Random Forest and LightGBM algorithms. Both ML algorithms can be used not only to develop robust classifiers but also to evaluate feature importance scores representing the contributions of each explanatory variable to the model prediction (Boutahir et al. 2022; Moorthy et al. 2020). The optimal ML algorithm for landslide susceptibility assessment was determined by comparing the classification performances of the prediction models.

2.1.2.1. Random forest

Random Forest is a bagging method based on decision trees (Breiman 2001). The final prediction result is determined by the voting results of multiple decision trees built using the bootstrap sampling method (Efron 1979).

In this study, the CART algorithm, using the Gini index as the attribute selection measure, was applied to construct a decision tree (Breiman et al. 1984). The Gini index is defined as follows: (3) $\mathrm{Gini}(v)=1-\sum _{j=1}^{C}p{(j|v)}^2$(3) where Gini(v) is the Gini index at node v, C is the number of classes ( =  2), and p(j|v) is the probability that the training sample belongs to class j at node v.

The node with the highest information gain was selected for further splitting when constructing the decision tree. The information gain is calculated as follows: (4) $\mathrm{IG}(v)=\mathrm{Gini}(v)-\{p_L\cdot \mathrm{Gini}(v_L)+p_R\cdot \mathrm{Gini}(v_R)\}$(4) where IG(v) is the information gain at node v, v_L and v_R are the left and right child nodes of node v, respectively, and p_L and p_R are the probabilities that the training sample belongs to v_L and v_R, respectively.

2.1.2.2. LightGBM

LightGBM is a tree-based learning algorithm proposed by Ke et al. (2017) that uses the gradient boosting decision trees (GBDT) algorithm (Friedman 2001). Histogram-based algorithm and leaf-wise growth strategy are applied to increase the training speed and reduce memory consumption in LightGBM (Liang et al. 2020).

For the given dataset with n samples and m features $D=\{({\mathbf{x}}_i,y_i)|i\in \{1,\dots ,n\}\}$, where ${\mathbf{x}}_i\in {\mathbb{R}}^m$ and $y_i\in \mathbb{R}$ are randomly selected from the training dataset in GBDT, the prediction result ${\hat{y}}_i$ for the i-th training sample x_i is expressed as follows: (5) ${\hat{y}}_i=\sum _{k=1}^K{f}_k({\mathbf{x}}_i)$(5) where K is the number of decision trees and f_k is a function of the k-th decision tree. The number of features (m) was nine in this study (i.e. rainfall index, elevation, slope angle, TRI, general curvature, vegetation, topography, soil and lithology), as shown in Figures 3 and 4.

The training goal of the GBDT model is to minimise objective function L, expressed as follows: (6) $L=\sum _{i=1}^{n}l(y_i,{\hat{y}}_i)+\sum _{k=1}^K\mathrm{\Omega }(f_k)$(6) where $l(y_i,{\hat{y}}_i)$ is the loss function representing the difference between the prediction result ${\hat{y}}_i$ and target y_i, and Ω(f_k) is the regularisation term that penalises the complexity of the model. In this study, lasso (L1) and ridge (L2) penalty terms were used to avoid overfitting and improve the accuracy of the model (Wang et al. 2021).

The objective function L^(t) at the t-th iteration is minimised by adding a new tree f_t(x_i). L^(t) is expressed using Equations (5(5) ${\hat{y}}_i=\sum _{k=1}^K{f}_k({\mathbf{x}}_i)$(5) ) and (6(6) $L=\sum _{i=1}^{n}l(y_i,{\hat{y}}_i)+\sum _{k=1}^K\mathrm{\Omega }(f_k)$(6) ), as follows: (7) $L^{(t)}=\sum _{i=1}^{n}l(y_i,{\hat{y}}_i^{(t-1)}+f_t({\mathbf{x}}_i))+\mathrm{\Omega }(f_t)$(7) where ${\hat{y}}_i^{(t-1)}$ is the prediction result for x_i at the t-1-th iteration.

In this study, Logarithmic Loss (logloss) was used as the loss function (Zhang and Gong 2020). logloss is calculated as follows: (8) $\mathrm{logloss}=-{\displaystyle \frac{1}{N}\sum _{i=1}^N\sum _{j=1}^C\{I_{i,j}\cdot \log (p_{i,j})\}}$(8) where N is the number of samples, I_i,j is the binary variable, which has a value of 1 when y_i is equal to j and a value of 0 otherwise, and p_i,j is the probability that ${\hat{y}}_i$ is equal to j.

2.1.3. Resampling imbalanced dataset

The dataset used for ML in this study is imbalanced data, where the amount of non-landslide data is significantly greater than that of landslide data. A classifier developed without considering the class imbalance problem can be overwhelmed by the major class and ignores the minor class (Ofek et al. 2017). In this study, the synthetic minority oversampling technique (SMOTE) proposed by Chawla et al. (2002) and random undersampling (RUS) were used to balance the numbers of landslide and non-landslide samples. Synthetic minority class samples (i.e. landslide data) can be generated by SMOTE based on the neighbours randomly selected from the k-nearest neighbours of the minority class, depending on the amount of oversampling required. The value of feature j in the new sample generated by SMOTE x_new(j) is expressed as follows: (9) ${\mathbf{x}}_{\mathrm{new}}(j)={\mathbf{x}}_i(j)+\lambda \{{\mathbf{x}}_z(j)-{\mathbf{x}}_i(j)\}$(9) where x_i(j) is the value of feature j in the sample randomly selected from the minority class samples, x_z(j) is the value of feature j in the sample randomly selected from the minority class k-nearest neighbours of x_i, and λ is a uniform random number in the range [0, 1].

Excessive resampling can cause performance deterioration of the model owing to overfitting, loss of important information, and misinterpretation of the data (Shahabadi et al. 2021). The accuracy of the prediction model depends on the number of training data points after resampling (Hahne et al. 2014). In this study, the optimal number of resampled training data and the ML algorithm for developing a rainfall-induced landslide prediction model were investigated based on performance metrics.

2.1.4. Performance metrics

The performance of the prediction models developed using the Random Forest and LightGBM algorithms was evaluated using accuracy, precision, recall, F-measure and area under the receiver operating characteristic (ROC) curve (AUC) as performance metrics. Accuracy, precision, recall and F-measure are defined as follows: (10) $\mathrm{Accuracy}={\displaystyle \frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{FP}+\mathrm{FN}+\mathrm{TN}}}$(10) (11) $\mathrm{Precision}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}}$(11) (12) $\mathrm{Recall}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}}$(12) (13) $F-\mathrm{measure}={\displaystyle \frac{2\mathrm{Precision}\cdot \mathrm{Recall}}{\mathrm{Precision}+\mathrm{Recall}}}$(13) where TP, TN, FP and FN are the numbers of positive samples correctly predicted as positive, negative samples correctly predicted as negative, negative samples incorrectly predicted as positive, and positive samples incorrectly predicted as negative (i.e. true positive, true negative, false positive, and false negative), respectively. In this study, positive and negative samples represent landslide and non-landslide data, respectively.

The ROC curve is defined as a plot of the false-positive rate (FPR) and true-positive rate (TPR) {FPR(c), TPR(c)}, where $c\in \mathbb{R}$. FPR is defined as follows: (14) $\mathrm{FPR}={\displaystyle \frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}}$(14) TPR is the same as recall defined in Equation (12(12) $\mathrm{Recall}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}}$(12) ). The AUC is calculated as follows: (15) $\mathrm{AUC}={\int }_{\mathrm{\infty }}^{-\mathrm{\infty }}\mathrm{TPR}(c)\mathrm{dFPR}(c)$(15) In this study, the hold-out method was used to calculate the performance metrics. The dataset was randomly divided into 70% for training and 30% for classifier testing. The training data were resampled to equalise the landslide and non-landslide data using SMOTE and RUS, before applying the ML algorithm. This process was repeated three times to reduce the effects of data selection bias on prediction performance. Subsequently, the means of accuracy, precision, recall, F-measure and AUC were calculated for various numbers of resampled training data, and compared to determine the appropriate number of training data after resampling and the optimal ML algorithm.

2.1.5. Determination of optimal prediction model

The performance metrics obtained without resampling are presented in Table 1. There were insignificant differences in the performances of the prediction models developed using Random Forest and LightGBM algorithms. The accuracy was remarkably high, while the recall and AUC were as low as approximately 0.13 and 0.56, respectively. Prediction models developed without considering the class imbalance problem can be biased towards the major class (i.e. the non-landslide class) and have poor classification performance for the minor class (i.e. the landslide class). Figure 5 shows the relationship between the performance metrics and number of resampled training data. Regardless of the ML algorithms, the recall and AUC were dramatically improved by resampling, whereas the precision decreased because of the significant increase in FP. This result demonstrates that the models developed using resampled training data can roughly predict landslide locations, although non-landslide class data can be misclassified as landslide class. It is acceptable that rainfall-induced landslide susceptibility would be overestimated when considering disaster management. The application of resampling techniques is valid for the datasets used in this study, considering that the prediction model is used to determine regions with high rainfall-induced landslide susceptibility.

Figure 5. Relationship between performance metrics and the number of resampled training data: (a) Random Forest and (b) LightGBM.

Table 1. Performance metrics obtained without resampling.

	Accuracy	Precision	Recall	F-measure	AUC
Random Forest	0.99468	0.58903	0.12228	0.20251	0.56090
LightGBM	0.99471	0.58751	0.13828	0.22387	0.56887

As shown in Figure 5, each performance metric converges to a constant value with an increase in the number of resampled training data. In the illustrative example described below, the prediction model developed using the training data resampled to 4,000,000 landslide and non-landslide class data was applied considering the performance metrics and computational cost. The LightGBM prediction model exhibited an equivalent AUC and higher accuracy and F-measure than the Random Forest prediction model; thus, it was determined to be the optimal ML algorithm. Table 2 lists the hyperparameter values for Random Forest and LightGBM algorithms used in this study. The hyperparameters were determined by parameter tuning based on the performance of prediction models developed using resampled training data with 4,000,000 samples for landslide and non-landslide classes.

Table 2. Hyperparameter values for Random Forest and LightGBM algorithms.

Hyperparameter	Value
Random Forest
Number of trees	3000
Maximum depth of a tree	30
Maximum number of features	3 (square root of total number of features)
LightGBM
Number of trees	5000
Maximum number of leaves in a tree	2047
Learning rate	0.01
Coefficient for L1 regularisation term	1.00
Coefficient for L2 regularisation term	1.00

The spatial distribution of the results obtained using the LightGBM prediction model in the data collection areas is shown in Figure 6. The FN that represents the landslide grids incorrectly predicted as non-landslide is distributed around the TP, which represents the correctly predicted landslide grids. The regions with densely populated FP, representing the non-landslide grids that are incorrectly predicted as landslides, are mainly distributed in mountainous areas with plutonic rock, as shown in Figures 3 and 4, where the landslide susceptibility is overestimated. These results demonstrate that the prediction model can avoid overestimation of landslide safety and can be used to determine regions susceptible to landslides triggered by rainfall.

Figure 6. Spatial distribution of results obtained using LightGBM prediction model in data collection areas: (a) Heavy Rain Event of July 2018 and (b) Typhoon Hagibis.

2.1.6. Feature importance

Feature importance is an indicator of the extent to which each feature contributes to model prediction. The feature importance score is calculated as follows: (16) $F{I}_n(j)={\displaystyle \frac{\sum _{v=1}^{N_j}\mathrm{Gain}(v,j)n_v}{\sum _{j=1}^m\sum _{v=1}^{N_j}\mathrm{Gain}(v,j)n_v}}$(16) where FI_n(j) is the normalised feature importance score of feature j, N_j is the set of nodes in which feature j is used as the splitting criterion, Gain(v,j) is the information gain, and n_v is the number of training data points in node v that belong to N_j.

Figure 7 shows the feature importance scores of the LightGBM prediction model. The feature importance can be used to select an appropriate set of features to improve the performance of the prediction model. In this study, feature selection was not performed, and all nine features were used in ML since features with a significantly low feature importance score were not identified.

Figure 7. Feature importance scores of LightGBM prediction model.

2.2. Rainfall-induced landslide susceptibility assessment

Landslide susceptibility can be assessed by combining the ML-based prediction model with the rainfall hazard curve. In this study, a generalised extreme value (GEV) distribution was applied to estimate the hazard curve associated with the annual maximum rainfall index. The GEV approach has been widely used to model extreme rainfall event values (Paola et al. 2018). The GEV distribution was applied to fit the annual maximum rainfall index calculated for each year of the observation period using Equations (1(1) $R_{\mathrm{fw}}=\sqrt{{(R_1-R_w)}^2+a^2{(r_1-r_w)}^2}$(1) ) and (2(2) $R^{\mathrm{\prime }}=R_{\mathrm{fw}0}-R_{\mathrm{fw}}$(2) ) based on historical rainfall data. The GEV cumulative distribution function is expressed as follows: (17) $F(x_r)=\exp \left\{-{\left[1+\gamma \left({\displaystyle \frac{x_r-\mu }{\theta }}\right)\right]}^{-1/\gamma }\right\}$(17) where x_r is the rainfall index, F(x_r) is the GEV cumulative distribution function, and γ, μ and θ are the shape, location and scale parameters, respectively. The return period associated with the rainfall index x_r, T(x_r), is calculated as follows: (18) $T(x_r)={\displaystyle \frac{1}{1-F(x_r)}}$(18) The denominator in Equation (18(18) $T(x_r)={\displaystyle \frac{1}{1-F(x_r)}}$(18) ), 1 − F(x_r), represents the annual exceedance probability associated with x_r (i.e. the rainfall hazard curve). The values of the rainfall indices with various return periods were generated using Equations (17(17) $F(x_r)=\exp \left\{-{\left[1+\gamma \left({\displaystyle \frac{x_r-\mu }{\theta }}\right)\right]}^{-1/\gamma }\right\}$(17) ) and (18(18) $T(x_r)={\displaystyle \frac{1}{1-F(x_r)}}$(18) ), and input into the LightGBM prediction model for rainfall-induced landslide susceptibility mapping.

2.3. Risk assessment of economic loss from structures damaged by landslides

The structures affected by rainfall-induced landslides can be determined by overlaying the landslide susceptibility and structure distribution maps in GIS. Landslides can cause damage to structures and associated economic losses (Hashimoto et al. 2020). The economic loss caused by a structure damaged by a landslide can be calculated based on Ishibashi et al. (2021b) as follows: (19) $Q_s=\sum _{\mathrm{ds}=1}^{N_{\mathrm{ds}}}{P}_{\mathrm{ds}}\cdot w_{\mathrm{ds}}\cdot \alpha \cdot A_s$(19) where subscript s and ds are the structure number in the analysed area and the damage state of the structure, respectively, Q_s is the economic loss from the structure affected by rainfall-induced landslides, P_ds is the failure probability that the damage state of the structure becomes ds due to landslides, w_ds is the damage ratio for the damage state, α is the monetary value per unit area of a structure (=186.1 thousand JPY/m² based on MLIT (2023)), and A_s is the building area of the structure. The damage states of individual structures are classified into four levels based on the disaster damage certification criteria provided by the Cabinet Office, Government of Japan (2020): slight, moderate, extensive and complete. Katano et al. (2020) statistically organised the characteristics of structure damages due to past landslides triggered by rainfall in Japan. Based on Katano et al. (2020), P_ds for the slight, moderate, extensive and complete damage states were set to 0.347, 0.224, 0.106 and 0.324, respectively. The damage ratio was used to estimate the economic loss from individual structures considering the extent of damage. w_ds for the slight, moderate, extensive and complete damage states were determined to be 0.1, 0.3, 0.75 and 1.0, respectively, in accordance with the definition of the disaster damage certification criteria (Cabinet Office, Government of Japan 2020) and the report on structural damage after natural disasters (Rogers et al. 2015).

The total economic loss caused by all damaged structures in the analysed area Q is calculated as follows: (20) $Q=\sum _{s=1}^{N_s}{Q}_s$(20) where N_s is the total number of structures affected by rainfall-induced landslides in the analysed area.

The risk curve associated with the economic loss from structures damaged by rainfall-induced landslides is estimated as follows: (21) $P(Q\ge q)=1-{\int }_0^q\left|-{\displaystyle \frac{d{F}_Q(q)}{\mathrm{dq}}}\right|\mathrm{dq}$(21) where P(Q ≥ q) is the probability that Q exceeds q, and F_Q(q) is the cumulative probability distribution of the economic loss from the structures damaged by rainfall-induced landslides.

3. Illustrative example

3.1. Analysed area

The risk of economic loss from structures affected by rainfall-induced landslides in the area located in Fukushima Prefecture, as shown in Figure 2, was analysed using the proposed methodology. The analysed area included two major cities (Koriyama and Iwaki Cities). Risk assessments in cities and surrounding mountainous areas can contribute to determining effective landslide mitigation countermeasures and enhancing disaster resilience.

Figure 8 shows the distribution maps of the dataset for risk assessment in the analysed area. The analysed area was divided into grid cells with a resolution of 10 m using GIS. As with the areas used to collect the ML dataset, the decomposed granite soil belonging to the plutonic rock is widely distributed in the analysed area, which can be vulnerable to landslides triggered by rainfall (Award 1988). The applicability of the LightGBM prediction model to the analysed area is ensured since the landslide conditioning factors, such as topography, soil and lithology, are similar for the data collection areas shown in Figures 3 and 4.

Figure 8. Distribution maps of dataset for risk assessment in the analysed area.

3.2. Rainfall hazard analysis

Rainfall hazard analysis was performed using historical rainfall data obtained from 1976 to 2021 at the meteorological station presented in Figure 8(a). Figure 9 shows the hazard curve associated with the rainfall index estimated using Equations (17(17) $F(x_r)=\exp \left\{-{\left[1+\gamma \left({\displaystyle \frac{x_r-\mu }{\theta }}\right)\right]}^{-1/\gamma }\right\}$(17) ) and (18(18) $T(x_r)={\displaystyle \frac{1}{1-F(x_r)}}$(18) ). Takasao, Takara, and Shimizu (1986) proposed the Standard Least-Squares Criterion (SLSC), which is an indicator of fit across probability density functions and used as a standard goodness-of-fit criterion in practice in Japan (Takara 2021). Lower SLSC values indicate better fits to the data, and SLSC=0.04 is generally considered a threshold for well-fitted models (Tanaka and Takara 1999). The SLSC value of the GEV distribution shown in Figure 9 was calculated to be 0.029, indicating a good fit with the rainfall index. In this illustrative example, the rainfall hazard curve presented in Figure 9 was applied to the entire analysed area considering the distribution of meteorological stations. XRAIN, which is a real-time and spatially dense rainfall observation system in Japan, has been in operation for approximately ten years. Continuous data collection from a rainfall observation system with a high spatiotemporal resolution can contribute to rainfall hazard analysis at each grid location.

Figure 9. Hazard curve associated with rainfall index.

3.3. Landslide susceptibility mapping

Landslide susceptibility was assessed by inputting the data on the landslide conditioning factors and the rainfall index values estimated based on the hazard curve into the LightGBM prediction model. Figure 10 shows the landslide susceptibility map based on the return period associated with the rainfall index and the locations of Koriyama and Iwaki Cities. Rainfall indices with return periods between 2 and 500 years were used to assess landslide susceptibility in the analysed area.

Figure 10. Landslide susceptibility map based on return period associated with rainfall index.

Figure 11 shows the relationship between the annual exceedance probability and the number of landslides in the entire analysed area as well as in Koriyama and Iwaki Cities. The number of landslides can dramatically increase when the annual exceedance probability of the rainfall index is less than approximately 0.01 (i.e. when the return period is greater than 100 years). Iwaki City has a larger number of landslides than Koriyama City, regardless of the annual exceedance probability. Since low landslide susceptibility areas, such as lowland and built-up areas, are widely distributed in the Koriyama and Iwaki Cities, the number of landslides in these cities is not significantly large compared to that in the entire analysed area.

Figure 11. Relationship between annual exceedance probability and number of landslides.

3.4. Results of risk assessment

The layers, including the grid data with landslide susceptibility and structure locations, were overlaid in the GIS to identify the structures exposed to rainfall-induced landslides. The structure distribution map of the analysed area is shown in Figure 12. Structures can be damaged not only by landslides occurring in the grid where the structure is located but also by nearby landslides. Structures potentially affected by landslide mass movements should be considered in economic risk assessment using Equations (19(19) $Q_s=\sum _{\mathrm{ds}=1}^{N_{\mathrm{ds}}}{P}_{\mathrm{ds}}\cdot w_{\mathrm{ds}}\cdot \alpha \cdot A_s$(19) ) and (20(20) $Q=\sum _{s=1}^{N_s}{Q}_s$(20) ). In this illustrative example, a 50 m buffer zone around each structure was used when determining the structures exposed to landslides by combining landslide susceptibility and structure distribution maps in GIS (Kamp et al. 2008).

Figure 12. Structure distribution map of the analysed area.

Figure 13 shows the risk curve associated with economic loss from structures damaged by rainfall-induced landslides. The economic risk can increase particularly when the annual exceedance probability is less than 0.01 since many structures would be affected by landslides triggered by rainfall. It is important to implement disaster mitigation countermeasures in the analysed area, considering the severe and extensive damage to structures with a low annual exceedance probability. In addition, Koriyama City has a higher risk than Iwaki City when the annual exceedance probability of economic loss ranges between 0.01 and 0.06 (i.e. the return period is between 15 and 100 years), although the number of landslides in Iwaki City is greater than that in Koriyama City in this range of annual exceedance probability, as shown in Figure 11. This result demonstrates that not only the landslide susceptibility but also the risk of economic loss from damaged structures needs to be considered to improve landslide disaster management. Countermeasures against rainfall-induced landslides in Koriyama City should be prioritised to efficiently mitigate the risk with a high annual exceedance probability, and countermeasures should be taken in Iwaki City to prepare for high-magnitude and low-frequency rainfall events. The risk curve is useful for developing an appropriate disaster mitigation strategy to reduce negative impacts on communities in areas potentially affected by rainfall-induced landslides.

Figure 13. Risk curve associated with economic loss from structures damaged by rainfall-induced landslides.

4. Conclusions

A framework for estimating the risk of economic loss from structures damaged by rainfall-induced landslides is proposed. The risk was assessed based on a landslide prediction model developed using ML, the rainfall hazard curve estimated by the GEV distribution using historical rainfall data, the spatial distribution of structures, and the monetary value of individual structures damaged by landslides. As an illustrative example, the risk curve representing the relationship between the economic loss from the structures affected by rainfall-induced landslides and its annual exceedance probability in an area located in Fukushima Prefecture of Japan was estimated. The following conclusions are drawn:

1. The grid data with a resolution of 10 m obtained in the areas affected by the Heavy Rain Event of July 2018 and Typhoon Hagibis were used for ML. The rainfall index was used as a feature associated with rainfall intensity, considering the effects of both long- and short-term rainfall on landslide occurrences.

2. Resampling techniques based on SMOTE and RUS were applied to improve the performance of ML-based landslide prediction models considering the class imbalance problem of the dataset. Random Forest and LightGBM algorithms were used to develop models for solving the binary classification problem of landslide/non-landslide classes. LightGBM was determined to be the optimal ML algorithm based on the performance metrics. Regions with high rainfall-induced landslide susceptibility can be determined using the LightGBM prediction model.

3. GEV distribution was used to estimate the hazard curve associated with the annual maximum rainfall index. Landslide susceptibility was assessed by inputting the values of the rainfall indices with various return periods generated based on the hazard curve and the data on landslide conditioning factors obtained in the analysed area into the ML-based prediction model.

4. The number of structures exposed to landslides was calculated by overlaying the landslide susceptibility and structure distribution maps in GIS. The risk of economic loss from structures damaged by rainfall-induced landslides in the area located in Fukushima Prefecture of Japan was estimated based on the number of damaged structures and the monetary value of individual structures. The economic loss from structures affected by landslides triggered by rainfall would be huge in the analysed area, especially when the annual exceedance probability is less than 0.01 (i.e. a return period of greater than 100 years). Even though the expected number of landslides is low, the risk can be higher because of the spatial distributions of landslides and structures, depending on the annual exceedance probability. The proposed methodology can contribute to decision-making on countermeasures against rainfall-induced landslides and the development of appropriate disaster mitigation strategies for enhancing resilience.

Data availability statement

Data that support the findings of this study are available from the corresponding author upon reasonable request.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by JSPS Grant-in-Aid for Scientific Research (KAKENHI): [Grant Number JP21K14235].