A Data-Driven Model on Google Earth Engine for Landslide Susceptibility Assessment in the Hengduan Mountains, the Qinghai–Tibetan Plateau

: Amplifying landslide hazards in the backdrop of warming climate and intensifying human activities calls for an integrated framework for accurately evaluating landslide susceptibility at ﬁne spatiotemporal resolutions, which is critical for the mitigation of increasingly high landslide disaster risks. Yet, dynamic landslide susceptibility mapping is still lacking. Using high-quality data, from 14,435 landslides and non-landslides, we developed an efﬁcient holistic framework for evaluating landslide susceptibility, considering landslide-relevant internal and external factors based on cloud computing platform and algorithmic models, which enables dynamic updating of a landslide susceptibility map at the regional scale, particularly in regions with highly complicated topographical features such as the Hengduan Mountains, as considered in this study. We compared Classiﬁcation and Regression Trees (CART), Support Vector Machines (SVM), and Random Forest (RF) classiﬁers to screen out the best portfolio model for landslide susceptibility mapping on the Google Earth Engine (GEE) platform. We found that the Random Forest (RF) classiﬁer integrated with synergy mode had the best modeling performance, with 90.48% and 89.24% accuracy and precision, respectively. We also found that forests and grasslands had the controlling effect on the occurrence of landslides, while human activities had a notable inducing effect on the occurrence of landslides within the Hengduan Mountains. This study highlights the performance of the holistic landslide susceptibility evaluation framework proposed in this study and provides a viable technique for landslide susceptibility evaluation in other regions of the globe. The validation results showed 0.9611 for the synergy model AUC, 0.9569 for the static model AUC, and 0.9419 for the dynamic model AUC. In addition, the accuracy indices calculated by the confusion matrix showed 90.48% in accuracy and 89.24% in precision for the synergy model, 89.74% in accuracy and 87.30% in precision for the static model, and 87.39% in accuracy and 83.80% in precision for the dynamic model. All these statistical metrics showed that the assessment model using synergy mode had higher evaluation accuracy and better performance. The results indicated that the paradigm of the leverage synergy mode combined with the machine learning algorithm based on the GEE platform for large-scale LS assessment was viable and feasible and had a satisfactory assessment performance. and Q.Z.; formal analysis, W.W.; investigation, W.W. and J.Z.; data curation, W.W.; writing—original preparation, W.W. and Q.Z.; writing—review and Q.Z., G.W., J.Z., Z.S. and S.S.; visualization, W.W.;


Introduction
A landslide is a natural disaster that may cause thousands of fatalities and inflict massive damage on infrastructure [1][2][3]. With the projected increase in extreme precipitation [4], land resources dwindling, and urban development spiraling, landslides are increasing in frequency, scope, and destructive capacity [5]. The International Disaster Database (EM-DAT) shows that between 2000 and 2020, landslides accounted for 4.51% of all natural disasters, while 13.05% of these landslide disasters occurred in China [6]. The actual landslide disasters and their risk may be higher than recorded in EM-DAT landslide catalogs, since many landslides go unreported [7,8]. Landslide susceptibility is regions with occurrence of landslides from a holistic viewpoint. In addition, this is needed for the mitigation of potential landslide hazards at the regional scale, especially in the Hengduan Mountains [60,61]. It should be noted that LS mapping is a challenging task, particularly for regions with complicated topographical features, since the driving factors triggering landslides are complicated, so an efficient assessment framework is needed. To develop a holistic LS evaluation framework based on GEE is a critical first step for the dynamic update of LS map and is also of practical significance for regional mitigation of landslide disasters in a warming climate. Therefore, development of an efficient LS evaluation framework combining static and dynamic explanatory factors, based on multisource datasets and cloud computing platform and machine learning algorithms, is urgently needed. This point constitutes the major motivation of this study.
Here, we present a synergy of static and dynamic explanatory variables based on GEE to evaluate LS at the regional scale and propose a holistic LS evaluation framework. The major objectives of this study, therefore, are to (1) evaluate the influence of the types of landslide-related factors (static, dynamic and both) on LS, (2) elucidate the controlling factors of LS across the Hengduan Mountains region, and (3) propose and develop an efficient and near-real-time LS evaluation framework at the regional scale. This study highlights an efficient LS evaluation framework and helps provide a technique framework for the mapping of LS in other regions of the globe.

Study Region
The Hengduan Mountains region is located in the southeastern part of the Qinghai-Tibetan Plateau and is a general term for the north-south-oriented parallel mountains in the western part of Sichuan and Yunnan provinces and the eastern part of the Tibet Autonomous Region, covering an area of more than 700,000 km 2 , with elevation ranging from 124 m to 7473 m ( Figure 1). Actually, the Hengduan Mountains region is a complete geomorphological partition, in which different subregions have similarities in geomorphological forms, tectonic features, and geomorphological development processes [62]. The Hengduan Mountains region has been recognized as being highly susceptible to landslides because of frequent earthquakes, extreme precipitation, and intensifying human activities [15,39,59,63].  The study area is located in the SSW of China, with an area of more than 700,000 square kilometers.

