Forecasting of Landslide Displacement Using a Probability-Scheme Combination Ensemble Prediction Technique

Data-driven models have been extensively employed in landslide displacement prediction. However, predictive uncertainty, which consists of input uncertainty, parameter uncertainty, and model uncertainty, is usually disregarded in deterministic data-driven modeling, and point estimates are separately presented. In this study, a probability-scheme combination ensemble prediction that employs quantile regression neural networks and kernel density estimation (QRNNs-KDE) is proposed for robust and accurate prediction and uncertainty quantification of landslide displacement. In the ensemble model, QRNNs serve as base learning algorithms to generate multiple base learners. Final ensemble prediction is obtained by integration of all base learners through a probability combination scheme based on KDE. The Fanjiaping landslide in the Three Gorges Reservoir area (TGRA) was selected as a case study to explore the performance of the ensemble prediction. Based on long-term (2006–2018) and near real-time monitoring data, a comprehensive analysis of the deformation characteristics was conducted for fully understanding the triggering factors. The experimental results indicate that the QRNNs-KDE approach can perform predictions with perfect performance and outperform the traditional backpropagation (BP), radial basis function (RBF), extreme learning machine (ELM), support vector machine (SVM) methods, bootstrap-extreme learning machine-artificial neural network (bootstrap-ELM-ANN), and Copula-kernel-based support vector machine quantile regression (Copula-KSVMQR). The proposed QRNNs-KDE approach has significant potential in medium-term to long-term horizon forecasting and quantification of uncertainty.


Introduction
As one of the most common natural hazards in the world, landslides pose a significant threat to public health and safety. According to statistics, landslides have affected 4.8 million people and caused 18,275 deaths during the period of 2009-2019 [1]. Landslide displacement prediction, which provides the necessary information to determine the extent of ongoing hazard, has proven to be the most cost-saving risk reduction measure [2][3][4]. However, landslide displacement prediction is complex and remains a key challenge in natural hazard research. This challenge arises because landslides are In this study, a probability-scheme combination ensemble prediction that employs quantile regression neural networks and kernel density estimation (QRNNs-KDE) was proposed for robust and accurate prediction and uncertainty quantification of landslide displacement. The Fanjiaping landslide with long-term and near real-time monitoring data was selected as a case study to explore the performance of the QRNNs-KDE approach. The deformation characteristics were clarified for fully understanding the triggering factors.

Description of Uncertainty Sources
Predictive uncertainty in data-driven models consists primarily of input uncertainty, parameter uncertainty, and model uncertainty [30][31][32].
The input uncertainty is related to the input data uncertainty and the input variable section uncertainty. The input data uncertainty is primarily due to measurement and sampling error and environmental noise. The input variable section uncertainty accounts for uncertainty inherent in the selection of input variables from the candidate data set. For physical models, the required inputs are pre-determined, being consistent with considered rheological models. However, for data-driven models, the selection of input variables is problem-dependent and cannot be determined in advance. Only major and relevant variables are selected as final inputs to train the data-driven model. The selection of the variables to include in a data-driven model from the original data set is inherently uncertain, especially when the input candidate pool is very large. For example, in data-driven models that utilize decomposition algorithms, only a portion of the decomposed sub-components are selected as input variables. The candidate input pool, which consists of sub-components, increases very quickly with the decomposition level and potentially increases the input variable selection uncertainty.
The parameter uncertainty refers to the uncertainty in the model parameter vector and mainly arises from the inability to identify a unique set of best parameters for the model [33].
Model uncertainty arises primarily from the model structure uncertainty and model error. Model structure uncertainty is associated with the specific model setting of learning algorithms, such as the polynomial order in polynomial regression models, the number of hidden nodes in an ANN or ELM, and the type of kernel function in an SVM. The input uncertainty may also account for model structure uncertainty, because different input variables "automatically" produce different model structures. Model error refers to the difference between two model estimates with respect to the corresponding target and is caused by the inability to reproduce the real processes.

Ensemble Prediction
Ensemble prediction is not a specific learning algorithm but a strategic combination of multiple predictions into a single output with a model combination process [21]. Based on the selection of the learning algorithm, ensemble prediction models can be further classified into homogeneous and heterogeneous ensemble models (Figure 1). A homogeneous ensemble model generates multiple learners with the same learning algorithm on different training datasets, which are produced by manipulating the original training data (schematic illustrated in Figure 1a). Bootstrap aggregation, also known as bagging for short, is the most straightforward and widely used method of manipulating the training dataset. By contrast, a heterogeneous ensemble model generates multiple learners with different learning algorithms on the same training data set (schematic illustrated in Figure 1b).

