基于NDVI时间序列的近12年中国植被覆盖的单调趋势分析

李 强，张 翀

（1.陕西师范大学 旅游与环境学院，陕西 西安710062；2.陕西学前师范学院 环境与资源管理系，陕西 西安710100）

植被作为陆地生态圈的重要组成，是气候系统的重要元素［1］。众所周知，全球气候正处于一个持续变暖的阶段，强烈影响陆地生物圈［2］，在这样的背景下，掌握陆地植被覆盖度年际间的变化规律，对评价陆地生态系统的环境质量、调节生态过程具有重要的理论和实际意义［3］。归一化植被指数（normalized differ－ence vegetation index，NDVI）由红波段与近红外波段的反射率计算而来［4］，和植物的生产力密切相关［5］，并且NDVI趋势可以用来衡量植被覆盖的改善与退化［6］。NDVI趋势被用在很多测算中，包括全球变暖的生态响应［7］、物候变化［8］、作物状况［9］、土地覆盖变化［10］以及沙漠化［11］。植被覆盖的趋势及年际变化会影响植被与大气之间能量的交换［12］。一系列研究表明，在全球，特别是北半球，存在生长季的开始时刻提前以及生长季长度增大的趋势［13］，生长季长度增大与气温的上升会加快地表水分蒸发，增大干旱胁迫与林火发生的可能［14］，并增加了固碳强度［15］。因此，生长季内植被覆盖相对于全年更能作为反映诸如土壤退化等的指示器。利用NDVI序列进行趋势分析时，由于数据集存在自相关会影响模型假设，线性模型应慎重使用。所以既要剔除序列的季节性，也需要进行非参数趋势检验［16］。本文基于时间序列谐波分析法、线性趋势和Kendall’sτ趋势等方法对中国植被生长季起止时刻进行了确定，并计算得到中国生长季植被覆盖的单调性趋势以及生长季长度趋势，并进一步解释不同植被类型植被覆盖趋势的差异以及生长季变化规律，以期为中国当前的生态建设、修复提供有用的空间信息和理论支撑。

1 数据与方法

1.1 数据来源

本文所用的资料包括黄土高原1999—2010年SPOT VEGETATION旬值NDVI数据、131个台站12a（1999—2010年）的年降水资料和年平均气温资料。NDVI数据来自于互联网（http：∥free.vgt.vito.be／home.php），空间分辨率为1 000m。植被类型数据来自寒区旱区科学数据中心。本文对旬值NDVI数据经过时间序列谐波分析重构其旬值序列，并得到旬值异常值序列，进而对逐年异常值序列求均值得到年际异常值序列；利用NDVI旬值重构序列提取植被物候特征，确定植被生长季始末期，从而计算得到逐年生长季内NDVI的均值时间序列。

1.2 NDVI时间序列的谐波分析

时间序列谐波分析法（harmonic analysis of time series，HANTS）是平滑和滤波两种方法的综合，它能够充分利用遥感图像存在时间性和空间性的特点，将其空间上的分布规律和时间上的变化规律联系起来。时间序列谐波分解法进行影像重构时充分考虑了植被生长周期性和数据本身的双重特点，能够用代表不同生长周期的植被频率曲线重新构建时序NDVI影像，真实反映植被的周期性变化规律。剔除季节影响几乎可以完全消除数据的季节性特征。剔除季节影响后，一期与另一期的对比将更有意义，而且可以帮助确定是否存在趋势。HANTS分析中单独年份的谐波分析与整个时间段的谐波分析之间的差异即为异常值序列，异常值序列剔除了数据的季节性特征［17］。

1.3 物候特征提取

多种处理方法被用来描述NDVI时间序列的生长季的起始时刻（start of growing season，SOS），如半极大值［18］、10% 阈 值 法［19］、拐 点 法［20］、最 大 曲 率法［21］、滑动 平 均 法［22］。 本 文 按 照 Whitetffu 等［8］的方法，采用HANTS平滑NDVI的一阶导数，一阶导数的最大值对应时刻为SOS（最大NDVI增长），生长季的结束时刻（end of growing season，EOS）为SOS过后首次NDVI下降到与SOS相等的时刻。该方法与其他方法比较相对可靠。本文主要采用各年NDVI年内变化的三阶导数，但是这种方法只适合于一年一熟作物区［17］，所以对于农耕区的一年多熟区域采用10%阈值法进行计算［19］。求得中国区域范围内每个栅格的每年的SOS（一阶导数的极大值）与EOS（一阶导数的极小值），SOS至EOS即为生长季，两者之间的时间长度即为生长季长度（length of growing season，LOS）。

