羅 爻，劉 碩，高俁晗，張美琦，宣立強(qiáng)

(哈爾濱師范大學(xué))

0 引言

近年來中國快速城市化和工業(yè)化進(jìn)程，造成了嚴(yán)重的大氣污染.大氣污染引發(fā)的一系列環(huán)境問題，嚴(yán)重危害了人類的身心健康[1-2].隨著對(duì)環(huán)境和健康問題的重視，大氣污染已經(jīng)成為中國需解決的重大民生問題和國家發(fā)展問題[3-4].探討哈爾濱市大氣污染時(shí)空變化特征，有利于引導(dǎo)城市實(shí)施可持續(xù)發(fā)展戰(zhàn)略.

隨著大數(shù)據(jù)時(shí)代的到來，傳統(tǒng)的統(tǒng)計(jì)分析軟件難以滿足海量數(shù)據(jù)的處理需求[5].因此出現(xiàn)了諸多新的技術(shù)與工具，其中開源統(tǒng)計(jì)分析語言R被廣泛應(yīng)用[6].R語言是一種由統(tǒng)計(jì)學(xué)家開發(fā)的統(tǒng)計(jì)計(jì)算和繪圖的語言和環(huán)境，具有大數(shù)據(jù)處理的能力[7].為響應(yīng)國家要求哈爾濱市于2015年底前完成了空氣質(zhì)量監(jiān)測(cè)系統(tǒng)的建立，對(duì)重點(diǎn)污染企業(yè)實(shí)施嚴(yán)格的在線監(jiān)控.要實(shí)現(xiàn)大氣污染物網(wǎng)絡(luò)化治理，必須充分發(fā)揮大數(shù)據(jù)的特性[8].大氣污染監(jiān)測(cè)數(shù)據(jù)量大，傳統(tǒng)的數(shù)據(jù)分析工具耗時(shí)費(fèi)力，不便于高效、快捷的處理和分析數(shù)據(jù)，難以發(fā)揮監(jiān)測(cè)數(shù)據(jù)在污染防治中的作用，R語言及其眾多工具包為R語言用于大量大氣數(shù)據(jù)分析和可視化提供了強(qiáng)有力的支持.

早期的時(shí)序分析使用直觀數(shù)據(jù)來尋找規(guī)律，隨著研究進(jìn)展，發(fā)現(xiàn)簡(jiǎn)單的描述時(shí)序分析具有很大的局限性[9].用統(tǒng)計(jì)學(xué)原理來分析時(shí)間序列，可以更準(zhǔn)確的估計(jì)隨機(jī)序列的演變[10].該文應(yīng)用R語言及時(shí)間序列分析技術(shù)對(duì)哈爾濱市2017年全年監(jiān)測(cè)數(shù)據(jù)進(jìn)行各種可視化分析，以探討R在大氣數(shù)據(jù)分析領(lǐng)域應(yīng)用的巨大潛力.并為相關(guān)部門制定大氣污染控制措施提供科學(xué)依據(jù).

1 R語言與所需工具包(package)的安裝與功能

R語言官方主頁https://www.r-project.org/ ，R語言開發(fā)環(huán)境下載安裝地址：https://cran.r-project.org/mirrors.html，輔助R的工具RStudio，下載地址： http://www.rstudio.com/ide進(jìn)入下載頁面后，有Desktop和Server 2個(gè)版本，選擇Desktop.下載安裝R與RStudio后，使用install.packages()安裝所需工具包，library()加載所需工具包.所需工具包與功能見表1.

表1 工具包功能

2 數(shù)據(jù)處理

2.1 數(shù)據(jù)來源

