2001—2020年青藏高原草地NDVI时空变化及驱动因子分析

吉珍霞，裴婷婷,,*，陈英,，侯青青，谢保鹏，吴华武

(1.甘肃农业大学资源与环境学院， 甘肃 兰州 730070； 2. 甘肃农业大学管理学院， 甘肃 兰州 730070； 3.甘肃农业大学草业学院，甘肃 兰州 730070； 4. 草业生态系统教育部重点实验室(甘肃农业大学)， 甘肃 兰州 730070； 5. 中国科学院南京地理与湖泊研究所,江苏 南京 210008)

全球变暖的背景下，陆地生态系统正在经受多环境变量的影响，而植被作为陆地生态系统的重要组成部分，可作为全球变化的“指示器”，对涵养水源、保持水土和调节气候有重要作用[1-2]。越来越多的证据表明，植被变化与环境要素变化严重影响植被的碳吸收，且植被发生变化时可能会通过一些生物-物理过程反馈区域气候[3-6]。因此，全面了解植被动态对环境变化的敏感性，对于生态系统碳平衡和草地资源的可持续利用至关重要。

草地是全球分布最广的生态系统之一，也是维持陆地生态系统稳定、保护生态环境和发展畜牧业经济的重要组成部分[7]。归一化差值植被指数(Normalized difference vegetation index，NDVI)是用于表征植被覆盖度和生长状况的度量指标之一，能够揭示环境的演化[8]。20世纪80年代以来，基于NDVI的植被动态变化已经被广泛的研究。如Eckert等[9]利用MODIS NDVI数据经时间序列分析发现NDVI适用于植被覆盖变化区域的监测和土地退化与再生的识别，且NDVI呈增加趋势和减少趋势的区域大多与植被增加或减少的土地覆盖等级变化的区域相吻合。田义超等[10]基于 2000—2011年SPOT-VEGETATION逐旬NDVI数据及逐日气温和降水数据，分析了北部湾沿海地区植被覆盖对温度和降水的旬响应特征，发现 2000—2011年植被处于恢复状态，且植被对降水和温度具有明显的阈值和滞后性。可见长时间序列植被数据中蕴藏了气候变化与环境的关系，开展大范围、长时间序列植被生长状况的研究具有重要意义。

全球水循环系统、水资源条件和生态系统发生了一系列的变化(如降水格局改变、草原灌丛化)[11]，这也导致草地生态系统植被的分布范围、生理生态特性、水分利用效率(Water use efficiency，WUE)发生了较大改变[12]。范广洲等[13]认为高原植被可以通过生理过程，产生净CO2吸收，从而降低大气CO2含量，起到缓解温室效应的作用。青藏高原是世界上海拔最高、中国面积最大的的高原，高寒草地大约占据了青藏高原面积的60%[14]。由于青藏高原有地形起伏较大、植被生长周期短和土壤贫瘠等特点，该地区的草地生态系统对环境变化格外敏感[15]。近些年研究表明，青藏高原年平均温度每年上升0.3～0.4℃，是全球同时期平均温度上升速率的2倍左右[16]。同时，青藏高原降水格局也在发生明显的改变，呈现北部增加、南部减少的趋势[17]。但在全球变暖的趋势中，气候变暖对青藏高原生态系统的影响可能存在差异。温度升高不仅会影响高原植被活动，而且还会促进表层土壤融化，造成土壤水分流失，从而加剧土壤干旱[18]。降水作为干旱半干旱地区草地生态系统水分补给的主要来源，有研究指出植被对水分的吸收在一定程度决定了自身的数量和分布格局[19]。当植被生产力受水热要素影响发生变化时，势必会造成整个生态系统功能的变化[20]。但目前对青藏高原草地的研究主要集中在草地生产力对平均气候变量(温度、降水)、二氧化碳和人类活动等的响应。相比以往对环境变量的研究，草地WUE和地表干旱程度直接作用于草地，对草地的生长也有着重要影响[21]。因此，探究温度、降水、草地WUE、地表干旱相互作用的过程，以及它们之间的交互作用对草地生长的影响，这对青藏高原地区植被变化机理的认识具有一定的生态价值和现实意义。

本研究基于2001—2020年植被NDVI数据，采用Theil-Sen中位数趋势分析、Mann-Kendall检验和变异系数分析法分析青藏高原草地NDVI的时空变化特征，运用岭回归分析和结构方程模型探究温度、降水、地表干旱和草地WUE对草地NDVI的影响，以描述该地区草地变化特征，并探明植被生长的主要驱动因素，旨在了解近年来青藏高原草地变化与环境要素变化的关系，为区域生态保护提供科学依据。

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

