Failure Mechanism and Long Short-Term Memory Neural Network Model for Landslide Risk Prediction

Abstract

Rockslides along a stepped failure surface have characteristics of stepped deformation characteristic and it is difficult to predict the failure time. In this study, the deformation characteristics and disaster prediction model of the Fengning granite rockslide were analyzed based on field surveys and monitoring data. To evaluate the stability, the shear strength parameters of the sliding surface were determined based on the back-propagation neural network and three-dimensional discrete element numerical method. Through the correlation analysis of deformation monitoring results with rainfall and blasting, it is shown that the landslide was triggered by excavation, rainfall, and blasting vibrations. The landslide displacement prediction model was established by using long short-term memory neural network (LSTM) based on the monitoring data, and the prediction results are compared with those using the BP model, SVM model and ARMA model. Results show that the LSTM model has strong advantages and good reliability for the stepped landslide deformation with short-term influence, and the predicted LSTM values were very consistent with the measured values, with a correlation coefficient of 0.977. Combined with the distribution characteristics of joints, the damage influence scope of the landslide was simulated by three-dimensional discrete element, which provides decision-making basis for disaster warning after slope instability. The method proposed in this paper can provide references for early warning and treatment of geological disasters.

1. Introduction

Excavation often causes slope damage in areas with complex geological conditions. This results in huge casualties, property losses, and construction delays [1]. Step-path rockslides are a common type of rock slope failure in engineering. The failure surface is often formed by a combination of steeply and gently inclined fissures. Therefore, the combination of steeply and gently inclined fissures has an important influence on the deformation, failure evolution, and stability of the slope [2]. In sedimentary rocks, stepped translational slides are mainly developed along bedding planes, whereas in the granite and metamorphic rocks, many relaxation joints are commonly developed as a result of stress release as the material is removed from the slope. Relaxation joints often have a gentle dip angle, consistent with the slope direction. Combined with steep dip fractures, bed-like rock mass structures readily form resulting in a landslide failure surface with a combination of steep and gentle fractures. This stepped failure mode is more common in western China and is often combined with steep structural fissures that form a step-path failure, such as the landslide of the Manwan hydropower station [3], the Songta hydropower station, and Maji hydropower station [4,5,6]. Generally, the failure of a stepped slope is characterized by multiple sliding zones, multistage deformation, long creep time, and displacement along with fractures [7]. It is difficult to predict the failure time because the evolution of the failure surface is uncertain [8]. Moreover, the deformation prediction model and failure mechanisms of step-path landslides, which are key factors in slope stability evaluation and reinforcement design [3], are still poorly understood.

The typical landslide caused by excavation is mainly due to excavation making the slope lose its anti-sliding force and blasting vibration. Expressway and hydropower projects are often associated with stepped failure along the relaxation joints, but there are few studies on this type of landslide [9]. To date, most of the existing research has been focused on the analysis of the mechanisms of combination landslides with steeply and gently inclined fissures [6,10,11]. Studies have also been focused on the stepped failure mode of different fracture combination forms by considering the joint persistence ratio (JPR) of fractures and the dip angle of the rock bridge [12,13]. Numerical methods have been used to simulate the development of stepped failures on rock slopes where the fracture surface was formed by brittle fracture propagation within the rock [12]. However, there is little research on the prediction of step-path landslide trends caused by excavation.

Landslides caused by rainfall have periodicity, but the displacement caused by the random vibration load of construction blasting is difficult to predict. Many researchers have discussed in detail landslides caused by rainfall or reservoir water level changes, and believe that the main factors of landslides are rainfall and groundwater level changes [14,15,16]. Numerical methods are generally used to geomorphologically assess the stability of slopes and evaluate risks associated with landslide hazards. Most are focused on soil or loose layer landslides [17,18,19]. Few studies on deep-seated landslides induced by rainfall are available [20].

Mechanical parameters are important factors in slope stability evaluation [21,22]. The commonly used methods in determining rock mass parameters are in-situ tests, laboratory tests, engineering analogy and numerical inversion. However, in-situ and laboratory tests do have limitations and uncertainties [23]. With the rapid development of computer technology, the displacement back analysis method based on field monitoring information has been applied in these years [24,25,26]. Combined with artificial intelligence and numerical analysis method, it is successfully applied to the estimation of slope rock mass strength parameters [24].

Many landslide prediction models require many expensive geotechnical characterizations of the materials involved in landslides, and therefore may be applied only in limited cases [27,28,29]. With the development of pattern recognition and artificial intelligence technologies, computational intelligence has been used in landslide deformation prediction. Intelligence models such as artificial neural networks (ANNs) [28,29,30,31] support vector machine (SVM) [9], extreme learning machine (ELM) models [32], Harmony Search (HS) algorithm [33] and random forests have been used for landslide prediction in recent years [34,35]. However, these traditional machine learning models have many drawbacks. For example, the inherent and deep features of input variables cannot be extracted [36], and they cannot adequately describe the internal processes of a landslide in motion [28]. Most of these models are only applicable to soil landslides with long-term creep characteristics, factors with periodic influence, and continuous sliding surfaces. The stepped sliding surface rock landslides prediction model has not been reported. Predictive recurrent neural networks (RNNs) can correct errors using backpropagation and gradient descent algorithms, but they do not solve the problem of gradient disappearance.

