Prediction of Air Pollution Concentration Based on mRMR and Echo State Network

Abstract

Air pollution has become a global environmental problem, because it has a great adverse impact on human health and the climate. One way to explore this problem is to monitor and predict air quality index in an economical way. Accurate monitoring and prediction of air quality index (AQI), e.g., PM_2.5 concentration, is a challenging task. In order to accurately predict the PM_2.5 time series, we propose a supplementary leaky integrator echo state network (SLI-ESN) in this paper. It adds the historical state term of the historical moment to the calculation of leaky integrator reservoir, which improves the influence of historical evolution state on the current state. Considering the redundancy and correlation between multivariable time series, minimum redundancy maximum relevance (mRMR) feature selection method is introduced to reduce redundant and irrelevant information, and increase computation speed. A variety of evaluation indicators are used to assess the overall performance of the proposed method. The effectiveness of the proposed model is verified by the experiment of Beijing PM_2.5 time series prediction. The comparison of learning time also shows the efficiency of the algorithm.

1. Introduction

With the rapid advancement of urbanization and industrialization, air quality has deteriorated severely, which has negatively affected the quality of the living environment and even hindered economic growth in some areas [1]. In particular, the inhalable particles produced by industrial pollution have small particle size, large diffusion area and strong activity, thus they can enter the human body through the respiratory tract, which has adverse effects on human health. Therefore, the prediction of air pollutants plays a crucial role in the early warning and control of environmental pollution [2]. Therefore, modeling and forecasting the air quality index (e.g., PM_2.5 concentration) has become an effective way to prevent and control air pollution, and it also provides a scientific basis for the development of effective measures [3]. The implementation of this idea can effectively reduce the health hazard of air pollution, thus achieving early warning and rational planning [4].

For a long time, many scholars have conducted in-depth research on air pollution. At the same time, a variety of predictive models have been proposed, such as autoregressive integrated moving average model [5], support vector machine [5], multiple linear regression model [6], neural networks [7,8], and so on [9]. All of them have been applied to predict air pollution concentration. Oprea et al. [10] applied artificial neural network and adaptive neuro-fuzzy inference system to predict PM_2.5 concentration. Deng et al. [11] proposed heterogeneous space-time artificial neural networks to deal with spatial heterogeneity, which has been applied to predict the concentration of fine particles in Beijing–Tianjin–Hebei. Ong et al. [12] proposed a deep recurrent neural network for time series prediction, which can accurately predict the PM_2.5 concentration. In order to solve the air pollution monitoring and prevention in Kunming and Yuxi, China, Li et al. [13] studied the method called co-integration flower-pollination-algorithm support vector machine for the prediction of PM₁₀ concentration time series. In order to study the attenuation effect of haze on solar radiation scattering, Yao et al. [14] proposed a new SVM-based method that increases the accuracy of global solar radiation models by increasing variables of daily global solar radiation, sunshine hours, temperature, relative humidity, and air quality index. Reid et al. [15] selected an optimal prediction model based on 10-fold cross-validation to estimate PM_2.5 concentrations during wildfires in Northern California in 2008, and reliably predicted major wildfire events. As can be seen, different types of predictive models can successfully solve different problems of air pollution time series prediction.

Because the causes of air pollution are very complicated, analyzing the main pollutants and influence variables of air quality index will lay the foundation for the establishment of predictive models. Although the methods described earlier are effective, they do not analyze the validity of the input variables, so the models are likely to contain irrelevant or redundant information. In order to build an accurate predictive model, Sun et al. [16] applied Pearson correlation coefficient to analyze the relationship between PM_2.5 and other variables and selected the appropriate input variables according to the correlation order. Zhang et al. [17] investigated cross-correlations between PM_2.5 and four meteorological factors based on multifractal detrended cross-correlation analysis method, which reveals the impact of meteorological variables on PM_2.5 concentration. Zhu et al. [18] proposed a novel graphical causality analysis approach and analyzed the impact of meteorological and traffic variables on air quality indexes. Chen et al. [19] applied the convergent cross-mapping method to analyze the causal relationship between meteorological factors and PM_2.5 concentration in Beijing–Tianjin–Hebei region, and obtained quantitative causality analysis results. Therefore, for complex air pollution problems, we also need to choose the appropriate correlation or causality analysis method to analyze the impact variables, thus achieving high-precision detection and forecasting.

