缺失數(shù)據(jù)下空氣質(zhì)量預(yù)測的ARIMA乘法季節(jié)模型

田文婷，宇世航

缺失數(shù)據(jù)下空氣質(zhì)量預(yù)測的ARIMA乘法季節(jié)模型

田文婷，宇世航

（齊齊哈爾大學(xué) 理學(xué)院，黑龍江 齊齊哈爾 161000）

以AQI指數(shù)作為評判空氣質(zhì)量狀況的標(biāo)準(zhǔn)進(jìn)行探究，對存在缺失值的AQI指數(shù)進(jìn)行填充，經(jīng)探究發(fā)現(xiàn)，填充后的AQI指數(shù)存在明顯季節(jié)特征，提出在原ARIMA模型上加入季節(jié)特征，建立ARIMA乘法季節(jié)模型。首先選東北四個(gè)城市進(jìn)行實(shí)證分析，用均值法和缺失森林法對缺失數(shù)據(jù)進(jìn)行填充，依據(jù)存在的季節(jié)特征分別建立適當(dāng)?shù)腁RIMA乘法季節(jié)模型，對未來12個(gè)月空氣質(zhì)量進(jìn)行預(yù)測并擬合未來變化趨勢，將預(yù)測的數(shù)據(jù)與真實(shí)數(shù)據(jù)進(jìn)行對比，得出ARIMA乘法季節(jié)模型在預(yù)測空氣質(zhì)量時(shí)優(yōu)于單純的ARIMA模型，使用缺失森林填充后建立的ARIMA乘法季節(jié)模型優(yōu)于用均值法填充。最后，確定基于缺失森林填充的缺失數(shù)據(jù)建立的ARIMA乘法季節(jié)模型為較優(yōu)的預(yù)測模型。

AQI指數(shù)；缺失森林；季節(jié)性；ARIMA乘法季節(jié)模型

目前，環(huán)境問題日益顯著，從國家長遠(yuǎn)發(fā)展來看，環(huán)境問題必須解決。近來，諸多城市出現(xiàn)霧霾天氣，空氣污染問題迫在眉睫，為解決環(huán)境問題，我國出臺了很多相應(yīng)政策。目前，有多位學(xué)者就環(huán)境問題、空氣質(zhì)量問題進(jìn)行了多方面的探究，以便提供合理的建議。余曉美[1]，孫銘[2]，孔朝莉[3]，王雪嬌[4]，王未來[5]等學(xué)者探究了空氣質(zhì)量的變化特征；嚴(yán)宙寧[6]，黃珊[7]等學(xué)者對空氣質(zhì)量進(jìn)行預(yù)測，肖悅[8]等對空氣質(zhì)量的空間分布特征進(jìn)行探究。先前學(xué)者均是基于數(shù)據(jù)完整的前提下提出的，但實(shí)際生活中，收集到的數(shù)據(jù)并非完整數(shù)據(jù)，可能存在缺失值。對于存在缺失值的數(shù)據(jù)，若直接刪除數(shù)據(jù)將會導(dǎo)致預(yù)測模型產(chǎn)生較大誤差，因此如何構(gòu)建預(yù)測模型至關(guān)重要。

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

本文數(shù)據(jù)主要來源于中華人民共和國生態(tài)環(huán)境部，部分?jǐn)?shù)據(jù)來源于http://www.tianqihoubao.com/。選取了2013年5月到2020年4月的大連市、沈陽市、哈爾濱市、長春市AQI指數(shù)的月度數(shù)據(jù)，對數(shù)據(jù)中存在的缺失值本文采取兩種處理方法。根據(jù)填充后的數(shù)據(jù)進(jìn)行分析，文章中進(jìn)行的數(shù)據(jù)分析均采用R語言[9-10]實(shí)現(xiàn)。

2 填補(bǔ)方法

