作物物候期是重要的农业生态系统特征，准确获取作物的物候信息是农业生产、田间精细管理、计划决策等的重要依据，对于监测作物长势、进行作物种植管理、预测作物产量等具有重要意义^[1-3]。遥感技术可用于监测植被对气候的响应和作物物候期^[4-5]，利用遥感时间序列数据集提取物候信息在监测植被物候动态变化趋势^[6]、评估季节条件和季节变化下的植被覆盖和植被响应^[7]、预测作物产量^[8]、优化作物和生态等动力模型的状态变量^[9]等诸多方面有广泛应用。尽管以高空间分辨率、高光谱分辨率遥感数据为主要信息源的农情精准监测研究已有大量文献报道，但对大尺度区域监测而言，无论是从数据获取难易程度，还是从数据成本、时间效率等方面看，低空间分辨率的多光谱卫星遥感数据仍然具有应用优势，其中以AVHRR和MODIS等数据应用最为广泛^[10-11]，生成以归一化植被指数和增强型植被指数为代表的植被指数数据^[12-13]，因其具有良好的时间序列特性，可以较好地描述不同类型植被的生长过程，是物候期遥感监测中使用最多的数据类型。虽然该类数据的空间分辨率就中国中东部地区的农田规模化程度和农田破碎度而言相对较粗，但在相似的生态区内作物物候期进程相对一致，低空间分辨率影像可反映区域作物物候信息。由于植被指数时间序列受大气污染等的影响在时间上和空间上不连续，在物候研究前需对表征植被生长过程的遥感时间序列数据进行平滑，再根据去除噪声后的遥感植被生长曲线的特征提取物候信息。早期研究中，Gallo等^[14]利用AVHRR归一化植被指数时序数据监测玉米吐丝期等。近年来，该领域受到持续关注，相关研究不断发展，Hmimina等^[15]利用地面测量归一化植被指数数据评价了MODIS归一化植被指数用于预测不同生态系统物候期时的潜力，Kandasamy等^[16]对SG(Savitzky-Golay)滤波、不对称高斯拟合、Logistic拟合等3种去噪平滑算法的性能进行比较，评估了地表物候产品的不确定性，Verger等^[17]使用遥感估算的叶面积指数时序数据对大尺度物候信息进行了提取。

目前，物候遥感监测研究尚存在两方面问题：一方面可用于平滑植被生长曲线的方法包括SG滤波、Logistic拟合、三次样条拟合、小波去噪、高斯函数拟合、傅里叶分解等在内的多种噪声去除算法和模型^{[16, 18]}，已有学者在不同方法间进行比较并加以选择^[19-21]，SG滤波是最受推荐的方法之一^{[16, 22]}，为了提高滤波平滑质量，在SG滤波前通常需要进行时序数据的等时间间距合成，而滤波后数据的时间分辨率不发生改变，多日合成降低的时间分辨率影响物候期遥感监测精度；另一方面，针对不同类型的植被，其生长曲线特征点位需根据植被生长规律加以定义，确定其与具体物候期的对应关系，自然植被生长曲线的特征点位与其季节变化的对应关系可以较好地直观定义^[23]，但作物的不同物候期往往彼此临近且不同物候阶段的差异不如自然植被明显，现有研究中作物生长曲线特征点位与其物候期的匹配关系尚需进一步明确^[24-25]。针对上述问题，本研究利用MODIS归一化植被指数时间序列开展作物物候期遥感监测研究，以期进一步提高作物物候期遥感监测的准确率：采用最大值合成、改进的SG迭代滤波方法和Logistic拟合提高对作物归一化植被指数生长曲线重构的质量和重构后数据的时间分辨率，提取作物生长曲线的多个特征点位，通过比较分析不同特征点位对应日期与作物进入不同物候期的实际日期之间的匹配情况，建立作物生长曲线特征点位与作物物候期之间的最佳匹配关系。

本文以河南省为研究区，以夏玉米为研究对象。河南省位于中国中部偏东的黄河中下游地区，西起太行山和豫西山地东麓，南至大别山北麓，与冀、晋、陕、鄂、皖、鲁六省毗邻，东西长约580 km，南北跨约550 km，空间范围为31°23′~36°22′N，110°21′~116°39′E，土地面积为1.67×10⁵ km²。河南省内中部、东部和北部平原属于由黄河、淮河和海河冲积而成的黄淮海平原，黄淮海平原面积广阔、土壤肥沃，是中国重要的农耕地区，也是中国夏播玉米种植面积最大的区域，夏玉米也是河南省主要秋收粮食作物。研究区内农业气象观测站分布如图 1所示。

图1

