三峡库区武隆段滑坡灾害驱动因子演变格局与人类活动的影响

张志兼，黄 勋,2，蔡雨微，傅镜羽，朱 悦，杨 锐，韩超群

（1.重庆师范大学地理与旅游学院,重庆 401331；2.三峡库区地表过程与环境遥感重庆市重点实验室,重庆 401331）

0 引言

一般来讲，滑坡的时空分布受限于内在控制因子和外在诱发因子[1]。对于三峡库区滑坡的内在控制因子，学界普遍认为滑坡受工程地质岩组影响较大，如侏罗系红层地区滑坡广泛发育[2−3]，也有学者分析了三峡库区滑坡与构造活动的关系[4−5]，以及斜坡结构类型与滑坡成因机制及空间分布的关系[6−7]。

外在诱发因子主要有降雨、地震、库水位升降以及人类工程活动等。大量案例已表明，随着人类活动对自然环境改造的加剧，滑坡发生频率日渐增高，人类活动已逐渐成为滑坡变形失稳的主要诱因[8]。例如，三峡库水位升降导致大量涉水老、古滑坡发生变形复活，代表性的有巴东县黄土坡滑坡[9]；大规模采矿导致的滑坡以2009年6月5日重庆武隆鸡尾山发生的大型滑坡-碎屑流为代表，体积7.00×106m3，造成74 人死亡[10]；还有部分滑坡是由人类工程建设切坡所致[11−12]，2020年7月13日发生在重庆武隆S529 省道旁的白马镇牛儿湾滑坡，体积超过1.80×106m3，造成4 栋房屋倒塌，93 人紧急撤离。

除了对单体滑坡的研究外，近年来大量学者从区域层面探讨了滑坡时空分布的驱动因子。Froude 等[13]对2004—2016年全球灾难性滑坡的分析发现，在不断增强的城市建设、非法采矿和肆意削坡的扰动下，滑坡数量与日俱增，人类活动极有可能超过降雨成为滑坡的第一诱因；Li 等[14]研究了2003—2014年全国28 省市土地城市化与滑坡的关系，发现滑坡数量随城市建成区面积和路网密度的增大而增加，尤其以我国西部地区最为显著；王新胜等[15]针对三峡库区重庆主城段1 703个历史滑坡点的研究发现，人口密度位居滑坡影响因子之首。这一观点在Li 等[16]的研究中也得到证实，除库水位升降外，距道路的距离和土地利用类型等人类活动因子，对三峡库区滑坡分布的影响也已超过降雨因子。

不单是针对三峡库区，国内外其他研究也已表明，在快速城镇化的背景下，不断增强的人类活动是控制滑坡分布的关键要素[17]。这一认识在经济发达的城市区域是适用的，但是在城镇化率较低、人类活动强度较小的农村县城是否可行？为此，文中选择三峡库区腹心的武隆区作为研究区，该区的城镇化率较低（仅为全市的67.4%，2018年），人类活动强度位列三峡库区（重庆段）倒数第二[18]。以重庆武隆1991—2015年330 处历史滑坡为样本点，运用重力模型和标准差椭圆模型描述滑坡的时空分布格局，利用地理探测器进行驱动力分析，试图解释该格局的演化机制，得出欠发达县域内的人类活动也是影响滑坡分布关键性因素的观点，其结论可为三峡库区滑坡灾害减灾防灾规划提供依据。

1 研究区概况与数据来源

1.1 研究区概况

武隆区地处重庆市东南部乌江下游、三峡库区腹心（图1）。地势东北高、西南低，以山地为主，相对高差700～1 000 m。气候属亚热带湿润季风气候，雨量丰沛、四季分明，多年平均降水量多在1 000 mm 以上，雨量在时间分配上不均匀，4—9月降雨占全年总量的70%以上。区内地质构造以NE-NNE 走向的褶皱为主，断层不发育，属新华夏构造体系和南北经向构造体系——川黔南北构造带。出露地层全为沉积岩，岩性以碳酸盐岩、砂岩、页岩和泥岩为主。

图1 研究区概况Fig.1 Research area

武隆境内地势起伏大、山高谷深、河流密布、降雨量大且时间集中，为滑坡提供了适宜的孕灾环境，武隆一直都是三峡库区地质灾害比较严重的地区，在重庆市地质灾害危险性分区中属于武隆—南川重灾带[19]。加之近年来，区内社会经济发展水平不断提高、城市建设规模逐步扩大，人类工程活动如道路开挖、矿山开采和水电站等建设活动日益频繁，人类对地表环境的干扰程度也不断加强，导致重大滑坡灾害时有发生。

