近紅外特征光譜的羊肉TVB-N濃度預測模型

張 旭，白雪冰，汪學沛，李新武，李志剛，張小栓，4*

1.中國農業大學信息與電氣工程學院，北京 100083 2.中國農業大學工學院，北京 100083 3.石河子大學信息科學與技術學院，新疆 石河子 832003 4.食品質量與安全北京實驗室，北京 100083

引 言

中國的羊肉產量及羊肉消費量均居世界首位。保證儲存期內羊肉的新鮮度和防范質量安全問題愈發緊要。在羊肉的貯藏過程中，在微生物和內外源酶的作用下，羊肉中的脂肪和蛋白質分解產生有毒的氨(NH3)和胺類(R-NH2)[1]，并與腐敗產生的有機酸結合生成揮發性鹽基氮(total volatile basic nitrogen，TVB-N)，因此TVB-N濃度是評估羊肉質量安全的關鍵參數。

食品安全國家標準中規定的TVB-N濃度的檢測方法包括自動凱氏定氮儀法、半微量定氮法、微量擴散法。這些化學檢測方法需要破壞樣品，操作過程復雜，耗時費力，且結果易受操作水平影響，不能滿足快速、非破壞的質量安全檢測要求。

近紅外光譜(near infrared spectroscopy，NIR)檢測技術具有分析快速、操作簡便、無破壞性等特點，在肉品質檢測領域已有大量研究。被應用于新鮮度分級[2]、摻假識別[3]、等級劃分[4]、鮮凍肉鑒別等定性分析，以及化學組成(包括膽固醇[5]、脂肪[6]、水分[7])分析、感官品質(包括肉色[8]、系水力、嫩度[8])評價等定量分析。

原始近紅外光譜數據中雖然包含與特定成分相關的有效信息，但也受到噪聲及散射等因素的干擾，這些干擾會降低光譜模型的預測性能。因此，應用適當的光譜預處理方法(散射校正、平滑處理、尺度縮放)和變量篩選方法可有效消除與被測指標無關的噪聲、散射等干擾，提高光譜與被測指標間的相關性。將競爭性自適應重加權法(competitive adaptive reweighted sampling，CARS)[9]、無信息變量消除法(uninformative variable elimination，UVE)[9]和連續投影算法(successive projections algorithm，SPA)等方法應用于簡化和優化近紅外光譜預測模型已有報道。

在回歸預測分析中，偏最小二乘(partial least squares，PLS)是被廣泛應用的線性回歸方法，當被測指標與光譜數據間存在非線性關系時，采用線性的回歸方法無法實現光譜信息的充分提取，影響模型準確性。為解決光譜數據的非線性問題，非線性算法包括支持向量機(support vector machine，SVM)和偏最小二乘支持向量機(least squares-support vector machine，LS-SVM)已被用于近紅外光譜預測模型的構建。為此，探討近紅外光譜的線性及非線性預測模型對穩定可靠的定量分析非常必要。

本研究以反映生鮮羊肉質量安全的TVB-N濃度為預測對象，采集樣品680～2 600 nm的近紅外光譜數據，經2種方法剔除離散程度大的異常樣本后，對比多種預處理方法對預測模型性能的影響，以不同變量篩選方法優選特征波長，探討線性及非線性建模方法的預測性能，建立優化的生鮮羊肉儲存期TVB-N預測模型，為實現快速無損檢測生鮮羊肉中的TVB-N濃度提供參考和技術支持。

1 實驗部分

1.1 材料與儀器

在江蘇省東臺市華東山羊市場購買當天屠宰的綿羊，取背最長肌，切除其脂肪和肌膜，整形切成3 cm×3 cm×1 cm(長×寬×高)的塊狀150個，置入無菌袋并以4 ℃冷鏈運輸車在24 h內運抵中國農業大學工學院實驗室，將樣品置于生化培養箱(LRH-250，上海一恒儀器公司)中，分組后在0，4，8和20 ℃溫度下儲藏。低溫試驗時每隔24 h各取出3個樣品進行測試，20 ℃下每隔12 h取2塊樣本進行測試，共測試11 d，共獲得121份樣本。

1.2 光譜數據獲取

采用SpectraStar 2600 XT-R型近紅外光譜儀(美國Unity Scientific公司)采集羊肉樣本光譜，掃描次數為12次，掃描范圍為680～2 600 nm，分辨率為1 nm。

1.3 揮發性鹽基氮(TVB-N)濃度測定

根據GB 5009.228—2016《食品安全國家標準食品中揮發性鹽基氮的測定》中的微量擴散法測定樣品的TVB-N濃度。

1.4 數據處理與建模方法

通過MATLAB R2018b軟件(美國Mathworks公司)完成數據處理和模型構建。

1.4.1 異常值剔除方法

利用蒙特卡洛采樣(Monte-Carlo sampling，MCS)法消除異常樣本。異常剔除閾值設置為式(1)和式(2)所示，超過閾值之一的即為異常樣本。

