CN108984930B - 一种用于近红外光谱分析的推进稀疏偏最小二乘方法 - Google Patents

Abstract

本发明提出一种用于近红外光谱分析的推进稀疏偏最小二乘方法，其包括使用蒙特卡洛交叉验证和10折交叉验证确定稀疏偏最小二乘因子数和稀疏参数；根据预测误差确定样品权重；使用稳健步骤减少异常样品影响；确定使用平方形式损失函数；对未知样品，将多个模型预测结果排序，根据模型权重，确定其最终结果。本发明可以用于近红外光谱分析，提高方法的预测精度和稳定性。

Description

一种用于近红外光谱分析的推进稀疏偏最小二乘方法

技术领域

本发明涉及分析化学领域的无损分析，特别涉及一种用于近红外光谱分析的推进稀疏偏最小二乘方法。

背景技术

近红外光谱是介于可见光和中红外光之间的电磁波谱，波数约为：10000～4000cm^-1。近红外光谱法是利用含有氢基团(X-H，X为：C，O，N，S等)化学键(X-H)伸缩振动倍频和合频，在近红外区的吸收光谱，通过选择适当的化学计量学多元校正方法，把校正样品的近红外吸收光谱与其成分浓度或性质数据进行关联，建立校正样品吸收光谱与其成分浓度或性质之间的关系-校正模型。在进行未知样品预测时，应用已建好的校正模型和未知样品的吸收光谱，就可定量预测其成分浓度或性质。另外，通过选择合适的化学计量学模式识别方法，也可分离提取样本的近红外吸收光谱特征信息，并建立相应的类模型。在进行未知样品的分类时，应用已建立的类模型和未知样品的吸收光谱，便可定性判别未知样品的归属。

建立一个稳健、可靠的校正模型在近红外光谱定量分析中是十分重要的，但是当建模样品较少，或者所研究问题很复杂时，单个模型经常得不到较好的分析结果。研究发现，多模型共识方法通过整合多个模型即子模型的结果可以克服单个模型存在的预测结果不准确、不可靠的缺点。同时，由于多模型共识方法在提高模型精度、降低过拟合方面具有优势，引起了广泛关注。因为多模型共识方法通过提高子模型的预测能力和提高子模型之间的差异来提高预测能力和稳定性，所以研究工作，如重采样集合(bagging)、推进(boosting)等，集中在如何选择建模样品构建模型。这两种方法都是通过随机取样技术从建模集中构建子模型，不同的是：推进方法通过整合多个预测结果、优化样品取样权重的方式降低模型的方差和偏差；而重采样集合方法对所有样品同等对待用于降低单个样品对模型的影响。

变量选择是近红外光谱分析中一个十分复杂的问题，对一组特定的数据，可以选择多种变量选择方法进行分析，如区间偏最小二乘方法(iPLS)可以用于选择一个或几个波段用于建模；无信息变量消除-偏最小二乘法(UVE-PLS)、随机检验-偏最小二乘法(RT-PLS)、基于信息向量的变量选择法(OPS)通过变量的重要性评价值进行变量选择；遗传算法-偏最小二乘(GA-PLS)、蚁群算法-偏最小二乘法(ACO-PLS)等通过优化方法进行变量选择；连续投影算法-多元线性回归(SPA-MLR)方法选择共线性最小的变量组进行变量选择；竞争性自适应权重取样(CARS)通过迭代竞争的方式逐步选择变量；弹性网-偏最小二乘法(EN-PLS)对信息相关变量进行分组的方法进行变量选择。虽然这些方法都已经成功应用于近红外光谱定量分析中，但是它们大多需要较长的计算时间，难以实现在线实时分析。

现行的很多建模方法大多只关注多模型方法带来的稳健性，如boostingPLS等，或者只关注变量选择方法带来的精简模型或者预测误差降低，如UVE、CARS、GA-PLS等，较少有结合二者来进行，尤其是在多模型过程中引入变量选择的方法很少。

发明内容

鉴于上述现行建模方法的不足，本发明的目的在于提供一种用于近红外光谱分析的推进稀疏偏最小二乘方法，通过推进(boosting)选样方式，可以提高模型的稳定性，从而减少模型过拟合；通过使用稀疏偏最小二乘方法，可以在建模过程中实现变量选择，精简模型，降低模型预测误差。

为达上述目的，本发明一种用于近红外光谱分析的推进稀疏偏最小二乘方法，其包括以下步骤：

