0.   引言

降雨是滑坡灾害的重要诱因^[1]。随着全球气候变化，极端降雨事件的频率与强度增大，诱发群发性浅层滑坡的问题在我国东南地区越来越突出，已造成巨大的人员伤亡和经济损失。例如，2021年6月10日，诸暨市遭遇突发性极端降雨天气，诱发了701处滑坡，造成4人遇难，经济损失高达7.5亿元。极端降雨诱发浅层滑坡的危险性评估已成为我国东南丘陵区地质灾害防控领域亟待解决的难题之一。

针对降雨型滑坡的危险性评估问题，目前国内外学者大多选取具有代表性的特定极端降雨事件为工况，计算安全系数（F_s）以评估斜坡稳定性^[2]，该方法仅能提供确定性的评估结果，未考虑降雨事件的随机过程特征。事实上，降雨强度、持续时间等降雨特征在时间分布上具有很强的不确定性，这种不确定性应被纳入滑坡危险性评估中。因此，有研究通过计算滑坡重现期以评估滑坡概率，此方法基于多年降雨数据，能够准确表征降雨特征的随机性。例如，Lu等^[3]采用降雨强度与持续时间的二维联合分布模拟年降雨量的不确定性，结合斜坡地质条件计算边坡年失效概率；He等^[4]提出了一种基于降雨事件交替随机更新过程的滑坡概率全寿命周期评价方法；Peres等^[5]利用蒙特卡罗方法将随机降雨模型与斜坡水文模型结合，估算了滑坡重现期。这类基于随机降雨过程的滑坡重现期计算方法，能够实现长时间尺度下降雨诱发滑坡的危险性定量评估^{[6 − 9]}。

本文在考虑降雨随机性的基础上重点刻画了极端降雨事件的特征，提出一种极端降雨诱发浅层滑坡的重现期计算方法：首先，对降雨数据进行统计分析，基于Copula函数构建极端降雨事件平均强度-持续时间（I-D）的二元联合分布；其次，考虑土体力学性质与降雨输入，结合非饱和渗流与极限平衡分析，模拟单体斜坡失稳的I-D阈值；最后，利用蒙特卡罗抽样方法模拟极端降雨事件I-D联合分布超过滑坡I-D阈值的概率，表征单次降雨导致特定斜坡失效的概率$ P_{F 1} $，能进一步计算斜坡年失效概率$ P_{F A} $及滑坡重现期。该方法考虑了斜坡岩土体的特定力学性能以及降雨的不确定性，重现期计算结果能动态评估滑坡危险性，为风险处置以及地质灾害综合防治工作提供参考。

1.   研究区概况

本文以浙江省诸暨市马剑镇某滑坡隐患点为例展开研究。马剑镇位于诸暨市西侧，年平均降水量为1801.36 mm。2021年6月10日诸暨市遭受极端降雨天气期间，马剑镇发生86处地质灾害，占总历史灾害数量的92.47%，因此该地区极端降雨条件下的滑坡危险性评估工作亟待开展。

研究选取的滑坡隐患点如图1所示。该斜坡属于侵蚀剥蚀丘陵地貌，地形较陡，整体呈北西-南东向延伸。坡高差约125 m，坡度34°。坡面坡脚无地表水，降雨主要沿斜坡汇入坡底。根据工程勘探，该斜坡覆盖层主要为含碎石砂土、含碎石粉质黏土、全风化～强风化下白垩统流纹斑岩（$ \lambda {\textit{π}} K_{1} $）。本研究选取a-a剖面进行二维边坡稳定性分析，具体见2.2节。

图  1  研究斜坡示意图

Figure  1.  Profile map of the of the investigated slope

下载: 全尺寸图片 幻灯片

2.   研究内容

文章提出了一个极端降雨诱发浅层滑坡的重现期计算框架，以浙江省马剑镇某隐患斜坡为例验证思路可行性。基于浙江省气象局提供的马剑站2005—2021年近17年的小时降雨数据，考虑该斜坡的土体力学性质，计算能够表征斜坡动态危险性的滑坡重现期。具体框架如图2所示。

图  2  研究方法

Figure  2.  Research methodology

下载: 全尺寸图片 幻灯片

2.1   极端降雨事件I-D联合分布

