CN105115938A - 一种饲料原料中三聚氰胺掺假定量检测方法 - Google Patents

Abstract

本发明公开了一种饲料原料中三聚氰胺掺假定量检测方法，通过采用成像光谱仪获取训练样本的近红外成像光谱，对光谱进行预处理；采用偏最小二乘判别分析(PLS_DA)方法得到训练样本的像元点个数比值，建立其与样本浓度之间的线性关系，得到分界点浓度值；波长变量选择后，采用平均光谱分阶段建立偏最小二乘(PLS)定量模型；然后利用光谱仪采集未知样本的光谱，经相同的光谱预处理及选择相同的波长变量，通过已建PLS_DA判别模型和PLS定量模型预测该未知样本的三聚氰胺浓度。采用本发明提供的方法，能够有效解决现有掺假检测方法的检测限高、难于实现痕量检测、光谱数据相关性高、信息冗余的问题，并能显著提高三聚氰胺浓度的预测精度及稳定性，降低检测限。

Description

一种饲料原料中三聚氰胺掺假定量检测方法

技术领域

本发明涉及农业中饲料掺假定量检测分析技术领域，尤其涉及一种饲料原料中三聚氰胺掺假定量检测方法。

背景技术

当前，我国作为世界上最大的畜禽养殖生产国之一，畜禽食品安全现状较为严峻，食品安全事件层出不穷。由于优质畜禽产品需要高效、安全的饲料原料作保障，因此，饲料原料安全成为畜禽产品安全的前提。蛋白类饲料原料蛋白质含量丰富，作为动物的主要营养来源，供应紧张且价格昂贵，导致一些不法企业将非蛋白氮类化合物三聚氰胺掺入其中，以提高粗蛋白含量，以次充好，降低成本，谋取利益。这种掺假饲料原料，不仅会对动物生长产生直接危害，而且给消费者的健康带来潜在威胁。采用传统的化学分析方法进行掺假检测时，需要繁琐复杂前处理操作过程，耗时耗力，成本高，人为因素影响大，并且对检测样品具有破坏性，无法实现实时检测，严重限制了所能够分析的样品数量，对于饲料产品质量检验工作来说这一缺点尤为显著。

光谱作为一个重要的光学分析手段，包含物质特有的分子振动信息，可用于分析测定与这些信息相关的化学性质和化学成分含量，具有快速准确、无损和实时检测的特点。目前对饲料原料中三聚氰胺掺假检测采用的最广泛的方法为近红外光谱法。然而，饲料原料掺假样本中存在三聚氰胺含量较低的现象，检测存在风险，其检测误差相对较大，样本的不均匀性和异质性及样品成分稀释的问题，严重影响了传统近红外光谱进行饲料原料掺假定量检测分析的精度。近几年来，新型的近红外显微成像光谱技术以其光谱和空间分辨率高，谱图合一的特点，克服了样本不均匀性的特点，但其分辨率的大幅度提高，大大的增加了光谱数据的维数，使得计算复杂，预测精度和速度降低。

目前所用的近红外显微成像光谱分析法通常为采集掺假饲料原料样本的微米级的近红外光谱曲线，存在采集的低浓度饲料原料掺假样本的光谱曲线中几乎不含三聚氰胺的光谱信息，或者掺假样本中大量饲料原料基质的光谱信息严重覆盖了其中三聚氰胺特征的现象，并且数据量大幅度的提高，增加了提取有效特征的难度，从而使得检测限较高，建立的定量模型预测结果不准确，误差大，定量检测精度低。

综上所述，利用常用的近红外光谱检测技术，存在检测限高、检测误差大、难以实现痕量检测的问题；同时，显微成像光谱技术存在光谱数据维数高、冗余信息多、高相关性的缺陷。

发明内容

本发明技术解决问题：克服现有技术的不足，提供一种简单、快速、高效的饲料原料中三聚氰胺掺假定量检测方法，以解决现有技术中难以实现痕量检测、定量检测精度低及成像光谱数据相关性高、信息冗余的问题。

本发明提供的饲料原料中三聚氰胺掺假定量检测方法以纯饲料原料近红外光谱为背景，采集掺假样本的近红外成像光谱数据，利用近红外成像光谱数据中像元点个数比值与样本三聚氰胺浓度之间的相关关系，并且采用波长变量选择方法降低成像光谱数据维数，分阶段建立掺假定量检测模型，使得检测限降低，实现痕量检测。

为达到上述目的，本发明的技术方案是这样实现的：

步骤101，以纯饲料原料的近红外光谱为背景，采用近红外显微成像光谱仪采集不同浓度掺假样本的近红外成像光谱数据，构成训练样本集；

步骤102，对所述近红外成像光谱数据进行预处理，以去除噪声、基线的干扰；

