用于位移预测的物理信息数据同化模型

刘勇 ^一个， 龙晶晶 ^b,c， 李昌东 ^b,c,*， 詹伟文 ^一个

^一个 中国 地质大学 机械工程与电子信息学院， 湖北 武汉 430074

^b 中国 地质大学 （武汉） 工程学院 ， 中国 武汉 430074

^c 中国 地质大学八东国家地质灾害观测研究 站 ， 武汉 430074， 中国

A R T I C L E I N F O

关键字：

水动力压力驱动滑坡预测模型

动力响应机制 滑动区劣化 颗粒过滤器

物理信息数据同化方法

A B S T R A C T

提出了一种物理，数学预测模型中滑坡意义的问题 PF 算法方程由于滑坡演化的复杂性和非线性特性，滑坡的变形规律总是难以预测。考虑到滑坡物理机理的动力学演化过程，信息数据同化 （PIDA） 方法首先 进行滑坡位移预测。 它可以 增强 预测能力 物理预测模型的并解决 忽略 物理。此外，粒子滤波器 （PF） 在非线性和非高斯系统中的优越性决定了其在复杂问题中的应用。针对流体动力压力驱动的滑坡， 了滑带劣化 （SZD） 模型 构建 作为 状态 的。观测方程是由滑坡变形与监测数据集之间的关系形成的。基于物理学的滑坡位移预测模型作为桥梁，由抗剪强度参数动态更新模型和物理机理分析模型组成，用于关键坡度反馈、修改和更新。将结果与传统的长短期记忆 （LSTM） 和反向传播神经网络 （BPNN） 方法进行比较，PIDA 方法结合了 触发机制和不稳定机制在准确性方面表现更好。此外，它还能够动态可靠地预测滑坡变形。

介绍

，造成了灾难性的影响（人 的以来物理力学方法复杂性的严重制约模拟的物理模型在估算边坡稳定性时山体滑坡是引起全世界关注的最重要问题之一因为它对人们的生活和财产安全Cruden，1991 年;Froude 和 Petley，2018 年;Liu 和 Wang，2021 年; Guo et al.， 2021）。位移监测和预测已被证明 是 不可或缺 的部分滑坡早期预警（Intrieri 等， 2012 年; Li et al.， 2021c）。 自 开创性工作 Saito （1965） 和 Fukuzono （1985），包括 基于物理 和现象学的模型已经被提出并进一步改进（Miao et al.， 2017）。 然而， 的发展 传统受到滑坡岩 土结构，而 多样性的失败机制（Huang et al.， 2017）。基于非线性动力学和物理考虑了地质因素、降水和 hydraulic 模型，但过于简化的模型几乎无法 表达 动力学演变 （山体滑坡Chen et al.， 2018）。 此外，大多数物理上有意义的预测模型

; Yan 。are mainly constructed under case-specific conditions (Jiang et al., 2010; Yan et al., 2020). Lack of continuous updates makes it difficult to study the landslide dynamic evolution process. Consequently, it is nearly impossible to achieve ideal analysis and prediction results.

数据采集和处理的发展，变形那么通用性存在固有的缺点数据驱动方法依赖于对监测数据的观察和解释，这已成为传统模型驱动方法的替代方案（Federico 等人，2011 年;He 和 Semnani，2023 年）。 随着监测 被认为是 一种有效的滑坡预测方法技术（Gili et al.， 2000; Casagli 等人，2010 年; Gong et al.， 2021;Wang et al.， 2022）。机器学习 （ML） 和深度学习 （DL） 方法 因其 非线性 （Mayoraz 和 Vulliet，2002 年; Zhang et al.， 2021）。 此外，已经提出了一些集成多种算法的改进集成模型，用于针对特色陆地类型以及考虑时间序列的预测（周 et al.， 2016;Han et al.， 2021; Lin et al.， 2022）。然而，如果 数据挖掘方法不 考虑，滑坡物理的

湖号 * 通讯作者：388鲁墨路洪山区武汉市北省邮编：430074， 中国。

