**1. Introduction**

Landslide geological disasters are a serious type of geological disaster that occur worldwide, inducing serious threats and losses to the development of human society. In recent years, under the influence of extreme global climate change, seismic activities, coupled with the rapid development of human engineering activities, have become more intense interferences to the natural environment, directly leading to geological disasters with greater intensity and higher frequency [1,2]. This increases the difficulty of developing landslide disaster reduction strategies [3]. Combined with land use change, population growth, uncontrolled urbanization, and so on in vulnerable areas, the landslide risk level continues to rise [4]. The Three Gorges Reservoir Area is one of the areas with a high incidence of landslide disasters in China [5,6]. Every year, landslides cause very large losses to local lives and property, so it is necessary to prevent and control such landslide disasters [7].

Landslide displacement refers to the distance that the soil or rock mass on the slope slides along the slope as a whole or separately under the influence of gravity under the influence of river erosion, groundwater activity, rain immersion, earthquakes, and artificial

**Citation:** Lin, Z.; Ji, Y.; Liang, W.; Sun, X. Landslide Displacement Prediction Based on Time-Frequency Analysis and LMD-BiLSTM Model. *Mathematics* **2022**, *10*, 2203. https:// doi.org/10.3390/math10132203

Academic Editors: Xiang Li, Shuo Zhang and Wei Zhang

Received: 24 May 2022 Accepted: 22 June 2022 Published: 24 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

slope cutting [8]. The generation of landslide displacement is influenced by both the internal geological conditions (geological structure, landform, lithology, etc.) and external influencing factors (rainfall, reservoir water level, etc.) of a given landslide location [9–11]. Landslide displacement prediction is a frontier problem in the international landslide research field [12]. Landslide displacement prediction is an important part of landslide disaster loss reduction [13]. It can be used to summarize the historical displacement of landslides and environmental conditions and also to analyze the potential relationship between geological and meteorological environmental changes and disasters. The combination of an accurate landslide displacement prediction model and landslide warning model can effectively improve people's ability to determine landslides in daily life, help decision-makers make more accurate decisions [14], and take active disaster reduction actions in advance [15] to achieve adaptive risk avoidance [16] and protect people's lives and health, and to achieve the purpose of improving people's livelihood [14]. In the Three Gorges Reservoir area, affected by rainfall and reservoir water level, landslide displacement usually shows the characteristic of a "step shape" that represents accelerated activity in the rainy season and remains almost steady in the dry season [17]. Therefore, defining methods for predicting the increase in landslide displacement has become the focus of scholars. However, landslides are very complex systems, and their deformation is affected by their own engineering geological conditions and externally induced factors [14]. The displacement curve is often highly nonlinear, which makes it difficult to accurately predict the landslide displacement [8].

At present, landslide displacement prediction models are mainly divided into physicsbased and data-based models [18]. Although both models can predict landslide displacement, physics-based models are complex, time-consuming, expensive [19], difficult to establish [5], and have strict application conditions [11], which can only be used in limited cases [20]. However, the data-based model has a simple process, accurate prediction, and low cost [19], and it is good at dealing with nonlinear relations [21]. Therefore, physicsbased models are not as popular as data-based models [10]. Thus, most of the landslide displacement prediction models in recent years are based on data models. The key factors for the occurrence of landslides can be roughly divided into two categories: the slope, lithology, and soil type are internal factors affecting landslides, while the rainfall, reservoir water level, and snowfall are external factors affecting landslides [22].

Based on the principle of time series analysis, many studies have decomposed landslide displacement into trend displacement and periodic displacement [5,8,10,14,15,23], which has well separated the nonlinear characteristics of landslide displacement. There have also been studies on the decomposition of landslide displacement data into several subsequences of different frequencies based on the time-frequency analysis method. Guo et al. [24] combined the variational modal decomposition (VMD) method with the WA-GWA-BP model, and Liu et al. [6] combined the VMD method with the periodic neural network (PNN) model to predict landslide displacement. The empirical mode decomposition (EMD) method was combined with the LSTM model, and a linear interpolation technique was used to increase the size of the training dataset to accurately predict landslide displacement [25]. According to the wavelet transform, multiscale analysis can be carried out on landslide displacement data through the operation functions of stretching and shifting [26–28]. To solve the problems of modal aliasing in EMD, the Ensemble Empirical Mode Decomposition (EEMD) algorithm was used to decompose the landslide data series, further improving the prediction ability of the model [29–31]. Taking into consideration the fact that EMD and EEMD have randomness and uncontrollability built in the decomposition times of landslide displacement, Xing et al. [32] used VMD to decompose landslide displacement.

Taking into consideration the lag fluctuation of the groundwater level, the SVC-PSO-SVR model was proposed to predict landslide displacement by Han et al. [33]. Deng et al. [34] used acoustic emission sound generation and rainfall as data inputs, and the equivalent reservoir water level function model was combined with Lasso-ELM to improve

the accuracy of landslide displacement prediction. Based on the traditional gray prediction model, L et al. improved and proposed a new gray prediction model [13]. SVR and the Hausdorff derivative operator were used to determine model parameters with the improved SALP group algorithm, further improving the prediction performance of the traditional gray model [35]. Based on the optimal weight allocation method, Li et al. [36] assigned different weights to the verhulst model and GM(1,1) model, and combined the advantages of the two models to form a new prediction model. Considering the importance of past experience, Hu et al. [37] used the verhulst inverse function to describe the motion characteristics of landslides, and constructed a displacement prediction model combined with a random forest algorithm. To predict the displacement more accurately, two new concepts, the trend sequence and sensitivity state, were proposed, and a new model was obtained by integrating the trend sequence and sensitivity state [38]. The cost function and the penalty mechanism were proposed in order to force the underestimated landslide displacement to be transferred to a higher estimate, and the ability of the model to avoid landslide risk is taken into full consideration while predicting landslide displacement [16].

