Earthquake-Induced Landslide Susceptibility Assessment Using a Novel Model Based on Gradient Boosting Machine Learning and Class Balancing Methods

: Predicting the susceptibility of a speciﬁc part of a landslide (SSPL) involves predicting the likelihood that the part of the landslide (e.g., the entire landslide, the source area, or the scarp) will form in a given area. When predicting SSPL, the landslide samples are far less than the non-landslide samples. This class imbalance makes it difﬁcult to predict the SSPL. This paper proposes an advanced artiﬁcial intelligence (AI) model based on the dice-cross entropy (DCE) loss function and XGBoost (XGB DCE ) or Light Gradient Boosting Machine (LGB DCE ) to ameliorate the class imbalance in the SSPL prediction. We select the earthquake-induced landslides from the 2018 Hokkaido earthquake as a case study to evaluate our proposed method. First, six different datasets with 24 landslide inﬂuencing factors and 10,422 samples of a speciﬁc part of the landslides are established using remote sensing and geographic information system technologies. Then, based on each of the six datasets, four landslide susceptibility algorithms (XGB, LGB, random-forest (RF) and linear discriminant analysis (LDA)) and four class balancing methods (non-balance (NB), equal-quantity sampling (EQS), inverse landslide-frequency weighting (ILW), and DCE loss) are applied to predict the SSPL. The results show that the non-balanced method underestimates landslide susceptibility, and the ILW or EQS methods overestimate the landslide susceptibility, while the DCE loss method produces more balanced results. The prediction performance of the XGB DCE (average area under the receiver operating characteristic curve (0.970) surpasses that of RF (0.956), LGB (0.962), and LDA (0.921). Our proposed methods produce more unbiased and precise results than the existing models, and have a great potential to produce accurate general (e.g., predicting the entire landslide) and detailed (e.g., combining the prediction of the landslide source area with the landslide run-out modeling) landslide susceptibility assessments, which can be further applied to landslide hazard and risk assessments.


Introduction
Regional landslide susceptibility is defined as an estimation of the probability that a landslide will occur in a given area under certain environmental conditions [1,2].A landslide susceptibility assessment provides valuable information that can be utilized in a landslide hazard and risk assessment and/or in land-use planning and traffic route decision-making in order to reduce the impacts of future landslide damage.Regional landslide susceptibility models usually focus on either predicting the susceptibility of an entire landslide [3][4][5][6] or predicting the susceptibility of one specific part of a landslide (SSPL) [7][8][9][10][11][12][13][14].When the model focuses on predicting the SSPL, only the zones that contain that specific part (e.g., the source area or a landslide scarp) are defined as the landslide class.Compared with the prediction of the entire landslide, the prediction of SSPL can make the landslide susceptibility assessments more precise and widely applied.For example, the prediction of the landslide source area can be combined with a physical-based landslide run-out model to calculate the landslide deposition and the kinetic parameters of landslides (i.e., velocity and landslide impact frequency or probability) in a way that is more precise and consistent with physical laws.In addition, these parameters can be further used for quantitative landslide hazard and risk assessment.Therefore, it is of great significance to study the prediction of SSPL.
However, by dividing a landslide into smaller areas and only classifying one or a few of these areas as the landslide class, an extreme class imbalance arises between the number of non-landslide samples and the number of landslide samples.This class imbalance makes it more difficult to predict the SSPL than it is to predict the susceptibility of an entire landslide.In order to ameliorate this class imbalance, it is necessary to apply a class balancing method.There are three types of class balancing methods in current landslide susceptibility assessments: under-sampling techniques, over-sampling techniques, and cost-sensitive techniques.Under-sampling methods synthetically or randomly collect the same number of landslide samples and non-landslide samples using algorithms such as the equal-quantity sampling (EQS) [12][13][14], fuzzy c-means [15], easy ensemble [16], density peak clustering [17] and self-organizing-map [18] algorithms.While under-sampling significantly reduces the sample size and improves the computational efficiency, model training may be affected by the loss of data from the original samples.In over-sampling methods, algorithms such as SMOTE (synthetic minority oversampling technique) [19] and ADASYN (adaptive synthetic sampling) [20] sample or synthetically generate landslide samples to ensure that the number of landslide samples is equal to the number of non-landslide samples.While over-sampling maintains the sample balance and increases the overall sample size, it may result in higher computational costs due to the introduction of incorrect and/or redundant information into the sample data.In cost-sensitive techniques such as the direct approach, the pre-processing approach, and the post-processing approach, the class balance is maintained by explicitly accounting for the cost of model mispredictions during the training process [21].In general, most landslide susceptibility studies employ undersampling methods such as EQS, while over-sampling, cost-sensitive methods, inverse landslide-frequency weighting, and loss correction are hardly ever used.The EQS method sets a landslide frequency (in this article the term frequency refers to the relative frequency) of 50% in the model training and validation sets.Because the actual landslide frequency is much lower than 50%, the EQS method is not a true representation of real-world conditions, and thus yields a biased prediction of landslide susceptibility.
The susceptibility algorithm is important in predicting the SSPL, or the susceptibility of an entire landslide.Early research on this topic focused mainly on susceptibility algorithms reliant on data-driven statistical methods such as linear discriminant analysis (LDA) [2,22], logistic regression [8,23], and the weight of evidence method [24,25].However, the landslide influencing factors used in these algorithms must be conditionally independent and linearly independent, and it is difficult to solve complicated nonlinear problems using such algorithms.In the past ten years, machine learning techniques such as decision tree (DT) [26,27], random forest (RF) [5,9,28], support vector machine (SVM) [29][30][31][32], artificial neural network (ANN) [33][34][35], and deep neural network [36] have been widely used in landslide susceptibility assessments.While these machine learning techniques can deal with nonlinear problems, they struggle with large datasets, which limits the sample size and resolution of the resulting landslide susceptibility assessment.Recently, several new breakthrough ensemble learning methods, such as the XGBOOST (XGB) [37] and LightGBM (LGB) [38] methods, have been invented.XGB and LGB can handle large datasets and have more flexible and accurate objective functions.These ensemble machine learning methods can provide new solutions for landslide susceptibility prediction.
Although a great deal of research has been carried out on prediction of the SSPL, several problems remain unsolved.
(1) The EQS class balancing method is commonly used in current research.However, this method sets a landslide frequency of 0.5 in the model training and validation sets, which is not a true representation of real-world conditions.At present, a more reasonable class balancing method has not been proposed.(2) No studies in the literature have systematically compared and evaluated the application of different class balancing methods in SSPL prediction.(3) The commonly used linear and nonlinear landslide susceptibility algorithms have many limitations.Despite the numerous advantages of using the XGB and LGB algorithms, no studies have employed these frameworks in prediction of the SSPL.
The goals of this study were to develop a landslide susceptibility model that combines the XGB, LGB, and dice cross-entropy (DCE) loss function to solve the class imbalance problem and achieve better model performance in SSPL prediction.We explored how the choice of class balancing methods (EQS, inverse landslide-frequency weighting (ILW), DCE loss, or non-balance (NB) methods) and susceptibility algorithms (XGB, LGB, RF, or LDA) affects the prediction of the SSPL.In this study, the landslides triggered by the 2018 Hokkaido earthquake (M W 6.6) were taken as a case study, and the prediction of the susceptibility of six different landslide parts was considered.

Materials and Methods
The main framework of this article is presented in Figure 1.In Part I, we describe our study area and construct our dataset.Our high-quality dataset consists of 10,422 earthquakeinduced landslides in the 2018 Hokkaido earthquake and 24 landslide influencing factors.We discuss our methodology in Parts II and III.Our technique automatically extracts six different parts of the landslide (e.g., the entire landslide, the landslide source area, and the landslide scarp) and compiles separate inventories for them.Each type of landslide inventory is combined with 24 landslide factors and then separated into six datasets.We created 13 landslide susceptibility models using combinations of four susceptibility algorithms (XGB, LGB, linear discriminant analysis (LDA), and random forest (RF)) and four class balancing methods (non-balance (NB), equal-quantity sampling (EQS), inverse landslide-frequency weighting (ILW), and DCE loss).We input all six datasets into each of the 13 models and generated 78 sets of results.In each case, the modified group crossvalidation method (MGCV) was used for the model training, prediction, and validation.Part IV includes the Results and Discussion sections.We used machine learning metrics and landslide statistics to evaluate and analyze the results of the landslide susceptibility.

Study Area
Our study area (Figure 2) is located in the eastern Iburi area of central Hokkaido (bounded by latitudes of 42.49  2a).The terrain in this area is low and hilly, with flatter profiles in the west and steeper topography in the east.The altitude ranges from 1 m to 873 m and the slope degree ranges from 0.1 • to 57.7 • .The minimum and maximum monthly average temperatures are −7 • C and 20.4 • C, respectively.Most (42%) of the yearly precipitation falls between July and September, with a monthly average precipitation of 142 mm.
The lithology of our study area (Figure 2a) is predominantly (63.4%) marine sedimentary rocks such as marine siliceous mudstones, marine conglomerates, marine mudstones, and marine sandstones.The igneous rocks (5.2%) in our study area mainly consist of rhyolite, andesite, basalt, and diorite, while metasandstone chlorite, pelitic schist, and mafic schist are the main metamorphic rocks (1.1%) found in this area.The remaining strata are mainly Quaternary deposits.In terms of its structural features, our study area is located in an active compression collision zone created by the tectonic action of the Amur Plate, the Philippine Sea Plate, the Okhotsk Plate, and the Pacific Plate.Per the USGS seismic database, a M W 6.6 earthquake (epicenter location of 42

Study Area
Our study area (Figure 2) is located in the eastern Iburi area of central Hokkaido (bounded by latitudes of 42.49°N and 43.02°N and longitudes of 141.69°E and 142.36°E; Figure 2a).The terrain in this area is low and hilly, with flatter profiles in the west and steeper topography in the east.The altitude ranges from 1 m to 873 m and the slope degree ranges from 0.1° to 57.7°.The minimum and maximum monthly average temperatures are −7 °C and 20.4 °C , respectively.Most (42%) of the yearly precipitation falls between July and September, with a monthly average precipitation of 142 mm.
The lithology of our study area (Figure 2a) is predominantly (63.4%) marine sedimentary rocks such as marine siliceous mudstones, marine conglomerates, marine mudstones, and marine sandstones.The igneous rocks (5.2%) in our study area mainly consist of rhyolite, andesite, basalt, and diorite, while metasandstone chlorite, pelitic schist, and mafic schist are the main metamorphic rocks (1.1%) found in this area.The remaining strata are mainly Quaternary deposits.In terms of its structural features, our study area is located in an active compression collision zone created by the tectonic action of the Amur Plate, the Philippine Sea Plate, the Okhotsk Plate, and the Pacific Plate.Per the USGS seismic database, a MW 6.6 earthquake (epicenter location of 42°41′10″N 141°55′44″E and focal depth of 35 km) struck central Hokkaido on 6 September 2018.According to our landslide inventory, this earthquake triggered at least 10,422 landslides.

Basic Data
Before carrying out a landslide susceptibility assessment, it is necessary to prepare the basic data such as the landslide inventory and the landslide influencing factors.These data were gained and processed through the remote sensing technology, geographic information system (GIS) technology, and field surveys.The dataset for preparing the landslide inventory and the landslide influencing factors are presented in Table 1.By comprehensively considering the data quality, sample size, and calculation cost, we selected the grid cells (referred to as pixels) as the mapping units and resampled all of the basic data to a resolution of 30 m.
We prepared a landslide inventory containing 10,422 co-seismic landslides from the Hokkaido earthquake on 6 September 2018.Our landslide inventory was created based on the landslide data (3307 landslides) obtained from the Geospatial Information Authority of Japan and the landslide inventory of [39] (5625 landslides, covering an area of 46.3 km 2 ).We fully expanded and revised the original landslide inventories through visual interpretation of high-resolution Google Earth images and UAV photos taken before and after the earthquake (Table 1) and through digitization.We interpreted the landslides by starting at the area with the highest seismic intensity and working our way outwards until we found no further landslides; the total area of the interpretation was greater than 2000 km 2 .The interpretation and digitization of each landslide covered the range from the source area to the deposition area (Figure 2d).
The landslides in this inventory have a total area of 49.97 km 2 , a total volume of 62 × 10 6 m 3 (based on the volume formula for soil landslides proposed by Larsen, et al. [40]), and a density of 5.2 landslides/km 2 .The quartiles and the average of the area of the individual landslides are 1378 m 2 , 2763 m 2 , 5575 m 2 , 170,255 m 2 , and 4795 m 2 .The landslides mainly occurred in the Quaternary deposits and the mudstone layer, and can be classified as flow-slides, debris avalanches, debris slides, and rock slides.The study area had experienced prolonged rainfall and strong earthquakes, and numerous landslides liquefied as a result.
In our landslide susceptibility assessment, we focused on 24 landslide influencing factors (Table 2 and Figure 3), including the seismology, geomorphology, hydrology, geology, vegetation, and rainfall of our study area.PGA, PGVA, and distance focus are directly related to both the amount of ground shaking and earthquake intensity.The distance ridge factor provides insight into the seismic site effects.Sites that are located close to mountain ridges generally experience seismic amplification, as after propagating through a mountain, seismic waves are reflected, diffracted, and channeled toward the mountain ridge, where the waves interfere constructively with one another [41].By summarizing the relationship between the direction of seismic wave propagation and the slope attitude [42], factors such as the epicentral direction, angle ES , angle FH , and angle FS influence the occurrence of landslides.
In terms of topography and geomorphology, slope MAX and the elevation difference are used to describe the steepness of the terrain.Landslides rarely develop in gentle terrain because it is difficult to generate effective free faces and the potential energy required to trigger a landslide.The elevation ratio describes the position of a point with respect to the slope (e.g., at the bottom of the slope, on the shoulder of the slope, etc.).Differences in the slope stress state and the thickness of the weathering layer or unloading layer affect the slope stability as well.In terms of the hydrological conditions, we considered river factors such as distance SR , distance MR , and SPI [44].We accounted for these variables because a river can alter the stress state (e.g., by undercutting unloading or via the seepage force) and/or the physical properties (e.g., by softening the rock or soil) of a slope.According to the influence a river of a given size exerts on a slope, that river is classified as either a major river or a minor river according to the Strahler stream order.The meteorological factors include the cumulative precipitation in the four days preceding the earthquake.Because the earthquake occurred during the rainy season, we inferred that abundant rainfall facilitated the onset of multiple landslides during the earthquake by saturating and softening the slope rocks and soils.
Other geological factors include the lithological and structural features of a given area.The lithology affects the material composition, strength, structure, and weathering resistance of the slope.For example, mudstone is characterized by low strength, low structural integrity, and weak intercalation.Structural factors include distance fault , distance fold , and fault density.Slopes with faults and folds often develop additional joints, fissures, weathered zones, and groundwater paths; naturally, these features increase the slope instability.The angle between epicentral direction and the slope aspect (angle ES ) DEM and QuiQuake 0.0-360.0 The angle between the horizontal and the line from calculated cell to focus (angle FH ).Angle FH represents the direction of seismic wave propagation at a location, which would influence the occurrence of landslides. 42.9-90.0 The sum of angle FH and the slope degree (angle FS ) 43.7-126.8 • 30 m

Geomorphologic
The maximum slope in the neighbourhood of the calculated cell (slope MAX ) DEM 0.0-57.7 •

m
The variation of the slope aspects in the neighbourhood of the calculated cell (aspect VAR ) 0.0-1.0-30 m The ratio of the elevation to the maximum elevation in the neighbourhood of the calculated cell (elevation ratio) Elevation difference 0.0-233.0m 30 m The percentage of convex cells in the neighbourhood of the calculated cell (surface convexity) [43].Surface convexity describes the shape (convex, concave, flat) of the slope, which affects the stability of the slope under earthquake shaking.In order to establish the landslide susceptibility model, we developed an algorithm that automatically and efficiently extracts the different parts of the landslide.This algorithm assumes that each part of a landslide has a different relative elevation [45].As is shown in Figure 4a, a landslide can be divided into two parts at a certain elevation E m , namely, the upper part (i.e., landslide scarp or landslide source area) as the landslide class and the lower part (i.e., landslide accumulation area) as the non-landslide class; this allows us to specifically extract the areas where the landslide initiates, such as the source area.To facilitate this process, we can define the relative elevation ratio R d (0-100%) as follows: where E max and E min (m) are the highest and lowest elevation of the landslide, respectively, and H d and H (m) represent the height difference of the accumulation zone and the height of the entire landslide, respectively.In order to accurately represent the landslide feature, we chose to define the landslide parts as different sections of a surface rather than as a single point (e.g., the centroid).The grid cells are used as the mapping units, with the landslide class and non-landslide class cells assigned values of one and zero, respectively.
Remote Sens. 2022, XX, x FOR PEER REVIEW 10 of 30 Using our algorithm (Figure 4b), we were able to automatically extract and identify 62,532 parts representing six different landslide features (Figure 4c-f).Instead of extracting the landslide accumulation area separately, we extracted the area containing the landslide source area for landslide susceptibility prediction, as the landslide accumulation area is not the location where the landslide initiates.We compared and analyzed the extraction results with satellite images of the landslides.An Rd value of 0% corresponds to an entire Using our algorithm (Figure 4b), we were able to automatically extract and identify 62,532 parts representing six different landslide features (Figure 4c-f).Instead of extracting the landslide accumulation area separately, we extracted the area containing the landslide source area for landslide susceptibility prediction, as the landslide accumulation area is not the location where the landslide initiates.We compared and analyzed the extraction results with satellite images of the landslides.An R d value of 0% corresponds to an entire landslide, from the source area to the accumulation area.An R d of 10-30% encompasses the parts of the landslide that are not the accumulation area, such as the source and the transport area of a rock avalanche.An R d value of 50-70% corresponds to the source area.An R d value of 90% corresponds to the landslide scarp.In the automatic landslide extraction algorithm, the accuracy of the landslide polygon boundaries is related to the output resolution, while the accuracy of landslide cell extraction and identification is affected by the DEM quality.A statistical analysis of each part of the landslide is shown in Table 3.Here, we use realistic landslide frequency values, rather than the ideal value of 0.5 that has often been used in previous studies.Landslide frequency: the ratio of the sum of the areas of the landslide parts to the total study area.

Validation Method
In this study, we devised and implemented the modified group cross-validation method (MGCV) for validation of the landslide susceptibility models.In this method, we treat each landslide or landslide part as a whole; in this way, we ensure that the cells in a given landslide can only be assigned to the training set or to the validation set.The process of the MGCV is as follows: (1) MGCV randomly divides all of the landslides (or landslide parts) in the inventory into K G groups (K G = 4 in this study), where each independent group is composed of the same number of samples.It then randomly selects non-landslide cells and distributes them evenly amongst the four groups.This ensures that the training set, validation set, and total sample set have the same realistic landslide frequency.(2) Then, taking the i-th group of samples as the model validation set, the remaining K G − 1 groups are combined into a model training set, the model landslide susceptibility prediction is compared to the observed landslides of the validation set, and the evaluation metric is calculated.(3) Finally, step 2 is repeated K G times to complete the cross-validation process, with each group of samples used as the validation set exactly once.The model performance is the average of the metric scores calculated in K G iterations.

Evaluation Metrics
The receiver operating characteristic curve (ROC) [46] is an evaluation metric frequently used in landslide susceptibility assessments.The area under the ROC (AUC) is used as a quantitative indicator of model performance; a higher AUC corresponds to a better model performance.This study employs the ROC as a baseline reference for comparison with other studies, not as a core indicator of model performance.In fact, the ROC only measures the relative ranks of landslide susceptibility, rather than the size of the landslide susceptibility value.
The Jaccard index (Jaccard), which ranges from 0 to 1, measures the similarity between the landslide susceptibility and the presence or absence of a landslide.Larger Jaccard values correspond to better model performance.In this paper, Jaccard is the geometric mean of Jaccard P and Jaccard N , which measure the similarity in negative samples (i.e., the non-landslide cells) and positive samples (i.e., the landslide cells), respectively: By equally considering the contributions of both the positive and negative samples to the final prediction result, the Jaccard value avoids placing too much influence on classes that have proportionally more samples in a given dataset.More traditional metrics (i.e., accuracy score) can be biased by such classes.The F1 score is ideal for tasks in which a class imbalance may arise, as it considers both the recall and precision: where Precision = TP/(TP + FP) and Recall = TP/(TP + FN).The true positive (TP), false positive (FP), true negative (TN), and false negative (FN) are calculated using the confusion matrix approach [46], where the terms positive and negative refer to the presence or absence of landslides, respectively.

Landslide Susceptibility Algorithms
XGBoost (XGB) [37] is a gradient boosting framework that supports GPU and parallel computing.Because XGB is flexible, efficient, and accurate, it can handle large datasets such as those found in landslide susceptibility assessments.XGB is a tree-based ensemble machine learning approach.The generation of an XGB model begins with a single classification and regression tree (CART); in each model iteration, a new CART is created using the current error, then that CART is then added to the output of the existing model.The XGB model grows iteratively in this fashion until the convergence criteria are satisfied.The output of an XGB model is the combination of the outputs of K CARTs: where i represents the i-th sample, x i is the landslide influencing factors, ŷi is the landslide susceptibility prediction, K is the number of iterations, and f t (x i ) is the landslide susceptibility prediction generated by CART T t added in the t-th iteration.Here, T t is determined mainly by minimizing L (t) , which is the objective function where ŷ(t) is the output of the model in the t-th iteration; y i represents whether landslides occur (1) or do not occur (0), l is the loss function, and R(f t ) is the regularization term of the model.By further expanding n ∑ i=1 l(y i , ŷ(t) i ) using the second-order Taylor expansion, we account for the structure of the CART in the objective function and rewrite the objective function as [37]: Where , V is the total number of leaves in T t , I j is a sample group belonging to leaf node j, and α and β are regularization coefficients.The optimal structure of T t and f t (x i ) can be determined based on Equation (7) and the specific growth mode of the CARTs.This ŷ(t) i calculation completes the tth iteration.
LightGBM (LGB) [38] is a gradient boosting framework that is similar to XGB.LGB was specifically designed to expedite the slower parts of XGB.Linear discriminant analysis (LDA) and random forest (RF) have been widely used in studies of debris flows and landslides.Detailed descriptions of LDA and RF can be found in [47] and the [48], respectively.None of the aforementioned algorithms considers hyperparameter optimization, because LDA, RF, XGB, and LGB all have good initial hyperparameters.The main hyperparameters are described as follows.In LDA, the solver is singular value decomposition and the threshold for a singular value is 1 × 10 −4 .In RF, the number of trees is 100, the split criterion is Gini impurity, and the depth of the tree is unlimited.In XGB, the number of trees is 100, the learning rate is 0.3, the subsample ratio is 1, and the depth of the tree is 6.In LGB, the number of trees is 100, the learning rate is 0.1, the subsample ratio is 1, and the depth of the tree is unlimited.The above algorithms and the DCE loss function (described below) were coded and calculated using Python's scientific computing libraries (i.e., NumPy, XGBoost and LightGBM).

DCE Loss Function
This study proposed two new XGB and LGB frameworks that employ the DCE loss function.We synthesize the DCE loss using the dice loss [49] and weighted cross-entropy loss algorithms, where the dice loss is actually equivalent to the F1 score.Let l(y, ŷ) be defined as an n-dimensional vector l(y, ŷ) = [l 1 , l 2 , . . . ,l n ] where l i is the DCE loss of the i-th sample.The relevant formulas are as follows: where l C i and l D i are the weighted cross-entropy loss and dice loss of the i-th sample, respectively, and a and (1 − a) are the weights of the non-landslide and landslide samples in l C i respectively.It is possible to address the class imbalance by using a smaller a weight in the training process; this choice reduces the contribution of non-landslide samples and increases the contribution of landslide samples.The parameters b and (1 − b) are the weights of l C i and l D i in l i .The parameters a and b must be tuned in order to identify the optimal value.The n l parameter is the coefficient that adjusts the relative sizes of l C i and l D i , while p i represents the landslide susceptibility prediction, where p i = 1/(1 + e − ŷi ).
The first-and second-order gradients of l(y, ŷ) with respect to ŷ are the Jacobi and Hessian matrices, respectively.From the diagonal elements of these matrices, we can identify the first-order gradient and its ith term: where In addition, we can identify the second-order gradient and its i-th term: The results of Equations ( 12) and ( 14) are substituted into the g i and h i terms of Equation (7), respectively, for the t-th iteration of model training.

Class Balancing Method
We explored four class balancing methods: the equal quantity sampling (EQS), inverse landslide frequency weighting (ILW), loss correction, and non-balance (NB) methods.EQS creates a class-balanced dataset by combining the landslide samples and an equal number of non-landslide samples obtained through random under-sampling.The landslide frequency of the dataset produced using this method is 0.5.In order to ensure that the landslide frequency in the validation set is realistic (rather than 0.5), the EQS in this study was used to generate only the training sets, instead of both training and validation sets as in previous studies.ILW addresses the class imbalance by using the landslide prior probabilities of the samples to amplify the weight of the landslide samples in the model training process.The weights of the landslide and non-landslide samples are set to (N l + N f )/N l and (N l + N f )/N f , respectively, where N l and N f are the number of landslide samples and non-landslide samples in the training set, respectively.
When using the loss correction method, loss functions are created or corrected in order to improve the model training when datasets are class imbalanced.In this study, we employed the DCE loss correction (described in Section 2.3.6), which is specifically applicable to the XGB and LGB.For comparison, we considered a model that does not include a class balancing method, which we name the non-balance (NB) model.
These four methods for balancing classes were used in the training of the MGCV method.The four landslide susceptibility algorithms were combined with four class balancing methods to create 13 landslide susceptibility models.We refer to a given landslide susceptibility model using nomenclature based on the algorithm and class balancing method used for that model; for example, XGB NB represents the model combining XGBoost (the susceptibility algorithm) and NB (the class balancing method).

Model Performance Evaluation
As discussed previously, we used 13 different models to predict the landslide susceptibility of various landslide parts.As is shown in Table 4, for every model, increasing R d values correspond to decreasing PRE, REC, F1, JAC P , JAC, and Jaccard values.This means that when landslide parts are reduced in size, the prediction performance of all models deteriorates.The AUC only decreased by up to 5% (i.e., RF ILW ) with increasing R d ; this means that this metric is insensitive to changes in the landslide part.The evaluation metric scores for each model varied considerably with the class balancing method, especially at higher values of R d .We grouped the models according to their balancing methods and calculated the average values of the evaluation metric scores (considering different R d ).The EQS (or ILW) method yielded the highest recall (mean recall: 0.92) and the lowest precision (mean precision: 0.21).The NB method produced the highest precision (0.52) and the lowest recall (0.18).The precision (0.38) and recall (0.60) values for the DCE loss (XGB DCE and LGB DCE ) models were more balanced.
The models, again classified by their class balancing methods, were comprehensively evaluated based on the F1 and Jaccard metrics.The DCE loss method produced the best average F1 (0.465) and Jaccard (0.540) scores, while the NB method produced the lowest average F1 (0.258) and Jaccard (0.365) scores.The EQS and ILW methods produced F1 (0.311 and 0.293) and Jaccard (0.406 and 0.381) scores that fell in between these two extremes.We conducted statistical analyses on the model results based on their primary algorithms.The average F1 scores for the algorithms, from best to worst, are XGB (0.369), LGB (0.337), RF (0.304), and LDA (0.163) (Table 4).
As mentioned previously, a and b in the XGB DCE and LGB DCE models represent specific model parameters.Here, we investigate how these values affect the models themselves and their scoring metrics (Figure 5).As is shown in Figure 5a,e, the highest AUC and F1 scores occurred when a = 0.2.A value of index a that is too high or too low proves detrimental to the performance of the model.When b = 0, the DCE loss function reduces into the dice loss function, resulting in very low F1 and AUC values.When b > 0, the AUC metric remains largely unchanged (Figure 5f).When b = 1, the DCE loss function reduces into the weighted cross-entropy loss function, causing a sharp reduction in the F1 score.This suggests that the model does not perform as well as the DCE loss function does when either the Dice loss or cross-entropy loss is used separately as an individual loss function.The best AUC and F1 scores correspond to a b value of 0.05.The LGB DCE model produced F1 and AUC trends (Figure 5c,d,g,h) that are similar to those in XGB DCE .
As mentioned previously, a and b in the XGBDCE and LGBDCE models represent specific model parameters.Here, we investigate how these values affect the models themselves and their scoring metrics (Figure 5).As is shown in Figure 5a,e, the highest AUC and F1 scores occurred when a = 0.2.A value of index a that is too high or too low proves detrimental to the performance of the model.When b = 0, the DCE loss function reduces into the dice loss function, resulting in very low F1 and AUC values.When b > 0, the AUC metric remains largely unchanged (Figure 5f).When b = 1, the DCE loss function reduces into the weighted cross-entropy loss function, causing a sharp reduction in the F1 score.This suggests that the model does not perform as well as the DCE loss function does when either the Dice loss or cross-entropy loss is used separately as an individual loss function.The best AUC and F1 scores correspond to a b value of 0.05.The LGBDCE model produced F1 and AUC trends (Figure 5c,d,g,h) that are similar to those in XGBDCE.LGBDCE models.As Rd increases, the part of the landslide that the model seeks to predict is reduced from the entire landslide to the landslide scarp.

Landslide Susceptibility Mapping
Using the MGCV method, we used four iterations of each model to generate cross predictions of landslide susceptibility for our entire study area (Figures 6 and 7).Because the results produced by the LGB and XGB algorithms are relatively similar, we chose to highlight the XGB and RF algorithm results in our figures.The statistical summary of the cells in each model (Figures 6 and 7) are summarized in Table 5.
The models using either the EQS or ILW class balancing methods included the largest number of cells in the high susceptibility class.However, the NB method models had the largest number of cells with a susceptibility of 0.0-0.1 (>92% of the total), indicating that the landslide susceptibility determined using the NB method is low overall.As the Rd value increase from 0.0 to 0.9, the high-susceptibility area (0.6-1.0) predicted by the XGBDCE model changed from the entire landslide area (Figure 6g) to the landslide source area (Figure 6h) and the landslide scarp (Figure 6i), demonstrating that the XGBDCE model is able to precisely adjust the spatial distribution of the landslide susceptibility according to the target landslide part.With the EQS balancing method, when Rd increases from 0.0 to 0.9, the number of cells in the high-susceptibility area (0.6-1.0)only decreases by 10% (Table 5).Figures 6d-f and 7d-f show that the high-susceptibility areas on the hillside are

Landslide Susceptibility Mapping
Using the MGCV method, we used four iterations of each model to generate cross predictions of landslide susceptibility for our entire study area (Figures 6 and 7).Because the results produced by the LGB and XGB algorithms are relatively similar, we chose to highlight the XGB and RF algorithm results in our figures.The statistical summary of the cells in each model (Figures 6 and 7) are summarized in Table 5.The models using either the EQS or ILW class balancing methods included the largest number of cells in the high susceptibility class.However, the NB method models had the largest number of cells with a susceptibility of 0.0-0.1 (>92% of the total), indicating that the landslide susceptibility determined using the NB method is low overall.As the R d value increase from 0.0 to 0.9, the high-susceptibility area (0.6-1.0) predicted by the XGB DCE model changed from the entire landslide area (Figure 6g) to the landslide source area (Figure 6h) and the landslide scarp (Figure 6i), demonstrating that the XGB DCE model is able to precisely adjust the spatial distribution of the landslide susceptibility according to the target landslide part.With the EQS balancing method, when R d increases from 0.0 to 0.9, the number of cells in the high-susceptibility area (0.6-1.0)only decreases by 10% (Table 5).Figures 6d-f and 7d-f show that the high-susceptibility areas on the hillside are much larger than the actual landslide areas; therefore, the EQS method cannot precisely predict a specific part of the landslides.For the NB method, Figures 6a-c and 7a-c show that the XGB NB and RF NB models are able to accurately adjust the spatial distribution of the landslide susceptibility according to the target landslide part.However, in these models, the susceptibility values within the dense landslide area were low.
We quantified the relationship between the proportion of landslides and the susceptibility on the landslide susceptibility maps (Figure 8).In the NB-based model, the proportion of landslides is negatively correlated with landslide susceptibility.However, in the EQS-based and DCE-based models the proportion of landslides increases with the landslide susceptibility, indicating that the landslide proportion in these models is more reasonable than that in the NB model.The proportion of landslides varies with increasing R d .In the NB-based model, the proportion of landslides in the 0.0-0.2susceptibility class increased very rapidly (average value of 55%).In the XGB DCE models, only when R d reaches 0.9 does the proportion of landslides decrease as the susceptibility value increases.In the EQS-based models, the proportion of landslides changed relatively little.The above results indicate that the proportion of landslides in EQS-based and DCE-based models is more stable than that in NB-based models when the target landslide part changes.Our findings on the relationship between the landslide frequency and the susceptibility on the landslide susceptibility maps are summarized in Figure 9.The landslide frequency differs significantly amongst the various models.The landslide frequency in the NB-based models is closest to the reference line (the actual landslide frequency), while that in the EQS-based models is very low and furthest away from the reference line.The landslide frequency of the DCE-based models falls in between the above two extremes.
As is shown in Figure 9, the landslide frequency decreased as R d increased in all of the models.In the EQS model, the largest landslide frequency reduction occurred in the 0.6-1.0susceptibility class (on average 6.5 times lower).In this susceptibility class, the model predicted that many landslide cells were actually non-landslide cells (i.e., false positives).For the NB-based model, when R d ≥ 0.7, the landslide frequency in the highsusceptibility class is even lower than that in the low susceptibility class.The frequency of landslides in the XGB DCE model decreased by a factor of only 1.5 on average.Therefore, the model was maximally stable when there were changes to the target landslide parts.0.6-1.0susceptibility class (on average 6.5 times lower).In this susceptibility class, the model predicted that many landslide cells were actually non-landslide cells (i.e., false positives).For the NB-based model, when Rd ≥ 0.7, the landslide frequency in the high-susceptibility class is even lower than that in the low susceptibility class.The frequency of landslides in the XGBDCE model decreased by a factor of only 1.5 on average.Therefore, the model was maximally stable when there were changes to the target landslide parts.

The Analysis of Landslide Influencing Factors
Next, we explored how various landslide influencing factors affect the prediction of different landslide parts in these models.We quantified the influence exerted by each factor using the reduction of the loss function in the XGB-based model (Figure 10).The rank, Mean, and STD of the landslide influencing factors vary greatly depending on the model (Figure 10).Because the XGBDCE model generated the most accurate predictions, we focus on the results for this model here.According to the ranked mean of the landslide influencing factors (Figure 10c), the PGA, curvatureSTD, angleFS, angleFH, precipitation,

The Analysis of Landslide Influencing Factors
Next, we explored how various landslide influencing factors affect the prediction of different landslide parts in these models.We quantified the influence exerted by each factor using the reduction of the loss function in the XGB-based model (Figure 10).The rank, Mean, and STD of the landslide influencing factors vary greatly depending on the model (Figure 10).Because the XGB DCE model generated the most accurate predictions, we focus on the results for this model here.According to the ranked mean of the landslide influencing factors (Figure 10c), the PGA, curvature STD , angle FS , angle FH , precipitation, distance focus , and lithology are the primary factors controlling landslide occurrence.It is worth noting that earthquakes and rainfall events are factors that can trigger landslides.As shown in Figure 11a, there is a significant positive correlation between the PGA and landslide frequency.However, contradicting known patterns, more precipitation corresponds to a lower landslide frequency, that is, the maximum landslide frequency (0.07) was observed in the low rainfall (5-17 mm) area (Figure 11f).Clearly, the landslides in this study area were more heavily influenced by ground motion than by the amount of rainfall.
Figure 10c shows that certain factors, such as PGA, curvature STD , angle FS , angle FH , precipitation, and distance focus , have relatively low-ranking standard deviations.Additionally, because there is no comparative change in the landslide frequency for a given factor (Figure 11) for increasing values of R d , we infer that these factors are not sensitive to changes in the parts of the landslide.Conversely, the ranked standard deviations of factors such as slope MAX , distance SR , distance ridge , elevation ratio, and stream power index (SPI) are significantly larger (Figure 10c).The relationships between these factors and landslide frequencies vary for the different parts of the landslide.For higher values of R d , the landslide frequency increases with the elevation ratio and decreases with the SPI and distance ridge metrics (Figure 11i,k,l).Therefore, the source areas of the landslides are most likely to be found on the upper slopes of the hills, closer to the ridges.The opposite trend is observed for smaller values of R d (Figure 11).
Remote Sens. 2022, XX, x FOR PEER REVIEW 22 of 30 distancefocus, and lithology are the primary factors controlling landslide occurrence.It is worth noting that earthquakes and rainfall events are factors that can trigger landslides.As shown in Figure 11a, there is a significant positive correlation between the PGA and landslide frequency.However, contradicting known patterns, more precipitation corresponds to a lower landslide frequency, that is, the maximum landslide frequency (0.07) was observed in the low rainfall (5-17 mm) area (Figure 11f).Clearly, the landslides in this study area were more heavily influenced by ground motion than by the amount of rainfall.Figure 10c shows that certain factors, such as PGA, curvatureSTD, angleFS, angleFH, precipitation, and distancefocus, have relatively low-ranking standard deviations.Additionally, because there is no comparative change in the landslide frequency for a given factor (Figure 11) for increasing values of Rd, we infer that these factors are not sensitive to changes in the parts of the landslide.Conversely, the ranked standard deviations of factors such as slopeMAX, distanceSR, distanceridge, elevation ratio, and stream power index (SPI) are significantly larger (Figure 10c).The relationships between these factors and landslide frequencies vary for the different parts of the landslide.For higher values of Rd, the landslide frequency increases with the elevation ratio and decreases with the SPI and distanceridge metrics (Figure 11i,k,l).Therefore, the source areas of the landslides are most likely to be found on the upper slopes of the hills, closer to the ridges.The opposite trend is observed for smaller values of Rd (Figure 11).
The mechanism of the earthquake-induced landslides in this study area was analyzed based on the main landslide influencing factors, as follows.The type of the soil and rock determines the strength of the slope.The landslide frequency of siliceous mudstone and marine conglomerate is significantly higher than that of other lithologies (Figure 11h), mainly because these two lithologies have low strength and weak regoliths.The slip surfaces of most landslides in the study area are shallow and nearly planar; therefore, the slip surface should mainly follow the interface of the bedrock and the overburden or interface between different layers in the overburden (mostly composed of the residual soil, unweathered volcanic ash, weathered pyroclastic soil, and pumice [50]).These weak interfaces are low in yield acceleration, and are more likely to deform under the seismic load.Precipitation plays a key role in landslide occurrence (Figure 10c).Prolonged rainfall before the earthquake (the cumulative precipitation in August and September was 200-300 mm) increased the water content of the soil and rock, softening them.The seismic factors are the most important landslide influencing factors (Figure 10c).Under strong ground motion, the pore pressure of the soil increased rapidly (most of the soils in the study area are fine grained with impeded drainage), leading to a significant decrease in the undrained shear strength (and the yield acceleration) of the soil [51], with soil layers even being liquefied.Each time the ground acceleration exceeded the yield acceleration, the slope mass slid along the aforementioned weak interfaces.In the process, the joints and cracks of the slope propagated, and the deformation modulus and the strength of the slope decreased [52,53].When the permanent displacement of the slope reaches a certain threshold (with higher PGA values associated with larger permanent displacement [54]), landslides are very likely to occur.Most landslides in the study area initiated and slid along the valley channel for a relatively long distance.The slip surfaces on the slope are smooth, without residual sliding mass [39].These phenomena indicate that the water content of the sliding mass was relatively high.The mechanism of the earthquake-induced landslides in this study area was analyzed based on the main landslide influencing factors, as follows.The type of the soil and rock determines the strength of the slope.The landslide frequency of siliceous mudstone and marine conglomerate is significantly higher than that of other lithologies (Figure 11h), mainly because these two lithologies have low strength and weak regoliths.The slip surfaces of most landslides in the study area are shallow and nearly planar; therefore, the slip surface should mainly follow the interface of the bedrock and the overburden or interface between different layers in the overburden (mostly composed of the residual soil, unweathered volcanic ash, weathered pyroclastic soil, and pumice [50]).These weak interfaces are low in yield acceleration, and are more likely to deform under the seismic load.Precipitation plays a key role in landslide occurrence (Figure 10c).Prolonged rainfall before the earthquake (the cumulative precipitation in August and September was 200-300 mm) increased the water content of the soil and rock, softening them.
The seismic factors are the most important landslide influencing factors (Figure 10c).Under strong ground motion, the pore pressure of the soil increased rapidly (most of the soils in the study area are fine grained with impeded drainage), leading to a significant decrease in the undrained shear strength (and the yield acceleration) of the soil [51], with soil layers even being liquefied.Each time the ground acceleration exceeded the yield acceleration, the slope mass slid along the aforementioned weak interfaces.In the process, the joints and cracks of the slope propagated, and the deformation modulus and the strength of the slope decreased [52,53].When the permanent displacement of the slope reaches a certain threshold (with higher PGA values associated with larger permanent displacement [54]), landslides are very to occur.Most landslides in the study area initiated and slid along the valley channel for a relatively long distance.The slip surfaces on the slope are smooth, without residual sliding mass [39].These phenomena indicate that the water content of the sliding mass was relatively high.

The Applicability of the Automatic Extraction Method
We devised an automatic extraction method to efficiently identify and extract the various landslide parts (Figure 4).This extraction method only requires DEM data and the polygonal landslide inventory.Using our technique, an ordinary computer can extract 360,000 landslide parts in 6 s.This automatic method is applicable to both the new landslide inventories and older landslide inventories that are missing data or are use low quality data.In contrast, manual extraction methods require experienced specialists, high-quality imagery, and DEMs.For large study areas, it may take months to prepare the data using manual extraction methods.While automatic landslide extraction methods that are based on relative elevation ratios (R d ) are empirical, their precision is sufficient for regional landslide studies, in which the resolution is usually >10 m.Within this resolution, landslides are represented by a small number of cells.Even when more precisely determining the boundaries of each part of the landslides using manual methods and then extracting these parts, the extraction results for most of the landslides change only slightly.Due to the inherent subjectivity and additional constraints present in the manual extraction method, automatic extraction methods should always match the accuracy of manual extraction methods at the regional scale.

Comparison and Prospect of Landslide Learning Algorithms
We mainly measured landslide susceptibility model performance using F1 score instead of the commonly used AUC of the ROC, as the number of landslide and nonlandslide samples in our dataset was extremely imbalanced.In such cases, the F1 score is far more stringent than the AUC of the ROC.Therefore, in this study the F1 score is far less than the AUC of the ROC (Table 4), which is greater than 0.9 (indicates an excellent model performance) [55] for all of the models.
In order to enhance the standards of landslide susceptibility evaluation, the landslide susceptibility algorithm should meet multiple criteria.Compared with more traditional methods, the XGB and LGB algorithms all enable the user to adjust the loss functions based on the specific application.In this study, we devise the XGB DCE and LGB DCE models, which are based on DCE loss algorithm, to address the extreme class imbalance that is inherent in landslide susceptibility modeling.On average, the F1 values of these models are 59% (XGB DCE ) and 49% (LGB DCE ) than that of the RF EQS (the optimum RF model).In terms of the computational efficiency, we experimented with RF and ANN as well.The computational costs of the RF and ANN algorithms are 66 and 114 times higher, respectively, than that of the XGB and LGB models.An increasing demand for large datasets is anticipated in the future.For example, it is necessary to analyze landslide susceptibility using rapidly updated data (i.e., meteorological data) or high-resolution data at larger regional scales.Because the XGB and LGB algorithms are more accurate in prediction and more efficient in computation, they have wide applicability in future landslide susceptibility assessment.

Influence and Applicability of Class Balancing Methods
We explored how landslide susceptibility prediction is affected by different class balancing methods.While the NB-based models, which did not utilize any class balancing method, yielded the highest precision, they had the lowest recall ratio (Table 4).Most of the landslides predicted by these NB models (approximately 76%) are located in the lowest susceptibility area (0.0-0.2) (Figure 8).As such, we conclude that these NB models underestimate landslide susceptibility.This underestimation is due to most of the samples being in the non-landslide class; however, the weights of the landslide and non-landslide class are the same in the objective function of NB method, and therefore predicting the samples as the non-landslide class (i.e., with a low landslide susceptibility) is beneficial to the optimization of the objective function of the NB method.Therefore, the NB model is suitable for small-scale earthquake events with low landslide frequencies, in which the underestimation of the landslide susceptibility does not result in serious consequences.Moreover, the more precise values found in these NB models can accurately predict small potential landslide areas.However, when landslides are triggered by large earthquakes (e.g., the Wenchuan earthquake [56] and Nepal earthquake [57]), the NB-based models may fail to predict most landslides.
The EQS and ILW class balancing methods, on the other hand, have the highest recall values and the lowest precision values (Table 4).The landslide frequency of the high landslide susceptibility area (0.8-1.0) in these models is very low.Therefore, we conclude that these methods result in significant overestimation of the landslide susceptibility.This overestimation occurs because these techniques involve an inherent amplification of the landslide frequency or the landslide-class weight in the training dataset.For example, when using the EQS-based model, the landslide frequency in the training dataset was set to 0.5, while the landslide frequency of historically strong earthquakes such as the 2018 Hokkaido earthquake (M w 6.6), the 2008 Wenchuan earthquake (M w 7.9) [56], and the 2015 Nepal earthquake (M w 7.8) [57] were only 0.047, 0.12, and 0.01, respectively.Therefore, we conclude that it is not appropriate to use the EQS and ILW class balancing methods in models used to predict landslide susceptibility in areas affected by small or medium scale earthquake events.Due to overestimation, regardless of which part of the landslide is used to create the model, these models can only predict the overall landslide susceptibility of entire hillsides, and are unable to predict the susceptibility of smaller landslide features (Figures 6 and 7).
The prediction performance of the DCE loss-based models is far better than that of the other class balancing methods (Table 4).In addition, by using a specific weighted F1 score as the optimization objective, the DCE loss method achieves more balanced precision and recall, thereby solving the problem of overestimating or underestimating the landslide susceptibility that occurs in other class balancing methods.In order to implement the DCE loss function, the user must determine the values of parameters a and b.Although the maximum number of non-landslide cells is 136 times that of landslide cells (Figure 5), the optimal weighted ratio of positive to negative is only four.Therefore, the a value does not need to refer to the prior landslide probability (i.e., 1/136 when R d = 0.9) of an event.Under different R d values, the pattern of the F1 score and AUC curve is consistent (Figure 5).These trends demonstrate that reasonable values of a and b can be applied to a dataset of different parts of a landslide as well as a dataset of different frequencies (e.g., earthquake events of different intensities).
Model performance decreased significantly when shifting the focus from predicting the susceptibility of an entire landslide to predicting the susceptibility of a specific landslide feature (Table 4).The main reason for these results is that predicting low-frequency events (a specific landslide feature) is more difficult than predicting high-frequency events (an entire landslide).Despite this, studies have found [13,23] landslide susceptibility models that perform better when focused solely on landslide scarps.These contradictory results can be explained by the fact that the frequencies of the entire landslide and the landslide scarp in these studies [13,58] were both set to 0.5, a value that is dozens of times larger than the actual landslide frequency.Therefore, for that specific study, it became easier to predict the susceptibility of landslide scarps and the prediction performance was significantly overestimated.
While standard models can only be applied to predict the susceptibility of entire landslides, the XGB DCE and LGB DCE algorithms can be applied to predict the susceptibility of smaller landslide features.The susceptibilities predicted for an entire landslide by the XGB DCE and LGB DCE models (Figure 6g) can be applied across large geographic regions.With regard to predicting a specific part of a landslide, such as the landslide source area (Figure 6h,i), these prediction results can be combined with a physical-based landslide runout model to calculate the landslide deposition area and the kinetic parameters of landslides in a way that is more precise and consistent with physical laws than the traditional methods (Figure 12).Because these algorithms are both flexible and accurate, they can be applied to large-scale regions to increase the precision of landslide hazard or risk assessments.As shown in Figure 12a,b, we calculated the landslide impact frequency and impact probability parameters for landslide hazard and risk assessment, respectively.The spatial distribution of these parameters is in agreement with the observed landslides.

Limitations
First, the predicted susceptibility of landslide scarps (R d = 0.9) was not ideal.An attempt to improve the model performance which involved the addition of 30 landslide influencing factors related to the microtopographical features and earthquakes was unsuccessful.There is a relationship between subpar model performance and the lack of high-quality data.At a resolution of 30 m, it is difficult to gather high-quality data that exactly describe the characteristics of a landslide scarp.However, data with resolution higher than 30m were unavailable in this study area.Earthquakes and rainfall have been proven to be the two most significant factors driving the triggering of landslides (Figure 10c).However, the available seismic data and rainfall data do not account for the localized characteristics of the rainfall or earthquakes (i.e., the site effects).
Second, certain factors that can improve the performance of landslide susceptibility models were not considered in the analysis of the earthquake-induced landslide mechanism in Section 3.3.(1) The undrained strength of the overburden and of the interface between the bedrock and overburden.Most of the landslides in the study area are shallow translational landslides, and therefore these strength factors are helpful in predicting this type of landslide.(2) The parameters that reflects the physical properties of the soil and rock under the rainfall conditions (e.g., degree of saturation, relative density, coefficient of permeability).The occurrence of earthquake-induced landslides in the study area is related to the prolonged rainfall before the earthquake.These factors could be considered in the model to better measure the effect of rainfall.(3) The profile characteristics of the soil and rock, such as the type, depth, and sequence.These factors influence the characteristics of the slip surface of landslides and the site effect of ground motions.However, at present, data for the above factors in the study area are unavailable.Due to the large spatial variability of these geotechnical factors [32], it is necessary to carry out a large number of field investigations and experiments in the future.
Remote Sens. 2022, XX, x FOR PEER REVIEW 26 of 30 traditional methods (Figure 12).Because these algorithms are both flexible and accurate, they can be applied to large-scale regions to increase the precision of landslide hazard or risk assessments.As shown in Figure 12a,b, we calculated the landslide impact frequency and impact probability parameters for landslide hazard and risk assessment, respectively.The spatial distribution of these parameters is in agreement with the observed landslides.

Limitations
First, the predicted susceptibility of landslide scarps (Rd = 0.9) was not ideal.An attempt to improve the model performance which involved the addition of 30 landslide influencing factors related to the microtopographical features and earthquakes was unsuccessful.There is a relationship between subpar model performance and the lack of high-quality data.At a resolution of 30 m, it is difficult to gather high-quality data that exactly describe the characteristics of a landslide scarp.However, data with resolution higher than 30m were unavailable in this study area.Earthquakes and rainfall have been proven to be the two most significant factors driving the triggering of landslides (Figure 10c).However, the available seismic data and rainfall data do not account for the localized characteristics of the rainfall or earthquakes (i.e., the site effects).
Second, certain factors that can improve the performance of landslide susceptibility models were not considered in the analysis of the earthquake-induced landslide mechanism in Section 3.3.(1) The undrained strength of the overburden and of the interface between the bedrock and overburden.Most of the landslides in the study area are shallow translational landslides, and therefore these strength factors are helpful in predicting this type of landslide.(2) The parameters that reflects the physical properties of the soil and rock under the rainfall conditions (e.g., degree of saturation, relative density, coefficient of permeability).The occurrence of earthquake-induced landslides in the study area is related to the prolonged rainfall before the earthquake.These factors could be considered in the model to better measure the effect of rainfall.(3) The profile characteristics of the Third, our method is not able to predict landslide occurrence prior to extreme triggering events such as strong earthquakes or storms.The spatial distribution of landslides largely depends on the intensity of the triggering event (i.e., the intensity or PGA of the earthquake).However, our method only establishes the relationship between the intensity of a single event and the occurrence of landslides.As such, this method is not yet suitable for landslide prediction related to future triggering events of different sizes.In future work, we intend to use more landslide inventories to thoroughly investigate the relationship between seismic event intensity and landslides in order to improve the generalization ability of our models.

Conclusions
This paper proposes a landslide susceptibility model which combines the XGB (or LGB) and the DCE loss function for predicting the susceptibility of a specific part of a landslide (SSPL).We systematically explored how different susceptibility algorithms (XGB, LGB, and RF) and class balancing methods (EQS, ILW, NB, and DCE) affect the prediction of the SSPL.For our study, we prepared six different datasets with 24 landslide influencing factors and 10,422 samples of a specific part of the landslides from the Hokkaido earthquake on 6 September 2018, and established samples of realistic landslide frequencies.Our results demonstrate that XGB and LGB outperform RF and LDA.
The class balancing method heavily influences the predictive performance of different landslide susceptibility models.The low recall of NB methods results in a tendency to underestimate the landslide susceptibility of a given area.NB models are suitable for low-intensity events (i.e., heavy rainfall and small earthquakes) that can trigger landslides.The EQS and ILW methods result in low precision values and significant overestimation of landslide susceptibility.Those two methods should only be applied in large-scale landslide assessments with low resolutions.Predicting a smaller part of a landslide is more difficult than predicting an entire landslide, and a model with a reasonable class balancing method is required.The proposed XGB DCE and LGB DCE models produce more balanced precision and recall values and predict the SSPL more accurately.When combined with the use of an automatic extraction algorithm on different parts of the landslide, these models are applicable for assessing landslide susceptibility, hazard, and risk at various scales and levels of landslide frequency.
In future research, higher quality data should be collected and applied to further improve the predictive performance of the models by focusing on smaller landslide features.A second possible improvement is to include more earthquake-induced landslide inventories and landslide influencing factors in future iterations.Ideally, our model should be able to evaluate landslide susceptibility either prior to a seismic event or in near real-time after the earthquake has begun.

Figure 1 .
Figure 1.The scheme of the present study.

Figure 1 .
Figure 1.The scheme of the present study.

Figure 2 .
Figure 2. Overview map of the study area.(a) Geologic map at a scale of 1:200,000 overlain by the distribution of co-seismic landslides (source: geologic map in Table 1.The classification of rock groups has been modified.).(b) The locations and satellite imagery (source: http://goto.arcgisonline.com/maps/World_Imageryaccessed on 1 September 2022) of the study area.(c) Distribution of peak ground acceleration in the Hokkaido region during the 6 September 2018, earthquake (source: QuiQuake in Table 1).(d) Satellite imagery (source: Maxar imagery acquired on 11 September 2018) of the co-seismic landslides in the local zone in (a).

Figure 2 .
Figure 2. Overview map of the study area.(a) Geologic map at a scale of 1:200,000 overlain by the distribution of co-seismic landslides (source: geologic map in Table 1.The classification of rock groups has been modified.).(b) The locations and satellite imagery (source: http://goto.arcgisonline.com/maps/World_Imagery accessed on 1 September 2022) of the study area.(c) Distribution of peak ground acceleration in the Hokkaido region during the 6 September 2018, earthquake (source: QuiQuake in Table 1).(d) Satellite imagery (source: Maxar imagery acquired on 11 September 2018) of the co-seismic landslides in the local zone in (a).

Remote
Sens. 2022, XX, x FOR PEER REVIEW 7 of 30 factors such as the epicentral direction, angleES, angleFH, and angleFS influence the occurrence of landslides.

Figure 3 .
Figure 3.The landslide influencing factors used in this study, obtained and processed by remote sensing and GIS technology.The yellow pentagram represents the epicenter.The names of the factors in (a-w) are below the color bar.

Figure 3 .
Figure 3.The landslide influencing factors used in this study, obtained and processed by remote sensing and GIS technology.The yellow pentagram represents the epicenter.The names of the factors in (a-w) are below the color bar.

Figure 4 .
Figure 4. Automatic extraction method for different parts of a landslide.(a) The division and extraction of the landslide parts.(b) A flowchart of our automatic extraction methodology.(c-f) The results of the landslide parts corresponding to Rd values of (c) 0.1, (d) 0.3, (e) 0.5, and (f) 0.9.The solid black lines represent the boundaries of the landslides.The colored cells represent the extracted parts of the landslide.

Figure 4 .
Figure 4. Automatic extraction method for different parts of a landslide.(a) The division and extraction of the landslide parts.(b) A flowchart of our automatic extraction methodology.(c-f) The results of the landslide parts corresponding to R d values of (c) 0.1, (d) 0.3, (e) 0.5, and (f) 0.9.The solid black lines represent the boundaries of the landslides.The colored cells represent the extracted parts of the landslide.
Gradient-based one-sided sampling reduces the training cost while preserving the accuracy of the model.The histogram algorithm divides the continuous landslide influencing factors into discrete bins and uses those bins to create new CARTs, which greatly reduces the costs of computation and storage in model training.Exclusive feature bundling binds mutually exclusive landslide factors to reduce the factors involved in model training, which significantly improves training speed without affecting model accuracy.

Figure 5 .
Figure 5.The influence that parameters a and b in the DCE loss function exert on landslide susceptibility prediction using the (a,b,e,f) XGBDCE and (c,d,g,h) LGBDCE models.As Rd increases, the part of the landslide that the model seeks to predict is reduced from the entire landslide to the landslide scarp.

Figure 5 .
Figure 5.The influence that parameters a and b in the DCE loss function exert on landslide susceptibility prediction using the (a,b,e,f) XGB DCE and (c,d,g,h) LGB DCE models.As R d increases, the part of the landslide that the model seeks to predict is reduced from the entire landslide to the landslide scarp.

Figure 6 .
Figure 6.Prediction of the susceptibilities of various landslide parts with combinations of the XGB algorithm with different class balancing methods: (a-c) XGBNB, (d-f) XGBEQS, (g-i) XGBDCE.As Rd increases, the part of the landslide that the model seeks to predict (the black polygon) is reduced from the entire landslide to the landslide scarp.

Figure 6 .
Figure 6.Prediction of the susceptibilities of various landslide parts with combinations of the XGB algorithm with different class balancing methods: (a-c) XGB NB , (d-f) XGB EQS , (g-i) XGB DCE .As R d increases, the part of the landslide that the model seeks to predict (the black polygon) is reduced from the entire landslide to the landslide scarp.

Figure 7 .
Figure 7. Prediction of the susceptibilities of various landslide parts with combinations of the RF algorithm with different class balancing methods: (a-c) RFNB, (d-f) RFEQS, (g-i) RFILW.As Rd increases, the part of the landslide that the model seeks to predict (the black polygon) is reduced from the entire landslide to the landslide scarp.

Figure 7 .
Figure 7. Prediction of the susceptibilities of various landslide parts with combinations of the RF algorithm with different class balancing methods: (a-c) RF NB , (d-f) RF EQS , (g-i) RF ILW .As R d increases, the part of the landslide that the model seeks to predict (the black polygon) is reduced from the entire landslide to the landslide scarp.

Figure 8 .
Figure 8.The relationship between the proportion of landslides and the landslide susceptibility.From subfigures (a) to (f), the Rd increases, and the part of the landslide that the model seeks to

Figure 8 .
Figure 8.The relationship between the proportion of landslides and the landslide susceptibility.From subfigures (a-f), the R d increases, and the part of the landslide that the model seeks to predict is reduced from the entire landslide to the landslide scarp.Because some of the plots almost overlap (XGB-based model versus LGB-based model and EQS-based model versus ILW-based model), for greater clarity we only show representative model plots.

Figure 9 .
Figure 9.The relationship between the landslide frequency and the landslide susceptibility.The models with the plots closest to the reference lines have predicted susceptibility values that are close to the actual landslide frequency.(a-f), the Rd increases, and the part of the landslide that the model seeks to predict is reduced from the entire landslide to the landslide scarp.Because some of the plots almost overlap (XGB-based model versus LGB-based model and EQS-based model versus ILWbased model), for clarity we show only representative model plots.

Figure 9 .
Figure 9.The relationship between the landslide frequency and the landslide susceptibility.The models with the plots closest to the reference lines have predicted susceptibility values that are close to the actual landslide frequency.(a-f), the R d increases, and the part of the landslide that the model seeks to predict is reduced from the entire landslide to the landslide scarp.Because some of the plots almost overlap (XGB-based model versus LGB-based model and EQS-based model versus ILW-based model), for clarity we show only representative model plots.

Figure 10 .
Figure 10.The ranking of the importance of different landslide influencing factors using the (a) XGBNB, (b) XGBEQS, (c) XGBDCE, and (d) XGBILW models.Different Rd (0.0-0.9) values correspond to predictions of different landslide parts in the models when landslide influencing factors are ranked by importance.Mean and STD refer to the mean and standard deviation, respectively, of the ranking of each landslide factor's importance for all Rd values.Lower rankings represent more significant factor contributions. EL = elevation, diff = difference, and DIR = direction.

Figure 10 .
Figure 10.The ranking of the importance of different landslide influencing factors using the (a) XGB NB , (b) XGB EQS , (c) XGB DCE , and (d) XGB ILW models.Different R d (0.0-0.9) values correspond to predictions of different landslide parts in the models when landslide influencing factors are ranked by importance.Mean and STD refer to the mean and standard deviation, respectively, of the ranking of each landslide factor's importance for all R d values.Lower rankings represent more significant factor contributions. EL = elevation, diff = difference, and DIR = direction.

Figure 12 .
Figure 12.A combination of landslide source area prediction and landslide run-out modelling portraying the respective distribution of (a) the landslide impact frequency (expressed as the cumulative distribution function) and (b) the landslide impact probability.The landslide run-out model uses r.randomwalk [59], in which the motion is controlled by the landslide angle of reach.The average angle of reach (25°) and the probability distribution function used for modeling were calculated during the statistical analysis of the landslides.

Figure 12 .
Figure 12.A combination of landslide source area prediction and landslide run-out modelling portraying the respective distribution of (a) the landslide impact frequency (expressed as the cumulative distribution function) and (b) the landslide impact probability.The landslide run-out model uses r.randomwalk [59], in which the motion is controlled by the landslide angle of reach.The average angle of reach (25 • ) and the probability distribution function used for modeling were calculated during the statistical analysis of the landslides.
• N and 43.02 • N and longitudes of 141.69 • E and 142.36 • E; Figure • 41 10"N 141 • 55 44"E and focal depth of 35 km) struck central Hokkaido on 6 September 2018.According to our landslide inventory, this earthquake triggered at least 10,422 landslides.

Table 1 .
Data set sources and meta-data.
Geospatial Authority Institute of Japan (https://fgd.gsi.go.jp/download/menu.phpaccessed on 1 September 2022) QuiQuake = Quick estimation system for earthquake maps triggered by observation records; CHIRPS = Climate Hazards Group InfraRed Precipitation with Station data.

Table 1 .
Data set sources and meta-data.
).(d) Satellite imagery (source: Maxar imagery acquired on 11 September 2018) of the co-seismic landslides in the local zone in (a).

Table 2 .
Description of factors influencing landslides.
F × W F /A F ). L F , W F , and A F are the total fault length, fault width, and area of the statistical zone, respectively.
O × W O /A O ).L O , W O ,and AO are the total fold length, fold width, and area of the statistical zone, respectively.

Table 3 .
Inventories of different landslide parts.

Table 4 .
Model prediction performance of different landslide features for every model.

Table 5 .
Statistical analysis of the number of cells (unit: thousand) in the generated landslide susceptibility maps for every model.