长江干流库岸滑坡受库水位变动的影响大，这已是学界的共识[1,16]。武隆位于长江二级支流乌江的下游（图1），三峡蓄水前武隆站水位变化区间约为169～185 m，蓄水后水位变化区间缩小为169～177 m，且季节间变动幅度小（图2），区内滑坡受库水位升降影响较小，同时文中重在揭示城镇化背景下人类生活生产活动对滑坡分布的影响，因此可排除库水位升降的干扰。

图2 武隆站与三峡大坝（坝前）的水位对比Fig.2 Water level comparison between Wulong Station and Three Gorges Dam (before the dam)

1.2 影响因子选取

结合武隆区自然环境条件及社会经济发展状况，文中将影响滑坡分布的多个指标分为降雨、地质条件、地形地貌和人类活动等4 个方面共14 种因子[20−21]（表1）。

表1 滑坡分布的影响影响因子分类及说明Table 1 Classification and description of impact factors affecting landslide distribution

（1）降雨：降雨往往是滑坡的直接诱发因素，由于强降雨，斜坡土体含水量增大，斜坡岩土体静、动水压力急剧增大，降低滑动面的抗剪强度[22]，增大了下滑力，从而诱发滑坡。

（2）地质条件：选取工程地质岩组、距背向斜轴区距离、距断层距离3 个因子来衡量该区的地质条件。地层岩性影响深部岩土体的物理属性，斜坡上的岩体往往在软岩、软层的影响下容易出现塑性形变，不同岩性接触带容易出现片理风化作用，形成强烈的错动带和风化带[23]；背向斜轴区构造应力相对集中[24]，裂隙分布密集，层间常有构造滑移面，且在局部地段存在顺向坡；断层附近岩土体结构不稳定，在外力扰动作用下易发生滑坡。

（3）地形地貌：地表起伏度可反映地面相对高差，影响着地表物质的侵蚀、搬运、堆积等过程，在很大程度上决定了滑坡灾害的易发程度；斜坡的坡度越大，土体下滑力越大，斜坡稳定性越低[1]；河流、沟谷的两侧冲刷侵蚀作用强烈，斜坡前缘土体被带走，形成临空面，容易发生滑坡。土壤类型、土壤侵蚀强度和NDVI 等6 个因子可以综合衡量该区地表土体的抗蚀性和抗冲性[25]。

（4）人类活动

武隆区破坏地质环境的人类工程活动主要体现在不合理切坡、填土、盲目采矿和矿山弃渣不合理堆放以及滥垦滥伐，因此选取距道路距离、距采矿点距离、人口密度和GDP 等4 个因子作为衡量该区人类活动强度的指标[16,26−27]。

1.3 数据处理

研究区内共有330 个历史滑坡灾害点，数据来自“重庆市武隆县地质灾害排查项目”调查成果，包括20世纪80年代以来武隆境内滑坡的类型、发生时间、位置、主要诱因、威胁人户及财产等信息。文中首先利用1991—2015年的330 处历史滑坡，分析该时段滑坡的时空分布格局。其次，考虑到GDP、人口、道路等相关数据的匹配性，选取2001—2015年的287 处历史滑坡分析该时间段滑坡空间分布的驱动因子。上述14 种因子中，对于连续型因子进行离散化，并采用自然断点法分成9 级[15,28]；对于离散型因子则保持其原有分级，部分因子见图3。

图3 滑坡影响因子Fig.3 Landslide impact factors

2 研究方法

采用趋势分析、核密度分析、重力模型和标准差椭圆模型分析滑坡的时空分布格局，利用地理探测器对不同时段滑坡分布的驱动因子进行分析，揭示其演化机制。

2.1 滑坡时空分布格局分析

2.1.1 重力模型

重心变化是研究海量时空数据演变趋势的有效手段[29]，“重心”可表征地理要素的时空分布特征，重心移动方向指示了空间现象的“高密度”部位，而重心移动的距离反映了要素变化幅度的空间差异。定义武隆滑坡的重心坐标分别为：

xi、yi——第i个滑坡的空间坐标；

n——滑坡总数。

2.1.2 标准差椭圆模型

