在全球气候变化、人类活动以及温室气体增加的背景下, 青藏高原作为气候变化的敏感区和生态脆弱带, 呈现出变暖趋势, 促使高原冰川退缩、湿地面积萎缩、积雪融化、冻土活动层厚度减薄, 草原冻融等问题愈发突出, 导致区域性的农牧区草地退化和沙漠化发生^[1-4]。高寒草原和草甸是当地牧民的主要牧草来源, 植被覆盖特征、生物量生产和分配、群落结构等的改变, 直接影响到农牧民的经济生活, 并且青藏高原生态环境变化对世界生态系统具有深刻的影响^[5-8]。因此, 对高寒植被的分布特征和变化规律及其对气候的响应进行研究, 能够加深对青藏高原气候和生态变化机理的认识, 具有十分重要的生态价值和现实意义。

以往针对青藏高原植被覆盖的研究主要包括植被覆盖的时空分布特征、影响植被生长的气候因子以及高寒草甸的物候信息等方面。高原植被时空分布特征研究主要采用不同来源卫星遥感监测的分析方法。朴世龙等^[9]对我国1982—1999年以来的植被覆盖动态变化进行了分析, 指出20世纪80年代和90年代青藏高原是植被覆盖度下降的主要分布区之一。对青藏高原局部地区植被变化而言, 卓嘎等^[10-11]分析了西藏地区植被的分布与变化趋势, 讨论了气候变化以及人类活动因素对不同区域植被变化特征的影响。边多、除多等^[12-14]认为阿里植被长势有所改善, 那曲草地退化出现先快后慢的现象, 拉萨植被与降水量的关系较气温密切。吕洋等^[15]基于趋势分析法和Hurst指数法, 指出2001—2012年雅鲁藏布江流域植被覆盖总体上呈现上升趋势。廖清飞等^[16]分析了青海省东部农业区植被覆盖度时空变化趋势及突变特征, 认为研究区降水量的变化以及退耕还林措施是植被覆盖度变化的重要影响因素。

在对高原整体植被覆盖特征研究中, Peng等^[17]基于Hurst指数指出青藏高原植被具有增长趋势, 而刘军会等^[18]利用生态模型揭示出植被变化趋势呈现整体升高局部退化的现象。不同时间尺度植被的时空变化表明^[19]青藏高原植被整体变化呈波动上升的趋势。李斌等^[20]利用半变异函数分析等方法探讨了青藏高原植被的时空变化特征。李辉霞等^[21]通过分离气候要素与人类活动对植被指数的贡献, 定量评估三江源地区生态保护与建设工程的实施效果。青藏高原降水量偏多能够促进植被的生长, 两者呈现较好的正相关^[22-23], 降水量与植被覆盖相关较好的区域主要位于高原东北部^[24]。气温升高有利于高山植物生长^[25], 而气温过高不利于高寒草地生长发育^{[24, 26]}。徐满厚等^[27]则对高寒草甸植被进行增温模拟试验, 认为适当的增温与降水均可极显著促进高寒草甸植被生长。针对青藏高原植被物候研究, 马晓芳等^[28-31]提取了青藏高原高寒草地的物候信息(植被返青期、枯黄期和生长季长度), 分析了不同高寒草地植被物候10年平均状况的变化趋势、空间分异特征及其对气候变化的响应规律。

综上所述, 以往对青藏高原植被覆盖研究多数采用美国大气海洋局(National Oceanic and Atmospheric Administration, NOAA)植被资料, 且侧重青藏高原不同区域植被覆盖分布特征或者高寒草甸物候信息的时空变化规律。研究很少涉及在全球变暖大背景下, 近期青藏高原整体的植被覆盖特征^[32-33], 尤其是不同时间尺度植被覆盖动态变化, 青藏高原降水量、气温与植被覆盖之间究竟呈正相关还是负相关等问题不甚明确。本文利用时空分辨率相对较高的MODIS(Moderate-resolution imaging spectrometer)NDVI卫星产品以及青藏高原区域101个气象站点的降水、气温资料, 采用最大合成法(Maximum Value Composites, MVC)、趋势性分析以及相关分析方法, 旨在分析在全球气候变暖背景下, 近17年来青藏高原植被覆盖在不同时间尺度上的空间分布和时间变化规律, 探讨降水和气温变化对植被覆盖的影响, 为高原地区环境保护、生态安全屏障建设提供理论依据。

