不同隨機(jī)場(chǎng)模擬對(duì)斜坡失穩(wěn)風(fēng)險(xiǎn)評(píng)估結(jié)果的影響

李丞，蘇立君,3

(1.中國(guó)科學(xué)院 山地災(zāi)害與地表過程重點(diǎn)實(shí)驗(yàn)室；中國(guó)科學(xué)院、水利部 成都山地災(zāi)害與環(huán)境研究所，成都 610041； 2.中國(guó)科學(xué)院大學(xué)，北京100049；3.中國(guó)科學(xué)院 青藏高原地球科學(xué)卓越創(chuàng)新中心，北京100101)

已有研究表明[1-3]，土體參數(shù)空間變異性對(duì)斜坡穩(wěn)定性影響顯著。通常，斜坡滑動(dòng)面傾向于尋找最不利的失穩(wěn)路徑，這些滑動(dòng)面不僅存在一定的空間變異性，而且相互之間具有相關(guān)性，因此，斜坡存在多種滑動(dòng)面形式[1-2]。從工程角度來看，不僅要考慮斜坡的穩(wěn)定性，還應(yīng)考慮斜坡的失穩(wěn)后果，若不考慮滑動(dòng)面間相關(guān)性的影響，就會(huì)高估斜坡的失穩(wěn)風(fēng)險(xiǎn)[3]。近年來，斜坡失穩(wěn)風(fēng)險(xiǎn)評(píng)價(jià)受到了較多的關(guān)注。Li等[4]將子集模擬算法引入隨機(jī)有限元算法，為小概率水平下的斜坡失穩(wěn)風(fēng)險(xiǎn)計(jì)算提供了一種有效方法。Li等[5]提出一種基于具有代表性滑動(dòng)面的極限平衡法，可為邊坡失穩(wěn)風(fēng)險(xiǎn)提供定量評(píng)價(jià)。Xiao等[6]提出輔助隨機(jī)有限元法，使得隨機(jī)有限元算法在三維斜坡風(fēng)險(xiǎn)評(píng)估中得到應(yīng)用。Jiang等[7]在極限平衡分析框架下，提出一種考慮土體二維空間變異性的邊坡失穩(wěn)風(fēng)險(xiǎn)定量評(píng)估方法。

上述研究中，均采用平穩(wěn)隨機(jī)場(chǎng)模型來描述土體性質(zhì)的空間變異性[8]，然而，土體受沉積條件、地質(zhì)和環(huán)境等作用影響，導(dǎo)致相應(yīng)的土體參數(shù)呈現(xiàn)沿埋深變化的趨勢(shì)[9]。值得注意的是，目前，已有學(xué)者[10-13]認(rèn)識(shí)到土體參數(shù)非平穩(wěn)分布特征對(duì)斜坡概率分析的影響。由文獻(xiàn)[10-13]可知，文獻(xiàn)[4-7]高估了斜坡失穩(wěn)概率，因此，在評(píng)價(jià)斜坡失穩(wěn)風(fēng)險(xiǎn)時(shí)需考慮土體抗剪強(qiáng)度參數(shù)的非平穩(wěn)分布特征。

學(xué)者們很少關(guān)注使用非平穩(wěn)隨機(jī)場(chǎng)表征土體抗剪強(qiáng)度參數(shù)來評(píng)價(jià)斜坡失穩(wěn)風(fēng)險(xiǎn)。針對(duì)這一問題，筆者采用隨機(jī)場(chǎng)模擬土體抗剪強(qiáng)度參數(shù)分布特征，結(jié)合有限元極限分析法、強(qiáng)度折減法和Monte Carlo模擬來計(jì)算斜坡安全系數(shù)和失效概率，并討論了使用非平穩(wěn)隨機(jī)場(chǎng)和平穩(wěn)隨機(jī)場(chǎng)模擬土體抗剪強(qiáng)度參數(shù)分布特征的差異對(duì)斜坡失穩(wěn)風(fēng)險(xiǎn)的影響。

1 不排水抗剪強(qiáng)度參數(shù)非平穩(wěn)隨機(jī)場(chǎng)模擬