Mthreshold=μM+3σM

(1)

Sthreshold=μS+3σS

(2)

其中，Mthreshold為各樣本預測誤差均值的閾值，μM和σM為各樣本預測誤差均值的均值和標準差；Sthreshold為各樣本預測誤差標準差的閾值，μS和σS為各樣本預測誤差標準差的均值和標準差。

同時利用馬氏距離法(Mahalanobis distance，MD)剔除MD過大的樣本，閾值設置為

MDthreshold=μ+3δ

(3)

式(3)中，MDthreshold為各樣本MD的閾值，μ和δ為各樣本MD的均值和標準差。

1.4.2 樣本集劃分方法

剔除異常樣本之后，運用光譜-理化值共生距離(sample set partitioning based on joint x-y distance，SPXY)算法劃分出75%的樣本為校正集，其余為驗證集樣本。

1.4.3 光譜預處理方法

根據處理效果預處理方法可分為散射校正、平滑處理、尺度縮放等。散射校正包括多元散射校正(multiple scattering correction，MSC)和標準正態變換(standard normal variate，SNV)，用于消除樣品顆粒尺寸差異和分布差異對漫反射光的影響。平滑處理可有效提高信噪比，常用的有Savitzky-Golay卷積平滑(S-G smoothing，SGS)、移動平均平滑(moving average smoothing，MAS)等方法。尺度縮放包括歸一化(Normalization)、中心化(Centering)、標準化(Autoscaling)等，用來消除數據尺度差異的影響。分別采用這7種方法對樣本進行光譜預處理，以找到最佳預處理方法。

1.4.4 特征波長篩選方法

所采集羊肉光譜共1 921個波長，存在冗余和多重共線性信息，篩選特征波長取代全光譜可提高模型簡潔性和計算效率。

CARS選擇波數的方法是基于回歸系數的權重，權重值越大則代表該變量對模型建立的貢獻越大，被選取的概率越大。

UVE以輸入變量及等量的隨機噪聲建立PLS模型得到回歸系數矩陣，計算各變量的穩定性并篩選穩定性大的變量，閾值設為隨機變量穩定性最大絕對值的0.99倍。UVE選出的波長呈局部連續分布，波段之間仍存在嚴重的多重共線性問題，嶺回歸(ridge regression，RR)法是解決此類問題的有效方法，但是鮮有研究將嶺回歸用于光譜檢測中的特征波長篩選，由此在UVE法的基礎上引入嶺回歸法得到改進的無信息變量消除(improved uninformative variable elimination，IUVE)法，以進一步簡化模型和提高預測精度。

SPA以正交投影分析全部波長變量，并保留對TVB-N敏感的特征變量，使變量之間共線性達到最小，降低模型輸入量。

1.4.5 預測建模方法

PLS在主成分分析的基礎上，對光譜和理化值同時進行分解，保留對光譜貢獻大的主成分，進而構建誤差最小化的最佳線性回歸模型。

SVM模型利用核函數將低維輸入映射到高維特征空間，并在高維特征空間進行線性回歸，適于處理小樣本、非線性以及高維數等問題。LS-SVM是SVM的改進方法，進一步降低計算復雜性和提高計算速度。設定徑向基函數(radial basis function，RBF)為SVM模型及LS-SVM模型的核函數，以具備交互驗證的網格搜索(grid-search)法對SVM模型的關鍵參數(C，γ)，以及LS-SVM模型的關鍵參數(γ，σ2)進行尋優。

2 結果與討論

2.1 剔除異常樣本

圖1為121個羊肉樣本的原始近紅外漫反射吸光度光譜。由于水分子中O—H鍵伸縮振動的二級倍頻和一級倍頻吸收，在980，1 440和1 940 nm附近呈現出吸收峰，在1 200和2 400 nm附近則是與C—H鍵拉伸和伸縮振動相關的波峰[8]。樣本的光譜曲線趨勢相似，但不同波段的上升和下降的趨勢不同，說明其內部化學成分存在差異。

圖1 羊肉樣本近紅外光譜圖

采用MCS法識別異常樣本，主成分個數設置為11，預處理方法為Centering，抽樣次數1000次。異常樣本檢測結果見圖2，得到閾值Mthreshold為6.4，Sthreshold為2.198，第47，89，94，103和121個樣本被判定為異常樣本。

圖2 蒙特卡洛采樣法異常值檢測結果

圖3是121個光譜樣本到平均光譜的MD分布圖，由MD的均值和標準偏差得到閾值為5.596，超出此閾值的異常樣本為47號，89號，94號，103號，與MCS法的檢測結果重合度較高。

圖3 羊肉樣本的馬氏距離分布

將剔除異常值后的116個樣本，采用SPXY算法劃分校正集和驗證集樣本。羊肉TVB-N濃度的統計分析結果如表1所示。

表1 羊肉樣品TVB-N濃度的統計結果

2.2 原始光譜及光譜預處理比較

