基于ARIMA模型的邊坡變形分析與預(yù)測

胡 波，譚 涵

(重慶市勘測院，重慶 400000)

由于受復(fù)雜的地質(zhì)條件的影響，邊坡穩(wěn)定問題一直是巖土工程界關(guān)注的焦點問題。為了保證工程施工和運行的安全，對邊坡監(jiān)測數(shù)據(jù)進行長期監(jiān)測并預(yù)測其變形趨勢至關(guān)重要。時間序列分析是一種認(rèn)識產(chǎn)生觀測序列的隨機機制及對序列未來的可能取值給出預(yù)測或預(yù)報的重要手段，廣泛應(yīng)用于商業(yè)、氣象學(xué)、農(nóng)業(yè)、形變監(jiān)測、醫(yī)學(xué)等領(lǐng)域[1-6]。如時序分析方法用于GNSS數(shù)據(jù)處理[7-8]，同時也可與灰色模型、小波分析、頻譜分析進行結(jié)合建模[9-11]，在大壩、基坑、橋梁等土木工程形變分析中取得了不錯的結(jié)果[12-15]。現(xiàn)有研究表明自回歸滑動平均求和(ARIMA)模型是一種常用時間序列分析模型。本文詳細(xì)論述建立ARIMA模型的關(guān)鍵步驟，并建立模型對某邊坡工程461 d的觀測數(shù)據(jù)進行分析和預(yù)測，驗證利用ARIMA模型對邊坡監(jiān)測數(shù)據(jù)進行分析與預(yù)測的可行性和有效性。

1 數(shù)據(jù)與方法

本文采用某邊坡工程從2016年6月29日—2017年10月4日期間的監(jiān)測數(shù)據(jù)成果，該數(shù)據(jù)包含了461 d中1202期觀測結(jié)果。由于時間序列分析的基礎(chǔ)數(shù)據(jù)一般為離散、等間隔的數(shù)據(jù)序列[2]，而原始觀測數(shù)據(jù)非等間隔觀測，試驗過程中采用內(nèi)插方法對數(shù)據(jù)進行重采樣，從而得到X方向和Y方向每天的平面位移的時間序列。利用R語言建立模型分析時間序列，利用前430個數(shù)據(jù)建立ARIMA預(yù)測模型，利用后31個數(shù)據(jù)對預(yù)測結(jié)果進行評估。

自回歸滑動平均(ARMA)模型是目前最常用的一種平穩(wěn)時間序列預(yù)測模型。一般來說，如果滿足

Yt=φ1Yt-1+φ2Yt-2+…+φpYt-p+et-

θ1et-1-θ2et-2-…-θqet-q

(1)

則稱{Yt}為自回歸滑動混合平均過程，階數(shù)分別為p和q，表示為ARMA(p,q)。

2 時序分析建模

2.1 數(shù)據(jù)平穩(wěn)性分析

根據(jù)觀測記錄對隨機過程的結(jié)構(gòu)進行統(tǒng)計推斷時，通常對其做出某些簡化的(大致合理的)假設(shè)，其中最重要的假設(shè)即是平穩(wěn)性。平穩(wěn)性的基本思想是：決定過程特性的統(tǒng)計規(guī)律不隨時間的變化而變化。一個隨機過程{Yt}稱為弱平穩(wěn)的條件是：均值函數(shù)在所有時間上恒為常數(shù)，且自協(xié)方差值只與滯后階數(shù)有關(guān)。通常非平穩(wěn)時間序列可以通過差分得到平穩(wěn)的時間序列。常用的時間序列平穩(wěn)性檢驗方法有圖示法和單位根檢驗法。

圖1為X方向位移和Y方向位移的時間序列圖，可以看出具有明顯的趨勢，均值隨時間不斷變化，因此可以判斷X方向和Y方向位移的時間序列為非平穩(wěn)時間序列。單位根檢驗法一般常用的有DF檢驗(Dickey-Fuller test)、ADF檢驗(augmented Dickey-Fuller test)和PP檢驗(Phillips-Perron test)。對于X方向的位移，ADF檢驗的P值為0.107 6，接受原序列非平穩(wěn)的假設(shè)；PP檢驗的P值小于0.01，拒絕原序列非平穩(wěn)的假設(shè)。對于Y方向的位移，ADF檢驗和PP檢驗的P值均小于0.01。結(jié)合時間序列圖像綜合判斷，X方向和Y方向的位移為非平穩(wěn)時間序列。

