0 引言

在实际的边坡工程中，岩土体的基本参数具有很强的复杂性和变异性，试验虽能获取其物理力学参数，但可获得的试验样本数量有限，不足以表明其分布类型和特征，而获取大批量的样本数据却因为环境或成本等因素的限制变得困难。此时，难以再用例如矩法等常规的可靠度方法来进行概率可靠性分析。尽管试验获得的原始样本数据较少，但它还是能够确定参数的一个大致区间，即表现出“未知而有界”特征。为解决此类问题，非概率可靠性的概念^[1]便应运而生，最早是在应用力学领域由Ben-Haim和Elishakoff提出。1994年，稳健可靠性理论^[2]被提出，作为稳健可靠性的数学基础，凸集合的概念被引入到非概率可靠性模型中，认为当已知的样本信息匮乏时，采用集合模型来描述不确定性的波动范围较合适。郭书祥等^[3]则将用来表述不确定性参数的变量区间化，从而定义了一种新的非概率可靠性指标。王晓军等^[4]用集合表示结构的不确定性信息，并将结构安全域与区间变量域体积之比用来度量结构的非概率可靠性。乔心州等^[5]在提出了一种基于椭球凸集的非概率可靠性模型的同时，还证明了概率模型与非概率模型的相容性。周凌等^[6]则针对已有的超椭球凸集非概率可靠性指标的片面性，提出一种综合指标并给出求解算法。

除了可靠度指标的求解途径外，边坡功能函数的建立也是边坡结构稳定可靠度研究的重点。众多不确定性因素导致边坡的极限状态方程一般都颇为复杂，呈隐含形式且非线性程度较高。因此，将边坡功能函数显式化后再进行相应的可靠度计算分析是本研究旨在解决的问题之一。运用代理模型不失为一种方法，响应面模型^[7]、支持向量机模型^[8]、神经网络模型^[9]等都是目前应用比较广泛的代理模型。穆雪峰等^[10]曾指出地质统计学中的Kriging模型对未知信息的预测能力要优于一些常用的代理模型。

本研究将按照如下思路展开研究：首先，由有限的样本数据确定一个大致区间，在参数的区间描述下构建边坡超椭球凸集模型；然后，将拉丁超立方试验获得的区间内样本信息代入边坡模型，得其响应值后利用Kriging代理模型实现边坡功能函数的显式化。最后，通过换算得到在标准正态空间内的非概率指标η。当η>1和0 < η < 1时，边坡稳定性评价分别采用非概率指标和概率指标。本研究将Kriging代理模型与已有的超椭球凸集模型非概率可靠性理论结合起来，为求解样本模糊特征下边坡可靠度提供一种方法。

1 边坡超椭球凸集模型的构建

目前在工程领域中，描述不确定性参数的凸集模型有不少，大体分为能量有界模型、包容有界模型、斜率有界模型和有界Fourier模型^[11]这4类。有界Fourier模型因其在几何上表示为多维超椭球也被称为椭球凸集模型。描述不确定性参数的凸集模型中，区间凸集模型和椭球凸集模型这两种凸集模型被应用得最为广泛，且区间凸集模型适用于各个参数互不相关的情况，而当参数具有相关性时则适用椭球凸集模型。因此本研究将采用超椭球模型来界定边坡参数的不确定量。

将结构的随机参数和功能函数分别设为X = (X₁, X₂, …, X_n)^T, Z=g_X (X)。根据有限的样本数据取X_i∈ [X_imin, X_imax]，则随机参数X_i的变差δ_i为：

(1)

式中，X为随机参数X_i在区间[X_imin, X_imax]分布的平均值。将随机参数X_i的区间波动程度定义为ΔX_i:

(2)

故有:

(3)

上式区间凸模型在几何空间中表示一超长方体，该长方体的外接等效单位椭球用式子表示如下:

(4)

式中θe_i为椭球半轴。

并且此时超长方体的顶点在球壳上，则有:

(5)

由此看来，如何确定式(3)的最小外接椭球问题就转换成了一个最值问题，即在已知式(5)的情况下求外接椭球的体积最小值。可利用拉格朗日乘数法求解，则半径为1的超椭球凸集模型^[12]表示为:

(6)

式中将W视为一各元素均大于0的对角矩阵，即意味着椭球各主轴分别平行于坐标轴，则有:

(7)

引进向量u_i:

(8)

联合式(8)和式(6)，原定超椭球凸集模型就转化成:

(9)

而E_c为空间U中的一个单位超球集合，由式(8)可得:

(10)

具体求解时，要实现原始空间向标准向量空间的转换^[12]，只需将式(10)代入结构功能函数即可。

2 基于凸集模型的非概率综合指标可靠度求解 2.1 边坡功能函数

在边坡稳定性分析中, 边坡的功能函数通常被定义为g (X)=g (X₁, X₂, …, X_n)。一般来说，当g (X)>0时, 边坡处于稳定状态；当g (X) < 0时, 边坡处于破坏状态；当g (X)=0时, 边坡处于极限平衡的临界状态。具体计算时, 通常利用滑动力S与抗滑力T之间的关系来表示边坡的极限状态方程^[13], 即

