基于高光譜圖像波段融合的獼猴桃軟腐病早期分類檢測

高宏盛,郭志強*,曾云流,丁 港,王逍遙,李 黎

1. 武漢理工大學信息工程學院,湖北 武漢 430070 2. 華中農業大學園藝植物生物學教育部重點實驗室,國家柑橘保鮮技術研發專業中心,湖北 武漢 430070 3. 中國科學院植物種質創新與特色農業重點實驗室,中國科學院獼猴桃產業技術工程實驗室,中國科學院種子創新研究院,中國科學院武漢植物園,湖北 武漢 430074

引 言

獼猴桃富含豐富的維生素和礦質元素,被譽為“水果之王”,深受消費者的喜愛。 然而由于獼猴桃果實質地柔軟,皮薄多汁,極易受到病原微生物的侵染而致病。 其中,軟腐病是目前獼猴桃采后最為嚴重的真菌病害,感病果實在貯藏前期完全無癥狀,在貯藏后期及銷售階段逐漸發病,發病部位表皮凹陷變軟,果皮上形成圓形或橢圓形的病斑,病斑邊緣呈現水漬狀,嚴重時會導致整個果實腐爛,并具傳染性導致整庫腐爛[1],造成極大的經濟損失。 由于軟腐病發病前期變化不明顯,導致人工篩選效率低下。 因此,開發一種快速軟腐病的早期無損檢測的方法成為果農、 銷售者、 消費者等共同關注的一個問題。

近年來,基于高光譜成像技術的無損檢測研究已經成為熱點。 高光譜成像技術融合了光譜技術和機器視覺技術。 機器視覺技術能夠通過紅、 綠、 藍三個波段觀測病害部位的大小、 紋理等來確定病害類型,但它無法探知病害的生化參數,而光譜技術可以獲取樣本的生物和化學參數,但它不能反映樣本的全局信息。 使用高光譜成像可以一次獲得光譜信息和空間信息,快速完成無損檢測工作[2],已廣泛應用于水果[3-5]、 蔬菜[6]、 肉類[7]等農副產品的相關無損檢測。 目前,桃[8]、 南果梨[9]、 草莓[10]、 庫爾勒梨[11]、 藍莓[12]、 柑橘[13]、 蘋果[14]等常見水果的高光譜成像病蟲害檢測均能達到快速、 無損的早期分類檢測效果。 這些研究表明,高光譜成像技術在針對水果病蟲害的早期無損檢測具有很大的應用潛力。

對于獼猴桃樣本,高光譜成像技術目前在可溶性固形物含量檢測[15]、 冷害分級[16]、 酸度檢測[17]、 成熟度檢測[18]等方面均有相關報道。 這些研究表明,高光譜成像技術能夠反映獼猴桃的內在成分,比如內在品質參數、 成熟度、 損傷等方面,而利用高光譜成像技術對獼猴桃軟腐病進行早期檢測鮮見報道。 結合前人相關研究,探索和驗證使用高光譜成像技術對獼猴桃軟腐病進行早期檢測的可行性。

上述使用高光譜成像對水果進行無損檢測的大部分文獻僅使用了光譜特征作為分類依據,而沒有利用高光譜圖像“空譜合一”的優點。 對于高光譜圖像特征的提取,主要是對其進行紋理特征的提取,Sun等[8]基于主成分分析(principal component analysis,PCA)的PC1圖像進行了紋理特征的提取,Yu等[9]基于PC1和PC2圖像進行紋理特征的提取,兩者研究結果均表明通過結合高光譜圖像的光譜特征與圖像特征能夠進一步提高分類的準確率。 然而,使用主成分圖像進行圖像特征提取時,容易存在較大的光譜失真,不能充分利用高光譜圖像多波段圖像的信息。 本工作以“云海一號”獼猴桃為研究對象,首先將獼猴桃的特征波段圖像融合得到融合圖像,并提取融合圖像的紋理特征。 同時,提取高光譜圖像的光譜特征,然后將光譜特征與紋理特征融合后使用支持向量機(support vector machine,SVM)進行分類,對獼猴桃感染軟腐病后進行了早期檢測與識別。

1 實驗部分

