基于高光谱成像技术的羊肉新鲜度预测

张 珏 田海清 王 轲 张丽娜 于 洋

(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含量方面的预测效果。