电子邮件地址：lichangdong@cug.edu.cn （C. Li）。

https://doi.org/10.1016/j.compgeo.2024.106085

在线2024 年 1 月 16 日提供

0266-352X/ © 2024 爱思唯尔有限公司 版权 所有。

机制。

模拟数据的出色的数据处理能力而性模型的 得到深入研究和广泛应用

的状态变量一，即时间 的概率仅 。 变量的噪声。

预测

可以通过 后验概率 p（x_t—1|y_1:t—1) 在

t —1. 公式 可以 表示如下：

p（x_t|y_1:t—1) = ∫ p（x_t|倍_t—1)p（x_t—1|y_1:t—1)DX系列_t—1 （3）

其中 p（x_t|倍_t—1) 是系统模型的状态 转换概率。

更新过程是使用观测变量 （y） 来校正先验概率密度 p（x|y） 获得后验概率密度 p（x|y） 在时间 t 处。 公式 可以 表示如下：_{tt1:t—1t1:t}

（SAR） 数据 。 Xue et al. （2018） 应用了 DA

p（xt |y t ） = p（y_t|倍_t)p（x_t|y_1:t—1)

(4)

数据的滑坡稳定性分析 仅限于浅研究 基于瞬态降雨渗透和基于网格的区域边坡稳定性 （TRIGRS） 模型以及全球定位系统 （GPS） 和干涉合成孔径雷达 （InSAR） 方法。 然而， （ InSAR 和 SAR 时间序列 数据采集 和处理江 等人，2016年）。 TRIGRS 模型 层、降雨诱发滑坡的（Baum 等人，2008年）。 此外， 应用于滑坡预测的 DA 方法主要关注地表位移，而忽略了深部位移，这更能反映滑坡变形机制（Shentu et al.， 2011）。

针对 具有独特地质结构演化规律和机械强度退化规律的流体动力压力驱动的滑坡，提出了 一种物理信息数据同化（PIDA）方法位移预测。 该文提出PIDA 方法 ，旨在提高物理预报模型的预测能力，解决了数学模型忽略了滑坡动力变形过程的物理影响的问题。

物理信息数据同化 （PIDA） 方法

p（yt |y1：t—1 )

根据 后验概率密度 p（x_t|y_1:t) ，迷你

mum 方差估计的 状态变量计算 如下：

E（x_t) = ∫ x_tp（x_t|y_1:t)DX系列_t （5）

在实际应用中，预测过程中涉及的积分运算难以实现，状态变量的后验概率密度 为 封闭的解析形式 很难表达。Monte Carlo 模拟通过提取一系列遵循状态后验概率分布 的样本，将积分运算问题转化为解决样本期望问题。 它

The PF algorithm is one of the most useful DA methods developed under the framework of the Bayesian estimate and Monte Carlo method,

from the prior probability distribution The posterior probability density

distribution of state variables is expressed as follows:

and Künsch, 2018). Bayesian filtering theory allows for prediction of the prior probability density of a state by utilizing the state transition equation of the system and the posterior probability density of the state from the previous moment (Nakano et al., 2007; Gao and Zhang, 2012; Dahlin and Lindsten, 2014). The state transition rule of a dynamic sys-

。

tem can be expressed by the state equation:

E（x ） = E（x ） = 1 ∑ x^我

(7)

x = f （x， ε） （1）_tt—1t

t t N

t

i=1

And observation equation:

y = g（x， η） （2）_ttt

In the PF method, according to the state empirical distribution of the system, a set of “particles” is randomly generated in the state apace to represent the state variables to be estimated. The PF algorithm effectively estimates the probability density distribution of state variables

and calculates optimal values by utilizing observations to adjust the weight of each particle (Green, 1995; Arulampalam et al., 2002). Theoretically, based on the Strong Law of Large Number, the actual posterior probability density distribution can be accurately simulated as long as the number of particles is large enough.

The realization process of PF algorithm is shown in Fig. 1, and

