間隔影響分析波長選擇算法在近紅外光譜鑒別貝類毒素中的應用

姜微，劉瑤，劉忠艷，曾紹庚，熊建芳，喬付

1(嶺南師范學院 計算機與智能教育學院，廣東 湛江，524048)2(嶺南師范學院 電子與電氣工程學院，廣東 湛江，524048)

在海洋生物中，貝類的主要食物為藻類、原生動物等一些浮游生物，在原甲藻和鰭藻等藻類中易產生海洋生物毒素——腹瀉性貝類毒素(diarrhetic shellfish poisons，DSP)[1-2]。貝類屬于非選擇濾食性動物，在生長過程中，難免濾食有毒藻類食物，由此引起DSP在其體內長期融積富集，當人們誤食染毒的貝類后就會對身體健康產生危害，引起中毒現象，甚至會增加患結腸癌等腫瘤的風險[3-5]。由于DSP中毒現象分布廣泛，人體損傷性大，應對食用貝類加強毒素檢測力度，確保貝類食用安全，因此研究高效高質的貝類毒素檢測方法是十分必要的。

目前腹瀉性貝類毒素檢測方法主要有小鼠生物測定法、酶聯免疫吸附法、高效液相色譜法等[6-8]。小鼠生物測定法是應用最多的貝類毒素檢測方法，但該方法受小鼠個體影響，檢測結果靈敏度不高，偏差較大。酶聯免疫吸附法以抗原與抗體特異性反應為基礎，交叉反應時產生較大偏差。高效液相色譜法具有靈敏度和準確度高等優勢，但該方法所用儀器昂貴，且前處理過程比較復雜，需要經過專業培訓的人員，比較適用于實驗室和較高精度需求的樣本檢測。這些方法均有各自的優勢，但都存在費時費力，對樣品具有破壞性等問題，無法實現快速無損分析。因此，開發高效、快速、無損的貝類毒素檢測技術具有重要的現實意義。

近紅外光譜屬于分子振動光譜，主要是對含氫基團(C—H、N—H、O—H、S—H等)振動的倍頻和合頻進行吸收。由于不同的物質含有不同的氫基團，這些基團對近紅外光吸收波長均有所不同，通過采集樣品的近紅外光譜，可以得到樣品中含氫基團的特征信息，從而得出樣品間的差異特征。近紅外光譜分析技術正是根據不同的光譜特征信息及傳統化學分析方法測定的樣品性質或數據，通過化學計量學方法，將光譜與數據進行關聯，從而實現對未知樣品進行定性和定量分析[9-10]。近紅外光譜技術因具有快速無損、綠色環保等優點，近年來在食品的真偽鑒定[11]、產地鑒別[12]、等級判別[13]以及微生物檢測[14]等領域得到了廣泛應用。在貝類檢測方面，黃冠明等[15]對牡蠣(Crassostreaangulata)樣本的蛋白質、水分、脂肪進行分析，LIU等[16]提出基于高光譜技術對菲律賓蛤仔(Ruditapesphilippinarum)重金屬污染進行檢測，劉忠艷等[17]提出基于近紅外光譜和多層感知機對腹瀉性貝毒的檢測方法。這些研究成果為貝類毒素高效、快速、無損檢測提供了技術可行性。

由于光譜采集過程中存在環境、樣本、操作人員等的影響，所以光譜往往存在譜峰重疊、基線漂移等干擾因素。為了消除這些干擾，需要結合化學計量學方法對光譜數據進行預處理，建立貝類毒素的準確鑒別模型。同時近紅外光譜波段數量較大，波段間的相關性較高，一定程度上存在數據冗余等問題，因此數據降維一直是近紅外光譜數據分析的流程之一。本文采用間隔影響分析(margin influence analysis， MIA)與連續投影算法(successive projections algorithm，SPA)相結合的方法對近紅外光譜數據進行降維，然后應用偏最小二乘線性判別分析算法(partial least squares linear discriminant analysis，PLS-LDA)對腹瀉性貝類毒素樣本和健康樣本的混合數據集進行分類，建立快速鑒別模型，以期為貝類毒素的快速無損鑒別提供一種新方法，為海產品污染物識別提供參考。

