The landslide inventory is an essential commission for GIS-based landslide susceptibility modeling (Titti et al., 2021). The quality and quantity of landslide locations have an impact on the outcome and accuracy assessment of the predictive model. Satellite images are valuable sources that support the landslide inventory, especially in mountainous regions that are difficult to access. However, dense tropical forests may present challenges and uncertainty in detecting landslide occurrences based on these remotely sensed data. Hence, taking advantage of satellite images for landslide remote detection and ground true validation is essential for enhancing the truthfulness of the landslide inventory database. In the present study, we identified 271 landslide points over the examined area, including historical destructive geological events and field investigation. The points were first sampled based on numerous geospatial resources, including high-resolution Earth Engine images, optical satellite images from Sentinel-2, and thematic maps, which were then verified in field campaigns during the rainstorm season in 2022. As illustrated in Figure 1, the dataset was split into 70% (196 points) for the training phase and 30% for validation (71 points).

In the domain of GIS-based landslide susceptibility, preparing conditioning factors in a consistent spatial database required outstanding knowledge of not only geographic accuracy but also of the various types of features to fully reflect the interest domain. To date, to select landslide-related model variables, it should be noted that there are no standard guidelines or regulations for choosing optimal factors, as these are well-defined based on the specific case and research areas. Among a plethora of causative factors relating to landslide susceptibility, we recognized 15 significant factors through state-of-the-art literature reviews (Pourghasemi et al., 2018; Reichenbach et al., 2018; Yong et al., 2022). We categorized all factors into five fundamental groups, 1) topography, 2) hydrology, 3) geology, 4) land cover conditions and anthropogenic activities, and 5) weather, as an effort to exhaustively consider the contribution of related datasets to the disaster. Continuous data layers were reclassified using Jenks natural breaks optimization, while discrete data layers were stored in unique values for further analyses. For the post-processing, the final raster layers were geometrically corrected according to the World Geodetic System (WGS), 1984, and the UTM (Universal Transverse Mercator) Zone 48°N (North), which then transformed into a 10-m resolution. Figure 3 displays the maps of all 15 factors after standardization, as mentioned previously. Detailed descriptions of each group-based conditioning factor are provided in the following sections.

The bivariate statistical FR method was applied to estimate the landslide densities in all subclasses in each factor. The goal of this method is to compute the percentage of landslide pixels located in subcategories of all factors, correcting the associated raters to the propensity of landslide occurrences (Lee and Talib, 2005). To perform the estimation, the landslide inventory training dataset and factor maps were used for obtaining the following equations: F R i , j = L S i , j   / ∑ L S i , j N i , j   / ∑ N i , j , (3) R F i , j = F R i , j ∑ F R i , j × 100 , (4)

where i is the subclass of the considered factor j; LS_i,j is the number of landslide pixels in each class; N_i,j is the number of class pixels; FR_i,j is the frequency ratio; RF_i,j is the relative frequency (%).

The subjective approach, AHP, was first introduced in Saaty, 1977 (1980), known as a widely used multi-criteria decision analysis that allows decision-makers to prioritize and evaluate alternative options based on multiple criteria. In the field of critical landslide area assessment, this method provides a robust, yet simple to handle, complex decision-making problems. It is based on the principle of pairwise comparison, where the relative importance of a criterion is assessed in relation to others. Here, the AHP was implemented by a subjective hierarchical structure that contains 15 landslide causative factors. Based on judgments of multidisciplinary experts who collaborated with local officers, pairwise comparisons are then made to determine their relative importance. These comparisons are repeated at each level of the hierarchy until an acceptance consistency ratio is obtained, indicating that the criteria are consistent with others. The final step involves combining the pairwise comparisons to calculate the overall weight of individual factors. The metrics in the process, including Consistency Index (CI) and Consistency Ratio (CR), were estimated according to the following equation: C I = λ max − n n − 1 , (5) C R = C I R I , (6)

where λ max is the maximum eigenvalue of the matrix, n is the number of considered criteria (n = 15); the random consistency index (RI) is 1.59, which was used for 15 criteria (Saaty, 1980). The obtained CR less than 0.1 implies the consistency and acceptance of the decision-makers’ pairwise comparison matrix.

