青藏高原高寒草地植被物候时空变化特征

黄文洁，曾桐瑶，黄晓东

（兰州大学草地农业科技学院 / 草地农业生态系统国家重点实验室，甘肃 兰州 730020）

物候学是研究生物圈动植物与气候、水文、土壤等环境因子之间相互关系的科学[1]。植物物候是指植物由于受气候和环境因素驱动而出现的以年为单位的周期现象[2-3]。物候变化直接影响生物生产，并在生态系统生产力和碳循环中发挥关键作用，许多学者形象地认为物候是季节性和内生性的最敏感的综合指标。近100年来，在变暖已成为事实的背景下[4]，植物物候因其对气候变化响应的敏感性已经成为生态系统对气候变化反馈的重要感应器和全球变化中重点研究的对象[5]。

传统的物候观测采用野外人工目视观察法，该方法耗费大量人力物力，并且不易开展大范围长时间序列的监测。遥感观测具有覆盖范围广、连续的空间、长时间序列的特点，能够定量反映植物的季节生长和发育及其年际间的年际变化等特征，而成为传统物候观测手段的有效补充[6]。越来越多的学者以遥感技术为手段深入研究植被物候变化，使其逐步发展为一门独立的学科并成为全球气候变化研究的重要组成部分。宋春桥等[7]利用2001-2010年MODIS EVI数据对藏北高原典型样区物候的时空变化进行了分析。田柳茜[8]选用GIMMS NDVI (1982-2006)、SPOT NDVI (1999-2013)和 MODIS NDVI(2001-2013)数据集探讨了青藏高原返青期的时空特征。李鹏[9]选择GIMMS NDVI(1982-2012)、 SPOT NDVI(1999-2012)和MODIS NDVI(2001-2012)3种植被指数数据来研究青藏高原植物枯黄期的时空转变。近年来，关于草地植被物候时空变化特征在其他地区已有大量的研究，但针对青藏高原地区高寒草地植被物候的研究，尚存在地面验证数据较少，时空尺度较小，以及不够深入等问题。因此，本研究基于地面观测数据以及遥感数据，分析青藏高原高寒草地的植被物候特征，以期为后期研究物候对气候的响应做好准备。

NDVI植被指数由于对叶绿素敏感，曲线变化幅度大，可消除部分太阳高度角和观测条件带来的噪声等特点，使其成为物候研究中广泛而有效的选择[10]。利用遥感数据提取植物物候期主要有阈值法、最大斜率法、滑动平均法和经验公式法等方法[11]。这些数学方法各有优劣，如丁明军等[12]利用1999-2009年SPOT VGT-NDVI数据采用最大比率法对高寒草地物候的时空变化进行了分析。司文才和刘峻明[13]选用2000-2009年NDVI数据，采用动态阈值法研究了河北省冬小麦(Triticum aestivum)的返青期，得出冬小麦返青期由南到北逐渐推迟的结论。侯学会等[14]利用Logistic模型对NDVI累积曲线进行拟合，依据曲率极值方法提取植被物候，研究了近10年东北地区的森林物候。

青藏高原是全球气候变化响应较为敏感的地区，其热液条件处于生物生长极限，生态系统极其脆弱[15]。本研究针对青藏高原高寒草地，选用2001-2015年MODIS-NDVI数据，基于TIMESAT软件提取植被返青期、枯黄期及生长季长度3个物候指标，分析其时空变化特征，探究其驱动因子，以期为青藏高原高寒草地物候监测打下基础。

1 研究区概况

图1 青藏高原草地类型及牧试台站空间分布Figure 1 Distributions of the grasslands and climate stations on the Tibetan Plateau

青藏高原位于我国西南部 73°19′-104°47′ E，26°00′-39°47′ N(图 1)，是世界上海拔最高的高原，有世界“第三极”的美誉[16]。青藏高原长约 2 800 km，宽 300～1 500 km，总面积达 257.24 万 km2，占我国陆地面积的26.8%[17]。青藏高原区域地形复杂多样，地势落差大，海拔由西北向东南逐渐降低。全年气温低、区域差异大，降水少且分布不均匀，东南多、西北部少，空气稀薄，具有很强的太阳辐射。由于海拔的原因使得青藏高原地区自然条件与其他中纬度暖温带和亚热带地区大有差异，青藏高原的草地类型主要有高寒草甸、高寒草原、山地草甸等17种。青藏高原的气候变化(水热、干湿)间接地影响了青藏高原植被的生长过程，而植被的生长反过来也会影响碳收支、气候变化、土地覆盖、农作物估产及放牧等经济活动[18]。研究青藏高原物候时空变化特征能够揭示全球变暖背景下高寒草地生态碳循环及其对气候变化的响应[19]。

