基于EMD-ARIMA模型的邊坡位移分析與預測

王江榮， 劉 碩， 靳存程， 劉靜芳

(蘭州石化職業技術學院信息處理與控制工程學院，甘肅 蘭州 730060)

對高速公路路基邊坡進行安全性監測并根據監測時間序列對邊坡變形準確預測是有效控制其變形、科學防災減災的重要舉措和保障。現有的變形預測方法主要有：灰色理論[1-2]、自回歸移動平均模型[3]、神經網絡[4]、支持向量機[5]、參數模型[6]和非參數模型[7]等，這些方法均取得了較好的預測效果，但也存在著對監測數據適應能力不強、對數據質量要求高、預測精度不高等問題；另外，這些方法在建模時未能體現監測數據所蕰含的物理地質特性，大多方法僅考慮數據本身，而當數據質量或數量不足及規律性較差時，往往預測精度不高。本文在已有研究成果的基礎上，提出一種基于經驗模態分解技術(empirical mode decomposition，EMD)[8]和自回歸移動平均(auto regressive integrated moving average model,ARIMA)[9-10]相結合的時間序列預測方法，該方法首先利用EMD技術將監測時間序列分解成頻率不等、地理物理特性不同的子序列(或稱模態分量)，這些子序列通常呈非線性及非平穩性，適合用ARIMA模型描述和預測，然后將ARIMA模型的擬合預測結果通過線性組合(線性回歸)生成路基邊坡的水平位移或沉降位移。實證分析表明EMD-ARIMA模型具有較高的預測精度，能夠滿足工程需要。

1 基于EMD的邊坡水平位移時間序列分解

1.1 EMD簡介

EMD技術是一種數據序列分解方法[11-12]，該技術能夠自適應地將復雜非線性、非平穩高速公路路基邊坡地表位移監測時間序列S(t)(t=1,2,…，m,m是總監測期數)分解成若干個不同頻率的子序列即本征模態函數IMF(intrinsic mode function,IMF)和一個殘差余量r(趨勢項序列)。記IMFi為IMF的第i(i=1,2,…,n)個分量(即S(t)的子序列)，則IMFi包含了原地表位移時間序列不同時間尺度的物理地質局部特征。將IMFi(i=1,2,…,n)和剩余分量r進行疊加便得到原時間序列的重構，即

(1)

EMD 得到的各分量IMFi(t)相對分解前的位移監測時間序列S(t)具有更簡單的波動規律，通過分析IMFi(t)(彼此獨立)和r(t)殘差分量的頻率特征，并對不同頻率的分量分別建立ARIMA模型，可提高集成模型的預測精度。

說明一點，EMD技術是根據數據列(信號)本身固有的特征而自適應分解，即從時間尺度出發，先把數據信號中特征時間尺度最小(即頻率最高)的模態函數IMF1分離出來，然后分離出特征時間尺度較大的模態函數IMF2(頻率較高)，依此下去，最后分離出特征時間尺度最大的分量即余量r(頻率最小)，而余量通常是單調函數(或不明顯存在極值)。不同的數據信號EMD分解出的模態個數不盡相同，模態個數取決于數據特征時間尺度，或受數據特征時間尺度驅動。

1.2 數據來源及工程簡介

數據來自重慶奉云高速公路K1360+500～+660段位于滑坡中部ZK2-3附近JC1-6監測點的累計水平位移。該路段小里程方向接鳳凰隧道，大里程方向接孫家溝特大橋，盡管路基邊坡已采取加固措施，但因周邊工程開挖及雨水浸蝕等影響，該段路塹路基面發生變形開裂，后緣裂縫已發育直達橋臺基礎地基范圍內，一旦裂縫出現持續擴大變形，必然會造成橋臺不穩定。圖1給出了2015年3月21日到2015年11月23日間所完成的37期監測數據。

圖1 監測點JC1-6處的累計水平位移時序圖

從圖1可看出，盡管在監測期間對剪出口采取填土反壓及抗滑樁施工治理措施，但因受周邊工程開挖、水文地理環境、土壤物理特性、大氣溫度、雨水等因素影響，監測出的位移時間序列出現了較大的波動性，也呈現出較明顯的非線性和非平穩性，挖掘數據序列中所蘊含的變形信息，對研究高速公路路基邊坡穩定性具有重要意義。

1.3 EMD分解結果分析

對原始邊坡位移監測時間序列進行EMD分解，得到了3個IMF分量和1個殘差余量r，如圖2所示。

圖2 地表水平位移EMD分解結果

