新疆地区草地植被物候时空变化

张仁平，郭靖，冯琦胜，梁天刚*

(1.绿洲生态教育部重点实验室,新疆大学干旱生态环境研究所，新疆 乌鲁木齐 830046； 2.兰州大学草地农业生态系统国家重点实验室,兰州大学草地农业科技学院，甘肃 兰州 730020；3.新疆林业科学院，新疆 乌鲁木齐 830000)

植被物候现象可以直观地反映植被年际与季节环境条件，不仅能有效地表达植物生长发育的季节特征，而且能直观表达植物群落和物种的生理学和形态学特征[1]。植被物候变化虽然是循序渐进的过程，但对环境的影响却极为深远，其会影响到植被净初级生产力、全球的气候和碳循环[2]。同时，草地牧草返青期和枯黄期的变化，会明显影响当地草地畜牧业的发展[3]。

随着全球变化研究的不断深入，植被物候已成为气候变化对陆地生态系统影响的研究热点。遥感数据可以大范围、长时间和空间连续地反映植被物候特征。近年来，越来越多的国内外学者选择遥感数据作为估算陆地植物物候变化的手段[4-5]，并且已取得了丰硕的成果[6-9]。自20世纪80年代以来，北半球植物生长季开始日期呈提前趋势，结束日期呈推迟趋势[10-12]；在国内，利用遥感数据对植物物候分析结果也和北半球植物变化趋势基本一致[13-14]。然而也有研究表明，北半球一些地区的植物生长季开始日期自20世纪90年代呈推迟趋势[15]。

目前，基于遥感数据提取物候信息的主要方法有：阈值法、滑动平均法、最大斜率法、倒数法、经验公式等方法[16-19]。这些数学模型和方法都有各自的优劣，如Reed等[11]基于AVHRR-NDVI的数据采用中值滤波法提取美国的植被物候信息，虽然大大地消除了云对物候信息的影响，但是损失了植被最大绿度的信息。Chen等[20]利用SG滤波法对中国区域的SPOT-VGT数据进行平滑处理，提取了高质量的物候信息数据。Zhang等[21]利用Logistic函数对美国东北地区的MODIS-EVI时间序列数据进行拟合，由于这种方法不必界定植物物候阈值，可以普遍适用于大范围的物候监测，美国NASA根据这种算法形成物候产品(MCD12Q2)，已被广泛用来研究全球地面的季节性物候学和年际变化[22-23]。

新疆地区具有独特的自然地理环境，属于典型的干旱和半干旱地区，生态系统极为脆弱，草地面积位居我国第三位，对气候变化的影响极为敏感。过去几十年，新疆地区正处于剧烈升温阶段，且明显高于全球的增温水平[24-25]。已有研究表明，新疆地区物候正发生显著变化[26]，但对新疆草地物候变化的研究仍处于相对薄弱状态。因此，准确掌握新疆地区草地植被物候信息，分析草地植被物候的变化趋势对保护和合理利用当地草地资源具有重要的意义。

1 材料与方法

1.1 研究区概况

新疆维吾尔自治区位于我国西北部，地理位置介于34°22′-49°33′ N，73°22′-96°21′ E，总面积为166×104km2，约占国土总面积的1/6(图1)。新疆地处欧亚大陆腹地，四面高山环抱，北有阿尔泰山，南有昆仑山系，中有横亘全境的天山，三山环抱中为广袤的准噶尔和塔里木盆地，“三山夹两盆”构成了新疆独特的地理环境特征。新疆属于典型的温带大陆干旱性气候，光热资源充足，年日照时数达2550～3500 h，年平均气温9～12 ℃，无霜期长达180～220 d，降水量稀少，北疆年降水为100～200 mm，南疆在100 mm以下。而蒸发量则相反，北疆为1500～2300 mm，南疆为2100～3400 mm。由于特殊的地理位置、地形条件和干旱气候的影响，新疆生态环境极为脆弱，植物种类稀少，覆盖度低，类型结构简单。新疆草地总面积居我国第三位，约57.26万km2，可利用草地面积约48万km2，占新疆国土面积的34.4%，新疆草地面积是耕地面积的15倍，是森林面积的22倍，占全区绿色植被面积的86%。

1.2 物候数据来源

