2000-2018年间青藏高原植被覆盖变化及其与气候因素的关系分析

韩炳宏， 周秉荣， 颜玉倩， 石明明， 苏淑兰， 赵恒和，牛得草， 傅 华

(1. 中国气象局成都高原气象研究所， 四川 成都， 610072； 2. 青海省防灾减灾重点实验室， 青海 西宁， 810001；3. 青海省气象科学研究所， 青海 西宁， 810001； 4. 青海省海南州气象局， 青海 共和， 813099；5. 兰州大学草地农业生态系统国家重点实验室， 甘肃 兰州， 730020)

植被作为全球陆地生态系统的重要组分，对全球物质能量循环、碳平衡调控及维持气候稳定等方面起着重要作用[1]。光照、气温和降水等气候要素是植被生长发育不可或缺的环境因子，对植被健康稳定的生长状况具有显著影响[2-3]。归一化植被指数(Normalized Difference Vegetation Index,NDVI)作为评价植被生长状况的重要因子，不仅可以反映全球或区域尺度上植被生长状况和分配格局，也可以表观长时间序列上的变化趋势[4]。通常与植被生产力和叶面积指数等重要植被特征参数密切相关[5]。有些研究结果表明，1982—1999年北半球中高纬度和中国大多数地区的NDVI表现出增加趋势[6-7]。而有些研究结果显示，2000—2009年北半球多个地区植被NDVI呈减弱趋势[8]。Piao等[9]对欧亚大陆温带和寒带生长季植被NDVI变化趋势进行了研究，认为1982—1997年显著增加，而1998—2006年呈下降趋势，尤其在春季和夏季表现最为明显；Mohammat等[10]认为，亚洲内陆地区生长季植被由于受春季变冷和夏季干旱的影响，其NDVI增加趋势于1990年停止；杜加强等[11]研究表明，1998—2012年新疆地区夏季植被NDVI由之前的极显著增加转变为显著减少。

青藏高原是中国面积最大、世界上海拔最高的高原，堪称地球“第三极”[12]，拥有独特的生态系统类型，被认为是亚洲乃至北半球气候变化的“感应器”[13]，不仅关系到高原周边地区数亿居民的供水安全，同时也对我国乃至东亚地区生态系统稳定起着屏障作用[12,14]。随着全球气候变暖和人类活动加剧，青藏高原正面临生态环境压力增加与生态安全屏障功能改变等风险[12]。由于受气候变化和人类活动的双重影响，青藏高原高寒植被变化受到越来越广泛的关注。近30a来，尤其是21世纪以来，青藏高原植被覆盖度呈不明显的绿度增加或变绿趋势[15-16]。赵紫薇[17]利用1982—2012年全球监测模型与制图研究归一化植被指数(Glaobal Inventory Modelling and Mapping Studies Normalized Difference Vegetation Index,GIMMS NDVI)和2001—2013年中尺度分辨率成像光谱归一化植被指数(Moderate Resolution Imaging Spectroradiometer Normalized Difference Vegetation Index,MODIS NDVI)对近30 a青藏高原植被动态变化时序进行了对比分析，发现青藏高原地区植被呈整体改善趋势；孟梦等[18]利用1982—2013年GIMMS NDVI数据研究了青藏高原NDVI变化趋势及其对气候变化的响应，指出近年来青藏高原植被长势逐渐变良好，覆盖度呈增加态势。以上研究虽得到了较为相近的结果，但由于受遥感数据时序的限制，并未涉及近5a的草地植被变化状况，加之研究方法和不同时段采用的遥感影像产品不同，获取的结果仅能代表过去某一时段草地植被生长状况的好坏，同时也缺少对不同时段草地植被与同期气象因子的关系分析。为此，本文基于遥感监测分析，利用2000—2018年MODIS NDVI逐旬数据对青藏高原植被覆盖时空变化特征及其驱动因素进行系统研究，以期为青藏高原植被在区域乃至全球气候变化过程中的调控机理提供科学依据。

1 材料与方法

1.1 研究区概况

青藏高原南起喜马拉雅山脉南缘，北至昆仑山、阿尔金山和祁连山北缘，西部为帕米尔高原和喀喇昆仑山脉，东及东北部与秦岭山脉西段和黄土高原相接，地处26°00′～39°47′ N，73°19′～104°47′ E，东西长约2.8×103km，南北宽约3.0×102～1.5×103km，总面积约2.5×106km2，地形上可分为藏北高原、藏南谷地、柴达木盆地、祁连山地、青海高原和川藏高山峡谷区6个部分[12-13]，主要包括中国西藏全部和青海、新疆、甘肃、四川、云南的部分地区(图1)。其自然历史发育极其年轻，受多种因素共同影响，东南部属于暖湿性气候，西北部属干冷性气候，具有太阳辐射强、气温低、气温日较差大、年变化小、干湿季节分明等特点；年降水量为486 mm，年均温在—5.6～8.6℃之间，年日照时数为2 300～3 600 h，太阳总辐射为140～180 C·m-2。高原上冻土广布，植被多为草地植被[19]。

