基于GDIV 模型的大渡河中游地区滑坡危险性评价与区划

阳清青，余秋兵，张廷斌，易桂花，张恺

（1.成都理工大学地球科学学院，四川 成都 610059；2.四川省冶金地质勘查局六〇六大队，四川成都 611730；3.成都理工大学旅游与城乡规划学院，四川 成都 610059）

0 引言

滑坡是极具破坏性的地质灾害之一，给人类社会造成了巨大的经济损失和人员伤亡，严重威胁人类生命财产安全，制约人类社会的可持续发展[1-2]。开展区域滑坡危险性评价，识别滑坡高危险区和阐明地质灾害的潜在分布对地质灾害防治具有重要意义，同时也是决策者预防和减少地质灾害损失的重要手段[3]。

滑坡危险性评价经历了从定性模型到半定量模型再到定量模型的发展过程[4]。定性模型主要是依据专家的知识和野外经验，受主观性限制较多；定量模型主要使用统计或数学建模技术，评价结果具有客观性[5-7]。随着科学技术的发展和定量风险管理的需要，滑坡定量评价模型成为主要的研究方向[8-10]。常用的定量模型主要包括信息量模型[11-12]、确定系数法[13]、逻辑回归模型[14-16]、地理加权回归模型[17]、支持向量机模型[18-19]、人工神经网络[20]和随机森林[21-22]等模型。上述模型为区域滑坡危险性评价提供了重要参考，已在省、市和县域尺度的地质灾害危险性评价中得到广泛应用[23-24]。但在评价过程中，多数研究在定性确定评价因子后，直接应用于危险性评价[25]，很少开展评价因子的定量优选工作；各评价因子分级量化主观性较强，评价因子权重赋值多依靠专家经验或先验知识，缺少客观依据[26]；运用常权栅格叠加各评价因子时，因子间存在抵消问题和极值区域埋没问题[27]。因此，结合区域滑坡成灾特征建立科学的评价指标体系，构建有效的评价模型是区域滑坡危险性评价的关键问题[28]。

大渡河中游地区位于青藏高原东缘的川西地区，区域地形地貌和地质条件复杂，生态环境较脆弱，气候复杂多变，是我国滑坡地质灾害高发地区[29]。由于地理环境恶劣，地质灾害调查和监测难度大，地质灾害综合防治能力较薄弱。综上，本文以大渡河中游地区为研究区，构建一种基于地理探测器、信息量法和变权栅格叠加的评价模型（geographic detector,information quantity and variable weight raster overlay model，GDIV）。通过地理探测器定量描述各因子对滑坡发生的影响力，基于影响力的显著程度实现评价因子的优选，同时以因子探测结果作为因子定权基础，增强权重的客观性，解决了定性评价因子选择主观性强的问题。采用变权栅格叠加评价因子信息量值，解决了评价因子叠加存在的抵消问题。基于GDIV 模型开展区内滑坡危险性评价，评价模型为区域滑坡地质灾害评价提供新参考，评价结果为研究区地质灾害防治提供决策支持。

1 研究区与数据源

1.1 研究区概况

大渡河中游地区地处四川省阿坝州和甘孜州交界处，区域面积约8 179.6 km2。该区地质构造复杂，包括滇藏“歹”字形构造体系和金汤弧形构造体系，出露的地层主要有三叠系、志留系和泥盆系。地处青藏高原东缘地势一级阶梯向二级阶梯过渡地带，区内海拔1 787～5 782 m，相对高差大。气候以高原型季风气候为主，年平均气温18.4℃，年均降雨量788 mm。区内水系发达，以大渡河为主流，包括革什扎河和东谷河等支流，河流汇集于丹巴县。在复杂的地质地貌条件、降雨和地震等因素影响下，地质灾害发生频繁。大渡河中游地区地质灾害类型主要包括滑坡、崩塌和泥石流，地质灾害详查数据显示，截至2021 年12 月，区内分布滑坡地质灾害隐患点313 处（图1）。滑坡发育规模以中小型为主，大型规模的地质灾害6 处，巨型规模的滑坡1 处。