要选择 重要性概率密度函数， 最小化重要性权重的方差至关重要 。为了简化计算，假设最简单的先验概率密度函数 p（x|x） 作为 重要性概率密度函数，则粒子重要性权重的递归公式可以表示如下：_tt—1

概括如下 （Djuric et al.， 2003）：

Step 1: The sample particles are initialized at first. And the weight of

ωt （xt ） = ωt—1

xt - 1

)p（y_t|倍_t ) （10）

each particle is set as 1/N.

第 3 步：根据观测值更新每个粒子的权重，并对权重值进行归一化。

后验概率分布是未知的，不能用粒子直接模拟。因此，PF 算法采用重要性抽样来解决这个问题 （Chen et al.， 2005）。 已知的重要性概率密度函数 q（x_t|y_1:t) ，则从中提取粒子，然后，实际的后验概率密度 p（x_t|y_1:t) 通过对粒子进行加权来模拟。状态变量的后验概率和最优估计值如下所示：

第 4 步：考虑到 “粒子降解”问题（Doucet et et al.， 2000），重新采样 方法用于保留具有较大权重值的颗粒并消除其他颗粒（Gordon等人，1993 年;Liu 和 Chen，1995 年）。将所有粒子的权重值重置为 1/N，然后开始迭代。

第 5 步：根据 概率密度分布， 计算时间 t 状态变量的最优估计值。

更新 t = t +1 并返回 步骤 2，直到 整个时间为

完成。

Construction of physics-informed landslide displacement prediction model

xiω_t 十一

t ω̃t t )

预测模型如图 1 所示。 2. 首先，根据 深部位移监测数据计算SZD 指数 ， 然后，

= p(y1:t |xt )p(xt ); ω (xi )

last moment and external inducing conditions. The SZD index is taken as

t t t

q(xt |y1:t )

̃_{t t}

the model parameter, which needs to be estimated during the PF algorithm operation. Secondly, the correlation between dynamic shear

i t

i=1

is the normalized importance weight.

strength parameters and the SZD index is constructed, and dynamic estimation of the shear strength parameters is updated accordingly.

无花果。 1. 粒子滤波 （PF） 算法的流程图和 重新采样的过程 。

Fig. 2. The construction of physics-informed landslide displacement prediction model.

Then, through mechanical analysis of landslide deformation and determination of the phreatic line, the physical mechanism model is constructed. Finally, the prediction model is constructed as the observation equation based on the internal force distribution and external deformation monitoring datasets. Considering the landslide dynamic evolution process, the differences between observed and predicted values of real-time surface displacement data are used to update the parameters of the state equation with the help of the PF algorithm, and realizes model prediction trajectory modification.

滑带劣化 （SZD） 模型

监测观察明确 表明，滑坡失败的特点是潜在滑动带的破坏（Li et al.， 2021a）。 提出SZD指数 是为了 反映滑坡滑带的结构完整性和对外部影响条件的动力响应。首先，应通过钻探研究和深位移曲线来确定滑动带位置，因为滑动带附近总是出现突变（Juang，2021 年;Li et al.， 2021b）。

然后，可以通过深位移-深度曲线计算 SZD 指数（图 D）。 3） （Li et al.， 2008）。 钻孔倾角计是

经常用于 深位移测量（Segalini et al.， 2014;Ha et al.， 2018）。累积变形 （W_j） 的钻井倾斜仪监测点 j（j = 0,1,2⋯n，其中 n 是指 数

深部位移监测站点）是指 水平位移

mation differences at time t relative to the deformation value of the fixed monitoring point (j = 0) at the initial moment. Due to the restriction of monitoring technology, only limited sensors are installed to obtain several deep displacement values at different depths, and four representative sensors are selected to quantitatively measure the damage degree of the sliding zone. As shown in Fig. 3a, sensor 3 is usually installed at the top of the drilling hole, which is at the same time above the upper sliding surface. Meanwhile, sensor 0 is a fixed point and installed at the bottom of the drilling hole, which is below the lower sliding surface. Sensors 1 and 2 are both installed in the slip zone, the deep displacement monitoring values of which are generally larger because the landslide sliding zone is prone to relatively more severe shear failure than that of the sliding mass and bedrock.