步骤103，采用偏最小二乘判别分析(PLS_DA)方法得到所述训练样本的近红外成像光谱数据中属于三聚氰胺的像元点，得到三聚氰胺像元点个数与总像元点个数的比值；

步骤104，建立不同浓度训练样本的像元点个数比值与三聚氰胺浓度之间的线性关系，得到分阶段定量模型的分界点浓度值；

步骤105，对训练样本的近红外成像光谱数据进行变量选择，根据所述的分界点浓度值，采用训练样本的平均光谱数据分阶段建立高浓度和低浓度的偏最小二乘(PLS)定量模型；

步骤106，用建立的PLS_DA模型判别未知样本的近红外成像光谱数据中属于三聚氰胺的像元点，得到像元点个数比值；

步骤107，根据所述的像元点个数比值确定定量模型阶段，对未知样本的近红外成像光谱数据进行变量选择后，取各个未知样本的平均光谱数据，用所构造的偏最小二乘定量模型，检测未知样本的三聚氰胺浓度。

其中，所述述步骤102中对所述近红外成像光谱数据进行预处理，具体包括SG平滑过程与自适应迭代权重-偏最小二乘方法进行基线校正。

其中，所述步骤103中采用PLS_DA方法判别训练样本的近红外成像光谱中属于三聚氰胺的像元点，获得三聚氰胺像元点个数与总像元点个数的比值，具体如下：

设已经存在训练样本集，包含成像光谱数据X(n×l)和类别矩阵Y(n×m)，其中n为样本，l为波长变量，m为类别，存在关系式记为：

Y＝BX+E

其中B为回归系数；E为噪声。训练样本集中，Y矩阵的类别值分别用“1”和“0”表示，其中1表示该像元为纯饲料原料像元点，0表示该像元为三聚氰胺像元点。选取纯饲料原料与三聚氰胺的平均光谱数据作为训练集，建立PLS-DA模型训练已知的纯饲料原料样本光谱数据，设定分类阈值使得纯饲料原料样本成像光谱数据中各像元点的类别均为1，将该阈值作为最终掺假鉴别中PLS_DA模型的分类阈值，进而逐一判别掺假样本的近红外成像光谱数据中属于三聚氰胺的像元点。

得到训练样本集的像元点个数比值，公式如下：

$p=\frac{{\stackrel{\‾}{x}}_{\begin{array}{cc}mela\min e& pixels\end{array}}}{{\stackrel{\‾}{X}}_{\begin{array}{cc}all& pixels\end{array}}}$

其中p为像元点个数比值；为训练样本的成像光谱中属于三聚氰胺的像元点个数；为所述样本的总像元点个数。

其中，所述步骤104中建立所述训练样本集的像元点个数比值与三聚氰胺浓度的线性关系，获得分阶段定量模型的分界点浓度值，具体如下：

将不同浓度训练样本的像元点个数比值与相应的样本中三聚氰胺浓度进行线性拟合，同时得到线性关系式与相关系数值，公式如下：

p＝aC+b

$R=\frac{\underset{i=1}{\overset{n}{\Σ}}(p_i-\stackrel{\‾}{p})(C_i-\stackrel{\‾}{C})}{\sqrt{\underset{i=1}{\overset{n}{\Σ}}{(p_i-\stackrel{\‾}{p})}^2\underset{i=1}{\overset{n}{\Σ}}{(C_i-\stackrel{\‾}{C})}^2}}$

其中p为像元点个数比值；C为样本中三聚氰胺的浓度；a和b为拟合系数；R为相关系数；n为样本个数；i＝1,2,…,n；p_i为第i个样本的像元点个数比值；为像元点个数比值平均值；C_i为第i个样本的三聚氰胺浓度；为样本浓度平均值。

观察所述拟合曲线，局部选取类似分界点浓度值，分别进行分阶段线性拟合，比较不同分界点浓度值构成的分阶段拟合曲线的相关系数值，选择相关系数最好的曲线相对应的分界点浓度值作为分阶段定量模型的分界点。

其中，所述步骤105中对训练样本的成像光谱数据进行变量选择，进一步包括最小角回归-最小绝对收缩选择(Lars_Lasso)方法，具体如下：

Lasso方法的算法原理为：

设有数据(Xⁱ,y_i),i＝1,2,...,n，n为样本个数，其中的X和y分别为第i个观测值的自变量和应变量，假设各观测值是彼此独立且X是标准化后的，j＝1,2,...,l为模型的回归系数，l为模型回归系数个数，则Lasso估计定义为：

