拉曼光谱结合光谱特征区间筛选算法快速定量鉴别植物调和油品质

吴升德，姜 鑫，李爱琴，郭志明，朱家骥,*

（1.盐城市产品质量监督检验所，江苏 盐城 224056；2.盐城工学院电气工程学院，江苏 盐城 224051；3.江苏大学食品与生物工程学院，江苏 镇江 212013）

植物调和油是食用植物油市场的重要产品之一，它由两种或两种以上纯植物油按照一定比例混合而成。一般以低价值植物油为主体，掺入一定比例的高价值植物油（如特级初榨橄榄油、山茶籽油等）[1]。然而，一些不法商家往往通过虚假宣传夸大植物调和油中高价值植物油的含量从而牟利。因此，明确植物调和油中高价值植物油的含量，对于保障消费者权益以及维护市场秩序具有十分重要的现实意义。目前，常规的植物调和油鉴定方法主要包括气相色谱-质谱（gas chromatographymass spectrometry，GC-MS）法[2-3]、液相色谱-质谱法[4-5]、高效液相色谱法[6-7]、核磁共振波谱技术[8-9]等。这些方法虽然具有灵敏度高、准确性好的优点，但是通常需要复杂的样品前处理，从而导致检测过程繁琐、检测时间长，无法满足现场快速鉴定的需求。

近年来，食品安全分子光谱检测技术成效显著，例如，拉曼光谱作为一种新兴的分子光谱分析技术已在食品质量与安全检测领域获得了广泛的应用。拉曼光谱不仅可以提供待测物质丰富的分子结构信息，而且具有破坏性小、检测速度快、操作简单、不受水分子干扰等优点[10]。然而，拉曼光谱数据作为一种高维数据矩阵，直接对其建模分析将面临过拟合的风险。为了克服该问题，已引入了偏最小二乘回归（partial least squares regression，PLSR）和主成分回归（principal component regression，PCR）建模方法，这两种方法通过提取隐变量实现了光谱数据的降维，且取得了较好的应用效果[11-12]。但是，越来越多的研究表明对高维光谱数据进行特征变量的筛选能够进一步提高PLSR或PCR模型的性能[13-15]。

近年来，随着人工智能技术的迅猛发展，群体智能优化算法也蓬勃兴起，粒子群优化（particle swarm optimization，PSO）算法是其中的典型代表。受鸟群捕食行为的启发，Kennedy等[16]于1995年首次提出了PSO算法，其基本思想为通过群体中粒子之间的协作和信息共享实现最优解的搜索。目前，基于PSO的光谱变量筛选算法已被大量提出。例如，Xue Long等[17]提出了一种基于可见-近红外光谱与PSO-PLSR算法的快速无损检测脐橙表面敌敌畏残留的方法，采用PSO算法对采集的脐橙表面敌敌畏残留的可见-近红外光谱进行特征变量的筛选，并在此基础上建立PLSR模型，与直接建立的PLSR模型相比，特征变量筛选后建立的PLSR模型性能得到了明显的提高。Zhao Jie等[18]提出了一种基于近红外光谱与机器学习算法快速定量检测无糖养胃颗粒中活性成分的方法，采用PSO算法对采集的无糖养胃颗粒的近红外光谱进行特征变量的筛选，并在此基础上建立支持向量机模型，实现了3 种活性成分——芍药内酯苷、芍药苷和苯甲酰芍药苷的高精度定量检测。尽管PSO算法具有鲁棒性好、易于实现、全局搜索能力强等优点，但是PSO容易陷入局部最优解，从而降低了其寻优的性能。2014年，Mirjalili等[19]受狼群等级制度及捕猎行为的启发提出了灰狼优化（grey wolf optimizer，GWO）算法。该算法收敛性能好、参数少、易于实现、局部搜索能力强，但其全局搜索能力一般。显然，GWO能够与PSO形成良好的互补，从而提高群体智能优化算法的性能。