1 研究区、资料和方法 1.1 研究区域

青藏高原(图 1)位于我国西南部, 在我国境内部分西起帕米尔高原、东至横断山脉、南自喜马拉雅山、北迄昆合山祁连山北侧; 东西长约2945km, 南北宽达1532km, 跨新疆、青海、甘肃、西藏等6个省区, 总面积为2.62×10⁶km², 占我国陆地总面积的26.8%^{[28, 34]}。青藏高原作为世界上平均海拔最高的高原, 在各种地理因素的共同作用下, 东南部属于暖湿性气候, 西北部属于干冷性气候。高原气候基本特点是太阳辐射强烈, 气温低, 气温日较差大, 年变化较小, 干湿季节分明, 气温和降水自东南到西北呈现递减的趋势。

1.2 资料来源

归一化植被指数(Normal Difference Vegetation Index, NDVI)通过红外与近红外波段的组合实现对植被信息状态的表达, 可以用来表征植被覆盖程度和植被生长状况、叶面积指数、生物量以及吸收的光合有效辐射等植被参数^[35-37]。NDVI可以对植被生长动态变化进行监测, 同时能够在较大时空尺度上客观反映植被覆盖程度和植被生长状况^[38-40], 其时序数据也已成为基于生物气候特征开展大区域植被的基本手段^[41], 归一化植被指数在全球及区域大尺度植被变化与气候相互作用中已被广泛应用。