图1 河南省农业气象观测站位置分布图 Fig.1 The distribution of agricultural meteorological observation stations in Henan Province

选用的遥感数据源为MODIS的地表反射率数据MOD09GA，数据级别为L2G(经过大气校正和几何校正)。使用空间分辨率为250 m×250 m的第1波段(中心波长为659 nm)数据和第2波段(中心波长为865 nm)数据进行归一化植被指数计算。MOD09GA产品使用的投影坐标系为正弦投影，应用前需进行坐标系转换。对2013，2014年和2015年共计1095幅MOD09GA数据进行投影转换及研究区域裁剪的批量处理操作，转换为等面积圆锥投影。研究区内夏玉米于5月下旬至6月上旬种植，9月下旬至10月上旬收获，为获取夏玉米全生育期生长曲线，选取每一年年积日为第121日—第304日(5—10月)的数据，经等面积圆锥投影、归一化植被指数计算和裁剪后，将其按时间顺序叠加在一起，构建2013—2015年夏玉米生长季的归一化植被指数时间序列数据集，并与基于多源数据利用决策树方法^[26]生成的河南省玉米种植区掩膜叠加。

夏玉米物候期地面观测数据取自河南省内的23个农业气象观测站，以2013年和2014年各站的物候期观测数据作为分析建模数据，以2015年各站物候期观测数据作为检验数据。每个观测站记录了其附近当年种植的夏玉米进入各个物候期的时间，观测的物候期包括播种、出苗、三叶、七叶、拔节、抽雄、吐丝、开花、乳熟和成熟期。

针对时间序列MODIS数据的特点，提出一种基于作物生长曲线重构和曲线特征点位匹配的夏玉米物候期监测方法。首先采用旬最大值合成、改进的SG迭代滤波和Logistic拟合对原始归一化植被指数时间序列数据进行平滑，而后使用动态阈值法和曲率法提取重构后作物生长曲线的多个特征点位，进而通过匹配特征点位与作物物候期监测夏玉米的关键物候期。

原始归一化植被指数时间序列数据存在大量噪声点，时间与空间不连续，无法直接使用，需对其进行去噪平滑。本研究选用旬最大值合成对原始归一化植被指数时间序列进行初步去噪，该方法对超过10 d的连续归一化植被指数低值噪声(通常由连续阴雨天气造成)只能减弱而无法消除。需要对最大值合成后的时序数据进行去噪处理。

使用SG滤波对最大值合成后的归一化植被指数时间序列进行平滑处理。SG滤波属于时间域内的去噪平滑算法，是一种利用滑动窗口进行局部多项式最小二乘法拟合的方法。这里使用一种考虑收割期限制的改进的SG迭代滤波方法以真实还原作物的归一化植被指数生长曲线，作者在前期研究中详细给出了该方法的处理流程^[27]。本研究利用改进的SG迭代滤波方法，在满足迭代终止条件时完成归一化植被指数曲线平滑。

经过SG迭代滤波后的归一化植被指数时间序列数据可用于描述夏玉米的生长过程，该序列的时间间隔为10 d，为了提高物候期提取的时间分辨率，对该时间序列进行Logistic曲线拟合^[28]。仅对夏玉米生长周期中归一化植被指数时间序列数值增加的阶段进行拟合，并将日期转换为计数值后输入Logistic函数计算得到作物生长阶段任意一日的归一化植被指数值。经Logistic拟合后，作物生长过程曲线的时间分辨率提高至1 d。

原始归一化植被指数时间序列数据经旬最大值合成、改进的SG迭代滤波和Logistic拟合后，原数据序列中的低值噪声点被去除，高值点则被保留，归一化植被指数曲线得到平滑且时间分辨率达到1 d。利用本研究方法平滑后的归一化植被指数时间序列数据与作物的真实生长过程更为吻合，为作物物候期的准确提取提供了更精准的数据支持。

完成夏玉米生长曲线重构后，使用动态阈值法和曲率法2种方法提取曲线的特征点位。使用动态阈值法时，参考Jonsson等^[29]的研究，根据作物生长规律，在归一化植被指数时间序列中数值增加的阶段，定义距离最小值为曲线增幅10%的位置为代表作物快速生长期开始的特征点位，记为动态阈值1；定义距离最小值为曲线增幅90%的位置为代表作物营养生长基本停止，转向以生殖生长为主的时期的特征点位，记为动态阈值2。使用曲率法时，定义第1个曲率最大值点位为代表作物快速生长期开始的特征点位，记为曲率最大值1；定义曲率最小值点位为代表作物生长最快的特征点位；定义第2个曲率最大值点位为代表作物营养生长基本停止，转向以生殖生长为主的时期的特征点位，记为曲率最大值2^[30]。以研究区内某一典型夏玉米像元为例，使用动态阈值法和曲率法提取曲线特征点位如图 2所示。

