园艺作物发育期和采收期模拟模型的最优模拟路径

程陈，李春，李文明，叶晨洋，王岳胜，赵承森，丁枫华，金志凤，冯利平，黎贞发※

（1. 丽水学院生态学院，丽水 323000；2. 中国农业大学资源与环境学院，北京 100193；3. 天津市气候中心，天津 300074；4. 丽水市气象局，丽水 323050；5. 浙江省气候中心，杭州 310020）

0 引言

作物模型以环境因子为驱动变量，运用计算机技术和物理数学方法，对作物生长发育和产量品质形成过程进行定量描述与预测，是一种面向作物生长发育过程、机理性强的数值模拟模型[1-2]。田间试验不仅需要收集大量农业资源数据，而且可能仍然无法提供有效管理决策，故作物模拟模型在提高作物生产力和经济效益等农业管理中发挥着越来越重要的作用。目前没有一个单独模型能确保在所有条件下比其他模型模拟效果更好[3-5]，因此需要探究模型的最优模拟路径，以提高模型的普适性和准确性。

根据发育期模拟模型算法原理的不同，主要分为4类建模方法：第1 种为温差法，昼夜温差缩短了果实发育周期，果实尺寸减小，开花数和果实数也减少，但是累积温差对发育期影响的定量研究较少[6-7]；第2 种为生长度日法，是简便实用的发育期预测方法，数据获取较为便利[8-9]；第3 种为生理发育时间法，综合考虑了日长和温度2个环境因子对作物发育的影响[10-11]；第4 种为辐热积法，综合考虑了辐射和温度2个环境因子对作物发育的影响[12-13]。在模型温度响应方面，有研究表示当高温或昼夜温差较大时，作物发育速率与温度存在非线性关系，会导致发育期模拟存在偏差[14-15]。在模型算法步长方面，随着观测仪器、计算机和通信技术的迅速发展，气象要素逐小时（甚至更高频次）采集和存储变得越来越容易，为精确模拟作物及其生长环境的定量关系提供了有力的数据支撑[16]，且已经应用于作物蒸散发校正方面[17]。PURCELL[18]比较了2 种线性模型小时步长和日步长在不同温度条件下热量累积的差异，还需进一步探究光照对不同步长模型稳定性及模拟精度的影响。在模型算法机理方面，WU等[19]分析了5 种常用小麦模型的不同算法，表示5 种模型均能较好地模拟作物发育期，但模拟精度有所差异。有学者利用贝叶斯方法提高了玉米[20]和水稻[21]模型的精确度。MARTINEZ等[22]结合AquaCrop 和MOPECO模型，提高了田间缺水状况的模拟精度。而针对不同发育期建模方法在不同类型园艺作物（叶菜、瓜果和花卉等）的精确度如何，还有待考证。关于采收期模拟模型研究中，黄瓜、番茄和茶叶均是无限生长型园艺作物[1,23]，采收期受到季节性生产能力和产量等因素影响，生产迫切需要探究采收期最优模拟路径，以期合理利用人物力资源。

本研究获取了9 a 58 茬分期播种试验观测数据，以气温和辐射为环境驱动量，构建基于4 类建模方法（温差法、积温法、生理发育时间法和辐热积法）的园艺作物发育期和采收期模拟模型，确定重要的模型参数，并以4 种方式（平均值、最值均值、中值和逐步回归）集成模拟结果，最终明确园艺作物发育期和采收期模拟模型的最优模拟路径，为园艺作物智慧生产管理和高效利用农业资源方面提供理论基础和技术支撑。

1 材料与方法

1.1 试验设计

黄瓜、番茄和芹菜试验地点设于天津市武清区农业科技创新基地园区（116°58'12''E，39°25'48''N，海拔8 m），郁金香试验地点设于北京市顺义区杨镇国际鲜花港日光温室基地（116°47'24''E，40°10'12''N，海拔38 m），芹菜、菠菜和香芹试验地点设于浙江省丽水市丽水学院玻璃温室基地（119°55'12''E，28°27'N，海拔61.8 m），茶叶试验地点设于浙江省丽水市松阳县新兴乡上安村（119°22'48''E，28°31'48''N，海拔166.2 m），表1 介绍了供试品种、试验地点、试验始末日期和数据介绍，每个播期设3个重复，采用随机区组设计。

表1 供试品种及试验设计Table 1 Crop varieties experimented and experimental design

1.2 数据获取

1.2.1 发育期数据

