韓非,胡超

(1.中航勘察設(shè)計(jì)研究院有限公司,北京 100098;2.四川輕化工大學(xué)土木工程學(xué)院,自貢 643000)

邊坡破壞(滑坡)是常見的自然災(zāi)害之一,頻繁地導(dǎo)致人員、財(cái)產(chǎn)損失以及生態(tài)環(huán)境的破壞[1-3]。邊坡穩(wěn)定性評價(jià)是防災(zāi)減災(zāi)工作的關(guān)鍵,為邊坡設(shè)計(jì)和治理提供了定量依據(jù)。近年來,基于概率論的可靠度分析方法越來越受到學(xué)術(shù)界的關(guān)注[4-6]。該方法將強(qiáng)度參數(shù)視為隨機(jī)變量,綜合考慮了強(qiáng)度參數(shù)的不確定性對邊坡穩(wěn)定性的影響,不僅能給出定量的評價(jià)指標(biāo)(安全系數(shù))還可以給出定性評價(jià)指標(biāo)(可靠度指標(biāo)β和失效概率Pf),為工程設(shè)計(jì)人員提供了更多有益參考。然而,可靠度分析方法的計(jì)算效率一直是推動其在實(shí)際工程中應(yīng)用的重要阻礙[7-8]。

邊坡可靠度分析方法是在概率論的框架下,利用傳統(tǒng)的邊坡穩(wěn)定性分析方法確定失效域范圍,計(jì)算功能函數(shù)在失效域上的概率積分[9-11]。當(dāng)前可靠度分析主要有解析法,近似法[如FORM(first-order reliability method)、SORM(second-order reliability method)和PEM(point estimate method)等],采樣法[如MCS(Monte Carlo simulation)、LHS(Latin hypercube sampling)和SS(subset simulation)等]以及響應(yīng)面法等。采用不同的可靠度分析方法,需要結(jié)合相適應(yīng)的邊坡穩(wěn)定性分析方法,如限平衡法的可靠度分析方法[12-14]和數(shù)值模擬方法[15-17]。極限平衡法需要假定一定數(shù)量的潛在滑動面,然后確定失效概率最大的滑動面,這無疑是個非常艱巨的計(jì)算工作[18-19]。數(shù)值模擬方法能考慮巖土材料的應(yīng)力-應(yīng)變關(guān)系,不需要搜索潛在滑動面[20-22],因此該方法具有非常好的潛力。然而,數(shù)值模擬技術(shù)通過將問題域離散為一定數(shù)量的計(jì)算單元來求解偏微分方程,這就導(dǎo)致其計(jì)算量特別巨大,這種情況對于三維邊坡穩(wěn)定性問題尤為明顯。

邊坡可靠度設(shè)計(jì)通常已知目標(biāo)失效概率,而需要確定擬設(shè)計(jì)的坡角和坡高[23-24]。近年來,基于FORM逆可靠度法在邊坡設(shè)計(jì)中得到了越來越多的關(guān)注。蔣水華等[23]基于一階逆可靠度方法提出了能考慮土體參數(shù)空間變異性的坡角設(shè)計(jì)方法,該方法可以在少量試驗(yàn)數(shù)據(jù)條件下獲得切合工程實(shí)際的邊坡設(shè)計(jì)方案,為坡角設(shè)計(jì)提供了一條新的有效途徑。蘇永華等[25]采用響應(yīng)面方法對邊坡穩(wěn)定性系數(shù)隱式表達(dá)式進(jìn)行近似處理,并與基于FORM逆可靠度算法相結(jié)合,構(gòu)建出一種具有可靠度與穩(wěn)定性系數(shù)雙重控制指標(biāo)的工程邊坡設(shè)計(jì)新方法?；贔ORM逆可靠度方法把坡角和坡高考慮為隨機(jī)變量,通過求解目標(biāo)失效概率或者可靠度指標(biāo),來獲得設(shè)計(jì)坡角和坡高。然而根據(jù)FORM法原理[26-27],計(jì)算得到的坡角和坡高只是概率空間空中到原點(diǎn)距離最近的設(shè)計(jì)點(diǎn),由于考慮了土體參數(shù)的不確定性,在該點(diǎn)的坡角和坡高可能并不是全局最小的坡角和坡高,這就使得設(shè)計(jì)的邊坡可能并不經(jīng)濟(jì)。此外,在逆可靠度方法的求解過程中,坡角和坡高的改變會導(dǎo)致計(jì)算模型幾何形狀發(fā)生變化,但是已有的研究并沒有介紹其實(shí)現(xiàn)的具體細(xì)節(jié)。這主要是由于當(dāng)前采用的邊坡穩(wěn)定性軟件大多是屬于第三方商業(yè)軟件,實(shí)際使用很大程度上受限于軟件提供的接口和功能等。

