一种基于地表形变和孕灾环境条件的滑坡易发性评价方法

1.本发明涉及地质灾害评价技术领域，具体涉及一种基于地表形变和孕灾环境条件的滑坡易发性评价方法。


背景技术：

2.滑坡是环境参数变化的结果。在极端气候事件的影响下，世界范围内发生的滑坡数量越来越多，造成了重大的经济和人员损失。因此开展滑坡易发性评价是十分有必要的，可为地质灾害风险评价提供依据。
3.目前地质灾害易发性评价的主要方法可以分为经验模型、统计模型和确定性模型三种基本类型，其中统计模型克服了经验模型中权重赋值的强烈主观性，并且避免了确定性模型中对地质灾害发生机制知识和详细地质环境数据的需求，因此适合并被广泛应用于区域尺度的地质灾害易发性定量评价的工作。统计模型在广义上主要包括频率比法、证据权法、逻辑回归法、层次分析法、人工神经网络法和随机森林法等。
4.其中统计模型的输入主要分为两部分内容，一部分为历史滑坡的编目数据，另一部分为孕灾环境条件因子。孕灾环境条件因子主要包括地层岩性、构造、坡度、坡向、高程、距水系距离、植被指数等，这些输入数据均为静态的，缺乏动态的数据控制。然而滑坡是一个动态变化的地质现象，缺乏动态因子的输入将严重影响滑坡易发性评价结果的精度。


技术实现要素：

