韩非,胡超

(1.中航勘察设计研究院有限公司,北京 100098;2.四川轻化工大学土木工程学院,自贡 643000)

边坡破坏(滑坡)是常见的自然灾害之一,频繁地导致人员、财产损失以及生态环境的破坏[1-3]。边坡稳定性评价是防灾减灾工作的关键,为边坡设计和治理提供了定量依据。近年来,基于概率论的可靠度分析方法越来越受到学术界的关注[4-6]。该方法将强度参数视为随机变量,综合考虑了强度参数的不确定性对边坡稳定性的影响,不仅能给出定量的评价指标(安全系数)还可以给出定性评价指标(可靠度指标β和失效概率Pf),为工程设计人员提供了更多有益参考。然而,可靠度分析方法的计算效率一直是推动其在实际工程中应用的重要阻碍[7-8]。

边坡可靠度分析方法是在概率论的框架下,利用传统的边坡稳定性分析方法确定失效域范围,计算功能函数在失效域上的概率积分[9-11]。当前可靠度分析主要有解析法,近似法[如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)等]以及响应面法等。采用不同的可靠度分析方法,需要结合相适应的边坡稳定性分析方法,如限平衡法的可靠度分析方法[12-14]和数值模拟方法[15-17]。极限平衡法需要假定一定数量的潜在滑动面,然后确定失效概率最大的滑动面,这无疑是个非常艰巨的计算工作[18-19]。数值模拟方法能考虑岩土材料的应力-应变关系,不需要搜索潜在滑动面[20-22],因此该方法具有非常好的潜力。然而,数值模拟技术通过将问题域离散为一定数量的计算单元来求解偏微分方程,这就导致其计算量特别巨大,这种情况对于三维边坡稳定性问题尤为明显。

边坡可靠度设计通常已知目标失效概率,而需要确定拟设计的坡角和坡高[23-24]。近年来,基于FORM逆可靠度法在边坡设计中得到了越来越多的关注。蒋水华等[23]基于一阶逆可靠度方法提出了能考虑土体参数空间变异性的坡角设计方法,该方法可以在少量试验数据条件下获得切合工程实际的边坡设计方案,为坡角设计提供了一条新的有效途径。苏永华等[25]采用响应面方法对边坡稳定性系数隐式表达式进行近似处理,并与基于FORM逆可靠度算法相结合,构建出一种具有可靠度与稳定性系数双重控制指标的工程边坡设计新方法。基于FORM逆可靠度方法把坡角和坡高考虑为随机变量,通过求解目标失效概率或者可靠度指标,来获得设计坡角和坡高。然而根据FORM法原理[26-27],计算得到的坡角和坡高只是概率空间空中到原点距离最近的设计点,由于考虑了土体参数的不确定性,在该点的坡角和坡高可能并不是全局最小的坡角和坡高,这就使得设计的边坡可能并不经济。此外,在逆可靠度方法的求解过程中,坡角和坡高的改变会导致计算模型几何形状发生变化,但是已有的研究并没有介绍其实现的具体细节。这主要是由于当前采用的边坡稳定性软件大多是属于第三方商业软件,实际使用很大程度上受限于软件提供的接口和功能等。

滑移线场理论避免了对任意滑动面的假设,并在给定的边界应力下产生同时满足平衡和塑性屈服的应力场。Jeldes等[28]将应用于边坡稳定性分析中,并提出了临界坡度的概念,基于该概念可判定边坡的稳定状态。Nandi等[29]基于改进的拟动力方法,推广了滑移线场理论,使其能考虑地震对边坡稳定性的影响,有效地考虑了土体的动力特性,如阻尼比和频率效应。然而,基于滑移线场理论的边坡可靠度分析方法研究鲜有报道[30]。

综上所述,为了推进基于可靠度的边坡优化设计方法的进一步发展,提出了一个基于滑移线场概率边坡优化设计方法。该方法拟采用更加高效的拟蒙特卡罗模拟(quasi-Monte Carlo simulation,QMCS)来确定边坡的失效概率,而在每次模拟过程中采用滑移线场理论来分析边坡的稳定性。为了减少计算时间,在所提方法中一个简单而有效地二分搜索方法用来计算目标失效概率的坡角。最后,通过案例研究论证所提方法的有效性。该方法可为基于可靠度边坡设计提供新的技术手段。

1 基于滑移线理论的边坡稳定性分析

1.1 滑移线理论

基于Mohr-Coulomb强度准则,边坡土体中的任意一点应力状态可表示为

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

(1)

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

(2)

τxy=σsinφsin2θ

(3)

式中:σx和σy分别为x和y向的正应力;τxy为xy平面内的剪应力;σ为与应力状态有关的特征应力;c为黏聚力;φ为内摩擦角;θ为最大主应力和x轴的夹角。

在平面应变条件下,均值各向同性土体的平衡微分方程为

