基于高光譜成像技術的寧夏枸杞產地溯源鑒別

袁偉東，姜洪喆,2，楊詩雨，張 聰，周 禹，周宏平,2,*

（1.南京林業大學機械電子工程學院，江蘇 南京 210037；2.南京林業大學 林業資源高效加工利用協同創新中心，江蘇 南京 210037）

寧夏枸杞（Lycium barbarumL.）是茄科枸杞屬的多年生灌木，營養豐富且具有藥用價值，被稱為“藥食同源”水果[1]。近代藥理研究表明，枸杞子具有延緩衰老、降血糖、降血脂、降血壓、免疫調節、益精明目、滋陰補腎等多方面功效[2-3]。我國枸杞種植面積廣泛，主要分布在寧夏、內蒙古、新疆、甘肅等地，其中寧夏枸杞是唯一被收錄于《中華人民共和國藥典》的枸杞品種[4]。隨著人們生活水平的不斷改善以及保健意識的增強，寧夏枸杞因具有肉質飽滿、甘甜、色艷、含糖量低等優點受到消費者的喜愛。但品質優良的寧夏枸杞產量有限，導致一部分不良商販濫用其他產區的枸杞冒充寧夏優質枸杞以假亂真，嚴重擾亂寧夏枸杞交易市場。

傳統的枸杞產地溯源檢測方法包括高效液相色譜法[5]、碳同位素分析[6]和電子鼻[7]，而這些技術成本高昂、耗時且具有破壞性。因此，迫切需要一種經濟、高效、快速且無損的檢測方法[8]。近紅外光譜作為一種快速、無損檢測食品化學質量屬性的技術已被廣泛使用[9-10]。然而，近紅外光譜是一種點測量技術，在某些情況下由于樣品的異質結構，獲得的光譜信息可能不足或不能代表整個樣品。高光譜成像技術將光譜和成像技術相融合，同時獲取光譜和空間信息的特殊優勢備受食品行業研究人員的青睞[11-12]。目前，高光譜成像技術已廣泛應用于農產品的產地鑒別[13]、摻假[14]、含油率[15]等多方面的檢測。然而高光譜數據包含大量的共線性和冗余信息，導致模型將面臨更多的計算負擔和所開發模型的不穩定和復雜性。特征波長的選擇對于便攜式儀器的開發具有重要意義，特征波長選擇有利于挑選出共線性最小、冗余最少且包含主要信息的變量。選擇方法通常是從原始波長中選擇少數關鍵波長，以便提高模型的運行效率和可解釋性，去除非信息變量波段，產生更簡單、更穩健的校正模型。

在本研究中，采用400～1 000 nm波段范圍內的高光譜成像系統結合化學計量學方法對不同產地枸杞進行檢測鑒別。應用多種特征波長選擇方法選取特征波長建立判別模型并對比結果，為消除干擾信號和光散射的影響還采用了多種光譜預處理方法，旨在提供一種能夠快速、無損鑒別枸杞產地的檢測方法。

1 材料與方法

1.1 材料

優質枸杞樣品（寧杞1號）分別購于寧夏銀川、甘肅靖遠、內蒙古巴彥淖爾、青海西寧和新疆烏魯木齊5 個具有代表性的生產基地。

1.2 儀器與設備

用于測試的高光譜成像檢測系統由暗箱、升降平臺、成像光譜儀（SOC710VP）、4 盞75 W鹵鎢燈和一臺計算機組成，如圖1所示。光譜范圍374.98～1 038.79 nm，分辨率為4.68 nm，在256 個波段上進行成像。

圖1 高光譜成像系統示意圖Fig.1 Schematic diagram of hyperspectral imaging system

1.3 方法

1.3.1 樣品準備

挑選大小、顏色均勻的枸杞，并剔除表面有明顯缺陷、有蟲眼、干枯的劣質枸杞。為了減少溫度對樣品鑒別的影響，所有的枸杞樣品在進行成像前均貯存于4 ℃的冰箱中，貯存時間為12 h。每種產地枸杞均采集450 粒樣品，5 種產地共收集了2 250 粒枸杞樣品高光譜圖像。每種產地枸杞隨機選取300 個樣品分配進入訓練集，其余150 個樣品進入預測集。