Deep learning does require model parameters to be set because of its autonomous processes, and it can effectively avoid local optima. Deep learning algorithms have received little attention for LSP [37]. Huang et al. proposed a novel deep learning algorithm, namely a fully connected sparse autoencoder (FC-SAE) LSP, which is suitable to predict potential landslide disaster points in the spatial prediction of landslide susceptibility areas [37]. There is little research on a prediction model for rock landslides coupled with excavation disturbances and rainfall. So, it is necessary to develop a novel machine learning method in landslide displacement prediction. Compared to the aforementioned artificial intelligence landslide prediction model, LSTM (long short-term memory) can effectively solve the problem of gradient disappearance or gradient explosion in a general recurrent neural network (RNN). It is a type of time recurrent neural network. Since it has the function of self-updating weight and significantly improves the accuracy of the network, LSTM can be used as a complex nonlinear element to construct larger deep neural networks [38]. The LSTM model does not need to separate the trend and periodic terms in the original deformation information and can compensate for the deformation trend prediction caused by the accidental interruption of monitoring data [39]. It is an important and effective technique for forecasting stepped rock landslide deformation.

Because the landslide treatment time and reinforcement method depend on the landslide deformation development trend, it is very important to predict the future deformation trend and damage influence area of landslides. Compared with other artificial intelligence methods, the LSTM model proposed in this paper can well reflect the shear process of sliding surface and the influence of rainfall, blasting and other short-time on the deformation trend, and greatly improve the accuracy of deformation prediction. Using displacement inversion analysis method coupling with LSTM of deformation monitoring data and discrete element numerical method, the joints persistence along sliding surface can be better determined, so as to further determine the shear strength of sliding surface. The novelty of this study is to provide an effective method of landslide disaster prediction.

2. Materials and Methods

2.1. Geological Setting

The Fengning Landslide is located at the excavation slope of the Fengning Pumped Storage Power Station, Fengning Manchu Autonomous County, Chengde City, Hebei Province (Figure 1). Rapid sliding first occurred on 6 March 2018. The boundary of the landslide mass formed an armchair shape, and the front edge formed a shear outlet of the excavated slope at an elevation of 1230 m. The sliding resulted in a stagger of 2.0 m at the rear edge. By 11 September 2019, the landslide shear outlet had moved down the excavated slope to an elevation of 1185 m. After October 2020, further accelerated deformation was noted following the rainy season. The landslide remains under constant deformation, and has not been completely stabilized. This seriously impacts the safety of construction workers, material yards, and other building equipment at the foot of the landslide. Moreover, it poses a direct threat to the safety of adjacent areas, and poses great difficulties for reinforcement and treatment (Figure 1c).

Field geological survey was carried out from 2018, and 3 drill holes were arranged to explore the landslide feature. The bedrock is a fine-grained biotite adamellite (T_h). The rock outcrops on the slopes were mostly weakly weathered granite. The vertical profile of a typical section of the slope is shown in Figure 2. Below the elevation of 1230 m is the excavation slope of the rockfill material, and the excavation slope ratio is 1:0.7. There are three groups of well-developed joints: (1) tectonic joints J1, described as having a northwest strike of 335–355° and a dip of 75–90° toward the northeast (NW335–355°/75–90° NE) and an average joint spacing of 35 cm; (2) tectonic joints J2, NE50–70°/60–80° SE and a general spacing of 30–80 cm; and (3) relaxation joints J3, NE40–60°/35–45° SE and a general joint spacing of 50–80 cm. Most of the joints were flat, rough, and slightly open. The drilling data show that the relaxation joints below 1185 m of elevation were not developed, and the spacing was only 5 m. By November 2017, the excavation elevation was 1185 m. There are three steeped sliding surface. Many tensile cracks occurred at the top area between YII1 and YII2, the sliding surface with buried depth of 10~20 m is composed of the gently inclined relaxation joint and the steeply inclined tectonic fissure with stepped features. Slope excavation engineering directly induced the slope failure.

The groundwater is mainly composed of granite fissure water, which occurs in bedrock and is recharged by atmospheric precipitation and discharged to the overburden and the surface. The depth of the groundwater table on the slope was 70–100 m. The landslide was not affected by groundwater table. The borehole core revealed that the relaxation zone below the lower elevation of 1185 m of the landslide mass was excavated. Therefore, the shear outlet of the lowest sliding surface of the Fengning stepped landslide was not lower than the elevation of 1185 m.

2.2. On-Site Monitoring

This slope was relatively stable before excavation of the rockfill material. Slope instability was triggered by excavation and rainfall. The sliding surface developed along the rock-joint system to form a sliding mass. Therefore, to ensure the safety of the dam rockfill material yard at the front of the landslide, surface displacement and borehole inclinometer monitoring systems were designed and installed in the landslide’s front, middle, and trailing edge sections (Figure 3). A Leica TCA2003 automatic total station, borehole fixed inclinometer, rain gauge, ground fissure meter, and unmanned aerial vehicle (UAV) were used to monitor the three-dimensional surface deformation and deep displacement. In 12 April 2018 surface displacement monitoring points were arranged on the landslide. Two monitoring stations (YII5-1 and YII5-2) were located on the excavated slope, seven monitoring stations (YII2, YII3-2–YII3-4, and YII4-2–YII4-4) were located at the landslide surface, and five monitoring stations were located outside the landslide boundary (YII1, YII3-1, YII3-5, YII4-1, and YII4-5). The locations of the monitoring stations are shown in Figure 3. Among them, the positive direction of X points to the north, the positive direction of Y points to the east, and the positive direction of H is vertical subsidence. Monitoring began in early April, 2018. As of 20 June 2021, 352 periods had been observed. Fixed borehole inclinometers were installed in the three boreholes. In August 2018, three 25-m inclinometer holes (CX01, CX02, and CX03) were installed at the front edge of the landslide excavated slope at elevations of 1215, 1215, and 1200 m, respectively.

2.3. LSTM Algorithms and Displacement Forecasting Model