1.1 材料

實驗樣本采摘自湖北省武漢市初陽奇異果園,選擇成熟度處于商業采收期、 可溶性固形物為8%～9%、 果重在80～110 g的“云海一號”獼猴桃作為樣本,并剔除了畸形果和疤痕果,然后使用咪鮮胺和抑酶唑進行處理,等待晾干后放入冷庫中保存。 用于樣品軟腐病接種的是擬莖點霉菌(Phomopsissp.),由華中農業大學園藝林學學院相關研究專業人員對獼猴桃進行接種。

由于獼猴桃的球形形狀會導致高光譜圖像中光譜反射的不均勻性,為了最大限度的減少這種影響,接種區域選擇在獼猴桃的中央赤道部位。 首先在無菌條件下使用滅菌后的移液槍進行接種。 實驗分為健康對照組和實驗組,健康對照組注射無菌生理鹽水來消除接種處傷口對分類識別的影響,實驗組使用兩種不同劑量的接種方法來模擬不同情況下的獼猴桃軟腐病發病情況,其中第一組接種6 μL孢子懸浮液,第二組接種10 μL孢子懸浮液。 接種完成后,對所有獼猴桃進行單果套袋并將接種面朝上,室溫20 ℃貯藏。 將接種當天記為第0天,并從第1天開始每天下午1點進行高光譜圖像的采集,同時關注軟腐病的發展狀況,直到獼猴桃軟腐病發病嚴重或出現白色菌絲后停止高光譜圖像的采集。

1.2 高光譜圖像采集系統

采用比利時IMEC公司開發的SNAPSCAN VNIR高光譜成像儀采集高光譜圖像,采集系統如圖1所示,該系統主要由相機機身、 鏡頭、 升降平臺、 鹵素燈光源及掃描控制器組成。 該系統配有官方的圖像采集軟件HSISnapscan,HSI(hyperspectral image); 傳感器有150多個可用的波段,經系統處理后,在470～900 nm之間有150個可用光譜波段,利用軟件可自動黑白校正得到校正后的光譜反射率數據。 根據獼猴桃大小選擇的成像分辨率為800像素×800像素,相機積分時間為4 ms,模擬增益為1.6 dB。 圖像采集時,整個系統放置于一個密閉的黑箱中,防止外部光源的干擾。

圖1 高光譜圖像采集系統(a): 示意圖; (b): 實物圖Fig.1 Hyperspectral image acquisition system(a): Schematic diagram; (b): Actual picture

1.3 方案設計

實驗方案由三部分組成: 圖像采集及處理、 特征提取及融合、 模型建立及評估。 總體實驗流程如圖2所示。

圖2 總體實驗流程Fig.2 Overall experimental process

所使用軟件除1.2節中提到的HSI Snapscan外,圖像感興趣區域的選擇使用了Labelme軟件,除此之外其他數據處理及建模分析均由Python3實現。

1.4 高光譜圖像處理及數據處理

1.4.1 獼猴桃圖像信息及光譜特征分析

圖3是不同類型的獼猴桃樣品,其中圖3(a)代表健康獼猴桃,圖3(b)、 (c)和(d)分別代表接種擬莖點霉菌(6 μL)后第四天、 第五天和第七天的獼猴桃。 由圖3可以看出,接種擬莖點霉菌后第四天的獼猴桃與健康獼猴桃并無明顯區別,而在第五天,獼猴桃接種部位變軟,說明從第五天起肉眼觀察可以識別出軟腐病發病。 在接種后的第七天獼猴桃軟腐病發展嚴重,有明顯的水漬狀。

圖3 獼猴桃樣品(a): 健康; (b): 接種第四天; (c): 接種第五天; (d): 接種第七天Fig.3 Kiwifruit samples(a): Healthy; (b): The fourth day of vaccination; (c): The fifth day of vaccination; (d): The seventh day of vaccination

