【论文阅读38】-结合应力预测位移

这篇论文提出了一种物理信息约束数据同化方法（PIDA，Physics-Informed Data Assimilation），用于预测受水动力压力驱动滑坡的变形位移。

论文将物理机制模型和数据驱动模型相结合，再通过粒子滤波（PF）算法进行动态参数更新和预测轨迹修正，从而提高滑坡位移预测的精度与可靠性。

1. 提出滑带劣化指数（SZD，Sliding Zone Deterioration），通过深部位移监测数据计算，刻画滑坡滑带的结构完整性与受外界条件驱动下的响应能力。

2. 基于 SZD 指数，动态更新滑带剪切强度参数（黏聚力、内摩擦角），使预测模型能随滑带劣化情况实时修正参数。

3. 构建物理机制模型：

   • Boussinesq 渗流方程求解滑坡潜水面位置；
   • 残余推力法+分块法分析滑坡内力；
   • 建立推力-位移关系模型作为观测方程。

4. 通过粒子滤波算法将 SZD 指数作为状态变量，将监测数据与物理模型预测值对比，实时反馈修正模型参数，动态更新预测轨迹。

• PIDA 方法预测精度优于 LSTM 和 BPNN，尤其在滑坡位移突变点预测效果显著提升。
• PIDA 模型兼具物理机理含义和数据驱动特性，能动态估计预测模型状态变量及 SZD 指数，减少模型结构误差和参数不确定性。
• 对粒子滤波中的粒子退化问题，论文指出了多次重采样可能降低粒子多样性，建议后续改进粒子滤波算法提升模型自适应性。
• 论文强调未来提高监测频率或插值补全数据，有助于增强模型时间精度。

作者

Yong Liu^a, Jingjing Long^b,c, Changdong Li^b,c,∗, Weiwen Zhan^a**

• ^a School of Mechanical Engineering and Electronic Information, China University of Geosciences, Wuhan 430074, China
  中国地质大学（武汉）机械与电子信息学院，武汉 430074

• ^b Faculty of Engineering, China University of Geosciences, Wuhan 430074, China
  中国地质大学（武汉）工程学院，武汉 430074

• ^c Badong National Observation and Research Station of Geohazards, China University of Geosciences, Wuhan 430074, China
  中国地质大学（武汉）地质灾害巴东国家野外科学观测研究站，武汉 430074

关键词：水动力压力驱动型滑坡；预测模型；动力响应机制；滑带劣化；粒子滤波；物理约束数据同化方法

摘要

由于滑坡演化过程复杂且具有非线性特征，其变形规律始终难以预测。针对滑坡物理机制的动态演化过程，本文首次提出了一种基于物理约束的数据同化（Physics-Informed Data Assimilation, PIDA）方法用于滑坡位移预测。该方法能够增强物理预测模型的预测能力，并有效解决传统数学预测模型中忽略滑坡物理意义的问题。此外，粒子滤波（Particle Filter, PF）在非线性、非高斯系统中的优越性也决定了其在复杂问题中的应用价值。
针对水动力压力驱动型滑坡，本文构建了滑带劣化（Sliding Zone Deterioration, SZD）模型作为PF算法的状态方程，同时通过滑坡变形与监测数据之间的关系构建观测方程。基于物理约束的滑坡位移预测模型作为关键参数反馈、修正与更新的桥梁，由剪切强度参数动态更新模型和物理机制分析模型组成。
将该方法与传统的长短期记忆网络（Long Short-Term Memory, LSTM）和反向传播神经网络（Back Propagation Neural Network, BPNN）方法进行对比，结果表明：结合触发与失稳机制的PIDA方法在预测精度方面表现更优。此外，该方法能够实现动态、可靠的滑坡变形预测。

1. 引言

滑坡因其对人民生命财产安全造成的灾难性影响，已成为全球广泛关注的重要地质灾害问题之一（Cruden, 1991；Froude and Petley, 2018；Liu and Wang, 2021；Guo et al., 2021）。滑坡位移监测与预测已被证明是滑坡早期预警体系中不可或缺的组成部分（Intrieri et al., 2012；Li et al., 2021c）。自 Saito（1965）和 Fukuzono（1985）开创性工作以来，众多基于物理和现象学的模型相继被提出并不断完善（Miao et al., 2017）。然而，受制于滑坡岩土结构的复杂性和失稳机制多样性的影响，传统物理力学方法的发展受到了严重限制（Huang et al., 2017）。虽然基于非线性动力学和数值模拟的物理模型在稳定性评价中考虑了地质因素、降雨和水动力学模型，但其过于简化的建模方式难以准确刻画滑坡的动态演化过程（Chen et al., 2018）。此外，大多数具有物理意义的预测模型通常是针对特定个例工况构建（Jiang et al., 2010；Yan et al., 2020），缺乏连续更新，导致难以动态研究滑坡演化过程，从而无法实现理想的分析与预测效果。

