近年来，中国的空气污染问题已成为当前面临的严峻的环境问题之一.空间统计学作为近年来逐渐热门的研究领域，在对空气污染物空间分布特征的研究中得到了越来越多的应用.空间自相关分析是空间统计学的重要内容之一，空间自相关性来源于空间地理学的重要假定，即“地理要素之间的距离越近，其属性值也就越接近”.常用的空间自相关的衡量方法有Moran′s I指数、Geary′s C指标等(Schabenberger et al., 2005).对空气污染的空间自相关特征进行研究，有助于掌握空气污染物的空间分布规律.许多学者对空气污染物浓度或者空气污染指数的空间自相关性进行了研究，例如，韩贵锋等(2006)对重庆市空气污染物的空间分布规律进行了研究，发现不同污染物有着不同的空间自相关特征；昌晶亮等(2015)对珠三角地区空气中PM_2.5浓度的空间分布规律进行了研究，探讨了其空间自相关特征的尺度效应；胡芳芳(2010)基于空气污染指数(API)分析了北京市及其周边区县空气污染程度的空间聚集特征，发现存在明显正的空间自相关.此外，Xu等(2016)对中国的京津冀、长三角和珠三角地区的PM_2.5浓度进行了空间自相关分析，也发现了区域尺度上的显著的空间自相关效应.

为了更加全面地评估中国的空气质量，从2013年之后，国家环保部停用了空气污染指数(API)，而开始使用新的空气质量指数(AQI).AQI是比API更为全面的衡量空气质量好坏的无量纲指标，其计算过程是首先计算单项污染物的分指数，然后取最大值计算AQI.目前已有多个对空气质量的空间分布特征的研究是基于AQI的(王冠岚等, 2016；柳笑盈, 2015).对于AQI而言，空间上接近的城市有着较为相近的AQI值也是正常的.首先，空间上接近的城市其经济发展水平和空气污染物排放水平一般也更加接近，因而导致邻近城市的AQI可能有着较为接近的均值水平.此外，空间上接近的城市更容易受到类似气象条件的影响，这也会导致其空气质量有着相类似的短期波动规律，即存在着空间自相关.

以往对于空气质量的空间自相关分析(昌晶亮等, 2015；韩贵锋等, 2006)，大都针对局部地区，少见关于中国整个东部地区城市群的分析.需要强调的是，全局空间自相关分析(如Moran′s I指数)与时间序列自相关分析类似，需要假定数据均值和方差的平稳性，也就是说计算Moran′s I指数需要假定空间变量的均值和方差是相同的(Schabenberger et al., 2005).空间上存在趋势项会干扰空间自相关分析，因为邻近空间点位样本值接近的原因可能不是由于空间自相关，而仅仅是因为它们来自于均值接近的相互独立的分布(Schabenberger et al., 2005).空间维度上存在低频变化趋势可以增大Moran′s I的值(Waller et al., 2004).在较大区域应用空间自相关分析时，空间平稳的假定很难被认为是满足的.这是因为存在气候条件和污染源分布的差异，中国不同区域城市的空气污染程度呈现出明显的区域性差异(段玉森等, 2008), 这种在空间上存在的均值非平稳性显然会影响对空间自相关性的判断(Schabenberger et al., 2005).

对于只有一次观测的空间变量，建立均值模型较为困难.在这种情况下，如果邻近位置的变量值类似，很难区分这是由空间趋势还是由空间自相关所导致的.但对于日AQI数据来说，由于其同时具有时间和空间的维度，因此，本文提出从时间维度来对均值过程进行建模，分别对每个城市的AQI序列建立单变量时间序列模型，进而以模型的残差项为基础进行空间自相关分析.也就是说，首先从时域上对AQI数据进行预置白(pre-whitening)，这样获得的残差序列在时域上就是独立的白噪声序列.从空间上看，残差的均值是相等的(都为0).这样，对残差的空间自相关检验就排除了在AQI原始数据中存在的空间趋势的干扰.此外，依据时序模型残差进行空间自相关分析有着较为明确的意义，叙述如下.

