Verhulst反函数预测模型的改进及滑坡时间概率预测

Model Modification of Verhulst Inverse-Function Forecasting Model and Probabilistic Forecast for Landslide Failure Time

摘要

滑坡时间预测是滑坡灾害防治的重要组成部分，但由于滑坡演化存在不确定性，准确预测滑坡的发生极为困难.Verhulst反函数模型是一种常用的滑坡时间预测模型，但其存在计算起始时刻选择不当会造成监测数据拟合质量差和预测精度低的问题.针对这一不足，提出了一种改进的Verhulst反函数模型（MVIF模型），并进行近实时概率预测分析.结果表明：（1）MVIF模型改善了原模型对计算起始时刻选择苛刻的问题；（2）MVIF模型预测精度较高，在滑坡进入中等加速变形阶段之后可进行可靠的预测；（3）预测滑坡时间与破坏概率结合提供了一种新的滑坡预报准则.该研究可为蠕滑型滑坡的预警预报提供有价值的参考.

关键词

滑坡 / 滑坡时间预测 / Verhulst反函数模型 / 模型改进 / 概率分析 / 预测决策 / 灾害

Key words

landslides / landslide time-of-failure forecast / Verhulst inverse-function model / model modification / probabilistic analysis / forecasting decision / hazards

0 引言

滑坡是我国主要的地质灾害形式之一，每年都会带来严重的生命财产损失（许强，2020）.可靠的滑坡时间预测是滑坡灾害预警和风险管理的关键组成部分.变形是表征滑坡演化最直观的指标，因此当前滑坡时间预测大多采用基于变形观测数据以及蠕变理论的“现象学”预测方法（Intrieri et al.，2019）.自Saito（1965）提出第一个实用的预测模型以来，国内外学者针对滑坡时间预测这一课题开展了广泛的研究，提出了数十种预测模型，其中基于速度监测数据的逆速度法（Fukuzono，1985）是现今最为流行的预测方法.然而，真实滑坡记录的速度监测数据受环境噪声的影响过大，因此采用位移监测数据进行滑坡预测被认为是更优的选择（Crosta and Agliardi，2003）.

Verhulst反函数模型（Verhulst inverse-function model，下文简称“VIF模型”）即为一种基于位移监测数据的滑坡时间预测模型（李天斌和陈明东，1996）.滑坡演化类似于生物孕育、生长、成熟和消亡的演变过程，为此晏同珍（1988）引入Verhulst生物生长模型预测滑坡时间并取得初步成功.李天斌和陈明东（1996）指出Verhulst生物生长模型反函数的反“S”曲线与典型滑坡位移曲线类似，建议采用预测理论和依据更为充分的VIF模型进行滑坡时间预测.当前，VIF模型不仅成为短期或临滑滑坡时间预测的常用工具之一，而且还可用于推断滑坡中长期位移发展趋势（Hu et al.，2021），这是逆速度法等只能用于加速变形阶段的预测模型所不具备的优势.此外，该模型已经被拓展至隧道围岩以及矿山地压灾害的预测预报中（贺小黑，2018）.然而，研究表明VIF模型有时对加速变形阶段位移监测数据的拟合效果很差，计算起始时刻的选择直接影响模型对监测数据的拟合质量，进而导致预测精度过低（龙万学等，2008）.因此，VIF模型的预测理论和适用性还有待进一步的研究.

在实际滑坡时间预测中，由于受滑坡演化的复杂特征、外界触发因素、观测误差、环境噪声以及预测模型假定的局限性等诸多不确定性因素的影响，精确预测滑坡的发生时刻极为困难（Zhang et al.，2020；唐辉明等，2022）.传统的滑坡时间预测基本是基于确定性分析的点估计预测，即计算一个固定、唯一的预测值代表未来最可能发生滑坡的时间，其预测结果往往与真实滑坡时间存在一定差异.因此，这种预测方法虽然计算简便且易于理解，但很难为滑坡预警和风险管理提供更多可靠的预测信息.近年，一些学者开始对滑坡时间预测的可靠性以及预测不确定性问题开展了研究.这些研究的思路大致可分为3类：其一，采用多种预测模型对滑坡时间进行综合分析，利用各模型预测结果的互补特性进行预测决策，从而提升预测的可靠性（任开瑀等，2020）；其二，采用区间估计替代点估计对滑坡可能发生失稳破坏的时间区间进行推断，并制定相应的预警准则（Carlà et al.，2017b）.其三，利用概率分析方法对滑坡时间进行不确定性量化，例如Zhang et al.（2020）采用极大似然估计以及一个滑坡数据库对滑坡时间进行概率预测，并综合考虑了观测不确性和模型不确定性.这些方法虽然在不同程度上量化了滑坡时间预测的不确定性，但是对实时或近实时预测可靠性的研究相对较少.滑坡演化具有阶段性，不同演化阶段的预测不确定性显然存在差异.

针对VIF模型存在的不足，本文提出了一种基于位移监测数据的改进VIF模型.此外，考虑到滑坡演化存在不确定性，本文还建立了一个滑坡案例库对提出模型的预测不确定性进行标定，并实施近实时概率预测分析，进而提出一种考虑预测不确定性的滑坡预报准则.提出的方法可为滑坡时间预测和灾害预警提供有价值的参考.

1 VIF模型及改进

Verhulst生物生长模型的微分方程为：

\frac{\mathrm{d}s}{\mathrm{d}t} = a s - b s^{\mathrm{2}}，(1)

式中，s为位移，t为观测时间，a、b为模型参数.假定计算起始时刻t ₀对应的位移为s ₀，则由上式可得位移s的表达式为（晏同珍，1988）：

s \left(t\right) = \frac{s_{\mathrm{0}}}{b{s}_{\mathrm{0}}+\left(a-s_{\mathrm{0}}\right){\mathrm{e}}^{-a\left(t-t_{\mathrm{0}}\right)}},(2)

