基于MODIS NDVI数据的郑州市植物物候监测研究

杨 阳 何瑞珍 田国行 任晓娟 李 岩 郑景飚

(1. 河南农业大学林学院，河南 郑州 450002；2. 上海工程技术大学电子电气工程学院，上海 201620；3. 信息工程大学地理空间信息学院，河南 郑州 450001)

植物由于受气候等环境因子影响形成的生长发育节律称为植物物候，因其能直接响应自然气候变化的季节、年际规律，被作为重要的生态指示器广泛应用于农业生产、全球气候变化、生态学实践等[1]。以往的物候监测一般采用目视观察法，直接在野外定点观察和测量生物物候的年内、年际变化[2]。该观测方法具有准确性和客观性的特点，但难以进行长周期、大尺度的物候时空变化监测和分析，同时会耗费大量的人力和物力。而日益发展的遥感技术则为研究物候提供了新的手段。

光谱曲线反应了物体对波的吸收和反射状况，其形状取决于物体内在的组织结构特征。因此，不同物体各具不同的光谱曲线，而同类物体由于内部构造相近形成了相似的光谱曲线形状，如植被特有的光谱曲线。遥感技术正是基于此原理实现对植物物候的监测。遥感监测覆盖范围广、时间序列较长，同时满足空间连续和多时相监测，对于指示植被动态变化与全球气候变化的关系具有重要意义[3]。

遥感植被监测是利用遥感植被指数与研究区植被参数的相互关系，以植被指数的变化来反映植被状况[4]。NDVI是基于MODIS数据的归一化植被指数，其利用绿色植物的波谱特性能在大范围覆盖区域内反映植被覆盖度、光合有效辐射等植被参数。Zhang等[5]采用MODIS 时间序列数据提取美国东北部植被的物候期，并分析了不同植被物候期随纬度变化的情况。刘玲玲等[3]利用2005年AVHRR和MODIS的NDVI数据研究了全球植被物候参数，并通过比较分析，证明2种传感器提取的结果在空间变化趋势上具有一致性。这些研究均表明NDVI能有效的评估植被生长状态[6]，可应用于监测植被、分析物候等。

本研究结合郑州市域植被类型的分布特点，基于MODIS NDVI时间序列数据进行分析，采用滑动平均法获取郑州市域植被关键物候期，以期实现郑州主要植被类型物候的实时监测，达到快速获取区域植被信息的目的，对评估中原地区的气候变化及其对植物的影响，指导当地的农业生产，监测生态环境变化等具有重要的理论和现实意义。

1 研究区概况

郑州市位于东经112°42′～114°14′，北纬34°16′～34°58′，地处中原，北临黄河，西依嵩山，东南为广阔的黄淮平原。全市总面积7 446.2 km2，市区面积1 010.3 km2，农用地、园林地和草地面积4 562.74 km2。郑州属于北温带大陆性季风气候，春季时间为每年的3月下旬至5月下旬；夏季为5月下旬至9月初；秋季为9月初至11月初；冬季为11月初至次年3月下旬。夏季炎热、降雨集中，秋高气爽、日照充足，冬季寒冷、降雪少。年平均气温15.6 ℃，年平均降雨量为542 mm，全年日照时间约1 869.7 h。

本研究根据郑州市域的植物分布情况，以林地、草地、农田3种代表性植被类型为对象进行研究。由于农田受气候、土壤等自然因素影响，不同地区的农田会进行各种作物轮种，出现了一年一熟、一年两熟，甚至于一年多熟的情况[7]。郑州的农田基本上为一年两熟，以旱地作物为主，包括玉米 (Zeamays)、小麦 (Triticumaestivum)、棉花 (Gossypiumherbaceum)、大豆 (Glycinemax) 等。本研究对于农田的复杂情况，以大面积种植的冬小麦和夏玉米为主要研究对象，进行植物物候的监测。

2 材料与方法

2.1 数据来源

本研究采用的MODIS陆地产品是搭载于EOS/Aqua卫星上的MODIS获取的数据，由NASA的EOS数据中心提供的16 d合成数据 (MYD13A1)。MYD13A1是陆地的标准数据产品，内容主要为NDVI，空间分辨率为500 m，时间范围为2011年全年。GlobCover为2009年全球陆地覆盖数据，分辨率为300 m。选取MERIS传感器在2009年1月1日至12月31日期间所接收的较高质量的影像数据，来进行图像合成，其采用 “美国食品和农业组织的地表覆盖分类系统” 作为图例生成标准。

