基于MODIS-LAI 数据的广西甘蔗物候期提取

谢鑫昌 ，杨云川,2,3*，田忆，廖丽萍,2,3，韦钧培，周津羽，陈立华,2,3

（1.广西大学土木建筑工程学院，广西 南宁 530004；2.工程防灾与结构安全教育部重点实验室，广西 南宁 530004；3.广西防灾减灾与工程安全重点实验室，广西 南宁 530004）

物候是指作物受气候变化及各种环境因素的影响而呈现出以年为周期的自然现象，当作物达到关键的生长发育状态所对应的时期即为作物物候期[1]。区域作物物候的动态变化是响应气候环境变化的重要指标[2]，即反映了作物的生长和发育规律，也影响着作物的产量和质量[3]。因此，及时、准确获取作物物候信息，对作物农情实时监测、田间精准管理、作物产量保障和气候响应研究等均具有重要意义[4-5]。

目前，地面监测和卫星遥感是及时获取作物物候信息的主要技术手段[6-7]。地面监测主要是在田间设立观测站点或手持光谱仪定时获取作物物候信息，该方式较为费时费力且局限于小范围田间管理[7]；相比而言，卫星遥感方式具有低成本、高时效、观测区域广及省时省力等多方面的显著优势[8-9]，已广泛应用于农业监测；遥感监测常通过植被指数或作物理化参数来反映大空间尺度上作物的光合作用状态及长势季节变化，最终实现作物物候期的过程识别[10]。

常用的遥感物候监测指标有归一化植被指数（NDVI，Normalized differential vegetation index）、增强植被指数（EVI，Enhanced vegetation index）及叶面积指数（LAI，Leaf area index）等。近年代表性的研究有：辛景峰等[11]、郑小波等[12]采用的条件时间内差法重构AVHRR、MODIS-NDVI 时序数据集，实现了冬小麦、水稻的关键物候期估测，其物候监测误差均小于1 旬。王尧等[13]、张喜旺等[7]、Chu等[14]、Liu 等[15]及Xu 等[5]在上述研究基础上采用S-G 滤波法重构MODIS-NDVI、HJ-NDVI 时序数据集，并通过动态阈值法提取了单双季稻、冬小麦、夏玉米、甜菜等多种作物的物候期，结果表明作物在空间上的物候分布与实地研究具有较好的一致性。上述方法不仅减少云雾对遥感数据质量的影响，而且降低了数据噪声，能有效提高作物物候信息监测的准确性。杨琳等[1]进一步采用非高斯对称、双Logistic 拟合法对MODIS-NDVI 时序曲线进行平滑降噪处理，实现了冬小麦关键物候期的提取，并通过对比论证了这两种方法与S-G 滤波方法的优劣。李艳等[3]则采用改进的最大值合成法重构MODIS-NDVI 时序数据，并进一步通过S-G 滤波降噪处理、Logistic 函数拟合，最后通过曲率法、动态阈值法实现夏玉米关键物候期提取。凌洋等[16]则采用HANTS、TSF 滤波法重构MODIS-EVI，提取了水稻的关键物候期，结果表明TSF 滤波法准确度高于HANTS 滤波，并将平滑和滤波处理综合，充分考虑作物生长的周期性规律和遥感数据自身特点。黄健熙等[4]基于S-G 滤波法重构MODISLAI，并结合双Logistic 函数拟合、历史积温法实现了冬小麦返青、拔节及开花期的提取，研究表明双Logistic 函数拟合精度较高。Huang 等[17]则基于重构的MODIS-NDVI 时序数据，采用有效积温法，实现了华北平原冬小麦抽穗期和开花期的提取，其中RMSE 均保持在4.5～6.5 d。上述研究表明，通过融入有效积温、太阳辐射等气象因素，能进一步提高作物不同关键物候期的识别精度。

综上所述，可通过多种时序重构法的对比分析来筛选具体区域、特定作物物候信息提取的最优方法；但已有研究主要应用于小麦、玉米、水稻等粮食作物，针对甘蔗等经济作物的物候期识别研究相对较少。广西属亚热带季风气候，水热条件充沛，是我国最大的甘蔗产区和蔗糖产业基地[18]。然而，广西多山地丘陵、云雾覆盖度高、农田地块破碎、插花种植普遍[19]，导致区域甘蔗种植面积虽广但平均单位面积产量不高、历年灾害损失严重。因此，及时、准确获取广西甘蔗的物候信息，是提升甘蔗品质和产量、优化甘蔗种植结构、统筹机械化收割等精准管理的迫切需求。本文即采用覆盖尺度广、时间分辨率高的MODIS-LAI 数据集[4,8]，选择S-G滤波函数、非对称高斯函数及双Logistic 函数等方法进行广西甘蔗物候的LAI 时序拟合，并通过动态阈值方法实现广西甘蔗关键物候期提取，对比论证各方法的识别精度差异，为广西甘蔗种植的智慧管理及蔗糖产业结构优化提供科学支撑。

