蒙古高原植被春季物候的时空动态及其驱动机制：基于XGBoost和SHAP的洞察 张宇、 程浩、 李富佳、 陈莉

《Land》：Spatiotemporal Dynamics and Driving Mechanisms of Vegetation Spring Phenology on the Mongolian Plateau: Insights from XGBoost and SHAP Yu Zhang, Hao Cheng, Fujia Li and Li Chen

　　**摘要**
干旱地区的植被春季物候对气候和人为压力非常敏感，然而不同植被类型在生长季节开始时间（SOS）上的非线性响应尚未得到充分量化。本研究从2001年至2020年获取了蒙古高原（MP）稳定植被区域的MODIS中分辨率成像光谱辐射计（MODIS）归一化差异植被指数（NDVI）时间序列，并应用极端梯度提升（XGBoost）模型结合Shapley加性解释（SHAP）分析法，研究了六个环境因素——降水量、温度、风速、牲畜密度、人口密度和海拔高度——对森林、灌木丛和草地的影响。研究结果显示，生长季节开始时间在空间上存在明显异质性：北部森林和灌木丛的开始时间较早，而南部干旱草地的开始时间则较晚。森林和灌木丛的开始时间平均每十年分别提前了6.8天和6.4天，而草地则没有显著变化。尽管温度在所有植被类型的SOS变化中起主导作用，但其他因素的相对重要性也有所不同：风速对森林的影响尤为显著，而降水量和海拔高度对草地和灌木丛尤为重要。SHAP分析揭示了强烈的非线性和阈值效应，例如草地存在一个350毫米降水量阈值，超过该阈值后SOS响应会发生显著变化。这些结果表明，干旱地区物候调节具有植被特异性和非线性特征，因此物候预测和生态系统监测应明确考虑植被类型和基于阈值的气候响应。

**1. 引言**
植被物候是生态系统响应气候变化的关键指标。在物候指标中，生长季节开始时间（SOS）被广泛用于反映植被活动和生态系统生产力[1]。SOS的变化会直接影响碳循环、水分平衡和能量交换，因此了解SOS的动态变化对于评估全球气候变化下的生态系统响应至关重要[2,3]。干旱和半干旱地区由于水资源有限和生态系统脆弱，对环境变化尤为敏感[4,5]。蒙古高原是世界上最大的干旱地区之一，以草地、灌木丛和森林生态系统为主[6]。近年来，该地区经历了显著的气候变化和人类活动的增加，如放牧和人口扩张，这些因素可能强烈影响植被物候[7,8]。然而，与湿润地区相比，关于干旱生态系统物候动态的研究仍然有限，尤其是在森林、灌木丛和草地生态系统对多种环境因素的不同响应方面[9]。以往的研究主要关注气候变量（如温度和降水量）与SOS之间的关系[10]，但仍存在一些局限性：首先，许多研究基于假设简单或单调响应的统计关系，可能无法捕捉到复杂的非线性行为[11]；其次，环境因素的阈值效应往往没有得到充分量化，尽管在干旱生态系统中植被响应很少是线性的[12,13]；第三，人类活动（如牲畜密度和人口压力）的影响相对较少受到关注[14,15]；另外，尽管森林、灌木丛和草地具有不同的生态特征，但它们之间的差异往往没有得到充分考虑[16,17]。基于这些不足，本研究旨在探讨蒙古高原以下问题：(1) 2001年至2020年间，SOS是否表现出显著的时间变化和空间异质性？(2) 控制SOS动态的主要驱动因素是什么，气候因素与人为因素的相对贡献如何？(3) SOS与其驱动因素之间是否存在非线性关系和阈值效应，这些关系是否因植被类型而异？为了解答这些问题，本研究首先从长期NDVI时间序列中提取SOS，以尽量减少土地利用变化的影响，然后将多种驱动因素（包括气候变量、地形条件和人类活动）整合到XGBoost模型中，并应用SHAP方法量化各因素的相对重要性及其非线性效应和阈值。这一框架有助于全面理解干旱生态系统中的植被物候动态及其相关影响因素。

