宋 威，劉開云，梁軍平，徐 沖

(1.北京交通大學(xué) 土木建筑工程學(xué)院, 北京 100044；2.清華大學(xué) 環(huán)境學(xué)院， 北京 100084；3.中鐵第一勘察設(shè)計(jì)院集團(tuán)有限公司， 陜西 西安 710043)

隧道圍巖位移的空間分布及隨時(shí)間的演化能夠表征圍巖的變形規(guī)律，在隧道工程的建設(shè)中，隧道位移的監(jiān)控量測具有十分重要的意義，可以利用隧道位移值的絕對大小及其變化速度對隧道圍巖穩(wěn)定程度進(jìn)行判斷，而通過數(shù)值模擬的方法，能夠?qū)λ淼赖淖冃芜^程進(jìn)行模擬復(fù)現(xiàn)，供設(shè)計(jì)施工參考。而數(shù)值模型中參數(shù)的取值決定了數(shù)值模擬的精度，由于巖土體自身的物理力學(xué)性質(zhì)復(fù)雜，受力變形過程具有高度非線性，所以準(zhǔn)確識別巖土體的物理力學(xué)參數(shù)顯得十分重要。相較于通過經(jīng)驗(yàn)或工程類比的方法來確定參數(shù)，位移反分析方法提供了一種科學(xué)的、可控的參數(shù)取值方法，在巖土工程中得到了廣泛應(yīng)用[1-3]。為解決正法反分析的計(jì)算效率問題，引入智能機(jī)器學(xué)習(xí)算法來對巖體復(fù)雜非線性系統(tǒng)進(jìn)行建模，以代替數(shù)值計(jì)算，對反分析問題進(jìn)行建模求解[4-7]。人工神經(jīng)網(wǎng)絡(luò)和單輸出支持向量機(jī)作為位移預(yù)測方法，被廣泛地采用，取得了良好的效果。

位移反分析中的數(shù)值直接解法的精度取決于正向位移預(yù)測模型的精度和優(yōu)化算法的問題尋優(yōu)能力，關(guān)于機(jī)器學(xué)習(xí)模型，之前大多數(shù)研究采用的支持向量回歸模型均為一維輸出，只能滿足單變量的預(yù)測，對于隧道變形而言，一般在同一斷面監(jiān)測得到的位移變量包括水平收斂和拱頂下沉，如果采用一維輸出SVR算法，就需要建立兩個(gè)單輸出的模型，在反分析過程中將水平收斂和拱頂沉降分開考慮，但這兩個(gè)變量實(shí)際上都是同一非線性系統(tǒng)的觀測值，兩者之間可能存在相互關(guān)系，而一維輸出SVR模型無法考慮這種相互之間的非獨(dú)立關(guān)系[8-9]，本文引入一種針對多輸出系統(tǒng)(Multiple Output System)改進(jìn)的支持向量回歸算法[10](MSVR)，能夠在同一個(gè)模型中考慮多個(gè)輸出，建立統(tǒng)一的優(yōu)化目標(biāo)，進(jìn)行模型訓(xùn)練。引入該方法考慮水平收斂和拱頂沉降的相互關(guān)系，建立統(tǒng)一的反演目標(biāo)。關(guān)于優(yōu)化算法，近年來，許多優(yōu)秀的集群智能優(yōu)化算法被用來對反分析中的優(yōu)化問題進(jìn)行求解，而人工免疫算法是模仿人體免疫機(jī)制提出的一種新的智能方法，并在近幾十年來被廣泛的應(yīng)用于非線性優(yōu)化、組合優(yōu)化、控制工程、機(jī)器人、故障診斷、圖像處理等諸多領(lǐng)域，并取得了一些成就[11-13]，本文采用免疫克隆選擇算法(Immune Clone Select Algorithm, ICSA)，對MSVR模型控制參數(shù)進(jìn)行優(yōu)化，使其對樣本的誤差最小，將其應(yīng)用于隧道位移反分析中，結(jié)合銅黃高速公路三口隧道工程進(jìn)行位移反分析，驗(yàn)證該方法在巖土工程反分析領(lǐng)域的應(yīng)用效果。

1 多維輸出支持向量回歸

1.1 一維輸出支持向量回歸

