基于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)能有效提高模型的预测精度，具有借鉴意义;另外，该模型还适用于围岩变形、路基沉降及大坝变形等领域。