劉子健，顧佳盛，周 聰，王游游，楊 健，黃 俊，王宏鵬,*，白瑞斌,*

（1.浙江科技大學生物與化學工程學院，浙江杭州 310023；2.中國中醫(yī)科學院中藥資源中心道地藥材品質(zhì)保障與資源持續(xù)利用全國重點實驗室，北京 100700；3.江西省道地藥材質(zhì)量評價研究中心，江西南昌 330000）

山楂（Crataegus pinnatifida）是薔薇科山楂屬植物，是典型的“藥食同源”植物[1-3]，在我國廣泛分布于吉林、遼寧、河北、河南、山東、山西等地區(qū)。我國山楂年產(chǎn)量超過150 萬噸，市場前景廣闊，但由于不同產(chǎn)地山楂的各類營養(yǎng)成分含量存在差異[4-5]，因此在價格上也有所區(qū)分，而當今山楂市場上產(chǎn)地混用、以次充好等現(xiàn)象屢見不鮮，使許多消費者上當受騙，這些現(xiàn)象嚴重破壞了市場秩序。因此，目前市場亟需一種能夠快速準確對山楂進行產(chǎn)地溯源的方法。

目前對于各類農(nóng)產(chǎn)品產(chǎn)地的溯源主要是依靠化學分析技術，如高效液相色譜技術（High Performance Liquid Chromatography，HPLC）[6]、氣相-質(zhì)譜聯(lián)用技術（Gas Chromatography-Mass Spectrometry，GC-MS）[7]、超高效液相色譜（Ultra Performance Liquid Chromatography，UPLC）[8]等，這些方法通常在測量前需要對樣本進行粉碎或勻漿處理，并使用有機溶劑對樣本中的化學成分進行萃取，這一過程不但會損壞樣本，同時有機溶劑還可能會對環(huán)境造成污染。與之相比，高光譜成像技術是一種基于非常多窄波段的影像數(shù)據(jù)技術，可以在樣本完好的情況下對其進行定性或定量分析[9-10]，具有快速、無損、無污染檢測的特點。目前，基于高光譜成像技術對各類水果進行產(chǎn)地識別的研究已有很多：張立欣等[11]基于高光譜成像技術，結合多種機器學習算法實現(xiàn)對不同產(chǎn)地的蘋果進行區(qū)分，其最佳識別模型為MSC-CARSSVM，預測集準確率達到100%；Sun 等[12]采集了不同產(chǎn)地水蜜桃的高光譜數(shù)據(jù)，并利用HPLC 分析解釋樣本的“高光譜指紋”，最終建立的分類器準確率達到99.3%。以上研究表明，高光譜成像技術在水果產(chǎn)地識別領域具有廣闊的應用前景。然而，目前使用高光譜成像技術對山楂進行產(chǎn)地溯源卻鮮有報道。

基于上述背景，為滿足市場需求，本文旨在探究高光譜成像技術在山楂產(chǎn)地識別中的應用及不同采樣方向?qū)τ谀Ｐ头诸愋阅艿挠绊?，利用高光譜成像系統(tǒng)（410～2500 nm），分別采集山楂樣本果梗面、側(cè)面及底面的光譜數(shù)據(jù)，結合多種機器學習算法分別建立產(chǎn)地識別模型，最終實現(xiàn)基于高光譜成像技術對山楂進行產(chǎn)地溯源的目的。

1 材料與方法

1.1 材料與儀器

山楂樣品 山楂樣品均采自2022 年10 月至12 月，其中山東省2 批、山西省3 批、遼寧省2 批、河北省3 批、河南省1 批。每個批次隨機選擇80～100 粒品相完好、大小相近的山楂，最終共采集900 粒樣品。使用干布擦拭樣品表面殘留泥土，然后于4 ℃環(huán)境中冷藏保存，便于后續(xù)圖像采集。

高光譜成像系統(tǒng) 高光譜成像系統(tǒng)由成像模塊、鹵鎢燈、水平移動平臺和Teflon 白板組成。其中成像模塊包含兩個鏡頭，分別為SN0605 VNIR、N3124 SWIR；鹵鎢燈 挪威Norsk Elektro Optikk公司；水平移動平臺 立陶宛Standa translation stage 公司；Teflon 白板 中國雙利合譜公司。

