Graph-Represented Broad Learning System for Landslide Susceptibility Mapping in Alpine-Canyon Region

Abstract

Zhouqu County is located at the intersection of two active structural belts in the east of the Qinghai-Tibet Plateau, which is a rare, high-incidence area of landslides, debris flow, and earthquakes on a global scale. The complex regional geological background, the fragile ecological environment, and the significant tectonic activities have caused great difficulties for the dynamic susceptibility assessment and prediction of landslides in the study area. Specifically, Zhouqu is a typical alpine-canyon region in geomorphology; currently there is still a lack of a landslide susceptibility assessment study for this particular type of area. Therefore, the development of landslide susceptibility mapping (LSM) in this area is of great significance for quickly grasping the regional landslide situation and formulating disaster reduction strategies. In this article, we propose a graph-represented learning algorithm named GBLS within a broad framework in order to better extract the spatially relevant characteristics of the geographical data and to quickly obtain the change pattern of landslide susceptibility according to the frequent variation (increase or decrease) of the data. Based on the broad structure, we construct a group of graph feature nodes through graph-represented learning to make better use of geometric correlation of data to upgrade the precision. The proposed method maintains the efficiency and effectiveness due to its broad structure, and even better, it is able to take advantage of incremental data to complete fast learning methodology without repeated calculation, thus avoiding time waste and massive computation consumption. Empirical results verify the excellent performance with high efficiency and generalization of GBLS on the 407 landslides in the study area inventoried by remote sensing interpretation and field investigation. Then, the landslide susceptibility map is drawn to visualize the landslide susceptibility assessment according to the result of GBLS with the highest AUC (0.982). The four most influential factors were ranked out as rainfall, NDVI, aspect, and Terrain Ruggedness Index. Our research provides a selection criterion that can be referenced for future research where GBLS is of great significance in LSM of the alpine-canyon region. It plays an important role in demonstrating and popularizing the research in the same type of landform environment. The LSM would help the government better prevent and confine the risk of landslide hazards in the alpine-canyon region of Zhouqu.

1. Introduction

Zhouqu County and its surrounding areas are located in the rapid deformation zone of the transition from the Qinghai-Tibet Plateau to the Loess Plateau. The landform is positioned in a typical alpine-canyon region [1]. Coupled with the complex geological environment conditions and intensive tectonic activities, this area is a rare high-incidence area for landslides, debris flows, and other geological disasters in China [2,3]. The large number and massive scale of geological hazards pose a huge threat to the local residents. In 2010, the “8.8” debris flow (SanYanyu catchment, Figure 1) disaster in Zhouqu County caused 1765 deaths along with huge damage to both the infrastructure and property [4]. In 2018, another large-scale landslide occurred in the lower reaches of Zhouqu County, causing a blockage of the Bailong River and an interruption of National Highway 345. In addition to the natural shearing movement of the slope caused by gravity, earthquake, rainfall, and other factors, the occurrence of landslides of Zhouqu County is also affected by geological changes in the surrounding area. In 2008, the Wenchuan earthquake induced or reactivated massive landslides in the area [5] and caused huge losses [6,7]. Therefore, estimating and measuring the landslide susceptibility in the Zhouqu area is critical to decrease or avoid disasters caused by landslides.

In this article, we aim to assess the landslide risk in the study area by visualizing the estimation through the landslide susceptibility mapping (LSM). Landslide susceptibility refers to the quantitative or qualitative assessment of the spatial distribution by category and quantity (or area) of existing or possible landslides in the area, which is a critical method for disaster reduction and mitigation of landslides and debris flows [8,9,10,11,12,13]. Several approaches have been proposed for landslide susceptibility, including geomorphological mapping, landslide inventories, heuristic terrain and susceptibility zoning, physically based numerical modeling, and statistical-based methods [14,15].

Statistical methods mainly focus on the correlation among the key factors, including topography, geological environment, climate, hydrology, vegetation, soil, and human activities, which affect the development and the spatial distribution of landslides [16]. However, the estimation and prediction via statistical methods rely on the strong hypothesis of data distribution, which may not hold true for real data. In recent years, machine learning is increasingly being used for the LSM with good results [17,18] attributed to the continuous development of AI technology and the improvement of computational hardware. As one of the generalized linear models, logistic regression (LR) [19,20,21] is a basic binary classifier estimating the possibility of occurrence, which is easy to understand and implement with low computational cost. However, LR performs poorly when there exists missing data, which is also severely sensitive to features and prone to fall into underfitting. Bivariate statistical-based kernel logistic regression (KLR) [22] is improved for higher accuracy, which is highly influenced by the fitness of the selected kernels, and the kernel parameters are also hard to optimize. The hybrid method may coordinate and balance the difficulty in choosing, but it brings more parameters to be decided upon. AdaBoost (AB) [23,24] shows remarkable performance in classification tasks, but it often falls into overfitting on the noisy sample sets and being sensitive to outliers. Stochastic gradient decent (SGD) [25] is a fast gradient decent in which only one training sample is used when the parameter is updated. Along the opposite direction of the corresponding gradient (or approximate gradient) of its current point, SGD iteratively searches for the local minimum of the objective function given the specified step size. If the sample is noisy, it easily falls into a local optimal solution and converges to an undesirable state. Random forest (RF) [26] is a classifier which is established in a random way and contains multiple decision trees. RF is one of the most popular classifiers that is applied well in many research areas; it can deal with large-scale data while assessing feature importance. Yet RF is not stable with randomly generated trees, so its performance is severely influenced by the noise. Support vector machine (SVM) [27] is a novel small sample learning method based on a solid statistical foundation. Its goal is to classify the feature space into the optimal hyperplane whereas the core of the method is the idea of maximizing the classification margin. SVM can be well applied to high-dimensional data, avoiding the problem of dimensional disasters, and it also can handle the interaction of nonlinear features. Because the computation complexity is $\mathcal{O}\left(N^3\right)$, the efficiency of the model will be rapidly reduced when the amount of observed data N is large. SVM suffers from overfitting because of noise data as well. The extreme learning machine (ELM) [28] is a single hidden layer feedforward neural network learning algorithm. This algorithm only needs to set the number of hidden layer nodes in the network, adjusting neither the input weights of the network nor the bias of hidden elements in the execution process. It can generate the unique optimal solution, so the learning speed is fast with the fair generalization. In most cases, good performance can be obtained for ELM, but the initial parameters of the hidden layer (the connection weight offset value and the number of nodes) still have a great influence on the classification accuracy. Inappropriate parameters will lead to relatively poor results, and the learning performance of a single ELM is unstable.

