Five typical models, namely, the multilayer perceptron (MLP) (You et al., 2017), ANFIS (Jang, 1993), LSTM (Muzaffar and Afshari, 2019), SVM (Brereton and Lloyd, 2010), and GRNN (Land and Schaffer, 2020), have been widely used for air pollutants forecasting because of their robustness, efficiency, and accuracy.

For the DE algorithm and gray wolf optimization (GWO) algorithm, the defects of prematurity, poor stability, and ease in falling into local optimum will occur when solving the optimization problem separately. Combining the advantages and disadvantages of the two algorithms, a more efficient hybrid optimization algorithm, DEGWO algorithm, is proposed to improve global search capabilities.

Firstly, in order to avoid the phenomenon in which the population is iteratively reduced to a certain area, the crossover and selection operations of the DE algorithm are used to maintain the diversity of the population. Then, as the initial population of the GWO algorithm, the objective function value of the individual is calculated. The optimal three individuals X _α , X _β , and X _δ are selected to update the positions of other gray wolf individuals. Then, the position of the gray wolf individual is updated by the intersection and selection operations of DE, and the iterative update is repeated until the optimal one is selected. The target function value is output.

The hybrid algorithm not only improves the global search ability but also effectively avoids the defects of early maturity stagnation and falling into local optimum. The specific implementation steps of the algorithm are as follows:

Step 1: Set the relevant parameters of the hybrid optimization algorithm, population size N, maximum iteration number t _max , crossover probability CR, search dimension D, search range ub, lb, and scaling factor range F.

