武漢空氣質(zhì)量預(yù)測及受周邊城市污染影響分析

張學(xué)新 周泳岑

摘 要 利用武漢及其周邊城市長沙、南昌、合肥、襄陽、孝感等城市空氣質(zhì)量指數(shù)（AQI）及相關(guān)污染物數(shù)據(jù)，探索了武漢市空氣質(zhì)量指數(shù)的統(tǒng)計(jì)分布規(guī)律及空氣污染治理效果的評價(jià)問題，給出了能較好預(yù)測空氣質(zhì)量指數(shù)等級的推理規(guī)則.用統(tǒng)計(jì)模型分析武漢及周邊大中城市的SO2等空氣污染物之間的傳播及相互影響.

關(guān)鍵詞 環(huán)境經(jīng)濟(jì)學(xué);空氣污染影響機(jī)制;規(guī)則模型預(yù)測;空氣質(zhì)量;因果檢驗(yàn);協(xié)整關(guān)系

中圖分類號 X823 ???????????文獻(xiàn)標(biāo)識碼 A

Prediction of Air Quality in Wuhan

and Analysis of the Influence on its Air Pollution

Exerted by Those of Surrounding Cities

ZHANG Xuexin1， ZHOU? Yongcen 2

（1. School of Mathematics and Statistics， Hubei Engineering University， Xiaogan， Hubei 432000， China;

2. Wuhan Foreign Languages School， Wuhan， Hubei 430022， China）

Abstract The data of air quality index（AQI） and related pollutants in Wuhan and in its surrounding cities such as Changsha， Nanchang， Hefei， Xiangyang and Xiaogan are used in this paper.Then， the statistical distribution of air quality （AQI） in Wuhan and a scientific evaluation of the effect of pollution control and management are explored. More importantly， some inference rules are obtained which can well predict the air quality （AQI） level. Finally， by applying some statistical regression mode， the spread and interaction between air pollutants such as SO2 in atmosphere of Wuhan and its surrounding cities has been analyzed.

Key words environmental economics; impact mechanism of air pollution; rulebased model forecast; air quality; testing for causality; cointegration relation

1 引 言

對城市空氣質(zhì)量狀況及其氣象誘因，已有許多研究.時(shí)連俊等（2015）[1]對成都市空氣質(zhì)量狀況、李剛（2017）[2]對克拉瑪依市空氣質(zhì)量特征進(jìn)行描述性統(tǒng)計(jì)分析.韓霄和張美根（2014）[3] 通過模擬華北平原氣象場及主要?dú)馊苣z粒子的時(shí)空分布分析重霾成因，趙金霞等（2017）[4] 探討了天津?yàn)I海新區(qū)灰霾的主要?dú)庀笳T因.劉超等（2017）[5]討論上海冬夏兩季大氣污染特征及其污染來源，蘇維等（2017） [6] 討論南昌市PM2.5和PM10的時(shí)空變異特征.在預(yù)測方面，叢琳等（2017） [7]對北京市PM2.5做回歸分析，用PM10等關(guān)聯(lián)指標(biāo)對PM2.5作靜態(tài)預(yù)測.對武漢市空氣質(zhì)量研究，岳巖裕等（2016）[8]研究空氣質(zhì)量狀況與氣象條件的關(guān)系;劉慧君（2014）[9]分析PM2.5污染的成因;郭浩天（2014） [10]分析PM2.5中有機(jī)酸的分布特征及來源.這些研究側(cè)重于分析氣象因素對武漢市空氣質(zhì)量指數(shù)的影響，主要關(guān)注PM2.5污染物成分的分解.

空氣質(zhì)量指數(shù)是研究空氣質(zhì)量常用的一個(gè)指標(biāo).利用當(dāng)天全部實(shí)體污染物信息預(yù)測該天空氣質(zhì)量等級很有意義.第一，擬合武漢市PM2.5，PM10，SO2，CO，NO2和O3-8h等實(shí)體污染物濃度的分布函數(shù)，確認(rèn)選擇非時(shí)間序列分析方法預(yù)測質(zhì)量指數(shù)等級，使用決策樹模型進(jìn)行預(yù)測.第二，分析空氣質(zhì)量治理效果.第三，鑒于城市群之間的空間相關(guān)性，應(yīng)用Granger因果檢驗(yàn)、協(xié)整檢驗(yàn)，聯(lián)系地、整體地研究武漢市及周邊大中城市的污染物間的相互影響及傳播機(jī)制.

2 武漢市空氣質(zhì)量指數(shù)分布特征

