**基于kNDVI与可解释机器学习XGBoost-SHAP的山西省植被覆盖时空演变、阶段性驱动机制及非线性响应阈值研究**

《Frontiers in Environmental Science》：Spatial and Temporal characteristics and driving force analysis of vegetation cover change in Shanxi Province, China based on kNDVI and XGBoost-SHAP model

　　

摘要：研究背景：系统阐明植被覆盖的时空演变特征并揭示其对气候波动和人类活动的响应机制，对于推进山西省"双碳"战略（碳达峰和碳中和）和优化国土空间开发保护格局具有重要现实意义。方法：基于核归一化植被指数（kNDVI，kernel Normalized Difference Vegetation Index）数据集以及多源气象、地形和社会经济数据集，本研究整合Theil-Sen斜率估计、Mann-Kendall显著性检验、Hurst指数（赫斯特指数，用于量化时间序列趋势持久性的指标）、标准差椭圆（SDE，Standard Deviational Ellipse）和重心迁移模型，系统表征了2000至2024年植被覆盖的时空格局并预测其未来演变的潜在持续性。同时，构建了可解释的XGBoost（极端梯度提升）-SHAP（Shapley加性解释）机器学习框架，量化了七个驱动因子——海拔、坡度、坡向、年均温度、年均降水量、人口密度和夜间灯光强度——的独立贡献、非线性边际响应及其时间演化特征。结果：（1）时间上，全省多年平均kNDVI约为0.1699，呈极显著波动上升趋势，年增长率为0.0041。年度kNDVI范围从2001年记录的0.1132到2024年的峰值0.2227。（2）空间上，植被覆盖呈现显著的"北低南高"分异格局。较高kNDVI值主要分布在太行山、吕梁山、中条山和太岳山等森林覆盖山地，而较低值分布在晋北沙丘丘陵区和汾河谷地城市群。研究区98.02%的区域观测到植被恢复，其中约92.18%呈极显著改善。（3）2000–2024年间，植被空间分布的主轴保持稳定的东北-西南方向，标准差椭圆呈现适度向外扩张，植被重心总体向西北方向迁移。（4）Hurst指数表明，区域98.03%的植被改善趋势将持续，而零星的工矿和城镇用地斑块（占1.97%）面临潜在的植被退化风险。（5）区域植被驱动系统经历了三个演化阶段：单一降水主导的自然驱动阶段（2000–2005年）、气候-人类耦合过渡阶段（2005–2010年）和多因子协同平衡阶段（2010–2020年）。降水和温度作为核心气候驱动因子，海拔在很大程度上决定了植被分布的垂直分异基线。人口密度和夜间灯光强度对植被生长产生持续性抑制效应。经验推导的主要预测因子初步阈值如下：海拔约1,200 m、坡度约8°、年均温度约7.5 °C、年均降水量约500 mm、人口密度500人/km²、夜间灯光强度约5。这些值可作为区分生态保育区与人类活动干扰区的参考边界。结论：本研究识别了黄土高原温带过渡带山西省植被动态的跨期分异与非线性耦合规律。提出了面向气候适应、地形分区和人为压力调控的差异化生态治理策略，为中国黄土高原生态安全屏障建设提供了观测科学支持。

**论文解读**

**一、研究背景与科学问题**

全球气候变暖与加速的人类土地利用活动正在重塑陆地植被生态系统的结构与功能，在半干旱生态过渡带引发了植被退化、土壤侵蚀加剧和生物多样性丧失等普遍性生态胁迫。植被调控着陆地碳循环、水循环以及地表能量平衡，其长期动态直接综合反映了气候背景条件和人为干扰信息，因而被视作全球生态环境监测的核心指示指标。

山西省地处黄土高原东部，属于暖温带落叶阔叶林与温带草原的生态过渡带（ecotone），为黄河中游水土保持、风力侵蚀防治和区域气候调节提供关键生态系统服务。该区域地形异质性强烈，形成了显著的南北水热梯度，煤矿开采产业和城市扩张带来的人为干扰与退耕还林工程（Grain-for-Green Project）、三北防护林工程（Three-North Shelter Forest Programme）等重大国家生态修复举措并存。这种气候变异性、地形异质性、人类开发与生态修复的交织共存，使山西成为揭示半干旱区气候-人类活动耦合作用下植被响应的天然实验场。

然而，既有研究面临三方面不足：第一，传统归一化植被指数（NDVI，Normalized Difference Vegetation Index）在山地森林区存在光谱饱和与混合像元偏差，可能导致植被覆盖定量评估的系统性误差；第二，多数驱动力研究依赖偏相关或偏最小二乘回归等线性统计方法，难以捕捉植被状态与环境协变量之间复杂的非线性关联；第三，缺乏多时期驱动力贡献的阶段对比分析，且尚未从经验层面定量识别多因子非线性响应阈值。围绕上述研究空白，本文聚焦三个核心问题：kNDVI植被条件在25年间如何时空演变？植被变化趋势的持续性如何？驱动植被动态的主导自然与人为因子是什么，关键预测因子是否表现出可识别的非线性响应阈值？

