A Novel Linear Time-Varying GM(1,N) Model for Forecasting Haze: A Case Study of Beijing, China

Abstract

Haze is the greatest challenge facing China’s sustainable development, and it seriously affects China’s economy, society, ecology and human health. Based on the uncertainty and suddenness of haze, this paper proposes a novel linear time-varying grey model (GM)(1,N) based on interval grey number sequences. Because the original GM(1,N) model based on interval grey number sequences has constant parameters, it neglects the dynamic change characteristics of parameters over time. Therefore, this novel linear time-varying GM(1,N) model, based on interval grey number sequences, is established on the basis of the original GM(1,N) model by introducing a linear time polynomial. To verify the validity and practicability of this model, this paper selects the data of PM₁₀, SO₂ and NO₂ concentrations in Beijing, China, from 2008 to 2018, to establish a linear time-varying GM(1,3) model based on interval grey number sequences, and the prediction results are compared with the original GM(1,3) model. The result indicates that the prediction effect of the novel model is better than that of the original model. Finally, this model is applied to forecast PM₁₀ concentration for 2019 to 2021 in Beijing, and the forecast is made to provide a reference for the government to carry out haze control.

1. Introduction

In recent years, the process of urbanization and industrialization has been accelerating in China, but air pollution is seriously increasingly. At present, haze is one of the greatest atmospheric, environmental pollution problems in China. Haze not only has adverse effects on the ecological environment [1,2,3] but also poses a major threat to human health [4,5]. As one of the largest developing countries in the world, China is facing the biggest challenge of sustainable development. In 2018, China’s government officially issued the ‘Three-Year Plan of Action to Win the Blue Sky Defense War’. Therefore, an accurate study of haze is of great importance.

Beijing is one of key cities for haze control in China. 16 haze events occurred in Beijing from November, 2012 to January, 2013, and the minimum visibility was even 667 m during this period [6]. Besides, the haze event in January, 2013 may have been the cause of 690 deaths in Beijing, which could lead to 253.8 million dollars losses [7]. To reduce the haze event, China’s government implemented the ‘Air Pollution Prevention and Control Action Plan’ in September, 2013. Since then, an increasing number of scholars started to research haze from different aspects. Current haze studies mainly focused on health impact [3,4,7], economic loss [7], chemical composition [8], statistical characteristics [9], trend prediction [10], formation mechanism [11], and the rest. Above researches mainly based on hourly or daily data, but the goal of Blue Sky Defense War is to reduce the average annual concentration of pollutants. Therefore, using annual data to forecast haze is necessary. Due to the uncertainty and suddenness of haze, the data of its related indices are uncertain. In the system research, grey number refers to the uncertain number in an interval or a general number set [12]. In other words, due to the constraints of data acquisition tools, acquisition conditions and errors by the acquisition personnel, the value of haze may contain inevitable measurement errors in a certain range, thus, haze has obvious grey number features. Therefore, Wu et al. established a grey prediction model to forecast annual pollutant concentrations effecting haze by limited data for the first time [13]. Grey prediction model is an important part of grey system. As an emerging uncertain system theory, grey system theory is characterized by small data modeling to obtain accurate results, and widely applied in all professions and trades. The grey system theory was established by Professor Deng, and assists with solving an uncertain system of small data and poor information, mainly through the deep mining of the inherent law in the existing information research system. Besides, grey number is the basic unit of grey system [12]. To summarize, grey prediction model is an effective tool for studying haze.

At present, the original grey model (GM)(1,N) based on interval grey number sequences is used to estimate structural parameters by the least squares method, and the result of its structural parameters are fixed values that are not related to time. Although this method is simple, it ignores the dynamic change characteristics of the parameters over time, which may lead the model to have a high precision fit, but not an ideal predictive effect. Therefore, this paper selects the data of PM₁₀ concentration (PM₁₀ refers to particles with aerodynamic equivalent diameter less than or equal to 10 microns in ambient air, called inhalable particulate matter), SO₂ concentration and NO₂ concentration in Beijing, China, from 2008 to 2018, and establishes a novel linear time-varying GM(1,3) model based on interval grey number sequences to perform simulations and make predictions. Meanwhile, we establish an original GM(1,3) model based on interval grey number sequences as the contrast model, and then we compare the results of these two models.

