基于频率比—层次分析法的金沙江巧家段滑坡灾害易发性评价

王 盈, 袁仁茂

0 引言

山体滑坡是一种非常严重的自然灾害, 常常造成人员伤亡和财产损失。 滑坡的发生受到多种因素的影响, 如地质构造、 岩性分布、 地貌条件、 土地利用、人类活动等, 其中极端强降雨、 地震是诱发滑坡的常见原因[1-6]。 例如: 2008 年汶川8.0 级地震诱发了数以万计的滑坡灾害, 在约48 678 km2的范围内, 滑坡数量高达480 007 个, 滑坡覆盖面积约711.8 km2[7],这些滑坡摧毁房屋建筑, 堵塞河道, 造成大约2 万人员伤亡[8], 同时造成了巨大的经济损失。 2017 年7月, 四川省茂县叠溪镇发生特大高位顺层岩质滑坡,造成10 人死亡, 73 人失踪, 整个新磨村被毁灭的重大损失。 滑坡源区山体在1933 年叠溪地震中产生拉张裂缝, 之后在多次地震、 长期重力以及降雨作用下, 最终整体失稳破坏[9-11]。 在构造活跃、 地形陡峭、 软弱岩性分布的山区, 滑坡易发生, 但是在构造不活跃、 地质条件稳定的区域, 由于人类工程的开挖及加载, 如高速公路的修建、 煤矿的开采, 通常也会造成严重的滑坡灾害。 因此对滑坡灾害的预测和治理对减轻灾害损失、 减少人员伤亡具有重要意义。

滑坡易发性评价是指在一定地形条件的前提下,在某一区域内发生滑坡的可能性程度[12]。 目前, 随着GIS 的发展, 有很多的滑坡易发性的评价方法, 一般可以分为定性评价和定量评价。 定性评价一般基于专家经验, 直接从现场调查进行判断或者根据各影响因素分区图进行判断, 如层次分析法[13-17]、 专家系统法[18,19]、 模糊综合评判法[20-24]等; 定量评价方法主要有两种, 一种是基于物理模型和工程地质的方法, 建立在滑坡发生的物理机制的基础上, 使用数学模型来定量描述斜坡稳定状态[25,26]。 第二种是基于统计分析理论, 它不需要收集滑坡工程地质方面的数据, 而是根据历史滑坡失稳的影响因素的特征来分析未来滑坡的发生概率[27], 如频率比法[28-30]、 信息量法[31-34]、 概率指数模型[35,36]、 证据权法[37-39]、 逻辑回归[40,41]、 人工神经网络[42,43]、 支持向量机[44-46]等。 应用不同的评价方法, 有着不同的评价效果, 因此在区域地质灾害评价时, 需要充分了解灾害情况,然后选择合适的影响因子以及评价方法, 这样才能获得更好的结果。

本文以金沙江云南巧家段为研究区, 基于野外调查和多源遥感解译, 在获得金沙江巧家段滑坡编目的前提下, 将层次分析法(Analytic Hierarchy Process,AHP) 和频率比法(Frequency Ratio, FR) 相结合,建立了滑坡易发性评价模型, 分析滑坡灾害的发生于各个影响因子的相关性, 利用层次分析法对各个影响因子赋予相对权重, 研究金沙江巧家段的滑坡灾害易发性程度, 以期待为研究区建设和规划减轻滑坡灾害风险。

1 研究区概况

金沙江巧家段位于四川和云南两省的交界地区,介于东经 102°72′—103°28′, 北纬26°00′—27°28′,总面积约4 037 km2, 涉及宁南县、 会东县、 巧家县、会泽县以及东川区(图1)。

图1 研究区位置图Fig.1 Location map of study area