Some deep learning methods, such as convolutional neural network (CNN) [29,30] and recurrent neural network (RNN) [31,32], have also been applied to the landslide susceptibility mapping [33,34] and gained excellent performance. Generative Adversarial Networks (GAN) is an unsupervised generative model which has recently been applied to correct imbalanced landslide data sets [35]. These methodologies showed great performance in prediction compared to the traditional methods, but the high computation complexity and parameter adjustment also caused much more time consumption compared with the limited performance improvement that deep learning methods could bring. Even worse, it was difficult to give a reasonable explanation for the complex structure of the model. Wang et al. [36] combined two transfer-learning strategies for landslide susceptibility modeling to make hazard-related information available to stakeholders in a timely manner, but transfer learning, whose performance is overly dependent on pretrained models, is not a general manner that requires a very large data base.

The broad learning system (BLS) [37,38,39,40] is a flat network considering the mapped features as the input of the random vector functional link neural network (RVFLNN) [41] to complete the process of feature selection and extraction. In contrast to other classifiers, BLS eliminates the shortcomings of the long training process and provides an outstanding generalization ability. BLS proves to be efficient and effective in the classification scenario and it has been widely used in various fields, such as control, engineering, bio-medicine, etc. However, the methodology is proposed generally to scope with all types of data without specific refinements for landslide susceptibility mapping, which causes less accuracy than expected and reduces its effectiveness. Therefore, we propose an improved model in the broad structure, called the graph-represented broad learning system (GBLS), which is modified specifically for landslide characteristics in the alpine-canyon region. The proposed model is supposed to improve the prediction accuracy of landslide sensitivity by embedding geometric structural features of the data to the network while maintaining the effectiveness and efficiency in the broad framework. In particular, GBLS can be quickly reconstructed in the incremental way without retraining the whole model for additional data due to its structure, which is demonstrated theoretically and empirically in the following section. In this way, the proposed method is functional in practice to deal with data in the alpine-canyon region or the same type of landform environment where the rapid and effective landslide susceptibility assessment is important for quickly grasping the regional landslide hazard situation and formulating disaster reduction strategies.

The rest of the article is organized as follows. In Section 2, the study area is described and analyzed with remote sensing interpretation and feature selection. Then, GBLS is constructed theoretically after the fundamental structure of the standard BLS model is introduced roughly; the incremental learning property is also illustrated. In Section 3, experimental results and comparisons among several state-of-the-art models are displayed, followed by the landslide susceptibility mapping analysis of Zhouqu County and its surrounding areas through the newly proposed GBLS. Detailed discussions are conducted in Section 4, including model comparison, feature ranking, resolution sensitive analysis, etc. Finally, Section 5 draws the conclusion.

2. Method

2.1. Study Area Overview

Zhouqu County lies on the Diebu-Wudu thrust nappe structural belt as the core of Bailong River compound anticline, adjoining the margin of the Qinghai-Tibet Plateau in the west and the West Qinling orogenic belt in the north. Significant tectonic activities in the area have resulted in the development of fault structures and the formation of various tectonic surfaces of varying scales, which provides sufficient habitat for landslides and collapses. The study area also has various lithologic strata with complex structures, including Silurian, Devonian, Carboniferous, Permian, Triassic, Jurassic, Cretaceous, and Quaternary materials (Figure 1). In particular, the Silurian phyllite, slate, and schist are highly sensitive to the development of landslides. Complex topographical factors and diverse geological structures promote the occurrence of landslides.

In terms of climate, the study area is located in the western Qinling subtropical northernmost area. It generally belongs to the humid warm temperate zone, exhibiting a notable monsoon climate, with an annual total precipitation of 420.6 mm, 362.8 mm of which occurs from May to October, accounting for 86.3% of the annual precipitation. The uneven distribution of rainfall makes it easy to form a wide range, long, and heavy rainstorm. The annual average temperature is 13.4 °C whereas the monthly average temperature ranges from the lowest 1.8 °C to the highest 23.3 °C. Zhouqu County is a typical alpine-canyon region [2] where most of the landscape is mountainous, exceeding 3200 m in elevation, and the relative discrepancy in elevation to the valley is greater than 800 m. Due to the large difference in altitude, the temperature reached an extreme maximum of 38.2 °C and a minimum of $-10.2$ °C in recent years.

Influenced by the long history of human activities, the ecological environment in the study area becomes more and more fragile. It has suffered frequent geological hazards in recent years, for instance, the most memorable one was the 2010 Zhouqu County “8.8” debris flow [5], and several large events occurred in the past consecutive years: the 2018 Zhouqu Jiangdingya landslide, the 2019 Zhouqu Yahuokou landslide [42], and the 2020 Southeast Gansu flash flood secondary geological disasters, resulting in huge casualties and property losses. In consequence, this is particularly urgent and important for the risk assessment of disasters.

2.2. Remote Sensing Interpretation of Landslide

The relevant landslide data in the study area are collected and interpreted visually based on multi-phase Google Earth satellite remote sensing image. According to the geomorphological characteristics of the landslide, the boundary of the landslide is determined by comparing images from different time periods, and then to better calibrate and validate the interpreted landslide inventory, several field explorations are carried out to verify the consistency against the survey data. Finally, feature information of 407 landslides is defined and applied to our model in the following research (Figure 2).

2.3. Acquisition and Description of Landslide Feature Factors

One of the difficulties of landslide susceptibility and risk assessment is the division of evaluation units. In this article, the grid division method is presented, which is a benefit for matching and analysis of remote sensing data with other base maps, such as elevation, land use, and lithology. The resolution of the grid size is chosen in the spirit of balancing the redundancy of the raster and computational efficiency, and the 50 m × 50 m grid is adopted.

The study area is located in a mountainous area that contains a major streamflow and extensive records of accumulated layer landslides. As we have been investigating the causes of a series of colluvial landslides in the study area in depth, such as the Jiangding Cliff landslide and the Yahuokou landslide in Zhouqu County [18], the landslide and its surrounding areas are always located in the low hills and gullies, with great topographic fluctuations and steep natural slopes. The old landslide area forms a relatively gentle platform, which is conducive to rainwater collection and leakage as well as the formation of landslides. Furthermore, the surface water in the gully erodes the slope foot, causing the edge area of the slope foot to be further reduced, and even forming an anticline angle, which has a disturbance effect on the slope foot of the leading edge of the landslide and reduces the stability. Generally, the special topography, geomorphology, and geotechnical structure create favorable conditions for landslide formation in Zhouqu County.

