目前，滑坡稳定性的定量计算方法可分为确定性计算方法和可靠度计算方法。传统的确定性分析方法主要采用极限平衡分析方法、强度折减法等计算坡体的安全系数，在实际工程中得到了广泛的应用。但由于坡体存在荷载环境、结构特征、岩土体特性参数等诸多不确定性因素，采用以概率统计为基础的可靠度分析方法计算滑坡的破坏概率和可靠度，比用安全系数在一定程度上能更客观、定量地反映滑坡的安全性^[1]，并可以为滑坡灾害的定量风险评价提供依据^[2]，近年来在滑坡工程中日益得到重视。

2018年10月10日和11月3日，西藏自治区江达县波罗乡白格村金沙江右岸先后两次发生大规模高位滑坡，堵塞金沙江，堰塞湖与溃坝洪水给金沙江上游沿岸居民及其生产和生活设施带来巨大灾害。白格滑坡及其导致的滑坡堵江–溃坝洪水灾害链受到了国内外的广泛关注，诸多学者对滑坡的地质环境条件、变形破坏机制及发展趋势进行了研究。许强等^[3]结合现场调查和历史遥感影像解译等手段，分析了滑坡区斜坡的变形历史和动态演化过程，以及两次滑坡–堰塞堵江的基本特征。邓建辉等^[4]基于现场调查对滑坡的形成机制与过程进行了系统分析。张永双等^[5]从地质角度剖析了白格滑坡的地质成因，提出了金沙江构造缝合带混杂岩体岸坡在持续重力作用下的失稳机理。自2018年连续两次滑坡发生之后，白格滑坡后缘出现3个裂缝区（图1），一旦失稳，再次堵江的危险性较大^[4]。因此，裂缝区的发展趋势和稳定性成为目前面临的关键难题。冯文凯等^[6]基于裂缝区形貌特征和变形迹象对其变形破坏特征和发展趋势进行了分析。陈菲等^[7]结合地质和监测资料对裂缝区的发展趋势进行了分析和预测。王立朝^[8]和Fan^[9]等分别采用2维离散元模拟软件UDEC和3维有限差分软件FLAC3D对裂缝区稳定性进行了定量数值模拟，由于缺乏现场试验参数，定量分析采用的均为经验参数，难以反映滑坡岩土体的真实性状。

图1(Fig. 1)

图1 白格滑坡平面图 Fig. 1 Overview of Baige landslide

鉴于白格滑坡地质结构的复杂性和岩土体物理力学参数的不确定性，采用可靠度方法对白格滑坡进行参数反演和裂缝区失稳概率分析，对于滑坡灾害的风险评价及处置具有重要意义。

白格滑坡地处青藏高原东部，属横断山区北部，芒康山与沙鲁里山间的金沙江河谷地带，为典型的构造高山峡谷侵蚀地貌。滑坡位于金沙江右岸，前缘为金沙江凹岸，平水位高程2 880 m；后缘为一走向N19°E的条形山脊，高程约3 730 m。滑坡坡向80°～100°，坡度上下陡、中间缓，相对高差达850 m，平均坡度约33°。

滑坡区位于金沙江缝合带上，由于受多期构造运动的影响，区域构造形迹较为复杂。主要断裂构造形迹为近北南向的江达—波罗金沙江断裂带及NW向断裂。根据现场调查及勘查报告^[10]，滑坡区出露岩性以片麻岩为主，中前部为片麻岩夹碳质板岩，岩层产状倾向（180°～220°）、倾角（∠36°～∠42°）；中部到后缘范围为蛇纹石化片麻岩，在斜坡顶部结晶良好，墨绿色；中下部呈碎粉岩状，灰绿色，绿泥石化，具有多期、多次变形与变质特点。滑坡体典型剖面如图2所示。

图2(Fig. 2)

图2 白格滑坡剖面图（1–1′剖面） Fig. 2 Cross section of Baige landslide (1–1′ section)

滑坡按高程划分为3区，即前缘的阻滑区、中部的主滑区和后缘的牵引区，分割高程大致为3 500 m和3 000 m。主滑区为楔形体，系两组发育良好的结构面切割形成；阻滑区为四面体，由两组发育较差的结构面切割形成；牵引区为完全风化的岩土体夹团块状碎裂岩体。白格“10·10”滑坡是相对完整的阻滑区岩体在主滑区重力驱动下的渐进破坏过程，方量约10⁷ m³。主滑区楔形体重力是滑坡的主要动力来源，滑坡的孕育过程中，白格“11·3”滑坡是牵引区的部分岩土体在起阻滑作用的碎裂岩体渐进解体后下滑的结果，方量约3×10⁶ m^3[2]。