**2. 材料与方法**
**2.1. 研究区域**
研究区域位于蒙古高原（37°22′–53°20′ N and 87°43′–126°04′ E），涵盖蒙古和中国内蒙古的部分地区（图1a）。该地区属于大陆性气候，冬季寒冷干燥，夏季温暖湿润。年平均气温介于-5°C至8°C之间，年降水量从北部的约400毫米逐渐减少到南部的不足100毫米，呈东北-西南方向递减。海拔高度从约500米到超过3000米不等（图1b）。西部和北部以阿尔泰山、杭爱山和肯特山为主，东部和南部则为低洼平原和高原[1]。图1展示了研究区域的位置（a）、蒙古高原的海拔高度（b）、基于MODIS土地覆盖产品（MCD12Q1, IGBP）得出的植被覆盖类型（c）以及植被物候监测区域的范围（d）。自然植被主要由森林、灌木丛和草地组成（图1c）。森林主要分布在北部山脉，以落叶松（Larix sibirica Ledeb）和白桦（Betula pendula subsp. mandshurica (Regel) Ashburner & McAll）为主；灌木丛则分布在森林和草地之间的过渡地带[18]；草地（包括草地草原、典型草原和沙漠草原）占据了研究区域的最大比例[19]。为排除非植被区域对物候检测的影响，分析仅包括那些土地覆盖类型在整个研究期间保持稳定且NDVI值持续高于0.1的区域（图1d）。

**2.2. 数据集与处理**
为便于后续研究，所有数据集均使用ArcGIS 10.2进行处理，并重采样至1000米的空间分辨率。
- **2.2.1. 植被覆盖类型数据**
我们使用了2001年和2020年的MODIS MCD12Q1 V6数据（https://modis.gsfc.nasa.gov/，访问日期：2026年3月25日），这些数据用于识别稳定植被区域。具体而言，只有2001年和2020年都被分类为植被的像素被视为稳定区域，排除在两个时间点之间存在任何土地覆盖变化的像素。这些稳定区域被确定为本研究的研究对象。该数据集定义了17种主要土地覆盖类型，其中1-5类归类为森林，6-9类归类为灌木丛，10类归类为草地。另外，11-17类（代表非植被覆盖）被排除，以确保不会干扰植被物候的监测（表1）。
- **2.2.2. 植被物候数据**
植被物候数据来自MODIS MOD13Q1 Collection 6产品，该产品提供了250米空间分辨率的16天复合NDVI（https://modis.gsfc.nasa.gov/，访问日期：2026年3月25日）。研究时间范围为2001年至2020年，所有可用图像均经过筛选和预处理。NDVI值低于0.1的像素被排除，以减少非植被表面的噪声。SOS是通过基于阈值的方法提取的[20]：对于每个像素和每年，计算最小和最大NDVI值，然后定义一个动态阈值，该阈值为主播NDVI值加上年度NDVI振幅的固定百分比（20%）[21]，并将SOS定义为NDVI首次超过该阈值的那一天（图2）。为确保数据质量，仅保留研究期间土地覆盖稳定的像素。
- **2.2.3. 驱动因素数据**
本研究选择了六个驱动因素（表2）：三个气候变量（降水量、温度和风速）、两个人类活动变量（人口密度和牲畜密度）以及一个地形变量（数字高程模型）。所有驱动因素数据集均重采样至一致的空间分辨率，并投影到同一坐标系中进行分析。

**2.3. 方法**
**2.3.1. SOS提取**
SOS数据来自2001年至2020年的MODIS NDVI时间序列（MOD13Q1），我们采用基于阈值的方法结合平滑技术来减少噪声。首先应用Savitzky–Golay滤波器平滑原始NDVI时间序列，该滤波器可减少云层污染和大气变化引起的噪声，同时保留季节性模式[23]。然后使用动态阈值提取每个像素和每年的SOS[21]：对于给定年份，设最小NDVI为min NDVI，最大NDVI为max NDVI，年度NDVI振幅为ΔNDVI = max NDVI - min NDVI。阈值θ定义为ΔNDVI的固定百分比（20%），SOS则被定义为NDVI首次超过θ的那一天（图2）。