This paper is organized as follows: the literature reviews of studying on haze and grey prediction model are presented in Section 2; the modeling algorithms and model testing methods of the linear time-varying GM(1,N) model based on interval grey number sequences are illustrated in Section 3; the linear time-varying GM(1,N) model based on interval grey number sequences is established to forecast PM₁₀ concentration of Beijing in Section 4; the main conclusions of this paper are summarized in Section 5.

2. Literature Reviews

2.1. Study of Haze

To control haze more scientifically and effectively, many scholars at home and abroad have carried out a series of extensive studies. PM_2.5 (PM_2.5 refers to particles with an aerodynamic equivalent diameter less than or equal to 2.5 microns in ambient air, called fine particulate matter) and PM₁₀ are main factors causing haze, which not only affect the quality of atmospheric environment and visibility, but also endanger human health. Voukantsis et al. used the neural network method of multilayer perceptive structure to predict the daily change of PM₁₀ concentration in Thessaloniki and Helsinki [14]. Kumar et al. used a principal component regression method to predict short-term air quality index days in Delhi [15]. Lang et al. generalized the season autoregressive integrated moving average (SARIMA) model to the small-scale time sequences and used this model to make a short-term prediction of PM_2.5 concentration at 10 stations in Hangzhou [16]. Mishra et al. studied the meteorological factors and atmospheric pollutants affecting PM_2.5 in Delingha, and predicted the PM_2.5 concentration under haze conditions in Delingha using a combination of a neural network and fuzzy logic [17]. Konovalov et al. predicted future PM₁₀ values in Europe by combining the deterministic prediction method and the linear regression method [18]. Based on the data of PM_2.5 concentration in 285 cities in China from 2001 to 2012, Cheng et al. used a dynamic space panel model to analyze the main driving factors affecting air pollution [19]. Wu et al. established an input-oriented zero sum gains-data envelopment analysis (ZSG-DEA) model, and studied the allocation efficiency of PM_2.5 emissions in China’s provinces under fixed target amount conditions [20]. The novel air quality early warning system was proposed by Xu et al., which included an evaluation module, a prediction module and a feature estimation module, and this literature used a novel dynamic fuzzy comprehensive evaluation method to determine the air quality level and major pollutants in the study area [21]. Han et al. predicted the spatial distribution of lung cancer in males induced by PM_2.5 in China from 2010 to 2015, and the spatial autocorrelation method was used to evaluate the spatial relationship between the incidence of lung cancer and the atmospheric level of satellite-derived PM_2.5 from 2006 to 2009 [22].

Nowadays, there are some scholars studying haze with grey system theory. Xiong et al. established the multivariate grey model (MGM) (1,M) based on interval grey number sequences to predict the visibility and relative humidity during a haze period in Nanjing [23]. Gong et al. combined the GM(1,1) model with both the Markoff chain model and the residual error correction model to establish a modified grey Markoff chain model, and used this model to predict the PM_2.5 concentration in Shanghai [24]. Wang studied the distribution characteristics of PM_2.5 concentration in Huaian with the non-parametric hypothesis test, and predicted the PM_2.5 concentration of Huaian in the next five years by using GM(1,1) model [25]. Chen et al. predicted the hourly PM_2.5 and PM₁₀ in Taichung’s Dali area by using several GM(1,1) models and back-propagation artificial neural network, and compared their predictive performances [26]. Wang et al. established a novel grey correlation degree model, and used this model to dynamically analyze the influencing factors of haze in southern China [27].

2.2. Study of Grey Prediction Model

The grey prediction model has been widely applied in many fields, and it occupies an important position in the grey system theory. The GM(1,1) model that has only one variable is the most widely used prediction model in the grey system, but it is most suitable for the time sequences of monotonic increasing or monotonic decreasing. Hence, many scholars have worked to improve it. Wu et al. proposed a new GM(1,1) model with the fractional order accumulation, and this model had a better predicted performance than the traditional model [28]. Focusing on why the discrete grey model simulated the constant value growth rate, Zhang et al. established a linear time-varying discrete grey model by introducing a linear time polynomial [29]. Wang introduced the time polynomial function into the GM(1,1) power model and optimized the power exponent of the model [30]. Zeng et al. built a self-adapting intelligent grey prediction model to predict the natural gas demand of China [31].

