隧洞圍巖損失位移估計(jì)的智能優(yōu)化反分析

張 研，蘇國(guó)韶，燕柳斌

（1. 廣西大學(xué) 土木建筑工程學(xué)院，南寧 530004；2. 廣西大學(xué) 工程防災(zāi)與結(jié)構(gòu)安全教育部重點(diǎn)實(shí)驗(yàn)室，南寧 530004；3. 廣西大學(xué) 廣西防災(zāi)減災(zāi)與工程安全重點(diǎn)實(shí)驗(yàn)室，南寧 530004）

1 引 言

隧洞施工過程中，監(jiān)測(cè)斷面的布置一般滯后于掌子面開挖，監(jiān)測(cè)儀器到位前圍巖發(fā)生的位移可稱為損失位移。損失位移的求取對(duì)于圍巖穩(wěn)定性的合理評(píng)價(jià)具有重要意義，特別是當(dāng)損失位移較大時(shí)，若將監(jiān)測(cè)到的少部分位移作為圍巖最終的收斂位移，必然給隧洞的安全性和正常工作性能的評(píng)價(jià)帶來較大的偏差。此外，自新奧法施工理念出現(xiàn)以來，利用圍巖位移量測(cè)信息進(jìn)行信息化施工已成為現(xiàn)代隧洞工程施工的一種趨勢(shì)，然而，在信息化施工過程中，如果忽視損失位移信息，將給施工安全帶來不利影響，并會(huì)導(dǎo)致圍巖反分析結(jié)果出現(xiàn)較大誤差[1]。因此，如何獲取損失位移成為隧洞圍巖穩(wěn)定性評(píng)價(jià)和安全施工的關(guān)鍵問題。

多年來，人們進(jìn)行損失位移求解時(shí)，常根據(jù)經(jīng)驗(yàn)引入荷載（或位移）釋放系數(shù)，進(jìn)而求得相應(yīng)的等效釋放荷載，再依據(jù)Boltzmann疊加原理得到監(jiān)測(cè)位移[2]，然而荷載（或位移）釋放系數(shù)一般依靠經(jīng)驗(yàn)來確定，存在一定的盲目性。近年來，楊志法等[3]考慮松動(dòng)圈的影響，提出將損失位移求取問題轉(zhuǎn)化為反分析圍巖多層分區(qū)彈性模量問題，對(duì)提高位移反分析精度起到了很大作用；張傳慶等[4]以Hoek經(jīng)驗(yàn)公式的通式為基礎(chǔ)，基于現(xiàn)場(chǎng)監(jiān)測(cè)收斂位移曲線，提出了求取損失位移的方法，取得了一定進(jìn)展。然而由于問題的高度復(fù)雜性，現(xiàn)有的理論和方法還難以完全滿足工程實(shí)踐的要求，因此，尋求更為合理有效的圍巖損失位移求取方法是很有必要的。

監(jiān)測(cè)收斂位移是隧洞開挖過程中獲得的最直接、最有效的信息，以監(jiān)測(cè)位移作為優(yōu)化反分析的依據(jù)是損失位移求取的有效途徑。本文在數(shù)值方法的基礎(chǔ)上，將損失位移的求取問題轉(zhuǎn)化為一種優(yōu)化反分析問題，然后再利用優(yōu)化方法進(jìn)行求取。

由于巖體結(jié)構(gòu)的復(fù)雜性，優(yōu)化反分析的目標(biāo)優(yōu)化函數(shù)常是一種復(fù)雜的非線性多峰函數(shù)，采用傳統(tǒng)的梯度優(yōu)化方法往往只能獲得局部最優(yōu)解。近年來，采用遺傳算法（genetic algorithm,GA）、粒子群優(yōu)化算法（particle swarm optimization,PSO）等隨機(jī)性全局優(yōu)化算法求解優(yōu)化反分析問題成為了一種趨勢(shì)[5－6]。但采用隨機(jī)性全局優(yōu)化算法時(shí)，為了評(píng)價(jià)隨機(jī)解的優(yōu)劣，往往需要對(duì)大量的隨機(jī)解進(jìn)行成千上萬次的適應(yīng)度評(píng)價(jià)。由于巖體結(jié)構(gòu)的復(fù)雜性，很難建立決策變量與適應(yīng)度函數(shù)的顯式表達(dá)式，往往需要借助數(shù)值計(jì)算建立這種隱含關(guān)系，而數(shù)值計(jì)算是相對(duì)耗時(shí)的，多次的適應(yīng)度評(píng)價(jià)使得總的計(jì)算耗時(shí)過于浩大，從而導(dǎo)致隨機(jī)性全局優(yōu)化算法在工程應(yīng)用中受到了極大限制。為了有效地減少優(yōu)化過程中適應(yīng)度評(píng)價(jià)次數(shù)以提高損失位移優(yōu)化反分析的求解效率，本文將 PSO[7－8]、高斯過程機(jī)器學(xué)習(xí)方法（Gaussian process,GP）[9－10]、FLAC3D數(shù)值計(jì)算方法[11]相結(jié)合，提出隧洞圍巖損失位移反分析的PSO-GP-FLAC3D方法，為隧洞圍巖損失位移求取提供一條新的途徑。