研究区气候属于亚热带季风气候, 夏季高温多雨, 冬季温和少雨, 在高海拔的山区还具有明显的垂直气候特征, 多年平均降雨量处于750 ～1 050 mm,而且由南向北, 降水量依次减少。 整个研究区背靠药山山脉, 整体地势东高西低, 地形以山地为主, 山高谷深, 海拔高差大, 金沙江贯穿其中, 主要岩性为冲洪积层、 板岩、 千枚岩、 灰岩、 粉砂岩、 砂岩、 白云岩、 页岩、 泥岩、 玄武岩等等, 地质构造比较复杂,主要受小江断裂带北段控制, 其次为则木河断裂和五莲峰断裂, 新构造活动强烈, 未来有发生强震的可能[47]。 该区域土地利用类型以林地、 草地和耕地为主, 大量分布红壤、 褐红土以及黄棕壤等多种类型的土壤。

研究区由于构造环境复杂活跃, 区域降水不均衡, 地形高差大, 岩性、 土壤及土地利用的差异, 导致崩塌、 滑坡、 泥石流灾害频发, 造成大量经济损失和人员伤亡。 例如巧家县城后山特大滑坡[48]及其他地质灾害[49], 还有罕见的蒋家沟泥石流, 每到雨季,多次冲毁农田道路, 堵塞江水, 淤埋村庄, 在野外崩塌、 滑坡灾害也非常常见。 因此, 很有必要对易发生地质灾害的研究区域进行评估, 采取适当的缓解措施。

2 研究方法

频率比法是一种常见的统计分析法, 将某一指标因子进行状态分类并计算各等级状态对滑坡的影响程度, 指的是在整个研究区内滑坡的面积比, 但是单一的频率比模型忽视了各指标因子的所占权重[50], 所以本研究采用频率比和层次分析法相结合的方法, 确定各个影响因子权重, 确定易发性指数并基于GIS 平台进行易发性分区。 一般步骤如下: 首先, 建立滑坡编目, 选择其中的70%作为训练样本, 并选择影响因子, 通过频率比方法, 计算滑坡的相对频率(mLRF), 然后利用层次分析法计算各个影响因子的相对权重, 接下来将各个因子的图层按照相对权重进行叠加, 得到易发性结果图, 最后, 利用剩余30%的滑坡, 再选择相同面积非滑坡数据, 通过ROC 曲线检验模型的准确性。

3 数据来源

滑坡数据库是进行滑坡灾害易发性评价的重要基础, 本次研究通过野外调查和借助高分辨率遥感影像解译获得, 遥感影像为2014 年的QuickBird 0.5 m 分辨率的正射影像图。 最终, 在野外调查和遥感影像解译的基础上, 得到了研究区的滑坡分布图, 共有滑坡3 573 个(图2), 滑坡面积达8.4 km2, 而滑坡以中小型浅层滑坡为主, 多沿着河谷及道路分布。 本次采用频率比—层次分析法进行易发性评价, 所选取的影响因子, 包括高程、 坡度、 坡向、 地形起伏度、 曲率、 多年平均降雨量, 土地利用等, 共13 个影响因子, 具体数据性质见表1。

表1 数据来源Table 1 The source and the attributes of geodatabase

图2 研究区滑坡分布图Fig.2 The distribution of landslide in the study area

4 易发性评价

4.1 滑坡影响因子的选取及滑坡相对频率(mLRF) 的计算

研究区内滑坡灾害共有3 573 个滑坡, 借助于ArcGIS 的地统计分析模块, 将滑坡随机分为两类, 一类为滑坡的70%, 共2 501 个, 面积约5.5 km2, 并作为该次模型的训练样本, 另一类滑坡1 072 个, 则作为模型的验证样本。 然后根据现有的数据, 选取影响研究区内滑坡发生的影响因子, 如高程、 坡度、 坡向、 地层岩性、 断裂、 土地利用类型等共13 个因子。确定好影响因子后, 在ArcGIS 平台上进行重分类,由此获得各个影响因子图层, 见图3, 然后通过公式(1) 来分析各个影响因子的分级因子对滑坡的相关性, 即滑坡的相对频率比, 见附表1。