土體參數(shù)的空間變異性由趨勢(shì)參數(shù)和隨機(jī)波動(dòng)參數(shù)組成，不同深度處的土體參數(shù)ξ(d)可表示為[9]

ξ(d)=t(d)+w(d)

(1)

式中：d為地面以下土體深度；t(d)為趨勢(shì)參數(shù)函數(shù)，即土體參數(shù)在埋深d處的均值；w(d)為隨機(jī)波動(dòng)參數(shù)函數(shù)，w(d)為均值和標(biāo)準(zhǔn)差不隨深度變化的統(tǒng)計(jì)平穩(wěn)隨機(jī)場(chǎng)[9]。t(d)與土體物質(zhì)組成、沉積條件和固結(jié)過程等有關(guān)[14]。對(duì)于正常固結(jié)粘土層，t(d)從地表處沿深度增加(起始值為0)；對(duì)于高度超固結(jié)粘土層，t(d)沿深度不發(fā)生變化；對(duì)于土層較厚的超固結(jié)粘土層，t(d)從地表處沿深度增加(起始值為某固定值)。盡管土體參數(shù)沿深度方向可能存在非線性變化的趨勢(shì)，為簡(jiǎn)單起見，t(d)一般選擇簡(jiǎn)單的線性函數(shù)[15]。為此，只研究超固結(jié)粘土層趨勢(shì)參數(shù)函數(shù)從地表處沿深度增加的情況。據(jù)此，土體不排水抗剪強(qiáng)度參數(shù)cu由某一定值cu0隨深度增加的表達(dá)式為[11,13]

cu(d)=cu0+ρσv=cu0+ργd

(2)

式中：cu0為地面處土體的不排水抗剪強(qiáng)度，kPa；ρ為土體不排水抗剪強(qiáng)度隨深度增加的速率，kPa/m；σv為垂直有效應(yīng)力，kPa，σv=γd；γ為土體重度，kN/m3，d為地面以下土體深度，m。采用文獻(xiàn)[11,13]的方法，模擬土體不排水抗剪強(qiáng)度的非平穩(wěn)分布特征：首先，將cu0模擬為均值為μcu0和標(biāo)準(zhǔn)差為σcu0的對(duì)數(shù)正態(tài)平穩(wěn)隨機(jī)場(chǎng)，然后，考慮cu隨深度線性變化的影響，就此得到cu的二維非平穩(wěn)隨機(jī)場(chǎng)為

(3)

根據(jù)式(3)，得cu的均值和標(biāo)準(zhǔn)差分別為

μcu=μcu+ργd

σcu(d)=COVcu0(μcu0+ργd)

(4)

模擬隨機(jī)場(chǎng)的數(shù)據(jù)為土體參數(shù)的均值、方差以及方差折減函數(shù)，而方差折減函數(shù)取決于相關(guān)函數(shù)及相關(guān)距離。方差折減函數(shù)用于描述空間平均后的方差和點(diǎn)方差的關(guān)系，如式(5)所示。

(5)

式中：D為平均間距；ρ(θ)為相關(guān)函數(shù)，可由協(xié)方差函數(shù)求得。Li等[16]研究了不同相關(guān)函數(shù)對(duì)邊坡可靠度的影響，發(fā)現(xiàn)相關(guān)函數(shù)對(duì)可靠度影響不大，因此，選用形式簡(jiǎn)單的指數(shù)型相關(guān)函數(shù)

(6)

式中：x1、x2、y1、y2是隨機(jī)場(chǎng)域內(nèi)的坐標(biāo)；lh和lv分別表示水平和豎直相關(guān)距離。通過Karhunen-Loeve級(jí)數(shù)展開法[17]生成隨機(jī)場(chǎng)。

2 斜坡安全系數(shù)計(jì)算和失穩(wěn)風(fēng)險(xiǎn)評(píng)估方法

