杨琳　高苹　居为民

摘要：农作物物候期是作物对气候变化响应的重要指标，也是提高区域农业管理水平的重要参数；遥感技术的发展为大面积提取农作物物候期提供了有效的技术途径。利用2010年MODIS NDVI序列提取江苏省冬小麦关键物候期的方法，首先采用非对称高斯函数、双Logistic函数对NDVI时序曲线进行拟合，进而采用动态阈值法确定冬小麦物候期，并利用观测数据对提取结果进行验证。结果表明，就冬小麦返青期、抽穗期、成熟期的提取而言，双Logistic函数拟合方法略优于非对称高斯函数拟合方法，利用前一种拟合方法提取的返青期、抽穗期、成熟期与观测数据比较的均方根误差分别为5.5、9.4、7.5 d，而利用后一种拟合方法提取的返青期、抽穗期、成熟期均方根误差分别为6.1、9.5、7.8 d；2010年全省的冬小麦普遍于第49天之前开始返青，抽穗期普遍开始于第105天至第113天，成熟期普遍开始于第145天至第161天，抽穗期、成熟期大体上表现出从南到北逐渐延迟的趋势，但是返青期的空间差异并不明显。

关键词：江苏省；冬小麦；物候期；生育期；MODIS NDVI；时间序列

中图分类号：S127 文献标志码：A 文章编号：1002—1302（2016）01—0315—06

物候是指自然界中的生物受气候和其他环境因素的影响而出现的以年为准周期的自然现象，农作物的物候期则是指农作物达到关键生育期所对应的日期。农作物物候的变化反映了气候变化对农作物生长发育的影响，是气候变化对农业生态系统影响的指标之一，因此成为全球气候变化研究的热点问题之一。

作物物候信息可以通过田间观测、积温预测和遥感监测等多种途径获得。但是田间观测方法不仅费时费力，而且代表的面积有限；积温预测方法需要播种日期、逐日气温以及某一生育期所需积温等信息。因此，这2种方法的大面积应用受到一定的限制。遥感技术的发展为大面积监测作物物候提供了崭新的手段，基于遥感信息计算的植被指数能够在大范围内较精确地反映植被绿度、光合作用和生长的季节和年际变化，因而被广泛应用于植被监测、分类和物候分析中。目前常用于物候监测的植被指数有归一化植被指数（normalized vegetation index，NDVI）、增強植被指数（enhanced vegetation index，EVI）等。有研究认为，NDVI作为指示地面生长季节开始和结束的遥感指标具有比较稳定的特点，因此国内外很多学者利用NDVI时序数据对作物物候及其时空变化进行了研究。如Sehgal等利用1981—2001年NOAA-AVHRR NDVI数据集对印度恒河平原的农作物物候信息进行提取并分析了其时空变化特征。常守志等基于MODISNDVI时序数据采用动态阈值法提取了三江平原2003年的农田物候信息。

冬小麦是江苏省主要的夏熟作物之一，其生长发育经常受到低温冻害、连续阴雨和病虫害的影响。这些灾害所造成的损失与灾害的强度、持续时间以及冬小麦所处的生育期有关，准确的冬小麦物候期分布数据有助于提高农业生产管理水平、减轻灾害的影响。但是目前利用遥感数据提取江苏省冬小麦物候期的研究还有待开展。因此，本研究分别利用非对称高斯函数、双Logistic函数方法对MODIS NDVI数据进行拟合，采用动态阈值法提取江苏省2010年冬小麦的关键物候期，比较利用2种拟合方法提取的冬小麦物候期差异，以期为利用遥感技术监测冬小麦的生长和发育进程提供参考。

1研究区概况

江苏省（30°45′～35°20′N，116°18′～121°57′E）位于长江、淮河下游，总面积10.26万km²，地形以平原为主。江苏省属东亚季风气候区，处在亚热带和暖温带气候过渡地带，全省年均降水量700～1200 mm，年均气温13～16℃，年日照时间2000～2 600 h，平均无霜期200～240 d。大多数地区实行一年两熟制，粮食作物以水稻、麦类为主，经济作物有棉花、花生、油菜、黄麻、蚕桑等。其中，冬小麦是主要的夏熟作物之一，一般10月份播种，翌年5月份成熟。