本研究旨在提出一种基于拉曼光谱与变量筛选算法对植物调和油中高价值植物油含量快速定量检测的方法，从而实现对植物调和油品质的定量鉴别。首先，针对PSO算法与GWO算法的弊端，将PSO与GWO融合构建混合智能优化算法，即PSOGWO算法；其次，将PSOGWO结合组合移动窗口（combined moving window，CMW）策略构建新型的光谱特征区间筛选算法，即PSOGWO-CMW算法；然后，将配制的玉米油（corn oil，CO）-特级初榨橄榄油（extra virgin olive oil，EVOO）植物调和油作为检测样本，并采集其拉曼光谱。为了评估PSOGWO-CMW模型的性能，将PLSR、PSO-CWM、GWO-CMW和PSOGWO-CMW分别用于检测自制CO-EVOO植物调和油样本中EVOO含量，并对检测结果进行对比分析。最后，将本方法与标准检测方法（GC-MS）分别用于检测真实CO-EVOO植物调和油样本中EVOO含量，并对结果进行对比分析。

1 材料与方法

1.1 材料

鲁花牌CO和欧丽薇兰牌EVOO购于当地永辉超市。

1.2 仪器与设备

RMS1000手持式拉曼光谱仪 上海如海光电科技有限公司；JP-010S超声波水浴振荡器 深圳洁盟清洗设备有限公司；6890N/5975N GC-MS联用仪 美国安捷伦公司。

1.3 方法

1.3.1 CO-EVOO植物调和油制备

将CO 与E VOO 按以下6 种比例（V（CO）∶V（EVOO）＝95∶5、90∶10、85∶15、80∶20、75∶25、70∶30）配制成CO-EVOO植物调和油。以95∶5为例，具体配制步骤如下：1）取9.5 mL的CO与0.5 mL的EVOO混合于干净的烧杯中；2）将上述混合物超声匀质5 min。其他比例的CO-EVOO植物调和油均按此方法配制。

1.3.2 光谱数据采集

采用RMS1000手持式拉曼光谱仪采集CO-EVOO植物调和油的主要参数设置如下：激光功率设置为100 mW，积分时间设置为2 s，扫描次数设置为3 次。对于每种比例的CO-EVOO植物调和油均采集10 条拉曼光谱，总共获得了60 条拉曼光谱。RMS1000手持式拉曼光谱仪采集的拉曼光谱范围为200～3 000 cm-1，光谱分辨率为2 cm-1。为了后续的定量分析，现将已获得的拉曼光谱数据集划分为校正集与预测集：1）对于每种比例的CO-EVOO植物调和油，从中随机挑选6 条拉曼光谱划入校正集，故校正集中拉曼光谱的数量为36；2）对于每种比例的COEVOO植物调和油，将剩余的4 条拉曼光谱划入预测集，故预测集中拉曼光谱的数量为24。

1.3.3 算法原理

1.3.3.1 CMW策略

许多研究表明筛选光谱特征区间比筛选离散的光谱特征变量更有意义，故本研究采用智能优化算法结合CMW策略筛选拉曼光谱的特征区间[20]。CMW策略的基本思想为在整个光谱范围内设置N个等宽的窗口（每个窗口代表一个光谱特征区间），这些窗口可以移动且可以相互覆盖。考虑到拉曼光谱的自身特性，每个窗口的宽度设置为5 个波数点[21]。此外，PSO-CMW、GWOCMW和PSOGWO-CMW算法各自所对应的最佳窗口数量将在本研究2.2节中进行优化。

1.3.3.2 PSO-CMW算法

本研究中，PSO-CMW算法的基本思路为通过PSO结合CMW筛选出最佳的光谱特征区间组合。PSO-CMW算法的具体实施步骤为：

1）假设pi代表第i个粒子，其可以表示为pi＝{wi1,wi2,…wij,…wiN}，其中wij表示第i个粒子的第j个搜索维度（即第j个窗口的中心位置），第i个粒子的速度vi＝{vi1,vi2,…vij,…viN}，初始化的粒子群由n个粒子构成。

2）为每个粒子建立PLSR模型（即目标函数），并计算5折交互验证均方根误差（root mean squared error of cross-validation，RMSECV）作为适应度值。记录每个粒子各自的RMSECV值、位置以及速度作为个体最优解，记为Pbesti。将最小RMSECV值所对应的粒子作为全局最优解，记为Gbestg。Pbesti和Gbestg将在迭代过程中更新。

