1 引言(Introduction)

大气污染是目前全球面临的主要环境问题之一, 其主要特点是污染物在流动时具有不可见性, 传播速度快、影响范围广、扩散能力强, 对环境、气候和人体健康产生重要影响(Mcconnell et al., 2010; Gurjar et al., 2010).据统计, 全球每年有高达300多万人因室外污染而过早死亡, 其中, 中国因空气污染而过早死亡的人数有近140万(WHO, 2014).随着人们对环境和健康问题的重视, 空气污染成为政府和人民普遍关注的热点问题(Yang et al., 2015).在2016年全国卫生健康大会上, 习近平总书记指出, 良好的生态环境是人类生存与健康的基础, 要按照绿色发展理念, 实行最严格的生态环境保护制度, 重点抓好空气污染的防治.华北地区是中国传统的政治中心和经济重心, 人口稠密, 近年来区域经济发展迅速, 城市化进程加快, 华北已经成为全国乃至全世界空气污染最为严重的地区( 张瑜等, 2015).2015年华北地区共发生1710次重污染及以上天气事件, 占全国的44.1%, 情势十分严重(新浪网, 2016), 因此, 对华北地区环境空气质量的分布特征及成因探究至关重要.

近年来, 国内外学者利用数理统计方法对空气质量展开了大量的研究, 主要集中于空气质量对人体健康和社会发展的影响(Lelieveld et al., 2015; Gauderman et al., 2015)、空气质量与气象因子之间的关系(肖建能等, 2016; 白鹤鸣等, 2015)、局部城市空气质量的成因(王占山等, 2015; 王冠岚等, 2016)、特定时段空气污染过程(潘竟虎等, 2016)、空气质量年际及季节变化(李名升等, 2013; 邓霞君等, 2013)等方面.也有学者通过立体观测和模型模拟对空气质量进行了定量研究, 例如, 利用遥感反演气溶胶对地面空气质量进行监测和分析(吴永红等, 2009); 利用第三代区域空气质量模式(Community Multiscale Air Quality, CMAQ)对京津冀、长三角、川渝地区的空气质量及相关污染物的浓度变化进行模拟, 并与地面观测数据进行了对比验证, 用于评价区域污染物的相互影响程度(程兴宏等, 2016; 黄晓波等, 2016; 李世广等, 2013).还有学者对中外空气质量评价标准的差异进行了研究并提出改进意见(高庆先等, 2015).虽然这些研究已经取得了一定的成果, 但将华北地区作为一个整体并分析其空气质量时空布局和综合成因的研究较少, 且空气污染作为一种空间现象, 现有研究多侧重于传统的数理统计, 缺乏对空间异质性和空间相关性的思考(刘永伟等, 2013); 同时, 空气质量是多种大气污染物的综合结果, 而现有成果一般侧重研究特定空气污染物, 如PM_2.5、PM₁₀等.

基于此, 本文以GIS空间分析技术为依托, 采用“点面结合”的思路, 在传统的数理统计基础上加入空间统计方法, 对华北地区空气质量指数(Air Quality Index, AQI)的时空分布规律进行探究, 并结合自然因子和社会经济因子进行成因分析, 以期为华北地区大气污染防治提供决策参考.

2 研究思路与方法(Research route and methods)

为探讨2015年华北地区AQI的空间分布特征及成因, 本文首先通过选取最合适的空间插值方法获得2015年全年及各季节的AQI空间分布数据, 然后分析降水、风速、风向及地形等对AQI的影响, 最后利用逐步回归分析法, 找出对AQI具有重要影响的社会经济因子, 并结合探索性空间数据分析(Exploratory Spatial Data Analysis, ESDA)方法和空间回归模型对AQI及社会经济因子进行空间相关分析.

2.1 GIS空间插值

