预测聚类树(predictive clustering tree, PCT)^[23]是一种基于决策树的学习方法, 样本x会根据节点上的规则落到决策树的某一叶子结点上, 叶子结点预测x对应的标记的概率值, 按照概率值可确定样本x属于哪一个类别。 考虑到单棵决策树的性能有限, 可将决策树集成得到森林应用于多标记分类问题中, 如基于PCT的随机森林(random forest of predictive clustering trees, RF-PCT)^[24], RF-PCT通过平均每棵树的预测结果输出一个概率向量, 按照其中概率值大小, 确定样本属于某个类别, 如图1所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 基于PCT的随机森林Fig.1 The random forest based on PCT

深度森林用于多标记分类问题, 需要考虑如何利用标记之间的相关性以提高分类器的性能。 将RF-PCT作为基分类器嵌入到深度森林中, 可以在逐层特征学习中挖掘多个标记之间的相关关系。

利用聚类树构造层森林, 每层包含2个完全随机森林和2个普通随机森林, 再用层森林构造layer-by-layer多层级联的深度森林模型, 模型结构如图2所示。 模型第一层随机森林的输入是经过预处理的冷鲜羊肉样本特征光谱数据, 输入数据经过第一层中各森林的计算输出不同的结果, 将结果拼接形成候选特征空间H_t, 特征筛选用于充分挖掘多个新鲜度评价指标的相关性, 将上一层的候选特征空间H_t经过度量指标的判定, 保留相关性较大的元素, 筛选后的特征空间为G_t。 由于模型训练到第一层时, 还未形成上一层的候选特征空间, 故在第一层G_t=H_t。 为保证原始特征, 每一层筛选后的特征空间G_t都与预处理后的特征光谱数据拼接共同作为下一层的输入。 层增长控制通过度量指标计算每一层的分类性能来判断模型是否充分学习, 进而确定层数来控制模型的复杂度。 最后模型的输出为概率值, 依据概率值确定所属类别。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 基于多标记深度森林的新鲜度评价模型结构Fig.2 The structure of freshness evaluation model based on multi-label deep forest

特征筛选根据不同的度量指标计算每一层随机森林输出的置信度, 通过比较置信度, 重新组成特征输入到下一层中。 特征筛选的核心是定义置信度计算方法。

假设评价模型中每一层的输出为H_t, 该值由若干森林的输出拼接得到, 通过森林数量对H_t求平均值得到该层的预测概率矩阵P, P的行数为样本数, 列数为标记数。 当度量指标基于实例时, 把矩阵P的每一行元素按照从大到小排序; 当度量指标基于标记时, 把矩阵P的每一列元素按照从大到小排序。

Hamming loss用来判断P上的分类是否正确。 假设阈值θ =0.5, p_ij> 0.5时预测结果为1, 该值越大, 证明预测为1的概率越大, 故置信度越大; p_ij≤ 0.5时预测结果为0, 该值越小, 证明预测为0的概率越大, 故置信度越大。 因此Hamming loss置信度可以定义为

αj=1N∑i=1Npij(1)

式(1)中, N为实例数, pij=pij, pij> 0.51-pij, pij≤0.5。

One-error用来判断相关标记中预测的最大概率, 故One-error置信度可以定义为预测最大概率值

αi=maxpij(2)

Ranking loss用来判断样本的所有标记的排列顺序, 由Ranking loss定义可知, Ranking loss为0时, 模型性能最佳, 故在定义置信度时, 需要列出Ranking loss为0时的各种组合。 若存在4个标记, 则可能的组合有{0000, 1000, 1100, 1110, 1111}五种, 通过计算这些组合的概率之和得到Ranking loss置信度。

αi=∑j=0Q∏k=0jpik∏k=j+1Q(1-pik)(3)

Macro-AUC用来判断标记上所有样本的排列顺序, 与Ranking loss类似, Macro-AUC为1 时, 模型性能最佳, 故在定义置信度时, 需要列出Macro-AUC为1 时的各种组合, 计算这些组合的概率之和, 得到Macro-AUC置信度。

αj=∑i=0N∏k=0ipkj∏k=i+1N(1-pkj)(4)

通过置信度对特征进行筛选, 表1为特征筛选过程。

表1

