基于多重二阶差商滤波的矿物油检测自动判别方法

关键词 矿物油测定；液相色谱-气相色谱联用法；鼓包峰轮廓；二阶差商；滤波

由饱和烃（MOSH）和芳香烃（MOAH）组成的矿物油污染物在食用油和相关加工食品中普遍存在[1-3]。污染源可能来自食品包装材料，也可能是食用油生产、运输过程中工业油残存导致的污染[4-8]。毒理学研究结果表明， MOSH 容易蓄积， MOAH 具有致癌和遗传毒性。近年来，欧盟对矿物油的监管日趋严格。2022 年，欧盟植物、动物、食品和饲料委员会SCoPAFF 发布了《新型食品和食物链的毒理学安全-2022/04/21》[9]，规定了不同脂肪含量食品中MOAH 的限量，即低脂肪含量（≤4%脂肪）的干性食品中MOAH 限量为0.5 mg/kg，较高脂肪含量（4%～50%脂肪）的食品中MOAH 限量为1.00 mg/kg，脂肪和油脂（≥50%脂肪）中MOAH 限量为2 mg/kg[9]。我国国家标准GB 2760—2024《食品添加剂使用标准》允许食品级矿物油（白油和液体石蜡）作为加工助剂使用，在糖果和鲜蛋中的最大使用量为5.0 g/kg[10]。然而，GB/T 21309—2007《涂渍油脂或石蜡大米检验》和GB/T 37514—2019《动植物油油脂矿物油的检测》两项检测标准对矿物油进行定性判断的检出限为百分级别，远不能满足欧盟限量标准要求[11-12]。

目前，矿物油定量分析主要采用液相-气相色谱联用技术（LC-GC-FID）[13-17]，分析得到的MOSH 和MOAH 谱图均呈不规则的驼峰形状。根据FID 的等碳响应，在驼峰外添加非矿物油组分的烃类-内标物，利用其峰面积即可实现矿物油的定量分析[18]。驼峰上方的尖峰或窄峰为天然烃类干扰物，通常不会影响驼峰轮廓的判别[1]。但是，当干扰物净化不完全时，会在驼峰上方形成一个或数个豁口，干扰驼峰轮廓的准确判断与勾勒，由于人们对驼峰的轮廓判断存在主观差异，导致矿物油分析结果的重现性较差。Bauwens 等[19]提出了一种去除窄尖峰、保留宽驼峰的方法，在每个二维（2D）点重建一维（1D）色谱图（类似于1D 轴上的投影），并根据设置的“基线宽度”消除窄峰。在保留宽驼峰的同时，作为驼峰顶部的峰值出现的信号会自动移除。但是，该方法可实现的功能较少，而且对某些较宽的尖峰去除效果不理想。

本研究针对现有问题设计了一套基于多重二阶差商滤波的驼峰曲线形态的自动判别算法。由于尖峰的下极值点高于驼峰曲线，而驼峰曲线较平滑，相对驼峰上的点，尖峰所处的下极值点的二阶差商通常为绝对值较大的负数。因此，通过过滤二阶差商较大的点可以达到去除尖峰的目的。此外，由于豁口是由尖峰簇组成的，通过多次过滤二阶差商较大的点，还可以去除豁口。结果表明，对驼峰曲线进行正向、反向、正向三次二阶差商过滤，即可去除色谱数据中的豁口。根据上述原理设计了相应的矿物油自动分析软件，在使用软件处理驼峰数据时，用户可以设置不同的过滤“阈值”，也可在自动形成的驼峰折线上手动删除不必要的数据点。在进行积分运算时，用户可以选择使用“自动基线”或“基线数据”，也可以选择“按驼峰积分”或“碳谱积分”以确定不同的积分区间。经实验验证，本算法及软件可以有效地处理矿物油检测的色谱数据，在保证计算精度的基础上大幅提高计算效率，缩短数据处理的时间。

1 实验部分

1.1 仪器与试剂

液相色谱-气相色谱联用仪（LC-GC，日本岛津公司）：配备二元泵和紫外检测器的LC20A 液相色谱仪和配备氢火焰离子化检测器（FID）的GC 2010 plus 气相色谱仪；LC-GC 接口（德国Axel Semrau 公司）和PAL 自动进样器（瑞士CTC 公司）；Biofuge primo R 低温高速离心机（美国赛默飞世尔科技公司）；R-100 旋转蒸发仪（瑞士步琦有限公司）；TUS-200P 恒温金属浴（上海一恒科技有限公司）。

