基于遥感和积温的冬小麦生育期提取方法

黄健熙 赵剑桥 汪雪淼 解智琨 卓 文 黄 然

(1.中国农业大学土地科学与技术学院， 北京 100083； 2.农业农村部农业灾害遥感重点实验室， 北京 100083；3.中国农业大学信息与电气工程学院， 北京 100083)

0 引言

生育期对农作物生长发育的动态监测、田间精细管理具有重要意义。生育期的准确提取有利于对作物的时空年际变化作出合理分析，为有效监测作物生长提供有力依据，进而反映气候变化对作物生长的影响[1-2]，并有利于估产模型的改进[3-8]。冬小麦是中国主要的粮食作物之一，研究冬小麦生育期的提取方法，精确监测冬小麦的关键生长阶段，对其产量预测具有重要意义[9]。遥感技术具有时效高、范围宽、成本低和时间序列连续等优点，能反映地面植被季节性生长发育的过程及其年际变化等特点，可为监测农作物生育期提供新的技术手段[10]。

已有较多学者利用Savitzky-Golay滤波(S-G滤波)平滑时间序列遥感数据，提取农作物生育期[11-12]。但在上述研究中，平滑NDVI时间序列时直接用了S-G滤波，使NDVI的值总是处于周围极大值和极小值之间。然而云雾和气溶胶的影响导致NDVI值偏低，因此时间序列上突降的点都应该作为噪声滤除，使用改进的S-G上包络线滤波能获得更高质量的时间序列[13]。

利用平滑模型函数拟合时间序列遥感数据及其产品，进而提取生育期，是近几年发展较快的一种方法，平滑模型函数包括Logistic函数法、非对称高斯函数法和谐波函数法[14]。SAKAMOTO等[15]利用MODIS/Terra数据，使用小波变换、傅里叶变换两种方法重构增强植被指数(Enhanced vegetation index，EVI)时序曲线。李铮等[16]以东北三省为研究区域,使用非对称性高斯函数拟合法平滑MODIS、CYCLOPES和GLASS叶面积指数时间序列,利用动态阈值法提取水稻的主要生育期。JONSSON等[17]提出基于非线性最小二乘拟合的非对称性高斯函数拟合AVHRR NDVI时序数据的方法，用于描绘地面植被的季节性生长和衰退曲线，并确定生育期参数。侯学会等[18]基于SPOT VGT NDVI数据，用5种方法提取冬小麦返青期，分析遥感监测结果与实测数据的均方根误差，iNDVI-Logistic提取误差为12.91 d，Logistic函数法提取误差为13.04 d，阈值法提取误差为17.48 d，导数法和DNA法提取误差大于35 d。

目前，很多研究综合使用上述两种时间序列遥感数据的处理方法，获得了较好效果[19-21]。但是前人成果主要集中于冬小麦一个或两个显著生育期(返青期和抽穗期)提取的研究，对于拔节期和开花期的研究很少，这是因为拔节期和开花期在LAI时间序列曲线上没有明显特征，因而，需要引入辅助数据。有效积温是作物基点温度之上日平均气温的积累[22]，有效积温与植物的生长速度和生育阶段有直接联系，是衡量作物生长发育过程热量条件的重要指示因子[23]。CHU等[11]利用MODIS数据提取了冬小麦返青期和抽穗期，发现积温每降低10℃·d，返青期延后4～5 d(R2=0.816,p<0.001)，抽穗期延后1～2 d(R2=0.401,p<0.001)。说明冬小麦生育期与积温具有显著的相关关系。胡乔玲等[24]在Logistic曲线拟合NDVI并提取返青期的基础上，结合积温进行拔节期推算研究，监测的冬小麦拔节期开始时间与观测值平均误差为4.3 d，最大误差为8 d。

MODIS LAI数据具有覆盖范围广、高时间分辨率的特征，在气候变化监测[25]、净初级生产力评估[25]、作物生育期监测[10]、作物产量预测[26-27]等方面有很广泛的应用。本文综合运用S-G上包络线滤波、Logistic曲线拟合、有效积温等，结合MODIS LAI数据和地面观测数据，提取并验证大范围冬小麦关键生育期，以实现大区域上冬小麦返青期、拔节期、抽穗期、开花期4个关键生育期的提取，并采用地面观测生育期数据定量评价提取精度。

1 研究区域与数据获取