1.3.2 高光譜圖像采集與校準

高光譜圖像采集軟件為Hyper Scanner_2.0.127，該系統由北京安洲科技有限公司提供。經過多次調試，樣品表面到高光譜鏡頭的距離為50 cm，積分時間設定為30 ms。為了捕捉到清晰的枸杞形狀，圖像分辨率設置為696×520像素，采集結束后得到一張696×520×256像素的三維圖像。

為消除光照不均勻和暗電流的影響，對原始圖像Ir進行反射率校正，根據式（1）計算Ic：

式中：Ic為校正后的高光譜圖像；Ir為原始高光譜圖像；Istd為校準后標準板的反射率；Is為標準板的反射率（約30%的反射率）。

1.3.3 反射光譜數據的提取和預處理

感興趣區域（region of interest，ROI）光譜信息提取的主要步驟如圖2所示。為了選擇單個枸杞樣品作為ROI，首先對原始高光譜圖像進行分水嶺分割、腐蝕、膨脹處理（圖2a），二值化生成掩膜并分割得到單個枸杞的原始光譜圖像（圖2b）。然后選取圖像中反射率最高和最低的兩個波段（885.89 nm和417.49 nm）進行波段運算，從而得到枸杞樣品與背景特征具有高對比度的灰度圖像（圖2c）。在生成的圖像中，通過將反射閾值設置為0.15得到二值化的掩膜圖像（圖2d），應用掩膜圖像將背景去除得到僅有枸杞的高光譜圖像（圖2e）。提取ROI的平均光譜作為對應樣品的反射光譜（圖2f）。由于336.2～400 nm和1 000～1 038.79 nm波段上的光譜強度較低，其數據波形抖動較大，包含很多噪聲信號，因此在后續分析中只保留400～1 000 nm之間的光譜數據，共計231 個波段。

圖2 高光譜信息提取過程Fig.2 Information extraction process from hyperspectral image

光譜數據中除了包含樣品的組分信息外，還夾雜著環境噪聲、雜散光、樣品背景等不相關信息。因此為了減小或消除這些無用信息的影響，提高檢測模型的精度和穩定性，建模前需要對原始光譜進行標準正態變量（standard normal variate，SNV）、SNV結合去趨勢、歸一化反射光譜（normalized reflectance spectrum，NR）、一階導數和二階導數等預處理。其中，一階導數和二階導數是基于二階多項式擬合和五點移動窗口的Savitzky-Golay平滑實現。

1.3.4 特征波長選擇方法

高光譜數據中包含大量的冗余信息，需要提取出最小冗余新變量消除數據間的共線性和重疊問題，提高模型的運行效率和模型的可解釋性。

連續投影算法（successive projections algorithm，SPA）是一種前向變量循環選擇方法。SPA在向量空間中使用簡單的投影操作，以選擇具有最小相關性的變量子集[16]。從一個波長開始，每次循環都計算它在未選中波長上的投影，最終將投影向量最大的波長引入循環后的波長組合。在本研究中，變量數的范圍設置為10～50，最小均方根誤差（root mean square error，RMSE）所對應的數量即為最終的選擇結果。

競爭性自適應重加權算法（competitive adaptive reweighted sampling，CARS）是一種將蒙特卡羅采樣與偏最小二乘回歸系數相結合的特征變量選擇方法，采用了“適者生存”的理論。其原理是通過自適應重加權采樣方法和指數衰減函數去除偏最小二乘回歸系數中絕對值最小的波長點，得到多個變量子集后，選擇交叉驗證均方根誤差最小的變量子集作為特征變量[17]。

