馮杭華
(1.浙江華東測繪與工程安全技術有限公司,浙江 杭州 310014)
目前監測數據分析的方法主要包括回歸分析法[1](雙曲線法、指數曲線法等)、灰色系統[2]、時間序列分析法[3]、遺傳算法[4]、人工神經網絡[5]、卡爾曼濾波[6]等方法。其中時間序列分析是通過建立合適的數學模型,擬合時間序列,并完成預測功能。ARIMA模型就是這一類的典型代表,被廣泛應用到建筑物沉降預測[7]、基坑沉降監測[8]和橋梁變形監測[9]等領域。胡波[10]等利用ARIMA模型對邊坡觀測數據組成的時間序列進行分析研究,驗證ARIMA模型預測邊坡變形的可行性。GM(1,1)模型作為灰色系統的核心,同樣被廣泛應用于各個領域[11-14],但該模型具有一定的局限性。張儀萍[15-16]等在GM(1,1)模型基礎上用多項式逼近模型中的參數,采用最小二乘法確定多項式中的待定系數,構建時變參數GM(1,1)模型。楊柳[17]等將時變參數GM(1,1)模型運用到邊坡位移預測分析并且取得很好的效果。
本文的主要內容利用權系數將ARIMA模型和時變參數GM(1,1)模型相結合形成組合模型,簡稱ARIMA-TGM組合模型。結合邊坡位移實測數據進行分析預測,將預測結果和ARIMA模型與時變參數GM(1,1)模型二者的預測結果對比,驗證ARIMA-TGM組合模型在邊坡位移變形監測中的可行性和優越性。
求和自回歸移動平均模型(ARIMA模型)是一種被廣泛應用的時間序列預測模型,該模型具有完備的數學理論基礎,且算法效率高,易于實現。ARIMA模型來自于自回歸滑動混合模型(ARMA模型),ARMA模型包括自回歸模型(AR模型)和滑動平均模型(MA模型)兩大部分。

運用ARMA模型的前提是原始序列必須為平穩的,當原始序列存在趨勢性和周期性時ARMA模型不再適用時,針對非平穩的時間序列ARIMA模型應運而生。ARIMA模型主要思想是先將一個非平穩序列時間序列{Yt},經過d次差分Wt=?dYt得到一個平穩的時間序列{Wt},然后運用ARMA模型對平穩的時間序列{Wt}進行分析預測。則稱{Yt}為求和自回歸滑動平均模型,即ARIMA(p,d,q)模型,其計算流程如圖1所示。

圖1 ARIMA模型計算流程圖
設x(0)=(x(0)(1),x(0)(2),…,x(0)(n))為已知的非負初始數據序列,一次累加生成序列為X(1)=(x(1)(1),x(1)(2),…,x(1)(n)),其中
GM(1,1)模型微分方程為:

式中,z(1)(t)=0.5×x(1)(t)+0.5×x(1)(t-1)。
傳統GM(1,1)模型的白化微分方程中參數a,b假定為固定常數,這樣不符合現實情況,王金柱[18]假定參數a,b也是關于時間的函數時,建立了時變參數GM(1,1)模型,具體建模步驟如下:
1)模型的建立,時變參數模型微分方程為:

白化方程為:

式中,a(t)和b(t)為所要求解的時變參數,可用以下多項式逼近。

則有:

為使s取最小值,對式(6)進行全微分。

將式(5)代入式(4)并展開為矩陣形式可得:

式中,Y=(x(0)(2),…,x(0)(n))T,B=(a0,ap,b0,bq)T。

根據最小二乘法,可得式(8)的解為:

2)時變參數模型的求解,將式(4)中的x(1)(t)用緊鄰均值生成代替,可解得微分方程的通解。

并且令初始條件x(1)(1)=x(1)(0)得到c=x(0)(1)。
根據復合梯形積分得式(10)的離散化形式。

式中,當t<n時,結果為模型的擬合值;當t>n時,
所求的結果為模型的預測值[17],對作一階累減得到擬合預測數據序列。
3)多項式次數的確定。時變參數模型中多項式次數p、q的最佳取值根據原始實測值與模型擬合值的平均擬合相對誤差達到最小值來確定[16]。

本文通過權系數將ARIMA模型和時變參數GM(1,1)模型相結合。采用取方差倒數的方法求解權系數的大小,它的主要思想是為了使組合預測模型的誤差平方和盡可能地小,對誤差平方和較小的模型賦予較高的權系數,而對誤差平方和較大的模型賦予較小的權系數[19]。

式中,ej為第j個模型的擬合誤差平方和。

式中,Xt為第t期實際觀測值;Xˉtj為第t期第j種預測方法的擬合值;先得到ARIMA模型和時變參數GM(1,1)模型的擬合值,根據式(12)和式(13)計算2種模型所對應的權系數,然后根據下式計算ARIMA-TGM組合模型的擬合預測值。