Initially proposed by Shannon (1948), SE is a well-established information theory-based method. The basic idea is to give more weight to events that have higher entropy, as they contain more uncertainty and information, and less weight to events with lower entropy. Mainly based on ready-to-use data, the SE formula measures the average amount of information contained in a dataset, and the weights are derived from these entropy values. In risk assessment, the entropy of a hazard event can be used to determine its likelihood and potential impact. Events with higher entropy carry more uncertainty and, therefore, have a higher risk associated with them. Here, we computed the entropy of landslide susceptibility based on the distribution of detected landslides in the training dataset to the contributing factors. SE weighting was estimated using the following equations (Roodposhti et al., 2016; Agrawal and Dixit, 2022): P i , j = F R i , j ÷ ∑ i = 1 m F R i j , (7) E j = − 1 log 2 m j × ∑ i = 1 m P i j log 2 P i j , (8) W S E j = 1 − E j ∑ j = 1 n 1 − E j , (9)

where i is the subclass of the considered factor j, m is the number of subclasses in each conditioning factor, n is the number of conditioning factors (n = 15), P_ij is the probability density, FR_ij is the frequency ratio, E_j is the entropy value, and W S E j is the entropy weight.

The synergy of both objective and subjective approaches may enhance the performance of the model and mitigate issues related to the ill-posedness. Therefore, the data-driven method SE and knowledge-based method AHP were integrated to derive the combined weights for each of 15 conditioning factors, as seen in the following equation (Wang and Zhang, 2018): W S E − A H P j = W S E j × W A H P j ∑ j = 1 n ( W S E j × W A H P j ) , (10) where W S E j and W A H P j are derived from the objective-based SE and the subjective-based AHP, respectively; j is the considered factor; n is the total number of conditioning factors of the landslide susceptibility model (n = 15).

To retrieve the landslide susceptibility index, a summation of the product of rating subclasses and respective weights of each conditioning factor is given by the following equation: L S I = ∑ j = 1 n W j × R j F R , (11)

where n is the number of landslide conditioning factors (n = 15); W_j is the factor weight obtained by the corresponding methods of AHP, SE, and AHP–SE; and R_j is the two-dimensional matrix of factor j that has been reassigned to the RF. The continuous range value of the landslide susceptibility index is divided into five probability categories, namely, very low, low, moderate, high, and very high, using the natural breaks Jenks function.

The performance of various models was evaluated using the receiver operating characteristic (ROC) curve. The primary statistical metrics employed to measure the model accuracy are true positive rate (TPR), false positive rate (FPR), and AUC. The TPR, also referred to as sensitivity, reflects the proportion of correctly classified positive cases among all positive cases, while the FPR, also known as specificity, measures the likelihood of a true negative case being classified as negative. The AUC spans from 0 to 1. A retrieval value of the AUC closer to 1 results in the better performance of the model: T P R = T P T P + F N = T P P , (12) F P R = T N T N + F P = T N N , (13) A U C = ∑ T P + ∑ T N P + N , (14)

where true positive (TP) and true negative (TN) are the values of correct landslide and non-landslide pixels, respectively; true negative and false positive (FP) are the values of incorrect landslide and non-landslide pixels, respectively; P and N are the corresponding total number of landslide and non-landslide pixels.

We also estimate the relative landslide density index to evaluate the performance of different models in each landslide-sensitive zone, using the following equation: R = n i N i   / ∑ n i N i × 100 , (15) where n_i and N_i are the numbers of landslide pixels and total pixels in each susceptibility class, respectively.

Table 2 showed the FR and RF calculated for subclasses in landslide-controlling factors. By considering the quantity relationship between training points (196 points) that appeared in each subcategory (total of 79,865,948 pixels), higher FR values indicate higher sensitivity of a class to landslide occurrence compared to the remaining classes within a predisposing factor. It should be noted that even if a particular subclass has a small number of landslide inventory points, the FR value may still be high. This can occur when the percentage of pixels in the subclass is low relative to the total number of pixels in the thematic layer. The calculated RF (%) values were then assigned to reclassify thematic maps of all 15 causative factors as input for the landslide susceptibility model.

Regarding topographical groups, the highest elevation ranging from 1,019 m to 1,971 m, slope degree of 17.5–26.8^o, east slope aspect, and convex plan curvature (>0.05) were determined as the most influential subclasses. In case of the group set by hydrological variables, distance to rive spanning 200–300 m, drainage density between 1.6 and 2.7 km²/km, and the range of the TWI of 0.6–4.7 become more sensitive compared to other corresponding subcategories. The soil of Epi Lithi Humic Acrisols (ACu-l1), wash slope and corroded slope, and 400–500 m proximity to the fault gave the most pressure in the geological group. The group of land cover conditions and anthropogenic activities witnessed the most significant susceptibility in subclasses of tree cover, distance to road of 400–500 m, and a range value of the NDVI, approximately 0.7 – 0.8. Last, weather-based components estimated the rainfall amount on average of 15.6–23.3 mm/day, and two climate patterns (IB1a and IIA2a) are the most impact subcategories contributing to the landslide probability.

