基于ICO-SPA特征提取的近紅外高光譜小麥赤霉病粒嘔吐毒素含量預測

丁 靜，梁 琨，韓東燊，徐劍宏，沈明霞

(1.南京農業大學工學院/江蘇省現代設施農業技術與裝備工程實驗室，江蘇南京 210031；2.江蘇省農業科學研究院農產品質量安全與營養研究所，江蘇南京 210014)

小麥在世界范圍內種植面積廣泛，既是人類主要的食物來源，又是重要的工業原料[1]。赤霉病是小麥的主要病害之一，在小麥抽穗揚花期氣候的溫、濕度適宜且雨水充沛時容易高發。赤霉病的病原菌主要是禾谷鐮刀菌，該菌在小麥籽粒中所產生的代謝產物脫氧雪腐鐮刀菌烯醇(DON)嚴重危害人畜健康，人畜誤食帶菌病粒會引起腹瀉、嘔吐等癥狀，因此DON又被稱作嘔吐毒素[2-5]。我國GB 2715-2005《糧食衛生標準》中規定小麥及全麥粉中DON限量標準為1 mg·kg-1。傳統的檢測方法主要是生物化學檢測法，包括皮膚毒性實驗、氣相色譜(GC)、液相色譜(LC)、氣(液)質聯用(GC-MS)、酶聯免疫吸附試驗(ELISA)等，這些檢測方法耗時、昂貴、重現性差、樣本處理過程復雜、需要專門的技術人員操作，難以進行大批量樣本的快速無損檢測。

高光譜成像技術融合了光譜技術與成像技術，不僅可根據提取的樣本光譜特征表征和識別給定的物質，還可通過對給定的物質進行成像，從而達到擴大檢測面積、提高檢測精確度的目的[6]，為進行復雜綜合的小麥赤霉病癥狀檢測提供保障。另外，由于樣本前處理簡單、高光譜圖像采集方便及獲得的樣本信息豐富，已有研究者將高光譜技術應用于小麥赤霉病檢測。Shahin等[7]將加拿大西部紅春小麥籽粒分為健康、鐮刀菌輕度感染和嚴重感染三個等級，并使用主成分分析(principal component analysis，PCA)在450～950 nm范圍內提取6個特征波段(484、567、684、817、900、950 nm)用于線性判別分析(linear discriminant analysis，LDA)建模，結果表明，基于6個特征波段建模與基于全波段建模準確率相當，均高于92%；Delwiche等[8]通過在400～1 000 nm和1 000～1 700 nm兩個波段范圍尋找能夠區分正常籽粒與赤霉病粒的特征波段對(分別為502、678 nm和1 198、1 496nm)建立二等級分類模型，可見光和近紅外波段分級準確率分別為94%和97%。梁 琨等[9]比較了連續投影算法(successive projections algorithm，SPA)和競爭性自適應重加權算法(competitive adaptive reweighted sampling，CARS)在小麥赤霉病籽粒識別上的應用效果，發現SPA選取的特征變量數量少，精度高，模型預測效果優于CARS，結合支持向量機(support vector machines，SVM)模型能夠達到94%以上的預測準確率。這些研究結果表明，高光譜技術結合特征波長挑選算法能夠通過少數幾個特征變量進行赤霉病籽粒識別。以上研究均是對小麥感染赤霉病的定性分析，使用高光譜成像技術進行小麥中嘔吐毒素定量分析的報道較少。Barbedo等[10-11]利用基于高光譜圖像的處理算法，采用數學形態學操作和光譜波段操作實現對小麥鐮刀菌毒素(FHB)的快速自動檢測，提出使用鐮刀菌指數(FI)計算小麥籽粒中FHB存在的可能性，進而估算DON濃度，準確率約為91%，但是FI值在0.4～0.6范圍內的正態分布曲線存在重疊區域，在這個區域FI值存在不確 定性。

