An Improved VMD–EEMD–LSTM Time Series Hybrid Prediction Model for Sea Surface Height

1. Introduction

In recent years, the continual rise in sea-level has had severe social impacts on coastal areas, including the degradation of freshwater resources, damage to infrastructure, and the depletion of agricultural resources [1,2,3,4,5]. The Sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) highlighted that under the influence of human activities, the rate of the rise in sea level has been steadily accelerating. Between 1901 and 1971, the average rate of the rise in sea level was 1.3 millimeters per year, which increased to 1.9 millimeters per year between 1971 and 2006, and rose further to 3.7 millimeters per year between 2006 and 2018 [6,7,8]. To address the threats posed by rising sea-level, accurate predictions of future changes in sea level are of paramount importance for the sustainable development of coastal regions [7,8,9].

Currently, time series forecasting methods for predicting sea-level primarily fall into two categories: statistical methods and machine learning methods [10,11,12,13,14,15,16]. Statistical methods are widely used in time series forecasting and are rooted in the core concept of conducting statistical analyses on historical data to capture the patterns and trends for predictive purposes [17,18,19,20]. Representative models in this category include ARIMA (auto-regressive integrated moving average) and exponential smoothing [21,22,23,24,25,26]. However, these models face challenges when dealing with complex nonlinear data, as they require manual selection of the features and adjustments of the model’s parameter, leading to some degree of systematic bias and limitations [27,28,29]. Traditional machine learning methods, while capable of automatically extracting features and adjusting the model’s parameters, may not adequately extract the features when dealing with large-scale and highly complex data, resulting in limited predictive performance [30,31]. In contrast, deep learning, as an emerging machine learning technology, has demonstrated a strong capacity for generalization and adaptability to complex data through specific neural network architectures and training methods. It is becoming increasingly prevalent in time series data of sea-level characterized by nonlinearity, non-stationarity, and multivariable attributes. Makarynskyy et al. (2004) utilized artificial neural networks (ANN) to perform multistep predictions based on measured sea level data from a tidal station in Australia, affirming the feasibility of using neural network methods in prediction sea-level [32]. Balogun and Adebisi (2021) conducted comprehensive predictions and comparisons of changes in sea level along the west coast of Peninsular Malaysia using three models: ARIMA, support vector regression (SVR), and long short-term memory (LSTM). Their research findings validated the superiority of the LSTM model in predicting sea-level [33].

With the continued advancement of deep learning methods in the field of time series forecasting, various domains of time series prediction have witnessed the superior predictive accuracy achieved by combining data decomposition techniques with deep learning, as supported by research findings [34,35,36]. Song et al. (2022) conducted multifaceted comparisons of various data decomposition methods, such as complementary ensemble empirical mode decomposition (CEEMD), time-varying filtering-based empirical mode decomposition (TFV–EMD), wavelet transform (WT), and the fusion of these methods with the Elman neural network (ENN) in minute-scale time series predictions of sea level. Their study confirmed that the TVF–EMD–ENN model exhibits the best predictive performance [37]. Wang et al. (2021) incorporated time series of wind speed that had been secondarily decomposed using CEEMD and wavelet packet decomposition (WPD) into a gated recurrent unit (GRU) for making predictions. Their experiments demonstrated that this hybrid model significantly improved the accuracy of short-term wind speed predictions [38]. Variational mode decomposition (VMD), as an emerging data decomposition technique, has stood out among the various decomposition methods due to its unique non-recursive variational approach and exceptional decomposition capability. It has been widely applied in the field of mixed deep learning time series prediction [39,40]. Wang et al. (2020) combined the VMD–LSTM model and used an improved particle swarm optimization algorithm (IPSO) to optimize model parameters, confirming the higher predictive accuracy of the VMD–LSTM model in photovoltaic short-term power time series prediction [41]. Huang et al. (2022) conducted a comparison between Empirical Mode Decomposition (EMD) and VMD, revealing that VMD exhibited stronger noise removal capabilities. They also verified the higher precision of the VMD–LSTM model in predicting variations in coal thickness [42]. Han et al. (2019) performed multifaceted comparisons of various prediction models, including VMD–LSTM, Persistence (PER), Wavelet (WT), and BP neural networks. Their research validated that the VMD–LSTM model exhibited higher accuracy in wind power prediction [43].