(11)

此时传统的边坡稳定性分析方法可能需要进行一些变换来适应这种表示形式。为此，以确定性分析中的安全系数计算方法为基础，一些专家学者重新定义了边坡稳定极限状态方程, 即

(12)

式中安全系数F_s是关于随机变量X的函数, X = (X₁, X₂, …, X_n)。

而确定性分析中安全系数的计算普遍采用极限平衡法，极限平衡法因其简便适用至今已形成多种计算方法，通常有瑞典圆弧法、Bishop法、Janbu法、萨尔玛法、摩根斯坦-普瑞斯法^[14]等。本研究基于极限平衡法中的简化Bishop法^[15]来建立边坡稳定功能函数。

简化Bishop法安全系数F_s计算公式如下：

(13)

边坡功能函数即为：

(14)

(15)

式中，W_i为计算参数；b_i为第i个条块重力、宽度、网弧底面倾角；u_i，N_i，T_i为第i个条块孔隙水压力、法向作用力、圆弧底面剪力；c_i为黏聚力；φ_i为滑面内摩擦角，见图 1。

显而易见，式(14)和(15)中都含有F_s和m_αi，同时F_s又是关于c_i，φ_i，W_i的函数，因此式(14)边坡的功能函数Z是一个高度非线性的隐式函数。需寻求与要求偏导的传统可靠度方法(如矩法)不同的方法以解决边坡的可靠度分析问题。

2.2 Kriging代理模型

Kriging方法通过部分已知信息去模拟某点的未知信息，是一种更具统计性的半参数化插值方法，可以在某点一定范围内求其线性最优无偏估计值。Kriging模型一般由回归部分加随机过程组成，其关系表示成如下形式^[16]:

(16)

式中，ω为回归系数；f (X)为回归模型，一般是关于变量X的多项式函数；Z (X)为一随机过程。

给定容量为m的训练样本S= [X⁽¹⁾, X⁽²⁾, …, X^(m)]与其真实响应Y= [Y⁽¹⁾, Y⁽²⁾, …, Y^(m)], 则可预测X_new的值为任意点。

(17)

式中，R为一对称矩阵，由R (λ: S)构成，大小为m×m，对角线上的元素为1；λ为相关参数；r (X_new)为待测点和训练样本间的相关向量；F是m维向量，由m个样本点处的回归模型组成；f (X_new)为回归多项式，视实际情况确定，一般可采用多项式的形式(2阶或以下)。

在正态随机过程的假设下，构造最优Kriging模型需要通过求解未知量λ，λ则通过求解式(19)所示优化问题获得。于是通过最大似然估计得:

(18)

(19)

根据上述公式推导，将式(17)作为边坡功能函数的近似显示表达式，将边坡稳定性的影响因素设为随机变量形式X = (X₁, X₂, …, X_n)，则边坡功能函数Z的近似表达式为：

(20)

首先将边坡参数训练样本S= [X⁽¹⁾, X⁽²⁾, …, X^(m)]代入简化Bishop模型中获得样本的真实响应值Y= [Y⁽¹⁾, Y⁽²⁾, …, Y^(m)]，然后由式(17)~式(19)得到功能函数近似表达式(20)，求得各待定系数后取回归多项式f (x)为0阶，得边坡功能函数的Kriging代理模型:

(21)

2.3 非概率综合指标可靠度求解原理

文献[6]所定义的一种多个超椭球凸集合描述时的可靠性综合指标k，其表达式为：

(22)

式中，R_set为非概率可靠度；η为非概率可靠性指标。

图 2为二维空间中椭球凸集模型与极限状态曲面的干涉关系，如图所示，不确定变量所对应的凸域由单位圆表示。将标准向量空间的坐标原点到极限状态曲面的最短距离定义为η。功能函数经变换由原始空间转化至标准向量空间U中，极限状态曲线将U划分为两部分，一个g (u) < 0的区域，一个g (u)>0的区域，分别对应失效域和可靠域。由此可得当η>1时，结构总是稳定的，此时可用非概率可靠性指标η，即坐标原点到极限状态曲面的最短距离来评估结构可靠度；而当0 < η < 1时，η指标便不再适用，因此对式(22)中的非概率可靠度R_set进行求解，以此作为可靠度评估指标，将功能函数变换到标准正态空间后再用蒙特卡洛法利用式(23)求解失效概率和可靠度指标^[17]。

(23)

式中，Z为边坡功能函数；N_MC为抽样总次数，N_Z≤0为功能函数Z≤0的样本数量。

综上，边坡凸集模型非概率综合指标可靠度求解方法的操作流程为:

(1) 根据已有的样本参数信息确定样本区间范围，设影响边坡稳定性的随机变量为X = (X₁, X₂, …, X_n)。

(2) 构建边坡超椭球凸集模型, 得到实现边坡功能函数由原始空间转换至标准正态空间的桥梁式(10)。