感染赤霉病的小麥產生的DON毒素會導致籽粒中的主要大分子有機物成分如淀粉、蛋白質、脂肪、纖維素的含量發生變化[12-13]，這種變化能夠反映在光譜上[11,14-15]，從而具備了使用光譜技術檢測DON的可能性。本研究基于近紅外高光譜成像技術進行DON定量檢測，將1 000～2 500 nm范圍內的光譜數據與標準方法測得的DON含量值進行關聯，將SPA以及區間組合優化算法(interval combination optimization algorithm，ICO)結合SPA提取的特征變量作為輸入建立偏最小二乘回歸(partial least squares regression，PLSR)、多元線性回歸(multiple linear regression，MLR)和最小二乘支持向量機回歸(least squares support vector regression，LS-SVR)模型預測樣本中DON含量，比較不同特征波長篩選算法并結合預測模型精度，獲得小麥赤霉病檢測的優選波長和優化模型，以期實現小麥感染赤霉病等級的高效、智能定量檢測，保證糧食食品 安全。

1 材料與方法

1.1 材料與儀器

試驗中用于磨粉的小麥籽粒均來自江蘇省農業科學研究院農產品質量安全與營養研究所，供試材料是2017年收獲的冬小麥，品種為煙農19。所用儀器包括PM8188A谷物水分測定儀(泰州市維科特儀器儀表有限公司)、XA-1型磨粉機(常州越新儀器制造有限公司)、Image-λ-N25E-HS型高光譜成像光譜儀(芬蘭 Specim公司)、HSIA-OLE23型鏡頭(德國施耐德公司)、4個HSIA-LS-T-H型鹵素燈(四川雙利合譜科技有限公司)、HSIA-T1000-IMS型電控位移平臺(四川雙利合譜科技有限公司)和GaiaSorter-Dual型暗箱(四川雙利合譜科技有限公司)。

1.2 試驗方法

1.2.1 樣本前處理

對獲取的小麥赤霉病籽粒樣本人工分揀除去麥稈、石子、土塊、塑料等異物，使用谷物水分測定儀檢測籽粒含水量，對于含水量過高的樣本進行通風干燥，含水量較低的樣本暴露在潮濕空氣中，最終將所有小麥赤霉病籽粒樣本含水量控制在12%～13%范圍[16-19]，以減少樣本含水量差異對光譜的影響，保證建模數據的可靠性。取約25 g小麥籽粒放入磨粉機中研磨30 s，將磨粉機杯蓋內壁附著的未磨碎的部分掃入磨粉倉，蓋上杯蓋繼續磨粉30 s，保證全部顆粒均可通過20目篩，以便后續使用標準方法檢測樣本中DON含量，取出研磨后的全麥粉，混勻，置于直徑60 mm、高10 mm的培養皿中，將稱量紙附在培養皿中的全麥粉表面起隔離作用，再用塑料平板將培養皿中的全麥粉壓平，使全麥粉表面高度與培養皿邊緣持平，用于光譜數據采集。

1.2.2 圖像獲取與標定

高光譜成像光譜儀的成像波段是1 000～ 2 500 nm，光譜分辨率為12 nm，采樣率為5.6 nm，像元數為384×288。在高光譜圖像數據獲取前，打開高光譜成像系統的光源預熱30 min。為了得到清晰無畸變的圖像，多次調整試驗系統參數后將曝光時間設置為10 ms，傳送帶移動速度為0.36 deg·s-1，樣本距鏡頭310 mm。采集光譜前，需要對儀器進行黑白校正，以減弱光譜儀中暗電流對圖像的影響。首先采集反射率為99%的白板漫反射圖像，記為W，然后蓋上鏡頭蓋采集暗圖像，記為B，根據式(1)計算原始漫反射圖像R的校正圖像Rt：

Rt=(R-B)/(W-B)×100%

(1)

1.2.3 ROI的選取與高光譜數據提取

將采集到的高光譜圖像用ENVI 4.8(美國 Exelis Visual Information Solutions公司)進行感興趣區域(region of interest，ROI)提取。沿圖像中每個培養皿的輪廓手動選取輪廓以內的全麥粉作為每個樣本的感興趣區域。將每個樣本ROI區域的平均反射率值作為該樣本的原始光譜數據，后續數據分析在MATLAB 2014a(美國 The MathWorks公司)軟件中完成。

1.2.4 DON毒素含量的測定

