基于增强回归树的海河平原小麦赤霉病预测模型构建与验证

陶晡，齐永志，屈赟，曹志艳，赵绪生，甄文超

1河北农业大学植物保护学院，河北保定 071001；2河北农业大学现代教育技术中心，河北保定 071001；3河北农业大学农学院/华北作物改良与调控国家重点实验室/河北省作物生长调控重点实验室，河北保定 071001

0 引言

【研究意义】小麦赤霉病（Fusarium head blight，FHB）是小麦生产上发生面积最广、危害程度最大的麦类病害之一[1]，该病是以禾谷镰孢（Fusarium graminearum）为主要致病菌的真菌性病害[2]。据报道，自 1990年以来美国小麦种植面积因赤霉病流行不断压缩，2018年小麦种植面积减少1 200万公顷[3]，2016年小麦赤霉病造成加拿大萨斯克彻温省经济损失约为10亿美元[4]。同时病菌产生的脱氧雪腐镰刀菌烯醇（DON毒素）和玉米赤霉烯酮（ZEN毒素）等危害人畜健康，对小麦品质和产量造成严重影响[5-6]。近年来，因气候条件变化、耕作制度改变，我国小麦赤霉病发生呈日趋严重的趋势，由长江中下游麦区逐渐向北扩展，淮河流域地区成为重发区，在黄淮北片麦区和北部冬麦区也成为常发病害[7]，2010年以来重发频率在50%以上，2015、2016、2018年发生面积均超过550万公顷[8]。自1995年以来，小麦赤霉病逐渐在海河平原（也称河北平原）蔓延，已由零星出现逐渐演变成连片发生，并由次要病害上升为主要病害之一，年均发生面积达26.7万公顷以上[7]。小麦赤霉病在流行年份具有短期内暴发快、面积大、损失重的特性，因此，明确海河平原影响赤霉病发生的关键气象因子，建立适宜该区域的病害预测模型，提供准确的预测预报信息，对有效防控病害蔓延具有重要意义。【前人研究进展】DE WOLF等[9-10]以小麦开花前7 d和开花后10 d的气象因子作为预测变量，用逻辑回归建立了小麦赤霉病测报模型；在此基础上，SHAH等以品种抗性、玉米残茬以及前期研究获得的4个气象因子作为变量，通过R语言，建立了基于Leaps and Bounds算法的Logistic回归模型[11]和增强回归树（boosted regression tree，BRT）模型[12]，结果表明，增强回归树模型误判率低于Logistic回归模型；HOOKER等[13]以抽穗前4—7 d降雨天数和温度为预测变量，建立了含有指数项的模型，预测小麦 DON毒素含量；DEL PONTE等[14]以空中孢子捕捉量和感病组织为基础建立测报模型；ROSSI等[15]以菌源量、小麦关键生育期为基础，综合考虑日产孢率、孢子分散率、侵染机率和小麦生育期等因素，预测小麦赤霉病发生风险；MUSA等[16]建立了基于web的瑞士小麦赤霉病预警系统FusaProg，预测小麦赤霉病发生、DON毒素含量并指导杀菌剂科学使用。国内专家学者从气候预测、菌量预测、气候菌量相结合预测等方面展开了研究，建立了长期预测、中期预测与短期预测模型，同时，借助神经网络[17-18]、支持向量机[19]、无人机高光谱图像[20]等技术，不断提高了预测预报的准确度。一般情况下，预测模型存在可移植性差、跨地区应用准确度下降等问题。【本研究切入点】增强回归树是以分类回归树（classification and regression tree，CART）算法为基础的一种自学方法，通过自我学习和随机选择生成多重回归树，提高模型稳定性和预测精度。YOU等利用该模型明确了环境变量与品种对牧草病害发生的影响，取得了较好的效果[21]，为评估小麦赤霉病主要影响因子重要性提供了一种新的思路。【拟解决的关键问题】根据影响小麦赤霉病流行的关键生育期，选择温度、湿度、降雨、日照、风速等气象因子为预测变量，筛选出重要预测变量，并分析其对病害发生的影响，以期提升模型预测准确度，为小麦赤霉病发生预测预报提供参考，同时也可为建立该病害综合、高效防控体系提供技术支撑。

