Landslide susceptibility prediction mapping with advanced ensemble models: Son La province, Vietnam

Landslide is a severe geohazard in many mountainous areas of Vietnam during the rainy season. They directly threaten human lives and properties every year. Landslide susceptibility maps are useful tools for risk mitigation, land-use planning, and early warning systems for local areas. It is necessary to update these maps continuously because of the complexity of landslide events. This fact requires further extending the approach techniques with practical implications. Therefore, this study aimed to develop landslide susceptibility prediction maps based on advanced machine learning (ML) techniques. Five state-of-the-art hybrid ML models were developed: bagging MLP, dagging MLP, decorate MLP, rotation forest MLP, and random subspace MLP with multilayer perceptron (MLP) as a base classifier. Sixteen causative factors were collected to build landslide susceptibility maps based on the relationship between historical landslide locations and specific local geo-environmental conditions. The model performance was verified using various statistical indexes. Based on the area under ROC curve (AUC) analysis results of the testing dataset, the rotation forest MLP model has the greatest predictive accuracy of AUC = 0.818. It is followed by the decorate MLP and bagging MLP (AUC = 0.804), the random subspace MLP model (AUC = 0.796), the dagging MLP (AUC = 0.789), and the single MLP (AUC = 0.698). The results of this study can be applied effectively to other mountainous regions to mitigate the risk of landslides.


Introduction
Landslides are among the most dangerous natural hazards in mountainous areas, generating more deaths and more economic damage than others (Wu et al. 1996;Zhu and Huang 2006). Landslides can wreak havoc on buildings, roads, and road infrastructure, damage agriculture, and threaten lives (Guzzetti et al. 2004;Harp and Jibson 1996).

3
The factors of precipitation, topographic, geological, geomorphic, and geoenvironmental conditions are closely related to landslide occurrences (Hang et al. 2021). Landslides may also be influenced by human activities such as deforestation, land-use change, transportation development, and plant degradation (Skilodimou et al. 2018).
In Vietnam, landslides often occur in mountain areas during the rainy season (Nguyen et al. 2020). Landslide susceptibility mapping can provide useful information for government agencies to implement land use and mitigation measure planning . Moreover, landslide spatial forecasting can significantly help managers determine the trend for landslide occurrences in a specific region (Aleotti and Chowdhury 1999). The susceptible landslide zoning is essential in identifying the areas where proper planning and infrastructure design could be developed and conducted considering landslide risks (Hang et al. 2021). In addition, the information on landslideprone sites needs to be updated continuously because of their complex geological processes (Chen and Li 2020;Van Westen et al. 2008).
Various methodologies and procedures have been utilized to map landslide susceptibility in hilly areas. Some statistical approaches used to construct landslide susceptibility maps, such as weights of evidence (Armaş 2012;Mahdadi et al. 2018;Pradhan et al. 2010), frequency ratio (Juliev et al. 2019), multicriteria decision analysis (Feizizadeh and Blaschke 2013), or a data-driven approach (Riaz et al. 2018). Recently, machine learning (ML) techniques have been explored and applied as a potential approach in landslide prediction studies based on geo-environmental factors (Bui et al. 2019). For example, landslide susceptible areas were forecasted based on logistic regression (Hemasinghe et al. 2018;Hong et al. 2016), artificial neural networks, support vector machine (Huang and Zhao 2018;Lee et al. 2017), random forest (Kim et al. 2018;Taalab et al. 2018), and naïve Bayes tree (Chen et al. 2018;Hu et al. 2021). Many advanced hybrid ML models have been constructed to improve the accuracy of predictive models (Achour and Pourghasemi 2020). Along with data availability, the selection of efficient predictive models is very important to receive reliable landslide susceptibility maps (Achour and Pourghasemi 2020;Pourghasemi et al. 2020).
Multilayer perceptron (MLP) as a base classifier method is one of the most popular artificial neural network techniques since it can bring a higher predictive accuracy for predictive models (Adnan et al. 2020;Hong et al. 2020;Li et al. 2019;Meghanadh et al. 2022;Zare et al. 2013). MLP can be applied as a single model in creating a landslide susceptibility map (Adnan et al. 2020;Meghanadh et al. 2022;Zare et al. 2013). MLP can also be combined with other ML models to forecast the spatial distribution of landslides, such as the particle swarm optimization model ) and the stochastic gradient descent model . These studies concluded that the ensemble ML models with MLP have higher predictive accuracy than other ML models (Achour and Pourghasemi 2020;Wang et al. 2020).
During the rainy season, large landslides frequently threaten human lives and properties in Son La province, which is located in Vietnam's northwest mountainous region (Ahlheim et al. 2008). However, it still lacks studies on evaluating landslide susceptibility for the area. This study aimed to create the landslide susceptibility map for the new case study of Son La province, using advanced ensemble ML models. We developed five hybrid models which combine complex ML models with the MLP model as a base classifier model. They are bagging MLP, dagging LP, decorate MLP, rotation forest MLP, and random subspace MLP. The landslide-prone areas are identified based on the research area's specific topographical, hydrological, geological, and environmental conditions. The resulting susceptibility map can supply useful information for local authorities in strengthening landslide disaster risk mitigation strategies at the regional scale.

