1. 引言
邛崃山系位于我国第一阶梯和第二阶梯的分界线上,地形地貌复杂,地质构造活跃[1]。该区域内分布着大量的滑坡灾害,加之地形、降水和人类活动等因素的耦合作用,研究区滑坡灾害将长期呈现多发、易发、高发态势,频繁发生的滑坡灾害不断造成生命财产损失、严重制约着山区基础设施建设的发展。
随机森林模型在广域灾害易发性评价中表现出了很好的预测精度,特别是在处理大数据量、高维度数据方面具有很强的非线性处理和泛化能力。然而,现有机器学习模型存在一些挑战。例如,难以有效选择非滑坡样本,从而降低了学习效率;模型的可解释性不足;传统的稀疏性分级方式可能导致信息丢失;超参数设置不合理可能降低准确率。因此,未来的机器学习模型需要解决这些问题,以提高滑坡易发性评价的准确性和可靠性。
针对上述机器学习模型存在的问题,本文利用随机森林模型对邛崃山系进行滑坡易发性评价。利用信息量模型进行滑坡易发性预评价,在低易发区提取与滑坡灾害点等量的非灾害点。在此基础上对邛崃山系进行滑坡易发性评价,并根据评价结果进行滑坡易发性分区,以期提高随机森林模型的预测能力和评价结果精度。
2. 研究区域概况
邛崃山系位于我国四川省中部,青藏高原向四川盆地过渡的山地地区,第一阶梯与第二阶梯的地理分界线和农业分界线上,同时也是大熊猫国家公园的第二大山系。邛崃山系地理位置如图1所示,该地区涵盖阿坝藏族羌族自治州、甘孜藏族自治州、雅安市和成都市的部分区域,其中包括宝兴县、汶川县、天全县、芦山县、崇州市、大邑县、理县、都江堰市、康定市、小金县和泸定县11个县市的部分区域。整个区域南北长约250 km,呈北宽南窄状。
Figure 1. Map of the geographical location of the Qionglai mountain range
图1. 邛崃山系地理位置图
邛崃山系位于中国最显著的地势陡变地带,主要由高山峡谷地貌类型构成。其地貌景观呈现河谷狭窄,谷壁陡峭的特点,岭谷高差达2000~3500 m [2]。如图1所示,该山系呈南北走向,地势从西北往东南递减。研究区域内的山峰海拔高度主要分布在3500~6000 m之间,其中四姑娘山、霸王山、巴朗山等海拔均超过5000 m,海拔5000 m以上的山峰积雪常年不化,有现代冰川分布和古冰川遗迹。研究区域海拔落差可达4000 m以上,地势起伏较大,沟谷深切,地形复杂多变。同时,研究区域由于受流水和冰川等外力作用的影响,地表受到强烈的侵蚀作用使其呈现出山河相间分布、纵列分布的地貌特征,该区域地形地貌为滑坡灾害的发育提供了有利的发育环境。
3. 研究方法
3.1. 信息量模型
信息量法是一种量化信息的价值和传输的统计方法,它的核心思想是通过度量信息的数量和质量,来评估事件的难易程度和不确定性。在地质灾害预测中,信息量法可以用于评估各个因子对地质灾害的影响,并计算它们对地质灾害易发性的贡献[3]。在信息量法中,通过计算信息熵来描述数据的不确定性和随机性。信息熵越高,数据的随机性和不确定性就越大,而信息熵越低,数据的确定性就越高。此外,可以利用信息量法来计算每个评价因子的信息量值,这个值可以用来表示每个评价因子对地质灾害发生的难易程度,信息量值越大,表示这个因素对地质灾害的影响越大。
在地质灾害易发性评价中,利用信息量法提取各个因子在地质灾害中提供的信息量,将各个因子对地质灾害的贡献量化。X表示成灾因子(i表示n个因子中的一个类别,j表示m个因子分级中的一个分级区间),则Xij对B点的致灾信息量表示为:
(1)
式中:P (B/Xij)表示在Xij区间发生地质灾害的概率,P (B)表示地质灾害发生的概率。在实际计算中很难直接计算某个区间的因子发生地质灾害的概率P (B/Xij),因此为方便计算采用样本的频率值来替代概率值,如下式所示:
(2)
式中:Nij/Sij表示在研究区内在因子所属区域发生地质灾害的数量与因子分布面积的比值。N/S表示研究区发生地质灾害数量与研究区面积的比值。
信息量法在地质灾害评价中可以将研究区域分成各个评价单元进行分析,每个评价单元可以单独计算其总信息量,并通过信息量大小进行比较,从而对可能发生地质灾害的区域进行有效评估。此外,信息量模型还可以准确地表达各种地质灾害影响因素的影响程度,这对于评价地质灾害易发性和采取相应的防控措施具有重要意义。因此,信息量模型是一种具有应用前景的地质灾害评价方法。
3.2. 随机森林模型
评价模型的选择直接决定评价结果的精确度。在地质灾害易发性评估模型中,随机森林模型(Random Forest, RF)基于多棵决策树(Decision Tree, DT)的预测优势,既能处理大样本数据,又能在二分类问题上表现突出,在大范围研究区域情况下表现良好[4]-[6],因此本文选取随机森林模型开展地质灾害易发性评价。
(1) 决策树
决策树是一种基于非参数化和非线性的监督学习模型。它通过对独立特征进行优化,将输入数据集分割成不同的子集,并支持类别型和连续型的输入和输出变量。在每个节点处做出二元决策,创建一个分裂来将一个或多个类别从全局数据集中分离出来[7]。在节点进行分裂的过程中,可以根据不同的基于规则,采用信息增益(IG)、信息增益率(IGR)或基尼指数(GI)等算法进行分裂,从而构建出不同的决策树模型,如ID3、C5.0和CART等[8]。以C5.0为例,它基于IGR计算最佳的分裂,IGR被认为是一种概率化的度量,用于计算不确定性减少的程度。通常,决策树通过计算具有最大IGR的分裂来生长,直到找到最佳解决方案。IGR的计算公式如下:
(3)
IGR是信息增益比,N是全局数据集,T是预测变量,而G (N, T)是原始节点和新节点之间熵差异,计算方式如下:
(4)
在这里,C是目标变量集合,t是C的类别数,K是T的类别数,即Ci (I = 1, 2, ..., t),Tj (j = 1, 2, ..., k)。
在机器学习中,决策树是一种常见的分类和回归算法。但是,单个决策树模型容易出现过拟合的问题。这是由于决策树的单一决策方式,使得树越深时决策规则越复杂,从而能够“完美”的拟合当前的数据集。但是,这也就使得模型泛化能力较差,不能很好地适应新的数据,预测结果不够准确。因此,为提高决策树模型的泛化能力,常常会采用集成学习的方法,如随机森林、梯度提升树等来组合多个决策树模型,以减小过拟合的风险,提高预测性能[9]。
(2) 随机森林(Random Forest, RF)
随机森林是一种有监督的分类算法,是一种集成方法,使用DT模型,使得每棵树都适合于独立采样的数据子集,使用bootstrapping方法[10]。随机森林模型是由许多单独的决策树模型组成,如图2。RF以高精度率而闻名,可以对预测中的异常值提供高精度率,这是由于在每个分割节点上根据两个数据对象进行随机选择,即Out-Of-Bag (OOB)和相似度,使用OOB数据获取变量重要性估计和内部无偏的OOB误差(即分类误差),当树被添加到森林中时,使用“装袋”(一种自助法)随机选择变量样本作为用于模型校准的训练数据集[11]。对于每个变量,该函数确定如果对该变量的值在OOB观察中进行排列,则该变量的建模预测误差。另一方面,相似性用于替换缺失数据,定位异常值,并且只能在每个树适用于每对案例后计算,并通过除以适用的总树数进行归一化[12]。
Figure 2. Random forest model comprising multiple decision trees
图2. 由多个决策树组成的随机森林模型
4. 滑坡易发性评价体系
4.1. 评价单元选取
在进行易发性评价之前需要确定评价单元,评价结果是在进行易发性评价时需要计算的最小单元。每个评价单元中包含各项评价因子数据,是构成研究区易发性评价的前提。地质灾害易发性评价中常用的评价单元有网格单元、栅格单元、斜坡单元、流域单元、地域单元等。其中网格单元、栅格单元形状较规则,斜坡单元、流域单元、地域单元形状不规则。不同类型的评价单元各自具有独特的优缺点和适用性[12]。
综合考虑评价单元划分方法在邛崃山系的适用性,本文主要选取栅格单元作为研究区进行地质灾害易发性评价的基本单元。栅格单元形状规整,边界清晰,并且在划分时快捷简单,利于ArcGIS进行空间分析和处理。根据经验公式[13]:
(5)
式中,𝐺𝑆为合适的栅格单元大小,S为原始基础数据精度的分母。本文基础数据比例尺为1:250,000,由公式5可得,本文使用30 m × 30 m的栅格单元作为研究区评价单元。
4.2. 评价因子相关性分析
Figure 3. Evaluation factors and spatial distribution of landslide hazards
图3. 评价因子与滑坡灾害空间分布图
本文根据邛崃山系具体情况,以邛崃山系历史滑坡灾害数据,综合考虑研究区域滑坡点的影响因素及相关地质资料,选择高程、坡度、坡向、降雨量、归一化植被指数、岩组、距断层距离、距河流距离、距道路距离9个评价因子构建评价指标体系。
考虑到不同因子栅格表达的物理意义以及量纲的不同,为保证数据在使用时的一致性,首先利用ArcGIS对9个评价因子进行归一化处理,然后将9个因子特征值提取到滑坡点属性表中。当相关系数绝对值小于0.5,可以认为两因子之间相互独立。通过个因子相关性检验,定量确定因子之间的相关性大小。在进行各因子相关系数检验时需要根据不同的数据类型选择合适的相关分析方法,常见的相关性分析方法包括:卡方检测,Spearman系数,Pearson系数和Eta系数,文中因子的数据类型既有标度型又有定序型,所以在选择相关性分析方法时应选择斯皮尔曼相关系数分析方法。斯皮尔曼相关系数通过对输入数据进行排序求解,进而确定各变量之间的相似程度,相似程度与斯皮尔曼相关系数值成正比关系,通常用R定量反映因子之间的相似程度,即相关系数。R的值域是[−1, 1],R > 0时,因子呈现正相关;R < 0时,因子呈现负相关。
Table 1. Spearman correlation coefficient of each factor
表1. 各因子斯皮尔曼相关系数表
因子 |
高程 |
坡度 |
坡向 |
NDVI |
降雨量 |
岩组 |
距断层距离 |
距道路距离 |
距水系距离 |
高程 |
1.00 |
−0.127 |
−0.035 |
−0.183 |
−0.074 |
0.103 |
0.010 |
0.305 |
−0.177 |
坡度 |
|
1.00 |
0.018 |
0.042 |
−0.144 |
−0.172 |
0.007 |
−0.031 |
0.090 |
坡向 |
|
|
1.00 |
0.179 |
−0.173 |
0.029 |
−0.091 |
0.050 |
0.104 |
NDVI |
|
|
|
1.00 |
−0.313 |
−0.007 |
−0.194 |
0.073 |
0.259 |
降雨量 |
|
|
|
|
1.00 |
0.021 |
0.098 |
−0.051 |
−0.091 |
岩组 |
|
|
|
|
|
1.00 |
0.070 |
−0.100 |
0.059 |
距断层距离 |
|
|
|
|
|
|
1.00 |
0.124 |
−0.243 |
距道路距离 |
|
|
|
|
|
|
|
1.00 |
−0.062 |
距水系距离 |
|
|
|
|
|
|
|
|
1.00 |
对评价因子进行相关系数分级,统计结果如表1所示。经过对各因子相关系数的绝对值进行分析,发现它们都小于0.4,说明这些因子之间的相关程度并不高,无需删去任何一个影响因子。因此,本文最终选定高程、坡度、坡向、归一化植被指数、降雨量、岩组、距断层距离、距道路距离和距水系距离这9个因子来构建研究区滑坡易发性评价体系。评价因子与滑坡灾害空间分布如图3所示。
4.3. 评价因子分级
Table 2. Landslide hazard evaluation factor grading scale
表2. 滑坡灾害评价因子分级表
影响因素 |
分级(分类) |
高程(m) |
① 417~1217 ② 1217~2201 ③ 2201~3116 ④ 3116~3962 ⑤ 3962~6250 |
坡度(˚) |
① 0~13.45 ② 13.45~23.02 ③ 23.02~32.46 ④ 32.46~43.14 ⑤ 43.14~77.85 |
坡向(˚) |
① 337.5~22.5 ② 22.5~67.5 ③ 67.5~112.5 ④ 112.5~157.5 ⑤ 157.5~202.5 ⑥ 202.5~247.5 ⑦ 247.5~292.5 ⑧ 292.5~337.5 |
归一化植被指数 |
① 0~0.38 ② 0.38~0.52 ③ 0.52~0.63 ④ 0.63~0.71 ⑤ 0.71~1 |
降雨量(mm) |
① 474~538 ② 538~576 ③ 576~612 ④ 612~648 ⑤ 648~791 |
岩组 |
① 坚硬岩组 ② 较硬 ③ 较硬夹较软岩组 ④ 较软夹较硬岩组 ⑤ 较软岩组 ⑥ 软岩组 ⑦ 土体 |
到断层距离(m) |
① 0~500 ② 500~1000 ③ 1000~1500 ④ 1500~2000 ⑤ 2000~2500 ⑥ 2500~3000 ⑦ >3000 |
距道路距离(m) |
① 0~100 ② 100~200 ③ 200~300 ④ 300~400 ⑤ 400~500 ⑥ 500~600 ⑦ >600 |
距水系距离(m) |
① 0~500 ② 500~1000 ③ 1000~1500 ④ 1500~2000 ⑤ 2000~2500 ⑥ 2500~3000 ⑦ >3000 |
自然间断点法是一种基于数据分布的分级方法,其主要思想是根据数据的分布情况划分连续属性为有限个离散区间,而无需人为指定分割点。具体来说,自然间断点法会寻找数据分布中的“自然间隔”,并将连续属性分级为有限个离散区间。这种方法的优点在于可以更好地反映数据的分布情况,同时能够保留数据的原始信息。但其缺点是可能因数据分布不均匀而导致某些区间内的样本数量过少。本文采用自然间断点法进行评价因子分级,滑坡灾害评价因子分级情况如表2所示。
5. 滑坡易发性分析结果
5.1. 信息量提取非灾害点
非灾害点数据集不能充分代表整个研究区中易发性低的无灾害区域的特征,那么可能会存在采样点不足以代表该区域的特征的问题。此外,如果选择的某些采样点正好位于高易发性级别,这些点记录的是易发生滑坡灾害的特征属性,这就会导致非灾害点数据集整体上呈现一种扰动。为保证非灾害点子数据的准确性,本文采用信息量法提取非灾害点,具体步骤如下:
(1) 首先将九个致灾因子进行重分类;
(2) 根据因子的重分类结果,统计各个因子分级上滑坡点的数量,并计算每个因子的分级信息量如表3所示;
(3) 在ArcGIS中叠加因子的信息量,获取信息量易发性栅格;
(4) 对易发性栅格进行重分类,从低易发性区域提取非灾害点。
Figure 4. Distribution of control datasets in the study area
图4. 研究区对照数据集分布情况图
采用信息量法提取出2373个非灾害点,非灾害点及滑坡点分布情况如图4所示。
Table 3. The magnitude of information for each evaluation factor
表3. 各评价因子信息量值
评价因子 |
分级 |
滑坡数量(处) |
信息量值 |
高程(/m) |
417~1217 |
1134 |
2.93 |
1217~2201 |
780 |
1.89 |
2201~3116 |
359 |
0.75 |
3116~3962 |
79 |
0.14 |
3962~6250 |
21 |
0.04 |
坡度(/˚) |
0~13.45 |
473 |
1.59 |
13.45~23.02 |
763 |
1.59 |
23.02~32.46 |
552 |
0.80 |
32.46~43.14 |
406 |
0.64 |
43.14~77.85 |
179 |
0.66 |
坡向 |
北 |
142 |
0.52 |
东北 |
354 |
1.35 |
东 |
319 |
0.96 |
东南 |
378 |
1.16 |
南 |
366 |
1.22 |
西南 |
460 |
1.66 |
西 |
248 |
0.80 |
西北 |
106 |
0.37 |
NDVI |
0~0.38 |
1 |
0.03 |
0.38~0.52 |
74 |
0.22 |
0.52~0.63 |
470 |
1.20 |
0.63~0.71 |
912 |
1.38 |
0.71~1 |
916 |
0.96 |
降雨量(/mm) |
474~538 |
499 |
0.65 |
538~576 |
518 |
0.58 |
576~612 |
590 |
1.82 |
612~648 |
466 |
1.95 |
648~791 |
285 |
2.14 |
岩组 |
坚硬岩组 |
163 |
0.84 |
较硬岩组 |
128 |
1.14 |
较硬夹较软岩组 |
9 |
0.26 |
较软夹较硬岩组 |
2 |
0.05 |
较软岩组 |
369 |
0.87 |
软岩组 |
1663 |
1.19 |
土体 |
38 |
0.23 |
距断层距离(/m) |
0~500 |
499 |
2.55 |
500~1000 |
317 |
1.83 |
1000~1500 |
238 |
1.64 |
1500~2000 |
205 |
1.63 |
2000~2500 |
145 |
1.33 |
2500~3000 |
102 |
1.07 |
距道路距离(/m) |
>3000 |
868 |
0.57 |
0~100 |
524 |
12.14 |
100~200 |
198 |
5.64 |
200~300 |
155 |
4.13 |
300~400 |
98 |
3.14 |
400~500 |
86 |
2.61 |
500~600 |
70 |
2.22 |
>600 |
1242 |
0.57 |
距水系距离(/m) |
0~500 |
643 |
3.46 |
500~1000 |
448 |
2.54 |
1000~1500 |
215 |
1.27 |
1500~2000 |
161 |
0.98 |
2000~2500 |
136 |
0.86 |
2500~3000 |
118 |
0.77 |
>3000 |
652 |
0.48 |
5.2. 信息量支持下随机森林滑坡易发性评价
对研究区内滑坡点和非灾害点进行分类,将滑坡点赋值为“1”,非灾害点赋值为“0”,将正样本和负样本一起用于随机森林模型训练。利用ArcGIS对9个评价因子的样本点进行提取,将样本中的70% (1661个)作为训练样本,将训练数据集输入预先编写好的代码进行训练,剩余的30% (712个)作为测试数据集,使用训练完成的模型对测试集进行验证。训练得出的随机森林模型ROC曲线如图5所示,该模型的AUC值为0.940。
接下来,将整个研究区域的因子数据输入模型,得出该区域的滑坡灾害易发性指数,并采用ArcGIS软件中的自然断点法,将其分为五个级别,分别为极低易发区、较低易发区、中易发区、较高易发区和极高易发区,最终生成滑坡易发性图,如图6所示。邛崃山系西部的极高易发区和高易发区呈线性分布,该区域沟谷交错,地势起伏大,内外营力作用强烈;东部呈带状分布,该区域属于山地、丘陵、平原建有地貌,水系发达,道路相对研究区西部较密集,断裂带较发育,且沉积了深厚的第四系冰水堆积、洪水堆积的松散沉积物,在雨季时受到冲刷,容易引发地质灾害。
低易发区和极低易发区占研究区总面积59.9%,呈不规则的片状分布。邛崃山系东部低易发区和极低易发区虽然河流、道路分布较密集,但该区域地势较平坦,绝大部分属于平原地貌,且该区域距离断层相对较远,因此,易发性等级较低。即使在邛崃山系海拔较高区域,仍然存在一定数量的低易发区和极低易发区。在高海拔区域,由于蓄水条件差、坡体含水率低等因素,地质灾害的发生可能性相对较低。此外,植被覆盖率较高的区域也存在较多的低易发区,因为植被的根系可以起到稳固土层的作用,从而在一定程度上减少了地质灾害的发生。
Figure 5. ROC curves of the random forest model
图5. 随机森林模型的ROC曲线
Figure 6. Zoning map of landslide susceptibility assessment in the Qionglai Mountains
图6. 邛崃山系滑坡易发性评价区划图
6. 滑坡灾害影响分析
6.1. 滑坡的影响分析
由图7所示,邛崃山系居民点和道路总体分布不均,主要分布于东部和山区河谷地带。东部居民点呈片状分布,交通路线密集,位于成都平原西部,该区居民点和道路分布密集,呈大范围聚集的特点,同时该区域滑坡易发概率低,属于滑坡极低、低易发区域。中部是山区向平原过渡地带,该区域滑坡易发性由较高和极高等级向中等及较低转变,居民点和道路分布最少。研究区西部和北部山区河谷地带居民点沿河流呈线状分布,主要集中在理县、小金县、康定市东北部、宝新县南部、泸定县东部、荥经县北部及石棉县东北部、汶川县东部河谷地带。道路、居民点主要位于滑坡高易发区和中等易发区,居民点呈小范围聚集的特点,道路呈现沿河谷分布的特点。
Figure 7. Landslide susceptibility and geographical element distribution map
图7. 滑坡易发性与地理要素分布图
居民点与滑坡易发等级的分布情况直接影响当地居民生命财产安全,如损毁房屋建筑甚至威胁到居民的生命安全等;道路与滑坡易发等级的分布情况影响着当地基础设施建设安全等问题,如滑坡会阻塞并破坏道路。因此,在考虑滑坡防治措施时,不仅需考虑滑坡本身,还需结合居民点、道路等要素的地理分布情况。
6.2. 邛崃山系滑坡灾害防治分区
滑坡防治分区与易发性分区有所不同,其考虑因素不仅包括滑坡发生的可能性,还关注对人类的威胁程度。具体而言,本研究将人口居住密度大有道路经过的极高易发区和高易发区作为重点防治区域;将居住密度较低的部分高易发区和中等易发区以及人口密度较大的低易发区作为次重点防治区域;而人口稀缺的地区以及滑坡极低易发区则作为一般防治区域。本文通过滑坡防治分区方式,可以更加有效地采取针对性的防治措施,确保人类在面临滑坡灾害时的安全。
综合考虑滑坡易发性分区、居民点、交通等因素得到的邛崃山系滑坡防治结果如图8所示。
Figure 8. Result map of landslide prevention and control zoning
图8. 滑坡防治分区结果图
图中明显可以看出重点防治区的居名点密度高,大部分居民点位于该区域。其中研究区最东部居民点同样很高,但是由于该区域位于极低易发区,因此被划分为一般防治区域。
7. 结论
以邛崃山系作为研究对象,通过分析滑坡空间分布规律与影响因子之间的关系,对影响因子进行相关性分析,以高程、坡度、坡向、归一化植被指数、降雨量、岩组、据断层距离、据道路距离和距水系距离这9个因子建立邛崃山系滑坡易发性评价指标体系。利用信息量支持下的随机森林模型开展滑坡易发性评价。
(1) 本文运用信息量支持下的随机森林算法,在该地区开展滑坡易发性评价研究。结果显示,滑坡极高易发区面积占研究区面积的7.1%,滑坡点密度为51.1个/100 km2。通过ROC曲线验证随机森林模型,ROC曲线的AUC值为0.940,说明信息量支持下随机森林模型具有较高的可靠性和准确性,本文滑坡易发性结果较理想。
(2) 研究区滑坡极高易发区和高易发区沿道路、河谷以及断层分布。邛崃山系西部的极高易发区和高易发区呈线性分布,该区域沟谷交错,地势起伏大,内外营力作用强烈;东部呈带状分布,该区域属于山地、丘陵、平原建有地貌,水系发达,道路相对研究区西部较密集,断裂带较发育,且沉积了深厚的第四系冰水堆积、洪水堆积的松散沉积物,在雨季时受到冲刷,容易引发滑坡。
基金项目
铁路病害智能检测技术与成套装备研发(KSNQ233011)。
NOTES
*第一作者。