本文對選取的數(shù)據(jù)集進(jìn)行缺失值處理，主要采用了均值填補(bǔ)法和缺失森林填充。均值填補(bǔ)是最常用的缺失值處理方法，其核心是用未缺失的數(shù)據(jù)均值代替其缺失值，因均值本身受異常值影響較大，使用該方法填充存在較大誤差。缺失森林[11]填補(bǔ)是利用隨機(jī)森林填充缺失值的非參數(shù)方法，其先是將缺失值從少到多進(jìn)行排序，用隨機(jī)森林回歸填補(bǔ)缺失值最少的，然后不斷進(jìn)行迭代即可得到填充后的數(shù)據(jù)，因缺失森林主要利用回歸進(jìn)行填補(bǔ)，所以該方法存在的缺點(diǎn)是可能存在過擬合。

3 模型介紹

4 實(shí)證分析

4.1 缺失值處理

為進(jìn)行建模，對選取的空氣質(zhì)量數(shù)據(jù)進(jìn)行缺失值處理，其缺失值情況如圖1所示。

圖1 四個(gè)城市數(shù)據(jù)缺失率

圖1中，大連市、哈爾濱市以及長春市的缺失率均為0.0138，而沈陽市的缺失率達(dá)到了0.0170。對缺失數(shù)據(jù)分別采用均值法和缺失森林進(jìn)行填補(bǔ)，后續(xù)建立ARIMA乘法季節(jié)均依據(jù)填補(bǔ)后的數(shù)據(jù)，下文中出現(xiàn)的繪制圖均依據(jù)于缺失森林填補(bǔ)的數(shù)據(jù)。

4.2 繪制各個(gè)城市的曲線圖

為探究AQI指數(shù)的變化特征，根據(jù)選取的大連市、沈陽市、哈爾濱市以及長春市2013年5月到2019年12月的AQI指數(shù)數(shù)據(jù)繪制趨勢圖如圖2所示。

圖2 AQI指數(shù)趨勢圖

如圖所示：（1）東北地區(qū)四個(gè)城市的AQI指數(shù)曲線呈現(xiàn)出季節(jié)變化特征。從2013年5月到2019年12月，每年呈現(xiàn)三個(gè)變化過程，第一個(gè)過程為下降階段，大約為每年一月到五月份，第二個(gè)過程為小幅波動過程，平穩(wěn)期大約為每年六月到九月左右，第三過程為急劇上升過程，大約為每年的十月到十二月。全年AQI的峰值基本出現(xiàn)在每年的冬季，出現(xiàn)這一現(xiàn)象的原因?yàn)闁|北地區(qū)依靠燃燒煤炭、秸稈等燃料進(jìn)行取暖，造成AQI指數(shù)增大，空氣污染嚴(yán)重。

（2）4個(gè)城市中AQI指數(shù)曲線存在差異。圖2中可知大連市的AQI指數(shù)曲線不同于其他城市，相比之下，大連市AQI指數(shù)在各個(gè)季節(jié)變化較小。相反，長春市在各個(gè)季節(jié)變化較大。出現(xiàn)該現(xiàn)象的原因與其城市所在的地理位置，經(jīng)濟(jì)發(fā)展水平及主要特征等因素有關(guān)。

4.3 平穩(wěn)性檢驗(yàn)及相關(guān)性檢驗(yàn)

4.3.1 平穩(wěn)性檢驗(yàn)

各個(gè)城市AQI指數(shù)數(shù)據(jù)呈現(xiàn)明顯周期性，且具有不顯著向下的趨勢性，為消除時(shí)間趨勢項(xiàng)，進(jìn)行一階差分得到新序列，對新序列進(jìn)行平穩(wěn)性檢驗(yàn)，檢驗(yàn)結(jié)果如表1所示。取顯著性水平為0.05時(shí)，ADF檢驗(yàn)的P值均小于0.05，說明拒絕原假設(shè)，接受備擇假設(shè)，即不存在單位根，認(rèn)為差分后序列是平穩(wěn)的，故本文中選取=1,=1, 構(gòu)造預(yù)測模型。