由式（2）推导出Verhulst生物生长模型的反函数：

s \left(t\right) = \frac{\mathrm{1}}{a} \mathrm{l} \mathrm{n} \frac{\left(a-b{t}_{\mathrm{0}}\right)t}{\left(a-bt\right)t_{\mathrm{0}}} + s_{\mathrm{0}} , {}_{}{}^{}t t_{\mathrm{0}} \ne \mathrm{0}，(3)

该式即为VIF模型（李天斌和陈明东，1996），其预测曲线为一反“S”曲线，可以反映典型滑坡演化的3个变形阶段，即初始变形、等速变形和加速变形阶段（图1a）.

假定滑坡发生时的位移趋于无穷大，则t → a/b，即预测滑坡时间t_f =a/b.将t_f 代入上式可得：

s \left(t\right) = \frac{\mathrm{1}}{b{t}_f} \mathrm{l} \mathrm{n} \frac{b{t}_f/t_{\mathrm{0}}-bt}{t_f-t} + s_{\mathrm{0}}，(4)

令A取代1/（bt_f ），B取代b和t ₀，则上式可改写成：

s \left(t\right) = A \mathrm{l} \mathrm{n} \frac{t_f-Bt}{t_f-t} + s_{\mathrm{0}}，(5)

该式即为改进的VIF模型（下文简称“MVIF模型”）.式中，s ₀是一个固定的监测位移值，具有一定噪音，直接输入预测模型，可能会限制模型的预测性能.因此，s ₀可替换成一模型参数C，则上式改成：

s \left(t\right) = A \mathrm{l} \mathrm{n} \frac{t_f-Bt}{t_f-t} + C .(6)

MVIF模型的位移预测曲线同样为反“S”曲线（图1b）.与VIF模型相比，MIVF模型中t_f 可由模型参数估计直接获取，无需再经由类似于VIF模型的模型参数转换，简化了计算流程.

2 滑坡时间预测不确定性

Zhang et al.（2020）将滑坡时间预测的不确定性分为两类，其一为由观测误差和环境噪声等因素引起位移监测数据出现离散性的观测不确定性，其二为模型假定造成的模型不确定性.对于MVIF模型的不确定性标定，本文参考Zhang et al.（2020）提出的方法进行分析，同时考虑观测不确定性和模型不确定性两种情况，并采用极大似然估计对滑坡时间进行估计.另外，本文还首次对滑坡加速变形阶段不同时期的模型不确定性进行量化，从而为滑坡实时或近实时预测分析提供参考.

滑坡时间的不确定性表述如下（Zhang et al.，2020）：

T_f = t_f + {\epsilon }_m，(7)

式中，T_f 为真实滑坡时间， ε_m 为T_f 与t_f 二者偏差的随机变量.假定t_f 和ε_m 服从正态分布，即t_f~N（μ_f， {\sigma }_f^{\mathrm{2}}）和ε_m~N（μ_m， {\sigma }_m^{\mathrm{2}}），则T_f 的均值μ_T 和标准差σ_T 分别为：

{\mu }_T = {\mu }_f + {\mu }_m，(8)

{\sigma }_T = \sqrt[]{{\sigma }_f^{\mathrm{2}}+{\sigma }_m^{\mathrm{2}}} .(9)

根据正态分布的假定，滑坡时间的概率密度函数和累积分布函数分别为：

f \left(t\right) = \varphi \left(\frac{t-{\mu }_f-{\mu }_m}{\sqrt[]{{\sigma }_f^{\mathrm{2}}+{\sigma }_m^{\mathrm{2}}}}\right)，(10)

.(11)

2.1　观测不确定性计算方法

受观测环境条件、仪器观测误差、人为误差等因素的影响，获取的位移监测数据往往携带一定的噪声.数据噪声使位移数据出现离散性，不可避免地对预测精度造成影响，因此滑坡预测需量化噪声带来的观测不确定性.观测不确定性采用极大似然估计进行计算.为量化观测不确定性，MVIF模型改写成：

s_{\epsilon } ( t ) = A \mathrm{l} \mathrm{n} \frac{t_f-Bt}{t_f-t} + C + {\epsilon }_{\mathrm{0}}，(12)

式中，ε ₀为一服从正态分布的随机变量，ε ₀ ~N（0， {\sigma }_{\mathrm{0}}^{\mathrm{2}}）.令 θ =｛A， B， C， t_f， ε ₀｝.当 θ 已知时，则时刻t_i 观测到位移值s_i 的概率为：

p \left(s_i\left|\mathit{\theta }\right.\right) = \varphi \left\{\frac{s_i-s\left(t_i\right)}{{\sigma }_{\mathrm{0}}}\right\} .(13)

假定各监测数据点相互独立，则 θ 的似然函数表示为：

l \left(\mathbf{\theta }\left|s_i\right.\right) = {\displaystyle \prod _{i=\mathrm{1}}^N} \varphi \left\{\frac{s_i-s\left(t_i\right)}{{\sigma }_{\mathrm{0}}}\right\}，(14)

依据极大似然估计， θ 的最优值为使该式取值最大时的情况.为了降低算法的收敛难度，该式可通过对数似然函数进行求解：

L \left(\mathit{\theta }\right) = \mathrm{l} \mathrm{n} l \left(\mathit{\theta }\left|s_i\right.\right) .(15)

θ 的概率密度分布可由上式获取的最优均值 θ*和协方差为 C_θ 的多元正态分布函数进行求解，其中 C_θ 表示如下：

{\mathit{C}}_{\mathit{\theta }} = {\left[-\frac{{\partial }^{\mathrm{2}}L\left(\mathit{\theta }\right)}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\mathrm{T}}}\right]}^{-\mathrm{1}} {\left|\right.}_{\mathit{\theta }={\mathit{\theta }}^{\mathrm{*}}} .(16)

2.2　模型不确定性计算方法

