克里格法對(duì)隧道圍巖穩(wěn)定性分析

薛洪濤

(中鐵十九局集團(tuán)第七工程有限公司，廣西 百色 533000)

?

克里格法對(duì)隧道圍巖穩(wěn)定性分析

薛洪濤

(中鐵十九局集團(tuán)第七工程有限公司，廣西 百色 533000)

隧道圍巖穩(wěn)定性一直是隧道施工控制的重點(diǎn)，采用克里格算法對(duì)隧道襯砌結(jié)構(gòu)穩(wěn)定性計(jì)算進(jìn)行了研究。根據(jù)不同的變異函數(shù)對(duì)隨機(jī)變量特征進(jìn)行表達(dá)，結(jié)合二次正交組合法，建立了隧道襯砌結(jié)構(gòu)功能函數(shù)的插值算法，避免了隱式函數(shù)無(wú)法直接求解的問(wèn)題；對(duì)隧道力學(xué)參數(shù)的取值,提出了基于超前地質(zhì)預(yù)報(bào)中的地震波法獲取方法，克服了隧道力學(xué)參數(shù)取值的模糊性。采用克里格法與采用傳統(tǒng)蒙特卡洛方法的對(duì)比結(jié)果表明：迭代次數(shù)減少,失效概率的絕對(duì)誤差僅為0.0152%，相對(duì)誤差為4.26%。證明了該方法的合理性。

克里格法；隧道；二次正交組合法；超前地質(zhì)預(yù)報(bào)

引 言

隧道作為很多公路或鐵路的控制性工程，其圍巖穩(wěn)定性是重點(diǎn)，一旦出現(xiàn)工程事將損失巨大。研究者采用理論方法、現(xiàn)場(chǎng)監(jiān)測(cè)和數(shù)值模擬等多種手段對(duì)隧道圍巖穩(wěn)定性進(jìn)行了大量研究，得到的成果頗多。對(duì)于理論算法，1951年，南非地質(zhì)學(xué)者Danie Krige提出了克里格模型[1-2]。該模型采用隨機(jī)過(guò)程算法，能夠得到方差估計(jì)較小的無(wú)偏估計(jì)，且具有半?yún)?shù)化的特點(diǎn)[3-4]，本文擬采用該方法應(yīng)用于隧道圍巖穩(wěn)定性分析。

對(duì)于隧道力學(xué)參數(shù)的取值，計(jì)算時(shí)往往采用室外取試樣后在室內(nèi)做試驗(yàn)確定，得到的參數(shù)往往存在誤差。而且隧道賦存的巖土環(huán)境本來(lái)就復(fù)雜多樣，隧道的力學(xué)參數(shù)在同一斷面的不同位置及不同斷面的相同位置，力學(xué)參數(shù)都是變化的。隨著隧道超前地質(zhì)預(yù)報(bào)的發(fā)展，出現(xiàn)了一系列新技術(shù)，比如地質(zhì)雷達(dá)、鉆孔攝像和地震波法等手段。特別是地震波法，即采用瑞士的TSP203 plus，可獲取測(cè)試掌子面前方巖體的泊松比μ、內(nèi)摩擦角φ、彈性模量E和粘聚力c等一系列力學(xué)參數(shù)，其他參數(shù)也可經(jīng)過(guò)計(jì)算得到[5]。這為隧道圍巖穩(wěn)定計(jì)算提供了新思路，即根據(jù)獲得的實(shí)測(cè)數(shù)據(jù)分析圍巖穩(wěn)定，有更高的精準(zhǔn)性與實(shí)時(shí)性[6]。

根據(jù)不同的變異函數(shù)對(duì)隨機(jī)變量特征進(jìn)行表達(dá)，結(jié)合二次正交組合法，建立了隧道襯砌結(jié)構(gòu)功能函數(shù)的插值算法，避免了隱式函數(shù)無(wú)法直接求解的問(wèn)題，超前地質(zhì)預(yù)報(bào)中的地震波法獲取隧道力學(xué)參數(shù)對(duì)隧道襯砌可靠度進(jìn)行計(jì)算。

1 襯砌結(jié)構(gòu)穩(wěn)定功能函數(shù)

隧道結(jié)構(gòu)主動(dòng)壓力經(jīng)歷彈性階段、彈塑性階段、塑性-分離階段和松動(dòng)階段四個(gè)變形階段，如圖1所示。實(shí)踐表明，D點(diǎn)是支護(hù)的最佳時(shí)刻，此時(shí)圍巖壓力PD達(dá)到最小，用Pamin表示。