依据文献将园艺作物生长发育过程分为多个关键发育期，其中，黄瓜分为6个关键发育期（移栽期、伸蔓期、初花期、坐果期、采收初期和采收盛期）[1]，番茄分为5个关键发育期（移栽期、三叶期、初花期、坐果期和成熟期）[23]，芹菜分为4个关键发育期（移栽期、外叶生长期、心叶肥大期和枯萎初期）[13]，香芹分为6个关键发育期（播种期、出苗期、一真叶期、三真叶期、五真叶期和枯萎初期）[24]，菠菜分为5个关键发育期（播种期、出苗期、三真叶期、五真叶期和枯萎初期）[24]，郁金香分为5个关键发育期（移栽期、现蕾期、始花期、盛花期和末花期）[9]，茶叶的采收标准为“两叶一心”[24]。当大于或等于50%的作物群体出现某一发育期特征时，即为该作物群体进入了某一发育期，并记录该日期。季节性生产能力、茬口搭配、光热资源利用率和人为因素等导致实际种植中作物拉秧期的随机性较大，故本研究不进行该发育期的模拟。

1.2.2 气象数据

温室内常规观测选用小气候观测仪（CAWS2000 型，北京华云尚通科技有限公司）自动记录空气温湿度、CO2浓度及太阳总辐射等气象要素，采集频率为10 min/次，加密辅助观测选用小气候观测仪（Hobo 型，美国ONSET 公司）自动记录空气温湿度和太阳总辐射等气象要素，采集频率为5 min/次。温室外小气候观测使用当地气象观测站气象要素，采集频率为1 h/次。

1.3 模型检验统计变量

2 模型的建立

2.1 发育进程算法

2.1.1 温差法

温差法（accumulated temperature difference,MTD）分为最大温差法（式（1））和昼夜温差法（式（2）），最大温差由最高气温与最低气温的差值计算出，昼夜温差由昼间均温与夜间均温的差值计算出。由此计算出累积最大温差（accumulated temperature difference between maximum and minimum,MMM）和累积昼夜温差（accumulated temperature difference between day and night,MDN），由于逐时气象数据没有昼夜之分，故本研究逐时算法中不考虑MDN。

式中Tmax为最高气温，℃；Tmin为最低气温，℃；Tday为昼间均温，℃；Tnight为夜间均温，℃，下述逐时算法均与逐日算法形式一致，不再赘述。

2.1.2 积温法

积温法（MAT）分为有效积温法和活动积温法（式（3）～（8）），其中有效积温由平均温度与生物学下限温度（表2）的差值计算出，活动积温由平均温度计算出。由此计算出累积有效积温（effective accumulated temperature,ME）和累积活动积温（active accumulated temperature,MA）[9]，本研究设置了3 种作物生长发育速率对温度响应的线性生长假设：

表2 不同园艺作物不同发育阶段三基点温度Table 2 Three base point temperature of different horticultural crops at different development stages ℃

1）假设A：仅考虑生物学下限温度Tb，生长发育速率随温度的增加而线性增长，计算见式（3）～（4）：

2）假设B：考虑生物学下限温度Tb和 生物学上限温度Tm，发育速率随温度升高而线性增长，当温度高于Tm时，发育停止，计算见式（5）～（6）：

3）假设C：考虑生物学下限温度Tb和生物学上限温度Tm，发育速率随温度升高而线性增长，当温度高于Tm时，发育随温度升高而保持恒定不变，计算见式（7）～（8）：

式中Tave为平均温度，℃，本研究采用气温和叶温2 种温度形式；Tb为生物学下限温度，℃；Tm为生物学上限温度，℃；MEA为基于假设A 的累积有效积温；MAA为基于假设A 的累积活动积温；MEB为基于假设B 的累积有效积温；MAB为基于假设B 的累积活动积温；MEC为基于假设C 的累积有效积温；MAC为基于假设C 的累积活动积温。通过整理文献资料，得各园艺作物不同发育阶段的三基点温度（表2）。

2.1.3 生理发育时间法

根据温度效应因子（temperature effect factor,FTE）的计算方式将生理发育时间（physiological developmental time,MPDT）分为线性（式（9））和正弦（式（10））2 种类型温度响应模式，再结合2 种温度形式，可将MPDT划分为4 种计算类型（式（11））。为了满足各园艺作物生长发育过程中对光长的需求，本研究试验过程中进行相应补光/遮光田间管理，即仅考虑FTE对园艺作物发育期和采收期的影响，以期探讨温度响应模式和温度形式对园艺作物发育期和采收期的影响。