At present, air pollution has become one of the major environmental problems. However, the causes of air pollution are very complicated. Taking the PM_2.5 concentration as an example, its concentration is not only affected by the air pollutants such as NO₂, CO, O₃, and SO₂, but also by meteorological variables such as temperature, pressure, humidity, wind speed and wind direction. How to choose effective information from a variety of variables for prediction is an important research topic. However, in the previous studies, the correlation between influencing factors and PM_2.5 was not considered in the establishment of most predictive models. In order to solve the above problem, this paper considers using the minimum redundancy maximum relevance (mRMR) [20] feature selection method to select the appropriate input variables, which can select the most relevant information and reduce redundant information. Furthermore, considering the chaotic characteristics of variables, we use phase space reconstruction to extract evolutionary information of relevant variables. At last, the new input variables are transferred to the supplemental leaky integrator echo state network (SLI-ESN) for prediction. The improved model not only enhances the feature extraction and memory ability of the reservoir for multivariate time series, but also improves the influence of historical evolution state on the current state. In practical applications, accurate prediction is conducive to monitoring air quality and making reasonable and scientific decision-making on air pollution prevention. In order to verify the validity of the proposed method, we select the dataset of air quality index and meteorological time series to predict PM_2.5 concentration in Beijing, China.

The rest of this paper are arranged as follows: Section 2 describes the preliminary knowledge of this paper, including feature selection method and echo state network; Section 3 introduces in detail the four aspects of the feature selection method, phase space reconstruction, PM_2.5 prediction model and algorithm flow; In Section 4, we analyzed the experimental results of the time series prediction of PM_2.5 concentration in Beijing, China. Section 5 gives conclusions and illustrates the challenges that will still be faced in the future.

2. Preliminaries

2.1. Feature Selection Method

Time series that affect air pollution are high-dimensional data, which not only contains rich information, but also has irrelevant or redundant factors. These irrelevant and redundant factors reduce the prediction accuracy and efficiency of the model. Thus, analyzing the relationship between variables and selecting valuable input variables are important for prediction.

Feature selection is the most typical data preprocessing method [21]. It consists of four parts: generation process, evaluation function, stop criterion and verification process. Currently, Common feature selection algorithms include random forest (RF), correlation feature selection (CFS), fast correlation-based filter (FCBF), mutual information (MI) [22], information gain (IG) [23], regularization models, relief-based algorithms, and genetic algorithm. However, these feature selection algorithms usually ignore the redundancy relationship between features. Random forests have over-fitting on noise regression problems and do not give continuous output. The CFS and FCBF feature selection methods are slow to calculate and cannot handle large-scale data efficiently. MI and IG do not require the type of data distribution but have high computational complexity for high-dimensional data. The regularization model has a better selection effect on high-dimensional data with a size much larger than the number of samples, but it does not perform satisfactorily in low-dimensional data selection. For the prediction of PM_2.5 concentration with high noise, how to select appropriate feature selection method is very important. In this paper, the minimum redundancy maximum relevance (mRMR) [20] is used for feature selection. After this, a set of purest features can be obtained, and the redundancy features are removed while guaranteeing the maximum relevance.

2.2. Echo State Network

The echo state network (ESN) proposed by Jaeger et al. [24] is a new type of recurrent neural network. It consists of an input layer, a reservoir and an output layer. The structure of ESN is shown in Figure 1. ${\mathit{W}}_{in}$ is the input weight matrix that connects the inputs and the reservoir, $\mathit{W}$ is the internal connection weight matrix of the reservoir, ${\mathit{W}}_{out}$ is the output weight matrix that connects the reservoir and the output, and ${\mathit{W}}_{back}$ is the output-to-reservoir feedback connection weight matrix, which is usually set as a zero vector.