支持向量機(jī)學(xué)習(xí)方法是由Vapnik等[14]根據(jù)統(tǒng)計(jì)學(xué)的理論提出的。

(1)

式中：w為擬合系數(shù)；ξi,ξi*為松弛因子；C為懲罰系數(shù)；k為樣本數(shù)量。

約束條件為

(2)

式中：xi為第i個(gè)樣本輸入值；yi為第i個(gè)樣本輸出值；f為預(yù)測函數(shù)；ε為損失計(jì)算閾值。

對這一凸二次優(yōu)化問題，采用Lagrange乘子法求解，Lagrange函數(shù)表達(dá)式為

L(w,b,ξi,ξi*,αi,αi*,γi,γi*)=

(3)

根據(jù)KKT條件，其對偶形式最大化函數(shù)表達(dá)式為

(4)

式中：xi為第i個(gè)樣本輸入值；xj為第j個(gè)樣本輸入值。

約束條件為

(5)

求解優(yōu)化問題得到支持向量回歸模型，表達(dá)式為

(6)

式中：x為需要預(yù)測的樣本輸入值。

上述分析適用于線性回歸問題，對于非線性回歸，需要引入特征空間，使用一個(gè)非線性映射把數(shù)據(jù)映射到一個(gè)高維特征空間進(jìn)行線性回歸，在高維特征空間用核函數(shù)來代替線性回歸中的內(nèi)積運(yùn)算，核函數(shù)的定義為

K(xi,xj)=φ(xi)·φ(xj)

(7)

式中：φ為非線性映射函數(shù)。

經(jīng)過與線性回歸相同的推導(dǎo)，最后得到的支持向量模型擬合函數(shù)為

(8)

核函數(shù)種類較多，徑向基函數(shù)核函數(shù)(Radical Basic Function,RBF)是一種常用的效果較好的核函數(shù)，表達(dá)式為

(9)

式中：x,y為兩樣本輸入值；σ為核函數(shù)參數(shù)。

1.2 多維輸出支持向量回歸

傳統(tǒng)的支持向量回歸的輸出變量為一維變量，該特點(diǎn)使其應(yīng)用場景受到限制，在一些復(fù)雜的系統(tǒng)里，需要建立多輸入-多輸出的映射系統(tǒng)，而一維SVR并不能完成此類任務(wù)，因此，文獻(xiàn)[10]在一維SVR的基礎(chǔ)上進(jìn)行擴(kuò)展，使其適用于多維輸出系統(tǒng)，以解決實(shí)際工程中更為復(fù)雜的問題。

將一維ε不敏感損失函數(shù)擴(kuò)展到多維，定義損失函數(shù)，表達(dá)式為

(10)

基于式(10)所示損失函數(shù)，構(gòu)造優(yōu)化目標(biāo)函數(shù)，表達(dá)式為

(11)

為解決MSVR模型的數(shù)學(xué)優(yōu)化問題，文獻(xiàn)[10]中提出了采用迭代加權(quán)最小二乘法(Iterative Reweighted Least Squares)來求解。

在式(11)的優(yōu)化目標(biāo)函數(shù)中，將損失函數(shù)用一階泰勒展開來近似替代，即

(12)

構(gòu)造式(12)的二次近似代替，文獻(xiàn)[10]中采用的近似公式為

(13)

式中：CT為不依賴于W和b的常數(shù)項(xiàng)。

采用近似公式(13)的原因是該公式中W和b解耦合，優(yōu)化求解不需迭代，直接取對W和b的偏導(dǎo)數(shù)等于0即可計(jì)算出W和b的近似解。

算法流程如下：

Step3計(jì)算下一迭代步的解，迭代計(jì)算公式為

(14)

式中：ηk為迭代步長。

采用回溯算法來確定迭代步長。每個(gè)迭代步的初始步長設(shè)置為1(初始步長的選取對于算法的收斂非常重要)，然后檢查LP(Wk+1,bk+1)

Step4計(jì)算uik+1和αi，設(shè)置k=k+1，回到Step2，直到滿足收斂條件，達(dá)到收斂。

優(yōu)化目標(biāo)求解完畢得到使樣本集合總體損失最小的W和b，多輸出支持向量回歸模型即建立完畢，多輸出支持向量回歸稱為MSVR。

2 免疫克隆選擇算法