图2

图2 利用两种方法提取曲线特征点位 (a)动态阈值法，(b)曲率法 Fig.2 Curve feature point extraction using two methods (a)dynamic threshold method, (b)curvature method

研究区内每个夏玉米像元重构后的生长曲线均可利用动态阈值法和曲率法提取出5个特征点位，这些特征点位蕴含着作物生长的农学规律，需进一步明确其与作物具体物候期的对应关系。筛选2013—2014年23个农业气象观测站的实测物候期数据，去除4个异常站(存在缺测、误测，导致数据不完整)，使用剩余19个站的观测数据进行特征点与物候期匹配。将遥感数据与观测数据对比时，选取离农业气象观测站的空间位置最近且落入玉米种植区掩模内的像元(即最临近玉米像元，若有多个像元则取平均光谱)。将每个站点位置提取出的5个特征点位时间分别与实地观测的进入各物候期时间进行均方根误差计算，计算结果见表 1。某个特征点与某个物候期的匹配关系越好，均方根误差就越小。由表 1可知，均方根误差最小时，动态阈值1和曲率最大值1均对应七叶期，动态阈值1与七叶期对应关系更好；曲率最小值对应拔节期；动态阈值2和曲率最大值2均对应抽雄期，动态阈值2与抽雄期对应关系更好。最优匹配关系为动态阈值1对应七叶期，均方根误差为5.4 d，曲率最小值对应拔节期，均方根误差为6.4 d，动态阈值2对应抽雄期，均方根误差为6.0 d。分析对应的农学规律，在非胁迫条件下，抽雄期前归一化植被指数表征的以叶面积指数为代表的玉米群体生长参量持续增长，其中七叶期以前增长缓慢，故七叶期对应邻近归一化植被指数最小值的动态阈值1位置，之后玉米群体生长参量开始快速增长，至拔节期前后增长速率达到最大，故拔节期对应作物生长曲线斜率最大值即曲率最小值位置，抽雄期后叶面积指数基本停止增长，生物量持续增加，故抽雄期对应邻近归一化植被指数最大值的动态阈值2位置。

表 1

表 1 2013—2014年曲线特征点位与各物候期的均方根误差(单位：d) Table 1 Root mean square error between feature points and maize phenology in 2013 and 2014(unit: d) 物候期 动态阈值1 动态阈值2 曲率最大值1 曲率最小值 曲率最大值2 播种期 23.0 54.3 29.0 38.6 48.2 出苗期 16.9 48.0 22.8 32.3 42.0 三叶期 12.6 43.6 18.4 27.9 37.5 七叶期 5.4 32.4 8.2 16.9 26.4 拔节期 15.0 19.2 9.6 6.4 13.5 抽雄期 32.0 6.0 25.9 16.5 8.2 开花期 33.4 6.4 27.3 17.9 9.4 吐丝期 34.5 6.9 28.4 19.0 10.4 乳熟期 57.0 26.3 50.9 41.4 32.0 成熟期 79.8 48.6 73.7 64.0 54.5

表 1 2013—2014年曲线特征点位与各物候期的均方根误差(单位：d) Table 1 Root mean square error between feature points and maize phenology in 2013 and 2014(unit: d)

根据各特征点位与各物候期的均方根误差大小，选择利用动态阈值1提取七叶期，利用曲率最小值提取拔节期，利用动态阈值2提取抽雄期。筛选2015年23个农业气象观测站的实测物候期数据，去除5个异常站，使用剩余18个站的观测数据对遥感提取的3个物候期进行精度验证。经验证，动态阈值1提取七叶期的均方根误差为5.9 d，动态阈值2提取抽雄期的均方根误差为4.9 d，曲率最小值提取拔节期的均方根误差为5.3 d。3个关键物候期的遥感监测误差在6 d以内，各验证站点的误差分布情况见图 3。

图3

图3 2015年遥感提取生育期与实测生育期的误差分布图 Fig.3 Error distribution between extracted and observed maize entered growth stages in each site in 2015