3）根据Pbesti和Gbestg，对每个粒子的位置及速度进行更新，其数学表达式如式（1）、（2）所示：

式中：t为当前迭代次数；ω为惯性权重；c1和c2为加速度系数；r1和r2为[0,1]之间的随机数。

4）该算法持续迭代（重复执行步骤2、3）直到满足预先设置的最大迭代次数，最终，该算法输出最佳的光谱特征区间组合。

根据参考文献[22]，PSO-CMW算法的主要参数设置如下：惯性权重ω＝2，加速度系数c1＝c2＝2，种群规模n＝100，最大迭代次数Imax＝1 000。PSO-CMW算法的流程如图1所示。

图1 PSO-CMW算法流程图Fig.1 Flow chart of PSO-CMW algorithm

1.3.3.3 GWO-CMW算法

GWO算法是基于狼群等级制度和捕猎行为的智能优化算法[23]。本研究中，GWO-CMW算法的基本思路为通过GWO结合CMW筛选出最佳的光谱特征区间组合。GWO-CMW算法的具体实施步骤为：

1）假设gi代表第i匹灰狼，其可以表示为gi＝{wi1,wi2,…wij,…wiN}，其中wij表示第i匹灰狼的第j个搜索维度（即第j个窗口的中心位置），初始化的灰狼群体由n匹灰狼构成。

2）为每匹灰狼建立PLSR模型（即目标函数），并计算5折RMSECV值作为适应度值。其次，将最小、次小和第三小RMSECV值所对应的灰狼分别作为最优、次优和第三优解，记为α狼、β狼和δ狼。α狼、β狼和δ狼的位置将在迭代过程中更新。

3）模拟灰狼群体包围猎物的数学表达式如式（3）、（4）所示：

式中：t为当前迭代次数；Xp为猎物的位置；X为灰狼个体的位置；D为灰狼个体与猎物之间的距离；A、C为系数向量，其数学表达式分别如式（5）、（6）所示：

式中：a为在迭代过程中，由2线性减小到0；Imax为最大迭代次数。

4）由α狼、β狼和δ狼带领ω狼进行追捕猎物的数学表达式如式（7）～（13）所示：

式中：C1、C2和C3可由式（6）得到；A1、A2和A3可由式（5）得到；Xα、Xβ和Xδ分别为α狼、β狼和δ狼的位置；Dα、Dβ和Dδ分别为灰狼个体与α狼、β狼和δ狼之间的距离；X1、X2和X3为ω狼朝向α狼、β狼和δ狼前进的距离。

5）该算法持续迭代（重复执行步骤2～4）直到满足预先设置的最大迭代次数，最终，该算法输出最佳的光谱特征区间组合。

根据参考文献[24]，GWO-CMW算法的主要参数设置如下：种群规模n＝20，最大迭代次数Imax＝1 000。GWO-CMW算法的流程如图2所示。

图2 GWO-CMW算法流程图Fig.2 Flow chart of GWO-CMW algorithm

1.3.3.4 PSOGWO-CMW算法

将PSO融入GWO得到的混合智能群体算法，可以平衡PSO算法在全局搜索和局部搜索中的性能，同时实现灰狼自身经验的信息交换，完善位置更新策略[25]。本研究中，PSOGWO-CMW算法的基本思路为通过PSOGWO结合CMW筛选出最佳的光谱特征区间组合。PSOGWOCMW算法具体实施步骤与GWO-CMW算法基本一致，其区别在于将PSO融入GWO，完善了位置更新策略，具体表现在α狼、β狼和δ狼带领ω狼进行追捕猎物的数学表达式（式（14）～（21））：

式中：c1、c2、c3为加速度系数；ω为惯性权重；r3为[0,1]之间的随机数；V为灰狼个体的速度。

根据参考文献[26]，PSOGWO-CMW算法的主要参数设置如下：惯性权重ω＝0.5＋rand(·)/2（rand(·)为[0,1]之间的随机数），加速度系数c1＝c2＝c3＝0.5，种群规模n＝20，最大迭代次数Imax＝1 000。PSOGWO-CMW算法的流程如图3所示。