标准差椭圆可用于度量地理要素的中心趋势、离散和方向趋势等空间分布特征[30]。标准差椭圆的中心点表示所有要素的中心位置，长半轴表示地理要素分布的方向性，短半轴表示地理要素分布的范围。因此通过绘制滑坡空间分布的标准差椭圆，比较不同时间段的椭圆变化特征，可以揭示滑坡分布的年代际空间动态过程和总体趋势。计算公式如下：

式中：C——标准差椭圆轴长；

n——滑坡总数；

xi、yi——第i起滑坡的坐标；

2.2 空间分布驱动因素探测

地理探测器是基于地理空间分异立论所提出的探测空间分异性、揭示其背后驱动因子的一种统计学方法[31]。其核心思想认为：若某个自变量对因变量有重要影响，那么二者在空间分布上应该具有相似性。地理探测器已广泛用于台风[32]、山洪[33]、滑坡[34]和泥石流[35]等自然灾害领域。地理探测器包括因子探测、风险探测、交互探测和生态探测。

文中采用因子探测和交互探测对滑坡空间分布的驱动因子进行分析，以降雨、地质条件、地形地貌和人类活动等4 个方面共14 种因子作为自变量X，以滑坡核密度值作为因变量Y，利用不同时间段的q值来反映驱动力的演化过程，通过探测各变量之间的空间分异程度和交互作用影响，从而找到影响滑坡空间分布的关键影响因子。

（1）因子探测

因子探测器通过因子解释力衡量自变量对因变量的贡献大小，进而判断某自变量是否为因变量空间分异性的驱动因子，以及所占权重。通过比较因变量在自变量各类别中的方差与全局方差的关系来计算因子解释力的大小，计算公式如下：

式中：q——解释力，为某影响因子在多大程度上解释了 滑坡的空间分布，q的值域为[0,1]，q值越大，表 示影响因子对滑坡空间分布的解释力越强；

h——变量因子分类数，h=1,···,L；

Nh、N——类型h、全区的样本数；

SSW——类型方差之和；

SST——全区总方差。

（2）交互探测

交互探测器分别计算和比较各单因子q值及两因子叠加后的q值，以判断两因子是否存在交互作用以及交互作用的强度等[36]。X1∩X2是将X1和X2进行空间叠加后形成的新的空间因子，利用表2 的判别依据，确定因子间的交互作用类型。

表2 交互探测关系对应表Table 2 Interaction detection relationship correspondence table

3 滑坡灾害的时空演变格局

3.1 时间变化特征

以武隆境内330 个历史滑坡点数据为基础，分析1991—2015年滑坡逐年变化趋势及其演化过程（图4）。武隆滑坡发生频次按累积曲线变化特征可划分为1991—1997年、1998—2007年和2008—2015年3 个阶段性变化过程。在1991—1997年处于低发阶段，滑坡发生频次为0～7 起/年，平均每年发生2.28 起；在1998—2007年进入快速增长阶段，发生滑坡277 起，占总数的83.94%，滑坡发生频次为4～84 起/年，平均每年发生27.7 起；在2008—2015年处于缓慢增长阶段，滑坡发生频次为2～8 起/年，平均每年发生4.63 起。3 个异常高值均出现在快速增长阶段，分别为1998年、2003年和2007年，对应的滑坡数量为18 起、84 起和53 起，共计155 起，占总数的46.97%。

图4 滑坡累积量与降雨量的对应关系Fig.4 Correspondence between landslide accumulation and rainfall

从滑坡累积量与年均降雨量的关系（图4）来看，1991—1997年的年均降雨量仅为973.38 mm，1998—2007年增至1 325.04 mm，增幅超过30%，滑坡累积曲线斜率较上一阶段明显变大，与降雨量的变化趋势基本一致，2008—2015年间年均降雨量为1 306.79 mm，与1998—2007年相比波动较小，然而滑坡累计曲线斜率却明显减缓。有研究表明，1 300 mm 是三峡库区地质灾害暴发的年降雨量临界值，在年均降雨量超过1 300 mm 后，地质灾害累积数与降雨量关系曲线斜率明显变缓[37]。因此，我们初步判断认为，降雨对于滑坡空间分布的影响在2007年以后存在一定程度的降低，对此还需结合滑坡发生机理和相关数据，进一步论证分析。

3.2 空间变化特征

3.2.1 滑坡分布密度分析

