边坡是一种广泛存在于自然和人工环境中的地质构造，其稳定性对于交通、水利、能源等工程的安全运行至关重要［1-3］。由于边坡岩土力学参数的不确定性和功能函数的高度非线性，边坡可靠度分析面临诸多挑战。边坡可靠度分析的核心挑战在于解决具有隐式的或高度非线性功能函数的边坡失效概率的计算精度较低的问题，当前应对此挑战的常用分析方法主要有蒙特卡罗法［4］、一次二阶矩法［5］、响应面法［6］等。其中，响应面法是一种基于近似模型的高效方法，它通过构建边坡安全系数与随机变量之间的函数关系（即响应面），利用数值积分或极值理论求解失效概率。响应面法具有计算速度快、精度高、适用范围广等优点，已被广泛应用于边坡可靠度分析［7］。
然而，响应面法在拟合高度非线性的功能函数时，仍面临计算精度不足、效率低下的固有局限，为克服上述局限，科研人员已提出多种改进策略与算法。其中，基于代理模型的方法常被用于优化响应面法的计算。苏国韶等［8］提出了一种高斯过程回归动态响应面法（Gaussian process regression dynamic response surface method，GPRDRSM），该方法利用高斯过程回归（Gaussian process regression，GPR）模型构建随机变量与功能函数响应值之间的关系，并进一步通过纳入新增训练样本点，动态地更新响应面函数。Hu等［9］利用GPR模型近似功能函数及其一阶偏导数，并考虑预测值的偏移量，动态地更新训练数据集。谭晓慧等［10］提出了一种基于径向基函数响应面法的岩质边坡可靠度分析方法，使用径向基函数神经网络模型近似岩质边坡的性能函数，并结合蒙特卡罗模拟（Monte Carlo simulation，MCS）估计边坡失效概率。Zhu等［11］提出了一种新颖的高斯过程智能响应面法（GPRSM-base MCS），使用GPR模型近似高度非线性的响应超曲面，并提出用最近邻策略动态地更新训练数据集。
代理模型的拟合精度高度依赖于训练点的选择策略，目前大多数采样策略倾向于优先选取最接近极限状态曲面的点作为训练点更新响应面模型，然而该策略的有效性依赖于训练点分布的均匀性，若训练点过度集中则易导致响应面模型过度拟合，进而影响失效概率的计算精度［12］。为此本文在现有的研究基础上提出一种新的采样函数——自适应搜索引导函数（adaptive search guidance function，ASGF），并结合GPR模型和MCS法，提出一种适用于边坡可靠度分析的自适应高斯过程响应面法（adaptive Gaussian process response surface method，AGPRSM）。
1 GPR模型
GPR模型是用高斯分布概率模型刻画被测数据之间内在联系的一种模型，该方法对数据的内在结构没有特定的限制，可以灵活处理各种非线性和线性关系。对于训练集T=[(X,Y)X∈Rn,Y∈Rn]，其中[X=[x1, x2, ?, xn]T]为[n]维的输入变量，[Y=[y1, y2, ?, yn]T]是相应的输出值。输入变量的随机过程状态集合[g(X)=g(x1), g(x2), ? , g(xn)T]服从联合高斯分布，因此[g]属于高斯过程，其概率函数用GP表示：
[gX ~ GPρX,CX, X’] （1）
式（1）中：[ρX]为均值函数，通常设置为0；[C(X, X’)]为输入变量[X]的[n×n]阶协方差函数矩阵。
[cx, x’=exp-i=1nθixi-xi’2] （2）
式（2）中：[cx, x’]为核函数，[θi]为[xi]的重要性程度。
当存在观测噪声[ε]时，即[y=hX+ε]，假设噪声满足高斯分布[ε~N(0,σ2)]，新数据的输出值为[y?=[y1?,y2?,?,yN?]T]，则其联合概率密度为：
[fy*;μ, σ2, θ=]
[12πσ2n2C12exp-12σ2(y?1)TC-1(y?μ1)] （3）
[μ=1TC-1y*1TC-11] （4）
[σ2=y*-1μTC-1y*-1μn] （5）
式（3）~式（5）中，[1]为元素全是1的矩阵。
将式（4）和式（5）代入式（3），通过极大化对数似然函数求得最优超参数[θ]，并利用式（4）和式（5）求得[μ和σ2]。根据[μ、σ2、θ]，在给定测试输入变量[x0]条件下，[y0]的条件分布为[N(μGP, σGP2)]，其均值和方差分别为：
[μGP=μ+rTC-1y?-μ1] （6）
[σGP2=σ21-rTC-1r+1-1TC-1r21TC-1r] （7）
[r= c(x0, x1), c(x0, x2),?, c(x0, xn)T] （8）
2 AGPRSM
2.1 基于GPR模型的蒙特卡洛模拟
假设[Z=hx]是极限状态功能函数，[x=(x1, x2,?, xn)]是一组随机变量，则系统的失效概率[(pf)]可表示为：
[pf=Z≤0hXxdx=]
[?Z≤0hXx1, x2,?, xndx1dx2?dxn] （9）
式（9）中：当[hx≤0]时，系统处于失效状态，当[hx>0]时，系统处于安全状态；[hXx]是随机变量的联合概率密度函数。
鉴于直接解析求解上述多重积分在计算上不可行，工程实践中通常采用MCS求解失效概率，其表达形式如下：
[pf=1NTi=1NTI[hXi]] （10）
式（10）中：[NT]为MCS样本数，[IhXi]为[hx]的指示函数，即当[hx>0]时，[IhXi=1]，否则[IhXi=0]。
但由于极限状态功能函数往往是隐式的或高度非线性的，因此采用GPR模型构建的近似响应面替代真实功能函数值，并结合MCS估计真实的失效概率。
2.2 自适应搜索引导函数
代理模型对响应面的拟合程度决定失效概率的预测精度，由于极限状态曲面附近的样本点对响应面的拟合精度影响最大，因此往往将这类样本点作为模型的训练点［7］。但如果加入的训练点过度集中在某一区域则易诱发近似响应面光滑性显著降低，从而影响失效概率的预测精度［11］，因此训练点的选择需要在逼近极限状态面和分布均匀之间保持平衡。
为了提高响应面的拟合精度，本文提出ASGF，通过引入一个指数项灵活地调节搜索方向，从而有效地平衡对当前已知区域的开发和对未知区域的探索，确保点的分布不会过于集中，其表达形式如下：
[FASGFx=exp-μhxT?Rimprx?Pimprx] （11）
[T=gσhx] （12）
式（11）和式（12）中：[T]是控制ASGF下降速度的参数；[μhx和σhx]是GPR模型对样本点[x]处的预测均值和标准差；[g]是真实极限状态功能函数的评估总次数。[Rimprx]是预测值[μhx]相较于[hmin]的改进量［13］，若[μhx]大于或等于[hmin]，则改进量为0；[Pimprx]是预测值[μhx]相较于[hmin]改进的可能性［13］。
[Rimprx=maxhmin-μhx,0hmin] （13）
[Pimprx=Φhmin-μhxσhx] （14）
式（13）和式（14）中：[hmin]是当前训练集中最接近极限状态曲面的函数值，当预测值为正时，[hmin]取其当前训练集中最小值，否则[hmin]取其绝对值的最小值；[Φ?]是标准正态分布的累积分布函数。
ASGF通过控制表达式中指数项函数值的大小动态调整新增训练点的位置。当训练点趋近于极限状态面时，指数项函数值趋近于1，此时ASGF倾向于开发已知区域；反之，函数值趋近于0，ASGF倾向于探索未知区域。
此外，参数[T]的值随评估总次数g的增加而变大，并随预测标准差[σhx]的减小而变大，这表明：当对真实极限状态功能函数有更多的信息或者模型更加精确时，ASGF倾向于选择改进程度较大的点加入训练集中，此时算法倾向于开发已知的区域；反之，当参数[T]的值减小时，指数项接近0，此时算法倾向于探索未知的区域。
2.3 实现步骤
AGPRSM算法进行可靠度分析的流程如图1所示，具体实现步骤如下。
步骤1：使用拉丁超立方抽样生成[5n]个（[n]是维度）空间上均匀分布的点，并根据真实极限状态功能函数计算出函数值，以此建立初始训练数据集[D={(Xi, h(Xi))|i=1, ?, 5n}]。
步骤2：基于训练数据集[D]，建立随机变量[Xi]与极限状态功能函数值[hXi]之间的非线性映射关系，由此建立极限状态功能函数的GPR模型。
步骤3：根据随机变量的分布特征，随机抽取[NT]个样本点，根据GPR模型预测得到近似的功能函数值[h(Xi)]。最后结合MCS模拟估计出真实的失效概率的值，为确保计算的精度，设置[NT]为[108]。
步骤4：根据收敛条件[|pf(g)-pf(g-1)|<ε]判断失效概率估计值是否收敛，若收敛则停止迭代，输出结果[pf]，否则继续步骤5。
步骤5：使用差分进化算法并结合ASGF进行搜索，将ASGF值最大的样本点视为最优样本点，并加入训练集[D]中，然后回到步骤2。
上述算法流程可由MATLAB R2020b软件实现。
<G:\武汉工程大学\2024\第3期\邝泽民-1.tif>[开始][拉丁超立方抽样生成
初始训练集D][用极限状态功能
函数计算真实值][用训练集D建立GPR模型][用GPR模型替代极
限状态功能函数
并计算失效概率][使用差分进化
算法优化][选择ASGF值最大的点加入训练集D][否][是][结束][输出失效概率[pf]][[|pf(g)-pf(g-1)|<ε]]
图1 AGPRSM算法流程图
Fig. 1 Flow chart of AGPRSM algorithm
3 数值实验
为了验证本文所提出的方法的有效性，选取2个案例进行实验。案例1是一个高度非线性二元多项式，其随机变量服从标准正态分布且独立。案例2是一个典型的边坡工程案例，其功能函数具有高度非线性且随机变量具有多种分布情况。
3.1 案例1
有一非线性功能函数［14］为
[h1x=exp0.4x1+7-exp0.3x2+5-200]
（15）
随机变量[x1]和[x2]均服从标准正态分布。
对真实极限状态函数进行[108]次MCS模拟，并将结果作为本案例的研究基准（[pf]=3.62×[10-3]）。运用所提出的方法估算失效概率。首先使用拉丁超立方抽样生成10个训练点，每次迭代添加1个新的样本点到训练数据集中，经过第7次迭代后GPR预测的失效概率达到收敛条件[(ε=0.01×10-3)]，此时[pf]=3.64×[10-3]，计算误差仅为0.5%。
图2显示了GPR估计的失效概率的迭代过程；图3显示了AGPRSM的搜索过程，新加入的训练点均位于极限状态函数及其领域，且整体空间分布均匀性良好，图中的序号表示新加入训练点的次序。
3.2 案例2
选取香港秀茂坪边坡［15］为研究对象（图4），假设其坡高[H]= 60 m，坡角[θ]=50°，滑面倾角[ψ]=35°，将黏聚力[c]、内摩擦角[φ]、张拉裂缝深度[z]、张裂隙中充水深度系数[iw]、水平地震加速度系数[α]视为随机变量，变量分布如表1所示，基于拉依达准则，正态分布变量的置信区间取为均值±3倍标准差。此外，假定[c]和φ之间以及[z]和[iw]之间存在负相关，相关系数为[ρc,φ=ρz,iw][=-0.5]。假设边坡的极限状态功能函数为[h2(c, ψ, iw, z, α)=FS(c, ψ, iw, z, α)-1]，则其安全系数[FS]的计算如下：
[FS=cA+[Wcosψ-αsinψ-U-Vsinψ]tanφWsinψ+αcosψ+Vcosψ] （16）
[A=H-zsinψ] （17）
[W=0.5γH21-zH2cotψ-cotθ] （18）
[U=0.5γwzwA] （19）
[V=0.5γwzw2] （20）
[zw=iw?z] （21）
其中：[γ]=26 kN/m3，为岩体重度；[γw]=10 kN/m3，为水的重度。
采用本文方法求解边坡失效概率，如图5所示，经过14次迭代后失效概率达到收敛条件（[ε=0.01×10-2]）。计算结果如表2所示，本文方法估计的失效概率相对误差为0.44%，而函数调用仅有39次，与GPRSM-base MCS方法［11］相比，本文方法的计算结果精度更优，函数调用次数却更少。
表1 随机变量分布参数
Tab. 1 Distribution parameters of random variables
[随机变量 概率分布类型 均值 标准差 相关系数 [c] 正态分布 100 kPa 20 kPa -0.5 [φ] 正态分布 35[°] 5[°] [z] 正态分布 14 m 3 m -0.5 [iw] 指数分布 0.5 截断到[0，1]区间 [α] 指数分布 0.08 截断到[0，0.016]区间 ]
<G:\武汉工程大学\2024\第3期\邝泽民-5.tif>[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
迭代次数 / 次][2.8
2.7
2.6
2.5
2.4
2.3
2.2
2.1
][GPR模型预测的失效概率 / %][（15，2.274）]
图5 案例2中GPR预测的失效概率迭代过程
Fig. 5 Iterative process of predicting failure probabilities
using GPR in Case 2
表2 案例2的可靠度分析结果
Tab. 2 Results of reliability analysis for case 2
[方法 失效
概率 / % 相对
误差 / % 函数调用
次数 / 次 MCS 2.26 - 108 GPRSM-base MCS 2.24 0.88 40 本文方法 2.27 0.44 39 ]
4 结 论
本文提出了AGPRSM，该方法采用GPR模型逼近极限状态函数，并根据自适应采样函数动态更新响应面模型，有效地提高了失效概率的估算精度。
（1）本文通过优化GPRSM-base MCS方法中的采样函数，提出一种用于搜索引导的函数ASGF，该函数通过充分考虑模型的预测信息以及真实极限状态功能函数的评价次数，动态调整搜索策略，实现样本点逼近极限状态曲面且分布均匀，有效地提高了近似响应面模型的拟合精度。
（2）将AGPRSM应用于具有独立标准正态分布随机变量的高度非线性功能函数和具有多种分布类型随机变量的高度非线性功能函数的2种案例。相较于GPRSM-base MCS方法，AGPRSM在求解可靠度问题时函数调用次数更少且计算结果精度更高。