Time series decomposition is a crucial preprocessing technique that enables the handling of nonlinear and nonstationary signal data in real-world applications, allowing models to effectively extract both high-frequency and low-frequency information from raw data, thereby enhancing the accuracy of predictive models (Talukdar et al., 2020; Youssef et al., 2016; Dong et al., 2023). Landslide displacement time series exhibit nonlinear and nonstationary characteristics, with general nonstationary data typically composed of components such as trend, seasonality, periodicity, or random noise (Min et al., 2019). By decomposing the original series, characteristic information at different frequencies can be extracted, facilitating the model’s analysis of each component in the raw data and uncovering underlying data patterns and developmental trends (Alabi et al., 2022). Applying predictive models to each decomposed component significantly improves the model’s predictive performance. This paper primarily employs the Variational Mode Decomposition (VMD) algorithm for this purpose.

Variational Mode Decomposition (VMD) is a novel method for decomposing complex signals, introduced by K. Dragomiretskiy and D. Zosso in 2014 as an advancement over Empirical Mode Decomposition (EMD) (Manibardo et al., 2022). Unlike EMD, VMD redefines the Intrinsic Mode Function (IMF) as an amplitude- and frequency-modulated (AM-FM) signal, denoted as u k t , with the following mathematical representation in Equation 1: u k t = A k t cos ϕ k t (1)

In this formulation, A k t represents the instantaneous amplitude, and ϕ k t denotes the phase.

In Variational Mode Decomposition (VMD), the signal sifting approach used in Empirical Mode Decomposition (EMD) is abandoned for decomposing Intrinsic Mode Functions (IMFs). Instead, after configuring parameters such as the number of modes K , penalty parameter a , and update step size τ , VMD decomposes the signal by solving an optimization problem within a variational framework. The mathematical formulation for this model is given in Equation 2: min u k , ω k ∑ k = 1 K ∂ t σ t + j π t u k t e − j ω k t 2 2 s . t .   ∑ k = 1 K u k = f t (2)

In the equations: u k represents the decomposed components obtained through the final decomposition process; u k = u 1 , . . . , u K ; ω k corresponds to the actual center frequencies of each Intrinsic Mode Function (IMF) component; ω k = ω 1 , . . . , ω K ; σ t + j / π t u k t denotes the analytic signals derived for each IMF component; e − j ω k t signifies the estimated center frequencies of these analytic signals; f t represents the original input signal.

To solve for the optimal solution of this constrained variational model, Equation 4 must be reformulated into an unconstrained optimization problem: a ∑ k ∂ t σ t + j π t u k t e − j ω k t 2 2 + f t − ∑ k u k t 2 2 + λ t , f t − ∑ k u k t (3)

In the equations: λ represents the Lagrange multiplier.

The Alternating Direction Method of Multipliers (ADMM) is utilized to compute the saddle point of the unconstrained model derived earlier, which corresponds to the optimal solution of the constrained variational model described in Equation 3. This process decomposes the original signal into K Intrinsic Mode Function (IMF) components. By employing this approach, each IMF component inherits distinct scale-dependent characteristics, thereby better preserving its inherent structural information. Consequently, this study adopts Variational Mode Decomposition (VMD) to analyze landslide displacement-time curves, leveraging its capacity to isolate meaningful components that reflect the underlying dynamics of slope deformation.

LSTM (Long Short-Term Memory) networks safeguard and regulate the information within their memory cells through three gates, with the manipulation of this information achieved via element-wise multiplication by activation functions. A set of parameters, trained through gradient descent, is employed to control the state of each gate (Zhang et al., 2024). Each gate in an LSTM serves a distinct and specific function. The forget gate, denoted as f , determines which pieces of information to discard from the previous state h t − 1   . After the input x t and h t − 1   are processed by the update gate u , they, in conjunction with the adjusted forget gate f , decide the weightage to be given to the candidate state h ∼ t for updating the current state h t . To generate the output y t , the network first filters its current state using a nonlinear function g 2 , and then combines it with the output of the output gate ο to produce y t . Here, a portion of the state y t is fed back as the next input y t − 1 . Each gate’s operation hinges on the current external input x t and the previous output (or more accurately, the previous hidden state y t − 1 , but again, adhering to the original text’s phrasing). Figure 1 provides a detailed illustration of the LSTM’s structural principles, and the equation governing the update of the LSTM’s state is presented in Equation 4. f o r g e t g a t e : f t = σ W f x t + R f y t − 1 + b f c a n d i d a t e   s t a t e : h ∼ t = g 1 W f x t + R f y t − 1 + b f u p d a t e   s t a t e : u t = σ W u x t + R u y t − 1 + b u c e l l   s t a t e : h t = u t ⊗ h ∼ t + f t ⊗ h t − 1 o u t p u t   g a t e : ο t = σ W ο x t + R ο y t − 1 + b ο o u t p u t : y t = ο t ⊗ g 2 h t (4)

