基于时间序列聚类和LSSVM的隧道拱顶位移预测

2019-03-14，，

公路工程 2019年1期

，，

(西南交通大学交通隧道工程教育部重点实验室，四川成都 610031)

1 概述

在隧道修建过程中，隧道的监控量测工作对隧道的安全施工具有指导意义。我国大部分公路隧道采用新奥法修建，隧道拱顶位移是监控量测的必测项目[1]。隧道的开挖会打破岩体的应力应变平衡状态，导致应力重分布和围岩的变形，拱顶位移便是其直观体现[2]。通过对拱顶位移的监控量测，判断围岩以及支护结构的状态并预测未来的发展趋势，将对支护设计提供重要参考，从而确保隧道的稳定性，避免施工中的安全事故[3]。

目前，隧道拱顶下沉预测一般有3种方法，即解析法、数值分析法和人工智能方法。在解析法中， Loganathan[4]等采用镜象法研究了粘土地层中隧道开挖导致的地表沉降。Bobet[5]等提出了饱和土中浅埋隧道的应力函数解答。Sulem[6]等提出隧道洞室的收敛函数可以通过距掌子面的位移和地层时变特性来表示。Y.-W. Pan[7]等将隧道收敛视为岩体应力状态和流变特性的函数。但由于解析法中大量的近似和简化，难以准确模拟岩体复杂的非线性特性和时变特性，也无法考虑施工过程的影响。在数值分析方法中，张顶立[8]等用FLAC-3D软件模拟地铁隧道的开挖过程，分析了不同地层条件下拱顶下沉与地表沉降的关系。汪小敏[9]等采用有限元程序对软弱围岩隧道施工过程进行了三维弹塑性分析，并比较了不同施工方法对控制围岩变形的作用。杨珺博[10]等用Ansys软件模拟了CD法开挖的穿江隧道的力学特性，总结了该工程条件下隧道拱顶下沉的特征。目前，不少计算软件已经能够考虑地层的各向异形、围岩的非线性以及围岩与隧道结构间的相互作用[11]，然而如何准确并有效地对实际工程进行数值模拟至今并没有得到很好的解决。围岩的力学特性十分复杂，得到准确全面的围岩信息比较困难，而且精确地建立模型并进行验证将耗费大量时间[12]。随着人工智能技术地发展，通过人工智能方法解决隧道工程中的问题吸引了广泛关注。苏道振[13]等针对大断面软弱地层中的隧道施工，通过BP神经网络预测了隧道拱顶下沉及周边收敛。Adoko[14]等采用多元自适应回归样条法和人工神经网络预测了隧道开挖过程中洞室的收敛位移。Mahdevari[15]等采用支持向量机算法，以围岩参数为输入值，隧道最终收敛为输出值，预测开挖工程中洞室的最终收敛位移。在解析法和数值方法中，实际工程的随机性和不确定性会对预测结果造成巨大影响，而人工智能方法基于已有的实测数据进行预测，因而被视为更有前景的方法。

神经网络作为一种典型的人工智能方法，可很好地反映输入、输出变量的非线性关系。然而对于隧道拱顶下沉的预测，神经网络算法要求收集大量与拱顶下沉相关的参数数据，包括隧道埋深、断面大小、围岩的各种特性以及含水情况[14]。此外，神经网络方法预测拱顶下沉存在隐含层数和隐含节点数难以确定，容易陷入局部极小值，收敛速度慢等问题[16]。而支持向量机在解决小样本，局部极小值等问题上有明显优点[17]。基于此，本文采用基于支持向量机的非线性时间序列方法对拱顶下沉进行预测。目前，对于隧道相关问题的时间序列研究，大部分是针对单一时间序列。然而，在实际工程中，一般每隔5 m或10 m就会对隧道断面的拱顶下沉进行实时监测，尤其是对于长大公路隧道，可以得到相当数量的断面沉降值的时间序列，这些数据理应得到充分利用。因而，本文先通过K-medoids算法对时间序列按时间序列相似性进行聚类[18]，然后对聚类中心通过LSSVM方法建立预测模型。实现了在断面监测初期对拱顶下沉的有效预测。