在连续两次滑坡事件发生之后，受失稳滑体牵引拖拽作用以及临空卸荷作用的影响，在滑坡体后缘及两侧出现多条深大裂缝，形成3个裂缝区（图1）。综合现场调查、监测及地质分析^[5]，各区域的范围及发展趋势为：C1区位于滑坡体下游侧边界，面积约2.3×10⁵ m²，厚度15～95 m，可分为3个亚区，变形速率最快的为C1–1区，体积约7×10⁵ m³；C2区位于滑坡体后缘，面积约7.3×10⁴ m²，厚度45～70 m，可分为2个亚区，潜在失稳范围为C2–1区，体积约3.2×10⁶ m³；C3区位于滑坡体上游侧边界，面积约1.2×10⁵ m²，厚度20～30 m，体积约1.9×10⁶ m³，可分为2个亚区，其中C3–1区体积约1.3×10⁶ m³。各裂缝区典型地质剖面如图3所示，ZK表示钻孔位置和编号。

图3(Fig. 3)

图3 裂缝区剖面图 Fig. 3 Cross section of crack zones

在滑坡可靠度分析中，滑坡的破坏概率与坡体抗剪强度参数的统计特征密切相关。受现场条件及试验方法限制，试验数量有限，试验成果往往难以准确地反映坡体抗剪强度参数的统计特征。贝叶斯分析方法能够将岩土体参数先验信息和有限的场地信息有机结合，在概率反演边坡抗剪强度参数和渗透系数等方面逐渐得到应用^[11]。Zhang等^[12]提出了基于极大似然估计的优化算法，优化过程可以在电子表格实现。Zhang等^[13]采用马尔可夫链蒙特卡洛（MCMC）方法对多层土体抗剪强度参数同时进行反分析，有效地解决了边坡反分析中多参数反演问题。Wang等^[14]分别采用MCMC方法和极大似然估计对台湾3号高速公路滑坡抗剪强度参数和锚固力参数进行了非确定性反演，表明上述两种方法的计算结果均与现场调查和勘测结果一致，基于极大似然估计的优化算法计算工作量小，更适用于滑坡事后调查和评估。基于极大似然估计的优化算法基本原理如下：

设 $g({{\theta }})$ 为滑坡安全系数的计算模型； ${{\theta }}$ 为随机参数向量，其先验分布满足均值 ${\;\mu _{{\theta}} }$ 、协方差 ${C_{{\theta}} }$ 的正态分布。由于计算模型的不确定性，模型预测值与滑坡实际安全系数之间的误差采用随机变量表示：

式中： $y$ 为滑坡实际安全系数； $\varepsilon $ 为表征模型不确定性的随机变量，为简化起见，设 $\varepsilon $ 服从正态分布 $N({\mu _\varepsilon },\sigma _\varepsilon ^2)$ 。

基于以上假设， ${{\theta}} $ 的后验分布也近似为均值 ${{{\;\mu}} _{{ \theta }\left| d \right.}}$ 、协方差 ${{{C}}_{{ \theta }\left| d \right.}}$ 的多元正态分布，滑坡稳定性概率反分析的目标即为根据滑坡的稳定性状态，确定均值矩阵 ${{{\;\mu}} _{{ \theta }\left| d \right.}}$ 和协方差矩阵 ${{{C}}_{{ \theta }\left| d \right.}}$ 。

对于出现失稳的滑坡， ${{{\;\mu}} _{{ \theta }\left| d \right.}}$ 也是导致滑坡破坏事件的最可能参数组合，其极大似然估计可通过最小化负对数似然函数来实现：

随机参数向量 ${{\theta}} $ 的后验协方差矩阵 ${{{C}}_{{{\theta}} \left| d \right.}}$ 表示各向量元素之间的相关性，可用式（3）和（4）确定。

Zhang等^[12]将上述方法称为基于优化的反分析方法，其优化过程可以在Excel中实现，计算简便。

边坡稳定性反分析一般是极限平衡分析方法。在极限平衡分析方法中，Morgenstern–Price法满足严格平衡条件，收敛性良好，并能对任意形状的滑动面进行分析，因此得到了广泛应用。在概率反分析优化过程中，需要重复调用安全系数的计算模型并对其进行求导，而Morgenstern–Price法的安全系数计算公式为隐函数，无法写出显式状态方程^[15]，其使用具有一定局限。为克服这一问题，采用响应面法（response surface method, RSM）对安全系数计算模型进行拟合。

现场调查及监测表明，白格“11·3”滑坡坡体内无明显地下水，影响滑坡稳定性最主要的参数是滑坡岩土体的抗剪强度指标黏聚力 $c$ 和内摩擦角 $\tan \; \phi $ 。为提高安全系数计算模型的线性度，本文选用 $c$ 和 $\tan \; \phi $ 作为基本随机参数。根据响应面法基本原理^[16-17]，将安全系数近似表达为岩土体抗剪强度参数的二次函数形式，有：