在《空氣質(zhì)量歷史數(shù)據(jù)查詢》網(wǎng)頁里采集武漢市2013年12月初到2018年6月中旬的日相關(guān)記錄，其中有空氣質(zhì)量指數(shù)值、質(zhì)量等級、污染物PM2.5、PM10、SO2、CO、NO2、O3-8h的濃度觀測值.對少量缺失記錄，用最鄰近的3個(gè)觀測值的均值替補(bǔ).

對武漢市空氣質(zhì)量指數(shù)序列的頻率圖及周期性模型擬合分析得到，6年內(nèi)武漢市空氣質(zhì)量指數(shù)序列沒有周期性.再考察2015年對2014年、2016年對2015年、2017年對2016年的每日差分序列，均值、方差及其他分布特征均有顯著性差異，印證空氣質(zhì)量指數(shù)序列沒有周期性.統(tǒng)計(jì)分布結(jié)果從一個(gè)側(cè)面說明空氣質(zhì)量治理取得了一定的效果.其中，2013年12月份有21天、2014年有36天、2015年有20天、2016年有7天、2017年有8天、2018年1月～6月份有4天分別是重度污染，這些日期多數(shù)分布在12月份、1月份、2月份，少數(shù)分布在3～6月、10～11月.2014年有29天、2015年有32天、2016年有52天、2017年有56天、2018年1月～6月份有14天分別是優(yōu)秀天氣，優(yōu)秀日期多數(shù)分布在7月份、8月份、10月份、9月份、5月份，少數(shù)分布在2月份、11月份.

從空氣質(zhì)量指數(shù)的月度平均值看，武漢市空氣質(zhì)量的同比沒有明顯下降，如圖1所示.

圖1 武漢市5年月度空氣質(zhì)量指數(shù)均值比較

3 空氣質(zhì)量治理效果分析

評價(jià)空氣質(zhì)量治理狀況，一般是統(tǒng)計(jì)某個(gè)時(shí)段的“優(yōu)良”等級天數(shù)的累計(jì)頻率.為了分析某個(gè)時(shí)段連續(xù)的“優(yōu)良”等級天數(shù)，需使用序列模式方法.該方法的步驟是：1）保持空氣質(zhì)量指數(shù)序列觀測值的時(shí)間順序，對序列進(jìn)行分割聚類;2）在每個(gè)聚類里，對空氣質(zhì)量指數(shù)求平均值，按“優(yōu)”、“良”、“重度污染”、“中度污染”、“輕度污染”5級給出該類的屬性標(biāo)簽.3）對每個(gè)聚類，計(jì)算不同屬性標(biāo)簽的空氣質(zhì)量指數(shù)日數(shù)的比例.

序列模式分析不打亂空氣質(zhì)量指標(biāo)值的觀測時(shí)間順序，它依據(jù)樣本點(diǎn)之間的相近程度，將性質(zhì)相近的樣本點(diǎn)聚為一類.設(shè)某一類Gi={ti，ti+1，…，ti+j-1}，j≥1，表示Gi包含j個(gè)樣本點(diǎn){xti，xti+1，…，xti+j-1}.該類的均值為i.=∑ti+j-1l=tixl/j，該類的直徑為D（ti，ti+j-1）=∑ti+j-1l=ti（xl-i.）2.記L（p（n，k））是把n個(gè)樣本點(diǎn)分為k類遭受的損失函數(shù)，并定義為全體分類的直徑總和.當(dāng)n和k固定時(shí)，最小的L（p（n，k））意味著分割聚類的離差平方總和最小，因而分割聚類是有效的.可證

L（p（n，2））=min 2≤j≤nD（1，j-1）+D（j，n）， L（p（n，k））=min k≤j≤n{L（P（j-1，k-1））+D（j，n）}，（k≥3）.? （1）

由式（1）可知，尋找把n個(gè)樣本點(diǎn)分成k類（k≥3）的最優(yōu)分割，需要在對j-1個(gè)樣本點(diǎn)做k-1類最優(yōu)分割（2≤j≤n）的基礎(chǔ)上進(jìn)行.

不失一般性，僅對武漢市2014年至2017年間四個(gè)冬季的空氣質(zhì)量指數(shù)序列做最優(yōu)分割聚類，結(jié)果見表1.從表1看，2014-2017年4個(gè)冬季的空氣質(zhì)量，前3年持續(xù)好轉(zhuǎn)，但是在第4年情況惡化.2014年冬季，兩次出現(xiàn)重度污染，共9天，良好一次，44天.2015年冬季，兩次出現(xiàn)重度污染，共6天，良好一次，100天.2016年冬季，零次出現(xiàn)重度污染，良好二次，57天，其余全是輕度污染.2017年冬季，一次出現(xiàn)重度污染，共3天，良好一次，51天.應(yīng)用序列模式方法分析和表述武漢市大氣環(huán)境質(zhì)量，結(jié)果更具體、更深刻.