圖1 圍巖形變與壓力發(fā)展過(guò)程

Pamin= γr0[Rmax/r0-1]

(1)

式中，Rmax為最大的圍巖允許松動(dòng)區(qū)半徑，γ為巖體重度，r0為隧道等代圓半徑。

根據(jù)卡斯特納方程，圍巖最大松動(dòng)圈半徑Rmax：

(2)

式中，c1、φ分別為已開(kāi)挖破裂區(qū)圍巖的粘聚力和內(nèi)摩擦角，σz為側(cè)壓系數(shù)λ=1的原巖應(yīng)力。

c1=c+τaAs/SaSb

(3)

式中：c為圍巖本身粘聚力，τa為錨桿抗剪強(qiáng)度，As為錨桿的橫截面積，Sa,Sb分別為錨桿的橫向和縱向間距。

根據(jù)可靠度有關(guān)定義，則極限狀態(tài)方程為：

Z=Pi-Pamin=0

(4)

式中，Pi為支護(hù)結(jié)構(gòu)即襯砌提供的支護(hù)阻力。

由式(1)～式(3)可知，Pamin為γ、r0、φ、c、σz、Sa、Sb、Pamin的函數(shù)，即

Pamin=f(γ,r0,φ,c,σz,Sa,Sb,Pamin)

(5)

Pamin為隱式函數(shù)，在確定性分析中，只能采用迭代方法求解。另外，由于式(5)不能表達(dá)為基本參數(shù)的明晰解析式，也就不可獲取Pamin的概率分布，從而無(wú)法采用直接概率積分及已有的可靠度計(jì)算方法，必須尋找等效計(jì)算方法。本文引入克里格插值模型解決這一問(wèn)題。

2 克里格插值算法

2.1 變異函數(shù)

設(shè)Z(x)為一隨機(jī)函數(shù)，x∈Rn,x1,x2, ... ,xm是n維空間中樣本點(diǎn)的位置，Z(x1)，Z(x2)…，Z(xm)為相應(yīng)的聯(lián)合分布的隨機(jī)變量?？死锔窀鶕?jù)隨機(jī)變量變異的空間位置特征，提出采用變異函數(shù)γ(h)進(jìn)行模擬。在有限樣本的情況下，γ(h)值可以計(jì)算：

(6)

Nh指在(xi+h,xi)之間用來(lái)計(jì)算n維空間中樣本的變異函數(shù)值的樣本對(duì)數(shù)，下標(biāo)h表示Nh是分離距離h的函數(shù)；Z(xi+h)指點(diǎn)xi偏離h處隨機(jī)變量的實(shí)現(xiàn)值；Z(xi) 指點(diǎn)xi處隨機(jī)變量的實(shí)現(xiàn)值；h為步長(zhǎng)，為樣本點(diǎn)的空間間隔距離。

變異函數(shù)可用于表征隨機(jī)變量的空間變異結(jié)構(gòu)或空間連續(xù)性，目前比較常用的有球形模型、指數(shù)模型、高斯模型和冪函數(shù)模型等。采用式(6)計(jì)算得出h及γ(h)值，做出h-γ(h)試驗(yàn)變異圖。在試驗(yàn)變異圖的基礎(chǔ)上，再對(duì)其擬合一個(gè)最優(yōu)的理論變異函數(shù)模型。

2.2 算法

普通克里格法估計(jì)的一般公式：

(7)

式中，Z(xi)為隨機(jī)函數(shù)在xi位置的已知值；Z*(xi)為在x0位置的隨機(jī)函數(shù)估計(jì)值；n為估計(jì)過(guò)程中已知值的個(gè)數(shù)；λi為分配給Z(xi)殘差的權(quán)重。

使用式(6)進(jìn)行插值估計(jì)時(shí)，克里格法需滿足2個(gè)要求：(1)選取Ki，使Z*(x0)的估計(jì)無(wú)偏；(2)估計(jì)誤差的方差σz最小。

根據(jù)要求(1)有：

E[Z*(x0)-Z(X0)]=0

(8)

將式(8)代入式(7)，經(jīng)整理可得：

(9)

根據(jù)要求(2)有：

σ2=Pamin=g(γ,r0,φ,c,σz,Sa,Sb)

(10)

將式(7)代入(8)，簡(jiǎn)化得：

(11)

式中，γ(h)為變異函數(shù)，由式(7)可得λi，根據(jù)拉格朗日原理和駐點(diǎn)的一階偏導(dǎo)數(shù)原理，對(duì)式(11)的極值運(yùn)算可得：