根據獼猴桃軟腐病的發展過程,將所采集的獼猴桃樣本分為三類: 第一類為健康對照組,第二類為軟腐病早期,第三類為軟腐病晚期。 軟腐病發病早期是根據病斑肉眼看不出明顯擴大,其大小為接種口的直徑為0～10 mm,晚期是指病斑肉眼可見的擴大,其直徑>10 mm。 實驗過程中發現,第二實驗組比第一組發病提早一天,因此,用于后續分類研究的樣本為: 第四天至第七天樣本(實驗組一),第三天至第七天樣本(實驗組二)。 最終用于實驗的獼猴桃樣本有295個,其中第一類樣本60個,第二類樣本128個,第三類樣本107個。

將所有樣本分類標定后,每個樣本以接種部位為中心,選擇對應的感興趣區域(region of interest,ROI)用于后續的圖像特征提取,其中ROI區域的平均反射率光譜代表樣本的原始光譜特征。

圖4顯示了不同類別獼猴桃的平均反射率光譜曲線,可以看出不同病害程度的獼猴桃平均光譜曲線總體趨勢上是相似的,在630 nm附近與670 nm均存在波峰和波谷,這主要是由于獼猴桃果肉和表皮中的葉綠素吸收引起的。 軟腐病早期和晚期的獼猴桃光譜曲線比較接近,但健康獼猴桃的光譜反射率比軟腐病獼猴桃低,這可能是真菌繁殖過程導致獼猴桃果實中的成分發生了變化。 這說明光譜曲線具備將健康果與軟腐病果初步區別的潛力,但是無法將發病早期和晚期果準確區分,需對光譜特征處理進一步分析。

圖4 不同類別獼猴桃光譜信息Fig.4 Spectra of different categories of kiwifruit

1.4.2 光譜特征提取

為了減少光譜數據中噪聲的影響,通常需要對原始光譜數據進行預處理,常用的光譜預處理方法有: SG卷積平滑(savitzky-golay)、 多元散射校正(multiplicative scatter correction,MSC)和變量標準化(standard normalize variate,SNV)。 MSC算法通過擬合樣本與理想光譜的關系來對光譜進行校正,而SNV算法則是通過對每條光譜做標準正態變換來對光譜進行校正。 兩者均是在已知樣本所屬類別的基礎上對該類別的平均光譜曲線進行校正,然而在之后的模型測試階段,樣本類別是未知的,故不采用這兩種預處理方法。 SG卷積平滑預處理可以減少隨機噪聲對光譜的影響,而本研究所采用的高光譜儀器在導出光譜數據時已進行過校正,有效減少了隨機噪聲的影響。 對于輸入特征只包括光譜特征時,分類器采用SVM,表1展示了使用SG卷積平滑與不使用預處理方法的實驗結果對比,可以看出,使用SG卷積平滑后,最好的準確率與不使用預處理相同,但在SPA或CARS光譜特征提取下,使用SG得到的結果要更差一點,這說明,對于本文實驗中采集的高光譜數據而言,已不需要再進行數據的平滑,使用SG平滑反而會造成過度平滑惡化結果。 綜合考慮,在獲取原始光譜特征后,不進行光譜預處理直接進行光譜特征提取。

表1 光譜數據未處理和經預處理的結果對比Table 1 Results of input data with different preprocessing methods

采集的高光譜圖像包含150個波段,使用全波段進行建模時,分析時間長,且各波段間相關性較強,很容易造成信息冗余,因此需要對原始光譜進行特征降維或者特征選擇。 采用PCA進行特征降維,選擇連續投影算法(successive projections algorithm,SPA)和競爭性自適應重加權算法(competitive adaptive reweighted sampling,CARS)進行特征波段的選取。

PCA降維通過對原始波段特征進行線性加權組合形成新的互不相關的主成分變量,根據主成分貢獻率的大小對主成分數量進行選擇,可以有效降低特征維度。 對原始光譜特征進行PCA降維后,PC1、 PC2、 PC3貢獻率分別為81.98%、 14.7%、 2.1%。 前三個主成分在三維空間中的分布情況如圖5所示,從圖中可以看出,三種類別的樣本之間有所交叉,說明僅根據前三個主成分還無法完全對樣本進行分類,后續還需借助分類器進行建模分析。 為了保證后續分類器的分類性能,選擇保留累計貢獻率達到99.9%的前10個主成分作為原始光譜數據降維后的光譜進行下一步分析。