Then the competition selection operation is performed according to Eq. 2 to generate. x i ( g + 1 ) = { v i ( g ) ,   f [ v i ( g + 1 ) ]   < f [ x i ( g + 1 ) ]   x i ( g ) ,   f [ v i ( g + 1 ) ]   ≥   f [ x i ( g + 1 ) ] (2)

Step 3: Calculate the objective function value of each gray wolf individual in the population, sort according to the size of the objective function value, and select the optimal first three individuals as X _α , X _β , and X _δ , respectively.

Finally, update the current position of each gray wolf individual according to Eqs 6–9. X 1 ( t + 1 ) = X α ( t ) − A 1 D α (6) X 2 ( t + 1 ) = X β ( t ) − A 2 D β (7) X 3 ( t + 1 ) = X δ ( t ) − A 3 D δ (8) X p ( t + 1 ) = X 1 + X 2 + X 3 3 (9)

Step 5: Update the values of a, A, and C in the algorithm, cross-operate the position of the individual population according to Eq. 1, retain the better components, then perform Eq. 2 to select new individuals and calculate the objective function values of all gray wolf individuals.

Fuzzy C-means clustering (FCM), known as fuzzy ISODATA, is a clustering algorithm that uses membership degrees to determine the extent to which each data point belongs to a certain cluster. In 1973, Bezdek proposed the algorithm as an improvement to the early hard C-means clustering (HCM) method (Gayen and Biswas, 2021). The clustering steps are as follows:

Step 1: Initialize the membership matrix U with a random number whose value is between 0 and 1, so that it satisfies the constraint in Eq. 11.

The process of establishing a fuzzy synthetic evaluation (FSE) system is as follows (Lu et al., 2011).

Step 1. The set of factors for the evaluation object is determined.

The selected factors should possess the traits of representativeness, feasibility, and system. Air quality evaluation relies on the concentration levels of the main air pollutants. Therefore, in this study, the indicators were chosen according to China's ambient air quality standards (AAQS: GB3095-2012). Moreover, different geographical areas have different topographic and economic characteristics, and consequently, the different key pollutants in the study areas should be also considered.

Step 2. The evaluation rank standard is determined.

The evaluation rank set is described as V = {v ₁, v ₂, ⋯ ,v _n }. In our study, the air pollution degrees were divided into five levels. The pollutants grading standard according to AAQS is shown in Supplementary Appendix S1.

Step 3. Index fuzzification.

In this step, the membership functions (MFs) corresponding to each index are obtained. The process of fuzzification constitutes the process of membership calculation by using MFs. In this study, we used the trapezoidal membership to calculate the membership value.

Step 4. The factor weight is calculated.

Weight reflects the importance of each factor in synthetic evaluation and directly affects the outcome of the evaluation. Many methods exist for determining the weight, such as weighted statistics, coefficient of variation method, the Delphi method, and entropy methods. In our study, the weight was calculated by fuzzy weighting method.

Step 5. The evaluation results are output.

The objective function of the DEGWO algorithm is based on stability and accuracy, in which MSE is the standard to measure accuracy and the variance of error is the standard to measure stability. Algorithm 1 briefly outlines the process of the MODEGWO.

Algorithm 1

MMODA

Input: Objective function Min fitness(x) = f ₁ + f ₂ m i n { f 1 ( x ) = M S E = 1 N ∑ i = 1 N ( x ^ i − x i ) ,   i = 1,2 , ... , N f 2 ( E ) = V a r ( E ) = 1 N ∑ i = 1 N ( E i − E ¯ ) 2 ,   i = 1,2 , ... , N

Note: E i is the test error; the calculation equation is E i = x ^ i − x i

  x ^ i and x i are the actual data and output data by each model

  Parameters of DEGWO

  CR is crossover probability: 0.2;

  MaxGen is the maximum number of the iteration: 500;

  F is the scaling factor: [0.2, 0.8];

  p_size is population size: 50.

Output: The optimal solution and the best objective function value.

Initialize a parent population, mutant population, and child population of gray wolf with a random position in a feasible region using equation; x p k = x p k ( l o w ) + ( x p k ( u p ) − x p k ( l o w ) ) × r a n d ( 0,1 )

Note: x p k ( l o w )   is the lower bound of the pth component of the kth individual.

  x p k ( u p ) is the upper bound of the pth component of the kth individual.

 rand(0, 1) represents a random number in [0, 1].

 p = 1, 2 ⋯ , d. k = 1, 2 ⋯ , p_size

Initialize crossover probability Pc and scaling factor F;

Initialize a, A, and C;

Evaluate f for all individuals in the parent population;

Sort the parent population in a non-decreasing order, according to the objective function value;

X _α is the best individual in the parent population of gray wolves;

X _β is the second individual in the parent population of gray wolves;

X _δ is the third individual in the parent population of gray wolves;

While (t < MaxGen)

for each individual in the parent population of gray wolves

Update the position using the following equation; X p ( t + 1 ) = X 1 + X 2 + X 3 3

end for

Obtain a mutant population of gray wolves using the following equation; V i ( g ) = X r 1 + F ⋅ ( X r 2 ( g ) + X r 3 ( g ) )

Note: g is the generation number,

 F is the scaling factor, g = 0, 1, 2, ⋯ , MaxGen,

 MaxGen is the maximum number of the iteration

Obtain a child population of gray wolves using the following equation; v i j ( g + 1 ) = { h i j ( g ) , r a n d ( 0,1 ) ≤ C R j = r a n d ( 1 , n ) x i j ( g ) , r a n d ( 0,1 ) > C R j = r a n d ( 1 , n )

CR represents the crossover probability

for each individual Parent _i in a parent population of gray wolves

If f(Child _i ) < f(Parent _i )

Replace Parent _i with Child _i

end if

Update A, C, and a;

Sort the parent population of gray wolves in a nondecreasing order;

Update X _α , X _β , X _δ ;

t = t + 1;

End while

Return Parentα and f(Parentα).

The hybrid AQI forecasting system in this paper is composed of the above three parts. A flow chart of the hybrid model is presented in Figure 1.

From the above, we can see that the AQI forecasting step using the hybrid forecasting system proceeds as follows:

Step 1: Data Mining

The AQI and six air pollution time series with missing points is filling processed by shape-preserving piecewise cubic spline interpolation, which maintains the continuity of each time series.

Step 2: Feature selection and data setting for each model

To ensure the forecasting performance, a modified optimization algorithm is used to further optimize the parameters of the best forecasting model (expect LSTM). Finally, the main air pollutants forecasting results are obtained and compared with those of different hybrid forecasting models.

The BJ-TJ-HE is the national capital region of the People's Republic of China. It is the biggest urbanized megalopolis region in Northern China, where Beijing, Tianjin, Baoding, and Langfang are the central core areas of BJ-TJ-HE. In this paper, the 13 cities of BJ-TJ-HE are evaluated to develop an early warning indicator for air quality. The datasets of hourly concentrations of the six major air pollutants used in this study are retrieved from the website of the China National Environment Monitoring Centre (http://www.cnemc.cn/sssj/). The first dataset includes AQI hourly concentrations collected from January 1, 2017, to December 31, 2018, in BJ-TJ-HE. Figure 2 shows the result of fuzzy C-mean cluster, which displays the construction of a fuzzy matrix based on the attributes of AQI in 13 cities and objectively and accurately cluster (Category I: Beijing, Baoding, Langfang; Category II: Shijiazhuang, Tangshan, Handan, Chengde, Hengshui, Xingtai; Category III: Tianjin, Qinhuangdao, Zhangjiakou, Cangzhou). The result of fuzzy comprehensive evaluation is shown in Table 1, which found that the main air pollutants are PM₁₀, PM₂, and NO₂ in 13 cities.

According to the analysis in Table 1, the main air pollutants from statistical analysis of BJ-TJ-HE are NO₂, PM_2.5, and PM₁₀ shown in Table 2, in which the average value of the main air pollutants shows obvious differences among the 13 cities. The average value of NO₂ in the different cities is between 22.2525 and 49.4348 μg/m³, in which the average value in Xingtai is higher than in the other cities. At the same time, the PM_2.5 and PM₁₀ average values in Xingtai are 69.6938 and 135.8368 μg/m³, which are also higher than in the other cities. The maximum values of NO₂, PM_2.5, and PM₁₀ were in Hengshui, Baoding, and Zhangjiakou, with values of 215, 402, and 1581 μg/m³, and the minimum values of the three main air pollutants were in Zhangjiakou, Beijing, and Zhangjiakou, with values of 1, 3, and 12 μg/m³.

In terms of skewness, all data sets are rightward, with values of skewness are greater than 0. For the values of kurtosis, only three data sets of NO₂ were less than 3, which meant that these three sets (Qinhuangdao, Shijiazhuang, and Xingtai) had a fat tail. At the same time, the other data sets had a thin tail.

This section focuses on the efficiency of the different forecasting model with respect to computational performance. Eight evaluation criteria are applied to estimate the forecasting performance, namely, mean absolute error (MAE), root mean square error (RMSE), mean absolute percentage error (MAPE), Theil U statistic 1 (U1), and Theil U statistic 2 (U2) were calculated for all the fits; the goodness of forecasting fit (R ²) and the standard of forecasting error (STDE) indicates the stability of the forecasting models; and the direction accuracy (DA) evaluates the optimal decision-making, often relying on correct forecasting directions or turning points between the actual and forecasting values. These performance metrics are defined in Table 3.

In the forecasting processing, there is no model that can be applied to all time series in the process of forecasting. Therefore, in this paper we developed a new metric, which measures accuracy of each hybrid model testing set and determines whether the model can provide the optimal forecasting value. The process of model selection is as follows:

Each model data is divided into 840 training samples, 168 testing samples, and one forecasting value. The accuracy of the testing sample is calculated by using the WIC. In order to eliminate the difference of the order of magnitude of forecasting metric, the MAE, MAE RMSE, MAPE, STDE, U1, and U2 are normalized. The calculation formula is as follows: WIC=NMAE+NRMSE+NMAPE+NSTDE+U1+U2+R2+ ( 1-DA )

For the first forecasting, the 1st to 840th samples are the training samples, the 841st to 1008th samples are the testing samples, and the 1009th sample is the forecasting value. At the end of the forecasting, the WIC value of the testing sample is calculated. If the WIC of the ith model is the smallest, the forecasting value of the ith model provides the optimal forecasting value.

For the second forecasting, the 2nd to 841st samples are the training samples, the 842nd to 1009th samples are the testing samples, and the 1010th sample is the forecasting value. At the end of the forecasting, the WIC value of the testing sample is calculated. If the WIC of the ith model is the smallest, the forecasting value of the ith model provides the optimal forecasting value.

For the kth forecasting, the kth to (840 + k − 1)th samples are the training samples, the (840 + k)th to (1008 + k − 1)th samples are the testing samples, and the (1008 + k)th sample is the forecasting value. At the end of the forecasting, the WIC value of the testing sample is calculated. If the WIC of the ith model is the smallest, the forecasting value of the ith model provides the optimal forecasting value.

In the whole forecasting process, 168 optimal forecasting values are generated, and 168*5 WIC values are generated at the same time. Table 4 only shows the optimal model and the percentage of optimal forecasting value for three main air pollutants.

From Table 4, it can be seen that SVM provides more optimal forecasting value for the three main pollutants at different times, especially in the PM₁₀ forecasting process; the optimal forecasting value for the first quarter and the third quarter is 82.14% (138 optimal forecasting value), and the other four models also provide corresponding optimal forecasting value.

In this portion, the hourly NO₂ time series for 13 cities in three categories were utilized as the testing data for the five hybrid models with one-step-ahead forecasting. Beyond that, with the purpose of comprehensively comparing the precision of the modeling forecasting, this experiment consisted of two parts: the multi-step forecasts demonstrated in Table 4 and, for the local analysis horizon, the local forecasts presented in Table ′5 and Figure 3, which focus on first season. Table 5 and Figure 3 demonstrate the following:

1) Focusing on Category I, the new proposed model based on model selection realizes excellent results on the eight evaluation indices in the first season forecasting. On the contrary, DEGWO-ANFIS has the lowest effectiveness. The maximum reduction of MAPE for the proposed model compared with the other hybrid models is approximately 71.18% in Beijing's NO₂ forecasting, 53.93% in Baoding's NO₂ forecasting, and 61.61% in Langfang's NO₂ forecasting, respectively. The reduction was about 62.39% and 76.49% for one-step forecasting and 2.79%, 6.10%, and 19.33% for the three cities at the hourly interval NO₂ forecasting in Category I. As for the R ², the proposed model has the best performance among the four single-hybrid models for hourly interval NO₂ time series.

In summary, from the analyses above, it can be concluded that the model selection forecasting system realizes the best forecasting results compared to the single hybrid model. Model selection also gives better forecasting performance in the other season with the results shown in Supplementary Appendices S2–S4 indicating better robustness of the model selection forecasting system.

This experiment mainly focused on the forecasting performance of each model for PM_2.5 of Category II in the first season, with the forecasting results of four different hybrid models (MODEGWO-SVM, MODEGWO-BPNN, MODEGWO-ANFIS, Adam-LSTM) and model selection represented in Table 6 and Figure 4

1) For first season PM_2.5 forecasting accuracy, the final forecast results of PM_2.5 for six cities in Category II are composed of four hybrid models, which include MODEGWO-SVM, MODEGWO-BPNN, MODEGWO-ANFIS, and Adam-LSTM. Among the four models, MODEGWO-SVM and MODEGWO-BPNN have better forecasting performance, with the MODEGWO-SVM obtaining 64.29% optimal forecasting points and the MODEGWO-BPNN obtaining 21.43% optimal points for six cities in Category II. The smallest MAPE values of MODEGWO-SVM are 0.92%, 0.94%, 1.36%, and 0.79% for Hengshui, Tangshan, Chengde, and Xingtai PM_2.5 forecasting, and the MODEGWO-BPNN obtains the best MAPE (1.08% and 0.85%) value for Shijiazhuang and Handan.

In summary, for the Category I and Category III PM_2.5 forecasting list in Supplementary Appendices S5–S7, the model selection forecasting system exhibits the best forecasting accuracy among the different hybrid models for four seasons. It is evident that the forecasting capacity of the model selection is robust when considering each forecasting metrics. The accuracy of model selection depends on the hybrid model, so it is necessary to increase the types of models in the modeling process which ensures that more forecasting results can be obtained, and the optimal forecasting value can be selected in the model selection process.

For the hourly PM₁₀ time series for three categories, it can be observed that the model selection forecasting system attains satisfactory results. Specifically, the lowest values of MAE are 0.4643, 0.4600, and 0.3869 and of RMSE are 0.7302, 0.7906, and 0.5561, corresponding to PM₁₀ forecasting in Category I in three cities, successively. The results of Category I indicate that the smaller the MAE and MSE, the smaller the deviation between the observations and forecasting, which verifies the forecasting performance. In addition, the average MAPE values of model selection for six cities of Category II in four seasons are lower than 1%. Compared with the optimal hybrid model, the model selection is approximately reduced by 10%. The analyses reveal the forecasting superiority of the model selection system. By parity of reasoning, a similar conclusion can be reached through the analyses of the hourly PM₁₀ forecasting results for Category III (the forecasting results are shown in Table 7 and Figure 5).

1) For PM₁₀ forecasting in the first season, the optimal hybrid models are DEGWO-SVM and DEGWO-BPNN, with which the MAPE values of the best hybrid model (DEGWO-SVM) for four cities of Category III are 0.71%, 0.81%, 1.09%, and 0.72%. The final forecasting values are obtained by model selection; based on the results of DEGWO-SVM and DEGWO-BPNN the MAPE values of model selection are 0.65%, 0.75%, 0.97%, and 0.67%. Additionally, the values of other forecasting metrics are at their best under the model selection.