正己烷、二氯甲烷和无水乙醇（农残级，德国默克公司）；异丙醇（色谱级，美国赛默飞世尔科技公司）；3-氯过氧苯甲酸（m-CPBA， 70%，日本东京化学工业公司）；硫代硫酸钠（分析纯，国药集团化学试剂有限公司）。MOSH/MOAH 标准溶液（北京坛墨质检科技有限公司），包括：正十一烷（n-C11）、正十三烷（n-C1₃）、环己基环已烷（Cycy）、5α-胆甾烷（Cho）、正戊基苯（5B）、1，3，5-三叔丁基苯（TBB）、1-甲基萘（1-MN）、2-甲基萘（2-MN）和苝（Per）。其中， n-C11、Cycy、5B、TBB、1-MN、2-MN 的浓度为300 mg/L， Cho 和Per 的浓度为600 mg/L， n-C₁₃ 的浓度为150mg/L， Cycy 为MOSH 定量内标， TBB 或2-MN 为MOAH 定量内标，其余内标和C₇～C₄₀的正构烷烃混合物（德国默克公司）用于评估LC-GC-FID 仪器性能[13]。硅胶60（0.063～0.200 mm，德国默克公司）使用前在400 ℃活化16 h，冷却后置于干燥器内，7 d 内使用。BAM-K009 润滑油标准品（德国柏林联邦材料研究与测试研究所）。

1.2 样品

2023年在电商平台购买了42 个食品样品，包括32 个植物油样品、6 个植物油基鱼罐头样品、2 个婴幼儿配方奶粉样品和2 个奶酪样品。所有样品在保质期内进行分析。参考样品RVEP233004 椰子油和RVEP200612 奶粉从德国标物局（DRRR）获得；参考样品葵花籽油和橄榄油为法国油脂研究所（ITERG）在ISO 20122 标准国际验证试验中提供的样品。

1.3 实验方法

1.3.1皂化和萃取

根据ISO 20122 标准对食用油样品进行前处理[20]。称取3 g 食用油样品至100 mL 离心管中，加入20 μL MOSH/MOAH 标准溶液和30mL 正己烷-乙醇（1∶1， V/V）混合溶剂，充分混合后移取10 mL 至另一个离心管中，加入3mL KOH 溶液（33%， m/V）， 60 ℃恒温振荡皂化30 min 直至溶液澄清；冷却至室温后，向皂化液中加入5 mL 正己烷和5 mL 乙醇-水（1∶1， V/V）混合溶剂，振荡，离心10 min；将正己烷层转移至梨形烧瓶中，再次用5 mL 正己烷萃取残留物；合并2 次正己烷提取液，转移至3 g 硅胶柱，用15 mL正己烷-二氯甲烷（7∶3， V/V）混合溶剂洗脱；收集所有洗脱液，浓缩至1～2 mL。

参照文献[21]的方法进行奶粉前处理。称取5 g 样品至离心管，加入10 μL MOSH/MOAH 标准溶液和10 mL 50 ℃温水，振荡2～3 min；加入10 mL KOH 溶液（0.6 g/mL）和20 mL 乙醇， 60 ℃恒温振荡皂化30 min；冷却至室温后，加入15 mL 正己烷，振荡5 min，加入10 mL 去离子水，振荡、离心，取上清液；向残留物中加入10 mL 正己烷重复提取，合并正己烷上清液，减压浓缩至2～3 mL；加入5 mL 乙醇-水（1∶1， V/V）混合溶剂，振荡、离心；将所有上清液转移到3 g 硅胶柱中，用15 mL 正己烷-二氯甲烷（7∶3，V/V）混合溶剂洗脱，收集所有淋洗液，浓缩至1～2 mL。

其它食品基质的前处理参照文献[22]的方法进行。将样品均质后，称取5 g 样品至离心管中，加入10 μL MOSH/MOAH 标准溶液和20 mL KOH-水-乙醇溶液（25∶50∶90， m/V/V）， 60 ℃恒温振荡皂化30 min后，提取步骤与奶粉样品一致。

1.3.2 环氧化