图3 PSOGWO-CMW算法流程图Fig.3 Flow chart of PSOGWO-CMW algorithm

1.3.4 模型评价指标

为了评估各模型的性能，一般采用校正集均方根误差（root mean squared error of calibration set，RMSEC），预测集均方根误差（root mean squared error of prediction set，RMSEP），校正集决定系数（coefficient of determination of calibration set，），预测集决定系数（coefficient of determination of prediction set，）和性能偏差比（ratio of performance to deviation，RPD）等模型评价指标，RMSE、R2和RPD的计算如式（22）～（24）所示：

式中：yi为真实值；为预测值；为真实值的平均值。当L表示校正集中样本的数量时，则式（22）对应为RMSEC，式（23）对应为；当L表示预测集中样本的数量时，则式（22）对应为RMSEP，式（23）对应为。本研究中，PLSR、PSO-CMW、GWO-CMW和PSOGWO-CMW算法由Matlab R2020a软件实现。

1.3.5 GC-MS检测CO-EVOO植物调和油中EVOO含量

为了验证本研究方法的准确性，采用标准方法（GCMS）检测CO-EVOO植物调和油中的EVOO含量，并与本研究方法的检测结果对比分析。GC-MS检测条件设置如下[27]。

1.3.5.1 色谱条件

色谱柱：HP-88毛细管柱（100 m×0.25 mm，0.20 μm）；升温程序：初始温度40 ℃，保持5.0 min，以5 ℃/min升温至245 ℃，保持5.0 min；进样口温度：250 ℃；进样模式：脉冲不分流进样，脉冲压力103.4 kPa，持续1.0 min；载气与流速：高纯氦气，1.0 mL/min。

1.3.5.2 质谱条件

传输线温度：280 ℃；电离模式：电子电离源；质量扫描范围：20～400 u；溶剂延迟时间：8.3 min。

1.4 数据处理

2 结果与分析

2.1 光谱预处理与光谱分析

光谱采集过程中带入的干扰信息（如噪声、基线漂移、光散射等）往往无法避免。因此，对光谱数据进行预处理有助于提高后续定性或定量分析的准确性。6 条代表性的CO-EVOO植物调和油原始拉曼光谱如图4A所示。本研究中，采用的光谱预处理策略主要包含以下3 个步骤：1）采用自适应迭代重加权惩罚最小二乘算法进行基线校正；2）采用多元散射校正算法进行光散射校正；3）采用卷积平滑算法进行光谱信号的平滑处理。预处理之后的拉曼光谱如图4B所示，对比图4A、B可以发现，预处理后拉曼光谱的基线漂移得到了抑制，光谱信号更加平滑，特征峰也愈加明显，为后续的定量分析奠定了良好的基础。由于CO-EVOO植物调和油中的主要成分为脂肪酸和甘油三酯，故主要的拉曼光谱特征峰及其振动归属为[28]：1 081 cm-1（C—C键拉伸）、1 260 cm-1（＝C—H键变形）、1 300 cm-1（C—H键变形）、1 434 cm-1（C—H键变形）、1 648 cm-1（C＝C键拉伸）和1 740 cm-1（C＝O键拉伸）。如图4C、D所示，随着EVOO含量的升高，特征峰1 260 cm-1和1 648 cm-1处的拉曼强度逐渐降低，这是能够定量分析CO-EVOO植物调和油中EVOO含量的重要因素。上述现象的主要原因在于CO-EVOO植物调和油中油酸与亚油酸比例的改变，油酸是EVOO中的主要脂肪酸，亚油酸是CO中的主要脂肪酸，随着油酸含量的升高以及亚油酸含量的降低，导致1 260 cm-1和1 648 cm-1处的拉曼强度降低[29]。最后，由于1 000～1 800 cm-1拉曼光谱区域信噪比高且包含了主要的特征峰，故本研究中拉曼光谱特征区间筛选及定量模型的构建均基于此区域。

