0.   引言

降雨诱发滑坡是世界上许多国家和地区面临的主要地质灾害之一，每年均会造成大量人员伤亡和财产损失^[1]。受沉积和后沉积、风化、不同荷载和应力历史条件等因素的影响，土体参数呈现出显著的空间变异性，这对边坡可靠度评估具有重要影响^[2]。客观描述土体参数空间变异性是评估边坡可靠度的重要一步。虽然随机场和地统计学等理论为量化土体参数空间变异性提供了有效的工具，为可靠度研究奠定了良好的基础，但是大多数研究主要是基于先验统计信息对土体参数空间变异性进行量化。土体参数先验信息一般是根据文献资料、工程类比和相关假设等途径统计得到^{[3 − 4]}，未充分利用边坡场地信息。因而，仅基于先验统计信息描述的土体参数空间变异特征会与工程实际存在偏差，进而造成对边坡可靠度评估的错误认识。为准确揭示降雨入渗下边坡失效概率演化规律，首先需要提前利用边坡场地多源信息对空间变异土体参数进行更新。目前从边坡场地不同阶段可收集的多源信息，主要包括降雨数据、室内外试验数据、压力水头和变形等监测数据、观测到的滑动面位置以及过去服役状态等实际观测信息^[5]。利用边坡场地多源信息能够更新空间变异土体参数概率分布，在较大程度上降低对土体参数不确定性的估计。进而，基于参数后验统计信息能够更加准确地评估边坡可靠度。

贝叶斯方法可以将土体参数先验信息与多源场地信息相结合，为土体参数概率反分析提供一种符合工程实践的有效工具^[6]。对于降雨入渗下确定性边坡渗流和稳定性分析，需要提前确定包括水力和抗剪强度参数在内的土体参数值。一旦考虑这些参数的空间变异性，需要同时反分析成千上万个随机变量。而常用的马尔科夫链蒙特卡洛（Markov chain Monte Carlo，MCMC）方法对于这个高维参数更新问题基本上无能为力，因为随着随机变量数量的增加，这类方法的计算复杂性也会急剧增加^[7]。幸运的是，Straub和Papaioannou^[8]提出的基于子集模拟的贝叶斯更新（BUS）方法为解决此类高维问题提供了有效工具。目前该方法已广泛应用于边坡、堤坝等岩土结构可靠度更新^[9]以及软土地层的高效识别^[10]。但是，在利用多源信息对空间变异土体参数进行反分析时，该方法存在计算精度不高、效果不佳等缺点^[7]。主要是因为利用多源信息进行贝叶斯更新时所构建的似然函数是一个高维函数，似然函数值非常小，甚至低于计算机浮点运算精度^[11]，如果对此处理不当便会影响土体参数更新的精度和效率。因此，如何有效处理考虑土体参数空间变异性造成的高维概率反分析和边坡可靠度预测问题仍是一个亟需深入研究的难题。

为有效解决极小似然函数值和高维随机变量的概率反分析问题，本文将改进的BUS方法（mBUS）与自适应条件抽样（aCS）算法^[12]相结合，并构建空间变异土体参数概率反分析和边坡可靠度预测框架。最后以某公路边坡为例说明所构建框架的有效性和合理性。

1.   基于多源信息融合的土体参数反分析

对于一个特定的边坡来说，土体参数往往具有一定的不确定性。在运行阶段，可以从边坡现场获取大量监测数据和现场观测信息。融合这些场地信息不仅可以缩减土体参数的不确定性，还可后续预测得到更准确的模型输出响应（如压力水头、安全系数），从而为工程决策提供依据。根据贝叶斯理论，融合监测数据或现场观测信息的参数后验概率密度函数（PDF）可表示为^[8]：

式中：x为不确定土体参数X的实现值，其中X = (X₁, X₂, ···, X_n)^T，n为随机变量数目；$f'({\boldsymbol{x}})$和$f''({\boldsymbol{x}})$分别是X的先验和后验概率密度函数；$L({\boldsymbol{x}})$为似然函数，表示当随机变量取某一特定实现值时某一场地信息事件发生的概率。构建恰当的似然函数对于空间变异土体参数概率反分析至关重要^[11]。对于间接信息（如压力水头或位移等监测数据）来说，由于模型误差和测量误差，边坡第i个位置处监测数据$ {y_i} $（i = 1, 2, ···, t，t为监测点数目）与模型预测h(x)之间的关系可表示为

式中：$ \varepsilon _i,(i = 1,2,\cdots,t) $,为第i个压力水头监测值和计算值之间的测量误差，假定其在统计上是相互独立的并服从均值μ_ε = 0，标准差σ_ε = 0.05的正态分布^[5]。相应地，基于等量信息建立的似然函数可表达成如下形式：

