基于TINDVI的伊犁河谷草地覆盖变化

2021-04-27闫俊杰沙吾丽达吾提拜刘海军

水土保持研究 2021年3期

闫俊杰, 陈晨, 赵阳, 沙吾丽·达吾提拜, 刘海军

(伊犁师范大学资源与生态研究所, 新疆伊宁 835000)

草地占据了全球将近20%的的陆表面积[1]，它不仅是陆地生态系统重要组成部分[2]，而且在生物多样性保护、水源涵养、大气调节及原材料供给等生态系统维护与调节方面占有重要地位[3]。全球变化背景下，草地植被动态变化及其对生态环境的影响被广泛关注[4-7]，而对草地植被生长科学而精确的量化表达是研究草地动态变化的基础。

遥感技术在数据获取的时效性、成本及数据丰富度等方面具有难以比拟的优势[8-9]，是目前区域植被动态监测的主要方法及数据来源。基于遥感技术，不同学者将植被反射光谱中的不同波段进行线性或非线性组合来突出地表植被信息，提出了多种植被指数来对植被的生长状况进行量化表达[10]，为植被动态变化研究提供了有效且计算简便的指标。在众多植被指数中，归一化植被指数(Normalized difference vegetation index，NDVI)应用最为广泛，被广泛应用于全球、区域及流域等不同尺度的植被动态研究[4-5,9,11-14]。MODIS NDVI，NOAA-AVHRR NDVI及SPOT NDVI时间序列数据是主要数据源[15]。在对NDVI时间序列数据预处理过程中，最大值合成法(maximum Value Composites，MVC)常被用于计算年内NDVI最大值(NDVImax)，来表征单个生长季内植被生长达到最好时的状况[10,16]，从而为植被的年际变化分析提供基础。然而植被的年际变化并不仅仅体现在不同年份植被生长可达到的最高水平的差异上，植被在不同年份的年内生长过程的差异也是植被动态变化的重要方面。时间累积NDVI(Time integered NDVI，TINDVI)被定义为年内NDVI变化曲线与时间轴围城的面积或一定时期内NDVI的累积值[17-19]。相对于NDVImax，TINDVI包含了植被生长过程的信息，且已有研究证明TINDVI与植被干物质量具有很好线性关系[20-21]，然而以TINDVI为衡量指标来分析植被变化的研究却鲜有报道。对比分析TNDVI与NDVImax在反映植被时空变化方面的比较优势，探讨TINDVI的应用价值，将为全面和深入认识植被的动态变化提供有利帮助。

西北地区是我国草地资源的重要分布区，而位于天山西段的伊犁河谷更是我国久负盛誉的优质牧场。已有学者以NDVImax或覆盖度为参考，指明伊犁河谷草地退化严重，且退化仍在持续恶化[22-25]。本研究利用包含生长过程信息的TINDVI作为表征草地植被的参数指标，分析伊犁河谷草地植被TINDVI时空变化，并探讨TINDVI与NDVImax的差异，以期进一步探明伊犁河谷草地植被的变化特征，并为相关研究提供借鉴。

1 研究区概况

伊犁河谷位于80°09′42″—84°56′50″E，42°14′16″—44°53′30″N，地处天山西段，整个河谷地势西低东高，地形呈向西敞开的“V”字形(图1)，西风湿润气流受此特殊地势地形影响，容易抬升凝结，给河谷带来丰沛降水，成就了其“塞外江南”的美誉。伊犁河谷气候虽属于温带大陆性气候，但却表现出明显的高山气候特征[26]，平原区与山区气候明显不同，平原区干旱少雨、蒸发能力强，而山区降水丰沛、天气冷凉[27]，平原区降水200 mm左右，但山区降水可达1 000 mm，气温变幅约为9.2～2.8 ℃[26]。河谷内植被类型以草地主，且类型丰富多样，是新疆乃至我国著名的优质牧场，分布有高寒草甸、山地草甸、温性草甸草原、温性草原、温性荒漠草原、温性荒漠及低平地草甸等多种类型草地[24,28]。伊犁河谷行政区划上属于伊犁哈萨克自治州，畜牧业在国民经济中占有重要比例，且畜牧方式仍然以传统放牧为主。相关研究表明，近年来受气候变化及人为影响，草地植被覆盖度持续降低[29]，草地生产力降低及毒害草蔓延等草地退化问题严重[22,30]。

图1 研究区示意图

2 材料与方法

2.1 数据来源与预处理

文章用到的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(NDVImax)作为年尺度上的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数据栅格大小一致。

2.2 研究方法

2.2.1 TINDVI计算 TINDVI通常被定义为生长季内NDVI变化曲线与时间轴围城的面积或NDVI的累积值[17-19]。累积值计算方便且简单，因此TINDVI利用生长季内NDVI的累积值表示，即：

(1)

式中:i即代表NDVI数据的日期。MOD13Q1 NDVI产品每年23期数据,从每年的第1天到第353天每16 d一期数据。根据研究区草地植被返青期和枯萎日期实际情况及其空间差异[31];TINDVI计算取3月初(第65天)到11月初(第305天)时段内共16期NDVI数据,并在计算过程中对3月初和11月初NDVI<0数据进行了剔除。

2.2.2 单调趋势Mann-Kendall非参数检验 Mann-Kendall非参数检验被广泛应用于时间序列数据趋势检验，其统计量Zc计算过程如下[32-33]：

(2)

(3)

(4)

(5)

式中:xi为xk为时间序列数据的集合。在α置信水平上,若|Zc|≥Z1-α/2,即时间序列在α置信水平上存在显著上升或下降趋势。Zc>0,则变化趋势为上升趋势;Zc<0,则变化趋势为下降趋势。

时间序列变化速率可用Kendall倾斜度来表示，其计算公式如下：

(6)

式中:1

3 结果与分析

3.1 TINDVI空间分布特征

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

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

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

经统计，整个河谷6.08.0的比例达到了21.67%；TINDVI<2.0的比例最少，为9.34%；2.0

3.2 2000-2018年伊犁河谷草地TINDVI时间动态

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

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

由图3可知，2000—2018年伊犁河谷草地TINDVI总体呈波动增加的变化趋势，但根据Mann—Kendall趋势检验，其Zc值为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年分别出现两个明显的波动增加时段。

3.3 2000-2018年伊犁河谷草地TINDVI空间动态

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

注：1.96为Mann—Kendall检验p=0.05显著水平上Zc临界值。

由图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%，分布也均较为离散。

3.4 TINDVI与NDVImax年际变化的时空差异

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

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

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

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

4 讨论

4.1 TINDVI与NDVImax年际变化时空差异的原因

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

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

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

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

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

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

4.2 TINDVI及NDVImax与降水和气温关系及其对比

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

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

注：UU为与降水及气温相关性均非显著；PU为与降水显著正相关；UP为与气温显著正相关；NU为与降水显著负相关；UN为与气温显著负相关；PP为与降水及气温均显著正相关；NN为与降水和气温均显著负相关；PN为与降水显著正相关，与气温显著负相关；NP为与降水显著负相关，与气温显著正相关。

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

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