基于冬小麦农业气候分区的WOFOST模型参数标定

李 颖 赵国强 陈怀亮 余卫东 苏 伟 程耀达

1)(中国气象局·河南省农业气象保障与应用技术重点实验室，郑州 450003) 2)(河南省气象科学研究所，郑州 450003) 3)(河南省气象局，郑州 450003) 4)(哈尔滨市气象局,哈尔滨 150028) 5)(中国农业大学土地科学与技术学院，北京100083) 6)(郑州大学生态与环境学院,郑州 450001)

引 言

作物生长模型起步于20世纪60年代，其实质是用数学方法表达作物生长过程[1]。伴随作物生长理论的完善、系统科学和计算机技术的发展，作物生长模型经历了从理论走向实用的发展历程[2]。在世界各国对粮食生产的持续关注下，以应用为目的的作物生长模型迅速发展，既可用于现代化和科学化的作物种植管理，也可辅助政策的分析制定[3]，已成为农业研究最有力的工具之一[4]。随着作物生长模型向应用多元化方向发展，其区域化大范围应用的要求也日益迫切[5]，其中典型应用包括区域作物生长监测和产量预测，准确进行大面积作物生长监测和产量预测对加强作物生长调控和保障粮食安全具有重要意义[6-7]。

现有的作物模型，如WOFOST，CERES，APSIM和WheatSM等，都具有较强的机理性和实用性，已在小麦、水稻、玉米等作物的长势监测和产量评估中得到广泛应用[8-9]。通过对作物模型进行参数的严格标定，在单点作物生育期的模拟误差普遍小于3 d，产量的估算误差一般不超过15%[10-13]。然而，随着作物模型在区域上的升尺度应用，其区域模拟效果受到来自模型品种参数、外部输入数据等因素不确定性增大的影响，模型的区域模拟精度较难满足实际应用需求[14]。

为提高作物模型在区域应用的精度，必须针对特定区域对其进行本地化，即完成模型参数的本地化校准[15]。当区域面积较大时，气候条件差异较大,作物品种不同，需对区域进行合理分区，结合科学高效的优化算法对模型参数分区标定。本文以WOFOST模型用于河南省大范围冬小麦生长模拟与产量预测为目标，通过划分不同农业气候生态区，各区分别对模型参数进行标定，将农业气候学知识与科学高效的优化算法相结合，通过合理、高效地对研究区域分区标定，提高作物模型在区域应用的精度，为作物模型区域应用和模型参数调整优化提供科学依据。河南全省农业气候生态区划分采用K均值聚类算法，模型参数的敏感性分析选择Sobol全局敏感性分析方法，敏感参数的标定采用差分进化马尔科夫链蒙特卡洛方法[16]。

1 研究区和研究方法

1.1 研究区概况

本文以河南省为研究区。河南省位于中国中部，黄河中下游地区，西起太行山和豫西山地东麓，南至大别山北麓，东西长约580 km，南北跨约550 km，全省土地面积为1.67×105km2，空间范围在110°21′～116°39′E，31°23′～36°22′N。河南省位于中国冬小麦核心主产区，平原面积广阔，土壤肥沃，其小麦种植面积占全国小麦种植面积的1/4以上，年产量约占全国的26%，为我国小麦产量第一大省[17]。

1.2 WOFOST模型

WOFOST模型由荷兰瓦赫宁根大学与世界粮食研究中心共同研发[18]。该模型是一个动态解释性模型，以逐日为步长对作物生长发育进行模拟，可描述光合作用、呼吸作用、蒸腾作用、干物质分配等作物基本生理过程，并可描述这些过程受环境影响情况，可模拟水分限制条件下的产量、营养限制条件下的产量和潜在产量[19]。

WOFOST模型以逐日气象数据为驱动，基于积温理论建立作物生长发育模块，并考虑作物对光长的敏感性和叶片的生长及衰老过程，对叶面积生长分两个阶段进行描述。模型中干物质的积累和分配则是经过呼吸作用的消耗后被分配到作物的各个器官。作物在每个时刻获得的总干物质重量，被按照一定的比例系数分配给地上部分和地下部分，地上部分再被分配到贮存器官和茎、叶等器官。现有研究表明，WOFOST模型在河南省及其相邻省份适用良好[19-20]。

1.3 数据来源与处理