针对武隆境内1991—2015年的历史滑坡点，使用ArcGIS 计算核密度估计值以展示其聚集区域与聚集程度。据核密度估计结果（图5），滑坡的空间分布存在明显差异，呈现出西、中、东3 个聚集区域：（1）西部：平桥镇—鸭江镇的带状聚集；（2）中部：羊角镇—白马镇—长坝镇的环形聚集；（3）东部：巷口镇—火炉镇的带状聚集。其中，白马镇西北部和巷口镇中部的灾害密度值>1.5，平桥镇、长坝镇、羊角镇、火炉镇西部和仙女山镇东南部等地的灾害密度值为1.0～1.5，鸭江镇、庙垭乡、白云乡、江口镇以及沧沟乡等地的灾害密度值为0.5～1，其余地区滑坡发生起数较少。

图5 历史滑坡点核密度估计图Fig.5 Estimated nuclear density map of historical landslide sites

3.2.2 滑坡重心及其变化趋势

将1991—2015年的武隆滑坡数据以5年为间隔，采用重力模型分析滑坡灾害重心的移动轨迹和移动距离（图6），利用标准差椭圆分析滑坡方向及分布（表3），以揭示滑坡空间分布的变化特征。

表3 1991—2015年武隆滑坡标准差椭圆参数变化Table 3 Variation of standard deviation ellipse parameters of Wulong landslide from 1991—2015

图6 历史滑坡灾害重心及方向分布特征Fig.6 Historical landslide hazard center of gravity and directional distribution characteristics

1991—2015年武隆滑坡的重心在29°21′N～29°23′N、107°32′E～107°43′E 之间移动。1996—2000年滑坡重心由白马镇向西北移动至羊角镇中部，2001—2005年向东北方向移动至羊角镇西部，2006—2010年由羊角镇向东南方向移动至巷口镇与仙女山镇交界处，2011—2015年由巷口镇向西北方向迁回至羊角镇南部。1991—2015年灾害重心移动7.6 km，重心整体向东迁移，其中1996—2000年、2001—2005年、2006—2010年和2011—2015年移动距离分别为6.2 km、8.0 km、11.4 km 和4.8 km，共计移动30.4 km，平均每年移动1.2 km。

1991—2015年的标准差椭圆方向在77.7°～101.4°之间，武隆滑坡的空间分布呈现出由西北—东南方向转向东北—西南方向的变化过程。标准差椭圆的长半轴整体呈先增后减的趋势，2001—2005年最长（27.3 km），2006—2010年最短（23.1 km）；标准差椭圆的短半轴整体呈增长趋势，1991—1995年最短（6.5 km），2006—2010年最长（13.4 km）。结合偏心率分析，武隆滑坡空间分布的方向性呈减弱趋势，武隆滑坡的空间分布由聚集走向分散。

4 滑坡时空分布的驱动力分析

4.1 单因子影响分析

考虑到人口密度、GDP 等数据存在一定滞后性及其与灾害数据的匹配性，仅选取2001—2015年的287个历史滑坡点为样本，将其按5年为间隔划分为3 个时间段，利用地理探测器定量评价因子对滑坡分布的影响，实现对武隆滑坡时空格局驱动力的演化分析，分析结果见图7。

图7 滑坡分布影响因子q 值的变化过程Fig.7 The process of changing q value of landslide distribution influence factors

从2001—2015年来看，各影响因子对滑坡空间分布的解释力均值最高的前5 位依次为：降雨（0.274）、工程地质岩组（0.220）、土壤类型（0.132）、河流（0.120）、GDP（0.116）；而背向斜轴区、断层、NDVI、采矿点、地表起伏度、坡度等因子的解释力均低于0.05，对滑坡空间分布的影响不显著。

从不同时间段来看，2001—2005年、2006—2010年和2011—2015年解释力最高的因子分别为降雨（0.382）、工程地质岩组（0.319）、工程地质岩组（0.204）。相比于2001—2010年，2011—2015年降雨、地质条件因子和地形地貌因子对滑坡空间分布的解释力均有不同程度降低，其中降雨和河流的降幅最为明显，而人类活动因子的解释力则逐渐增强。根据前文（图3、图6）不难看出，滑坡高密度区往往也是人类活动密集区，武隆区2001—2015年间各项社会经济指标均呈现不同幅度增长（表4），表明随着经济的不断发展，人类对地表环境的改造越来越强烈，人类不合理的工程建设活动如切坡、填土、采矿和矿山弃渣不合理堆放以及由于经济发展和人口增长，所带来的人均耕地面积的减少，植被破坏、水土流失等问题，使得人类活动对于滑坡空间分布的影响也越来越显著。

表4 武隆区2005—2015年社会经济相关统计数据Table 4 Relevant socio-economic statistics of Wulong from 2005 to 2015