Study area
Son La is a mountainous province in Vietnam's northwest region. It has 250 km of borderline with the Lao People's Democratic Republic. The province covers an area of 14,174 km 2 and is situated within 20° 39ʹ and 22° 02ʹ North latitude and between 103° 11ʹ and 105° 02ʹ East longitude (Fig. 1). The population of Son La was 1,252,650 people in 2019, of which 13.85% are in urban areas and 86.15% in rural areas. The Son La province has 12 ethnic groups, of which 53.2% Thai group, 17.6% Kinh group, 14.6% H'Mong group, 7.6% Muong group, and 7.2% other ethnic minority groups (available at https:// www. sonla. gov. vn/). This province has many mountains, hills, and rivers and is surrounded by primitive forests. In particular, the terrain of Son La province is strongly divided by the Da River, Ma River, and high mountain ranges (Thach and Canh 2011).
In the study area, landslides are influenced by extreme weather, prolonged heavy rainfall, and human activities. (Van Hoang et al. 2019). The rainfall averages 300-700 mm per month during the rainy season, with maximum daily rainfall exceeding 100 mm  (Hoang and Tien Bui 2018). In addition, the other factors of topography, geology, fault density, land use, weathering crush, type of soil, and soil thickness also affect landslide susceptibility in this province (Thach and Canh 2011). Landslide hazard has caused serious damage to humans, houses, road networks, agriculture, road infrastructure, and drainage system in the province (IFRC 2021). For example, there were 10 killed people, 4 injured, and 258 damaged houses by landslides in the Muong La district in 2017 (Dung et al. 2021).

Methodology flowchart
The methodological framework of this study is shown in Fig. 2. The framework consists of six main steps: (1) data collection and preparation, (2) checking the multicollinearity of the landslide-related variables, (3) landslide inventory mapping, (4) landslide modeling process using a hybrid ML ensemble framework, (5) landslide susceptibility mapping, and (6) model comparison and validation. First, the topographical, geological, hydrological, environmental data and historical landslide locations were collected and prepared for modeling samples. Second, these data were checked for multicollinearity to avoid computational instability in model assessment. Third, the landslide inventory map was created. Using the Sample tool in ArcGIS Pro, the sample data were randomly split into training (70%) and validating (30%) datasets. Fourth, ensemble ML techniques were developed for the modeling, including bagging MLP, dagging MLP, decorate MLP, rotation forest MLP, and random subspace MLP. Fifth, the above hybrid ML models were used to create landslide susceptibility maps. These susceptibility maps were divided into five susceptibility classes: very low, low, moderate, high, and very high susceptibility, based on the Natural Break technique in ArcGIS Pro. Finally, the susceptibility models were compared and verified using the cross-validation approach.

