Data-Driven Combination-Interval Prediction for Landslide Displacement Based on Copula and VMD-WOA-KELM Method

0. INTRODUCTION

Landslide is a common natural disaster which poses a huge threat to human life and property (Liu et al., 2016; Li et al., 2010). Recently, thousands of landslides have occurred in the Three Gorges Reservoir in China, which are triggered by the strong rainfall and massive change of the reservoir water level (Ma et al., 2017; Zhou et al., 2016; Yin et al., 2010). A bunch of studies have been carried out on landslide prediction and treatment (Huang et al., 2017; Intrieri et al., 2013; Li et al., 2012). Thereinto, predicting landslide displacement reasonably is significant and essential for early warning of such hazard.

The evolution of landslide is a nonlinear complex dynamic process, while the related deformation is usually affected by multiple factors, such as geological settings, external environment change and human activities, etc. Some scholars (Herrera et al., 2009; Corominas et al., 2005) have attempted to establish physical models by building the relationship between deformation and the physical parameters of geotechnical materials in some given area. Although the established model is physically straightforward, the values of the related physical parameters are usually hard to obtain and always change with the transformation of landslide mass, especially in some complex large landslides. Therefore, the so-called data-driven models have been put forward from the initial deterministic model to the disposition-time statistical analysis model, subsequently to the nonlinear machine learning model during recent years (Zhou et al., 2016; Lian et al., 2014; Yin and Yan, 1996; Saito, 1965). These models are mostly established based on the deformation monitoring data and the deformation-related impact factors. Compared with the physical model, the data required by the data-driven model is easier to obtain. Nevertheless, the frequently-used statistical analysis model, nonlinear regression or polynomial fitting methods have been verified to be hard to perform deformation prediction efficiently. More recently, a number of scholars have resorted to the signal decomposition techniques (for instance, the empirical mode decomposition (EMD), the ensemble empirical mode decomposition (EEMD), the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN)) (Li L Q et al., 2021; Deng et al., 2017; Lin et al., 2011), combining with machine learning methods (for instance, the support vector machine (SVM), the extreme learning machine (ELM), the artificial neural network (ANN)) (Cao et al., 2016; Du et al., 2013; Liu et al., 2016; Wu et al., 2016; Zhou et al., 2016; Lian et al., 2014), to deal with this problem. In this measure, the landslide displacement is decomposed into a series of components with different characteristics, whose predicted values obtained by the machine learning models are added together to get the total predicted value. Compared with the traditional statistical prediction method, this measure has been verified to have a higher prediction accuracy and stronger generalization ability. However, the above studies mostly focus on the deterministic point prediction, and did not consider how to evaluate the uncertainties exist in the prediction process and quantify the reliability of the predicted result. For a complex dynamic landslide system, there are a variety of uncertainties involving the geological settings, measurement technology and prediction model exist during the displacement prediction process (Phoon and Kulhawy, 1999a, b), which can be quantified by constructing a PI (Wan et al., 2014). In recent years, although interval prediction has been widely applied in different fields (Lins et al., 2015; Shrivastava and Panigrahi, 2013; Khosravi et al., 2010; Shrestha and Solomatine, 2006), there is only a little work focusing on the landslide displacement PI. For instance, Wang et al. (2019) used the PSO-ELM of two-node output layer to directly obtain the upper and lower limits of the PI. Ma et al. (2018) adopted the Bootstrap combining with machine learning method to calculate the system error and random error caused by the uncertainties during constructing the PI. However, the former has a relatively low stability and higher calculation cost, while the latter has a higher requirement for parameter selection. Once the parameters fall into the local optimum during the calculation process, the prediction accuracy will be greatly reduced.

In this paper, we propose a combination-interval prediction method based on Copula and VMD-WOA-KELM model to construct the PI. Firstly, the real monitoring displacement is decomposed by VMD to some displacement components, whose key impact factors are respectively determined by establishing their individual Copula model. The WOA-KELM model of each displacement component is then built using these key impact factors to conduct the dynamic point prediction. The point prediction error of each displacement component is calculated simultaneously. By comparing the error fitting degrees of the PDFs fitted by Gaussian distribution and Gaussian-t mixture distribution for each displacement component, whose optimum prior distribution form can be selected according to the highest fitting degree to construct the parametric-method-based PI. On the other hand, the optimal window width kernel density estimation is used to obtain the PDF of the prediction error for each displacement component, further to construct the non-parametric-method-based PI. Finally, we combine the PIs of each displacement component calculated by the above two methods through DE to fulfill the better performance. The PI of total displacement is finally obtained by superimposing the PI of each displacement component.

1. METHODOLOGY

1.1 VMD Decomposition of Landslide Displacement

The VMD is a new adaptive signal decomposition algorithm proposed by Dragomiretskiy and Zosso (2014). It can decompose the complex non-stationary signal into different modal component μ_k with its central frequency. It can not only effectively overcome the modal aliasing phenomenon exist in the EMD process, but also deal with the endpoint effect through mirroring continuation, and filter out the noise while retain the useful components in the original signal (Li et al., 2018; Xu and Niu, 2018).

