基于高光譜成像技術的羊肉新鮮度預測

張 玨 田海清 王 軻 張麗娜 于 洋

(1.內蒙古農業大學機電工程學院，呼和浩特 010018； 2.內蒙古師范大學物理與電子信息學院，呼和浩特 010022)

隨著生活水平的不斷提高，人類膳食結構逐漸發生變化，消費者對畜禽生鮮肉的需求量也在迅速增加。羊肉味道鮮美，肉質細膩，但在存儲、加工及物流等環節易受酶、微生物等作用發生腐敗變質，準確、快速地對生鮮羊肉進行品質檢測和安全評定尤為重要。肉品新鮮度是衡量肉品品質主要指標，新鮮度的評定主要采用感官檢測與實驗室檢測分析相結合的方法[1]。感官檢測主要通過專業人員對肉的色澤、黏度、彈性和氣味等進行評價，存在易受主觀因素影響、測量誤差大且結果不宜量化等缺陷[2]。實驗室檢測通常以揮發性鹽基氮(Total volatile basic nitrogen，TVB-N)、pH、硫代巴比妥酸值(Thiobarbituric acid reactive substances，TBARS)、顏色參數和菌落總數(Total viable counts，TVC)等為主要檢測指標，該方法檢測精度較高，但通常需要經過復雜的樣品前處理過程，存在步驟繁瑣、耗時費力等缺點，且無法滿足快速、無損檢測肉品新鮮度的要求[3]。

高光譜成像技術具有波段多、光譜分辨率高和圖譜合一等特點[4]，現已成為遙感領域的前沿技術。目前，國內外研究者應用該技術在肉品新鮮度無損檢測方面已經取得不少成果。朱榮光等[5]以全波段羊肉反射光譜信息作為輸入量，利用逐步多元線性回歸(Stepwise multiple linear regression，SMLR)、偏最小二乘回歸(Partial least squares regression，PLSR)和主成分回歸(Principal component regression，PCR)3種建模方法，建立并驗證了羊肉TVB-N含量的預測模型。張晶晶等[6]運用標準變量變換等5種光譜預處理方法變換原始光譜，利用競爭性自適應重加權算法(Competitive adaptive reweighting sampling，CARS)、連續投影算法(Successive projection algorithm，SPA)提取特征波長，建立了灘羊肉TVB-N含量的偏PLSR預測模型，并對灘羊肉貯藏時間進行判別分析。上述研究表明，應用高光譜成像技術進行肉品新鮮度檢測具有很好的應用前景。但SMLR和PLSR等常規模型較多擬合了肉品生化組分與高光譜變量之間的線性關系，而肉品腐敗過程復雜，具有明顯的時空分異和非線性特征，且因變量的變化并非只受單個自變量的影響，而是多個變量綜合影響的結果，因此該類方法對肉品品質變化的解釋力較弱，建立模型的預測精度較低。

隨機森林回歸(Random forest regression，RFR)算法有較強的數據擬合能力，能有效地分析非線性、共線性和具有交互作用的數據[7]，解決利用多變量預測因變量的問題，在土壤養分診斷[8]、植被重金屬含量分析[9]和空氣質量檢測[10]等方面有較好的應用效果，但該算法在肉品TVB-N含量預測方面的應用鮮有報道。另外，RFR建模參數通常通過試數法確定或根據其應用對象的經驗值設定[11-12]，而依據上述方法建立的RFR模型算法很難取得最優解，因此，控制參數實現智能尋優對提高RFR建模算法在羊肉TVB-N含量方面的預測能力十分必要。借助二進制粒子群優化算法(Binary particle swarm optimization，BPSO)在高光譜特征波長提取方面的應用優勢[13-14]，本研究采用BPSO法和相關系數分析法(Correlation coefficient，CC)優選表征羊肉TVB-N含量的高光譜特征變量，并利用袋外均方根誤差RMSEOOB對RFR模型最佳回歸子樹和分裂特征2個重要參數進行尋優，建立羊肉TVB-N含量的最佳預測模型，旨在為高光譜遙感技術在肉品新鮮度檢測方面提供理論參考。

1 材料與方法

1.1 試驗材料與TVB-N測定

