一种区域滑坡易发性更新评价方法

技术领域

本发明属于地质灾害防治预警模拟技术领域，为一种区域滑坡易发性更新评价方法。

背景技术

近年来，国内外研究者针对区域尺度潜在滑坡早期识别的难题进行了积极探索，并提出了多种方法。应用最广泛的主要有基于统计模型的区域滑坡易发性评价和基于物理模型的区域斜坡稳定性评价。其中，区域滑坡易发性是指区域内斜坡在地形地貌、地质条件、水文条件和地表覆被等环境因子作用下发生滑坡的概率，主要通过统计分析和机器学习模型研究滑坡与各类环境因子之间的非线性关系，构建区域滑坡易发性评价模型。虽然该方法具有建模简单、效率高、可考虑多种环境因子与滑坡之间的关系等优势，但是缺乏对滑坡发生物理机制的解释，非常依赖滑坡样本数据的质量。另外，当前滑坡易发性结果主要表现为静态的空间分布，忽略了滑坡数据库的多期次性和滑坡易发性结果的时效性，致使滑坡易发性结果实际应用效果较差。基于物理模型的区域斜坡稳定性评价方法，一般是基于无限长边坡模型，采用水文模型和岩土模型对区域斜坡进行稳定性分析，能够模拟降雨过程中地下水的变化，可以更好地解释滑坡形成的物理机制。但是，该方法没有充分融合历史发生滑坡的信息，同时依赖区域尺度相关参数(内聚力、内摩擦角、渗透系数等)的精度，而实际上区域尺度相关参数获取困难、成本高，并且存在较大的不确定性，严重制约了该方法的推广应用。

发明内容

本发明的目的是针对当前区域滑坡易发性评价和基于物理模型区域斜坡稳定性评价方法中存在的不足，提供一种区域滑坡易发性更新评价方法，本发明充分考虑历史滑坡数据的多期次性和滑坡的失稳机制，利用最新滑坡信息，结合贝叶斯更新理论对研究区相关参数进行反演和更新，采用物理模型计算区域斜坡的稳定性和失效概率，进而融合滑坡易发性结果和失效概率结果，实现区域滑坡易发性结果的更新。本发明具有物理含义明确、针对性强、成本低、易操作、结果准确和可视化好等优点。

本发明是通过以下技术方案实现的。

本发明所述的一种区域滑坡易发性更新评价方法，包括以下步骤：

步骤S1、基础数据获取：包括研究区的数字高程数据(DEM)、滑坡数据、初始滑坡易发性分区图、诱发因素数据、初始地下水位、土体参数；

所述的滑坡数据包括历史滑坡数据和最新滑坡数据；所述的土体参数包括有粘聚力c、内摩擦角φ、土体厚度H