The occurrence of landslides is usually accompanied by abrupt displacement variation. If the original displacement data sequence is directly used in displacement prediction, the trained model may misidentify the abrupt segment as abnormal data (Li et al., 2019) to increase difficulties of displacement prediction. Meantime, estimating error at the first sharp point according to the prediction error PDF in the earlier stable deformation stage will result a great result bias, and the poor point prediction effect at this point will lead to an error spurt, which may result in accumulated effect to the prediction error estimation for the following point and reduce reliability of the interval prediction model. We use the VMD to decompose displacement into a series of displacement component to extract the displacement local features. The prediction error PDF of each displacement component is relatively stable, and the work mainly focuses on how to find the specific error PDF for each component, which can be used to estimate their error at the next point effectively.

1.2 Copula Correlation Analysis

Selecting the effective input data can improve the prediction accuracy and give full play to the prediction performance of the model, which usually relies on evaluating the correlation between the input and output data correctly. The traditional linear correlation analysis has limitations for random variables that satisfy the non-normal distribution characteristic, which is difficult to deal with the variables with "thick tails". The coupla method gives a new insight to cope with this problem. Sklar (1996) suggested that for multiple marginal distribution functions ${\mathrm{F}}_{1}\left({x}_{1}\right), {\mathrm{F}}_{2}\left({x}_{2}\right), ..., {\mathrm{F}}_{\mathrm{N}}\left({x}_{N}\right)$, there is a copula function C that satisfies the following formula.

where, $\mathrm{F}({x}_{1}, {x}_{2}, ..., {x}_{N})$ is the joint distribution function of the marginal distribution ${\mathrm{F}}_{1}\left({x}_{1}\right), {\mathrm{F}}_{2}\left({x}_{2}\right), ..., {\mathrm{F}}_{\mathrm{N}}\left({x}_{N}\right)$.

Typical copula functions are mainly categorized into two groups (Reboredo, 2011). One is the elliptic copula function including Gaussian copula and Student t copula. The other is the Archimedes copula function including Gumbel copula, Clayton copula, and Frank copula. In practical application, the optimal copula function for target data is usually determined by Akaike information criterion (AIC) and Bayesian information criterion (BIC) (Gao, 2019; Posada and Buckley, 2004).

where $p{\mathrm{e}}_{i}$ is the empirical probability; ${p}_{i}$ is the theoretic frequency; m is the number of the model parameter; n is the number of data samples, and L is the maximum log-likelihood function. In this study, we choose the copula function with the minimum AIC and BIC values as the optimal function.

The correlation between the impact factors and displacement can be evaluated by the Kendall rank correlation coefficient τ and the upper tail correlation coefficient. When their values are larger than 0.6, there is a strong correlation between the two variables (Liao et al., 2019; Deng et al., 2017). The detailed calculation process can refer to the literature by Tang (2018).

1.3 Point Prediction Method WOA-KELM

According to the previous study (Sebbar et al., 2020), compared with the traditional neural network and SVM methods, the kernel extreme learning machine (KELM), which is an optimized method of the ELM, can effectively save the training time and avoid the overfitting as well as the local optimum problems through reducing the number of relative parameters. The ELM is a feedforward neural network training model with maximum likelihood of single hidden layer, which has great generalization ability and fast learning speed (Huang et al., 2006a, b; Huang, 2003). The basic network structure is shown in Figure 1.

Figure  1.  KELM network structure diagram

DownLoad: Full-Size Img PowerPoint

For sample $({x}_{i}, {y}_{i})$, ${x}_{i}=[{x}_{i1}, {x}_{i2}, ..., {x}_{im}{]}^{T}\in {R}^{m}$, ${y}_{i}\in R$, the output value of ELM is expressed as follows.

where N is the number of hidden node; ${\beta }_{i}=[{\beta }_{1}, {\beta }_{2}, ..., {\beta }_{N}]$ is the weight vector connecting the ith hidden node and the output nodes; ${\omega }_{i}=[{\omega }_{1i}, {\omega }_{2i}, ..., {\omega }_{mi}]$ is the weight vector connecting the ith hidden node and input nodes; ${b}_{i}=[{b}_{1}, {b}_{2}, ..., {b}_{2}]$ is the threshold of the ith hidden node; $h({\omega }_{i}\cdot {x}_{i}+{b}_{i})$ is the output function of the hidden node. The training goal of ELM is to calculate the output weight matrix $\beta$ that can achieve the best prediction effect. This training process can be fulfilled by least square linear regression using the hidden layer output matrix and real output data to implement. The least square solution $\stackrel{-}{\beta }$ is shown in Eq. (6).

where ${H}^{+}$ is the Moore-Penrose generalized inverse matrix of the hidden layer output matrix H.