$\begin{array}{ccc}\widehat{\mathrm{\β}}=\arg \min \left\{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(y_i-\mathrm{\α}-\underset{j=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\β}}_j{x}_{ij}\right)}^2\right\},& S.T.& \underset{j=1}{\overset{l}{\mathrm{\Σ}}}\left|{\mathrm{\β}}_j\right|\≤t\end{array}$

式中t为调和参数，α为回归模型参数。对一切的t，有α的估计 为应变量的平均值，不失一般性，假定当t的取值比较小时，某些相关度较低的系数就被压缩为0，从而达到特征选择的目的。然而，参数t不易确定，本文采用Lars方法解决Lasso思想的计算问题。

Lars方法的算法原理为：

Step1：假设估计的模型则当前的残差为y；

Step2：采用逐步向前阶梯法找出和应变量y相关性最强的自变量x₁，作为选择出的第一个特征量；

Step3：寻找出自变量x₂，使得当前残差R与x₁和x₂之间的相关系数相等，作为选择出的第二个特征量；

Step4：选择新的寻找路径，使得选入模型的变量和当前残差R的系数均相等，即自变量x₁和x₂的角平分线；

Step5：寻找符合条件的相关性最强的自变量x₃作为第三个特征量，同时寻找新的前进路径，依次类推，直至所有的变量都选择出来。

其中，所述步骤105中分阶段建立偏最小二乘(PLS)定量模型，具体如下：

根据上述步骤104中获得分界点浓度值，将训练样本集划分为高浓度与低浓度两个样本集，分别采用相应训练样本集的平均光谱数据建立高浓度与低浓度PLS定量模型。

其中，所述步骤106中判别未知样本的近红外成像光谱数据中属于三聚氰胺的像元点，具体如下：

Step1：用近红外显微成像光谱仪采集未知样本的近红外成像光谱数据；

Step2：对于未知样本的近红外成像光谱数据，采用与训练集样本相同的预处理方法进行光谱处理，去除噪声、基线的干扰；

Step3：对于经预处理后的未知样本的光谱数据进行PLS_DA判别分析，得到未知样本成像光谱中属于三聚氰胺的像元点，获得所述样本的像元点个数比值。

其中，所述步骤107中检测未知样本的三聚氰胺浓度，具体如下：

Step1：根据上述步骤106中得到未知样本的像元点个数比值与步骤104中训练样本集中像元点个数比值与三聚氰胺浓度的线性关系曲线，确定定量模型阶段；

Step2：对于未知样本的光谱，采用与训练集样本相同的变量选择方法进行波长变量选择；

Step3：通过未知样本的平均光谱数据代入已建好的相应阶段的定量模型获得未知样本的三聚氰胺浓度的预测值。

本发明与现有技术相比的优点在于：本发明所提供的饲料原料中三聚氰胺掺假定量检测方法以纯饲料原料近红外光谱为背景光谱，在突出三聚氰胺光谱信息特征，减弱饲料原料纯物质光谱信息影响的同时，也保证了掺假样本的成像光谱数据均为纯像元，提高了定性鉴别三聚氰胺的可靠性和准确性；成像光谱技术可以实现微米级显微尺度上的光谱成像，对微小目标物的分析极为有利，实现了痕量检测；结合分类方法与成像光谱数据中微米级像元点的光谱信息，获得掺假饲料原料样本建立定量模型的分布规律，建立不同浓度训练样本的像元点个数比值与三聚氰胺浓度的线性关系，得到分阶段定量模型的分界点浓度值，采用各个样本的平均光谱数据进一步分阶段建立PLS定量模型，提高预测精度，降低检测限；采用波长变量选择方法提高定量分析精度，解决过拟合和高相关性问题，提高检测速度。

附图说明

图1为本发明所述的饲料原料中三聚氰胺掺假定量检测方法的流程图；

图2为本发明实施例的纯豆粕与三聚氰胺的原始近红外光谱曲线的示意图；

图3为本发明实施例的经预处理后的不同浓度掺假样本的近红外光谱曲线的示意图。其中：(a)为以白板近红外光谱为背景光谱；(b)为以纯豆粕近红外光谱为背景光谱；

图4为本发明实施例的PLS_DA判别模型的分类阈值的示意图；

图5为本发明实施例的不同浓度掺假样本经PLS_DA模型判别后的像元点个数比值与三聚氰胺浓度的线性关系的示意图；

图6为本发明实施例的以三聚氰胺浓度0.1％为分界点浓度值，分阶段的像元点个数比值与三聚氰胺浓度的线性关系示意图。其中：(a)为高浓度阶段的线性关系图；(b)为低浓度阶段的线性关系图；

图7为本发明实施例的以三聚氰胺浓度0.075％为分界点浓度值，分阶段的像元点个数比值与三聚氰胺浓度的线性关系示意图。其中：(a)为高浓度阶段的线性关系图；(b)为低浓度阶段的线性关系图；