第一，时间序列模型的残差也就是模型在样本内的一步预测误差.时序模型的本质是依据时间序列的“惯性”(即时域自相关性)，根据历史数据来预测未来值.不考虑时序模型的具体形式，对于时间序列Y₁, Y₂, …, Y_t, …，时序模型可以表达为式(1), 其中, e_t为t时刻的随机项.而残差就是实际值与预测值之差，即式(2).

(1)

(2)

因此，t时刻的残差实际上就是无法被历史数据Y_t－1, …, Y₁解释的那一部分信息.在某种意义上，对于AQI预测(无论是基于统计模型还是数值模型)来说，对时序模型残差的解释能力应该是评估模型的标准.也就是说，任何有效的AQI预测模型应该能进一步解释历史AQI数据无法解释的那一部分信息.残差中存在的空间自相关特征可能对预测模型提供指导信息.

第二，大量研究表明，影响空气质量的因素可以归结为两类，即污染物排放量和气象因素(王冠岚等, 2016；王红丽, 2015；张国琏等, 2010；谈建国等, 2007).对每个城市来说，污染源的排放在相邻2 d之内一般不会有大的变化，所以在相邻2 d，AQI偏离其时序“惯性”的变化(即时序模型没有解释的残差项)则主要来自于有利于(或者不利于)污染物累积的气象条件(如风向、风速、降水等)及气象条件导致的跨区域传输的影响.这样，不同城市残差序列的空间自相关结构实际上就反映了邻近城市空气质量在气象条件驱动下的同步响应特征.

基于此，本文以AQI时序模型残差作为空间统计分析对象，探讨其空间自相关特征，以期更深刻地了解中国东部城市群空气污染的空间分布规律.

本研究所使用的AQI数据来自国家环保部网站的全国城市空气质量日报，共包括366个城市的数据.首先剔除了9个缺测天数在10 d以上的城市；同时，考虑到站点的分布密度，本文只分析100°E东部的城市.这样，本研究一共使用了327个城市的数据，具体如图 1所示.在这些城市中，有156个城市的AQI指数记录从2014年1月1日开始，另有169个城市的数据从2014年12月31日或者2015年1月1日开始有记录；此外，还有2个城市的数据记录时间从2015年1月23日开始.在下文中，将观测值开始于2014年1月1日的156个城市称作“第一批次城市”.本文使用的数据截至日期为2017年1月1日.

图 1(Fig. 1)

图 1 本文使用的AQI监测城市分布(实心点表示从2014年1月1日即有记录，空心点表示记录开始时间为2014年12月31日之后) Fig. 1 Locations of cities used for AQI analysis in this study

全局空间自相关最常用的检验方法为Moran′s I指数检验(张松林等, 2007).Moran′s I指数的定义类似于皮尔逊相关系数，其计算公式见式(3)(Bivand et al., 2008；Negreiros et al., 2009).

(3)

式中，x为样本的平均值，w_ij为第i个和第j个样本点的距离权重(即权重矩阵W第i行第j列的元素)，n为样本点的个数，S₀为所有权重之和，即S₀.如果对于所有的样本点，其本身数值大小与其邻域样本点数值大小都呈现出“高-高”或者“低-低”的分布类型，则式(3) 分子的乘积项为正值，所以I值也是正值，即说明空间变量之间呈现正相关性.而如果样本点其本身数值大小与其邻域样本点数值大小都呈现出“高-低”或者“低-高”的分布类型，则I值为负值，说明空间变量之间呈现负相关性.与皮尔逊相关系数相同，Moran′s I指数的值介于-1~1之间.

在用式(3) 计算Moran′s I指数时，需要确定权重矩阵W.在计算w_ij时首先确定第i个样本点的邻域，在这里将邻域半径定义为500 km.如果第j个样本点落在邻域中，则使用其与第i个样本点距离的反比作为w_ij的值；如果第j个样本点落在邻域之外，则w_ij值为0.最后，将每个样本点对应的所有近邻点的权重之和标准化为1，即权重矩阵W的每行之和标准化为1.

