CN105092519A - 基于增量偏最小二乘法的样品成份测定方法 - Google Patents

Abstract

本发明公开了一种基于增量偏最小二乘法的样品成份测定方法，包括以下步骤：S1，采集待测样品的近红外光谱数据；S2，通过增量偏最小二乘模型，获得该近红外光谱数据所对应的样品中各成分的含量。本发明通过利用增量偏最小二乘模型对待测样品的近红外光谱数据进行处理，从而即可获得该近红外光谱数据所对应的样品中各成分的含量，与采用传统的偏最小二乘模型进行数据处理相比，节约了时间和空间，所得回归系数基本相同，预测均方根误差却更小，可见，本发明中的增量偏最小二乘模型具有更高的预测精度和建模效率。此外，本发明采用梯度学习的方法寻找优化的回归系数，从而可以更快速的更新模型，提高模型适应新数据的能力。

Description

基于增量偏最小二乘法的样品成份测定方法

技术领域

本发明涉及一种成份测定方法，尤其是一种基于增量偏最小二乘法的样品成份测定方法。

背景技术

在食品、药品及石油化工产品的生产中，按照美国食品药品协会提出的过程分析技术(ProcessAnalysisTechnology，PAT)相关标准，需要对生产过程中的中间产品进行分析和检测，以明确其各物质含量、中间生成物状态及其变化规律，从而满足产品质量设计要求，生产出可靠的最终产品。近红外光谱检测技术以其快速、无损、低成本的特性成为过程分析技术中的通用方法，得到业界的广泛应用。其原理是：当近红外光照射(穿过或反射)待测样品时，样品分子中的基团吸收红外光产生振动，使偶极矩发生变化，在不同的波长处，吸光率不同，从而获得红外吸收光谱。对于溶液，按照朗伯比尔定理，吸光率和溶液的浓度呈线性关系。

在过程分析中，随着生产的进行，样本是逐个获得的。设采集的第i个样本点＜x_(i),y_(i)＞，x_(i)∈R^1×p是观测值，为红外光谱数据，y_(i)∈R^1×q是因变量，为x_(i)对应的物质浓度，前n个样本点构成校正集＜X⁽ⁿ⁾,Y⁽ⁿ⁾＞，自变量X⁽ⁿ⁾∈R^n×p为包含n个校正样品的光谱数据，Y⁽ⁿ⁾∈R^n×q代表物质浓度，如果用＜X,Y＞代表训练样本全集，随着＜X,Y＞样本数量不断增加，X与Y回归模型预测精度应该逐渐提高。

目前，过程分析技术主要采用化学计量学中的已有方法建立过程分析模型。但是，在PAT分析中，参与训练的样本具是逐步获得的，具有增量特性，与普通的化学计量学建模方法有所不同。其样本空间随着不同批次样本点的采集会发生改变，当新样本到达后，需要用这个新样本更新模型。因此，尽管偏最小二乘法(PartialLeastSquares，PLS)在化学计量学中得到了广泛应用，但是由于PLS采用批量学习模式，调整PLS参数时，需要抛弃已有模型，重新训练全部数据，并采用交叉验证等方法选择新参数，建立新模型，因而存在训练时间长、空间消耗大、建模效率低等问题；另外，应用PLS模型对样品中的成分含量进行预测时，预测精度有待进一步提高。

发明内容

本发明的目的在于，提供一种基于增量偏最小二乘法(IncrementalPartialLeastSquares，IPLS)的样品成份测定方法，它可以有效解决现有的PLS模型用于实际生产产品在线检测中面临的实际问题，尤其是采用PLS模型进行产品的增量数据检测时，其需要抛弃已有模型，重新训练全部数据，建立新模型，从而导致训练时间长、空间消耗大、建模效率低的问题，以及应用PLS模型对样品中的成分含量进行预测时，预测精度有待进一步提高的问题。

为解决上述技术问题，本发明采用如下的技术方案：一种基于增量偏最小二乘法的样品成份测定方法，包括以下步骤：

S1，采集待测样品的近红外光谱数据；

S2，通过增量偏最小二乘模型，获得该近红外光谱数据所对应的样品中各成分的含量。优选的，步骤S2中所述的增量偏最小二乘模型通过以下方法建立：

a.收集n个样本的近红外光谱数据及其对应的样品中各成分的含量数据作为初始校正集；

b.将初始校正集中的数据进行中心化处理，得初始中心化样本数据；

c.对所述的初始中心化样本数据进行偏最小二乘(PLS)回归，得初始回归系数及偏最小二乘初始回归模型；

d.当第i个样本点加入该校正集时，对其进行中心化处理，并利用偏最小二乘初始回归模型获得相应的预测值；其中，i＝n+1,n+2，……；

e.计算所述的预测值与真实值之间的误差；若该误差小于等于阀值δ，则令i＝i+1，转至d；否则对回归系数进行更新后，令i＝i+1，转至d，直至得到最终的回归系数B_i，即得增量偏最小二乘模型。

本发明的步骤e中，通过以下方式对回归系数进行更新，得最终的回归系数B_i：

其中，i＝n+1,n+2，……，B_i为第i个样本点＜x_(i),y_i＞加入校正集时得到的更新的回归系数，B_i-1为对所述的初始校正集中的n个训练样本数据进行偏最小二乘法回归时获得的初始回归系数，η为正常数，其大小由沿着梯度反方向搜索时的步长决定；为第i个样本点的中心化数据，为第i个样本点的预测值，为的转置。

