CN112288169A - 二维滑坡动力学参数概率反分析与滑程超越概率评估方法 - Google Patents

二维滑坡动力学参数概率反分析与滑程超越概率评估方法 Download PDF

Info

Publication number
CN112288169A
CN112288169A CN202011192379.2A CN202011192379A CN112288169A CN 112288169 A CN112288169 A CN 112288169A CN 202011192379 A CN202011192379 A CN 202011192379A CN 112288169 A CN112288169 A CN 112288169A
Authority
CN
China
Prior art keywords
landslide
model
probability
past
distribution
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011192379.2A
Other languages
English (en)
Other versions
CN112288169B (zh
Inventor
曾鹏
孙小平
李天斌
许强
王升
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chengdu Univeristy of Technology
Original Assignee
Chengdu Univeristy of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chengdu Univeristy of Technology filed Critical Chengdu Univeristy of Technology
Priority to CN202011192379.2A priority Critical patent/CN112288169B/zh
Publication of CN112288169A publication Critical patent/CN112288169A/zh
Application granted granted Critical
Publication of CN112288169B publication Critical patent/CN112288169B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0635—Risk analysis of enterprise or organisation activities

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Strategic Management (AREA)
  • Economics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Development Economics (AREA)
  • General Business, Economics & Management (AREA)
  • Tourism & Hospitality (AREA)
  • Quality & Reliability (AREA)
  • Operations Research (AREA)
  • Marketing (AREA)
  • Game Theory and Decision Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Probability & Statistics with Applications (AREA)
  • Artificial Intelligence (AREA)
  • Software Systems (AREA)
  • Educational Administration (AREA)
  • Mathematical Physics (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Data Mining & Analysis (AREA)
  • Pit Excavations, Shoring, Fill Or Stabilisation Of Slopes (AREA)
  • Testing Or Calibration Of Command Recording Devices (AREA)

Abstract

本发明实施例提供了二维滑坡动力学参数概率反分析与滑程超越概率评估方法，涉及滑坡危险性预测技术领域。该方法将动态数值模型的输入参数作为服从某种分布的随机变量，即先验分布，并通过基于马尔可夫链蒙特卡罗的贝叶斯推理，利用过去事件的运动距离作为观察特征，从而有效地改善其先验分布，得到后验分布。随后，如果一个潜在的滑坡与这个过去滑坡相似，那么即可利用改进的后验分布进行潜在滑坡的危险性率预测，估计了滑坡路径上不同位置点的超越概率，获得运动距离超越曲线。本发明实施例通过选取不同的超越概率阈值，从而生成潜在滑坡的运动危险性区划图，为二维滑坡的风险评价和管理提供了重要依据。

Description

二维滑坡动力学参数概率反分析与滑程超越概率评估方法
技术领域
本发明实施例涉及滑坡危险性预测技术领域，特别是涉及二维滑坡动力学参数概率反分析与滑程超越概率评估方法。
背景技术
滑坡是最大的自然灾害之一，在全世界范围内造成许多人死亡和严重的经济损失。为了减少滑坡的社会影响，有必要对二维潜在滑坡进行风险评估并实施缓解策略。运移强度(例如，运动距离、堆积厚度、速度和冲击压力)的预测是风险评估的关键步骤，尤其是涉及极其快速的流状滑坡。
基于物理的动态数值模型已成为一种越来越受欢迎的流状滑坡运动分析工具，包括过去事件的反向分析和潜在未来事件的预测。对于滑坡定量风险评估，运动预测是必要的；但是，由于大多数这些模型使用简化的流变关系来描述二维滑坡复杂的运动过程，因此无法直接和准确地测量模型输入参数，导致模型存在不确定性，为了全面的进行定量风险评估，有必须考虑并量化这些不确定性。
为了评估这些不确定性的影响，一些研究者将动态数值模型中的输入参数视为随机变量，随后采用MCS(蒙特卡罗模拟，Monte Carlo method)进行概率运移预测(也称概率运动距离预测)。例如，使用950MCSs获得关于芬兰峡湾滑坡的流出距离、最大沉积深度和最大速度的统计信息；或采用5000MCSs估算泥石流多个控制点的最大速度和最大高度的概率分布；或应用MCS评估冲积扇各点冲击力的概率密度分布。不幸的是，概率运移预测的准确性很大程度上依赖于不确定参数的统计信息，而这在实验室或现场是很难获得的。
如何选择输入参数以将这些模型用于可靠的滑坡运动预测问题是当前最大的研究挑战之一。通常，可以通过对(相对)有据可查的过去事件进行反向分析来获得合理的输入预测参数。对于在潜在滑坡地区中已经发生了类似的且有据可查的过去滑坡，那么就可通过对过去滑坡的反分析校准确定模型的输入参数，从而对潜在滑坡进行正向预测，即从反向分析到运动预测。
概率反分析的方法可以分为两类：(i)确定性反分析；(ii)概率反分析。确定性反分析旨在确定一组唯一的参数，这些参数将通过反复试验计算出的运动距离、速度或堆积厚度与给定的真实事件相匹配。然而，最佳拟合参数集通常不是唯一的，因为多个参数集可能产生高度相似的模型输出；因此，它们忽略了参数的非唯一性。相反，概率反分析方法能够通过考虑输入参数的不确定性来确定多组参数。这对于概率运移预测是有用的，因为概率反分析的结果可以告知校准参数的不确定性。
最大似然法和贝叶斯方法可能是最常用的概率反分析方法。使用最大似然法对滑坡稳定性问题进行反分析，具有良好的计算效率。尽管概率反分析方法比确定性方法具有多个优势，但很少有研究集中在滑坡运动过程的概率反分析上。一些研究者提出了雪崩危害建模的贝叶斯概率框架，其中动态雪崩模型输入参数的后验分布是通过贝叶斯推论使用先前雪崩的观测结果确定的。还有一些研究者运用贝叶斯定律校准了滑坡运动距离模型并评估了滑坡参数的后验分布。但是，校准后的模型预测未来事件的性能仍然不清楚。因此，需要一种有效、可靠的概率过程来对滑坡运动过程进行反分析，以便能够对滑坡失稳时的运动距离(运动距离为滑坡顶部与堆积坡脚之间的距离)进行可靠的概率正向预测。
发明内容
本发明实施例提供了二维滑坡动力学参数概率反分析与滑程超越概率评估方法，以克服上述的问题。
为了解决上述问题，本发明实施例公开了二维滑坡动力学参数概率反分析与滑程超越概率评估方法，所述方法包括：
分别构建潜在滑坡的克里金代理模型和过去滑坡的克里金代理模型，其中，所述过去滑坡与所述潜在滑坡的物质性质和滑坡类型相似；
定义随机变量，确定所述随机变量的先验分布，并基于马尔可夫链蒙特卡罗的贝叶斯定理，利用所述过去滑坡的克里金代理模型对所述随机变量的先验分布进行更新，得到所述随机变量的后验分布；
利用所述潜在滑坡的克里金代理模型对所述后验分布中的T个样本点分别进行计算，得到每个样本点的运动距离，其中，所述T是蒙特卡罗模拟的样本数；
针对所述T个样本点，建立相应样本点的功能函数，以计算每个样本点的运动距离超过预先设定的运动距离阈值的超越概率，针对多个不同的所述运动距离阈值重复计算所述超越概率，生成运动距离超越概率曲线；
针对所述运动距离超越概率曲线，根据预设的概率分类标准评估所述潜在滑坡的不同区域的危险等级。
在本发明一实施例中，分别构建潜在滑坡的克里金代理模型和过去滑坡的克里金代理模型，包括：
分别构建所述潜在滑坡的动态数值模型和所述过去滑坡的动态数值模型；
在预先设定的上限和下限范围内，通过均匀设计UD方法对所述随机变量进行采样，得到输入数据集，所述输入数据集包括多个采样点；
利用所述潜在滑坡的动态数值模型对每个所述采样点对应的运动距离进行计算，得到第一输出数据集，所述第一输出数据集包括多个所述采样点的运动距离；
利用所述过去滑坡的动态数值模型对每个所述采样点对应的运动距离进行计算，得到第二输出数据集，所述第二输出数据集包括多个所述采样点的运动距离；
基于所述输入数据集和所述第一输出数据集，对克里金Kriging模型进行训练，得到所述潜在滑坡的克里金代理模型；
基于所述输入数据集和所述第二输出数据集，对克里金Kriging模型进行训练，得到所述过去滑坡的克里金代理模型。
在本发明一实施例中，分别构建所述潜在滑坡的动态数值模型和所述过去滑坡的动态数值模型，包括：
分别获得所述潜在滑坡与所述过去滑坡滑动前的数字高程模型和滑动后的数字高程模型；
基于所述过去滑坡滑动前的数字高程模型和滑动后的数字高程模型，建立所述过去滑坡的动态数值模型；
基于所述潜在滑坡滑动前的数字高程模型和滑动后的数字高程模型，建立所述潜在滑坡的动态数值模型。
在本发明一实施例中，所述方法还包括：
在所述输入数据集中随机选择多个附加采样点；
针对每个所述附加采样点，使用动态数值模型和克里金代理模型分别计算所述附加采样点的运动距离；其中，当所述动态数值模型为所述过去滑坡的动态数值模型时，所述克里金代理模型为所述过去滑坡的克里金代理模型；当所述动态数值模型为所述潜在滑坡的动态数值模型时，所述克里金代理模型为所述潜在滑坡的克里金代理模型；
将所述动态数值模型计算获得的运动距离与所述克里金代理模型计算获得的运动距离进行比较，以对所述克里金代理模型的替代性进行评价。
在本发明一实施例中，所述方法还包括：
确定模型偏差系数ε；
滑坡运移距离模型表示为：
y＝ε·D(θ) (1)；
其中，y表示所述过去滑坡的实际运动距离，θ表示所述随机变量的样本点，D(θ)表示由校准前的所述过去滑坡的克里金代理模型计算获得的运动距离；D(θ)通过以下公式计算：
D(θ)＝F(θ)+z(θ) (2)；
其中，F(θ)是通过回归分析获得的趋势函数，z(θ)是一个均值为零的平稳高斯过程；
假设ε服从正态分布，平均值为1，标准偏差为σε，对所述过去滑坡的克里金代理模型进行校准，得到y与D(θ)之间的表达式：
L(θ|y＝Dobs)＝N[Dobs/D(θ)] (3)；
其中，Dobs表示由校准后的所述过去滑坡的克里金代理模型计算获得的运动距离，N为标准的概率密度函数PDF，以平均值1和标准偏差σε来表征；
利用所述过去滑坡的克里金代理模型对所述随机变量的先验分布进行更新，得到所述随机变量的后验分布，包括：
根据所述随机变量的先验分布fprior(θ)，确定一个起始值θ(t)，设置t＝1；
令t＝t+1，从建议分布J(θ*|θ(t-1))生成提案θ*；其中，所述建议分布为多元正态分布，均值为θ(t-1)，协方差矩阵为ξ·Cθ，其中ξ是比例因子，Cθ是θ的先验分布的协方差矩阵；
使用下述公式计算密度比r：
Figure BDA0002753132430000041
其中，q(θ|y)是得到的非标准化后验概率密度，表达式如下：
q(θ|y)＝N[Dobs/D(θ)]·fprior(θ) (5)；
从均匀分布(0,1)生成随机数U，判断U是否≤min(1,r)；
如果U≤min(1,r)，则接受建议并设置θ(t)＝θ*；否则，拒绝建议并设置θ(t)＝θ(t-1)；
迭代上述过程，直到t＝T，根据贝叶斯定理，得到所述随机变量的后验分布fpost(θ|y＝Dobs)：
fpost(θ|y＝Dobs)＝m·N[Dobs/D(θ)]·fprior(θ) (6)；
其中，m是使所述PDF有效的归一化常数。
在本发明一实施例中，针对所述T个样本点，建立相应样本点的功能函数，以计算每个样本点的运动距离超过预先设定的运动距离阈值的超越概率，针对多个不同的所述运动距离阈值重复计算所述超越概率，生成运动距离超越概率曲线，包括：
针对所述T个样本点，将每个样本点带入以下功能函数G(θ)：
G(θ)＝Dtv-D(θ)' (7)；
其中，θ表示所述随机变量的样本点，D(θ)'是所述潜在滑坡的克里金代理模型对应所述θ计算得到的运动距离，Dtv是预先设定的滑坡运动距离阈值；
在G(θ)≤0时，计算每个样本点的D(θ)'超过Dtv的超越概率Pe：
Pe＝P[G(θ)≤0]＝∫G(θ)≤0fpost(θ|Fobs)dθ (8)；
其中，fpost(θ|Fobs)与fpost(θ|y＝Dobs)等同；
用多个不同的所述Dtv带入公式(7)、(8)中进行计算，获得每个样本点的多个所述Pe；
基于所述T个样本点的多个所述Pe，生成运动距离超越概率曲线。
在本发明一实施例中，所述概率分类标准包括：
超越概率≥50.0％，危险等级为极高；
10.0％≤超越概率<50.0％，危险等级为高；
1.0％≤超越概率<10.0％，危险等级为中；
0.1％≤超越概率<1.0％，危险等级为低；
超越概率<0.1％，危险等级为极低。
与现有技术相比，本发明实施例包括以下优点：
本发明实施例将概率反分析与概率运动预测相结合，对过去滑坡事件进行概率反分析，然后将反分析结果用于潜在滑坡的概率运动预测，具体在动态数值模型的输入参数中选择随机变量，确定所述随机变量的先验分布，并基于马尔可夫链蒙特卡罗的贝叶斯定理，利用所述过去滑坡的克里金代理模型对所述随机变量的先验分布进行更新，得到所述随机变量的后验分布；再利用所述潜在滑坡的克里金代理模型对所述后验分布中的蒙特卡罗模拟的样本点分别进行计算，得到每个样本点的运动距离；并针对不同的样本点和不同的运动距离阈值，建立相应样本点的功能函数，以计算每个样本点的运动距离超过该运动距离阈值的超越概率，从而生成运动距离超越概率曲线；最终针对所述运动距离超越概率曲线，根据预设的概率分类标准评估所述潜在滑坡的不同区域的危险等级，达到指导二维滑坡风险评价和风险管理的目的。
本发明实施例通过试验和误差来校准一组最佳参数，该方法考虑了嵌入在输入参数中的不确定性，并得到其后验分布，以此校准了输入参数中不确定性的信息，可以直接用于概率运移预测。
本发明实施例利用UD建立了一个克里金代理模型去近似动态数值模型，代理模型的计算成本低，可以有效地进行后续的概率反分析和概率预测，大大提高了计算效率。
附图说明
图1是本发明实施例提出的滑坡失稳运动模型的一般框架示意图；
图2是本发明实施例二维滑坡动力学参数概率反分析与滑程超越概率评估方法的步骤流程图；
图3是本发明实施例图2的实施程序图；
图4a是黑方台阶地位置示意图；
图4b是DC#2和DC#3发生滑坡后的卫星图；
图5是DC#2滑坡和DC#3滑坡滑动前和滑动后的数字高程模型示意图；
图6是DC#2和DC#3滑坡的动态数值模型示意图；
图7是拟合分布下有效摩擦角
Figure BDA0002753132430000061
的概率密度直方图；
图8是基于克里金代理模型计算的DC#2滑坡的运动距离与基于动态数值模型计算的DC#2滑坡的运动距离的比较示意图；
图9是接受率与比例因子ξ之间的关系示意图；
图10是ξ＝1.5时马尔可夫链样本的性能示意图；
图11a是有效摩擦角
Figure BDA0002753132430000062
的正态分布和毛刺分布示意图；
图11b是孔隙压力比ru的正态分布和毛刺分布示意图；
图12是基于克里金代理模型计算的DC#3滑坡的运动距离与基于动态数值模型计算的DC#3滑坡的运动距离的比较示意图；
图13是拟合分布下DC#3滑坡运动距离的概率密度直方图；
图14是基于超越概率曲线的DC#3滑坡的危险区划示意图；
图15a是模型偏差系数的标准差对随机变量后验均值的影响示意图；
图15b是模型偏差系数的标准差对随机变量后验标准差的影响示意图。
具体实施方式
为使本发明实施例的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明实施例作进一步详细的说明。
已有技术中，过去滑坡的概率运动反分析是通过以下方法实现的：使用(i)动态数值模型计算出的运动距离、堆积厚度或速度以及(ii)滑坡的观察特征Fobs对过去滑坡进行反向分析。使用反向分析校准动态数值模型时，通常会从工程判断，实验室测试或先前的反向分析中得知一些有关输入参数可能值的先验信息。该先验信息可以通过具有概率密度函数(PDF)的先验分布fprior(θ)来描述。随后，贝叶斯定理使用观测值Fobs更新输入参数θ的先验分布，更新得到具有概率密度函数的后验分布fpost(θ|Fobs)。贝叶斯定理的一般形式是：
fpost(θ|Fobs)∝L(θ|Fobs)fprior(θ)；其中，L(θ|Fobs)是似然函数。
然后通过使用校准的输入参数进行动态数值模拟，计算出运移强度(例如，运动距离，堆积厚度以及速度)，可以实现潜在滑坡的运移预测。
但已有技术中滑坡失稳动态数值模型的输入参数不能精确测量，导致预测不确定性较大。为此，本发明实施例提出了二维滑坡动力学参数概率反分析与滑程超越概率评估方法来标定这些参数，并量化它们的不确定性，然后将标定的参数用于前向预测滑坡运动概率。与试错法不同，本文提出的反分析方法采用概率方法，将动态数值模型的输入参数作为服从某种分布的随机变量(先验分布)，并通过基于马尔可夫链蒙特卡罗的贝叶斯推理，利用过去事件的运动距离作为观察特征，从而有效地改善其先验分布，得到后验分布。随后，如果一个潜在的滑坡与这个过去滑坡相似，如物质性质和滑坡类型相似，那么即可利用改进的后验分布进行潜在滑坡的危险性率预测，估计了滑坡路径上不同位置点的超越概率，获得运动距离超越曲线。本发明实施例通过选取不同的超越概率阈值，从而生成潜在滑坡的运动危险性区划图，为二维滑坡风险评价和管理提供了重要依据。此外，在实现过程中，本发明实施例采用基于Kriging的代理模型来近似动态数值模型，显著提高了该方法的计算效率。最后，以甘肃黑方台两个连续的二维滑坡为例，验证了该方法的有效性。
参照图1，描述了本发明实施例提出的滑坡失稳运动模型的一般框架——从滑坡动力学参数概率反分析与滑程超越概率预测。在说明所提出的概率框架之前，必须强调的是，由于材料性质和滑坡类型广泛影响滑坡的失稳行为，在对过去事件的参数进行反向分析校准时，只有当潜在事件发生在已经发生过类似记录的过去事件的区域时，才可以直接将这些参数传输到潜在事件的概率预测中。因此，只有满足上述前提，才能建立提出的框架。
参照图2，示出了本发明实施例二维滑坡动力学参数概率反分析与滑程超越概率评估方法的步骤流程图，该方法包括准备阶段和实施阶段，所述方法可以包括以下步骤：
准备阶段：
步骤S201，分别构建潜在滑坡的克里金代理模型和过去滑坡的克里金代理模型，其中，所述过去滑坡与所述潜在滑坡的物质性质和滑坡类型相似；
实施阶段：
步骤S202，定义随机变量，确定所述随机变量的先验分布，并基于马尔可夫链蒙特卡罗的贝叶斯定理，利用所述过去滑坡的克里金代理模型对所述随机变量的先验分布进行更新，得到所述随机变量的后验分布；
步骤S203，利用所述潜在滑坡的克里金代理模型对所述后验分布中的T个样本点分别进行计算，得到每个样本点的运动距离，其中，所述T是蒙特卡罗模拟的样本数；
步骤S204，针对所述T个样本点，建立相应样本点的功能函数，以计算每个样本点的运动距离超过预先设定的运动距离阈值的超越概率，针对多个不同的所述运动距离阈值重复计算所述超越概率，生成运动距离超越概率曲线；
步骤S205，针对所述运动距离超越概率曲线，根据预设的概率分类标准评估所述潜在滑坡的不同区域的危险等级。
由于本发明实施例基于贝叶斯(Bayes)定理，采用马尔可夫链蒙特卡洛{MCMC，MCMC由两个MC组成，即蒙特卡罗方法(Monte Carlo Simulation，简称MC)和马尔科夫链(Markov Chain，也简称MC}模拟，以对过去滑坡的运动概率进行反分析。但MCMC模拟的一个可能的局限性是它需要对动态数值模型进行迭代评估，这在计算上可能是不允许的。为了解决这一问题，本发明实施例采用克里金代理模型来替代传统的动态数值模型，以此可以有效地进行MCMC模拟。
具体而言，在准备阶段，步骤S201可通过以下子步骤构建潜在滑坡和过去滑坡的克里金代理模型：
子步骤S201-1，分别构建所述潜在滑坡的动态数值模型和所述过去滑坡的动态数值模型；子步骤S201-2，在预先设定的上限和下限范围内，通过均匀设计UD(UniformDesign)方法对所述随机变量进行采样，得到输入数据集，所述输入数据集包括多个采样点；子步骤S201-3，利用所述潜在滑坡的动态数值模型对每个所述采样点对应的运动距离进行计算，得到第一输出数据集，所述第一输出数据集包括多个所述采样点的运动距离；子步骤S201-4，利用所述过去滑坡的动态数值模型对每个所述采样点对应的运动距离进行计算，得到第二输出数据集，所述第二输出数据集包括多个所述采样点的运动距离；子步骤S201-5，基于所述输入数据集和所述第一输出数据集，对克里金Kriging模型进行训练，得到所述潜在滑坡的克里金代理模型；子步骤S201-6，基于所述输入数据集和所述第二输出数据集，对克里金Kriging模型进行训练，得到所述过去滑坡的克里金代理模型。
在本发明实施例中，动态数值模型是一种连续介质力学模型，滑坡的失稳行为由深度积分质量守恒方程和动量守恒方程控制。该模型使用基于麦科马克(MacCormack)全变分递减有限差分格式的数值方法来求解方程组，它的特点是计算时间效率高，在时间和空间都具有二阶精度。在动态数值模型中，基底阻力τb在滑坡动态运移模拟中起着关键作用，三种常见的流变学模型用于计算基底阻力：摩擦模型、Volley模型和宾汉模型。
在子步骤S201-1中，分别构建所述潜在滑坡的动态数值模型和所述过去滑坡的动态数值模型，可以通过以下步骤实现：
分别获得所述潜在滑坡与所述过去滑坡滑动前的数字高程模型和滑动后的数字高程模型；基于所述过去滑坡滑动前的数字高程模型和滑动后的数字高程模型，建立所述过去滑坡的动态数值模型；基于所述潜在滑坡滑动前的数字高程模型和滑动后的数字高程模型，建立所述潜在滑坡的动态数值模型。具体构建潜在滑坡和过去滑坡的数字高程模型的方法，本发明实施例在后续详细叙述，在此不多赘述。在本发明实施例中，分别构建所述潜在滑坡的动态数值模型和所述过去滑坡的动态数值模型，用于模拟潜在滑坡和过去滑坡在二维地形上的崩塌行为。
克里金(Kriging)模型是一种基于统计假设的插值方法，关于训练基于克里金模型的代理模型的细节可参照现有技术，在此不多赘述，如采用由Lophaven等人开发的DACE工具箱来构建基于Kriging的代理模型。此外，为了有效地建立基于Kriging的代理模型，开发了MATLAB与动态数值模型之间数据自动传输的接口代码。
实际中，代理模型的准确性不仅与代理模型的类型有关，而且与用于选择采样点的方法有关。目前，随机抽样方法很流行；但是，它不能保证样本在感兴趣区域中的均匀分布。为了确保随机样本充分覆盖整个感兴趣区域，通常需要更多的采样点，这意味着计算成本增加。因此，在本发明实施例中，在动态数值模型的基础上，采用UD方法(一种确保采样点尽可能均匀分布在感兴趣区域的空间填充技术)来构建基于克里金模型的代理模型的采样点，以此训练得到了潜在滑坡的克里金代理模型和过去滑坡的克里金代理模型。具体实现步骤参照子步骤S201-1～子步骤S201-6。
为保证本发明实施例的克里金代理模型替代动态数值模型进行使用的可靠性，具体使用前，在本发明一实施例中，还提出了对该克里金代理模型的替代性进行评价的方法，可以包括以下步骤：
在所述输入数据集中随机选择多个附加采样点；针对每个所述附加采样点，使用动态数值模型和克里金代理模型分别计算所述附加采样点的运动距离；其中，当所述动态数值模型为所述过去滑坡的动态数值模型时，所述克里金代理模型为所述过去滑坡的克里金代理模型；当所述动态数值模型为所述潜在滑坡的动态数值模型时，所述克里金代理模型为所述潜在滑坡的克里金代理模型；将所述动态数值模型计算获得的运动距离与所述克里金代理模型计算获得的运动距离进行比较，以对所述克里金代理模型的替代性进行评价。
在实施阶段中，如步骤S202所述，本发明实施例采用动态数值模型中考虑的输入参数的向量作为随机变量，确定其随机性，并确定其先验分布，如果没有可用的测试数据，则还可以根据专家的判断，规格或文献数据来推断该随机变量的先验分布。
本发明实施例使用运动距离作为滑坡概率运移反分析中的观察特征，其被定义为滑坡顶部与堆积坡脚之间的距离。但是，一些不确定因素可能会影响反分析，如由于实际条件过于简化而导致的模型不确定性，测量不确定性以及相应的模型输入参数通常不符合给定实际滑坡的“真实”参数。因为这些不确定性，在此，本发明实施例用模型偏差系数ε来表征实际运动距离与模型输出得到的运动距离的比值，来表征它们的影响。滑坡运动距离模型表示为：
y＝ε·D(θ) (1)；
其中，y表示所述过去滑坡的实际运动距离，θ表示所述随机变量的样本点，D(θ)表示由校准前的所述过去滑坡的克里金代理模型计算获得的运动距离；D(θ)通过以下公式计算：
D(θ)＝F(θ)+z(θ) (2)；
其中，F(θ)是通过回归分析获得的趋势函数，z(θ)是一个均值为零的平稳高斯过程，协方差由高斯相关函数计算。
假设ε服从正态分布，平均值为1，标准偏差为σε，对所述过去滑坡的克里金代理模型进行校准，得到y与D(θ)之间的表达式：L(θ|y＝Dobs)＝N[Dobs/D(θ)] (3)；
其中，Dobs表示由校准后的所述过去滑坡的克里金代理模型计算获得的运动距离，N为标准的概率密度函数PDF，以平均值1和标准偏差σε来表征；
接着，本发明实施例采用MCMC模拟从任意分布中抽取样本，然后对其进行校正，以更好地逼近(并最终收敛到)后验分布。使用MCMC模拟的一个优点是它可以考虑任何类型的先验分布，另一个优点是MCMC样本可以直接用于潜在滑坡的概率运移预测。在本发明实施例中，采用MCMC-Metropolis算法方法对过去滑坡进行运动概率反演分析。Metropolis算法的具体实现过程可以概括如下：
根据所述随机变量的先验分布fprior(θ)，确定一个起始值θ(t)，设置t＝1；起始值θ(t)可以从先验分布中随机选择，也可以为其平均值；
迭代：
令t＝t+1，从建议分布J(θ*|θ(t-1))生成提案θ*；其中，所述建议分布为多元正态分布，均值为θ(t-1)，协方差矩阵为ξ·Cθ，其中ξ是比例因子，Cθ是θ的先验分布的协方差矩阵；
使用下述公式计算密度比r：
Figure BDA0002753132430000111
其中，q(θ|y)是得到的非标准化后验概率密度，表达式如下：
q(θ|y)＝N[Dobs/D(θ)]·fprior(θ) (5)；
从均匀分布(0,1)生成随机数U，判断U是否≤min(1,r)；
如果U≤min(1,r)，则接受建议并设置θ(t)＝θ*；否则，拒绝建议并设置θ(t)＝θ(t-1)；
迭代上述过程，直到t＝T，根据贝叶斯定理，得到所述随机变量的后验分布fpost(θ|y＝Dobs)：
fpost(θ|y＝Dobs)＝m·N[Dobs/D(θ)]·fprior(θ) (6)；
其中，m是使所述PDF有效的归一化常数，T是蒙特卡罗模拟的样本数，即迭代的最大次数。
随机变量的后验分布中有多个样本点，本发明实施例利用所述潜在滑坡的克里金代理模型对所述后验分布中的T个样本点分别进行计算，得到每个样本点的运动距离，计算公式可参照式(1)～式(6)。在这一过程中，由于T非常大，本发明实施例因采用克里金代理模型进行计算，相比直接使用动态数值模型进行计算，极大地节省了计算工作量，节约了计算成本。
接着，执行步骤S204，本发明实施例针对所述T个样本点，将每个样本点带入以下功能函数G(θ)：G(θ)＝Dtv-D(θ)'(7)；
其中，θ表示所述随机变量的样本点，D(θ)'是所述潜在滑坡的克里金代理模型对应所述θ计算得到的观察运动距离，Dtv是预先设定的滑坡运动距离阈值；
在G(θ)≤0时，计算每个样本点的D(θ)'超过Dtv的超越概率Pe：
Pe＝P[G(θ)≤0]＝∫G(θ)≤0fpost(θ|Fobs)dθ (8)；
其中，fpost(θ|Fobs)与fpost(θ|y＝Dobs)等同，具体计算时，将式(8)中的fpost(θ|Fobs)替换为fpost(θ|y＝Dobs)进行计算；
用多个不同的Dtv带入公式(7)、(8)中进行计算，获得每个样本点的多个所述Pe；
基于所述T个样本点的多个所述Pe，生成运动距离超越概率曲线。
在本发明实施例中，超越概率可以反映滑坡路径上某一点受到威胁的可能性。因此，本发明实施例针对所述运动距离超越概率曲线，根据预设的概率分类标准评估所述潜在滑坡的不同区域的危险等级，其中：所述概率分类标准包括：超越概率≥50.0％，危险等级为极高；10.0％≤超越概率<50.0％，危险等级为高；1.0％≤超越概率<10.0％，危险等级为中；0.1％≤超越概率<1.0％，危险等级为低；超越概率<0.1％，危险等级为极低。
本发明实施例采用不同超越概率等级进行潜在滑坡的危险区域划分，有助于识别需要立即采取行动的脆弱(失稳)区域，以此，划分后得到的灾害区划图可用于指导二维滑坡风险评价和风险管理。参照图3，示出了本发明实施例图2的建议实施程序图。
为验证本发明实施例的有效性，下面采用一研究案例详细说明。
该研究案例位于距中国甘肃省兰州市45公里的黑方台阶地(见图4a)，属于干旱黄土地区。由于强烈的构造活动，黄河在侵蚀作用下形成了多级阶地，其中黑方台阶地是第四阶地。岩性剖面可依次分为四层：黄土，粘土，砾石和基岩。顶层是马兰黄土，厚度约25-48m。下层粘土层的厚度为3-19m。粘土层覆盖砾石层，厚度为1-6m。基岩由泥岩和少量砂岩构成。黄土下面的粘土起到不可渗透层的作用。自1960年代以来的长期灌溉已大大提高了黑方台阶地的地下水位。由于特殊的地质条件，阶地边缘发生了不同类型的滑坡。发明人通过对黑方台地区69处滑坡的地质构造和特征的研究，将这些滑坡分为黄土滑坡和黄土基岩滑坡。根据滑坡的二维运移特性，将黄土滑坡进一步分为黄土流滑、黄土滑坡和黄土流。在这些黄土滑坡中，位于当川村附近的DC#2和DC#3黄土滑坡(图4b)分别发生在2015年和2016年。在本研究中将它们视为研究案例。
DC#2和DC#3滑坡的滑动前和滑动后的数字高程模型(DEM)如图5所示。DC#2滑坡的源区沿主滑动方向的长度为150m，宽度175m，面积2×104m2。根据滑前DEM(见图5a，DC#2滑前的DEM示意图)与滑后DEM(见图5b，DC#2滑后的DEM示意图)对比得到DC#2滑坡滑前到滑后高程变化图(见图5c)，源区最大厚度为35m，平均厚度约为20m，总体积估计为4×105m3。滑坡经过一个缓坡地带，冲出距离(L)为730米，造成14栋房屋、3家工厂和大片农田受损。滑坡最大堆积厚度约20m，平均厚度约5m，滑坡顶部与堆积坡脚之间的垂直高差(H)为122m，即DC#2滑坡的流动性指数H/L为0.167，具有明显的长距离滑坡的特征。
DC#3滑坡顶高程1715m，堆积坡脚高程1600m，高差115m，滑坡宽度120m，长度60m，滑距520m，因此，DC#3滑坡流动性指数H/L为0.221，同时具有长距离滑坡的特征。根据DC#3滑前的DEM(见图5d)和DC#3滑后的DEM(见图5e)，得到DC#3滑前到滑后高程变化图(见图5f)，总源区面积约为7.2×103m2，最大厚度约为25m，平均厚度为18m，总体积约为1.3×105m3。滑坡最大堆积厚度为10m，平均堆积厚度为2.7m。
显而易见的，黑方台地区黄土滑坡具有相似的物理性质(如岩性、二维形态特征相同)和滑坡类型(如均具有长距离滑坡的特征)。DC#2滑坡发生在靠近DC#3滑坡的区域，如图4b所示。此外，先前的研究表明，这两个滑坡都是黄土滑坡，并且它们到地下水位的深度相似。因此，可以合理地假设DC#2和DC#3滑坡在物质性质和滑坡类型上的相似性。为了说明本发明实施例的方法，将DC#2滑坡视为过去事件，即过去滑坡，用于概率运移反分析，而DC#3滑坡被假定为潜在的未来事件，即潜在滑坡，用于概率运移预测。通过对DC#2滑坡的概率反分析标定输入参数，对DC#3滑坡进行了概率运移预测。
基于滑动前和滑动后DEMs，本研究建立了DC#2和DC#3滑坡的动态数值模型(见图6)。为了给滑坡的最终沉积提供足够的空间，动态数值模型进行了水平扩展，由若干控制参数组成，包括地形网格分辨率、网格数目、总时间和时间步长。表1中列出了用于模拟中的地形网格分辨率、网格数目、总时间和时间步长的值，这些值由发明人经过反复试验确定。
表1：DC#2和DC#3滑坡的质量流控制参数值
Figure BDA0002753132430000131
以往的研究表明，超孔隙水压力是黑方台滑坡长期存在的主要原因之一，孔隙水压力过高是黑方台滑坡长距离运动的主要原因之一。因此，由于其简单性和基于物理的参数，本研究采用了考虑过超孔隙水压力影响的摩擦模型。基础阻力剪切应力τb表示为：
Figure BDA0002753132430000132
式中，σ＝ρgh为滑坡的基底正应力，ρ为质量密度，h为质量高度，g为重力，
Figure BDA0002753132430000133
为质量有效摩擦角；ru为孔隙压力比，表示滑坡的液化程度。其中，ru是由孔隙水压力除以基底正应力得出的，其值在0到1之间。通过对模型参数的敏感性分析发现密度对运动距离的影响较小，Ρ通常是一常数，可设置为1.5g/cm3。孔隙水压力比ru虽然不能准确地反映滑坡路径上的孔隙水压力分布，但它是一个与复杂的未知孔隙水压力相当的指标。由于有效摩擦角
Figure BDA0002753132430000134
的不确定性和孔隙水压力的未知性，本研究将有效摩擦角
Figure BDA0002753132430000135
和孔隙水压力比ru作为随机变量。
具体实现时，考虑了DC#3滑坡的失稳特性，采用黄土的有效残余摩擦角。鉴于有限的测试数据，不可能拟合出先验分布；因此，进行了广泛的文献综述以获得先验分布。表2给出了甘肃34组黄土样品的有效残余摩擦角，这些样品是通过固结排水或不排水剪切试验获得的。
表2：甘肃黄土的有效残余摩擦角
Figure BDA0002753132430000141
根据文献资料，有效摩擦角
Figure BDA0002753132430000142
的平均值、标准差及变异系数分别为34.60°、3.45°以及0.10。图7显示了拟合分布下有效摩擦角
Figure BDA0002753132430000143
的概率密度直方图(有六个区间)；还包括四种常用概率模型(即正态分布、对数正态分布、伽马分布和威布尔分布)的理论分布以供比较。通过K-S(斯米尔诺夫，Kolmogorov-Smirnov)检验，发现正态分布是更合适的分布。然后将有效摩擦角
Figure BDA0002753132430000144
建模为一个正态分布的随机变量，平均值为34.60°，标准偏差为3.45°。
为了考虑过高的孔隙水压力对滑坡运动距离的影响，还需要有关孔隙压力比ru的信息。尽管孔隙压力比ru是基于物理的，但由于受颗粒扩容和水力扩散系数的影响，在滑动运动过程中对超孔隙水压力的产生或消散是非常复杂。因此，平均孔隙压力比ru的先验分布主要通过文献数据进行评估。根据相关文献，ru的参考范围为0.2-0.8，在滑坡稳定性的概率分析中，假设ru服从正态分布。因此，在这项研究中，假设孔隙压力比ru使用3σ法则，服从，平均值为0.5，标准差为0.1的正太分布。表3总结了用于进行概率运移反分析的随机变量的先验分布。此外，为简单起见，假设有效摩擦角
Figure BDA0002753132430000145
和孔隙压力比ru在统计上是独立的，进行了敏感性分析，以检查先验分布对概率反向分析结果的影响。
表3：随机变量的先验分布
Figure BDA0002753132430000151
在进行后验分析之前，本研究构建了克里金代理模型以有效地逼近动态数值模型。首先，通过UD使用(θ+5σ)的上限和(θ-5σ)的下限生成300个训练数据θj(j＝1,2,…,300和
Figure BDA0002753132430000153
)，得到输入数据集。随后，使用过去滑坡的动态数值模型在每个采样点评估DC#2滑坡的运动距离DDC#2(θj)，以生成训练数据(第二输出数据集)的输出。最后，根据θj和DDC#2(θj)的输入和输出数据集，对克里金模型进行了训练以建立了DC#2滑坡的克里金代理模型。同理，DC#3滑坡的克里金代理模型的训练方法参照DC#2滑坡的克里金代理模型的训练方法。
为了检查克里金代理模型的充分性，本研究从输入数据集中随机选择了50个附加采样点，并使用动态数值模型和克里金代理模型计算了它们相应的运动距离。如图8所示，两个模型的计算结果吻合良好。因此，可以使用克里金代理模型
Figure BDA0002753132430000152
进行概率反分析。
用于后验分析的模型偏差系数的标准偏差σε选择为0.15。此外，选择合适的比例因子ξ是进行有效MCMC仿真的关键。对于DC#2滑坡，接受率与比例因子ξ之间的关系如图9所示，接受率随比例因子ξ的增加而减小。当接受率约为0.2-0.4时，马尔可夫链效率最高。分析发现，当缩放ξ太小或太大时，MCMC仿真的效率较低，而当ξ＝1.5时，相应的接受率为0.28。如图10所示，ξ＝1.5时马尔可夫链样本主动移动，这被认为更有效，发现ξ对马氏链中
Figure BDA0002753132430000154
和ru的影响是一致的。因此，对于DC#2滑坡的概率运移反分析，比例因子ξ可设置为1.5。
本发明实施例在这项研究中，使用了具有100,000个样本的马尔可夫链，并且由于马尔可夫链在开始阶段可能未达到稳定状态，因此丢弃了马尔可夫链中的前500个样本。使用选定比例因子(ξ＝1.5)进行MCMC模拟，有效摩擦角
Figure BDA0002753132430000158
和孔隙压力比ru通过DC#2滑坡观测到的运动距离(Dobs＝730m)进行更新。马尔可夫链样本的概率密度直方图如图11所示。有效摩擦角
Figure BDA0002753132430000159
和孔隙压力比ru的后验分布分别近似服从正态分布和毛刺分布。根据这些马尔可夫链样本，计算出随机变量的后均值和标准差：有效摩擦角
Figure BDA0002753132430000157
分别为32.06和3.26。孔隙压力比ru分别为0.63和0.06。图11比较了随机变量的先验分布和后验分布。结果表明，有效摩擦角
Figure BDA0002753132430000155
的后验均值减小，说明先验分布高估了有效摩擦角
Figure BDA0002753132430000156
孔隙压力比ru的后验均值增加，说明先验分布低估了孔隙压力比ru。此外，更新后有效摩擦角
Figure BDA0002753132430000163
和孔隙压力比ru的标准偏差都减小了，这表明在这种情况下可以通过反分析来减小输入参数的不确定性。
按照前述步骤，本研究基于300个输入的训练数据，构造了一个克里金代理模型
Figure BDA0002753132430000161
来近似DC#3滑坡的动态数值模型DDC#3(θ)。如图12所示，
Figure BDA0002753132430000162
为DC#3滑坡的动态数值模型提供了很好的近似值。随后，使用代理模型计算出对应于每个后验样本的运动距离fpost(θ|Fobs)。图13显示了计算得到的DC#3滑坡运动距离的概率密度直方图，其中包括五个常用的理论分布进行比较。通过K–S检验，选择了广义极值分布作为数据的最佳拟合。发现运动距离的平均值为574.84m，标准偏差为74.22m，在95％置信度下的置信区间为[457m，747m]。DC#3滑坡观测到的运动距离为520m，在计算出的置信区间内，表明预测结果是合理的。
通过构造沿滑坡二维路径上任意点的性能函数，并使用带有100,000个MCMC样本，可以获得沿滑坡路径上任意点的超越概率。然后，如图14所示，绘制了滑坡运动距离的超越概率曲线。根据提出的概率分类标准可知，DC#3滑坡的潜在影响区域分为五类：极高、高、中等、低以及极低，如图14所示。在实际工程中，极高危险区应最令人担忧的，因为影响该区域的可能性大于不影响区域的可能性。灾害分区的结果表明，DC#3滑坡的观测积累区位于极高的危险区，表明灾害分区的结果是合理的，可以在一定程度上指导二维滑坡的风险评估和管理。
还需要说明的是，潜在滑坡的概率运移预测的可靠性不仅取决于与过去滑坡的相似性，还取决于运动概率返回中涉及的一些用户定义参数(例如，模型偏差系数ε和先验分布)。通过对用户定义参数进行了敏感性研究以检查其作用。如公式(1)中引入了模型偏差系数ε。如前所述，模型偏差系数ε的概率分布是未知的，假设它是平均值为1、标准偏差为σε的正态分布。利用当前研究案例中给出的先验分布和10万个样本的马尔可夫链，对DC#2滑坡进行了在不同σε值的概率反分析。图15给出了随机变量σε不同值的后验统计量，表明σε显著影响概率反分析的结果，并且随着σε增加，后验统计量逐渐接近先验统计。这一结果的原因是，在σε较大(即残差较大)的情况下，反分析的结果更多地依赖于先验分布而非观察特征。然而，概率反分析的目的是利用观察特征更新先验分布。因此，如果σε不能准确估计，可使用较小的σε值来获得更依赖于观测值的后验结果。
如前所述，概率运移反分析需要指定随机变量的先验分布。为了检查先验分布对后验统计量的影响，表4给出了四个具有不同假设的先验分布：本研究案例中使用的先验分布表示为先验1；先验2假设随机变量遵循多元正态分布，且均值不同。先验3假设随机变量服从具有不同标准偏差的多元正态分布；先验4假设随机变量遵循与先验1不同的对数分布，且分布参数不同。对DC#2滑坡的概率运移反分析在不同先验分布(先验1～先验4)下重复进行，计算结果列于表5，表明后验统计对先验分布较为敏感，说明先验分布对概率运移反分析的结果有重要影响。
表4：四种不同假设的先验分布
Figure BDA0002753132430000171
表5：先验分布对后验统计量的影响
Figure BDA0002753132430000172
DC#2滑坡的结果表明，概率反分析有效地更新了先验均值，表明它们被低估或高估。为了证明后验更新的作用，利用先验分布和后验分布预测了DC#3滑坡的运动距离，其预测均值分别为374.90和574.84m。与观测到的运动距离(520m)相比，使用先验分布的预测结果证明是不准确的，因为对运动距离的低估会对生命和基础设施造成极大的威胁。然而，使用后验分布的预测提供了更为合理的运动距离预测，并且也有可能通过使用连续贝叶斯更新(从几个类似的滑坡事件中观察到的运动距离)进一步改进该预测。此外，图11显示，通过概率反分析，输入参数的标准偏差减小，表明其不确定性正在降低。
综上，本发明实施例提出的滑坡运动危险性评价的方法，即将概率反分析与概率运动预测相结合。从上述研究案例中看出，本发明实施例将DC#2滑坡作为过去事件进行概率反分析，DC#3滑坡(靠近DC#2滑坡区域)假设为潜在事件进行概率运动预测，利用改进的动态数值模型(即克里金代理模型)来计算运动距离，采用基于Bayes定理和MCMC模拟的概率反分析方法，对DC#2滑坡的模型输入参数进行更新，并将其结果用于DC#3滑坡的概率运移预测。从这项研究中可以得出几个结论：
1)与反复试错的确定反分析方法不同，该方法可以考虑嵌入在输入参数中的不确定性，并得到其后验分布，以此校准了输入参数中不确定性的信息，可以直接用于概率运移预测。
2)利用DC#2滑坡的概率运移反分析结果，对DC#3滑坡进行了概率运移预测，计算出滑坡路径上各点运动距离的超越概率。计算出的超越概率可直接用于单体滑坡的定量风险评估。此外，根据产生的运动距离-超越概率曲线，使用不同的超越概率分级值将受DC#3滑坡影响的区域划分为五个危险区。该危险性区划表明，DC#3滑坡的观测堆积区处于极高危险区，表明结果是合理的，对二维滑坡风险评价和治理具有一定的指导意义。
3)由于需要进行成千上万次的动态数值模拟，因此无论是概率运移反分析还是概率运移预测都需要大量的计算成本。为了解决这一问题，本发明实施例利用UD建立了一个基于Kriging的代理模型(克里金代理模型)去近似动态数值模型。模型精度检验表明，用300个训练数据(即300个动态数值模拟)训练代理模型就足够，训练工作量较低，使用代理模型可以有效地进行后续的概率反分析和概率预测，大大提高了计算效率。
以上对本发明实施例所提供的二维滑坡动力学参数概率反分析与滑程超越概率评估方法进行了详细介绍，本文中应用了具体个例对本发明实施例的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明实施例的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明实施例的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明实施例的限制。