1.4 趋势分析

1.4.1 线性趋势 运用线性趋势线分析植被覆盖的变化趋势，即以时间为自变量，NDVI为因变最，利用最小二乘法计算斜率值K。K值的符号反映上升或下降的变化趋势，K＜0表示在计算时段内呈下降趋势，K＞0表示呈上升趋势。K值绝对值的大小可以度量其演变趋势上升、下降的程度［23］。

1.4.2 Kendall’sτ趋势 Mann等［24］首次建议用Mann—Kendall法来检验时间趋势Kendall’sτ的显著性，该方法已经被应用在季节数据上，主要是水文分析方面，近年应用在 NDVI数据上［6，25］。为了更好显示出Kendall’sτ的空间分布，本文将其分为8个等级：极显著持续退化（－1～－0.75）、较显著持续退化（－0.75～－0.5）、显著持续退化（－0.5～－0.25）、不显著持续退化（－0.25～0）、不显著持续改善（0～0.25）、显著持续改善（0.25～0.5）、较显著持续改善（0.5～0.75）、极显著持续改善（0.75～1）。

2 结果与分析

2.1 剔除季节性的NDVI趋势

图1为中国范围内内剔除季节性NDVI数据年际变化曲线，植被覆盖呈明显增加趋势（0.03／10a），低于中国 1998—2007年 植 被 覆 盖 速 度 0.048／10a［26］。（1）1999—2000年中国植被覆盖减小，2000—2004年植被覆盖增加，2000年植被覆盖最低，2001年开始上升，主要是由于2000和2001年属特大干旱年，且2000年干旱程度较2001年严重，该时间中国大范围降水 偏 少，发 生 了 建 国 以 来 最 严 重 的 旱 灾［27－28］；（2）2004年植被覆盖达到峰值，随之在2005年下降，2005年以后植被覆盖呈稳定上升趋势。

图1 中国1999－2010年剔除季节性NDVI数据的年际变化

1999—2010年剔除季节性旬NDVI数据小波变化系数如图2所示，根据其规则交替，显示出在80～110频率尺度上存在5～7a（180～252旬）的周期，在1～30频率尺度上存在1～2a（36～72旬）的周期。图1也显示出中国植被覆盖变化呈现出一定的周期性特点：4a的稳定增长期（2000—2004与2005—2009年）、1a的突然下降期（1999—2000，2004—2005与2009—2010年），即存在5a的周期。

图2 中国1999－2010年剔除季节性NDVI旬数据的小波变换系数

利用线性趋势原理，得到中国剔除季节性NDVI的线性趋势空间分布（如附图8所示）。植被明显改善区集中在黄土高原—坝上高原以南，横断山脉—黄土高原西缘以东的大面积区域，另外新疆绿洲区也呈改善趋势；明显退化区，在中国东部及东南部，主要分布于黄河三角洲平原、长江三角洲平原、珠江三角洲平原与腾格里沙漠以南—陇中高原以北的区域，其中黄三角植被结构简单，生态系统年轻化的特点，并且是国内外少有的资源富集区，开发潜力大，脆弱的生态环境在近12a来经济快速发展中呈退化趋势，而长江三角与珠江三角以及腾格里沙漠—陇中高原之间的区域，主要是由于快速经济发展、开发利用以及城市化等要素导致植被覆盖的下降，朴世龙等［29］发现长三角与珠三角在过去18a间也呈下降趋势。另外在乌兰察布高原—浑善达克沙地南缘—锡林郭勒高原东南部一线的条带地区退化较为明显，与穆少杰等［3］研究结果一致。

图3反映出，中国各种植被类型均呈改善趋势。（1）常绿针叶林、常绿阔叶林、稀疏灌丛、低山草原、耕地为明显改善区，趋势值在0.004～0.007，年际波动性最强。常绿针叶林与常绿阔叶林主要分布在中国华东南部、华南和西南地区的山地丘陵区，尤以横断山脉、大巴山改善明显；稀疏灌丛分布在两广南部的地势平坦区；低山草原主要分布在黄土高原及其周边低山地区；耕地分布较为广泛，但陕北、河南—安徽的黄淮平原的植被改善趋势明显高于其他地区，其中陕北自1999年退耕还林还草开始以来，林草替换原有耕地，植被覆盖明显增加［30］，黄淮平原近50a来气温变化不明显，降水呈增加趋势［31］，降水量的增加促使该地区植被覆盖的增加。（2）落叶针叶林、落叶阔叶林、密集灌丛、高山亚高山牧场草地、典型草原、牧场草地和沼泽湿地的改善趋势次之，趋势值在0.001～0.004，波动性相对较弱。落叶针叶林集中分布在大小兴安岭以及新疆的阿尔泰山与天山山脉；落叶阔叶林分布于秦岭北麓、小兴安岭以及长白山脉与坝上高原；密集灌丛分布于青藏高原东部、云贵高原、罗霄山等山地区。（3）海岸湿地、荒漠草原与温带高山—亚高山高寒草原改善趋势不明显，趋势值在0～0.001，波动性最小，为生态脆弱区。