1 数据来源与研究方法

1.1 数据来源

本研究采用覆盖广西区（图1（a））的MODISMOD15A2H 叶面积指数（LAI）数据产品，来源于美国航天航空局（NASA，https://ladsweb.modaps.eosdis.nasa.gov/），空间分辨率500 m，时间分辨率8 d，时段为2014—2018 年1—12 月，共230 期影像数据。文中使用的甘蔗种植空间分布数据来源于本课题组采用随机森林分类方法提取2014—2018 年的LANDSAT8 OIL（30 m 空间分辨率）遥感影像获得（图1（b）—（f）），该数据总体识别精度达到90%以上，Kappa 系数均超过0.8。此外，采用广西6 个农业气象站统计的2014—2018 年甘蔗物候期数据作为精度验证（表1）。

表1 广西甘蔗物候期统计结果（d）Table 1 Statistical results of sugarcane phenology in Guangxi (d)

1.2 数据处理与函数拟合方法

1.2.1 MODIS-LAI 数据预处理 MOD15A2H 数据的预处理主要采用NASA 官网软件MRT（MODIS Reprojection Tool）工具；首先对该数据进行重投影、拼接及镶嵌、数据重采样等处理，然后通过ENVIIDL 软件实现广西研究区及甘蔗种植区数据的裁剪。

1.2.2 函数拟合方法 MODIS-LAI 数据产品虽然经过8 d最大值合成，减弱了受大气中云和气溶胶的影响，但是由于卫星传感器角度、太阳高度角、水蒸气等因素的变化产生影像噪声[1]，致使LAI 数据产品的时序特征存在显著的不规则波动，难以真实地反映作物的时序特征。为此，文中对MODIS-LAI 时序数据进行必要的滤波处理，通过平滑降噪去空值处理的方法来削弱这种不规则波动的影响；具体选择了S-G 滤波函数法（Savitzky-Golay）、非对称高斯拟合法（Asymmetric Gauss）及双Logistic 函数拟合法（Double Logistic）等。

1）S-G 滤波函数法。该方法是由Savitzky[20]等提出的一种移动窗口式加权平均算法，其加权系数主要通过最小二乘法拟合得出，计算公式如下：

式中：Y指像元对应的LAI 原始值；Y*指对应像元通过滤波处理后的LAI 值；Ci为第i期LAI 值滤波过程的系数；j指原始LAI 数组的系数。

2）非对称高斯拟合法。Jonsson 和Eklundh[21]提出了一种通过局部时序曲线拟合拓展到整体拟合的平滑函数处理法，即非对称高斯拟合法，主要是通过不同的分段高斯函数模拟植被或作物的生长情况，最后采用平滑连接的方法实现整个时序曲线拟合，完成时间序列重建。局部拟合函数如下式（2）和（3）：

式中，g(t,a1,…a5)为高斯函数；a1为决定最大值和最小值所在整条时序曲线的位置；c1和c2为控制整条曲线的基准及振幅；a4、a5和a2、a3分别为控制曲线左、右部分的宽度及斜率。

式中：[tL,tR]表示时序曲线中待拟合部分的取值区间范围；fL(t)、fC(t)和fR(t)分别为该区间[tL,tR]内左边最小值、中间最大值及右边最小值所对应的局部拟合函数；α(t)和β(t)分别为介于0 到1 之间的剪切系数。

3）双Logistic 函数拟合法。双Logistic 函数拟合法是由Beck 等[22]提出，其拟合过程与上述非对称高斯拟合类似，都需要进行局部拟合再到整体拟合，其相比基于傅里叶变换的滤波函数拟合法，对植被或作物生长过程的拟合更为准确，对于高纬度地区的作物生长模拟也普遍适用，其拟合结果与非对称高斯拟合相似，不同在于所采用的局部拟合函数为双逻辑形式，且公式中缺少一个参数，仅有在大量统计对比中非对称高斯拟合对比该方法才会产生较小优势。其局部拟合公式如下：