Table 3 represented the AHP order matrix for 15 conditioning factors, with a range of positive integer values from 1 to 7 applied for comparison among these aforementioned criteria. The achieved metrics λ max , CI, and CR are equal to 15.64, 0.045, and 0.028, respectively, strongly confirmed for the consistency of the matrix. Relative weights for all factors in the method of the AHP, SE, and AHP–SE were determined and illustrated in Figure 4. From the charts, subjective judgment strictly obeys the AHP process, with the most important causative factors being precipitation and slope. In contrast, data-based statistics SE accounting for soil and climate affect much of the landslide susceptibility model. The final ranking, which incorporated both AHP and SE, identified slope and soil as the most crucial factors for detection of landslide-prone regions. When considering all three methods, soil, climate, slope, and precipitation were recognized as the top important factors, with significantly higher weights than those of the rest. For example, the lowest weight value for slope in SE was 0.138, while the highest weight value for climate in AHP–SE was 0.194. The plan curvature was found to contribute the least to the model, with the lowest weight assigned to it in all three methods (AHP: 0.009, SE: 0.021, and AHP–SE: 0.002). Notably, AHP–SE showed a significant difference between the most crucial factor, climate, and the least important one, plan curvature, with weight values ranging from 0.002 to 0.194.

Figure 5 illustrated the distribution of landslide sensitivity using the AHP, SE, and the synergy of these two methods. The maps clearly demonstrated that the most sensitive landslides are located along the Truong Son mountain range in Vietnam, stretching from the northwest to the southeast. In contrast, coastal areas showed a very low to low proportion of landslide occurrences. The polar chart (Figure 5D) summarized the statistics area of the five levels of landslide susceptibility based on AHP, SE, and the integrated approach. Overall, the major landslide-prone areas belong to the high level in all three methods, on average, accounting for nearly a quarter of the total study area. In particular, the AHP method showed the highest percentage (29.83%) of the high landslide-prone category, while the SE method exhibited slightly different proportions for high and very high levels, i.e., 28.76% and 26.05%, respectively. The combined approach AHP–SE resulted in 28.97% for high risk and 21.93% for very high risk zones. The figures for moderate risk of landslide trigger in the three methods appeared to be similar, with negligible differences in area percentage, ranging from the lowest 21.32% (AHP) to the highest 22.70% (SE).

Looking into the classes of low and very low susceptibility, it is evident that these classes resulted in the smallest percentage of the susceptibility area compared to the three remaining categories. In the SE, approximately 11% of the natural area was identified as having the least landslide probability, and the number remains unchanged for the level of low sensitivity. Using the AHP model, the proportions of very low and low levels were found to be 11.71% and 15.65% of the total area, respectively, whereas 14% and 13.30% of the area are shown in very low and low landslide susceptibility degrees according to the combined method AHP-SE. In comparison to categories of the high landslide sensitivity class, the statistical data for the level of very low susceptibility revealed significantly lower values, i.e., only one-third of the values obtained through the AHP and the SE method and half the number of the high landslide sensitivity class identified through the AHP–SE method. By comparing these three methods, the statistics strongly suggest that the examined area exhibited landslide predisposition mostly at high and very high levels.

Figure 6 summarized the performance of the SE, AHP, and the AHP–SE method for landslide susceptibility mapping based on ROC curves with AUC values and R-index for landslide probability levels. Overall, all three models were reliable, with the AUC higher than the random guess (0.5) and classified as a very good prediction (AUC ranging from 0.8–0.9). Interestingly, the integrated AHP–SE method outperformed both the two individual methods, with the AUC reaching the highest score of well over 0.87. Objective-based SE was better than subjective-based AHP, i.e., AUC values are nearly 0.85 and 0.81, respectively. Compared to the lowest AUC derived from the AHP, the synergy of the objective and subjective weighting method, AHP–SE, led to a notable improvement in the landslide susceptibility model by approximately 7%. The R-index derived for very high landslide sensitivity levels in AHP, SE, and AHP–SE maps is ̉63.35%, 71.09%, and 78.41%, respectively, indicating the best prediction of AHP–SE models.