图4 6 条代表性的CO-EVOO植物调和油的拉曼光谱Fig.4 Six representative Raman spectra of CO-EVOO blends

2.2 各模型窗口数量的优化

为了获得各模型最佳的窗口数量，将PSO-CMW、GWO-CMW和PSOGWO-CMW在CO-EVOO植物调和油拉曼光谱数据上各自独立运行30 次，并记录下RMSEP值进行对比分析。同时，考虑到计算量与模型的复杂度，窗口数量不宜过多，故将窗口数量的上限设置为100。图5所示为不同的窗口数量对各模型性能的影响。例如，当PSO-CMW的窗口数量为10时，RMSEP值很大，这是由于模型中包含的光谱变量过少，从而导致模型性能较差。当窗口数量由10增加到30时，RMSEP值明显下降。当窗口数量由40增加到100时，RMSEP值再次变大，这是由于模型中包含的光谱变量过多，从而导致模型出现了过拟合现象。因此，对于PSO-CMW而言，其最佳的窗口数量为30。同理，对于GWO-CMW和PSOGWO-CMW而言，其最佳的窗口数量均为40。

2.3 各模型检测结果

将CO-EVOO植物调和油拉曼光谱分别输入PLSR、PSO-CMW、GWO-CMW和PSOGWO-CMW模型筛选最佳的光谱区间组合，并在此基础上构建定量模型预测EVOO的含量。各模型建模过程中，最佳隐变量数（latent variables，LVs）由5折RMSECV值确定。各模型的检测结果如表1所示。

表1 各模型定量检测EVOO含量的结果Table 1 Performance parameters of different predictive models for EVOO content

2.3.1 PLSR模型

由表1可知，PLSR模型定量检测结果如下：RMSEC＝1.046 8，＝0.984 7，RMSEP＝1.847 2，＝0.952 2，RPD＝4.574 6。图6A为不同的LVs对应的RMSECV值，其中最小RMSECV值对应的最佳LVs为5。CO-EVOO植物调和油中EVOO含量PLSR预测值与真实值之间的关系如图6B所示。

图6 PLSR模型定量检测结果Fig.6 Quantitative results of PLSR model

2.3.2 PSO-CMW模型

由表1可知，最优的PSO-CMW模型定量检测结果如下：RMSEC＝0.647 2，＝0.990 1，RMSEP＝1.094 3，＝0.983 6，RPD＝7.806 4。PSO-CMW筛选的拉曼光谱特征区间如图7A所示，其中部分区间，如1 078～1 087、1 254～1 262、1 292～1 305、1 425～1 445、1 638～1 654 cm-1和1 735～1 748 cm-1分别覆盖了拉曼特征峰1 081、1 260、1 300、1 434、1 648 cm-1和1 740 cm-1。CO-EVOO植物调和油中EVOO含量PSOCMW预测值与真实值之间的关系如图7B所示。

图7 PSO-CMW模型定量检测结果Fig.7 Quantitative results of PSO-CMW model

2.3.3 GWO-CMW模型

由表1可知，最优的GWO-CMW模型定量检测结果如下：RMSEC＝0.553 1，＝0.992 0，RMSEP＝1.026 1，＝0.985 2，RPD＝8.216 9。GWO-CMW筛选的拉曼光谱特征区间如图8A所示，其中部分区间，如1 073～1 084、1 254～1 267、1 290～1 305、1 425～1 441、1 642～1 655 cm-1和1 733～1 742 cm-1分别覆盖了拉曼特征峰1 081、1 260、1 300、1 434、1 648 cm-1和1 740 cm-1。CO-EVOO植物调和油中EVOO含量GWOCMW预测值与真实值之间的关系如图8B所示。

图8 GWO-CMW模型定量检测结果Fig.8 Quantitative results of GWO-CMW model

2.3.4 PSOGWO-CMW模型