综上所述，对于趋势及波动性：常绿林＞落叶林、稀疏灌丛＞密集灌丛、低山草原＞高山亚高山牧场草地、典型草原与牧场草地＞荒漠草原与温带高山—亚高山高寒草原、耕地＞沼泽湿地＞海岸湿地。反映出中国植被覆盖受纬度地带性影响较为明显。纬度较低，生物多样性程度较高，则植被覆盖改善趋势明显，相反则生物多样性程度下降，改善趋势次之。轻微改善的地区多为生态环境脆弱区，其生物多样性很低，易受人类活动等因素的影响，趋势值处于0值附近。

2.2 生长季NDVI的单调趋势

附图8b为中国1999—2010年生长季植被覆盖的Kendall’sτ趋势空间分布。Kendall’sτ趋势平均值为0.2873，持续改善的面积比重为85.85%，持续退化的比重占14.15%，说明中国生长季植被覆盖总体上呈持续改善趋势。

空间分布上（附图8b）：（1）中国生长季植被覆盖显著、较显著与极显著持续退化区比重为3.91%，主要分布于长江三角、珠江三角、腾格里沙漠—陇中高原之间的区域、柴达木盆地东南部与乌兰察布高原—浑善达克沙地南缘—锡林郭勒高原东南部一线的条带地区［3］，植被覆盖退化区还零散分布在青藏高原整个地区，尤以南部为主的区域。（2）不显著持续性变化区为易变化区域，比重为38.25%，集中分布与新疆—内蒙古的沙漠区，该地区NDVI小于0.1，变化可能与到达传感器的能量的微小差异引起，这种微小差异与传感器本身以及大气的影响有关；四川盆地、东南丘陵以及退化区周边地区为不显著持续性变化区。（3）持续改善区分布范围遍及全国，持续性较高的区域（极显著与较显著区域）主要分布在105°E以东，30°—40°N 的区域，比重为23.52%，该区域北界大致为中温带与暖温带过渡处，南界大致为中亚热带与北亚热带的过渡处，西界至高原温带边缘；另外在河西走廊、横断山区、新疆绿洲区以及青藏高原东北部也有分布。

由图3b发现，各种植被类型的Kendall’sτ趋势与线性趋势规律具有很强的相似性，即纬度地带性明显，低纬度地区植被多样性程度高，持续性强，反之持续性弱，从而也可以表明趋势的持续性与趋势的大小存在很强的相关性，即趋势高则持续性大，反之持续性弱。

图3 中国不同植被类型的统计值注：1.落叶针叶林；2.常绿针叶林；3.常绿阔叶林；4.落叶阔叶林；5.密集灌丛；6.稀疏灌丛；7.海岸湿地；8.高山—亚高山牧场草地；9.低山草原；10.典型草原；11.荒漠草原；12.牧场草地；13：沼泽湿地；14.耕地；15.温带高山—亚高山高寒草原

2.3 生长季长度趋势

受物候变化影响的生长季Kendall’sτ趋势由SOS与LOS的变化决定。如果植物生长季稳定，则年际SOS与EOS基本不变。附图8c表明LOS的线性趋势与其变异系数的空间变化，中国植被生长季长度平均趋势为0.763 6，总体上呈增加趋势。

附图8c—d显示出，LOS趋势值主要集中在－1～3.5，减小与增加趋势的比重相当，生长季长度减小的区域比重占49.88%。（1）连片减小区域主要分布在青海东部及南部与天山山脉，波动性较小（图3d）；陕北高原—吕梁山—太行山以及坝上高原大兴安岭两侧、小兴安岭、长白山周围的低山区，波动性很小（图3d）；内蒙古的呼伦贝尔高原波动性较小（图3d），以及山东丘陵，主要分布着耕地，波动性较高（图3d），其生长季长度也呈减小趋势。（2）生长季长度增加的比重为50.12%，连片增加区域主要分布在高度较高的山地区，包括阿尔泰山脉、秦巴山区、大小兴安岭与长白山脉，该地区除常绿林以外的其他植被类型的生态环境较为脆弱，波动性很高，值大于1.5（图3d）。

