高光谱技术结合迭代决策树的香肠菌落总数预测

郭培源，徐 盼，董小栋，许晶晶

（北京工商大学计算机与信息工程学院，食品安全大数据技术北京市重点实验室，北京 100048）

随着人们生活水平的提高以及愈来愈多食品问题的发生，食品安全逐渐成为人们重点关注的问题，因而食品检测的方法及其性能近年来成了研究热点之一。香肠作为种类丰富且深受大家喜爱的食品，其检测研究具有重要意义。受加工过程中不达标的卫生条件以及运输贮藏过程中环境因素的影响，香肠品质会有所下降，而带来很多安全问题和隐患[1]。评价香肠品质常用的指标是挥发性盐基氮含量、蛋白质分解产生、过氧化值和酸价含量[1-2]，以及亚硝酸盐含量，少量的添加可给香肠上色且提供独特的风味，过量使用，则会对身体造成危害。除此之外，香肠的菌落总数也是一项重要的评价指标，因为细菌会加速香肠的腐败，从而给人体的健康带来不利影响[3-4]。

香肠中菌落总数的测定，传统是采用理化实验的方法，即通过培养皿计数获得，但是用理化实验获取香肠菌落总数周期长、耗试剂、操作繁琐，且对样品具有破坏性[5]。而近年来兴起的高光谱成像技术是一种无损检测技术，与理化实验相比，它具有在线实时、对样品无破坏性、准确便捷等优点，现已广泛应用于食品检测领域[6-7]。王莉等[8]采用波长范围400～1 000 nm可见近红外高光谱对冷鲜羊肉的菌落总数和挥发性盐基氮含量进行新鲜度的检测研究，其中，采用偏最小二乘回归得到了最佳预测结果，建立的菌落总数和挥发性盐基氮含量预测模型的校正集相关系数分别为0.906 7和0.914 7，预测集相关系数分别为0.874 3和0.880 2。刘善梅等[9]采用高光谱成像技术对生鲜猪肉的含水率进行无损检测，建立偏最小二乘回归预测模型，交叉验证和预测相关系数分别为0.926和0.924，均方根误差（root mean square error，RMSE）分别为0.467%和0.438%。张雷蕾等[10]在470～1 000 nm波长范围内，从高光谱图像中提取反射光谱，对光谱进行多元散射校正（multiplicative scatter correction，MSC）处理，并采用偏最小二乘建模分析，实现对猪肉的新鲜度评价。Jin Huali等[11]利用偏最小二乘方法分别在400～1 000 nm全波段上和1 000～2 500 nm中选取的6 个特征波长上进行建模预测花生油中的成分含量，两种方法的效果都很好，但是后者的效果优于前者。Xiong Zhenjie等[12]采用偏最小二乘-连续投影算法的方法实现了红肉中色素含量的定量检测，并采用图像处理的方法将色素在红肉中的分布进行可视化研究。

虽然高光谱成像技术已广泛应用于食品检测领域，但利用高光谱技术检测香肠内化学物质的含量以及对香肠进行分级的相关研究与应用非常少。本实验采用400～1 000 nm高光谱仪采集香肠的高光谱数据，并分别采用迭代决策树（gradient boosting decision tree，GBDT）和支持向量回归（support vector regression，SVR）方法建立香肠菌落总数的预测模型，以期为香肠菌落总数的快速定量和品质控制提供参考。

1 材料与方法

1.1 材料与仪器

广式香肠，购于北京永辉超市，将香肠切块，每块香肠厚2 cm，获取50 份样品，其中每份样品有200 g，每天取一份样本，将样品平放于电移台上，采用“推扫式”成像的方法获取香肠的光谱值，然后进行菌落数测定。

肉制品光谱检测仪购自北京卓立汉光公司，波段范围400～1 000 nm，采样分辨率5 nm，共有128 个波段。高光谱成像系统硬件由装有图像采集卡的计算机、CCD相机、成像光谱仪、光源等系统组成。

1.2 高光谱检测原理

高光谱成像技术既可以获取含有物质内外理化信息的光谱值，同时也能通过成像设备获取样品各个波段的图像数据，这种图谱合一的三维数据称为“数据立方体”[13]，如图1所示。其中，图像代表两维的空间维度，而波长代表一维的光谱维度。“数据立方体”中每个波段可获取一幅二维图像，样品的每个像素可以获取一条光谱曲线[14-15]。图像信息能够全面反映物体的外在特征，而光谱信息则能够反映物体的内在物理结构和化学成分等信息[16]。

图1 高光谱成像技术检测原理Fig. 1 Detection principle of hyperspectral imaging technology

1.3 方法

1.3.1 高光谱数据采集

由于高光谱图像采集系统获得的原始高光谱图像在各个波段范围内的光源强度、光源亮度分布不均匀，并且暗电流和噪声等因素会对光谱信息产生影响[17]，因而需要对采集到的高光谱图像进行黑白板校正处理[18]，得到样品的光谱反射值，具体如式（1）所示：

