基于历史气象资料和WOFOST模型的区域产量集合预报

马鸿元 黄健熙,2 黄 海 张晓东,2 朱德海,2

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

0 引言

作物生长模型是采用数学模型方法描述作物光合、呼吸、蒸腾、营养等机理过程，以特定时间步长来动态模拟作物生长发育期间的生理生化参数的数学模型，能较为准确地描述光、温、水、肥等因子以及田间栽培和管理措施对作物生长和发育的影响[1-2]。自20世纪60年代作物模型研究领域创始以来，作物模型逐渐成为农学、土壤学、植物生理学、气象学等学科的重要研究工具，尤其在作物产量模拟方面，作物模型能够动态、定量地模拟作物生长发育和产量形成，因此在产量预报[3]、农业气象服务[4]、产量的气候变化响应模式[5]等领域发挥了重要作用。如欧盟委员会联合研究中心(JRC)很早就将作物模型应用于作物监测和产量预报，在农业遥感监测项目(MARS)中建立了专门的作物生长监测系统(CGMS)，定时发布基于作物模型的产量月报[6]；美国[7]、印度[8]、巴西[9]等农业大国也都很早开展了基于作物模型的产量预报，我国李振海[10]、帅细强等[11]、陈上等[12]以预报或历史气象资料为基础初步开展了DSSAT、Oyrza2000及CERES-Maize模型在作物产量预报中的应用研究。

随着应用领域的不断深入，目前作物模型已从最早的田间尺度模拟扩展到了大区域尺度的模拟[13-14]，在气候变化情景模拟和遥感数据同化方面应用尤为广泛，但由于作物模型大多按天积分运算，需要输入全生育期的完整气象要素才能模拟得到最终的产量结果。基于作物模型的产量预报方法主要还是受制于预报期内的气象资料的获取，因此根据不同的研究需要，形成了各种气象资料的获取方法。

在月度、季度预报中，一些研究直接利用区域气候模型生成的数据作为作物模型的输入，进而预测农作物产量[15]，由于生成的气象数据存在一定的误差，还需要对气候模型的输出进行纠偏和订正[16-17]，同时低分辨率的气候模型模拟的逐日天气数据与气象站点实测数据有较大的尺度差异，直接使用气候模型产生的逐日气象数据时不可避免地存在时间和空间尺度不匹配的问题[18]。

为了实现气候模型与作物模型的高效耦合，必须将气候/天气模型的数据进行降尺度，常用的方法一是统计降尺度[19]和随机天气发生器法等，其中天气发生器在作物模型中应用较广，多项研究通过天气生成器对各气候模型的月平均数据降尺度生成逐日气象要素，作为作物模型的驱动数据[20-21]；二是通过动力学模型降尺度[22]，该方法利用复杂的动力学方程把气候模型的模拟结果降尺度到更小的网格和时间尺度上。动力降尺度模型考虑了局地的地形因子，从而能更好地反映出局地气候特征，对极端天气状况有较好的模拟效果，但此方法计算复杂费时，应用在区域化的作物模型中成本较高。

而在中短期预报中，数值天气模式则可以和作物模型相结合，BAIGORRIA等[16]使用COAPS的区域模型生成了20个数值天气预报集合成员，耦合陆面模型后输入到CERES-Maize模型，其中缺失的辐射量要素使用了天气发生器生成，最终使用主成分回归分析进行了产量预报。LEE等[23]则开发了一种气象数据服务端，利用站点观测数据和数值预报数据直接生成ORYZA 2000作物模型的气象驱动数据，大幅提高了作物模型的部署效率。DUMONT等[24-25]使用STICS作物模型进行了大量的预报气象驱动数据比较工作，验证了使用历史平均资料和天气发生器等不同气象驱动策略对当年小麦产量预报的影响，结果显示使用两种策略没有明显差异，使用天气发生器精度稍好，但使用历史平均资料的计算效率则更高。

