基于MODIS遥感数据的贵州省森林物候监测

摘要 以贵州省森林植被为研究对象，基于对MODIS卫星遥感影像的分析，探讨了遥感监测手段在区域尺度上提高森林植被物候监测精细程度的可能性。结果表明：通过分析比较归一化植被指数（NDVI）、增强植被指数（EVI）、相对绿度指数（Gcc）、绝对绿度指数（ExG）等不同植被指数，发现增强型植被指数（EVI）更能表达森林植被物候变化趋势；采用S-G滤波及动态阈值法提取物候关键期，发现贵州省森林植被物候变化趋势呈现西部地区植被生长开始期晚于东部地区，结束期早于东部地区，整个植被生长季长度西部短于东部；基于EVI的贵州省森林物候空间分异规律总体上不明显，仅纬度与生长开始期（SOS）和生长季长度（LOS）相关性显著，经度海陆地带性及海拔垂直地带性差异不显著。

关键词 森林物候；MODIS；植被指数；NDVI；EVI

中图分类号 S771.8  文献标识码 A

文章编号 0517-6611（2023）14-0114-06

作者简介 董奎（1990—），男，贵州织金人，工程师，硕士，从事林业资源调查研究。*通信作者，工程师，从事石漠化检查防治研究。

植物的萌发、展叶、开花、叶变色、落叶、飞禽迁徙、昆虫活动等物候现象，既反映了当地和当时的环境状况，也反映了过去一段时期环境状况的积累［ 1］。植物物候学是以年为周期，研究植物受气候、水文、土壤等生物因子和非生物因子的影响而产生自然现象的学科［ 2］。植物物候研究旨在认识自然季节现象变化规律，实现农林业、园林园艺业等经营管理的自动化、智能化、精细化，以及森林观花（果品）赏叶旅游的准确预报，服务于农（牧）林业生产、森林观花（果品）赏叶旅游、生态应用与全球气候变化研究、人类健康与科普教育等［ 3］。

植物物候现象是直观和敏感的环境条件季节和年际变化的生物指示器，其发生时间可反映陆地生态系统对气候变化的快速反应和适应，被认为是“自然界的语言”，是全球变化的“诊断指纹”［ 4-5］。精细的植物物候测定，通过植物物候相的变化节律，精细地反映植物生长节律和生态效应的波动变化，可以准确监测各种植物在每年生长季的起始时间、终止时间及生长季中的物候节律的精细变化。

传统的物候学研究方法主要是以野外观测为基础，即以眼观察和手记录的方式，通过定点定时记录生物物候现象的周年庆变化［ 6-7］，这种方法简单易行，操作方便，但耗时耗力，而且由于监测站点有限，在空间上很难实现由点到面的空间尺度转换，因此，监测站点有限，经常造成数据缺失，难以实现由点向面的空间尺度转换［ 7］。近40余年来，随着计算机技术和遥感技术的飞速发展，植物物候观测手段逐渐发展出一套从地面到空间卫星，广泛应用于植物物候监测及其与气候变化有关的研究，以获取植物物候连续、区域尺度信息的一套比较成熟的遥感监测技术和方法，从而使传统的物候监控手段的缺失得到了有效弥补。近几年，国内和国际植被物候研究以大田观测为主，以气候物候为基础，还包括以遥感为基础的大尺度物候研究［ 8-10］。笔者以贵州省森林植被为研究对象，探讨不同植被指标下的森林植被物候监测方法和途径，利用MODIS时间序列遥感数据进行遥感物候监测。

1 材料与方法

1.1 研究区概况