1.2 實驗方法

1.2.1 高光譜數(shù)據(jù)采集 樣本圖像采集前，關閉環(huán)境燈光，打開鹵鎢燈并對高光譜成像系統(tǒng)進行預熱。為探究擺放方式對山楂產(chǎn)地識別模型的影響，將山楂樣本以果梗朝上（G）、側(cè)面朝上（C）和底面朝上（D）三種方式擺放（圖1），分別拍攝圖像。采集圖像時，將15～20 粒樣本放置在水平移動平臺上，在樣本排列末端放置Teflon 白板，分別采集三個方向的圖像數(shù)據(jù)，光譜儀鏡頭與樣品的距離為25 cm，平臺移動速度為1.5 mm/s，SN0605 VNIR 鏡頭的積分時間為3500 μs，幀時間為18000，光譜范圍410～990 nm，共108 個波段；N3124 SWIR 鏡頭積分時間為4500 μs，幀時間46928，光譜范圍950～2500 nm，共288 個波段，兩個鏡頭的光譜分辨率均為6 nm。

圖1 三種樣本擺放方式Fig.1 Placement methods of three samples

為減小環(huán)境以及儀器對圖像數(shù)據(jù)的影響，在圖像采集完成后使用軟件（HySpex RAD，Norsk Elektro Optikk，挪威）對原始光譜數(shù)據(jù)進行RAD 校正。隨后進行黑白板校正以消除空氣等外界因素對圖像的影響并得到相對反射率，相對反射率計算公式如下：

式中：R 表示相對反射率；Rraw表示原始反射率；Rw表示白板反射率，即Teflon 白板反射率（反射率接近1）；Rd表示黑板反射率（反射率接近0）。

校正完成后，使用軟件ENVI 5.3 在圖像中手動選取感興趣區(qū)域（ROI），對于不同拍攝方向的樣本圖像，分別取其相應部位（即果梗面、側(cè)面和底面）作為ROI，以ROI 平均反射率作為樣本的光譜值。手動合并兩個鏡頭得到的光譜數(shù)據(jù)，最終得到包含396個波段反射率的數(shù)據(jù)集。將樣本按照7:3 的比例隨機劃分為訓練集和預測集，用于后續(xù)分類建模。

1.2.2 主成分分析 主成分分析（Principal Components Analysis，PCA）是一種常用的聚類分析方法[13]，其基本原理是通過線性變換的方式，將原始數(shù)據(jù)轉(zhuǎn)換成一組線性無關的“特征”，而每個“特征”稱為“主成分”，是一種通用的統(tǒng)計方法。本研究在得到樣本光譜原始數(shù)據(jù)后，首先利用PCA 方法，對樣本數(shù)據(jù)進行初步的可視化分析。

1.2.3 光譜數(shù)據(jù)預處理方法 受儀器和光散射等因素的影響，樣本的原始光譜數(shù)據(jù)中存在很多噪聲[14]，會影響后續(xù)的建模分析，因此對原始光譜數(shù)據(jù)進行預處理可以有效提高后續(xù)建模的準確率[15]。本研究為消除噪聲的影響，分別采用多元散射校正（Multiplicative Scatter Correction，MSC）、一階導數(shù)（Derivative，D1）、SG 平滑（Savitzky-Golay，SG）和標準正態(tài)變量變換（Standard Normal Variate Transformation，SNV）四種方式對原始光譜數(shù)據(jù)進行預處理，再使用預處理后的數(shù)據(jù)進行分類建模。

1.2.4 特征波長提取方法 由于高光譜圖像包含波段多，數(shù)據(jù)維度高，在分類建模時經(jīng)常會面臨Hughes 現(xiàn)象和維數(shù)災難問題[16]，而特征波長提取是一種常見的降維方式，可以有效降低模型復雜度，從而提升模型運算速度[17]。本研究在建立全波段分類模型后，為降低模型復雜度，分別采用連續(xù)投影算法（Successive Projections Algorithm，SPA）和競爭性自適應重加權采樣算法（Competitive Adaptive Reweighted Sampling Algorithm，CARS）對原始光譜數(shù)據(jù)進行特征波長提取，然后基于特征波長數(shù)據(jù)建立分類模型，為山楂專屬小型高光譜設備的開發(fā)提供參考。