样品中的生物烯烃会干扰MOAH 分析，必须通过环氧化将其去除[23]。移取1.3.1 节所得提取物1 mL 至离心管中，加入0.5 mL m-CPBA-乙醇溶液（200 g/L， m/V）， 40 ℃振荡15 min 后，加入0.5 mL 乙醇和2 mL 100 g/L 硫代硫酸钠溶液，继续振荡30 s，终止反应，离心，取上清液进行后续分析。

1.4 LC-GC分析

MOSH 和MOAH 的分离分析通过LC-GC-FID 联用系统进行。液相色谱的色谱柱为Allure Si HPLC色谱柱（250 mm×2.1 mm， 5 μm， 60 Å，美国瑞思泰克公司），流动相为正己烷和二氯甲烷。将50～100 μL样品注入液相色谱仪，在进样1.5 min 内，以300 μL/min 的流速将流动相从最初的100%正己烷转变为70%正己烷，保持至6.2 min；采用500 μL/min 的流速仪用100%二氯甲烷反冲9 min；切换至100%正己烷，以500 μL/min 的流速正向平衡色谱柱10 min，降低流速至300 μL/min，直至程序结束。MOSH 馏分的洗脱时间为2.0～3.5 min， MOAH 馏分的洗脱时间为4.5～6.0 min，通过紫外检测器（230 nm）监测MOSH 和MOAH 的流出情况， MOSH 和MOAH 馏分通过HPLC-GC 接口从HPLC 全部转移至GC，随后，通过保留间隙技术和部分并行溶剂蒸发技术将大量的溶剂蒸汽排出，溶质保留在分析柱头。保留间隙系统包括MXT Siltek 无涂层的预柱（10 m×0.53 mm，美国瑞斯泰克公司）、MXT-1 100%二甲基聚硅氧烷分析柱（15 m×0.25 mm×0.25 μm，美国瑞斯泰克公司）和溶剂排除出口（SVE）。GC 在MOSH 和MOAH 馏分转移时压力分别为84 和90 kPa（氢气），大量溶剂蒸汽从SVE 出口排出， SVE 出口在馏分转移前0.5 min 打开，在转移结束0.3 min 后关闭。SVE 出口关闭后， MOSH 和MOAH 通道的分析压力升至150 kPa 并保持恒定。柱温箱升温程序设定为：初始温度60 ℃保持6 min，以15 ℃/min 升至120 ℃，再以25 ℃/min 升至370 ℃，保持6 min。FID 温度为380 ℃。氮气、氢气（燃烧气）和空气（助燃气）的流速分别设置为30、40 和400 mL/min。

1.5 数据处理

人工积分使用Lab solutions 软件（Version 5.92，日本岛津公司）进行数据处理。分别以Cycy、TBB 或2-MN 作为MOSH 或MOAH 的定量内标对MOSH 和MOAH 进行定量分析。通过积分计算鼓包峰与基线之间的空白面积与内标峰的面积比获得矿物油的含量。

基于多重二阶差商滤波的驼峰曲线形态的自动判别算法，将原始谱图数据导入软件后，自动积分得到MOSH 或MOAH 含量。

2 结果与讨论

2.1 矿物油含量自动计算算法

在矿物油检测中，由于仪器配套的数据分析软件不具备自动建立矿物油谱线的功能，需要检测人员手工建立矿物油谱线，这给矿物油检测工作增加了难度和工作量。Bauwens等[19]提出通过在每个2D点重建1D色谱图（类似于1D轴上的投影）来消除尖峰的方法。该方法可以提取出粗略的矿物油色谱图，但对有些较宽的尖峰去除效果不佳。例如，在文献[19]中，在1D 时间1300s附近，原始色谱数据（红色曲线）包含一个较宽的尖峰，应被平滑去除，但该方法得到的矿物油色谱图（绿色曲线）未能将其去除。

本研究设计了一种思路简单、功能完备且能准确去除多种尖峰的算法。首先，通过对矿物油色谱数据进行谷值点搜索建立驼峰的轮廓；然后，对谷值点进行三次二阶差商过滤，去除无效的谷值点（包括尖峰和豁口），得到色谱数据正确的下包络线，即矿物油谱线；最后，通过对矿物油谱线和基线之间进行面积积分，根据用户输入的内标质量、样品质量和矿物油谱线面积即可以计算出样品中的矿物油的含量。本方法不仅能够自动建立矿物油谱线，也能自动建立基线，还能根据需要对自动建立的谱线对参数进行人工微调。实验结果表明，本方法能很好地提高矿物油含量分析的效率和精确度。

