岷江上游生态脆弱性时空分异及驱动因子交互效应分析——以杂谷脑河流域为例

肖成志，计扬，李建忠，张志，巴仁基，曹亚廷

1. 中国地质调查局军民融合地质调查中心，四川 成都 610036；2. 中国地质大学（武汉）地球物理与空间信息学院，湖北 武汉 430074；3. 中国地质大学（武汉）工程学院，湖北 武汉 430074

岷江上游杂谷脑河流域是典型的山地型生态脆弱区（刘铁冬，2014），是中国重要的生态功能区，生态屏障战略地位十分突出，对维护长江流域生态安全，促进长江流域经济社会可持续发展具有举足轻重的作用（师嘉祺等，2021；汪璇等，2022）。该区域地形起伏大，易发生崩塌、滑坡泥石流等地质灾害（余恩旭等，2023），特别是干旱河谷区地质条件、自然条件非常复杂，生态十分脆弱，经济发展非常落后（陈小平等，2014），加上2008 年5•12 汶川地震影响，使得该区域生态脆弱性更加典型。因此了解岷江上游地震前后生态环境状况对该区域生态修复与重建至关重要。

生态系统脆弱性是一个生态系统在时间和空间上承受压力的能力（Williams et al.，2000）。早期的研究主要基于单一风险因子，以单维度、单要素的评价为主，随着脆弱性评价一步研究，逐渐发展为多维度、多要素的综合评价（李平星等，2014），其过程主要是选取生态脆弱性指标，在此基础上建立区域生态脆弱性的评估指标体系，然后进行全面评估。评价方法主要有专家打分法、层次分析法、空间主成分分析法、熵权法、物元法、人工神经网络法等（Niu et al.，2020；Qiao et al.，2020；Rezaeisabzevar et al.，2020）。生态系统作为一个动态系统，具有一定的自我调节机制，即适应能力（肖桐等，2010），不同特征生态系统受到的自然过程（如降水、极端事件等）和人类活动影响（如森林砍伐、绿化、土地使用变化等）不同（Wanyama et al.，2021）。如Polsky et al.（2007）认为人类和自然环境不是独立的系统是紧密耦合的，并且对威胁的暴露、敏感和适应能力不同，“脆弱性”的认知需要对人类与环境相互作用的生物物理、认知和社会维度进行调查，并在此基础上提出了暴露度-敏感性-适应能力（Vulnerability Scoping Diagram，VSD）模型。Wang et al.（2008）将生态脆弱性指标分为“土地资源条件”、“水热气象条件”和“地形条件”三大类。Kang et al.（2018）将生态脆弱性指标分为“自然环境”和“社会经济”两个代表性类别。以上研究大多以自然环境和社会经济方面来研究生态脆弱性，但对地质条件在生态脆弱性研究中起到的作用研究较少，未单独进行讨论。因此本文基于生态地质学理论和杂谷脑河流域生态本底条件，从“自然环境-地质地貌-地质灾害-社会环境”4 个维度选取评价指标，以探究复杂地质条件对研究区生态脆弱性的作用，为该区域绿色发展提供支撑。

1 研究区及数据

1.1 研究区概况

杂谷脑河流域（图1）位于青藏高原东缘，阿坝藏族羌族自治州东南部，岷江干流西部。地理位置31°26′—31°33′N，103°9—103°36′E，东西长约94 km，南北宽约56 km，流域面积为4 629 km2。区内山峦绵延、层岭叠嶂、山高谷深、地势西北高、东南低，区内坡度较陡多在20°—60°之间，流域内地貌类型为：侵蚀、深切河谷；龙门山剥蚀、侵蚀高山、中山；岷山、邛崃山剥蚀、侵蚀高山；岷山、邛崃山冰蚀、冰缘极高山、高山。气候相对复杂，垂向气候特征显著，水热环境具有明显的区域性和垂向变异特征，特别是干旱河谷区降水量较少。植被类型主要有针叶林、阔叶林、灌草，及樱桃树、李子树等。土壤类型主要有暗棕壤、棕壤、褐土、高山草甸土、亚高山草甸土、高山寒漠土等。

图1 研究区概况示意图Figure 1 Overview of the study area

1.2 数据源及处理