用羊肉樣本為察哈爾羊，取羊酮體里脊肉置于低溫冷藏箱運至實驗室。在無菌操作臺上將鮮羊肉剔除表面脂肪和肌膜，盡量保持樣本表面平整，用無菌刀分割成84塊，尺寸大小約為45 mm×45 mm×20 mm，自封保鮮袋密封后逐個編號，整齊無擠壓地擺放在貯藏溫度為4 ℃的冰箱環境中貯藏1～12 d。每隔24 h取出7個樣本，于室溫下靜置30 min，用濾紙吸收表面水分后對樣本進行光譜采集，隨后將樣本送至氨基酸實驗室，立即通過半微量凱氏定氮法[15]測定樣本TVB-N含量，作為衡量羊肉新鮮度的標準。參照國標GB/2707—2016《食品安全國家標準-鮮(凍)畜、禽產品》與前人研究成果[16-17]，以TVB-N含量為主要依據將羊肉鮮度劃分為3個等級：一級鮮度(TVB-N≤15 mg/100 g)，二級鮮度(15 mg/100 g 25 mg/100 g)。

1.2 高光譜圖像采集

采用五鈴光學(ISUZU OPTICS)高光譜成像系統，主要部件包括成像光譜儀、鹵素燈、電控位移平臺、光源控制器和計算機。整個系統置于黑箱內，系統主要軟件參數見表1。為保證獲取的羊肉樣本圖像清晰且不失真，預實驗反復調試后，將成像光譜儀的曝光時間設為2.1 ms，物鏡高度為40 mm，電控位移平臺速度為22.9 mm/s，起點和終點位置分別為165和235 mm。圖像分辨率選擇800×428像素，通過高光譜圖像采集軟件得到935～2 539 nm 的樣本高光譜圖像。

表1 高光譜成像系統主要元件信息Table 1 Information of main components of hyperspectral imaging system

1.3 光譜圖像黑白校正及感興趣區域選取

為消除光強的變化和鏡頭中暗電流對光譜信息的影響，在數據處理前對高光譜圖像進行黑白校正[18]：

(1)

式中：R為樣本黑白校正后光譜反射率；Is為樣本原始光譜反射率；ID為黑板校正光譜反射率；IW為白板校正光譜反射率。

使用ENVI 5.3軟件提取羊肉高光譜圖像信息，避開羊肉結締和筋腱部位，選取圖1所示的左上、左下、右上、右下和中間5個大小為20×20像素點作為感興趣區域(Region of interesting，ROI)，并計算ROI內所有像素點的平均值得到樣本平均反射光譜。

圖1 羊肉高光譜圖像感興趣區域
Fig.1 Hyperspectral image region of interest for lamb

1.4 BPSO法提取特征波長

為降低高光譜數據維數和冗余度，利用BPSO法[19]提取特征波長。二進制編碼將粒子位置向量中的每一位取1或0，1表示相應波長被選中，0則表明該波長未被選中。

算法執行過程如下：

1)初始化粒子群，設定粒子的初始位置xi=(xi1,xi2,…,xid)和初始飛行速度vi=(vi1,vi2,…,vid)；

2)計算適應度函數值Fitness；

(2)

(3)

(4)

(5)

圖2 粒子群提取特征變量流程
Fig.2 Flow chart of particle swarm optimization algorithm

1.5 隨機森林算法

隨機森林算法是Breiman提出的一種基于決策樹的集成學習算法[20]。利用自助抽樣法(bootstrap)從原數據集有放回地隨機抽取多個不同的訓練數據集，且每個訓練數據集的樣本數量與原數據集相同。利用隨機子空間法對每個bootstrap數據集分別構建決策樹模型，分別利用投票法和平均值法確定模型預測的分類和回歸結果。RFR模型中每顆決策樹都是回歸樹，這些樹并行建立多個預測子模型，構建過程如圖3所示。

圖3 RFR模型構建流程圖
Fig.3 Flow chart of RFR model construction

定理1 假設S為原始樣本，N為S中的樣本數，則S中每個樣本沒有被抽取的概率為：

(6)

根據定理1，RFR模型利用bootstrap法隨機抽取的自助訓練集中，每次約有36.8%的樣本未被抽取，這些樣本稱為袋外數據(Out of bag，OOB)。由于OOB誤差估計近似等于交叉驗證結果，計算森林中每棵回歸子樹的OOB估計誤差即可得到RFR模型的泛化誤差，因此，利用OOB對RFR算法的建模結果進行泛化誤差估計。回歸子樹個數和候選分裂特征個數是影響RFR模型預測精度的主要離散型調節參數，依據袋外數據均方根誤差RMSEOOB對上述2個參數尋優以提高RFR模型預測精度，尋優過程如圖4所示。

具體步驟如下：

1)設定回歸子樹個數k與候選分裂特征個數mtry初值，從M個輸入特征中隨機選擇mtry個特征建立RFR訓練模型，并計算RMSEOOB。