由圖2可看出，EMD技術將原始監測時間序列自適應地分解成4 個頻率不等的IMF 分量IMF1～IMF3和余量r，揭示了不同因素影響下公路路基邊坡變形的波動特征，殘余序列r屬于低頻成分，代表了邊坡水平變形隨時間變化的總趨勢。其中，分量IMF1～IMF3的波動周期性特點較明顯，但波動周期非均勻且不穩定，這種非均勻和不穩定性反映了邊坡土體受土壤的物理性質、地質環境、雨水、溫度及周邊施工等因素影響其內部結構所發生的變化情況，非線性特點較明顯。可見，各模態IMF分量揭示了蘊含在邊坡變形數據序列中的物理地質特性(或特征)，且除了序列IMF1曲線外，其他分量的變化曲線較原序列曲線( 如圖1所示) 要光滑和平穩，這有助于提高ARIMA模型對邊坡變形分析與預報的準確度。

2 ARIMA模型介紹

ARIMA(p,d,q)模型是建立在差分運算基礎上的自回歸移動平均ARMA 模型(其中p,q分別是自回歸模型AR和移動平均模型MA的階數，d為差分階數)。因ARMA模型適用于平穩非白噪聲時間序列的擬合預測分析，所以當時間序列非平穩時應進行差分運算使其平穩化。對于平穩時間序列{yt}(t表示時間，yt表示序列值)可建立自回歸移動平均ARMA(p,q)模型：

yt=φ0+φ1yt-1+φ2yt-2+…+φpyt-p+

ξt-θ1ξt-1-θ2ξt-2-…-θpξt-p

(2)

式中：φ0是模型的常數項；φ1,…φp是自回歸模型AR(p)的系數；θ1，θ2，…，θp是移動平均模型MA(q)的系數；ξt,ξt-1,…,ξt-p是時間t,t-1,…,t-q的隨機誤差。

當θ1=θ2=…=θp=0時，模型(2)成為自回歸模型AR(p)：

yt=φ0+φ1yt-1+φ2yt-2+…+φpyt-p+ξt

(3)

當φ0=φ1=…=φp=0時，模型(2)為移動平均模型MA(q)：

yt=ξt-θ1ξt-1-θ2ξt-2-…-θpξt-p

(4)

模型的具體形式，可以借助模型統計特性來識別(稱作人為識別法)。

借助AIC或BIC信息準則法識別模型，即分別計算AR(p)、MA(q)及ARMA(p,q)的AIC值(p,q=0,1,2,…,5)，以最小的AIC值或BIC值所對應的模型為選定模型[13]。考慮到AR(p)模型是線性方程估計，而MA和ARMA模型是非線性估計(參數估計較因難)，故在實際建模時可用高階(取較大的p值)AR模型替代MA或ARMA模型[14]。

3 EMD-ARIMA模型預測

3.1 模型預測流程

由圖2知，IMF1是高頻且伴有突變現象的時間序列，這與邊坡位移的漸變性相悖,可將其視為噪聲信號，在建模時可不予考慮，而其他分量即IMF2、IMF3及余量r則漸趨平緩且光滑，可視為反映邊坡位移的真實信號;其中序列r可看作理想的邊坡位移量，它較好地反映了邊坡位移變化的趨勢。將去掉高頻后的IMF2+IMF3+r看成逼近實際邊坡位移的信號，為了提高逼近能力和逼近水平，采用IMF2、IMF3及r的線性組合即下面的表達式：

(5)

圖3 EMD-ARIMA模型預測流程

3.2 預測過程

基本思路：選用前32期的原始監測水平位移構建EMD-ARIMA模型，然后再對后5期的水平位移進行預測。具體如下：

步驟一，分別建立模態分量IMF2、IMF3及余量r的擬合模型AMMA(p，q)，對模型的階數p和q按相對最優模型識別法確定，即計算模型ARMA(p,q)所有p和q(p，q均小于等于5)組合的AIC信息量，取其中使信息量AIC達到最小的p和q為模型的階數。本著易操作和盡可能使模型結構簡單化，在不影響精度的前提下，經試算IMF2、IMF3及余量r的擬合模型結構設定成ARMA(4,0)時均取得了理想效果(擬合精度很高)，此時模型變成自回歸模型AR(4)。

步驟二，將IMF2、IMF3及余量r時間序列導入EViews軟件[15]，在“Equation Specification”命令窗口按模型結構輸入待估參數(按規定格式輸入)，然后選擇“LS-Least Squares(NLS and ARMA)”并執行運算，輸出模型參數(含常數項)及模型性能指標。由參數估計得到3個自回歸模型 AR(4)如下：

IMF2的擬合模型，