在对Moran′s I指数进行显著性检验时，首先需要将Moran′s I指数进行标准化(式(4)).

(4)

式中，E(I)是Moran′s I指数在原假设(即不存在空间自相关)条件下的数学期望，而SD(I)是标准差.在无空间自相关的原假设条件下，局域Moran′s I指数的数学期望和方差的计算方法可以参考相关文献(Anselin, 1995).当Z(I)过大或者过小时可以拒绝原假设，而单边假设检验的p值为p=P(Z＞Z_观测).后文使用的Moran′s I指数都是标准化的Z(I).本文进行的Moran′s I检验基于R语言扩展包spdep.

全局Moran′s I指数从整体上描述了样本点之间自相关的大小，但这可能会掩盖空间自相关特征的空间异质性，即不同的区域空间自相关特征是不同的.也就是说，全局Moran′s I指数只能检验空间自相关性，不能对空间自相关进行定位.为了描述局部的自相关特征，Anselin(1995)提出了局域Moran′s I指数，如式(5) 所示.

(5)

对某个样本点，其局域Moran′s I指数反映了空间变量在此位置附近的自相关性质.某一点的局域Moran′s I指数高，说明这一点与周围的点倾向于同方向变动.所有的I_i相加之和正比于全局Moran′s I指数.与全局Moran′s I指数一样，后文所使用的局域Moran′s I指数值也是标准化之后的Z(I)值.

经典的时间序列建模流程来自于Box和Jenkins的分析方法(Box et al., 1976).建模的一般过程可以归结为模型识别、参数估计和模型诊断(Cryer et al., 2008).时间序列的自相关函数可以用来判断时间序列的平稳性，对于时间序列Y₁, Y₂, ..., Y_n，其时滞为k的样本的自相关系数见式(6).

(6)

如果自相关函数r_k随着时滞k的增加迅速衰减到0附近，则表明序列可以被看作是平稳的.平稳时间序列可以用自回归滑动平均模型(ARMA)来建模，ARMA模型的形式见式(7).

(7)

式中，e_t是白噪声序列，p和q是ARMA模型的阶数，φ₁, φ₂, …, φ_p, θ₁, θ₂, …, θ_q是待估计的参数，可以用极大似然法来对其进行估计.

ARMA模型只能用来对平稳时间序列进行建模，如果时滞k很大而样本自相关函数r_k仍然不衰减到0，则表明序列不平稳.对于非平稳时间序列，可以通过差分来使时间序列平稳化，进而对差分序列建立平稳时间序列模型，这种建模方式被称为自回归滑动平均求和模型(ARIMA).

尽管可以通过自相关函数的形态来判断ARIMA模型的阶数，但更加通用的方法是通过最小化来实现：

(8)

式中，L是极大化的似然函数值，k=p+q+1.这被称为赤池信息准则，其原理是惩罚阶数较大的模型，来确保选择尽可能简单的模型.

通过以上的ARIMA模型，可以获得在空间上均值平稳的残差项.但计算Moran′s I指数需要变量在空间上的方差也是平稳的，而由于不同城市的残差项波动程度不相同，所以需要对其进行标准化.在分析每个城市的残差序列时，发现普遍存在波动集群效应，即存在着波动和平静交替的模式，这很大程度上是受到不同季节气象条件的影响.本文使用广义自回归条件异方差模型(GARCH)来对这种方差前后的相依性进行建模，其定义见式(9)(Bollerslevb, 1986).

(9)

式中，e_t为ARIMA模型的残差，ε_t为具有零均值和单位方差的独立同分布随机变量序列，e_t的条件方差为σ_t|t-1².式(9) 表明，条件方差具有记忆性，即大的波动会跟随大的波动.通过条件方差的估计值，获得标准残差，具体见式(10).由于标准残差的方差都是1，所以在空间上标准残差是方差平稳的.本文以标准残差序列为基础分析其空间自相关特征.

(10)