2)判斷RMSEOOB是否收斂，設定RMSEOOB收斂精度為10-3。如果RMSEOOB變化差值大于收斂精度，則進行模型參數優化，bootstrap重新采樣訓練RFR模型；當RMSEOOB變化差值達到設定精度時，確定模型回歸子樹數量k值。

3)設定分裂特征個數mtry的取值范圍為[1,M]，步長為1，對mtry進行全局迭代尋優，以RMSEOOB最小化原則確定最佳分裂特征個數mtry的值。

4)依據確定的調優參數生成回歸子樹，構建羊肉TVB-N含量的RFR預測模型。

圖4 RFR模型參數尋優流程圖
Fig.4 Flow chart of parameter optimization process for RFR model

2 結果與討論

2.1 羊肉TVB-N隨時間的變化趨勢

分析羊肉TVB-N隨儲存時間的變化關系(圖5)可知，TVB-N含量隨儲存時間的增加呈遞增趨勢。樣本存儲1 d后TVB-N平均含量為10.31 mg/100 g，從第1～ 3d，樣本TVB-N均值從10.31 mg/100 g緩慢增加至13.37 mg/100 g。當TVB-N含量積累達到一定值，假單細菌等有氧細菌的增加在一定程度上有效抑制了蛋白質代謝，進而延緩了TVB-N的增加速率[21]。從第4 d起TVB-N迅速增加至16.12 mg/100 g，超過一級新鮮度閾值15 mg/100 g。隨著保鮮袋內現有空氣的消耗，有氧細菌減少，乳酸菌等無氧細菌此時達到高峰，蛋白質分解速度進一步加快[22]，肉樣存儲9 d后TVB-N顯著增加到27.32 mg/100 g，樣本已經開始腐敗。

圖5 TVB-N含量變化趨勢
Fig.5 Variation trend of TVB-N content

2.2 羊肉樣本近紅外反射光譜分析

蛋白質主要含C-H、O-H和N-H等基團，近紅外光譜與食品分子含氫基團振動合頻與各級倍頻的吸收有關，因此通過分析近紅外光譜特征則可獲取樣品有機分子含氫基團特征信息[23]。圖6為第2、6和11 d樣本平均反射光譜曲線，其中每條特征曲線為當天樣本平均反射光譜。分析樣本整個近紅外譜區反射光譜發現，光譜在974、1 211及1 440 nm處存在3個明顯的特征吸收峰。974和1 440 nm附近存在強吸收峰，分別為水分子O-H伸縮振動二級和一級倍頻，表明水分子在該波長下對近紅外輻射存在強吸收。1 211 nm處的相對弱峰為C-H伸縮振動二級倍頻，被認為與脂肪含量相關。1 074 nm附近存在N-H 伸縮振動二級倍頻，而N-H一級倍頻存在于1 500 nm附近，在此波長下能獲得大量蛋白質相關信息[24-25]。分析圖6可知，不同樣本TVB-N光譜反射曲線整體趨勢一致，但光譜反射率會隨TVB-N含量高低發生變化，TVB-N高的樣本的光譜相對反射率較高，這應該與肉品新鮮度變化過程中樣本的脂肪、蛋白質等化學成分的變化有關，上述光譜特征為利用近紅外高光譜技術分析肉品新鮮度提供了理論依據。

11.86、18.93和35.17 mg/g分別為樣本第2、6和11 d的TVB-N值
圖6 第2、6和11 d樣本平均反射光譜曲線
Fig.6 Mean reflection spectrum on the 2nd, 6th and 11th day

2.3 特征波長提取

2.3.1BPSO法提取特征波長

試驗獨立尋優20次得到最優適應度收斂曲線如圖7所示，算法迭代至60次時，適應度函數值基本接近最優值。利用BPSO法優選出1 055、1 112、1 238、1 257、1 307、1 345、1 446、1 490、1 692、1 723、1 748、1 781、1 836、1 868、2 069和2 157 nm共16個特征波長。

圖7 適應度函數收斂曲線
Fig.7 Fitness function convergence curve

2.3.2CC法提取特征波長

對樣本平均反射光譜與其TVB-N測定值進行相關性分析，結果如圖8所示。不同波長下光譜反射率與TVB-N均為正相關關系，選擇相關系數值高于0.5的極值點波長變量，將其作為特征波長變量。最終通過CC算法篩選出特征波長個數為13個，分別為1 049、1 099、1 125、1 339、1 364、1 673、1 723、2 144、2 251、2 263、2 307、2 332和2 351 nm。

圖8 相關系數分析法提取特征波長
Fig.8 Correlation coefficient method for extracting characteristic wavelengths

2.4 預測模型建立

2.4.1數據集劃分

