北疆地区草地TI-NDVI与NDVImax时空异质性评价

焦阿永， 陈伏龙， 闫俊杰， 凌红波， 申瑞华

（1.石河子大学水利建筑工程学院，新疆 石河子 832003；2.伊犁师范大学资源与生态研究所，新疆 伊宁 835000；3.中国科学院新疆生态与地理研究所，新疆 乌鲁木齐 830011；4.中国科学院大学，北京 100049）

气候变化和人类活动带来了全球或区域生态环境的改变。草地是陆地生态系统的重要组成部分，草地变化构成了全球变化的重要组成部分［1］。草地是人类生存的重要自然资源和环境基础，通过对草地变化研究能为人类有效应对全球变化带来的环境问题提供重要参考和支持。而准确且敏感的参考指标是草地变化研究的基础［2］。

异质性是地表生态系统的基本特征，也是评价地表覆被动态要回答的基本问题。为充分理解草地动态的特征或规律，学者们在时间和空间维度上揭示了关于草地异质性的更多细节［3-5］。传统的实地调查的草地监测方法由于耗时和劳动密集型，很难在大尺度上实现对草地的时空监测［6-7］。由于遥感技术在获取数据的时效性、空间性及成本等方面具有无法比拟的优势，是目前区域植被动态监测的主要方法及数据来源［8-9］。借助于遥感数据计算得到的植被指数是量化地表植被状况的有效和实用方法。各植被指数中归一化植被指数（NDVI）与光合有效辐射、生物量、植被生产力等紧密相关，是目前表征植被生长状态及植被覆盖度中应用最多的植被指数［10］。且多家机构或研究学者已制备了NOAA-AVHRR NDVI，SPOT-VGT NDVI，TERRA-MODIS NDVI 和GIMMS NDVI 等多种覆盖全球的NDVI成品数据，以方便使用。MODIS数据因其较高的空间分辨率，并且改进了在红光波段和近红外波段的波幅及辐射定标等技术，使其更详尽地反映地表植被的空间差异，已被广泛应用于区域植被的动态监测［11-14］。

在对NDVI 时间序列数据预处理过程中，年最大值合成法（MVC）常被用于计算年内NDVI最大值（NDVImax），来表征单个生长季内植被生长达到最好时的状况，并作为植被在年时间尺度上的状态指标［15-17］。然而，NDVI在浓密的植被冠层中容易达到饱和，敏感性降低，削弱了NDVI对植被时空异质性的表达能力。同时植被的动态变化不仅存在于植被能达到最佳时的状况上，也体现在植被生长过程的年际差异上，如植被生长季长度的变化等。时间累积NDVI（Time Integrated NDVI，TI-NDVI）被定义为年内NDVI变化曲线与时间轴围成的面积或一定时期内NDVI 的累积值［18-22］。基于地面监测和卫星遥感的多项研究表明，在农作物、灌木、草地和森林等多种生态系统中TI-NDVI与植物生产力具有很强的相关性［23-28］。基于此，TI-NDVI 常被作为生产力的替代指标广泛用于植被动态、生态功能及环境脆弱性等研究之中［29-32］。相对于NDVImax，TI-NDVI不仅积累了植被的NDVImax，还积累了植被在返青期和枯萎期的NDVI，揭示更多关于植被异质性的信息。

鉴于此，本研究选择草地分布广泛且类型丰富的新疆北部（北疆）作为研究区，在16 d 合成的MODIS NDVI 时间序列数据的支持下，借助于GIS 空间分析技术以及Mann-Kendall非参数趋势检验和变异系数（CV）等数理统计方法，研究北疆草地时空动态，分析NDVImax 和TI-NDVI 在表达草地异质性方面的优势比较，以期为北疆地区乃至全球草地有效评估和管理提供参考。

1 研究区概况

