文章用到的NDVI数据为美国NASA EOS数据中心发布的MODIS MOD13Q1产品。MOD13Q1数据为16 d合成的时间序列数据,空间分辨率为250 m,每年23期,本研究用到2000年2月—2018年12月共19 a数据。在数据预处理过程中,对NDVI时间序列数据进行了镶嵌、转投影、重采样、研究区裁剪以及Savitzky-Golay滤波处理。用于年际变化分析的NDVI年数据通过对每年23期数据的最大值合成(MVC)处理获得^[10,16],即用年内植被生长达到最好状况时的NDVI(NDVI_max)作为年尺度上的NDVI。

草地分布边界的矢量数据是通过对2018年6—7月份sentinel-2假彩色合成影像目视解疑获得。对解疑获得的矢量数据进行了栅格化处理,作为NDVI数据研究区裁剪的掩膜。为保障数据的空间匹配,草地分布数据的栅格化处理及NDVI数据重采样处理过程中,像元大小均设置为250 m×250 m。

气温和降水数据来自中国气象科学数据共享服务网(http:∥data.cma.gov.cn/),包括伊犁河谷内伊宁市、霍尔果斯市、霍城等10个气象站和其周边巴音布鲁克在内的11个气象站的月数据。在11个气象站点数据的基础上,借助于Anusplin4.2软件及DEM数据,进行空间插值,获得伊犁河谷月气温和月降水空间分布数据,之后合成生长季(3—10月份)内的平均气温和累积降水空间分布数据,插值过程中栅格大小设置为250 m,与NDVI数据栅格大小一致。

利用2000—2018年伊犁河谷草地TINDVI的年空间分布数据,计算其多年平均值。基于TINDVI的多年平均数据分析伊犁河谷草地TINDVI的空间分布特征(图2)。

由图2可知,伊犁河谷草地多年平均TINDVI的值域约为0～11.0。空间上伊犁河谷草地TINDVI呈现明显的垂直分异特征。具体表现为TINDVI先随海拔的升高而逐步增加,河谷西部海拔较低的平原及山前洪积冲积扇区TINDVI较低(2.0<TINDVI<4.0),到中山区域后TINDVI达到最大(TINDVI>8.0),但空间上主要集中在河谷中南部及东部的中山区域,之后随海拔升高TINDVI变化转为逐步减小,河谷周边积雪时间较长的高山区TINDVI最小值(TINDVI<2.0)。

经统计,整个河谷6.0<TINDVI<8.0的面积比例最高,为34.24%,空间上主要分布在河谷西北和乌孙山南部的低山区域、昭苏盆地—特克斯谷低、河谷东部和南部中山向高山过度区以及河谷南部高山区的山谷地带; TINDVI>8.0的比例达到了21.67%; TINDVI<2.0的比例最少,为9.34%; 2.0<TINDVI<4.0的比例为14.97%; 4.0<TINDVI<6.0的比例为19.78%,空间上主要分布在特克斯河下游两侧的低海拔区以及伊犁河北部的低山丘陵区。

图2 草地TINDVI空间分布(A)及其直方图(B)

对2000—2018年伊犁河谷草地TINDVI的年空间分布数据求逐年平均,绘制TINDVI的年际变化曲线(图3),分析伊犁河谷草地TINDVI总体的年际变化。

图3 2000-2018年TINDVI年际变化曲线

由图3可知,2000—2018年伊犁河谷草地TINDVI总体呈波动增加的变化趋势,但根据Mann—Kendall趋势检验,其Z_c值为1.05,未达到显著水平,其年变化率β值为0.015。在年际波动上,TINDVI于2008年和2014年出现较大谷值,值分别为5.34,5.36; 于2007年、2013年和2016年出现较大峰值,值分别为6.26,6.40,6.39; 2008年谷值前后的2000—2007年及20008—2013年分别出现两个明显的波动增加时段。

基于2000—2018年伊犁河谷草地TINDVI的逐年空间分布数据,根据Mann-Kendall检验的计算步骤,进行逐像元计算,获得TINDVI变化趋势Z_c值和变化率β值的空间分布数据(图4),基于此分析TINDVI变化的空间特征。

图4 TINDVI变化趋势Zc值(A)和变化率β值(B)的空间分布

由图4可知,伊犁河谷绝大部分草地TINDVI呈增加趋势,经统计,其面积比例达到了73.95%,但其中仅有11.56%达到了显著水平,在空间上离散分布在河谷南部、东部及中西部的部分区域; 整个区域有26.05%的区域TINDVI呈降低变化趋势,空间上主要分布在昭苏盆地、特克斯河下游周围及河谷北部的中山区域,而TINDVI降低达到显著水平的比例仅为0.67%。