Claims (7)

1.二维滑坡动力学参数概率反分析与滑程超越概率评估方法，其特征在于，所述方法包括：
分别构建潜在滑坡的克里金代理模型和过去滑坡的克里金代理模型，其中，所述过去滑坡与所述潜在滑坡的物质性质和滑坡类型相似；
定义随机变量，确定所述随机变量的先验分布，并基于马尔可夫链蒙特卡罗的贝叶斯定理，利用所述过去滑坡的克里金代理模型对所述随机变量的先验分布进行更新，得到所述随机变量的后验分布；
利用所述潜在滑坡的克里金代理模型对所述后验分布中的T个样本点分别进行计算，得到每个样本点的运动距离，其中，所述T是蒙特卡罗模拟的样本数；
针对所述T个样本点，建立相应样本点的功能函数，以计算每个样本点的运动距离超过预先设定的运动距离阈值的超越概率，针对多个不同的所述运动距离阈值重复计算所述超越概率，生成运动距离超越概率曲线；
针对所述运动距离超越概率曲线，根据预设的概率分类标准评估所述潜在滑坡的不同区域的危险等级。
2.根据权利要求1所述的方法，其特征在于，分别构建潜在滑坡的克里金代理模型和过去滑坡的克里金代理模型，包括：
分别构建所述潜在滑坡的动态数值模型和所述过去滑坡的动态数值模型；
在预先设定的上限和下限范围内，通过均匀设计UD方法对所述随机变量进行采样，得到输入数据集，所述输入数据集包括多个采样点；
利用所述潜在滑坡的动态数值模型对每个所述采样点对应的运动距离进行计算，得到第一输出数据集，所述第一输出数据集包括多个所述采样点的运动距离；
利用所述过去滑坡的动态数值模型对每个所述采样点对应的运动距离进行计算，得到第二输出数据集，所述第二输出数据集包括多个所述采样点的运动距离；
基于所述输入数据集和所述第一输出数据集，对克里金Kriging模型进行训练，得到所述潜在滑坡的克里金代理模型；
基于所述输入数据集和所述第二输出数据集，对克里金Kriging模型进行训练，得到所述过去滑坡的克里金代理模型。
3.根据权利要求2所述的方法，其特征在于，分别构建所述潜在滑坡的动态数值模型和所述过去滑坡的动态数值模型，包括：
分别获得所述潜在滑坡与所述过去滑坡滑动前的数字高程模型和滑动后的数字高程模型；
基于所述过去滑坡滑动前的数字高程模型和滑动后的数字高程模型，建立所述过去滑坡的动态数值模型；
基于所述潜在滑坡滑动前的数字高程模型和滑动后的数字高程模型，建立所述潜在滑坡的动态数值模型。
4.根据权利要求2所述的方法，其特征在于，所述方法还包括：
在所述输入数据集中随机选择多个附加采样点；
针对每个所述附加采样点，使用动态数值模型和克里金代理模型分别计算所述附加采样点的运动距离；其中，当所述动态数值模型为所述过去滑坡的动态数值模型时，所述克里金代理模型为所述过去滑坡的克里金代理模型；当所述动态数值模型为所述潜在滑坡的动态数值模型时，所述克里金代理模型为所述潜在滑坡的克里金代理模型；
将所述动态数值模型计算获得的运动距离与所述克里金代理模型计算获得的运动距离进行比较，以对所述克里金代理模型的替代性进行评价。
5.根据权利要求1或2所述的方法，其特征在于，所述方法还包括：
确定模型偏差系数ε；
滑坡运移距离模型表示为：
y＝ε·D(θ) (1)；
其中，y表示所述过去滑坡的实际运动距离，θ表示所述随机变量的样本点，D(θ)表示由校准前的所述过去滑坡的克里金代理模型计算获得的运动距离；D(θ)通过以下公式计算：
D(θ)＝F(θ)+z(θ) (2)；
其中，F(θ)是通过回归分析获得的趋势函数，z(θ)是一个均值为零的平稳高斯过程；
假设ε服从正态分布，平均值为1，标准偏差为σε，对所述过去滑坡的克里金代理模型进行校准，得到y与D(θ)之间的表达式：
L(θ|y＝Dobs)＝N[Dobs/D(θ)] (3)；
其中，Dobs表示由校准后的所述过去滑坡的克里金代理模型计算获得的运动距离，N为标准的概率密度函数PDF，以平均值1和标准偏差σε来表征；
利用所述过去滑坡的克里金代理模型对所述随机变量的先验分布进行更新，得到所述随机变量的后验分布，包括：
根据所述随机变量的先验分布fprior(θ)，确定一个起始值θ(t)，设置t＝1；
令t＝t+1，从建议分布J(θ*|θ(t-1))生成提案θ*；其中，所述建议分布为多元正态分布，均值为θ(t-1)，协方差矩阵为ξ·Cθ，其中ξ是比例因子，Cθ是θ的先验分布的协方差矩阵；
使用下述公式计算密度比r：
Figure FDA0002753132420000031
其中，q(θ|y)是得到的非标准化后验概率密度，表达式如下：
q(θ|y)＝N[Dobs/D(θ)]·fprior(θ) (5)；
从均匀分布(0,1)生成随机数U，判断U是否≤min(1,r)；
如果U≤min(1,r)，则接受建议并设置θ(t)＝θ*；否则，拒绝建议并设置θ(t)＝θ(t-1)；
迭代上述过程，直到t＝T，根据贝叶斯定理，得到所述随机变量的后验分布fpost(θ|y＝Dobs)：
fpost(θ|y＝Dobs)＝m·N[Dobs/D(θ)]·fprior(θ) (6)；
其中，m是使所述PDF有效的归一化常数。
6.根据权利要求1所述的方法，其特征在于，针对所述T个样本点，建立相应样本点的功能函数，以计算每个样本点的运动距离超过预先设定的运动距离阈值的超越概率，针对多个不同的所述运动距离阈值重复计算所述超越概率，生成运动距离超越概率曲线，包括：
针对所述T个样本点，将每个样本点带入以下功能函数G(θ)：
G(θ)＝Dtv-D(θ)' (7)；
其中，θ表示所述随机变量的样本点，D(θ)'是所述潜在滑坡的克里金代理模型对应所述θ计算得到的运动距离，Dtv是预先设定的滑坡运动距离阈值；
在G(θ)≤0时，计算每个样本点的D(θ)'超过Dtv的超越概率Pe：
Pe＝P[G(θ)≤0]＝∫G(θ)≤0fpost(θ|Fobs)dθ (8)；
其中，fpost(θ|Fobs)与fpost(θ|y＝Dobs)等同；
用多个不同的所述Dtv带入公式(7)、(8)中进行计算，获得每个样本点的多个所述Pe；
基于所述T个样本点的多个所述Pe，生成运动距离超越概率曲线。
7.根据权利要求1或6所述的方法，其特征在于，所述概率分类标准包括：
超越概率≥50.0％，危险等级为极高；
10.0％≤超越概率<50.0％，危险等级为高；
1.0％≤超越概率<10.0％，危险等级为中；
0.1％≤超越概率<1.0％，危险等级为低；
超越概率<0.1％，危险等级为极低。
CN202011192379.2A 2020-10-30 2020-10-30 二维滑坡动力学参数概率反分析与滑程超越概率评估方法 Active CN112288169B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011192379.2A CN112288169B (zh) 2020-10-30 2020-10-30 二维滑坡动力学参数概率反分析与滑程超越概率评估方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011192379.2A CN112288169B (zh) 2020-10-30 2020-10-30 二维滑坡动力学参数概率反分析与滑程超越概率评估方法

