結(jié)合MLR和ARIMA模型的時空建模及預(yù)測

李 莎，林 暉

1.湖北第二師范學(xué)院 物理與機(jī)電工程學(xué)院，武漢430205

2.武漢數(shù)字工程研究所，武漢430205

隨著全球氣候變暖，氣候變化已成為全球熱點關(guān)注的問題之一。中國氣候變化趨勢與全球變化基本保持一致，1951—2009年中國全國年平均氣溫上升1.38 ℃，變暖速率達(dá)到0.23 ℃/10a[1]。近年來，國內(nèi)外有大量關(guān)于氣候變化分析和預(yù)測方面的研究。比如Hengl等[2]運用回歸克里金法（regression-kriging）對日均氣溫進(jìn)行時空預(yù)測，其中回歸處理過程使用觀測站點經(jīng)緯度、離海距離、海拔、時間、MODIS 地表氣溫圖像等信息對氣溫進(jìn)行回歸模擬。Meyer等[3]分別運用線性回歸法與隨機(jī)森林（RF）、廣義增強(qiáng)回歸模型（GBM）、Cubist[4]等3種機(jī)器學(xué)習(xí)算法，對南極洲地區(qū)的日均氣溫進(jìn)行空間分辨率為1 km 的插值預(yù)測。Kilibarda 等[4]針對全球陸地日均氣溫進(jìn)行分辨率為1 km 的時空插值研究，其中運用多元線性回歸模型（MLR）通過站點緯度和測量時間等協(xié)變量模擬氣溫的變化趨勢。Li 等[5]將時空MLR 與克里金模型結(jié)合對東北三省月均氣溫進(jìn)行插值研究，研究結(jié)果表明引入MLR 的插值精度要好于普通時空克里金。朱文博和朱連奇[1]運用AUNSPLINF方法對甘肅省50年的月均氣溫進(jìn)行插值處理，獲得空間分布圖并進(jìn)行趨勢特征分析。李均力等[6]以阿牙克庫木湖為研究對象，利用Corona、Landsat、HJ1A/1B 等形成的長時間序列遙感數(shù)據(jù)，分析其50年的年度和季節(jié)性變化特征，并結(jié)合降水、氣溫等數(shù)據(jù)分析與氣候變化的響應(yīng)。任婧宇等[7]采用距平、Mann-Kendall 趨勢檢驗和Sen’s 斜率估計方法分析在黃土高原地區(qū)未來時期四季氣候變化的時空分布特征。

ARIMA（差分自回歸移動平均）模型，即Box-Jenkins法，是一種常見的時間序列預(yù)測分析方法，廣泛應(yīng)用于醫(yī)療衛(wèi)生預(yù)測[8-9]、國民宏觀經(jīng)濟(jì)控制[10]、氣象預(yù)測[11]、交通運輸預(yù)測分析[12-13]等各個領(lǐng)域。傳統(tǒng)的氣象數(shù)據(jù)插值預(yù)測多為時間域的預(yù)測或空間域的插值（比如常見的克里金插值法[14-15]）。但是氣象變量具有典型的時空分布特點，比如地表氣溫的變化，不僅受時間的影響，還與氣象站點的地理位置、海拔高度等多種因素有關(guān)，因此，為了更好地擬合氣象變化規(guī)律，近年來時空多元線性回歸（MLR）模型被廣泛應(yīng)用于時空數(shù)據(jù)相關(guān)結(jié)構(gòu)的建模和分析。比如在氣象和環(huán)境等領(lǐng)域，時空MLR 模型與克里金插值法相結(jié)合稱作時空回歸克里金法，它被認(rèn)為是一種有效的時空插值方法而常被使用[5，16]。相比于普通克里金法，回歸克里金法可以大幅度提高解釋變量變化趨勢的比例[2]。此方法中，回歸建模實際上與克里金法是分離的。引入回歸建模，通過多個已知協(xié)變量模擬被研究變量的全局變化趨勢（空間變化趨勢、時間變化趨勢或者時空變化趨勢），繼而去除全局趨勢以保證數(shù)據(jù)的平穩(wěn)性[17-18]。為了進(jìn)一步描述氣溫與其他變量的相關(guān)性，提高時空預(yù)測精度，本文構(gòu)建了一個基于ARIMA模型和時空MLR 模型的預(yù)測方法：（1）通過時空MLR 建模，獲得地表氣溫關(guān)于時間、空間位置以及海拔等因素的時空變化趨勢，并分析和評價各因素對氣溫的影響程度；（2）運用ARIMA 模型對MLR 擬合殘差進(jìn)行時序擬合及預(yù)測，繪制預(yù)測曲線和誤差分布圖，并評價預(yù)測精度。