Landslide Inventories
A landslide inventory map contains valuable information about the spatial pattern of landslide events in any given region, enhancing the understanding of landslide behavior Remote Sens. 2022, 14, 4662 4 of 24 and the LS evaluation. It is a critical step in LS studies to make a landslide inventory map [23]. Here, we obtained historical landslide data from Natural Resources Bureaus of eight administrative divisions, including Tibetan Qiang Autonomous Prefecture of Ngawa, Diqing Tibetan Autonomous Prefecture, Tibetan Autonomous Prefecture of Garzê, Nujiang of the Lisu Autonomous Prefecture, Liangshan Yi Autonomous Prefecture, Lijiang City, Nyingchi, and Ya'an, covering 74 districts and counties accounting for 48.59% of the Hengduan Mountains region area, with a total of 21,112 geological disaster points in three categories, i.e., landslide, rockfall, and debris flow. The period of landslides covers a period of 2000-2020. Here, we use the term "landslide" to describe rotational or translational mass movement slide type. However, a portion of the raw landslide data lacks geographic coordinates, time information of the landslides, and description information. Therefore, we exclude this part of landslide data with missing information, for reliability of the assessment results. After cleaning and eliminating data absent of attributes, a total of 7217 landslide locations (centroid) were obtained (Figure 2), and we named this landslide points dataset as Expedition Data. In addition, we collected landslide data from the Resource and Environmental Science and Data Center of the Chinese Academy of Sciences (https://www.resdc.cn/ accessed on 10 July 2020), which covers the whole range of the study area and is mainly used for visual comparative verification of LS assessment results. We named this dataset as Investigation Data. The statistical time of expedition data and the investigation data were up to 2020. Due to different sources of landslide data, landslide surveys vary in detail and landslide inventories. Here, we use expedition data as training data for LS evaluation and investigation data as visual comparison data for evaluating performance.

Landslide Inventories
A landslide inventory map contains valuable information about the spatial pattern of landslide events in any given region, enhancing the understanding of landslide behavior and the LS evaluation. It is a critical step in LS studies to make a landslide inventory map [23]. Here, we obtained historical landslide data from Natural Resources Bureaus of eight administrative divisions, including Tibetan Qiang Autonomous Prefecture of Ngawa, Diqing Tibetan Autonomous Prefecture, Tibetan Autonomous Prefecture of Garzê, Nujiang of the Lisu Autonomous Prefecture, Liangshan Yi Autonomous Prefecture, Lijiang City, Nyingchi, and Ya'an, covering 74 districts and counties accounting for 48.59% of the Hengduan Mountains region area, with a total of 21,112 geological disaster points in three categories, i.e., landslide, rockfall, and debris flow. The period of landslides covers a period of 2000-2020. Here, we use the term "landslide" to describe rotational or translational mass movement slide type. However, a portion of the raw landslide data lacks geographic coordinates, time information of the landslides, and description information. Therefore, we exclude this part of landslide data with missing information, for reliability of the assessment results. After cleaning and eliminating data absent of attributes, a total of 7217 landslide locations (centroid) were obtained (Figure 2), and we named this landslide points dataset as Expedition Data. In addition, we collected landslide data from the Resource and Environmental Science and Data Center of the Chinese Academy of Sciences (https://www.resdc.cn/ accessed on 10 July 2020), which covers the whole range of the study area and is mainly used for visual comparative verification of LS assessment results. We named this dataset as Investigation Data. The statistical time of expedition data and the investigation data were up to 2020. Due to different sources of landslide data, landslide surveys vary in detail and landslide inventories. Here, we use expedition data as training data for LS evaluation and investigation data as visual comparison data for evaluating performance.

