楔形体边坡的系统可靠度分析

谭晓慧，宋传中，侯晓亮

（合肥工业大学 资源与环境工程学院，合肥230009）

边坡的稳定性分析是岩土工程中十分重要的研究课题。由于岩土体的地质条件、环境条件、荷载条件等因素的不确定性，使得边坡工程中存在大量的变异性［1－3］。尤其是对于楔形体岩质边坡，不仅存在材料参数的变异性，还存在结构面几何参数的变异性。而且，楔形体边坡有多种失稳模式，随着参数均值及变异性的不同将可能发生不同模式的失稳。因此，对于楔形体岩质边坡需要进行边坡稳定的系统可靠度分析。

文献［4］将多滑面边坡看成一个串联系统，进行了多滑面岩质边坡的系统可靠度分析，并研究了参数的敏感性。文献［5］视边坡为4种失稳模式组成的串联系统，采用电子表格法求解了每一种失稳模式及整个边坡体系的可靠度。文献［6］对平面滑动及楔形体滑动边坡采用蒙特卡罗模拟法进行了系统可靠度分析。文献［7］进行了平面滑动岩质边坡的系统可靠度分析，考虑了拉张裂缝位置的对边坡稳定性的影响，并且视整个边坡的失稳为串并联体系，采用线性编程法求解了系统失效概率的范围值。文献［8］在进行边坡稳定的有限元可靠度分析时，视单一失稳模式是由滑动面上的单元组成的并联体系，采用逐步等效线性化求交法来计算单一失稳模式的可靠指标，再采用界限法估算了由各失稳模式组成的串联体系的可靠指标。

在系统可靠度分析方面，文献［9－10］采用系统可靠度分析的界限法，进行了挡土墙的系统可靠度分析。文献［11］采用逐步等效平面法分析了挡土墙的抗倾覆及抗滑移可靠度。文献［12］对系统可靠度分析的条件余量法进行了改进，提高了计算精度。

本文将楔形体边坡视为由4个并联子系统组成的串联系统，提出了在进行系统可靠度分析时对非线性功能函数进行线性化的新方法，通过算例分析验证了该方法的精度，并采用条件余量法及界限法进行了楔形体边坡的系统可靠度分析。

1 楔形体边坡稳定的定值法分析

如图1所示某楔形体边坡，坡顶倾角为Ω，坡面倾角为α，二者倾向相同，其交线TX 与水平面的距离为h。记节理面及分别为滑面1及滑面2，其走向与坡面走向之夹角分别为θ1、θ2；倾角分别为δ1、δ2；内聚力及内摩擦角分别为c1、c2、φ1及φ2；节理面交线与水平面的夹角为ε。当Ω＜ε＜α时，楔形体可能会沿着滑面1或滑面2滑动。文献［4－5］对此问题进行了深入研究，求得了各种情况下边坡安全系数的解析解，具体如下。

图1 楔形体边坡示意图

1.1 楔形体同时沿滑面1及滑面2滑动（模式1）

式（1）成立的前提是式中tanφ1及tanφ2前面的系数项为正，即：

式中：Gw1及Gw2分别为归一化的水压力参数，分别为节理面1、2上的水压力；Sγ为岩石的比重，及γw分别为岩石及水的容重。其它符号的含义如下：

1.2 楔形体只沿节理面1滑动（模式2）

若作用在滑面2（面B′CO）上的净法向力是上浮力时，楔形体B′DOC只沿着滑面1（面B′DO）滑动，其安全系数表达式如下：

式中：

式（10）成立的前提是：

1.3 楔形体只沿节理2滑动（模式3）

此式成立的前提条件：

1.4 楔形体浮起（模式4）

若作用在滑面1、2上的净法向力都是上浮力时，则楔形体B′DOC由于水压力的作用而浮起，失去与两个节理面的接触，此时，记边坡的安全系数Fs4＝0。此情况出现的条件是：