图8为本发明实施例的高浓度样本的PLS定量模型的预测结果示意图。其中：(a)为基于近红外全波段光谱的预测结果；(b)为基于特征波长变量的预测结果；

图9为本发明实施例的低浓度样本的PLS定量模型的预测结果示意图。其中：(a)为基于近红外全波段光谱的预测结果；(b)为基于特征波长变量的预测结果。

具体实施方式

下面结合附图和具体实施例对本发明的技术方案进一步详细阐述。

本发明所提供的饲料原料中三聚氰胺掺假定量检测方法的流程图，如图1所示，主要包括以下步骤：

步骤101，饲料原料采用豆粕，以纯豆粕的近红外光谱为背景，采用近红外显微成像光谱仪采集不同浓度掺假样本的近红外成像光谱数据，构成训练样本集。

这里，本发明采用英国PerkinElmer公司的FITRMicroscopeSpotlight400显微近红外成像光谱仪，选取市场上饲用的饲料原料豆粕与纯度为99％的三聚氰胺，分别配制纯豆粕样本、纯三聚氰胺样本及混合比例为2％-0.001％的掺假样本，对样本进行压片处理。采集被测样本的近红外成像光谱数据，设置实验参数为：成像模式下采集反射光谱、空间分辨率为25×25μm、光谱分辨率为8cm^-1、扫描速度为1.0cm/s、扫描波段为4000-7800cm^-1，波长变量个数为951个，光谱矩阵中每一行各元素代表一个样本在各个波长下的吸光度值。本实施例中训练样本集为20个，其原始光谱的光谱矩阵X的图形如图2所示，其中分别为纯豆粕与三聚氰胺的原始近红外光谱曲线，可看出波数为6816、5006、4627和4491cm^-1处的三聚氰胺光谱曲线吸光度较高，呈尖峰，在6713和6559cm^-1波数处的光谱吸光度较低，呈明显的小峰，而纯豆粕的近红外光谱图较宽，无尖峰，其中6850-6500cm^-1之间的吸收峰为三聚氰胺分子中的-NH2键的一倍频引起的，5150-4100cm^-1之间的多个吸收峰是由分子中含氮杂环和-NH2键的合频引起的，为准确鉴别出豆粕中的三聚氰胺提供了依据。

步骤102，对所述近红外成像光谱数据进行预处理，得到预处理后的训练样本集。

对训练样本集光谱数据进行预处理的目的是主要是用来去除上述样本光谱测量数据中由于仪器、环境等造成的噪声及基线漂移，从而提高光谱质量。预处理的方法为SG平滑方法与自适应迭代权重-偏最小二乘方法。

其中，所用的自适应迭代权重-偏最小二乘(air-PLS)方法的原理及步骤如下：

该方法是通过迭代改变拟合基线与原始信号之间的总体方差权重实现的，而该权重由当前拟合基线和原始信号之间的差异得到，因此它是个不断拟合和逼近的过程。air-PLS方法主要为采用惩罚最小二乘算法将信号平滑的原理应用到基线拟合，然后自适应迭代将惩罚过程转化为一个基线逼近的过程。

惩罚最小二乘(penalizedleastsquares)算法原理为：

假设x是要分析的光谱信号，z是平滑后信号，它们的长度均为n，i＝1,2,…,n，z对于x的保真度用二者间的误差平方和表示为：

$F=\underset{i=1}{\overset{n}{\Σ}}{(x_i-z_i)}^2---\left(1\right)$

z的粗糙度用其差分平方和表示为：

$R=\underset{i=2}{\overset{n}{\Σ}}{(z_i-z_{i-1})}^2=\underset{i=1}{\overset{n-1}{\Σ}}{\left({\mathrm{\Δz}}_i\right)}^2---\left(2\right)$

保真度和粗糙度之间的平衡可以用准确度加上惩罚的粗糙度，λ为惩罚系数，式子表示为：

Q＝F+λR＝||x-z||²+λ||Dz||²(3)

对其求偏导数且令其为零，可得一个易解的线性方程为：

(I+λD′D)z＝x(4)

其中I为单位矩阵；D为差分矩阵；D′为D的转置。

要实现最小二乘算法进行基线校正，引入保真度的权重向量w，并将其在有峰的位置设置为零，其中m为x与z的长度，i＝1,2,…,m，则z对x的保真度变为：

$F=\underset{i=1}{\overset{m}{\Σ}}{w}_i{(x_i-z_i)}^2=(X-Z{)}^{\′}W(X-Z)---\left(5\right)$

自适应迭代重加权原理：

