Spatiotemporal Landslide Susceptibility Mapping Incorporating the Effects of Heavy Rainfall: A Case Study of the Heavy Rainfall in August 2021 in Kitakyushu, Fukuoka, Japan

Abstract

Landslide represents an increasing menace causing huge casualties and economic losses, and rainfall is a predominant factor inducing landslides. Landslide susceptibility assessment (LSA) is a commonly used and effective method to prevent landslide risk, however, the LSA does not analyze the impact of the rainfall on landslides which is significant and non-negligible. Therefore, the spatiotemporal LSA considering the inducing effect of rainfall is proposed to improve accuracy and applicability. In this study, the influencing factors are selected using the chi-square test, out-of-bag error and multicollinearity test. The spatial LSA are thus obtained using the random forest (RF) model, deep belief networks model and support vector machine, and compared using receiver operating characteristic curve and seed cell area index to determine the optimal assessment result. According to the heavy rainfall characteristics in the study area, the rainfall period is divided into four stages, and the effective rainfall model is employed to generate the rainfall impact (RI) maps of the four stages. The spatiotemporal LSAs are obtained by coupling the optimal spatial LSA and various RI maps and verified using the landslide warning map. The results demonstrate that the optimal spatiotemporal LSA is obtained using the spatial LSA of the RF model and temporal LSA of the rainfall data in the peak stage. It can predict the area where rainfall-induced landslides are likely to occur and prevent landslide risk.

1. Introduction

Landslides are one of the most common geo-hazards on Earth and result in considerable damage to life and property [1,2]. It is critical to accurately evaluate landslide-prone areas to effectively determine the locations and times of landslides [3]. Landslide susceptibility assessment (LSA) is thus an effective risk management measure to prevent landslides and reduce landslide risk [4].

Rainfall-induced landslides are a typical cascading geo-hazard and are attracting the increasing attention of many researchers [5,6]. Rainfall-induced landslides commonly occur after heavy rainfall and may develop into potentially catastrophic movements [7,8]. Heavy rainfall or prolonged rainfall may result in the increase of pore-water pressure and the significant decrease of friction force, which causes slope instability [9,10,11]. Hence, in addition to analyzing the spatial LSA using the environmental condition factors, it is also necessary to analyze the temporal LSA using the rainfall data to illustrate the impact of heavy rainfall on landslides, which ensures the availability and accuracy of LSA [12,13,14]. Spatiotemporal landslide susceptibility mapping (LSM) is thus studied as an excellent method of susceptibility assessment and provides information on potential landslides under the impact of heavy rainfall [15,16].

The analysis of spatial LSA generally assumes that future landslides are more likely to occur in similar environmental conditions to previous landslides [17]. The LSA thus can predict the potential risk areas using the historical landslide data and environmental condition factors [18]. The selection of influencing factors is an essential part of LSA, and most current research selects the influencing factors using correlation test and multicollinearity test to prove the correlation between influencing factors and landslides and avoid a redundancy of influencing factors [19,20,21].

Meanwhile, the LSA results are affected by evaluation models, and it is very important to select the appropriate models for LSA [22,23,24]. There are various models to analyze the LSA, and these are generally classification models [25]. Early studies map the landslide susceptibility using traditional statistical analysis models, such as linear discriminant analysis, weights of evidence, frequency ratio models and regression models [26,27,28,29]. Most current studies propose various machine learning algorithms, such as support vector machines (SVM), naive Bayesian model, random forest (RF) models, and artificial neural networks, to evaluate the landslide susceptibility [28,30,31,32,33]. Because of the complex environmental conditional factors of landslides in study areas, the accuracy and applicability of LSA obtained using various models are very different. It is thus essential to select the evaluation models that are more applicable to the study area.

On the other hand, rainfall is one of the most critical inducing factors, especially heavy and prolonged rainfall, and rainfall analysis is thus one of the core steps in the temporal LSA [34,35]. The impact of rainfall intensity and duration on landslides has long been an issue of great interest for researchers in LSA [36,37]. The effective rainfall model (ERM) is a common and effective method and employs the rainfall data of some days before the study time to indicate the impact of rainfall, which can illustrate the temporal landslide susceptibility of the study area [38,39].