文中NDVI数据来自美国LP DAAC数据中心(Land Processes Distributed Active Archive Center, http://lp-daac.usgs.gov/main.asp)MODIS仪器提供的2000—2017年2月逐月NDVI数据, 分辨率为1km, 数据格式为HDF-EOS, 投影方式为Sinusoidal。MODIS数据是植被遥感监测领域中重要的数据源之一。文中对遥感数据的合成、再投影处理是通过投影软件MRT(MODIS Reprojection Tool)、遥感影像处理软件ENVI(The environmental for visualizing Images)实现。把下载的数据依据青藏高原范围进行拼接, 使其投影为标准的经纬度网格坐标, 将软件输出的MODIS资料保存为GEOTIF格式; 在ENVI里完成青藏高原矢量边界裁剪, 获得研究区的NDVI数据。为反映研究区不同时间尺度植被覆盖状况, 计算月、季节、生长季和年的NDVI平均值, 将数据转存为ASCII格式, 以便用于植被指数与气象要素的相关计算。文中使用GIS(Geographic Information System)绘制植被空间分布和变化趋势图, IDL绘制曲线图。

气象数据来源于中国气象科学数据共享服务网(http://cdc.cma.gov.cn), 选取2000年1月—2016年12月青藏高原地区101个气象站的月降水和温度数据, 气象站包括青海和西藏的全部站点; 新疆的塔什库尔干、甘肃的乌哨岭、玛曲和合作; 云南的德钦、贡山和中甸; 四川的石渠、若尔盖、德格、甘孜、色达、道孚、阿坝、马尔康、红原、小金、松潘、巴塘、新龙、理塘、稻城、康定、得荣、木里和九龙。

1.3 分析方法

文中将研究时段划分为年际、生长季、季节、月等4种时间尺度。季节划分采用气象学上的标准：春季3—5月, 夏季6—8月, 秋季9—11月, 冬季从12月到次年2月, 生长季选取4—9月。采用最大合成法(Maximum value Composites, MVC)获得各研究时段(月、季节、生长季和年)内最大NDVI值作为不同时段的NDVI值, 该方法在一定程度上能够消除云对植被指数的影响。由于文中关注植被覆盖及其变化趋势, 没有单独考虑云、积雪、水体和无数据区域的植被状况。对植被覆盖的等级划分采用等间隔等分方法^[32], 认为植被不同等级之间的变化速率一致, 依据植被变化速率大小划分为严重退化、中度退化、轻度退化、基本不变、轻度改善、中度改善以及明显改善7个等级。

NDVI随时间的变化趋势利用趋势线分析^[42], 即对每个栅格的NDVI数值进行一元线性回归分析。slope为NDVI回归方程的变化速率, 即斜率, 若slope＞0, 则该像元17年间的植被覆盖变化趋势增加, 且数值愈大植被覆盖度增加趋势愈明显; 反之表明植被变化趋势减少; slope=0表示此像元的植被覆盖无变化。采用最小二乘法法拟合NDVI和气候要素随时间的变化速率, 即NDVI和气候要素x随时间t的线性回归。在分析植被指数和气候要素关系时, 先通过插值计算出各观测站点上的NDVI, 再计算年NDVI和年降水量、年平均气温的相关系数。在变化趋势或相关系数分析中, 如果相关系数大于显著性检验表中自由度N对应的临界值, 说明变化速率或相关系数通过指定显著性水平P值的检验。

2 结果分析 2.1 NDVI空间变化

图 2为2000—2016年青藏高原不同季节、生长季和年合成平均植被指数的空间分布图。可以看出, 无论季节、生长季或全年NDVI, 青藏高原东南部植被状况明显好于高原西北部。由于受到东亚夏季风的影响, 东南部能够得到较好的水热条件, 分布着常绿阔叶林以及针叶林, 而青藏高原西北部基本属于干旱、半干旱区域, 水热条件较差, 分布着草原及荒漠, 表明植被覆盖的分布格局与该区域水热条件的时空分布保持了较好的一致性。

春季高原地区植被覆盖整体状况比较差, 植被状况最好的区域(NDVI＞0.8)位于西藏东南部的林芝、昌都等地, 这些区域分布着比较茂密的森林, 相对于其他地表类型具有更高的NDVI值; 植被覆盖较好的区域(NDVI＞0.4)位于青藏高原东部和东南部, 即青海省东部、川西高原以及西藏东南部; 高原大部分区域的NDVI均小于0.4, 尤其在高原西部和北部柴达木盆地, NDVI甚至小于0.2, 基本对应着青海西部和藏北辽阔的荒漠和半荒漠区域。

随着夏季和雨季来临, 高原植被状况有所好转, NDVI＞0.4的区域由高原东南部向高原西北部逐渐增加, 青海省东部和南部、川西高原以及西藏中东部NDVI＞0.6, 尤其是川西高原和西藏东南部NDVI＞0.8。秋季高原雨季结束, 气温逐渐下降, 植被覆盖与夏季植被空间分布比较相似, 但是随着温度降低, 植被趋于枯黄, 高原东南部的NDVI较夏季有所降低, 基本处于0.4—0.8之间, NDVI大于0.8的区域明显减少。冬季高原地区异常寒冷, NDVI大于0.4的区域急剧南缩, 主要位于川西高原东部和南部以及西藏东南部, NDVI＞0.8的区域基本位于西藏东南部。由于青藏高原植被的生长期较平原地区短, 植被生长最快时段基本集中在夏季。因此, 高原生长季和全年植被覆盖的空间分布与夏季植被覆盖特征比较接近, 仅在个别区域存在着差别, 比如生长季和年植被覆盖NDVI＞0.8的区域较夏季范围偏大。

图 3为青藏高原2000—2016年不同研究时段植被的变化趋势。由图可见, 高原植被既有增长又有退化趋势, 植被覆盖变化具有明显的季节和区域差异。春季高原大部分区域植被呈增加趋势, 增加地区分布在青海东部和南部、川西高原北部、西藏东南部, 而西藏那曲中东部、川西高原南部退化比较严重。夏季植被增加区域位于高原北部以及青海省大部, 植被覆盖减少区域位于西藏中西部, 尤其以藏北最为明显。秋季高原植被减少最为明显, 减少区域集中在高原东南部以及青海省南部。冬季高原大部植被覆盖呈增加趋势, 减少区域位于青海省东北部和南部。生长季和年植被覆盖变化空间分布与夏季特征接近。整体上青藏高原植被覆盖改善区域大于退化区域, 改善面积达到高原面积的72%—75%, 退化面积占25%—28%。

根据年平均NDVI变化斜率数据及其分布规律, 采用等间隔等分法^[32]将变化斜率分为7个等级(表 1)。结果表明, 严重退化的区域主要位于青藏高原西南部, 尤其是西藏那曲中东部, 退化面积为9.09km², 占总面积的0.0003%;中度退化区域分布于西藏中西部以及青海南部, 退化面积为136.37km², 占总面积的0.0045%;轻度退化区域范围相对较广, 几乎分布于除西藏东南部之外的大部分区域, 退化面积为3424.32km², 占总面积的0.113%;基本不变的区域面积为2711859.93km², 占总面积89.4894%;轻度改善面积319561.61km², 占总面积10.5453%;中度改善面积2239.44km², 占总面积0.0739%;明显改善面积139.40km², 仅占总面积0.0046%。由此可见, 近年来高原地区退耕还林、草原生态补偿机制等生态环境保护措施的实施, 使得高原地区植被中度和轻度改善面积增大, 青藏高原大部分区域的植被状况有了明显的改善, 但仍然存在植被的不同程度退化区域。

2.2 NDVI时间变化

依据青藏高原月植被指数统计出的近17年NDVI的月平均值(图 4)可以看出, 随着季节更替NDVI相应地发生改变。由于高原冬半年时段长, 降水较少且天气寒冷, 1—3月NDVI值较小, 基本维持在0.15左右, 4月以后随着气温升高和降水增加, 植被开始返青, 尤其以高原东南部较为显著, NDVI值也迅速增大, 不同类型植被生长旺盛, 在8月NDVI值达到最大值(0.3393), 9月后随着农作物的成熟, 大部分植被开始落叶, 植被覆盖逐渐减少, NDVI值急剧下降, 至次年2月下降为全年中的最低值(0.1572)。2000—2016年多年平均的标准误差为0.0021, 最大误差值出现在6月, 仅达到0.0038, 说明随着夏季到来, 植被覆盖呈现明显变化, NDVI反映的植被状况较实际情况略微偏大, 其他月份误差较小, 其中最小误差出现在1月(0.0014)。由此可见, 高原植被指数月际变化的标准误差均非常小, NDVI能很好地反映呈“单峰型”的年变化规律。

图 5显示了青藏高原地区不同季节、生长季和年NDVI随时间的变化趋势。从年均值曲线及其线性趋势上看, 植被指数呈显著的上升趋势, 增幅为0.007/10a, 相关系数为0.619, 通过显著性水平检验(P＜0.01)。2000年的NDVI值为0.337, 为植被指数最低值, 而后NDVI持续上升, 2005年达到0.3641, 以后2006—2012年经历了幅度较小的波动后逐渐上升, 2012年达到近17年来NDVI最大值(0.3726), 随后呈现微弱的下降趋势, 2015年NDVI值为0.3559。不同季节的NDVI值变化曲线表明, 植被指数均呈明显上升趋势, 春、夏、秋、冬季和生长季的NDVI增幅分别为0.01/10a、0.008/10a、0.003/10a、0.007/10a和0.008/10a, 对应的相关系数分别为0.686, 0.629, 0.242, 0.672和0.621。青藏高原地区不同时段NDVI值呈幅度较小的增加趋势, 除秋季外, 其他季节和生长季植被变化均通过显著性水平检验(P＜0.01)。

2.3 气候因子变化特征

多年月平均降水量和平均气温(表 2)表明降水量和气温呈现基本一致的单峰型变化, 即1—3月降水量较小, 气温较低, 最低温度出现在1月, 12月降水最少; 4月以后无论是气温还是降水量均呈显著增加, 降水量和气温都在7月份达到一年中的最大值。结合NDVI月际变化和表 2可以看出, 降水量、平均气温和NDVI具有类似的年变化特征, NDVI值的最大值略微滞后于降水量和气温的峰值。

近17年来不同时段的降水量和气温变化趋势表明, 冬季降水量和温度最低, 夏季温度最高, 生长季降水量最大。由降水量的变化趋势可见, 除夏季降水量呈减少趋势(4.09mm/10a)外, 其他季节降水量均呈增加趋势, 以秋季和春季变化较为明显, 增幅分别为2.51mm/10a和1.99mm/10a, 但是总体上降水量增加的变化趋势不显著, 没有通过显著性水平检验(P＜0.01)。冬季、春季和秋季降水量的年际波动较小, 夏季和生长季降水量年际波动较大。不同研究时段的气温均呈增加趋势, 尤其以春季增加最为显著, 可以达到0.63℃/10a;秋季、夏季和生长季次之, 增温幅度分别达到0.46、0.43℃/10a和0.43℃/10a;冬季增温幅度最小, 仅为0.14℃/10a。除冬季气温外, 其他季节气温变化趋势均通过显著性水平检验(P＜0.10), 尤其是春季和生长季平均气温通过P＜0.01显著性水平检验。不同时段平均气温的年际波动总体较小, 其中, 冬季和秋季的变化相对较大。

对青藏高原地区2000—2016年101个气象站观测的降水量和平均气温实测资料分析(图 6)发现, 近17年青藏高原地区年平均气温变化幅度较大, 多年年平均气温为4.56℃; 气温从2000年至2009年呈在波动之中显著增加, 2009年之后到2012年持续下降, 随后急剧增加; 平均气温最低值出现在2000年(3.79℃), 气温最高值(5.15℃)出现在2016年; 整体上年平均气温呈显著上升趋势, 增幅为0.45/10a, 相关系数为0.613, 通过显著性水平检验(P＜0.01)。在全球气候变暖的大背景下, 区域气候呈现显著的变暖趋势。由年均降水量可见, 近17年青藏高原地区降水变化幅度不大, 多年平均年降水量为485.83mm; 降水量从2000年至2006年呈减少—增加—减少的变化趋势, 并在2006年达到近17年来的最小值(432.11mm), 随后呈现波动变化, 2015年达到次低值437.13mm, 2016年达到最大值(537.37mm); 降水量总体上呈缓慢增加趋势, 增幅为2.04mm/10a, 相关系数为0.032, 未通过显著性水平检验(P＜0.01)。综上所述, 青藏高原气候变化呈现“暖湿化”趋势。

2.4 植被覆盖与气候因子的响应

植被指数与气候因子关系密切, 并且气候条件能够调节和改变植被覆盖的分布状况^[43-44]。从2000—2016年多年平均的高原地区年最大合成NDVI、年降水量、年平均气温的突变规律来看, 降水量的变化不存在显著突变现象, 平均气温和NDVI则分别在2003年和2004年表现出了显著增加的突变现象, NDVI的突变时间滞后于平均气温的增加, 表明近年来平均气温与植被覆盖呈基本一致的增加趋势。

将青藏高原多年平均的各像元的年最大合成NDVI插值到各气象站点上, 得到了NDVI和降水量、气温的散点分布图。无论年降水量还是年平均气温, 总体上均与年最大合成NDVI呈较好的正相关, 通过P＜0.10显著性水平检验。年降水量较年平均气温在空间尺度上与植被指数有更好的对应, 降水量增加可极显著地促进高寒草甸植被的生长。温度升高在一定程度上可以满足植被正常生长所需的热量供应, 适当增温可促进高原植被生长, 但是温度过高时, 过高的温度会改变植物群落的小气候环境, 直接或间接影响到植被的生长发育和生物量生产, 对植被覆盖起到抑制作用^{[24, 26]}。

为了解植被指数对降水量、平均气温在青藏高原不同区域的响应, 分别计算了各高原观测站年降水量和年平均气温与对应植被指数的相关系数(图 7)。结果表明, 年降水量和植被指数在青藏高原大部分区域呈正相关, 尤其是青藏高原东北部, 包括青海大部、西藏大部以及川西高原北部, 体现了降水对草原植被生长有利^[22-23]; 负相关区主要位于青海东南部、西藏北部和东南部以及川西高原南部。年平均气温和植被指数的正相关区主要位于青海大部、西藏西北部和川西高原南部; 负相关区位于高原西南部, 尤其是西藏中东部, 说明气温对植被的生长在特定区域具有不利的影响。

由此可见, 对高原西部降水较少和气温较低的区域(以西藏阿里为例), 降水量和热量条件都会影响植被生长, 并且与植被均呈正相关; 对于高原东南部降水相对充沛的部分区域, 降水量和热量条件对植被生长影响不显著, 甚至降水偏多, 可能还会抑制植被生长^[24]。从不同类型下垫面来看, 草原和荒漠覆盖区域的NDVI与降水量、气温表现为较密切的正相关^[11], 尤其以草原更为显著, 通过P＜0.10显著性水平检验; 森林区NDVI与降水量、气温关系不显著。总体来说, 青藏高原降水量、气温分别与植被指数既存在正相关, 也存在负相关, 降水量对植被覆盖的影响大于平均气温对植被覆盖的影响。

3 讨论

综上所述, 无论不同季节、生长季或年最大合成的植被状况, 青藏高原植被均呈现逐渐增加的趋势。但是, 有研究^[45]应用增强的植被指数(Enhanced Vegetation Index, EVI)分析2000—2012年青藏高原生长期植被覆盖变化, 认为青藏高原植被呈现不显著的减少趋势, 因为植被减少的速率大于植被增加的速率, 并且气候变化的初始状态会导致不同植被类型对气候变化的敏感性的差异。文中研究结果与此结论存在着研究时段以及植被增加、减少速率定量化分析上的差异, 没有涉及气候初始状态对植被覆盖变化的影响。降水量、气温与植被覆盖存在着正、负相关的区域, 气温与植被覆盖的正相关区主要位于高原东北部, 负相关区位于西南部, 降水量对植被覆盖的影响较气温明显, 这与以往研究^{[24, 46]}结论接近一致, 并且指出了植被指数和气候要素的相关存在着区域性差异。由于青藏高原观测站点主要位于高原中东部, 高原西部大多数区域没有站点或者站点很少, 研究结论的适用性局限于高原中东部。文中关注的气候因子主要限于降水量和气温的站点资料, 今后需要采用格点数据、卫星估算降水以及相对湿度等其他资料, 进一步开展青藏高原不同区域、尤其是高原东、西部植被类型差异特征及其对气候变化的响应和互馈作用。

4 结论

本文分析了2000—2016年青藏高原地区植被的时空特征和变化趋势, 得出以下结论：

(1) 青藏高原东南部植被状况明显好于高原西北部, 植被覆盖的分布格局与区域水热条件的时空分布保持了较好的一致性。近年来青藏高原大部分区域的植被状况有明显改善, 严重退化区域主要位于青藏高原西南部。除秋季外, 青藏高原地区NDVI值在2000—2016年呈幅度较小的增加趋势, 各时段植被变化均通过显著性水平检验(P＜0.01)。

(2) 除夏季降水量呈减少趋势(减幅为4.09mm/10a)外, 其他季节降水量呈增加趋势; 不同研究时段的气温均呈增加趋势, 整体上青藏高原气候变化呈现“暖湿化”趋势。降水量和平均气温总体上与年最大合成NDVI呈较好的正相关, 但也存在负相关区域, 青藏高原上降水量对植被覆盖的影响较平均气温明显。