模型不确定性为预测模型的固有偏差，反映模型假定偏离真实滑坡演化特征的程度.模型不确定性可由ε_m 的均值μ_m 和标准差σ_m 进行计算（Zhang et al.，2020）.假定存在n个滑坡案例，T_fj 为第j案例的真实滑坡时间， T =｛T_{f 1}， T_{f 2}， …， T_fn ｝，则由极大似然估计计算的第j个案例t_f 的均值和标准差分别为μ_fj 、σ_fj .当参数 θ' =｛μ_m， σ_m ｝已知，观测到N_f 个滑坡案例滑坡时间的似然函数表示为：

l \left(\mathit{\theta }\left|\mathit{T}\right.\right) = {\displaystyle \prod _{i=\mathrm{1}}^n} \varphi \left(\frac{T_{fj}-{\mu }_{fj}-{\mu }_m}{\sqrt[]{{\sigma }_{fj}^{\mathrm{2}}+{\sigma }_m^{\mathrm{2}}}}\right)，(17)

该式同样可采用极大似然估计获得 θ'的最优估计值.需要指出的是，通常认为加速变形行为是滑坡发生失稳破坏的前兆，因此本文只对加速变形阶段的模型不确定性进行标定，不考虑初始变形和等速变形阶段的情况.

当前的研究基本只采用滑坡之前最后一次观测的数据对模型不确定性进行标定（Zhang et al.，2020）.实际上，由于滑坡演化机制的复杂性，不同演化阶段的模型不确定性显然会有所区别.这里引入改进切线角α作为划分加速变形阶段不同时期的量化尺度.α的计算方法详见许强等（2015）的研究.依据α的定义，当45° < α < 80°时，滑坡处于初始加速变形阶段；当80° ≤ α < 85°时，滑坡处于中加速变形阶段；当α ≥ 85°时，滑坡处于加加速变形阶段；而当α趋于90°时，意味着滑坡即将发生失稳破坏.研究表明，当α ≥ 60°时，滑坡的加速变形趋势趋于明显且α的波动性降低（方汕澳等，2021）.为此，本研究以α=60°作为计算起点，并5°为计算步长，分别量化滑坡在加速变形阶段α ≤ 60°、< 65°、< 70°、< 75°、< 80°、< 85°和 < 90°七个时期的模型不确定性.考虑到原始监测数据存在噪声，加速变形阶段的位移监测数据采用MVIF模型预先校准和平滑处理后再求α.

3 研究结果

3.1　预测不确定性标定结果

为量化预测不确定性，本文共搜集了20个具有相对完整位移时间序列和给出真实滑坡时间T_f 的滑坡案例（表1）.首先，采用图像数字化识别工具对文献中滑坡案例的监测数据进行数字化识别.其次，按5°改进切线角的步长计算滑坡的观测不确定性，采用极大似然估计计算每个滑坡在加速变形阶段不同时期滑坡时间预测的均值μ_f 和标准差σ_f，其中μ_f 也作为滑坡失稳破坏的最可能发生时间.最后，基于观测不确定性的结果，计算加速变形阶段不同时期的模型不确定性.需要指出的是，20个滑坡案例破坏时α均超过85°，且在划分的七个变形时期均具有对应的位移数据点.

图2展示了各个案例在加速变形阶段不同时期预测的μ_f，其中虚线为各案例的真实滑坡时间T_f .由图可知，μ_f 和T_f 之间存在一定误差，且误差存在不确定性.总体而言，当α < 70°时，预测误差偏大；当α > 75°时，预测误差减小；而当α > 80°时，预测的μ_f 基本在T_f 附近波动.这一结果表明，随着滑坡变形加剧，滑坡失稳破坏趋势趋于明显，滑坡时间的可预测性逐步提升.另外，值得关注的是对于20个案例的位移监测数据，VIF模型在加速变形阶段不同时期的预测均具有较高的拟合质量.

表1汇总了各个案例最后一次观测的预测滑坡时间均值μ_f 和标准差σ_f，并计算了各案例仅考虑观测不确定性的预测滑坡时间的95%置信区间（CI）.由表可知，20个滑坡案例最后一次观测获取的t_f 与T_f 之间存在0.01~3 d的误差，其中6个案例的预测误差绝对值大于1 d，且只有2个案例的T_f 位于95%CI之内.这一结果表明，滑坡时间的预测误差除了受监测数据噪声的影响，还受到模型假定等因素的影响.

图3为MVIF模型在加速变形阶段不同时期模型不确定性｛μ_m， σ_m ｝的量化结果，由式（17）计算得到，未叠加观测不确定性的均值和标准差.总体趋势上，随着α增大，模型不确定性降低.当α < 70°时，虽然模型不确定性的μ_m 为3.53~4.55 d，但σ_m 已经大于19 d，表明在初始加速变形阶段的前期，滑坡失稳破坏的趋势还未明晰，难以可靠地预测滑坡时间.滑坡进入中加速变形阶段（α ≥ 80°）之后，模型不确定性的μ_m 和σ_m 均比初始加速变形阶段的要低.在加加速变形阶段，模型预测不确定性的μ_m 和σ_m 分别为0.29 d和0.95 d，此时滑坡时间的可预测性较高.表1罗列了考虑模型不确定性最后一次观测的预测结果.这里的“考虑模型不确定性”是指同时考虑观测不确定性和模型不确定性，下文同.预测结果显示，考虑模型不确定性后，共有18个案例的T_f 位于95%CI之内；只有Mt. Beni滑坡和New Tredegar滑坡2个案例的T_f 位于95%CI之外，但接近于95%CI.因此，在滑坡时间预测中，必须同时考虑观测不确定性和模型不确定性.

3.2　滑坡案例分析