2 楔形体边坡稳定的可靠度分析

2.1 楔形体边坡的系统组成

由上述分析可见，图1所示的楔形体边坡失稳模式有4种，其中任一种情况发生时边坡都将失稳。因此，整个边坡可以视为由这4种失稳模式组成的串联系统。

对于楔形体边坡的每一种失稳模式，它又分别需要满足2个限制条件，并且当Fsi≤0（i＝1，2，3，4）时才会失稳。亦即当3个条件同时满足时，该失稳模式才可能发生。因此，楔形体边坡的每一种失稳模式又可以视为3个元件组成的并联系统。

在进行可靠度分析时，设基本变量为X，分别记各元件的功能函数如下：

则整个楔形体边坡的系统组成可由图2表示。

由图2可见，整个楔形体边坡可视为由4种失稳模式组成的串联系统，而每一种失稳模式又可视为由3个元件组成的并联子系统。因此，为了求解整个楔形体边坡的可靠度，需要先对每一元件进行可靠度分析。

图2 楔形体边坡失稳模式的系统组成

2.2 元件的可靠度分析

元件的可靠度分析是系统可靠度分析的基础，有关元件可靠度分析的研究成果较为丰富，分析方法也较多，如：一阶可靠度分析方法、响应面法、蒙特卡罗模拟法等［13－15］。一阶可靠度分析方法的实质是采用通过功能函数Z＝0上某一点（验算点）的超切平面来代替原来的非线性功能函数，使非线性功能函数线性化，其验算点坐标x＊及可靠指标β可通过迭代法求解：

式中：μX、σX及CX分别为随机变量X的均值向量、均方差矩阵及协方差矩阵；▽g（x＊）为验算点x＊处功能函数对基本变量的偏导数向量；αX为敏感性系数（方向余弦）向量。

由于可靠指标是标准正态空间中坐标原点到极限状态曲面的最短距离，因此亦可采用最优化方法求解可靠指标。

蒙持卡罗模拟法（MCSM）通过随机模拟和统计试验来求解结构的可靠度。在目前的结构可靠度分析中，它被认为是一种相对精确的方法［13］。

可靠指标β与失效概率Pf具有如下对应关系：

式中：Ф为标准正态分布函数。

2.3 系统的可靠度分析

设整个结构存在m个失稳模式，包含n个基本变量，则并联系统和串联系统的失效事件分别是各失稳模式失效事件的交事件及和事件：

准确求解系统的失效概率需要高维积分，因此一般采用简化的近似方法求解系统的失效概率及可靠指标，即：先对各功能函数分别进行线性化，再用超切平面围成的失效域代替原极限状态曲面围成的失效域进行求解。经过近似处理后，并联系统及串联系统的失效概率分别为

式中：ZLi为Zi的线性化函数；Фm为m维标准正态分布函数，其简化计算方法很多，由于条件余量法简单、精度高［12］，本文采用该方法来求解Фm；β为可靠指标向量；ρ为各失稳模式间的相关系数矩阵，它可由各失稳模式线性化功能函数的敏感性系数向量求得。对于第i及第j失稳模式，其相关系数为：

式中：α为验算点处的敏感性系数，可由一阶可靠度分析方法求得。

对于串联系统，亦可采用界限法来近似求解系统失效概率的界限值：

该不等式左右端分别对应于各失稳模式完全相关及完全不相关的情况。由于式（30）中不包含各失稳模式间的相关系数矩阵ρ，因而可大大简化系统可靠度分析的过程。

2.4 非线性功能函数的线性化

为了求解各线性化功能函数间的相关系数矩阵ρ，需要对系统中各元件的非线性功能函数进行线性化。对功能函数进行线性化时，线性化点位的选择非常重要。对于并联系统，在不同的点进行线性化，得到的线性化函数围成的新失效域的位置及大小也有所不同；而对于串联系统，线性化点位的选择对计算结果的影响相对较小［15］。因此，本文主要研究并联系统中非线性功能函数的线性化方法。

