胡 波,谭 涵
(重庆市勘测院,重庆 400000)
由于受复杂的地质条件的影响,边坡稳定问题一直是岩土工程界关注的焦点问题。为了保证工程施工和运行的安全,对边坡监测数据进行长期监测并预测其变形趋势至关重要。时间序列分析是一种认识产生观测序列的随机机制及对序列未来的可能取值给出预测或预报的重要手段,广泛应用于商业、气象学、农业、形变监测、医学等领域[1-6]。如时序分析方法用于GNSS数据处理[7-8],同时也可与灰色模型、小波分析、频谱分析进行结合建模[9-11],在大坝、基坑、桥梁等土木工程形变分析中取得了不错的结果[12-15]。现有研究表明自回归滑动平均求和(ARIMA)模型是一种常用时间序列分析模型。本文详细论述建立ARIMA模型的关键步骤,并建立模型对某边坡工程461 d的观测数据进行分析和预测,验证利用ARIMA模型对边坡监测数据进行分析与预测的可行性和有效性。
本文采用某边坡工程从2016年6月29日—2017年10月4日期间的监测数据成果,该数据包含了461 d中1202期观测结果。由于时间序列分析的基础数据一般为离散、等间隔的数据序列[2],而原始观测数据非等间隔观测,试验过程中采用内插方法对数据进行重采样,从而得到X方向和Y方向每天的平面位移的时间序列。利用R语言建立模型分析时间序列,利用前430个数据建立ARIMA预测模型,利用后31个数据对预测结果进行评估。
自回归滑动平均(ARMA)模型是目前最常用的一种平稳时间序列预测模型。一般来说,如果满足
Yt=φ1Yt-1+φ2Yt-2+…+φpYt-p+et-
θ1et-1-θ2et-2-…-θqet-q
(1)
则称{Yt}为自回归滑动混合平均过程,阶数分别为p和q,表示为ARMA(p,q)。
根据观测记录对随机过程的结构进行统计推断时,通常对其做出某些简化的(大致合理的)假设,其中最重要的假设即是平稳性。平稳性的基本思想是:决定过程特性的统计规律不随时间的变化而变化。一个随机过程{Yt}称为弱平稳的条件是:均值函数在所有时间上恒为常数,且自协方差值只与滞后阶数有关。通常非平稳时间序列可以通过差分得到平稳的时间序列。常用的时间序列平稳性检验方法有图示法和单位根检验法。
图1为X方向位移和Y方向位移的时间序列图,可以看出具有明显的趋势,均值随时间不断变化,因此可以判断X方向和Y方向位移的时间序列为非平稳时间序列。单位根检验法一般常用的有DF检验(Dickey-Fuller test)、ADF检验(augmented Dickey-Fuller test)和PP检验(Phillips-Perron test)。对于X方向的位移,ADF检验的P值为0.107 6,接受原序列非平稳的假设;PP检验的P值小于0.01,拒绝原序列非平稳的假设。对于Y方向的位移,ADF检验和PP检验的P值均小于0.01。结合时间序列图像综合判断,X方向和Y方向的位移为非平稳时间序列。
对X方向和Y方向的位移进行一阶差分后的时间序列图像如图2所示。从图2中可以发现X方向和Y方向位移一阶差分后的结果基本上在0附近波动,并且ADF检验和PP检验的结果均表明,两个一阶差分后的时间序列是平稳的。因此,可以考虑使用一阶差分的ARIMA模型对这两个时间序列建模。
要确定ARIMA(p,d,q)模型的参数p、d、q的值,一般先通过样本自相关函数(ACF)图和样本偏自相关函数(PACF)图来判断。对于AR(p)模型,其ACF图一般呈现拖尾特征,而PACF图则在滞后p阶后截尾;对于MA(q)模型,其ACF图一般在滞后q阶后截尾,而PACF图呈现拖尾特征;对于ARMA(p,q)模型,ACF图和PACF图都呈现拖尾特征。从图3可以看出,一阶差分后的X方向位移时间序列的1阶滞后自相关函数值非常显著且在1阶以后截尾,其偏自相关函数值只在1、2、3阶较为显著。一阶差分后的Y方向位移时间序列的1阶滞后自相关函数值在滞后1、2、3、13阶显著,其偏自相关函数同样在1、2、3、13阶显著。判断这两个时间序列可能是ARMA混合模型,需要进一步判断其阶数。
样本ACF和PACF可以有效地识别纯AR(p)或MA(q)模型,但是对于混合ARMA模型来说,它的理论ACF和PACF有着无限多的非零值,使得根据样本ACF和PACF来识别混合模型非常困难。因此,扩展自相关法(EACF)常常被用于确定ARMA模型的阶数,并且EACF法对于适度大的样本容量具有较好的样本性质。图4为两个一阶差分时间序列的EACF图,可以看出X方向的EACF的零三角区域非常清楚地显示出p=1、q=1的ARMA模型比较合适。而Y方向的EACF图的零三角区域为p=2、q=2,还需进一步分析,某些系数可能是因为偶然因素在统计上显著不为零。
为了得到一些可供深入研究的有用的初步模型,使用BIC检查几个最佳子集ARMA模型,汇总在图5中。图中的每一行对应着一个子集ARMA模型,模型所选变量的单元格用阴影表示,根据BIC值将模型分类,较好的模型(有较低的BIC值)处于较高的行中,并且颜色也较深。对于X方向的时间序列,最上面一行说明具有最小BIC值的子集ARMA(7,7)模型只包括观测时间序列误差过程的1阶滞后;其次最好的模型包括时间序列1阶滞后及误差过程的2阶滞后;第3个较好的模型包括时间序列1、2阶滞后和误差过程的2阶滞后。在不同的子集模型中,时间序列的1阶滞后和误差过程的2阶滞后是出现最频繁的变量。而对于Y方向的时间序列,最好的子集包括观测时间序列2阶滞后和误差过程的1阶滞后;其次最好的模型只包括误差过程的1、2、6阶滞后;第3个较好的模型只包括误差过程的1阶滞后。在不同的子集模型中,时间序列的2阶滞后和误差过程的1阶滞后是出现最频繁的变量。综上所述,笔者选择ARIMA(1,1,1)模型为X方向位移建立时间序列模型,选择ARI(2,1)模型为Y方向位移建立时间序列模型。
在识别X方向和Y方向位移的时间序列模型为ARIMA(1,1,1)和ARI(2,1)后,接下来就需要估计模型的参数。常用的参数估计方法有矩估计、最小二乘估计、极大似然估计和无条件最小二乘估计。本文使用R语言中的arima函数进行参数估计,其中的method参数可以选择参数估计的方法为“CSS”(最小二乘估计)、“ML”(极大似然估计)或“CSS-ML”,默认的方法是“CSS-ML”,即先通过最小二乘估计计算初值,再用极大似然估计方法计算。计算结果见表1。
表1 平面位移的时间序列模型参数的极大似然估计值
注:P为似然对数值。
为了判断模型的拟合优度,并给出适当的调整建议,往往需要分析拟合模型的残差。残差等于实际值减去预测值,如果模型正确识别,残差应当具有以下两条性质:①参数估计充分接近真值,则残差应近似具有白噪声的特性;②呈现独立、同分布、零均值和相同标准差的正态变量。与这些性质的偏离有助于发现更合适的模型。
过度拟合是指在识别并拟合出我们认为合适的模型之后,拟合一个更一般的模型,即一个“接近”的模型。该模型以原始模型为特例包容原始模型,如果额外参数的估计不显著地不为零并且两个模型共同参数的估计与原始估计相比没有显著的改变,则认为原始模型合理。针对拟合出的ARIMA(1,1,1)和ARIMA(2,1,0)模型,本文选取包容它们的ARIMA(2,1,1)模型作为一个“接近”的模型,并进行参数估计,结果见表2。
方向系数估计值标准差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为似然对数值。
从表2可以看出,对于X方向位移的时间序列,新增的参数ar2估值为-0.040 4,并不显著,并且过度拟合的参数标准差都有所增大,AIC值也明显增大。因此,认为ARIMA(2,1,1)模型对X方向位移是过度参数化的,ARIMA(1,1,1)模型更加合理。对于Y方向位移的时间序列,新增的参数ma1估值为-0.926 5,非常显著,AIC值也明显变小,因此判断ARIMA(2,1,1)模型是更合理的。其残差检核结果如图7所示,可以看出结果是优于原ARIMA(2,1,0)模型的。
利用前430个X和Y方向的位移数据分别建立了合理的ARIMA模型,然后向后预测20步,与实测数据的结果对比,对比结果如图8所示。图8中黑色实心点为预测值,十字标记为实测数据值,上下虚线为95%预测极限。可以看出两个时间序列的预测值变化都很平滑并逐渐趋近于一个均值,预测值的变化趋势和实测数据大致一致,且实测数据值都在预测值的95%预测极限之内,预测结果较好。
本文使用2016年6月29日—2017年10月4日共461 d的边坡监测数据进行时间序列分析。首先利用图示法、ADF检验和PP检验分析了数据的平稳性,发现X方向和Y方向位移的时间序列并不平稳,然后通过一阶差分得到了平稳的时间序列,再利用ACF、PACF、EACF和BIC的结果来识别ARIMA模型的阶数,对X方向位移得到了ARIMA(1,1,1)模型,对Y方向位移得到了ARIMA(2,1,0)模型,并使用CSS-ML方法进行了参数估计。随后利用残差分析和过度拟合的方法进行了模型评价,发现X方向位移的ARIMA模型符合标准,通过检验,而Y方向位移的ARIMA模型则应该扩展为ARIMA(2,1,1)。最后利用前430个数据的模型预测了之后的31个数据,与实测数据符合较为理想。因此,通过ARIMA模型对边坡监测数据进行时间序列建模是一种有效的分析和预测手段,对工程施工和防灾减灾具有重大意义。