CN108802002B - 一种快速无损鉴别解除滞育的蚕卵拉曼光谱模型构建方法 - Google Patents

Abstract

本发明属于拉曼光谱检测技术领域，本发明公开了一种快速无损鉴别解除滞育的蚕卵拉曼光谱模型构建方法。该方法通过拉曼光谱检测技术和PLS_DA分类算法，主要包括甄选样品、数据采集、数据挖掘、结论分析等。所述甄选样品包括甄选大造和9芙×7湘两个品种滞育卵与解除滞育卵的筛选，所述数据采集通过微区激光拉曼光谱仪进行光谱数据采集，并对拉曼光谱数据的进行预处理，进行PLS_DA定性分析拉曼光谱模型的构建，该方法能快速无损地识别出滞育或解除滞育的蚕卵，最终实现了快速无损家蚕卵解除滞育的鉴别分析。

Description

一种快速无损鉴别解除滞育的蚕卵拉曼光谱模型构建方法

技术领域

本发明属于拉曼光谱检测技术领域，更具体地，涉及一种快速无损鉴别解除滞育蚕卵的拉曼光谱模型构建方法。

背景技术

拉曼光谱分析法是基于分子振动、转动方面信息，并应用于分子结构研究的一种分析方法。对于各种物质的化学基团，拉曼光谱具有尖锐的特征谱峰，谱峰的位置和强弱可以灵敏地反映出有关物质的结构及其变化信息。因此，拉曼光谱分析法是一种对物质进行定性鉴别与定量分析的有效手段。

滞育性家蚕卵的物质构成及其构成比例都和解除滞育的家蚕卵有着非常大的不同。家蚕卵的主要物质是水65％，蛋白质20％，脂质10％，炭水化物4％，无机盐类等。昆虫进入滞育前都有一个积蓄营养物质和能量的过程。家蚕在蛹期卵巢发育过程中积蓄大量的糖原，使得刚产下的滞育卵的糖原含量就高于解除滞育的蚕卵。同样，脂类也是昆虫体内重要的能源和营养物质，在滞育卵内脂质的含量远高于糖原，滞育卵脂质含量高于解除滞育的蚕卵。

综上所述，滞育性卵中的糖原和脂质在含量上均远远大于解除滞育蚕卵，且构成脂质的成分也不同。从而导致两者的拉曼光谱图的谱峰位置和峰强都有差别，经过计算机的数据处理，分类算法可以很好地将滞育卵、解除滞育卵分辨出来。利用拉曼光谱检测技术来分辨滞育的蚕卵和解除滞育的蚕卵，有理论基础和技术支持，具有可行性。

家蚕是一种以卵繁殖的卵滞育昆虫，一般在卵产下后25℃的保护下，约经过7天蚕卵即进入到完全滞育状态，滞育卵须在低温中保存数日并结合相应的刺激才能解除滞育并孵化，从而满足生产养蚕用种需要。滞育卵与解除滞育卵在物质构成及构成比例上有着明显差异，但滞育卵与解除滞育卵在外观上用肉眼难以区别，目前生产上主要是运用蚕卵胚子解剖方法，并结合显微镜观察而作出判断，但解剖过的蚕卵不能再孵化养蚕，且需要技术人员操作熟练，掌握胚子解剖关键技术要领，视力要好。

在目前的家蚕行业上，尚未发现有人将拉曼光谱检测技术应用于家蚕滞育性的区分上。

发明内容

为了解决上述现有技术存在的不足和缺点，提供一种快速无损鉴别解除滞育的蚕卵拉曼光谱模型构建方法。该检测方法将分类算法与拉曼光谱进行有机融合，以期能得到较好的应用前景。

本发明的目的通过下述技术方案来实现：

一种快速无损鉴别解除滞育的蚕卵的拉曼光谱模型构建方法，包括如下具体步骤：

S1.拉曼光谱数据的采集：取不同样品蚕卵的滞育卵与滞育解除卵各100粒，将拉曼光谱仪的入射光对准蚕卵中部进行检测，得到清晰的图像后，采集“波数-光谱强度”的拉曼光谱数据；

