， ，

(西南交通大學 交通隧道工程教育部重點實驗室， 四川 成都 610031)

1 概述

在隧道修建過程中，隧道的監(jiān)控量測工作對隧道的安全施工具有指導意義。我國大部分公路隧道采用新奧法修建，隧道拱頂位移是監(jiān)控量測的必測項目[1]。隧道的開挖會打破巖體的應力應變平衡狀態(tài)，導致應力重分布和圍巖的變形，拱頂位移便是其直觀體現(xiàn)[2]。通過對拱頂位移的監(jiān)控量測，判斷圍巖以及支護結構的狀態(tài)并預測未來的發(fā)展趨勢，將對支護設計提供重要參考，從而確保隧道的穩(wěn)定性，避免施工中的安全事故[3]。

目前，隧道拱頂下沉預測一般有3種方法，即解析法、數(shù)值分析法和人工智能方法。在解析法中， Loganathan[4]等采用鏡象法研究了粘土地層中隧道開挖導致的地表沉降。Bobet[5]等提出了飽和土中淺埋隧道的應力函數(shù)解答。Sulem[6]等提出隧道洞室的收斂函數(shù)可以通過距掌子面的位移和地層時變特性來表示。Y.-W. Pan[7]等將隧道收斂視為巖體應力狀態(tài)和流變特性的函數(shù)。但由于解析法中大量的近似和簡化，難以準確模擬巖體復雜的非線性特性和時變特性，也無法考慮施工過程的影響。在數(shù)值分析方法中，張頂立[8]等用FLAC-3D軟件模擬地鐵隧道的開挖過程，分析了不同地層條件下拱頂下沉與地表沉降的關系。汪小敏[9]等采用有限元程序對軟弱圍巖隧道施工過程進行了三維彈塑性分析，并比較了不同施工方法對控制圍巖變形的作用。楊珺博[10]等用Ansys軟件模擬了CD法開挖的穿江隧道的力學特性，總結了該工程條件下隧道拱頂下沉的特征。目前，不少計算軟件已經能夠考慮地層的各向異形、圍巖的非線性以及圍巖與隧道結構間的相互作用[11]，然而如何準確并有效地對實際工程進行數(shù)值模擬至今并沒有得到很好的解決。圍巖的力學特性十分復雜，得到準確全面的圍巖信息比較困難，而且精確地建立模型并進行驗證將耗費大量時間[12]。隨著人工智能技術地發(fā)展，通過人工智能方法解決隧道工程中的問題吸引了廣泛關注。蘇道振[13]等針對大斷面軟弱地層中的隧道施工，通過BP神經網絡預測了隧道拱頂下沉及周邊收斂。Adoko[14]等采用多元自適應回歸樣條法和人工神經網絡預測了隧道開挖過程中洞室的收斂位移。Mahdevari[15]等采用支持向量機算法，以圍巖參數(shù)為輸入值，隧道最終收斂為輸出值，預測開挖工程中洞室的最終收斂位移。在解析法和數(shù)值方法中，實際工程的隨機性和不確定性會對預測結果造成巨大影響，而人工智能方法基于已有的實測數(shù)據(jù)進行預測，因而被視為更有前景的方法。

神經網絡作為一種典型的人工智能方法，可很好地反映輸入、輸出變量的非線性關系。然而對于隧道拱頂下沉的預測，神經網絡算法要求收集大量與拱頂下沉相關的參數(shù)數(shù)據(jù)，包括隧道埋深、斷面大小、圍巖的各種特性以及含水情況[14]。此外，神經網絡方法預測拱頂下沉存在隱含層數(shù)和隱含節(jié)點數(shù)難以確定，容易陷入局部極小值，收斂速度慢等問題[16]。而支持向量機在解決小樣本，局部極小值等問題上有明顯優(yōu)點[17]?；诖?，本文采用基于支持向量機的非線性時間序列方法對拱頂下沉進行預測。目前，對于隧道相關問題的時間序列研究，大部分是針對單一時間序列。然而，在實際工程中，一般每隔5 m或10 m就會對隧道斷面的拱頂下沉進行實時監(jiān)測，尤其是對于長大公路隧道，可以得到相當數(shù)量的斷面沉降值的時間序列，這些數(shù)據(jù)理應得到充分利用。因而，本文先通過K-medoids算法對時間序列按時間序列相似性進行聚類[18]，然后對聚類中心通過LSSVM方法建立預測模型。實現(xiàn)了在斷面監(jiān)測初期對拱頂下沉的有效預測。

2 算法

2.1 動態(tài)時間彎曲距離