2.2 时间序列谐波分析法

云和大气能够降低NDVI数值，是影响遥感影像最大的噪声[8]。本研究采用时间序列谐波分析法 (HANTS) 对NDVI数据进行滤波处理，其核心算法是傅里叶变换和最小二乘法拟合，即把时间序列波谱数据的复杂函数分解成若干正弦、余弦函数，通过叠加其中具有时间序列特征的多条曲线，对时序数据进行平滑。该方法在保持原时间序列周期性的同时，又能使异常像元恢复正常[9-11]。在NDVI时间序列数据中，受大气干扰及处理误差影响的谐波频率一般较高，而反映植被相关信息的谐波频率一般较低。因此平滑重建较低频率的谐波即可减弱噪声影响，重构NDVI时间序列曲线[12-13]，并真实反映曲线的周期性变化规律。

2.3 物候划分

将郑州植物物候分为3个关键期：开始期、成熟期 (抽穗期)、结束期。开始期指植物开始生长或绿度开始增加的时间，为NDVI值开始增大，NDVI时间序列曲线上升速率迅速增大处；成熟期为植物绿度和光合作用最大化的时间，NDVI值达到最大；结束期为植物绿度和光合作用迅速下降的时间，NDVI时间序列曲线下降速率变化最大处；生长季的长度为结束期和开始期之间的差值。对于农作物一年两熟的情况，根据对农作物的辨别和轮作的规律，分时段进行研究。

2.4 滑动平均法

采用滑动平均法进行植物物候关键期的提取。其基本原理是利用滑动平均后的曲线与实际NDVI时间序列曲线相交来获取植物物候。把实际NDVI曲线第1次超过滑动平均曲线的时间作为开始期，第1次低于滑动平均曲线的时间作为结束期[14]，实际曲线最大值处对应的日期即为成熟期。该方法对NDVI时间序列数据的计算更为稳定、可靠，但其准确性与滑动平均时间间隔的选取密切相关，而滑动平均时间间隔又与不同植被类型特征有关，需依研究植被实际状况而定。

3 结果与分析

3.1 郑州市域土地分类结果与验证

对GlobCover数据进行裁剪得到郑州市域地表分类数据，通过地表覆盖分类系统中的唯一值解译出农田、林地、草地3种主要的土地类型，可得到郑州市域主要植被类型的分布图 (图1)。由图1可知，郑州市域内土地利用以农田为主，分布最广。林地次之，与地形相结合，集中在西部、南部山区。草地较少，主要在北部黄河沿岸分布，同时在西部山区中草地与林地交错分布。利用Google地球高分辨率影像数据对其进行抽样验证，结果表明其能满足郑州市域尺度的分类精度。另外，采用2009年全球陆地覆盖数据与2011年MODIS数据进行综合分析，由于两者的获取时间不一致，在采样中同时借助Google地球高分辨率影像数据进行数据比对调整，以避免随着城市的发展，草地、农田和林地的分布变化对研究结果造成影响。

图1郑州主要植被类型分布
Fig.1 Distribution of main plant types in Zhengzhou

3.2 滤波效果比较

频率个数、曲线匹配阈值 (FET)、剩余点个数 (DOD) 等参数是影响曲线平滑效果的重要因子，其设置没有客观标准。本研究经过反复试验，取频率个数为2，FET为1 000，DOD为8时，拟合效果最好。郑州林地采样点NDVI滤波前后的对比结果见图2。由图2可以看出，滤波前NDVI时间序列曲线呈无规则剧烈波动，难以进行分析。滤波重构的NDVI时间序列曲线排除了干扰值同时保持了曲线特征。经HANTS滤波的NDVI时间序列曲线可直接用于分析植被物候期。

图2HANTS滤波前后NDVI曲线对比
Fig.2 NDVI curve comparison before and after HANTS filtering

3.3 物候关键期的提取