图1 大渡河中游地区滑坡分布图和地质条件背景图Fig.1 Map of landslide distribution and geological conditions in the middle reach area of Dadu River

1.2 数据来源及处理

本研究所需的基础数据包括大渡河中游地区滑坡点、地层岩性、土壤类型、数字高程模型和地震动参数等。其中滑坡隐患点数据来源于四川省国土空间生态修复与地质灾害防治研究院1∶50 000 地质灾害详查数据；地层岩性和地质构造数据来源于1∶25 万炉霍县幅、康定幅和马尔康幅地质图和相关地质报告；降雨量数据来自国家气象科学数据中心—中国气象数据网（http://data.cma.cn/）；土壤类型数据来源于中国科学院资源环境科学数据中心（http://www.resdc.cn），基于2010—2020 年降雨量插值生成年平均降雨量；数字高程模型（digital elevation model，DEM）和2020 年8 月两景Landsat8 OLI_TIRS 数据来源于地理空间数据云（http://www.gscloud.cn/），空间分辨率30 m。利用DEM生成坡度、坡向、地形湿度指数（TWI）和径流强度指数（SPI），基于Landsat8 影像计算归一化植被指数（NDVI）。地震动参数源自中国地震动参数区划图GB 18306—2015（http://www.gb18306.net/）。同时将年平均降雨量、地震动参数等栅格数据统一空间分辨率为30 m，地理参考为WGS_1984_UTM_zone_47N。

2 研究方法

2.1 滑坡危险性评价初级指标构建

滑坡地质灾害是由内部因素和外部因素共同作用的结果[30]，本研究基于区域滑坡发育规律，从地质特征、地震、地形地貌、地表水系、降雨和人类活动等6 个方面，选取13 个影响因子构建大渡河中游地区滑坡危险性初级评价指标体系（图2）。其中工程地质岩组分为4 类，分别为坚硬岩、较坚硬岩、较软岩、松散土。坚硬岩主要为花岗岩、闪长岩、石英砂岩和石英岩等；较坚硬岩主要为板岩、石英岩状砂岩和硅质岩等；较软岩主要为页岩、千枚岩和粉砂岩等；松散土主要是第四系河漫滩等。

图2 大渡河中游地区滑坡危险性初级评价指标体系分级图Fig.2 Grading chart of the primary hazard assessment index system for landslides in the middle reach area of Dadu River Basin

2.2 地理探测器

将滑坡点密度作为因变量，高程、坡度、坡向、地震动参数、土壤类型、与道路距离和与断层距离等13 个因子作为自变量，运用地理探测器的因子探测器和交互探测器对因子进行计算分析。为消除因子类型差异的影响和优化各级因子对地质灾害影响程度的体现，先利用信息量法对自变量进行处理[31]：

式中：Xi——成灾因素X中的第i区间；

Y——成灾因素X中第i区间滑坡发生的信息量值；

B——滑坡事件；

Ni——研究区内包含评价因素Xi的单元数且存在滑坡的单元数；

Si——研究区内包含评价因素Xi的单元数；

N——研究区存在滑坡单元总数；

S——研究区中评价单元的总数。因子探测器的计算如下[32-33]：

式中：h=1,2,···,L——变量Y或因子X的分层；

Nh、N——层h和全区的单元数；

SSW、SST——层内方差之和和全区总方差；

q——因子对滑坡发育的解释力，值域为[0，1]，值越大说明对滑坡的影响力越大，因子对滑坡发育的解释力越强。

交互探测器是评估因子X1和X2共同作用时是否会增加或减弱对因变量Y的解释力，或这些因子对Y的影响是相互独立的。因子间若存在交互作用便能通过对比q值体现，交互作用方式的判断依据如表1 所示。

表1 交互作用探测器因子关系Table 1 Factor relationships of interaction detectors

2.3 评价模型

2.3.1 LRI 模型