本研究数据包括数字高程模型（Digital Elevation Model，DEM）、Landsat 数据、土地利用数据、潜在蒸发数据、降雨数据、人口密度数据、地层岩性数据、地貌分区数据、地质灾害点数据、干旱河谷区边界数据、道路矢量数据和土壤类型数据。其中，DEM 数据为美国航空航天局（National Aeronautics and Space Administration，NASA）和日本宇宙航空研究开发机构（Japan Aerospace Exploration Agency，JAXA）共同开发的数字高程模型，可在地理空间数据云（https://www.gscloud.cn/）下载。Landsat 数据来源于美国地质勘探局（https://earthexplorer.usgs.gov/），考虑到数据质量、植被生长旺盛时期、1998 年洪水及2008 年5•12 地震影响，震前选择了2001 年6 月13 日数据，震后选择了2009 年6 月9 日数据，现状选择了2020 年7 月19 日数据，道路矢量数据由Google Earth 软件解译得到，其他数据来源如表1 所示。所有数据空间参考均统一为WGS_1984_UTM_48N，其空间分辨率统一至30 m×30 m 大小计算。

表1 数据源Table 1 Data Sources

2 研究方法

2.1 评价指标

根据指标体系构建原则，本文从“自然环境-地质地貌-地质灾害-社会环境”4 个维度选取了14 个生态脆弱性指标，指标计算和分级方法如表2 所示，生态脆弱性指标量纲不同，为统一量纲本文将所有指标分为5 级（图2、3）。有参考标准的参考相关标准和前人研究成果，无参考标准的根据现场实际调查情况或自然间断点法分级如表3 所示。

表2 杂谷脑河流域生态脆弱性评价指标体系Table 2 Ecological vulnerability assessment index system in the Zagunao River basin

表3 杂谷脑河流域生态脆弱性评价指标分级标准Table 3 Classification criteria of ecological vulnerability indexes in the Zagunao River basin

图3 稳定指标空间分布特征Figure 3 Spatial distribution characteristics of stability indicators

2.2 空间主成分分析

空间主成分分析可以通过对原始空间轴的旋转，将相关空间变量转换为少数几个互不相关的综合指标，从而实现对数据结构简化的目的（Abson et al.，2012；李志学等，2021），因此该方法可有效解决生态脆弱性评价指标之间相关性过大和指标冗余问题（孙宇晴等，2021）。本文选取累积贡献率在85%以上的主要成份，以其为替代因子，来评估生态脆弱程度，其计算公式为：

式中：

Wj——第j个替代指标（主成分）；

Yi——第i个初始指标的标准化值；

Zij——第i个初始指标在第j个主成分所对应的特征向量；

n——初始指标个数。

主成分分析需要对相关指标进行中心化和归一化处理。一般可将指标分为正面指标和负面指标两类。对于正向指标，指标值越大，脆弱性就越强烈，而对于负向指标，指标值越大，脆弱性就越低。最后生态脆弱性指数计算公式如下：

式中：

E——生态脆弱性指数；

Wj——第j个主成分；

Qj——第j个主成分的贡献率；

m——累计贡献率大于85%的前几个主成分个数。生态指数指数值的大小表示生态脆弱性程度，其值越大生态脆弱性越高生态环境越差，其值越小生态脆弱性越低，生态环境越好。为进行多年生态脆弱性的比较分析，对EI 进行归一化处理，归一化公式如下：

式中：

V——标准化值；

Ei——生态脆弱性指数第i个区域内的值；

Emin和Emax——生态脆弱性指数在研究区的最小值和最大值。

参考已有分级标准采用等间距法将生态脆弱性指数分为5 级：潜在脆弱（0—0.2）、微度脆弱（0.2—0.4）、轻度脆弱（0.4—0.6）、中度脆弱（0.6—0.8）、重度脆弱（0.8—1）。生态脆弱性值越大表示生态脆弱性越强，生态环境质量越差，生态脆弱性值越小表示生态脆弱性越低，生态环境质量越好。

2.3 地理探测器模型

地理探测器是一种基于数据挖掘和统计学方法的空间分析工具，旨在分析和解释不同因素对空间分布现象的影响，包括因素探测、相互作用探测、风险探测和生态探测4 个模块（Ren et al.，2022）。本文基于因子探测模块和交互探测模块对杂谷脑河流域的驱动因子进行探讨。因子检测的基本目的是识别地理差异的原因，即自变量X对生态脆弱性变化属性Y的解释力，并用下式进行计算：

