Prediction of sea ice area based on the CEEMDAN-SO-BiLSTM model

Complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN)

Due to the “mode mixing” caused by EMD and the noise residual caused by EEMD, this article introduces complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN), which overcomes the defects of EEMD decomposition in terms of loss of completeness and mode mixing by adaptively adding white noise. In the algorithm, E_i(⋅) is defined as the i-th mode generated by EMD decomposition, $\overline{Ci\left(t\right)}$represents the ith mode generated by CEEMDAN decomposition, ɛ is the standard deviation of the noise, v^j follows N (0,1) and j =1,2, ….,N denotes the number of times white noise is added, while r represents the residue. The specific steps of the CEEMDAN algorithm are as follows:

1. Add Gaussian white noise to the original signal y(t) to obtain a new signal y(t) + (−1)^qɛv^j, where q = 1,2. EMD is performed on the signal to obtain the first-stage intrinsic mode component C₁ (1)${\displaystyle E\left(y\left(t\right)+{\left(-1\right)}^q\varepsilon v^j\left(t\right)\right)=C_1^j\left(t\right)+r^j}$

2. Taking the overall average of the N generated mode components produces the first intrinsic mode function of the CEEMDAN decomposition. (2)${\displaystyle \overline{C_1\left(t\right)}=\frac{1}{N}\sum _{j=1}^N{C}_1^j\left(t\right)}$

3. Calculate the residue after removing the first modal component (3)${\displaystyle r_1\left(t\right)=y\left(t\right)-\overline{C_1\left(t\right)}}$

4. Add paired positive and negative Gaussian white noise to r ₁(t) to obtain a new signal, and perform EMD on the new signal to obtain the first-order modal component D ₁. Then, the second intrinsic mode component of the CEEMDAN decomposition can be obtained. (4)${\displaystyle \overline{C_2\left(t\right)}=\frac{1}{N}\sum _{j=1}^N{D}_1^j\left(t\right)}$

5. Calculate the residue after removing the second modal component (5)${\displaystyle r_2\left(t\right)=r_1\left(t\right)-\overline{C_2\left(t\right)}}$

6. The above steps are repeated until the obtained residual signal is a monotonic function and cannot be further decomposed, and the algorithm ends. The number of intrinsic mode functions obtained at this time is denoted as K, and the original signal y(t) can be decomposed as: (6)${\displaystyle y\left(t\right)=\sum _{k=1}^K\overline{C_k\left(t\right)}+r_k\left(t\right)}$

Bidirectional long short-term memory neural network (BiLSTM) and optimization

BiLSTM

The BiLSTM allows for the transmission and feedback of past and future states of the hidden layers through a bidirectional network (Fig. 1).