The present study proposes a methodology to evaluate spatiotemporal landslide susceptibility incorporating the effects of heavy rainfall, and the methodology is applied to the Kitakyushu, Fukuoka, Japan. Multi-source data such as the condition factors, historical landslides, heavy rainfall data and landslide warning map are considered. In the spatial LSA, the chi-square (${\chi }^2$) test, out-of-bag (OOB) error and multicollinearity test are employed to select the influencing factors, and the spatial LSMs are thus obtained using the RF model, deep belief networks (DBN) model and SVM model. The accuracies of the spatial LSMs are proved using receiver operating characteristic (ROC) curve and seed cell area index (SCAI) to obtain the optimal spatial LSM used for spatiotemporal LSM. Based on the heavy rainfall data in August of the study area, the heavy rainfall period is divided into four stages, and the rainfall impact (RI) maps of the four stages are thus generated using ERM. The spatiotemporal LSMs are obtained by coupling the optimal spatial LSM and temporal LSMs and verified using the landslide warning map to determine the optimal spatiotemporal LSM, which provides a reference for landslide prevention.

2. Materials and Methods

2.1. Methodological Flow

The steps and processes of the LSM procedures involved in the present study (Figure 1) are as follows:

The methodology framework of the proposed LSM can be divided into seven steps: (a) data collection: the data of various condition factors of landslides and the historical landslides are collected, and heavy rainfall data in August and landslide warning map are prepared. (b) Influencing factor selection: the influencing factors are selected as the LSA indices using the ${\chi }^2$ test, OOB error and multicollinearity test. (c) LSA: the LSMs are obtained using the RF model, DBN model and SVM model. (d) Result verification: the results of LSA are verified using ROC curve, area under the curve (AUC) and SCAI to determine the optimal spatial LSA. (e) Rainfall map generation: the heavy rainfall period in the study area is divided into four stages, namely the development stage, peak stage, stagnation stage and decline stage, and the rainfall maps of various stages are generated using the interpolation method, and the RI maps of various stages are obtained and normalized by the effective rainfall model. (f) Spatiotemporal LSA: various RI maps are coupled with the optimal spatial assessment result to obtain the spatiotemporal susceptibility maps under the influence of heavy rainfall in August. (g) Spatiotemporal assessment verification: the various spatiotemporal LSMs are verified using the landslide warning map to determine the optimal spatiotemporal LSM.

2.2. Factor Selection Methods

The ${\chi }^2$ test is one of the most common methods to determine the correlation between influencing factors and landslides, and the variance inflation factor (VIF) is a commonly used index to test multicollinearity which refers to the linear correlation between independent variables [19,20,21]. The ${\chi }^2$ test determines the accuracy of the hypothesis by calculating the deviation of the actual value from the theoretical value, and the ${\chi }^2$ can be obtained using Equation (1). The hypothesis is verified by comparing the ${\chi }^2$ with the critical value of a significance test.

$${\chi }^2={{\displaystyle \sum }}_{i=1}^n\frac{{\left(A_i-T\right)}^2}{T},$$

(1)

where $A_i$ is the $i^{th}$ actual value and $T$ is the theoretical value.

Multicollinearity of influencing factors is verified using VIF. The smaller the VIF value, the less significant the multicollinearity, and the greater the VIF value, the more significant the multicollinearity. The VIF value is thus closer to 10 indicating stronger multicollinearity between the influencing factors [26,40,41]. The VIF can be obtained using Equation (2).

$$VIF=\frac{1}{1-R_i^2},$$

(2)

where $R_i$ is the correlation coefficient of the $i^{th}$ independent variable to other independent variables.