**二、研究目标与总体思路**

本研究以山西省为研究区，基于MOD13Q1遥感产品生成NDVI时间序列并重建kNDVI，联合应用Theil-Sen斜率趋势分析、Mann-Kendall显著性检验和Hurst指数法，全面表征植被年际波动、空间分异、重心迁移及未来演变可持续性；在此基础上，构建XGBoost-SHAP可解释机器学习框架，定量识别影响植被覆盖变化的主导因子，阐明各因子影响强度、贡献机制与边际效应。该研究旨在为黄土高原山西段差异化生态修复规划与国土空间治理提供经验证据。

**三、主要技术方法**

数据方面，植被数据采用国家青藏高原科学数据中心发布的MODIS MOD13Q1 NDVI数据集（2000–2024年，250 m分辨率），经Savitzky–Golay（SG）滤波和最大值合成法（MVC，Maximum Value Composite）生成月尺度时序；气象数据来自彭守璋团队研制的1 km月尺度气温和降水格网产品；地形数据为30 m分辨率GDEM数字高程模型（DEM，Digital Elevation Model），并派生出坡度和坡向；社会经济数据采用LandScan人口密度数据（250 m）和Harvard Dataverse夜间灯光数据（500 m）。全部数据经统一空间坐标转换和最近邻重采样，实现多源数据时空精准匹配。

方法上，首先依据kNDVI = tanh(NDVI²)公式重建kNDVI序列；其次采用Theil-Sen中位数斜率估计联合Mann-Kendall显著性检验量化kNDVI年际趋势及其显著性；运用Hurst指数评估趋势持续性，并叠合趋势显著性分类未来植被变化轨迹；利用标准差椭圆（SDE）和重心迁移模型刻画植被分布空间方向性、离散程度及重心迁移路径；最后构建XGBoost-SHAP可解释机器学习框架，以海拔、坡度、坡向、年均温度、年均降水量、人口密度和夜间灯光强度七项因子为预测变量，以5年为间隔选取2000、2005、2010、2015和2020年五个时点切片建模，通过SHAP均值柱状图、蜂群散点图和依赖曲线图（SHAP dependence plot），系统揭示因子贡献时序演化、非线性响应关系和关键阈值。所有样本按70%与30%随机划分为训练集和独立测试集，以R²和RMSE评估模型拟合与泛化性能。

**四、主要研究结果**

**4.1 植被覆盖时空演变特征（2000–2024年）**

时间趋势方面，全省多年平均kNDVI为0.1699，呈极显著波动上升趋势（p < 0.01），年增长率为0.0041 yr^-1，最低值0.1132出现在2001年，最高值0.2227出现在2024年。2010年后年际kNDVI波动幅度收窄，植被生态系统稳定性持续增强。

空间分异方面，多年度平均kNDVI空间分布呈显著"北低南高"带状分异。高值区集中于太行山、吕梁山、中条山和太岳山等森林覆盖山地；中等值过渡区分布于山间盆地和黄土塬面的农牧镶嵌景观；低值区集中于晋北风沙黄土丘陵区（大同、朔州、忻州）和汾河谷地城市密集带。Theil-Sen斜率均值为0.0043，98.02%的区域植被呈增加趋势，其中92.18%呈极显著改善，植被退化区零星分布且面积占比不足1%。

**4.2 未来趋势持续性**

Hurst指数（H）范围为0.2761–0.7691，H > 0.5的持续型区域占98.03%，其中极显著持续改善占92.17%，构成主导趋势类型。退化型区域仅占1.97%，以持续性退化为主（1.87%），工矿和城镇建成区是退化风险热点。非持续性改善和稳定状态合计仅占0.44%，表明当前植被恢复态势具有高度长期稳定性。

**4.3 空间分布格局迁移**

SDE分析显示，植被空间分布主轴始终维持稳定的东北-西南方向，椭圆面积从2000年的77,998.52 km²波动扩大至2024年的80,540.66 km²，旋转角θ从11.25°微调至11.15°。半短轴（东西向）持续扩张，半长轴（南北向）基本稳定，表明植被空间覆盖向外扩张且以东西向为主。重心迁移呈阶段性振荡：2000–2005年向东北迁移6.56 km，2005–2010年向西迁移3.55 km，2010–2015年向南迁移4.42 km，2015–2020年向西北迁移6.63 km，2020–2024年向北迁移10.59 km。整体而言，重心长期呈明显西北向迁移趋势——经度从112.33°西移至112.27°，纬度从37.23°北移至37.31°。