4.2 因子交互影响分析

根据地理探测器交互探测结果，各时间段解释力最高的前10 位如表5所示，各因子叠加后其解释力明显提高，表明滑坡的空间分布是多因素共同作用的结果[34]。2001—2015年间，主要以降雨和工程地质岩组与其它因子的交互作用最为显著，其中2001—2005年降雨与土壤侵蚀的交互结果以及2006—2010年工程地质岩组与降雨和河流的交互解释力均超过0.5，表明其对研究区50%以上的滑坡空间分布具有解释力。

表5 交互探测后的q 值Table 5 q-value after interaction detection

该区岩土体均受水力侵蚀影响，2001—2005年有29%的滑坡分布于强度水力侵蚀区和剧烈水力侵蚀区，50%的滑坡分布于中度水力侵蚀区，在土壤侵蚀作用下，地表岩、土体等被剥蚀、破坏和搬运的过程加剧。随着降雨的入渗，土体发生软化形成软滑面，造成土体失稳进而引发滑坡。同时，从工程地质岩组分布来看，区内滑坡点多沿河谷呈带状分布，河谷地带岩性组合以灰岩与泥、页岩岩组为主，此类岩性组合往往存在软弱夹层，因差异风化使斜坡成陡—缓—陡形式，形成凹腔[38]，在后缘加载与强降雨触发双重影响下极易发生滑坡，因而滑坡灾害沿河流走向呈带状密集发育。

需要指出的是人类活动因子与其它因子的交互解释力在2001—2015年间比重逐步上升，其中2006—2010年，工程地质岩组与距道路距离的交互解释力高达0.458，同时，2011—2015年解释力较高的交互结果也多以人类活动因子为主，且所有因子交互后的q值都超过原因子之和，其解释力均呈非线性增长，进一步论证了上述人类活动对于滑坡空间分布的影响越来越明显的观点。

5 讨论

尽管人们对影响滑坡发生的人为作用并不陌生，但大多属于感性认识，文中通过对武隆历史滑坡中人类活动的影响过程进行量化，有助于理解灾害系统中“人”的地位和作用，有利于形成合理的人地协调观。基于文中的认识，在武隆区未来的滑坡减灾防灾工作中，应当重视人类活动的影响，重点关注城镇建设、道路开挖等活动对自然坡体的扰动，同时在这些人类活动密集的区域内，人员、建筑物等承灾体也相对集中，因此，灾害造成的损失往往都很大。这也提醒我们在滑坡灾害的风险评价，尤其是危险性评价中，需要着重考虑人类活动的相关指标。

同时文中在以下两个方面尚存在不足之处：（1）本文降雨因子仅采用了年平均降雨数据，且由于无法获取研究区地层产状数据，仅以工程地质岩组、断层和背斜向斜轴区作为地质影响因子稍有欠缺，同时人类活动因子的选取也有一定的局限性。（2）仅以滑坡密度表征灾害强度，同时以网格作为评价单元，忽视了滑坡面的整体性和连续性，滑坡与各因子间的内在空间关系和影响机制也未得到良好体现。文中仅对武隆区滑坡驱动因子的演变格局和人类活动影响进行了初步的探讨，上述不足之处将在今后的研究中进一步完善。

6 结论

（1）武隆1991—2015年的滑坡累计曲线在时间分布上呈现出“缓-陡-缓”的特征，以2008年为界，之前的滑坡发生速率随降雨量的增加而增长，2008—2015年，在降雨量保持稳定的情况下，滑坡发生速率明显减缓。

（2）在空间分布上，武隆历史滑坡呈现西、中、东3 个聚集区域；滑坡分布的重心集中于中部的白马镇、羊角镇和东部的巷口镇等地；标准差椭圆分析表明：武隆滑坡的空间分布呈现出由西北—东南方向转向东北—西南方向的变化过程，同时表现出方向性减弱和离散化的趋势。

（3）针对上述滑坡的时空分布格局，利用地理探测器找寻了2001—2005、2006—2010 和2011—2015 三个时段滑坡分布的驱动因子，结果显示，降雨和“软-硬”互层结构岩组与其它因子的交互作用是影响武隆区滑坡分布的主要影响因素；同时降雨、地质因子和地形地貌等因子对滑坡空间分布的解释力均有不同程度降低，而人类活动因子的解释力不断增强，已逐渐成为影响滑坡分布的关键影响因素之一，是未来减灾防灾重点关注的对象。