The OOB error is a widely employed method to select the influencing factors and determine their significance. The OOB error is calculated using the ratio of misclassification number to the total number. In the present study, based on the recursive feature elimination (RFE), the OOB error is thus used to eliminate the factor with the least importance using variable importance measures (VIM) and to determine the optimal number of landslide influencing factors [42].

2.3. Landslide Susceptibility Assessment (LSA) Methods

The RF model and SVM model are currently commonly used models [28,31], while the DBN model is a relatively novel model used for the LSA [43]. The applicability of the three models for the LSA is proved, and the present study thus selects the three models, namely the RF model, DBN model and SVM model to analyze the LSA.

2.3.1. Random Forest (RF) Model

The RF model was proposed as an advanced integrated learning algorithm using the ensemble of unpruned classification trees in 1995. It is a classifier training samples by multiple trees [44,45]. The trees are constructed according to the following four steps: (a) N denotes the number of training samples, while M denotes the number of characteristics. (b) The number of characteristics is inputted to determine the decision result of the previous node on the decision tree. (c) The training set is generated using the bootstrap sampling method, and the unselected samples are used to predict and evaluate the errors. (d) The characteristics are selected randomly to determine the decision of each node on the decision tree, and the optimal splitting method is calculated [46].

2.3.2. Deep Belief Networks (DBN) Model

The DBN model is an efficient unsupervised learning algorithm and is also a probabilistic generative model [47]. Compared with the traditional neural network, it generates a joint distribution between the observation data and the label. The DBN model is composed of restricted Boltzmann machines (RBM) which consists of the visual layer $\upsilon$ and hidden layer $h$. There is no connection between the neurons in a layer, while there are connections between layers, and every layer can be considered as a logistic regression model. The hierarchical structure is very close to the structure of the human brain. Several RBMs structure a DBN, and the hidden layer of the previous RBM is the visual layer of the current RBM. The structure of the DBN model is presented in Figure 2.

In RBM, each neuron has a bias coefficient $b$ for visual neuron and $c$ for hidden neuron, and the weight $W$ between two connected neurons is the connection weight. The energy of the RBM can be obtained using Equation (3).

$$E\left(\upsilon ,h\right)=-{{\displaystyle \sum }}_{i=1}^{N_{\upsilon }}{b}_i{\upsilon }_i-{{\displaystyle \sum }}_{j=1}^{N_h}{c}_j{h}_j-{{\displaystyle \sum }}_{i,j=1}^{N_{\upsilon },N_h}{W}_{ij}{\upsilon }_i{h}_j,$$

(3)

where $N_{\upsilon }$ is the number of visual neurons; $N_h$ is the number of hidden neurons.

2.3.3. Support Vector Machines (SVM) Model

The SVM model is a classifier developed from the generalized portrait algorithm in pattern recognition. It can be used for classification, regression and outlier detection based on the principle of minimizing statistical and structural risk. The sample data are normalized and mapped using kernel functions to build the SVM model, and the hyperparameters of the SVM model are optimized to train the model. The SVM model can be performed using various kernel functions, such as linear kernel (Equation (5)), polynomial kernel (Equation (6)), radial basis function kernel (Equation (7)) and sigmoid kernel (Equation (8)). The main control hyperparameters are the penalty parameter and kernel function parameter. The penalty parameter is introduced, and the objective function can be obtained using Equation (4):

$$f\left(x\right)=min\frac{1}{2}\Vert w{\Vert }^2+C{{\displaystyle \sum }}_{i=1}^n{\xi }_i,$$

(4)

where $\frac{1}{2}\Vert w{\Vert }^2$ is the geometric interval; $C$ is the penalty parameter; and ${\xi }_i$ is the slack variable. And the various kernel functions can be obtained using Equations (5)–(8):

$$K\left(x_i,x\right)=x_i·x,$$

(5)

$$K\left(x_i,x\right)={[(x_i·x)+1]}^m,$$

(6)

$$K\left(x_i,x\right)=\exp \left(-\frac{\Vert x_i-x,{\Vert }^2}{2{\delta }^2}\right),$$

(7)