本发明采用梯度学习的方法寻找优化的回归系数，一方面采取了增量学习的方法，另一方面通过上述的回归系数优化公式可以更快速的更新模型，提高模型适应新数据的能力。

上述方法的步骤e中，所述的增量偏最小二乘模型为：

y＝XB_i+E

其中，所述的y为待测样品中各成分的含量，X为待测样品的近红外光谱数据，B_i为增量偏最小二乘模型最终的回归系数，E为残差。

前述的基于增量偏最小二乘法的样品成份测定方法中，步骤e中所述的阀值δ采用K-折交叉验证的方式确定。

优选的，本发明具体通过以下方法确定阀值δ：

(1)将初始校正集中的n个样本数据平均分为k份，以第j份作为增量训练集，剩余的k-1份作为初始训练集进行交叉验证，其中，1≤j≤k(k≥4)；

(2)对所述的初始训练集中的数据进行中心化处理，得初始中心化样本数据；

(3)对初始中心化样本数据进行PLS回归，得初始回归系数及PLS初始回归模型；

(4)当第j份作为增量训练集中的数据加入该初始训练集时，对增量训练集中的数据进行中心化处理，并利用PLS初始回归模型获得相应的预测值；

(5)设δ＝2^-i，1≤i≤n_top，其中n_top∈N；

(6)计算所述的预测值与真实值之间的误差；若该误差小于等于δ，则转至(7)；否则对所述的回归系数进行更新，并将该更新后的回归系数应用到作为增量训练集的第j份数据上，获得相应的预测值及交叉验证均方差MSE(k,i)；将所述的交叉验证均方差MSE(k,i)放入均方差矩阵中；

(7)令i＝i+1，若i＜n_top转至(5)；否则，则令j＝j+1，转至(1)；若j＞k，得k行n_top列的均方差矩阵；

(8)对所述的k行n_top列的均方差矩阵的每一列求均值，获得均值矩阵；

(9)查找均值矩阵中的均方差最小值，该均方差最小值在矩阵中所对应的列号为i_best，相应的即为最优阀值δ。

采用本发明的上述方法确定的最优阀值δ，可以使得增量偏最小二乘模型的预测精度最高。

本发明中所述的待测样品为谷物、土壤、草或复方中药药物。

优选的，所述的待测样品为含有芍药苷的复方中药药物。

与现有技术相比，本发明具有以下优点：

1、通过利用增量偏最小二乘模型对待测样品的近红外光谱数据进行处理，从而即可获得该近红外光谱数据所对应的样品中各成分的含量，与采用传统的偏最小二乘模型进行数据处理相比，节约了时间和空间，所得回归系数基本相同，预测均方根误差却更小，可见，本发明中的增量偏最小二乘模型具有更高的预测精度和建模效率；

2、本发明采用梯度学习的方法寻找优化的回归系数，一方面采取了增量学习的方法，另一方面通过上述的回归系数优化公式可以更快速的更新模型，提高模型适应新数据的能力；

3、本发明中，引入了模型更新的阀值，对样本起到筛选作用，降低部分高密度样本反复取样对模型的影响，同时可以有效增加稀疏样本到模型，降低模型复杂度的同时，提高了模型的精度。本发明中通过采用K-折交叉验证的方式确定最佳阀值，可以使得增量偏最小二乘模型的预测精度更高，预测误差更小，同时节约了时间和空间，可以更好的进行增量光谱数据的定量分析；

4、本发明中，随着建模数据的增多，会比PLS节约更多的建模时间，同时，基于阀值的样本筛选机制，使得建立的IPLS模型变得更准确，预测精度更高；

5、本发明中的基于增量偏最小二乘模型的样品中成分含量的测定方法，对复方中药药物中芍药苷含量的测定非常有效，相对于PLS的预测值的改善程度达9.18％，说明了IPLS的提出，在复方中药药物数据集表现出了极强的预测能力；

6、本发明中的基于增量偏最小二乘模型的样品中成分含量的测定方法，对草中碳、氮、硫元素含量的测定非常有效，相对于PLS的预测值的改善程度达0.57％，说明了IPLS的提出，在草数据集表现出了极强的预测能力；

7、本发明中的基于增量偏最小二乘模型的样品中成分含量的测定方法，对土壤中有机质含量的测定非常有效，相对于PLS的预测值的改善程度达2.78％，说明了IPLS的提出，在土壤数据集表现出了极强的预测能力；

8、本发明中的基于增量偏最小二乘模型的样品中成分含量的测定方法，对谷物中水分、油脂、蛋白质、淀粉含量的测定非常有效，其中，测定谷物中的水分含量时，相对于PLS的预测值精度的改善程度为17.38％；测定谷物中的油脂含量时，相对于PLS的预测值的精度改善程度为9.29％；测定谷物中的蛋白质含量时，相对于PLS的预测值的精度改善程度为6.10％；测定谷物中的淀粉含量时，相对于PLS的预测值的精度改善程度为6.42％，说明了IPLS的提出，在谷物数据集表现出了极强的预测能力。

附图说明

图1是实验例1中前四个批次联合建模回归系数；

图2是实验例1中前二个批次阀值与交叉验证预测均方根误差的关系示意图；

图3是实验例1中前六个批次123条近红外光谱；

图4是实验例1中前六个批次联合建模所得回归系数；