本研究中冬小麦农业气候分区数据来源于1981—2010年河南全省113个气象观测站冬小麦全生育期逐日气象数据，WOFOST模型运行所需气象数据来源于研究区内20个气象观测站2013—2019年的逐日观测，包括日最低气温、日最高气温、日平均水汽压、日平均瞬时风速和日降水量。由于WOFOST模型的气象驱动数据包括太阳辐射量，利用Hargreaves辐射公式[21]由大气温差计算得到(式(1))，将6种气象要素按照WOFOST作物模型需要的格式建立对应数据库文件。

(1)

式(1)中，Ra为天顶辐射(单位：M·J·m-2·d-1)，Tmax为最高气温(单位：℃)，Tmin为最低气温(单位：℃)，KRs为调整系数。在内陆地区，陆地占主体，气团不会受到广阔水体的强烈影响，KRs约为0.16[22]，河南属于内陆地区，取0.16。

农业气象观测数据包括作物生长数据、土壤数据和田间管理信息，其中作物生长数据主要包括冬小麦播种期、关键生育期、产量等；土壤数据包括土壤重量含水率、土壤相对湿度、降水量与灌溉量、地下水位深度、土壤水分总贮存量、土壤有效水分贮存量、田间持水量、土壤容重、凋萎湿度等；田间管理信息主要为灌溉时间和灌溉量、施肥时间和施肥量。本研究所需农业气象观测数据、土壤数据和田间管理信息来源于研究区内20个农业气象试验站和农业气象观测站2013—2019年观测数据以及2019年田间调查的叶面积指数，剔除异常值，整理为标准数据格式。

1.4 研究方法

针对WOFOST模型用于大范围冬小麦生长模拟与产量预测目的，首先根据30年气象观测数据，利用K均值聚类算法将研究区河南省划分为不同的农业气候生态区，在每一个区间利用Sobol全局敏感性分析方法筛选出WOFOST模型的敏感参数，采用差分进化马尔科夫链蒙特卡洛方法对模型的敏感参数进行自动标定，并利用实测数据对模型分区标定结果进行验证。

1.4.1 K均值聚类算法

K均值聚类算法属于聚类分析方法，具有快速收敛的优势，其基本思想是以距离为相似度进行分类。给定聚类数目K后，首先随机把数据集划分为K个簇，计算每个簇平均值得到K个初始聚类中心，计算每簇中每个样本与每个聚类中心的距离，将其重新分配到距离最近的聚类中心所在的簇，再计算每个新簇的中心。不断重复该过程，直至聚类中心不再发生变化，则聚类结束[23]。

1.4.2 Sobol方法

Sobol方法是一种典型的全局敏感性分析方法[24-25]，经验证是最为有效的模型参数敏感性分析方法之一[26]。Sobol方法基于方差分解的思想，将参数化模型[16]表示为

y=f(X，θ)。

(2)

式(2)中，y为目标函数；X为驱动数据；θ为参数向量。

将模型总方差分解为单个参数的影响和参数之间组合的影响：

(3)

式(3)中，Di为参数θi对目标函数y的一阶偏方差，Dij为参数θi和θj之间相互作用的二阶偏方差，k为参数的维数。

对式(3)进行归一化，通过偏方差和总方差之比计算参数敏感性指标。

一阶敏感性指标

(4)

二阶敏感性指标

(5)

总敏感性指标

SOi=∑S(i)。

(6)

其中，SOi表示包含所有参数θi的敏感性，即参数θi单独对目标函数y的影响程度及与其他参数相互作用对y影响程度的总和[27]。

1.4.3 差分进化马尔科夫链蒙特卡洛方法

马尔科夫链蒙特卡洛方法是一种贝叶斯后验采样方法，可以提供趋近于真实后验分布情况的样本序列[28]。该方法的基本原理见文献[29]，在求解随机变量的期望或随机事件的发生概率时，为了克服蒙特卡洛随机抽样所需样本量随维数增加而指数增大，导致计算速度过慢的问题，通过构造马尔科夫链的极限平稳分布情况模拟计算积分，即利用先验知识从建议概率分布抽样产生候选样本值，建立状态转移规则并判断状态是否转移，从而产生马尔科夫链，马尔科夫链的各个状态为随机变量的样本值，在进行足够多次的迭代后，马尔科夫链的转移概率将趋于平稳分布，得到目标的后验分布[30]。由于马尔科夫链蒙特卡洛方法收敛效率较低，本研究使用Braak[31]提出的差分进化马尔科夫链蒙特卡洛方法对待优化参数的后验分布进行估算。该方法解决了传统马尔科夫链蒙特卡洛方法中使用单一马尔科夫链，易发生局部收敛的问题，可以同时构建多条马尔科夫链估算参数的后验分布，并在不同链之间进行相互学习，大大提高了计算效率[16]。