1 材料與方法

1.1 試驗樣品制備

樣品購買自廣東省湛江市東風海鮮市場，品種為翡翠貽貝。首先選取大小相似的貽貝，將其放置于塑料容器中進行馴化，以適應實驗環境。馴化3 d后，選擇生命力強的貽貝進行后續實驗。預先準備好2個大小為119 cm×108 cm×32 cm塑料養殖箱，分別飼養健康樣本(對照組)和毒素污染樣本(實驗組)。每個養殖箱中承裝海水80 L，鹽度為3%，溫度為26 ℃。實驗組養殖箱中加入濃度為7.3×109細胞/L的利馬原甲藻來模擬受DSP污染的海洋環境。對照組養殖箱中每天用0.5 L光合細菌喂養貽貝樣本。

在實驗過程中，為了保持貽貝良好的生理狀態，利用氧氣泵對養殖箱中的海水不斷充氣，海水每24 h更換1次，以保持貽貝的生活環境清潔。將貽貝樣本在養殖箱中連續喂養6 d，使得DSP在貽貝樣本中積累。剔除死的貽貝，共收集240個樣本(每組120個樣本)進行光譜采集。

1.2 光譜采集與預處理

貽貝從海水中取出后，利用近紅外光譜系統(圖1)獲得貽貝的光譜信息，該系統由近紅外光譜儀、鹵素光源、光纖、計算機和可調位移平臺組成。近紅外光譜儀的型號為SW2520-050-NIRA(中國臺灣OTO光電子有限公司)。采集貽貝的反射光譜范圍在950 nm到1 700 nm之間，包括114個波段。光譜采集前，預先進行黑白校正以降低噪聲[18]。

圖1 近紅外光譜系統Fig.1 Near-infrared spectroscopy system

光譜采集后，將實驗組DSP污染貽貝的樣本內殼肉取出并冷凍，用于檢測DSP含量，研究中采用GB 5009.212—2016《食品安全國家標準 貝類中腹瀉性貝類毒素的測定》中的LC-MS/MS法進行檢測。

采集到的原始光譜數據使用Savitzky-Golay卷積平滑求導結合標準正態變量變換用以消除光譜數據的噪聲以及基線平移。對于表征樣本是否為健康樣本矩陣采用平均中值法預處理，去除變量中沒有信息含量的均值部分。

1.3 MIA數據降維方法

MIA是一種基于模型集群分析思想，并以支持向量機(Support vector machines， SVM)的內在工作機制為基礎的波段選擇算法[19]。MIA方法源于間隔是反映SVM模型預測能力的一個重要指標，SVM模型的間隔越大，其結構風險就越小，其模型的泛化性能越好。也就是說，能夠增加SVM模型間隔的變量是有信息變量，反之則是無信息變量甚至是干擾變量。MIA算法原理如圖2所示。

圖2 MIA方法的主要原理Fig.2 Main principle of MIA method

假設有光譜訓練樣本集T={xi,yi|i=1，2，…，n}∈(X×Y)n，其中，n是樣本數量；xi是訓練樣本集X內的第i個樣本(p維向量)，xi∈X=Rn；yi是訓練樣本集Y內的第i個樣本的分類標簽，yi∈Y={-1，1}。

(1)基于變量空間的蒙特卡洛抽樣。圖2中行為采樣次數用N表示(一般較大)，列為變量空間大小用p表示。采樣過程如下：預先定義待采樣變量的數量，用Q表示；每次采用無放回地從p個變量中隨機選取Q個變量，得到n×Q的子數據集。重復N次，得到N個子數據集。圖2每一行表示一次抽樣即在變量空間隨機選擇Q個變量(黑色方塊)。

(2)建立SVM子模型。給定一個懲罰因子C，對于(1)中每個隨機抽樣的子數據集，建立SVM分類模型。本研究中采用交叉驗證來選擇C。由此計算出其對應的N個間隔值，記為Mi，i=1，2，…，N。由圖2可知每個子模型建模所用變量及其間隔。

