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

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

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

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

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

1 材料與方法

1.1 試驗材料與TVB-N測定

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

1.2 高光譜圖像采集

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

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

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

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

(1)

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

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

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

1.4 BPSO法提取特征波長

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

算法執(zhí)行過程如下：

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

2)計算適應(yīng)度函數(shù)值Fitness；

(2)

(3)

(4)

(5)

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

1.5 隨機森林算法

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

圖3 RFR模型構(gòu)建流程圖
Fig.3 Flow chart of RFR model construction

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

(6)

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

具體步驟如下：

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

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

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

4)依據(jù)確定的調(diào)優(yōu)參數(shù)生成回歸子樹，構(gòu)建羊肉TVB-N含量的RFR預(yù)測模型。

圖4 RFR模型參數(shù)尋優(yōu)流程圖
Fig.4 Flow chart of parameter optimization process for RFR model

2 結(jié)果與討論

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

分析羊肉TVB-N隨儲存時間的變化關(guā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。當(dāng)TVB-N含量積累達到一定值，假單細菌等有氧細菌的增加在一定程度上有效抑制了蛋白質(zhì)代謝，進而延緩了TVB-N的增加速率[21]。從第4 d起TVB-N迅速增加至16.12 mg/100 g，超過一級新鮮度閾值15 mg/100 g。隨著保鮮袋內(nèi)現(xiàn)有空氣的消耗，有氧細菌減少，乳酸菌等無氧細菌此時達到高峰，蛋白質(zhì)分解速度進一步加快[22]，肉樣存儲9 d后TVB-N顯著增加到27.32 mg/100 g，樣本已經(jīng)開始腐敗。

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

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

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

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法提取特征波長

試驗獨立尋優(yōu)20次得到最優(yōu)適應(yīng)度收斂曲線如圖7所示，算法迭代至60次時，適應(yīng)度函數(shù)值基本接近最優(yōu)值。利用BPSO法優(yōu)選出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 適應(yīng)度函數(shù)收斂曲線
Fig.7 Fitness function convergence curve

2.3.2CC法提取特征波長

對樣本平均反射光譜與其TVB-N測定值進行相關(guān)性分析，結(jié)果如圖8所示。不同波長下光譜反射率與TVB-N均為正相關(guān)關(guān)系，選擇相關(guān)系數(shù)值高于0.5的極值點波長變量，將其作為特征波長變量。最終通過CC算法篩選出特征波長個數(shù)為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 相關(guān)系數(shù)分析法提取特征波長
Fig.8 Correlation coefficient method for extracting characteristic wavelengths

2.4 預(yù)測模型建立

2.4.1數(shù)據(jù)集劃分

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

表2 羊肉樣本TVB-N統(tǒng)計Table 2 TVB-N statistics for lamb samples

2.4.2模型建立

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

圖9 BPSO-RFR預(yù)測模型參數(shù)尋優(yōu)過程
Fig.9 Parameter optimization process for BPSO-RFR prediction model

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

2.5 預(yù)測模型結(jié)果分析

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

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

3 結(jié) 論

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

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