环境科学  2026, Vol. 47 Issue (10): 6933-6944   PDF    
耦合SHAP与地理探测器的青藏高原生态驱动归因
陆艺杰1,2, 张震2, 翁国明1, 胡克宏2     
1. 江苏省地质局大数据中心,南京 210007;
2. 安徽理工大学空间信息与测绘工程学院,淮南 232001
摘要: 为系统揭示青藏高原县域生态环境状况的空间分异规律与复杂驱动机制,以《生态环境状况评价技术规范》(HJ 192—2015)为依据,对各县域单元的生态环境状况指数(EI)进行测算. 进而耦合SHAP解释器与地理探测器模型,从特征贡献与因子交互等多维度解析其核心驱动因子. 结果表明:①青藏高原的生态环境状况呈现出由东南向西北递减的显著空间分异规律,EI优良等级区域集中于湿润的东南部林草过渡带,而较差与差等级区域则广泛分布于干旱高寒的西北腹地. ② EI的空间格局是多因子交互作用下非线性增强效应的直接体现. 其中,植被覆盖状况(NDVI和FVC)是居于主导地位的自然驱动力,其后依次为年均降水量(Pre)和土地利用类型(Land). ③交互探测证实,植被因子(NDVI和FVC)是关键的放大器,任一因子与之交互均能显著增强对EI空间分异的整体解释力. SHAP模型进一步将此非线性关系进行量化归因,精准揭示了各驱动因子在不同情景下的贡献方向与强度. 研究结论指出,青藏高原生态环境状况是一个对多重驱动因子产生复杂非线性响应的综合体,其宏观格局由植被和气候系统主导,并受到土地利用格局的显著调制.
关键词: 青藏高原      生态环境状况指数      SHAP解释器      地理探测器      空间分异规律      非线性增强效应     
Coupled Attribution of Driving Mechanisms for the Spatial Heterogeneity of the Ecological Environment on the Qinghai-Xizang Plateau via SHAP and Geodetector
LU Yi-jie1,2 , ZHANG Zhen2 , WENG Guo-ming1 , HU Ke-hong2     
1. Jiangsu Geological Bureau Big Data Center, Nanjing 210007, China;
2. School of Geomatics, Anhui University of Science and Technology, Huainan 232001, China
Abstract: To uncover the spatial patterns and complex drivers of the county-level ecological environment status across the Qinghai-Tibet Plateau, we calculated the ecological index (EI) for each unit, following the "Technical Criterion for Ecosystem Status Evaluation"(HJ 192—2015). A coupled SHapley Additive exPlanations (SHAP) explainer and Geodetector were then used to analyze the key driving factors through multi-dimensional analysis, including feature contribution and factor interaction. The results indicate that: ① The ecological environment status of the Qinghai-Tibet Plateau exhibited a significant spatial gradient, decreasing from southeast to northwest. Areas with "Excellent" and "Good" EI grades were concentrated in the humid southeastern forest-grassland transition zone, whereas those graded as "Poor" and "Bad" were widely distributed in the arid and alpine northwestern hinterland. ② The spatial pattern of the EI was a direct reflection of the nonlinear enhancement effect driven by multi-factor interactions. Among these, vegetation condition, characterized by normalized difference vegetation index (NDVI) and fractional vegetation coverage (FVC), stood as the predominant natural driving force, its importance far surpassing that of secondary factors such as mean annual precipitation (Pre) and land use type (Land). ③ Factor interaction detection confirmed that vegetation factors (NDVI and FVC) acted as a key amplifier, as their interaction with any other factor significantly enhanced the overall explanatory power for the spatial differentiation of the EI. Furthermore, the SHAP explainer quantitatively attributed this nonlinear relationship, precisely revealing the direction and magnitude of each driving factor's contribution under various contexts. This study profoundly concludes that the ecological environmental status of the Qinghai-Tibet Plateau is a composite system resulting from complex, nonlinear responses to multiple drivers, whose macro-scale pattern is dominated by the vegetation-climate system and significantly modulated by land use patterns.
Key words: Qinghai-Xizang Plateau      ecological index      SHAP explainer      Geodetector      spatial heterogeneity      nonlinear enhancement effect     

区域生态环境状况是生态治理目标的核心表征与解释变量[1]. 基于遥感物理量的评价方法,其局限性在于过度依赖单一指标,难以全面刻画生态系统的复杂机制与多维属性,当前动态耦合指标与行业标准化规范指数已协同构成区域生态环境评价的核心基准[2]. 《生态环境状况评价技术规范》(HJ 192—2015)[3]不仅构建了一套官方权威的评价指标体系与标准化核算方法,更标志性地推动了评估范式的深刻转型,即从早期以像元为基本分析单元、侧重于单一遥感物理参数(植被指数和地表温度等)的评价模式,演进为以行政区或流域为评价实体,系统性地耦合了自然生态系统状况、环境质量水平以及社会经济活动压力等多维要素的综合性诊断框架. 具体来说,HJ 192—2015通过衡量物种多样性计算生物丰度指数、基于遥感植被数据计算植被覆盖指数、整合河流和湖泊等水体数据计算水网密度指数、模拟土地侵蚀类型计算土地胁迫指数和量化工业污染物排放计算污染负荷指数,推动了县级、省级和生态区尺度统一评价[1],并已应用在全国多个重点生态县.

传统的线性回归等模型难以处理生态系统中普遍存在的空间异质性问题[4],地理探测器模型[5](Geodetector)基于统计学原理可解释地理空间分布规律以及驱动因素. 宋媛等[6]首次使用地理探测器对生态环境状况的多因素交互作用进行了分析. 胡克宏等[7]定量分析自然特征的驱动机制,增强了模型对空间异质性[4]的解释能力. 然而,地理探测器在精细刻画变量间非线性关系方面存在局限. 与此同时,以机器学习为代表的非线性模型虽因其强大的预测拟合能力在生态学研究中得到广泛应用,但其“不可解释”特性制约了对驱动机制的解析. SHAP解释器的出现解决了这一问题[8]. Lundberg等[9]首次将SHAP线性模型的系数估计Shapley值引入机器学习模型解释领域,用于分析特征对预测结果的贡献. Bae等[10]通过SHAP结合LightGBM模型,提取了空间效应,展现了SHAP在空间数据分析中的应用价值. 地理探测器侧重于从宏观层面识别空间格局分异的主导因素并判断因子交互的类型[5,11];而SHAP则精于从微观层面精准量化各因子复杂的非线性贡献,并能实现对具体分析单元的个体化归因[12,13],二者通过空间分异探测与非线性归因的耦合,构建起双重互补的驱动力分析框架,其结论的交叉验证与相互印证可显著提升驱动机制解析结果的可靠性与鲁棒性.