步骤S1、用X₀表示光谱数据矩阵，y表示光谱对应的测量数据向量，设参与建模的样品总数为m，其对应的光谱变量数为p，初始化每个样品权重WS为1/m，并设定推进取样量、迭代次数T，稀疏偏最小二乘模型因子数n和稀疏参数γ，其中稀疏偏最小二乘模型因子数n设定与偏最小二乘方法一致，用蒙特卡罗交叉验证确定，稀疏参数通过10折交叉验证确定，样品权重是指被选择参与建模的样品的权重，推进取样量是建模样品数量占样品总数的比例；

步骤S2、将参与建模的样品总数m乘以推进取样量得选取样品数，按照轮盘赌方式选取与选取样品数一致的选取样品；

步骤S3、使用选取样品建立迭代次数为i时的稀疏偏最小二乘回归模型M_i，并计算各个光谱变量数p对应的回归系数，其具体包括：

S3-1：初始化回归系数b为p×1的0向量，令k＝1，X＝X₀；

S3-2：计算

其中T表示矩阵转置；

S3-3：

S3-4：更新变量集S＝{j,1≤j≤p,w_j≠0}U{j,1≤j≤p,b_j≠0}，其中w_j表示第j个变量对应的权重，b_j表示第j个变量对应的回归系数；

S3-5：以现有变量集S，使用偏最小二乘算法，建立回归模型，计算出回归系数b_s，X的光谱数据权重P_s；

S3-6：计算

S3-7：升级X，X_i＝X_s,j，i∈S，i与j要对应；

S3-8：升级b，b_i＝b_s,j，i∈S，i与j要对应；

S3-9：令k＝k+1，重复S3-3-S3-8，直至k>n时终止；

步骤S4、根据稀疏偏最小二乘回归模型M_i计算迭代次数为i时所有m个参与建模的样品的预测误差

其中，y_k,i为迭代次数为i时第k个样品的化学成分的测量值，

为迭代次数为i时第k个样品的采用稀疏偏最小二乘回归模型M_i得到的化学成分的预测值；

步骤S5、根据公式

计算稀疏偏最小二乘回归模型M_i的模型权重，其中

为模型损失函数，WS_k,i为迭代次数为i时第k个样品的样品权重；

步骤S6、根据公式

计算迭代次数为i+1时所有m个参与建模的样品的样品权重；

步骤S7、i＝i+1，重复S2-S6，直至完成T次迭代，得到T个稀疏偏最小二乘回归模型及对应的模型权重；

步骤S8、对一个未知样品，利用T个稀疏偏最小二乘回归模型预测未知样品的化学成分，得到T个样品的化学成分的预测值y_i，(i＝1,2,3,…,T)，并对y_i进行升序排列，满足下式的第r个预测值作为最终的预测结果：

即WM需要按照预测值y_i的升序排序，依次增加

直至累计和刚好大于

时对应的样品预测值，为该样品的化学成分的预测值。

上述方法中，设定推进取样量为50％，迭代次数T为100。

上述方法中，步骤S4中如果e_k,i≥median(e_i)+γ_r.MAD(e_i)，则令e_k,i＝0，其中MAD是测量值与中位值绝对偏差的中位值，γ_r是稳健参数，优选的取值4。

上述方法中，步骤S3的S3-5中使用SIMPLS偏最小二乘算法计算回归稀疏和光谱数据权重。

上述方法中，所述蒙特卡洛交叉验证，10折交叉验证之前对样品的光谱变量进行标准化，所述标准化就是将每个光谱变量的值减去所有光谱变量的平均值，然后再除以其标准偏差。

附图说明

图1显示为本发明的方法整体流程图；

图2显示为本发明的稀疏偏最小二乘方法流程图；

图3显示为本发明的实施例中交叉验证均方跟误差随稀疏参数变化图；

图4显示为本发明稀疏偏最小二乘回归系数变化图；

图5显示为本发明的增加稳健步骤对比图；

图6显示为本发明的未知样品预测示意图。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

实施例

本实施例是应用于近红外光谱定量分析，对小麦样本中的蛋白质成分含量值进行建模分析。具体的步骤如下：

数据由Williams提供，包含884个加拿大小麦样本的可见-近红外光谱，可从网址：http://www.idrc-chambersburg.org/shootout2008.html下载。样本的近红外光谱通过Foss Model 6500近红外光谱仪器测得，波长范围为400-2498nm，采样间隔2nm，包括1050个波长点。参考值在温尼伯华的谷物实验室测定。其中，第680和681号样本没有有效的蛋白质含量，为奇异样本。本研究考察蛋白质含量，并采用网站上对数据集划分的描述，775个样本用作建模集，107个样本用作预测集。建模集用于建立模型，预测集用于测试方法。

步骤S1、使用蒙特卡罗交叉验证，确定偏最小二乘因子数为12，即确定稀疏偏最小二乘模型因子数为12。图3为使用10折交叉验证确定稀疏参数图，可以看出稀疏参数取值0.95时，模型的交叉验证均方根误差(RMSECV)最小。所以稀疏参数γ定为0.95。设定迭代次数T为100，推进取样量为50％，共775个样品，每个样品权重WS初始化为1/775。