表1 ADF檢驗(yàn)表

4.3.2 相關(guān)性檢驗(yàn)

根據(jù)各個(gè)城市的AQI指數(shù)繪制其自相關(guān)圖3所示。

圖3 原始數(shù)據(jù)自相關(guān)圖

如圖3所示：各個(gè)城市的AQI指數(shù)存在明顯季節(jié)特征，根據(jù)自相關(guān)圖顯示，自相關(guān)函數(shù)之間存在較強(qiáng)的相關(guān)關(guān)系，故進(jìn)行季節(jié)差分，消除季節(jié)趨勢，使序列變?yōu)槠椒€(wěn)時(shí)間序列。

4.4 預(yù)測模型建立

4.4.1 大連市預(yù)測模型

根據(jù)大連市2013年5月到2019年12月的AQI指數(shù)數(shù)據(jù)建立預(yù)測模型，為建立較優(yōu)的模型，建立了多個(gè)模型對比得表2。

表2 大連市預(yù)測模型表

如表2所示：模型一為傳統(tǒng)的時(shí)間序列模型，其建立依據(jù)于均值處理缺失值數(shù)據(jù)，不含季節(jié)性，此時(shí)其AIC的值為635.01，當(dāng)將數(shù)據(jù)的季節(jié)性考慮進(jìn)去時(shí)，建立模型二，其AIC的值明顯減小，當(dāng)采用缺失森林填充數(shù)據(jù)時(shí)，建立了模型三時(shí)，其AIC的值變成552.02，略有減小，故此時(shí)選擇模型三。對模型三殘差序列繪制自相關(guān)圖和偏自相關(guān)圖，其并無顯著自相關(guān)，故該模型是充分的，模型的表達(dá)式為

將=12帶入上式即可得模型表達(dá)式，進(jìn)一步對該模型的殘差序列進(jìn)行Box-Ljung檢驗(yàn)，檢驗(yàn)結(jié)果中P值為0.9123大于顯著性水平0.05，故接受原假設(shè)，認(rèn)為殘差序列為白噪聲。與傳統(tǒng)模型相比，加入季節(jié)乘法模型后其AIC顯著減小。

4.4.2 沈陽市預(yù)測模型

根據(jù)沈陽市2013年5月到2019年12月的AQI指數(shù)數(shù)據(jù)建立預(yù)測模型，為建立較優(yōu)的模型，建立多個(gè)模型對比得表3。

表3 沈陽市預(yù)測模型表

如表3所示，模型一是傳統(tǒng)的模型預(yù)測方法，即不含季節(jié)乘法模型的ARIMA模型，其數(shù)據(jù)依據(jù)于均值填充的數(shù)據(jù)，模型的AIC為699.26，在此模型的基礎(chǔ)上加入季節(jié)乘法項(xiàng)，其AIC減小為592.04，AIC的值顯著減小，將原數(shù)據(jù)使用缺失森林填充，再依據(jù)該數(shù)據(jù)建立ARIMA乘法季節(jié)模型三，其AIC的值繼續(xù)減小。為檢驗(yàn)?zāi)Ｐ腿欠癯浞?，繪制殘差序列的自相關(guān)圖和偏自相關(guān)圖，其并無顯著關(guān)系中系數(shù)顯著，所以該模型是充分的，模型的表達(dá)式為

將=12帶入上述模型即得到預(yù)測模型，進(jìn)一步對該模型的殘差序列進(jìn)行Box-Ljung檢驗(yàn)，檢驗(yàn)結(jié)果中P值為0.9068大于顯著性水平0.05，故接受原假設(shè)，認(rèn)為殘差序列為白噪聲。

4.4.3 長春市預(yù)測模型

根據(jù)長春市2013年5月到2019年12月的AQI指數(shù)數(shù)據(jù)建立預(yù)測模型，為建立較優(yōu)的模型，建立多個(gè)模型對比得表4。