滑移線場理論避免了對任意滑動面的假設(shè),并在給定的邊界應(yīng)力下產(chǎn)生同時滿足平衡和塑性屈服的應(yīng)力場。Jeldes等[28]將應(yīng)用于邊坡穩(wěn)定性分析中,并提出了臨界坡度的概念,基于該概念可判定邊坡的穩(wěn)定狀態(tài)。Nandi等[29]基于改進(jìn)的擬動力方法,推廣了滑移線場理論,使其能考慮地震對邊坡穩(wěn)定性的影響,有效地考慮了土體的動力特性,如阻尼比和頻率效應(yīng)。然而,基于滑移線場理論的邊坡可靠度分析方法研究鮮有報(bào)道[30]。

綜上所述,為了推進(jìn)基于可靠度的邊坡優(yōu)化設(shè)計(jì)方法的進(jìn)一步發(fā)展,提出了一個基于滑移線場概率邊坡優(yōu)化設(shè)計(jì)方法。該方法擬采用更加高效的擬蒙特卡羅模擬(quasi-Monte Carlo simulation,QMCS)來確定邊坡的失效概率,而在每次模擬過程中采用滑移線場理論來分析邊坡的穩(wěn)定性。為了減少計(jì)算時間,在所提方法中一個簡單而有效地二分搜索方法用來計(jì)算目標(biāo)失效概率的坡角。最后,通過案例研究論證所提方法的有效性。該方法可為基于可靠度邊坡設(shè)計(jì)提供新的技術(shù)手段。

1 基于滑移線理論的邊坡穩(wěn)定性分析

1.1 滑移線理論

基于Mohr-Coulomb強(qiáng)度準(zhǔn)則,邊坡土體中的任意一點(diǎn)應(yīng)力狀態(tài)可表示為

σx=σ(1+sinφcos2θ)-ccotφ

(1)

σy=σ(1-sinφcos2θ)-ccotφ

(2)

τxy=σsinφsin2θ

(3)

式中:σx和σy分別為x和y向的正應(yīng)力;τxy為xy平面內(nèi)的剪應(yīng)力;σ為與應(yīng)力狀態(tài)有關(guān)的特征應(yīng)力;c為黏聚力;φ為內(nèi)摩擦角;θ為最大主應(yīng)力和x軸的夾角。

在平面應(yīng)變條件下,均值各向同性土體的平衡微分方程為

(4)

(5)

式(5)中:γ為土的重度。

將式(1)～式(3)代入式(4)和式(5),可得

(6)

(7)

方程(6)和(7)為雙曲線型方程組,存在兩簇滑移線(α簇和β簇,如圖1所示),可分別表示為

圖1 典型的滑移線場示意圖Fig.1 A typical net of characteristic lines

(8)

式(8)中:2μ為α簇和β簇滑移線之間的夾角,且有μ=π/4-φ/2,其中φ為內(nèi)摩擦角。

(9)

滑移線上每一個點(diǎn)可以由x、y、θ和σ惟一確定。

Mα(xα,yα,θα,σα)為滑移線α簇上的一點(diǎn),Mβ(xβ,yβ,θβ,σβ)為滑移線β簇上的一點(diǎn)。

(10)

(11)

式中:xα和xβ分別為Mα和Mβ點(diǎn)的x向坐標(biāo);yα和yβ分別為Mα和Mβ點(diǎn)的y向坐標(biāo);θα和θβ分別為在Mα和Mβ點(diǎn)處最大主應(yīng)力和x軸的夾角;σα和σβ分別為在Mα和Mβ點(diǎn)處的特征應(yīng)力。

為了計(jì)算邊坡的初始應(yīng)力狀態(tài),在邊坡頂部的點(diǎn)Mα和Mβ可通過式(12)進(jìn)行確定。

(12)

式(12)中:Δx為坡頂?shù)挠?jì)算步長;ω為坡角;H為邊坡高度;P為作用在坡頂?shù)暮奢d。

當(dāng)邊坡頂部沒有外部荷載P時,或者外部荷載P小于式(13)給出的值時,P由式(13)進(jìn)行確定,即

(13)

M(x,y,θ,σ)是滑移線α簇和β簇的交點(diǎn),它可以由Mα和Mβ確定,x、y、θ和σ可分別表示為

(14)

(15)

