The Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) is an advanced signal processing method designed for analyzing nonlinear and non-stationary time series. It extends the original Empirical Mode Decomposition (EMD) by improving the decomposition stability and mitigating mode mixing—an issue in traditional EMD where intrinsic mode functions (IMFs) with overlapping frequency content compromise the interpretability of results (Guo et al., 2023).

CEEMDAN enhances EMD by introducing Gaussian white noise into the decomposition process. Initially, Gaussian noise is added to the original signal to generate multiple noisy versions. Each noisy signal is decomposed using EMD to extract its IMFs. Subsequently, the corresponding IMFs from each realization are averaged to obtain a more robust and stable decomposition. The incorporation of adaptive noise further ensures that each IMF retains high independence and spectral purity, thereby improving decomposition accuracy. This process enables effective noise reduction while preserving the underlying characteristics of the original signal (Hu et al., 2021). In this study, the IMFs and residual components extracted via CEEMDAN serve as refined input features for subsequent spatial-temporal modeling.

Figure 1 illustrates the basic workflow of Graph Neural Networks (GNN), whose main objective is to extract spatial features from multivariate data at each time step and model the dependencies between variables. Specifically, a graph structure is used to represent the relationships between variables:

Nodes represent the 18 variables, including AQI and major pollutant concentrations (e.g., PM_2.5, PM₁₀, SO₂, CO, NO₂, O₃_8 h), as well as the IMFs (IMF_1 to IMF_10) and Residue from CEEMDAN decomposition.

Edges define the connections between the variables through an adjacency matrix.

Graph Convolutional Networks (GCN), a common variant of GNN, are used to update the spatial features, and the update process can be represented by the following equations, as shown in Equations 1–3 (Zhou et al., 2020): H t = σ A ∼ X t W + b (1)

Where X t ∈ R N × F represents the input feature at time t , A ∼ ∈ R N × N is the normalized adjacency matrix. W ∈ R F × D is the learnable weight matrix, where D is the hidden layer dimension. b ∈ R D represents the bias term. σ is the activation function, representing spatial features at time t .

To enhance expressiveness, multiple layers can be stacked, allowing the network to capture more complex relationships. H t l = σ A ∼ H t l − 1 W l + b l (2) H t l = X t is the initial input, and H t l is the output of the L th layer, serving as the final spatial feature. The GNN is applied across all time steps, producing the output as: H = H 1 , H 2 , … , H T ∈ R T × N × F (3)

This process is applied at every time step, yielding spatial features across time.

The spatial features extracted by GNN need to be fused with the original input features to generate a unified input for the Transformer model.

To fully leverage the spatial feature extraction capabilities of GNN, the spatial features are fused with the original input features to form a combined feature suitable for Transformer modeling. The fusion process is carried out as shown in Equations 4–6:

Concatenation: The features are concatenated. Z t = X t , H t ∈ R N × F + D (4)

Linear Transformation: A linear transformation is applied. Z t = W f X t , H t + b f (5)

Attention Fusion: Multi-head attention mechanisms dynamically fuse the features. Z t = M u l t i H e a d A t t e n t i o n Q = X t , K = H t , V = H t (6)

