基于高光譜技術及SPXY和SPA的玉米毒素檢測模型建立

于慧春，婁 楠，殷 勇*，劉云宏

新鮮玉米在存貯過程中，由于其胚部大、水分含量高、帶菌量多，高溫高濕環境下極易霉變[1]，不僅給經濟造成重大損失，而且霉變玉米在代謝過程中會產生多種對人體具有極強致病性、致癌性的毒素[2-3]，危害人畜健康。黃曲霉毒素B1（aflatoxin B1，AFB1）和玉米赤霉烯酮（zearalenone，ZEN）是其中2 種比較穩定的代謝產物，容易在霉變玉米中積累，導致含量升高，AFB1和ZEN的多少與玉米霉變情況密切相關[4]。因此，可以通過監測玉米中AFB1和ZEN含量變化表征玉米的霉變情況，實現對玉米霉變程度的準確評價[5]。但是常規的AFB1和ZEN檢測方法存在樣品處理繁瑣、費時、對樣品有破壞性等缺點，難以實現簡單、快速、無損檢測，無法滿足實際生產的需要[6-8]。

高光譜技術通過提取被測對象的圖像和光譜信息特征并與其內外品質指標間建立聯系，從而實現對其內外綜合品質的評價，因而在果蔬品質[9-11]、禽肉品質[12-14]、蛋[15]和煙葉[16]等農產品檢測領域獲得廣泛應用研究。在霉變玉米的高光譜檢測方面，也有學者進行了一些探索[17-19]，但目前研究主要集中在高光譜技術對不同霉變程度玉米的定性判別方面，而對霉變玉米中AFB1和ZEN等的定量檢測分析方面尚鮮有報道。袁瑩[20]和褚璇[21]等將不同濃度的AFB1溶液滴在玉米籽粒上，然后利用高光譜技術進行檢測，但其最終的分析結果仍然只是對表面含有不同濃度AFB1的玉米樣品進行分類識別，也未嘗試通過高光譜技術對測試樣本在特定條件下自身代謝生成的AFB1含量進行預測。此外，上述霉變玉米高光譜技術檢測研究中，對所使用的校正集樣本質量也未進行考察，大量的樣本數據之間可能存在相同或者差異過小的情況，從而影響分析結果精度和耗費建模時間。

本研究以玉米霉變過程中產生的AFB1和ZEN 2 種代謝產物含量為玉米霉變程度的表征指標，通過測定不同霉變程度玉米樣本的光譜信息和相應的AFB1和ZEN含量，建立基于較少校正集樣本和特征波長光譜信息的AFB1和ZEN含量偏最小二乘回歸（partial least squares regression，PLSR）預測模型，從而實現對玉米霉變程度的準確分析，以期為實現玉米霉變的在線、快速、精確檢測提供借鑒。

1 材料與方法

1.1 材料

玉米來自于2016年10月在洛陽當地收獲，品種為中單909，購于洛陽市中原農貿城，由實驗室自行培育出不同霉變程度的玉米。

1.2 儀器與設備

LHS-HC-100恒溫恒濕培養箱 上海資一儀器設備有限公司；IST 50-3810高光譜成像光譜儀 德國Inno-spec公司；康標達鏡頭 日本Computar公司；RK90000420108線性鹵素燈（500 W） 德國Esylux公司。

高光譜圖像采集系統如圖1所示，該系統主要包括高光譜成像光譜儀、康標達鏡頭、4 個500 W線性鹵素燈、實驗室自制傳送裝置和計算機等主要部件。其中高光譜成像光譜儀由攝像機和光譜儀2 部分組成，攝像機為CCD相機，光譜攝像儀為可見-近紅外光譜儀，光譜波長范圍為371.05～1 023.82 nm，可采集到的光譜數據點數為1 288 個，光譜分辨率為2.8 nm，傳送帶移動速率1.25 mm/s，曝光時間100 ms，樣本與鏡頭的距離為285 mm。高光譜成像儀通過USB 2.0接口數據線連接計算機，通過SICap-STVR V1.0.x軟件平臺驅動控制成像儀，及時記錄和存貯高光譜數據。

圖1 高光譜圖像采集系統Fig. 1 Schematic presentation of the hyperspectral image acquisition system

1.3 方法

1.3.1 玉米樣本制備

將新鮮玉米放入恒溫恒濕培養箱，在培養溫度30 ℃、相對濕度90%[22]的條件下進行培養。培養時間不同，玉米霉變程度也不同，本實驗用培養初始（新鮮玉米）、第2天、第4天、第6天、第8天標記玉米的不同霉變程度，5 個霉變程度的玉米看作5 個霉變等級，每個等級樣品各取50 個測試樣本，每個樣本含60 g（±0.1 g）玉米。