θ=(σβ-σα)+2(σβθβ+σαθα)tanφ+

γ[(yα-yβ)+(2x-xα-xβ)tanφ]×

[2(σβ+σα)tanφ]-1

(16)

(17)

對于臨界邊坡線的點(diǎn)應(yīng)滿足:

(18)

因此,對于已知的分別位于臨界邊坡線和滑移線β簇上的點(diǎn)Mb(xb,yb,θb,σb)和Mc(xc,yc,θc,σc),依據(jù)差分法,其交點(diǎn)Mij(xij,yij,θij,σij)可表示為

(19)

(20)

式中:xb和xc分別為Mb和Mc點(diǎn)的x向坐標(biāo);yb和yc分別為Mb和Mc點(diǎn)的y向坐標(biāo);θb和θc分別為在Mb和Mc點(diǎn)處最大主應(yīng)力和x軸的夾角;σb和σc分別為在Mb和Mc點(diǎn)處的特征應(yīng)力;xij為Mij點(diǎn)的x向坐標(biāo);yij為Mij點(diǎn)的y向坐標(biāo);θij為在Mij點(diǎn)處最大主應(yīng)力和x軸的夾角;σij為在Mij點(diǎn)處的特征應(yīng)力。

進(jìn)一步結(jié)合式(16)和式(17),則交點(diǎn)Mij(xij,yij,θij,σij)可通過式(21)～式(24)計(jì)算得到。

xij=xbtan(θb-μ)-xctan(θc+μ)-

(yc-yb)[tanθb-tan(θc+μ)]-1

(21)

(22)

θij=(σc-σb)+2(σcθc+σbθb)tanφ+

γ[(yb-yc)+(2xij-xb-xc)tanφ]×

[2(σc+σb)tanφ]-1

(23)

(24)

在邊坡頂點(diǎn)O處,θb和σb可由式(25)、式(26)確定。

(25)

(26)

1.2 破壞準(zhǔn)則及安全系數(shù)

前文已給出了基于滑移線場理論的邊坡穩(wěn)定性分析方法,然而對于如何判定滑坡的極限狀態(tài)仍需要進(jìn)一步闡述。如果已知γ、c、φ、ω和H,基于前述的滑移線場理論可以計(jì)算得到滑坡的滑移線,計(jì)算過程直至臨界邊坡線上的最低點(diǎn)的y*坐標(biāo)值小于0截止(圖1)。由Jeldes等[28]提出的臨界坡度的概念可知,當(dāng)臨界邊坡線通過坡腳的點(diǎn)時,也就是邊坡處于極限平衡狀態(tài),如圖2中曲線b所示;當(dāng)臨界邊坡線處于坡腳右邊時,滑坡處于穩(wěn)定狀態(tài),如圖2中曲線c所示;當(dāng)臨界邊坡線處于坡腳左邊時,滑坡處于不穩(wěn)定狀態(tài),如圖2中曲線a所示。由圖2可知,當(dāng)滑坡處于穩(wěn)定時,計(jì)算過程應(yīng)直至臨界邊坡線上的最低點(diǎn)的y*坐標(biāo)值小于0;而當(dāng)滑坡處于不穩(wěn)定狀態(tài)時,臨界邊坡線處于坡腳左邊,因此計(jì)算過程應(yīng)直至臨界邊坡線上的最低點(diǎn)的x*坐標(biāo)處于等效坡面左側(cè)時,計(jì)算即可停止,從而減小不必要的計(jì)算(圖1)。

圖2 滑移線法確定的臨界邊坡線Fig.2 Determination of the critical slope contour by slip line field theory

基于臨界坡度的概念并結(jié)合傳統(tǒng)的強(qiáng)度折減法,通過不斷折減邊坡強(qiáng)度參數(shù)使得臨界邊坡線通過坡腳,即可獲得邊坡的安全系數(shù)?；趶?qiáng)度折減概念,土體強(qiáng)度參數(shù)可以采用式(27)進(jìn)行折減。

(27)

式(27)中:R為折減因子;cr為折減后的黏聚力;φr為折減后的內(nèi)摩擦角。

當(dāng)邊坡處于極限平衡狀態(tài)時,有安全系數(shù)Fs=R。

2 基于概率的邊坡設(shè)計(jì)

2.1 擬蒙特卡洛模擬

