李 莎，林 暉

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

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

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

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

1 算法描述

1.1 時(shí)空MLR模型

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

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

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

1.2 ARIMA模型

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)數(shù)據(jù)

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

圖1 氣象觀測(cè)站點(diǎn)分布圖

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

3 實(shí)驗(yàn)分析

3.1 時(shí)序分解

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

圖3 沈陽(yáng)站點(diǎn)月均氣溫時(shí)序分解

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

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

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

3.2 時(shí)空趨勢(shì)擬合

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

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

圖4 SRD全子集回歸圖

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

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

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

3.3 ARIMA建模

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

表2 各站點(diǎn)隨機(jī)項(xiàng)時(shí)間序列ARIMA模型

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

圖7 沈陽(yáng)站隨機(jī)項(xiàng)ACF圖和PACF圖

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

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

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

表4 時(shí)空預(yù)測(cè)統(tǒng)計(jì)指標(biāo)

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

圖9 2006年12個(gè)月各站點(diǎn)時(shí)空預(yù)測(cè)誤差分布圖

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

4 結(jié)語(yǔ)

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