为了保证计算的数值稳定性，通常将式（3）似然函数进行如下对数化处理：

此外，边坡在过去降雨事件下的服役状态（安全系数 FS>1.0 或FS≤1.0），也可为土体参数反分析和边坡安全评估提供有价值的信息^[13]。利用观测信息（安全系数y>1.0）可建立基于不等量信息的似然函数如下：

式中：Φ（·）为标准正态变量的累积分布函数，FS(x)为通过极限平衡等方法或代理模型估算的边坡安全系数，ζ为真实安全系数与估算安全系数之间存在的模型偏差，一般假定ζ服从均值为$ {\mu _\zeta } $，标准差为$ {\sigma _\zeta } $的正态分布^[13]。如果使用监测到的压力水头作为监测信息来构建$ {L_1}(\boldsymbol{x}) $，只能反分析水力参数。如果使用边坡服役状态作为现场观测信息来构建$ {L_2}(\boldsymbol{x}) $，可以同时反分析水力参数和抗剪强度参数。

求解式（1）便可得到土体参数后验概率密度函数。Straub和Papaioannou^[8]提出的BUS方法可以有效地数值求解式（1）得到高维$f''(\boldsymbol{x})$的近似解，通过引入一个与似然函数L(x)有关的常数c，构建一个等效的场地信息事件失效区域$\varOmega $如式（6）:

式中：u为服从[0, 1]的均匀分布，c为似然函数乘子，对于所有的x，均需要满足：

据此，基于等量或不等量信息的贝叶斯更新问题就转换为一个结构可靠度问题，之后便可采用子集模拟方法求解结构可靠度问题^[14]。相应的子集模拟驱动变量为

然而，在采用子集模拟进行求解时，需要提前确定似然函数乘子c，c的确定对于确定驱动变量g的分布和产生后验样本有不可忽视的影响。为了避免直接估算c，Diazdelao等^[7]提出了改进贝叶斯更新方法（mBUS）方法，将式（6）中的失效区域改写为：

在此基础上，结合自适应条件抽样（aCS）概率反分析空间变异土体水力参数和抗剪强度参数的后验统计特征（均值、标准差和概率分布），具体计算步骤详见文献^{[7, 12]}。

在参数反分析过程中一般需要调用有限元等数值模型进行多次确定性边坡渗流及稳定性分析。例如对于一个接受概率为10⁻⁶量级的反分析事件来说，需要进行10⁵～10⁶次模型调用，这对于复杂模型来说计算量非常可观。为提高计算效率，本文采用基于人工神经网络（ANN）的代理模型提前建立输入参数和输出响应量（边坡压力水头和安全系数）之间的显式函数关系，代替数值模型进行参数反分析。

2.   算例分析

2.1   模型构建和参数设定

2004年9月10—11日台风“海马”席卷台湾，导致9月12日台湾桃园市桃115公路桩号1k+470处发生重大山体滑坡^[15]。根据Wang等^[15]的研究，本文选择了台风“海马”边坡的潜在不稳定区，采用GeoStudio 2012商业软件建立了某公路边坡的二维数值模型。边坡几何模型和有限元网格如图1所示，该边坡长100 m，宽40 m。共划分了3294个节点，3173个有限元单元，单元尺寸为1 m×1 m。模型边界条件概述如下：总水头边界设置在边坡左侧（AB）和右侧（GH），地下水位以上（BC、FG）设为不透水边界。坡顶（DE）使用单位流量边界条件，将降雨数据输入到水力边界函数中，作为模型的降雨边界。由于现场缺乏雨量计，台风“海马”的降雨数据（2004年9月10—11日）取自距离边坡3.8 km的复兴气象站C0C460，图2为2004年9月10—11日的每小时降雨监测数据。最大降雨强度为25 mm/h，累计降雨量为145 mm。

图  1  边坡数值模型示意图

Figure  1.  Illustration diagram of slope numerical model

下载: 全尺寸图片 幻灯片

图  2  2004年9月10—11日降雨数据

Figure  2.  Rainfall data for September 10−11, 2004

下载: 全尺寸图片 幻灯片