Explanatory Variables
Selection of appropriate explanatory variables is a critical step in the LS assessment, but no commonly accepted selection criterion is available [40]. Based on landslide behavior across the Hengduan Mountains region [34,51,64,65], 50 explanatory variables were accepted as input variables for LS assessment. For the LS evaluation across the large-scale and spatially heterogeneous region, we subdivided these explanatory variables into three categories, i.e., static, dynamic, and triggering factors [39,41,66]. The underlying surface-related explanatory variables that were related to sliding behaviors of landslides were taken as static factors [67,68]. Here, we defined dynamic factors as the explanatory variables derived from remote sensing images associated with dynamic changes in the underlying surface. Explanatory variables that trigger mass movements were called triggering factors [69], such as rainfall and earthquakes considered in this study. The explanatory variables are listed in Table A1.
Static factors involved 9 explanatory variables from DEM, road net, river net, and geology fault vector data. Elevation has been evidenced to be substantially viable in depicting spatial distribution of landslides [70]. Slope is an important factor for landslide processes, since it controls the shear forces acting on mountain hill slopes [19,71]. Aspect is associated with diurnal anisotropic heat, root development, weathering, and other geomorphological processes [72,73]. Curvature shows topographic relief that is considered to be a causative factor for landslide, which controls the convergence or divergence of landslide materials and water in the direction of landslide movement [23]. The terrain ruggedness index is an indicator of the variation of surface relief and the degree of erosion [25]. The topographic wetness index (TWI) reflects the effects of topography and soil characteristics on the spatial distribution of soil moisture [29]. The distance from fault indicates the likelihood of landslides induced by ruptures and fault surfaces easily become sliding surfaces because the stress on the rock surrounding a fault is unstable [74]. The distance from river and LS are in close interaction, and this distance indirectly describes the erosion power of streams, which plays an important role in landslides [24]. The distance from the road net is used as an indicator to measure landslides derived from or related to road construction [75].
Dynamic factors included 33 explanatory variables, which were obtained or calculated from Sentinel-1/Sentinel-2 and Landsat satellite data (Table A1). Satellite remote sensing data (SAR or optical data) have been widely used in landslide identification and spatiotemporal analysis of landslide activity [40,43,46,47,76]. In general, LS are identified by analyzing vegetation cover changes and relief-oriented parameters using different vegetation indices [77]. The texture of soil represents the relative proportion of sand, silt, and clay content, which is highly correlated with landslide susceptibility. In addition, landslide initiation shows close dependency to high soil moisture levels and low vegetation density [78][79][80][81]. Here, we use various vegetation indices to characterize the fingerprint characteristics of land cover changes and spectral-texture information relevant to landslides. The amplitude component of an SAR image is related to surface roughness, surface material, soil water, and radar polarization, which has a strong relation with landslide occurrence [76,82]. In order to better characterize the dynamic changes of surface material composition and soil water content, we used the amplitude data of SAR with different polarization methods and calculated the spatial correlation of information after data smoothing [83,84]. In this study, we integrated optical and SAR satellite data incorporated with extracted index information and Geary's C spatial analysis information to delineate the dynamic change process of the underlying surface (Table A1). Considering the diversity of dynamic factors and remarkable value ranges of the dynamic factors, we normalized the dynamic factors before further analysis.
Triggering factors are mainly rainfall and earthquakes in this study area, according to the hazard property information of landslide inventories. Quantile level analysis of precipitation data can reflect the precipitation characteristics of the study area [85]. Here, we used the distance from the epicenter and precipitation percentiles to determine the spatial clustering of landslides [8,86].

LS Assessment Framework
The LS assessment framework proposed in this study is based on binary classification over the study area through pixel-by-pixel calculation, which was separated into four phases: (a) preparation of training and validation sample data, (b) construction of data features, (c) comparison of the selected three LS assessment models, and (d) mapping of the LS across the study region ( Figure 3). used the distance from the epicenter and precipitation percentiles to determine the spatial clustering of landslides [8,86].

LS Assessment Framework
The LS assessment framework proposed in this study is based on binary classification over the study area through pixel-by-pixel calculation, which was separated into four phases: (a) preparation of training and validation sample data, (b) construction of data features, (c) comparison of the selected three LS assessment models, and (d) mapping of the LS across the study region ( Figure 3).  Table A1.  Table A1.

Sample Data
The LS assessment model requires both landslides, denoted as 1, and non-landslides, denoted as 0 [87]. The expedition landslide inventory data of 7217 landslide sites were taken as positive samples. Here, we used the Create Random Point module in ArcGIS to generate non-landslide samples for further modeling [16]. Two strategies were adopted to avoid a blind selection of data points introducing uncertainty into the model training procedure which may produce an unreliable LS assessment. At first, the spatial range of generated random points should be in consistency with these eight municipalities where expedition data were collected. Second, random points within a 1 km buffer zone from landslide positive samples were screened out. At last, 7218 non-landslide points were generated as negative samples. The ratio of landslide and non-landslide samples is determined as 1:1 [19]. Thus, a total of 14,435 samples were produced for LS assessment. Sample data were randomly split into dataset for model training (70% of the total dataset) and model verification (30% of the total dataset) [88].

Development of Feature Modes
We grouped explanatory variables into three feature modes for LS assessment, i.e., static mode, dynamic mode, and synergy mode. Static factors and triggering factors were combined to form the static mode; combination of dynamic factors and triggering factors formed the dynamic mode; and the integration of static factors, dynamic factors and triggering factors formed the synergy mode. Triggering factors are those factors that cause mass movements, so we added triggering factors for different data modes. A set of performance metrics for different feature mode testing was recorded during mode validation using the RF method. Then, the feature mode with the highest modeling performance was used to map LS. The preparation of the feature modes includes the calculation of static factors and dynamic factors. The static factors were calculated through ArcGIS and SAGA GIS software and then uploaded to GEE LEGACY ASSETS. The calculation of dynamic factors is implemented in GEE by coding and then stored in GEE LEGACY ASSETS.