QMCS是一種非常有效的基于抽樣的可靠性分析方法。計(jì)算效率高于蒙特卡羅模擬(MCS),平均收斂速度為O(N-1)[31]。QMCS在概念上類似于MCS,不同之處在于其使用低偏差序列樣本而不是偽隨機(jī)數(shù)。序列樣本是通過低偏差采樣(LDS)生成,并且比MCS中的偽隨機(jī)數(shù)生成的序列樣本更均勻[32]。Halton序列是常用的LDS采樣方法[33],采用Halton數(shù)列來估計(jì)邊坡的失效概率。

擬蒙特卡洛模擬的計(jì)算過程可分為3個步驟:首先,依據(jù)xi的概率分布生成N個隨機(jī)樣本;然后,將每個隨機(jī)樣本輸入到功能函數(shù)中計(jì)算對應(yīng)的值。最后,統(tǒng)計(jì)失效樣本的數(shù)量,并采用式(28)來確定邊坡的失效概率Pf。

(28)

式(28)中:xi為第i個樣本;I(·)為指標(biāo)函數(shù),如果g(xi)≤0,則I[g(xi)]=1,否則I[g(xi)]=0。采用滑移線場理論分析邊坡穩(wěn)定性時,如果g(xi)≤0時,即對應(yīng)的是圖2中曲線a和曲線b情況;如果g(x)>0時,即對應(yīng)的是圖2中曲線c情況。由式(28)可知,失效概率Pf的精度與樣本數(shù)量直接相關(guān)。為了保證計(jì)算精度以及計(jì)算效率,采用Pf的變異系數(shù)COVPf來估計(jì)樣本數(shù)量。

(29)

2.2 設(shè)計(jì)邊坡概率估計(jì)

圖3 不同坡角、坡高邊坡的安全系數(shù)Fig.3 Safety factor of slopes under different slope angles and slopes

(30)

圖4 計(jì)算流程圖Fig.4 Flowchart of the proposed method

為了實(shí)現(xiàn)所提邊坡設(shè)計(jì)方法,基于Python平臺編制了相應(yīng)的計(jì)算程序。計(jì)算程序包括3個腳本文件Main.py、Slip_line_Slope.py、Cal_Pf.py以及Cal_Fs.py。Main.py腳本用于設(shè)置初始參數(shù),展示最終計(jì)算結(jié)果;Slip_line_Slope.py腳本用于分析邊坡穩(wěn)定性;Cal_Pf.py腳本用于計(jì)算邊坡失效概率;Cal_Fs.py腳本用于計(jì)算邊坡安全系數(shù)。

3 算例研究

表1 土體剪切強(qiáng)度參數(shù)Table 1 Shear strength parameters of soils

圖5 設(shè)計(jì)邊坡橫截面圖Fig.5 Cross section of design slope

3.1 邊坡參數(shù)統(tǒng)計(jì)特征

表2 剪切強(qiáng)度參數(shù)分布類型驗(yàn)證(α=0.05)Table 2 Verification of probabilistic distribution of shear strength parameters(α=0.05)

圖6 強(qiáng)度參數(shù)頻數(shù)直方圖Fig.6 Histogram of shear strength parameters

相關(guān)系數(shù)ρ是隨機(jī)場的重要參數(shù),描述兩個隨機(jī)變量的聯(lián)合分布結(jié)構(gòu)。對于服從Lognormal分布的黏聚力和內(nèi)摩擦角,在不同相關(guān)系數(shù)下,基于蒙特卡洛模擬生成的樣本如圖7所示。圖7表明盡管黏聚力和內(nèi)摩擦角都具有相同邊緣分布,但聯(lián)合分布是不同的。因此,有必要研究隨機(jī)場的相關(guān)性對隨機(jī)分析結(jié)果的影響。采用Spearman法確定了黏聚力和內(nèi)摩擦角的相關(guān)系數(shù)ρc,φ=-0.13。

圖7 不同互相關(guān)系數(shù)下c和φ的聯(lián)合概率分布Fig.7 Joint probability distribution of c and φ under different cross-correlation coefficients

為了驗(yàn)證基于滑移線場理論邊坡穩(wěn)定性分析的有效性,以及自編腳本的可靠性。在邊坡未開挖時采用土體參數(shù)的平均值,計(jì)算得到的邊坡安全系數(shù)為1.27,該計(jì)算結(jié)果與GEOSTUDIO軟件采用Morgenstern-Price法和Spencer法計(jì)算得到的1.288結(jié)果相近。為了從分驗(yàn)證基于滑移線場理論邊坡穩(wěn)定性分析的有效性,計(jì)算了在不同坡角下的邊坡安全系數(shù),計(jì)算結(jié)果如表3所示?？梢钥闯?本文計(jì)算結(jié)果與GEOSTUDIO軟件計(jì)算結(jié)果幾乎一致。