图3 滑坡影响因子分级图层Fig.3 Landslide triggering factors and their classes

其中, mLRF 表示滑坡的相对频率比; LF 表示各个影响因子的各分级内的滑坡面积; CA 表示各个影响因子的各分级的总面积。 各个影响因子分类及统计结果如下:

1) 高程

在山区, 高程影响灾害的分布, 主要原因在于高程与降水量之间有着密切的相关性。 研究区高程范围在527 ～3 675 m 之间, 金沙江沿岸是地势最低的区域, 两侧沟壑发育, 地势增高。 在本文中, 将研究区高程按照每200 m 分为一级, 分为527 ～700 m,700～900 m, ＞900～1 100 m, …, ＞2 900 ～3 100 m,＞3 100 m 等14 级。 其中在＞700～1 900 m 范围内的滑坡面积较多, 在1 100～1 300 m 范围内的滑坡面积达到最大, ＞2 900 m 的滑坡面积最小。 滑坡的相对频率比随着高程的增加, 不断增加, ＞700～900 m 为0.162 0, ＞900～1 100 m 达到最大(0.163 8), 之后则随着高程的增加又不断下降, 到＞3 100 m 基本为0。

2) 坡度

坡度也是非常重要的影响因子, 表现在坡度影响地表水径流、 地下水补给和排泄、 物质的搬运和堆积、 斜坡体的应力分布等[51]。 研究区坡度多在0°～82°之间, 按照每10°划分为7 类, 分别为0 ～10°,＞10° ～20°, ＞20° ～30°, ＞30° ～40°, ＞40° ～50°,＞50°～60°, ＞60°。 在＞20°～50°分类内的滑坡面积较多, 在＞30°～40°分类内面积最大, 较大滑坡频率比分布在＞30°～60°, 这与其他研究也比较符合。

3) 坡向

坡向对滑坡的影响主要表现为山坡的小气候和水热比的规律性差异[52], 间接影响其他因素, 如降水,土壤湿度, 植被覆盖等。 研究区坡向分为5 类, 分别为平坦、 北、 东北、 东、 东南、 南、 西南、 西、 西北。 统计结果表明, 在西南、 南、 东南3 个方向的滑坡面积最大, 相对频率比也最大, 与前人所认为阳坡的滑坡易发性更高具有一致性。

4) 地形起伏度

地形起伏度代表了研究区的地形高差, 一般情况下, 地形起伏度越大, 边坡的重力势能越大, 越容易发生滑坡灾害。 本次采用变点分析法确定了最佳起伏度[53], 然后以50 m 为间距, 将研究区地形起伏度划分为8 类, 结果发现近80%的滑坡面积在100～200 m的分类范围内。 另外, 随着起伏度的增加, 滑坡密度不断增大, 到＞200 ～250 m 达到最大(0.320 2), 之后随着起伏度增加, 滑坡密度不断下降。

5) 曲率

曲率代表了坡面的形状, 曲率为正, 为凸形坡;为负, 为凹形坡; 为0, 为平坡。 研究区曲率在-2 ～ 2的范围内, 滑坡面积所占比重大, 另外比起凸形坡,在凹形坡山的滑坡面积较大, 且滑坡密度也相对较大。

6) 多年平均降雨量

降水是重要触发因素, 经常导致重大的滑坡灾害。 而它对边坡的影响主要表现为, 降水会进入破体裂隙内, 会软化物质, 减小接触面的摩擦力, 增加坡体重力, 动水压力和超静孔隙水压力的存在, 都可能会导致边坡失稳破坏。 因为研究区位于川滇两省的交界处, 且只有宁南、 巧家、 会泽和东川4 个雨量观测站, 因此本文收集了川滇全部雨量观测站自1981—2010 年的降水资料, 采用协同克里金差值, 以50 mm为间距, 获得研究区的多年平均降水分布图。 然后统计各个分类内的滑坡面积, 得到研究区降水多分布在700～1 100 mm 的范围内, 在800 ～900 mm 分类内的滑坡面积及滑坡密度最大。