Publications (2)

Publication Number Publication Date
CN112288169A true CN112288169A (zh) 2021-01-29
CN112288169B CN112288169B (zh) 2022-03-18

Family

ID=74353625

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011192379.2A Active CN112288169B (zh) 2020-10-30 2020-10-30 二维滑坡动力学参数概率反分析与滑程超越概率评估方法

Country Status (1)

Country Link
CN (1) CN112288169B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113343504A (zh) * 2021-08-02 2021-09-03 成都理工大学 一种三维滑坡运动危险性概率评价方法
CN113420165A (zh) * 2021-06-11 2021-09-21 北京达佳互联信息技术有限公司 二分类模型的训练、多媒体数据的分类方法及装置
CN114781731A (zh) * 2022-04-26 2022-07-22 成都理工大学 基于贝叶斯理论的滑坡运动距离超越概率预测方法和系统
CN115510578A (zh) * 2022-09-26 2022-12-23 成都理工大学 基于InSAR近实时监测的滑坡失稳时间概率预测方法及产品
CN117172368A (zh) * 2023-09-04 2023-12-05 中国地质大学（武汉） 滑坡-涌浪最大高度超越概率预测方法、设备及存储设备

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110054869A1 (en) * 2008-05-05 2011-03-03 Dachang Li Modeling Dynamic Systems By Visualizing and Narrowing A Parameter Space
CN104820750A (zh) * 2015-05-11 2015-08-05 广西大学 一种基于判别分析的结构可靠度动态响应面方法
CN111445120A (zh) * 2020-03-24 2020-07-24 成都理工大学 一种滑坡运动距离超越概率计算方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110054869A1 (en) * 2008-05-05 2011-03-03 Dachang Li Modeling Dynamic Systems By Visualizing and Narrowing A Parameter Space
CN104820750A (zh) * 2015-05-11 2015-08-05 广西大学 一种基于判别分析的结构可靠度动态响应面方法
CN111445120A (zh) * 2020-03-24 2020-07-24 成都理工大学 一种滑坡运动距离超越概率计算方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
XIAOPING SUN等: "Run-out distance exceedance probability evaluation and hazard zoning of an individual landslide", 《LANDSLIDES》 *
滕斌 等: "基于Kriging模型的土质边坡可靠度计算方法", 《公路工程》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113420165A (zh) * 2021-06-11 2021-09-21 北京达佳互联信息技术有限公司 二分类模型的训练、多媒体数据的分类方法及装置
CN113420165B (zh) * 2021-06-11 2024-03-05 北京达佳互联信息技术有限公司 二分类模型的训练、多媒体数据的分类方法及装置
CN113343504A (zh) * 2021-08-02 2021-09-03 成都理工大学 一种三维滑坡运动危险性概率评价方法
CN114781731A (zh) * 2022-04-26 2022-07-22 成都理工大学 基于贝叶斯理论的滑坡运动距离超越概率预测方法和系统
CN114781731B (zh) * 2022-04-26 2023-04-18 成都理工大学 基于贝叶斯理论的滑坡运动距离超越概率预测方法和系统
CN115510578A (zh) * 2022-09-26 2022-12-23 成都理工大学 基于InSAR近实时监测的滑坡失稳时间概率预测方法及产品
CN115510578B (zh) * 2022-09-26 2023-07-14 成都理工大学 基于InSAR近实时监测的滑坡失稳时间概率预测方法及产品
CN117172368A (zh) * 2023-09-04 2023-12-05 中国地质大学（武汉） 滑坡-涌浪最大高度超越概率预测方法、设备及存储设备
CN117172368B (zh) * 2023-09-04 2024-02-27 中国地质大学（武汉） 滑坡-涌浪最大高度超越概率预测方法、设备及存储设备