圖5 前三個主成分分布Fig.5 The first three principal components distributions

SPA算法根據最大投影向量不斷循環,在篩選特征波段的同時降低了變量組之間的共線性關系。 設置波段選擇過程中的最少波段數為2,最多波段數為50,發現在選擇波段數為3時,均方根誤差波動幅度減小,均方根誤差(RMSE)為0.49。

CARS算法通過將蒙特卡洛采樣與PLS回歸系數相結合,在每次自適應加權采樣過程中,保留回歸系數絕對值大的波段,剔除絕對值較小的波段,經過多次采樣選擇交叉驗證均方根誤差(RMSECV)最小的波段子集。 設置蒙特卡洛采樣次數為50,5組交叉驗證,發現在采樣次數為15時,RMSECV最低,此時篩選的特征波段數時39個。 表2展示了SPA和CARS算法所選擇的具體波段,SPA選取的波長處于可見光區,而CARS選取的波長在可見光和近紅外區均有分布。

表2 特征波段選擇結果Table 2 Feature band selection results

1.4.3 圖像特征提取

高光譜圖像各波段的圖像可以清晰的反映軟腐病獼猴桃的圖像特征,因此,在提取光譜特征的基礎上加入波段圖像特征,可以進一步強化獼猴桃的病害信息,增強健康組與不同程度病害組獼猴桃之間的區分度。 由于采集的高光譜圖像波段數多達150個,如果對每個波段進行圖像特征提取會造成極大的計算量,并且其中存在一些對分類貢獻較小的波段,因此,在提取圖像特征之前,需要進行特征波段的選擇。

選擇的特征波段是在SPA過程中產生的一個波段范圍子集,共有8個波段,分別是第0、 11、 27、 62、 77、 87、 132、 149波段,對應波長為470、 502、 548、 649、 692、 721、 851和900 nm。 沒有使用SPA最終選擇的3個波段的原因是: 對于光譜特征而言,這3個波段的平均反射率在光譜信息上具有很好的代表性,但對于圖像特征而言,不同波長圖像下樣本對光進行吸收反射后對圖像特征的反映也不同,所篩選的8個波段能夠很好的覆蓋原始150個波段中不同的波段范圍。

采用灰度共生矩陣(gray-level co-occurrence matrix,GLCM)進行圖像紋理特征的描述。 設置不同的距離(2、 4、 8、 16)以及不同的角度(0°、 45°、 90°、 135°)進行灰度共生矩陣的提取,并計算灰度共生矩陣的對比度(contrast)、 同質性(homogeneity)、 相關性(correlation)、 能量(ASM)、 熵(entropy)。 將4個距離以及4個角度的平均值作為圖像的紋理信息,將每幅圖片計算得到的5個紋理參數代表該圖像的紋理特征。

對于各子帶系數的融合規則,通常對低頻子帶進行加權平均,對高頻子帶進行絕對值取最大的方法來求解。 經過實驗發現,在使用絕對值最大法對高頻子帶進行融合時,由于每個像素值只用到了一個波段圖像信息,最后所得到的高頻子帶系數將存在一定的信息丟失。 因此,為了在簡化融合規則的同時考慮融合更多的細節信息,對于高頻子帶系數也采取了加權平均的方法。 式(1)和式(2)分別描述了低頻子帶和高頻子帶的融合規則。

(1)

(2)

圖6 高頻子帶圖像(a): 第132波段; (b): 絕對值最大融合; (c): 加權平均融合Fig.6 High frequency subband image(a): Band 132; (b): Absolute maximum fusion; (c): Weighted average fusion

1.4.4 模型建立與評估

采用Kennard-Stone算法將采集的實驗樣本按照訓練集∶測試集=7∶3進行劃分,共獲得訓練集207個,測試集88個。 由于樣本的光譜特征與圖像特征量綱不同,在將兩種特征融合之前需要進行歸一化處理。 采用最大最小值歸一化方法將兩組特征歸一化至[0,1]之間,計算過程如式(3)

(3)

式(3)中,X表示原始特征,Xmin、Xmax分別代表原始特征中的最小值和最大值,Xnorm為歸一化至[0,1]之間的特征,作為后續分類器的輸入。