土体参数表现出的固有空间变异性，对边坡可靠度会产生重要影响^[16]，因此合理确定空间变异土体参数的统计特征就显得尤为重要。本文土体参数先验统计信息根据文献^{[2, 15 − 16]}确定。由于缺乏非饱和土的水力参数试验数据，本文根据Wang等^[15]的研究，使用GeoStudio 2012中内置的粉砂样本含水量函数来表征边坡的土水特征曲线（SWCC）。需要指出的是，内置函数只需输入饱和体积含水量和相应的土体类型即可，不需要输入其他非饱和渗流参数值。水力传导函数采用Fredlund-Xing模型^[17]。表1给出了土体参数的先验信息。

表  1  回填层岩土体参数先验信息

Table  1.  Prior statistical characteristics of soil parameters for the backfill layer

土体类型	土体参数	μ	σ	COV	概率分布
回填层	饱和渗透系数/（m·h−1）	0.108	0.0648	0.6	对数正态分布
	黏聚力/kPa	10	3	0.3	对数正态分布
	内摩擦角/（°）	15	3	0.2	对数正态分布
	饱和体积含水量	0.315	—	—	—
	饱和重度/（kN·m−3）	20	—	—	—

下载: 导出CSV 

| 显示表格

本文将饱和渗透系数（k_s）、有效黏聚力（c'）和有效内摩擦角（φ'）模拟成平稳对数正态随机场。整个边坡共离散了800个随机场单元，随机场单元尺寸为2 m×2 m，1个随机场单元覆盖了4个尺寸为1 m×1 m的有限元单元，使用高斯型自相关函数来描述参数空间自相关性：

式中：(x₁, y₁)、(x₂, y₂)是随机场空间内任意两点坐标，l_h和l_v分别为水平和垂直自相关距离。水平相关距离和竖直相关距离分别取50 m和10 m^[18]。随机场单元与自相关距离的比值小于0.2，可以满足计算精度要求^[19]。采用Karhunen-Loève（KL）展开方法离散三个平稳对数正态随机场^[16]，当展开项数为12时，随机场的期望能比率因子为0.95，可满足离散精度要求。因此，总共需要36个随机变量来离散k_s、c'、φ'这三个平稳对数正态随机场。

2.2   土体参数概率反分析结果

大多数现场监测数据都十分稀疏、不完整。例如香港岩土工程办公室1999年在香港大屿山东涌进行了一项监测工程。在一个面积大约25000 m²的区域内，仅在1～3 m的浅层区域安装了10个孔压计，很难完整地描述整个区域的压力水头分布情况^[20]。文中公路边坡同样缺失压力水头实地监测数据，因此本文采用一种虚拟监测方案来获取监测数据^[17]。首先基于参数先验信息随机产生一个k_s典型实现，以此作为真实的参数参考场，见图3（a），将参考场输入GeoStudio 2012计算得到降雨历时第48 小时压力水头分布（即2004年9月11日末，见图3（b））。接着根据Zhang等^[21]的研究，在计算的压力水头分布中加入5%的高斯白噪声，得到压力水头虚拟监测数据。此外，将边坡在2004年9月11日末保持稳定作为现场观测信息。然后将虚拟监测数据和实际现场观测信息融合到概率反分析中，推断空间变异土体参数的概率分布，进而预测后续的边坡可靠度。

图  3  饱和渗透系数随机场的典型实现以及相应的渗流分析结果

Figure  3.  A typical realization of the random fields of saturated hydraulic conductivities and corresponding seepage analysis results

下载: 全尺寸图片 幻灯片

本节使用MATLAB神经网络工具箱中的反向传播算法构建了21个基于人工神经网络（ANN）模型的边坡压力水头和安全系数代理模型，以近似计算20个监测点（见图1）边坡压力水头与安全系数和36个独立标准正态变量之间的关系。将1500个随机样本组合以及2004年9月11日末20个监测点处边坡压力水头和安全系数作为训练数据集。其中压力水头主要受k_s影响，因此在建立压力水头代理模型时，只需输入与20个监测点处压力水头有关的12个输入随机变量。由于FS同时受3个参数（k_s、c'、φ'）的影响，因此在建立FS的代理模型时需要输入36个随机变量。此外，额外生成了500个随机样本和对应的有限元输出响应量作为验证数据集，以验证所构建代理模型的有效性。图4比较了使用代理模型和有限元模型计算得出的20个监测点处边坡压力水头值和第48小时安全系数。从图4（a）和（b）可以看出，基于ANN的代理模型得出的20个监测点压力水头值和第48 小时安全系数与有限元计算结果一致。此外，21个代理模型的决定系数（R²）均大于97%。这表明所构建的代理模型具有良好的预测性能，可用于后续的参数概率反分析和边坡可靠度预测，进而评估边坡压力水头和安全系数，从而提高计算效率。

图  4  基于代理模型和有限元方法计算的压力水头和安全系数的比较