研究区共收集到71个主要城市空气质量监测站点数据, 由于站点的数目有限, 本文通过空间插值方法获取华北地区的空气质量分布情况.空间插值常用的有反距离权重法、规则样条函数法和克里金法, 其中, 克里金法有多种经验半变异函数, 常用的有球面函数、圆形函数、高斯函数、指数函数和线性函数(Teegavarapu et al., 2012).为了选取较为合适的插值方法, 本文采用交叉验证对以上4种插值方法进行对比.具体方法为：总监测站点的75%(53个)作为训练样本进行插值, 用剩余的25%(18个)进行验证, 采用的评价指标包括平均误差(ME)、标准化平均偏差(NMB)、标准化平均误差(NME)及均方根误差(RMSE).上述4个指标的绝对值越小表明插值效果越好, 由于篇幅有限, 本文只列出年均插值验证结果(表 1).

验证结果表明, 克里金法的ME、NMB、NME和RMSE值普遍小于反距离权重法和样条函数法, 并以指数函数插值得出的结果综合最优.因此, 本文选用克里金法中的指数函数作为空气质量面状数据的生成方法.

2.2 ESDA技术

ESDA是以空间关联性为核心, 描述并揭示地理数据的空间分布, 分析其空间联系、集聚及异质性的空间分析方法(Kies et al., 2010).通过ESDA可以有效地检验地理数据的空间相关性, 一般采用莫兰指数(Moran′s I)进行度量, 判断出地理数据的属性分布是否具有统计上显著的集聚或分散, 以此来探讨数据的空间关联和空间差异程度(熊昌盛等, 2014), Moran′s I的数学计算公式如式(1) 所示.

(1)

式中, X_i和X_j分别表示城市i和城市j的AQI, 并且i≠j, ω_{i, j}表示城市i和城市j之间的空间权重, n等于要素总数.Moran′s I指数取值范围为.在给定显著性水平时, 若Moran′s I显著为正值, 表示AQI高值或低值的区域在空间上为显著集聚, 值越接近于1, 总体空间差异越小; 若Moran′s I显著为负, 则表示区域与周边地区的AQI具有显著的空间差异, 值越接近-1, 总体空间差异越大; 若Moran′s I为零值, 表示空间不相关.本文利用ESDA技术来进行AQI与相关社会经济因素的空间自相关分析.

2.3 空间回归模型

本文利用地理加权回归(Geographical Weighted Regression, GWR)模型对华北地区市域AQI进行空间回归分析.空间回归分析是在传统回归分析的基础上考虑了空间事物的自相关性, 通过空间关系把属性数据与空间关系结合起来, 以更好地解释地理事物的空间关系(张新峰, 2009).当市域AQI的影响因素的所有回归参数不受空间位置的影响, 保持固定不变, 则可使用最小二乘法(OLS)回归模型, 计算公式见式(2).

(2)

式中, y为因变量, x为自变量, β₀为回归常数, β₁, β₂, β₃, …, β_n为回归系数, ε为残差, n为研究城市总个数.

GWR模型是对最小二乘法(OLS)模型的改进, GWR模型充分考虑了影响因素的空间位置的变化, 将数据的地理位置嵌入到回归参数之中, 计算公式见式(3).

(3)

式中, u_i, v_i为第i个城市的空间地理位置坐标, β_j(u_i, v_i)是第i个城市上的第j个回归参数, 是地理位置函数(Sassi, 2010; 胡艳兴等, 2015).本文从OLS模型出发, 对华北地区各地级市的AQI空间结构特征进行检验, 然后构建GWR模型对影响AQI的社会经济因子进行模拟, 探讨区域差异及成因.

3 数据来源与处理(Data sources and processing)