4 空氣質(zhì)量等級預(yù)測

設(shè)樣本集S的大小為|S|，樣本點(diǎn)共有m個(gè)不同的類別，其中屬于第i類Ci（i=1，2，…，m）的樣本點(diǎn)構(gòu)成集合Si，樣本點(diǎn)落入Ci的概率是pi.再設(shè)屬性A將S劃分為S=∪vj=1Sv，當(dāng)A取值aj時(shí)，落入類Ci的樣本點(diǎn)構(gòu)成集合Sij，發(fā)生的概率是pij=|Sij|/|S|，則屬性A的信息增益定義為 Gains（C，A）=（∑vj=1∑mk=1Skj/|S|）（∑mi=1pijlog 2pij）-∑mi=1pilog 2pi，信息增益率定義為

GainsR（C，A）=Gains（C，A）/[（∑vj=1∑mk=1Skj/|S|）（∑mi=1pijlog 2pij）].?? （2）

現(xiàn)在構(gòu)建分類與回歸樹（C&R tree）預(yù)測模型.這是一種基于樹的分類方法，它以遞歸劃分的方式將訓(xùn)練記錄分割為具有相似輸出變量值的若干個(gè)子集.C&R樹從根節(jié)點(diǎn)開始，每次利用信息增益率選擇一個(gè)當(dāng)前最佳的屬性進(jìn)行分枝，采用成本復(fù)雜性的修剪策略去控制樹的生長，最后產(chǎn)生基本形式是A→B的推理規(guī)則來預(yù)測新樣本點(diǎn)的類別.

為了檢驗(yàn)C&R樹模型的預(yù)測效果，把所有樣本按7：3的比例隨機(jī)分成兩個(gè)數(shù)據(jù)集，一個(gè)用于訓(xùn)練模型，一個(gè)用于測試模型的性能，結(jié)果如表2所示.

由C&R樹得到的推理規(guī)則有6條，預(yù)測準(zhǔn)確率在93%以上.變量的重要性依次是PM2.5，PM10，O3-8h，Co，NO2和SO2.

規(guī)則用于 優(yōu) - 包含 1 個(gè)規(guī)則： 如果 PM2.5 <= 75.5且 PM10 <= 46.5 則 優(yōu);

規(guī)則用于 良 - 包含 1 個(gè)規(guī)則：如果 PM2.5 <= 75.5且 PM10 > 46.5 且 O3-8h <= 160.5 則 良;

規(guī)則用于 輕度污染 - 包含 2 個(gè)規(guī)則：

規(guī)則 1如果 PM2.5 <= 75.5且 PM10 > 46.5且 O3-8h > 160.5則 輕度污染;

規(guī)則 2如果75.5 < PM2.5 <= 115.5則 輕度污染;

規(guī)則用于 中度污染 - 包含 1 個(gè)規(guī)則：如果 115.5 < PM2.5 <= 149.5則 中度污染;

規(guī)則用于 重度污染 - 包含 1 個(gè)規(guī)則：如果PM2.5 > 149.5則 重度污染;

缺?。毫?

給出一個(gè)上述推理規(guī)則的使用和驗(yàn)證實(shí)例.查閱2019年5月1日武漢空氣質(zhì)量指數(shù)日歷史數(shù)據(jù)（來源：https：//www.aqistudy.cn/historydata/daydata.php？city=武漢）知：AQI=79，質(zhì)量等級=良，PM2.5=28，PM10=71，SO2=8，CO=0.8，NO2=35，O3-8h =134.按照第二條規(guī)則，滿足條件PM2.5 =28<75.5且PM10=71>46.5且O3-8h =134<160.5，因此它被規(guī)則預(yù)測為良，的確如此.

5 武漢市及周邊大城市的空氣質(zhì)量相互影響

為了分析武漢市空氣質(zhì)量是否受周邊大城市長沙、南昌、合肥、襄陽、宜昌等地空氣質(zhì)量的影響，現(xiàn)對PM2.5、PM10、SO2、CO、NO2、O3-8h等實(shí)體污染物觀測值序列做統(tǒng)計(jì)建模分析.

5.1 長期穩(wěn)定關(guān)系