研究区域位于新疆北部（北疆，42°8′24″～48°3′N、80°36′36″～95°55′12″E），地形高差悬殊，海拔高度介于-157～7077 m，主要以山地和盆地为主（图1），年均降水量200 mm以下，年均气温-4～9 ℃，属于温带大陆性干旱半干旱气候。区域中部的准噶尔盆地被天山、阿尔泰山以及准噶尔西部山地环绕，西风气流从西南部的伊犁河谷地带进入盆地，后被东部天山阻隔，无法继续东进，水汽停留在原地，给伊犁河谷带来丰富的降水。北疆地区降水量空间分布极不均匀，山区降水量远高于平原和盆地，西南部多于东北部。伊犁河谷、天山和阿尔泰山年平均降水量超过200 mm，属于半干旱区，其他区域属于干旱区。研究区域内除了额尔齐斯河为外流河，其他河流均为内流河，补给水源主要为大气降水、冰雪融水和地下水。

图1 研究区位置及地形Fig.1 Location and topography of the study area

2 材料与方法

2.1 数据来源与预处理

研究用到的NDVI 数据为美国NASA EOS 数据中心（https://ladsweb.modaps.eosdis.nasa.gov）发布的MODIS MOD13Q1 产 品。MOD13Q1 数 据 为16 d 合成的时间序列数据，空间分辨率为250 m，每年23期，时间序列是2000—2019年。在数据预处理过程中，对NDVI 时间序列数据进行了镶嵌、转投影、重采样、研究区裁剪以及Savitzky-Golay 滤波处理［33］。用于年际变化分析的NDVImax 数据是通过对每年23期数据的最大值合成（MVC）处理获得，即用年内植被生长达到最好状况时的NDVI作为年尺度上的NDVI。

北疆草地分布边界的矢量数据是通过对Sentinel-2假彩色合成影像的目视解译获得。Sentinel-2影像选取2018 年6—9 月植被覆盖达到最大，且云量＜10%的高质量影像。假彩色合成影像像元大小为10 m。对解译获得的草地分布矢量数据进行了栅格化处理，作为裁剪NDVI 数据的掩膜。草地分布数据的栅格化处理及NDVI数据重采样处理过程中，像元大小均设置为250 m×250 m，以确保其与NDVI图像的空间匹配。同时NDVImax＜0.1的像元被作为无植被区进行了剔除。

2.2 计算TI-NDVI

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

式中：i代表NDVI数据的日期。MOD13Q1产品每年23期数据，从每年的第1 d到第353 d每16 d一期数据，并在计算过程中对NDVI＜0 数据进行了剔除，避免植被尚未变绿或已经变黄时的NDVI［34］。由于北疆植被生长季通常开始于3 月初，结束于11 月初，为了进一步降低非植被的噪声，消除一年中早于65 d、晚于305 d的NDVI。TI-NDVI计算取3月初（第65 d）到11 月初（第305 d）时段内共16 期NDVI数据。由于NDVI 的取值范围为0～1.0，因此，理论上TI-NDVI的取值范围为［0，16］。

2.3 NDVImax与TI-NDVI的比较

分别对2000—2019 年NDVImax 和TI-NDVI 时间序列影像求平均值，获得NDVImax 和TI-NDVI 多年平均影像，对NDVImax 和TI-NDVI 所揭示的草地空间格局进行差异性比较。并提取NDVImax 多年平均影像所有像素值以及其所对应的TI-NDVI 值，以此计算NDVImax 与TI-NDVI 的相关系数，评估两者相关性。

利用Mann-Kendall 非参数趋势检验方法，对NDVImax 和TI-NDVI 时间序列影像进行逐像元计算，获得2000—2019年NDVImax和TI-NDVI变化趋势空间分布图，分析NDVImax 和TI-NDVI 变化趋势的空间分布。同时，对NDVImax 变化趋势图与TINDVI 变化趋势图进行叠置操作，制作两者变化趋势差异的空间分布图，分析NDVImax 与TI-NDVI 变化趋势的空间差异。