**2.3.2. 趋势分析**
2001年至2020年间SOS开始时间的时间趋势采用非参数Theil–Sen斜率估计器结合Mann–Kendall（MK）显著性检验进行逐像素分析[25,26,27]。这种组合方法在植被物候研究中得到广泛应用，因为它对异常值具有鲁棒性，且不需要数据正态性的假设[28]。Sen斜率（S）用于量化时间序列中的单调趋势强度。对于给定年份，令y_i表示前i年的SOS值，S = Σ(y_i - y_(i-1))/n，其中n为年份数。任意两个年份之间的斜率表示为：
S = (y_i - y_(i-1))/(y_n - y_1)
整体Sen斜率估计值为所有成对斜率的中位数[4]。正的S值表示SOS的延迟趋势（更晚的绿化），负的S值表示提前趋势（更早的绿化）。S的绝对值表示年均变化天数。MK检验用于评估检测到的趋势的统计显著性，检验统计量T计算为：
T = (S - μ) / (σ√n)，其中μ为均值，σ为标准差。根据预先建立的显著性阈值，将S的趋势分为七类（表3）。

**2.3.3. 多重共线性检验**
在构建XGBoost模型之前，进行了多重共线性检验，以避免解释变量之间的冗余信息。方差膨胀因子（VIF）用于量化变量之间的线性相关程度，较高的VIF表示较强的多重共线性，其倒数（1/VIF）表示每个变量贡献的唯一信息占比[29]。根据常见阈值，VIF大于10的变量被认为具有严重多重共线性并被排除。本研究所有变量的VIF值均低于10，表明多重共线性程度可接受。由于森林、灌木丛和草地生态系统具有不同的环境条件，分别进行了VIF分析，结果见表4（表4）。所有三种植被类型的平均VIF值均低于3，个别VIF值均小于6，说明所选驱动因素之间的多重共线性较低。因此，所有六个变量均被保留用于后续的XGBoost建模和SHAP分析。

**2.3.4. XGBoost模型**
使用XGBoost算法研究了SOS对多种环境因素的非线性响应。XGBoost是一种集成学习方法，通过顺序构建决策树序列，每棵新树纠正前一棵树的错误[30]。由于其能够处理复杂交互作用和异构输入数据，这种方法在生态学应用中越来越受欢迎[31]。基于VIF诊断结果，所有六个候选预测变量均被保留用于模型构建。这些预测变量包括三个气候变量（PRE、TEM和WIN）、两个人为指标（DENSITY和POP）以及一个地形描述符（DEM），SOS作为响应变量。为了考虑森林、灌木丛和草地生态系统的独特特征，我们分别开发了XGBoost模型。XGBoost算法最小化一个正则化目标函数，该函数结合了损失项和惩罚项：（9）其中测量了观测到的SOS（）与预测值（）之间的差异。该惩罚项是一个正则化函数，用于惩罚树木的数量和每棵树的深度。这种正则化可以减少过拟合并提高模型的泛化能力。与传统线性回归不同，XGBoost不对数据施加线性或正态性的假设，这使得它特别适合检测物候指标与环境驱动因素之间的非线性关系和阈值行为。对于每种植被类型，建模数据集包含了相应植被类别内的所有有效像素-年份样本，从而获得了足够大的样本量以进行模型训练。数据集被随机划分为训练子集（70%）和验证子集（30%），使用固定的随机种子（123）以确保可重复性。XGBoost模型使用梯度提升框架实现，并采用平方误差损失函数。关键超参数的设置旨在平衡模型复杂性和泛化能力，包括最大树深度为8，学习率为0.09。为了允许灵活的树构建，最小子节点权重被设置为10，并最小化了正则化项。

除了决定系数（R2）之外，还使用均方根误差（RMSE）和平均绝对误差（MAE）进一步评估了模型性能，以提供更全面的预测准确性评估。

2.3.5 SHAP分析
为了阐明各个预测因子的贡献，我们应用了基于合作博弈论的SHAP方法。在SHAP框架中，模型生成的预测被视为集体收益，每个解释变量被视为一个玩家，其边际贡献通过所有可能的特征组合进行评估[32]。
对于给定的变量，其Shapley值（）是其在所有可能变量子集上的边际贡献的加权平均值：（10）其中是所有解释变量的完整集合，是不包含变量 的任何子集，是基于子集 的模型预测。正的Shapley值意味着该变量使SOS预测时间延迟；负值则意味着它使SOS预测时间提前。

我们同时追求全局和局部的可解释性。全局上，通过平均所有观测值上的绝对Shapley值来评估变量重要性——较大的平均绝对Shapley值表示对SOS预测的更大影响。局部上，构建了ShapSH dependence图，以可视化在考虑所有其他协变量平均效应的情况下，预测的SOS如何沿着单个驱动因子的梯度变化。通过绘制特征值与其对应的Shapley值并用局部加权回归（LOESS）平滑处理来构建ShapSH dependence曲线，以捕捉潜在的非线性趋势。阈值被定义为平滑后的ShapSH曲线穿过零点或表现出明显拐点的位置，这表明了延迟和提前效应之间的转换。