The aforementioned studies across various time series domains have consistently demonstrated the superior predictive accuracy of the VMD–LSTM model. However, in practical applications, due to variations in VMD parameter settings and data characteristics, incomplete VMD decomposition may occur. This results in residual components that still contain a certain level of fluctuations and nonwhite noise elements. This inadequately processed information can potentially have a detrimental impact on the predictive accuracy of the VMD–LSTM model, particularly in forecasting complex nonlinear and irregular time series. In light of these challenges, this study proposed a deep learning hybrid model based on the VMD–LSTM model, known as the VMD–EEMD–LSTM model, which incorporates both VMD and ensemble empirical mode decomposition (EEMD). This model aims to enhance the predictive accuracy of the VMD–LSTM model for sea level time series by further processing the residual sequences obtained from VMD decomposition. Evaluation of the VMD–EEMD–LSTM model's predictive performance was conducted with multiple models, various perspectives, and multiple monitoring stations for a comprehensive validation.

2. Principles and Methods

2.1. Signal Processing Methods

VMD, EMD, EEMD, and CEEMDAN (Conformal Empirical Mode Decomposition with Adaptive Noise) are all widely used adaptive methods of data decomposition in the fields of signal processing and data analysis [44,45,46]. Among them, VMD is a completely non-recursive modal decomposition method. Its core idea involves modeling a signal as a variational problem and subsequently seeking the optimal solution through iterative transformations. Ultimately, this process decomposes nonstationary signals into a series of standard orthogonal modal functions. The corresponding principles are as follows [47,48].

(1) With the objective of minimizing the summation of the estimated bandwidths for each modal component ${\mu }_K\left(t\right)$, a constrained variational problem model aimed at identifying the optimal solution. The specific formulation of the constrained variational problem model is provided below.

$$\{\begin{array}{c}\underset{\left\{{\mu }_K\right\},\left\{{\omega }_K\right\}}{{\displaystyle \min }}\left\{{\displaystyle \sum _K{\Vert d_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_K\left(t\right)\right]e^{-j{\omega }_{K}t}\Vert }_2^2}\right\}\\ s.t.{\displaystyle \sum _K{\mu }_K=f}\end{array}$$

(1)

In Equation 1, $j^2=-1$, $\delta (t)$ represents the Dirac function, $\left\{{\mu }_K\right\}$ corresponds to the modal functions obtained after decomposition, $\left\{{\omega }_K\right\}$ denotes the central frequencies associated with each mode, and $f$ represents the original signal.

To attain the optimal solution for the constrained variational problem, the introduction of quadratic penalty factors $\alpha$ and Lagrange multiplier operators ${\lambda }_t$ transforms the problem into an unconstrained variational problem.

$$L\left(\left\{{\mu }_K\right\},\left\{{\omega }_K\right\},\lambda \right)=\alpha {\displaystyle \sum _K{\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast {\mu }_K\left(t\right)\right]e^{-j{\omega }_{K}t}\Vert }_2^2+}{\Vert f\left(t\right)-{\displaystyle \sum _K{\mu }_K\left(t\right)}\Vert }_2^2+\langle \lambda \left(t\right),f\left(t\right)-{\displaystyle \sum _K{\mu }_K\left(t\right)\rangle }$$

(2)

where $L(*)$ represents the augmented Lagrangian function, $\Vert f\left(t\right)-{\displaystyle \sum _K{\mu }_K\left(t\right)}\Vert$ is the quadratic penalty term. Subsequently, an alternatively direction method with multiplier operators is used to solve the unconstrained variational problem, and the optimal solution is obtained by alternating updating ${{\mu }_K}^{n+1}$, ${{\omega }_K}^{n+1}$, and ${\lambda }^{n+1}$.

As a recursive method, the EMD method decomposes the data into a finite number of intrinsic mode functions (IMF) that reflect the inherent properties of the time series signals, along with a "residual sequence" [49,50]. However, the EMD may suffer from the mode-mixing problem in the IMF sequences. To overcome this challenge, this study introduced the EEMD method. EEMD gradually introduces normally distributed white noise into the original signal and then offsets this noise through multiple averaging calculations. This process leads to more precise decomposition of the signal and effectively avoids the mode mixing phenomenon that can occur during the EMD decomposition process [51,52,53]. The specific process is as follows:

(1) Initially, white noise denoted as $\omega (t)$ is introduced into the original signal$x(t)$.

$$x_i(t)=x(t)+{\omega }_i(t),i=1,2,\dots ,m$$

(3)

(2) Subsequently, the EMD method is employed to decompose the initial noisy signal, resulting in n IMF, represented as $C_i(t)$, and a residual sequence represented as $r_i(t)$.

$$x_i(t)={\displaystyle {\sum }_{j=1}^n{C}_{ij}(t)}+r_i(t)$$

(4)

(3) Steps (1) and (2) are iteratively executed for a total of $m$ times, in which white noise is added and IMF components are obtained through decomposition in each iteration. Finally, all the IMF components thus obtained are integrated and averaged to obtain the ultimate result of EEMD signal decomposition.

CEEMDAN introduces an adaptive noise complete set to automatically construct noise components, enabling more effective extraction of modal components in the signal compared to EEMD. This enhances the accuracy and robustness of data decomposition [54].

2.2. Long Short-Term Memory

LSTM is an improved type of recurrent neural network (RNN). Its distinctive memory module is beneficial for handling long-term dependencies and mitigating the challenges related to vanishing and explosion gradients [55]. Compared with traditional neural networks, LSTM networks exhibit pronounced advantages when addressing tasks pertaining to the prediction of lengthy time series data. Consequently, LSTM networks find extensive applications in domains such as time series forecasting [56,57].

The architectural framework of LSTM network comprises an input layer, intermediate hidden layers, and an output layer. Each hidden layer controls storage and access of the data through input gates, forget gates, and output gates, as visually depicted in Figure 1.

As illustrated in the figure, LSTM processes the input of high-temporal data related to sea surface elevation and the previous moment's hidden state output using three gates. The primary process is as follows:

(1) LSTM, through the forget gate (denoted as $f_t$), determines whether to discard or retain information related to $X_t$ and $h_{t-1}$ is governed by the activation function $\sigma$ of the forget gate.

$$f_t=\sigma (W_f\cdot [h_{t-1},X_t]+b_f)$$

(5)

In the equations, $W$ and $b$ represent the weight matrices and biases, respectively. $f_t$ is a vector with values in the range of 0 to 1, where the values within the vector indicate whether information in the cell state $C_{t-1}$ is preserved. A value of 0 implies no preservation, while 1 implies full preservation.

(2) The cell state is updated through the input gate by passing $X_t$ and $h_{t-1}$ to the activation function $\sigma$ to determine the information update. $X_t$ and $h_{t-1}$ are passed to the tanh function to create a new candidate value vector $C{\text{'}}_t$ (where $C{\text{'}}_t$ is a vector in the range of -1 to 1), and the tanh output is multiplied by $\sigma$ output.

$$i_t=\sigma (W_i\cdot [h_{t-1},X_t]+b_i)$$

(6)

