青藏高原NDVI变化趋势及其对气候的响应

孟 梦, 牛 铮, 马 超, 田海峰, 裴 杰

(1.中国科学院 遥感与数字地球研究所, 遥感科学国家重点实验室, 北京 100101; 2.中国科学院大学, 北京 100049; 3.河南理工大学 测绘与国土信息工程学院, 河南 焦作 454000)

青藏高原总面积约250万km2，约占中国陆地总面积的1/4，平均海拔在4 000 m以上，是中国面积最大、世界上海拔最高的高原，被称为世界的“第三极”[1]。其特殊的中低纬度地理位置、地貌条件以及巨大的海拔高度，形成了湿润、干旱、热带和寒带等多种气候类型，几乎囊括了东南亚所有类型的生态系统[2]。青藏高原具有辽阔的生态脆弱区，生态系统敏感而脆弱，且受人类影响较小，这也成为其受广大生态学研究者钟爱的一个重要原因[3]。

近年来，在全球变化和人为因素的共同作用下，青藏高原生态环境产生了一系列负面变化，具体表现为湿地草地面积减小、冻土退化、土地沙漠化严重、水域减少、多种动植物濒临灭绝以及自然灾害增多等现象[4]。青藏高原作为中国乃至亚洲地区的生态安全屏障，其“一举一动”都备受社会和学术界关注，因此科学有效的评判青藏高原生态环境变化以及探讨其变化规律成为许多学者研究的重要课题[5]。植被作为高原生态系统中重要的组成部分，植被覆盖变化在一定程度上能够反映生态环境的变化，因此植被覆盖变化研究对青藏高原生态环境变化规律的探索具有深刻意义。

归一化植被指数(Normalized Difference Vegetation Index，NDVI)对植被的生长量和长势十分敏感，是反映地表植被生长状况的重要指标，因而在定性和定量评价植被覆盖及其生长活力时被作为重要的信息源[6-7]。其定义为近红外波段与可见光红外波段波谱数值之差与两波段波谱数值之和的比值，计算公式为NDVI=(NIR-IR)/(NIR+IR)，NDVI能较好的反映植被生长状况，且可以消除地形、部分太阳角以及云等因素的影响[8-9]，其次，NDVI对干旱和半干旱区地区的低值被覆盖信息十分敏感[10]，因此在低值被覆盖区用NDVI反映植被信息得到了广泛应用。

近年来，很多学者利用NDVI研究青藏高原植被的动态变化及植被与温度、降水等气候因子的关系。李艳芳等[11]基于SPOT-Vegetation植被指数数据，得出了青藏高原NDVI总体呈上升趋势，季节变化趋势不明显的结论。向波等[12]基于全球归一化植被指数(Glob-NDVI)发现近年来青藏高原植被向好的方面发展的现象。杨元合等[13]利用NOAA-AVHRR NDVI数据，发现青藏高原植被指数在生长季阶段显著增加，春季对NDVI增长率贡献最大，且不同植被类型在不同的生长季NDVI增长率也有所区别。王青霞等[14]基于1982—2006年GIMMS NDVI数据，研究表明青藏高原植被与温度降水等气候因子具有相关性，不同等级的植被覆盖度对温度降水的响应程度是不同的，植被生长较好的区域对气候响应更为明显。杨保等[15]利用冰芯中的花粉以及氧同位素记录，恢复了青藏高原过去500多年的植被演变与气候变化趋势，研究表明当高原为暖干气候时，植被生长范围扩大，高原气候变冷变湿则植被生长范围缩小。陆晴等[16]基于GIMMS NDVI和地面气象站台观测数据，发现青藏高原植被与温度降水的相关关系具有空间差异性，在温度较高、热量条件较好的高原东北部地区，降水为植被NDVI变化的主导因素；在降水充足的高原东中部地区，温度为植被NDVI变化的主导因子。

尽管前人做了大量的研究，但多以行政边界划分研究区域，并且研究区中包括了西藏自治区东南部、四川省西南部、云南省西北部和青海省中部植被十分茂盛的地区，且多数研究所用遥感数据和气象数据时序短，因此不能突出地反映青藏高原高原植被的时空变化。鉴于此，本文以地理界线划分研究区，突破人为划定的行政边界，选取青藏高原海拔4 000 m以上的生态极端脆弱敏感区作为本次的研究区域，此外，选用长时间序列1982—2012年GIMMS NDVI数据、1982—2012年温度降水数据为数据源，结合多种统计方法，分析青藏高原植被NDVI的变化趋势及NDVI对气候的响应，旨在揭示全球变化背景下青藏高原植被变化趋势以及气候对高原植被生长的影响机制。

