0.   引言

边坡失稳会造成道路掩埋、交通中断、建筑物坍塌等,给人民生命财产安全造成重大的损失^[1]。受地形地质条件的影响,在水利水电、高速公路、铁路、矿山等基础工程建设中,常常需要进行边坡工程开挖,从而形成了大量不同规模的工程边坡。边坡稳定性分析及设计是岩土工程中的经典问题,与中国国民经济基础建设发展密切相关。传统边坡设计主要是基于单一安全系数的确定性设计方法,该方法通过试算确定满足边坡规范允许安全系数要求的边坡设计方案。如熊爽等^[2]采用赤平投影图解法和极限平衡法,通过确定性分析对矿山岩质边坡坡角进行设计。该方法虽然计算简单、直观,但是不能有效地考虑边坡工程中各种不确定性的影响。此外,工作应力法^[3]、荷载和抗力系数设计方法^[4]以及偏安全系数法^[5]是目前欧洲和北美较为流行的工程设计方法,但由于依然属于确定性分析的范畴,仍未拓展至不确定性边坡设计领域。

边坡工程由于自身条件的复杂性及试验条件的限制,其土体参数往往存在一定的不确定性。即使对于相同土体,土体特性在空间上也呈现明显的不均匀分布,这就要求在可靠度设计中考虑土体参数空间变异性的影响。近年来,考虑边坡工程多种不确定性的可靠度设计方法受到了学者和工程师的青睐。如彭兴等^[6]提出一种基于蒙特卡洛模拟的岩质边坡可靠度设计方法。邓志平等^[7]同时考虑地层变异性和土体参数空间变异性,提出了基于广义耦合马尔可夫链与随机场模型耦合的边坡可靠度全概率设计方法。陈祖煜等^[8]基于分项系数设计方法,通过改变承德机场工程加筋土坡坡高、坡度和其他材料参数,标定边坡安全系数和分项系数取值标准。李典庆等^[9]基于安全系数累积分布函数提出了边坡广义可靠指标相对安全率设计方法。综上,基于概率分析（蒙特卡洛法、广义可靠指标、分项系数设计方法等）的边坡可靠度设计方法能够更为客观、真实地反映边坡安全性能。然而,因较少解释土体参数空间变异性作用导致边坡设计结果与工程实际存在一定的偏差,同时需要进行多次可靠度试算来满足目标失效概率要求,设计工作量较大。

边坡可靠度设计通常已知目标失效概率设计坡角、坡高等,但是如何在少量实际试验数据条件下考虑土体参数空间变异性的影响,高效地进行边坡可靠度设计仍是一个关键技术难题。幸运的是,在可靠度分析基础上发展起来的逆可靠度方法为解决这一难题提供了一个有效的工具。可以根据边坡工程实际及重要级别,在已知目标失效概率的条件下寻求最优边坡参数（如坡角、坡高等）。如苏永华等^[10]将响应面方法与基于一次二阶矩的逆可靠度算法有机结合,提出了基于风险控制的边坡稳定逆可靠度设计方法。Papaioannou等^[11]结合一阶逆可靠度方法和有限元强度折减法设计边坡最优坡角,并探讨了土体参数空间变异性对坡角设计的影响。此外,Ji等^[12]基于（Hasofer-Lind-Rackwitz-Fiessler,HLRF）算法发展了简化正一阶和逆一阶可靠度迭代算法可用于边坡可靠度设计。但是,少量的实际数据条件下空间变异土坡可靠度设计问题仍研究较少。

为此,本文提出了适合任意坡角的空间变异土坡模型重构方法,基于一阶逆可靠度方法和极限平衡分析,综合考虑3种不确定性（试验数据测量不确定性、模型转换不确定性和模型修正系数不确定性）的影响,从可靠度逆分析角度对少量实际试验数据条件下空间变异砂土边坡进行设计,从而获得切合工程实际的边坡设计方案。

1.   空间变异土坡可靠度设计

1.1   空间变异土坡稳定性分析