案例1：陕西省周至 G108 路段滑坡.陕西省周至G108路段滑坡（下文简称“周至滑坡”）位于陕西省周至县的南部山区，滑坡区域内地质构造强烈，发育多条张拉裂缝，岩性主要为含炭绢云母石英片岩（喻小等，2019）.滑坡的长期变形对国道G108的路面和路肩挡土墙造成了破坏.2017 年 10月，滑坡因降雨作用而发生局部滑塌.2017年12月15日，相关单位在国道的上边坡和下边坡路基共布置了4个GNSS监测点并开展监测活动.2018年2月19日凌晨04： 30左右，滑坡的左侧区域发生滑塌.根据喻小等（2019）提供的原始监测数据，本案例主要以滑塌区域监测点JC02对提出的预测方法进行分析（图4）.起始观测时间定义为t ₁=1 d，则T_f 约为67.19 d.

MVIF模型和VIF模型利用最后一次观测数据的预测曲线见图4.以不同计算起始时刻获取MVIF模型和VIF模型的模型参数见表2.MVIF模型选取不同计算起始时刻进行预测均表现出高质量的拟合效果（拟合优度R ² > 0.99）.由表2 可知，本例中以计算起始时刻1（MVIF-1）和计算起始时刻2（MVIF-2）为计算起点获取的MVIF模型的模型参数较为接近，表明该滑坡在加速阶段变形趋势相对平稳.MVIF-1和MVIF-2的t_f 均约为68.05 d，与真实滑坡时间T_f 相差约1 d.对于VIF模型，以计算起始时刻1（VIF-1）为计算起始点的拟合效果较差，预测曲线没有沿着原始监测数据的趋势发展，无法进行有效的预测；以计算起始时刻2（VIF-2）为计算起始点的拟合效果较好（R ²=0.992 8），t_f 为67.14 d，较为准确地预测滑坡时间.由分析可知，虽然VIF-2的预测结果最接近于T_f，但VIF模型计算起始时刻选择不当会严重影响预测质量，而MVIF模型可有效地改善这一问题，并且具有较高的预测精度.

图5为仅考虑观测不确定性以及考虑模型不确定性的近实时预测结果.模型不确定性采用图3的结果进行标定，依据α，加速变形阶段不同时期采用对应的模型不确定性.由图5可知，当α < 60°时，预测结果波动性较大，难以可靠地预测滑坡时间；当α > 70°时，t_f 基本在T_f 上下波动，提出的模型可较为可靠地预测滑坡时间.

值得关注的是，对于仅考虑观测不确定性的情况，除了极个别预测点，预测滑坡时间的95%CI几乎没有将T_f 容纳其中；对于考虑模型不确定性的情况，除了加速初期的个别预测点，预测滑坡时间的95%CI基本将T_f 容纳其中.造成这一结果的主要原因是该滑坡加速变形演化趋势的变异性低且MVIF模型拟合精度极高，使得滑坡观测不确定性远小于模型不确定性.以最后一次观测为例，仅考虑观测不确定性时，预测滑坡时间的σ_f 为0.043 d，而这一时期模型不确定性的σ_m 为0.95 d.图6为最后一次观测仅考虑观测不确定性以及考虑模型不确定性两种情况的预测滑坡时间概率密度函数和累积分布函数.两种情况的函数曲线差异明显，仅考虑观测不确定性的曲线形态极为“高瘦”，T_f 位于95%CI之外；而考虑模型不确定性的函数曲线则相对“矮胖”，预测滑坡时间的95%CI为［66.78 d， 69.91 d］，将T_f 容纳其中.显然，若只考虑观测不确定性，滑坡时间预测的不确定性将被低估，考虑模型不确定性之后更符合实际预测结果.

案例2：茂县滑坡.茂县滑坡位于四川省阿坝州茂县叠溪镇新磨村背部的富贵山，该区域新构造运动活跃，松坪沟断裂带从其间穿插而过，地质灾害频发（闫斌，2021）.滑坡体主要由变质砂岩、大理化灰岩以及千枚岩组成.据调查，滑源区在2003年已出现多条明显的地表裂缝.滑坡发生前，滑坡区域出现持续性降雨，累积降雨量约80 mm.持续性降雨是滑坡失稳破坏的直接诱发原因.2017年6月24日凌晨，富贵山突发特大高位滑坡事件并形成碎屑流，滑坡体积约450万m³，水平滑动距离超过 2 500 m.滑坡体摧毁了山下的新磨村，造成83人死亡.滑坡发生之后，很多学者采用干涉合成孔径雷达技术（InSAR）对滑坡地表历史变形进行反演（任开瑀等，2020；闫斌，2021）.本案例选取闫斌（2021）提供的InSAR位移监测数据（监测点P2）进行预测分析.起始观测时间t ₁=1 d，真实滑坡时间T_f 约为426 d.

MVIF模型和VIF模型利用最后一次观测数据的预测曲线见图7.以不同计算起始时刻获取MVIF模型和VIF模型的模型参数见表3.与周至滑坡类似，MVIF模型以不同时刻作为计算起点进行预测均表现出良好的拟合质量（R ² > 0.99），以计算起始时刻1和2为计算起点获取的t_f 分别约为427.03 d和 425.54 d，接近于T_f .在VIF模型的预测中，同样出现计算起始时刻选择不恰当（计算起始时刻1）导致拟合质量差的情况，从而无法有效地进行滑坡时间预测.VIF模型以计算起始时刻2为计算起点可获得较为理想的拟合结果（R ²=0.9951），预测的t_f 为416.50 d.

图8为仅考虑观测不确定性以及考虑模型不确定性的近实时预测结果，预测的模型不确定性同样采用图3的结果进行标定.由图8可知，当α < 70°时，t_f 与T_f 的误差偏大，难以可靠地预测滑坡时间；当α > 70°时，t_f 开始逐步向真实滑坡时间收敛；当α > 75°时，t_f 接近于T_f .