2 结果与分析

2.1 河南省冬小麦农业气候区划

本研究采用K均值聚类算法，根据河南省气候特点和冬小麦生产品种布局将全省划分为不同的农业气候生态区[32-33]。依据气候条件对冬小麦全生育期影响情况，利用专家知识法，选取1981—2010年全省113个气象观测站冬小麦全生育期降水(单位：mm)、全生育期大于等于0 ℃积温(单位：℃·d)、全生育期日照时数(单位：h)、冬前积温(单位：℃·d)、3—4月日照时数(单位：h)、5月降水量(单位：mm)、3—4月雨日(单位：d)、5月平均气温日较差(单位：℃)、3—4月降水量(单位：mm)、5月日照时数(单位：h)和5月雨日(单位：d)为区划指标进行K均值聚类，结合冬小麦生产品种的布局情况，完成河南省冬小麦农业气候分区(图1)。

如图1所示，河南省冬小麦生产分为5个农业气候生态区，基本呈纬向分布。其中Ⅰ区分布在河南省北部，包括安阳、濮阳等地区，该区日照充足，但大于等于0℃积温较少，冬季温度较低，常有冻害发生。春旱机率较高，多数年份小麦生长受到一定水分胁迫。Ⅱ区主要分布在河南省中北部，包括新乡、焦作、郑州、开封、商丘、洛阳等地区，这些地区春季降水变率大，自然降水偏少，春季低温霜冻对小麦有一定影响。Ⅲ区以河南省中南部为主，包括漯河、平顶山、周口、驻马店中北部、南阳中东部等地区。该区正常年份降水量可以满足小麦需要，但春季连阴雨天气较多，光照不足。Ⅳ区主要包括河南省南部的信阳、驻马店南部和南阳西部等地区，该区小麦全生育期大于等于0℃积温在2500℃·d以上，冬前气温高，春季多雨，灌浆期间高温、多雨、日较差较小。Ⅴ区主要分布在河南省西部的三门峡等地，该区干旱少雨，小麦生育期内降水量不足300 mm，春季日照充足，小麦生育后期常出现干热风。该分区结果一方面可体现纬度对冬小麦种植品种和生长发育的重要影响，另一方面也反映出河南省复杂地形影响气候、耕作方式等，进而对冬小麦种植品种和生长发育造成影响。河南省北部、西部、南部由太行山、伏牛山、桐柏山和大别山沿省界环绕，西南部为南阳盆地、中东部为黄淮海平原，图1中的分区良好地吻合该地形特点，显示本研究对河南省冬小麦农业气候分区结果的合理性。

图1 河南省冬小麦农业气候分区Fig.1 Agro-climatic division of winter wheat in Henan Province

作物模型中的积温参数可通过气象观测数据准确计算，在完成河南省冬小麦农业气候分区后，利用30年气象资料计算各区小麦播种至出苗、出苗至开花和开花至成熟各生育阶段的积温(表1)。

表1 不同区域冬小麦积温(单位：℃·d)Table 1 Accumulated temperature parameters of winter wheat in different zones(unit:℃·d)

2.2 敏感参数分析

应用Python中的SALib灵敏度分析库，排除不适合自动标定的参数，对WOFOST模型的21种小麦生长发育参数进行敏感性分析，分析过程包括待分析参数输入、参数样本生成、模型模拟过程以及分析结果输出。本研究通过设置待分析参数围绕默认值变化的比例确定其取值范围(表2)。

使用Saltelli采样器对表2中的参数进行随机采样，设置采样数为10000，产生(21+2)×10000=230000组样本输入进WOFOST模型运转，输出冬小麦模拟产量，从而建立从输入参数空间到输出结果空间的映射关系。利用Sobol方法筛选出对冬小麦模拟产量变化影响较大的参数。参数的实际采样值由如下公式计算：

表2 待分析参数及其围绕默认值的变化比例Table 2 The proportion of changes of parameters to be analyzed around the default value

