怒江州滑坡灾害易发性区划研究

简小婷，赵康，左小清，朱琪，朱文

(1.云南省地矿测绘院有限公司，云南 昆明 650218； 2.云南省自然资源厅，云南 昆明 650224；3.昆明理工大学国土资源工程学院，云南 昆明 650093； 4.中国电建集团昆明勘测设计研究院有限公司，云南 昆明 650051；5.云南省基础地理信息中心，云南 昆明 650034)

1 引 言

怒江傈僳族自治州(简称“怒江州”)位于云南省西北部，山高坡陡、峡谷深切、地质环境脆弱、气候复杂，属于滑坡灾害高易发区。近年来，由于自然环境变化和人类工程活动影响，怒江州滑坡灾害频发。迄今为止已爆发了“上帕镇6·30”“福贡县4·10”“兰坪县7·31”等大型滑坡灾害事件，对当地造成了重大人员伤亡和财产经济损失。因此，对滑坡灾害易发性的分析能为怒江州滑坡灾害预测及防治提供可靠的依据，对防灾减灾工作具有重要意义。

滑坡一般可分为单体滑坡和区域滑坡[1]。单体滑坡主要评估灾害个体属性和周围环境因子对滑坡灾害造成的影响[2～4]；区域滑坡则是根据区域地质环境背景、成灾诱发因素以及人类活动状况等对区域滑坡灾害做出评价[5～7]。由于单体滑坡不是区域性地质灾害防治规划重点，也不能支撑从宏观上分析地质灾害的分布规律[8]。因此，本文从区域滑坡视角综合评判怒江州滑坡灾害易发性。滑坡灾害易发性区划描述了滑坡发生概率的空间分布情况，是支撑滑坡防灾减灾的通用工具[9]。针对滑坡易发性区划研究，一般可采用定性或定量方法。定性方法主要结合滑坡发展规律和专业人员知识经验展开评价分析，如专家经验[10]、层次分析[11]和加权线性组合法等；定量方法则是在统计数据基础上建立数学模型进行定量评价，如确定性方法、人工智能法和多元统计法等[12～17]。定性方法主要依靠先验知识，受主观影响较大，有较大局限性。因此，本文拟结合GIS和Logistic回归模型，从定量视角综合评判怒江州滑坡灾害易发性程度。

2 研究区概况

怒江州坐落于滇西北部，辖泸水市、福贡县、贡山独龙族怒族自治县和兰坪白族普米族自治县，包含29个乡镇，总面积 14 703 km2，其行政区划见图1(a)。怒江州地处青藏高原东南部横断山脉峡谷地带，山高坡陡，水系发育密集，立体气候突出，地质灾害分布广、突发性强，是云南省滑坡易发、多发地区[18]。

图1 怒江州乡镇区划图(a)和滑坡点分布图(b)

根据云南省地质调查局数据显示，怒江州共有两百余处滑坡隐患点。如图1(b)所示，滑坡灾害在怒江州下辖的4个县市中均有分布， 并且多分布于怒江河谷及其各级支流沿岸。其中，福贡县和泸水市的滑坡发育密度最大，贡山县分布较多，兰坪县分布较少。

3 滑坡灾害易发性因子分析

3.1 数据来源

滑坡灾害易发性区划研究是一个综合性评价的过程，滑坡影响因子选取的正确与否直接关系到评价结果的可靠度[19]。本文结合怒江州滑坡灾害形成机理和发育规律的已有研究[20]，并收集了怒江州自2008年起至2020年历史滑坡灾害点数据和相关的环境因子数据(如表1所示)作为研究数据。

表1 数据来源表

3.2 初始影响因子提取

滑坡灾害是多种内外因素共同影响作用的产物[21]。本文通过对怒江州地形地貌和滑坡隐患点分布情况的研究，发现主要影响因素是高程、坡向、地形起伏度、降水量、植被指数、距河流距离、河流密度、土地利用类型、坡度，如图2所示。