LS Assessment Model
Three classifiers were chosen for LS assessment, i.e., Random Forest (RF), Classification and Regression Tree (CART), and Support Vector Machines (SVM). Our analysis based on these three classifiers was done via GEE, allowing large-scale predictions at the pixel scale [57].
RF classifier is an efficient and optimal in computation performance in dealing with high dimensionality of data [89], which has been widely used in satellite-based applications [28].
CART is a rule-based algorithm that splits the dataset subsets using all predictor variables to create two child nodes repeatedly, and the predicted value of a "terminal" node is the average of the response values in that node [90]. Predictor variables can be of any type (numeric, binary, categorical, etc.), and model outputs cannot be influenced by monotone transformations and different scales of measurement among predictors [91][92][93].
SVM, a kernel-based algorithm for classification and regression issues, has been widely used in LS mapping [19,94]. For nonlinear feature datasets, as in LS modeling, SVM is used to project features into high dimensional space with kernel functions, allowing classification within a plane [95].
We verified these three classifiers using the same feature mode. A set of performance metrics for testing were recorded for comparison of these classifiers. Then, a model selected with the highest modeling performance and the well-chosen feature mode were used to map the LS over the Hengduan Mountains region.

LS Mapping
The landslide identification model developed here was used to analyze the spatial pattern of landslides at the pixel scale with a spatial resolution of 30 m. LS is strictly a classification problem, with the binary outcome of presence or absence of a landslide [96]. The predicted value of 1 means that landslides will occur within the pixel range, and the predicted value of 0 means no landslides within the pixel range. The high spatial resolution of 30 m of the predicted result allows us to better identify the fingerprint of landslides [13]. Besides, whether or not landslide occurs at a given location is closely related to the triggering factors within and in the proximity of the regions with the occurrence Remote Sens. 2022, 14, 4662 8 of 24 of landslides. To construct the landslide susceptibility map, two main steps should be followed: first generate the landslide susceptibility indexes (LSIs), and then reclassify the LSIs. Here, we aggregated the predicted landslides and summed the pixel values of landslides based on the 10 × 10 domain range. So, we obtained a 300 m resolution spatially aggregated raster image, at which the size of the image element values indicates the occurrence possibility of landslides within a 300 m × 300 m range (Equation (1)). The abovementioned processing procedure can well avoid the subjective judgment error of the model-based assessment results and can better reflect the trend of landslide occurrences through statistical characteristics of landslides within a specific region. Finally, we used the Jenks Natural Breakpoint method to classify the LS into five levels, seeking to minimize the differences within levels and maximize the differences between levels [11,29]. Here, we classified the LSI into very low LS, low LS, medium LS, high LS, and very high LS: where m and n are the row and column number of LS assessment unit, respectively, with spatially aggregated 10 × 10 domain range. V ij is the predicted pixel value corresponding to row i and column j.

Validation of LS Assessment Framework
The expedition landslide data collected through field survey and the investigation landslide data mapped by Chinese Academy of Science (CAS) were used to validate the LS assessment framework proposed in this study. Expedition landslide data was divided into training data and validation data. The training data were used for the development of the model and the validation data were used for the selection of the model. Here, we used the receiver operating characteristic curve (ROC) and area under curve (AUC) to evaluate and compare the effects of different feature modes on the LS map [97,98]. The algorithm used in this study included binary classification, confusion matrix, and statistical indicators (i.e., precision, accuracy), which were adopted to quantitatively compare the accuracy of LS assessment models [99]. In binary classification, the accuracy is the proportion of correct predictions (both true positives and true negatives) among the total number of cases, while precision is the ratio of true positives to the total of the true positives and false positives. The formula is: Accuracy = (TP + TN)/(TP + TN + FP + FN) where TP = true positive; FP = false positive; TN = true negative; and FN = false negative. Investigation landslide data mapped by CAS delivers valuable information on the spatial pattern of landslides in the landslide susceptible zone. Then, we used the investigation landslide data to compare the modeling results from a macro perspective and to qualitatively measure the spatial consistency of the LS maps.

