敖希琴， 郑 阳， 虞月芬， 汪金婷， 李 凡

(安徽新华学院，合肥 230088)

0 引 言

随着国家环境污染防治政策的出台，环境污染问题总体有所改善，但局部地区仍日趋恶化，多地爆发持续性雾霾天气，严重影响了居民身体健康。因此急需加大空气质量研究，特别是监控和预测大气污染。PM2.5[1]指的是每立方米空气中粒径小于或等于2.5 μm的颗粒物的含量，其值大小就表示了空气污染的严重程度，值越大污染就越严重。人们经常将其作为衡量空气质量优劣的重要指标之一。

当前预测空气污染的统计学方法主要有两种，分别为多元回归和时间序列模型等[2]。例如，沈劲[3]等将聚类分析与多元回归进行组合，建立了空气质量预报模型。张艺耀等[4]利用多元统计方法分析了PM2.5的影响因素，并建立了PM2.5的预测模型。喻彩丽[5]使用ARMA、GARCH时间序列算法对空气污染物进行预测分析。彭斯俊等[6]选用ARIMA模型建立了PM2.5预测模型。多元线性回归主要是从污染物的影响因素角度进行分析，而时间序列从污染物时间上的关联性角度进行分析，两种方法各有优势和不足。基于此，拟将二者结合起来，试图建立基于多元时间序列的PM2.5预测模型，并将其预测结果分别与多元线性回归和时间序列模型的预测结果进行比较，以比较3种方法在预测PM2.5预测时的精度。

1 理论基础

1.1 多元线性回归模型

多元线性回归分析是研究一个因变量与多个自变量之间线性关系的方法[7]，其形式一般为

Y=β0+β1X1+β2X2+β3X3+…+βjXj+…+βkXk+μ

(1)

式(1)中，β0为常数项，βj(j=1,2,,k)为回归系数，k为自变量的个数，μ为除去所有自变量对因变量影响后的随机误差。

1.2 时间序列模型

在统计研究中，一组按时间顺序排列的随机变量称之为时间序列，可以用一个数学模型来阐述或模拟该序列，当模型结构确定后，就可用该序列的历史数据来预测未来值。常用的随机时间序列模型是自回归移动平均模型(ARMA)，其结构可表示为

xt-φ1xt-1-…-φpxt-p=at-Φ1at-1-…-Φqat-q

(2)

式(2)中：{at}是残差序列；当q为0时，为模型AR(p)；当p为0时，为模型MA(q)；引入延迟算子B，当前时刻的序列值乘以一个延迟算子B，等于前一时刻的序列值，即xt-1=Bxt，xt-2=B2xt，…，xt-p=Bpxt。因此模型AR(p)模型、MA(q)模型、ARMA(p,q)模型可分别记为

φ(B)xt=at
xt=Φ(B)at
φ(B)xt=Φ(B)at

只有平稳的时间序列才能够直接建立ARMA模型，若是非平稳序列，需对数据作逐步差分，使序列满足平稳性的要求后再建模。

ARMA(p，q)模型建模的基本思路：首先可通过时序图和ADF检验对序列的平稳性进行验证，接着通过序列的自相关图和偏相关图进行模型识别，并确定模型阶数，然后根据得到的模型及阶数进行模型参数估计和模型检验，最后利用所建立的模型进行未来值预测。

1.3 多元时间序列模型

多元时间序列模型(ARMAX)是多元线性回归和ARMA模型的结合，称为传递函数模型动态回归模型[8]，其本质是指带回归项的ARMA模型，假设响应序列{yt}和输入变量序列(即自变量序列){x1t}，{x2t}，…，{xkt}均平稳，ARMAX模型的结构为

其中，Φi(B)为第i个自变量的自回归系数多项式；Θi(B)为第i个自变量的移动平均系数多项式；li为第i个自变量的延迟阶数；{εt}为回归残差序列。

由于{yt}和{x1t}，{x2t}，…，{xkt}均平稳，二者的线性组合也是平稳的，因此残差序列{εt}是平稳序列，其结构为

继续对{εt}序列中的相关信息进行取值，得到ARMAX模型为

其中，Φ(B)、Θ(B)分别是残差序列的自回归系数多项式和移动平均系数多项式；at是零均值白噪声序列。

2 数据来源及分析

2.1 数据的来源及预处理

通过“天气后报网”[9]采集了合肥地区2017-08-01到2018-07-30-的各污染物共2 184(364×6)组的监测数据。对收集到的数据整理时发现有少量缺失值和异常值。对于缺失值，采用简单插补法进行填充，即用缺失值相邻两个值的平均值进行替代。对于异常值，首先确定其位置并将其置零，然后按照缺失值处理方法对其进行处理。

2.2 相关性分析

相关性分析是用于判断变量与变量之间关系的密切程度。为考察PM2.5各影响因素间是否存在相关性，借助于R软件[10]，可得PM2.5和各个影响因素之间的散点图矩阵和相关系数表，分别如图1和表1所示。

图1 综合散点图Fig.1 The integrated scatter plot