平均降雨强度I和持续时间D是反映极端降雨事件的最主要特征。为完整刻画极端降雨事件，本文考虑两者的统计相关性，采用Copula函数构建马剑镇极端降雨的I-D联合分布模型。

2.1.1   极端降雨事件提取

极端降雨事件的选样方法多样^[10]，此研究基于以下规则从马剑站的小时降雨数据中提取降雨事件并进行极端降雨选样：①保证事件间间隔（IETD）大于等于6小时，且每场降雨事件的强度大于等于强度阈值$ I_{0}=1 \mathrm{～mm} \cdot h^{-1} $，提取每场降雨事件的持续时间D与总雨量E，并计算平均降雨强度I；②利用峰值超过阈值法（POT）法进行极端降雨事件取样，确保每年至少有一场降雨事件被纳入分析^[11],使用累计降雨量的90%分位数作为阈值，确保数据拟合的准确性。

2.1.2   边际分布的构建

（1）边际分布函数的优选

初选表1中的5个分布函数用于拟合I和D。计算5种分布下I和D的赤池信息准则（AIC）值及K-S拟合优度指标，优选分布函数。其中，AIC的计算公式如式（1）：

表  1  马剑站边际分布拟合优度检验结果

Table  1.  Results of the goodness-of-fit test for the marginal distribution of the Majian station

	降雨强度（I）
	Lognormal	GEV	GPD	Exponential	Gumbel
AIC	834.9	832.4	850.5	856.5	920.2
统计量D	0.0678	0.0673	0.0961	0.1239	0.1517
P-value	0.4239	0.4322	0.0921	0.0122	0.0010
	降雨持续时间（D）
	Lognormal	GEV	GPD	Exponential	Gumbel
AIC	1231.4	1233.9	1244.9	1250.3	1234.6
统计量D	0.0538	0.0544	0.1253	0.1562	0.0583
P-value	0.7117	0.7001	0.0109	0.0006	0.6163

下载: 导出CSV 

| 显示表格

式中：$ k $——所拟合分布中参数的数量；

$ L $——对数似然值。

AIC值越小及 K-S检验中的D值越小且P-value越接近1，表明边际分布的拟合效果越好，检验结果如表1所示。

分析表中结果，优选GEV分布和Lognormal分布分别用于拟合I和D的边际分布。

（2）基于MCMC方法的边际分布参数估计

常见的参数估计方法有矩量法、概率加权矩量法、I矩量法和最大似然估计，然而，这些方法仅给出了参数的单点估计，无法传达参数估计中持续存在的不确定性，且在小样本的参数估计中准确性不高。贝叶斯推理提供了有力的框架来估计非平稳模型的后验分布，并量化预测的不确定性，是分析稀疏数据的首选方法之一^[12]。在贝叶斯推理中，马尔可夫链蒙特卡罗算法（MCMC）通常用于计算后验概率分布的矩和可信区间。本研究优选MCMC算法对边际分布参数的后验分布进行稳健估计。

在拟合I的GEV分布参数估计中，三个参数的先验均选用正态分布，正态分布的均值设置为利用极大似然法计算的参数值，方差设置为2，即：

在拟合D的Lognormal分布参数估计中，参考Kwon H及Sarhadi等的研究^{[12 − 13]}，形状参数$ s $先验选用Gamma分布，另外两个参数先验选用正态分布，方差设置为2，即：

式中，$ { mle {\mathrm{params}} } $——极大似然法计算的参数值。

边际分布的参数估计结果见表2。

表  2  边际分布与Copula函数参数值

Table  2.  Marginal distribution and Copula function parameter values

I（GEV）		D（Lognormal）		Copula函数单参数$ \theta $
$ \text { E } $	−0.52	$ 5 $	0.54	−7.8903
$ \mu $	3.34	$ \operatorname{loc} $	−4.57
$ \sigma $	2.01	$ \text { scale } $	19.32

下载: 导出CSV 

| 显示表格

2.1.3   Copula函数的构建

构建Copula函数之前，采用Kendall秩相关和Spearman秩相关系数检验极端降雨事件样本中I与D的相关性，结果显示两者间存在较强的负相关，因此优选Frank Copula构建I-D联合分布函数，其分布见式（2）：