物候数据来源于美国NASA的MODIS土地覆盖动态产品(MCD12Q2)，通过网站(https://wist.echo.nasa.gov/)获得2001-2014年MCD12Q2产品，空间分辨率为500 m，每年6幅影像。MCD12Q2产品利用每年8 d增强性植被指数(enhanced vegetation index，EVI)生成，EVI通过地表反射率(nadir bidirectional reflectance distribution function-adjusted reflectance)数据计算，空间分辨率为500 m，由于EVI数据比归一化植被指数数据(normalized difference vegetation index，NDVI)提供了更大的动态范围，因此选择EVI生成土地覆盖动态产品。

图1 研究区草地类型及地面返青观测站点空间分布Fig.1 Location of grassland types and grassland ground spring phenology observation sites

MCD12Q2产品已被广泛用来研究全球地面的季节性物候学和年际变化[22]，陆地覆盖动态产品是基于Zhang等[21]的算法(Logistic函数拟合法)，该算法不需要预先定义阈值并进行数据平滑，在一定程度上减少了主观因素的影响。根据MCD12Q2产品提供的周期，是基于MODIS时间序列数据的EVI曲率变化率过渡数据：返青为EVI增长的节点日期，特指季节开始的日期(start of the growing season，SOG)；枯黄为EVI最小的节点日期，特指季节结束的日期(end of the growing season，EOG)；生长季长度指季节结束日期与季节开始日期的差值(length of the growing season，LOG)。

利用MRT(MODIS ReProjection Tools)软件对MCD12Q2产品数据进行格式转换和投影转换，把HDF格式转换为Geotiff格式，投影转换为WGS84。

1.3 精度验证

在2012-2014年2-5月期间，每隔5～7 d对采样点完成一次观测，共获取123个样地，每个样地大小为500 m×500 m，每个样地里面设置5～8个样方，当50%的样地里草地返青时即认为该样地返青[1]。采用皮尔森相关系数(R)，平均误差(Bias)和均方根误差(root mean square error，RMSE)对MCD12Q2物候产品进行精度验证。

1.4 空间趋势分析

Theil-Sen Median 趋势分析和Mann-Kendall(MK)检验方法结合起来，已成为分析长时间序列数据趋势变化的重要方法之一，并且已逐渐应用到植被长时间序列分析中[27-29]。其中，Theil-Sen Median趋势分析可以减少数据异常值的影响，是一种可靠的非参数统计的趋势计算方法[30-31]。

2 结果与分析

2.1 产品验证

图2 MCD12Q2物候产品反演与地面观测的返青期对比与验证Fig.2 The comparison and verification between phenology observed in the ground stations and phenology based on MCD12Q2 retrieval

图2给出了3年(2012-2014年)123个MCD12 Q2识别的返青与地面物候观测值的精度评价结果。可以看出，2012-2014年地面实测点返青期主要集中在每年的第68～122天，而MCD12Q2物候产品的返青期主要集中在第57～128天。遥感识别的返青期与地面观测值的R=0.76，显著性水平P<0.05，Bias=-2.84 d，RMSE值为16.44 d。基于MCD12Q2物候产品的返青期验证结果表明，遥感识别比地面监测的物候期提前，MCD12Q2物候产品基本可以反映新疆草地的返青期。

2.2 物候多年平均值的空间分布特征

分析新疆草地植被各物候参数2001-2014年平均值的空间分布特征(图3)，可以看出新疆草地各物候参数明显呈区域性和垂向差异，由低山区到高山区，SOG逐渐推迟。SOG主要出现在第90～170天的区域占新疆草地的74.1%，SOG最早的区域主要位于准噶尔盆地、伊犁河谷和塔里木盆地，在一些海拔较高的区域，如阿尔泰山、天山以及昆仑山区域，SOG晚于第170天。在北疆地区的EOG较早，大部分区域草地的EOG早于第290天，南疆地区的EOG较迟，也表现出海拔上的差异，草地EOG出现在第250～300天的区域占新疆草地面积的89.9%。LOG明显表现出海拔上的差异，在阿尔泰山和天山以及昆仑山草地的LOG较长，塔里木盆地周边区域的LOG较短，草地LOG出现在100～200 d的区域占83.3%。分析草地物候标准差可知，北疆草地SOG、EOG和LOG的标准差较南疆草地的小，而山区草地较平原草地的小。

新疆草地物候与海拔关系密切(图4)。SOG随海拔的升高而推迟，EOG随海拔的升高具有提前的趋势，LOG与EOG的分布规律基本一致。海拔每升高1000 m，SOG推迟13 d，EOG提前7 d，LOG缩短20 d。从图中可以看出，在海拔1500 m以下区域，SOG和LOG随海拔变化的规律不明显，波动比较大；在海拔3000 m以下区域，EOG随海拔变化的规律不明显，在海拔4500 m以上，EOG随海拔变化的关系特别密切。此外，新疆不同草地类型的物候明显不同(图5)，新疆不同草地类型的SOG在第100～136天，EOG在第256～291天，LOG为122～190 d。其中：温性荒漠类、温性荒漠草原类、低地草甸类和温性草原化荒漠类的SOG较早，在第100～103天；温性草原类、温性草甸草原类和山地草甸类的SOG次之，在第107～116天；高寒草原类、高寒草甸类和高寒荒漠类的SOG最晚，在第127～136天。高寒荒漠类的EOG最早(第256天)，低地草甸的EOG最迟(第291天)。高寒荒漠类和高寒草甸类的LOG较短，为122～129 d；而温性草原化荒漠类、温性荒漠类和温性荒漠草原类的为173～177 d，低地草甸类的LOG最长，为190 d。

图3 2001-2014年新疆草地物候多年均值和标准差的空间分布Fig.3 Spatial pattern of mean and standard deviations of Xinjiang grassland phenology during 2001-2014

2.3 物候的年际变化特征分析

2.3.1物候区域尺度上的年际变化 为了从整体上把握新疆草地物候提前和延迟的动态变化趋势，结合Theil-Sen趋势分析和Mann-Kendall检验方法，借助于Matlab2012a软件进行新疆草地物候逐像元分析并绘制研究区像元尺度上的Theil-Sen分布图和Mann-Kendall中值检验图(图6)。结果表明，新疆草地SOG提前的面积比例(57.5%)大于推迟的(42.5%)，其中SOG显著提前的区域占研究区总面积的4.6%，SOG提前但不显著的区域占52.9%；SOG显著推迟的区域占2.6%，推迟但不显著的区域占39.9%。SOG提前的区域主要分布在阿尔泰山中山带和伊犁河谷中山带，SOG推迟的区域零星分布在准噶尔盆地边缘地带以及天山北坡的区域。EOG提前的面积比例(47.8%)略小于推迟的(52.2%)，其中显著提前的区域占2.6%，提前但不显著的区域占45.2%；EOG显著推迟的区域占研究区总面积的2.9%，推迟但不显著的区域占49.3%。草地枯黄期提前的区域零星分布在伊犁河谷中山带以及准噶尔盆地。新疆草地生长季长度延长的面积比例(52.3%)略大于缩短的(47.7%)，其中LOG显著延长的区域占研究区总面积的3.4%，延长但不显著的区域占48.9%；显著缩短的区域占3.1%，缩短但不显著的区域占44.6%。

图4 2001-2014年新疆草地物候年平均值与海拔的关系Fig.4 The relations between grassland phenology and altitude of the Xinjiang during 2001-2014

2.3.2物候年际变化在草地类型的分异特征 图7反映了新疆主要草地类型的物候变化速率。新疆大部分草地类型返青期呈提前趋势，平均提前速率为0.11 d·yr-1，其中温性草甸草原类和山地草甸类的返青提前得最早，分别达0.56和0.53 d·yr-1，而温性草原化荒漠类和温性荒漠类的返青期呈推迟趋势，分别达0.02和0.04 d·yr-1。绝大多数草地类型的枯黄期呈推迟趋势，平均推迟速率为0.14 d·yr-1，低地草甸类的枯黄推迟趋势最明显，达0.70 d·yr-1，温性草甸草原类枯黄期呈提前趋势，达0.22 d·yr-1。所有草地类型的生长季长度呈延长趋势，平均延长速率为0.25 d·yr-1，低地草甸类的生长季长度延长最长，达0.75 d·yr-1。

图5 2001-2014年间新疆不同草地类型物候年平均值Fig.5 Phenological mean values of different grassland types in Xinjiang during 2001-2014 1：温性草甸草原类Temperate meadow steppe；2：温性草原类Temperate steppe；3：温性荒漠草原类Temperate desert steppe；4:高寒草原类Alpine steppe；5:温性草原化荒漠类Temperate steppe desert；6:温性荒漠类Temperate desert；7:高寒荒漠类Alpine desert；8:低地草甸类Lowland meadow；9:山地草甸类Mountain meadow；10:高寒草甸类Alpine meadow；11:所有类型All.下同 The same below.

图6 2001-2014年新疆地区草地植被物候变化特征Fig.6 Spatial distribution on inter-annual variation of phenology in the Xinjiang during 2001-2014

3 讨论

北半球温带是全球气候变化最显著的区域之一，研究该区域的植被物候特征，对深入认识气候变化对陆地生态系统的影响具有典型性。目前已有许多学者利用不同遥感数据对北半球的植物物候进行了研究。Zeng等[32]利用MODIS数据分析表明，北半球植被返青期在2000-2010年呈提前趋势，枯黄期呈推迟趋势。Kaufmann等[33]利用AVHRR数据分析了欧亚大陆40°-70°地区植被物候变化特征，得出返青期呈提前趋势，枯黄期有推迟趋势，本研究得出的植被返青期提前，枯黄期推迟的结论与此一致。就生长季变化的幅度而言，Delbart等[34]利用1982-2004年AVHRR和SPOT-VGT数据分析欧亚大陆北纬50°-72° N植被物候的变化，结果也显示返青期提前。已有的研究结果得出SOG提前幅度在0.28～9.00 d·10 yr-1，EOG推迟幅度在0.6～14.0 d·10 yr-1[33-34]。本研究得出返青期提前0.11 d·yr-1，枯黄期推迟为0.14 d·yr-1，这跟前人的研究结果基本一致。不同的研究结果得出物候变化幅度不同，可能与研究区域、时间以及遥感数据不同有关。

IPCC第五次(2013年)评估报告指出，近30年全球气温持续增高[35]，而气候变化会导致陆地植被改变[36]。近年来，新疆的气温在冬季呈降低趋势，而在春季、夏季和秋季呈增加趋势，冬季温度的降低有利于植被的返青，由于植物的返青必须经过一定的低温过程，在温度最低的时期，温度越低，越有利于植物的返青[37]。像大多数温带和寒冷地区的植物，新疆草地在冬季经历一个休眠的过程，这可以避免霜冻对植物敏感生长组织的伤害。然而，春季的高温会加速热量的积累，这有助于草地返青的提前，植被的枯黄期与温度关系也很大，夏季和秋季的高温会推迟草地枯黄期的到来[38-40]。降水也对植物物候期产生一定程度的影响，但不同植被类型物候期对降水的响应明显不同[37]。

图7 2001-2014年间新疆草地生长季物候变化趋势Fig.7 Trend of grassland phenology in Xinjiang during 2001-2014 years

到目前为止，植被物候的遥感分析主要采用AVHRR，MODIS和SPOT等遥感数据，由于每种数据波段宽度、空间分辨率、传感器、最大值合成等不同，都会导致结果存在一定的差异，而且在大气校正方面：AHVRR数据不仅考虑云污染及水汽的影响，同时考虑了火山气溶胶散射的影响，而MODIS和SPOT数据在此基础上还考虑了气体的吸收，例如，MODIS数据还考虑了臭氧吸收，SPOT数据不仅考虑了臭氧吸收，还考虑了水蒸气和其他气体的吸收。Fontana等[41]利用AVHRR、SPOT和MODIS数据对高寒草地物候做了分析，表明SPOT和MODIS数据要优于AVHRR数据。Zhang等[21]对MODIS-EVI时间序列数据进行了季节性拟合，提取了植物物候变化趋势，这在一定程度上降低了噪声水平。

目前，由于遥感数据的分辨率较低，且受遥感数据本身的质量、数据预处理方法和物候识别方法影响，使得基于遥感数据获取的植物物候期的可靠性评价尤为必要。目前，大部分研究采用地面观测的物候数据对遥感识别的物候数据进行验证。例如，Fisher等[42]利用地面物候观测进行验证，量化评价了两种尺度下的植物物候遥感监测精度，发现植物物候遥感监测结果的均值可以反映精细尺度到粗尺度的统计转换。Yu等[39]利用NOAA AVHRR NDVI数据研究1982-2006年青藏高原植被春季物候时，采用均方根误差与平均绝对误差(mean absolute error，MAE)评价了地面观测数据与遥感监测结果之间的误差。Chen等[1]利用AVHRR NDVI数据研究1982-1993年中国温带的植被生长季开始日期和结束日期时，采用中国科学院的地面数据进行了验证。本研究通过2012-2014年的物候观测数据对MCD12Q2产品进行了验证，这可以让更多的学者掌握MCD12Q2物候产品的相关信息。

4 结论

本研究基于2001-2014年的MCD12Q2时间序列影像，分析了新疆草地物候SOG、EOG和LOG关键物候参数的空间分布格局，同时，结合Theil-Sen Median 趋势分析和Mann-Kendall检验方法，分析了该地区草地物候在近年来的时空变化状况。

新疆草地物候多年均值由低海拔到高海拔呈明显的区域性差异。随着海拔的不断升高，SOG逐渐推迟，EOG逐渐提前，LOG逐渐缩短，海拔每升高1000 m，SOG推迟13 d，EOG提前7 d，LOG缩短20 d。此外，不同草地类型的物候明显不一致。新疆不同草地类型的平均SOG在第100～136天,平均EOG在第256～291天，其中高寒荒漠类的EOG最早(第256天)，平均LOG为122～190 d。

2001-2014年间，新疆草地SOG整体上呈提前趋势，提前的面积比例(59.5%)大于推迟的(40.5%)，提前速率为0.11 d·yr-1。EOG呈推迟趋势，推迟的面积比例(52.7%)略大于提前的(47.3%)，推迟速率为0.14 d·yr-1。LOG呈延长趋势，延长的面积比例(52.9%)略大于缩短的(47.1%)，延长速率为0.25 d·yr-1。