To improve generalization ability and stability of ELM algorithm, Huang et al. (2012) introduced the kernel function into the ELM to form the KELM. Based on the orthogonal projection method and ridge regression theory, the value of $\beta$ can be calculated by adding the coefficient 1/C. The calculation formula is expressed as follows.

Therefore, the output function of ELM is expressed as follows.

The kernel function is used instead of H(x), and the output function of KELM is expressed as follows.

Herein, we select the most widely used radial basis function as the kernel function $K(x, {x}_{i})$ of the KELM algorithm.

where $\lambda$ is the kernel function parameter, and ${\lambda }_{0}$; ${‖x-{x}_{i}‖}^{2}$ is the square Euclidean distance between input feature vector $x$ and ${x}_{i}$.

In this study, we take the whale optimization algorithm (WOA) to optimize the penalty coefficient C and kernel function parameter $\lambda$, which has a great impact on the prediction performance to make the effect of point prediction better. The WOA is a novel nature-inspired meta-heuristic optimization algorithm, which mimics the social behavior of humpback whales and whose specific principles can refer to related literature (Mirjalili et al., 2014).

1.4 Prediction Interval Construction Method

There is a great deal of uncertainties that bring error to the landslide displacement point prediction, mainly including the inherent variability of the geotechnical mass, the monitoring system and the prediction model. In this study, we use a straightforward method to quantify overall effect of the three uncertainties. The landslide displacement prediction error ε is defined as the displacement deviation. Herein, $y$ represents the actual measured value and $\widehat{y}$ represents the predicted value of the landslide displacement.

The prediction error of historical data is used to estimate error PDF, and then the displacement PI under the 95% confidence level can be calculated. The PDF of prediction error is usually evaluated by PM or NPM (Yang et al., 2020). In PM, it is estimated by making an ensemble distribution hypothesis. However, the NPM can directly fit it according to the real error data. The former method has some demerits in fitting accuracy and the latter method is prone to reach local optimization. In this study, we construct a more flexible CI by combining the PIs based on the above two methods.

1.4.1   Parametric interval prediction method (PIPM)

In PIPM, the parameters of the hypothetical prior distribution model are determined by fitting the historical real error PDF, and the PI can be obtained according to the chosen optimal prior distribution. The Gaussian distribution is frequently used as prior distribution model under most cases (Yang et al., 2020), the PDF of which is as follows.

where $x$ is the number of sample points; $\mu$ stands for the sample mean; ${\sigma }^{2}$ presents sample variance.

However, some practice reveals that the error distribution of target data may not follow the Gaussian distribution, where some unfavorable features like peak, thick tail, skewness, and multi-peak usually can be observed. The t-distribution model can fit these features to some extent, whose PDF is expressed as

where v is the degrees of freedom; $\mu, \sigma$ are expectation and standard deviation respectively, and y obeys t distribution with ν degrees of freedom if $y=\frac{x-\mu }{\sigma }$.

The mixed Gaussian-t distribution method has integrated the advantage of the above two models. It has been verified to perform better than a single distribution model in terms of the previous study (Yang et al., 2019). The PDF of the Gaussian distribution and t-distribution can be mixed as

where ${x}_{i}$ is the ith sample point; ${f}_{i}\left({x}_{i}\right)$ is the PDF of the original distribution; ${\omega }_{i}$ is the proportion of the ith distribution in the mixed distribution.

The least square method is used to fit the mixed PDF of the prediction error for each displacement component. According to the principle of the minimum residual error between the fitted and the real PDF, the parameters of the optimal mixed distribution model of each displacement component error are determined.

Following the above methods, we match each displacement component with corresponding optimal mixed prior distribution, and subsequently calculate the cumulative probability distribution function (CDF) according to its PDF (see Eq. (16)).

The confidence interval at the moment of t + 1 is correspondingly calculated as follows.

where ${\widehat{Q}}_{M}(\cdot)$ is the inverse function of ${Q}_{M}\left(\varepsilon \right)$; M represents the sequential order, assigning a value of 1, 2, 3, 4;$\alpha$ is the confidence level, ${\alpha }_{1}=\alpha /2$, ${\alpha }_{2}=1-\alpha /2$.

1.4.2   Non-parametric interval prediction method (NPIPM)

In the NPIPM, the PDF is directly fitted based on the real error data and avoids assuming the prior error distribution by means of the kernel density estimation (Yang et al., 2020). In this method, the accuracy is mainly determined by a moving window, of which a small width will lead to extremely unstable distribution result while a large one will make some information loss.

If we define the PDF of the random variable x as f(x), it can be obtained according to the kernel density estimation.

where $n\to +\mathrm{\infty }$, $h\to 0$, $n\to +\mathrm{\infty }$and $nh\to +\mathrm{\infty }$; n is the total number of samples; h is the window width; ${X}_{i}$ is the given sample; $K(\cdot)$ is the kernel function which is assigned the Gaussian function whose expression is shown in Eq. (19). The value of the coefficient $h$ is determined by an improved interpolation method, as the measure of the previous study (Ye et al., 2017).