对于TINDVI的变化速率,全区大部分区域的变化率为-0.025～0.050(图4B),其比例达到了85.91%,其中TINDVI变化率为0～0.025的分布比较广泛,比例达到了41.29%,TINDVI变化率为-0.025～0的比例较小,为19.95%,空间上主要离散分布于河谷中部及北部的中山和低山区域,TINDVI变化率为0.025～0.050的比例为24.67%,空间上离散分布于河谷的不同区域; TINDVI变化率<-0.025和>0.050的比例均小,比例分别为5.99%和8.10%,分布也均较为离散。

相关研究表明,近年来伊犁河谷草地的NDVI_max及覆盖度变化总体均呈现减小趋势^[24-25,29],但根据本文的分析,伊犁河谷草地TINDVI却总体有所增大,空间上73.95%草地的TINDVI呈现不同程度增加。鉴于本文研究与前人研究得出相反的结论,有必要对TINDVI及NDVI_max变化的差异性进行探讨。

对比TINDVI与NDVI_max年际变化曲线(图5A),2000—2018年NDVI_max变化总体呈减小趋势,据Mann-Kendall趋势检验,其Z_c值为-0.28,未达到显著水平,变化速率为-0.000 1; 而TINDVI总体呈增加趋势,Z_c值为1.05,也未达到显著水平。在年际波动上,TINDVI以及NDVI_max均在2008年及2014年出现较大谷值,在2011年出现较大峰值,但在2005—2006年、2009—2010年及2016—2017年,NDVI_max呈现增大变化,而TINDVI则呈现减小变化。

变化的空间分异方面,利用TINDVI变化趋势Z_c值图与NDVI_max变化趋势Z_c图进行叠加,对比TINDVI及NDVI_max变化趋势的空间差异(图5B)。由图5B可知,伊犁河谷57.28%的草地NDVI_max有所减小,42.58%草地NDVI_max和TINDVI的变化趋势相反,其中NDVI_max减小而TINDVI增加的比例达到了36.91%,空间上主要位于NDVI_max及TINDVI值较高的中山区域,而NDVI_max和TINDVI呈相同变化趋势的面积比例虽然也达到了57.42%,但在空间上主要分布在NDVI_max及TINDVI值相对较低的山前洪积冲积扇、低山丘陵及积雪时间较长的高山区域。

图5 伊犁河谷草地TINDVI与NDVImax年际变化曲线(A)及变化趋势的叠加图(B)

由上文分析可知,2005—2006年、2009—2010年以及2016—2017年,全区平均TINDVI和NDVI_max出现相反的波动变化(图5A),同时不同区域内,TINDVI和NDVI_max的年际变化也不尽相同(图5B)。NDVI_max代表了生长季内植被覆盖到达最大时的NDVI值,而TINDVI代表生长季内NDVI的累积值。生长季内NDVI变化曲线即能直接展示NDVI_max值的大小,也能间接展示TINDVI值的大小,因此,对比不同年份及不同区生长季内NDVI变化曲线,能为解答NDVI_max及TINDVI时空变化产生差异的原因提供答案。

以16 d为间隔,计算生长季内(每年第65—305天)TINDVI与NDVI_max变化趋势相反年份(2005—2006年、2009—2010年以及2016—2017年),以及研究期前期(2000—2005年)和末期(2014—2018年)TINDVI与NDVI_max变化趋势不同叠加类型区的区域NDVI平均值,制作生长季内TINDVI与NDVI_max变化趋势相反年份(图6A)及研究期前期和末期两者变化趋势不同叠加类型(图6B)的NDVI变化曲线,探讨NDVI_max及TINDVI时空变化产生差异的原因。

对于在2005—2006年、2009—2010年及2016—2017年NDVI_max与TINDVI呈现截然相反的变化,由图6A可知,2005年草地NDVI在达到最大值之前春夏季的第65—161天和之后夏秋季的第209—289天内其值均小于2006年同时间的值,致使NDVI_max虽有所增加,但TINDVI却发生减小; 2009—2010年主要是在春夏季的第65—161天内2009年NDVI值小于2010年; 2016—2017年则主要是在夏秋季的第193—257天内2016年NDVI值小于2017年。

图6 TINDVI与NDVImax变化趋势相反年份及两者变化趋势不同叠加类型在研究期前期和末期NDVI生长季内的变化