S2.拉曼光谱数据的预处理Ⅰ：将步骤S1中采集的“波数-光谱强度”的拉曼光谱数据进行插值处理，提取拉曼强度在800～1800cm^-1波数范围内的所有整数点上的近似值，使每个样本的波数取值相同；

S3.拉曼光谱数据的预处理Ⅱ：对插值后样本数据依次进行窗口移动多项式最小二乘平滑和自适应迭代惩罚最小二乘法去除荧光背景；一阶求导操作以除去存在噪声和荧光背景的干扰，提高分类效果准确率；

S4.拉曼光谱模型构建：用PLS_DA算法对预处理后样本数据进行拉曼光谱模型的构建；

S5.对测试集进行预测：将已测滞育和解除滞育的蚕卵的拉曼光谱剩余的数据作为测试集，用所建拉曼光谱模型对测试集进行预测；

S6.求取滞育蚕卵和解除滞育蚕卵拉曼光谱的均值的差，再乘以对应的回归系数矩阵B，可得到上述两种蚕卵光谱在PLS_DA模型里的数学距离ΔY，即有：ΔY＝B×(X_滞育-X_解除滞育)，其中，X_滞育是滞育蚕卵的光谱强度，X_解除滞育是解除滞育蚕卵的光谱强度，B是PLS_DA算法构建的建模的回归系数矩阵，提取出ΔY的绝对值大于0.05且波数在拉曼谱峰附近所对应的波数，再次利用PLS_DA进行数学建模，结合建模效果和拉曼谱峰的位置，可以提取家蚕卵的拉曼光谱的特征峰，即可鉴别解除滞育的蚕卵。

优选地，步骤S1中所述的拉曼光谱仪的激光波长设为780～790nm，激光功率40mw；中心波长设为1400cm^-1；光谱检测最低波数为400cm^-1，最高波数为1800cm^-1；积分时间为10秒，累积次数为1次，使用石英显微镜头镜聚焦。

优选地，步骤S2中所述自适应迭代惩罚最小二乘法的惩罚因子为5000。

本发明中采用的偏最小二乘判别分析(partial least squares discriminationanalysis，PLS_DA)，是一种用于判别分析的多变量统计分析方法。PLS-DA是基于为了分类目的而搜索最优潜变量数据集的监督方法。优点是适用于对象数量少于变量数量的情况下。它融合了偏最小二乘和判别分析。偏最小二乘，是一种数学优化技术，它通过最小化误差的平方和找到一组数据的最佳函数匹配。用最简的方法求得一些绝对不可知的真值，而令误差平方之和为最小，它主要是为了解决回归中的共线性问题，在光谱分析中得到了广泛的应用。判别分析是一种根据观察或测量到的若干变量值，来判断研究对象如何分类的常用统计分析方法。其原理是对不同处理样本(如观测样本、对照样本)的特性分别进行训练，产生训练集，并检验训练集的可信度。

作为一个多元线性回归方法，偏最小二乘判别分析的主要目的是要建立一个线性模型：Y＝XB+E，其中Y是具有m个变量、n个样本点的输出矩阵，X是具有p个变量、n个样本点的输入矩阵，B是回归系数矩阵，E为噪音校正模型，与Y具有相同的维数。

求导可以消除基线漂移、强化谱带特征、克服谱带重叠，是常见的光谱预处理方法。一阶微分可以去除和波长无关的基线漂移；二阶微分可以去除同波长线性相关的漂移。但求导在消除基线漂移的过程中会放大噪声，因此求导前需对原始谱图进行平滑处理，常釆用卷积求导法。一阶求导后处理的目的在于锐化光谱特征，并减少光谱的偏移，以便后续的算法尽可能地起到更好的分辨效果。当然，与此同时不可避免地放大了噪声，但这仍在可接受的范围。

较之可见光和红外线等常用光谱，拉曼效应十分微弱，加上检测器的背景噪声、样品放置方向及其他因素的影响，获得的拉曼光谱信号往往叠加了各种噪声。这些噪声不仅影响光谱有用信息的获得，同时影响后续分类模型的建立以及对样品的预测效果。主要对影响较大的背景噪声进行处理。而激光引起的荧光又是拉曼光谱技术中最普遍最主要的背景噪声来源。荧光光谱外观通常比拉曼光谱峰要宽得多，看起来像拉曼光谱缓慢变化的基线，可以使用平滑、求导、曲线拟合、小波变换、归一化等的算法给予去除。

平滑：平滑可以提高分析信号的信噪比，能有效地消除谱图中的高频成分而保留低频成分，最简单的平滑方法是窗口移动平均法，但常用的窗口移动多项式最小二乘平滑(Savitzky-Golay Smoothing)。其优点是算法相对简单，计算速度较快，缺点是窗口过大可能会失去有用的高频信号。

求导：求导可以消除基线漂移、强化谱带特征、克服谱带重叠，是常见的光谱预处理方法。一阶微分可以去除和波长无关的基线漂移；二阶微分可以去除同波长线性相关的漂移。但求导在消除基线漂移的过程中会放大噪声，因此求导前需对原始谱图进行平滑处理，常釆用卷积求导法。

自适应迭代惩罚最小二乘法：它可以非常快速地完成大量光谱数据的背景扣除，且对于低信噪比的情况下可以保持一个有效的信号。它的做法是：利用迭代重加权惩罚最小二乘算法，逐步逼近背景，引入参数来调整曲线的平滑程度，将原始信号减去匹配的背景，完成背景扣除。

小波变换：小波变换是空间(时间)和频率的局部变换，因而能有效地从信号中提取信息。通过伸缩和平移等运算功能可对函数或信号进行多尺度的细化分析。小波变换可以有效消除背景和噪声。小波变换也被成功运用红外光谱、紫外光谱、拉曼光谱等的预处理中，由于拉曼光谱是由尖锐的拉曼信号叠加在平缓的荧光光谱上组成的，因此利用小波变换可以检测到拉曼光谱中尖锐的拉曼谱峰信号。尽管小波变换用于光谱预处理取得了很好的效果，但是对小波变换中某些参数的选取，还没有一定的标准，参数设置的正确与否，直接影响到预处理的效果，进而影响到校正模型的预测精度。

归一化：目前使用的归一化方法有外标法和内标法。外标法在测量样品的同时通过分光检测激光强度由此可以克服激光器功率波动的影响此方法要求采样位置固定，控制积分时间，同时需要使用多路光谱仪，适用于在线分析。使用外标法的另一个优势是可以获得激光器波长漂移的信息，有利于在线仪器的连续稳定运行。但由于荧光物质会产生强烈的荧光，从而抑制拉曼信号，使得通过外标法无法获得这些信息。因此，外标法通常只适用于荧光较小的试样。内标法是指利用光谱内部特定谱峰的强度，进行归一化。在复杂有机物分析中，可以将处的拉曼峰作为参考峰来进行归一化。

与现有技术相比，本发明具有以下有益效果：

1.本发明的拉曼光谱模型构建方法无需复杂的解剖操作，可快速无损鉴别解除滞育蚕卵，对操作者的技术熟练程度的要求大大降低，可依据用需要，灵活快速鉴别蚕卵解除滞育与否，其准确率为高，可达95％以上。

2.本发明将分类算法与拉曼光谱进行有机融合，以期能得到较好的应用前景。

附图说明

图1为滞育蚕卵和解除滞育的蚕卵的照片。

图2为大造和9芙×7湘品种滞育与解除滞育蚕卵样本拉曼光谱图。

图3为拉曼光谱数据的预处理时进行窗口移动多项式最小二乘平滑前后的处理效果的对比图。

图4为拉曼光谱数据的预处理时进行自适应迭代惩罚最小二乘法去除荧光背景前后的处理效果对比图。

图5为拉曼光谱数据的预处理时进行一阶求导处理前后的效果对比图。

图6为大造和9芙×7湘各品种滞育与解除滞育蚕卵的平均拉曼峰。

图7为对滞育和解除滞育蚕卵(大造滞育、大造非滞育、9芙×7湘滞育、9芙×7湘非滞育)求得的波数-ΔY折线图。

具体实施方式

下面结合具体实施例进一步说明本发明的内容，但不应理解为对本发明的限制。若未特别指明，实施例中所用的技术手段为本领域技术人员所熟知的常规手段。除非特别说明，本发明采用的试剂、方法和设备为本技术领域常规试剂、方法和设备。