Development Environment of Landslides
The spatial pattern of triggering factors behind landslides can decide the spatial patterns of landslides. Therefore, we depicted the spatial pattern of triggering factors behind landslides. Static factors included nine explanatory variables ( Figure 4). The topography of the Hengduan Mountains region is generally high in the north and low in the south, with remarkable differences in terrain, which cause unstable mountain slopes. The mountains run generally in a north-south direction between 25 • and 30 • N, with a slight tilt toward the west at their northern end and a tilt toward the southeast at their southern end (Figure 4a). The slope, curvature, and roughness delineate the topographical changes, which are consistent with the undulations of the mountains (Figure 4b-d). The aspect has no distinctive spatial characteristics (Figure 4e). TWI is higher at the relatively low-lying terrain (Figure 4f). The distance from the road net indicates a higher density of road networks in the southern parts of the study region and, hence, a higher intensity of human activities (Figure 4h). Large rivers are all developed parallel to the faults and deep major fractures and are densely distributed throughout the study area (Figure 4g,i).
landslides. Static factors included nine explanatory variables (Figure 4). The topography of the Hengduan Mountains region is generally high in the north and low in the south, with remarkable differences in terrain, which cause unstable mountain slopes. The mountains run generally in a north-south direction between 25° and 30°N, with a slight tilt toward the west at their northern end and a tilt toward the southeast at their southern end (Figure 4a). The slope, curvature, and roughness delineate the topographical changes, which are consistent with the undulations of the mountains (Figure 4b-d). The aspect has no distinctive spatial characteristics (Figure 4e). TWI is higher at the relatively low-lying terrain (Figure 4f). The distance from the road net indicates a higher density of road networks in the southern parts of the study region and, hence, a higher intensity of human activities (Figure 4h). Large rivers are all developed parallel to the faults and deep major fractures and are densely distributed throughout the study area (Figure 4g,i).   (Table A1), and we analyzed the dynamic factors of the landslide and non-landslide sample points ( Figure 5). Before further analysis, we normalized the dynamic factors considering the diversity of dynamic factors and remarkable value ranges of the dynamic factors. The quartiles of dynamic factors of the landslides are larger than those of the non-landslides, implying large differences in the eigenvalues of the samples throughout the study region. For Geary's C spatial analysis, the values of both landslide and non-landslide factors mostly ranged from 0 to 0.1, which are much smaller than 0.5, indicating a strong spatial correlation of landslides and homogeneous features within the domain of landslides. Therefore, it is an ideal protocol for remote-sensing-based identification of landslide hazards in the Hengduan Mountains region.
tors of the landslides are larger than those of the non-landslides, implying large d ences in the eigenvalues of the samples throughout the study region. For Geary's C s analysis, the values of both landslide and non-landslide factors mostly ranged from 0.1, which are much smaller than 0.5, indicating a strong spatial correlation of land and homogeneous features within the domain of landslides. Therefore, it is an ideal pr for remote-sensing-based identification of landslide hazards in the Hengduan Mounta gion.

Performance Evaluation of LS Assessment Models
Results of the machine learning method are heavily reliant on the quality of the selected feature data and the merits of the algorithms. The evaluation of LS assessment models included two steps, i.e., comparison of various feature modes and comparison of different algorithm performances.

Feature Modes Comparison
The RF model was chosen to compare and analyze the landslide assessment results by different feature modes. To assess the importance of each feature mode, given a prescribed set of hyperparameters, three feature modes were used separately with the RF model to identify the landslides: (1) static mode, (2) dynamic mode, and (3) synergy mode. We used the hyperparameters numOfTrees = 200, variablesPerSplit = 'sqrt', and min_samples_leaf = 1 when implementing the model in GEE [100]. Figure 6 shows the accuracy and performance metrics of the LS assessment obtained with three different feature modes, i.e., static mode, dynamic mode, and synergy mode.
The validation results showed 0.9611 for the synergy model AUC, 0.9569 for the static model AUC, and 0.9419 for the dynamic model AUC. In addition, the accuracy indices calculated by the confusion matrix showed 90.48% in accuracy and 89.24% in precision for the synergy model, 89.74% in accuracy and 87.30% in precision for the static model, and 87.39% in accuracy and 83.80% in precision for the dynamic model. All these statistical metrics showed that the assessment model using synergy mode had higher evaluation accuracy and better performance. The results indicated that the paradigm of the leverage synergy mode combined with the machine learning algorithm based on the GEE platform for large-scale LS assessment was viable and feasible and had a satisfactory assessment performance.
scribed set of hyperparameters, three feature modes were used separately with the RF model to identify the landslides: (1) static mode, (2) dynamic mode, and (3) synergy mode. We used the hyperparameters numOfTrees = 200, variablesPerSplit = 'sqrt', and min_sam-ples_leaf = 1 when implementing the model in GEE [100]. Figure 6 shows the accuracy and performance metrics of the LS assessment obtained with three different feature modes, i.e., static mode, dynamic mode, and synergy mode. The validation results showed 0.9611 for the synergy model AUC, 0.9569 for the static model AUC, and 0.9419 for the dynamic model AUC. In addition, the accuracy indices calculated by the confusion matrix showed 90.48% in accuracy and 89.24% in precision for the synergy model, 89.74% in accuracy and 87.30% in precision for the static model, and 87.39% in accuracy and 83.80% in precision for the dynamic model. All these statistical metrics showed that the assessment model using synergy mode had higher evaluation accuracy and better performance. The results indicated that the paradigm of the leverage synergy mode combined with the machine learning algorithm based on the GEE platform for large-scale LS assessment was viable and feasible and had a satisfactory assessment performance.