1 算法描述

1.1 時空MLR模型

在回歸分析中，如果有兩個或兩個以上的自變量，則稱為多元回歸，而如果自變量既包含空間變量又包含時間變量，則是時空多元回歸。事實上，一種現(xiàn)象常常是與多個因素相聯(lián)系的，由多個自變量的最優(yōu)組合共同來預(yù)測或估計因變量，比只用一個自變量進(jìn)行預(yù)測或估計更有效，更符合實際。氣溫就是一種常見的時空變量，同時用其空間和時間的分布特征來進(jìn)行預(yù)測，理論上應(yīng)該會具有更高的預(yù)測精度。

假設(shè)數(shù)據(jù)樣本為(Yi,X1i,…,Xki)，i=1,2,…,n，那么其多元回歸總體回歸函數(shù)為：

模型擬合的優(yōu)劣可通過判斷殘差的正態(tài)性、獨立性等來診斷，而模型中自變量的選取常采用逐步回歸法和全子集回歸法[19]。

1.2 ARIMA模型

常用的時間序列模型有自回歸模型AR(p)、移動平均模型MA(q)、自回歸移動平均模型ARMA(p,q)以及差分自回歸移動平均模型ARIMA(p,d,q)。其中，前3種常用于分析平穩(wěn)時間序列，而ARIMA通過差分可以用于處理非平穩(wěn)時間序列，因此ARIMA 更具有廣泛性和通用性。非平穩(wěn)時間序列是指包含趨勢性、季節(jié)性或周期性等特性的序列，它可能只含有其中的一種成分，也可能是幾種成分的組合[20]，且平穩(wěn)時間序列可以看作是非平穩(wěn)序列的一種特例。在消去局部水平或者趨勢之后，非平穩(wěn)時間序列顯示出一定的同質(zhì)性，經(jīng)過差分處理后可以轉(zhuǎn)換為平穩(wěn)時間序列，這樣的時間序列稱為齊次非平穩(wěn)時間序列，其中差分的次數(shù)就是齊次的階。

假設(shè)Yt是一個時間序列，那么其ARIMA(p,d,q)模型可表示為：

其中，p為自回歸階數(shù)，d為使之成為平穩(wěn)序列所做的差分次數(shù)（階數(shù)），q為移動平均階數(shù)，L是滯后算子，?i為自回歸部分的系數(shù)，θi為移動平均部分的系數(shù)，εi為殘差，是一個零均值白噪聲序列。

建立ARIMA模型的步驟通常如下：

（1）驗證序列的平穩(wěn)性。

（2）擬合模型，確定模型參數(shù)(p,d,q)。

（3）模型評價。一般來說，若模型合適，則模型的殘差應(yīng)滿足獨立正態(tài)分布；若不滿足，則返回至第（2）步重新確定模型參數(shù)。

（4）用模型進(jìn)行預(yù)測。

2 實驗數(shù)據(jù)

本文實驗數(shù)據(jù)由中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)提供，該數(shù)據(jù)包含若干氣象要素：風(fēng)速、氣溫和降水。選取1972 年1 月至2008 年12 月遼寧省月均氣溫數(shù)據(jù)為研究對象，應(yīng)用ARIMA 和MLR 模型進(jìn)行時空分析和預(yù)測。數(shù)據(jù)共涉及444 個月，27 個氣象觀測站點，站點分布如圖1 所示，其中4 個站點（新民、草河口、皮口、長海）的部分氣溫數(shù)據(jù)缺失，缺失情況如圖2 所示。存在缺失數(shù)據(jù)的站點，其站點的地理位置以及有效采樣的氣溫數(shù)據(jù)依然能夠為樣本的建模分析提供有效的時空分布信息，因此實驗中沒有刪除。實驗以1972—2005年共408 個月數(shù)據(jù)進(jìn)行建模，對2006—2008 年各月各站點進(jìn)行預(yù)測。

圖1 氣象觀測站點分布圖

圖2 數(shù)據(jù)缺失站點的月均氣溫

3 實驗分析

3.1 時序分解