表4 長春市預(yù)測模型表

如表4所示：模型一使用了均值填充缺失值的數(shù)據(jù)，建立了不含季節(jié)項(xiàng)的(1,1,1)，其AIC的值為724.01，當(dāng)給模型一加入季節(jié)乘法項(xiàng)時(shí)，得到模型二，其AIC顯著降低，降低到602.99，進(jìn)一步探究，采用缺失森林填補(bǔ)缺失值時(shí)，建立了模型三，與前兩個(gè)模型相比，AIC的值降低，其值為602.96。在模型三中取顯著性水平為0.05時(shí)，所有系數(shù)均顯著，該模型的殘差序列的自相關(guān)圖和偏自相關(guān)圖都沒有顯著的自相關(guān)，所以該模型是充分，模型的表達(dá)式為

將=12帶入上述模型得到長春市預(yù)測模型，進(jìn)一步對該模型的殘差序列進(jìn)行Box-Ljung檢驗(yàn)，檢驗(yàn)結(jié)果中P值為0.7649大于顯著性水平0.05，故接受原假設(shè)，認(rèn)為殘差序列為白噪聲。

4.4.4 哈爾濱市預(yù)測模型

根據(jù)哈爾濱市2013年5月到2019年12月的AQI指數(shù)數(shù)據(jù)建立預(yù)測模型，為建立較優(yōu)的模型，建立多個(gè)模型對比得表5。

表5 哈爾濱市預(yù)測模型表

如表5所示：模型一為不含有季節(jié)效應(yīng)的模型，其AIC的值為777.51，在模型一的基礎(chǔ)上增加季節(jié)模型并進(jìn)行修正，得到模型二，模型二的AIC顯著降低。同理，將原始數(shù)據(jù)使用缺失森林進(jìn)行填充，然后建立模型三，繪制模型三的自相關(guān)和偏自相關(guān)圖可知，該模型的殘差序列的自相關(guān)圖和偏自相關(guān)圖不存在顯著的相關(guān)性所以該模型是充分，模型的表達(dá)式為

進(jìn)一步對該模型的殘差序列進(jìn)行Box-Ljung檢驗(yàn)，檢驗(yàn)結(jié)果中P值為0.6437大于顯著性水平0.05，故接受原假設(shè)，認(rèn)為殘差序列為白噪聲。

4.5 殘差序列相關(guān)圖檢驗(yàn)

根據(jù)上述模型中，繪制相應(yīng)的自相關(guān)圖與偏自相關(guān)圖如圖4所示。

圖4 殘差相關(guān)性圖

根據(jù)圖4各個(gè)模型殘差序列的相關(guān)圖顯示，殘差的自相關(guān)系數(shù)和偏自相關(guān)系數(shù)均落在2倍的標(biāo)準(zhǔn)差以內(nèi)，其無顯著相關(guān)性，故上述模型均是充分的，且其殘差序列為白噪聲序列。

4.6 繪制預(yù)測圖并得出預(yù)測數(shù)據(jù)

4.6.1 繪制預(yù)測圖

根據(jù)四個(gè)城市建立的預(yù)測模型分別預(yù)測未來12個(gè)月的AQI指數(shù)，得到預(yù)測圖如圖5所示。

圖5 預(yù)測趨勢圖

圖5中，沈陽市、哈爾濱市與長春市預(yù)測情況相似，各個(gè)城市AQI指數(shù)大致趨勢相同，根據(jù)預(yù)測圖可知，大連市未來一年的AQI指數(shù)出現(xiàn)波動，但整體呈下降趨勢，其他三個(gè)城市仍呈現(xiàn)先減后增的趨勢，符合對原有數(shù)據(jù)探索的規(guī)律。總而言之，其AQI指數(shù)在未來一年中有所降低。為檢驗(yàn)預(yù)測模型是否能較好地預(yù)測數(shù)據(jù)，對各個(gè)城市的預(yù)測數(shù)據(jù)提取四個(gè)月作為參照，將其與真實(shí)數(shù)據(jù)進(jìn)行對比。