文献［7］在进行平面滑动边坡的并联子系统可靠度分析时，分别对各元件的功能函数求解各自的验算点，再在这些验算点处对相应的功能函数分别进行线性化。此时，线性化函数围成的新失效域与原失效域可能有一定差别；并联系统中各失稳模式的非线性化程度越大，新失效域与原失效域的差别也越大［15］。

文献［5］在进行楔形体滑动边坡的并联子系统可靠度分析时，将可靠指标的求解视为非线性有约束最优化问题，通过Excel电子表格的优化求解功能同时求解了各功能函数的联合验算点及相应的可靠指标，并将此可靠指标视为整个并联系统的可靠指标。该方法计算过程较为简单，但是，由可靠指标的几何意义可知，该方法的实质是将标准正态空间中与坐标原点至验算点连线相垂直的超切平面视为线性化平面，因此新失效域与原失效域可能相差较大。

文献［15］的方法如下：先求解并联子系统中各功能函数的联合验算点，再在联合验算点处对各元件的功能函数进行线性化。若联合验算点只由并联子系统中的几个元件决定，则只对这几个元件的功能函数进行线性化。

由于文献［15］只对经过联合验算点的功能函数进行线性化而忽略了其它功能函数的影响，因此线性化后的新失效域与原失效域可能有一定的差别。本文对此进行改进，提出一种新的线性化方法，该方法分两步：1）求解联合验算点；2）当功能函数经过联合验算点时，在联合验算点对该函数进行线性化；否则，在功能函数各自的验算点进行函数的线性化。这样既可考虑联合验算点对失效域的影响，又能考虑所有功能函数对失效域的影响，从而能使新失效域与原失效域尽可能一致。

为了研究文献［5］的方法（记为方法1）、文献［15］的方法（记为方法2）及本文方法的优缺点，分别采用这些方法进行功能函数的线性化，再采用条件余量法进行并联子系统的可靠度分析。其中，联合验算点的求解由如下的非线性有约束最优化问题确定：

式中：n为基本变量的个数；Np为并联子系统中元件的个数（Np＝3）；U 代表标准正态空间；满足式（31）的点uj（j＝1，2，…，n）即为验算点。

3 算例分析

3.1 边坡的几何外形及材料参数

某高速公路挖方边坡由中风化粉砂岩组成，边坡高24 m，拟开挖坡面产状为255°∠53°。坡体发育2组节理，其产状分别为220°∠50°及340°∠68°。这两组节理面形成楔形体，对边坡稳定性影响较大。

在进行楔形体边坡的可靠度分析时，取2个节理面的水压参数及强度参数分别相等；视节理面走向与坡面走向之夹角 （θ1、θ2）、节理面倾角（δ1、δ2）、节理面的水压力参数及强度参数（Gw1、tanφ1及c1／γh）为基本变量，其统计参数见表1，表中N及LN分别代表正态及对数正态分布。

表1 楔形体边坡的几何及材料参数

为了形象化地分析系统中各元件功能函数间的关系，先假设基本变量而将其它变量视为定值。在此基础上，再研究基本变量为7个，即的情况。以下分别简记这两种分析为2变量情况及7变量情况。在进行可靠度分析时，单一元件的可靠指标由一阶可靠度分析方法求解；并联系统的失效概率由条件余量法求解；联合验算点由科学计算语言MATLAB的非线性有约束最优化函数fmincon求解；各种情况下可靠指标及失效概率的精确值由蒙特卡罗模拟法（MCSM）求得。所有计算程序均由MATLAB编写。

3.2 2变量情况

3.2.1 并联子系统的可靠度分析 由前述分析可知：失稳模式1～4分别是由3个元件组成的并联子系统，而求解并联子系统的失效概率时对功能函数进行线性化的方法又有多种，因此，表2列出了4种失稳模式下采用不同方法求得的可靠指标β及失效概率Pf。为了与可靠指标的精确值进行对比，表中亦列出了106次蒙特卡罗模拟法的计算结果。