通过结合XGBoost和SHAP的分析框架，我们实现了高预测准确性和透明的模型解释。这种方法使得可以量化预测因子的 importance，并表征春季物候对气候和人为驱动因素的非线性、特定于生态系统的响应。

3. 结果
3.1 SOS的空间分布
2001年至2020年，整个研究区域（MP）SOS的空间分布如图2所示。观察到明显的空间异质性。总体而言，SOS在研究区域的西南部出现较晚，在北部和东部出现较早。图2显示了2001年、2005年、2010年、2015年和2020年SOS的空间分布和时间趋势（a-e）。年度平均SOS的时间变化（f）。这种模式与植被分布密切相关。东北部主要由森林和灌木丛组成，SOS倾向于更早出现。相比之下，西南部主要被草地覆盖，SOS相对延迟，这可能是由于水分限制更强以及植被对早期温度升高的反应较慢。此外，位于东南部的草地显示出比其他草地地区更早的SOS，这可能与放牧或土地利用强度等人为因素有关。在整个研究期间，平均SOS大约为125日期（DOY），在区域尺度上没有观察到显著的整体时间趋势。然而，空间上的趋势差异是明显的。特别是，从2001年到2020年，以草地为主的南部地区显示出SOS延迟的趋势。

3.2 SOS的时间趋势
图3展示了2001年至2020年Sen斜率（Sen’s slope）、MK Z值及其衍生的SOS趋势分类的空间模式。Sen斜率估计器（）量化了物候变化的幅度，而MK测试（MK test）评估了检测到的趋势的统计显著性。图3显示了2001年至2020年SOS斜率估计器β值（a）、SOS MK测试Z值（b）和SOS趋势（c）的空间分布。Sen斜率图（图3a）揭示了MP范围内SOS变化的显著空间异质性。负的β值表明SOS提前，主要出现在研究区域的北部和东北部，特别是在肯蒂（Khentii）和杭盖（Khangai）山脉的森林和山区。在这些区域，SOS的年提前率在当地超过-0.8天。正的β值表明SOS延迟，主要集中在南部和西南部的干旱草原，其中最显著的延迟超过0.8天/年，特别是在沙漠-草原过渡带。变化幅度向中央高原减小，那里β值接近零。

MK Z值的空间分布（图3b）进一步明确了这些趋势的统计稳健性。北部森林-草原生态带的大片区域显示出高于2.58的Z值，表明对观察到的SOS提前的信心很高。相比之下，中部和南部大片草原的Z值低于1.645，表明这些水限制区域的延迟或微小波动在统计上无法与年际変动区分开来。综合趋势分类（图3c）综合了SOS变化的方向和显著性。ES?和S?类别主要集中在北部山区森林和相邻的灌木丛中，这里变暖可能对春季物候产生了主导和一致的影响。相反，ES+和S+类别分布稀疏，主要局限于西南部的沙漠草原，表明局部干旱压力或其他环境限制系统地推迟了植被的绿化开始。在排除了这些非显著区域后，剩余区域（即有显著趋势的区域）约34.56%显示出春季物候的延迟。在这个延迟子集中，28.99%、37.85%和33.16%分别对应于极显著延迟、显著延迟和轻微显著延迟，这些类别主要集中在南部草原区域。相比之下，其余65.44%的显著变化区域显示出SOS的提前趋势。在这个提前子集中，31.03%、39.68%和29.29%分别对应于极显著提前、显著提前和轻微显著提前。这种模式突出了不同植被类型的物候敏感性差异，森林显示出更强的提前趋势，而草地则表现出较弱且大部分不显著的响应。