采集完高光譜圖像后，參考SN/T3137-2012標準中規定的液相色譜-質譜法定量檢測每個樣本中的DON含量，作為每個樣本的理化參考值。提取液由乙腈和水按 84∶16的比例配置，離心時轉速為2 500 r·min-1，活化時使用3 mL提取液，提取的樣本過柱時的流速為1 mL·min-1。利用型號為3500QTRAP色譜儀-液相色譜質譜聯用儀(ABSCIEX公司)進行毒素含量測定，流動相A為5 mmol·L-1的醋酸銨水，按照時間0、3、7、13、13.1、16 min時，流動相A和流動相B分別按照85∶15、30∶70、20∶80、10∶90、85∶15、 85∶15的濃度梯度進行操作，流速為0.6 mL·min-1，進樣量為5 μL。毒素檢出限為20 μg·kg-1。

1.3 特征變量提取

SPA采用連續投影策略進行變量排序產生一系列變量子集，然后通過比較變量子集模型的預測能力篩選出最優變量子集，能最大限度降低被選變量間的共線性[20]。但對于近紅外高光譜而言，有效變量間的投影距離并不一定最大，SPA篩選出的變量子集中可能包含一些無信息甚至是干擾變量[21]。ICO是基于模型集群分析策略(MPA)開發出的變量選擇方法，MPA強調從數據集中挖掘有效信息時，應分析大量隨機產生的子集模型信息而非單一模型信息[21]。高光譜包含豐富的樣本數據，同時也不可避免地包含較多噪聲，ICO作為基于柔性收縮策略進行波長區間組合優化的算法，采用波長區間替換波長點作為優化對象，一些與待測指標可能存在偶然相關性的噪聲波長點會因為所在光譜區間的重要性不足而被排除，從而能夠減少偶然誤差，更穩妥地保留反映真實樣本信息的有效波長區間。區間挑選算法選擇的波長數量一般較多，不便于快速預測，最終選中的變量區間可以采用單一變量選擇算法進行精簡[22]。將ICO選擇的優化組合區間結合SPA進行特征波長提取，既能有效地鎖定特征波段區間，又能夠減少有效區間內特征變量的冗余，形成優勢互補。

1.3.1 連續投影算法(SPA)

設光譜矩陣為XM×J(M為樣本數，J為波長數)，記xk(0)為初始波長迭代向量，N為需要提取的波長變量個數，該算法從一個波長開始，計算它在未選入波長上的投影，將投影向量最大的波長引入到波長組合，直到循環N次結束計算。具體步驟[23]如下：

(1)第1次迭代之前(n=1)，在訓練集光譜矩陣X中任選一個列向量xj，記為xk(0)；

(2)令S為所有未被選入的波長變量的集合，S={j,1≤j≤J,j?{k(0),k(1),…k(n-1)}}；

(3)利用式(2)計算xj對S中向量的投影：

(2)

(4)根據式(3)獲得投影值范數最大的波長 位置：

k(n)=arg(max|Pxj|,j∈S)

(3)

(5)根據將最大投影向量作為下輪的投影向量：xj=Pxj,j∈S；

(6)令n=n+1，若n

1.3.2 區間組合優化算法(ICO)

ICO算法首先采用軟收縮的方式對波長區間組合進行優化，然后采用局部捜索的方式對最終入選區間的寬度進行自動優化，運算過程[24]如下：

(1)將樣本的近紅外光譜劃分為寬度大體相當的N等份，每一份即為一個光譜區間；

(2)采用加權自舉采樣(WBS)生成M個不同波長區間隨機組合形成的子集，每個波長點初始采樣權重都默認為1；

(3)采用PLS算法和5折交互檢驗的方式計算每個區間組合子集對應的RMSECV值；

(4)從全部區間組合中提取一定比例(記作α)的最優區間組合子集(α通常可設置為0.1，也可根據數據的實際情況進行優化)，并計算出這一部分區間組合子集對應RMSECV的平均值；

(5)統計每個區間在那一部分最優區間組合中出現的數目，下一次迭代中第i個區間對應的采樣權重可以通過式(4)計算出來：

(4)

其中fi是第i個區間在最優區間組合中出現的頻次，kbest是提取的最優區間組合的數目；

(6)重復步驟2～5形成迭代循環，直到RMSECV平均值出現上升，迭代終止；

(7)在最后一次迭代中，RMSECV值最小的那一組波長區間被視為最終選中的波長區間；

