群體小麥條銹病發病動態無人機遙感監測方法

蘇寶峰，劉昱麟，黃彥川，蔚 睿，曹曉峰，韓德俊

（1. 西北農林科技大學機械與電子工程學院，楊凌 712100； 2. 農業農村部農業物聯網重點實驗室，楊凌 712100；3. 陜西省農業信息感知與智能服務重點實驗室；4. 西北農林科技大學農學院，楊凌 712100；5. 西北農林科技大學旱區作物逆境生物學國家重點實驗室，楊凌 712100）

0 引 言

無人機成像多光譜遙感在作物病害檢測方面已經取得了一定的成果。Liu 等[13]通過無人機搭載傳感器采集感染赤霉病的小麥灌漿期高光譜圖像，利用誤差反向傳播（Back Propagation，BP）算法建立了赤霉病監測模型，總體準確率達98%。Rodríguez 等[14]對比了5 種機器學習算法在馬鈴薯晚疫病檢測和評估中的性能表現，證明了可以通過無人機多光譜技術采集馬鈴薯冠層反射率，并建立線性支持向量機（Liner Support Vector Machine，LSVM）和隨機森林分類模型來檢測和評估馬鈴薯晚疫病發病狀況。Guo 等[15]通過無人機搭載光譜成像傳感器采集小麥條銹病冠層反射率，提取其植被指數和紋理特征，建立了不同感染期的偏最小二乘回歸的小麥條銹病監測模型，病害指數與所提取特征的決定系數可達0.82。Su等[16]利用無人機采集了接種條銹病的小麥生長季的冠層多光譜數據，利用不同的冠層植被指數在不同發病期識別條銹病發病區域。Bhandari 等[17]通過無人機獲取小麥冠層RGB 圖像計算植被指數，證明了其與條銹病侵染系數之間具有較好的相關性。然而上述針對作物病害識別監測的研究中，只對少數品種進行分析，所得出的結論在多基因型育種群體中的應用還有待進一步研究。Chivasa 等[18]基于無人機多光譜技術檢測玉米條紋病，建立不同品種抗病等級的隨機森林分類模型，表明無人機成像光譜能夠提高作物表型鑒定效率。上述研究表明無人機低空遙感技術不但能夠提高對作物病害監測效率，還能為育種群體的表型分析提供支持[19-21]，然而僅通過單一時間點的冠層數據無法實現對作物發病情況的動態分析[22]。

因此，本文以田間感染條銹病的育種小麥自然群體為研究對象，獲取其無人機成像多光譜數據與條銹病嚴重度人工調查數據，提取冠層光譜植被指數，并利用特征選擇和分類算法建立群體小麥田間條銹病發病階段和嚴重度的分類模型，并篩選出對發病階段和病害嚴重度分類敏感的特征，以探究通過多時相的光譜植被指數進行群體小麥田間條銹病動態分析的可行性；并通過指數響應時間序列，量化并分析田間條銹病發病狀況，以期為群體小麥條銹病的發病動態分析提供一種客觀、高效的高通量方法。

1 材料和方法

1.1 試驗區概況及試驗材料

試驗區域位于陜西省咸陽市楊陵區曹新莊試驗農場（108°5′34″E，34°8′17″N），海拔約480 m，屬于溫帶季風性氣候，年均氣溫約13 ℃，年均降水量635 mm。試驗期間的日平均空氣溫濕度情況如圖1 所示，2021 年4月19—26 日，由于降水，日最大空氣濕度保持在較高的水平，且空氣溫度較為適宜，利于條銹病夏孢子繁殖與傳播，使試驗區域內大部分小麥品種受到感染，產生條銹病夏孢子堆，出現春季流行癥狀。

試驗材料為中國小麥聯合攻關和小麥產業技術體系各單位選送參加區試的510 份小麥新品種（系）與重復種植15 次的6 個參考品種（濟麥22、百農207、周麥18、西農511、偃展4110、周麥36）構成的自然群體，按育種試驗增廣設計實行小區種植，共600 個小區，小區大小為1 m2。試驗區域的正射影像如圖2 所示。