Assuming that the system obtains a time series $\mathit{U}(n)={\left[u_1(n),u_2(n),\cdots ,u_K(n)\right]}^T$ at time n, which is the input of the reservoir at this time. The output of the reservoir well be ${\mathit{y}}_{in}(n)=\mathit{U}(n+\lambda )$, and the prediction of different time steps is realized according to the adjustment of $\lambda$. The reservoir receives two inputs. The one is $\mathit{U}(n)$ from the input layer and the other is $\mathit{x}(n-1)$ from the previous state of the reservoir. The matrixes $\mathit{W}$ and ${\mathit{W}}_{in}$ are randomly generated and remain unchanged during the training. The weight $\mathit{W}$ is a large-scale sparse matrix, in which non-zero elements indicate activated neurons in the reservoir. The update formula of current state of reservoir is as follows:

$$\mathit{x}(n+1)=\tanh \left({\mathit{W}}_{in}\times \left[1;\mathit{U}(n+1)\right]+\mathit{W}\times \mathit{x}(n)\right)$$

(1)

where $\mathit{x}(n)\in {\mathbf{R}}^{L\times 1}$ is the state of the reservoir at time n, and its initial state is a zero vector. The operator $\tanh (\cdot )$ is the hyperbolic tangent activation function. With the infinite growth of time, the dependence of the current state affected by the initial state of reservoir gradually decreases or even disappears [25]. The information in reservoir is transferred to the output layer with a linear connection, and the output of the network is:

$$\mathit{y}(n)={\mathit{W}}_{out}\mathit{x}(n)$$

(2)

Afterwards, Jaeger et al. [26] improved the ESN and proposed the leaky integrator ESN (LI-ESN) by introducing leaky integrator neurons into the reservoir. LI-ESN is a variant of ESN, but its reservoir contains the leak integrator neurons. The update formula of the state of reservoir in LI-ESN is:

$$\mathit{x}(n+1)=a\times \tanh \left({\mathit{W}}_{in}\times \left[1;\mathit{U}(n+1)\right]+\mathit{W}\times \mathit{x}(n)\right)+(1-a)\times \mathit{x}(n)$$

(3)

where $a\in [0,1]$ is the leaking rate. The leaky neuron in reservoir has a leaky integration of its activation, partially remembers its previous activation, and retains the effect of the previous moment on the current state.

3. PM_2.5 Time Series Prediction Model

In order to accurately predict the concentration of PM_2.5, this paper proposes a PM_2.5 time series prediction model. To handle high-dimensional input variables, the mRMR method is utilized to select optimal subset. Then, the phase space is reconstructed from the subset to obtain evolution information of time series. Finally, we predict PM_2.5 time series based on supplemental leaky integrator echo state network (SLI-ESN) model, which is proposed in this paper.

3.1. mRMR Feature Selection Method

In this paper, mRMR method [20] is used for feature selection. Its basic idea is to maximize the correlation between input features and output, and to minimize the correlation between input features. The goal is to find the most representative feature subset. The evaluation of correlation between features is performed by mutual information method, and the appropriate feature subset is selected from original feature set. The maximum relevance and minimum redundancy of mRMR are calculated as follows:

$$\max D\left(S,C\right)=\frac{1}{\left|S\right|}{\displaystyle \sum _{f_i\in S}I\left(f_i;C\right)}$$

(4)

$$\min R\left(S\right)=\frac{1}{{\left|S\right|}^2}{\displaystyle \sum _{f_i,f_k\in S}I\left(f_i;f_k\right)}$$

(5)

where $f_i$ and $f_k$ represent the ith and kth feature in set S, respectively. C denotes the target output. $I\left(f_i;C\right)$ denotes the correlation function of the ith feature and output C. $I\left(f_i;f_k\right)$ represents the correlation function of the ith and kth feature. Mutual information is defined as follows:

$$I\left(x;y\right)={\displaystyle \iint P\left(x,y\right)\log \frac{P\left(x,y\right)}{P\left(x\right)P\left(y\right)}dxdy}$$

(6)

Therefore, the expression of mRMR is as follows:

$$\underset{S}{\max }\left[\frac{1}{\left|S\right|}{\displaystyle \sum _{f_i\in S}I\left(f_i;C\right)-\frac{1}{{\left|S\right|}^2}{\displaystyle \sum _{f_i,f_k\in S}I\left(f_i;f_k\right)}}\right]$$

(7)