選用三種分類模型作為驗證,分別是最近鄰算法(K-nearest neighbor,KNN)、 隨機森林(random forest,RF)和支持向量機(SVM)。 模型在訓練過程中采用了十折交叉驗證,并通過測試集的分類準確率(accuracy)和混淆矩陣對模型進行評價,分類準確率的計算過程如式(4)

(4)

式(4)中,TP、TN分別表示模型預測正確的正、 負樣本,FP表示模型預測為正類的負樣本,FN表示模型預測為負類的正樣本。

2 結果與討論

2.1 識別準確率分析

基于不同處理組合對獼猴桃樣本數據進行分析,實驗結果如表3所示,其中前六種方式利用的是獼猴桃樣本的高光譜圖像數據或RGB圖像數據,分為僅使用光譜特征、 僅使用紋理特征和使用光譜、 紋理融合特征。 最后一種是將RGB光譜特征結合HSI波段融合圖像的紋理特征進行分析,即融合了RGB光譜數據和高光譜圖像數據。

表3 不同組合下測試集分類結果Table 3 Classification results of test sets under different combinations

從表3可以看出,(1)KNN是三種分類模型中分類效果最差的,RF與SVM兩種模型的測試集分類準確率較為接近,但使用RF分類器時,訓練集和驗證集的分類準確率大部分為100%,而測試集準確率基本處于80%～88%之間,存在較嚴重的過擬合現象。 使用SVM分類器時,訓練集、 驗證集和測試集的分類準確率基本相近,有不錯的泛化能力,因此,SVM是最適合識別獼猴桃軟腐病的分類器。 分類器模型的參數選擇是使用十折交叉驗證結合網格尋優進行確定的,表3在SVM-acc列的括號處給出了本實驗中SVM的具體參數,核函數均使用RBF。 (2)單獨使用光譜特征和紋理特征的分類結果都不夠理想,而將兩種特征進行融合可以獲得更好的分類效果,其中將波段降維后的光譜特征和紋理特征進行融合的分類結果要優于全波段光譜與紋理特征的融合,證明了光譜降維的有效性。 (3)和高光譜圖像相比,RGB圖像的分類準確率偏低,證明了高光譜圖像在獼猴桃軟腐病分類中的有效性。 此外,RGB圖像光譜值與本文提出的波段融合圖像的紋理特征融合后也能得到85.23%的準確率,說明本文通過提取特征波段圖像融合后的紋理信息能夠有效強化獼猴桃軟腐病的病害特征,從而提高了模型的預測能力。

使用光譜特征、 RGB圖像光譜特征以及融合特征進行SVM分類的混淆矩陣結果如圖7所示。 可以看出,單獨使用光譜特征對健康獼猴桃和軟腐病早期獼猴桃進行分類的準確率不高,并且不能區分發病早期和發病晚期的樣本,這與前文的光譜分析結果一致,發病早期和發病晚期樣本光譜曲線比較接近,僅使用光譜特征對兩類獼猴桃難以區分。 而使用RGB圖像的光譜特征進行分類時,不同類別樣本的分類結果均不理想,進一步說明RGB圖像不適合分類研究。 當使用融合特征進行分類時,各個類別的誤判率最低,大多數的獼猴桃軟腐病發病早期樣本可以得到正確識別,說明本方法能夠實現獼猴桃軟腐病的早期無損檢測,同時在獼猴桃感染軟腐病第三天或第四天病斑還未顯現之時,通過提取光譜特征與融合圖像的紋理特征的方法能夠將軟腐病早期獼猴桃與健康獼猴桃分離出來,有效減少了成本損失。

圖7 模型測試結果(a): 光譜特征; (b): RGB光譜特征; (c): 融合特征Fig.7 Model test results(a): Spectral features; (b): RGB spectral features; (c): Fused features

