黄河流域是连接青藏高原、黄土高原、华北平原的生态廊道，在中国经济社会发展和生态安全方面具有十分重要的地位。植被作为陆地生态系统的主体，是碳循环、水循环与能量循环的重要纽带，也是气候变化最敏感的指示标志^［1-2］，被认为是实现流域经济社会可持续发展的基础。随着退耕还林（草）、“三北”防护林和天然林保护工程等一系列生态恢复工程的实施与推进，黄河流域植被覆盖呈显著增长趋势^［3-4］。与此同时，全球变暖进程的持续加快使得干旱和热浪等极端天气事件并发频率与强度不断增加，给粮食生产和生态环境带来巨大而深远的影响，极大加剧了植被生态系统的生态风险^［5-6］。黄河流域位于东亚季风区，近年来面临着水土流失、土地荒漠化、人地矛盾日益突出等生态环境问题，且其植被生长对于全球气候变化有着较为敏感的响应过程^［7-8］。因此，阐明黄河流域植被覆盖变化及其影响因素是实现生态保护与高质量发展的关键环节。

气候变化和人类活动的双重扰动加剧了全球水分循环过程，由此导致的全球植被覆盖变化及其互馈机理成为学界高度关注的科学问题^［9-12］。当前国内外学者已从不同时空尺度上开展了植被覆盖变化的研究工作，其中Kaufmann等^［13］指出气温与CO₂浓度的升高驱动了北半球中高纬度陆地地区的植被生长不断增强；Chen等^［14］则从全球尺度探讨了植被覆盖变化特征，并指出中国和印度两国植被覆盖增加趋势最为显著，其中中国增长趋势最显著的地区位于中国北方黄河流域-黄土高原一带。区域尺度上，赵安周等^［15］和周德成等^［16］分析了2000—2014年黄土高原植被指数（NDVI）呈现显著增加趋势，提出了人类活动对于植被同时具有积极、消极的双重影响作用，并论证了它在生态环境效应与社会经济效益方面均有着较大的正面作用。应当指出的是，已有研究多探讨了生态工程对黄河流域植被变化的贡献，而就多重气候因子对植被变化的共同作用及贡献报道较少，且近20年黄河流域植被时空变化的最新趋势尚不明晰。因此，阐明黄河流域植被覆盖变化最新趋势及其对气候因子的综合响应机制有助于开展植被覆盖变化的预测研究。

本文基于MODIS NDVI数据、土地利用数据和气象数据集，以集合经验模态分解、趋势分析和随机森林回归等方法，全面分析黄河流域2001—2020年植被覆盖时空变化特征，并通过构建统计模型方式，定量分析植被覆盖变化对不同气候影响因子的响应特征，阐明植被覆盖变化的驱动因素与作用机理，旨在提高学界对黄河流域植被覆盖变化特征的最新认识，为科学开展流域生态可持续治理提供参考依据。

黄河西起巴颜喀拉山，在山东省注入渤海，流域面积约74万km²，由西至东依次流径青藏高原、内蒙古高原、黄土高原、黄淮海平原四大地理区域。该流域主要位于半干旱区-半湿润区气候过渡带，具有太阳辐射较强、降水集中且分布不均匀、蒸发量大等特点。上游地区主要面临荒漠化、凌汛等问题；中游水土流失较为严重；下游地区人口稠密，城市建设用地、农垦区密集分布，环境承载压力加大，面临着地上悬河化、用水紧张等问题。为改善生态环境，国家层面已实施了一系列生态恢复工程措施，显著改变了流域土地利用格局，提高了植被覆盖度（图1）。

遥感数据为2001—2020年的MOD13A3卫星产品数据中的NDVI数据，来源于美国航空航天局（https：//ladsweb.modaps.eosdis.nasa.gov/），空间分辨率为1 km，时间分辨率为30 d，本研究利用MCTK（MODIS Conversion Toolkit）完成NDVI数据提取、几何校正、重投影等工作。MODIS数据具有光谱信息丰富、重返周期短、宽视域等明显优势，被广泛应用于全球及区域植被变化、生态环境监测、自然灾害与极端气候事件响应等方面研究^［17］。土地利用数据来源于GlobeLand30全球地表覆盖数据（http：//www.globallandcover.com/）。

气象数据主要包括气温、降水、辐射、大气CO₂浓度等，来源于国家科技基础条件平台、国家地球系统科学数据共享服务平台——黄土高原科学数据中心（http：//loess.geodata.cn）^［18］。本研究中主要选取2001—2020年5—9月作为黄河流域植被生长季，并分别对降水、气温和辐射等逐月数据集进行年数据合成。其中大气饱和水汽压差（VPD）、辐射量（RAD）和CO₂浓度等年合成数据则分别来源于ERA5再分析数据集（https：//cds.climate. copernicus.eu/）和GOSAT卫星近地面逐月温室气体观测数据（https：//data2.gosat.nies.go.jp/ index_en.html）^［19-20］。需要说明的是，本文在NDVI趋势分析和NDVI与降水、气温的偏相关分析中统一采用1 km分辨率进行分析，而在植被变化驱动因素建模时将所有数据统一重采样为0.25°空间分辨率。

