Landslide Susceptibility Assessment in Hong Kong with Consideration of Spatio-Temporal Consistency

Abstract

Landslide susceptibility is crucial for assessing the probability and severity of landslide disasters in a region. Previous studies have focused on static landslide susceptibility, using landslide assessment factor data from varying years, making it difficult to estimate spatio-temporal consistency and resulting in low prediction accuracy. Taking Hong Kong, China, as the study region, this study proposes a framework to estimate spatio-temporally consistent landslide susceptibility. The landslide assessment factors are divided into static and dynamic factors, with a temporal resolution of 5 years. Specifically, the dynamic assessment of landslide susceptibility is conducted for the periods 2000–2004, 2005–2009, 2010–2014, and 2015–2019, covering a total span from 2000 to 2019. Results show that the accuracy of the proposed model, defined as the proportion of correctly classified samples relative to the total number of samples, exceeds 0.7 across these four time periods. Both the F1-Score and the receiver operating characteristic (ROC) curve indicate that the proposed research framework exhibits good accuracy and practicality in susceptibility assessment. The proposed framework could capture temporal variations in landslide occurrence, allowing for a more accurate prediction of landslide susceptibility. The findings provide valuable insights for landslide disaster prevention and mitigation in Hong Kong and would also be applicable in other countries or regions.

1. Introduction

Landslide refers to the geological phenomenon of soil or rock sliding down a slope under the action of gravity [1]. Landslides often occur on steep slopes, and the movement can be either the overall sliding of the slope body or the local soil flow or detachment of rock strata. Landslides are typically triggered by factors such as geological structures, lithology, vegetation cover, rainfall, and earthquakes. Landslides can have severe impacts on nearby residents, transportation, infrastructure, and the environment, making the prevention, monitoring, and management of landslides crucial. The assessment of the probability of landslide occurrence from a spatial perspective, also known as landslide susceptibility assessment (LSA), aims to systematically evaluate various geological, topographical, climatic, vegetation, and anthropogenic factors to determine the extent of landslide likelihood in a specific region [2].

The calculation methods for landslide susceptibility can be primarily divided into knowledge-driven approaches and data-driven approaches [3,4]. Knowledge-driven methods can make use of the experiences of the experts [5]. Knowledge-driven assessment methods focus on utilizing experts’ subjective judgments and domain knowledge to explain the mechanisms of landslide occurrence and apply them in the assessment of landslide susceptibility. However, it is subjective, and the opinions of experts may be incorrect. Data-driven models are primarily based on the statistical analysis theory [6]. Data-driven methods refer to the utilization of large amounts of observational data and data analysis techniques such as machine learning. Through processing and analyzing the data, patterns and characteristics of landslide occurrence are directly extracted from the data to predict and evaluate the likelihood and hazard of landslide occurrence. Compared to traditional assessment methods based on expert knowledge and experience, data-driven methods place greater emphasis on obtaining information from actual data and can adapt to more complex geological environments and climatic conditions [7]. Commonly used data-driven methods include logistic regression (LR) and decision trees (DTs), which have been widely used in early assessments, along with kernel density estimation (KDE), which identifies high-risk areas based on historical data [8]. As computing technology has evolved, data-driven methods have advanced into more sophisticated modern machine learning techniques, such as support vector machines (SVMs), random forests (RFs), and artificial neural networks (ANNs) [9,10]. These methods have significantly improved LSA in regions like the Himalayas, Iran, and China, offering more accurate predictions and better hazard management solutions [11,12,13,14,15,16,17].

