結合MLR和ARIMA模型的時空建模及預測

李 莎，林 暉

1.湖北第二師范學院 物理與機電工程學院，武漢430205

2.武漢數字工程研究所，武漢430205

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

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

1 算法描述

1.1 時空MLR模型

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

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

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

1.2 ARIMA模型

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

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

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

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

（1）驗證序列的平穩性。

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

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

（4）用模型進行預測。

2 實驗數據

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

圖1 氣象觀測站點分布圖

圖2 數據缺失站點的月均氣溫

3 實驗分析

3.1 時序分解

氣溫的時空分布并不平穩，尤其是時間分布上具有顯著的季節變化特點[21]（如圖3 所示）。ARIMA 模型雖然可以通過差分處理將非平穩時間序列轉換為平穩時間序列，但是為了優化模型擬合效果，通常會預處理其中顯著影響平穩的成分。而一個含有周期（季節）變化的時空變量通常可以分解為三個部分：時空趨勢項、時空季節項和時空隨機項。

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

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

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

式中的季節項Sea可以通過對每個站點時序分解進行逐一提取（圖3所示為沈陽站點月均氣溫的時序分解），而去除Sea之后的數據記為去季節項數據SRD。

3.2 時空趨勢擬合

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

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

圖4 SRD全子集回歸圖

圖5 MLR模型各變量相對權重點圖

圖6 MLR模型學生化殘差分布圖

表1 時空MLR模型參數表

3.3 ARIMA建模

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

表2 各站點隨機項時間序列ARIMA模型

表3 ARIMA（2，0，0）模型系數

圖7 沈陽站隨機項ACF圖和PACF圖

3.4 預測結果分析

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

圖8 時空預測整體結果

表4 時空預測統計指標

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

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

圖10 氣溫觀測值與對應氣溫均值對比圖

4 結語

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