In summary, whether for Category III or the other categories (the results are shown in Supplementary Appendix S8 and Supplementary Appendix S9) PM₁₀ forecasting, the model selection system attained the best performance for 13 cities. In the comparison of various hybrid models, the forecasting performance of MODEGWO-SVM is better than other hybrid models. Additionally, it can be observed that some models have low accuracy, which can still provide some optimal forecasting values for the model selection system for hourly PM₁₀ forecasting. Based on the above analysis, it can be seen that none of the models has been playing the best forecasting performance in the forecasting process, and various hybrid models are needed to make up for the shortcomings of the single hybrid model.

According to the results of each experiment, SVM provided the more optimal forecast values for the three main pollutants in the four quarters of 13 cities. The reason for the favorable score produced by SVM is that SVM provides a way to avoid the complexity of high-dimensional space by directly using the inner product function of the space (which is the kernel function) and then directly solving the corresponding decision-making problem in high-dimensional space by using the solution method under the condition of linear separability. When the kernel function is known, it can simplify the difficulty of solving the problem in high-dimensional space. Meanwhile, SVM is based on the small sample statistical theory, which conforms to machine learning. In addition, support vector machine has better generalization ability than neural network. Although the time series of the three main air pollutants are neither regular nor seasonal, SVM can also effectively capture future changes of the three main air pollutants.

