①该作者生于1999年。

②该作者生于1997年， 与①处作者同名不同人。

全球环境变化是当前国内外学者关注的热点问题， 其中全球气候变暖问题是焦点之一。受全球及区域气候变暖影响， 青藏高原作为除南极、 北极外的“第三极”， 其地表自然地理要素正在经历复杂的变化过程^［1-3］。植被作为地表物质与能量循环的首要环节， 深入探讨气候变暖背景下青藏高原植被时空变化及响应过程， 对于反映和认识全球及区域气候变化的陆地植被生态系统响应具有重要意义。

归一化植被指数（Normalized Difference Vegetation Index， NDVI）是利用遥感技术获取的具有高时间分辨率、 多空间分辨率的数据集， 为开展全球及区域植被生态系统变化研究提供了有力支撑^［4-5］。根据NDVI的已有研究， 主要围绕4个方面开展， 分别是基于NDVI的植被生态系统时空变化过程及其对气候变化的响应、 植被生态系统对气候变化及人类活动响应的差异研究、 基于NDVI的植被生态系统功能变化尤其是植被碳汇功能、 基于NDVI的青藏高原植被物候变化研究。第一个研究方面， 主要基于GIMMS NDVI（1982 - 2006年或1982 - 2013年）、 MODIS NDVI（2000年以来）、 SPOT-NDVI（1999年以来）等开展了青藏高原^［6-8］及三江源地区^［9］、 雅鲁藏布江地区^［10］、 西藏自治区^［11］和青海省^［12］的植被时空变化研究。该类研究中， 由于使用的数据源具有一定的时、 空尺度差异， 导致同时段的研究结论相似性高于差异性， 而不同时段研究结论差异性较大。由于数据更新的先后， 同类数据产生了不同的版本， 如GIMMS NDVIg和NDVI3g数据， 杜加强等^［13］对比分析认为NDVI3g数据NDVI均值较NDVIg、 MODIS NDVI和Landsat NDVI高。第二个研究方面， 已有研究借助NDVI、 气象数据、 生态保护政策及重大工程实施等内容， 定量区分了青藏高原植被生态系统变化的气候变化和人类活动因素^［14-17］， 但在结论方面存在差异， 部分研究认为气候变化是青藏高原植被生态系统的主控因素， 也有学者认为气候变化和人类活动对植被生态系统的控制存在此消彼长的过程， 具有时间变化阶段性和空间分布的异质性等特点^［18-20］。第三个研究方面， 以NDVI为基础， 探讨了青藏高原不同地区、 不同植被类型净初级生产力（Net Primary Productivity， NPP）及碳储量/碳汇的时空变化及其与气候、 人类活动的关系。研究认为随着全球气候变暖、 生态安全屏障保护与建设工程的实施， 青藏高原植被生态系统总体趋好， 局部变差， 植被碳储量/碳汇水平提高， 植被生态系统服务功能得到提升^［21-24］。第四个研究方面， 青藏高原植被物候的返青期提前， 生长期延长， 对于高原植被生态而言， 有利于植物的生长和发育以及植被碳储量的增加^［25-27］。

已有研究总体表明， 青藏高原自20世纪80年代以来， 植被状况总体向好， 局部变差。气候变暖引起的冰川融水增加， 湖泊面积扩大， 冻土面积萎缩， 土壤温度升高等， 均为高原植被生长和发育提供了有利条件， 植被生态系统服务功能、 植被碳储量/碳汇进一步增强等。尽管青藏高原植被生态系统研究较为丰富和全面， 但在植被生态系统尤其是植被类型变化特征及气候响应的空间差异方面还有不完善之处， 如已有研究多使用1982 - 2006年GIMMS NDVIg或2000年以来 MODIS NDVI及二者相结合方法进行分析^［28-30］， 而最新推出的1982 - 2015年GIMMS NDVI3g数据直接提供了长时间序列数据。本研究即以1982 - 2015年GIMMS NDVI、 气温、 降水、 植被类型数据为基础， 利用一元线性回归模型、 相关系数、 偏相关系数及t检验方法， 分析青藏高原植被生态系统及不同植被类型NDVI变化特征及其气候响应， 以期为青藏高原生态环境保护与建设提供科学依据。

青藏高原绝大部分位于位于中国西部， 介于74°27′ ~ 104°59′ E和27°12′ ~ 38°20′ N之间， 总面积约259.82 × 10⁴ km²。青藏高原南起中国边境线的喜马拉雅山南缘， 北至昆仑山、 阿尔金山、 祁连山北缘， 西部为帕米尔高原和喀喇昆仑山， 东与秦岭山脉西段和黄土高原相连， 海拔介于3 000 ~ 5 000 m之间， 平均海拔4 000 m以上， 被称为“世界屋脊”。青藏高原可划分为藏北高原、 藏南谷地、 柴达木盆地、 祁连山地、 青海高原、 川藏高山峡谷区6个部分。