本研究采用集合经验模态分解（Ensemble Empirical Mode Decomposition，EEMD）和非参数化趋势度模型进行趋势分析，其中EEMD是针对经验模式分解（EMD）模型存在的不足而提出的一种噪声辅助数据分析方法，其原理是在EMD基础上加入一组或多组白噪声信号，从而对原始时间序列数据进行平稳化处理。EEMD模型具有完全的自适应性、非线性等固有优势，并且其不受数据的线性限制，可更加准确地描述数据的变化趋势^［21-22］。

式中：imf_i（t）是EMD分解得到的第i个IMF；r_n（t）是分解筛除n个IMF后的信号残余分量，常常代表信号的直流分量或信号的趋势。

SEN趋势分析方法的优点是不需要样本服从一定的分布，并且不受异常值的干扰，对测量误差或离群数据具有较强的规避能力^［23-24］。

式中：β为植被变化趋势；x_i，x_j分别为i，j时刻植被指数的值。当β>0时，表示植被覆盖呈增加趋势，而当β<0时，表示植被覆盖呈减小趋势。

偏相关分析指在多个要素所构成的模型或系统中，在单独研究两个要素之间的相互关系的密切程度时，将其他要素的影响作用视为保持不变的恒定常量，即只分析两个变量之间相关程度且不考虑其他要素影响过程，所得结果为偏相关系数^［25］。

式中：r_xy·z为变量z固定后变量x和y的偏相关系数，即xy相关中剔除z的影响；r_xy为变量x与变量y的相关系数；r_xz为变量x与变量z的相关系数；r_yz为变量y与变量z的相关系数。在开展偏相关分析前，本研究对时间序列数据进行了去趋势处理，以消除不同时间序列数据变化趋势对于数据内部波动变化的影响。

随机森林回归原理是从n个原始训练数据集中，有放回地随机抽取k个样本，并生成k个决策树建模结果，依据所有决策树的建模结果，以其平均值来获取最终的预测结果^［26］。本研究在探究NDVI对于气象因子敏感性时，选取了降水、气温、饱和水汽压差、辐射量、CO₂浓度等因子构建统计模型。由于不同因子数值范围相差较大，难以同时进行比较，因此在建模前对不同因子进行了0—1的标准化处理。

近20年黄河流域植被覆盖呈显著增加趋势（图2），平均增速为 0.055/10a（P<0.05），且2010之前增速（0.067/10a）大于2010年之后的增速（0.051/10a）。通过对2001年与2020年土地利用与覆被类型转移矩阵统计可以看出，退耕还林还草、农用地开垦、农用地保持不变区域的植被覆盖增速相对较高，分别达到了0.082/10a、0.076/10a、0.067/10a，均通过了显著性水平0.05检验；其他区域的植被覆盖增速相对较低，说明在人类活动干预程度相对较高的退耕还林还草与农用地区域，植被生长同时受到了自然因素与人为干预因素的双重影响，因而植被增长趋势达到了一个相对较高的水平；而自然植被下垫面、非植被类型向植被类型过渡等区域内，植被增长趋势相对平缓。

空间趋势上，黄河流域有90.36%（Sen趋势结果）或91.22%（EEMD趋势结果）的区域植被覆盖呈增加趋势，且呈极显著增长的区域（P<0.01）占黄河流域总面积的51.03%（Sen趋势结果）或81.18%（EEMD趋势结果，图3，表1）。此外，植被覆盖增加趋势超过0.1/10a的区域集中在延安市、榆林东部、河套平原、宁夏回族自治区南部以及甘肃省东南部等地区，占整个研究区域的17.54%（Sen趋势结果）或21.76%（EEMD趋势结果）；而植被覆盖减小的区域占研究区的9.64%（Sen趋势结果）或8.78%（EEMD趋势结果），且集中在黄河源、银川-兰州-包头-西安-郑州城市周边、黄河下游黄淮海平原等地区，说明黄河上游青藏高原高海拔地区的植被退化、黄河中游各主要城市建设用地扩张、黄河下游-黄淮海平原地区农用地开垦等是造成黄河流域植被覆盖下降的主要原因。