ae=a0×f。

(7)

式(7)中，ae表示参数实际采样值，f表示本次采样比率，a0表示该参数默认值。

对于某些默认值为0的参数，采用如下公式计算参数实际采样值：

(8)

式(8)中，amax表示该参数最大值，amin表示该参数最小值。

在河南省5个农业气候生态区分别进行21种作物参数的敏感性分析，每个分区选择5个观测站点，使用2013—2017年的气象数据驱动模型。经对比，各分区总敏感性指数与一阶敏感性指数和二阶敏感性指数的筛选结果基本一致。以总敏感性指数为代表，河南5个农业气候生态区的参数敏感性分析结果见图2。

由图2可知，尽管各区敏感性分析结果中待分析参数的总敏感性指数不同，参数的敏感程度排序也存在差异，但在5个农业生态区，达到0.05的显著性水平，且总敏感性指数大于0.01的参数筛选结果一致，均为FL，FR，FO，SP，SL05，AM00，AM10，AM13和SL00。具体地，总敏感性指数大于 0.05的参数相同，共4个：FL，FR，FO 和 SL05。总敏感性指数小于0.05且大于0.01 的参数相同，共5个(不同区顺序不同)：AM00，AM10，AM13，SL00和 SP。

图2 5个农业气候生态区参数的总敏感性指数Fig.2 Total sensitivity index of parameters in five agro-climatic ecological zones

2.3 参数标定与精度验证

对敏感性分析筛选出的9种敏感参数进行本地化标定，使用差分进化马尔科夫链蒙特卡洛方法对其后验分布进行估算。差分进化马尔科夫链蒙特卡洛种群(并行链)数目设置为5条，每次采样都进行一次链更新，并通过方差法计算各个参数当前的收敛指标R，通过R判断是否达到收敛。当存在某个参数的R>1.1时，认为链未收敛；当所有参数的R<1.1 时，认为链达到收敛。收敛后再次进行1000次采样，由此得到参数的后验样本和分布。差分进化马尔科夫链蒙特卡洛方法中使用的似然函数为

lgL=lgLlai+lgLyield。

(9)

续图2

式(9)中，L表示总似然函数。Llai表示叶面积指数的似然函数，Lyield表示产量的似然函数，具体如下：

lgLlai=-0.5(p-pobs)T·V-1(p-pobs)-

0.5K·lg(2π)-lg(detV)。

(10)

式(10)中，p和pobs分别表示不同生育期叶面积指数的模型模拟值和观测值，V表示叶面积指数观测值的协方差矩阵，K表示空间维数，即叶面积指数观测值数量。

(11)

式(11)中，z和zobs分别表示成熟期产量的模型模拟值和观测值，σ表示产量观测值的标准差。

对各区使用2013—2017年多站点实测数据联合标定，本研究中参数标定所用建模数据和验证数据见表3。5个农业生态区待优化参数的后验分布见表4～表8，可以看到，多种作物参数在各区之间差异明显，可能体现不同区域适宜种植品种的不同，这也表明在大区域进行分区调参的必要性。将河南省作为一个整体区域使用2013—2017年多站点实测数据联合标定，整个研究区使用同一套参数的后验分布见表9。

表3 参数标定所用建模数据和验证数据Table 3 Modeling data and verification data used for parameter calibration

9种参数后验分布的平均值、中值以及最大似然值优化结果(代价函数取最优解时的一组参数值)见表4～表9，分别代入WOFOST模型进行模拟，使用2018—2019年观测数据进行验证，模拟及验证结果见图3。由图3可知，分区标定后，在各子区模拟的叶面积指数与实际观测值均吻合较好，综合5个分区的模拟结果，使用分区调参后验平均值模拟关键生育期叶面积指数的总均方根误差为0.655，其中，模拟返青期叶面积指数的均方根误差为0.428，模拟拔节期叶面积指数的均方根误差为0.753，模拟抽穗期叶面积指数的均方根误差为0.796，模拟乳熟期叶面积指数的均方根误差为0.578。使用默认参数模拟关键生育期叶面积指数的总均方根误差为2.897，整个研究区使用同一套优化参数(后验平均值)模拟关键生育期叶面积指数的总均方根误差为1.277，对比可知，经分区调参，模型动态模拟冬小麦多个关键生育期叶面积指数的精度均明显提高，表明基于农业气候分区标定WOFOST模型可准确模拟作物生长过程。在估产精度方面，综合5个分区模拟结果，使用默认参数模拟产量的均方根误差为 2282.192 kg·hm-2，整个研究区使用同一套优化参数(后验平均值)模拟产量的均方根误差为 1311.303 kg·hm-2，使用分区调参后验平均值模拟产量的均方根误差为 672.016 kg·hm-2，使用分区调参后验中值模拟产量的均方根误差为 684.622 kg·hm-2，使用分区最大似然值参数优化结果模拟产量的均方根误差为 796.436 kg·hm-2。对比可知，使用本研究提出的基于冬小麦农业气候分区的差分进化马尔科夫链蒙特卡洛参数标定方法，显著优于默认参数以及不分区调参的产量模拟精度，其中以使用分区调参后验平均值模拟产量的精度最优，较使用默认参数时产量模拟误差减小70.55%，较全区使用同一套优化参数(平均值)时产量模拟误差减小48.75%。

