楊麗麗，王振鵬，吳才聰

中國農業(yè)大學信息與電氣工程學院，北京 100083

引 言

新疆作為我國最重要的產棉基地，保障棉花產業(yè)健康發(fā)展，對促進新疆經濟具有重要意義[1]。病蟲害的大范圍爆發(fā)嚴重影響了棉花的產量和品質，其中棉葉螨因其隱蔽性強、個體微小的特點，導致監(jiān)測難度大，危害程度深[2]。傳統(tǒng)的螨害監(jiān)測方式主要是由當?shù)刂脖Ｈ藛T定期在棉田內憑經驗進行巡查，手工標記棉葉螨受害中心植株位置。這種監(jiān)測方式范圍小、人力消耗大且時效性差，不能準確評估大范圍棉田的受害情況。

國內外學者在作物病蟲害遙感監(jiān)測方面進行了大量研究。蘭玉彬等[3]通過獲取柑橘果園的無人機高光譜影像建立柑橘黃龍病判別模型，實現(xiàn)了低空高光譜遙感對柑橘黃龍病的大面積監(jiān)測。黃煥華等[4]利用小型固定翼無人機同時采集枯死松樹的可見光和近紅外航攝影像，對松材線蟲病的監(jiān)測準確率可達到80%以上。田野[5]利用多時相HJ小衛(wèi)星遙感影像實現(xiàn)棉葉螨害監(jiān)測，但因其光譜分辨率難以達到要求還需結合地面光譜數(shù)據(jù)進一步研究。王守會等[6]基于棉花冠層高光譜數(shù)據(jù)和Landsat8衛(wèi)星遙感數(shù)據(jù)對區(qū)域螨害發(fā)生進行預測研究，模型預測效果較好，但衛(wèi)星遙感影像分辨率過低使得監(jiān)測樣點之間的光譜差異大大降低。崔美娜等[7]基于無人機多光譜影像建立棉田螨害監(jiān)測識別模型，實現(xiàn)區(qū)域范圍內螨害識別，但并未考慮棉田氣象和土壤環(huán)境對螨害的影響。趙亮等[8]通過獲取不同螨害程度葉片的地面高光譜數(shù)據(jù)，篩選出螨害敏感波段，建立了棉葉螨葉片光譜識別的最佳模型。

本研究利用無人機遙感的高時效性、低成本和高時空分辨率，彌補了地面光譜難以實現(xiàn)以點帶面以及衛(wèi)星遙感空間分辨率低、時效性差的問題，并綜合考慮不同環(huán)境因素對于螨害發(fā)生的影響，建立了基于多源數(shù)據(jù)的棉田大范圍螨害監(jiān)測模型和螨害面積預測模型，可在實際應用環(huán)節(jié)完成快速大范圍的監(jiān)測需求，為棉田變量施藥提供依據(jù)，有效減少損失。

1 實驗部分

1.1 研究區(qū)概況

研究區(qū)位于新疆生產建設兵團第八師145團六分場六連，中心位置為86°3′30.85″E，44°32′30.23″N，屬中高緯度的典型大陸性氣候，光照充足，雨量稀少，氣候干燥，冬季嚴寒，夏季嚴熱，研究區(qū)氣候條件適宜棉葉螨正常生長[9]。研究區(qū)棉田面積34.7萬m2，約520畝，種植品種為新陸早166號和新陸早163號，采用1膜3管6行種植模式，蕾期一般處于當?shù)氐?月上中旬至7月上旬。

1.2 數(shù)據(jù)獲取與處理

1.2.1 棉花冠層高光譜數(shù)據(jù)采集與處理

2018年于北京時間7月12日12:00—14:00對蕾期棉花進行不同螨害程度的冠層光譜采集，采集設備為ASD FieldSpec Pro FR2500型背掛式野外高光譜輻射儀，光譜波長范圍為350～2 500 nm，采樣間隔為1 nm，光譜儀視場角為25°。采集時保持探頭垂直向下，光譜儀距離棉花冠層頂垂直高度約為0.5 m，每次測量前后均進行標準白板校正，每個采樣點采集5組光譜數(shù)據(jù)，取平均值作為該采樣點的光譜值。

1.2.2 無人機多光譜數(shù)據(jù)采集與處理

2019年棉花蕾期內，選擇晴朗無云、無風或風力較小的天氣，分別于6月22日、27日、29日和7月3日、9日、13日在研究區(qū)內進行6次多光譜數(shù)據(jù)采集，采集設備為瑞士eBeeSQ固定翼無人機搭載ParrotSequoia五通道多光譜相機。該相機鏡頭焦距為3.98 mm，分辨率為1 280×960，可同時獲取4個多光譜波段和可見光數(shù)據(jù)。4個多光譜波段分別為綠光(中心波長550 nm，波段寬度40 nm)、紅光(660 nm，40 nm)、紅邊光(735 nm，10 nm)和近紅外光(790 nm，40 nm)。

無人機起飛前用Sequoia多光譜相機配備的參考板進行輻射定標處理，根據(jù)預設航線垂直拍攝，影像范圍覆蓋整個研究地塊。設置無人機飛行航線為S型，飛行高度為120 m，飛行速度為13 m·s-1，旁向重疊度為60%，航向重疊度為75%，影像空間分辨率為0.12 m。數(shù)據(jù)采集后，將影像數(shù)據(jù)導入Pix4Dmapper軟件進行輻射定標和拼接處理。

1.2.3 螨害實地調查數(shù)據(jù)采集

在每期無人機多光譜影像采集的時間點附近由人工進行實地葉片螨害抽樣調查。在研究區(qū)棉田內使用VIVO X20手機拍攝健康和螨害葉片圖像，采樣比例1∶1，共135個采樣點。采樣時記錄葉片受害程度并根據(jù)手機自帶的GPS精確定位功能記錄采樣點經緯度信息。螨害等級參照國家標準《棉花葉螨測報技術規(guī)范》進行劃分為4級： 0級為完害； 1級為葉面有零星黃斑； 2級為紅斑占葉面積1/3以下； 3級為紅斑占葉面積1/3以上。4種不同螨害程度棉葉示例圖像如圖1。

圖1 不同螨害程度棉葉圖像(a): 0級; (b): 1級; (c): 2級; (d): 3級Fig.1 Images of cotton leaves with different degree of mite damage(a): Level 0; (b): Level 1; (c): Level 2; (d): Level 3

1.2.4 環(huán)境數(shù)據(jù)采集與處理

無人機多光譜數(shù)據(jù)采集期內，通過棉田研究區(qū)搭建的氣象監(jiān)測系統(tǒng)和土壤墑情設備，獲取2019年6月—2019年7月期間逐日每小時的氣象數(shù)據(jù)和土壤溫濕度數(shù)據(jù)共13種田間環(huán)境數(shù)據(jù)。再以4 d作為一個周期進行統(tǒng)計，計算出9種氣象數(shù)據(jù)，即每周期最高溫度、平均溫度、最低溫度、平均濕度、累計降雨量、平均風力、溫濕系數(shù)、溫雨系數(shù)、積溫(逐日平均溫度≥25 ℃)； 4種土壤數(shù)據(jù)，即10 cm土壤最高溫度、10 cm土壤最低溫度、10 cm土壤平均溫度和10 cm土壤平均濕度。同時考慮到農戶施藥對螨害的影響，在施藥日期(6月24日和7月8日)節(jié)點附近將周期長度細化為2 d。

2 研究方法

首先，通過采集棉花冠層高光譜數(shù)據(jù)及無人機多光譜影像數(shù)據(jù)，篩選敏感波段位置，構建螨害監(jiān)測植被指數(shù)； 其次，通過相關性分析篩選相關性顯著的螨害監(jiān)測植被指數(shù)和田間環(huán)境數(shù)據(jù)，建立基于支持向量機(support vector machine，SVM)的棉田螨害發(fā)生監(jiān)測模型； 最后，基于最優(yōu)模型分類結果生成螨害發(fā)生空間分布圖，通過統(tǒng)計健康點和螨害點的像元數(shù)計算螨害面積占比值，并結合田間環(huán)境數(shù)據(jù)采用多元逐步回歸分析方法建立棉田螨害面積預測模型?？傮w技術路線如圖2所示。

圖2 棉田螨害監(jiān)測技術路線Fig.2 Technical route of cotton field spider mite monitoring

2.1 螨害監(jiān)測植被指數(shù)計算

植被指數(shù)(vegetation index，VI)通過不同波段反射率的線性或非線性組合變化，在一定程度上可以突顯植被的不同特征，削弱背景等因素對植被光譜特征的干擾。工作中根據(jù)地面尺度的棉花冠層高光譜數(shù)據(jù)提取出多個螨害敏感波段，從中選取與低空尺度的無人機多光譜波段位置相對應的4個敏感波段，初步構建23種螨害監(jiān)測植被指數(shù)。植被指數(shù)及其計算公式如表1所示。

表1 植被指數(shù)及其計算公式Table 1 Calculation formulas for vegetation index

2.2 模型建立

通過相關性分析篩選反映螨害最敏感的植被指數(shù)和田間環(huán)境數(shù)據(jù)作為建模特征值，分別用0(無螨害)和1(有螨害)作為目標值，將135個樣本數(shù)據(jù)集輸入SVM分類器進行訓練和分類，采用單一環(huán)境數(shù)據(jù)、單一植被指數(shù)、環(huán)境數(shù)據(jù)與植被指數(shù)相結合的3種方法建立基于SVM棉田螨害發(fā)生監(jiān)測模型。模型訓練和測試使用Matlab R2019b和LIBSVM 3.24工具箱，并選取準確率(Accuracy)、精確率(Precision)、召回率(Recall)和對精確率和準確率整體評價的F1值作為模型精度的評價指標，確定出最優(yōu)螨害發(fā)生監(jiān)測模型。

使用田間環(huán)境數(shù)據(jù)分別與螨害面積占比值作相關性分析，篩選出顯著特征值，再基于多元逐步回歸分析法選出與螨害面積值關系最密切的環(huán)境特征建立棉田螨害面積預測模型，選取決定系數(shù)R2作為模型評價指標。

3 結果與討論

3.1 地面冠層光譜數(shù)據(jù)反射率分析

棉花葉片在遭受不同程度的棉葉螨危害后，棉花冠層光譜反射率曲線呈現(xiàn)規(guī)律性變化，呈現(xiàn)明顯的波峰和波谷等特征。圖3為不同螨害程度的棉花冠層光譜反射率曲線。

圖3 不同螨害程度地面冠層光譜反射率曲線Fig.3 Spectral reflectance curves of ground canopies with different degrees of mite damage

由圖3可見，在綠光波段范圍內553 nm波長附近出現(xiàn)第一個反射波峰后反射率逐漸下降，至紅光波段范圍內680 nm波長附近出現(xiàn)反射波谷-紅谷，此時具有較低的反射率。反射率于705 nm波長處開始呈現(xiàn)加速上升趨勢，且在近紅外波段760～1 350 nm范圍內反射率較高，同時在此范圍內光譜反射率隨著螨害程度的逐漸加重表現(xiàn)出較為明顯的下降趨勢，在中紅外波段1 550～1 800 nm范圍內的光譜反射率也表現(xiàn)出了同樣的規(guī)律。

紅邊是綠色植物最明顯的光譜特征之一。為分析不同螨害程度棉葉冠層光譜的紅邊特征，對紅邊(680～750 nm)范圍內的光譜反射率求一階微分，如圖4所示?？梢钥闯?80～750 nm范圍內的一階微分光譜曲線圖像呈現(xiàn)“雙峰”現(xiàn)象，且隨著螨害加重，反射率一階微分值逐漸減小，表明紅邊波段對不同程度的螨害較為敏感。