一阶可靠度方法^[12-13]（first-order reliability method, FORM）已广泛应用于考虑岩土工程各种来源不确定性的边坡失效概率计算中。为了表征土体参数的空间变异性,通常将参数随机场离散为一组随机向量X = (X₁, X₂, …, X_N)^T,其中N为随机场单元数目。目前随机场离散方法很多,如Karhunen-Loève展开方法、中点法和局部平均法。本文采用基于乔列斯基分解的中点法离散相关非高斯参数随机场,该方法简便且易编程实现^[14]。对于已知坡角为α的土坡而言,基于极限平衡方法（LEM）的边坡可靠度分析极限状态函数可表示为

式中：x = (x₁, x₂,…, x_N)^T为原始空间随机向量X的实现值；w为模型修正系数,用来表征采用LEM计算F_S存在的模型误差。在FORM方法中,需要将原始空间随机向量x转换为独立标准正态随机向量u,x = T(u),这可以通过Nataf逆变换^[15]得到,具体可通过以下两个公式实现：

式中 y = (y₁, y₂,…, y_N)^T为相关标准正态随机向量；$F_{X_i}{}^{-1}(\cdot )$为随机变量X_i的边缘累积分布函数的逆函数；$\Phi (\cdot )$为标准正态变量累积分布函数；L为下三角矩阵,可通过对自相关系数矩阵R进行乔列斯基分解得到

根据式（2）～（4）,可将式（1）的极限状态函数表示为

在基于FORM的边坡可靠度分析中,关键一步是求解独立标准正态空间中极限状态面上概率密度最大的点^[16],即验算点${\text{u}}^{\ast }$,可通过求解以下有约束的二阶非线性方程得到

Ji等^[12]利用HLRF算法将独立标准正态U空间中的极限状态函数表示为点${\text{u}}_k$处的泰勒级数展开：

进而u的一阶近似可表示为

式中 ${\text{u}}_k$和${\alpha }_k$均为独立标准正态U空间中u和$\alpha$的第k次迭代点；$g({\text{u}}_k,{\alpha }_k)$和${\nabla }_{u}g({\text{u}}_k,{\alpha }_k)$分别为迭代点${\text{u}}_k$处的功能函数值及梯度向量。为获得验算点u^*,需要对式（8）进行迭代求解。为提高迭代计算效率,定义如下搜索向量：

另外,为了稳定迭代计算过程,引入一个步长${\lambda }_k$,进而迭代计算表达式^[16]可表示为

式中,${\lambda }_k$和${\text{d}}_k$分别为U空间中的搜索步长和方向。一旦获得验算点${\text{u}}^{\ast }$,便可根据验算点与原点之间的距离计算边坡可靠度指标：

式中,e为验算点处的单位法向量,$e=-\nabla \text{g}(u^{\ast },\alpha )/$ $\Vert \nabla \text{g}(u^{\ast },\alpha )\Vert$,其中$\nabla \text{g}(u^{\ast },\alpha )$是极限状态函数在验算点处的梯度向量。进而可以求解得到边坡失效概率为

1.2   一阶逆可靠度坡角设计方法

可靠度逆分析的目的是在考虑不确定性因素条件下反算满足目标失效概率的设计参数,是对经典可靠度分析的拓展。逆可靠度方法需要定义一个约束性条件：$p_f$=$p_{f,t}$,其中$p_{f,t}$为目标失效概率。本文采用在FORM基础上发展起来的一阶逆可靠度方法设计一个满足目标失效概率的坡角α,运用验算点法对极限状态函数进行逆分析,基于目标失效概率确定边坡设计参数,无需进行额外的数值计算。本文逆可靠度问题定义为搜索最优坡角α,使得边坡失效概率等于目标失效概率$p_{f,t}$,并且极限状态函数值等于0,即

Zhang等^[16]和Ji等^[12]给出了U空间逆可靠度分析中的验算点u和设计参数$\alpha$的迭代计算表达式如下：

若采用式（14）和（15）进行迭代计算,每次迭代计算都需要将参数由U空间转换至原始X空间来计算边坡安全系数,而后又需要进行逆变换推求验算点,且偏导$\partial g({\text{u}}_k,{\alpha }_k)/\partial \alpha$的求解也更加困难,这大大增加了编程的难度。为便于边坡安全系数计算,Ji等^[17]对原有的HLRF算法^[18]进行了改进,将式（14）,（15）拓展到了原始相关非正态X空间。式（14）,（15）中U空间的梯度向量可以表示为