表3 不同方法計(jì)算的安全系數(shù)對比Table 3 Comparison of safety factors calculated by different methods

在邊坡未開挖時采用擬蒙特卡洛模擬,計(jì)算得到的邊坡安全系數(shù)的平均值為1.27,失效概率為3.56×10-2。盡管安全系數(shù)的平均值表明邊坡在為開挖狀態(tài)下邊坡處于穩(wěn)定狀態(tài),但是此時邊坡的坡角較大,邊坡失效概率并不滿足目標(biāo)失效概率要求,因此需進(jìn)行邊坡挖開處理。

3.2 概率邊坡優(yōu)化設(shè)計(jì)

圖8 計(jì)算過程Fig.8 Calculation process

為了驗(yàn)證所提方法的有效性,在現(xiàn)有的土體參數(shù)以及計(jì)算得到的坡角值的基礎(chǔ)上,利用直接蒙特卡洛模擬(MCS),重新計(jì)算的邊坡失效概率分別為1×10-2、1×10-3和1×10-4,計(jì)算結(jié)果與目標(biāo)失效概率一致,驗(yàn)證了所提方法的有效性。

4 參數(shù)研究

由圖9(a)、圖9(b)可知,隨著黏聚力和內(nèi)摩擦角平均值的增加設(shè)計(jì)坡角將逐步增大,而隨著黏聚力和內(nèi)摩擦角標(biāo)準(zhǔn)差的增加設(shè)計(jì)坡角逐步減小。隨著黏聚力和內(nèi)摩擦角平均值的增加,土體抗剪強(qiáng)度逐步增大,邊坡穩(wěn)定性逐步增加,因此設(shè)計(jì)坡角也隨之逐漸增大。當(dāng)黏聚力和內(nèi)摩擦角的標(biāo)準(zhǔn)差增加時,盡管黏聚力和摩擦角的平均值未發(fā)生變化,但是邊坡失效概率也隨之增大,因此要滿足目標(biāo)失效概率的設(shè)計(jì)坡角也應(yīng)減小。由圖9(c)可知,隨著重度平均值和標(biāo)準(zhǔn)差的增加設(shè)計(jì)坡角逐步增小。由圖9(d)可知,隨著相關(guān)系數(shù)的減小,設(shè)計(jì)坡角逐步增大。隨著相關(guān)系數(shù)的減小,黏聚力和內(nèi)摩擦角逐步呈負(fù)相關(guān)性,由圖7(c)可知,此時在失效域范圍的隨機(jī)樣本數(shù)量越少,邊坡失效概率較小,因此在相對較大的設(shè)計(jì)坡角下即可滿足目標(biāo)失效概率。

5 結(jié)論

將滑移線場理論應(yīng)用于邊坡設(shè)計(jì),提出了基于失效概率的邊坡設(shè)計(jì)新方法。新方法采用了二分搜索算法來求解滿足目標(biāo)失效概率解,并且在每次搜索過程中使用了能更加快速收斂的擬蒙特卡洛模擬方法來計(jì)算邊坡失效概率。以一個實(shí)際邊坡為例進(jìn)行了邊坡概率優(yōu)化設(shè)計(jì)。得出如下主要結(jié)論。

(1)所提方法計(jì)算效率高,最多只需36.29 min就能得到滿足目標(biāo)失效概率的結(jié)果。傳統(tǒng)的基于極限平衡法和數(shù)值模擬方法的邊坡穩(wěn)定性分析,在搜索過程中邊坡的坡角和坡高會不斷變化,從而導(dǎo)致邊坡幾何模型的變化,這使得計(jì)算過程的模型重構(gòu)工作十分復(fù)雜。所提出的邊坡概率優(yōu)化設(shè)計(jì)方法可簡化這一問題的難度,提高了計(jì)算效率,為一般的隱式邊坡可靠度設(shè)計(jì)奠定基礎(chǔ)。

(2)參數(shù)研究表明,土體參數(shù)的平均值、標(biāo)準(zhǔn)差以及相關(guān)系數(shù)與設(shè)計(jì)坡角的規(guī)律,研究成果為基于概率論的邊坡優(yōu)化設(shè)計(jì)提供了新的方法。

(3)僅將土體參數(shù)考慮為隨機(jī)變量,忽略了土體參數(shù)空間變異性,具有一定的局限性。如何將所提方法進(jìn)一步推廣能更加全面的考慮土體參數(shù)的空間異性將是需要進(jìn)一步研究的內(nèi)容。