SOS与EOS的年际变化以及线性趋势显示，SOS呈减小趋势，而EOS成增加趋势，并且EOS趋势的增大程度大于SOS减小程度（0.589 6＞0.021 1），SOS与EOS的年际变化特征基本相反，峰值对应谷值，即SOS提前—EOS推迟或SOS推迟—EOS提前，并且后者情况只有在2007年表现明显，所以中国LOS总体上的增长趋势由SOS提前和EOS推迟造成，主要是因为在气候增温［30］的趋势下，春季气温会提前达到植被生长所需的适宜温度，而在秋季，气温的下降趋势得到延缓。

3 结论

（1）中国植被覆盖呈明显增加趋势，各种植被类型年际变化规律基本相同。中国植被覆盖变化呈现出一定的周期性特点：4a的稳定增长期，1a的突然下降期。

（2）中国生长季植被覆盖Kendall’sτ趋势平均值为0.287 3，总体上呈持续改善趋势，持续改善的面积比重为85.85%，持续退化的比重占14.15%，持续性较高的改善区主要分布在105°E以东，30°—40°N之间的区域，比重为23.52%，该区域北界大致为中温带与暖温带过渡处，南界大致为中亚热带与北亚热带的过渡处，西界至高原温带边缘。

（3）中国植被覆盖受纬度地带性影响较为明显，纬度低，生物多样性程度较高，则植被覆盖改善趋势明显，持续性强，反之改善趋势及持续性较低；黄三角和风沙区等生态脆弱区，长江三角、珠江三角与兰州市周边地区等经济快速发展区，植被覆盖退化严重。

（4）中国植被生长季长度平均趋势为0.763 6，总体上呈增加趋势，主要由EOS的增大趋势引起。本文之所以采用2000年土地覆盖数据，就是因为退耕还林从1999年开始，而1999—2010年耕地的植被覆盖增大趋势明显，而且呈持续性相对较高，这意味着在此期间有很大面积的耕地转化为了林地［30］；从1978年“三北防护林”建设以来，西北、华北北部、东北西部植被覆盖的改善趋势及持续性较弱，且退化区面积仍然较大，但是总体上是改善的，与王强等［32］的结论非常接近，并且沙漠以及荒漠草原均呈持续改善趋势，表明1999—2010年沙漠与荒漠草原均有转化为植被覆盖度较高的林草地的区域，植树种草、禁牧轮牧和防沙治沙生态恢复措施的广泛实施起到了重要作用［3］。

［1］ Foley J A，Levis S，Costa M H，et al.Incorporating dynamic vegetation cover within global climate models［J］.Ecological Applications，2000，10（6）：1620－1632.

［2］ IPCC.4th Assessment report of the intergovernmental panel on climate change［M］.Synthesis Report，2007：52.

［3］ 穆少杰，李建龙，陈奕兆.2001—2010年内蒙古植被覆盖度时空变化特征［J］.地理学报，2012，67（9）：1255－1268.

［4］ Tucker C J.Red and photographic infrared linear combinations for monitoring vegetation［J］.Remote Sensing of Environment，1979（8）：27－150.

［5］ Prince S D，Tucker C J.Satellite remote sensing of rangelands in Botswana（Ⅱ）：NOAA AVHRR and herbaceous vegetation［J］.International Journal of Remote Sensing，1986，7（11）：1555－1570.

［6］ Alcaraz－Segura D，Chuvieco E，Epstein H E，et al.Debating the greening vs.browning of the North American boreal forest：Differences between satellite datasets［J］.Global Change Biology，2009，16（2）：760－770.

［7］ Pettorelli N，Vik J O，Mysterud A，et al.Using the satellite－derived NDVI to assess ecological responses to environmental change［J］.Trends in Ecology ＆ Evolution，2005，20（9）：503－510.

［8］ White M A，de Beurs K M，Didan K，et al.Intercomparison，interpretation，and assessment of spring phenology in North America estimated from remote sensing for 1982—2006［J］.Global Change Biology，2009，15（10）：2335－2359.

［9］ Tottrup C，Rasmussen M S.Mapping long－term changes in savannah crop productivity in Senegal through trend analysis of time series of remote sensing data［J］.Agriculture，Ecosystems ＆ Environments，2004，103（3）：545－560.