式中,${\nabla }_{x}f({\text{x}}_k,{\alpha }_k)$为X空间的梯度向量。进而可以通过HLRF-x算法计算X空间逆可靠度分析中的验算点x^*和设计参数$\alpha$如下：

式中,${\text{μ}}_k^N$为第k次迭代点处的等效正态均值向量；T_k为在第k次迭代点计算的修正自相关系数矩阵,其第i行和第j列的元素$T_{k,ij}={\sigma }_{k,i}^N{R}_{ij}{\sigma }_{k,j}^N$,其中$R_{ij}$为自相关系数矩阵R中的元素；${\sigma }_{k,i}^N$和${\sigma }_{k,j}^N$为随机变量等效正态标准差对角矩阵$[{\sigma }_k^N]$中的元素。根据式（18）便可设计出满足目标失效概率$p_{f,t}$的最优坡角α。由于边坡极限状态函数通常是一个隐式非线性函数,梯度向量${\nabla }_{u}g({\text{u}}_k,{\alpha }_k)$和偏导$\partial f({\text{x}}_k,{\alpha }_k)/\partial \alpha$没有解析解,故本文采用有限差分法进行计算。此外,在坡角设计过程中,土体参数随机场值${\text{x}}_k$在每次迭代计算时通过式（17）进行更新,更新后的${\text{x}}_{k+1}$将逐渐靠近于验算点${\text{x}}^{\ast }$,坡角α值通过式（18）进行更新。需要指出的是,坡角α变化会改变边坡模型的几何形状,相应的边坡破坏模式也会发生明显的变化。因此,针对每个不同的坡角均需要重新构建空间变异土坡模型。

2.   适合任意坡角的空间变异土坡模型重构方法

当考虑土体参数空间变异性时,土坡模型的不断重构会给边坡可靠度设计带来很大的难度。为解决这一难题,本文提出了适合任意坡角的空间变异土坡模型重构方法。以仅考虑土体参数垂直方向一维空间变异性的土坡为例,该方法主要计算步骤简要介绍如下：

（1）假设该边坡在垂直方向上划分为N个随机场单元,每个随机场单元与斜坡的交点为m_i（i = 1, 2,…, N）。已知坡高为H和坡角为α,可以得到N个交点的坐标分别为m₁（H/Ntanα,H/N）,m₂（2H/Ntanα,2H/N）,…,m_N（H/tanα,H）。图1中取N = 12来说明参数随机场划分情况。

图  1  边坡模型重构示意图

Figure  1.  Schematic diagram of slope model reconstruction

下载: 全尺寸图片 幻灯片

（2）找到某一代表性滑动面S₁与斜坡和坡顶的两个端点（见图1中的点A和B）。将端点A与B之间的水平距离等分成N_l个土条（见图1中分割实线）,通过垂直投影将对应的滑动面S₁划分成N_l段。一般来说,滑动面分割的数目越多,则对边坡安全系数的估计就越准确,当划分数目为100时可满足计算精度要求^[19]。图1中仅取N_l = 10来说明滑动面S₁的划分情况。

（3）为了保证每段滑动面只落在一个随机场单元内,需要找到随机场单元与滑动面S₁的交点（如图1中的点C和D）,根据这些点的垂直投影,进一步将滑动面划分为额外的N_s段（见图1中的分割虚线）。

（4）根据步骤（2）和（3）,可以将滑动面S₁划分为（N_l+N_s）段。由图1可见,每段滑动面只落在一个随机场单元内。提取每段滑动面所在的随机场单元的中心点坐标,产生土体参数（内摩擦角和黏聚力等）随机场实现值,并赋给对应的那段滑动面^[14]。

（5）一旦获得每段滑动面上的土体参数取值,便可采用极限平衡法计算滑动面S₁所对应的边坡安全系数。

（6）如果边坡中存在N_p条潜在滑动面,重复步骤（2）～（5）N_p次,便可获得N_p条潜在滑动面所对应的安全系数,取其中最小的安全系数视作边坡安全系数FS,对应的滑动面即为最危险滑动面（破坏模式）。