1.3.2 AFB1和ZEN含量的測定

為了減少AFB1和ZEN的損失，制備的樣本在完成光譜數據采集后立即按照GB/T 18979—2003《食品中黃曲霉毒素B1的測定》[23]和GB/T 5009.209—2008《食品中玉米赤霉烯酮的測定》[24]分別測定AFB1和ZEN含量，并將其測定結果作為建模的標準參考值。

1.3.3 高光譜圖像數據采集與校正

在進行高光譜圖像采集時，為了滿足實際玉米霉變檢測的需要，將玉米樣本散亂均勻平鋪在規格為10 cm×10 cm的白色載物盒中，然后將裝有玉米的載物盒放置在移動載物平臺上，依次對玉米樣本進行高光譜圖像采集，設定高光譜的掃描圖像大小為550×600像素。本實驗共采集到250 幅玉米樣本高光譜圖像。

同時，為了消除光源強度分布不均以及光譜暗電流噪聲的影響，對每個樣本按校正公式進行黑白校正，計算公式如下：

式中：R為校正后的樣本高光譜圖像；R0為樣本原始高光譜圖像；B為全黑標定圖像；W為全白標定圖像[25]。樣本的高光譜圖像采集在SICap-STVR V1.0.x軟件平臺上完成，實驗數據分析與處理在MATLAB 2014a（美國The Math Works公司）和遙感圖像處理平臺ENVI5.1（美國Boulder公司）兩個軟件上完成。

1.4 數據處理

1.4.1 原始光譜數據進行預處理

為了減少背景噪聲、雜散光等無用信息對原始光譜數據的干擾，需要對原始光譜數據進行預處理。本實驗對比標準化、變量標準化、多元散射校正（multiplicative scatter correction，MSC）、標準正態變量校正（standard normal variate，SNV）和卷積平滑（Savitzky-Golay smoothing，S-G）5 種預處理方法。為保證模型校正集最大程度地表征樣本均勻分布，提高模型穩定性，本實驗采用光譜-理化值共生距離（sample set partitioning based on joint x-y distance，SPXY）算法[26]劃分模型的初始校正集和預測集。在此基礎上，為降低或者消除校正集樣本間的共線性[27]，簡化模型運算量，采用SPXY算法結合PLSR法分析不同校正集樣本子集預測AFB1和ZEN含量的差異，從而進一步對劃分的初始校正集樣本進行優選。在采用均勻光譜間隔（uniform spectral spacing，USS）法對原始光譜變量進行初步篩選的基礎上對比連續投影算法（successive projections algorithm，SPA）[28]和競爭性自適應重加權算（competitive adaptive reweighted sampling，CARS）法[29]2 種特征波長提取方法，最大程度地剔除原始光譜矩陣中的冗余信息，提高模型精度，減少模型運算量。

1.4.2 模型建立與評價

霉變玉米的光譜信息與其AFB1和ZEN含量之間的關系屬于非確定性問題，采用回歸分析構造變量間的數理統計模型可用于該類問題的研究分析[30]。PLSR法集主成分分析、普通多元線性回歸和相關分析于一體，它在描述光譜矩陣X變量的同時也描述了指標矩陣Y變量，較好地解決了自變量的多重共線性問題，在光譜分析領域得到了廣泛應用。

2 結果與分析

2.1 不同霉變等級的玉米AFB1和ZEN含量

在進行AFB1和ZEN含量檢測時，對每類玉米樣品的50 個平行樣本進行測定，以獨立測定結果的絕對值差不超過算術平均值的10%為準，并將其平均值作為此類玉米樣品的實際指標值，結果如表1所示。由表1可知，第4、5個等級的玉米中，這2 類毒素含量已超過國家標準（AFB1≤20 μg/kg，ZEN≤60 μg/kg）。

表1 5 個霉變等級玉米樣品中AFB1和ZEN含量Table 1 AFB1 and ZEN values in 5 grades of moldy maize samples

2.2 光譜反射值曲線

提取每個玉米樣本在371.05～1 023.83 nm波長范圍內的高光譜圖像的平均光譜反射值，結果如圖2所示，371.05～480.55 nm和999.66～1 023.83 nm兩個波段受噪聲的影響較大，因此剔除了這2 個波段，最終采用的波段范圍為481.06～999.15 nm。由圖3可以看出，5 個等級玉米樣品光譜曲線變化趨勢基本相似，但不同霉變等級玉米的光譜值有差別，總體看來，5 個等級玉米樣本具有一定的可分性，這為預測建模提供了可能。

圖2 250 個玉米樣本的光譜反射值曲線Fig. 2 Spectral reflectance curves of 250 maize samples