To simplify the calculation for the SZD index, the deep dis-placement–depth curve is generalized to a broken line (Fig. 3b), which can better reflect the geometric information of shear deformation in the sliding zone. As shown in Fig. 3c, the landslide deformation geometric

无花果。 3. 滑坡深位移-深度曲线 示意图 。

模式中，在滑坡动力演化过程中，相对减振值 ΔW₁₂ (ΔW₁₂ = W₂ —W₁） 的 Grad 逐渐增加，而 相对变形值 ΔW₂₃ (ΔW₂₃ = W₃ —W₂) 在滑动区上方的滑移质量中，滑动区相应减少。因此，SZD 指数 （S_t ） 定义如下：

ΔW₁₂

遵循：

c_t = c_f — c_f — c_r) × S_t （15）

在实际条件下，c 、 φ 、 c 和 φ 是 饱和峰内聚力、饱和 峰值内耗角、饱和残差_ffrr

+ ΔW₂₃

(11)

当滑动区位于 phreatic 线。 换句话说 ，

当 ΔW₁₂ = 0、S_t = 0， 这意味着

滑移区， 完整性最高。 边坡内部结构 完整，变形对环境变化不太敏感。当滑动区发生严重剪切破坏时（当 ΔW₂₃ = 0、S_t = 1）， 则滑动区的 完整性最低。 边坡失去稳定性， 滑移区受损。 在这种情况下 ，内部因素对滑坡的发生影响很大，并且边坡对外部诱发因素很敏感，例如 CTGR 地区的降雨和库水位波动（Long et al.， 2020;Yao等人，2022 年）。SZD 索引具有明确的物理含义，S 的取值范围 _t 是 [0,1]。 随着滑动区的 不断破坏，S_t 显示逐渐增加的趋势。

由于 监控设备和技术的限制， 要获得 Long Time Series 监控数据集。 因此，基于 水动力条件下滑带结构退化的响应机制，构建了序贯 SZD 指数预测模型 来确定 由于数据不足而导致的后续 SZD 指数。拟合 PF 算法中的状态函数，随时计算和提高 SZD 指数的预测精度。在外部诱因因素作用下状态函数的一般形式如下所示：

S_t = f_t(S_t—1、u_t) + ε_t （12）

物理机理模型

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

岩体属于自然状态，c_f 、φ_f 、c_r和 φ_r 是相应的自然剪切强度参数。

水线的位置

计算流体动力压力驱动的滑坡的裂缝线对于分析边坡稳定性非常有用，但也非常复杂（Yan et al.， 2010）。 水库水位条件下 水库线位置的简化方程 降雨量可以根据 Boussinesq 非稳态渗流微分方程使用 Laplacian 矩阵及其 逆变换来确定。 基于一些基本假设：a） 含水层是均质的、各向同性的，具有无限的横向延伸，而 含水层是水平的;

湾） 储层两侧 的潜水流 被视为 一维流;c） 平行于坡面的喷发流是由库水位的变化引起的，降雨渗透导致 喷发流 垂直于 坡面; d） 水库水位以 恒定速度增加或减少 v₀; e） 水库岸被认为是垂直斜坡（Sun et al.， 2018），潜水非定常运动的微分方程可以表示如下：

其中 μ 是 比产; π 是 降水强度; h 和 x 是坐标的垂直和水平方向; t 是时间; T 是含水层的传递系数，T = Kh，其中 K 是渗透系数。

在滑坡演化过程中，滑坡的物质组成和地质结构不断变化 更改 （Xue et al.，

2020 年）， 岩石和土壤的参数也 随

∂h

∂²h

∂x2

π

+ μ （18）