2 算法

2.1 动态时间弯曲距离

在许多研究中，时间序列之间的相似性度量都通过距离来表示，其中，最常用的是欧式距离。然而，对于曲线形状相似，但时间轴不一致的时间序列，欧式距离并不适用。动态时间弯曲方法通过对时间轴的变化可以有效解决这一问题[19]。

假设时间序列Q，C，长度分别为m，n：

Q=q1,q2,…,qi,…,qm

(1)

C=c1,c2,…,cj,…,cn

(2)

定义n×m阶矩阵D=[d(qi,cj)]，其中，d(qi,cj)为qi，ci两点间的距离，一般取d(qi,cj)=(qi-cj)2。定义集合W=w1,w2,…,wk,…,wK，wK为矩阵D中的元素，该集合即为弯曲路径。弯曲路径必须满足以下3个条件：

a. 边界条件：w1=D(1,1)，wK=D(m,n)；

b. 连续性：集合中元素wk=D(a,b)，w(k-1)=D(a′,b′)应满足a-a′≤1及b-b′≤1；

c. 单调性：集合中元素wk=D(a,b)，w(k-1)=D(a′,b′)应满足a-a′≥0及b-b′≥0。

在诸多满足上述条件的弯曲路径中，一般而言只关心最小的路径长度。路径长度采用迭代算法计算，定义路径至元素D(i,j)的最小路径长度为γ(i,j)：

γ(i,j)=D(i,j)+min{γ(i-1,j-1),

γ(i-1,j),

γ(i,j-1) }

(3)

则:

(4)

为了降低计算的时间复杂度，弯曲路径将被限制在对角线附近的一定范围内[20]。

2.2 K-medoids聚类算法

最常用的K-means聚类算法用欧式距离表示数据之间的差异，通过直接求均值来计算聚类的中心点，该方法仅适用于数据对象都在同一个欧式空间中[21]。因而本文选用以数据对象作为聚类中心的K-medoids算法，数据间的差异性采用上述动态时间弯曲距离来度量。K-medoids算法的流程如下：

a.随机选择初始聚类中心点m1,m2,…,mk,…,mK；

b.将其他对象分配到与其距离最小的中心点所在的聚类：

(5)

c. 按照动态时间弯曲距离的平均值最小的原则更新每个聚类的中心点：

∑xi∈CKDTW(xl,xi)

(6)；

d. 重复第b.c.步直至算法收敛。

(7)

2.3 小波包变换

(8)

(9)

(10)

式中:i是频率参数；j是尺度参数；k是平移参数；h(k)，g(k)是一组共轭镜像滤波器函数。

小波包分解的过程就是将离散信号通过一个低通滤波器H和一个高通滤波器G进行滤波[25]，分解过程如图1。第j层和第j+1层有如下关系：

(11)

(12)

(13)

(14)

(15)

(16)

图1 三层小波包分解过程Figure 1 The process of three-level WPT decomposition

2.4 LSSVM算法

和标准支持向量机不同，最小二乘支持向量机(LSSVM) 在目标函数中使用二范数，用等式约束代替不等式约束，将求解优化问题变为求解线性方程组，提高了收敛速度和求解效率[26]。

给定训练集(xi,yi),i=1,2,…,N，xi∈RN，xi∈R。用映射函数φ(·)将输入数据转化到高维空间。则目标函数为：

(17)

约束条件为：

图2 算法流程图Figure 2 Workflow of the algorithm

yi=ωTφ(xi)+b+ei,i=1,2,…,N

(18)

式中:ω权向量；γ为正则化参数；e为误差变量；b为偏置量。

建立拉格朗日函数，求解后得到最小二乘支持向量机回归模型：

(19)

式中:K为满足Mercer条件的核函数；xi为训练集；x为新输入的数据。

本文采用径向基函数作为核函数，其具有较好的泛化能力[26]，表达式为：