7) 土地利用类型

土地利用主要反映人类活动的程度, 表示未利用土地、 林草地、 农业用地、 城镇用地等土地利用程度的高低[54]。 根据研究区土地利用的实际情况, 国家土地利用分类方法, 结合刘纪远等在建设“中国20世纪LUCC 时空平台” 建立的LUCC 分类系统, 划分为18 种二级分类土地利用类型, 其中林地和草地类型的滑坡面积最大, 但是裸岩石质地和其他建设用地类型的滑坡密度最大, 分别为0.358 9 和0.148 1。

8) 土壤类型

土壤的颗粒组成、 冻胀性、 透水程度和风化程度等都对滑坡有着重要的影响。 根据三级分类标准, 分为2 类, 其中褐红土和红壤分类的滑坡面积最大, 占到了滑坡面积的60%以上。 其他的土壤虽然也有分类, 但是滑坡面积相对较少。 尤其在棕壤性土(0.175 8)、 褐红土(0.137 1)、 红壤(0.113 1) 这3类土壤具有较高的滑坡密度。

9) 地层岩性

地层岩性是非常重要的影响因子[55], 因为不同的岩性有着不同的硬度和抗风化能力, 这对边坡的稳定性有着决定作用。 本文依据《工程岩体分级标准》(GB/T 50218—2014) 以及研究区岩性分布特征, 将研究区分为5 类, 分别为坚硬岩(花岗岩、 玄武岩、辉长岩、 石英砂岩、 辉绿岩、 硅质白云岩等), 较坚硬岩(白云岩、 灰岩、 板岩、 碳硅质板岩等), 较软弱岩(泥灰岩、 页岩、 泥砂质白云岩、 粉砂岩、 千枚岩、 泥质粉砂岩等), 软弱岩(泥岩、 泥页岩等),极软岩(残坡积、 冲洪积、 冰川、 湖河沉积等成因的砂、 砾石、 黏土等)。 对各级分类内的滑坡面积进行统计, 结果显示, 以较坚硬岩和较软弱岩分类的滑坡分布密度最大。

10) 河流

河流由于水的影响, 会浸润、 冲蚀、 掏蚀河岸,使得坡度变陡, 促进滑坡发育。 研究区滑坡有好多发育在河流两岸, 在ArcGIS 平台中, 以河流为中心进行缓冲区分析, 划分为5 个缓冲区, 分别为0～50 m,＞50～100 m, ＞100 ～150 m, ＞150 ～200 m, ＞200 m。结果发现, 随着与河流距离的增加, 滑坡面积逐渐减少, 但是滑坡密度在＞100～200 m 的范围内达到最大 。在＜200 m 的缓冲区内, 滑坡密度相对较大, 一旦超过200 m, 滑坡密度迅速下降, 这很有可能是超过200 m的范围, 河流的影响减小, 受到了其他因素较多。

11) 公路

在山区, 公路的开挖, 破坏了岩体完整性, 导致高陡边坡的增加, 破坏了原有边坡的应力平衡, 在卸荷作用下, 易发生崩塌、 滑坡灾害。 在ArcGIS 平台中, 以道路为中心进行缓冲区分析, 也划分为5 个缓冲区, 与河流缓冲区分类一致。 结果发现, 随着与道路距离的增加, 滑坡面积也逐渐减少, 而滑坡密度也具有递减趋势。

12) 断裂

断裂活动影响地形和岩体结构, 在断裂带沿岸,岩性一般较为破碎, 这是滑坡的物质来源, 而断裂活动导致的地震, 会在短时间内诱发大量的地质灾害,造成重大损失。 研究区内小江断裂带纵贯南北, 是主要的控制断裂。 通过以断裂为中心的缓冲区分析, 以1 km 为间距, 共划分为15 类, 总体趋势为随距离断裂越远, 滑坡面积迅速减少, 滑坡密度虽有起伏, 但总体上仍随距离越远而也不断减小。 由此说明, 距离断裂越远, 则受影响越小。