Applying the forward selection algorithm to solve the objective function (7), we can get the sort result of the input features. Finally, the optimal input feature subset can be obtained by cross-validation.

3.2. Phase Space Reconstruction

Phase space reconstruction is a very important step in time series analysis and prediction. In order to extract useful information from time series, Takens [27] proposed embedding theorem of phase space reconstruction. In the actual calculation, since the numerical differentiation is sensitive to the error. Therefore, phase space reconstruction of time series generally adopts coordinate delay method. The essence is to construct m-dimensional vector by different delay times of one-dimensional time series. Its calculation formula is as follows:

$${\mathit{U}}_i(n)=\left[u_i(n),u_i(n+{\tau }_i),\cdots ,u_i\left(n+(m_i-1){\tau }_i\right)\right]$$

(8)

where $\left\{u_i(n)|i=1,2,\cdots ,K\right\}$ are input time series, $m_i$ is the embedding dimension, and ${\tau }_i$ is the delay time. The reconstructed phase space contains all the evolution information of original system.

The embedding theorem shows that for one-dimensional time series of infinitely long, noise-free m-dimensional chaotic attractors, an m-dimensional phase space can be found, as long as the embedding dimension satisfies $m\ge 2d+1$, where d is dimension of the dynamic system [28]. However, the existing time series are finite-length sequences with noise. The embedding dimension and delay time cannot be arbitrarily selected, otherwise the phase space quality of the reconstruction will be seriously affected.

3.3. Supplementary Leaky Integrator Echo State Network

The echo state network, including the model proposed in this paper, has the ability to memorize historical information. However, due to the unpredictability of chaotic time series, it is difficult to accurately retain the information before the current state of the reservoir for a long time. Moreover, due to the complexity of PM_2.5 time series, its current concentration information may be more dependent on the most recent state. Inspired by the leaked integrator reservoir, the update of the reservoir state is strongly dependent on the historical state. In this paper, considering the impact of previous state on current state, the state update formula of reservoir neurons is improved, and SLI-ESN model is proposed. The state update formula of reservoir is as follows:

$$\mathit{x}(n+1)=a\times tanh\left({\mathit{W}}_{in}\times [1;\mathit{U}(n+1)]+\mathit{W}\times \mathit{x}(n)\right)+b\times \mathit{x}(n)+(1-a-b)\times \mathit{x}(n-1)$$

(9)

As with LI-ESN, the meaning of the attenuation parameter $a$ has not changed, and its satisfy $0<a<1$. The parameter $b$ is the supplement factor, and it also need to satisfy $0<b<1$. More importantly, like the side length of a triangle, each coefficient must be greater than 0, and satisfy $a+b<1$. Because the current state of reservoir is more dependent on the state of neighboring moment, this paper considers influence of the first two historical moments. When $n\ge 2$, supplement $(1-a-b)\times \mathit{x}(n-1)$ is established. Set ${\mathit{W}}_{back}=0$, the output layer is calculated as follows:

$$\mathit{y}(n)={\mathit{W}}_{out}\mathit{x}(n)$$

(10)

The calculation from reservoir to output layer is solved by ridge regression [29]. It is a biased estimation regression method for collinear data analysis. It overcomes ill-posed problem in least squares solution and prevents over-fitting problem. The expression is as follows:

$${\mathit{W}}_{out}={({\mathit{X}}^T\mathit{X}+k\mathit{I})}^{-1}{\mathit{X}}^T\mathit{Y}$$

(11)

where $k$ is regularization parameter, $\mathit{I}$ is the identity matrix, $\mathit{x}(n)$ and $\mathit{y}(n)$ are column vectors of $\mathit{X}$ and $\mathit{Y}$, respectively.

3.4. Algorithm Flow

In order to accurately predict the PM_2.5 time series, this paper proposes the prediction structure based on feature extraction and improved echo state network model. The process is summarized as follows.

• Feature selection: select the optimal subset from original dataset based on mRMR feature selection method.
• Phase space reconstruction: reconstruct phase space of the selected optimal subset based on Takens’ theorem and form a new set of input features.
• Data division: divide training set and testing set according to a certain proportion.
• Model training: train SLI-ESN model using ridge regression algorithm on training set.
• Prediction: predict PM_2.5 time series using SLI-ESN model on testing set.