Thus, several features are selected to capture the geomorphological characteristics and the interactive trait of water runoff and terrain traits. Slope has a controlling effect on the landslide and it significantly affects the stress distribution of the rock and soil [43]. Aspect directly affects the distribution of solar radiation and water vapor and causes regular differentiation of factors such as precipitation, vegetation, erosion, and weathering [44]. Plane curvature (Cpl) and Profile curvature (Cpr) reflect the unevenness of the terrain, which directly affect the collection of surface runoff and then affect the distribution of pore water pressure; both features can be used to understand the process of erosion and runoff formation. Roughness (R) and Terrain Ruggedness Index (TRI) measures the smoothness of terrain surface texture by reflecting the largest and mean value of the inter-cell differences in elevation [45], as well as the Topographic Position Index (TPI) to describe the difference between the average value of the center pixel and its surrounding cells [45]. Topographic Wetness Index (TWI) reflects the impact of local topography on runoff flow direction and accumulation [46]. Stream Power Index (SPI) measures water erosivity [47] and Sediment Transport Index (STI), comprehensively characterizing the transport and deposition of surface materials with flow.

Rainfall data measured by governmental meteorological stations and hydrological stations from 1964 to 1997 is gathered and processed, the average rainfall values of 10 min (R10m), 1 h (R1h), 6 h (R6h), 1 day (R1d), 3 days (R3d), and 1 year (R1y) are used for this analysis. Different lengths of the moving average on precipitation features not only enable the analysis to capture the short-term, midterm, and long-term moisture level and possible runoffs but also can provide a clear indicator of landslide early warning if any of the precipitation features prevail as a significant feature in the analysis results. By combining the geomorphological features and the precipitation information, we hope to capture how precipitation and terrain impact the susceptibility of landslides, especially for the accumulated layer landslides, and we provide results that can be interpreted without an investigation on the triggering rate of precipitation.

Other features we used include geology, meteorology, and hydrology, such as NDVI, land use, lithology, and distant to road and river. These features are chosen due to their possible impacts on slope stability [13,48,49,50]. NDVI is calculated from the image data of Gaofen-1 (8m resolution). Land use implemented a 5-level indicator based on 1:250,000 map published by China Academy of Sciences (CAS), to measure the extent of human activities and inhabitants. Similarly lithology is divided into 5-level labels from the 1:200,000 China official geological map according to the hardness of the rock [11]. Distance from roads and rivers is an important driving factor affecting slope stability and it reflects human activities as well.

Afterwards, we consider a large variety of factors that are associated with the development and occurrence of landslides and specify 21 features from the aspect of topography, geological environment, meteorological, and hydrological conditions.

Table 1. Feature Factors Description and Data Resources.

No.	Features	Abbr.	Unit	Data Resource
1	Slope	S	°	Calculated by DEM with resolution of 12 m (JAXA/METI ALOS PALSAR L1.0 2011).
2	Aspect	A	°
3	Plane curvature	Cpl	/
4	Profile curvature	Cpr	/
5	Roughness	R	/
6	Terrain Ruggedness Index	TRI	/
7	Topographic Position Index	TPI	/
8	Topographic Wetness Index	TWI	/
9	Stream Power Index	SPI	/
10	Sediment Transport Index	STI	/
11	Average rainfall of l0 min	R10m	mm	Obtained by interpolation of meteorological station data from 1964 to 1997.
12	Average rainfall of l h	R1h	mm
13	Average rainfall of 6 h	R6h	mm
14	Average rainfall of 1 day	R1d	mm
15	Average rainfall of 3 days	R3d	mm
16	Average rainfall of 1 year	R1y	mm
17	Normalized Difference Vegetation Index	NDVI	/	Calculated by Gaofen-1 image with resolution of 8 m (2013).
18	Land use	/	/	1:250,000 map published by China Academy of Sciences (2000).
19	Lithology	/	/	1:200,000 China official geological map (2002).
20	Distant to road	/	km	Obtained by image interpretation Gaofen-1 image with resolution of 8 m (2013).
21	Distant to river	/	km	Calculated by DEM with resolution of 12 m (2011).

details factor characteristics and data sources used in this article.

Illustrations on distribution of some major features are drawn by ArcGIS (Figure 3).

2.4. Model

2.4.1. Preliminary of Broad Learning System

In brief, the structure of standard BLS is composed of three parts. Given data $\{\mathbf{X},\mathbf{Y}\}\in {\mathbb{R}}^{N\times (M+C)}$, the ith group of mapped features is

$${\mathbf{Z}}_i\stackrel{\Delta }{=}{\varphi }_i({\mathbf{XW}}_{e_i}+{\mathbf{\beta }}_{e_i}),i=1,\dots ,n.$$

(1)

For simplicity, we denote the collection of n groups of mapped features as ${\mathbf{Z}}^n\stackrel{\Delta }{=}\left[{\mathbf{Z}}_1,{\mathbf{Z}}_2,\dots ,{\mathbf{Z}}_n\right]$, which also forms the input of enhancement layer. Similarly, the outputs of the jth group of enhancement nodes are denoted by

$${\mathbf{H}}_j\stackrel{\Delta }{=}{\xi }_j({\mathbf{Z}}^n{\mathbf{W}}_{h_j}+{\mathbf{\beta }}_{h_j}),j=1,\dots ,m.$$

(2)

and set ${\mathbf{H}}^m\stackrel{\Delta }{=}\left[{\mathbf{H}}_1,{\mathbf{H}}_2,\dots ,{\mathbf{H}}_m\right]$ in the same way. The weights ${\mathbf{W}}_{e_i},{\mathbf{W}}_{h_j}$ and their corresponding bias ${\mathbf{\beta }}_{e_i},{\mathbf{\beta }}_{h_j}$ are randomly generated with matching dimensions. Note that ${\xi }_j$ are activation functions connecting feature nodes to enhancement nodes to introduce nonlinear correlation into the system. Let $\mathbf{A}=\left[{\mathbf{Z}}^n,{\mathbf{H}}^m\right]$; then, the output $\mathbf{Y}=\mathbf{AW}$ shows the predicted classification results with the output-layer weight $\mathbf{W}$ computed by the ridge regression approximation of pseudoinverse [51] as

$$\mathbf{W}={(\lambda \mathbf{I}+\mathbf{A}{\mathbf{A}}^T)}^{-1}{\mathbf{A}}^T\mathbf{Y},$$

(3)

where $\lambda$ denotes the ${\ell }_2$-norm regulator constraining the sum of squared weights. Specifically,

$${\mathbf{A}}^{+}=\underset{\lambda \to 0}{lim}{(\lambda \mathbf{I}+\mathbf{A}{\mathbf{A}}^T)}^{-1}{\mathbf{A}}^T.$$

(4)

Complete deduction and principles can be checked in [37].