In the equation, x t represents the input vector at time step t; W f , W h , W u and W ο are weight matrices associated with the input units and connections to the hidden layer; R f , R h , R u and R ο are weight matrices for the connections within the hidden layer; b f , b h , b u and b ο are bias vectors; the activation function σ is either the sigmoid or tanh function; g 1 and g 2 are nonlinear activation functions; and ⊗ denotes element-wise multiplication.

The effectiveness of predictive models is evaluated and compared using the Root Mean Square Error (RMSE) and the coefficient of determination (R²) (Praveen et al., 2025; Zhang et al., 2021; Huang et al., 2017). The primary formulas for these metrics are as follows in Equations 5, 6:

(1) RMSE

Where, x ^ represents the predicted values, x denotes the actual values, x ¯ is the mean of the actual samples, and N indicates the total number of samples. A smaller RMSE value indicates stronger predictive capability of the model, while an R² value closer to 1 signifies better model fitting performance.

Normalization is a critical preprocessing step in data analysis, particularly for time series prediction tasks like landslide displacement forecasting. The primary objective of normalization is to rescale the input data to a standardized range (typically [0, 1] or [-1, 1]), thereby eliminating the influence of varying magnitudes among different features or samples. This process ensures that all variables contribute equally to the model’s training and prevents dominant features (e.g., large-scale displacement values) from overshadowing others (e.g., subtle rate-of-change indicators) (Naemitabar and Asadi, 2021).

In the context of landslide displacement time series, raw displacement data often exhibit significant variations in scale due to factors such as sensor precision, environmental noise, or cumulative displacement effects. Without normalization, models like LSTM or hybrid VMD-SVR-LSTM architectures may struggle to converge efficiently or could produce biased predictions. For instance, if displacement values range from 0 to 1,000 m while rate-of-change values range from 0 to 1 m/day, the model might disproportionately prioritize displacement magnitude over temporal trends. The most commonly used normalization technique is Min-Max Scaling, defined as in Equation 7: y = x i − x i m i n x i m a x − x i m i n (7)

In the formula: x i m a x represents the maximum value of variable i ; x i m i n represents the minimum value of variable i ; x i denotes the original value; y denotes the normalized value.

The Baishuihe landslide is located on the southern bank of the Yangtze River, approximately 56 km downstream from the Three Gorges Dam site, within Baishuihe Village, Shazhenxi Township (Figure 2). Its geographical coordinates are 110°32′09″ longitude and 31°01′34″ latitude (Dou et al., 2019). The landslide mass is situated in a broad river valley section of the Yangtze, characterized by a monoclinic bedding slope that dips in the same direction as the slope face. The terrain slopes downward from south to north, exhibiting a stepped profile toward the river. The trailing edge of the landslide reaches an elevation of 410 m, bounded by the geotechnical interface, while the leading edge extends to the Yangtze River. The eastern and western boundaries are defined by bedrock ridges, with an overall slope gradient of approximately 30°. The landslide measures 600 m in north-south length and 700 m in east-west width, with an average slide mass thickness of about 30 m and a total volume of 1,260 × 10⁴ m³. Classified as a colluvial landslide, the slope exhibits a dip-slope orientation.

Since June 2003, professional monitoring of the Baishuihe landslide has been conducted, with numerous monitoring stations deployed across the site. Among these, the ZG118 station, located in the central portion of the landslide mass, is capable of capturing the entire evolutionary process of the landslide. This station also boasts relatively comprehensive monitoring data, facilitating data-driven modeling efforts. Given these two advantages—central location and robust datasets—the monitoring data from the ZG118 station were utilized to validate the effectiveness and superiority of the proposed methodology in this study. Figure 3 illustrates the cumulative displacement, reservoir water level, and rainfall curves monitored at the ZG118 station from January 2007 to December 2012. Notably, rapid movement episodes are observed at the onset of the rainy season (May to September annually) and toward the end of reservoir water level drawdowns (June and July). Moreover, these rapid movement periods typically conclude before the rainy season’s end. It is evident that fluctuations in both reservoir water levels and rainfall significantly influence cumulative landslide displacement, suggesting that precipitation and reservoir water level variations are primary factors driving deformation and failure at the Baishuihe landslide.

Based on prior research, nine influential factors affecting cyclical displacement were extracted: maximum monthly rainfall, cumulative monthly rainfall, cumulative rainfall over a two-month period, monthly reservoir water level, monthly change in reservoir water level, two-month change in reservoir water level, monthly cumulative displacement increment, two-month cumulative displacement increment, and three-month cumulative displacement increment. These nine factors have demonstrated high effectiveness in predicting cyclical displacement. For the experiments, data from January 2007 to December 2011 were used as the training set, while data from January 2012 to December 2012 served as the testing set. Predictive models were developed using both Support Vector Regression (SVR) and Long Short-Term Memory (LSTM) networks.

