基于非參數(shù)檢驗和ARIMA模型的蕪湖市空氣質(zhì)量指數(shù)研究

余婉風(fēng) 呂科 劉洋 朱偉杰

摘要：在近幾年全國空氣質(zhì)量總體有所好轉(zhuǎn)的大環(huán)境下，本文通過分析蕪湖市空氣質(zhì)量指數(shù)數(shù)據(jù)，探究蕪湖市空氣質(zhì)量現(xiàn)狀，并構(gòu)建AQI短期預(yù)測模型，為蕪湖市政府控制環(huán)境污染和有效地治理提供科學(xué)的依據(jù)。2013年12月1日—2020年10月31日近8年蕪湖市空氣質(zhì)量指數(shù)（AQI）數(shù)據(jù)作為研究對象，R語言為實現(xiàn)工具。首先，分析AQI數(shù)據(jù)曲線圖，采用非參數(shù)檢驗Kruskal-Wallis檢驗比較這8年AQI數(shù)據(jù)是否具有顯著性差異;其次，根據(jù)對AQI時間序列平穩(wěn)性分析結(jié)果，選擇合理的時間序列模型—ARIMA模型，估計模型參數(shù)，建立擬合模型，并評價模型有效性;最后，利用模型預(yù)測未來幾個月AQI。

關(guān)鍵詞：AQI;非參數(shù)檢驗;時間序列分析;ARIMA預(yù)測模型; R語言

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

文章編號：1009-3044（2021）11-0239-03

Study on Air Quality Index（AQI） of? Wuhu City Based on Nonparametric Test and ARIMA Model

YU Wan-feng， LV Ke， LIU Yang， ZHU Wei-jie

（College of Big Data and Artificial Intelligence， Anhui Institute of Information Technology， Wuhu 241003，China）

Abstract： In recent years， air conditions nationwide has been improving， through analyzing the air quality index data of Wuhu city， to explore the current situation of air conditions of Wuhu City， and model for AQI prediction ， which provides scientific basis for Wuhu city government to control environmental pollution effectively. From December 1， 2013 to October 31， 2020， air quality Index （AQI） data of Wuhu city in recent 8 years were taken as the research object， and R language was used as the implementation tool. Firstly， the AQI data graph was analyzed， and the kruskal-Wallis test was used to compare the significant differences of the AQI data over the past 8 years. Secondly， according to the stationary analysis of AQI time series， an reasonable time series model -- ARIMA model is selected， to estimate fitting model parameters and evaluate the effectiveness of the model. Finally， the ARIMA model is used to predict AQI in the coming months.

Key words： AQI; nonparametric tests; Time-Series analysis; ARIMA prediction model; R softwre

隨著城市工業(yè)的成長壯大，空氣污染嚴(yán)重，空氣質(zhì)量惡化，不僅影響到人們的正常生活，而且威脅著人們的身心健康[1]。環(huán)保作為國家戰(zhàn)略性重點產(chǎn)業(yè)，全國各級政府對本省市的環(huán)保工作高度重視，紛紛積極推動節(jié)能減排和環(huán)境治理工作，截至目前，全國幾乎所有的省市，均已出臺生態(tài)保護(hù)相關(guān)政策、資金支持或項目管理方案，為我國全面推進(jìn)環(huán)保事業(yè)提供有力的支持[2-3]。蕪湖市政府以科學(xué)發(fā)展觀為指導(dǎo)，緊緊圍繞國家環(huán)?？偩痔岢龅钠邆€方面整改要求，開展集中整治行動，推進(jìn)環(huán)保執(zhí)法，推進(jìn)節(jié)能減排，推進(jìn)全市經(jīng)濟(jì)發(fā)展模式的轉(zhuǎn)變[4]。

