张学新， 周聿瑾

(1.湖北工程学院数学与统计学院， 湖北 孝感 432000； 2.武汉英中高级中学， 武汉 430022)

引 言

城市空气污染治理是一个复杂的、系统性的工作，需要联动周边城市一起治理。对于大都市而言，它的空气污染除了自身的因素之外，也与周边城市的空气污染密切相关，因此需要探索大都市及其周边城市空气污染物之间的影响机理，不能孤立地分析一个城市的空气污染问题。

关于城市空气污染物的影响分析，较多的是研究空气污染物与气象因素、地理因素、人为因素等非空气污染物因素的关系。李琛等[1]以汉中市城区2015年空气质量监测数据和气象资料为基础，分析温度、露点温度、平均海平面大气压等气象因素对空气质量的影响及其空气污染特征。Tao等[2]利用长江三角洲的350个城市的有关数据研究城市形态对城市中PM2.5浓度的影响，城市越紧凑，其PM2.5浓度越高。安玉琴等[3]探讨2015-2016年石家庄市、张家口市空气污染物与气象要素的相关性，认为温度、相对湿度、风速为影响污染物浓度的主要气象要素。王莉等[4]采用Spearman相关系数分析泰州市2016年度各种空气污染物和气象因素间的关系，结论是不同气象因素对空气污染物的浓度变化有较大影响，PM2.5与其他常规空气污染物间存在一定相关性。孟昭伟等[5]应用多元回归分析西安市莲湖区和雁塔区PM2.5质量浓度的变化特征及其与气象条件的关系。Xiao等[6]应用多元线性回归揭示京津冀及其周边地区气象因素(尤其是风和降水)与空气质量之间关系。结果表明，从2015年到2018年，京津冀及周边地区的空气质量表现出明显的季节特征，气象因素与空气质量负相关。李惠娟等[7]根据环保部门的污染源解析结果，建立面板计量经济学模型探讨空气污染人为影响因素，结果表明煤炭消耗、汽车拥有量是城市空气污染的显著性影响因素。高正旭等[8]利用武汉市2014-2017年大气污染物和气象要素的观测数据，使用HYSPLIT、PSCF和CWT分析方法揭示武汉市不同大气污染物的潜在源区分布及其贡献特性，认为武汉市本地及其周边地区污染对武汉市大气污染物的贡献较大。陈静[9]基于GIS的变异性分析，采用普通克里金算法对江浙沪16个城市的PM2.5月均浓度数据插值，分析了该地区PM2.5月均浓度的空间变异性，分析了人为源VOCS排放污染源与当地的风向、风速和降雨量对空气污染指数空间分布的影响。白杨等[10]利用空间插值、趋势分析和相关性分析方法，研究了西安市一次重空气污染期间 6 大污染物的质量浓度时空变化及彼此间的相关关系，未涉及西安周边城市的大污染物。何英华等[11]收集A市、B市连续三年每日的空气主要污染物环保监测数据，通过描述性统计对比分析A、B两市各种空气污染物在时间、空间上的分布规律。

关于一个大城市及其周边地区的空气污染物相互影响与治理质量的深度研究比较少见。张学新等[12]利用武汉及其周边城市空气质量指数(AQI)及相关污染物数据，探索了武汉市空气质量指数的统计分布规律及空气污染治理效果的评价问题，并用统计模型分析武汉及周边大中城市的SO2等空气污染物之间的传播及相互影响机理。总体看,对某个城市空气污染物影响分析，鲜有从空间关联上、联系地、整体地研究城市群之间的空气质量及影响。这里把北京及周边大中城市作为一个系统，应用向量自回归理论整体地探讨各地实体污染物间的相互影响、应用变点检测技术实证北京市空气质量治理效果。

1 污染物平稳特性

数据来源及记号。从周边八个方位看，北京有廊坊、天津，保定、石家庄、大同、呼和浩特，唐山、秦皇岛，张家口、乌兰察布、承德、沈阳共12座大中城市。对每一个城市，采集2013-12-02到2019-01-29期间的所有空气质量指数日历史数据(仅缺乌兰察布2014年的数据)。