1　研究区与数据

1.1　研究区概况

青藏高原雄踞亚洲东南，西起帕米尔高原，东到横断山脉，北以昆仑山、阿尔金山和祁连山为邻，南抵喜马拉雅山脉，西南边缘部分分属印度、巴基斯坦、尼泊尔、锡金、不丹及缅甸等国，其地理位置大致为27°17′—40°19′N，68°52′—102°43′E(图1)。青藏高原气温随纬度和海拔的升高而降低，气温日差大，年平均气温和降水量由东南向西北递减。该区域的气候变化不仅对北半球有巨大的影响，甚至对全球气候变化也具有显著的超前性和敏感性。

1.2　数据及数据预处理

1.2.1数 据文中所用主要数据：GIMMS NDVI数据，来源于美国国家航空航天局(https:∥ecocast.arc.nasa.gov/data/pub/)，时间分辨率为15 d，空间分辨率为8 km×8 km，时间跨度为1982—2012年，31 a的数据共包括了744期半月合成图像。GIMMS NDVI数据经过辐射定标、大气校正、除云、除坏线等处理，是目前研究干旱、半干旱区植被变化质量比较可靠的NDVI数据集[17]。气象数据来源于中国气象科学数据共享服务网(http:∥cdc.nmic.cn)提供的气候资料年值数据集。高程数据来源于GTOP30项目所提供的DEM数据，分辨率为90 m。

1.2.2数据预处理首先利用青藏高原矢量数据对GIMMS NDVI数据批量裁剪，其次采用最大值合成法(Maximum Value Composite，MVC)获取青藏高原月度、年度NDVI数据，最后对得到的月度、年度NDVI值进行计算，获取研究区月度、年度NDVI均值并统计入Excel表格。

位于青藏高原地区的气象站点数目极少，本文选取该数据集青藏高原地区102个地面气象站台的气象数据(1982—2012年)进行克里金插值，获取分辨率为8 km×8 km的年平均气温和年平均降水量图像，从而获取青藏高原地区年气温、年降水量信息。

使用ArcGIS软件中的双线性插值法将DEM数据重采样为和NDVI相同的分辨率，即8 km×8 km，其次在Global Mapper中生成海拔为4 000 m的等高线，确定本次研究区范围，即青藏高原海拔高于4 000 m的区域。

图1研究区域及气象站点分布

2　研究方法

(1) 最大值合成法。为了研究青藏高原地区植被在空间上的整体变化特征，采用最大值合成法获取研究区月NDVI值和年NDVI值，最大值合成法可以有效减少大气中气溶胶、云阴影、太阳高度角等因素的影响[18]。计算公式如下：

(1)

式中：MNDVIi为第i月NDVI最大值；i=1,2，…，12为月序号；NDVI1为第i月上半月NDVI值；NDVI2为第i月下半月NDVI值。以求取的月NDVI值为基础求取年NDVI。

(2) 植被变化分析。结合植被绿度变化率(Greenness Rate of Change,GRC)[19]分析青藏高原植被变化趋势。GRC定义为在一定时间范围内采用最小二乘法拟合年度NDVI均值的斜率，它能够反映出每个栅格的变化趋势[20]。计算公式如下：

(2)

式中：变量i为年份序号，取值范围为1～31；n为研究的时间序列长度，n=31；YMNDVIi为第i年的年度NDVI平均值。GRC>0说明NDVI在31 a间的变化趋势是增加的，反之则是减少的[21]。

(3) 植被覆盖度分析。计算公式如下[22]：

(3)

(4) NDVI变化率分析。由下式计算：

(4)

式中：直线斜率表示31 a(1982—2012)间NDVI年平均值与年份之间求直线回归所得直线的斜率，均值为31 a的年度NDVI平均值，该变化率在数值上相当于初期值除以研究期间年平均NDVI的末期值与初期值之差[22]。