2012年环保部规定采用AQI替代原有的空气污染指数(Air Pollution Index, API).AQI是定量描述空气质量状况的无量纲指数, 相比于API可以更好地表征环境空气质量状况, 参与AQI评价的主要污染物为细颗粒物(PM₁₀)、可吸入颗粒物(PM_2.5)、二氧化硫(SO₂)、二氧化氮(NO₂)、臭氧(O₃)、一氧化碳(CO)等6项.AQI划分为6个等级, 数值越大、级别越高, 说明空气质量越差; 数值越小、级别越低, 空气质量就越好(王艳琴, 2012)(表 2).本文的AQI数据来源于中华人民共和国环境保护数据中心2015年公布数据, 共30031个样本数据; 行政区划数据来源于1 : 400万中国基础地理信息数据库; DEM来源于国际农业研究磋商组织(CGIAR)官网发布的SRTM 90 m分辨率的全球数字高程数据; 气象资料数据来源于中国气象科学数据共享服务网、中央气象台官方网站; 社会经济数据来源于2016年中国统计年鉴.根据分析需要, 对研究区2015年1—12月的日AQI求月平均值, 然后计算出各季度平均值, 对于缺失值的日期, 用其相邻日期数据的平均值进行代替.运用ArcGIS软件将各月和各季度的AQI平均值作为对应城市点的属性值, 建立GIS数据库.

4 结果与分析(Results and analysis) 4.1 华北地区空气质量总体特征

对研究区71个城市的AQI等级日数分布进行多组独立样本Kruskal-Wallis H检验, 可知p=0.000 < 0.05, 均具有统计学意义(表 3).从表 3中可以看出, 在华北地区AQI等级平均日数构成中, 空气质量主要为良, 占整体的49.2%;71个城市的AQI值范围为14~500, 年均值为106.52, 已超过国家AQI二级标准(AQI < 100) 的临界值.有70个城市发生过中度污染、重度污染和严重污染天气共3959 d, 即平均每个城市57 d, 日数占整体的15.4%, 其中, 北京市(30.4%) > 河北省(19.9%) > 河南省(18.2%) > 天津市(14.9%) > 山东省(14.4%) > 山西省(8%).对比2015年全国其他片区监测站点的AQI大于150(即中度污染以上天气)的日数比例：华北地区(15.4%) > 西北地区(8.3%) > 华中地区(8.2%) > 东北地区(7.1%) > 华东地区(5.3%) > 西南地区(2.7%) > 华南地区(1.5%), 可知华北地区的空气污染在全国居首并远高于其他片区.

4.2 AQI年内时空分布特征

通过已建立的GIS数据库中的各月监测点均值数据求得季度均值和年均值, 利用GIS中的克里金指数模型插值方法对华北地区71个监测点进行空间插值, 得到华北地区2015年年均AQI空间分布图和4个季度的AQI空间分布图(图 1).从全年尺度看, 华北地区空气污染分布范围较广, 重污染区域集中在北京、河北南部、河南北部和山东西部, 山西西部、河北北部和山东半岛空气质量较好; 季节尺度上, AQI均值为70~200, 并以70~150为主, AQI均值为150~200的区域只出现在冬季, 且分布范围广泛.同时, 空气污染呈现出的季节变化趋势较为明显, 冬季污染范围分布最广, 高值区域集中成片, 只有山西西北部、河北北部和山东东部地区环境较好; 春季AQI均值降低, 污染范围收缩变小, 高值出现分离; 夏季是全年空气质量最好的季节, 大部分地区的AQI均值都未超国家二级标准; 秋季污染范围又开始扩张, 高值在河北南部和山东西部聚集.

4.3 AQI影响因素分析 4.3.1 AQI与气象要素

气象要素是制约污染物在大气中稀释、扩散、迁移和沉降的重要因素(赵普生等, 2012).根据上文分析结果可知, 华北地区AQI随季节温度、湿度变化产生的波动性较强并在空间分布上存在明显的差异.为进一步分析气象要素中的风向、风速、降水对AQI产生的影响, 本文采用迭代自组织数据分析技术(Iterative Self-Organizing Data Analysis Technique, ISODATA), 以年均AQI克里金插值数据作为分类变量, 找出华北地区空气核心污染区(图 2), 并以核心污染区中的城市作为空气污染较高的典型对象进行分析.