在許多研究中，時間序列之間的相似性度量都通過距離來表示，其中，最常用的是歐式距離。然而，對于曲線形狀相似，但時間軸不一致的時間序列，歐式距離并不適用。動態(tài)時間彎曲方法通過對時間軸的變化可以有效解決這一問題[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. 連續(xù)性：集合中元素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聚類算法用歐式距離表示數(shù)據(jù)之間的差異，通過直接求均值來計算聚類的中心點，該方法僅適用于數(shù)據(jù)對象都在同一個歐式空間中[21]。因而本文選用以數(shù)據(jù)對象作為聚類中心的K-medoids算法，數(shù)據(jù)間的差異性采用上述動態(tài)時間彎曲距離來度量。K-medoids算法的流程如下：

a.隨機選擇初始聚類中心點m1,m2,…,mk,…,mK；

b.將其他對象分配到與其距離最小的中心點所在的聚類：

(5)

c. 按照動態(tài)時間彎曲距離的平均值最小的原則更新每個聚類的中心點：

∑xi∈CKDTW(xl,xi)

(6)；

d. 重復第b.c.步直至算法收斂。

(7)

2.3 小波包變換

(8)

(9)

(10)

式中:i是頻率參數(shù)；j是尺度參數(shù)；k是平移參數(shù)；h(k)，g(k)是一組共軛鏡像濾波器函數(shù)。

小波包分解的過程就是將離散信號通過一個低通濾波器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) 在目標函數(shù)中使用二范數(shù)，用等式約束代替不等式約束，將求解優(yōu)化問題變?yōu)榍蠼饩€性方程組，提高了收斂速度和求解效率[26]。

給定訓練集(xi,yi),i=1,2,…,N，xi∈RN，xi∈R。用映射函數(shù)φ(·)將輸入數(shù)據(jù)轉化到高維空間。則目標函數(shù)為：

(17)

約束條件為：

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

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

(18)

式中:ω權向量；γ為正則化參數(shù)；e為誤差變量；b為偏置量。

建立拉格朗日函數(shù)，求解后得到最小二乘支持向量機回歸模型：

(19)

式中:K為滿足Mercer條件的核函數(shù)；xi為訓練集；x為新輸入的數(shù)據(jù)。

本文采用徑向基函數(shù)作為核函數(shù)，其具有較好的泛化能力[26]，表達式為：

(20)

正則化參數(shù)γ和核函數(shù)參數(shù)σ的優(yōu)化采用網格搜索算法。

在對時間序列通過LSSVM算法建立預測模型時，輸入變量和輸出變量的長度分別為N和1，即輸出值由其前面的N個值決定。

算法流程圖如圖2所示。

3 工程實例分析

為了驗證算法的有效性，采用汶馬高速鷓鴣山隧道左線拱頂下沉數(shù)據(jù)進行分析。四川省汶川至馬爾康高速公路項目起于汶川縣城以南鳳坪壩，止于卓克基，路線全長173 km。鷓鴣山隧道為本項目全線控制性工程。右幅隧道起訖樁號為YK179+730～YK188+496 m，長8766 m；左幅隧道起訖樁號為ZK179+696～ZK188+486，長8790 m，最大埋深約1350 m。隧道區(qū)出露的地層主要為第四系全新統(tǒng)，三疊系上統(tǒng)新都橋組，三疊系上統(tǒng)侏倭組及三疊系中統(tǒng)雜谷腦組，巖性主要為變質砂巖、板巖、千枚巖組成。該隧道為雙向四車道高速公路，凈寬10.25 m，凈高5 m。隧道內輪廓拱高715 cm，采用上半圓半徑為553 cm的三心圓斷面。采用新奧法施工，二臺階開挖。

鷓鴣山隧道施工過程中的拱頂下沉采用全站儀進行監(jiān)測。圍巖主要為IV、V級圍巖，其中，IV級圍巖測量間距為10 m左右，V級圍巖測量間距為5 m左右，在圍巖變化出應增設監(jiān)測點。在保證爆破作業(yè)不會破壞測點的基礎上，測點應靠近工作面埋設，通常距工作面0.5～2 m。在下一循環(huán)開挖前，必須獲得拱頂下沉的初始讀數(shù)，一般在12 h內測得[27]。

3.1 K-medoids聚類結果分析

選取鷓鴣山隧道左線ZK185+030至ZK185+780共100個斷面(斷面間間距為10 m或5 m)進行基于動態(tài)時間彎曲距離(DTW)的K-medoids聚類分析。在實際監(jiān)控量測過程中，拱頂下沉的量測頻率與圍巖級別以及下沉速度等因素有關，因此拱頂下沉時間序列是非等間距的，應首先對時間序列進行預處理。本案例采用3次樣條插值的方法間非等間距的時間序列轉化為等間距的時間序列。預處理后，計算這100個時間序列兩兩之間的DTW距離，衡量時間序列之間的相似程度。以斷面ZK185+740和ZK185+750為算例(見圖3)，距離矩陣：

圖3 算例斷面的拱頂位移曲線Figure 3 Vault displacement curves of examples

動態(tài)時間彎曲距離：DTW=10°