表2 并联子系统的可靠指标及失效概率（2变量）

由表2可见：除了模式1以外，方法1、2的可靠指标及失效概率与其对应的精确解有一定差别，而本文方法的结果则与精确解非常接近。前两种方法的误差可以由图3～6预以说明。

图3显示了模式1的3个元件在标准正态空间中的极限状态曲线。线AB、CD、EF分别代表g1＝0、g5＝0及对应的失效域分别在这3条极限状态曲线的上部、下部及下部，这3个区域无交集，即无失效域，因此破坏概率Pf＝0。对于此特殊情况，各种方法的求解结果相同。

图3 极限状态曲线（模式1）

图4显示了模式2的3个元件在标准正态空间中的极限状态曲线，它们在原点附近近似于直线。线AB、CD、EB分别代表g2＝0、g6＝0及g7＝0；g2≤0、g6≤0及－g7≤0对应的失效域分别在这3条极限状态曲线的上部、上部及下部，它们围成的失效域为一个近似于三角形的区域PFB。图中P点为联合验算点；该点到坐标原点的距离OP即为方法1的可靠指标。连接点OP，过P点作线OP的垂线GH。由可靠指标及失效概率的几何意义可知：方法1对应的失效域对应于线GH的右上部分区域，它大于真实失效域PFB。因此，方法1的失效概率Pf偏大，可靠指标β偏小。

由于联合验算点P是极限状态线g6＝0及g7＝0的交点，因此，方法2分别取过P点的g6及g7的线性化函数（线PD、线PB）所围成的区域BPD作为失效域。可见，该失效域也大于真实的失效域PFB，但小于方法1的失效域，故方法2的失效概率Pf大于真实值但小于方法1的相应值，而可靠指标β值则小于真实值而大于方法1的相应值。

图4 极限状态曲线（模式2）

本文方法与方法2的区别是：对于g6及g7分别在联合验算点进行线性化；对于g2在其自己的验算点（A点）进行线性化，因而本文方法得到的失效域为图中PFB，它与真实失效域几乎重合，故其对应的失效概率及可靠指标与真实解亦非常接近。

图5显示了模式3在标准状态空间中的极限状态曲线，AB、CD、EB分别代表g3＝0、g5＝0及g8＝0。真实失效域为这三者围成的区域FGDKJ，该失效域边界上与坐标原点距离最近的点P为联合验算点，它亦是g5＝0的验算点。另两个功能函数的验算点分别为B及K。功能函数g3、g5及g8在其各自的验算点处线性化的极限状态曲线分别为AB、CD及ED。

方法1视联合验算点与坐标原点的距离OP为可靠指标，切线CPD以上的区域为失效域。它大于真实失效域FGDKJ。因此，该方法的失效概率Pf偏大，可靠指标β偏小。

由于联合验算点P只满足g5＝0，因此，方法2只取过P点的g5的线性化函数（CPD）以上的区域作为失效域。因此，方法2与方法1的失效域相同，故二者的可靠指标及失效概率亦相同。

图5 极限状态曲线（模式3）

本文方法视功能函数g3、g5及g8在其各自的验算点处线性化极限状态曲线AB、CD及ED围成的区域AHDE为新失效域，它与原失效域FGDKJ比较一致，故所求失效概率及可靠指标也与真实解比较接近。当然，由于功能函数的非线性性，新失效域与原失效域不完全一致。

图6显示的是模式4在标准状态空间中的极限状态曲线，由于Z4恒小于0，因此图中无该功能函数的极限状态曲线。图中AB、CE分别代表g7＝0及g8＝0。真实失效域为CDB以上的区域。P点为该失效域边界上离原点距离最近的点，即联合验算点，它亦是的g8＝0的验算点；G点是g7＝0的验算点。由前述分析可知，方法1及方法2的失效域均为FDE以上的区域，它远大于真实的失效域；而本文方法的失效域为FDB以上的区域，它与真实失效域非常接近。因此，本文方法求得的失效概率及可靠指标优于方法1及方法2的相应解。

