基于R语言的哈尔滨市大气污染数据时间序列分析*

罗 爻，刘 硕，高俣晗，张美琦，宣立强

(哈尔滨师范大学)

0 引言

近年来中国快速城市化和工业化进程，造成了严重的大气污染.大气污染引发的一系列环境问题，严重危害了人类的身心健康[1-2].随着对环境和健康问题的重视，大气污染已经成为中国需解决的重大民生问题和国家发展问题[3-4].探讨哈尔滨市大气污染时空变化特征，有利于引导城市实施可持续发展战略.

随着大数据时代的到来，传统的统计分析软件难以满足海量数据的处理需求[5].因此出现了诸多新的技术与工具，其中开源统计分析语言R被广泛应用[6].R语言是一种由统计学家开发的统计计算和绘图的语言和环境，具有大数据处理的能力[7].为响应国家要求哈尔滨市于2015年底前完成了空气质量监测系统的建立，对重点污染企业实施严格的在线监控.要实现大气污染物网络化治理，必须充分发挥大数据的特性[8].大气污染监测数据量大，传统的数据分析工具耗时费力，不便于高效、快捷的处理和分析数据，难以发挥监测数据在污染防治中的作用，R语言及其众多工具包为R语言用于大量大气数据分析和可视化提供了强有力的支持.

早期的时序分析使用直观数据来寻找规律，随着研究进展，发现简单的描述时序分析具有很大的局限性[9].用统计学原理来分析时间序列，可以更准确的估计随机序列的演变[10].该文应用R语言及时间序列分析技术对哈尔滨市2017年全年监测数据进行各种可视化分析，以探讨R在大气数据分析领域应用的巨大潜力.并为相关部门制定大气污染控制措施提供科学依据.

1 R语言与所需工具包(package)的安装与功能

R语言官方主页https://www.r-project.org/ ，R语言开发环境下载安装地址：https://cran.r-project.org/mirrors.html，辅助R的工具RStudio，下载地址： http://www.rstudio.com/ide进入下载页面后，有Desktop和Server 2个版本，选择Desktop.下载安装R与RStudio后，使用install.packages()安装所需工具包，library()加载所需工具包.所需工具包与功能见表1.

表1 工具包功能

2 数据处理

2.1 数据来源