(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,其中φ为内摩擦角。

(9)

滑移线上每一个点可以由x、y、θ和σ惟一确定。

Mα(xα,yα,θα,σα)为滑移线α簇上的一点,Mβ(xβ,yβ,θβ,σβ)为滑移线β簇上的一点。

(10)

(11)

式中:xα和xβ分别为Mα和Mβ点的x向坐标;yα和yβ分别为Mα和Mβ点的y向坐标;θα和θβ分别为在Mα和Mβ点处最大主应力和x轴的夹角;σα和σβ分别为在Mα和Mβ点处的特征应力。

为了计算边坡的初始应力状态,在边坡顶部的点Mα和Mβ可通过式(12)进行确定。

(12)

式(12)中:Δx为坡顶的计算步长;ω为坡角;H为边坡高度;P为作用在坡顶的荷载。

当边坡顶部没有外部荷载P时,或者外部荷载P小于式(13)给出的值时,P由式(13)进行确定,即

(13)

M(x,y,θ,σ)是滑移线α簇和β簇的交点,它可以由Mα和Mβ确定,x、y、θ和σ可分别表示为

(14)

(15)

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

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

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

(16)

(17)

对于临界边坡线的点应满足:

(18)

因此,对于已知的分别位于临界边坡线和滑移线β簇上的点Mb(xb,yb,θb,σb)和Mc(xc,yc,θc,σc),依据差分法,其交点Mij(xij,yij,θij,σij)可表示为

(19)

(20)

式中:xb和xc分别为Mb和Mc点的x向坐标;yb和yc分别为Mb和Mc点的y向坐标;θb和θc分别为在Mb和Mc点处最大主应力和x轴的夹角;σb和σc分别为在Mb和Mc点处的特征应力;xij为Mij点的x向坐标;yij为Mij点的y向坐标;θij为在Mij点处最大主应力和x轴的夹角;σij为在Mij点处的特征应力。

进一步结合式(16)和式(17),则交点Mij(xij,yij,θij,σij)可通过式(21)～式(24)计算得到。

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)

在边坡顶点O处,θb和σb可由式(25)、式(26)确定。

(25)

(26)

1.2 破坏准则及安全系数

前文已给出了基于滑移线场理论的边坡稳定性分析方法,然而对于如何判定滑坡的极限状态仍需要进一步阐述。如果已知γ、c、φ、ω和H,基于前述的滑移线场理论可以计算得到滑坡的滑移线,计算过程直至临界边坡线上的最低点的y*坐标值小于0截止(图1)。由Jeldes等[28]提出的临界坡度的概念可知,当临界边坡线通过坡脚的点时,也就是边坡处于极限平衡状态,如图2中曲线b所示;当临界边坡线处于坡脚右边时,滑坡处于稳定状态,如图2中曲线c所示;当临界边坡线处于坡脚左边时,滑坡处于不稳定状态,如图2中曲线a所示。由图2可知,当滑坡处于稳定时,计算过程应直至临界边坡线上的最低点的y*坐标值小于0;而当滑坡处于不稳定状态时,临界边坡线处于坡脚左边,因此计算过程应直至临界边坡线上的最低点的x*坐标处于等效坡面左侧时,计算即可停止,从而减小不必要的计算(图1)。

图2 滑移线法确定的临界边坡线Fig.2 Determination of the critical slope contour by slip line field theory

基于临界坡度的概念并结合传统的强度折减法,通过不断折减边坡强度参数使得临界边坡线通过坡脚,即可获得边坡的安全系数。基于强度折减概念,土体强度参数可以采用式(27)进行折减。

(27)

式(27)中:R为折减因子;cr为折减后的黏聚力;φr为折减后的内摩擦角。

当边坡处于极限平衡状态时,有安全系数Fs=R。

2 基于概率的边坡设计

2.1 拟蒙特卡洛模拟

QMCS是一种非常有效的基于抽样的可靠性分析方法。计算效率高于蒙特卡罗模拟(MCS),平均收敛速度为O(N-1)[31]。QMCS在概念上类似于MCS,不同之处在于其使用低偏差序列样本而不是伪随机数。序列样本是通过低偏差采样(LDS)生成,并且比MCS中的伪随机数生成的序列样本更均匀[32]。Halton序列是常用的LDS采样方法[33],采用Halton数列来估计边坡的失效概率。

拟蒙特卡洛模拟的计算过程可分为3个步骤:首先,依据xi的概率分布生成N个随机样本;然后,将每个随机样本输入到功能函数中计算对应的值。最后,统计失效样本的数量,并采用式(28)来确定边坡的失效概率Pf。

(28)

式(28)中:xi为第i个样本;I(·)为指标函数,如果g(xi)≤0,则I[g(xi)]=1,否则I[g(xi)]=0。采用滑移线场理论分析边坡稳定性时,如果g(xi)≤0时,即对应的是图2中曲线a和曲线b情况;如果g(x)>0时,即对应的是图2中曲线c情况。由式(28)可知,失效概率Pf的精度与样本数量直接相关。为了保证计算精度以及计算效率,采用Pf的变异系数COVPf来估计样本数量。

(29)

2.2 设计边坡概率估计

图3 不同坡角、坡高边坡的安全系数Fig.3 Safety factor of slopes under different slope angles and slopes

(30)

图4 计算流程图Fig.4 Flowchart of the proposed method