4.6.2 預(yù)測數(shù)據(jù)與真實(shí)數(shù)據(jù)進(jìn)行比較

根據(jù)各個(gè)城市預(yù)測模型得到的數(shù)據(jù)與真實(shí)數(shù)據(jù)進(jìn)行對比得到各個(gè)城市的誤差表6。

表6 誤差表

如上述表所示：上述表為東北四個(gè)城市的AQI指數(shù)誤差表，根據(jù)誤差表顯示，各大城市的預(yù)測模型預(yù)測效果一般，存在預(yù)測誤差，但與傳統(tǒng)不含季節(jié)乘法模型的ARIMA模型相比，其預(yù)測效果較優(yōu)。

5 結(jié)論

研究結(jié)果表明，在針對含有缺失值的AQI指數(shù)數(shù)據(jù)時(shí)，采用缺失森林法填充數(shù)據(jù)要優(yōu)于采用均值法填充數(shù)據(jù)建立的模型，而針對AQI指數(shù)的預(yù)測模型中ARIMA乘法季節(jié)模型優(yōu)于ARIMA模型。對于東北這四個(gè)城市的AQI指數(shù)。該四個(gè)城市的AQI指數(shù)呈現(xiàn)季節(jié)變化特征，根據(jù)這一變化特征對四大城市分別建立了基于缺失森林填補(bǔ)缺失值數(shù)據(jù)下的預(yù)測模型，預(yù)測結(jié)果顯示各個(gè)城市預(yù)測效果一般，但與傳統(tǒng)不含季節(jié)乘法模型的ARIMA模型相比，其預(yù)測效果較優(yōu)。在實(shí)際生活中存在一些不可抗因素導(dǎo)致AQI指數(shù)的浮動較大，無法做到精準(zhǔn)預(yù)測，但該模型能比較好地揭示AQI值的發(fā)展規(guī)律，與未加入季節(jié)效應(yīng)的ARIMA模型相比，效果較為理想。

研究結(jié)果顯示，在對數(shù)據(jù)進(jìn)行隨機(jī)森林填充后建立的ARIMA季節(jié)乘法模型擁有比傳統(tǒng)模型較優(yōu)的預(yù)測效果，但仍存在誤差，其主要原因如下：

（1）各個(gè)城市空氣質(zhì)量變化趨勢相似，但具體變化幅度存在差異。根據(jù)研究結(jié)果表明，針對不同城市，其變化幅度不同，該模型相對而言對變化幅度較小的城市，預(yù)測效果較優(yōu)。

（2）數(shù)據(jù)受疫情影響較大。本文選取的預(yù)測的數(shù)據(jù)為2020年1月到4月的數(shù)據(jù)，但2020年初正值疫情爆發(fā)時(shí)期，突如其來的疫情使得各項(xiàng)生產(chǎn)活動不能如期進(jìn)行，而選擇用疫情之前的數(shù)據(jù)預(yù)測疫情期間的數(shù)據(jù)，使得預(yù)測存在一定誤差，現(xiàn)疫情已經(jīng)得到有效控制，后續(xù)會繼續(xù)收集近期數(shù)據(jù)進(jìn)行預(yù)測，可以有更高預(yù)測精度的預(yù)測結(jié)果。

[1] 余曉美，沈永昌. 中國環(huán)境保護(hù)重點(diǎn)城市空氣質(zhì)量的動態(tài)特征分析[J]. 統(tǒng)計(jì)與決策，2019, 35(11): 91-94.

[2] 孫銘，許夢婷，謝佩錦，等. 基于ARMA模型的哈爾濱市空氣質(zhì)量研究[J]. 哈爾濱師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2018, 34(04): 21-25.