本文在确定ARIMA阶数时，使用了R语言forecast包中的auto.arima函数.参数估计时，首先对原始AQI序列进行差分处理，然后基于R语言rugarch包对ARMA-GARCH过程进行参数估计.GARCH模型的阶数取为GARCH(1, 1).对于时间序列模型的自相关分析、定阶、参数估计等方法可以参考相关文献(Cryer et al., 2008)，本文不再赘述.

在分析城市逐日AQI序列时，需要注意的是其时间非平稳性.非平稳的一个最重要的原因是不同的季节其气象特征不同，污染物扩散水平也不相同，导致空气污染水平呈现出大致的以年为周期的变化特征.此外，在长期时间范围内，污染物的排放量也会有变化趋势，所以日AQI序列不是平稳的时间序列.

对城市群AQI序列的时间自相关分析如图 2所示.图 2a给出了所有城市原始AQI序列从滞后1 d到滞后4 d的自相关函数值的分布情况.从图 2a可以看出，对于大部分城市的AQI序列来说，其一阶自相关函数都大于0.5，最高者可以大于0.8.大多数城市的自相关函数随着阶数的增加衰减缓慢(图略)，表明原始的AQI序列是非平稳的，这与上面的分析是一致的.

图 2(Fig. 2)

图 2 中国城市群AQI指数序列(a)、差分序列(b)及时序模型标准残差序列(c)的自相关系数的频次分布特征 Fig. 2 Frequency distribution of autocorrelation coefficients of the original daily AQI series(a), the difference series(b) and the residuals series from time series models(c)

对日AQI序列差分之后，再次计算自相关函数(图 2b)，发现当滞后达到2~3阶时，自相关函数已经很接近于0，表明一阶差分序列是平稳的.从图 2b中可见，对于差分序列，前几阶的自相关函数都是负值，表明从统计上来说，日AQI序列在增高后的几天内略有回落的倾向.根据上面的结果，本文对每个城市的AQI序列分别建立ARIMA(p, d, q)模型，其中，d设置为1(即一阶差分).使用2.3节叙述的时间序列建模方法，获得每个城市AQI序列的残差序列和标准残差序列.本文计算了标准残差序列的自相关函数，结果见图 2c.可见时间序列模型已经充分提取了序列的自相关信息，标准残差的时滞自相关都在0附近.

为了看出残差反映信息的物理意义，选取北京市2014年10月1—15日这个时间窗口来进行分析.在这段时间内，出现了一次典型的污染过程(从7日持续到11日)，此次污染过程已有研究者进行了详尽的讨论(王占山等，2015).图 3a显示了这个时间窗口的AQI观测值和时序模型的拟合值(即样本内一步预测值).从图 3a中可以发现，时序模型的预测值变化轨迹滞后于观测值，这是由于时序模型未考虑气象因素的影响，所以只能捕捉持续性的变化，而无法充分反映AQI序列相邻2 d之间的剧烈变动.图 3b给出了这个时间窗口对应的标准化残差序列.从图 3b可以看出，在10月7—10日，标准化残差序列为较高的正值，反映了不利于污染物扩散的气象因素的影响.而在10月12日，标准化残差为绝对值较大的负值，反映了有利于污染物扩散的气象因素的影响.根据王占山等(2015)的分析，从7日持续到10日，主要的天气形势为弱气压和静风天气，还有偏南风导致的跨区域污染物传输，这些因素导致了污染物的累积.而从11日晚间开始的北风造成的污染物清除过程导致了空气质量迅速好转.这些过程都在图 3b中的标准化残差序列中得到了反映.

图 3(Fig. 3)

图 3 北京市2014年10月一次污染过程的日AQI序列的观测值和拟合值(a)及标准化残差序列(b) Fig. 3 Observations and fitted values of AQI series in a pollution process(a) in October, 2014 in Beijing and the corresponding standardized residuals(b)