以往对青藏高原的生态评价多依赖单一的遥感物理量[14~16],本研究采用国家规范的生态环境状况指数(EI)作为评价因变量,评价结果不仅在科学上更全面,而且具有更强的现实指导意义. 具体而言,本研究借助EI指数表征呈现青藏高原地区的生态环境空间状况,并耦合基于机器学习的SHAP解释器与地理探测器两种前沿方法,构建一个优势互补的驱动力分析框架,以捕捉青藏高原生态环境复杂的驱动机制. 基于此,以期为青藏高原的生态保护和管理提供系统性和整体性的思维.

1 材料与方法 1.1 研究区概况

青藏高原,被誉为“世界屋脊”与“亚洲水塔”[17],是全球气候变化的敏感区和生态环境的脆弱区[18,19]. 从图 1海拔的空间分布可知,青藏高原整体呈现出西北高、东南低的总趋势,独特的地形地貌决定其气候[20]、水文[21,22]、植被[23,24]和生态系统[25]分布,进而导致了其脆弱的生态系统. 青藏高原是一个自然地理单元,其范围横跨中国西部的6个省级行政区,涵盖我国甘肃省、四川省、云南省、青海省、新疆维吾尔自治区和西藏自治区的县级行政区划220个. 从图 1各级行政边界可知,青藏高原的自然边界与省、市、县的行政边界并不完全重合,多个县级行政区被高原边界线分割,这为跨区域的生态保护和协同管理带来挑战.

图 1 青藏高原区划和自然地理概况 Fig. 1 Administrative divisions and physical geographical conditions of the Qinghai-Tibet Plateau

1.2 数据来源

本研究实验数据以2020年为基准,针对青藏高原220个县域构建统一的地理空间分析数据库. 为确保数据处理的一致性,所有栅格数据均被统一至30 m空间分辨率,并在此基础上采用针对性方法整合至县域尺度:首先,对土地利用和植被覆盖指数等栅格变量,运用分区统计聚合其像素级信息;其次,对源于统计公报的宏观数据,采用面积权重法进行空间降尺度;最后,对年鉴中的县域面板数据,经行政区划代码标准化后与GIS矢量底图进行属性连接. 详细数据来源如下所示.