(8)采用局部捜索策略對步驟7中選出各個波長區間的寬度進行優化。在此策略中，各個被選中波長區間邊緣相鄰的波長點逐一被納入或排除建模，并根據該波長點對模型RMSECV值的影響情況進行評價。如果一個相鄰波長點的納入會使模型的RMSECV值下降，則該波長點被選入，反之則被剔除。反復在每次優化后區間上進行局部捜索，直到沒有新波長點能夠對模型的RMSECV值產生影響，這個時候優化出的區間則可以被視為ICO算法最終選中的波長區間。

1.4 建模方法與模型評價

1.4.1 異常樣本的剔除與光譜數據集劃分

異常樣本的存在影響定標模型的預測性能和適配性。為避免異常樣本對模型精度的干擾，通過基于XY變量聯合的異常數據檢測算法(outlier sample eliminating algorithm based on joint XY distances，ODXY)[25]剔除異常樣本3個。采用基于光譜-理化值共生距離算法(sample set partitioning based on joint X-Y distances，SPXY)[26-27]，從樣本集中篩選校正集樣本。對117個樣本按照SPXY算法以3∶1比例劃分校正集和預測集，得到校正集樣本89個，預測集樣本28個，并且DON含量最大值和最小值對應的樣本均落在校正集中。

1.4.2 光譜預處理

為減少電噪音、光散射、基線漂移等干擾對光譜數據的影響，本研究采用Savitzky-Golay(SG)平滑、標準正態變量(standard normal variable，SNV)、矢量歸一化(normalize)、移動平均算法(moving average，MA)、SG-SNV、SG-Normalize、SG-Normalize-MA等7種算法分別對原始光譜進行預處理，再根據PLSR模型的效果尋找最合適的光譜預處理方法。

1.4.3 模型建立與評價

分別將SPA和ICO-SPA提取的特征變量作為輸入建立PLSR、MLR、LS-SVR模型預測樣本中DON含量，比較模型性能選取最優模型。模型性能利用以下指標進行評價：校正集相關系數(rc)、校正集均方根誤差(root mean square error of calibration，RMSEc)、預測集相關系數(rp)、預測集均方根誤差(root mean square error of prediction，RMSEP)、交叉驗證均方根誤差(RMSECV)和預測相對分析誤差(residual predictive deviation，RPD)。相關系數越大，均方根誤差越小，所建模型性能越好。當RPD值為1.5～2.0時僅能夠區分響應變量的低值和高值；當RPD值為2.0～2.5時能夠粗略地定量預測；當RPD值為2.5～3.0或者更高時模型預測精度較高[28]。

2 結果與分析

2.1 樣本光譜特征

試驗所采用的光譜儀采集的波長范圍為 920～2 528 nm，由于存在系統誤差，光譜曲線在首尾端可能存在較大噪聲，會直接影響試驗的準確性，因此數據處理時取1 000～2 500 nm波段的268個波長的光譜數據作為后續數據處理區域。在圖1的樣本光譜曲線上，在1 470、1 940 nm處出現明顯吸收峰，反映麥粉樣本中水分對光譜的影響[8]；1 201 nm與脂肪的C-H彎曲振動2倍頻帶1 195 nm接近[29]；1 582 nm與淀粉的特征峰 1 550 nm接近[30]；1 767 nm與纖維素C-H彎曲振動1倍頻帶1 780 nm接近[29]； 2 103 nm與多糖的O-H伸縮和C-O伸縮振動組合頻帶 2 100 nm接近[29]；2 299 nm與蛋白酰胺的C-H伸縮振動二倍頻帶2 300 nm接近[29]。這說明樣本反射光譜曲線反映了試驗樣本中的淀粉、蛋白質、脂肪、纖維素等大分子成分的含量 信息。

圖1 近紅外波段下全部小麥赤霉病籽粒樣本反射光譜曲線

2.2 光譜預處理效果比較

從光譜經過不同預處理后建立的PLSR模型(表1)看，SG-Normalize-MA的PLSR模型最好，建模的相關系數和均方根誤差分別為0.927和0.402 mg·kg-1，模型擬合性較好；預測集的相關系數和均方根誤差分別為 0.926和0.354 mg·kg-1，相對分析誤差為 2.949，說明模型有較高預測精度。因此，最終選取SG-Normalize-MA處理后的光譜數據進行后續處理。

2.3 樣本特征波長提取