变异系数（CV）用于衡量数据系列中数据点的离散程度。CV越大，说明数据点相对平均值离散程度越大。在用NDVImax 或TI-NDVI 等指标衡量植被动态时，所选指标的CV越大，能揭示的植被波动就越大，能揭示植被波动的细节就更多，对植被动态也更加敏感。因此，分别在时间和空间维度上对NDVImax 和TI-NDVI 的CV进行了计算。时间维度上，分别对NDVImax 和TI-NDVI 时间序列影像进行逐像元计算，获得时间CV（CVt）的空间分布图。空间维度上，分别在NDVImax 和TI-NDVI 多年平均值影像数据的基础上，利用15×15像元的滑动窗口（即移动窗口中包含225个像元）进行逐像元计算，获得空间CV（CVs）的空间分布图。利用CVt和CVs空间分布图，在时间和空间维度上，对比分析NDVImax 和TINDVI 对草地异质性表达的敏感性及两者的比较优势。在计算CVs时，笔者也尝试了5×5像元和9×9像元的移动窗口，其CVs空间分异与15×15像元移动窗口的结果基本一致。

3 结果与分析

3.1 北疆地区NDVImax和TI-NDVI空间格局

在像元尺度上计算2000—2019 年北疆NDVImax 和TI-NDVI 的年平均值，表征NDVImax 和TINDVI的空间格局（图2）。

从图2 可以看出，北疆地区NDVImax 和TI-NDVI的空间格局基本相同，均表现出明显的海拔分异。北疆地处欧亚大陆干旱半干旱区腹地，水是制约植被生长的关键因素，高程控制着水的空间分布［35-38］。阿尔泰山山区、天山山脉和塔城盆地周边山区降水相对丰沛，草本植物生长良好，草地NDVImax 和TINDVI 分别＞0.5 和＞5.0。而在海拔3500 m 以上的山区，由于冬季较长，气温较低，植被的生长受到制约，草地NDVImax和TI-NDVI相对较低［39-40］，一般分别＜0.4和＜3.0，而在准噶尔盆地和伊犁河谷的低海拔地区，由于降水稀少，气候炎热，干旱制约了植被的生长，草地NDVImax 和TI-NDVI 最低，NDVImax＜0.2，TINDVI＜2.0。

图2 北疆年平均NDVImax（a）和TI-NDVI（b）的空间分布Fig.2 Spatial distribution of yearly averaged NDVImax(a)and TI-NDVI(b)in northern Xinjiang

对比图2a和图2b可以看出，阿尔泰山北部和天山西南部的高山区域，TI-NDVI的低值区域（TI-NDVI＜2.0）明显多于NDVImax 的低值区域（NDVImax＜0.3）。这主要是因为此区域生长季短，即使其NDVImax与其他山区差异较小，其累积NDVI却相差很大。表明TI-NDVI 比NDVImax 更能表达高山区域草地空间分布的异质性。

3.2 北疆地区NDVImax和TI-NDVI年际变化特征

图3显示了北疆整个草地NDVImax 和TI-NDVI的年际变化。Mann-Kendall 趋势检验表明，2000—2019 年北疆草地NDVImax 和TI-NDVI 均呈现显著增加趋势（P＜0.01），Mann-Kendall 统计量Z分别为2.75 和2.69。年 际 波 动 上，除2000—2001 年、2005—2006 年、2014—2015 年和2018—2019 年外，NDVImax和TI-NDVI的波动变化也基本一致。

图3 2000—2019年北疆草地区域平均NDVImax和TINDVI的年际变化Fig.3 Inter-annual variations of the regional averaged NDVImax and TI-NDVI in the grassland of northern Xinjiang from 2000 to 2019

在空间上，NDVImax 和TI-NDVI 的变化趋势存在明显差异，但总体上均表现出明显的海拔分异（图4a，图4b）。NDVImax 或TI-NDVI 呈上升趋势的区域主要分布在准噶尔盆地、伊犁河谷和塔城盆地的洼地。下降趋势主要分布在天山山脉、阿尔泰山山脉和塔城盆地周围山区。经计算，北疆20.87%草原的NDVImax 呈下降趋势，其中2.31%达到显著水平（Z＜-1.96）。NDVImax 呈上升趋势的草地面积占79.13%，达到显著水平（Z＞1.96）的草地面积占19.24%。TI-NDVI呈上升趋势和下降趋势的比例分别为13.80%和86.20%，达到显著水平（|Z|＞1.96）的比例分别为0.48%和30.22%（图5a）。