图6 极限状态曲线（模式4）

3.2.2 整个边坡系统的可靠度分析

整个楔形体边坡的失稳可视为4个失稳模式的串联系统，因此可以采用串联系统的可靠度分析方法来求解边坡的失效概率及可靠指标。由图4～6可知：各并联子系统的失效域是由各元件的极限状态曲线围成的区域；尽管各元件的极限状态曲线可能是线性的，但由这些元件的极限状态曲线所围区域的边界是高度非线性的；而且，在联合验算点处的并联子系统的极限状态曲线还可能会出现明显的转折点（见图4）。因此，无论采用何种方法对并联子系统的极限状态函数进行线性化，都将引起较大的误差。为了避开线性化各并联子系统的功能函数所带来的误差，本文采用宽限法来求解串联系统的失效概率及可靠指标。以上述各种方法求得的并联子系统的失效概率为基础，由宽限法求得的整个系统的失效概率及可靠指标见表3。与MCSM的结果相比可知：基于本文方法求解并联子系统的失效概率后，再采用宽限法求得的系统失效概率较为准确，而基于方法1及方法2求得的失效概率偏大，可靠指标偏小。

表3 串联系统的失效概率及可靠指标（2变量）

3.3 7变量情况

上节进行了2变量情况下楔形体边坡的可靠度分析，并用图形的形式直观地显示了分析结果。为了更真实地反映边坡中各变量对可靠度分析结果的影响，再对7变量的情况进行分析，此时，。由于基本变量数较多，无法用图形直观地表示各极限状态曲线间的关系，只能列表分析。表4～5为7变量情况下的可靠度分析结果。

表4 并联子系统的失效概率及可靠指标（7变量）

续表4

表5 串联系统的失效概率及可靠指标（7变量）

根据式（1）～（15）可以预测：当基本变量数较多时，各元件极限状态曲线的非线性化程度将更高，由这些极限状态曲线围成的失效域将更复杂。此时，由表4～5可见：方法1～2的计算结果与精确解有较大的误差，但本文方法的计算结果仍与精确解一致。

由表4还可以看出：在楔形体边坡的4种失稳模式中，模式4的可靠指标最小，失效概率最大；而且，该最小可靠指标小于目标可靠指标［1］，因此，该边坡的失稳模式是楔形体同时与两个滑面失去接触，由于水的浮托力作用而发生破坏。但是，定值法分析结果表明：该边坡的安全系数为Fs＝2.070 6，可见该边坡相对稳定；而且，在4种失稳模式中，模式1的安全系数最小，即边坡最有可能沿着滑面1及滑面2同时滑动。二种分析方法不一致的原因是定值法分析时未考虑基本变量的不确定性，因此定值法分析结果偏于安全。由表3及表5亦可发现在相同的情况下，分析时考虑的基本变量数越多，得到的可靠指标就越小，这是因为考虑了更多的不确定性对边坡稳定性的影响。

4 结 论

将楔形体边坡的失稳事件视为由4个失稳模式组成的串联系统，再将这4个失稳分别视为由3个元件组成的并联子系统，采用系统可靠度分析的方法进行了各失稳模式及整个边坡的可靠度分析。在求解并联子系统的失效概率时，提出了对非线性功能函数进行线性化的方法。