2.4.2. Graph-Represented Broad Learning System

In the flat single structure of BLS, the mapped features ${\mathbf{Z}}^n$ are extended to the enhancement nodes ${\mathbf{H}}^m$. Both mapped feature nodes and enhancement nodes are directly linked to the last layer of model and the output weight is estimated by fast pseudo-inverse. Due to its structure, BLS is a kind of neural network that does not depend on the depth structure, which is suitable for the system with few features but has a high requirement for real-time prediction. However, standard BLS is not capable of extracting and processing correlation of geometric structure between samples, which is very important for landslide data. Therefore, we aim to build a new model based on the broad framework that inherits the efficient computational performance and incremental learning capability of BLS while exploiting the spatial geometric correlations of the data.

It is reasonable to assume that the input sample space is composed of multiple low-dimensional manifolds on which all samples lie and data points lying on the same manifold have the same label from a novel geometric perspective. This manifold assumption is popular and wildly used in many semi-supervised learning algorithms. Demonstrated by the manifold learning [52] and spectral graph theory [53], the data manifold can be effectively approximated through the graph Laplacian matrix $\mathcal{L}$, which is constructed in an essential way by

$$\mathcal{L}=\mathcal{D}-\mathcal{W},$$

(5)

where $\mathcal{D}=diag(d_1,\dots ,d_N)$ is the degree matrix,