(12)

式中，μ為拉格朗日乘數(shù)，λi滿足：

(13)

由式(12)、式(13)得：

(14)

式中，γij=γ(xi-xj)即距離為xi和xj之間的變異函數(shù)值，式(14)可簡(jiǎn)化為

[F][L]=[B]

(15)

式中，[F]為式(14)左邊的系數(shù)矩陣，[L]為式(14)左邊的待求解的向量，[B]為式(14)右邊的已知向量。解這些線性方程組，即可得到所有的權(quán)重1,..., λn和拉格朗日乘數(shù)μ。

由計(jì)算所得的權(quán)重和拉格朗日乘數(shù)，克里格法的估計(jì)方差S可通過(guò)式(7)求得估計(jì)值Z*計(jì)算：

(16)

3 結(jié)合二次正交法的克里格插值可靠度計(jì)算

3.1 試驗(yàn)設(shè)計(jì)抽樣方法

二次回歸方程的一般形式為：

k=1，2，…，m-1(j≠k)

(17)

其中，a,{bj},{bkj},{bjj}為回歸系數(shù)，可以看出該方程共有 (m+1)(m+2)/2項(xiàng)，因此試驗(yàn)次數(shù)n≥(m+1)(m+2)/2。

二次回歸正交組合設(shè)計(jì)的總試驗(yàn)次數(shù)為：

n=mc+mγ+m0

(18)

根據(jù)正交性的要求，可以推導(dǎo)出星號(hào)臂長(zhǎng)度γ必須滿足關(guān)系式：

(19)

可見(jiàn)，星號(hào)臂長(zhǎng)度γ表征試驗(yàn)區(qū)間，與因素?cái)?shù)m、零水平試驗(yàn)次數(shù)m0及二水平試驗(yàn)數(shù)mc有關(guān)。

3.2 確定合適的二次回歸正交組合設(shè)計(jì)

根據(jù)因素m選擇合適的正交表進(jìn)行交換，明確二水平試驗(yàn)方案，確定二水平試驗(yàn)次數(shù)mc和星號(hào)試驗(yàn)次數(shù)mγ。

根據(jù)二次回歸正交組合設(shè)計(jì)表確定的試驗(yàn)方案，共進(jìn)行n次試驗(yàn)，得到n個(gè)試驗(yàn)指標(biāo)。有多少個(gè)未知數(shù)，便至少需進(jìn)行多少次試驗(yàn)計(jì)算回歸系數(shù)，建立含規(guī)范變量的回歸方程。回歸系數(shù)的計(jì)算：

常數(shù)項(xiàng)：

一次項(xiàng)偏回歸系數(shù)：

交互項(xiàng)偏回歸系數(shù)：

二次項(xiàng)偏回歸系數(shù)：

總平方和：

自由度：dfT=n-1

一次項(xiàng)偏回歸平方和：

交互項(xiàng)偏回歸平方和：

二次項(xiàng)偏回歸平方和：

各種偏回歸平方和的自由度都為1。

殘差平方和：

SSe=SST-SSR

自由度：

dfe=dfT-dfR

4 可靠性計(jì)算

4.1 TSP軟件獲取巖石力學(xué)參數(shù)

TSP203探測(cè)是根據(jù)人工制造一系列地震波的回波原理，在以開(kāi)挖的隧道邊墻布置若干個(gè)炮眼并置炸藥半截與孔內(nèi)，通過(guò)瞬發(fā)電雷管引爆每個(gè)藥孔形成一系列的輕微震源，地震波經(jīng)過(guò)掌子面前方返回被接收元件接受并傳輸于電腦。經(jīng)處理分析軟件可分析出掌子面后未開(kāi)挖不良地質(zhì)體，如：軟弱破碎帶、溶洞和溝槽等不良地質(zhì)構(gòu)造，同樣可測(cè)掌子面為參考，根據(jù)波速、動(dòng)態(tài)模量等一系列巖體力學(xué)參數(shù)判斷將開(kāi)挖里程段的圍巖級(jí)別[5-6]。

