孔繁昌,劉煥軍,,于滋洋,孟祥添,韓 雨,,張新樂※,宋少忠,羅 沖
(1.東北農業大學公共管理與法學院,哈爾濱 150030;2.中國科學院東北地理與農業生態研究所,長春 130012;3.吉林工程技術師范學院信息工程學院,長春 130052)
水稻(OryzasativaL.)作為中國重要的糧食作物,種植面積穩定在3 000萬hm2以上,稻谷產出量占全國糧食總產量的1/4左右,因此水稻在中國的糧食結構中發揮著至關重要的作用。稻瘟病是水稻生產過程中最為嚴重的病害之一,尤以穗頸瘟因其傳播速度快,防控難度大,破壞性強,對產量影響最大[1]。隨著精準農業的不斷推進,對農情預測預報、病害防治的時效和準確性提出了提出更高層次的需求,傳統的“以點代面”的稻瘟病監測手段已不能滿足其要求。
高光譜成像技術是結合了高光譜豐富的波段信息與計算機成像技術二者優勢的一種新的技術,在農產品品質檢測[2-4],農作物長勢監測[5-6],病蟲害的識別[7-8]等方面廣泛應用。高光譜成像技術作為一種無損、精準、高效的檢測技術手段已經在稻瘟病的研究中發揮了舉足輕重作用,Huang等[9]提出了“光譜詞袋”模型分析方法,采用卡方-支持向量機算法,實現了稻瘟病程度的分級模型;袁建清等[10]篩選不同的光譜輸入量,使用偏最小二乘判別分析和主成分加支持向量機的方法構建起水稻葉瘟病與缺氮葉片的判別與診斷模型;黃雙萍等[11]以田間獲取的離體穗株為試驗研究樣本,測定每株稻穗的高光譜影像,基于變種深度卷積神經網絡方法,結合植保專家分級意見建立起了水稻穗瘟的檢測模型;Zhang等[12]提出了一種感染稻瘟病葉片與健康葉片光譜反射率比的數據構建方法,建立起支持向量機水稻葉瘟評價模型,對水稻營養生長后期的水稻葉瘟實現分級檢測。在已有的對稻瘟病成像光譜反演研究中,多數處在實驗室理論階段,尚沒有考慮實際的生產過程中會存在著諸多干擾因素對稻瘟病識別的影響。在其他病蟲害高光譜研究中,部分學者已經開始嘗試在室外環境條件下進行觀測,喬洪波等[13]使用近地成像光譜儀對小麥全蝕病等級展開監測,證明了應用的可行性;Li等[14]使用無人機高光譜影像數據對柑橘的黃龍病進行識別,為黃龍病的大面積監測提供技術基礎。因此,實驗室條件下的稻瘟病高光譜識別已經具備一定的基礎,但對于稻瘟病的外業監測尚需進行更深入的研究[15]。
無人機高光譜遙感不僅能夠實現更大范圍內、更高空間分辨率的病蟲害精準監測,而且能夠快速完成田塊尺度下目標信息的傳遞,獲得目標地物與周圍環境背景的相互關系。本研究采用無人機高光譜遙感數據,分別以光譜預處理變換數據、光譜特征參數(Spectral Characteristic Parameters,SCPs)、植被指數組合(Combination ofVegetation Indices,CVIs)作為輸入量,通過隨機森林(Random Forest,RF)分類算法構建水稻穗頸瘟判別模型。結合最佳輸入量,實現病害等級反演,并解釋穗頸瘟發病期間水稻的生理變化與光譜特征關系,以期為無人機高光譜遙感實現穗頸瘟病定量遙感監測與預警分級提供支持。
本研究數據來自吉林市永吉縣開展的水稻稻瘟病試驗研究,試驗地位于萬昌鎮(43°44′49″N,125°54′11″E),研究區地理位置及無人機高光譜影像真彩色合成圖如圖1所示。該地區位于長白山向松遼平原過渡地帶,屬于溫帶大陸性干寒季風型氣候,夏季雨量充沛,適宜水稻種植,農業熟制為一年一熟。選取的水稻品種為吉玉粳(代號吉90-g4),在每個試驗小區中央插播蒙古稻(易感病,對稻瘟病抗性很差)作為誘發株,最大限度模擬自然發病,并設有不接種對照組,試驗區總面積約5 000 m2。于2019年4—9月開展水稻稻瘟病試驗,2019年 7月 9日吉玉粳進入拔節期,將稻梨孢(Pyricularia oryzae)真菌原液經稀釋后均勻噴灑到蒙古稻葉片表面,稻梨孢侵入蒙古稻后形成誘發株。2019年8月15日,水稻生長進入乳熟期,吉玉粳出現穗頸瘟,隨即請植保專家對不同的采樣區域(1 m×1 m)按照國家標準[16]并結合實地條件對該區域發病程度進行連續的觀測、定級,按照發病程度不同分為無病、輕度、中度、重度4個等級;采樣定級區域共計30個(表1)。

