CN104237158B - 一种具有普适性的近红外光谱定性分析方法 - Google Patents

Abstract

一种具有普适性的近红外光谱定性分析方法，获取训练集和预测集样本近红外光谱，构造余弦字典，求解训练集和预测集样本在余弦字典下的表示系数，然后结合表示系数的稀疏度完成光谱重构，再根据重构后的训练集近红外光谱及训练集样本类别编号寻找定性分类的有效投影方向，并完成重构后训练集和预测集近红外光谱的投影变换，接下来计算变换后的训练集和预测集近红外光谱的图拉普拉斯矩阵，并将分类问题转化成优化问题，最后完成分类函数的求解及执行预测集近红外光谱的定性分析。本发明提供了一种精度较高、准确性良好、有效适用于小样本复杂体系的具有普适性的近红外光谱定性分析方法。

Description

一种具有普适性的近红外光谱定性分析方法

技术领域

本发明涉及化工、食品、农业、环境、生物等行业的近红外仪器分析和检测领域，特别是一种近红外光谱定性分析方法。

背景技术

近红外光谱分析技术是利用物质中敏感基团对不同波长近红外光的吸收强度来进行定性分析的方法，近年来在农业、食品、药物、医学、石油化工等领域有着广泛的应用。各个领域应用近红外分析技术后，优势明显，主要表现如下几个方面：分析测试速度有较大幅度地提高，测试时间通常在几分钟内；测试效率有较大提升，通过一次光谱测量，可对样本多个性质同时进行测量；测试成本进一步降低，近红外技术是无损检测，除仪器供电外无其他试剂、预处理损耗。近红外分析技术被越来越多的分析工作者认可和使用，近红外光谱技术作为在线快速分析技术在生产过程和工业控制中发挥了的巨大作用，应用潜力巨大。将近红外光谱与待检测参数建立定性分析模型是近红外法检测的关键技术。

在定性模型建立方面，由于光谱中不同波长变量间存在复共线性关系，即某个波长点的信息可以用其他波长点线性表示，这些波长点的存在对现有的定性分析方法存在干扰作用，降低建模精度。此外，由于光谱波长点过多，而不同波长点之间的信息差别较小，在计算过程中易引起散度、方差矩阵奇异。

为解决上述问题及完成近红外光谱的定性分析，光谱分析工作者主要使用两类方法，一类是无监督的方法，以根据距离测度分类的聚类分析为典型代表，之后出现了K均值聚类、模糊K均值聚类等动态聚类算法。采用了迭代分级聚类策略思想的系统聚类分析，神经网络中的自组织神经网络也逐步在近红外光谱聚类分析中应用。另一类是有监督的方法，根据已知样本距离判断的最小距离判别法、K最邻近法，根据样本之间的相关性判断的线性判别分析法，在主成分分析基础上进行拟合判别的SIMCA方法等模式识别方法在近红外光谱分析文献中较为常见。同时对定量分析方法进行适当改造，将定量值定义为特定的类别编号，就可将定量分析方法用于定性分析，比如判别偏最小二乘法(PLS-DA)，支持向量机，BP神经网络等。

现有的神经网络、支持向量机等定性分析方法难以拟合样本光谱和类别特征存在的非线性关系，容易造成预测失败。光谱各吸收峰之间存在交叠现象，加上在光谱采集过程中由于光谱仪器本身或者外界因素的干扰，不可避免地存在噪声，这仍旧给光谱的定性析造成很大困难。通常近红外法应用的检测对象，如农业、食品的产品，由于种类的多样性和不确定性，在定性模型建立过程中，在特定波长下，朗伯比尔定律约束具有多样性、交叠性，光谱建模就从单一体系问题演变成为复杂体系问题，加大分类预测的难度。

样本数量对于复杂体系下近红外分类模型精度有较大的影响，常规近红外定性分析方法需通过收集大量样品建立分类模型，使模型尽可能地涵盖样品信息，从而具有较好的稳定性和适用性。在工程实际中，每次取得的样本不一定处于同一体系，则大多数近红外光谱定性分析成为一小样本建模问题，其预测精度和稳定性尚待提升。同时获得样本物理化学类别特性信息的过程将会消耗大量的时间、人力和财力，在极端情况下，部分样本光谱的物理化学类别信息可能无法获得。而获得未知物理化学类别信息的样本光谱较为容易，仅需一次或多次的光谱采集。

如上所述，现用于提高近红外光谱定性分析结果精度的方法和途径存在各种局限性和不足，制约了近红外光谱分析技术在复杂体系中的应用范围。目前大多应用实践尚不能满足近红外光谱定性分析中需大量已知物理化学类别特性信息样本的要求，且针对复杂体系的应用对象，无法提供标准物理化学类别特性信息的样本。那么为了提高近红外光谱定性分析的精度和稳定性，必须发展新的近红外光谱定性分析方法，克服噪声、吸收峰混叠、复共线性的干扰，从未知分类信息的样本光谱提取更多信息，实现复杂混合体系近红外光谱的准确定性分析，对于当前近红外光谱定性分析的研究和应用具有十分重要的意义。

