基于连续时间序列NDVI数据的中国生态状况时空变化特征

宗加权, 白淑英, 冯朝阳, 蒲 阳

(1.南京信息工程大学, 南京 210044 ; 2.中国环境科学研究院,国家环境保护区域生态过程与功能评估重点实验室, 北京 100012)

植被是陆地生态系统中的重要组成部分，在陆地表层能量交换、水分循环和生物地球化学循环过程中扮演着重要角色[1]，植被覆盖度变化是生态环境变化的直接结果，很大程度上表示了生态环境总体状况[2]。相对于传统实地调查植被覆盖特征费时费力，通过遥感手段获取数据，优势有空间范围覆盖面积大、时间序列长、方便获取等。归一化差值植被指数(NDVI)是目前常用的研究地表植被覆盖特征的重要指标之一[3]，长时间序列监测植被覆盖特征变化对更好地模拟陆地生态系统的动态变化特征具有重要意义[4]。AVHRR NDVI是目前时间序列最长的全球连续数据集，具有很高的科研价值和实际意义[5-6];MODIS NDVI是从2000年开始的数据集，提高了叶绿素敏感度和空间分辨率、排除了大气水汽的干扰、改善了合成方法，但时间序列长度有限。探寻两者间的关系，通过数据拟合构建长时间序列的NDVI数据集，是实现对植被长期观测的关键[7-8]。一些学者基于不同NDVI数据集重构NDVI数据，并结合气象数据等分析了植被状况的时空变化特征，如:李净等[9]基于重构后的数据产品和植被类型数据集，分析了不同时期植被NDVI变化特征，发现西北地区NDVI整体呈上升趋势;李飞等[10]根据MODIS NDVI与AVHRR数据构建了较长时间序列NDVI数据集，采用相关性分析法分析了全国农业、森林、草地与稀疏植被区植被活动的变化特征，结果指出近30 a来中国陆地植被活动整体增强，这与方精云等[11]利用AVHRR NDVI研究中国地区植被活动的结果相同。赵玉萍等[12]利用1982—2003年GIMMS NDVI数据集和逐月气象资料，用相关分析法研究了藏北高原草地生态系统NDVI与气候因子的相关性及滞后性，得到NDVI与月平均气温等呈高度正相关(p<0.001)。刘可等[13]利用1982—2012年GIMMS NDVI3g和中国陆地生态系统类型数据，采用一元线性经验模态分解和相关分析方法，研究了近30 a中国各生态系统NDVI的时空变化特征，结果表明近30 a中国植被活动显著上升，农田、森林、草地和水体与湿地生态系统的NDVI总体为非稳定上升趋势。中国陆地生态系统中的森林、草原、荒漠、水体与湿地和其他生态系统类型是对土地利用/覆盖的分类判别汇总，对生态系统状况有指示意义。水体与湿地在生态环境评价中起到重要作用，其水体与湿地状况好坏和面积变化影响着生态环境状况，水域面积的减少会导致生态环境恶化[14]。水体和湿地因为水源丰富、土壤肥沃，很容易被开垦为良田和进行水产的养殖，水体与湿地生态系统因此被严重破坏[15]，牛振国等[16]研究表明，近30 a，中国湿地类型由湿地转化为非湿地主要是气候变化和农业活动引起的。

以上研究存在一些问题，首先采用的AVHRR数据为一代数据产品，数据时间跨度不长，数据的空间分辨率不高(8 km×8 km);其次大多数研究多从行政区划和地理单元探究植被活动特征，未考虑到不同生态系统间的镶嵌性，因此难以详细地探究不同生态系统间NDVI的差异。对此，本研究基于重构的1981—2018年长时间序列植被指数数据集(GIMMS NDVI3g和MODSI NDVI)，采用分区统计、NDVI趋势分析法，分析中国大陆范围(包含香港、澳门，未包含中国台湾省)非耕地、非聚落和湿地生态系统NDVI的空间分布和时间动态变化特征，试图为中国生态环境保护和可持续发展提供决策依据。

1 数据来源及预处理

1.1 数据来源

本研究利用2000年2月—2018年12月MODIS13 A3数据集植被指数产品的月合成数据(https:∥ladsweb.modaps.eosdis.nasa.gov/search/)，空间分辨率为1 km。还有美国国家航天航空局NASA Goddard Space Flight Center的全球监测与模型研究组(GIMMS)采用最大值合成法(MVC)合成的15 d、空间分辨率为8 km的最新GIMMS NDVI 3g.v1数据，时间范围为1981年7月—2015年12月(https:∥ecocast.arc.nasa.gov/data/pub/gimms/3g.v1)。