青藏高原气候为独特的高原气候， 辐射强烈， 日照时间长、 平均气温低、 年积温少， 气温日较差大， 干湿分明， 冬季干冷漫长， 夏季温凉多雨。年平均气温由南而北由20 ℃降至 -6 ℃以下， 年降水量由2 000 mm递减至50 mm以下。青藏高原是东亚、 东南亚和南亚众多河流的发源地， 如黄河、 长江、 澜沧江、 雅鲁藏布江等。此外， 青藏高原也是中国湖泊分布较为密集的区域， 主要湖泊有青海湖、 色林错、 纳木错、 鄂陵湖、 班公湖等。

青藏高原植被类型多样， 包括针叶林、 针阔混交林、 阔叶林、 灌丛、 荒漠、 草原、 草丛、 草甸、 沼泽、 高山植被、 栽培植被及其他植被， 其中以草原植被覆盖面积最广， 占青藏高原区域总面积的27.76%， 其次为草甸， 占25.32%， 高山植被和灌丛分别占11.26%和10.35%， 其他植被类型覆盖面积占比均低于10%（图1）。随着全球气候变暖， 冰川加速融化， 高原湖泊水量增加， 湖水位上升， 青藏高原脆弱的陆地植被态系统如何变化受到广泛关注。

研究使用数据包括： （1）GIMMS NDVI3g（1982 - 2015年）， GIMMS （Global Inventory Modelling And Mapping Studies）NDVI数据是美国国家航天航空局（NASA）提供的全球植被指数数据， 下载自https： //nex.nasa.gov/nex网站， 为3g.v1版本， 数据为ncf格式压缩文件， 通过MATLAB对其进行解压、 投影变换、 重采样等， 得到时间分辨率为15 d、 空间分辨率为5 km的GIMMS NDVI数据。对15 d数据采用平均值方法， 计算得到月值数据， 月值数据采用最大值合成（Maximum Value Composition， MVC）计算得到1982 - 2015年年值数据。年值数据重采样为1 km分辨率， 用于时空变化过程及其与气温降水相关关系分析。（2）气温、 降水空间插值数据（1982 - 2015年）， 该数据是基于全国2 400多个气象站点日观测数据， 利用澳大利亚ANUSPLIN插值软件插值处理生成， 数据分辨率为1 km。（3）植被类型数据， 基于2001年植被生态学家侯学煜院士主编的《1∶1 000 000中国植被图集》数字化后得到， 分辨率为1 km。（4）青藏高原矢量边界数据。以上数据中， （2）、 （3）和（4）数据来源于中国科学院资源环境科学数据中心（http： //www.resdc.cn）， 利用青藏高原矢量边界， 裁切得到1982 - 2015年GIMMS NDVI逐年数据、 气温和降水逐年插值数据、 植被类型数据。

以NDVI为因变量y， 年份为自变量x， 利用最小二乘法计算得到NDVI与年份之间的一元线性模型， 计算公式为^［31］

式中： a为系数（变化速率）； b为常数； y¯和x¯分别为1982 - 2015年NDVI和年份的平均值。一元线性模型得到后， 利用F检验进行显著性检验^［31］。

式中： y'为利用一元线性回归模型计算得到的y的拟合值。计算得到F值后， 在F分布表中查找α=0.05时， F的临界值4.15， 若F>4.15， 则回归方程是显著的， 相反， 则不显著。

相关分析用于分析NDVI与气温、 降水关系。其中相关系数分析两个独立变量的相关关系， 如NDVI与气温、 NDVI与降水的相关关系， 计算公式为^［31］

式中： r_xy 为相关系数， 值介于-1 ~ 1之间， 绝对值越接近于1， 相关性越强， 反之， 相关性越弱； x¯和y¯分别为变量x和y的平均值。计算得到相关系数后， 查找相关系数临界值表， α=0.05时， 临界值为0.3494， |r_xy |>0.3494， 相关关系显著， 反之则不显著。

式中： r12∙3为假定第3个变量不变， 变量1和变量2的偏相关系数； r₁₂、 r₁₃和r₂₃分别为变量1和变量2、 变量1和变量3、 变量2和变量3的相关系数。偏相关系数的检验采用t检验方法， 计算公式为^［31］