Researchers considered that it is normal that the landslide displacement will fluctuate within the normal range in the future, and the prediction interval method was adopted instead of point prediction; this method can obtain clear data, and the resulting model was presented as an interval [15,29,39,40]. Interval prediction can not only predict the future variation trend of data but also obtain the variation range of landslide displacement, providing an important basis for decision-making regarding landslide disaster prevention and mitigation [21]. In recent years, with the development of technology, the local mean decomposition (LMD) algorithm has been increasingly applied, and an increasing number of studies have begun to use the LMD algorithm to address nonlinear problems in various fields [41–45]. Inspired by previous studies, this paper presents a new theory that intends to apply the LMD algorithm to process typical nonlinear landslide displacement data on the basis of previous research results and time-frequency analysis methods and to propose a new local mean decomposition-bidirectional long short-term memory (LMD-BiLSTM) model.

The LMD-BiLSTM model can decompose nonlinear and nonstationary data series well, it can solve the problem of ignoring random displacement in time series analysis, and the LMD algorithm is used to decompose the original displacement data and influencing factors into several sub-time-series data of different frequencies. Due to the complex relationship between influencing displacement factors and landslide displacement, the correlation between them cannot be well quantified. Therefore, to improve the accuracy of prediction, the maximum information coefficient (MIC) algorithm is introduced in this paper.

MIC correlation calculations are carried out between the subsequence of influencing displacement factors and each of the subsequences of landslide displacement, and the data with high correlation are selected as the input variable of each subsequence prediction, which improves the validity and reliability of the input data. Considering that the rainfall and reservoir water level of landslides have a similar change trend in each time period, the bidirectional long short-term memory (BiLSTM) model is used to predict each sub-series and the final predicted displacement of the landslide is obtained by adding the predicted results. The case of a Baishuihe landslide in the Three Gorges region of China is taken to verify the prediction performance and advantages of the model.

The main contributions of this study are as follows.

1. In order to solve the problem where the time series analysis method ignores the random factors, which would affect the accuracy of prediction, based on the timefrequency analysis method, the LMD algorithm is used for the first time to decompose landslide displacement data and factors affecting displacement into multiple instantaneous frequency subsequences with physical significance. The disadvantages of EMD or EEMD mode aliasing and endpoint effects are solved, and the integrity of the signal is better preserved.


#### **2. Materials and Methods**

#### *2.1. Local Mean Decomposition*

The LMD algorithm was first proposed by Jonathan S. Smith [46], and its advantages over EMD in EEG signal processing were discussed. In recent years, the LMD algorithm, a decomposition method for nonlinear nonstationary signals, has been applied to tool chatter areas [41], solar radiation prediction [42], underwater acoustic signal processing [43], daily natural gas load forecasting [44], time series compressor stall processes [45], and many other fields.

The core idea of the LMD algorithm is to adaptively decompose a complex nonstationary multicomponent signal into the sum of several product functions (*PF*) with the physical significance of instantaneous frequency, for which each *PF* component is directly calculated by an envelope signal and a pure frequency modulation (FM) signal. The envelope signal is the instantaneous amplitude of the *PF* component, and the instantaneous frequency of the *PF* component can be directly calculated from the pure frequency modulation signal. Furthermore, by combining the instantaneous amplitude and instantaneous frequency of all *PF* components, the complete time-frequency distribution of the original signal can be obtained. The LMD algorithm can decompose complex nonlinear landslide displacement data into several *PFs* with physical significance, and its decomposition steps are as follows [43,46]:

$$m\_i = \frac{n\_i + n\_{i+1}}{2} \tag{1}$$

$$\alpha\_{i} = \frac{|n\_{i} - n\_{i+1}|}{2} \tag{2}$$

where *m<sup>i</sup>* is the *i*-th average value of two consecutive extreme values *n<sup>i</sup>* ; *ni*+1*α<sup>i</sup>* is the local magnitude of each halfwave oscillation; *i* = 1, 2, . . . , *m* − 1 (where m is the number of extrema). The mean *m<sup>i</sup>* of all continuous extreme values must form a straight line. The moving average method is used to smooth the mean *m<sup>i</sup>* and *α<sup>i</sup>* , and the local mean function *m*11(*t*) and the envelope estimation function *α*11(*t*) are obtained. Then, *m*11(*t*) is separated from the original signal *X*(*t*) to obtain the residual signal *h*11(*t*):

$$h\_{11} = X(t) - m\_{11}(t) \tag{3}$$

Then, *h*11(*t*) is divided by *α*11(*t*) to conduct amplitude modulation and yield:

$$s\_{11}(t) = \frac{h\_{11}(t)}{\mathfrak{a}\_{11}}\tag{4}$$

where *s*11(*t*) is a pure frequency modulation signal, and the process is repeated q times until a pure FM signal *s*1*q*(*t*) whose envelope function meets *α*1(*q*+1) (*t*) = 1 is obtained. The corresponding envelope *α*1(*t*) is shown in the formula:

$$a\_1(t) = a\_{11}(t)a\_{12}(t) \cdots a\_{1n}(t) = \prod\_{q=1}^{n} a\_{1q}(t) \tag{5}$$

where *q* is the number of iterations and the first component *PF*1(*t*) is:

$$PF\_1(t) = s\_{1q}(t) \times a\_1(t) \tag{6}$$

The *PF*1(*t*) is subtracted from the original signal *X*(*t*) to obtain the new signal. The above steps are repeated to obtain *PF*2(*t*). These steps are repeated until the last signal becomes a constant or contains no more oscillations, and then the residual signal *u<sup>k</sup>* (*t*) can be obtained. Therefore, the original signal *X*(*t*) can be decomposed into the sum of the *PF* component and *u<sup>k</sup>* (*t*), as shown in the formula below:

$$X(t) = \sum\_{p=1}^{k} PF\_p(t) + \mu\_k(t) \tag{7}$$

#### *2.2. Maximal Information Coefficient*

