ARIMA模型和Holt-Winters方法在邊坡土體深層水平位移預測中的比較

晏 凱， 吳佳曄， 張應遷， 楊 露

(1.四川輕化工大學管理學院，自貢 643000；2.四川升拓檢測技術股份有限公司，自貢 643000；3.四川輕化工大學土木工程學院，自貢 643000；4.東北林業大學土木工程學院，哈爾濱 150000)

邊坡穩定性問題一直是邊坡工程的一項重要研究內容[1]。而邊坡工程是一個不確定的、非線性的復雜的系統[2]，目前國內外存在的邊坡安全監測技術也是多種多樣。如2010年Pieraccini等[3]、Barla等[4]、Gischig等[5]最早采用了GB-In SAR對邊坡的變形進行監測研究；2012年阮志新等[6]利用測斜儀對高速公路邊坡進行了監測并證明了監測的正確性；2013年王志勇等[7]利用In SAR技術對北京房山地區的滑坡進行了監測；以及2015年徐茂林等[8]利用TM30對鞍山市某露天礦邊坡位移進行了監測實驗，實現了自動化監測，等等。這些邊坡工程監測方法大多都用于邊坡變形分析，由此推斷研究邊坡變形對分析邊坡穩定性意義重大。

土體深層水平位移是邊坡變形監測中的重要指標之一。本文嘗試運用ARIMA模型[9]和霍爾特-溫特斯方法[10]，對土體深層水平位移的歷史數據進行短期預測和實證研究，并對結果加以比較，分析兩種方法的優劣。

1 問題描述

1.1 方法選取的原則

在長期的邊坡變形過程中，監測的土體深層水平位移值是隨季節周期變化，并存在一定增長(或降低)趨勢的非平穩時間序列。傳統的預測方法如回歸模型、k-最近鄰分類的IBK模型等適用于平穩時間序列，對該類值的預測結果并不理想[11]。而由Box和Jenkins提出來的自回歸求和滑動平均模型(ARIMA模型)，是將非平穩時間序列先經過d階差分平穩化，再對得到的平穩時間序列利用自回歸[autoregressive, AR(p) process]和滑動平均過程[moving average, MA(q) process]，并通過樣本自相關系數(autocorrelation coefficient, ACF)和偏自相關系數(partial autocorrelation coefficient, PCF)等數據對建立的模型進行辨識，同時還提出了一整套的建模、估計、檢驗和控制方法[12]；其次Holt-Winters預測方法，由于它可以同時處理趨勢和季節性變化[13]，并能適當地過濾隨機波動的影響。因此選用ARIMA模型和Holt-Winters預測方法來預測具有季節性波動和趨勢性增長(或減少)的非平穩時間序列問題，并對結果進行比較分析。

1.2 數據集的構建

實驗的邊坡土體深層水平位移值數據來源于四川升拓檢測股份有限公司的SQLsever數據庫，與之相關邊坡工程項目為茂縣石大關鄉拴馬村梯子槽滑坡監測項目。SQLsever數據庫是美國Microsoft公司設計的一個為分布式服務器計算的、高性能的、可擴展的關系型數據庫管理系統。

由于數據采集設備故障和外界環境等原因，導致了部分時間點的采集數據出現缺失。因此為降低預測誤差，故選取監測數據值連續而穩定的時間段來進行預測實驗。

隨機選取監測項目中某個土體深層水平位移監測點的500～650個數據。之所以選500～650個有兩個原因：一是如果數據取的太少就無法充分挖掘原序列中的信息，而選取數據過多又會因序列前后間隔太長，使前期數據與后期預測數據關聯性變弱，而造成不必要的計算量；二是查閱諸多相關的時間序列預測文獻，若對季度數據(m=4，按1年4個季度循環)實證分析的原始序列個數在50～100，月度數據(m=12，按1年12個月循環)的個數在250～350，因此對時度數據(m=24，按1天24小時循環)實證分析的原始序列個數可取在500～650。

最終選取從2018年8月15日0：00～2018年9月10日23：00的“X軸”方向的所有整點數據構成“H”數據集及序列{Ht}，數據量為648，見圖1。從圖1中可以發現，“序列{Ht}是明顯的有確定性趨勢又有周期性的時間序列。

