基于极限学习机的低温液态乳制品中金黄色葡萄球菌污染程度预测方法

王建新，张懿文，连梦雪，王晔茹

(1.北京林业大学信息学院，北京100083；2.国家食品安全风险评估中心，北京100021)

金黄色葡萄球菌较易污染蛋白质含量较高的食品，比如鱼、肉、蛋、乳及其制品[1]。与其他制品不同，乳及乳制品在制作完成后一般不再进行二次烹饪，更容易受到金黄色葡萄球菌及其所产肠毒素的威胁，由其引发的食物中毒事件在全球范围内普遍存在[1-2]。我国巴氏杀菌乳国家标准(GB 19645—2010)中规定，金黄色葡萄球菌在同一批次采集5个样品，均不得检出。因此，对金黄色葡萄球菌进行风险评估具有重要的现实意义。

我国目前所进行的微生物风险评估多数聚焦在零售和消费环节，较少从源头上即从原料及生产加工环节开展风险评估[3]。此外，我国微生物风险评估的基础数据大多来源于食品监测和专项调研数据，在食品生产加工、运输贮藏等多个环节中的关键数据普遍存在缺失情况[4]。关键环节数据的获取成为微生物风险评估的难点。在无法获取充足数据的情况下，对食品风险评估系统进行多环节定量预测研究成为研究者关注的重点。

近年来，低温液态乳制品已成为居民日常高频消费的食品[5]，其运输多为冷链物流，对温度要求较为严格[6]。我国的冷链运输尚处于起步阶段，冷链运输过程中的监控不够严密或监控缺失，对于制冷机关闭、未提前预冷等情况无法做到准确把握[7-8]，这使低温液态乳制品容易受到金黄色葡萄球菌等食源性致病菌的污染，也会在一定程度上限制食品相关商贸供应链发展[6,8]。建立并且完善整个冷链链条，实现数据和信息的一体化是解决问题的关键[9]。虽然已有试点的一体化监控冷链[10-11]，但由于技术和成本等原因，大面积的推广应用还需要相当长的时间。因此，需要基于不完善的低温液态乳制品冷链监控数据，充分利用智能算法，达到或接近获取完善的冷链监控信息的目的。

1 预测方法和模型的选择

影响低温液态乳制品中金黄色葡萄球菌繁殖的因素主要集中在初始污染浓度、温度、湿度、pH、繁殖时间和容器类型上。这些因素对应的变量是原因变量。在历年国家食品污染和有害因素风险监测以及体系检查过程中，管理部门积累了原因变量和对应结果变量(即时段末端的最终污染浓度值)的大量数据，研究人员可充分利用这些数据，找到原因和结果之间的数量关系，在允许存在部分数据缺失的情况下，根据数据中蕴含的规律补全缺失值，从而实现预测和诊断的目的。

线性回归的方法是建立原因变量和结果变量之间关系的传统方法，但由于该方法仅限于线性关系，不能建立复杂的、高度非线性的关系，因此，对于从多原因变量数据预测结果变量的微生物浓度预测需求，线性回归方法是不合适的。BP(back propagation)神经网络，能够反映高度非线性的映射关系，在乳制品质量安全风险预警中起到了显著作用[12]。但是，BP神经网络的训练过程复杂、耗时较长，当数据量较大时，这个问题更为严重。极限学习机是神经网络的一种，已在很多领域得到了广泛的应用，取得了良好的效果[13]。极限学习机有诸多优势：①极限学习机的结构简单，只有三层网络结构，一层输入层、一层隐藏层、一层输出层；②传统的神经网络模型训练一般是迭代执行的，需要经历大量的迭代步骤，因而耗费大量的训练时间，一般是几十分钟甚至更长，而极限学习机的训练是基于最小二乘法的，训练时间一般能控制在1 s之内；③极限学习机的预测准确度通常也比较高[14]，这是因为所有的训练数据一次性使用，避免了因训练数据集中的部分噪声数据而使得待估计参数远离中心点。综上，本研究把极限学习机作为预测低温液态乳制品中金黄色葡萄球菌浓度和进行条件推断的计算工具。

与其他神经网络模型类似，极限学习机并不能很好地拟合指数型映射关系，但微生物繁殖过程大致遵循指数增长函数，因此，本研究把极限学习机模型与指数增长函数相结合，也就是将微生物繁殖机制嵌入极限学习机模型中，更好地拟合原因变量和结果变量之间的复杂映射关系，以构建低温液态乳制品多环节金黄色葡萄球菌浓度超标的预警模型，试图通过训练、计算和预测，改善现阶段有关数据匮乏的现状，提供针对每个环节预测的工具。

2 数据来源与优化方法

2.1 数据来源

本研究基于2017年低温液态乳制品中金黄色葡萄球菌的专项监测及相关的环境数据(合作工厂提供)，其中包含有初始污染浓度、温度、湿度和繁殖时间等变量，数据涵盖了生产加工、运输、存储和消费等多个环节，其中，最终浓度值是结果变量，其他全部是原因变量。因真实数据保密要求，笔者在真实数据基础上通过函数映射和增加随机扰动，生成了模拟数据。具体模拟数据格式见表1。