圖4 不同螨害程度地面冠層一階微分光譜曲線Fig.4 First derivative spectral curve of ground canopy with the different degrees of spider mites damage

分析結果表明，棉花葉片在遭受不同程度的棉葉螨為害后，可通過棉花冠層光譜反射率及反射率一階微分的變化進行表征，綠光波段553 nm附近、紅光波段680 nm附近、紅邊波段680～750 nm范圍、近紅外波段760～1 350 nm范圍都是與螨害發(fā)生相關的敏感波段，與文獻[8]研究結果一致。

3.2 無人機多光譜數(shù)據(jù)反射率分析

3.1節(jié)中高光譜的螨害敏感波段分析結果顯示，本實驗中無人機攜帶的多光譜相機波段范圍屬于螨害敏感波段，為進一步確定無人機攜帶多光譜相機進行蟲害監(jiān)測的可行性，對無人機多光譜遙感影像進行了不同螨害程度采樣點的反射率提取，圖5為反射率曲線。

由圖5可見，綠光波段550 nm附近、550～660 nm波段范圍、660～735 nm波段范圍處不同等級螨害發(fā)生點的光譜反射率基本重疊在一起，可能的原因是無人機的飛行高度影響影像分辨率，造成不同螨害程度的采樣點在光譜響應上的變化差異較小，光譜差異很難區(qū)分，不能進行受害等級區(qū)分，但有害和無害的光譜反射率差異比較明顯。故利用無人機多光譜遙感進行螨害監(jiān)測模型構建時，只構建基于SVM的二分類模型，將螨害1級及以上均視為螨害發(fā)生，即未發(fā)生螨害記為0，發(fā)生螨害記為1。

圖5 不同螨害程度的多光譜影像反射率曲線Fig.5 Reflectance curve of multispectral image with different degrees of spidermite damage

3.3 螨害監(jiān)測模型特征篩選

將選取的23種植被指數(shù)和預處理后的13種田間環(huán)境數(shù)據(jù)分別與地面采樣點的螨害發(fā)生數(shù)據(jù)做相關性分析，用顯著性(sig)作為評價指標，當sig≤0.05時，植被指數(shù)與螨害發(fā)生達到顯著相關水平，當sig≤0.01時，植被指數(shù)與螨害發(fā)生達到極顯著相關水平。從表2看出，所選的植被指數(shù)中，SAVI，OSAVI，TVI和NDGI與螨害發(fā)生呈極顯著相關，RDVI，RVI和MSR與螨害發(fā)生呈顯著相關。從表3可以看出，平均濕度、溫濕系數(shù)和10 cm土壤平均溫度與螨害發(fā)生呈極顯著相關，最高溫度、平均溫度、積溫、10 cm土壤最高溫度和10 cm土壤平均濕度與螨害發(fā)生呈顯著相關。研究中選取sig≤0.05的植被指數(shù)及環(huán)境特征作為螨害監(jiān)測建模的特征值。

表2 植被指數(shù)與螨害發(fā)生之間相關性Table 2 Correlation between the occurrence of spider mite and vegetation indices

表3 環(huán)境數(shù)據(jù)與螨害發(fā)生之間相關性Table 3 Correlation between the occurrence of spider mite and environmental data

3.4 基于SVM算法的棉田螨害監(jiān)測模型

選取表2及表3中的7種植被指數(shù)和8種田間環(huán)境數(shù)據(jù)作為建模特征值，在135個地面采樣點中，隨機選取出總樣本數(shù)的2/3作為訓練集，其余的1/3作為測試集進行模型構建，分別基于單一環(huán)境數(shù)據(jù)、單一植被指數(shù)、環(huán)境數(shù)據(jù)與植被指數(shù)相結合的3種方法建立SVM棉田螨害發(fā)生監(jiān)測模型M1，M2及M3。SVM模型的精度主要與核函數(shù)和SVM內的參數(shù)選取有關，研究中選擇徑向基函數(shù)(radial basis function，RBF)作為SVM訓練模型的核函數(shù)，使用交叉驗證的方法調整SVM內的懲罰參數(shù)c和gamma參數(shù)，得到各類評價指標效果最優(yōu)的監(jiān)測模型。模型結果如表4所示。