图1 研究区地势及高程Fig.1 Terrin and elevation of study area

1.2 数据来源及处理

NDVI是通过红外与近红外波段的组合实现对植被信息状态的表达，通常用来表征植被覆盖程度、生长状况、叶面积指数、生物量以及光合有效辐射等植被参数[19]，已被广泛应用到全球及区域尺度上植被变化与气候互作研究中。

文中NDVI数据来自美国LP DAAC数据中心(http://lp-daac.usgs.gov/main.asp)MODIS仪器提供的2000年1月—2018年12月逐旬NDVI数据，分辨率为1km，数据格式为HDF-EOS，投影方式为Sinusoidal。数据的合成、再投影处理通过MRT和投影软件ENVI实现。根据青藏高原范围将下载的MODIS数据进行拼接，使其投影为标准的经纬度网格坐标，将软件输出的MODIS资料保存为GEOTIF格式；在ENVI里完成青藏高原矢量边界裁剪，获得研究区的NDVI数据。为更好地反映研究区不同时间尺度植被覆盖状况，将数据另存为ASCII格式，以便用于植被指数与气象要素的相关分析。

1.3 研究方法

1.3.1最大值合成法 最大值合成法(Maximum Synthesis Method,MVC)是目前国际上通用的最大合成法，采用该方法取一个月每旬的数据最大值为月值，进一步消除云、大气和太阳高度角的干扰[20-21]。此法假设一个月每旬中NDVI值最大的那一天天气是晴朗的，不受云层的影响，就取这个最大值作为月NDVI值，计算公式如下：

NDVIi=Max(NDVIij)

(1)

式中NDVIi为第i月的NDVI值，NDVIij为第i月第j旬的NDVI值，可以认为NDVIi是一月内植被最丰盛时期NDVI值[21]。利用(1)方法依次获得年和季节NDVI值，然后采用公式NDVI=0.004DN-0.1将其转化为真实NDVI值。NDVI的取值范围为—1.0～1.0，一般认为NDVI达到0.1以上表示有植被覆盖，值越大植被覆盖度越大；0.1以下则表示地表无植被覆盖，如裸土、沙漠、戈壁、水体、冰雪和云[22-23]。

1.3.2均值法 均值处理是将某时间间隔内的NDVI数据求均值，以消除或降低由时间段端点年份气候异常对植被生长状况的影响；每年的NDVI值由各月最大NDVIi求均值获得[24]，计算公式如下：

(2)

式中MNDVI,i表示第i年的NDVI，n表示月份，NDVIi是第i月的最大NDVI值。

1.3.3差值法 差值法是用于2个时期NDVI值的变化，即后一时期所有像元NDVI值与前一时期所有像元的差。宋怡等[25]利用该方法计算了年份间NDVI值的差；徐慧等[24]也同样采用该方法计算了年份间NDVI值的差。不同时期NDVI差值的计算公式如下：

ΔNDVI=MNDVI,i-MNDVI,j

(3)

式中NDVI表示NDVI差值，MNDVI,i表示后一年NDVI值，MNDVI,j表示前一年NDVI值。

1.3.4NDVI变化趋势(Slope)分析 一元线性回归分析可以模拟每个栅格的变化趋势，Stow等[20]就用该方法来模拟植被的绿度变化率(GreennessRateofChange,GRC)，GRC被定义为某时间段内的季节合成归一化植被指数年际变化的最小次方线性回归方程的斜率。此处同样用该方法来模拟多年NDVI的变化趋势，计算公式为：

(4)

式中变量i为1～19序号，MNDVI,i表示第i年的NDVI值。变化趋势图则反映了近19a青藏高原地区的NDVI的变化趋势。某格点的趋势线是这个格点19a的NDVI值用一元线性回归模拟出来的一个总的变化趋势，即θSlope为这条趋势线的斜率。这个趋势线并不是简单的最后一年与第一年的连线。若θSlope>0，则说明NDVI在19a间的变化趋势是增加的；若θSlope<0则是减少。另外，利用ArcGIS空间分析工具Reclass命令，将θSlope值进行重分类，依次划分为严重退化[—0.030，—0.010]、中度退化[—0.010，—0.005]、轻微退化[—0.005，—0.003]、保持不变[—0.003，0.003]、轻微改善[0.003，0.005]、中度改善[0.005，0.010]和明显改善[0.010，0.030]7个等级[19]；再通过Zonal Satistics命令，得出7个等级所对应的像元数、面积、均值、标准差和面积百分比。

2 结果与分析

2.1 2000—2018年年际NDVI变化

通过2000—2018年的年平均NDVI进行逐年差值分析，从变化幅度上看，逐年NDVI的差值结果显示2001，2006和2012年分别在2000，2005和2011年的基础上显著增加，增加幅度分别为0.126，0.106和0.11；2005，2007，2008和2009年青藏高原年均NDVI较2004，2006，2007和2008年分别减少了0.02，0.032，0.052和0.046。另外，除2001和2006年外，年平均NDVI在2000—2009年呈减少趋势；2010—2018年呈不同程度的增加态势。整体来看，青藏高原地区近19a年平均NDVI呈逐渐增加趋势(图2)。

图2 2000—2018年年平均NDVI差值Fig.2 Difference of annual average NDVI from 2000 to 2018

2.2 2000—2018年NDVI季节变化

本文采用气候学上公认的季节划分方法，3—5月为春季、6—8月为夏季、9—11为秋季、12—翌年2月为冬季，分别将其NDVI平均获得季节NDVI，再将19个年份同一季度NDVI求均值，从而获得4张季节NDVI影像，以分析青藏高原2000—2018年NDVI平均生长状况在季节上的空间分布及差异(图3)。

整体来看，2000—2018年青藏高原NDVI自西北向东南依次增加，且四季NDVI空间分布差异较大(图3)。从地理位置来看，低植被覆盖区(NDVI<0.2)主要分布于西藏大部、新疆和甘肃局部以及青海西北部地区；中植被覆盖区(0.2≤NDVI≤0.3)主要位于青海与甘肃、西藏、四川和云南的交界区；高植被覆盖区(NDVI>0.3)主要分布在四川和云南大部、青海和甘肃以及西藏东南局部地区。从草地类型分布来看，低植被覆盖区主要以荒漠、戈壁和裸地为主，中植被覆盖区以高寒草原和温性草原为主，高植被覆盖区以高寒草甸为主。

图3 2000～2018年NDVI季节变化Fig.3 NDVI seasonal changes from 2000 to 2018

2.3 趋势分析

从NDVI趋势分析结果(图4)可以看出，在四川和云南局部、青海中东部及南部、西藏东南部及北部部分地区，NDVI的变化趋势呈明显的上升趋势，说明近19a中这些区域的植被情况得到比较好的改善。而在新疆和甘肃局地、青海东南部及东北部、西藏东南部及其西部的边远地区，NDVI呈下降趋势。尤其在和田南部局地、林芝大部、祁连山部分地区及果洛藏族自治州大部，降低的趋势比较明显，即在这些地区植被退化情况较为严重。另外，由表1可知19a年来，草地植被退化总面积为60.03 km2，恢复改善总面积达69.13 km2。故整个青藏高原地区植被生长状况得到改善的区域面积大于植被退化的面积。其中，植被改善区域面积占整个青藏高原总面积的27.35%，退化区域占23.75%。因此，青藏高原植被虽局部恶化，但整体仍处于恢复状态。这也进一步验证了前面的分析结果。

图4 2000—2018年MNDVI变化趋势Fig.4 Change trends of MNDVI from 2000 to 2018

表1 2000—2018年MNDVI变化趋势结果统计Table 1 Statistics of MNDVI change trend from 2000 to 2018

θSlope植被变化趋势Vegetation variation trend像元数Pixelelements面积Areas/ km2平均值Average标准差Standard deviation百分比Percentage %—0.030～—0.010严重退化Serious degradation48 8073.89×1040.170.141.54—0.010～—0.005中度退化Moderate degradation145 01711.00×1040.150.154.35—0.005～—0.003轻微退化Slight degradation566 20945.14×1040.240.1617.86—0.003～0.003保持不变Remain unchanged1 551 107123.65×1040.140.1248.910.003～0.005轻微改善Slight improvement616 03649.11×1040.180.1419.430.005～0.010中度改善Moderate improvement193 92016.03×1040.530.136.340.010～0.030明显改善Serious improvement50 1423.99×1040.670.061.58

2.4 2000—2018年NDVI与气候因子的关系

全球气候变暖有可能加速土壤有机质的分解而促进较为寒冷的地区植被净初级生产力的提高已基本形成共识[1,12-13]。在大部分相对干旱的区域，植被生长与降水呈正相关；在湿润地区，则为负相关。本研究发现19a中青藏高原地区NDVI与气温和降水的关系较好(图5a，b)。其中，NDVI与气温的相关性系数达0.909；NDVI与降水的相关性虽不及气温，但相关系数也能达到0.793。水在植物的营养物质输送、结构维持和各种生理生化过程中起着十分重要的作用，因此植被的生长与降水量的多少具有明显的正相关关系；植物叶片的净光合速率随温度的升高而增加，气温过高或过低都会影响植物的生理生化过程。从两者的线性关系可以看出，当气温每升高1℃，NDVI增加0.128；降水每增加100 mm，NDVI相应增加0.172。

图5 NDVI与气象因子的关系Fig.5 Relationship between NDVI and meteorologic factors

3 讨论

青藏高原作为亚洲生态安全屏障和全球气候变化敏感区，其植被生长状况对区域或全球尺度上气候变化的响应具有重要的调控作用。为此，本文利用2000—2018年MODIS NDVI数据从季节、年际水平对青藏高原不同地区、不同时间尺度上植被覆盖变化及其驱动因素进行了分析。结果表明，近19a来，青藏高原地区NDVI虽局部有所恶化，但整体仍呈稳定恢复趋势；与卓嘎等[19]的研究结论一致。然而，Shen等[26]利用增强型植被指数(Enhanced Vegetation Index,EVI)研究了2000—2012年生长季青藏高原植被覆盖变化，认为青藏高原植被覆盖变化呈不显著减少趋势，主要表现为植被减少速率大于植被增加速率；且不同植被类型对气候变化的敏感性不同。其研究结论与本文不一致，这主要是因为Shen等[26]的研究时段以及研究方法不同。另外，我们发现自21世纪以来，特别是2000—2009年期间，青藏高原植被覆盖呈现不显著的减少趋势，即植被退化面积和速率大于植被恢复，2010—2018年间植被覆盖不同程度的增加，这可能与近10年青藏高原地区气温和降水的明显增多有关[19]。宋怡等[25]基于1998—2004年SPOT_VGT NDVI数据对中国西北地区植被变化特征进行了系统研究，研究表明：西北地区植被覆盖在1998—2004年间出现退化趋势，尤其以2000—2002年间的变化幅度较大；青藏高原地区植被虽局部区域有改善，但总的改善幅度小于退化幅度。由此可见，青藏高原植被覆盖在2010年之后出现了较大幅度的改善，且逐年恢复的速率和面积大于退化速率和面积。

植物对气温和降水的响应机制主要表现在其对植物生理生化过程中酶活性的影响方面，水是植物代谢途径中的主要生理生化过程(糖酵解、三羧酸循环、磷酸戊糖、乙醛酸循环和乙醇酸氧化等)的参与者，而气温则对植物代谢途径中的一系列酶(氧化酶、还原酶和脱氢酶等)的活性产生重要影响。以往研究表明，青藏高原气候条件与NDVI的相关存在明显的区域性差异。青藏高原东北部气温和降水与NDVI具有较好的正相关关系，西南部呈负相关关系[19]。孟梦等[18]研究发现青藏高原植被覆盖与空气中的水汽含量显著相关，对最低气温的响应最为敏感。我们发现青藏高原地区NDVI与气温和降水的相关性较好。其中，NDVI与气温间的相关性大于其与降水间的相关性，这与以往研究[18-19,26-28]结论一致；这也说明在全球气候变暖背景下，青藏高原地区气候在暖干化向暖湿化演变的过程中，青藏高原地区植被覆盖对气温响应的敏感性高于降水。

由于青藏高原下垫面类型分布多样，地貌特征复杂，降水格局不均等特点，青藏高原地区植被覆盖和类型表现出明显的区域性差异。另外，青藏高原地区气象站点主要集中于高原中东部，高原西北部地区站点相对较少，甚至个别区域无观测站点，故利用本文所得结论对整个青藏高原地区植被覆盖进行描述仍存在一定差异。因此，今后需采用气候格点数据进一步探讨青藏高原不同地区植被覆盖差异及其对气候变化的响应及互馈机制。

4 结论

近19a来，青藏高原地区植被NDVI主要经历了2个显著阶段，即前10a呈减少趋势，后9a呈不同程度的增加态势。NDVI整体呈逐渐增加趋势。从地理位置来看，低植被覆盖区主要分布于西藏大部、新疆南部和甘肃局部以及青海西北部地区；中植被覆盖区主要位于青海与甘肃、西藏、四川和云南的交汇区域；高植被覆盖区主要分布在四川和云南西北部大部分地区、青海和甘南以及西藏东南局部地区。NDVI趋势分析结果显示，除高原局部地区植被有所退化，但大部地区植被生长状况明显改善，且草地植被改善的面积大于退化。因此，青藏高原植被虽局部恶化，但整体仍处于恢复状态。青藏高原地区NDVI与同期气温和降水间均有较好的正相关关系，但气温要好于降水。且气温每升高1℃，NDVI增加0.128；降水每增加100 mm，NDVI相应增加0.172。