圖1 研究區地理位置及地面采樣區域Fig.1 Geographical location of the study area and ground sampling area

表1 采樣區穗頸瘟病害程度等級描述Table 1 Description of panicle blast disease levels in the sampling areas
1.2.1 無人機高光譜影像采集系統
試驗采用的無人機遙感平臺主要包含:大疆 M600 Pro六旋翼高性能無人機(最大起飛質量6 kg;續航時間16~20 min)、三軸增穩云臺、地面控制站。無人機高光譜成像系統包括Cubert S185機載成像光譜儀和超微型計算機,其中Cubert S185是德國Cubert公司生產的一款全畫幅式機載成像高光譜儀,光譜采集范圍 450~950 nm,光譜采樣間隔為4 nm,擁有125通道。Cubert S185在執行任務時,能夠在0.001 s內獲取1幅1 000×1 000(像素)的全色影像和1幅50×50(像素)的高光譜圖像立方體,同時無人機內置的定位定向系統(Position and Orientation System,POS)記錄每張影像的經緯度信息。
1.2.2 無人機高光譜影像的采集與拼接
無人機高光譜影像的采集時間為2019年8月20日,無人機高光譜遙感平臺對地觀測時間段為10:00—12:00,一次直接獲取整個試驗區所有的正射高光譜影像,為盡可能地避免光線和風速的影響,無人機采集數據時選擇天氣晴朗無云、風速低于5級、控制飛行時間在40 min內進行;起飛前對設備進行白板、暗電流校正;飛行旁向重疊率75%,航向重疊率90%,飛行高度50 m,對應空間分辨率為 3.0 cm;地面采樣、觀測和定級工作同步進行。
無人機高光譜遙感平臺采集的是多幅單景影像,需要對影像進行拼接。采用設備供應商提供的拼接方案,使用光譜儀自帶軟件Cube-Pilot 1.5.5對高光譜影像進行反射率影像的輸出;將已經輸出的全色影像和 POS點數據導入Agisoft PhotoScan 1.5.0 軟件[17-18],Agisoft PhotoScan 1.5.0軟件會根據 POS點位置并結合臨近影像的同名點對影像進行匹配,生成三維點云數據;最后對全色影像與高光譜立方體影像進行影像融合與正射校正,得到本研究區包含地理坐標的無人機正射高光譜影像[19]。
使用專業遙感影像處理軟件ENVI 5.3對影像進行幾何校正,校正后影像能夠覆蓋整個試驗區。隨后采用軟件中的感興趣區(Region Of Interest,ROI)工具,按照地面采樣定級時使用載波相位差分系統獲取的精確地理坐標,將每個對應的采樣區按照30×30(像素)大小矩形區域提取1個ROI,實際對應地面空間分辨率為0.9 m×0.9 m,與實地采樣大小較為吻合。將該區域的所有像素進行計算,取區域像素點反射率平均值作為該采樣區反射光譜。
對ROI的光譜進行9點平滑、一階微分、主成分分析(Principal Component Analysis,PCA)、連續統去除(Continuum Removal,CR)等光譜預處理變換,植被指數計算、組合和光譜特征參數的提取。選取的CVIs是通過文獻查閱[20],選擇與植物感染病害后能夠表征其生理參數發生變化的相關植被指數,并結合本次實際采樣點數據進行相關性計算得出。該組合由反映葉綠素B和葉綠素A比值變化的葉綠素指數(Chlorophyll Index,CI)、作物氮素水平變化的作物氮反映指數(Nitrogen Reflectance Index,NRI)、生物量變化的歸一化植被指數(Normalized Difference Vegetation Index,NDVI)、Vogelman 紅邊指數(Vogelman red edge index,VOG)、冠層溫度變化的光化學植被指數(Photochemical Reflectance Index,PRI)組成,其數學表達式如表2所示。