发明内容

为了克服已有近红外光谱定性分析方法的精度较低、准确性较差、不能适用于小样本复杂体系的不足,本发明提供了一种精度较高、准确性良好、有效适用于小样本复杂体系的具有普适性的近红外光谱定性分析方法。

本发明解决其技术问题所采用的技术方案是：

一种具有普适性的近红外光谱定性分析方法，所述方法包括以下步骤：

(1)获取训练集和预测集样本近红外光谱，将训练集近红外光谱表示为X(n×p)，其中n为训练集样本数量，p为光谱数据维度；将预测集近红外光谱表示为X'(m×p)，其中m为预测集样本数量；X(i,j)表示训练集近红外光谱矩阵中坐标位置为(i,j)的吸光度值，X'(i,j)表示预测集近红外光谱矩阵中坐标位置为(i,j)的吸光度值；构造训练集近红外光谱的类属性矩阵Y(n×K)，Y(i,j)表示类属性矩阵中坐标位置为(i,j)的类别属性，其中K为校正集样本中类别数；Y矩阵的每一行对应一个训练集样本，若第t号训练集样本的类别编号为z，则Y(t,j＝z)＝1，

(2)用离散余弦函数g_a[b]＝cos(bπa/(10*n*m))逼近近红外光谱中余弦信号特征的吸收峰并构造字典中原子g_a，每个原子g_a为一p维列向量，其中b＝0,1,...,p-1，下标a＝0,1,...,10*n*m-1代表频率参量；在10*n*m个原子形成后，使用原子g_a组成具有过完备特性的字典矩阵G，字典矩阵G的计算式为：

G(p×(10*n*m))＝[g₀,g₁,g₂...g₁₀*_n*_m-1]；

(3)逐一提取训练集近红外光谱矩阵X(n×p)的每个样本光谱，即光谱矩阵X(n×p)的p维行向量，向量记为x_c，c＝1,2...n；逐一提取预测集近红外光谱矩阵X'(m×p)的每个样本光谱，即光谱矩阵X'(m×p)的p维行向量，向量记为x'_d，d＝1,2...m；计算每个近红外光谱x_c和x'_d在字典矩阵G下的稀疏表示系数α_c和α_d，c＝1,2...n，d＝1,2...m，其中α_c和α_d均为10*n*m维行向量；

(4)使用稀疏表示系数向量α_c重构训练集近红外光谱矩阵中的行向量x'_c＝(G*α_c)^T，c＝1,2...n，其中，上标T为向量或矩阵转置运算；使用稀疏表示系数向量α_d重构预测集近红外光谱矩阵中的行向量x″_d＝(G*α_d)^T，d＝1,2...m，其中上标T为向量或矩阵转置运算；使用重构后的x'_c，c＝1,2...n，逐行替换训练集近红外光谱矩阵X(n×p)中的行向量，形成新的训练集近红外光谱矩阵X_new(n×p)；使用重构后的x″_d，d＝1,2...m逐行替换预测集近红外光谱矩阵X'(m×p)中的行向量，形成新的预测集近红外光谱矩阵X'_new(m×p)；

(5)利用训练集近红外光谱数据，寻找实施定性判别的q个投影方向w_k，w_k为p维向量，k＝1,2...q该方向能对光谱类别信息数据进行有效分类；

(6)将新的训练集近红外光谱矩阵X_new(n×p)和新的预测集近红外光谱矩阵X'_new(m×p)向向量w_k，k＝1,2,3...q，方向进行投影变换X_cal＝(X_neww₁,X_neww₂...X_neww_q)、X_tst＝(X'_neww₁,X'_neww₂...X'_neww_q)，得到变换后的训练集近红外光谱矩阵和预测集近红外光谱矩阵，分别记作X_cal(n×q)和X_tst(m×q)；

(7)将变换后的训练集近红外光谱矩阵X_cal(n×q)和预测集近红外光谱矩阵X_tst(m×q)按行合并为一全体光谱矩阵X_all((n+m)×q)，逐一提取全体光谱矩阵X_all((n+m)×q)的每个样本光谱，即光谱矩阵X_all((n+m)×q)的q维行向量，向量记为x_ii或x_jj，ii,jj＝1,2...n+m；计算关联矩阵W，W_ii,jj表示关联矩阵中坐标位置为(ii,jj)的值，ii,jj＝1,2...n+m且W_ii,jj＝exp(||x_ii-x_jj||²/0.48)，其中||||为向量2范数；计算对角矩阵D，D_ii,jj表示对角矩阵中坐标位置为(ii,jj)的值，且D_ii,jj＝0(ii≠jj)ii,jj＝1,2...n+m；计算n+m阶图拉普拉斯方阵L，L＝D-W；