生物免疫系統(tǒng)是一個(gè)復(fù)雜的自適應(yīng)系統(tǒng)。人體免疫系統(tǒng)能夠識別病原體并做出應(yīng)答反應(yīng)，從而具有一定的學(xué)習(xí)、記憶和模式識別的能力?？梢詫⑵湫畔⑻幚淼脑?、機(jī)制采用計(jì)算機(jī)算法進(jìn)行描述，來解決科學(xué)和工程問題。Castro[15]最早提出了克隆選擇算法 (ICSA)，該算法是受人體免疫系統(tǒng)啟發(fā)，模擬生物免疫系統(tǒng)功能和作用機(jī)理對復(fù)雜問題進(jìn)行求解的智能方法，它保留了生物免疫系統(tǒng)所具有的若干特點(diǎn)，包括全局搜索能力；多樣性保持機(jī)制；魯棒性強(qiáng)；并行的求解搜索過程。將其引入對優(yōu)化問題進(jìn)行求解。

在ICSA算法中加入了種群抑制過程，來對種群的平均濃度進(jìn)行控制，避免算法過早的收斂到局部最優(yōu)解，增加全局優(yōu)化能力。ICSA算法的詳細(xì)過程見圖1。

圖1 ICSA算法流程圖

采用典型多峰函數(shù)Rosenbrock函數(shù)(香蕉函數(shù))檢驗(yàn)改進(jìn)后的ICSA算法的優(yōu)化能力，其函數(shù)表達(dá)式為

(15)

式中：X=[x1,x2,…,xn]∈RN

Rosenbrock函數(shù)全局最小值點(diǎn)在所有自變量取值為1時(shí)取得，函數(shù)值為0。用10元Rosenbrock函數(shù)檢驗(yàn)ICSA算法對于多元函數(shù)的優(yōu)化效果。自變量的搜索區(qū)間為(-10,10)，算法參數(shù)設(shè)置見表1。

表1 優(yōu)化算法控制參數(shù)取值

運(yùn)行10次優(yōu)化算法，大約4次找到函數(shù)的最小值點(diǎn)，這10次的優(yōu)化結(jié)果見表2。由表2可知，ICSA算法具有良好的多維多峰值函數(shù)的尋優(yōu)能力，可以應(yīng)用于求解MSVR模型的優(yōu)化問題以及基于ICSA-MSVR耦合算法的巖體物理力學(xué)參數(shù)位移反分析。

表2 Rosenbrock函數(shù)優(yōu)化結(jié)果

3 基于ICSA-MSVR耦合算法的巖體物理力學(xué)參數(shù)位移反分析

MSVR模型的建立過程中，控制參數(shù)(懲罰系數(shù)C，敏感系數(shù)ε，核函數(shù)參數(shù)σ)的取值需要人為指定，為了控制參數(shù)取值來達(dá)到最小的樣本訓(xùn)練誤差及最好的MSVR模型泛化精度，采用ICSA算法對參數(shù)進(jìn)行優(yōu)化求解。

在模型訓(xùn)練階段，定義訓(xùn)練樣本集合的總體誤差函數(shù)為優(yōu)化目標(biāo)，單個(gè)樣本的誤差同樣采用ε不敏感損失函數(shù)，見式(10)，將訓(xùn)練樣本分為學(xué)習(xí)樣本和測試樣本，采用K折交叉驗(yàn)證的方法來計(jì)算樣本整體誤差，優(yōu)化目標(biāo)函數(shù)表達(dá)式為

(16)

式中：Lall(C,ε,σ)表示整體的訓(xùn)練損失函數(shù)；k表示樣本等分的份數(shù)；km表示訓(xùn)練樣本k等份后每份的數(shù)量；上標(biāo)星號表示求得的最優(yōu)參數(shù)。

待MSVR模型訓(xùn)練完畢，即通過ICSA算法優(yōu)化求得最優(yōu)的MSVR模型參數(shù)后，即完成正法反分析過程中的位移預(yù)測模型的建立過程。

MSVR模型建立完成后，可以進(jìn)行巖體的物理力學(xué)參數(shù)辨識過程。

給定待反演巖體物理力學(xué)參數(shù)的取值范圍，優(yōu)化目標(biāo)函數(shù)為預(yù)測位移(拱頂沉降和水平收斂)的誤差函數(shù)，其表達(dá)式為