表1 PM2.5与各影响因素间相关系数Table 1 Correlation coefficients between PM2.5 and influencing factors

由图1可知PM2.5浓度与PM10、SO2、NO2、CO、O3有一定的线性关系，与O3的关系不明显。由表1可得，PM2.5与PM10、CO呈现较强的正相关性；与SO2、NO2呈一般正相关；与O3呈负相关趋势。

3 多元回归模型的建立

3.1 模型的建立

根据相关性分析结果，课题将PM10、SO2、NO2、CO、O3作为影响因子引入模型，对PM2.5的浓度进行预测。建立多元线性回归方程，得到模型系数表如表2所示。

表2 系数表一Table 2 Coefficient table

根据系数表一可建立回归模型：

PM2.5=-20.17 1+0.42 7PM10-1.21SO2-
0.19NO2+73.401CO-0.092O3

3.2 模型的检验

分别采用了多重共线性检验、系数显著性检验、方程显著性检验、拟合优度检验对模型进行检验和评价，以达到验证模型有效性的目的。由系数的显著性检验结果可知除NO2外，剩余3个自变量均通过参数显著性检验。由多重共线性检验结果可得所有因素的方差膨胀因子(VIF)均小于10，表明自变量彼此相关程度不大，该多元线性回归方程是合理的。方程的显著性检验结果表明方程满足显著性要求。由拟合优度检验结果可知模型的R2为0.790 3，表明模型拟合优度较好。

3.3 模型的优化

由于在进行模型检验时，影响因子SO2和O3未通过系数的显著性检验，说明模型仍需优化，因此选用逐步回归分析法对模型进行筛选和优化。

逐步回归就是在多个自变量共同影响着一个因变量的关系中，从可供选择的所有变量中选出对因变量有显著影响的变量建立方程。课题中逐步回归分析结果如表3所示。

表3 系数表二Table 3 Coefficient table

根据系数表二可建立回归模型：

PM2.5=-18.825+0.402PM10-
1.607SO2+68.617CO-0.093O3

结合逐步回归分析结果，相较于原来的模型，新模型去除了NO2，为鉴别新模型是否能够应用于课题，同理采用拟合优度检验、方程显著性检验、参数显著性检验等指标验证模型，经过检验模型符合要求。并选用了AIC准则对两个模型的拟合优度进行了比较，结果表明经过逐步回归优化后的模型拟合度更好。

4 时间序列模型的建立

4.1 序列平稳性检验

根据研究需要并且结合数据的特征，选用了时间序列模型中ARMA模型进行建模，由于ARMA模型适用于平稳的时间序列，因此，建模前选用了单位根(ADF)检验对序列的平稳性进行检验。由检验结果可知PM2.5序列为平稳序列。

4.2 模型识别与定阶

对于平稳时间序列，模型的识别与定阶可通过自相关图和偏相关图的截尾和拖尾情况进行判定[10]。通过观察样本的自相关图和偏相关图，均呈拖尾趋势，分别如图2和图3所示。因此适用于ARMA(p,q)模型，通过图中滞后阶数的特征，预选模型ARMA(1,1)和ARMA(1,2)进行比较和评价。根据AIC、SC最小以及拟合度最优准则，最终选定了ARMA(1,2)作为预测模型。

图2 自相关图Fig.2 Autocorrelation graph

图3 偏相关图Fig.3 Partial correlation graph

4.3 参数估计以及模型的检验

对预测模型ARMA(1,2)模型进行参数估计，结果如表4所示：

表4 参数估计表Table 4 Parameter estimation table

进一步对预测模型的残差序列进行白噪声检验，结果表明残差序列的参数在5%的显著水平下，说明残差序列中信息提取充分，预测模型ARMA(1,2)建立合理。根据参数估计的结果，可得模型具体形式为

(1+0.701 1B)yt=(1+0.003 7B+0.065 1B2)at

5 多元时间序列模型的建立

5.1 协整关系检验

协整关系检验是用于确定输入序列与输出序列之间是否具有协整关系，可通过输入与输出之间的回归残差序列是否平稳进行判断[11]。若残差序列平稳，说明输入、输出序列存在长期均衡的关系，可建立多元时间序列模型。根据实验数据建立输入序列与输出序列的回归模型：

PM2.5=-22.735 2+0.396 8PM10-1.143 4SO2-
0.205 9NO2+76.195 1CO-0.066 6O3

对上述回归模型提取的残差进行平稳性检验，检验结果表明残差序列符合平稳性要求，可建立多元时间序列模型。

5.2 平稳性检验

分别对输入序列和输出序列进行平稳性检验，ADF检验结果表明输入序列{PM10}、{ SO2}、{NO2}、{CO }、{O3}序列为非平稳序列，输出序列{PM2.5}为平稳序列。对所有的非平稳序列进行差分，在经过一阶差分后，各输入序列均平稳。

5.3 各输入变量延迟阶数的确定

为了合理描述输入序列和输出序列的延迟关系，需计算输入序列和输出序列的延迟相关系数，找出延迟相关系数显著非零的延迟阶数，从而确定延迟效应。延迟相关系数[12]计算公式为