步骤S2、参数敏感性分析：基于研究区DEM和土体参数等数据，采用TRIGRS物理模型进行区域斜坡稳定性分析。考虑TRIGRS模型输入参数集X＝[H

步骤S3、敏感参数贝叶斯更新：基于贝叶斯更新理论，将最新滑坡数据作为观测信息，构建滑坡区域安全系数的似然函数；利用未发生滑坡区域的数据作为观测信息，构建稳定区域安全系数的似然函数；最后通过多维联合滑坡区域和稳定区域似然函数对整个区域进行表征，实现对敏感参数的更新；

步骤S4、失效概率计算：选取栅格单元为制图单元，通过将步骤S3中的研究区域联合似然函数应用到整个研究区域，进行多次模拟计算得到各栅格单元输入参数的后验分布：待更新输入参数的均值、标准差的空间分布及其概率密度函数曲线；并对参数先验分布和后验分布进行比较分析，选取合理的输入参数分布对研究区进行稳定性分析和失效概率计算；

步骤S5、失效概率分级制图：根据研究区域的滑坡失效概率计算结果，在ArcGIS10.2软件中通过自然间断法划分为极高、高、中等、低和极低五个等级，并绘制滑坡失效概率等级图；

步骤S6、滑坡灾害预警等级制图：通过矩阵分析法对研究区滑坡失效概率等级图与历史滑坡易发性分区图进行融合，获得研究区预警等级表和区域滑坡易发性更新后的滑坡灾害预警等级分布图，并将整个区域划分为安全区、预警区、危险区。

由上述技术方案可见，本发明提供一种区域滑坡易发性更新评价方法，利用最新滑坡信息对历史区域滑坡易发性等级区划图进行动态更新，首先根据资料收集和现场勘测等方法获得研究区域中基础DEM数据、降雨等诱发因素数据、土体参数等信息。选取栅格单元作为基础制图单元，采用TRIGRS物理模型计算研究区栅格单元的稳定系数，并对模型中输入参数进行敏感性分析。再基于贝叶斯更新方法，将最新滑坡信息作为似然函数中的观测信息分别构建滑坡区域和稳定区域的似然函数，反演并更新更符合工程实际的相关参数，从而更新研究区域的滑坡失效概率分布图，最后利用ArcGIS软件将历史滑坡空间易发性分区结果与滑坡失效概率分布图进行融合得到更新后的区域滑坡灾害预警等级分布图。具有物理含义明确、针对性强、成本低、易操作、结果准确和可视化好等优点，同时克服了基于机器学习模型滑坡易发性评价研究中无法考虑滑坡力学机制和滑坡数据多期次性的缺点，适用性更为广泛。该模型方法能够呈现滑坡灾害的动态分布情况，对区域滑坡灾害治理具有重要的参考价值。

所述步骤S2的区域斜坡稳定性分析为：

采用公式(1)建立区域栅格单元的稳定分析模型，分析计算得到栅格单元不同时段的安全系数，以此为依据判定栅格单元的稳定情况，当栅格单元安全系数≤1处于不稳定状态，栅格单元安全系数>1处于稳定状态。

上式中：F

所述步骤S2中的敏感性分析为：

考虑TRIGRS模型中输入参数集X＝[H

所述步骤S3的具体包括：

依据实际测量值和模拟实现值之间的误差构建似然函数的具体表达式为公式(2)，通过将参数一般正态分布特征转化成标准正态分布特征，公式(2)可转换为公式(3)，如果当误差项大于模拟值与设定阈值的差值时，即y

上式中：f

上式中：

上式中：Φ为标准正态累积分布函数。

采用公式(5)来表示研究区域中栅格单元在诱发因素作用下由稳定状态变为失稳状态的似然函数，即滑坡区域似然函数；采用公式(6)来表示研究区域中栅格单元在诱发因素作用下保持稳定状态的似然函数，即稳定区域似然函数；最后通过公式(7)多维联合滑坡区域和稳定区域中构建的似然函数对整个研究区域进行表征；

式中：L

式中：n

所述步骤S5具体包括：

统计各栅格单元中安全系数≤1的次数，用公式(9)来表征结构可靠度问题，用公式(10)得出各栅格单元的滑坡失效概率；

g

P

式中：F

根据研究区域中最新计算出来的滑坡失效概率在ArcGIS10.2软件中通过自然间断法划分为五类：极高、高、中等、低和极低等级；同理将历史滑坡易发性结果分为极高、高、中等、低和极低五个等级。

所述步骤S6将研究区域融合最新滑坡信息的滑坡灾害预警等级分布图划分为安全区、预警区、危险区，如表1。

表1预警分级表

与现有技术相比，本发明的有益效果是：本发明提供一种区域滑坡易发性更新评价方法，充分考虑滑坡易发性的时效性和滑坡的物理力学机制，基于贝叶斯更新理论和TRIGRS物理模型，将最新滑坡信息作为观测信息构建了研究区的似然函数，对研究区内物理模型中相关参数进行反演更新，进而实现区域内斜坡的安全系数更新和失效概率计算。通过融合区域失效概率分布图和区域滑坡易发性分级图实现区域滑坡易发性实时动态更新的目的。本发明实现了基于数据驱动模型的区域滑坡易发性评价和基于物理驱动的区域斜坡失效概率评价方法的有机融合，具有物理含义明确、针对性强、成本低、易操作、结果准确和可视化好等优点，同时克服了基于机器学习模型滑坡易发性评价研究中无法考虑滑坡力学机制和滑坡数据多期次性的缺点，适用性更为广泛。该模型方法能够呈现滑坡灾害的动态分布情况，对区域滑坡灾害治理具有重要的参考价值。

附图说明

附图1为本发明的模拟方法流程图。

附图2为寻乌县城南滑坡所处的地理位置示意图。

附图3为研究区2006～2020滑坡位置示意图。

附图4为无限长斜坡模型示意图。

附图5为研究区域的历史滑坡易发性示意图。

附图6为研究区域后验失效概率分布示意图。

附图7为研究区域滑坡易发性更新结果示意图。

具体实施方式

为使本发明的目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。

需要说明的是，除非另有定义，本发明所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。在本发明的说明书中所使用的术语只是为了描述具体的实施方式的目的，不是旨在于限制本发明。本发明所使用的术语“及/或”包括一个或多个相关的所列项目的任意的和所有的组合。

本发明所述的一种区域滑坡易发性更新评价方法，具体过程如图1所示，下面结合附图对本发明进行进一步详细说明。

步骤S1、获取基础数据资料；所述基础数据资料包括以下资料：

(1)研究区数字高程模型(DEM)文件；

(2)最新的滑坡发生位置资料，参见附图2和附图3；

(3)TRIGRS模型所需的初始输入参数，包括有粘聚力、内摩擦角、土体厚度、初始地下水位、饱和土体渗透系数和水力扩散系数的均值、标准差和波动范围，所有的输入参数取值都是依据江西省赣州市寻乌县长宁镇土体情况来定义的，随机模型参数的取值和模型误差如下表2。

表2物理模型输入参数分布

(4)本次实例只模拟研究区24小时内的降雨情况，将降雨序列分为两段，即前半段(0～12h)降雨强度为61.72mm/day和后半段(12～24h)为大暴雨工况下的125mm/day。

(5)研究区历史滑坡易发性分区结果图，参见附图5。

步骤S2、在TRIGRS模型中进行栅格单元稳定性和参数敏感性分析。

根据研究区的DEM数据、土体参数等信息，结合TRIGRS模型，参见附图4，采用公式(1)建立区域栅格单元的稳定分析模型，分析计算得到栅格单元不同时段的安全系数。以此为依据判定栅格单元的稳定情况，当计算得到的安全系数≤1，表明栅格单元处于不稳定状态，当计算得到的安全系数>1，表明栅格单元处于稳定状态。

步骤S3、将最新滑坡数据作为观测信息构建滑坡区域安全系数的似然函数，利用未发生滑坡区域的数据作为观测信息构建稳定区域安全系数的似然函数，最后通过多维联合滑坡区域和稳定区域似然函数对整片研究区域进行表征，实现对敏感参数的贝叶斯更新。

具体的说，将最新发生的滑坡作为似然函数中的观测信息，公式(2)是通过参数实际测量值和模拟实现值之间的误差来构建似然函数，将公式(2)转换为公式(3)，如果当误差项大于模拟值与实际测量值的差值时，即y

与以往在区域斜坡中所使用的贝叶斯方法不同的是，本发明在区域滑坡易发性更新研究中构建了两种不同的似然函数形式，即栅格单元在降雨作用下由稳定变为不稳定状态的滑坡区域似然函数；降雨作用前后斜坡都保持稳定状态的稳定区域似然函数，分别用公式(5)和(6)来表示滑坡区域和稳定区域相对应的似然函数。最后通过公式(7)多维联合滑坡区域和稳定区域中构建的似然函数对整个研究区域进行表征。

步骤S4、选取栅格单元作为基本的制图单元，并进行区域斜坡机理分析。每一个栅格单元假设为无限长边坡，代表一个真实的斜坡，通过贝叶斯更新方法获得各栅格单元敏感性较高参数的后验分布，并对参数先验分布和后验分布进行分析比较，选取更加合理的参数分布对研究区进行稳定性分析并计算失效概率。

需要说明的是，在区域滑坡更新研究中，所研究的区域可能包含几万个栅格单元，将每个栅格单元视为一个无限长斜坡模型，都有对应的安全系数，根据栅格单元的精度，将整个研究区域划分为K行L列，获取K×L个栅格单元，通过将贝叶斯更新方法应用到整个研究区域，进行W次模拟计算得到各栅格单元M个敏感性较高参数的后验分布，通过MATLAB将各栅格单元中反演得到的M×K×L×W维后验分布值导入TRIGRS模型中，获得整个研究区栅格单位的安全系数，其中每个栅格单元生成W个安全系数。

在本实施例中，寻乌县城南滑坡作为研究区域，以2006～2020年间发生的6处滑坡作为似然函数中融入的观测信息，栅格单元的精度是30m，整个研究区域为54行53列，一共包含2862个栅格单元，每一个栅格单元代表一个真实的斜坡。对表2中的初始输入参数进行关于安全系数的敏感性分析，剔除部分非敏感的参数，得到敏感性较高的不确定性参数集Z＝[c,φ,H

步骤S5、统计各栅格单元安全系数，得出各栅格单元发生的滑坡失效概率，并根据研究区域的滑坡失效概率在ArcGIS10.2中通过自然间断法划分为五个等级，得到滑坡失效概率等级图；

具体的，统计各栅格单元中安全系数小于或等于1的次数，用公式(9)来表征结构可靠度问题，用公式(10)得出栅格单元发生滑坡的失效概率(图6)。

进一步的，将研究区域中最新计算出来的每个栅格单元的滑坡失效概率在ArcGIS中通过自然间断法划分为五类：极高危险、高危险、中危险、低危险和极低危险，得到滑坡失效概率等级图。

步骤S6、通过矩阵分析法对滑坡失效概率等级图与历史滑坡易发性分区结果图进行融合，从而得到指定研究区域融合最新滑坡信息的滑坡灾害预警等级分布图。

具体的说，通过预警分级表(表1)在ArcGIS将研究区域的滑坡失效概率等级图与历史的滑坡易发性区划图进行融合，得到融合最新滑坡发生信息的指定研究区域的滑坡灾害预警等级分布图(图7)。将其划分为安全区、预警区、危险区。因此当斜坡危险指标达到预警区以上，需要对其重点关注，结合斜坡实际情况进行后续的滑坡灾害治理措施。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。

显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。在本发明中提及实施例意味着，结合实施例描述的特定特征、结构或者特性可以包含在本实施例申请的至少一个实施例中。在说明书中的各个位置出现该短语并不一定均是相同的实施例，也不是与其它实施例互斥的独立的或是备选的实施例。本领域技术人员可以显式地和隐式地理解的是，本发明所描述的实施例可以与其它实施例相结合。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。