Quantile regression
Quantile regression is a common statistical technique for conducting inferences concerning conditional quantile functions [34,35]. More formally, any real-valued random variable Y may be characterized by its distribution function as follows: whereas for any 0 1 τ < < , Given a data set ( ( ), ( )) i x t Y t for 1, 2, , i I =  and 1, 2, , t N =  , the linear quantile regression can be expressed as follows: where 0 1 τ < < is the quantile, and b is an error with zero expectation.
The estimated parameters i θ can be approximated by minimizing a sum of the asymmetrically weighted absolute residual cost functions, which are expressed as follows: where ( ) Y t is the observation at time t and τ ρ is the check function, which is also known as the pinball loss function and is defined as follows: The base learner combination is the main step in the ensemble prediction model. Summation and averaging are simple combination schemes. A more general approach involves assigning a weight to each base learner. In the present study, a heterogeneous ensemble model was built based on QRNNs and KDE. QRNNs serve as base learning algorithms to produce multiple base learners, and the probability combination scheme based on KDE is used to combine the base learners into the final ensemble prediction.

Quantile Regression
Quantile regression is a common statistical technique for conducting inferences concerning conditional quantile functions [34,35]. More formally, any real-valued random variable Y may be characterized by its distribution function as follows: whereas for any 0 < τ < 1, is called the τth quantile of Y. Given a data set (x i (t), Y(t)) for i = 1, 2, · · · , I and t = 1, 2, · · · , N, the linear quantile regression can be expressed as follows:Ŷ where 0 < τ < 1 is the quantile, and b is an error with zero expectation. The estimated parameters θ i can be approximated by minimizing a sum of the asymmetrically weighted absolute residual cost functions, which are expressed as follows: where Y(t) is the observation at time t and ρ τ is the check function, which is also known as the pinball loss function and is defined as follows:

Quantile Regression Neural Network
Given inputs x i (t) and an output Y(t), the output from a QRNN is calculated as follows: Consider a hidden-layer transfer function h(·); the output from the j-th hidden-layer node g j (t) is given by applying the hidden-layer transfer function to the inner product between x i (t) and hidden-layer weights w (h) ij plus the hidden-layer bias b j (h) , which can be calculated as follows: An estimate of the conditional τ-quantileŷ τ (t) iŝ where w j are the output-layer weights, b (o) is the output-layer bias, and f (·) is the output-layer transfer function. The transfer function h(·) and f (·) are usually set as the hyperbolic tangent sigmoidal and linear function, respectively [36].
As an alternative method to prevent overfitting, weight delay regularization for the magnitude of the input-hidden layer weight can be applied by setting a penalty with a nonzero value.

Kernel Density Estimation (KDE)
Nonparametric density estimation is the process of fitting a parametric density model of a random variable without making the assumption that the density belongs to a particular parametric family [37,38]. Various methods have been proposed for nonparametric density estimation, e.g., k-nearest neighbors method, Parzen windows, histogram, and KDE [38]. In the domain of nonparametric density estimation, the K-nearest neighbors method has a very limited scope of practical applications due to its very poor performance. The Parzen windows method presents slightly better performance but also produces discontinuities (stair-like curves) that are quite annoying in practice [38]. A histogram is a simple form of the nonparametric density estimation. However, it suffers serious and noticeable drawbacks. First, the resulting visualization strongly depends on the choice of binning. Second, the natural feature of the histogram is discontinuity, which causes extreme difficulty if derivatives of the estimates are required.
Fortunately, those abovementioned drawbacks can be easily eliminated by using KDE [38,39]. In fact, KDE has been extensively studied and has become the most popular method in nonparametric 6 of 23 density estimation. Given a random sample Y 1 , Y 2 , · · · , Y m , the value of the density at the point y estimated by the KDE method is given by the following: where h is the bandwidth with positive real value and K(·) is the kernel function. In this study, the most effective Epanechnikov kernel [38] was adopted and expressed as where R(·) is the indicator function, that is, R(y ∈ A) = 1 for y ∈ A and R(y ∈ A) = 0 for y A.
The selection of bandwidth parameter is a crucial issue in KDE. The bandwidth parameter influences the smoothness of the KDE curve and also determines the tradeoff between the bias and variance. In general, the smaller the bandwidth, the smaller the bias, and the larger the variance. A number of methods have been proposed to find the optimal bandwidth, such as Silverman's rule of thumb and the Sheather-Jones method. Silverman's rule of thumb bandwidth with a Gaussian kernel and Epanechnikov kernel can be computed as follows: h optimal ≈ 1.06σn − 1 5 (10) h optimal ≈ 2.34σn − 1 5 (11) whereσ is the estimation of σ (standard deviation of the input data) [38].

Ensemble Prediction Employing QRNNs and KDE
The proposed ensemble prediction employing QRNNs and KDE is shown in Figure 2. The QRNNs-KDE approach consists of four stages: (1) data splitting and normalization, (2) QRNN modelling, (3) probability density function (PDF) estimation by KDE, and (4) final ensemble prediction. thumb and the Sheather-Jones method. Silverman's rule of thumb bandwidth with a Gaussian kernel and Epanechnikov kernel can be computed as follows: where σ is the estimation of σ (standard deviation of the input data) [38].