对比仅考虑观测不确定性以及考虑模型不确定性两种情况，可发现：考虑模型不确定性的情况在α > 75°之后，预测结果的95%CI将T_f 容纳其中；而仅考虑观测不确定性的情况需在α > 80°之后，预测结果的95%CI才能包含T_f .以最后一次观测为例（图9），仅考虑观测不确定性时，预测滑坡时间的μ_f 和σ_f 分别为427.03 d和1.14 d，预测滑坡时间的95%CI为［425.15 d， 428.90 d］；对于考虑模型不确定性的情况，预测滑坡时间的95%CI为［424.87 d， 429.76 d］.两种情况的概率分布曲线较为接近，这是由于σ_f 略大于σ_m 造成的，即观测不确定性对预测不确定性的影响大于模型不确定性.

4 讨 论

4.1　MVIF模型的适用性

上述分析表明，MVIF模型具有较高的预测精度，并且有效地改善了VIF模型对计算起始时刻选择苛刻的问题，因而具有一定的实时或近实时预测能力，但在使用时仍应重视其适用性和局限性.

（1）MVIF模型更适用于蠕滑型滑坡，对于脆性破坏型滑坡的适用性仍需进一步探讨.通常，脆性破坏型滑坡前期变形很小或可忽略不计，从出现加速至最终破坏的持续时间非常短暂，具有一定突发性，因此预测此类滑坡的失稳破坏时间更为困难.Carlà et al.（2017a）指出高频的监测数据可以显著地提升此类滑坡的可预测性.这里以Carlà et al.（2017a）提供的两个露天矿脆性破坏案例Failure#1和Failure#3来研究MVIF模型在脆性破坏型滑坡预测的可行性（图10）.由图可知，MVIF模型对两个案例的位移监测数据均具有较高的拟合质量，预测的滑坡时间也接近于T_f .虽然如此，但由于急剧发展的脆性破坏行为留给预测活动的时间过短，因此建议MVIF模型在脆性破坏预测中不宜单独使用或不使用.

（2）MVIF模型完整的预测曲线具有三阶段蠕变曲线特征，因此具有预测滑坡中长期位移发展趋势以及识别滑坡演化阶段的潜力.以呷爬滑坡为例（Hu et al.，2021），该滑坡位于锦屏I级水电站上游11.5 km处， 受库水调度以及降雨等因素的影响而发生缓慢变形（图11）.由图可知，MVIF模型展现出了对该滑坡位移监测数据的良好拟合效果（R ²=0.992 8）.依据预测曲线的发展趋势，滑坡在未来一段时间内仍处于匀速变形阶段.另外，从图中不难发现MVIF模型更多的是对滑坡总体位移趋势的预测，难以精细刻画如出现多期加速变形行为的复杂位移模式.在这一方面，MVIF模型的位移预测能力逊于现今较为流行的机器学习方法.机器学习方法以一种“黑匣子”的方式构建滑坡位移与外部触发因素之间的关系，使其可以对滑坡长期演化的复杂位移模式进行良好的预测，但它们中的大多数预测方法很难解释滑坡的变形机制以及缺乏一定的泛化能力（Intrieri et al.，2019；郭子正等，2022）.Hu et al.（2021）将VIF模型与随机森林方法结合，建立了一个能够充分考虑滑坡运动学和外部触发因素的综合滑坡运动预测模型，提升了机器学习方法预测滑坡位移的泛化能力.

（3）MVIF模型在预测时可降低对加速变形起始点（OOA）识别的要求.通常认为，加速变形行为是滑坡发生破坏的前兆，因此大多数预测模型是基于加速蠕变理论建立的，如逆速度法 （Fukuzono，1985）.这些预测模型大多数只能用于加速变形阶段，因此在使用时必须识别OOA，但在实时预测中准确、及时地识别OOA是一项具有挑战性的工作（Carlà et al.，2017b）.在前文茂县滑坡的分析中，MVIF模型选取了匀速变形阶段和加速变形阶段两个阶段的位移监测数据进行分析，依然获得较好的预测结果（图7）.

（4）MVIF模型可利用非位移监测数据进行预测分析.除了传统的位移监测外，近年一些新的监测技术和手段被引入滑坡监测之中，从而衍生出多种观测变量，如微震信号、斜坡变形倾角等.图12a、12b分别为采用微震监测（Mesnil-Val崩塌）以及采用斜坡变形倾角监测（黑方台滑坡）的两个失稳破坏案例（Amitrano et al.，2005；Wang et al.，2022）.图中的预测结果验证了MVIF模型利用非位移监测数据进行预测的可行性.广义上，MVIF模型可利用那些在临滑前表现出加速变形特征的观测变量进行预测分析，但其预测精度仍需要更深入的研究.

（5）滑坡体材料具有空间变异性，同一滑坡不同滑动区域的变形特征也可能存在差异性，因此采用MVIF模型进行滑坡时间预测时，应选用最能代表滑坡整体变形特征的监测点或者综合多源监测点进行分析.

（6）受文章篇幅限制，当前用于标定模型不确定性的案例相对较少，模型不确定性的大小可能在一定程度上受到选用案例的影响.同时，这些案例的监测数据大多采用图像数字化识别软件从文献图片中提出，因而不可避免地出现一定偏差.因此，需要一个案例数量更多、滑坡类型更完备并具有原始监测数据的大型数据库对MVIF模型的预测能力进行下一步的验证.

4.2　基于MVIF模型的滑坡预报准则

模型不确定性的标定与两个滑坡案例的分析表明，在初始加速变形阶段，t_f 与T_f 差异较大，但在改进切线角α > 70°后开始向真实滑坡时间收敛；在滑坡进入中等加速变形阶段（α ≥ 80°）之后，预测模型可提供较为可靠的预测结果，t_f 基本接近于T_f .因此，为充分利用预测信息以及增强预测的可靠性，建议依据滑坡的不同演化阶段，对MVIF模型的预测结果进行分阶段决策.阶段划分的判据采用改进切线角α和稳定性系数K，划分标准参考许强等（2015）的研究.分段决策的示意图及对应的判据具体见图13，具体描述如下：

