李丞，蘇立君,3

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

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

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

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

1 不排水抗剪強度參數(shù)非平穩(wěn)隨機場模擬

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

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

(1)

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

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

(2)

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

(3)

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

μcu=μcu+ργd

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

(4)

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

(5)

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

(6)

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

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

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

使用強度折減法，計算邊坡安全系數(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

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

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

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

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

(8)

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

R=C×Pf

(9)

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

(10)

3 算例分析

3.1 土體不排水強度參數(shù)隨機場模擬

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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