GIS支持下CF与Logistic回归模型耦合的九寨沟景区滑坡易发性评价*

罗路广 裴向军 黄润秋 裴 钻 朱 凌

(地质灾害防治与地质环境保护国家重点实验室(成都理工大学)， 成都 610059， 中国)

0 引 言

2017年九寨沟“8·8”地震是四川省继汶川MS8.0级地震、芦山MS7.0地震以来的又一次强震(Fan et al.，2018)，诱发了大量的崩滑地质灾害，边坡岩土体结构受到不同程度的震裂损伤，加之区内降雨量大、人类活动频繁等特点使失稳概率急剧增高。震后九寨沟景区经过短暂的开放后，受强降雨诱发崩塌、滑坡、泥石流等地质灾害的威胁，再次临时闭园。因此，基于GIS平台的滑坡易发性评价可以为危险范围预测及重点防治区域区划等方面提供理论依据与技术支撑。

滑坡易发性评价为滑坡危险性和风险评价的基础，它是指特定区域内多个影响因子组合导致滑坡发生的可能性。国内外学者利用某种单一数学方法及模型进行地质灾害易发性分析和评价(胡瑞林等， 2013)，如信息量模型、确定性系数法(CF)、支持向量机法(SVM)、证据权法及逻辑回归法(Logistic)、层次分析法(AHP)、频率比模型(FR)等，均取得了较好的应用效果(许冲等， 2009； Pradhan et al.，2010； 张艳玲等， 2012； 刘丽娜等， 2014； 韩振华等， 2015； 许英姿等， 2016； 戴岚欣等， 2017)，然而仍存在易发性影响因子的量化、各因子之间及同一因子不同特征值之间对地质灾害易发性的影响值确定的人为因素干扰等缺陷。因此，两种方法甚至多方法、多模型相互耦合是近年来新的探究方向，如基于地理加权回归模型的云南省地质灾害易发性评价(饶品增等， 2017)、采用聚类分析和最大似然法的汶川极重灾区震后地质灾害易发性评价(胡凯衡等， 2012)、基于CF和SVM建立的CF-SVM耦合易发性评估模型(李远远等， 2018)、加权信息量模型(范林峰等， 2012)、信息量与Logistic回归耦合模型(樊芷吟等， 2018)等。沈玲玲等(2016)应用模糊逻辑法、信息量模型及Shannon熵改进的信息量模型，对岷县MS6.6级地震诱发滑坡进行了易发性评价； 黄发明等(2018)为保证SVM等机器学习模型所选的非滑坡栅格单元是真正的“非滑坡”，提出了基于自组织映射(SOM)的聚类分析和SVM耦合的滑坡易发性评价模型，结果显示SOM-SVM模型具有比单独SVM模型更高的成功率和预测率； 郭子正等(2019)基于逻辑回归-模糊层次分析方法(LR-FAHP)的加权频率比模型(WFR)，对三峡库区万州区内滑坡进行易发性等级预测，与单一的LR、FAHP和FR模型相比，WFR模型能将滑坡易发性评价精度提升4%～9%。

诸多研究结果表明，相比单一模型，多模型耦合在评价精度、合理性及成功率方面具有明显的优越性。传统确定性系数模型可客观地反映评价因子内部不同级别对滑坡易发性程度的贡献值，但是忽略了各因子之间关系的复杂性及其对滑坡易发性影响的差异性，不能很好地确定因子间的相对权重，而逻辑回归方法根据影响因子与历史灾害点的关系基于统计结果可较好地确定滑坡易发性影响因子之间的相互权重大小。因此，本文基于遥感影像解译和野外调查成果，分析地质灾害空间发育分布规律及其影响因素，选取地震峰值加速度(PGA)、距断层距离、高程、坡度、坡向、坡面曲率、地层岩性、公路及水系等9个因子，将确定性系数模型与逻辑回归模型相结合开展九寨沟景区滑坡易发性评价，为次生灾害防灾减灾提供参考。

1 滑坡易发性评价方法

1.1 确定性系数模型

确定性系数模型是由Shortliffe et al.(1975)提出的一个概率函数，由于评价流程相对简单，精度较高，在地质灾害易发性评价中应用较为广泛，其前提是假设已发生的地质灾害和未来将发生的灾害处于相同的地质环境条件下。计算公式为：

(1)

式中:Pa为滑坡在因子分类数据a中发生的条件概率，即a条件下滑坡面积与数据分类a的面积比值比表示；Ps为整个研究区的滑坡发生概率，即滑坡面积与研究区总面积的比值(本文Ps=0.0059)。