图4c 是NDVImax 和TI-NDVI 的变化趋势的叠加图。北疆大面积草地NDVImax 和TI-NDVI 呈相反变化趋势（图4c 中I-D 和D-I），占比达到17.55%（图5b），且多位于植被覆盖度较高的山区。尤其是阿尔泰山与伊犁河谷周边的山区，草地NDVImax和TI-NDVI变化趋势的差异性明显不同。阿尔泰山山区草地变化主要表现为NDVImax 增加而TI-NDVI减小，而伊犁河谷周边山区则主要表现为NDVImax减小而TI-NDVI 增加。整个北疆地区，2000—2019年，8.56%的区域NDVImax 和TI-NDVI 同时下降（DD），空间上主要分布在天山南部的部分地区；73.89%的区域NDVImax 和TI-NDVI 同时增加（I-I），主要分布在植被覆盖度较低的准噶尔盆地；5.24%的草地NDVImax增加而TI-NDVI减少（I-D），空间上主要分布在在阿尔泰山脉以及准噶尔盆地中部和伊犁河谷的部分区域；12.31%的草地NDVImax减少而TI-NDVI增加（D-I），空间上集中分布在伊犁河谷和塔城盆地周边山区。

图4 NDVImax（a）和TI-NDVI（b）Mann-Kendall非参数检验统计量Z值的空间分布及变化趋势叠加（c）Fig.4 Spatial distribution of the statistic Z of the Mann-Kendall non-parametric test for NDVImax(a)and TI-NDVI(b),and the superimpose d trend image of NDVImax and TI-NDVI(c)

图5 Mann-Kendall检验（a）和叠加趋势（b）的Z范围比例Fig.5 Proportions of Z ranges of Mann-Kendall test(a)and superimposed trends(b)

3.3 NDVImax 与TI-NDVI 在草地时空异质性表达上的对比

时空异质性是评价植被动态的基本问题之一。通过计算NDVImax 和TI-NDVI 的时间和空间CV，来对比两者在表达草地时空异质性的差异［41-42］。

图6a 显示了NDVImax（图6a1）和TI-NDVI（图6a2）的CVs的空间分布及其差异（图6a3）。NDVImax 和TI-NDVI 的CVs的空间格局几乎相同。CVs（CVs＞0.3）较大的区域，多位于天山和阿尔泰的高海拔区，这些区域NDVImax 和TI-NDVI 均随着海拔的升高而快速增大或减小（图2），因而CVs均较大。CVs较小的区域，主要分布在准噶尔盆地低地或天山、阿尔泰山半山腰，这些区域NDVImax 和TI-NDVI 在空间上的变化不大（图2）。然而，54.18%的区域TINDVI的CVs大于NDVImax 的CVs（图6a3），空间上主要分布在NDVImax 通常较大（NDVI＞0.6）或沿海拔梯度变化较大的山区，以及小部分植被覆盖度最低的地区（NDVImax＜0.2 和TI-NDVI＜1.0）（图6a3，图2），这说明TI-NDVI 更能表达山区植被覆盖度较高地区的草地空间异质性。相反，在准噶尔盆地地区的低地和伊犁河谷地区，NDVImax 和TI-NDVI 都很小（NDVI＜0.4 和TI-NDVI＜4.0），TI-NDVI 的CVs多小于NDVImax 的CVs（图6a），这意味着NDVImax 在表达低覆盖草地的空间异质性方面更敏感。