Comparison of Landslide Identification Performances of Models
We determined that synergy mode with the machine learning algorithm was the right choice for further LS assessment practice. To do LS assessment with higher accuracy, we compared three types of classifiers for modeling performance, and the classifier with the highest modeling performance would be adopted for LS mapping over the study region.

Comparison of Landslide Identification Performances of Models
We determined that synergy mode with the machine learning algorithm was the right choice for further LS assessment practice. To do LS assessment with higher accuracy, we compared three types of classifiers for modeling performance, and the classifier with the highest modeling performance would be adopted for LS mapping over the study region. Figure 7 showed the accuracy and performance metrics of the LS assessment using three different classifiers, i.e., RF, CART, and SVM. The accuracy indices calculated by the confusion matrix showed 90.48% in accuracy and 89.24% in precision for the RF classifier, 84.80% in accuracy and 84.75% in precision for the CART classifier, and 72.20% in accuracy and 66.83% in precision for the SVM classifier. Evaluation metrics showed that the RF classifier had the best performance and highest accuracy when compared to CART and SVM. In addition, we perform LS mapping using CART and SVM classifiers to get the spatial performance. Through visual interpretation, it can be seen from LSM that there is a common overfitting phenomenon in CART and SVM ( Figure A1). In this case, we developed a LS assessment paradigm at the regional scale, i.e., the RF classifier coupled with synergy feature mode using the GEE platform for LS assessment. 84.80% in accuracy and 84.75% in precision for the CART classifier, and 72.20% in accuracy and 66.83% in precision for the SVM classifier. Evaluation metrics showed that the RF classifier had the best performance and highest accuracy when compared to CART and SVM. In addition, we perform LS mapping using CART and SVM classifiers to get the spatial performance. Through visual interpretation, it can be seen from LSM that there is a common overfitting phenomenon in CART and SVM ( Figure A1). In this case, we developed a LS assessment paradigm at the regional scale, i.e., the RF classifier coupled with synergy feature mode using the GEE platform for LS assessment.

Spatial Pattern of LS
The LS map obtained by the assessment methodology ( Figure 8a) indicated free LS in the northern part of the Hengduan Mountains region. The areas without LS colors mean the aggregation result value is zero. In addition, we take a zero LS value as the background, which indicates that landslides in these areas are unlikely to occur. These areas include, but are not limited to, flat land and bodies of water such as lakes ( Figure A3). Most high or very high LS levels can be attributed to unstable slope bodies due to weathering and/downcutting of the river channel, while downcutting processes can also be motivated by road construction. Besides, intensifying human activities also intensify the weathering processes of rocks and soils. The LS level is usually higher along the road than the surrounding areas, and potential human activities, mainly the construction of roads, have modified the stability of slope and induced landslides (Figures 9-14). We found higher LS in the southern parts of the study region when compared to that in the northwestern region, while population density and economic activities in the southern region had the dominant contribution to higher LS. The density of road network (Figure 4) also showed that the roads in the southern part of the study area were subject to higher density, and these regions matched well the regions with high LS, implying a considerable fractional contribution by human activities to LS in the study region.

