金 秀,齊海軍,李紹穩*
(安徽農業大學信息與計算機學院,安徽 合肥 230036)
馬鈴薯是全球第四大重要的糧食作物,其地位僅次于小麥、水稻和玉米。儲藏期中馬鈴薯干腐病是最常見的真菌類病害,病害的發生會嚴重影響馬鈴薯品質,降低其食用價值和經濟價值[1-5]。干腐病的誘因主要是鐮刀真菌侵染,當真菌孢子侵染到馬鈴薯外表后,其不斷繁殖使得馬鈴薯外部和內部發生病變。植物病理學上將病害過程分為侵入期、潛育期和發病期,當病害發展到發病期階段時,病癥從隱性變成顯性,危害性遞增[6-11],所以病害發病期診斷對于病害防治和預測是非常重要的。
病害監測方法主要有圖像識別[12-17]、環境因子的預測[18-21]、機器視覺[22-23]等,圖像技術在病害病癥顯現時(即病害發病期)有效,但對隱藏類病癥尤其在病害潛育期階段無法診斷。鐮刀真菌孢子侵染馬鈴薯時,馬鈴薯內部色素含量、水分和細胞間隙首先發生變化,光合作用和養分水分吸收、運輸、轉化等生理機能衰退,隨后發病部位出現相應的病癥,馬鈴薯內部成分形態和結構發生改變。因此通過對比不同病害階段的光譜曲線差別,可以判斷染病作物內部結構變化情況。大量研究表明,在光譜的紅外和近紅外波段范圍內,染病作物曲線會發生明顯變化[24-25]。高光譜技術由于其光譜分辨率高,與普通光譜對比,其觀察到的對象特征更加豐富;因此可對病害的病癥外觀進行監測,從而也可對病害造成的內部變化進行分析[26-29]。在病害潛育期時可利用光譜對病斑進行診斷,在病害發病期時可利用圖像進行診斷,所以基于時序高光譜成像的數據分析,可獲得潛育期到發病期時序關鍵點,即診斷馬鈴薯干腐病發病期。
近年來,國內外專家學者使用了多種技術來對病害進行檢測、分析。路陽等[30]利用圖像處理技術對水稻的病害進行智能識別,通過利用灰度值、邊緣檢測、幾何特征提取等方法進行識別。此外,近紅外高光譜技術也被大量用于馬鈴薯檢測上[31]。Dacal-Nieto等[32]應用高光譜技術檢測馬鈴薯黑星病,采用支持向量機(support vector machine,SVM)和隨機分類對映射在馬鈴薯塊莖中的每個像素進行分類,準確度達97.1%,同時利用高光譜技術對馬鈴薯的空心病進行了檢測,其病癥檢測的正確率達89.1%。周竹等[33]利用高光譜技術對馬鈴薯的損傷進行檢測和比較,識別馬鈴薯損傷的準確率高達97.39%。蘇文浩等[34]利用近紅外高光譜技術選擇了7 個特征波長,通過第二主成分識別馬鈴薯外部缺陷,正確率達97.08%。黃濤等[35]針對馬鈴薯的空心病利用近紅外高光譜技術提取8 個光譜特征,再基于SVM算法建立的模型識別率達100%。謝傳奇等[36]利用高光譜中的光譜和紋理特征對番茄早疫病進行早期檢測,通過主成分分析、連續投影算法結合最小二乘SVM構建了早期鑒別模型,其識別率達95%左右。Ray等[37]使用了高光譜技術對馬鈴薯葉片的晚疫病進行研究分析,發現了病葉片在710、720 nm和750 nm特征波長處具有明顯特征。Sugiura等[38]使用RGB成像系統對馬鈴薯早疫病進行評估,發現病害的發生與時序圖像具有一定的相關性,其相關系數為0.77。胡耀華等[39]利用高光譜技術對馬鈴薯葉片的晚疫病進行了診斷研究,其建立的最小二乘SVM的識別模型對診斷的準確率高達94.87%,結果表明高光譜技術可以實現晚疫病脅迫下馬鈴薯病害程度的有效區分。López-Maestresalas等[40]利用短波近紅外高光譜對馬鈴薯表皮的損傷進行無損檢測,其模型的準確率可達98.56%。從國內外的研究可知,高光譜技術既可利用圖像進行病癥分析,又可通過特征波長來進行病害診斷;因此對于RGB圖像等病害識別技術方法而言,高光譜技術更能利用其圖像和光譜特征的融合,充分地采集病害發生變化的信息特征,從而精確地進行病害診斷和識別。
因此,圖譜合一的高光譜技術因其具有大量數據信息而被廣泛用于成分分析中,其中在病害診斷應用上的研究較多;但國內外研究普遍未將不同時序階段的高光譜進行分析和對比,即病害時序高光譜圖像分析,從而丟失了大量有效的診斷信息,無法進行有效的病害預測和判定。本實驗利用時序高光譜數據對馬鈴薯干腐病發病期進行了分析研究。實驗環境一致性、實驗數據量龐大等問題會導致時序高光譜成像實驗和分析方法復雜化[41-42];而結合時序算法與高光譜數據的交叉研究相對較少,其中主要的研究也僅局限于利用時序高光譜技術對化學成分和病害變化進行監測,僅有部分人做了時序高光譜的基礎研究。Stuckens等[43]使用時序高光譜技術對果園的生長狀態進行檢測。Wu Di等[44]利用時序高光譜技術分析了牛肉中水分狀態的變化。Xie Chuanqi等[45]應用時序高光譜技術監測了茶葉水分含量變化。Kim等[46]利用時序高光譜技術分析了松材線蟲病的發生過程。時序高光譜在數據結構上將時序信息與高光譜信息相結合,使得觀測對象的信息量遞增,雖然數據處理難度也隨之增加,但其在病害判定和預測上的前景會更加廣闊。
本研究以馬鈴薯干腐病對象,利用時序高光譜技術對不同階段病癥進行分析,再基于時序數據分析算法,判定出馬鈴薯干腐病潛育期到發病期的時序關鍵點,最終實現馬鈴薯干腐病發病期的診斷。基于時序高光譜的發病期初始點判定能為馬鈴薯干腐病的早期預測和機理研究提供參考。
實驗所用試材為安徽地區的沿淮河春馬鈴薯,購自安徽合肥本地超市。
高光譜成像儀器 美國Opto-Knowledge Systems公司。技術參數如下:相機類型為Scientific Complementary Metal Oxide Semiconductor(SCMOS),相機尺寸大小為226 mm×220 mm×113 mm,質量為13.6 kg,波長范圍為400~1 000 nm,光譜分辨率為1.45 nm/pixel,空間分辨率最大為2 048×2 048 像素,掃描角度為-30°~30°。高光譜成像儀由2 個鹵素燈光源和1 臺進行數據接收的筆記本電腦組成,將樣本放置于觀察臺上進行成像信息采集,圖1為實驗中高光譜儀器的組件架構。