式中：R为校正后图像；IR为原始图像；ID为黑板校正图像；Iw为白板校正图像。

使用高光谱分析处理软件ENVI5.1，在每个样本的高光谱图像上选取感兴趣区域（region of interest，ROI），对ROI采用N维可视化工具获取平均光谱曲线[19]，如图2所示。对50 个样本中每个样本选取10 个ROI，共得到500 个光谱数据。

图2 香肠样本的ROI（a）及其平均光谱曲线（b）Fig. 2 Region of interest in sausage samples (a) and its average spectral curve (b)

1.3.2 光谱预处理

在采样过程中，由于样品的不均匀性、高频随机噪声、基线漂移、光散射等因素会对建模效果产生负面影响[20]，所以为了减少或消除此类因素的影响，需对采集的原始高光谱数据进行不同的预处理，本研究采用MSC的预处理方法。MSC是高光谱建模最常用的预处理方法，分析结果较好[21]。它可以有效地消除样品颗粒分布不均匀或者样品大小不同等情况造成的散射误差。

首先计算样品得到的所有高光谱的平均光谱，将得到的平均光谱作为基准光谱。每个光谱与基准光谱进行一元线性回归运算，求得各光谱相对于基准光谱的线性平移量（回归常数）和倾斜偏移量（回归系数）。在每个原始光谱中减去回归常数且除以回归系数后，每个光谱的基线平移和偏移都得到了修正，而样品成分含量对应的光谱信息在数据处理的过程中没有受到影响，进而提高原始光谱的信噪比。平均光谱、回归方程和MSC运算的算法过程如（2）～（4）：

式中：Ai,j为香肠样品的平均光谱曲线；Ai（MSC）为经多元散射校正后的光谱；A为n×p维定标光谱数据矩阵；n为光谱数量；p为波长点数；A为原始光谱的平均矢量；Ai为1×p维矩阵，表示每个光谱矢量；mi为Ai与A线性回归得到的相对偏移系数，Bi为Ai与A线性回归获得的平移量。

1.3.3 主成分分析

高光谱数据的数据量较为庞大，且相邻波段的图像相互重叠，具有很大的关联性，因此高光谱数据降维处理的效果将影响后续的建模效果[22]，而主成分分析（principal component analysis，PCA）在数据降维方面具有独特的优势，所得的主成分分量之间相互独立，可以有效地消除高光谱数据中的冗余信息[23]。一般情况下，PC1包含波段中80%的方差信息，前3 个主成分包含了所有波段中90%以上的信息量[24]。本实验对光谱的128 个波段进行PCA，选取前5 个主成分建立模型，这5 个主成分的方差累计贡献率达到95%以上。

1.3.4 SVR

SVR是将复杂实际问题通过非线性变换转换到高维特征空间，在高维空间中构造线性决策函数来实现原空间中的非线性决策[24-27]。它是找到一个超平面，使到超平面最远的样本点的“距离”最小。本研究选用SVR方法实现香肠菌落总数的预测，所用的高斯核函数[28]如式（5）所示：

式中：σ为核函数的宽度参数；x、xi分别为超平面最远的样本点及中心点。

高斯核函数对数据中存在的噪声有着较好的抗干扰能力，由于其具有很强的局部性，其参数决定了函数作用范围，随着参数σ的增大而减弱。

SVR高斯核函数的c和g参数，常用粒子群（particle swarm optimization，PSO）算法、网格搜索、遗传算法（genetic algorithm，GA）3 种方法进行寻优。其中，c为惩罚系数，即对误差的宽容度。c越高，说明越不能容忍出现误差，容易出现过拟合；c越小，容易欠拟合。c过大或过小，泛化能力都差。g为选择高斯核函数后自带的一个参数，隐含地决定了数据映射到新的特征空间后的分布。g越大，支持向量越少；g越小，支持向量越多，支持向量的个数影响训练与预测的速度。

1.3.5 迭代决策树

迭代决策树（gradient boosting decision tree，GBDT）作为回归树的一种，相对于一般决策树算法具有防止过拟合、泛化能力较强等优点。模型训练的时候，对于输入的一个样本，首先会赋予一个初值，然后会遍历每一棵决策树，每棵树都会对预测值进行调整修正，每一次迭代是为了改进上一次结果，减少上一次模型的残差，并且在残差减少的梯度方向上建立新的组合模型[29]。其基本思想是通过构建M个弱分类器，经过多次迭代最终组合而成一个强分类器。

GBDT又被称为提升树，其可以表示为决策树的加法模型：

式中：T（xi,θm）为决策树；θm为决策树的参数；M为树的个数。

针对不同问题的GBDT，其主要区别在于使用的损失函数不同，包括用平方误差损失函数的回归问题，用指数损失函数的分类问题，以及用一般损失函数的一般决策问题。本研究使用平方误差损失函数实现回归。