研究时段内，黄河流域降水、气温、VPD和CO₂浓度均呈上升趋势，而辐射量整体呈现下降趋势，且上述气候因子均存在较为明显的年际波动特征（图4）。降水、气温、CO₂与NDVI呈现正相关关系的区域分别占研究区总面积的85.11%、63.14%和72.72%，说明这些因子整体上会对植被生长起到促进作用；而VPD、辐射量与NDVI主要呈负相关关系，二者分别占研究区总面积的59.69%和69.52%，说明VPD和辐射量的增加可能会对植被生长起到抑制作用。同时，VPD、辐射量与NDVI呈现强负相关的区域集中在甘肃东南部、宁夏南部、陕西西部等地区，说明VPD和辐射量对上述地区植被生长具有更为显著的抑制作用，而近年来辐射量的下降可能缓解对植被生长的抑制作用。

为进一步探究NDVI对不同气候影响因子的敏感性，选取上述与植被生长密切相关的气候因子，结合随机森林多元回归方法，构建了黄河流域NDVI模拟模型，模拟了黄河流域NDVI时序数据（受GOSAT-CO₂数据可获取时间范围限制，选取2010—2019时间范围内的各时序数据，构建此模型）。将观测时序数据与模拟时序数据逐像元对比分析可以看出，两者在空间格局上具有较好的一致性，且相关系数在0.8以上的区域占研究区的60.47%，而相关系数在0.5以下的像元呈零星分布状态（图5A、B）。模拟数据与观测数据总体回归系数解释率为0.925，可以满足气候因子对植被覆盖变化影响的定量分析。

基于上述构建的NDVI模拟模型，本研究分别控制其他变量保持不变，逐个将不同气候因子的数据放大1.5倍与2倍进行模拟实验，进而将模拟结果与原始模拟的数据做差运算，得到不同实验条件下的结果（图5C）。可以看出，降水量、气温、CO₂浓度在两种升高情景下均表现出正向促进作用，即NDVI均较原始模拟数据高，其中降水量和CO₂浓度表现为持续升高，而气温则在1.5倍情景下大于2倍情景下的模拟结果，说明植被生长有着适宜的气温范围，气温数值超过此范围后，可能会对植被生长起到抑制作用。VPD和辐射量的增加对黄河流域NDVI表现出一致且持续的下降趋势，两者下降幅度均超过降水量、气温和CO₂浓度的正向促进幅度。其原因可能是VPD的提高会加剧大气干旱程度、降低植物气孔导度，限制植被的光合作用并使得植被长势、生产力持续下降；而辐射量的增加可能会进一步增加蒸散，进而造成植被水分散失增多，降低气孔导度、进一步导致植被长势的下降。上述分析结果表明，在生态恢复工程等人类活动干预下，黄河流域植被覆盖呈显著增加趋势，然而气候变化对植被长势的潜在风险却在增加。因此，加强多气候因子对植被生长影响的协同作用分析有助于全面理解植被变化趋势。

2001—2020年黄河流域植被覆盖呈显著增加趋势，平均增速为0.055/10a；且黄河流域植被恢复呈现出阶段性，2001—2010年的平均增速（0.067/10a）高于2011—2020年的平均增速（0.051/10a），整体上均通过置信度为0.05的显著性检验。

空间格局上，Sen趋势和EEMD方法检测出的黄河流域植被覆盖呈极显著增长区域（P<0.01）分别占整个研究区总面积的51.03%与81.18%，且平均增速大于0.1/10a的地区分别占研究区总面积的17.54%与21.76%。

回归分析结果表明，降水、气温和CO₂浓度的提高会促进植被生长，且黄河流域内气温对促进植被生长有着最高的贡献率；而VPD与辐射量的提高则对于植被生长有着持续的抑制作用。

黄河流域生态工程的实施使得原先通过砍伐森林扩大耕地面积的现象已得到根本逆转，自然植被呈现显著恢复的趋势，但植被恢复趋势是否具有可持续性、生态工程的实施能否与地区经济发展相协同等方面，仍存在问题^［27］。人类活动对于黄河流域植被生长积极/消极作用的良好区分、预测未来干旱等极端气候事件对于植被所产生的胁迫作用等问题^［28］，同样有待于进一步的研究改进与提升。本研究从宏观格局上揭示了植被对于气象因子变化的敏感程度，对于深入认识植被对气候变化的响应过程及其规律具有一定的意义。本研究存在的不确定性主要表现在：①NDVI模拟时序数据的像元值区间，整体略低于真实NDVI时序数据，由于本研究中参与模拟NDVI的数据均为自然气象要素因子，说明退耕还林还草政策、生态工程等人工植被恢复措施的实施，使得NDVI以超出自然条件下的增长速率进行增长，最终促使真实NDVI的平均数值，高于模型模拟的平均数值。②本研究在完成探讨不同气候因子对于植被生长的作用及定量化贡献率后，尚未将自然、人为驱动力分离开来，并与全球变化背景下频度和强度不断增加的极端天气气候事件（干旱、洪涝、高温热浪等）相结合，充分探讨上述事件对于研究区内的植被长势与时空变化分布格局所造成的影响。