步骤S2、根据样品权重，使用轮盘赌方法，从建模集取出50％(推进取样量)样品，共387个。样品权重越大，越容易被选上。

步骤S3、使用选出的样品，建立稀疏偏最小二乘模型：

S3-1、用X₀表示光谱数据矩阵，包含1050个波长点数据(变量)，用p表示波长点数。

y表示光谱对应的测量数据向量。设定稀疏参数γ＝0.95，因子数12，初始化回归系数b

为1050×1的零向量，X＝X₀，k＝1；

S3-2、计算

其中T表示矩阵转置；

S3-3、根据Z计算w，

S3-4、更新选入变量集S＝{j,1≤j≤p,w_j≠0}U{j,1≤j≤p,b_j≠0}，其中w_j表示第j个变量

对应的权重，b_j表示第j个变量对应的回归系数；

S3-5、以现有变量集S，使用偏最小二乘算法，建立模型，计算出回归系数b_s，X权重

P_s；

S3-6、计算

S3-7、升级X，X_i＝X_s,j，i∈S，i与j要对应。

S3-8、升级b，b_i＝b_s,j，i∈S，i与j要对应。

S3-9、令k＝k+1，重复S3-3-S3-8，直至k＝12。此时的b即为计算的回归系数。

图4显示了因子数从1变化到12，对应的回归系数的变化。可以看出随着因子数的增大，选择的光谱变量逐渐增多。稀疏偏最小二乘方法在建模过程中实现了变量选择。

步骤S4、根据稀疏偏最小二乘回归模型M_i计算迭代次数为i时所有建模样品的预测误差

其中，y_k,i为迭代次数为i时第k个样品的蛋白质的测量值，

为迭代次数为i时第k个小麦样品的采用稀疏偏最小二乘回归模型M_i得到的化学成分的预测值。为了降低异常样品获取较大权重影响模型，进行了如下限制：

e_k,i＝0，如果e_k,i≥median(e_i)+γ_r.MAD(e_i)。

MAD是测量值与中位值绝对偏差的中位值。γ_r是稳健参数，这里取4。具体的，通过设置γ_r可以使预测误差特别大的样品的权重降低，有助于减少异常样品对模型的影响。本文后续也称该操作为稳健步骤。

图5为推进稀疏偏最小二乘不加稳健步骤和增加稳健步骤，经过连续100次迭代计算后，样品取样权重的箱须图。从图中看出，如果不加稳健步骤，在终止运行时部分样品权重会非常大，将近25％，如果该样品是异常样品将干扰模型，得到异常结果。增加稳健步骤后，样品的极端权重得到有效去除，所有样品权重都在5％以下。

步骤S5、根据公式

计算稀疏偏最小二乘回归模型M_i的模型权重，其中

WS_k,i为迭代次数为i时第k个样品的样品权重。其中，

为模型损失函数。

步骤S6、根据公式

计算迭代次数为i+1时所有建模样品的样品权重。

步骤S7、i＝i+1，重复S2-S6，直至完成T次迭代，得到T个稀疏偏最小二乘回归模型及对应的模型权重。

步骤S8、对一个未知样品，利用T个稀疏偏最小二乘回归模型预测未知样品的化学成分，得到T个样品的化学成分的预测值y_i，(i＝1,2,3,…,T)。对y_i进行升序排列，满足下式的第r个预测值作为最终的预测结果：

具体的，WM需要按照预测值y_i的升序排序，依次增加

直至累计和刚好大于

时对应的样品预测值，为该样品的化学成分的预测值。以预测集样品1为例，100个模型预测其光谱得到100个预测值，见图6中的(a)。将100个预测值按升序排列，见图6中的(b)，得到排序后序号；将

按照预测值升序序号排列，并计算

其中r为序号，r≤T，以

对序号r作图，见图6中的(c)。图6中的(c)标识的第一个

大于

位置(由竖线标出)，为最终预测值对应的序号，该值对应的图6中的(b)竖线和排序后预测值的交点，即为最终样品1的预测值。

下表1列出了推进稀疏偏最小二乘方法与偏最小二乘、稀疏偏最小二乘、推进偏最小二乘方法结果比较。从表1可以看出，四种方法预测效果排序为推进稀疏偏最小二乘>稀疏偏最小二乘>推进偏最小二乘>偏最小二乘。本发明提出的方法结合了推进算法和稀疏偏最小二乘算法的优点。

表1：推进稀疏偏最小二乘法和其它方法的比较

*^a偏最小二乘；