需要指出的是，邻近城市AQI对气象因素的响应可能是区域性同步的.为了看出这种特征，本文给出了以北京市为中心的部分城市群在这个污染过程中的标准化残差的空间分布(图 4).从图 4可以看出，在从7日开始的污染增强阶段，大面积的城市有着同样符号为正的标准化残差；而在12日的污染物清除阶段，华北城市群的残差变为负值.这种在残差中存在的空间聚集模式在下文将使用Moran′s I指数来进行量化.

图 4(Fig. 4)

图 4 北京市2014年10月一次污染过程对应的标准化残差的空间分布 Fig. 4 Spatial distribution of standardized residuals corresponding to the pollution process in October, 2014 in Beijing

Moran′s I指数是表征全局空间自相关程度的指标.通过逐日的标准残差值，计算并给出了逐日的Moran′s I指数变化(图 5).考虑到2015年之前和之后可以使用的城市数量有差别(使用不同数量城市计算的结果会有差异)，且为了更充分利用数据，本文分别给出了以第一批次城市数据(从2014年1月1日开始)和以所有城市数据(从2015年1月1日开始)为基础的计算结果.

图 5(Fig. 5)

图 5 标准残差序列的逐日标准化Moran′s I指数(Z(I)值及其平滑曲线, 虚线为Z(I)=1.96的水平) Fig. 5 Daily Moran′s I series of the standardized residuals and corresponding smooth series(the dot lines mean Z(I)=1.96 level)

从图 5中可以看出，基于标准化残差的标准化Moran′s I指数(Z(I))序列存在着明显的以年为周期的变化规律.通过局部加权线性回归获得的低频变化曲线可以看出，每年冬季，Z(I)会达到峰值；而在每年夏季，Z(I)达到谷值.这种以年为周期的变化特征支持了本文的观点，即残差反映了气象因素对空气质量的影响.冬季较高的Z(I)指数反映了邻近城市AQI值对气象条件变化较高的同步响应程度；在夏季，这种同步响应的程度减弱.但需要指出的是，在一年四季中，这种同步响应都是显著的(Z(I)＞1.96).

由于全局的Moran′s I指数可能掩盖空间自相关特征的空间异质性，所以本文又计算了以标准化残差为基础的逐日局域Moran′s I指数.为了看清其季节平均特征，将逐日的标准化局域Moran′s I指数(Z(I))使用分位数进行季节合成.图 6给出了每个季节不同城市Z(I)值的累积频率P分别为0.7和0.6的样本分位数Q_0.7和Q_0.6.从样本分位数的含义可知，Z(I)小于Q_0.7的频率是0.7，而小于Q_0.6的频率是0.6.在图 6中，只标出了Z(I)大于1.64的城市(即在α=0.05的水平上呈显著为正的空间自相关).如果对于某个季节，某个城市的Q_0.6出现在图 6中，说明这个城市的局域Moran′s I指数在这个季节大于40%的天数中是显著的；而如果某个城市的Q_0.7出现在图 6中，则说明大于30%的天数是显著的.

图 6(Fig. 6)

图 6 标准残差序列的逐日标准化局地Moran′s I指数的P=0.7和P=0.6的分位数Q0.7和Q0.6(基于2014年12月31日之后的数据) Fig. 6 P=0.7 and P=0.6 quantiles(Q0.7和Q0.6) of the standardized local Moran′s I index of the standardized residuals

从图 6可知，在4个季节中，冬季的自相关水平较高.从Q_0.7来看，冬季局域Z(I)超过1.64的区域在冬季最大，几乎充满了整个中国东部地区.Z(I)值在4个季节中都有着一个类似的规律，即在包括长三角到京津冀的大片东部区域都有着较大的Z(I)值.相比而言，Q_0.6的空间分布特征呈现出更加明显的季节变化特征.如果将Q_0.6＞1.64的区域作为高自相关区，则其在春季的区域很小，集中在长三角附近.在夏季，这个高自相关区向北延伸到山东半岛，并明显向内陆延伸.秋季，京津冀地区开始显现高空间自相关特征.在冬季，高自相关区域的面积最大，其范围从京津冀地区延伸到长三角地区，并覆盖山西、河南和安徽等地.