如果k個(gè)城市的某污染物濃度序列PC1t、PC2t，… PCkt不是平穩(wěn)的，但其d階差分是平穩(wěn)的，而且存在不全為零的常數(shù)c1，…，ck使得∑ki=1ciPCit是平穩(wěn)的，則稱這k個(gè)序列有協(xié)整關(guān)系（長期穩(wěn)定的關(guān)系）.把平穩(wěn)序列視為d=0的情形.對于具有協(xié)整關(guān)系PC*2=k0+k1PC*1的兩個(gè)序列，其動態(tài)關(guān)系用誤差修正模型給出：

ΔPC2t=α（PC2-k0-k1PC1）t-1+

β2ΔPC1t+εt，?? （3）

這里， 差分ΔPCt=PCt-PCt-1，假定殘差εt為白噪聲.

首先通過單位根檢驗(yàn)判斷武漢PM2.5W與長沙PM2.5C兩個(gè)時(shí)間序列之間的平穩(wěn)性.各單位根檢驗(yàn)的原假設(shè)都是：“序列沒有單位根”.細(xì)微差別是，Levin方法假設(shè)“兩個(gè)序列有共同的單位根”，Pesaran方法假設(shè)“每個(gè)序列各自有單位根”.在計(jì)算檢驗(yàn)的概率時(shí)，ADFFisher、Fisher、PPFisher方法使用漸近卡方分布，其它檢驗(yàn)方法則使用漸近正態(tài)分布.各個(gè)檢驗(yàn)結(jié)果都顯著性地拒接了原假設(shè)，說明兩個(gè)序列都是平穩(wěn)的，見表3.

其次，分析它們之間的協(xié)整關(guān)系.表4是武漢PM2.5W與長沙PM2.5C的誤差修正模型擬合結(jié)果.

表4對應(yīng)的誤差修正模型是：

ΔPM2.5Wt=-0.355（PM2.5W-16.225-0.85PM2.5C）t-1+0.507ΔPM2.5Ct+εt，??? （4）

它表明武漢的PM2.5W與長沙的PM2.5C有長期均衡關(guān)系，PM2.5W=16.225+0.85PM2.5C，同樣的，武漢與南昌、武漢與合肥的PM2.5都有長期均衡關(guān)系.值得注意的是，兩地之間的影響，不局限于滯后一天，經(jīng)初步分析，滯后可達(dá)7天.

對有關(guān)城市的其它污染物PM10、SO2、CO、NO2、O3-8h的時(shí)間序列分析可以得到類似的結(jié)論.

5.2 Granger因果關(guān)系

判斷一個(gè)城市的空氣污染是否是由另一個(gè)附近城市漂移過來的，即一個(gè)城市的空氣質(zhì)量指數(shù)是否與另一個(gè)城市的空氣質(zhì)量指數(shù)具有Granger因果關(guān)系，需要進(jìn)行相關(guān)的檢驗(yàn).把一個(gè)城市的某種污染物濃度PC1的歷史信息加入到關(guān)于另一個(gè)城市的同樣污染物PC2的預(yù)測模型中，如果該模型的預(yù)測能力得到顯著提高，就說PC1是PC2的“Granger原因”.

對模型

檢驗(yàn)H0：βj=0，（j=0，1，…l），其含義是假設(shè)“一地PM2.5不是另一地PM2.5的“Granger原因”.統(tǒng)計(jì)量F=（SSE0-SSE1）/mSSE1/（n-l-m-1）～F（m，n-l-m-1），這里n是觀測值個(gè)數(shù)，SSE1、SSE0是模型及零假設(shè)下的模型的殘差平方和.

5.3 武漢、長沙、南昌、合肥各PM2.5序列之間的Granger因果檢驗(yàn)

武漢、長沙、南昌、合肥各PM2.5序列PM2.5W、PM2.5C、PM2.5N、PM2.5H之間的Granger因果檢驗(yàn)結(jié)果見表5.在0.05、0.01的顯著性水平下，拒接“一地PM2.5不是另一地PM2.5的Granger原因”的假設(shè).即武漢、長沙、南昌、合肥四城市的PM2.5污染物是互為相互影響的，一地至少受另外一地的污染物提前1～2天的影響.

5.4 湖北省武漢市周邊城市大氣污染物的相互影響

對湖北省內(nèi)武漢市周邊大中城市空氣質(zhì)量狀況對武漢市的影響分析，僅選擇SO2污染物這一個(gè)指標(biāo)做實(shí)證分析.

武漢市某天的SO2污染物不僅受省內(nèi)其周邊大中城市當(dāng)天SO2污染物的影響，而且也受這些周邊大中城市前若干天SO2污染物的影響、還受武漢市自身前若干天SO2污染物的影響，因此，可以建立一個(gè)帶有分布滯后項(xiàng)的多元線性回歸模型.模型的概要見表6.

