土木建筑与环境工程  2013, Vol. 35 Issue (6): 33-39   PDF    
伪蒙特卡罗法及其在边坡可靠度分析中的应用
褚雪松, 李亮    
青岛理工大学 土木工程学院, 山东 青岛 266033
收稿日期:2013-04-18
基金项目:国家自然科学基金(51274126、51008167);山东省高等学校科技计划项目(J10LE07);教育部高等学校博士学科点专项科研基金(20103721120001);大连理工大学海岸与近海工程国家重点实验室开放基金项目(LP1214)
作者简介:褚雪松(1977-), 女, 博士, 主要从事岩土工程防灾减灾研究, (E-mail)celldl@126.com
摘要:在确定具有最小可靠度指标的滑动面(即临界可靠度滑动面)时,由于常规的蒙特卡罗法抽样耗时巨大,临界可靠度滑动面的获得较为耗时。对于均质边坡,利用简化Bishop法构建了可靠度分析的功能函数,设计了6种随机变量的标准差组合,假定了随机变量的4种抽样范围,利用抽样次数较小的蒙特卡罗法即伪蒙特卡罗法对随机生成的132组可行滑动面进行了伪可靠度指标的计算并与蒙特卡罗法计算得到的可靠度指标进行了比较分析,研究发现:只有一个随机变量的前提下,滑动面的伪可靠度指标与蒙特卡罗法计算的可靠度指标呈完全线性关系,在其它条件相同的情况下,伪蒙特卡罗法抽样范围越大,伪可靠度指标与蒙特卡罗法计算的可靠度指标之间的拟合直线斜率越小,反之亦然;伪蒙特卡罗法抽样次数越大,伪可靠度指标与蒙特卡罗法计算的可靠度指标之间的拟合直线斜率越大,反之亦然。对均质边坡,可应用伪蒙特卡罗法快速计算其临界可靠度指标。
关键词边坡稳定    可靠度分析    极限平衡法    蒙特卡罗法    
Analysis on the Quasi Monte-Carlo Method and its Application in the Slope Reliability
Chu Xuesong, Li Liang    
School of Civil Engineering, Qingdao Technological University, Qingdao 266033, Shandong, P. R. China
Received: 2013-04-18
Abstract: The critical reliability slip surface with minimum value of reliability index is difficult to be obtained by the Monte-Carlo method. The simplified Bishop method was adopted to calculate the factor of safety for given circular slip surface thereby defining the performance function in reliability analysis. Six variance values of random variables were designed and four sample ranges were assumed. A quasi-Monte Carlo method was defined with fow samples. One hundred and thirty-two potential slip surfaces were randomly generated. The reliability indexes regarding those slip surfaces were calculated by Monte-Carlo method and quasi-Monte Carlo method, respectively. The comparative results show that in the case of one random variable, the reliability index from Monte Carlo method is found to be linear with that from quasi-Monte Carlo method. In the similar conditions, the wider the sample range, the smaller the tangent value of line, and Vice versa. In addition, the bigger the sample number, the bigger the tangent value of line. The Quasi-Monte Carlo method is approved to be practicable for the determination of critical reliability slip surface for homogeneous slope.
Key Words: slope stability    reliability analysis    limit equilibrium method    Monte-Carlo method    

