An explainable hybrid deep learning framework for offshore wind speed forecasting

2.1. VMD and fuzzy entropy-based signal processing

To capture the complex, non-stationary, and nonlinear characteristics present in offshore wind speed data, this study adopts a two-stage preprocessing technique combining Variational Mode Decomposition (VMD) and Fuzzy Entropy (FE), hereafter referred to as VFE.

By applying VMD, $K$ intrinsic mode function (IMF) $\left\{u_k\right(t){\}}_{k=1}^K$ can be derived from the original wind speed signal $f\left(t\right)$ through adaptive decomposition, each with compact bandwidth around an adaptive center frequency $w_k$. To achieve accurate signal reconstruction, the variational formulation Eq. (1) aims to constrain the overall bandwidth of the extracted modes as much as possible [40]:

where, ${\partial }_t$ is the time derivative, $*$ denotes convolution, and $\delta \left(t\right)$ is Dirac function. To address this optimization task under constraints, a Lagrange multiplier $\lambda \left(t\right)$ and penalty term $\alpha$ are introduced, yielding the augmented Lagrangian:

The decomposition quality being highly sensitive to the penalty factor $\alpha$ and mode number $K$: inappropriate settings may lead to under- or over-decomposition, impacting model performance.

After decomposition, FE is employed to evaluate the complexity of each IMF. Unlike sample entropy, FE uses an exponential membership function to quantify the similarity between phase space vectors [41]. Given a time series $X=\{x_1,x_2,...,x_N\}$, the fuzzy similarity is defined as Eq. (3):

where $d_{ij}$ is the Chebyshev distance, $r$ is the similarity tolerance, and $n$ is the fuzzy power. The FE is computed by comparing the average similarity of vectors of dimension $m$ and $m+1$:

A higher FE value indicates greater irregularity and complexity. In this work, IMFs with higher FE are retained for prediction, while low-entropy components are discarded to reduce noise and computation burden.

By combining the adaptive decomposition capability of VMD and the complexity filtering power of FE, the VFE module enhances the signal quality and preserves meaningful temporal features for downstream prediction.

2.2. Ivy algorithm (IVYA)

Inspired by the adaptive growth and climbing patterns of ivy plants, the Ivy Algorithm (IVYA) was introduced as a meta-heuristic optimization strategy by Ghasemi et al. (2024) [42]. By simulating the biological mechanisms of ivy propagation – including coordinated expansion, adaptive search, and competition for sunlight – IVYA offers an effective strategy for solving complex, high-dimensional optimization problems.

IVYA models the ivy’s growth process using discrete-time dynamic equations and evolutionary rules. Its optimization framework consists of five major stages: (1) population initialization, (2) coordinated growth, (3) light-seeking movement, (4) expansion and adaptation, and (5) survivor selection. These components work together to guide the population toward global optima through both exploration and exploitation mechanisms. The IVYA’s process flow is depicted in Fig. 2.

Population initialization. Eq. (5) is employed to initialize the IVYA population with random positions distributed throughout the search space:

where, a vector with dimension $D$ of uniformly distributed random numbers is represented by $rand(1,D)$. $Npop$ is the number of populations.

Coordinated growth. It is postulated that the growth rate $Gv$ of ivy plants varies as a function of time. Drawing upon extensive experimental data, a corresponding difference equation describing the growth rate $G{v}_i\left(t\right)$ of individual members $I_i$ is derived:

where the vectors $\mathrm{\Delta }G{v}_i\left(t\right)$ and $\mathrm{\Delta }G{v}_i(t+1)$ represent the growth rate of a discrete-time system, $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}$ denotes a uniformly distributed real number within the range [0, 1], while $N(1,D)$ indicates a normally distributed random vector of dimension $D$.

Light-seeking movement. Young ivy is often guided to grow toward nearby trees or established ivy that has secured a support structure, allowing access to sunlight and promoting population sustainability.

The Eq. (7-8) describes how Ivy $I_i$ uses $I_{ii}$ to move logically in the direction of the radiant source:

Expansion and adaptation. Following a global search, member $I_i$ locates the closest and most influential neighbor $I_{ii}$ within the solution space, and $I_i$ tries to follow the best member of the population $I_{Best}$ to find the best candidate solution. This phase is represented by the Eq. (9):