外部影响条件的变化 （Springman et al.， 2003; Zuo et al.， 2020）。 滑带抗剪强度参数 的计算 对 稳定性影响 很大 山体滑坡的分析和预测。 但是， 几乎不可能 获得 实际值 通过现场或室内实验进行数值模拟。实际上，从初始状态到完全剪切破坏状态的过程被认为是滑带抗剪强度 逐渐降低的宏观结果。 因此， 抗剪强度参数的动态变化可以表示为 SZD 指数的函数。

假设 S_t 是 SZD 索引，c_t 和 φ_t 是时间 t 处对应的内聚力和内摩擦角。与 SZD 指数相关的剪切强度参数动态更新模型如下所示

其中 σ = T = KHm ，这被认为是 压力传导率，并且

当水库水位升高时（图 D）。 4a）， 位置 x 在时间 t 的 phreatic 线 的变化 可以表示如下：

u（x， t） = h（x， t） — h（0， 0） （19）_向上向上

其中 h（0， 0） 是指 不透水层与储水库水位之间的 距离，位于 初始时间， h（x， t） 是 _向上

时间 t 位置 x 处 的 phreatic 线。

遵循：

⎧⎪ ∂h = σ

∂h²

∂x2

π

⎪⎨ u（0， t） = h（0， t） — h（0， 0） = h（0， 0） + v t — h（0， 0） = v t， t > 0

⎪ u（x， 0） = h（x， 0） — h（0， 0） = h（0， 0） — h（0， 0） = 0， 0 < x < ∞

t f S_r — S_f

⎩⎪ u（∞， t） = h（∞， t） — h（0， 0） = h（0， 0） + πt — h（0， 0） = 0， t > 0

其中 v₀ 表示水库水上下移动的恒定速度。

当水库水位下降时（图 D）。 4b）， 位置 x 在时间 t 的 phreatic 线 的变化 可以表示如下：

无花果。 4. 水库水位上升时 滑坡水流线计算模型 （a） 和 减少 （b）; （c） 切片法和残差推力法的示意图 。

u（x， t） = h（0， 0） — h（x， t） （21）_下下

The unsteady seepage flow mathematical model of the landslide phreatic line is calculated as follows:

according to the equations mentioned above.

迭代计算

边坡稳定性一直是岩土工程中非常重要的课题（Liao and Ji， 2021; Zhang et al.， 2022）。 但是，

∂²h

∂x2

π

总体安全系数和失稳性分析很难反映滑坡的 DY Namic 演化过程（Miao et al.， 2017）。 限制

基于 Mohr-Coulomb 准则的平衡法 （LEM） 是扩展的

在应用拉普拉斯变换及其逆变换求解数学模型后，使用回归分析来拟合系数（Zheng et al.， 2004）。因此，当储层水位升高和降低时 ，水滴线的解析解可以表示如下：

rock， 而 slice 方法是 实现 极限平衡的 一种方法

载玻片质量分析（Xu等人，2021 年; Su et al.， 2022）。 采用残余推力法来分析物理机制，因为它很容易 实现（Song 和 Xu，2012 年; Zou等人，2020a）。 通过 滑坡结构特性分析，可以得到任意形状的滑面在复杂荷载作用下的滑坡推力和残余推力计算。将滑动体大致均匀地分成 M 个切片， 并在 i 上 进行 力分析 _th 切片 （图 D） 4c）。 考虑到 地下水的影响 ，计算

⎪⎧ 高 x 吨

(v π） a λ⁴ 一个 λ³ 一个 λ² 一个 λ

)t π t h

公式如下：

1 + — （0， 0)

F = [（W + W ）sinα

+ Di cos（αi — β ）]— c_我l_我 + [(W_1我 + W_2我)COSα_我 — D_我sin（α_我 — β_我)]棕褐色φ_我 + Fi ψ

其中 λ =

(23)

(24)

拟合系数分别在水库水位增加和减少折痕时 进行拟合。 可以 确定每个月的温度 线

S

其中 F_我 是切片 i 的残余滑动力; F_s 是稳定性因子。