自适应迭代重加权方法与加权最小二乘以及迭代惩罚最小二乘相似，不同的是采用了不同的方法计算权重，并增加了一个惩罚项来控制拟合基线的平滑度，优化的目标函数为：

$Q^t=\underset{i=1}{\overset{m}{\Σ}}{w_i}^t|x_i-{z_i}^t|+\mathrm{\λ}\underset{j=2}{\overset{m}{\Σ}}|z_j^t-z_{j-1}^t{|}^2---\left(6\right)$

其中i＝1,2,…,m；j＝1,2,…,m；t为迭代次数；λ为惩罚系数。

迭代的过程如下：设置初始权重采用下式进行迭代：

$w_i^t=\left\{\begin{array}{cc}0& x_i\&GreaterEqual;z_i^{t-1}\\ e^{\frac{t(x_i-z_i^{t-1})}{\left|d^t\right|}}& x_i<z_i^{t-1}\end{array}\right.---\left(7\right)$

其中向量d^t包含了原始光谱(x)与迭代过程中的拟合矢量差值为负的所有元素。在前t-1次的迭代中，拟合矢量是背景估计的一个候选值。如果当前计算大于此背景估计的一个候选值，则被视为处于峰的位置之中，此时权重会重置为零，使其不影响下一次迭代。在air-PLS算法过程中，迭代和重加权不断的自动执行，就可自动地、逐渐地消除处于峰位置之中的数据点，将背景点在权重矢量w中保留下来。迭代终止的条件可以是设定的迭代次数或者是一个收敛条件。

本次实施例中，SG平滑处理过程采用的平滑窗口点数为7，其预处理后不同浓度掺假样本分别以白板近红外光谱和纯豆粕近红外光谱为背景光谱的光谱曲线图如图3所示，其中图(a)为基于白板光谱背景的不同三聚氰胺浓度的掺假豆粕样本经过预处理后的近红外光谱曲线，图(b)为基于纯豆粕光谱背景的不同三聚氰胺浓度的掺假豆粕样本经过预处理后的近红外光谱曲线，可看出以纯豆粕的近红外光谱作为背景光谱采集的掺假豆粕样本的光谱曲线随着三聚氰胺浓度的变化呈现规律性的变化(如(b)图中波数5006cm^-1与6816cm^-1处)，并且不同光谱背景下同浓度掺假样本的光谱曲线特征波长处的强度明显增强，而以白板的近红外光谱为背景采集的掺假样本的光谱曲线相似，特征峰强度的变化不明显。

步骤103，采用偏最小二乘判别分析(PLS_DA)方法得到所述训练样本的近红外成像光谱数据中属于三聚氰胺的像元点，得到三聚氰胺像元点个数与总像元点个数的比值。

这里，采用纯豆粕与三聚氰胺的平均光谱数据建立PLS_DA模型的训练集X，分别对应类别值1与0，构成类别矩阵Y，得到关系式为：

Y＝BX+E(8)

其中B为回归系数；E为噪声。首先预测纯豆粕样本的成像数据数据中各像元点的类别值为1，保证纯豆粕样本的判别准确性，从而确定分类阈值，如图4所示，其中采用PLS_DA判别模型鉴别纯豆粕与三聚氰胺样本的成像光谱数据中的各个像元点信息，得到各个像元点的类别值，可看出选取的分类阈值0.3正确地判别了纯豆粕与三聚氰胺样本的类别，保证了PLS_DA模型判别的正确性。

其次，得到分类阈值后，采用PLS_DA模型判别分析训练样本集中不同掺假样本的成像光谱中属于三聚氰胺的像元点，分类结果如表1所示；同时根据下述公式得到三聚氰胺像元点个数与总像元点个数的比值：

$p=\frac{{\stackrel{\&OverBar;}{x}}_{\begin{array}{cc}mela\min e& pixels\end{array}}}{{\stackrel{\&OverBar;}{X}}_{\begin{array}{cc}all& pixels\end{array}}}---\left(9\right)$

其中p为像元点个数比值；为训练样本的成像光谱数据中属于三聚氰胺的像元点个数；为所述样本的总像元点个数。

表1PLS-DA模型的分类结果

其中n为同浓度的掺假样本个数；C为掺假样本中三聚氰胺的浓度值；为n个掺假样本的成像光谱数据中属于三聚氰胺像元点的平均个数；SD为三聚氰胺像元点个数的标准偏差；为n个掺假样本的成像光谱数据中总像元点的平均个数；p为与的比值。

步骤104，建立所述训练样本集的像元点个数比值与三聚氰胺浓度的线性关系，获得分阶段定量模型的分界点浓度值。

将不同浓度训练样本的像元点个数比值与相应的样本中三聚氰胺浓度进行线性拟合，同时得到线性关系式与相关系数值，公式如下：

p＝aC+b(10)

$R^2=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}(p_i-\stackrel{\&OverBar;}{p})(C_i-\stackrel{\&OverBar;}{C})}{\sqrt{\underset{i=1}{\overset{n}{\&Sigma;}}{(p_i-\stackrel{\&OverBar;}{p})}^2\underset{i=1}{\overset{n}{\&Sigma;}}{(C_i-\stackrel{\&OverBar;}{C})}^2}}---\left(11\right)$