(5) 相关性分析。此外，考虑到植被生长与地形有密切的关系，借助GTOP 30DEM数据，分析了青藏高原归一化植被指数随地形的变化情况。

3　结果与分析

3.1　青藏高原NDVI时间序列变化特征

图2　青藏高原月度NDVI变化趋势

3.1.2年度NDVI变化将利用最大值合成法获取的的年NDVI进行计算得到年NDVI均值，将其与年份线性回归，得到青藏高原31 a NDVI均值变化趋势(图3)。可知31 a来青藏高原总体年NDVI均值呈缓慢增加趋势(slope=0.000 2)，利用公式(4)可计算得到青藏高原31 a年NDVI变化率为0.029，表明该研究区植被向好的方向发展，与王敏等[25]得出的结论一致。

图3　青藏高原年度NDVI变化趋势

植被覆盖度是单位面积内植被(包括叶、茎、枝)垂直投影面积所占的百分比[26]，它是陆面过程模式中最重要的参数[27]。植被覆盖度不仅可以表示植被进行光合作用面积的大小，它在一定程度上对植被的茂密程度、生长态势以及区域生态系统环境变化也有所指示[28]。由公式(3)计算得到青藏高原31 a植被覆盖度(图4)，将其与年份线性回归，可知青藏高原31 a间植被覆盖度呈0.000 3/a 的速率缓慢增长，其中在1988年、1994年、2000年、2010年有大幅度的跃升。认为全球变化对青藏高原起到了正面影响，植被长势逐渐向好的方向发展，植被覆盖度相应增加[29]。

3.2　青藏高原NDVI对气候因子的响应

陆气与气候间相互作用的变化，通过降水机制改变和气温变化等对植被的生长和分布产生重要影响。图5A、图5B分别展示了青藏高原1982—2012年31 a年平均气温以及年平均降水量的变化趋势。

图4　青藏高原31 a植被覆盖度变化趋势

图5A展示了青藏高原1982—2012年共31 a的年平均降水量变化趋势，年最低降水量为255.2 mm(1994年)，最高降水量为343.7 mm(2010年)，31 a平均降水量为309.9 mm。1982—1995年，年降水量呈下降趋势，1995—2012年，年降水量呈上升趋势，31 a青藏高原地区总体降水量呈逐渐增加趋势，变化率为0.223 mm/10 a。

由图5B可以看出，近31 a青藏高原海拔高于4 000 m的区域平均气温呈明显增加趋势，年最低气温为4.82℃(1983年)，最高气温为6.99℃(2009年)，增长速率为0.018℃/10 a。气温增长速率在1995年前后存在转折点，1982—1995年年平均气温增长速率为0.037℃/10 a，1995—2012年年平均气温增长速率为0.029℃/10 a，宋辞等[30]认为青藏高原气温发生突变，整体呈波动上升的趋势。

图5　青藏高原31 a降水、温度变化趋势

对青藏高原地区NDVI与降水、温度做相关性分析，见图6A,6B所示。研究区年均NDVI与降水、温度的皮尔逊相关系数分别为0.183，0.281，均呈现不显著的弱相关。相关系数表明该地区温度对植被的影响大于降水对植被的影响。此外，对降水、温度与NDVI做超前滞后分析，发现降水对研究区NDVI有三年的滞后效应，与刘雅勤等[31]的研究结果保持一致，当降水向后移动三年时，降水与NDVI的相关系数达到最大值0.387。温度对植被指数滞后性影响不明显。

由图6A，6B可以看出，青藏高原地区年NDVI均值在2000—2008年呈现明显下降趋势。对此，分别对2000—2008年的年NDVI均值、年降水量、年气温与年份做一元线性回归分析(图7)。由图7可以看出，2000—2008年年降水量呈下降趋势，年气温呈升高趋势，温度的升高及降水量的减少有可能打破了青藏高原地区降水量与蒸发量之间的平衡状态，使青藏高原干旱地区更加干旱。已有研究表明[32]，干旱对青藏高原高山草甸植被的生长具有抑制作用，加剧了植被的退化，因此对总体区域植被NDVI产生了重要影响。

此外，在2008—2012年年NDVI均值呈增长趋势，对于这一变化趋势，考虑到时间序列短、所用气象数据和NDVI数据集分辨率低、以及青藏高原自然环境复杂等因素，因此尚不能对此情况做出合理的解释。