(20)

正则化参数γ和核函数参数σ的优化采用网格搜索算法。

在对时间序列通过LSSVM算法建立预测模型时，输入变量和输出变量的长度分别为N和1，即输出值由其前面的N个值决定。

算法流程图如图2所示。

3 工程实例分析

为了验证算法的有效性，采用汶马高速鹧鸪山隧道左线拱顶下沉数据进行分析。四川省汶川至马尔康高速公路项目起于汶川县城以南凤坪坝，止于卓克基，路线全长173 km。鹧鸪山隧道为本项目全线控制性工程。右幅隧道起讫桩号为YK179+730～YK188+496 m，长8766 m；左幅隧道起讫桩号为ZK179+696～ZK188+486，长8790 m，最大埋深约1350 m。隧道区出露的地层主要为第四系全新统，三叠系上统新都桥组，三叠系上统侏倭组及三叠系中统杂谷脑组，岩性主要为变质砂岩、板岩、千枚岩组成。该隧道为双向四车道高速公路，净宽10.25 m，净高5 m。隧道内轮廓拱高715 cm，采用上半圆半径为553 cm的三心圆断面。采用新奥法施工，二台阶开挖。

鹧鸪山隧道施工过程中的拱顶下沉采用全站仪进行监测。围岩主要为IV、V级围岩，其中，IV级围岩测量间距为10 m左右，V级围岩测量间距为5 m左右，在围岩变化出应增设监测点。在保证爆破作业不会破坏测点的基础上，测点应靠近工作面埋设，通常距工作面0.5～2 m。在下一循环开挖前，必须获得拱顶下沉的初始读数，一般在12 h内测得[27]。

3.1 K-medoids聚类结果分析

选取鹧鸪山隧道左线ZK185+030至ZK185+780共100个断面(断面间间距为10 m或5 m)进行基于动态时间弯曲距离(DTW)的K-medoids聚类分析。在实际监控量测过程中，拱顶下沉的量测频率与围岩级别以及下沉速度等因素有关，因此拱顶下沉时间序列是非等间距的，应首先对时间序列进行预处理。本案例采用3次样条插值的方法间非等间距的时间序列转化为等间距的时间序列。预处理后，计算这100个时间序列两两之间的DTW距离，衡量时间序列之间的相似程度。以断面ZK185+740和ZK185+750为算例(见图3)，距离矩阵：

图3 算例断面的拱顶位移曲线Figure 3 Vault displacement curves of examples

动态时间弯曲距离：DTW=10°

构造一个100×100的相似度矩阵，以聚类数为2～6分别聚类200次，统计不同聚类数出现频次最高的聚类中心，并对聚类结果用DB指标进行评价(见表1)。为了确定各聚类数对应的聚类中心，对同一聚类数不同聚类中心的DB值进行比较(见图4)。由表1可见，DB值在聚类数为4时取最小值。由图4可见，对于同一聚类数，出现频次越高的聚类中心对应的DB值越小。因此将拱顶下沉时间序列分为4类，聚类中心取聚类数为4时出现频次最高的聚类中心。聚类中心的拱顶下沉位移时间序列如图5所示。

表1 不同聚类数频次最高的聚类中心Table 1 The highest-frequency clustering centers of different clustering number聚类数最高频次聚类中心DB值2ZK185+705,ZK185+5801.243ZK185+380,ZK185+640,ZK185+7001.664ZK185+055,ZK185+115,ZK185+600,ZK185+7300.985ZK185+075,ZK185+125,ZK185+140,ZK185+350,ZK185+7201.256ZK185+045,ZK185+070,ZK185+125,ZK185+300,ZK185+375,ZK185+7601.87

图4 不同聚类中心的DB值Figure 4 DB value of different clustering centers

图5 聚类中心的拱顶位移曲线Figure 5 vault displacement curves of clustering centers

3.2 LSSVM预测结果分析