表4 不同螨害監(jiān)測模型分類結果對比Table 4 Comparison of classification results of the different degrees of spider mite monitoring models

表4中，M3效果最好，測試集分類準確率為80%，精確率為75%，召回率為70.59%，F(xiàn)1值為72.73%； M2次之，M1分類效果最差。表明，同時使用環(huán)境數(shù)據(jù)和植被指數(shù)綜合考慮了影響螨害發(fā)生的主要因素，從而使模型具有更高的精度和監(jiān)測效果。

3.5 棉田螨害嚴重度空間分布

基于棉田螨害監(jiān)測模型M3，分別繪制6月22日、27日、29日和7月3日、9日、13日共計6景無人機遙感螨害監(jiān)測空間分布圖[圖6(a—f)]。

圖6 不同時期(6月22日，27日，29日，7月3日，9日，13日)棉田螨害監(jiān)測空間分布圖(a, b, c, d, e, f)Fig.6 Spatial distribution map (a, b, c, d, e, f) of cotton field spider mite monitoring in (June 22, 27, 28, Jule 3, 9, 13) different periods

圖6中可見，6月22日開始，研究區(qū)螨害成片發(fā)生且基本分布在棉田東部以及中路以南區(qū)域[圖6(a)]，棉田實際調查發(fā)現(xiàn)棉田東部及南部區(qū)域靠近路邊，地邊雜草較多且農戶活動多聚集于兩地，同時發(fā)現(xiàn)監(jiān)測期間研究區(qū)中部以南地區(qū)因棉田灌溉設施漏水，棉花長勢較弱，這些因素綜合導致了棉田邊緣易發(fā)生的螨害擴散到了棉田中部及以南。根據(jù)6月22日棉田螨害監(jiān)測情況，農戶于6月24日使用25%阿維乙螨唑懸浮劑對葉螨進行施藥防治，從圖6(b)可以看出，6月27日分布圖中螨害面積逐漸減小，葉螨數(shù)量較之前有所下降。但由于棉田內施藥時間不同及葉螨抗藥性強，觀察6月29日棉田螨害監(jiān)測圖[圖6(c)]發(fā)現(xiàn)棉田北部螨害擴散較其他地區(qū)略微嚴重。

隨著7月新疆地區(qū)溫度逐漸升高，高溫條件對棉葉螨擴散極其有利，同時藥劑對葉螨的防治效果下降，導致7月3日螨害逐漸加重，呈現(xiàn)局部棉田葉螨暴發(fā)狀態(tài)且多集中于東、西部地塊邊緣區(qū)域[圖6(d)]。因此，農戶于7月8日再次對棉葉螨進行施藥防治，但隨著植保人員和農戶在田中不斷移動進行實地調查等操作，導致葉螨隨之遷移擴散，通過圖6(e)可以看出，中部以北區(qū)域出現(xiàn)局部田塊爆發(fā)成災現(xiàn)象，防治難度增加，易由逐漸擴散轉變成迅速擴散。7月13日再次對棉田進行監(jiān)測，發(fā)現(xiàn)螨害中部以北區(qū)域已達到大面積爆發(fā)狀態(tài)，而中部以南區(qū)域螨害防治程度較好，無局部爆發(fā)狀態(tài)[圖6(f)]。

監(jiān)測圖整體趨勢與實際調查情況較為相符，通過基于環(huán)境數(shù)據(jù)與植被指數(shù)相結合建立的SVM棉田螨害監(jiān)測模型可以反映田間螨害變化，基于無人機的低空遙感可以輔助用來實現(xiàn)棉田大范圍的螨害發(fā)生監(jiān)測。

3.6 棉田螨害面積預測模型的建立

首先基于模型分類結果計算螨害面積占比值，考慮到實際應用環(huán)節(jié)，環(huán)境數(shù)據(jù)更易獲取，故只將初選的13種田間環(huán)境數(shù)據(jù)與棉田螨害面積值作相關性分析，選取sig≤0.05的田間環(huán)境數(shù)據(jù)，結果如表5所示。