图5是实验例1中前四个批次阀值与交叉验证预测均方根误差的关系示意图；

图6是实验例1中前四个批次联合建模对五、六批次的预测结果；

图7是实验例1中前六个批次联合建模对第七个批次的预测结果；

图8是实验例2中谷物数据的近红外光谱示意图；

图9是实验例2中土壤数据的近红外光谱示意图；

图10是实验例2中草数据的近红外光谱示意图；

图11是实验例2中corn(moisture)的回归系数对比示意图；

图12是实验例2中corn(oil)的回归系数对比示意图；

图13是实验例2中corn(protein)的回归系数对比示意图；

图14是实验例2中corn(starch)的回归系数对比图；

图15是实验例2中grass的回归系数对比示意图；

图16是corn(moisture)阀值与交叉验证预测均方根误差的的关系示意图；

图17是corn(oil)阀值与交叉验证预测均方根误差的关系示意图；

图18是corn(protein)阀值与交叉验证预测均方根误差的关系示意图；

图19是corn(starch)阀值与交叉验证预测均方根误差的关系示意图；

图20是grass阀值与交叉验证预测均方根误差的关系示意图；

图21是soil阀值与交叉验证预测均方根误差的关系示意图。

图22是实验例2中corn(moisture)预测值与真实值的比较结果示意图；

图23是实验例2中corn(oil)预测值与真实值的比较结果示意图；

图24是实验例2中corn(protein)预测值与真实值的比较结果示意图；

图25是实验例2中corn(starch)预测值与真实值的比较结果示意图；

图26是实验例2中grass预测值与真实值的比较结果示意图；

图27是实验例2中soil预测值与真实值的比较结果示意图；

图28为增量偏最小二乘模型的建模方法流程图；

图29为确定最佳阀值的方法流程图。

下面结合附图和具体实施方式对本发明作进一步的说明。

具体实施方式

本发明的实施例：一种基于增量偏最小二乘法的样品成份测定方法，如图28所示，包括以下步骤：

S1，采集待测样品的近红外光谱数据；

S2，通过增量偏最小二乘模型，获得该近红外光谱数据所对应的样品中各成分的含量；其中，所述的增量偏最小二乘模型通过以下方法建立：

a.收集n个样本的近红外光谱数据及其对应的样品中各成分的含量数据作为初始校正集；

b.将初始校正集中的数据进行中心化处理，得初始中心化样本数据；

c.对所述的初始中心化样本数据进行偏最小二乘(PLS)回归，得初始回归系数及偏最小二乘初始回归模型；

d.当第i个样本点加入该校正集时，对其进行中心化处理，并利用偏最小二乘初始回归模型获得相应的预测值；其中，i＝n+1,n+2，……；

e.计算所述的预测值与真实值之间的误差；若该误差小于等于阀值δ，则令i＝i+1，转至d；否则对回归系数进行更新后，令i＝i+1，转至d，直至得到最终的回归系数B_i，即得增量偏最小二乘模型。

步骤e中，通过以下方式对回归系数进行更新，得最终的回归系数B_i：

其中，i＝n+1,n+2，……，B_i为第i个样本点＜x_(i),y_i＞加入校正集时得到的更新的回归系数，B_i-1为对所述的初始校正集中的n个训练样本数据进行偏最小二乘法回归时获得的初始回归系数，η为正常数，其大小由沿着梯度反方向搜索时的步长决定；为第i个样本点的中心化数据，为第i个样本点的预测值，为的转置。

步骤e中，所述的增量偏最小二乘模型为：

y＝XB_i+E

其中，所述的y为待测样品中各成分的含量，X为待测样品的近红外光谱数据，B_i为增量偏最小二乘模型最终的回归系数，E为残差。

步骤e中所述的阀值δ采用K-折交叉验证的方式确定，具体包括以下步骤(如图29所示)：

(1)将初始校正集中的n个样本数据平均分为k份，以第j份作为增量训练集，剩余的k-1份作为初始训练集进行交叉验证，其中，1≤j≤k(k≥4)；

(2)对所述的初始训练集中的数据进行中心化处理，得初始中心化样本数据；

(3)对初始中心化样本数据进行PLS回归，得初始回归系数及PLS初始回归模型；

(4)当第j份作为增量训练集中的数据加入该初始训练集时，对增量训练集中的数据进行中心化处理，并利用PLS初始回归模型获得相应的预测值；

(5)设δ＝2^-i，1≤i≤n_top，其中n_top∈N；

(6)计算所述的预测值与真实值之间的误差；若该误差小于等于δ，则转至(7)；否则对所述的回归系数进行更新，并将该更新后的回归系数应用到作为增量训练集的第j份数据上，获得相应的预测值及交叉验证均方差MSE(k,i)；将所述的交叉验证均方差MSE(k,i)放入均方差矩阵中；

(7)令i＝i+1，若i＜n_top转至(5)；否则，则令j＝j+1，转至(1)；若j＞k，得k行n_top列的均方差矩阵；

(8)对所述的k行n_top列的均方差矩阵的每一列求均值，获得均值矩阵；

(9)查找均值矩阵中的均方差最小值，该均方差最小值在矩阵中所对应的列号为i_best，相应的即为最优阀值δ。

所述待测样品为谷物、土壤、草或复方中药药物，且对于复方中药药物芍药苷含量的测定尤为精确。

本发明中的回归系数的更新公式——的推导过程如下：