2.3.1. LSTM Algorithms

Long short-term memory (LSTM) belongs to a time-recursive recurrent neural network, which can effectively solve the problem of gradient disappearance or gradient explosion in general recurrent neural networks (RNNs). RNNs have good learning ability for information within short distances [40]. Therefore, it is difficult to capture long-term time-related information. To solve this problem, Hochreiter and Schmidhuber proposed a long-term and short-term memory network in 1997 [41]. LSTM is a new type of deep machine-learning neural network. Based on RNN, this model adds a memory cell, input gate, forget gate, and output gate to control the transmission of the information at separate times. This can improve the ability of the RNN to process long sequence data. After building the LSTM, the gradient disappearance problem can be solved, and the network accuracy can be improved significantly. In this paper, the discrete element method was used to solve the joint persistence ratio of sliding surface under the action of excavation, rainfall and blasting through displacement back analysis. Then, the LSTM model was used to predict the joint persistence ratio of sliding surfaces under different displacements stages. The results can evaluate the landslide influence range.

The LSTM memory unit has a forgetting gate ($f_t$), input gate ($i_t$), and output gate $o_t$. The most important of these gates is the forget gate [42]. The main inputs and outputs of the LSTM structure are illustrated in Figure 4.

The step-by-step LSTM algorithm is described as follows [43,44].

Firstly, to decide what information we’re going to throw away from the cell state. The decision is made by a sigmoid layer named as forgetting gate layer (Equation (1)). It looks at $h_{t-1}$ and $x_t$, and outputs a number between 0 and 1 for each number in the cell state $C_{t-1}$. If the output of the function sigmoid is close to 0, then the previously stored information $C_{t-1}$ will be “forgotten”. If the output of the function sigmoid is close to 1, the previously stored information $C_{t-1}$ will be completely accepted.

$$f_t=\sigma \left(W_f·\left[h_{t-1},x_t\right]+b_f\right)$$

(1)

where $x_t$ is the input series, $h_{t-1}$ is the output of the previous memory block, $W_f$ is the weight vector, $b_f$ is the bias vector, and $\mathsf{\sigma }$ is the sigmoid function.

Step 2: To decide what new information we’re going to store in the cell state. Firstly, to decide which values we’ll update by using a sigmoid layer named the “input gate layer” (Equation (2)). Then, a $\tan h$ layer creates a vector of new candidate values $\widetilde{C_t}$, that could be added to the state (Equation (3)).

$$i_t=\sigma (W_i·\left[h_{t-1},x_t\right]+b_i$$

(2)

$$\widetilde{C_t}=\tan h(W_c·\left[h_{t-1},x_t\right]+b_{c)}$$

(3)

Step 3: To create an update state of combining with Equations (2) and (3). Add the gender of the new subject to the cell state, to replace the old one we’re forgetting, update the old cell state $C_{t-1}$ into the new cell state $C_t$.

Step 4: Multiply the old state by $f_t$, drop the information which decided to forget earlier. Then add the new information $i_t\times C_{t-1}$ as decided in the previous steps.

$$C_t=f_t\times C_{t-1}+i_t\times \widetilde{C_t}$$

(4)

Finally, it is necessary to decide what we’re going to output. This output will be based on cell state, but will be a filtered version. Firstly, run a sigmoid layer which decides what parts of the cell state we’re going to output (Equation (5)). Then, put the cell state through $\tan h$ (to push the values to be between −1 and 1) and multiply it by the output of the sigmoid gate, it only output the parts we decided to (Equation (6)).

$$O_t=\sigma \left(W_0·\left[h_{t-1},x_t\right]+b_o\right)$$

(5)

$$h_t=\tan h\left(C_t\right)\times O_t,$$

(6)

where $O_t$ is the output status of the input gate, $h_t$ is the output status of the output gate and is the output of the LSTM core unit.

The approach used to establish the forecasting model in this study is shown in Figure 5.

2.3.2. Data Preprocessing

The component response mode of the cumulative displacement S(t) of a landslide is usually represented by the following model:

$$S\left(t\right)=A\left(t\right)+B\left(t\right)+C\left(t\right),$$

(7)

where A(t) is the trend component of landslide displacement, which is controlled by its own geological condition, rock mass deformation capacity, potential energy, and evolves with time; B(t) is the component of the fluctuation term, which is usually caused by seasonal rainfall; and C(t) is a random disturbance component, which is mainly affected by blasting, earthquakes, temperature, humidity, air pressure, and other meteorological information. In this study, the displacement data were extracted from the raw observation data using the moving average method (MAM). In the LSTM model, it is not necessary to separate trend items and periodic items from the monitoring information.

2.3.3. Loss Function

The loss function is used to evaluate the deviation between the value predicted by the machine learning model and the observed value. The mean square error (MSE) loss function was used in this study.

$$MSE=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{\left(x_i-{\widehat{x}}_i\right)}^2,$$

(8)

where N is the number of samples, ${\widehat{x}}_i$ and $x_i$ are the predicted and real values, respectively.

2.3.4. Machine Learning Model

The gradient descent method is the most used method to train machine learning models. It can be divided into batch gradient descent, random gradient descent, and small batch gradient descent. In this study, an extended algorithm of the random gradient descent method was adopted. It is a very popular algorithm in the deep learning field and is the adaptive moment estimation algorithm (Adam), which adaptively adjusts the learning rate for each parameter. The advantage of the Adam algorithm is that it calculates the exponential weighted average of gradient square $g_t^2$, and it calculates the exponential weighted average of gradient $g_t$ [43].

$$M_t={\beta }_1{M}_{t=1}+\left(1-{\beta }_1\right)g_t$$

(9)

$$G_t={\beta }_2{G}_{t=1}+\left(1-{\beta }_2\right)g_t☉g_t,$$

(10)

where β₁ and β₂ are the decay rates of two moving averages; usually given as β₁ = 0.95 and β₂ = 0.99. When β₁ and β₂ were close to 1.00, they were corrected according to the following formulas:

$$\widehat{M_t}=\frac{M_t}{1-{\beta }_1^t}$$

(11)

$$\widehat{G_t}=\frac{G_t}{1-{\beta }_2^t}$$

(12)

Therefore, the Adam parameter update difference is:

$$∆{\theta }_t=-\frac{\alpha }{\sqrt{\widehat{G_t}}}\widehat{M_t}$$

(13)

The learning rate, α, is usually set as 0.001.

2.3.5. Model Forecast Accuracy Evaluation

The back-propagation algorithm is used in the LSTM model parameter training. In this study, the root mean square error (RMSE), correlation coefficient (R), mean absolute error (MAE), and mean absolute percentage error (MAPE) were selected for the evaluation index of the LSTM effect. The RMSE and R can be calculated using Equations (14) and (15).

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left(x_i-\widehat{x_i}\right)}^2}{N}}$$