(3)對SVM間隔進行統計分析。對于波長變量i，將N個子模型分為A類(包含波長變量i的模型)和B類(不包含波長變量i的模型)兩類。基于這兩類間隔的數據，可以計算得到波長變量i對應的2個分布，兩類的間隔均值分別記為Mi,A和Mi,B，則波長變量i間隔分布均值的差可以表示為公式：

Dmeani=Mi,A-Mi,B

若Dmeani>0，表明SVM模型中加入波長變量i可以提高模型性能，反之亦然。對各波長變量的間隔分布，進行非參數Mann-Whitney U檢驗[20]，計算差異顯著性的P值。選擇Dmeani>0且P<0.05的波長變量作為影響波段。MIA的最終目標是選擇能夠顯著增加間隔的波長變量。

為了降低模型的復雜性，提高模型的性能，采用SPA進一步減少最終模型中包含的波段數量。

1.4 參數調優

MIA方法中有3個調優參數，分別是N、C和Q。為了考察這3個參數對模型的影響，在一系列取值下進行研究，并對結果進行了分析，以期得到一組最優化的參數。對于蒙特卡洛抽樣次數，已有研究表明，N越大，得到的結果越好，但計算成本較高。考慮到計算成本和再現性，本工作中選擇N=10 000。SVM的懲罰因子C用于控制對分類錯誤的懲罰力度，取值越大表示越不能容忍錯分現象，但取值過大，容易出現過擬合。因此，選擇合適的懲罰因子C是至關重要的。Q是每次采樣劃分子數據集的變量個數，本研究通過檢驗MIA識別的信息變量的再現性和相應的預測誤差來確定Q的最佳取值。本研究利用交叉驗證進行C和Q的優化。

為了選擇合適的C和Q，本研究將貽貝的近紅外數據集隨機劃分為訓練集(96個健康樣本和96個DSP樣本)和測試集(24個健康樣本和24個DSP樣本)20次。樣本的近紅外光譜數據在訓練集和測試集之間沒有重復。本研究中，懲罰因子C在1～12取值，步長為1；子數據集變量個數Q的取值范圍為10～50，步長為10。圖3給出了對DSP貽貝和健康貽貝樣本混合數據集執行MIA算法后得到的預測誤差隨參數C與Q的變化趨勢。由圖3可知，隨著C值的增加，預測誤差的變化會因Q的取值不同而不同，說明C和Q的取值均會影響MIA算法的結果，應選擇預測誤差最小結果對應的參數作為最優參數值。由圖3-a亦可知，在采樣劃分子數據集的變量個數Q為30時，對于不同的C值，模型均取得較低的預測誤差(圖3-a中紅色實線表示)。當Q=30且C=6時，算法取得了最小預測分類誤差0.044 5(圖3中五角星表示)。優化參數時，當對應C和Q的步長分別取值更小時，參數優化效果不明顯，對模型結果影響不大。故此研究中，所選最優參數N為10 000，C為6，Q為30。

1.5 分類模型及評價標準

采用PLS-LDA，SVM和隨機森林(random forest，RF)3種方法分別建立分類模型，比較不同模型對DSP貽貝、健康貽貝的識別準確率。同時，為了減少因樣本集劃分而導致的結果差異，提高模型穩定性，采用20次5折交叉驗證法進行訓練集、測試集的劃分和模型訓練，并將其平均值作為最終分類結果。

PLS-LDA是基于PLS算法改編的用來建立自變量與響應變量之間映射關系的一種監督分類方法。SVM是一種利用超平面對樣本進行分類的方法，本研究中SVM使用的RBF作為核函數。RF是一種基于裝袋和隨機子空間的決策樹集成方法，其輸出的類別是由個別決策樹輸出類別的眾樹來決定的，可以用來解決高維數據分類等問題。