表4 河南省Ⅰ区敏感参数的后验分布Table 4 The posteriori distribution of sensitive parameters in Zone Ⅰ of Henan Province

表5 河南省Ⅱ区敏感参数的后验分布Table 5 The posteriori distribution of sensitive parameters in Zone Ⅱ of Henan Province

表6 河南省Ⅲ区敏感参数的后验分布Table 6 The posteriori distribution of sensitive parameters in Zone Ⅲ of Henan Province

表7 河南省Ⅳ区敏感参数的后验分布Table 7 The posteriori distribution of sensitive parameters in Zone Ⅳ of Henan Province

表8 河南省Ⅴ区敏感参数的后验分布Table 8 The posteriori distribution of sensitive parameters in Zone Ⅴ of Henan Province

表9 河南省同一套敏感参数的后验分布Table 9 The posteriori distribution of sensitive parameters for the whole Henan Province

图3 2019年产量与叶面积指数验证结果Fig.3 Verification results of yield and leaf area index in 2019

续图3

3 结论与讨论

本研究采用K均值聚类算法根据1981—2010年气象数据将河南省划分为5个农业气候生态区，在各区分别计算作物模型中的积温参数，再根据2013—2017年的观测数据利用Sobol全局敏感性分析方法选择其他敏感性参数，随后采用差分进化马尔科夫链蒙特卡洛方法对敏感参数进行标定，并用2018—2019年的观测数据进行验证，得到结论如下：

1) 河南省5个农业气候生态区，整体呈纬向分布，良好吻合研究区的地形特点。冬小麦农业气候分区结果可为作物模型分区标定提供科学的分区依据。

2) 对5个农业气候生态区分区进行WOFOST模型参数的敏感性分析，结果表明：不同冬小麦品种筛选出的敏感参数具有较高一致性。

3) 在农业气候分区的基础上，使用差分进化马尔科夫链蒙特卡洛方法优化WOFOST模型参数，模型模拟冬小麦多个关键生育期叶面积指数的精度和产量预测精度均显著优于使用默认参数或整个研究区使用同一套优化参数时的模拟精度。

本研究以河南省为研究区，将农业气候学知识与科学高效的优化算法相结合，有效实现作物模型区域应用时模型参数的优化调整，但存在以下不足：一方面，本研究局限在省域范围内，若区域面积进一步扩大，冬小麦种植品种有更大变化时，尚需进一步研究模型的敏感参数是否发生变化以及模型输出的不确定性将如何变化等问题；另一方面，本研究使用差分进化马尔科夫链蒙特卡洛方法通过最优化代价函数自动为模型敏感参数赋值，不同参数在不同区间差异较大，有待根据不同区间冬小麦实际种植品种进行测量分析，进一步验证采用差分进化马尔科夫链蒙特卡洛方法为作物参数自动赋值的准确性。下一步用于整个黄淮海地区可进一步优化本研究提出的基于农业气候分区的参数标定方法，实现WOFOST模型在大区域应用的本地化工作。

此外，本研究以近30年的气象观测数据为依据对研究区进行分区，但各区内模型模拟效果仍受到来自模型品种参数、外部输入数据等因素的不确定性影响。为进一步提高作物模型的区域化应用精度，可采用网格单元运转模型，并在该过程中将遥感信息与作物模型耦合解决模型空间扩展过程中格点参数获取困难的问题。