Out of the 15 selected factors, the weighting approach subjectively and objectively resulted in different relative importance levels. For knowledge-based weighting using AHP, confirmation for the most important factor, slope, was also found in different case studies (Kayastha et al., 2013; Mondal and Maiti, 2013; Moragues et al., 2021; Agrawal and Dixit, 2022; Khalil et al., 2022). Slope is a widely used variable in the field of landslide hazard zonation (Pourghasemi et al., 2018; Reichenbach et al., 2018). Notwithstanding, its importance in relation to landslide sensitivity may vary depending on the specific circumstance. Recent research settings on different characteristics of the area found other factors to be more influential, such as land use (Guo et al., 2023) and lithology (Yalcin, 2008; Pourghasemi et al., 2013). Particularly in the tropical climate region, precipitation was recognized as the most significant component related to landslide sensitivity (Shahabi and Hashim, 2015). Related to the data-driven methods, the results indicated that calculation based on SE gave the most important factors of soil. It was also confirmed by the same method (Devkota et al., 2013; Agrawal and Dixit, 2022). For the integrated model, climate was the most influential factor in the landslide susceptibility model, especially in tropical monsoon areas. The factor was rarely used due to data availability and it is unnecessary in some general case studies. Here, we aimed to highlight the impact of the typical climate form under tropical monsoon areas as one of the key variables leading to the high degree of landslide triggers. An updated global map of landslide sensitivity strongly suggested that Southeast Asia, with its classical tropical environment, has one of the most frequent landslide occurrences (Stanley et al., 2021).

We also have a deep insight into the contribution of different landslide conditioning factors, which allows us to suggest the necessary datasets in regions with the same environmental conditions. Based on the AHP method, our subjective knowledge indicated that the top five important factors are soil, climate, land cover, TWI, and elevation. The AHP-based results suggested that each fundamental group contains at least one factor with a significant contribution to the landslide hazard predictive model. However, approaching data-based statistics, the SE method has shown that precipitation, slope, distance to fault, climate, and soil were estimated as the highest weight factors. For the combined AHP–SE, climate, soil, slope, elevation, and geomorphology become the highest weight factors of the landslide-prone predictive model. Considering both individual and integrated methods, two groups of land cover conditions and anthropogenic activities and hydrology contributed the least to the landslide susceptibility model. This is a confirmation for the most essential group of weather, topography, and geology compared to the two remaining groups of hydrology and natural–artificial conditions. The topography group also contained the most significant factor, slope, and the least important factor, plan curvature, in all the three methods.

It should be noted that subclasses in conditioning factors significantly affect the results of landslide susceptibility modeling, especially for high-weight factors. Figure 7 represented the percentage of landslide inventory validation recorded in each specific subclass, which analyzed the relationship between landslide occurrence and different subclasses in causal factors. In summary, the sensitivity of landslide over the case study area was observed mostly in hilly areas (98–268 m), steep slopes (18–27^o), slope aspects of east and northeast, both concave and convex of planform curvature, distance to river under 100 m, drainage density below 0.58 km²/km, TWI ranging from 0.64 to 6.26, distance to fault higher than 500 m, soil class of Epi Lithi Ferralic Acrisols, geomorphology landform of slow gravity slope, land cover of tree cover, distance to road higher than 500 m, high NDVI values, daily rainfall spanning between 11.1 and 13.3 mm/day, and the climate form of tropical monsoon climate Ia1b separated by dry seasons (high temperature and moderator drought) and rainy seasons (from May to November with the total amount of precipitation 2,500–2,800 mm).

Among the followed methods, AHP is one of the most common techniques due to its easy approach and implementation (Tyagi et al., 2022). Even with the subjective limitation, AHP was still conducted in many applications of case studies (Pourghasemi et al., 2018; Yong et al., 2022). The method applied in our study was less accurate than the SE, with no significant difference, however. Previous experiences also confirmed the slightly lower performance than that of data-driven methods such as statistical analysis (Panchal and Shrivastava, 2021) and machine learning models (Khalil et al., 2022). Nonetheless questions about the requirement of expert knowledge contributing to the landslide susceptibility model still remained, with the continuous proposal for different combined approaches. Various studies demonstrated a better performance by integration of AHP with different strategies, i.e., WLC (Hung et al., 2016), FR (Mondal and Maiti, 2013), fuzzy (Agrawal and Dixit, 2022), and evidential belief function (EBF) (Althuwaynee et al., 2014). Furthermore, incalculable fluctuation of changing climate leads to ill-posed problems, while depending individually on data-driven methods (Gariano and Guzzetti, 2016; Zêzere et al., 2017), thus emphasizing the necessity of knowledge-based contribution in diverse landslide susceptibility models.