式中：g(t,a1,…,a5)为高斯函数，a1为决定最大值和最小值所在整条时序曲线的位置，a2、a3和a4、a5分别为控制整条曲线左右半部分的宽度及斜率，最后通过整体拟合函数实现MODIS-LAI 时间序列重构。

1.3 物候期提取与精度检验方法

1.3.1 物候期提取 TIMESAT 是一款由Jonsson 和Eklundh[23]基于Matlab 程序开发的专门用来处理时间序列卫星遥感数据的软件，能很好地实现植被或作物的物候期动态特性分析，并获取短期甚至长期的植被或作物生长发育信息。该软件中能够实现S-G滤波函数法（Savitzky-Golay）、非对称高斯拟合法（Asymmetric Gauss）及双Logistic函数拟合法（Double Logistic）等多种时间序列降噪重构的平滑处理，本文主要采用动态阈值方法[23]并结合重构的作物生长曲线来实现广西甘蔗播种-萌芽期、茎伸期及成熟期的提取。其中本研究主要将20%阈值位置处确定为甘蔗的播种-萌芽期，也就是LAI 曲线上升阶段最大值的20%位置处所对应的时间即为甘蔗的播种-萌芽期开始时间；下降段最大值的20%处所对应的时间即为甘蔗成熟期开始的时间[23]；根据FAO作物生育期的规划方法[18]，其中作物生长中期对应甘蔗生育的茎伸期，即甘蔗LAI 曲线达到峰值所对应的时间即为茎伸期开始时间。

1.3.2 精度检验 本文主要通过图1（a）中6 个主要农业气象站点的2014—2018 年甘蔗物候期地面观测数据对广西甘蔗生育期遥感识别结果进行精度验证，分别采用最大误差、最小误差、平均误差及均方根误差等指标验证甘蔗物候期识别的准确度。均方根误差计算公式如下：

式中：di为广西甘蔗物候期提取的反演日期，为对应的广西甘蔗物候期的观测日期，d；n为样本总数。

2 结果与分析

2.1 拟合结果比较

选取广西崇左市扶绥县甘蔗种植区域内的某甘蔗像元，并查看该像元2014—2018 年LAI 时序重构曲线结果；由图2 可知，主要受到卫星传感器角度及水蒸气等不确定性因素影响，不论是哪年的LAI 原始时序曲线均呈现出显著的锯齿状波动，不能准确地反映整个甘蔗物候期的长势变化特征。因此，文中通过平滑函数处理方法可以很好地保持原始时序曲线的时间变化形态，并校正极低的离异值，使整个时序曲线拟合更接近甘蔗生长发育的真实情况。从图2 可以看出，S-G 滤波拟合效果较低，虽然消除了曲线不规则波动的影响，但局部曲线仍存在不规则波动变化，主要由于S-G 滤波处理的局限性导致的；根据杨琳等[1]、边金虎等[24]研究可知，S-G 滤波窗口大小及函数阶数的选取需要根据不同的条件进行主观人为地设定，如果选取的滤波窗口过大则容易遗漏掉较多的作物生长发育信息；反之偏小，则导致数据冗余，不能及时、准确地获取到整个物候期的长期长势变化趋势；S-G 函数的阶数使用过高，则容易出现曲线过拟合，且局部曲线数值就会产生不规则波动的情况。而非高斯对称及双Logistic 函数拟合的效果较好，消除了S-G 拟合存在的锯齿状波动情况，较好地模拟了整个甘蔗生育期的生长发育，并保持了S-G 滤波及原始数据长势季节性变化的趋势。

由图2（a）—（e）结合广西甘蔗的实际生长发育概况可知，广西甘蔗的播种-萌芽期一般在2—6 月份，因此前期甘蔗的LAI 值增长缓慢；6—8 月份广西甘蔗进入分蘖期-茎伸期后，其生长速度加快，生长力旺盛，LAI 指数变化加快，曲线呈现出陡峭上升趋势，直到8 月中旬左右达到最高值；9月下旬后进入甘蔗茎伸期末，生长缓慢，曲线呈缓慢下降趋势；11 月份前后进入成熟期，曲线下降速度加快，该甘蔗物候曲线变化趋势与丁美花[25-26]等在广西甘蔗长势监测及长势信息提取研究中的结论相一致，表明采用平滑函数降噪取空值的拟合方法可以较好地模拟广西甘蔗长势变化情况。

综上可知，S-G 滤波窗口及函数阶数是S-G 函数准确拟合及获取作物长势变化信息的关键因素；而非高斯对称及双Logistic 函数拟合的效果较好地弥补了S-G滤波处理中局部曲线值动荡变化的缺点，并保持了原始数据中作物长势季节性变化的趋势，能更准确地模拟广西甘蔗生长过程。