其中p为像元点个数比值；C为样本中三聚氰胺的浓度；a和b为拟合系数；R为相关系数；n为样本个数；i＝1,2,…,n；p_i为第i个样本的像元点个数比值；为像元点个数比值平均值；C_i为第i个样本的三聚氰胺浓度；为样本浓度平均值。

本次实施例中，采用所述步骤103中PLS_DA模型得到的像元点个数比值与三聚氰胺浓度建立线性拟合曲线，如图5所示，得到二者的线性函数关系，相关系数为0.64。观察所述拟合曲线，局部选取类似分界点浓度值0.1％与0.075％，分别进行分阶段线性拟合，得到以分界点浓度值0.1％建立的分阶段拟合曲线如图6所示，图中(a)为三聚氰胺浓度值大于0.1％的训练集样本建立的线性拟合曲线，(b)为三聚氰胺浓度小于0.1％的训练集样本建立的线性拟合曲线，二者的相关系数均达到0.901；而以分界点浓度值0.075％建立的分阶段拟合曲线如图7所示，图中(a)为三聚氰胺浓度值大于0.075％的训练集样本建立的线性拟合曲线，(b)为三聚氰胺浓度小于0.075％的训练集样本建立的线性拟合曲线，二者的相关系数分别为0.69与0.974。比较图6与图7中的相关系数值，选取分界点浓度值为0.1％。

步骤105，对训练样本的近红外成像光谱数据进行变量选择。

对训练样本集光谱数据进行变量选择的目的是主要是用来降低光谱数据的维数和消除冗余信息，从而提高定量分析速度与精度。波长变量选择的方法为最小角回归-最小绝对收缩选择(Lars_Lasso)方法。

Lasso方法的算法原理为：

设有数据(Xⁱ,y_i),i＝1,2,...,n，n为样本个数，其中的X和y分别为第i个观测值的自变量和应变量，假设各观测值是彼此独立且X是标准化后的，j＝1,2,...,l为模型的回归系数，l为模型回归系数个数，则Lasso估计定义为：

$\begin{array}{ccc}\widehat{\mathrm{\&beta;}}=\arg \min \left\{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(y_i-\mathrm{\&alpha;}-\underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_j{x}_{ij}\right)}^2\right\},& S.T.& \underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}\left|{\mathrm{\&beta;}}_j\right|\&le;t\end{array}---\left(12\right)$

式中t为调和参数，α为回归模型参数。对一切的t，有α的估计 为应变量的平均值，不失一般性，假定当t的取值比较小时，某些相关度较低的系数就被压缩为0，从而达到特征选择的目的。然而，参数t不易确定，本发明采用Lars方法解决Lasso思想的计算问题。

Lars方法的算法原理为：

Step1：假设估计的模型则当前的残差为y；

Step2：采用逐步向前阶梯法找出和应变量y相关性最强的自变量x₁，作为选择出的第一个特征量；

Step3：寻找出自变量x₂，使得当前残差R与x₁和x₂之间的相关系数相等，作为选择出的第二个特征量；

Step4：选择新的寻找路径，使得选入模型的变量和当前残差R的系数均相等，即自变量x₁和x₂的角平分线；

Step5：寻找符合条件的相关性最强的自变量x₃作为第三个特征量，同时寻找新的前进路径，依次类推，直至所有的变量都选择出来。

本次实施例中，将训练样本集的近红外光谱的波长变量由951个减少为15个，实现了显著性变量的选择和对应参数的估计，模型更易于理解，较好地解决了最小二乘建模中过度拟合和高相关性的问题。

步骤105，分阶段建立偏最小二乘(PLS)定量模型，根据上述步骤104中获得分界点浓度值0.1％，将训练样本集划分为高浓度与低浓度两个样本集，分别采用相应训练样本集的平均光谱数据建立高浓度与低浓度PLS定量模型。

步骤106，判别未知样本的近红外成像光谱数据中属于三聚氰胺的像元点。

进一步地，所述判别未知样本的近红外成像光谱数据中属于三聚氰胺的像元点，包括如下步骤：

Step1：用近红外显微成像光谱仪采集未知样本的近红外成像光谱数据；

Step2：对于未知样本的近红外成像光谱数据，采用与训练集样本相同的预处理方法进行光谱处理，去除噪声、基线的干扰；