構造一個100×100的相似度矩陣，以聚類數(shù)為2～6分別聚類200次，統(tǒng)計不同聚類數(shù)出現(xiàn)頻次最高的聚類中心，并對聚類結果用DB指標進行評價(見表1)。為了確定各聚類數(shù)對應的聚類中心，對同一聚類數(shù)不同聚類中心的DB值進行比較(見圖4)。由表1可見，DB值在聚類數(shù)為4時取最小值。由圖4可見，對于同一聚類數(shù)，出現(xiàn)頻次越高的聚類中心對應的DB值越小。因此將拱頂下沉時間序列分為4類，聚類中心取聚類數(shù)為4時出現(xiàn)頻次最高的聚類中心。聚類中心的拱頂下沉位移時間序列如圖5所示。

表1 不同聚類數(shù)頻次最高的聚類中心Table 1 The highest-frequency clustering centers of different clustering number聚類數(shù)最高頻次聚類中心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小波基函數(shù)將原始時間序列進行2層分解，產生4個分離的時間序列。以斷面ZK185+055為例，分解過后的時間序列如圖6所示。對于分解后的這4個時間序列，用LSSVM算法建立預測模型，這4個時間序列的時間長度為38，輸入變量的長度為6，輸出變量的長度為1，則訓練樣本個數(shù)為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

為了判斷模型的準確性，對實際監(jiān)測值和模型預測值進行對比，以聚類中心斷面ZK185+055自身和該類中斷面ZK185+070為例，數(shù)據(jù)點的對比和擬合線如圖8所示。均方根誤差(RMSE)分別為0.056 mm和0.293 mm，平均絕對百分比誤差(MAPE)分別為0.18%和1.08%。如果不對原始數(shù)據(jù)進行小波包分解，直接用LSSVM算法建立預測模型，監(jiān)測值和預測值的對比如圖9所示。均方根誤差(RMSE)分別為0.202 mm和1.266 mm，平均絕對百分比誤差(MAPE)分別為0.59%和3.89%。分別采用WPT-LSSVM和LSSVM建模所得的預測結果如圖10所示。

按照相同的方法得到其他各類的預測模型，統(tǒng)計每一類中所有斷面的RMSE值和MAPE值，將結果與直接用最小二乘支持向量機建立模型的結果進行比較(見表2)。

圖8 監(jiān)測值與預測值(WPT-LSSVM)對比Figure 8 Comparison between monitored values and predicted

圖9 監(jiān)測值與預測值(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聚類中心聚類個數(shù)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 拱頂位移外推預測

針對當前正在監(jiān)測的斷面ZK185+020，對其前10 d的拱頂下沉位移進行監(jiān)測，得到一個時間長度為10的時間序列，將該時間序列與上述4個聚類中心的前10個點構成的時間序列計算DTW距離，斷面ZK185+020與聚類中心斷面ZK185+055的DTW值最小。以斷面ZK185+055所建立的預測模型預測ZK185+020的發(fā)展趨勢，并與之后的實測值進行對比(見圖11)，RMSE為0.515 mm， MAPE為1.29%。由圖可見，對于當前正在監(jiān)測的斷面，基于K-medoids聚類和LSSVM結合的算法具有較高的準確性，且該算法僅監(jiān)測當前斷面的前幾個時間點就可以進行預測，為可能出現(xiàn)較大下沉值的斷面的處置措施提供更為及時的參考依據(jù)。而基于該斷面拱頂下沉時間序列的自身特征進行預測則需要足夠多的時間點的監(jiān)測數(shù)據(jù)才可以建模。但該算法的預測結果在下臺階開挖時間點誤差較大，由于鷓鴣山隧道采用二臺階開挖法，造成的最終累積誤差較小。對于多步開挖的隧道，該算法需要進一步優(yōu)化。

圖11 斷面ZK185+020拱頂位移外推預測Figure 11 Extrapolating prediction of vault displacements at section ZK185+020

4 結論

對隧道拱頂下沉的實時監(jiān)測以及準確預測是隧道施工安全的重要保障。本文采用基于時間序列聚類和最小二乘支持向量機結合的智能算法預測隧道拱頂下沉，得到了如下結論：

a. 基于拱頂下沉時間序列相似度進行聚類，以聚類中心為對象，通過WPT-LSSVM算法建立的預測模型在預測同一聚類中其他對象的發(fā)展趨勢上具有較高的準確性，初步驗證了同一類時間序列用同一模型預測的可行性。

b. 隧道拱頂下沉是不平穩(wěn)時間序列，趨勢項和隨機項具有不同的變化特性，相比直接用LSSVM建立模型，采用WPT方法和LSSVM結合的算法具有更高的預測精度，能夠更真實地反映拱頂下沉時間曲線。

c. 對于當前正在監(jiān)測的斷面，如果僅利用該斷面拱頂下沉時間序列的自身特征進行預測，隨著訓練點減小預測點增多，誤差將不斷累積。而基于聚類中心的預測可以充分利用已監(jiān)測斷面的數(shù)據(jù)，在當前斷面的監(jiān)測初期就可以進行有效預測。