同時，隨著R語言中各類工具包（package）的不斷涌現(xiàn)，數(shù)理統(tǒng)計分析工具和可視化方案的R程序包得以開發(fā)應(yīng)用開來。本研究數(shù)據(jù)來自2013年-2020年日報AQI數(shù)據(jù)，非參數(shù)檢驗Kruskal-Wallis檢驗可以對多組獨立樣本的多重比較，安裝、引用R語言pgirmess程序包，調(diào)用kruskalmc函數(shù)實現(xiàn)Kruskal-Wallis檢驗[5];利用R語言軟件包tseries進(jìn)行時間序列分析;調(diào)用軟件包forecast的auto.arima（）函數(shù)進(jìn)行系統(tǒng)定階，擬合最優(yōu)ARIMA模型。本文中，運(yùn)用非參數(shù)檢驗對比近幾年蕪湖市空氣質(zhì)量指數(shù)有無明顯改善;應(yīng)用時間序列分析建立預(yù)測模型，對模型進(jìn)行參數(shù)估計，診斷和評價，確定最優(yōu)模型，并檢驗?zāi)Ｐ陀行訹6]，最后運(yùn)用模型進(jìn)行預(yù)測。

1 非參數(shù)檢驗—Kruskal-Wallis檢驗

1.1抽樣數(shù)據(jù)

研究對象是蕪湖市2013年12月到2020年7月的空氣質(zhì)量指數(shù)（AQI），取每個月AQI平均值。

為了更直觀地分析2013年到2020年AQI的變化趨勢，描繪數(shù)據(jù)趨勢圖。如圖1所示，2013年12月到2020年10月蕪湖市空氣質(zhì)量指數(shù)（AQI）有逐年降低趨勢，但不明顯，且可能呈周期性變化。僅從圖1無法確定AQI是否有明顯的逐年降低趨勢，且考慮到AQI分布未知，樣本量少，采用非參數(shù)檢驗進(jìn)行8組獨立樣本的多重比較，判斷不同年份的AQI數(shù)據(jù)是否具有統(tǒng)計學(xué)差異。

1.2 R語言實現(xiàn)Kruskal-Wallis秩和檢驗

利用抽樣技術(shù)，從2013年12月到2020年10月中抽取每個月1號的數(shù)據(jù)，即83個樣本數(shù)據(jù)，不同年份作為組別，共有8個獨立組別。安裝并加載程序包pgirmess，調(diào)用Kruskal-Wallis秩和檢驗函數(shù)計算結(jié)果p-value=0.02139，顯著性水平α=0.05，p-value小于顯著性水平，說明拒絕原假設(shè)（原假設(shè)是8組樣本數(shù)據(jù)無顯著性差異），即這8年的AQI數(shù)據(jù)有顯著性差異。結(jié)合圖1分析，可以認(rèn)為近幾年蕪湖市的空氣質(zhì)量指數(shù)有明顯的降低趨勢，空氣質(zhì)量明顯好轉(zhuǎn)，環(huán)境治理效果顯著。為了給蕪湖市政府控制空氣污染和有效地治理提供科學(xué)的依據(jù)，分析數(shù)據(jù)特征，建立有效的AQI預(yù)測模型。

2 AQI時間序列預(yù)測模型

從2013年1月到2020年10月蕪湖市空氣質(zhì)量指數(shù)（AQI）總趨勢雖無明顯降低，但AQI不是雜亂無章的，如圖2所示，分解了數(shù)據(jù)的變化趨勢、季節(jié)性和不確定性因素，分析圖中“seasonal”曲線， AQI有明顯的季節(jié)性、周期性趨勢?；贏QI的數(shù)據(jù)特征，采用時間序列分析建立AQI短期預(yù)測模型，在預(yù)測模型有效的前提下利用模型預(yù)測2020年未來幾個月的空氣質(zhì)量指數(shù)，并繪制擬合圖形。

2.1 平穩(wěn)性分析

（1）平穩(wěn)性分析

常見的時間序列模型包括ARAM和ARIMA模型等，根據(jù)時間序列的平穩(wěn)性來選擇合適的預(yù)測模型，常見的時間序列平穩(wěn)性檢驗方法有PP檢驗法和ADF檢驗法。