圖3 5 類玉米樣本的平均光譜反射值曲線Fig. 3 Average spectral reflectance curves of five grades of maize samples

2.3 光譜數據預處理

隨機選取200 個樣本作為校正集，剩余的50 個樣本作為預測集。采用不同預處理方法對原始光譜數據進行預處理，并基于各預處理后的光譜數據建立PLSR模型，結果如表2所示。由表2可知，與基于原始光譜數據的預測結果相比，5 種預處理方法中，基于SNV建立的PLSR模型對這2 種毒素含量的預測效果最好，對應AFB1和ZEN含量的預測集相關系數和均方根誤差（R2pre，RMSEP）分別為（0.994 4，0.984 6）和（0.991 6，2.320 9），因此確定SNV為預處理方法。

表2 基于不同預處理方法的PLSR建模結果Table 2 Results of PLSR models based on different pretreatments

2.4 校正集樣本的優選

2.4.1 SPXY算法劃分樣本集數據

校正集樣本的劃分在一定程度上能夠決定所建模型的預測性能，本實驗采用SPXY算法對樣本集進行劃分。首先設定初始校正集樣本個數為200 個，剩余的50 個樣本為預測集樣本，SPXY算法的詳細步驟參考文獻[26]，劃分結果如表3所示，不同霉變等級的樣品間，校正集和預測集的樣本個數存在較明顯的差別。

表3 SPXY劃分樣本集結果Table 3 Calibration and prediction set partitioned by SPXY

2.4.2 SPXY算法優選校正集樣本

針對表3中初分的初始校正集樣本，采用SPXY算法進行進一步優選，以盡量有效地降低樣本間的共線性。確定樣本數N的范圍為80～200，步長為10，校正集樣本優選過程中，分別基于優選出的校正集樣本子集建立PLSR模型，其預測集的相關系數R2pre和RMSEP的變化如圖4所示。

圖4 R2pre和RMSEP隨校正集樣本數量變化曲線Fig. 4 Plots of R2pre and RMSEP values versus number of calibration set samples

由圖4可以看出，R2pre隨樣本數N的增加呈遞增趨勢，但總體變化不明顯（0.995～0.999），RMSEP曲線呈遞減趨勢，總體變化較明顯。圖4a中，對于RMSEP變化曲線，N在80～130范圍內取值時，RMSEP變化差異明顯，N在130～200范圍內取值時，RMSEP的變化趨于平緩，N取值為130 為該曲線的拐點，該點對應的R2pre和RMSEP為0.997 4和0.672 0，當N取值為190時，R2pre達到極大值0.997 6，RMSEP降為極小值0.641 4，與N取130 時對應的R2pre和RMSEP相比，R2pre僅增加0.000 2，RMSEP僅減少0.030 6，從數值上看，兩者差異較小，但校正集樣本數增加了60 個，因此，為簡化模型運算量，綜合圖4a中R2pre和RMSEP曲線的變化趨勢，對于AFB1，最終從初始校正集中優選出130 個樣本組成模型校正集。同樣根據圖4b，對于ZEN，最終從初始校正集中優選出140 個樣本組成模型校正集，對應的R2pre和RMSEP為0.998 7和0.862 1。

經SPXY算法劃分的預測集樣本和初始校正集樣本以及經SPXY算法優選的校正集樣本的AFB1和ZEN含量變化情況如表4所示。由表4可知，無論是對AFB1還是對ZEN，SPXY法優選出的校正集樣本的各指標值含量變化范圍在初始校正集范圍內，且標準差相近，證明SPXY算法優選后的樣本具有一定的代表性。

表4 不同樣本集玉米AFB1和ZEN含量分布Table 4 Distribution profiles of AFB1 and ZEN values in different sample sets

2.5 特征光譜的選擇

2.5.1 光譜數據的降維

高光譜圖像具有較高的光譜分辨率，容易對被測物進行分辨，但也會導致數據量的劇增和數據冗余。相鄰波段間差異很小，直接進行特征波段的提取可能會漏掉某些有用信息，因此，本實驗首先對光譜變量進行初降維。USS法通過控制步長選出間隔均勻、相關性低、信息量大的少量波段，從而可對光譜數據進行有效降維。在481.06～999.15 nm范圍內，分別以正整數2、3、4、5、6……20為間隔，從1 023 個經預處理后的光譜波段中抽取光譜組成新的光譜矩陣并建立PLSR模型，根據模型的和RMSEP，確定最佳的步長間隔。表5為利用不同步長間隔下抽取的光譜建立的PLSR模型對AFB1和ZEN含量的預測結果。