青藏高原地处中国西南部(26°～40° N和 73°～105° E)，总面积约257.24×104km2，被称为“世界屋脊”[22]。该地区属于高山高原气候区，其特点是温度低、日较差较大(图1C)。降水量自东南向西北呈递减趋势(图1E)，超过80.0%的年降水量发生在生长季(5月至9月)，且大部分地区WUE较低(图1D)。受印度洋季风和地形影响，青藏高原大部分地区属于干旱和半干旱区，仅中部和南部部分地区以半湿润和湿润为主[23](图1B)。青藏高原具有典型的水平地带性，受海拔、高原季风气候和高大山脉等的影响，其植被类型主要以草地为主(图1A)。

图1 青藏高原地理位置图Fig.1 Geographical location map of the Qinghai-Tibet Plateau注：TVDI为温度植被干旱指数(B)、TEMP为温度(C)、WUE为水分利用效率(D)、PRE为降水(E)，下同Note:TVDI is temperature vegetation drought index (B),TEMP is temperature (C),WUE is water use efficiency (D),and PRE is precipitation (E). The same as below

1.2 数据来源与处理

研究采用的遥感数据主要来源于美国国家航空航天局(National Aeronautics and Space Administration,NASA)提供的MOD13A2 NDVI，MOD11A2 LST(Land surface temperature),MOD17A2 GPP(Gross Primary Productivity)和MOD16A2 ET(Evapotranspiration)产品(https://lpdaacsvc.cr.usgs.gov)，时间分辨率分别为16 d和8 d，空间分辨率分别为500 m和1 km，时间尺度为2001—2020年。后经重采样和计算将所有数据空间分辨率统一为1 km，时间分辨率统一为生长季(每年的5月至9月)。为了减少云和大气对数据的影响，本文采取Savitzky-Golay(S-G)滤波对数据进行平滑处理[24]。另外，月数据采用了最大合成法，该方法可以进一步消除云和大气等的干扰，提高数据的精度。草地生长季年均WUE由生长季总GPP和生长季总ET比值计算，具体公式如下[25]：

(1)

式中：WUE为水分利用效率(g·m-2·mm-1)，GPP为总初级生产力(kg C·m-2)，ET为蒸散发(mm)。

温度植被干旱指数(Temperature vegetation drought index，TVDI)可根据LST和NDVI来计算，具体计算公式如下[26]：

(2)

Tsmax=a1+b1NDVI

(3)

Tsmin=a2+b2NDVI

(4)

式中：Ts为地表温度；Tsmax为NDVI对应的最高地表温度，即干边；Tsmin为NDVI对应的最小地表温度，即湿边；a1，b1，a2，b2分别为干、湿边线性拟合方程的系数。

90 m分辨率的数字高程模型(Digital elevation model，DEM)和1∶100万植被类型数据来源于中国科学院资源环境科学数据平台(http://www.resdc.cn)。在ArcGIS 10.2中经拼接、裁剪和重采样得到青藏高原1 km的DEM。植被类型数据经重分类得到森林、灌木、草地、农田和其他五种植被类型(图1A)，本文主要研究草地植被类型。

气候数据来源于中国气象数据网的日值数据集(http://data.cma.cn/data/)，时间跨度为2001—2020年。使用MATLAB R2020b软件计算温度、降水在生长季的时间序列，使用ANUSPLINE 4.3软件将该尺度数据插值为空间分辨率1 km的数据集。

1.3 研究方法

1.3.1趋势分析法 本文采用Theil-Sen中位数趋势分析和Mann-Kendall检验来研究生长季草地NDVI的空间趋势特征。Theil-Sen中位数趋势分析是一种稳健的非参数统计的趋势计算方法，对测量误差和离群数据不敏感，可用于生长季草地NDVI的时间序列趋势分析。而Mann-Kendall趋势检验作为非参数统计检验的方法，用于评估Theil-Sen斜率估计的显著性，即检验生长季草地NDVI趋势的显著性[25]。具体公式如下：

(5)

式中：β为斜率；i和j代表年份；如果β>0，NDVI数据集时间序列具有正趋势；如果β<0，NDVI数据集序列具有负趋势。

(6)

(7)

(8)

式中：xj和xi为NDVI时间序列数据；n为NDVI序列中数据个数；sgn()为符号函数；如果Zc≥Z1-α/2，则原假设被拒绝，NDVI序列呈显著变化趋势；如果Zc>0，则表明NDVI序列呈上升趋势；如果Zc<0，则表明NDVI序列呈下降趋势；α为给定的置信水平(显著性检验水平)，Zc≥1.96时表示NDVI序列通过了置信度为95%的显著性检验。

1.3.2变异系数 变异系数用来表示地理数据的波动程度[27]，在一定程度上可用于指示草地生长状况。其计算公式如下:

(9)

1.3.3岭回归分析法 由于岭回归可以消除自变量之间的共线性，又摒弃了最小二乘法的不偏性，因此作为一种被改进的最小二乘估计法被广泛运用于实际的回归过程[28]。本文利用岭回归来探究植被对环境因子的敏感性，以下是它的原理：

(10)

1.3.4结构方程模型 采用AMOS 23软件构建结构方程模型(Structural equation model，SEM)，该模型基于研究者的先验知识预先设定各因素之间的依赖关系，以此判断各因素之间关系强弱，可获得自变量(TVDI,WUE、温度和降水)对因变量(草地NDVI)影响的直接效果、间接效果和总效果，并能够拟合和判断整体模型，可更全面的了解研究系统[29]。本研究使用SEM模型中极大似然估计来研究草地NDVI,TVDI,WUE、温度和降水之间的直接影响和间接影响，构建影响草地生长的整体模型，对整体模型进行分析可得到各因子之间的路径系数[29-30]。通过SPSS 24中的因子分析方法对NDVI,WUE,TVDI、温度和降水进行KMO度量(KMO值域为0～1，越接近1，变量效果越好)，其中KMO值均大于0.75，表明各数据均呈正态分布，可进行建立结构方程模型。根据Amos 23软件构建SEM模型，计算模型适应性参数值，其中卡方自由度比为1.00(χ2/df∈[0，3])，近似误差均方根为0.04(RMSEA∈[0，0.05])，这些参数表明该模型配适度较好。

2 结果与分析2.1 青藏高原生长季草地NDVI的时空变化特征

由2001—2020年青藏高原生长季草地NDVI空间分布来看(图2A)，草地NDVI多年均值空间异质性明显，整体从东南到西北呈现逐渐下降趋势，多年均值为0.26。低值占比面积较大，主要分布在青藏高原西北地区(0.60)。2001—2020年生长季草地NDVI变化趋势如图2B所示(通过P<0.05的像元面积占比63.0%)，草地NDVI多年呈现稳定上升趋势，仅有西南部部分区域呈现下降趋势(面积占比19.0%)。将草地NDVI变异系数分为5个等级，分别为低波动区(0.3)。草地生长季NDVI的空间稳定性如图2C所示，青藏高原大部分地区的草地NDVI波动较低，在低波动区的面积占比66.2%，表明该地区的植被处于稳定生长的状态。高波动区和较高波动区零散的分布在整个区域，变异系数大于2的面积总占比为3.3%，表明这些地区的草地相比其他地区草地而言较为脆弱。2001—2020年青藏高原草地生长季NDVI年际变化如图3所示，区域内草地NDVI多年呈上升趋势，但多年斜率约为0.000 9。多年均值始终保持在0.25～0.28之间，草地NDVI最小年份为2003年，最大年份为2018年。

图2 2001—2020年青藏高原草地NDVI空间变化特征Fig.2 Spatial variation of grassland NDVI in the Qinghai-Tibet Plateau from 2001 to 2020注：A为草地NDVI多年均值，B为草地NDVI多年变化趋势，C为草地NDVI的变异系数Note:A is the multi-year mean value of grassland NDVI,B is the multi-year trend of grassland NDVI,and C is the variation coefficient of grassland NDVI

图3 2001—2020年青藏高原草地NDVI时间变化特征Fig.3 Temporal variation of grassland NDVI in the Qinghai-Tibet Plateau from 2001 to 2020

2.2 青藏高原地形对生长季草地NDVI的影响

将研究区海拔分为6个等级，分别为低海拔区(5 000 m)。植被NDVI在不同高程梯度中存在差异(图4A)，低海拔区、较高海拔区和极高海拔区的草地NDVI值主要集中在0.00～0.20之间，而较低海拔区、中海拔区和高海拔区的草地NDVI主要集中在0.20～0.60之间。将研究区坡度分为5个等级，分别为平坡地(0～5°)、缓坡地(5～8°)、斜坡地(8～15°)、陡坡地(15～25°)、急坡地(25～35°)。如图4B所示，草地NDVI值的分布随着坡度增大而越发离散，在坡度为0～25°区间内草地NDVI值主要集中在0～0.4。将研究区坡向分为5个等级，分别为平地(-1°)、阴坡(0～67.5°,337～360°)、半阴坡(67.5～112.5°,292.5～337.5°)、半阳坡(112.5～157.5°,247.5～292.5°)、阳坡(157.5～247.5°)。不同坡向草地NDVI存在显著差异(图4C)，草地NDVI值主要集中在0.30～0.60之间。青藏高原阴坡草地NDVI最集中，其次是阳坡、半阴坡和半阳坡。结果表明，海拔在3 000～4 000 m区间内的草地生长情况要好于其他海拔区间内的草地，坡度较小的区域植被分布较集中，阴坡草地生长情况较好。

图4 青藏高原草地NDVI与地形因子的相关性Fig.4 Correlation between grassland NDVI and topographic factors in the Qinghai-Tibet Plateau注：A为高程，B为坡度，C为坡向Note:A is elevation,B is slope,and C is aspect

2.3 青藏高原生长季草地NDVI对环境因子的敏感性响应

2001—2020年青藏高原草地生长季NDVI对WUE,TVDI、温度和降水敏感性空间分布如图5所示。草地NDVI对WUE的敏感性具有明显的空间异质性，且整个区域普遍呈现显著正敏感性，正敏感性高的区域主要分布在青藏高原东部(图5A)。草地NDVI对TVDI的敏感性正负差异明显，零散分布在整个青藏高原(图5B)。草地NDVI对温度的敏感性整体呈现出正敏感性，仅在西北部分地区呈现出高负敏感性(图5C)。草地NDVI对降水的敏感性空间分布规律较为明显，青藏高原西部敏感性较低，而东部敏感性较高，尤其在三江源地区(图5D)。结果表明，青藏高原整体草地受WUE升高和温度升高而生长情况变好，但在西南部部分地区的植被会随干旱程度加剧而退化。

图5 2001—2020年青藏高原植被NDVI对环境因子的敏感性Fig.5 Sensitivity of vegetation NDVI to environmental factors in the Qinghai-Tibet Plateau from 2001 to 2020

SEM模型表明青藏高原生长季草地与不同环境因素存在复杂的交互作用(图6)。不同因素均与草地NDVI呈显著正相关关系，其中草地WUE的正向影响最大(0.64)，其次是温度(0.35)、降水(0.07)和干旱(0.09)。此外，在整个相互作用路径中，温度(0.57)、降水(0.46)和干旱(0.29)对草地WUE的影响系数均高于温度、降水和干旱对草地NDVI的影响。由此可知，青藏高原地区草地生长的主导因子是WUE，而温度、降水和干旱作为草地NDVI变化的间接因子，主要通过影响草地WUE来影响草地生长。

对青藏高原生长季草地NDVI和WUE进行分段线性拟合，用于明确WUE对草地NDVI的影响，结果如图7所示。草地NDVI在WUE为3 g·m-2·mm-1时出现明显的转折，WUE在0～3 g·m-2·mm-1区间时草地NDVI呈现显著的上升趋势，上升斜率为0.28(P<0.01)。当WUE大于3 g·m-2·mm-1时，草地NDVI呈不显著下降趋势，下降斜率为0.03。该结果表明青藏高原草地在生长过程中对WUE存在阈值效应，即草地NDVI会随WUE增大而增大，但是超过了这一阈值，草地NDVI会随WUE增大而变小。

图6 青藏高原草地生长季NDVI与环境要素的结构方程模型Fig.6 Structural equation model of grassland NDVI and environmental factors during the growing season in the Qinghai-Tibet Plateau注：*代表显著性0.001。单向箭头表示由自变量指向因变量。双向箭头表示两因子之间的相关关系。线上的数值表示标准化路径系数，表示自变量对因变量影响的大小。数值正负表示相关关系的正负Note:* indicate significant difference at the 0.001 level. A one-way arrow points from an independent variable to a dependent variable. The bidirectional arrows indicate the correlation between two factors. The value on the line represents the standardized path coefficient and represents the influence of the independent variable on the dependent variable. Numerical values are positive or negative of the correlation

图7 青藏高原生长季草地NDVI与草地WUE之间的相关关系Fig.7 Correlation between grassland NDVI and grassland WUE in the growing season of Qinghai-Tibet Plateau注：*表示0.05显著性差异Note:* Significant difference at the 0.05 level

3 讨论

本研究发现青藏高原整体草地呈稳定生长状态。其中，张戈丽等[31]利用GIMMS NDVI数据研究1982—2006年青藏高原植被覆盖变化时指出青藏高原地区植被覆盖整体表现出“整体升高，局部退化”的趋势。李达等[32]基于MOD09A1数据反演青藏高原生长季植被NDVI分析其时空特征和空间差异时，发现从2001—2018年青藏高原生长季NDVI呈缓慢上升趋势，且空间差异明显。Hao等[33]通过耦合季节趋势模型和增强回归树模型分析青藏高原植被的变化趋势，结果发现青藏高原草地生长季NDVI整体呈现上升的趋势，且多年绿化趋势显著，这些结果均与本文的研究一致。另外，本文从地形角度分析了青藏高原生长季草地的分布特征，结果发现不同海拔、不同坡度和不同坡向间植被分布存在显著差异。基于以往研究发现，青藏高原海拔可以调节高寒草地物候变化、草地生产力及其对气候变暖的响应[34-35]，且因不同海拔、坡度和坡向内气候、太阳辐射、日照时数等的差异，造成不同地形条件下草地生长的差异[36]。

青藏高原是一个由连绵的山脉和丘陵组成的大型山地系统，山坡和海拔创造了一个明显的温度和水分梯度(图1)。而青藏高原草地是对全球气候变化最敏感的生态系统之一，气候因子的空间变化趋势，以及青藏高原及其周边地区生态系统服务功能和水资源的变化，可能会导致高寒草地的空间异质性[37]。这一结果经本文证实，即温度、降水、WUE和TVDI被认为是影响青藏高原草地NDVI时空分异的重要因素。且经岭回归分析发现2001—2020年青藏高原大部分区域草地会受WUE升高和温度升高的影响而生长情况变好，而部分地区植被会受干旱程度加剧而退化。可能的原因是温度升高可以通过增强植物光合作用或促进矿质营养物质的吸收来促进有机质的产生，从而有利于植物的生长发育[34]。但是随着温度增加，草地干旱程度加剧，水分亏缺，也会限制部分植被的生长，况且青藏高原属于干旱半干旱地区，植被的有效用水主要来源于降水。Zhao等[38]通过研究干旱对植被生产力的影响时得出植被生产力对温度的响应随时间增长是由水分主导的，这也充分证实了水分和温度对高寒干旱半干旱地区草地生长的积极促进作用。

本文与以往研究不同的是，将WUE作为衡量草地生长的指标之一，通过结构方程模型对草地生长过程进行归因分析后发现，温度、降水和TVDI主要是通过影响植被WUE来影响植被生长。生态系统碳、水循环过程是一个极其复杂的反馈过程，在植被内部生理生态特征尺度上，植物光合作用和蒸腾作用主要受气孔导度的影响[39]。孙双峰等[40]指出植被会在不同的生长阶段根据外界环境来调整气孔导度，调整自身的水分利用效率，保证植被的正常生长。如干旱时植物自身能够缓解蒸腾耗水，增加碳吸收量[41]。此外，本文发现青藏高原草地NDVI对WUE存在阈值效应，原因可能是温度升高可以使叶片气孔导度加大，而光合作用增加幅度大于蒸散作用幅度时会促使WUE升高，但当温度过高，土壤发生水分亏缺时植被WUE下降，尤其是干旱区一年生的植被[42]。此外，本研究不足的是未考虑人类活动对青藏高原草地变化的影响，在未来的研究过程中，需进一步考虑政府政策和人类活动对青藏高原植被的影响。

4 结论

本研究通过详细分析青藏高原草地的时空变化特征及其对多影响因素的响应，发现青藏高原草地NDVI在坡度较小的区域植被分布较集中，阴坡草地生长情况较好，且青藏高原草地整体呈稳定生长的状态，其受WUE升高和温度升高而生长情况变好。影响青藏高原草地NDVI变化的主导因子是WUE，草地NDVI对WUE存在阈值效应，即草地NDVI会随WUE增大而增大，当WUE大于3 g·m-2·mm-1时，NDVI会随WUE增大而变小，该结论可为全球变暖背景下植被变化提供新的见解。