1.2 數據采集

1.2.1 無人機成像多光譜數據獲取

本研究中無人機成像多光譜數據采集系統采用深圳大疆創新科技有限公司的經緯Matrice 100 四旋翼無人機作為平臺，搭載美國Micasense公司的RedEdge成像系統，如圖3 所示。

設B層中與Aj相關的因素成對比較判斷矩陣經過了一致性檢驗，求得單排序一致性指標為CI(j)，(j=1，…，m)，相應的，平均隨機一致性指標為RI(j)，CI(j)、RI(j)已經在層次單排序時求得，那么B層總排序隨機一致性比例為：

RedEdge 傳感器具有5 個通道，其名稱、中心波長/半峰全寬分別為：藍（Blue，B）475 nm /20 nm、綠（Green，G）560 nm/20 nm、紅（Red，R）668 nm/10nm、紅邊（Red Edge，RE）717 nm/10nm、近紅外（Near Infrared，NIR）840 nm/40 nm，每個通道空間分辨率為1 080×756 像素；成像系統下行光照傳感器（Downwelling Light Sensor，DLS）模塊用于數據采集時環境光線與太陽角度的校正；平臺集成的全球定位系統（Global Positioning system，GPS）用于記錄圖像的位置信息；反射率校正板用于校正各個通道的反射率。使用該平臺在2020 年3 月31 日至5月14 日采集15 次小麥冠層成像多光譜數據，數據采集時間為能見度良好的10:00－15:00，無人機飛行高度為20～21 m，前向與側向重疊率為80%～90%，飛行速度保持在1.4～1.6 m/s。

1.2.2 小麥條銹病人工鑒定數據采集

按照小麥葉片上條銹病夏孢子所占整個葉片面積的百分比，將嚴重度劃分為13 個等級（0、1%、5%、10%、20%、30%、40%、50%、60%、70%、80%、90%、100%）[23]。在2021 年5 月12 日，按照上述標準對研究區域內600 個小區進行病害嚴重度鑒定數據采集，將其劃分為3 個嚴重度等級，即病害嚴重度在0～30%為高抗反應型（Resistance，R），在40%～60%為中抗反應型（Moderate Resistance，MR），在70%～100%的為感病反應型（Susceptible，S）[24]。

依照發病后時間順序將條銹病發病劃分為5 個階段，對應數據采集日期為：4 月22 日、4 月26 日、5 月1 日、5 月6 日、5 月12 日，分別記錄為S0、S1、S2、S3 和S4，其中S0 為田間第一次出現明顯夏孢子堆的狀態。（Optimized Soil Adjusted Vegetation Index，OSAVI）圖像分割小麥冠層與背景土壤，通過統計比較確定每次多光譜數據的閾值生成二值化掩膜圖像。最后，通過圖像裁剪處理得到試驗群體600 個種植小區的反射率指數圖像和掩膜，并計算和提取各小區光譜植被指數。本文中所用到的光譜植被指數的定義及計算公式如表1 所示（人工選擇的閾值會影響OSAVI 指數的提取結果，故該指數，僅用作閾值分割，不參與后續建模與分析）。

表1 本研究中使用植被指數公式Table 1 Spectral vegetation indices used in this study

1.3 研究方法

首先對無人機多光譜圖像進行預處理。由于光譜植被指數擁有明確的物理意義，能夠探測病害引起的生理生化過程[25]，故構建由光譜植被指數、條銹病發病階段和條銹病嚴重度等級數據組成的數據集；通過特征篩選和機器學習算法建立條銹病發病階段和嚴重度的分類模型，同時篩選出對分類結果敏感的特征；并基于篩選出的光譜植被指數響應的時間序列，分析部分試驗材料在不同階段的田間條銹病發病動態。本研究整體技術路線如圖4 所示。

1.3.1 無人機圖像預處理與植被指數提取