2.2 矿物油谱线的自动生成算法

矿物油谱线的生成主要分为2步：（1）寻找有效区域内的谷值点，建立矿物油驼峰的基本轮廓；（2）对基本轮廓中的无效点进行过滤，确定准确的矿物油谱线。

2.2.2 去除无效的谷值点

虽然谷值点基本反映了矿物油驼峰的形状轮廓，但有一些数据点明显不属于矿物油在色谱图中形成的驼峰。如图3 中远高于大部分绿色数据点的尖峰数据点，单个的峰以尖峰形式出现，簇状的峰则以豁口的形式出现，这些都不属于矿物油在色谱图中的驼峰。观察这些无效的谷值点可发现一个共同特征，即二阶差商是绝对值较大的负数，说明这些无效谷值点在谷值点序列中既是局部极大值点，又是数据系列的拐点。因此，可以采用对谷值点序列进行二阶差商过滤的方法去除这些无效的谷值点。序列的二阶差商计算方法如下：

其中， D₀是根据经验设置的判别阈值，为负数。当谷值点序列中一个点的二阶差商小于D₀时，即可认为该点为无效的谷值点，将其从序列中删除。通过调节阈值D₀可设置算法对豁口的屏蔽程度：阈值的绝对值越小，屏蔽程度越大；阈值的绝对值越大，屏蔽程度越小。

从谷值点序列中过滤二阶差商绝对值较大的点，可以去除单个的尖峰。但是，豁口是由一簇尖峰形成的，一次差商过滤难以完全去掉整个豁口。实验发现，对谷值点序列进行正向、反向、正向3次差商过滤，即可基本去除色谱数据中的所有豁口。对于一些特殊的豁口，如果采用默认参数去除效果不佳，还可通过适度调节阈值D0，增大或减小对豁口的去除程度，优化矿物油谱线。

图3中紫色数据点及其形成的折线即为矿物油气相色谱图中的驼峰，也是色谱图数据对应的正确的矿物油谱线。对于较宽的尖峰，本算法也能很好地去除豁口，获得准确的驼峰谱图。因此，本方法能够有效地去除多种尖峰。

在实际处理数据时，若通过调整阈值得到的矿物油谱线未达到最佳，本程序还提供手动删除无效数据点的功能，从而使矿物油谱线更符合实际。

2.3 基线的确定和积分计算

计算矿物油谱线对应的积分数值，还需确定积分使用的基线，通过计算矿物油谱线和基线之间包络的面积换算出矿物油检测的含量。

2.3.1 完全自动生成基线

大部分情况下，积分使用的基线就是包络线序列的起始点和终止点之间连接的直线。但是，矿物油谱线驼峰对应的中间谷值点有时会低于包络线序列的起始点和终止点之间连接的直线，这时基线将变为包络线序列的起始点、终止点以及驼峰对应的中间谷值点3 点构成的折线，如图4所示。

2.3.2 处理加载的基线数据

如图5A 所示，实验数据形成的基线存在干扰，不利于积分计算，为了确保检测结果的可靠性，有必要对基线数据进行曲线拟合。通过观察可知，二次曲线拟合基本可以满足形状上的要求（图5B）。

在计算数值积分时，由于矿物油谱线和自动生成的基线数据均表现为折线形状，需要在具体的积分位置计算折线上对应的y 值。

其中， j 是折线上对应的数据点的编号。对于加载的数据基线，可以记录二次曲线拟合的3 个参数，直接采用二次曲线公式计算数值积分对应的y值：

2.3.3 计算积分面积

确定了矿物油谱线和基线即确定了数值积分需要的上、下边界，就可以计算矿物油谱线与基线间包含的面积，依据积分面积可换算出矿物油检测的含量。通常， MOSH 选择Cycy 作为定量内标， MOAH 选择2MN 或TBB作为定量内标。在加载检测数据时，软件可直接读入色谱数据中的内标峰面积。在处理色谱数据时，手动输入内标质量M_ISTD（ng）和样品质量M_smp（g）。根据矿物油谱线面积A_smp 和标准物质的峰面积A_ISTD，利用式（8）计算出矿物油在样品中的含量。

根据上述原理自动生成基线，获得矿物油谱线和基线数据，对二者之间包含的区域进行自动数值积分，再根据所得到的积分面积换算出矿物油的含量。与文献[19]的方法相比，本方法功能更完备，更利于提高用户操作效率。