分別使用原始光譜及經過7種方法預處理的光譜建立全波段的PLS預測模型，建模結果見表2。與原始光譜所建立的PLS模型相比，由SNV，MSC，Autoscaling處理的光譜數據建立的模型性能均下降，經Normalization，Centering，MAS處理的光譜數據所建立模型的性能沒有明顯改善，而以SGS預處理的數據建模效果最好，確定SGS為最優預處理方法。

表2 不同預處理方式的PLS預測模型比較

2.3 特征波長的選取

2.3.1 應用CARS篩選特征波長

采用CARS提取特征波長，設置蒙特卡洛采樣次數為50，采用7折交叉驗證計算。隨著采樣次數增加，圖4(a)曲線呈指數衰減，在運行次數1～5次，變量選擇個數曲線快速下降，對應粗選過程，之后進入緩慢遞減的細選過程。圖4(b)為交互驗證均方根誤差的變化趨勢圖，在運行次數1～36，交互驗證均方根誤差緩慢波動降低，隨后逐漸升高。從圖4(c)回歸系數曲線中的“*”標出了交互驗證誤差的最低點，在采樣次數為36次時，達到最小值2.206。此時變量篩選個數為14個，分別為720，725，823，834，925，1 162，1 230，1 278，1 441，1 473，1 867，1 981，2 484和2 554 nm。

圖4 基于CARS的變量選擇過程

2.3.2 應用IUVE篩選特征波長

運行IUVE算法計算1 921個光譜波長和等量隨機噪聲的穩定性，閾值設為±11.76，將超過閾值的輸入變量進行嶺回歸分析，以嶺跡法確定嶺回歸參數k值為0.2，變量篩選的原則參考文獻[10]。嶺回歸分析的結果如圖5所示，可以看出各回歸系數的嶺估計在k=0.2時基本穩定，根據回歸系數隨k值的變化趨勢結合誤差選擇變量。IUVE與未改進的UVE的波長選擇結果如圖6所示，藍色的曲線為未改進的UVE法計算出的各變量穩定性曲線，超過閾值的波長達到703個，占總波長的36.60%，洋紅色豎線為結合嶺回歸分析的IUVE法最終選擇出的144個有效變量，占總波長的比值降至7.50%，因此改進后的UVE可有效地消除各波長變量間的共線性。

圖5 嶺回歸分析嶺跡圖

圖6 基于IUVE的特征波長篩選

2.3.3 應用SPA篩選特征波長

在變量個數1～28的范圍內優選波長，SPA算法以RMSE的大小為依據確定特征波長數量。隨著特征波長數量的增加，RMSE的變化過程如圖7所示。當波長數量由1增加到7時，RMSE迅速下降，表明此類波長變量為與羊肉TVB-N相關的重要波長變量。當波長數量由7個增加到15個時RMSE呈波動式下降，此類波長為有用信息變量。隨著波長數量由15個繼續增加，RMSE繼續緩慢下降。因此以15個特征波長作為輸入的特征變量。圖8為15個特征波長在全光譜中的位置分布，分別為680，798，1 067，1 266，1 497，1 498，1 901，1 920，1 936，2 009，2 263，2 386，2 391，2 575和2 583 nm。

圖7 基于SPA的波長個數選擇

圖8 基于SPA的特征波長選擇

2.4 模型建立

2.4.1 PLSR模型

圖9 CARS-PLS模型對羊肉TVB-N濃度的預測結果

表3 不同波長提取方法的PLS預測模型比較

2.4.2 SVM模型

圖10 CARS-SVM模型對羊肉TVB-N濃度的預測結果

2.4.3 LS-SVM模型

表4 不同波長提取方法的SVM預測模型比較

圖11 CARS-LS-SVM模型對羊肉TVB-N濃度的預測結果

表5 不同波長提取方法的LS-SVM預測模型比較

3 結 論

利用近紅外光譜對羊肉TVB-N濃度進行預測，主要結論如下：

(1)以MCS法和MD法剔除了羊肉的光譜數據的5個異常值，且2種方法的檢測結果重合度和可信度較高。

(2)以原始光譜和6種不同方法預處理的光譜建立了PLS預測模型，SGS處理的光譜建模效果最好，平滑處理、尺度縮放方法的建模效果整體上好于散射校正方法。

(3)利用CARS，UVE，IUVE，SPA提取特征光譜得到的波長個數分別為14，703，144，15，占全光譜1921個波長的0.73%，36.60%，7.50%，0.78%。對比全光譜和各方法提取的特征波長所建立的預測模型，CARS提取的波長建立的模型性能最優。對比UVE法，IUVE可消除波長間的共線性和提高模型性能。

(4)對提取的特征波長建立了儲存期生鮮羊肉TVB-N的PLS，SVM和LS-SVM預測模型，最好的校正集預測結果由SVM模型取得，最好的驗證集預測效果由LS-SVM模型得到，這兩種方法的建模效果與模型參數和建模樣本密切相關。