（7）在迭代计算过程中,一旦土体参数随机场值${\text{x}}_k$和坡角α_k发生变化,便需要更新步骤（1）中的N个交点的坐标进而获得新的边坡几何模型,接着重复步骤（2）～（6）,计算新坡角α_k对应的边坡安全系数和最危险滑动面。

需要说明的是,参数随机场值${\text{x}}_k$和坡角α_k的更新也将导致偏导$\partial f({\text{x}}_k,{\alpha }_k)/\partial \alpha$、自相关系数矩阵T_k及等效正态均值${\text{μ}}_k^N$等参数发生变化,据此便可借助式（17）和（18）进行下一步迭代计算,从而获得一个稳定的设计坡角α。

综上,基于一阶逆可靠度方法进行空间变异土坡坡角设计的主要步骤如下：

（1）收集静力触探试验（CPT）、标准贯入试验（SPT）等试验数据,包括锥尖阻力$q_c$、标准贯入击数N_SPT等。本文选用国际土力学及基础工程协会提供的一些相关边坡工程现场试验数据,具体详见https://www.issmge.org/news/are-we-overdesigning-a-survey-of-international-practice。

（2）采用相应的模型转换公式获得土体内摩擦角及不排水抗剪强度等参数数据,并据此计算各自的统计特征（包括均值、标准差、概率分布）,在模型转换过程中需要考虑试验数据的测量确定性及模型转换不确定性的影响。

（3）根据土体参数统计特征和目标边坡失效概率$p_{f,t}$,考虑模型修正系数w的不确定性,采用一阶逆可靠度方法设计满足目标边坡失效概率$p_{f,t}$的最优坡角。

3.   算例分析

下面以一个代表性的砂土边坡为例,基于一阶逆可靠度方法进行坡角设计。边坡几何形状和基本条件如图2所示。对于砂土边坡而言,其破裂面大多为浅层滑动面^[20-21],故采用图3所示的砂土边坡模型,砂土边坡的黏聚力为0,内摩擦角为${\phi }^{\prime }$,砂土重度取20 kN/m³。边坡高为6 m,坡角为α。以坡趾为原点建立坐标轴,其坡顶坐标为（6/tanα m,6 m）。设计满足3种目标失效概率$p_{f,t}$（1×10^-4,1×10^-3和1×10^-2）的坡角。

图  2  砂土边坡示意图

Figure  2.  Schematic diagram of a sandy slope

下载: 全尺寸图片 幻灯片

图  3  砂土边坡随机场离散模型

Figure  3.  Random field discretization model for sandy slope

下载: 全尺寸图片 幻灯片

为了表征土体参数固有的空间变异性,将砂土边坡在垂直方向上分为12个随机场单元,每个随机场单元厚度为0.5 m。采用指数型自相关函数模拟土体参数空间自相关性,计算表达式为^[21]

式中 ${\delta }_z$为垂直方向上的波动范围,z_i和z_j分别为第i个和第j个随机场单元的中心点坐标。在此基础上,采用基于乔列斯基分解的中点法生成土体参数随机场实现值。

3.1   砂土边坡参数统计特征

首先,选用现场静力触探CPT试验获得的5组锥尖阻力$q_c$数据,如图4所示。

图  4  砂土边坡的5组砂土CPT试验数据

Figure  4.  Five sets of CPT test data for sandy slope

下载: 全尺寸图片 幻灯片

根据文献[22],采用如下模型转换公式可获得砂土有效内摩擦角${\phi }^{\prime }$为

式中 $q_{\text{c}}$为锥尖阻力；${{\sigma }^{\prime }}_{\text{v0}}$为垂直有效应力,${{\sigma }^{\prime }}_{\text{v0}}=$ $\gamma h$,其中$\gamma$和h分别为土体的重度和埋深；$p_{\text{a}}$为标准大气压,取100 kPa；$\epsilon$为模型转换误差,通常将其模拟为均值为0,标准差为一常数的正态随机变量。由此可得到5组随埋深变化的砂土内摩擦角数据,如图5所示。

图  5  由CPT数据估计的砂土内摩擦角随埋深变化曲线

Figure  5.  Variation of effective friction angles for sandy slope estimated from CPT data with depth

下载: 全尺寸图片 幻灯片