采用有限元極限分析下限法進(jìn)行斜坡穩(wěn)定性計(jì)算。有限元極限分析下限法指出：在任何靜力許可應(yīng)力場(chǎng)內(nèi)可以計(jì)算真實(shí)極限荷載的下限值(或“安全值”)。靜力許可應(yīng)力場(chǎng)需滿足平衡條件、應(yīng)力邊界條件和屈服條件(應(yīng)力必須位于應(yīng)力空間的屈服面內(nèi)部或之上)。下限解是求解滿足靜力許可條件的最大破壞荷載[18]。

使用強(qiáng)度折減法，計(jì)算邊坡安全系數(shù)[18]

(7)

圖1給出了有限元極限分析下限法的斜坡安全系數(shù)分析流程[18]。TOL為收斂精度，文中取0.001。

圖1 有限元極限分析下限解的斜坡安全系數(shù)分析流程Fig.1 Analysis flow of lower limit solution of finite element limit analysis of slope safety factor

在強(qiáng)度折減過程中，采用K聚類方法(K-means clustering method)[17]搜索滑面的形狀和位置，據(jù)此滑面計(jì)算滑坡體積。

斜坡每一次計(jì)算都會(huì)生成一個(gè)安全系數(shù)。對(duì)于一系列隨機(jī)場(chǎng)實(shí)現(xiàn)，使用Monte Carlo方法可獲得全部安全系數(shù)，Monte Carlo方法的具體思路如圖2所示。因此，Monte Carlo方法被用于產(chǎn)生足夠數(shù)量的不排水抗剪強(qiáng)度隨機(jī)場(chǎng)，并進(jìn)行斜坡穩(wěn)定性分析。

圖2 利用Monte Carlo模擬計(jì)算斜坡失效概率和失效風(fēng)險(xiǎn) 的流程圖Fig.2 The flow chart of calculating slope failure probability and failure risk by Monte Carlo simulation

用式(8)計(jì)算斜坡失效概率。

(8)

在得到斜坡失效概率及失效后果后，斜坡失穩(wěn)風(fēng)險(xiǎn)R可以被評(píng)價(jià)為

R=C×Pf

(9)

式中：C為斜坡失穩(wěn)后的量化結(jié)果[17]。式(9)僅適用于斜坡滑動(dòng)面形式唯一時(shí)的情況。然而，在評(píng)價(jià)空間變異性的斜坡中可能存在多種滑動(dòng)面形式，所以，每種斜坡滑動(dòng)面形式相關(guān)的后果都應(yīng)該單獨(dú)評(píng)估，于是，采用文獻(xiàn)[17,19]中的公式評(píng)估斜坡失穩(wěn)風(fēng)險(xiǎn)。

(10)

3 算例分析

3.1 土體不排水強(qiáng)度參數(shù)隨機(jī)場(chǎng)模擬

以文獻(xiàn)[13]的一個(gè)不排水飽和粘土斜坡為例進(jìn)行分析，如圖3所示。根據(jù)文獻(xiàn)[13]可知，斜坡高度10 m，坡比為1∶2，土體重度γ=20 kN/m3。

圖3 不排水飽和粘土斜坡Fig.3 The undrained saturated clay slope

斜坡模型劃分為1 000個(gè)三角形單元，Monte Carlo模擬次數(shù)在計(jì)算結(jié)果滿足精度的條件下，計(jì)算次數(shù)為1 000次。為描述cu沿深度方向的非平穩(wěn)分布特征，即cu均值和標(biāo)準(zhǔn)差隨埋深的增加而增大的特性，將cu模擬為μcu0=14.669 kPa和σcu0=4.034 kPa的對(duì)數(shù)正態(tài)平穩(wěn)隨機(jī)場(chǎng)，水平相關(guān)距離δh=19 m，豎直相關(guān)距離δv=1.9 m，參數(shù)ρ為0.15，變異系數(shù)為25%。同時(shí)，為了比較采用非平穩(wěn)隨機(jī)場(chǎng)和平穩(wěn)隨機(jī)場(chǎng)模擬土體抗剪強(qiáng)度參數(shù)空間變異性的差異對(duì)斜坡可靠度和失穩(wěn)風(fēng)險(xiǎn)的影響，參照文獻(xiàn)[10,20]的做法，取斜坡中部(z=-10 m)處的不排水強(qiáng)度均值和變異系數(shù)分別為44.67 kPa和27.5%，模擬對(duì)應(yīng)的cu對(duì)數(shù)正態(tài)平穩(wěn)隨機(jī)場(chǎng)實(shí)現(xiàn)值，也取δh=19 m和δv=1.9 m。