图6　青藏高原地区NDVI、降水相关曲线NDVI(A)、温度相关曲线(B)

图7　青藏高原地区(2000-2008)年NDVI、降水量、气温变化趋势

3.3　青藏高原NDVI与海拔的相关性分析

高原植被在全球变化中发挥着重要的指示作用，因此备受关注。高原植被不仅对温度、降水等气候因子响应敏感，而且受海拔高度的影响也很大[33]。在研究区2000年的NDVI影像中均匀选取100个样本点，分别统计它们所对应的NDVI值和高程值，将NDVI值和高程值映射到以NDVI值为纵坐标，海拔高度为横坐标的的坐标系中(图8)。结果表明整体NDVI与海拔高度呈现相关程度不是很强的负相关，二者皮尔逊相关系数为-0.309，究其原因可能是植被类型过多，且植被指数受植被覆盖度、海拔高度、气候因子等综合影响，因此使海拔高度的影响变的不是太明显[34]。

在ArcGIS软件中用2000年的NDVI影像生成间距为0.1的等值线(图9)，并将等值线与青藏高原DEM影响叠加，可知在青藏高原中部、东南部、以及西北部海拔较低的区域NDVI等值线分布比较密集，印证了NDVI与海拔相关的事实。

图8　青藏高原地区NDVI与DEM相关曲线

图9　青藏高原地区NDVI等值线与DEM叠加

4　结 论

(1) 在全球变化影响下，1982—2012年各研究区NDVI值呈现以年为周期的变化特点，年NDVI值表明，植被的生长旺季处于每年的7月、8月、9月份，且年NDVI值总体呈增长趋势，年NDVI变化率为0.029，表明青藏高原地区植被生长状况发展良好。

(2) 青藏高原地区1982—2012年植被覆盖度呈0.000 3/10 a的速率增长，认为全球变化对青藏高原起到了正面影响，植被长势逐渐向好的方向发展，植被覆盖度相应增加。

(3) 青藏高原地区温度、降水都呈增长趋势，在1995年之前，年降水量呈下降趋势，1995年之后，年降水量呈增长趋势。温度在1995年之后增长趋势比1995年之前慢，总体增长趋势为前高后低。

(4) 研究区NDVI与温度降水的皮尔逊相关系数分别为0.281，0.183，均呈现不显著的弱相关，降水对归一化植被指数的影响有三年的滞后效应，温度对植被指数滞后性影响不明显。

(5) 植被生长受海拔高度影响明显，NDVI与海拔高度皮尔逊相关系数为-0.309，青藏高原中部、东南部、以及西北部海拔较低的区域植被长势相对较好，且NDVI值相对较高。

此外，青藏高原地区地面气象站台分布稀疏，即使我们采用了国际相对认可的克里金插值方法，但气象站台稀疏地区插值数据的不稳定性仍然是存在的，这种不稳定性对于后期的趋势分析会产生影响，因此，采用更多的气象观测数据来提高青藏高原地区气候插值的结果精度是十分必要的。

参考文献：

[1]郑度,赵东升.青藏高原的自然环境特征[J].科技导报,2017,35(6):13-22.

[2]莫申国,张百平,程维明,等.青藏高原的主要环境效应[J].地理科学进展,2004,23(2):88-96.

[3]梁四海,陈江,金晓媚,等.近21年青藏高原植被覆盖变化规律[J].地球科学进展,2007,22(1):33-40.

[4]姚檀栋,朱立平.青藏高原环境变化对全球变化的响应及其适应对策[J].地球科学进展,2006,21(5):459-464.

[5]孙鸿烈,郑度,姚檀栋,等.青藏高原国家生态安全屏障保护与建设[J].地理学报,2012,67(1):3-12.

[6]田庆久,闵祥军.植被指数研究进展[J].地球科学进展,1998,13(4):327-333.

[7]孙红雨,王长耀,牛铮,等.中国地表植被覆盖变化及其与气候因子关系：基于NOAA时间序列数据分析[J].遥感学报,1998,2(3):204-210.

[8]韩瑞梅,刘培,马超,等.鄂尔多斯植被的NDVI3g动态及气候响应[J].水土保持通报,2016,36(5):28-33.