The BiLSTM is interpreted through Eq. (7). (7)${\displaystyle \left\{\begin{array}{c}{h}_f=\text{LSTM}\left(x_i,h_{f-1}\right)\hfill \\ h_b=\text{LSTM}\left(x_t,h_{b-1}\right)\hfill \\ h_t=w_t{h}_f+v_t{h}_b+b_t\hfill \end{array}\right.}$ where x_i is the input, h_f is the forward-pass implicit layer state, h_b is the reverse-pass implicit layer state, h_t is the implicit layer state, w_t is the forward-pass implicit layer output weight, v_t is the reverse-pass implicit layer output weight, and b_t is the error value.

BiLSTM for snake optimizer optimization (SO-BiLSTM)

Based on the special mating behavior of snakes, Hashim & Hussien (2022) proposed the snake optimizer (SO). The algorithm is divided into two stages: global exploration when there is no food and local exploitation when there is food. When food is scarce (Q<0.25), the snakes search for food and update the location by choosing any random location to achieve global optimization. When there is sufficient food (Q ≥0.25), snakes move toward food and update their position if Temp>0.6; if the temperature Temp ≤0.6, the snakes enter combat mode at a random number Rand<0.6 taking the value of [0,1] and enter mating mode at Rand ≥0.6 and replace the individual snake with the worst fitness value by laying and hatching eggs.

The prediction accuracy of the BiLSTM neural network is affected by the number of neurons in the hidden layer, the learning rate and the L2 regularization coefficient. Considering the randomness of artificially set parameters, the root mean square error RMSE is selected as the fitness function in this article, and the parameters are optimized using the SO algorithm (Fig. 2).

CEEMDAN-SO-BiLSTM

Through the above analysis, this article introduces the CEEMDAN module for noise reduction of the original data based on the advantages of the optimization algorithm and deep learning and proposes a combined CEEMDAN-SO-BiLSTM prediction model (Fig. 3).

As seen from the above figure, the prediction model in this article is divided into three parts. The first part is CEEMDAN decomposition, which decomposes the time series into K modes for noise reduction; the second part is the SO-BiLSTM neural network model to train the prediction of K modes; and the third part superimposes and reconstructs the prediction results of the second part to obtain the prediction results of the original data.

Data sources

The data set used in this article was obtained from the sea ice data set on the official website of the National Environmental Information Center (https://www.ncei.noaa.gov/) and the Grambling Sea Ice area time series data was used for the study in the model analysis. The time range is from day 272 of 2017 to day 271 of 2022, with a total of 1,824 data points.

Model evaluation criteria

To evaluate the prediction performance of the model, this article selects four indicators: mean absolute error (MAE), root mean square error (RMSE), determinability coefficient R², and symmetric mean absolute percentage error (SMAPE). Equations for calculating these indicators are provided in Eqs. (9)–(12).

These indicators quantify the accuracy of the model’s predictions. MAE and RMSE represent the average and root of the squared errors between the predicted and actual values, respectively, while R² measures the correlation between the predicted and actual values. SMAPE calculates the symmetric percentage difference between the predicted and actual values. (9)${\displaystyle \mathrm{MAE}=\frac{1}{n}\sum _{i=1}^n\left|y_i-{\hat{y}}_i\right|}$ (10)${\displaystyle \text{RMSE}=\sqrt{\frac{1}{n}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left(y_i-{\hat{y}}_i\right)}^2}}$ (11)${\displaystyle {\mathrm{R}}^2=\frac{\sum _{i=1}^n{\left(y_i-{\hat{y}}_i\right)}^2}{\sum _{\mathrm{i}=1}^{\mathrm{n}}{\left(y_i-\overline{y}\right)}^2}}$ (12)${\displaystyle \text{SMAPE}=\frac{100\text{\%}}{n}\sum _{i=1}^n\frac{\left|{\hat{y}}_i-y_i\right|}{\left(\left|{\hat{y}}_i\right|+\left|y_i\right|\right)/2}}$ where y_i is the i-th observation, ${\hat{y}}_i$ is the i-th prediction, $\overline{y}$ is the mean, and n is the number of samples.

Experimental analysis

According to the characteristics of the data sets, this article first partitions the original data into a training set comprising 70% of the data and a testing set consisting of the remaining 30%. Then, using the rolling window method with CEEMDAN decomposition, it predicts the values within each roll-out window. The noise ratio is set at 0.2, with 500 iterations of adding noise to the signal and a maximum span of 2000. For example, it shows the different modalities obtained by decomposing the first prediction window (Fig. 4), including 10 modes representing various frequency dimensions of the ice concentration time series. The high correlation coefficients of IMF8 and IMF9 with the original sequence, reaching 0.91 and 0.77, respectively, indicate that these two modes contribute significantly to the periodic trend of the original signal. IMF1 to IMF3 are considered as noise signals, while IMF10 represents a short-term trend component (Fig. 5). As seen from the figure, compared to the original sequence, the decomposed modalities are more stable and smooth and exhibit clear information features, providing a solid foundation for predictions.

The specific process of making predictions using the rolling window method is as follows: assuming that the model training and test time series data are D, with length T + k, where the first T data points are used for model training and the last k data points are used for model testing. It uses the previous 30 data points to predict the next data point in the testing set. To obtain each predicted value from the testing set, follow these steps:

1. CEEMDAN decomposition is performed on the first T data points to extract multiple features and train separate models for each feature.

2. Following completion of model training, utilize data from D[T-29:T] to obtain predicted outputs from each feature model and sum them to acquire the prediction value for D[T+1].

3. Execute CEEMDAN1 decomposition on data D[2:T+1] to extract multiple features once again and train individual models for every feature. Utilizing 30 data points from D[T-28:T+1], it receives the predicted data from each feature model and aggregate them to attain the forecast for D[T+2].

Repeat actions 1 through 3, sliding the window of 30 data points and performing decomposition and model training on the enclosed data before making predictions, ultimately arriving at projections for all future time intervals D[T+1:T+k].

After applying SO-BiLSTM training for predictive outcomes from the preceding splitting apart, SO adjusts BiLSTM’s primary acquisition pace, opaque node contingent, and L2 regulation factor. Following numerous trials, the versatile extent characteristic ultimately slopes toward equilibrium (Fig. 6), and the versions for every mode within the initial estimation aperture are improved as per Table 1.

For each modality, the subsequent predictive value is recreated so that the concluding forecast can be generated for the initial information. By persistently pushing the prognostication casement ahead and producing supplementary projections, up to thirty percent of all estimations have been accomplished (Fig. 7). Evaluation parameters for the objective model are determined using Table 2.

Model comparison

To verify the superiority of the model proposed in this article, models such as BiGRU and BiLSTM were tested and compared with the model in this article, and the test results of each model are shown in Table 3.

Based on our findings (Fig. 8), it appears that the performance of the comparison model deteriorates during specific periods, particularly around May 2021 and January 2022, characterized by substantial fluctuations and increased prediction bias. Notably, the ARIMA model demonstrates a remarkable shortcoming in fitting the original time series during these periods.

On the other hand, the proposed CEEMDAN-SO-BiLSTM model yields a superior fitting capability, effectively capturing the variability of the time series. As evidenced in Tables 2 and 3, the proposed method exhibits remarkable advantages over both the single-model counterparts and the benchmark SVR model. Specifically, the combined CEEMDAN-SO-BiLSTM model reduces MAE, RMSE, and SMAPE by approximately 81.8%, 81.5%, and 9.150%, respectively, and raises R² by approximately 3.9%.

Furthermore, our experimental findings reveal that applying decomposition and optimization strategies consistently enhances prediction accuracy, with the CEEMDAN-SO-BiLSTM model outperforming its VMD counterpart. Finally, among the different variants of the LSTM and BiLSTM models, CEEMDAN-SO-BiLSTM provides the best performance. These results underscore the efficacy and robustness of our proposed framework in accurately predicting sea ice areas in the Greenland Sea.