Reshef et al. [47] proposed the MIC concept that can be used to measure the nonlinear correlation between two variables based on mutual information theory in information theory. The main idea is as follows: if there is some correlation between the two variables, then after some sort of meshing on the scatter plot of the two variables, the mutual information of these two variables can be calculated according to the approximate probability density distribution of the variables in the grid. After regularization, this value can be used to measure the correlation between the two variables. When the MIC is 0, it indicates that the pairs of variables are completely independent; when the MIC is 1, it indicates that there is some functional relationship between the pairs of variables. The larger the MIC is, the stronger the correlation between variables is. For a set of ordered pair datasets *D* = {(*x<sup>i</sup>* , *yi*), *i* = 1, 2, . . . , *n*}, if the X-axis is divided into X cells and the Y-axis into Y cells, the result is an *x* × *y* grid *G*. The points in the dataset *D* land on *G* based on the proportion of approximate *D*|*G*, representing its probability distribution. Therefore, the maximum mutual information can be defined as follows:

$$I^\*(D, \mathfrak{x}, y) = \max I(D|G) \tag{8}$$

where the maximum value in the above formula is the maximum value of mutual information on *G* of all possible networks in which *D* divides the X-axis into X grids and the Y-axis into Y grids. *I*(*D*|*G*) indicates mutual information in the case of probability distribution *D*|*G*. The elements of the *x*-th row and *y*-th column of the eigenmatrix *M*(*D*) on the ordered pair dataset *D* are shown in the formula:

$$M(D)\_{\mathbf{x},\mathbf{y}} = \frac{I^\*(D,\mathbf{x},\mathbf{y})}{\log(\min\{\mathbf{x},\mathbf{y}\})} \tag{9}$$

Then, the ordered pair dataset *D* divided by data scale n and the number of grids is less than or equal to *B*(*n*). Then, the MIC of dataset *D* is defined as follows:

$$MIC(D) = \max\_{\mathbf{x}y < B(n)} \{ M(D)\_{\mathbf{x}, \mathbf{y}} \} \tag{10}$$

According to the actual application [47–49], *B*(*n*) = *n* 0.6 is recommended, so this paper also chooses this value.

#### *2.3. Bidirectional Long Short-Term Memory Model*

Landslide displacement is affected by external factors such as rainfall and reservoir level, and the collected data series have certain fluctuations and periodicities. To analyze and predict the time series data of landslide displacement, this paper uses an improved recurrent neural network (RNN) model to predict landslide displacement. A long shortterm memory network (LSTM) was proposed by Hochreiter et al. [50]; it is a special RNN

that, because of the added gate mechanism that is added, can solve the RNN gradient explosion and gradient disappearance problems to a certain extent. that, because of the added gate mechanism that is added, can solve the RNN gradient explosion and gradient disappearance problems to a certain extent.

Landslide displacement is affected by external factors such as rainfall and reservoir level, and the collected data series have certain fluctuations and periodicities. To analyze and predict the time series data of landslide displacement, this paper uses an improved recurrent neural network (RNN) model to predict landslide displacement. A long shortterm memory network (LSTM) was proposed by Hochreiter et al. [50]; it is a special RNN

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 6 of 20

*2.3. Bidirectional Long Short-Term Memory Model* 

On the basis of a traditional RNN, the LSTM model introduces *c<sup>t</sup>* to store long-term memory. The model can be made to selectively update the memory unit so that the current node learns the characteristics of the input data from a long time ago. There are three control gates in each LSTM unit, namely the input gate *i<sup>t</sup>* , forget gate *f<sup>t</sup>* , and output gate *o<sup>t</sup>* . The structure of the LSTM model is shown in Figure 1. On the basis of a traditional RNN, the LSTM model introduces ௧ to store long-term memory. The model can be made to selectively update the memory unit so that the current node learns the characteristics of the input data from a long time ago. There are three control gates in each LSTM unit, namely the input gate ௧, forget gate ௧, and output gate ௧. The structure of the LSTM model is shown in Figure 1.

**Figure 1.** The LSTM model structure. **Figure 1.** The LSTM model structure.

̃

Assuming the input sequence = {ଵ, ଶ,…௧}, ௧ = {௧ଵ, ௧ଶ,…௧} represents the real vector data of the k-dimension under time step t, the LSTM model is used to construct Assuming the input sequence *x* = {*x*1, *x*2, . . . *xt*}, *x<sup>t</sup>* = {*xt*1, *xt*2, . . . *xtk*} represents the real vector data of the k-dimension under time step t, the LSTM model is used to construct landslide displacement data, and the updating formula of each internal unit is as follows.

landslide displacement data, and the updating formula of each internal unit is as follows. The forget gate ௧ is used to forget the information of the state ௧ିଵ of the upper memory unit, and its formula is The forget gate *f<sup>t</sup>* is used to forget the information of the state *ct*−<sup>1</sup> of the upper memory unit, and its formula is

$$f\_t = \sigma \left(\mathcal{W}\_f \mathbf{x}\_t + \mathcal{U}\_f h\_{t-1} + b\_f\right) \tag{11}$$

where is the weight matrix of the forget gate, is the bias of the forget gate, and is the sigmoid function. where *W<sup>f</sup>* is the weight matrix of the forget gate, *b<sup>f</sup>* is the bias of the forget gate, and *σ* is the sigmoid function.

The calculation of the candidate state of memory unit ̃ <sup>௧</sup> is shown in the formula, and the input gate ௧ determines the reserved information of the candidate state in the current unit. The calculation of the candidate state of memory unit *<sup>c</sup>*e*<sup>t</sup>* is shown in the formula, and the input gate *i<sup>t</sup>* determines the reserved information of the candidate state in the current unit.

$$\widetilde{\boldsymbol{c}}\_{t} = \tanh(\mathcal{W}\_{\mathcal{C}}\mathbf{x}\_{t} + \mathcal{U}\_{\mathcal{C}}h\_{t-1} + b\_{\mathcal{C}}) \tag{12}$$