(1)土地利用数据源自中国科学院资源环境科学与数据中心(https://www.resdc.cn/Default.aspx),采用专家目视判读方法解译获得(2000年国家大地坐标系),空间分辨率30 m,数据覆盖耕地、林地、草地、水域、建设用地和未利用土地等25个土地利用二级类型.

(2)植被覆盖指数是基于美国NASA数据中心MODIS 16 d 250 m连续时间序列NDVI数据产品(https://ladsweb.modaps.eosdis.nasa.gov/),采用最大值合成法生成的年度植被指数数据.

(3)数字地形高程模型(DEM)数据源自美国奋进号航天飞机的雷达地形测绘的SRTM V4.1数据(https://search.earthdata.nasa.gov/search?q=SRTM).

(4)河流长度、河流面积、水资源量和年均降水量数据源自甘肃省、四川省、云南省、青海省、新疆维吾尔自治区和西藏自治区水利厅发布的水土保持公报和水资源公报.

(5)以县区为单元的二氧化硫、化学需氧量、氨氮、氮氧化物及烟(粉)尘等的排放量和固体废物丢弃量数据源自甘肃省、四川省、云南省、青海省、新疆维吾尔自治区和西藏自治区统计局发布的统计年鉴数据(https://www.stats.gov.cn/sj/ndsj/index.html).

为清晰界定本研究的驱动因子并支撑后续的区域驱动机制解释,本研究将所有参与计算的数据定义为具体的特征,其特征名称和类型如表 1所示.

表 1 参与生态环境状况指数计算的各特征变量 Table 1 Characteristic variables involved in the calculation of ecological index

2 研究方法 2.1 生态环境状况指数

本文的EI计算方法依托《生态环境状况评价技术规范》(HJ 192—2015)确定. 评价指标中4类指数(生物丰度指数、植被覆盖指数、水网密度指数和污染负荷指数)根据前述的依赖数据并参考《生态环境状况评价技术规范》(HJ 192—2015)方法[3]完成计算:[EI=0.35×生物丰度指数+0.25×植被覆盖指数+0.15×水网密度指数+0.15×(100-土地胁迫指数)+0.10×(100-污染负荷指数)+环境限制指数]. 由于土地胁迫指数缺乏直接的来源数据,因此只能通过土地侵蚀类型结合土地利用类型间接获取. 土地侵蚀类型是基于DEM-FVC协同机制[26]推算所得,推算依据如表 2所示.

表 2 土壤侵蚀强度分级指标 Table 2 Soil erosion intensity classification index

土地侵蚀强度的计算见式(1).

$ \begin{aligned} & \text {ESI_values}=\alpha \cdot \text { Slope_values }+\beta \cdot(100-\text {FVC_values})+ \\ & \quad \gamma \cdot \text {Slope_values} \cdot\left(100-\mathrm{FVC}_{-} \text {value}\right) \end{aligned} $ (1)

式中,ESI_values表示土地侵蚀强度指数,值越高表示侵蚀越严重;Slope_values表示坡度,单位为度(°);FVC_values表示植被覆盖度,单位为百分比(%);α表示坡度权重系数,反映坡度对侵蚀的独立贡献;β表示植被缺失权重系数,反映植被缺失对侵蚀的独立贡献;γ表示交互作用系数,反映坡度与植被缺失之间的协同效应. 其中,FVC产品是根据年度植被指数数据通过像元二分模型反演形成[27].

将上述得到的生物丰度指数、植被覆盖指数、水网密度指数、土地胁迫指数和污染负荷指数,参考《生态环境状况评价技术规范》(HJ 192—2015)进行权重累加[3],得到青藏高原地区的EI;然后以青藏高原地区的县级行政区为单元获取数据,通过数据映射计算出区域统计范围内所有像元的平均值,以此反映该区域以县级行政单元为基础的生态环境态势.

根据EI,将生态环境分为5级,即:优、良、一般、较差和差[3]. 同样地,将生物丰度指数、植被覆盖指数、水网密度指数、土地胁迫指数和污染负荷指数也进行分级,具体见表 3.

表 3 生态环境态势分级 Table 3 Ecological environment status classification

2.2 生态质量特征贡献

基于ARIMA模型、LASSO回归模型、Ridge回归模型、Random Forest模型、Extra-Trees模型、AdaBoost模型、KNN模型和LightGBM这8类回归模型构建全样本的SHAP解释器,结合最优模型的基线预测值可视化EI的特征贡献度.

SHAP值是一种将模型对单个样本的预测结果归因于其输入特征的加性解释方法[28]. 其理论基础在于,将单个样本的预测值与全体样本的平均预测值(即基线)之间的差值,作为需要解释的总效应. 依据Shapley值原理将总效应“公平”分配给每个特征. 因此,每个特征的SHAP值精确地表示综合考虑所有其他特征的交互作用下,该特征将模型预测结果从基线值推向最终预测值的贡献量[16,29]. SHAP值的计算[30]见式(2).

$ \begin{aligned} \text { SHAP_values }= & \sum\limits_{S \subseteq\left\{x_1, \ldots, x_M\right\} \backslash\{i\}} \frac{|S|!(M-|S|-1)!}{M!} \times \\ & {[f(S \cup\{i\})-f(S)] } \end{aligned} $ (2)

式中,SHAP_values表示模型预测值;M表示特征的总数;S表示特征的一个子集,不包含当前特征{i};fS⋃i表示在特征子集S的基础上加入特征{i}后的模型预测值;fS表示仅基于特征子集S的模型预测值. 权重S!M-S-1!M!表示基于组合数学中的Shapley值分配原理,对所有可能的特征子集都进行公平且全面地考虑. 本研究将青藏高原地区涵盖的220个县级行政区划对应的生态环境态势作为特征子集进行模型训练.

2.3 生态质量驱动因素

地理探测器模型运用于驱动因素分析时,需将自变量离散化为类型量[7]. 对应生态环境状况指数的5类分级,本研究采用自然断点分类算法[31]将土地利用类型、植被覆盖指数、数字高程模型、植被覆盖度、水资源量、水域面积、河流长度、年均降水量、化学需氧排放量、氨氮排放量、氮氧化物排放量、二氧化硫排放量、固体废物丢弃量和烟(粉)尘排放量的数据值重分类为5类类型量,分别标记为类别1、类别2、类别3、类别4和类别5,将生态环境状况指数和所有影响特征输入地理探测器软件并进行因子探测、生态探测及交互探测,进而分析各个影响特征对生态环境质量的解释力及其交互作用.

3 结果与分析 3.1 生态质量现状分析

青藏高原整体生态环境状况呈现出极其显著的从东南向西北逐渐变差的空间分异规律. 生态环境优良区(蓝色区域)主要集中在青藏高原的东南部,包括四川西部的甘孜州、阿坝州南部,云南西北部的迪庆州,以及西藏东南部的林芝和昌都等区域. 生态环境一般区(黄色区域)主要分布在优良区和较差区之间的过渡地带,呈弧形带状分布,横贯青海南部、西藏中部等区域. 生态环境较差/差区(橙色/红色区域)广泛分布于青藏高原的西北部和腹地,包括西藏北部的羌塘高原、青海省的大部分地区(特别是柴达木盆地周边)以及新疆南部区域.

青藏高原EI的空间分异格局[图 2(a)],是其各分项指标[图2(b)~2(f)]综合作用的宏观体现,主要由自然生态本底与外部环境压力两大维度共同塑造.

图 2 青藏高原生态环境状况指数及其空间分布 Fig. 2 Spatial distribution of ecological index and sub-indexes in Qinghai-Tibet Plateau

(1)自然生态本底格局呈由东南向西北的梯度分异  由生物丰度[图 2(b)]、植被覆盖[图 2(c)]及水网密度[图 2(d)]等正向指标所表征的自然生态本底,共同呈现出由东南向西北显著递减的梯度特征. 这一格局主要受控于水热条件的宏观地理分布. 高值区集中分布于受印度洋暖湿气流影响显著的东南部,该区域降水充沛,水系发达[图 2(d)],孕育了以森林和草甸为主的高覆盖度植被群落[图 2(c)],进而支撑了较高的生物丰度[图 2(b)],使得整体EI等级评定为“优”或“良”. 低值区则广泛分布于气候高寒干旱的西北内陆,水分成为主要的生态限制因子,导致该区域水网稀疏[图 2(d)]和植被以荒漠和高寒草原为主[图 2(c)],生态系统结构相对简单,生物多样性维持在较低水平[图 2(b)].

(2)环境压力格局是自然胁迫与人为扰动的空间解耦  土地胁迫指数[图 2(e)]与污染负荷指数[图 2(f)]作为负向指标,揭示了影响生态质量的环境压力源,且二者的空间分布逻辑存在显著差异. 土地胁迫指数的空间格局与自然本底高度耦合. 高胁迫区[图 2(e)]主要位于自然禀赋恶劣的西北部,反映了该区域生态系统的内在脆弱性与低承载力. 相反,东南部生态系统稳定性及恢复力较强,土地胁迫程度较低. 污染负荷指数的空间分布则与自然本底呈现空间解耦特征,其格局主要由人类活动的强度与空间分布决定. 高负荷区[图 2(f)]并未出现在生态最脆弱的西北部,反而集中在人口相对密集、城镇化与工矿业发展水平较高的区域,如青海东部的河湟谷地及西藏中部的“一江两河”流域. 而广大的西部和北部地区,尽管生态脆弱,但因人类活动强度低,所承受的污染负荷也最低.

3.2 生态质量特征贡献分析

为保证训练/测试数据的真实性,本研究预测数据均不进行事先归一化处理,最后模型预测结果如表 4所示. 根据目视结果,该预测任务的真实数据具有高度非线性和剧烈波动的特点,传统的线性模型或简单的时间序列模型(ARIMA)无法胜任. 集成学习模型,特别是基于树的集成模型(LightGBM、Random Forest、Extra-Trees和AdaBoost),整体表现优于线性模型(LASSO和Ridge)和简单模型(KNN)[32],这表明目标变量与特征之间存在复杂的非线性关系,其强大的非线性拟合能力能够精准地预测目标值. 因此所有模型的性能排序为:LightGBM>Random Forest>Extra-Trees > KNN > LASSO>Ridge > AdaBoost > ARIMA,即LightGBM在所有模型中综合表现最佳,尤其是在预测精度(MAE,MAPE)方面具有明显优势,是完成此预测任务的首选模型.

表 4 8类回归模型预测结果的对比1) Table 4 Comparison of prediction results of eight regression models

基于LightGBM回归模型构建SHAP解释器,并对获取所有特征SHAP值的平均绝对值进行排序,表征模型预测影响最大的特征. SHAP解释器不仅验证了地理探测器的宏观发现,更从样本层面提供了对驱动因子贡献方向、非线性关系及数值依赖性的精细化、多维度归因. 图 3为SHAP依赖图对特征交互效应的可视化解构,揭示了NDVI的主效应呈显著正单调性,是关键驱动因子. 通过颜色编码量化特征间的交互效应可直观感知,高交互特征(SO2)能持续放大NDVI的边际贡献(SHAP值),这证实了二者存在强烈的正向协同作用. SHAP特征依赖图成功捕捉了特征间非简单线性叠加的条件依赖复杂关系,极大地增强了模型的可解释性与可信度.

横坐标为特征“NDVI”的特征值;纵坐标为特征“NDVI”的SHAP值;色柱为特征交互性,红色表示高交互性,蓝色表示低交互性 图 3 青藏高原生态环境状况指数各特征依赖图 Fig. 3 Dependence of each characteristic to the ecological index of Qinghai-Tibet Plateau

图 4分析揭示了青藏高原生态环境(EI)的形成是由自然要素强力主导、人为压力局部加剧的复杂过程.

上横坐标为每个特征的平均SHAP值,值越大,表示该特征对模型输出的平均影响越大;下横坐标为每个特征对单个预测的具体贡献,负值表示降低模型输出,正值表示提高模型输出;纵坐标为影响因素的重要性,百分数为特征全局重要性的量化指标;色柱为特征排名,红色表示高特征值,蓝色表示低特征值 图 4 青藏高原生态环境状况指数主要特征贡献度 Fig. 4 Contribution of important characteristics to the ecological index of Qinghai-Tibet Plateau

(1)植被状况与水热条件构成EI的正向驱动核心  SHAP特征贡献度图显示,植被状况与水热条件累积贡献度高达64.758%. 其中,NDVI(32.760%)和FVC(17.430%)合计贡献逾50%,确立了植被作为青藏高原生态环境状况的主要表征地位. 植被高特征值(红色散点)与正向SHAP贡献值呈现强耦合关系,定量证实了植被是构成高EI评分的主导因素. Pre(14.568%)作为关键水分来源,其贡献模式与植被因子高度一致,三者共同构成了控制高原生态基底的“水-植被”正反馈系统[25],是EI宏观格局的决定性力量.

(2)环境胁迫因子构成EI的主要负向驱动力,且可区分为自然源与人为源  SHAP特征贡献度图显示,自然源胁迫以Land(10.874%)和DEM(7.617%)为主要表征. 值得注意的是,Land和DEM高特征值(红色散点)对应负向贡献,这有力地表明不合理的人类活动和极端海拔等因素对区域生态环境状况构成了显著的负面压力[33,34]. 与之相对,如COD(2.618%)和NOx(1.805%)等人为源胁迫污染指标,虽全局贡献较小,但其作用方向清晰地指向环境压力,进一步补充和验证了人为胁迫在生态系统中的负面角色.

(3)驱动因子贡献分异揭示机制的层次性与局部性  SHAP特征贡献度图显示,模型驱动机制呈现出鲜明的层次性,特征重要性从NDVI(32.760%)到NOx(1.805%)的剧烈衰减,表明EI的形成并非由众多因子均等驱动,而是由少数全局性主导因子和多数次级影响因子构成. 其次,贡献度较低的因子揭示了驱动机制的局部性,如水体和污染物等平均SHAP值较低的因子,虽全局影响较弱,但在特定样本(地理单元)上可产生极大的正向或负向贡献,是局部生态异质性的关键驱动者. DEM的复杂双向贡献更是佐证了驱动机制的“全局主导与局地制约”特征,深刻反映了青藏高原生态系统驱动力在空间尺度上的异质性与复杂性.

为进一步解构模型对单个样本的预测逻辑,SHAP热力图(图 5)通过对所有样本进行相似性聚类排序,直观地揭示了不同生态状况区域的主导驱动模式. 分析结果清晰地将青藏高原的生态环境状况归纳为两大典型的驱动模式.

纵坐标为特征名;f(x)曲线为每个样本经过模型计算后最终预测值的变化情况,黑实线表示模型对每个样本的预测值,灰虚线表示模型输出的期望值(基准线),黑线高于虚线表示正贡献主导;低于虚线表示负贡献主导;黑色柱状图为每个特征对模型输出的总体贡献强度;色柱为特征重要性,颜色越深表示特征对该预测结果的贡献越大,颜色的正负表示贡献值的方向 图 5 生态环境状况指数各特征的基线预测热力图 Fig. 5 Thermal map of baseline prediction of each characteristic of ecological index

(1)高生态质量区的正向协同驱动模式主导样本集的主要部分(如图 5中样本序列0~55)  在这一模式下,模型预测输出的f(x)曲线总体处于相对较高水平,对应着生态环境状况“优”或“良”的区域. 其内在驱动机制表现为:首先,NDVI、FVC与Pre的SHAP值在此区间内呈现持续且强烈的正贡献,其数值普遍稳定在0.4以上的深红色区域. 这定量地表明,优越的植被覆盖、旺盛的植被活力以及充沛的降水,是维持高水平EI的协同正向驱动力. 其次,Land和DEM等固有负向因子的潜在威胁被有效抑制,其SHAP值多在中性区域(白色)徘徊,表明在这些样本中,较低的土地脆弱性和相对适宜的海拔条件未构成胁迫. 最后,污染负荷指标(COD和NOx)的SHAP值普遍接近于零,贡献微乎其微,进一步印证了该模式下的高EI值主要归因于卓越的自然生态本底.

(2)生态脆弱区的多重胁迫叠加模式解释生态环境的急剧恶化(如图 5中样本序列55~65)  与前一模式形成鲜明对比,该模式下的模型预测输出值急剧下降,骤降至0.3以下,对应生态环境脆弱或较差的区域. 其驱动机制呈现为一种负面效应的累积与放大:首先,NDVI、FVC和Pre的SHAP值急剧转为强烈的负贡献,其数值常常跌破-0.5,清晰地表明植被退化与水分亏缺是导致生态环境降级的根本原因. 更重要的是,Land与高海拔DEM的负向贡献在此模式下被显著放大,其SHAP值从近0水平加强至-0.4以下. 这一现象深刻揭示了胁迫的叠加效应,即在生态本底已然脆弱的条件下,高土地胁迫度和严酷高程环境的负面影响不再是线性叠加,而是被非线性地放大,共同将生态系统推向一个更低质的状态.

3.3 生态质量驱动因素分析 3.3.1 因子探测

表 5呈现了基于地理探测器因子探测模块的分析结果,该模块旨在量化各驱动因子对青藏高原生态环境状况指数(EI)空间分异的解释力,其解释力大小由q统计量衡量.

表 5 各特征因子探测 Table 5 Factor detector of feature factors

(1)自然地理要素是塑造EI空间格局的主导力量  分析结果表明,FVC的q统计量位列第一,表明其空间分层异质性与EI的空间分布具有最高的耦合度,是决定EI宏观格局的首要自然因子. Land和Pre紧随其后,同样展现出强大的解释力. 这印证了土地利用格局及其退化状态,以及由降水决定的水热条件,是构成高原生态环境本底的核心基石. WaterTotal亦表现出显著的解释力,说明水资源的宏观空间分布不均衡性是导致区域生态环境差异的另一重要原因[34].

(2)污染负荷因子和水系形态因子是塑造EI空间格局的次要及特殊因子  分析结果表明,尽管污染负荷相关指标的q统计量相对较低,但其P值均通过了显著性检验(P < 0.050). 这表明,虽然其独立解释的EI空间变异有限,但作为一种真实的环境胁迫,其影响在统计上是成立的. WaterArea与WaterLeng的q统计量为0,且P值远大于0.050,表明这两个因子与EI的空间分布不存在显著的统计学关联,其原因在于指标反映的是局部的水体形态,其生态效应在宏观尺度上被更能代表区域水分状况的主导变量Pre与WaterTotal的效应所掩盖或涵盖[35,36].

(3)关于DEM因子的方法论辨析  一个值得关注的现象是,DEM的q统计量相对较低(q=0.142),这与前述SHAP分析中揭示的其较高重要性似乎存在矛盾. 这种差异源于两种分析方法的内在机制不同:地理探测器通过比较因子分层后各层内的方差与总方差来衡量该因子对因变量空间分布的线性解释力. SHAP值则量化了在机器学习模型内部,一个特征对单个预测的边际贡献,能有效捕捉复杂的非线性与交互效应. 因此,DEM的低q值表明其简单的分层并不能很好地解释EI的空间格局,但其在SHAP中的高重要性则揭示了,在模型复杂的函数关系中,DEM通过与其他因子的非线性交互,对最终的EI预测值产生了关键影响. 这凸显了多模型、多视角分析对于全面理解复杂地理现象驱动机制的必要性.

3.3.2 生态因子

表 6展示了基于地理探测器生态探测模块的分析结果. 该模块通过比较不同驱动因子对生态环境状况指数(EI)空间分布影响的q统计量,判断两个因子的影响模式是否存在显著差异.

表 6 各特征生态探测1) Table 6 Ecological detector of feature factors

(1)驱动机制相似的因子对(无显著差异)  空间影响模式上无显著差异的因子通常存在强烈的空间耦合、共线性或同源性. FVC(表征“量”)与NDVI(表征“质”)之间无显著差异. 这表明在宏观尺度上,二者高度协同,植被覆盖广阔的区域其生长活力通常也较高,因此对EI的空间影响模式具有同质性. DEM与WaterArea、WaterLeng之间无显著差异. 这反映了青藏高原水系分布受地形地貌的强控制作用[22],河流多沿河谷发育,湖泊则分布于构造盆地,导致水系的空间格局与高程分异高度耦合,因此其生态效应在空间上呈现出相似的模式. 多种污染物指标之间无显著差异,暗示了指标间具有共同的来源或相似的扩散路径. Pre与Land之间无显著差异,这揭示了降水格局对土地利用类型的决定性制约. 青藏高原从东南湿润区到西北干旱区的降水梯度,决定了从森林、草甸、农田到荒漠的土地覆被梯度,从而使两者的空间影响模式高度重叠.

(2)驱动机制异质性的因子对(存在显著差异)  NDVI、FVC、Pre和Land这四大核心因子,与绝大多数其他因子之间均存在显著差异. 这证明了核心因子作为独特且不可替代的主导驱动力的地位,即分别从地表、生物和气候等关键维度塑造了生态格局,其作用机制与次要因子(污染物)存在本质区别. 所有污染物指标与NDVI、FVC、Pre、Land、DEM和WaterTotal等主要自然因子之间均存在显著差异,揭示了自然生态过程与人为污染过程在空间影响模式上的根本性分异. 自然要素决定了区域宏观尺度、呈梯度或面状分布的生态背景;而污染物则更多地反映了在特定区域(城市、工矿区)呈点状或廊道式分布的局部人类活动胁迫. 二者的作用尺度、强度与空间模式截然不同,构成了影响生态环境的两种异质性力量.

3.3.3 交互探测

图 6呈现了基于地理探测器交互探测模块的分析结果,旨在揭示不同驱动因子对青藏高原生态环境状况(EI)空间分异的综合影响. 结果普遍表明,任意两个因子的交互作用解释力均大于其单独作用,证实了EI的空间格局是由多因子协同驱动而非单一因子主导,反映了研究区生态系统的内在复杂性.

1.Land,2.NDVI,3.DEM,4.FVC,5.WaterTotal,6.WaterArea,7.WaterLeng,8.Pre,9.COD,10.NH3-N,11.NOx,12.SO2,13.SOL,14.YFC;下三角为热力q值(颜色越“红”,q越大),上三角为交互q值,数值越大,交互越强;★为Top10交互组合 图 6 各特征交互探测热力图 Fig. 6 Thermal map of interactive detection of each feature

(1)交互探测结果凸显了植被指标(NDVI和FVC)的核心中介作用  当植被与几乎所有其他因子交互时,其联合解释力(q值)均达到极高水平. 这有力证明,地形地貌(DEM)和水热条件(Pre)[37]等环境因子的生态效应,在很大程度上通过影响植被的生长活力与覆盖格局而最终得以表达. 更为关键的是,交互作用揭示了次要因子在特定情景下的效应放大机制. 例如,单独解释力相对较弱的DEM和NOx,一旦与植被因子叠加,其联合解释力便显著增强. 这种情景依赖性体现在,植被本底脆弱(NDVI值低)的高海拔区域,低温和缺氧等高程带来的负面胁迫效应会被显著放大;植被覆盖较好但受人类活动干扰的区域,氮氧化物(NOx)对生态系统的影响会变得更为复杂和关键[25].

(2)基于交互作用的解释力排序可识别出塑造EI空间格局的主导驱动组合  其中,“植被(NDVI/FVC)⋂降水(Pre)”“植被(NDVI/FVC)⋂海拔(DEM)”以及“植被(NDVI/FVC)⋂氮氧化物(NOx)”等组合表现出最强的解释力. 强交互组合进一步证实了青藏高原的生态环境状况是由自然本底条件(水、热和地形)与植被响应这两大系统的耦合所决定,同时在局部地区,以污染物为表征的人为压力已经成为不可忽视的关键调制变量. 因此,对该区域生态环境的评估与管理,必须从单一因子分析转向对这些关键协同作用的综合考量.

4 讨论

本研究揭示了青藏高原生态环境格局的形成,是一个由自然地理梯度奠定宏观背景,并由多因子协同作用塑造局部异质性的复杂过程. 其根本驱动机制体现为以植被为核心,气候和地表等多重要素相互耦合、显著增强的复杂网络.

4.1 宏观梯度与局部异质性耦合的生态环境空间格局

青藏高原的生态环境状况(EI)在空间上并非均质,而是呈现出宏观梯度与局部异质性并存的复杂格局. 研究结果清晰地展示了由东南向西北递减的宏观梯度,这与区域内巨大的自然地理环境差异存在直接因果联系. 高原东南部边缘区,受惠于季风带来的充沛水热[36,38],生态系统服务功能强大,EI整体较高. 与之形成鲜明对比的是,广袤的高原腹地与西北部,在严酷高寒和极度干旱的气候胁迫下[38],生态系统结构单一且稳定性差,EI整体偏低.

更为重要的是,本研究在宏观梯度之上,识别出了由人类活动主导的局部异质性. 通过横向对比发现,不同区域的生态压力源存在本质差异:在生态本底较好的东南部,压力主要表现为人口经济活动带来的附加型环境负荷(污染物排放);而在生态本底脆弱的西北部,压力则更多体现为系统固有的脆弱性(土地高胁迫和低承载力). 这种格局的分异,印证了本研究构建的EI评价体系能够有效地区分由自然梯度决定的背景质量与由人类活动造成的局部压力[38].

4.2 “主导⁃协同”的生态环境双重驱动机制

青藏高原的生态环境演变并非由单因子主导,而是遵循核心因子主导与多因子协同增强的复杂运行规律.

4.2.1 植被、降水与土地利用的核心主导作用

因子探测结果(表 5)定量辨识出植被、降水和土地利用是驱动EI空间分异的3个核心因子,且三者扮演的角色各不相同. 植被(FVC,q=0.87;NDVI,q=0.86)作为最敏感的生物指示器,是生态系统健康状况的最终体现. 降水(Pre,q=0.66)则是决定这一格局的首要限制性自然因子,其空间梯度深刻地塑造了植被与生态系统的宏观分布. 而土地利用(Land,q=0.70)作为一个综合性因子,同时耦合了自然禀赋(地形、土壤决定的利用潜力)与人类扰动(农牧业、城镇化),成为连接自然过程与社会经济活动的关键纽带. 三者之间明确的解释力差异与角色分工,为理解高原生态系统的主导链路提供了清晰的证据.

4.2.2 多驱动因子的综合效应与交互放大机制

本研究的科学性不仅在于识别主导因子,更在于揭示它们之间复杂的非线性关系. 交互探测(表 6)与SHAP分析的结论相互印证,为驱动机制的分析提供了极高的可靠性与鲁棒性. 地理探测器从宏观空间关联上证实了普适的协同增强规律:任意两因子的交互作用对EI的解释力均显著超越其单因子效应. 这证明,高原生态系统的演变是多重因素耦合叠加的非线性结果,而非简单的线性相加.

在此基础上,SHAP分析从微观归因层面进一步揭示了这一增强效应的核心机制——植被的中介与放大作用. 植被在因子交互网络中扮演了核心中介变量,其他因子(无论是自然背景还是人为压力)的生态效应多需通过作用于植被才能传导至整个系统. 植被作为“效应放大器”,使得独立解释力较弱的因子(DEM和NOx),在与植被的交互中其边际效应被显著放大,成为在特定情景下不可忽视的关键驱动力. 深刻揭示了高原生态系统对多重压力响应的非线性与情景依赖性特征.

5 结论

(1)青藏高原生态环境空间格局由宏观地理梯度与局域胁迫耦合决定.  青藏高原生态环境整体上呈现由东南向西北递减的地理特征,与水热条件主导的自然分区高度一致,构成区域生态质量的基础格局. 在此宏观背景之下,人类活动的胁迫效应具有明显空间差异,东南部生态禀赋优越的区域面临的主要风险源自社会经济活动带来的局部环境压力,而在西北部生态脆弱区,核心胁迫则更多表现为土地退化与高寒环境等内生性系统退化压力.

(2)植被是该生态系统的核心指示器与调节枢纽.  植被覆盖度与生长活力(FVC和NDVI)可有效反映生态系统健康状况,同时在气候、人类活动等外部驱动因子的作用传导中发挥关键调节作用. 因此,植被的动态监管是理解并调控生态系统响应机制的重要路径.

(3)非线性协同增强是青藏高原生态系统响应外部驱动的基本机制.  青藏高原生态系统对多重驱动力的响应表现出非线性协同增强特征,而非单因子效应的线性叠加. 交互探测明确证实,不同驱动因子组合普遍存在协同增强作用,揭示了生态退化与恢复过程均源于多因子耦合作用下的非线性机制. 该机制也意味着,即便贡献度较低的因子,在与主导因子(尤其是植被)交互时,其影响也被放大,成为影响局部生态演替方向的关键控制因子.

参考文献
[1] 董贵华, 王业耀, 于洋, 等. "十三五"以来我国生态质量状况时空变化分析[J]. 中国环境监测, 2023, 39(1): 1-9.
Dong G H, Wang Y Y, Yu Y, et al. Analysis on spatio-temporal changes of ecological quality status in China since the 13th five year plan period[J]. Environmental Monitoring in China, 2023, 39(1): 1-9.
[2] 李艳翠, 袁金国, 刘博涵, 等. 滹沱河流域生态环境动态遥感评价[J]. 环境科学, 2024, 45(5): 2757-2766.
Li Y C, Yuan J G, Liu B H, et al. Ecological environment dynamical evaluation of Hutuo River Basin using remote sensing[J]. Environmental Science, 2024, 45(5): 2757-2766. DOI:10.13227/j.hjkx.202305251
[3] HJ 192-2015, 生态环境状况评价技术规范[S].
[4] 郭浩, 董磊, 邬伦, 等. 空间异质性建模方法[J]. 地理学报, 2025, 80(3): 567-585.
Guo H, Dong L, Wu L, et al. Tackling spatial heterogeneity in geographical analysis: an overview[J]. Acta Geographica Sinica, 2025, 80(3): 567-585.
[5] 王劲峰, 徐成东. 地理探测器: 原理与展望[J]. 地理学报, 2017, 72(1): 116-134.
Wang J F, Xu C D. Geodetector: principle and prospective[J]. Acta Geographica Sinica, 2017, 72(1): 116-134.
[6] 宋媛, 石惠春, 谢敏慧, 等. 2000—2017年甘肃省生态环境质量时空演变格局及其影响因素[J]. 生态学杂志, 2019, 38(12): 3800-3808.
Song Y, Shi H C, Xie M H, et al. Spatiotemporal evolution pattern and influencing factors of eco-environmental quality in Gansu from 2000 to 2017[J]. Chinese Journal of Ecology, 2019, 38(12): 3800-3808.
[7] 胡克宏, 张震, 郜敏, 等. 中国丝绸之路经济带沿线植被覆盖变化及自然影响因素分析[J]. 农业工程学报, 2020, 36(17): 149-157.
Hu K H, Zhang Z, Gao M, et al. Variations in vegetation cover and natural factors of provinces in China along silk road economic belt during 2000-2018[J]. Transactions of the Chinese Society of Agricultural Engineering, 2020, 36(17): 149-157.
[8] 宋沛林, 解吉波, 杨腾飞, 等. 基于生态敏感性的区域环境时空变化及其驱动因素探究——以定安县为例[J]. 测绘通报, 2023(7): 18-24.
Song P L, Xie J B, Yang T F, et al. Exploring the spatio-temporal variation of the regional environment and driving factors based on ecological sensitivity: taking Ding'an County as an example[J]. Bulletin of Surveying and Mapping, 2023(7): 18-24.
[9] Lundberg S M, Lee S I. A unified approach to interpreting model predictions[A]. In: Proceedings of the 31st International Conference on Neural Information Processing Systems[C]. Long Beach, CA, USA: Curran Associates Inc., 2017. 4768-4777.
[10] Bae S, Son B, Sung T, et al. Advancing hourly gross primary productivity mapping over East Asia using Himawari-8 AHI and artificial intelligence: unveiling the impact of aerosol-induced radiation dynamics[J]. Remote Sensing of Environment, 2025, 323. DOI:10.1016/j.rse.2025.114735
[11] Song Y Z, Wang J F, Ge Y, et al. An optimal parameters-based geographical detector model enhances geographic characteristics of explanatory variables for spatial heterogeneity analysis: cases with different types of spatial data[J]. GIScience & Remote Sensing, 2020, 57(5): 593-610.
[12] Werther M, Odermatt D, Simis S G H, et al. Characterising retrieval uncertainty of chlorophyll-a algorithms in oligotrophic and mesotrophic lakes and reservoirs[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2022, 190: 279-300. DOI:10.1016/j.isprsjprs.2022.06.015
[13] Li C Q, Guo S C, Huang Q, et al. Dynamic monitoring of fine-grained ecological vulnerability in dryland urban agglomeration integrating novel remote sensing index and explainable machine learning[J]. GIScience & Remote Sensing, 2025, 62(1). DOI:10.1080/15481603.2025.2528302
[14] Jin W, Li H Z, Wang J Z, et al. Continuous remote sensing ecological index (CRSEI): a novel approach for multitemporal monitoring of eco-environmental changes on large scale[J]. Ecological Indicators, 2023, 154. DOI:10.1016/j.ecolind.2023.110739
[15] Gu Z L, Gong J, Wang Y. Construction and evaluation of ecological networks among natural protected areas based on "quality-structure-function": a case study of the Qinghai-Tibet area[J]. Ecological Indicators, 2023, 151. DOI:10.1016/j.ecolind.2023.110228
[16] 刘慧文, 刘欢, 胡鹏, 等. 基于可解释机器学习的青藏高原草地物候变化多因素影响分析[J]. 环境科学, 2024, 45(6): 3375-3388.
Liu H W, Liu H, Hu P, et al. Multi-factor impact analysis of grassland phenology changes on the Qinghai-Xizang Plateau based on interpretable machine learning[J]. Environmental Science, 2024, 45(6): 3375-3388. DOI:10.13227/j.hjkx.202306185
[17] Lu Y J, Zhang Z, Kong Y R, et al. Integration of optical, SAR and DEM data for automated detection of debris-covered glaciers over the western Nyainqentanglha using a random forest classifier[J]. Cold Regions Science and Technology, 2022, 193. DOI:10.1016/j.coldregions.2021.103421
[18] 王晓峰, 郑媛元, 孙泽冲, 等. 青藏高原自然保护区生态脆弱性演化与模拟[J]. 环境科学, 2025, 46(3): 1633-1644.
Wang X F, Zheng Y Y, Sun Z C, et al. Evolution and simulation of ecological vulnerability in Qinghai-Xizang Plateau nature reserve[J]. Environmental Science, 2025, 46(3): 1633-1644. DOI:10.13227/j.hjkx.202403096
[19] 李平星, 高晨真, 罗艳华, 等. 青藏高原生态资产空间差异及其价值化潜力[J]. 生态学报, 2024, 44(5): 1808-1821.
Li P X, Gao C Z, Luo Y H, et al. Spatial differences and value conversion potential of ecological assets on the Qinghai-Tibet Plateau, China[J]. Acta Ecologica Sinica, 2024, 44(5): 1808-1821.
[20] Ding B J, Feng L, Ba S, et al. Temperature drives elevational diversity patterns of different types of organisms in Qinghai-Tibetan Plateau wetlands[J]. iScience, 2023, 26(8). DOI:10.1016/j.isci.2023.107252
[21] Li Y G, Liu F, Zhou Y D, et al. Large-scale geographic patterns and environmental and anthropogenic drivers of wetland plant diversity in the Qinghai-Tibet Plateau[J]. BMC Ecology and Evolution, 2024, 24(1): 74. DOI:10.1186/s12862-024-02263-w
[22] 徐祥德, 赵天良, Chungu L, 等. 青藏高原大气水分循环特征[J]. 气象学报, 2014, 72(6): 1079-1095.
Xu X D, Zhao T L, Chungu L, et al. Characteristics of the water cycle in the atmosphere over the Tibetan Plateau[J]. Acta Meteorologica Sinica, 2014, 72(6): 1079-1095.
[23] Yu H B, Yang M, Lu Z X, et al. A phylogenetic approach identifies patterns of beta diversity and floristic subregions of the Qinghai-Tibet Plateau[J]. Plant Diversity, 2024, 46(1): 59-69. DOI:10.1016/j.pld.2023.07.006
[24] Liu M X, Xu L, Mu R L, et al. Plant community assembly of alpine meadow at different altitudes in northeast Qinghai-Tibet Plateau[J]. Ecosphere, 2023, 14(1): e4354. DOI:10.1002/ecs2.4354
[25] Duan H C, Huang B Y, Liu S L, et al. Impact of extreme climate indices on vegetation dynamics in the Qinghai–Tibet Plateau: a comprehensive analysis utilizing long-term dataset[J]. ISPRS International Journal of Geo-Information, 2024, 13(12): 457. DOI:10.3390/ijgi13120457
[26] SL 190-2007, 土壤侵蚀分类分级标准[S].
[27] Geng X, Wang X M, Fang H L, et al. Vegetation coverage of desert ecosystems in the Qinghai-Tibet Plateau is underestimated[J]. Ecological Indicators, 2022, 137. DOI:10.1016/j.ecolind.2022.108780
[28] Gröschler K C, Martens T, Schrautzer J, et al. Data-driven identification of high-nature value grasslands using harmonized Landsat sentinel-2 time series data[J]. Remote Sensing Applications: Society and Environment, 2025, 37. DOI:10.1016/j.rsase.2024.101427
[29] 杨昌兰. 基于深度学习及SHAP解释的城市空间扩展识别与驱动分析[D]. 武汉: 武汉大学, 2022.
Yang C L. Identification and driving analysis of urban spatial expansion through intergrating deep learning models and SHapley Additive exPlanations[D]. Wuhan: Wuhan University, 2022.
[30] 周亚男. 基于多源遥感数据和机器学习算法的土壤质地类型预测研究[D]. 重庆: 西南大学, 2023.
Zhou Y N. Prediction of soil texture class based on multi-source remote sensing data and machine learning algorithms[D]. Chongqing: Southwest University, 2023.
[31] Cao F, Ge Y, Wang J F. Optimal discretization for geographical detectors-based risk assessment[J]. GIScience & Remote Sensing, 2013, 50(1): 78-92.
[32] 向兴华. 树集成预测模型规则提取应用研究[D]. 湘潭: 湖南科技大学, 2021.
Xiang X H. Application research of tree ensemble prediction model rule extraction[D]. Xiangtan: Hunan University of Science and Technology, 2021.
[33] 旦增桑姆, 单增尼玛. 高原气象变化对生态环境的影响分析与解读[J]. 资源节约与环保, 2019(2): 146.
[34] 牛振国, 景雨航, 张东启, 等. 气候变化背景下青藏高原湿地生态系统响应特征: 回顾与展望[J]. 气候变化研究进展, 2024, 20(5): 509-518.
Niu Z G, Jing Y H, Zhang D Q, et al. An overview and the outlook for wetland ecosystems in the Qinghai-Tibetan Plateau under climate change[J]. Climate Change Research, 2024, 20(5): 509-518.
[35] Härkönen L H, Lepistö A, Sarkkola S, et al. Reviewing peatland forestry: implications and mitigation measures for freshwater ecosystem browning[J]. Forest Ecology and Management, 2023, 531. DOI:10.1016/j.foreco.2023.120776
[36] Liu X, Hu L T, Sun K N, et al. Improved understanding of groundwater storage changes under the influence of river basin governance in Northwestern China using GRACE data[J]. Remote Sensing, 2021, 13(14). DOI:10.3390/rs13142672
[37] 陈卫东, 取宗, 李晓童. 青藏高原植被指数时空特征及其影响因素[J]. 经济地理, 2025, 45(1): 193-203.
Chen W D, Qu Z, Li X T. Spatiotemporal evolution of vegetation index and its influencing factors in the Qinghai-Xizang Plateau[J]. Economic Geography, 2025, 45(1): 193-203.
[38] 张云霞, 汪仕美, 李焱, 等. 2000—2020年青藏高原生态质量时空变化[J]. 生态学杂志, 2023, 42(6): 1464-1473.
Zhang Y X, Wang S M, Li Y, et al. Spatiotemporal patterns of ecological quality across the Qinghai-Tibet Plateau during 2000-2020[J]. Chinese Journal of Ecology, 2023, 42(6): 1464-1473.