$$K\left(x_i,x\right)=\tan \mathrm{h}(\eta <x_i,x>+\theta ).$$

(8)

2.4. Effective Rainfall Model (ERM)

The ERM reflects the rainfall accumulation that stays in the rock and contributes to the development of landslides. The model was proposed by Kohler and Linsley in 1951 as an antecedent precipitation index to predict the surface runoff based on hydrological data [48]. The effective rainfall can be obtained using Equation (9).

$$R_e={{\displaystyle \sum }}_{i=0}^n{K}^i{R}_i,$$

(9)

where $K$ is the attenuation coefficient and considered as 0.8; and $R_i$ is the rainfall of the previous $i^{th}$ day.

In the present study, the rainfall data are collected from meteorological observation stations in the study area. The maps of rainfall during the study time are thus generated using the kriging interpolation method, and the effective rainfall map can be obtained using the ERM.

2.5. Validation Methods

Verifying the accuracy of susceptibility assessment results obtained using various models is a key task of susceptibility research. In the present study, the SCAI and ROC curve methods were selected to verify the accuracy of susceptibility assessment. The SCAI can be obtained using Equation (10).

$$SCAI=\frac{P_i}{P{’}_i},$$

(10)

where $P_i$ is the areal percentage of the susceptibility level $i$; and ${P^{\prime }}_i$ is the areal percentage of the landslides at level $i$.

The ROC curve is a commonly used method to verify the accuracy of assessment results [49]. The ROC curve was obtained using two variables, namely the specificity and sensitivity. The variables can be obtained using Equations (11) and (12):

$$Specificity=\frac{TN}{FP+TN},$$

(11)

$$Sensitivity=\frac{TP}{TP+FN},$$

(12)

where $FP$ is the false positive; $TN$ is the true negative; $TP$ is the true positive; and $FN$ is the false negative.

The threshold for AUC values is 0.5 to 1, and the closer it is to 1, the more accurate the model is [50].

3. Study Area and Dataset

3.1. Study Area

In the present study, Kitakyushu, belonging to Fukuoka in Japan, was selected as the study area (Figure 3). Geographically, the area is located within 33°58′00″ N to 33°43′00″ N and 130°40′00″ E to 131°01′00″ E, with an area of 488.78 km². Kitakyushu is on the northern Kyushu Island, and the southern side of the area is the hibikinada of the Japan Sea, while its eastern side is Suwou of SetoNaikai. Geologically, the whole area is characterized by rugged topography and tilting terrain from north to south. The highest and lowest altitudes of the study area are about 900 m and –68 m, respectively. Meanwhile, the area can be mainly divided into four regions: southern mountain region, central plain region, northeastern mountain region and northwestern hilly region. In accordance with the Ministry of Land, Infrastructure, Transport and Tourism of Japan, the study area has active geological tectonic movements with complex geological conditions. The lithology of the study area dominates sedimentary rock and igneous rocks, with a landfill area of more than 5% of the study area.

The annual average precipitation of the study area is about 1265 mm, which is the main inducing factor for landslides. The area is warm and humid throughout the year with an average annual temperature of 16.2 °C. Under the influence of the geographical environment, the northern region has a typical Japan Sea climate, while the eastern region has a Seto Inland Sea climate. The precipitation significantly varies, concentrated in the rainy season and typhoon season. Because the study area is located in the Pacific Rim Volcanic Seismic Zone at the junction of Eurasian and Pacific plates, it has frequent crustal movements. The heavy rainfall and crustal movement are thus factors inducing landslides.

Kitakyushu and Fukuoka are considered as the center of the metropolitan area and the hub of all administrative activities of the Kyushu area. According to the census of Japan, 2017, Kitakyushu has a population of 935,000. According to the Bureau of Land Policy in the Ministry of Land, Infrastructure, Transport and Tourism of Japan, the historical landslide areas can be obtained in Figure 3. The frequent landslides thus restrict the development of the industry and the socio-economic growth of the study area.

3.2. Condition Factors