(3) 对样本区间采用拉丁超立方试验来构造样本点，代入简化Bishop模型式(14)获得其响应值。

(4) 通过Kriging代理模型使边坡功能函数显式化。

(5) 引入基于超椭球凸集的非概率可靠性综合指标η，当η>1时，通过迭代计算得到坐标原点到极限状态曲面的最短距离，以此作为评估标准；当0 < η < 1时，将式(10)代入Kriging代理模型，运用蒙特卡洛法根据式(23)求解失效概率P_f。

3 实例分析

一双层边坡^[18]，坡比为1∶2，高15.24 m。其土层参数信息见表 1。

其中黏聚力c₁，c₂和内摩擦角φ₂未知其具体数值与服从的分布类型，根据试验经验，把c₁，c₂，φ₂作为随机参数，表示为X = [X₁, X₂, X₃]。

根据上述的可靠度求解步骤，所建立的半径为1的超椭球凸集模型为：

(24)

计算得X，ΔX_i:

求得加权矩阵W :

(25)

对区间采用拉丁超立方试验构造样本点，代入边坡功能函数式(14)得响应值后，由Kriging代理模型拟合获得边坡功能函数。

然后由式(10)，进行标准化变换的矩阵：

(26)

代入Kriging模型将功能函数转化至标准正态空间。最后通过判断，该边坡有失稳的可能性，为0 < η < 1的情况，再运用蒙特卡洛法求得失效概率P_f和可靠度指标β，结果如表 3所示。

由表可知，当取30组或40组样本数据时求得的失效概率和可靠度指标相对误差较小，且40组数据抽样结果：P_f=1.47%，β=2.177 5，与文献[18]中已知随机变量的具体数值和已知分布类型的情况下求得的结果：P_f=1.46%，β=2.180 0非常接近。由此证明本研究提出的模糊样本特征条件下的基于超椭球凸集模型和Kriging代理模型的非概率边坡稳定可靠度求解方法是可行的，并且精确度很高，较于蒙特卡洛法其效率也明显提高。

4 工程案例分析

湖南省某高速公路是国家规划“七纵九横”的第4纵南益高速公路的重要一段，全长86.724 km，从G56杭瑞高速程家山枢纽开始到S7101益阳绕城高速迎丰桥枢纽止。其中一段路堑边坡高度约17.5 m，坡比为1∶3，层状构造，其岩土层主要由(1)杂填土(2)粉质黏土(3)强风化泥质砂岩组成。

实地勘探得到其有限的物理力学参数样本数据如表 4所示。

表中黏聚力c₁，c₂，c₃，内摩擦角φ₁，φ₂，φ₃和重度γ₁，γ₂，γ₃都只知其大致范围，根据试验经验，把c₁，c₂，c₃，φ₁，φ₂，φ₃，γ₁，γ₂，γ₃都作为随机变量，表示为X = [X₁, X₂, …, X₉]。

根据2.3节可靠度求解步骤, 首先构建半径为1的超椭球凸集模型：

(27)

其中X，ΔX_i计算结果如表 5所示。

从而求得加权矩阵

(28)

然后通过区间拉丁超立方试验得到40组样本，代入边坡功能函数式(14)得到响应值，如表 6所示，将表 6中数据代入Kriging代理模型计算得功能函数。

再由式(10)，进行标准化变换的矩阵：

(29)

代入Kriging模型将功能函数转化至标准正态空间。

最后判断该边坡属于0 < η < 1的情况，运用蒙特卡洛法求得失效概率和可靠度指标P_f=4.0×10^-5，β=3.944 4，参考《公路工程结构可靠性设计统一标准》 (JTG 2120—2020)中高速公路路基目标可靠度指标：3~2.5，结果表明该边坡失稳概率非常低，边坡比较稳定，同时对照表 6中40组数据，由简化Bishop法计算得到的40组安全系数平均值F_S=2.033，远大于边坡安全系数临界值1，进一步验证了该边坡的稳定性。

5 结论

研究表明，本研究主要取得以下几个方面的成果：

(1) 由于边坡岩土体参数信息的不完备性，利用其有限的原始样本数据构建了边坡超椭球凸集模型，为边坡功能函数的空间转换提供了便利。

(2) 利用LHS抽样和具备超强学习能力的Kriging代理模型实现边坡功能函数显式化，然后将其变换至标准正态空间，求解失效概率和可靠度指标。

(3) 求解过程引入基于超椭球凸集的非概率可靠性综合指标η，当η>1时，通过坐标原点到极限状态曲面的最短距离来衡量边坡稳定性；当0 < η < 1时，运用蒙特卡洛法即可求得边坡失效概率P_f与可靠度指标β。

(4) 算例分析表明本研究方法可行且选取40组样本数据时精确度较高；对照《公路工程结构可靠性设计统一标准》 (JTG 2120—2020)，采用本研究方法分析工程案例得出的结论与运用简化Bishop法的计算结果一致。算例分析与工程案例分析充分证明了本研究方法的可操作性、计算精度及工程实用性。