分類準確率可以直觀評價模型好壞，但當使用不平衡數據集時，準確率并不能反映全面情況。本研究中引入受試者工作特征曲線(receiver operating characteristic，ROC)來衡量模型性能，它能夠準確反映模型特異性和敏感性的關系，是實驗準確性的綜合代表，一般以ROC曲線下面積(area under the curve，AUC)作為模型評價指標，其值越大代表模型的鑒別效果越好，其值最大為1，分類器的性能與AUC成正比[21]。

a-不同Q值下預測誤差隨C變化二維圖； b-不同Q值下預測誤差隨C變化三維圖圖3 懲罰因子C與采樣變量數Q對預測誤差的影響Fig.3 Effect of penalty factor C and sampling variable number Q on prediction error

2 結果與分析

2.1 健康樣本和DSP污染樣本的原始光譜分析

混合樣品集(120個DSP樣品和120個健康樣品)的原始近紅外光譜如圖4-a所示。圖中光譜曲線走向趨勢相似，說明二者的內部含有的化學成分大致相似。在1 050 nm附近有明顯的波峰，為N—H的三級倍頻[22]，可能與貽貝中蛋白質的近紅外吸收有關。圖4-b顯示了混合樣本集120個DSP貽貝和健康貽貝樣品平均光譜曲線，二者表現出一定的差異，在950 nm至1 080 nm波長范圍內，健康樣品的光譜反射率值高于DSP污染樣品的光譜反射率值；在1 150～1 700 nm，DSP污染貽貝反射強度高于健康樣本，在其他波長范圍內，2種樣品的平均光譜曲線幾乎重疊。光譜采集后，對實驗組樣本采用LC-MS/MS法檢測DSP含量為35 μg/kg。

貝類樣品基質復雜，內源性化合物的存在可能會改變酶的自然環境，影響其活性。當貽貝受到海洋毒素污染時，會產生蛋白質、酶和脂質等組織成分的變化，這些變化反映在光譜曲線上。這些光譜的差異為區分DSP污染樣品和健康樣品提供了可行性。

a-原始光譜圖；b-平均光譜圖圖4 樣本近紅外光譜Fig.4 Near-infrared spectra of samples

由于波長和重要背景的重疊，很難明確識別特定化學成分的波長。為了探索DSP污染和健康樣本的分布，采用主成分分析(principal component analysis，PCA)提取新變量，以發現樣品之間可能存在的聚類現象。第一個主成分(PC1)、第二個主成分(PC2)和第三個主成分(PC3)分別解釋了總方差的95.60%、2.08%和1.32%，累計方差貢獻度在99%以上。圖5所示為DSP污染貽貝和健康貽貝的PCA得分圖和載荷圖。可以注意到，圖5-a中，健康樣品比DSP污染樣品更緊密地聚集在一起。然而，DSP污染樣品和健康樣品之間仍然存在重疊。明顯的重疊意味著僅僅用簡單的分類方法不能夠區分它們。因此，需要化學計量學方法和模式識別方法對DSP污染和健康樣品進行分類。如圖5-b所示，對PC1貢獻最大的波段為1 030 nm，對PC2貢獻最大的波段為1 470 nm和1 050 nm，對PC3貢獻最大的波段為1 690 nm和1 010 nm。根據相關文獻，波段1 470 nm主要對應的是與蛋白質相關的譜帶，1 690 nm歸屬于脂質信號[23]。

在獲取近紅外光譜數據時，由于近紅外光譜系統的工作條件和樣品結構特性的變化可能會導致光譜中的基線漂移、隨機噪聲和多元散射效應。為了改進近紅外光譜數據，通常采用光譜預處理方法來減少這些影響[24]。本研究采用Savitzky-Golay卷積平滑求導(多項式階數為2，滑動窗口寬度為13，求導階數為1)結合標準正態變量變換對原始光譜進行預處理，以消除光譜數據的噪聲及基線漂移，可以增強鑒別模型的魯棒性。

a-得分圖；b-載荷圖圖5 DSP樣本與健康樣本數據集PCA得分散點圖及載荷圖Fig.5 PCA score scatter diagram and load diagram of data set for DSP sample and health samples

2.2 選取特征波長