1.2.5 分類模型的建立 對原始數(shù)據(jù)進行預處理或特征波長提取后，基于處理得到的數(shù)據(jù)，分別采用以下方法建立分類模型，并綜合對比各項評估指標以篩選出最優(yōu)模型。偏最小二乘判別分析（Partial Least Squares Discriminant Analysis，PLS-DA）是一種常見的線性判別模型。偏最小二乘法的原理是通過協(xié)方差極大化準則，分解自變量數(shù)據(jù)和因變量數(shù)據(jù)，建立相互對應的回歸關系方程[18]。而PLS-DA 則是基于偏最小二乘法建立自變量和分類變量之間的回歸模型，提取出與分類相關的特征變量，實現(xiàn)樣本的類別預測[19]。本研究將特征變量的個數(shù)控制在10～15 個，以防止模型過擬合。

支持向量機（Support Vector Machine，SVM）的基本思想是將樣本特征數(shù)據(jù)映射到n 維空間中，n 的大小取決于核函數(shù)和樣本特征維數(shù)，然后在空間中構造最優(yōu)的分類超平面[20]。SVM 在樣本數(shù)量較少的情況下也能取得較好結果，具有優(yōu)秀的泛化能力。核函數(shù)是支持向量機映射數(shù)據(jù)的重要手段[21]，本研究使用徑向基核函數(shù)的支持向量機進行分類建模，并通過網(wǎng)格搜索篩選最優(yōu)模型參數(shù)（懲罰系數(shù)和核寬度）。

隨機森林（Random Forests，RF）由預先設定數(shù)量的分類決策樹組成。決策樹可以表示樣本屬性與其特征值之間的映射關系，樹中每一個節(jié)點表示對象屬性的判斷條件[22]，隨機森林通過對所有決策樹的預測值進行平均或投票得到最終結果[23]。本研究通過網(wǎng)格搜索找出最佳的決策樹數(shù)量和深度。

1.2.6 模型評估標準 模型建立完成后，分別通過以下指標篩選出最優(yōu)模型：準確率（Accuracy）是分類問題最常用的評價指標；精確率（Precision）和召回率（Recall）則反映了模型對于正例的敏感程度，三個指標計算公式如下：

式中：TP 表示真陽性樣本個數(shù)；FP 表示假陽性樣本個數(shù)；TN 表示真陰性樣本個數(shù)；FN 表示假陰性樣本個數(shù)。

混淆矩陣可以直觀呈現(xiàn)模型預測結果，通常由m 行m 列組成（m 取決于樣本類別總數(shù)），每一列代表模型預測值，每一行則代表真實值。本研究通過建立混淆矩陣，綜合對比模型指標，篩選出最優(yōu)分類模型。

1.3 數(shù)據(jù)處理

高光譜圖像的采集和校正分別使用儀器自帶軟件HySpex GROUND 和HySpex RAD 實現(xiàn)；樣本反射率數(shù)據(jù)使用軟件ENVI5.3 進行提取；PCA 分析、原始數(shù)據(jù)預處理、特征波長提取及分類模型的建立均基于Pycharm（Python 3.10）軟件實現(xiàn)；光譜曲線繪圖基于Origin 2023b 軟件完成。

2 結果與分析

2.1 原始光譜曲線分析

在進行分類建模前，首先對各產(chǎn)區(qū)樣本的光譜特征進行分析并探究部分特征峰的成因，不同產(chǎn)區(qū)樣本的平均光譜曲線如圖2 所示。對比發(fā)現(xiàn)不同產(chǎn)區(qū)山楂樣品的平均反射率總體趨勢相似；但是同產(chǎn)區(qū)山楂平均反射率在不同數(shù)據(jù)集（C、G 和D）上有所不同，這可以歸因于樣品表面信息的差異。另外，不同產(chǎn)區(qū)山楂樣品的反射率數(shù)值存在一定差異，這些差異主要與樣品的表面信息（如果皮、果斑顏色）和品質(zhì)特性有關，其中山東產(chǎn)區(qū)的山楂在400～800 nm 波段下的反射率明顯高于其他產(chǎn)區(qū)，區(qū)域特征較為明顯，根據(jù)楊曉寧等[4]的研究報道：相比于其他產(chǎn)區(qū)，山東產(chǎn)區(qū)山楂的有機酸含量較高，這與上述現(xiàn)象相吻合。不同產(chǎn)區(qū)山楂在600～700 nm 處的吸收峰略有不同，但總體趨勢相似；對于短波紅外波段（short wave infrared，SWIR），各產(chǎn)地反射率曲線趨勢相近，但在1000～1200 nm 處的吸收峰有所區(qū)分。對不同波段下的吸收峰進行分析，700～800 nm 處的吸收峰可歸因于樣本中的葉綠素[24]；970 nm 附近的吸收峰可能是水中O-H 鍵的伸縮振動造成[25-26]；1200 nm 附近的吸收峰可能與C-H 的第二拉伸泛音有關，可歸因于碳水化合物和脂肪[27]，總體而言，各產(chǎn)地樣本所含化學成分種類相似，但具體含量存在差異，這與張悅等[28]報道的不同產(chǎn)地陳皮光譜曲線規(guī)律一致。