1.1 研究区域

为了得到普适性的研究结果，选取河北、河南、山东三省为研究区域(图1)，进行冬小麦的生育期预测研究。该区域位于31°23′～42°40′N，110°21′～122°42′E，地处黄淮海平原，温带季风气候，夏季降水集中，雨热同期，冬春干旱少雨。年降水量500～900 mm，年均温11～14℃。

图1 研究区及农气站点Fig.1 Study area and agrometeorological stations

1.2 数据获取

1.2.1遥感数据

遥感数据使用MODIS LAI标准产品中的MCD15A3H陆地4级数据产品，该产品的时间分辨率为4 d，空间分辨率为500 m(https:∥ladsweb.modaps.eosdis.nasa.gov/search/)。本文使用的MODIS LAI产品时间为2015年1—7月，轨道号为H26V4、H26V5、H27V4、H27V5。该数据已经进行了几何校正和辐射校正,本文根据研究区域对该数据进行了拼接、裁剪。

1.2.2气象数据

气象驱动数据来源为中国区域高时空分辨率地面气象要素驱动数据[28-29](http:∥westdc.westgis.ac.cn/data/7a35329c-c53f-4267-aa07-e0037d913a21)，其时间分辨率为3 h，水平空间分辨率0.1°。本文使用的时间范围为2012—2015年，经过数据预处理，将其转换为366或365个波段的栅格文件，一个波段代表一天的日平均气温。

1.2.3观测数据

地面观测数据来自河北、河南、山东三省2012—2015年的农气站点观测记录，包含作物生育状况观测记录和土壤水分观测记录。其中作物生育状况观测记录包括台站编号、作物品种、栽培方式、冬小麦各生育阶段日期、生长状况、生长高度、生长密度、产量、产量因素和产量构成，以及主要田间管理措施等，在农气站点上还有关键生育期叶面积指数和生物量等观测值。本文使用了该系列数据中来自研究区域的64个农气站点的生育期数据。其中，2012—2014年的观测数据用于计算返青-拔节、抽穗-开花2个阶段冬小麦所需的平均积温，2015年的观测数据用于检验各生育期的提取精度。

由于本研究的空间跨度较大，返青期、拔节期、抽穗期、开花期几个关键生育期在不同地面观测站点呈现明显的南北差异，时间差别较大。如表1所示，从南到北各生育期的日期大致呈现逐渐推迟的规律，该现象符合从南到北辐射和降雨的空间分布规律。

表1 2015年冬小麦生育期观测数据Tab.1 Observed data of winter wheat growth stages in 2015 DOY

表1是对冬小麦生育期观测数据的统计，单位为一年中的天数顺序(Day of year, DOY)。4个关键生育期观测数据最小值均出现在河南省，最大值均出现在河北省。河南省整体生育期靠前，山东省居中，河北省普遍较后。拔节期的南北差异较返青期明显，河北省拔节期基本都在第100天之后，只有最南端的一个站点拔节期在第100天之前，为第87天；河南省主要在第90天之前，只有最北边的两个站点拔节期在第90天之后，分别为第94天和第104天；山东省拔节期分布在第85～100天，其中最南边的站点对应第85天，最北边的站点对应第100天。

抽穗期和开花期对温度的响应更加敏感，因此这2个生育期的南北差异更明显。河北省抽穗期主要在第120天之后，只有离河南省最近的站点在第120天之前，为第116天；河南省则均匀分布在第98～120天；山东省抽穗期分布在第112～126天。河北省开花期均匀分布在第120～131天，最北边的站点对应第131天，最南边对应第120天；河南省的开花期均匀分布在第104～126天；山东省分布在第119～131天，最北边的站点对应第131天，最南边对应第119天。

2 研究方法

2.1 技术路线

本研究利用MODIS数据重构LAI时间序列，根据时间序列的特征提取返青期和抽穗期。在此基础上利用气象数据，计算积温阈值并提取冬小麦的开花期、拔节期。通过LAI时间序列得到抽穗期日期后，从抽穗期开始计算有效积温，有效积温一旦达到近三年抽穗期到开花期有效积温的平均值，则当天为开花期。同理，在返青期的基础上得到拔节期。技术路线见图2。

图2 冬小麦生育期提取技术路线图Fig.2 Flow chart of winter wheat growth stages extraction

2.2 LAI时间序列重构