圖 1 高光譜成像儀Fig. 1 Hyperspectral imaging instrument
1.3.1 樣品處理
挑選表皮無任何明顯病害、光滑的馬鈴薯,人工進行鐮刀菌接種。先在馬鈴薯表皮上制作2 個小孔,然后鐮刀菌的孢子懸浮液注射入孔內,進行真菌侵染,從而產生馬鈴薯干腐病。為觀察發病的情況,以打孔未接種樣本、未打孔樣本進行參照對比[47-49]。其中為保證馬鈴薯干腐病發生,將制作的樣本放置到人工氣候箱中儲藏,其環境始終保持為24~25 ℃,相對濕度為90%,通過樣本間的對比發現,經過18 d后,接種的馬鈴薯樣本完全發病。本研究在18 d的發病過程中,每天都對所有的接種樣本進行高光譜數據采集,由于部分接種孔并未發病,導致病斑基本無任何變化,即干腐病并未發生;因此最后在30 個馬鈴薯干腐病病斑中除去未有效發病的樣本,優選了18 個有效樣本進行高光譜數據分析,并將其組合為時序高光譜數據集,因此分析的樣本總數為324(18×18=324)。
時序高光譜圖像采集方法是以病原菌接種當天為初始點,每隔24 h對樣本進行高光譜成像實驗,直至完全發病后(18 d)結束實驗。實驗環境中固定了每天采集時間點(10:00~12:00)和儀器光強(50 W鹵鎢燈)、相機垂直角度(90°)、相機到樣本距離(35 cm)等因素,對同一對象進行數據采集,可確保實驗對象數據的時序一致性。實驗采集的數據為時序高光譜數據,即采集信息是所有時序點上高光譜圖像集成而組成的時序高光譜圖像,此數據由三部分組成:光譜、圖像、時間軸,如圖2所示,不同時間點采集的高光譜信息可構成一個完整的時序高光譜數據集。