(1)设初始校正集的均值为计算公式如下：

(2)初始中心化样本数据可表示为：

(3)对初始中心化样本数据和进行PLS1回归，得初始回归系数B_n；

(4)当有新的样本点＜x_(n+1),y_n+1＞加入到初始校正集＜X⁽ⁿ⁾,Y⁽ⁿ⁾＞中时，新的校正集用＜X,Y＞表示，则新的均值可用为：

(5)对第n+1个新增样本点＜x_(n+1),y_n+1＞进行中心化处理，得：

(6)将代入所得的回归模型中，求得预测值

(7)计算所述的预测值与真实值之间的误差，若误差不在设定的阀值δ(为2^-i，i可取任意正整数)范围内，则对回归系数进行更新，此时将代入新的回归模型中，得新的预测值

(8)将预测值的改变量用表示：

因此，当增加一个新样本时，回归系数的更新过程可以表示为：

B_n+1←B_n+ΔB

(9)由于回归系数的更新思想是使得每个新增单独样本点的真实值与预测值的误差e达到最小，那么对于第n+1个样本点＜x_(n+1),y_n+1＞，其真实值与预测值的误差记为：

(10)采用最小二乘的思想进行处理，得：

(11)为了使E达到最小，计算E对于B_n的梯度，其梯度可用▽E_n+1(B_n)表示为：

(12)由于梯度的反方向可以使回归系数向E最陡峭下降的方向移动，因此，沿着梯度的反方向即可找到E的最小值；那么，ΔB可表示为：

ΔB＝-η▽E_n+1(B_n)

其中，η是一个正常数，它决定了沿着梯度反方向搜索时的步长；因此，公式B_n+1←B_n+ΔB还可表示为：

B_n+1←B_n-η▽E_n+1(B_n)

(13)将公式B_n+1←B_n-η▽E_n+1(B_n)表示成分量的形式：

将公式代入梯度向量的分量得：

(14)将公式代入公式得回归系数每个分量的更新公式：

因此，回归系数B_n+1的更新模型为：

(15)当第i个样本点＜x_(i),y_i＞，i＝n+2,n+3...加入到校正集＜X,Y＞中时，新的均值可用表示为：

(16)第i个样本点的中心化数据可表示为：

(17)第i个样本点的预测值为：

(18)计算预测值与之间的误差，若该误差不在设定的阀值δ范围内，则对回归系数进行更新；

(19)重复上述步骤继续处理下一个样本点，直到得到最终的回归系数向量。回归系数的更新过程如下所示：

.

.

.

.

.

.

对上述过程进行归纳，可得回归系数的更新公式为：

其中，i＝n+1,n+2,...。

实验例1：对含有芍药苷的复方中药药物进行测定

1.实验目的

(1)在中药煎煮中，能够确定合理的时间，得到芍药苷浓度满足要求的原药；

(2)验证增量偏最小二乘法(即IPLS)是否比传统偏最小二乘法(即PLS)更能有效的处理增量数据。

2.样品处方和工艺：

【处方】柴胡180g延胡索(炙)200g白芍240g炙甘草100g

【制法】以上四味加水煎煮二次，第一次8倍量，煎煮2小时，第二次6倍量煎煮1小时。

1.实验准备与取样规则：

(1)实验仪器及药材以10L的三颈瓶进行提取，每批次720g药材，有含量检查项的药材均准备不同含量的三个批次，用于完成6个批次的建模提取，1个批次进行验证采用含量不同的三个批次药材随机混合投料。

(2)样品数量样品数量按提取工艺时间安排取样间隔，均从沸腾开始取样，共取了7个批次，每个批次都煎煮两次，所以每个批次都是由一煎和二煎数据组成，共取样144个，其中第一批次得到18个样本数据(一煎12个，二煎6个)，其余六批次均得到21个样本数据(一煎15个，二煎6个)。本次实验采用其它近红外厂商实验后的存样进行，样品已经过筛网过滤和冰箱统一冷藏密封存放。

4.实验仪器：

硬件设备：Bruker公司生产的MATRIX-F型傅里叶变换近红外光谱仪，配有光纤长度2m的2mm固定光程手持式透射测量实验室光纤探头。

测试方法：光谱扫描范围(12000～4000cm-1)，光谱分辨率(8cm-1)，扫描次数(128次)，干涉仪扫描频率10kHz，双边扫描得到吸光度光谱；每个样品测量1次。

5.增量模型的建立：

按照处方药品配方分批次，对其进行煎煮，测量到的样本数据，按其到达时间顺序，用IPLS建立样品浓度y关于光谱X的模型，并将所得结果与PLS的结果进行比较。

建模方法：IPLS；对比方法：PLS。

预处理方法：中心化。

参数选择方法：通过K-fold和F检验确定主成分数PC，基于K-折交叉验的阀值确定方法选择δ。

建模的软件：Matlab软件。

实验环境：处理器：Intel(R)GMA4500MHD

内存：2.00GB

操作系统：window7

建模过程：

模型一：首先选取前两个批次的数据，使其作为初始校正集，得到初始校正模型，然后，用三，四批次的数据进行增量模型的修正，得到前四个批次的赠量模型；

模型二：将前四个批次的数据作为初始校正集数据，得到初始校正模型，用五，六批的数据进行增量模型的修正，得到前六个批次的增量模型。

5.1IPLS建立芍药苷浓度模型