Sample selection in landslide susceptibility research is crucial. Huang et al. [18] used a non-supervised classification method to cluster landslide susceptibility areas and select non-landslide grid cells from extremely low-susceptibility zones. Liu et al. [7] employed fuzzy c-means clustering to improve the quality of “non-landslide” samples in the LSA model. Meanwhile, Zeng et al. [19] addressed the issue of unbalanced positive and negative samples that can arise from global random negative sampling strategies by utilizing graph neural networks. In terms of model training, various methods have been proposed. Jiping et al. [20] and Shuai et al. [21] introduced a multi-kernel SVM approach that considers optimized sample selection for LSA. Wang et al. investigated the application of convolutional neural networks (CNNs) in landslide susceptibility mapping in Yanshan County, China [22]. Wang et al. compared conventional recurrent neural networks with three of their variants, long short-term memory, gated recurrent unit, and simple recurrent unit, in a landslide susceptibility case study [23,24]. Wang et al. utilized a deep neural network architecture and interpreted the results of geological hazard susceptibility models through the spatial patterns of SHAP values [25]. Wei et al. combined infinite slope stability analysis with classic statistical algorithms such as logistic regression and CNN to create a hybrid model that enhances the practicality, interpretability, and cross-regional generalization capabilities of regional landslide susceptibility models [26].

Hong Kong is a region that frequently experiences landslide disasters. Numerous studies have assessed its landslide susceptibility. These studies rely on historical landslide points as samples and collect multi-source geospatial data as landslide assessment factors [27,28,29]. However, most of these studies employ relatively static LSA models, resulting in a loss of spatial and temporal correlations in landslide susceptibility predictions and insufficient understanding of landslide evolution patterns. Therefore, it is crucial to dynamically quantify regional landslide susceptibility across multiple time periods. A robust landslide susceptibility model can provide decision-making support for landslide disaster prevention and control in the region [30], while also enabling people to clearly understand the development patterns and future trends of landslides [31,32,33].

To address the aforementioned issues, this study takes Hong Kong, China, as the research area and proposes a dynamic LSA model that takes into account spatial and temporal consistency. The aim is to address the lack of spatial and temporal correlation in landslide susceptibility results and the insufficient understanding of landslide evolution patterns. The main contributions of this study are as follows: ① The results of LSA in Hong Kong from 2000 to 2019 were calculated with a time resolution of 5 years. This approach enables the model to capture temporal variations in landslide occurrence, allowing for a more accurate representation of the dynamic nature of landslide susceptibility. ② This study analyzes the 20-year evolution patterns of landslides in the region of Hong Kong. By examining the spatial distribution and temporal trends of landslide occurrences, this study provides insights into the factors that influence landslide susceptibility and how they change over time. This analysis not only enhances the understanding of landslide mechanisms but also contributes to the development of more effective landslide prevention and mitigation strategies.

Overall, the proposed dynamic LSA model offers a comprehensive approach to assessing landslide risk in the region of Hong Kong. By incorporating spatial and temporal consistency, this model provides a more accurate representation of landslide susceptibility and enables a deeper understanding of landslide evolution patterns. This information can be used to inform decision-making in landslide disaster prevention and control, ultimately contributing to the safety and resilience of the region [34,35].

2. Study Area and Dataset

2.1. Study Area

Hong Kong is located in the south of China, east of the Pearl River Estuary, with longitude and latitude ranging from 22.15° to 22.56° N and 113.84° to 114.44° E (as shown in Figure 1). It neighbors Shenzhen toward the north, faces Macao across the sea toward the west, connects to the Zhuhai Wanshan Archipelago toward the south, and faces the South China Sea toward the east. Hong Kong, Macao, and Guangdong jointly constitute the Guangdong–Hong Kong–Macao Greater Bay Area. The elevation range of the research area is 0–925 m, and the slope range is 0–60° (with an average slope of 12.12°). The terrain is mainly characterized by low-altitude mountainous areas with medium and small fluctuations, underwater deltas, and low-altitude hills. The vegetation coverage is generally high. Hong Kong has a subtropical monsoon climate with long summers and short winters. Summers are hot and rainy, while winters are warm and dry with less rain. This climate type makes Hong Kong’s temperature pleasant throughout the year, but there are also significant seasonal differences. Summer temperatures usually range from 26 to 30 °C, and Hong Kong experiences frequent and sometimes heavy rainfall during this season, accompanied by extreme weather such as typhoons. Winters (from December to February) are relatively cool, with temperatures generally ranging from 10 to 20 °C. It rarely drops below 5 °C, and the weather is relatively dry. Hong Kong is also influenced by a marine climate, resulting in high humidity throughout the year, with an average relative humidity of approximately 78%. Additionally, Hong Kong’s geographic location in a typhoon-prone area exposes it to potential typhoon impacts during the summer and autumn seasons. The climatic characteristics of Hong Kong determine that the recharge of groundwater mainly comes from atmospheric precipitation. A large amount of rainwater infiltrates through the surface and underground runoff, providing a rich source of recharge for groundwater. The stratigraphic lithology includes Quaternary deposits, Tertiary red clay sandstone and fan-shaped sandstone, as well as Permian granite and andesite [36].