本研究還對波段圖像融合及紋理特征提取過程中的其他情況進行了對比實驗,分別是(1)高頻子帶系數融合采用絕對值取最大方法,記作max-fusion。 (2)紋理特征提取時,分別計算不同角度的特征,每幅圖片提取20個紋理參數代表該圖像的紋理特征,記作four-angle。 (3)對不進行融合的8個特征波段提取紋理特征得到40維的向量作為圖像特征,記作non-fusion。 (4)使用SPA選擇的三個波段進行融合,記作SPAband-fusion。 實驗結果如圖8所示。 max-fusion實驗準確率均偏低,與1.4.3節分析一致,采用絕對值最大法融合后提取的紋理特征,高頻子帶存在一些圖像失真以及信息丟失,該紋理特征無論是和哪種光譜特征進行融合,結果均劣于加權平均法融合后所提取的紋理特征。 在four-angle與non-fusion兩個實驗中,紋理特征與SPA光譜特征融合的分類結果也偏低,這是因為SPA所得光譜特征僅有3維,而兩種方法所得紋理特征分別為20、 40維,紋理特征存在冗余,特征融合后削弱了光譜特征在分類中的作用。 SPAband-fusion實驗則表明,本工作所選的8個特征波段包含了更多的圖像特征,因此結果要優于SPA三個波段的融合。

圖8 對比實驗結果Fig.8 Comparative experimental results

2.2 與其他方法對比

為了更全面的說明本工作在波段融合后提取紋理特征的有效性,本節與一些相關文獻中的圖像特征提取方法進行了對比,使用的數據集均為本工作的獼猴桃軟腐病數據集,統一使用SVM分類器進行分類。 結果如表4所示,可以看出本方法優于其他方法。 其中文獻[8]分別提取了6維顏色特征,16維GLCM特征以及32維直方圖特征; 文獻[9]使用光譜特征與主成分圖像提取的紋理特征融合后進行分類; 文獻[19]根據PC1權重曲線波峰波谷獲取特征波段,并將特征波段圖像的平均灰度值作為圖像特征; 文獻[20]在SPA獲得的波段中選擇了一個最為清晰、 亮度均勻的波段圖像,利用此圖像提取了共16維GLCM特征。

表4 本方法與其他方法的對比Table 4 Comparison between this method and other methods

在進行圖像特征提取時,類比的眾多文獻大多選擇主成分圖像或特征波段圖像進行提取,本工作與文獻[9]雖然使用的特征變量都是8維,但主成分圖像由于造成了光譜失真導致主成分圖像提取的紋理特征對分類的效果要劣于圖像融合后提取的紋理特征。 文獻[20]雖然使用的特征變量是19維,但是其中有16維是圖像特征,削弱了光譜特征的分類能力,此外,僅使用一個波段圖像提取的紋理特征也不夠充分,故分類準確率較低。 而其他文獻中輸入分類器的特征維度都較高,特征中可能存在著冗余特征,因此分類結果較差。

3 結 論

以湖北省武漢市“云海一號”獼猴桃為研究對象,通過對健康獼猴桃及感染軟腐病的不同時期獼猴桃進行高光譜圖像采集,提出了一種特征波段圖像融合的獼猴桃軟腐病早期分類檢測方法。 首先,對樣本原始光譜分別使用PCA,SPA,CARS進行光譜特征提取得到光譜特征,然后對SPA求解過程中所得的8個特征波段圖像使用NSCT進行圖像融合并提取融合圖像的紋理特征,最后通過融合SPA光譜特征與紋理特征建立SVM模型對軟腐病早期分類效果達到最優,主要得出以下結論:

(1)健康獼猴桃和軟腐病獼猴桃光譜差異明顯,可以利用光譜特征將兩者區分,但是獼猴桃在軟腐病早期和晚期的光譜差異較小,僅利用光譜信息難以區分。

(2)相比于單一光譜特征和圖像特征,將光譜特征與圖像特征進行融合可以顯著提高分類準確率,說明圖像的紋理特征可以在一定程度上彌補獼猴桃軟腐病在光譜特征上的分類缺陷。 通過圖像融合后再提取紋理特征,可以充分利用高光譜圖像“空譜合一”的優點,實現在降低紋理特征數據維度的同時,進一步提高分類準確率。

(3)本研究使用高光譜成像技術能夠在獼猴桃感染軟腐病3～4 d時將染病果與健康果成功區分,實現了獼猴桃軟腐病的早期無損檢測,為獼猴桃的銷售分級提供了一定的指導意義。