中国陆地生态系统数据来源于中国科学院地理科学与资源研究所、中国科学院资源环境科学数据中心的“资源环境数据云平台”。包括1980年、1990年、1995年、2000年、2005年、2010年和2015年，共7时期数据。该数据是在遥感解译获取的1∶10万比例尺土地利用/土地覆盖数据的基础上，通过对各生态系统类型进行辨识和研究，经过分类处理形成多期中国陆地生态系统类型空间分布数据集，矢量数据栅格化生成的1 km栅格数据。生态系统类型依据土地利用/土地覆盖遥感分类系统，归分为7大生态系统类型:(1) 农田生态系统，主要包括水田和旱地;(2) 森林生态系统，主要包括密林地(有林地)、灌丛、疏林地、其他林地;(3) 草地生态系统，主要包括高覆盖度草地、中覆盖度草地、低覆盖度草地;(4) 水体与湿地生态系统，主要包括沼泽地、河渠、湖泊、水库、冰川与永久积雪、滩地;(5) 聚落生态系统，主要包括城镇、农村居民地、工矿;(6) 荒漠生态系统，主要包括沙地、戈壁、盐碱地;高寒荒漠;(7) 其他生态系统，主要包括裸土地和裸岩砾石地。

1.2 数据预处理

采用MRT软件对MODIS NDVI数据进行拼接，投影转换，数据格式转换。采用Matlab软件对GIMMS NDVI数据进行逐月最大值合成(MVC)，进一步减少云的影响，并降低月内物候循环的影响[17]，还原真值，格式转换。由于GIMMS NDVI数据的空间分辨率约为8 km，将GIMMS NDVI由8 km重采样为1 km的数据，不损失影像信息，以匹配MODIS NDVI。本文采用的地理坐标系是WGS-84,Albers投影坐标系，随后统一采用IDL程序对两套数据进行逐月裁切、掩膜、还原真值等操作。

2 研究方法

MODIS NDVI数据被称为GIMMS NDVI数据的延续[18]，两套数据集有一段时间的重叠部分，即2000年2月份-2015年12月份。已有学者[17,19]的研究表明可以利用两套数据的重叠时间部分，以此建立逐月线性转换方程，但未能从文章中找出适合本研究区域的转换方程，本研究中建立的最佳转换方程为y=0.970034x+0.016953，拟合方程的判定系数R2为0.88,p<0.01，拟合效果较好。利用转换方程将1981年7月—2001年12月的GIMMS NDVI进行修正拟合，为了保留MODIS数据1 km的空间分辨率，将GIMMS NDVI数据也重采样为1 km，并连接MODIS NDVI 2002年1月至2018年12月数据，最终获得中国大陆地区的1981年7月至2018年12月长时间序列空间分辨率为1 km的NDVI数据集。

2.1 NDVI时空变化趋势分析法

2.1.1 NDVI时间变化趋势分析法 趋势分析是对一组随时间变化的变量进行回归分析，得到其变化趋势的方法。对多年植被指数数据逐像元NDVI值进行一元线性回归，得到研究区内植被随时间的变化趋势，即NDVI的年际变化。计算公式为[20-21]:

(1)

式中:k为研究时间段内的年数,取1-n；YNDVIk是第k年某季相的NDVI平均值,如夏季NDVI平均值采用当年6月、7月、8月份NDVI值求得平均；slope>0表示植被覆盖呈增加趋势;slope小于0则相反;slope=0表示植被覆盖没有明显的变化。

2.1.2 NDVI空间变化分析法 将大陆地区分为六大区域(东北、华北、华东、中南、西北、西南)，分区统计分析NDVI的总体变化特征，能整体上体现NDVI均值和趋势变化情况。

2.2 主要生态系统类型的NDVI时空变化分析法

采用分区统计的方法来统计不同生态系统NDVI的变化情况，及其空间分异特征。由于生态系统数据为多年一幅的栅格数据，为充分利用生态系统数据，将连续时间序列的NDVI与最邻近年的生态系统数据进行对应统计，表1是分区统计对应表。

表1 分区统计年份对应

3 结果与分析

3.1 中国植被NDVI时空变化分析

3.1.1 中国植被NDVI时间变化分析 图1是整个中国大陆范围内NDVI大于0值，2000年2月至2015年12月GIMMS NDVI平均值修正前后与MODIS NDVI平均值变化曲线，由图1可知，原始GIMMS NDVI的周期性变化趋势与MODIS NDVI相似，每年冬季原始GIMMS NDVI都小于MODIS NDVI，在植被生长季，两者相差不大，修正后GIMMS NDVI与MODIS NDVI值相接近，具有一致性。值得注意的是MODIS NDVI数据在2000年初和2001年冬季，其值明显低于修正前后GIMMS NDVI数据，可见MODIS数据在刚获取的阶段，数据质量可能不稳定，所以对于中国大陆地区，在使用MODIS NDVI数据时，可以从2001年7月份开始选用。