W_1我 和 W_2我 分别是层 I 的 phreatic 线上方的重力和 phreatic 线以下的浮重; α_我 是切片 i 的滑动面倾角线; c_我 和 φ_我 分别是切片 I 的内聚力和内摩擦角。 D_我 = γ_w一个_我罪恶β_我 是 OSMotic 压力，γ_w 是水的体积密度，A _我 是切片 I 在裂缝线下方的面积 。β_我 是水力坡度方向与水平方向之间的夹角。 ψ_我 是残余推力的传递系数。

方程的右侧。 （24） 由切片 I 的滑动力、抗滑力和残余推力组成。采用迭代计算求解方程。 开始时， 将第一个切片的 残余推力设置为 0，然后计算每个滑块的残余滑动力。 当最后一个切片的 残余推力为 0 时， 迭代将停止，此时每个切片的残余推力可以反映滑坡的受力变形。

推力-位移模型

在总结大量实验结果和现场调查的基础上，Lu （2015） 提出了一种新的岩土节理本构模型，用于描述不同法向应力下滑动面剪切应力-剪切应变的全过程特征 （Sun et al.， 2022）。对于特定的岩石和土壤块，其机械行为仅随法向应力而变化 。 在这种情况下 ，有四个参数可以描述其机械行为。根据这些特性，提出了以下本构等式：

τ = Gγ[1 + γ^m/S]ρ （26）

其中 τ 表示剪切应力，γ 表示剪切应变，G 表示初始剪切模量，m、s、ρ 表示不同法向应力下的常数系数。

与严格理论推导的物理模型不同，上述本构模型是引用数学拟合过程的数学模型。它可以拟合岩石和土体的剪切应力--应变关系。此外，本构模型也可以用于拟合所有物理或化学过程，只要它们能够满足函数曲线特性即可。本文构建包括力和监测数据的方程，将剪切应变替换为残余推力，将剪应力替换为 滑坡月累积位移。 该模型改进如下：

F = GY[1 + Y^m/s]ρ （27）

其中 F 表示残余推力，Y 表示月位移，G 表示初始剪切模量，m、s、ρ 表示 不同的法向应力。

因此，滑坡逆冲和地表监测数据的观测方程一般可以表示如下：

Y_t = g_t(F_t) + η_t （28）

其中 Y_t 是指在时间 t、F 处观测到的累积表面位移 _t 是时间 t 处滑动带 的 推力， g_t（•） 是 反映 力分布和表面位移之间 关系的数据驱动函数 。

物理信息滑坡位移预测模型的步骤 如下所示：

SZD 指数首先 根据方程根据 深位移数据计算。（11）。然后，通过 方程（12）中的 SZD 指数回归模型计算状态方程 的初始化参数。 也可以得出没有深位移监测数据的后续 SZD 指数。

动态剪切强度参数根据 SZD 指数和剪切强度参数之间的关系进行更新 ，并与方程一致。（15） 和 （16），这是状态和观测方程之间的重要联系。

为了进一步分析湿润动力学条件影响下的边坡稳定性， 可以根据方程分别确定滑坡喷发线的位置 。（23） 当水库水位上升和下降时。

滑动面上方的 土壤被 分成几个垂直切片。 根据 切片和滑动面的 phreatic 线的 相对位置， 剪切强度参数 分为 滑动面 高于 切片的 phreatic 线和 Saturated 状态（当它低于 phreatic 线时）。

每个切片的残余推力以 Eqs 为单位获得。 （24） 和 （25） 通过残差推力法的迭代分析，然后，可以通过推力-位移回归模型 Eq 获得预测的位移值。 （28） 称为观测方程。

基于 PF 算法的 PIDA 方法的流程

PIDA 方法结合了物理知情滑坡位移预测模型和 PF 算法，以提供更好的地质参数估计。基于 PF 算法的 PIDA 方法的步骤说明如下（图 1）：

在 PF 过程中，SZD 指数被视为 颗粒。 首先， SZD 指数的初始 值 （S_t) 以平均值为 平均值，初始化一组满足高斯分布的粒子，使这些粒子在 SZD 指数范围内具有较宽的分布。