Table 1

表1(Table 1)

表1 特征筛选过程 Table 1 The process of feature screening Input: Measure M, forests' output Ht, previous representation Gt-1 1. Initialize matrix Gt=Ht 2. If Measure M is label-based then 3. For j=1 to l do 4. Computer confidence αjt on H·jt 5. Update G·jt to G·jt-1 when αjt< θ t 6. End for 7. End if 8. If Measure M is instance-based then 9. For i=1 to m do 10. Compute confidence αit on Hi·t 11. Update Gi·t to Gi·t-1 when αit< θ t 12. End for 13. End if Output: Representation Gt

表1 特征筛选过程 Table 1 The process of feature screening

假设训练集为X, 标记集为Y, 指定一个度量指标M, 第t层的随机森林输出为H^t, 度量值为 mjt, 第t层的阈值为θ _t。 将H^t的值赋给G^t, G^t用来表示一组新的特征, 若度量是基于标记的, 则对于某一标记的所有实例, 计算其置信度 αjt; 若度量是基于实例的, 则对于某一实例的所有标记, 计算其置信度 αit, 如果置信度小于阈值θ _t, 则将G^t中对应的元素替换成上一层的特征表示G^t-1中的对应元素, 如图3。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 特征筛选Fig.3 Feature screening

根据森林的输出H^t和标记集Y计算当前层的度量值, 在计算度量值时也需要将基于实例的度量和基于标记的度量分开计算。 每一层的阈值将由该层的度量值和置信度决定, 若当前层特征上的度量值小于前一层特征上的度量值, 则将当前层特征上的置信度存储到集合S中, 最后对S中的置信度求平均值作为当前层的阈值。

为了降低模型过拟合的风险, 达到控制模型复杂度的目的。 提出的评价模型采用K折交叉验证将训练数据分为K个组, 对于每一组数据, 使用其他所有组的数据进行训练, 并对当前组数据进行预测。 表2为层增长控制流程, 假设模型最大深度为T, 训练集为X, 标记集为Y, 度量指标为M, 以及包含每层度量值的数组q, 其中性能最好的值定义为q_best。 当评价模型通过特征筛选得到第t层的输出H^t之后, 根据度量指标M计算该层的度量值q[t], 若q[t]大于q_best, 则更新q_best值; 若q[t]连续三次小于q_best并且t在T的范围内, 则停止层的增长, 同时保留包括q_best所在的层与前面的所有层, 删除后面所有层。

表2

Table 2

表2(Table 2)

表2 层增长控制流程 Table 2 The process of layer growth control Input: measure M, maximal depth T, training data X, label data Y 1. Initialize pasameters: performance array in each layer q[1: T] the best performance qbest the first layer threshold θ 1 the best performance layer index L=1 the model set C=ϕ 2. For t=1 to T do 3. Get classifier ht by training the t-layer forests 4. Ht=ht ([X, Gt-1]) 5. When t> 1 6. Compute θ t 7. Gt=filtered characteristics 8. Compute performance q[t] on measure M with Gt 9. if q[t] is better than qbest, 10. Let qbest=q[t], L=t 11. else if qbest is not updated in recent 3 layers then 12. break 13. end if 14. Add the t-layer to model set C=C∪ layert 15. end for 16. Keep {layter1, …, layerL} in model set C and drop others {layerL+1, …} Output: model set C, the best performance layer index L

表2 层增长控制流程 Table 2 The process of layer growth control

试验所用冷鲜羊肉样本取自内蒙古锡林郭勒盟苏尼特右旗农贸市场, 选择屠宰后经过排酸的5只羊酮体里脊部位, 去除脂肪和结缔组织, 均匀分割成6 cm× 6 cm× 1 cm的肉片, 用保鲜袋分3组密封包装、 并编号, 无挤压放置在温度为4 ℃的冰箱中, 存放14 d。 每隔24 h取一次样本, 在室内放置25 min, 挥发掉样本表面的水分, 用于挥发性盐基氮(total volatile basic nitrogen, TVB-N)、 酸碱度(pH值)、 菌落总数(total aerobic count, TAC)、 大肠菌群近似数(approximate number of coliforms, ANC)测定和光谱反射率采集, 试验样本覆盖了新鲜、 次新鲜和不新鲜三个冷鲜羊肉新鲜度等级, 新鲜羊肉表面有光泽, 肉细而紧密; 不新鲜羊肉表面无光泽, 肉色深暗, 肉质松弛无弹性; 次新鲜为过渡阶段, 无论表面情况或触摸手感均处于新鲜与不新鲜之间。

