城市空气污染监测数据的统计建模与探究

发布者:陆敏发布时间:2021-05-25浏览次数:16


一、摘要


       空气质量是近年来备受关注的话题,PM2.5作为衡量环境空气质量的重要指标,具有重要的研究意义。本文选取我国31个城市的空气污染监测数据,利用R软件分别进行了相关性分析、回归分析、聚类分析以及过程监控。首先,文章分别利用普通最小二乘法、主成分分析法、偏最小二乘法和BP神经网络四种分析方法对所选取的不同城市六项污染物浓度进行分析,探究了空气质量指数AQI中其他五项污染物浓度对PM2.5浓度的影响。其次,基于弱模式聚类法,对各个污染程度不同的城市进行了分类,为研究不同区域的污染状况以及采取差异化的空气污染治理措施提供了依据。最后,选取各类代表性城市的空气污染数据建立时间序列模型,并用带自相关数据的CUSUM控制图方法进行过程监控,对同一污染物的浓度在不同时间的数据进行时间序列分析、污染趋势分析以及过程监控,为掌握环境空气质量指数变化趋势和及时调整环保策略提供决策支持。


关键词:PM2.5;偏最小二乘法;主成分分析;BP神经网络;弱模式聚类;AR模型


二、研究主要流程


      1. 研究背景

      2. 数据与模型介绍

      3. 建立模型

      4. 结论与建议


三、分析及主要结论


(一)背景概要


       PM2.5 是指直径小于或等于 2.5 微米的空气颗粒,也称可入肺颗粒物,入肺一词直截了当的点明其危害,由于其体积比较小、质量比较轻,从而能在空气中长时间停留,增大进入人体的几率,从而成为研究空气污染颗粒物中的重点对象。不仅如此,相比其他大气颗粒,PM2.5 更容易富集空气中的有毒重金属、酸性氧化物、有机污染物、病毒和细菌。并且颗粒物的粒子半径越小,其化学成分就越复杂,对人体健康的危害也就越大。此外,PM2.5 对城市能见度也有着非常大的影响。

       传统衡量空气质量的指标是空气污染指数(API),是根据近地面的几种主要空气污染物浓度及其持续时间确定的,目前计入 API 的指标项目有 SO2、NO2、PM10,然而灰霾的主要原因-- PM2.5 并未进行考量。2012 年 2 月 29 日,中国国家环境保护部公布了最新修订的《环境空气质量标准》(GB3095—2012)和《环境空气质量指数(AQI)技术规定(试行)》(HJ633-2012), 启用了空气质量指数(AQI)作为空气质量监测指标并于 2016 年 1 月 1 日在全国实施。该指标代替了API,是无量纲指数,其参与评价的污染物有细颗粒物(PM2.5)、可吸入颗粒物(PM10)、二氧化(SO2)、二氧化氮(NO2)、臭氧(O3)、一氧化碳(CO)等六项。PM2.5 进入公众视线的时间不长,但由于 AQI 采用了更严格的标准,发布频率更高,且首次将产生灰霾的主要因素--对人类健康危害极大颗粒物PM2.5 的浓度指标作为空气质量监测指标,反映了短期的空气质量状况,可以为公众的出行和健康提供指导,因此研究 AQI的六项污染物指标,尤其是 PM2.5,具有重要意义。



(二)研究的问题


      1. 对31个省会城市和直辖市(港澳台除外)的空气质量数据,通过普通最小二乘法、主成分分析法、偏最小二乘法、BP神经网络四种分析方法对AQI中的六项污染物浓度进行分析,探究了AQI中其他五项污染物浓度对PM2.5浓度的影响。

      2. 对选定的以上城市的空气质量数据,采用弱模式聚类法对污染程度不同的城市进行了分类,为研究不同区域的污染状况以及采取差异化的空气污染治理措施提供了依据。

      3. 选取各类代表性城市的空气污染数据建立时间序列模型,并用带自相关数据的CUSUM控制图方法进行过程监控,为实时应对空气质量变化提供了依据。



(三)数据与模型概要


      1. 数据简介

       本文所有数据来自中国环境监测总站官方网站,数据为 2015 年 10 月 1 日到 2016 年 3月 31 日我国27个省会城市(港、澳、台除外)和4个直辖市的日AQI 指标数据,包括可入肺颗(PM2.5)、可吸入颗粒物(PM10)、二氧化硫(SO2)、二氧化氮(NO2)、臭氧(O3)、一氧化碳(CO)。AQI 中六项指标说明和原始数据示例如下:


AQI 各项指标说明