2材料与方法

2.1数据来源与预处理

2.1.1数据来源 研究采用的覆盖江苏省的MODISMOD09Q1反射率数据产品来源于NASA（htp：∥ladsweb.na-scom.nasa.go.index.html），空间分辨率为250 m，时间分辨率为8 d，时间范围覆盖了2010年1—12月，共46期数据。研究使用的辅助数据有江苏省行政边界数据、土地利用数据，以及用于结果验证的6个农业气象观测站的2009—2010年的冬小麦生育期数据（图1）。其中，土地利用数据由中国科学院南京地理与湖泊研究所通过融合2006—2007年SPOT2/4卫星遥感影像（全色，10 m空间分辨率）和IJandsat TM卫星遥感影像（多光谱，30 m空间分辨率）生成，比例尺为1：50 000，土地利用类型共有32类，研究使用的是种植农业空间数据层。

2.1.2数据预处理 首先使用NASA提供的MRT（MODISReprojection Tools）软件对下载的MODIS反射率数据进行拼接、投影转换和格式转换。然后利用江苏省行政边界数据对转换后的影像进行裁剪处理，使用土地利用数据中的种植农业空间数据层进行掩膜处理。最终得到覆盖江苏省农作物种植区的反射率数据，生成研究区的NDVI序列数据：

农田的NDVI值应大于0，然而个别农田像元的NDVI值由于冬季或早春时积雪或其他噪声的影响出现了负值，将这些不合理的值用0代替。为了提高数据存储效率，将浮点型的NDVI数据扩大100倍，保存为无符号的短整型二进制文件。

每8 d的MOD09Q1反射率数据是基于每天的反射率数据、采用最大值合成法（MVC）以消除云和气溶胶影响而生成的影像，但是传感器观测角度、太阳高度角、水蒸气等引起的噪声仍然存在，会导致计算的NDVI出现不合理的波动。为此采用Savizky-Golay（简称S-G）滤波法对计算的NDVI序列数据进行3次平滑处理。参照梁守真等的成功经验，本研究选取的滤波窗口大小为7，平滑多项式的阶数为2。

2.2非对称性高斯函数和双Logistic函数拟合

为了重建冬小麦NDVI的变化动态，分别利用非对称性高斯函数和双Logistic函数对经过平滑处理的NDVI时序数据进行拟合。拟合时首先寻找冬小麦NDVI时序曲线上的极大值、极小值，然后用局部模型函数拟合极大值、极小值之间的数据。局部模型函数的通用表达式为：

2.3物候期提取

在分别利用非对称高斯函数和双Logistic函数对NDVI数据进行拟合后，采用动态阈值法提取冬小麦的返青期、抽穗期、成熟期（时间范围覆盖33～169 d）。根据司文才等的研究，本试验选择20%的阈值确定返青期、成熟期，即在NDVI的上升阶段，当NDVI值达到最大值的20%时对应的时间为冬小麦的返青期开始时间；在NDVI的下降阶段，当NDVI值减小到最大值的20%时对应的时间为冬小麦的成熟期开始时间。根据鹿琳琳等对抽穗期开始的定义，确定NDVI达到最大的时间点为抽穗期开始时间。

2.4物候期监测结果的验证

基于遥感数据的区域物候期监测结果的实地验证是目前物候研究的关键问题。本研究利用江苏省冬小麦主产区内的淮安市盱眙县、淮安市淮安区、盐城市滨海县、沭阳县、徐州市铜山区、连云港市赣榆区6个农业气象站的2009—2010年冬小麦生育期观测数据对冬小麦物候期的遥感提取结果进行验证。假设农业气象站的观测数据为真值，分别采用最大误差、最小误差、平均误差及均方根误差（root mean squared error，RMSE）为指标表征物候期提取的精度。其中均方根误差（RMSE）的计算公式为：