In order to briefly explain prediction process of the proposed model, the schematic diagram of the method is shown in Figure 2.

4. Results and Discussion

In order to prove the validity and practicability of the proposed method, the SLI-ESN model was applied to predict the actual observed PM_2.5 time series of Beijing. At the same time, this paper also carried out comparative experiments of ESN, LI-ESN [26], extreme learning machine (ELM) [30], hierarchical ELM (H-ELM) [31] and stacked auto-encoder (SAE) [32]. The experimental environment is Windows 7 system, MATLAB 2016a programming software. The computer memory is 6 GB, clocked at 3.50 GHz, and Intel-i3 CPU.

Firstly, the AQI dataset and evaluation indicators are explained. Secondly, the feature selection experiment based on mRMR is introduced in detail, and the optimal subset is obtained. Finally, the reconstructed features of optimal subset is used for training the SLI-ESN model, and the validity of the proposed method in PM_2.5 time series prediction is verified by predictive indicators.

4.1. Data Description

This paper selects 8759 samples of hourly air pollution data from January to December 2016 in Haidian District, Beijing. The dataset is from the US Embassy (Harvard University Geographic Analysis Center Dataverse), including the average concentration of PM_2.5, PM₁₀, NO₂, CO, O₃, and SO₂ per hour, and hourly temperature (T), pressure (P), humidity (H), wind speed (WS) and wind direction (WD).

Because PM_2.5 is a particulate matter that can enter the lungs, it has a great impact on human health and the quality of atmospheric environment. Therefore, predicting PM_2.5 concentration can not only make an estimate of environmental quality, but also have a greater impact on better environmental governance and human health. This paper uses five predictive indicators to evaluate the prediction results of PM_2.5 time series, namely root mean square error (RMSE), normalized root mean square error (NRMSE), mean absolute error (MAE), symmetric mean absolute percentage error (SMAPE) and Pearson correlation coefficient (R).

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{t=1}^N{\left[\widehat{y}\left(t\right)-y\left(t\right)\right]}^2}}$$

(12)

$$\mathrm{NRMSE}=\frac{1}{{\widehat{y}}_{\max }-{\widehat{y}}_{\min }}\sqrt{\frac{1}{N}{\displaystyle \sum _{t=1}^N{\left[\widehat{y}\left(t\right)-y\left(t\right)\right]}^2}}$$

(13)

$$\mathrm{MAE}=\frac{1}{N}{\displaystyle \sum _{t=1}^N\left|\widehat{y}\left(t\right)-y\left(t\right)\right|}$$

(14)

$$\mathrm{SMAPE}=\frac{1}{N}{\displaystyle \sum _{t=1}^N\frac{\left|y\left(t\right)-\widehat{y}\left(t\right)\right|}{\left(\left|\widehat{y}\left(t\right)\right|+\left|y\left(t\right)\right|\right)/2}}$$

(15)

$$\mathrm{R}=\frac{{\displaystyle \sum _{t=1}^N\left(y\left(t\right)-\overline{y}\left(t\right)\right)\left(\widehat{y}\left(t\right)-\overline{\widehat{y}}\left(t\right)\right)}}{\sqrt{{\displaystyle \sum _{t=1}^N{\left(y\left(t\right)-\overline{y}\left(t\right)\right)}^2}{\displaystyle \sum _{t=1}^N{\left(\widehat{y}\left(t\right)-\overline{\widehat{y}}\left(t\right)\right)}^2}}}$$

(16)

where N is the number of samples, $\widehat{y}\left(t\right)$ is the predicted output, and $y\left(t\right)$ is the target value.

In the above evaluation indicators, the smaller the values of the estimated indicators of RMSE, NRMSE, MAE, and SMAPE are, the better the prediction results of the model are. $\mathrm{R}=1$ indicates that $\widehat{y}\left(t\right)$ and $y\left(t\right)$ are linear correlation, and $\mathrm{R}=0$ indicates there is no correlation. When $\mathrm{R}\in \left(0,1\right)$, it indicates that there is a correlation, and the larger value is, the stronger the linear correlation is.

4.2. Data Processing