5.本发明的目的之一在于提供一种基于地表形变和孕灾环境条件的滑坡易发性评价方法，该评价方法能够提高目标区域的滑坡易发性评价结果的精度。
6.为实现上述目的，本发明的技术方案为：
7.一种基于地表形变和孕灾环境条件的滑坡易发性评价方法，利用动态因子进行滑坡易发性评价，该评价方法具体包括以下步骤：
8.获取滑坡易发性评价目标区域的滑坡静态因子和动态因子；
9.根据所述滑坡静态因子，建立随机森林模型进行目标区域的滑坡易发性评价，获得第一易发性评价结果；
10.建立第一易发性评价结果与动态因子的列联矩阵，获得第二易发性评价结果；
11.通过gis对第一易发性评价结果和第二易发性评价结果进行对比整合。
12.进一步地，通过接受者操作特征曲线(roc)和曲线下方面积(auc)对所述第一易发性评价结果进行精度分析。
13.进一步地，获取滑坡易发性评价目标区域的历史滑坡数据、遥感解释和现场调查数据；
14.根据所述历史滑坡数据、遥感解释和现场调查数据获取滑坡易发性评价目标区域的滑坡分布信息；
15.将所述滑坡分布信息分为滑坡训练数据集和滑坡验证数据集；
16.将所述滑坡训练数据集和滑坡验证数据集输入随机森立模型中进行滑坡易发性评价。
17.进一步地，所述滑坡易发性评价目标区域的滑坡静态因子包括坡度、坡向、剖面曲率、高程、距断裂距离、地形湿度指数(twi)、归一化植被指数(ndvi)和岩性；
18.通过多重共线性分析对所述滑坡易发性评价目标区域的滑坡静态因子进行相关性分析。
19.进一步地，所述滑坡易发性评价目标区域的动态因子为毫米级别的地表形变速率。
20.进一步地，所述地表形变速率的获取步骤包括：
21.获取sentinel-1和radarsat-2的合成孔径雷达(sar)数据；
22.利用差分干涉测量小基线集时序分析技术(sbas-insar)，获取所述地表形变速率。
23.进一步地，对所述地表形变速率进行校正，获得多源数据的los向形变速率和最优结果的斜坡向形变速率vslop。
24.进一步地，根据雷达数据和斜坡数据对所述los向形变速率进行转换，获得斜坡向形变速率vslop。
25.进一步地，建立第一易发性评价结果与动态因子的列联矩阵时，根据斜坡向形变速率的标准偏差(δ＝13mm/yr)确定地表形变速率区间；
26.确定地表形变速率区间后，将易发程度从1级增加到5级。
27.进一步地，利用差分干涉测量小基线集时序分析技术(sbas-insar)获取目标区域毫米级别的地表形变速率时，基于所述滑坡易发性评价目标区域的斜坡形变特征进行分级。
28.与现有技术相比，本发明具有如下优点：
29.本发明通过sbas-insar技术获取动态因子即地表形变速率，再使用传统的滑坡易发性评价方法获取滑坡易发性评价结果；使用列联矩阵将地表形变速率与滑坡易发性评价结果整合到一起，得到精度更高的易发性结果；最后对比分析整合前后的易发性评价结果，能够更加清楚地显示出本发明中地表形变速率对滑坡易发性评价结果的影响。
附图说明
30.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍。显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。
31.图1为本发明中的滑坡易发性评价方法流程图；
32.图2为本发明中滑坡静态因子中的高程图；
33.图3为本发明中滑坡静态因子中的坡度图；
34.图4为本发明中滑坡静态因子中的坡向图；
35.图5为本发明中滑坡静态因子中的剖面曲率图；
36.图6为本发明中滑坡静态因子中的twi图；
37.图7为本发明中滑坡静态因子中的ndvi图；
38.图8为本发明中滑坡静态因子中的距断层的距离图；
39.图9为本发明中滑坡静态因子中的岩性图；
40.图10为本发明中的第一易发性评价结果示意图；
41.图11为本发明中的目标区域的地表形变速率示意图，其中a为sentinel-1数据获取的卫星los向形变速率，b为radarsat-2数据获取的卫星los向形变速率，c为对比后选取的斜坡向形变速率；
42.图12为本发明中的目标区域第一评价结果与第二评价结果对比示意图，其中a为应用校正矩阵后的第二滑坡易发性评价结果，b为每个单元的原始滑坡易发性等级和新的滑坡易发性等级之间的差异；标号为1的方框区域对应日通区域，标号为2的方框区域对应雪如区域；
43.图13为本发明中通过接受者特征曲线(roc)和曲线下方面积(auc)对第一易发性评价结果进行精度分析的示意图。
具体实施方式
44.为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
45.所举实施例是为了更好地对本发明进行说明，但并不是本发明的内容仅局限于所举实施例。所以熟悉本领域的技术人员根据上述发明内容对实施方案进行非本质的改进和调整，仍属于本发明的保护范围。
46.实施例1
47.广域的滑坡易发性目前最常用的评价模型主要为统计型模型，以随机森林模型为例，主要通过统计现有的滑坡数据，将有利于滑坡发生的数据用于计算预测区域的滑坡易发性。随机森林(rf)是一种结合装袋法(bagging)生成多个独立的样本集和多棵分类回归树来进行预测的集成学习方法，结果由投票得分最多或取平均决定。其主要思想在于多个弱分类器经过一定策略进行组合，形成一个较单一分类器预测性能更加优越的集成模型。随机森林算法能够处理大维数、大数据量的数据集，具有很高的泛化能力。与其他统计学习方法相比，该方法不易过拟合，在计算无明显改善的前提下，提高了预测精度。此外，在数据缺失和不平衡的情况下，结果是稳健的。它在数据分类和管理中得到了广泛的应用，并在滑坡易发性评价中得到了很好的应用。在本实施例中，我们使用r语言环境中的随机森林模块来评价滑坡的易发性。
48.本实施例通过毫米级地表变形速率，有效解决滑坡易发性评价动态因子输入的问题，有效提高滑坡易发性评价结果的精度，本实施例中的滑坡易发性评价方法流程参见图1。
49.insar技术作为监测地质灾害，特别是滑坡的工具，在过去的二十年里，受到了越来越多的关注。利用星载雷达手段进行滑坡易发性评价，是改善滑坡监测与识别的可靠途径，尤其是在区域尺度上识别滑坡的能力，提供了滑坡现今活动性的宝贵信息。由于滑坡易
发性的现场测绘是一项复杂而耗时的任务，涉及收集现场数据和分析，定期重复的现场调查来获取滑坡新近的变化被认为是不切实际的，特别是大规模开展的情况下。为了解决这个问题，本实施例中利用小基线子集insar(sbas insar)来获取滑坡的活动性，作为易发性评价的指标之一。该技术克服了时间非相干的局限性，避免了传统干涉测量方法中的长时间分离、空间非相干和大气效应，产生了在时间和空间上更加连续的陆地形变结果。该方法适用于长序列的缓慢线性和非线性变形监测，广泛应用于地面沉降、地震、活动断层，以及蠕变和滑坡等边坡失稳的识别与监测。
50.本实施中针对西藏昌都市澜沧江卡若至卡贡段的区域进行滑坡易发性评价，主要为以下几个步骤：
51.s1、获取西藏昌都市澜沧江卡若至卡贡段区域的滑坡静态因子
52.滑坡静态因子为目标区域的坡度、坡向、剖面曲率、高程、距断裂距离、地形湿度指数(twi)、归一化植被指数(ndvi)和岩性，数据来源为区域地质图、数字高程模型、遥感影像。提取的高程图参见图2，坡度图参见图3，坡向图参见图4，剖面图参见图5，twi图参见图6，ndvi图参见图7，距断层的距离图参见图8，岩性图参见图9。在图9中，1代表层理发育的砾岩、砂岩和泥岩(岩组2)，2代表灰岩、板岩、白云岩和大理岩(岩组3)，3为片麻岩、变粒岩和石英片岩(岩组4)，4为块状结构的花岗岩、闪长岩以及岩浆岩岩脉(岩组5)，5代表松散结构的砾石、砂土、碎石土(岩组12)。
53.再通过多重共线性分析对所述滑坡易发性评价目标区域的滑坡静态因子进行相关性分析，分析结果参见表1。
54.表1滑坡孕灾环境静态因子的多重共线性分析
[0055][0056]
由表1中的相关性分析结果可知，vif值最高为2.842，最低tol值为0.985，表明8个滑坡条件因子之间不存在多重共线性。利用多重共线性对静态因子进行相关性分析，能够进一步提升评价结果的精确性。
[0057]
s2、获取第一易发性评价结果
[0058]
获取滑坡易发性评价目标区域的历史滑坡数据、遥感解释和现场调查数据；然后根据历史滑坡数据、遥感解释和现场调查数据获取滑坡易发性评价目标区域的滑坡分布信息；将滑坡分布信息分为滑坡训练数据集和滑坡验证数据集。
[0059]
将滑坡训练数据集、滑坡验证数据集、步骤s1中的滑坡静态因子输入随机森立模型中，进行目标区域的滑坡易发性评价，该评价结果为第一易发性评价结果，本实施例中的第一易发性评价结果图参见图10。通过接受者特征曲线(roc)和曲线下方面积(auc)对第一易发性评价结果进行精度分析，参见图13，可得曲线下方面积(auc)值为0.897，满足评价精度需求。
[0060]
s3、获得滑坡动态因子
[0061]
选取sentinel-1和radarsat-2的合成孔径雷达(sar)数据，利用差分干涉测量小基线集时序分析技术(sbas-insar)，获取目标区域毫米级别的地表形变速率。对地表形变
速率进行校正，获得多源数据的los向形变速率和最优结果的斜坡向形变速率vslop。
[0062]
对los向的形变速率进行转换，获得斜坡向形变速率vslop，转换方式为:其中index＝nlos
·
nslope；nlos＝(sinθsinαs，-sinθcosαs，-cosθ)；在上述公式中：vslop为斜坡向形变速率；vlos为los向的形变速率；α为斜坡的坡向，单位为
°
；为斜坡的坡度，单位为
°
；θ为雷达入射角，单位为
°
；αs为卫星轨道方向与正北方向的夹角，即雷达卫星飞行方向，升轨数据为负，降轨数据为正。
[0063]
本实施例中目标区域的地表形变速率分布参见图11，其中a为sentinel-1数据获取的卫星los向形变速率，b为radarsat-2数据获取的卫星los向形变速率，c为对比后选取的斜坡向形变速率。
[0064]
s4、获取第二易发性评价结果
[0065]
建立滑坡易发性与地表形变速率之间的列联矩阵来整合第一易发性评价结果，最终获取第二易发性评价结果，列联矩阵整合对照表参见表2，第二易发形评价结果见图12a。
[0066]
表2滑坡易发性和形变速率整合对照表
[0067][0068]
其中形变速率区间是根据斜坡向形变速率的标准偏差(δ＝13mm/yr)确定的。确定了速率区间，将易发程度从1级增加到5级。速率值越大，滑坡易发程度越高。易发程度不能超过5。这种整合降低了错误分类的概率，并可以将易发程度低但即将发生滑坡的区域的易发性修正为高易发性。
[0069]
其中insar获取的形变量等级分区，基于野外现场调查的斜坡形变特征进行分级。
[0070]
s5、对比第一易发性评价结果与第二易发性评价结果
[0071]
通过gis对比整合前后的易发性评价结果，对比结果参见图12，其中a为应用校正矩阵后获得的第二滑坡易发性评价结果，b为每个单元的原始滑坡易发性等级和新的滑坡易发性等级之间的差异；标号为1的方框区域对应日通区域，标号为2的区域对应雪如区域。由图12中可以看出，滑坡易发性明显增强的地区有多个区域，其中日通和雪如区域变化最为明显。新的易发性评价图中各级单元百分比变化如表3所示，1类栅格单元占比31.01％，2类栅格单元占比26.55％，3类栅格单元占比19.17％，4类栅格单元占比15.64％，5类栅格单元占比7.63％。通过对比新老评价结果的差异，发现有563430个栅格单元的滑坡易发性发生了改变，其中有412563个增加了1个易发性等级，113775个栅格单元增加了2个易发性等
级，36985个栅格单元增加了3个易发性等级，106个栅格单元增加了4个易发性等级。
[0072]
表3新老易发形评价结果的对比结果
[0073][0074][0075]
上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。