Taking the appropriate relevant factors into consideration is essential for evaluating landslide susceptibility. Based on the previous survey study of the study area, the following 14 condition factors were collected and prepared: soil thickness (ST), elevation, slope, topographic wetness index (TWI), topography, cumulative runoff (CR), lithology, curvature, plane curvature, profile curvature, distance to the road (DRO), distance to the railway (DRA), distance to a river (DRI) and distance to a structure line (DSL).

The data of the condition factor were obtained by the field environmental survey and related literature research. The digital elevation model (DEM) data were provided by the Geospatial Information Authority of Japan, and the maps of elevation, slope, TWI, curvature, plane curvature and profile curvature in the study area were thus prepared. The geology factors, namely the soil thickness, topography, cumulative runoff, structure line and lithology, were provided by the Ministry of Land, Infrastructure, Transport and Tourism of Japan. Information about the roads, railways and rivers was provided by the Ministry of Land, Infrastructure, Transport and Tourism Kyushu Regional Development Bureau. The factors involve both quantitative and qualitative elements. The quantitative factors are divided into five levels using the natural breakpoint method, and Figure 4 shows the maps of the 14 factors.

3.3. Heavy Rainfall Data and Landslide Warning Map

Heavy rainfall is one of the most significant landslide-inducing factors, and because of the climate characteristic of the rainy season in the study area, there always is heavy rainfall in August. The daily rainfalls in August at various monitoring stations are collected from the Japan Meteorological Agency of Ministry of Land, Infrastructure, Transport and Tourism (Figure 5).

As can be seen from Figure 5, the heavy rainfall started on 7 August and ended on 27 August. According to the rainfall variability and trends in Figure 5, the heavy rainfall can be divided into four stages, namely the development stage (from 7 August to 10 August), peak stage (from 11 August to 15 August), stagnation stage (from 16 August to 20 August) and decline stage (from 21 August to 27 August). Meanwhile, the geo-hazard warning map after the heavy rainfall also can be collected from the Japan Meteorological Agency of Ministry of Land, and the map is shown in Figure 6.

4. Results

4.1. Selection of Spatial Influencing Factors

The spatial influencing factors are selected using the ${\chi }^2$ test, OOB (out-of-bag) error and multicollinearity test. The ${\chi }^2$ values of the 14 factors can be obtained using Statistical Product and Service Solutions (SPSS) in

Table 1. ${\chi }^2$ values of the 14 factors.

Factors	${\mathit{\chi }}^2\mathbf{Values}$	Factors	${\mathit{\chi }}^2\mathbf{Values}$
Soil thickness	36.77	Curvature	5.85
Elevation	41.27	Plane curvature	8.15
Slope	31.71	Profile curvature	3.10
TWI	9.56	DRO	3.10
Topography	27.77	DRA	3.22
Cumulative runoff	0.73	DRI	10.79
Lithology	14.73	DSL	4.60

.

The ${\chi }^2$ values of various factors are compared with the critical value (k = 3.84), and the factors whose ${\chi }^2$ values are smaller than the critical value are eliminated. Therefore, 10 factors, namely the soil thickness, elevation, slope, TWI, topography, lithology, curvature, plane curvature, DRI and DSL, are selected as evaluation indices. According to the RFE, the OOB errors of various factor selection and the VIMs of selected factors are calculated, and the factor with the least VIM is eliminated. The process is repeated, and the OOB errors and VIMs are thus shown in Figure 7.

As can be seen from Figure 7a, the number of factors in the optimal factor selection is 10, and the VIMs of the 10 factors are shown in Figure 7b. The VIFs of the 10 factors are then obtained to verify the multicollinearity in

Table 2. Variance inflation factors (VIFs) of the 10 factors.

Factors	VIFs	Factors	VIFs
Slope	1.327	Lithology	1.105
Soil thickness	1.186	Plane curvature	1.109
Elevation	1.446	Curvature	1.103
Topography	1.384	TWI	1.173
DRI	1.078	DSL	1.033

.

