局部偏最小二乘法結合可見-近紅外光譜預測豬肉揮發性鹽基氮

王文秀，彭彥昆*，王 凡，馬 營

（1.中國農業大學工學院，國家農產品加工技術裝備研發分中心，北京 100083；2.河北農業大學食品科技學院，河北 保定 071000）

我國是豬肉生產和消費大國，2017年豬肉產量達5 340萬 t。揮發性鹽基氮（total volatile basic nitrogen，TVB-N）是豬肉在貯藏期間蛋白質的分解產物，是反映豬肉是否具有可食用性的重要指標[1]。GB 5009.228—2016《食品中揮發性鹽基氮的測定》規定了TVB-N的標準測定方法，但是該方法存在前處理繁瑣、耗時長、破壞樣品等缺點，無法滿足快速、實時、無損傷的檢測需要[2]。

近年來，可見-近紅外光譜技術與化學計量學結合，已經應用到肉品品質[3-6]、安全[7]、摻假[8-9]、種類識別[10]等研究上，并取得了令人滿意的結果。Prevolnik[11]、Prieto[12]、Weeranantanaphan[13]等綜述了光譜技術在肉品檢測中的應用，表明該技術能夠無損檢測肉品主要參數。在預測TVB-N方面，Cai Jianrong等[14]基于1 000～2 500 nm之間的近紅外光譜，建立了豬肉中TVB-N的偏最小二乘模型，預測相關系數為0.808 4。馬世榜等[15]通過特征波長篩選建立了牛肉TVB-N的預測模型，相關系數為0.925 0。上述研究基于校正集樣品建立的模型直接預測驗證集樣品，證實了光譜技術預測TVB-N的可行性。

由于近紅外模型建立過程也包括了樣品狀態、環境變化、儀器條件等因素，因此會出現建立的TVB-N模型在預測不同批次樣品時效果不佳的問題，限制了近紅外光譜技術在肉品行業的實際應用。為了提高模型的預測能力，Naes等[16]提出了局部回歸方法，根據“相似樣品產生相似輸出”的原理，依據待預測樣品的光譜特征，通過某種相似度依據，在數據集樣品中尋找部分與之相似的樣品，重新建立局部校正模型。該方法可以充分利用原始數據集樣品的信息，在土壤全氮含量反演[17]、蘋果糖度預測[18]、甜味劑濃度測定[19]等研究中已有應用，但是在豬肉TVB-N檢測方面鮮見報道。此外，在建立局部模型時，相似樣品選擇的依據和數量是至關重要的因素，一些學者對此展開了研究。Dambergs等[20]以待測樣品與數據庫樣品光譜的相關系數作為依據選擇局部建模樣品，鄢悅等[21]通過光譜信息散度選擇樣品建立局部校正模型，張紅光等[22]利用凈信號分析結合歐式距離構建局部模型。這些研究證實了局部建模策略的優勢，但是多從單個角度評價不同樣品的相似程度，判斷能力有待提高。相似度準則的選擇仍然是目前需要進一步研究和解決的問題。

為了解決豬肉TVB-N預測模型在應用中面臨的上述具體問題，本實驗針對兩個批次的豬肉實驗樣品，基于350～2 500 nm波段的可見-近紅外反射光譜矩陣，提出了基于距離、信息測度和投影的相似性度量方法，通過構建相似度函數和相似度因子，建立TVB-N的局部偏最小二乘模型，提高模型的預測能力。

1 材料與方法

1.1 材料與試劑

豬肉背最長肌部位購買于北京順鑫農業股份有限公司鵬程食品分公司、雙匯冷鮮肉專賣店以及北京二商大紅門肉類食品有限公司。

0.01 mol/L標準鹽酸滴定溶液 廈門海標科技有限公司；氧化鎂、硼酸、95%乙醇（均為分析純） 西隴化工股份有限公司；指示劑甲基紅和溴甲酚綠 廣東光華科技股份有限公司。

1.2 儀器與設備

KDY-9820半自動凱氏定氮儀 北京瑞邦興業科技有限公司；MYP11-2A磁力攪拌器 上海梅穎浦儀器儀表制造有限公司；5 mL移液槍 德國Eppendorf公司；可見-近紅外反射光譜采集系統如圖1所示，包括可見-短波近紅外光譜儀AvaSpec-2048x14（荷蘭Avantes公司）、長波近紅外光譜儀AvaSpec-NIR256-2.5（荷蘭Avantes公司）、Y型光纖、環形光導、鹵鎢燈光源等。兩臺光譜儀的波長范圍分別為350～1 100 nm和1 000～2 500 nm，以下分別簡稱為前波段光譜和后波段光譜。