(14)

$$R=\frac{{{\displaystyle \sum }}_{i=1}^N\left(x_i-\overline{x}\right)\left(\widehat{x_i}-\overline{\widehat{x}}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^N{\left(x_i-\overline{x}\right)}^2}\sqrt{{{\displaystyle \sum }}_{i=1}^N{\left(\widehat{x_i}-\overline{\widehat{x_i}}\right)}^2}},$$

(15)

where $x_i$ and $\widehat{x_i}$ are the monitoring and prediction values, respectively, $\overline{x}$ and $\overline{\widehat{x}}$ are the average values of the monitoring value and prediction value, respectively, and $N$ is the number of samples.

2.4. Discrete Element Numerical Method

The discrete element numerical method was used to determine the sliding surface parameters and the dynamical action. In the parameter inversion of displacement back analysis, we can calculate the displacement corresponding to different shear strength parameters in numerical analysis [25,45]. This method provided sufficient samples for the inversion of the shear strength parameter. But in the previous work, since two indexes of shear strength (c, φ) were obtained, the inversion accuracy is limited. In this paper, only one index of joint persistence ratio on the sliding surface is obtained by displacement back analysis to determine the strength parameters, so the calculation results are more reasonable.

The shear strength of fractured rock mass is mainly determined by the joint persistence ratio (JPR). The universal function between slope displacement and JPR is established in numerical simulation analysis. The inversion result is unique and can make the accuracy higher. The shear strength parameters of stepped sliding surface can be expressed as following formulas.

$$S=f\left(k\right)$$

(16)

$$c=\left(1-k\right)c_r+k{c}_j$$

(17)

$$\phi =\arctan \left[\left(1-k\right)\tan {\phi }_r+k\tan {\phi }_j\right]$$

(18)

where $S$ is displacement of slope (m), $c_r$ is the cohesion of rock bridge (MPa), $c_j$ is the cohesion of joints (MPa), ${\phi }_r$ is the internal friction angle of rock bridge (°); ${\phi }_j$ is the internal friction angle of joints (°). $k$ is the joint persistence ratio (JPR).

For the granite slope cut by structural plane, the discrete element method [46] can effectively analyze the slope stability under blasting dynamic action, and is suitable for simulating the overall dynamic failure process of rock mass. Shi Chong took the steep rock slope as an example, established the seismic wave propagation model by discrete element method, and then analyzed the stability characteristics of steep rock slope under seismic dynamic action [47]. In the discrete element numerical analysis in this paper, we suppose the normal force and tangential force exists between elements [48,49]. The normal force and tangential force are expressed by normal spring force (F_n) and tangential spring force (F_s) respectively [50].The damping force is changed by adjusting the damping coefficient. Damping force (F_v) is directly proportional to the unit velocity ($x^{\prime }$) and the direction is opposite to the unit motion direction:

$$F_{\mathrm{v}}=-\eta x^{\prime }$$

(19)

where η is the damping coefficient.

For the application of blasting load, set the sinusoidal loading wave through the following formula and input it into the numerical model.

$$V=A\sin (2\pi t/T)$$

(20)

where V is the amplitude value corresponding to different t times. A is the loading amplitude and T is the sine wave period. According to the common blasting parameters, the incident wave frequency is 0.25 MHz and the loading displacement amplitude is 1 × 10⁻⁷ m. The dynamic simulation of slope stability by blasting is carried out with damping coefficient of 0.1 as the parameter value.

On considering the influence of rainfall on slope stability, the seepage pressure is equivalent to node force.

3. Results

3.1. Displacement Monitoring Results

3.1.1. Surface Displacement

By 20 June 2021, a total of 352 periods of displacement had been monitored. The displacement of seven surface monitoring stations within the landslide mass changed significantly, especially in the main sliding direction (close to the Y direction of the coordinate). From 30 April 2018 to 20 June 2021, the deformation in the main sliding area is much greater than the monitoring points outside the boundary of the landslide mass. The maximum cumulative displacement was 580.6 mm for monitoring station YII5-2. The downward displacement trend of the landslide mass along the vertical contour line is notable and continues (Figure 6a). The cumulative displacement of the five monitoring stations outside the boundary of the landslide mass did not change significantly (Figure 6b). Since 6 March 2018, when the first sliding occurred, the monitoring stations were displaced approximately 130–250 mm within the first four months. The deformation rate was fast and represented the first accelerated deformation stage. Then, the slope entered a slow deformation stage, with a displacement over the next seven months of only 50–100 mm. The deformation of the landslide mass increased again in October 2018, when it entered the second accelerated deformation stage. It entered the third rapid deformation stage in September 2019 and in October 2020, entered the fourth accelerated deformation stage, which is still developing. All monitoring points in the landslide area show that the deformation trend has the stage of acceleration-creep–acceleration-creep, and up to now, the deformation is still increasing.