As shown in Table 2, the VIFs are far less than 10, and this indicates there is no strong correlation between the factors, which proves that there is no redundancy in the factor selection.

4.2. Spatial Landslide Susceptibility Mapping (LSMs) and Validation

The spatial landslide susceptibility of the study area is evaluated using the RF model, DBN model and SVM model. The landslide region is set to 1, while the non-landslide region is set to 0, and the map of historical landslide region is coupled with the maps of selected influencing factors. We selected 70% of grid data as the training data to predict the landslide susceptibility of each grid in the study area. The spatial LSMs were thus divided into five categories—very low susceptibility, low susceptibility, moderate susceptibility, high susceptibility and very high susceptibility—using a natural breakpoint method and shown in Figure 8, and the area percentages of various susceptibility levels are thus presented in

Table 3. Area percentages of various susceptibility levels.

Models	Very Low	Low	Moderate	High	Very High
RF model	15.25%	15.30%	40.05%	27.45%	1.95%
DBN model	12.27%	9.48%	53.62%	19.13%	5.50%
SVM model	20.14%	1.16%	71.37%	6.24%	1.09%

.

As shown in the figure and table above, the distribution of the LSMs is similar, and the very low susceptibility area is mainly concentrated in the northern part of the study area. The area with the moderate susceptibility level is the largest area of the study area, while the area with the very high susceptibility level is the smallest. Then, the SCAI and ROC curve are used to verify the accuracy of the three LSMs using the RF model, DBN model and SVM model. The SCAI is an effective method to verify the accuracy of the three assessment results, and the SCAIs of various models are shown in

Table 4. Seed cell area indexes (SCAIs) of RF model, DBN model and SVM model.

Model	Level	Area Percentage	Landslide Area (km2)	Landslide Percentage	SCAI
RF model	Very low	15.25%	0.208	0.24%	63.40
	Low	15.30%	12.888	14.90%	1.03
	Moderate	40.05%	40.603	46.95%	0.85
	High	27.45%	26.342	30.46%	0.90
	Very high	1.95%	6.443	7.45%	0.26
DBN model	Very low	12.27%	0.367	0.42%	28.91
	Low	9.48%	0.482	0.56%	17.01
	Moderate	53.62%	57.368	66.34%	0.81
	High	19.13%	18.391	21.27%	0.90
	Very high	5.50%	9.876	11.42%	0.48
SVM model	Very low	20.14%	1.668	1.93%	10.44
	Low	1.16%	2.396	2.77%	0.42
	Moderate	71.37%	73.135	84.57%	0.84
	High	6.24%	4.955	5.73%	1.09
	Very high	1.09%	4.329	5.01%	0.22

.

Table 4 shows that the very high susceptibility has the smallest SCAI value, while the very low susceptibility has the greatest SCAI value. It indicates that the three models are applicable to the LSA in the study area. Meanwhile, the ROC curve and AUC are methods commonly used to evaluate quantitatively the accuracy of the assessment result, and they are thus drawn in Figure 9.

The ROC curves and AUCs illustrate that the RF model presents the best LSA performance. The RF model can be considered as the optimal model in the three models for the LSM in the study area.

4.3. Spatiotemporal LSMs and Validation

There is always heavy rainfall in August in the study area, and the heavy rainfall is divided into four stages, namely the development stage, peak stage, stagnation stage and decline stage. The effective rainfall in the four stages is obtained using ERM to analyze the spatial impact of rainfall in the various stages on the landslide susceptibility. According to the rainfall data collected from various monitoring stations, the daily rainfall maps during the heavy rainfall period are thus obtained using kriging interpolation. The effective rainfall maps of the four stages are then generated, and the RI is normalized and calculated using Equation (13).

$$RI=\frac{ER-E{R}_{min}}{E{R}_{max}-E{R}_{min}},$$

(13)

where $ER$ is the effective rainfall of the grid; $E{R}_{min}$ is the minimum effective rainfall in the study area; and $E{R}_{max}$ is the maximum effective rainfall in the study area. The RI maps are thus obtained and presented in Figure 10.