数据源自《空气质量历史数据查询》(https://www.aqistudy.cn/historydata/index.php)、《天气后报》(http://www.tianqihoubao.com/aqi/)。对少量缺失记录，通过2-邻近值方法插补。以CO_BJ表示北京CO指标序列，CO_BJ(-1)表示北京CO指标序列的滞后1天序列,其他同理。

经初步分析，北京及周边城市空气污染物的统计分布很复杂，需要把它们作为一个向量系统深刻分析其相互关系。其中取Epanechnikov核函数K(u)=0.75(1-u2)I(|u|≤1)时，各城市CO污染物的核密度估计如图1所示,它们不属于任何常见的概率密度曲线。

图1 北京及周边城市CO核密度曲线

相关系数分析显示，北京的CO与呼和浩特、石家庄、天津、保定、承德当天的CO及滞后一天的CO(-1)，以及与沈阳、廊坊、唐山、秦皇岛、张家口当天的CO有较强的关联(相关系数均大于0.5)。北京的CO与大同、乌兰察布CO的相关性不明显。为避免伪相关，还需检验两个独立序列的平稳性。假设污染物序列{PCt}满足：

(1)

基于SIC准则、取最大滞后长度为25、视外生变量为常数，对北京、天津、石家庄、呼和浩特、沈阳、廊坊、保定、唐山、张家口、承德的CO序列CO_BJ、CO_TJ、CO_SJ、CO_HT 、CO_SY、CO_LF、CO_BD、CO_TS、CO_ZJ、CO_CD进行ADF检验，结果见表1。在0.05的显著性水平下，分别拒绝以上序列“有一个单位根(α=0)”的零假设，表明它们是平稳序列。

表1 CO序列的ADF检验

同理，北京及周边城市的SO2序列、NO2序列、PM2.5序列、PM10序列都是平稳的。对于臭氧8小时(O3_8h)序列，除了北京西面的大同、呼和浩特是一阶差分平稳序列外，其余都是平稳序列。

2 北京污染物的影响因子

使用Granger因果检验可以判断一个城市的空气污染是否由另一个城市的渗透漂移过来。把一个城市的某种污染物浓度PC1的历史信息加入到另一个城市同名污染物PC2的预测模型，如果预测能力得到显著性提高，就说PC1是PC2的(单方面)Granger原因，如果PC2也是PC1的(单方面)Granger原因，则称它们互为Granger因果关系。对模型

(2)

设观测值个数n，SSE1、SSE0分别是模型及H0下模型的残差平方和，则模型的统计量为H0:βj=0(j=0,1,…l)，其整体显著性检验为：

F=m-1(n-l-m-1)(SSE0/SSE1-1)～F(m,n-l-m-1)

(3)

在0.05的显著性水平下，各污染物序列的成对Granger因果检验结果如下：

考察l(l≥3)期滞后，对于PM10、SO2、O3_8h污染物，北京与周边城市均互为Granger因果关系。对于PM2.5，大同、乌兰察布、沈阳是北京单方面Granger原因，北京是承德单方面Granger原因，北京与其他周边城市均互为Granger因果关系。对于CO，张家口、沈阳是北京的单方面Granger原因，北京与周边其他城市均互为Granger因果关系。对于NO2，张家口、沈阳是北京单方面Granger原因，北京是唐山单方面Granger原因，北京与其他周边城市均互为Granger因果关系。

考察2期滞后，对于PM2.5、O3_8h污染物，北京与周边城市均互为Granger因果关系。对于PM10，呼和浩特是北京单方面Granger原因，北京是承德单方面Granger原因，张家口、承德、沈阳是北京单方面Granger原因，北京与其他互为Granger原因。对于SO2，张家口、沈阳是北京单方面Granger原因，北京与其他互为Granger原因。对于CO，廊坊、沈阳是北京单方面Granger原因，北京与其他互为Granger原因。对于NO2，张家口、承德、沈阳是北京的单方面Granger原因，北京与其他互为Granger原因。

总体来看，沈阳单方面污染了北京，其他三个大城市与北京相互污染。在2天内，主要是周边中等城市污染了北京。

分析发现北京PM10与北京等地SO2、PM2.5互为Granger因果关系。例如，在任意滞后期上，北京PM2.5与北京PM10均互为Granger因果关系，1～2期的结果见表2。

表2 北京PM2.5与PM10的成对Granger因果检验

下面结合Granger因果关系，给出各污染物的回归方程。

对PM2.5分析结果见表3,其回归方程估计式为：

表3 北京及周边城市PM2.5回归估计摘要

PM2.5_BJ=2.99+0.349PM2.5_BJ(-1)+0.536PM2.5_LF-0.083PM2.5_TJ-0.077PM2.5_BD(-1)-0.121PM2.5_ZJ(-1)+0.954PM2.5_CD-0.533PM2.5_CD(-1)+0.029PM10_BJ(-1)

(4)

R-squared=0.850547，Adjusted R-squared=0.849909，Durbin-Watson stat=1.981030，Prob(F-statistic)=0.0000

对误差项νt的ADF检验结果是t-Statistic=-18.181 35，Prob.=0.0000，再对其作Breush-Godfrey LM序列相关性检验，结论是2期内不存在自相关，但有异方差。这些表明虽然残差序列νt的分布比较复杂，但不违背模型的假定，它是平稳随机项。

这表明北京PM2.5主要受到北京前一天PM2.5及承德、廊坊当天PM2.5、承德、张家口前一天PM2.5影响。同理，有：

PM10_BJ=-6.206+0.109PM10_BJ(-1)+0.619PM10_LF+0.138PM10_TJ+0.125PM10_DT(-1)+0.039PM10_QH(-1)+0.911PM10_ZJ-0.093PM10_ZJ(-1)-0.041PM10_SY(-2)+0.869PM2.5_BJ-0.208PM2.5_BJ(-1)-0.750PM2.5_LF+0.074PM2.5_TJ(-2)-0.321PM2.5_CD+0.575SO2_BJ-0.148SO2_TJ(-2)-0.097SO2_DT(-1)-0.111SO2_TS+0.294SO2_CD(-1)

(5)

方程(5)的R-squared=0.864 170，模型整体显著性Prob(F-statistic)= 0.000 000，系数显著性检验的Prob.值除一个0.0177、一个0.0006、一个0.0003外，余下部分几乎全是0.0000，Durbin-Watson stat=1.945 032。这表明北京PM10主要受到北京前一天的PM10及廊坊、张家口、天津当天的PM10影响，受北京、廊坊、承德的PM2.5影响，受北京、承德、唐山的SO2影响。

SO2_BJ=-2.99+0.407SO2_BJ(-1)+0.268SO2_LF-0.063SO2_LF(-1)-0.070SO2_TJ+0.046SO2_SJ+0.056SO2_TS-0.015SO2_TS(-1)+0.066SO2_ZJ+0.219SO2_CD-0.171SO2_CD(-1)+0.022PM2.5_BJ-0.965CO_BJ(-1)

(6)

方程(6)的R-squared=0.884 466，模型整体显著性Prob(F-statistic)= 0.000 000，系数显著性检验的Prob.值除一个0.0151外，其余全是0.0000，Durbin-Watson stat= 1.944 419。它表明北京市的SO2污染物主要受北京前一天的SO2、廊坊、承德、张家口、唐山当天的SO2、承德前一天的SO2及北京前一天的CO影响。

CO_BJ=0.458CO_BJ(-1)+0.401CO_LF-0.086CO_LF(-1)-0.087CO_TJ+0.115CO_TJ(-1)+0.048CO_BD(-1)-0.032CO_BD(-2)-0.046CO_DT+0.147CO_SJ-0.107CO_SJ(-1)+0.053CO_TS-0.034CO_TS(-1)+0.258CO_ZJ-0.147CO_ZJ(-1)+0.426CO_CD-0.345CO_CD(-1)-0.061CO_SY(-1)

(7)

方程(7)的R-squared=0.846 731，模型整体显著性Prob(F-statistic)=0.000 000，系数显著性检验的Prob.值除一个0.0044、一个0.0020、0.0012外，其余全是0.0000， Durbin-Watson stat= 1.960 725。它表明北京的CO污染物主要受北京前一天的CO、承德、廊坊、张家口、石家庄当天的CO、承德、张家口、石家庄、天津前一天的CO影响。

NO2_BJ=-2.87+0.495NO2_BJ(-1)+0.382NO2_LF-0.279NO2_LF(-1)+0.040NO2_TJ(-1)+0.028NO2_SJ+0.123NO2_QH+0.162NO2_ZJ+0.697NO2_CD-0.445NO2_CD(-1)

(8)

方程(8)的R-squared=0.879 152，系数显著性检验的Prob.值除一个0.0197、一个0.0153外，其余全是0.0000，Durbin-Watson stat=2.072 212。它表明北京市NO2污染物主要受北京前一天NO2、承德、廊坊当天及它们前一天的NO2影响。

O3_8h_BJ=-5.01+0.243O3_8h_BJ(-1)+0.520O3_8h_LF-0.151O3_8h_LF(-1)-0.107O3_8h_TJ+0.171O3_8h_BD-0.099O3_8h_SJ(-1)+0.082O3_8h_DT+0.115O3_8h_HT+0.062O3_8h_TS+0.104O3_8h_TS(-1)-0.082O3_8h_Qh(-1)-0.061O3_8h_WL+0.374O3_8h_CD-0.114O3_8h_CD(-1)

(9)

方程(9)的R-squared=0.938 044，系数显著性检验的Prob.值除一个0.0038、一个0.0016外，其余是0.0000，Durbin-Watson stat= 2.049 539。它表明北京市的O3_8h污染物主要受北京前一天的O3_8h、廊坊、承德、保定、呼和浩特当天的O3_8h、廊坊、承德前一天的O3_8h影响。

以上Granger因果检验方程中有些系数估计值是负的，可理解为互为Granger原因的体现，一方的污染加重(减轻)是因为吸入(排放给)另一方的污染物，导致另一方的污染物浓度降低(增加)。虽然可以做很多额外的努力去修正模型，比如改善因地理空间相邻而出现的复共线现象，但是用修正模型解释北京周边城市污染物对北京污染物的影响将更晦涩。总体来看，各种回归模型良好。

3 协整检验

取α=0.05，对北京PM2.5与其周边城市PM2.5的协整检验均不显著，对其他PM2.5、PM10、SO2、CO、NO2、O3-8h污染物的协整检验也没有通过，表明在目前它们已经没有明显的稳定关系，北京周边城市的大气污染防治工作均取得一些效果。

4 影响因子的异动性

为分析北京及周边城市空气污染物的变动情况，了解哪些回归因子的作用发生了显著性变化，需要对前述线性模型作稳定性诊断。这里提出一个线性模型多个未知变点的检验方法。设有m个变点的多元线性回归模型：

(10)

其中：T0=1，Tm+1=T，xt∈p与zt∈q是回归因子向量，x变量的系数不在体制内变化，而z变量的系数是体制变化决定的，误差项ut～i·i·d。

对任意一个m-分割{Tj}=(T1···,Tm)，Tj=[λjT],0<λ1<···<λj<···<λm<1，最小二乘法通过

(11)

(12)

对预先指定的变点个数k>0，检验假设是：

H0:m=0，H1:m=k

(13)

supFT(k;q)=sup(λ1,···,λk)FT(λ1,···,λk;q)

(14)

其极限分布依赖于回归因子以及残差序列相关性和异方差性，对于各种修剪参数Trimming，回归因子个数和变点个数，有相应的临界值及响应曲面的计算，当supFT(k;q)超过临界值时，拒绝原假设，认为存在k个变点。

前述分析知道，北京前一天的CO，廊坊、天津、保定、石家庄、唐山、张家口、承德当天或前一天的CO对北京当天的CO影响较大。现考虑变量CO_BJ(-1)、CO_LF(-1)、CO_TJ、 CO_TJ(-1)、CO_BD(-1)、CO_SJ、CO_TS、CO_ZJ、CO_CD是否因体制变化而变化，做多变点检验。有关参数设置及检验结果见表4。这里体制变化是指北京市历次空气重污染应急预案的发布施行。对比北京市最近几次发布《北京市空气重污染应急预案》修订版的时间，考虑政策实施前的试行及政策起作用的滞后效应，可见估计的变点日期比较符合实际。

表4 北京CO浓度分布因应急预案发布变异性检验

设置Trimming=0.15,最大变点数4, 显著性水平α=0.05。

模型稳定性诊断表明北京、廊坊、天津、保定、石家庄、唐山、张家口、承德等CO对北京CO的污染发生了变化，影响力度实际变小。为更简单一些，把在2017年9月15日至2019年1月29日的子样本上拟合的模型与在全段上的模型比较，北京周边城市CO的影响明显减弱。CO_BJ(-1)、CO_LF(-1)、CO_TJ、CO_TJ(-1)、CO_BD(-1)、CO_BD(-2)、CO_SJ、CO_TS、CO_ZJ、CO_CD的估计系数显著下降，此外对CO_BD(-1)、CO_DT、CO_SJ(-1)、CO_TS(-1)、CO_ZJ、CO_ZJ(-1)、CO_SY(-1)的系数已经没有统计上的显著性，不宜再保留在原模型中。这说明，经过多方治理，保定、大同、张家口当天的CO、北京、石家庄、唐山、张家口、沈阳前一天的CO对北京当天CO的污染显著性降低。

5 结束语

北京市各实体空气污染物浓度的统计分布复杂，不宜使用回归分析与时间序列分析建模处理。当缺乏相关理论动态分析两个变量之间的作用方向时，Granger检验提供了参考。它是通过考察在模型中加入一个变量的历史信息能否提升模型的预测能力来指明两个变量之间的先后影响方向。

把Granger检验应用于北京及其周边城市的两个污染物浓度序列，发现北京各空气污染物的主要影响是：承德、廊坊、张家口三地的PM2.5浓度；廊坊、张家口两地的PM10浓度及承德的SO2浓度；廊坊、承德、张家口三地的SO2浓度；承德、廊坊两地的CO浓度；承德、廊坊、张家口三地的NO2浓度；廊坊、承德、保定、呼和浩特四地的O3_8h浓度。北京当天污染物也受其前一天污染物的持续影响。

基于当前可用数据(最新数据)的协整检验没有发现北京与各周边城市污染物两者之间存在协整关系，表明北京与周边城市污染物之间过去高度相关的关系已经消逝，说明北京与各周边城市两个污染物序列之间不具有长期稳定的关系，证实北京及周边城市大气污染防治工作系统治理取得了一定成效。进一步通过结构多变点检验，实证了2012年至2017年间北京实施空气重污染应急预案治理空气质量效果显著。这种动态监控一个城市与其周边地区污染物流动的研究方法，给区域城市群协同、系统地治理大气污染提供了精准治理方向。