首先，使用Pix4Dmapper（Pix4DInc., Switzerland,http://www.pix4d.com/）軟件對無人機采集的多光譜圖像進行拼接和輻射校正，最后生成5 個通道的單波段反射率指數地圖；為了后續批量處理多次遙感數據，利用QGIS 軟件（https://www.qgis.org/）結合地面控制點對拼接完成的單波段反射地圖進行圖像幾何校正；而后，采用閾值法進行圖像分割，以優化土壤調節植被指數

1.3.2 光譜響應分析

在判斷條銹病嚴重度時，由于群體小麥不同品種（系）間的初始狀態、對病害的敏感程度和染病速率不同，單一時間點的光譜特征無法評估不同品種（系）小麥對條銹病脅迫的響應，為了消除單一時間點的誤差，根據公式（1）計算小麥冠層22 個光譜植被指數的平均變化率作為條銹病脅迫下小麥的光譜響應[37]，并根據響應的時間序列作為構建嚴重度等級分類模型的特征和條銹病的發病動態分析指標。

式中response 為某一植被指數的響應，S0為該植被指數的初始狀態，本文中取4 月22 日的各個光譜植被指數作為初始狀態；Si為S0后第i個時期的光譜植被指數，本文中取i取1、2、3、4，對應1.2.2 中介紹的S1 ～ S4。

1.3.3 特征篩選與分類算法

隨機森林是一種由多決策樹集成的機器學習算法，其通過重采樣方式多次隨機抽取原始訓練集一部分樣本，進行多個決策樹的建模，并通過投票方法輸出最后的分類結果。隨機森林算法在特征篩選上的隨機性使其不容易出現過擬合，抗噪能力較好，且能對特征的重要性進行排序[38]。在對條銹病發病階段的分類模型構建和敏感特征篩選中，利用試驗群體600 個小區5 個發病階段的22 個光譜植被指數（表2 中除OSAVI）構建建模數據集。數據集共3 000 個樣本，每個樣本具有22 個特征，分為S0、S1、S2、S3、S4 共5 個類別，并按照3:1 的比例劃分訓練集與測試集，設置隨機森林算法的決策樹個數為300，建立發病階段的隨機森林分類模型，同時給出22 個特征的重要性排序；在對條銹病病害嚴重度等級的分類模型構建和敏感特征篩選中，以5 月12 日采集的無人機多光譜數據所提取的22 個植被指數的響應作為特征，構建條銹病嚴重度等級隨機森林分類模型（參數設置和數據集劃分與發病階段分類模型相同），并給出22個特征的重要性排序。

為了對比不同算法在篩選條銹病發病階段和病害嚴重度敏感特征的結果，使用隨機蛙跳算法對特征進行篩選。隨機蛙跳是一種能夠利用少量的變量迭代進行建模的高維數據變量選擇方法，它能夠輸出每個變量選擇的可能性，從而進行變量的篩選[39]。支持向量機（Support Vector Machine，SVM）算法在解決小樣本、非線性及高維模式識別中具有一定的優勢[40]，它首先尋找一個最大邊際超平面，并將低維數據通過核函數映射到高維空間中，從而使線性不可分的樣本變為線性可分的樣本，并通過引入模型懲罰因子，提高分類模型的泛化性，達到更好地分類效果。通過隨機蛙跳算法（設置族群中青蛙的數量N為10 000，子族群中青蛙的數量Q設為2）實現對條銹病發病階段和病害嚴重度敏感特征的篩選；為了驗證篩選特征對于分類問題的有效性，將篩選后的被選擇概率排在約前1/3 的（前7 個）植被指數作為SVM的輸入特征，構建條銹病發病階段和病害嚴重度等級分類模型。

2 結果與討論

2.1 群體發育階段的動態分析