孔口邊坡的位移監測數據,監測數據自2010-07-05~2013-02-05,以10 d為一周期,一共32期監測數據,分別使用ARIMA模型、時變參數GM(1,1)模型和ARIMA-TGM組合模型對前20期的位移監測數據進行擬合建模,并預測后12期的邊坡位移變化量。前20周期的實測值見表1。

表1 邊坡位移監測的前20周期實測值
ARIMA模型預測過程主要分為以下4步:
1)原始序列平穩性檢驗。本文先通過觀察位移變化曲線,然后再采用ADF單位根檢驗法進行平穩性檢驗。由圖2所得原始時間序列位移值存在明顯的上升趨勢,再經過ADF檢驗發現檢驗統計量遠大于1%、5%、10%顯著水平下的臨界值,并且p為0.938遠大于0.05,所以可以確定原始序列為非平穩序列,需要進行平穩化處理。首先取d=1進行一階差分,圖3為一階差分后的時間序列圖,從圖3可以看出樣本數據在0附近波動,并且一階差分后序列的ADF檢驗統計量為-12.386,遠小于1%、5%、10%顯著水平下的臨界值,p的值為5.836e-08無限接近于0,因而經過一階差分后為平穩時間序列,故d取值為1。

圖2 實測位移的時間序列圖

圖3 一階差分后的時間序列圖
2)模型識別與定階。通過觀察差分后平穩序列的自相關函數(ACF)圖和偏自相關函數(PACF)圖截尾性等特征確定模型形式,并且根據自相關函數和偏自相關函數初步確定p和q的取值范圍,然后再結合BIC準則進一步確定ARIMA模型的最佳階數p、q,通過反復實驗,最終確定當p=0,d=1,q=1時BIC最小值為-3.174,所以ARIMA(0,1,1)為最優模型。
3)模型參數估計與模型檢驗。使用極大似然估計進行參數估計,得到模型參數θ1為0.389。模型確定后,模型擬合統計R方為0.966,說明擬合效果較好。還需要對殘差進行白噪聲檢驗。檢驗結果Ljung-BoxQ為11.803,P值為0.812,殘差檢驗為白噪聲,自相關系數均在95%可信區間,進一步說明ARIMA(0,1,1)模型可以進行預測。
4)模型預測。利用以上確定的最優ARIMA(0,1,1)模型對后12期邊坡位移進行預測,具體預測結果如表2所示。
同樣使用時變參數GM(1,1)模型對邊坡位移前20期的監測數據進行擬合建模,然后預測后12期的邊坡位移值。時變參數GM(1,1)模型的最佳多項式次數p、q的取值不能取太大值,否則會導致矩陣出現病態現象,所以本文設定p、q的取值在0~4之間,利用matlab平臺結合平均擬合相對誤差達到最小原則設計搜索算法確定最佳p、q的取值,通過搜索最終確定當p值取0、q值取1時平均擬合相對誤差達到最小,時變參數GM(1,1)模型可以取得很好的擬合效果。利用建立好的時變參數GM(1,1)模型對邊坡后12期位移值進行預測,具體預測結果見表2。
根據邊坡位移的實測值和ARIMA模型與時變參數GM(1,1)模型前20期邊坡位移的擬合結果,利用式(12)、(13)求2種模型所對應的權系數,求出賦予ARIMA模型的權系數為0.738,賦予時變參數GM(1,1)模型的權系數為0.262。分別將ARIMA模型與時變參數GM(1,1)模型后12期邊坡位移預測值乘以各自對應得權系數然后相加就得到ARIMA-TGM組合模型的邊坡位移預測值,具體預測結果見表2。

表2 3種模型的預測結果對比
由表2可以看出,ARIMA模型的預測最小絕對誤差為0.001,最大為0.828;時變參數GM(1,1)模型的預測最小絕對誤差為0.004,最大為1.188;ARIMA-TGM組合模型的預測最小絕對誤差為0.043,最大為0.741,3種模型的預測值和實際觀測值都很接近,都能很好地描述邊坡位移的變化趨勢。
為了更加客觀地對比3種預測模型的預測精度,通過計算平均絕對誤差(MAE)、平均絕對百分比誤差(MAPE)、和均方根誤差(RMSE)3種評價指標進行模型評價,具體計算結果見表3。

表3 3種模型3種指標的計算結果
由表3可得,本文的ARIMA-TGM組合模型的MAE、MAPE和RMSE 3種精度評價指標均小于ARIMA模型和時變參數GM(1,1)模型所對應的MAE、MAPE和RMSE值,所以ARIMA-TGM組合模型的預測精度更高。
本文使用權系數將ARIMA模型和時變參數GM(1,1)模型結合形成ARIMA-TGM組合模型,分別運用3種模型對邊坡位移進行擬合預測,通過計算MAE、MAPE和RMSE精度評價指標對3種模型的預測精度進行對比,結果表明ARIMA-TGM組合模型的預測精度高于單一模型的預測精度,能更加有效的對邊坡監測數據進行分析和預測,對邊坡穩定性評估和工程的防災減災具有一定的參考價值。