2 損失位移優(yōu)化反分析的計(jì)算模型

2.1 損失位移優(yōu)化反分析的總體思路

假設(shè)有一隧洞，借助數(shù)值計(jì)算模型，容易獲得開挖后圍巖某測(cè)點(diǎn)的計(jì)算位移u隨時(shí)間t變化的全過程曲線，即位移全收斂計(jì)算曲線。保持計(jì)算模型和本構(gòu)關(guān)系不變，給定不同的圍巖力學(xué)參數(shù)能夠計(jì)算獲得不同的位移全收斂計(jì)算曲線。同時(shí)，將實(shí)測(cè)位移依照對(duì)應(yīng)的監(jiān)測(cè)時(shí)間描繪于同一坐標(biāo)系中，該測(cè)點(diǎn)的兩條位移全收斂計(jì)算曲線和一條位移實(shí)測(cè)曲線見圖1。以實(shí)際位移監(jiān)測(cè)開始時(shí)刻為分隔點(diǎn)，圖1中的兩條位移全收斂計(jì)算曲線均可分為兩部分:一部分是損失位移，另一部分是監(jiān)測(cè)位移。從圖中可看出，與計(jì)算曲線1相比較，計(jì)算曲線2的監(jiān)測(cè)位移部分與位移實(shí)測(cè)曲線較接近，由此可認(rèn)為，計(jì)算曲線2的可信度較高，進(jìn)而可以推斷計(jì)算曲線2包含的損失位移即為可能的損失位移。如何獲取計(jì)算曲線 2，就成為了問題的關(guān)鍵?？梢圆捎脙?yōu)化反分析方法來推求曲線 2，這時(shí)損失位移的求解就轉(zhuǎn)化成一種優(yōu)化反分析問題。

圖1 圍巖位移變形曲線Fig.1 Displacement curves of surrounding rocks

2.2 優(yōu)化反分析的目標(biāo)函數(shù)與決策變量

優(yōu)化反分析以計(jì)算位移曲線與實(shí)際監(jiān)測(cè)位移曲線的誤差最小化為優(yōu)化目標(biāo)，以圍巖力學(xué)參數(shù)為決策變量，通過優(yōu)化參數(shù)，搜索出計(jì)算位移曲線與實(shí)際監(jiān)測(cè)位移曲線最為接近的位移全收斂計(jì)算曲線，該曲線中包含的損失位移部分即為相對(duì)可信的損失位移。

為了評(píng)價(jià)計(jì)算收斂位移曲線與實(shí)際監(jiān)測(cè)收斂位移曲線的接近程度，依據(jù)監(jiān)測(cè)時(shí)間確定用于比較的樣本點(diǎn)。由此，優(yōu)化反分析的優(yōu)化目標(biāo)函數(shù)可設(shè)為

式中:f(x)為目標(biāo)函數(shù)；x為決策變量；m為監(jiān)測(cè)總時(shí)間；n為位移測(cè)點(diǎn)總數(shù)；dit為第i個(gè)位移監(jiān)測(cè)點(diǎn)在t時(shí)刻的位移計(jì)算值；為第i個(gè)位移監(jiān)測(cè)點(diǎn)在t時(shí)刻的位移實(shí)測(cè)值。

約束條件為

式中:常數(shù)a、b分別為決策變量取值區(qū)間的下限和上限。

建立上述優(yōu)化問題的優(yōu)化目標(biāo)函數(shù)與約束條件后，就可采用優(yōu)化方法進(jìn)行求解。

3 損失位移優(yōu)化反分析智能優(yōu)化方法

3.1 粒子群優(yōu)化方法

PSO算法最早由Eberhart和Kennedy博士首先提出[7]，它是模擬鳥類群體捕食行為的一種仿生智能優(yōu)化算法。由于其具有算法簡(jiǎn)單、尋優(yōu)能力強(qiáng)等優(yōu)點(diǎn)，近年來受到各個(gè)領(lǐng)域?qū)W者的關(guān)注。

算法首先隨機(jī)初始化一群粒子，然后調(diào)整速度、位置，通過多次迭代來尋找最優(yōu)值。在每次迭代中，粒子通過跟蹤本身迄今為止尋找到的最優(yōu)值，即個(gè)體極值Pid；當(dāng)前整個(gè)粒子群找到的最優(yōu)值，即全局極值 Pgd，以及粒子前一次迭代的狀態(tài)來調(diào)整當(dāng)前迭代步的位置和速度。粒子速度和位置的調(diào)整方程為

式中:vid（k）為第i粒子在第k次迭代中第d維度上的速度；xid（k）為第i粒子在第k次迭代中第d維度上的位置；c1和c2均為學(xué)習(xí)因子，通常c1=c2=1.8～2.0；r1、r2為（0，1）之間的均勻分布的隨機(jī)數(shù)；w為慣性權(quán)重。為了避免算法后期在全局最優(yōu)解附近難以收斂這個(gè)問題，w值應(yīng)隨算法迭代的進(jìn)行而減小，一般依據(jù)下式進(jìn)行線性減小:

式中:wmax、wmin分別為最大、最小慣性權(quán)重；t為當(dāng)前迭代數(shù)；tmax為最大迭代數(shù)，一般wmax=0.9，wmin=0.4。

3.2 高斯過程機(jī)器學(xué)習(xí)方法

GP是基于高斯隨機(jī)過程與貝葉斯學(xué)習(xí)理論發(fā)展起來的一種新的機(jī)器學(xué)習(xí)方法。由于其具有良好的適應(yīng)性，能夠很好地處理高維數(shù)、小樣本、非線性等復(fù)雜問題，近年來得到廣泛應(yīng)用[12－13]。在保證性能的前提下，與神經(jīng)網(wǎng)絡(luò)(ANN)和支持向量機(jī)(SVM)相比，高斯過程具有靈活的非參數(shù)推斷、參數(shù)自適應(yīng)獲取、實(shí)現(xiàn)過程簡(jiǎn)單等突出優(yōu)點(diǎn)。

以統(tǒng)計(jì)學(xué)理論的觀點(diǎn)分析，GP是隨機(jī)變量的1種集合。其任意有限隨機(jī)變量集合均服從聯(lián)合高斯分布，即對(duì)任意一組n維隨機(jī)變量X，與其對(duì)應(yīng)的t時(shí)刻的過程狀態(tài)()fX 的聯(lián)合概率分布是n維高斯分布。當(dāng)給定它的均值m(t)和協(xié)方差函數(shù) k(t，t′)時(shí)，()fX 在t時(shí)刻的過程狀態(tài)()f t完全確定，GP定義式表示如下:

假設(shè)有 n組觀察數(shù)據(jù)的集合 D={(xi,yi) |i=1,...,n}作為GP模型的訓(xùn)練樣本集合，其中:xi為d維輸入矢量，yi∈R為觀察目標(biāo)值。如果X為d×n維輸入矩陣，y表示輸出矢量，那么訓(xùn)練樣本集合可表示為D=(X,y)。

對(duì)于新的d維輸入矢量 X*，學(xué)習(xí)過的GP模型可以依據(jù)先驗(yàn)知識(shí)預(yù)測(cè)出與X*相對(duì)應(yīng)的輸出矢量y*。

實(shí)際建模過程中通常要考慮噪聲的影響，則帶高斯噪聲的標(biāo)準(zhǔn)線性回歸模型為

式中:f(X)為回歸函數(shù)值；ε為噪聲，且符合均值為0、方差為的高斯分布。

根據(jù)貝葉斯學(xué)習(xí)理論，觀察目標(biāo)值y的先驗(yàn)分布為

式中:(=K K X,)X 為n×n階對(duì)稱正定的協(xié)方差矩陣，矩陣中的任一項(xiàng)kij度量了xi和xj的相關(guān)性；I為單位矩陣。

n個(gè)訓(xùn)練樣本的觀察目標(biāo)矢量y和1個(gè)測(cè)試樣本的回歸函數(shù)輸出值f*所形成的聯(lián)合高斯先驗(yàn)分布為

式中:K(X,X*)為測(cè)試點(diǎn) X*與訓(xùn)練樣本集合的所有輸入點(diǎn)X的n×1階協(xié)方差矩陣；k(X*,X*)為測(cè)試點(diǎn) X*本身的協(xié)方差矩陣，可分別簡(jiǎn)寫為K(X*)、k(X*)。常用的協(xié)方差函數(shù)為

式中:超參數(shù)l、σf、σn的選取對(duì)學(xué)習(xí)與預(yù)測(cè)結(jié)果的影響較大。在GP回歸模型中，最優(yōu)超參數(shù)可采用極大似然法自適應(yīng)獲得[9]；xp、xq根據(jù)情況可以是訓(xùn)練樣本之間、測(cè)試樣本之間或者訓(xùn)練樣本與測(cè)試樣本之間的變量組合；δpq是符號(hào)函數(shù)，當(dāng)p=q時(shí)，δpq=0，反之 δpq=1。

GP學(xué)習(xí)完成后，將根據(jù)貝葉斯原理在訓(xùn)練集先驗(yàn)分布的基礎(chǔ)上預(yù)測(cè)出與X*對(duì)應(yīng)的最可能的輸出值。貝葉斯原理能夠利用觀察到的真實(shí)數(shù)據(jù)不斷更新概率預(yù)測(cè)分布，當(dāng)給定新的輸入 X*、訓(xùn)練集的輸入值X和觀察目標(biāo)值y的條件下，推斷出f*預(yù)測(cè)值的均值和方差為

3.3 PSO-GP-FLAC3D優(yōu)化反分析方法

3.3.1 方法構(gòu)建的思路

采用粒子群優(yōu)化-高斯過程協(xié)同優(yōu)化方法[14]，結(jié)合FLAC3D數(shù)值計(jì)算軟件進(jìn)行損失位移優(yōu)化反分析。其中，將式（1）作為適應(yīng)度函數(shù)。