式中FTE有4 种计算方式，其中FTE_1表示线性温度响应模式和气温形式，FTE_2表示正弦温度响应模式和气温形式，FTE_3表示线性温度响应模式和叶温形式，FTE_4表示正弦温度响应模式和叶温形式。

2.1.4 辐热积法

相对辐热积（relative product of thermal effectiveness and photosynthetically active radiation,RTEP）是将温度效应因子与相应光合有效辐射的乘积，与MPDT算法一致，可将RTEP划分为4 种计算类型（式（12）），然后累加得到累积辐热积（accumulated product of thermal effectiveness and photosynthetically active radiation,MTEP）（式（13））。

式中R为光合有效辐射，MJ/（m2·d）；Q为太阳总辐射，MJ/（m2·d）；K为光合有效辐射占比，一般取0.47[13]。

2.2 算法集成

集成法（methods integration,MI）是将上述4 类建模方法的模拟结果进行4 种方法的集成，计算见式（14）～（17）。

式中MI_ave为平均值法集成法模拟值，d；MI_max为最值均值法集成法模拟值，d；MI_med为中值法集成法模拟值，d；MI_step为逐步回归法集成法，d；N为建模方法数量，其中逐时尺度有21 种方法，逐日尺度有22 种方法（表3）；SN为第N种建模方法模拟值，d；SN_max为N种建模方法模拟值最大值，d；SN_min为N种建模方法模拟值最小值，d；SN_med为N种建模方法模拟值中值，d；aN为第N种建模方法逐步回归系数，c为逐步回归常数，均为模型待定参数。

表3 4 类建模方法算法介绍Table 3 Introduction to algorithms for four types of modeling methods

3 模型检验与结果分析

3.1 模型参数的确定

3.1.1 发育期及采收期模拟模型参数的确定

根据上述4 类发育进程算法（式（1）～（12）），按最小二乘法原理求解出发育期模型初始参数，并通过试错法调试参数，最终得到发育期模拟模型参数。由于黄瓜、番茄和茶叶是无限生长型作物，在采收初期-拉秧期阶段内会经历多次采收过程[1,23-24]，其季节性生产能力和产量与采收期的关系巨大，根据每次采收阶段所需积累的发育进程与产量间的定量关系，按最小二乘法原理求解出不同作物采收期模型参数，其中式（18）适用于逐时尺度下JY 品种采收期模拟，式（19）适用于其他情况下采收期模拟。

式中Y为作物产量，kg/hm2；Y0为作物初次采收产量，kg/hm2；DN为第N种建模方法发育期进程累积量；D0为第N种建模方法初次采收发育期进程累积量；a和b为模型待定参数。

3.1.2 发育期及采收期模拟模型集成参数的确定

由于部分试验没有获取辐射数据，故本研究逐步回归参数分为2 种情况，即有辐射数据（包含辐热积算法）和没有辐射数据（不包含辐热积算法），最终按最小二乘法原理求解出模型逐步回归系数aN和逐步回归常数c（表4）。

表4 发育期和采收期模拟模型逐步回归集成参数Table 4 Stepwise regression integrated parameters of simulation models for development and harvest periods

3.2 模型的检验

3.2.1 发育期模拟模型的检验