3.1.2. Deep Seated Displacement

Currently, the effective monitoring stations for landslide inclinometer monitoring are fixed inclinometers CX03 and CX01 at elevations of 1200 and 1215 m, respectively. The deformation process curve of CX03 from 12 August 2018, to 4 October 2019, is shown in Figure 7a. The deformation process curve of CX01 from 3 October 2018, to 12 September 2019, is shown in Figure 7b. The inclinometer observation data shows that there is a shear displacement at the front edge of the landslide mass at elevations of 1185–1190 m. The displacement of monitoring point CX03 (1200 m elevation) reached 1341 mm on 4 October 2019. The shear displacement of monitoring point CX01 (1215 m elevation) reached 972 mm on the same date.

3.1.3. Displacement with Rainfall and Blasting

According to the displacement characteristic curves of fixed inclinometer monitoring points CX01 and CX03, the displacement can be roughly divided into four stages. The first stage was a uniformly accelerated deformation stage. From 12 August to 20 October 2018, the displacement reached 363.8 mm and the deformation rate was 5.35 mm/d. The second stage was a slow deformation stage. From 21 October 2018, to 7 June 2019, the displacement reached 933.1 mm and the deformation rate was 2.73 mm/d. The third stage was an accelerated deformation stage. From 8 June to 11 August 2019, the displacement reached 1340.8 mm and the deformation rate was 6.39 mm/d. The fourth stage was a slow deformation stage. From 12 August to 4 October 2019, the displacement reached 1343.7 mm. The change trend for these four stages of deep displacement was consistent with the deformation trend of the landslide surface.

It can be seen from Figure 8 that after the rapid moving on 6 March 2018, the deformation of landslide mass can be divided into two stages: accelerated deformation and creep deformation. Except for monitoring station YII1, all accelerated deformation stages are closely related to the rainfall season, which occur between July and October every year, with a certain time lag. Instead of the rainy season, the landslide mass enters the creep deformation stage. A comparative analysis of rainfall and landslide displacement in the rainy seasons in 2018 and 2019 shows a significant correlation between rainfall and landslide displacement. The displacement speed is closely related to the rainy season (Figure 8). In the first rainy season, from 15 June to 14 September 2018, with a cumulative rainfall of 347.0 mm, the deformation increased the most. Deformation at station YII5-1 increased from 46.4 mm to 196.8 mm, and the deformation speed was 1.671 mm/d. At monitoring point YII3-3 in the middle of the landslide mass, the deformation increased from 26.9 mm to 120.2 mm. The deformation speed was 1.037 mm/d. The deformation of YII2 at the upper part of the landslide increased from 29.9 mm to 135.3 mm, at a rate of 1.171 mm/d. From 15 June to 14 September 2019, in the second rainy season, the cumulative rainfall was 254.0 mm, and the deformation rate was smaller than that in the first rainy season (

Table 1. Deformation with rainfall.

Rainy Season	Date	Monitoring Station					Cumulative Rainfall, mm
		YII1	YII2	YII3-3	YII4-3	YII5-1
The first	16 June 2018	2.2	29.9	26.9	44.9	46.4	0.5
	14 September 2018	12.4	135.3	120.2	186.9	196.8	347.0
	Accelerated displacement (mm/d)	0.113	1.171	1.037	1.578	1.671
The second	13 June 2019	33.8	216.6	199	303.4	323.9	0.5
	15 September 2019	52.9	267.5	251.1	372.1	405.2	254.0
	Accelerated displacement (mm/d)	0.212	0.566	0.579	0.763	0.903

). The spatial characteristics of deformation show that YII5-1, located at the shear outlet at the slope toe, was still the most deformed part of the slope. The acceleration was greatest during the rainy seasons as water infiltrating the cracks of the upper slope gradually increased the pressure from top to bottom. This shows that fissure water caused by vertical rainfall infiltration plays an important role in slope stability.

The impact of blasting vibration on the landslide was mainly in the area between the front edge and the middle part of the landslide (Figure 9a). There were five blasting cycles from July to the end of August 2018, and the blasting charge each time was 3480 kg. During this period, the displacement change near the front edge of the landslide increased at a uniform speed, and the sliding deformation increased significantly. In December 2018, there were 16 blasts 1140 m from the main monitoring profile, with an average charge of 1992 kg. The maximum charge was 2496 kg on 15 December. However, the displacement of the landslide mass did not increase significantly at this stage. The relationship between YII4-1 and YII4-5 displacement outside the boundary of the landslide mass and blasting was not obvious (Figure 9b). Excavation by blasting damaged and loosened the strongly weathered rock mass. This served to promote the expansion of tension cracks that were generated at the top of the slope, providing a more convenient channel for rainwater infiltration and adversely effecting slope stability.

3.2. Comprehensive Prediction Results of Landslide Risk

3.2.1. The LSTM Model for Deformation Process Prediction

The monitoring data obtained from surface monitoring stations II5-1, II4-3, II3-3, and II2 and the deep-seated displacement monitoring data of CX01 at different elevations were employed to establish the forecast LSTM model for the Fengning landslide. In the LSTM model, the monitoring sequence is divided into two parts. One part of the data is used for model training to obtain the time series displacement prediction model. Then the model is used to test the other part of the data, and the reliability of the model is evaluated through the correlation analysis between the predicted value of the model and the measured value.