2 数据与方法

2.1 数据来源

本研究采用2001-2015年青藏高原归一化植被指 数 (normalized difference vegetation index, NDVI)分析青藏高原高寒草地物候变化趋势。NDVI能够很好地反映植被的生物量、初级生产力以及植被覆盖状况，在植被生长和趋势变化的研究中得到广泛应用[20-21]。本研究所采用的MODIS(moderate resolution imaging spectroradiometer)V006 MOD13A1 NDVI数据集，时间为2001-2015年，空间分辨率为500 m，时间分辨率为16 d，每年23期。DEM数据从地理空间数据云(http://www.gscloud.cn/)中下载得到，空间分辨率为500 m。

研究采用了青海省气象科学研究所提供的共计20个牧试站点的地面牧草物候观测资料，用于验证遥感提取的物候参数(表1)。站点观测记录表明，2003-2015年，青藏高原牧草返青期主要集中在每年4月初至6月中旬(4月10日-6月12日)，返青期在空间上表现出显著的差异，不同区域不同草地类型的返青期相差最大达到2个月。枯黄期主要集中在8月末至10月初(8月29日-10月14)，相比返青期，青藏高原牧草进入枯黄期的时段相对比较集中。

2.2 时间序列重建

16 d 最大合成法 (maximum value composite, MVC)虽然消除了部分云、大气等对NDVI的影响，但NDVI时间序列仍然存在噪声，气溶胶、云层等噪声很可能影响提取植被物候的准确性[22]。因此，在提取物候参数之前需要对NDVI数据进行平滑处理。TIMESAT 3.2平台提供了3种时间序列重建的方法：非对称高斯函数拟合法(AG)[23]、双Logistic函数拟合方法(DL)[24]和Savizky-Glolay过滤算法(S-G)[25]。S-G算法最符合原始NDVI曲线的细节，在降低噪声的同时保持曲线的范围与形状不发生偏移[26]，能够真实地反映出青藏高原NDVI数据的周期性变化[27]，因此，本研究选用S-G过滤算法拟合NDVI时间序列数据，S-G滤波窗口的值设为4。

表1 牧试台站经纬度与海拔Table 1 Latitude, longitude, and altitude of the pasture test stations

2.3 物候参数提取

植被的 3 个参数，生长季开始 (the start of growth season, SOG)、生长季结束 (the end of growth season,EOG)和生长季长度 (length of growth season, LOG)被用来反映植物的物候变化。 采用儒略日 (Julian days)换算方法，将物候出现日期转化为距1月1日的实际日数，即年序列累积日数，得到各物候期的时间序列。基于遥感技术的SOG是指大多数植物开始扩大叶片并正常生长的日期，此后植被进入快速生长阶段。EOG指大多数植物不能正常生长并且叶子开始变黄的日期，之后植物进入休眠或死亡阶段。LOG等于EOG减去SOG，这与传统概念上生长季节是指“一年中植物生长的天数”不同。目前，针对时间序列的遥感数据发展了不同的物候期提取方法，包括阈值法[28-29]、曲线拟合法[30]、最大变化速率法[31]和经验回归法[32]等，但未对最佳提取方法形成普遍共识[33]。结合该区域的前期研究，本研究采用阈值法提取青藏高原高寒草地植被物候期。其中，当NDVI增加至拟合函数左半部分振幅20%时定义为植物返青期，当NDVI降低至拟合函数右半部分振幅60%时定义为植物枯黄期[9]，生长季的长度则为返青期与枯黄期时间之差，据此利用TIMESAT软件逐年提取青藏高原高寒草地植被物候参数。

2.4 趋势分析方法

为了研究青藏高原高寒草地近15年植被物候时空变化特征，本研究采用趋势分析法来拟合2001-2015年物候的年际变化[34-35]。年变化率(Slope)大于零表示返青期推迟、枯黄期推迟或生长季延长的趋势，年变化率小于零表示返青期提前、枯黄期提前或生长季缩短的趋势。公式如下：

式中：y为第i年的物候；i为年变量，从1年到n年。

3 结果与分析

3.1 精度验证

利用20个牧试站提供的青藏高原牧草物候期观测数据作为真值，采用回归系数(R2)、平均误差(Bias)和均方根误差(RMSE)对MOD13A1产品提取的物候期进行精度验证。2003-2015年20个验证站点MOD13A1识别的返青期和枯黄期与地面观测返青期和枯黄期的精度评价结果 (图2)表明，2003-2015年，青藏高原高寒草地牧试站观测返青期主要集中在第110-160天，而MOD13A1遥感数据提取的返青期主要集中在第100-161天。遥感识别的返青期与地面观测值的R2= 0.61，Bias=-1.23，RMSE= 8.31，误差超过 10 d 的占 19.6%。2003-2015年，青藏高原高寒草地牧试站观测枯黄期主要集中在第250-290天，而 MOD13A1 遥感数据 提取的枯黄期主要集中在第240-288天。遥感识别的枯黄期与牧试站观测值的R2= 0.64，Bias=5.64，RMSE=8.39，误差超过 10 d 的占 22.2%。返青期和枯黄期的验证结果表明，基于MOD13A1产品提取的物候参数可以较为准确地反映青藏高原高寒草地返青期与枯黄期的时空变化特征。

图2 MOD13A1产品反演与地面观测的返青期和枯黄期对比Figure 2 Comparison between the phenology observed from the ground stations and phenology based on the MOD13A1 retrieval

3.2 物候期的空间分布

通过栅格处理计算得出青藏高原高寒草地2001-2015年的物候空间分布特征图(图3)。青藏高原高寒草地返青期存在明显的东早西迟的趋势，从东南到西北，返青期开始时间逐渐推迟。多年平均返青期主要集中在儒略年第110-170天，这与牧试站地面观测数据基本一致。青藏高原高寒草地平均返青期的时间为第137天，东部返青期主要在第130天左右，而西部则主要在第150天左右开始。返青期在高原东南部地势较低的河谷地区最早，在纬度或海拔较高的区域，牧草地返青期较晚。

青藏高原高寒草地枯黄期在空间分布上也存在明显的差异，由中部向西南和东南方向逐渐推迟，高原中部枯黄期在第260-280天，西部边缘枯黄期迟于东部，约在10月中下旬枯黄。从枯黄期空间分布特征图来看，整个青藏高原高寒草地的平均枯黄期时间为第286天左右，草地集中在9月初到10月底进入枯黄期。在返青期和枯黄期的共同作用下，牧草生长季长度自东南至西北逐渐缩短，这一趋势符合研究区地势和气候的特征。生长季长度集中在第100-170天，平均生长季长度为 149 d。

3.3 物候期的年际变化

通过线性拟合青藏高原2001-2015年高寒草地物候期的年际变化趋势(图4)。返青期的趋势并不是固定的，在不同的时间段内趋势不尽相同。2001-2009年间，返青期有明显的提前趋势，返青期提前 18 d；2009-2015 年间，返青期推迟 4 d。总体而言，青藏高原高寒草地的返青期呈现出提前的趋势，2009年草地的返青期最早。青藏高原高寒草地植被枯黄期呈现波动提前的趋势，这一趋势较之返青期而言并不明显，2002、2008和2013年牧草进入枯黄期较早，2001-2015年间青藏高原高寒草地植被枯黄期整体提前5 d。将2001-2015共15年分为两个时间段，在2001-2008年期间枯黄期以每年0.857 d的速度提前，在随后的2008-2015年又出现了推迟的趋势，每年推迟0.286 d(图4)。从整个研究期来看，2001-2015年枯黄期呈现出微弱的提前趋势。在青藏高原高寒草地返青期和枯黄期年际变化的共同影响下，青藏高原高寒草地植被生长期呈现波动延长的趋势，平均生长季长度延长了9 d，其中2007年的平均生长季长度最长，达到158 d。

图4 2001-2015年青藏高原高寒草地植被物候年际变化Figure 4 Inter-annual variation of the vegetation phenology in the alpine grasslands of the Tibetan Plateau from 2001 to 2015

选取青藏高原山地草甸类、高寒草原类、高寒草甸类3种典型的高寒草地类型，分析其物候期差异。青藏高原不同草地类型的物候差异显著(图5)，结果表明，山地草甸类植物返青期较早，返青期集中在第124-134天；高寒草甸类的植物返青期次之，集中在第132-141天；高寒草原类的植物返青期最晚，在第142-152天。不同草地类型的枯黄期相比返青期时间相对比较集中，其中高寒草原类与高寒草甸类在第270-282天进入枯黄期，山地草甸类在第278-288天进入枯黄期。在返青期与枯黄期的共同影响下，山地草甸类的生长季长度最长，范围在145～162 d，高寒草原类的生长季长度最短，范围在123～140 d，高寒草甸类草地的生长季长度在133～148 d。2001-2015年间，不同草地类型的返青期都有提前的趋势，但趋势并不显著，而枯黄期高寒草原类与高寒草甸类表现出一定的提前趋势，山地草甸类基本没有表现明显的提前或推迟趋势。

图5 青藏高原不同草地类型物候年际变化特征Figure 5 Inter-annual variation of the vegetation phenology in the different grassland types of the Tibetan Plateau

2001-2015年，青藏高原高寒草地物候年际变化在空间分布上也呈现出显著的差异(图6)。牧草返青期中呈提前趋势的占71.71%，呈推迟趋势的占28.29%，在青藏高原高寒草地东部和西北地区返青期呈现出推迟的趋势。其中，返青期显著提前的占5.12%，显著推迟的占2.11%。牧草枯黄期中呈提前趋势的占65.18%，呈推迟趋势的占34.82%，枯黄期提前的区域集中在青藏高原高寒草地中部。枯黄期显著提前的占4.05%，显著推迟的占2.81%。从牧草生长季长度变化来看，具有增长趋势的占57.81%，具有缩短趋势的占42.19%，整体呈现延长趋势。生长季长度中约2.67%的区域显著提前，约6.15%的区域显著推迟。

3.4 物候期的海拔差异

以100 m间隔为海拔梯度提取不同草地类型的物候参数，对物候与海拔建立线性拟合关系，反映出2001-2015年青藏高原高寒草地物候受海拔因素影响的趋势。结果表明，青藏高原2001-2015年高寒草地植被物候受海拔因素影响(图7)。在海拔2 000-3 500 m，返青期呈现不规则波动，海拔对不同草地类型的返青基本没有任何影响；在3 500-5 100 m，返青期随海拔的升高而推迟，在 5 100 m以上，随着海拔的升高返青期反而提前。枯黄期以5 200 m为界，出现两种不同的趋势，在2 000-5 200 m范围内，随着海拔升高枯黄期呈现波动提前的趋势；在5 200 m以上，枯黄随着海拔的升高而推迟。对于生长季长度而言，以3 500和5 000 m为节点，小于3 500 m生长季长度变化不规则，海拔在 3 500-5 000 m，生长季长度缩短。海拔大于5 000 m 时，生长季长度反而延长。在海拔 5 000 m以上出现了植物物候特征与低海拔相反的变化趋势，可能是由于受积雪覆盖地表的影响，造成草地植被物候提取出现偏差。

综上所述，海拔因素很大程度上影响了青藏高原高寒草地物候期。当海拔在3 500-5 000 m时，海拔每升高500 m，牧草返青期推迟4.3 d，枯黄期提前 1.3 d，生长季长度缩短 5.6 d。在海拔 3 500 m以下，牧草返青期和生长季长度与海拔并无明显的相关关系(图7)，这可能是受到人类活动的影响[13]。在海拔5 000 m以上，植被物候随海拔变化的关系呈现出相反的趋势，因为在5 000 m以上的区域，主要以寒漠或苔原为主，该区域植被生命活动短暂，对物候参数进行提取基本是不可行的，造成分析结果可信度较差。

青藏高原不同草地类型对海拔变化的响应有所不同，本研究选取了高寒草甸类、高寒草原类、山地草甸类3种草地类型来分析其在不同海拔高度植物的返青期、枯黄期和生长季长度的变化。分析结果表明，在 3 500-5 000 m 海拔区间，3 种草地类型的返青期都随着海拔的升高而推迟。其中高寒草原类的推迟幅度最大，达到9.41 d。枯黄期随海拔变化的趋势不明显，波动较大，山地草甸类与高寒草甸类与海拔呈负相关的关系，说明随着海拔的升高，枯黄期逐渐提前，但高寒草原类随着海拔的升高却呈现枯黄期推迟的趋势，R2达到 0.513 6。生长季长度在海拔 3 000-5 000 m 区间，都表现出随着海拔的升高而缩短的趋势，且趋势显著(图8)。

4 讨论与结论

图7 2001-2015年青藏高原高寒草地物候期与海拔的关系Figure 7 Relationship between the phenology and altitude of the alpine grasslands of the Tibetan Plateau from 2001 to 2015

图8 青藏高原不同草地类型物候期与海拔的关系Figure 8 Relationship between the phenology and altitude of the different grassland types of the Tibetan Plateau

青藏高原高寒草地植被物候的变化趋势不是单纯的线性变化，在不同的时间段，有不同的趋势，呈现出的变化非常复杂。因此，推测这是由于物候提取方法不同，即物候期是由多种因素共同决定导致的，如水热因素和人为因素等。丁明军等[12]认为，1999-2009年间，青藏高原高寒草地返青期整体上呈提前趋势，枯黄期呈推迟趋势，生长期长度延长，在重合的时间段本研究也得出相同的结果。田柳茜[8]选用GIMMS (1982-2006)、SPOT (1999-2013)、MODIS (2001-2013)3种NDVI数据提取青藏高原植被返青期，尽管在绝对数值上不相同，均得出返青期提前的趋势。马晓芳等[36]选用1982-2005年GIMMS NDVI数据提取青藏高原高寒草地物候，得出返青期和枯黄期逐渐提前，生长季长度缩短的结论。吕灿宾[37]认为1980-2010年间生长季结束时间呈提前趋势。受到时间尺度和不同遥感数据源精度的影响各研究得出的物候变化存在一定的差异，整体上，本研究结论与前人的研究基本吻合[7-9]，但在分析物候期与海拔之间关系方面，本研究得出了一些更有意义的结论。物候期并非随着海拔的增加而呈线性变化，以5 000 m为分界点，在海拔3 500-5 000 m 间，草地物候与海拔之间有明显的相关性，随着海拔的升高牧草返青期推迟，枯黄期提前，生长期长度缩短。当海拔大于5 000 m之后，寒漠或苔原受极低温度的影响，不具有明显的物候期[14]。

不同的研究方法将影响草地关键物候参数的提取，本研究采用阈值法，提取青藏高原高寒草地植被物候期。由于遥感数据分辨率较低和存在云层、大气的噪声等问题，并且数据预处理方法与物候参数提取方法都将影响研究结果，所以很有必要结合地面实测数据对遥感数据提取的物候期进行精度评价。例如，Chen等[38]利用中国科学院的地面实测数据对1982-1993年中国东部温带植被物候期进行了验证。本研究采用2003-2015年的牧试站物候观测值对遥感数据进行均方根误差与偏差评价。验证结果表明，地面观测数据与遥感监测结果之间的误差在可接受范围内，本研究提取的青藏高原高寒草地植被物候期参数的精度较高，可以反映青藏高原的草地植被物候的时空变化。

本研究利用MOD13A1产品分析了青藏高原高寒草地2001-2015年植被物候时空变化特征并探究了海拔对物候期影响。结果表明：1)青藏高原高寒草地物候在空间上存在较大差异性。由于受地形抬升和水热条件恶化的影响，呈现出由东南向西北返青期逐渐推迟，枯黄期逐渐提前，生长季长度逐渐缩短的趋势。2)不同草地类型物候期存在明显差异，山地草甸类植物返青期较早，高寒草原类植物返青期较晚。高寒草甸草原类的枯黄期最早，山地草甸类的枯黄期最迟。3)青藏高原高寒草地植被物候年际变化的趋势表现为返青期和枯黄期提前，整体上，生长季长度增长。具体分析结果表明，2001-2015年间，植被返青期提前了 14 d，枯黄期提前 5 d，生长季长度延长了 9 d。4)青藏高原高寒草地物候空间异质性与海拔关系密切。存在着 3 500 和 5 000 m 的分界线，3 500 m以下物候随海拔变化规律不明显，在3 500-5 000 m之间，二者关系显著，随着海拔的升高，返青期推迟，枯黄期提前，生长季长度缩短。