基于观测数据的预测方法近年来逐渐成为传统模型驱动方法的重要补充（Federico et al., 2011；He and Semnani, 2023）。随着数据采集与处理技术的发展，变形监测逐渐被视为滑坡预测的有效手段（Gili et al., 2000；Casagli et al., 2010；Gong et al., 2021；Wang et al., 2022）。依托非线性建模能力，机器学习（ML）和深度学习（DL）方法已广泛应用于滑坡预测领域（Mayoraz and Vulliet, 2002；Zhang et al., 2021）。同时，部分集成多算法的改进型集成预测模型也被提出，针对不同类型滑坡和时间序列特性提升预测性能（Zhou et al., 2016；Han et al., 2021；Lin et al., 2022）。但若数据挖掘方法未能融合滑坡的物理机理，则其泛化能力仍存在一定局限性。总体而言，常规数值模拟与物理预测模型多基于静态模型，在特定工况下实现预测，难以动态表征滑坡演化过程及其触发因素（Zhang et al., 2019）。数据驱动模型凭借强大的数据处理能力，得益于监测数据量的激增而快速发展，但大多未涉及物理机制分析，其数学模型存在“黑箱”问题，无法获取内部变量及参数，缺乏可解释性（Bathaee, 2018）。此外，目前仍缺乏将监测数据反馈用于修正预测模型、提升预测精度的研究。

近年来，结合滑坡触发机制与失稳机制，构建物理约束方法以实现滑坡变形动态与可靠预测的趋势愈加显著（Yan et al., 2022）。鉴于当前研究中存在的问题，综合考虑滑坡内部与外部动态条件，数据同化（Data Assimilation, DA）方法凭借其模型参数持续优化能力，成为理想的候选方案。DA方法能够将机制模型与多源观测数据融合，对模型参数与预测轨迹进行更新修正（Rodell et al., 2004），目前已在气象学、海洋学及陆表过程等领域得到深入研究与广泛应用（Robinson and Lermusiaux, 2000；Evensen, 2003；Brezzi et al., 2015）。然而，该方法在滑坡灾害领域的应用仍较为有限。例如，Jiang 等（2016）基于序贯DA框架，耦合合成孔径雷达（SAR）数据与水文因子，开展了滑坡响应机制分析；Xue 等（2018）将DA方法应用于基于降雨入渗与区域格网坡面稳定性（TRIGRS）模型、GPS与InSAR数据的滑坡稳定性分析。然而，InSAR和SAR时序数据在获取与处理方面存在一定局限（Jiang et al., 2016），且TRIGRS模型主要针对浅层、降雨型滑坡（Baum et al., 2008）。此外，现有滑坡预测中DA方法多关注表面位移，忽略了能够更好反映滑坡变形机制的深部位移数据（Shentu et al., 2011）。

针对具有独特地质结构演化与力学强度劣化规律的水动力压力驱动型滑坡，本文提出了一种基于物理约束的数据同化（Physics-Informed Data Assimilation, PIDA）方法用于位移预测。该方法旨在提升物理预测模型的预测能力，解决数学模型忽略滑坡动态变形物理意义的问题。

2. 物理约束数据同化（PIDA）方法

2.1 粒子滤波（Particle Filter, PF）算法

粒子滤波（PF）算法是一种基于贝叶斯估计和蒙特卡洛方法（Monte Carlo method）框架发展而来的常用数据同化方法，因其适用于非线性、非高斯系统而受到广泛关注（Chopin, 2002；Han and Li, 2008；Fearnhead and Künsch, 2018）。贝叶斯滤波理论利用系统的状态转移方程与前一时刻的后验概率密度，预测系统状态的先验概率密度（Nakano et al., 2007；Gao and Zhang, 2012；Dahlin and Lindsten, 2014）。动态系统的状态转移规则可由状态方程表示为：

$$x_t=f(x_{t-1},{\epsilon }_t)$$

以及观测方程：

$$y_t=g(x_t,{\eta }_t)$$

其中，$x\_t$ 和 $x\_t-1$ 分别表示当前时刻 $t$ 和上一时刻 $t-1$ 的状态变量，假设 $x\_t$ 满足一阶马尔可夫过程，即当前状态仅与前一时刻状态相关。$y\_t$ 表示当前时刻 $t$ 的观测变量，$f(\cdot )$ 为状态转移方程，$g(\cdot )$ 为观测方程，$\epsilon \_t$ 和 $\eta \_t$ 分别表示状态预测过程噪声与观测过程噪声。贝叶斯滤波通过观测变量 $y\_1t$ 递推状态变量 $x\_t$ 的后验概率密度，获取最优状态估计（Doucet et al., 2001）。该过程包括预测与更新两部分，前者根据上一时刻的状态与系统模型预测下一时刻状态的先验概率密度，后者在获得当前观测值后修正先验概率密度，得到状态的后验概率密度。

预测阶段，通过上一时刻 $t-1$ 的后验概率 $p(x\_t-1|y\_1t-1)$，预测当前时刻 $t$ 的先验概率 $p(x\_t|y\_1t-1)$，公式如下：

$$p(x_t|y_{1:t-1})=\int p(x_t|x_{t-1})p(x_{t-1}|y_{1:t-1})d{x}_{t-1}$$