（1）滑坡处于初始变形和匀速变形阶段（α ≤ 45°，1.0 ≤ K < 1.15），变形发展相对缓慢，短期内发生失稳破坏的迹象不明显，准确预测滑坡时间的可能性较低，因此在这两个阶段应侧重于识别滑坡的演化阶段，利用MVIF模型预测滑坡的未来变形发展趋势，并及早地发现可能出现的加速变形事件.由于受降雨、库水涨落、人类活动等因素的影响，这两个变形阶段可能会出现多次加速变形事件.因此，对于识别出的加速变形事件，必须依据滑坡长期演化趋势特征以及K的大小，综合判断滑坡是否已经进入最终的加速变形阶段.

（2）在滑坡进入初始加速变形阶段之后（45° < α < 80°，0.95 ≤ K < 1.0），变形速度开始逐渐增加，可以开始开展滑坡时间预测工作.然而，这一阶段的前期距离真实滑坡时间相对较长，失稳破坏趋势仍较为模糊，微小的扰动也可能会造成较大的预测误差（亓星等，2020）.因此，在初始加速变形阶段前期，建议侧重于变形破坏趋势的预测，判断滑坡是否存在短期内发生失稳破坏的可能性.当α > 70°之后，可以逐步开展具体的滑坡时间预测工作.

（3）在滑坡进入中等加速变形阶段之后（80° ≤ α < 85°，K < 0.95），变形速度出现较快的增长趋势，失稳破坏趋势趋于明显，这一阶段可进行较为可靠的滑坡时间预测.

（4）在滑坡进入加加速变形阶段之后（α ≥ 85°，K < 0.95），变形骤然增加且对外界的扰动异常敏感，位移曲线与时间坐标轴近似垂直，任何微小的扰动都可能触发滑坡发生，此时可对滑坡时间进行确定性预测预报.

据前文分析可知，由于滑坡演化复杂多变以及监测数据存在误差，即使滑坡进入加加速变形阶段，t_f 与T_f 之间仍存在一定偏差，因此引入预测不确定性进行预测决策是极为关键的.预测的均值和概率区间可以为评估滑坡最可能发生失稳破坏的时刻和时间区间提供有价值的参考.此外，这里提供一种将预测滑坡时间t_f 与破坏概率P_f 结合的滑坡预测预报决策方法：

，(18)

