植被是全球陆地生态系统的重要组成部分[1], 在调节气候变化、维持陆地生态系统功能等方面发挥着重要作用[2-3]。青藏高原作为全球平均海拔最高的独特自然地理单元, 是我国乃至亚洲重要的生态安全屏障[4], 素有“世界第三极”之称, 也是全球气候变化的敏感区和生态脆弱带[5-6]。在过去几十年, 随着全球气候变化与人类活动加剧, 青藏高原地表自然地理要素经历了复杂的变化过程[7-8], 这种变化不仅显著影响了当地的植被生态系统[1, 9], 也对毗邻周边、乃至全球变化产生深刻的影响[10]。从20世纪80年代以来, 青藏高原一直是气候-植被关系及生态系统变化研究的热点和难点区域[11-12], 被视为我国乃至全球气候变化的“天然实验室”[13]。因此, 深入探讨全球变化背景下高原植被的时空演变机制, 理解高原生态系统对气候变化的响应和适应机理, 对于应对区域和全球气候变化、提升和优化高原生态功能和全球生态建设都具有重要意义。

一般而言, 区域植被的短期变化多是受人类活动和自然干扰的影响, 而长期变化则取决于温度、降水、辐射等水热条件的影响[14]。在全球气候变化的背景下, 植被生长受水热条件胁迫的影响日趋明显[15], 然而, 植被对水热变化的响应由于区域生态系统的异质性而呈现明显的空间差异[16]。目前国内外学者主要基于AVHRR、MODIS、Landsat等卫星遥感数据对整个青藏高原或高原局地区域开展了一系列植被动态研究[17-19], 主要包括植被时空变化格局[20]、演变规律[9]及其驱动机制[21]等内容。研究结果显示, 由于研究区域、时空尺度的不同, 青藏高原的植被对水热条件变化的响应机制也明显的不同[15, 22]。近40年来, 青藏高原植被呈现“整体增加、局部减少”的趋势[14-15], 但是不同子区域植被动态机制存在显著差异[14, 23], 水热条件缺乏区域植被显著增加, 而水热条件充足地区植被却表现为减少趋势[22]。大部分研究认为, 青藏高原植被的变化主要受气温影响[24-25], 三江源等局部区域植被对气温的响应尤为明显[26]。也有学者认为, 青藏高原植被变化受降水[18]或辐射[27]的影响较气温密切。而最新研究结果显示, 不同类型植被受不同水热条件的影响, 如高寒草甸主要受温度影响, 高寒草原受降水和温度双重影响[28], 温性草地植被生长受地表温度的影响[21]。

综上所述, 以往对青藏高原植被动态的研究基于不同数据源、时间尺度和地理区域, 在植被时空变化规律、植被对水热条件变化的响应方向和响应程度, 尤其是长时间尺度上的响应过程还存在着不确定性, 甚至争议。为此, 本文基于近40年的长时间序列植被数据, 从不同时空尺度上探究不同类型植被的变化规律及其对水热条件的响应, 以期对青藏高原植被-气候变化之间的关系有更深一步的了解, 为气候变化背景下青藏高原植被动态监测、生态修复、生态资源的合理利用和保护等方面提供参考。

青藏高原面积约为257万km2[29], 平均海拔超过4000 m, 年平均气温为-6—20℃, 年降水量50—2000 mm, 太阳辐射量达0.5861—0.7954 MJ cm-1 a-1, 是我国太阳总辐射值最高的地区[30]。由于水热条件在水平和垂直方向上的区域差异, 不同气候梯度的热量和水分约束下, 青藏高原的植被群落由东南的亚热带雨林过渡到西北的高寒荒漠[23, 31]。中国植被区划将青藏高原划分为11个不同的植被类型区(图 1)。植被类型分区在空间上是呈现出完整的、连续的植被带状分异性[33], 可以指示植被地理分布的经纬度和垂直地带性规律及其与水热条件的关系。