其中，$p(x\_t|x\_t-1)$ 为系统模型的状态转移概率。
更新阶段，利用观测变量 $y\_t$ 修正先验概率 $p(x\_t|y\_1t-1)$，得到当前时刻的后验概率 $p(x\_t|y\_1t)$，表达式如下：

$$p(x_t|y_{1:t})=\frac{p(y_t|x_t)p(x_t|y_{1:t-1})}{p(y_t|y_{1:t-1})}$$

其中，$p(y\_t|x\_t)$ 表示当前状态变量 $x\_t$ 与观测值 $y\_t$ 的似然概率，$p(y\_t|y\_1t-1)$ 称为证据（evidence），即对观测值的一步预测，通常为归一化常数。最终得到的后验概率将作为下一预测时刻的输入，完成递推。

根据状态变量的后验概率密度 $p(x\_t|y\_1t)$，可计算状态变量的最小方差估计：

$$E(x_t)=\int x_{t}p(x_t|y_{1:t})d{x}_t$$

实际应用中，预测过程中涉及的积分运算较难实现，且状态变量后验概率密度一般难以解析表达。蒙特卡洛仿真方法通过抽取服从后验分布的一系列样本，将积分运算问题转化为样本期望求解。设定从先验概率分布中独立抽取 $N$ 个随机粒子 $x\_t^i\_{i=1}^N$，则状态变量的后验概率密度近似表示为：

$$p(x_t|y_{1:t})\approx \hat{p}(x_t|y_{1:t})=\frac{1}{N}\sum _{i=1}^N\delta (x_t-x_t^i)$$

其中，$\delta (\cdot )$ 为狄拉克函数，$\hat{p}(x\_t|y\_1t)$ 表示 $p(x\_t|y\_1t)$ 的近似值。

状态变量 $x\_t$ 的估计值则可表示为：

$$E(x_t)=\frac{1}{N}\sum _{i=1}^N{x}_t^i$$

在PF方法中，根据系统的状态经验分布，在状态空间内随机生成一组“粒子”来表示待估计的状态变量。该方法通过观测值修正粒子权重，进而有效估计状态变量的概率密度分布及其最优值（Green, 1995；Arulampalam et al., 2002）。理论上，依据大数定律，只要粒子数量足够多，便可较准确地模拟实际后验概率密度分布。PF算法的实现过程如图1所示，并可归纳为以下步骤（Djuric et al., 2003）：

步骤1：初始化粒子样本，设置每个粒子的初始权重为 $1/N$。

步骤2：根据状态方程（式(1)）预测粒子状态 $x\_t^i,,i=1,2,\cdots ,N$。

步骤3：依据观测值更新粒子权重，并对权重进行归一化。由于后验概率分布未知，无法直接通过粒子模拟，因此PF方法采用重要性采样（importance sampling）方法（Chen et al., 2005）。选择已知的重要性概率密度函数 $q(x\_t|y\_1t)$，从中抽取粒子，再通过对粒子加权模拟真实后验概率密度。后验概率与状态变量最优估计表达如下：

$$p(x_t|y_{1:t})=\frac{p(y_{1:t}|x_t)p(x_t)}{q(x_t|y_{1:t})p(y_{1:t})}$$

粒子的估计值为：

$$E(x_t)=\frac{{\sum }_{i=1}^N{x}_t^i{\omega }_t(x_t^i)}{{\sum }_{i=1}^N{\omega }_t(x_t^i)}=\sum _{i=1}^N{x}_t^i{\tilde{\omega }}_t(x_t^i)$$

其中，$\omega \_t(x\_t)$ 为重要性权重，$\omega \_t=\frac{p(y\_1t|x\_t)p(x\_t)}{q(x\_t|y\_1t)}$，$\tilde{\omega }\ast t(x\_t^i)=\frac{\omega \_t(x\_t^i)}{\sum \ast {i=1}^N\omega \_t(x\_t^i)}$ 为归一化后的权重值。

重要性概率密度函数的选取需尽可能减小权重方差。为简化计算，通常将先验概率密度函数 $p(x\_t|x\_t-1)$ 作为重要性概率密度函数，此时粒子权重递推公式可表示为：

$${\omega }_t(x_t)={\omega }_{t-1}(x_{t-1})p(y_t|x_t)$$

步骤4：考虑到“粒子退化”（particle degradation）问题（Doucet et al., 2000），采用重采样（resampling）方法保留权重较大的粒子，舍弃权重较小的粒子（Gordon et al., 1993；Liu and Chen, 1995），并将所有粒子的权重重置为 $1/N$，随后进入迭代。

步骤5：依据概率密度分布，计算当前时刻 $t$ 的状态变量最优估计值，更新 $t=t+1$，返回步骤2，直至全部时刻预测完成。

2.2. 构建物理信息驱动的滑坡位移预测模型

结合滑坡的物理机制模型与数学预测模型，本文构建了物理信息驱动的滑坡位移预测模型，整体框架如图 2 所示。首先，依据深部位移监测数据计算滑带劣化指标（SZD），并构建包含前一时刻 SZD 和外部诱发因子的粒子滤波（PF）算法状态方程，将 SZD 指标作为需要在 PF 迭代中估计的模型参数。其次，建立 SZD 指标与剪切强度参数之间的相关关系，从而动态更新剪切强度参数。