［10］ Hüttich C，Herold M，Schmullius C，et al.Indicators of Northern Eurasia’s land－cover change trends from SPOT－VEG ETATION time－series analysis 1998—2005［J］.International Journal of Remote Sensing，2007，248（18）：4199－4206.

［11］ Symeonakis E，Drake N.Monitoring desertification and land degradation over sub－Saharan Africa［J］.International Journal of Remote Sensing，2004，25（3）：573－592.

［12］ Baldocchi D，Falge E，Gu L，et al.FLUXNET：A new tool to study the temporal and spatial variability of ecosystem－scale carbon dioxide，water vapor，and energy flux densities［J］.Bulletin of the American Meteorological Society，2001，82（11）：2415－2434.

［13］ Sparks T H，Aasa A，Huber K，et al.Changes and patterns in biologically relevant temperatures in Europe 1941—2000［J］.Climate Research，2009，39（3）：191－207.

［14］ Zhang Ke，Kimball J S，Mu Qiaozhen，et al.Satellite based analysis of northern ET trends and associated changes in the regional water balance from 1983to 2005［J］.Journal of Hydrology，2009，379（1／2）：92－110.

［15］ White M A，Running S W，Thornton P E.The impact of growing－season length variability on carbon assimila－tion and evapotranspiration over 88years in the eastern US deciduous forest［J］.International Journal of Biometeorology，1999，42（3）：139－145.

［16］ De Beurs K M，Henebry G M.Trend analysis of the Path finder AVHRR Land（PAL）NDVI data for the deserts of central Asia［J］.Geoscience and Remote Sensing Letters，IEEE，2004，1（4）：282－286.

［17］ Rogier de Jong，Sytze de Bruin，Allard de Wit.Analysis of monotonic greening and browning trends from global NDVI time－series［J］.Remote Sensing of Environment，2011，115（2）：692－702.

［18］ White M A，Thornton P E，Running S W.A continental phenology model for monitoring vegetation responses to interannual climatic variability［J］.Global Biogeochemical Cycles，1997，11（2）：217－234.

［19］ Jones M O，Kimball J S，Jones L A，Satellite passive microwave detection of North America start of season［J］.Remote Sensing of Environment，2012，123：324－333.

［20］ Moulin S，Kergoat L，Viovy N，et al.Global－scale assessment of vegetation phenology using NOAA／AV HRR satellite measurements［J］.Journal of Climate，1997，10（6）：1154－1170.

［21］ Zhang Xiaoyang，Friedl M A，Schaaf C B，et al.Monitoring vegetation phenology using MODIS［J］.Remote Sensing of Environment，2003，84（3）：471－475.

［22］ Reed B C，White M，Brown J F.Remote sensing phenology［M］∥Schwartz M D.Phenology：An Integrative Environmental Science.Dordrecht，The Nether－lands：Kluwer Academic Publishing，2003.

［23］ 徐建华.现代地理学中的数学方法［M］.北京：高等教育出版社，2002.

［24］ Hirsch R M，Slack J R，Smith R A.Techniques of trend analysis for monthly water quality data［J］.Water Resources Research，1982，18（1）：107－121.

［25］ Alcaraz－Segura D，Liras E，Tabik S，et al.Evaluating the consistency of the 1982—1999NDVI trends in the Iberian Peninsula across four time－series derived from the AVHRR Sensor：LTDR，GIMM S，FASIR，and PAL（Ⅱ）［J］.Sensors，2010，10（2）：1291－1314.

［26］ 邱海军，曹明明.基于SPOT VEGETATION数据的中国植被覆盖时空变化分析［J］.资源科学，2011，33（2）：335－340.

［27］ 程殿龙.2000年旱灾与抗旱工作［J］.防汛与抗旱，2001（1）：30－35.

［28］ 民政部救灾处.2001年全国自然灾害和救灾工作情况［J］.中国减灾，2002（1）：30－33.

［29］ 朴世龙，方精云.最近18年来中国植被覆盖的动态变化［J］.第四纪研究，2001，21（4）：294－302.

［30］ 李双双，延军平，万佳.近10年陕甘宁黄土高原区植被覆盖时空变化特征［J］.地理学报，2012，67（7）：960－970.

［31］ 何太蓉，庄红娟，刘存东.秦岭—黄淮平原交界带中东部近50年气候变化特征与趋势［J］.安徽农业科学，2009，37（14）：6532－6534.

［32］ 王强，张勃，戴声佩.基于GIMMS AVHRR NDVI数据的三北防护林工程区植被覆盖动态变化［J］.资源科学，2011，33（8）：1613－1620.