2.2.1S-G上包络线滤波

S-G滤波最早于1964年由SAVITZKY和GOLAY[30]提出。它可以理解为一种权重滑动平均滤波，其权重取决于在一个滤波窗口范围内做多项式最小二乘拟合的多项式次数[31]。S-G滤波过程为

(1)

其中

N=2m+1

i——时间索引

Ci——第i个LAI值的滤波系数

Yj+i——时间j处第i个LAI的原始值

m——窗口半径

如果只使用S-G滤波，滤波后时间序列上每一点的值是窗口内各点值的均值，不能将窗口内的极大值包含在内，存在部分突降无法消除的问题。因此，使用了CHEN等[13]提出的S-G上包络线滤波法来重构冬小麦LAI时间序列。该方法的处理步骤为：

(1)对原始LAI时间序列进行S-G滤波，分别存储滤波后和滤波前的时间序列。

(2)对比步骤(1)中存储的2个时间序列，得到新的时间序列，并将其作为原始序列。

(2)

式中O——原始的LAI值

L——滤波后的LAI值

T——迭代次数

采用S-G上包络线滤波算法对遥感数据进行去噪处理。由于云污染和水汽等的影响，遥感图像存在数据质量偏低甚至缺失的情况，因此滤波前的LAI时间序列曲线噪声较多，存在尖锐拐点，不够平滑，难以直接用于提取冬小麦的生育期；由图3可以看出，滤波后得到了外包络原始序列的平滑曲线，消除了原始数据的云污染和缺失数据造成的误差，更准确地反映了冬小麦的生长变化规律，便于之后的Logistic曲线拟合。

图3 2015年原始MODIS LAI时序曲线与S-G上包络线滤波结果对比Fig.3 Comparison of MODIS LAI and S-G upper-envelope LAI in 2015

2.2.2Logistic曲线

Logistic模型是由比利时数学兼生物学家VERHULST于1838年首先提出的。其特点是一开始缓慢增长，而在以后的某一范围内迅速增长，到达一定限度后，增长再度缓慢下来[32]，公式为

(3)

式中t——LAI时间序列中的时间

y(t)——t时间的LAI值

a、b、c、d——拟合参数

对拟合后的Logistic曲线方程求二阶导数，得到

(4)

原始MODIS LAI时序曲线、S-G上包络线滤波结果、Logistic曲线拟合结果及曲线二阶导数的对比见图4。

图4 2015年MODIS LAI、S-G上包络线滤波LAI、Logistic拟合LAI、Logistic曲线二阶导数对比Fig.4 Comparison of MODIS LAI, S-G upper-envelope LAI, fitting Logistic LAI and second derivative of fitting Logistic LAI in 2015

2.3 冬小麦种植区域提取

在研究区域内提取2015年LAI时间序列曲线受冬小麦控制的像元，根据冬小麦生长期内LAI时间序列的特征，给出提取较纯像素的条件：由于冬小麦在抽穗期生长旺盛，LAI值较高，因此要求第31波段(121 DOY)的LAI值大于1。由于冬小麦LAI于抽穗期达到峰值后持续下降，且本研究的空间跨度较大，各区域抽穗期DOY差异大，因此要求在第22～36波段(85～141 DOY)内，LAI有一个极大值，且该极大值大于第36波段(141 DOY)的LAI值。以上2个限制条件，能滤除MODIS LAI时序曲线不符合冬小麦发育情况的像素。最后得到的研究区域内2015年冬小麦分布如图5所示。用地面采样的方法验证得到冬小麦种植区域提取总体精度为90.75%，Kappa系数为0.86。

图5 2015年研究区冬小麦分布图Fig.5 Winter wheat distribution map in study area in 2015

在研究区域内，具有有效地面观测数据的农气站点共64个。但有些站点附近(3×3的栅格)混合像元问题严重，以农气站点周围3×3栅格中至少有5个像元种植冬小麦为筛选条件，最后保留了35个站点，可用于冬小麦生育期提取后的验证。

2.4 冬小麦关键生育期提取

返青期是指早春麦田半数以上的麦苗心叶长度达到1～2 cm的时期。冬季麦苗停止生长，在该时期突然开始生长，LAI时间序列曲线表现为突然上升。抽穗是禾谷类作物发育完全的穗，随着茎秆的伸长而伸出顶部叶的现象。全田50%植株抽穗为抽穗期，抽穗期处于冬小麦营养生长和生殖生长并进阶段，且LAI在该时期前后达到最大值。

