CN113343504A - 一种三维滑坡运动危险性概率评价方法 - Google Patents

一种三维滑坡运动危险性概率评价方法 Download PDF

Info

Publication number
CN113343504A
CN113343504A CN202110878125.4A CN202110878125A CN113343504A CN 113343504 A CN113343504 A CN 113343504A CN 202110878125 A CN202110878125 A CN 202110878125A CN 113343504 A CN113343504 A CN 113343504A
Authority
CN
China
Prior art keywords
landslide
model
motion
maximum
observation point
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.)
Pending
Application number
CN202110878125.4A
Other languages
English (en)
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 CN202110878125.4A priority Critical patent/CN113343504A/zh
Publication of CN113343504A publication Critical patent/CN113343504A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Strategic Management (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Economics (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Development Economics (AREA)
  • Geometry (AREA)
  • Educational Administration (AREA)
  • Evolutionary Computation (AREA)
  • Game Theory and Decision Science (AREA)
  • Computer Hardware Design (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Testing Or Calibration Of Command Recording Devices (AREA)

Abstract

本发明公开了一种三维滑坡运动危险性概率评价方法，涉及滑坡运动危险性预测技术领域。该方法将动态数值模型用于滑坡运动分析，并以最终堆积深度为观测信息，利用基于MCMC模拟的多观测贝叶斯反分析方法对动态数值模型输入参数的不确定性进行校准，得到了后验分布。然后，使用后验分布作为潜在滑坡的输入，以在三维地形上估计潜在滑坡最大运动高度和最大运动速度的超越概率，并生成用于单体滑坡定量风险评价的概率运动危险性图。本发明采用多响应全局克里金代理模型近似模型输入和输出之间的关系，使其完全替代动态数值模型进行贝叶斯反分析和概率运动危险性评估，在不损失精度的前提下大大提高了计算效率。

Description

一种三维滑坡运动危险性概率评价方法
技术领域
本发明涉及滑坡运动危险性预测技术领域，特别是涉及一种三维滑坡运动危险性概率评价方法。
背景技术
山体滑坡造成了严重的经济损失和人员伤亡，但它们的识别、监测和控制措施在技术上或经济上并不总是可行的。因此，考虑灾害情景及其潜在后果的滑坡QRA(Quantitative Risk Analysis，定量风险评估)正成为一个研究热点。为了评估单体滑坡造成的潜在风险，必须定量评估这几个概率，包括：(i)在给定回归期内以某种规模失稳的概率；(ii)滑坡物质以某种运动强度到达某位置的概率，通常称为运动危险性评估；(iii)承灾体在网格上出现在该位置的概率；(iv)承灾体损坏或人员伤亡的概率，亦称易损性。风险R是这些概率和成本C的乘积。网格影响概率
Figure 996947DEST_PATH_IMAGE001
是连接滑坡发生和潜在后果的桥梁；因此，滑坡运动危险性评价是滑坡风险评价的关键问题。
考虑到预测模型中的不确定性，计算
Figure 153122DEST_PATH_IMAGE002
(即运动危险性概率评估)相当于评估滑坡运动强度超越概率(即运动强度超过给定阈值的概率，P e )。这个问题可以用功能函数G(θ)来表示：
Figure 137127DEST_PATH_IMAGE003
；
上式中，
Figure 745963DEST_PATH_IMAGE004
是承灾体的给定阈值；θ 是动态数值模型所涉及的不确定输入变量；
Figure 653876DEST_PATH_IMAGE005
是计算的运动强度。因此，滑坡运动强度超过阈值的概率（即P e ）即为
Figure 664558DEST_PATH_IMAGE006
的概率，并可使用以下概率积分确定：
Figure 632514DEST_PATH_IMAGE007
；
上式中，
Figure 728646DEST_PATH_IMAGE008
是θ的联合概率密度函数。为了对上式中的P e 进行良好的估计，至少应适当解决三个主要问题：(i)在给定输入参数的情况下，一个能合理预测任意位置的运动强度的动态数值模型；(ii)一种有效且准确的可靠性方法，该方法可解释输入参数中的不确定性，并可通过动态数值模型将这些不确定性传播到计算的运动强度中；(iii)利用有效且廉价的方法获得滑坡输入参数可靠的统计信息。
由于基于物理的动态数值模型能够在不规则地形上进行运动分析，并且计算结果可以直接与QRA的易损性耦合，因此它们已越来越多地用于概率运动危险评估。目前，滑坡运动危险性预测通常是确定性的，没有考虑动态数值模型中的不确定性。然而，在滑坡QRA的框架下，考虑不确定性的运动危险概率评估更为合适。
近年来，滑坡运动危险性引起了人们的广泛关注。但目前关于滑坡运动(Run-out)分析的研究主要集中在二维模型上。众所周知，二维模型不能考虑滑坡的横向扩展；因此，它无法提供潜在滑坡的网格影响范围。另外，二维运动预测需要预先假设滑坡的运动路径，建立二维运动模型，这对于预测未来的事件并不总是可行的。因此，有必要对滑坡进行三维运动危险性概率评价。
但相关技术中很少有人尝试进行三维运动危险性概率预测。一方面是成本问题，包括计算成本和研究成本，这是滑坡运动危险性概率评价实际应用的关键；另一方面是输入参数统计信息的可靠性对概率运动预测的性能也至关重要。然而，如何正确地确定动态数值模型中输入参数的统计信息，目前的研究还很有限。
发明内容
本发明实施例提供一种三维滑坡运动危险性概率评价方法，可以克服上述技术问题。
为了解决上述问题，本发明实施例公开了一种三维滑坡运动危险性概率评价方法，所述方法包括：
步骤S1：确定潜在滑坡对应的过去滑坡，并在所述过去滑坡的滑坡堆积区内选择多个观测点；
步骤S2：对所述多个观测点中的每一个观测点构造克里金代理模型，并基于每一个观测点的克里金代理模型、贝叶斯定理以及差分进化适应大都会市算法，对预先确定的随机变量的先验信息进行反分析，得到所述随机变量的后验分布；
步骤S3：将所述潜在滑坡预估的滑坡堆积区网格化，选取所述潜在滑坡整个运动过程中各网格位置的最大运动高度和最大运动速度作为评价指标，建立所述潜在滑坡在每一个网格的最大运动高度功能函数和最大运动速度功能函数，并计算得到所述潜在滑坡的最大运动高度克里金代理模型矩阵和最大运动速度克里金代理模型矩阵；
步骤S4：将所述随机变量的后验分布输入到所述潜在滑坡的最大运动高度克里金代理模型矩阵中，输出得到每一个网格的最大运动高度，以及将所述随机变量的后验分布输入到所述最大运动速度克里金代理模型矩阵中，输出得到每一个网格的最大运动速度，并对每一个网格的最大运动高度超越概率以及最大运动速度超越概率分别进行计算，获得所述潜在滑坡的最大运动高度超越概率矩阵以及所述潜在滑坡的最大运动速度超越概率矩阵；
步骤S5：在三维地理信息系统平台上，根据所述潜在滑坡的最大运动高度超越概率矩阵以及所述潜在滑坡的最大运动速度超越概率矩阵，生成用于单体滑坡定量风险评价的概率运动危险性图。
进一步的，所述多个观测点的个数为K个；所述步骤S2包括：
子步骤S2-1：对所述K个观测点中的每一个观测点构造动态数值模型，并建立每一个观测点的克里金代理模型，将所述动态数值模型的输入参数作为随机变量θ，利用同一个观测点的克里金代理模型计算该观测点处所述随机变量对应的最终堆积深度模拟值，获得所述K个观测点的最终堆积深度模拟值
Figure 377933DEST_PATH_IMAGE009
，其中，
Figure 243121DEST_PATH_IMAGE010
表示第K个观测点的最终堆积深度模拟值；
子步骤S2-2：利用设定的模型偏差系数
Figure 381978DEST_PATH_IMAGE011
和所述K个观测点的最终堆积深度模拟值
Figure 965406DEST_PATH_IMAGE012
，计算得到所述K个观测点的最终堆积深度观测值
Figure 480701DEST_PATH_IMAGE013
，其中，
Figure 200395DEST_PATH_IMAGE014
；
子步骤S2-3：基于所述过去滑坡的K个观测点对应的最终堆积深度观测值，通过贝叶斯定理对所述随机变量的先验分布
Figure 697105DEST_PATH_IMAGE015
进行反分析，利用差分进化适应大都会市算法对反分析结果进行采样，得到所述随机变量的后验分布
Figure 767829DEST_PATH_IMAGE016
；
Figure 821235DEST_PATH_IMAGE017
；
上式中，
Figure 661015DEST_PATH_IMAGE018
是似然函数，用于表征现场测量数据与动态数值模型模拟结果之间的匹配程度，表示为所述随机变量的联合条件概率密度函数；其中，在假设不同网格位置之间的模型偏差系数
Figure 141675DEST_PATH_IMAGE019
是不相关的情况下，所述
Figure 699696DEST_PATH_IMAGE020
表示为：
Figure 228897DEST_PATH_IMAGE021
；
上式中，
Figure 923184DEST_PATH_IMAGE022
表示
Figure 840324DEST_PATH_IMAGE023
的概率密度函数在
Figure 620061DEST_PATH_IMAGE024
对应的概率密度。
进一步的，所述子步骤S2-1中的对所述K个观测点中的每一个观测点构造动态数值模型，并建立每一个观测点的克里金代理模型，包括：
子步骤S2-101：使用摩擦模型计算基底阻力
Figure 15270DEST_PATH_IMAGE025
，即
Figure 829643DEST_PATH_IMAGE026
；
上式中，𝜎 = ρgh，𝜎表示基底正应力，ρ表示滑坡物质的密度，h表示滑坡物质的运动高度，g代表重力；
Figure 652105DEST_PATH_IMAGE027
是有效残余摩擦角；r u 是孔隙压力比系数，反映了基底材料的液化程度，取值范围为0~1；
子步骤S2-102：基于所述摩擦模型，利用质量流模型对所述过去滑坡进行运动模拟，对所述K个观测点中的每一个观测点构造动态数值模型；其中，所述质量流模型是基于三维纳维-斯托克斯方程沿深度积分的连续介质模型；
子步骤S2-103：将所述动态数值模型的输入参数
Figure 106089DEST_PATH_IMAGE028
和r u 定义为所述随机变量
Figure 304989DEST_PATH_IMAGE029
，并通过均匀设计对所述随机变量进行采样，得到输入数据集，所述输入数据集包括多个初始采样点；
子步骤S2-104：将所述输入数据集的每个初始采样点代入每一个观测点的动态数值模型中，以分别输出每一个观测点位置对应的响应值，基于所述输入数据集以及每一个观测点位置对应的响应值，对克里金模型进行训练，得到每一个观测点的初始克里金代理模型；
子步骤S2-105：从候选样本池中选择一个新样本点，其中，所述候选样本池
Figure 973868DEST_PATH_IMAGE030
中的候选样本点通过拉丁超立方体抽样生成；j=1，2，⋯, s；其中s是候选样本点的数量；
子步骤S2-106：将所述新样本点添加到所述输入数据集中，以对每一个观测点的初始克里金代理模型重新训练，持续迭代，直到该观测点的克里金代理模型满足预先设定的收敛准则，得到最终该观测点的克里金代理模型。
进一步的，所述子步骤S2-105中的从候选样本池中选择一个新样本点，包括：
通过最大化质量参数
Figure 967232DEST_PATH_IMAGE031
从候选样本池中选择一个新样本点，其中，
Figure 721561DEST_PATH_IMAGE032
的表达式如下：
Figure 458573DEST_PATH_IMAGE033
；
上式中，
Figure 919641DEST_PATH_IMAGE034
是第j个候选样本点与当前设计样本点之间的最小欧氏距离，其中，当前设计样本点包括所述初始采样点和所述新样本点；
Figure 349486DEST_PATH_IMAGE035
是第j个候选样本点处的预测方差，由每一个观测点的第i个克里金代理模型进行估计获得；
w i 为权重因子，表示每一个观测点的第i个克里金代理模型对总预测方差的贡献，其中：
Figure 591111DEST_PATH_IMAGE036
；
上式中，
Figure 866235DEST_PATH_IMAGE037
是每一个观测点的第i个克里金代理模型的预测方差的平均值，
Figure 509706DEST_PATH_IMAGE038
是每一个观测点的第i个克里金代理模型的最大预测方差，p是每一个观测点的克里金代理模型的数目。
进一步的，所述子步骤S2-106包括：
在α≥0.001时，将所述新样本点添加到所述输入数据集中，以对每一个观测点的初始克里金代理模型重新训练，持续迭代，直到该观测点的克里金代理模型满足序贯抽样算法的收敛准则，得到最终该观测点的克里金代理模型；
其中，α定义为第d - 49次迭代到第d次迭代的确定性系数R 2的变化程度，d为自适应序贯抽样算法的迭代次数且d ≥ 50。
进一步的，在得到最终该观测点的克里金代理模型之后，所述方法还包括：
利用拉丁超立方体抽样生成多个验证样本点；
针对每一个观测点，分别使用为该观测点构造的动态数值模型和最终该观测点的克里金代理模型计算每个验证样本点对应的最终堆积深度；
将该观测点的动态数值模型计算得到的每个验证样本点对应的最终堆积深度与最终该观测点的克里金代理模型计算得到的每个验证样本点对应的最终堆积深度进行比较，利用所述R 2评价最终该观测点的克里金代理模型的预测质量；其中，所述R 2用α来衡量最终该观测点的克里金代理模型的预测准确度随迭代次数的增加而改善的程度。
进一步的，所述步骤S3包括：
子步骤S3-1：将所述潜在滑坡预估的滑坡堆积区网格化，选取所述潜在滑坡整个运动过程中各网格位置的最大运动高度和最大运动速度作为评价指标，建立所述潜在滑坡在每一个网格的最大运动高度功能函数
Figure 844872DEST_PATH_IMAGE039
和最大运动速度功能函数
Figure 573794DEST_PATH_IMAGE040
：
Figure 839559DEST_PATH_IMAGE041
；
Figure 337536DEST_PATH_IMAGE042
；
Figure 109183DEST_PATH_IMAGE043
；
上式中，n和m分别是网格在所在位置的列数和行数，N和M分别是所述潜在滑坡预估的滑坡堆积区网格化后的列总数和行总数；
Figure 59822DEST_PATH_IMAGE044
为所述潜在滑坡在（n, m）网格处最大运动高度克里金代理模型输出的最大运动高度，
Figure 942327DEST_PATH_IMAGE045
为所述潜在滑坡在（n, m）网格处最大运动速度克里金代理模型输出的最大运动速度；
Figure 294811DEST_PATH_IMAGE046
和
Figure 237359DEST_PATH_IMAGE047
分别是最大运动高度的阈值和最大运动速度的阈值；
子步骤S3-2：基于所述潜在滑坡在每一个网格的最大运动高度功能函数
Figure 612977DEST_PATH_IMAGE048
，计算得到所述潜在滑坡的最大运动高度克里金代理模型矩阵
Figure 33594DEST_PATH_IMAGE049
：
Figure 506164DEST_PATH_IMAGE050
；
子步骤S3-3：基于所述潜在滑坡在每一个网格的最大运动速度功能函数
Figure 619613DEST_PATH_IMAGE051
，计算得到所述潜在滑坡的最大运动速度克里金代理模型矩阵
Figure 544844DEST_PATH_IMAGE052
：
Figure 769152DEST_PATH_IMAGE053
。
进一步的，在通过所述子步骤S3-2计算得到所述潜在滑坡的最大运动高度克里金代理模型矩阵
Figure 17599DEST_PATH_IMAGE054
以及通过所述子步骤S3-3计算得到所述潜在滑坡的最大运动速度克里金代理模型矩阵
Figure 301950DEST_PATH_IMAGE055
之后，所述步骤S4包括：
将所述随机变量的后验分布输入到所述潜在滑坡的最大运动高度克里金代理模型矩阵
Figure 714477DEST_PATH_IMAGE056
，输出得到每一个网格的最大运动高度；
每个网格位置（n，m）处的最大运动高度超越概率计算为
Figure 742476DEST_PATH_IMAGE057
：
Figure 924058DEST_PATH_IMAGE058
；
上式中，T为后验分布中的后验样本数，θ t 是由马尔可夫链蒙特卡罗模拟采样获得的第t个后验样本；当
Figure 113731DEST_PATH_IMAGE059
时，表示
Figure 13554DEST_PATH_IMAGE060
发生；
所述潜在滑坡的最大运动高度超越概率矩阵为：
Figure 517348DEST_PATH_IMAGE061
；
以及，
将所述随机变量的后验分布输入到所述潜在滑坡的最大运动速度克里金代理模型矩阵
Figure 819016DEST_PATH_IMAGE062
，输出得到每一个网格的最大运动速度；
每个网格位置（n，m）处的最大运动速度超越概率计算为
Figure 179590DEST_PATH_IMAGE063
：
Figure 566709DEST_PATH_IMAGE064
；
上式中，T为后验分布中的后验样本数，θ t 是由马尔可夫链蒙特卡罗模拟采样获得的第t个后验样本；当
Figure 936511DEST_PATH_IMAGE065
时，表示
Figure 827106DEST_PATH_IMAGE066
发生；
所述潜在滑坡的最大运动速度超越概率矩阵为：
Figure 624161DEST_PATH_IMAGE067
。
进一步的，所述步骤S5包括：
利用三维地理信息系统平台将得到的所述潜在滑坡的最大运动高度超越概率矩阵以及所述潜在滑坡的最大运动速度超越概率矩阵投影到栅格文件上；其中，所述栅格文件用于在所述三维地理信息系统平台上可视化；
将所述栅格文件叠加在所述潜在滑坡的数字高程模型或卫星图像地图上，得到用于单体滑坡定量风险评价的概率运动危险性图。
进一步的，所述多个观测点的选取来源包括：
临近所述过去滑坡的滑坡堆积区底边界的点；
临近所述过去滑坡的滑坡堆积区最顶部边界的点；
临近所述过去滑坡的滑坡堆积区左右边界的点。
本发明实施例包括以下优点：
本发明基于现有提出的二维运动预测基本框架，提出了一种改进的三维滑坡运动危险性概率评价的总体框架，在该框架下，所提出的三维滑坡运动危险性概率评价方法在类似于潜在滑坡的过去滑坡的滑坡堆积区内选择多个观测点，对多个观测点中的每一个观测点构造动态数值模型，动态数值模型用于滑坡运动分析，并以动态数值模型输出的响应值即最终堆积深度作为观测信息，利用基于马尔可夫链蒙特卡罗模拟的多观测贝叶斯反分析方法对为每一个观测点构造的动态数值模型的输入参数（即随机变量）的不确定性进行校准，得到了后验分布。然后，使用后验分布作为潜在滑坡的输入，以在三维地形上计算潜在滑坡预估的滑坡堆积区中每一个网格的最大运动高度超越概率以及最大运动速度超越概率，并生成用于单体滑坡定量风险评价的概率运动危险性图。为了在不损失精度的前提下提高计算效率，本发明还提出了一种多响应克里金代理模型，该模型采用了序贯自适应采样方法，以全局近似模型输入和输出之间的关系，将其完全替代原有的动态数值模型用于贝叶斯反分析和概率运动危险性评估，在不损失精度的前提下大大提高了计算效率。
附图说明
图1是本发明三维滑坡运动危险性概率评价的总体框架图；
图2是本发明一种三维滑坡运动危险性概率评价方法的步骤流程图；
图3是CJ8号滑坡位置示意图；
图4是CJ8号滑坡发生两次滑动的示意图；
图5(a)和图5(b)分别是CJ8号滑坡第一次滑动前后的影像，图5(c)和图5(d)分别是CJ8号滑坡第二次滑动前后的影像。
图6是本实例选择用于贝叶斯反分析的观测点网格位置和最终堆积深度示意图；
图7是本实例拟合的模型偏差系数的概率密度分布直方图；
图8是本实例7个观测点处最终堆积深度代理模型的R 2随迭代次数增加的变化示意图；
图9(a)是有效残余摩擦角
Figure 419948DEST_PATH_IMAGE068
的先验分布和后验分布示意图，图9(b)是有孔隙压力比r u 的先验分布和后验分布示意图，图9(c)是后验样本点的示意图；
图10是通过后验均值模拟的CJ8号滑坡第一次滑动的堆积区示意图；
图11是本实例对17017个最大运动高度克里金代理模型和17017个最大运动速度克里金代理模型的确定性系数R 2的统计图；
图12是通过后验均值模拟的CJ8号滑坡第二次滑动堆积区示意图；
图13(a)~图13 (d)分别是图12中A点的最大运动高度和最大运动速度的直方图和超越概率曲线；
图14(a)~图14(d) 分别是图12中B点的最大运动高度和最大运动速度的直方图和超越概率曲线；
图15是基于最大运动高度超越概率的滑坡运动危险性区划图(
Figure 593440DEST_PATH_IMAGE069
=0.0m)；
图16(a)~图16(d)是
Figure 604121DEST_PATH_IMAGE070
=0.75,
Figure 572077DEST_PATH_IMAGE071
=1.5m,
Figure 668209DEST_PATH_IMAGE072
=3m/s以及
Figure 379813DEST_PATH_IMAGE073
=6m/s条件下的滑坡运动危险性概率评价图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
针对本发明的技术问题，参照图1和图2，本发明公开了一种三维滑坡运动危险性概率评价方法，该方法可以包括以下步骤：
步骤S1：确定潜在滑坡对应的过去滑坡，并在所述过去滑坡的滑坡堆积区内选择多个观测点；
在本发明中，由于滑坡物理性质和滑坡类型广泛地影响着滑坡的运动特性，只有当潜在事件在以前有记录良好的类似事件时，通过对过去事件的反分析校准的输入参数才能直接应用到潜在事件的运动预测中。因此，本发明根据潜在滑坡的物质性质和滑坡类型，确定同样物质性质和滑坡类型的过去滑坡。其中，物质性质可以为地质环境、岩土性质等特征，滑坡类型可以为长距离滑坡特征、短距离滑坡特征等。
在本发明中，可以通过机载激光探测与测距(LiDAR)、地面激光扫描(TLS)和无人机(UAV)等对过去滑坡在滑动前后地形进行测绘，如此可以获得过去滑坡在滑动前后的高分辨率数字高程模型，进而得到过去滑坡的滑坡堆积区。
在本发明中，通过将过去滑坡的滑坡堆积区网格化，可以得到大量的网格，这些网格分布在不同的网格位置。本发明可以从网格化后的滑坡堆积区中选择不同网格位置的网格作为观测点(也可称为观察点)。观测点的选取原则是：1）观测点的选取足以控制整个过去滑坡的滑坡堆积区；2）根据不同实际案例的滑坡堆积区实际形状选择观测点，主要原则是在滑坡堆积区边界附近选点；3）在满足观测点的选取足以控制整个过去滑坡的滑坡堆积区的前提下，观测点不宜过多。
在本发明一实施例中，可以按照以下选取来源进行多个观测点的选取：①临近过去滑坡的滑坡堆积区底边界的点；临近过去滑坡的滑坡堆积区最顶部边界的点；临近过去滑坡的滑坡堆积区左右边界的点。
步骤S2：对所述多个观测点中的每一个观测点构造克里金代理模型，并基于每一个观测点的克里金代理模型、贝叶斯定理以及差分进化适应大都会市(differentialevolution adaptive metropolis，简称DREAM)算法，对预先确定的随机变量的先验信息进行反分析，得到所述随机变量的后验分布；
在本发明中，采用基于贝叶斯定理进行反分析时，观测信息可以采用运动距离、不同位置的最终堆积深度、运动速度或冲击力等。
在本发明一实施例中，提供了至少以最终堆积深度为观测信息进行反分析，假设所选择的观测点的个数为K个；步骤S2可以通过以下步骤实现：
子步骤S2-1：对所述K个观测点中的每一个观测点构造动态数值模型，并建立每一个观测点的克里金代理模型，将所述动态数值模型的输入参数作为随机变量θ，利用同一个观测点的克里金代理模型计算该观测点处所述随机变量对应的最终堆积深度模拟值，获得所述K个观测点的最终堆积深度模拟值
Figure 245001DEST_PATH_IMAGE074
，其中，
Figure 321542DEST_PATH_IMAGE075
表示第K个观测点的最终堆积深度模拟值；
在本发明一实施例中，由于基底阻力
Figure 904970DEST_PATH_IMAGE076
在滑坡动态运动模拟中起着关键作用，因此子步骤S2-1中的对所述K个观测点中的每一个观测点构造动态数值模型，并建立每一个观测点的克里金代理模型，具体可以采用以下方法实现：
子步骤S2-101：使用摩擦模型计算基底阻力
Figure 154685DEST_PATH_IMAGE077
：
Figure 139959DEST_PATH_IMAGE078
；
上式中，𝜎 = ρgh，𝜎表示基底正应力，ρ表示滑坡物质的密度，h表示滑坡物质的运动高度，g代表重力；
Figure 449718DEST_PATH_IMAGE079
是有效残余摩擦角；r u 是孔隙压力比系数，反映了基底材料的液化程度，取值范围为0~1；
子步骤S2-102：基于摩擦模型，利用质量流模型（Massflow）对所述过去滑坡进行运动模拟，对K个观测点中的每一个观测点构造动态数值模型；其中，质量流模型是基于三维纳维-斯托克斯（3D Navier–Stokes）方程沿深度积分的连续介质模型；滑坡的运动行为受深度积分质量守恒方程和动量守恒方程的控制，该模型可以采用基于Mac-Cormack有限差分格式的数值方法求解方程。如此，本发明完成了对滑坡三维运动模拟的动态数值模型的构造。
考虑到概率反分析和概率运动预测都涉及数万次重复运动模拟。因此，对于复杂且计算成本高的动态数值模型在计算成本上是不可行的，如Massflow。因此，本发明所提出的多响应克里金代理模型被用来近似输入随机变量
Figure 520442DEST_PATH_IMAGE080
和不同网格位置的输出强度之间的关系，基于拟合的克里金代理模型，可以有效地进行概率反分析和概率运动预测。代理模型的目标是用最少的训练样本建立最精确的代理模型；因此，需要有效的样本设计策略。混合自适应采样策略是网格填充和自适应方法的结合，其针对已有的单个响应系统开发。在此，本发明将混合自适应采样策略修改为适合多响应克里金代理模型。具体而言：
子步骤S2-103：将所述动态数值模型的输入参数
Figure 760799DEST_PATH_IMAGE081
和r u 定义为所述随机变量
Figure 600579DEST_PATH_IMAGE082
，并通过均匀设计（Uniform Design，简称UD）对所述随机变量进行采样，得到输入数据集，所述输入数据集包括多个初始采样点。在本发明中，将动态数值模型的输入参数
Figure 81239DEST_PATH_IMAGE083
和r u 定义为随机变量，而ρ定义为确定性变量，因为它对运动模拟结果的影响有限。初始采样的主要目的是在随机变量的整个参数网格中生成采样点，以捕捉上述所构造的动态数值模型的全局行为。
子步骤S2-104：将所述输入数据集的每个初始采样点代入每一个观测点的动态数值模型中，以分别输出每一个观测点位置对应的响应值，基于所述输入数据集以及每一个观测点位置对应的响应值，对克里金模型(Kriging)进行训练，得到每一个观测点的初始克里金代理模型。在本发明中，每个初始采样点的响应值为最终堆积深度的模拟值（即最终堆积深度模拟值），克里金模型可以采用Lophaven等人(2002)开发的DACE（Design andAnalysis of Computer Experiments)工具箱构建。本发明通过将输入数据集中的每个初始采样点和该初始采样点对应的响应值作为训练数据，可以实现对克里金模型的训练，从而得到该观测点的初始克里金代理模型。
子步骤S2-105：从候选样本池中选择一个新样本点，其中，所述候选样本池
Figure 373680DEST_PATH_IMAGE084
中的候选样本点通过拉丁超立方体抽样(Latin hypercube sampling，缩写LHS)生成；j=1，2，⋯,s；其中s是候选样本点的数量；
在本发明实施例中，考虑本发明的动态数值模型的多响应性质，即同时考虑所有近似代理模型的行为，子步骤S2-105中的从候选样本池中选择一个新样本点，可以通过以下步骤实现：
通过最大化质量参数
Figure 230778DEST_PATH_IMAGE085
从候选样本池中选择一个新样本点，其中，
Figure 925064DEST_PATH_IMAGE086
的表达式如下：
Figure 576625DEST_PATH_IMAGE087
；
上式中，
Figure 559625DEST_PATH_IMAGE088
是第j个候选样本点与当前设计样本点之间的最小欧氏距离，其中，当前设计样本点包括所述初始采样点和所述新样本点；
Figure 954834DEST_PATH_IMAGE089
是第j个候选样本点处的预测方差，由每一个观测点的第i个克里金代理模型进行估计获得；
w i 为权重因子，表示每一个观测点的第i个克里金代理模型对总预测方差的贡献，其中：
Figure 503627DEST_PATH_IMAGE036
；
上式中，
Figure 591669DEST_PATH_IMAGE090
是每一个观测点的第i个克里金代理模型的预测方差的平均值，
Figure 858702DEST_PATH_IMAGE091
是每一个观测点的第i个克里金代理模型的最大预测方差，p是每一个观测点的克里金代理模型的数目。
子步骤S2-106：将所述新样本点添加到所述输入数据集中，以对每一个观测点的初始克里金代理模型重新训练，持续迭代，直到该观测点的克里金代理模型满足预先设定的收敛准则，得到最终该观测点的克里金代理模型。
可选的，子步骤S2-106的通过以下条件实现：
在α≥0.001时，将所述新样本点添加到所述输入数据集中，以对每一个观测点的初始克里金代理模型重新训练，持续迭代，直到该观测点的克里金代理模型满足序贯抽样算法的收敛准则，得到最终该观测点的克里金代理模型；
其中，α定义为第d - 49次迭代到第d次迭代的确定性系数R 2的变化程度，d为自适应序贯抽样算法的迭代次数且d ≥ 50。
综上步骤，本发明可以获得所选取的多个观测点的克里金代理模型，多个观测点的克里金代理模型能够输出不同位置处的最终堆积深度，如此得到上述所提及的基于多响应的克里金代理模型。
在得到最终该观测点的克里金代理模型之后，为保证本发明实施例的克里金代理模型替代动态数值模型进行使用的可靠性，具体使用前，在本发明另一实施例中，还提出了对该克里金代理模型的近似精度进行检验的方法，可以包括以下步骤：
利用拉丁超立方体抽样生成多个验证样本点；
针对每一个观测点，分别使用为该观测点构造的动态数值模型和最终该观测点的克里金代理模型计算每个验证样本点对应的最终堆积深度；
将该观测点的动态数值模型计算得到的每个验证样本点对应的最终堆积深度与最终该观测点的克里金代理模型计算得到的每个验证样本点对应的最终堆积深度进行比较，利用所述R 2评价最终该观测点的克里金代理模型的预测质量；其中，所述R 2用α来衡量最终该观测点的克里金代理模型的预测准确度随迭代次数的增加而改善的程度。
子步骤S2-2：利用设定的模型偏差系数
Figure 57602DEST_PATH_IMAGE092
和所述K个观测点的最终堆积深度模拟值
Figure 460902DEST_PATH_IMAGE093
，计算得到所述K个观测点的最终堆积深度观测值
Figure 906796DEST_PATH_IMAGE094
，其中，
Figure 661125DEST_PATH_IMAGE095
；
虽然通过上述方法本发明计算得出了K个观测点的最终堆积深度模拟值，但由于模型误差和测量误差，因此模拟的最终堆积深度模拟值与观测得到的最终堆积深度观测值之间总是存在偏差，偏差可以表示为模型偏差系数
Figure 398137DEST_PATH_IMAGE096
。因此，本发明结合模型偏差系数
Figure 921522DEST_PATH_IMAGE097
与K个观测点的最终堆积深度模拟值
Figure 351366DEST_PATH_IMAGE098
，可以得出K个观测点的最终堆积深度观测值
Figure 327413DEST_PATH_IMAGE099
。实际中，模型偏差系数
Figure 868115DEST_PATH_IMAGE100
可以通过过去滑坡数据计算获得，即可以通过其他滑坡的数字高程模型（Digital Elevation Model，简称DEM）滑动前后观察到的最终堆积深度观测值与模型输出最佳模拟结果（最终堆积深度模拟值）进行比较计算获得。通过计算模型偏差系数，可以确保其适用于滑坡运动的概率反分析。
子步骤S2-3：基于所述过去滑坡的K个观测点对应的最终堆积深度观测值，通过贝叶斯定理对所述随机变量的先验分布
Figure 449269DEST_PATH_IMAGE101
进行反分析，利用DREAM算法对反分析结果进行采样，得到所述随机变量的后验分布
Figure 784436DEST_PATH_IMAGE102
；
Figure 513357DEST_PATH_IMAGE103
；
上式中，
Figure 592172DEST_PATH_IMAGE104
是似然函数，用于表征现场测量数据与动态数值模型模拟结果之间的匹配程度，表示为所述随机变量θ的联合条件概率密度函数。需要说明的是，先验分布
Figure 90149DEST_PATH_IMAGE105
也称为随机变量θ的先验概率密度函数，后验分布
Figure 596217DEST_PATH_IMAGE106
也称为随机变量θ的后验概率密度函数；其中，在假设不同网格位置之间的模型偏差系数
Figure 812435DEST_PATH_IMAGE107
是不相关的情况下，
Figure 881891DEST_PATH_IMAGE108
可以表示为：
Figure 234375DEST_PATH_IMAGE021
；
上式中，
Figure 176923DEST_PATH_IMAGE109
表示
Figure 614857DEST_PATH_IMAGE023
的概率密度函数在
Figure 35475DEST_PATH_IMAGE110
对应的概率密度。
马尔可夫链蒙特卡罗(Markov Chain Monte Carlo，简称MCMC)模拟的基本思想是从随机变量的先验分布中采样，然后对这些样本进行修正，以更好地逼近并最终收敛到目标后验分布。DREAM算法是一种自适应算法，对复杂的高维后验分布具有很高的采样效率，在采样过程中能够自适应地更新建议分布的尺度和方向。因此，本发明采用DREAM算法进行多链MCMC仿真，然后对随机变量的先验分布反分析获得的后验分布(反分析结果)进行采样，可以得到随机变量θ 的后验分布。
步骤S3：将所述潜在滑坡预估的滑坡堆积区网格化，选取所述潜在滑坡整个运动过程中各网格位置的最大运动高度和最大运动速度作为评价指标，建立所述潜在滑坡在每一个网格的最大运动高度功能函数和最大运动速度功能函数，并计算得到所述潜在滑坡的最大运动高度克里金代理模型矩阵和最大运动速度克里金代理模型矩阵；
在本发明一实施例中，步骤S3可以通过以下步骤实现：
子步骤S3-1：将所述潜在滑坡预估的滑坡堆积区网格化，选取所述潜在滑坡整个运动过程中各网格位置的最大运动高度和最大运动速度作为评价指标，建立所述潜在滑坡在每一个网格的最大运动高度功能函数
Figure 508044DEST_PATH_IMAGE111
和最大运动速度功能函数
Figure 293598DEST_PATH_IMAGE112
：
Figure 218828DEST_PATH_IMAGE113
；
Figure 443136DEST_PATH_IMAGE042
；
Figure 770212DEST_PATH_IMAGE114
；
上式中，n和m分别是网格在所在位置的列数和行数，N和M分别是所述潜在滑坡预估的滑坡堆积区网格化后的列总数和行总数；
Figure 54563DEST_PATH_IMAGE115
为所述潜在滑坡在（n, m）网格处最大运动高度克里金代理模型输出的最大运动高度，
Figure 467090DEST_PATH_IMAGE116
为所述潜在滑坡在（n,m）网格处最大运动速度克里金代理模型输出的最大运动速度；
Figure 229510DEST_PATH_IMAGE117
和
Figure 621480DEST_PATH_IMAGE118
分别是最大运动高度的阈值和最大运动速度的阈值；
子步骤S3-2：基于所述潜在滑坡在每一个网格的最大运动高度功能函数
Figure 76733DEST_PATH_IMAGE119
，计算得到所述潜在滑坡的最大运动高度克里金代理模型矩阵
Figure 976555DEST_PATH_IMAGE120
：
Figure 542666DEST_PATH_IMAGE050
；
子步骤S3-3：基于所述潜在滑坡在每一个网格的最大运动速度功能函数
Figure 578755DEST_PATH_IMAGE051
，计算得到所述潜在滑坡的最大运动速度克里金代理模型矩阵
Figure 204908DEST_PATH_IMAGE121
：
Figure 592027DEST_PATH_IMAGE053
。
在本发明中，潜在滑坡预估的（即可能的）滑坡堆积区可参照过去滑坡的滑坡堆积区进行划分，但一般比过去滑坡的滑坡堆积区面积大一些。将潜在滑坡的滑坡堆积区网格化的方法与过去滑坡类似，在此不多赘述。由于整个运动过程中的最大运动强度可以反映滑坡的破坏力，因此，本发明选取潜在滑坡整个运动过程中各网格位置的最大运动高度和最大运动速度作为评价指标。其中，
Figure 633933DEST_PATH_IMAGE122
和
Figure 790108DEST_PATH_IMAGE123
的具体值可以根据相应承灾体的承受值进行确定。
由于一个网格具有基于最大运动高度克里金代理模型以及基于最大运动速度克里金代理模型，因此，本发明的潜在滑坡一共需建立2 × N × M个克里金代理模型，这些克里金代理模型共享样本点，因此，实施时只需要一个自适应采样过程。每个网格位置处的最大运动高度克里金代理模型和最大运动速度克里金代理模型可参照过去滑坡中每个观测点的克里金代理模型的构造方法，在此不多赘述。
步骤S4：将所述随机变量的后验分布输入到所述潜在滑坡的最大运动高度克里金代理模型矩阵中，输出得到每一个网格的最大运动高度，以及将所述随机变量的后验分布输入到所述最大运动速度克里金代理模型矩阵中，输出得到每一个网格的最大运动速度，并对每一个网格的最大运动高度超越概率以及最大运动速度超越概率分别进行计算，获得所述潜在滑坡的最大运动高度超越概率矩阵以及所述潜在滑坡的最大运动速度超越概率矩阵；
在本发明一实施例中，在通过所述子步骤S3-2计算得到所述潜在滑坡的最大运动高度克里金代理模型矩阵
Figure 587162DEST_PATH_IMAGE124
以及通过所述子步骤S3-3计算得到所述潜在滑坡的最大运动速度克里金代理模型矩阵
Figure 195998DEST_PATH_IMAGE125
之后，所述步骤S4包括：
将利用DREAM算法获得的随机变量的后验分布输入到潜在滑坡的最大运动高度克里金代理模型矩阵
Figure 369491DEST_PATH_IMAGE126
，输出得到每一个网格的最大运动高度；
当后验分布中的后验样本数T足够大时（比如后验样本数的数量达到10000以上），每个网格位置（n，m）处的最大运动高度超越概率计算为
Figure 380172DEST_PATH_IMAGE127
：
Figure 535079DEST_PATH_IMAGE058
；
上式中，T为后验分布中的后验样本数，θ t 是由MCMC模拟采样获得的第t个后验样本；当
Figure 631211DEST_PATH_IMAGE128
时，表示
Figure 77235DEST_PATH_IMAGE129
发生；
所述潜在滑坡的最大运动高度超越概率矩阵为：
Figure 208003DEST_PATH_IMAGE061
。
在本发明中，每个网格位置（n，m）处的最大运动速度超越概率的计算方法与最大运动高度超越概率的计算方法相同，仅需把潜在滑坡的最大运动高度超越概率矩阵中“h”替换为“v”即可，即：
将所述随机变量的后验分布输入到所述潜在滑坡的最大运动速度克里金代理模型矩阵
Figure 346860DEST_PATH_IMAGE130
，输出得到每一个网格的最大运动速度；
当后验分布中的后验样本数T足够大时（比如后验样本数的数量达到10000以上），每个网格位置（n，m）处的最大运动速度超越概率计算为
Figure 930288DEST_PATH_IMAGE131
：
Figure 180004DEST_PATH_IMAGE064
；
上式中，T为后验分布中的后验样本数，θ t 是由MCMC模拟采样获得的第t个后验样本；当
Figure 102960DEST_PATH_IMAGE132
时，表示
Figure 412719DEST_PATH_IMAGE133
发生；
所述潜在滑坡的最大运动速度超越概率矩阵为：
Figure 483443DEST_PATH_IMAGE067
。。
步骤S5：在三维地理信息系统平台(Geographic Information System或 Geo－Information system，GIS)上，根据所述潜在滑坡的最大运动高度超越概率矩阵以及所述潜在滑坡的最大运动速度超越概率矩阵，生成用于单体滑坡定量风险评价的概率运动危险性图。
在本发明中，通过使用不同的临界阈值(
Figure 536850DEST_PATH_IMAGE134
和
Figure 376630DEST_PATH_IMAGE135
)，可获得最大堆高和最大速度的三维超越概率分布图(即运动危险性图)，可直接用于滑坡易损性评估。
本发明生成了用于单体滑坡定量风险评价的概率运动危险性图，概率运动危险性图明确显示了滑坡路径上每个网格位置的运动强度超过一定水平的概率的网格分布，反映了某位置受到威胁的程度。因此，危险性区划有助于确定更脆弱的地区，利用这些灾危险性区划图为土地规划和管理提供可能性。此外，由于潜在滑坡运动危险性概率评价图中使用的强度定义为潜在破坏力的度量，因此可以直接在易损性曲线中使用，以进行概率风险分析。在防灾结构的设计中，运动强度超越概率也可以为基于可靠性的防灾结构设计提供信息。
在本发明一实施例中，步骤S5可以通过以下步骤实现：利用三维地理信息系统平台将得到的所述潜在滑坡的最大运动高度超越概率矩阵以及所述潜在滑坡的最大运动速度超越概率矩阵投影到栅格文件上；其中，所述栅格文件用于在三维地理信息系统平台上可视化；
将所述栅格文件叠加在所述潜在滑坡的数字高程模型或卫星图像地图上，得到用于单体滑坡定量风险评价的概率运动危险性图。在本发明中，关于栅格文件如何可以在三维地理信息系统平台上可视化的技术以及叠加在DEM或卫星图像地图上的技术，可参见相关技术，在此不多赘述。将栅格文件叠加在所述潜在滑坡的DEM或卫星图像地图上，用于单体滑坡定量风险评价的概率运动危险性图的以DEM图或卫星图像地图的形式显示，如此，可以便于决策者从该图中了解可能受到滑坡影响的区域，简单直观。
为验证本发明实施例的有效性，下面采用一研究实例详细说明。
研究区位于黑方台，距中国甘肃省兰州市45公里（见图3），阶地被厚度为26–48 m的马兰黄土覆盖。由于灌溉始于20世纪60年代，阶地边缘经常发生地面沉降和滑坡。黑方台阶地滑坡一般可分为两类：浅层黄土滑坡和黄土基岩滑坡。由于黄土层基底的饱和导致黄土的应变软化，因此浅层黄土滑坡的运动距离相对较长。图3中的CJ8号滑坡位于黑方台阶地东侧。参见图4，2015年共发生过两次滑动，利用高分辨率DEM在Massflow中获得了两次滑动的前后影像，详见图5(a)~图5(d)。从图5(b)可知，第一次滑动发生在2015年3月29日，总体积为2.14× 104 m3，运动距离375 m。第一次滑坡的流动指数H/L为0.160，具有明显的长距离滑坡特征。从图5(d)可知，第二次滑动发生在2015年9月20日，总体积为1.22× 104 m3，运动距离320 m。第二次滑坡的流动指数H/L为0.156，也具有长距离滑坡的特征。
由于两个滑坡在物质性质和滑坡类型上有明显的相似性。因此，为了说明概率运动危害评估方法，在本发明中，将第一次滑动视为过去事件（参见图5(a)和图5(b)的影像），用于多观测贝叶斯反分析，将第二次滑动滑假定为潜在的未来事件（参见图5(c)和图5(d)的影像），通过对过去滑坡（过去事件）进行反分析得到的校准输入参数可用于潜在滑坡（未来事件）的正向运动预测，以评估其滑坡运动危险性。
一、CJ8号滑坡第一次滑动（过去滑坡）的反分析
通过差分CJ8#滑坡在2015年两次滑动的滑前和滑后DEM，可以确定两次滑动的原始地形高程、初始滑动位置、滑体厚度和滑坡堆积区（简称堆积区）范围。在此基础上，以地形高程和初始体积（近似为初始滑动面积与滑体厚度的乘积）为基础，建立了CJ8号滑坡体在质量流中第一次滑动的数值模型。在DEM中，网格尺寸设定为3m，网格数为90× 143，运动模拟总时间为100s，时间步长设置为可变时间步长（初始时间步长为0.05s，courant数为0.25）。在实例中，采用摩擦模型
Figure 857290DEST_PATH_IMAGE136
进行计算，其中有效残余摩擦角
Figure 149731DEST_PATH_IMAGE137
和孔隙压力比r u 被定义为随机变量。如表1所示，示出了CJ8号滑坡的有效残余摩擦角
Figure 193779DEST_PATH_IMAGE138
和孔隙压力比r u 这两个输入参数的先验分布，并假设
Figure 888066DEST_PATH_IMAGE139
和r u 在统计上是独立的。
表1：输入参数的先验分布
输入参数 有效残余摩擦角
Figure 539627DEST_PATH_IMAGE140
(°)		 孔隙压力比<i>r</i><sub><i>u</i></sub>
方法 34.60 0.50
标准差 3.45 0.10
分配类型 正态分布 正态分布
二、观测集和模型偏差系数
选择四类最终堆积深度观测值来控制滑坡堆积区的主要几何特征（见图6）：1）临近滑坡堆积区底边界的点（例如，控制滑坡运动距离的点1）；2）临近滑坡堆积区最顶部边界的点（例如，观测点6）；3）临近滑坡堆积区左右边界的点（例如，控制滑坡横向扩展的观测点2、3、4、5和7）；4）根据不同案例的实际堆积区形状选择其他观测点。在本实施例中，最后选取7个观测点进行贝叶斯反分析，对输入参数进行校正，可以控制滑坡堆积区的几何形状和最终堆积深度的网格分布。从滑动前后DEM之间的差异获得每个观测点的最终堆积深度。值得注意的是，尽管根据滑坡前后的DEM可以获得大量的最终堆积深度观测值，但本发明仅使用有限的观测值（即7个观测点）即可进行贝叶斯校正。主要原因是：(i)滑坡堆积区几何观测点的选取足以控制整个堆积区；(ii)选择更多的观测点意味着需要建立更多的克里金代理模型，这增加了计算成本。
通过将观察到的最终堆积深度与最佳模拟结果进行比较，本发明计算了两者的差异，并计算了模型偏差系数，如图7所示。对数正态分布被发现是一个可接受的适合模型偏差系数的分布。分布的平均值和标准差分别为1.027和0.797。因此，假设每个观测点的模型偏差系数服从对数正态分布，平均值
Figure 584943DEST_PATH_IMAGE141
为1.027，标准差
Figure 980152DEST_PATH_IMAGE142
为0.797，并且假设不同观测点的模型偏差系数是独立的。然后，利用模型偏差系数的概率分布建立
Figure 528945DEST_PATH_IMAGE143
中的似然函数。
为了降低贝叶斯反分析的计算成本，构造了计算简单的过去滑坡的观测点的克里金代理模型来模拟随参数值变化的模型输出（即每个观测点的最终堆积深度模拟值）。首先，使用均匀设计生成20个初始采样点，初始采样的主要目的是在整个参数网格中生成采样点，以捕捉动态数值模型的全局行为。此后，通过拉丁超立方体抽样生成10000个候选样本点作为新样本的候选样本池。因为5σ的区间概率（大约0.99999943）包含了大部分信息，θ= [
Figure 616987DEST_PATH_IMAGE144
, r u ] 的取样下限为[17.35,0]，上限为[51.85,1]。然后，利用动态数值模型计算每个观测点对应于每个初始采样点的模型响应（即最终堆积深度模拟值），利用初始采样点和该初始采样点对应的最终堆积深度模拟值作为训练数据，可以对克里金模型进行训练，从而构建了该观测点的初始克里金代理模型。此外，通过求解
Figure 821703DEST_PATH_IMAGE145
的表达式中的优化问题，从候选样本池中选择一个新的样本点。新的样本点被添加到初始样本集中，以重新训练该观测点的初始克里金代理模型。迭代继续进行，直到满足预先设定的收敛准则。
为了验证所构建的克里金代理模型的准确性，利用拉丁超立方体抽样生成了100个额外的测试样品，并分别使用Massflow和最终克里金代理模型计算了相应的最终堆积深度。然后，利用确定性系数(R 2)来评价多响应克里格模型的预测质量。在本实例中，可使用R 2在迭代过程中的变化程度α（α定义为第d-49次迭代到第d次迭代的R 2的变异系数，d为自适应序贯抽样算法的迭代次数且d≥50）来衡量随着迭代次数增加代理模型精度的改善程度，当α的值很小时，即说明随着样本的继续增加而代理模型精度的改善程度并不大，此次迭代即可终止。因此，自适应序贯采样算法的收敛准则可定义为α低于某一目标值，在本实例中该目标值设为0.001。同时，为了避免收敛准则过于严苛，设置最大500次迭代作为附加的收敛准则。当d=288时，收敛准则α≤0.001满足。图8显示了随着迭代次数的增加，每个克里金代理模型中的R 2的变化。研究发现，随着自适应序贯抽样样本数的增加，克里金代理模型的预测质量逐渐提高并趋于稳定。因此，验证了本发明根据R 2的相对变化来确定收敛准则是可行的。各克里金代理模型的最终确定性系数R 2分别为0.9988、0.9914、0.9806、0.9937、0.9947、0.9988和0.9899，表明克里金代理模型具有足够的预测精度。因此，本发明的克里金代理模型可以完全替代动态数值模型来进行贝叶斯反分析。
三、贝叶斯反分析结果
采用DREAM算法进行多链MCMC仿真，对随机变量的后验分布进行采样。在这项研究中，马尔可夫链的数量是四个（通常定义为随机变量数量的两倍），生成的样本数是60000个。这里，Rubin收敛诊断R stat 准则（R stat ≤1.2）验证了算法的收敛性。通过测试MCMC模拟，发现R stat 值均小于1.2，表明60000的链长度足以确保链内收敛以实现稳定的后验分布。此外，在本发明中，前10000个样本被丢弃，最后50000个样本被用作后验分布的平稳样本，如此可以克服马尔可夫链在开始阶段不可能达到平衡状态的问题。
图9(a)和图9(b)分别显示
Figure 755024DEST_PATH_IMAGE146
和 r u 的后直方图。后验统计量
Figure 423903DEST_PATH_IMAGE147
和 r u 也进行了评估，平均值分别为31.324和0.769，标准差分别为3.532和0.039。输入参数的后验均值与先验均值相比有明显的调整，其中，r u 的标准差显著降低，说明贝叶斯反分析后该参数的不确定度降低。图9(c)显示了后验分布的样本点，可以清楚地观察到有效残余摩擦角
Figure 682846DEST_PATH_IMAGE148
与孔压比r u ；之间存在明显的相关性，相关系数为0.879。后验平均值的模拟堆积区与观测堆积区非常接近（见图10），表明多观测贝叶斯反分析是有效的。
四、用于运动预测的潜在滑坡的多响应克里金代理模型的建立
将第一次滑坡（过去滑坡）的校准结果应用于CJ8号滑坡第二次滑动（潜在滑坡）的概率运动预测。以地形高程和初始体积为基础，建立了CJ8号滑坡第二次滑动的动态数值模型。在动态数值模型中，网格尺寸设定为3m；网格数量为17017（143×119); 运动模拟总时间为100s；将时间步长设为可变时间步长（初始时间步长为0.05s，courant为0.25）。对于每个模拟，Massflow提供了栅格输出，这些输出由网格单元组成。为了生成滑坡运动危险性概率评价图，需要使用从第一次滑动的贝叶斯模型校准获得的后验样本，获得每个网格单元的最大运动高度和最大运动速度的概率密度。由于后验样本量（50000个）过大，无法通过直接重复评估动态数值模型来获得概率密度。因此，为每个网格单元（简称网格）构造了计算节省的克里金代理模型。最后，建立了17017个最大运动高度克里金代理模型和17017个最大运动速度克里金代理模型。因为α≤0.001过于严格，直到达到最大500次迭代，迭代才停止。为了进一步检验所构建的克里金代理模型的准确性，本实例还对17017个最大运动高度克里金代理模型和17017个最大运动速度克里金代理模型的决定系数R 2进行了统计，如图11所示。结果发现，对于最大运动高度的克里金代理模型，96.38%模型的R 2大于0.95，98.29%模型的R 2大于0.85，99.04%模型的R 2大于0.70。对于最大运动速度的克里金代理模型，发现94.06%的模型R 2大于0.95，97.26%的模型R 2大于0.85，98.63%的模型R 2大于0.70，这表明本发明所构建的克里金代理模型的近似精度是较高的。
五、潜在滑坡失稳概率危险评估与制图
图12显示了使用后验平均值和Massflow模拟的CJ8号滑坡第二次滑动模拟的最终堆积区。结果表明，模拟的最终堆积范围与实测堆积范围基本一致，表明反演结果具有较好的预测效果。然后，采用MCMC仿真和基于克里金代理模型方法进行了概率运动预测。基于每个网格单元上的克里金代理模型，计算了MCMC每个后验样本对应的最大运动高度和最大运动速度。
在此基础上，进一步研究了各网格单元的最大运动高度和最大运动速度的分布，得到了各网格单元的最大运动高度和最大运动速度的超越概率曲线。为了说明结果，选择了两个网格单元（图12中的A点和B点）进行说明。点A(
Figure 171596DEST_PATH_IMAGE149
)的最大运动高度、点A(
Figure 174187DEST_PATH_IMAGE150
)的最大运动速度的直方图和超越概率曲线如图13(a)~图13 (d)所示；点B(
Figure 697573DEST_PATH_IMAGE151
)的最大运动高度和点B(
Figure 48788DEST_PATH_IMAGE152
)的最大运动速度的直方图和超越概率曲线如图14(a)~图14 (d)所示。正态分布被选为
Figure 290414DEST_PATH_IMAGE153
和
Figure 831117DEST_PATH_IMAGE154
的最佳拟合，而广义极值分布和极值分布分别被选为
Figure 209008DEST_PATH_IMAGE155
和
Figure 809754DEST_PATH_IMAGE156
的最佳拟合。此外，在95%置信水平下
Figure 273096DEST_PATH_IMAGE157
,
Figure 617490DEST_PATH_IMAGE158
,
Figure 53151DEST_PATH_IMAGE159
和
Figure 559218DEST_PATH_IMAGE160
的置信区间分别为[1.14 m，1.57 m]，[12.01 m/s，16.12 m/s]，[0.31 m，0.95 m]和[1.49 m/s，7.41 m/s]。置信区间改进了对给定位置处最大运动高度和最大运动速度的估计。在本实例中，通过互补累积分布函数(CCDF)获得了任何网格单元的最大运动高度和最大运动速度的超越概率曲线，如图13~图14所示。每个网格单元的最大运动高度和最大运动速度的超越概率曲线有助于基于概率的防灾减灾设计。
通过选择一个统一的给定阈值，可以计算出每个网格单元的最大运动高度和最大运动速度的超越概率。图15为允许生成最大运动高度和最大运动速度的超越概率的网格分布图，图15显示了
Figure 775436DEST_PATH_IMAGE161
=0.0m下的概率运动危险性图。CJ8号滑坡可能影响的区域被划分为五个危险等级：极高(P e ≥ 50.0%)，高(10.0% ≤ P e < 50.0%)，中等 (1.0% ≤ P e <10.0%)，低（0.1%≤ P e <1.0%和极低（P e <0.1%）。其中，概率大于0.5的区域（即极高危险区）是防灾减灾工作中最需要关注的区域，因为事件发生的概率大于不发生的概率。如图15所示结果，可验证出极高危险区与本发明观测到的最终堆积区吻合良好，这再次验证了基于反分析的概率运动预测的准确性。值得注意的是，阈值可以根据用户的要求而变化。图16(a)~图16(d)显示了
Figure 657941DEST_PATH_IMAGE162
=0.75,
Figure 10425DEST_PATH_IMAGE163
=1.5 m,
Figure 687394DEST_PATH_IMAGE164
=3 m/s以及
Figure 312280DEST_PATH_IMAGE164
=6 m/s条件下的运动危险性概率图。
综上，本发明所提出的一种三维滑坡运动危险性概率评价方法，以量化运动强度超过不同阈值的概率，从而生成用于单体滑坡定量风险评估（QRA）的运动危险概率图。利用CJ8号滑坡第一次滑动所测得的多个最终堆积深度，采用基于MCMC的贝叶斯反分析方法对动态数值模型的输入参数进行了校准。应用校准后的参数（即后验分布）估计CJ8号滑坡第二次滑动各网格单元的最大运动高度和最大运动速度的超越概率，进行运动危险性概率图绘制。本发明的主要结论如下：1）该方法有效地综合了先验信息（即输入参数的先验分布）、观测信息（即不同网格位置的最终堆积深度）和模型偏差系数，得到了输入参数的后验分布，充分利用了现有信息。后验分布提供了有关校准输入参数不确定性的信息，MCMC的后验样本可直接用于概率运动危险评价。2）利用三维动态数值模型Massflow将输入参数中的不确定性传播到预测的运动强度中。生成了超越不同最大运动高度阈值（0 m，0.75 m，1.5 m）和不同的最大运动速度阈值（3 m/s和6 m/s）的概率运动危险性图，可应用于QRA中的概率易损性评估。3)基于序贯自适应采样策略的多响应克里金代理模型，大大降低了反演和正演预测的计算量，从149d（基于动态数值模型）降低到1.3d，并且精度可接受，具有工程应用价值。4）概率方法的一个普遍局限性是，它继续承受着相对较大的计算成本。CJ8号的第一次和第二次滑动的单次质量流正向运行都需要在工作站上至少2分钟[Intel（R）Xeon（R）E5-2630 v4 CPU@2.2 GHz和32 GB RAM]。然而，为了获得后验统计量的稳健估计，通常需要60000次MCMC模拟，这意味着需要60000次运动模拟。因此，直接基于动态数值模型的计算时间约为120000min（约83d），这显然是不可接受的。此外，概率运动预测还需要50000次数值模型评估，以获得超越概率的稳定估计；其基于动态数值模型的计算成本（约69d）也是不可接受的。而在本发明中，克里金代理模型解决了动态数值模型的局限性，将概率反分析和概率运动危险性预测的时间分别减少到0.5d和0.8d。
本说明书中的各个实施例均采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似的部分互相参见即可。
本领域内的技术人员应明白，本发明实施例的实施例可提供为方法、装置、或计算机程序产品。因此，本发明实施例可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明实施例是参照根据本发明实施例的方法、终端设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理终端设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理终端设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理终端设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理终端设备上，使得在计算机或其他可编程终端设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程终端设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
尽管已描述了本发明实施例的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例做出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明实施例范围的所有变更和修改。
最后，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者终端设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者终端设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者终端设备中还存在另外的相同要素。
以上对本发明所提供的一种三维滑坡运动危险性概率评价方法，进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种三维滑坡运动危险性概率评价方法，其特征在于，所述方法包括：
步骤S1：确定潜在滑坡对应的过去滑坡，并在所述过去滑坡的滑坡堆积区内选择多个观测点；
步骤S2：对所述多个观测点中的每一个观测点构造克里金代理模型，并基于每一个观测点的克里金代理模型、贝叶斯定理以及差分进化适应大都会市算法，对预先确定的随机变量的先验信息进行反分析，得到所述随机变量的后验分布；
步骤S3：将所述潜在滑坡预估的滑坡堆积区网格化，选取所述潜在滑坡整个运动过程中各网格位置的最大运动高度和最大运动速度作为评价指标，建立所述潜在滑坡在每一个网格的最大运动高度功能函数和最大运动速度功能函数，并计算得到所述潜在滑坡的最大运动高度克里金代理模型矩阵和最大运动速度克里金代理模型矩阵；
步骤S4：将所述随机变量的后验分布输入到所述潜在滑坡的最大运动高度克里金代理模型矩阵中，输出得到每一个网格的最大运动高度，以及将所述随机变量的后验分布输入到所述最大运动速度克里金代理模型矩阵中，输出得到每一个网格的最大运动速度，并对每一个网格的最大运动高度超越概率以及最大运动速度超越概率分别进行计算，获得所述潜在滑坡的最大运动高度超越概率矩阵以及所述潜在滑坡的最大运动速度超越概率矩阵；
步骤S5：在三维地理信息系统平台上，根据所述潜在滑坡的最大运动高度超越概率矩阵以及所述潜在滑坡的最大运动速度超越概率矩阵，生成用于单体滑坡定量风险评价的概率运动危险性图。
2.根据权利要求1所述的三维滑坡运动危险性概率评价方法，其特征在于，所述多个观测点的个数为K个；所述步骤S2包括：
子步骤S2-1：对所述K个观测点中的每一个观测点构造动态数值模型，并建立每一个观测点的克里金代理模型，将所述动态数值模型的输入参数作为随机变量θ，利用同一个观测点的克里金代理模型计算该观测点处所述随机变量对应的最终堆积深度模拟值，获得所述K个观测点的最终堆积深度模拟值
Figure 761171DEST_PATH_IMAGE001
，其中，
Figure 968162DEST_PATH_IMAGE002
表示第K个观测点的最终堆积深度模拟值；
子步骤S2-2：利用设定的模型偏差系数
Figure 81611DEST_PATH_IMAGE003
和所述K个观测点的最终堆积深度模拟值
Figure 6842DEST_PATH_IMAGE004
，计算得到所述K个观测点的最终堆积深度观测值
Figure 168833DEST_PATH_IMAGE005
，其中，
Figure 495909DEST_PATH_IMAGE006
；
子步骤S2-3：基于所述过去滑坡的K个观测点对应的最终堆积深度观测值，通过贝叶斯定理对所述随机变量的先验分布
Figure 780260DEST_PATH_IMAGE007
进行反分析，利用差分进化适应大都会市算法对反分析结果进行采样，得到所述随机变量的后验分布
Figure 927207DEST_PATH_IMAGE008
；
Figure 955206DEST_PATH_IMAGE009
；
上式中，
Figure 136789DEST_PATH_IMAGE010
是似然函数，用于表征现场测量数据与动态数值模型模拟结果之间的匹配程度，表示为所述随机变量的联合条件概率密度函数；其中，在假设不同网格位置之间的模型偏差系数
Figure 778992DEST_PATH_IMAGE011
是不相关的情况下，所述
Figure 678814DEST_PATH_IMAGE012
表示为：
Figure 244925DEST_PATH_IMAGE013
；
上式中，
Figure 281014DEST_PATH_IMAGE014
表示
Figure 907168DEST_PATH_IMAGE015
的概率密度函数在
Figure 28707DEST_PATH_IMAGE016
对应的概率密度。
3.根据权利要求2所述的三维滑坡运动危险性概率评价方法，其特征在于，所述子步骤S2-1中的对所述K个观测点中的每一个观测点构造动态数值模型，并建立每一个观测点的克里金代理模型，包括：
子步骤S2-101：使用摩擦模型计算基底阻力
Figure 336192DEST_PATH_IMAGE017
，即
Figure 492367DEST_PATH_IMAGE018
；
上式中，𝜎 = ρgh，𝜎表示基底正应力，ρ表示滑坡物质的密度，h表示滑坡物质的运动高度，g代表重力；
Figure 289421DEST_PATH_IMAGE019
是有效残余摩擦角；r u 是孔隙压力比系数，反映了基底材料的液化程度，取值范围为0~1；
子步骤S2-102：基于所述摩擦模型，利用质量流模型对所述过去滑坡进行运动模拟，对所述K个观测点中的每一个观测点构造动态数值模型；其中，所述质量流模型是基于三维纳维-斯托克斯方程沿深度积分的连续介质模型；
子步骤S2-103：将所述动态数值模型的输入参数
Figure 898257DEST_PATH_IMAGE020
和r u 定义为所述随机变量
Figure 806170DEST_PATH_IMAGE021
，并通过均匀设计对所述随机变量进行采样，得到输入数据集，所述输入数据集包括多个初始采样点；
子步骤S2-104：将所述输入数据集的每个初始采样点代入每一个观测点的动态数值模型中，以分别输出每一个观测点位置对应的响应值，基于所述输入数据集以及每一个观测点位置对应的响应值，对克里金模型进行训练，得到每一个观测点的初始克里金代理模型；
子步骤S2-105：从候选样本池中选择一个新样本点，其中，所述候选样本池
Figure 816852DEST_PATH_IMAGE022
中的候选样本点通过拉丁超立方体抽样生成；j=1，2，⋯, s；其中s是候选样本点的数量；
子步骤S2-106：将所述新样本点添加到所述输入数据集中，以对每一个观测点的初始克里金代理模型重新训练，持续迭代，直到该观测点的克里金代理模型满足预先设定的收敛准则，得到最终该观测点的克里金代理模型。
4.根据权利要求3所述的三维滑坡运动危险性概率评价方法，其特征在于，所述子步骤S2-105中的从候选样本池中选择一个新样本点，包括：
通过最大化质量参数
Figure 784808DEST_PATH_IMAGE023
从候选样本池中选择一个新样本点，其中，
Figure 67890DEST_PATH_IMAGE024
的表达式如下：
Figure 779495DEST_PATH_IMAGE025
；
上式中，
Figure 910262DEST_PATH_IMAGE026
是第j个候选样本点与当前设计样本点之间的最小欧氏距离，其中，当前设计样本点包括所述初始采样点和所述新样本点；
Figure 49119DEST_PATH_IMAGE027
是第j个候选样本点处的预测方差，由每一个观测点的第i个克里金代理模型进行估计获得；
w i 为权重因子，表示每一个观测点的第i个克里金代理模型对总预测方差的贡献，其中：
Figure 366968DEST_PATH_IMAGE028
；
上式中，
Figure 882263DEST_PATH_IMAGE029
是每一个观测点的第i个克里金代理模型的预测方差的平均值，
Figure 805219DEST_PATH_IMAGE030
是每一个观测点的第i个克里金代理模型的最大预测方差，p是每一个观测点的克里金代理模型的数目。
5.根据权利要求3所述的三维滑坡运动危险性概率评价方法，其特征在于，所述子步骤S2-106包括：
在α≥0.001时，将所述新样本点添加到所述输入数据集中，以对每一个观测点的初始克里金代理模型重新训练，持续迭代，直到该观测点的克里金代理模型满足序贯抽样算法的收敛准则，得到最终该观测点的克里金代理模型；
其中，α定义为第d - 49次迭代到第d次迭代的确定性系数R 2的变化程度，d为自适应序贯抽样算法的迭代次数且d ≥ 50。
6.根据权利要求5所述的三维滑坡运动危险性概率评价方法，其特征在于，在得到最终该观测点的克里金代理模型之后，所述方法还包括：
利用拉丁超立方体抽样生成多个验证样本点；
针对每一个观测点，分别使用为该观测点构造的动态数值模型和最终该观测点的克里金代理模型计算每个验证样本点对应的最终堆积深度；
将该观测点的动态数值模型计算得到的每个验证样本点对应的最终堆积深度与最终该观测点的克里金代理模型计算得到的每个验证样本点对应的最终堆积深度进行比较，利用所述R 2评价最终该观测点的克里金代理模型的预测质量；其中，所述R 2用α来衡量最终该观测点的克里金代理模型的预测准确度随迭代次数的增加而改善的程度。
7.根据权利要求1所述的三维滑坡运动危险性概率评价方法，其特征在于，所述步骤S3包括：
子步骤S3-1：将所述潜在滑坡预估的滑坡堆积区网格化，选取所述潜在滑坡整个运动过程中各网格位置的最大运动高度和最大运动速度作为评价指标，建立所述潜在滑坡在每一个网格的最大运动高度功能函数
Figure 114978DEST_PATH_IMAGE031
和最大运动速度功能函数
Figure 920123DEST_PATH_IMAGE032
：
Figure 973530DEST_PATH_IMAGE033
；
Figure 813310DEST_PATH_IMAGE034
；
Figure 293970DEST_PATH_IMAGE035
；
上式中，n和m分别是网格在所在位置的列数和行数，N和M分别是所述潜在滑坡预估的滑坡堆积区网格化后的列总数和行总数；
Figure 851990DEST_PATH_IMAGE036
为所述潜在滑坡在（n, m）网格处最大运动高度克里金代理模型输出的最大运动高度，
Figure 630459DEST_PATH_IMAGE037
为所述潜在滑坡在（n, m）网格处最大运动速度克里金代理模型输出的最大运动速度；
Figure 590325DEST_PATH_IMAGE038
和
Figure 241886DEST_PATH_IMAGE039
分别是最大运动高度的阈值和最大运动速度的阈值；
子步骤S3-2：基于所述潜在滑坡在每一个网格的最大运动高度功能函数
Figure 21623DEST_PATH_IMAGE040
，计算得到所述潜在滑坡的最大运动高度克里金代理模型矩阵
Figure 416832DEST_PATH_IMAGE041
：
Figure 231204DEST_PATH_IMAGE042
；
子步骤S3-3：基于所述潜在滑坡在每一个网格的最大运动速度功能函数
Figure 991350DEST_PATH_IMAGE043
，计算得到所述潜在滑坡的最大运动速度克里金代理模型矩阵
Figure 258383DEST_PATH_IMAGE044
：
Figure 457283DEST_PATH_IMAGE045
。
8.根据权利要求7所述的三维滑坡运动危险性概率评价方法，其特征在于，在通过所述子步骤S3-2计算得到所述潜在滑坡的最大运动高度克里金代理模型矩阵
Figure 126162DEST_PATH_IMAGE046
以及通过所述子步骤S3-3计算得到所述潜在滑坡的最大运动速度克里金代理模型矩阵
Figure 119526DEST_PATH_IMAGE047
之后，所述步骤S4包括：
将所述随机变量的后验分布输入到所述潜在滑坡的最大运动高度克里金代理模型矩阵
Figure 873855DEST_PATH_IMAGE048
，输出得到每一个网格的最大运动高度；
每个网格位置（n，m）处的最大运动高度超越概率计算为
Figure 63397DEST_PATH_IMAGE049
：
Figure 586782DEST_PATH_IMAGE050
；
上式中，T为后验分布中的后验样本数，θ t 是由马尔可夫链蒙特卡罗模拟采样获得的第t个后验样本；当
Figure 751047DEST_PATH_IMAGE051
时，表示
Figure 992673DEST_PATH_IMAGE052
发生；
所述潜在滑坡的最大运动高度超越概率矩阵为：
Figure 533376DEST_PATH_IMAGE053
；
以及，
将所述随机变量的后验分布输入到所述潜在滑坡的最大运动速度克里金代理模型矩阵
Figure 911267DEST_PATH_IMAGE054
，输出得到每一个网格的最大运动速度；
每个网格位置（n，m）处的最大运动速度超越概率计算为
Figure 512013DEST_PATH_IMAGE055
：
Figure 913039DEST_PATH_IMAGE056
；
上式中，T为后验分布中的后验样本数，θ t 是由马尔可夫链蒙特卡罗模拟采样获得的第t个后验样本；当
Figure 257432DEST_PATH_IMAGE057
时，表示
Figure 489830DEST_PATH_IMAGE058
发生；
所述潜在滑坡的最大运动速度超越概率矩阵为：
Figure 261477DEST_PATH_IMAGE059
。
9.根据权利要求1所述的三维滑坡运动危险性概率评价方法，其特征在于，所述步骤S5包括：
利用三维地理信息系统平台将得到的所述潜在滑坡的最大运动高度超越概率矩阵以及所述潜在滑坡的最大运动速度超越概率矩阵投影到栅格文件上；其中，所述栅格文件用于在所述三维地理信息系统平台上可视化；
将所述栅格文件叠加在所述潜在滑坡的数字高程模型或卫星图像地图上，得到用于单体滑坡定量风险评价的概率运动危险性图。
10.根据权利要求1所述的三维滑坡运动危险性概率评价方法，其特征在于，所述多个观测点的选取来源包括：
临近所述过去滑坡的滑坡堆积区底边界的点；
临近所述过去滑坡的滑坡堆积区最顶部边界的点；
临近所述过去滑坡的滑坡堆积区左右边界的点。
CN202110878125.4A 2021-08-02 2021-08-02 一种三维滑坡运动危险性概率评价方法 Pending CN113343504A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110878125.4A CN113343504A (zh) 2021-08-02 2021-08-02 一种三维滑坡运动危险性概率评价方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110878125.4A CN113343504A (zh) 2021-08-02 2021-08-02 一种三维滑坡运动危险性概率评价方法

