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

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

0 引 言

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

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

1 材料和方法

1.1 試驗區(qū)概況及試驗材料

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

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

1.2 數(shù)據(jù)采集

1.2.1 無人機成像多光譜數(shù)據(jù)獲取

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

設(shè)B層中與Aj相關(guān)的因素成對比較判斷矩陣經(jīng)過了一致性檢驗，求得單排序一致性指標為CI(j)，(j=1，…，m)，相應(yīng)的，平均隨機一致性指標為RI(j)，CI(j)、RI(j)已經(jīng)在層次單排序時求得，那么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 像素；成像系統(tǒng)下行光照傳感器（Downwelling Light Sensor，DLS）模塊用于數(shù)據(jù)采集時環(huán)境光線與太陽角度的校正；平臺集成的全球定位系統(tǒng)（Global Positioning system，GPS）用于記錄圖像的位置信息；反射率校正板用于校正各個通道的反射率。使用該平臺在2020 年3 月31 日至5月14 日采集15 次小麥冠層成像多光譜數(shù)據(jù)，數(shù)據(jù)采集時間為能見度良好的10:00－15:00，無人機飛行高度為20～21 m，前向與側(cè)向重疊率為80%～90%，飛行速度保持在1.4～1.6 m/s。

1.2.2 小麥條銹病人工鑒定數(shù)據(jù)采集

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

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

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

1.3 研究方法

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

1.3.1 無人機圖像預(yù)處理與植被指數(shù)提取

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

1.3.2 光譜響應(yīng)分析

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

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

1.3.3 特征篩選與分類算法

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

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

2 結(jié)果與討論

2.1 群體發(fā)育階段的動態(tài)分析

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

2.2 條銹病發(fā)病階段分類

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

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

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

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

2.3 條銹病嚴重度等級分類

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

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

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

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

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

2.4 條銹病脅迫下光譜響應(yīng)動態(tài)分析

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

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

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

3 結(jié) 論

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

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

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

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