基于研究对象特点及分布，采用样点统计分析方法。将郑州市域主要植被类型分布图与处理后的研究区NDVI影像叠加，在草地、农田分布区域内分别随机选取20个采样点。草地采样点植物主要包括：白头翁 (Anemonechinensis)、东风菜 (Doellingeriascaber)、卷丹 (Liliumlancifolium)、小糠草 (Agrostisgigantea)、马唐 (Digitariasanguinalis)、胡枝子 (Lespedezabicolor)、复盆子 (Rubusidaeus) 等。考虑到林地在城市附近因人为的引种、选育、栽培各种林木而物候期差别甚大。为满足本研究要求，故在以乡土树种为主的郑州西部山区天然林地内随机选取非常绿植被采样点20个。林地采样点乔木主要包括：枫杨 (Pterocaryastenoptera)、刺槐 (Robiniapseudoacacia)、枣树 (Zizyphusjujuba)、黄栌 (Cotinuscoggygria)、栎属 (Quercus)、硬阔、软阔等。逐一分析各类型植物的物候关键期。

经过反复试验，取林地、草地的滑动平均时间间隔为5个16 d合成周期 (图3)。因小麦、玉米的生长季长度较短，农田滑动平均时间间隔取3个16 d合成周期时效果最好 (图4)。根据各植被类型的最适滑动平均时间，利用滑动平滑法对各植被的物候关键期进行提取，提取结果见图5。

图3滑动平均前后林地NDVI曲线对比及物候参数提取
Fig.3 Comparison of NDVI curves and extraction of phenological parameters of forest land before and after the moving average

根据各植被物候期分布，可知郑州西部山区林地大部分在第85～110天开始生长，在第215～240天生长达到最大程度，在第305～325天停止生长，林地生长季长度为200～220 d，生长时间较长 (图5a)。大部分草地在第85～110天开始生长，在第220～240天生长旺盛，在第300～325天停止生长、枯萎，草地生长季长度为190～220 d (图5b)。大部分冬小麦开始期在第40～60天，处于拔节期。在第100～120天生长达到极限，在第135～150天逐渐停止生长，进入收获期，小麦生长季长度集中在90～110 d (图5c)。玉米在小麦收割后进行播种，与小麦采用轮作模式。大部分玉米生长期集中在第165～180天，在第220～235天进入抽穗期，在第260～270天生长结束，进入收获期，生长季长度最短，玉米生长季长度为90～105 d (图5d)。

图4滑动平均前后农田NDVI曲线对比及物候参数提取
Fig.4 Comparison of cropland NDVI Curves and extraction of phenological parameters before and after the moving average

在郑州市域内不同地点，同类型植被的物候关键期略有波动，但总体保持一致。天然林地、草地形成了与郑州自然气候变化相适宜的物候期。农作物小麦、玉米在人为作用下依自然季节的更替进行有规律的轮作，在各气候区均形成了相对稳定的物候期。

3.4 物候关键期的验证

验证是基于遥感数据提取大面积物候信息的重要工作。柳晶等[15]的研究中，春季物候主要为前1～3个月的影响，平均温度每升高1 ℃，物候期提前2.8～7.9 d、2.7～9.0 d；对于生长季，年平均温度每升高1 ℃，物候期延长9.5～18.6 d，结束期推迟。

把刘玲玲等[3]对该区域的物候期研究作为真实值与本研究得出的物候期进行对比分析，结果见表1。经查询郑州2011年前3个月的均温比2005年高1.17 ℃，对应开始期提前12 d左右。2011年的年平均气温相较2005年升高0.2 ℃，对应结束期推迟20 d左右，生长季长度延长35 d左右。可以发现，本研究中各项温度均有所升高，物候期提前，结束期推迟，生长季长度延长，与上述研究结论基本是一致的。但由于环境条件改变及存在的可能误差，物候期的延长天数相较柳晶的研究有所增加。另外，根据对玉米和小麦的实地调查发现，小麦返青期基本从2月开始,收割期为5月底。玉米播种基本上为6月初小麦收割后，玉米收获通常在9月底、10月初进行。通过对比以上资料，看出基于MODIS NDVI数据获得的植物物候期与研究和调查资料具有一定的可比性，表明可利用遥感来监测植物关键物候期。

图5各植被物候期分布
Fig.5 Distribution of vegetation phenology

表1 2005年与2011年物候期对比Table 1 Contrast of phenology between 2005 and 2011

4 结论与讨论

本研究通过时间序列谐波分析法 (HANTS) 平滑重构NDVI时间序列曲线，有效减弱了噪声的干扰，能直接用于提取植物的物候信息。同时结合全球陆地覆盖数据，解译出郑州市域主要植被类型分布图，并针对林地、草地和农田的不同特征，运用滑动平均法分别进行各类型植物采样点的滑动平均分析，实现快速、有效监测郑州市域植物关键物候期，包括开始期、成熟期、结束期以及生长季长度。通过与前人相关研究及实际调查的物候资料比较验证，结果基本吻合，表明该方法获得的物候期具有一定科学性。