Fig. 2Flowchart of the IVYA

Eq. (10) is applied to compute the present growth rate of member $I_i^{new}$, $\mathrm{\Delta }G{v}_i^{new}$ (which is exactly similar to the formula used to calculate $\mathrm{\Delta }G{v}_i$):

Survivor selection. In order to simulate the alternating stages, namely “climb” and “expand”, the method is used. If the objective value $f\left(I_i\right)$ for member $I_i$ falls below the product of $f\left(I_{best}\right)$ and parameter $\beta =\frac{(2+rand)}{2}$, Eq. (7) is used to expand the branch and leaf width of the ivy tree. Otherwise, the upward growth and climbing behavior of the ivy is guided by Eq. (9).

2.3. Convolutional neural network (CNN)

Due to their effectiveness in capturing local temporal patterns from structured inputs, Convolutional Neural Networks (CNNs) have become a prevalent choice for feature extraction tasks. The CNN architecture typically consists of three primary layers: the input layer, a series of hidden layers (including convolution and pooling operations), and the output layer [43, 44]. The CNN’s structure is shown as Fig. 3.

The core component of a CNN is its convolutional layer, which utilizes learnable filters to derive significant local characteristics from the input. This process is mathematically described in Eq. (11):

where $x_i^{(l-1)}$ denotes the feature map fed from the preceding layer, $w_{ij}^{\left(l\right)}$ is convolutional kernel, $b_j^{\left(l\right)}$ is the bias term, $*$ denotes convolution operation, $\sigma (\cdot )$ is the nonlinear activation function (commonly ReLU).

Following convolution, a pooling layer – specifically, max pooling – is applied to reduce dimensionality and suppress noise, preserving dominant features while mitigating overfitting. The pooling operation is defined as:

where $down(\cdot )$ represents the pooling function applied to the output feature map.

The final layer is typically fully connected, which aggregates extracted features and maps them to the target prediction.

This study uses CNN primarily to capture short-term temporal dependencies in preprocessed wind speed components, providing an effective input representation for subsequent temporal modeling.

Fig. 3Architecture of the CNN

2.4. Bidirectional gated recurrent unit (BiGRU)

To effectively model temporal relationships in sequential data, the BiGRU is employed. Compared to LSTM, the GRU architecture achieves similar performance while simplifying the gating mechanism by using fewer components [45, 46]. Its structure revolves around two essential components – the reset and update gates – which are visualized in Fig. 4(a).

The mathematical formulation of a GRU unit is as follows:

Eqs. (13-16) are related to each other and cannot be used alone. $r_t$ is reset gate. $z_t$ is update gate. $\widetilde{h_t}$ is candidate hidden layer state, reflecting the input information. $h_t$ is the output of the hidden layer. $\sigma$ and $\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}$ are the Sigmoid function and activation function, respectively; $W_r$, $U_r$, $W_z$, $U_z$, $W$, $U$ are all training parameter matrices.

BiGRU extends the GRU by processing the input sequence in both forward and backward directions, capturing richer temporal dependencies. As shown in Fig. 4(b), it consists of two parallel GRUs: one traversing the sequence from past to future, the other in reverse.

The output at time $t$ is computed as:

where $h_t^{\left(f\right)}$ and $h_t^{\left(b\right)}$ are the forward and backward hidden states, $W_f$, $W_b$ are the corresponding output weights, $b$is the bias term.

Compared with LSTM, BiGRU offers similar accuracy with fewer parameters and faster training convergence, thereby rendering it suitable for applications involving real-time forecasting. In this work, BiGRU is employed after CNN to model both forward- and backward-influenced temporal patterns in the offshore wind speed sequence, enhancing prediction robustness.

Fig. 4Architectures of a) GRU and b) BiGRU

a)

b)

2.5. SHAP Theory

SHAP (SHapley Additive exPlanations) is an explainable AI method based on cooperative game theory, designed to quantify the impact of each feature input on the prediction results produced by the model. It attributes the prediction of a model $f\left(x\right)$ to the sum of feature contributions, where each feature’s effect is measured by its Shapley value [47, 48].