本发明中的增量偏最小二乘(即IPLS)法采用了近似梯度下降的想法。该方法的原理是，对于初始校正集＜X⁽ⁿ⁾,y⁽ⁿ⁾＞，首先对其进行PLS1回归，得到初始回归系数B_n，当有新的样本点＜x_(i),y_i＞加入到校正集＜X⁽ⁿ⁾,y⁽ⁿ⁾＞中时，计算新增样本点x_(i)的预测值及预测误差，如果该误差不在设定的阀值δ范围内，则更新回归系数，否则继续处理下一个新增样本点。本方法具体过程如下；

输入：初始校正集＜X⁽ⁿ⁾,y⁽ⁿ⁾＞，潜变量数A

输出：回归系数B

1：

2：

3：fori＝n+1,n+2,...

4：

5：

6：

7：

8：

9：ifΔ＞δ

10：

本发明采用K-fold交叉验证确定阀值δ，具体方法如下：

将初始校正集＜X⁽ⁿ⁾,y⁽ⁿ⁾＞平均分为k份，每份所含样本数为n/k。轮流从这k份数据中抽取一份作为测试集X_test∈R^(n/k)×p和y_test∈R^(n/k)×1，记测试样本个数为n_test，剩下的k-1份作为训练集X_tran∈R^{(n×(n-1)/k)×p}和y_tran∈R^{(n×(n-1)/k)×1}。令δ＝2^-i，则对于不同的指数i＝1:n_top，其中n_top为指数i的上限，分别将训练集X_tran，y_tran作为初始训练集，对应的测试集X_test，y_test作为增量训练集代入算法2中进行IPLS1回归，得到回归模型B，然后用如下公式预测X_test，得预测值

用如下公式计算交叉验证的均方差MSE：

其中，y_i为y_test的第i个分量，为的第i个分量，计算不同的指数i下均方差MSE的均值mMSE，找出最小的mMSE，记为minMSE，其对应的下标为阀值的最佳指数i_best，对应的阀值为最佳阀值。

用K-fold交叉验证确定阀值δ的具体步骤如下：

输入：初始校正样本集＜X⁽ⁿ⁾,y⁽ⁿ⁾＞，潜变量数A，学习率η，阀值的指数上限n

输出：阀值最佳指数i_best，阀值δ

1：计算校验集＜X⁽ⁿ⁾,y⁽ⁿ⁾＞的样本数，设为n；

2：将校验集随机分为SEG_calib(≥4)份；

3：fork＝1:SEG_calib

4：第k份代表测试集X_test，y_test，样本数为n_test；

5：剩余的为训练集X_tran，y_tran，样本数为n-n_test；

6：fori＝1:n_top

7：δ＝2^-i

8：B＝IPLS1(X_tran,y_tran,A,δ,η,X_test,y_test)

9：将回归模型应用到测试集X_test上，得到交叉验证预测值；

10：计算第k份数据在不同指数i下的均方差MSE(k,i)；

11：求得不同指数i下的均方差MSE(k,i)的均值mMSE；

12：找出最小的mMSE，记为minMSE及对应的下标i_best；

13：i_best为最终阀值指数，最佳阀值作为结果返回。

5.1.1前四个批次联合建模(IPLS1-C1,PLS1-C1)

首先提取前四个批次的样本数据，按提取工艺时间安排取样间隔，均从沸腾开始取样，每个批次数据都由一煎和二煎数据组成。选择一、二批次39个样本数据作为初始训练样本集；三、四批次42个样本数据作为增量训练样本集。建立前四个批次的增量模型，这里我们用单因变量增量偏最小二乘(IPLS1)所得结果与单因变量偏最小二乘(PLS1)的结果进行比较。

IPLS1建模结果：lv数量＝3，δ＝2^-10，RMSECV＝0.0662，RMSEC＝0.0748(如表1所示)。

PLS1与IPLS1回归系数比较：前四个批次联合建模的回归系数如图1所示，图1中共有五条曲线，其中黑色实线代表PLS1所得回归系数，从图中可以看出，IPLS1算法的回归系数一共更新了三十次，图中仅表示了IPLS1第一次，第十次，第二十次，第三十次回归系数的更新模型，可知，尽管两种方法对回归系数的更新方法不同，两者所得的回归系数走势却相同。

阀值与交叉验证预测均方根误差的关系：对中药前二个批次的初始校正集实施K-fold交叉验证的阀值确定方法，设阀值指数的取值上限为10，通过阀值确定方法，得到该数据的最佳阀值，如图2所示，可知，最佳阀值为2^-10。

5.1.2前六个批次联合建模(IPLS1-C2,PLS1-C2)

继续提取第五、六批次的样本数据，按提取工艺时间安排取样间隔，同样从沸腾开始取样，将已测好的前四个批次的81个样本数据作为初始训练样本集；五、六批次42样本数据个作为增量训练样本集，建立前六个批次的增量模型，所得结果与PLS1的结果进行比较。

IPLS1建模结果：lv数量＝4，δ＝2^-9，RMSECV＝0.0760，RMSEC＝0.0761(如表1所示)。

PLS1与IPLS1回归系数比较：前六个批次联合建模回归系数比较示意图如图4所示，图4中共有四条曲线，其中黑色实线代表PLS1所得回归系数，从图中可以看出，IPLS1算法的回归系数一共更新了二十九次，图中仅表示了IPLS1第一次，第十九次，第二十九次回归系数的更新模型，并将波数在800到900的回归模型进行放大。可知，尽管两种方法对回归系数的更新方法不同，但两者所得的回归系数走势基本相同，同时可知随着模型不断的更新，IPLS1所得的回归系数越来越接近PLS1的回归系数。其中前六个批次123条近红外光谱如图3所示。