圖1 “H”數據集的序列圖Fig.1 The sequence diagram of the “H” dataset

2 ARIMA模型的建立與預測過程

2.1 模型的建立

設{Ht}為平穩時間序列且均值為零，則ARMA(1,1)的1階自回歸1階滑動平均模型公式可表述為

Ht=φHt-1+et-θet-1

(1)

故ARMA(p,q)公式[14]表述為

Ht-φ1Ht-1-…-φpHt-p=et-θ1et-1-…-θqet-q

(2)

式(1)可簡寫為

φ(B)Ht=θ(B)et

(3)

式(3)中：et表示白噪聲參數，t=1, 2, …,p；φ表示自回歸系數；θ表示回歸參數。在模型ARIMA(p,d,q)中，d代表差分階數。因為邊坡監測項目中，土體深層水平位移的時間序列數據普遍都是不平穩的，使序列平穩的途徑有很多，差分就是其中之一，且ARIMA模型差分后的建模過程幾乎與ARMA相同。

表1 一階差分序列的單位根檢驗

注：AIC全稱AIC信息準則，于1973年由日本統計學家赤池弘次提出，又稱赤池(Akaike)信息準則[15]; SC指的是施瓦茨(Schwartz)準則，又稱貝葉斯信息準則(BIC); HQ指的是漢南-奎因(Hannan-Quinn)準則[16]。

依據上述“三準則”，AIC、SC、HQ小的較優，則由表1的單位根檢驗結果，判斷該序列為存在既無截距項又無趨勢項的單位根過程。

由表2可以看出一階差分的ADF檢驗的T統計值的絕對值為18.577 4，遠大于顯著性水平為0.05時的T統計值1.942，則存在一個單位根的原假設成立；其次P也很小。這均說明原序列的一階差分序列是平穩的，因此d=1。

表2 一階差分序列(無截距、無趨勢)的ADF檢驗

從圖2可以清晰地看到，DH的時序圖圍繞一個常數值上下波動，且波動區間不大；又從圖3中可看出，自相關系數迅速衰減趨于0。綜合所述，一階差分后原序列是平穩的，不再進行二階差分檢驗。

圖2 “H”數據集的一階差分時序圖Fig.2 The first-order difference sequence diagram of the “H” dataset

圖3 一階差分自相關圖Fig.3 The first-order differential autocorrelation graph

從圖3中可清晰地看到，P都小于0.05且接近于0，則原序列是非白噪聲的，存在一定相關性，并有規律可循。在滯后階數(lag)為4以后，偏自相關系數迅速趨近于0，初步判斷一階差分序列為偏自相關的4階截尾，可試著擬合AR(4)模型；在lag=5時，偏自相關系數處于2倍標準差的置信帶邊緣，可嘗試擬合AR(5)模型。而在lag=1時，自相關系數顯著不為零，在lag=2、8時，它又處于2倍標準差的置信帶邊緣，則可試著去擬合MA(1)、MA(2)和MA(8)模型。因此可以初步擬合的模型有ARIMA(4,1,1)、ARIMA(5,1,1)等。最后由模型定階發現，p可能等于4，q可能等于1、2、8，表3、表4分別是各個擬合模型的估計結果和預測參數對比情況。

依據P小于顯著性水平0.05且越小結果越顯著(模型較優)的判斷原則，比較表3中c、ar(p)、ma(q)的P大小，綜合得出ARIMA(3,1,1)模型最優。

先根據“三準則”，AIC、SC、HQ小的較優，對表4中各模型的檢驗統計量進行比較，得出ARIMA(3,1,1)有兩項最小(SC、HQ)；又依據F分布原則，F-statistics大的較優， ARIMA(3,1,1)的F統計量顯然最大，P也小于顯著性水平0.05。因此綜合判斷ARIMA(3,1,1)為最優擬合模型。

表3 方程的P比較

注：c表示常數；ar表示自回歸過程，括號內表示滯后階數；ma表示滑動平均過程，括號內表示滯后階數。

表4 各種模型的精度指標對比

模型ARIMA(3,1,1)的公式推導[9，14]如下：

因為：

(4)

又：

(5)

(6)

在式(1)兩邊同時乘以以(6)得：

(7)

系數代入式(7)得出ARMA(3,1)模型的公式為