Akkaya等^[23]采用指数型自相关模型对德克萨斯农工大学国家岩土工程试验和现场数据进行分析,发现砂土参数的垂直波动范围${\delta }_z$约为1.75 m,故本文采用${\delta }_z$=1.75 m进行坡角设计。同时根据Phoon等^[24]可知,通过CPT数据估算空间平均有效内摩擦角${{\phi }^{\prime }}_{\text{a}}$的变异系数（COV）的计算表达式为

式中,$S{D}_{\epsilon }$为式（20）中模型转换误差εe的标准差,取$S{D}_{\epsilon }$= 2.8°^[24]；$CO{V}_{\text{e}}$表示为锥尖阻力$q_c$测量误差的变异系数,取值范围一般为5%～15%^[24],本文按照文献[25]的做法,取$CO{V}_{\text{e}}$=15%；$CO{V}_{\text{w}}$为锥尖阻力$q_{\text{c}}$的变异系数,可通过对图4的5组CPT试验数据计算获得；当随机场单元尺寸较大,需要考虑局部平均效应,用表征局部平均效应的COV缩减因子$\Gamma (L)$表示,对于指数型自相关函数,$\Gamma (L)$的计算表达式为^[26]

式中,L表示为随机场网格尺寸,L = 0.5 m。由此可通过式（21）计算得到12个随机场单元的$CO{V}_{{{\phi }^{\prime }}_{\text{a}}}$,如表1所示。另外,每个随机场单元的均值等于图5中落在每个随机场区域内的有效内摩擦角数据的样本均值。综上,12个土层的有效摩擦角的均值,见表1,并假设其服从对数正态分布。

表  1  砂土边坡不确定参数统计特征

Table  1.  Statistics of uncertain input parameters for sandy slope

${{\phi }^{\prime }}_i$ /(°)	埋深/m	${\mu }_{{\phi }^{\prime }}$ /(°)	$CO{V}_{\text{w}}$	$CO{V}_{{{\phi }^{\prime }}_{\text{a}}}$	概率分布
${{\phi }^{\prime }}_1$	0.25	37.228	0.391	0.090		对数正态
${{\phi }^{\prime }}_2$	0.75	40.791	0.463	0.086		对数正态
${{\phi }^{\prime }}_3$	1.25	41.405	0.310	0.077		对数正态
${{\phi }^{\prime }}_4$	1.75	40.918	0.335	0.079		对数正态
${{\phi }^{\prime }}_5$	2.25	39.412	0.469	0.090		对数正态
${{\phi }^{\prime }}_6$	2.75	36.287	0.501	0.100		对数正态
${{\phi }^{\prime }}_7$	3.25	37.248	0.472	0.095		对数正态
${{\phi }^{\prime }}_8$	3.75	39.315	0.306	0.081		对数正态
${{\phi }^{\prime }}_9$	4.25	39.727	0.203	0.076		对数正态
${{\phi }^{\prime }}_{10}$	4.75	40.262	0.252	0.077		对数正态
${{\phi }^{\prime }}_{11}$	5.25	40.059	0.309	0.080		对数正态
${{\phi }^{\prime }}_{12}$	5.75	37.787	0.420	0.091		对数正态

下载: 导出CSV 

| 显示表格

3.2   砂土边坡坡角可靠度设计结果

根据工程实际,不考虑地下水位及裂缝对砂土边坡的影响。对于每个坡角α及参数随机场实现值,采用简化毕肖普方法计算边坡安全系数。图6为一个坡角为α = 34°的砂土边坡,其中随机生成了1636条潜在滑动面。基于表1中参数均值进行确定性边坡稳定分析,计算的边坡安全系数为1.204,其最危险滑面如图6所示。另外,图7给出了基于表1中参数均值的砂土边坡安全系数随坡角的变化关系曲线。可以看出,在容许安全系数为1.3的情况下,只要设计坡角小于32°,边坡便不会发生失稳破坏。

图  6  含1636条潜在滑动面的砂土边坡

Figure  6.  Sandy slope with 1636 potential slip surfaces

下载: 全尺寸图片 幻灯片

图  7  砂土边坡安全系数随坡角变化曲线关系

Figure  7.  Variation of factor of safety for sandy slope with slope angle

下载: 全尺寸图片 幻灯片