$$\dot{\mathbf{u}}\_{l} = \sigma(\mathbf{W}\_{l}\mathbf{x}\_{l} + \mathbf{U}\_{l}\mathbf{h}\_{l-1} + \mathbf{b}\_{l}) \tag{13}$$

௧ = (௧ + ℎ௧ିଵ + ) (13) where and represent the weight matrices of the input gate ௧ and candidate state <sup>௧</sup> , respectively; and represent the corresponding bias quantities of ௧ and ̃ <sup>௧</sup> . ௧ where *W<sup>i</sup>* and *W<sup>c</sup>* represent the weight matrices of the input gate *i<sup>t</sup>* and candidate state *c*e*t* , respectively; *<sup>b</sup><sup>i</sup>* and *<sup>b</sup><sup>c</sup>* represent the corresponding bias quantities of *<sup>i</sup><sup>t</sup>* and *<sup>c</sup>*e*<sup>t</sup>* . *i<sup>t</sup>* and *f<sup>t</sup>* combine with the previous memory state *<sup>c</sup>t*−<sup>1</sup> and the current candidate state *<sup>c</sup>*e*<sup>t</sup>* to update the current memory unit state *c<sup>t</sup>* .

$$\mathcal{c}\_{t} = f\_{t} \odot \mathcal{c}\_{t-1} + i\_{t} \odot \widetilde{\mathcal{c}}\_{t} \tag{14}$$

where  indicates multiplication. The input gate *o<sup>t</sup>* is mainly used to control the output of the memory unit state value.

$$\rho\_t = \sigma(\mathcal{W}\_o \mathbf{x}\_t + \mathcal{U}\_o h\_{t-1} + b\_o) \tag{15}$$

The weight and bias of each unit in the above types are dynamic and can be updated through data training to predict landslide displacement in time series. In traditional LSTM models, the information is one-way, and the model can use past information but not future information. To adapt to the variation characteristics of landslide displacement, precipitation, and reservoir water level, the BiLSTM model was selected to construct the prediction model. to update the current memory unit state ௧. ௧ = ௧ ⊙ ௧ିଵ + ௧ ⊙ ̃ ௧ (14) where ⊙ indicates multiplication. The input gate ௧ is mainly used to control the output of the memory unit state value.

and ௧ combine with the previous memory state ௧ିଵ and the current candidate state ̃

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 7 of 20

The BiLSTM model is formed by the combination of positive and reverse LSTM [51], and its structure is shown in Figure 2. Forward LSTM can obtain the past data information of the input sequence, and backward LSTM can obtain the future data information of the input sequence [52]. The forward and backward LSTM training processes of time series data can further improve the global integrity of feature extraction. At time t, the output value *H<sup>t</sup>* of the hidden layer of BiLSTM is composed of forward → *h <sup>t</sup>* and backward ← *h <sup>t</sup>* : ௧ = (௧ + ℎ௧ିଵ + ) (15) The weight and bias of each unit in the above types are dynamic and can be updated through data training to predict landslide displacement in time series. In traditional LSTM models, the information is one-way, and the model can use past information but not future information. To adapt to the variation characteristics of landslide displacement, precipitation, and reservoir water level, the BiLSTM model was selected to construct the pre-

$$\stackrel{\rightarrow}{h}\_{t} = \begin{array}{c} \overrightarrow{\text{LSTM}} \ (h\_{t-1}, \mathbf{x}\_{t}, c\_{t-1}), t \in [1, T] \end{array} \tag{16}$$

$$\stackrel{\leftarrow}{h}\_{t} = \overleftarrow{\text{LSTM}} \ (h\_{t+1}, \mathbf{x}\_{t}, c\_{t+1}), t \in [T, 1] \tag{17}$$

$$H\_t = \begin{bmatrix} \stackrel{\rightarrow}{h}\_{t\prime} \stackrel{\leftarrow}{h}\_{t\prime} \end{bmatrix} \tag{18}$$

<sup>௧</sup> = ⃖ሬሬሬሬሬሬሬሬሬሬሬ(ℎ௧ାଵ, ௧, ௧ାଵ), ∈ ሾ, 1ሿ (17)

௧

**Figure 2.** The BiLSTM model structure. **Figure 2.** The BiLSTM model structure.

ℎ → <sup>௧</sup> = ሬሬሬሬሬሬሬሬሬሬሬ⃗(ℎ௧ିଵ, ௧, ௧ିଵ), ∈ ሾ1, ሿ (16) After obtaining *H<sup>t</sup>* , the landslide displacement is obtained by a fully connected layer as the prediction output of the model.

#### **3. Results**

**3. Results**  *3.1. A Real Case* 

#### *3.1. A Real Case*

diction model.

௧ = ሾℎ→ ௧, ℎ ← ௧ሿ (18) After obtaining ௧, the landslide displacement is obtained by a fully connected layer as the prediction output of the model. The Baishuihe landslide is located on the south bank of the main trunk road of the Yangtze River between Zigui County and Badong County, 56 km away from the three Gorges Dam sites, with a longitude of 1100320009<sup>00</sup> and latitude of 31001003400. The specificgeographical location of the Baishuihe landslide is shown in Figure 3.

The Baishuihe landslide is located on the south bank of the main trunk road of the Yangtze River between Zigui County and Badong County, 56 km away from the three Gorges Dam sites, with a longitude of 110′32″09″ and latitude of 31′01″34″. The specific

geographical location of the Baishuihe landslide is shown in Figure 3.

ℎ ←

**Figure 3.** The geographical location of Baishuihe landslide. **Figure 3.** The geographical location of Baishuihe landslide. **Figure** The geographical location of Baishuihe