样品的TVB-N含量依据GB/5009.228— 2016《食品安全国家标准食品中挥发性盐基氮的测定》中的半微量凯氏定氮法测定^[25]; pH值依据GB/5009.237— 2016《食品安全国家标准食品pH值的测定》中的非均值化试样测定法测定^[26]; TAC含量依据GB/4789.2— 2022《食品安全国家标准食品微生物学检验菌落总数测定》中单位质量菌落总数标准值进行检测试验^[27]。 ANC含量依据GB/4789.3— 2016《食品安全国家标准 食品微生物学检验 大肠菌群计数》中单位质量大肠菌群数标准值进行检测试验^[28]。

图4为冷鲜羊肉样本在14天内四项指标的变化趋势, 由图可知, 贮藏初期样本中各指标含量较低, 随时间增加, 各指标含量逐渐增加, 其中TVB-N在第5天、 第10天和第14天增幅较大; pH值在第4天和第10天增幅较大; TAC含量在第5天增幅较大; ANC含量在第4天和第9天增幅较大。 可以初步判断冷鲜羊肉样本分别在第4~5天、 第9~14天内新鲜度等级发生变化。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 各新鲜度指标变化趋势 (a): TVB-N变化趋势; (b): pH变化趋势; (c): TAC变化趋势; (d): ANC变化趋势Fig.4 The trends of each freshness indicator (a): The trend of TVB-N; (b): The trend of pH; (c): The trend of TAC; (d): The trend of ANC

根据国家食品卫生监测标准和以往的研究成果, 当TVB-N含量≤ 15 mg· (100 g)^-1时, 为新鲜肉; 当15 mg· (100 g)^-1< TVB-N含量≤ 25 mg· (100 g)^-1时, 为次鲜肉; 当TVB-N含量> 25 mg· (100 g)^-1时, 为不新鲜肉。 目前, 我国食品卫生检测标准中尚未制定冷鲜羊肉的菌落总数、 pH和大肠菌群近似数指标, 故以TVB-N含量作为主要指标, 通过比较TAC、 pH值和ANC指标测量值共同确定样本的新鲜度等级。

高光谱采集系统包括照明设备、 机械扫描平台、 高光谱成像仪(Hyperspec VNIR N-series)、 反射参考板和图像采集软件, 光谱仪可采集的波长范围在400~1 000 nm, 共有750个光谱通道, 分辨率为2.8 nm。

每次试验时, 提前30 min打开光谱仪预热, 将样本放置在距光谱仪镜头约40 cm处, 设置像元混合次数为6次, 光谱仪曝光时间为3 ms, 调节光谱像元亮度(DN)值小于8 500。 试验时, 利用调焦板调节光谱仪, 设置光谱仪扫描方向、 次数和移动速度, 采集黑白校正光谱图像, 用于获取样本校正光谱数据。 利用ENVI软件, 从样本的每一个光谱图像中随机选取20个感兴趣点作为特征提取和校正模型建立试验数据。 图5为实验样本部分感兴趣区域原始光谱DN值。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 部分感兴趣区域原始光谱反射DN值Fig.5 Original spectral reflection DN values in the rigin of interest from 400 to 1 000 nm

黑白校正后的光谱数据仍存在部分噪声且光谱强度有差异, 需要对其进行预处理, 因此, 采用卷积平滑法(Savitzky-Golay)进行平滑滤波, 再采用多元散射校正对平滑滤波后的光谱数据进行处理, 消除光谱中的基线平移或偏移现象, 提高光谱信噪比, 用于后期的特征提取和分类识别。 经过平滑滤波和散射校正后的光谱曲线如图6所示。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 预处理后光谱反射DN值Fig.6 The spectral reflection DN values after pretreatment

试验采用连续投影算法对样本光谱数据进行特征提取, 假设数据集X中的样本数为M, 原始特征数为J, 选择的第一个波段为i(0), 算法在每次迭代中合并新的波段, 直到集合中存在N个波段, 算法流程如表3。