图6b 显示了2000—2019 年期间NDVImax（图6b1）和TI-NDVI（图6b2）CVt的空间分布及其差异（图6b3）。相较于NDVImax 和TI-NDVI 的CVs空间格局，NDVImax和TI-NDVI的CVt的空间格局差异很大（图6b1，图6b2）。对于NDVImax 而言，在准噶尔盆地、塔城盆地和伊犁河谷的洼地的CVt明显大于山区。相反，除准噶尔盆地植被覆盖度最低的部分地区外，低地和山区TI-NDVI 的CVt差异不大，均处于较低水平（0.05＜CVt＜0.1）（图6b2）。低地的植被覆盖率通常比山区低的多。植被覆盖度高的地区NDVI 容易饱和，导致高覆盖度草地NDVI 的异质性被低估，如天山和阿尔泰山区的草地。图6b3 为TINDVI 的CVt减去NDVImax 的CVt的差值图，从图上可以看出，植被覆盖度较高的山区TI-NDVI的CVt多大于NDVImax 的CVt。这意味着TI-NDVI 可以弥补NDVImax 饱和度的不足，更能表达高覆盖草地的时间异质性。然而，低地区域TI-NDVI的CVs比NDVImax的CVs要小。

图6 NDVImax和TI-NDVI时空CV的空间分布及其差异Fig.6 Spatial distribution of the spatial and temporal CV of both NDVImax and TI-NDVI,and their difference

4 讨论

4.1 NDVImax与TI-NDVI的关系

NDVImax 是植被在一个生长季所能达到的最大NDVI，而TI-NDVI定义了一个生长季植被的累积NDVI，这是NDVImax 和TI-NDVI 的根本区别。然而，北疆NDVImax 和TI-NDVI 的空间变化基本一致(图2)，图7a也显示TI-NDVI随NDVImax的增加而增加，这表明NDVImax 高的地区TI-NDVI 总体也较大。NDVImax 与TI-NDVI 的Pearson 相关系数达到0.89（P＜0.01）。因此，NDVImax 与TI-NDVI 在表示植被空间分布方面有很多共同之处。

从图7a 中可以看出，随着NDVImax 从0.1 变化到0.7，TI-NDVI离散分布在2个明显不同的范围内，这表明具有相同NDVImax 的区域通常可以分为两部分，即TI-NDVI 较小的部分和TI-NDVI 较大的部分（图7b1）。同样，TI-NDVI相同的区域也可以用同样的方式分成两部分（图7b2）。表明在相同NDVImax 或TI-NDVI 的区域可能有完全不同的TI-NDVI或NDVImax。TI-NDVI 离散分布的2 个不同范围分别对应了研究区内的高山区域和低地平原区。高海拔地区气温低、冬季长、生长季节短，如图7b1 和图7b2中的A2和B2区，因此，相较于生长季节较长的区域，即使其NDVImax可能相同或相近，TI-NDVI却要小的多（图7b1）。然而，山区具有降水丰富的优势，夏季温度升高，低温限制被解除后，NDVImax则远远高于干旱的低地荒漠区（图7b2）。

图7 NDVImax和TI-NDVI的散点图（a）及NDVImax或TI-NDVI相同但TI-NDVI或NDVImax不同的区域的NDVI（b）Fig.7 Scatter plot of NDVImax and TI-NDVI(a)and NDVI curves of pixels with same NDVImax but different TI-NDVI,NDVI curves of pixels with same TI-NDVI but different NDVImax(b)

NDVImax 和TI-NDVI 均呈现明显的海拔分异。水和热是控制植被空间分布的基本因素。然而水和热沿海拔梯度变化。北疆地处亚欧大陆干旱区腹地，准噶尔盆地、伊犁河谷和塔城盆地等海拔较低的区域植被受缺水和高温制约，所以这些区域的NDVImax 和TI-NDVI 均很小。随海拔的升高，降水增加，气温降低，对植被生长的限制逐渐减缓，NDVImax 和TI-NDVI 均逐步增大。当海拔达到一定高度且水和热均最为适宜时，NDVImax 和TI-NDVI 达到最大，如天山和阿尔泰山的半山腰等地区。随着海拔的升高，植被受低温和饱和土壤水的限制，NDVImax 和TI-NDVI 又开始下降，如天山和阿尔泰山的高海拔地区。海拔对水热条件的再分配是草地NDVImax和TI-NDVI空间分异的根本原因，这也是北疆草地NDVImax和TI-NDVI均呈海拔分异的原因。

4.2 NDVImax与TI-NDVI变化差异的解释