CF的取值范围为[-1，1]， 0～1表示滑坡发生的确定性高，该地质环境条件下滑坡易发生，CF值越接近于1，表示该单元对滑坡易发性越敏感； 反之，- 1～0表示滑坡发生的确定性低，该地质环境条件下滑坡不易发生； “0”表示该地质环境下不能确定是否会发生滑坡。

1.2 逻辑回归模型

逻辑回归模型是二项分类因变量常用的统计分析模型，它描述的是因变量滑坡是否发生(0代表不发生， 1代表发生)和多个致灾因子(x1，x2，x3，…xn)之间的关系。该模型中自变量可以是连续的也可以是离散的，不需要满足正态的频率分布。逻辑回归函数表达式为：

P=1/[1+e-(a+β1x1+…+βixi)]

(2)

式中：P为滑坡发生概率，范围为0～1；α为逻辑回归计算出的一个常数项；β为逻辑回归计算而得出的回归系数；i为评价因子种类数目。

将式(2)两边取自然对数，In(P/(1-P))作为因变量，将影响因子xi(i=1，2，…，n)作为自变量，得：

In(P/(1-P))=α+β1x1+β2x2…+

βnxn=α+βx

(3)

1.3 确定性系数与逻辑回归耦合模型

鉴于确定性系数模型未能考虑各评价因子对滑坡易发性影响的差异性问题，为提高评价精度，本文将确定性系数法与逻辑回归模型相耦合，利用逻辑回归法根据影响因子与历史灾害点的关系统计得到各评价因子的权重，采用确定性系数法计算各因子不同特征变量对应的CF值，并引入权重值将所有影响因子的CF值进行加权求和，滑坡易发性指数I的表达式为：

(4)

式中：I为滑坡加权易发性指数；ωi为第i个评价因子的权重；Ii为第i个评价因子不同级别对应的CF值。

2 研究区概况

九寨沟国家级自然保护区(图1)地处青藏高原东缘的深切割高山峡谷地带，是以高山湖泊群、瀑布群和钙华滩流为主体的国家重点风景名胜区及世界自然文化遗产中心，属于高原寒温带-亚寒带季风气候，年平均大气降雨量为704.3mm，降水多集中在5～9月。九寨沟沟域面积为655.49 km2，为嘉陵江支流白水河的主要支流，由则查洼沟(长约32 km，流域面积221.10 km2)、日则沟(长约28.55 km，流域面积245.88 km2)及树正沟(长约13.4 km)构成，3条支沟在平面上呈“Y”字型展布，形成一个完整的树枝状水系。地貌以高寒山地与峡谷为主，地势南高北低，最高海拔4828 m，最低海拔2000 m，平均海拔3000 m以上，平均坡度超过30°。出露地层主要有：泥盆系、石炭系、二叠系、三叠系灰岩、白云岩及第四系堆积物。

图1 研究区位置与滑坡分布图

研究区附近主要有岷江断裂、塔藏断裂和虎牙断裂的北西段(图1所示“8·8”地震推测发震断层； 李渝生等， 2017)，美国地质调查局(USGS)公布“8·8”地震峰值加速度(PGA)为0.08～0.26 g。区内构造作用复杂，断层发育(图1)，在构造应力作用下，岩体节理裂隙发育，为滑坡的产生提供了条件。2017年九寨沟“8·8”地震诱发了大量滑坡，根据震前、震后高清遥感影像对比及现场复核，目视解译滑坡1047处，滑坡总面积为3.88 km2。

3 九寨沟景区滑坡空间发育分布规律影响因子分析

滑坡的发生是斜坡自身内部基础地质条件与外界环境因素共同作用的结果。前者主要包括地层岩性、地形地貌及地质构造等，后者为水文地质环境、人类工程活动等触发因素。基于前人的研究成果及研究区地质条件，选取地质构造因子(地震峰值加速度、距断层距离)、地形地貌因子(高程、坡度、坡向、坡面曲率)、地层岩性、其他因子(公路及水系)等分析其对滑坡发育分布的影响(图2)。滑坡面密度是指滑坡面积与各级别因子总面积的比值，滑坡数密度是指单位面积内的滑坡数量。

图2 滑坡分布与影响因子关系

3.1 地质构造因子

3.1.1 地震峰值加速度

地震峰值加速度(PGA)是反映地震时地表震动强度的重要参数，震中附近值较高，随距震中距离增加而逐渐衰减。图2a显示滑坡主要分布在0.24 g和0.26 g 地震峰值加速度区域，占滑坡总面积的94.11%，随着PGA由0.2 g增加至0.26 g，滑坡面密度及数密度迅速增加。