where ${k}_{2}=\int {t}^{2}K\left(t\right)\mathrm{d}t$. For the given sample data, it is divided into two subsets, one of which is used to estimate the CDF of displacement prediction error, and the other of which is used to get the confidence interval that satisfies the confidence level of (1 $-$ α).

where ${\widehat{L}}_{2}$ and ${\widehat{U}}_{2}$ are the lower limit and upper limit of the confidence interval, respectively; two coefficients ${\alpha }_{1}$ and ${\alpha }_{2}$ can be calculated as ${\alpha }_{1}=\alpha /2$, ${\alpha }_{2}=1-\alpha /2$; $\widehat{F}(\cdot)$ is the inverse function of $F(\cdot)$, which is the integration of $f\left(x\right)$. In this study, the indicator x is the error between the predicted displacement and the actual displacement.

1.4.3   Combination interval prediction method (CIPM) based on the PIPM and NPIPM

In this section, the CIPM based on the PIs of PIPM and NPIPM is built. Thereinto, the choice of the combination weight coefficient plays a decisive role in the model performance. The DE algorithm, a heuristic random search algorithm based on population difference, is introduced because its good robustness and low calculation cost (Wang, 2018). Following the principle that the PI should realize the highest coverage while its width is the narrowest, the DE algorithm is used to combine the PIs based on the PIPM and NPIPM, and the following results are obtained.

where $\lambda$ and $(1-\lambda)$ are the combined weights of displacement PI based on PIPM and NPIPM, respectively. We take the interval evaluation comprehensive index CWC as fitness to optimize $\lambda$ to maximize the interval prediction performance.

1.5 Evaluation Index of Point Prediction and Interval Prediction

The prediction results of the model should be evaluated from two aspects. Firstly, the point prediction usually refers to the traditional evaluation indicators such as the degree of fitting (R²), the root mean square error (RMSE) and the mean absolute percentage error (MAPE) (Xiang, 2012). Secondly, three indicators including the prediction interval coverage probability (PICP), prediction interval normalized average width (PINAW), and coverage width-based criterion (CWC) are selected to evaluate the performance of the PI under a given confidence level (1 - α) (Wang et al., 2020; Shrivastava et al., 2016; Khosravi et al., 2011).

where ${N}_{test}$ is the sample size of test data set; $\xi$ is the number of actual displacement values that fall within the PI for the confidence level.

where R is the difference between the maximum and the minimum value in the predicted samples, and the index is normalized. $\widehat{U}\left({x}_{i}\right)$ is the upper limit of prediction. $\widehat{L}\left({x}_{i}\right)$ is the lower limit of prediction.

where $\gamma$ is Boolean variable whose value equals to zero when $PICP\mu$, which means eliminating the exponential terms. On the contrary, the setting for $\gamma =1$ means the exponential term remained. To enlarge the difference of PICP and the confidence interval$\mu$, the control parameter $\eta$ generally takes a larger value, which is 30 referring to a previous study (Khosravi et al., 2011).

1.6 Interval Prediction Process

The overall framework of PIs construction based on the above methods is shown in Figure 2.

Figure  2.  The overall framework of interval prediction.

DownLoad: Full-Size Img PowerPoint

2. CASE STUDY OF THE BAIJIABAO LANDSLIDE

2.1 Geological Settings

The Baijiabao landslide is located at the right bank of Xiangxi River in Zigui County, Hubei Province (see Figure 3), is a typical reservior bank deposit landslide (Li C D et al., 2021). The rear edge of the landslide is bounded by bedrock with an elevation of 275 m, and the front edge directly reaches the Xiangxi River. The elevation of the shear outlet is between 125 and 135 m. The left side is bounded by the bedrock and the right side is bounded by the mountain ridge. The landslide shows morphology of dustpan with a front edge width of 500 m, a rear edge width of 300 m, a longitudinal length of about 550 m and an area of 2.2 × 10⁵ m². The average thickness of the sliding body is 45 m. The overall terrain is high in the west with an inclination of 35°–50° and low in the east with an inclination of 0°–25°. The landslide mass is mainly Quaternary deposits composed of silty clay with crushed stone and crushed stone, and the sliding zone is mainly composed of silty clay, the thickness of which is generally 0.2–0.3 m; the sliding bed is composed of Lower Jurassic arkose sandstone, silty mudstone and argillaceous siltstone with dip of 260°–285° and dip angle of 30°–40°. The plan view and the engineering geological profile of the landslide are shown in Figure 4 and Figure 5. Four GPS stations including ZG323 to ZG326 were positioned on the landslide surface for periodic displacement monitoring.

Figure  3.  Location of the Baijiabao landslide in Three Gorges Reservoir area of China.

DownLoad: Full-Size Img PowerPoint

Figure  4.  Plan view of the Baijiabao landslide.

DownLoad: Full-Size Img PowerPoint

Figure  5.  Engineering geological profile of Baijiabao landslide.

DownLoad: Full-Size Img PowerPoint

2.2 Deformation Characteristic Analysis