Located in the hilly region of the southeast coast, Hong Kong has complex geological conditions with numerous mountains and hills. Coupled with the influence of monsoon winds, frequent summer rainstorms provide conditions for the occurrence of landslide disasters. As one of the most densely populated regions globally, Hong Kong faces a scarcity of land resources, leading to the construction of residential buildings and public infrastructure, often in hilly areas. This, in turn, increases the risk of landslide disasters. Landslides in Hong Kong are characterized by frequent occurrence, destructiveness, and unpredictability. Once a landslide occurs, it not only causes casualties and property damage but also severely impacts local transportation, communication, and other infrastructure. Typical landslide-prone areas include the hills of the New Territories in Hong Kong, as well as some residential areas and public facilities near mountainous regions. Take the landslide disaster on Po Shan Road as an example. On 18 June 1972, after continuous heavy rainstorms hit Hong Kong Island, large-scale landslides occurred consecutively on the slopes behind Po Shan Road, resulting in a significant landslide deposit and causing a severe accident with 67 fatalities, 20 injuries, the disruption of roads at the lower edge of the landslide, and the complete destruction of two high-rise buildings. The economic losses amounted to tens of millions of Hong Kong dollars. To prevent and control landslide disasters, the Hong Kong government has taken a series of measures, such as establishing the Geotechnical Engineering Office, specifically responsible for the design and maintenance of slopes; strengthening geological surveys and monitoring and early warning systems; and raising public awareness and prevention consciousness regarding landslide disasters. These measures have reduced the occurrence and losses caused by landslide disasters to some extent, but continuous reinforcement and improvement are still necessary.

2.2. Dataset

In this study, a five-year temporal resolution landslide susceptibility database covering the period from 2000 to 2019 has been constructed. The landslide samples were sourced from the Hong Kong Spatial Data Sharing Platform. The original data are provided as geographic coordinates (latitude and longitude), accompanied by attribute information, including disaster ID, year of occurrence, and the number of individuals affected, as illustrated in Figure 2. Specifically, there were 1535 historical landslide points from 2000 to 2004, 3206 from 2005 to 2010, 781 from 2010 to 2014, and 1399 from 2015 to 2019. Using the landslide samples (historical landslide points) as spatial references, an equal number of random points were sampled outside a 5 km radius and defined as non-landslide samples. The landslide samples and non-landslide samples collectively compose the landslide dataset.

Geospatial data relevant to landslide disasters within the study area were selected, including fixed factors such as elevation, slope, aspect, curvature, lithology, and topographic wetness index (TWI), as well as time-varying factors such as average rainfall, max normalized difference vegetation index (NDVI) values, and land cover [37]. Please refer to Figure 3, Figure 4, Figure 5 and Figure 6 for further details. Administrative boundary data for the Hong Kong study area were sourced from the CSDI Geoportal [38]. Among these, elevation, slope, aspect, curvature, TWI, and lithology data serve as static landslide assessment factors, while average annual rainfall (AAR), average annual maximum NDVI (AMN), and annual land cover (LC) are considered dynamic landslide assessment factors. Elevation data were sourced from the geospatial data cloud [39]. Slope, aspect, and curvature data were derived from the elevation data. Lithology data were obtained from the Hong Kong Civil Engineering and Development Department [40], and the legend also referred to this website. The data format is in vector polygons, but in this study, they were rasterized. Hydrological analysis was conducted on elevation and slope data to generate landslide assessment factors such as TWI. Rainfall data were sourced from the geospatial resource data cloud (GRDC) [41]. NDVI data were obtained from the National Ecosystem Survey Data Center (NESDC) [42], and the data were stored after being multiplied by 10,000 [43]. Land cover data were sourced from Yang Jie’s team at Wuhan University [44]. Among these time-varying landslide assessment factors, average rainfall and NDVI values were calculated for a five-year period. Land cover data were fine-tuned based on the real-world scenario, with the value taken as the most frequent occurrence within the five-year period. All the aforementioned data were resampled to a 30 m spatial resolution raster format and normalized for analysis.

