边坡可靠度分析的随机响应面法及程序实现

第 29 卷 第 8 期 岩石力学与工程学报 Vol.29 No.8 2010 年 8 月 Chinese Journal of Rock Mechanics and Engineering Aug.，2010 收稿日期：2010–01–31；修回日期：2010–04–27 基金项目：国家自然科学基金重点项目(50839004)；国家杰出青年科学基金项目(50725931)；“十一五”国家科技支撑 计划 项目进度计划表范例计划下载计划下载计划下载课程教学计划下载 重大项目(2008BAB29B01)； 教育部新世纪优秀人才计划项目(NCET–08–0415) 作者简介：李典庆(1975–)，男，博士，1998 年毕业于河海大学机电工程学院机械设计与制造专业，现任教授、博士生导师，主要从事岩土工程 可靠度及风险调控理论、高坝大库系统寿命评估理论及风险分析等方面的教学与研究工作。E-mail：dianqing@whu.edu.cn 边坡可靠度分析的随机响应面法及程序实现 李典庆 1，2，周创兵 1，2，陈益峰 1，2，姜清辉 1，2，荣 冠 1，2 (1. 武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉 430072；2. 武汉大学 水工岩石力学教育部重点实验室，湖北 武汉 430072) 摘要：提出分析相关非正态变量可靠度计算问题的随机响应面法，采用 Nataf 变换成功地解决输入变量相关时随 机响应面法的配点问题及可靠度计算问题。推导 4～6 阶 Hermite 随机多项式展开的解析 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 达式，并编写基于 C# 语言的随机响应面法计算程序。以岩质边坡平面滑动破坏模式为例证明随机响应面法在边坡可靠度分析中的有效 性。研究结果表明，基于 Nataf 变换的随机响应面法能够有效分析含有相关非正态变量的边坡可靠度问题。随机 响应面法的计算精度优于传统的 FORM 方法，其计算效率高于传统的蒙特卡罗模拟方法，其收敛性在数学意义上 是有保证的。随机多项式展开的阶数几乎对边坡安全系数均值的估计没有影响，但是在边坡失效概率的计算中要 选择适当的随机多项式展开的阶数。在基于随机响应面法的可靠度分析框架内，边坡安全系数计算和可靠度分析 2 个过程分开独立进行，同时计算安全系数和失效概率能够更加系统地进行边坡稳定性分析。研究成果为拓展随 机响应面法在边坡可靠度分析中的应用奠定了一定的基础。 关键词：边坡工程；边坡可靠度；随机响应面法；随机多项式；概率配点；相关非正态变量 中图分类号：P 642.22 文献标识码：A 文章编号：1000–6915(2010)08–1513–11 RELIABILITY ANALYSIS OF SLOPE USING STOCHASTIC RESPONSE SURFACE METHOD AND CODE IMPLEMENTATION LI Dianqing1，2，ZHOU Chuangbing1，2，CHEN Yifeng1，2，JIANG Qinghui1，2，RONG Guan1，2 (1. State Key Laboratory of Water Resources and Hydropower Engineering Science，Wuhan University，Wuhan，Hubei 430072， China；2. Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering，Ministry of Education，Wuhan University， Wuhan，Hubei 430072，China) Abstract：This paper aims at proposing a stochastic response surface method(SRSM) for reliability analysis involving correlated random variables. The Nataf transformation is adopted to effectively transform the correlated nonnormal random variables into the independent standard normal variables，which facilitates the collocation points associated with the correlated random variables and reliability computation using the SRSM. Explicit polynomials are derived for fourth-order to sixth-order Hermite polynomial chaos expansions of any number of random variables. A C#-language based computer program WHUSRSM (Wuhan University SRSM) is developed. An example of reliability analysis of rock slope stability with plane failure is presented to demonstrate the validity and capability of the proposed SRSM. The results indicate that the proposed SRSM can evaluate the reliability of rock slope stability involving correlated random variables efficiently. The proposed SRSM has a higher accuracy than the conventional first-order reliability method；and its efficiency is higher than Monte Carlo simulations. • 1514 • 岩石力学与工程学报 2010年 Moreover，the convergence in Hilbert space of the SRSM can be ensured in any case. The mean of factor of safety for slope stability can be accurately estimated by the proposed SRSM with different orders. However，the orders of the SRSM should be selected carefully to accurately estimate the probability of rock slope failure. The calculation of factor of safety and the reliability analysis can be conducted separately within the framework of the SRSM based on reliability analysis. Slope stability analysis would be investigated systematically with the results associated with factor of safety as well as probability of failure. These results can provide a basis for extending the application of the SRSM to reliability analysis of rock slope stability. Key words：slope engineering；slope reliability；stochastic response surface method；polynomial chaos；probabilistic collocation；correlated nonnormal random variables 1 引 言 随着可靠度理论的迅速发展和日趋完善，它在 岩土工程中日益得到重视并已经得到广泛的应用[1]。 如 2010 年 1 月 15～16 日在宜昌三峡工地召开的全 国性水利水电工程风险分析及可靠度设计方法研讨 会又涌现出一批有关水利水电工程结构可靠度分析 的 论文 政研论文下载论文大学下载论文大学下载关于长拳的论文浙大论文封面下载 。在国内，陈祖煜等[2]认为，岩土工程中应 该采用考虑不确定性因素的可靠度分析方法。在国 外，以 H. T. Wilson 等为代表的学者也在积极推进 岩土工程可靠度分析方法的研究[3]。规范[4]明确建 议，对于一级边坡有条件时应进行可靠度分析，可 见边坡采用可靠度分析方法是大势所趋。 目前常用的边坡可靠度计算方法有一次二阶矩 法(FORM)[5]、响应面法(RSM)[6]、蒙特卡罗模拟方 法(MCS)[7]、自适应重要抽样方法[8]、J. M. Duncan[9] 提出的简易分析方法等。本文目的在于提出一种新 的边坡可靠度计算方法——随机响应面法(stochastic response surface method，SRSM)。S. S. Isukapalli 等[10，11]提出了随机响应面法，将其应用于环境和生 物系统的不确定性分析，并开发了 Web 版本的随机 响应面法计算程序 WebSRSM。S. P. Huang 等[12]采 用 KL 正交级数展开将随机场进行离散，研究了基 于概率配点的随机响应面法在结构可靠度分析中的 应用。S. P. Huang 等[13]编写了基于 EXCEL 的随机 响应面法计算程序 CSRSM，并探讨了随机响应面法 在边坡可靠度分析中的应用。王 卿和黄淑萍[14]采用 基于概率配点的随机响应面法分析了单桩沉降可靠 度。为了考虑随机变量间的相关性，S. S. Isukapalli[10] 将基于乔列斯基分解的相关变量独立化方法引入到 随机响应面法中，但并没有给出算例来证明所提方 法的有效性。经过分析发现，Isukapalli 方法只适用 于相关正态变量的情况，这也是许多可靠度计算方 法称可以解决相关变量问题最容易隐藏的缺点之 一。 尽管随机响应面法在可靠度分析中得到了一定 的应用，也取得了一些可喜的进展，然而综合国内 外研究可以看出，目前基于随机响应面法的可靠度 研究主要存在以下 3 个方面不足：首先，现有随机 响应面法均不能分析含有相关非正态变量的可靠度 问题。众所周知，岩体力学参数如黏聚力和内摩擦 角间存在明显的统计负相关性，岩质边坡张裂缝中 充水深度百分比和张裂缝深度也存在明显的负相关 关系，而且张裂缝中充水深度百分比服从截尾指数 分布[5]。因此，如果不解决相关非正态变量的可靠 度分析问题，随机响应面法在边坡可靠度分析中的 应用将会大受限制。其次，目前随机响应面法的 Hermite 随机多项式展开的解析表达式仅限于 2 和 3 阶，当然它们对于一些简单的情况如输入随机变量 近似服从正态分布、失效概率大于 10－3时是有效的。 但是当输入随机变量或输出随机响应量是强非正态 分布变量时如均匀分布变量，常用的 2 和 3 阶随机 多项式的计算精度远远不够。另外，当所分析问题 的失效概率较小，如失效概率量级为 10－3～10－4时， 2 和 3 阶随机多项式的误差也将很大。因此，有必 要采用更高阶的 Hermite 随机多项式展开，遗憾的 是目前还未有 3 阶以上 Hermite 随机多项式展开的 解析表达式，这严重妨碍了高阶 Hermite 随机多项 式展开在可靠度计算中的应用。最后，现有随机响 应面法考虑的变量分布类型也非常有限，如常用的 截尾指数分布与 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 正态随机变量的映射关系未见 报道。针对上述 3 个问题，本文提出了考虑相关非 正态变量的随机响应面法，采用 Nataf 变换方法[15] 第 29 卷 第 8 期 李典庆，等. 边坡可靠度分析的随机响应面法及程序实现 • 1515 • 分析随机响应面法中输入非正态变量间的相关性问 题。推导了 4～6 阶 Hermite 随机多项式展开的解析 表达式。最后采用算例证明了所提随机响应面法在 边坡可靠度分析中的有效性，并定量地研究了边坡 稳定性分析的安全系数和可靠度之间的关系。 2 随机响应面法基本原理 随机响应面法最早是由 Isukapalli 于 1998 年在 研究环境和生物系统的不确定性问题中提出来的。 从数学理论背景来讲，随机响应面法最早可追溯 到 N. Wiener[16]对布朗运动的研究，他在对布朗运 动的随机过程理论研究中提出了齐次随机函数的 概念。随机响应面法分析的第一步是将输入随机变 量向量 X 表示为标准正态分布随机变量 U 的函数， 即 1[ ( )]F X U (1) 式中： 1[ ]F   为 X 累计概率分布函数的反函数， ( )  为标准正态分布的累计概率分布函数，U 为标准正 态分布随机变量。如对于正态分布的随机变量 X 来 说，它和标准正态分布随机变量之间的映射关系[11] 如下： X X  X U (2) 式中： X ， X 分别为 X 的均值和标准差。输入随 机变量的个数(随机空间的自由度)根据定义问题的 需要确定，但为了减少计算工作量，通常选择最敏 感的随机变量作为输入随机变量。 在此基础上，把输出响应也表示为相同的标准 正态随机变量的函数。输出响应量 Y 通常采用如下 的 Hermite 随机多项式[17]展开来表示： 1 1 1 1 2 1 2 1 1 2 1 2 1 2 3 1 2 3 1 2 3 0 0 1 2 1 1 1 3 1 1 1 ( ) ( ) ( ) i i i i i i i i i i i i i i i i i i i i i a a U a U U a U U U                       Y  ， ， ， 11 2 1 2 1 2 1 2 3 1 1 1 1 ( ) n n n n ii i i i i n i i i i i i i a U U U        ， ， ， ， ， (3) 式中： 0a ， 1ia ， 1 2i ia ， 1 2 3i i ia ， 1 2 ni i ia ， ， 均为待定系 数；n 为标准正态随机变量的个数； 1 (n iU ， 2iU ， ， ) ni U 为 n 阶 Hermite 多项式，其计算公式如下： 1 2 ( ) nn i i i U U U ， ， ， T T 1 2 1 1 2 2( 1) e e n n n i i iU U U        U U U U (4) 为了简化，式(3)通常表示为 0 ( ) ( ) P j j j Y c   U U (5) 这里， jc 和 ia ， ( )j U 和 ( )p U 之间存在一一 对应的关系。 ( )j U 和 ( )p U 均构成 Hilbert 空间的 完备正交基，且有 2 i j i ij   ， (6) 式中： ij 为 Kronecker delta 符号；  ， 为内积运 算，定义为 ( ) ( ) ( ) ( ) ( )df g f g W   U U U U U U， (7) 式中：W(U)为权重函数，可表示 T1 2 1( ) (2 ) e n W   U U U (8) 根据式(3)可得，输出随机响应在随机空间自由 度为 M 时，p 阶 Hermite 随机多项式展开的项数 Na 的计算公式[17]为 a ( )! ! ! M pN M p  (9) 在式(2)的基础上，S. S. Isukapalli 等[11]给出了 2 和 3 阶 Hermite 随机多项式展开表达式。对于可靠 度计算来说，当失效概率较小(如小于 10－4)时，2 和 3 阶随机多项式展开可能不能满足计算精度的需 要[1]。此外，对于均匀和  分布等变量来说，2 和 3 阶 Hermite 随机多项式的误差会很大。为了说明， 图 1 给出了均匀分布的概率密度函数曲线的比较。 从图 1 可以看出，2 阶随机多项式展开完全不能逼 近理论概率密度函数曲线，4 阶随机多项式展开可 以近似理论概率密度函数曲线，但是在曲线拐角处 仍然没有办法逼近。相对来说，6 阶随机多项式展 开可以很好地逼近原概率密度函数曲线。由此可 见，发展高阶 Hermite 随机多项式展开的解析表达 式是非常有必要的。为了便于高阶 Hermite 多项式 展开的推广应用，本文推导了 4～6 阶 Hermite 随机 多项式展开的解析表达式 Y4～Y6，其展开表达式分 别为 • 1516 • 岩石力学与工程学报 2010年 x 图 1 均匀分布的概率密度函数曲线的比较 Fig.1 Comparison of probability density functions curves for uniform distribution random variables 2 4 0 1 1 3 4 2 1 1 ( 1) ( 3 ) ( 6 3) n n i i ii i i i n n iii i i iiii i i i i Y a a U a U a U U a U U                  ≈ 2 3 1 1 1 1 ( ) ( 3 ) n n n n ijj i j i ijjj i j i j i j i j j i j i a U U U a U U U U           1 2 2 2 2 1 1 1 2 1 1 1 ( 1) ( ) n n iijj i j i j i j i n n n n n iijk i j k j k ij i j i j k j i j i j i k i a U U U U a U U U U U a U U                      2 1 3 2 1 1 1 ( ) n n n n n n n ijk i j k ijkl i j k l i j i k j i j i k j l k a U U U a U U U U               (10) 5 3 5 4 1 4 2 1 1 2 3 2 3 1 1 ( 10 15 ) ( 6 3 ) ( 3 3 ) n iiiii i i i i n n ijjjj i j i j i i j j i n n iijjj i j i j j j i j j i Y Y a U U U a U U U U U a U U U U U U                    ≈ 2 1 2 1 1 1 3 1 1 ( ) ( 3 ) n n n n iijkl i j k l j k l i j k j l k j i k i l i l j n n n iiijk i j k i j k i j k j j i k i a U U U U U U U a U U U U U U                    1 2 2 2 2 1 1 ( ) n n n ijjkk i j k i j i k i i j k j j i k i a U U U U U U U U          4 3 2 1 1 ( ) n n n n n ijklm i j k l m i j i k j l k m l a U U U U U           (11) 6 4 2 6 5 1 5 3 1 1 ( 15 45 15) ( 10 15 ) n iiiiii i i i i n n ijjjjj i j i j i j i j j i Y Y a U U U a U U U U U U              ≈ 2 4 2 2 2 4 2 1 1 1 3 3 3 3 1 ( 6 3 6 3) ( 3 3 9 ) n n iijjjj i j i j i j j i j j i n n iiijjj i j i j j i i j i j i a U U U U U U U a U U U U U U U U                  1 4 2 1 1 3 2 3 2 1 1 1 ( 6 3 ) ( 3 3 ) n n n iiiijk i j k i j k j k i j k j j i k i n n n iiijjk i j k i k i j k i k i j k j i k j k i a U U U U U U U U a U U U U U U U U U U                   1 2 2 2 2 2 1 2 2 2 2 2 2 2 2 1 3 1 1 ( 1) ( 3 ) n n n iijjkk i j k i j i j i k j i k j k i j k n n n n iiijkl i j k l i j k l i j k j l k j i k i l i a U U U U U U U U U U U U a U U U U U U U U                       ＞ ＞ 1 1 2 2 2 2 1 1 3 2 1 2 1 1 ( ) ( ) n n n n iijjkl i j k l i k l j k l i j i k l k k j l j k i l i n n n n n k l iijklm i j k l m j k l m i j k j l k m l j i k i l i m i a U U U U U U U U U U U U a U U U U U U U U U                       ＞ ＞ ＞ ＞ ＞ 5 4 3 2 1 1 ( ) n n n n n n ijklmo i j k l m o i j i k j l k m l o m a U U U U U U        ＞ ＞ ＞ ＞ ＞ (12) 由式(10)～(12)可以看出，4～6 阶 Hermite 随机 多项式的计算较 2 和 3 阶要复杂得多。 为了求解上述随机多项式的待定系数，常用的 2 阶随机多项式展开 4 阶随机多项式展开 6 阶随机多项式展开 理论概率密度函数 概 率 密 度 函 数 PD F 第 29 卷 第 8 期 李典庆，等. 边坡可靠度分析的随机响应面法及程序实现 • 1517 • 方法是概率配点法。配点法的结果通常取决于配点 的选择，一般来说 p 阶随机多项式展开的配点是通 过 p+1 阶 Hermite 多项式的根来确定输入随机变量 的取值，如 3 阶 Hermite 多项式 33 3H U U  的根 为 0， 3 ，4 阶 Hermite 多项式 4 24 6 3H U U   的根为 3 6  。按照上述方法进行配点，可供 选择的配点数目为 c ( 1)MN p  (13) 一般来说， cN >> aN ，因此选取的配点数只是 所有配点数 cN 中的一部分。选取配点的原则是优先 选择高概率区域内的点[18]，而且在布置配点时应尽 量关于原点对称。选取配点的个数一般要多于待定 系数的个数，S. S. Isukapalli 等[11]建议配点数目取为 待定系数的 2 倍，这样可以平衡每个配点的影响。 上述独立标准正态随机变量的配点经过 Nataf 变换 后得到原始随机变量相应的配点，再将其代入确定 性分析程序可得相应配点对应的响应值向量 Y = 1 2{ }ny y y， ， ， 。最后将该配点和 Y 一起代入式(3)， 得到一个以待定系数 ia 为未知数的线性代数方程。 同理，根据 m 组配点可得 m 个线性代数方程构成的 线性代数方程组，求解这个线性代数方程组可得待 定系数。 上述随机响应面法只能考虑变量相互独立的情 况，对于输入变量是相关非正态变量来说，S. S. Isukapalli 等[11]建议采用如下的方法：设输入随机变 量向量 X 的协方差矩阵为 x ， ijy 与 ijx 的关系为 / ( ) i iij ij x y y x     。再将 y 进行 Cholskey 分解可 得下三角矩阵 L，进而得到相关标准正态随机变量 向量 Y = LZ(Z 为独立标准正态随机变量向量)。最 后可得原始随机变量计算公式 ix  i ix x iy  。值得 指出的是，S. S. Isukapalli 等[11]还特别说明其并没有 采用算例验证上述方法的正确性，经过简单分析就 可以发现上述方法只适用于相关正态变量，对于相 关非正态变量来说，上述方法在大多数情况下误差 会很大。为了将随机响应面法拓展应用于相关非正 态变量的可靠度分析问题，本文提出采用 Nataf 变 换方法来建立原始的相关非正态变量和变换后的独 立标准正态随机变量间的映射关系，从而有效地解 决了随机响应面法在分析输入变量相关时的可靠 度分析问题，下面简要介绍 Nataf 变换方法基本原 理。 3 Nataf 变换方法 Nataf 变换[15]利用变量的边缘概率密度函数、 相关系数矩阵以及 n 维标准正态分布概率密度函 数，通过 Gaussian Copula 函数[19]构造变量联合概 率密度函数，具体步骤如下：定义 n 维随机变量向 量 T1 2{ }nX X X ， ， ，X ，其相关系数矩阵  ( )ij n n  。根据等概率变换原则可得 X 和 Y 中的各个 变量有如下函数关系： 1 ( ) ( ) ( 1 2 ) [ ( )] i i i X i i X i y F x i n y F x        ， ， ， (14) 式中： ( ) iX F  为任意变量的边缘概率密度函数， 1( )   为 ( )  的反函数。 根据 Nataf 变换理论，利用隐函数求导法则可 推导出变量 X 的联合概率密度函数为 1 21 2 0 1 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) nX X X n n n f x f x f x f x y y y y           ，X (15) 式中： ( ) iX f x 为变量 iX 的边缘概率密度函数， (n y ， 0 ) 为 n 维标准正态变量向量的概率密度函数。一 般地，将变量 X 的概率分布模型称为 Nataf 分布。 相关非正态变量当量正态化过程一般会改变原始变 量间的相关性，也就是说当量正态化后的标准正态 随机变量的相关系数 0ij 和原相关非正态变量的相 关系数 ij 一般都是不相等的。根据相关系数的定义 及式(14)和(15)可得变量 X 的相关系数 ij 与等效标 准正态随机变量 Y 的相关系数 0ij 有以下关系： 1 [ ( )] i i i i x ij x F y              1 2 0 [ ( )] ( )d dj j i i x i j ij i j x F y y y y y          ， ， (16) 式中： 2 0( )i j ijy y ， ， 为相关系数是 0ij 的二维标准 正态分布随机变量的联合密度概率函数。 当 iX 和 jX 的边缘概率分布函数及其相关系数 ij 已知时，通过求解式(16)所示的非线性方程就可 以确定等效相关系数 0ij 。A. Der Kiureghian 和 P. L. Liu[15]已经给出了 2 类共计 10 种概率分布间等效相 关系数经验计算公式，但并没有给出本文要用到的 截尾指数分布和正态分布间的经验公式。为此，本 文编写了基于二维 Gaussian-Hermite 积分的计算程 • 1518 • 岩石力学与工程学报 2010年 序 CoefNataf 来求解式(16)中等效相关系数 0ij 。 对标准正态随机变量向量 Y 的相关系数矩阵 0 进行 Cholesky 分解或特征值分解可得 T 0 0 0 Γ Γ (17) 式中： 0Γ 为 0 的分解矩阵。利用 0Γ 可将相关标准 正态随机变量向量 Y 转换为独立的标准正态随机变 量向量 U： 1 0 U Γ Y (18) 至此完成了Nataf变换的正变换过程。 因为随机响应面法的输入变量是标准正态随机 变量，首先对标准正态随机变量进行配点，然后再 通过一定的转换关系映射到原始随机变量，也即将 独立标准正态随机变量向量转化为相关非正态随机 变量向量，显然这个过程对应的变换就是 Nataf 变 换的逆变换过程。将式(18)进行变换可得 0Y Γ U (19) 然后，采用式(14)可得 1[ ( )] i i ix F y (20) 这样就可以实现标准正态随机变量配点到原始 相关非正态变量配点的转换过程。 4 随机响应面法计算程序WHUSRSM 由上述介绍可以看出，采用随机响应面法计算 可靠度比较复杂，必须借助一定的计算程序才能进 行。为此，本文编写了基于 C#语言的随机响应面法 计算程序 WHUSRSM[20]，其结构如图 2 所示。该程 序主要包括定义问题(含输入随机变量个数和 Hermite 随机多项式展开阶数)、定义随机变量(变量均值、 标准差以及其他参数)、配点(独立标准正态变量 U 空间、原始的相关非正态变量 X 空间配点及相应的 输出响应值)、求解和结果显示(随机多项式的待定 系数、拟合的随机响应面、输出响应的统计参数如 均值、标准差、变异系数、峰度系数、偏度系数等) 显示 5 个部分组成。国内外也开发了类似的计算程 序如 WebSRSM[11]和 CSRSM[13]。与上述 2 个计算程 序相比，WHUSRSM 具有以下 3 个优点： (1) WHUSRSM 采用 Nataf 变换成功地解决了相关非正 态变量的配点及其可靠度计算问题，从而极大地拓 展了随机响应面法的应用范围；(2) 输入和输出随 机变量最多可达 20 个，而且 Hermite 随机多项式可 以展开到 4～6 阶，在变量个数和随机多项式展开的 阶数方面具有明显的改进；(3) WHUSRSM 可以计 算 12 种概率分布的随机变量，如截尾指数分布、 威布尔分布等，变量分布类型更加广泛。 图 2 WHUSRSM 分析软件结构示意图 Fig.2 WHUSRSM outline structure guide 5 算 例 本文考虑的岩质边坡示意图如图 3 所示[5]，该 模型具有一定的代表性，也是最典型的。然而天然 边坡并非如此规则，为了简化边坡安全系数的计算， 计算中采用的假定与 E. Hoek 和 J. Bray[21]中的一样， 如：计算中假定滑面及张裂缝的走向平行于坡面， 仅有一条张裂缝且张裂缝垂直，水沿张裂缝底进入 滑动面渗漏，张裂缝底与坡趾间长度内的水压力按 线性变化至 0(三角形分布)等。算例中确定性变量 的取值如下：岩石和水的容重分别为   2.6×104 N/m3， w  1×104 N/m3， f  50°， p  35°，T = 0，  0，H = 60 m。随机变量及分布参数如表 1 所示。表 1 中：c， ，z 服从正态分布；r， 服从 截尾指数分布，r 的下界为 0，上界为 1， 的下界 为 0，上界为 0.16。边坡安全系数的计算公式[5]为 p p p tan (sin cos ) cos sin cA NFS W V T           (21) 显示区域 状态栏 主菜单栏 快捷菜单栏 导航菜单栏  Z W W T U V b B H 水压力分布 拉裂面 底滑面 Z w p f 第 29 卷 第 8 期 李典庆，等. 边坡可靠度分析的随机响应面法及程序实现 • 1519 • 图 3 边坡示意图[5] Fig.3 Illustration of slope[5] 表 1 随机变量及分布参数 Table 1 Random variables and distribution parameters 取值 类型 黏聚力 c /kPa 内摩擦 角(°) 单宽拉裂面 面积 z/m2 张缝中充水深 度系数 r 水平地震加速 度系数 均值 100 35 14 0.5 0.08 标准差 20 5 3 相应的功能函数为 1G FS  (22) 式(21)，(22)中：A 为底滑面的面积(单宽)， A  (H  p) / sinz  ；z 为拉裂面的面积(单宽)， z   f p1 cos tanH   ； N 为作用在底滑面上的正压 力，N   p p p(cos sin ) sinW U V       cosT  ； W 为滑块自重， 2 2 p0.5 [(1 ( / ) )cotW H z H    fcot ] ；U 为作用在底滑面上的水压力合力， w0.5U rzA ；r 为张缝中充水深度系数 w( / )r z z ； V 为作用在拉裂面上的水压力合力， 2 2w0.5V r z 。 采用 WHUSRSM 计算程序可得算例的计算结 果，图 4 给出了不同计算方法的功能函数累积分布 函数曲线的比较，相关系数 0.5c z r   ， ， 。这 里比较了 Isukapalli 方法和本文提出的方法，其中相 关系数矩阵相关性分解又分别考虑了乔列斯基分解 和特征值分解，这样共有 4 种计算方法。这里需要 指出的是4阶随机多项式展开(378个配点)的计算精 度已经满足要求，没有必要采用更高阶的随机多项 式展开。从图 4 可以看出，本文方法的 CDF 曲线 和 100 万次 MCS 的 CDF 曲线几乎完全重合，而 Isukapalli 方法的 CDF 曲线和 100 万次 MCS的 CDF 曲线有明显的差别。将 100 万次 MCS 结果视为精 确解，本文方法和 Isukapalli 方法的最大相对误差分 别为 33.3%和 19 621%，但是本文关注的只是功能 函数小于 0 的失效概率，由 CDF 曲线可得本文方 法、Isukapalli 方法以及 100 万次 MCS 方法的失效 概率分别为 0.023 0，0.179 0，0.022 4，本文方法和 Isukapalli 方法的相对误差分别为 2.7%和 699%，显 然本文方法的计算精度很高，而 Isukapalli 方法的结 果是值得商榷的。此外，无论是本文方法还是 Isukapalli 方法，乔列斯基分解和特征值分解的结果 几乎一样，可见相关系数矩阵的相关性分解方法对 随机响应面法的可靠度计算结果几乎没有影响，因 图 4 不同计算方法的功能函数累积分布函数曲线的比较 Fig.4 Comparison among cumulative distribution function curves of performance functions using different calculation methods 此下面计算都采用乔列斯基分解方法。 图 5，6 分别给出了概率密度函数(PDF)和累计 分布函数(CDF)曲线的比较。为了比较，配点数目 均取为随机多项式待定系数数目的 3 倍。由图中可 以看出，2 阶随机多项式展开 63 个配点的 PDF 和 CDF 曲线与 100 万次 MCS 结果有一些差别，3 阶 随机多项式展开 168 个配点的 PDF 和 CDF 曲线与 100 万次 MCS 结果也存在非常明显的差别，4 阶随 机多项式展开 378个配点的PDF和CDF曲线与 100 万次 MCS 的 PDF 和 CDF 曲线基本重合。此外，虽 然图 5中不同阶数随机响应面法PDF曲线存在一定 的差别，但是功能函数均值的估计值几乎完全一样， 均值都是 0.353。这说明随机响应面法在估计输出响 应的均值方面具有一定的稳健性。由此可得，2，3， 4 阶随机多项式展开以及 100 万次 MCS的失效概率 分别为 0.026 0，0.178 0，0.023 0，0.022 4，可见失 效概率存在明显的差别，尤其是 3 阶随机多项式展 开的失效概率和其他 3 个失效概率值根本不在一个 数量级上。这进一步说明虽然不同阶数的随机响应 面法能够精确地估计功能函数的均值，但是它们并 不能准确地估计功能函数小于 0 的概率(失效概率)。 －1.0 －0.5 0.0 0.5 1.0 1.5 1×10－3 功能函数 G 累 计 分 布 函 数 CD F 4 阶随机多项式展开(378 个配点) 乔列斯基分解特征值基分解 Nataf+乔列斯基分解 Nataf+特征值基分解 MCS(100 万次)模拟 1×10－4 1×10－1 1×10－2 1×10－0 • 1520 • 岩石力学与工程学报 2010年 原因在于随机变量 PDF 曲线尾部的形状对失效概 率计算有很大的影响，图 5 中 3 阶随机多项式展开 的PDF曲线尾部和MCS结果存在非常明显的差别， 从而导致了图 6 中 CDF 曲线尾部和 MCS 结果的差 别更加明显。 功能函数 G 图 5 概率密度函数曲线的比较 Fig.5 Comparison among probability density functions curves 功能函数 G 图 6 累积分布函数曲线的比较 Fig.6 Comparison among cumulative distribution function curves 为了进一步研究随机响应面法的计算精度，图 7 给出了不同计算方法可靠指标的比较。可以看出，4 阶随机多项式展开的结果和 100 万次 MCS 结果保 持较好的一致性，尤其是当变量间相关性较低时。 传统的 FORM 方法可靠指标和 MCS 方法可靠指标 存在明显的差别，而且这种差别随着变量间负相关 性的增加而显著增大。以图 7(b)中同时考虑 c  ， 和 z r ， 为例，当负相关系数大于 0.6 时，4 阶随机多 项式展开的可靠指标和 MCS 可靠指标最大绝对误 差为 0.023 2，而 FORM 方法可靠指标和 MCS 可靠 指标最大绝对误差为 0.173。可见功能函数(式(22)) 只计算 378 次(4 阶随机多项式展开)的效果和计算 100 万次(MCS 方法)的效果几乎相同。但是当 c  ， = z r ， = 0.96 时，与 100 万次 MCS 结果相比， 4 阶随机多项式展开和 FORM 方法可靠指标绝对误 差分别为 0.386 和 0.645。可见，当随机变量间负相 关性非常高时，4 阶随机多项式误差较大，此时采 相关系数 (a) 仅考虑 c  ，  相关系数 (b) 同时考虑 c  ， 和 z r ， 图 7 不同计算方法可靠指标的比较 Fig.7 Comparison among reliability indices using different calculation methods 用更高阶的随机多项式展开或者更多的配点也 许可 商标使用许可商标使用许可商标使用许可商标使用许可商标使用许可 以减小上述误差。总体来说，基于 4 阶随机多项式 展开的随机响应面法的计算精度要优于传统的 FORM 方法，同时其计算效率要优于传统的 MCS 方法，而且在大多数情况下其结果和 MCS 保持较 好的一致性。 传统的边坡可靠度计算的 FORM 方法不能分 析输入变量和输出响应量间的相关关系，然而这在 随机响应面法中是非常容易实现的。图 8 给出了安 全系数和输入变量间的相关性。可以看出，边坡安 全系数与黏聚力和内摩擦角都呈正相关关系，这是 概 率 密 度 函 数 PD F 2 阶随机多项式展开(63 个配点) 3 阶随机多项式展开(168 个配点) 4 阶随机多项式展开(378 个配点) MCS(100 万次)模拟 c， = z，r =－0.5 累 积 分 布 函 数 CD F c， = z，r =－0.5 2 阶随机多项式展开(63 个配点) 3 阶随机多项式展开(168 个配点) 4 阶随机多项式展开(378 个配点) MCS(100 万次)模拟 可 靠 指 标 4 阶随机多项式展开 FORM MCS 可 靠 指 标 4 阶随机多项式展开 FORM MCS c， = z，r 第 29 卷 第 8 期 李典庆，等. 边坡可靠度分析的随机响应面法及程序实现 • 1521 • 因为抗剪强度参数是抗力作用的变量，其量值的增 加对于边坡稳定性是有利的。此外还可以看出，内 摩擦角对边坡稳定性的影响大于黏聚力的影响。相 反，边坡安全系数与张裂缝中充水深度系数和水平 地震加速度系数都呈明显的负相关关系，而且两者 的负相关性均较大，这是因为边坡张裂缝中水压力 黏聚力/(10 kPa) (a) 黏聚力 内摩擦角/(°) (b) 内摩擦角 张裂缝中充水深度系数 (c) 张裂缝中充水深度系数 水平地震加速度系数 (d) 水平地震加速度系数 图 8 安全系数和输入变量间的相关性 Fig.8 Correlation between factors of safety and input variables 和水平地震力都是荷载作用的变量，其量值的增加 对于边坡稳定性是非常不利的。 为了研究滑面抗剪强度参数变化对边坡安全系 数的影响，图 9 给出了不同强度折减系数对安全系数 累积分布函数曲线的影响。这里采用 4 阶随机多项 式展开 378 个配点和 c  ， =－0.5。由图 9 可见，滑面 抗剪强度对边坡安全系数有非常明显的影响，如当  ，，，，，时，相应的失效概 率分别为 1.000，0.960，0.690，0.320，0.120，0.038， 这意味着当抗剪强度降为原来一半时，该边坡肯定 会失稳；当抗剪强度降为原来 80%时，边坡失效概 率将变为原来的 8.4 倍。因此，在岩质边坡地质勘 测中要尽可能获得诸如断层、节理等准确的资料， 从而提高边坡可靠度分析结果的准确性。 安全系数 FS 图 9 不同强度折减系数对安全系数累积分布函数曲线的影响 Fig.9 Effects of different strength reduction factors on cumulative distribution function curves of factors of safety 为了定量地研究边坡稳定性分析的安全系数和 可靠度的关系，有必要分析安全系数和失效概率或 安 全 系 数 FS 相关系数 = 0.17 抽样数目 N = 5 000 安 全 系 数 FS 相关系数 = 0.48 抽样数目 N = 5 000 安 全 系 数 FS 相关系数 =－0.48 抽样数目 N = 5 000 安 全 系 数 FS 相关系数 =－0.49 抽样数目 N = 5 000 累 积 分 布 函 数 CD F  = 1.0  = 0.9  = 0.8  = 0.7  = 0.6  = 0.5 • 1522 • 岩石力学与工程学报