LRI（logistic regression and informative）模型是逻辑回归模型和信息量模型的组合模型，逻辑回归模型是一种对二元分类变量进行回归分析的统计方法，定量描述二元因变量和一系列自变量之间的关系。在危险性评价时，选取评价因子为自变量，以信息量模型计算出的各个指标等级的信息量值表征各评价因子，以区内是否存在滑坡为因变量（1 代表存在，0 代表不存在）。构建逻辑回归函数：

式中：P——滑坡发生的概率，范围为0～1；

α——截距；

β——回归系数；

i——评价因子种类数目。

将式（5）两边取自然对数ln[P/1-P]，将评价因子Xi(i=1,2,···,n)作为自变量，建立线性回归方程：

2.3.2 变权栅格叠加

变权栅格叠加是变权理论在栅格叠加分析中的应用，该方法根据各叠加图层的初始权重及叠加栅格的像元值，通过一定规则对极值像元进行惩罚型、激励型或混合型变权，可以凸显个别极值栅格的效应[28]，实现层内因子重要区间的最大化应用。变权栅格叠加的实现主要包括状态变权函数确定和归一化变权权重计算，其过程如图3 所示。图3 中，W为各评价因子的初始权重，S(X)是根据状态变权函数获得各图层每个像元的状态变权值，W(X)是各个栅格像元的状态变权值和原始权重值乘积的归一化值，M是所有叠加图层各个栅格像元值和变权权重乘积的综合值。本文的状态变权函数采用局部激励型变权，当信息量值＞0.6 时，相应权重翻倍，即调权幅度为200%；当信息量值≤0.6 时，权重保持不变，即调权幅度为100%。

图3 变权栅格叠加过程Fig.3 The variational raster overlay process

2.3.3 GDIV 模型

GDIV 模型是本文提出的地理探测器、信息量法和变权栅格叠加的组合模型，该模型定量描述各因子对滑坡发生的影响力，并以定量结果作为因子定权基础，增强了权重的客观性；同时以变权栅格叠加因子值，消除了因子叠加存在的抵消问题。该模型的实现主要有四步：（1）基于信息量法，对初选因子分级和量化，增强因子各级的差异性；（2）基于地理探测器的因子探测和交互探测结果，选取影响力较强的因子作为评价因子；（3）构建权重方程，见式（6），归一化各评价因子q值以获得权重；（4）确定状态变权函数，生成各因子的变权权重，利用变权栅格叠加实现危险性评价。采用自然间断法将评价结果划分为地质灾害极高、高、中、低和极低5 个等级的危险区，实现滑坡危险性评价区划。具体流程如图4 所示。

图4 GDIV 模型计算流程图Fig.4 The flowchart of GDIV model calculation process

式中：ωi——归一化各因子q值所得的权重值；

qi——各因子因子探测q值。

3 结果与分析

3.1 影响因子分析

3.1.1 影响因子显著性分析

因子探测结果表明，地形地貌、地质特征、人类活动、地表水系、地震和降雨6 个方面的因子对地质灾害发生的影响力较大，但各因子对地质灾害发生的解释力存在差异（表2）。高程、土壤类型、年平均降雨量、地震动参数、与河流距离、工程地质岩组和与道路距离对滑坡发生的影响力较显著，这7 个因子是滑坡发育的重要因子，其q值均大于0.1，其中高程的平均解释力度为58.3%，土壤类型的平均解释力为41.5%，则高程和土壤类型是研究区滑坡的重要背景因子。与断层距离、归一化植被指数、坡向、径流强度指数和坡度对滑坡发生的解释力度相对一般，q值仅介于0～0.1。地形湿度指数p值大于0.05，未通过显著性检验，该因子对大渡河中游地区滑坡发生的影响程度较小，在本研究中为次要因子。

表2 滑坡初级评价指标q 值统计Table 2 Statistical analysis of primary evaluation index qvalues for landslides

3.1.2 影响因子交互作用