根据相关研究, 降水在任何季节都对空气中的颗粒物质具有冲刷作用, 能够改善空气质量(韩燕等, 2012; 周国兵等, 2013).但由于降水对污染物的作用比较缓慢, 本文以降水日当天的AQI和前一日AQI的变化率来分析降水对空气中颗粒物的清除作用.利用城市降水量的逐日数据及对应AQI较前日的变化率, 通过散点图分析降水量对AQI变化率的影响(图 3).从图 3中可以看出, 日降水量在0~10 mm分布居多, 对应的AQI变化率中负值占56.1%, 略多于正值数, 说明降水对颗粒物有一定的冲刷作用, 但不太明显; 当降水量超过10 mm时, AQI变化率绝大部分为负值, 占83.1%, 并且变化率较大, 说明降水量达到10 mm以上, 对空气的净化作用较明显.

风是影响污染物稀释扩散的重要因子, 风速是造成快速水平输送或平流的主要原因(李宗恺等, 1985).对2015年各个季节的逐日风速数据与其对应的AQI数据及AQI评价污染物浓度进行相关性分析(表 4).结果显示, 在空气质量较好的夏季, 风速与污染物之间的相关性不太明显, 只有部分污染物通过了显著性检验, 而在春季、秋季和冬季, 风速与众多污染物之间呈现显著的负相关; 从全年尺度上看, 除了O₃未通过显著性检验外, 其余参与AQI评价的污染物与风速都呈显著的负相关, AQI与风速的相关系数为-0.201, PM_2.5、PM₁₀、SO₂、CO、NO₂与风速的相关系数分别为-0.255、-0.150、-0.156、-0.216和-0.280, 说明风速的增大使得单位时间内通过烟波断面的空气量增大和湍流扩散增强, 对华北地区日常发电、工业生产、汽车尾气等经济生产过程中的污染排放物起着稀释作用, 而O₃只在夏、秋两季与风速具有负相关性, 主要原因是臭氧是大气中的氮氧化物和挥发性有机物经光化学反应形成的光化学烟雾的主要成分, 主要集中在日照充足的夏季, 因此, 在夏、秋季节, 风速的增大有利于氮氧化物和挥发性有机物扩散, 而就全年尺度, 风速对其影响不太明显.值得注意的是, 春季风速对PM₁₀的抑制作用不够显著, 这与华北地区春季回暖和强风天气有关, 风速过大容易引起扬尘天气, 增加空气中的PM₁₀浓度, 降低了相关系数.

以空气污染最严重的12月为例, 分析风向在空气污染过程中产生的影响.12月各地的AQI波动较大, 存在3次空气污染过程, 分别是12月3—12日、12月12—16日、12月16—27日(图 4).以污染持续最长的12月16—27日为时间段绘制城市玫瑰风向图(图 5).时段平均风速为8.1 km·h^-1, 主要盛行风向为北风和偏北风, 出现频率合计为41.9%.

结合该时间段中核心污染区主要城市的AQI级别日变化(图 6)可知, 12月16日各城市空气质量较好; 17日北部城市出现污染天气, 污染集中在北京和保定附近, 而西南和南部城市空气质量较好; 18日北部城市污染进一步加重并开始向东南扩散, 保定出现重度污染, 衡水、德州和聊城也出现中度污染天气, 只有西南部的枣庄、临汾、莱芜和东营空气质量为良好; 19日污染继续向南部和东南部扩散, 山东西部城市大部分都出现了重度污染天气, 空气质量为良好的城市也转为轻微污染; 20日开始各城市的空气污染按由北至南的顺序不断加重; 到12月24日, 区域内绝大部分城市空气质量都为严重污染; 25日北部和西北部城市空气污染有所降低, 严重污染集中在山东西部地区; 26日各城市的污染开始消退; 27日空气污染消退, 各城市空气质量转为良好.综合空气污染由北至南的动态变化过程和偏北盛行风传输方向可以看出, 污染物随盛行风向引起的大气流动进行传输, 再叠加本地的源污染, 形成空气污染的空间移动, 造成了大范围的空气污染天气.

4.3.2 AQI与地形