式中：

q值——影响因素解释力；

h=1, 2, 3...,L——影响因子x分层，即分类或分区；

Gh和G——h层和整个区域的样本数；

W和T——层内方差和全区的方差。q的范围是0—1，q的值越大，因变量的解释力就越大。在极端情况下，q值为1 表示自变量的空间完全分异，而q值为0 则表示自变量随机分布。

交互检测可以定量表示两个驱动因素对生态脆弱性的分布格局的相互作用，表4 中概述的交互检测类型可用于识别各驱动因素之间如何交互作用。

表4 交互探测类型Table 4 Types of interactive detection

2.4 交互效应变异比

驱动因子相互作用效应可以看作是间接效应，是指两个驱动因子相互作用时对因变量的综合作用。交互效应变异比（Interactive Effect Variation Ratio，IEVR）是一种量化驱动因子的个体效应（直接影响）与交互效应（间接影响）之间差异的模型（Deng et al.，2021）。驱动因子IEVR 值越大驱动因子对因变量间接影响越强，IEVR 值越小驱动因子对因变量直接影响越强。其原理是计算某影响因子独立影响和与其他因子交互影响值的差值，然后计算该因子独立影响值和交互影响值所有差值的均值，表达式如下：

式中：

Rxi——变量xi与其他自变量的交互效应相对于其独立影响力的平均变异率；

n——因子个数；

P(xi)——xi的独立影响力；

P(xi∩xj)——变量xi和xj的交互影响力。

3 结果分析

3.1 生态脆弱性时空分布及演变

本文基于14 个生态脆弱性评价指标进行空间主成分分析，选取累计贡献率大于85%的主成分如表5 所示，基于公式（2）、公式（3）计算EVI，按照生态脆弱性分级标准得到生态脆弱性等级图（图4），由图4 可知研究区生态脆弱性变化活跃的区域主要为域干旱河谷区以及西部、北部高山极高山区域，生态脆弱性分布特征为西部和北部高于中部和东部，生态脆弱性随海拔升高呈现出先减弱后增强趋势。计算2001、2009 和2020 年3 个时期生态脆弱性面积占比（图5），3 个时期研究区均为潜在脆弱区面积占比最大面积占比最大分别为46.4%、41.6%和43.5%，重度脆弱区面积占比最小分别为7.1%、7.5%和5.3%，微度脆弱、轻度脆弱和中度脆弱区面积占比接近在13%—21.8%之间。图6 为各时期生态脆弱性转移桑基图，显示出2001—2009 年生态脆弱性主要变化类型为潜在脆弱转微度脆弱、微度脆弱转轻度脆弱、轻度脆弱转中度脆弱和中度脆弱转重度脆弱。2009—2020 生态脆弱性主要变化类型为，潜在脆弱转微度脆弱、微度脆弱转潜在脆弱、轻度脆弱转微度脆弱和中度脆弱转轻度脆弱。

表5 各主成分的特征值、贡献率及累计贡献率表Table 5 Eigenvalue, contribution rate and cumulative contribution rate of principal component

图4 2001、2009 和2020 年杂谷脑河流域生态脆弱性分级图Figure 4 Ecological vulnerability classification diagram of the Zagunao River Basin in 2001, 2009, 2020

图6 不同时期生态脆弱性转移桑基图Figure 6 Sankey diagram of ecological vulnerability transfer in different periods

图7a 为干旱河谷（a1），侵蚀、深切河谷（a2），龙门山剥蚀、侵蚀高山、中山（a3），岷山、邛崃山剥蚀、侵蚀高山（a4）和岷山、邛崃山冰蚀、冰缘极高山、高山（a5）在3 个时期的生态脆弱性均值，结果显示3 个时期差距较小但不同区域之间生态脆弱性具有明显差异。其中干旱河谷区和岷山、邛崃山冰蚀、冰缘极高山、高山生态脆弱性较强（生态脆弱性多年均值分别为2.27 和3.29），侵蚀深切河谷、龙门山、剥蚀、侵蚀高山、中山和岷山、邛崃山剥蚀、侵蚀高山生态脆弱性弱（生态脆弱性多年均值分别为1.28、1.15 和1.19）。

图7 干旱河谷区、不同地貌类型区生态脆弱性均值及面积比例Figure 7 Mean value and area proportion of ecological vulnerability in arid valley and different geomorphic types