對X方向和Y方向的位移進行一階差分后的時間序列圖像如圖2所示。從圖2中可以發(fā)現(xiàn)X方向和Y方向位移一階差分后的結(jié)果基本上在0附近波動，并且ADF檢驗和PP檢驗的結(jié)果均表明，兩個一階差分后的時間序列是平穩(wěn)的。因此，可以考慮使用一階差分的ARIMA模型對這兩個時間序列建模。

2.2 時序分析建模

要確定ARIMA(p,d,q)模型的參數(shù)p、d、q的值，一般先通過樣本自相關(guān)函數(shù)(ACF)圖和樣本偏自相關(guān)函數(shù)(PACF)圖來判斷。對于AR(p)模型，其ACF圖一般呈現(xiàn)拖尾特征，而PACF圖則在滯后p階后截尾；對于MA(q)模型，其ACF圖一般在滯后q階后截尾，而PACF圖呈現(xiàn)拖尾特征；對于ARMA(p,q)模型，ACF圖和PACF圖都呈現(xiàn)拖尾特征。從圖3可以看出，一階差分后的X方向位移時間序列的1階滯后自相關(guān)函數(shù)值非常顯著且在1階以后截尾，其偏自相關(guān)函數(shù)值只在1、2、3階較為顯著。一階差分后的Y方向位移時間序列的1階滯后自相關(guān)函數(shù)值在滯后1、2、3、13階顯著，其偏自相關(guān)函數(shù)同樣在1、2、3、13階顯著。判斷這兩個時間序列可能是ARMA混合模型，需要進一步判斷其階數(shù)。

樣本ACF和PACF可以有效地識別純AR(p)或MA(q)模型，但是對于混合ARMA模型來說，它的理論ACF和PACF有著無限多的非零值，使得根據(jù)樣本ACF和PACF來識別混合模型非常困難。因此，擴展自相關(guān)法(EACF)常常被用于確定ARMA模型的階數(shù)，并且EACF法對于適度大的樣本容量具有較好的樣本性質(zhì)。圖4為兩個一階差分時間序列的EACF圖，可以看出X方向的EACF的零三角區(qū)域非常清楚地顯示出p=1、q=1的ARMA模型比較合適。而Y方向的EACF圖的零三角區(qū)域為p=2、q=2，還需進一步分析，某些系數(shù)可能是因為偶然因素在統(tǒng)計上顯著不為零。

為了得到一些可供深入研究的有用的初步模型，使用BIC檢查幾個最佳子集ARMA模型，匯總在圖5中。圖中的每一行對應(yīng)著一個子集ARMA模型，模型所選變量的單元格用陰影表示，根據(jù)BIC值將模型分類，較好的模型(有較低的BIC值)處于較高的行中，并且顏色也較深。對于X方向的時間序列，最上面一行說明具有最小BIC值的子集ARMA(7,7)模型只包括觀測時間序列誤差過程的1階滯后；其次最好的模型包括時間序列1階滯后及誤差過程的2階滯后；第3個較好的模型包括時間序列1、2階滯后和誤差過程的2階滯后。在不同的子集模型中，時間序列的1階滯后和誤差過程的2階滯后是出現(xiàn)最頻繁的變量。而對于Y方向的時間序列，最好的子集包括觀測時間序列2階滯后和誤差過程的1階滯后；其次最好的模型只包括誤差過程的1、2、6階滯后；第3個較好的模型只包括誤差過程的1階滯后。在不同的子集模型中，時間序列的2階滯后和誤差過程的1階滯后是出現(xiàn)最頻繁的變量。綜上所述，筆者選擇ARIMA(1,1,1)模型為X方向位移建立時間序列模型，選擇ARI(2,1)模型為Y方向位移建立時間序列模型。

在識別X方向和Y方向位移的時間序列模型為ARIMA(1,1,1)和ARI(2,1)后，接下來就需要估計模型的參數(shù)。常用的參數(shù)估計方法有矩估計、最小二乘估計、極大似然估計和無條件最小二乘估計。本文使用R語言中的arima函數(shù)進行參數(shù)估計，其中的method參數(shù)可以選擇參數(shù)估計的方法為“CSS”(最小二乘估計)、“ML”(極大似然估計)或“CSS-ML”，默認(rèn)的方法是“CSS-ML”，即先通過最小二乘估計計算初值，再用極大似然估計方法計算。計算結(jié)果見表1。

表1 平面位移的時間序列模型參數(shù)的極大似然估計值

注：P為似然對數(shù)值。

3 結(jié)果分析

3.1 模型殘差