Step3：对于经预处理后的未知样本的光谱进行PLS_DA判别分析，得到未知样本成像光谱中属于三聚氰胺的像元点，获得所述样本的像元点个数比值。

本实施例中，采集34个未知的样本光谱，采用与训练集样本相同的SG平滑与air_PLS方法进行预处理，平滑窗口点数为7，代入PLS_DA模型判别分析未知样本的像元点个数比值。

表2未知样本的定量模型阶段判定结果

其中C为掺假样本中三聚氰胺的浓度；p为掺假样本的像元点个数比值。

步骤107，检测未知样本的三聚氰胺浓度。

这里，所述的检测未知样本的三聚氰胺浓度，是指将未知样本通过已建立的定量模型，得到该未知样本的三聚氰胺浓度的预测值，包括如下步骤：

Step1：根据上述步骤106中得到未知样本的像元点个数比值与步骤104中训练样本集中像元点个数比值与三聚氰胺浓度的线性关系曲线，确定定量模型阶段；

Step2：对于未知样本的光谱，采用与训练集样本相同的变量选择方法进行波长变量选择；

Step3：通过未知样本的平均光谱数据代入已建好的相应阶段的定量模型获得未知样本的三聚氰胺浓度的预测值。

本实施例中，对于经PLS_DA模型判别分析后的未知样本的像元点个数比值选择定量模型，采用与训练模型相同的Lars_Lasso方法进行波长变量选择，并且获得未知样本的平均光谱数据，通过相应定量模型预测得到该36个未知样本的三聚氰胺浓度，得到采用PLS方法建立的高浓度与低浓度样本定量模型的预测值与性能参数表。

为与未进行波长变量选择的PLS方法进行比较，采用PLS和PLS_Lars_Lasso方法建立的高浓度样本定量模型的实验结果如表3所示，预测结果如图8所示，图中(a)为基于全波段光谱的PLS模型的预测结果，(b)为基于特征波长变量的PLS模型的预测结果。

表3高浓度样本的PLS定量模型的试验结果

采用PLS和PLS_Lars_Lasso方法建立的低浓度样本定量模型的实验结果如表4所示，预测结果如图9所示，图中(a)为基于全波段光谱的PLS模型的预测结果，(b)为基于特征波长变量的PLS模型的预测结果。

表4低浓度样本的PLS定量模型的试验结果

高浓度与低浓度样本的PLS定量模型中，训练样本集的近红外光谱全谱的波长变量为951个，经过Lars-Lasso方法变量选择后，用于PLS回归建模的波长变量为15个，定量预测结果中PLS-Lars-Lasso定量模型的均方根误差位于0.009-0.105％之间，决定系数位于0.96-0.99之间，残余预测误差为1.77和2.99，均方根误差较小，且决定系数高，实现了显著性变量的选择和对应参数的估计，优于未进行波长变量选择的PLS定量模型。同时，通过PLS-Lars-Lasso分阶段定量模型的预测结果图8(b)和图9(b)中可看出，分阶段定量模型的预测值稳定、可靠。

采用本发明中的方法，对掺假豆粕样本进行了精确定量，预测结果的相关性达到0.96以上，预测误差低于0.105，且定量检测限达到了0.0025％。在已发表的他人文献中，采用像元点个数比值与样本浓度的拟合函数进行浓度预测，仅实现了半定量检测，定量检测限为0.05％，本发明的技术效果明显优于现有方法。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (8)

1.一种饲料原料中三聚氰胺掺假定量检测方法，其特征在于包括以下步骤：

步骤101，以纯饲料原料的近红外光谱为背景，采用近红外显微成像光谱仪采集不同浓度掺假样本的近红外成像光谱数据，构成训练样本集；

步骤102，对所述近红外成像光谱数据进行预处理，以去除噪声、基线的干扰；

步骤103，采用偏最小二乘判别分析(PLS_DA)方法得到所述训练样本的近红外成像光谱数据中属于三聚氰胺的像元点，得到三聚氰胺像元点个数与总像元点个数的比值；

步骤104，建立不同浓度训练样本的像元点个数比值与三聚氰胺浓度之间的线性关系，得到分阶段定量模型的分界点浓度值；

步骤105，对训练样本的近红外成像光谱数据进行变量选择，根据所述的分界点浓度值，采用训练样本的平均光谱数据分阶段建立高浓度和低浓度的偏最小二乘(PLS)定量模型；

步骤106，用建立的PLS_DA模型判别未知样本的成像光谱中属于三聚氰胺的像元点，得到像元点个数比值；

步骤107，根据所述的像元点个数比值确定定量模型阶段，对未知样本的近红外成像光谱数据进行变量选择后，取各个未知样本的平均光谱数据，用所构造的偏最小二乘定量模型，检测未知样本的三聚氰胺浓度。