$$C{\text{'}}_t=\tanh (W_c\cdot [{\mathrm{h}}_{t-1},X_t]+b_c)$$

(7)

(3) The cell state from the previous layer is element-wise multiplied with the forget vector, and then this value is element-wise added to the output of the input gate, resulting in the updated cell state.

$$C_t=f_t\ast C_{t-1}+i_t\ast C{\text{'}}_t$$

(8)

In the equations, $f_t\ast C_{t-1}$ determines the forgetting of information in $C_{t-1}$, while $i_t*C{\text{'}}_t$ determines the addition of information in $C{\text{'}}_t$ to the new memory cell state $C_t$.

(4) Through the output gate $O_t$, the value of the next hidden state $h_t$ is determined, and this hidden state contains information from previous inputs.

$$O_t=\sigma (W_O\cdot [h_{t-1},X_t]+b_O)$$

(9)

$$h_t=O_t\ast \tanh (C_t)$$

(10)

2.3. The VMD–EEMD–LSTM Hybrid Second-Order Decomposition Prediction Model

VMD and EEMD, as two classical data processing methods, have been widely applied in hybrid modeling. Their effectiveness in enhancing the predictive accuracy of deep learning models has been well-established [58,59,60,61]. The VMD–LSTM model, as a prevalent hybrid deep learning approach, has been widely employed in the realm of time series forecasting. Its applications encompass load forecasting and wind speed prediction, where it has showcased remarkable performance [62,63]. The VMD–LSTM model leverages VMD to perform decomposition of the initial data into a sequence of IMF sequence and a residual sequence. Subsequently, the model individually forecasts each IMF sequence and the residual sequence using the LSTM model. Ultimately, the predicted outcomes of each sequence are aggregated to derive the final model prediction. During the prediction process, as the standard normal mode functions obtained through VMD decomposition are stationary signals, predicting each IMF separately can achieve higher prediction accuracy. However, in practical VMD decomposition, the residual sequence still contains some fluctuating characteristics and high-frequency noise, and their values are relatively large. If these parts of the data are not appropriately processed, they will adversely affect the overall predictive accuracy of the model [64,65,66]. In contrast, the EEMD–LSTM model is a recursive decomposition method, and its main predictive errors are concentrated in the IMF components, which perform well in predicting the residual sequence and overall data. Based on this, this study proposed a deep learning hybrid model called VMD–EEMD–LSTM. This model employs VMD for the initial data decomposition and then utilizes EEMD to further break down the residual components with lower prediction accuracy resulting from the VMD decomposition. Subsequently, each IMF obtained through both VMD and EEMD decomposition is used as a feature used as input into the LSTM model for making predictions. Ultimately, the forecasted outcomes of each IMF are aggregated to yield the model's comprehensive prediction. This approach augments the overall predictive precision of the model by handling the residual components produced by VMD. The detailed procedure is elucidated in Figure 2.

The specific prediction process of the mixed VMD–EEMD–LSTM second-order decomposition model is as follows.

Step 1: Preprocess the time series data on sea level from each station and then input them into the VMD model (with K as the number of components in the model) for decomposition.

Step 2: Take the residual sequence "Residual ₁" obtained from the VMD decomposition, and input it into the EEMD model for further decomposition. This will yield various model components as well as "Residual ₂".

Step 3: Through extensive experiments, it has been determined that among the IMF components obtained through EEMD decomposition, the IMFs after IMF_K (IMF _K+1 to IMF _n) and "Residual ₂" have smaller prediction errors. To mitigate experimental intricacies and guarantee the precision of the model's predictions, the IMF components beyond IMF _K and 'Residual ₂' are combined and utilized as input features for the LSTM model to facilitate the prediction process.

Step 4: Utilizing the distinct IMF components acquired from both VMD decomposition and the EEMD decomposition as distinct features, these components are fed into the LSTM model for prediction purposes. This process yields a total of 2K+1 predictions.

Step 5: Aggregate and amalgamate the 2K+1 predictions to derive the ultimate prediction generated by the VMD–EEMD–LSTM model.

2.4. Evaluation index

To evaluate the precision and dependability of the diverse deep learning models in predicting performance, this study employs the subsequent assessment metrics: RMSE (root mean square error), MAE (mean absolute error), and R² (Coefficient of determination). The definitions of these three-evaluation metrics are elaborated as per references [67,68,69]:

(1)

RMSE

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}}$$

(11)

(2)

MAE

$$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|\left(y_i-{\widehat{y}}_i\right)\right|}$$

(12)

(3)

R²

$$R^2=1-\frac{{\displaystyle \sum _{i=1}^{\mathrm{n}}(y_i-{\widehat{y}}_i}{)}^2}{{\displaystyle \sum _{i=1}^{\mathrm{n}}(y_i-\overline{y}}{)}^2}$$

(13)

where $y_i$ represents the actual values of sea level, ${\widehat{y}}_i$ represents the values predicted by each model, $\overline{y}$ is the mean of the actual values of sea level, and n denotes the total number of data points related to sea level. For RMSE and MAE, smaller values indicate higher predictive accuracy, while for R², values closer to 1 indicate accurate predictions, and values closer to 0 suggest that the model has weaker explanatory power.

To visually assess the enhanced performance of the VMD–EEMD–LSTM model in comparison to other hybrid models across diverse accuracy evaluation metrics, this study introduces the concept of an improvement ratio (I). Through the computation of I, the degree of enhancement achieved by the VMD–EEMD–LSTM model in terms of accuracy can be precisely quantified. The formula for calculating I is defined as follows:

$$I_{y\widehat{y}}=\frac{y-\widehat{y}}{y}$$

(14)

where y and $\widehat{y}$ signify diverse evaluation metrics, y represents the evaluation metric of the hybrid models compared against the VMD–EEMD–LSTM model, while $\widehat{y}$ represents the evaluation metric of the VMD–EEMD–LSTM model. If $I_{y\widehat{y}}$ is greater than 0, it indicates a decreasing trend in the accuracy. If $I_{y\widehat{y}}$ is less than 0, it indicates an increasing trend. The greater the absolute value of $\left|I_{y\widehat{y}}\right|$, the greater the improvement in that evaluation metric for the hybrid model and vice versa.

3. Data and experiments

3.1. Data Preprocessing

The satellite altimetry grid data used in this study were obtained from the European Union's Copernicus Earth Observation Program, specifically from the GLORYS12V1 product (GLOBAL_MULTIYEAR_PHY_001_030). The data have a spatial resolution of 0.083° × 0.083° and a temporal resolution of 1 day [70]. The GLORYS12V1 product is a reanalysis of the global ocean with a 1/12° horizontal resolution and 50 vertical levels, covering sea level measurements from 1993 onwards. It has undergone the necessary standard corrections. To ensure the fairness of the experiments, this study utilized data from six satellite altimetry grid points near the coast of the Netherlands. All the selected data have a consistent temporal coverage, spanning from early 1993 to the end of 2020, totaling 28 years. The distribution of the data is presented in Table 1.

3.2. Experimental pretreatment

3.2.1. Parameter Settings of VMD

Unlike EMD and EEMD, VMD allows for the autonomous selection of the number of mode components obtained during decomposition.

Therefore, in the context of utilizing VMD for data decomposition, the choice of an appropriate number of mode components, denoted as K, is of paramount importance to attain high-quality decomposition outcomes. Opting for a K value that is excessively large may lead to over-decomposition, whereas selecting one that is overly small may result in under-decomposition. To ascertain the optimal K value for the sea level height time series post-decomposition, this investigation employs the signal-to-noise ratio (SNR) as an evaluative criterion for decomposition quality. A higher SNR corresponds to a more distinct signal decomposition and improved noise removal. Following comprehensive experimental inquiry and empirical observations, this study confines the selection of K values to the range of 2 to 10, and identifies the K value within this range that yields the highest SNR as the optimal K value for each individual time series [71,72].

$$\mathrm{SNR}=10\lg \frac{{\displaystyle \sum _{i=1}^N{m}^2(i)}}{{\displaystyle \sum _{i=1}^N{[m(i)-n(i)]}^2}}$$

(15)

where $m(i)$ represents the original signal, and $n(i)$ represents the reconstructed signal. In VMD, the penalty factor α also exerts a certain influence on the decomposition outcomes. Given that the optimal range for the penalty factor α is typically between 1.5 and 2 times the size of the decomposed data [73], and to ensure experimental consistency while considering the size of the decomposed data in the experiments, this study set the penalty factor to 15,000 for all decomposition processes.

Because the range of sites covered in this study was relatively small, the frequency of fluctuation and the amplitude of the sequences of sea level height were quite similar. Therefore, the optimal parameters obtained in the experiments were consistent, all indicating that K=5 was the best number of components for decomposition (Figure 4 in Section 4.2 shows the results VMD decomposition for K=5). To reduce the complexity of the subsequent experiments and ensure experimental consistency, this study combined the data with a K greater than 5 from the IMF obtained by EMD and EEMD decomposition with the residual term for a better predictive analysis.

To further validate the reliability of the selected K value, the LSTM model was employed to conduct comparative experiments for sea level data prediction at the Maassluis station. The experimental results are presented in Table 2.

As presented in Table 2, distinct values of K in VMD decomposition produce residual sequences that manifest substantial predictive errors, constituting the primary source of discrepancies within the VMD-LSTM model. A comparative analysis of predictive outcomes across varying K values reveals that with the escalation of K, the R² for residual sequence predictions gradually diminishes, while the cumulative errors for each IMF increase. This observation implies that the selection of an excessively diminutive K value may result in an inadequate decomposition of the signal, ultimately yielding inferior predictive performance. Conversely, opting for an excessively large K value may lead to an exorbitant decomposition of the signal, which is also not conducive to model prediction.

When K is set at 5, the VMD-LSTM model attains the highest level of predictive accuracy. This reaffirms that, in the context of time series prediction for sea level data, K=5 represents the optimal number of decompositions for VMD.

3.2.2. Parameter Settings of the Model

In deep learning prediction models, a variety of different parameters are involved, and the sizes of the parameters have different degrees of influence on the model’s predictive accuracy. In order to ensure the reliability, this study conducted an experiment by setting the same model parameters. The configuration details of each model are presented in Table 3. In this experimental setup, the parameters for the LSTM model and the hybrid models were set to identical sizes.

According to Table 3, all models employed a uniform data partitioning scheme in this research: the data spanning from early 1993 to the end of 2012 constituted the training set, data from the end of 2012 to the end of 2015 served as the validation set, and data from the end of 2015 to the end of 2020 comprised the test set. This data splitting approach was adopted with the intention of ensuring that the models had access to a sufficient volume of training data to thoroughly grasp the underlying data characteristics. Additionally, employing a substantial dataset for testing purposes enabled a more comprehensive assessment of the model's predictive performance.

4. Results and analysis

4.1. Analysis of the predictions of a single deep learning model

In this section, the predictive performance of four different models, namely ANN [74,75], RNN [76,77], gated recurrent units (GRU) [78,79], and LSTM, was comprehensively evaluated and compared using three different sequences of sea level height. The goal was to determine the best-performing model in terms of time series forecasting, providing a reliable foundation for constructing the subsequent hybrid models. The precise evaluation metrics for the predictions of each model are shown in Figure 3.

As shown in Figure 3, for the three different time series datasets of sea level height, the ANN model exhibited the poorest predictive performance, with an average RMSE of 150.85 mm, an average MAE of 114.06 mm, and an average R² of 0.28 across the different monitoring stations. In contrast, the LSTM model performed the best, with an average RMSE of 137.92 mm, an average MAE of 100.13 mm, and an average R² of 0.40 across the different monitoring stations. LSTM outperformed ANN, RNN, and GRU, demonstrating its superiority. However, since LSTM is a single model, it failed to fully extract the features of the data during training, resulting in a relatively high RMSE and MAE and a relatively low R² for the predictions. This phenomenon highlights the challenge that single models face in accurately capturing all the fluctuations and trends in time series data, especially in complex time series forecasting tasks. Therefore, in the subsequent work of constructing the hybrid models, it is necessary to combine the characteristics of the data decomposition methods to further improve the predictive accuracy of the models.

4.2. Analysis of the hybrid deep learning first-order decomposition model

In response to the issue of insufficient extraction of the features of the data by single models in complex time series forecasting, this study introduced and compared four different data decomposition methods: VMD, EMD, EEMD, and CEEMDAN. Taking the original sea level data from the MAASSLUIS station as an example, these methods decomposed the data into multiple IMFs and a residual sequence. Subsequently, the decomposed sequences were used as the model’s features and individually fed into the LSTM model for making predictions. The results for each IMF and residual sequence are shown in Figure 4 and Figure 5. This experiment aimed to gain a deeper understanding of how the different data decomposition methods impact the performance of the LSTM model, and evaluated their potential for improving the accuracy of time series predictions.

In Figure 4 and Figure 5, the “Residual” presented for EMD, EEMD, and CEEMDAN refers to the results obtained by adding up the various IMFs after IMF5 and the residual sequence. From the figures, it can be observed that the IMFs obtained after VMD decomposition have well-defined frequency signals and waveform characteristics. Therefore, the LSTM model produced excellent predictions for each IMF. However, the residual sequence generated after VMD decomposition was relatively large and contained a significant amount of white noise. Consequently, even though there were some waveform features and patterns in the residual sequence, they were challenging for the LSTM model to capture, resulting in less accurate predictions, subsequently affecting the overall accuracy of the VMD–LSTM model’s predictions. In contrast, the EMD, EEMD and CEEMDAN methods, while not performing as well as VMD for predicting the various IMFs, yielded better prediction results for the residual sequence. To further analyze the accuracy of the predictions, this study summarized the evaluation metrics of each hybrid model's results, as shown in Table 4.

From the data in Table 4, it becomes apparent that the EEMD–LSTM model achieved the highest overall predictive accuracy, followed by the EMD–LSTM model and the CEEMDAN model, while the VMD–LSTM model exhibited the lowest predictive accuracy. However, it is noteworthy that a significant portion of the prediction errors in the VMD–LSTM model stemmed from the predictions of the residual sequence, and the prediction errors for the various IMFs were notably lower than those of the EMD–LSTM, EEMD–LSTM and CEEMDAN models.

Although the EEMD-LSTM model may have lower predictive accuracy for the residual sequence compared to the EMD-LSTM and CEEMDAN models, it excels in IMF prediction accuracy and overall accuracy. The CEEMDAN decomposition method, despite its enhanced robustness and applicability compared to EMD, yields predictive accuracy similar to that of the EMD-LSTM model. This indicates that the CEEMDAN-LSTM model does not significantly improve predictive performance in high-resolution sea level data compared to the EMD-LSTM model.

While the EEMD–LSTM model did not perform as strongly as the EMD–LSTM model in forecasting the residual sequence, it outperformed the EMD–LSTM model in forecasting the IMFs. As a result, the VMD–LSTM model excelled in IMFs prediction, whereas the EEMD–LSTM model exhibited the highest overall predictive accuracy. Building upon these insights, this study introduced the VMD–EEMD–LSTM model, which enhances overall predictive accuracy by reprocessing the residual components obtained from VMD decomposition with EEMD in addition to the VMD–LSTM model.

5. Discussion

5.1. Analysis of the predictions of the mixed VMD–EEMD–LSTM second-order decomposition model

To comprehensively evaluate the predictive performance of the VMD–EEMD–LSTM model relative to the VMD–LSTM and EEMD–LSTM models, this study conducted comparative experiments using sea level data from six different monitoring stations (Maasluis, Vlissingen, Hoek Van Holland, Delfzijl, Harlingen, Ijmuiden). In this section, Maaluis station is taken as an example to analyze the differences in the predictions of the hybrid models. To distinguish the models’ results more clearly, this section introduces the prediction error R and analyzes the differences between the original data and the predictions of each hybrid model. The comparative results are shown in Figure 6.

As depicted in Figure 6, the VMD–LSTM model, while reasonably aligning with the overall trend of sea level fluctuations, exhibited suboptimal performance near extreme points, particularly in proximity to local maxima. This observation suggests that the VMD–LSTM model struggled to capture the nuanced characteristics of sea level fluctuations, leading to notable prediction errors. In contrast, the EEMD-LSTM model's predictions closely match the original data, notably in capturing the amplitude of fluctuations, which significantly outperformed those of the VMD–LSTM model. Nevertheless, on a comprehensive scale, the results achieved by the EEMD–LSTM model still lagged behind those of the VMD–EEMD–LSTM model. This indicates that the VMD–EEMD–LSTM model not only represents an enhancement over the VMD–LSTM model but also surpasses the EEMD–LSTM model in predictive accuracy. It underscores the effectiveness of this hybrid model in combining the predictive strengths of the VMD–LSTM and EEMD–LSTM models, resulting in superior outcomes and overall improved predictive performance.

5.2. Analysis of the accuracy of the predictions of the mixed VMD–EEMD–LSTM second-order decomposition model

To gain a more precise insight into the enhancement achieved by the VMD–EEMD–LSTM model in comparison to the VMD–LSTM and EEMD–LSTM models across diverse time series, this section scrutinizes the RMSE, MAE, and R² of the predictions made by the three hybrid models for sea level time series data collected from six different stations. The evaluation metrics used to gauge the precision of each hybrid model's predictions are depicted in Figure 7, while the improvement ratios of the VMD–EEMD–LSTM model relative to the VMD–LSTM and EEMD–LSTM models for various evaluation metrics are detailed in Table 4.

Figure 7 clearly demonstrates that both the VMD-EEMD-LSTM model and the EEMD-LSTM model exhibit markedly superior predictive accuracy in comparison to the VMD-LSTM model. Furthermore, the VMD-EEMD-LSTM model showcases a noticeable degree of enhancement over the EEMD-LSTM model. The three hybrid models consistently demonstrated similar performance when predicting accuracy across various stations. This suggests that the sea level heights observed at the selected stations in the Netherlands displayed a degree of consistency, resulting in relatively minor variations in prediction accuracy. However, in comparison to the VMD–EEMD–LSTM model, the EEMD–LSTM model exhibited some fluctuations in the evaluation metrics across different time series predictions. This signifies that the stability and accuracy of the EEMD–LSTM model in forecasting results for diverse time series are not as robust as those of the VMD–EEMD–LSTM model. This result underlines the superiority of the VMD–EEMD–LSTM model in handling time series from different stations and, to some extent, validates its ability to adapt more stably to various requirements and scenarios of prediction.

To further evaluate the degree of improvement of the VMD-EEMD-LSTM model over the EEMD-LSTM model and the VMD-LSTM model across various accuracy assessment metrics, this study introduces the concept of an improvement ratio (I) for in-depth analysis. The results are presented in Table 5.

$I_{R^2}$$I_{R^2}$Table 5 provides clear evidence that the VMD–EEMD–LSTM model significantly improved prediction accuracy when compared to the EEMD–LSTM model across various stations. On average, it achieved a remarkable 26.95% reduction in RMSE, a 28.00% reduction in MAE, and a 6.53% increase in R². The EEMD–LSTM model showed a relatively modest increase of only 6.53% in the R², indicating that it could fit the actual distribution of the data well. The limited improvement in R² for the EEMD–LSTM model also indirectly confirmed the high predictive accuracy and superior performance of the VMD–EEMD–LSTM model.

Compared with the VMD–LSTM model, the VMD–EEMD–LSTM model exhibited even more significant improvements in the accuracy of its prediction, with an average reduction of 58.68% in the RMSE, an average reduction of 59.96% in the MAE, and an average increase of 49.85% in the R². This demonstrates that in practical VMD–LSTM predictions, there is significant room for improvement due to the incomplete decomposition of VMD.

In summary, the VMD–EEMD–LSTM model not only leverages the advantages of the LSTM model in handling long-term time series but also optimizes the variational decomposition of VMD and the adaptive iterative nature of EEMD. This results in the model demonstrating superior performance and producing better predictions in the field of time series forecasting of sea level height.

In summary, this study utilizes satellite altimetry data to estimate and forecast sea surface height. The findings indicate that the VMD-EEMD-LSTM model, which leverages the strengths of both hybrid prediction models, substantially enhances both predictive accuracy and the overall performance of sea surface height forecasts. Notably, it leads to significant improvements in forecasting the GSMSL, as evidenced by tests conducted along the Dutch coast.

6. Conclusion

This article discusses a new method for the high-precision time series forecasting of sea level height based on VMD–LSTM, named VMD–EEMD–LSTM. It addresses the limitations in the VMD–LSTM model, such as the insufficient decomposition of VMD, and enhances the robustness compared with the EEMD–LSTM model. The method's reliability was validated using multiple experiments involving Dutch coastal satellite altimetry data. The key findings are as follows.

(1) By comparing the predictions of different individual models, it is evident that the LSTM model exhibits the best predictive performance. However, the average RMSE remains high at 137.92 mm, the average MAE is 100.13 mm, and the average R² is only 0.40 across different measurement stations. This indicates that single deep learning predictive models often suffer from insufficient feature extraction when dealing with complex time series data, resulting in generally lower predictive accuracy.

(2) Comparing the four hybrid prediction models, VMD-LSTM, EMD-LSTM, EEMD-LSTM, and CEEMDAN-LSTM, the VMD-LSTM model has the lowest predictive accuracy across different measurement stations, with an average RMSE of 111.35 mm, an average MAE of 80.98 mm, and an average R² of 0.61. In contrast, the EEMD-LSTM model demonstrates the highest predictive accuracy, with an average RMSE of 63.82 mm, an average MAE of 45.71 mm, and an average R² of 0.87. Although the VMD-LSTM model lags EMD-LSTM EEMD-LSTM and CEEMDAN-LSTM models in overall predictive accuracy, its individual IMF components exhibit exceptionally high predictive accuracy within the LSTM model. While the IMF components of the EEMD-LSTM model may not match the VMD-LSTM model in predictive accuracy, the overall predictive accuracy of EEMD-LSTM surpasses that of VMD-LSTM.

(3) The VMD-EEMD-LSTM model, compared to the EEMD-LSTM model, achieves an average reduction of 26.95% in RMSE, an average reduction of 28.00% in MAE, and an average increase of 6.53% in R². Compared to the VMD-LSTM model, it achieves an average reduction of 58.68% in RMSE, an average reduction of 59.96% in MAE, and an average increase of 49.85% in R². These results illustrate that the VMD-EEMD-LSTM model, through the synergistic combination of the strengths from both hybrid prediction models, markedly improves both predictive accuracy and the overall performance of SSH forecasts.