图2 滑坡灾害易发性影响因子分布图

(1)高程

高程表征了研究区的宏观地貌，研究表明地质灾害与高程分布具有明显的区域规律[22]。怒江州地势北高南低，以怒江为中心沿东西两侧延伸，受构造抬升，两侧地势逐渐增高。如图2(a)所示，滑坡灾害主要分布于高程较低的地区。

(2)坡向

坡向是某一地面点处高程变化量最大的方向[23]，对太阳辐射面影响较大。向阳一面坡体易导致岩体裂隙发育破碎，而阴坡土层多易于累积堆积。本文使用ArcGIS软件的坡向工具从DEM中提取坡向值如图2(b)所示，可见灾害点多分布于怒江流域两岸的不同坡向上，其中斜坡上分布的滑坡灾害更多。

(3)地形起伏度

地形起伏度反映了一个区域海拔最高点与最低点之间的差值，与滑坡灾害具有较强相关性。本文采用ArcGIS软件的Spatial Analyst工具，基于原始DEM数据分别计算出邻域内的最大值和最小值，将其相减得到如图2(c)所示的地形起伏度图。

(4)降水量

结合怒江州历年的滑坡灾害情况，强降水是造成滑坡灾害的主要致灾因子之一。怒江州降雨主要集中在每年4月～9月期间，其中7月～8月是强降雨时期，滑坡灾害多发。本文利用ArcGIS软件的栅格计算器，计算出2011年～2015年的平均降水量栅格数据(图2(d))。

(5)植被指数

植被对地质灾害发育和稳定性具有深刻影响，植被指数(Normalized Difference Vegetation Index，NDVI)反映了一个地区的植被覆盖情况，植被指数越大，表示植被覆盖程度越高。本文计算出2018年～2020年的怒江州平均NDVI空间分布情况如图2(e)所示。

(6)距河流距离

由于很多滑坡隐患点均分布于怒江及其支流沿岸，因此水系对怒江州滑坡灾害有潜在影响。以怒江、澜沧江、独龙江三大干流及通甸河、老窝河的水域边界为基础，依次建立 600 m、1 200 m、1 800 m、2 400 m和3 000 m的缓冲区，如图2(f)所示。

(7)河网密度

同时本文基于DEM数据，进行“填洼-流向-流量-栅格矢量化-线密度分析”等一系列空间分析，计算出单位面积内的河网密度图，如图2(g)所示。

(8)土地利用分类

土地利用分类情况在宏观上表征了怒江州不同土地类型的分布情况。本文利用中国科学院空天信息创新研究院发布的2020年30m地表覆盖精细分类产品，将怒江州的土地利用类型分为建筑用地、耕地、裸地、水域、林地5种类型。由图2(h)可见，滑坡隐患点多分布于建设用地区域。

(9)坡度

坡度与滑坡灾害的发生有着紧密联系，坡度越大的地方越易发生滑坡灾害。基于DEM数据，利用表面分析工具，计算出坡度图(图2(i))。

3.3 影响因子筛选

由于滑坡影响因子错综复杂且具有非线性特点，极大地影响了易发性分析的精度。因此，本文采用主成分分析方法筛选对比成灾因子，筛选出有效的影响因子进行易发性分析，提高准确度。

(1)将原始数据标准化，以消除量纲影响

本文为统一影响因子的数据类型和单位，采用极差标准化方法进行归一化处理，计算公式如下：

(1)

式中：xi表示各影响因子的值，xmax和xmin分别表示各影响因子的最大值和最小值。

(2)建立变量之间的相关系数矩阵R

R=(rij)m×m

(2)

(3)

(3)计算相关系数矩阵R的特征值λj(j=1，2，…，m)的信息贡献率和累积贡献率

(4)

(5)