表2 植被指數組合及數學表達式Table 2 Combination of vegetation indices and mathematical expression
PCA是通過正交變換將較多的變量轉換成一組線性不相關的變量,利用降維的思想對數據量進行壓縮,高光譜數據由于波段多,數據量大,常采用PCA進行處理。
在植被遙感應用中,紅邊一般定義為 670~760 nm的波段范圍,此范圍處在由紅光波段向近紅外波段過渡,植被光譜反射率快速增高的銜接處,當植物受到病害脅迫或葉綠素發生變化等因素影響時,紅邊往往會產生向短波方向移動的“藍移”現象[20]。此外,在450~950 nm整條原始光譜曲線范圍內會有多處特殊的反射、吸收位置,常把450~510 nm(藍光范圍內)的強吸收谷稱為藍谷,520~580 nm(綠光范圍內)的光譜強反射峰稱為綠峰,610~700 nm的(紅光范圍內)的光譜強吸收谷稱為紅谷。
連續統去除(CR)又被稱為去包絡線處理,其最早是用于礦石光譜化學分析,主要的作用是去除環境背景值的影響,有效突出目標物體光譜的差異性,Kokaly等[26]在對植物的氮元素、木質素在葉片含量的分析中首次將該方法引入到植被光譜中。植被在不同的生理狀態下,光譜曲線在某一波段范圍強吸收或強反射會形成吸收谷(反射峰),光譜曲線兩點間斜率變化等是識別植物理化成分變化的重要理論依據。通過分析 CR處理后光譜曲線,提升了稻瘟病在不同光譜范圍內的反射(吸收)差異,提取更多的有關更多病害生理變化特征參數。
本研究使用CR處理后光譜曲線計算SCPs包括谷深(峰值)、谷(峰)面積、光譜曲線波段間斜率。將 CR處理后整條光譜曲線中的450~530 nm波段范圍內形成的吸收谷最大深度(1減去某一對應波長去包絡線后曲線值)記為H1;500~670 nm波段范圍內形成的反射峰最大值(某一對應波長去包絡線后曲線值)記為H2;530~760 nm波段范圍內形成的吸收谷最大深度記為 H3。將450-530 nm波段范圍的連續統去除后光譜曲線與連續統去除最大反射值所在直線圍成的部分進行積分,所得到的吸收谷面積記為BA,以H1為界限,左半部分積分面積記為BA1,右半部分記為BA2;500~670 nm波段范圍內連續統去除光譜曲線與連續統去除最小反射率值所在直線圍成的部分進行積分,所得到的反射峰面積記為GA,以H2為界限,左半部分積分面積記為GA1,右半部分記為GA2;530~760 nm波段范圍連續統去除光譜曲線與連續統去除反射率最大值所在直線圍成的部分進行積分,所得到的吸收谷面積記為RA,以H3為界限,左半部分積分面積記為RA1,右半部分記為RA2。500~530 nm范圍內連續統去除光譜曲線兩端點間連線斜率記為S1;530~670 nm范圍內連續統去除光譜曲線兩端點間連線斜率記為S2;670~750 nm范圍連續統去除光譜曲線兩端點間連線斜率記為 S3。本研究以上操作均在函數繪圖軟件Origin 2017中批量完成,使用以上提取的參數構成稻瘟病識別的SCPs輸入量(圖2)。
采用RF的方法進行建模,RF是一種統計學習理論,具有很高的預測準確率,因其優越的性能,可以摒除大部分噪聲,避免過擬合問題[27]。將不同輸入量代入RF算法進行建模,構建水稻穗頸瘟的分類模型,并應用預測集生成混淆矩陣(Confusion Matrix,CM),CM中針對整個模型準確率(Accuracy,ACC,%)的計算如式(1)所示:

式中TC為正確分類的樣本數;FC為錯誤分類的樣本數。

圖2 連續統去除后光譜特征參數Fig.2 Spectral characteristic parameters after continuum removal
Kappa系數(K)用于對精度進行判斷,可以避免樣本數量不均衡帶來的“偏向性”,其計算如式(2)所示:

式中N為CM中所有觀測樣本數量,r為CM中列的數量,Xii為CM中對角線上的觀測樣本數量,Xi+為CM中第i行的觀測樣本數量,X+i為第i列的觀測樣本數量。
本研究數據集共計30個樣本,20個樣本用于建模,10個用于驗證。數據預處理和模型構建分別在 EXCEL 2010、Origin 2017、ENVI 5.3、RStudio 1.2.5軟件中實現和完成。
在大田環境下獲取到的高光譜影像中提取不同發病區域的反射光譜曲線,隨后對其進行9點平滑處理,9點平滑可以很好地去除光譜曲線隨機噪聲,使曲線變化更加清晰。選取其中具有代表性的健康、輕度、中度和重度等級樣本區域繪制光譜反射曲線(圖3)。由圖3a可知,水稻感染穗頸瘟后,從健康到重病,光譜整體反射率呈現依次下降的趨勢,其中出現 5處明顯的差異性:在530 nm的綠峰波段處,綠峰反射率隨病害加深逐漸降低;在610~700 nm紅谷波段范圍內,隨著病害加深,形成重度、中度、輕度、健康依次由小到大明顯的吸收差異;綠峰與紅谷受病害程度的加深表現出相反的變化,在二者變化共同的作用下,550~670 nm范圍光譜曲線斜率的絕對值逐漸減小;在植物特有紅邊范圍內,由于穗頸瘟病害的脅迫出現紅邊“藍移”現象;在近紅外高反射區域,健康到重度病害反射率逐漸降低的趨勢明顯,而且下降的幅度在不斷的減小。
由圖3b可知,對于野外獲取的水稻光譜數據進行CR處理,水稻反射率變化趨勢進一步拉大。從健康到重病的反射率依次呈現上升趨勢,在466~750 nm之間出現了3個明顯反射吸收的極值,分別以498、534和666 nm為中心點,分別對應藍光波段的吸收谷,綠光波段的反射峰和紅光波段的吸收谷,隨著病害等級升高,極值處出現增加或者減少10%~20%的反射率。因此,計算每一個吸收谷和反射峰圍成的面積,吸收谷深深度和反射峰值,2個極值之間的斜率等光譜曲線特征參數可以很好地反映光譜的差異變化。但在近紅外高反射區光譜變化差異變得不再顯著,這與原始光譜曲線表現出的差異有明顯不同。

圖3 大田環境下不同病害程度的水稻光譜曲線Fig.3 Rice spectrum curves of different disease levels in field environment
使用光譜預處理變換數據、SCPs、CVIs分別作為輸入量代入RF分類模型中(表3),其中基于CVIs的建模方法精度最高,建模精度達到了85%,驗證精度為90%,Kappa系數為0.86,驗證集只有一個重度病害樣本被分到了中度樣本里面,Kappa系數也達到了幾乎完全一致水平,說明CVIs作為輸入量建??梢赃_到對穗頸瘟早期識別達到了理想的效果,針對于處于發病后期的穗頸瘟識別模型需要加強;9點平滑、光譜曲線特征參數、去包絡線也取得不錯的效果,建模精度范圍在70%~80%之間,驗證精度范圍在 70%~80%之間,Kappa系數在 0.59~0.70之間,主要問題表現在將健康、輕度病害錯分為重度病害,模型加重了病害程度分級,相對于CVIs的精度來看,使用這些輸入量建模仍有待優化;PCA作為一種常用的光譜處理方法作為輸入量建模在本研究中建模精度僅為60%,驗證精度為70%,Kappa系數為0.60,并未表現出較好的效果,這可能由于主成分在進行正交變換時,將無用的地物光譜信息混入了正交向量中,導致精度較實驗室理論結果下降;而對于一階微分來看,精度明顯偏低,驗證精度僅有60%,Kappa系數也僅有0.47,這種數據處理方法沒有凸顯穗頸瘟在特征波段的光譜變化,結合RF進行建模精度也就最低。