(8)将预测集样本的定性分析转换成 $\underset{f\∈H_k}{\min }\frac{1}{n}\underset{\mathrm{ii}=1}{\overset{n}{\mathrm{\Σ}}}\max (\mathrm{0,1}-\mathrm{yf}\left(x_{\mathrm{ii}}\right))+{\mathrm{\γ}}_K{\left|\right|f\left|\right|}_K^2+{\mathrm{\γ}}_I{f}^T\mathrm{Lf}$ 优化问题，其中；γ_K为控制空间复杂度参数，γ_I为控制结构复杂度参数，H_k为希尔伯特空间，x_ii为光谱矩阵X_all((n+m)×q)的q维行向量，ii＝1,2...n+m，x_ii的前n个样本由训练集近红外光谱组成，y为前n个x_ii光谱向量对应样本的类别编号向量；f(x_ii)为x_ii的分类函数，该函数在输入后m个预测集光谱向量x_ii时，预测x_ii光谱对应样本的类别编号；

(9)使用拉格朗日乘子法，计算 $\underset{f\∈H_k}{\min }\frac{1}{n}\underset{\mathrm{ii}=1}{\overset{n}{\mathrm{\Σ}}}\max (\mathrm{0,1}-\mathrm{yf}\left(x_{\mathrm{ii}}\right))+{\mathrm{\γ}}_K{\left|\right|f\left|\right|}_K^2+{\mathrm{\γ}}_I{f}^T\mathrm{Lf}$ 优化问题的解，得到ii＝1,2...n+m，其中x_ii为光谱矩阵X_all((n+m)×q)的q维行向量，α_ii为核表示系数，e为截距，均由拉格朗日乘子法求得，K'(,)为核函数；

(10)逐一提取预测集近红外光谱矩阵X_tst(m×q)的每个样本光谱，即光谱矩阵X_tst(m×q)的q维行向量，向量记为x_u，u＝1,2...m；将m个预测集光谱x_u代入分类函数中，即可获得预测集光谱的类别编号，完成定性分析。

进一步，所述的步骤(3)中，x_c和x'_d在字典下G的稀疏表示系数α_c和α_d计算具体过程为：

(3-1)x_c在字典矩阵G下的稀疏表示系数α_c求解过程等价于优化问题，其中||α_c||₀表示α_c的l₀范数，即α_c中非零元素的个数，l₀范数上限T为稀疏度；初始化残差r＝x_c、清空索引集和重组字典

(3-2)计算稀疏表示系数α_c；

(3-2-1)：找出残差r与字典中内积最大原子g_a所对应的标号其中(,)为向量内积运算，并更新索引集Λ＝Λ∪{λ}和重组字典G'＝[G',g_λ]；

(3-2-2)：由最小二乘法计算当前x_c向量的稀疏表示系数α_c＝argmin₂||x_c-G'α_c||₂并更新残差r＝x_c-G'α_c；

(3-2-3)：重复(3-2-1)和(3-2-2)T次后停止迭代；

(3-3)x'_d在字典矩阵G下的稀疏表示系数α_d求解过程等价于优化问题，其中||α_d||₀表示α_d的l₀范数，即α_d中非零元素的个数；初始化残差r＝x'_d、清空索引集和重组字典

(3-4)计算稀疏表示系数α_j；

(3-4-1)：找出残差r与字典中内积最大原子g_a所对应的标号其中(,)为向量内积运算，并更新索引集Λ＝Λ∪{λ}和重组字典G'＝[G',g_λ]；

(3-4-2)：由最小二乘法计算当前x'_d向量的稀疏表示系数α_d＝argmin₂||x'_d-G'α_d||₂并更新残差r＝x'_d-G'α_d；

(3-4-3)：重复(3-4-1)和(3-4-2)T次后停止迭代。

再进一步，所述的步骤(5)的具体过程为：

(5-1)计算K方阵D'＝1/nY^TY；定义K行k列矩阵Q_k(K×k)，k＝1,2,3...q，Q₁作为Q_k的初始矩阵为全1矩阵；q为投影方向个数，同时也代表Q_k矩阵更新过程的次数；随机生成一个K阶过渡方阵θ_*，计算K阶方阵

(5-2)按照k＝1,2,3...q的取值顺序，执行q次两步循环法求解p维投影列向量w_k；

(5-2-1)利用弹性网眼回归(ElasticNet)优化求解方法求解优化问题

$\underset{w_k}{\arg \min }\{{\left|\right|Y{\mathrm{\θ}}_k-X_{\mathrm{new}}{w}_k\left|\right|}^2+0.1*w_k^{T}I{w}_k+0.2*{\left|\right|w_k\left|\right|}_1\},$ 其中I为p阶单位矩阵，符号“||||₁”表示求向量1范数运算；优化求解结果为此次循环的w_k：