上面采用的简化毕肖普方法是一种极限平衡分析方法,由于条间力假设导致计算的边坡安全系数与实际边坡安全系数$F_{\text{S}}{}_{\text{A}}$之间可能会存在着一定的偏差^[27-28],为此引入模型修正系数w表征这一偏差,进而边坡安全系数可表示为

式中,${\phi }^{\prime }$为参数随机场一次实现,${\phi }^{\prime }={({{\phi }^{\prime }}_1,{{\phi }^{\prime }}_2,\cdots ,{{\phi }^{\prime }}_{12})}^T$；w为模型修正系数。根据Bahsan等^[29],可将w模拟为均值为1.066,标准差为0.296的对数正态随机变量。由此,可以得到砂土边坡的极限状态函数为

任取初始坡角α = 20°,基于表1的参数统计特征,采用FORM方法^[30]和式（24）的极限状态函数可以算出该状态下的边坡失效概率为5.02×10^-3。对于3个不同的目标失效概率（1×10^-4,1×10^-3和1×10^-2）,采用一阶逆可靠度方法可快速设计出对应的坡角,如图8所示。

图  8  不同目标失效概率下砂土边坡坡角可靠度设计结果

Figure  8.  Reliability-based designed slope angle for sandy slope under three scenarios of target probability of failure

下载: 全尺寸图片 幻灯片

由图8可见,只需经过4,5次迭代计算,得到的可靠度设计坡角变化不大。当目标失效概率分别为1×10^-4,1×10^-3和1×10^-2时,基于CPT试验数据设计的边坡坡角分别为14.13°,17.22°和21.56°。显然,设计的坡角随着目标失效概率的增大而增大。同时,为了验证提出方法的有效性,本文重新建立了坡角分别为14.13°,17.22°和21.56°的3个砂土边坡模型,基于表1的统计特征,采用FORM方法^[30]计算的边坡失效概率分别为9.909×10^-5,1.083×10^-3和1.062×10^-2,计算结果与目标失效概率吻合,说明了本文提出方法的有效性。因此,对于目标失效概率1×10^-4,应保证设计的坡角小于14.13°才能满足边坡稳定性要求。值得注意的是,迭代过程中坡角的变化会导致边坡模型几何形状发生变化,采用直接蒙特卡洛模拟等常规可靠度方法进行坡角试算非常耗时。相反,本文提出方法计算效率高,对于常见的隐式边坡可靠度设计问题尤其有效。

4.   结论

本文提出了适合任意坡角的空间变异土坡模型重构方法,并定量考虑多种不确定性因素,基于一阶逆可靠度方法对少量实际试验数据条件下空间变异土坡坡角进行可靠度设计。并以一个砂土边坡为例对边坡坡角进行可靠度设计。主要结论如下：

（1）提出方法可以考虑试验数据测量不确定性、模型转换等不确定性以及模型修正系数不确定性的影响,较好地解决了参数空间变异性条件下土坡坡角设计问题。相比之下,确定性分析方法因不能定量解释边坡工程多种不确定性因素的作用,获得的边坡设计方案与工程实际相差较大。就本文算例而言,对于目标失效概率1×10^-4,采用提出方法设计的砂土边坡坡角应小于14.13°,而采用确定性分析方法设计的边坡坡角与之相差较大。表明十分必要定量考虑多种不确定性进行边坡可靠度设计,其计算结果可作为边坡确定性设计的有益补充。

（2）提出方法计算效率高,仅需进行4,5次迭代计算便可设计出满足目标失效概率的最优坡角。在迭代计算中,坡角变化势必导致边坡模型几何形状发生改变,相应的边坡破坏模式也会发生明显的变化。尤其当考虑土体参数空间变异性时,边坡模型的不断重构会给可靠度设计带来很大的难度。而本文提出的适合任意坡角的空间变异土坡模型重构方法可有效解决这个难题,并极大地提高了计算效率,可为常见的隐式边坡可靠度设计奠定基础。

（3）由于边坡工程实际试验数据通常较少,本文仅考虑了土体参数垂直方向上一维空间变异性的影响,并且是通过相关文献统计资料确定了参数波动范围取值。如何基于少量的实际试验数据对高维空间变异土坡进行可靠度设计仍需要进一步研究。