表3 不同輸入量的建模精度Table 3 Modeling accuracy of different inputs
目前,針對于稻瘟病的分級,使用PCA后數據作為輸入量建模是一種常用的方式。由于本次試驗為無人機平臺下的大田試驗,在室外受到光線以及周圍環境背景的影響,光譜數據進行主成分時可能忽略了有效的光譜信息和夾雜了環境背景噪聲,這導致本次研究結果與實驗室條件下結果有不同。同時,采用主成分后的光譜數據缺乏相應的物理意義的解釋。
稻瘟病的發生和發展帶來一系列水稻生理參數的變化,這些生理參數的變化無論是在光譜曲線還是水稻生長表征都能夠體現出來。在光譜曲線上,稻瘟病等級與光譜曲線的最大正相關出現在 666~670 nm 處,這與Kobayashi等[28]對日本水稻穗頸瘟研究得出的敏感波段范圍結果是一致的,此處屬于紅谷吸收波段,與葉綠素A含量有關;當作物受到脅迫或者接近成熟時,葉綠素 A的含量會出現下降,從而導致植物細胞傳遞光能的能力下降,穗頸細胞遭受到稻梨孢侵染后出現發黃、干枯現象與此相符合。最大負相關出現在750~754 nm處,該處雖然不是葉綠素的強吸收帶,但是隨著葉綠素的減少,吸收帶變窄,引起紅邊“藍移”,因此,750 nm處反射率仍受葉綠素含量的間接影響,同時紅邊位置與水分弱吸收帶(720~740 nm)臨近,紅邊位置變化與穗頸細胞間持水率降低導致的枯黃也有關聯。
CR處理后的光譜曲線,在466~534 nm的位置出現一個明顯的吸收谷,這是較之前原始光譜曲線新出現的特征,計算該處的谷面積、谷深度后發現與稻瘟病病害等級相關性較強。結合水稻生理變化來看,此處與葉綠素 A、類胡蘿卜素強吸收有關系,當穗頸瘟發生時,穗頸細胞組織結構病變,葉綠素 A、類胡蘿卜素降低,因此在此處的吸收谷逐漸收縮。同時,將CR處理后光譜谷深、谷面積、斜率等構建一個能夠反映稻瘟病脅迫條件下水稻植株葉綠素、類胡蘿卜素、細胞結構、紅邊變化等植物生理參數組合作為輸入量進行建模,得到較好的效果。
不同波段組合、變換構建的植被指數可以從多角度反映當作物生長受到脅迫時所發生的生理變化。在已有的病蟲害研究中,很多學者已經指出病害能夠帶來植物本身葉綠素、LAI、氮、作物冠層溫度等生理指標改變。因此,可以將植被指數作為水稻穗頸瘟發病程度與生理指標變化之間的“紐帶”,構建起CVIs來對稻瘟病進行分級。本研究采用的CI由640和674 nm波段組成,分別代表葉綠素B和葉綠素A的吸收峰,穗頸受到脅迫后對藍紫光和紅光吸收率出現差異,從而導致能夠反映穗頸葉綠素B和葉綠素A相對含量變化的CI發生變化;當水稻抽穗期后,水稻稻穗和穗頸部分占據了水稻冠層的優勢,穗頸瘟發生時,穗頸營養運輸發生障礙,大量營養無法正常運輸到稻穗部分,因而導致稻穗內部的氮元素含量較正常稻穗不足,選擇能夠反映氮變化的NRI作為參考可以考察發生穗頸瘟時水稻植株表層氮含量;在530 nm處存在葉黃質敏感的波段,將其與570 nm處反射率組合可以構成冠層光化學指數PRI,該指數指出當作物光合功能出現問題,更多的光能會轉換成熱量散失,因此與冠層溫度密不可分,由于病害發生程度的差異,光合功能受到干擾程度不同,能夠使用該指數來反演水稻冠層溫度;植物受到病害脅迫時,使近紅外反射率下降而可見光部分亮度增高,會出現紅邊“藍移”現象,VOG紅邊指數可以判斷紅邊的變化,借助該指數輔以判斷穗頸瘟發生時紅邊的移動以及紅邊幅值的變化;NDVI作為預測生物量的指標已被廣泛采用,穗頸瘟的發生帶來生物量的改變,在對植被指數與生物量相關性分析過程中,NDVI與生物量的相關性在置信度為 0.01水平下相關性達到了0.699,可以用來估算生物量的變化??偨Y以上選擇特定指數所關注因素、用途,形成了一個特定的CVIs,這個指數組合綜合考慮了水稻穗頸瘟發生不同程度的生理參數的變化,并且成為判斷稻瘟病發展程度的中間參量,建模時獲得了最高的精度,說明具有建模的可行性。
水稻穗頸瘟病室內的識別與分級已有很多學者進行了相關的研究,這類試驗大多都是采樣離體帶回到實驗室,在相對穩定條件下進行試驗,取得了比較高的精度,采用的各種方法均達到 90%以上。但是,在水稻實際生產中,這種方式不僅破壞作物本身的生理結構,而且沒有擺脫“以點代面”的困境。無人機遙感平臺以實際應用為出發點,在室外大田環境下進行數據采集,環境中的光線、隨機噪聲、地物間反射的相互作用等干擾因素不可避免。這導致大田數據建模驗證精度較室內數據低,如何進一步降低數據的噪聲,選擇合適的尺度數據作為輸入量提高建模識別精度應為主要的研究方向,可以為大尺度稻瘟病的監測提供思路、為后續研究開展儲備相應的關鍵技術。
1)分析4類不同穗頸瘟受災程度葉片的反射光譜曲線規律,在反射光譜曲線上有 3個明顯特征:在近紅外波段,隨著葉片感染程度加深,反射率明顯呈現下降趨勢;在550~670 nm之間曲線,健康樣本呈現一個明顯的斜率,而染病樣本斜率趨近于0;670~760 nm的紅邊范圍隨著病害等級升高出現“藍移”現象。連續統去除(Continuum Removal,CR)處理可以增強光譜在466~750 nm的差異性,提升不同病害等級之間可辨識性,從而更有效提取光譜曲線特征參數(Spectral Characteristic Parameters,SCPs)。
2)由不同波段組合變換構成的植被指數能夠反演出水稻在穗頸瘟脅迫條件下的不同生理參數變化,組合生成一個特定的植被指數組合(Combination of Vegetation Indices,CVIs),搭建生理參數與穗頸瘟發病程度之間的“紐帶”,提高了模型準確性。
3)使用隨機森林(Random Forest,RF)算法對穗頸瘟分類建模具有可行性,在所有輸入參數中,使用多種植被指數構成的CVIs作為輸入參量的建模驗證精度最高為90%,Kappa系數為0.86,建模結果可以用來識別大田穗頸瘟。
相對于在室內環境下對稻瘟病的光譜理論研究,無人機高光譜遙感監測是理論研究走向大田應用的關鍵一環,本研究可以為水稻抽穗期發生在穗頸節部位的穗頸瘟病害的大范圍監測、及時指導農戶精準施肥施藥、減少產量損失提供新的支持。