(5-2-2)计算矩阵 ${\mathrm{\θ}}_k^{\′}=(I-Q_k{Q}_k^T{D}^{\′})D^{\′-1}{Y}^T{X}_{\mathrm{new}}{w}_k;$ 更新矩阵 ${\mathrm{\θ}}_k={\mathrm{\θ}}_k^{\′}/\sqrt{{\mathrm{\θ}}_k^{\′T}{D}^{\′}{\mathrm{\θ}}_k^{\′}}$ 和矩阵Q_k+1＝(Q_k:θ_k)，其中符号“：”表示向量或矩阵的联结运算。

所述的步骤(8)的控制空间复杂度参数γ_k＝0.9，控制结构复杂度参数γ_I＝1；

所述的步骤(9-10)中，核函数为K'(yy,zz)＝exp(-||yy-zz||²/2σ²)，yy、zz表示核函数的输入参数，σ为平滑度，σ＝0.35。

与现有技术相比，本发明首先构造了余弦冗余字典，并求取了近红外光谱数据在该字典下的稀疏表示系数，通过控制稀疏度，完成光谱信号的重构，滤除了近红外光谱仪器本身或者外界因素的噪声，再通过构造一弹性网眼优化问题，确定了利于定性分类的最优投影方向，随即完成近红外光谱数据投影变换，克服了复共线性、光谱混叠、光谱非线性的干扰，最后使用核函数表示理论，将压缩后的训练集和预测集样本共同形成定性分类优化问题并求解该问题完成定性分类，由于本发明提出的方法，让预测集样本在分类优化问题中参与了模型的辅助训练，提取了其内部结构和空间信息，丰富了定性分析过程中的已知信息量，因此本发明的方法能取得相对较优的分类准确率；此外，由于预测集样本的使用，本发明方法需要的训练集样本数量可缩减，适合在近红外光谱小样本复杂体系建模时使用。

附图说明

图1为本发明所述一种具有普适性的近红外光谱定性分析方法的流程图；

图2为所采集的199个油品样本近红外光谱图；

图3为参数γ_K和参数γ_I在区间搜索时定性分类平均准确率的等高线图。

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1～图3，一种具有普适性的近红外光谱定性分析方法，由字典构造、表示系数求解、光谱重构、投影方向搜索、投影变换、图拉普拉斯矩阵计算、分类优化问题构造及分类函数求解等部分组成，图1给出了本发明方法的总体实现框图，其具体包括以下步骤：

(1)获取训练集和预测集样本近红外光谱，将训练集近红外光谱表示为X(n×p)，其中n＝20为训练集样本数量，p＝2051为光谱数据维度；将预测集近红外光谱表示为X'(m×p)，其中m＝179为预测集样本数量，p＝2051为光谱数据维度；X(i,j)表示训练集近红外光谱矩阵中坐标位置为(i,j)的吸光度值，X'(i,j)表示预测集近红外光谱矩阵中坐标位置为(i,j)的吸光度值；构造训练集近红外光谱的类属性矩阵Y(n×K)，Y(i,j)表示类属性矩阵中坐标位置为(i,j)的类别属性，其中K＝2为校正集样本中类别数；类属性矩阵Y的每一行对应一个训练集样本，若第t号训练集样本的类别编号为z＝0,1，则Y(t,j＝z)＝1，其中c＝0代表食用油样本，c＝1代表混有地沟油的食用油样本；图2为所采集的199个油品样本近红外光谱；

(2)用离散余弦函数g_a[b]＝cos(bπj/(10*n*m))逼近近红外光谱中余弦信号特征的吸收峰并构造字典中原子g_a，每个原子g_a为一p维列向量，其中b＝0,1,...,2050，下标a＝0,1,...,10*n*m-1代表频率参量；在10*n*m个原子形成后，使用原子g_a组成具有过完备特性的字典矩阵G，字典矩阵G的计算式为：

G(p×(10*n*m))＝[g₀,g₁,g₂...g₁₀*_n*_m-1]；

(3)逐一提取训练集近红外光谱矩阵X(n×p)的每个样本光谱，即光谱矩阵X(n×p)的p＝2051维行向量，向量记为x_c，c＝1,2...20；逐一提取预测集近红外光谱矩阵X'(m×p)的每个样本光谱，即光谱矩阵X'(m×p)的p＝2051维行向量，向量记为x'_d，d＝1,2...179；计算每个近红外光谱x_c和x'_d在字典下G的稀疏表示系数α_c和α_d，c＝1,2...20，d＝1,2...179，其中α_c和α_d均为10*20*179维行向量；在此具体实施例中，x_c和x'_d在字典矩阵G下的稀疏表示系数α_c和α_d计算具体过程为：