粒子群優化算法（particle swarm optimization，PSO）是一種基于種群的全局智能優化算法。在PSO中，每個優化問題的解都被抽象作為多維空間一個沒有質量和體積的粒子。所有粒子都有一個評價函數決定的適應值（以F＝R2作為適應度函數，其中R2表示決定系數），粒子們通過相互作用發現復雜搜索空間中的最優區域[15]。本研究中，粒子種群數量大小設置為20，迭代次數設置為1 000。

迭代保留信息變量算法（iteratively retaining informative variables，IRIV）是一種基于二進制矩陣重排過濾器提出的特征波長選擇算法，將所有變量分為強信息變量、弱信息變量、無信息變量和干擾變量4 類[18]。通過多次迭代的方式，對于出現比例較多的變量賦予較高的權重，并去除無信息變量和干擾變量燈對模型無關的變量，最后保留有用的信息變量，經過反向消除獲得特征變量。隨著變量選擇方法的不斷發展，變量方法聯用組合使用也逐漸受到重視。梁琨等[19]通過CARS-IRIV篩選高光譜數據特征變量建立庫爾勒香梨可溶性固形物含量預測模型，在提高預測精度的同時簡化了模型的運算。CARS組合IRIV充分利用不同算法之間的互補性，構建更少、更有效的變量模型，對設備的開發和基礎建設具有重要的指導意義。

1.3.5 模型的建立與評價

偏最小二乘判別分析（partial leasts quares discriminant analysis，PLS-DA）是光譜分析中常用的定性建模方法，基于最小二乘算法開發的監督線性算法，適用于解決變量較多且存在多重線性的問題[20]。PLS-DA將光譜變量視為獨立輸入（X變量），并將不同的類別信息分配為1、2、3、4、5等（即Y變量）。為了校準最優PLS-DA模型，本研究采用10折交叉驗證法并依據最小均方根誤差確定PLS-DA模型的最佳潛在變量（latent variables，LVs），在全光譜模型和簡化模型中將LVs最大數量分別設置為40和30。

為了驗證上述建立模型的性能，本研究中采用訓練集、交叉驗證集和預測集的分類準確率（correct classification rate，CCR）評估分類模型的性能，CCR計算如式（2）所示：

式中：CCR為分類準確率（正確分類的樣本占全部樣本的比率）/%；TP為真陽性；TN為真陰性；FN為假陰性；FP為假陽性。

為了進一步評估所選模型的結果，基于混淆矩陣的靈敏度、特異性和Kappa系數評價模型的分類性能，計算分別如式（3）～（6）所示：

式中：Pe為預期的隨機準確率；P0為CCR；M為類別總數；i為類別數。

具體來說，靈敏度（又稱真正率或召回率）量化了模型識別正樣本的能力，特異性（又稱真負率）則是模型在正確識別偽樣本方面的效能。Kappa系數用于衡量分類模型的性能與隨機分類的性能之間的差異。當Kappa系數位于0.8～1之間，表示模型的預測性能與真實值之間幾乎完全一致，表明模型預測性能具有極高的可靠性和穩健性。

2 結果與分析

2.1 原始光譜

圖3顯示了從校準高光譜圖像中提取的5 種不同產地枸杞樣品的光譜曲線以及平均光譜曲線。5 種不同產地枸杞樣品之間的光譜曲線表現出相似的輪廓，主要源于枸杞內部組織共性，而反射率的強度差異主要受到其內部化學成分含量影響。多糖、黃酮、總糖和多酚是枸杞主要的活性成分，它們也是衡量其內部質量的主要特征指標，這些化學成分在不同產地的枸杞樣品中含量不同[21]。從圖3b可以看出，在可見光區域內的400～550 nm波段，不同產地枸杞的平均光譜反射強度幾乎相同，表明它們在藍色和綠色的顏色分量上幾乎沒有差異。然而，在550～700 nm的波段區域，各產地枸杞的平均光譜反射強度開始有所不同，且呈現出明顯的上升趨勢，表明紅色波段內存在的顏色分量差異較為明顯。在波長450 nm附近的反射峰主要與枸杞中的酚類物質阿魏酸相關[22]。550 nm波長處的波谷是總糖的有效波長[21]，而在波長560 nm附近的吸收峰與枸杞表面葉綠素和類胡蘿卜素的吸收帶相關[23]。在近紅外區域，黃酮的有效波段為800～900 nm，而在860 nm附近吸收峰為枸杞黃酮的有效波長[18]。在910～960 nm波段可歸因于水或碳水化合物的O—H拉伸模式的第二泛音[24]。由于5 個不同產地枸杞光譜之間存在很多交叉和重疊，因此將光譜與化學計量學方法相結合進行深度分析，并作出準確判別。