Ensemble Prediction Employing QRNNs and KDE
The proposed ensemble prediction employing QRNNs and KDE is shown in Figure 2. The QRNNs-KDE approach consists of four stages: (1) data splitting and normalization, (2) QRNN modelling, (3) probability density function (PDF) estimation by KDE, and (4) final ensemble prediction.
Data splitting and normalization: The original landslide monitoring dataset is divided into training data and testing data. The training data are used for model construction, and the testing data are used to evaluate the performance of the constructed model. To eliminate the influence of dimensional data, the training data and testing data are first normalized in the range of 0 to 1.
QRNNs modelling: QRNNs serve as base learning algorithms to generate multiple base learners by applying a finite number of conditional quantities 1 2 m τ τ τ ≤ ≤ ≤  within the domain 0 1 τ < < , e.g., τ = 0.01, 0.02, …, 0.98, 0.99. The base learners of landslide displacement are obtained after renormalizing the outputs from the QRNNs approach. To avoiding overfitting in QRNNs modelling, a penalty parameter with nonzero value is applied. PDF estimation by KDE: Multiple base learners from the QRNNs base model are treated as the input for KDE to estimate the probability density function (PDF) of the base learners. The kernel function and bandwidth influence the shape of the KDE curve. An appropriate kernel function and an optimal bandwidth should be chosen to best match the features of the original dataset.
Final ensemble prediction: In the present study, the final ensemble prediction was obtained through a probability combination scheme as follows: where ( )  Data splitting and normalization: The original landslide monitoring dataset is divided into training data and testing data. The training data are used for model construction, and the testing data are used to evaluate the performance of the constructed model. To eliminate the influence of dimensional data, the training data and testing data are first normalized in the range of 0 to 1.
PDF estimation by KDE: Multiple base learners from the QRNNs base model are treated as the input for KDE to estimate the probability density function (PDF) of the base learners. The kernel function and bandwidth influence the shape of the KDE curve. An appropriate kernel function and an optimal bandwidth should be chosen to best match the features of the original dataset.
Final ensemble prediction: In the present study, the final ensemble prediction was obtained through a probability combination scheme as follows: where p i (t) is the probability value of the i-th base learner and Y i (t) is obtained from the KDE for monitoring period t.

Evaluation Metrics and Uncertainty Quantification
In this study, five indices-coefficient of determination (R 2 ) MSE, RMSE, NRMSE, and MAPEwere applied to assess the performance of point prediction. R 2 , MSE, RMSE, NRMSE, and MAPE are defined as whereû t and u t denote the t-th predictive value and observation, respectively, and u andû denote the mean of the observation and the mean of the predictive value, respectively. In the present study, the associated predictive uncertainties were quantified with PIs. After the above procedures, full PDFs of the future landslide displacement were achieved. An interval prediction with a (1 − α) × 100% confidence interval can be obtained from the α/2 and 1 − α/2 quantiles of the obtained PDF. The α level, also called the significance level, ranges from 0 to 1 and is the probability of not capturing the value of the parameter. The predictive values of the α/2 quantity and 1 − α/2 quantity are set as the upper bound (U 1−α t ) and lower bound (L 1−α t ), respectively. For example, a 90% central PI can be obtained from the 0.05 and 0.95 quantiles of the PDF. The upper bound and lower bound of the 90% confidence level correspond to the predictive values of the 0.95 and 0.05 quantiles of the obtained PDF.
The prediction interval coverage probability (PICP), normalized mean PI width (NMPIW), and coverage width-based criterion (CWC) are three indices for evaluating the correctness of the approximated PIs. The PICP reflects the degree of reliability of PIs and is defined as where I 1−α t is defined as follows: NMPIW measures the width of the PI; it is defined as where ς is the range of the underlying targets.
For high-quality PIs, narrow PIs (smaller NMPIW) with a high coverage probability (large PICP close to 100%) have great value [40,41]. Theoretically, NMPIW and PICP are conflicting. Therefore, CWC, which is a new balance criterion between PICP and NMPIW [42], is proposed to give a comprehensive assessment of PIs. CWC is defined as where ψ is a small positive value within the range of (0.1%, 0.5%), µ corresponds to the nominal confidence level associated with PIs that is usually set to 1 − α, and δ is a small positive value less than 1. γ is set to 1 during the training process; for testing, it is defined by the following step function:

Features of the Fanjiaping Landslide
The Fanjiaping landslide is located on the southern bank of the Yangtze River and upstream of the Baishuihe landslide and downstream of the well-known Huangtupo landslide, which is approximately 56 km northwest of the Three Gorges Reservoir Dam (see Figure 3 for location). The Fanjiaping landslide is an ancient landslide [43,44]   The site-specific investigation shows that the landslide materials are arranged in two different layers: a colluvial deposit at the upper surface and highly disturbed sandstone at the lower surface. The cataclastic sandstone is underlaid by sandstone and mudstone of the Jurassic Xiangxi formation (J1x) with an average dip direction of 10-25° and a dip angle of 27-36° (Figure 4b,d). Soft coal layers are prevalent in the J1x formation, and many landslides have developed along the soft coal layers. The borehole data indicates that the landslide mass of the Muyubao and Tanjiahe landslide slide along a soft coal layer with a thickness ranging from 0.1 to 0.3 m. According to laboratory testing of sliding zone soil obtained from the borehole, the natural moisture content of the soil is 12.6%, and the natural density is 1.9 g/cm 3 . The site-specific investigation shows that the landslide materials are arranged in two different layers: a colluvial deposit at the upper surface and highly disturbed sandstone at the lower surface. The cataclastic sandstone is underlaid by sandstone and mudstone of the Jurassic Xiangxi formation (J1x) with an average dip direction of 10-25 • and a dip angle of 27-36 • (Figure 4b,d). Soft coal layers are prevalent in the J1x formation, and many landslides have developed along the soft coal layers. The borehole data indicates that the landslide mass of the Muyubao and Tanjiahe landslide slide along a soft coal layer with a thickness ranging from 0.1 to 0.3 m. According to laboratory testing of sliding zone soil obtained from the borehole, the natural moisture content of the soil is 12.6%, and the natural density is 1.9 g/cm 3 .

Input Data
A total of sixteen GPS beacons were installed on the landslide mass to monitor the landslide movements in September 2006 (see Figure 4 for the GPS locations): four on the Tanjiahe landslide and twelve on the Muyubao landslide. The GPS monuments were manually surveyed once a month. In April 2016, four GPS monitoring points, ZG295, ZG296, ZG297, and ZG298, were updated to near real-time monitoring. At most, thirteen years' worth of monitoring data were obtained. Figure 5 shows the monthly rainfall intensity obtained from the Shazhenxi Meteorological Station near the Fanjiaping landslide, the reservoir water level, and the displacement from GPS survey monuments

Input Data
A total of sixteen GPS beacons were installed on the landslide mass to monitor the landslide movements in September 2006 (see Figure 4 for the GPS locations): four on the Tanjiahe landslide and twelve on the Muyubao landslide. The GPS monuments were manually surveyed once a month. In April 2016, four GPS monitoring points, ZG295, ZG296, ZG297, and ZG298, were updated to near real-time monitoring. At most, thirteen years' worth of monitoring data were obtained. Figure 5 shows the monthly rainfall intensity obtained from the Shazhenxi Meteorological Station near the Fanjiaping landslide, the reservoir water level, and the displacement from GPS survey monuments over the thirteen-year period from October 2006 to March 2018. The available data indicate that the landslide was unstable and continuously deforming during the entire monitoring period. The landslide exhibits a step-like deformation behavior because of the periodic fluctuations in the reservoir water level and heavy precipitation. The monitoring data from both Muyubao and Tanjiahe show that larger displacements occurred in the upper middle part of the landslide mass. From the sequence of the surface cracks and displacement magnitude, we speculate that the movement occurred first at the rear part and progressed downslope. Based on a previous study on the relations between slip-surface geometry, material structures, and deformational structures [45,46], the observed kinematic behaviors are expected independent of the characteristics of the landslide material. However, more work is needed to confirm these findings.  [45,46], the observed kinematic behaviors are expected independent of the characteristics of the landslide material. However, more work is needed to confirm these findings.

Triggering Factors of the Landslide Movements
Although the Fanjiaping landslide is one of the largest landslides in the TGRA, very few publications have reported detailed information on the triggering factors of the landslide movements. Fully understanding the triggering factors is critical for landslide mitigation and early warning. In this study, long-term and near real-time monitoring data were used to comprehensively analyze the landslide movements.  (Figure 6f). The results of the above analysis suggest that landslide deformation occurred at the preliminary operation phase, and more significant movement is likely to occur when the reservoir water level reaches a new higher level.
(2) A large deformation occurred when the reservoir water level slightly dropped from 175 m to 170 m in November to February (I in Figure 6). During 2009 to 2015, the monthly deformation rate during this drawdown period was greater than 20 mm per month (Figure 6a-d). For example, when the reservoir water level dropped from 174 m to 172.21 m in January 2012 to February 2012, the