R軟件安裝加載tseries包，進(jìn)行PP檢驗和ADF檢驗。PP檢驗結(jié)果中出現(xiàn)警告信息，ADF檢驗p值=0.04803< 0.05，默認(rèn)顯著性水平[α=0.05]，拒絕原假設(shè)（原假設(shè)認(rèn)為時間序列是非平穩(wěn)的），可以認(rèn)為AQI數(shù)據(jù)是平穩(wěn)的時間序列。

（2）白噪聲檢驗

對于平穩(wěn)的時間序列需要進(jìn)行白噪聲檢驗，因為白噪聲是純隨機(jī)序列，對純隨機(jī)序列建模毫無意義。對AQI時間序列白噪聲檢驗結(jié)果p值=4.6637e-07小于顯著性水平0.05，拒絕原假設(shè)，可以認(rèn)為AQI數(shù)據(jù)為非白噪聲時間序列，對該平穩(wěn)時間序列建模有意義。

考慮到AQI呈季節(jié)性、周期性變化，選擇目前最常用的擬合平穩(wěn)序列的模型ARAM（p， q）模型（自回歸移動平均模型）作為AQI時間序列預(yù)測模型。ARMA模型（自回歸移動平均模型）又可以細(xì)分為AR（p）模型、MA（q）模型和ARMA（p， q）模型三大類，確定AQI預(yù)測模型屬于哪一類這都取決于階數(shù)p和q的值，其一般準(zhǔn)則如表1所示。

2.2 確定ARAM（p， q）模型階數(shù)p和q

對于ARAM（p， q）模型，一般通過分析ACF自相關(guān)圖和PACF偏相關(guān)圖來估計兩個未知參數(shù)p和q。在R軟件中，擬合線性回歸模型，再調(diào)用acf（）和pacf（）兩個函數(shù)繪制ACF自相關(guān)圖和PACF偏相關(guān)圖。

圖3中，自相關(guān)系數(shù)拖尾，且從1開始控制在置信區(qū)間之內(nèi)，p= 1;圖4中，偏相關(guān)系數(shù)拖尾，且從2開始控制在置信區(qū)間之內(nèi)，q= 2。根據(jù)表1確定ARAM模型準(zhǔn)則，初步推測AQI預(yù)測最優(yōu)模型為ARMA（1， 2）。

2.3 系統(tǒng)自動定階

為了驗證通過時間序列穩(wěn)定性、ACF自相關(guān)圖和PACF偏相關(guān)圖分析確定的ARAM（1， 2）模型是否比較合理，利用R語言軟件包forecast的ARIMA模型函數(shù)，默認(rèn)預(yù)測模型為ARIMA（p， d， q）模型（差分整合移動平均自回歸模型）進(jìn)行系統(tǒng)自動定階，與通過平穩(wěn)時間序列ACF自相關(guān)圖和PACF偏相關(guān)圖分析確定的ARAM（1， 2）模型進(jìn)行比較，找到最優(yōu)模型。安裝、加載forecast包，時間序列ARIMA建模，結(jié)果如表2，ARIMA（1，0，0）（2，0，0）12 是擬合AQI時間序列的最佳模型，其中參數(shù)p= 1，d= 0，q= 2。參數(shù)d表示差分階數(shù)，當(dāng)d為0時，ARIMA模型就等同于ARMA模型，即ARIMA（1， 0， 2）模型與ARMA（1， 2）模型是等價的。由此，通過分析自相關(guān)系數(shù)圖和偏相關(guān)系數(shù)圖確定的最優(yōu)模型與R語言進(jìn)行系統(tǒng)自動定階結(jié)果一致。

在確定了最優(yōu)模型ARIMA（1，0，2）模型之后，并采用極大似然估計思想進(jìn)行參數(shù)估計，如表2，ARIMA（1，0，2）擬合模型：

[yt=0.6101yt-1+0.3725εt-1+0.2043εt-2+84.4825]