我国大部地区为典型的季风气候，年际变率大，且季节分布不均[26]，因此，基于气象资料多年平均值作为预报期的驱动数据并不能很好地体现气象要素的波动性，而且灾害性的异常天气会被平滑，这对模拟潜在的产量损失十分不利；而如果将历史同期所有年份的气象要素都作为集合成员进行驱动模拟，则可以重现本地历史天气情景下的产量[27]，通过集合预报的技术手段得出产量预报的不确定性，将以往单一数值的预报提升为概率预报[28-29]。

本文以多年气象资料为依托，以Python平台的WOFOST模型(Python crop simulation environment，PCSE)为基础，构建使用历史资料作为预报期气象驱动数据的系统框架，区域化分布式模拟研究区内多年冬小麦产量，并系统分析产量集合的统计特征，结合农业气象观测资料与产量统计数据分析历史集合预报方法的预报能力和精度，以期为区域化产量概率预报提供科学依据。

1 数据与方法

1.1 研究区概况

选取河北省保定市、衡水市冬小麦主产区为研究区(图1)，该区地处黄淮海平原，属暖温带半湿润季风气候，多年平均气温3～15℃，年降水量500～800 mm，四季分明，日照充足。冬小麦在本区粮食组成中占有较大比重。研究获取了衡水市2009年29个观测点的冬小麦田间观测数据[30]，主要包括播种与收获日期、物候期、叶面积指数(Leaf area index，LAI)和产量等，试验点的分布充分考虑了该地区冬小麦的生长分布状况，具有很好的代表性。

图1 研究区示意图Fig.1 Location of study region

1.2 气象资料

所用的气象驱动数据是中国科学院青藏高原研究所发布的中国区域高时空分辨率地面气象要素驱动数据集(China meteorological forcing dataset)01.06.0014版本[31-32]，该驱动数据集包括近地面气温、地表气压、近地面空气比湿、近地面全风速、降水率和向下长、短波辐射7个要素，空间分辨率为0.1°，时间分辨率为3 h。该数据集是以国际上现有的Princeton再分析资料、GLDAS资料、GEWEX-SRB辐射资料以及TRMM降水资料为背景场，融合了中国气象局常规气象观测数据制作而成的。经过前期处理，将这些要素转换为与作物模型量纲、格式一致的气象输入文件，包括日照辐射量、日最高温、日最低温、清晨水汽压、风速和降水6个要素。

1.3 WOFOST模型

WOFOST(World food studies)模型是由荷兰瓦赫宁大学和世界粮食研究中心共同开发的一种机理性作物生长模型[33]，该模型通过利用光截获量和CO2同化量作为生长驱动过程，利用作物物候发育过程来描述作物生长，模拟的生长过程主要包括作物的碳同化、呼吸、蒸腾、干物质的分配等。模型主要有两种应用模式，潜在模式下作物生长仅仅由温度和太阳辐射决定，没有考虑其他生长限制因素；水分限制模式下作物生长受水分限制。WOFOST模型是一种通用模型，通过定制不同的参数可以模拟不同种类的作物，具有很强的适用性，是大范围环境条件下的各种作物模拟的有利工具。模型最初由Fortran语言开发，而今WOFOST研究组已经开发了适合Python平台的版本PCSE，通过配置丰富的第三方库，大大加强了模型应用的灵活性。

1.4 产量预报框架

图2 基于历史气象资料的作物模型产量集合预报技术流程Fig.2 Flow chart of research technology route

传统的作物模型运行过程只是输入数据到输出数据的单一映射，无法提供模型模拟的不确定性信息，但作物模型本身，以及输入参数、气象驱动存在着不确定性。如果以集合的思想将作物模型放入一个外部框架内，用不同年份的气象驱动模拟的结果作为集合成员，可以使模拟结果的集合代表其概率分布，最终将单一数值的模拟输出转换为概率分布。基于上述考虑，本文引入多年(1979—2015年)历史气象资料，开发出以从某一预报时刻模拟该年作物生长和产量状况的概率预报的集合预报框架。该框架以Python平台下WOFOST(PCSE)模型为基础，作物品种参数和区域化由本课题组历史标定成果[34-39]确定，通过外部程序框架的搭建，完成多年气象输入文件生成、预报时刻气象集合成员文件生成、分布式PCSE模拟、集合预报成员后处理等主要功能，从而实现区域化集合预报作物长势和产量的主要功能(图2)。利用该预报框架开展产量集合预报时，首先进行前期数据准备，预先收集作物品种参数、田间管理信息及播种至预报时刻的气象资料，然后根据预报时刻生成模型的气象输入文件集合，该集合的每个成员文件由播种日至预报时刻的实际气象和预报时刻至成熟日的历史资料拼接而成，利用该输入集合进行区域化的分布式模拟，最终得到模拟产量的集合并进行相应的后处理过程。

1.5 验证方法与误差分析

以1979—2015年气象资料为历史集合进行集合预报，以2009年29个观测站的田间实测产量为观测值，以2009年每月1日为预报起始时刻，同一日期下历史集合模拟结果的平均值作为预报值，观测值与预报值的差值即为预报误差。

集合预报给出的预报是一个概率分布或者是一个包含多个成员的集合，若需要对集合预报进行确定性验证，一般是对集合成员的平均值或集合中位数进行验证[40]。集合成员的平均值计算公式为

(1)

Mi——集合成员值

n——集合成员数量

当有了成员集合后，还可进行某一产量水平的概率预报，假设集合中所有成员的权值相等，计算某一事件的成员数与集合成员总数的比例，即为该事件发生的概率，计算公式为

(2)

式中pi——模拟结果中某一事件(如产量高于7 000 kg/hm2)出现与否的标识，当该事件发生时pi=1，否则为0

平均差会由于差值的正负抵消，因此本文采用平均绝对误差(Mean absolute error，MAE)表示误差，另外还将均方根误差(Root mean square error，RMSE)与皮尔逊相关系数(Pearson correlation coefficient，PCC)作为精度评价指标，相关计算式为

(3)

(4)

(5)

式中m——验证样本数量

i——验证样本序号

Fi——第i实测点模型预报值

Oi——第i实测点观测值

其中，MAE越小说明预报越精确，PCC越接近于1说明预报模型的预报能力越强。通过对每个观测点进行1979—2008年的集合模拟，对预报的误差进行分布拟合，最后统计3种指标作为预报能力。

2 结果分析

2.1 单点预报分析

图3为深州市2009年3月1日(图3a)、4月1日(图3b)、5月1日(图3c)、6月1日(图3d)的叶面积与籽粒干物质量预报结果，其中不同颜色的线条表示不同集合成员的全生育期模拟值，左上方蓝色直方图分别表示6月1日的LAI与最终产量的集合成员分布频率。红色虚线标记的是该点的实测产量(6 786 kg/hm2)，蓝色竖线指示预报日期，该时刻之前的气象驱动资料为2009年实际气象数据，该天及之后的气象驱动为1979—2008年的历史集合成员。结果显示预报开始后集合成员的模拟结果逐渐出现分化，LAI的集合分布在5.0～6.5之间，产量的集合分布在6 000～8 300 kg/hm2之间，不同年份的模拟结果有较大的差异，但随着预报时刻的不断更新，LAI的集合成员模拟值变化幅度减小，而且基于历史气象的模拟集合其变化区间能够较好地涵盖实测值，说明历史气象数据的集合化表达可以用来部分表达模型模拟的不确定性，在小麦拔节后期即可得到对小麦长势、LAI较为准确的集合预报。

模型模拟的产量即存贮器官干物质质量(Total dry weight of storage organs，TWSO)的集合平均值与实测值相差始终在500 kg/hm2左右，而直到5月初产量模拟结果仍有很大变异性(图4)，这是因为冬小麦生长前期是以营养生长为主，同化的干物质主要分配到叶、根、茎中；生育中期是小麦营养生长和生殖生长并进时期；在小麦生长后期同化物已停止向营养器官中分配，主要分配给贮藏器官；同时根、茎、叶等营养器官中的干物质会有部分转移到穗部形成产量，因此在抽穗后小麦籽粒开始形成、积累大量干物质并进行器官间的养分转移直至成熟，所以在冬小麦的产量形成中抽穗—灌浆的生殖生长阶段十分关键，这一时期作物遭受气象胁迫作用会十分明显地影响产量，这也说明在后续实际业务预报中，在冬小麦生长早期集合成员模拟的变异性较大，此时进行预报效果较差，而最终产量预报结果对抽穗—灌浆阶段的气象预报最为敏感，从生育后期即5月开始进行华北地区产量预报有一定的预报价值。这段时间气象数据的不确定性将决定产量预报结果的成败，因此在生长早期进行业务产量预报时还需要结合中长期的数值天气预报产品。

图3 深州市观测点4个预报时期的LAI和TWSO集合模拟值与实测值对比Fig.3 Comparison of simulated and measured values of LAI and TWSO in four forecast periods of Shenzhou City

图4 2009年深州市实测点不同起报时刻的冬小麦产量滚动集合预报及其箱形图Fig.4 Rolling ensemble forecasts of winter wheat yields in 2009 and their box plots at different forecasting times in Shenzhou City

2.2 区域预报分析

通过分布式区域化模拟预报，得到研究区内区域化的产量集合预报结果。为了对产量预报结果进行验证，本文以预报的集合中值产量和实测点产量进行对比验证，可以看出随着预报日期的推进，研究区南部和北部的产量差异加大，南部预报的产量有所提高，北部产量逐渐降低，说明2009年该区域南部整体天气条件较好，实际产量优于历史条件下的模拟(图5)，而北部区域天气条件差于历史同期。从该研究区产量大于7 500 kg/hm2的区域化概率预报(图6)，可以看出在2月初开始预报时，根据历史同期气象资料，大部分地区尤其是保定地区的产量高于7 500 kg/hm2的概率低于50%，然而随着预报时刻的推进，绝大多数地区产量高于7 500 kg/hm2的概率大于70%，说明2009年气象条件适宜冬小麦生长，天气条件优于历史同期。

2.3 预报性能定量分析

图5 2009年8个预报时期下区域化集合产量预报(集合中位数)Fig.5 Regionalized yield ensemble forecast (median of ensemble) for eight forecast periods in 2009

图6 2009年8个预报时期下产量高于7 500 kg/hm2的概率预报Fig.6 Probability forecast for event of production above 7 500 kg/hm2 on eight forecast dates in 2009

为了进一步评估集合预报的性能和误差，本文对所有田间实测点的产量集合预报进行了统计，根据实测点的实测产量计算16个预报时期集合中位数的验证指标MAE、RMSE、PCC。16个预报时期的MAE为458～578 kg/hm2, RMSE在598～758 kg/hm2之间，PCC在0.503～0.563之间(图7)，集合平均数的情况与之类似，但由于篇幅所限未列出。总体而言基于历史气象的集合预报对于本年的产量有一定的预报能力，相对误差在6.6%左右，集合中位数与实测值相关性也较好，说明历史气象集合的变异性很好地表达了气象输入参数的不确定性，有助于概率预报和数据等应用。但随着预报期的推进集合中位数的误差指标收敛较慢，直到5月中旬左右预报的性能达到最大值，说明气象数据的不确定性仍然不足以完全表现实际产量的不确定性，尤其是潜在模式下模型对于水分胁迫和土壤特性的影响没有得到表达，而且受制于气象驱动数据空间分辨率的限制，对于精细尺度的产量无法仅凭气象数据进行预报。

图7 衡水市实测点不同起报时刻下产量滚动预报的集合中位数验证结果Fig.7 Verification results of rolling ensemble forecast median of 29 observation points in Hengshui City at different forecasting times

3 讨论

气象数据是作物模型的主要驱动数据，利用作物模型模拟作物生长需要完整描述整个生育期的气象数据，缺乏预报期气象驱动数据是制约作物模型应用于产量预报的主要瓶颈，因而目前作物模型的应用主要局限于历史产量估测、单点模拟验证等，而作物模型机理明确、可扩展性强的优势没有得到发挥。气象历史资料对于本地的气候模式和变异性有很好的描述能力，因此可以使用历史资料作为预报期的输入资料。本文采用潜在模式下的WOFOST模型作为预报模型，对2009年衡水保定地区的冬小麦产量进行了模拟预报，并使用集合预报的方法进行了后处理，结果发现基于历史气象的集合预报可以进行本年产量的预报，相对误差在6.6%左右，预报精度较好，表明该方法可以作为缺乏预报期气象资料时的补充手段，而且集合成员的分布可以表现产量的不确定性，可以根据集合成员提供产量水平的概率预报。