(3-1)x_c在字典矩阵G下的稀疏表示系数α_c求解过程等价于优化问题，其中||α_c||₀表示α_c的l₀范数，即α_c中非零元素的个数，l₀范数上限T为稀疏度；初始化残差r＝x_c、清空索引集和重组字典

(3-2)计算稀疏表示系数α_c；

(3-2-1)：找出残差r与字典中内积最大原子g_a所对应的标号其中(,)为向量内积运算，并更新索引集Λ＝Λ∪{λ}和重组字典G'＝[G',g_λ]G'＝[G',g_λ]，λ表示残差r与字典中内积最大原子g_a所对应的标号；

(3-2-2)：由最小二乘法计算当前x_c向量的稀疏表示系数α_c＝argmin₂||x_c-G'α_c||₂并更新残差r＝x_c-G'α_c；

(3-2-3)：重复(3-2-1)和(3-2-2)T次后停止迭代；

(3-3)x'_d在字典矩阵G下的稀疏表示系数α_d求解过程等价于优化问题，其中||α_d||₀表示α_d的l₀范数，即α_d中非零元素的个数；初始化残差r＝x'_d、清空索引集和重组字典

(3-4)计算稀疏表示系数α_j；

(3-4-1)：找出残差r与字典中内积最大原子g_a所对应的标号其中(,)为向量内积运算，并更新索引集Λ＝Λ∪{λ}和重组字典G'＝[G',g_λ]；

(3-4-2)：由最小二乘法计算当前x'_d向量的稀疏表示系数α_d＝argmin₂||x'_d-G'α_d||₂并更新残差r＝x'_d-G'α_d；

(3-4-3)：重复(3-4-1)和(3-4-2)T次后停止迭代。

(3-5)：所述的稀疏度T＝6；

(4)使用稀疏表示系数向量α_c重构训练集近红外光谱矩阵中的行向量x'_c＝(G*a_c)^T，c＝1,2...n，其中上标T为向量或矩阵转置运算；使用稀疏表示系数向量α_d重构预测集近红外光谱矩阵中的行向量x″_d＝(G*α_d)^T，d＝1,2...179，其中上标T为向量或矩阵转置运算；使用重构后的x'_c，c＝1,2...20逐行替换训练集近红外光谱矩阵X(n×p)中的行向量，形成新的训练集近红外光谱矩阵X_new(n×p)；使用重构后的x″_d，d＝1,2...179逐行替换预测集近红外光谱矩阵X'(m×p)中的行向量，形成新的预测集近红外光谱矩阵X'_new(m×p)；

(5)利用训练集近红外光谱数据，寻实施找定性判别的q个投影方向w_{k＝1,2,3...q},w_k为p维向量,该方向能对光谱类别信息数据进行有效分类；

(5-1)计算K方阵D'＝1/nY^TY；定义K行k列矩阵Q_k(K×k)，k＝1,2,3...q，Q₁作为Q_k的初始矩阵为全1矩阵；q为投影方向个数，同时也代表Q_k矩阵更新过程的次数；随机生成一个K阶过渡方阵θ_*，计算K阶方阵

(5-2)按照k＝1,2,3...q的取值顺序，执行q次两步循环法求解p维投影列向量w_k；

(5-2-1)利用弹性网眼回归(ElasticNet)优化求解方法求解优化问题

$\underset{w_k}{\arg \min }\{{\left|\right|Y{\mathrm{\θ}}_k-X_{\mathrm{new}}{w}_k\left|\right|}^2+0.1*w_k^{T}I{w}_k+0.2*{\left|\right|w_k\left|\right|}_1\},$ 其中I为p阶单位矩阵，符号“||||₁”表示求向量1范数运算；优化求解结果为此次循环的w_k：

(5-2-2)计算矩阵 ${\mathrm{\θ}}_k^{\′}=(I-Q_k{Q}_k^T{D}^{\′})D^{\′-1}{Y}^T{X}_{\mathrm{new}}{w}_k;$ 更新矩阵 ${\mathrm{\θ}}_k={\mathrm{\θ}}_k^{\′}/\sqrt{{\mathrm{\θ}}_k^{\′T}{D}^{\′}{\mathrm{\θ}}_k^{\′}}$ 和矩阵Q_k+1＝(Q_k:θ_k)，其中符号“：”表示向量或矩阵的联结运算。

(5-3)所述的投影方向个数q＝8；

(6)将新的训练集近红外光谱矩阵X_new(n×p)和新的预测集近红外光谱矩阵X'_new(m×p)向向量w_k，k＝1,2,3...8，方向进行投影变换X_cal＝(X_neww₁,X_neww₂...X_neww₈)、X_tst＝(X'_neww₁,X'_neww₂...X'_neww₈)，得到变换后的训练集近红外光谱矩阵和预测集近红外光谱矩阵，分别记作X_cal(20×8)和X_tst(179×8)；