根據人工調查群體中600 個小區的生育期結果，該群體4 月16 日—22 日基本完成抽穗，4 月28 日—5 月4日期間群體小麥揚花，揚花授粉后灌漿緩慢啟動，大約兩周內籽粒含水率迅速增加，而干物質的積累很少，群體尚未明顯衰老，期間條銹病是影響小麥長勢的主導因素。這與利用群體冠層NDVI 指數時間序列分析[41]的結果一致：如圖5 所示，NDVI 在拔節期后的飽和效應提前結束。

2.2 條銹病發病階段分類

使用22 個光譜植被指數作為分類特征構建的群體小麥條銹病發病階段分類模型在測試集中的混淆矩陣如圖6a 所示，篩選出的敏感特征如圖7a 所示。通過隨機蛙跳對條銹病發病階段敏感特征篩選的結果如圖7b 所示。將被選擇概率排名約前1/3 的特征作為輸入特征構建的SVM 分類模型在測試集中的混淆矩陣如圖6b 所示。

由圖6 可知，通過隨機蛙跳篩選特征建立的SVM 分類模型和利用22 個植被指數構建隨機森林分類模型對條銹病發病階段的分類都取得了較好的效果，表明可以通過光譜植被指數來描述群體小麥田間條銹病發病階段，模型的評價結果如表2 所示。SVM 算法的F1 分數為0.985，隨機森林算法的F1 分數為0.970，略低于前者，表明隨機蛙跳算法和隨機森林算法都可以用來篩選對于條銹病發病階段更為敏感的指數。S2 和S3 時期的精度和召回率相對較低，從混淆矩陣中也可以看出兩種分類模型在S2 和S3 時期出現相對較為明顯的混淆，可能是由于在兩個階段前后病害表現變化較小，導致光譜植被指數差異較小。

表2 發病階段分類模型評價Table 2 Accuracy evaluation of classification models of disease stages

DVIRE、NDVIrededge、GVI 指數在隨機森林和隨機蛙跳算法中具有相似的排名，對條銹病發病階段較為敏感，表明這些光譜植被指數能夠更好的描述群體小麥田間條銹病發病動態。

2.3 條銹病嚴重度等級分類

使用隨機森林算法建立的條銹病病害嚴重度等級分類模型在測試集中的混淆矩陣如圖8a 所示，篩選出的病害嚴重度敏感特征如圖9a 所示。通過隨機蛙跳對條銹病病害嚴重度敏感特征篩選的結果如圖9b 所示。將被選擇概率排名約前1/3 的特征作為輸入特征構建SVM 分類模型在測試集中的混淆矩陣如圖8b 所示。

如圖8 所示，兩種分類算法分類精度較之對條銹病發病階段的分類精度都有所下降。測試集中錯誤的樣本主要集中在嚴重度相鄰的兩級之間（“R”與“MR”、“MR”與“S”）。由表3 可以看出，通過隨機蛙跳篩選特征建立的SVM 分類模型的F1 分數為0.780，而利用22 個原始特征構建的隨機森林模型的F1 分數為0.741。由于數據集中類別“R”的數量較多，因此在兩種分類模型中類別“R”的精度和召回率最高；兩種分類模型中“MR”的精度和召回率較低，易出現與嚴重度相鄰級別的混淆。由于人工鑒定具有一定的主觀性，給出的分類標簽可能會與實際條銹病嚴重度情況有所差異，造成分類結果并不十分準確，但總體上看，人工鑒定數據與機器學習算法構建的分類模型對條銹病嚴重度等級判斷結果基本一致。由此可知，能夠通過由無人機多光譜數據提取的冠層光譜植被指數響應量化群體中不同小麥品種（系）條銹病病害嚴重程度差異。

表3 病害嚴重度分類模型精度評價Table 3 Accuracy evaluation of classification models of disease severity

如圖9 所示，通過隨機蛙跳算法篩選出較為重要的光譜植被指數與隨機森林算法的分類貢獻性結果相似程度不高，其中WI、NPCI、NDVI 在兩種敏感特征篩選算法中都具有較高的重要性。