圖 2 時序高光譜圖像Fig. 2 Time series hyperspectral images
時序高光譜可以充分描述樣本在接種病原菌后的變化過程,而這種變化過程也可分成線性或非線性。當病害變化過程為線性時,可用簡單聚類算法來劃分病害程度;而當病害變化過程為非線性時,必須針對時序樣本的共性特征進行聚類分析,從而總結出樣本在不同時間點上的變化規律。因此時序聚類方法是通過時序點的樣本聚類分析來劃分病害程度,由于其聚類的算法更符合時序類的非線性特征,所以可精確地判斷病害潛育期到發病期的時序關鍵點,從而達到對馬鈴薯進行干腐病發病期判定的目的。
時序高光譜數據采集的一致性是整個實驗中重點需要解決的問題,本次實驗為保證信息采集的一致性,采取了以下措施:1)為阻隔其他光線干擾、保證拍攝環境的一致性,在搭建的黑房中進行實驗;并且在拍攝高光譜圖像前燈光需預熱約10 min,保證燈光光線強度一致再進行光譜數據采集;同時為保證在拍攝過程中樣本的方向一致性,實驗在樣本的背面進行標注定位;2)為保證拍攝過程中對樣本不進行干擾,本次實驗過程要求非常熟練,須快速將樣本在人工氣候箱和樣本臺之間進行轉移,從而避免樣品污染等問題;3)設備進行預熱,保證輻射的光譜強度基本一致;4)使用了專業拍攝的黑布材料作為背景臺,減少反射等干擾。
1.3.2 數據處理
對采集的數據進行分析,關鍵要建立病癥與高光譜數據之間的關系。通過對馬鈴薯干腐病發病過程的觀察可知,在病害的潛育期無法通過圖像特征確定病癥所在位置,因此必須實現病斑感興趣區域(region of interest,ROI)的時序動態匹配。同時為了提高病害時序關鍵點判定的精確度,利用異常樣本剔除、閾值分割、特征提取等方法對從接種到完全發病的病害病斑時序高光譜進行特征提取和驗證。病斑光譜預處理算法中在圖像的ROI提取上使用了最大類間方差法(OTSU)[50-52],其方法是在判決分析或最小二乘原理的基礎上進行計算的,其算法在目標域背景面積相差不大時能有效地進行分割;在異常值剔除中利用常用的概率密度比法;而在特征提取算法中使用了核函數主成分法(kernel principal component analysis,KPCA)[53],其通過核函數映射結果代替原來的數據,從而可以解決高維數據的線性不可分問題。
在樣本分析過程中,本研究使用了模糊聚類(fuzzy C-means,FCM)算法和動態時間彎曲(dynamic time warping,DTW)算法,FCM是基于傳統的C均值聚類使用了模糊算法,它使劃分到同一集合的對象之間相似度最高,是一種柔性的模糊劃分[54-55];但是聚類算法對時序數據的分析效果較差,主要因為在聚類過程中并未考慮時序具有的相似性,所以發病后的病斑光譜與發病前的病斑光譜容易混淆。為解決該問題,本研究利用DTW算法對時序光譜數據進行聚類分析,并進行時序關鍵點提取,該算法能更準確地對時序數據進行聚類分析并得到結果。時序數據具有高維度、有序性、屬性多、噪聲雜等特征,因此可通過不同向量值之間的距離來決定類別歸屬,DTW算法最初是用于語音識別,其優點在于:允許2 個時間序列非等長比較以及不同步的點對應計算,可以對時間序列的同步問題進行處理,將兩個序列的錯位點關聯起來。因此DTW算法(式(1)~(3))可對時序序列進行劃分[56-57],其算法定義為:設兩條序列Q={q1,q2,…,qn},C={c1,c2,…,cm},長度分別為n和m,式(1)表示了序列Q和序列C中的一條最短路徑,即DTW距離。