该研究数据来源于黑龙江省生态环境厅(http://www.hljdep.gov.cn/)，通过Python编写接口程序，自动获取哈尔滨市12个监测点位CO、NO2、O3、PM2.5、PM10、SO2的小时数据.该文应用R语言以及相关工具包对哈尔滨市2017年12个测站点的小时数据进行时间序列分析，取12个监测站均值.以探讨该软件在空气质量数据分析领域应用的巨大潜力[11].

2.2 数据处理

使用Excel表格对监测数据进行筛选排查，对数据进行简单的分列处理，并导出为CSV格式文件以便R读取.通过read.CSV()命令导入R中，利用TimePoint函数对检测数据进行时间格式转化，sqldf包对检测数据进行数据汇总与计算平均值，运用summaryPlot函数快速概览数据整体情况，利用ggplot2函数的分面功能展示各污染物的时间变化，最终进行数据分析.

3 时间序列分析技术

时间序列分析是大气监测数据的常用分析方法，包括对季度、月份、周-日、小时变化等特征进行分析，以揭示污染物的时间变化规律以及预测变化趋势.时间序列分析的目的一般有两个方面：一是认识产生观测序列的随机机制，即建立数据生成模型；二是基于序列的历史数据，也许还要考虑其他相关序列或因素，对序列未来的可能取值给出预测或预报[12].该研究中主要用到滑动平均过程.一阶滑动平均过程公式：

Yt=et-θet-1:

E(Yt)=0

(1)

ρ1=(-θ)/(1+θ2)

γk=ρk=0k≥2

方程式中θ和ρ1的一些数值可以帮助说明各种可能性.需要注意的是，负的θ对应的ρ可以通过简单地取正的θ所对应的ρ的负数得到.

3.1 时间序列分析

为了解哈尔滨市6项污染物的分布情况，对2017年监测数据运用summaryPlot函数快速概览数据整体情况，绘制出图1.对污染物分布情况进行概述统计分析，并参照国家空气质量标准(GB3095 2012)[13]，对各项污染物超标情况进行描述.从图1、表1中可看出，除SO2外，PM2.5、PM10、NO2年均值全都超过了国家空气质量2级标准，其中PM2.5超标最严重.将数据集进一步整合成适用于程序语言处理的形式(由宽变长)，利用ggplot2函数的分面功能展示各污染物的时间变化序列如图2所示.

图1 6项污染物的时间序列图

表1 环境空气污染物基本项目浓度限值

图2 各污染物分列的时间序列变化

SO2、CO全年排放量不超标，但秋冬季排放量高于春夏季.NO2排放量超标率为7.8%，污染物浓度起伏不大，峰值出现在冬季.PM10与PM2.5污染物浓度起伏规律相似，在10月后污染物浓度达到峰值，秋冬污染物浓度值高于春夏.PM10、PM2.5超标率分别为23.1%和69%，从图2可看出，SO2、CO、PM2.5、PM10、NO2这几项污染物浓度值都是秋冬季大于春夏季，其主要原因秋冬季为哈尔滨市的采暖期，气温O℃以下，早晚为燃煤供暖和出行的高峰期.寒冷天气路面易结冰，机动车行驶速率降低，在外滞留时间延长，污染物排放也变高[14].哈尔滨市冬季受内蒙古-西伯利亚高压控制，下沉逆温易出现并时间长、范围广、势力强，不利于污染物扩散.且正值秋季农忙结束，农民大量燃烧秸秆等生物质，导致空气严重污染.O3呈锯齿状对称分布，1～4月逐渐升高，在5月达到最大值，8月出现一个次峰，此后浓度值逐渐下降.

从图2中可以看出PM10与PM2.5为哈尔滨市大气污染主要成分，timeVariation函数提供了将不同量级的数据进行标准化处理的方法，该函数还提供了计算充值的功能，利用difference这一参数计算粗颗粒物浓度并展示其时间变化如图3所示.

图3可反应出粗颗粒物的变化特征，月变化中显示4、5月份浓度较高，说明春期扬尘污染较严重，日变化中8：00-9：00与19：00-20：00出现峰值，说明道路交通扬尘是该点位粗颗粒物的重要来源.周变化中周末平均小时浓度低于工作日浓度，表现出了明显的周末效应.

3.2 所需代码

读取数据

raw_data=read.csv('2017.csv', head = T)

时间格式转化

raw_data$TimePoint=as.POSIXct(strptime(raw_data$TimePoint,format='%d/%m/%Y %T'))

计算站点的平均值作为最终画图的数据

图3 PM10、PM2.5以及其插值的时间变化

raw_data3=sqldf("select TimePoint,avg(SO2_value) as SO2_value,avg(NO2_value) as NO2_value,avg(O3_value) as O3_value,avg(CO_value) as CO_value,avg(PM10_value) as PM10_value,avg(PM2.5_value) as PM2.5_value from raw_data2group by TimePoint")

画图开始

6项污染物的时间序列图

plot.ts(subset(raw_data_1001,select=-TimePoint),col="red")

colnames(raw_data_1001)=c('date','SO2','NO2','O3','CO','PM10','PM2.5')

summaryPlot(raw_data_1001)

粗颗粒物浓度及差值变化

timeVariation(raw_data_1001,pollutantc("PM2.5","PM10"),difference = TRUE)

4 结论

(1)案例分析结果表明：2017年哈尔滨市大气污染物中PM2.5、NO2、PM10、O3为超标项目，其中PM2.5与PM10为主要污染物.秋冬季污染物浓度值高于春夏季，早晚出行时也出现了污染物浓度峰值，主要原因是取暖期对燃煤需求量大，且气温低，路面易结冰，机动车在外滞留时间长，导致污染物排放量增大.且哈尔滨市冬季受内蒙古-西伯利亚高压控制，下沉逆温易出现并时间长、范围广、势力强，不利于污染物扩散.

(2)R语言具有大数据处理的能力，可高效、快捷的处理和分析数据量庞大的大气污染数据.时间序列分析可直观的揭示污染物的时间变化规律以及预测变化趋势，为相关部分制定大气污染措施提供大数据支持.R语言用户也可根据分析需求，调用不同的函数编制自己的程序.