基于增強回歸樹的海河平原小麥赤霉病預測模型構建與驗證

陶晡，齊永志，屈赟，曹志艷，趙緒生，甄文超

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%，病穗率預測值與實際觀測值的波動趨勢基本一致。研究結果不僅為海河平原小麥赤霉病預測預報提供技術支撐，也為小麥赤霉病預測模型優化和改進提供了參考。