原始数据示例


      2. 数据相关性分析

为了解变量间的相关关系以保证模型的科学性,在做回归分析之前,本案例首先对所有变量进行相关性分析。这里首先通过变量的散点图矩阵对各变量进行初步判断,结果见下图。

全国 AQI 六项指标的散点图


中可以大致看出,除了和O3SO2以外,PM2.5AQI其他各项指标之间大都呈现出较高的相关性。与此同时,AQI各项指标两两之间的相关性大多较低,但是PM10NO2PM10SO2之间的相关程度较高,因此若直接对各项指标进行回归,则变量之间可能存在多重共线性。当自变量的集合内部存在较高程度相关性时,称自变量存在多重共线性,多重共线性的存在常常会扩大模型系数的估计误差。为了进一步确定对数据相关性的判断,我们做出了各项指标之间的相关系数表进行定量判断,见下表


AQI 各项指标之间的相关系数

中可以清楚地看到,除了O3SO2以外,PM2.5浓度和其他各项指标浓度的相关系数都超过了0.7,具有高度相关性,认为可以利用AQI中的各项指标对PM2.5浓度存在较密切的关系。值得注意的是,PM2.5O3呈现负相关关系,有待进一步的验证。此外,PM10NO2PM10CO之间的相关系数都在0.7以上,说明变量间存在多重共线性。

由于空气污染物在不同城市的分布可能存在较大差异,因此本案例以南京市为例做了单个城市的AQI各项污染物浓度的相关性分析,下图给出了南京市2015101日到2016331AQI各项指标的相关性散点图矩阵:

南京市 AQI 各项指标的相关性散点图


上图种看出,相比于31个城市的整体数据,南京市PM2.5浓度与其他五项污染物浓度呈现出更加明显的相关关系,尤其是与PM10近乎线性相关,并且PM2.5浓度与O3浓度仍然呈现负相关关系。此外,除PM2.5之外的五项指标都呈现出明显的正相关关系,说明如果进行普通的回归,则变量之间会存在多重共线性。类似地,我们通过相关系数表进一步确定变量间的相关性,见下表


AQI 各项指标之间的相关系数(南京市)

PM2.5PM10的相关系数高达0.93,和CO的相关系数高达0.74,和NO2SO2也呈现出较高的相关性,而和O3相关程度不是很高,但是也呈现出负相关关系。由于不同城市的经济发达程度、工业发达程度、人口数量等因素各不相同,所以会造成每个城市的各项空气污染物的含量也存在差异,从而导致这些污染物浓度之间的相关关系不同。总体而言,PM2.5O3存在负相关关系,和其他污染物存在不同程度的正相关关系。


      3. 预测模型建立与分析


(1)基于普通最小二乘法的回归

      我们首先运用最小二乘法对 2015 年 10 月 1 日到 2016 年 3 月 31 日的日 AQI 的所有指

标数据进行多元线性回归,数据标准化后的回归结果如下表:


普通最小二乘回归结果

       模型整体显著,且修正后的可决系数也较高,但是变量O3在 0.05 的显著性水平下并不显著,我们推测这可能是多重共线性导致,因此我们进一步采用可以在很大程度上消除多重共线性的模型—主成分回归模型、偏最小二乘回归模型进行拟合,最后再采用机器学习算法—BP 神经网络模型进行拟合。


(2)基于主成分分析法的回归


       对标准化后的数据计算各主成分的贡献率和累计贡献率,其中,PV 表示方差贡献率,CP 表示累计贡献率。


主成分分析结果(OLS)

       从上表的主成分分析可知,前两个主成分的累计贡献率达 75.85%,前三个主成分的累计贡献率达89.79%,前四个主成分的累计贡献率已经达 95.88%。

       结合累计方差贡献率和碎石图,分别选择对三个成分(累计贡献率 89.79%)、四个成分(累计贡献率 95.88%)进行回归,然后比较回归模型的优劣。

碎石图


三个成分回归(OLS)


四个成分回归(OLS)

最后对还原主成分回归模型的变量系数,得原始变量关于 PM2.5 浓度的回归模型如下:


       PM2.5 与 SO2、CO 都呈现出正相关关系,这与现实中 PM2.5 主要由人为排放,在空气中转化成 PM2. 5 的气体污染物主要有二氧化硫、氮氧化物、氨气和挥发性有机物等的事实也符合。


 (3)基于偏最小二乘法的回归


偏最小二乘回归结果