Landslide inventory map
It is necessary to rely on detailed information on previous landslide events in studies on the formation mechanism of landslides, landslide susceptibility mapping, and developing landslide risk mitigation strategies (Mirus et al. 2020). Thus, landslide inventory mapping is often the first and the most crucial step in data preparation for modeling (Bui et al. 2019). In this study, 1,689 landslide locations were collected from many different sources, in which 1,225 landslide positions were explored from the website of the Institute of Geosciences and Minerals of Vietnam (available at http:// canhb aotru otlo. vn/ phanv ungca ctinh. html), and 464 landslide locations were added through field survey combined with the interpretation of Google Earth images.
There are three main causes leading to landslide formation in Son La province, including construction structures on steep slopes, thick weathering crust, and high groundwater levels (Thach and Canh 2011). In addition, heavy prolonged rainfall events are considered the primary cause of landslides in the study area (Hoang and Tien Bui 2018). The statistics have recorded many landslide developments in Quynh Nhai, Sop Cop, Song Ma, Moc Chau, Van Ho, and Thuan Chau districts (Hoang and Tien Bui 2018;IFRC 2021;Thach and Canh 2011). The landslide curves stretched from 10 to 100 m, and the landslide regions were bigger than 30 m 2 (Thach and Canh 2011).
The landslide inventory data, which are historical landslide positions, need to be randomly divided into two parts for training and validating predictive models. The training and testing datasets ratio does not need to be a fixed number; it depends on the predictive methods and the number of landslide locations within the study area. For example, it could be 80/20, 70/30, 60/40, or 50/50 (Chen and Chen 2021). Most analyses used the ratio of 70/30 to divide the samples for hazard susceptibility modeling (Huang and Zhao 2018). Thus, in this study, the landslide inventory data were randomly split into two samples: 70% training dataset (1183 landslide sites) and 30% validating dataset (506 landslide sites) using the Sample tool of ArcGIS Pro software.

Landslide causative factors
Modeling work begins with the identification of landslide causative factors. The factors are referred to in previous studies and are based on the available data in the research area (Kavzoglu et al. 2019). Precipitation, topography, hydrology, geology, geomorphology, and geo-environment are factors that significantly impact landslide formation in mountain areas (Bui et al. 2019;Meghanadh et al. 2022). In this study, 16 landslide causative factors were taken for the modeling in the following subsections ( Fig. 3 and Table 1). We used frequency ratio analysis to evaluate the relationship between the landslide sites and the causative factors. The frequency ratio is the ratio between the number of landslide sites in each classification divided by the total number of landslide positions and the number of pixels in each classification divided by the total number of pixels (Khan et al. 2019;Pham et al. 2016). The analysis result in Fig. 4 shows the percentage of landslide locations in each classification of each factor. In all causative factor maps, the higher frequency ratio value of the classification correlates with the higher probability of landslide occurrences and vice versa. The trend indicates that such a classification is reasonable.

Elevation
In investigations of landslide susceptibility mapping, elevation is a crucial factor in studies on landslide susceptibility mapping (Goetz et al. 2015;Tehrany and Kumar 2018). The higher-elevation zones correspond to the higher landslide frequencies (Myronidis et al. 2016). The elevation map was constructed using the digital elevation model (DEM). The DEM was created using an ALOS picture with a spatial resolution of 30 m obtained from https:// www. eorc. jaxa. jp/ ALOS/ en/ aw3d30/ data/ index. htm in March 2021. The research area's elevation ranges from 70 to 2884 m and is separated into nine groups (Fig. 3a).

Slope
The slope is a critical factor in landslide susceptibility evaluations since it can control the landslide creation and movement in tropical mountainous areas (Dai et al. 2002;Guns and Vanacker 2013). In terms of magnitude and frequency, the steep slope is closely associated with landslide development (Aditian et al. 2018). However, moderate slopes also cause landslides, particularly in colluvial or deeply weathered slopes (Ghimire 2011). In this study, the slope angle map was created using a DEM with a spatial resolution of 30 m in ArcGIS Pro software. The map is separated into six classes, ranging from 0° to 78.3° (Fig. 3b).

Aspect, curvature
Other topographical characteristics, aspect and curvature, are often used as primary input data in landslide prediction (Arabameri et al. 2020;Hang et al. 2021;Kavzoglu et al. 2019). The aspect factor can affect the formation of shallow landslides on hill areas (Shao et al. 2021). In particular, northerly aspects (Northeast to Northwest) often impact landslide occurrences more than in other directions (Douglas et al. 2013). However, the influence of the aspect factor on landslide formation can change due to the vegetation cover type (Karsli The curvature factor reflects the morphology of the topography (Lee and Talib 2005). The curvature indicates the convex surface, concave surface, or flat surface of the topographical morphology (corresponding to positive, negative, or zero values), which is related to water holding capacity after prolonged heavy rainfall (Al-Najjar et al. 2019). The surfaces are more convex or more concave, and the probability of landslide formation is higher and higher (Ohlmacher 2007). The aspect map (Fig. 3c) was divided into nine classes, and the curvature map ( Fig. 3d) was categorized into five levels.

Elevation difference
Elevation difference represents the altitude difference of all points on the Earth's surface (Corsini et al. 2005). This factor also reflects the exiguous separation of elevations where water redistribution is very significant to landslide formation (Van Westen et al. 2003). The elevation difference is often associated with the number and the impacted area of the landslide events (Lineback Gritzner et al. 2001). In particular, if the elevation difference increases, the mass movement and the erosion due to potential energy will be higher (Ghimire 2011). In this study, the elevation difference factor on the 1:50,000 topographical map was derived using relative altitude (meters) in each square grid (1 km 2 ). It was categorized into nine classes (Fig. 3e).

TWI
The Topographic Wetness Index (TWI) measures the impact of topography on the location and amount of saturated runoff source zones (Pourghasemi et al. 2012). TWI was built from a 30 m DEM and was divided into seven categories in ArcGIS Pro software Fig. 3f. TWI is calculated as follows: where A S denotes the specific basin area (m 2 /m) and β denotes a sloped angle in degrees.

Rainfall
Heavy prolonged rainfall is considered the primary cause of landslides in mountain areas (Singh et al. 2021). It might trigger unexpected landslides depending on the topographical and geological characteristics of the ground/rock mass (Dung et al. 2021). In this study, the rainfall data were derived from 25 rain gauge stations in Son La province and neighboring provinces from 2010 to 2021. The rainfall map was created by interpolating the rainfall distribution based on DEM using the geostatistical Kriging technique (Fig. 3h).

Drainage density
Drainage density describes the drainage availability in the short term in response to changes in environmental conditions (Mezughi et al. 2011). Drainage density is directly associated with landslide formation in mountain places (Arabameri et al. 2020). The drainage density is calculated by dividing the total drainage length (km) in each square grid by (1) TWI = ln A S tan the number of square grids (Fig. 3i). The drainage density map of the Son La province was provided by the provincial Department of Natural Resources and Environment.

Road density and distance from the road
The road network is often associated with an increase in landslide events (Skilodimou et al. 2018). Meanwhile, road density is often used to measure the effect of development on landslide formation and distribution (Simon et al. 2015). Distance from the road represents a negative relation with landslide events (Akgün and Bulut 2007). The shorter road distances are, the higher the landslide occurrences are (Skilodimou et al. 2018). The road density map was grouped into five levels (Fig. 3l), and the distance from the road map was classified into six classes (Fig. 3k). The road network map within the study area was extracted from the national road network map provided by the Department of Survey, Mapping, and Geographic Information.

Distance from the rivers
Many previous studies have proved that landslide has a close relationship with distance from the rivers (Arabameri et al. 2020;Bui et al. 2019). Landslides often happen along the sides of the valleys where the groundwater flows toward rivers and streams (Raja et al. 2017). In this study, the distance from the river map was calculated in ArcGIS Pro software and consisted of seven classes (Fig. 3n).

Hydrogeology
Hydrology and hydrogeology significantly affect landslide formation in the hilly areas (Kayastha et al. 2012). Many studies have looked at the importance and complexity of hydrogeology in landslide susceptibility evaluations (Frodella et al. 2021;Sujatha 2021). The hydrogeology map of Son La province was obtained on a scale of 50,000 in 2020 from the Vietnamese Ministry of Natural Resources and Environment. The map is covered by four hydrogeological units (Fig. 3m).

NDVI
The Normalized Difference Vegetation Index (NDVI) measures the development of vegetation on the Earth's surface (Jaafari et al. 2014). The NDVI explains the link between vegetation density and the occurrence and distribution of landslides (W. Chen et al. 2019). The NDVI is expressed as below: where NIR denotes the infrared reflectance of the electromagnetic spectrum, and Red denotes the red reflectance of the electromagnetic spectrum.
The NDVI value in this research ranged from -0.643 to 0.694 and was separated into five groups (Fig. 3g).

Geology and geomorphology
Although most of the geological and geomorphological factors change over relatively long periods, their characteristics have a substantial influence on the evolution of erosional and landslide processes in mountainous areas (Bui et al. 2019;Pisano et al. 2017). In 2020, the Vietnamese Ministry of Natural Resources and Environment released geological and geomorphological maps of Son La province on a scale of 50,000. The geological map is covered by nine geological groups (Fig. 3s). The detail of the nine geological groups is shown in Table 1. The geomorphological map is covered by fifteen geomorphological types (Fig. 3p).

Land cover
Anthropogenic activities can change land cover due to the transformation of the natural landscape (Promper et al. 2014). Therefore, this factor is often considered an important triggering factor for landslide occurrences in mountain regions (Van Westen et al. 2003). Land cover can affect the landslide frequency and distribution quickly (Hong et al. 2007). In this study, the land-cover map was developed by ESRI in 2020 using deep learning models and satellite images (downloaded at https:// livin gatlas. arcgis. com/ landc over/). The land-cover map of the study area consists of eight groups, including bare ground (0.02%), built area (2.57%), crops (6.59%), flooded vegetation (0.01%), grass (0.53%), scrub/shrub (28.86%), trees (59.55%), and water bodies (1.87%) (Fig. 3m).

Multilayer perceptron (MLP)
The MLP method is one of the most common artificial neural network approaches and is widely applied in landslide susceptibility assessment (Gómez and Kavzoglu 2005;Li et al. 2019). The backpropagation algorithm is the training rule of MLP (Gómez and Kavzoglu 2005). The goal function's minimal value and the method's optimal weight values may be modified and calculated . This model has three main components: input, hidden, and output layers. Landslide influencing factors are considered the input layers. The output layers are made up of the categorized findings that are used to classify landslides and non-landslides. The hidden layers are applied to transform inputs into outputs (Gómez and Kavzoglu 2005). The number of neurons in the input layer of X = x 1 , x 2 , … , x m0 , hidden layer, and output layer are designated as m 0 , m 1 , and m 2 , respectively, if the MLP model comprises multi-input variables and multioutput variables. The input and output of neurons in the hidden layer are calculated as follows: where h j , θ j , and z j denote the input, the threshold, and the output of the jth neuron in the hidden layer; w ij denotes the weight value of the ith input and the jth in the hidden layer neurons; and f h j denotes the activating function. Afterward, the inputs and outputs of output layer neurons are represented as below: where h k , θ k , and z k denote the kth input, the threshold, and the output in the output layer neurons; and w jk denotes the weight value between the jth hidden layer neuron and the kth output layer neuron.

Bagging (BA)
Bootstrap aggregating or Bagging algorithm was proposed by Breiman (1996). This algorithm was applied to achieve an aggregated predictor based on different bootstrap samples (Breiman 1996). This algorithm uses a training dataset T (x k ,y k ), where x k ∈ Q, y k ∈ (landslide; non-landslide), k = 1-M, and M is the number of bootstrap samples. Next, a bootstrap sample T k is generated from the initial training dataset based on the replacement method. Therefore, the model is formed through a base classifier B k using the bootstrap sample T k and a classifier B k (x) is developed from each bootstrap sample T k . Finally, the classifier B* is synthesized from B 1 , B 2 , …, B n and calculated as below (Bauer and Kohavi 1999): where B i (x) represents a classifier that is generated from each bootstrap sample T k .

Dagging (DA)
Dagging was first proposed by Ting and Witten (1997). This method determined the final prediction based on the majority vote. In this technique, the training dataset was divided into many separate classified parts, and each part of the data corresponded with a basic learner (Ting and Witten 1997). If the input training dataset D has K samples, the dagging algorithm creates N datasets from the input training dataset. Each dataset consists of k samples (kN < k), and other datasets do not contain a similar sample. After that, each dataset would be trained by a basic classifier to build a classification model. Thus, N classification models can be formed from the N original datasets. These models make their prediction classes based on the given query samples. Finally, the prediction class of the dagging method has the most votes (W. Chen and Li 2020).

Decorate (DE)
The decorate was developed by Melville and Mooney (2003). Decorate is an ensemble meta-learner to create Diverse Ensemble Creation by Oppositional Relabeling of Artificial Training Examples (DECORATE). This algorithm is applied to create a new classifier based on the combination of original training datasets. The content of the decorate algorithm can be understood as follows (Melville and Mooney 2003): 1. I n p u t a t r a i n i n g d a t a s e t D = x 1 , y 1 , x 2 , y 2 , … , x p , y p w i t h x i ∈ R b , y i ∈ Y = l 1 , l 2 , … , l q , i = 1, 2, … , p;j = 1, 2, … , q ; in which p denotes the number of training sub-datasets (p > 1), b denotes the number of attributes in each subdataset (b > 1), m denotes the number of class labels in a classification model (q > 1). 2. The base learning algorithm-BaseLearn-is used to train a classifier C 1 from the original input dataset D, and the first ensemble C * = C 1 is got. 3. Decorate algorithm generates classifiers iteratively.

Rotation forest (RF)
Rotation forest (RF) was an ensemble learning method proposed by Rodríguez et al. (2006). This method developed an ensemble of decision trees based on the random subspace and bagging methodology with principal component analysis (PCA). In this technique, the input training datasets were divided into many sub-datasets to create the classifiers, and PCA was used for each sub-dataset (Kuncheva and Rodríguez 2007).
The input training dataset is E, the class label of the input training dataset is F, and the feature set is G. If the quantity of training times is T with t features, then E is T × t matrix. Assuming that class labels F are formed based on the class set of H = (h 1 , h 2 , … , h j ) ; the feature set is allocated into M sub-datasets, P-decision tree in a forest rotation of D = (D 1 , D 2 , … , D P ) . According to that, two indices-M and P-need to be pre-calculated. The procedure to determine the input training dataset consists of four steps (Sahin et al. 2020): (1) Divide the feature set G into M feature subset, and the attribute of each feature subset of N = n/M; (2) G i,j is the attribute subset to train the classifier D i , and Y i,j is the dataset in G i,j ; (3) Construct a rotation matrix R i with the nominated ratio in the matrix F i,j . The coefficients of the matrix F i,j include: ij , that is determined based on a linear transformation; (4) The rotation matrix R i can be calculated by: In addition, the classification step can be evaluated for a given case y; then, the classifier Q i = q ij yR b i will be used to classify the probability of this case. In this way, the confidence of a class can be determined as the following equation: Finally, y is assigned to a class which has the largest confidence determined.

Random Subspace (RSS)
The random subspace (RSS) was first proposed by Ho (1998). The training dataset was created in the modified feature space in the RSS method to build a higher number of training variables (Ho 1998). The RSS algorithm can be expressed as follows: input a training dataset. D(x i ,y i ); where x i ∈ T and is a m-dimensional vector x i = (x i1 , x i2 , …, x im ), y i ∈ (landslide; non-landslide). First, m * features are randomly selected from the training dataset, where m * < m. In this way, the m *-dimensional random subspace of the original m-dimensional feature space is generated. Second, the modified training dataset D* consists of m * -dimensional training features x i = (x i1 , x i2 , …, x im* ). Finally, a final classifier is developed based on the combination of primary classifiers according to a voting scheme (Lai et al. 2006).

Model validation and comparison
It is required to validate the prediction models to assess the applicable ability. We used a variety of statistical metrics to assess the performance of five suggested ensemble models in this research. They are positive predictive value (PPV), negative predictive value (NPV), sensitivity, specificity, accuracy (ACC), F-measure, Jaccard, mean squared error (MSE), root mean squared error (RMSE), receiver operating characteristic (ROC) curve, and area under the ROC curve (AUC). Sensitivity and specificity denote the true positive rate and true negative rate of landslide and non-landslide locations (Bui et al. 2019). Accuracy is the mathematics average of sensitivity and specificity or the scale between the number of landslide and non-landslide pixels (Arabameri et al. 2020). F-measure is expressed as a weighted harmonious average of the accuracy and revocation in binary classification (Ha et al. 2021). Jaccard is the Jaccard sameness coefficient factor (Pham et al. 2017). The ACC, kappa, MSE, and RMSE statistical indicators are determined as in the following equations: in which TP denotes true positive, TN denotes true negative, FP denotes false positive, and FN denotes false negative.
RMSE represents the level of dispersion of predictive values from actual values (Nsengiyumva and Valentino 2020). The less the RMSE value is, the higher the accuracy of the prediction model is (Bui et al. 2019). Meanwhile, the AUC indicator is considered the main value in assessing the performance of the predictive models (Ye et al. 2016).

Multicollinearity and factor selection
Variance inflation factors (VIF) and tolerance parameters were used for checking the multicollinearity of the landslide causative factors in order to select optimal input factors. The tolerance and VIF determine the transformation in the standard errors of the attended factors (Bui et al. 2019).). The multicollinearity problem of used factors will happen if VIF is greater than 10 or tolerance is less than 0.1 (Lin et al. 2017). Figure 5 shows the results of VIF and tolerance calculated for sixteen landslide influencing factors. All factors have VIF > 10 and tolerance < 0.1; therefore, they can be applied to build landslide susceptibility modeling in this research.

Validation and comparison of the susceptibility models
In this study, five advanced hybrid ML models (BAMLP, DAMLP, DEMLP, RFMLP, and RSSMLP) and the single model (MLP) were applied to generate landslide susceptibility models. The cross-validation method was applied to validate these used models. Several standard quantitative indexes were tested on training and validating datasets (Table 2 and Fig. 6).

Landslide susceptibility mapping
Five advanced hybrid ML models and a single MLP model were used to create the landslide susceptibility maps based on the training dataset. These resulting maps were RSSMLP, and f MLP divided into five susceptibility classes: very low, low, moderate, high, and very high susceptibility (Fig. 8).
Based on the statistical indicators in Table 2, the RFMLP model is the best prediction model compared to the other models (DEMLP, BAMLP, DAMLP, RSSMLP, and MLP). As a result, the RFMLP model was chosen to create the most accurate landslide susceptibility map for the research region. The final results represented 217,377 ha (15.45%) in very high-susceptibility areas; 303,376 ha (21.56%) in high-susceptibility areas, 271,515 ha (19.30%) in moderate-susceptibility areas, 215,139 ha (15.29%) in low-susceptibility areas, and 399,679 ha (28.40%) in very low-susceptibility areas (Fig. 9).

Discussion
Landslides are dangerous and devastating geological procedures in mountainous regions (Šilhán 2020). They have caused thousands of deaths and injuries annually and are the main natural hazards that cause loss of life in Asia (Highland and Bobrowsky 2008). Therefore, the information on landslide susceptible areas needs to be updated continuously because of their complex geological processes (Chen and Li 2020;Van Westen et al. 2008). The present study provided a feasible approach for susceptible landslide evaluation using Fig. 9 Areas of landslide susceptibility and percentages by classification advanced ensemble ML models. The five advanced ensemble ML models were also developed to estimate predictive landslide susceptibility maps.
ML approaches have been considered an effective and potential method in natural hazard assessment because of their flexibility and predictive accuracy (Bui et al. 2019). Hybrid ML models, in particular, do not rely on statistical assumptions and may quantify the relevance and effect of landslide-related factors (Achour and Pourghasemi 2020). Previously, single models, such as the MLP model, were often applied to build landslide susceptibility maps (Adnan et al. 2020;Meghanadh et al. 2022;Zare et al. 2013). Table 2 and Fig. 7 also show that ensemble ML models had better predictive accuracy than a single ML model. Several studies have currently developed hybrid models that combine MLP and other ML models; these hybrid models showed better predictive accuracy. Li et al. (2019) employed an ensemble of PSO-MLP to create a landslide susceptibility map for Shicheng County in China with the AUC = 0.881. In another study, Hong et al. (2020) applied an ensemble of MLP-SGD to predict the landslide susceptible areas for Yanshan County in Jiangxi province in China with the AUC = 0.822. In this study, we explored advanced hybrid ML models, namely BAMLP, DAMLP, DEMLP, RFMLP, and RSSMLP, with the MLP as a base classifier model. The validation results showed that the RFMLP model has the highest predictive accuracy of AUC = 0.818. Compared with the previous studies, the RFMLP ensemble model with the case study of Son La province was found to be appropriate. This study also provided a new and efficient hybrid model for assessing landslide susceptibility in mountainous areas.
Son La is situated in the center of the northwest area of Vietnam, and extensive landslides often destroy many regions during the rainy season (Ahlheim et al. 2008). The increasing landslide events have threatened human lives and properties in the study area (IFRC 2021). The landslide susceptibility maps can supply relevant information for government agencies and local authorities to implement planning, land-use sustainable management, or build early warning systems (Bălteanu et al. 2020). The result of the susceptibility map may supply helpful information for analyzing and evaluating the landslide risk with practical meanings at the local scale. This present study generally added a new approach to assessing landslide susceptibility based on ML techniques.

Conclusion
We used state-of-the-art hybrid ML models to create landslide susceptibility prediction maps in this work. Historical landslide placements, as well as precise geo-environmental, meteorological, hydrogeological, geological, and topographical data, are included in the input data. Several statistical indices were utilized to verify and compare the forecast models. It can be seen that the predictive accuracy values of the used models are good. The RFMLP model has the best AUC (0.818), followed by the DEMLP and BAMLP (0.804), the RSSMLP (0.796), the DAMLP (0.789), and the MLP (0.698). The RFMLP ensemble model was selected to build the landslide susceptibility prediction map for the research region. This resulting map can supply useful information for policy-makers and decisionmakers in preventing landslide risks in the future. This research also supports the use of ML approach to analyze and manage natural risks in mountainous locations.