Publications (1)

Publication Number Publication Date
CN113343504A true CN113343504A (zh) 2021-09-03

Family

ID=77480526

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110878125.4A Pending CN113343504A (zh) 2021-08-02 2021-08-02 一种三维滑坡运动危险性概率评价方法

Country Status (1)

Country Link
CN (1) CN113343504A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114626008A (zh) * 2022-03-15 2022-06-14 中铁二院工程集团有限责任公司 一种基于幂相关随机过程的铁路路基沉降预测方法和装置
CN114781731A (zh) * 2022-04-26 2022-07-22 成都理工大学 基于贝叶斯理论的滑坡运动距离超越概率预测方法和系统
CN115162376A (zh) * 2022-07-28 2022-10-11 河海大学 基于gis平台进行滑坡过程概率分析的滑坡危险性区划方法
CN115510578A (zh) * 2022-09-26 2022-12-23 成都理工大学 基于InSAR近实时监测的滑坡失稳时间概率预测方法及产品

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101770027A (zh) * 2010-02-05 2010-07-07 河海大学 基于InSAR与GPS数据融合的地表三维形变监测方法
CN112288169A (zh) * 2020-10-30 2021-01-29 成都理工大学 二维滑坡动力学参数概率反分析与滑程超越概率评估方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101770027A (zh) * 2010-02-05 2010-07-07 河海大学 基于InSAR与GPS数据融合的地表三维形变监测方法
CN112288169A (zh) * 2020-10-30 2021-01-29 成都理工大学 二维滑坡动力学参数概率反分析与滑程超越概率评估方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
PENG ZENG 等: "3D probabilistic landslide run-out hazard evaluation for quantitative risk assessment purposes", 《ENGINEERING GEOLOGY》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114626008A (zh) * 2022-03-15 2022-06-14 中铁二院工程集团有限责任公司 一种基于幂相关随机过程的铁路路基沉降预测方法和装置
CN114626008B (zh) * 2022-03-15 2023-03-21 中铁二院工程集团有限责任公司 一种基于幂相关随机过程的铁路路基沉降预测方法和装置
CN114781731A (zh) * 2022-04-26 2022-07-22 成都理工大学 基于贝叶斯理论的滑坡运动距离超越概率预测方法和系统
CN114781731B (zh) * 2022-04-26 2023-04-18 成都理工大学 基于贝叶斯理论的滑坡运动距离超越概率预测方法和系统
CN115162376A (zh) * 2022-07-28 2022-10-11 河海大学 基于gis平台进行滑坡过程概率分析的滑坡危险性区划方法
CN115510578A (zh) * 2022-09-26 2022-12-23 成都理工大学 基于InSAR近实时监测的滑坡失稳时间概率预测方法及产品
CN115510578B (zh) * 2022-09-26 2023-07-14 成都理工大学 基于InSAR近实时监测的滑坡失稳时间概率预测方法及产品