為了進(jìn)一步提高優(yōu)化反分析的效率，本文對(duì)粒子群優(yōu)化-高斯過程協(xié)同優(yōu)化方法做了如下的改進(jìn)。

（1）對(duì)于全局最優(yōu)解的超前預(yù)測(cè)

原方法在每次循環(huán)中僅利用一次GP機(jī)器學(xué)習(xí)模型在全局空間下對(duì)可能最優(yōu)解進(jìn)行預(yù)測(cè)，本文通過構(gòu)建內(nèi)層循環(huán)，循環(huán)產(chǎn)生q（q=1,2,…,q）代粒子（q為內(nèi)循環(huán)次數(shù)），多次采用GP機(jī)器學(xué)習(xí)模型對(duì)粒子適應(yīng)度進(jìn)行評(píng)價(jià)，并尋找q代粒子中適應(yīng)度最小的粒子即全局空間下預(yù)測(cè)可能最優(yōu)解，作為預(yù)測(cè)全局最優(yōu)解。這樣就能充分發(fā)揮GP機(jī)器學(xué)習(xí)模型對(duì)粒子適應(yīng)度進(jìn)行評(píng)價(jià)，無需調(diào)用數(shù)值計(jì)算這一優(yōu)點(diǎn)，提高尋找全局最優(yōu)解的效率。

（2）對(duì)于尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)的動(dòng)態(tài)更新

原方法是將當(dāng)前迭代步的最優(yōu)粒子信息作為新的知識(shí)源添入知識(shí)庫(kù)，而本文對(duì)尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)的更新從以下兩個(gè)途徑實(shí)現(xiàn):①對(duì)q代粒子中適應(yīng)度最小的粒子進(jìn)行了一次真實(shí)適應(yīng)度評(píng)價(jià)，不浪費(fèi)此次適應(yīng)度評(píng)價(jià)，將該粒子個(gè)體信息作為新的知識(shí)源替換尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)中最差粒子，更新尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)；②當(dāng)不滿足精度要求進(jìn)入外循環(huán)時(shí)，在最優(yōu)粒子指引下主粒子群結(jié)構(gòu)產(chǎn)生第2+p（p=0,1,…,p）代粒子（p為外循環(huán)次數(shù)），經(jīng)過適應(yīng)度評(píng)價(jià)，將此代粒子信息與原尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)的所有樣本信息按照適應(yīng)度由小到大順序進(jìn)行排列，選擇與原知識(shí)庫(kù)數(shù)量相等的適應(yīng)度較小的粒子組建新的尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)。與原方法相比，通過保持知識(shí)庫(kù)的容量不變并提高知識(shí)庫(kù)質(zhì)量，可以明顯提高GP模型的預(yù)測(cè)能力與精度。

（3）此外，本文方法采用雙層嵌套的粒子群結(jié)構(gòu)，充分利用粒子群自身尋優(yōu)技術(shù)，降低對(duì)粒子群尋優(yōu)結(jié)構(gòu)的破壞，大大提高粒子群獲得最優(yōu)解的機(jī)會(huì)，而原方法采用的搜索空間動(dòng)態(tài)壓縮、全局最優(yōu)解隨機(jī)擾動(dòng)等措施，需要干擾粒子群本身尋優(yōu)，增加了算法的難度。

3.3.2 方法的實(shí)現(xiàn)步驟

本文方法的流程圖見圖 2，其框架主要由兩個(gè)嵌套的粒子群算法組成，主粒子群結(jié)構(gòu)由外層循環(huán)構(gòu)建，內(nèi)層循環(huán)構(gòu)建子粒子群結(jié)構(gòu)。具體實(shí)現(xiàn)步驟如下:

（1）應(yīng)用FLAC3D數(shù)值計(jì)算軟件建立隧洞數(shù)值計(jì)算模型，采用某種巖體流變本構(gòu)模型，將其參數(shù)作為決策變量。編制Fish程序，選取測(cè)點(diǎn)、監(jiān)測(cè)時(shí)間，依據(jù)式（1）建立適應(yīng)度函數(shù)。

（2）啟動(dòng) MATLAB環(huán)境下的優(yōu)化程序，根據(jù)決策變量的個(gè)數(shù)設(shè)置PSO算法的維數(shù)。為了保證隨機(jī)解的互異性，依照粒子群規(guī)則隨機(jī)初始化第1代和第2+p代粒子作為隨機(jī)解（p為動(dòng)態(tài)更新最優(yōu)粒子的循環(huán)次數(shù)）。調(diào)取FLAC3D程序，將隨機(jī)解作為巖體的計(jì)算力學(xué)參數(shù)，通過FLAC3D獲取各測(cè)點(diǎn)的計(jì)算位移，將各測(cè)點(diǎn)的計(jì)算位移代入適應(yīng)度函數(shù)，獲得每個(gè)粒子的適應(yīng)度，然后返回 MATLAB環(huán)境下的優(yōu)化程序，對(duì)隨機(jī)解的適應(yīng)度進(jìn)行比選。如果最小適應(yīng)度滿足精度要求，則退出計(jì)算，最小適應(yīng)度對(duì)應(yīng)粒子即為最優(yōu)解；否則，轉(zhuǎn)入下一步。

（3）當(dāng)步驟（2）中p=0時(shí)，基于第1代和第2代粒子建立尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)；當(dāng)p≠0時(shí)，更新尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)。接著，將知識(shí)庫(kù)中適應(yīng)度最小粒子作為歷史最優(yōu)解。

（4）根據(jù)PSO規(guī)則在內(nèi)循環(huán)產(chǎn)生第q代粒子，采用GP機(jī)器學(xué)習(xí)模型對(duì)知識(shí)庫(kù)中樣本進(jìn)行學(xué)習(xí)，初步建立自變量與函數(shù)值的非線性映射關(guān)系。借助學(xué)習(xí)后的GP模型取代FLAC3D計(jì)算對(duì)此代粒子進(jìn)行適應(yīng)度預(yù)測(cè)，通過q次循環(huán)選取預(yù)測(cè)適應(yīng)度最小粒子作為基于GP模型的最優(yōu)解的預(yù)測(cè)值。

（5）考慮到最優(yōu)解的預(yù)測(cè)值與數(shù)值計(jì)算解之間必然存在一定的誤差，將最優(yōu)解的預(yù)測(cè)值代入FLAC3D進(jìn)行正向計(jì)算，將此解獲得對(duì)應(yīng)的真實(shí)的適應(yīng)度函數(shù)值返回 MATLAB環(huán)境下，并與當(dāng)前迭代步歷史最優(yōu)解的適應(yīng)度進(jìn)行比較，較小者為當(dāng)前迭代步對(duì)應(yīng)的當(dāng)前最優(yōu)解。將當(dāng)前最優(yōu)解返回尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)，動(dòng)態(tài)更新尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)。

（6）對(duì)當(dāng)前最優(yōu)解進(jìn)行精度判斷，如果滿足精度要求，則退出計(jì)算；反之，則動(dòng)態(tài)更新最優(yōu)粒子，循環(huán)次數(shù)p增加一次，指導(dǎo)主粒子群結(jié)構(gòu)產(chǎn)生下一代粒子，引領(lǐng)粒子群飛向最有前景的尋優(yōu)方向。同時(shí)，將此新一代粒子經(jīng)過適應(yīng)度評(píng)價(jià)產(chǎn)生的新信息與原尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)的所有樣本信息，按照適應(yīng)度由小到大順序進(jìn)行篩選，選取與原知識(shí)庫(kù)中同樣數(shù)目的樣本更新知識(shí)庫(kù)。按照上述兩種策略動(dòng)態(tài)更新尋優(yōu)經(jīng)驗(yàn)知識(shí)庫(kù)，克服隨機(jī)優(yōu)化算法過度依賴初始學(xué)習(xí)樣本這一局限性，保證算法的全局性。

依據(jù)上述步驟編制了MATLAB計(jì)算程序，程序中采用接口程序直接調(diào)用FLAC3D數(shù)值計(jì)算程序。

圖2 PSO-GP-FLAC3D方法流程圖Fig.2 Flow chart of PSO-GP-FLAC3D method

4 算例研究

為了驗(yàn)證上述方法的可行性，分別采用標(biāo)準(zhǔn)PSO-FLAC3D與PSO-GP-FLAC3D方法求取隧洞圍巖損失位移。取1/4圓形隧洞如圖3所示，隧洞直徑為20 m，隧洞沿坐標(biāo)軸x、z方向圍巖厚度均為90 m，沿 y方向取單位厚度，均質(zhì)巖性，初始地應(yīng)力為σx=σy=σz=30 MPa。本構(gòu)模型采用經(jīng)典黏彈性本構(gòu)模型——Viscous模型，決策變量為泊松比μ、彈性模量E和黏滯系數(shù)η共3個(gè)變量，則PSO的維數(shù)為3維。作為例證，沿環(huán)向分別在 30°、45°、60°處選取3個(gè)測(cè)點(diǎn)a、b、c作為監(jiān)測(cè)點(diǎn)（見圖4），真解設(shè)為:μ=0.2，E=2.0 GPa，η=10 GPa·d，總的計(jì)算時(shí)間為10 d，滯后監(jiān)測(cè)為4 d，即實(shí)際監(jiān)測(cè)總時(shí)間為6 d。將后6 d測(cè)點(diǎn)計(jì)算位移作為實(shí)測(cè)位移，采用本文方法預(yù)測(cè)前4 d測(cè)點(diǎn)已發(fā)生的損失位移。

目標(biāo)函數(shù)即適應(yīng)度函數(shù)參見式（1），其中:監(jiān)測(cè)總時(shí)間m=6 d；位移測(cè)點(diǎn)總數(shù)n=3；兩種方法的搜索范圍設(shè)為:0.1≤μ≤0.5，1.0 GPa≤E≤5.0 GPa，8.0 GPa·d≤η≤12.0 GPa·d。