1 材料与方法

试验数据采集于2014—2016年完成，数据处理、模型构建及检验于2017—2019年完成。

1.1 病害数据和气象信息采集

收集整理 2001—2013年海河平原小麦主产区安新、望都、定州、新乐、正定、无极、栾城、辛集、平山、行唐、灵寿、阜城、武邑、景县、临西、宁晋、磁县、馆陶、曲周、永年、大名共21个县（市）定点监测小麦赤霉病发病情况基本数据，来源于河北省植保植检总站。

2014—2016年在上述21县（市）田间调查小麦赤霉病发生情况，每县选择10个调查点，每个调查点随机取样500穗，并计算病穗率。依据国家《小麦赤霉病测报技术规范》GB/T15796—2011将小麦赤霉病划分为不发生、轻度流行、中度流行、重度流行4个等级：不发生（0级，病穗率＜0.1%，对小麦生产未造成减产）、轻度流行（1级，0.1%≤病穗率＜5%，对小麦生产造成局部减产）、中度流行（2级，5%≤病穗率＜10%，对小麦生产造成部分减产）、重度流行（3级，病穗率≥10%，对小麦生产造成明显减产）。

小麦生育期内气象资料：来源于河北省气象局2001—2016年21个县（市）逐日最高温度、最低温度、平均温度、日照数据、平均风速、平均相对湿度、总降雨量等。

1.2 数据处理

河北省南部和北部冬小麦生育进程存在一定的时间差[22]，南部麦区比中部麦区的播种期、越冬期晚 3—5 d，其他生育期早3—7 d；北部麦区比中部麦区的播种期、越冬期早5—7 d，其他生育期晚5—7 d（表1）。按照河北中部麦区常年小麦生育进程分别计算21个县（市）小麦抽穗期初始日期，每5 d编为一组，分别以当地的小麦抽穗期初始日为起点，小麦抽穗期初始日向前选择6组，即抽穗期初始日向前26—30、21—25、16—20、11—15、6—10 和 1—5 d，抽穗期初始日向后选择2组，即抽穗期初始日向后1—5、6—10 d。以中部麦区的正定、栾城等地最高温度为例，选择小麦抽穗期初始日（5月1日）为起点，HT-65、HT-55、HT-45、HT-35、HT-25、HT-15分别代表小麦抽穗期初始日之前 26—30、21—25、16—20、11—15、6—10和1—5 d的最高温度平均值，HT15、HT25分别代表小麦抽穗期初始日之后1—5、6—10 d的最高温度平均值。主要包括9个气象因子：最高温度（HT）、最低温度（LT）、平均温度（MT）、平均风速（MWS）、平均相对湿度（MRH）、相对湿度≥85%天数（RH85）、降雨天数（DRain）、总降雨量（Rain）、总日照时数（SD）。

表1 河北省中部麦区常年小麦生育进程Table 1 Growth process of wheat in middle wheat region in Hebei Province[22]

1.3 预测模型构建及参数选择

1.3.1 预测模型构建 研究以不同气象因子为预测变量、小麦赤霉病病穗率为响应变量，采用增强回归树模型建模，模型拟合使用R语言（3.6.1版本，R核心开发组，2019）gbm包和 Elith的函数包[23]。增强回归树结合了提升（boosting）和分类回归树（CART）两种技术，通过组合大量相对简单的决策树的方式以优化模型的预测性能，模型可写成M棵分类回归树相加的形式。

式中，Tm(X, γm)为第m棵分类回归树，X为气象因子预测变量，γm为其参数，是该决策树分裂点和叶子结点的赋值，求解 γm的过程即为单棵决策树的学习过程。