Also Published As

Publication number Publication date
CN112288169B (zh) 2022-03-18

Similar Documents

Publication Publication Date Title
CN112288169B (zh) 二维滑坡动力学参数概率反分析与滑程超越概率评估方法
Di Stefano et al. Flow resistance equation for rills
Domeneghetti et al. Assessing rating-curve uncertainty and its effects on hydraulic model calibration
US20090164186A1 (en) Method for determining improved estimates of properties of a model
Sun et al. Data-space approaches for uncertainty quantification of CO2 plume location in geological carbon storage
Sun et al. From probabilistic back analyses to probabilistic run-out predictions of landslides: A case study of Heifangtai terrace, Gansu Province, China
Wang et al. Bayesian stochastic soil modeling framework using Gaussian Markov random fields
CN113610945A (zh) 一种基于混合神经网络的地应力曲线预测方法
CN113343504A (zh) 一种三维滑坡运动危险性概率评价方法
Kouzehgar et al. Physical modeling into outflow hydrographs and breach characteristics of homogeneous earthfill dams failure due to overtopping
Wang et al. Data-driven analysis of soil consolidation with prefabricated vertical drains considering stratigraphic variation
Ghoshal et al. Distribution of sediment concentration in debris flow using Rényi entropy
US11927717B2 (en) Optimized methodology for automatic history matching of a petroleum reservoir model with Ensemble Kalman Filter (EnKF)
Pektaş Determining the essential parameters of bed load and suspended sediment load
Muduli et al. Evaluation of liquefaction potential of soil based on shear wave velocity using multi-gene genetic programming
Kovacevic et al. The use of neural networks to develop CPT correlations for soils in northern Croatia
Spacagna et al. Data-driven soil profile characterization using statistical methods and artificial intelligence algorithms
Kheiry et al. Uncertainty quantification of steady-state seepage through earth-fill dams by random finite element method and multivariate adaptive regression splines
Gómez-Hernández Geostatistics
Barnhart et al. Evaluation of debris-flow building damage forecasts
CN117933103B (zh) 一种基于贝叶斯深度学习的碳封存模型不确定性分析方法
Witty et al. Comparison of Gaussian and Indicator Based Sequential Simulation Methods for 3D Spatial Uncertainty Quantification in Subsoil Modeling Using Cone Penetration Tests
Sheinin Use of the lognormal distribution for processing the results of mechanical testing of soils
Kovacevic et al. Developing correlations between the soil fines content and CPT results using neural networks
Zhou et al. Lateral variations of shear wave velocity (VS) profile and VS30 over gentle terrain

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant