基于MODIS与WOFOST模型同化的区域冬小麦成熟期预测

2019-10-10黄健熙高欣然马鸿元朱德海

农业机械学报 2019年9期

黄健熙高欣然黄海马鸿元苏伟朱德海

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

0 引言

作物收获期对作物产量与品质影响显著，收获过早或过晚都不利于作物的丰产增收[1]。对于冬小麦而言，其主产区多采用冬小麦-夏玉米轮作的种植制度，收获过晚将直接影响夏玉米的及时耕种。此外，随着农村土地制度改革的不断推进，传统的以家庭为单位的小农经济正向大规模的机械化农业转型[2]，优化收割次序、提前调度农机，以实现作物的准时收获是机械化农业发展带来的新问题，因此，在实际生产中迫切需要对作物成熟期进行大区域精准预报。

在农学意义上来说，成熟期是指作物的营养存贮器官成长到可收获的阶段。影响成熟期变化的因素较为多样，形成了不同的研究方法。文献[3-4]着重分析气象因素(如温度、光周期、降水等)与作物成熟期的相关性，通过在单点尺度上建立气象因子与成熟期的统计模型估测实验年份的成熟期。段金省[5]利用多年的气象数据和某年份作物的成熟期数据建立查找表，认为气候条件越相近的年份生长的作物成熟期也越接近，并以此为理论基础估测成熟期。基于气象因子的估测方法易于操作，在使用气象预报数据或区域气候模型[6]的情况下能够较快进行成熟期的预报，但是该类统计方法难以在实验区域外推广，且无法应对年际间种植品种和田间管理措施等因素的变化。

另一类方法是通过分析生育期内作物生化参量的变化特征提取成熟期，这类方法多以遥感数据作为数据源，主要分为两种技术路线：①分析生育期内植被指数的时间序列曲线，计算指数时序曲线的曲率、波峰等特征，提取成熟期[7-8]，主要方法有：阈值法[9]、Logistic曲线拟合法[10]和滑动平均法[11]等。②监测生育期内作物水分[12]、叶绿素[13]、氮素[14]、干物质[15]等指标的变化，使用回归分析等统计学方法建立成熟期估测模型，进而实时监测作物是否达到成熟所需指标。基于遥感数据监测作生化参量提取成熟期的方法能够在大区域上提取作物的成熟期，但是无法进行预报。此外，由于云、雨天气的影响，作物生育期内不可避免地会遇到遥感数据缺失或数据质量较差的问题，将给时间序列曲线和生化指标监测带来误差，影响成熟期提取精度，而且这类方法基于回归分析构建模型，对于品种调整和气候变化缺乏适应性，难以在实验区域外推广应用。

作物成熟期受到气候条件、作物品种特性和耕种时间的共同影响[16]，现有研究方法大多基于历史数据对成熟期建模，不能很好地反映当年作物品种及耕种时间的变化。如果根据当年作物实际生长发育情况动态调整模型中与作物品种相关的参数，将能极大提高模型的适应性。因此，引入作物生长模型，并将遥感数据作为观测数据，构建代价函数优化作物模型参数，不仅可以利用遥感数据大范围、实时性强的优点，而且又能以气象预报数据驱动机理过程模型，实现区域尺度的成熟期预测。

1 研究区与数据

1.1 研究区概况

选择河南省冬小麦种植区作为研究区，其覆盖范围为31.38°～36.37°N, 110.35°～116.64°E，总面积约5.425×106hm2(中国国家统计局，2015年)，地处黄淮海平原，属暖温带季风气候区，年平均降水量约650 mm，年平均日照时数约2 300 h，四季分明，雨热同期，光照充足。冬小麦是该地区种植的主要粮食作物之一，通常在10—11月播种，2—3月返青，5—6月成熟。由于气象条件、播种时间和种植品种的差异，研究区内冬小麦的成熟期在空间上和年际间也存在一定差异，适用于进行冬小麦成熟期的预测与验证。本文所用冬小麦种植区通过随机森林方法提取得到[17]，空间分辨率为500 m，在研究区内包含35个国家级农业气象观测站(图1)。

图1 研究区与农业气象站点分布Fig.1 Map of study area and agrometeorological sites

1.2 研究数据

1.2.1遥感数据

遥感数据来自2017—2018年MODIS LAI产品MCD15A3H (https:∥ladsweb.modaps.eosdis.nasa.gov/)，其时间分辨率为4 d，空间分辨率为500 m，通过使用ArcGIS工具对影像进行投影转换、镶嵌和裁剪等处理，输出覆盖整个研究区的MODIS LAI数据。

1.2.2气象数据

气象数据来自2017—2018年欧洲中期天气预报中心提供的气象要素数据集以及集合预报格点数据集TIGGE(https:∥apps.ecmwf.int/datasets/data/tigge/levtype=sfc/type=cf/)，包括近地面气温、近地面气压、露点温度、近地面全风速、太阳辐射、降水量等要素，其时间分辨率为6 h，空间分辨率为0.25°，其中预报数据的时限为16 d。对原始气象数据进行预处理得到时间分辨率为1 d，空间分辨率为500 m的数据，包括：日最高温、日最低温、日照辐射量、水汽压、平均风速与降水量6个要素，目前在多个领域得到广泛应用[18-20]。超出预报数据时限的气象数据参照本课题组之前的研究方法[21]，用历史气象数据代替。

1.2.3作物与土壤数据

作物和土壤数据分别来自农业气象站记录的《作物生育状况观测记录年报表》和《土壤水分观测记录年报表》，作物数据具体包括站点的区站号、台站名称、作物名称、品种名称、品种类型、作物关键生育期、台站经纬度、台站海拔、株高、密度、田间管理措施等；土壤数据具体包括土壤质量含水率、土壤相对湿度、土壤容重、田间持水量、凋萎湿度等。部分站点还包括关键生育期叶面积、植株干质量等数据。

我国工程建设监理制度的推行虽然走过了20多年，但有序竞争的市场环境还没有形成。2010年，中国水利工程协会组织开展了首批水利建设市场主体信用评价工作，但采取的是工程施工和监理市场主体自愿申请参加的方式，仅有137家水利建设市场主体自愿申报进行信用评级。其中监理市场主体只有67家自愿申报进行信用评级，信用评级是在市场主体申报材料的基础上，由全国性行业自律组织——中国水利工程协会进行初审、复审和组织评审委员会终审，信用评价缺乏强制性和约束性，信用评价结果缺乏公信力，实际的信用监管效果还不够理想，优胜劣汰的市场准入和清出机制没有形成，制约监理活动良性发展的因素依然存在。

2 研究方法

2.1 WOFOST模型的本地化

WOFOST模型是由瓦赫宁根大学和瓦赫宁根研究中心共同开发的一种机理性作物生长模型，能够模拟生育期内作物对气象和其他环境因子的反应，实现不同类型作物发育阶段的量化、生物量和作物产量的模拟[22]，目前已经广泛应用于指导农业生产和模拟作物长势等多个方面[23-26]。为了使模型能够准确地模拟研究区内冬小麦的生长发育过程，首先需要对模型参数进行标定，实现WOFOST模型的本地化，使用郑州农业气象试验站的历史观测数据或查阅文献[27-28]来确定模型的参数取值，部分参数校准值见表1。

2.2 MODIS LAI数据预处理

本文将研究区的MODIS LAI数据按照时间顺序叠加，在各像元上生成对应的LAI时间序列曲线，使用Savitzky-Golay(S-G)滤波方法剔除MODIS LAI数据中由云、水汽、气溶胶等带来的噪声。由于500 m分辨率影像受到混合像元因素的影响，且MODIS LAI数据的尺度低于WOFOST模型模拟LAI的尺度，为了减小系统误差给同化造成的影响，采用归一化的方法处理LAI数据，即把MODIS LAI和模型模拟LAI均映射到0～1内，在保留MODIS LAI趋势信息的情况下使两种LAI数据源尺度相同，公式为

表1 WOFOST模型中主要参数校准值Tab.1 Calibrated values of main parameters in WOFOST

注：表中DVS表示作物的发育阶段。

LAIn=(LAIo-LAImin)/(LAImax-LAImin)

(1)

式中LAIo、LAIn——归一化前、后叶面积指数

LAImax、LAImin——该时间序列上叶面积指数的最大值和最小值

2.3 待优化参数的选择

由文献[29-31]可得，对作物成熟期影响较大的作物参数有IDEM(出苗日期)、TBASEM(出苗下限温度)、TSUMEM(从播种到出苗期有效积温)、TSUM1(出苗期至开花期有效积温)、TSUM2(开花期至成熟期有效积温)等。由于TBASEM和TSUMEM均反映作物从播种期到出苗期的特性，而本文选择LAI作为耦合变量，且遥感LAI数据只能反映作物出苗之后的生长状况，因此无法有效地优化播种期相关参数。分析历史数据发现，在研究区TSUM1、TSUM2、IDEM 3个参数没有表现出明显的空间规律，且年际间变化较大，因此将三者设为待优化参数，根据文献[28-29]及历史实测数据规定其初始值与上下限如表2所示。

2.4 SCE-UA算法与代价函数构建

优化算法采用复合形混合演化算法(Shuffled complex evolution-University of Arizona, SCE-UA)。该算法是DUAN等[32-33]在1993年提出的引入了复杂形分割、混合思想的全局优化算法，在样本空间的搜索效率、计算速率和全局搜索最优的能力上表现突出且对待优化参数初始值不敏感[34]，为遥感与作物生长模型同化的实际应用提供了可行的方法。