84個羊肉樣本中，其中2個樣本反光嚴重，2個樣本TVB-N測定值異常，去掉上述4個明顯離群樣本，共得到80個有效樣本。為提高模型的預測精度，確保校正集樣本所含信息的代表性，采用K-S(Kennard-Stone)算法基于樣本間的歐氏距離對樣本進行篩選以劃分校正集和預測集，56個樣本為校正集，24個樣本為預測集。K-S算法劃分羊肉樣本的TVB-N統計值如表2所示。

表2 羊肉樣本TVB-N統計Table 2 TVB-N statistics for lamb samples

2.4.2模型建立

以BPSO法優選特征波長為自變量，建立羊肉TVB-N含量的BPSO-RFR預測模型。首先確定回歸子樹數量k值，試驗預設k=500，將20次獨立運行試驗的平均值作為測試結果。輸入量M為BPSO法優選的16個特征變量，當候選分裂特征個數mtry分別取M，M/2，M/3時，RFR模型的RMSEOOB隨k的變化曲線如圖9(a)所示，由圖可知，所有曲線RMSEOOB均隨k的增加而降低，并總體呈現迅變-緩變-平穩的變化趨勢。變化初期k<50時，RMSEOOB隨k增加而迅速減小，k值超過50后，曲線逐漸趨于平緩，直至k增加到300，RMSEOOB值接近收斂，因此k值設定為300。遍歷分裂特征個數mtry在區間[1，16]的全部數值，步長設定為1，算法尋優結果如圖9(b)所示。當mtry=9時，BPSO-RFR預測模型的RMSEOOB取得最小收斂值4.85，模型訓練時間為0.243 s。綜上分析，當k=300，mtry=9時，可獲得最優的BPSO-RFR預測模型。

圖9 BPSO-RFR預測模型參數尋優過程
Fig.9 Parameter optimization process for BPSO-RFR prediction model

以CC分析法優選的13個特征波長為自變量，依照上述參數尋優過程建立羊肉TVB-N含量的CC-RFR預測模型。設定分裂特征個數mtry值分別為13，7和4，k值設定為300時，對應模型的RMSEOOB值已經全部收斂，因此k值設定為300。再遍歷分裂特征個數mtry在區間[1，13]的全部數值，步長設定為1，當mtry=7時，CC-RFR模型的RMSEOOB取得最小收斂值5.47，模型訓練時間為0.186 s。綜上分析，當k=300，mtry=7時，可獲得最優的CC-RFR預測模型。

2.5 預測模型結果分析

表3 不同特征變量和算法結合模型的檢驗結果Table 3 Results of the models combined with different characteristic variables and algorithms

BPSO-RFR和BPSO-PLSR模型的預測效果均優于CC-PLSR、CC-RFR模型，分析認為，CC分析法以樣本實測值與光譜信息的相關性最大原則，優選了表征力較強的16個特征波長，“開環”方式下特征變量只單向影響模型預測精度，預測結果的優劣信息無法反饋給模型，因此限制了光譜信息的有效提取，導致模型反演精度降低。BPSO法將模型預測精度作為特征波長的提取標準，“閉環”方式下智能地提取特征波長，將模型預測精度作為特征變量的提取標準，既提升了模型計算效率，又提高了模型穩定性和準確性。BPSO-RFR和BPSO-PLSR預測模型精度均略高于朱榮光等[5]利用全波段光譜信息建立的TVB-N含量預測模型，一定程度上表明BPSO算法提取特征波長的有效性，另外，該算法的優越性在水質預測[19]、烤煙煙葉圖像特征提取[26]及青貯飼料含水率檢測[27]等方面也得到了驗證。

3 結 論

TVB-N含量是衡量肉品新鮮度的關鍵指標。為快速、無損檢測羊肉新鮮度，本研究采集不同存儲天數羊肉樣本935～2 539 nm的近紅外高光譜圖像，分別以BPSO法和CC分析法優選特征波長為自變量，建立羊肉TVB-N含量的BPSO-RFR、BPSO-PLSR、CC-PLSR和CC-RFR預測模型，比較各模型預測精度并確定最佳預測模型，探索應用高光譜技術預測羊肉TVB-N含量的可行性，主要結論如下：

1)回歸子樹個數k為300，候選分裂特征個數mtry分別取9和7時，BPSO-RFR、CC-RFR預測模型獲得袋外均方根誤差RMSEOOB最小收斂值。表明實現最佳回歸子樹和分裂特征2個參數的智能尋優不僅可提高RFR建模算法的計算效率，還能改善BPSO-RFR、CC-RFR預測模型在羊肉TVB-N含量方面的預測效果。