^b稀疏偏最小二乘；

^c推进偏最小二乘；

^d推进稀疏偏最小二乘；

^e相关系数；

^f预测均方根误差，其值越小，模型越好；

^g相对分析误差Residual Prediction Deviation，其值越大，模型越好；

^h相对标准偏差(100次随机计算统计结果)，其值越小，越好。

Claims (6)

1.一种用于近红外光谱分析的推进稀疏偏最小二乘方法，其特征在于包括以下步骤：

步骤S1、用X₀表示光谱数据矩阵，y表示光谱对应的测量数据向量，设参与建模的样品总数为m，其对应的光谱变量数为p，初始化每个样品权重WS为1/m，并设定推进取样量、迭代次数T，稀疏偏最小二乘模型因子数n和稀疏参数γ，其中稀疏偏最小二乘模型因子数n设定与偏最小二乘方法一致，用蒙特卡罗交叉验证确定，稀疏参数通过10折交叉验证确定，样品权重是指被选择参与建模的样品的权重，推进取样量是建模样品数量占样品总数的比例；

步骤S2、将参与建模的样品总数m乘以推进取样量得选取样品数，按照轮盘赌方式选取与选取样品数一致的选取样品；

步骤S3、使用选取样品建立迭代次数为i时的稀疏偏最小二乘回归模型M_i，并计算各个光谱变量数p对应的回归系数，其具体包括：

S3-1：初始化回归系数b为p×1的0向量，令k＝1，X＝X₀；

S3-2：计算

其中T表示矩阵转置；

S3-3：

S3-4：更新变量集S＝{j,1≤j≤p,w_j≠0}U{j,1≤j≤p,b_j≠0}，其中w_j表示第j个变量对应的权重，b_j表示第j个变量对应的回归系数；

S3-5：以现有变量集S，使用偏最小二乘算法，建立回归模型，计算出回归系数b_s，X的光谱数据权重P_s；

S3-6：计算X_s＝X_s(I-P_s(P_s ^TP_s)^-1P_s ^T)；

S3-7：升级X，X_i＝X_s,j，i∈S，i与j要对应；

S3-8：升级b，b_i＝b_s,j，i∈S，i与j要对应；

S3-9：令k＝k+1，重复S3-3-S3-8，直至k>n时执行S4步骤；

步骤S4、根据稀疏偏最小二乘回归模型M_i计算迭代次数为i时所有m个参与建模的样品的预测误差

其中，y_k,i为迭代次数为i时第k个样品的化学成分的测量值，

为迭代次数为i时第k个样品的采用稀疏偏最小二乘回归模型M_i得到的化学成分的预测值；

步骤S5、根据公式

计算稀疏偏最小二乘回归模型M_i的模型权重，其中

为模型损失函数，WS_k,i为迭代次数为i时第k个样品的样品权重；

步骤S6、根据公式

计算迭代次数为i+1时所有m个建模样品的样品权重；

步骤S7、i＝i+1，重复S2-S6，直至完成T次迭代，得到T个稀疏偏最小二乘回归模型及对应的模型权重；

步骤S8、对一个未知样品，利用T个稀疏偏最小二乘回归模型预测未知样品的化学成分，得到T个样品的化学成分的预测值y_i，(i＝1,2,3,…,T)，并对y_i进行升序排列，满足下式的第r个预测值作为最终的预测结果：

即WM需要按照预测值y_i的升序排序，依次增加

直至累计和刚好大于

时对应的样品预测值，为该样品的化学成分的预测值。

2.根据权利要求1所述的一种用于近红外光谱分析的推进稀疏偏最小二乘方法，其特征在于：设定推进取样量为50％，迭代次数T为100。

3.根据权利要求1所述的一种用于近红外光谱分析的推进稀疏偏最小二乘方法，其特征在于：步骤S4中如果e_k,i≥median(e_i)+γ_r.MAD(e_i)，则令e_k,i＝0，其中MAD是测量值与中位值绝对偏差的中位值，γ_r是稳健参数。

4.根据权利要求1所述的一种用于近红外光谱分析的推进稀疏偏最小二乘方法，其特征在于：所述稳健参数γ_r为4。

5.根据权利要求1所述的一种用于近红外光谱分析的推进稀疏偏最小二乘方法，其特征在于：步骤S3的S3-5中使用SIMPLS偏最小二乘算法计算回归稀疏和光谱数据权重。

6.根据权利要求1所述的一种用于近红外光谱分析的推进稀疏偏最小二乘方法，其特征在于：所述蒙特卡洛交叉验证，10折交叉验证之前对样品的光谱变量进行标准化，所述标准化就是将每个光谱变量的值减去所有光谱变量的平均值，然后再除以其标准偏差。

Images