有些系數(shù)估計(jì)值是負(fù)的，比如-0.0805，意指前2天黃石市SO2污染物濃度的減少將引致武漢市今天SO2污染物濃度0.0805個(gè)單位的增加，或者說前2天黃石市SO2污染物濃度的增加將引致武漢市今天SO2污染物濃度0.0805個(gè)單位的減少.說明武漢市及周邊SO2污染物有關(guān)聯(lián)，污染途徑可能取決于風(fēng)向（因?yàn)闆]有取得風(fēng)向的歷史數(shù)據(jù)，不能確定其具體關(guān)系），它引致了SO2污染物的流動.對宜昌沒有檢測到這種關(guān)系.

Durbin統(tǒng)計(jì)量的值是2.077950，說明殘差序列非自相關(guān)，圖2顯示殘差分布近似于正態(tài)分布.模型總體效果很好.

類似的，對其它污染物濃度的分析也表明武漢市及省內(nèi)部分周邊大中城市的污染物濃度是相互影響的.

6 結(jié)論與討論

武漢市空氣質(zhì)量指數(shù)的統(tǒng)計(jì)分布難以得到，使得當(dāng)前評價(jià)空氣質(zhì)量治理效果的方法比較簡單，序列模式方法嘗試更加科學(xué)化的評價(jià)，但是帶來一定的時(shí)間復(fù)雜度.規(guī)則模型不僅可準(zhǔn)確地預(yù)測空氣質(zhì)量等級，而且能給出等級對應(yīng)的各種污染物濃度的臨界值，給人具體的空氣質(zhì)量認(rèn)識.

武漢市與周邊大中城市空氣質(zhì)量污染物之間的影響機(jī)制很復(fù)雜.Granger因果檢驗(yàn)及其它統(tǒng)計(jì)分析表明這種影響機(jī)制是相互的，武漢市空氣污染物與長沙、南昌、合肥的空氣污染物有著長期均衡關(guān)系，可以有比較穩(wěn)定的2～7天的前期影響.武漢當(dāng)天空氣污染物還受孝感、黃石及武漢市自身前1～3天污染物影響.建立帶有分布滯后項(xiàng)的多元線性回歸可以定量地描述武漢市與周邊大中城市空氣質(zhì)量污染物之間的影響力度.模型（6）只用PM2.5、PM10、SO2、CO、NO2、O3-8h污染物的濃度就較好地預(yù)測了武漢市未來一天SO2污染物濃度.如果能加入其它氣象因素，將極大提升模型的預(yù)測能力.

參考文獻(xiàn)

[1] 時(shí)連俊，徐建，王變芳，等.成都市空氣質(zhì)量狀況研究[J].資源與環(huán)境，2015，31（8）： 986-989.

[2] 李剛.克拉瑪依市空氣質(zhì)量特征分析[J].干旱環(huán)境監(jiān)測，2017，31（2）： 75-79.

[3] 韓霄，張美根.2013年1月華北平原重霾成因模擬分析[J].氣候與環(huán)境研究，2014，19（2）：127-139.

[4] 趙金霞，沈岳峰， 范蘇丹.天津市濱海新區(qū)持續(xù)性重度霧霾成因分析[J].沙漠與綠洲氣象，2017， 11（6）： 69-74.

[5] 劉超，花叢，康志明. 2014-2015年上海地區(qū)冬夏季大氣污染特征及其污染源分析[J].氣象，2017，43（7）：823-830.

[6] 蘇維，賴新云，賴勝男，等.南昌市城市空氣PM2.5和PM10時(shí)空變異特征及其與景觀格局的關(guān)系[J].環(huán)境科學(xué)學(xué)報(bào).2017，37（7）： 72432-2439.

[7] 叢琳，孫德山，鄒存利，等.北京市PM2.5的相關(guān)因素研究[J].經(jīng)濟(jì)數(shù)學(xué)，2017，34（4）： 26-29.

[8] 岳巖裕，王曉玲，張蒙晰，等.武漢市空氣質(zhì)量狀況與氣象條件的關(guān)系[J].暴雨災(zāi)害，2016，35（3）：271-278.

[9] 劉慧君.武漢市PM2.5污染的演變預(yù)測及成因分析和仿真[D].長沙：湖南大學(xué)數(shù)學(xué)與計(jì)量經(jīng)濟(jì)學(xué)院，2014.

[10]郭浩天.武漢市大氣PM2.5中有機(jī)酸的時(shí)空分布特征及來源解析[D].武漢：武漢理工大學(xué)資源與環(huán)境工程學(xué)院，2014.