SVM has two very important parameters: c and g. c is the penalty coefficient, that is, tolerance of errors. The higher the c, the more intolerable the errors and easy to over-fit. The smaller the c, the less easy fitting is. If c is too large or too small, the generalization ability becomes worse. g is a parameter that comes with RBF function when it is selected as a kernel. Implicitly it determines the distribution of data after mapping to a new feature space. The bigger the g, the less support vectors; it will only act near the support vector samples. For the unknown samples, the classification effect is very poor. There is a possibility that the training accuracy can be very high, and the test accuracy is not high, that is, over-fitting. The smaller the g value is, the more support vectors there are, and the greater the smoothing effect will be; the higher accuracy of the training processing cannot be obtained, and the accuracy of the testing processing will also be affected. This paper used DEGWO algorithm and GWO algorithm to optimize the parameters of SVM (c and g). The results of two types of hybrid SVMs are shown in Table 8, which displays that the optimum penalty coefficients of SVM corresponding to pollutant forecasting in different cities vary widely. For example, in the forecasting processing of NO₂, the variation range of parameters is [2, 99]. However, the fluctuation range of g is small, with most variations ranging from 0 to 1. In general, the air pollutants forecasting performances of support vector machine are very dependent on the penalty coefficient. In the whole experiment it can be observed that the support vector machine has good forecasting accuracy for three main air pollutants forecasting, but it cannot provide the best forecasting value in each point. It indicates that the support vector machine is suitable for hourly air pollutants forecasting.

