作品介绍

1 研究目的

基于MODIS NDVI植被指数、土地利用数据和气象数据集，辅以趋势分析、偏相关分析、马尔科夫转移矩阵变化分析、多元回归分析等方法，全面分析黄河流域2001-2020年植被时空变化特征，并通过构建统计模型方式，定量分析植被覆盖变化对不同影响因子的响应特征，阐明植被覆盖变化的驱动因素与作用机理，旨在提高公众和学界对黄河流域植被覆盖变化特征的最新认识，为科学开展流域生态可持续治理提供参考依据。

2 研究区数据与方法

2.1 研究区概况

黄河流域西起巴颜喀拉山，东至山东省注入渤海，流域面积约74万km2，由西至东依次横跨青藏高原、内蒙古高原、黄土高原、黄淮海平原四大地理区域。该流域主要位于半干旱区-半湿润区气候过渡带，气候上具有太阳辐射较强、降水集中且分布不均匀、以及蒸发量大等特点。上游地区主要面临荒漠化、凌汛等问题；中游水土流失较为严重；下游地区人口稠密，城市建设用地、农垦区密集分布，环境承载压力加大，面临着地上悬河化、用水紧张等问题。为改善生态环境，国家层面已实施了一系列生态恢复工程措施，显著改变了流域土地利用格局，提高了植被覆盖度。本项目重点研究各植被下垫面类型范围内植被变化的时空分布特征，并对于2001-2020年黄河流域非植被区域进行掩膜处理。

图1 2020年黄河流域土地利用图

2.2 研究方法

2.2.1 趋势分析

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

      (1)

式中，β为植被变化趋势，

分别为i,j时刻植被指数的值，当β>0时，表示植被覆盖呈增加趋势，而当β<0时，表示植被覆盖呈减小趋势。

本作品中，利用ENVI Task中的Slope K Calculator扩展工具，估算黄河流域NDVI时序数据的逐像元线性趋势。与此同时，编写SEN趋势分析的自定义IDL代码，对于NDVI线性趋势结果进行逐像元的显著性检验，详细过程请参见3.3章节部分。

2.2.2 偏相关分析

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

   (2)

式中，

为表示变量z固定后变量x和y的偏相关系数，即xy相关中剔除z的影响；为变量x与变量y的相关系数；为变量x与变量z的相关系数；为变量y与变量z的相关系数。

在开展偏相关分析前，本研究对时间序列数据进行了去趋势处理，以消除不同时间序列数据变化趋势对于数据内部波动变化的影响。具体实现方式为：首先逐像元地计算并得到线性趋势分析的斜率数值，再对原始数据减去一个以斜率为等差的数列即可得到了消除趋势后的新数据序列。本作品中，首先编写去趋势分析的IDL代码，对于NDVI/降水/气温时序数据进行去趋势波动分析。完成上述操作之后，利用ENVI Task中的Partial Correlation Coefficient扩展工具，对各个数据进行偏相关分析，详细过程请参见3.4章节部分。

2.3 技术流程设计

本作品技术流程大致划分为数据预处理、指标提取与判别分析、评价模型构建与技术应用三部分，具体介绍如下：

数据预处理：在ENVI-IDL、ArcGIS、HDF Explorer、NC Panoply等软件的支持下，对于MODIS数据、NC/HDF5气象数据、土地利用产品数据进行定标参数读取、重投影、提取目标专题信息等遥感图像预处理工作，并进行植被指数年最大值合成、降水/气温数据生长季合成、气象数据标准化处理等操作。

指标提取与判别分析：在各输入数据及趋势分析、偏相关分析、马尔科夫转移矩阵分析的支持下，分别从土地利用类型变化、植被动态驱动力因素、气象要素多元回归分析等方面对于黄河流域植被时空变化特征及其驱动因素，进行全面的评估与分析。

评价模型构建与技术应用：获得上述各评估模型的结果后，从不同角度对于黄河流域植被动态与生态环境效应进行归因分析、综合评估、专题制图、应用分析。

图2 项目技术流程设计

3 技术流程步骤

3.1 数据预处理

3.1.1 MOD13A3产品数据预处理