從表5可以看出，當間隔為8，波段為128 個時，模型對AFB1含量的預測效果較好，R2pre和RMSEP為0.997 4和0.673 4；當間隔為7，波段為147 個時，模型對ZEN含量的預測效果較好，R2pre和RMSEP為0.998 8和0.835 4。對比只采用SPXY算法優選校正集樣本而未對光譜變量進行篩選所建的PLSR模型結果可知，通過USS法進行變量篩選后所建模型對AFB1含量的預測精度變化不大，對ZEN含量的預測精度有一定的提高，但參與建模的變量數能顯著減少。

表5 基于不同間隔光譜數據建立PLSR模型預測結果Table 5 Predictive results of PLSR models with different numbers of spectral intervals

2.5.2 光譜特征波段的選取

圖5 基于SPA的變量篩選Fig. 5 Wavelengths selected by SPA

為進一步實現霉變玉米在線檢測的檢測時間短、準確率高的要求，本實驗在USS法篩選波長變量的基礎上，對比了SPA和CARS 2 種特征波長選擇方法。

在SPA中，最小的交互驗證RMSE對應的波長變量個數即為最終的選擇結果，指定波段數N的范圍為2～25，SPA的具體運算步驟可參考文獻[28]。由圖5a可以看出，從128 個波長變量中提取17 個波長時RMSE達到最小；由圖5b可以看出，從147 個波長變量中提取17 個波長時RMSE達到最小。表6為基于SPA提取的特征波長建立的AFB1和ZEN含量的PLSR模型預測結果。

在CARS算法中，設定蒙特卡羅采樣次數為50，CARS算法的具體運算步驟參考文獻[29]，如圖6所示，基于CARS算法針對AFB1和ZEN含量的變量篩選過程，對這2 種毒素指標值，當采樣次數分別為14 次和19 次時，對應的特征波長個數分別為18 個和19 個。

圖6 基于CARS的變量篩選Fig. 6 Wavelengths selected by CARS

表6 基于SPA特征波長提取方法的PLSR模型結果Table 6 Parameters of PLSR models with characteristic wavelengths extracted by SPA

表7 基于CARS特征波長提取方法的PLSR模型結果Table 7 Parameters of PLSR models with characteristic wavelengths extracted by CARS

對比表6和表7可知，基于SPA篩選出的波長所建立的PLSR模型對AFB1和ZEN含量的預測效果較CARS算法好。

2.6 模型分析

原始光譜數據經過SNV預處理后，采用SPXY算法對校正集樣本進行劃分與優選，運用USS法結合SPA對光譜數據進行篩選，然后建立基于優選后的校正集樣本及特征波長的PLSR模型。如圖7所示，波長優選后，PLSR模型的預測能力依舊表現良好，AFB1含量的預測精度（，RMSEP）為（0.997 3，0.681 5），與優選之前的預測精度（0.997 4，0.672 0）相近；ZEN含量的預測精度（，RMSEP）為（0.997 7，1.144 1），略低于優選前的精度（0.998 7，0.862 1）。

圖7 基于特征波段建立PLSR模型預測結果Fig. 7 Predictive results of PLSR models based on characteristic wavelengths

3 結 論

本實驗利用不同霉變時期的玉米光譜實驗數據，研究了高光譜技術用于霉變玉米中AFB1和ZEN含量檢測模型的構建，并得到以下結論：

1）對于AFB1和ZEN含量，經SPXY法分別優選出的130 個和140 個校正集樣本的指標值變化范圍在初始校正集范圍內，且標準差相近，說明SPXY算法優選的樣本具有一定的代表性。

2）基于SPXY法劃分的初始校正集樣本建立的PLSR模型預測精度明顯高于劃分之前的預測精度，且通過SPXY法對初始校正集樣本進一步優選，PLSR模型精度基本不變，但校正集樣本數減少為初始校正集的65%和70%，在保證模型穩健性的前提下，有效降低了校正模型的復雜性。

3）利用USS法結合SPA進行特征波段選擇，最終從1 023 個波段中分別篩選出17 個波段。分別建立2 種毒素值基于各自特征波段的PLSR模型，結果顯示AFB1含量預測精度（0.997 3，0.681 5）與優選前的預測精度（0.997 4，0.672 0）相差不大；ZEN含量預測精度（0.997 7，1.144 1）略低于優選之前的預測精度（0.998 7，0.862 1），R2pre從0.998 7降低到0.997 7，RMSEP從0.862 1上升到1.144 1，但波長減少為17 個，在保證模型精度的前提下，實現光譜數據的壓縮。

綜合上述研究結果，基于SPXY算法和SPA建立的高光譜檢測模型，可以實現對霉變玉米中AFB1和ZEN含量的準確預測。