To address the challenge of insufficient data volume for direct modeling of the raw landslide displacement series, this study first employed Variational Mode Decomposition (VMD) to decompose the displacement series into several components, each with a unique central frequency. This decomposition indirectly resolved the data scarcity issue by enabling modeling of individual components rather than the raw series. The cumulative landslide displacement is decomposed into trend and periodic displacements, as shown in Figure 4. From Figure 4, it can be observed that the trend displacement obtained using VMD effectively reflects the changing trend of the cumulative landslide displacement. The periodic displacement is primarily influenced by cyclical precipitation and reservoir water level fluctuations. Therefore, before establishing a predictive model for periodic displacement, it is necessary to extract the factors affecting the periodic displacement of the landslide.

Subsequently, the study leveraged the LSTM network’s strength in extracting temporal sequence features and the SVR’s capability in nonlinear regression to create a hybrid modeling approach, compensating for the limitations of single-model regression performance. This combined strategy effectively tackled the data-driven landslide displacement prediction problem.

The experiments were conducted on the deep learning framework PyTorch-CPU, with the model optimized using the Adam optimizer. Based on the size of the dataset, the number of network training epochs was set to 300, and the batch size was set to 8. The LSTM network had an input layer size of 4, a time window size of 5, a step size of 1, and a hidden layer size of 50. The optimal hyperparameters for the SVR model were determined to be C = 1.015 and g = 0.01. The VMD-LSTM-SVR modeling workflow is illustrated in Figure 5, where the blue box represents the hybrid modeling process and the red box denotes the individual component modeling phase. As shown in Figure 5, the VMD-LSTM-SVR modeling approach can be broken down into the following steps.

In the first step, the landslide displacement cycle data are proportionally divided into training, validation, and testing sets in sequence. The training set is used to train the landslide displacement prediction model, the validation set is employed to tune hyper parameters and validate model performance (specifically during LSTM model training), and the testing set is reserved for forecasting the landslide’s future displacement behavior. The modeling data collectively refers to the datasets utilized during both training and validation phases. The specific proportional splits are illustrated in Figure 6.

In the second step, each sequential component of the training and validation sets is normalized and structured into time-series vectors to enhance convergence speed and prediction accuracy during model training.

In the third step, the degradation feature sequences, standardized using Variational Mode Decomposition (VMD), are decomposed into multiple components, each with a limited bandwidth and a unique central frequency.

In the fourth step, a hybrid LSTM-SVR predictive model is constructed for each decomposed component. Since the LSTM and SVR models operate in parallel, their respective parameters are optimized independently.

In the fifth step, the weights of the hybrid model are optimized based on the principle of minimizing the sum of squared prediction errors (SSPE), expressed mathematically as in Equation 8: min ⁡ ∑ t = 1 N y t − β 1 y ^ L S T M , t − β 2 y ^ S V R , t 2 s . t . β 1 + β 2 = 1   β 1 , β 2 ≥ 0 (8) where β 1 and β 2 represent the predictive weights assigned to the LSTM and SVR models, respectively. The predicted output y ^ t of the hybrid model at time t is calculated as in Equation 9: y ^ t = β 1 y ^ L S T M , t + β 2 y ^ S V R , t (9)

In the sixth step, the LSTM-SVR predictive results for all decomposed components are aggregated, and inverse normalization is applied to obtain the final forecasted landslide displacement time series. The predictive performance is then evaluated by comparing the forecasted series with the ground truth (actual displacement data) on the testing set.

Variational Mode Decomposition (VMD) requires prior specification of the number of decomposed components (Zhu et al., 2025). An excessively high number of components risks introducing modal aliasing and noise interference, while an insufficient number may fail to fully decompose the landslide displacement sequence, thereby hindering the extraction of meaningful information. Thus, an optimal number of VMD components must be determined. This study leverages the property that VMD components possess unique central frequencies, which progressively converge as the number of components increases, to guide the selection of the component count. Table 1 presents the central frequencies of VMD components for K = 2 to 6, where γ quantifies the distinctness (or separation) of the components under different values of K.

γ can be expressed as in Equation 10: γ = 1 k − 1 ∑ k = 1 k − 1 ω k + 1 − ω k (10)