(17)

aff(x)=[fu(x)-u]2+[fv(x)-v]2

式中：x為待反演參數(shù)向量，即圍巖的物理力學(xué)參數(shù)，選取的7個(gè)參數(shù)為重度、變形模量、黏聚力、內(nèi)摩擦角、水平側(cè)壓力系數(shù)、泊松比、剪脹角；fu(x)為MSVR模型預(yù)測水平收斂；fv(x)為MSVR模型預(yù)測拱頂下沉；u為現(xiàn)場實(shí)測水平收斂；v為現(xiàn)場實(shí)測拱頂下沉；aff為親和度函數(shù)(即優(yōu)化目標(biāo)函數(shù)，最小化問題)。

優(yōu)化完成后輸出優(yōu)化結(jié)果，完成參數(shù)反演，得到巖體力學(xué)參數(shù)的辨識結(jié)果。

將辨識參數(shù)輸入到MSVR模型中，對隧道開挖的后續(xù)位移進(jìn)行預(yù)測。以驗(yàn)證辨識得到的參數(shù)是否符合工程實(shí)際。

ICSA-MSVR耦合算法的流程見圖2，在定義誤差函數(shù)和結(jié)果輸出部分模型訓(xùn)練過程和參數(shù)辨識過程有所不同。

圖2 ICSA-MSVR耦合算法流程

在ICSA算法的優(yōu)化過程中，為了擴(kuò)大參數(shù)的搜索范圍和搜索效率，將待優(yōu)化參數(shù)進(jìn)行了指數(shù)映射，即在種群中的參數(shù)的取值范圍是實(shí)際的取值范圍的自然對數(shù)，在實(shí)際親和度計(jì)算中，將抗體個(gè)體進(jìn)行指數(shù)映射，計(jì)算公式為

P′=exp(P)

(18)

式中：P為每個(gè)個(gè)體的取值(對于模型訓(xùn)練過程P為三個(gè)支持向量回歸控制參數(shù)C,ε,σ，對于力學(xué)參數(shù)反演過程P為待反演的7個(gè)巖體的物理力學(xué)參數(shù))。親和度計(jì)算采用指數(shù)映射后的個(gè)體取值P′。

4 工程實(shí)例

4.1 三口隧道工程概況

銅(陵)—黃(山)高速公路銅湯段三口隧道為雙向分離式四車道高速公路隧道，左、右線隧道間距約為46 m，左線隧道超前右線隧道施工，超前距離約為30 m，因而可以忽略右線隧道施工對左線隧道的影響。本文考察其左線出口段。該段為Ⅲ級圍巖段，里程樁標(biāo)號為ZK252+270—ZK252+511段，軸向掘進(jìn)總長度約為241 m。

三口隧道施工方法為爆破掘進(jìn)開挖，隧道研究段為全斷面一次開挖，以3 m為循環(huán)進(jìn)度，每天掘進(jìn)6 m。在隧道施工過程中設(shè)置監(jiān)測斷面，對隧道收斂和沉降變形進(jìn)行監(jiān)測，ZK252+510監(jiān)測面測點(diǎn)位移見表3。

4.2 訓(xùn)練樣本的獲取——數(shù)值試驗(yàn)

三口隧道橫斷面見圖3，通過Midas建立網(wǎng)絡(luò)模型，然后導(dǎo)入FLAC3D中進(jìn)行計(jì)算，整體數(shù)值計(jì)算模型見圖4。采用均質(zhì)材料的假設(shè)，沿隧道軸線方向?yàn)閥軸方向，垂直向上為z軸方向，x軸水平向右為正方向，模型從左至右寬為60 m；從坐標(biāo)原點(diǎn)往下長為50 m；計(jì)算模型沿縱向長80 m，單元共42 500個(gè)，節(jié)點(diǎn)共43 768個(gè)。

圖3 三口隧道斷面(單位：m)

圖4 三口隧道數(shù)值計(jì)算模型

對模型前后左右及下邊界施加法向變形約束，上邊界為自由邊界。由于沒有實(shí)測地應(yīng)力資料，垂直地應(yīng)力等于上覆巖層自重。兩個(gè)水平方向地應(yīng)力相等，等于豎向地應(yīng)力乘以水平側(cè)壓力系數(shù)。