(7)将变换后的训练集近红外光谱矩阵X_cal(20×8)和预测集近红外光谱矩阵X_tst(179×8)按行合并为一全体光谱矩阵X_all(199×8)，逐一提取全体光谱矩阵X_all(199×8)的每个样本光谱，即光谱矩阵X_all(199×8)的q＝8维行向量，向量记为x_ii或x_jj，ii,jj＝1,2...199；计算关联矩阵W，W_ii,jj表示关联矩阵中坐标位置为(ii,jj)的值，ii,jj＝1,2...199且W_ii,jj＝exp(||x_ii-x_jj||²/0.48)，其中||||为向量2范数；计算对角矩阵D，D_ii,jj表示对角矩阵中坐标位置为(ii,jj)的值，且D_ii,jj＝0(ii≠jj)ii,jj＝1,2...n+m；计算n+m阶图拉普拉斯方阵L，L＝D-W；

(8)将预测集样本的定性分析转换成 $\underset{f\∈H_k}{\min }\frac{1}{n}\underset{\mathrm{ii}=1}{\overset{n}{\mathrm{\Σ}}}\max (\mathrm{0,1}-\mathrm{yf}\left(x_{\mathrm{ii}}\right))+{\mathrm{\γ}}_K{\left|\right|f\left|\right|}_K^2+{\mathrm{\γ}}_I{f}^T\mathrm{Lf}$ 优化问题，其中；γ_K为控制空间复杂度参数，γ_I为控制结构复杂度参数，H_k为希尔伯特空间，x_ii为光谱矩阵X_all(199×8)的q＝8维行向量，ii＝1,2...199，x_ii的前n＝20个样本由训练集近红外光谱组成，y为前n＝20个x_ii光谱向量对应样本的类别编号向量；f(x_ii)为x_ii的分类函数，该函数在输入后m＝179个预测集光谱向量x_ii时，预测x_ii光谱对应样本的类别编号；所述的步骤(8)的控制空间复杂度参数γ_k＝0.9，控制结构复杂度参数γ_I＝1；

(9)使用拉格朗日乘子法，计算 $\underset{f\∈H_k}{\min }\frac{1}{n}\underset{\mathrm{ii}=1}{\overset{n=20}{\mathrm{\Σ}}}\max (\mathrm{0,1}-\mathrm{yf}\left(x_{\mathrm{ii}}\right))+{\mathrm{\γ}}_K{\left|\right|f\left|\right|}_K^2+{\mathrm{\γ}}_I{f}^T\mathrm{Lf}$ 优化问题的解，得到ii＝1,2...199，其中x_ii为光谱矩阵X_all((n+m)×q)的q＝8维行向量，α_ii为核表示系数，e为截距，均由拉格朗日乘子法求得，K'(,)为核函数；

(10)逐一提取预测集近红外光谱矩阵X_tst(179×8)的每个样本光谱，即光谱矩阵X_tst(179×8)的q＝8维行向量，向量记为x_u，u＝1,2...179；将179个预测集光谱x_u代入分类函数中，即可获得预测集光谱的类别编号，完成定性分析。

进一步，所述的步骤(9-10)的核函数为K'(yy,zz)＝exp(-||yy-zz||²/2σ²)，yy、zz表示核函数的输入参数，σ为平滑度，σ＝0.35。

本实施实例收集了花生油、大豆油、菜籽油、葵花仁油、芝麻油五个品种的食用油样本各20个和由绍兴市质量技术监督局查获含地沟油的食用油样本99个，累计199个样本分别装入100ml烧杯中，使用MPA傅立叶变换型近红外光谱仪(德国Bruker)及自带光程2mm的光纤附件伸入烧杯中测量每个样本的近红外光谱。光谱测量范围为2500cm^-1到4000cm^-1，分辨率设定为4cm^-1,每个光谱扫描16次并取平均值。该光谱数据被划分为训练集样本20个，验证集样本179个。训练集样本中，随机选取了食用油样本10个和地沟油样本10个，构造了一复杂体系的小样本问题。

比较本发明提出的近红外光谱定性分析方法与使用偏最小二乘判别分析(PLSDA)和支持向量机(SVM)传统方法在分类准确度的性能差异。偏最小二乘判别分析和支持向量机只使用了训练集20个样本建立了定性分析模型，而本发明提出的方法，在步骤(7)-(10)，让预测集样本参与了模型的辅助训练，提取了其内部结构和空间信息，丰富了定性分析过程中的已知信息量。表1给出了不同方法的预测结果，结果表明本发明的方法取得了相对最优分类准确率，说明了本发明方法的有效性和优越性。