該研究數(shù)據(jù)來源于黑龍江省生態(tài)環(huán)境廳(http://www.hljdep.gov.cn/)，通過Python編寫接口程序，自動(dòng)獲取哈爾濱市12個(gè)監(jiān)測(cè)點(diǎn)位CO、NO2、O3、PM2.5、PM10、SO2的小時(shí)數(shù)據(jù).該文應(yīng)用R語言以及相關(guān)工具包對(duì)哈爾濱市2017年12個(gè)測(cè)站點(diǎn)的小時(shí)數(shù)據(jù)進(jìn)行時(shí)間序列分析，取12個(gè)監(jiān)測(cè)站均值.以探討該軟件在空氣質(zhì)量數(shù)據(jù)分析領(lǐng)域應(yīng)用的巨大潛力[11].

2.2 數(shù)據(jù)處理

使用Excel表格對(duì)監(jiān)測(cè)數(shù)據(jù)進(jìn)行篩選排查，對(duì)數(shù)據(jù)進(jìn)行簡(jiǎn)單的分列處理，并導(dǎo)出為CSV格式文件以便R讀取.通過read.CSV()命令導(dǎo)入R中，利用TimePoint函數(shù)對(duì)檢測(cè)數(shù)據(jù)進(jìn)行時(shí)間格式轉(zhuǎn)化，sqldf包對(duì)檢測(cè)數(shù)據(jù)進(jìn)行數(shù)據(jù)匯總與計(jì)算平均值，運(yùn)用summaryPlot函數(shù)快速概覽數(shù)據(jù)整體情況，利用ggplot2函數(shù)的分面功能展示各污染物的時(shí)間變化，最終進(jìn)行數(shù)據(jù)分析.

3 時(shí)間序列分析技術(shù)

時(shí)間序列分析是大氣監(jiān)測(cè)數(shù)據(jù)的常用分析方法，包括對(duì)季度、月份、周-日、小時(shí)變化等特征進(jìn)行分析，以揭示污染物的時(shí)間變化規(guī)律以及預(yù)測(cè)變化趨勢(shì).時(shí)間序列分析的目的一般有兩個(gè)方面：一是認(rèn)識(shí)產(chǎn)生觀測(cè)序列的隨機(jī)機(jī)制，即建立數(shù)據(jù)生成模型；二是基于序列的歷史數(shù)據(jù)，也許還要考慮其他相關(guān)序列或因素，對(duì)序列未來的可能取值給出預(yù)測(cè)或預(yù)報(bào)[12].該研究中主要用到滑動(dòng)平均過程.一階滑動(dòng)平均過程公式：

Yt=et-θet-1:

E(Yt)=0

(1)

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

γk=ρk=0k≥2

方程式中θ和ρ1的一些數(shù)值可以幫助說明各種可能性.需要注意的是，負(fù)的θ對(duì)應(yīng)的ρ可以通過簡(jiǎn)單地取正的θ所對(duì)應(yīng)的ρ的負(fù)數(shù)得到.

3.1 時(shí)間序列分析

為了解哈爾濱市6項(xiàng)污染物的分布情況，對(duì)2017年監(jiān)測(cè)數(shù)據(jù)運(yùn)用summaryPlot函數(shù)快速概覽數(shù)據(jù)整體情況，繪制出圖1.對(duì)污染物分布情況進(jìn)行概述統(tǒng)計(jì)分析，并參照國家空氣質(zhì)量標(biāo)準(zhǔn)(GB3095 2012)[13]，對(duì)各項(xiàng)污染物超標(biāo)情況進(jìn)行描述.從圖1、表1中可看出，除SO2外，PM2.5、PM10、NO2年均值全都超過了國家空氣質(zhì)量2級(jí)標(biāo)準(zhǔn)，其中PM2.5超標(biāo)最嚴(yán)重.將數(shù)據(jù)集進(jìn)一步整合成適用于程序語言處理的形式(由寬變長)，利用ggplot2函數(shù)的分面功能展示各污染物的時(shí)間變化序列如圖2所示.

圖1 6項(xiàng)污染物的時(shí)間序列圖

表1 環(huán)境空氣污染物基本項(xiàng)目濃度限值

圖2 各污染物分列的時(shí)間序列變化

SO2、CO全年排放量不超標(biāo)，但秋冬季排放量高于春夏季.NO2排放量超標(biāo)率為7.8%，污染物濃度起伏不大，峰值出現(xiàn)在冬季.PM10與PM2.5污染物濃度起伏規(guī)律相似，在10月后污染物濃度達(dá)到峰值，秋冬污染物濃度值高于春夏.PM10、PM2.5超標(biāo)率分別為23.1%和69%，從圖2可看出，SO2、CO、PM2.5、PM10、NO2這幾項(xiàng)污染物濃度值都是秋冬季大于春夏季，其主要原因秋冬季為哈爾濱市的采暖期，氣溫O℃以下，早晚為燃煤供暖和出行的高峰期.寒冷天氣路面易結(jié)冰，機(jī)動(dòng)車行駛速率降低，在外滯留時(shí)間延長，污染物排放也變高[14].哈爾濱市冬季受內(nèi)蒙古-西伯利亞高壓控制，下沉逆溫易出現(xiàn)并時(shí)間長、范圍廣、勢(shì)力強(qiáng)，不利于污染物擴(kuò)散.且正值秋季農(nóng)忙結(jié)束，農(nóng)民大量燃燒秸稈等生物質(zhì)，導(dǎo)致空氣嚴(yán)重污染.O3呈鋸齒狀對(duì)稱分布，1～4月逐漸升高，在5月達(dá)到最大值，8月出現(xiàn)一個(gè)次峰，此后濃度值逐漸下降.

從圖2中可以看出PM10與PM2.5為哈爾濱市大氣污染主要成分，timeVariation函數(shù)提供了將不同量級(jí)的數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理的方法，該函數(shù)還提供了計(jì)算充值的功能，利用difference這一參數(shù)計(jì)算粗顆粒物濃度并展示其時(shí)間變化如圖3所示.

圖3可反應(yīng)出粗顆粒物的變化特征，月變化中顯示4、5月份濃度較高，說明春期揚(yáng)塵污染較嚴(yán)重，日變化中8：00-9：00與19：00-20：00出現(xiàn)峰值，說明道路交通揚(yáng)塵是該點(diǎn)位粗顆粒物的重要來源.周變化中周末平均小時(shí)濃度低于工作日濃度，表現(xiàn)出了明顯的周末效應(yīng).

3.2 所需代碼

讀取數(shù)據(jù)

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

時(shí)間格式轉(zhuǎn)化

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

計(jì)算站點(diǎn)的平均值作為最終畫圖的數(shù)據(jù)

圖3 PM10、PM2.5以及其插值的時(shí)間變化

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項(xiàng)污染物的時(shí)間序列圖

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 結(jié)論

(1)案例分析結(jié)果表明：2017年哈爾濱市大氣污染物中PM2.5、NO2、PM10、O3為超標(biāo)項(xiàng)目，其中PM2.5與PM10為主要污染物.秋冬季污染物濃度值高于春夏季，早晚出行時(shí)也出現(xiàn)了污染物濃度峰值，主要原因是取暖期對(duì)燃煤需求量大，且氣溫低，路面易結(jié)冰，機(jī)動(dòng)車在外滯留時(shí)間長，導(dǎo)致污染物排放量增大.且哈爾濱市冬季受內(nèi)蒙古-西伯利亞高壓控制，下沉逆溫易出現(xiàn)并時(shí)間長、范圍廣、勢(shì)力強(qiáng)，不利于污染物擴(kuò)散.

(2)R語言具有大數(shù)據(jù)處理的能力，可高效、快捷的處理和分析數(shù)據(jù)量龐大的大氣污染數(shù)據(jù).時(shí)間序列分析可直觀的揭示污染物的時(shí)間變化規(guī)律以及預(yù)測(cè)變化趨勢(shì)，為相關(guān)部分制定大氣污染措施提供大數(shù)據(jù)支持.R語言用戶也可根據(jù)分析需求，調(diào)用不同的函數(shù)編制自己的程序.