The overall topography of the landslide area is high in the south and low in the north. The relative elevation difference of the terrain in the landslide area is approximately 300 m, the frontal elevation is approximately 70 m, the north–south length is 600 m, the east– west width is 700 m, and the average thickness of the slide body is approximately 30 m [6,10]. The longitudinal sliding surface of the Baishuihe landslide area is a folded line, which is steep at the back and shallow at the front, and the middle slide surface is between the two sliding surfaces. The Baishuihe landslide area has two slip zone layers, the shallow slide zone is the interface of gravel soil and cataclasite, the thickness is approximately 0.9~3.13 m, and the buried depth is 12.4~20.3 m [19,27]. The deep slide zone is the contact surface between cataclastic rock and carbonaceous silt mudstone, with a thickness of 0.6~1.5 m and burial depth of 18.9~34.1 m [31]. According to the classification of Hungr et al. [53], the Baishuihe landslide belongs to the clay/silt planar slide. Its deformation is slow. A topographical map of the Baishuihe landslide area is shown in Figure 4. The overall topography of the landslide area is high in the south and low in the north. The relative elevation difference of the terrain in the landslide area is approximately 300 m, the frontal elevation is approximately 70 m, the north–south length is 600 m, the east–west width is 700 m, and the average thickness of the slide body is approximately 30 m [6,10]. The longitudinal sliding surface of the Baishuihe landslide area is a folded line, which is steep at the back and shallow at the front, and the middle slide surface is between the two sliding surfaces. The Baishuihe landslide area has two slip zone layers, the shallow slide zone is the interface of gravel soil and cataclasite, the thickness is approximately 0.9~3.13 m, and the buried depth is 12.4~20.3 m [19,27]. The deep slide zone is the contact surface between cataclastic rock and carbonaceous silt mudstone, with a thickness of 0.6~1.5 m and burial depth of 18.9~34.1 m [31]. According to the classification of Hungr et al. [53], the Baishuihe landslide belongs to the clay/silt planar slide. Its deformation is slow. A topographical map of the Baishuihe landslide area is shown in Figure 4. The overall topography of the is high the south and low the north. The relative elevation difference of the terrain in the landslide area is approximately 300 m, the frontal elevation is approximately 70 m, the north–south length is 600 m, the east– west is 700 m, and the average thickness the slide is approximately 30 m [6,10]. sliding surface the Baishuihe area is a folded line, which is steep the back and at the and the middle slide is the two The Baishuihe landslide area has two slip zone layers, the shallow slide zone is the interface of gravel soil and cataclasite, the thickness is approximately 0.9~3.13 m, and the buried depth is 12.4~20.3 m [19,27]. The deep slide zone is the contact surface between cataclastic and silt mudstone, with a of 0.6~1.5 depth of 18.9~34.1 [31]. According to the of Hungr et al. [53], the landslide the planar slide. deformation slow. A topographical map of the Baishuihe landslide area is shown in Figure 4.

**Figure 4.** The GPS monitoring location of the Baishuihe landslide area. **Figure 4.** The GPS monitoring location of the Baishuihe landslide area.

Eleven global positioning system (GPS) monitoring sites have been installed in the Baishuihe landslide area. The GPSs are installed by the National Cryosphere Desert Data **Figure 4.** The GPS monitoring location of the Baishuihe landslide area. Eleven global positioning system (GPS) monitoring sites have been installed in the Baishuihe landslide area. The GPSs are installed by the National Cryosphere Desert Data Eleven global positioning system (GPS) monitoring sites have been installed in the Baishuihe landslide area. The GPSs are installed by the National Cryosphere Desert Data Center. The purpose is to observe the impact of the Three Gorges Dam on landslides

Center. The purpose is to observe the impact of the Three Gorges Dam on landslides in

Center. The purpose to the impact of the Gorges Dam landslides in

in this area. The ZG118 monitoring station is located in the center of the landslide area, and it can well reflect the whole situation of the Baishuihe landslide. It is used by most research institutes [11,19,54], so the data from the ZG118 monitoring station are also used as representatives of Baishuihe landslide data in this study. In this study, the displacement data of the Baishuihe landslide area for 108 months are used, and the variations in rainfall and reservoir water level within the landslide range in the same period are monitored. The data-collection time period was from January 2004 to December 2012, and one data point was collected every month, as shown in Figure 5. This dataset and the topographical map are provided by the National Cryosphere Desert Data Center/National Service Center for Speciality Environmental Observation Stations. this area. The ZG118 monitoring station is located in the center of the landslide area, and it can well reflect the whole situation of the Baishuihe landslide. It is used by most research institutes [11,19,54], so the data from the ZG118 monitoring station are also used as representatives of Baishuihe landslide data in this study. In this study, the displacement data of the Baishuihe landslide area for 108 months are used, and the variations in rainfall and reservoir water level within the landslide range in the same period are monitored. The data-collection time period was from January 2004 to December 2012, and one data point was collected every month, as shown in Figure 5. This dataset and the topographical map are provided by the National Cryosphere Desert Data Center/National Service Center for Speciality Environmental Observation Stations.

*Mathematics* **2022**, *10*, x FOR PEER REVIEW 9 of 20

**Figure 5.** The variations in the displacement, rainfall, and reservoir water level of the Baishuihe landslide area. **Figure 5.** The variations in the displacement, rainfall, and reservoir water level of the Baishuihe landslide area.

#### *3.2. Analysis of Influencing Factors 3.2. Analysis of Influencing Factors*

Figure 5 shows that from April to August every year, the rainy season of the Baishuihe landslide area occurs, and the amount of rainfall continues to rise. Every year, Baishuihe's landslide displacement increases most rapidly during the height of the rainy season, and the largest increase in landslide displacement occurred during the rainy season in 2007. The dry season in the Baishuihe landslide area is from January to March. There is less rainfall in these months, and the landslide displacement hardly increases. As rainfall increases, the amount of water entering the slope increases, making the overall weight of the slope rise, thus increasing the speed of the landslide. After rainwater enters the slope, the interior of the slope becomes wet, reducing the friction force and increasing the landslide probability. Rainfall is usually used as one of the inputs to predict landslide Figure 5 shows that from April to August every year, the rainy season of the Baishuihe landslide area occurs, and the amount of rainfall continues to rise. Every year, Baishuihe's landslide displacement increases most rapidly during the height of the rainy season, and the largest increase in landslide displacement occurred during the rainy season in 2007. The dry season in the Baishuihe landslide area is from January to March. There is less rainfall in these months, and the landslide displacement hardly increases. As rainfall increases, the amount of water entering the slope increases, making the overall weight of the slope rise, thus increasing the speed of the landslide. After rainwater enters the slope, the interior of the slope becomes wet, reducing the friction force and increasing the landslide probability. Rainfall is usually used as one of the inputs to predict landslide displacement models [55–58].