氣溫的時空分布并不平穩(wěn)，尤其是時間分布上具有顯著的季節(jié)變化特點[21]（如圖3 所示）。ARIMA 模型雖然可以通過差分處理將非平穩(wěn)時間序列轉(zhuǎn)換為平穩(wěn)時間序列，但是為了優(yōu)化模型擬合效果，通常會預(yù)處理其中顯著影響平穩(wěn)的成分。而一個含有周期（季節(jié)）變化的時空變量通?？梢苑纸鉃槿齻€部分：時空趨勢項、時空季節(jié)項和時空隨機(jī)項。

圖3 沈陽站點月均氣溫時序分解

若將遼寧省各站點的月均氣溫觀測值定義為時空變量MAT={MAT(s,t)|s∈S,t∈T}，(s,t) 為時空域(S,T)內(nèi)的某一點，則：

其中，Sea(s,t)為時空季節(jié)項，Tr(s,t)為時空趨勢項，R(s,t)為時空隨機(jī)項。對應(yīng)于實驗數(shù)據(jù)的時間范圍，1972 年1 月記為t=1，2006 年1 月記為t=409。由于氣溫的變化周期為12 個月，則同一個站點的氣溫在不同年度相同月份的季節(jié)項一致，即：

式中的季節(jié)項Sea可以通過對每個站點時序分解進(jìn)行逐一提?。▓D3所示為沈陽站點月均氣溫的時序分解），而去除Sea之后的數(shù)據(jù)記為去季節(jié)項數(shù)據(jù)SRD。

3.2 時空趨勢擬合

除了顯著的季節(jié)變化之外，月均氣溫還與站點的地理位置、所處的海拔高度以及一年當(dāng)中所處月份有關(guān)，這種變化趨勢對應(yīng)的就是式（5）中的時空趨勢項Tr(s,t)，它可以通過MLR模型進(jìn)行擬合。在去季節(jié)項數(shù)據(jù)的基礎(chǔ)上，如果將Tr進(jìn)行提取分離，理論上剩下的數(shù)據(jù)應(yīng)該是平穩(wěn)的隨機(jī)項R(s,t)。

時空趨勢擬合可能涉及站點的東西坐標(biāo)x，南北坐標(biāo)y，海拔高度d，時間坐標(biāo)m等變量，本文通過全子集回歸法來確定變量的選擇。圖4所示，基于調(diào)整R平方（adjr2），同時含有截距項（INT）、d、y、x和m變量時，adjr2 達(dá)到最大值0.38，即此時的MLR 模型最佳，具體參數(shù)見表1。圖5 所示，站點的南北坐標(biāo)對SRD 的影響最大（權(quán)重達(dá)到45.5%左右），其次是東西坐標(biāo)（權(quán)重約24.7%）。這說明在消除季節(jié)性影響因素之后，氣溫受地理位置的影響較大。圖6所示，MLR擬合誤差的核密度曲線與正態(tài)曲線基本一致，這說明誤差很好的服從了正態(tài)分布，模型擬合效果良好。

圖4 SRD全子集回歸圖

圖5 MLR模型各變量相對權(quán)重點圖

圖6 MLR模型學(xué)生化殘差分布圖

表1 時空MLR模型參數(shù)表

3.3 ARIMA建模

進(jìn)行ARIMA 預(yù)測之前，通常要對時間序列進(jìn)行平穩(wěn)性檢驗。實驗依次對每個站點的R(s,t)時間序列進(jìn)行ADF 檢驗，均是平穩(wěn)的。以沈陽站點為例，其自相關(guān)圖和偏相關(guān)圖如圖7 所示，可見該站點R(s,t)是平穩(wěn)的，基于AIC 信息準(zhǔn)則，其ARIMA 模型參數(shù)分別為p=2，d=0，q=0，連同其他站點R(s,t)時間序列的ARIMA模型見表2。表2中顯示有多個站點的ARIMA模型參數(shù)相同，比如開原、清原等14 個站點的模型都為ARIMA（2，0，0），彰武等3 個站點的模型都為ARIMA（1，0，3），章黨等4 個站點的模型都為ARIMA（1，0，0），阜新和黑山站的模型為ARIMA（0，0，3）。由于各模型的自回歸系數(shù)（AR）、移動平均系數(shù)（MA）、截距（INT）等不同，所以各站點的ARIMA 模型實際互不相同，表3 是部分站點模型為ARIMA（2，0，0）的系數(shù)列表。