研究区位于贵州省103°36′～109°35′E，24°37′～29°13′N，位于我国西南部云贵高原东端（图1）。贵州省地势西高东低，自中部向北、东、南三面倾斜，平均海拔1 100 m，地貌属高原山地，其中92.5%为山地和丘陵。全省气候属亚热带湿润季风气候，年平均气温9～19 ℃，年降水量1 100～1 500 mm，降水量分布东、南偏多，西、北偏少，相对湿度在70%以上，全年日照时数1 300 h，无霜期270 d左右。降水占45%～52%。全省植被类型复杂多样，可分为针叶林、阔叶林、竹林、灌丛与草丛、沼泽植被和水草植被五大类，全省森林覆盖率在2021年达到62.12%。针叶林以杉木（Cunminghamia lanceolata）、马尾松（Pinus massoniana）、云南松（P.yunnanensis）、柏木（Cupressus funebris）等为主，阔叶林以壳斗科（Fagaceae）、樟科（Lauraceae）、木兰科（Magnoliaceae）、山茶科（Theaceae）等为主。

1.2 数据预处理

采用的卫星遥感数据为MODIS获取的遥感影像数据，该数据搭载于EOS/TERRA卫星上的中等分辨率成像光谱仪。精选贵州省2021年全年8 d最大值合成地表反射率产品（MODIS/TERRA8-DAYL3GLOBAL500M）MOD09A1，数据下载地址：https：//ladsweb.modaps.eosdis.nasa.gov/。贵州板块编号为：H27V06，数据包含7个波段，空间分辨率500 M。MOD09A1数据格式为HDF（HierarchicalDataFormat分级数据格式），投影方式为等面积正弦曲线投影，每张图像覆盖面积为1 200 km×1 200 km。资料已经辐射定标，几何精细校正，大气及气溶胶严格校正。全年共46幅（景物）影像，于2021年完成。

利用ENVI/IDL编写程序，利用ENVI/IDL函数File_Search和ENVI_CONVERT_File_MAP_Projection对MODIS数据进行批量的数据格式转换和投影转换，将正弦投影转换为经纬度坐标。利用贵州省行政边界矢量数据，投影转换后的数据采用ENVI_MASK_APPLY_DOIT函数进行批量裁剪。最后用IDL编写程序，按（1）～（4）对各种植被指数进行批量计算。

为了方便数据分析与存储，在IDL中利用fix（）函数转换各植被指数，将其扩大10 000倍，转换成整型。

森林植被数据来源于Globe Land 30 m空间分辨率全球地表覆盖数据（2020版），由国家基础地理信息中心发布［ 11］（http：//www.globallandcover.com/），空间分辨率为30 m，包括耕地、林地、草地等10个一级类型。其中，林地定义为乔木覆盖，树冠盖度在30%以上的土地，包括落叶阔叶林、常绿阔叶林、落叶针叶林、常绿针叶林、混交林，以及树冠盖度在10%～30%的疏林地［ 11］。

由于使用的MODIS数据空间分辨率为500 m，需要对地表覆盖数据进行重分类，使其与MODIS数据空间分辨率一致，再将其重分类为森林与非森林，以掩膜的方式处理为森林植被范围，得出研究区森林植被指标，在ENVI软件中完成这一过程。

1.3 植被物候确定方法

传统的物候被定义为一个间断性、不连续的过程，在特定的生命周期（如萌芽、抽枝、叶的展开、开花、结果和落叶、休眠等）的物种或个体尺度下的规律性现象。物候是指以卫星遥感为基础，以群落尺度的动态连续过程为基础，主要表现为随季节变化而变化的群落结构和外貌，即季相变化的物候节律。

森林物候参数提取是由Jnsson等［ 12］提出的动态阀值方法，又称比例阀值方法，即植被生长的起始期和结束期为20%，植被指数生长曲线的左右振幅。该方法所采用的阀值并不是某一特定的植被指标，而是采用了动态的比值形式［ 13］。物候关键期参数定义如下：生长季开始阶段（SOS）、生长季结束阶段（EOS）、生长季长度（LOS）、生长季最大值阶段（MOE）。