圖3 5 種不同產地的枸杞樣品光譜曲線Fig.3 Spectral curves of L.barbarum samples from five different geographical origins

2.2 基于全光譜的判別模型

對原始全光譜應用不同的預處理方法，并采用原始和預處理全光譜建立PLS-DA枸杞產地判別模型。為了確定最佳的預處理方法，本研究基于5 個產地的枸杞數據建立模型，結果如表1所示。結果表明，無論采用原始光譜還是預處理光譜PLS-DA模型，分類準確率均大于90%，表明所建立的模型可以輕松區分不同產地的枸杞樣本。在PLS-DA模型中，LVs的選擇會嚴重影響最終結果。當LVs過多時，會造成過擬合。相反，如果LVs的數量太少，會丟失一些有用的信息，降低模型的準確性。因此，通過交叉驗證計算最小預測殘差平方和以確定最優LVs[25]。圖4顯示了不同LVs對原始光譜PLS-DA判別模型性能的影響。還可以觀察到，NR和SNV預處理均提升了模型的性能，而SNV＋去趨勢、一階導數和二階導數略微降低了性能。其中，基于NR預處理全光譜的PLSDA模型表現最佳，最佳LV為30，訓練集分類準確率為95.5%，交叉驗證集分類準確率為91.9%，預測集分類準確率為93.1%。NR預處理可有效抑制部分光照差異的影響和消除無關信息，這一結論與徐新剛等[26]的研究結果一致。因此，在隨后的分析中最終選擇NR預處理光譜用于枸杞產地溯源鑒別。

表1 基于不同預處理方法全光譜PLS-DA模型性能Table 1 Performance of full-spectrum PLS-DA models based on different pre-processing methods

圖4 基于不同LVs原始光譜PLS-DA判別模型性能Fig.4 Performance of the raw spectral PLS-DA discriminant model based on different LVs

本研究對枸杞產地多元化鑒別需求進行深入分析，采用原始光譜和NR預處理光譜構建PLS-DA產地溯源判別模型，研究結果如表2所示?？梢园l現，隨著枸杞樣本產地的增加，無論采用原始光譜還是預處理光譜PLS-DA模型分類準確率總體呈下降趨勢。當模型輸入枸杞產地數量從2增加到5時，經預處理后模型的訓練集分類準確率從99.2%下降到95.5%，交叉驗證集分類準確率從94.8%下降到91.9%，測試集分類準確率從98.3%下降到93.1%。

表2 基于PLS-DA構建枸杞產地多元化判別模型Table 2 Multivariate identification model based on PLS-DA for geographical origins of L.barbarum

2.3 特征波長選擇

全光譜數據包含大量的冗余信息，特征波長的選擇有利于降低數據維度，這對于開發實時多光譜檢測系統非常需要[27]。在本研究中，SPA、CARS、PSO、IRIV和CARS＋IRIV對4 組光譜集篩選出的特征波長分布如圖5所示（為了便于觀察，同時將5 種特征選擇分布情況置于一張圖中，其中縱坐標值按比例增加）。當枸杞產地數量為2時，相較于其他產地數量所篩選的特征波段較少，其中5 種方法選擇的特征波長都相對分散不連續。當模型輸入為5 個產地枸杞樣本數據時所選的特征波長相對較多，但僅占全光譜的14.3%～42.4%?？梢园l現不同的特征變量選擇方法將選取不同數量的特征波長，因此確定最優的變量選擇方法對于構建高質量判別模型至關重要。枸杞產地溯源的進一步分類將使用基于選擇的特征波長進行PLS-DA建模。