By constantly selecting the size of the training set, the monitoring data set in this study was divided into the training set and the test set according to the proportions of 78% and 22%, which are 352 × 0.78 ≈ 272 and 352 × 0.22 ≈ 80, respectively. The window period was set at 30 (that is, a sequence of 30 consecutive time points was used to predict the value of the next time point). The model accuracy evaluation index also adopted the RMSE, R, MAE, and MAPE. The obtained model training and prediction curves are shown in Figure 10. For example, using surface monitoring station YII5-1 and the CX01 inclinometer at an elevation of 1205 m, the LSTM model predicted that the surface displacement of YII5-1 would reach 568.92 mm by 25 February 2022 (Figure 10e). The horizontal displacement of the CX01 inclinometer was predicted to reach 1202 mm by 10 August 2019 (Figure 10f).

3.2.2. Prediction of Influence Range of Landslide Mass by Using Discrete Element Method

Displacement Back Analysis of the Joint Persistence Ratio

In fact, the sliding surface parameters of landslide mass are variables with the displacement. At present, they are assumed to be constant in the stability evaluation of landslide, which brings unrealistic results to the stability prediction of landslide. Therefore, based on the displacement back analysis method and discrete element numerical simulation method, the strength index of sliding surface can be obtained through the displacement prediction model, which provides a reasonable basis for stability evaluation.

The three-dimensional discrete element numerical model of the excavated slope was established using the software 3DEC (Figure 11a,b). Select 10 groups of profiles along the strike direction of 110° according to the actual terrain, establish a slope model similar to the actual terrain in surfer software Surfer 14.0, and then import it into 3DEC software for modeling. The selected area elevation is from 1100 m to 1337 m, and the model height is 237 m; model width 144 m. The length of the model is 275.2 m. In the model, the x-axis of the Cartesian coordinate system points horizontally into the slope, the y-axis points horizontally along the slope to the downstream side, and the z-axis points vertically upward. There are 33,798 units and 6647 nodes in total. Complete granite block and joint parameters are selected for calculation, and the parameters are shown in

Table 2. Physical and mechanical parameters of rock mass.

Rock Mass	Φ (°)	c (kPa)	Natural Density (kN/m3)	Saturated Density (kN/m3)	Poisson Ratio (ν)	Elastic Modulus (×104 MPa)
Strongly weathered granite	27	50	24.1	24.6	0.38	0.50
Weakly weathered granite	31	350	25.0	25.3	0.23	0.60
Slightly weathered granite	35	1050	26.0	26.2	0.18	0.75

and

Table 3. Joint parameters.

Joint Type	Cohesive c (MPa)	Internal Friction Angle φ (°)	Normal Stiffness (Pa/m)	Tangential Stiffness (Pa/m)
Fault	0.024	16	4 × 109	9.5 × 108
joint	0.048	20	4 × 109	1.1 × 109

. The yield criterion followed the Mohr–Coulomb criterion, and the fault and structural plane are set as the Coulomb slip model. On the four side surfaces, displacements along the directions normal to the surfaces were also fixed. The top surface was the free boundary. The stability coefficient is solved by strength reduction method in 3DEC. In the numerical simulation, the slope instability criteria are mainly as follows: (1) take the convergence of numerical calculation as the instability criterion; (2) The catastrophe of the displacement of the characteristic part is taken as the instability criterion; (3) The continuity of plastic zone is used as the instability criterion.

Meanwhile, considering the layout of the monitoring points in the sliding stage, YII5-1 and YII4-3 were chosen as the basis for the displacement back analysis. In order to establish the sample for displacement back analysis, according to the on-site joint investigation and persistence calculation, the minimum JPR of rock mass on the excavation surface is 0.56. Therefore, 40 groups of JPR value ranging from 0.60 to 1.00 and interval of 0.01 were substituted into the numerical model to calculate the displacement of the landslide under rainfall and blasting. The influence of rainfall on slope stability is mainly simulated from the change of void water pressure caused by rainfall infiltration, and the influence of blasting on slope deformation is simulated according to the method discussed in Section 2.4. The two monitoring points are YII4-3 and YII5-1 respectively. The displacement changes of the two points in X, Y and Z directions are monitored respectively. The location of the monitoring points is shown in Figure 11c. The numerical calculation model is used to obtain the sufficient training samples in LSTM model.

The numerical simulation results of the sliding process are shown in Figure 12. The displacement under action of rainfall and blasting could be obtained using numerical analysis. After 4000 steps of calculation iteration of point YII4-3, the total displacement in X direction is about 160 mm and that in Z direction is about 160 mm, which is consistent with the slope rapid deformation stage. After iteration to 8000 steps, the displacement in X direction is about 250 mm and the displacement in Z direction is about 210 mm, which corresponds to the slope deformation in the second stage. After iterating to 10,000 steps (Figure 12a,b), the displacement in X direction reaches about 320 mm, and the displacement in Z direction is 360 mm, which is consistent with the slope deformation in the third stage, that is, it enters the rapid deformation stage again.

At point YII5-1, after 4000 steps of calculation iteration, the total displacement in X direction is about 200 mm and that in Z direction is about 120 mm, which is also consistent with the slope deformation in the first stage. After iteration to 8000 steps, the displacement in X direction is about 300 mm and the displacement in Z direction is about 210 mm. After iterating to 10,000 steps (Figure 12c,d), the maximum displacement in X direction reaches about 450 mm, and the displacement in Z direction is 260 mm, which is consistent with the slope deformation in the third stage, that is, it enters the rapid deformation stage again.