提取的反射層可以顯示小范圍內(nèi)巖石構(gòu)造的顯著變化。巖石特性的大幅度變化(與地震的主波長(zhǎng)相當(dāng)，主波長(zhǎng)即波速除以主頻)將通過(guò)速度函數(shù)，以較簡(jiǎn)單、可靠的方式進(jìn)行評(píng)估。TSPwin通過(guò)縱橫波的速度計(jì)算巖石力學(xué)參數(shù)，對(duì)于某些參數(shù)，則必須借助于經(jīng)驗(yàn)關(guān)系式[7-8]?？蓮乃俣绕拭嬷苯佑?jì)算出縱橫波速度比和泊松比，TSP win中的巖石密度ρ通過(guò)縱橫波速度的經(jīng)驗(yàn)公式計(jì)算，求得密度后，可以計(jì)算出體積模量k、剪切模量μ、拉梅參數(shù)λ以及動(dòng)態(tài)楊氏模量Edyn，縱橫波速度比和泊松比是無(wú)量剛的，其他參數(shù)可由公式推導(dǎo)[6-8]。對(duì)于變質(zhì)巖、火山巖、火成巖、沉積巖四大巖類，不同類型的巖石使用不同的經(jīng)驗(yàn)公式。對(duì)于所有的由經(jīng)驗(yàn)公式計(jì)算出的彈性參數(shù)，取由兩種公式所求得結(jié)果的平均值。

計(jì)算過(guò)程中，不能由TSP獲得的參數(shù)，即粘聚力、內(nèi)摩擦角、錨桿的橫向和縱向間距，需要根據(jù)地質(zhì)資料、圍巖分類和反分析來(lái)確定力學(xué)參數(shù)。等代圓半徑r0變異性不大，取隧道等代圓半徑即可，而σz原巖應(yīng)力由隧道埋深h計(jì)算，參數(shù)變化見(jiàn)表1。

表1 圍巖基本參數(shù)統(tǒng)計(jì)特征

4.2 克里格法的實(shí)現(xiàn)過(guò)程

計(jì)算方法1：

計(jì)算方法2：

選取普通克里格方法建立功能函數(shù)Pamin的插值模型；隨機(jī)從隨機(jī)變量(x1,x2, ... ,xm)中抽取一組變量代入Pamin的插值模型中計(jì)算對(duì)應(yīng)的Pamin值，再將Pamin的值代入極限狀態(tài)方程中求出功能函數(shù)Z值，重復(fù)n次，即可獲得ζ1,ζ2...ζn的n個(gè)值計(jì)算出其均值和方差。

(20)

(21)

由可靠度定義式求出可靠度指標(biāo)。

因?yàn)镻amin是關(guān)于γ、r0、φ、c、σz、Sa、Sb、Pamin的函數(shù)，一共含有7個(gè)隨機(jī)變量。若在計(jì)算時(shí)同時(shí)考慮每一個(gè)變量，則計(jì)算工作量較大，過(guò)程繁瑣。在滿足實(shí)際工程需要的前提下，可以從中選取若干因素，將其視為隨機(jī)變量進(jìn)行分析。在利用TSP系統(tǒng)進(jìn)行地質(zhì)預(yù)報(bào)的過(guò)程中，可以獲得巖體密度、泊松比以及靜和動(dòng)態(tài)彈性模量等參數(shù)，部分其他參數(shù)取自地質(zhì)資料、公式計(jì)算和反分析，其中可以用作可靠度計(jì)算的參數(shù)為巖體密度γ(kN/m3)。

運(yùn)用二次回歸正交試驗(yàn)設(shè)計(jì)方法和正交試驗(yàn)設(shè)計(jì)方法，對(duì)將巖體密度視為隨機(jī)變量γ取60個(gè)水平進(jìn)行抽樣做可靠度分析[11]，并將結(jié)果進(jìn)行比選。

二次回歸正交方法抽樣的樣本方差為：

Var(f(x))=0.4743

正交試驗(yàn)設(shè)計(jì)方法抽樣的樣本方差為：

Var(f(x))=1.3743

二次回歸正交方法能夠取得較小的樣本方差，因此更適合本文分析。

用二次回歸正交方法分別獲取60、80、100、120個(gè)樣本點(diǎn)計(jì)算結(jié)果。

由表2可以看出，當(dāng)樣本點(diǎn)分別為20、30、40、50與60時(shí)，其失效概率及可靠度的計(jì)算結(jié)果大體相近。為了驗(yàn)證克里格法的正確性，將60個(gè)樣本點(diǎn)的計(jì)算結(jié)果與蒙特卡洛算法的結(jié)果進(jìn)行比較，其中克里格法與蒙特卡洛法的誤差為：

(22)

相對(duì)誤差為：

(23)

其中，Pfk為克里格法的誤差，Pfm為蒙特卡洛法的誤差。由計(jì)算結(jié)果可知選用克里格法進(jìn)行可靠度分析與已有方法結(jié)果較為接近。