提升树的训练流程如下：

输入：训练数据集T＝{（x1,y2）, （x2,y2）, …, （xN,输出：提升树fM（x）。

1）初始化f0（x）＝0；2）对m＝1, 2, …,M；a）按式（6）计算残差rmi＝yi-fm-1（xi），i＝1, 2, …,N；b）拟合残差rmi学习一个回归树，得到T（x,θm）；c）更新fm（x）＝fm-1（x）+T（x,θm）；3）得到回归问题提升

1.3.6 模型评价指标

1.3.6.1 RMSE

RMSE是观测值与真值偏差的平方和与观测次数n比值的平方根，它能很好地反映测量的精密度，RMSE越小，模型的预测效果越好。具体如式（7）所示：

式中：Xobs,i为观测值；Xmodei,i为真实值。

1.3.6.2 决定系数R2

决定系数用于判断回归方程的拟合程度，R2越接近1，模型的预测效果越好。具体如式（8）所示：

式中：yi为真实值；y为均值；为估计值。

1.4 数据处理

实验采用10折交叉验证的方法对原始数据进行处理，即每次用9 个子集的并集作为训练集，余下的1 个子集作为测试集，这样总共获得10 组训练/测试集，从而进行10 组训练和测试，最终训练与测试的结果返回的是10 组的均值。因而本实验中每次有450 个香肠光谱样本作为训练集，50 个作为测试集。对经MSC预处理和PCA降维处理后的光谱分别采用SVR和GBDT方法建立香肠菌落总数的预测模型，并验证模型的预测效果。采用RMSE和R2作为评价模型预测效果的指标，获得了较好的实验结果，其训练集和测试集均方根误差分别为0.001和0.003，决定系数R2分别为0.998和0.996。

2 结果与分析

2.1 光谱预处理的结果

原始光谱以及经MSC后的光谱图如图3所示。

图3 原始（A）和经MSC预处理后（B）的光谱图Fig. 3 Original spectra (A) and spectra preprocessed by MSC (B)

2.2 SVR建模结果

表1 不同参数寻优算法对应的SVR建模结果Table 1 Modeling results of SVR with different parameter optimization algorithms

3 种寻优算法得到的最优c和g参数及其对应模型的预测结果如表1所示。SVR方法采用PSO算法进行c和g参数寻优得到的预测模型最优。PSO算法的收敛速度快，受问题维数的变化影响较小，使得求解过程更容易计算[30]。采用PSO算法进行c和g参数寻优的SVR方法对香肠光谱训练集和测试集样本的菌落总数预测结果如图4所示。

图4 SVR模型的训练集（a）与测试集（b）的预测结果Fig. 4 Prediction results of training (a) and test (b) sets with SVR model

2.3 GBDT建模结果

迭代1 000、1 500 次和2 000 次的GBDT建模结果如表2所示。迭代2 000 次得到的建模结果最好，并且迭代过程很快。迭代2 000 次的GBDT方法对香肠光谱样本的训练和测试结果如图5所示。

表2 不同迭代次数的GBDT建模结果Table 2 Modeling results of GBDT with different iterations

图5 GBDT模型的训练集（a）与测试集（b）的预测结果Fig. 5 Prediction results of training (a) and test (b) sets with GBDT model

2.4 两种方法建模结果的比较

由SVR和GBDT的建模结果，比较采用PSO算法进行参数寻优的SVR建模结果与迭代2 000 次的GBDT建模结果可知，GBDT建模结果要远优于SVR的，GBDT所得的RMSE非常小，比SVR所得的要小一个数量级，并且GBDT所得的R2几乎为1。除此之外，SVR建模所需的训练时间很长，GBDT训练时间则很短。因而基于高光谱成像技术利用GBDT方法预测香肠菌落总数的方法可行且有效。

3 结 论

本实验通过高光谱成像系统采集50 个香肠样本的高光谱数据，并利用高光谱分析处理软件ENVI5.1，在每个香肠样本的高光谱图像中选择10 个ROI，从而获得500 个香肠样本的平均光谱数据。实验采用MSC方法对光谱预处理，并采用PCA方法从128 个光谱波段中选择5 个特征波段，从而提高了模型的预测精度。以处理过的光谱数据作为输入，理化实验所得的香肠菌落总数值作为输出，分别采用SVR和GBDT方法建立香肠菌落总数的预测模型。实验结果可知，迭代2 000 次的GBDT建模结果最优。本实验中，GBDT模型迭代2 000 次时，训练集和测试集的RMSE都很小，R2也都接近1，但是当迭代次数多于2 000 次时，是否会产生过拟合、建模效果需要进一步论证。除此之外，进一步需要探索研究地是，将每个像素点下预测出的菌落总数定量反演到香肠样本表面图像上，生成可视化分布图，使香肠新鲜度的动态变化趋势更加直观、形象地呈现。