The calculation iteration reaches 12,000 steps, and the slope displacement in each direction increases to varying degrees compared with 10,000 steps. The maximum displacement is about 1m, and most of the displacement is about 500 mm–600 mm.

In this study, the displacements of the monitoring point calculated numerically were selected as the input layer in LSTM model, and the JPR parameter was regarded as the output layer to train the LSTM model. The LSTM is helpful in developing the relationship between input parameters and the outcome of the model. The calculation results are shown in

Table 4. Three dimensional discrete element numerical inversion results.

Route of Failure Surface	Cohesion c (MPa)	Internal Friction Angle φ (°)	JPR (k)
Shear outlet at EL.1230	0.218	23.976	0.85
Shear outlet at EL.1215	0.158	22.598	0.90
Shear outlet at EL.1185	0.108	21.429	0.95

.

According to the JPR of the sliding surface at the deepest 1185 shear outlet obtained from the above inversion, the factor of safety program is compiled by using the fish function in the three-dimensional discrete element, and the slope stability coefficient is calculated to be 0.975. The slope has entered the stage of long-term creep deformation now.

Rockslide Damage Influence Scope Prediction by Using Discrete Element Method

On the basis of considering the size and distribution position of granite blocks, divided the blocks into two basic forms: sliding and rolling. The contact between blocks is edge to edge. The damage range of instability along 1185 m sliding shear outlet under rainfall and blasting is analyzed by three-dimensional discrete element method. The top elevation of the numerical model is 1340 m and the bottom elevation is 1050 m (Figure 13a). The length of the model in X direction is 488 m, and the length in Y direction is 220 m. The horizontal direction is constrained at both sides of the boundary, the bottom boundary is constrained in the horizontal and vertical directions, and the upper boundary of the model is a free surface. The blocks were assumed as rigid, the size of the block is determined by the spacing of joints, and the block does not disintegrate during sliding. At the same time, considering the change of the JPR of the sliding surface with the displacement, the calculation results of different time steps were shown in

Table 5. Calculation results of damage influence range of 1185 m shear outlet.

Calculation Time Step T	Sliding Distance in Sliding Direction (m)	Accumulation Width of Landslide Mass (m)	Front Elevation of Landslide Accumulation Mass (m)
7000	Tension slide start		1185
10,000	46	126	1170
13,000	94	168	1140
16,000	132	183	1100
19,000	135	190	1100
22,000	138	198	1100
25,000	146	211	1100
28,000	152	216	1100

and Figure 13.

At T = 28,000, the sliding failure has stopped, the fragmented blocks in front of the landslide accumulation are distributed in a near fan shape, slightly inclined to the upstream side, the farthest block slides to the 1100 m platform, the maximum sliding distance along the slope is 152 m, and the nearest place to the river is 129 m, so it has no impact on the river and Highway No.9. The maximum width scattered to both sides is 216 m, the maximum distance of sliding to the upstream side is 59 m, and 159 m away from the buildings on the upstream side, which has no impact on the buildings. The widest part is located on the 1185m platform. At 1100 m, the dispersion width of the platform is 161 m, and the ratio of slope height to the migration distance of landslide accumulation is 1:1.04.

Through the comparative analysis of the influence range of the landslide mass sliding along the three sliding surfaces (Figure 14), the results show that when sliding along the sliding surface in the middle, the influence range is the largest due to the large volume and potential energy of the landslide mass. When sliding occurs on the lower sliding surface, the sliding distance is the smallest, the damage influence range is the smallest and the accumulation thickness is the largest due to the influence of landslide potential energy and excavated slope platform. Therefore, when carrying out landslide safety warning, the warning shall be carried out according to the influence range of sliding along the 1200 sliding surface.

4. Discussion

4.1. Mechanisms of the Stepped Rock Mass Landslide

The relaxation joints along the slope simulated a bedding-plane along the base of the rock mass structure, readily forming a steep-gentle stepped sliding zone. JPR plays a controlling role in the stability of a rockslide. The sliding surface of the step-path landslide was composed of three sections. As discussed in Section 3.1.3, the landslide was mainly triggered by rainfall and excavation. Movement of the lower part was caused by excavation relaxation and rainfall infiltration. The displacement was fast, which reflects the relaxation traction effect and tensile effect on the landslide mass. At the trailing (upper) edge, after the first stage of sliding, the rock mass was broken, which resulted in a pushing effect on the lower part of the landslide mass and acted on the sliding surface with compressive stress (Figure 15a). The sliding surface in the middle was stepped, and the existence of steep cracks weakened the tensile strength or shear resistance of the rock bridge. The different rainfall conditions induced deformation growth, resulting in the formation of a stepped sliding route.

As a landslide moves, the gently inclined joints produce shear sliding under compressive shear stress, the steeply inclined joints produce shear sliding under tensile shear stress, and the rock bridge produces shear or tensile fractures under tensile shear stress. The stepped failure of a rock slope is the result of multiple sliding shears of a potentially deformed slope along steeply and gently inclined joints in stepwise fashion. Owing to the combined stepped distribution of steeply and gently inclined fractures, there are multiple sliding shear planes, having the failure characteristics of multi-stage shear and step-by-step tensile fracture of the rock bridge under the action of shear force (Figure 15b). After the rock bridge is pulled step-by-step, a landslide mass with a stepped fracture surface is formed. Because the sliding surface still has a certain shear strength and the rock bridge closest to the slope top is not completely connected, the entire sliding mass will not shear down and slide out immediately (Figure 15c). During the sliding shear process of the landslide mass, the slope mass moved stepwise from top to bottom along the steeply and gently inclined fracture surfaces (Figure 15d,e).