圖1 可見-近紅外光譜采集系統示意圖Fig. 1 Schematic diagram of the visible and near-infrared spectral acquisition system

1.3 方法

1.3.1 樣品處理及實驗設置

本研究包括兩個批次的實驗。在第1批實驗中，選取屠宰后經過排酸后熟的豬肉背最長肌部位作為樣品，去除外層筋膜后，將樣品分割為大小約8 cm×5 cm×2.5 cm（長×寬×高）的肉塊，放置在自封袋中包好并編號，然后用蓄冷箱運送至中國農業大學工學院無損檢測實驗室，無積壓放置在4 ℃冰箱中保存。共獲得有效實驗樣品108 個，包括鵬程肉40 個，雙匯肉38 個，大紅門肉30 個。為了增大TVB-N標準參考值的范圍，每天從冰箱中隨機取出9 個樣品用于光譜的采集和TVB-N含量的測定，實驗持續12 d。

在第2批實驗中，同樣以排酸后的豬肉背最長肌部位為樣品，包括鵬程肉15 個，雙匯肉15 個，大紅門肉10 個。樣品的前處理方法與第1批實驗保持一致，分割包裝后置于4 ℃冰箱中保存待用。為了獲得不同新鮮程度的樣品，每天隨機取出4 個樣品用于實驗，共持續10 d。

上述兩個批次的實驗在不同的時間進行，第1批次樣品用來建立包含多個品種豬肉的校正模型并形成建?；A數據集，第2批次樣品用來驗證所建模型對不同批次樣品的適用性，同時驗證提出的局部偏最小二乘建模方法的可靠性。

1.3.2 光譜采集

實驗開始前，將樣品從冰箱取出并在室溫下靜置，同時打開儀器進行預熱。實驗時，首先調節環形光導與標準硫酸鋇白板的距離至形成光強均勻無暗影的光斑，確定距離為4 cm。然后，依次采集參比光譜和暗背景光譜進行校準，并在樣品表面選取5 個不同位置采集反射光譜信息，平均后作為該樣品的最終光譜。

1.3.3 TVB-N理化值測定

參照GB 5009.228—2016方法對TVB-N的標準理化值進行測定。將豬肉絞碎后準確稱?。?0±0.1）g置于錐形瓶中，加入100 mL蒸餾水攪拌30 min并過濾。然后，準確量取10 mL濾液和10 mL質量濃度為10 g/L的氧化鎂溶液，加入到消化管中蒸餾5 min，硼酸吸收液用0.01 mol/L的鹽酸標準液進行滴定，根據消耗的鹽酸體積計算TVB-N含量。每個樣品做3 個平行實驗，取平均值作為該樣品最終TVB-N值。

1.4 數據分析

1.4.1 雙波段數據融合方法

為了更加充分的利用兩個光譜儀采集的信息，需要對雙波段光譜數據進行融合，以得到連續覆蓋整個可見光及近紅外區域的光譜曲線。融合方法如下：首先，截掉后波段光譜在1 000～1 074 nm和2 279～2 500 nm范圍內噪音較大的光譜數據；然后，保持前波段光譜在350～1 074 nm范圍及后波段光譜在1 369～2 279 nm范圍的數據不變，利用式（1）對后波段光譜在1 074～1 369 nm范圍的數據進行拋物線擬合校正，得到校正后的光譜數據；最后，將350～1 074 nm范圍的原始光譜數據、1 074～1 369 nm范圍的校正光譜數據以及1 369～2 279 nm范圍的原始光譜數據組合形成融合后的光譜矩陣。

式中：Rbefore（λ）和Rafter（λ）分別為校正前后波長λ位置處的反射率/%；a、b、c分別為拋物線擬合方程的3 個參數。由于1 369 nm為拋物線方程的頂點，因此b和c分別是1 369和1，且為固定值。參數a通過將待融合的前后兩個波段光譜在1 074 nm波長處反射率比值代入到式（1）求得。

1.4.2 光譜預處理及建模方法