式中： n为样本数； m为自变量个数； 其他同前。t值计算出后， 查t分布表中α=0.05时的临界值为2.042， 若t>2.042， 表示偏相关关系显著， 反之， 则不显著。

青藏高原植被类型总体以草原、 草甸、 高山植被、 灌丛为主， 表现为NDVI均值总体较低， 1982 - 2015年多年平均值仅为0.3568， 研究时段内年均NDVI呈不显著的增加过程［图2（a）］， 2003年年均NDVI最低， 为0.3466， 2012年最高， 为0.3627， 二者相差较小， 仅为0.0161。据青藏高原NDVI距平值变化过程［图2（b）］， 可将NDVI变化过程划分为4个阶段， 第一阶段为1982 - 1987年， 该阶段NDVI多年平均值为0.3521， 相对于1982 - 2015年NDVI平均值而言， 该阶段NDVI以负距平为主， 即每年年均NDVI值均低于1982 - 2015年均值。第二阶段为1988 - 2000年， 该阶段NDVI多年平均值为0.3589， NDVI以正距平为主， 有8个年份距平值为正值， 5个年份距平值为负值， 这一阶段NDVI的波动较为复杂。第三阶段为2001 - 2008年， NDVI多年平均值为0.3537， 是4个阶段中NDVI平均值最低的时期。除2002年和2005年NDVI为正距平外， 其他年份均为负距平。第四阶段为2009 - 2015年， NDVI多年平均值为0.3603， 除2014年NDVI为负距平外， 其他年份均为正距平。该阶段NDVI的提升可能与2008年开始实施的青藏高原生态安全屏障工程建设有关。总体来看， 第一阶段和第三阶段NDVI低于1982 - 2015年青藏高原NDVI平均水平， 为植被水平相对较差的两个阶段。第二和第四阶段NDVI水平上升， 但第二阶段表现出NDVI升高的往复性， 第四阶段则表现出近期青藏高原NDVI呈现出增加过程， 且总体NDVI水平相对较高。

1982 - 2015年NDVI线性变化速率空间分布结果表明［图3（a）］， 青藏高原NDVI以线性增加为主， 占总面积的85.87%（不含地表水体、 裸地等无数据的栅格）。NDVI年平均变化速率集中分布在-0.002 ~ 0.004范围内， 共占总面积的97.80%， 其中-0.002 ~ 0类型占13.07%， 0 ~ 0.002类型占64.40%， 0.002 ~ 0.004类型占20.33%。从空间分布看， 1982 - 2015年NDVI增加速率较高的区域集中分布在青海省柴达木盆地的南部、 东部和东北部边缘地带， 西藏南部， 而NDVI呈减少速率， 且减少速率较高的区域大部分集中在四川省境内的青藏高原东部边缘、 小部分分布在西藏东部。

利用青藏高原1982 - 2015年NDVI线性速率及其F检验结果， 将NDVI变化趋势划分为不显著减少、 显著减少、 不显著增加、 显著增加。图3（b）表明， 青藏高原大面积区域NDVI以显著增加为主， 分布在高原北部、 西部和南部， 如青海省西部和北部、 新疆南部、 西藏西部和南部， 占总面积的63.26%。通过青藏高原省会城市、 地级城市和县区城市分部情况可知， 显著增加分布区域， 各级城市分布数量较少， 可以认为该区域人类活动相对较弱， 植被变化主要受到气候影响， 即水热条件的改善有利于植被的生长和发育， NDVI数值提升， 反之则限制植被生长和发育， NDVI数值下降。NDVI显著减少区域集中在各级城市分布密集的高原东南部， 如西藏东部、 青海省南部、 四川省西部及三省区交界区域。这些区域显著减少与不显著减少， 以及不显著增加交叉分布， 其中不显著增加面积比例占22.61%、 不显著减少占10.68%， 显著减少区域仅占3.45%。随着全球气候变暖， 青藏高原平均气温升高， 为植被生长和发育提供了有利条件， 促使NDVI在1982 - 2015年出现大面积的显著增加过程。

青藏高原植被类型以草原、 草甸、 高山植被、 灌丛、 荒漠、 针叶林和阔叶林为主， 共占总面积的91.45%。利用以上植被类型空间分布数据， 提取1982 - 2015年与植被类型相对应区域的NDVI， 探讨其时间变化过程和空间变化显著性结果， 有利于进一步认识青藏高原植被动态变化情况， 结果见图4。