表1 低温液态乳制品中金黄色葡萄球菌相关数据格式

从表1可以看出：原因变量和结果变量之间没有直观的映射关系和规律，下一步需用采用非线性的映射模型。

2.2 哑变量变换、数据集分割和数据归一化

在原因变量中，“容器类型”是一个类型变量，其值包括A、B和C共3种，但这些值是类型，不是数值，不能直接参与极限学习机的训练和预测，需要把它转变为数值型变量。转变的规则是：如果类型变量有几个值，就转变为几个数值型变量，每个变量只取0、1。这样相对应的数值型变量称为哑变量。如果类型变量的值为A，则哑变量a的值为1，其他哑变量b和c的值为0。转变前后的变量值如表2所示。

a、b、c这3个哑变量并不是独立的，由其中任意2个哑变量的值可以得到另外1个哑变量的值[15]。这是因为3个变量中有且仅有1个变量的值为1，其余变量的值为0。例如，如果a和b的值分别为0和1，则可以推断c的值为0；如果a和b的值都是0，则可以推断c的值为1。因此，笔者只保留2个哑变量a和b，舍弃哑变量c。

表2 类型变量转变为哑变量

为了验证算法的准确率和性能，笔者把数据集分为2个部分：训练集和测试集，其中，训练集约占2/3，用于训练模型；其余约1/3用于模型测试和比较。

对所有的输入数据，不在[0，1]区间的属性值要通过缩放和平移，转化为[0，1]范围内属性值，并作为输入值输入极限学习机算法中去，对训练集和测试集中的数据必须做统一的缩放和平移操作才能保证算法的正确性。

2.3 总体框架

检测低温液态制乳品中金黄色葡萄球菌浓度的具体环节是分别检测加工前浓度、生产加工后/运输前的浓度、运输后/消费前的浓度、食用前的浓度。金黄色葡萄球菌的风险过程描述如图1所示。

图1 低温液态乳制品中金黄色葡萄球菌污染的风险过程描述Fig.1 Risk profile of S.aureus in low temperature dairy products

加工前浓度指的是生牛乳中金黄色葡萄球菌的浓度。我国生牛乳中携带金黄色葡萄球菌的比例较高[1]。生牛乳的加工过程，包括对生牛乳进行预热、加热灭菌、冷却、包装等操作，在这一系列过程中，生产加工工具带菌、操作人员带菌、产品加热灭菌环节不充分、产品包装不严密等情况均可能导致金黄色葡萄球菌的污染[16]。

在运输贮藏过程中，对于温度的要求较高。在运输过程中，运输车辆如果无法达到规定要求的温度，那么温度的升高易导致金黄色葡萄球菌的较快繁殖。消费者购买乳制品并进行贮藏时，如无法确保处于低温状态，也可能导致金黄色葡萄球菌的大量繁殖。

从生产加工阶段起，通过分析与金黄色葡萄球菌增殖相关的关键因素，即初始污染菌浓度、时间、温度、湿度、pH及容器类型等，获得金黄色葡萄球菌在不同影响因素数值组合下产生的不同预测结果，其中，每环节最后的输出结果值作为下一环节的初始输入值，采用既定的算法及其组合，利用已有的数据进行训练得到网络结构，构建预测评估模型。如果需要预测评估多环节中的金黄色葡萄球菌及肠毒素浓度，则需要进行两个方面的选择：一是选取待预测数据集文件；二是选取相对应的网络结构。其中，待预测数据集文件中包含的数据缺少最终的结果，而网络结构中包含训练学习得到的极限学习机的大量连接权重，二者结合才能够得到相应的预测数据集。

功能之相反，从菌浓度等因素倒推时间、温度等有时未知的因素，可以判定具体环节的具体原因变量出现问题，从而能够查找问题根源，划清责任。所用逆向的诊断方法与正向的预测方法完全一致，因此，本文重点阐述预测方法。

2.4 算法选择和改进方法

从表1可以看出，相关变量之间是典型的多个原因变量影响一个结果变量的关系。虽然极限学习机模型能够拟合这种数量关系映射，但是现有的模型存在固有的缺陷：预测过程是一个黑盒，对外不可观察，只考虑输入和输出，并不考虑问题本身的显式数据规律。所以，针对具体的问题，极限学习机在预测准确率提升方面还有改进的空间。

金黄色葡萄球菌繁殖呈指数增长，时间是影响金黄色葡萄球菌繁殖的最直接因素[17]，但被极限学习机作为普通因素对待；指数增长规律也是该研究中的一条影响变化的深刻规律，被极限学习机视为数据间的一般规律。

为此，本研究改变极限学习机的结构，把时间变量和指数增长规律进行单独计算后,再参与到极限学习机的运行过程，使得微生物繁殖的基本规律与极限学习机有机地结合起来。结合的主要原则是减少时序因素对极限学习机的影响，直接处理指数关系映射问题，从而得到更精准的预测结果。