2.2 广西甘蔗关键物候期提取结果精度验证

文中采用6 个主要农业气象站点的2014—2018年甘蔗生育期地面观测数据对上述三种方法提取的广西甘蔗物候期结果进行了对比和精度验证（图3和表2）。

由图3 可知，2014—2018 年的甘蔗播种-萌芽期、茎伸期及成熟期的提取误差普遍都在15 d 以内，只有少数是超过15 d；其中2014—2018 年，基于S-G滤波拟合提取的甘蔗物候期误差超过15 d 的样本数为播种-萌芽期3 个（10%），茎伸期7 个（23.33%），成熟期9 个（30%）；2014—2018 年基于A-G 函数拟合提取的甘蔗物候期误差超过15 d 的样本数为播种-萌芽期0 个，茎伸期2 个（6.67%），成熟期6个（20%）；2014—2018 年基于D-L 函数拟合提取的甘蔗物候期误差超过15 d 的样本数为播种-萌芽期2 个（6.67%），茎伸期3 个（10%），成熟期8 个（26.67%）。

由表2 可知，基于S-G 滤波拟合提取的播种-萌芽期、茎伸期、成熟期开始时间的均方根误差分别为10.43 d、11.67 d、12.46 d，而平均误差分别为9.16 d、10.17 d、11.28 d；基于A-G 函数拟合提取的三个关键物候期开始时间的均方根误差分别为7.63 d、8.43 d、10.48 d，其中播种-萌芽期及成熟期精度较高，平均误差分别为6.69 d、6.91 d、9.13 d；基于D-L 函数拟合提取结果的均方根误差分别为9.02 d、10.01 d、11.76 d，平均误差为7.75 d、8.42 d、10.28 d。

综合图3 和表2 结果表明，S-G、A-G 和D-L对广西甘蔗的三个关键物候期提取结果的均方根误差及平均误差均在15 d 内，其中A-G、D-L 方法的提取结果优于S-G 方法结果，而A-G 结果则较D-L结果优，也即A-G 函数拟合方法对广西甘蔗的物候期识别效果最好，在广西甘蔗物候监测识别方面具有良好的应用价值。

2.3 广西甘蔗关键物候期空间分布格局

由3.2 节的甘蔗物候期提取精度分析可知，非对称高斯函数（A-G）拟合法对广西甘蔗物候期的提取效果最优，其中图5 为基于该方法对MODISLAI 时序曲线重构，并采用动态阈值法提取的2014—2018 年广西甘蔗播种-萌芽期、茎伸期和成熟期开始时间的空间差异分布图。

由图4 可知：2014 年广西甘蔗播种-萌芽期的开始时间呈“桂东向桂西推迟”的空间特征，桂东开始时间主要集中在第121～153 d（即4 月1 日到5 月2 日），桂西的开始时间则推迟至第153 d 后；2015—2018 年广西甘蔗播种-萌芽期开始时间的空间格局相似，整体趋于5 月2 日后。

2014、2016 年广西甘蔗茎伸期开始时间呈“桂西南向桂东北提前”的空间分布格局，由第281 d后向217 d前过渡（即10月8日后提前至8月5日前）；2015 年甘蔗茎伸期开始时间的空间格局则整体推迟至10 月8 日后；2017 年甘蔗茎伸期开始时间呈“桂东向桂西提前”的空间分布特征，即时间提前至8月5 日前；2018 年茎伸期开始时间则整体集中在8月5 日到10 月8 日。

2014 年广西甘蔗成熟期开始时间的空间特征与茎伸期相似，即由2014 年12 月11 日—2015 年1月12 日提前至2014 年11 月9 日—12 月10 日（即第346～377 d 提前至第313～345 d）；2015 年成熟期开始时间则由11 月9 日前（桂东北）推迟至12 月10 日后（桂西南）；2016、2018 年成熟期开始时间则呈“桂西向桂东提前”的空间变化特征，由2016年12 月10 日后提前至11 月9 日前。

综上可知，通过遥感技术提取的广西甘蔗物候信息能够很好地反映甘蔗生长发育的时空变化过程，对于广西甘蔗物候监测的年际、年内变化分析具有较好的应用价值。

2.4 物候提取结果的影响因素分析

表2 2014—2018 年广西甘蔗物候期提取精度评价（d）Table 2 Accuracy evaluation of sugarcane phenology inversion in Guangxi from 2014 to 2018 (d)