表3

Table 3

表3(Table 3)

表3 连续投影算法 Table 3 Continuous projection algorithm Input: Data X(M× J), the first wavelength i(0), the number of wavelengths N 1. Let xj=jth column of X; j=1, …, J 2. Let a set S to store unselected wavelengths S={j, 1≤ j≤ J, j∉{i(0), …, i(n-1)}} 3. Take the first selected wavelength i(0) as the normal line and establish the orthogonal subspace. Calculate the projection of xj on the subspace orthogonal as Pxj=xj-(xTjxk(n-1))xk(n-1)(xTk(n-1)xk(n-1))-1, j∈S 4. Let i(n)=arg(max‖ Pxj‖ , j∈ S) 5. Let xj=Pxj, j∈ S 6. Let n=n+1. If n< N, go back to Step 2 Output: {i(n); n=0, …, N-1}

表3 连续投影算法 Table 3 Continuous projection algorithm

若N和i(0)是未知的, 则对N定义一个范围N_min≤ N≤ N_max, 对每个N, 需要考虑初始波段i(0)从1到J的每一种情况, 进行步骤2— 步骤6计算, 根据输出结果i(n)建立多元线性回归分析模型, 以i(n)对应的光谱数据作为测试集, 以TVB-N、 pH、 TAC和ANC含量为标记, 计算均方根误差(root mean square error, RMSE), 其最小值对应的i(0)和N即为最优初始波段和选择的波段数。 试验设置特征波长个数范围为5~30, 共提取了18个特征波段, 如图7所示。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 提取特征波段Fig.7 Extracting characteristic bands

采用本文提出的深度森林建立冷鲜羊肉新鲜度评价模型, 将280个光谱样本按照3∶ 1的比例划分为训练集和测试集, 训练集和测试集样本个数分别为196个和84个, 表4为训练集和测试集的新鲜、 次新鲜和不新鲜样本数统计结果。

表4

Table 4

表4(Table 4)

表4 冷鲜羊肉不同新鲜度样本数统计 Table 4 Sample number statistics of different freshness of chilled mutton 样本集 新鲜 次新鲜 不新鲜 总计 训练集 47 92 71 210 测试集 13 28 29 70 总计 60 120 100 280

表4 冷鲜羊肉不同新鲜度样本数统计 Table 4 Sample number statistics of different freshness of chilled mutton

将理化微生物实验方法测得的TVB-N、 pH、 TAC和AVC四项新鲜度评价指标按照国家食品卫生监测标准和前人研究成果, 分别划分出新鲜、 次新鲜和不新鲜的区间, 并组成标记, 图8为建立新鲜度分类模型所用的冷鲜羊肉样本标记。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 某冷鲜羊肉样本标记Fig.8 The label of chilled mutton sample

其中三种不同的颜色对应三个新鲜度等级, 绿色表示新鲜区间, 黄色表示次新鲜区间, 红色表示不新鲜区间。 新鲜度分类模型中的所有标记均按照此规则产生, 若将其视为一个长度为12的数组Z。 当样本新鲜度为新鲜时, 绿色单元格所在的标记为属于该样本的标记, 其数组为[1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0]; 当样本新鲜度为次新鲜时, 橙色单元格所在的标记为属于该样本的标记, 其数组为[0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0]; 当样本新鲜度为不新鲜时, 红色单元格所在的标记为属于该样本的标记, 其数组为[0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1]。

在新鲜度等级预测过程中, 若四项新鲜度评价指标含量都处于新鲜区间, 则表明样本为新鲜; 若其中一项或多项处于次新鲜区间, 则表明样本为次新鲜; 若其中一项或多项处于不新鲜区间, 则表明样本为不新鲜。

试验建立的深度森林参数设置如下, 将最大层数设置为10, 每层的森林数为2, 分别为PCT组成的一个随机森林和一个极端随机森林, 其中每一个森林设置5棵树, 后面的每一层都比前一层多5棵树, 该方法可以保证模型在每一层都能学习到不同的表示。 同理, 森林的最大深度为3, 后面的每一层都比前一层的森林最大深度多3个单位。 最后设置5折交叉验证以防止过拟合。