1.3.2 模型参数选择 在模型运行过程中，需要优化迭代次数（the number of trees，nt）、树的复杂度（tree complexity，tc）、学习效率（learning rate，lr）、抽样比率（bag fraction，bf）、函数损失形式（distribution）、交叉验证折数（cv.folds）等参数[12]。树的复杂度即为单棵决策树的叶节点数量，它是模型拟合环境因子间交互作用的阶数。增强回归树中所有决策树的叶节点数量相同，训练过程中叶节点达到一定数量时则停止生长，不需要剪枝[23-24]。学习效率决定了模型达到最优所需训练的时间，lr值过小，则收敛速度慢、训练时间越长；lr值过大，抽样时容易产生噪音，导致函数平滑性降低、稳定性差[25]。通常情况下，迭代次数（nt）要达到1 000以上模型才趋于稳定，树的复杂度（tc）1—16，学习效率（lr）0.001—0.1，抽样比率（bf）为 0.75，函数损失形式为“gaussian”。由于 tc和 lr的取值影响模型的预测准确度，随机选择70%训练集数据用于构建模型，剩余30%的数据用于计算模型的预测偏差，根据模型预测偏差大小选择最优的tc值和lr值。

随着决策树数量的增加，增强回归树模型的拟合效果会越来越好，但决策树数量过大会出现过拟合，导致预测精度降低。本研究以 10倍交叉验证法（10-fold cross-validation）确定最优决策树的数量。

1.3.3 预测因子相对重要性计算 在分类回归树模型中，FRIEDMAN[26]提出了用I2(T)j 作为第j个预测因子Xj的相关性的度量，该度量基于选择Xj变量进行决策树的节点分裂时平方误差加权改进，该度量比其对应的单个分类树更加可靠。

式中，预测因子Xj的相对重要性的平方即为平方误差加权改进在模型中M棵分类回归树上的平均[27]。通常情况下预测因子的相对重要性以百分数形式表示，所有预测因子的相对重要性之和为100。

1.4 模型预测准确度检验

采用最大误差参照法计算预测准确度[28]：

其中，R为模型预测准确度，Fi为模型预测病害流行等级，Ai为实际病害流行等级，Mi为第i次预测的最大参照误差。

2 结果

2.1 模型预测变量确定

通过对预测集正态性判断、误差独立性判断、线性判断、同方差性判断以及多元回归模型综合判断，剔除离群值和强影响点，筛选出影响海河平原小麦赤霉病发生的 8个关键气象因子，即 LT-65、MWS-55、MRH-55、Rain-35、MT-25、SD15、MRH15、DRain15。同时，构建了小麦赤霉病多元线性回归预测模型（multiple linear regression model，MLRⅠ）：y=-13.2427+0.3145LT-65-0.9824MWS-55+0.1209MRH-55+0.1377Rain-35-0.4184MT-25+0.0814SD15+0.28024MRH15-0.8832DRain15，该模型R2=0.8158，矫正R2=0.8018，P＜2.2×10-16。

其中，y为小麦赤霉病病穗率，LT-65为抽穗期初始日之前26—30 d最低温度，MWS-55为抽穗期初始日之前21—25 d平均风速，MRH-55为抽穗期初始日之前21—25 d相对平均湿度，Rain-35为抽穗期初始日之前11—15 d总降雨量，MT-25为抽穗期初始日之前6—10 d平均温度，SD15为抽穗期初始日之后1—5 d总日照时数，MRH15为抽穗期初始日之后1—5 d相对平均湿度，DRain15为抽穗期初始日之后 1—5 d降雨天数。

2.2 预测模型拟合的关键影响因素