For a given instance $x$, the model output $f\left(x\right)$ is defined using Eq. (18):

where ${\varphi }_0$ represents the expected output over the entire dataset, and ${\varphi }_i$ denotes the marginal contribution of the $i$th feature.

The Shapley value ${\varphi }_i$ of feature $i$ is defined as Eq. (19):

In this context, $F$ is the complete set of the features, while $S$ refers to a subset of $F$ that excludes the specific feature $i$. The $f\left(S\right)$ represents the model’s output based solely on the features contained in $S$.

This formula computes the average marginal contribution of feature $i$ over all possible subsets $S$, ensuring fairness and consistency in feature attribution. In practical implementation, SHAP approximates these values through sampling and model-specific explainers, such as TreeSHAP or KernelSHAP.

The key advantages of SHAP lie in its solid theoretical foundation and practical utility. It ensures consistency, meaning that if a feature has a greater impact on the output, its SHAP value will not be lower than that of a less influential feature. SHAP also provides local accuracy by ensuring that the sum of all feature attributions equals the difference between the model prediction and the expected output. Furthermore, by aggregating individual explanations across multiple instances, SHAP enables global interpretability, offering insights into overall feature importance and interaction effects, which is valuable for high-stakes applications like wind power forecasting.

3.1. Architecture of the CNN-BiGRU network

The CNN model alternates between convolutional and pooling layers, utilizing local connectivity and weight sharing to efficiently extract temporal features from the original signal. This process enables the network to automatically learn local patterns and form a dense, informative feature representation. Given its ability for automatic feature extraction and compatibility with sequential data, CNN is applied to capture temporal dependencies in offshore wind speed time series.

In the BiGRU component, the input data is processed bidirectionally using Gated Recurrent Units (GRUs). This allows the network to get temporal dependencies from both future and past contexts, enhancing its capacity to model the entire sequence structure. The reuse of weights across time steps enhances model expressiveness without increasing the number of parameters, which helps mitigate underfitting. Considering the volatility and uncertainty of offshore wind speed, combining CNN and BiGRU allows the model to achieve both efficient feature extraction and effective temporal modeling. The structure of the CNN-BiGRU model is illustrated in Fig. 5. The structure can be described as follows:

Input Layer: The input is structured as a two-dimensional matrix of size, where corresponds to the temporal dimension of the sequence, and indicates the feature dimensionality per time step. The CNN module first processes this data to extract temporal patterns.

Convolutional Layer: 1D convolutional layers are used to extract temporal features. Each has a kernel size of 2 with a stride of 1, and the number of filters is set to 32 and 64 respectively. ReLU is used as the activation function, introducing non-linearity and accelerating convergence by zeroing out negative activations.

Pooling Layer: A max-pooling layer follows each convolutional layer to reduce dimensionality and eliminate less informative features. The pooling size is 2×1 with a stride of 2. The output of the pooling layers is flattened into a time sequence format suitable for BiGRU processing.

BiGRU Layer: The extracted features are then input into the BiGRU network, which processes the sequence bidirectionally to capture comprehensive temporal dependencies. A Dropout layer is added within the BiGRU to prevent overfitting by randomly deactivating some units during training. The final output is passed through a fully connected layer to generate the predicted wind speed.

Optimization and Training: The Adam optimizer is employed for training due to its adaptive learning rate capabilities. The number of training epochs is set between 1 and 200, and the initial learning rate ranges from 0.001 to 0.01. A learning rate decay factor of 0.1 is applied to ensure stable convergence. The configurations are chosen to achieve a trade-off between training speed and the model’s ability to generalize.

The specific configurations of the CNN-BiGRU model, including convolutional kernel sizes, number of filters, dropout rate, and optimization settings, are detailed in Table 1. These settings are based on standard practices in deep learning for time series modeling and have been refined to ensure effective temporal feature extraction, training stability, and model generalization. Particular attention has been paid to adapting the convolution and pooling parameters to the 1D structure of time series data, as well as to selecting appropriate ranges for hyperparameter tuning.

Fig. 5Architecture of the CNN-BiGRU