目前边坡稳定性分析中,虽然基于极限平衡方法(诸如简化Bishop法[1]、SARMA法[2]、摩根斯坦普莱斯法[3]以及不平衡推力法[4])的确定性分析仍然在工程中广泛应用,但边坡稳定可靠度分析也越来越受到岩土工程界的关注。与确定性分析仅能得到关于滑动面的一个安全系数不同,可靠度分析能得到关于滑动面的安全系数均值及标准差,所以能更合理地评价边坡的稳定性。边坡可靠度分析与传统的极限平衡方法是密切相关的,当前可靠度分析中所需要的功能函数一般是基于极限平衡方法构建的,对于给定的某个滑动面,可采用蒙特卡罗法[5-10]、一次二阶矩法[11-13]、响应面法[14]等来计算其可靠度指标(或者失效概率),几种方法各有优缺点,譬如蒙特卡罗法可避免求偏导数等复杂的数学运算,但是抽样次数一般较大,与滑动面的失效概率有关,失效概率越小,抽样次数越大[5];一次二阶矩法仅需要随机变量的均值和方差即可进行可靠度指标的求取,然而却需要多次与多层的迭代;响应面法是一种近似模拟功能函数的方法,或者说是一种便于功能函数求导运算的途径。在边坡可靠度分析中,最终需要确定具有最小可靠度指标的滑动面(称之为临界可靠度滑动面,下同),临界可靠度滑动面的确定步骤与具有最小安全系数滑动面(称之为临界滑动面,下同)的确定类似,适用于临界滑动面确定的许多搜索方法(譬如传统的网格法、穷举法以及近期发展的智能方法等)均可以用于临界可靠度滑动面的搜索中[5]。两者的区别在于,临界滑动面的获取仅需要几分钟甚至更小的时间,然而临界可靠度滑动面的确定却非常耗时,其原因是显而易见的,譬如蒙特卡罗法的上千、万次的抽样(每一次抽样意味着一次安全系数的计算),以及一次二阶矩法的多层迭代等。为了快速地对边坡进行可靠度分析,往往先得到临界滑动面,然后再对临界滑动面进行一次可靠度指标的计算作为评价之用,这就导致临界可靠度滑动面与临界滑动面是一致的,但显然这是不符合工程实际的,本文针对均质边坡探索一种快速确定临界可靠度滑动面的方法。

1 临界可靠度滑动面的确定

临界可靠度滑动面确定的第1步是如何从数学上描述潜在滑动面,作为一般情况,可假定潜在滑动面为任意形状滑动面,笔者曾对数种可行的构造方法进行了比较分析[15]。本文针对均质边坡进行,因此可假定滑动面为圆弧,需要3个变量确定一条滑动面;

第2步就是采用一种方法计算可行滑动面的可靠度指标,譬如蒙特卡罗法、一次二阶矩法等。本文假定随机变量符合正态分布,利用蒙特卡罗法计算其可靠度指标,功能函数为g=Fs-1,Fs可以是任意一种极限平衡方法得到的安全系数,本文用简化Bishop法[1]计算功能函数中的安全系数;

第3步就是变化潜在滑动面寻求临界可靠度滑动面的策略,即搜索算法。关于搜索算法,目前算法种类繁多,本文选用简单易行的和声搜索算法[16]进行,关于和声搜索算法在边坡稳定分析中的具体应用可参见文献[17]。

2 伪蒙特卡罗法抽样
2.1 蒙特卡罗法抽样

蒙特卡罗法适用于随机变量的概率密度分布形式已知或符合假定的情况,它避免了功能函数对随机变量的求导运算从而可方便地与任意一种极限平衡方法结合。实际上,蒙特卡罗法直接求解的是滑动面的失效概率,其主要思想为:首先对随机变量进行大量的随机抽样,然后将这些抽样值逐个代入功能函数并判断功能函数是否小于零(即失效),最终将失效的抽样次数与总的抽样次数之比定义为失效概率,当然也可以根据功能函数的均值和标准差来计算与此失效概率相对应的可靠度指标。

蒙特卡罗法的最大缺陷就是抽样次数巨大,一般可用下式(1)近似确定最小的抽样次数[5]

$ T \ge \frac{{100}}{{{p_{\rm{f}}}}} $ (1)

式中:T为抽样次数,pf为滑动面的失效概率。

对于随机生成的可行滑动面而言,其失效概率一般较小,因此T值非常大。实际上,如果考虑最终结果,即找到临界可靠度滑动面而言,能不能采用次数较少的抽样来计算滑动面的“伪”可靠度指标呢?因为抽样次数减少后所得到的可靠度指标肯定不是滑动面真正的可靠度指标,所以此处称之为伪可靠度指标,相应地抽样称之为伪蒙特卡罗法。该伪可靠度指标虽然不准确,但是与其真正的可靠度指标之间如果存在对应关系的话,就可以用于临界可靠度滑动面的搜素中。

2.2 伪蒙特卡罗法抽样

对于伪蒙特卡罗法抽样来说,第一步就是要确定随机变量的变化范围,研究表明:对于一个具有正态分布的随机变量,99.73%的数据介于[μ-3σ, μ+3σ]之间,此处μ为随机变量的均值,σ为其标准差。为便于本文计算之用,给定了4种伪蒙特卡罗法抽样的范围,分别是[μσ, μ+σ](范围1,余者类推),[μ-1.5σ, μ+1.5σ],[μ-2σ, μ+2σ],[μ-3σ, μ+3σ]等。设抽样次数为Q,可采用下式(2)来确定抽样的值。