表2 待优化参数初始值及取值范围Tab.2 Initial value and range of parameters for optimization

由于之前的研究发现，与顺序同化算法相比，SCE-UA算法更易受到观测数据误差的影响，当观测LAI偏低时同化后的LAI偏低较严重[35]，因此在构建代价函数时，先将遥感观测LAI和WOFOST模拟LAI归一化，得到代价函数

(2)

代价函数J为时间序列遥感观测LAI与模型模拟LAI的误差平方和。SCE-UA算法通过全局搜索，当满足邻近5个代价函数值均达到设定阈值或者迭代次数达到设定次数时停止运算并得到最优参数集，将所得最优参数集连同其他参数输入WOFOST模型，实现模型模拟冬小麦生长发育过程和成熟期的预测。

3 结果与分析

3.1 WOFOST模型的标定

在本研究中，本地化的WOFOST模型是准确预测冬小麦成熟期的必要条件。耦合变量LAI是表征植被冠层结构的基本参数，在WOFOST模型中是直接影响作物的光合、呼吸等重要过程的状态变量[36]。因此，采用2016年郑州和泛区农业气象站全生育期实测LAI数据对标定后的WOFOST模型模拟结果进行验证。图2是WOFOST模型模拟冬小麦LAI时间序列与实测LAI对比。

图2 WOFOST模型模拟冬小麦LAI时间序列与实测LAI的对比Fig.2 Comparisons of WOFOST simulated time series LAI and field-measured LAI

冬小麦LAI在出苗期至越冬期前增长较缓慢，越冬期停止增长，返青期开始增长迅速，在开花期前达到最大值，开花期至成熟期逐渐降低，成熟期之后迅速降低。从图2可以看出，标定后的WOFOST模型模拟的郑州站和泛区站的LAI时间序列曲线的趋势特征与实测值基本符合，经计算，模拟LAI和实测LAI间的RMSE为0.92 m2/m2，说明模拟值与实测值具有良好的一致性，WOFOST模型模拟的冬小麦LAI精度较高。此外，部分点的模拟产生了一些偏差，可能的原因有：①年际间冬小麦生育期的变化导致WOFOST模型用2018年的作物参数无法模拟出较为准确的开花期和成熟期，影响了LAI的模拟精度。②模型标定是以郑州站的观测数据为准，因此不一定适用于其他地区冬小麦的模拟。

3.2 单点同化结果分析

开花期和成熟期分别是冬小麦营养生长和生殖生长阶段的结束日期，是评价模型模拟发育期合理的关键指标。因此本文使用研究区内农业气象站的开花期和成熟期观测日期共同检验实验结果的精度，图3为2018年滑县站同化前后模型模拟LAI和MODIS LAI的时间序列曲线。

图3 滑县站同化前后模型输出LAI和MODIS LAI的时间序列曲线Fig.3 Changing curves of WOFOST simulated LAI, assimilated LAI and MODIS LAI in Huaxian station

从图3可以看出，滑县站同化后WOFOST模型模拟的出苗期、开花期与成熟期略晚于同化前。从LAI时序曲线的趋势分析，同化后LAI曲线的变化趋势更接近于MODIS LAI，表明WOFOST模型优化TSUM1等参数能够使其输出的LAI曲线趋势更接近遥感数据，而且并没有受到遥感数据尺度较低的影响。相对应地，同化后LAI明显高于同化前LAI和MODIS LAI。通过分析滑县气象数据和模型的逐日模拟结果，发现原因是同化前模型模拟的开花期在DOY110附近，由前文可知，冬小麦的快速生长时期发生在返青期至开花期之间，同化前返青期至开花期气温较低，光照不足，不利于冬小麦的快速生长；而同化后模拟的开花期较晚，为DOY122，该阶段气温回升，使同化后的冬小麦在返青期至开花期期间长势优于同化前。由此看出，准确地模拟开花期不仅是成熟期模拟的基础，更对WOFOST模型生物量、产量等多方面的模拟精度有很大的影响。

图4 开花期、成熟期的预测值与实测值的回归分析Fig.4 Regression analysis diagrams of predicted and measured anthesis and maturity date

在表3中，同化前开花期、成熟期取值范围分别为DOY97～108、DOY135～152，同化后两者取值范围变为DOY108～120、DOY141～156。可以看出，多数点同化前的开花期和成熟期均早于同化后，且同化后开花期、成熟期的时间变异性大于同化前。对比各站点同化前后参数值，发现部分开花、成熟较晚的站点同化前开花期、成熟期均提前，说明标定的作物参数并不符合这些站点的实际情况，导致样本内部差异较小。同化前开花期、成熟期误差较大也是因为在标定WOFOST模型时研究区内的TSUM1、TSUM2等作物参数是基于历史数据计算取值，没有考虑到年际间品种与播种时间等因素的变化，同化后误差明显减小，说明同化MODIS LAI和WOFOST模型的方法可以提高研究区冬小麦开花期、成熟期的预测精度。