Where Q (queries) are derived from original features, and K ((keys) and V (values) are obtained from GNN features. The output Z t ∈ R N × D ` is reshaped into Z ∈ R T × N × D ` , facilitating dynamic weighted fusion and improving sensitivity to key variables.

Figure 2 illustrates the main structure of the Transformer model. The Transformer model receives the fused feature Z , which is used to model the temporal dependencies in the time series for future prediction. The input processing first reshapes Z to match the input format required by the Transformer. Then, position encoding is applied to add temporal information to the input sequence, as Transformer models do not inherently account for sequence order (Vaswani et al., 2017). Position encoding typically uses sinusoidal functions, as shown in Equations 7, 8: PE t , 2 i = s in t 1000 2 i d model (7) PE t , 2 i + 1 = cos t 1000 2 i d model (8)

In the attention mechanism, the input sequence is projected into queries ( Q ), keys ( K ), and values ( V ), where W Q , W K , and W V are trainable projection matrices. The attention score is computed by taking the dot product of the query and key, followed by scaling and applying the Softmax function, as shown in Equation 9: Attention Q , K , V = Softmax Q · K t d k · V (9)

Multi-head attention involves parallelizing multiple self-attention heads to capture different feature subspaces. For each head h , the computation can be represented as shown in Equation 10: h ead h = Attention Q h , K h , V h (10)

The outputs of all heads are concatenated and passed through a linear transformation, as shown in Equation 11: MultiHead Q , K , V = Concat h ead 1 , … , h ead h · W O (11)

The attention output is then passed through a Feed-Forward Network (FFN) consisting of two linear layers with a ReLU activation, as shown in Equation 12: FFN x = max 0 , x · W 1 + b 1 · W 2 + b 2 (12)

Residual connections and Layer Normalization are used after each sub-layer to ensure training stability, as shown in Equation 13: LayerNorm x + SubLayer x (13)

For time series forecasting, the output of the final Transformer layer is typically the predicted sequence, representing future values at each time step. If Z out is the output of the final Transformer layer, it can be mapped to the desired prediction dimensions through a linear layer, as shown in Equation 14: y ^ = Z out · W y + b y (14)

Figure 3 presents the overall architecture of the proposed CEEMDAN-GNN-Transformer hybrid framework. The process begins with data collection and preprocessing. The target air quality variable (AQI) is decomposed using CEEMDAN into IMFs and a residual component to reduce non-stationarity and noise.

An ablation study is conducted to verify the effectiveness of each model component, and hyperparameters are tuned to achieve optimal performance.

Step 1

The GNN component extracts spatial dependencies among multiple pollutants and decomposed features at each time step, revealing inter-variable relationships.

Step 2

These spatial features are fused with the original features through a designed fusion mechanism, yielding representations suitable for sequence modeling.

Step 3

The Transformer module captures long-term temporal dependencies and nonlinear patterns within the air quality time series, ultimately producing forecasts of future AQI values.

This hybrid framework effectively combines multi-scale signal decomposition, spatial graph modeling, and temporal sequence learning, offering a robust and interpretable approach to air quality forecasting, particularly under the complex conditions present in industrialized urban environments like Chang’an Town.

This study selects Chang’an Town in Dongguan City as the research area, as shown in Figure 4, for the empirical analysis of air quality time series forecasting. Chang’an Town is one of the most economically developed areas in Dongguan, covering a wide range of industries, including manufacturing, technology, and others. The area has a high concentration of manufacturing enterprises, which results in significant air pollution, thus necessitating precise air quality prediction methods to support environmental management and pollution control.

The dataset spans from 2015 to 2024 and includes 3,647 daily records of major air pollutants: PM_2.5, PM₁₀, NO₂, SO₂, CO, and O₃, alongside the AQI. PM_2.5, PM₁₀, NO₂, SO₂, CO, and O₃ are used as predictor variables, while AQI serves as the target variable. The data are sourced from official environmental monitoring platforms with high reliability and accuracy, ensuring the robustness of subsequent modeling tasks.

Situated in the core of the Greater Bay Area, Chang’an is surrounded by several industrial cities, and its air quality is influenced by both local emissions and regional pollutant transport, making it highly representative for studying industrial pollution dynamics. The results obtained here are expected to be generalizable to similar urban-industrial zones.

Figure 10 presents the CEEMDAN decomposition results for the AQI series. CEEMDAN is an improved version of the EMD method. By introducing adaptive noise and performing ensemble averaging, CEEMDAN enhances the stability of the decomposition and improves the separation of modes. The first row (red) shows the original signal, which contains multiscale information, including high-frequency noise and low-frequency trends. The middle green subplots represent different IMFs, arranged from high to low frequencies. The higher-order IMFs capture high-frequency noise and fast oscillations, while the lower-order IMFs gradually reveal low-frequency trends. The final row represents the residual trend component, reflecting the long-term variations of the signal. CEEMDAN provides better decomposition performance than traditional EMD and EEMD (Ensemble Empirical Mode Decomposition), effectively suppressing mode mixing and yielding IMFs with clearer physical significance, making it possible to decompose the target prediction data into interpretable features. With feature engineering completed, the next step involves building the forecasting model and conducting experimental analysis.

This study implemented the proposed CEEMDAN-GNN-Transformer and baseline models within a Python environment using the PyTorch deep learning framework. All experiments were conducted on a workstation configured with an NVIDIA GeForce RTX 3060 GPU, ensuring consistent training efficiency and reproducibility across experiments. Unless stated otherwise, all experiments were conducted under the same hardware and software conditions.

During the data preprocessing phase, the dataset was divided into training (80%) and test (20%) sets. The feature variables and prediction target were standardized using the StandardScaler method, yielding a mean of 0 and a standard deviation of 1 across both sets. This normalization mitigates the effects of differing variable scales, improves model training efficiency, and promotes stable convergence. The standardized data were then transformed into PyTorch tensors to enable high-performance batch processing within the deep learning pipelines.

To comprehensively assess the performance of the proposed model, this study adopts four widely used regression evaluation metrics: Mean Absolute Error (MAE), Mean Squared Error (MSE), Coefficient of Determination (R ²), and Explained Variance Score (EVS). These metrics assess model performance from multiple perspectives, including error size, model fit, and the extent of data variance explanation (Chai and Draxler, 2014; Chicco et al., 2021).

MAE is the average of the absolute differences between actual observations and predicted values.

MSE is the mean of the squared differences between actual and predicted values.

R ² quantifies the model’s ability to explain variance in the data, with values closer to 1 indicating better explanatory power.

EVS is similar to R ² but focuses more on the model’s ability to explain data variability, with higher EVS values indicating better fit.

To comprehensively assess the performance of the models in air quality forecasting, four common regression metrics were selected: Mean Absolute Error (MAE), Mean Squared Error (MSE), Coefficient of Determination (R ²), and Explained Variance Score (EVS). These metrics evaluate the model’s prediction performance from multiple dimensions, including error magnitude, goodness of fit, and its ability to capture the data’s variability. MAE measures the average absolute difference between the actual observations and the model’s predictions, while MSE captures the average squared error between actual and predicted values. The R ² coefficient quantifies the proportion of data variance explained by the model, ranging from 0 to 1, with values closer to 1 indicating stronger explanatory capacity. EVS operates similarly to R ² but focuses more explicitly on the proportion of data variance explained by the prediction, making it especially valuable for assessing the model’s ability to capture variability across the dataset.

To validate the effectiveness of the proposed CEEMDAN-GNN-Transformer model, this study selects nine representative predictive models for comparison. These models include deep learning approaches (e.g., CNN-Transformer), ensemble learning methods (e.g., XGBoost, GBDT, Random Forest), traditional machine learning methods (e.g., KNN, SVM), and classical time series models (e.g., ARIMA). These models exhibit significant differences in terms of time pattern modeling ability, variable relationship modeling strength, IMF component utilization, implementation complexity, and computational resource requirements, as summarized in Table 4 (Cui et al., 2023; Ye et al., 2025; Lei et al., 2023).

Through comparative analysis, it is evident that CNN-Transformer combines CNN’s local feature extraction ability and Transformer’s global dependency modeling capability, making it highly suitable for time series modeling. Compared to traditional time series models, CNN-Transformer efficiently extracts short-term patterns while modeling long-range dependencies, with relatively low computational complexity. By comparing it with GNN-Transformer, this study highlights the necessity of the GNN structure for temporal tasks and examines whether it enhances the Transformer’s performance in complex time series modeling. The comparison validates the advantage and applicability of GNN-Transformer in handling complex time series data and relational modeling tasks.

To explore the contribution of each key module in the proposed model and verify the combination effect, this study designs a systematic ablation experiment. Keeping hyperparameters consistent, the experiment sequentially removes or combines the signal decomposition module (CEEMDAN), the GNN, and the Transformer module, with a Multi-Layer Perceptron (MLP) used as the baseline model for comparison. Each experiment is repeated 20 times, with the average result taken to ensure robust and fair evaluation. The results are shown in Table 5.

The experimental results indicate that both GNN and Transformer significantly improve the model’s performance. Transformer, in particular, demonstrates a more substantial effect (MAE reduced from 16.97 to 14.032, R ² improved to 0.54). The combination of GNN and Transformer further enhances prediction capability (MAE reduced to 10.021, R ² improved to 0.79), indicating a synergistic effect. Moreover, the introduction of CEEMDAN signal decomposition significantly optimizes the Transformer’s prediction performance (MAE reduced to 8.9022, R ² improved to 0.81). When CEEMDAN is added on top of the GNN-Transformer combination, the performance improvement is limited (MAE decreases from 10.021 to 9.2495, R ² improves to 0.87). In summary, the CEEMDAN-GNN-Transformer architecture demonstrates superior modeling capability in time series prediction tasks, with the complementary strengths of the three components achieving the best performance in complex temporal modeling tasks.

This study employs the Optuna framework for hyperparameter optimization, aiming to minimize the MSE of the GNN and Transformer models on the validation set. Optuna is a Bayesian optimization-based automated hyperparameter search tool that intelligently explores the search space to find the optimal combination of hyperparameters with the fewest trials. It supports various machine learning and deep learning frameworks, including Scikit-learn, PyTorch, and TensorFlow. By defining the search space, the optimizer selects appropriate hyperparameters such as the number of attention heads, hidden layer dimensions, model depth, learning rate, dropout rate, and number of neighbors. Through 50 trials, Optuna automatically tunes the hyperparameters to find the best combination, improving model prediction accuracy, reducing computational cost, and enhancing generalization ability.

Table 6 lists various hyperparameters affecting model performance, along with their respective search ranges. These include structural parameters (scaling factor for hidden layers, attention heads, and encoder layers), graph structure parameters (number of neighboring nodes in the graph), regularization parameters (dropout rate), and optimizer parameters (learning rate). The tuning of these hyperparameters directly impacts model training performance and generalization ability. After 20 iterations, the optimal combination of hyperparameters was obtained. A partial example of the optimization results is shown in Table 7.

As shown in Table 7, the optimal hyperparameter combination found through the optimization process includes 8 attention heads, a hidden dimension scaling factor of 47, 1 attention layer, a learning rate of 0.000153, a dropout rate of 0.197, and 5 neighbors. This combination resulted in the best performance, with a loss of 0.1929, demonstrating the model’s best performance under these settings. The actual hidden layer dimension is calculated by multiplying the hidden_dim_factor by num_heads.

The final model was built using the best parameters obtained from Optuna’s optimization. The following steps were undertaken: the graph structure was constructed, with adjacency matrices calculated based on feature similarities for both the training and testing sets, where the k nearest neighbors were identified for each data point. The edge index for GNN computation was then built. The same k_neighbors value was used for both training and testing sets, but separate graph structures were created for each. The model was trained for 200 epochs. Each epoch included a forward pass through the GNN layers, Transformer layers, and a fully connected layer, followed by calculating the mean squared error loss between predicted and true values. Backpropagation was performed to compute gradients, and the Adam optimizer was used to update the model parameters.

From the training logs, as shown in Table 8, the loss decreases progressively over the training process, indicating that the model is converging and improving its prediction accuracy. The loss decreases rapidly during the initial stages (Epochs 20–40), from 0.3592 to 0.2892, reflecting the larger optimization steps in the early training phase. In subsequent stages, the loss stabilizes, indicating that the training process is becoming more stable. By Epoch 200, the loss converges to 0.1635. Overall, the decreasing trend in loss indicates the model’s effective training and suggests that it achieves an optimal performance level.

Table 9 provides a summary of the final model’s parameter counts, highlighting its complexity and capacity.

The model developed in this study integrates a GNN and a Transformer architecture to process graph-structured data and capture global dependencies. The input feature dimension is 17, and after two layers of GNN processing, it is mapped to a 376-dimensional hidden representation. Subsequently, a single-layer Transformer encoder with 8 attention heads is used to further model global feature interactions. Finally, a fully connected layer maps the features to a 1-dimensional output for pollutant prediction. The total number of parameters is approximately 1.85 million, making it a relatively complex model for pollutant time series forecasting.

To comprehensively evaluate the forecasting performance of the proposed CEEMDAN–GNN–Transformer model for air quality prediction, a comparative analysis was conducted against a range of classical machine learning and deep learning approaches. All models were tested on the same dataset preprocessed via CEEMDAN signal decomposition, and their performance was assessed using four commonly used metrics: MAE, MSE, R ², and EVS. The results are presented in Table 10.

As shown in Table 10, the forecasting performances of the compared models vary significantly. The GNN–Transformer approach achieves the best results across all metrics, with the lowest MAE (7.6495) and MSE (98.2009), along with the highest R ² (0.9192) and EVS (0.9491). This highlights its superior ability to capture complex temporal dynamics and to accurately explain the variance within the data. Meanwhile, XGBoost also delivers strong performance, yielding an R ² of 0.8890 and an EVS of 0.9044, underscoring its effectiveness in capturing feature interactions and nonlinear relationships. In comparison, Random Forest, SVM, and CNN-Transformer perform slightly worse than XGBoost and GNN-Transformer in terms of MSE and MAE, while Adaboost and KNN show relatively weaker performance in all metrics. Notably, KNN has the highest MSE (422.7499) and the lowest R ² (0.5454), which indicates its limitations in modeling complex dynamic relationships, especially in high-dimensional regression tasks. Overall, the GNN-Transformer and XGBoost models demonstrate superior prediction performance in this study, while traditional machine learning methods (e.g., KNN, Adaboost) show certain limitations in modeling complex dynamic changes.

Figure 11 illustrates the forecasting results of the GNN-Transformer model compared with observed AQI values. The predicted time series closely tracks the actual observations, demonstrating the model’s strong temporal modeling capability and its ability to capture long-term trends in air quality fluctuations. However, slight discrepancies are observed during certain extreme events, such as pollution peaks or sharp declines, where the predicted values tend to be smoother than the actual data. This limitation may arise from the signal decomposition and denoising process of CEEMDAN, which can attenuate high-frequency variations. Despite this, the proposed GNN-Transformer approach shows considerable promise in capturing the overall dynamics of air quality. In future work, incorporating meteorological covariates or further refining the model architecture could help better capture these abrupt variations and further improve prediction precision.

Although CEEMDAN is effective in reducing noise and extracting multi-scale features, its decomposition inherently introduces a smoothing bias, diminishing the retention of high-frequency signals. In this study, extreme pollution episodes (such as sharp AQI surges or declines) appear attenuated in the fitted curve. While beneficial for overall stability and trend extraction, this effect may be suboptimal for real-time air quality management, where timely detection of extreme events is critical for public health and emergency response. Future work could consider augmenting CEEMDAN with high-frequency component weighting or integrating auxiliary submodels optimized for peak detection to enhance responsiveness to short-term pollution shocks.

In environmental air quality forecasting, beyond accuracy, it is equally important to quantify the uncertainty of predictions, as this information supports scientific decision-making and risk management. To this end, the Monte Carlo (MC) Dropout method—a Bayesian approximation technique—was applied to the CEEMDAN–GNN–Transformer model for uncertainty analysis.

Originally proposed by Gal and Ghahramani. (2016), MC Dropout involves keeping the dropout mechanism active during inference, performing multiple stochastic forward passes to simulate sampling from the posterior weight distribution. This enables approximate estimation of the predictive distribution. For each of T stochastic passes, predictions y t ^ are collected, from which the predictive mean and variance are computed, as shown in Equations 15, 16: E y * ^ ≈ 1 T ∑ t = 1 T y t ^ (15) Var y * ^ ≈ 1 T ∑ t = 1 T y t ^ − E y * ^ 2 (16)

Using these results, two uncertainty evaluation metrics were calculated: Prediction Interval Coverage Probability (PICP) and Mean Prediction Interval Width (MPIW). The results are shown in Table 11.

The results show that the PICP (0.916) is close to the ideal 95% coverage, indicating that the prediction intervals successfully cover the majority of observed values. The MPIW (22.001) is substantially narrower than the typical AQI category width (≈50 for mild pollution), implying high-precision prediction intervals.