表2 各站點隨機(jī)項時間序列ARIMA模型

表3 ARIMA（2，0，0）模型系數(shù)

圖7 沈陽站隨機(jī)項ACF圖和PACF圖

3.4 預(yù)測結(jié)果分析

以27個站點前408個月的月均氣溫數(shù)據(jù)為訓(xùn)練集，依次進(jìn)行時序分解、MLR 擬合以及ARIMA 預(yù)測，將各階段值按照式（5）合成各站點2006—2008 年共36 個月的月均氣溫預(yù)測值，并與同時段原觀測值進(jìn)行比較，預(yù)測結(jié)果如圖8所示。圖8（a）中觀測值對預(yù)測值的點絕大多數(shù)落在或靠近y=x直線，且相關(guān)系數(shù)COR=0.993 4，非常接近于1；其時空預(yù)測誤差近似服從正態(tài)分布，如圖8（b）所示，約60%的誤差在±1℃以內(nèi)。進(jìn)一步引入統(tǒng)計指標(biāo)評價基于ARMIA 和MLR 模型的時空預(yù)測結(jié)果（如表4 所示），包括誤差中位數(shù)（MED）、誤差均值（ME）、誤差絕對值均值（MAE）以及均方根（RMSE），其中RMSE按式（7）得到，式中MAT*(si,ti)為月均氣溫的預(yù)測值，N為預(yù)測點總數(shù)。

圖8 時空預(yù)測整體結(jié)果

表4 時空預(yù)測統(tǒng)計指標(biāo)

由表4 可見，2006—2008 年的整體預(yù)測誤差MAE為0.927 7 ℃，RMSE 為1.151 0 ℃，說明時空預(yù)測整體效果良好。若按年份分別預(yù)測，MAE雖逐年增大，但相差不大，都在1 ℃以內(nèi)，而RMSE 均略高于1 ℃。預(yù)測誤差的空間分布以2006 年為例，圖9 分別顯示了2006年12 個月各站點的預(yù)測誤差分布圖，圖中氣泡大小對應(yīng)該站點的誤差絕對值大小，紅色代表預(yù)測值小于觀測值，綠色則相反。12個月中，大部分站點誤差較小，2、3、5、6、8、9、11、12月共計8個月的MAE小于1 ℃，其中5、8、9、11月MAE甚至小于0.5 ℃。另外，2、4、6、7月各站點的月均氣溫預(yù)測值普遍大于實際觀測值，而10 月的預(yù)測值普遍偏小，這主要是因為該月全省的月均氣溫較往年同期偏低或者偏高。如圖10所示，2006年4月遼寧省全部站點（其中26、27 號站點觀測值缺失）的月均氣溫值（紅色）明顯低于1972—2006 年4 月的月均氣溫平均值（藍(lán)色），因此2006年4月的各站點氣溫預(yù)測值普遍偏高；而2006年10月的情況恰恰相反，由于該月全省各站點月均氣溫（紅色）明顯高于35 年的平均值（藍(lán)色），導(dǎo)致當(dāng)月的時空預(yù)測值整體偏低。由此可見，當(dāng)某區(qū)域的氣溫出現(xiàn)較大幅度升降甚至出現(xiàn)極端天氣時，將對氣溫預(yù)測增加難度。

圖9 2006年12個月各站點時空預(yù)測誤差分布圖

圖10 氣溫觀測值與對應(yīng)氣溫均值對比圖

4 結(jié)語

區(qū)域氣溫是時空分布的變量，就某一站點而言，氣溫是具有顯著季節(jié)變化的時間序列，但就研究區(qū)域而言，它受地理位置、海拔高度等多種因素影響，且氣溫在時間和空間域的分布變化是交互的。本文將時空MLR模型和ARIMA模型應(yīng)用于遼寧省的月均氣溫預(yù)測分析中，前者擬合氣溫的時空變化趨勢，后者進(jìn)行隨機(jī)項預(yù)測，同時考慮了氣溫受空間、時間因素的影響。實驗以1972 年至2005 年的數(shù)據(jù)為訓(xùn)練集，對2006 年各月進(jìn)行整體預(yù)測，結(jié)果說明本文的方法能夠?qū)崿F(xiàn)較高精度的時空預(yù)測，同時針對少數(shù)預(yù)測誤差較大的點，給出了進(jìn)一步的原因分析。