Triggering Factors of the Landslide Movements
Although the Fanjiaping landslide is one of the largest landslides in the TGRA, very few publications have reported detailed information on the triggering factors of the landslide movements. Fully understanding the triggering factors is critical for landslide mitigation and early warning. In this study, long-term and near real-time monitoring data were used to comprehensively analyze the landslide movements. The cumulative displacement at monitoring point ZG295, monthly rainfall intensity, and reservoir water levels in 2009, 2011, 2012, 2015 are shown in Figure 6a-d. The available data shows the following trends: reservoir dropped from 174.45 to 171.18 m on November 9 2016 to January 16 2017, deformations of 24.07 and 26.46 mm occurred at ZG296 and ZG297, respectively. The corresponding monthly deformation rates were 12.8 and 13.2 mm/month, respectively.
From the above analysis we can conclude that landslide movement was especially pronounced under prolonged periods of dropping reservoir levels, especially during periods of slight dropdown at the highest reservoir level, and the minimum triggering threshold consisted of episodes lasting one month, with cumulative rainfall exceeding 158 mm.  (Figure 6f). The results of the above analysis suggest that landslide deformation occurred at the preliminary operation phase, and more significant movement is likely to occur when the reservoir water level reaches a new higher level.
(2) A large deformation occurred when the reservoir water level slightly dropped from 175 m to 170 m in November to February (I in Figure 6). During 2009 to 2015, the monthly deformation rate during this drawdown period was greater than 20 mm per month (Figure 6a-d). For example, when the reservoir water level dropped from 174 m to 172.21 m in January 2012 to February 2012, the displacements measured at monitoring points ZG295, ZG296, ZG297, and ZG298 were 38.98, 26.39, 35.12, and 41.78 mm, respectively (Figure 6c). However, when the reservoir water level significantly dropped from 170 m to 145 m in February to June (II in Figure 6), the monthly deformation rate decreased to less than 20 mm per month.
(3) When the reservoir level remained at 145 m in July to September (III in Figure 6) and the landslide area suffered a heavy rainfall event, the landslide deformation was likely to be suspended except for 2012. The maximum monthly rainfall intensity during those suspended activities was 158 mm. In July 2012, the landslide area suffered from a heavy rainfall event with a monthly rainfall intensity of 208 mm, and the monitoring point deformed at a high rate. From those comparative analyses, we can speculate than the minimum triggering threshold consists of episodes lasting one month with cumulative rainfall exceeding 158 mm.
(4) When the reservoir rose from 145 m to 175 m in September to November (IV in Figure 6), monthly deformation rate decreased to a small positive (less than 10 mm per month) or even negative value.
(5) The near real-time monitoring data also showed the abovementioned trends: when the reservoir water level rose from 147. 25  From the above analysis we can conclude that landslide movement was especially pronounced under prolonged periods of dropping reservoir levels, especially during periods of slight dropdown at the highest reservoir level, and the minimum triggering threshold consisted of episodes lasting one month, with cumulative rainfall exceeding 158 mm.

Data Splitting and Normalization
The available data ( Figure 5) indicate that for the two active blocks, the largest displacements were observed at monitoring points ZG289 and ZG291, respectively. Therefore, monitoring points ZG289 and ZG291 were selected to establish a prediction model for the Fanjiaping landslide.
Previous correlation analysis in [3] revealed that weak to very strong correlations exist between landslide displacement and triggering and state variables. Therefore, based on triggering factor analysis and previous work on correlation analysis in [3], seven variables including four trigger variables and three state variables were selected as the inputs: rainfall intensity over the past month (x 1 (t)), rainfall intensity over the past two months (x 2 (t)), average reservoir water level in the current month (x 3 (t)), variation in the reservoir water level in the current month (x 4 (t)), displacement over the past one month (x 5 (t)), displacement over the past two months (x 6 (t)), and displacement over the past three months (x 7 (t)). In addition, the displacement in the current month (Y(t)) was selected as the output. A data set (x i (t), Y(t)), i = 1, 2, · · · , 7 was generated based on the inputs and corresponding outputs. For the Tanjiahe landslide, the data from October 2006 to January 2015 with a size of 100 were treated as the training set, and the data from February 2015 to June 2015 with a size of 5 were used as the testing set. For the Muyubao landslide, the data from October 2006 to January 2015 with a size of 133 were treated as the training set, and the data from November 2017 to October 2018 with a size of 12 were used as the testing set.

QRNN Modelling
Two nonlinear models with a sigmoidal transfer function and linear transfer function for τ = 0.01, 0.02, . . . , 0.98, 0.99 with an interval of 0.01 were trained for monitoring points ZG289 and ZG291 to generate multiple base learners. The number of hidden nodes in the QRNNs model was set to 5. The penalty for weight delay regularization was set to 1 to prevent overfitting in QRNNs model construction. For each monitoring period, a total of 99 base learners were obtained at conditional quantities ranging from 0.01 to 0.99 based on the QRNNs. The main parameters applied in the modelling of QRNNs are shown in Table 1.