The landslide cumulative displacement, reservoir water level fluctuation and rainfall in the reservoir area from January 2007 to December 2019 are illustrated in Figure 6. To be more clearly, the detailed data from January 2010 to December 2012 are listed in the Table 1. Referring to the previous study (Miao et al., 2018), the evolution models of landslides can be divided into four categories including steady-type, exponential-type, step-like-type and convergent-type, while the Baijiabao landslide belongs to the step-like type. Like most landslides in the Three Gorges area, the step-like deformation of Baijiabao landslide is mainly affected by the internal geological factors, the fluctuation of reservoir water level and the heavy rainfall (Li L Q et al., 2021; Yao W M et al., 2019; Yao W et al., 2015; Du et al., 2013). It can be seen that during the period of March to July in each year, the reservoir water level drops from the highest level to the lowest level, while the large accumulative displacement changes usually presented from May. From August to the following April, although the reservoir water level has increased significantly, the variation of cumulative displacement tends to be gentle. It indicates that the landslide deformation is positively correlated with the decrease of reservoir water level. Since the water level starts to decrease in March every year, while the displacement begins to increase significantly in May every year. Due to the two-month time lag here, we cannot directly identify the water level change of which month is closely related to the sharp landslide deformation. Therefore, the decrease of reservoir water level (DRWL) in the current month K1, the DRWL in the past one month K2 and the DRWL in the past two months K3 are usually taken as the primary impact factors when identifying the key impact factors. Moreover, the largest landslide displacement change usually occurs in July, while the heavy rainfall in this region appears from May to September in each year. The rapid movement period ends before the end of rainfall reason, indicating that the change of landslide deformation is less sensitive to the heavy rainfall than the decrease of reservoir water level. Similarly, the rainy season begins from May in each year while the largest displacement changes occurred in July. It is also difficult to determine the rainfall in which month is greatly contributed to the displacement change, so we usually choose the accumulative rainfall (AR) in the current month J1, the AR in the past one month J2 and the AR in the past two months J3 as the other primary impact factors.

Figure  6.  Landslide cumulative displacement, rainfall and reservoir water level from January 2007 to December 2019.

DownLoad: Full-Size Img PowerPoint

Table  1.  The detailed data of landslide cumulative displacement, rainfall and reservoir water level from January 2011 to December 2012

Monitoring date	Rainfall (mm)	Reservoir water level (m)	ZG323 (m)	ZG324 (m)	ZG325 (m)	ZG326 (m)
2011/1/1	12.20	174.545 45	368.200 8	414.062 5	401.673 6	510.460 3
2011/2/1	24.39	174.710 74	359.832 6	414.062 5	401.673 6	510.460 3
2011/3/1	36.59	166.280 99	359.832 6	414.062 5	401.673 6	510.460 3
2011/4/1	95.12	157.851 24	368.200 8	417.968 8	401.673 6	510.460 3
2011/5/1	182.93	154.380 17	368.200 8	417.968 8	401.673 6	518.828 5
2011/6/1	73.17	145.289 26	376.569 0	421.875 0	401.673 6	518.828 5
2011/7/1	185.37	149.586 78	460.251 1	476.562 5	401.673 6	644.351 5
2011/8/1	70.73	145.289 26	476.987 5	527.343 8	502.092 1	661.087 9
2011/9/1	114.63	153.719 01	485.355 7	546.875 0	527.196 7	669.456 1
2011/10/1	114.63	145.785 12	476.987 5	550.781 3	527.196 7	669.456 1
2011/11/1	109.76	167.603 31	476.987 5	550.781 3	527.196 7	669.456 1
2011/12/1	14.63	174.876 03	485.355 7	542.968 8	527.196 7	669.456 1
2012/1/1	7.32	174.710 74	485.355 7	539.062 5	527.196 7	669.456 1
2012/2/1	4.88	171.570 25	476.987 5	542.968 8	527.196 7	669.456 1
2012/3/1	48.78	164.462 81	476.987 5	546.875 0	527.196 7	669.456 1
2012/4/1	53.66	163.471 00	476.987 5	546.875 0	527.196 7	669.456 1
2012/5/1	114.63	162.809 92	485.355 7	550.781 3	527.196 7	677.824 3
2012/6/1	163.41	145.223 97	518.828 5	574.218 8	535.564 9	702.928 9
2012/7/1	31.71	158.843 00	644.351 5	734.375 0	552.301 3	870.292 9
2012/8/1	112.20	162.975 21	635.983 3	734.375 0	702.928 9	887.029 3
2012/9/1	60.98	146.280 99	652.719 7	742.187 5	702.928 9	887.029 3
2012/10/1	31.71	166.115 70	661.087 9	746.093 8	711.297 1	903.765 7
2012/11/1	0.00	173.388 43	661.087 9	746.093 8	711.297 1	912.133 9
2012/12/1	4.88	174.545 45	652.719 7	750.000 0	711.297 1	903.765 7

 | Show Table

DownLoad: CSV

2.3 Decomposition of Cumulative Displacement