在冬小麦生育期内，其LAI变化曲线近似于抛物线。进入抽穗期时，冬小麦的长势较好，叶片的生长状况在整个生育期中属于最好时期，冬小麦LAI在整个生育期中处于峰值[33]。因此，从滤波后的LAI时间序列中提取LAI最大值所对应的天数顺序，即为当年冬小麦的抽穗期。

从返青期到抽穗期，冬小麦的LAI呈单调递增，在抽穗期达到极大值后，从抽穗期到开花期，处于下降状态。返青期到抽穗期这一增长过程，由Logisitc曲线较为准确地拟合出来。从拟合的曲线中提取出二阶导数的最大值，最大值所对应的天数顺序即为2015年冬小麦的返青期。

根据2012—2014年的气温格网数据以及农气站点记录的返青期、拔节期，计算出返青期-拔节期的平均积温。将此平均积温作为阈值，结合2015年提取的返青期、气温格网数据，提取2015年的拔节期。同理，根据2015年提取的抽穗期，基于3年平均积温，提取2015年的开花期。

2.5 验证方法

根据各农气站点的经纬度，提取遥感影像中对应的像元，以该像元为中心，扩展至3×3像元区域，即1.5 km×1.5 km的区域。取该9个像元中的冬小麦区域生育期的平均值作为提取值，与对应的农气站点观测的生育期对比。假设农气站点的观测值为真值，分别采用最大误差、最小误差、平均误差及均方根误差(Root mean square error, RMSE)作为冬小麦生育期提取精度的评价指标。

MODIS LAI的时间分辨率为4 d，再考虑到混合像元等因素的影响，当生育期RMSE小于6 d时，认为该生育期提取精度较高。

3 结果分析

3.1 冬小麦关键生育期提取结果与分析

获得2015年研究区冬小麦生育期提取结果如图6所示。提取结果具有显著的空间变异性，与观测数据的时空变异规律基本吻合。分析具有地面观测数据的提取值与观测值，可得：返青期提取值的范围是29～91 DOY，观测值的范围是39～69 DOY。如图7a所示，剔除2个异常样本，返青期误差在0～5 d内的样本数为17个(51.5%)，误差在6～10 d内的样本数为6个(18.2%)，误差超过10 d的样本数为10个(30.3%)。

图6 2015年冬小麦生育期提取结果Fig.6 Results of extracted winter wheat growth stages in 2015

LAI数据源的空间分辨率为500 m，因此混合像元是引起各生育期误差的主要原因之一。由图7a可看出，混合像元的效应造成提取的返青期日期偏后，即提取日期大于观测日期。各生育期中，返青期的提取对混合像元非常敏感。由于夏季作物LAI快速增长的时期晚于冬小麦，LAI时序曲线二阶导数最大的天数顺序在像元内夏季作物的影响下产生延迟。此外，MODIS LAI能够较好地反映自然植被与林地的LAI，然而对于农作物而言，往往低于其实测LAI[34]。MODIS LAI的这一特性，降低了对返青期的提取精度。

另一影响返青期提取的因素是Logistic曲线拟合的精度。返青期在LAI上体现的特征非常细微，本文采用求二阶导数的方法提取返青期。因此，如果Logistic曲线不能准确刻画LAI快速增长期变化趋势的特性，就会导致提取的返青期产生一定误差。

根据提取的2015年返青期，利用2012—2014年3年历史积温平均值作为积温阈值，提取出当年拔节期。拔节期提取值的范围是52～133 DOY，观测值的范围是69～108 DOY。如图7b所示，剔除2个异常样本，拔节期误差在0～4 d内的样本数为20个(60.6%)，误差在5～8 d内的样本数为9个(27.3%)，误差超过8 d的样本数为4个(12.1%)。

图7 2015年生育期提取值与观测值的对比Fig.7 Comparison of extracted and observed growth stages in 2015

与站点观测值相比，返青期提取日期以延迟居多，导致开始计算积温的日期也产生延迟。研究区域内冬小麦返青期一般处于2月中下旬，在返青期实测日期—返青期推测日期这段时间内，各天的日平均气温多数都小于基点温度0℃，即这段时间的有效积温接近于0 ℃·d。因此，返青期提取值的延迟对积温计算的影响较小，不会使拔节期的提取产生较大误差。