其中，cov(yt,xt-k)为延迟k阶的协方差函数，D(yt)为输入序列的方差，D(xt-k)为输出序列的方差，k为延迟阶数，即相关系数显著非零时的阶数，取值[0∶10]。通过计算，结果如表5所示。

表5 PM2.5与各影响因素间延迟相关系数Table 5 Delay correlation coefficient between PM2.5 and influencing factors

由表5可知，当k=0时，PM10、SO2、NO2、CO均显著非零，k=4时，O3显著非零，因此PM10、SO2、NO2、CO的取值从2017-08-01到2018-07-30，而O3取值从2017-07-28到2018-07-26。

5.4 ARMAX模型参数确定

首先对输入变量▽PM10(PM10的一阶差分序列)、▽SO2(SO2的一阶差分序列)、▽NO2(NO2的一阶差分序列)、▽CO(CO的一阶差分序列)、▽O3(O3的一阶差分序列)分别建立时间序列模型。使用相对最优法则定阶，确定▽PM10模型为ARMA(2,2)，▽SO2模型为ARIMA(1,1)，▽NO2模型为ARIMA(1,2)，▽CO模型为ARIMA(2,2)，▽O3模型为ARIMA(1,2)。另外各个输入变量的残差自相关检验结果表明各个拟合模型显著成立，具体形式如下：

▽PM10拟合模型的具体形式：

(1-0.305 2B+0.407 8B2)▽PM10=
(1-0.077 3B-0.820 7B2)aPM10

▽SO2拟合模型的具体形式：

(1+0.597 6 B)▽SO2=(1-0.954 6B)aSO2

▽NO2拟合模型的具体形式:

(1+0.479 9B)▽NO2=
(1-0.707B-0.209 7B2)aNO2

▽CO拟合模型的具体形式:

(1-0.354 4B+0.540 5B2)▽CO=
(1+0.009 7B-0.917 1B2)aCO

▽O3拟合模型的具体形式：

(1-0.345 8B)▽O3=(1-0.721 0B-0.170 3B2)aO3

根据延迟阶数并经实验可得如下结构模型：

拟合该ARMAX模型，依据条件最小二乘估计方法，模型可表示为

对回归残差序列{εt}进行白噪声检验，检验结果表明{εt}为白噪声序列，说明建模型有效。通过观察{εt}的自相关图、偏相关图，并根据AIC准则为残差构建ARMA模型。最终残差序列的模型为ARMA(2,2)。对参数进行最大似然估计后得到模型：

6 模型的验证

为了验证多元时间序列模型拟合效果是否优于单一的多元线性回归模型和时间序列模型，选用了均方根误差(RMSE)、平均绝对误差(MAE)和Theil不相等系数[13]3个指标来检验模型的预测效果。

其中yi表示真实值，gi表示预测值，ARMSE值和AMAE值都是越小，表示预测值与真实值越接近，预测准确度越高；U取值在0到1之间，U越接近0，模型预测越准确。

根据已建立的回归模型、时间序列模型、多元时间序列模型，预测合肥市2018-07-31至2018-08-13合肥市的PM2.5浓度。各评价指标结果如表6所示。

表6 模型指标评价结果Table 6 Model index evaluation results

由表6可知，3个评价指标均表明ARMAX模型的预测精度最高，模型拟合度最好，多元线性回归预测效果次之，预测精度最低的为时间序列模型。

7 总 结

以合肥市2017年和2018年PM2.5和其他污染物相关数据为分析基础，构建了用于预测PM2.5浓度的多元时间序列模型(ARMAX)，相较于多元线性回归模型和时间序列模型，模型预测精度更高。首先，通过散点图矩阵和相关系数分析了各影响因素与PM2.5之间的相关性，在此基础上，分别构建了多元线性回归模型和时间序列模型。继而对输入、输出序列进行了协整关系以及平稳性检验，并根据延迟相关系数确定了各影响因素延迟阶数，建立了各影响因素的时间序列模型，之后将时间序列模型与多元线性回归模型进行整合，得到最终的多元时间序列模型。最后选用RMSE、MAE和Theil不相等系数对3个模型的预测精度进行评价，得到多元时间模型的预测精度明显优于单一的多元线性回归模型和时间序列模型。

虽然提出的ARMAX模型可以在一定范围内提高PM2.5浓度的预测能力，但依然有以下不足之处：

(1) 由于实验数据中存在少量缺失值和异常值，数据波动性较大，会导致建立的模型预测效果不理想，而课题在建模前仅选用了简单插补法进行处理，在后期的研究中是否可以采用一些其他更为科学合理的方法进行缺失值和异常值的处理，从而提高模型的预测精度。

(2) ARMAX模型会随着预测期的延长，而产生预测精度下降的情况，是否可以建立一个动态的PM2.5预测模型，即当预测精度超出所设定的精度时，模型就会自动进行修正，以提高预测的精准性，这方面的工作值得进一步研究。