In recent years, the displacement of ZG326 varies more significantly than that of other points. Taking the monitoring displacement of this point as base data, we establish a VMD model by repeat trails to determine the optimal values of model relevant parameters, like the mode number K of 4, penalty factor $\alpha$ of 2 000 and time step of 0.2 referring to the measure of previous literature (Li et al., 2018). On this basis, the original displacement of ZG326 is decomposed into three IMF components and one residual component with different frequencies, as shown in Figure 7. The four displacement components have their own corresponding physical meaning. Referring to the study by Zhang et al. (2021), we also use the K-means clustering theory to analyze the features of these IMFs. It is found that the IMF2 and IMF3 components can be classified as one group as the periodic displacement, and the residual component can be treated as the trend displacement while the IMF1 represents the random displacement. This measure has an advantage of describing the periodic displacement more detailed.

Figure  7.  VMD decomposition results of displacement at ZG326 from 2007 to 2019.

DownLoad: Full-Size Img PowerPoint

The 156 groups of displacement data at ZG326 are divided into training data set and back-up data set by a ratio of 8 : 4. In the back-up set, the last 12 groups of data are selected as a testing set, and the rest data are used to generate the prediction error. The back-up data set is predicted for only one point stepwise. The monitoring data at the current moment is incorporated into the training set before the data at the next moment being predicted, which makes the training set be updated in real time. The updated training data set is learned by the WOA-KELM model to regain the values of the optimal relevant parameters. The above procedure is repeated at each prediction moment. With more and more real monitoring data being obtained, the training set is extended by the new monitoring data sequentially, making the prediction accuracy constantly increased.

2.4 Selection of Key Impact Factors Based on Copula Model

Referring to previous studies (Cao et al., 2016; Zhou et al., 2016), the trend displacement over the past 1, 2, and 3 months should be input into the WOA-KELM model to predict the trend displacement of the current month. Besides, for the random and periodic displacements, it can be known from the section of "Deformation characteristic analysis" that both the heavy rainfall and the decrease of the reservoir water level are the primary impact factors triggering the displacement fluctuation. It is hard to directly determine the key impact factors for each IMF displacement component. Therefore, we build Copula model between each IMF component and six impact factors, including K1, K2, K3, J1, J2 and J3, in which the above-mentioned five typical copula functions are used to get individual joint-distribution function for different marginal distributions. The values of the parameters in each marginal distribution are estimated using the maximum likelihood estimation method, and then the copula function with minimum values of AIC and BIC is selected for each copula model. It is found that the Gumbel copula function performs best for each model, maintaining a strong link between the impact factors and the IMF displacement component at the upper tip of the joint distribution, while being asymptotically independent at the lower tip. This process is carried out on the platform of MATLAB R2020b running on a Core I7-4700MQ @2.4 GHz CPU with 32-GB RAM. The density function of the joint distribution of K2-IMF3 and J2-IMF3 based on Gumbel Copula function is individually shown in Figures 8a and 8b. The correlation coefficients of the models of K1-IMF3, K2-IMF3, K3-IMF3 and J1-IMF3, J2-IMF3, J3-IMF3 returned by this function are shown in Table 2. Note that although some researchers have resorted to the traditional grey correlation degree to conduct correlation analysis (Miao et al., 2017), it can hardly deal with the thick tail feature between the impact factors and landslide deformation. The Copula theory takes inherent strengths to tackle this issue.

Figure  8.  The density function between the impact factor and IMF3 displacement component.

DownLoad: Full-Size Img PowerPoint

Table  2.  The correlation coefficient of Gumbel copula model

Correlation coefficient	K1-IMF3	K2-IMF3	K3-IMF3	J1-IMF3	J2-IMF3	J3-IMF3
Kendall rank correlation coefficient τ	0.603 4	0.738 7	0.694 4	0.696 4	0.653 6	0.608 5
Upper-tail correlation coefficient	0.747 5	0.883 4	0.855 3	0.754 6	0.775 6	0.624 6

 | Show Table

DownLoad: CSV

In Table 2, the kendall rank correlation coefficient $\tau$ is used to describe in quantity whether there is harmonious between the impact factor and the IMF displacement component. The upper-tail correlation coefficient is used to weigh if there is notable influence of a larger independent variable on the related variable. It can be seen that there is a harmony and close relationship exist between the six individual impact factors and the IMF3. For the reservoir water level, the related coefficient values for both K2-IMF3 and K3-IMF3 are larger than those for K1-IMF3. On the other hand, the coefficient values of J1-IMF3 and J2-IMF3 are larger than J3-IMF3. It suggests that the factors K2, K3 and J1, J2 have greater influence on IMF3, despite the close relationships exist between all the given six impact factors and IMF3. The same phenomenon can be also observed in the relationship between those impact factors and the other two IMF displacement components. Therefore, it is reasonable to adopt K2, K3 and J1, J2 as the input data set to establish the following WOA-KELM model.

2.5 Point Prediction of Displacement