3.3 不同植被类型之间的差异
尽管研究区域的整体SOS趋势不显著，但在单独分析植被类型时出现了明显的差异。森林、灌木丛和草地的春季物候的时间和趋势有显著差异（图4）。图4显示了2001年至2020年不同植被覆盖区域的SOS趋势。森林的SOS最早，平均值为107.72 DOY，并显示出约每十年6.82天的显著提前趋势（n = 20, t = ?2.98, p < 0.05）。灌木丛也显示出类似的趋势，平均SOS为107.59 DOY，每十年显著提前约6.36天（n = 20, t = ?2.55, p < 0.05）。森林和灌木丛的趋势高度一致。相比之下，草地的SOS最晚，平均值为137.82 DOY，并且随时间没有显著趋势，仅有约每十年0.75天的轻微延迟（n = 20, t = 0.16, p = 0.787）。为了进一步研究年际变化，分析了每种植被类型的SOS时间动态（图5）。从2001年到2020年，草地的SOS始终比森林和灌木丛晚。森林和灌木丛显示出类似的时间模式，有明显的提前趋势。两种植被类型的平均SOS从2001年的大约106 DOY下降到2020年的约94–95 DOY。然而，草地显示出不同的模式。从2001年到2015年，SOS相对稳定，介于133到136 DOY之间。然而，2020年出现明显的延迟，平均SOS增加到约143 DOY。这种差异表明，与森林和灌木丛生态系统相比，草地物候可能对短期环境波动或外部干扰更加敏感。

3.4 SOS的驱动因素
XGBoost模型对不同植被类型的性能通过R2值来展示。该模型对所有植被类型的SOS都具有良好的预测能力。对于森林、灌木丛和草地，R2值的范围从最低的0.62到最高的0.71不等。这些结果表明，模型能够有效捕捉SOS与环境驱动因素之间的关系，在草地生态系统中具有相对更好的泛化性能。不同植被类型中驱动因素的相对重要性各不相同（图6）。总体而言，温度（TEM）成为所有植被类型中的主要驱动因素，表明热条件在控制MP内的SOS方面起着根本性作用。此外，与人类相关的因素通常比气候相关因素贡献更多（图6）。对于森林（图6a），人口密度（POP）的贡献最低（2.5%），而温度（TEM）、降水（WIN）和密度（DENSITY）的作用更为重要，分别贡献了24.0%、23.2%和19.8%。这表明森林物候主要受气候条件调节，同时也受到放牧相关干扰的影响。对于灌木丛（图6b），温度（TEM，28.8%）、海拔（DEM，21.9%）和降水（PRE，15.9%）的贡献相对显著，表明温度条件和地形因素在塑造SOS方面起着关键作用。与森林相比，灌木丛的SOS似乎受到气候和地表条件共同的影响。对于草地（图6c），温度（TEM，26.5%）被认为是最重要的因素，其次是降水（PRE，24.1%）和海拔（DEM，20.6%），而人口密度（POP）的贡献最低（3.0%）。温度（TEM）和降水（PRE）的更强影响表明草地物候更直接受水热条件控制，反映了草地生态系统对水和热量可用性的高敏感性。

总体而言，这些结果突出了不同植被类型中与SOS相关因素的共性和差异。虽然温度在所有生态系统中都起着主导作用，但其他因素的相对重要性有所不同，表明在分析干旱地区物候动态时应考虑植被特定的响应。

3.5 非线性效应和阈值
ShAP dependence图（图7）揭示了不同植被类型中SOS与其相关因素之间的显著非线性关系和类似阈值的响应。尽管温度（TEM）在三种植被类型中对SOS的影响始终最显著，但SOS响应曲线有所不同。在森林和灌木丛中观察到线性趋势，而在草地中则是非线性的。对于森林（图7a–d），温度（TEM）与SOS显著负相关。当TEM低于约-1°C时，Shap值迅速增加，表明在低温条件下TEM对SOS有显著的延迟效应。当TEM超过-1°C时，这种延迟效应逐渐转变为提前效应。WIN通常表现出积极的影响，在风速达到约0.8米/秒后，SHAP值显著增加，表明较高的风速会延迟SOS的到达。DENSITY与SOS呈正相关，但其增长率在每公顷30个Head之后减缓，过高的DENSITY会延迟SOS。DEM和SOS显示出非线性的负相关，波动范围在600-100米之间。超过这个范围后，海拔的升高会略微提前SOS。图7显示了森林中SOS驱动因素的SHAP依赖性（a-d）、灌木丛中SOS驱动因素的SHAP依赖性（e-h）以及草地中SOS驱动因素的SHAP依赖性（i-l）。对于灌木丛（图7e-h），TEM也表现出强烈的负效应，与森林类似，但响应范围更广。观察到一个明显的阈值在-2°C左右；低于这个阈值时，SOS会显著延迟。DEM显示出U形关系，在600-2000米的范围内SHAP值较低，表明中高度地区可能不利于SOS的提前。PRE通常表现出负相关，阈值在550毫米左右。超过这个阈值后，PRE的增加通常会略微提前SOS。相比之下，DENSITY表现出正的非线性效应；适度的放牧会提前SOS，而过度放牧则会延迟SOS。对于草地（图7i-l），非线性响应更为复杂。TEM显示出明显的U形模式，其SHAP值在-3°C附近达到最小，表明极低和极高的温度都会延迟SOS，而中等温度则可以提前SOS。PRE总体上表现出负面影响，阈值接近350毫米，表明PRE不足可能会抑制早期生长。DEM表现出正相关，表明海拔的升高会持续延迟SOS。WIN显示出非线性的正相关，具有明显的阈值在0.3°C和1.3°C左右；在这些范围之外，风速会显著延迟SOS。这些结果突显了SOS对环境驱动因素的非线性和植被依赖性响应。温度仍然是一个主要因素，但其影响的形式和强度因植被类型而异。建立的阈值进一步表明，一旦环境条件超过某些临界值，植被响应可能会发生突变，强调了考虑非线性因素的重要性。