2.4 算法流程

（1）加载数据。首先按固定格式将色谱数据转存为ASCII 码，采用软件加载ASCII 码数据文件。在加载了色谱数据文件的基础上，继续加载基线数据和碳数谱图的数据。（2）计算有效区间。利用式（2）和式（3）计算矿物油检测数据的有效区间，提取对应的数据编号。（3）计算谷值点。利用式（4）计算矿物油色谱数据在有效区间内包含的谷值点。（4）计算矿物油谱线。设置合适的阈值，利用式（5）删除谷值点序列中的无效数据点，必要时手动删除谷值点序列中无效数据点。（5）确定基线数据。可以根据矿物油谱线的起始点、终止点和中间谷值点构建基线，也可由用户直接加载基线数据。基线数据需要在有效区间内采用最小二乘法进行二次曲线拟合，同时保留3 个参数。（6）计算数值积分。缺省状态下，积分区间为矿物油谱线的全部有效区间，可根据碳数谱图数据选择相应的积分区间。数值积分时，利用式（6）和式（7）计算y 值。（7）计算含量。根据用户输入的内标质量M_ISTD（ng）、样品质量M_smp（g）、矿物油谱线面积A_smp 和标准物质的峰面积A_ISTD，利用式（8）计算样品中的矿物油含量。（8）显示并保存数据。计算结果直接在窗口中显示，用户可以选择保存计算结果到数据文件。本方法已全部用软件实现并投入使用中，软件操作界面如图6所示。

2.5 算法验证与应用

测定了德国标物局（DRRR）国际能力验证的食用油（RVEP233004）和ISO 20122 国际标准验证的葵花籽油和橄榄油样品，相关数据分别使用传统的色谱工作站人工处理和本研究建立的计算机软件自动计算，结果见表1。采用z′比分数衡量实验室测定值与参考值的偏离程度，人工处理和软件自动计算的z′比分数的绝对值均小于2，表明两种数据处理方式均能满足能力验证要求。

以市售米粉样品为实际样品，采用润滑油标准品（MOSH∶MOAH=89∶11）制备2 个浓度水平的加标样品，用于加标回收实验。按照1.3节方法处理后，注入HPLC-GC-FID 进行测定，采用所建立的计算机软件自动计算，结果见表2，MOSH 和MOAH 的加标回收率分别为82.4%～108.7%和92.9%～95.7%，RSD≤3.5%。

选择食用油、奶粉、奶酪和鱼罐头等真实样品基质共42 个样品，进行相应的前处理操作，将供试液注入HPLC-GC 联用仪器，分别进行人工积分和软件自动积分，再剔除低于定量限的数据，得到42 组MOSH 数据和28 组MOAH 数据。采用SPSS 软件对人工积分和软件积分的数据进行处理，结果见表3，MOSH 和MOAH 不同处理方式的差值均符合正态分布，采用配对t 检验方法进行显著性差异分析，结果见表4， pgt;0.05，表明两种数据处理方式不存在显著性差异。

检测结果表明，本算法适用于食用油、乳酪和奶粉等基质样品中定量限以上的MOSH和MOAH检测。

3结论

建立了一种矿物油检测色谱数据分析的自动判别方法。分析结果表明，本方法不仅原理简单，而且比现有方法更准确，功能更完备，在很大程度上提高了矿物油色谱数据分析的准确性和工作效率。（1）构建的驼峰生成方法符合矿物油数据的实际形态，可准确地描绘矿物油在色谱数据中的实际区域。通过构建矿物油色谱数据的“谷点”序列，在“谷点”序列中排除二阶差商较小的“凸点”，有效消除了干扰数据，最终建立了符合实际的矿物油谱线。（2）开发的数据处理软件工具可以在矿物油色谱数据中自动形成积分计算的有效区间和基线，还可以针对加载的基线数据，通过最小二乘法在有效区域局部形成基线的二次曲线。本软件可由使用者根据碳数谱图确定积分区间，方便确定碳数区间对应的物质含量。本方法是按照传统数据分析算法设计的，还存在少量需要手动删除部分干扰峰的情况。在积累大量的实验数据后，通过机器学习等人工智能方法，开发更准确高效的数据分析模型，为矿物油的数据分析提供更科学高效的分析工具。