采用MIA方法進行波長變量篩選，對于預處理后的DSP貽貝與健康貽貝混合集光譜數據，分別選出一個有信息的波長與一個無信息的波長，其對應的間隔的分布如圖6所示。圖6-a和圖6-b所對應的分布分別為強信息變量(波段1 473.77 nm，P=6.733 2×10-12)和無信息變量(波段989.78 nm，P=1.730 2)的分布。其中，峰“1”表示含有相應波長變量的SVM模型的間隔分布，峰“0”表示不包含此波長變量的SVM模型的間隔分布。很明顯，圖6-a中包含波長1 473.77 nm與不包含波長1 473.77 nm變量建立模型的間隔分布存在明顯差異。包含波長1 473.77 nm 變量建立的模型間隔分布發生了右移，其間隔平均值大于不包含該波長變量的，意味著該變量蘊藏著對DSP貽貝和健康貽貝進行判別的有用信息。這意味著將波長1 473.77 nm變量加入SVM模型中，可以提高模型的泛化性能。相比之下，圖6-b中，包含波長989.78 nm與不包含波長989.78 nm變量所建立的SVM模型的間隔分布沒有明顯差異，且包含波長989.78 nm變量的SVM模型的間隔分布位于左側，說明此變量對DSP貽貝和健康貽貝鑒別是無用的信息變量，將其加入到SVM模型中，不但不會增加，反而會減小模型的間隔，降低模型的預測能力，需要從模型中剔除。

為了建立DSP貽貝和健康貽貝混合樣本數據集預測的分類模型，應該確定一個波段子集。經過MIA方法分析后，對于950～1 700 nm波段范圍內的114個波長變量，各波長變量Dmean值的分布如圖7所示。根據算法原理，Dmean<0為無信息變量，由圖7可知，位于0線以下的均是此類變量，有25個變量被剔除。利用非參數Mann-Whitney U和Holm-Bonferroni方法進行統計檢驗，對剩余的89個變量按計算獲得的每個波長變量P值進行排序(Dmean>0)，然后按順序選出一定數量的波長變量建模，并采用交叉驗證對預測性能進行評價。由于MIA中使用了蒙特卡洛策略，所建立的SVM分類器的預測誤差不能準確地再現，因此，通過運行20次MIA程序來研究分類誤差隨變量數量的變化。光譜數據集的平均5折交叉驗證誤差及標準差如圖8所示。由圖8可知，當模型輸入變量少于35個時，平均5折交叉驗證誤差隨變量數的增加都是先逐漸減小，然后達到最小。超過35個變量時，誤差隨變量數的增加逐漸增大。最終確定35個特征變量能夠顯著增加SVM間隔，即對DSP貽貝和健康貽貝鑒別有用的變量。

a-有信息波長1 473.77nm；b-無信息波長989.78nm圖6 兩類信息變量SVM模型的間隔分布Fig.6 Margin distribution of SVM model with two types of information variables

為了進一步減少波段數量，簡化最終模型，對MIA選擇的35個特征變量采用SPA獲取具有最小冗余信息波段，并依據模型的交叉驗證均方根誤差的最小值來確定特征波長變量個數，最終篩選出來的特征波長共有7個，分別為1 029.56、1 473.77、1 480.40、1 520.18、1 526.81、1 672.67、1 692.56 nm。最后選擇的特征波段如圖9所示。

圖7 運行MIA方法后每個波長變量的Dmean值Fig.7 Dmean value of each wavelength variable after running the MIA method

圖8 選擇不同變量數目時的5折交叉驗證誤差及標準差Fig.8 Error and standard deviation of 5-fold cross-validation when choosing different numbers of variables

圖9 MIA-SPA方法選擇的波長變量的分布Fig.9 Distribution of wavelength variables selected by the MIA-SPA method

2.3 判別模型及評價