实例计算表明：将楔形体边坡的失稳事件视为串并联系统、采用系统可靠度分析的方法进行楔形体边坡的可靠度分析是可行的。本文提出的线性化方法能较好地模拟并联系统的真实失效域，具有较高的求解精度。并联子系统形成的失效域边界是高度非线性的，在进行串联系统的可靠度分析时若对各并联子系统的功能函数再次进行线性化会引起较大的误差，此时适于用宽限法求解整个系统的失效概率及可靠指标。可靠度分析结果表明：该楔形体边坡的最可能失稳模式是：由于水的浮托力的作用，滑体与2个滑面失去接触而发生滑动；而定值法的分析结果则表明滑体将同时沿着两个滑面发生滑动，这是因为定值法的分析结果未考虑基本变量的变异性。因此，在进行楔形体边坡的设计及治理时必须采用可靠度分析的方法，以便全面考虑各种参数的变异性对失稳模式及边坡稳定性的影响。

本文以楔形体边坡为例进行了系统可靠度分析，但该方法不仅适用于楔形体边坡的系统可靠度分析，还适用于其它结构或工程的系统可靠度分析。

［1］Srivastava A，Babu G L S，Haldar S.Influence of spatial variability of permeability property on steady state seepage flow and slope stability analysis ［J］.Engineering Geology，2010（110）：93－101.

［2］Suchomel R，Mašín D.Comparison of different probabilistic methods for predicting stability of a slope in spatially variable cφsoil［J］.Computers and Geotechnics，2010，（37）：132－140.

［3］Wang Y，Cao Z J，Au S K.Efficient Monte Carlo Simulation of parameter sensitivity in probabilistic slope stability analysis［J］.Computers and Geotechnics，2010（37）：1015－1022.

［4］谭晓慧.多滑面边坡的可靠性分析［J］.岩石力学与工程学报.2001，20（6）：822－825.Tan X H.Reliability analysis on a slope with several slide surfaces ［J］.Chinese Journal of Rock Mechanics and Engineering，2001，20（6）：822－825.

［5］Low B K.Reliability analysis of rock wedges［J］.Journal of Geotechnical and Geoenvironmental Engineering，1997，123（6）：498－505.

［6］Park H J，West T R，Woo I.Probabilistic analysis of rock slope stability and random properties of discontinuity parameters，Interstate Highway 40，Western North Carolina，USA ［J］.Engineering Geology，2005，79：230－250.

［7］Jimenez－Rodrigueza R，Sitar N.Rock wedge stability analysis using system reliability methods［J］.Rock Mechanics and Rock Engineering，2007，40（4）：419－427.

［8］吴震宇，陈建康，许唯临，等.岩质边坡稳定的体系可靠度分析及工程应用［J］.四川大学学报：工程科学版，2008，40（2）：32－37.Wi Z Y，Chen J K，Xu W L，et al.Systematic reliability analysis of rock slope stability and its engineering application［J］.Journal of Sichuan University：Engineering Science Edition，2008，40（2）：32－37.

［9］Low B K，Zhang J，Tang W H.Efficient system reliability analysis illustrated for a retaining wall and a soil slope ［J］.Computers and Geotechnics，2011，38：196－204.

［10］Zevgolis I E，Bourdeau P L.Probabilistic analysis of retaining walls［J］.Computers and Geotechnics，2010，（37）：359－373.

［11］杜永峰，余钰，李慧.重力式挡土墙稳定性的结构体系可靠度分析［J］.岩土工程学报，2008，30（3）：349－353.Du Y F，Yu Y，Li H.Analysis of reliability of structural systems for stability of gravity retaining walls ［J］.Chinese Journal of Geotechnical Engineering，2008，30（3）：349－353.

［12］Yuan X X，Pandey M D.Analysis of approximations for multinomial integration in system reliability computation.Structural Safety，2006（28）：361－377.

［13］Wang Yu，Cao Zijun，Au Siu－Kui.Efficient Monte Carlo Simulation of parameter sensitivity in probabilistic slope stability analysis［J］.Computers and Geotechnics，2010（37）：1015－1022.

［14］Reliability analysis using radial basis function networks and support vector machines［J］.Computer and Geotechnics，2011，38（2）：178－186.

［15］贡金鑫.工程结构可靠度分析方法［M］.大连：大连理工大学出版社，2003.