3.2. Forecasting model formulation

The modeling process of the proposed VFE-IVYA-CNN-BiGRU framework for offshore wind speed forecasting consists of three main stages, as illustrated in Fig. 6. The steps are detailed as follows:

(1) Data Acquisition and Preprocessing.

Historical offshore wind speed data are collected on a daily basis. The raw dataset is first cleaned by removing zero-value records and interpolating missing entries [49]. To avoid information leakage, the dataset is then chronologically partitioned into training and testing subsets, after which all preprocessing procedures are performed. Subsequently, min-max normalization is fitted on the training subset, and the learned scaling parameters are applied to the testing data to ensure proper time-ordered validation, as shown in Eq. (20):

where $x$ denotes the original data; $x_{max}$ and $x_{min}$ represent the maximum and minimum values of $x$, respectively, while $x^{\text{'}}$ refers to the normalized result after transformation.

In addition to normalization, Pearson Correlation (PCs) is used to measure linear relationships among features. This step assists in filtering relevant meteorological variables before entering the VFE module while minimizing the risk of collinearity-driven leakage across target-related predictors.

Fig. 6Flowchart of the proposed VFE-IVYA-CNN-BiGRU framework

(2) Signal Decomposition and Feature Processing.

Step 2.1: For the training subset, the normalized wind speed series is decomposed using the VFE module, which combines VMD and FE. VMD adaptively separates the signal into several IMFs, each capturing specific frequency content and a residual (RES), helping to reduce signal non-stationarity. The decomposition parameters derived from the training set are then applied to process the testing subset to ensure consistency and prevent data leakage.

Step 2.2: FE is computed for each IMF component in the training subset to assess its complexity. Based on the entropy values, IMFs are reconstructed to form a new sequence that emphasizes informative structures and suppresses noise-dominated components. The reconstruction scheme is applied consistently to the testing subset using the thresholds determined from the training data.

Step 2.3: CNN is employed to extract localized temporal features, producing a high-dimensional feature matrix. This matrix is then input into the BiGRU layer to capture temporal dependencies from both forward and backward directions.

(3) Model Training and Prediction.

Step 3.1: Key hyperparameters of the CNN-BiGRU architecture are optimized using the Ivy Algorithm (IVYA), including the number of convolutional filters, learning rate, batch size, and number of hidden units.

Step 3.2: With the optimized hyperparameters, the feature-enriched dataset is used to train the CNN-BiGRU architecture. The trained model is then used to predict wind speed values on the test set. Each reconstructed component is predicted independently.

Step 3.3: Finally, the predicted sub-sequences are aggregated and then inverse normalized to restore the original scaleto yield the complete offshore wind speed forecast. The model performance is evaluated using standard error metrics, and results are compared against benchmark models.

3.3. Evaluation metrics

To assess how closely the VFE-IVYA-CNN-BiGRU model’s predictions align with actual wind speed values, this study applies four commonly used evaluation metrics: R², RMSE, MAE, and MAPE [50]. These indicators help in analyzing the predictive reliability and error magnitude of the model. The expression is shown as Eq. (21-24):

where $p_i$ is the true value of the $i$th sample; $p_i^{\text{'}}$ is the predicted value of the $i$th sample; ${\bar{p}}_l$ is the average of all the samples; $N$ represents the number of predicted samples.

4.1. Dataset description

This study employs data collected from an offshore wind farm in Guangdong Province, China, during the one-year period spanning from the beginning to the end of 2022 (January 1 to December 31). Wind speed and associated meteorological variables were recorded at 30-minute intervals, forming the basis for seasonal prediction analysis. The dataset includes ten meteorological factors: wind direction (WDIR), gust speed (GSP), wave height (WVHT), dew point temperature (DPT), dominant wave period (DWP), sea level pressure (SLPR), air temperature (ATMP), sea surface temperature (SST), average wave period (AWP), and actual wind speed (AWS).

For seasonal evaluation, the full-year dataset is divided into four subsets corresponding to distinct time frames: Winter (January 1-31), Spring (April 1-30), Summer (July 1-31), Autumn (October 1-31). Each seasonal subset is further split into training (70 %) and testing (30 %) sets. The number of samples per season varies slightly due to the differing lengths of each month: 1488 (Winter), 1436 (Spring), 1474 (Summer), and 1484 (Autumn). Data preprocessing includes zero-value removal, interpolation for missing values, and min-max normalization.

The seasonal variations in offshore wind speed are illustrated in Fig. 7, where the temporal trends across the four seasons show distinct fluctuations and amplitudes, reflecting meteorological seasonality and variability. Table 2 provides a statistical overview of wind speed across different seasons, reporting values such as the mean, median, range (maximum and minimum), standard deviation, skewness, and kurtosis.

All experiments were conducted using a computing platform operating under Windows 11, powered by an Intel Core i9-12900H CPU and supported with 16 GB of RAM. The implementation of the model was carried out using Python 3.10, with TensorFlow 2.14 serving as the underlying deep neural modeling toolkit. All baseline and proposed models were trained under identical training budgets, input window sizes, batch settings, optimizer schedules, and stopping criteria to ensure a fair comparison. Computational metrics such as convergence curves and latency are not included to maintain focus on forecasting performance and interpretability.

Fig. 7Original offshore wind speed data in four seasons

4.2. Meteorological feature analysis

Offshore wind speed is influenced by nine marine meteorological parameters, including WDIR, GSP, WVHT, DWP, AWP, SLPR, ATMP, SST and DPT. To quantitatively assess the relationships between these factors and wind speed, Pearson Correlation (PCs) analysis is employed.

The PCs between each meteorological variable and the actual offshore wind speed is calculated by Eq. (25):

where the correlation coefficient $r$ ranges from [–1, 1], when it is greater than 0, the stronger the correlation, the larger the value; $n$ is the total number of dataset; $x$ is the meteorological factor data; $y$ is the offshore wind speed data [51].

Fig. 8Feature correlation heatmap (Winter)

The correlation heatmap of the winter dataset is visualized in Fig. 8, clearly showing that gust speed (GSP) exhibits the highest positive correlation with wind speed, followed by SLPR and DPT. To generalize this analysis, Table 3 presents the PCs values between wind speed and all nine features across the four seasonal datasets.

Based on the PCs results, the top six variables with the highest absolute correlation values were selected for each seasonal subset as input features for the prediction model. These selected variables are listed in Table 4. Notably, GSP consistently ranks highest across all seasons, indicating its dominant influence on wind speed variation. Other features, such as SLPR, DPT, and WDIR, also show moderate and seasonally consistent correlations.

4.3. Wind speed decomposition via VMD

In order to improve the stability and predictability of the offshore wind speed time series, VMD serves to break down the raw input sequence into a set of IMFs, each representing distinct oscillatory patterns embedded in the original data. A critical parameter in VMD is the decomposition level $K$, which determines the number of modes. If $K$ is set too low, under-decomposition may occur, failing to capture the full signal complexity; conversely, an excessively high $K$ may result in mode aliasing.

To determine an appropriate value of $K$, the center frequency method is adopted. Winter dataset is exemplified, Table 5 presents the center frequencies of the IMFs for $K=\text{2}$ to $K=\text{9}$. When $K=\text{8}$ and $K=\text{9}$, the highest center frequency stabilizes at 0.463, indicating convergence. Hence, the optimal decomposition level is set to $K=\text{8}$.

The results of VMD decomposition for winter, spring, summer, and autumn are shown in Fig. 9. The offshore wind speed signals are effectively decomposed into low- and high-frequency components, which serve as the basis for further reconstruction and prediction.

Fig. 9VMD-based decomposition of wind speed series in four seasons

4.4. Wind speed reconstruction via fuzzy entropy

After VMD decomposition, the offshore wind speed time series is represented as a set of IMF components. Directly predicting all components may lead to increased computational complexity and error accumulation. To address this, the Fuzzy Entropy (FE) of each component is calculated to quantify its complexity. Components with similar FE values are grouped and reconstructed, reducing redundancy while preserving key dynamic patterns.

Based on the FE evaluation, IMFs with similar entropy levels were grouped to form low-, mid- and high-complexity components. Lower-entropy IMFs typically exhibited more regular oscillations, whereas higher-entropy IMFs contained richer dynamical variations and contributed more predictive information. This FE-guided reconstruction ensures that the retained components reflect distinct temporal behaviors while suppressing noise-dominated signals.

Fig. 10Fuzzy entropy analysis results for seasonal datasets

The FE values and seasonal differences are visualized in Fig. 10, which shows that complexity varies by season and IMF component. Taking the winter dataset as an example, this reconstruction result illustrated in Fig. 11, enables the generation of three clear sub-sequences with distinguishable frequency characteristics. Compared with the original full IMF set, the reconstructed sequences offer lower dimensionality and improved interpretability, which contributes to more robust and accurate prediction. It should be noted that the fuzzy entropy (FE) values of the IMFs differ across the four seasonal datasets. This variation arises because wind-speed time series exhibit distinct turbulence intensities and fluctuation characteristics under different seasonal meteorological conditions. Since FE measures the irregularity and complexity of a signal, the resulting entropy values naturally reflect these seasonal dynamics rather than remaining uniform across all datasets.

4.5. Comparison of forecasting techniques

To thoroughly evaluate the proposed VFE-IVYA-CNN-BiGRU model, this study conducts a comparative assessment against five benchmark forecasting models: BiGRU, TCN, Transformer, CNN-BiGRU, IVYA-CNN-BiGRU, VMD-IVYA-CNN-BiGRU and VFE-IVYA-CNN-Transformer. The models are tested on four seasonal offshore wind speed datasets – winter, spring, summer, and autumn – and evaluated using MAE, RMSE, MAPE, and the coefficient of determination $R^2$. Quantitative results are reported in Table 6 and visualized in Figs. 12–19.

Fig. 11Reconstructed wind speed signal after FE-based IMF filtering (Winter)

Among the four seasonal datasets, the winter results most clearly demonstrate the performance advantages of the proposed VFE-IVYA-CNN-BiGRU model. Winter wind speed series exhibit strong periodicity coupled with moderate low-frequency fluctuations, making this season an essential benchmark for evaluating multi-scale feature extraction and temporal-dependency modeling. On this dataset, the proposed model achieves an $R^2$ exceeding 0.97, significantly outperforming BiGRU (≈ 0.8860), TCN (≈ 0.8803), Transformer (≈ 0.8753), and CNN-BiGRU (≈ 0.9037), corresponding to a relative improvement of approximately 8.9-11.3 %. The reduction in error metrics is even more pronounced: compared with the basic BiGRU, the proposed architecture reduces MAE by more than 55 % and RMSE by nearly 60 %, indicating substantially enhanced accuracy in capturing both the amplitude variation and temporal evolution of winter wind speed. Even against stronger convolution-based baselines such as CNN-BiGRU and IVYA-CNN-BiGRU, MAE and RMSE still decrease by 40-50 % and 45-55%, respectively, demonstrating the synergistic benefits introduced by VFE decomposition and IVYA-enhanced convolution for representing localized structures and high-frequency fluctuations. Furthermore, relative to more advanced hybrid approaches – such as VMD-IVYA-CNN-BiGRU and VFE-IVYA-CNN-Transformer – the proposed model maintains clear superiority, achieving 15-25 % reductions in RMSE and 20-35 % reductions in MAPE. These improvements further highlight the decisive role of the BiGRU component in modeling long-range dependencies within winter’s quasi-stationary wind patterns.

The winter ablation results show a clear stepwise improvement as each module is added. Starting from the BiGRU baseline ($R^2\approx$ 0.8860), introducing CNN reduces MAE and RMSE by roughly 20-30 %, while replacing CNN with IVYA brings a further 10-15 % reduction. Adding VMD provides another 15–20% improvement in RMSE through more effective mode separation. With VFE incorporated, the model reaches $R^2>$ 0.97 and achieves over 55 % and 60 % decreases in MAE and RMSE compared with BiGRU, confirming that each component contributes incremental and complementary gains.

By contrast, the spring, summer, and autumn datasets exhibit stronger irregularity and transitional seasonal characteristics, resulting in relatively narrower performance gaps among models. Nevertheless, the proposed method remains the top performer across all evaluation metrics, typically reducing RMSE by 10-20 % and MAPE by 10-25 % compared with the next-best model. In particular, for the highly non-stationary and fluctuation-intensive summer dataset, the proposed model improves $R^2$ by more than 30 % relative to BiGRU, confirming its strong adaptability under chaotic wind conditions. Overall, the winter experiments provide the clearest quantitative evidence of the method’ s core advantages: VFE ensures high-fidelity mode separation and mitigates modal mixing; IVYA-CNN enhances the extraction of local and high-frequency structures; and BiGRU delivers robust long-term dependency modeling. The synergistic integration of these components enables the proposed model to achieve the highest forecasting accuracy across all seasons, with winter exhibiting particularly substantial improvements.

Fig. 12Results from six forecasting approaches (Winter)

Fig. 13Results from six forecasting approaches (Spring)

These statistical gains are corroborated by the time-series prediction plots (Figs. 12-15), where the proposed model demonstrates significantly better temporal alignment with ground truth, especially around peaks, troughs, and inflection zones. In winter and autumn, where daily amplitude varies widely, the model accurately follows trend direction without lag – whereas other models, especially BiGRU and CNN-BiGRU, show smoothing effects or phase shift. In summer, the proposed model shows reduced overfitting and higher responsiveness to random fluctuations, which is particularly evident in fast-changing segments.

Fig. 14Results from six forecasting approaches (Summer)

Fig. 15Results from six forecasting approaches (Autumn)

Visual evaluations from radar plots (Figs. 16-19) further emphasize the model's advantage. The VFE-IVYA-CNN-BiGRU consistently occupies the innermost region, indicating minimum error across all axes (MAE, RMSE, MAPE). In bar plots, the percentage decrease in MAE reaches up to 57.5 %, and MAPE up to 55.1 %, compared to the weakest models. Moreover, scatter-fitting plots reveal denser clustering near the identity line with a steeper slope closer to unity and reduced spread, reflecting both reduced bias and variance. The residual distribution is significantly narrower, especially in high-variance seasons, further confirming better model reliability.

The consistent improvement of the proposed model is attributed to the synergistic integration of three key components. The original signal is decomposed into a set of intrinsic mode functions (IMFs) using VMD, isolating multiscale oscillations and reducing non-stationarity. Fuzzy Entropy quantifies the complexity of each IMF, allowing the model to suppress low-informative or noise-prone components. This process produces a reconstructed sequence with a higher signal-to-noise ratio, which enhances learning stability. Finally, IVYA adaptively optimizes the model’s architecture – including convolutional depth, recurrent capacity, and learning parameters – resulting in efficient convergence and enhanced generalization across seasonal regimes.

In summary, the VFE-IVYA-CNN-BiGRU model exhibits not only lower average forecasting error but also superior temporal alignment, reduced sensitivity to extreme values, and better fit in noisy or high-variance conditions. Its performance advantage is both statistically significant and structurally grounded, enabling robust offshore wind speed forecasting across all seasonal contexts.

Fig. 16Comparative forecasting accuracy of six models (Winter)

Fig. 17Comparative forecasting accuracy of six models (Spring)

Fig. 18Comparative forecasting accuracy of six models (Summer)

Fig. 19Comparative forecasting accuracy of six models (Autumn)

To further rule out the possibility of target leakage or artificially inflated performance caused by highly collinear meteorological variables, an additional ablation experiment was conducted on the winter dataset using the proposed VFE-IVYA-CNN-BiGRU model. In this test, the gust-speed variable (GSP), which exhibits the strongest correlation with wind speed, was removed from the input set. Excluding GSP led to a clear degradation in forecasting accuracy, with $R^2$ decreasing to 0.9322, MAE increasing to 0.4729 m/s, RMSE rising to 0.6536 m/s, and MAPE increasing to 7.864 %. This consistent deterioration confirms that the gust-related variable provides genuine predictive information rather than leaking future values or duplicating the target signal. The results further demonstrate that the proposed framework effectively captures physically meaningful meteorological drivers, and that the performance gains attributed to gust-related features arise from true model learning rather than from spurious correlations or unintended data leakage.