根据增强回归树（BRT）模型拟合曲线（图 1）可知，在不同学习效率（lr）和树的复杂度（tc）下，当lr为0.1和0.05时，模型的最小预测偏差与其他学习效率相比偏差相对较大，在不同树的复杂度下，模型会较早的发生过度拟合。当lr为0.001时，模型迭代次数一般在2 000左右达到最小预测偏差，当lr为0.01时，模型迭代次数一般在500—800范围内达到最小预测偏差。当lr为0.005时，模型迭代次数在900—1 800范围内达到最小预测偏差。

设置lr为0.01、0.005，由不同tc的残差标准误（residual standard error，图2）可知，在lr为0.01和0.005的学习效率条件下，当tc为6时增强回归树模型的残差标准误分别为0.01004和0.006311，随着tc值的增加，增强回归树模型的预测偏差相对变化不大。综合考虑不同lr和tc下模型预测偏差，选择模型的lr为0.005，tc为6。

2.3 预测变量重要性及对病害发生预测的影响

设置lr为0.005，在不同tc下预测了变量的重要性（图3），随着tc增大，各预测变量重要性排名未发生太大变化，MRH15、Rain-35是相对重要的两个预测变量；其次是 MRH-55、SD15、LT-65、MT-25、MWS-55、DRain15。

当tc=6时，预测变量MRH15、Rain-35、MRH-55、SD15、LT-65、MWS-55、MT-25、DRain15的重要性由高到低依次为69.62%、14.08%、4.89%、4.34%、3.35%、2.02%、1.20%、0.50%。

设置lr为0.005、tc为6、抽样比率（bf）为0.75、函数损失形式为“gaussian”，交叉验证折数为10次、n.trees为5 000，确定拟合最终的拟合模型，模型各预测变量的反应曲线见图4。

（1）平均相对湿度对小麦赤霉病发生风险的影响预测变量 MRH15对赤霉病发生风险的重要性最高，为69.62%。当其小于46%时，其变化对赤霉病发生的影响较小；当其在46%—67%时，随其增加，赤霉病发生风险迅速上升；当其高于67%时，其对赤霉病发生的促进作用趋于平稳。预测变量MRH-55对赤霉病发生风险的重要性居第3位，为4.89%。其对赤霉病发生风险是非线性关系，当其在35%—48%时，其对赤霉病发生风险的作用效果为先升后降；当其在48%—58%时，随其增加，赤霉病发生的风险迅速上升；当其高于58%时，其对赤霉病发生的促进作用趋于平稳。

（2）总降雨量对小麦赤霉病发生风险的影响 预测变量Rain-35对赤霉病发生风险的重要性居第2位，为14.08%。赤霉病的发生风险随其增加而增加，当其高于16 mm时，其对赤霉病发生风险的促进作用趋于平稳。

（3）日照时数对小麦赤霉病发生风险的影响 预测变量SD15对赤霉病发生风险的重要性居第4位，为4.34%。模型显示当该时期日照时数在25—49 h时，赤霉病发生风险迅速降低，当其大于50 h时，日照时数对赤霉病发生风险的抑制作用趋于平稳。

2.4 模型的优化

根据增强回归树模型筛选的 4个重要的预测变量，即 MRH15、Rain-35、MRH-55、SD15，简化多元线性回归模型（MLRⅡ）：y=-19.45376+0.11689MRH15+0.17346Rain-35+0.04185SD15+0.26592MRH-55，该模型R2=0.7575，矫正 R2=0.7468，P＜2.2×10-16。

2.5 模型的验证

以2008、2010、2012年“同年份多点”的部分地区历史数据为测试集，对多元线性回归模型、增强回归树模型预测结果进行验证，观测值与预测值对比发现，预测结果与实际观测值基本相符，但个别地区预测结果略有出入（图5）。预测变量由8个简化为4个时，多元线性回归模型预测准确度由 88.43%降至85.90%，增强回归树模型预测准确度由 87.72%升至91.23%。