3结果与分析

3.1 S-G滤波结果

随机选取盐城市射阳县境内的1个耕地像元（行列号分别为593、1372）以查看S-G滤波的效果（图2）。可以看出，滤波前NDVI时序曲线存在明显的短期锯齿状波动，经S-G滤波后的NDVI时序曲线基本上保持了原有曲线的形态，并且变得较为平滑，极低的离异值得到很好的纠正。据调查，该像元的农作物为一年两熟制，夏熟为冬小麦，秋熟为水稻。滤波后的NDVI时序曲线呈明显的双峰形态，2个波峰分别对应于冬小麦、水稻的抽穗开花期。在冬小麦生长过程中，返青期之后随着气温的升高，冬小麦生长加速，叶片的叶绿素含量上升，对红光的吸收和近红外光的反射增强，NDVI值随之相应上升，到抽穗期达到峰值；抽穗期过后，随着冬小麦逐渐成熟，叶片开始变黄，叶绿素含量降低，对红光波段的反射率开始回升，近红外波段的反射率逐渐降低，NDVI值随之相应减少，直到冬小麦成熟收获。

但S-G滤波也具有局限性：滤波窗口的大小和平滑多项式阶数的选择均需要根据不同的条件人为进行设定，主观性影响较大。边金虎等研究表明，如果滤波窗口大小设置得偏小，容易产生大量的冗余数据，不易获取数据集的长期趋势；相反则容易遗漏一些细节所描述的正确信息。梁守真等研究表明，在一般情况下，平滑多项式阶数越小，处理后曲线越平滑，但拟合效果差。阶数越高，拟合效果越好，但有可能出现过拟合现象，而且在个别数据点会有振荡。因此，在进行S-G滤波处理时，滤波窗口的大小及平滑多项式的阶数是保证NDVI数据拟合准确性的关键。

图3是分别采用非对称性高斯函数、双Logistic函数拟合后的上述耕地像元NDVI时序曲线的对比效果，可以看出：2条拟合曲线均与滤波后NDVI值保持了相同的季节趋势；除到达波谷和第2个波峰阶段外，2条曲线上的NDVI值基本保持一致；在波谷和第2个波峰周围，利用双Logistic函数拟合的曲线斜率变化要比非对称高斯函数大，并且用前者拟合的NDVI值在波峰处比后者大，在波谷处比后者小。

3.2冬小麦关键物候期的空间分布

图4、图5、图6分别为采用非对称高斯函数（图4-a、图5-a、图6-a）和双Logistic函数（图4-b、图5-b、图6-b）拟合并利用动态阈值法得到的2010年江苏省冬小麦返青期开始时间、抽穗期开始时间、成熟期开始时间的空间分布。总体来看，基于2种拟合方法提取的物候期结果相似：返青期普遍开始于第49天之前，即2月18日之前；抽穗期普遍开始于第105天至第113天，即4月15—23日；成熟期普遍开始于第145天至第161天，即5月25日至6月10日。在空間分布上，研究区冬小麦返青期开始时间并没有随着纬度的变化而发生规律性的改变，大部分地区冬小麦开始返青的时间较集中，但在阜宁县、建湖县的大部分地区以及滨海县西南部、海门市东部地区返青期开始的时间较其他冬小麦区晚。抽穗期开始时间和成熟期开始时间的空间分布总体上表现出从南到北逐渐延迟的趋势：苏南地区抽穗期普遍开始于第105天之前，成熟期普遍开始于第153天之前；苏中地区抽穗期普遍开始于第113天之前，成熟期普遍开始于第145天至第161天；苏北地区抽穗期普遍开始于第113天之后，成熟期普遍开始于第153天之后。但是抽穗期、成熟期开始最晚的地区却没有集中在研究区的北部，在响水县、滨海县的大部分地区以及阜宁县北部、海门市东部和启东市北部地区，冬小麦抽穗期较其他地区晚；在泗洪县北部、阜宁县西南部、楚州区东部以及涟水县西北部地区，冬小麦成熟期开始得比其他地区晚。