為了判斷模型的擬合優(yōu)度，并給出適當(dāng)?shù)恼{(diào)整建議，往往需要分析擬合模型的殘差。殘差等于實際值減去預(yù)測值，如果模型正確識別，殘差應(yīng)當(dāng)具有以下兩條性質(zhì)：①參數(shù)估計充分接近真值，則殘差應(yīng)近似具有白噪聲的特性；②呈現(xiàn)獨立、同分布、零均值和相同標(biāo)準(zhǔn)差的正態(tài)變量。與這些性質(zhì)的偏離有助于發(fā)現(xiàn)更合適的模型。

3.2 過度擬合

過度擬合是指在識別并擬合出我們認(rèn)為合適的模型之后，擬合一個更一般的模型，即一個“接近”的模型。該模型以原始模型為特例包容原始模型，如果額外參數(shù)的估計不顯著地不為零并且兩個模型共同參數(shù)的估計與原始估計相比沒有顯著的改變，則認(rèn)為原始模型合理。針對擬合出的ARIMA(1,1,1)和ARIMA(2,1,0)模型，本文選取包容它們的ARIMA(2,1,1)模型作為一個“接近”的模型，并進行參數(shù)估計，結(jié)果見表2。

方向系數(shù)估計值標(biāo)準(zhǔn)差X方向位移ar1 0.63140.1262ar2-0.04040.0602ma1-0.77200.1177^σ2=0.4435,P=-434.36,AIC=874.72Y方向位移ar1 0.45680.0600ar20.10490.0560ma1-0.92650.0340^σ2=0.5138,P=-466.32,AIC=938.65

注：P為似然對數(shù)值。

從表2可以看出，對于X方向位移的時間序列，新增的參數(shù)ar2估值為-0.040 4，并不顯著，并且過度擬合的參數(shù)標(biāo)準(zhǔn)差都有所增大，AIC值也明顯增大。因此，認(rèn)為ARIMA(2,1,1)模型對X方向位移是過度參數(shù)化的，ARIMA(1,1,1)模型更加合理。對于Y方向位移的時間序列，新增的參數(shù)ma1估值為-0.926 5，非常顯著，AIC值也明顯變小，因此判斷ARIMA(2,1,1)模型是更合理的。其殘差檢核結(jié)果如圖7所示，可以看出結(jié)果是優(yōu)于原ARIMA(2,1,0)模型的。

3.3 模型預(yù)測

利用前430個X和Y方向的位移數(shù)據(jù)分別建立了合理的ARIMA模型，然后向后預(yù)測20步，與實測數(shù)據(jù)的結(jié)果對比，對比結(jié)果如圖8所示。圖8中黑色實心點為預(yù)測值，十字標(biāo)記為實測數(shù)據(jù)值，上下虛線為95%預(yù)測極限?？梢钥闯鰞蓚€時間序列的預(yù)測值變化都很平滑并逐漸趨近于一個均值，預(yù)測值的變化趨勢和實測數(shù)據(jù)大致一致，且實測數(shù)據(jù)值都在預(yù)測值的95%預(yù)測極限之內(nèi)，預(yù)測結(jié)果較好。

4 結(jié) 語

本文使用2016年6月29日—2017年10月4日共461 d的邊坡監(jiān)測數(shù)據(jù)進行時間序列分析。首先利用圖示法、ADF檢驗和PP檢驗分析了數(shù)據(jù)的平穩(wěn)性，發(fā)現(xiàn)X方向和Y方向位移的時間序列并不平穩(wěn)，然后通過一階差分得到了平穩(wěn)的時間序列，再利用ACF、PACF、EACF和BIC的結(jié)果來識別ARIMA模型的階數(shù)，對X方向位移得到了ARIMA(1,1,1)模型，對Y方向位移得到了ARIMA(2,1,0)模型，并使用CSS-ML方法進行了參數(shù)估計。隨后利用殘差分析和過度擬合的方法進行了模型評價，發(fā)現(xiàn)X方向位移的ARIMA模型符合標(biāo)準(zhǔn)，通過檢驗，而Y方向位移的ARIMA模型則應(yīng)該擴展為ARIMA(2,1,1)。最后利用前430個數(shù)據(jù)的模型預(yù)測了之后的31個數(shù)據(jù)，與實測數(shù)據(jù)符合較為理想。因此，通過ARIMA模型對邊坡監(jiān)測數(shù)據(jù)進行時間序列建模是一種有效的分析和預(yù)測手段，對工程施工和防災(zāi)減災(zāi)具有重大意義。