本发明采用的实验样品是广东常见蚕种，由动物与科学学院的王叶元副教授培育家蚕卵(包括滞育卵与解除滞育卵)作为实验样品，人工随机挑选出大造品种滞育蚕卵和解除滞育蚕卵各100枚，9芙×7湘品种滞育蚕卵和解除滞育蚕卵各100枚进行实验。实验样本如图1所示，图1为滞育蚕卵和解除滞育的蚕卵的照片，其中：(a)为解除滞育蚕卵，(b)为滞育蚕卵。从图1中可以观察到滞育与解除滞育蚕卵为芝麻大小，仅从形态和颜色上，肉眼不能看出二者明显的区别。

拉曼光谱仪采用北京卓立汉光仪器有限公司生产的FinderOne激光共聚焦显微拉曼光谱仪，其采用优化的光学显微成像光路设计，搭配高性能影像校正光谱仪和LDC-DD技术的CCD探测器。激光器激发波长可选：325nm、532nm、785nm等。其光谱检测范围为：100^-3至3000cm^-1，光谱分辨率：3cm^-1。探测器：TE深制冷型背感光CCD阵列探测器。

实施例1家蚕卵拉曼光谱模型的构建

1.拉曼光谱数据的采集：拉曼光谱仪的激光波长设为785nm，激光功率40mw；使用FinderOne的连续检测功能，中心波长设为1400cm^-1；光谱检测最低波数为400cm^-1，最高波数为1800cm^-1；积分时间为10秒，累积次数为1次；使用石英显微镜头镜聚焦；测量室温为20℃。每次检测时将入射光对准蚕卵中部，调整焦距，使之得到清晰的图像，再采集得到“波数-光谱强度”格式的拉曼光谱数据。图2为大造和9芙×7湘品种滞育与解除滞育蚕卵样本拉曼光谱图，其中，(a)为大造滞育，(b)为9芙×7湘滞育，(c)为大造解除滞育，(d)为9芙×7湘解除滞育。从图2中可以看出，大造和9芙×7湘的拉曼强度平均值相差不大，但大造的标准差要明显大于9芙×7湘。蚕卵的拉曼信号特征峰比较明显，但存在噪声和荧光背景的干扰，故需对原始拉曼谱图进行数学方法的预处理以及算法的分析。

2.拉曼光谱数据的预处理Ⅰ：由于使用的拉曼光谱仪所采集的光谱数据自变量(波数)并不刚刚在规律的整数点上，且在每次采集时都有些许的偏移，这必然对数据分析造成了负面影响。为使各样本的数据点的自变量相同，所有光谱数据首先被进行了插值处理。

3.拉曼光谱数据的预处理Ⅱ：再进行窗口移动多项式最小二乘平滑、自适应迭代惩罚最小二乘方法、一阶求导的算法等预处理，以除去存在噪声和荧光背景的干扰，提高分类效果准确率。

4.拉曼光谱模型构建：取预处理后样本数据的60％作为训练集，利用计算机分析软件Matlab2014b的强大数据处理功能，进行偏最小二乘判别分析(PLS_DA)算法的拉曼光谱模型的构建。

5.家蚕卵拉曼特征峰的提取：求取滞育蚕卵和解除滞育蚕卵拉曼光谱的均值的差，再乘以对应的回归系数矩阵B，可得到上述两种蚕卵光谱在PLS_DA模型里的数学距离ΔY，即有：ΔY＝B×(X_滞育-X_解除滞育)，其中，X_滞育是滞育蚕卵的光谱强度，X_解除滞育是解除滞育蚕卵的光谱强度，B是PLS_DA算法构建的建模的回归系数矩阵，提取出ΔY的绝对值大于0.05且波数在拉曼谱峰附近所对应的波数，再次利用PLS_DA进行数学建模，结合建模效果和拉曼谱峰的位置，可以提取家蚕卵的拉曼光谱的特征峰，即可鉴别解除滞育的蚕卵。

如果ΔY越大，则说明对应的波数更为重要。当提取ΔY的绝对值中最大的某几项，再次利用PLS_DA进行数学建模，可以简化模型，提高预测稳定性。根据PLS_DA的原理，其回归系数矩阵中，波数对应的系数越大，则说明该波数在建模中的重要性越大。结合光谱数据的分布特点，可以对家蚕卵的拉曼光谱进行特征峰的识别。