为了实现本研究的算法，实验用计算机的配置是台式机，Windows10操作系统，CPU为英特尔I7处理器，DDR3内存16G，编程语言为Java1.8。

3 结果与讨论

3.1 改进的计算算法

时段末端的最终污染菌浓度用式(1)估算。

(1)

式中：i为初始污染菌浓度，c为最终污染菌浓度，T为繁殖倍增周期，t为繁殖时长。

式(1)得到的值作为极限学习机的一个输入，但要在原极限学习机中删除时长这个输入变量，因为该变量已经在指数函数中体现了。

为了阐述清楚二者结合的机制，示例如下。假设金黄色葡萄球菌的初始污染菌浓度为1 000 CFU/g，时间为78 h，繁殖倍增周期为36 h，上述数值都纳入极限学习机的输入中。如果使用微生物繁殖机制与极限学习机相结合，则这些值都要从极限学习机的输入变量中删除，继而替换为通过式(1)计算获得的数据，由此得到

二者结合前后，极限学习机的结构变化如图2所示。

图2 极限学习机与微生物机制相结合前后的结构变化Fig.2 The working structure change before and after combination of microbial mechanism with the extreme learning machine

极限学习机与微生物机制结合的食品风险预测模型构建步骤如下：

(2)

2)设置隐藏层的结点个数为N，设置输入层与隐藏层之间的随机权重为wi=[wi1,wi2,…,win]T∈Rn，设置偏置值bi,i=1,2,...,N。

运用g(x)公式计算出隐藏层的输出矩阵H和权值β(β=H+Y)，其中，H+是H的摩尔-彭罗斯广义逆矩阵。极限学习机中输入节点和隐藏节点之间的权重被设置成随机值，在[0,1]之间。所有隐藏节点与输出节点之间的权重则是通过对训练数据进行机器学习后获得。

3)构建极限学习机与微生物机制结合的食品风险评估模型，得到预测数值。基于极限学习机的预测值为

(3)

其中，

(4)

3.2 储存环节极限学习机和改进算法计算结果对比

根据金黄色葡萄球菌检测数据，对单独使用极限学习机进行运算和使用极限学习机与微生物机制结合的改进算法运算进行对比实验。以储存环节为例，展示储存环节的检测和预测数据，结果如图3所示。图3比较了极限学习机结合微生物繁殖原理前后对金黄色葡萄球菌浓度的预测准确度。图3的纵坐标分别表示实验检测的真实值、把储存条件数据输入极限学习机而计算获得的值、把储存条件数据输入改进的极限学习机而计算获得的值。

图3 储存环节极限学习机与微生物繁殖机制结合前后的预测结果比较Fig.3 Comparison of prediction results between pre-and post-combination of the extreme learning machine and microbial reproduction mechanism during storage

由图3可知：2种算法都有比较好的预测性能，比较接近实验检测值。

由于各数据组中的数据并不是由相同条件产生的，也就是初始污染浓度、温度、湿度和容器类型等条件不完全一致，本研究分别计算不同条件下相对于不同基准值的预测值百分比误差，然后取平均值得到预测的平均百分比误差(μpy)，具体计算见式(5)。

(5)

式中：xi和yi分别表示真实浓度值和预测浓度值，i=1,2,…,n；其中，n是数据的组数。

由对比实验结果可知，2种方法都能够比较准确地预测各环节的浓度值。但是，结合微生物繁殖机制之后，虽然有些环节预测准确性基本保持不变，但有些环节的预测准确度有了明显的提升，μpy更小，从14.77%降至2.85%。

极限学习机的工作原理是非线性加权，会削弱明确的指数关系(菌群繁殖与时长的关系)的影响，因此将极限学习机与微生物机制结合，将时间作为一个输入属性，代入生物繁殖规律函数中，进而把极限学习机输入值中初始污染和时间2个变量转换为由微生物繁殖规律所产生的函数计算值。将明确的关系不参与到极限学习机的学习和计算中，而是直接获取，去除了时间因素对极限学习机的影响，成功解决了极限学习机和其他工具不能很好地处理指数关系映射的问题。

研究结果也表明：数据训练集的数量会很大程度上影响预测和诊断结果的质量；而且，训练集中涉及的条件范围越广、组合越复杂，预测效果就会越好。因此，后续需要在加强各个环节数据的收集和融合基础上，对算法进行更进一步的完善和验证。

4 结论

本文将极限学习机与微生物机制相结合，应用于食品安全风险评估中，以金黄色葡萄球菌多环节风险预警为例，利用极限学习机和微生物繁殖机制结合的改进算法对金黄色葡萄球菌在各个环节、不同影响因子值下的生长繁殖情况进行拟合、预测和推断，从而构建了不同情境下金黄色葡萄球菌浓度的预测模型。该模型可以对未检测或缺失的变量值做出比较准确的估计，训练执行效率和预测执行效率都非常高。本研究可为低温液态乳制品多环节中金黄色葡萄球菌的预测及诊断提供快速的预筛功能，迅速定位问题环节，对实验检测起到便捷的辅助作用。