Spatial Pattern of LS
The LS map obtained by the assessment methodology (Figure 8a) indicated free LS in the northern part of the Hengduan Mountains region. The areas without LS colors mean the aggregation result value is zero. In addition, we take a zero LS value as the background, which indicates that landslides in these areas are unlikely to occur. These areas include, but are not limited to, flat land and bodies of water such as lakes ( Figure A3). Most high or very high LS levels can be attributed to unstable slope bodies due to weathering and/downcutting of the river channel, while downcutting processes can also be motivated by road construction. Besides, intensifying human activities also intensify the weathering processes of rocks and soils. The LS level is usually higher along the road than the surrounding areas, and potential human activities, mainly the construction of roads, have modified the stability of slope and induced landslides (Figures 9-14). We found higher LS in the southern parts of the study region when compared to that in the northwestern region, while population density and economic activities in the southern region had the dominant contribution to higher LS. The density of road network (Figure 4) also showed that the roads in the southern part of the study area were subject to higher density, and these regions matched well the regions with high LS, implying a considerable fractional contribution by human activities to LS in the study region.            -d) show the similarities and differences of LS assessment results in static feature mode, dynamic feature mode, and synergy feature mode. (e-g) display details using Ersi and Google satellite image. Dynamic feature mode can sense detailed information of the pregnant environment, but the identification of floodplain fans and moraine landscapes can be easily confused. However, the synergy feature mode is beneficial to improve recognition accuracy. shows an recently constructed artificial slope protection building. Results show that using the dynamic feature mode can identify spatial texture structure information, and, when combined with static feature mode, it can reduce the false alarm rate and ensure the reliability of results. From the results of LS assessment in the study area, regions with high LS were in agreement with the spatial pattern of the landslide hazard density map produced from the investigation dataset by the Resource and Environment Science and Data Center of the Chinese Academy of Sciences (Figure 8b). It should be noted that the landslide poin density map could not locate every landslide hazard point. Here, we show a couple o possibilities. (1) The landslide hazard occurred a long time ago, and now we cannot iden From the results of LS assessment in the study area, regions with high LS were in agreement with the spatial pattern of the landslide hazard density map produced from the investigation dataset by the Resource and Environment Science and Data Center of the Chinese Academy of Sciences (Figure 8b). It should be noted that the landslide point density map could not locate every landslide hazard point. Here, we show a couple of possibilities. (1) The landslide hazard occurred a long time ago, and now we cannot identify it due because the traces of landslides disappeared due to long-term weathering and vegetation coverage changes. In this case, the landslide point density map might omit these obliterated landslide disaster relics. (2) The landslide/landslides occurred but this/these landslide event/events were not recorded. It could be because these landslide locations were inaccessible. Thus, the landslide locations recorded by the landslide point density map can be modeled by the LS assessment framework proposed in this study, implying that the modeling performance of the framework is acceptable. The proposed framework can locate other landslide locations that were not recorded in the landslide point density map. The density map of the training data used in the proposed framework is shown in Appendix B Figure A2. Hence, the modeling results can provide unique and irreplaceable information for enhancing mitigation of landslide disasters over the study region and can provide the reference information for the evaluation of landslide susceptibility in other regions of the globe.