ANN as a nonlinear mapping model is used to solve the problem of time series forecasting, because the ANN model can find the optimal solution of a complex problem with the help of high-speed computing ability of the computer. In order to ensure the forecasting accuracy of the ANN model, parameters of ANN need to be elaborately configured. However, there is no effective rule for establishing the values of these parameters on air pollutants forecasting. Although there are many studies on the tuning of the parameters of the neural network, it is obvious that the selection of the whole parameter space is beyond the scope of this study. Therefore, the parameters of the neural network are set by means of simulation experiment and optimization algorithm: the experimental design is as follows:

This processing configured various input layers and a number of hidden layers to find out the influence of the usage of recent data on the performance of different ANN models. The number of input layers from 1 to 10 increases for three main air pollutants, which means there are 1,008 pieces of sample data on NO₂, PM_2.5, and PM₁₀; the train-to-verify ratio 5:1 means that 840 pieces of sample data were used as training data for building the ANN model, while 168 pieces of sample data were used as testing data for finding the training-to-testing ratio and parameter of each ANN model (the optimal number of input layers of each model and the number of hidden layers of LSTM and BPNN). Figure 6 shows the forecasting performance with the different configurations of the optimal number of input layers and the number of hidden layers of LSTM and BPNN, in which it is difficult to find a regular correlation between the forecasting performances and the parameters. Consequently, it is difficult to find an optimal combination of ANN's parameters that brings the model to the best performance in the practical air pollutant forecasting where MAPE and R ² are unknown.