式中，t′=t_p + ∆t，t_p 为当前观测时刻，∆t为某一假定可能发生滑坡时间 T_f^{\text{'}}的前∆t d.该式可表述为未来∆t d内，滑坡发生失稳破坏的概率为Φ（∙）.∆t可依据滑坡的实际情况和决策者的风险态度进行选取.根据式（18）以及滑坡加速变形阶段的不同时期，本研究制定了如下滑坡时间预报准则：

\left\{\begin{array}{c}{P}_f\ge \mathrm{50}\mathrm{\%}{,}_{}{}_{}\alpha >\mathrm{70}°\mathrm{且}\mathrm{\Delta }t=\mathrm{15}{}_{}{}^{}\mathrm{d},\\ P_f\ge \mathrm{30}\mathrm{\%}{,}_{}{}_{}\alpha \ge \mathrm{80}°\mathrm{且}\mathrm{\Delta }t=\mathrm{10}{}_{}{}^{}\mathrm{d},\\ P_f\ge \mathrm{5}\mathrm{\%}{,}_{}{}_{}\alpha \ge \mathrm{85}°\mathrm{且}\mathrm{\Delta }t=\mathrm{5}{}_{}{}^{}\mathrm{d},\end{array}\right. (19)

该式的含义为当改进切线角α>α′（如70°）时，未来∆t d（如15 d）内滑坡破坏概率P_f ≥ P_f^{\text{'}}（如50%）时可发出滑坡时间预报.下文利用前文的周至滑坡和茂县滑坡对提出预报准则的可行性进行验证，计算结果见图14.

图14a为周至滑坡的滑坡时间预报结果.当α > 75°时，该滑坡未来15 d发生滑坡的概率为57.65%，此后滑坡发生破坏的概率均不低于50%.根据式（19），意味着α > 75°之后可以发出滑坡在未来15 d内可能发生失稳破坏的预报，对照真实滑坡时间，可提前9.2 d进行预报.

图14b为茂县滑坡的滑坡时间预报结果.由图可知，只有α > 85°之后，即滑坡进入加加速之后才能发出预报.造成这一结果的主要原因有2个：其一，与该滑坡的成因机制和变形特征相关，从而影响滑坡的可预测性；其二，该滑坡的位移采用InSAR技术进行观测，数据的采样频率以及精度均低于常规的地面监测技术和手段，因此可能会导致预测误差偏大和预报延迟.虽然发出预报的时间节点晚于周至滑坡，但对照真实滑坡时间，依然可提前5 d发出预报.

上述两个案例验证了提出预报准则的可行性，但P_f 和∆t的阈值标准仍需要更多的案例进行标定和修正，同时还要考虑工程实例的具体情况，如成因机制、变形特征、演化阶段、监测技术和数据处理方法等.此外，考虑到滑坡演化的复杂性，可靠的滑坡时间预测预报还需要结合多种预测方法、多种观测变量（如地下水位、斜坡变形倾角和微震信号）以及融合多源监测数据进行综合分析和时空预测.提出的模型和方法如何有机地结合这些要素进行滑坡预测预报还有待于进一步研究.

5 结论

本文针对Verhulst反函数滑坡预测模型（VIF模型）对计算起始时刻选择不当会造成监测数据拟合质量差及预测精度低的问题，提出了改进的VIF模型（MVIF模型），并进行近实时概率预测分析以及构建预测预报决策方法.主要得出以下结论：

（1）提出的MVIF模型与VIF模型类似，预测曲线同样为反“S”曲线，与典型滑坡三阶段位移演化曲线类似，适用于蠕滑型滑坡.MVIF模型改善了VIF模型对计算起始时刻选择苛刻的问题，并且具有高质量的拟合效果和较高的预测精度.

（2）利用20个滑坡案例对MVIF模型的预测不确定性进行标定，预测不确定性分为观测不确定性和模型不确定性两类.由于模型拟合精度较高，观测不确定性相对较小.模型不确定性的标定结果显示，加速变形阶段不同时期的模型不确定性存在差异.整体趋势上，模型不确定性的均值和标准差随滑坡变形发展而逐步减少.在滑坡进入中等加速变形阶段后，MVIF模型可进行可靠的滑坡时间预测.在滑坡进入加加速变形阶段后，模型不确定性的均值和标准差分别为0.29 d和0.95 d，此时预测不确定性较小，可进行确定性预测预报.

（3）预测滑坡时间与破坏概率结合提供了一种新的滑坡预报量化准则.案例分析验证了该方法的可行性，为滑坡灾害的预警预报提供了有价值的参考.但是考虑到滑坡演化的复杂性，建议提出的方法不宜单独使用，且破坏概率的阈值标准仍需更多的案例进行校准.

参考文献

[1] Amitrano, D., Grasso, J. R., Senfaute, G., 2005. Seismic Precursory Patterns before a Cliff Collapse and Critical Point Phenomena. Geophysical Research Letters, 32(8): L08314. https://doi.org/10.1029/2004gl022270

[2] Carey, J.M., Moore, R., Petley, D., et al., 2007. Pre- Failure Behaviour of Slope Materials and Their Significance in the Progressive Failure of Landslides. Landslides and Climate Change: Challenges and Solutions, London, 207-215.

[3] Carlà, T., Farina, P., Intrieri, E., et al., 2017a. On the Monitoring and Early-Warning of Brittle Slope Failures in Hard Rock Masses: Examples from an Open-Pit Mine. Engineering Geology, 228: 71-81. https://doi.org/10.1016/j.enggeo.2017.08.007

[4] Carlà, T., Intrieri, E., Di Traglia, F., et al., 2017b. Guidelines on the Use of Inverse Velocity Method as a Tool for Setting Alarm Thresholds and Forecasting Landslides and Structure Collapses. Landslides, 14(2): 517-534. https://doi.org/10.1007/s10346-016-0731-5

[5] Carlà, T., Farina, P., Intrieri, E., et al., 2018. Integration of Ground-Based Radar and Satellite InSAR Data for the Analysis of an Unexpected Slope Failure in an Open-Pit Mine. Engineering Geology, 235: 39-52. https://doi.org/10.1016/j.enggeo.2018.01.021

[6] Chen, H., Tang, H., Ge, X.R., et al., 2019. Study on Early Warning Index and Early Warning Forecast of Creep Landslide Based on Deep Displacement. Chinese Journal of Rock Mechanics and Engineering, 38(S1): 3015-3024 (in Chinese with English abstract).

[7] Crosta, G. B., Agliardi, F., 2003. Failure Forecast for Large Rock Slides by Surface Displacement Measurements. Canadian Geotechnical Journal, 40(1): 176-191. https://doi.org/10.1139/t02-085

[8] Fang, S.A., Xu, Q., Xiu, D.H., et al., 2021. A Study of the Predicted Instability Time of Sudden Loess Landslides Based on the SLO Model. Hydrogeology & Engineering Geology, 48(4): 169-179 (in Chinese with English abstract).

[9] Fukuzono, T., 1985. A New Method for Predicting the Failure Time of a Slope. The Fourth International Conference and Field Workshop on Landslides, Tokyo, 145-150.

[10] Gigli, G., Fanti, R., Canuti, P., et al., 2011. Integration of Advanced Monitoring and Numerical Modeling Techniques for the Complete Risk Scenario Analysis of Rockslides: The Case of Mt. Beni (Florence, Italy). Engineering Geology, 120(1/2/3/4): 48-59. https://doi.org/10.1016/j.enggeo.2011.03.017

[11] Guo, Z.Z., Yang, Y.F., He, J., et al., 2022. A Novel Deep Learning Model to Predict Landslide Displacement Considering Attention Mechanism. Earth Science, 1-21 (in Chinese with English abstract).

[12] Harries, N., Noon, D., Rowley, K., 2006. Case Studies of Slope Stability Radar Used in Open Cut Mines. In: Stability of Rock Slopes in Open Pit Mining and Civil Engineering Situations, Proceedings International Symposium on Stability of Rock Slopes in Open Pit Mining and Civil Engineering. The South African Institute of Mining and Metallurgy, Johannesburg, 335-342.

[13] He, M.C., Han, X., Zhang, B., et al., 2012. Real-Time Remote Monitoring and Forecasting Technology for Landslide Disasters Based on Sliding Force Variation and Its Engineering Application. Journal of Heilongjiang Institute of Science and Technology, 22(4): 337-342 (in Chinese with English abstract).

[14] He, X.H., 2018. Modification of Verhulst Inverse Function Model of Landslide Forecast. Geology and Mineral Resources of South China, 34(4): 289-293 (in Chinese with English abstract).

[15] Hu, X. L., Wu, S. S., Zhang, G. C., et al., 2021. Landslide Displacement Prediction Using Kinematics-Based Random Forests Method: A Case Study in Jinping Reservoir Area, China. Engineering Geology, 283: 105975. https://doi.org/10.1016/j.enggeo.2020.105975

[16] Intrieri, E., Carlà, T., Gigli, G., 2019. Forecasting the Time of Failure of Landslides at Slope-Scale: A Literature Review. Earth-Science Reviews, 193: 333-349. https://doi.org/10.1016/j.earscirev.2019.03.019

[17] Kieffer, D. S., Valentin, G., Unterberger, K., 2016. Continuous Real-Time Slope Monitoring of the Ingelsberg in Bad Hofgastein, Austria. Geomechanics and Tunnelling, 9(1): 37-44. https://doi.org/10.1002/geot.201500047

[18] Leng, C.Q., Zhang, Y., 2016. Prediction of Slope Instability Time in Jinlonggou Yard Based on Deformation Monitoring. Design of Hydroelectric Power Station, 32(1): 84-88 (in Chinese).

[19] Li, T.B., Chen, M.D., 1996. Felhaas Inverse Function Model Method for Landslide Time Prediction. Journal of Geological Hazards and Environment Preservation, 7(3):13-17, 29 (in Chinese with English abstract).

[20] Loew, S., Gschwind, S., Gischig, V., et al., 2017. Monitoring and Early Warning of the 2012 Preonzo Catastrophic Rockslope Failure. Landslides, 14(1): 141-154. https://doi.org/10.1007/s10346-016-0701-y

[21] Long, W.X., Lin, J., Xu, X.H., et al., 2008. Selection of Initial Prediction Time of Landslide Based on Verhulst Inverse Function Model. Chinese Journal of Rock Mechanics and Engineering, 27(S1): 3298-3304 (in Chinese with English abstract).

[22] Naismith, W.A., Wessels, S.D.N., 2005. Management of a Major Slope Failure at Nchanga Open Pit, Chingola, Zambia. Journal of the Southern African Institute of Mining and Metallurgy, 105(9):619-626. https://hdl.handle.net/10520/AJA0038223X_3086

[23] Qi, X., Zhu, X., Xu, Q., et al., 2020. Improvement and Application of Landslide Proximity Time Prediction Method Based on Saito Model. Journal of Engineering Geology, 28(4): 832-839 (in Chinese with English abstract).

[24] Ren, K.Y., Yao, X., Zhao, X.M., et al., 2020. Study on Landslide Instability and Failure Prediction Based on Time Series InSAR, GPS and Image Migration Measurement. Chinese Journal of Rock Mechanics and Engineering, 39(S2): 3421-3431 (in Chinese with English abstract).

[25] Saito, M., 1965. Forecasting the Time of Occurrence of a Slope Failure. In: Proceedings of the 6th International Conference on Soil Mechanics and Foundation Engineering. Montreal, 2: 537-541.

[26] Saito, M., 1979. Evidential Study on Forecasting Occurrence of Slope Failure. OYO Technical Report, 1:1-23.

[27] Seguí, C., Rattez, H., Veveakis, M., 2020. On the Stability of Deep-Seated Landslides. The Cases of Vaiont (Italy) and Shuping (Three Gorges Dam, China). Journal of Geophysical Research: Earth Surface, 125(7): e2019JF005203. https://doi.org/10.1029/2019JF005203

[28] Suwa, H., Mizuno, T., Ishii, T., 2010. Prediction of a Landslide and Analysis of Slide Motion with Reference to the 2004 Ohto Slide in Nara, Japan. Geomorphology, 124(3/4): 157-163. https://doi.org/10.1016/j.geomorph.2010.05.003

[29] Tang, H.M., Li, C.D., Hu, W., et al., 2022. What is the Physical Mechanism of Major Landslides?. Earth Science, 47(10): 3902-3903 (in Chinese with English abstract).

[30] Wang, D., Du, H., Wang, Q.L., 2020. Research and Application for Early Warning of Landslides Based on Hierarchical Clustering Coupling Weighting Markov Chain Model. Journal of China Coal Society, 45(5): 1783-1794 (in Chinese with English abstract).

[31] Wang, H. H., Zhong, P., Xiu, D. H., et al., 2022. Monitoring Tilting Angle of the Slope Surface to Predict Loess Fall Landslide: An on-Site Evidence from Heifangtai Loess Fall Landslide in Gansu Province, China. Landslides, 19(3): 719-729. https://doi.org/10.1007/s10346-021-01727-0

[32] Wang, Z.W., Wang, L., Han, Q.Q., et al., 2020. Data Decoding Method of the Displacement Sensor and Its Application in Landslide Monitoring. Journal of Geodesy and Geodynamics, 40(4): 436-440 (in Chinese with English abstract).

[33] Xu, Q., 2020. Understanding the Landslide Monitoring and Early Warning: Consideration to Practical Issues. Journal of Engineering Geology, 28(2): 360-374 (in Chinese with English abstract).

[34] Xu, Q., Tang, M.G., Huang, R.Q., et al., 2015. Monitoring, Early Warning, and Emergency Disposal of Large-Scale Landslides. Science Press, Beijing (in Chinese).

[35] Xu, Q., Zeng, Y.P., 2009. Research on Acceleration Variation Characteristics of Creep Landslide and Early- Warning Prediction Indicator of Critical Sliding. Journal of Rock Mechanics and Engineering, 28(6): 1099-1106 (in Chinese with English abstract).

[36] Yamaguchi, U., Shimotani, T., 1986. A Case Study of Slope Failure in a Limestone Quarry. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 23(1): 95-104. https://doi.org/10.1016/0148-9062(86)91670-0

[37] Yan, B., 2021. MT-InSAR Technology Combined with Inverse-Velocity Method for High Level Landslide Monitoring and Early Warning. Southwest Jiaotong University, Chengdu (in Chinese with English abstract).

[38] Yan, T.Z., 1998. Landslide Dynamic Law and Forecasting Application. The 3th National Congress of Engineering Geology, Chengdu, 37-43 (in Chinese).

[39] Yu, X., Zhao, Q.H., Zhang, C.H., et al., 2019. Application of GNSS Real-Time Monitoring in Landslide Early Warning: Case of Landslide of G108 Section in Zhouzhi County, Shannxi Province. Yangtze River, 50(10): 126-130, 142 (in Chinese with English abstract).

[40] Zhang, J., Wang, Z. P., Zhang, G. D., et al., 2020. Probabilistic Prediction of Slope Failure Time. Engineering Geology, 271: 105586. https://doi.org/10.1016/j.enggeo.2020.105586

基金资助

广西科技基地和人才专项项目(AD23026111)

广西自然科学基金项目(2024GXNSFBA010226)

广西科技重大专项(AA23023016)