华北地区西部是南北走向的太行山山脉, 北部是东西走向的燕山山脉, 地形抬升使得西部和北部地形起伏明显高于东部平原.通过对研究区的DEM进行起伏度分析, 将其重分类为10级, 并统计各级的AQI均值水平(图 7).发现地形起伏度为1级(0~13 m)的AQI均值最高, 2级之后的AQI均值波动不大, 说明华北地区空气污染主要集中在平原地区, 在起伏度较大的山地和高原地区污染较轻.

为探讨山脉对AQI的影响, 以太行山和燕山山脉为中心向两侧作100 km的缓冲区, 分别对AQI年均栅格进行分区统计(表 5).太行山山脉东西部及燕山山脉的南北部AQI均存在较大差异：太行山山脉东侧 > 西侧, 燕山山脉南侧 > 北侧, 并且太行山山脉东侧的最小值大于西侧的最大值, 燕山山脉南侧的最小值大于北侧的最大值.主要原因是太行山山脉东部和燕山山脉南部位于平原地区, 城市密集, 人口基数大, 内部污染物排放量大(张殷俊等, 2015), 而污染物在向太行山和燕山山脉方向的扩散过程中受到地形抬升的阻塞, 容易形成闭合环流风, 对气流产生阻挡效应, 污染物难以越过起伏较大的山脉, 因此, 空气质量出现了较大分化.

4.3.3 AQI与社会经济因素

社会经济与空气污染之间的关系较为复杂, 需要进行多因素综合分析.本文选取的社会经济因子包括总人口(万人)、人口密度(人·km^-2)、地区GDP(万元)、第二产业产业产值比重、煤气及液化石油气供应量(万m³)、民用汽车拥有量(辆)、城市建成区绿化率、森林覆盖率.由于这8个指标是按各省市的地级市行政区统计得来, 为方便分析社会经济因素和AQI的关系, 本文利用GIS中栅格数据的空间统计方法, 对AQI年均值进行地级市行政区空间统计, 通过经济因子与AQI进行逐步多元线性回归分析可消除多变量间的多重共线问题.回归结果表明, 人口密度(p=0.003)、第二产业产值比重(p=0.000)、民用汽车保有量(p=0.000)、森林覆盖率(p=0.001)4个因子通过了显著性检, 因此, 剔除其余4个因子.为探讨AQI与社会经济因子的空间自相关性, 计算出对应的Moran′s I (表 6).

由表 6可知, 所有变量的Moran′s I指数均大于0, 且均通过了显著性检验(p < 0.05, Z > 1.96), 说明AQI、人口密度、第二产值比重、民用汽车保有量和森林覆盖率并不是随机分布, 而是存在空间自相关性, 且表现出较为显著的空间集聚, 因此, 可以利用空间回归模型定量解释不同因子对AQI产生的影响.为消除量纲的影响, 将因变量AQI和经济因子进行标准化处理, 然后利用最小二乘法(OLS)模型对华北地区AQI进行全局回归分析(表 7).

由表 7可知, 联合卡方统计量p值小于0.05, 表明模型具有统计显著性; 影响AQI的4个经济解释变量中, 除了地区GDP外, 其余3个因子的p值均小于0.01, 与AQI呈现正相关关系, 4个因子的VIF值均小于7.5, 因此, 都可以作为模型的解释变量; 同时还得到OLS回归模型的Adjusted R²为0.688026, AICc值为105.454746, 根据稳态评估系数Koenker(BP)统计量的p值小于0.05可知, 模型具有统计学上的显著性非稳态, 因此, 适合进行地理加权回归分析.通过GWR模型得出的R²为0.89219, Adjusted R²为0.83373, 大于OLS模型, AICc值为85.48035, 与OLS模型的AICc值差距大于3, 表明GWR的拟合结果比OLS更加准确; 通过对OLS和GWR模型结果的残差进行空间自相关性分析, GWR模型残差的Moran′s I指数为0.046392, 小于OLS模型残差的Moran′s I指数0.211412, 表明GWR的残差自相关性比OLS小, 再次证明用GWR模型拟合AQI和经济因子之间的关系优于OLS模型.