為了構造最優鑒別模型，本項目對MIA-SPA特征波長提取方法得到的特征變量，分別應用PLS-LDA、SVM和RF建立DSP樣本和健康樣本的無損鑒別模型，并與MIA、SPA和全波段光譜建立的模型進行對比分析。不同降維方法下各模型的鑒別結果如表1所示。由表1可知，PLS-LDA對DSP樣本和健康樣本均具有很高的鑒別能力，其中DSP樣本的分類準確率達到了100%；RF對DSP樣本和健康樣本的分類準確率偏低，不能達到分類要求；而SVM對DSP樣本的分類準確率較高，達到了96.67%，但是對健康樣本的分類準確率在85%以下，低于PLS-LDA的分類性能。由表1可知采用全波段建立的PLS-LDA、SVM和RF模型的分類準確率較采用降維后的變量建立的模型低，這是由于采集的樣本光譜數據中包含了對表征樣本特征相關性較差的波段，這些波段會降低模型的分類效果，而數據降維可篩選出更能表征樣本的特征波段并將無關的波段予以刪除。與MIA、SPA相比，采用MIA-SPA方法篩選的特征波長建立的模型具有較好的分類準確率，特別是MIA-SPA-PLS-LDA模型的分類效果最優，對DSP污染樣本鑒別準確率達到100%。

表1 不同降維方法下各模型的分類效果Table 1 The classification effect of each model under different dimensionality reduction methods

為了進一步評價各種分類模型的性能，分別作出了上述分類模型ROC曲線。圖10所示為不同降維方法下的3種鑒別模型的ROC曲線對比結果圖。由圖10可知，相同的降維方法得到的輸入變量所建立的3類模型中，PLS-LDA模型的AUC相對SVM與RF模型的要大，說明本研究中PLS-LDA模型具有更好的性能；由圖10-a可知，采用MIA-SPA降維后建立的PLS-LDA模型的AUC為0.989，均高于其他模型，說明MIA-SPA-PLS-LDA模型對DSP樣本和健康樣本均具有很高的鑒別能力。

2.4 討論

本研究采用MIA-SPA-PLS-LDA結合近紅外光譜技術對腹瀉性貝類毒素進行檢測。貽貝被DSP污染后體內發生極其復雜的化學轉化和酶轉化機制，將導致受污染貽貝和健康貽貝的化學成分存在差異。已有研究表明，這些差異可以被近紅外光譜捕獲[17,25]，因此，本文采用近紅外光譜檢測DSP污染是有依據支撐的。本文進一步驗證了近紅外光譜結合波段選擇方法以及判別模型，無損檢測DSP污染的貽貝是可行的。

a-MIA-SPA;b-全波段;c-MIA;d-SPA圖10 不同降維方法下各分類模型的ROC曲線Fig.10 ROC curves of each classification model under different dimensionality reduction methods

貽貝因形態及非光滑外殼可能會導致光譜變異性和光散射，因此，本文采用Savitzky-Golay卷積平滑求導結合標準正態變量變換預處理，可有效消除光譜數據的噪聲及基線漂移。特征選擇可以篩選出與DSP檢測最相關的波段信息，同時提高分類模型的性能，本文提出采用MIA-SPA方法獲取最優特征波段，據報道，在1 470 nm處發現了與蛋白質相關的譜帶，在1 690 nm發現了脂類信號，說明MIA-SPA算法選擇的波段可以反映光譜差異[23]。結果表明，從光譜中提取重要信息，可以實現對含DSP貽貝的快速檢測。

3 結論

本文采用近紅外光譜技術和化學計量學方法，對DSP污染和健康貽貝的混合樣本集進行無損鑒別研究。采用MIA-SPA方法能夠有效地對光譜數據進行降維(降維后變量為7個)，同時結合PLS-LDA方法能夠較好地對樣本進行無損鑒別(DSP分類準確率為100%)，通過ROC曲線分析，MIA-SPA-PLS-LDA方法的AUC最大，其模型鑒別效果最優。結果表明，MIA-SPA-PLS-LDA模型用于DSP污染快速無損鑒別是可行的，對DSP污染貽貝的檢測準確率達到100%，DSP含量檢出限為35 μg/kg。該研究結果可為后續各種海洋貝類毒素鑒別分析提供理論依據。