The WOA-KELM model for each displacement component based on the training data set is established respectively, and then is used to predict the back-up data set. Referring to previous studies (Zhang et al., 2021; Mirjalili et al., 2014), for the WOA in each model, we set the number of variables as 2, the number of search agents as 30 and the maximum iteration number as 200. The upper limit and lower limit of the variable vector to be optimized (C, $\lambda$) are respectively set to (1 000, 10) and (10, 0). The predicted value and the true value are cross-validated to obtain the fitness value. The optimal C and $\lambda$ of each model are obtained by WOA optimization, as shown in Table 3. By introducing the optimal C and $\lambda$ into each KELM model, the predicted displacement of each displacement component was obtained.

Table  3.  The optimal parameters value of KELM for each IMF displacement component

KELM parameter	IMF1	IMF2	IMF3	R
C	433.893 2	284.894 8	165.438 0	40.538 5
λ	6.237 4	5.850 3	3.582 2	0.785 0

 | Show Table

DownLoad: CSV

The predicted cumulative displacement is obtained by superimposing each predicted displacement component. To verify the effectiveness and superiority of the constructed WOA-KELM model, the same data are adopted to the WOA-ELM, WOA-SVM and WOA-BPNN to conduct the prediction work. The evaluation indexes values of prediction accuracy are shown in Table 4. It indicates that the WOA-KELM achieves the best prediction performance among the four models with the R², RMSE and MAPE values of 0.946, 17.328 7, and 0.327 2, respectively. It also reflects that the WOA-KELM has stronger generalization and greater learning ability for high-dimensional nonlinear data.

Table  4.  Comparison of prediction results of different models

Evaluation index	WOA-KELM	WOA-ELM	WOA-SVM	WOA-BPNN
R2	0.946	0.915	0.889	0.864
RMSE	17.328 7	21.594 2	27.432 4	29.786 2
MAPE	0.327 2	0.398 7	0.483 5	0.608 9

 | Show Table

DownLoad: CSV

2.6 Interval Prediction of Displacement

There may be an underfitting phenomenon existing during the process of generating the error data set for each displacement component because only 40 groups of data can be used in this study. It may introduce difficulties in choosing a reasonable prior distribution that reflects the real error distribution. To tackle this problem, we choose the displacement error data at other three monitoring points as the complementary data for that at the monitoring point ZG326, because these four locations share the same internal and external conditions.

The Gaussian distribution and Gaussian-t mixed distribution are used to fit the PDF of each displacement component error which has been supplemented by the error of other three points using the least square method. The parameters of the prior distribution model are determined to optimize the fitting degree for each displacement component error. The fitting results are shown in Figure 9. The curves fitted by the Gaussian-t mixed distribution are closer to the true distribution than those fitted by the Gaussian distribution for all four different displacement components. Therefore, we choose the former distribution to conduct the following PIPM work.

Figure  9.  Fitting curve of prediction error probability density. (a) IMF1 component prediction error probability density fitting; (b) IMF2 component prediction error probability density fitting; (c) IMF3 component prediction error probability density fitting; (d) R component prediction error probability density fitting.

DownLoad: Full-Size Img PowerPoint

Before predicting the displacement interval at a moment, the real error at the last moment should be added into the error generation set. Simultaneously, the parameters of the Gaussian-t mixed distribution model and its PDF should be re-determined according to the updated real error distribution. By repeatedly using the method proposed at Section 1.4.1, each displacement component PI under the 95% confidence level can be obtained, then be added up to attain the total displacement PI based on the PIPM, as shown in Figure 10a.

Figure  10.  Prediction interval of different methods. (a) Prediction interval of PIPM; (b) prediction interval of NPIPM; (c) prediction interval of CIPM.

DownLoad: Full-Size Img PowerPoint

On the other hand, according to the real data in the error generation set, the PDF of each displacement component prediction error is estimated. The estimated result is used to calculate the PI based on the NPIPM at each moment in the test set. Before the PI at a moment is estimated, the actual error data at the last moment is included in the error generation set, and the real PDF of the error is simultaneously updated. The PI of the NPIPM is shown in Figure 10b. Moreover, the population size, the magnification factor and the maximum iteration number of the DE is set to 40, 1 and 100, respectively. The optimal weight coefficient λ is obtained as 0.448 after 100 iterations. According to the Eq. (22), the PIs constructed according to the PIPM and the NPIPM are combined to obtain the CI of the displacement component (see Figure 10c).

From the Figure 10, it can be seen that although the PI constructed by the PIPM has a high coverage rate, the bandwidth is relatively large and accuracy is relatively low. The bandwidth of the PI constructed by the NPIPM is obviously narrower, but its coverage rate for the actual displacement is lower, meaning a poor prediction reliability. Neither of the above two PIs can guarantee the prediction accuracy while achieve a high prediction reliability. The PI constructed by the CIPM has higher accuracy than that constructed by the PIPM and higher prediction reliability than that constructed by the NPIPM, which presents that the CIPM is the most effective means to predict the displacement interval of Baijiabao landslide among the above three methods.