圖2 不同產(chǎn)地在VNIR 和SWIR 波段下的平均反射率曲線Fig.2 Average reflectance curves of different origins in VNIR and SWIR bands

對比各數(shù)據(jù)集的平均反射率曲線（圖2），發(fā)現(xiàn)G 數(shù)據(jù)集在700～1000 nm 處反射率略高于其他數(shù)據(jù)集，而此波段反射率與樣品水分及葉綠素含量密切相關，因此推測山楂樣本不同部位所含成分略有不同。山東與遼寧產(chǎn)區(qū)樣品的平均反射率在三個數(shù)據(jù)集上都表現(xiàn)出了較大差異（山東產(chǎn)區(qū)樣品反射率較高，而遼寧產(chǎn)區(qū)樣品則偏低），說明兩組樣品差異明顯。光譜平均反射率曲線雖然展現(xiàn)出樣本的部分差異，但是僅憑這些特征很難對樣本進行產(chǎn)地溯源。綜上所述，有必要建立分類模型以挖掘樣品光譜數(shù)據(jù)的潛在特征。

2.2 樣本數(shù)據(jù)PCA 分析

使用主成分分析（PCA）對三個數(shù)據(jù)集進行初步的可視化分析，繪制的PCA 得分圖見圖3，保留了前兩個主成分。初步分析發(fā)現(xiàn)，無監(jiān)督模型分類效果并不好，三個數(shù)據(jù)集前兩個主成分能解釋的方差占比之和在75%左右。山東與遼寧產(chǎn)地的樣本區(qū)分相對較好，這與原始光譜分析時得出的結論相符。對于大部分樣本，使用無監(jiān)督算法進行分類的效果并不理想，因此后續(xù)還需要采用PLS-DA、SVM 和RF 方法進行有監(jiān)督分類建模。

圖3 原始數(shù)據(jù)PCA 得分圖Fig.3 Raw data PCA score plot

2.3 基于全波段的建模分析

2.3.1 預處理及分類建模方法篩選 為篩選出最佳預處理和分類建模方法，分別采用4 種預處理方法和3 種分類建模方法建立模型，以樣本底面數(shù)據(jù)為代表，各模型分類準確率見表1。對比四種預處理數(shù)據(jù)分類模型準確率可以發(fā)現(xiàn)，引入預處理方法之后，大部分模型的分類精度得到了提高，而D1 對于三種分類模型（PLS-DA、SVM 和RF）均為最優(yōu)預處理方式。對比三種不同模型（PLS-DA、SVM 和RF）分類準確率，發(fā)現(xiàn)無論采用哪種預處理方式，采用RF 建立的分類模型雖然有較高的訓練集準確率，但是預測集準確率一般；采用PLS-DA 和SVM 建立的分類模型訓練集和預測集準確率良好，其中以SVM 模型分類準確率最高。綜上所述，對于底面數(shù)據(jù)，D1 為最佳預處理方式，采用SVM 建立的分類模型分類準確率高，且具有優(yōu)秀的穩(wěn)定性和泛化能力。為進一步驗證結論，分別使用C 和G 數(shù)據(jù)集進行建模對比，均呈現(xiàn)相同的規(guī)律，故判斷D1 為最優(yōu)預處理方式，SVM 為最佳分類建模算法，后續(xù)均采用D1-SVM（經(jīng)D1 預處理后建立的SVM 模型）方式進行分類建模。

表1 不同預處理分類模型準確率Table 1 Accuracy of different preprocessing classification models

2.3.2 不同采樣方式分類建模分析 本研究為探究不同采樣方向?qū)δＰ头诸惤Y果的影響，分別收集了樣本側(cè)面朝上（C）、果梗面朝上（G）和底面朝上（D）的高光譜圖像。同時為模擬實際應用時隨機拍攝到的高光譜數(shù)據(jù)，將三個數(shù)據(jù)集進行等比混合建立一個新數(shù)據(jù)集（R），使用四個數(shù)據(jù)集分別進行分類建模，建模方法均采用D1-SVM，綜合對比各項指標篩選出最優(yōu)模型。各模型分類準確率結果見表2。

表2 不同方向數(shù)據(jù)分類模型準確率Table 2 Accuracy rate of data classification model in different directions

對于使用R 數(shù)據(jù)集建立的分類模型，其準確率較高（100%，96.7%），根據(jù)圖4 并由公式（3）和公式（4）計算得出，不同產(chǎn)區(qū)的精確率和召回率均超過90%。對比四個數(shù)據(jù)集模型的準確率可以發(fā)現(xiàn)，三種單面數(shù)據(jù)集（C、D 和G）模型準確率均高于使用R 數(shù)據(jù)集建立的模型，這說明對于山楂樣本，在高光譜數(shù)據(jù)采集時保持樣品方向一致可以有效提高分類模型準確率，這一規(guī)律與Mansuri 等[29]在玉米真菌感染檢測中的發(fā)現(xiàn)一致。橫向?qū)Ρ菴、G 和D 三個模型，其中使用D 數(shù)據(jù)集建立的分類模型準確率最高，訓練集和預測集準確率均達到100%，各產(chǎn)區(qū)樣本全部預測正確。為避免過擬合現(xiàn)象，對D-D1-SVM 模型進行十折交叉驗證，其平均準確率為98.8%。綜上所述，D-D1-SVM 模型對于不同產(chǎn)區(qū)山楂的分類效果最優(yōu)。

圖4 全波段模型混淆矩陣Fig.4 Confusion matrix of full band model

2.4 基于特征波長的建模分析

2.4.1 特征波長的選擇 為篩選出最佳特征提取方法，分別使用2 種提取方式提取4 個數(shù)據(jù)集的特征波長，最終得到的波長見表3 及圖5。對比兩種方法提取得到的特征波長數(shù)量發(fā)現(xiàn)，使用SPA 提取出的特征波長數(shù)量明顯少于CARS，進一步觀察特征波長分布（圖5），發(fā)現(xiàn)使用SPA 提取出的特征波長分布均勻，各個波段均有涉及；而CARS 提取的特征波長分布較為集中，主要分布于750、2000 及2250 nm處的三個特征峰。觀察各組特征波長重合的部分，發(fā)現(xiàn)750、1700 和2200 nm 附近的重合波長較多，說明這三處吸收峰可能包含不同產(chǎn)區(qū)樣本的差異信息。對這些特征峰進行深入分析，700～800 nm 處的吸收峰來自于樣品內(nèi)部的葉綠素，也受樣品的外部顏色特征影響；1700 nm 附近的吸收峰可歸因于酰胺基團[30]；2200 nm 處的吸收峰為C-H 和C-O 的聯(lián)合吸收峰[31]。

表3 不同方法提取特征波長數(shù)量Table 3 Number of feature bands extracted by different methods

圖5 不同數(shù)據(jù)集特征波長Fig.5 Feature bands of different datasets

2.4.2 特征波長建模分析 使用4 個數(shù)據(jù)集的特征波長分別建立SVM 模型，其準確率見表4。觀察發(fā)現(xiàn)使用SPA 篩選特征波長建立的模型分類準確率優(yōu)于CARS，這一現(xiàn)象在G 和D 數(shù)據(jù)集上尤為明顯。綜合考慮波長數(shù)量和模型準確率，SPA 篩選的波長數(shù)量更少，模型復雜度較低，且準確率更高。與本研究得到的結果不同，李濤等[32]在基于特征波段建立紅景天分類模型時，發(fā)現(xiàn)CARS 為最佳特征波段提取方法，這說明對于不同的檢測對象，應當選用不同的特征提取方法，而對于山楂樣本，SPA 相比于CARS 特征波長提取效果更好。

