基于地形单元的滑坡易发性评判方法
技术领域
本发明涉及一种基于地形单元的滑坡易发性评判方法，属于地球系统科学中的区域地质灾害评价领域。
背景技术
滑坡易发性评价是实现滑坡地质灾害防控的基础，而滑坡致灾因子及滑坡评价单元划分则是滑坡易发性评价中两个最为重要的组成部分。滑坡致灾因子分析是区域滑坡灾害风险评估和防控的首要工作，同时针对滑坡易发性评价单元的划分目前仍存在一些问题，没有较为统一的划分标准，导致滑坡易发性评判方法的可靠性差。
滑坡发生本身是一个较为复杂的过程，构建科学的易发性评价模型对于滑坡易发性评价结果至关重要。
发明内容
针对现有滑坡易发性评判方法可靠性差的问题，本发明提供一种基于地形单元的滑坡易发性评判方法。
本发明的一种基于地形单元的滑坡易发性评判方法，所述方法包括：
S1、确定滑坡点位置及相关的致灾因子数据；
S2、将致灾因子数据进行离散化，获得致灾因子数据的离散区间，每个离散区间对应一个分级类别，获取各个致灾因子在各离散区间分布面积占比AR；
S3、根据各个致灾因子在各离散区间分布面积占比AR，获取权重比指数曲线，采用面积分类比的曲线与权重比指数曲线交点及拐点作为判断标准，将权重比指数曲线与面积分类比的曲线的拐点或交点作为区间的端点值，获得优化后的区间；权重比指数QW为：
其中，H为每一分级类别内的滑坡点数，H0为滑坡点总数；
S4、统计优化后的各离散区间面积及对应滑坡点数，获取滑倾度指数为：
其中，TTs表示滑坡点总数与致灾因子数据所占总体面积的比值，TTa为相应致灾因子对应的优化后的区间中存在的滑坡点总数与区间a面积的比值：
若OS>0，则表示在优化后的区间范围内滑倾度指数OS为正值，滑坡发生的可能性较大；
若OS<0，则表示在优化后的区间范围内滑倾度指数OS为负值，滑坡发生的可能性较小；
若OS≈0，则表示在优化后的区间范围内滑倾度指数OS接近于零，不能确定滑坡是否会发生；
S4、利用地形表面模型DEM构建河网；
S5、计算河网中的连接点得到某一汇水流域的出水点，结合水流方向数据，分别提取出正负地形下的汇水流域，确定汇水流域所有栅格单元位置，将正负地形下的汇水流域合并，得到地形单元；
S6、将每种致灾因子的各个优化区间的滑倾度指数赋给每个地形单元；
S7、构建的随机森林模型，其输入为每个地形单元中的数据，输出为每个地形单元的滑坡易发性概率值。
作为优选，所述S2中，运用等间距方法将连续型的致灾因子数据进行离散化，所使用面积分类比AR为：
其中，Yi为致灾因子数据的第i个分级类别对应面积，Ys为致灾因子数据所占总体面积。
作为优选，所述S4中，利用地形表面模型DEM构建河网，包括：
S41、根据地形表面模型DEM，对河流区域内因地形特点设定不同的模拟洼地深度值，利用模拟洼地深度值对地形表面模型DEM中洼地进行洼地填充；
S42、在地形表面模型DEM中，模拟水流方向：
采用流向均指数模型作为真实水流方向判断依据，流向均指数公式为：
其中，Gi为地形表面模型DEM中两个栅格单元中心之间的直线距离，Gz为地形表面模型DEM中两个栅格单元之间的高程差：
若Max(FI)>0，则表明两个栅格单元之间的流向未定；
若Max(FI)≤0，并且最大值只有一个，则将对应此最大值的方向作为两个栅格单元之间的流向；
若Max(FI)＝0，并且有一个以上的0值，则以这些0值所对应的方向值相加的结果作为两个栅格单元之间的流向；
若Max(FI)<0，并且有一个以上的最大值，则任选一个最大值所对应的方向作为两个栅格单元之间流向；
S43、在地形表面模型DEM中，模拟汇流量：在地形表面模型DEM的每个网格节点处赋予不同深度的水位，按照水流从高处流向低处的自然规律，依据区域地形和水流方向计算各个栅格单元所有方向的水量，获得汇流多向累积量；
S43、根据S43获得的各个栅格单元的汇流多向累积量，设置不同的汇水面积值以模拟出不同条件下的河网，与真实河网进行比较，找出最接近真实河网的汇水面积值，利用该汇水面积值生成河网。
作为优选，汇流多向累积量为：
DU＝D1+D2+D3+...+DN；其中，Dn为某一中心栅格单元第n个流向流来的水流量，n＝1,2,3,...,N,N表示中心栅格单元的水流流向的数量。
作为优选，Dn＝DnM，其中M为1或0，若流向指向目前需要计算汇流多向累积量的栅格单元，则M＝1，否则M＝0。
作为优选，S1中，致灾因子数据包括连续型数据和离散型数据，连续型数据为高程、坡向、距河流距离、距道路距离、剖面曲率、NDVI，离散型数据为岩土体类型、土地利用类型以及斜坡结构。
本发明的有益效果：本发明方法构建了一种基于地形单元的滑坡易发性评判模型，具有如下特点：首先，本发明方法充分考虑多方面致灾影响要素，一定程度上保证发明内容的全面充实。其次，本发明方法能够对滑坡致灾影响要素进行量化分析，从而为滑坡防治工作提供基础。另外，本发明方法在充分考虑地形地貌的基础上选择基于地形的评价单元，保证研究结果的相对可靠。
附图说明
图1是本发明的流程图；
图2是实施例区正地形模拟汇水面积值与对应的沟谷密度曲线图；
图3是基于地形单元的实施例区随机森林模型ROC曲线图；
图4是基于地形单元的实施例区滑坡易发性等级图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
下面结合附图和具体实施例对本发明作进一步说明，但不作为本发明的限定。
本实施方式的基于地形单元的滑坡易发性评判方法，所述方法包括：
步骤一、确定滑坡点位置及相关的致灾因子数据；
步骤二、所选致灾因子分为连续型数据和离散型数据，统一使用离散型数据代入运算。运用等间距方法将致灾因子数据进行离散化，获得致灾因子数据的离散区间，每个离散区间对应一个分级类别，获取各个致灾因子在各离散区间分布面积占比AR；
面积分类比AR为：
其中，Yi为致灾因子数据的第i个分级类别对应面积，Ys为致灾因子数据所占总体面积。
步骤三、根据各个致灾因子在各离散区间分布面积占比AR，获取权重比指数曲线，采用面积分类比曲线与权重比指数曲线交点及拐点作为判断标准，将权重比指数曲线与面积分类比曲线的拐点或交点作为区间的端点值，获得优化后的区间；
权重比指数QW为：
其中，H为每一分级类别内的滑坡点数，H0为滑坡点总数；
S4、统计优化后的各离散区间面积及对应滑坡点数，获取滑倾度指数为：
其中，TTs表示滑坡点总数与致灾因子数据所占总体面积的比值，TTa为相应致灾因子对应的优化后的区间中存在的滑坡点总数与区间a面积的比值：
若OS>0，则表示在优化后的区间范围内滑倾度指数OS为正值，滑坡发生的可能性较大；
若OS<0，则表示在优化后的区间范围内滑倾度指数OS为负值，滑坡发生的可能性较小；
若OS≈0，则表示在优化后的区间范围内滑倾度指数OS接近于零，不能确定滑坡是否会发生；
步骤四、利用地形表面模型DEM构建河网；
步骤五、计算河网中的连接点得到某一汇水流域的出水点，结合水流方向数据，分别提取出正负地形下的汇水流域，确定汇水流域所有栅格单元位置，将正负地形下的汇水流域合并，得到地形单元；
本步骤中，通过计算河流连接点可得到某一汇水流域的出水点即该汇水流域的最低点，结合水流方向数据，分析搜索出该出水点上游所有流过该出水口的栅格，直到该汇水流域所有的栅格都确定了位置，代表已经搜索到流域的边界，分水岭的位置。输入水流方向数据及河流连接点数据输出汇水流域栅格，利用ArcGIS中的栅矢转换功能获得汇水流域的矢量文件。将DEM进行反转可得到负地形再进行上述相同的步骤提取汇水流域，将两次分析所得汇水流域合并即可得到地形单元。
步骤六、将每种致灾因子的各个优化区间的滑倾度指数赋给每个地形单元；
步骤七、构建的随机森林模型，其输入为地形单元中每个栅格单元的数据，输出为每个栅格单元的滑坡易发性概率值。并利用ROC(Receiver Operating CharacteristicCurve)曲线下面积AUC值对模型分类效果进行评价，具体表现为AUC值越大，分类效果越好。
本实施方式的步骤四中，利用地形表面模型DEM构建河网，包括：
步骤四一、根据地形表面模型DEM，对河流区域内因地形特点设定不同的模拟洼地深度值，利用模拟洼地深度值对地形表面模型DEM中洼地进行洼地填充；
DEM是一种较光滑的地形表面模型，由于DEM误差以及一些真实地形的存在，使DEM表面存在一些凹陷的区域，在进行水流方向计算时往往会导致不合理的甚至错误的水流方向，因此在计算水流方向前应先对原始DEM数据进行洼地填充，得到无洼地的DEM。采用模拟洼地深度值对原始DEM数据进行洼地填充，模拟洼地深度值判断公式为
其中，B计为通过考虑研究区地形特点设定不同的洼地深度而得到的模拟洼地分布图，B真为研究区真实地形的洼地分布图。当N的值趋近于1表示模拟洼地深度最接近实际的洼地分布，由此得出可将洼地填平的模拟洼地深度值；
步骤四二、在地形表面模型DEM中，模拟水流方向：
采用流向均指数模型作为真实水流方向判断依据，流向均指数公式为：
其中，Gi为地形表面模型DEM中两个栅格单元中心之间的直线距离，Gz为地形表面模型DEM中两个栅格单元之间的高程差：
若Max(FI)>0，则表明两个栅格单元之间的流向未定；
若Max(FI)≤0，并且最大值只有一个，则将对应此最大值的方向作为两个栅格单元之间的流向；
若Max(FI)＝0，并且有一个以上的0值，则以这些0值所对应的方向值相加的结果作为两个栅格单元之间的流向；
若Max(FI)<0，并且有一个以上的最大值，则任选一个最大值所对应的方向作为两个栅格单元之间流向；
步骤四三、在地形表面模型DEM中，模拟汇流量：
在地形表面模型DEM的每个网格节点处赋予不同深度的水位，按照水流从高处流向低处的自然规律，依据区域地形和水流方向计算各个栅格单元所有方向的水量，获得汇流多向累积量为：
DU＝D1+D2+D3+...+DN；其中，Dn为某一中心栅格单元第n个流向流来的水流量，n＝1,2,3,...,N,N表示中心栅格单元的水流流向的数量。
Dn＝DnM，其中M为1或0，若流向指向目前需要计算汇流多向累积量的栅格单元，则M＝1，否则M＝0。
步骤四四、根据获得的各个栅格单元的汇流多向累积量，设置不同的汇水面积值以模拟出不同条件下的河网，与真实河网进行比较，找出最接近真实河网的汇水面积值，利用该汇水面积值生成河网。
其中，C计为通过设定不同的汇水面积值而获得的模拟数字河网数据，C真为真实数字河网数据。当Q的值趋近于1表示此汇水面积值下对应的数字河网接近于真实数字河网，即可将此时的汇水面积值作为生成河网的输入数据。
本实施方式步骤一中，致灾因子数据包括连续型数据和离散型数据，连续型数据为高程、坡向、距河流距离、距道路距离、剖面曲率、NDVI，离散型数据为岩土体类型、土地利用类型以及斜坡结构。
具体实施例：为了实现上述的预测和规划目标，提供了一种基于地形单元的滑坡易发性评判方法。通过同时考虑评价单元以及致灾因子，采用随机森林模型作为滑坡易发性评价模型，提供未来滑坡易发性地段的合理预测。下面结合附图和具体实例对本发明采用的计算方案分几个步骤作进一步说明：
1、数据获取：
本实施例使用的数据包括某省某市某县1:1万县级二调数据、1：650000地质灾害分布图，地理空间数据云的GDEMDEM 30m分辨率数字高程数据、GF-1遥感影像以及湖北省气象局雨量站点数据。
2、所选致灾因子分为连续型数据和离散型数据，统一使用离散型数据代入运算。为方便之后统一代入模型，首先运用等间距方法将上述所选致灾因子中的连续型致灾因子进行离散化，并将其离散化数据代入面积分类比公式统计各个致灾因子在各离散区间分布面积占比，所使用面积分类比AR；
3、采用面积分类比曲线与权重比指数曲线交点及拐点作为判断标准将上述离散化区间进一步优化，将滑坡数据以及各致灾因子计算所得的面积分类比值数据代入面积分类比公式；
4、统计经上述面积分类比及权重比指数优化后各离散区间面积及对应滑坡点数，计算滑倾度指数以代入随机森林评价模型，从而计算各致灾因子的致灾倾向。
5、首先依据实施例区相关资料得到此区域相对高差位于500-1300m之间。然后利用在ArcGIS平台中计算得到的实施例区洼地深度分布图与真实地形的洼地分布图进行对比，从而计算得出正地形的洼地深度。用0减去正地形的高程值即可得到负地形，参考正地形洼地填充阈值计算方法，计算出负地形洼地填充阈值。最终将实施例区正地形的洼地充填阈值设置为500，负地形的洼地充填阈值设置为8000。
6、将实施例区GDEMDEM 30m分辨率数字高程数据代入流向均指数模型进行计算
从而确定各栅格单元的水流方向，模拟真实水流方向。
其中，Gi为两个栅格单元中心之间的直线距离，Gz为两个栅格单元之间的高程差：
若Max(FI)>0，则表明此栅格单元流向未定；
若Max(FI)≤0并且最大值只有一个，则将对应此最大值的方向值作为中心栅格单元处的流向值；
若Max(FI)＝0并且有一个以上的0值，则以这些0值所对应的方向值相加的结果作为流向值。
7、为模拟真实汇流量，采用汇流多向累积量并利用上述所得水流方向统计各个栅格单元所有方向的水量，即将每个中心单元邻近栅格对应的水流量输入汇流多向累积量；
基于以上所得的汇流多向累积量，设置不同的汇水面积值以模拟出不同条件下的河网，与真实河网进行比较，当汇水面积值判断公式Q的值趋近于1表示此汇水面积值下对应的数字河网接近于真实数字河网，即可将此时的汇水面积值作为生成河网的输入数据。
本实施例首先基于ArcGIS平台并结合上述操作得到充填好的无洼DEM，分别取正地形汇水面积阈值为500、1000、1500…12000等9个值；取反地形汇水面积阈值500、1000、1500…12000等9个值，计算相应的沟谷密度并与真实地形河网密度进行比较，其具体计算结果如图2所示。最终本实例将正地形的汇水面积阈值设置为12000。同理，选取负地形的汇水面积阈值为8000。将上述计算得到的正负地形模拟汇水面积值输入ArcGIS中对应模块分别提取正负地形的汇水流域，将二者合并即可得到对应于实施例区的地形单元。
以上述得到的地形单元作为滑坡易发性评价的基本单元，将每种致灾因子各个优化后区间的敏感性值赋给每个地形单元，带入构建的随机森林模型，并利用ROC曲线下面积AUC值对模型分类效果进行评价。AUC值介于0.5到1.0之间，具体表现为AUC值越大，分类效果越好。本实施例模型AUC值为0.82(图2)，该值满足空间分布验证精度要求。
本实施例基于县内的9个滑坡致灾因子，采用随机森林模型经过反复调整关键模型参数对实施例区进行滑坡易发性评价，从而确定实施例区每个栅格的易发性概率值，据此利用自然间断法绘制出实施例区滑坡易发性分布图，其结果如表1、图4所示。结果表明，实施例区易发性等级为中易发以下的区域占77.84％，滑坡发生在极高和较高易发区的占55.92％。较高易发区和高易发区大多分布于长江及其支流两岸所在的沙镇溪镇、郭家坝镇、香溪镇、泄滩乡等区域。经实地考察及历史滑坡资料进行比对可知，基于地形单元的滑坡易发性评价结果是合理可靠的，与秭归县历史滑坡分布特点吻合程度较高。
表1实施例区不同易发性等级对应的滑坡灾害点数及分级单元数