笔者所采用的8 d合成MOD09产品，经过严格的前期预处理，将8 d中最佳观测值作为影像像素值，采用最大值合成的方式进行处理。然而，经过处理的资料仍会受到云雾等影响，尤其是多云雾地区，如西南部，可能会在降低影像近红外线与红外线波段反射数值，导致资料在时间序列上突变，使植被物候分析合成期内持续出现杂音。因此，要想获得流畅的时间序列图像，重构时间序列植被指数是必不可少的一步。近年来，有更多的算法针对时间序列曲线去噪，如最大值合成法（MVC）、Savitzky-Golay滤波法（S-G）、双Logistic函数拟合法（D-L），AG（非对称高斯函数）、傅里叶分析法、小波分析法、中值迭代滤波法和时间序列谐波分析法（harmonic analysis of time series，HANTS）等。每种滤波方法都有优缺点，分别采用时代卫星软件模块中主要的S-G、D-L、AG 3种滤波方法和HANTS方法，对2015年EVI时间序列样本数据进行滤波，对原数据重建的效果进行不同方法的分析比较后发现［ 14］。原始曲线的细节保留在S-G滤镜的原始数据上。因此，笔者用Savizky-Golay（S-G）滤波法重构植被指数时间序列数据。S-G滤波法是应用最小二乘卷积拟合方法［ 15］，由Savizky和Golay于1964年提出。笔者采用ENVI/IDL编写程序进行植被指标重构，需要确定滤波窗的大小及在S-G滤波过程中的滤波窗多项式系数，并且滤波窗越大，时间序列曲线越平滑，但会导致细节信息的曲线变化而丢失。经连续测试，确定滤镜视窗为5的大小，滤镜多项式系数为3的大小。重构结果如图2所示（27°03′46.11″N，104°27′07.80″E）。

为更加精细地反映植被生长变化过程，获取一个生长季内植被各生长阶段的物候信息，对各植被指数在春季树木萌发和展叶生长期内分别按10%、30%、50%、70%、90%的速度递增进行测算。并记录各植被指数在秋季树木叶片变色衰老期下降10%、30%、50%、70%、90%和植被指数达到最大值时的天数。

2 结果与分析

2.1 植被指数选择

由各植被指数的时间序列变化曲线（图3）可知，MODIS卫星图像中各植被指数总体上可以反映出植被物候在年内的变化趋势。随着植被生长季的到来，植被指数（VI）曲线逐渐上升，当植被生长旺盛时，VI上升速度加快，曲线变化迅速，之后达到最大值，秋季植被逐渐进入凋落期，植物叶片开始变色，VI开始下降，且变化速度加快。但经比较分析可知，增强型植被指数EVI更难准确表征植被生长和衰老过程。

图3 不同植被指数曲线

Fig.3 Different vegetation index curves

2.2 植被物候关键期确定

利用TimeSAT软件，计算得出不同植被指数物候关键期（儒略日）（表1）（27°03′46.11″N，104°27′07.80″E）。由表1可知，植被指数EVI、NDVI、Gcc等指标与ExG植被物候关键期的指标存在差异。SOS在生长季的起始阶段为110 d左右，生长季结束期328 d左右，生长季长度218 d，生长季最大值阶段232 d。NDVI得出的SOS最小为96 d，Gcc最大为128 d；NDVI和ExG得出的EOS均为最大，为344 d，Gcc最小为304 d；生长季最大值阶段NDVI最大达248 d，EVI最小，为216 d。

设置不同植被指数变化率阈值，利用ENVI/IDL编写程序，计算不同变化率植被指数发生的天数（表2）。生长季植被指数分别提高10%、30%、50%、70%、90%时，衰老季植被指数下降了10%、30%、50%、70%、90%。从表2和图4可以看出，不同生长速度的日数在各植被指标上也出现了不同程度的背离。生长季节植被指数上升10%的天数，除NDVI异常值外，其余为80～100 d，而此时森林植物叶片尚未开始生长；植被指数在112～136 d上升到30%，树叶开始返青，也是草木群落物候返青期的开端；达到200 d以后，植被指数提高90%，这时大部分草木叶片展开定型，进入草木成熟期；至220 d左右时，草木指数达到最大值，此后一直处于稳定状态，并持续至260 d左右，草木叶片开始老化、变色，每片草木指数仅减少30%；到290 d以后，叶片开始大量变黄，植被指数下降70%左右，开始枯萎；至344 d以后，草木指数降至90%，几乎所有叶片枯萎，树木开始进入休眠状态。