从各区域生态脆弱性百分比来看（图7b），重度脆弱区只分布在岷山、邛崃山冰蚀、冰缘极高山、高山且轻度和重度脆弱区分布最广面积百分比为25.8%和36.3%。该区域内植被覆盖度低，山顶主要为雪被和裸地雪线以下主要为高山草甸，海拔高自然环境恶劣环境，生态环境极其脆弱，生态受损后不易修复，因此该区域需要将强保护，防止区域环境继续恶化，在保护的同时适地适树植树造林可能起到生态修复辅助作用。干旱河谷微度脆弱和轻度脆弱区分布最广面积百分比为46.6%和33.9%，主要位于汶川县城威州镇至理县县城杂谷脑镇。该区域人口密集，人类工程活动强烈，降雨较少，降雨时空分布不均，植被覆盖度较低，植被类型主要为灌草以及果树，岩性较为软弱易发生崩塌滑坡泥石流等地质灾害，多种因素综合导致该区域生态脆弱性变化强烈，该区域在农业活动、工程建造、资源开发时应更加注重生态环境保护。侵蚀深切河谷、龙门山剥蚀侵蚀高山中山和岷山邛崃山剥蚀侵蚀高山主潜在脆弱区分布最广面积百分比分别为74.5%、87.6 和83.8%生态环境较好。

3.2 生态脆弱性关键驱动因子分析

表6为因子探测结果表，结果显示高程（0.862）、生物丰度指数（0.77）、植被覆盖度（0.739）、水源涵养指数（0.724）和干旱指数（0.655）解释力均处于较高水平是生态脆弱性关键驱动因子，其余因子解释力均较低对研究区生态脆弱性空间分布影响较小。图8 为因子交互探测结果，各驱动因子交互大多为双因子增强和非线性增强不存在非线性减弱或者独立的关系，双因子交互q值均大于单因子q值，说明两因子交互作用解释力均大于单一因子解释力，其中植被覆盖度、干旱指数、水源涵养指数、生物丰度指数和高程与其他因子交互解释力最高，进一步说明植被覆盖度、干旱指数、水源涵养指数、生物丰度指数和高程对研究区生态脆弱性空间分布影响较大。

表6 生态脆弱性驱动因子探测表Table 6 Detection table of driving factors of ecological vulnerability

图8 驱动因子交互探测热力图Figure 8 Thermodynamic diagram of driving factor interaction detection

3.3 关键驱动因子交互效应分析

图9 为3 个时期生态脆弱性主要驱动因子交互效应变异比的变化，显示出2001、2009 和2020 年植被覆盖度IEVR 值别为2.41%、5.29%和18.73%，呈现逐渐增大趋势，一方面是2001—2020 年区域内中等植被覆盖度和较高植被覆盖度区域面积不断减小，另一方面是区域内土壤轻度侵蚀区域不断增加，地质灾害环境影响力不断加强。干旱指数在2001 年和2020 年IEVR 值相差不大分别为15.2%和14.43%，但在2009 年呈现出较低水准（8.09%），这是因为2009年相较于2001 年和2020 年气候变化更大，气候较为干燥，特别是干旱河谷区降雨较少。水源涵养指数IEVR 值分别为8.05%、12.82%和29.98%，IEVR 值逐渐增大，这是因为区域内中等植被覆盖度和较高植被覆盖度区域面积不断减小，而具有高质量生态系统的地区，如森林、灌木、草原和湿地等区域，地表径流较少，保水能力相对较高（Duan et al.，2022）。生物丰度指数IEVR 值在2001 年和2020 年较小分别为4.71%和3.87%，但在2009 年显著提升（8.17%），这是由于5•12 地震影响2009 年生物丰度和大大降低，地质灾害等因素影响力大大加强，而生物多样性决定了生态系统的抗变化能力，生物多样性越单一生态就越脆弱。高程IEVR 值在3 个时期均最低分别为1.48%、1.96%和2.12%，这是因为山区生态环境质量与地质环境息息相关（Liu et al.，2023），研究区垂直分异明显土壤类型和厚度、植被类型、降雨量、温度均随海拔的变化而呈现出不同特征。

图9 不同时期主要驱动因子交互效应比变化图Figure 9 Variation diagram of main driver IEVR in different periods

4 讨论与结论

4.1 讨论