式中:W表示變量i在序列Q上、變量j在序列C上進行遞歸計算的結果,d表示兩個相同長度序列的距離函數(曼哈坦距離、歐幾里得距離等)。
本研究基于DTW的聚類分析,先通過時序光譜之間距離計算不同病斑光譜數據之間的關聯度,然后再通過時序點的模式形態來判定時序關鍵點,圖3為基于DTW的高光譜數據處理流程圖。
在高光譜數據處理過程中,H值代表了DTW的時序點百分比,其中3 個時序鄰接點的H值可被描述為單調序列和非單調序列變化兩種基本的模式[58-59],圖4的a、b為3 個序列點的非單調序列變化,圖4的c、d為3 個序列點的單調序列變化,本研究基于DTW的算法將具有H值非單調序列變化的第一個時序點判定為時序關鍵點,圖4的時序點2、3、4號構成的第一個H值非單調點,因此時序點3號為病害時序關鍵點,即為發病期的初始點。

圖 3 基于DTW的高光譜數據處理流程圖Fig. 3 Flow chart of hyperspectral data processing on the basis of DTW

圖 4 1號樣本的H值Fig. 4 H value of sample 1
本研究中,時序高光譜數據以cub的格式導入到代碼中進行數據批量處理和分析,并實現以上算法。本實驗使用的數據處理平臺為envi 4.2、IDL 6.2、Matlab 7.0。
1.3.3 正確率評價標準
病害發病期初始點的人工診斷方法是以病害表現形態的差異化來定義的,馬鈴薯干腐病從接種真菌開始進行病斑狀態的觀察,當病斑發生擴大或顏色加深時,則判定當天的病斑為發病期時序關鍵點,但人工診斷發病期初始點也具有小范圍的誤差。在本實驗中,因為時序高光譜數據是以1 d為單位進行采集的,為消除人工診斷結果和實驗過程造成的正確率偏差,病害時序關鍵點的判定方法以人工診斷判定為基準點,允許判定結果有前后1 d以內的誤差;但誤差范圍超過時,則代表該病害時序關鍵點判定錯誤。
馬鈴薯干腐病的高光譜圖像數據量大,每個像素點包含339 個波長原始輻射值,因此必須劃分出病斑的ROI。從發病過程可知,病害后期病癥較為明顯,光譜提取較易,而病害前期病癥圖像幾乎無法分辨,無法通過閾值分割進行光譜提取。因此本實驗利用明顯病斑區域進行OTSU閾值分割后,再基于分割后的ROI來映射病害前期的病斑(圖5A),最終可提取出時序高光譜數據(圖5B)。其中病斑高光譜的時序變化并不具有規律性,下一步將對病斑光譜特征進行分析及提取。

圖 5 1號馬鈴薯干腐病的病斑光譜數據集Fig. 5 Spectral data of dry rot lesions for sample 1
為分析樣本高光譜的時序性特點,先對數據進行主成分分析,以第1主成分的兩個波峰538、794 nm作為特征進行時序圖的繪制,發現馬鈴薯干腐病的病斑光譜隨著病害的發生呈現明顯的非單調性(圖6)。因此為匹配病斑時序高光譜的非單調性特點,利用高斯KPCA進行光譜特征的提取,通過前3 個主成分的加權重系數曲線,提取了14 個波峰作為病斑高光譜的特征波長(圖7)。

圖 6 馬鈴薯在儲藏18 d過程中病斑區域于538、794 nm波長處的光譜時序圖Fig. 6 Time series spectral plot of dry rot lesions during 18 days of storage at 538 and 794 nm

圖 7 KPCA的權重系數Fig. 7 Weight coeff i cient of KPCA
通過對馬鈴薯進行18 d連續的高光譜拍攝獲得了時序高光譜圖像,該數據描述了干腐病從接種到發病的整個過程。在實驗的過程中,馬鈴薯干腐病的病斑變化呈現出兩個特性:1)病斑包含了腐爛區域和菌絲區域兩部分,首先是接種位置周邊腐爛,然后腐爛區域加速了真菌的生長,當菌絲擴散到一定范圍時,腐爛區域會繼續變大、加深,這種病斑階段性的轉變導致了光譜反射率呈非線性時序變化;2)干腐病發病困難,與未接種馬鈴薯對比發現,大部分樣本在接種后的第5~6天才發生明顯變化,第18天確定干腐病已產生,切開馬鈴薯表皮進行觀察,其內部明顯有發病癥狀,所以可確定該馬鈴薯進入了干腐病的發病期(圖8)。