对比于信芳等[1]与刘玲玲等[3]的研究，均表明通过HANTS滤波处理的MODIS NDVI时间序列数据可直接用于植物关键物候期的监测。而与本研究采用滑动平均法不同，两者均采用动态阈值法进行物候关键期的提取，该方法通过设定阈值条件，将植物物候期限制在合理生育期内，但阈值本身的确定易受主观影响。针对本研究，滑动平均法对物候的计算则更为稳定、可靠，但容易受春季雪融的影响，所得结果可能早于植物实际返青期。研究中已查询郑州地区2011年历史天气，发现不存在大量积雪状况，可排除NDVI数据可能受春季雪融的影响。另外，两者均通过逐个像元运算得到区域内的关键物候期分布格局，相较本研究的随机采样统计分析，更加直观的反应了整个区域关键物候期的分布特点，可为本研究的深入提供借鉴。

本研究通过将NDVI时间序列数据与全球陆地覆盖遥感数据相结合，实现了不同植被类型物候的分类研究、比较。但受遥感影像时间分辨率、空间分辨率的限制，仅研究了植物的3个物候关键期，下一步可利用高精度影像详细研究植物各阶段的物候特征。另外，随机采样可能会对提取的植被物候精确度产生影响，今后也可通过提取多年的物候平均值来做相关性分析将会更加科学。

[参 考 文 献]

[1] 于信芳, 庄大方. 基于MODIS NDVI数据的东北森林物候期监测[J]. 资源科学, 2006, 28(4): 111-117.

[2] Wang Q, Tenhunen J D. Vegetation mapping with multitemporal NDVI in North Eastern China Transect (NECT) [J]. International Journal of Applied Earth Observation and Geoinformation, 2004, 6(1): 17-31.

[3] 刘玲玲, 刘良云, 胡勇. 基于AVHRR和MODIS数据的全球植被物候比较分析[J]. 遥感技术与应用, 2012, 27(5): 754-762.

[4] 强建华, 赵鹏祥, 陈国领. 基于NDVI的油松天然林生长状况的遥感监测研究[J]. 西北林学院学报, 2007, 22(1): 149-151, 172.

[5] Zhang X Y, Friedl M A, Schaaf C B, et al. Monitoring vegetation phenology using MODIS[J]. Remote Sensing of Environment, 2003, 84(3): 471-475.

[6] 朱源, 彭光雄, 王志, 等. 西藏林芝地区近30 a来的NDVI变化趋势研究[J]. 西北林学院学报, 2011, 26(4): 69-74.

[7] 张峰, 吴炳方, 刘成林, 等. 利用时序植被指数监测作物物候的方法研究[J]. 农业工程学报, 2004, 20(1): 155-159.

[8] 侯学会, 牛铮, 高帅, 等. 基于SPOT-VGT NDVI时间序列的农牧交错带植被物候监测[J]. 农业工程学报, 2013, 29(1): 142-150, 294.

[9] Julien Y, Sobrino J A. Comparison of cloud-reconstruction methods for time series of composite NDVI data[J]. Remote Sensing of Environment,2010,114(3): 618-625.

[10] Menenti M, Azzali S, Verhoef W, et al. Mapping agroecological zones and time lag in vegetation growth by means of fourier analysis of time series of NDVI images[J]. Advances in Space Research, 1993, 13(5): 233-237.

[11] Roerink G J, Su Z, Menenti M. S-SEBI: a simple remote sensing algorithm to estimate the surface energy balance[J]. Physics and Chemistry of the Earth, Part B: Hydrology, Oceans and Atmosphere, 2000, 25(2): 147-157.

[12] 李正国, 杨鹏, 周清波, 等. 基于时序植被指数的华北地区作物物候期/种植制度的时空格局特征[J]. 生态学报, 2009, 29(11): 6216-6226.

[13] Julien Y, Sobrino J A, Verhoef W. Changes in land surface temperatures and NDVI values over Europe between 1982 and 1999[J]. Remote Sensing of Environment, 2006, 103(1): 43-55.

[14] 崔凯. 基于遥感技术的作物物候监测方法及动态变化分析研究[D]. 长沙: 中南大学, 2012.

[15] 柳晶, 郑有飞, 赵国强, 等. 郑州植物物候对气候变化的响应[J]. 生态学报, 2007, 27(4): 1471-1479.