圖10 為人工鑒定的條銹病嚴重度等級與分類效果較好的SVM 分類模型結果的可視化。可以更加直觀地看出構建的SVM分類模型基本可以完成對田間小麥條銹病病害嚴重度的分類，個別混淆主要發生在嚴重度等級相鄰的樣本中（“R”與“MR”、“MR”與“S”）。

2.4 條銹病脅迫下光譜響應動態分析

根據兩種分類算法對條銹病發病階段和嚴重度的分類結果，篩選出6 個對條銹病發病階段和病害嚴重度都相對敏感的指數，分別為 DVIRE、GVI、NDVI、NDVIrededge、NPCI 和WI，計算6 個指數相對于4 月22 日的響應，并在群體重復種植的15 組參考品種中隨機選擇一組，分析其冠層植被指數響應時間序列，如圖11 所示。

在6 個參考品種中，濟麥22 和百農207 病害嚴重度等級為“S”，偃展4110、周麥18、西農511 的嚴重度等級為“MR”，周麥36 的嚴重度等級為“R”。在5 月12 日，6 個指數基本能夠區分3 類病害嚴重度等級，其中DVIRE、NDVI、NDVIrededge、NPCI 和WI 響應的時間序列中，效果更為明顯；然而在上述5 個指數中NPCI與WI 對條銹病發病階段的分類并不敏感；對發病階段分類貢獻率較大的DVIRE、GVI、NDVI 和NDVIrededge中，GVI 對病害嚴重度等級也不夠敏感；而DVIRE、NDVI和NDVIrededge 同時對發病階段和病害嚴重度的分類有較大的貢獻。在條銹病脅迫早期，小麥冠層DVIRE 的響應最快，其他植被指數響應較慢，GVI、NPCI 和WI、在條銹病脅迫早期出現了與最終結果相反方向的變化。

如圖11a 所示，從DVIRE 的響應時間序列來看，在小麥條銹病脅迫早期，發病速率較快的是偃展4110，而不是最終病害嚴重度等級最高的百農207；偃展4110 在4 月29 日之后發病速率存在明顯的降低，并最終處于“MR”級別，而百農207 的DVIRE 的響應基本呈線性趨勢下降；而最終嚴重度等級最低的周麥18 在4 月26 日之后就已經發生類似于偃展4110 的變化，發病速率明顯減緩，并在5 月12 日人工鑒定時呈現出“R”等級。

3 結 論

本文提出了一種群體小麥條銹病發病動態無人機遙感分析方法，并通過所建立的機器學習模型驗證了其可靠性，主要結論如下：

1）利用22 個原始特征構建的隨機森林分類模型與通過隨機蛙跳篩選特征建立的SVM分類模型對群體小麥田間條銹病發病階段的分類都極為準確，證明了利用無人機遙感采集并提取的光譜植被指數能夠描述群體小麥田間條銹病發病隨時間的動態變化；篩選出的NDVIrededge、DVIRE 和GVI 對模型分類結果更加敏感，可以更好地實現對群體小麥田間條銹病發病階段的分類。

2）利用22 個植被指數的響應構建的隨機森林分類模型與通過隨機蛙跳篩選特征建立的SVM分類模型都能夠較好的完成對條銹病嚴重度等級的分類，證明了無人機采集并提取的光譜植被指數可以描述小麥群體中不同品種（系）之間的田間條銹病嚴重度差異；篩選出的WI、NPCI、NDVI 在兩種特征選擇算法中取得相似的重要性，對條銹病病害嚴重度更為敏感。

3）通過多時相的光譜指數響應的時間序列能夠在病害發生的不同階段，精細、客觀地量化并分析育種群體中不同品種（系）的發病狀態，實現對群體小麥條銹病表型的動態分析。在全部提取的植被指數中，DVIRE、NDVI 和NDVIrededge 可以同時描述群體小麥條銹病發病階段和病害嚴重度，其中小麥冠層DVIRE 可以作為最佳指標精細量化群體小麥條銹病發病動態情況。