[3] 孔朝莉，何文煒. 我國環(huán)保重點(diǎn)城市空氣質(zhì)量統(tǒng)計(jì)分析和評價(jià)[J]. 統(tǒng)計(jì)與管理，2020, 35(05): 28-36.

[4] 王雪嬌. 黑龍江省不同類型城市空氣質(zhì)量變化特征研究[D]. 哈爾濱：哈爾濱師范大學(xué)，2019.

[5] 王未來，李明艷，張美航，等. 我國31個(gè)重點(diǎn)城市2014-2018年空氣質(zhì)量變化特征分析[J]. 中國公共衛(wèi)生管理，2020, 36(04): 441-445.

[6] 嚴(yán)宙寧，牟敬鋒，趙星，等. 基于ARIMA模型的深圳市大氣PM_(2.5)濃度時(shí)間序列預(yù)測分析[J]. 現(xiàn)代預(yù)防醫(yī)學(xué)，2018, 45(02): 220-223, 242.

[7] 黃珊. 數(shù)據(jù)驅(qū)動模型下的鄭州大氣污染預(yù)測[D]. 鄭州：華北水利水電大學(xué)，2020.

[8] 肖悅. 中國空氣質(zhì)量時(shí)空分布特征及影響因素分析[D]. 重慶：西南大學(xué)，2018.

[9] 劉強(qiáng)，裴艷波，張貝貝. R語言與現(xiàn)代統(tǒng)計(jì)方法[M]. 北京：清華大學(xué)出版社，2016: 299-307.

[10] 李亞偉，劉玲，宋士勛，等. ARIMA乘法季節(jié)模型的R軟件實(shí)現(xiàn)[J]. 環(huán)境衛(wèi)生學(xué)雜志，2018, 8(04): 345-349.

[11] 劉佳星，張宏烈，劉艷菊，等. 基于缺失率的不完整數(shù)據(jù)填補(bǔ)算法[J]. 統(tǒng)計(jì)與決策，2021, 37(02): 39-41.

[12] 唐繼強(qiáng)，鐘鑫偉，劉健，等. 基于時(shí)間序列季節(jié)分類模型的軌道交通客流短期預(yù)測[J/OL]. 重慶交通大學(xué)學(xué)報(bào)（自然科學(xué)版）:1-9[2021-04-11].

Air quality prediction based on ARIMA multiplicative seasonal model

TIAN Wen-ting，YU Shi-hang

(College of Science, Qiqihar University, Heilongjiang Qiqihar 161006, China)

In this paper, we use AQI as a standard to evaluate the air quality, and fill in the missing AQI, it is found that the filled AQI has obvious seasonal characteristics, and seasonal features are added to the original Arima model to establish a multiplicative seasonal Arima model. Firstly, four cities in northeast China were selected for empirical analysis, and the missing data were filled in with the mean method and the missing forest method, and then the proper Arima multiplicative seasonal models were established according to their seasonal characteristics, the air quality in the next 12 months is predicted and the future trend is fitted. Finally, the forecast data is compared with the real data, and the Arima model is better than the Arima model in the forecast of air quality, the Arima multiplicative seasonal model with missing forest is better than that with mean value method. All in all, the Arima multiplicative seasonal model based on missing data of missing forest filling is determined to be a better prediction model.

AQI；missing forest；seasonality；arima multiplicative seasonal model

2021-09-11

黑龍江省自然科學(xué)基金（LH2019A027）；黑龍江省教育廳科學(xué)技術(shù)與研究項(xiàng)目（135309477）

田文婷（1998-），女，甘肅定西人，在讀碩士，主要從事時(shí)間序列應(yīng)用研究，1549358806@qq.com。

宇世航（1971-），女，黑龍江齊齊哈爾人，教授，博士，主要從事整值時(shí)間序列和復(fù)雜數(shù)據(jù)統(tǒng)計(jì)推斷研究，qqhrysh@163.com。

X831

A

1007-984X(2022)04-0083-07