基于因子探测器中对滑坡发生解释力较强的7 个因子，采用交互探测器分析两因子交互作用后对滑坡发育的解释力（表3）。结果显示，21 对交互因子中，双因子增强16 对，非线性增强5 对，不存在相互独立或者非线性减弱的关系，与单因子对滑坡发育的解释力相比，两因子交互作用明显增强了对滑坡发育的解释力。交互作用解释力最高的是高程与工程地质岩组（q=0.736），其次是高程与地震动参数（q=0.676），解释力均高于65%，两者均是滑坡发生的重要因子组合。高程与土壤类型、工程地质岩组、与道路距离、地震动参数、年平均降雨量、与河流距离均呈现双因子增强作用类型，解释力均在59.6%以上，说明高程增强了土壤类型、工程地质岩组、与道路距离、与河流距离、地震动参数、年平均降雨量对滑坡发生的影响力。土壤类型与地震动参数、工程地质岩组的解释力均大于53%，是滑坡发育的较重要因子组合。从各因子的交互作用结果来看，高程、土壤类型和工程地质岩组等7 个因子的协同作用是大渡河中游地区滑坡发生的主导交互作用方式，也是滑坡形成的重要影响因子。

表3 部分滑坡初级评价指标交互作用Table 3 Interactions of primary evaluation indicators for landslides

3.2 评价因子优化与分析

通过地理探测器分析了研究区内高程、土壤类型和工程地质岩组等13 个因子与滑坡发生频率变化关系,结合因子探测器和交互探测器结果，最终选定影响力较强和交互作用显著的7 个因子作为大渡河中游地区滑坡危险性评价因子，即高程、年平均降雨量、土壤类型、地震动参数、与道路距离、工程地质岩组和与河流距离。

通过式（1）得到各评价因子信息量值（表4），其中高程值＜2 700 m 的区间信息量值大于2，表明区内滑坡主要发生在高程值＜2 700 m 的区间；土壤类型为淋溶土和人为土的区间信息量值大于1，表明区内滑坡主要发生在这两类土壤中。工程地质岩组通过基岩和堆积体的物理力学性质来影响滑坡发生，区内较软岩的信息量值大于1.8，意味着较软岩分布区更有利于滑坡的发生。归一化植被指数、与河流距离、地震动参数和与道路距离是滑坡形成的重要因子，年平均降雨量在750～775 时，信息量值最大；道路在建设中因剥离地表土壤和松散岩石对周边环境产生影响，区内滑坡发生概率与道路距离成负相关关系。区内与道路距离小于100 m 的缓冲区信息量值最高，随着与道路距离的增加信息量值依次递减。

表4 危险性评价因子分级与信息量值Table 4 Grading and information value of hazard evaluation factors

3.3 滑坡危险性评价

3.3.1 LRI 模型滑坡危险性评价

基于大渡河中游地区各评价因子的信息量值，运用LRI 模型分析得到各因子的逻辑回归系数值（表5）。结果表明，各评价因子的显著性水平值均小于0.05，通过显著性检验。各评价因子的回归系数均为正值，表明各评价因子对滑坡发生均具有促进作用。高程的回归系数最大（B=4.992），区内地质灾害的发生与高程相关程度最大。其次是土壤类型（B=3.001）、工程地质岩组（B=1.606）、年平均降雨量（B=1.103）、与道路距离（B=0.995）和地震动参数（B=0.802），回归系数均大于0.8，对滑坡的发生具有较强的影响。基于各因子的回归系数构建线性回归方程，实现区内滑坡危险性评价，见图5（a）。评价结果表明，随着危险等级的降低对应区内的

表5 滑坡危险性评价因子逻辑回归分析结果Table 5 Results of logistic regression analysis for landslide hazard evaluation factors

图5 滑坡危险性区划图Fig.5 Landslide hazard zoning map

面积逐步增加，其中极高和高危险的面积占比为18.1%，分布295 处滑坡，随着危险等级的降低对应区内的滑坡点也逐渐减少。

3.3.2 GDIV 模型滑坡危险性评价