3.3冬小麦关键物候期结果验证

表1、图7是采用非对称高斯函数、双Logistic函数分别拟合NDVI数据进而采用动态阈值法提取的不同站点物候期与观测数据的比较，可以看出遥感提取的物候期与观测的差普遍小于8 d。基于非对称高斯函数拟合提取的物候期误差超过8 d的样本数为：返青期1个（16.7%），抽穗期2个（33.3%），成熟期1个（16.7%）；基于双Logistic函数拟合提取的物候期误差超过8 d的样本数为：返青期1个（16.7%），抽穗期为1个（16.7%），成熟期为1个（16.7%）。结果表明，利用MODIS NDVI数据提取江苏省冬小麦物候期是可行的。

表2是基于MODIS NDVI数据提取的冬小麦生育期与观测数据比较的统计。可以看出，基于双Logistic函数拟合提取的冬小麦生育期精度稍好于基于非对称高斯函数拟合提取的结果，抽穗期的误差大于返青期、成熟期。利用双Logistic函数拟合方法提取的返青期、抽穗期、成熟期与观测数据比较的均方根误差分别为5.5、9.4、7.5 d，平均误差分别为4.5、7.8、5.5 d；利用另一种拟合方法提取的返青期、抽穗期、成熟期均方根误差分别为6.1、9.5、7.8 d，平均误差分别为5.3、8.5、6.3 d。

综上所述，利用非对称高斯函数、双Logistic函数拟合NDVI时序数据并采用动态阈值法进行判别可以有效提取研究区的冬小麦生育期，但提取的农作物物候期与站点观测数据尚存在一定的差异。主要原因有以下几个方面：本研究采用的是MODIS反射率的8 d合成产品，合成数据虽然在一定程度上保证了数据的质量，但明显降低了数据的时间分辩率，会导致物候提取的一定误差；此外，本研究使用的MODIS数据空间分辨率为250 m，遥感提取的物候期是像元内的平均值，而观测数据代表的是站点附近的物候，遥感数据与观测数据空间尺度差异也会影响遥感提取的物候与观测数据的一致性。当然，观测数据的误差也会影响遥感提取的物候的验证结果，将实地观测的和遥感提取的物候信息相结合是一种有效的物候提取途径。

4结论

快速、准确获得区域冬小麦物候信息，对于实施田问精细管理以及产量预测等具有重要意义。本研究基于MODISNDVI时间序列数据，分别利用非对称高斯函数、双Logistic函数進行拟合，并采用动态阈值法提取了江苏省2010年冬小麦物候信息。主要得出以下结论：

（1）S-G滤波处理可以有效地纠正时序数据中的奇异值，使NDVI时序曲线在保持原有形状的基础上变得较为平滑。

（2）非对称高斯函数、双Logistic函数2种拟合方法都可以使NDVI时序曲线更符合作物的生长过程。验证结果表明，就江苏省的冬小麦物候提取而言，双Logistic函数拟合方法略优于非对称高斯函数拟合方法。

（3）基于MODIS NDVI提取的2010年江苏省冬小麦返青期普遍开始于第49天之前，抽穗期普遍开始于第105天至第113天，成熟期普遍开始于第145天至第161天。抽穗期和成熟期的空间分布总体上表现出从南到北逐渐延迟的趋势，但返青期的空间分布特征不明显。同时也存在着一些不足，研究中采用动态阈值法提取物候期需要人为设定阈值，并且没有考虑阈值的空间变化，会影响物候提取的精度；由于观测数据的限制，仅使用了淮北地区6个点的数据对物候提取结果进行了验证，在其他地区，物候提取的误差尚不清楚；影响冬小麦物候的因素有很多，比如气候、地形、冬小麦品种、耕作制度差异等，它们的综合作用使其存在明显的空间变化。因此在今后的研究中，要综合考虑这些因素的影响，实现物候期提取的更细致化。