山东科技大学测绘与空间信息学院刘雅婷：顾及空间异质性和特征优选的滑坡易发性评价方法|《测绘学报》2024年53卷第7期

本文内容来源于《测绘学报》2024年第7期（审图号GS京(2024)1261号）

刘雅婷,1,2, 陈传法,1,2

1.山东科技大学测绘与空间信息学院，山东 青岛

2.

山东省基础地理信息与数字化技术重点实验室，山东 青岛

基金项目

国家自然科学基金
(


);

山东省自然科学基金
(

ZR2020YQ26
);

山东省高等学校青创科技支持计划
(

2019KJH007）

作者简介

第一作者：刘雅婷（1998—），女，硕士生，研究方向为地质灾害早期识别和风险分析。E-mail：

chlyt2017@163.com


通讯作者: 陈传法 E-mail：chencf@sdust.edu.cn
摘要
高效、精准、可靠的滑坡易发性评价方法是灾前科学预警和全面防治的关键手段。然而，传统滑坡易发性评价方法未能有效解决空间异质性和冗余特征造成的预测偏差。针对该问题，本文提出了一种顾及空间异质性和特征优选的滑坡易发性评价方法（spatial feature optimized Stacking，SF-Stacking）。该方法首先使用AGNES聚类（agglomerative nesting）将研究区全局栅格单元分成若干个局部子区，然后采用一种特征优选策略为每个子区选择最优致灾因子组合，最后采用Stacking集成技术耦合多种机器学习算法实现滑坡易发性评价。以宜宾市为研究区，基于滑坡灾害易发性分区图和统计学指标，将SF-Stacking方法与7种传统方法对比表明，SF-Stacking方法的准确性最优，稳健性最强，可解释性最高。
关键词


滑坡易发性
;

空间异质性
;

特征优选
;

Stacking集成学习

本文引用格式

刘雅婷, 陈传法.
顾及空间异质性和特征优选的滑坡易发性评价方法
[J]. 测绘学报, 2024, 53(7): 1417-1428 doi:10.11947/j.AGCS.2024.

LIU Yating, CHEN Chuanfa.
Landslide susceptibility evaluation method considering spatial heterogeneity and feature selection
[J].

Acta Geodaetica et Cartographica Sinica
, 2024, 53(7): 1417-1428 doi:10.11947/j.AGCS.2024.

阅读全文

http://xb.chinasmp.com/article/2024/1001-1595/1001-1595-2024-07-1417.shtml

我国地质环境复杂，地质灾害呈高发态势，防治形势严峻[1]。在众多地质灾害中，滑坡灾害发生最频繁，影响范围最广，严重威胁着人民生命财产安全、生态环境安全及社会经济建设[2]。然而，山体滑坡具有隐蔽性、突发性和复杂性等特点，并且不同地质环境孕育的滑坡成灾机理各异，预警难度非常大[3]。因此，如何实现区域精细化滑坡易发性准确评价，有效提高滑坡灾害预报预警能力，一直是滑坡灾害防治领域的研究热点和难点[4]。

滑坡易发性评价是通过分析滑坡特征因子与滑坡空间分布的潜在关系来预测区域内滑坡发生概率。滑坡灾害的发育过程复杂，在内部孕灾因子（如地形地貌、地层岩性、地质构造、交通水系等）与外部诱发因子（如降雨、地震、人类工程活动等）的作用下，不同灾害环境发生滑坡的空间概率各不相同[5]。当前常用的滑坡易发性评价方法主要包括数理统计法[6-7]和机器学习法[8-9]。数理统计法旨在构建明确的数学表达式以预测滑坡发生概率。该类方法简单易行，可解释性高，却难以准确表达特征因子与滑坡之间的复杂关系，预测性能不高[10]。区别于数理统计法，机器学习法不需要在特征因子与滑坡间建立可解释的数学公式，具有较好的预测性能，但当训练数据不足时容易出现欠拟合或过拟合现象[11]。因此，为提高机器学习模型的泛化能力，集成学习方法被广泛采用[12]。然而，上述研究均采用全局模型对研究区滑坡进行整体预测，忽略了不同致灾因子作用差异带来的空间异质性，由此导致滑坡预测偏差。为了解决空间异质性问题，文献[13—14]提出了基于空间分区的局部栅格单元法，但是未能考虑冗余特征的影响。此外，文献[15—16]采用特征选择方法筛选出相对重要的特征因子，在一定程度上增加了模型的可靠性，然而仅使用单一特征选择方法难以平衡特征子集最小化和性能最大化，导致精度提升有限[17]。

鉴于此，本文提出一种顾及空间异质性和特征优选的滑坡易发性评价方法（spatial feature optimized Stacking，SF-Stacking）。该方法的创新性主要包括：①采用聚类分析将研究区分成若干个均质子区，消除空间特征异质性造成的预测偏差；②耦合多种特征选择方法筛选各个子区最优空间特征，提升特征优选的可靠性；③使用集成学习算法拟合重要特征与滑坡之间的复杂非线性关系，提高模型的普适性。

1 研究区概况及数据源

1.1 研究区概况

宜宾市（图1）位于四川省南部（103°20’E—103°36’E，27°50’N—29°16’N），总面积约13 298 km2，地处川、滇、黔三省交界地带，岷江、金沙江、长江三江交汇口。该区域东北部与四川盆地盆底丘陵相接，西北侧属盆中方山丘陵区，西部为大小凉山余脉，南部为云贵高原北坡，地势西高东低，地形复杂，地貌多样，最高海拔达1995 m。宜宾市属于亚热带湿润季风气候，年平均降雨量为844~1068 mm，每年7~9月是暴雨频发的季节，加之地震、边坡开挖和修建道路等人类工程活动，使得该市滑坡灾害频发[18]。

图1

图1 研究区滑坡灾害点分布

Fig.1 Distribution of landslide disaster points in the study area

1.2 数据源

用于滑坡易发性评价的数据源包括：①滑坡灾害数据来源于四川省自然资源厅，共686例历史滑坡数据（2008—2018年）（图1）；②年平均降雨量、岩性数据来源于中国科学院资源环境科学数据中心；③DEM数据来源于美国NASA地球数据网站的30 m SRTM，主要用于高程、坡度、坡向和地形湿度指数（terrain wetness index，TWI）等地形地貌基础环境因子的获取；④30 m空间分辨率的Landsat-8遥感影像数据来源于美国地质勘探局，用于归一化植被指数（normalized difference vegetation index，NDVI）的提取；⑤道路数据、水系数据、兴趣点（point of interest，POI）密度数据是由OpenStreetMap提供的公开数据集；⑥地震动峰值加速度（peak ground acceleration，PGA）来源于中国地震局地球物理研究所）；⑦断层数据来源于地震活动断层探察数据中心；⑧土地利用数据来源于欧洲航空局制作的10 m空间分辨率World Cover数据。

1.3 特征因子

根据研究区滑坡致灾机理，本文从地形地貌、水文地质、人类工程活动、降雨等方面选取了高程、坡度、坡向、平面曲率、剖面曲率、NDVI、TWI、年降雨量、断层距离、河流距离、道路距离、POI、PGA、岩性及土地利用15个特征因子。为保证空间一致性，将所有因子的空间分辨率均重采样至30 m，借助自然间断法对连续型特征因子进行分级[19]（图2）。

图2

图2 特征因子

Fig.2 Feature factors

2 研究方法

为了降低滑坡易发性建模的不确定性，本文构建了一种顾及空间异质性和特征优选的滑坡易发性评价方法（图3），主要包括数据预处理、模型构建和模型精度评价3部分。具体而言，数据预处理主要是根据研究区地理概况选取基础特征因子，并利用信息量法生成非滑坡点，与滑坡点共同组成样本数据集。在模型构建中，首先利用AGNES聚类方法将研究区划分成不同子区，然后在各子区中采用特征优选方法（VIF-Boruta-GeoDetector，VBG）剔除冗余特征因子，最后基于各子区最优特征因子构建相应的Stacking集成模型。模型精度评价主要借助统计学指标和易发性结果图来评估模型性能。

图3

图3 本文方法技术流程

Fig.3 Flowchart of the proposed method

2.1 数据预处理

为降低传统滑坡易发性评价方法随机选取非滑坡点带来的不确定性，本文利用信息量法选取各区非滑坡点。信息量法是一种统计评价方法，广泛应用于地质灾害预测中的非滑坡点选取[20]。数据预处理流程为：首先，通过式（1）计算上述连续型特征因子各分级状态下的滑坡易发性，进而求和得到研究区全部栅格的易发性；然后，在滑坡极低易发性区随机抽取与滑坡灾害点相同数目的单元作为非滑坡点；最后，合并滑坡点与非滑坡点作为样本数据集，并将原始样本集按照7∶3的比例划分为训练集和验证集。其中，训练集用于训练模型的参数，验证集用于模型精度检验

(1)

式中，x为致灾因子；y为滑坡灾害事件；I为因子x在各分级状态下提供的滑坡灾害信息量；S为研究区栅格单元总数；Sx为研究区x各分级状态的单元数；N为研究区滑坡总数量；Nx为x各分级状态含有滑坡灾害点的栅格数量。当I>0时，表示有利于滑坡灾害发生；当I<0时，表示不利于滑坡灾害发生。

2.2 模型构建

2.2.1 特征空间异质性分区

由于滑坡孕灾机理的复杂性，滑坡特征对滑坡内部构造的影响会随着地理环境不同而存在差异。因此，在构建滑坡易发性评价模型时，需先解决滑坡特征因子空间异质性问题[13]。为此，本文采用AGNES[21]聚类滑坡点致灾因子来弱化空间异质性的影响。具体而言，首先根据CH指数（Calinski-Harabasz）[22]确定簇间距离度量准则及最优聚类簇个数，然后利用AGNES对滑坡点进行聚类，最后对聚类簇结果构建泰森多边形以划分各个子区。其中，簇间距离的度量准则可分为3种：最小距离法（dmin）、最大距离法（dmax）及平均距离法（davg）。

2.2.2 特征选择

不同子区的滑坡发育机理及诱发因素存在明显差异，若将所有特征无差别地输入机器学习模型中会降低数据的有效性及影响模型准确性[23]。因此，训练模型构建前必须进行特征选择。以往研究表明，在进行空间特征选择时，为保证特征间相互独立以及提高模型稳定性，需要利用方差膨胀系数（variance inflation factor，VIF）先对特征因子进行多重共线性检验[24]。不同于VIF，Boruta算法采用基于随机森林分类算法（random forest，RF）的思路对特征选择，可以有效识别重要特征，并且在特征选择的过程中使用随机化策略，避免了因数据变化导致的不稳定性问题。相较于基于机器学习的特征选择方法，地理探测器（GeoDetector）通过对样本空间的划分和对局部线性模型的拟合来挖掘空间异质性，能够筛选出最具有区分性和解释性的特征[25]。鉴于此，本文提出了基于VIF-Boruta-GeoDetector组合（VBG）的特征优选方法。该方法具体流程如下。

步骤1 计算特征因子之间的VIF值。若VIF大于阈值，删除该特征，否则保留。计算公式为

(2)

式中，表示将第i个特征因子当作因变量，其余因子当作自变量进行线性回归分析所得到的决定系数。VIFi值越大表明该特征共线性越强，通常情况下，当VIFi小于10时，表明特征间多重共线性较低。

步骤2 采用Boruta算法对步骤1中的保留特征进行重要性排序，进一步保留被标记为“重要”的特征。

步骤3 对步骤2中保留下来的特征，利用地理探测器中的因子探测器筛选显著性较高的特征，从而形成最优特征子集。其计算公式为

(3)

式中，h=1，2，…，L为评价因子X的分层；Nh和N分别为层h和全区的单元数；σh和σ分别是层h和全区的Y值的方差；q为因子对滑坡易发性空间分异的影响力；取值范围为[0，1]。q值越大，表示因子X对滑坡易发性Y空间分异的影响越大；反之则越小。地理探测器使用P值检验结果的显著性，P<0.05表明结果具有显著性。

2.2.3 模型集成

为了克服单个机器学习算法性能瓶颈，本文采用Stacking集成策略耦合多个机器学习算法的优势来提高模型可靠性及降低过拟合风险。Stacking集成学习算法由基学习器和元学习器组成，通过组合异质学习器构成多层次强学习器，以提升模型的泛化能力[12]。其中，为保证模型多元化和预测性能最大化，本文构建的Stacking集成学习策略采用支持向量机（SVM）、BP神经网络（BPNN）、K最近邻算法（KNN）作为基学习器，逻辑回归（LR）作为元学习器。具体而言，首先利用基学习器对训练集进行k折交叉训练（如k=5），并堆叠每个基学习器每折的输出结果作为新的训练集；然后通过每一次训练对验证集进行预测，将k次预测平均值作为新的验证集；最后将新的训练集和验证集输入到元学习器中得到最终的预测值。

2.3 模型精度评价

为了验证本文方法的有效性，与7种方法进行滑坡易发性评价对比，包括顾及空间异质性的Stacking方法（spatial Stacking，S-Stacking）、顾及特征优选的Stacking方法（feature optimized Stacking，F-Stacking）、全局Stacking集成学习方法（Stacking）和全局机器学习方法（SVM、BPNN、KNN、LR）。其中，全局Stacking集成学习方法和全局机器学习方法均未考虑空间异质性和冗余特征的影响。本文选用准确率（ACC）、F1值、ROC曲线、Kappa系数和Sridevi Jadi试验概率精度（SJ）5个指标对模型性能定量评估。其中，各指标取值范围均为0-1，其值越接近1，表示预测准确率越高[26]。各精度指标公式为

(4)

(5)

(6)

(7)

(8)

式中，TP为滑坡点的正确分类数；TN为非滑坡点的正确分类数；FP为非滑坡点分为滑坡点的数量；FN为滑坡点分为非滑坡点的数量；n表示滑坡点和非滑坡点的栅格数；N为研究区总栅格单元数；S为存在滑坡单元总数；M为极高易发性和高易发性的单元总数；MS为极高易发性和高易发性中存在滑坡的单元总数。

3 结果与分析

3.1 空间分区

为获取最优聚类数目，本文测试了不同聚类数的CH值。图4展示了不同簇间距离度量准则下各聚类数的CH值。可以看出，簇间距离为dmin时，聚类效果最优，最佳聚类数为3，对应的CH值为75；簇间距离为dmax时，聚类效果次之，最佳聚类数为3，对应的CH值为56；簇间距离为dmax时，聚类效果最差，最佳聚类数为11，对应的CH值为19。综上，本文选择dmin法确定滑坡聚类数（即为3），并使用AGNES算法对滑坡点进行分簇。

图4

图4 不同聚类数目下的CH值

Fig.4 CH values under different number of clusters

基于滑坡点聚类结果构建泰森多边形，以此将研究区划分为3个子区（图5）。其中，Ⅰ区有137个滑坡点，分布于宜宾市西北部，该区属于山区，地势陡峭，岩体多为松散岩石，降雨量较大，采矿和建设公路等人类工程活动较多[27]，易发生滑坡灾害。Ⅱ区有239个滑坡点，分布于宜宾市东北部，该区以丘陵为主，地势相对平缓，岩性较为坚硬，存在许多褶皱带和断层，在长期降雨作用下易发生滑坡灾害。Ⅲ区有310个滑坡点，分布于宜宾市南部，该区位于青藏高原和盆地之间的过渡地带，地形起伏较大，地质构造复杂，常有大量的松散土层或岩石碎屑堆积，易受降雨等因素的影响而发生滑坡（图2）。

图5

图5 聚类分区情况

Fig.5 Cluster partition situation

3.2 特征优选

3.2.1 各区特征优选结果

图6是各子区滑坡特征因子多重共线性检验结果。可以看出，Ⅰ区中VIF值最高的特征为高程（4.1），最小为坡向（1.1）；Ⅱ区中VIF值最高的特征为坡度（2.8），最小为坡向（1.1）；Ⅲ区中VIF值最高的特征为坡度（2.4），最小为坡向（1.1）。不同子区下各初始特征因子的VIF值均小于5，说明各特征因子间共线性程度低，具有较好的独立性。

图6

图6 各子区特征多重共线性检验结果

Fig.6 Multi-collinearity test results of characteristic factors in each subregion

图7是各子区特征因子的Boruta算法分析结果。结果表明，不同特征因子对不同子区滑坡的贡献度不同。其中，对于Ⅰ区滑坡具有显著影响的特征为高程，重要性达25%，Ⅱ区为坡度，重要性达33%，Ⅲ区为坡度，重要性达38%。经Boruta算法特征筛选后，Ⅰ区保留高程、坡度、道路距离、岩性、年降雨量、平面曲率、POI、PGA和剖面曲率9个特征；Ⅱ区保留坡度、平面曲率、TWI、剖面曲率、高程、断层距离和土地利用7个特征；Ⅲ区保留坡度、高程、年降雨量、平面曲率、TWI、剖面曲率、断层距离、坡向和岩性9个特征。

图7

图7 各子区特征Boruta算法分析结果

Fig.7 Analysis results of Boruta algorithm for characteristic factors of each subregion

为了得到最小优选特征子集及提高模型解释力，进一步利用GeoDetector来度量各分区经Boruta算法筛选后保留的特征对滑坡空间分布的影响力（图8）。由图8可知，不同子区特征因子的q值有所不同。对于Ⅰ区的滑坡而言，高程的影响最为显著，其解释力高达55%。高程可以反映不同海拔地区的气候、岩石风化程度对滑坡发育的影响，Ⅰ区内滑坡大多分布于高程1000 m以下的区域，并且随着高程增高，滑坡数量逐渐减少。对于Ⅱ区和Ⅲ区的滑坡而言，坡度的影响最为显著，其解释力分别为53%和40%。坡度影响着堆积物的形成与斜坡稳定性，Ⅱ区、Ⅲ区内滑坡集中分布于坡度为5°~30°的区域，较小和较大的坡度不利于滑坡发育。此外，Ⅰ区中的剖面曲率、Ⅲ区中的岩性特征均未通过5%的显著性检验，且解释力较弱，因此在Boruta算法保留的特征基础上进一步剔除了Ⅰ区中的剖面曲率、Ⅲ区中的岩性。综上所述，经VBG方法筛选后，Ⅰ区的最优特征子集为高程、岩性、年降雨量、道路距离、POI、坡度、平面曲率和PGA；Ⅱ区为坡度、平面曲率、TWI、剖面曲率、高程、断层距离和土地利用；Ⅲ区为坡度、高程、年降雨量、TWI、剖面曲率、平面曲率、断层距离和坡向。

图8

图8 各子区特征探测结果

Fig.8 Detection results of characteristic factors in each subregion

3.2.2 特征因子组合分析

为验证VBG特征优选方法的有效性，本文测试了各子区不同特征因子数量的AUC值。首先将15个特征因子按照GeoDetector计算的q值排序，然后逐一剔除q值排序低的特征因子，最后计算不同特征因子组合的AUC值，结果如图9所示。由图9可知，各子区不同特征因子组合下的AUC值不同；当特征因子数量小于6时，AUC值相对较低，说明特征因子较少时难以准确描述滑坡灾害的空间分布。当特征因子组合数量分别为8、7、8时3个子区的AUC值最高，分别为0.992、0.910、0.951，这与VBG特征优选后的各子区特征因子组合相同。当特征因子组合数量进一步增加时，AUC值有降低趋势，说明本文特征优选方法在保证特征因子数量最少的同时还保持预测性能最优。

图9

图9 不同特征因子组合的AUC值

Fig.9 AUC values of different feature factor combinations

3.3 模型精度评价

SF-Stacking方法与7种对比方法的精度如图10所示。结果表明，SF-Stacking的滑坡易发性预测精度明显优于其他方法。从各精度指标来看，SF-Stacking的ACC（0.913）、F1值（0.909）、Kappa系数（0.825）、SJ（0.707）和AUC（0.936）最高，S-Stacking次之，LR最低，说明本文方法准确性最高，稳健性最强。此外，4种全局机器学习模型中，SVM的精度最高。相较于SVM，Stacking的ACC值提升了7%，F1值提升了9%，Kappa系数提升了3%，SJ提升了6%，AUC提升了1.9%；相较于Stacking，S-Stacking和F-Stacking的ACC分别提升了2.5%和1.5%，F1值分别提升了2.1%和1.4%，Kappa系数分别提升了4.8%和2.9%，SJ分别提升了3.5%和3.2%，AUC分别提升了1.7%和1.2%。因此，集成学习技术可以改善单个机器学习模型性能，而且顾及空间异质性和特征优选的策略可以有效提高滑坡易发性精度。

图10

图10 精度评价指标

Fig.10 Precision index

3.4 易发性分区结果

利用自然间断法分别将各方法计算的滑坡易发性指数划分为5个等级，即极低易发区、低易发区、中易发区、高易发区和极高易发区（图11）。由图11可知，8种方法得到的易发性结果图具有相似的空间分布：极高易发区主要位于研究区西北部的丘陵和南部的高原，这是由于该区域地形起伏大，地势较为陡峭，降雨量较大，加上日益增多的人类工程活动影响（如煤矿开采）[18,27]，众多因素共同加剧了该区域滑坡灾害风险（图2）。同时，不同方法下的分区图也存在些许差异。例如，BPNN的高易发性区分布最多，SVM、KNN和LR的中易发区分布最多，SF-Stacking、S-Stacking、F-Stacking、Stacking的各易发区相对均衡，且SF-Stacking的低易发区分布较多。此外，由滑坡机理可知，水系流域两岸的岩体长期受动水压力和水流侵蚀的影响，易发生滑坡灾害，而流域内部不存在发生滑坡的可能。结合图11中的典型区域（图11中的放大窗口）可知，SF-Stacking方法的预测结果能够将流域两岸的滑坡易发性划分为极高，而流域内部的易发性划分为极低，而其他7种对比方法的流域内部均有中、高易发性的存在。相较之下，SF-Stacking的预测结果更符合实际滑坡分布状况。

图11

图11 各方法构建的研究区易发性区划

Fig.11 Vulnerability zoning map of different methods

为进一步检验各方法获取的易发性分区图的准确性，本文通过对比2019—2023年多时相谷歌卫星影像（米级）的光学特征，借助目视解译和实地核查手段识别出了23个新的滑坡区域（图11），并将其与易发性分区图叠加分析。其中，23个外部验证点全部落入SF-Stacking的极高易发区，有18个外部验证点落入S-Stacking和F-Stacking的极高易发区，5个外部验证点落入高易发区，其余方法的极高易发区内分别落入了10~20个外部验证点，但中和低易发区内分别落入了2~6个外部验证点。由此可见，相比于其他7种方法，SF-Stacking方法表现出更高的预测性能。

4 结论

为提高滑坡灾害预测能力，本文提出了一种顾及空间异质性和特征优选的滑坡易发性评价方法。该方法首先使用AGNES将滑坡点聚类，解决滑坡致灾因子空间异质性问题；然后采用VBG对各个子区域特征优选，避免冗余特征对预测模型干扰；最后将4种机器学习方法（SVM、BPNN、KNN和LR）通过Stacking集成技术构建局部易发性评价模型，实现高精度滑坡预测。借助统计学指标和易发性分区图将SF-Stacking方法与7种传统方法对比分析表明：①空间异质性影响滑坡易发性预测精度；②本文提出的特征优选策略提高了模型计算效率和准确率；③SF-Stacking方法预测精度明显优于其他7种方法的预测精度，具有较高的稳健性、可靠性。

综上所述，本文方法有效解决了空间异质性和冗余特征造成的滑坡预测偏差缺陷，实现了滑坡易发性高精度制图，进而可为自然灾害风险评估、预警和防治等国家重大工程提供理论方法和技术支撑。但本文方法也存在一定的局限性，例如该方法需要研究区有足够的滑坡灾害点，且滑坡点与非滑坡点的比例不同可能会造成滑坡易发性评价的不确定性。因此，后续研究中本文将基于先验知识[28]扩大训练样本数量，并寻求滑坡点与非滑坡点间的最优比例来增强SF-Stacking模型的泛化能力。

初审：张 琳

复审：宋启凡

终审：金 君