3. Methodology

The workflow of LSA is illustrated in Figure 7. Multi-source geospatial data were collected as samples and assessment factors (consisting of dynamic and static factors) for LSA. The collinearity of the factors was evaluated using variance inflation factor (VIF), and RF and SVM were employed to predict the probability of landslide occurrence in the study area. The prediction results of the models were then mapped and analyzed. An enumeration of parameter combinations was used, as this approach is more refined and ensures that the optimal parameter combination for the study is achieved. The experimental environment was CPU: Intel i7-14700KF, GPU: RTX 4070 (Intel, Santa Clara, CA, USA), machine learning framework scikit-learn 1.3.1, and Python version 3.7.

3.1. Multicollinearity Analysis

When using landslide assessment factors as model inputs, the presence of high correlation or linear dependency among the input features can lead to inaccurate parameter estimation of the model, reduce its stability, and pose the risk of overfitting [7]. In this study, the VIF was employed to assess the presence of collinearity among the input features, as specified in Equation (1):

$$VI{F}_i={\displaystyle \frac{1}{1-R_i^2}}$$

(1)

where R² represents the correlation between a given variable among the independent variables and the remaining independent variables through multiple linear regression. The resulting model accuracy is used as a measure of the correlation between this variable and the remaining independent variables. When the VIF is close to 1, it indicates that there is no multicollinearity among the independent variables, and each variable independently explains the dependent variable, resulting in good model stability. When the VIF is greater than 1 but less than 10, it suggests a certain degree of multicollinearity, which may affect the interpretability of the model but not significantly impact its performance. However, when the VIF is greater than 10, it indicates severe multicollinearity, leading to unstable parameter estimation and potentially severe impacts on both the interpretability and predictive performance of the model.

3.2. Machine Learning

3.2.1. Random Forest

RF consists of a collection of multiple decision trees, each of which is generated by training on a random subsample of the data [45]. In RF, each decision tree is trained independently, meaning that different subsets of training data and features are used in the construction of each tree. This randomness helps reduce overfitting and improves the generalization ability of the model. The training process of RF involves randomness in two main aspects: data sampling and feature selection. Data sampling is achieved by randomly selecting a subsample of the training data with replacement to ensure that each decision tree is trained on a different subset of data. Feature selection involves randomly selecting a subset of features at each node of the decision tree for splitting. This randomness helps reduce the correlation between decision trees, thereby enhancing the performance of the overall model. During prediction, RF adopts a voting mechanism, where each decision tree predicts the outcome for a new sample, and the final prediction result is obtained by voting or averaging the predictions of all decision trees. This integrated approach effectively reduces the variance of predictions and improves prediction accuracy. RF has the advantage of being less sensitive to noise and overfitting and has the ability to handle high data dimensionality and multicollinearity. However, RF also has disadvantages such as high computational complexity and slow computation speed.

3.2.2. Support Vector Machine

SVM is a supervised learning algorithm primarily used for classification problems, but it can also be applied to regression problems [46]. In classification tasks, the main objective of SVM is to find a hyperplane to segment the samples, with the principle of maximizing the margin. Support vector classification (SVC) is an important application of SVM in classification problems. When the input labels are categorical, SVC can be used for classification. Its basic model is defined as a linear classifier with the largest margin in the feature space. In other words, the learning strategy of SVM is to maximize the margin, which ultimately can be transformed into solving a convex quadratic programming problem. SVM has several notable advantages, such as strong robustness for small samples with a large number of feature attributes, strong learning capabilities for both simple and complex classification models, and the ability to avoid overfitting when adopting complex mathematical models. Additionally, SVM has strong convergence performance, a simple structure, few adjustable parameters, and is easy to implement. However, SVM also has some drawbacks, such as being sensitive to parameter selection, lacking interpretability, having low convergence accuracy, and having insufficient convergence speed in some complex problems. In practical applications, SVM is widely used in areas such as image recognition and text classification and is particularly suitable for small samples, high-dimensional features, and nonlinear problems.