利用文中方法对研究区2015年夏玉米物候期进行监测，图 4为利用遥感提取夏玉米进入七叶期、拔节期和抽雄期的时间分布情况。监测结果显示，研究区内夏玉米进入七叶期的日期以6月25日—6月27日为最多，该时段进入七叶期的像元数占夏玉米种植总像元数的45.6%；夏玉米像元进入拔节期的日期多分布在7月10日—7月17日，该时段进入拔节期的像元数占夏玉米种植总像元数的66.8%；夏玉米像元进入抽雄期的日期以7月27日—8月2日为最多，该时段进入抽雄期的像元数接近夏玉米种植总像元数的71%。从监测得到的夏玉米进入3个关键物候期的时间分布情况和空间分布情况看，七叶期和抽雄期的时空一致度较高，而进入拔节期的时间和空间变异度相对较大，这可能与拔节期前玉米生长速率受环境因素影响较大，且研究区夏玉米苗期经常遭遇初夏旱等灾害有关，由于各田块灌溉条件不同，不同的土壤墒情下进入拔节期的时间早晚不同。夏玉米花期则主要受光周期和积温影响，故同一七叶期对应的花期相对一致。本研究利用遥感监测研究区夏玉米进入关键物候期的时空分布结果符合农学规律。

图4

图4 2015年利用MODIS NDVI时序数据监测夏玉米进入关键物候期的时间 (a)七叶期，(b)拔节期，(c)抽雄期 Fig.4 Dates of summer maize entered key growth stages monitored by using time-series MODIS NDVI data in 2015 (a)seven-leaf, (b)jointing, (c)tasseling

本研究以MODIS为数据源对夏玉米的归一化植被指数生长曲线进行重建，提取夏玉米3个关键物候期结果的精度经验证在6 d以内，比权文婷等^[31]使用SG滤波法对小麦物候期监测的精度提高了4 d；与Sakamoto等的系列同类研究相比，比其采用小波变换方法对水稻物候期监测的精度提高了6 d^[1]，比其采用二步滤波法^[32]对玉米和大豆物候期监测的精度提高了1 d，比其采用形状模型拟合^[8]对美国35个州的玉米物候期监测的精度提高了3 d。

分析本研究方法较国内外同类研究对物候期监测精度改进的原因包括：①本研究综合采用旬最大值合成、改进的SG迭代滤波方法和Logistic拟合优化重构夏玉米归一化植被指数生长曲线，重构后的作物生长曲线既消除了噪声又符合作物的生长规律，同时数据的时间分辨率提高至1 d，而仅使用SG滤波、小波去噪等方法校正作物生长过程线时，虽然可以有效消除噪声，但未考虑作物生长规律进行曲线拟合，曲线形态的农学意义不够明确，并且使用离散点滤波方法校正后数据的时间分辨率不会发生改变，通常保持在MODIS植被指数标准产品的周期8 d左右^[1]，以此时间分辨率提取作物物候期精度受到制约。②不同于前人研究中往往直接定义作物生长曲线上的特征位置与作物物候期的匹配关系^[24-25]，本研究根据各特征点位与各物候期的均方根误差大小匹配特征点位与物候期，更具有客观性。

本文所使用MODIS数据的空间分辨率为250 m×250 m，在此空间分辨率上混合像元普遍存在，而混合像元的植被指数过程线不仅包含作物的生长信息，还包含有地表其他地物信息，这导致了作物生长信息提取中可能存在一定的不确定性；且混合像元中目标作物端元组分的丰度越低，像元光谱体现出的作物生长信息越少，不确定性越大。相比之下，遥感数据的空间分辨率越高，耕地像元中的混合信息越少，相同方法对物候期监测的精度越高，如杨浩等^[33]采用30 m×30 m分辨率的环境卫星CCD(charge coupled device camera, 电荷耦合器件相机)数据监测水稻的播种期、抽穗期、收获期，均方根误差减少为3.4 d。但在大尺度区域监测应用时，高空间分辨率遥感数据存在时空分辨率矛盾这一难题。后续研究中，拟应用多源遥感数据提高作物物候期监测的精度。此外，由于Logistics函数是单增函数，难以模拟花期后植被指数值下降的发育进程，故该函数不能满足对夏玉米全生育期的准确监测。在后续研究中，还需进一步对Logistics模型进行扩充或使用分段函数，以实现对作物全生育期的准确监测。

1) 综合采用最大值合成、改进的SG迭代滤波方法和Logistic拟合可提高对作物归一化植被指数生长曲线重构的质量，并将重构后数据的时间分辨率提高至1 d，为准确监测作物物候期奠定基础。

2) 采用动态阈值法和曲率法提取重建后作物生长曲线的多个特征点位，通过比较分析不同特征点位对应日期与作物进入不同物候期的实际日期之间的匹配程度，可建立作物生长曲线特征点位与作物物候期之间的最佳匹配关系。

3) 本文采用的方法基于时间序列MODIS数据监测夏玉米关键物候期的精度经验证在6 d以内，适用于大面积作物物候期遥感监测。