Figure  4.  Comparison of pressure heads and factors of safety calculated using the surrogate model and finite element method

下载: 全尺寸图片 幻灯片

本文采用贝叶斯序贯更新方法进行两阶段概率反分析。在第一阶段，利用k_s的先验信息（表1）和虚拟监测数据反分析k_s的后验分布。在第二阶段，将第一阶段的后验分布视为先验分布，并融合现场观测信息进一步推断k_s、c'、φ'的后验分布。为了兼顾计算精度和效率，外层子集模拟的样本数为N_outer = 2000；内层子集模拟的样本数为N_inner = 2000，条件概率p₀ = 0.1。重复进行10次独立的子集模拟计算，取平均值得到最终结果。第一阶段的后验样本总数N_p1 = 20000，第二阶段的后验样本总数N_p2 = 20000，基于这些后验样本推断土体参数的后验统计特征。如此量大的样本可用于对参数后验分布进行全面探索，并对相关统计特征进行稳健估计。

为说明mBUS_aCS方法相较于BUS_aCS方法在解决高维反分析问题方面的优越性，以边坡任意两个空间位置A和B，见图3（a），为例。图5比较了第一阶段概率反分析中由mBUS_aCS和BUS_aCS方法分别计算的不同位置处k_s的后验PDF。BUS_aCS方法与本文使用的mBUS_aCS方法具有相同的外层子集模拟样本数和条件概率。可以发现，使用BUS_aCS方法基于大量监测数据推断出的后验PDF不够平滑，且存在多个峰值。这是因为BUS_aCS方法在处理高维似然函数时，在生成条件样本时会产生很多重复样本，导致样本接受率较小。相比之下，mBUS_aCS方法减少了重复样本的产生，样本接受率高，因此推断得到的后验PDF更平滑。

图  5  边坡两个任意位置的k_s后验PDF比较

Figure  5.  Comparison of the posterior PDFs of k_s at two arbitrary locations on the slope

下载: 全尺寸图片 幻灯片

与先验PDF相比，mBUS_aCS方法推断得到的k_s后验PDF更高瘦，表明参数不确定性显著降低。在位置A处，COV从0.60降至0.38，在位置B处，COV从0.60降至0.40。这表明压力水头对k_s的分布有重要影响。因本文提供了真实的参考场，见图3（a），故可以将推断得出的参数后验分布均值和相应的参考场值（ks_A = 0.055 m/h, ks_B = 0.044 m/h）进行比较。由图5可知，mBUS_aCS方法对应的后验均值与参考场值基本吻合。这表明mBUS_aCS方法可以有效利用大量监测数据推断空间变异土体参数的后验统计特征，并缩减其不确定性。因此，后续进一步利用mBUS_aCS方法对土体水力和抗剪强度参数进行概率反分析。

图6比较了4个监测点位压力水头先验分布和基于第一次更新后的20000个k_s后验样本计算得到的压力水头后验分布。由图6可知，基于后验样本计算得到的4个监测点处压力水头后验均值与真实监测值基本吻合，且压力水头的后验分布均集中在狭窄区域内。表明采用概率反分析方法推断的后验样本计算的响应值与真实监测值吻合较好。

图  6  4个监测点位先验及后验压力水头分布

Figure  6.  The prior and posterior distributions of pressure heads at the 4 monitoring points

下载: 全尺寸图片 幻灯片

在融合监测数据反分析水力参数的基础上，利用观测信息对土体水力参数及抗剪强度参数同时进行第二阶段概率反分析。根据文献^[21]，假定模型偏差系数ζ 服从均值为−0.002，标准差为0.044的正态分布，利用FS>1的现场观测信息构建似然函数，如式（5）所示。将第一阶段融合压力水头监测数据获得的后验分布作为第二阶段概率反分析的先验信息。

图7分别展示了c'和φ'的后验均值和标准差的空间分布。其中颜色的深和浅分别代表土体参数统计特征值的大和小。从图7（c）中可以看出，临界滑移面附近φ'的后验均值最大，约为16.4°，大于先验均值15°。这是因为φ'值越大，边坡越稳定，进而与边坡保持稳定的观测信息一致。在远离临界滑移面的区域，φ'的后验均值几乎等于先验均值。如图7（d）所示，融合现场观测信息后，φ'的标准差从先验的3.0°减小至2.5°。临界滑移面附近的φ'标准差降幅最大。图7（a）和（b）中的c'也有类似的规律。但是，c'的统计量变化明显小于φ'，这表明该边坡稳定性对φ'更为敏感。基于k_s、c'、φ'的后验均值，采用简化Morgenstern-Price法计算的2004年9月11日末边坡FS为1.15，也与现场边坡保持稳定这一观测信息保持一致。