圖3 隧洞網(wǎng)格劃分圖Fig.3 Mesh of tunnel

圖4 隧洞開挖面的測(cè)點(diǎn)布置圖Fig.4 Monitoring points of excavation section for tunnel

以精度滿足f＜1×10-4作為終止條件，標(biāo)準(zhǔn)PSO算法參數(shù)設(shè)置:維數(shù) D=3，種群個(gè)數(shù) NP=30，c1=c2=2.0，p=500，最大飛行速度Vmax= [0.01,0.01,0.01]。PSO-GP-FLAC3D模型中PSO的參數(shù)設(shè)置同標(biāo)準(zhǔn)PSO算法，q =10，GP模型初始超參數(shù)lnl =[-1,-1,-1]，lσnf=-1，σn=1×10-6。

為了保證兩種算法對(duì)比的公平性，每種算法都獨(dú)立運(yùn)行5次，取其適應(yīng)度評(píng)價(jià)次數(shù)的平均值作為該算法的最終適應(yīng)度評(píng)價(jià)次數(shù)。計(jì)算結(jié)果見表 1，兩種方法的相對(duì)誤差令人滿意。但與標(biāo)準(zhǔn)PSO方法相比，采用PSO-GP-FLAC3D方法時(shí)計(jì)算代價(jià)得到了大幅度削減，適應(yīng)度評(píng)價(jià)次數(shù)和計(jì)算耗時(shí)僅為標(biāo)準(zhǔn)PSO-FLAC3D的1/5。

表1 兩種方法計(jì)算結(jié)果比較Table 1 Results comparison with different methods

PSO-GP-FLAC3D方法在3個(gè)測(cè)點(diǎn)處與實(shí)測(cè)位移曲線的擬合情況見圖5。從圖中可見，在3個(gè)測(cè)點(diǎn)處計(jì)算監(jiān)測(cè)收斂位移曲線與實(shí)際監(jiān)測(cè)收斂位移曲線基本吻合。最終計(jì)算損失位移與實(shí)測(cè)損失位移大小比較如表 2所示，3個(gè)測(cè)點(diǎn)中最大相對(duì)誤差為0.37%，平均相對(duì)誤差僅有0.32%，說明本文方法是可行的，能夠通過對(duì)實(shí)際監(jiān)測(cè)收斂位移曲線的準(zhǔn)確擬合，可獲取合理的損失位移。

圖5 各測(cè)點(diǎn)的計(jì)算位移與實(shí)測(cè)位移對(duì)比Fig.5 Comparison between calculated displacement and measured displacement of monitoring points

表2 計(jì)算損失位移與實(shí)測(cè)損失位移比較Table 2 Comparison between calculated loss displacement and measured loss displacement

為了探討滯后監(jiān)測(cè)時(shí)間對(duì)損失位移預(yù)測(cè)精度的影響，本文還對(duì)滯后監(jiān)測(cè)天數(shù)為6、8 d情況下不同測(cè)點(diǎn)的損失位移進(jìn)行預(yù)測(cè)，計(jì)算結(jié)果如表3所示。從中可見，滯后監(jiān)測(cè)天數(shù)越大，反分析獲得的損失位移精度越低。當(dāng)滯后監(jiān)測(cè)天數(shù)達(dá)到8 d時(shí)，3個(gè)測(cè)點(diǎn)獲取的損失位移相對(duì)誤差均超過了 5%，計(jì)算結(jié)果可信度相對(duì)較低。

5 工程應(yīng)用

錦屏二級(jí)水電站位于四川省涼山彝族自治州境內(nèi)的雅礱江錦屏大河彎處雅礱江干流上，系利用雅礱江錦屏150 km長(zhǎng)大河彎的天然落差，裁彎取直、鑿洞引水建成的。錦屏二級(jí)水電站共有引水隧洞 4條，最大埋深為2550 m，為超深埋長(zhǎng)隧洞特大型地下水電工程。為了保證隧洞施工的安全性，開挖輔助洞A、B以超前探明地質(zhì)情況[15]。輔助洞大部分洞段處高地應(yīng)力區(qū)，開挖后圍巖應(yīng)力釋放劇烈。下面以輔助洞B的BK14+599斷面為例，利用實(shí)測(cè)位移資料，采用本文方法進(jìn)行圍巖損失位移的優(yōu)化反分析。

表3 不同監(jiān)測(cè)天數(shù)損失位移計(jì)算結(jié)果比較Table 3 Results comparison of calculated loss displacement for different monitoring days

圖6 BK14+599斷面的AC測(cè)線布置圖Fig.6 AC monitoring line of section BK14+599

數(shù)值計(jì)算網(wǎng)格劃分如圖7所示，模型的水平x方向與豎直z方向均取60 m，沿隧洞軸線y方向取單位厚度。

圖7 BK14+599斷面的網(wǎng)格劃分圖Fig.7 Mesh of section BK14+599