首先对聚类中心ZK185+055、ZK185+115、ZK185+600、ZK185+730这4个断面的拱顶下沉时间序列进行小波包分解，采用DB16小波基函数将原始时间序列进行2层分解，产生4个分离的时间序列。以断面ZK185+055为例，分解过后的时间序列如图6所示。对于分解后的这4个时间序列，用LSSVM算法建立预测模型，这4个时间序列的时间长度为38，输入变量的长度为6，输出变量的长度为1，则训练样本个数为38-6=32，通过LSSVM预测后的时间序列如图7所示。将这4个新的时间序列对应的预测模型组合即得到原始时间序列的预测模型。然后，将该模型应用于断面ZK185+055所在的类的其他所有序列。

图6 分解后的拱顶位移时间序列Figure 6 Decomposed time series of original vault displacements

图7 预测的拱顶位移时间序列Figure 7 Predicted time series of vault displacements values

为了判断模型的准确性，对实际监测值和模型预测值进行对比，以聚类中心断面ZK185+055自身和该类中断面ZK185+070为例，数据点的对比和拟合线如图8所示。均方根误差(RMSE)分别为0.056 mm和0.293 mm，平均绝对百分比误差(MAPE)分别为0.18%和1.08%。如果不对原始数据进行小波包分解，直接用LSSVM算法建立预测模型，监测值和预测值的对比如图9所示。均方根误差(RMSE)分别为0.202 mm和1.266 mm，平均绝对百分比误差(MAPE)分别为0.59%和3.89%。分别采用WPT-LSSVM和LSSVM建模所得的预测结果如图10所示。

按照相同的方法得到其他各类的预测模型，统计每一类中所有断面的RMSE值和MAPE值，将结果与直接用最小二乘支持向量机建立模型的结果进行比较(见表2)。

图8 监测值与预测值(WPT-LSSVM)对比Figure 8 Comparison between monitored values and predicted

图9 监测值与预测值(LSSVM)对比Figure 9 Comparison between monitored values and predicted values

图10 WPT-LSSVM与LSSVM预测结果对比Figure 10 Comparison of predicted results between WPT-LSSVM and LSSVM

由图10和表2可知，预测模型在预测同一聚类中其他对象的拱顶下沉位移值上具有较高的预测精度，在实际工程中具有参考意义。相比直接采用最小二乘支持向量机建模，先通过小波包变换进行分解再采用最小二乘支持向量机建模的方法在预测隧道拱顶下沉位移值时具有更高的准确性和可靠性。

表2 WPT-LSSVM和LSSVM模型预测结果对比Table 2 The comparison of predicting results between WPT-LSSVM and LSSVM聚类中心聚类个数RMSE/mmMAPE/%WPT-LSSVMLSSVMWPT-LSSVMLSSVMZK185+055210.2121.3020.613.91ZK185+115170.2211.3180.664.10ZK185+600380.2011.2240.583.68ZK185+730240.1971.1890.553.55

3.3 拱顶位移外推预测

针对当前正在监测的断面ZK185+020，对其前10 d的拱顶下沉位移进行监测，得到一个时间长度为10的时间序列，将该时间序列与上述4个聚类中心的前10个点构成的时间序列计算DTW距离，断面ZK185+020与聚类中心断面ZK185+055的DTW值最小。以断面ZK185+055所建立的预测模型预测ZK185+020的发展趋势，并与之后的实测值进行对比(见图11)，RMSE为0.515 mm， MAPE为1.29%。由图可见，对于当前正在监测的断面，基于K-medoids聚类和LSSVM结合的算法具有较高的准确性，且该算法仅监测当前断面的前几个时间点就可以进行预测，为可能出现较大下沉值的断面的处置措施提供更为及时的参考依据。而基于该断面拱顶下沉时间序列的自身特征进行预测则需要足够多的时间点的监测数据才可以建模。但该算法的预测结果在下台阶开挖时间点误差较大，由于鹧鸪山隧道采用二台阶开挖法，造成的最终累积误差较小。对于多步开挖的隧道，该算法需要进一步优化。