阀值与交叉验证预测均方根误差的关系：对中药前四个批次的初始校正集实施K-fold交叉验证的阀值确定方法，设阀值指数的取值上限为10，通过阀值确定方法，得到该数据的最佳阀值，如图5所示，可知，最佳阀值为2^-10。

5.2用IPLS预测后续批次芍药苷浓度

(1)对于前四个批次建立的增量模型，我们选择五，六批次的数据作为测试数据，得到预测结果；

(2)对于前六个批次建立的增量模型，我们选择第七个批次的数据作为测试数据，得到预测结果。

结果比较：对比本方法与使用相同数据集的PLS建模方法的精度与其他性能指标。

5.2.1五、六批次数据的预测分析(IPLS1-P1,PLS1-P1)

预测：用前四个批次已建好的模型对第五、六批次样本的浓度进行预测(如表1所示)：

预测结果：IPLS1:RMSEP＝0.0624，SEP＝0.0643

PLS1:RMSEP＝0.0653，SEP＝0.0654

PLS1与IPLS1测试集真实值与预测值的比较：由预测结果和图6可知，测试集的真实值与预测值的数值对均在直线y＝x附近，并向直线y＝x靠拢，IPLS1几乎所有的数值对都比PLS1的数值对更靠近直线y＝x，说明IPLS1预测精度相比PLS1得到了提高。

5.2.2第七个批次数据预测分析(IPLS1-P2,PLS1-P2)

预测：用前六个批次已建好的模型对第七个批次样本的浓度进行预测(如表1所示)：

预测结果：IPLS1：RMSEP＝0.0465，SEP＝0.0472

PLS1：RMSEP＝0.0512，SEP＝0.0524

PLS1与IPLS1测试集真实值与预测值的比较：由预测结果和图7可知，测试集的真实值与预测值的数值对均在直线y＝x附近，并向直线y＝x靠拢，IPLS1几乎所有的数值对都比PLS1的数值对更靠近直线y＝x，说明IPLS1预测精度相比PLS1得到了提高。

5.3预测结果的比较与分析

如表1所示，RMSECV为交叉验证均方根误差；RMSEC为校正均方根误差，由训练样本得到；RMSEP为预测均方根误差，由测试样本得到；SEP为预测标准误差；h_improve为IPLS1对模型预测精度的改进程度；No.为训练所用样本点数量。

表1不同模型的预测结果比较

具体的说，RMSEP为模型的均方根误差，通过模型的均方根误差来衡量模型的预测精度，RMSEP的计算公式为：

其中，y_pre为待测组分的含量预测值，y_ref为待测组分的含量真实值，N为测试样本数。

SEP为预测标准误差(standarderrorofprediction)，SEP的计算公式为：

其中，bias为平均偏差，n_tst为预测样本点数量，y_pre为待测组分的含量预测值，y_ref为待测组分的含量真实值。

h_improve为IPLS1对模型预测精度的改进程度，h_improve的计算公式为：

其中，RMSEP_IPLS、RMSEP_PLS分别代表IPLS1和PLS1方法的预测均方根误差。

time为PLS1算法及IPLS1算法为对初始校正集的训练时间，随着新增样本点的加入，所述的time还包括对回归系数更新的时间以及用最终的回归系数对预测集进行预测的时间，具体计算公式如下：

time＝time(plsr)+time(B_update)+time(predict)

由表1可知：IPLS1所得的RMSEP均比PLS1所得的RMSEP要小，可知，IPLS1的预测精度相比PLS1都有所提高；另外，从这两种方法所用的时间来看，IPLS1算法相比传统的PLS1算法节约了更多的时间；此外，从表1中还可以看出，前六个批次的建模效果及预测效果都比前四个批次的建模效果及预测效果要好，说明随着建模数据的增多，本发明所建立的IPLS模型越来越准确，预测精度也越来越高。

6实验结论

实验结果表明，PLS1算法及IPLS1算法对回归系数的更新过程不同，二者所获得的回归系数基本相同，但IPLS1算法的预测均方根误差却更小且运行时间更短，可见，本发明中的增量偏最小二乘模型具有更高的预测精度和建模效率；此外，本发明中的IPLS1算法通过确定最佳阀值，从而进一步降低了IPLS1模型的预测误差。

实验例2：通过谷物，土壤及草样本数据集的仿真实验验证本发明方法的有效性

1.实验数据来源

谷物数据集由EigenvectorResearch公司提供。该数据含有80个谷物样本的近红外光谱和对应的物质含量(水分，油脂，蛋白质和淀粉)。光谱采集范围为1100-2498nm，采集间隔为2nm。

在试验中，使用Kennard-Stone(KS)方法将数据集分为60个训练样本和20个测试样本。其中训练集同样使用KS抽取30个作为IPLS1初始训练集，30个作为其增量训练集。

土壤数据集为土壤样本中有机质含量的吸光率。土壤样本来自于位于瑞典北部的阿比斯库市的长期田间试验。该数据含有108个土壤样本，光谱采集范围为400-2498nm，采集间隔为2nm。在试验中，使用Kennard-Stone(KS)方法将数据集分为80个训练样本和28个测试样本。其中训练集同样使用KS抽取40个作为IPLS1初始训练集，40个作为其增量训练集。