表2 不同樣本點(diǎn)時(shí)的可靠度比較

5 結(jié)束語(yǔ)

在隧道結(jié)構(gòu)穩(wěn)定功能函數(shù)的基礎(chǔ)上，引入了克里格插值算法，并結(jié)合二次回歸正交組合統(tǒng)計(jì)，建立了基于克里格算法的隧道圍巖穩(wěn)定性評(píng)價(jià)方法。

利用TSP獲得的圍巖力學(xué)參數(shù)，對(duì)實(shí)例圍巖穩(wěn)定性進(jìn)行了可靠度分析。結(jié)果表明，當(dāng)樣本點(diǎn)個(gè)數(shù)為20、30、40、50、60時(shí)，其圍巖失效概率與可靠度結(jié)果較為接近。同蒙特卡洛法進(jìn)行比較時(shí)，其誤差在5%之內(nèi)，表明采用克里格法進(jìn)行圍巖可靠度計(jì)算的方法是可靠的。

[1] 蘇永華,羅正東,張盼鳳,等.基于Kriging的邊坡穩(wěn)定可靠度主動(dòng)搜索法[J].巖土工程學(xué)報(bào),2013(10):1863-1869.

[2] 靳國(guó)棟,劉衍聰,牛文杰.距離加權(quán)反比插值法和克里格插值法的比較[J].長(zhǎng)春工業(yè)大學(xué)學(xué)報(bào),2003,24(3):53-57.

[3] 張團(tuán)峰,王家華.試論克里格估計(jì)與隨機(jī)模擬的本質(zhì)區(qū)別[J].西安石油學(xué)院學(xué)報(bào),1997,12(2):52-55.

[4] 張崎,李興斯.基于Kriging模型的結(jié)構(gòu)可靠性分析[J].計(jì)算力學(xué)學(xué)報(bào),2006,23(2):175-179.

[5] 張平松,吳健生.中國(guó)隧道及井巷地震波法超前探測(cè)技術(shù)研究分析[J].地球科學(xué)進(jìn)展,2006,21(10):1033-1038.

[6] 楊正華,張文波,王衛(wèi)東.淺析地震波法用于隧道病害的診斷與預(yù)測(cè)[J].災(zāi)害學(xué),2008(1):27-30.

[7] 李堅(jiān).TSP 法在鐵路客運(yùn)專線隧道超前地質(zhì)預(yù)報(bào)工作中的應(yīng)用前景[J].鐵道勘察,2006,31(6):45-49.

[8] 李華,李富,魯光銀.TSP 法與探地雷達(dá)相結(jié)合在隧道超前地質(zhì)預(yù)報(bào)中的應(yīng)用研究[J].工程勘察,2009,37(7):86-90.

[9] BOUCNEAU G,MEIRVENNE V M,Thas O,et al.Integrating properties of soil map delineations into ordinary kriging[J].European Journal of Soil Science,1998,49:213-229.

[10] 張仁鐸.空間變異理論及應(yīng)用[M].北京:科學(xué)出版社,2005.

[11] REZAEE H,ASGHARI O,YAMAMOTO J K.On the reduction of the ordinary kriging smoothing effect[J].Journal of Mining and Environment,2011,2(2):102-117.

Analysis of the Stability of Surrounding Rock in Tunnel Based on Kriging Method

XUEHongtao

(No.7 Engineering China Rail Way 19th Bureau Group, Baise 533000, China)

The stability of surrounding rock in tunnel is always the key point during construction, which was studied through the application of Kriging algorithm. The characteristic of random variable was expressed according to different variation function. In combination with quadratic orthogonal combination method, the interpolation algorithm of the function of the tunnel lining structure was established, which avoid the problem that the implicit function cannot be solved directly. For the value of tunnel mechanical parameters, earthquake wave method based on advance geology forecast was put forward, which overcome the fuzziness of tunnel mechanical parameter values. At last, compared Monte Carlo method with kriging method, the results showed that the number of iterations was reduced, and the absolute error of the failure probability was only 0.0152%, the relative error was 4.26%. This method was proved to be rational.

Kriging method; tunnel; quadratic orthogonal combination method; advance geological prediction

2015-07-24

國(guó)家自然科學(xué)基金(41130742)

薛洪濤(1975-)，男，遼寧沈陽(yáng)人，工程師，主要從事隧道施工技術(shù)方面的研究，(E-mail)1412887537@qq.com

1673-1549(2016)01-0071-05

10.11863/j.suse.2016.01.15

TU94

A