MOD13A3卫星产品数据中的NDVI数据，来源为美国航空航天局(https://ladsweb.modaps.eosdis.nasa.gov/)，空间分辨率为1km*1km，时间分辨率为30d，本研究利用MCTK(MODIS Conversion Toolkit)完成NDVI数据提取、几何校正、重投影等批处理工作。MODIS数据具有光谱信息丰富、重返周期短、宽视域等明显优势，在于全球及区域植被变化研究、生态环境监测、自然灾害与极端气候事件响应等方面，有着广泛的应用前景。

四景MODIS影像可完全覆盖黄河流域研究区（图3），本项目选取2001-2020时间范围内（20年，累计240个月）共计960景MOD13A3产品数据，首先选用HDF Explorer软件读取相关参数，在ENVI+IDL开发环境下调用MCTK,一键化批量完成960景数据的相应的预处理工作（图4）。

基于MCTK函数接口完成预处理工作后，借助ENVI Modeler建模工具，一次输入一年共计48景数据，构建模型一键化完成图像镶嵌、特征波段提取、时间顺序图层叠加、逐月植被指数最大值合成为年数据、空间裁剪等图像处理过程。在实现过程中，同时提取MOD13A3产品数据中的NDVI、EVI、NIRV（近红外与红外波段相乘）植被指数，进行对比分析。最终，结合指数提取效果与研究区实际状况（植被指数数值普遍低于0.8），选择NDVI参与后续的建模过程。

图3 NASA_MODIS检索研究区范围

图4 HDF参数读取及MOD13A3批处理展示

图5 ENVI Modeler合成植被指数年数据

3.1.2 MCD12Q1 产品数据预处理

MODIS土地覆盖类型产品数据包含多个分类方案，这些分类方案描述了一年来输入的Terra和Aqua传感器的观测数据得出的土地覆盖属性（该影像数据重返周期为一年）。初步的土地覆盖计划确定了国际地圈生物圈计划（IGBP）定义的17种土地覆盖类别，其中包括11种自然植被类别，3种衍生合成土地类别以及3种非植被土地类别。MODIS Terra + Aqua土地覆盖类型年度L3全球500 m 正弦投影网格产品（MCD12Q1）结合了五种不同的土地覆盖分类方案，这些方案是通过监督决策树分类方法得出的：

表2 MCD12Q1产品数据所含土地覆盖分类方案

编号	土地覆盖类型分类方案
1	IGBP全球植被分类方案
2	马里兰大学（UMD）计划
3	MODIS衍生的LAI / FPAR方案
4	MODIS衍生的植被第一净初级生产力（NPP）方案
5	植物功能类型（PFT）方案

同MOD13A3批处理方式，借助HDF Explorer、ENVI+IDL调用MCTK、ENVI Modeler建模的方式，一键化完成MCD12Q1数据预处理、IGBP全球植被分类方案提取、图像镶嵌、时序数据构建、空间裁剪等图像处理过程。本项目中，引入MCD12Q1数据的主要目的为提取植被区域（第一景、最后一景均为非植被类型的像元判定为非植被区）、同GlobeLand30土地利用利用数据提取结果进行对比分析等。

图6 HDF参数读取及MCD12Q1批处理展示

图7 ENVI Modeler合成IGBP年数据

表3为IGBP分类方案，根据项目实际需求，利用Band Math工具编写判别函数，将地物类别归类为林地、灌木、草地、湿地、农用地、城市建设用地、裸地、水域/冰雪八大类，各自的DN值分别赋值为1-8，判别函数如下：

(b1 gt 0.9 and b1 lt 5.1)*1+(b1 gt 5.9 and b1 lt 7.1)*2+(b1 gt 7.9 and b1 lt 10.1)*3+(b1 eq 11)*4+(b1 eq 12 or b1 eq 14 )*5+(b1 eq 13)*6+(b1 eq 16)*7+(b1 eq 17 or b1 eq 15)*8

完成归类操作后，以最前、最后两景数据提取植被区掩膜文件，判别函数如下（此公式代入趋势分析等ENVI Modeler建模工程文件中以流程化的形式实现）：

(B1 GE 1 AND B1 LE 5) OR (B2 GE 1 AND B2 LE 5) *1+(B1 GT 5 AND B2GT 5) *0

表3 IGBP全球植被分类方案

DN值	分类方案	DN值	分类方案	DN值	分类方案
0	背景值	6	稠密灌丛	12	农用地
1	常绿针叶林	7	稀疏灌丛	13	城市和建筑区
2	常绿阔叶林	8	木本热带稀疏草原	14	农用地/自然植被
3	落叶针叶林	9	热带稀疏草原	15	冰雪
4	落叶阔叶林	10	草地	16	稀疏植被
5	混交林	11	永久湿地	17	水

3.1.3 GlobeLand30产品数据预处理

30米全球地表覆盖数据GlobeLand30是中国研制的30米空间分辨率全球地表覆盖数据，采用WGS-84坐标系，以GEO-TIFF的格式进行存储与发布，每十年提供一期产品数据，分类方案如下：

表4 GlobeLand30分类方案

类型	代码	类型	代码
耕地	10	水体	60
林地	20	苔原	70
草地	30	人造地表	80
灌木地	40	裸地	90
湿地	50	冰川和永久积雪	100

利用ENVI Modeler工具，可同时对于2000、2020两期数据,一键化完成图像镶嵌、空间裁剪、基于Color Slice Classification工具将其转化为ENVI标准格式分类图像，以用于后续的专题制图与马尔科夫转移矩阵分析所需的输入数据。

图8 ENVI Modeler处理GlobeLand30数据

3.1.4 气温、降水气象数据预处理

降水、气温逐月数据来源于国家科技基础条件平台—国家地球系统科学数据共享服务平台-黄土高原科学数据中心(http://loess.geodata.cn)，空间分辨率为1km*1km，以NC格式进行发布。另外，根据官方文档说明，降水、气温的单位分别为0.1mm与0.1℃,需要在ENVI软件中进行除以10的波段运算，确保将数据格式转化为浮点型，以及计量单位的正确。

本项目选取5-9月作为黄河流域的生长季。借助ArcMap软件将NC数据转化为多波段ENVI标准格式时序数据之后，借助ENVI Modeler建模工具，构建模型一键化完成2001-2020全部数据读取、各年份5-9月波段光谱裁剪、逐像元特征统计（降水数据取累加值、气温数据取平均值）、时序数据构建、空间范围裁剪操作。

图9 ArcMap读取解译NC数据

图10 降水、气温栅格数据提取结果

图11  ENVI Modeler 合成生长季降水-气温数据

3.1.5 ERA5再分析数据集预处理

对于NC格式的ERA5再分析数据集，预处理方式同3.1.4中的对应部分。本项目所需的辐射量(RAD)数据可以直接从数据集中获取，代表大气干旱程度的大气饱和水气压差(VPD)需要通过数据集中的温度、露点温度，以及ENVI自带的DEM高程数据估算得到，计算公式如下：

  (3)

    (4)

   (5)

   (6)

    (7)

其中，

为ERA5数据集的气温数据；为ERA5数据中的露点温度数据；Z为ENVI自带的全球DEM高程数据。在IDL中编写波段运算扩展函数，将上述公式转化为IDL语言，调用ENVI Classic界面，在Band Math工具界面中将运算波段设置为输入时序数据的多波段文件，在将开尔文温度转化为摄氏度与消除图像异常值的基础上，一键化完成对于VPD的估算。

图12 IDL波段运算二次开发编程及调用ENVI Classic

得到VPD、RAD的估算结果之后，同样在ENVI Modeler中，利用迭代器(Iterator)+聚合器(Aggregator)的组合方式，一键化完成逐年生长季均值估算、时序数据构建、空间裁剪、数据标准化等图像处理过程。

图13 ENVI Modeler合成生长季VPD、RAD数据

3.1.6 GOSAT-CO2数据预处理

GOSAT以HDF5(h5)的格式储存数据，可利用ENVI5.6自带的HDF Dataset Browser，选取平均二氧化碳栅格图像，赋予数据集中的经纬度信息至此图像栅格，即可输出预处理结果。

根据图15，CO2提取结果存在空间分辨率低、数据缺失严重、有效数据仅覆盖2010年至2019年等问题。根据该数据在时间范围上的变化，远大于空间范围上的变化的特点，本项目以此数据逐月全球CO2浓度平均值，赋值到空间分辨率为2.5°*2.5°的ERA5数据的全部像元上，并在此基础上完成生长季均值计算、数据标准化等图像处理过程。

图14 ENVI HDF Dataset Browser读取GOSAT数据

图15 GOSAT全球平均二氧化碳浓度提取结果

3.2 土地利用与覆被类型变化探究分析

图16为MCD12Q1与GlobeLand30数据的土地利用与地物转移矩阵空间分布图,其中转移矩阵空间地图中的各类别由下方表5中各地物相互转换的情形整合得到。由结果可得，GlobeLand30数据可以有效提取得到众多的农用地开垦（第一时期的林-灌-草-湿地转化为第二时期的农用地）、退耕还林还草（第一时期的农用地转化为第二时期的林-灌-草-湿地）等转移类型的细碎斑块，而空间分辨率为500m的MCD12Q1数据对于评估景观破碎化程度较高的黄河流域内呈斑块状分布的退耕还林还草区域精度较差，故在接下来的研究分析过程中采用GlobeLand30数据。

另外，结合GlobeLand30前后两期分类图像与表5转移矩阵可得，黄河流域内草地下垫面占比由49.02%下降到45.65%；林地下垫面由40.00%下降到30.79%；耕地下垫面由11.3%上升到12.69%；人造地表由1.72%上升到3.65%。空间格局上，渭河流域、黄河流域以及黄河下游地区呈现出较为明显的城市建设用地与农用地扩张现象，初步说明黄河流域内主要自然植被占比有所下降、城市建设用地和农用地占比有所上升，即人类活动对于黄河流域植被时空分布格局，有着较大的影响作用。

表5 2001-2020年土地利用类型变化转移矩阵

2020地物类型占比(%)	2001地物类型占比(%)
	耕地	林地	草地	灌木	湿地	水体	人造地表	裸地	冰川/积雪
耕地	87.03	2.85	4.9	7.09	7.46	16.68	20.37	3.04	0
林地	1.47	88.49	4.04	8.79	0.26	0.65	0.3	0.17	0
草地	4.91	8.04	85.53	53.42	12.88	7.11	1.71	28.42	4.64
灌木	0.11	0.24	1.43	26.13	0.17	0.23	0.14	0.32	0
湿地	0.2	0.02	0.48	0.07	73.37	5.42	0.02	0.05	0
水体	0.59	0.08	0.24	0.49	4.91	67.89	0.22	0.29	0
人造地表	5.59	0.2	0.93	1.05	0.37	1.43	77.16	0.75	0
裸地	0.1	0.08	2.12	2.96	0.58	0.59	0.08	66.82	6.2
冰川/积雪	0	0	0.32	0	0	0	0	0.13	89.15

3.3 植被指数时空变化特征

在利用Slope K Calculator扩展工具进行NDVI趋势分析的基础之上，编写IDL自定义MK Trend程序，在关键函数中设置P=0.05、0.01的显著性阈值（程序截图中为双侧显著性检验，故实际输入的参数为0.025、0.005，未通过显著性检验的像元将被忽略为背景值），将输出结果同扩展工具输出的NDVI趋势分析结果整合，输出显著性检验等级结果，详细判别函数为：

(b1 gt 0 and b2 ne 0)*1+(b1 gt 0 and b2 eq 0 and b3 ne 0)*2+(b1 gt 0 and b2 eq 0 and b3 eq 0)*3+(b1 lt 0 and b2 eq 0 and b3 eq 0)*4+(b1 lt 0 and b2 eq 0 and b3 ne 0)*5+(b1 lt 0 and b2 ne 0)*6

其中，b1,b2,b3分别为NDVI趋势结果、P=0.05显著检验下的有效值部分、P=0.01下的有效值部分，1-6重分类结果值依次对应图19B中的各重分类类别。

图17 IDL显著性检验参数设置

图18 各土地利用转移类型下NDVI趋势分区统计

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

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

图19 NDVI趋势时空分布特征

表1 2001-2020年植被覆盖变化趋势显著性检验统计表

变化趋势	变化程度	分级标准	像元个数	百分比/%
增长	极显著增长**	Slope>0且p<0.01	460736	51.03
	显著增长*	Slope>0且0.01≤p<0.05	118097	13.08
	不显著增长	Slope>0且p≥0.05	245733	27.21
减少	不显著减少	Slope<0且p≥0.05	63395	7.02
	显著减少*	Slope<0且0.01≤p<0.05	7307	0.81
	极显著减少**	Slope<0且p<0.01	7671	0.85

注：显著性检验方式为t检验，*表示p通过置信度为0.05的检验，**表示p通过置信度为0.01的检验

3.4 植被指数与气象因子偏相关分析

在基于NDVI、气温、降水多波段时序数据进行偏相关分析之前，结合其趋势分析结果，编写IDL自定义去趋势程序，输出各因子去趋势后的多波段时序数据，导入Partial Correlation Coefficient扩展工具中进行偏相关分析，IDL程序核心部分与偏相关分析输出结果如下所示：

图20 IDL去趋势波段分析实现过程 

图21 NDVI时序数据去趋势效果对比

图22A-B分别为NDVI与降水、气温因子的偏相关分析空间专题地图。由结果可得，二者分别在约85.11%与63.14%的空间上呈现正相关关系，初步说明黄河流域大部分区域内，降水量与气温的升高均对植被长势有着正向促进作用。再者，二者在空间格局上正相关程度最高的区域均主要集中于400mm等降水线、半干旱区-半湿润区分界线以西、以北地区，说明水热条件为制约此区域植被生长的决定性主导因素。与此同时，根据二者的相关性分级结果，NDVI与降水的偏相关结果中的较强正相关性部分（相关系数大于0.3）所占比重较多，整体样式为单峰偏态分布；NDVI与气温的偏相关结果中弱正相关性部分(相关系数介于0与0.3之间)所占比重较多，整体样式为单峰正态分布。上述结果表明，NDVI分别与降水、气温的相关性程度存在明显的差异性。另外，由于植被长势对于水热条件变化的响应存在一定的滞后性，因此基于年尺度的NDVI与气象要素的偏相关分析，无法有效获取部分现象背后的机理与规律。因此，尚需进一步分析探究气象要素对于植被长势的影响作用。

由图22C-D可得，黄河流域5-9月生长季内，逐年累计降水量与平均气温整体均呈现出上升的“暖湿化”趋势；且二者在黄河流域植被区内约75.16%的区域内呈现负相关关系，同时，二者相关系数绝对值小于0.2的极弱相关性部分的占比达到了约59.99%，说明降水-气温的耦合性较弱，尚需引入更多的气象因子要素，对于植被长势的促进/抑制作用与贡献率的高低，进行定量化的评估与分析。

图22 NDVI与气象因子偏相关分析结果

3.5 植被动态驱动力分析

以上述部分NDVI的趋势分析、显著性检验、分别与降水/气温偏相关分析结果作为输入数据，依据2.3部分图2的植被动态驱动力分析判别模型的流程，编写Band Math公式进行掩膜提取、函数判别等过程，并用ENVI Modeler流程化、一键化输出植被动态驱动力分析的重分类结果，建模过程以及输出结果如下所示：

图23 植被动态驱动力分析建模过程

通过定量分析气候变化和人类活动对植被活动的影响，具体结论如下：

（1）受复杂因素驱动退化面积占比为0.41%，受气象因子驱动退化区域为1.25%，说明上述两类植被退化并非黄河流域植被动态的主导模态。在空间上，植被退化区域主要分布于关中平原城市群、黄淮海平原农垦区以及青藏高原地区，且气象因子和复杂因素影响的退化区域相互交织分布。

（2）受生态建设驱动绿化区域占比为22.64%，受多因素驱动绿化区域占比为2.27%，受气候因素驱动绿化区域占比为39.18%。说明黄河流域植被恢复以气候因素驱动为主，生态建设驱动为辅，且二者的地理分界线大致与半湿润区/半干旱区、400mm等降水线大致相同。空间上，NDVI受人类活动驱动绿化区域，连片分布于甘肃省东南部、陕西省中部、山西省西南部、毛乌素-库布齐沙漠北部以及黄河下游部分地区，基本与NDVI趋势呈现极显著增长的范围吻合；受气候因子驱动绿化区域，主要分布于陕西省中部、宁夏回族自治区南部、青海省中部等。

（3）无显著相关区域表征含义为：植被变化稳定。空间上主要分布于2类区域：一是受人类活动影响的耕地、城市区，如黄河下游沿岸地区、关中-中原城市群等；二是分布于植被相对稳定的森林、草地区，如青藏高原区、鄂尔多斯南部、宁夏回族自治区中部等。

图24 植被动态驱动力分析空间分布图

3.6 气候因素对植被指数影响综合分析

为进一步探究NDVI对不同气候影响因子的敏感性，本研究在降水、气温数据的基础上、添加对植被生长密切相关的VPD、辐射、CO2气象因子，首先将上述所有因子统一为空间范围一致、空间分辨率为2.5°的多波段时序数据、进行0-1之间的数据标准化（Raster Normalization扩展工具实现），具体实现过程参见图13。再编写IDL程序构建多元回归模型，模拟了黄河流域NDVI时序数据（受GOSAT-CO2数据可获取时间范围限制，选取2010-2019时间范围内的各时序数据，构建此模型）。

通过将原始时序数据与模拟时序数据逐像元对比分析可以看出(将两组NDVI时序数据导入Pearson Correlation Coefficient扩展工具实现线性相关性分析)，两者在空间格局上具有较好的一致性，且相关系数在0.8以上的区域占研究区的60.47%，而相关系数在0.5以下的像元呈零星分布状态（图25A-B）。总体而言，模拟数据与原始数据总体回归系数解释率为0.925，可以满足气候因子对植被覆盖变化影响的定量分析。

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

4 结论与讨论

基于MODIS NDVI数据和气象数据集，本文分析了2001-2020年黄河流域植被覆盖时空变化特征，并对其气候影响因子进行定量探究，主要得到以下结论： 

（1）基于高空间分辨率的土地利用与覆被产品数据，可对黄河流域呈现细碎斑块化分布的退耕还林还草、农用地侵蚀破坏等区域，进行有效地识别与提取。

（2）2001-2020年黄河流域植被覆盖呈显著增加趋势，平均增速为0.055·10 yr-1；且黄河流域植被恢复呈现出阶段性，2001-2010年的平均增速（0.067·10 yr-1）高于2011-2020年的平均增速（0.051·10 yr-1），整体上均通过置信度为0.05的显著性检验。空间格局上，NDVI趋势分析检测出的黄河流域植被覆盖呈极显著增长区域（p < 0.01）占整个研究区的51.03%的面积，且平均增速大于0.1·10 yr-1的地区占据研究区总面积的17.54%。

（3）黄河流域植被恢复过程受气候变化和人类活动共同作用。其中，受复杂因素驱动绿化面积占比为2.27%、受气候因子驱动绿化区域为39.18%、受生态建设驱动绿化区域占比为22.64%，其中NDVI受生态工程驱动绿化的区域与受气候因子驱动绿化区域，大体呈现以半湿润区/半干旱区为分界线的空间分布格局。

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

研究表明，黄河流域生态工程的实施已使得原先通过砍伐森林扩大耕地面积的现象已得到根本逆转，自然植被呈现显著恢复的趋势，但植被恢复趋势是否具有可持续性、生态工程的实施能否与地区经济发展相协同等方面，仍存在问题。与此同时，人类活动对于黄河流域植被生长积极/消极作用的良好区分、预测未来干旱等极端气候事件对于植被所产生的胁迫作用等问题，同样有待于进一步的研究改进与提升。本研究从宏观格局上揭示了植被对于气象因子变化的敏感性程度，对于深入认识植被对气候变化的响应过程及其规律，具有一定的意义。本研究存在的不确定性主要表现在：（1）空间分辨率为1km的NDVI数据对于评估黄河流域这一景观破碎化程度较高且区域内呈斑块状分布的草地精度较差，之后可尝试通过更高空间分辨率的NDVI数据进行评价。

（2）本研究在完成探讨不同气候因子对于植被生长的作用及定量化贡献率后，尚未将自然、人为驱动力分离开来，并与全球变化背景下频度和强度不断增加的极端天气气候事件（干旱、洪涝、高温热浪等）相结合，充分探讨上述事件对于研究区内的植被长势与时空变化分布格局所造成的影响。