从前面的局域空间自相关分析(图 6)，结合实际空气污染的程度，本研究认为有2个跨季节的面积较大的高空间自相关区域值得重视，其一是以京津冀地区为主的华北地区，另一个是长三角地区.这2个区域的空间自相关特征有着不同的季节变化特征，反映了不同季节气象条件对空气质量影响的特征.长三角地区AQI变化的空间自相关性在冬季和夏季最强；而京津冀地区AQI变化的空间自相关性在秋季和冬季最强.需要指出的是，高的空间自相关性并不代表高的AQI值本身，而是指邻近的城市AQI指数倾向于在气象条件的驱动之下同方向变动.

从首要污染物来看，京津冀地区在冬季的首要污染物是PM_2.5(图略).根据王冠岚等(2016)的研究，京津冀地区各地的污染物主要为原地生成.而根据王占山等(2015)的研究，跨区域的传输也是重要的影响因素.陈朝晖等(2007)研究了北京市一次重度污染过程中PM₁₀的变化规律，发现了空气污染程度与天气形势存在的对应关系，并强调并不是最高PM₁₀浓度的时刻对应最不利的天气形势，而是最高PM₁₀浓度增量的时刻对应着最不利的天气形势.这说明空气污染物是受到不利气象条件的影响而逐渐累积的，且天气形势影响的是空气污染物质浓度的增量.需要注意的是，残差(未标准化的)与增量有着很好的关联.从图 2a可见，对大多数城市的AQI序列来说，差分序列的自相关系数很低(大部分在-0.2~0之间).这说明差分序列中所包含的前后关联信息很有限，以差分序列进行时序建模，向前外推一步的预测结果必然十分接近0.所以，对于原始AQI序列来说，向前外推一步的预测必然是接近于前一日的AQI值，这就说明差分序列和残差序列之间一定有着很好的相关性.因此，本文研究AQI时序模型的残差与对API增量的分析有着内在的联系，本研究的思路可以看作是对增量分析的拓展.

从长三角区域来看，空间自相关性在冬、夏两季均较高(图 6).从首要污染物来看，长三角地区的冬季和夏季有所不同，冬季的首要污染物是PM_2.5, 而夏季的首要污染物是臭氧(王红丽, 2015)，长三角地区的臭氧浓度呈现出冬季低而夏季高的现象(刘芷君等, 2016).以往的研究指出，近地面臭氧浓度变化的支配因素是大气的光化学过程(谈建国等, 2007).根据对上海近地面臭氧污染的研究(谈建国等, 2007)，影响臭氧浓度的气象因素主要有温度、湿度、风向和风速等.高温和湿度较低的天气下，偏北风和较小的风速容易产生较高的臭氧污染.易睿等(2015)发现，长三角城市群各城市臭氧月浓度有着类似的变化趋势.本研究表明，从AQI指数来看，长三角城市群的日变化也有着较高的空间自相关.

对于长三角的冬季来说，以往研究发现，长三角地区的PM_2.5浓度有着明显的季节分布特征，呈现出冬季高而夏季低的规律(戴昭鑫等, 2016).对苏州和昆山PM_2.5的研究发现(江瑶, 2014)，气层稳定、气温逐渐升高、风速小和湿度低的综合作用，导致了空气污染物不易扩散，空气污染加重；而随着冷空气的到来，空气质量好转.这也为冬季长三角地区AQI残差的空间自相关性提供了一个解释.

1) 全局Moran′s I指数表明，中国城市群日AQI时间序列模型标准残差的空间自相关呈现出以年为周期的变化，冬季空间自相关程度增大，而夏季空间自相关程度减小.

2) 局地Moran′s I指数表明，有2个高空间自相关区值得重视，即京津冀地区和长三角地区.在春、夏两季，长三角地区是代表性的高自相关区域; 在秋季，京津冀地区成为代表性的高自相关区; 在冬季，从京津冀地区到长三角地区形成了一个联通的高空间自相关区.残差的高空间自相关体现了邻近城市空气质量对气象条件的较强的同步响应特征.