综上可知，基于三种方法重构MODIS-LAI 时序数据，并通过动态阈值法提取的广西甘蔗物候期存在显著差异，这与三种序列重构法本身的特点密切相关。首先，A-G 和D-L 拟合方法均属于函数拟合法，通过局部拟合扩展至整体，受到拟合函数自身形式的影响较大，可获得平滑的重构曲线；而S-G滤波属于时域滤波法，属于局部处理方法，时序曲线的平滑程度与滤波窗口及函数阶数的选取有关，人为主观性较大，受原始数据质量的影响大。其次，A-G 函数和D-L 函数拟合的基本原理相同，仅是局部拟合函数不同，导致两种方法无论是时间还是空间上得到的甘蔗物候期模拟结果非常相似；S-G 滤波法的原理和前者完全不同，因此，造成其与A-G和D-L 函数拟合法提取甘蔗物候期的时空信息存在显著差异，而该曲线上升及下降阶段的左右偏移也是造成物候期出现提前或推迟的主要原因。

虽然通过AG 函数提取的广西甘蔗物候期效果最优，但其估测结果仍与地面观测值存在误差，而导致该误差的主要原因有：

1）数据方面，MODIS-LAI 数据受空间分辨率限制（空间分辨率500 m），其值仅代表像元平均值，而地面观测数据来自农业站点的物候观测值，二者空间坐标匹配存在明显差异，影响了遥感提取的物候期与观测值空间分辨率的一致性；但农业气象站点并不是在甘蔗田间的单点设置，而是在当地甘蔗种植区范围内布置多个观测站点，通过多站点的共同协调作用实现当地甘蔗物候监测，所以站点布置的观测范围与MODIS-LAI 的空间分辨率接近，采用该遥感数据进行广西甘蔗物候信息监测具有一定的合理性。

2）广西属于亚热带季风区，降雨充沛，在甘蔗等作物生育季期间云雨较多，对于受云雾影响的区域，甘蔗物候期内LAI 值普遍偏低，时序曲线波动变化显著，不能反映真实的甘蔗生长发育状况；而时序曲线重构法具有较好的鲁棒性，可更好地提高影像数据质量，并获得连续完整的甘蔗区LAI 时序曲线，使其更接近甘蔗的物候生长过程，具有一定的合理性；但其LAI 模拟值仍与甘蔗的实际值存在一定误差，因此，导致遥感提取的甘蔗物候期与真实情况对比存在高估或低估现象。

3）广西山地丘陵偏多，地势复杂，农田地块破碎，插花种植现象严重，导致甘蔗种植区解译结果偏少或偏多，且影像获取的蔗区位置相对偏移，与真实的甘蔗空间种植结构分布存在一定偏差，进而影响甘蔗的物候提取结果；且山地丘陵区也属于云雨及大雾天气易发区，严重影响遥感数据质量，与2）中影响效应相同。

随着遥感技术的发展，GF-1、GF-6、Sentinel-2、ZY-3 等中高分辨率影像的日益丰富，有效积温、太阳辐射等有关作物生长因素的融入，有助于实现更加精准的作物物候信息提取。因此，采用遥感提取的作物物候信息与地面观测相结合是实现作物物候信息监测的有效途径，也是未来农业资源遥感监测的研究热点之一。

3 结论

本文基于MODIS-LAI 数据，通过三种时序重构方法实现了广西甘蔗关键物候期提取，得到以下主要结论：

1）针对MODIS-LAI 数据时序曲线存在锯齿状波动而难以准确反映广西甘蔗生育期长势变化的问题，S-G 滤波函数法、非对称高斯拟合法和双Logistic 函数拟合法等均可有效消除该序列的不稳定波动及其奇异值，但后两种方法的拟合结果较好地弥补了S-G 滤波处理中局部曲线值动荡变化的缺点，能更准确地获取广西甘蔗生长过程信息。

2）在广西甘蔗物候期提取精度方面，三种方法提取结果的均方根误差及平均误差普遍在15 d内；其中，A-G 方法的精度最高，其次为D-L 方法，S-G 方法的提取结果精度相对最差。

3）由于S-G 函数法的滤波窗口及函数阶数选取的主观性和A-G、D-L 拟合法自身函数形式的限制，导致三者在时空上获取的广西甘蔗物候期模拟结果存在明显的差异性；地面观测数据与遥感数据的空间一致性，云雨及大雾天气导致遥感数据质量的下降及广西复杂地势条件对甘蔗种植分布解译结果的干扰等均是造成广西甘蔗物候期估测结果与地面观测值存在误差的主要原因。