圖 8 馬鈴薯1號的干腐病時序圖Fig. 8 Time series images of potato sample 1 infected with dry rot
根據以上的病斑觀察情況,利用人工專家知識可以確定出馬鈴薯干腐病發病期,第4~5天由于馬鈴薯1號病斑出現了明顯的菌絲和腐爛癥狀,基本可以確定該病斑點進入發病期,其時序關鍵點(發病期初始點)即為第4~5天。因此利用人工觀察的病害程度可確定出馬鈴薯干腐病的發病期,其18 個樣本的時序關鍵點見表1。

表 1 人工觀察判定的時序關鍵點Table 1 Time series key points by artif i cial judgment
本研究使用了FCM進行時序關鍵點的判定,該算法是通過隸屬度來確定每個數據點的類別。因為本次實驗目的為判定病斑的發病期,即發病的時序關鍵點,所以將FCM設定為2 類,且將第1天的病斑光譜固定為第1類。通過分析了18 個樣本的18 d光譜數據后,得到的結果見表2。與人工診斷結果(表1)進行對比分析可知,7、12號樣本發病期的時序關鍵點判定完全正確,1、4、11號樣本正確率偏差了0.5 d,2、3、5、8、13、16、18號的樣本正確率偏差了1 d,且其判定的關鍵點在實際發病后,因此,12 個樣本的結果在規定范圍內,故為正確的時序關鍵點;其中,18 個樣本中,第6、9、10、15、17號樣本的判定結果相對提前,雖然其偏差也為1 d,但從病害的時序發病狀態下進行判斷,此類結果也屬于嚴重偏差;同時第14號樣本判定結果離準確值較遠,因此該類結果偏差更為嚴重,所以本次FCM進行時序關鍵點判定的正確率為12/18=66.7%。基于FCM的發病期判定結果可知,可以通過時序高光譜的聚類分析進行時序關鍵點判定病害,但與正確的時序關鍵點相比偏差較大,因此判定精確性和準確度都較差。

表 2 干腐病病斑的時序關鍵點Table 2 Time series key points for dry rot disease lesions

圖 9 基于DTW算法的18 個樣本H值Fig. 9 H values of 18 samples based on DTW
馬鈴薯干腐病的時序高光譜數據通過DTW算法來進行聚類分析,可得到時序點的DTW距離百分比(圖9)。利用算法進行聚類分析后可知,發病期時序關鍵點如表3所示。將基于DTW的時序關鍵點與表1相對比可知:第12號病害樣本出現半天的偏差,第11號病害樣本出現了完全錯誤,所以判定出整體的發病期時序關鍵點判定的正確率為17/18=94.4%,而且其結果僅有1 個病害樣本出現了半天的偏差,其準確度較高。研究結果表明基于DTW距離的聚類分析能判定馬鈴薯干腐病的潛育期與發病期的時序關鍵點,且與FCM聚類相比其判定結果精確性和準確度都有很大的提升。

表 3 基于DTW算法的時序關鍵點判定Table 3 Time series key points based on DTW
馬鈴薯的干腐病病癥發生過程中除了病斑的變化,同時也伴隨鐮刀菌的菌絲生長,因此病害的潛育期相對于發病期更加隱蔽,監測和診斷難度更大,所以發病期時序關鍵點判定非常困難。本實驗基于病癥時序高光譜數據進行模糊聚類分析后發現,FCM的時序關鍵點判定的錯誤率高達33.3%,表明了復雜病癥時序數據情況下,僅用非時序性的聚類分析進行發病期時序關鍵點的判定是具有一定缺陷的。因此,為匹配時序數據的特性,本實驗在時序關鍵點判定上,使用了基于DTW的聚類方法,分析結果發現該算法能準確預判馬鈴薯干腐病發病期的時序關鍵點。其病害的時序關鍵點診斷的結果中,18 個病斑樣本完全錯誤的時序關鍵點和半天偏差的時序關鍵點都僅1 個,整體的正確率高達94.4%,明顯優于FCM聚類結果。研究結果表明,時序高光譜技術可以診斷病害的潛育期到發病期的時序關鍵點,其中基于DTW的時序高光譜聚類方法對病害發病期的診斷準確度較高,其研究結果也為后期的病害預測模型提供了參考。