$ {x_i} = {L_x} + \left( {i - 1} \right)\frac{{{U_x} - {L_x}}}{{Q - 1}}, \mathit{i = }{\rm{1, 2, }} \cdots \mathit{Q} $ (2)

式中,xi是随机变量的第i次抽样值,Lx是抽样范围的下限,Ux是上限,如上介绍的4种范围,如果采用范围1的话,Lx=μσUx=μ+σ,余者类推。

第2步就是确定Q,为便于临界可靠度滑动面搜索,Q当然越小越好,但是Q不能等于1,因为抽样一次,意味着仅得到一个安全系数,安全系数和功能函数的标准差均为零,无法计算伪可靠度指标,因此Q≥2。

2.3 滑动面可靠度指标计算

前面介绍的蒙特卡罗法以及伪蒙特卡罗法抽样得到的是TQ个随机变量的值,对于给定的滑动面,将这些随机变量值代入简化Bishop法[1]公式,即可得到TQ个安全系数,以伪蒙特卡罗法为例,其安全系数系列为F1, F2, …, FQ,可用式(3)来计算其伪可靠度指标。

$ {\beta _q} = \frac{{{\mu _{\rm{g}}}}}{{{\sigma _{\rm{g}}}}} $ (3)

式中${\mu _{\rm{g}}} = \frac{{\sum\limits_{i = 1}^Q {{g_i}} }}{Q}, {\sigma _{\rm{g}}} = \sqrt {\frac{{\sum\limits_{i = 1}^Q {{{({g_i} -{\mu _{\rm{g}}})}^2}} }}{{Q -1}}}。$

类似地,如果滑动面的可靠度指标是由蒙特卡罗法计算的话,同样利用式(3)计算,只不过公式中的Q要用T来代替,不再细叙。在下文的算例分析中,利用设定的4种范围,假定Q=2、3来比较随机生成的一系列可行滑动面的真正可靠度指标与其相应的伪可靠度指标之间是否存在对应的关系。

3 算例分析

考虑一均质土坡,其剖面如图 1所示,不考虑土体强度参数中容重γ的随机性,只考虑粘聚力c与内摩擦角$φ$之间的一个随机变量(表 1μcσc分别是粘聚力的均值与标准差,类似地,${\mu _\varphi }$${\sigma _\varphi }$分别是内摩擦角的均值与标准差)。为了比较参数变异程度(标准差的大小)对计算结果的影响,设计了如表 1所示的6种情况,表 1中第2行称之为情况1,依此类推。对每种情况下均可得到其临界可靠度滑动面,相应地称之为临界可靠度滑动面1、2、…、6。

图 1 均质土坡剖面

表 1 土层强度参数及标准差

3.1 伪蒙特卡罗法抽样试验

设计的数值试验如下:

1) 首先随机生成150组滑动面,经过可行性验证之后,剩下132组可行的滑动面;然后利用蒙特卡罗法计算132组可行滑动面的可靠度指标,抽样次数应该与每一个滑动面的失效概率相关,为方便计,取T=200 000次;最后利用伪蒙特卡罗法计算132组可行滑动面的伪可靠度指标,采用上文介绍的4种范围,抽样次数Q=2、3,以滑动面的可靠度指标作为横坐标,以滑动面的伪可靠度指标为纵坐标,绘制散点图如下所示。

2) 鉴于Q=1时无法进行伪蒙特卡罗法抽样,对132组可行滑动面利用简化Bishop法计算其安全系数,同样以滑动面的可靠度指标为横坐标,以滑动面的安全系数为纵坐标,绘制散点图如下所示。

注:图 35791113中,图例范围1、2、3、4即2.2节定义的4个范围,并且给出了纵坐标y对横坐标x拟合直线,R为相关系数,下同。

图 3 情况1下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=2)

图 5 情况2下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=2)

图 7 情况3下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=2)

图 9 情况4下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=2)

图 11 情况5下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=2)

图 13 情况6下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=2)

3.2 伪蒙特卡罗法抽样试验分析

图 24681012所示的滑动面可靠度指标与滑动面的安全系数散点图,由图可见,两者之间的关系较为杂乱,从垂直线上看,许多滑动面具有同一个可靠度指标,然而其安全系数却是不同的,也就是说不存在一一对应的关系,如此一来,根据安全系数为导向搜索得到的临界滑动面肯定与临界可靠度滑动面是不同的。从水平线上来看,同样存在多对一的关系。