基于第1批次的實驗樣品建立TVB-N的預測模型，并利用該模型直接預測第2批次樣品，驗證模型在不同批次樣品之間的適用性。建模前，采用標準正態變量變換（standard normal variate transformation，SNV）對光譜進行預處理，消除散射對光譜的影響[23]。由于兩臺光譜儀獲取的數據間隔不同（分別為0.59 nm和6.80 nm），導致前后兩個波段范圍的光譜數據個數不同，前波段數據明顯多于后波段數據。為了在光譜預處理和建模過程中，兩個波段的光譜信息能具有相同的權重，采用“cublic”插值法對融合后光譜以2 nm為間隔進行數據重排。建模采用偏最小二乘（partial least square，PLS）法分析，利用校正集相關系數（correlation coefficient in the calibration set，Rc）、預測集相關系數（correlation coefficient in the prediction set，Rp）、校正集標準分析誤差（standard error of calibration，SEC）和驗證集標準分析誤差（standard error of prediction，SEP）對模型進行評價。

1.4.3 局部偏最小二乘回歸模型建立

以第1批次樣品為建模基礎數據集，建立第2批次樣品的局部偏最小二乘回歸（locally partial least square，LPLS）模型。建模的步驟包括：首先，確定相似度函數來評價不同樣品之間的相似程度；然后，根據相似度函數，在建?；A數據集樣品中尋找與待測樣品相似的樣品，確定局部回歸建模鄰域窗口；最后，利用局部空間內樣品建立待測樣品的LPLS模型。各個步驟具體實施過程如下。

1.4.3.1 相似度函數選擇

為更加準確全面地衡量光譜之間的相似度，提出基于距離、信息測度和投影的相似性度量方法，分別采用歐式距離、光譜信息散度（spectra information divergence，SID）和光譜角（spectra angle metric，SAM）計算不同樣品的相似程度。目前，SID和SAM在高光譜分析中有所應用，但是在近紅外分析中應用較少[24-27]。將SID和SAM結合使用，可以從光譜形狀和光譜信息熵角度判斷相似程度，其計算公式如式（2）所示。同時定義了相似度函數S，對歐氏距離和光譜信息散度-光譜角（spectral information divergence-spectral angle，SID-SAM）進行加權求和，計算公式如（3）所示：

式中：x和y為待分析的兩條光譜；d（x, y）為光譜之間的歐式距離；λ為權重因子，取值范圍為0～1，當λ為0或1時，表示只計算光譜形狀信息或距離。不同樣品之間的相似程度越高，d（x, y）和SID-SAM值越小。在計算歐式距離時，將原始光譜矩陣均值中心化后利用主成分分析進行降維，利用得分矩陣代替原始矩陣計算歐式距離。

1.4.3.2 鄰域確定

利用式（2）、（3）計算出待測樣品與建?；A集中各樣品的相似度后，將樣品按相似性從大到小排序，依次選擇樣品組成鄰域窗口h。為了提高確定窗口h的效率，定義相似度因子SM，其計算公式如式（4）所示，其中Mmax為建立LPLS模型的最大樣品數，為經驗值，M為實際建模的樣品數，介于1和Mmax之間。相似度因子SM可理解為窗口h內樣品的相似度總和在最大LPLS建模樣品相似度總和中的比重，可通過代價函數（式（5））計算確定，代價函數最小時的SM為LPLS建模時的SM。

式中：Si為待測樣品與排序后建模集第i個樣品的相似度；X為建立LPLS模型的樣品組成的子集；MSE為對外部驗證樣品的預測誤差；y（k）為第k個樣品的標準參考值；f（X, k）為LPLS模型對該樣品的預測值。

1.4.3.3 LPLS模型建立

結合實驗數據對參數λ和SM尋優，確定最佳取值后，可確定每個待測樣品LPLS模型的建模樣品。然后，利用競爭性自適應重加權采樣（competitive adaptive reweighted sampling，CARS）算法實時優選特征波長，并基于特征波長下的光譜建立LPLS模型。最佳潛變量數通過留一法交叉驗證確定，分析流程如圖2所示，分析過程在Matlab2012a中進行。

圖2 LPLS模型分析流程圖Fig. 2 Flow chart of local partial least squares analysis

2 結果與分析

2.1 原始光譜和TVB-N參考值分析

圖3 實驗樣品的原始光譜Fig. 3 Original spectra of experimental samples