bj为主成分yj的信息贡献率；αp为主成分y1，y2，…，yp的累积贡献率。当αp接近于1(αp=0.85，0.90，0.95)时，则选择前P个指标变量y1，y2，…，yp作为P个主成分进行综合分析。滑坡灾害影响因子特征值及主成分分析贡献率如表2所示，表中可以得出前6个主成分的累积贡献率达到93.307%，包括了9个因子的整体信息，所以本文选取这6个因子进行易发性分析。

表2 滑坡灾害影响因子特征值及主成分贡献率

4 滑坡灾害易发性区划

4.1 Logistic回归模型建立

Logistic回归模型，是一种因变量为二分类变量的回归分析。在滑坡灾害的易发性评价中，将各评价指标数据作为自变量，而灾害发生与否可用0(滑坡灾害不发生)和1(滑坡灾害发生)表征，是典型的二分类变量[24，25]。由于滑坡灾害影响因子为非线性变量，不适合用线性回归推导，因此，本文选用Logistic回归模型分析滑坡灾害易发性，其表达式为：

(6)

(7)

(8)

式中P(y=1|x1，…，xi)是发生滑坡的概率；xi为影响因子；εi表示滑坡影响因子的线性函数；α表示在没有其他因子影响下，发生滑坡与不发生滑坡之比的对数值；βk是逻辑回归系数，表示改变影响因子时发生滑坡灾害与不发生概率之比的变化值；p表示滑坡发生的概率。根据逻辑回归模型，本文假设滑坡灾害发生的概率为P，取值范围为[0，1]。以滑坡灾害发生概率为因变量，各影响因子x1，…，xi为自变量，建立Logistic回归方程，则滑坡灾害发生的概率为：

(9)

本文利用公式(5)计算滑坡灾害发生的可能性，数值越大，则发生滑坡的可能性越大，反之则越小。

4.2 二元Logistic回归分析

本文建立了矩形渔网覆盖整个研究区域，并基于怒江州边界提取了 13 958个网格，并将该网格连接归一化后的影响因子属性表及灾害隐患点，得到每个网格内的影响因子和滑坡灾害发生情况，最后以表格形式导出统计结果。由于Logistic回归模型采用的是最大似然估计参数法，为了保证结果的准确性，样本规模需要大于100，但是样本数量过大会令任何多元相关都会出现统计显著[26]。因此，本文采用随机抽样的方法抽取 6 000组作为分析样本，其中滑坡样本 1 425个，非滑坡样本 4 575个。本文采用SPSS数据分析软件进行二元Logistic回归分析，所得结果如表3所示。

表3 Logistic回归分析结果输出表

根据最终的拟合结果可知，P值<0.05。“B”为偏回归系数，“S.E”为标准误差，“Wald”是一个统计量，用以检验自变量对因变量是否有影响；“df”是自由度；“EXP(B)”为相应变量的OR值(又叫优势比，比值比)，表示在其他条件不变的情况下，自变量每改变一个单位，事件的发生比“Odds”的变化率。

根据表3影响因子回归系数，高程和地形起伏度的回归系数为负值，表示高程值及地形起伏度与滑坡灾害的发生呈反比关系，则在海拔较低、地形平坦的区域更容易发生滑坡灾害。除此之外，研究区多数滑坡灾害的发生与降水量和坡度大小有关，而表中降水量和坡度的逻辑回归系数最大，分别是3.304和3.376。而河网密度的逻辑回归系数达到3.206，表征了滑坡隐患点多分布于怒江流域两岸及支流的现状。综上所述，表明了二元逻辑回归模型的模拟结果与实际情况相符。

本文为进一步定性分析模拟结果的准确性，以回归分析的预测值P为自变量、滑坡灾害发生情况为因变量，在SPSS中进行受试者工作特征曲线(Receiver Operating Characteristic Curve，ROC)分析，得到相应的结果和ROC曲线分别如表4和图3所示。

表4 ROC分析曲线输出表