4. 讨论
4.1. 植被类型对环境驱动因素的特定响应
尽管相同的环境变量影响所有植被类型的SOS，但它们的响应模式存在显著差异，反映了它们对气候变化的敏感性差异[1]。森林和灌木丛在春季物候上表现出一致且显著的提前趋势，而草地则表现出弱且不显著的方向性变化。这种差异表明木质植被对年际气候变化的响应更为一致，而草地物候则更受限制性环境条件的强烈调控[33]。森林和灌木丛主要分布在相对湿润且地形有利的北部地区，那里的湿度减少了春季物候发展的限制。相比之下，草地分布在更干旱的南部地区，水资源的限制增加了降水变化对温度驱动的物候响应的缓冲作用[16]。可以从两个地理角度理解这种差异。首先，植被类型在MP上的分布遵循明显的气候梯度，森林和灌木丛集中在相对湿润的北部山区，而草地则主导干燥的南部平原。这种空间格局意味着草地天生受到更强的水资源限制，使其物候对PRE的变化更加敏感[34]。由于PRE在调节干旱/半干旱地区物候的温度依赖性中起着关键作用，森林和灌木丛的缓冲响应可能部分反映了它们占据更有利的水热生态位[35]。其次，最近研究中报道的不同气候敏感性——灌木丛显示出最高但呈下降趋势的敏感性，草地表现出中等且呈上升趋势的敏感性，而森林保持最低且最稳定的敏感性——表明不同的植被类型具有不同的缓冲或放大外部气候信号的能力[36]。人类活动也进一步加剧了这些差异。放牧压力在高原上广泛存在，并通过改变植被结构和土壤暴露来调节物候时间[37]。然而，在区域尺度上，其影响可能次于气候控制，特别是在气候更加受限的草地生态系统中[38]。适度的放牧可以通过清除残留物和增加土壤对太阳辐射的暴露来促进SOS的提前，而过度放牧则会通过土壤压实和降低植被活力来延迟SOS[39,40]。不同植被类型的水资源可用性的不同基线条件可能进一步调节放牧效应的表现——已经处于较大水分压力的草地可能会对综合的气候和人为压力表现出更复杂的响应[41]。总之，这些发现强调了在评估干旱地区生态系统的物候动态和制定生态管理策略时，纳入植被特定响应和区域差异的重要性。

4.2. 干旱地区物候的独特特征
在MP上观察到的物候响应与其他干旱和高纬度生态系统的发现大体一致，但也突出了干旱地区物候的一些独特特征。森林和灌木丛中的SOS提前与北半球温带生态系统中报告的春季物候加速的一般模式相符，表明木质植被对近期变暖具有广泛的敏感性[16,42,43]。然而，草地中没有显著的趋势，这与更湿润的草地系统中观察到的模式明显不同，在那里变暖通常会导致一致的物候提前[44,45,46]。这种差异可能反映了干旱和半干旱生态系统中水资源限制的主导作用[47]。在这样的环境中，温度驱动的物候加速可能被不足的水资源可用性所抵消甚至抑制[48,49]。因此，草地物候受到热力和水文限制之间的更复杂相互作用的影响，使其对温度升高的响应较低[50,51,52]。降水和风相关的因素的重要性进一步强调了干旱生态系统多压力因素的本质。与能量受限的系统相比，MP中的物候受到水资源可用性、大气需求和表面条件的共同调节，这些因素共同塑造了植被对气候变化的特定响应[53,54,55]。非线性响应和阈值行为进一步支持了这一解释。在水资源有限的生态系统中，植被响应很少遵循线性的气候-物候关系，而是当环境条件超过临界阈值时表现出突然的变化[56]。这突出了在解释干旱地区物候动态时考虑非线性机制的重要性。最后，模型中人类活动的有限独立贡献，特别是在草地中，并不意味着放牧和人口压力在生态上是无关紧要的。相反，正如MP上的先前研究所指出的，土地利用强度的信号经常与区域尺度上的气候变化混淆在一起[57,58,59]。