选择的青藏高原7种植被类型中， 1982 - 2015年阔叶林植被NDVI平均值最高， 达到0.8628， 其次为针叶林植被， 为0.7197， 再次分别为灌丛、 草甸、 高山植被， 平均值分别为0.6335、 0.4899、 0.2524， 草原和荒漠植被类型NDVI平均值最低， 分别仅为0.1972和0.1032。从各植被类型NDVI时间变化过程看， 草原、 草甸、 高山植被和荒漠植被类型的NDVI呈增加过程， 其中草原、 高山植被、 荒漠植被类型NDVI呈显著增加过程， NDVI年平均线性增加速率分别达到0.0003、 0.0003、 0.0002［图4（a）、 4（b）、 4（c）、 4（e）］。随着全球平均气温升高， 高原平均气温也逐步升高， 高原植被生长所需的热量条件逐步改善， 使得植被总体水平提高。灌丛、 针叶林、 阔叶林植被类型NDVI呈减少过程， 但均不显著［图4（d）、 4（f）、 4（g）］。F检验显著性统计结果也表明［图4（h）］， 荒漠、 草原和高山植被的NDVI显著增加面积分别达到85.36%、 76.63%和65.43%， 草甸、 阔叶林、 针叶林、 灌丛植被类型NDVI显著增加的面积比例也较高， 但远低于高山植被、 草原和荒漠植被类型， 分别为51.93%、 50.69%、 42.80%和39.78%。相反， 灌丛是显著减少面积最高的， 其次为针叶林、 草甸和阔叶林植被。

结合各植被类型多年平均值可知， 1982 - 2015年青藏高原NDVI变化过程中， 平均值较高的植被类型， 包括阔叶林、 针叶林、 灌丛、 草甸， 表现出NDVI下降或小幅上升（草甸）的过程， 而平均值较低的植被类型， 如高山植被、 草原、 荒漠， 其NDVI表现出显著上升过程。反映出研究时段内， 青藏高原NDVI水平较低的植被类型增幅或增量高于NDVI水平较高的植被类型。可能反映出气候变化对青藏高原不同植被生态系统产生了差异化的影响， 如相对脆弱的高山植被、 草原、 荒漠植被生态系统得到改善， 而阔叶林、 针叶林、 灌丛植被生态系统出现了变差的迹象， 需要引起重视。

气温和降水是表征气候的重要指标， 探讨青藏高原NDVI对气温和降水变化的响应， 可为认识气候变化对青藏高原陆地植被生态系统影响提供帮助^［32］。以相关系数及偏相关系数t检验临界值为阈值， 将青藏高原NDVI与气温、 降水的相关系数和偏相关系数t检验结果划分为4种类型， 显著负相关、 负相关、 正相关和显著正相关。具体结果见图5。

气温是热量的表征指标， 随着全球气候变暖， 青藏高原气温也呈上升过程， 为植被生长提供了更多热量资源。青藏高原NDVI与气温相关系数空间分布呈纬度地带性特征， 正相关系数多集中在高原北部［图5（a）］； NDVI与降水相关系数空间分布呈东西向的干湿地带性分布， 南北向集中在区域中东部和中西部两个区域［图5（b）］。青藏高原1982 - 2015年NDVI与气温的相关系数和偏相关系数t检验结果均以正相关为主， 且空间分布较为一致， 呈南北向带状分布， 显著正相关均分布在青藏高原北部的青海省境内及青海省、 新疆、 西藏交界区域， 西藏南部地区也有集中分布， 分别占总面积的15.80%和15.15%； 显著负相关自西向东集中分布在青藏高原中南部地区， 主要位于西藏和四川省境内， 青海省东南与四川省交界处也有小面积分布， 分别占总面积的6.46%和7.59%［图5（c）、 5（e）］。青藏高原热量资源较少， 气温升高， 热量资源增加， 有助于植被生态系统的稳定和健康发展， 处于有利的一面。但NDVI与气温的负相关关系反映出， 随着气温的升高， 呈负相关关系地区的NDVI水平未实现增加， 反而出现了减少过程， 这一过程的驱动因素还不得而知， 需要另行研究。

青藏高原NDVI与降水的相关关系和偏相关关系均较NDVI与气温数值低， 但NDVI与降水的相关系数和偏相关系数t检验结果仍以正相关为主， 且空间分布较为一致， 呈东西向带状分布［图5（d）、 5（f）］。其中， 显著正相关主要分布在青海省北部、 西藏中南部和东部， 以及西藏和新疆的交界区域， 分别占总面积的9.12%和9.09%； 显著负相关的分布区域相对较小， 青海省东南部、 西藏西部和中南部均有小范围分布， 分别占总面积的1.80%和1.92%。