以2001—2016年“多年份定点”的正定、栾城两地小麦赤霉病发生的历史数据为测试集，对多元线性回归模型、增强回归树模型的预测结果进行验证，预测结果与实际观测值曲线基本一致（图6）。预测变量由8个简化为4个时，两个模型预测准确度无显著变化，多元线性回归模型预测准确度由87.53%变为87.42%，增强回归树模型预测准确度由89.20%变为89.21%。

3 讨论

小麦赤霉病是典型的气候性病害[1]，其预测模型受品种抗性、气象因素、田间菌源量等多因素影响，具有典型的地域特异性。GIROUX等评价了 9种不同模型在加拿大魁北克地区预测小麦赤霉病发生或DON毒素含量的效果，美国的两个模型（DE WOLF等开发[9-10]）和阿根廷的模型（MOSCHINI等开发[29-30]）预测效果优于其他模型[31]，表明预测模型

应用存在一定的地域特异性。SHAH等[12]利用增强回归树模型预测了严重度大于10%的小麦赤霉病发生概率，测试数据误识率与 logistic回归模型相比下降31%，但该模型未能准确反映出病穗率和病害发生等级；LANDSCHOOT等[32]基于比例优势模型，明确了比利时小麦赤霉病病情指数和 DON毒素含量的关键参数及参数值，并开发了web系统，预测效果良好，但该系统无法在其他地区应用；XU等[33]开发了基于logistic回归的欧洲小麦赤霉病DON毒素含量预测模型，利用重采样和全子集回归分析表明，通过气象因子预测DON毒素含量的“最佳”模型存在不唯一性；基于气象因子的 DON毒素含量预测模型具有可移植性或可替代性，其通用性较差，预测结果易出现假阳性。本实验室前期研究[34]对安徽桐城小麦赤霉病预测模型中的气象因子进行了物候期数据的本地化修正，构建了河北南部麦区小麦赤霉病预测模型，历史数据

验证模型准确度为70.00%，与安徽桐城市小麦赤霉病病穗率的预测结果相比，其准确度相对较低。

本研究综合借鉴本领域专家的研究经验[9,12,35-40]，以小麦物候期的抽穗期初始日为起点，每5 d编为一组，抽穗期初始日向前选择编6组，向后选择编2组，以各组内最高温度、最低温度等9个气象因子为自变量，赤霉病病穗率为因变量，通过逐步回归分析，筛选出显著的变量，得到最优子集[41]。观测样本中异常值对预测结果产生影响，利用R语言car包中的outlierTest( )函数查找离群点、Cooke距离判断强影响点，通过消除预测数据集中的离群点和强影响点，提升模型预测准确度。同时，通过对预测集正态性判断、误差独立性判断、线性判断、同方差性判断以及多元回归模型综合判断，证实了气象因子与病穗率的多元线性回归假设。SHAH等[12]研究表明，增强回归树模型预测效果优于logistic回归模型；本研究以“同年份多点”和“多年份定点”历史数据验证了多元线性回归模型和增强回归树模型预测结果，当预测变量由8个调减至4个时，多元线性回归模型预测准确度呈下降趋势，而增强回归树模型预测准确度呈上升趋势。利用历史数据验证模型预测效果时，2008年阜城地区和 2010年曲周地区病穗率的预测值与实际观测值偏差较大，其原因可能与本研究仅选择气象因子有关，赤霉病的发生与发展还与小麦品种抗性、田间管理措施等因素有关。在今后的研究中，还需进一步考虑小麦品种抗性、田间菌源量等因素，以期进一步提高预测准确度。

4 结论

研究构建了基于增强回归树的海河平原小麦赤霉病病穗率预测模型，该模型含有4个预测变量。经两地16年历史数据验证，模型预测准确度为89.21%，病穗率预测值与实际观测值的波动趋势基本一致。研究结果不仅为海河平原小麦赤霉病预测预报提供技术支撑，也为小麦赤霉病预测模型优化和改进提供了参考。