本文从3月开始进行滚动预报，3月初模型即有一定预报能力(PCC为0.505, MAE 为504 kg/hm2)，说明在本地的历史多年气象资料集合可以代表当年的气象要素变率；而且预报误差在后期逐渐减小，说明随着气象输入数据不断被更新为实测值，气象数据的不确定在下降；但产量集合的收敛过程较慢，直到灌浆中后期达到最高预报精度，这可能与模型模拟的空间尺度较粗而且模型其他参数未能区域化有关，另外在潜在模式下作物生长的模拟未考虑水分等胁迫作用的影响，区域上的产量差异仅仅来自于气象输入的差异，因而使得在空间上产量的分布梯度较不明显。

由于作物品种、土壤特性、管理措施对于产量的变异性有很强的影响，各种输入参数的变异性共同构成了产量的变异性，为了进一步提高产量区域化预报的精度，需要对敏感参数进行区域化，然而目前作物模型区域化调参受制于数据缺乏因而难度较大，将来借助农业气象观测、土壤数据库等数据进行参数区域化是进行模型产量预报的重要前期过程。另外卫星遥感数据是栅格类型的观测数据源，具有覆盖面积广、空间连续、周期观测、相对成本低的优势，通过遥感观测和作物模型的数据同化可以获得大区域范围的模型参数分布[41]，因此数据同化将会在模型产量预报方面有很大的潜力。

作物生长模型是一个非线性的复杂系统，通过模拟可以给出产量和其他生长参数的确定性模拟值，而集合预报的技术手段可以得到预报值不确定性信息，本文采用集合预报的思想进行了作物模型模拟，集合成员的生成方式是历史多年的气象输入，而作物模型的不确定性来源除了气象输入外，还有模型本身、输入参数的不确定性，因此多模式的作物模型、作物参数扰动都可以用来生成预报集合[42-43]，这样的超级集合更有助于表达整个系统的不确定性，为提高模拟精度、概率预报和数据同化[44]提供了一条很好的途径。和数值天气预报的技术发展历史相似，从单值的模拟拓展为概率预报是产量预报的发展方向，因此使用集合预报的思想进行作物模型模拟具有很好的应用前景。

本文初步探讨集合预报手段在作物模型模拟中的技术流程和应用前景，因此仅通过较小规模研究区一年的模拟预报对技术方案的可行性和误差进行了实验，而本文的技术方案还需要更多验证，多年大区域的验证实验还有待进行，而且解决预报期天气输入的根本方法在于结合数值天气预报，对于集合数值预报和作物模型的结合需要更深一步的研究。此外，在气候变化的大背景下，近年来的极端、异常天气事件增多，需要采用其他方法得到预报期更高精度的气象输入集合，如历史相似年份替代、天气发生器等技术手段。

4 结论

(1)使用多年气象资料作为驱动建立预报集合，集合的平均值与实测产量PCC最高为0.563，MAE最低为458 kg/hm2,表明历史气象集合对于本年度的作物产量具有一定的预报能力，能够作为预报期的气象数据补充手段，而集合的平均值或中位数能否代表当年实际气候变异下的模拟产量，取决于在历史气象集合中当年气象条件的一般性。

(2)抽穗至灌浆时期是作物产量形成的关键期，冬小麦产量对于这一阶段的气象条件十分敏感，在抽穗后开始产量预报有最高的预报精度。因此为进一步提高基于作物模型的产量预报精度，需要准确的抽穗至灌浆时期的气象参数，作物模型与中长期数值预报结合是提高预报时效的必然要求。

(3)现阶段研究在潜在模式下模拟预报，产量的变异性只与气象数据中的温度和辐射量有关，因此驱动数据的低空间分辨率和输入参数无区域差异是制约预报精度的关键因素，使得模拟产量的空间梯度较小，预报精度难以进一步提高，为了解决输入参数的估计问题，与遥感观测进行数据同化是可能的解决途径。