图1 2000年2月至2015年12月MODIS NDVI、原始与修正后GIMMS NDVI时间序列曲线对比

由表2看出，修正后GIMMS NDVI在3种数据中的平均值和最小值都为最大，平均值最大为0.373，最小值最大为0.273;除了2000年和2001年，修正后GIMMS NDVI与MODIS NDVI的年平均值环比增长率和最大值有相似的变化规律，且有相同的NDVI线性趋势(0.000 8/a)，最大值最小仅相差0.002。这说明，除了2000年、2001年，修正后GIMMS NDVI在年尺度上与MODIS NDVI相似的变化规律。

表2 不同数据年平均值、最大值、最小值和年平均值环比增长率

从表3中看出，MODIS NDVI和修正后GIMMS NDVI在2月至5月NDVI平均值环比增长率逐渐增加，而从5月至6月NDVI值月平均环比增长率有所下降，这是因为冬小麦成熟收割，新播种的作物未生长成熟造成的NDVI值下降现象，原始GIMMS NDVI未能捕捉到这个现象;修正后GIMMS NDVI平均值和最大值与MODIS NDVI有相似的变化规律，GIMMS NDVI月均值环比增长率最大为0.266%,MODIS NDVI为0.273%;修正后GIMMS NDVI最小值大于MODIS NDVI,2月份最大差值为0.096，这可能是由于2000年2月份MODIS NDVI数据刚刚获取，数据质量不稳定导致NDVI值偏低，10月份相差值最小为0.008。

表3 不同数据月平均值、最大值、最小值和月平均值环比增长率

3.1.2 中国植被NDVI空间变化分析 从图2中可以看出中国大陆夏季的NDVI的空间分布整体呈现东南和东北高、西北低的格局。在沿海的华北平原、长江三角洲和珠江三角洲等经济发达地区呈现低NDVI值的现象，这是由于人类自身发展的需要，NDVI值随着城镇化的发展而降低;在新疆的阿尔泰山南坡、天山山脉北坡的NDVI值却反常高，最高可达到0.896，是因为在准格尔盆地西部有缺口，大西洋和北冰洋的水汽能够影响到两侧的山脉，形成迎风坡，对于植被生长有利;另外，在西北新疆塔里木盆地、准格尔盆地、甘肃西北部、青海西北部、内蒙古西部有沙漠、荒漠戈壁等，西藏自治区内的高原有高原荒漠和裸岩，气候恶劣，不利于植物生长的环境，所以NDVI值普遍偏低，平均值为0.275。华北平原主要有农田，冬小麦收割之后，夏季播种农作物还未完全长成，掺杂土壤背景信息，因此NDVI值较低，平均值为0.471。位于华东地区的浙江、福建有丘陵，贵州、重庆、湖南和湖北西部有山地适宜森林生长，NDVI值较高，平均值为0.668。内蒙古东北部和东北三省内有大兴安岭、小兴安岭和长白山山脉，夏季森林和农作物生长旺盛NDVI值较高，平均值为0.775。

图2 2018年夏季中国大陆NDVI平均值

表4 2018年夏季中国大陆六大分区NDVI均值

由图3和表5，6可以看出过去38 a中，有49.923%的NDVI稳定不变，49.846%的NDVI轻微改善，中国大陆地区整体植被没有发生剧烈变化。西北地区NDVI趋势整体无明显变化(0.000 098/a)，但新疆天山山脉北坡和塔里木河流域处发现NDVI的变化较为剧烈，最大为0.030 272/a，呈上升趋势;内蒙古阴山南、甘肃省内祁连山、宁夏贺兰山等NDVI变化为上升趋势;在经济发达的长江三角洲、珠江三角洲、京津翼地区NDVI趋势呈现下降现象。在西藏藏南地区、内蒙古大兴安岭、东北三省内的小兴安岭和长白山，夏季NDVI值较高，而年际变化趋势中呈现下降趋势，说明植被覆盖浓郁，多年以来的植被覆盖有所下降。

表5 1981-2018年六大分区NDVI趋势均值

表6 1981-2018年NDVI变化趋势分类级别

图3 1981-2018年中国大陆夏季NDVI变化趋势

3.2 主要生态系统类型的NDVI变化分析

中国陆地生态系统具体分为农业、森林、草地、水体与湿地、聚落、荒漠和其他共7种生态系统类型。农田和聚落受人类活动影响较大，近30 a来农作物总产和单产都有所提高[22]，对于全国范围内植被NDVI增加是正向作用，但想要更详细地分析中国大陆的生态状况则需要中国陆地生态系统分类数据的支持。自然生态系统有森林、草地、水体与湿地、荒漠或者其他生态系统，这些生态系统类型受人类影响较小，对生态系统状况有指示意义，按表1对上述生态系统分区统计其NDVI长时间序列变化，有助于中国生物多样性保护和开展生态环境状况评价。