During the experiment on configurations of ANN's parameter, the optimal parameter setting trained the networks models for each ANN model in 13 cities. The forecasting values with the best performance of each ANN (the best forecasting metric in Table 9) were selected to forecast three main air pollutants in each experiment. However, it can be found that there were no giant differences on forecasting performance among the networks trained with the same configuration, even if the neural network has the randomness and probability mechanism inside the training processing. A large sample of the times series is another reason that the training stability of the neural network can be ensured. For example, the best MAPE of BPNN for the forecasting of NO₂ in Beijing shown in Figure 6 is 5.82%, and the worst one is 8.94%. Most of the MAPEs are between 6% and 7%. It is practical to use ANN in real air pollutants forecasting application where forecasting the changing air pollutant time series is suitable.

In summary, with the rapid development of ANN, it has become a powerful tool to solve prediction problems. Neural network is used in the field of air pollution to solve the problem of non-linear forecasting which cannot be solved by statistical models. Its non-linear mapping is especially suitable for the application of air pollutant forecasting. The main reason is that the ANN has the following advantages:

1) Non-linear mapping ability: ANN realizes a mapping function from input to output in essence. Mathematical theory proves that three-layer neural network can approximate any non-linear continuous function with arbitrary precision. This makes it especially suitable for solving complex internal mechanism problems, that is, ANN has strong non-linear mapping ability.

In order to improve the computing efficiency and save the computing time, training and forecasting processing of all the models for the main air pollutants time series with parallel computing by central processing unit (CPU) and graphics processing unit (GPU). The computing times of every independent hybrid model in each experiment are shown in Table 10, from which we can further research the computational efficiency of the developed model selection forecasting system for the main air pollutants. Specifically, the average computation time of the model selection forecasting system ranges from 330.6037 to 363.2167 s for NO₂ forecasting in Category II, with the longest computing time appearing in the different categories. Notably, Adam-LSTM with complex model structure has longer computing time than the other hybrid models, taking more time in the iterative optimization process. Moreover, this paper establishes multiple hybrid models and uses the model selection method to find the best forecasting value, in which the final forecasting accuracy is improved but needs more computing time.