图 2 情况1下滑动面的可靠度指标与安全系数之间的关系

图 4 情况2下滑动面的可靠度指标与安全系数之间的关系

图 6 情况3下滑动面的可靠度指标与安全系数之间的关系

图 8 情况4下滑动面的可靠度指标与安全系数之间的关系

图 10 情况5下滑动面的可靠度指标与安全系数之间的关系

图 12 情况6下滑动面的可靠度指标与安全系数之间的关系

图 35791113以及14~19所示的滑动面可靠度指标与滑动面的伪可靠度指标的散点图,由图可见:

图 14 情况1下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=3)

图 15 情况2下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=3)

图 16 情况3下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=3)

图 17 情况4下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=3)

图 18 情况5下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=3)

图 19 情况6下滑动面的可靠度指标与伪可靠度指标之间的关系(Q=3)

1) 假定的4种范围,对表 1设计的每种情况而言,散点图均呈现出完全的线性关系,这可以由相关系数R=1看出。

2)Q一定的话,对表 1设计的每种情况而言,相同的范围得到的拟合直线斜率基本一致,譬如当Q=2时,范围1为0.7左右,范围2为0.47,范围3为0.353,范围4为0.235。

3) 对每种情况,其它条件一样的话,范围越大,拟合的直线斜率越小,反之亦然。

4) 其它条件一样的话,Q越大,拟合得到的直线斜率就越大,反之亦然。

由以上4条结论可知,对于表 1设计的6种情况,可以应用任何一种范围来计算滑动面的伪可靠度指标,然后再换算成可靠度指标用于临界可靠度滑动面搜索就变得快速、可行。

3.3 伪可靠度指标的应用

图 1所示的土坡,在表 1所设计的6种情况下,临界滑动面是一致的,因为确定性分析用到的参数均值是一样的,临界滑动面示于图 20中,在表 1设计的6种情况下,对临界滑动面进行伪可靠度指标的计算,然后利用图 35791113所得到的拟合直线公式换算得到其这正的可靠度指标列于表 2中。此外,对表 1设计的6种情况直接应用伪蒙特卡罗法进行临界可靠度指标的搜索(范围1,Q=2),得到的临界滑动面示于图 20中,相应的可靠度指标示于表 2中。由图 20可见,临界滑动面位于临界可靠度滑动面4、5、6与临界可靠度滑动面1、2、3之间,临界可靠度滑动面1、2、3基本重合,临界可靠度滑动面4、5、6也基本重合。再从表 2可知,对临界滑动面进行蒙特卡罗法抽样得到的可靠度指标均比临界可靠度滑动面对应的可靠度指标偏大,最大误差达25%左右。

图 20 均质土坡临界滑动面与临界可靠度滑动面比较

表 2 临界滑动面及临界可靠度滑动面可靠度指标比较

4 结论

鉴于蒙特卡罗法抽样次数巨大导致的临界可靠度滑动面搜索困难的问题,针对一均质边坡,在只有一个随机变量的前提下,提出了一种伪蒙特卡罗抽样方法,比较了伪蒙特卡罗法抽样次数以及抽样范围对结果的影响,初步研究得到以下结论:

1) 对均质边坡而言,滑动面的伪可靠度指标与蒙特卡罗法计算的可靠度指标呈完全线性关系,伪蒙特卡罗法抽样范围以及抽样次数Q都没有影响。

2) 在其它条件相同的情况下,伪蒙特卡罗法抽样范围越大,伪可靠度指标与蒙特卡罗法计算的可靠度指标之间的拟合直线斜率越小,反之亦然。

3) 在其它条件相同的情况下,伪蒙特卡罗法抽样次数越大,伪可靠度指标与蒙特卡罗法计算的可靠度指标之间的拟合直线斜率越大,反之亦然。

4) 对均质边坡,可应用Q=2的伪蒙特卡罗法进行临界可靠度指标的搜索,然后换算为蒙特卡罗法的可靠度指标。