圖5 不同方法選擇的特征波長分布Fig.5 Distribution of characteristic wavelengths selected by different methods

2.4 簡化模型

基于不同特征波長選擇方法選取的特征波長構建多元化枸杞產地溯源PLS-DA判別模型，結果如表3所示。與全光譜PLS-DA模型性能相比，結果有所下降。由于特征波長極大簡化了分類模型，可能造成這種性能下降。整體來看，基于SPA篩選特征波長建立的模型性能表現較差，主要原因是所篩選的波長主要集中在900～1 000 nm，無法提供全部的有效信息。雖然使用CARS和IRIV選擇的波長建模取得了優異的結果，由于選擇的波長較多不利于光譜檢測系統的開發，對于PSO-PLS-DA和CARS＋IRIV-PLS-DA模型表現出令人滿意的結果，所選的特征波長僅占全波長的15.6%～27.7%，在二元分類模型中預測集分類準確率分別為96.0%和97.7%，在三元分類模型中預測集分類準確率分別為90.0%和90.9%，在四元分類模型中預測集分類準確率分別為86.7%和89.2%，在五元分類模型中仍取得84.1%和87.1%的分類結果。

表3 基于不同波長選擇方法的最優PLS-DA建模結果Table 3 Results of optimal PLS-DA models based on different wavelength selection methods

為了進一步探索CARS＋IRIV-PLS-DA模型鑒別枸杞產地溯源的能力，圖6給出了該簡化模型預測集的混淆矩陣、靈敏度和特異性以及Kappa系數計算結果。在混淆矩陣中，縱坐標表示實際類，橫坐標表示預測類。主對角線內的值表示正確分類的樣本，主對角線外的值表示錯誤分類的樣本。結果表明，內蒙古和新疆枸杞樣本發生錯誤分類數量最少（靈敏度均大于94%），可能由于內蒙古和新疆枸杞相較于其他產地枸杞含糖量高，具有較好的區分性。還可以觀察到寧夏枸杞易被錯誤分類成內蒙古枸杞，在二元分類模型中由于模型簡單寧夏枸杞的識別率高達96.7%。隨著輸入產地數量的增加，寧夏枸杞識別率和Kappa系數整體呈下降趨勢，在五元分類模型中寧夏枸杞仍取得了82.7%的識別率。4 組簡化模型的Kappa系數均超過0.83，說明分類模型具有較強的穩定性和魯棒性。這些結果表明，CARS＋IRIV-PLS-DA模型在沒有任何化學或物理信息的情況下識別寧夏枸杞產地溯源具有巨大的潛力。

3 結論

鑒別枸杞產地的傳統方法耗時且具有破壞性和主觀性，因此，本研究旨在開發一種基于高光譜成像技術，從而快速、無損識別寧夏枸杞產地的檢測方法。通過高光譜成像采集寧夏、甘肅、內蒙古、青海和新疆共2 250 個枸杞樣本，并從ROI提取光譜數據。研究發現相比于SNV、SNV結合去趨勢、一階導數和二階導數，NR可以更好地降低光譜噪聲和散射效應。為了降低數據維度并進一步減少建模時間，采用SPA、CARS、PSO、IRIV和CARS＋IRIV選擇特征波長，并基于特征波長建立PLS-DA判別模型。隨著枸杞產地數量的增加，模型性能呈下降趨勢。當僅輸入兩個枸杞產地數量時，全光譜模型分類準確率高達98.3%，鑒于實用性最佳簡化模型CARS＋IRIV-PLS-DA分類準確率高達97.7%。當輸入為5 個枸杞產地數量，簡化模型CARS＋IRIV-PLS-DA仍能獲得87.1%的分類準確率和0.839的Kappa系數。綜合分析表明，高光譜成像技術（400～1 000 nm）結合化學計量學方法可以作為一種快速、無損的檢測方法鑒別寧夏枸杞的真偽性。