獲得兩個波段的光譜數據后，首先利用拋物線擬合法對雙波段數據進行融合，對某個樣品連接前后的光譜對比圖如圖3A所示。從融合前的雙波段光譜（前波段為350～1 100 nm，后波段為1 074～2 279 nm）可以看出，在1 074～1 100 nm波段范圍內，由于儀器響應不同，兩臺光譜儀獲取的光譜數據并不一致。而利用拋物線擬合法進行校正后，兩段光譜可以有效的融合為一條連續無陡變的光譜。融合后的兩個批次樣品的原始光譜圖如圖3B和3C所示，其中545 nm為肌紅蛋白的吸收峰，980 nm為O—H鍵的吸收峰，1 280 nm與N—H鍵振動有關，2 200 nm與C—H鍵振動有關[28-30]。對比圖3B和C，可以看到不同批次實驗樣品的光譜趨勢基本一致，僅在光譜強度上有所差異。針對一個批次樣品的光譜來看，在630～925 nm波段范圍之間，不同樣品的光譜存在差異，這可能與樣品品種和來源不同有關。

2 個批次實驗樣品的標準參考值統計信息如表1所示，TVB-N含量范圍分別為7.49～44.17 mg/100 g和8.17～39.54 mg/100 g。根據GB 2707—2016《鮮（凍）畜、禽產品》的規定，TVB-N含量小于15 mg/100 g時樣品為新鮮肉，可見2 次實驗均涵蓋了新鮮和腐敗的樣品，且第1次實驗的理化值范圍大于第2次實驗。

表1 2 個批次實驗樣品的TVB-N含量統計信息Table 1 Descriptive statistics of TVB-N contents in the two sets of samples mg/100 g

2.2 第1批次樣品建模結果分析

對第1批實驗樣品的雙波段光譜進行數據融合且重排后，利用SNV對光譜進行預處理。然后以3∶1的比例對108 個樣品進行分組，則有80 個樣品組成校正集建立模型，28 個樣品組成驗證集對模型進行驗證。建立的PLS模型預測性能較好，Rc和Rp分別為0.952 0和0.940 8，SEC和SEP分別為2.286 5 mg/100 g和2.442 4 mg/100 g。

2.3 直接預測第2批次樣品結果分析

利用2.2節中建立的PLS模型直接預測第2批實驗樣品，驗證所建立的模型對不同批次樣品的適用性。在進行模型驗證時，同樣需要對第2批樣品的光譜進行雙波段數據融合，并以2 nm為間隔進行數據重排。進行SNV預處理后，將光譜矩陣與2.2節建立的模型系數矩陣相乘可以得到預測值，第2批40 個樣品的真實值與預測值散點圖如圖4所示。

從圖4可以直觀地看出，一些樣品的預測效果不理想，真實值和預測值的偏差較大，總體相關系數R為0.845 6，SEP為4.581 0 mg/100 g。這表明不同批次實驗樣品之間存在差異，利用建立的模型直接預測不同樣品時容易產生較大的誤差，模型的適用性需要進一步增強。

圖4 直接預測第2批次樣品時真實值和預測值散點圖Fig. 4 Comparison between actual values and predicted values of the second batch of samples using the established model

2.4 LPLS模型預測第2批次樣品結果分析

由于直接利用第1批樣品建立的模型預測第2批樣品時，預測結果有待提高，因此進一步考察LPLS模型的預測效果。以第1批樣品為建模基礎數據集，計算第2批樣品與其中每個樣品的歐式距離和SID-SAM值。從圖3可以看出，不同來源的樣品光譜存在差異，為了消除這種差異對求取歐式距離和SID-SAM值的影響，在計算之前對兩批樣品進行極差歸一化預處理，使光譜處在一個相同的數據范圍內。具體方法為：將樣品的光譜數據看作是p元行向量（p為變量數），將該行的每一個數據與最小值求差，然后除以該行數據中最大值和最小值的差。經過處理后的兩批實驗樣品光譜曲線如圖5所示，紅色和藍色分別為第1批和第2批實驗樣品的光譜。可見數據都映射到0～1之間，各個變量和平均值分布更加均衡，光譜間無明顯差異，有利于后續的計算和分析。

圖5 極差歸一化處理后兩個批次樣品的光譜Fig. 5 Normalization pretreated spectra of the two batches of samples