参考文献
[1] Bishop A W. The use of the slip circle in the stability analysis of slopes[J]. Geotechnique, 1955, 5(1): 7–17. DOI:10.1680/geot.1955.5.1.7
[2] Sarma S K. Stability analusis of embankments and slopes[J]. Geotechnique, 1973, 23(3): 423–433. DOI:10.1680/geot.1973.23.3.423
[3] Morgenstern N R, Price V E. The analysis of the stability of general slip surfaces[J]. Geotechnique, 1965, 15(1): 79–93. DOI:10.1680/geot.1965.15.1.79
[4] 时卫民, 郑颖人, 唐伯明, 等. 土坡稳定不平衡推力法的精度分析及其使用条件[J]. 岩土工程学报, 2004, 26(3): 313–317.
Shi W M, Zheng Y R, Tang B M, et al. Accuracy and application range of imbalance thrust force method for slope stability analysis[J]. Chinese Journal of Geotechnical Engineering, 2004, 26(3): 313–317. (in Chinese)
[5] 陈祖煜. 土质边坡稳定分析-原理·方法·程序[M]. 北京: 中国水利水电出版社, 2003: 372-373.
[6] 吴振君, 王水林, 汤华, 等. 边坡可靠度分析的一种新的优化求解方法[J]. 岩土力学, 2010, 31(3): 713–718.
Wu Z J, Wang S L, Tang H, et al. A new optimization approach for slope reliability analysis[J]. Rock and Soil Mechanics, 2010, 31(3): 713–718. (in Chinese)
[7] 吴振君, 王水林, 葛修润. LHS方法在边坡可靠度分析中的应用[J]. 岩土力学, 2010, 31(4): 1047–1054.
Wu Z J, Wang S L, Ge X R. Application of latin hypercube sampling technique to slope reliability analysis[J]. Rock and Soil Mechanics, 2010, 31(4): 1047–1054. (in Chinese)
[8] Au S K, Cao Z J, Wang Y. Implementing advanced Monte Carlo simulation under spreadsheet environment[J]. Structural Safety, 2010, 32(5): 281–292. DOI:10.1016/j.strusafe.2010.03.004
[9] Wang Y, Cao Z, Au S K. Practical reliability analysis of slope stability by advanced Monte Carlo simulations in a spreadsheet[J]. Canadian Geotechnical Journal, 2011, 48(1): 162–172. DOI:10.1139/T10-044
[10] Zhang J, Zhang L M, Tang W H. New methods for system reliability analysis of soil slopes[J]. Canadian Geotechnical Journal, 2011, 48: 1138–1148. DOI:10.1139/t11-009
[11] el Ramly H, Morgenstern N R, Cruden D M. Probabilistic slope stability analysis for practice[J]. Canadian Geotechnical Journal, 2002, 39(3): 665–683. DOI:10.1139/t02-034
[12] Duncan J M. Factors of safety and reliability in geotechnical engineering[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2000, 126(4): 307–316. DOI:10.1061/(ASCE)1090-0241(2000)126:4(307)
[13] Li K S, Lumb P. Probabilistic design of slopes[J]. Canadian Geotechnical Journal, 1987, 24(4): 520–535. DOI:10.1139/t87-068
[14] 吕大刚, 贾明明, 李刚. 结构可靠度分析的均匀设计响应面法[J]. 工程力学, 2011, 28(7): 109–116.
Lyu D G, Jia M M, Li G. Uniform design response surface method for structural reliability analysis[J]. Engineering Mechanics, 2011, 28(7): 109–116. DOI:10.6052/j.issn.1000-4750.2009.10.0748 (in Chinese)
[15] 李亮, 王玉杰, 王秋生, 等. 土坡稳定分析中模拟任意滑动面的新策略及其效率分析[J]. 水利学报, 2008, 39(5): 535–541.
Li L, Wang Y J, Wang Q S, et al. New procedure for simulating arbitrary slip surface of soil slope in stability analysis[J]. Journal of Hydraulic Engineering, 2008, 39(5): 535–541. (in Chinese)
[16] Geem Z W, Kim J H, Loganathan G V. Harmony search[J]. Simulation, 2001, 76(2): 60–68. DOI:10.1177/003754970107600201
[17] 李亮, 于广明, 褚雪松. 边坡临界滑动场方法的优化算法实现[J]. 岩土工程学报, 2010, 31(6): 827–831.
Li L, Yu G M, Chu X S. The Simulation of critical slip field method of slopes based on optimization algorithm[J]. Chinese Journal of Geotechnical Engineering, 2010, 31(6): 827–831. (in Chinese)
    伪蒙特卡罗法及其在边坡可靠度分析中的应用
    褚雪松, 李亮