4.3. 限制和未来工作
尽管这项研究获得了许多见解，但仍需承认一些限制，这也指出了未来研究的有希望的方向。首先，从NDVI时间序列中提取SOS存在不确定性。MOD13Q1的16天合成间隔虽然适合区域尺度监测，但在植被稀疏的地区可能会掩盖微妙的物候转变，尤其是在土壤背景效应明显的地方[60,61]。此外，这里使用的固定20%动态阈值虽然在物候研究中广泛使用，但并未考虑NDVI季节性曲线形状的年际变化，这在极端气候条件下可能会发生变化[62]。其次，某些驱动因素数据集的空间分辨率限制了分析。ERA5-Land气候变量的原生分辨率约为11公里，可能无法充分解析细尺度的地形和微气候对物候的影响，尤其是在山区。第三，虽然本研究中使用的建模框架在检测非线性关系方面非常强大，但在因果推断方面存在固有局限性。XGBoost量化了预测因子的的重要性，但并未确定因果关系。尽管SHAP分析通过隔离边际贡献减轻了这一问题，但残余的混杂因素无法完全消除。草地中识别出的U形温度响应虽然在生态上合理，但仍需要通过控制实验或基于过程的模拟来进一步验证，以明确代表这种阈值行为背后的生理机制。第四，分析仅限于2001-2020年期间土地覆盖保持稳定的地区。这种设计有效地将物候信号与土地利用变化的混杂效应分离，但也排除了正在经历积极转变的地区，如农作物扩张、城市化或植被退化。因此，这里记录的物候响应可能是一个相对保守的信号，MP上人类引起的物候变化的全部幅度可能被低估了。基于这些限制，未来的研究方向有几个：更高时间分辨率的遥感数据，结合自适应阈值技术，可以提高干旱景观中SOS检测的准确性。整合基于现场的放牧强度调查和高分辨率动物追踪数据有助于区分真正的生态作用和与粗糙空间数据集相关的统计伪影。此外，将机器学习方法与基于过程的物候模型相结合，为超越预测重要性、迈向机制理解提供了有希望的途径，特别是对于本研究中揭示的非线性阈值行为。最后，扩展分析范围以包括过渡性和人为改造的景观，将提供更全面的评估，了解气候强迫和地表变化如何共同塑造MP上的植被物候。

5. 结论
基于对2001年至2020年间MP稳定植被区域的MODIS衍生SOS的分析，这项研究揭示了春季物候趋势的明显空间异质性。北部和东北部山区的森林和灌木丛生态系统显示出显著的提前趋势，平均SOS分别在大约107天和101天发生，并且每十年分别提前6.8天和6.4天。相比之下，主导南部和西南部干旱地区的草地生态系统显示出较晚的平均SOS（138天），并且没有显著的时间趋势，在沙漠-草原过渡带局部地区延迟超过每天0.8天。排除非显著区域后，剩余区域中有65.4%显示出SOS的提前趋势，而34.6%显示出延迟趋势。温度始终是所有植被类型的主要驱动因素，但其相对重要性和响应形式存在显著差异。SHAP分析进一步识别出强烈的非线性，包括草地的U形温度响应和350毫米附近的降水阈值，以及森林中的风速和牲畜密度阈值。这些发现强调了干旱地区物候调节的植被特定性和非线性本质，为预测生态系统对气候变化的响应提供了更机械的基础。它们还强调了将区域和植被特定的阈值纳入生态管理和气候适应策略的重要性，特别是在水资源受限的草地和风吹暴露的森林中，不同的物候敏感性可能对碳循环、饲料可用性和生态系统稳定性产生连锁影响。