在建立LPLS預測模型，不同的λ和SM取值會獲得不同的鄰域，進而得到不同的預測結果，因此需要對這兩個參數進行尋優以確定最佳取值。在計算相似度時，由于距離對光譜差異的影響大于形狀和角度，因此賦予距離更大的權重，設定λ的取值分別為0.5、0.6、0.7、0.8、0.9和1.0。SM的取值直接影響建立局部模型的窗口大小，取值太大，會導致建立LPLS模型的樣品數太多，影響模型的準確性及實時建模的速度和效率；若取值太小，會導致選擇的樣品數過少，不能準確反映光譜特征與標準參考值之間的關系。在計算出待測樣品與建模基礎數據集中每個樣品的相似度后，將樣品按相似性從大到小排序。根據經驗值設定Mmax為40，即建立LPLS模型的最大樣品數量為40，計算出前40 個樣品的相似度總和。然后設定SM分別為0.50、0.55、0.60、0.65、0.70、0.75、0.80、0.85、0.90、0.95和1.00，共11 個取值。結合λ的6 個取值，可以得到66 種λ和SM組合，因此對于第2批實驗中的每一個樣品，可以建立66 個LPLS預測模型，進而得到66 個預測誤差。對每種組合下第2批40 個樣品的預測誤差求平均，比較66 種組合下的平均誤差值，通過最小化該誤差值，確定λ和SM的最佳組合。經過比較分析，在歐式距離和SID-SAM的權重分別為0.8和0.2，相似度因子SM為0.8時，第2批40 個樣品的平均預測誤差最小，因此確定0.8為λ和SM兩個參數的最優值。

根據圖2的數據分析流程，在λ和SM值確定后，分別得到第2批40 個樣品的建模鄰域窗口，利用CARS算法實時優選特征變量，并建立LPLS預測模型，真實值與預測值的散點圖如圖6所示。與圖5中直接預測的結果相比，模型效果有了明顯的改善，預測相關系數R上升至0.948 1，SEP下降至2.650 8 mg/100 g。利用顯著性檢驗方法對真實值和預測值進行統計學分析可知，P＞0.05，表明兩組數據無顯著性差異。這也說明利用LPLS方法建模，能有效提高對外部驗證樣品的預測能力。

圖6 采用局部建模方法對第2批次樣品預測結果Fig. 6 Prediction results for the second batch of samples using LPLS

本研究結合歐式距離和SID-SAM，可以更加綜合全面地評價不同光譜之間的相似性，更有利于找尋與待測樣品光譜特征相似的建模樣品。通過優化SM值，對每一個待測樣品，均能從建模基礎數據集樣品中動態選擇最佳建模鄰域窗口。相比在建立LPLS模型時，根據待測樣品的光譜特征，在基礎數據中以待測樣品為中心搜索固定距離空間內的樣品，或搜尋與之相似性較高的固定個數的樣品，該方法具有更高的靈活性，可根據樣品特性實時優選出最佳建模用樣品子集。在整個分析過程中，確立鄰域空間和建立模型都是根據樣品“需要”判斷后自動進行，具有較強的自適應學習能力。與利用模型直接預測不同批次樣品的結果相比，LPLS方法具有更高的靈活性和適應性，尤其對于樣品差異較大引起非線性問題時，LPLS方法比全局PLS方法對樣品的預測能力更佳。

3 結 論

針對光譜法檢測豬肉TVB-N時，建立的校正模型對不同批次實驗樣品預測效果不佳的問題，提出了局部偏最小二乘法預測TVB-N含量的方法。計算外部驗證樣品與建?；A數據集中所有樣品的歐式距離和SID-SAM，以權重分別為0.8和0.2對二者進行加權求和，計算出相似度函數。以相似度因子為0.8選擇建模鄰域，對外部驗證的每一個樣品構建LPLS預測模型。與利用校正模型直接預測外部驗證樣品時的結果相比，LPLS模型具有更大的預測能力，相關系數R從0.845 6上升至0.948 1，SEP從4.581 0 mg/100 g下降至2.650 8 mg/100 g。在今后的研究中，收集更多品種和來源的樣品組成建模基礎數據集，結合本研究的方法，可以實現對更多樣品的無損預測。此外，也可將該方法應用于利用可見-近紅外光譜法預測豬肉其他參數的研究中。