displacement models [55–58]. Although the largest increase in landslide displacement occurred during the rainy season of 2007, there were three years with more rainfall than that in 2007, suggesting that rainfall is not the only factor affecting landslide displacement. Because of the complexity of the internal geological structure of a landslide, landslides usually contain several different states, and landslide stability differs between such states. When a landslide is relatively stable, it is difficult to produce a large displacement when strong external factors are encountered. When a landslide is in an unstable state, a large displacement may occur even if external factors of general strength are encountered. Therefore, the maximum rain-Although the largest increase in landslide displacement occurred during the rainy season of 2007, there were three years with more rainfall than that in 2007, suggesting that rainfall is not the only factor affecting landslide displacement. Because of the complexity of the internal geological structure of a landslide, landslides usually contain several different states, and landslide stability differs between such states. When a landslide is relatively stable, it is difficult to produce a large displacement when strong external factors are encountered. When a landslide is in an unstable state, a large displacement may occur even if external factors of general strength are encountered. Therefore, the maximum rainfall is not the factor that corresponds to the maximum landslide displacement; this is consistent with the actual displacement in the Baishuihe landslide area. To some extent, the displacement before the landslide can represent the state of the landslide at that time and the stability of the landslide [14,17,28,59]. Therefore, we take the displacement of the landslide in the previous month as the state of the landslide as one of the inputs of the prediction model.

Since 2003, the Three Gorges Dam has released water during the rainy season to ensure that the dam is safe; subsequently, the reservoir level of the dam has dropped significantly. It can be observed from Figure 5 that the rate of landslide displacement increases at the end of the decrease in the reservoir water level every year, indicating that the decrease in the reservoir water level has a certain lag effect on landslide displacement. When the water level of the reservoir decreases for a certain period of time, the resistance of the landslide surface will be reduced. When the reservoir releases more water, the impact force of the landslide also increases with the increase in water flow, making the structure of the landslide more easily affected such that a landslide is more likely to occur. Therefore, the reservoir water level is also considered one of the influencing factors of landslide displacement [60–63].

In this paper, it is speculated that landslide displacement is the result of the comprehensive action of rainfall, reservoir level, and landslide state. Therefore, rainfall, reservoir water level, and landslide state are selected as the input for the prediction model in this paper.

#### *3.3. Decomposition Data Using the LMD Algorithm*

In this experiment, we used 108 months of data for simulation verification. The data of the first 96 months were used as a training set for the training model, and the data from the last 12 months were used as a test set for model prediction to verify the performance of the landslide displacement prediction model. According to the timefrequency analysis method, the LMD algorithm was used to decompose the landslide displacement, precipitation, reservoir water level, and landslide state from high frequency to low frequency into six subsequences, as shown in Figure 6.

The landslide displacement, rainfall, reservoir water level, and landslide state were decomposed by the LMD algorithm, and five *PF<sup>s</sup>* with different frequencies and one *u<sup>k</sup>* residual component were obtained.

#### *3.4. Calculating the Correlation between Landslide Displacement and Influencing Factors*

Due to the complexity of the landslide attributes, there are many factors affecting landslides, and different landslide states cause different landslide consequences when the same influence conditions are met. After the landslide displacement and influencing factor data are decomposed into subsequence data of different frequencies by the LMD algorithm, the subsequences of landslide displacement for each frequency correspond to multiple influencing factors of different frequencies. However, any influencing factors with too little correlation reduce the prediction performance of the model; additionally, using too many irrelevant factors for prediction leads to an increasing number of deviations in the results. Therefore, the MIC algorithm is adopted in this paper to calculate the correlation degree between each of the subsequences of the landslide displacement and the different frequency subsequences of the various influencing factors. Table 1 shows the MIC correlation results for each set of landslide displacement subsequence with other subsequences.

The selection of an appropriate MIC value is particularly critical for the selection of influencing factors and the final prediction results of the model. MIC values that are too small result in factors with low correlation participating in the training and prediction of the prediction model. A large range of MIC values leads to fewer datasets for the model to use. According to the results of the study [17] and several experiments, this paper selects the influential factors with MIC > 0.3 as the input of the prediction model. The final results show that there are 16, 13, 16, 17, 14, and 17 input variables for *PF*1, *PF*2, *PF*3, *PF*4, *PF*5, and *uk* in the landslide displacement subsequences, respectively.

**Figure 6.** Subsequence decomposition of the landslide displacement, precipitation, reservoir water level, and landslide state. **Figure 6.** Subsequence decomposition of the landslide displacement, precipitation, reservoir water level, and landslide state.


**Table 1.** MIC results for each landslide displacement subsequence with other subsequences.

#### *3.5. Prediction Using the BiLSTM Model*

After the MIC algorithm calculation, the input variables of each subsequence of landslide displacement are obtained. The forecasting process of the BiLSTM model is introduced with *PF*<sup>1</sup> as an example. *PF*<sup>1</sup> has 16 input variables, which are combined into two 16-dimensional vectors. A vector of the model's training set contains data for the first 96 months, and the other vector contains data for the next 12 months as a test set for the model. Training sets are used to train the BiLSTM model, and the model adjusts the parameters adaptively in the training process. After the training, the test set is input into the trained BiLSTM model to obtain the prediction results. The process of *PF*<sup>1</sup> prediction is similar to those of the other five subsequences, and the prediction results obtained are shown in Figure 7.