图3 受试者工作特征曲线

由图3可知，ROC曲线下方与坐标轴所围面积(Area Under Curve，AUC)达到0.766，根据ROC曲线特征，当0.7

4.3 滑坡灾害综合易发性区划

将表3每个影响因子的逻辑回归系数代入式(5)，可知滑坡灾害易发性模型的表达式为：

(10)

基于ArcGIS软件，根据式(6)计算得到研究区所有网格的滑坡灾害易发性概率，并采用自然间断点分类法将评价结果划分为低易发区、较低易发区、中易发区、较高易发区和高易发区五个等级，从而生成研究区域的滑坡灾害易发性区划图和各县市分布图，如图4所示。

图4 怒江州乡镇易发性区划图(a)及滑坡易发性区划图(b)

4.4 滑坡灾害易发性区划分析

根据滑坡灾害易发性区划结果，本文利用空间统计工具，分别计算5个易发性分区面积与灾害点数量，结果如表5所示。

表5 怒江州滑坡灾害易发性分区统计表

综合图4(a)和表5，怒江州滑坡灾害点与高易发区都集中在怒江流域两岸，并且高易发区和较高易发区共包含了196处滑坡灾害点，占滑坡灾害总数的78.088%；而低易发区和较低易发区主要分布在怒江州海拔较高地区，仅包含灾害点29处，只占灾害总数的11.554%。由此说明，本文基于逻辑回归模型计算得到的滑坡灾害易发性区划结果与实际灾害情况相吻合，图4(a)的易发性区划结果有较好的准确性。

本文为进一步研究怒江州滑坡灾害易发性情况，将灾害区划与研究区乡镇级行政区划边界相叠加，并统计出每个市县各滑坡灾害易发性等级所占面积百分比如图4(b)所示和表6所示。

表6 县市内各滑坡灾害易发性面积百分比(%)

综合表6和图4(b)分析结果可知，贡山县高易发区和较高易发区面积占比分别是13.713%、25.885%，主要集中在独龙江乡、茨开镇和普拉底乡，分布于怒江流域和独龙江流域两岸；而低易发区占比3.133%，主要集中于丙中洛镇和独龙江乡北部区域。怒江流域自北向南贯穿整个福贡县，滑坡灾害多发，高易发区和较高易发区占比达到65.596%。泸水市高易发区占比为18.584%，主要分布于怒江流域两岸的乡镇；而较高易发区面积占比达到38.682%，上江镇、鲁掌镇、片马镇等均有涉及。兰坪县中易发区面积占比达到43.798%，主要集中分布于澜沧江流域两岸的中排乡、石登乡、营盘镇、兔峨乡，而较低易发区主要分布于通甸镇、啦井镇和金顶镇。

5 结 论

本文以云南省怒江州为研究区域，根据研究区自然环境和滑坡隐患点分布情况，选取了高程、坡度、地形起伏度、植被覆盖等多个影响因子。由于影响因子错综复杂，容易影响易发性分析精度，本文利用主成分分析方法筛选对比成灾因子，最终选取了高程、坡度、地形起伏度、降水量、河网密度、距河流距离6个影响因子构建滑坡灾害易发性评价指标体系。采用Logistic回归方法，建立了滑坡灾害易发分区模型，将怒江州滑坡灾害易发区分为低易发、较低易发、中易发、较高易发、高易发五个等级，绘制了怒江州滑坡灾害易发区划图。研究结果表明：高程、坡度、地形起伏度、降水量、河流密度、距河流距离6个影响因子对怒江州滑坡灾害影响显著，本文划分的怒江州滑坡易发分区准确地反映了滑坡高隐患区域，可为怒江州滑坡灾害防治提供科学指导。后续研究工作将深入研究怒江流域滑坡成灾机理和特点，进一步完善滑坡易发分区评价模型，构建整个怒江流域滑坡易发分区评价模型。