图3显示了空间复杂度参数γ_K和结构复杂度参数γ_I在区间{0.5e-6,1e-6,0.5e-5,1e-5,0.5e-4,1e-4,0.5e-3,1e-3,0.5e-2,1e-2,0.5e-1,1e-1,0.5151020304050}搜索时，利用所述的一种具有普适性的近红外光谱定性分析方法对177个预测集食用油和地沟油样本执行100次定性分类平均准确率的等高线图。在控制空间复杂度参数γ_k＝0.9，控制结构复杂度参数γ_I＝1区间内，等高线图显示本发明方法最优的分类准确率为96.67％。

表1为本发明方法最终分类准确率与常规方法的比较结果

表1。

Claims (5)

1.一种具有普适性的近红外光谱定性分析方法，其特征在于：所述方法包括以下步骤：

(1)获取训练集和预测集样本近红外光谱，将训练集近红外光谱矩阵表示为X(n×p)，其中n为训练集样本数量，p为光谱数据维度；将预测集近红外光谱矩阵表示为X'(m×p)，其中m为预测集样本数量；X(i,j)表示训练集近红外光谱矩阵中坐标位置为(i,j)的吸光度值，X'(i,j)表示预测集近红外光谱矩阵中坐标位置为(i,j)的吸光度值；构造训练集近红外光谱的类属性矩阵Y(n×K)，Y(i,j)表示类属性矩阵中坐标位置为(i,j)的类别属性，其中K为校正集样本中类别数；Y矩阵的每一行对应一个训练集样本，若第t号训练集样本的类别编号为z，则Y(t,j＝z)＝1，

(2)用离散余弦函数g_a[b]＝cos(bπa/(10*n*m))逼近近红外光谱中余弦信号特征的吸收峰并构造字典中原子g_a，每个原子g_a为一p维列向量，其中b＝0,1,...,p-1，a＝0,1,...,10*n*m-1代表频率参量；在10*n*m个原子形成后，使用原子g_a组成具有过完备特性的字典矩阵G，字典矩阵G的计算式为：

G(p×(10*n*m))＝[g₀,g₁,g₂...g_10*n*m-1]；

(3)逐一提取训练集近红外光谱矩阵X(n×p)的每个样本光谱，即光谱矩阵X(n×p)的p维行向量，向量记为x_c，c＝1,2...n；逐一提取预测集近红外光谱矩阵X'(m×p)的每个样本光谱，即光谱矩阵X'(m×p)的p维行向量，向量记为x'_d，d＝1,2...m；计算每个近红外光谱x_c和x'_d在字典矩阵G下的稀疏表示系数α_c和α_d，c＝1,2...n，d＝1,2...m，其中α_c和α_d均为10*n*m维行向量；

(4)使用稀疏表示系数α_c重构训练集近红外光谱矩阵中的行向量x'_c＝(G*α_c)^T，c＝1,2...n，其中，上标T为向量或矩阵转置运算；使用稀疏表示系数α_d重构预测集近红外光谱矩阵中的行向量x”_d＝(G*α_d)^T，d＝1,2...m，其中上标T为向量或矩阵转置运算；使用重构后的x'_c，c＝1,2...n，逐行替换训练集近红外光谱矩阵X(n×p)中的行向量，形成新的训练集近红外光谱矩阵X_new(n×p)；使用重构后的x″_d，d＝1,2...m逐行替换预测集近红外光谱矩阵X'(m×p)中的行向量，形成新的预测集近红外光谱矩阵X'_new(m×p)；

(5)利用训练集近红外光谱数据，寻找实施定性判别的q个投影方向w_k，w_k为p维向量，k＝1,2...q，该方向能对光谱类别信息数据进行有效分类；

(6)将新的训练集近红外光谱矩阵X_new(n×p)和新的预测集近红外光谱矩阵X'_new(m×p)向向量w_k，k＝1,2,3...q，方向进行投影变换X_cal＝(X_neww₁,X_neww₂...X_neww_q)、X_tst＝(X'_neww₁,X'_neww₂...X'_neww_q)，得到变换后的训练集近红外光谱矩阵和预测集近红外光谱矩阵，分别记作X_cal(n×q)和X_tst(m×q)；

(7)将变换后的训练集近红外光谱矩阵X_cal(n×q)和预测集近红外光谱矩阵X_tst(m×q)按行合并为一全体光谱矩阵X_all((n+m)×q)，逐一提取全体光谱矩阵X_all((n+m)×q)的每个样本光谱，即光谱矩阵X_all((n+m)×q)的q维行向量，向量记为x_ii或x_jj，ii,jj＝1,2...n+m；计算关联矩阵W，W_ii,jj表示关联矩阵中坐标位置为(ii,jj)的值，ii,jj＝1,2...n+m且W_ii,jj＝exp(||x_ii-x_jj||²/0.48)，其中||||为向量2范数；计算对角矩阵D，D_ii,jj表示对角矩阵中坐标位置为(ii,jj)的值，且D_ii,jj＝0(ii≠jj)ii,jj＝1,2...n+m；计算n+m阶图拉普拉斯方阵L，L＝D-W；