PDF Estimation by KDE
The multiple base learners from QRNNs were employed as inputs of Epanechnikov KDE to estimate the PDF. The optimal bandwidths for PDF estimation were calculated based on Silverman's rule of thumb. The optimal bandwidths for PDF estimation of testing data at ZG289 were set to 7.98, 8.05, 5.66, 7.05, and 9.26. The optimal bandwidths for PDF estimation of testing data at ZG291 were set to 5.28, 8.61, 7.77, 5.62, and 5.30.

Final Ensemble Prediction
Final ensemble predictions for the Fanjiaping landslide were generated through a probability combination scheme. PIs were constructed from the obtained PDF to estimate the predictive uncertainty. For the purpose of aiding decision-making, it is preferable to have prediction information with high confidence levels to reduce risks. Therefore, PIs at a high PINC value of 90% were obtained and analyzed in the study.

Results
PDFs: The PDFs of predictive displacement at ZG289 and ZG291 constructed by the proposed QRNNs-KDE approach are shown in Figures 7 and 8. The fast movement is the main concern in landslide displacement prediction. Here, only a portion of the prediction describing the fast landslide is selected and shown. Figures 7 and 8 show that rather than a single estimate, the range and complete PDF of the predictive displacement are provided by the proposed approach. All landslide displacement observations are distributed in the middle of the PDFs with high probability in addition to the observations of May and June at ZG289, which appear at the tail of the probability density curve. The small fraction falling into the right tail follows the increase in the prediction period; here, there are more uncertainties associated with longer-term landslide predictions. Final ensemble predictions for the Fanjiaping landslide were generated through a probability combination scheme. PIs were constructed from the obtained PDF to estimate the predictive uncertainty. For the purpose of aiding decision-making, it is preferable to have prediction information with high confidence levels to reduce risks. Therefore, PIs at a high PINC value of 90% were obtained and analyzed in the study.

Results
PDFs: The PDFs of predictive displacement at ZG289 and ZG291 constructed by the proposed QRNNs-KDE approach are shown in Figures 7 and 8. The fast movement is the main concern in landslide displacement prediction. Here, only a portion of the prediction describing the fast landslide is selected and shown. Figures 7 and 8 show that rather than a single estimate, the range and complete PDF of the predictive displacement are provided by the proposed approach. All landslide displacement observations are distributed in the middle of the PDFs with high probability in addition to the observations of May and June at ZG289, which appear at the tail of the probability density curve. The small fraction falling into the right tail follows the increase in the prediction period; here, there are more uncertainties associated with longer-term landslide predictions.   Final ensemble prediction: Figure 9 shows the final ensemble predictions. As shown in Figure 9, the ensemble predictions obtained via the probability combination scheme showed a high degree of consistency in the landslide displacement observations, with coefficient of determination values of 0.999932 and 0.999944. To further evaluate the prediction performances of ensemble prediction based on the QRNNs-KDE, the evaluation metrics of the BP, RBF, ELM, and SVM approaches are shown in Table 2. As shown in Table 2, the final ensemble predictions using the QRNNs-KDE approach outperformed the persistence methods with the smallest MSE, RMSE, NRMSE, and MAPE and the largest R 2 . Moreover, compared with predictions at monitoring point ZG289 using the Copula-KSVMQR approach in [3], the QRNNs-KDE approach provided more accurate prediction with smaller MAPE and RMSE.  Final ensemble prediction: Figure 9 shows the final ensemble predictions. As shown in Figure 9, the ensemble predictions obtained via the probability combination scheme showed a high degree of consistency in the landslide displacement observations, with coefficient of determination values of 0.999932 and 0.999944. To further evaluate the prediction performances of ensemble prediction based on the QRNNs-KDE, the evaluation metrics of the BP, RBF, ELM, and SVM approaches are shown in Table 2. As shown in Table 2, the final ensemble predictions using the QRNNs-KDE approach outperformed the persistence methods with the smallest MSE, RMSE, NRMSE, and MAPE and the largest R 2 . Moreover, compared with predictions at monitoring point ZG289 using the Copula-KSVMQR approach in [3], the QRNNs-KDE approach provided more accurate prediction with smaller MAPE and RMSE. Uncertainty quantification: Based on the PDFs shown in Figures 7 and 8, PIs at a high confidence level (90%) were constructed for ZG289 and ZG291 (Figure 10a,c, respectively). To evaluate the prediction performances based on the QRNNs-KDE approach, 90% PIs were constructed based on the bootstrap-ELM-ANN approach (Figure 10b,d). The corresponding evaluation metrics are shown in Table 3. As shown in Figure 10 and Table 3, the constructed PIs based on the QRNNs-KDE  Uncertainty quantification: Based on the PDFs shown in Figures 7 and 8, PIs at a high confidence level (90%) were constructed for ZG289 and ZG291 (Figure 10a,c, respectively). To evaluate the prediction performances based on the QRNNs-KDE approach, 90% PIs were constructed based on the bootstrap-ELM-ANN approach (Figure 10b,d). The corresponding evaluation metrics are shown in Table 3. As shown in Figure 10 and Table 3, the constructed PIs based on the QRNNs-KDE approach perfectly covered the observations with a high percentage, and the QRNNs-KDE approach outperformed the bootstrap-ELM-ANN approach with smaller NMPIW and CWC. For example, the performance indices NPIW and CWC of 90% PIs at ZG289 were 0.0215 and 0.1661, respectively, which were lower than those obtained using the bootstrap-ELM-ANN approach. The normalized mean PI width using the QRNNs-KDE approach was approximately 90% narrower than that for the bootstrap-ELM-ANN approach.
The experimental results show that the final ensemble predictions based on the QRNNs-KDE approach outperformed the traditional BP, RBF, ELM, SVM, and Copula-KSVMQR algorithms with regard to deterministic point prediction. The QRNNs-KDE approach was more informative than traditional algorithms because it provided the likely range of landslide displacement. The landslide observations were distributed in the middle of the prediction range with high probability. Moreover, regarding the aspect of uncertainty quantification, the QRNNs-KDE provided more satisfactory PIs than the bootstrap-ELM-ANN approach. Therefore, we believe that the final ensemble predictions based on the QRNNs-KDE approach have the advantages of accurate prediction and uncertainty quantification of landslide displacement. Note: The prediction results with a narrower PI range are shown in bold italics.