利用表1 中21 组验证数据集对发育期模型进行验证，结果表明，发育期模拟模型整体的观测值（45.52±34.02 d）与模拟值（45.39±34.44 d）的RMSE 为6.29 d，NRMSE 为13.82%，R2为0.97，D值为0.99，表明模型具有较高的模拟精度。由表5 可知，不同时间尺度发育期模拟模型的RMSE 为4.85～7.42 d，NRMSE 为10.65%～16.31%，且不同时间尺度发育期模型模拟精度大小依次是逐时尺度、逐日尺度。不同作物发育期模拟模型的RMSE 为0.50～7.83 d，NRMSE 为4.33%～20.24%，且不同园艺作物发育期模型模拟精度大小依次是郁金香、番茄、香芹、黄瓜、芹菜、菠菜。不同发育阶段发育期模拟模型的RMSE 为3.59～9.79 d，NRMSE 为9.99%～18.93%，且不同发育阶段发育期模型模拟精度大小依次是4th、5th、3rd、2nd、1st发育阶段。不同模拟方法发育期模拟模型的RMSE 为0.08～24.37 d，NRMSE 为0.18%～54.81%，其中，不同建模方法发育期模型模拟精度大小依次是集成法MI（均值RMSE 为1.46 d，均值NRMSE 为3.28%，模拟效果优）、积温法MAT（均值RMSE 为4.63 d，均值NRMSE 为10.40%，模拟效果良好）、生理发育时间法MPDT（均值RMSE 为5.51 d，均值NRMSE 为12.40%，模拟效果良好）、辐热积法MTEP（均值RMSE 为8.68 d，均值NRMSE 为15.96%，模拟效果良好）、温差法MTD（均值RMSE 为15.46 d，均值NRMSE 为34.76%，模拟效果较差）；不同温度响应模式发育期模型模拟精度大小依次是正弦式FTE_2（均值RMSE 为6.83 d，均值NRMSE 为13.67%，模拟效果良好）、线性式FTE_1（均值RMSE 为7.36 d，均值NRMSE 为14.68%，模拟效果良好）；不同线性温度响应模式发育期模型模拟精度大小依次是假设B（均值RMSE 为4.25 d，均值NRMSE 为9.57%，模拟效果优）、假设C（均值RMSE 为4.74 d，均值NRMSE 为10.65%，模拟效果良好）、假设A（均值RMSE 为5.18 d，均值NRMSE 为11.64%，模拟效果良好）；不同温度形式发育期模型模拟精度大小依次是叶温形式（均值RMSE 为5.25 d，均值NRMSE 为11.14%，模拟效果良好）、气温形式（均值RMSE 为5.98 d，均值NRMSE 为12.69%，模拟效果良好）。综上所述，最优园艺作物发育期模拟路径为先选择逐时尺度、考虑生物学下限和上限温度的正弦温度响应模式、叶温温度形式构建发育期模型，再选择中值集成法优化模型。

表5 基于不同发育进程算法的发育期模拟模型验证统计量Table 5 Verification statistics of development period simulation model based on different developmental process algorithms

3.2.2 采收期模拟模型的检验

利用表1 中12 组验证数据集对采收期模拟模型进行验证，结果表明，采收期模拟模型整体的观测值（114.28±37.41 d）与模拟值（114.43±34.10 d）的RMSE 为15.96 d，NRMSE 为13.97%，R2为0.82，D值为0.95，表明模型具有较高的模拟精度。由表6 可知，不同尺度采收期模拟模型的RMSE 为14.79～17.01 d，NRMSE 为12.94%～14.89%，且不同时间尺度采收期模型模拟精度大小依次是逐时尺度、逐日尺度，与发育期模拟模型的结果一致。不同作物采收期模拟模型的RMSE 为8.33～17.08 d，NRMSE 为5.44%～15.74%，且不同园艺作物采收期模型模拟精度大小依次是番茄、茶、黄瓜。不同模拟方法采收期模拟模型的RMSE 为11.92～20.09 d，NRMSE 为10.45%～17.60%，其中，不同建模方法采收期模型模拟精度大小依次是MI（均值RMSE 为14.10 d，均值NRMSE为12.35%，模拟效果良好）、MTEP（均值RMSE 为15.54 d，均值NRMSE 为13.51%，模拟效果良好）、MAT（均值RMSE 为16.17 d，均值NRMSE 为14.17%，模拟效果良好）、MTD（均 值RMSE 为16.38 d，均值NRMSE 为14.35%，模拟效果良好）、MPDT（均值RMSE 为17.02 d，均值NRMSE 为14.91%，模拟效果良好）；不同温度响应模式采收期模型模拟精度大小依次是FTE_2（均值RMSE 为16.02 d，均值NRMSE 为13.99%，模拟效果良好）、FTE_1（均值RMSE 为16.54 d，均值NRMSE 为14.43%，模拟效果良好）；不同线性温度响应模式发育期模型模拟精度大小依次是假设C（均值RMSE 为14.89 d，均值NRMSE 为13.05%，模拟效果良好）、假设A（均值RMSE 为16.61 d，均值NRMSE 为14.55%，模拟效果良好）、假设B（均值RMSE 为17.39 d，均值NRMSE 为15.23%，模拟效果良好）；不同温度形式发育期模型模拟精度大小依次是叶温形式（均值RMSE为15.22 d，均值NRMSE 为13.32%，模拟效果良好）、气温形式（均值RMSE 为17.21 d，均值NRMSE 为15.05%，模拟效果良好）。综上所述，最优园艺作物采收期模拟路径为先选择逐时尺度、考虑生物学下限和上限温度的正弦温度响应模式和叶温温度形式构建采收期模型，再选择逐步回归集成法优化模型。

表6 基于不同发育进程算法的采收期模拟模型验证统计量Table 6 Verification statistics of harvest period simulation model based on different developmental process algorithms