選取AC測(cè)線的實(shí)測(cè)收斂位移，依據(jù)式（1）建立適應(yīng)度函數(shù)，其中:監(jiān)測(cè)總時(shí)間m=45 d；位移測(cè)點(diǎn)總數(shù) n=2；決策變量的搜索范圍設(shè)為:0.1≤μ≤0.3，1.0 GPa≤E≤12.0 GPa，1.0×104GPa·d≤η≤12.0×104GPa·d。采用本文方法進(jìn)行求解，相關(guān)參數(shù)設(shè)置與以上算例相同。

經(jīng)優(yōu)化反分析，獲得的全局最優(yōu)解是:圍巖泊松比μ=0.263、圍巖彈性模量E =10.295 GPa、圍巖滯系數(shù)η=79567.5 GPa·d。測(cè)線AC的損失位移達(dá)10.45 mm（見圖8），占實(shí)際總位移的比例高達(dá)36%。如果僅依據(jù)遠(yuǎn)小于實(shí)際總位移的監(jiān)測(cè)位移來指導(dǎo)隧洞開挖或評(píng)價(jià)圍巖穩(wěn)定性，將給該隧洞的施工和運(yùn)行安全帶來潛在隱患。因此，該工程由于埋深較大，開挖后的彈性變形大，其引水隧洞及輔助洞的圍巖損失位移不可忽視，在圍巖穩(wěn)定性評(píng)價(jià)中應(yīng)給予足夠重視。

圖8 測(cè)線AC的計(jì)算位移與實(shí)測(cè)位移對(duì)比Fig.8 Comparison between calculated displacement and measured displacement of AC monitoring line

6 結(jié) 論

（1）損失位移的準(zhǔn)確預(yù)測(cè)對(duì)于地下工程圍巖穩(wěn)定性的合理評(píng)價(jià)具有重要意義。本文提出了1種損失位移求解的新思路，即根據(jù)滯后監(jiān)測(cè)的圍巖位移信息，通過優(yōu)化反分析方法預(yù)測(cè)監(jiān)測(cè)時(shí)刻前已發(fā)生的損失位移，并可同時(shí)獲取合理的圍巖計(jì)算力學(xué)參數(shù)。研究表明，該思路是可行的，為圍巖損失位移問題的求解提供了一條新的途徑。

（2）針對(duì)位移全局優(yōu)化反分析過程計(jì)算量較大的問題，本文結(jié)合數(shù)值計(jì)算模型，提出了一種基于高斯過程機(jī)器學(xué)習(xí)與粒子群優(yōu)化算法的損失位移智能優(yōu)化反分析方法，算例研究與工程應(yīng)用表明，該方法是可行的，與傳統(tǒng)優(yōu)化反分析方法和以往智能優(yōu)化反分析方法相比較，該方法采用動(dòng)態(tài)更新學(xué)習(xí)樣本的策略，具有不依賴初始學(xué)習(xí)樣本的優(yōu)點(diǎn)，在獲取全局最優(yōu)解的前提下，計(jì)算效率具有一定的優(yōu)越性。

（3）通過算例驗(yàn)證，隨著監(jiān)測(cè)天數(shù)的減少，獲取的損失位移精度降低。在實(shí)際應(yīng)用中，應(yīng)盡可能早的布置監(jiān)測(cè)儀器，盡可能多的獲取監(jiān)測(cè)位移數(shù)據(jù)，將有利于進(jìn)一步提高所獲取的損失位移的可靠性。算例研究表明，隨著滯后監(jiān)測(cè)天數(shù)的增加，采用本文方法獲取的損失位移的誤差也會(huì)隨著增高。因此，在本文方法的工程應(yīng)用中，應(yīng)注意監(jiān)測(cè)滯后天數(shù)對(duì)預(yù)測(cè)精度的影響，滯后天數(shù)較大時(shí)，應(yīng)適當(dāng)調(diào)整對(duì)預(yù)測(cè)結(jié)果的可靠度的預(yù)期。另一方面，在開挖后及時(shí)布置監(jiān)測(cè)儀器，盡量減小損失位移，將有利于提高圍巖損失位移與力學(xué)參數(shù)反分析結(jié)果的可靠性。

（4）錦屏二級(jí)水電站輔助洞某斷面的損失位移反分析結(jié)果表明，損失位移占實(shí)際發(fā)生總位移的比例高達(dá)36%，說明深部地下巖體工程中，受高地應(yīng)力的影響，開挖后圍巖在短時(shí)間內(nèi)彈性變形大，圍巖損失位移不可忽視，在圍巖穩(wěn)定性評(píng)價(jià)與圍巖參數(shù)反分析中應(yīng)給予足夠重視。

（5）為便于工程應(yīng)用，本文僅對(duì)基于經(jīng)典黏彈性本構(gòu)模型——Viscous模型的二維空間圍巖損失位移求取方法進(jìn)行了探討。在工程應(yīng)用中，采用本文方法時(shí)應(yīng)考慮工程圍巖特點(diǎn)，選擇合適的巖體本構(gòu)模型。此外，隧道圍巖變形還受到開挖方式、開挖進(jìn)尺、支護(hù)等多種因素的綜合影響。關(guān)于多種因素影響下三維空間圍巖損失位移反分析問題的研究將另文闡述。