而对于NDVI_max及TINDVI变化趋势在空间上的差异,图6B分别展示了NDVI_max及TINDVI变化趋势不同交互类型(图5B)中2000—2018年的前5 a(2000—2004年)及末端5 a(2014—2018年)的NDVI平均年内变化曲线。由图可知,对于NDVI_max及TINDVI均呈减小趋势的区域,主要是由于NDVI_max减小的同时,在夏秋季的绝大部分时间内(第129—289天)NDVI均有所减小(图6B1); 对于NDVI_max减小而TINDVI增大的区域,主要是由于在春夏季的第65—145天内及秋季的第257—305天的NDVI均有所增大(图6B2); 对于NDVI_max增大TINDVI减小的区域以及NDVI_max增大TINDVI增大的区域,分别是由于夏秋季的第209—305天的NDVI均有所减小(图6B3)或整个生长季NDVI几乎都有所增大(图6B4)。

通过上述讨论和分析可以看出,生长季内NDVI随植被从返青到枯萎的变化过程也随区域的不同和年份的不同而存在明显差异,单利用植被生长达到最大时的NDVI_max作为基础指标来评估草地的年际变化明显存在一定局限,用TINDVI作为参考指标来衡量草地植被变化具有一定优势。然而,虽然相关研究证实了TINDVI与玉米和大麦等农作的干物质量具有线性关系^[20-21],且TINDVI也常被作为植被生物量的相似或替代指标^[34-35],但产草量及净初等生产力是衡量草地光合作用所积累有机质或干物质的关键指标,也是评价草地退化与否的重要参考^[4-5,28],而TINDVI与其关系如何仍需深入探究。

气候变化是驱动植被变化的关键因素。分别计算TINDVI及NDVI_max与降水和气温的偏相关系数,探讨TINDVI及NDVI_max与降水和气温相关性及其差异(图7—8)。

由图7可知,无论是对于TINDVI还是NDVI_max,降水与草地植被变化的相关性明显高于气温,且与降水以正相关为主,这与前人关于天山及中亚区域的研究结果一致^[13-15]。而由图8可知,2000—2018年生长季内伊犁河谷草地降水以减少为主,而全区气温均有所升高,气温升高促使水分蒸发增强,加剧了水分散失,而降水的减少又减少了水分来源,加剧水分对植被生长的胁迫,这可能是大面积草地NDVI_max减小的原因。但在水分胁迫加剧的条件下,TINDVI却发生增加,而由上文分析可知,TINDVI增加主要得益于春季(3—4月份)植被NDVI较大幅度的增加(图6B2),经分析,2000—2018年3—4月平均气温和累积降水均有所增加(气温和降水Mann—Kendall检验Z_c值分别为0.28及0.98),为春季(3—4月份)NDVI的增加提供了有利条件。

图7 2000-2018年TINDVI与降水和气温(A)及NDVImax与降水和气温(B)相关性空间分布

图8 2000-2018年变化趋势Zc值空间分布

气候是驱动植被变化的关键因素,而海拔又使气候对植被的影响产生明显的海拔分异,因而植被与气候变化的相关性也常存在明显的海拔分异^[4-5]。根据图7A,在伊犁河谷中部低海拔区和中海拔区,植被与降水呈正相关而与气温相关性较弱,而在河谷周边的高海拔区,植被气温呈正相关而与降水呈负相关。伊犁河谷地处亚欧大陆干旱与半干旱区的腹地,水分是植被生长的关键限值因子,尤其在海拔较低的区域,因而低海拔区植被与降水变化呈正相关^[14]; 而高海拔区,降水虽然相对丰沛,但冬季漫长,全年气温明显低于平原区,气温则是限制植被生长的关键因子,因而植被变化气温呈正相关,与降水呈负相关^[14,36]。植被与气候相关性的海拔分异不仅在伊犁河谷有所表现,在整个天山、祁连山和青藏高原等区域均有所表现^[14,4-5,37]。

对比TINDVI及NDVI_max与降水和气温的相关性(图7),经计算,TINDVI与降水和气温呈显著相关的比例分别达到了38.51%和12.22%,而NDVI_max的该比例则分别为22.33%和5.27%,TINDVI与降水和气温变化相关性明显高于NDVI_max,尤其是与气温的相关性,可见TINDVI对气候变化响应更为敏感; 同时,在空间上,相对与NDVI_max,TINDVI与降水和气温相关性的海拔分异规律也更为明显,这也间接表明TINDVI对气候变化响应更为敏感。相对于植被生长达到覆盖最高时的NDVI_max,TINDVI包含了整个生长季内植被生长过程中NDVI变化的信息,而累积降水和平均气温也包含有整个生长季内各自信息,这应该是TINDVI对气温和降水变化响应更敏感的原因。