表4 特征波長建模準確率Table 4 Accuracy of feature band model

采用SPA 提取特征波長的分類模型預測集混淆矩陣見圖6，對比四個數(shù)據(jù)集的準確率（表4）看出，R-SPA 模型預測集準確率為87.8%，根據(jù)其混淆矩陣（圖6d）并由公式（3）和公式（4）計算得出，模型對于河北產(chǎn)區(qū)的精確率和召回率僅為79.2%和82.4%，分類能力一般。而C-SPA、G-SPA 和D-SPA 三個模型準確率均超過90%（分別為90.3%、91.5%和93%），這一現(xiàn)象再次證明在高光譜數(shù)據(jù)采集時，保持樣品方向一致可以有效提高分類模型準確率。綜合對比所有模型，D-SPA 模型擁有最高的分類準確率，訓練集和預測集準確率分別為95.2%和93%，根據(jù)其混淆矩陣（圖6c）并由公式（3）和公式（4）計算得出，模型對于各產(chǎn)區(qū)的精確率和召回率均超過90%（其中山東產(chǎn)區(qū)精確率和召回率最低，分別為91.6%和90%）；且這一模型涉及的特征波長數(shù)量最少，在保證分類準確率的情況下?lián)碛休^低的模型復雜度。

圖6 特征波長模型混淆矩陣Fig.6 Confusion matrix of feature band model

綜上所述，采集高光譜數(shù)據(jù)時保持樣品擺放方式一致有助于提高模型分類準確率。采用SPA 提取特征波長建立的產(chǎn)地分類模型復雜度較低且準確率良好?？梢栽诓ㄩL數(shù)量有限的情況下對山楂產(chǎn)地進行判別，為后續(xù)山楂專屬小型化高光譜設備的開發(fā)提供了方法參考。

綜合考慮全波段模型和特征波長模型的分類結果，發(fā)現(xiàn)采集樣本光譜數(shù)據(jù)時，樣本的擺放方式會影響后續(xù)分類建模準確率。無論全波段還是特征波長模型，使用D 數(shù)據(jù)集建模分類效果都明顯優(yōu)于R 數(shù)據(jù)集（提高了約5%），相對于C 和G 數(shù)據(jù)集也有所提高。觀察山楂樣品的外部特征，發(fā)現(xiàn)樣品底面存在萼片部位，結合寧素云等[33]的研究報道：山楂不同部位的化學成分含量存在差異，推測不同產(chǎn)地山楂其萼片部位各成分含量的差異相比于其他部位更大，進而導致分類特征更加明顯。

3 結論

本研究基于高光譜成像技術建立了山楂產(chǎn)地識別模型。為探究樣本拍攝方向?qū)Ψ诸惤Y果的影響，采集了山楂樣本三個不同方向（C、G 和D）的光譜數(shù)據(jù)，分別使用偏最小二乘判別分析（PLS-DA）、支持向量機（SVM）和隨機森林（RF）三種方法建立模型，通過對比模型分類準確率得到最優(yōu)建模方法，最終成功區(qū)分了5 個不同省級產(chǎn)區(qū)的山楂，為山楂無損檢測設備的開發(fā)提供了參考。經(jīng)過對比篩選發(fā)現(xiàn)，一階導數(shù)（D1）為最優(yōu)預處理方式，SVM 為最優(yōu)建模算法；使用連續(xù)投影算法（SPA）提取特征波長數(shù)量少且分類模型準確率高。全波段最優(yōu)建模方法為D-D1-SVM，訓練集和預測集準確率均達到100%；特征波長最優(yōu)建模方法為D-SPA-SVM，訓練集和預測集準確率分別為95.2%和93%。本研究證明基于高光譜成像技術對山楂產(chǎn)地進行溯源是可行的，為維護山楂市場秩序提供一種新的識別方式；同時驗證高光譜圖像采集方向會對檢測結果產(chǎn)生影響，為后續(xù)開發(fā)山楂專屬高光譜檢測設備提供理論依據(jù)和參考。

? The Author(s) 2024.This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by-nc-nd/4.0/).