A higher value of γ indicates greater distinctness among the decomposed components. As shown in Table 1, γ peaks when K = 3 and gradually stabilizes around 0.0339 as K increases to 6. For example, when K = 6, the differences between the components are minimal and significantly smaller than the disparities observed at K = 3. Consequently, K = 3 is determined as the optimal number of components for decomposition. The zero-crossing rate (ZCR), defined as the ratio of zero-crossings to the sequence length, is used to distinguish high-frequency and low-frequency components at K = 3. Empirically, a ZCR threshold of 5% is applied to classify components as either high-frequency or low-frequency. The ZCR values for Components 1–3 are 0%, 5.98%, and 16.43%, respectively, clearly identifying Component 1 as low-frequency and Components 2–3 as high-frequency. Low-frequency components capture the overall trend of the data, while high-frequency components reflect short-term fluctuations. Thus, by applying VMD to the original frequency-aliased data, the prediction task is simplified into forecasting the trend and fluctuation sequences separately, enabling the model to better capture underlying data patterns. The IMF1 component can be regarded as reflecting a long-term deformation trend influenced by internal control factors such as the landform, geological structure, etc., of the landslide itself—that is, the trend displacement. The IMF2 component can be considered to represent periodic fluctuations in landslide displacement caused by cyclical external influencing factors like rainfall and reservoir water level fluctuations—that is, the periodic displacement. The IMF3 component is typically induced by sudden, short-term influencing factors, including random factors such as human activities and microseisms. This type of displacement manifests as relatively irregular fluctuations in the time series—that is, the random displacement.

Model ablation experiments in data-driven modeling involve systematically removing certain components or factors to assess their impact on the overall model performance. To evaluate the effectiveness of the proposed model, comparative performance analyses were conducted across six models: SVR, LSTM, LSTM-SVR, VMD-SVR, VMD-LSTM, and VMD-LSTM-SVR. The experimental results are summarized in Table 2.

When comparing the three models (SVR, LSTM, and LSTM-SVR) without VMD integration, LSTM demonstrates superior feature extraction capabilities for time-series data compared to SVR, which does not inherently account for temporal dependencies. By combining LSTM and SVR in a weighted ensemble (LSTM-SVR), the strengths of both models are synergized, leading to improved predictive performance. After incorporating VMD, all models exhibit substantial performance enhancements, with RMSE values dropping from 0.05 to the range of 0.03–0.04. Notably, the R² score for the VMD-LSTM-SVR model surpasses 0.6, indicating strong alignment between predicted and actual trends. At this stage, the VMD-SVR and VMD-LSTM-SVR models outperform the VMD-LSTM model significantly, suggesting that VMD enables SVR to leverage its nonlinear regression strengths effectively for time-series forecasting. Figure 7 illustrates the ablation experiment results, highlighting stark differences in predictive performance between models with and without VMD. This disparity arises because VMD decomposes the original sequence into trend and fluctuation components, which are modeled separately and then aggregated. Consequently, VMD-enhanced models capture both the overarching trends (via trend component learning) and fine-grained fluctuations (via fluctuation component learning), yielding substantial performance improvements across all tested architectures.

The improvement effects are illustrated in Figure 8. Compared to the single LSTM and SVR models, the LSTM-SVR model demonstrates a reduced residual range and a residual mean closer to 0, indicating improved prediction performance. After applying VMD to the original sequence, the residual range is significantly narrowed down to −0.04 to 0.1. The quartiles in the boxplot and the residual mean also notably approach 0. Therefore, the VMD-LSTM-SVR model not only captures the underlying trend effectively but also exhibits smaller prediction errors and a narrower residual range, resulting in more robust prediction outcomes.

Both EMD (Empirical Mode Decomposition) and EEMD (Ensemble Empirical Mode Decomposition) are algorithms for sequence decomposition, sharing a modeling and prediction workflow similar to VMD (Variational Mode Decomposition). However, they fundamentally differ from VMD by relying on envelope and extremum-based principles, which do not inherently guarantee modal or spectral separation between decomposed components. This limitation often leads to modal/frequency aliasing, hindering improvements in model predictive performance. To effectively demonstrate the advantages of the VMD algorithm, this study establishes comparative models EMD-LSTM-SVR and EEMD-LSTM-SVR for performance evaluation. The experimental results are summarized in Table 3.

A comparative analysis of Tables 2, 3 reveals that, compared to the R² value of 0.3061 achieved by the LSTM-SVR model in the ablation experiments, the LSTM-SVR models integrated with sequence decomposition algorithms (e.g., EMD, EEMD, VMD) demonstrate varying degrees of improvement in R². This enhancement stems from the inherent properties of sequence decomposition algorithms, which isolate components representing trend sequences to varying extents, thereby boosting predictive performance. However, the remaining components capturing fluctuation sequences contribute differentially to the model’s predictive capacity. For instance, the EMD-LSTM-SVR model shows no improvement or even a slight decline in RMSE compared to the standalone LSTM-SVR, suggesting that the fluctuation components extracted by EMD may not provide as much predictive value in this context.