采用极大似然法进行参数估计，计算得到Frank Copula的唯一参数$ \theta $值为−7.8903，边际分布与Copula函数的参数计算结果见表2，Copula的概率密度函数见图3。

图  3  Copula函数的概率密度函数PDF图

Figure  3.  Probability Density Function of Copula Function

下载: 全尺寸图片 幻灯片

由图3可得平均降雨强度I与持续时间D之间存在显著的负相关关系，强降雨往往以短历时形式出现。同时，随着I或D的增大，即降雨事件特征趋于极端时，联合概率密度逐渐减小，对应的降雨事件发生概率显著降低。

2.1.4   I-D联合分布的构建

结合边际分布与Copula函数框架，能构建I-D二元联合分布，作为马剑镇极端降雨事件特征的概化模型，如式（3）所示：

式中，$ c\left(u_{1}, u_{2}\right) $——Copula的概率密度函数；

$ u_{1}=F_{I}(I) $——I的边际累积分布函数；

$ u_{2}=F_{D}(D) $——D的边际累积分布函数；

$ f_{I}(I) $——I的边际概率密度函数;

$ f_{D}(D) $——D的边际概率密度函数。

2.2   滑坡I-D阈值的确定

参考工程地质勘察结果，在GeoStudio中对研究斜坡的a-a剖面建模，模型如图4所示。

图  4  斜坡剖面模型

Figure  4.  Profile modeling of the investigated slope

下载: 全尺寸图片 幻灯片

在模型中输入斜坡土体参数。土体参数包括抗剪强度参数和水力参数，通过试验获得。其中水力参数包括饱和渗透系数以及土水特征曲线。对于非饱和土体，其基质吸力的大小是土体含水量的函数，含水量和基质吸力的关系曲线称为土水特征曲线（SWCC）。土体的渗透系数为饱和度或体积含水量的函数，因此渗透系数也是基质吸力的函数。计算得到SWCC曲线及渗透系数函数如图5所示。

图  5  斜坡SWCC曲线及渗透系数函数

Figure  5.  SWCC curve and permeability coefficient function of the investigated slope

下载: 全尺寸图片 幻灯片

汇总表征斜坡岩土体水力和力学特性的参数，如表3所示。为简化分析，下覆风化基岩的饱和渗透系数取土体饱和渗透系数的1/100，即2×10⁻⁷ m·s⁻¹，其余参数与表中土体参数相同。

表  3  斜坡土体参数表

Table  3.  Soil parameters of the investigated slope

岩土体水力参数			岩土体抗剪强度参数
饱和体积含水率	50%		饱和重度	21 （kN·m−3）
残余体积含水率	5%		有效粘聚力	13 （kPa）
饱和渗透系数	2×10−5 （m·s−1）		有效内摩擦角	18 （°）

下载: 导出CSV 

| 显示表格

考虑平均降雨强度I为5、10、20、30、40、50、60 mm·h⁻¹的均匀型降雨工况，分别设置于入渗边界。采用基于Richards理论的二维非饱和渗流场模拟和Morgenstern-Price条分法计算斜坡的时变安全系数。记录安全系数到达指定阈值F_S=1所需的降雨持续时间D，拟合7组模拟结果，得到滑坡的I-D阈值曲线如图6所示。

图  6  I-D降雨阈值曲线

Figure  6.  I-D rainfall threshold curve

下载: 全尺寸图片 幻灯片

2.3   滑坡重现期的计算

当降雨事件的平均雨强I和持续时间D联合坐标（I，D）位于降雨阈值曲线的上方时，表明斜坡会发生失稳，其发生的概率等价于单次降雨导致斜坡失效的概率$ P_{F 1} $^[14]，文章采用蒙特卡罗模拟方法计算$ P_{F 1} $，根据I-D二元联合概率密度函数式（3），随机生成$ N $次降雨事件$ X_{i}=(i=1,2, \cdots, k) $，其中每一个$ \boldsymbol{X}_{i} $包含一组$ \left(I_{i}, D_{i}\right) $变量，与降雨阈值图6比对，若模拟中有$ n $次降雨事件超出阈值，则$ P_{F 1} $可以表示为式（5）：