隧道初期支護(hù)采用噴錨支護(hù)，噴射混凝土采用C25混凝土，厚度10 cm，錨桿布置長度為2.5 m，環(huán)向間距1.5 m，錨桿采用全長粘接水泥砂漿錨桿，鋼材為HRB335，直徑為25 mm。錨桿鋼材和C25噴射混凝土的本構(gòu)均為彈性本構(gòu)，參數(shù)見表4。錨桿數(shù)值模型見圖5。

表4 隧道初期支護(hù)材料力學(xué)參數(shù)

隧道采用全斷面開挖，數(shù)值模擬中循環(huán)進(jìn)尺為3 m，初期支護(hù)緊跟掌子面。變形監(jiān)測斷面設(shè)置在距離隧道出口處1 m的位置。變形監(jiān)測點(diǎn)見圖6。

圖6 變形監(jiān)測點(diǎn)布置

參考JTG D70—2014 《公路隧道設(shè)計(jì)規(guī)范》[16]，圍巖物理力學(xué)指標(biāo)見表5。

表5 圍巖物理力學(xué)參數(shù)取值范圍

我國地應(yīng)力實(shí)測資料表明水平側(cè)壓力系數(shù)范圍為0.55～2.00。剪脹角對計(jì)算結(jié)果的影響不大，剪脹角取值范圍是3°～16.5°。采用7因素28水平的均勻試驗(yàn)劃分，因素選擇為重度、變形模量、黏聚力、內(nèi)摩擦角、水平側(cè)壓力系數(shù)、泊松比、剪脹角7個(gè)因素，應(yīng)用均勻試驗(yàn)法設(shè)計(jì)28個(gè)數(shù)值試驗(yàn)方案，采用 FLAC3D 軟件進(jìn)行三口隧道左線施工的三維數(shù)值模擬，考慮不同的掌子面與 ZK252+510 監(jiān)測面的距離，從數(shù)值實(shí)驗(yàn)的結(jié)果中隨機(jī)抽取35次數(shù)值計(jì)算結(jié)果，組成供MSVR 算法進(jìn)行網(wǎng)絡(luò)訓(xùn)練的 30個(gè)訓(xùn)練樣本見表6，5個(gè)檢驗(yàn)樣本見表7。訓(xùn)練樣本采用K折交叉驗(yàn)證法進(jìn)行MSVR模型的學(xué)習(xí)及建立。表7所示的檢驗(yàn)樣本完全不參與MSVR模型的訓(xùn)練過程，只用來檢驗(yàn)訓(xùn)練得到的MSVR模型的泛化性能。MSVR采用徑向基核函數(shù)，并對樣本進(jìn)行歸一化處理，將數(shù)據(jù)線性映射到[0,1]。將拱頂下沉和水平收斂作為MSVR模型的輸出變量建立MSVR模型，對模型參數(shù)(C,ε,σ)進(jìn)行優(yōu)化。三個(gè)參數(shù)的搜索區(qū)間依次為(0,10 000)、(0,100)、(0,10 000)，經(jīng)過ICSA算法優(yōu)化后，MSVR模型最優(yōu)參數(shù)C,ε,σ的數(shù)值分別為19.4、0.004 59、2.10。

表6 MSVR模型訓(xùn)練樣本集

表7 檢驗(yàn)MSVR訓(xùn)練效果的檢驗(yàn)樣本集

為檢驗(yàn)?zāi)Ｐ头夯阅埽瑢Ρ?所示檢驗(yàn)樣本，用訓(xùn)練好的MSVR模型進(jìn)行位移預(yù)測，結(jié)果見表8。

表8 檢驗(yàn)樣本的MSVR模型預(yù)測效果

由表8可知，MSVR模型對于拱頂下沉和水平收斂的預(yù)測精度較好，表明映射圍巖物理力學(xué)參數(shù)與位移之間非線性關(guān)系的MSVR模型已經(jīng)建立，可以用于巖體物理力學(xué)參數(shù)的辨識。

4.3 巖體參數(shù)的反演及后續(xù)位移預(yù)測