抽穗期提取值的范围是85～137 DOY，观测值的范围是98～127 DOY。如图7c所示，抽穗期误差在0～4 d内的样本数为20个(57.14%)，误差在5～8 d内的样本数为11个(31.43%)，误差超过8 d的样本数为4个(11.43%)。

抽穗期提取受混合像元的影响较小。混合像元对冬小麦生育期提取的影响主要来自于各种不同生长周期的夏季作物，其中树木的LAI远高于冬小麦，使冬小麦LAI的峰值显著升高；大棚、园艺作物、蔬菜作物等使冬小麦LAI峰值下降，曲线在达到峰值之后的下降趋势不明显。

但是混合像元仅对LAI的数值产生影响，对LAI峰值出现的时间影响较小，即峰值出现的时间主要由冬小麦决定。因此，抽穗期提取结果与地面观测结果基本一致。

根据提取的2015年抽穗期，利用2012—2014年3年历史积温平均值作为积温阈值，提取出当年开花期。开花期提取值的范围是87～150 DOY，观测值的范围是104～134 DOY。如图7d所示，开花期误差在0～4 d内的样本数为21个(60.0%)，误差在5～8 d内的样本数为11个(31.4%)，误差超过8 d的样本数为3个(8.6%)。经过对比，开花期提取值与站点观测值吻合情况良好。

此外，提取的生育期还受到以下因素的影响：原始数据存在误差，提取的生育期需要精确到1 d，而MCD15A3H的时间分辨率为4 d。对空间尺度的差异而言，农气站点的观测数据是点上数据，而通过遥感数据提取的生育期是3×3栅格内的平均值。农气站点人工观测生育期数据存在误差。

3.2 精度评价

由图7可知，提取的冬小麦生育期时间与观测值之间的误差有提前和延迟的现象，也存在基本一致的情况，其中多数返青期提取值存在延迟情况，其他生育期误差正负分布较为均衡。表2为与农气站点观测数据对比后，基于遥感与气象数据提取的冬小麦生育期的精度评价。由表2可知，提取的返青期、拔节期、抽穗期、开花期与农气站点观测数据比较，其平均误差分别为7.4、4.5、4.4、3.8 d，均方根误差分别为9.5、5.5、5.2、4.9 d。

表2 2015年冬小麦生育期提取精度评价Tab.2 Evaluation of extraction accuracy of winter wheat stages in 2015 d

与已有研究相比，本研究中返青期平均误差与均方根误差显著小于WANG等[21]的研究结果，均方根误差显著小于侯学会等[18]用5种方法研究得出的均方根误差；拔节期平均误差与胡乔玲等[24]研究结论比较一致；而目前关于冬小麦开花期时间的研究尚无数据加以对比验证。总体来看，本研究对冬小麦生育期的提取精度较高，达到了较以往的生育期提取方法更符合实际的提取结果。

4 结论

(1)以河北、河南、山东三省为研究区域，MCD15A3H产品的空间分辨率为500 m、时间分辨率为4 d，通过S-G上包络线滤波重构的LAI 时间序列，结合时序曲线特征和积温方法，提取出较为准确的冬小麦关键生育期。

(2)S-G上包络线滤波方法可以将滤波窗口内的极大值包含在内，解决了直接使用S-G滤波时部分突降无法消除的问题，更准确地反映了冬小麦的生长变化情况。结果分析表明，各生育期开始时间由南到北逐渐推迟，空间变异性符合实际的辐射和降雨的空间分布规律。提取的生育期与农气站点观测日期较为吻合，返青期的平均误差在8 d之内，拔节期、抽穗期、开花期的平均误差都在5 d之内。这是由于求二阶导数的方法对混合像元及Logistic函数拟合准确度敏感，返青期的提取结果出现延后现象，而拔节期、抽穗期、开花期的提取精度较高。

(3)由于冬小麦种植区提取存在误差，研究中MODIS LAI时序曲线可能包含非作物信息，且积温模型所使用的气象插值产品精度有待验证，对生育期提取精度有一定影响。其次，农气站点在研究区内分布不均，疏密程度和代表性不同，可能影响结果验证的准确性。此外，拔节期、开花期的积温计算采用历史年份积温均值，其建模精度受年际间气候和作物品种差异程度影响。