相关系数和偏相关系数分别反映2个和3个独立变量之间的相关关系， 通过青藏高原NDVI与气温、 降水的显著相关和显著偏相关面积比例可知， 在气温或降水交叉作用影响下， NDVI与降水、 与气温的显著负相关、 显著偏负相关面积比例均表现为增加， 而显著正相关、 显著偏正相关面积比例均表现为减少， 反映出气温、 降水分别在NDVI与降水、 NDVI与气温的相关关系中呈现出负面交叉作用， 导致显著负相关的比例升高， 而正相关的比例下降。

NDVI在反映全球气候变化和表征地表植被生态及其生产力状况方面具有重要作用。已有对青藏高原NDVI研究结果与本研究结果基本一致， 表明受全球气候变暖影响， 20世纪80年代以来青藏高原NDVI及模拟的净初级生产力等均呈增加过程^{［19， 33］}， 但由于数据时序长度不同， 导致NDVI变化的空间分布具有一定的差异， 如本文使用的GIMMS NDVI数据具有时间序列长的特点（1982年以来）， 但分辨率较低， 而MODIS NDVI数据分辨率较高， 但时间序列较短（2000年以来）^［6］， 二者结合使用可以满足长时间序列和刻画空间细节等方面的需求。

姚檀栋等^［34］提出了青藏高原现代西风与印度季风相互作用表现的3种模态， 即印度季风模态、 西风模态和过渡模态。印度季风影响下高原南部雅鲁藏布江流域湖泊面积普遍萎缩， 西风控制区的北部羌塘高原湖泊普遍强烈扩张^［35］。本研究表明1982 - 2015年青藏高原植被NDVI显著增加的区域主要分布在受西风控制的高原北部和西部， 而受印度季风控制的高原东南部植被呈减少态势， 印证了青藏高原印度季风模态和西风模态的在近30多年的表现特征。全球气候变暖在青藏高原不同区域NDVI变化方面产生了明显的差异， 即水热条件丰富的东南部地区NDVI下降， 而水热条件相对较差的高原北部等区域出现了NDVI的升高迹象。

郑然等^［36］研究认为青藏高原气温升高和降水增加， 导致高原干湿界线向西北方向移动。1982 - 2015年青藏高原年平均气温、 年降水量变化速率空间分布结果反映出高原总体年平均气温升高、 年降水量增加， 仅有高原南部和西部小范围区域表现为下降特征（图6）， 即青藏高原大部分区域表现为暖湿化趋势， 为植被生长提供了有利条件， 原本植被生态系统脆弱、 NDVI值较低的青藏高原高山植被、 草原、 荒漠植被NDVI呈显著增加趋势^{［6-7，28，30］}， 但青藏高原不同植被类型NDVI的增减过程是否与气候变化的周期性波动有关， 还需要通过设定NDVI周期性波动阈值， 深入分析20世纪80年代以来青藏高原的NDVI变化趋势是否在气候变化引起的植被波动范围内， 如果仅是NDVI随气候变化的周期性波动， 则无需过度关注， 但若超过了周期性波动范围， 则需引起重视。

（1）青藏高原1982 - 2015年NDVI总体上呈不显著的增加过程， 可划分为1982 - 1987年和2001 - 2008年NDVI以负距平为主， 1988 - 2000年和2009 - 2015年NDVI以正距平为主的变化过程。NDVI空间变化以显著增加为主， 分布在青藏高原北部、 西部和南部， 如青海省西部和北部、 新疆南部、 西藏西部和南部， 占总面积的63.26%。显著减少集中分布在高原东南部， 如西藏东部、 青海省南部、 四川省西部及三省交界区域， 仅占总面积的3.45%。

（2）青藏高原主要植被类型NDVI平均值表现为： 阔叶林>针叶林>灌丛>草甸>高山植被>草原>荒漠， 其中草原、 高山植被和荒漠植被的NDVI呈显著线性增加过程， 草甸植被呈不显著增加过程， 而灌丛、 针叶林和阔叶林植被的NDVI呈不显著的减少过程。这一过程可能反映了气候变化对青藏高原植被生态系统影响的空间差异性， 即青藏高原北部水热条件缺乏区域NDVI出现显著增加趋势， 而高原东南部水热条件充足地区呈现出显著减少趋势。

（3）青藏高原NDVI与气温相关系数空间上呈南北向分布， 具有纬度地带性特征， 显著正相关分布在高原中北部， 占15.80%， 显著负相关分布在高原中南部， 占6.46%； NDVI与降水的相关系数呈东西向分布， 具有干湿度地带性特征， 显著正相关分布在高原中部， 占9.12%， 显著负相关分布在高原东西两侧， 占1.80%。

①②