3.3. Model Validation

In machine learning classification tasks, accuracy measures the overall correctness of the classifier, while precision and recall focus more on the classifier’s performance on different categories. The F1-Score, which combines precision and recall, is especially useful for evaluating models in scenarios with different numbers of categories and unbalanced data, as referenced in Equations (2)–(5).

$$Accuracy={\displaystyle \frac{TP+TN}{TP+FP+TN+FN}}$$

(2)

$$Precision={\displaystyle \frac{TP}{TP+FP}}$$

(3)

$$Recall={\displaystyle \frac{TP}{TP+FN}}$$

(4)

$$F1-score=2\times {\displaystyle \frac{Precision\times Recall}{Precision+Recall}}$$

(5)

where true positive (TP) represents correctly predicted positive instances, true negative (TN) represents correctly predicted negative instances, false positive (FP) represents falsely predicted positive instances, and false negative (FN) represents falsely predicted negative instances. Additionally, the receiver operating characteristic (ROC) curve is used to demonstrate the relationship between the true positive rate (also known as the recall rate) and the false positive rate of the classifier at different thresholds. The area under the ROC curve, also known as the area under the curve (AUC) value, is a commonly used measure of a classifier’s ability to distinguish between positive and negative samples.

4. Results

4.1. Results of Multicollinearity Analysis

The results of the collinearity assessment for landslide assessment factors across four time ranges in this study are presented in

Table 1. VIF result.

	2000–2004	2005–2019	2010–2014	2015–2019
Elevation	1.614	1.454	1.742	2.027
Slope	3.337	3.574	3.287	3.585
Aspect	1.081	1.150	1.147	1.305
Curvature	1.032	1.160	1.120	1.038
TWI	3.957	4.001	3.962	3.782
Lithology	1.095	1.248	1.494	1.201
LC	1.768	2.434	1.356	1.497
AAR	1.402	1.289	2.089	2.085
AMN	1.873	2.681	1.504	1.515

. As introduced in the Section 3, when the VIF result of an assessment factor approaches 10, it indicates a certain degree of collinearity with other factors. In this study, the VIF values of the landslide assessment factors are all less than 5, indicating that there are no significant issues of multicollinearity among the assessment factors. This also suggests that the selection of landslide assessment factors in this study is relatively robust, which is conducive to improving the performance of the model.

4.2. Model Performance

In this study, a comprehensive evaluation of the overall performance of the models was conducted using accuracy, precision, recall, F1-Score, and ROC curve (AUC value). The landslide dataset, composed of both landslide and non-landslide samples, was divided into 70% for model training and 30% for model testing. Using the time range as the dividing line, the performance of the RF and SVM models on the test set was compared, as detailed in

Table 2. The 2000–2004 model performance.

	Accuracy	Precision	Recall	F1-Score	Parameters
SVM	0.7351	0.7154	0.7719	0.7426	{ kernal = ‘rbf’, C = 25, gamma = 20 }
RF	0.7980	0.7987	0.7917	0.7952	{ Best Estimators = 102, Best Max Depth = 15 }

,

Table 3. The 2005–2009 model performance.

	Accuracy	Precision	Recall	F1-Score	Parameters
SVM	0.7734	0.7623	0.7785	0.7703	{ kernal = ‘rbf’, C = 28, gamma = 18 }
RF	0.8046	0.8110	0.7817	0.7961	{ Best Estimators = 82, Best Max Depth = 20 }

,

Table 4. The 2010–2014 model performance.

	Accuracy	Precision	Recall	F1-Score	Parameters
SVM	0.7100	0.6774	0.8043	0.7354	{ kernal = ‘rbf’, C = 39, gamma = 29 }
RF	0.7740	0.7792	0.7660	0.7660	{ Best Estimators = 66, Best Max Depth = 8 }

and

Table 5. The 2015–2019 model performance.

	Accuracy	Precision	Recall	F1-Score	Parameters