[9]李春晖,杨志峰.黄河流域NDVI时空变化及其与降水/径流关系[J].地理研究,2004,23(6):753-759.

[10]刁鸣军,夏朝宗.1982—2013年准噶尔盆地植被长势变化分析[J].林业资源管理,2016(5):39-46.

[11]李艳芳,孙建.青藏高原NDVI时空变化特征研究(1982—2008)[J].云南农业大学学报,2015,30(5):790-798.

[12]向波,缪启龙,高庆先.青藏高原气候变化与植被指数的关系研究[J].高原山地气象研究,2001,21(1):29-36.

[13]杨元合,朴世龙.青藏高原草地植被覆盖变化及其与气候因子的关系[J].植物生态学报,2006,30(1):1-8.

[14]王青霞,吕世华,鲍艳,等.青藏高原不同时间尺度植被变化特征及其与气候因子的关系分析[J].高原气象,2014,33(2):301-312.

[15]杨保,唐领余,李春海,等.近550年来青藏高原中部植被演化与气候变化研究[J].科学通报,2010,55(21):2139-2147.

[16]陆晴,吴绍洪,赵东升.1982—2013年青藏高原高寒草地覆盖变化及与气候之间的关系[J].地理科学,2017,37(2):292-300.

[17]张戈丽,徐兴良,周才平,等.近30年来呼伦贝尔地区草地植被变化对气候变化的响应[J].地理学报,2011,66(1):47-58.

[18]郭鹏,徐丽萍.基于GIMMS—NDVI的新疆植被覆盖时空变化[J].水土保持研究,2014,21(4):97-100.

[19]陈志彪,朱鹤健.不同水土流失治理模式下的土壤理化特征[J].福建师大学报:自然科学版,2006,22(4):5-9.

[20]仝莉棉,曾彪,王鑫.山西省不同生态区NDVI时空变化及其影响因素[J].水土保持研究,2016,23(3):71-76.

[21]张戈丽,欧阳华,张宪洲,等.基于生态地理分区的青藏高原植被覆被变化及其对气候变化的响应[J].地理研究,2010,29(11):2004-2016.

[22]方精云,朴世龙,贺金生,等.近20年来中国植被活动在增强[J].中国科学,2003,33(6):554-565.

[23]刘双俞,张丽,王翠珍,等.基于MODIS数据的青藏高原植被物候变化趋势研究(2000年—2010年)[J].遥感信息,2014，29(6):25-30.

[24]李娜,谢虹.青藏高原NDVI对ENSO事件的响应[J].遥感技术与应用,2014,29(5):818-822.

[25]王敏,周才平,吴良,等.近10a青藏高原干湿状况及其与植被变化的关系研究[J].干旱区地理,2013,36(1):49-56.

[26]郭敏杰,张亭亭,张建军,等.1982—2006年黄土高原地区植被覆盖度对气候变化的响应[J].水土保持研究,2014,21(5):35-40.

[27]何建军,余晔,陈晋北,等.植被覆盖度对兰州地区气象场影响的模拟研究[J].高原气象,2012,31(6):1611-1621.

[28]穆少杰,李建龙,陈奕兆,等.2001—2010年内蒙古植被覆盖度时空变化特征[J].地理学报,2012,67(9):1255-1268.

[29]张宪洲,杨永平,朴世龙,等.青藏高原生态变化[J].科学通报,2015,60(32):3048-3056.

[30]宋辞,裴韬,周成虎.1960年以来青藏高原气温变化研究进展[J].地理科学进展,2012,31(11):1503-1509.

[31]刘雅勤,范广洲,周定文,等.青藏高原冬、春植被归一化指数变化特征及其与高原夏季降水的联系[J].气象学报,2007,65(6):959-967.

[32]王敏,周才平,吴良,等.近10a青藏高原干湿状况及其与植被变化的关系研究[J].干旱区地理,2013,36(1):49-56.

[33]魏伟,赵军,王旭峰.天祝高寒草原区NDVI, DEM与地表覆盖的空间关系[J].干旱区研究,2008,25(3):394-401.

[34]万红,李希灿,万剑华,等.青藏高原地区不同植被类型MODIS植被指数与海拔高度的关系研究[J].干旱区资源与环境,2013,27(11):39-44.