1）本文选取14 个生态脆弱性评价指标，利用主成分分析法对杂谷脑河流域2001、2009 和2020年3 个时期生态脆弱性进行评价，并对3 个时期生态脆弱性变化进行分析。研究结果表明，受2008 年5•12 地震影响该区域植被大面积损毁地质灾害多发，加剧了该区域水土流失损害了原有山地生态系统功能（吴宁等，2008），相对于震前震后区域内生态脆弱性增强。自5•12 地震以来政府部门实施了许多生态保护修复措施，从而促进了杂谷脑河流域植被恢复（姚昆等，2020），现今区域内生态脆弱性降低。研究结果与杨渺等（2020）对5•12 地震极重灾区生态恢复进行评估结果类似。其中干旱河谷生态脆弱性较高和高寒地区生态脆弱性最高印证了前人干旱河谷区和高寒地区易形成生态脆弱去的结论。但本文生态脆弱性评价也存在局限性，一是指标分级存在主观性，如高程分级是按照现场调查情况和植被垂直分异划分，该标准只适用于杂谷脑河流域。二是本文评价指标较多，本文仅在评价方法上去除了评价指标的相关性，未对评价指标的相关性具体分析。因此，如何构建适用于不同区域的指标分级标准，深入分析各指标之间复杂关系避免评价指标冗余是下一步工作的重点方向。

2）研究区除主要影响因子外土壤侵蚀因素具有一定的影响力且不断加强，这与生态环境因素中植被覆盖度影响力逐渐减小相关，植被可以防止水土流失和土壤侵蚀，保持土壤湿度和养分促进水循环，如果土壤易于流失和侵蚀，水资源的净化能力下降，生态系统服务能力和生产力都会受到影响，生态脆弱性就会增加。地质灾害密度的影响力较小2009 和2020 年较2001 年其影响力稍有提高，这是因为地质灾害点主要区域位于杂谷脑河沿线特别是汶川县至理县县城地质灾害点较多，但对于整个流域生态脆弱性而言影响较弱。人口密度和离道路距离在所有时期的影响力均处低水平且变化不大，主要因为研究区人口分布密度较大的区域主要为杂谷脑河河谷区域，特别是威州镇至杂谷脑镇区域人口最为集中，区域内城镇建设和农业活动基本位于杂谷脑河河流阶地以及缓坡区域，对于研究区其他区域影响力较小。地质灾害与社会环境数据不够丰富，在今后的研究中需进进一步补充和分析。

4.2 结论

基于杂谷脑河流域环境背景和人文特征选取指标，利用空间主成分分析法得出5•12 地震震前、震后以及现状的生态脆弱性，利用地理探测器和交互效应变异比，得出了主要驱动因子及其交互效应。主要有以下结论：

1）杂谷脑河流域总体生态脆弱性较低以潜在脆弱区为主，但局部区域生态脆弱性较高。震后（2009 年）相对于震前（2001 年）杂谷脑河流域潜在和轻度脆弱区面积减少，微度、中度和重度脆弱区面积增加，生态脆弱性增强生态环境呈变差趋势，现今（2020 年）相对于震后（2009 年）潜在脆弱、微度脆弱和轻度脆弱区域面积增加，中度脆弱、重度脆弱区域面积有不同程度的减少，生态脆弱性降低生态环境得到改善。生态脆弱性变化强列区域主要位于研究区东部干旱河谷和西部北部及高山区域。

2）杂谷脑河流域不同区域生态脆弱性分布特征明显，东部干旱河谷区生态脆弱性相对较高，主要为微度脆弱和轻度脆弱区。重度脆弱区只分布在西部北部岷山、邛崃山冰蚀、冰缘极高山、高山，同时该区域轻度和中度脆弱区面积占比最大，生态脆弱性最高。其他区域潜在脆弱区面积占比最大生态脆弱性较低。

3）杂谷脑河生态脆弱性空间分异主要驱动因子为植被覆盖度、干旱指数、水源涵养指数、生物丰度指数和高程。对比主要驱动因子各时期IEVR值，结果表明高程是研究区生态脆弱性最直接影响因子，植被覆盖度和水源涵养指数对研究区生态脆弱性间接影响力不断加强，干旱指数在气候干燥时期对研究区生态脆弱性直接影响力增强，生物丰度指数对研究区生态脆弱性主要为直接影响，在震后对研究区生态脆弱性间接影响增强。