利用式（6）获得各评价因子权重（表6），经变权栅格叠加实现基于GDIV 模型的滑坡危险性评价（图5b）。结果表明，各评价因子的权重大小顺序为：高程（0.329）＞土壤类型（0.234）＞年平均降雨量（0.103）＞地震动参数（0.092）＞与河流距离（0.089）＞工程地质岩组（0.088）＞与道路距离（0.065）。研究区极高和高危险区分布滑坡共302 处，占总地质灾害的96.5%。其中极高危险区滑坡数量占比为65.5%，高危险区滑坡数量占比为31%，中等及以下危险区滑坡发育较少。随滑坡危险等级的降低，对应区内的滑坡点数量也依次减少，评价结果与实际情况相符。滑坡危险性等级从低到高的各区面积占比分别是25.2%（极低）、40%（低）、13.2%（中）、11.2%（高）和10.5%（极高）。GDIV 模型充分考虑因子间的抵消效应和极值埋没问题，变权栅格叠加评价因子的综合信息量值在极值区间占比较小，极低危险区面积占比小于低危险区。极高危险区和高危险区主要沿东谷河、大渡河和革什扎河两侧分布，区内海拔小于3 200 m，以淋溶土为主，人类活动频繁。低危险区和极低危险区主要分布在水系发育程度较低、海拔相对较高和距交通干线相对较远的地区，集中分布在大渡河中游地区的阿科里乡、毛日乡、丹东乡、银恩乡、七美乡和协德乡等地。

表6 滑坡危险性评价因子权重值Table 6 Weight values of landslide hazard assessment factors

3.4 精度验证

基于ROC 曲线对GDIV 模型和LRI 模型滑坡危险性评价结果进行精度验证和对比分析。研究区共分布滑坡点313 处，在研究区生成313 处随机点作为非灾害样本点合并绘制ROC 曲线（图6）。LRI 模型的AUC值为0.903，标准误差为0.013，在满足显著性条件下AUC值最高可大0.929。GDIV 模型的AUC值为0.917，标准误差为0.012，在满足显著性条件下AUC值最高可大0.941。两种模型的评价效果均较好，但GDIV 模型的评价精度相对较高。由此表明基于GDIV 模型的滑坡危险性评价方法在大渡河中游地区的实用性更好，能较好地评估滑坡危险情况。

图6 滑坡危险性评价结果ROC 曲线Fig.6 ROC curve of landslide hazard evaluation results

4 结论

本文以大渡河中游地区为研究区，初选高程、土壤类型和与河流距离等13 个因子构建初级评价指标体系，利用地理探测器和信息量法实现因子的筛选和优化，构建GDIV 模型开展研究区滑坡危险性评价，为验证模型的精度引入LRI 模型，基于ROC 曲线对比分析和验证两种模型的评价精度，得出以下结论。

（1）地理探测器结果表明，高程（q=0.583）和土壤类型（q=0.415）是区内滑坡危险性的重要背景因素，高程与工程地质岩组（q=0.736）、高程与地震动参数（q=0.676）、高程与年平均降雨量（q=0.672）的协同作用对滑坡的发生具有重要影响。

（2）基于GDIV 模型的大渡河中游地区滑坡危险性评价结果显示，区内大部分地区为中度及以下危险区，面积占比为78.3%；而极高和高危险区面积占比为21.7%。极高危险区主要集中分布在河谷地区，工程地质岩组为松散土类均分布在该区内，土壤类型为淋溶土；高危险区要分布在东谷河、大渡河和革什扎河两侧极高危险区以外的区域；中危险区主要集中在海拔3 200～4 000 m范围内，极低和低危险区主要集中在人口稀少的高海拔区域。

（3）GDIV 模型（AUC=0.917）评价精度高于LRI 模型（AUC=0.903），在极高危险区和高危险区内GDIV 模型所含滑坡点总量较LRI 模型多7 处。GDIV 模型评价结果更符合区内滑坡分布情况，该模型适用于大渡河中游地区滑坡危险性评价。