[1]楊林德. 巖土工程問題的反演理論與工程實(shí)踐[M]. 北京:科學(xué)出版社,1996.

[2]蔡美峰,何滿潮,劉東燕. 巖石力學(xué)與工程[M]. 北京:科學(xué)出版社,2002.

[3]楊志法,熊順成,王存玉,等. 關(guān)于位移反分析的某些考慮[J]. 巖石力學(xué)與工程學(xué)報(bào),1995,14(1):11－16.YANG Zhi-fa,XIONG Shun-cheng,WANG Cun-yu,et al.Some consideration of the back-analysis from displacements[J]. Chinese Journal of Rock Mechanics and Engineering,1995,14(1):11－16.

[4]張傳慶,馮夏庭,周輝,等. 隧洞圍巖收斂損失位移的求取方法及應(yīng)用[J]. 巖土力學(xué),2009,30(4):997－1003.ZHANG Chuan-qing,FENG Xia-ting,ZHOU Hui,et al.Method of obtaining loss convergence displacement and its application to tunnel engineering[J]. Rock and Soil Mechanics,2009,30(4):997－1003.

[5]高瑋,鄭穎人. 采用快速遺傳算法進(jìn)行巖土工程反分析[J]. 巖土工程學(xué)報(bào),2001,23(1):120－122.GAO Wei,ZHENG Ying-ren. Back analysis in geotechnical engineering based on fast-convergent genetic algorithm[J]. Chinese Journal of Geotechnical Engineering,2001,23(1):120－122.

[6]趙洪波. 基于微粒群優(yōu)化的智能位移反分析研究[J].巖土工程學(xué)報(bào),2006,28(11):2035－2038.ZHAO Hong-bo. Back analysis of intelligent displacement based on particle swarm optimization[J].Chinese Journal of Geotechnical Engineering,2006,28(11):2035－2038.

[7]李麗,牛奔. 粒子群優(yōu)化算法[M]. 北京:冶金工業(yè)出版社,2009.

[8]KENNEDY J,EBERHART R. Particle swarm optimization[C]//Proceedings of IEEE International Conference on Neural Networks. Perth:IEEE Press,1995.

[9]SEEGER M. Gaussian processes for machine learning[J].International Journal of Neural System,2004,14(4):69－106.

[10]SOFIANE B B,AMINE B. Gaussian process for nonstationary times series prediction[J]. Computational Statistics &Data Analysis,2004,47(4):705－712.

[11]彭文斌. FLAC3D實(shí)用教程[M]. 北京:機(jī)械工業(yè)出版社,2008.

[12]蘇國(guó)韶,宋詠春,燕柳斌. 巖體爆破效應(yīng)預(yù)測(cè)的一種新方法[J]. 巖石力學(xué)與工程學(xué)報(bào),2007,26(增刊):3509－3514.SU Guo-shao,SONG Yong-chun,YAN Liu-bin. A new method for forecasting of blasting effect in rock mass[J].Chinese Journal of Rock Mechanics and Engineering,2007,26(Supp.):3509－3514.

[13]蘇國(guó)韶,燕柳斌,張小飛,等. 基坑位移時(shí)間序列預(yù)測(cè)的高斯過程方法[J]. 廣西大學(xué)學(xué)報(bào)(自然科學(xué)版),2007,32(2):223－226.SU Guo-shao,YAN Liu-bin,ZHANG Xiao-fei,et al.Time series prediction of foundation pit displacement using Gaussian process method[J]. Journal of Guangxi University(Natural Science Edition),2007,32(2):223－226.

[14]蘇國(guó)韶,張克實(shí),呂海波. 位移反分析的粒子群優(yōu)化-高斯過程協(xié)同優(yōu)化方法[J]. 巖土力學(xué),2011,32(2):510－515,524.SU Guo-shao,ZHANG Ke-shi,Lü Hai-bo. A cooperative optimization method based on particle swarm optimization and Gaussian process for displacement back analysis[J]. Rock and Soil Mechanics,2011,32(2):510－515,524.

[15]劉建友,伍法權(quán),盧丙清,等. 雅礱江錦屏二級(jí)水電站皮帶輸送隧洞施工中的地質(zhì)問題分析及其處理措施[J].工程地質(zhì)學(xué)報(bào),2009,17(5):591－596.LIU Jian-you,WU Fa-quan,LU Bing-qing,et al.Geological problems encountered in rock tunnel at Jinping Ⅱ hydroelectric power station and their treatment measures[J]. Journal of Engineering Geology,2009,17(5):591－596.

[16]陳炳瑞,馮夏庭,黃書嶺,等. 基于快速拉格朗日分析——并行粒子群算法的黏彈塑性參數(shù)反演及其應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報(bào),2007,26(12):2517－2525.CHEN Bing-rui,FENG Xia-ting,HUANG Shu-ling,et al.Inversion of viscoelasto-plastic parameters based on fast Lagrangian analysis of continuum—parallel particle swarm algorithm and its application[J]. Chinese Journal of Rock Mechanics and Engineering,2007,26(12):2517－2525.