**4.4 XGBoost-SHAP驱动机制分析**

五期SHAP均值贡献显示，区域植被驱动系统经历了三个清晰演化阶段：

**阶段一（2000–2005年）**：降水单一自然驱动阶段。降水（SHAP均值0.0331–0.0350）和温度（0.0323→0.0114）主导植被空间分异，人为因子贡献微弱。

**阶段二（2005–2010年）**：气候-人类耦合过渡阶段。温度取代降水成为首要驱动因子（SHAP均值0.0406），海拔跃升至第二位（0.0337），人口密度贡献持续攀升（0.0186），退耕还林工程增强了对海拔和坡度的正面植被效应。

**阶段三（2010–2020年）**：多因子协同平衡阶段。温度保持主导地位，海拔、降水和人口密度四项核心因子SHAP贡献差距大幅缩小，气候、地形与人为活动协同调控植被空间演化。2020年排名为：温度（0.0291）> 海拔（0.0277）> 降水（0.0257）> 人口密度（0.0235）> 坡度（0.0108）> 夜间灯光（0.0060）> 坡向（0.0044）。

SHAP依赖曲线定量识别了各因子非线性响应阈值：海拔约1,200 m为盆地压制与山地促进的分界，呈U型响应；坡度约8°为人类占用与生态保育的地形分界线；年均温度约7.5 °C以下热量不足抑制植被，7.5–12 °C为促进区间，高于12 °C边际增益递减；年降水量500 mm为水分限制阈值，超过500 mm后植被响应陡增，呈S型单调递增，未出现降水饱和抑制拐点；人口密度500人/km²为人为干扰突变区间，超过该阈值后SHAP值持续为负且不可逆；夜间灯光强度约5为轻度干扰与严重退化的分界值。坡向全值域SHAP波动幅度不足0.01，属宏观非显著性因子，仅对山地微生境产生局部调节。

**五、讨论与结论**

**讨论方面**，研究人员指出，山西省25年间广泛的植被恢复与黄土高原整体生态改善趋势一致，其核心驱动力来自退耕还林、天然林保护工程和三北防护林等数十年持续生态修复政策的叠加效应。基于XGBoost-SHAP识别的非线性阈值和阶段性驱动规律，提出了差异化生态治理框架：以海拔1,200 m和坡度8°为地形分界，划分核心生态保育区、生态-农业过渡区和人为密集干扰区，实施三级地形分区管理；以年均降水500 mm为界，对干旱北部推行封育保护与灌草修复为主、避免盲目大规模造林的被动修复策略，对中高山湿润林区实施水源涵养林建设；对超过人口密度和夜间灯光阈值的汾河城市群和煤矿区，采取紧凑发展模式、矿区全周期生态治理和生态缓冲带建设等差异化人为扰动调控措施。研究同时指出，Hurst持续性推断以气候和人类活动背景稳定为前提，无法反映未来政策调整或极端气候事件引发的突变；SHAP阈值来源于全省尺度全像元建模，为省级平均响应，尚未开展子区域独立建模以获取地域特异性阈值；多源数据空间分辨率不匹配（kNDVI 250 m与气象数据1 km）可能引入偏差；阈值拐点属观测经验值，未进行显著性检验。

**研究结论**：（1）25年间全省kNDVI呈极显著波动上升，年增长率为0.00406，2010年后生态系统稳定性逐步增强，表明数十年级联式生态修复工程显著提升了陆地植被系统的功能韧性。（2）植被覆盖呈南北空间异质性，北低南高，山地森林为高覆盖核心热点，低值区集中于晋北风沙丘陵和汾河谷地人为扰动走廊带，98.02%区域植被改善，植被空间分布主轴稳定为东北-西南方向，重心持续向西北迁移。（3）在现有气候和生态治理条件稳定的前提下，98.03%区域植被改善趋势将维持长期持续性，仅1.97%的零散工矿城镇斑块面临退化风险。（4）植被驱动系统经历三个阶段演替——降水主导期（2000–2005年）、气候-人类耦合过渡期（2005–2010年）和多因子协同平衡期（2010–2020年），人为活动从次级背景干扰演变为与气候地形等同重要的核心驱动力。（5）SHAP依赖曲线揭示了主要因子的分段非线性响应及经验临界阈值：海拔1,200 m、坡度8°、年均温7.5 °C、年降水500 mm、人口密度500人/km²和夜间灯光强度5，这些阈值可作为黄土高原山西段生态安全屏障建设的观测经验参考基准，支撑从定性经验管理向阈值导向的差异化空间精准治理范式转变。