实施例2

基于PLS_DA算法的回归系数矩阵的家蚕滞育性分析：根据PLS_DA的原理，其回归系数矩阵中，波数对应的系数越大，则说明该波数在建模中的重要性越大。结合光谱数据的分布特点，可以对家蚕卵的拉曼光谱进行特征峰的识别。操作如下：在去除荧光之后，求取滞育蚕卵和解除滞育蚕卵拉曼光谱的均值的差，再乘以对应的回归系数矩阵B，可得到两种蚕卵光谱在PLS_DA模型里的数学距离ΔY，即有：ΔY＝B×(X_滞育-X_解除滞育)。

图3为拉曼光谱数据的预处理时进行窗口移动多项式最小二乘平滑前后的处理效果的对比图。窗长为11，从图3可以看到，该平滑方式能较好地去除高频噪声，使光谱曲线更平滑，同时保留光谱的特征信息，为后续的处理扫除障碍。

图4为拉曼光谱数据的预处理时进行自适应迭代惩罚最小二乘方法去除荧光背景前后的处理效果对比图。从图4中可知，其中惩罚因子为5000，该方法拟合出了背景荧光，其提取出的拉曼谱峰较明显，适合做进一步的处理。图5为拉曼光谱数据的预处理时进行一阶求导处理前后的效果对比图。从图5中可知，一阶求导后分辨效果并不理想，不可避免地放大了噪声，但这仍在可接受的范围。其目的在于锐化光谱特征，并减少光谱的偏移，以便后续的算法尽可能地起到更好的分辨效果。图6为大造和9芙×7湘各品种滞育与解除滞育蚕卵的平均拉曼峰。从图6中看出，四种蚕卵(大造滞育、大造非滞育、9芙×7湘滞育、9芙×7湘非滞育)1274、1444、1506、1605、1671cm^-1等波数的拉曼特征峰较明显，其中在1274、1506cm^-1波数处，滞育的拉曼峰都要明显比解除滞育的拉曼峰高。在1562cm^-1波数处，解除滞育的蚕卵有拉曼峰，而滞育的蚕卵在1573cm^-1波数处形成拉曼峰，二者有11个波数的偏移。在800-1800cm^-1波数范围内，主要是C-H面外弯曲振动和C-N的伸缩振动所引起的特征峰。在1400-1800cm^-1波数，该范围内的拉曼谱峰主要由C＝C伸缩振动引起。其中，1274cm^-1处是C-N的伸缩振动和C-H面内弯曲振动的耦合；1444cm^-1附近代表的是包涵体蛋白，1506cm^-1附近代表的是蛋白质中酰胺Ⅱ谱线，1671cm^-1是蛋白质的酰胺Ⅰ谱线，1605cm^-1附近是邻苯二甲酸酯类。说明1274、1444、1506、1605、1671cm^-1波数对应物质在滞育与解除滞育蚕卵存在差异。

实施例3家蚕卵拉曼光谱模型性能的评判

使用混淆矩阵及其衍生指标和受试者工作特征曲线等进行模型好坏的界定。混淆矩阵(Confusion Matrix)用于展示模型预测精度的一种较为常见的方法。它是一个N*N列联表(N为分类的类别数)。以二分类为例，分类结果的混淆矩阵如下表1所示。

表1混淆矩阵

根据表1中的四个基本指标可以衍生出多个分类器评价指标，常用指标如下：

准确率(accuracy)：它针对的是测试集中所有识别正确的数量。计算公式为：accuracy＝(TP+TN)/(TP+TN+FN+FP)。

真正率(true positive rate，TPR)，它表现出在所有真正是正样本中，分类模型能召回数量的多少。计算公式为：TPR＝TP/(TP+FN)。

真负率(true negative rate，TNR)，与真正率相反，它表现出在所有真正是负样本中，分类模型能正确识别出的数量的多少。计算公式为：TNR＝TN/(FP+TN)。

假正率(False positive rate，FPR)，它表现出被错预测为正例的真实负例占所有真实负例的比例，它与真负率之和为1。计算公式为FPR＝FP/(FP+TN)；