On the other hand, the objective-based weighting approach, SE, predicted better landslide sensitivity than the AHP with the confirmation of all the mentioned accuracy metrics. The method was applied in numerous works as demonstrated for the best model while being compared to LR (Devkota et al., 2013), Bayesian conditional probability model (Pourghasemi et al., 2012), AHP (Panchal and Shrivastava, 2021), statistical information value (Singh et al., 2021), and FR (Jaafari et al., 2014). These techniques employ mathematical models to determine criteria weights without incorporating the subjective preferences of decision-makers, thereby avoiding the potential influence of personal bias on the final decision outcome. Despite the fact that weights assigned by the relationship of history landslide records and thematic layers are explicit and objective, the method is affected by the resolution of spatial explanatory factors and both the quality and quantity of landslide inventory data (Zêzere et al., 2017).

In the light of various combined methods, we deployed the synergy of AHP and SE, as noticed for the lack of a combination in the field of landslide susceptibility assessment. Generally, both of the proposal weighting approaches have self-advantages and disadvantages. The SE is based on objective data, while the field dataset of landslide occurrences and the detail or resolution of the thematic layers affect the performance of the predictive process. The AHP is based on subjective decisions that may have the tendency of following rigid perspectives or epistemic uncertainty. To overcome the disadvantages and yield advantages in both the methods, we integrated subjective and objective approaches to deliver an optimal weight for each corresponding landslide-controlling factor. Statistical metrics of the AHP–SE method were revealed in this research. The performance of the combined approach was demonstrated to be more accurate compared to that of AHP and SE, thus increasing the reliability of the predictive model.

In the present study, we involved 15 factors in an effort to consider the entire prevalent landslide-predisposing aspects of topography, geology, hydrology, land cover environment and anthropogenic activities, and weather. However, it is hard to confirm the best selection of conditioning factors as the task is a major challenge in hazard prediction modeling. Among hundreds of factors (with their original name) used in the landslide susceptibility, more than 23 factors were mostly used (Reichenbach et al., 2018). Analysis of studies during 2005–2016 also indicated the same number of factors used more than 30 times (Pourghasemi et al., 2018). Due to the unavailability of standard guidelines to select the most effective factors, the optimal factors are mainly based on landslide type, examined area conditions, methods applied in available data, and scale requirements (Pourghasemi et al., 2018). The accuracy of models by adopting different factors was demonstrated by Gaidzik and Ramírez-Herrera (2021). Therefore, future implementation should involve quantitative analysis to suggest the most significant factors and to improve our profound understanding of how predisposing factors affected landslide susceptibility.

The geospatial data DEM is the main source used to generate different landslide-related spatial products such as slope, aspect, and curvature. Although various studies confirmed the unimpacted topographic product resolution on the performance of landslide susceptibility models (Tian et al., 2008; Chen et al., 2020), suggestions for the most congruous raster pixel size of the DEM were also discussed in detail. Gaidzik and Ramírez-Herrera (2021) indicated that a finer resolution of topographic data leads to better prediction; however, the quantity and quality of input data seem to be important with lower resolution. Lee et al. (2004) concluded that the proper resolution of 30 m, while constructing the map with scales ranging from 1:5,000 to 1:50,000 (Tian et al., 2008), implied that for a changed study area, the size decided the meaningful spatial resolution, while also mentioning that flat, ridge, and slope foot terrain shape is more difficult to predict than landslide probability in the case of lower DEM resolution. In addition, the choice of suitable methods is also mentioned, for instance, the FR applied in low-resolution products seemed to be better than the entropy-based method and WoE (Chen et al., 2020), or the AHP is appropriate in the medium-resolution scale (1:250,000–1:25,000) (Yong et al., 2022). Therefore, further analysis needs to be considered for calculating the spatial resolution of topographic products in specific conditions of the area and for selecting the most suitable input data, methodologies, and map scale results.

Apart from the proposed methods already applied in this study, further experiments should attempt to continue accessing different methods including statistics, expert-based methods, artificial intelligence solutions, and a synergy of these methods to choose the most suitable approach (Chakrabortty et al., 2022). Big data also leads to the uncertainty of traditional methods due to the complex analysis process. Combining advanced high-performance machine learning and deep learning models with big geospatial data may help in resolving complicated issues and significantly improving the performance of the landslide model. This could be a potential solution for the transferability of accurate landslide susceptibility models over different areas and for temporal analysis.