1 000～2 500 nm采集的光譜共包含268個波段。為簡化建模過程，加快計算速度，提高模型的預測能力，需要提取光譜的特征波段，壓縮光譜數據，以實現對DON含量的快速定量檢測。本研究分別采用SPA及ICO-SPA算法篩選獲得特征波段，再根據反映的物質信息結合所建模型的預測準確率判斷所提特征波長的合理性。

2.3.1 基于SPA算法的特征波長篩選

設定波長數N的范圍為1～25，根據交叉驗證均方根誤差RMSECV值確定最佳的波長變量個數。從圖2a可以看出，隨著波長個數的增加，RMSECV值逐漸減小；當波長個數大于16時，RMSECV值變化不再顯著，此時RMSECV為 0.544。波長過多時易增加模型的運算量及復雜度，因此本研究選取16個波長作為最終特征波長，依次為1 223、1 285、1 347、1 414、1 526、 1 587、1 660、1 700、1 772、1 812、1 907、1 946、 2 047、2 092、2 137、2 215 nm，所選波長位置如圖2b，這些特征波段能夠反映小麥蛋白質 (1 526、 1 587、2 047)、脂肪(1 223、1 414、2 137 nm)、纖維素(1 772、1 812 nm)和水(1 946、2 092 nm)等大分子有機物的含量信息[29]。

表1 不同預處理方法的PLSR預測模型Table 1 PLSR prediction models for different pre-processing methods

圖2 SPA變量選擇

2.3.2 基于ICO-SPA算法的特征波長篩選

基于ICO算法，利用5折交叉驗證方法建立PLS模型選擇特征變量。PLS模型中最大潛變量數為10，加權自舉采樣次數為1 000，所選子模型的比例為0.01，等間距光譜間隔的數量為40。圖3a為ICO算法迭代過程中每個波長區間的采樣權重隨迭代次數增加的變化的情況。其中，顏色越紅，采樣權重值越接近于1；顏色越藍，采樣權重值越接近于0；如果顏色介于深藍色和深紅色之間，則采樣權重值處于0和1之間。

從圖3a可以看出，第9個波長區間在第1～5次迭代過程中顏色偏綠藍，權重系數約為0.3～0.5，但是在第6次迭代過程中權重系數仍有機會變大，最終被選中。該區間包含淀粉的O-H伸縮振動1倍頻帶(1 490 nm)[29]，是進行DON定量檢測重要依據。第14個波長區間的采樣權重值在第1次迭代過程中已經變為1，其采樣權重值在第2次迭代過程中依然有可能變得小于1，最終該區間由于權重過小被淘汰。由圖3b看出，經過8次迭代后均方根誤差趨于穩定，此時其值 0.520。最終，在ICO挑選的特征波長區間內，共挑選了22個特征波長 (1 083、1 089、1 100、 1 117、1 137、1 145、1 156、 1 319、1 324、1 335、 1 347、1 363、1 375、1 380、 1 475、1 487、1 498、 1 711、1 722、1 733、1 739、 2 445 nm)(圖4)，這些基于ICO-SPA提取出的特征波波段也可反映小麥赤霉病籽粒中淀粉、蛋白質、脂肪、纖維素等有機物大分子含量的特征[29,31]。

圖3 ICO變量選擇

圖4 ICO-SPA變量選擇

Fig.4 Variable selection by ICO-SPA

2.4 建模方法優選

將SPA挑選出的16個敏感波長和ICO-SPA優選出的22個敏感波長作為輸入變量，分別使用PLSR、MLR、LS-SVR 3種方法對DON含量進行建模，結果(表2)表明，基于兩種特征波長的MLR模型的RPD均達到2.5以上，建模效果比LS-SVR(RPD低于2)好，且不存在PLSR模型中出現的欠擬合現象。PLSR模型出現的欠擬合現象可能是由于提取的特征變量過少，后期可以考慮增加新的特征或特征組合用于增大假設空間，減小欠擬合現象。所建模型中，ICO-SPA-MLR模型預測效果最優，預測集相關系數rp為0.921，相對分析誤差為2.789，大于 2.5，能實現精準預測，且預測效果與全波段的PLSR模型建模效果相當(rp=0.926，RPD=2.949)。基于ICO-SPA提取的22個敏感波長建立的MLR模型預測效果與基于全波段268個波長建立的PLSR模型相當，說明ICO-SPA是一種有效的特征波長挑選方法，結合MLR模型能夠接近全波段的預測效果。另外，基于ICO-SPA建立的PLSR、MLR、LS-SVR模型相較于基于SPA建立的PLSR、MLR、LS-SVR模型，其rp分別提高了 2.3%、0.7%、1.6%，RMSEp分別下降了 0.070、0.017、0.018 mg·kg-1，RPD分別提高了 0.506、0.123、0.059，表明ICO-SPA篩選出的變量有效性比SPA高，在ICO選定的有效區間上進行SPA特征波長挑選，最終獲得的特征變量既有效又精簡。ICO-SPA比SPA多提取了淀粉含量的信息，對于控制在同一水平的水分含量信息，ICO-SPA沒有提取，表明淀粉含量信息的加入和水分含量信息的減少有助于DON含量 預測。