根据GWR模型各因子的综合回归系数, 它们对AQI的影响程度大小依次为：第二产业产值比重(0.34) > 民用汽车保有量(0.31) > 人口密度(0.28) > 森林覆盖率(-0.22).采用Jenks自然断点法(Natural Breaks)对华北地区GWR模型的局部系数进行可视化(图 8), 可以看出, 不同的社会经济因子对AQI的影响存在空间差异.在几个解释因子中, 第二产业产值比重对AQI的影响最为显著, 表明工业经济活动是影响华北地区空气质量的主要原因, 回归系数高值主要集中在华北中部地区, 包括河北南部、河南北部和山东西部和山西北部, 该区域是华北地区的重工业基地, 经济增长在很大程度上依托第二产业产值的贡献, 因此, 第二产业产值比重回归系数呈现高值集聚.民用汽车保有量对空气质量的影响仅次于第二产业产值比重, 说明汽车尾气的排放对华北地区的空气质量也造成了较强的影响, 回归系数高值集中在华北中西部地区, 主要原因是该地区近年来随着经济发展, 汽车保有量增多, 并且地处季风气候和大陆性气候的交界处, 汽车排放的尾气难以沉降, 因此, 对AQI的影响较大.人口密度与AQI也呈现出空间正相关关系, 表明人口的集聚会加重空气污染, 回归系数的高值集中在华北南部地区, 主要原因是该地区人口基数多、密度大, 大量的人口集聚对空气质量产生影响.森林覆盖率是唯一一个对AQI起抑制作用的因子, 说明森林覆盖率的提高有助于缓解空气污染, 其回归系数整体呈现西高东低, 最低值出现在山东半岛区域, 主要原因是华北东部地区森林覆盖率普遍高于西部地区, 因此, 东部地区森林覆盖率对空气质量的抑制作用也大于华北西部地区.

5 结论(Conclusions)

1) 华北地区AQI分布存在较大的差异, 且呈现出明显的季节变化趋势：冬季>春季>秋季>夏季.冬季污染范围分布最广, 高值区域集中连片, 夏季是全年空气质量最好的季节, 大部分地区的AQI未超国家二级标准; 华北地区空气污染整体分布范围较广, 重污染区域集中在京津冀经济圈、河北南部、河南北部、山东西部, 而山西西部、河北北部和山东半岛空气质量较好.

2) AQI与较为典型的气象要素如温度、气压、降水、风向和风速的关系较为密切.当降水量大于10 mm时, 对AQI具有明显的抑制作用, 有利于空气质量的提升; AQI与风速呈显著的负相关, 风速的增大有利于污染物的稀释, 风向在空气污染过程中影响污染物的传输方向, 打破了污染源的孤立性, 从而扩大了空气污染的影响范围; 华北地区AQI与地形之间也存在一定的相关性, 空气污染主要集中在平原地区, 并且太行山脉和燕山山脉对气流产生阻挡, 阻塞污染物的扩散, 使太行山脉和燕山山脉两侧地区的AQI出现了较大分化.

3) AQI与社会经济因子均具有空间自相关性, 且表现出较强的空间集聚效应, 说明一个地区的AQI与相邻地区的AQI不是孤立的, 而是相互联系的.GWR模型中各社会经济因子的回归系数表明, AQI与第二产业产值比重、民用汽车保有量、人口密度呈正相关, 与森林覆盖率呈负相关, 其中, 第二产业产值对空气污染的贡献最大, 其次是汽车保有量.

6 建议(Suggestions)

根据本文研究成果, 建议华北地区有关部门需要重点防控冬、春季节的空气污染, 并在城市规划中重视风道的建设, 在拥有成雨条件的雾霾天气时可考虑采取人工降水的方式加快空气污染物的稀释; 其次是贯彻绿色发展的理念, 转变华北地区的经济增长方式, 加快第二产业的经济转型和模式创新(关成华等, 2016), 加强对工业企业的排放监管, 对超标排放的企业进行整治; 然后需要控制机动车的数量, 对机动车的排放标准进行严格的监测; 最后推进海绵城市的建设, 提高森林覆盖率, 并加大对空气污染防范和治理费用的投入.