(8)将预测集样本的定性分析转换成

优化问题，其中；γ_K为控制空间复杂度参数，γ_I为控制结构复杂度参数，H_k为希尔伯特空间，x_ii为光谱矩阵X_all((n+m)×q)的q维行向量，ii＝1,2...n+m，x_ii的前n个样本由训练集近红外光谱组成，y为前n个x_ii光谱向量对应样本的类别编号向量；f(x_ii)为x_ii的分类函数，该函数在输入后m个预测集光谱向量x_ii时，预测x_ii光谱对应样本的类别编号；

(9)使用拉格朗日乘子法，计算

优化问题的解，得到ii＝1,2...n+m，其中x_ii为光谱矩阵X_all((n+m)×q)的q维行向量，α_ii为核表示系数，e为截距，均由拉格朗日乘子法求得，K'(,)为核函数；

(10)逐一提取变换后的预测集近红外光谱矩阵X_tst(m×q)的每个样本光谱，即光谱矩阵X_tst(m×q)的q维行向量，向量记为x_u，u＝1,2...m；将m个预测集光谱x_u代入分类函数中，即可获得预测集光谱的类别编号，完成定性分析。

2.如权利要求1所述的一种具有普适性的近红外光谱定性分析方法，其特征在于：所述的步骤(3)中，x_c和x'_d在字典矩阵G下的稀疏表示系数α_c和α_d计算具体过程为：

(3-1)x_c在字典矩阵G下的稀疏表示系数α_c求解过程等价于优化问题，其中||α_c||₀表示α_c的l₀范数，即α_c中非零元素的个数，l₀范数上限T为稀疏度；初始化残差r＝x_c、清空索引集和重组字典

(3-2)计算稀疏表示系数α_c；

(3-2-1)：找出残差r与字典中内积最大原子g_a所对应的标号其中(,)为向量内积运算，并更新索引集Λ＝Λ∪{λ}和重组字典G'＝[G',g_λ]；

(3-2-2)：由最小二乘法计算当前x_c向量的稀疏表示系数α_c＝argmin₂||x_c-G'α_c||₂并更新残差r＝x_c-G'α_c；

(3-2-3)：重复(3-2-1)和(3-2-2)T次后停止迭代；

(3-3)x'_d在字典矩阵G下的稀疏表示系数α_d求解过程等价于优化问题，其中||α_d||₀表示α_d的l₀范数，即α_d中非零元素的个数；初始化残差r＝x'_d、清空索引集和重组字典

(3-4)计算稀疏表示系数α_d；

(3-4-1)：找出残差r与字典中内积最大原子g_a所对应的标号其中(,)为向量内积运算，并更新索引集Λ＝Λ∪{λ}和重组字典G'＝[G',g_λ]；

(3-4-2)：由最小二乘法计算当前x'_d向量的稀疏表示系数α_d＝argmin₂||x'_d-G'α_d||₂并更新残差r＝x'_d-G'α_d；

(3-4-3)：重复(3-4-1)和(3-4-2)T次后停止迭代。

3.如权利要求1或2所述的一种具有普适性的近红外光谱定性分析方法，其特征在于：所述的步骤(5)的具体过程为：

(5-1)计算K方阵D'＝1/nY^TY；定义K行k列矩阵Q_k(K×k)，k＝1,2,3...q，Q₁作为Q_k的初始矩阵为全1矩阵；q为投影方向个数，同时也代表Q_k矩阵更新过程的次数；随机生成一个K阶过渡方阵θ_*，计算K阶方阵

(5-2)按照k＝1,2,3...q的取值顺序，执行q次两步循环法求解p维投影列向量w_k；

(5-2-1)利用弹性网眼回归优化求解方法求解优化问题

其中I为p阶单位矩阵，符号“||||₁”表示求向量1范数运算；优化求解结果为此次循环的w_k：

(5-2-2)计算矩阵更新矩阵和矩阵Q_k+1＝(Q_k:θ_k)，其中符号“：”表示向量或矩阵的联结运算。

4.如权利要求1或2所述的一种具有普适性的近红外光谱定性分析方法，其特征在于：所述的步骤(8)的控制空间复杂度参数γ_k＝0.9，控制结构复杂度参数γ_I＝1。

5.如权利要求1或2所述的一种具有普适性的近红外光谱定性分析方法，其特征在于：所述的步骤(9)-(10)中，核函数为K'(yy,zz)＝exp(-||yy-zz||²/2σ²)，yy、zz表示核函数的输入参数，σ为平滑度，σ＝0.35。