表5 環(huán)境數(shù)據(jù)與棉田螨害面積之間相關性Table 5 Correlation coefficient between environmental data and cotton field spider mite area

分別使用SVM回歸和多元逐步回歸方法構建螨害面積模型，以決定系數(shù)R2作為指標來衡量模型的性能。在Matlab下調用LIBSVM工具包建立SVM螨害面積預測模型M1，并通過交叉驗證方式尋優(yōu)參數(shù)，得到最終模型R2=0.796。采用SPSS進行多元逐步回歸分析，從表5中逐步選出與螨害面積值關系最密切的田間環(huán)境數(shù)據(jù)作為建立預測模型M2的特征值，最終得到多元逐步回歸方程為

y=110.826-0.193x1-0.089x2-0.5x3+

0.423x4-4.022x5

(1)

式(1)中，x1是最高溫度，x2是平均濕度，x3是溫濕系數(shù)，x4是10 cm土壤最高溫度，x5是10 cm土壤平均溫度，y是預測的螨害面積占比值。模型評價結果見表6。

表6 模型評價結果Table 6 Model evaluation results

由表6可知，多元逐步回歸預測模型的決定系數(shù)R2=0.848，說明該模型擬合度優(yōu)于SVM回歸預測模型，并對多元逐步回歸預測模型進行顯著性檢驗，計算出P值為0.001，小于0.05。

本研究最終使用多元逐步回歸模型對測試樣本進行預測，將螨害面積占比值單位轉換成畝數(shù)，結果如表7所示。

表7 棉田螨害面積預測模型預測結果對比Table 7 Comparison of prediction results of cotton field spider mite area prediction models

結果表明模型對棉田螨害發(fā)生面積的預測效果較好。

4 結 論

對棉田大范圍螨害監(jiān)測和預測方法進行研究，通過地面尺度的棉花冠層高光譜數(shù)據(jù)結合低空尺度下的無人機多光譜遙感影像數(shù)據(jù)進行敏感波段篩選，計算螨害相關植被指數(shù)，并結合試驗區(qū)的氣象和土壤等田間環(huán)境數(shù)據(jù)，使用單一環(huán)境數(shù)據(jù)、單一植被指數(shù)以及環(huán)境數(shù)據(jù)與植被相結合的3種方式建立SVM螨害發(fā)生監(jiān)測模型并進行對比分析。結果表明，同時使用環(huán)境數(shù)據(jù)和植被指數(shù)建立的模型具有更高的準確度?；谝呀⒌腟VM棉田螨害監(jiān)測模型，繪制不同時期的螨害監(jiān)測空間分布圖，并分別統(tǒng)計健康點和螨害點像元數(shù)計算螨害面積占總面積比值，結合田間環(huán)境數(shù)據(jù)建立棉田螨害面積預測模型，模型的決定系數(shù)R2=0.848，可以用來對棉田螨害面積發(fā)生進行預測。

研究所建立的模型基于田間環(huán)境數(shù)據(jù)和植被指數(shù)，這兩個因素本身就具有時間和空間的變異性，對一段時間內的棉田進行監(jiān)測和面積預測，驗證了時間尺度的可移植性，對于區(qū)域上的驗證，需要采集不同區(qū)域的田間環(huán)境數(shù)據(jù)和植被指數(shù)，理論上與時間上的驗證同理。目前只針對棉田是否發(fā)生螨害建立了監(jiān)測模型，對不同螨害等級的發(fā)生監(jiān)測未進行研究，且模型雖然考慮了田間環(huán)境數(shù)據(jù)和施藥防治對螨害的影響，但是越冬基數(shù)、天敵等因素同樣會影響螨害的發(fā)生。本研究結果可為大范圍棉田螨害監(jiān)測和面積預測提供方法參考，后續(xù)研究中應融合更多相關因素，提升螨害監(jiān)測與預測的精度。