$${\mathcal{D}}_{ij}=\left\{\begin{array}{cc}\sum _j{\mathcal{W}}_{ij}\hfill & \hfill \mathrm{if}i=j\\ 0\hfill & \hfill \mathrm{if}i\ne j\end{array}\right.$$

(6)

The adjacent graph matrix $\mathcal{W}$ is obtained by a k-nearest neighbor graph with Gaussian kernel as

$${\mathcal{W}}_{ij}=\left\{\begin{array}{cc}{𝓈}_{ij}\hfill & \hfill \mathrm{if}{\mathbf{x}}_i\in \mathrm{KNN}\left({\mathbf{x}}_j\right)\mathrm{and}{\mathbf{x}}_j\in \mathrm{KNN}\left({\mathbf{x}}_i\right)\\ 0\hfill & \hfill \mathrm{if}{\mathbf{x}}_i\notin \mathrm{KNN}\left({\mathbf{x}}_j\right)\mathrm{or}{\mathbf{x}}_j\notin \mathrm{KNN}\left({\mathbf{x}}_i\right)\end{array}\right.$$

(7)

where ${𝓈}_{ij}=exp\{-\frac{\parallel {\mathbf{x}}_i-{\mathbf{x}}_j{\parallel }^2}{h^2}\}$ measures the similarity of sample ${\mathbf{x}}_i$ and its k nearest neighbor $\mathrm{KNN}\left({\mathbf{x}}_i\right)$. The Gaussian bandwidth h is calculated by the average distance between ${\mathbf{x}}_i$ and its kth neighbors [54], whereas k is resolved by $min\{6,logN\}$ where N is the number of samples [55].

Usually $\mathcal{L}\in {\mathbb{R}}^{N\times N}$; then, the computation of $\mathcal{L}$ might run out of memory given a large sample size. Note that adjacent samples are more likely to have feature correlation with each other. It is rational to divide the data from all study areas into several groups, each of which contains equal sizes of samples. If the sample size is not divisible, we add adequate zero vectors to the undersized group to balance the dimensions. Then, we conduct the graph Laplacian matrix on each group of samples to gain the geometric features of data, and then the graph feature nodes are combined with those block Laplacian matrices. In this way, we manage to extract manifold features and avoid too much computational complexity on large-scale data set.

To maintain the effectiveness and efficiency, we bring in the graph feature nodes to describe geometric features of data while capturing the association between different samples based on the broad structure. In this article, we utilize the symmetric normalized Laplacian $\mathcal{L}={\mathcal{D}}^{-\frac{\mathbf{1}}{\mathbf{2}}}(\mathcal{D}-\mathcal{W}){\mathcal{D}}^{-\frac{\mathbf{1}}{\mathbf{2}}}$ to prevent the gradient disappearance in the training process. Consequently, we construct a group of graph feature nodes to capture the inner data manifold within the broad framework to estimate the landslide susceptibility efficiently and precisely. Let $\mathbf{A}=[{\mathbf{Z}}^n,{\mathbf{H}}^m,\mathcal{L}]$; thus, the output-layer weight $\mathbf{W}$ is computed by Equation (3). The proposed model is called the graph-represented broad learning system (GBLS) whose brief construction is illustrated in Figure 4. So far, we have built the basic framework of GBLS that has the same flat single-layer structure as BLS but with the ability to obtain more spatially relevant information through graph nodes. In the following section, the incremental learning principle of GBLS is presented and demonstrated.

2.4.3. Incremental Learning with Additional Data

In landslide susceptibility mapping procedure, it costs massive time consumption to retrain the network from the beginning once the grid size is changed. Hence, the incremental learning capability of GBLS is demonstrated to reduce the unnecessary waste of sources and reveal the association between existing data and additional data.

Denote ${\mathbf{X}}_{old}$ as the initial input data whose mapped feature nodes ${\mathbf{Z}}_{old}^n$, enhancement nodes ${\mathbf{H}}_{old}^m$, and graph feature nodes ${\mathcal{L}}_{old}$ are calculated by Equations (1), (2) and (5), respectively. Then, ${\mathbf{A}}_{old}=[{\mathbf{Z}}_{old}^n,{\mathbf{H}}_{old}^m,{\mathcal{L}}_{old}]$. Given new input data ${\mathbf{X}}_{new}$, the new mapped feature nodes are generated by

$${\mathbf{Z}}_{new}^n=[{\mathbf{Z}}_{new1},{\mathbf{Z}}_{new2},\dots ,{\mathbf{Z}}_{newn}]=[{\varphi }_1({\mathbf{X}}_{new}{\mathbf{W}}_{e_1}+{\mathbf{\beta }}_{e_1}),{\varphi }_2({\mathbf{X}}_{new}{\mathbf{W}}_{e_2}+{\mathbf{\beta }}_{e_2}),\dots ,{\varphi }_n({\mathbf{X}}_{new}{\mathbf{W}}_{e_n}+{\mathbf{\beta }}_{e_n})]$$

(8)

the new enhancement nodes are widely expanded as

$${\mathbf{H}}_{new}^m=[{\mathbf{H}}_{new1},{\mathbf{H}}_{new2},\dots ,{\mathbf{H}}_{newm}]=[{\xi }_1({\mathbf{Z}}_{new}^n{\mathbf{W}}_{h_1}+{\mathbf{\beta }}_{h_1}),{\xi }_2({\mathbf{Z}}_{new}^n{\mathbf{W}}_{h_2}+{\mathbf{\beta }}_{h_2}),\dots ,{\xi }_m({\mathbf{Z}}_{new}^n{\mathbf{W}}_{h_m}+{\mathbf{\beta }}_{h_m})]$$

(9)

where the ${\mathbf{W}}_{e_i}$, ${\mathbf{W}}_{h_j}$ and ${\mathbf{\beta }}_{e_i}$, ${\mathbf{\beta }}_{h_j}$ are randomly generated in the initial preprocessing step. Then, newly added graph feature nodes ${\mathcal{L}}_{new}$ are calculated the same way as block Laplacian matrices in Equation (5), involving only new input data. Hence, the new combined feature matrix ${\mathbf{A}}_{new}=[{\mathbf{Z}}_{new}^n,{\mathbf{H}}_{new}^m,{\mathcal{L}}_{new}]$. In summary, we have the updating feature matrix

$${\mathbf{A}}_{update}=\left[\begin{array}{c}{\mathbf{A}}_{old}\\ {\mathbf{A}}_{new}^T\end{array}\right]$$

(10)

The corresponding pseudoinverse updating algorithm is deduced by

$${\mathbf{A}}_{update}^{+}=[{\mathbf{A}}_{old}^{+}-{\mathbf{PQ}}^T|\mathbf{P}]$$

(11)

in which ${\mathbf{Q}}^T={\mathbf{A}}_{new}^T{\mathbf{A}}_{old}^{+}$, and

$${\mathbf{P}}^T=\left\{\begin{array}{cc}{\mathbf{G}}^{+}\hfill & \hfill \mathrm{if}\mathbf{G}\ne 0\\ {(\mathbf{I}+{\mathbf{Q}}^T\mathbf{Q})}^{-1}{\mathbf{A}}_{old}^{+}\mathbf{Q}\hfill & \hfill \mathrm{if}\mathbf{G}=0\end{array}\right.$$

and $\mathbf{G}={\mathbf{A}}_{new}^T-{\mathbf{Q}}^T{\mathbf{A}}_{old}$.

Accordingly, the updated weights are inferred by

$${\mathbf{W}}_{update}={\mathbf{W}}_{old}+({\mathbf{Y}}_{new}^T-{\mathbf{A}}_{new}^T{\mathbf{W}}_{old})\mathbf{P}$$

(12)

where ${\mathbf{Y}}_{new}$ is the respective labels of new input data ${\mathbf{X}}_{new}$. In this manner, we only need to calculate partial required pseudoinverse rather than starting the the algorithm over with all samples. In general, the proposed model is constructed within the broad structure, exploiting the geometrical features of data by the graph-represented nodes to gain spatial correlation between regions. The specific framework of GBLS is shown in Figure 5. The performance of GBLS will be testified empirically in the next section.

3. Result

In this section, we compare the performance of the proposed GBLS with another seven popular machine learning methods, including AB, BLS, ELM, LR, RF, SVM, and SGD; then, we sort out the best method to evaluate the landslide susceptibility. Note that because we have already figured out several specific feature factors, the deep learning methods, such as CNN, are unable to demonstrate their superiority in feature extraction for image data, which will lead to low classification accuracy. Therefore, only the fully connected network has been utilized for the performance test.

We brief the experimental settings and then discuss the detailed results in the following.

3.1. Hyperparameter Setting

The hyperparameter determines in which hypothesis space the model is learning; hence, hyperparameter tuning is actually selecting the hypothesis space. For fair comparison, the hyperparameters in each methodology are determined by cross validation to gain their best performance on the testing sets while taking the complexity into consideration. Specifically, the number of ensemble learning circles in AB is picked from $\{5,10,15\}$. In SGD, the learning rate is set to 1. In RF, the number of trees in the forest is selected from $\{10,20,\dots ,100\}$; the max depth of tree is in $\{2,4,\dots ,10\}$. In LR, the inverse of regularization strength $C\in \{{10}^{-2},{10}^{-1},\dots ,{10}^2\}$. In SVM, Gaussian kernel is adopted with regularization parameter $C\in \{{10}^{-2},{10}^{-1},\dots ,{10}^2\}$ and $\gamma \in \{{10}^{-2},{10}^{-1},\dots ,{10}^2\}$. In ELM, the size of hidden nodes is searched for in the range of $\{1000,2000,\dots ,10,000\}$. In BLS and GBLS, the number of mapped features, mapping groups, and enhancement nodes is determined on $[1,2,\dots ,20]\times [1,2,\dots ,20]\times [1000,\dots ,10,000]$. The other hyperparameters use the default setting. Once the optimal hyperparameters are settled, the models run on this hyperparameter setting with their performance recorded. All experiments are conducted on the Matlab R2018a software platform, running on Intel-i7 2.4 GHz CPU with 16GB RAM PC.

3.2. Evaluation Indicators

In order to evaluate the performance of the model more comprehensively, we choose three commonly used indicators in the classification and briefly explain their role.

(1)

Test Accuracy

Test accuracy is defined as the percentage of correct predictions for the classification in the test data set, which is an intuitive indicator to assess the classifiers and which is easy to calculate. However, it is easy to obtain a high accuracy score by simply classifying all observations as the majority class when there is highly imbalanced data.

(2)

AUC-ROC

Receiver Operating Characteristics (ROC) is one of the principal evaluation metrics for evaluating a classification model’s performance. It shows the probability distribution where the positive case outranks the negative case by the classifier. Here, rank is decided according to the predicted order. Area Under the Curve (AUC) represents the degree or measure of separability. It tells how much the model is capable of distinguishing between classes. From an interpretation standpoint, the AUC-ROC score shows how good at ranking predictions the model is, and it is not sensitive to whether the sample categories are balanced. We use this visible indicator to compare the performance of multiple classifiers.

(3)

Sensitivity and Specificity

Sensitivity is the number of items correctly identified as positive out of the total true positives, whereas Specificity is the number of items correctly identified as negative out of the total negatives. Sensitivity is also termed as Recall. We always aim to have both sensitivity and specificity as high as possible. For the landslide susceptibility, if the sensitivity is high, that means the model is good at predicting the true landslide correctly. On the other hand, if the sensitivity is low, then the model will mark the landslide as a non-landslide area because it implies that the true positive rate is low. Furthermore, if the specificity is high, that means the true negative rate is high and as a result, the model is good at identifying non-landslide area correctly. Conversely, if the specificity is low, then the model falsely classifies the non-landslide as a landslide hazard area. Therefore, it is more important to have a model with high sensitivity for the test data.

3.3. Performance Comparison on Imbalanced Data

We aim to verify the predictive ability of the models to find out the best-fitting one for further research. Furthermore, the landslide susceptibility of the study area is required to be estimated as the model output. After grid division, there are a total of 129,883 samples in the study area in which the positive ratio is 27.36%. In this section, we compare the incremental learning efficiency along with the performance against several state-of-the-art machine learning methodologies listed above on imbalanced data. The landslide data extracted from the study area with 21 feature factors is roughly divided into five regions according to the river catchments, which might influence the geographical condition. Detailed information is listed in

Table 2. Description of Regional Data Sets.

Region	Samples	Positive Ratio	Catchment Range
I	25,008	28.25%	No.0∼No.9
II	22,897	29.25%	No.10∼No.19
III	22,085	26.40%	No.20∼No.29
IV	29,227	35.46%	No.30∼No.39
V	30,666	18.19%	No.40∼No.49

with a graphical illustration shown in Figure 6. Note that the distribution of landslides and non-landslides in those five regions is uneven, which matches the real hazard scenario. Then, we simulate the incremental data input by adding regions to the training process, that is, firstly the data in Region I is considered as the training data whereas the rest is used for testing, then the data in Region II is converted from the test set to the training set. This continues until only the data from Region V remains in the test set. During this process, the quantity of landslide and non-landslide data is kept asymmetrical. The performance and efficiency of selected models on those imbalanced data sets are testified with test accuracy and training time recorded in

Table 3. Test accuracy (%) and training time (s) of models on imbalanced data sets. In the first experiment, the training set contains 25,008 data points from Region I whereas the remainder is considered as test data. Similarly, Experiment 2 contains 47,905 training data points from Region I and II, Experiment 3 contains 69,990 training data points from Region I to III, and Experiment 4 contains 99,217 training data points from Region I to IV.

Model	Experiment 1		Experiment 2		Experiment 3		Experiment 4
	TestAcc	TrainTime	TestAcc	TrainTime	TestAcc	TrainTime	TestAcc	TrainTime
GBLS	72.43	7.46	72.64	10.39	73.17	14.80	81.81	17.11
AB	63.14	1.54	60.24	3.17	64.63	4.74	75.16	6.81
BLS	72.38	7.20	67.39	8.55	69.90	11.33	75.61	14.25
ELM	60.07	411.39	60.83	548.93	52.12	720.33	65.65	942.05
LR	40.87	3.11	53.89	5.52	42.22	8.89	53.19	16.09
RF	67.38	0.18	61.66	0.38	64.91	0.60	72.05	0.86
SGD	37.74	0.01	66.78	0.03	69.60	0.05	69.67	0.09
SVM	71.86	153.91	70.50	911.60	69.15	2503.31	80.52	5880.90

. Note that of all the models, only BLS and GBLS possess and utilize the incremental learning properties.

GBLS is proved to achieve greater performance on the imbalanced data than other methods based on the experimental results in Table 3. Note that massive data may fail to offer useful knowledge as more regions are added for training, which reduces the accuracy in some models. For example, the test accuracy of RF unexpectedly decreases by 5.72% as the amount of training data is increased from 25,008 to 47,905. Our algorithm exploits the training data more effectively through the broadly expanded feature nodes with graph-represented nodes to make better usage of the geometric correlation on data structures. This is demonstrated by the fact that the classification performance of GBLS gradually becomes better as the training data increases. Moreover, the capability of incremental learning inherited from the broad framework ensures low computational time consumption of GBLS as possible. In conclusion, GBLS is able to balance the classification performance and operational efficiency very well on imbalanced data.

To verify the performance of GBLS on the whole, the Friedman test [56] is exploited to compare numerical differences among multiple models from a statistical perspective. The p-value of the Friedman test is 0.0013. Therefore, we can reject the null hypothesis that the performances of these models are similar at significance level $\alpha =0.05$. The Nemenyi test [57] is immediately taken after this conclusion as a nonparametric post-hoc test to specify the discrepancy of model ranks. Once the difference between the average ranks of the two models is greater than the critical distance, the hypothesis that the two models have the same computational capability will be rejected at the corresponding significance level. The critical distance (CD) is calculated by

$$\mathrm{CD}=q_{\alpha }\sqrt{\frac{m(m+1)}{6n}},$$

(13)

where $q_{\alpha }$ is the critical value of studentized range distribution [58] based on the values of $\alpha$, the number of models m and the number of data sets n. Considering total 4 data sets and 8 models here, the critical distance for $\alpha =0.05$ is 2.407, and the ranking result in Figure 7 confirms excellent performance for GBLS. Therefore, the competitive capability of GBLS has been proved statistically.

3.4. Landslide Susceptibility Estimation

In the previous section, GBLS proves its efficiency and effectiveness on imbalanced data with incremental learning property. In the following, we will conduct the landslide susceptibility estimation in Zhouqu County. We choose the split of 70% training and 30% testing in order to balance the accuracy and generalization of the model. The spatial random sampling tool of GIS is used in sampling to avoid the reduction in representativeness of results due to the spatial correlation of data. In addition, stratified sampling is considered to predict the landslide susceptibility more precisely on the imbalanced samples. The cross validation is avoidable given large-scale data.

The experimental results are listed in

Table 4. Experimental results for landslide susceptibility estimation.

Model	Accuracy	Sensitivity	Specificity
GBLS	93.66%	86.30%	96.43%
AB	79.98%	51.73%	90.62%
BLS	91.25%	81.08%	95.09%
ELM	87.45%	72.70%	93.01%
LR	76.19%	32.54%	92.62%
RF	92.55%	83.37%	96.01%
SGD	68.83%	35.66%	81.32%
SVM	91.27%	80.52%	95.32%

, including the accuracy, sensitivity, and specificity for testing data.

The ROC curves for all models are shown in Figure 8. On the scale of test accuracy, the three best performers are GBLS (93.66%), RF (92.55%), and SVM (91.27%). On the scale of test sensitivity, the top three models are GBLS (86.30%), RF (83.37%), and BLS (81.08%). With AUC score as the standard, the highest three are GBLS (0.982), RF (0.975), and BLS (0.958). In total, GBLS presents its outstanding performance with the highest accuracy, sensitivity, and AUC score compared to other classifiers. Hence, we utilize the result of GBLS to draw the LSM.

Define the probability of the output category of 0–0.2 as a very low landslide risk area, 0.2–0.5 as a low landslide risk area, 0.5–0.7 as a moderate risk area, 0.7–0.9 as a high risk area, and 0.9–1 as a very high risk or landslide area. The LSM is drawn by estimated landslide susceptibility in GBLS (Figure 9), which well reflects the distribution status of regional landslides intuitively. According to the LSM, landslides in the study area are highly concentrated along the tectonic line, and the strike of the landslide is basically the same as that of the regional tectonic line. These highly susceptibility areas predicted by the model deserve more attention for landslide hazard identification and emergency response plans. If necessary, field exploration can be carried out to determine whether there is a hidden danger of landslides at the relevant location.

4. Discussion

4.1. Sensitivity Analysis on Resolution of Grid Units

The grid resolution is one of the important influencing factors of LSM. Too low a resolution does not guarantee the required accuracy of LSM. For instance, some research shows that the performance of landslide susceptibility prediction decreases when the resolution of grid units grows from 10 m to 100 m [59]. On the other hand, too high a resolution will greatly increase the complexity of calculation, which will cause a lot of time wastage and influence the observation of trends and patterns in the alpine-canyon region where geological conditions change frequently. Moreover, too high accuracy would easily lead to an overfitting problem because the model is always trained on the past remote sensing data. In this article, the 50 m × 50 m grid is adopted for study. This resolution can reduce the noise and false positives in the susceptibility analysis and in the meantime keep the accuracy in the susceptibility mapping process [60].

4.2. Relative Importance of Feature Factors in GBLS

Based on the empirical comparison of eight popular machine learning classifiers, GBLS is ultimately utilized as the best predictive methodology for landslide susceptibility. Though the extracted factors are not excessive, it is important to understand how the model works and which features play a key role. At the same time, the interpretable model reveals the decision-making process more clearly. Therefore, the importance of 21 relevant feature factors in GBLS is evaluated in this section. Because the origin feature factors have been transformed to the feature nodes and graph-represented feature nodes in GBLS, the Recursive Feature Elimination (RFE) [61] is applied to rank the feature importance more appropriately.

RFE is basically a wrapper-based method starting from building a model on the entire set of features and computing an importance score for each factor. The least important factor is removed; then, the model is rebuilt with the remaining factors. This procedure is recursively repeated until the desired quantity of features is eventually obtained. After the whole iteration, feature factors are ranked based on their time of elimination. RFE is an applied time-consuming method to obtain top features and scikit-learn makes it possible to implement the whole ranking procedure. Even better, scikit-learn provides RFECV to process recursive feature elimination with built-in cross-validated selection of the best number of features. So, we run RFE on Python using the 5-fold cross-validation for more generalization and the ranking of importance is shown in Figure 10. The heat plot shows that the most influential feature factors to the performance of GBLS are rainfall, NDVI, aspect, and TRI.

(1)

Rainfall

In the analysis of factor importance, all the factors related to the rainfall have be marked out as the most crucial to the landslide prediction, which meets our expectation. Note that among six indicators of rainfall, R1y, R3d, and R1d top the ranking list, followed by R6h and R1h. R10m, the lowest ranking in importance, is also in sixth place, well ahead of more than two-thirds of the indicators. As the main inductive factor of landslides, rainfall changes the structure of the soil through the infiltration of water, which ultimately causes the instability of the slope. In particular, short-term heavy rainfall causes the saturation zone on the surface of the soil and the peak of rainfall infiltration. It is prone to shallow rainfall landslides.

(2)

NDVI

Vegetation plays an important role in regulating climate change as one of the main components of the terrestrial ecosystem. Vegetation coverage is an important index to measure the surface vegetation condition which affects soil erosion; this is of great significance for regional environmental changes and monitoring. The soil in the Zhouqu area is relatively loose, and the roots of vegetation can effectively fix the soil and reduce the frequency of landslides. This explains the head position of NDVI in the feature importance ranking.

(3)

Aspect

Zhouqu County is a typical alpine-canyon landform located at the edge of the West Qinling and Qinghai-Tibet Plateau. Due to the influence of the monsoon climate, the vertical difference of the climate is obvious. The difference in aspect leads to a discrepancy in solar radiation, temperature, evaporation, and moisture, which have a profound impact on vegetation, soil, and hydrology. Hence, the influence of aspect on slope stability and landslide is significant.

(4)

TRI

From the perspective of origin, the mechanical properties of slope rock and soil are the direct and fundamental cause of landslides. Slope is not the only measure of the terrain roughness because a large slope in a highly inclined terrain does not represent a more rugged terrain. TRI implies the variation of slope in a terrain which is indicated by curvature. Though we could find that slope ranks relatively high as well, TRI is more crucial to the landslide susceptibility estimation in GBLS.

Though we have sorted all feature factors by their influence on the model, it is worth mentioning that RFE is a greedy wrapper-based algorithm coupled with the machine learning methodology in use, i.e., GBLS, for finding the most important features of the model. High ranking demonstrates that the feature’s attribute to model is great; nevertheless, low ranking does not necessarily mean the feature factor has nothing to do with the landslide. Some factors, such as slope, have a controlling effect on the landslide and significantly affect the stress distribution of the rock and soil. Those factors also play an important role in landslides today. As landslides are caused by comprehensive topographic factors, we should not underestimate or even ignore the actual role of factors due to their low contribution to the model.

4.3. Comparison of Models

(1)

GBLS vs. BLS

Basically, GBLS can take better use of spatial features by adding graph-represented nodes within the broad framework, which improves accuracy while inheriting the incremental capability of BLS. For landslide data, there exists some peculiarity that fits the usage of the broad structure. First of all, most of the feature factors related to landslide hazards are geological and meteorological data, which are characterized by numerous sparse features in disparate degrees of collinearity. The more universal the study area, the larger the sample size, and the more obvious the representation of nonlinear characteristics. Statistical methods generally assume data independence in some specific distribution. However, landslide hazards have certain swarm due to the same triggering factors, such as earthquakes, heavy precipitation processes, etc., which would lead to potential interactions and massive correlations between features. If we directly apply the original data indiscriminately, the latent correlation between feature factors will be ignored and missed, and the feature information cannot be fully utilized, resulting in the decline of accuracy. The broad structure is capable of improving this dilemma by generating feature nodes and enhancement nodes from the original data. Furthermore, the volume of landslide data keeps growing over time. The incremental learning ability of the methodology is one of the factors that needs to be premeditated, not only to reduce the redundant cost or computation consumption but also the correlation between the current landslide and the future occurrence.

(2)

GBLS vs. Single-layer methods

In Section 3.3 and Section 3.4, we have compared the performance of GBLS to other seven models including boosting, statistical learning, and single-layer machine learning methods. The results show that GBLS has superiority over other methods in handling imbalanced data while restricting the time consumption in a low range through incremental learning algorithm. Although GBLS is not the lowest in terms of time consumption, it is still very efficient in handling increasing amounts of data and avoiding repetitive computations among algorithms with comparable performance, such as SVM.

(3)

GBLS vs. Deep learning methods

Compared to the deep structure, GBLS is very simple because there is no cascade of layers. Similarly, the broad network does not need to use gradient descent to update the weights because there is no multi-layer connection, so the computational speed is much faster than deep learning. Moreover, when the accuracy of the network does not meet the requirements, the performance can be improved by increasing the “width” of the network. GBLS is suitable for learning with few data features but with high requirements for real-time prediction.

4.4. Contribution and Limitation

The GBLS is proposed as a graph-represented method aiming to balance the effectiveness and efficiency of landslide susceptibility learning. It is a useful method of the classification problem, which performs well on the imbalanced data as shown in Section 3. Specifically, it should deal with the data in the alpine-canyon region where the ecological and geological environment is fragile. The remarkable tectonic activities might cause hazards frequently; thus, the regional data also varies from time to time. In this case, the incremental property of GBLS comes in handy, avoiding time consumption while maintaining model performance and mining patterns of variation through a comparison with past data.

In addition, LSM has a great potential to support and advise the local government on land-use planning, disaster mitigation, etc. GBLS with incremental learning capability and high performance can fit the ever-changing needs of emergency management. The LSM might implement regular updates before each monsoon season or emergency updates for disaster mitigation after a sudden earthquake. The high susceptibility areas predicted by GBLS deserve more attention for landslide hazard identification and emergency response plans. If necessary, field exploration can be carried out to determine whether there is a hidden danger of landslides at the relevant neighborhood.

Even though GBLS has many great competencies, it inevitably has some pressing problems that need to be solved.

In the algorithm, the construction of GBLS relies basically on the transferred features, including feature nodes, enhancement nodes, and graph nodes, but we do not currently know the clear relationship between the number of nodes and the prediction accuracy. In practice, we simply utilize grid search to find the optimal value of parameters in which the searching range is determined empirically without effective theoretical assistance. Accordingly, the nodes’ combination obtained by grid search can only be considered as locally optimal rather than globally optimal in theory. Due to this status quo, we are currently seeking the regular pattern of performance as the number of nodes changes in order to generalize a feasible algorithm.

Theoretically, we tend to constantly expand the searching range of nodes in order to achieve higher precision. However, it is not guaranteed that more nodes indicates higher accuracy, thus simply increasing the number of nodes does not necessarily improve model performance when the accuracy does not meet expectations. Furthermore, a large volume of nodes will enlarge the transferred data matrix in the first layer, resulting in heavy computation of pseudo inverse. If the spatial resolution of grid units in the study area is chosen to be high at the same time, the algorithm might run out of computational memory and cannot be computed due to the computation complexity. We are currently considering using the drop-out method to solve this problem, but it is difficult to avoid the loss of accuracy.

In addition, the selection of conditioning feature factors is adjustable. Some potential influencing factors, such as wind aspect, which affects local changes in precipitation, can also be taken into consideration as wind directions might cause higher local rainfall. GBLS can play an essential role in alpine-canyon regions due to its incremental learning capability, yet this property may not always be a benefit when dealing with other less frequently changing geological situations. Hence, our proposed algorithm is a general algorithm which is more advantageous and applicable in handling geographic data with more significant changes, such as alpine-canyon regional data.

Zhouqu County is one of the most iconic areas for the alpine-canyon accumulation landslide, and the characteristics of landslides in this area represent well other regions with similar geo-conditions. GBLS is supposed to be able to provide solutions under the climate change conditions, for example, incorporating the precipitation projection in the IPCC reports. In this study, we find that rainfall is the most important feature overall whereas rainfall data for each smoothing window belongs to a set of constantly changing time series data. By its incremental learning capability of GBLS, it is possible to update LSM according to different IPCC precipitation scenarios whenever new precipitation data are available. In the further research, we plan to explore the feasible combination with longitudinal data for multidimensional modeling analysis.

5. Conclusions

As a typical alpine-canyon region coupled with complex geological environments and significant tectonic activities, Zhouqu County is a rare, high incidence area of landslides, debris flows, and earthquakes in China, which brings a huge threat to the residents. The frequent occurrence of a large number of geological disasters in Zhouqu and the surrounding areas in recent years has attracted many scholars to conduct relevant studies on environmental issues. In this article, we focus on the landslide susceptibility assessment in the study area so as to take precautions in advance. The landslides that occurred around Zhouqu County have some characteristics in common, such as high frequency of occurrence with fast data update. If the landslide susceptibility needs to be recalculated every time from the first beginning, it will inevitably lead to redundant time waste, and it is impossible to observe the changing pattern of disasters in time. Moreover, the distribution of landslides is chaotic and without obvious rules, and there exist complex correlations between influencing factors, which will reduce the prediction accuracy and fail to obtain a good explanation if the internal rules of features are not deeply explored. To solve these problems, we propose an algorithm named graph-represented broad learning system (GBLS) to explore the manifold structure and geometric correlation of landslide features. The framework of the standard broad learning system is utilized to maintain its efficiency of classification while graph feature nodes are generalized simultaneously to extract the manifold attributes. In this way, GBLS can balance the estimating performance as well as the efficiency in processing incremental data.

Then, we compare the performance of GBLS along with seven other state-of-the-art methods, in which GBLS shows its excellent capability both in efficiency and effectiveness in landslide susceptibility estimation on either imbalanced data or stratified samples. By sorting the importance of factors, the most important factors affecting landslide susceptibility in Zhouqu County and its surrounding area are rainfall, NDVI, aspect, and Terrain Ruggedness Index. The result is consistent with the actual situation because the 2010 mudslide was caused by torrential rain. The LSM could help the government to better carry out the risk prevention and monitor landslide dynamics in the alpine-canyon region with more validity.

As a new algorithm in LSM and geology, GBLS is also an incremental learning algorithm that is capable of updating the training process without recalculating the whole model given new data. Note that the performance of GBLS is influenced by the feature node combination which is determined by grid search in practice, and we expect to find a more rigorous theoretical solution to figure out the correlation between nodes and accuracy. On the other hand, the selection of feature factors greatly affects the manifold representation in the graph nodes because the distance between regional data is calculated based on the features. Inappropriate features will cause incorrect manifold correlation while reducing the accuracy of the model. In further study, we are planning to select more landslide-related features to more deeply explore landslide influencing factors for better disaster prevention and control. Moreover, we will take advantage of incremental property to establish the dynamic landslide susceptibility mapping model containing time series data, which is beneficial to fuse transverse data features with longitudinal time-varying features to obtain better research outcomes.