The optimal subset selection of original data was performed using mRMR. In this paper, the first 75% of the dataset is used for training set, and the last 25% is used for testing set. The optimal subset is selected based on training set, and each model parameter in the simulation is obtained only from the training set.

Firstly, PM_2.5 time series is chosen as reference variable. Other 10-dimensional variables are used as comparison variables. The data was quantitatively analyzed using mRMR to obtain sorting results: PM_2.5, CO, WS, PM₁₀, H, WD, SO₂, NO₂, P, T, and O₃. According to the influence of different factors on PM_2.5 concentration, the correlation between different variables and PM_2.5 is shown in Figure 3.

According to the mRMR method, the irrelevant and redundant variables in the original dataset are reduced. The selection result of obtained optimal subset is shown in Figure 4.

According to the prediction results in Figure 4, when the predictor is with 5 dimensions, the prediction error is the smallest, namely 9.199. Therefore, the optimal subset is PM_2.5, CO, WS, PM₁₀, and H. The phase space reconstruction is then performed on the selected optimal subset.

The delay time τ and embedding dimension m calculated by C-C method [33] are shown in

Table 1. Phase space reconstruction parameters of variables. Abbreviations: hourly temperature (T), pressure (P), humidity (H), wind speed (WS), and wind direction (WD).

Variables	PM2.5	PM10	NO2	CO	O3	SO2	T	P	H	WS	WD
$\tau$	8	8	4	6	4	6	4	12	4	4	6
m	2	2	3	2	4	2	6	2	4	4	3

. The 5 variables in the optimal subset are expressed as bold in the table. As shown in Table 1, the delay time for the obtained optimal subset is [8,8,6,4,4], and the embedding dimension is [2,2,2,4,4] for PM_2.5, PM₁₀, CO, H and WS respectively.

4.3. Experimental Results and Analysis

Phase space reconstruction is performed on the optimal subset to obtain a 14-dimensional reconstructed time series, which are used for inputs of prediction model. For the LI-ESN model, the reasonable range of the leaking rate a is (0, 1). For the SLI-ESN model, we use the cross-validation to select the two parameters a and b. According to experience, the feasible and effective range of the supplementy factor b is mainly within (0, 0.1). In this paper, ESN, LI-ESN, ELM, H-ELM and SAE are selected as comparison methods. The specific one-step (1 hour) prediction results are shown in

Table 2. Comparison of one-step (1 hour) prediction results. Abbreviations: root mean square error (RMSE), normalized root mean square error (NRMSE), mean absolute error (MAE), symmetric mean absolute percentage error (SMAPE), Pearson correlation coefficient (R), echo state network (ESN), leaky integrator ESN (LI-ESN), extreme learning machine (ELM), hierarchical ELM (H-ELM), stacked auto-encoder (SAE), supplementary leaky integrator echo state network (SLI-ESN).

Methods	RMSE	NRMSE	MAE	SMAPE	R
ESN	10.2020	0.0160	6.6948	0.1053	0.9936
LI-ESN	9.7063	0.0152	6.2322	0.1001	0.9943
ELM	11.7330	0.0184	7.1902	0.1082	0.9914
H-ELM	14.1520	0.0222	8.0575	0.1102	0.9876
SAE	32.1700	0.0505	20.1840	0.2764	0.9448
SLI-ESN	9.3953	0.0147	5.8447	0.0894	0.9945

.

It can be seen from Table 2 that the proposed method achieves better prediction results in one-step (1 hour) prediction. The one-step (1 hour) prediction result of PM_2.5 concentration is shown in Figure 5. And Figure 6 plots the fit of predicted value to actual data. It can be seen from the figures that the prediction has a good linear relationship with actual value. SLI-ESN performs satisfactorily at peaks and undulating moments, which mainly depends on the full utilization of the historical state of the reservoir and the effective information obtained by phase space reconstruction.

At the same time, the simulation results of the five-step (5 hours) prediction are given in Figure 7. The prediction curve can better track original input, and the medium-term prediction effect is also good.

Table 3. Comparison of five-step (5 hours) prediction results.