從圖5可以看出，DON含量的實測值與預測值的相關性很好，數據集的預測樣本點基本均勻分布在y=x線兩側，說明ICO-SPA-MLR預 測模型能夠實現對樣本中DON含量的定量 檢測。

表2 三種定量模型預測結果Table 2 Prediction results with the three quantitative models

圖5 小麥赤霉病籽粒中DON含量的ICO-SPA-MLR預測模型

3 討 論

本研究是在小麥赤霉病病害發生后，通過近紅外高光譜技術預測小麥赤霉病粒中的DON含量，與國家標準中的相關限量比對，確認該批小麥是否安全，以減小病粒對人畜健康的威脅，保證糧食食品安全；另外，面粉廠等大型小麥原糧收購單位一般會依據原糧的DON含量進行等級劃分，不同等級的原糧存入不同的谷倉，傳統的理化方法由于耗時較長，會嚴重延長原糧入庫的時間，使用近紅外高光譜的方法能夠快速獲得小麥籽粒的DON含量結果，從而指導原糧及時入庫，減少對運糧車或運糧船的長時間占用，提高原糧收購單位的經濟效益。

高光譜技術對光譜范圍內的每一個波段進行成像，在獲得樣品外部特征信息的同時，還提供反映樣品內部物理結構和化學成分的光譜信息，但因此產生的巨大的信息量一方面會影響檢測速度，另一方面由于大量噪聲或冗余信息的存在也會影響檢測精度。為了簡化近紅外高光譜分析模型，提高檢測的實時性和有效性，本研究分別采用連續投影算法(SPA)、區間組合優化結合連續投影算法(ICO-SPA)進行特征波長提取，結合偏最小二乘回歸(PLSR)、多元線性回歸(MLR)和最小二乘支持向量機回歸(LS-SVR)模型，比較兩種特征波長提取方法輸入三種模型后的預測效果，結果表明，ICO算法的引入有助于提高SPA算法所選波長的有效性，基于ICO-SPA挑選出的特征波長作為輸入變量建立的模型可以用來進行小麥赤霉病籽粒樣本的DON含量預測，其中最優模型為基于ICO-SPA提取的特征波段建立的MLR模型，預測集相關系數為0.921，預測集均方根誤差為0.375 mg·kg-1，相對分析誤差為2.789，大于2.5，預測精度略高于前人的研究結果[10-11]。這種將光譜區間挑選算法與單一波長挑選算法聯合使用的策略充分利用了算法之間優勢的互補性，因此所建模型的預測效果比算法單獨使用時有所提升。

在特征波長挑選之前一般要對光譜進行預處理，光譜預處理方法不僅會影響模型的預測性能，還會影響變量選擇算法的結果[22,32]。對SG-Normal-ize-MA預處理后的光譜進行ICO-SPA特征波長提取，所得特征波長能較充分地反映試驗樣本中各有機物成分含量差異，但是光譜預處理方法的組合及運算順序仍需要配合變量選擇算法作進一步優化，這一問題有待進一步探索。本研究針對小麥品種煙農19進行小麥赤霉病粒DON含量預測試驗，ICO-SPA的特征波長挑選算法獲得了讓人滿意的結果，若換成其他小麥品種，ICO-SPA的效果如何有待進一步的研究和驗證，后期還可以進行多個品種混合的小麥赤霉病籽粒DON含量預測，從而提高先進技術和算法對于糧食食品安全以及原糧收購單位經濟效益的應用價值。