接下來，對ARIMA（1，0，2）模型的3個系數(shù)和1個截距進(jìn)行假設(shè)檢驗，用極大似然估計思想估計的系數(shù)的絕對值除以其標(biāo)準(zhǔn)差（s.e.）得到的商與t檢驗5%的臨界值1.96比較，商的絕對值大于1.96，拒絕原假設(shè)（原假設(shè)為參數(shù)影響不顯著），否則認(rèn)為參數(shù)影響顯著。t（ar1）= 6.532> 1.96， t（sar1）= 3.242> 1.96，t（sar2）= 1.625< 1.96， t（mean）= 9.621> 1.96，所以參數(shù)ar1、sar1和截距對該時間序列模型都具有顯著影響。

2.4 預(yù)測模型有效性檢驗

常用的檢驗ARIMA模型有效性的方法是純隨機(jī)性檢驗方法Box-Ljung檢驗，R語言提供了Box-Ljung檢驗的tsdiag（）函數(shù)。調(diào)用tsdiag（）函數(shù)，檢驗結(jié)果如圖5所示。

第二行的ACF檢驗說明殘差沒有明顯的自相關(guān)性;第三行的Box-Ljung檢驗顯示所有的p值都大于顯著性水平0.01，殘差序列不能拒絕純隨機(jī)的原假設(shè)，說明殘差為白噪聲，所以此預(yù)測模型是有效模型。利用ARIMA模型預(yù)測過去8年的AQI數(shù)據(jù)，將預(yù)測值與真實值對比，如表3所示，該預(yù)測模型在誤差允許的范圍內(nèi)能反映AQI數(shù)據(jù)的變化趨勢和規(guī)律，有參考價值。

3 ARIMA模型預(yù)測未來空氣質(zhì)量指數(shù)

在檢驗ARIMA（1， 0， 2） 預(yù)測模型是有效的前提下，利用預(yù)測2020年未來幾個月的空氣質(zhì)量。調(diào)用forecast（）函數(shù)繪制預(yù)測模型曲線，如圖6所示，兩條曲線分別表示AQI時間序列和ARIMA（1， 0， 2）模型，右邊粗線表示預(yù)測模型預(yù)測的未來幾個月AQI序列。

調(diào)用forecast（）函數(shù)預(yù)測11月和12月每月平均的空氣質(zhì)量指數(shù)，有：

[predict（AQI—11）=68.67];

[predict（AQI—12）=78.73]。

4 總結(jié)

時間序列模型的缺點是精確度不夠高，為了提高精確度，可以在模型中增加與空氣質(zhì)量相關(guān)的變量如PM2.5、CO和N02含量等相關(guān)變量。考慮到目前蕪湖市政府正在逐漸加大環(huán)保力度，從長期來看，由于外部因素干預(yù)加強(qiáng)，導(dǎo)致模型的預(yù)測能力下降，所以此模型較適合短期內(nèi)預(yù)測。

參考文獻(xiàn)：

[1] 史美義.淺析當(dāng)前我國城市環(huán)境污染的現(xiàn)狀及原因[J].科技信息，2012（18）：79.

[2] 國務(wù)院.國務(wù)院關(guān)于落實科學(xué)發(fā)展觀加強(qiáng)環(huán)境保護(hù)的決定[J].中國環(huán)境監(jiān)測，2006，22（1）：1-6.

[3] 安徽省人民政府.安徽省人民政府貫徹國務(wù)院關(guān)于落實科學(xué)發(fā)展觀加強(qiáng)環(huán)境保護(hù)決定的實施意見[J].安徽省人民政府公報，2006（21）：21-26.

[4] 付偉，司紅君，盧堯，等.蕪湖市空氣質(zhì)量特征及其受氣象要素的影響分析[C]//第35屆中國氣象學(xué)會年會論文集.合肥，2018：128-129.

[5] 金英良，趙華碩，孫桂香，等.基于R軟件的多組獨立樣本秩和檢驗的多重比較[J].預(yù)防醫(yī)學(xué)論壇，2016，22（11）：805-806，809.

[6] 牟敬鋒，趙星，樊靜潔，等.基于ARIMA模型的深圳市空氣質(zhì)量指數(shù)時間序列預(yù)測研究[J].環(huán)境衛(wèi)生學(xué)雜志，2017，7（2）：102-107，117.

【通聯(lián)編輯：唐一東】