其中，真正例(True Positive，TP)：指被分类器正确分类的正元组。令TP为真正例的个数。真负例(True Negative，TN)：指被分类器正确分类的负元组。令TN为真负例的个数。假正例(False Positive，FP)：指被分类器错误标记为正元组的负元组。令FP为假正例的个数。假负例(False Negative，FN)：指被分类器错误标记为负元组的正元组，令FN为假负例的个数。

表2样本数据分组

PLS_DA模型的预测结果的混淆矩阵分别如表3所示。

表3PLS_DA预测结果的混淆矩阵

为客观评价各个模型，根据各个混淆矩阵，计算出各个模型的准确率、真负率和真正率等参数，如表4所示，可以得出以下结果：若从不同品种上看，三种品种的各项参数相差不大，以PLS_DA的准确率为例，大造品种最高，达到了0.987，而9芙×7湘和混合模型为0.9750，稍逊于大造。

表4两种模型的评价参数

以混合品种的PLS_DA模型为对象，绘制“波数-ΔY”的折线图，如图7所示。图7为对滞育和解除滞育蚕卵(大造滞育、大造非滞育、9芙×7湘滞育、9芙×7湘非滞育)求得的波数-ΔY折线图。从图7中可知，提取出ΔY的绝对值大于0.05所对应的波数，即：1274、1283、1437、1516、1585cm^-1共五个波数。但与家蚕的拉曼谱峰(图6)进行对比后，剔除1283cm^-1和1585cm^-1并不在拉曼谱峰附近的拉曼谱峰。再次提取出剩下的3个波数，仍以原来的测试集和预测集进行PLS_DA建模和预测验证，准确率为0.9625。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合和简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (3)

1.一种快速无损鉴别解除滞育的蚕卵的拉曼光谱模型构建方法，其特征在于，包括如下具体步骤：

S1.拉曼光谱数据的采集：取不同样品蚕卵的滞育卵与滞育解除卵各100粒，将拉曼光谱仪的入射光对准蚕卵中部进行检测，得到清晰的图像后，采集“波数-光谱强度”的拉曼光谱数据；

S2.拉曼光谱数据的预处理Ⅰ：将步骤S1中采集的“波数-光谱强度”的拉曼光谱数据进行插值处理，提取拉曼强度在800～1800cm^-1波数范围内的所有整数点上的近似值，使每个样本的波数取值相同；

S3.拉曼光谱数据的预处理Ⅱ：对插值后样本数据依次进行窗口移动多项式最小二乘平滑和自适应迭代惩罚最小二乘法去除荧光背景；一阶求导操作以除去存在噪声和荧光背景的干扰，提高分类效果准确率；

S4.拉曼光谱模型构建：用PLS_DA算法对预处理后样本数据进行拉曼光谱模型的构建；

S5.对测试集进行预测：将已测滞育和解除滞育的蚕卵的拉曼光谱剩余的数据作为测试集，用所建拉曼光谱模型对测试集进行预测；

S6.求取滞育蚕卵和解除滞育蚕卵拉曼光谱的均值的差，再乘以对应的回归系数矩阵B，可得到上述两种蚕卵光谱在PLS_DA模型里的数学距离ΔY，即有：ΔY＝B×(X_滞育-X_解除滞育)，其中，X_滞育是滞育蚕卵的光谱强度，X_解除滞育是解除滞育蚕卵的光谱强度，B是PLS_DA算法构建的建模的回归系数矩阵，提取出ΔY的绝对值大于0.05且波数在拉曼谱峰附近所对应的波数，再次利用PLS_DA进行数学建模，结合建模效果和拉曼谱峰的位置，可以提取家蚕卵的拉曼光谱的特征峰，即可鉴别解除滞育的蚕卵。

2.根据权利要求1所述的快速无损鉴别解除滞育蚕卵的拉曼光谱模型构建方法，其特征在于，步骤S1中所述的拉曼光谱仪的激光波长设为780～790nm，激光功率40mw；中心波长设为1400cm^-1；光谱检测最低波数为400cm^-1，最高波数为1800cm^-1；积分时间为10秒，累积次数为1次，使用石英显微镜头聚焦。

3.根据权利要求1所述的快速无损鉴别解除滞育蚕卵的拉曼光谱模型构建方法，其特征在于，步骤S2中所述自适应迭代惩罚最小二乘法的惩罚因子为5000。

Images