However, the GM(1,1) model ignores the effect of related factors on the system behavior data, and then, Professor Deng proposed the GM(1,N) model, which has one system behavior sequence and N-1 related factors sequences [12]. The GM(1,N) model, as an extension of the GM(1,1) model, fully considers the impact of related factors on the system behavior data. Therefore, the GM(1,N) model and its optimization models are increasing in researches and applications. Ding et al. established a novel GM(1,N) model combined with the changing trend of the driving term, and this model was applied to predict CO₂ emissions from fuel combustion in China [32]. Zeng et al. proposed the optimal background-value GM(1, N) model through optimizing the background-value coefficient with the particle swarm optimization algorithm [33]. Wang et al. constructed a nonlinear GM(1,N) model by introducing the power exponent of the related factors, and this model was used to predict the carbon emissions of fossil energy consumption in China [34]. The above studies on grey prediction models were mainly based on real number sequence. Recently, increasingly more scholars have begun to explore the modeling problem of grey prediction models based on interval grey number sequences. Ye et al. fully explored and expanded the axiom of generalized non-decreasing grey degrees and established a prediction model for interval grey number sequences [35]. Luo et al. established a discrete GM(1,1) model for kernel and measure sequences, and then restored the predicted value of the interval grey number [36]. Yang et al. established a prediction model for the normal distribution based on interval grey number sequences in the context of the normal distribution of uncertain information [37]. A new nonlinear GM(1,N) model was established based on interval grey number sequences by Xiong et al., which was applied to the prediction of air quality index in haze period [38].

3. Methodology

3.1. The Introduction of Interval Grey Number

In this section, we will mainly introduce the basic concepts of interval grey number, including the definitions of interval grey number, kernel and grey radius.

Definition 1. 

[12] grey number that has both lower bound $a_k$ and upper bound $b_k$ is called interval grey number, denoted as ${\otimes }_k\in \left[a_k,b_k\right]$.

Definition 2. 

[12] suppose that interval grey number ${\otimes }_k$ is a continuous function, $\tilde{\otimes }=(a_k+b_k)/2$ is called a kernel of interval grey number ${\otimes }_k$.

Definition 3. 

[12] when ${\otimes }_k$ is a continuous interval grey number, $r\left(k\right)=(b_k-a_k)/2$ is called a grey radius of interval grey number ${\otimes }_k$.

3.2. Linear Time-Varying GM(1,N) Model Based on Interval Grey Number Sequences

In this section, the modeling mechanism of the linear time-varying GM(1,N) model based on interval grey number sequences will be introduced. Besides, this model will be constructed from the kernel and grey radius sequences, respectively. Finally, the kernel and grey radius sequences will be restored to the interval grey number sequences. To illustrate the model more clearly, we will show the modeling steps of the linear time-varying GM(1,N) model based on interval grey number sequences in Figure 1.

Then, this section will mainly introduce the linear time-varying GM(1,N) model based on kernel sequences as follows:

Definition 4. 

Suppose that the system behavior sequence is ${\tilde{\otimes }}_1^{\left(0\right)}=\left({\tilde{\otimes }}_1^{\left(0\right)}\left(1\right),{\tilde{\otimes }}_1^{\left(0\right)}\left(2\right),\cdots ,{\tilde{\otimes }}_1^{\left(0\right)}\left(n\right)\right)$, and the related factor sequences is ${\tilde{\otimes }}_i^{\left(0\right)}=\left({\tilde{\otimes }}_i^{\left(0\right)}\left(1\right),{\tilde{\otimes }}_i^{\left(0\right)}\left(2\right),\cdots ,{\tilde{\otimes }}_i^{\left(0\right)}\left(n\right)\right),i=2,3,\cdots ,N,$ where ${\tilde{\otimes }}_i^{\left(1\right)}$ is the first order accumulating generation sequence of ${\tilde{\otimes }}_i^{\left(0\right)}, i=1,2,\cdots ,N$; $z_1^{\left(1\right)}\left(k\right)=0.5\left({\tilde{\otimes }}_1^{\left(1\right)}\left(k\right)+{\tilde{\otimes }}_1^{\left(1\right)}\left(k-1\right)\right)$ is the mean sequence generated by consecutive neighbors of ${\tilde{\otimes }}_1^{\left(1\right)}$; thus, the linear time-varying GM(1,N) model based on kernel sequences is shown as follows:

$${\tilde{\otimes }}_1^{\left(0\right)}\left(k\right)+a{z}_1^{\left(1\right)}\left(k\right)={\displaystyle \sum }_{i=2}^N\left(b_{i1}+b_{i2}k\right){\tilde{\otimes }}_i^{\left(1\right)}\left(k\right).$$

(1)

Also, the whitening equation of the linear time-varying model based on kernel sequences is shown as follows:

$$\frac{d{\tilde{\otimes }}_1^{\left(1\right)}}{dt}+a{\tilde{\otimes }}_1^{\left(1\right)}={\displaystyle \sum }_{i=2}^N(b_{i1}+b_{i2}t){\tilde{\otimes }}_i^{\left(1\right)}.$$

(2)

Theorem 1. 

The undetermined coefficient vector of the linear time-varying GM(1,N) model based on kernel sequences is $\widehat{\mathit{a}}= {\left[a,b_{21},···,b_{N1},b_{22},···,b_{N2}\right]}^T$, which can be estimated by utilizing the least squares method as follows:

(1) When $n-1=2\mathrm{N}-1$, that is $n=2\mathrm{N}$, and $\left|\mathit{B}\right|\ne 0$, $\widehat{\mathit{a}}={\mathit{B}}^{-1}\mathit{Y}$;

(2) When $n-1>2\mathrm{N}-1$, that is $n>2\mathrm{N}$, and $\left|{\mathit{B}}^T\mathit{B}\right|\ne 0$, $\widehat{\mathit{a}}={\left({\mathit{B}}^T\mathit{B}\right)}^{-1}{\mathit{B}}^T\mathit{Y}$;

(3) When $n-1<2\mathrm{N}-1$, that is $n<2\mathrm{N}$, and $\left|{\mathit{B}}^T\mathit{B}\right|\ne 0$, $\widehat{\mathit{a}}={\mathit{B}}^T{\left(\mathit{B}{\mathit{B}}^T\right)}^{-1}\mathit{Y}$; where $\mathit{B}=\left[\begin{array}{ccccccc}-z_1^{\left(1\right)}\left(2\right)& {\tilde{\otimes }}_2^{\left(1\right)}\left(2\right)& \dots & {\tilde{\otimes }}_N^{\left(1\right)}\left(2\right)& 2{\tilde{\otimes }}_2^{\left(1\right)}\left(2\right)& \dots & 2{\tilde{\otimes }}_N^{\left(1\right)}\left(2\right)\\ -z_1^{\left(1\right)}\left(3\right)& {\tilde{\otimes }}_2^{\left(1\right)}\left(3\right)& \dots & {\tilde{\otimes }}_N^{\left(1\right)}\left(3\right)& 3{\tilde{\otimes }}_2^{\left(1\right)}\left(3\right)& \dots & 3{\tilde{\otimes }}_N^{\left(1\right)}\left(3\right)\\ ⋮& ⋮& & ⋮& ⋮& & ⋮\\ -z_1^{\left(1\right)}\left(n\right)& {\tilde{\otimes }}_2^{\left(1\right)}\left(n\right)& \dots & {\tilde{\otimes }}_N^{\left(1\right)}\left(n\right)& n{\tilde{\otimes }}_2^{\left(1\right)}\left(n\right)& \dots & n{\tilde{\otimes }}_N^{\left(1\right)}\left(n\right)\end{array}\right],\mathit{Y}=\left[\begin{array}{c}{\tilde{\otimes }}_1^{\left(0\right)}\left(2\right)\\ {\tilde{\otimes }}_1^{\left(0\right)}\left(3\right)\\ ⋮\\ {\tilde{\otimes }}_1^{\left(0\right)}\left(n\right)\end{array}\right]$