接着，通过滑坡变形的力学分析与潜水面位置确定，构建物理机制模型。最后，基于滑坡内力分布与变形监测数据建立观测方程，形成完整的预测模型。该模型考虑滑坡的动态演化过程，利用 PF 算法，通过实时地表位移观测值与预测值的差异不断更新状态方程的参数，从而修正预测轨迹。

2.2.1 滑带劣化（SZD）模型

监测数据显示，滑坡破坏通常伴随着潜在滑动带的结构破坏（Li 等, 2021a）。因此，提出滑带劣化指标（SZD）以反映滑动带的结构完整性及其对外界诱发因素的动态响应能力。

首先，滑动带位置可通过钻探及深部位移曲线确定，因为滑动带附近常伴有突变（Juang, 2021；Li 等, 2021b）。然后，依据深部位移–深度曲线（见图 3），可计算 SZD 指标（Li 等, 2008）。其中，钻孔倾斜仪是常见的深部位移测量设备（Segalini 等, 2014；Ha 等, 2018）。设第 $j$ 个监测点的累积位移为 $W\_j$（$j=0,1,2,\dots ,n$，$n$ 为监测点数），定义为该点在时刻 $t$ 相对于初始时刻基准点（$j=0$）的水平位移差。

由于受传感器数量限制，通常选取 4 个代表性监测点进行 SZD 定量分析。图 3a 所示，传感器 3 安装于钻孔顶部，位于上滑面之上；传感器 0 安装于钻孔底部，作为固定基准点，位于下滑面之下；传感器 1 和 2 均安装于滑动带内，通常具有较大位移响应，表征剪切破坏。

为简化 SZD 计算，深部位移–深度曲线被折线近似处理（图 3b），以更好地刻画滑动带剪切变形的几何信息。如图 3c 所示，滑坡动态演化过程中，滑动带内的相对变形 $\Delta W\_12=W\_2-W\_1$ 随时间增加，而滑带上方滑体的相对变形 $\Delta W\_23=W\_3-W\_2$ 相对减小。因此，定义 SZD 指标 $S\_t$ 为：

$$S_t=\frac{\Delta W_{12}}{\Delta W_{12}+\Delta W_{23}}\text{\hspace{1em}(11)}$$

当 $\Delta W\_12=0$ 时，$S\_t=0$，表示滑动带未发生变形，结构完整性最高；当 $\Delta W\_23=0$ 时，$S\_t=1$，说明滑带发生严重剪切破坏，滑坡失稳，结构完整性最低。此时，滑坡对外部诱因（如降雨、水位波动等）响应增强。