為了便于說明，定義使用非平穩(wěn)隨機(jī)場(chǎng)模擬土體抗剪強(qiáng)度參數(shù)時(shí)為工況1，使用平穩(wěn)隨機(jī)場(chǎng)模擬土體抗剪強(qiáng)度參數(shù)時(shí)為工況2。圖4和圖5給出了cu(x,d)沿深度方向(圖3中的AB截面)的3次典型實(shí)現(xiàn)值。由圖4可知，工況1的cu均值隨著埋深的增加而增大，通過分析cu典型實(shí)現(xiàn)值沿埋深方向的變化趨勢(shì)可知，σcu同樣沿埋深增加，說明此模型能夠模擬cu的非平穩(wěn)分布特征。由圖5可知，工況2由于沒有考慮趨勢(shì)分量的影響[11]，cu值的大小與深度沒有關(guān)系，其相應(yīng)的cu典型實(shí)現(xiàn)值呈現(xiàn)雜亂無章的變化，導(dǎo)致土體內(nèi)部隨機(jī)場(chǎng)實(shí)現(xiàn)值變化較大。此外，由工況1獲得的隨機(jī)場(chǎng)實(shí)現(xiàn)值，大多分布在cu均值的右側(cè)，由工況2獲得的隨機(jī)場(chǎng)實(shí)現(xiàn)值在cu均值兩側(cè)較均勻的波動(dòng)。

圖4 工況1的不同cu實(shí)現(xiàn)值Fig.4 Different cu for condition 1

圖5 工況2的不同cu實(shí)現(xiàn)值Fig.5 Different cu for condition 2

將工況1和工況2獲得的二維非平穩(wěn)隨機(jī)場(chǎng)和二維平穩(wěn)隨機(jī)場(chǎng)cu(x,d)典型實(shí)現(xiàn)值分別賦予給斜坡模型，如圖6和圖7所示。由圖6和圖7可知，紅色部分為cu值較大區(qū)域，藍(lán)色部分為cu值較小區(qū)域。此外，土性參數(shù)空間變異性導(dǎo)致多種滑動(dòng)面形式，滑動(dòng)面位置和形狀及其安全系數(shù)均發(fā)生明顯變化，即：斜坡安全系數(shù)與滑動(dòng)面位置和形狀沒有直接聯(lián)系。

圖6 工況1某次典型實(shí)現(xiàn)值及滑動(dòng)面形式Fig.6 A typical realization value and failure mode for condition 1

圖7 工況2某次典型實(shí)現(xiàn)值及滑動(dòng)面形式Fig.7 A typical realization value and failure mode for condition 2

3.2 斜坡可靠度分析和失穩(wěn)風(fēng)險(xiǎn)定量評(píng)估

獲得cu隨機(jī)場(chǎng)實(shí)現(xiàn)值后進(jìn)行斜坡可靠度分析與失穩(wěn)風(fēng)險(xiǎn)定量評(píng)估。采用Monte Carlo進(jìn)行非平穩(wěn)隨機(jī)場(chǎng)和平穩(wěn)隨機(jī)場(chǎng)斜坡可靠度計(jì)算得到的Pf分別為0.044和0.062，與文獻(xiàn)[13]計(jì)算的非平穩(wěn)隨機(jī)場(chǎng)Pf=0.044 1基本一致，而且采用非平穩(wěn)隨機(jī)場(chǎng)計(jì)算得到的Pf小于采用平穩(wěn)隨機(jī)場(chǎng)計(jì)算得到的Pf，這與文獻(xiàn)[20-21]得出的結(jié)論吻合。工況1和工況2情況下混合滑動(dòng)體積的直方圖如圖8和圖9所示。從圖8和圖9可知，當(dāng)考慮土體抗剪強(qiáng)度參數(shù)非平穩(wěn)分布特征時(shí)，斜坡滑動(dòng)體積的峰值約在250 m3左右，而不考慮土體抗剪強(qiáng)度參數(shù)非平穩(wěn)分布特征時(shí)，斜坡滑動(dòng)體積的峰值約在750 m3左右。