Discussion
In this study, we developed a framework for assessing LS at the regional scale, allowing the efficient assessment of spatiotemporal variations in LS at a finer spatial resolution with the timely update of data, which is a kind of breakthrough into timely and systematic landslide susceptibility evaluation in regions with enormous difference in terrains. Based on the GEE platform and a synergy feature mode for comprehensively utilizing static and dynamic explanatory variables, our approach provides discernment of potential landslide locations, where landslides are subject to higher occurrence probability in a timely way. In this sense, the LS assessment framework can provide a viable technique for mitigation of landslide risk.
For machine learning methods for assessments of landslide susceptibility, the reliability of the model is closely related to the training sample, that is, the landslide inventory [14]. However, a landslide inventory cannot include all landslides occurred over a given region. Here, we used two inventories of landslide from different agencies to ensure acceptable representation of spatiotemporal pattern of landslides over the study region, i.e., expedition data were used for training samples and investigation data were used for comparison of the LS assessment map.
Availability of high-quality multisource and extensive data is a must to achieve an acceptable spatial assessment of LS [101]. Therefore, a synergy of static factors representing topography proxies and dynamic factors related to Earth surface properties and dynamic changes was used to identify landslide prone area over large-scale regions dominated by spatially heterogeneous geographical and geological features. Besides basic topographical and hydrometeorological features (Table A1), we further integrated various vegetation indices (i.e., NDVI, CRC, etc.) and spatial structure elements (Geary's C test) to derive landslide-related vegetation changes and diagnostic features [43,58]. Synergy feature mode combines the advantages of both types of data features, which reduces the noise error introduced by static feature mode or dynamic feature mode and can better represent the hazard-pregnant environment. Moreover, application of the proposed feature mode revealed that the LS assessment model using synergy mode achieved the highest modeling accuracies, when compared with the LS assessment model with static or dynamic feature mode only ( Figure 6). Then, we compared machine learning algorithms based on the synergy feature mode and found that the RF algorithm had a better performance and generalization ability and the modeling accuracy reached 90.48% and precision reached 89.24% (Figure 7). An improved framework integrating multiple explanatory variables with the fine spatial resolution ensured a reliable LS assessment with respect to the occurrence of landslide (Figure 9). Aforementioned procedures can help to produce a LS map using synergy feature mode, which can be put into practice (Figure 9a).
The LS assessment paradigm was designed to facilitate landslide research with a finer spatial resolution at the regional scale. We found that our modeling results were in agreement with the existing findings [7,15,59]. However, we can better identify the subtle spatial differences of the LS within the study area in a near-real-time way based on the GEE platform. Due to a large spatial range of the study region, the study region is dominated by discernible spatial heterogeneity in geomorphology, hydrometeorology, and triggering factors and the resulting distinctly variable LS in both space and time. Regions with a higher LS across the northern parts of the study region are spread along rivers, roads, and large fractures (Figure 8a), while a relatively high LS can be found over the southern part of the study area, which is attributed mainly to unstable slopes due to interferences of human activities such as an interwoven road network, a widespread residential area, and croplands ( Figure 10).
In high-latitude mountains, the stability of colluvial slope bodies is highly sensitive to external interferences or external forces, such as road construction, land reclamation, rainstorms, and so on. Different external forces may cause different degrees of instability of the colluvial slope bodies. In this case, we quantified percentages of landslide and non-landslide areas over regions with different land-use types. We found that nonlandslide areas were 4.41 times and 10.38 times more than the landslide areas in the forest and grassland land types, respectively, while the landslide areas were 3.66 times and 4.97 times higher than the non-landslide areas in the cropland and impervious surface land types, respectively ( Figure 11). These findings evidenced close relations between land use changes and LS in both space and time, illustrating a large spatial heterogeneity in LS variations that is highly influenced by forest cover types, cropland shifts, and settlement patterns [102][103][104].
Aiming to further verify and evidence LS assessment results, we did field surveys across the Hengduan Mountains region with a focus on typical landslides. Besides, we also attempted to obtain evidence from remote sensing images for additional evidence (Figures 9 and 11-14). Specifically, Figures 9, 12 and 13 describe the impact of use different feature mode on the LS results. Dynamic feature mode can sense detailed information of the pregnant environment and identify spatial texture structure information, but the identification of floodplain fans and moraine landscapes can be easily confused. However, the synergy feature mode is beneficial to improve recognition accuracy, reduce the false alarm rate, and ensure the reliability of results (Figures 12 and 13). Furthermore, we compare LS assessment results with landslide disaster information obtained from the internet media. The LS level of the area where the landslide occurred had high susceptibility, indicating that the quality of the LS assessment result is reliable ( Figure 14). We carefully compared evidence from field surveys and remote sensing images and found that the LS assessment results of the proposed LS evaluation framework predicted and/or simulated the landslide locations with satisfactory accuracy and high reliability. We also compared the LS map with the landslide inventory data density map (Figure 8a). Since the landslide inventory data were collected from different sources, they provided a great opportunity for cross validation. It was interesting to find that the LS was well consistent with the density map of landslide point in the spatial pattern. Ongoing LS studies at the regional scale will benefit from improved confidence in the LS assessment framework.

Conclusions
With the consideration of complicated driving factors behind landslide disasters and particularly the damaging impacts of landslides on the infrastructure, human settlements, and socioeconomy in the mountainous regions, we proposed a novel holistic landslide susceptibility evaluation framework to map the LS by a synergy of static, dynamic, and triggering factors using the GEE platform and RF model. Here, we considered the Hengduan Mountains region as a case study region, since it is dominated by remarkable differences in terrain (>7000 m in difference of terrain), different land-use patterns, and high sensitivity to weather extremes. Frequent landslide disasters and relevant disastrous consequences, par-ticularly, call for mitigation of landslide disasters in mountainous regions in the backdrop of warming climate. The proposed LS assessment framework involving landslide-related external and internal factors has been shown to be an effective technique for LS mapping in a near-real-time way, by routinely updating satellite data based on GEE's high-performance computing and abundant satellite data resources.
We mapped the LS over the landslide-prone Hengduan Mountains region with an area of~700,000 km 2 in Southwest China and found that this LS assessment framework worked well in mapping LS at a finer spatial scale and hourly time scale. We also showed reliability and accuracy of the LS map over the Hengduan Mountains region via comparison between multisource landslides inventory data, remote sensing images, and the LS map by the proposed framework, which showed a match with these data. Besides, landslides by field surveys also corroborated the modeling results. The results indicated that the spatial patterns of forests and grasslands had a significant controlling effect on the spatial patterns of the occurrence of landslides, while human interferences due to road construction, land reclamation, and pasturing further triggered instability of the colluvial slope bodies, which could easily trigger the occurrence of landslides within the Hengduan Mountains region. The LS map showed a complete picture of the susceptibility of landslide disasters, providing valuable information of the location of the landslide-prone regions, and is, therefore, beneficial for preparedness for and mitigation of landslide risk.
The LS assessment framework is transferrable and can be used in other landslide-prone regions in the globe. Landslide hazards in mountainous regions can be expected to occur with higher frequency due to intensifying weather processes, intensifying human interferences, and accelerating melting processes of the permafrost layers in high-altitude regions, such as the Hengduan Mountains region. Reliable and timely LS mapping is necessary for the preparedness for and mitigation of landslide disasters in a warming climate.

Acknowledgments:
We are thankful to the Bureau of Natural Resources in Yunnan, Sichuan, Qinghai, and Tibet for their assistance in the collection of the data analyzed in this study.