13) 地震峰值加速度(PGA)

PGA 是影响滑坡分布的地震因素之一, 一般峰值加速度值越大, 越容易发生滑坡, 影响越大。 研究区PGA 共分为3 类, 其中近98%的滑坡发生在0.2 ～0.3 g 范围内, 无明显规律性但是在0.3 g 的分类范围内的滑坡密度最大。

确定好影响因子后, 在ArcGIS 平台上进行重分类, 由此获得各个影响因子图层, 见图3, 然后通过公式(2) 来分析各个影响因子的分级因子对滑坡的相关性, 即滑坡的相对频率比, 见附表1。

其中, mLRF 表示滑坡的相对频率比; LF 表示各个影响因子的各分级内的滑坡面积; CA 表示各个影响因子的各分级的总面积。

4.2 层次分析法确定影响因子的相对权重

在确定了高程、 坡度、 坡向、 地形起伏度、 土地利用、 土壤等13 个影响因子后, 通过计算各个因子的滑坡相对频率比, 可以看到不同的影响因子以及各因子的分级因子, 都与滑坡的发育有着不同的影响程度。 那么为了比较各个影响因子对滑坡的影响程度,就需要确定各个因子的相对权重。 而层次分析法是一种多指标决策方法, 一般是借助专家经验对同层指标进行两两对比, 对相对重要性进行量化表示, 然后根据这些量化值构建判断矩阵; 之后, 根据判断矩阵计算各个指标的相对重要性, 并对结果做一致性检验,若检验结果合格, 则指标的权重分配合理, 否则重新进行比较, 最后得到最终评价结果。

原始层次分析法中运用专家经验法, 具有一定的主观性。 而频率比的方法则有效的减少了主观性的影响。 前面, 我们已经将13 个影响因子分级并计算出了各级影响因子的相对频率比, 然后基于滑坡的面积分布, 计算求得各个因子的相对频率的面积加权平均值(表3) 即各影响因子的相对频率比。 得到各影响因子的相对频率比后, 将结果两两相除, 就可以得到他们之间的相对重要性。 然后利用得到的结果构建判断矩阵(表4)。

表3 各个影响因子的滑坡相对频率比(面积加权平均值)Table 3 Relative frequency ratio of landslides for various influencing factors (area weighted average)

表4 影响因子的判断矩阵结果Table 4 Pairwise comparison matrix of all factors

由两两比较的结果看, PGA 与断裂距离的结果最大。 然后根据层次分析法计算影响因子的相对权重,即需要计算判断矩阵的特征向量(表5), 使得权重值和为1。 计算结果表明, 土地利用因子是最重要的影响因子, 权重值为0.153 9, 其次为与断裂距离、曲率、 高程, 权重值都在0.1 以上, 地震峰值加速度(PGA) 是权重值最小的影响因子。

表5 影响因子相对权重值Table 5 Relative weight value of influence factors

接下来, 对各因子的权重分配结果进行一致性检验(公式3、 4)。

其中, CR 为一致性指数, 一般要求结果小于0.1, 表示该矩阵通过一致性检验, 否则需要对判断矩阵进行修正; CI 为一致性指标, λmax为判断矩阵的最大特征值, n 表示判断矩阵的阶数; RI 为平均随机一致性指标, 取值如表6 所示(许树伯, 1998)。 经计算, 本文CR 值等于0, 远小于0.1, 因此各个影响因子的权重分配合理。

表6 平均随机一致性指标RI 取值表Table 6 Values of the random index (RI)

4.3 易发性评价

基于频率比—层次分析法的易发性评价, 主要借助ArcGIS 的空间分析功能, 将得到的各个影响因子图层以及计算求得的因子相对权重值, 代入到评价模型公式中(公式5), 利用栅格计算器功能得到易发性结果图。 最后按照自然间距分类方法, 将易发性指数值划分为5 类, 即低易发区、 中等易发区、 高易发区、 极高易发区(图4)。