式中： ${{{\theta}} _i}$ 为随机参数向量，即岩土体的抗剪强度指标 $c$ 、 $\tan \; \phi $ ； ${a_0}$ 、 ${b_i}$ 、 ${c_i}$ 为待定系数。针对白格“11·3”滑坡，采用“轴向点+中心点”的取样组合进行取样，在取样点利用GeoStudio软件^[18]，采用Morgenstern–Price极限平衡法进行安全系数计算。极限平衡计算模型见图4，模型中滑坡岩土体容重为确定性参数，取为 $\gamma = 24\; {\rm{kN}}/{{\rm{m}}^3}$ 。取样点及相应的安全系数见表1。

表1(Tab. 1)

表1 中心复合设计取样点及安全系数 Tab. 1 Sampling points by central composite design and corresponding safety factor

表1 中心复合设计取样点及安全系数 Tab. 1 Sampling points by central composite design and corresponding safety factor 序号 $c/{\rm{kPa}}$ $\tan \; \phi $ $g({ \theta} )$ 1 15.00 0.65 0.952 2 10.76 0.65 0.943 3 15.00 0.51 0.752 4 19.24 0.65 0.961 5 15.00 0.79 1.152

图4(Fig. 4)

图4 白格 “11·3”滑坡安全系数计算模型 Fig. 4 LEM model of Baige “11·3” landslide

将表1中的数据代入式（5），进行响应面求解，得到滑坡安全系数计算模型的响应面近似方程为：

式中， $c$ 和 $\tan \; \phi $ 的二次项系数均为0，表明该滑坡安全系数和岩土体抗剪强度参数近似为线性关系。同时，由于该滑坡的平均厚度较大， $\tan \; \phi $ 对该滑坡安全系数的贡献较大，因此 $\tan \; \phi $ 的一次项系数远大于 $c$ 的一次项系数。

根据坡体强风化蛇纹岩的室内直剪试验结果（表2）及相关文献[19]，随机参数向量的先验分布的黏聚力和内摩擦角分别可以取为 $c{\text{～}} N(15,{3^2})$ ， $\tan \; \phi {\text{～}} N (0.65,{0.1^2})$ 。设随机参数向量满足独立对数正态分布，协方差矩阵为：

表2(Tab. 2)

表2 直剪试验结果 Tab. 2 Results of direct shear test

表2 直剪试验结果 Tab. 2 Results of direct shear test 序号 干燥固结快剪 饱和固结快剪 c/kPa $ \tan \; \phi $ c/kPa $ \tan \; \phi $ 1 31.6 0.86 4.2 0.79 2 29.3 0.84 3.2 0.81 平均值 30.5 0.85 3.7 0.80

根据Christian等^[20]的研究，简化Bishop法的模型误差服从均值0.05、标准差0.07的正态分布。考虑到Morgenstern–Price法与简化Bishop法的计算结果一般比较接近，且Morgenstern–Price法满足力的严格平衡条件，因此假设安全系数计算模型 $g({{\theta}} )$ 的误差满足正态分布，均值、标准差分别取为 ${\;\mu _\varepsilon } = 0.05$ ， ${\sigma _\varepsilon } = 0.07$ 。

根据第2.2.1和2.2.2节计算模型和参数，采用Excel建立式（2）～（4）的电子计算表格；并采用Excel的规划求解功能，通过使式（2）的 $S\!({{\theta }})$ 最小化得到 ${{\theta}} $ 的后验均值 $\;{{{\mu}} _{{{\theta}} \left| d \right.}}$ ，通过式（3）和（4）计算后验分布的协方差 ${{{C}}_{{{\theta }}\left| d \right.}}$ 、标准差 ${{{\sigma}} _{{{\theta}} \left| d \right.}}$ 。

经过优化求解，求得随机参数向量的后验分布黏聚力和内摩擦角可分别取为 $c{\text{～}} N(13.05,{3^2})$ ， $\tan \; \phi {\text{～}} N(0.64,{0.04^2})$ 。图5为随机参数向量的先验分布与后验分布的对比。由图5可见：抗剪强度黏聚力、内摩擦角参数通过反分析可以降低先验参数的不确定性。由于 $\tan \; \phi $ 对边坡安全系数的贡献大，其概率密度的改变程度也相应较大。由于试验方法为室内直剪试验，且试验数据较少，反演参数更能反映坡体材料参数的统计特征。

图5(Fig. 5)

图5 随机参数向量先验分布与后验分布对比 Fig. 5 Comparison of prior and posterior distribution