3.1.2 距断层距离

区域断层构造对于地质体中节理裂隙的发育程度起着控制性作用，将岩体切割成块状或碎裂状，破坏了斜坡的稳定性。滑坡分布与距断层距离分布统计结果如图2b所示，滑坡主要分布在距断层2 km范围内，随着距断层距离的增加，滑坡数密度衰减迅速，表明岩体越破碎，滑坡越容易发生。

3.2 地形地貌因子

地形地貌条件是影响滑坡发育的重要因素，因此，在前人研究成果基础上，选择高程、坡度、坡向、坡面曲率等因子(由10 m×10 m的DEM数据提取获得)作为分析滑坡发育分布的影响因素。

3.2.1 高 程

高程是坡体应力值大小的重要影响因素，应力会随着坡高的增加而增加，影响着滑坡的势能。图2c为滑坡分布与海拔的关系，由图2c可知滑坡主要分布在2600～3200 m范围，该高程范围的面积仅占研究区面积的19.67%，却包含了69.21%的滑坡面积，其他高程范围滑坡面积及数量显著减少。其中高程2600～2800 m范围面密度最高达到2.786%，数密度为6.97个·km-2，高程2800～3000 m范围面密度次之， 3800 m以上基本没有滑坡发生。

3.2.2 坡 度

坡度是滑坡易发性评价的重要因子，平缓的山坡由于低剪应力使其发生滑坡概率变小，坡度越高，应力也随之增加，发生灾害的概率也将变大。图2d可知，滑坡主要分布于30°～60°范围，尤其在40°～50°达到峰值，随着坡度增大，灾害面密度呈现迅速增加，当坡度> 60°时，灾害面密度高达4.043%，数密度为8.79个·km-2，符合滑坡易发生在高坡度范围的规律。

3.2.3 坡 向

地震滑坡坡向上的分布情况与发震断层位置存在相关性，许强等(2010)将其受到地震波传播和断裂错动的影响作用归结为“背坡面效应”与“断层错动方向效应”。由图2e可知，滑坡主要发生在北东至南坡向范围内，滑坡面积占滑坡总面积的73.68%。推测发震断层(图1)走向为北西—南东，北东向至南范围与断层走向垂直或大角度相交且多处于背坡面，与汶川地震滑坡的分布呈现出相似的规律性。

3.2.4 曲 率

坡面曲率是对坡表一点扭曲变化程度的定量度量指标，正值表示斜坡为凸坡，负值表示凹坡，曲率为0或接近于0表示坡面越平坦。由斜坡曲率统计关系图(图2f)可见，坡面曲率越大，滑坡面积及其面密度、数密度越大，虽表现得不显著，但仍表明凹凸不平的斜坡较相对平坦的斜坡更有利于滑坡发生。

3.3 地质因子

地层是滑坡发生的物质基础，不同性质的岩石因坚硬程度和岩体结构的差异，直接制约滑坡发育类型与规模，对斜坡的变形破坏起着重要的作用。基于1︰5万地质图，根据地质年代划分为泥盆系(D)、石炭系(C)、石炭二叠系(Cp)、二叠系(P)、三叠系(T)及第四系(Q)。统计结果见图2g，石炭系岩组在所有地层中灾害最为发育，滑坡面积占总滑坡面积的68.6%，滑坡数密度为2.21个·km-2，第四系与石炭二叠系滑坡面密度次之，其次为二叠系、泥盆系及三叠系。

3.4 其他因子

3.4.1 公 路

景区内观光线路进行大规模的切坡、开挖等工程改变了斜坡的应力状态，加剧了滑坡的发生。在GIS平台中将研究区每隔300 m建立多环缓冲区，统计结果见图2h，滑坡面密度与数密度均随着与公路距离的增加迅速衰减，表明人类活动如公路建设对滑坡的发育具有一定的影响。

3.4.2 水 系

河流水系对两岸存在不同程度的冲刷和浸润现象，导致岸坡形成高陡的临空面，从而诱发滑坡。由图2i可知，滑坡主要分布在距水系600 m范围内，滑坡面密度随距水系距离增加由0～300 m范围处的1.19%降低至大于2100 m范围处的0.04%，规律性较为明显。

4 GIS支持下的滑坡易发性评价

4.1 评价因子选取与分级

根据前文中对各个因子对研究区滑坡发生的影响作用分析，共确定4类9个指标因子，构成本次滑坡易发性评价的指标体系，据前人经验与统计结果在GIS环境下对各个因子进行分级，结果见表1和图3。

图3 滑坡易发性评价因子分级图