表3 部分站点同化前后开花期、成熟期日期对比Tab.3 Comparison of simulated and assimilated anthesis and maturity date DOY

开花期、成熟期的预测值与实测值的回归分析见图4。从图4可以看出，同化后的开花期、成熟期与对应观测日期基本符合，所成散点的趋势线斜率均在1左右，说明同化后的开花期、成熟期与观测日期具有良好的一致性，且RMSE分别为2.10 d和2.48 d，预测精度较高。误差的主要来源有：①遥感数据空间分辨率较低，同化混合像元可能导致开花期模拟出现误差，随着模型的运行误差逐渐累积，进而影响成熟期的预测精度。②由于研究区内冬小麦品种和田间管理等方面存在差异，而大量参数在空间上的分布无法获得，只能采用标定值，因此模型运行时其他参数取值相同也会带来一定的误差；使用SCE-UA算法最小化代价函数的过程中，可能很难准确地找到最接近实际情况的全局最优参数集，因此优化后的参数集也可能存在一定误差。

3.3 开花期与成熟期空间分布

图5 研究区同化前后开花期空间分布对比Fig.5 Comparison of anthesis date distribution before and after assimilation in study area

图5、6分别为同化前后开花期、成熟期的空间分布图。从图5可见，由于只受气象输入数据的影响，同化前开花期模拟值呈较为明显的区块状分布，相邻的像元开花期时间十分接近。西北部和东部地区由于冬小麦播种时间较早且气候适宜，因此开花期最早，南部、中部次之，而河南省西部海拔较高，北部纬度较高，因此气温相对较低，开花期最晚。同化后的开花期总体趋势为从西南部向东北部逐渐推后，西部高海拔地区仍然最晚。同化前后开花期的空间分布变化主要表现在：①同化后研究区西南部开花期整体提前、研究区西北部和东部略有延后，原因是实际的IDEM、TSUM1数值与标定值相比发生了较大变化，导致同化后相应区域开花期普遍延后或提前。②同化后的开花期表现出了更详细的空间变异性，在表现出开花期整体趋势的前提下，可以看出相邻冬小麦像元的开花期存在一定差异，说明用MODIS LAI优化WOFOST模型的参数后可以在一定程度上表现出较细空间尺度上种植品种等方面的差异。

对比图5a和图6a可以看出，研究区中部、东部和西南部地区成熟期相对提前(位于研究区内的先后顺序)，说明这些地区当年度夏季气温较高，光照充足，因此生育期提早结束。从图6可以看出，同化后成熟期的空间分布除了增加空间细节之外，还可以看出，中北部地区成熟期大幅延后，研究区内成熟期整体趋势为从西南部向东北部逐渐延迟，主要原因是优化后的TSUM2(开花期至成熟期有效积温)大于标定值，更符合该地区的实际情况，提高了模型的预测精度。

图6 研究区同化前后成熟期空间分布对比Fig.6 Comparison of maturity date distribution before and after assimilation in study area

4 讨论

通过数据同化框架耦合遥感数据和作物生长模型，并结合气象预报数据实现区域尺度作物成熟期的预报，虽然总体精度较高，但仍存在一定的误差，误差产生的原因有：①预报延伸期气象驱动数据存在误差，由于5月16日(DOY136)后的气象数据是使用历史相似年份的同期气象数据替代的，因此与实际气象数据有一定偏差，以后可以尝试天气发生器[37-38]、区域气候模型[6]等多种中长期气象预报方法。②本方法仅优化了WOFOST模型的3个输入参数，其他参数取值在研究区内相同，在一定程度上也限制了模型的预测精度，可以考虑增加对LAI敏感性较高、在生育期不同阶段波动较大的参数作为待优化参数，如AMAXTB(叶片最大CO2同化速率)、SLATB(比叶面积)等。

冬小麦开花期的模拟精度会直接影响其成熟期、LAI模拟的精度，误差过大时甚至会影响LAI模拟的尺度，因此通过数据同化方法优化WOFOST模型参数进而提高开花期模拟精度对于成熟期预测具有重要意义。而本文使用的MODIS LAI数据存在混合像元效应且容易受到云的影响，使用滤波的方法平滑MODIS LAI也会导致一定的误差，因此为进一步提高成熟期预测精度，可以考虑增加高分辨率遥感数据或非光学遥感数据(如SAR数据)进行多遥感数据同化。此外，可以考虑将量化开花期模拟的不确定性量化，在开花期模拟集合的基础上进行同化，也可以加入作物参数扰动、多来源气象预报数据分别构建预报集合[39-40]，这些集合更有助于表征整个系统的不确定性，进而为决策者提供多维度的信息。