草数据集来自于第98届钱伯斯堡会议国际开发研究中心，该数据集包含141个羊茅草粉末样品近红外反射光谱数据，羊茅草粉末样品包含规定的碳、氮、硫成分的范围分别为29.6％～40.9％，1.1％～6.6％，0.3％～1.7％。相关的化学值是未知平行样品的均值，该平行样品由力可CNS-2000碳、氮和硫分析仪测得。在试验中，使用Kennard-Stone(KS)方法将数据集分为100个训练样本和41个测试样本。其中训练集同样使用KS抽取50个作为IPLS1初始训练集，50个作为其增量训练集。

2.处理方法

对于相同的训练集和测试集，分别采用PLS1方法和本发明的IPLS1方法处理上述样本数据。两种方法潜变量数的选择方式都是通过K-折交叉验证和F检验来确定的，IPLS1方法阀值参数的确定如实施例1。

首先比较两种方法对于上述样本数据所得到的回归系数，观察通过不同方法对回归系数更新，二者所得的回归系数是否存在显著差异；其次，通过计算模型的均方根误差来衡量模型的预测精度，以及IPLS1对模型预测精度的改进程度。

其中，谷物数据的近红外光谱如图8所示；土壤数据的近红外光谱如图9所示；草数据的近红外光谱如图10所示；光谱图的纵坐标表示样本的吸光率，横坐标表示光谱的波长。

阀值δ的选取：采用本发明中的方法进行确定。

3.实验结果与讨论

3.1回归系数的比较

用PLS1与IPLS1法分别处理谷物、土壤和草样本数据，如图11～图15所示，红色曲线代表PLS1算法得到的回归模型，其他颜色的曲线代表IPLS1算法的回归模型。

图11为corn(moisture)的回归系数对比图，图11中一共有四条曲线，除黑色实线外的曲线表示IPLS1算法所得的回归系数的更新过程，从图11中可以看出，IPLS1算法的回归系数更新了六次，图11中仅画出IPLS1第一次，第三次，第六次的回归系数。可知，尽管两种方法对回归系数的更新方法不同，但两者所得的回归系数却基本相同。

图12为corn(oil)的回归系数对比示意图，图12中一共有四条曲线，除黑色实线外的曲线表示IPLS1算法所得的回归系数的更新过程，从图中可以看出，IPLS1算法的回归系数一共更新了十四次，图中仅画出IPLS1第一次，第七次，第十四次的回归系数。可见，尽管两种方法对回归系数的更新方法不同，但两者所得的回归系数却基本相同。

图13是corn(protein)的回归系数对比示意图；图13中一共有五条曲线附图说明，除黑色实线外的曲线表示IPLS1算法所得的回归系数的更新过程，从图中可以看出，IPLS1算法的回归系数一共更新了二十四次，图中仅画出IPLS1第一次，第六次，第十二次，第二十四次的回归系数。可知，尽管两种方法对回归系数的更新方法不同，但两者所得的回归系数却基本相同。

图14是corn(starch)的回归系数对比图；图14中一共有四条曲线，除黑色实线外的曲线表示IPLS1算法所得的回归系数的更新过程，从图中可以看出，IPLS1算法的回归系数一共更新了二十四次，图中仅画出IPLS1第一次，第十二次，第二十次的回归系数。可知，尽管两种方法对回归系数的更新方法不同，但两者所得的回归系数却基本相同。

图15是grass的回归系数对比示意图；图15中一共有七条曲线，除黑色实线外的曲线表示IPLS1算法所得的回归系数的更新过程，从图中可以看出，IPLS1算法的回归系数一共更新了五十次，图中仅画出IPLS1第一次，第十次，第二十次，第三十次，第四十次，第五十次的回归系数。可知，尽管两种方法对回归系数的更新方法不同，但两者所得的回归系数却基本相同。

3.2阀值δ的选取：

分别对谷物、草和土壤的初始校正集实施K-fold交叉验证的阀值确定方法，设阀值指数的最大取值为10，通过阀值确定方法，得到每组数据的最佳阀值，如图16～21所示，可知，谷物中四种成分的最佳阀值分别为2^-10、2^-8、2^-8和2^-5，草的最佳阀值为2^-10，土壤的最佳阀值为2^-6。通过本发明的方法所得的阀值，对应的IPLS模型的预测均方根误差最小，可知，通过本发明的方法选出来的阀值可以进一步降低模型的预测误差，提高模型的预测精度。

3.3预测结果的比较与分析

分别采用PLS1和IPLS1法对上述样品数据集进行处理，获得的各个样品成分含量的预测结果如表2所示：

表2PLS1与IPLS1预测结果比较

由表2可知：对于不同的数据集，IPLS1所得的RMSEP均比PLS1所得的RMSEP要小，可见，IPLS1的预测精度相比PLS1均有所提高；同时从两种方法所用的时间来看，IPLS1算法相比传统的PLS1算法节约了更多的时间，因而建模效率也更高。

另外，从图22～图27中可以看出，两种方法所得的预测值与真实值的比较点均在直线y＝x周围，并向直线y＝x靠拢，说明IPLS1算法对谷物、草和土壤这三组数据的预测效果较好。其中：

图22是corn(moisture)预测值与真实值的比较结果示意图；