Similar Documents

Publication Publication Date Title
CN113343504A (zh) 一种三维滑坡运动危险性概率评价方法
Moghaddam et al. Developing comparative mathematic models, BN and ANN for forecasting of groundwater levels
Callaghan et al. Statistical simulation of wave climate and extreme beach erosion
Charvet et al. A multivariate generalized linear tsunami fragility model for Kesennuma City based on maximum flow depths, velocities and debris impact, with evaluation of predictive accuracy
Moradkhani et al. Evolution of ensemble data assimilation for uncertainty quantification using the particle filter‐Markov chain Monte Carlo method
CN110837669B (zh) 基于多源异构数据融合的滑坡不确定模型动态构建方法
Vergara et al. Estimating a-priori kinematic wave model parameters based on regionalization for flash flood forecasting in the Conterminous United States
He et al. Assessing hydrological model predictive uncertainty using stochastically generated geological models
Dawson et al. Adaptive importance sampling for risk analysis of complex infrastructure systems
CN108563837B (zh) 一种冲积河流水沙模型的模型参数实时校正方法和系统
Carreau et al. Stochastic downscaling of precipitation with neural network conditional mixture models
Sun et al. From probabilistic back analyses to probabilistic run-out predictions of landslides: A case study of Heifangtai terrace, Gansu Province, China
KR102061500B1 (ko) 인공신경망 기반 지진해일 예측 시스템
CN112288169B (zh) 二维滑坡动力学参数概率反分析与滑程超越概率评估方法
Zeng et al. 3D probabilistic landslide run-out hazard evaluation for quantitative risk assessment purposes
Yoshida et al. Optimal sampling placement in a Gaussian random field based on value of information
Oumeraci et al. XtremRisK—integrated flood risk analysis for extreme storm surges at open coasts and in estuaries: methodology, key results and lessons learned
CN110569477A (zh) 一种基于粒子群优化算法的滑坡区间稳定性分析方法
Schmidinger et al. Validation of uncertainty predictions in digital soil mapping
Svalova et al. Emulating computer experiments of transport infrastructure slope stability using Gaussian processes and Bayesian inference
Falk et al. Estimating uncertainty in the revised universal soil loss equation via Bayesian melding
CN117077554B (zh) 一种基于ConvGRU的三维咸潮预报方法
Gong et al. The influence of seismic displacement models on spatial prediction of regional earthquake-induced landslides
Fu et al. Reliability of the prediction model for landslide displacement with step-like behavior
Amaya et al. Hydrogeological multiple-point statistics inversion by adaptive sequential Monte Carlo

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
RJ01 Rejection of invention patent application after publication

Application publication date: 20210903

RJ01 Rejection of invention patent application after publication