图4 研究区易发性分区图Fig.4 Landslide susceptibility zones of study area

其中, OLSI 为易发性指数; Weight 为影响因子相对权重值; mLRF 为滑坡相对频率比。

研究区易发性结果图显示, 低易发区、 中等易发区、 高易发区、 极高易发区的面积分别占研究区的24.63%、 36.14%、 27.68%、 11.55%。 总体上, 研究区的60%以上属于中低易发区, 高等易发区主要分布在研究区南侧的金沙江及河谷两侧。

4.4 评价结果检验

易发性模型的检验和模型的构建一样重要, 若是模型的准确性不够, 则同样不能用于研究区的易发性评价。 目前, 对于评价模型的检验, 一般采用受试者工作特征曲线, 即ROC 曲线, 是根据一系列不同的二分类方式(分界值或决定阈), 以真阳性率(灵敏度) 为纵坐标, 假阳性率(1-特异度) 为横坐标绘制的曲线。 然后, 通过计算曲线下面积(AUC) 来测量模型的准确性。 若在AUC＞0.5 的情况下, AUC 越接近于1, 说明效果越好。 AUC 在 0.5～0.7 时有较低准确性, AUC 在0.7 ～0.9 时有一定准确性, AUC 在0.9 以上时有较高准确性。 AUC 越大, 利用滑坡模型预测的效果越好。

在本文中, 基于滑坡数量30%, 共1 072 个滑坡,对该易发性模型进行验证。 首先, 在ArcGIS 平台上,将这些滑坡转化为5 m×5 m 栅格的栅格单元, 共有104 559 个。 然后, 在研究区的非滑坡区域随机选取等量的栅格单元, 共同作为该模型的验证样本。 然后, 提取各个栅格点处的易发性指数值, 接下来, 数值导入到SPSS 软件中, 绘制ROC 曲线(图5), 并得到AUC 数值, 为0.827。 由此可见, 该易发性评价模型准确性较好。

图5 研究区易发性评价模型的ROC 曲线图Fig.5 Receiver operating characteristic (ROC) curve of susceptibility evaluation model in the research area

5 结论

金沙江巧家段滑坡灾害广泛分布, 经野外调查和遥感影像解译, 共有滑坡3 573 个。 本文将层次分析法和频率比法相结合, 建立了金沙江巧家县滑坡灾害易发性评价模型, 选取滑坡数量的70%作为训练样本, 其余的作为测试样本。 另外选取了高程、 坡度、起伏度、 多年平均降雨量、 土地利用、 土壤类型等13个影响因子, 统计这13 个影响因子与滑坡的分布关系表明: 在高程为700 ～1 100 m, 坡度为30°～60°,向阳坡, 地形起伏度为200～250 m, 在凹形坡的山谷地带, 多年平均降水量为800 ～900 mm, 裸岩石质地和其他建设用地类型, 以棕壤性土、 褐红土、 红壤为主的土壤, 较坚硬岩和较软弱岩分布, 地震峰值加速度(PGA) 为0.3 g 的区域的滑坡密度较大, 且随着与河流、 道路以及断裂的距离越大, 滑坡的分布密度不断降低。

通过层次分析法确定了各因子的相对权重, 其中, 土地利用、 与断裂的距离以及高程因子是权重较大的因子, 而PGA 是权重最小的影响因子。 基于ArcGIS 空间分析功能可得到研究区滑坡灾害易发性图, 表明: 研究区60%以上属于中低易发性, 而只有约12%区域为极高易发性, 并多位于研究区南侧金沙江河谷两侧。 最后利用相等面积的滑坡和非滑坡像元, 通过统计ROC 曲线的线下面积(AUC) 检验模型的准确性, 结果为0.827, 该模型是可靠的。