图23是corn(oil)预测值与真实值的比较结果示意图；

图24是corn(protein)预测值与真实值的比较结果示意图；

图25是corn(starch)预测值与真实值的比较结果示意图；

图26是grass预测值与真实值的比较结果示意图；

图27是soil预测值与真实值的比较结果示意图。

4实验结论

通过PLS1算法和IPLS1算法分别对谷物、土壤和羊茅草三种数据建立校正模型，并用每组数据的测试集进行测试，求得预测误差。实验结果表明，虽然两种方法对回归系数的更新过程不同，但是，IPLS1算法与传统的PLS1算法获得的回归系数基本相同，且IPLS1算法的预测均方根误差却更小且运行时间更短，因此，本发明中的增量偏最小二乘模型具有更高的预测精度和建模效率；此外，本发明中的IPLS1算法通过确定最佳阀值，从而进一步降低了IPLS1模型的预测误差。可见，本发明提出的增量偏最小二乘算法可以更好的进行增量光谱数据的定量分析，在时间和空间上提高了建模的效率。

Claims (8)

1.一种基于增量偏最小二乘法的样品成份测定方法，其特征在于，包括以下步骤：

S1，采集待测样品的近红外光谱数据；

S2，通过增量偏最小二乘模型，获得该近红外光谱数据所对应的样品中各成分的含量。

2.根据权利要求1所述的基于增量偏最小二乘法的样品成份测定方法，其特征在于，步骤S2中所述的增量偏最小二乘模型通过以下方法建立：

a.收集n个样本的近红外光谱数据及其对应的样品中各成分的含量数据作为初始校正集；

b.将初始校正集中的数据进行中心化处理，得初始中心化样本数据；

c.对所述的初始中心化样本数据进行偏最小二乘回归，得初始回归系数及偏最小二乘初始回归模型；

d.当第i个样本点加入该校正集时，对其进行中心化处理，并利用偏最小二乘初始回归模型获得相应的预测值；其中，i＝n+1,n+2，……；

e.计算所述的预测值与真实值之间的误差；若该误差小于等于阀值δ，则令i＝i+1，转至d；否则对回归系数进行更新后，令i＝i+1，转至d，直至得到最终的回归系数B_i，即得增量偏最小二乘模型。

3.根据权利要求2所述的基于增量偏最小二乘法的样品成份测定方法，其特征在于，步骤e中，通过以下方式对回归系数进行更新，得最终的回归系数B_i：

$B_i=B_{i-1}+\mathrm{\η}(y_i^{cen}-{\hat{y}}_i)x_{\left(i\right)}^{cenT}$

其中，i＝n+1,n+2，……，B_i为第i个样本点＜x_(i),y_i＞加入校正集时得到的更新的回归系数，B_i-1为对所述的初始校正集中的n个训练样本数据进行偏最小二乘法回归时获得的初始回归系数，η为正常数，其大小由沿着梯度反方向搜索时的步长决定；为第i个样本点的中心化数据，为第i个样本点的预测值，为的转置。

4.根据权利要求3所述的基于增量偏最小二乘法的样品成份测定方法，其特征在于步骤e中，所述的增量偏最小二乘模型为：

y＝XB_i+E

其中，所述的y为待测样品中各成分的含量，X为待测样品的近红外光谱数据，B_i为增量偏最小二乘模型最终的回归系数，E为残差。

5.根据权利要求2～4任一项所述的基于增量偏最小二乘法的样品成份测定方法，其特征在于，步骤e中所述的阀值δ采用K-折交叉验证的方式确定。

6.根据权利要求5所述的基于增量偏最小二乘法的样品成份测定方法，其特征在于，具体通过以下方法确定阀值δ：

(1)将初始校正集中的n个样本数据平均分为k份，以第j份作为增量训练集，剩余的k-1份作为初始训练集进行交叉验证，其中，1≤j≤k(k≥4)；

(2)对所述的初始训练集中的数据进行中心化处理，得初始中心化样本数据；

(3)对初始中心化样本数据进行PLS回归，得初始回归系数及PLS初始回归模型；

(4)当第j份作为增量训练集中的数据加入该初始训练集时，对增量训练集中的数据进行中心化处理，并利用PLS初始回归模型获得相应的预测值；

(5)设δ＝2^-i，1≤i≤n_top，其中n_top∈N；

(6)计算所述的预测值与真实值之间的误差；若该误差小于等于δ，则转至(7)；否则对所述的回归系数进行更新，并将该更新后的回归系数应用到作为增量训练集的第j份数据上，获得相应的预测值及交叉验证均方差MSE(k,i)；将所述的交叉验证均方差MSE(k,i)放入均方差矩阵中；

(7)令i＝i+1，若i＜n_top转至(5)；否则，则令j＝j+1，转至(1)；若j＞k，得k行n_top列的均方差矩阵；

(8)对所述的k行n_top列的均方差矩阵的每一列求均值，获得均值矩阵；

(9)查找均值矩阵中的均方差最小值，该均方差最小值在矩阵中所对应的列号为i_best，相应的即为最优阀值δ。

7.根据权利要求1或6所述的基于增量偏最小二乘法的样品成份测定方法，其特征在于，所述的待测样品为谷物、土壤、草或复方中药药物。

8.根据权利要求7所述的基于增量偏最小二乘法的样品成份测定方法，其特征在于，所述的待测样品为含有芍药苷的复方中药药物。