由于白格“11·3”滑坡滑源区与后缘裂缝区位于同一高程，岩性基本一致。因此，采用第2.2.3节白格“11·3”滑坡的参数反演分析结果，对后缘各裂缝区的失稳概率进行可靠度分析。计算方法为蒙特卡洛法，采用GeoStudio软件进行计算。计算时，针对每个典型剖面（图3）均考虑两种破坏模式：一是，整体破坏模式，以勘探及监测确定的潜在滑动面作为确定性滑面进行计算；二是，局部破坏模式，由软件自动搜索最危险滑动面，如图6所示。

图6(Fig. 6)

图6 白格滑坡后缘裂缝区潜在破坏模式 Fig. 6 Failure mode of potential failure mode in the rear fracture zones of Baige landslide

由图6可知：裂缝区C1自动搜索的最危险滑动区范围与监测得到的变形速率最大的C1–1区范围基本一致。同样，裂缝区C2、C3自动搜索的最危险滑动区范围也分别对应裂缝区前部的潜在失稳区C2–1和C3–1。

表3为各裂缝区的破坏模式及相应的失稳概率计算结果。

表3(Tab. 3)

表3 白格“11·3”滑坡后缘裂缝区可靠度计算结果 Tab. 3 Reliability calculation results of the fracture zone at the rear edge of Baige “11·3” landslide

表3 白格“11·3”滑坡后缘裂缝区可靠度计算结果 Tab. 3 Reliability calculation results of the fracture zone at the rear edge of Baige “11·3” landslide 编号 滑动模式 可靠度指标 $\;\beta $ 失稳概率pf/% 滑坡危险性 C1 整体 2.02 1.80 高 局部 0.79 21.70 极高 C2 整体 1.78 3.05 高 局部 0.44 33.90 极高 C3 整体 1.19 10.20 高 局部 0.64 27.30 极高

由表3可知：各裂缝区在局部破坏模式下，其破坏概率均明显大于整体破坏模式，表明裂缝区破坏以渐进牵引破坏为主，与实际的裂缝发展规律及破坏模式一致。根据文献[9]的滑坡半定量标准，可通过失稳概率判断滑坡的危险性： ${p_{\rm f}} \ge 5.00{\text{%}}$ 为极高危险性， $1.00{\text{%}} \le {p_{\rm f}} < 5。00{\text{%}}$ 为高危险性， $0.20\% \le {p_{\rm f}} < 1.00{\text{%}}$ 为中危险性， $0.04{\text{%}} \le {p_{\rm f}} < 0.20{\text{%}}$ 为低危险性， ${p_{\rm f}} < 0.04{\text{%}}$ 为极低危险性。据此和表3结果可以判断，3个裂缝区局部均为极高危险性，整体为高危险性。特别是方量较大、失稳概率最高的C2–1区如再次启动，堵江的危险性较大。

针对白格“11·3”滑坡，采用基于极大似然估计的优化算法对滑坡岩土体的抗剪强度参数进行了反演，并利用反演参数，对滑坡后缘裂缝区进行了可靠度分析，主要结果如下：

1）本文将响应面法与基于极大似然的概率反分析方法结合起来，对白格“11·3”滑坡滑坡体不确定参数进行反演，计算过程简便，计算方法适应性强。

2）通过对白格“11·3”滑坡抗剪强度参数进行反演，得到滑坡岩土体抗剪强度后验分布的黏聚力和内摩擦角分别为 $c{\text{～}} N(13.05{\rm{\;kPa}},{3^2})$ 、 $\tan \; \phi {\text{～}} N\left(0.64, \right.$ $\left.{0.04^2}\right)$ ，计算结果表明，通过反演可以降低先验参数的不确定性。

3）采用反演参数对白格“11·3”滑坡后缘3个裂缝区的整体和局部稳定性进行了可靠度分析。各裂缝区在局部破坏模式下的最危险滑动区范围与现场调查及监测得到的失稳区范围基本一致，且局部破坏模式下的破坏概率均明显大于整体破坏模式，表明裂缝区破坏以渐进牵引破坏为主，与实际的裂缝发展规律及破坏模式一致。

4）白格“11·3”滑坡后缘3个裂缝区的潜在失稳区域C1–1、C2–1和C3–1的体积分别为7×10⁵ m³、3.2×10⁶ m³和1.3×10⁶ m³，对应的失稳概率分别为21.70%、33.90%和27.30%，均为极高危险性。类比白格“11·3” 滑坡堵江事件，裂缝区如再次启动，特别是方量较大、失稳概率最高的C2–1区如再次启动，堵江的危险性较大。有必要根据可靠度计算结果，对滑坡堵江–溃坝灾害链进行风险评价，并提出合理的处置措施和应急预案。

5）鉴于白格“11·3”滑坡后缘3个裂缝区已逐步连通，裂缝区之间的变形和稳定性存在相互影响。如何考虑裂缝区之间的相互作用机制，从系统可靠度的角度对其总体稳定性进行评价，有待进一步进行研究。