SVM	0.7024	0.6747	0.7896	0.7277	{ kernal = ‘poly’, C = 13, gamma = 6 }
RF	0.8466	0.7786	0.7447	0.8015	{ Best Estimators = 16, Best Max Depth = 19 }

and Figure 8.

Essentially, this study is a binary classification task. The spatial probability analysis of landslide occurrence is a predictive analysis of the spatial scope of landslide occurrence, aiming to solve the problem of where landslides will occur. In this task, landslide areas are regarded as the positive class, while non-landslide areas are regarded as the negative class. Accuracy represents the proportion of correctly classified landslide and non-landslide areas by the model. Precision measures the proportion of areas predicted by the model as landslide areas that actually experience landslides. Recall focuses on the proportion of all areas that actually experience landslides that are successfully predicted by the model, reflecting the model’s ability to identify landslides. F1-Score is the harmonic mean of precision and recall, used to comprehensively evaluate the model’s performance. A higher F1-Score in this study indicates that the model performs well in both precision and recall.

Specifically, the RF model outperforms the SVM model in terms of accuracy, precision, and F1-Score across all four time ranges. In terms of recall, the RF model leads the SVM model in the time ranges of 2000–2004 and 2005–2010. Additionally, ROC curves were plotted for the four study time ranges. Consistent with the results shown by the above-mentioned machine learning evaluation metrics, the ROC curves of the RF model are generally higher than those of the SVM model. The AUC values of the RF model lead the SVM model by between 0.4 and 0.8 across all time ranges. Overall, the model accuracy across the four time ranges is higher than 0.7, and the F1-Score and ROC curves both indicate that the models and parameters selected in this study are suitable for the actual situation and exhibit good performance.

4.3. Landslide Susceptibility Mapping

The landslide susceptibility of the entire study area was evaluated using the aforementioned models. The probability of landslide occurrence at each spatial location was predicted pixel by pixel based on spatial indices and calculated the likelihood of it belonging to the landslide class, followed by spatial visualization. In ArcGIS 10.8 software, the Jenks statistical method was employed to classify the landslide susceptibility of the entire study area into five levels: very low, low, middle, high, and very high.

In the 2000–2004 time range (Figure 9a and Figure 10a), both the RF and SVM models predicted two very high areas in the central and western regions of the study area. Combined with Figure 2a, it can be observed that these areas have a relatively dense distribution of landslide points, which is consistent with the spatial distribution of the prediction results. Additionally, in the northeastern part of Lantau Island, the SVM model predicted a significantly higher landslide susceptibility level than the RF model. When compared with the distribution of historical landslide points, it indirectly reflects that the RF model is more consistent with the actual situation.

In the 2005–2009 time range (Figure 9b and Figure 10b), a severe landslide disaster occurred in Lantau Island in 2008, with more than 3000 recorded landslide points. Therefore, both the RF and SVM models predicted a very high landslide susceptibility in the Lantau Island area. In other regions, the prediction results of the RF and SVM models were relatively similar in spatial distribution. Landslide disasters were mainly concentrated on Lantau Island during this time period, while the probability of landslide occurrence in other areas of Hong Kong was moderately low.

In the 2010–2014 time range (Figure 9c and Figure 10c), the prediction results of the RF and SVM models again showed similar patterns in the northern area around the High Island Reservoir in the northeastern part of the study area. This region has a relatively high number of historical landslide points, which is consistent with the results of landslide susceptibility. In the northern area of the Plover Cove Reservoir, there are also a large number of historical landslide points, and the prediction results of both models are also similar to the spatial distribution pattern of historical landslide points.

In the 2015–2019 time range (Figure 9d and Figure 10d), historical landslide points were clustered in two clusters. The first cluster is located in the northeastern part of the study area, similar to the distribution of landslides in the 2010–2014 time range. The other cluster is located in the central–northern part of the study area, where the elevation is significantly higher than the surrounding areas. However, the prediction results of both models were not very accurate. The RF model overlooked the historical landslide points in the central–northern part of the study area, while the SVM model predicted most of the study area as having moderate to high susceptibility, which is unreasonable. It is conceivable that during the selection phase of non-landslide samples, the commonly used random sampling method was applied. This may have resulted in the inclusion of samples from areas with high susceptibility, potentially influencing the model’s accuracy. While the model’s performance is not without limitations, the integration with geospatial visualization data revealed a strong correlation between the predicted results and established geological knowledge. This correspondence indirectly supports the reliability and applicability of the research findings.