Discussion
In this study, with regard to point prediction, the probability-scheme combination ensemble prediction, which employs QRNNs-KDE, provided the best prediction. The fundamental reasons behind this can be explained from statistical, computational, and representational perspectives [47]. From a statistical perspective, the available training data set may not be able to provide sufficient information for training the true model (h * in Figure 11). Constructing an ensemble model (h ' in Figure   Figure 10. Comparisons of the observations and the constructed PIs at a 90% confidence level for the Fanjiaping landslide at ZG289 and ZG291 using QRNNs-KDE and bootstrap-ELM-ANN. (a) 90% PIs at ZG291 using QRNNs-KDE; (b) 90% PIs at ZG291 using bootstrap-ELM-ANN; (c) 90% PIs at ZG289 using QRNNs-KDE, (d) 90% PIs at ZG289 using bootstrap-ELM-ANN.

Discussion
In this study, with regard to point prediction, the probability-scheme combination ensemble prediction, which employs QRNNs-KDE, provided the best prediction. The fundamental reasons behind this can be explained from statistical, computational, and representational perspectives [47].
From a statistical perspective, the available training data set may not be able to provide sufficient information for training the true model (h * in Figure 11). Constructing an ensemble model (h ' in Figure 11) might not be better than the single best prediction model h * , but it does reduce the risk of choosing a bad learner with poor generalizability (schematic in Figure 11a). From a computational perspective, in a single model the training algorithms might get stuck in lock optima by only performing a local search. Constructing an ensemble model by searching from different starting positions might be a better alternative (schematic in Figure 11b). From a representational perspective, it is possible that the searched hypothesis space might not contain the true model h * . Constructing an ensemble model might expand the representable space (schematic in Figure 11c). 11) might not be better than the single best prediction model h * , but it does reduce the risk of choosing a bad learner with poor generalizability (schematic in Figure 11a). From a computational perspective, in a single model the training algorithms might get stuck in lock optima by only performing a local search. Constructing an ensemble model by searching from different starting positions might be a better alternative (schematic in Figure 11b). From a representational perspective, it is possible that the searched hypothesis space might not contain the true model h * . Constructing an ensemble model might expand the representable space (schematic in Figure 11c). Figure 11. Schematic that shows the fundamental benefits of the ensemble prediction model from statistical (a), computational (b), and representational (c) perspectives. h * is the true prediction model; h1, h2, and h3 are single prediction models; and h ' is the ensemble prediction model obtained by combining the single prediction models h1, h2, and h3. The outer black curve is the hypothesis space of all possible models. The inner blue curve denotes the subset of hypotheses that give reasonable accuracy with the available training data (modified from [47]).
In the proposed QRNNs-KDE approach, the probability combination scheme is employed to combine 99 base learners into one final ensemble to improve the model performance. However, a concern about computational time may be associated with this ensemble strategy. The required computational time is highly related to the number of base learners. For the case of ZG 291, the required computation time is 191.85 s to train 99 base learners in RStudio Version 1.2.5042 on an Intel(R) Xeon(R) E-2176M @ 2.70 GHz CPU with 64 GB RAM. Thus, we believe that the proposed approach is computationally efficient.
Nevertheless, the probability-scheme combination ensemble prediction, which employ QRNNs and KDE, also holds inherent limitations associated with data-driven models, such as the lack of an explicit input-output relationship, and the requirement of large training data to maintain the model performance.
In practical applications, the main motivation for the construction of predictive range and complete PDF is to quantify the likely predictive uncertainty in the deterministic point predictions. Availability of range and complete PDF of the predictive displacement allows the researchers and practitioners to efficiently quantify the level of predictive uncertainty with the deterministic point predictions and to consider a multiple of solutions/scenarios for the best and worst conditions. Wide ranges are an indication of presence of a high level of uncertainty in the operation. This information can guide the researchers and practitioners to avoid the selection of risky actions under uncertain conditions. In contrast, narrow range means that decisions can be made more confidently with less chance of confronting an unexpected condition in the future, for example, if a sharp displacement increment with a wider range was predicted for the further. An alert should be carefully determined whether reaching tertiary creep stage by researchers and practitioners through comprehensive analysis. Under this circumstance, time-of-failure forecasting should be run in parallel, and a multiple of solutions/scenarios should be considered until either failure precursors are identified or the movements suspended.
The proposed QRNNs-KDE approach is suitable for medium-term to long-term horizon forecasting. Results from previous studies [2,48] have shown that the performance of data-driven Figure 11. Schematic that shows the fundamental benefits of the ensemble prediction model from statistical (a), computational (b), and representational (c) perspectives. h * is the true prediction model; h 1 , h 2 , and h 3 are single prediction models; and h ' is the ensemble prediction model obtained by combining the single prediction models h 1 , h 2 , and h 3 . The outer black curve is the hypothesis space of all possible models. The inner blue curve denotes the subset of hypotheses that give reasonable accuracy with the available training data (modified from [47]).
In the proposed QRNNs-KDE approach, the probability combination scheme is employed to combine 99 base learners into one final ensemble to improve the model performance. However, a concern about computational time may be associated with this ensemble strategy. The required computational time is highly related to the number of base learners. For the case of ZG 291, the required computation time is 191.85 s to train 99 base learners in RStudio Version 1.2.5042 on an Intel(R) Xeon(R) E-2176M @ 2.70 GHz CPU with 64 GB RAM. Thus, we believe that the proposed approach is computationally efficient.
Nevertheless, the probability-scheme combination ensemble prediction, which employ QRNNs and KDE, also holds inherent limitations associated with data-driven models, such as the lack of an explicit input-output relationship, and the requirement of large training data to maintain the model performance.
In practical applications, the main motivation for the construction of predictive range and complete PDF is to quantify the likely predictive uncertainty in the deterministic point predictions. Availability of range and complete PDF of the predictive displacement allows the researchers and practitioners to efficiently quantify the level of predictive uncertainty with the deterministic point predictions and to consider a multiple of solutions/scenarios for the best and worst conditions. Wide ranges are an indication of presence of a high level of uncertainty in the operation. This information can guide the researchers and practitioners to avoid the selection of risky actions under uncertain conditions. In contrast, narrow range means that decisions can be made more confidently with less chance of confronting an unexpected condition in the future, for example, if a sharp displacement increment with a wider range was predicted for the further. An alert should be carefully determined whether reaching tertiary creep stage by researchers and practitioners through comprehensive analysis. Under this circumstance, time-of-failure forecasting should be run in parallel, and a multiple of solutions/scenarios should be considered until either failure precursors are identified or the movements suspended.
The proposed QRNNs-KDE approach is suitable for medium-term to long-term horizon forecasting. Results from previous studies [2,48] have shown that the performance of data-driven models varies for landslides with different deformation behaviors. Usually, for landslides with drastic step-like deformation, the prediction accuracy is lower, and the corresponding prediction error is larger. Therefore, in practical applications of medium-term to long-term horizon forecasting, when predicting landslides with drastic deformation, the proposed QRNNs-KDE approach should be applied with caution. To achieve excellent performance, sufficient data are recommended and needed for model training.

Conclusions
In this study, a QRNNs-KDE approach was proposed to improve the prediction accuracy and uncertainty quantification of landslide displacement. The Fanjiaping landslide in the TGRA was selected as a case study to explore the performance of the QRNNs-KDE approach. The following conclusions from the study were obtained: The movements of the Fanjiaping landslide was especially pronounced under prolonged periods of dropping reservoir levels, especially during periods of slight dropdown at the highest reservoir level, and the minimum triggering threshold consists of episodes lasting one month, with cumulative rainfall exceeding 158 mm.
The QRNNs-KDE approach achieves perfect performance and outperforms the traditional BP, RBF, ELM, SVM, bootstrap-ELM-ANN, and Copula-KSVMQR methods. Additionally, the proposed approach is more informative by providing the likely range and complete PDFs of landslide displacement. The landslide displacement observations are distributed in the middle of the prediction range with high probability.
In practical application, the proposed QRNNs-KDE approach is suitable for medium-term to long-term horizon forecasting. The range and complete PDF of the predictive displacement can supplement final point predictions for decision making.