森林和农田生态系统在过去38 a中NDVI呈现波动式上升，上升趋势分别为0.001 4/a，0.001 1/a，森林生态系统NDVI最低值为1982年0.692，最高值为2017年0.766。1987年5月6日至6月2日大兴安岭发生特大森林火灾;2000年6月17，18日黑龙江大兴安岭呼中区呼源、松岭区发生特大森林火灾，同年8月6日新疆裕民县巴旦杏林保护区发生重大森林火灾;2002年7月28日和2004年6月22日内蒙古大兴安岭北部原始林区分别发生特大、重大森林大火;2018年6月1日和2日内蒙古大兴安岭汗马国家级自然保护区和北部原始林区发生森林火灾等，这些森林大火使NDVI值在时间序列上呈现出偏低或者谷。在过去38 a中草地生态系统的NDVI在0.358至0.402间变化，草地生态系统NDVI变化较为剧烈，说明NDVI对于不同覆盖度草地的变化敏感，1995年、1997年锡林浩特市草原分别发生3次和1次特大火灾，该年份NDVI值有所下降，草地生态环境状况变差。荒漠生态系统的NDVI在过去38 a中呈下降趋势(-0.000 6/a)，最低为0.094，最高为0.126，荒漠生态系统主要包含沙地、戈壁、盐碱地、高寒荒漠，这些土地类型环境恶劣，十分脆弱，NDVI能够反映出其变化状况，在1981—2001年NDVI几乎没有大的变化，处于相对稳定的状态，而在2002—2018年NDVI值的上升趋势为0.000 9/a，说明近17 a的生态环境呈现向好的态势发展。其他生态系统的NDVI在过去38 a中呈现下降趋势(-0.000 7/a)，最低为0.138，最高为0.184，与荒漠生态系统类似，在1981—2001年，NDVI也变化不大，2002—2018年，其NDVI的上升趋势与荒漠生态系统一样，生态系统呈现向好态势。水体与湿地生态系统NDVI在过去38 a中呈现下降趋势(-0.001 7/a)，最低为0.311，最高为0.401，在1981—2001年NDVI呈现较稳定状态，在2002—2018年中呈现上升趋势(0.001 6/a)。

图4 1981-2018年森林和农田生态系统NDVI变化

4 讨 论

本研究利用NDVI探究生态系统状况，但可能有所欠缺，另如气候、地形等多种因素可以加以综合考虑;趋势线分析方法能很好地体现NDVI的变化趋势，但是如何寻找合适的NDVI变化趋势分级区间，需要进一步研究。刘可等[13]仅使用2010年一时期的中国陆地生态系统类型数据统计近30 a GIMMS NDVI3g数据，指出森林、农田、草地和水体与湿地总体NDVI非稳定上升，本文利用生态系统类型数据统计重构后的NDVI数据，指出森林和农田生态系统NDVI整体呈增强趋势，草地生态系统NDVI呈非稳定上升趋势，这可能是因为本文利用的是多期生态系统类型数据。利用生态系统分区统计不同年份的NDVI数据时，虽然尽可能利用生态系统数据，但最终统计的结果可能与实际的有差别，荒漠、水体与湿地和其他生态系统的NDVI曲线都在2000—2001年处有下降的陡坡，这可能是因为对GIMMS NDVI数据进行重采样，仍无法解决其分辨率粗的问题，但若将MODIS NDVI数据进行降分辨率，则会损失大量信息。

5 结 论

(1) 构建数据集时，修正后GIMMS NDVI和MODIS NDVI相似的年份和月份变化规律;中国大陆地区MODIS NDVI在2000年、2001年的冬季数据质量不佳，可以从2001年7月开始选用数据。中国大陆地区夏季NDVI整体分布特征呈现西北低，东南和东北高，在经济发达区域，由于人类经济社会活动因素出现低NDVI值现象。

(2) 在过去38 a中，对受人类活动影响较小的生态系统，森林生态系统植被呈上升趋势，森林火灾对森林的影响在NDVI序列变化能够反映出来;草地生态系统呈现稳定趋势，NDVI变化较为剧烈，NDVI序列可以作为监测草地生态状况变化的指示;荒漠、水体与湿地、其他生态系统的都呈下降趋势，三者在1981—2001年NDVI变化都较为稳定，且2002年以后NDVI呈现向好态势发展，这与MODIS传感器分辨率和波段宽度设置优于AVHRR传感器有关。