中看出,协变量前三个因子的累计贡献率已经达到82.62%,前四个因子的累计贡献率已经达到94.1%CV表示的是交叉验证结果中的PRESS。这里通过5个因子(Z1-Z5)和PM2.5X1)两两之间的相关系数散点图见下图,发现5个因子之间已经没有明显的相关关系,观察第一列可知X1与其他因子的线性关系比较明显,即被提取协变量的因子与响应变量的相关性较高,说明偏最小二乘法因子提取效果较好。

各因子和 PM2.5 相关关系散点图



这里选用Z1Z2Z3Z4四个因子进行回归。


四个因子回归(PLS)


       最后,用 R 语言“pls”程序包对因子回归模型进行系数还原,得到原始变量的回归系数和回归方程:


模型中,PM2.5PM10NO2CO都呈现出正相关关系,这与现实中PM2.5主要由人为排放,以及在空气中转化成PM2.5的气体污染物主要有氮氧化物、氨气和挥发性有机物等的事实也符合。与O3SO2呈现负相关关系。


4基于 BP 神经网络的回归

       在此处使用 R 语言 nnet 包中 nnet()方法,该方法是常见的前馈反向传播神经网络算法,参数设置隐层单元个数为 10,权重调整速度为 0.01,最大迭代次数为 10000 次,采用线性输入。


(5)选用不同模型的结果比较分析


四种方法的 CV 值

结果分析:

       仅仅从 CV 值的大小而言,BP 神经网络模型的拟合效果在四种方法中最优,更适合进行数据拟合的模型是偏最小二乘回归模型,因为它可以消除或在很大程度上降低变量间的多重共线性,从而使拟合值更加可信。


      4. 城市空气污染程度的聚类分析

聚类结果:


类 1(空气质量一般):北京、成都、广州、哈尔滨、杭州、合肥、呼和浩特、兰州、南昌、南京、南宁、上海、沈阳、太原、天津、武汉、重庆、西安、西宁、银川、长春、长沙

类 2(空气质量较差):济南、石家庄、乌鲁木齐、郑州

类 3(空气质量较好):福州、贵州、海口、昆明、拉萨


      5. 南京市 PM2.5 的过程监控

       由于污染物浓度存在较大的自相关性,且大部分污染物浓度数据不满足正态性假定,这使得传统的统计过程控制(SPC)方法无法适用。为了解决这一问题,本案例尝试用 AR 时间序列模型对污染物浓度进行建模,其次采用自相关数据下的 CUSUM 统计过程控制图方法对其模型进行监控。

       在第一类城市中以南京市为例,对 PM2.5 浓度数据进行过程监控。数据取自 2015年10月1日到2016年3月31日共 183 组观测。根据我们的分析发现该数据在前 40个时间点未发生均值变化,于是我们采用前 40 个数据进行第一阶段的模型估计,得到模型参数的估计为(0.7753,-0.2823),然后对模型进行 Shapiro-Wilk 正态性检验,结果如下:


Shapiro-Wilk 正态性检验(南京)


       通过 AR(2)模型建模使得残差通过了正态性检验。此时残差的 CUCUM 控制图调整参数选为 k=0.5,控制线h1=3.502。这是正态分布下使平均运行时间 ARL0=200 的参数设定。另一方面,原始观测数据不满足正态性假定,则不能使用这个调整参数,为此我们考虑用 MC方法将模型参数代入 AR 模型计算 CUSUM 控制图的平均运行时间 ARL0,模拟得到原始观测的参数为 k=0.5,h2=9.86。最终的过程控制图如下:

(a)南京市 PM2.5 观测数据,(b)AR 模型残差图,

(c)观测数据的 CUSUM 控制图,(d)残差的 CUSUM 控制图


       可以看出该组观测数据在第 80 天左右检测到均值和残差的变化,在 160 天左右检测到残差的变化。从监控的整个变化趋势来看,整体反映出 PM2.5 在冬季浓度较高,特别是在温度比较寒冷的 12 月份、1 月份和 2 月份均值变大较大。环境相关部门可根据监控的信号及时做出积极的应对措施,从而降低 PM2.5 浓度,改善空气质量指数。


(四)建议


      1. 削减硫氧化物和氮氧化物的排放总量--限制高硫煤的开采和使用、防治工业生产过程中产生的二氧化硫、二氧化氮污染、限制汽车尾气排放等。

      2. 限制碳排放量--规范燃料的使用、改善能源结构,改进技术并提高燃料的转化率和利用率、限制森林砍伐、防范火灾等。

      3. 环境相关部门可根据监控过程的信号及时做出相应改善措施,从而降低 PM2.5 浓度,改善空气质量指数。