NDVImax 和TI-NDVI 在年际波动（图3）和变化趋势的空间格局（图4c）存在差异。用生长季内不同时间的NDVI 绘制的曲线不仅显示了NDVImax（NDVI 曲线的顶点），而且还显示了TI-NDVI（曲线上累积的NDVI）的大小（图8）。因此，通过比较不同生长季或不同地区的NDVI曲线，可以揭示NDVImax与TI-NDVI变化差异的原因。

对于整个研究区域，2000—2001 年、2005—2006 年、2014—2015 年和2018—2019 年的区域平均NDVImax 和TI-NDVI呈现相反的波动（图3）。从图8a1 可以看出，2001 年NDVI 曲线的峰值较低，但2001 年第209～305 d 的NDVI 比2000 年大得多，表明2001 年仲夏到秋季的草地状况较好。这是2000—2001 年NDVImax 下降而TI-NDVI 增加的原因。从图8a2、图8a3和图8a4中可以看出，2005年、2015年和2019年NDVI曲线的峰值均大于2004年、2014 年和2018 年，但春秋季的NDVIs 均较小。2004—2005 年、2014—2015 年 和2018—2019 年，NDVImax增加，TI-NDVI减少。

图8 TI-NDVI与NDVImax变化趋势相反年份及两者变化趋势不同叠加类型在研究期前期和末期NDVI生长季内的变化Fig.8 Changes of TI-NDVI and NDVImax in the years with opposite trends and different superimposed types of the two trends in the NDVI growing season at the beginning and end of the study period

对于NDVImax 和TI-NDVI 变化趋势的空间差异。天山山脉和塔城盆地周边的大部分地区NDVImax 降低而TI-NDVI 增加。相反，阿尔泰山大部分地区NDVImax 增加而TI-NDVI 降低（图4c）。对比研究期前5 a（2000—2005 年）与后5 a（2015—2019年）的平均曲线（图b）可知，在天山及塔城盆地周边山区，2000—2004 年NDVI 曲线的峰值比2015—2019年略大，但2015—2019年曲线中第65～145 d的所有NDVI 均大于2000—2004 年（图8b1）的值，这表明尽管NDVI 峰值有所下降，但春季NDVIs 则有所改善，提高了这些地区的TI-NDVI。对于阿尔泰山地区，夏季NDVI虽有增加，但春季和秋季减少更多（图8b2），导致NDVImax 增加而TI-NDVI 减少。图8b3 为准噶尔盆地广阔洼地的平均曲线，其中NDVImax 和TI-NDVI 在2000—2019 年间均呈现增加趋势（图4c）。可以看出，整个生长季的NDVI 均有所提高（图8b3）。

单一生长季的NDVI曲线代表了草地的年生长过程。这种一年一度的生长过程通常在时间和空间上都有所不同。而年生长过程的年际变化也是植被动态的重要特征之一，但这一点在过去的草地动态监测中则很少被考虑。与NDVImax 相比，TINDVI 积累了草地整个生长季的NDVI，这使得TINDVI 成为在表征草地年生长过程方面具有价值的量化指标。

5 结论

（1）北疆地区NDVImax和TI-NDVI的空间变化基本相同，均表现出明显的海拔分异。TI-NDVI 总体随NDVImax的增大而增大，但NDVImax或TI-NDVI相同的区域，其TI-NDVI或NDVImax 却存在较大差异。

（2）2000—2019 年北疆地区草地TI-NDVI 和NDVImax 的总体呈现增加趋势（P＜0.01），但空间上，整个北疆17.55%草地的TI-NDVI和NDVImax呈现相反的变化趋势，且多位于草地覆盖度较高的山区。尤其是阿尔泰山与伊犁河谷山区，阿尔泰山草地主要表现为NDVImax 增加而TI-NDVI 减小，伊犁河谷山区草地则主要表现为NDVImax 减小而TINDVI增加。

（3）在时空异质性上，北疆阿尔泰山、伊犁河谷及塔城盆地周边山区等草地高覆盖区，TI-NDVI 的CVs和CVt均高于NDVImax，TI-NDVI 能弥补草地动态检测中NDVI 光饱和缺陷的影响，能更准确地表达高覆盖草地的时间异质性。