图  7  有效黏聚力和内摩擦角后验均值和标准差空间分布

Figure  7.  Spatial distribution of posterior mean and standard deviation of effective cohesion and internal friction angle

下载: 全尺寸图片 幻灯片

图8为2004年9月11日末融合不同信息量计算的安全系数的PDF图，可以看出，未融合任何信息的安全系数均值偏小（FS = 1.03），且标准差较大，对应的失效概率也更大。融合了压力水头监测数据后安全系数PDF没有明显变化。然而一旦同时融合监测数据和现场观测信息等多源信息后，安全系数的均值（FS = 1.13）向右出现明显偏移，安全系数的PDF更加高瘦，相对应的失效概率也得到了极大的缩减。可以看出，前文中经过贝叶斯概率反分析后土体参数不确定性的缩减进一步导致了安全系数不确定性的降低，进而获得更为可靠的边坡失效概率预测。考虑土体参数空间变异性时先验失效概率为43.7%，此时边坡存在极大可能出现失稳破坏，但是一旦融合边坡保持稳定的现场观测信息，边坡后验失效概率急剧减小为9.3%，在允许模型偏差范围内与现场边坡稳定这一观测信息吻合，也证明了本方法的有效性。

图  8  融入不同信息后的2004年9月11日末安全系数PDF图

Figure  8.  PDFs of FS at the end of September 11, 2004 with incorporating different information

下载: 全尺寸图片 幻灯片

2.3   目标降雨失效概率预测

一旦获得了空间变异土体参数的后验信息，便可预测未来降雨事件下边坡可靠度。图9显示了研究区2004年9月12日暴雨期间的每小时降雨数据。降雨持续时间为24 h，累计降雨量达287 mm。这样的暴雨导致桃园市公路1k+470发生了山体滑坡。利用9月12日的每小时降雨数据和更新后的参数，对2004年暴雨期间的边坡渗流进行了分析，进而预测边坡稳定性。

图  9  边坡先验和后验失效概率随降雨历时的变化

Figure  9.  Variation of prior and posterior probability of slope failure with rainfall duration

下载: 全尺寸图片 幻灯片

图9为2004年暴雨期间边坡先验和后验失效概率随降雨历时的变化图。融合监测数据和观测信息之后，边坡失效概率由9月12日初的9.3%急剧增大至9月12日末的23.1%，如此高的失效概率被认为是极其危险的，需要采取紧急措施来减轻灾害。这些结果表明，在该地区极有可能发生边坡失稳。这与2004年9月12日在该边坡附近及周边地区观测到的边坡失效情况吻合^[15]。综上，如果仍使用参数先验信息进行边坡失效概率预测，就会高估边坡失稳的概率，进而会错误预测发生滑坡的时间。

3.   结论与讨论

本文基于mBUS_aCS方法构建了融合边坡监测数据和现场观测信息等多源信息的土体参数概率反分析及边坡可靠度预测的计算框架。以某公路边坡为例，说明所构建框架的有效性和合理性。主要结论如下：

（1）mBUS_aCS方法可以充分利用监测数据、观测信息等多源场地信息，推断空间变异土体参数后验分布并预测降雨条件下边坡可靠度。与BUS_aCS方法相比，mBUS_aCS方法生成的后验样本分布更均匀，计算精度更高。因此，该方法为解决低接受概率水平的挑战和克服概率反分析和可靠性预测中遇到的高维似然函数问题提供了有价值的研究思路。

（2）融合监测数据和观测信息后，土体参数的不确定性明显降低，其中受压力水头影响的渗透系数（k_s）不确定性降低最为明显，而由现场观测信息反分析的有效黏聚力（c'）和有效内摩擦角（φ'）不确定性仅略有降低。更新后的土体参数可用于准确预测2004年9月12日因暴雨引发的山体滑坡事件。

（3）值得注意的是，本文采用的是虚拟的监测数据，如果有真实的压力水头监测数据，则参数概率反分析结果更具有说服力。并且本文融合的是当前阶段的所有信息来概率反分析土体参数，关于贝叶斯更新中不同时段的更新策略问题有待进一步的研究。

（4）本文并未考虑边坡失稳机制对概率反分析的影响，如未将后缘拉张裂隙对边坡稳定性的影响考虑入内，这也是本文下一步的研究方向。后续研究需要采用更能反映边坡实际情况、更具说服力的边坡稳定性模型融入到参数反分析中，进一步拓宽本文提出方法的适用范围。