4.2. Applicability and Accuracy of LSTM Method

The root mean square error (RMSE), mean absolute error (MAE), and regression coefficient (R) were used to compare the model reliability. Figure 16 and Table 5 show that the proposed LSTM model provided a very good prediction of the measured behavior of the landslide surface displacement, with a correlation coefficient of 0.977 along with average and maximum relative errors of 1.91% and 0.87%, respectively. These are satisfactory results for critical early warning before landslide treatment. To test the validity and reliability of the LSTM model proposed in this paper, a comparative study with BP (BP neural network), ARMA (autoregressive moving average), and SVM (support vector machine) models were conducted. The results are shown in Figure 16.

Figure 16 shows that after training the LSTM model, the MAE and RMSE converged to within 2.00 mm. The RMSE and region coefficient R were better than those of the other two artificial intelligence methods and the traditional time series model, ARMA. However, the BP neural network showed unstable results for different deformation monitoring sequences.

Table 6. Evaluation indexes of the testing set of surface monitoring points for the LSTM model.

Monitoring Station	RMSE (mm)	R	MAE (mm)	MAPE (%)
YII1	1.20	0.956	1.29	0.87
YII2	0.52	0.977	0.82	0.19
YII3-3	0.67	0.976	0.65	0.16
YII4-3	1.35	0.936	1.24	0.22
YII5-1	1.09	0.962	1.91	0.35

and Figure 16 shows that the LSTM model has a stronger learning and fitting ability for training set samples than traditional machine learning methods, especially for the prediction results of the test set. LSTM trained in the input sample (training set) showed better adaptability in the fresh sample (test set), and its prediction accuracy and prediction ability were better than those of ARMA and the traditional machine learning models, BPNN and SVM. Thus, LSTM is a suitable and reliable method for predicting rockslide displacement. The predicted values showed good agreement with the simulated values. This comparison indicates that the LSTM can be used to predict the future displacement of a stepped failure landslide.

From Figure 10 and Figure 16, the prediction accuracy of the LSTM model is much higher than that of other models for both surface and deep displacements. The accuracy of the LSTM prediction model depends on the parameters, mainly the epoch (number of training rounds), batch_size, and time window.

This study trained different epochs and determined the influence of different epochs on the simulation results (Figure 17a). When the epoch was 200, the resulting model error was the smallest. Because the epoch is difficult to determine, the number of iterations is usually determined by the value of the loss function (Figure 17b).

The batch_size refers to the number of samples (n) selected in a specific iteration. One iteration in the network first randomly selects n samples to form one batch, and then inputs the batch into the LSTM model to obtain the output results. Different batch sizes affect the accuracy of the model. If epoch remains unchanged (e.g., the epoch is fixed at 200), the LSTM model error MAE generally shows an upward trend with the increase in batch_size, which indicates that the prediction accuracy gradually decreases. Usually, batch_size is set to a power of 2. Figure 15c shows that different batch sizes affected the prediction accuracy of the YII3-3 monitoring point.

This study adopted the “sequence to sequence” method in LSTM prediction. That is, using the displacement data of n consecutive time series as the input data, the model will output m displacement data (where n is the size of the time window). In the LSTM model, n can be different from m. Figure 15d shows the variation law of model accuracy by changing the length of the time window of monitoring station YII3-3. When the window was in the interval [26,36], the model error MAE was the lowest. When the window was greater than 36, the MAE error increased gradually. A window of 30 was selected as the best choice for predicting the YII3-3 point displacement time series. Similar laws were obtained for other surface displacement monitoring stations for the landslide mass. The optimal LSTM prediction model was when batch_size was 16, epoch was 200, and n was 30. One disadvantage of the LSTM model is that requires multiple training cycles to obtain the best model parameters.

5. Conclusions

Predictions of landslide deformation trends and influence scope are important bases for landslide treatment and disaster prediction. In this study, we presented a modeling approach for landslide disaster prediction using LSTM and numerical method on considering rainfall and explosive action. From this study, the following conclusions can be drawn:

(1) For a stepped sliding surface landslide formed by a combination of steeply–gently inclined joints near the front edge of the landslide mass, the deformation was large owing to the influence of slope toe excavation, rainfall, and blasting vibration (monitoring stations YII5-1, YII5-2, and YII4-3). However, the deformation was small near the middle (YII3-3) and upper part (YII2) of the landslide mass. The deformation of various parts of the landslide mass was controlled by the JPR in the stepped sliding surface.

(2) The failure mechanisms of the Fengning landslide were: excavation and relaxation at the slope toe, blasting vibration, and rainfall infiltration. The joint persistence ratio of sliding surface can be predicted by LSTM model combined with the back analysis of discrete element displacement. The predicted results can reflect the external load and different deformation stages, and can be used to determine the strength parameters of sliding surface in different sliding time.

(3) Compared with other deep learning models, the LSTM model does not need to separate the periodic term from the monitoring sequence and can effectively solve the problem of gradient appearance or gradient expansion in a general RNN. The simulation accuracy is higher than that of the other three types because of the great advantages of LSTM in data processing. This shows that the LSTM is more suitable for the displacement prediction of rock landslides with stepped sliding surfaces under the action of external uncertain interference factors. The LSTM method is adopted to predict future phased deformation of the Fengning landslide. One disadvantage of the LSTM method is that is requires multiple training cycles to obtain the best model parameters.

(4) The discrete element numerical analysis method shows that the potential energy, volume, connection condition of sliding surface and free condition are important factors affecting the distribution of landslide accumulation. Sliding along the 1200 m shear outlet has the greatest influence range.