The RI maps of the four stages are coupled to the spatial assessment result using the optimal model. The temporal landslide probability, namely the RI, in the study area is thus accumulated with the spatial landslide probability obtained using the RF model to obtain the spatiotemporal LSA. The accumulation results are then divided into five categories: very low susceptibility, low susceptibility, moderate susceptibility, high susceptibility and very high susceptibility using the natural breakpoint method, and the spatiotemporal LSMs of the four stages are generated and shown in Figure 11.

Figure 11 shows that the spatiotemporal LSMs are tightly associated with the temporal LSAs, namely the RI maps. The area with the very high susceptibility is highly concentrated in the central-south region in the study area, while the area with the very low susceptibility is mainly located in the northern region. Meanwhile, the changes of the susceptibility assessment are fluctuating with the volatility of changing heavy rainfall. According to the geological hazard warning area after the heavy rainfall, the spatiotemporal LSMs are then validated using the ROC curve to determine the optimal stage used for spatiotemporal LSA. The ROC curves and AUC of the various stages are generated in Figure 12.

The above ROC curves indicate that the rainfall data in the development stage and peak stage show better applicability of the LSMs, and the spatiotemporal LSM result using the effective rainfall of the peak stage is the most accurate. The spatiotemporal LSM in the peak stage can be considered as the optimal spatiotemporal assessment result incorporating the effects of heavy rainfall.

5. Discussion

In previous studies, the spatial LSA has been investigated sufficiently, and the spatial LSMs are optimized from various aspects, such as the evaluation units, evaluation indices and evaluation models [6,51,52]. The spatial models perform well in evaluating spatial LSMs using the condition factors and historical landslides, but the temporal factors, such as climate change and heavy rainfall, may rapidly reshape the environment in the study area and thus increase the incidence of landslides. This further increases the landslide risk in the study area and reshapes the distribution of landslide susceptibility levels [53,54].

The present study considers the impact of condition factors and temporal factors based on the variability and trends of heavy rainfall in the study area. Although this study analyzes the fluctuating variation of the spatiotemporal LSA which is affected by volatility changes of heavy rainfall, the analysis may only be applied to the area where heavy rainfall is the main landslide-inducing factor. On the other hand, the spatial LSA in the previous study is validated using the historical landslide data, because they considered that the areas under the same spatial conditions as the landslide area are more prone to landslides. However, real-time rainfall data need to be employed to analyze the temporal LSA, and it is thus inaccurate to validate the accuracy of spatiotemporal LSA using historical landslide data. The future landslide data or landslide warning area can be used for the spatiotemporal LSA validation, however, it is difficult to implement if there is no landslide warning map because it is hard to collect the future landslide data.

6. Conclusions

In heavy rainfall areas, rainfall-induced landslides are always an important threat. The present study thus proposes a methodology to generate a spatiotemporal LSM, and the proposed method is applied to Kitakyushu, Fukuoka, Japan, under the impact of heavy rainfall in August 2021. The condition factors are selected as evaluation indices using the ${\chi }^2$ test, OOB error and multicollinearity test. The spatial LSMs are then generated using the RF model, DBN model and SVM model, and the results are compared and verified using ROC curves, AUCs and SCAIs to determine that the RF model is the optimal model for the spatial LSM in the study area. According to the rainfall variability and trends, the heavy rainfall period is divided into four stages, and the RI maps of various stages are obtained using ERM. The spatiotemporal LSMs under the influence of heavy rainfall in August are generated by coupling the four RI maps with the landslide probability map obtained using the RF model, and the changes of the susceptibility assessment are fluctuating over time. The spatiotemporal LSMs of various stages are verified using the landslide warning map, which illustrates that the four spatiotemporal results are accurate and the evaluation performance of spatiotemporal LSMs with the spatial LSM obtained using the RF model and RI map of the peak stage is the highest. The proposed method is scientific and feasible in spatiotemporal LSA and useful for areas with rainfall-induced landslides threatening human lives and property.