The proof is similar to literature [39].

Theorem 2. 

After calculating the coefficient vector $\widehat{\mathit{a}}$, the solution of the linear time-varying GM(1,N) model based on kernel sequences is shown as follows:

$${\widehat{\tilde{\otimes }}}_1^{\left(1\right)}\left(k+1\right)=\frac{1}{a}{\displaystyle \sum }_{i=2}^N(b_{i1}+b_{i2}k){\tilde{\otimes }}_i^{\left(1\right)}\left(k+1\right)+e^{-ak}\left[{\tilde{\otimes }}_1^{\left(1\right)}\left(0\right)-\frac{1}{a}{\displaystyle \sum }_{i=2}^N(b_{i1}+b_{i2}k){\tilde{\otimes }}_i^{\left(1\right)}\left(k+1\right)\right],$$

(3)

where the solution can be obtained using the initial condition ${\tilde{\otimes }}_1^{\left(1\right)}\left(0\right)= {\tilde{\otimes }}_1^{\left(0\right)}\left(1\right)$. Also, the inverse accumulating reduction equation is shown as follows:

$${\widehat{\tilde{\otimes }}}_1^{\left(0\right)}\left(k+1\right)={\widehat{\tilde{\otimes }}}_1^{\left(1\right)}\left(k+1\right)-{\widehat{\tilde{\otimes }}}_1^{\left(1\right)}\left(k\right).$$

(4)

The modeling mechanism of the linear time-varying GM(1,N) model based on grey radius sequences are the same as the linear time-varying GM(1,N) model based on kernel sequences, thus we will not present the linear time-varying GM(1,N) model based on grey radius sequences repeatedly.

After obtaining the values of kernel and grey radius respectively, we will calculate the upper and lower bounds of the interval grey number sequences as follows [40]:

$$\{\begin{array}{c}{\widehat{a}}_{k+1}={\widehat{\tilde{\otimes }}}_1^{\left(0\right)}\left(k+1\right)-{\widehat{r}}_1^{\left(0\right)}\left(k+1\right)\\ {\widehat{b}}_{k+1}={\widehat{\tilde{\otimes }}}_1^{\left(0\right)}\left(k+1\right)+{\widehat{r}}_1^{\left(0\right)}\left(k+1\right)\end{array}.$$

(5)

3.3. Model Evaluation Criterion

To analyze the reliability and credibility of the prediction model, we will show the model evaluation criterion for testing the model accuracy in this section. By comparing the relative error and average relative error of the upper and lower bounds of the interval grey number sequences to test the prediction model, the testing equations are shown as follows:

The relative error of the upper and lower bounds of the interval grey number sequences are shown as follows:

$$\{\begin{array}{c}\underset{\_}{{\mathsf{\Delta }}_k}=\frac{\left|{\widehat{a}}_k-a_k\right|}{a_k}\times 100\%\\ \overline{{\mathsf{\Delta }}_k}=\frac{\left|{\widehat{b}}_k-b_k\right|}{b_k}\times 100\%\end{array},k=1,2,\cdots ,n.$$

(6)

The average relative error of the upper and lower bounds of the interval grey number sequences are shown as follows:

$$\{\begin{array}{c}\overline{{\widehat{a}}_k}=\frac{1}{n}{\displaystyle {\displaystyle \sum }_{k=1}^n}\underset{\_}{{\mathsf{\Delta }}_k}\\ \overline{{\widehat{b}}_k}=\frac{1}{n}{\displaystyle {\displaystyle \sum }_{k=1}^n}\overline{{\mathsf{\Delta }}_k} \end{array},k=1,2,\cdots ,n.$$

(7)

Prediction accuracy is an important criterion for measuring the reliability of the prediction model. Therefore, this paper provides the prediction accuracy corresponding to the average relative error in

Table 1. The average relative error criterion for testing model [41].

Average Relative Error	Prediction Accuracy
$<$10%	High prediction
10–20%	Good prediction
20–50%	Reasonable prediction
$>$50%	Weak prediction

.

4. Empirical Results and Discussion

In this section, the linear time-varying GM(1,N) model and the original GM (1,N) model, based on interval grey number sequences, will be established to simulate the development trend of haze in Beijing, and the model with high prediction accuracy will be selected to forecast the haze situation in Beijing from 2019 to 2021.

4.1. Data Selection and Processing

PM₁₀ can be formed by the interaction of sulfur oxides, nitrogen oxides and other compounds in the ambient air. Additionally, PM₁₀ is highly correlated with SO₂ and NO₂ in 31 cities of China [42]. Therefore, PM₁₀ concentration in Beijing from 2008 to 2018 is selected as the system behavior sequence, and SO₂ concentration and NO₂ concentration are selected as related factor sequences. In addition, the data are from the annual report of air quality in Beijing. In data processing, the maximum and minimum observed values over the previous three values are the upper and lower bounds of the third interval grey number. Moreover, the interval grey number sequences for 2008 to 2014 are used as the modeling data, and the interval grey number sequences for 2015 to 2018 are used as the prediction data. Besides, we denote the interval grey number sequences of PM₁₀ concentration, SO₂ concentration and NO₂ concentration as $X_1(\otimes )$, $X_2(\otimes )$, $X_3(\otimes )$ respectively. The original data is shown in

Table 2. The interval grey number sequences of PM₁₀ concentration, SO₂ concentration and NO₂ concentration

$\mathit{k}$	Year	${\mathit{X}}_1(\otimes )(\mathsf{\mu }\mathbf{g}/{\mathbf{m}}^3)$	${\mathit{X}}_2(\otimes )(\mathsf{\mu }\mathbf{g}/{\mathbf{m}}^3)$	${\mathit{X}}_3(\otimes )(\mathsf{\mu }\mathbf{g}/{\mathbf{m}}^3)$
1	2008	[122,161]	[36,53]	[49,66]
2	2009	[121,148]	[34,47]	[49,66]
3	2010	[114,122]	[32,36]	[49,55]
4	2011	[114,121]	[28,34]	[53,55]
5	2012	[109,121]	[28,32]	[52,55]
6	2013	[108,114]	[27,28]	[52,56]
7	2014	[108,116]	[22,28]	[52,59]
8	2015	[102,116]	[14,27]	[50,59]
9	2016	[100,116]	[10,22]	[48,59]
10	2017	[95,108]	[8,14]	[46,50]
11	2018	[91,100]	[6,10]	[42,48]

.

4.2. Establishment and Comparison of Model

Step 1: on the basis of the data from 2008 to 2014 in Table 2, the kernel and grey radius sequences of PM₁₀ concentration, SO₂ concentration and NO₂ concentration are calculated according to Definition 2 and 3, and the results are shown in

Table 3. The kernel and grey radius sequences of PM₁₀ concentration, SO₂ concentration and NO₂ concentration

$\mathit{k}$	${\tilde{\otimes }}_1\left(\mathit{k}\right)$	${\tilde{\otimes }}_2\left(\mathit{k}\right)$	${\tilde{\otimes }}_3\left(\mathit{k}\right)$	${\mathit{r}}_1\left(\mathit{k}\right)$	${\mathit{r}}_2\left(\mathit{k}\right)$	${\mathit{r}}_3\left(\mathit{k}\right)$
1	141.5	44.5	57.5	19.5	8.5	8.5
2	134.5	40.5	57.5	13.5	6.5	8.5
3	118	34	52	4	2	3
4	117.5	31	54	3.5	3	1
5	115	30	53.5	6	2	1.5
6	111	27.5	54	3	0.5	2
7	112	25	55.5	4	3	3.5

.

Step 2: after calculating model parameters by least squares method, a linear time-varying GM(1,3) model for the PM₁₀ concentration kernel sequence is established as follows:

$${\tilde{\otimes }}_1^{\left(0\right)}\left(k\right)+1.9972z_1^{\left(1\right)}\left(k\right)=\left(1.5419-0.0591k\right){\tilde{\otimes }}_2^{\left(1\right)}\left(k\right)+\left(3.7449-0.0014k\right){\tilde{\otimes }}_3^{\left(1\right)}\left(k\right)$$

(8)

Similarly, a linear time-varying GM(1,3) model for the PM₁₀ concentration grey radius sequence is established as follows:

$$r_1^{\left(0\right)}\left(k\right)+1.6900z_1^{\left(1\right)}\left(k\right)=\left(-6.9833+2.2742k\right)r_2^{\left(1\right)}\left(k\right)+\left(9.5542-2.0577k\right)r_3^{\left(1\right)}\left(k\right).$$

(9)

According to the linear time-varying GM(1,3) model of PM₁₀ concentration kernel and grey radius sequences, the simulated values of PM₁₀ concentration kernel and grey radius sequences can be respectively obtained. The SO₂ concentration and NO₂ concentration sequences that are predicted by the GM(1,1) model are used as the related factor sequences for 2015 to 2018. Then, a linear time-varying GM(1,3) model is used to obtain the predicted values of PM₁₀ concentration kernel and grey radius sequences.

Step 3: according to Equation (5), the simulated and predicted values of the upper and lower bounds over the PM₁₀ concentration interval grey number sequence are obtained by restoring the simulated and predicted values of kernel and grey radius sequences. The results are shown in

Table 4. The simulated and predicted values and the average relative error of the PM₁₀ concentration upper and lower bounds.

		Linear Time-Varying GM(1,3) Model			GM(1,3) Model
Year	Actual Value (μg/m3)	Simulated Value (μg/m3)	Lower Bound Relative Error (%)	Upper Bound Relative Error (%)	Simulated Value (μg/m3)	Lower Bound Relative Error (%)	Upper Bound Relative Error (%)
2008	[122,161]	[122.00,161.00]	0.00	0.00	[122.00,161.00]	0.00	0.00
2009	[121,148]	[107.05,130.05]	11.54	12.13	[105.05,126.89]	13.18	14.26
2010	[114,122]	[124.80,141.14]	9.47	15.69	[124.73,138.16]	9.41	13.25
2011	[114,121]	[121.66,123.27]	6.72	1.88	[115.38,125.57]	1.21	3.78
2012	[109,121]	[110.83,121.07]	1.68	0.06	[111.76,119.48]	2.53	1.26
2013	[108,114]	[110.18,115.51]	2.02	1.32	[110.10,114.59]	1.95	0.52
2014	[108,116]	[105.85,118.84]	1.99	2.45	[104.30,117.56]	3.43	1.35
Average simulated relative error (%)			4.78	4.79		4.53	4.92
Year	Actual Value (μg/m3)	Predicted Value (μg/m3)	Lower Bound Relative Error (%)	Upper Bound Relative Error (%)	Predicted Value (μg/m3)	Lower Bound Relative Error (%)	Upper Bound Relative Error (%)
2015	[102,116]	[103.77,109.13]	1.74	5.92	[103.88,106.75]	1.85	7.98
2016	[100,116]	[100.99,106.51]	0.99	8.18	[101.46,103.48]	1.46	10.80
2017	[95,108]	[98.52,104.11]	3.70	3.60	[99.14,100.56]	4.36	6.89
2018	[91,100]	[96.31,101.91]	5.84	1.91	[96.95,97.95]	6.53	2.05
Average predicted relative error (%)			3.07	4.90		3.55	6.93

.

Step 4: the relative error and average relative error of the PM₁₀ concentration upper and lower bounds are calculated according to Equation (6) and Equation (7). The results are shown in Table 4.

At the same time, the original GM(1,N) model based on interval grey number sequences is selected as the contrast model. Establishing the GM(1,3) model based on interval grey number sequences to simulate and forecast PM₁₀ concentration in Beijing from 2008 to 2018, and the relative error and average relative error of upper and lower bounds are calculated. The results are shown in Table 4.

This paper will compare the simulated and predicted values, relative errors of upper and lower bounds of the PM₁₀ concentration interval grey number sequence. According to Table 4, the average simulation relative error of upper and lower bounds of linear time-varying GM(1,3) model are 4.79% and 4.78%, respectively, and the average simulation relative error of upper and lower bounds of GM(1,3) model are 4.92% and 4.53%, respectively. The average simulation relative errors of those two models are very close, and both belong to high simulation accuracy. However, the average predicted relative error of upper and lower bounds of linear time-varying GM(1,3) model are 4.90% and 3.07% respectively, which are less than those of GM(1,3) model.

For ease of comparison, we make comparison maps as shown in Figure 2, Figure 3, Figure 4 and Figure 5. For simplicity, we briefly note the linear time-varying GM(1,N) model based on interval grey number sequences as LTVGM(1,N), and briefly note the GM(1,N) model based on interval grey number sequences as GM(1,N). From Figure 2 to Figure 5, the novel linear time-varying GM(1,3) model based on interval grey number sequences can more accurately describe the upper and lower bounds of PM₁₀ concentration in Beijing, which especially has the better prediction effect.

4.3. Forecast Results and Discussion

To understand the trend of haze in Beijing in the future, the linear time-varying GM(1,3) model proposed in this paper is used to forecast PM₁₀ concentration in Beijing for 2019 to 2021, the forecast results are shown in

Table 5. The predicted value of PM₁₀ concentration in Beijing for 2019 to 2021.

Year	2019	2020	2021
PM10 concentration (μg/m3)	[94.36,99.88]	[92.62,98.07]	[91.06,96.29]

. According to Table 5, PM₁₀ concentration in Beijing will decrease slowly from 2019 to 2021, but still exceed the China’s environmental air quality standard [13].

According to all the results considered, we will perform a discussion as follows: the average predicted relative errors of the linear time-varying GM(1,3) model are less than those of the original GM(1,3) model, which is attributed to the improved adaptability of the model to the dynamic change characteristics data. The model proposed in this paper is simple in the modeling method and convenient for calculation and applications. In addition, this novel model expands the range of predicted values from the real number to interval grey number, which can broaden the range of applications. To a certain extent, it can make up for the errors caused by data acquisition tools, acquisition conditions and acquisition personnel. However, this novel model still has certain limitations. It only considers that the related factor sequences are air pollutants, and it neglects the influence of meteorological factors on haze. Therefore, the future research will select both air pollutants and meteorological factors as the related factor sequences to study haze.

5. Conclusions

Aiming at forecasting haze of Beijing more accurately, this paper introduced a linear time polynomial into the GM(1,N) model based on interval grey number sequences, and established a novel linear time-varying GM(1,N) model based on interval grey number sequences. The data of PM₁₀ concentration, SO₂ concentration and NO₂ concentration in Beijing, China, were selected as modeling data to establish a linear time-varying GM(1,3) model from 2008 to 2018. The results indicate that the prediction model proposed in this paper has a higher prediction accuracy than the original model, and both the corresponding prediction errors are less than 5%, which proves the validity and practicability of this model. When using this novel model to make a forecast of PM₁₀ concentration in Beijing, the forecast shows that PM₁₀ concentration will be a downward trend for 2019 to 2021 in Beijing. This is because PM₁₀ concentration is not only determined by air pollutants, but also related to meteorological factors, such as wind speed, precipitation, and so forth. Specifically, when wind speed or precipitation increases, haze will decrease in Beijing [43]. Meteorological factors vary from year to year, and the extreme weather event such as El Niño cannot be accurately forecasted a year in advance. When the extreme weather event occurs, the actual PM₁₀ concentration in Beijing may be different with our forecast for 2019 to 2021. If there is no extreme weather event in that year, the haze predicted value of this paper will be accurate. Therefore, this model can provide decision support for the government when working toward greater haze control.

The model proposed in this paper had a high accuracy of prediction, but there is much room for improvement. This paper only selected two pollutants as related factors, which were SO₂ and NO₂, in order to forecast PM₁₀ concentration in Beijing. In future research, we will consider a few meteorological factors as related factors to construct model, such as wind speed, relative humidity, air pressure, temperature, precipitation, etc. Besides, we will apply this novel model to several key cities in the Beijing-Tianjin-Hebei region of China. On this basis, we will use data mining technology to compare and analyze the haze of different cities in the Beijing-Tianjin-Hebei region.