表1 因子分级和CF、WCF值

4.2 评价因子权重的确定

将是否发生滑坡作为因变量(0表示滑坡不发生， 1表示滑坡发生)，将解译滑坡导入SPSS软件进行二元逻辑回归分析，因子权重由大到小依次为坡度(0.2316)、坡向(0.2245)、高程(0.1872)、地层(0.1263)、距断层距离(0.0742)、曲率(0.0414)、地震峰值加速度(0.0395)、距公路距离(0.0384)及距水系距离(0.0365)。

4.3 评价结果

根据式(1)计算各因子分类级别的CF值，计算结果见表1，将各因子的CF值进行叠加，得出基于确定性系数模型的滑坡易发性分区图(图4)。二元逻辑回归模型得到的滑坡易发性分区结果见图5。将CF值与权重值W代入式(4)，得到加权易发性指数WCF(表1)，在ArcGIS中将计算结果利用自然间断法把滑坡易发性划分为极高易发区、高易发区、中度易发区和低易发区(图6)，分区面积占比见表2。其中低易发区面积为397.99 km2，占总面积的60.72%； 中度易发区面积为158.51 km2，占总面积的24.18%。高易发区和极高易发区面积分别为64.83 km2和34.17 km2，占比为9.89%和5.21%。

图4 确定性系数模型评价结果

图5 Logistic回归模型评价结果

图6 耦合模型评价结果

表2 不同评价模型下的易发性分区统计结果

将确定性系数模型(图4)、逻辑回归模型(图5)和耦合模型(图6)获得的易发性评价结果进行对比发现，由3种方法所计算出的滑坡易发性评价结果与滑坡点分布规律可知分区结果主要有以下特点：

(1)滑坡高-极高易发区主要沿沟谷分布，日则沟比树正沟、则查洼沟滑坡更为发育，尤其集中在熊猫海附近，说明地震峰值加速度对滑坡的分布有明显的控制作用。

(2)基于耦合模型相比确定性系数与逻辑回归模型获得的易发性分区图中，更少的滑坡分布在低易发区，更多的滑坡分布在高-极高易发区，具有更好的实用性。

4.4 模型精度评价及检验

滑坡易发性评价结果是否准确直接关系到模型的可靠性，因此有必要对评价结果进行精度检验。受试者工作特征曲线(receiver operating charac-teristic curve，简称ROC曲线)是滑坡易发性评价模型性能检验的常用方法，以未发生滑坡的单元被正确预测的比例(假阳性率)为横坐标，以发生滑坡的单元被正确预测的比例(真阳性率)

为纵坐标绘制曲线，曲线越靠近左上角，说明模型分类的准确率越高。如图7所示，AUC值为ROC曲线以下至横坐标的面积，是评价模型准确度的指标，取值范围为0.5～1。由表2可以看出， 3种模型下滑坡分别有95.35%、93.37%及96.65%在中及以上易发区，仍有极少数滑坡分布在低易发区，这是由于滑坡受众多因素影响，本文模型仅考虑了最主要影响因素。AUC值分别达到了0.832、0.813和0.847，说明耦合模型比单一模型评价及预测能力更好。

图7 不同评价模型下的ROC曲线

5 结 论

(1)本文以九寨沟国家级自然保护区为研究对象，通过遥感解译和野外调查，发现九寨沟景区滑坡受断层控制作用明显，并沿沟谷及公路、水系相对发育。基于GIS平台，采用确定性系数法、逻辑回归及其耦合模型进行了滑坡易发性区划研究。

(2)基于前人研究成果及灾害发育分布规律，选取地质构造、地形地貌、地层岩组、其他因子等4类共9个因子采用3种模型建立滑坡易发性评价模型，其中坡度、坡向、高程、地层岩性4个因子对景区滑坡易发性影响相对较大。

(3)根据耦合模型下滑坡易发程度分区结果，将研究区划分为低易发区(60.72%)、中度易发区(24.18%)、高易发区(9.89%)和极高易发区(5.21%)，采用耦合模型得到的易发区划结果与解译及野外调查成果存在较好的一致性。结果表明采用耦合模型AUC值达到了0.847，比确定性系数模型、逻辑回归模型等单一模型评价结果更加合理、精度更高。

(4)由于地震滑坡具有长期效应，根据易发性评价结果，熊猫海及老虎海周边为滑坡高-极高易发区，地震产生的大量松散固体碎屑物质在强降雨条件下极易转化为新的滑坡或泥石流形成灾害链，为保障游客安全，建议针对重点防治区域进行有效的工程治理、监测及风险管控，降低成灾风险后开放此景点。