AVHRR (Advanced Very High Resolution Radiometer) NDVI (Normalized Difference Vegetation Index) 数据来源于美国国家航空和航天局(NASA)的全球监测与模型研究组GIMMS (Global Inventor Modeling and Mapping Studies) 发布的半月最大值合成产品GIMMS NDVI3g(以下简称GIMMS NDVI), 空间分辨率为0.083°, 时间跨度为1982—2015年。MODIS (Moderate Resolution Imaging Spectroradiometer) NDVI数据亦来源于NASA陆地产品组MODIS发布的16天的合成植被指数产品MOD13C1, 其空间分辨率为0.05°, 时间跨度为2000—2020年。气温、降水和太阳辐射数据(1979—2018年)来源于中国气象科学数据网(http://data.cma.cn), 该数据是基于全国2400多个气象站点日观测数据, 利用澳大利亚ANUSPLIN插值软件插值处理生成, 空间分辨率0.01°。数据高程模型(DEM)来源于中国科学院资源环境科学数据中心(http://www.resdc.cn), 空间分辨率为1km。植被类型数据为1 ∶ 100万中国植被类型图。

由于AVHRR和MODIS两种数据的传感器及时空分辨率的不同, 造成两种数据的NDVI值之间存在一定的差异[34]。为降低这种差异, 同时保持数据在空间上的一致性, 首先利用Savitzky-Golay滤波对青藏高原1982—2020年GIMMS NDVI、MODIS NDVI数据进行滤波去噪, 采用最大值合成(Maximum Value Composition, MVC)计算生长季(5—9月)NDVI, 并将MODIS NDVI采样到0.083°；对两组数据集重叠观测年份(2000—2015年)的NDVI逐像元建立数据的映射关系；最终, 通过逐像元回归分析将NDVI时序数据从1982年延伸至2020年(图 2), 用其代表研究区植被的年际变化。GIMMS NDVI和MODIS NDVI数据在重叠年份的年际变化趋势是一致的, 相关系数介于0.88—0.93之间, 均通过了0.001的置信度检验, 适用于分析青藏高原植被覆盖的时空变化格局。为了消除非植被因素对研究结果的影响, 本文取NDVI≥0.1的像元进行分析[35]。

Mann-Kendall(MK)检验法，是一种非参数的统计方法。由于该方法不需要数据服从特定的分布特征，且受异常值的影响较小[36-37]，能更好地监测轻微的变化区域，变化趋势更为准确[38]。因此本文利用MK趋势检验来计算长时间序列的植被指数变化趋势。计算方法如公式(1)—(4)：

式中，Zmk标准化检验统计量，S是由1982—2020年的NDVI计算得到的近似服从正态分布的趋势统计量，sgn(θ)为逻辑判别函数，var(S)为S统计量的方差。本文中的检验取显著性水平α=0.01，在Z1-α=2.58上进行双边趋势检验。

植被生长通常受到多个水热因子共同作用的影响，相关分析法已广泛应用于探索植被生长与水热因子之间的联系[16]。本文分别采用偏相关分析和复相关分析法逐像元对青藏高原生长季最大合成NDVI与生长季月平均温度、累计降水量、太阳辐射量进行相关分析，其中偏相关系数用以表示植被与单个水热因子之间的相关程度，复相关系数表示植被与多个水热因子之间的复合关系程度。偏向关系数和复相关系数计算方法为公式(5)—(7)：

式中，Rxy为植被与单个水热因子间的简单相关系数；Ry1·23为表示固定水热因子2、3之后NDVI(y)与水热因子1的偏相关系数；Ry·123为NDVI(y)与水热因子1、2、3之间的复相关系数。最后，采用t检验方法完成偏相关系数的显著性检验，采用F检验完成复相关系数的显著性检验。

如图 3所示青藏高原1982—2020年近40年NDVI均值在空间分布上具有较强的规律性, 在区域水热条件的影响下, 青藏高原植被整体上呈现由东南-西北梯度递减的分布特征。青藏高原植被NDVI最高的区域位于东南部等, 以草甸、草原、森林为主, NDVI最高达0.92, 青藏高原主要的农耕地及人工园林也分布于此(图 1)。植被NDVI最低的区域则集中于西北部和北部边界的高寒荒漠地带, 这与高原上的地形、气候分异和植被类型相关。

由海拔梯度分布特征统计结果(图 3)看, 青藏高原的植被覆盖呈现明显的垂直地带性, 随着海拔的升高, 植被NDVI逐渐减小。海拔 < 2500m的高原东南缘, 面积仅占整个高原的2.37%, NDVI均值达到0.90；2500—3500m的中海拔带既有植被覆盖极低(NDVI < 0.1)的北部沙漠区域, 也有以山区林地和农耕地为主的东部主要人类活动区(NDVI>0.8), 面积占14.37%；3500—4500m的中高海拔带, 是青藏高原主要的牧场所在地, 植被覆盖较高的高寒草甸、草原均分布于此, 面积占比27.57%, NDVI均值达0.55；4500—5500m的高海拔带主要为高原西北部的高寒草原及高山植被区, 面积超过高原总面积的一半以上(51.79%), 而NDVI均值仅为0.33；>5000m的极高海拔带则基本无植被分布。

1982—2020年青藏高原植被生长季NDVI以平均0.006/10a(P < 0.01)的速度呈增加趋势(图 4)。MK趋势分析结果表明, 植被NDVI的变化在不同的植被类型分区中有所不同, 从高寒荒漠区的0.014/10a降低到高寒草甸区的-0.007/10a。其中显著增加的区域主要集中分布在高原北部、西部和南部, 面积约为183.47万km2, 占总面积的73.97%；显著减少的区域主要分布在高原腹地及海拔相对较低、人类活动频繁的高原东南部地区, 面积约为45.59万km2, 占比18.38%；植被基本保持不变的区域分散在高原各处, 东南部居多。此外, 在海拔梯度上, 除 < 2500m的海拔带外, 植被显著增加的区域面积均大于减少的区域面积。

青藏高原植被以草原、草甸、高山植被、灌丛、荒漠、针叶林和阔叶林为主, 占青藏高原总面积的90%以上[11]。1982—2020年间, 青藏高原所有类型的植被生长季NDVI均呈现增加的趋势(图 5), 这与Zhou等[2]的研究结果一致。其中荒漠、草原、高山植被NDVI年际变化为极显著增加(P < 0.01), 增加面积分别达到19.2万km2(7.74%)、56.60万km2(22.92%)、18.65万km2(7.52%)；森林和灌草丛的NDVI年际变化在2000年前后呈现先波动后趋于平稳的增加趋势, 这与NDVI数值在高植被密度区敏感性降低有较大关系, 而在2000年后不同传感器的数据融合弱化了这种波动使得年际变化趋于平稳增加, 增速分别达到0.009/10a(P < 0.01)和0.005/10a(P < 0.01), 显著增加面积分别为13.88万km2(5.6%)和17.94万km2(7.23%)；年际变化上只有高寒草甸NDVI呈现不显著地增加(0.0002/10a), 但增加的面积达到42.31万km2(17.06%)。

从1982—2020年的四个年代际尺度上分析高原植被变化的变化特征(图 6)。1982—1990年, 青藏高原88.86%的区域植被NDVI呈现不显著的增加, 仅有高原西南边缘和柴达木盆地西南约9.01%的区域呈显著增加趋势。显著减少的区域分散在一江两河、河湟谷地等主要人类开发聚集区域。1991—2000年间, 高原大部分区域(84.1%)的植被变化依然为不显著增加。显著增加的区域主要聚集于黄河、长江、雅鲁藏布江等水源的区域及川西高原, 而显著减少的区域则分布于高原北部荒漠地带, 说明该阶段植被主要受水分影响。2001—2010年, 随着三江源生态保护、禁牧、退耕还草等生态工程的实施, 生态正效应明显。三江源、祁连山及高原东部区域植被显著增加, 面积占比达到21.45%。显著减少的区域主要在高原西南、西藏境内, 这与卓嘎等[18]、张镱锂等[39]和王敏等[40]的研究结果一致。2011—2020年, 整个高原面上64.79%的区域植被NDVI表现为不显著的增加, 明显增加的区域集中于高原北部的祁连山、长江与澜沧江水源区；显著减少的区域零星分布在柴达木盆地、横断山脉、喜马拉雅山脉沿线。

以8km的方格网为单元, 建立1982—2018年青藏高原生长季平均NDVI与平均气温、累计降水量、以及辐射量之间的相关性散点图(图 7), 从静态角度探讨植被对水热因子的响应。发现NDVI值与气温、降水和辐射量之间均为非线性关系。根据NDVI与气温二次函数关系(R2=0.202, P < 0.001), 取气温阈值区间2—11℃, 在该气温区间内, NDVI随气温升高而迅速增加, 增速达0.026/℃, 当气温超过11℃时, NDVI增加不再显著。从空间分布上看, 青藏高原大部分植被均分布2—11℃的温度区域内。同样, 根据NDVI与降水的二次函数关系(R2=0.693, P＜0.001), 取阈值区间250—650mm, 在该降水量区间内, NDVI随降水量增加而迅速升高, 增速为0.001/(mm), 高原腹地的草原、草甸等植被均分布在该降水区间内, 而在阈值区间外的NDVI与降水量关系则不显著。根据NDVI与太阳辐射的二次函数关系(R2=0.649, P < 0.001), 取阈值区间2900—3750 MJ/m2, 在该太阳辐射区间内, NDVI随着辐射量的增加而迅速减少, 减少速率为0.0007/(MJ/m2), 表明过高的太阳辐射抑制植被生长, 在辐射量>3350 MJ/m2高原西部高海拔地带, 较少有植被分布。

为探讨植被对单个水热条件因子变化的动态响应, 逐像元计算了生长季NDVI与气温、降水、辐射量之间的偏相关关系(图 8)。从空间分布看, 1982—2018年青藏高原生长季气温在大部分区域为增加趋势, NDVI与气温的偏相关关系由东南向西北呈现中部高、外围低的“夹心”状特征, 其中显著正相关(P < 0.01)区域集中于青海境内的柴达木盆地、羌塘高原、及西藏西南部；显著负相关(P < 0.01)的区域在高原西北部、东南部的横断山脉等均有分布。降水量变化虽在整个高原上呈现不显著增加趋势, 但高原南部的大部分区域表现为降水减少的趋势。NDVI与降水偏相关系数呈东西向的干湿地带性分布, 并由北向南递减, 除森林、灌草丛植被区域外, 高原上大部分植被类型区的NDVI都与降水呈正相关, 其中北部的高山植被、草原、草甸等植被类型区NDVI与降水显著正相关(P < 0.01)。太阳辐射除了在喀喇昆仑山西北部、横断山脉西南部的分布地区外, 其他大部分区域为减少趋势, 而辐射显著减少的区域也正是与NDVI呈现显著的负相关关系(P < 0.01)的区域。

从不同类型植被(图 9)及区划(图 10)与水热因子的相关程度看, 青藏高原上除森林和灌草丛区(IVA、IVB)外, 各类型植被与降水的相关程度要好于气温和辐射。降水对高原不同植被类型的正效应大小排序为：荒漠＞草原＞草甸>高山植被＞森林＞灌丛＞栽培植被, 说明降水是荒漠、草原、草甸的主要促进因子。温度对高原不同植被类型负效应大小排序是：栽培植被>高寒植被＞荒漠＞灌丛>森林>草甸＞草原, 气温是荒漠、高寒植被的主要抑制因子。辐射量与不同植被类型的相关系数大小排序是：森林＞灌丛＞高山植被＞栽培植被＞草甸>荒漠＞草原, 辐射对森林、灌丛和高山植被的影响为正效应, 而对荒漠、草原影响为负效应。植被类型相同的不同分区间, 由于地形地貌、海拔、经纬度、土壤类型等自然环境的不同存在显著的差异, 如IIIA1柴达木盆地高寒荒漠区、IA1昆仑山高寒荒漠区、及IA2阿里高寒荒漠区都是荒漠类型, 但由于三个区域地处不同的纬度和海拔带, 地形地貌也显著不同, 因此植被与水热条件的相关关系有明显的差异。

根据以上的分析发现, 青藏高原地表植被变化的水热因子驱动存在着明显的区域分异, 大部分区域植被变化受单一因子(气温或降水或辐射)的驱动较为显著, 部分区域受多个因子的共同作用。对植被变化的驱动因子进行分区, 才可以进一步的概括植被在区域尺度上的变化规律。因此, 本文通过分析青藏高原多年NDVI与气温、降水、以及辐射量的偏相关和复相关关系(图 11), 同时为保证变化特征的最大一致性提取其显著性水平为0.01的区域, 按表 1规则进行水热驱动因子分区研究。

整体而言, 在全球变暖和青藏高原向暖湿化演变的过程中, 大部分区域植被变化受气候因素影响, 呈现气候驱动型, 不同类型植被生长的主要驱动因子与植被的生态特征及气候环境密切相关。其中, 高原北部高寒植被生长主要受气温和太阳辐射影响, 高原腹地干旱半干旱的草原主要受降水影响, 高原中南部草甸、灌丛区则主要受气温驱动, 对于高原南部边缘的森林生长, 受辐射和温度的影响多一点, 但则没有明显的规律性。这一结果与Wu D等[16]、Li等[31]和Jong等[41]的研究具有较好的一致性。

基于1982—2020的时序NDVI数据发现, 在过去的近40年中, 青藏高原73.97%的区域呈现绿化趋势, 且在大多数的植被类型区, 生长季NDVI的变化速率都是正值, 平均增速达0.006/10a。其中2000年以前增速为0.005/10a, 而2000年以后增速达0.011/10a, 植被迅速增加, 这与国家在地区实行“退耕还林”与“退牧还草”等生态政策密切相关。植被变化的这一增长趋势与杨元合等[30](1982—1999年)、与Huang等[27](2001—2017年)的结果近似或一致, 而与韩炳宏等[24]认为的2000—2009年呈减少态势, 2010—2018年呈增加态势有所不同, 这种差异可能与趋势分时所选取的数据起始年份有一定的关系。从海拔梯度上的变化(图 3)可以看出, 整个高原除 < 2500m的海拔带外, 其他海拔带上植被NDVI变化速率也均为正值, 其中2500—3500m的海拔带为青藏高原主要人类活动区, 植被类型以草甸、灌丛、栽培植被为主, 植被变化受人类活动(耕作制度、灌溉及施肥方式、退耕还林等相关政策)影响较大；在3500—5500m海拔带, 植被NDVI增加的面积随海拔升高而不断增加, 这一海拔高度的植被主要受气温、降水等水热条件影响, 受人为干扰较少, 在全球变暖的大背景下, 气候变化对该区植被产生正效应。

在气候变化影响下, 植被总是在不断地区适应外界条件的改变使其自身生长更为有利, 该过程是动态的、非线性的, 植被活动随着水热气候因子的变化过程存在明显的阈值[42]。当温度、降水等水热条件的改变超过植被的适应能力时, 植被将会在结构功能等多方面受到抑制作用[46]。如图 8所示, 柴达木盆地、羌塘高原、及西藏西南部等高海拔区域, 植被增加与气候变暖有显著的正相关；而高原西北部、东南部的横断山脉等区域为稀疏的高山植被或荒漠, 升温加速土壤水分蒸发, 抑制了植被生长；东南部的灌木及森林地带的升温加速植被蒸腾, 光合作用减弱[46], 温度成为影响这些区域植被生长的限制性因子。NDVI与降水量显著正相关的高原中北部为高山植被、草原、草甸等低矮植被类型, 降水量增加有利土壤水分和养分的运转、植被吸收和生长。虽然西藏东部与横断山脉的森林和灌木区降水量较大, 但这些区域的NDVI与降水量呈现不显著负相关关系, 而与辐射量显著正相关, 说明降水量不构成植被生长的限制因子, 但由于降水量大, 云量较大, 导致辐射减弱, 不利于植被蒸腾和光合作用, 辐射成为影响这些区域植被生长的限制因子。NDVI与辐射负相关的地区如柴达木盆地, 地处干旱区, 海拔低、辐射强、地下水出露、加之人类活动等各种因素, 植被的生长规律与其他荒漠区截然不同[47], 但这一过程的机制机理还不得而知, 需要另行研究。

长时间序列的植被变化是一个复杂的过程, 不同的植被类型对水热条件和人类活动的干扰存在不同的响应, 而本文只研究了青藏高原生长季植被在长时间尺度上对的水热条件响应的规律, 并未深入植被对水热条件变化的时滞性不同。同时, 对于如柴达木盆地、青海南部高原等特殊高寒区脆弱自然环境下植被对水热响应表现出的不同现象并未深入研究。由于青藏高原的地理特殊性, 植被的生长不仅与气温、降水、辐射等水热条件有密切的关系, 同时也与冰川冻土、土壤水热、植被蒸散等多种因素有复杂的相互关系, 这些都将是未来探索的方向。

(1) 1982—2020年青藏高原植被NDVI在年际上呈波动增长趋势(0.006/10a)。植被显著增加的区域占高原总面积的73.97%, 主要分布在青藏高原北部、西部及南部边缘, 植被减少的区域主要为高原东部及南部的人类活动区。在年代际上, 2001—2010年三江源等区域的生态保护工程正效应明显。

(2) 青藏高原植被对水热气候条件的响应具有显著的空间异质性。静态上, 植被对气温的显著响应区间为2—11℃, 对降水的响应区间为250—650mm, 对辐射的响应区间为2900—3750MJ/m2, 空间上表现为高原腹地响应明显, 东南部四川盆地、西北部昆仑山脉和柴达木盆地无明显响应的一致性；动态上, 植被对气温、降水和辐射的响应在不同植被类型中，其响应关系、方向、程度均存在差异, 而除森林和灌丛类型外, 植被与降水的偏相关程度要优于温度和辐射。

(3) 青藏高原植被生长受水热因子驱动的区域占高原总面的55.95%, 其中42.72%的区域由气温、降水和辐射互补驱动, 首要驱动因子是温度, 占比26.54%, 13.23%的区域由多因子联合驱动；非气候驱动区域占高原总面积的44.05%。