$S\_t$ 的取值范围为 KaTeX parse error: Undefined control sequence: \[ at position 1: \̲[̲0,1]，具有明确物理意义，并随滑动带劣化呈递增趋势。鉴于监测数据稀缺，依据滑动带结构受水动力作用劣化的响应机制，构建了 SZD 序列预测模型作为 PF 状态函数，以提升其任意时刻的预测能力。其形式为：

$$S_t=f_t(S_{t-1},u_t)+{\epsilon }_t\text{\hspace{1em}(12)}$$

其中，$S\_t$ 与 $S\_t-1$ 分别表示当前与前一时刻的 SZD 指标，$u\_t$ 表示相关外部诱发因子时间序列，$\epsilon \_t$ 表示预测过程中的扰动项。

2.2.2 物理机制模型

I. 剪切强度参数动态更新模型

滑坡演化过程中，材料组成与结构特征不断演变（Xue 等, 2020），滑动带岩土的强度参数也随外部条件变化而动态变化（Springman 等, 2003；Zuo 等, 2020）。然而，现场或实验手段难以获取其真实数值。实际上，从初始状态至完全破坏状态，可视为剪切强度逐渐下降的宏观结果。因此，将剪切强度参数建模为 SZD 指标的函数。

设 $S\_t$ 为 SZD 指标，$c\_t$ 和 $\varphi \_t$ 为时刻 $t$ 的黏聚力和内摩擦角，模型如下：

$$c_t=c_f-\left(\frac{c_f-c_r}{S_r-S_f}\right)\times (S_t-S_f)\text{\hspace{1em}(13)}$$

$$\tan {\varphi }_t=\tan {\varphi }_f-\left(\frac{\tan {\varphi }_f-\tan {\varphi }_r}{S_r-S_f}\right)\times (S_t-S_f)\text{\hspace{1em}(14)}$$

其中，$S\_f$ 表示完整结构下的 SZD 指标（一般取 0），$S\_r$ 表示完全破坏时的 SZD 指标（一般取 1）。代入简化得：

$$c_t=c_f-(c_f-c_r)\times S_t\text{\hspace{1em}(15)}$$

$${\varphi }_t=\arctan [\tan {\varphi }_r+(\tan {\varphi }_f-\tan {\varphi }_r)(1-S_t)]\text{\hspace{1em}(16)}$$

其中，$c\_f$、$\varphi \_f$、$c\_r$、$\varphi \_r$ 分别为滑带处于饱和状态或自然状态下的峰值与残余强度参数。

II. 潜水面位置计算模型

水动力作用下的潜水面位置对滑坡稳定性影响显著，但计算复杂（Yan 等, 2010）。根据 Boussinesq 非稳定渗流方程，并结合 Laplace 变换与反变换，可推导在水库水位变化和降雨作用下的简化模型。

假设如下：

• 含水层均质各向同性，水平延展无限；
• 两侧水库渗流为一维流；
• 水位变化导致平行坡面渗流，降雨引起垂直渗流；
• 水库水位以恒速 $v\_0$ 升降；
• 岸坡为垂直边坡。

非稳定渗流控制方程：

$$\mu \frac{\partial h}{\partial t}=\frac{\partial }{\partial x}\left(T\frac{\partial h}{\partial x}\right)+\varpi \text{\hspace{1em}(17)}$$

其中，$\mu$ 为有效含水率，$\varpi$ 为降雨强度，$h$ 为水头，$T=Kh$ 为渗透系数（$K$ 为渗透率）。该二阶非线性方程可简化为：

$$\frac{\partial h}{\partial t}=\sigma \frac{{\partial }^{2}h}{\partial x^2}+\frac{\varpi }{\mu }\text{\hspace{1em}(18)}$$

其中，$\sigma =\frac{KH\_m}{\mu }$ 为渗透导率，$H\_m$ 为平均含水层厚度。

水位上升时：

$$u(x,t{)}_{\text{up}}=h(x,t{)}_{\text{up}}-h(0,0)\text{\hspace{1em}(19)}$$

其中，$h(0,0)$ 表示初始时刻不透水层与水库水位之间的距离，$h(x,t)\_up$ 是在时间 $t$、位置 $x$ 处的潜水面高度。滑坡潜水面的非稳定渗流数学模型表达如下：

$$\{\begin{array}{rl}& \frac{\partial h}{\partial t}=\sigma \left(\frac{{\partial }^{2}h}{\partial x^2}\right)+\frac{\varpi }{\mu },0<x<\infty ,t>0\\ & u(0,t)=h(0,t)-h(0,0)=h(0,0)+v_{0}t-h(0,0)=v_{0}t,t>0\\ & u(x,0)=h(x,0)-h(0,0)=0,0<x<\infty \\ & u(\infty ,t)=h(\infty ,t)-h(0,0)=h(0,0)+\frac{\varpi t}{\mu }-h(0,0)=\frac{\varpi t}{\mu },t>0\end{array}\text{\hspace{1em}(20)}$$
其中，$v\_0$ 表示水库水位上下波动的恒定速度。 当水库水位下降时（见图 4b），位置 $x$ 处在时刻 $t$ 的潜水面变化可表示为：

当水库水位下降（见图4b）时，位置 $x$ 处在时刻 $t$ 的潜水面变化可表示为：

$$u(x,t{)}_{down}=h(0,0)-h(x,t{)}_{down}\text{\hspace{1em}(21)}$$
滑坡潜水面的非稳定渗流数学模型计算如下：

$$\{\begin{array}{rl}& \frac{\partial h}{\partial t}=\sigma \left(\frac{{\partial }^{2}h}{\partial x^2}\right)+\frac{\varpi }{\mu },0<x<\infty ,t>0\\ & u(0,t)=h(0,t)-h(0,0)=v_{0}t,t>0\\ & u(x,0)=h(x,0)-h(0,0)=0,0<x<\infty \\ & u(\infty ,t)=h(\infty ,t)-h(0,0)=\frac{\varpi t}{\mu },t>0\end{array}\text{\hspace{1em}(22)}$$

求解该 PDE 模型后，应用回归分析拟合参数，最终得到水位上升与下降时的解析解：

$$h(x,t{)}_{\text{up}}=\left(\frac{v_0\mu }{\pi }\right)(a_1{\lambda }^4+a_2{\lambda }^3+a_3{\lambda }^2+a_4\lambda +1)t+\frac{\pi }{\mu }t\text{\hspace{1em}(23a)}$$

$$h(x,t{)}_{\text{down}}=h(0,0)-\left(\frac{v_0\mu }{\pi }\right)(a_1^{\prime }{\lambda }^4+a_2^{\prime }{\lambda }^3+a_3^{\prime }{\lambda }^2+a_4^{\prime }\lambda +1)t-\frac{\pi }{\mu }t\text{\hspace{1em}(23b)}$$

其中，$\lambda =\frac{x}{2\sqrt{KH\_mt/\mu }}$，$a\_i$ 和 $a\_i^{\prime }$ 为拟合系数。

2.3. 迭代计算过程

边坡稳定性一直是岩土工程领域的重要研究课题（Liao and Ji, 2021；Zhang et al., 2022）。然而，传统的整体安全系数和失稳分析方法难以反映滑坡的动态演化过程（Miao et al., 2017）。以 Mohr–Coulomb 强度准则为基础的极限平衡法（Limit Equilibrium Method, LEM）被广泛应用于滑坡稳定性分析，其通常假设岩体为刚体，且通过条分法（slice method）实现滑体的极限平衡分析（Xu et al., 2021；Su et al., 2022）。其中，残余推力法（Residual Thrust Method）因其实现简便，常用于分析滑坡物理机制（Song and Xu, 2012；Zou et al., 2020a）。

在复杂荷载条件下，可通过滑坡结构特征分析与残余推力计算得到任意形状滑动面的滑坡推力。将滑体沿滑面方向等分为 $M$ 个条带，对第 $i$ 个条带进行受力分析（如图 4c 所示），并考虑地下水的影响，其计算公式如下：

$$F_i=\frac{\left[(W_{1i}+W_{2i})\sin {\alpha }_i+D_i\cos ({\alpha }_i-{\beta }_i)\right]c_i{l}_i+\left[(W_{1i}+W_{2i})\cos {\alpha }_i-D_i\sin ({\alpha }_i-{\beta }_i)\right]\tan {\varphi }_i}{FS}+F_{i-1}{\psi }_i\text{\hspace{1em}(24)}$$

$${\psi }_i=\frac{\cos ({\alpha }_{i-1}-{\alpha }_i)}{\sin ({\alpha }_{i-1}-{\alpha }_i)\tan {\varphi }_i/FS}\text{\hspace{1em}(25)}$$

其中，$F\_i$ 表示第 $i$ 个条带的残余滑动力，$FS$ 为边坡稳定系数；$W\_1i$ 和 $W\_2i$ 分别表示条带 $i$ 上部饱和线以上部分的重力与饱和线以下部分的浮重；$\alpha \_i$ 为第 $i$ 条带的滑动面倾角；$c\_i$ 和 $\varphi \_i$ 分别表示条带的黏聚力和内摩擦角；$D\_i=\gamma \_wA\_i\sin \beta \_i$ 表示渗透水压力，$\gamma \_w$ 为水的单位重，$A\_i$ 为条带 $i$ 位于饱和线以下的面积，$\beta \_i$ 为水力坡降方向与水平线的夹角；$\psi \_i$ 为条带间残余推力的传递系数。式（24）右侧依次表示滑动力、抗滑力与残余推力。

采用迭代计算方法求解该方程：初始时将第一个条带的残余推力设为 0，逐条带计算残余滑动力。当最后一条带的残余推力为 0 时，迭代终止，此时各条带的残余推力可反映滑坡的受力变形过程。

2.3.1 推力–位移模型

基于大量实验与现场调查结果，Lu（2015）提出了一种用于描述不同法向应力条件下滑面剪应力–剪应变全过程特性的岩土接触面本构模型（Sun et al., 2022）。对于特定岩土体，其力学行为仅随法向应力变化。在此基础上，用四个参数即可描述其剪切行为，其本构关系可表示为：

$$\tau =G\gamma {\left[1+\frac{\gamma m}{s}\right]}^{\rho }\text{\hspace{1em}(26)}$$

其中，$\tau$ 表示剪应力，$\gamma$ 为剪应变，$G$ 为初始剪切模量，$m$、$s$ 和 $\rho$ 为在不同法向应力下的常数系数。

该模型并非严格理论推导的物理模型，而是基于数学拟合过程的数学模型，可拟合岩土体的剪应力–应变关系；只要满足函数曲线特性，也可用于拟合其他物理或化学过程。

为了构建同时包含力学参数与监测数据的观测方程，本文将剪应力替换为残余推力 $F$，将剪应力替换为滑坡的月累积位移 $Y$，并改写模型如下：

$$F=GY{\left[1+\frac{Ym}{s}\right]}^{\rho }\text{\hspace{1em}(27)}$$

其中，$F$ 表示残余推力，$Y$ 表示月位移，$G$ 为初始剪切模量，$m$、$s$ 和 $\rho$ 为常数系数。

因此，滑坡推力与表面监测数据之间的观测方程可表示为：

$$Y_t=g_t(F_t)+{\eta }_t\text{\hspace{1em}(28)}$$

其中，$Y\_t$ 表示时刻 $t$ 的观测累积表面位移，$F\_t$ 为时刻 $t$ 滑带的推力，$g\_t(\cdot )$ 为反映力分布与表面位移关系的数据驱动函数。

2.3.2 PIDA 方法的迭代步骤

PIDA 方法结合了物理约束滑坡位移预测模型与粒子滤波算法，旨在提升地质参数估计的精度。其实现流程如下（见图1）：

1. 初始化粒子：将 SZD 指标（滑带劣化指标）视为粒子。首先，取 SZD 指标初值 $S\_t$ 的均值，初始化一组服从高斯分布的粒子，使其在 SZD 指标范围内分布广泛。

2. 预测位移与加权更新：根据每个粒子的参数，使用物理约束滑坡预测模型计算对应的月表面位移预测值。将预测值与实测值进行对比，根据重要性概率密度函数 $q(x\_t|y\_1t)$ 计算对应粒子的权重 $\omega \_t(x\_t)$。权重越大，表示该粒子所代表的 SZD 指标越接近实际滑坡状态。

3. 加权估计最优值：根据粒子权重对 SZD 指标加权求和，获得其最优估计值。在更新粒子权重过程中，保留与实测月位移最接近的粒子预测结果，作为最终预测值。

4. 重采样操作：采用重采样方法保留高权重粒子，剔除低权重粒子，避免无效计算。根据当前时刻的最优 SZD 估计值与重采样后的粒子，递推生成下一时刻的 SZD 指标预测值 $S\_t+1$ 与对应粒子。

通过上述过程，PIDA 模型基于 PF 算法实现滑坡位移的连续预测。在每次递推中利用粒子对模型参数不断优化，从而提升滑坡预测模型的精度。

3. 应用与验证

3.1. 案例研究

白水河滑坡位于三峡库区，是典型的水动力压力驱动型滑坡，具有阶梯式变形特征（Zou 等，2020b）（图5a）。自2003年7月起，陆续布设了包括11个GPS监测点和3个钻孔测斜仪在内的多种监测设备，用于地表和深部位移监测（Liu 等，2020；Long 等，2022）。白水河滑坡的监测点布置及剖面示意分别见图5b和图5c。滑坡从2003年6月至2016年12月的累计位移、降雨和水库水位监测曲线见图6a。2007年7月，GPS监测显示白水河滑坡发生了剧烈变形，立即发布预警并关闭长江航道，造成每日损失逾亿元及较大社会负面影响（Li 等，2009）。随后监测显示，滑坡变形减缓并逐渐稳定持续十余年。这对经验模型或数据驱动模型的滑坡预测提出了新挑战。因此，以白水河滑坡为例，结合监测系统与物理机制的PIDA方法为急需解决的问题提供了潜在方案。

3.2. 位移预测

SZD指标是连接粒子滤波（PF）算法与物理机制模型的桥梁。选取2013年5月至2015年5月的深部位移数据计算SZD指标（图6b），计算结果见图7（蓝线）。同时，利用灰色关联度法（Yang 等，2019）筛选出的两大外部因素——月降雨量和水库水位下降量（2013年5月至2016年6月）用于构建序列SZD指标预测模型，采用多源回归分析，初始最优参数由粒子群优化（PSO）算法确定。状态函数如下：

$$S_t={\alpha }_0{S}_{t-1}+{\alpha }_1{R}_t+{\alpha }_2{R}_t{K}_t+{\alpha }_3{R}_t+{\alpha }_4+{\epsilon }_t$$

其中，Rt和Kt分别为时刻t的月降雨量和水库水位下降量；α0、α1、α2、α3和α4为对应回归系数。计算初值为α0=0.2224，α1=2.490×10⁴，α2=0.046×10⁴，α3=0.0006，α4=0.4874。εt为系统过程噪声协方差。当预测位移与观测位移误差大于0.01时，PF算法过程将更新这些系数。

选取白水河滑坡1-1′剖面作为横断面（图5c），按照条分法将坡体划分为84个条分，平均长度9.66米，最后一个条分9.06米（图8a）。条分编号自左至右为1至84号，ZG93监测点对应第39号条分。根据岩土勘察报告和实验（Sassa 等，2004），白水河滑坡峰值和残余剪切强度参数的平均值见表1。滑动带土体的自然和饱和单重分别为22 KN/m³和22.5 KN/m³。渗透系数K=0.1，比排系数μ=0.25。根据式（23）通过回归分析拟合方程系数，可得2013年5月至2016年6月的月度潜水线位置。其中，2013年10月为水库水位上升期典型阶段，2016年2月为下降期。潜水线位置分别见图8b和图8c。潜水线计算结果与相同条件下数值模拟吻合良好。滑动面可根据其与潜水线的相对位置判断为自然状态或饱和状态。

结合滑坡SZD指标模型及剪切强度参数类型，利用式（15）和（16）可获得每个条分的动态剪切强度参数。将动态剪切强度参数代入式（24）和（25）迭代计算，得各条分的残余推力。基于2013年5月至2015年5月ZG93点月度地表位移监测数据，采用数据驱动方法建立推力-位移数学模型，参数通过PSO方法优化。观测方程为：

$$G{Y}_t{\left[1+\frac{Y_t^m}{s}\right]}^{\rho }=g_t(F_t)+{\eta }_t\text{\hspace{1em}(30)}$$

计算得m=1，s=7.15×10⁶，ρ=0.53。进一步计算2015年6月至2016年6月第39条分的残余推力用于滑坡累计位移预测。

PIDA方法基于PF算法，SZD指标动态范围为0~1，粒子数量设为100，初始化采用均值为0.4的高斯分布。根据观测月度位移与预测值差异，动态更新SZD指标预测模型参数实现优化。最终，基于PIDA方法实现了参数动态更新与调整的准确滑坡位移预测。

3.3. 结果验证

人工神经网络（ANNs）因其非线性优势和卓越预测能力，常用于滑坡位移预测。ANN通常由输入层、隐藏层和输出层组成。反向传播神经网络（BPNN）是一种ANN，信息沿网络正向流动，误差通过反向传播调整权重（Li 等，2014）。长短时记忆网络（LSTM）设计考虑了时间序列，含门控机制控制信息流入与流出，实现选择性记忆与遗忘（Yang 等，2019）。

为验证基于PF算法的PIDA方法的准确性与可靠性，采用LSTM和BPNN方法对白水河滑坡位移进行对比预测。监测数据（2013年5月至2015年5月）按70%训练、30%测试划分，2015年6月至2016年6月的位移值由训练模型预测。图9a显示，PIDA方法预测位移绝对误差不超过10毫米，大部分集中于5毫米以下。

先前研究表明，白水河滑坡的变形受季节性降雨和水库水位变化驱动（Tang 等，2019）。图6a中，滑坡累计位移曲线呈阶梯式递增趋势，每年5月至7月出现明显阶梯，对应水库水位下降期及低水位期（或雨季）。通常，6月或7月暴雨后或水库水位快速下降时出现变形“跳变”，属变异状态事件；而其他如水位上升期和高水位期，累计位移曲线趋于平缓，滑坡处于蠕变状态。因此，2015年6月位移显著增加，且滑坡突变变形状态较难预测且极易失稳。

值得注意的是，PIDA方法在2015年6月“跳变点”处的绝对误差极小，甚至接近零，表现出极佳预测性能。PIDA方法、LSTM和BPNN预测结果对比表明，PIDA在滑坡突变状态预测中表现最优，因其能利用同化算法快速调整参数应对实际突变。图9a和9b还显示，PIDA预测值与观测数据最为一致，LSTM次之，BPNN在随后相对稳定阶段表现最差。

传统上，均方根误差（RMSE）和平均绝对误差（MAE）用于评价模型整体预测精度，计算公式如下：

• 均方根误差 (RMSE)：

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}\sum _{t=0}^n{\left({\hat{Y}}_t-Y_t\right)}^2}$$

• 平均绝对误差 (MAE)：

$$\mathrm{MAE}=\frac{1}{n}\sum _{t=0}^n\left|{\hat{Y}}_t-Y_t\right|$$

其中，Yt为观测值，Ŷt为预测值。PIDA、LSTM和BPNN方法的RMSE分别为2.810、27.271和47.215，MAE分别为1.419、16.495和19.007（图9c）。结果证明，考虑滑坡动态物理机制及监测数据的PIDA方法在滑坡预测中明显优于纯数据驱动方法，同时验证了参数动态调整的有效性。

4. 讨论

所提出的PIDA方法不仅考虑了滑坡的物理机制，还实现了滑坡预测参数的反馈、修正与动态更新，能够准确刻画滑坡的动态变形响应。在PIDA模型构建过程中，SZD指标由深部位移数据计算并通过序列SZD模型推导得出，是连接状态方程与观测方程的关键桥梁。一方面，状态方程中的SZD指标反映了滑坡结构完整性及其对外部诱发条件的动态响应；另一方面，也可以选择其他代表滑坡结构退化的指标，基于粒子滤波算法构建状态方程。通过建立剪切强度参数与SZD指标之间的关系，实现剪切强度参数的动态更新，进而构建物理机制模型。

本文采用剩余推力法分析滑坡稳定性并计算各条分单元的动态推力，因其具有迭代过程的优势。毫无疑问，更可靠或改进的稳定性分析方法也可以应用于PIDA方法以提升性能。由本构模型推导的推力-位移模型被视为滑坡位移预测的观测方程。然而，未来可进一步研究力分布与地表位移的关系，建立更可靠的观测方程以提升计算精度。

此外，值得关注的是PIDA方法在滑坡预测中参数更新与修正方面的显著优势。同时，粒子滤波算法中存在粒子退化问题——部分粒子具有很高的适应度权重，而大部分粒子权重较低甚至接近零，导致计算资源浪费且难以准确反映真实后验概率分布。重采样是减少无效样本、增加有效样本的有效手段，但多次重采样可能削弱粒子多样性，影响滑坡物理预测模型的适应性。因此，改进粒子滤波算法以解决粒子退化问题仍需更多努力。

本研究所用PIDA预测模型采用了月度降雨量、水库水位、深部及地表位移监测数据进行更新与预测，然而较粗的采样间隔无法充分捕捉连续的变形过程，进而降低预测精度。为提升模型时间精度，增加监测数据采集频率或结合数值方法采用适当插值技术补充数据将更为有效。

5. 结论

本文提出了一种考虑物理机制演化过程与数据驱动预测的动态PIDA方法。该模型能够根据位移偏差实时更新与修正关键参数，从而提升滑坡预测精度。在滑坡物理机制模型中，基于深部位移数据和相关外部诱发因素计算SZD指标；利用Boussinesq方程识别滑动带潜水面位置；通过SZD指标动态获取剪切强度参数；采用条分法与剩余推力法计算滑坡内部推力。基于粒子滤波算法的PIDA方法实现了剪切强度的动态更新，利用自适应动态参数增强了模型的适应能力，使预测结果更符合滑坡滑动前的实际变化。

对比结果表明，基于粒子滤波算法的PIDA方法在滑坡预测中优于LSTM和BPNN方法，尤其在变形突变点的预测上表现突出。PIDA模型不仅具备物理机制意义，还结合了监测数据与外部诱发因子，能够估计预测模型的状态变量与SZD指标，减少剪切强度的不确定性与偏差，有效降低模型结构误差，提升预测准确性。总体而言，该模型为滑坡分析与预测提供了全新的思路与方法。