Wt=0.000 061+0.381 209wt-1+0.147 211wt-2+

(8)

又令Wt=Ht-Ht-1,得ARIMA(3,1,1)模型的公式為

(Ht-Ht-1)+0.381 209(Ht-2-Ht-1)+

0.147 211(Ht-3-Ht-2)+0.127 094(Ht-4-

(9)

2.2 模型的診斷檢驗

在確定ARMA模型參數的估計值后，為了驗證擬合模型的適應性，還需進行診斷檢驗。運用Eviews時間序列分析軟件建立殘差的自相關圖，對殘差進行純隨機性檢驗，圖4是ARIMA(3,1,1)的殘差自相關檢驗結果。

圖4 ARMA(3,1,1)模型的殘差相關圖Fig.4 Residual correlation graph of the ARMA(3,1,1) model

從圖4中可以看到，自相關函數值幾乎都在0.95的置信區域內；其次，P都大于0.05，則表明殘差為白噪聲且是純隨機性的；又從ARIMA(3,1,1)模型擬合圖(圖5)中看出，擬合效果好。綜合判斷，擬合的ARIMA(3,1,1) 模型有效。

圖5 ARMA(3,1,1)模型的擬合圖Fig.5 The fitting graph of the ARMA(3,1,1) model

2.3 模型的應用

為驗證擬合模型的適用性，現用模型ARIMA(3,1,1)進行預測，預測結果見圖6，可以看出該模型的預測效果較理想。

圖6 序列靜態預測圖Fig.6 The sequence static forecasting graph

又將預測值存放在DH序列中，作出DH和DHF預測擬合圖，見圖7。從圖7中可以看出，該模型的預測效果較好。

圖7 靜態預測擬合圖Fig.7 The static forecasting fitting graph

現由在DHF中已有的預測值，利用模型ARIMA(3,1,1)的公式可反推出第一期預測值，又將預測值整合到原序列，同樣的過程可得出第二期預測值，以此類推直到得出m(m=24)期，具體結果見表6。

3 基于阻尼模型的Holt-Winters方法預測

3.1 Holt-Winters預測方法

Holt-Winters預測方法的基本思想是把具有線性趨勢、季節變動和隨機波動的時間序列進行分析研究，并與指數平滑法結合，分別對長期趨勢、趨勢的增量和季節波動做出估計，然后建立預測模型，外推預測值[17]。Holt-Winters季節模型包含加法模型、乘法模型[18]。

Holt-Winters乘法模型的預測外推公式為

y′t+h=(at+hbt)st+h-m

(10)

a、b、s的計算公式如下：

at=α(yt/st-m)+(1-α)(at-1+bt-1)

(11)

bt=β(at-at-1)+(1-β)bt-1

(12)

st=γ[yt/(at-1-bt-1)]+(1-γ)st-m

(13)

為了預測更小時間周期的數據集，在乘法模型的基礎上，添加延伸系數φ；再通過指數平滑法變換得到阻尼Holt-Winters模型的預測外推公式，如下：

y′t+h=[at+(φ+φ2+…+φh)bt]st+h-m

(14)

式中：a表示截距；b表示趨勢；s表示阻尼模型的季節因子；t表示預測時間，且t=1,2,…,m；h表示平滑期數，h>0；m表示周期長度，例如，對于季度數據取4，月度數據取12。因實驗數據間隔為小時，循環一次要24個采集數據，故取m=24。

a、b、s的計算公式如下：

at=α(yt/st-m)+(1-α)(at-1+φbt-1)

(15)

bt=β(at-at-1)+(1-β)φbt-1

(16)

st=γ[yt/(at-1-φbt-1)]+(1-γ)st-m

(17)

式中：α、β、γ為阻尼因子，也稱平滑系數，在0～1之間取值。

3.2 可行性分析

為了檢驗數據進行預測研究的可行性，利用OriginPro，對“H”數據集進行離散分析，得到如圖8的相鄰時刻數據散點圖。

圖8 “H”數據集的散點圖Fig.8 The distribution point diagram of the “H” dataset

3.3 預測與結果分析

利用WEKA數據挖掘軟件，讓α、β、γ在[0,1]的區間內進行精度為小數點后兩位的等值迭代運算，且預測誤差=|實際值-預測值|/實際值×100%，得出：當α=0.03、β=0.03、γ=0.03時，“H”數據集的平均預測誤差最小，為1.69%。

再分析平滑系數α、β、γ對預測誤差的影響大小(α、β、γ的系統默認值都為0.2)。由表5可以發現，對于“H”數據集，當β為默認值時，平均預測誤差變化最大，約為0.75%，因此推斷改變β最有可能提高預測精度。

表5 平滑系數的影響分析

根據以上分析結果，對“H”數據集，擬定α、γ為不變量，β為變量，分析預測結果得出：當β=0.04時，最小預測誤差約為1.68%，見圖9。

圖9 “H”數據集的預測誤差變化曲線Fig.9 The change curve in forecasting error of the “H” dataset

綜上所述，當取α=0.02、β=0.03、γ=0.03時，“H”數據集的平均預測誤差最小，為1.53%，具體預測結果見表6。

4 比較兩種方法的預測結果

因為ARIMA模型和Holt-Winters方法對長期預測存在局限性，這里只給出24期真實值與預測值的比較情況，如表6所示。

表6 兩種方法的預測結果及誤差

注:預測時間為“2018-09-11：00～2018-09-11：23”；方法1=ARIMA模型、方法2=Holt-Winters方法；預測誤差=(|實際值-預測值|/實際值)×100%；實際值數據來源于四川升拓檢測股份有限公司的SQLsever數據庫。

(18)

式(18)中：h1i和h′1i分別表示方法1的實際值和預測值；h2i和h′2i分別別表示方法2的實際值和預測值；i表示預測期數，i=1, 2, …, 24；n=24。

從表6可以看出，兩種方法預測誤差都偏小，且平均預測誤差只有1.73%和1.53%；又從圖11中不難看出兩種方法的預測結果曲線與實際值曲線的變化趨勢基本吻合；這都足以說明兩種方法對“H”數據集的短期預測精度都很高。

現將兩種方法的預測誤差比較發現：雖然Holt-Winters方法的最大預測誤差(5.05%)大于ARIMA模型(4.69%)，但在總體上，ARIMA模型有11期預測誤差大于2.22%，Holt-Winters方法只有5期，因此Holt-Winters方法短期預測效果更好，如圖10所示。

圖10 比較兩種方法的預測相對誤差Fig.10 Comparing relative errors of the two methods

然而，根據實際值的95%置信度，計算預測區間，可以得到兩種方法的預測曲線效果圖，見圖11。

圖11 兩種方法的預測結果曲線Fig.11 The curvilinear in forecasting results of the two methods

在圖11中，Holt-Winters方法盡管平均預測誤差更小，但預測區間所包含真實值的概率小于ARIMA模型(圖11中綠色矩形框標記的Holt-Winters方法預測點在預測區間外)，因此ARIMA模型短期預測更穩定。

綜上所述，對于土體深層水平位移的短期預測，兩種方法各有特點。

5 結論

邊坡土體深層水平位移預測是一個充滿挑戰性的問題，但時間序列預測一直被認為是對具有趨勢性增長(降低)和時變性波動序列進行統計預測的有效手段。利用ARIMA模型和Holt-Winters方法對土體深層水平位移進行了預測，可以為邊坡變形監測與預警提供一定的依據。

ARIMA模型和Holt-Winters模型都屬于短期預測模型，隨著時間的推移，它們的預測效果也會逐漸變差[19]。從兩種方法的預測結果來看，短期預測效果都很好。其次通過分析其預測誤差和結果曲線發現，兩種方法各具特點：Holt-Winters方法的短期預測效果更好；ARIMA模型的短期預測結果更穩定。又從實際操作來講，Holt-Winters方法簡便，而ARIMA模型的建模過程較為煩瑣，并且在p和q的選擇存在很多評判標準。

總之，通過兩種時間序列短期預測方法的比較研究，希望其結果對邊坡治理人員加強邊坡安全治理和預防邊坡事故發生有幫助。另外也可考慮融合其他預測方法，加強對各種影響因素的關注，如邊坡自身結構、天氣變化、工程施工等，常常這些不穩定的因素對邊坡土體深層水平位移的長期走勢也有著重要的應用價值。