假设一年内多次降雨事件的发生是独立的，极端降雨事件的年发生次数遵循泊松分布^[15]：

式中：$ P(N=k) $——一年内发生$ k $次降雨事件的概率， $ k $是非负整数；

$ \lambda $——泊松分布的正参数，等于泊松分布的均值和方差。

泊松分布的$ \lambda $可视作年平均极端降雨事件发生次数，本案例中$ \lambda $的拟合结果为9.56次。

参考Liu等^[3]的研究可知，年失效概率$ P_{F A} $与单次降雨失效概率$ P_{F 1} $间的关系如式（7）所示：

对应的滑坡重现期见式（8）：

计算出研究斜坡的年失效概率$ P_{F A} $为5.85%，滑坡重现期为17年。

3.   结果与讨论

3.1   案例斜坡滑坡危险性分析

（1）地质条件与斜坡稳定性分析

本文选取的隐患斜坡整体剖面形态为凹形坡，可能导致降雨集中汇流，增加滑坡风险；钻探揭露基岩节理裂隙发育，岩芯较破碎，表明岩体强度较低，抗剪强度较差；切坡建房导致斜坡局部应力集中，增加了斜坡失稳的可能性，可能成为滑坡的潜在触发点。

（2）极端降雨诱发滑坡的重现期分析

本文以滑坡重现期作为评估滑坡危险性的指标。滑坡重现期是指在特定区域内，某一特定强度或规模的滑坡事件发生的平均时间间隔。其估计结果并不意味着在重现期内斜坡必须保持稳定或在重现期之后斜坡就会失效，斜坡失稳在任何年份都可能发生。

通过本文提出的滑坡重现期计算框架，得到极端降雨诱发该斜坡滑坡的重现期为17年，意味着该斜坡在未来17年内有较高的概率发生一次极端降雨诱发的滑坡事件。

（3）滑坡危险性综合评价

案例斜坡坡脚有公路与居民区，屋前有便道连接村道，坡脚为公路，滑坡发生会对居民区造成直接威胁，引发一定程度上的人员伤亡、财产损失和环境破坏。根据重现期计算结果和地质条件分析，参考相关风险评价调查规范，该斜坡在极端降雨条件下具有中等滑坡风险。宜采取监测预警、工程治理等相关的防治措施以降低灾害可能造成的损失。

3.2   研究方法局限性

（1）本研究基于历史小时降雨数据构建极端降雨事件I-D联合分布。降雨数据长度为17年，不足以表征长时间尺度下降雨事件的非平稳特征。因此预测结果在未来的外推精度与适用时间范围有限，保守估计其仅适用于未来20年以内的滑坡危险性评估工作。

（2）斜坡土体力学性质理论上会随着降雨入渗而动态折减，但由于模型与计算的复杂性，本研究未考虑土壤力学参数的不确定性。

4.   结论

本文提出了一种计算斜坡年失效概率$ P_{P A} $与滑坡重现期的方法，该方法基于历史降雨数据，考虑了极端降雨事件的不确定性，实现了斜坡长期稳定性的定量评估，主要结论如下：

（1）重点考虑对极端降雨事件的定义与完整刻画，提出基于Copula函数构建平均降雨强度-持续时间（I-D）二元联合分布的方法，作为表征极端降雨事件特征的概化模型。

（2）综合考虑斜坡岩土体力学特性与降雨条件，结合降雨输入与斜坡稳定性评价，提出了单个斜坡I-D阈值曲线的构建方法。

（3）基于获得的I-D联合分布函数与单个斜坡I-D阈值曲线，采用蒙特卡罗模拟方法，建立了滑坡重现期的动态估计框架。

（4）以浙江省诸暨市马剑镇某滑坡隐患点为例，验证了本文方法的可行性。计算得到该斜坡极端降雨条件下的失效概率为0.63%、年失效概率为5.85%、滑坡重现期为17年，上述结果可为降雨型滑坡的风险防控提供参考。

本文提出的极端降雨事件特征概化模型考虑了降雨过程的随机性，能反映一定时间尺度下当地的降雨特征。而全球气候变化导致极端降雨事件发生频率与强度呈现非平稳特征，如何量化气候变化对极端降雨的影响，实现更长时间尺度下的滑坡重现期评估，值得进一步探讨。