5. Discussion

5.1. Data Distribution

As shown in Figure 11, a statistical analysis of the Jenks classification results for all models was conducted. This analysis serves two purposes: firstly, to understand the data distribution and subsequently infer the overall occurrence probability of landslide disasters in the study area; secondly, to compare the Jenks classification results among different models, which facilitates the cross-validation between them.

Throughout the various stages of this study, the Jenks classification results of the prediction outcomes from the two models did not seem to align perfectly. Referring to the results of model performance evaluation, during the 2000–2004 period, SVM’s predictions had a higher proportion of low and moderate susceptibility areas, while RF’s predictions tended to have a higher proportion of the three middle categories. Given that RF outperformed SVM in overall performance during this period, combined with the analysis of the actual situation in Figure 9 and Figure 10, it is believed that the data distribution of RF’s prediction results is closer to the true values. For the intermediate ten years from 2005 to 2014, the distribution of prediction results from the two models was relatively similar. However, it is worth noting that in the 2015–2019 period, there was a significant difference in data distribution between SVM’s prediction results and RF’s. Based on the findings in Section 4.3, SVM’s predictions for this period cannot be considered reliable for assessing landslide susceptibility in the area.

5.2. Assessment Factor Importance Analysis

Each decision tree in the RF is trained using different samples and subsets of features to avoid overfitting issues. After the prediction results of all decision trees are integrated through methods such as voting or averaging, the final prediction is obtained. During this process, the frequency of use and prediction effectiveness of each feature across multiple decision trees are recorded. By summing up the feature importance scores of all trees, the importance of each feature within the entire RF model is determined.

Based on the discussion in Section 5.1, the performance of the RF model was specifically analyzed. The importance levels of assessment factors for different study time ranges were calculated. Figure 12 demonstrates that the lithology factor is the most significant driving factor in all tasks, which is consistent with the conclusion of Ermini et al. [47] that lithology is a classic variable that controls landslide susceptibility. Additionally, AMN and curvature play a negligible role in driving landslides in this study, and AMN’s importance level is even zero during the 2010–2014 period.

By comparing the assessment factors with the prediction results of the RF model in Figure 4, a certain degree of overlap is observed between the spatial distribution of middle-to-high susceptibility areas and regions with high rainfall. In other words, high susceptibility areas often correlate with greater rainfall, which is consistent with the work by Turner and Schuster [48], Wong [49], Froude and Petley [50], and Ip et al. [51]. But this does not necessarily mean that places with high rainfall will always have middle-to-high susceptibility. Regarding the lithology factor, taking the high-susceptibility area of Lantau Island as an example, its stratigraphy mainly comprises Jurassic volcanic rocks with an age ranging from approximately 147.5 ± 0.2 to 146.6 ± 200,000 years. The lithology is primarily composed of fine-grained glassy tuff and flow-banded rhyolitic lava, with a small amount of euhedral coarse-grained crystal tuff. Given the analysis results from the past 20 years, Lantau Island has consistently been a high-risk area for landslides. Therefore, there is reason to believe that Lantau Island will continue to be a prime focus for landslide disaster prevention and control in Hong Kong for the foreseeable future.

5.3. Typical Region Analysis

Correlation analysis was conducted between the prediction results of this study and the three-dimensional landscapes in the real world, aiming to achieve two purposes. The first is to verify the validity, reliability, and usability of the research outcomes in the real world. The second is to observe the dynamic changes in landslide susceptibility by comparing data from different time periods in the same region.

5.3.1. Comparison of Boundary Areas of Different Terrains