圖8 工況1混合滑動(dòng)體積直方圖Fig.8 Mixed sliding volume histogram for condition 1

淺層滑動(dòng)面位置和深層滑動(dòng)面位置的定義如文獻(xiàn)[17]所述，本文不再描述。不同工況下斜坡發(fā)生不同滑動(dòng)面位置的直方圖如圖10和圖11所示，從圖10和圖11可知，深層滑動(dòng)面是主要的滑動(dòng)面形式，在工況2情況下，深層滑動(dòng)面的比例遠(yuǎn)大于工況1。工況1和工況2的斜坡失穩(wěn)風(fēng)險(xiǎn)分別為9.8、35.69 m3。

圖9 工況2混合滑動(dòng)體積直方圖Fig.9 Mixed sliding volume histogram condition 2

圖10 工況1分解滑動(dòng)體積直方圖Fig.10 Decomposition sliding volume histogram for condition 1

3.3 空間變異性對(duì)失效概率和失穩(wěn)風(fēng)險(xiǎn)的影響

圖11 工況2分解滑動(dòng)體積直方圖Fig.11 Decomposition sliding volume histogram for condition 2

圖12 豎直相關(guān)距離對(duì)失效概率和失穩(wěn)風(fēng)險(xiǎn)的影響Fig.12 Effectofvertical correlation distance on failure probability and failure risk

圖13 豎直相關(guān)距離對(duì)滑動(dòng)面形式的影響Fig.13 Effectofvertical correlation distance on sliding surface types

4 結(jié)論

結(jié)合有限元極限分析、強(qiáng)度折減法和蒙特卡洛模擬，比較了使用非平穩(wěn)隨機(jī)場(chǎng)模型和平穩(wěn)隨機(jī)場(chǎng)模型對(duì)計(jì)算斜坡失穩(wěn)風(fēng)險(xiǎn)和滑動(dòng)面位置和形狀的影響，得到以下結(jié)論：

1)當(dāng)土體參數(shù)存在空間變異性時(shí)，僅由斜坡安全系數(shù)是無法準(zhǔn)確評(píng)價(jià)滑坡體積大小和位置及所產(chǎn)生的失穩(wěn)風(fēng)險(xiǎn)，故需與具體滑動(dòng)面形式及滑坡體積結(jié)合，進(jìn)行滑坡風(fēng)險(xiǎn)綜合分析。

2)采用平穩(wěn)隨機(jī)場(chǎng)計(jì)算得到的斜坡深層滑動(dòng)面的數(shù)量遠(yuǎn)遠(yuǎn)大于采用非平穩(wěn)隨機(jī)場(chǎng)。

3)斜坡豎向空間變異性對(duì)采用非平穩(wěn)隨機(jī)場(chǎng)計(jì)算其失穩(wěn)風(fēng)險(xiǎn)的影響較小，對(duì)失效概率的影響較大，相比之下，斜坡豎向空間變異性對(duì)采用平穩(wěn)隨機(jī)場(chǎng)計(jì)算其失穩(wěn)風(fēng)險(xiǎn)和失效概率的影響顯著。

4)在邊坡穩(wěn)定性理論分析中(如極限分析法和極限平衡法)，需預(yù)先設(shè)定滑動(dòng)面位置和形狀。在這種情況下，精確的滑動(dòng)面位置和形狀可以大大提高理論方法的分析精度。不同隨機(jī)場(chǎng)模擬得到的邊坡滑動(dòng)面位置和形狀及失穩(wěn)風(fēng)險(xiǎn)可以幫助工程師預(yù)測(cè)和解釋實(shí)際場(chǎng)地的邊坡滑動(dòng)面位置和形狀，從而提高滑坡危險(xiǎn)性評(píng)價(jià)的準(zhǔn)確性。