yt=0.419 038+2.708 109yt-1-3.416 578yt-2

+2.319 444yt-3-0.704 687yt-4

(6)

IMF3的擬合模型，

yt=-0.006 059+3.402 036yt-1-4.631 039yt-2

+2.982 214yt-3-0.774 002yt-4

(7)

r余量的擬合模型，

yt=21.033 59+3.854 744yt-1-5.621 552yt-2

+3.676 619yt-3-0.910 097yt-4

(8)

相關檢驗結果見表1。數據表明3個模型的擬合效果是極顯著的。

表1 3個模型檢驗結果

步驟三，構建EMD-ARIMA預測模型。即

(9)

1.014 526·Yimf3(t)+1.022 983·Yr(t)

(10)

模型檢驗:決定系數(擬合優度)R2=0.973 7，調整后的決定系數R′2=0.971 3，均方根誤差RMSE=1.867 6，殘差平方和SSE=129.06，MAPE=0.195 6，表明該模型擬合效果是極顯著的，可以用于后期邊坡水平位移預測。

利用步驟二中的式(6)、式(7)及式(8)對后5期即33～37期的模態分量IMF2、IMF3和余量r分別預測，將預測值代入式(9)得后5期的邊坡水平位移量，預測值見表2。EMD-ARIMA模型的擬合預測效果及殘差見圖4。

表2 預測結果及比較

圖4 EMD-ARMA模型的擬合預測效果及殘差圖

由表2知預測值的最大絕對誤差為1.778 2 mm，最小絕對誤差為0.361 3 mm,平均絕對誤差為1.156 0 mm，說明模型具有很高的預測精度，即外推能力強，能夠滿足工程需要。

作為對比，再給出下面兩種形式的預測模型：

(1)不采取EMD分解技術，直接利用前32期監測數據建立ARIMA模型，借助AIC信息準則識別法將模型結構確定為ARIMA(2,1,2),即有:

yt=13.051 16+1.885 423yt-1-0.912 305yt-2

+ξt+0.969 266ξt-1-0.288 398ξt-2

(11)

式(11)對后5期的預測結果為(29.515 4，35.150 1，39.451 8，41.166 7，42.139 7)，預測值的最大絕對誤差為4.694 6 mm，最小絕對誤差為0.358 2

mm,平均絕對誤差為1.937 2 mm。

(2)不采取EMD分解技術，直接利用前32期監測數據建立AR(4)模型(AR(p)，p≥5時模型精度沒有得到改善)，即有:

yt=13.739 48+1.091 022yt-1-0.179 053yt-2

+0.626 871t-3-0.605 78t-4

(12)

式(12)對后5期的預測結果為(30.524 1，36.026 0，37.229 5，42.666 5，42.314 7)，預測值的最大絕對誤差為3.685 9 mm，最小絕對誤差為1.076 5 mm,平均絕對誤差為2.100 5 mm。

綜上分析知EMD-ARIMA預測模型高于ARIMA(2,1,2)模型和AR(4)模型，為路基邊坡變形預測分析提供了一種新思路和新方法。另外，重慶奉云高速公路K1360+500～+660路段邊坡在2015年8月之前，累積水平位移為34.21 mm，變形不斷加劇。而在2015年8月之后，由于采取了剪出口填土反壓及抗滑樁施工治理，截止2015年11月23日，累積水平位移為43.42 mm，這段期間的水平變形增加量9.21 mm，變形趨于穩定，滑坡處于基本穩定狀態。

為進一步驗證EMD-ARIMA預測模型的有效性，以1～25期的原始監測水平位移量建模并對26～37期的水平位移量進行預測，結果見表3。

表3 EMD-ARIMA模型對后12期的預測結果及誤差分析

續表

由表3知預測值的最大絕對誤差為3.547 3 mm，最小絕對誤差為0.070 3 mm,平均絕對誤差為1.493 4 mm，說明模型具有較高的預測精度，即外推能力強，能夠滿足工程需要。

說明一點，本文雖然僅以重慶奉云高速公路K1360+500～+660段監測點JC1-6的時間序列驗證了EMD-ARIMA模型的有效性，但對其他監測點同樣適用，限于篇幅不再討論。

4 結束語

路基邊坡土體易受土壤物理地質現象、雨水、氣溫及周邊環境等因素影響會發生不同程度的變形，而對邊坡土體實施監測并以此預測未來的變形是預防和有效控制變形的重要手段。實例分析表明“分解-建模-預測-合成”的建模方法(例如EMD-ARIMA)能有效提高模型的預測精度，具有借鑒意義;另外，該模型還適用于圍巖變形、路基沉降及大壩變形等領域。