2.根据权利要求1所述饲料原料中三聚氰胺掺假定量检测方法，其特征在于：所述步骤102中对所述近红外成像光谱数据进行预处理采用SG平滑过程与自适应迭代权重-偏最小二乘方法进行基线校正。

3.根据权利要求1所述饲料原料中三聚氰胺掺假定量检测方法，其特征在于：所述步骤103是采用PLS_DA方法判别训练样本的近红外成像光谱数据中属于三聚氰胺的像元点，获得三聚氰胺像元点个数与总像元点个数的比值，具体如下：

设已经存在训练样本集，包含近红外成像光谱数据X(n×l)和类别矩阵Y(n×m)，其中n为样本，l为波长变量，m为类别，存在关系式记为：

Y＝BX+E

其中B为回归系数；E为噪声，训练样本集中，Y矩阵的类别值分别用“1”和“0”表示，其中1表示该像元为纯饲料原料像元点，0表示该像元为三聚氰胺像元点；选取纯饲料原料与三聚氰胺的平均光谱数据作为训练集，建立PLS-DA模型训练已知的纯饲料原料样本光谱数据，设定分类阈值使得纯饲料原料样本近红外成像光谱数据中各像元点的类别均为1，将该阈值作为最终掺假鉴别中PLS_DA模型的分类阈值，进而逐一判别掺假样本的近红外成像光谱数据中属于三聚氰胺的像元点，得到训练样本集的像元点个数比值，公式如下：

其中p为像元点个数比值；为训练样本的近红外成像光谱数据中属于三聚氰胺的像元点个数；为所述样本的总像元点个数。

4.根据权利要求1所述饲料原料中三聚氰胺掺假定量检测方法，其特征在于：所述步骤104中建立所述训练样本集的像元点个数比值与三聚氰胺浓度的线性关系，获得分阶段定量模型的分界点浓度值，具体如下：

将不同浓度训练样本的像元点个数比值与相应的样本中三聚氰胺浓度进行线性拟合，同时得到线性关系式与相关系数值，公式如下：

p＝aC+b

其中p为像元点个数比值；C为样本中三聚氰胺的浓度；a和b为拟合系数；R为相关系数；n为样本个数；i∈(1,2,…,n)；p_i为第i个样本的像元点个数比值；为像元点个数比值平均值；C_i为第i个样本的三聚氰胺浓度；为样本浓度平均值；

观察所述拟合曲线，局部选取类似分界点浓度值，分别进行分阶段线性拟合，比较不同分界点浓度值构成的分阶段拟合曲线的相关系数值，选择相关系数最好的曲线相对应的分界点浓度值作为分阶段定量模型的分界点。

5.根据权利要求1所述饲料原料中三聚氰胺掺假定量检测方法，其特征在于：所述步骤105中对训练样本的成像光谱数据进行变量选择采用最小角回归-最小绝对收缩选择(Lars_Lasso)方法。

6.根据权利1要求所述饲料原料中三聚氰胺掺假定量检测方法，其特征在于：所述步骤105中分阶段建立偏最小二乘(PLS)定量模型，具体如下：

根据上述步骤104中获得分界点浓度值，将训练样本集划分为高浓度与低浓度两个样本集，分别采用相应训练样本集的平均光谱数据建立高浓度与低浓度PLS定量模型。

7.根据权利要求1所述饲料原料中三聚氰胺掺假定量检测方法，其特征在于：所述步骤106中判别未知样本的近红外成像光谱数据中属于三聚氰胺的像元点，具体如下：

Step1：用近红外显微成像光谱仪采集未知样本的近红外成像光谱数据；

Step2：对于未知样本的近红外成像光谱数据，采用与训练集样本相同的预处理方法进行光谱处理，去除噪声、基线的干扰；

Step3：对于经预处理后的未知样本的近红外成像光谱数据进行PLS_DA判别分析，得到未知样本成像光谱数据中属于三聚氰胺的像元点，获得所述样本的像元点个数比值。

8.根据权利要求1所述饲料原料中三聚氰胺掺假定量检测方法，其特征在于：所述步骤107中检测未知样本的三聚氰胺浓度，具体如下：

Step1：根据上述步骤106中得到未知样本的像元点个数比值与步骤104中训练样本集中像元点个数比值与三聚氰胺浓度的线性关系曲线，确定定量模型阶段；

Step2：对于未知样本的近红外成像光谱数据，采用与训练集样本相同的变量选择方法进行波长变量选择；

Step3：通过未知样本的平均光谱数据代入已建好的相应阶段的定量模型获得未知样本的三聚氰胺浓度的预测值。