Kowloon, located in the central part of Hong Kong, faces the sea in the south and is bordered by Beacon Hill and Lion Rock in the north. It is a typical area where mountains and cities intersect. Common sense tells us that the plain areas of the city are not affected by landslides, while mountainous areas with steep slopes are prone to landslide disasters. The prediction results were superimposed with the three-dimensional perspective of this region, as shown in Figure 13 (the large image on the left is the real image of the region, and the four sub-images on the right show the superimposed results for different time ranges. The composition of all images in Section 5.3 follows the same idea). The legend is consistent with Figure 9, indicating that the red areas are moderately to highly susceptible to landslides. As can be seen from the figure, the results are in good agreement with the actual terrain of the city and are fully applicable to this region.

5.3.2. Comparison of Interface Areas of Different Land Covers

As shown in Figure 14, the results reflect the differences between the reservoir shoreline and the reservoir itself, demonstrating the authenticity and accuracy of the findings. That is to say, landslides may occur at the intersection of water and land, but never in the middle of the water body (only the geographical surface is studied here, and submarine landslides are not considered in this study).

5.3.3. Areas with Significant Changes in Landslide Susceptibility

Figure 15a,b depicts the western part of the study area, where a significant number of historical landslide points were clustered between 2000 and 2004, forming a contiguous high-susceptibility area for landslides. Similar to Figure 13, the prediction results accurately match the local conditions in the plain residential areas between the two mountains. Over time, the susceptibility level gradually decreases in the mountainous areas on both sides of the residential areas, which is also consistent with local conditions.

Figure 16 shows Lantau Island in Hong Kong [52], which has been mentioned multiple times before. The overall trend here differs from that in Figure 15. In 2008, a severe landslide disaster occurred here. Although landslides occurred in other years, none were as severe as the one in 2008. Since 2008 falls within the time range of 2005–2009, the results effectively reflect this temporal difference.

Overall, the superiority of the research results was verified through the comparative analysis of susceptibility outcomes from different perspectives, including various terrain interfaces, surface cover areas, and time ranges.

5.4. Limitations and Future Work

This study employed SVM and RF to assess landslide susceptibility in Hong Kong across four time periods. A limitation of this study is the random sampling method used to select non-landslide samples, which may have included areas with higher susceptibility, potentially affecting model accuracy. Additionally, the performance of SVM for the 2015–2019 period was slightly hindered by a limited number of landslide samples in certain years, impacting its ability to generate precise temporal susceptibility maps. However, visual inspection of the results showed good consistency with geological knowledge and existing datasets, supporting the model’s reliability.

To improve model performance, future research will integrate interferometric synthetic aperture radar (InSAR) technology to monitor ground subsidence and enhance the landslide sample database. Additionally, more advanced evaluation techniques, such as root mean square error (RMSE), mean square error (MSE), and confusion matrices, will be employed [53]. These methods will help assess model performance more robustly, optimize results, and facilitate the broader application of machine learning models in landslide susceptibility mapping.

6. Conclusions

In this study, multi-source geospatial data were collected to assess the landslide susceptibility in Hong Kong over a period of 20 years. Specifically, historical landslide data from 2000 to 2019 in Hong Kong were used as labels, while elevation, slope, aspect, curvature, TWI, and lithology were considered as static assessment factors for landslide risk. Additionally, AAR, AMN, and LC were identified as dynamic assessment factors. VIF was used to determine the collinearity among the assessment factors, and RF and SVM models were employed to predict the landslide susceptibility in the study area over a 20-year period. By conducting a comparative analysis of the results and discussing them in the context of local conditions, the accuracy of this study was fully validated. The main conclusions of this study are as follows:

① The prediction accuracy of the RF model is superior to that of the SVM model, and the stability and generalization ability of the RF model are also better than those of the SVM. Based on the research results obtained using the RF model, the findings can serve as a reference for landslide disaster prevention and control in Hong Kong.

② Lithology, AAR, and elevation were the primary driving factors for landslide disasters in Hong Kong from 2000 to 2019, and the type of future landslide disasters in Hong Kong is still expected to be rainfall-induced.

③ The LSA system proposed in this study, which takes into account spatio-temporal consistency, has been well validated in terms of accuracy and applicability.

The research findings focus on LSA and can provide critical spatial insights for disaster prevention, government decision-making, and public safety [54]. They help develop early warning systems and mitigation strategies, such as guiding land use planning and infrastructure development in stable areas.