利用訓(xùn)練完成的MSVR模型，對巖體物理力學(xué)參數(shù)進(jìn)行反演，7個(gè)待反演參數(shù)的反演區(qū)間見表9。將拱頂沉降和水平收斂統(tǒng)一反演，兩者權(quán)重相等，反演目標(biāo)函數(shù)見式(17)，取監(jiān)測面距離掌子面5 m處時(shí)的位移值作為反演基礎(chǔ)信息。反演得到的巖體物理力學(xué)參數(shù)見表9。

表9 待反演、反演參數(shù)取值范圍

利用反演得到的巖體力學(xué)參數(shù)，輸入MSVR模型對后續(xù)位移(d=8、11、14 m)進(jìn)行預(yù)測，預(yù)測效果見表10。由表10可知，利用ICSA-MSVR算法對隧道進(jìn)行位移反分析，其后續(xù)位移預(yù)測中，拱頂下沉和水平收斂的預(yù)測具有較好的精度，誤差最大在15%左右,說明基于免疫多輸出支持向量回歸算法的位移反分析方法具有較高的精度。

表10 反演參數(shù)預(yù)測后續(xù)位移

為進(jìn)一步說明本文所述ICSA-MSVR算法對位移反分析問題的求解效果，與文獻(xiàn)[7]中方法進(jìn)行對比，文獻(xiàn)[7]中采用遺傳-廣義回歸神經(jīng)元網(wǎng)絡(luò)算法(GRNN)對塢石隧道進(jìn)行位移反分析，基于塢石隧道的實(shí)測位移值，對隧道圍巖的物理力學(xué)參數(shù)進(jìn)行辨識。文獻(xiàn)[7]采用數(shù)值試驗(yàn)、機(jī)器學(xué)習(xí)與人工智能相結(jié)合的智能巖土工程反分析方法，與本文具有可對比性。

根據(jù)文獻(xiàn)[7]中表6所述的隧道實(shí)際監(jiān)測位移，可以確定如式(17)所示的反演目標(biāo)函數(shù)，按照本節(jié)所示工程實(shí)例的相同步驟，對塢石隧道的圍巖參數(shù)進(jìn)行辨識，本文方法和文獻(xiàn)[7]方法反分析得到的圍巖力學(xué)參數(shù)對比見表11。

表11 反演參數(shù)對比

將反演得到的巖土參數(shù)輸入到MSVR模型，進(jìn)行后續(xù)的開挖步的隧道變形的預(yù)測，對比結(jié)果見表12。由表12可知，本文ICSA-MSVR算法對隧道變形的預(yù)測精度高于單一的BP方法，與遺傳-廣義回歸神經(jīng)元網(wǎng)絡(luò)算法(GRNN)對比，精度略有提高，本文方法對拱頂沉降和水平收斂的統(tǒng)一反演效果較好。也體現(xiàn)出多輸出支持向量回歸算法與免疫克隆選擇算法相結(jié)合的良好的位移反分析應(yīng)用前景。

表12 位移預(yù)測結(jié)果對比

本文所述算法，均編制C++程序，并借助Intel OpenMP并行算法庫實(shí)現(xiàn)了算法的并行，提高其計(jì)算效率。算法計(jì)算時(shí)間見表13。

表13 算法花費(fèi)時(shí)間表 s

5 結(jié)論

(1) 本文提出了改進(jìn)的免疫克隆選擇算法，加入了種群抑制過程來改善收斂于局部極值以及算法的早熟問題，經(jīng)過算例檢驗(yàn)，該算法對多維優(yōu)化問題具有良好的求解能力，且收斂速度較快。

(2) 引入多輸出的支持向量回歸算法，用來解決輸出變量之間的依賴性以及計(jì)算效率的問題，增強(qiáng)了位移反分析方法的工程應(yīng)用價(jià)值。

(3) 將本文方法應(yīng)用于銅黃高速公路三口隧道的位移反分析中，結(jié)果顯示本文方法的參數(shù)辨識效果較好，辨識的參數(shù)較為準(zhǔn)確，對后續(xù)位移的預(yù)測精度較好，具有一定的工程應(yīng)用價(jià)值。并將本文所述ICSA-MSVR算法與遺傳-廣義回歸神經(jīng)元網(wǎng)絡(luò)算法進(jìn)行了參數(shù)辨識的對比，驗(yàn)證了本文方法的優(yōu)勢。