3. DISCUSSION

3.1 Model Performance Evaluation

In this section, we use two measures to evaluate the prediction performance of the CIPM. One is to compare the PIs obtained by different interval construction methods mentioned in Section 2.6, and the other is to compare the results acquired by the CIPM with those based on other methods mentioned in the previous studies.

For the former measure, the index values of the PICP, PINAW, and CWC based on the three methods are extracted in Figure 11. Although the PI constructed by the CIPM exhibits a medium performance in the PICP and PINAW with the values of 91.67 and 16.21, the CWC value is the minimum for 56.55, which 15.57 smaller than the PI constructed by the PIPM and 303.7 smaller than that constructed by the NPIPM. That reflects the superiority of the CIPM for interval prediction of landslide displacement.

Figure  11.  The interval evaluation index of different prediction methods.

DownLoad: Full-Size Img PowerPoint

For the latter measure, we apply the CIPM to the cases of the Shuping landslide and Baishuihe landslide reported in the previous studies by Wang et al. (2019) and Ma et al. (2018). The prediction performance of three methods is compared in Table 5. The PI obtained by the CIPM shares the same PICP value with those obtained by DES-PSO-ELM and the Bootstrap-ELM-ANN method for both landslide cases, indicating that the three methods are able to cover the same enough information of real displacement data. Nonetheless, the PIs obtained by the CIPM exhibit a smaller bandwidth than those obtained by the other two methods, which leads to the most excellent prediction performance for the corresponding CWC is the minimum.

Table  5.  Performance comparison of PI constructed by CIPM and other methods

Case		Shuping landslide			Baishuihe landslide
Method		CIPM	DES-PSO-ELM	Bootstrap-ELM-ANN	CIPM	DES-PSO-ELM	Bootstrap-ELM-ANN
Evaluation index	PICP	0.916 7	0.916 7	0.916 7	0.916 7	0.916 7	0.916 7
	PINAW	30.984 0	42.038 5	38.956 4	29.551 7	38.006 9	39.430 6
	CWC	108.114 3	146.687 4	135.932 8	103.116 5	132.619 7	137.493 3

 | Show Table

DownLoad: CSV

3.2 The Limitations and Further Work of CIPM

The Baijiabao landslide is a typical deposit landslide in the Three Gorges Reservoir area with the notable step-like deformation behavoir. The CIPM is proposed to make displacement interval prediction for this kind of landslide, while its feasibility for other types of landslides need to be verified in further study.

Due to that the scale of data obtained by GPS is relatively small in this study, the error generation set may be insufficient to fit an accurate PDF. As the measure we have used in the previous context, the error data of the near monitoring point in the same landslide can be added as the complementary to those of the target monitoring point ZG326. Subsequently, the increasing real monitoring data will be accessed with the time elapsing, and the magnitude of the error generation set will finally satisfy the fitting requirements for the PDF. Meanwhile, the monitoring data of the target point can be accessed directly to fit the PDF. Our method has some flexibility on different conditions that the existing data is relatively small or big, thus presenting a wide applicability.

In this study, the displacement detected by GPS is used as the data set to calculate the prediction error. Although the GPS monitoring technology have been greatly improved compared with the traditional manual monitoring means, it can still have some monitoring limitation. It is highly recommended that the multi-source data from GPS, Insar and TDR technology are integrated to excavate more real information for the target slope.

Another issue is that the known displacement is assumed to be unknown and used as the test data set during the process of prediction. Note that this measure is necessary to give a way to verify and evaluate the validity of the established prediction model and provide a reliable support for the follow-up real unknown displacement prediction. In practical application, we just need to construct the unknown displacement PI relying on the pre-monitored real displacement data and relative impact factor data obtained from the future rainfall forecast and the regulation plan for the future reservoir water level, which can provide the engineers with a valid reference.

4. CONCLUSION

Interval prediction of landslide displacement, which overcomes the limit of point prediction and quantities the uncertainties existing in the process of displacement prediction, is a significant component of landslide early warning systems. This paper proposes a novel CIPM based on Copula and VMD-WOA-KELM using landslide data from Baijiabao landslide. Firstly, the major impact factors for different displacement IMF component is determined by Copula model, and their nonlinear relationships are built via the WOA-KELM. Through this process, the point prediction model is constructed with a greatly improved performance. Secondly, the prediction error is obtained to estimate the error PDF by a DE algorithm in combination with PM and NPM. The final displacement CI is formed by adding a bunch of IMF PIs.

To verify the superiority of the suggested CIPM, the CI is used to compare with the PI constructed by the PIPM and NPIPM, which indicates that the CI can not only satisfy the prediction reliability but also take consider of the prediction accuracy. Subsequently, the prediction results acquired by our method are compared with those from other two interval prediction method for Shuping landslide and Baishuihe landslide cases. The result shows that the CIPM preforms best and is robust to construct the landslide displacement PI. This method can provide a good reference for the assessment and mitigation of landslide disasters.