由表1可知，最优的PSOGWO-CMW模型定量检测结果如下：RMSEC＝0.499 2，＝0.993 7，RMSEP＝0.978 4，＝0.988 3，RPD＝9.242 1。PSOGWOCMW筛选的拉曼光谱特征区间如图9A所示，其中部分区间，如1 075～1 100、1 252～1 269、1 290～1 307、1 427～1 441、1 640～1 655 cm-1和1 731～1 748 cm-1分别覆盖了拉曼特征峰1 081、1 260、1 300、1 434、1 648、1 740 cm-1。CO-EVOO植物调和油中EVOO含量PSOGWOCMW预测值与真实值之间的关系如图9B所示。

图9 PSOGWO-CMW模型定量检测结果Fig.9 Quantitative results of PSOGWO-CMW model

2.4 各模型检测结果对比分析

各模型评价指标的变化趋势如图10A所示，显然，PLSR模型的性能较差，主要在于PLSR模型中包含了所有的光谱变量，其中既包含有信息的变量，也包含无信息的变量甚至干扰性的变量。相较于PLSR模型，PSOCMW、GWO-CMW和PSOGWO-CMW通过筛选光谱特征区间，各模型的预测性能均获得了提高。以PLSR模型的预测性能作为基准，PSO-CMW、GWO-CMW和PSOGWO-CMW模型的RMSEC分别降低了38.2%、47.2%和52.3%；PSO-CMW、GWO-CMW和PSOGWO-CMW模型的分别提高了0.5%、0.7%和0.9%；PSO-CMW、GWO-CMW和PSOGWO-CMW模型的RMSEP分别降低了40.8%、44.5%和47.0%；PSO-CMW、GWO-CMW和PSOGWO-CMW模型的分别提高了3.3%、3.5%和3.8%；PSO-CMW、GWO-CMW和PSOGWO-CMW模型的RPD分别提高了70.6%、79.6%和102.0%。各模型迭代收敛曲线如图10B所示，显然，PSOGWOC M W 模型的收敛速度更快，并拥有更好的寻优性能。综上所述，相较于PSO-CMW和GWO-CMW模型，PSOGWO-CMW模型具有最优的性能，主要在于PSOGWO充分利用了PSO与GWO各自的优势，有效地平衡了局部搜索和全局开发的能力，从而提升了模型的整体性能。

图10 各模型检测性能对比Fig.10 Comparisons of prediction performance of different models

2.5 真实样本检测分析

3 种品牌的植物调和油：金龙鱼（V（C O）∶V（EVOO）＝90∶10）、长寿花（V（CO）∶V（EVOO）＝94∶6）和鑫榄源（V（CO）∶V（EVOO）＝90∶10）分别购于当地的永辉、大润发和雅家乐超市。每种品牌的植物调和油分别准备10 个样本，每个样本为10 mL。每个样本分别采用本方法与标准方法（GC-MS）[27]检测EVOO含量，检测结果如表2所示。将两种方法的检测结果做双侧配对t检验，结果表明两者无显著性差异（P＝0.38＞0.05）。根据公式检测限＝3S0/K（S0为多个空白样本响应值标准差，K为校正曲线的斜率）[30]，可估算得到本方法对EVOO含量的检测限为1.25%。由于市场上植物调和油中高价值植物油的含量一般大于等于5%，故本方法可实现对真实植物调和油样本品质的定量鉴别。

表2 本方法与标准方法检测真实植物调和油样本中EVOO含量的结果Table 2 Comparison of the results of the proposed method and the standard method for EVOO contents in real BEVO samples

3 结论

本研究提出了一种基于拉曼光谱与PSOGWO-CMW算法实现植物调和油中高价值植物油含量快速定量检测的方法。以配制的CO-EVOO植物调和油为检测对象，与PLSR、PSO-CWM和GWO-CMW模型相比，PSOGWO-CMW模型具有最佳的建模性能，RMSEC＝0.499 2，＝0.993 7，RMSEP＝0.978 4，＝0.988 3，RPD＝9.242 1。同时，将本方法与标准方法分别检测真实的CO-EVOO植物调和油样本中EVOO含量，结果表明两者的检测性能无显著性差异（P＝0.38＞0.05）。综上所述，通过本研究结果，验证了本方法快速定量检测CO-EVOO植物调和油中EVOO含量的有效性与可行性。同时，本方法也适用于其他植物调和油中高价值植物油含量的快速定量检测。