After the six components of landslide displacement are predicted, they are summed to obtain the final prediction result, as shown in Figure 8.

**Figure 7.** The components of the landslide displacement prediction results. **Figure 7.** The components of the landslide displacement prediction results.

**Figure 8.** Final prediction results of landslide displacement. **Figure 8.** Final prediction results of landslide displacement.

#### **4. Discussion**

and accuracy of the model.

sults are shown in Figure 9.

**4. Discussion**  Landslide displacement prediction is a typical regression problem, so to describe the prediction performance of the LMD-BiLSTM model more accurately, this study selects Landslide displacement prediction is a typical regression problem, so to describe the prediction performance of the LMD-BiLSTM model more accurately, this study selects four

four performance indicators to evaluate the prediction effects of various models. The mean absolute percentage error (MAPE), mean absolute error (MAE), root-mean-square

> 1 <sup>∑</sup> ୀଵ

> > <sup>∑</sup> ୀଵ

> > > ො −

(ො − )ଶ

(ത − )ଶ

100% <sup>∑</sup> ୀଵ |

> ∑ ୀଵ

> ∑ ୀଵ

where indicates the number of landslide data points, ො = {ොଵ, ොଶ,…,ො} is the predicted value of the model, = {ଵ, ଶ,…} is the actual value of the Baishuihe landslide, and ത = {തଵ, തଶ,…ത} is the average value of the actual Baishuihe landslide. The three evaluation indexes MAE, RMSE, and MAPE all indicate that the smaller the value is, the better the performance and accuracy of the model. R2 also evaluates the model according to the numerical value, and the higher the evaluation value is, the better the performance

To further verify the effectiveness and predictive performance of the LMD-BiLSTM model, this study uses the LMD-BiLSTM without MIC and the LMD-LSTM model to simulate 108 data points of the Baishuihe landslide simultaneously. Four models are used to simultaneously predict five components and components, and the prediction re-


(ො − )ଶ (20)


(22)

MAE =

RMSE = <sup>ඨ</sup><sup>1</sup>

R<sup>ଶ</sup> =1−

MAPE =

performance indicators to evaluate the prediction effects of various models. The mean absolute percentage error (MAPE), mean absolute error (MAE), root-mean-square error (RMSE), and determination coefficient R-squared (R<sup>2</sup> ) are calculated. Each evaluation index is defined as follows:

$$\text{MAE} = \frac{1}{n} \sum\_{i=1}^{n} |\hat{y}\_i - y\_i| \tag{19}$$

$$\text{RMSE} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} (\mathcal{y}\_i - y\_i)^2} \tag{20}$$

$$\text{MAPE} = \frac{100\%}{n} \sum\_{i=1}^{n} \left| \frac{\mathcal{Y}\_i - y\_i}{y\_i} \right| \tag{21}$$

$$\mathbb{R}^2 = 1 - \frac{\sum\_{i=1}^n (\hat{y}\_i - y\_i)^2}{\sum\_{i=1}^n (\overline{y}\_i - y\_i)^2} \tag{22}$$

where *n* indicates the number of landslide data points, *y*ˆ = {*y*ˆ1, *y*ˆ2, . . . , *y*ˆ*n*} is the predicted value of the model, *y* = {*y*1, *y*2, . . . *yn*} is the actual value of the Baishuihe landslide, and *y* = {*y*<sup>1</sup> , *y*<sup>2</sup> , . . . *y<sup>n</sup>* } is the average value of the actual Baishuihe landslide. The three evaluation indexes MAE, RMSE, and MAPE all indicate that the smaller the value is, the better the performance and accuracy of the model. R<sup>2</sup> also evaluates the model according to the numerical value, and the higher the evaluation value is, the better the performance and accuracy of the model.

To further verify the effectiveness and predictive performance of the LMD-BiLSTM model, this study uses the LMD-BiLSTM without MIC and the LMD-LSTM model to simulate 108 data points of the Baishuihe landslide simultaneously. Four models are used to simultaneously predict five *PF* components and *u<sup>k</sup>* components, and the prediction results are shown in Figure 9.

After obtaining five *PF* components and *u<sup>k</sup>* components, the final prediction result is obtained by adding the subsequence prediction results of these models. To make a better comparison, the BiLSTM model is added to compare the final landslide displacement prediction results. The comparison of the final landslide displacement prediction results is shown in Figure 10.

We can see from Figures 9 and 10 that the BiLSTM model can also be used to predict landslide displacement without the processing of the LMD algorithm and MIC method, but it can only roughly predict the overall trend, and the predicted fluctuations are large. The prediction result of the LMD-BiLSTM model is similar to the curve of the LMD-BiLSTM without the MIC model, but it can be observed that the prediction deviation of some points is too large, which is the result of selecting too many input factors with low correlation without the MIC algorithm. The prediction results of the LMD-LSTM model are not as accurate as those of the LMD-BiLSTM model either on the whole or at a certain point, due to the lack of information about the future during model training and prediction. The experimental results also verify this idea.

To intuitively compare the performance of these prediction models, this paper uses MAE, MAPE, RMSE, and R<sup>2</sup> to evaluate the models. The results are shown in Table 2.

**Figure 9.** Comparison of the subsequence prediction results for landslide displacement based on four models. **Figure 9.** Comparison of the subsequence prediction results for landslide displacement based on four models. *Mathematics* **2022**, *10*, x FOR PEER REVIEW 16 of 20

**Figure 10.** Final prediction results for the landslide displacement of four models. **Figure 10.** Final prediction results for the landslide displacement of four models.

We can see from Figures 9 and 10 that the BiLSTM model can also be used to predict landslide displacement without the processing of the LMD algorithm and MIC method, but it can only roughly predict the overall trend, and the predicted fluctuations are large. The prediction result of the LMD-BiLSTM model is similar to the curve of the LMD-BiLSTM without the MIC model, but it can be observed that the prediction deviation of