为了实现所提边坡设计方法,基于Python平台编制了相应的计算程序。计算程序包括3个脚本文件Main.py、Slip_line_Slope.py、Cal_Pf.py以及Cal_Fs.py。Main.py脚本用于设置初始参数,展示最终计算结果;Slip_line_Slope.py脚本用于分析边坡稳定性;Cal_Pf.py脚本用于计算边坡失效概率;Cal_Fs.py脚本用于计算边坡安全系数。

3 算例研究

表1 土体剪切强度参数Table 1 Shear strength parameters of soils

图5 设计边坡横截面图Fig.5 Cross section of design slope

3.1 边坡参数统计特征

表2 剪切强度参数分布类型验证(α=0.05)Table 2 Verification of probabilistic distribution of shear strength parameters(α=0.05)

图6 强度参数频数直方图Fig.6 Histogram of shear strength parameters

相关系数ρ是随机场的重要参数,描述两个随机变量的联合分布结构。对于服从Lognormal分布的黏聚力和内摩擦角,在不同相关系数下,基于蒙特卡洛模拟生成的样本如图7所示。图7表明尽管黏聚力和内摩擦角都具有相同边缘分布,但联合分布是不同的。因此,有必要研究随机场的相关性对随机分析结果的影响。采用Spearman法确定了黏聚力和内摩擦角的相关系数ρc,φ=-0.13。

图7 不同互相关系数下c和φ的联合概率分布Fig.7 Joint probability distribution of c and φ under different cross-correlation coefficients

为了验证基于滑移线场理论边坡稳定性分析的有效性,以及自编脚本的可靠性。在边坡未开挖时采用土体参数的平均值,计算得到的边坡安全系数为1.27,该计算结果与GEOSTUDIO软件采用Morgenstern-Price法和Spencer法计算得到的1.288结果相近。为了从分验证基于滑移线场理论边坡稳定性分析的有效性,计算了在不同坡角下的边坡安全系数,计算结果如表3所示。可以看出,本文计算结果与GEOSTUDIO软件计算结果几乎一致。

表3 不同方法计算的安全系数对比Table 3 Comparison of safety factors calculated by different methods

在边坡未开挖时采用拟蒙特卡洛模拟,计算得到的边坡安全系数的平均值为1.27,失效概率为3.56×10-2。尽管安全系数的平均值表明边坡在为开挖状态下边坡处于稳定状态,但是此时边坡的坡角较大,边坡失效概率并不满足目标失效概率要求,因此需进行边坡挖开处理。

3.2 概率边坡优化设计

图8 计算过程Fig.8 Calculation process

为了验证所提方法的有效性,在现有的土体参数以及计算得到的坡角值的基础上,利用直接蒙特卡洛模拟(MCS),重新计算的边坡失效概率分别为1×10-2、1×10-3和1×10-4,计算结果与目标失效概率一致,验证了所提方法的有效性。

4 参数研究

由图9(a)、图9(b)可知,随着黏聚力和内摩擦角平均值的增加设计坡角将逐步增大,而随着黏聚力和内摩擦角标准差的增加设计坡角逐步减小。随着黏聚力和内摩擦角平均值的增加,土体抗剪强度逐步增大,边坡稳定性逐步增加,因此设计坡角也随之逐渐增大。当黏聚力和内摩擦角的标准差增加时,尽管黏聚力和摩擦角的平均值未发生变化,但是边坡失效概率也随之增大,因此要满足目标失效概率的设计坡角也应减小。由图9(c)可知,随着重度平均值和标准差的增加设计坡角逐步增小。由图9(d)可知,随着相关系数的减小,设计坡角逐步增大。随着相关系数的减小,黏聚力和内摩擦角逐步呈负相关性,由图7(c)可知,此时在失效域范围的随机样本数量越少,边坡失效概率较小,因此在相对较大的设计坡角下即可满足目标失效概率。

5 结论

将滑移线场理论应用于边坡设计,提出了基于失效概率的边坡设计新方法。新方法采用了二分搜索算法来求解满足目标失效概率解,并且在每次搜索过程中使用了能更加快速收敛的拟蒙特卡洛模拟方法来计算边坡失效概率。以一个实际边坡为例进行了边坡概率优化设计。得出如下主要结论。

(1)所提方法计算效率高,最多只需36.29 min就能得到满足目标失效概率的结果。传统的基于极限平衡法和数值模拟方法的边坡稳定性分析,在搜索过程中边坡的坡角和坡高会不断变化,从而导致边坡几何模型的变化,这使得计算过程的模型重构工作十分复杂。所提出的边坡概率优化设计方法可简化这一问题的难度,提高了计算效率,为一般的隐式边坡可靠度设计奠定基础。

(2)参数研究表明,土体参数的平均值、标准差以及相关系数与设计坡角的规律,研究成果为基于概率论的边坡优化设计提供了新的方法。

(3)仅将土体参数考虑为随机变量,忽略了土体参数空间变异性,具有一定的局限性。如何将所提方法进一步推广能更加全面的考虑土体参数的空间异性将是需要进一步研究的内容。