Methods	RMSE	NRMSE	MAE	SMAPE	R
ESN	43.1243	0.0677	29.3879	0.3777	0.8907
LI-ESN	41.1574	0.0646	27.6611	0.3577	0.9027
ELM	45.7617	0.0718	29.7892	0.3726	0.8678
H-ELM	49.9159	0.0783	32.1977	0.3974	0.8368
SAE	51.9847	0.0816	34.5487	0.4039	0.8394
SLI-ESN	37.6874	0.0591	25.5871	0.3392	0.9108

gives the five-step (5 hours) prediction results. The ten-step (10 hours) prediction result of SLI-ESN is shown in Figure 8. As seen in Figure 8, in some of the peaks, the prediction curve is still able to roughly fit the fluctuation trend of the original data, which is precisely because SLI-ESN makes full use of the role of historical information. Its comparison with other algorithms is shown in

Table 4. Comparison of ten-step (10 hours) prediction results.

Methods	RMSE	NRMSE	MAE	SMAPE	R
ESN	72.3723	0.2299	50.5928	0.5620	0.6858
LI-ESN	70.5700	0.2944	49.3628	0.5569	0.6895
ELM	71.6619	0.2429	47.1085	0.5385	0.6582
H-ELM	66.7932	0.2760	46.1529	0.5003	0.7053
SAE	71.7758	0.3631	49.7844	0.5617	0.6623
SLI-ESN	65.7108	0.1966	46.3633	0.5443	0.7314

. The proposed algorithm achieves the optimal value on all four parameter indicators, except for SMAPE, which fully demonstrates the effectiveness of SLI-ESN in long-term prediction.

Longitudinal comparison between Table 2, Table 3 and Table 4 shows that the longer the prediction time is, the larger the error will be, which is consistent with the basic characteristics of the chaotic time series. Comparing Table 2, Table 3 and Table 4, it can be found that the recurrent neural networks in three kinds of networks (SAE as the representative of deep learning, ELM and H-ELM as the representative of feedforward neural network and ESN as the representative of recurrent neural network) have more satisfactory prediction performance. This exactly demonstrates the validity of the structure of the reservoir in time series prediction.

In order to further illustrate the performance of the proposed method, the running time of all comparison methods are shown in

Table 5. Comparison of running time on six methods (second).

Methods	One-step (1 hour)		Five-step (5 hours)		Ten-step (10 hours)
	Training Time	Testing Time	Training Time	Testing Time	Training Time	Testing Time
ESN	0.1145	0.0238	0.1138	0.0213	0.1281	0.0226
LI-ESN	0.1139	0.0225	0.1648	0.0307	0.5086	0.0418
ELM	0.0624	0.0312	0.0624	0.0312	0.0624	0.0312
H-ELM	0.4292	0.1039	0.4894	0.1160	0.1748	0.0696
SAE	143.1321	0.0292	139.4019	0.0268	145.4450	0.0816
SLI-ESN	3.1061	0.1730	3.3148	0.1669	4.3705	0.2193

. The results show that the training time based on the ESN and ELM model is much smaller than the deep learning model. This is mainly because the training process of deep learning consumes a lot of time, and this time-consuming process is not required by the other neural networks. Moreover, SLI-ESN can complete training and testing within an acceptable time frame. It indicates that SLI-ESN can achieve good prediction results in both prediction accuracy and time consuming.

5. Conclusions

PM_2.5 is the main source of air pollution, and predicting PM_2.5 concentration is of great significance for protecting the environment. In order to improve prediction accuracy and reliability, it is important to preprocess the data to eliminate irrelevant and redundant variables before prediction. In this paper, the mRMR is used to screen original dataset to obtain optimal subset. Phase space reconstruction is performed on the optimal subset. And the reconstructed data is used as new input time series of SLI-ESN model for prediction. Experiments show the validity of the SLI-ESN model, which has high prediction accuracy in medium and long term projects, good generalization performance and good application prospects.

Although this paper has achieved the desired results, there are still some issues that need to be addressed in future work. First of all, long-term predictions are not satisfactory. We want to extend the model to longer prediction interval, such as one day, one week or one month. In addition, optimal subset selection and model optimization take a lot of time. In the future, we expect to simultaneously implement input variable selection and model optimization based on an optimization algorithm. Optimization objects include, but are not limited to, input variables, model structure, and model parameters.