综上可知，园艺作物发育期和采收期模型逐时优于逐日尺度（RMSE 相差2.40 d，NRMSE 相差3.81%）；集成法模拟优于单一模拟（RMSE 相差3.59 d，NRMSE相差6.28%）；FTE_2优于FTE_1（RMSE相差0.52 d，NRMSE 相差0.73%）；叶温优于气温温度形式（RMSE相差1.64 d，NRMSE 相差1.58%）；考虑生物学下限和上限温度的温度响应模式优于仅考虑生物学下限温度的温度响应模式（RMSE 分别相差0.10 d 和0.39 d，NRMSE 分别相差0.13%和0.82%）。

4 讨论

本研究表明时间尺度对发育期和采收期模拟模型的模拟精度有影响，小时步长模拟精度要高于日步长，这与PURCELL等[18,30]研究结果一致。作物类型对模型模拟精度也有影响，发育期模拟模型的模拟精度依次为花卉类、瓜果类、叶菜类，而灌木类和瓜果类采收期模拟模型的模拟精度相近，未来研究还需要增加试验样本以明确该结论的普适性。本研究中采用的非线性模型是正弦式温度响应模式，与线性模型相比，非线性模型更能较好地提高模型模拟精度，这与WANG等[31]研究结果一致。彭世彰等[32]探究了叶片气孔温度差对气孔导度、蒸腾速率和光合作用的影响，故本研究分别采用了叶温和气温的温度形式，结果表明叶温可以提高模型模拟精度。CHEN等[33-34]表示由于气候变暖，作物发育期将提前且生长期将缩短，说明极端温度对作物发育期和采收期有影响，本研究结果表明考虑下限和上限温度的温度响应模式可以提高模型模拟精度。本研究中集成法模型要显著优于单一模型的精度，这与SYED等[3-5]研究结果一致，但是发育期模块集成法模拟与单一法模拟精度差值较大，而采收期模块集成法模拟与单一法模拟精度差值较小。

通常用试错法和步长加速法等调试至较适宜但不是最适宜的作物模型参数，即模型可能存在多参同效的情况[35]，还需进一步给合试验观测资料进行论证，提高模型参数的普适性和合理性。本研究仅考虑4 种集成方法，未来可探究多种集成方法优化模型的模拟效果[36]，提高模型的精确度。当温度处于极端温度时，模型的参数不确定性增加，导致模型模拟结果出现较大偏差，本研究需要通过基于环境控制试验数据开发温度响应函数来降低模型的参数不确定性，对模型进行校正和检验，提高模型的普适性和精确度。

5 结论

本研究建立了基于4 类建模方法（温差法、积温法、生理发育时间法和辐热积法）的发育期和采收期模拟模型，并以4 种方式（平均值、最值均值、中值和逐步回归）集成模拟结果，为园艺作物智慧生产管理和高效利用农业资源方面提供理论基础和技术支撑。

1）不同时间尺度模拟模型的均方根误差（root mean square error，RMSE）为4.85～17.01 d，归一化均方根误差（normalized root mean square error，NRMSE）为10.65%～16.31%；不同作物模拟模型的RMSE 为0.50～17.08 d，NRMSE 为4.33%～20.24%，最优发育期模拟模型的园艺作物是郁金香，最优采收期模拟模型的园艺作物是黄瓜；不同模拟方法模拟模型的RMSE 为0.08～24.37 d，NRMSE 为0.18%～54.81%。

2）逐时优于逐日尺 度（RMSE相差2.40 d，NRMSE 相差3.81%）；集成法模拟优于单一模拟（RMSE 相差3.59 d，NRMSE 相差6.28%）；FTE_2优于FTE_1（RMSE 相差0.52 d，NRMSE 相差0.73%）；叶温优于气温温度形式（RMSE 相差1.64 d，NRMSE 相差1.58%）；考虑生物学下限和上限温度的温度响应模式优于仅考虑生物学下限温度的温度响应模式（RMSE分别相差0.10 d 和0.39 d，NRMSE 分别相差0.13%和0.82 %）。

3）最优园艺作物发育期模拟路径为先选择逐时尺度、考虑生物学下限和上限温度的正弦温度响应模式和叶温温度形式构建模型，再选择中值集成法优化模型。最优园艺作物采收期模拟路径为先选择逐时尺度、考虑生物学下限和上限温度的正弦温度响应模式和叶温温度形式构建模型，再选择逐步回归集成法优化模型。