correlation without the MIC algorithm. The prediction results of the LMD-LSTM model are not as accurate as those of the LMD-BiLSTM model either on the whole or at a certain point, due to the lack of information about the future during model training and predic-

MAE, MAPE, RMSE, and R2 to evaluate the models. The results are shown in Table 2.

LMD-BiLSTM without MIC 17.061 0.749 21.156 0.853

**Table 3.** Accuracy of the predicted landslide displacement based on the four models.

**Error** 

To intuitively compare the performance of these prediction models, this paper uses

**Models MAE MAPE RMSE R2** LMD-BiLSTM 7.276 0.322 8.511 0.976

LMD-LSTM 15.147 0.670 18.215 0.891 BiLSTM 24.501 1.074 29.524 0.713

In addition to using the model prediction performance evaluation index to evaluate the prediction results of these models, this paper also uses the minimum error, maximum error, relative error, and mean absolute error to intuitively judge the prediction perfor-

LMD-BiLSTM 0.518 16.544 3.784 7.109 LMD-BiLSTM without MIC 1.386 39.061 9.016 17.062

LMD-LSTM 1.680 30.583 7.983 15.168 BiLSTM 1.987 49.078 12.861 24.489

**Maximum Error** 

**Relative Error** 

**Mean Absolute Error** 

**Table 2.** Comparison of four models for landslide displacement prediction.

mance of the model, and the results are shown in Table 3.

**Models Minimum** 

tion. The experimental results also verify this idea.


**Table 2.** Comparison of four models for landslide displacement prediction.

In addition to using the model prediction performance evaluation index to evaluate the prediction results of these models, this paper also uses the minimum error, maximum error, relative error, and mean absolute error to intuitively judge the prediction performance of the model, and the results are shown in Table 3.


**Table 3.** Accuracy of the predicted landslide displacement based on the four models.

Tables 2 and 3 show that the model processed by the LMD algorithm and MIC algorithm is superior to the pure neural network BiLSTM model in terms of the comprehensive prediction level at all time points, the stability of the whole prediction, and the fluctuation change at a single time point. This superiority is because the hybrid model exploits the advantages of the LMD algorithm, which is good at analyzing the characteristics of data signals and reasonably reflecting the time and frequency distribution of data in various spaces and scales, which are advantages of the use of the MIC algorithm, which can calculate the degree of nonlinear association between two variables and also use the advantages of the BiLSTM model, which is good at processing time series data. Therefore, the hybrid model has the advantages of these algorithms and models. Because LMD-BiLSTM without the MIC model does not have the MIC algorithm to calculate the correlation of influencing factors, its prediction performance is between the LMD-BiLSTM model and BiLSTM model. Although the BiLSTM model takes the information sharing each cycle into account in the prediction, its prediction performance degrades because it does not remove redundant influencing factors, with the final prediction result being inferior to LMD-LSTM.

#### **5. Conclusions**

Landslide displacement prediction has been studied for a long time but is still a challenging research topic. In order to solve the disadvantage of ignoring random displacement in most time analysis methods, this paper proposes an LMD-BiLSTM model based on the time-frequency analysis method for landslide displacement prediction. The LMD algorithm is used for the first time to deconstruct the nonlinear and nonstationary data series of landslide displacement and influencing factors into multiple subseries. To improve the prediction accuracy, the MIC algorithm is used to quantify the correlation between the subsequences of landslide displacement and the subsequences of factors affecting the displacement; moreover, the factors with greater correlation are selected as the input variables of the model. According to the constantly changing characteristics of landslide displacement, precipitation, and reservoir water level in multiple time periods, the BiLSTM model is used to predict the subsequence components of landslide displacement, and the

final landslide predicted displacement is obtained by adding the predicted results. The real landslide dataset of the Baishuihe landslide in the Three Gorges Reservoir area of China is used in the experiment, and excellent results are obtained. The final results show that the new LMD-BiLSTM model proposed in this study can predict landslide displacement smoothly and accurately. In the future, the LMD-BiLSTM model will be improved according to the different characteristics of each landslide and then popularized and applied to the displacement prediction of other landslides. In addition, the LMD-BiLSTM model could also be applied to other forecasting fields, such as rainfall prediction and power generation prediction, to assist decision-makers in continuously improving the process of making reasonable judgments. The decision makers can add the landslide displacement prediction results into the landslide warning system to make their own judgment according to the landslide displacement, and the scientific community can build on these findings and apply these methods to other areas of prediction. The results presented in this paper may not be applicable to the early warning of landslides, because prediction models require a certain amount of monitoring data, which is not suitable for early warning.

**Author Contributions:** Conceptualization, X.S., Y.J. and Z.L.; methodology, Z.L. and X.S.; formal analysis, Y.J. and X.S.; resources, Z.L., Y.J., W.L. and X.S.; writing—original draft preparation, Z.L., Y.J. and X.S.; writing—review and editing, Z.L., Y.J. and X.S.; funding acquisition, X.S. and W.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant No. 62161007, Grant No. 62061010, and Grant No. 61861008), Technology Major Project of Nanning Qingxiu District (Grant No. 2018001), Department of Science and Technology of Guangxi Zhuang Autonomous Region (Grant No. AB21196041, Grant No. AA20302022, Grant No. AA19182007, and Grant No. AA19254029), Natural Science Foundation of Guangxi Province of China (Grant No. 2019GXNSFBA245072 and Grant No. 2018GXNSFAA294054), and Young Teachers Promotion Project of Guangxi Universities (Grant No. 2022KY0182).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Restrictions apply to the availability of these data. Data were obtained from the National Cryosphere Desert Data Center/National Service Center for Speciality Environmental Observation Stations and are available from the http://data.casnw.net/portal/ (accessed on 19 January 2021) with the permission of the National Cryosphere Desert Data Center/National Service Center for Speciality Environmental Observation Stations.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

