本发明属于近红外光谱技术在线分析领域,具体涉及一种基于GT-KF-PLS(Gamma test-Kalman filter-Partial least square,伽马测试的卡尔曼滤波偏最小二乘)近红外光谱动态演化模型校正方法及系统。
背景技术:
KF-PLS(Kalman FilterPartial Least Squares)已被证明是提高近红外光谱多元校正模型适应性的一种新方法。采用该方法建立的校正模型具有适应设备陈旧化、环境变化和模型界外样品等优点。
但是,KF-PLS噪声方差方法难以得到观测噪声不确定的动态近红外光谱仪测量过程的噪声方差值,因而常将观测噪声方差值置零,并引入了遗忘因子以进行KF-PLS模型校正。这在一定程度上虽然抑制了模型的不确定性发生,但是由于校正模型误差累积作用,在被测样品数量不断增加的情况下,校正模型不可避免地会出现发散现象,甚至完全失效。极大地限制了KF-PLS在观测噪声不确定近红外光谱模型校正中的应用。因此,如何保证自适应校正模型的稳定性成为了难点。
因此,为解决上述问题,本发明提出了一种基于GT-KF-PLS近红外光谱自适应模型校正方法及系统。
技术实现要素:
为了实现上述目的,本发明提供一种基于GT-KF-PLS(Gamma test-Kalman filter-Partial least square,伽马测试的卡尔曼滤波偏最小二乘)近红外光谱自适应模型校正方法及系统,以解决观测噪声不确定引起的KF-PLS校正模型发散问题。
本发明提供一种基于GT-KF-PLS近红外光谱自适应模型校正方法,包括:
S1:利用K/S算法从标准样品中选择有代表性的建模样品;
S2:采用PLS法对所述建模样品建立近红外光谱数据与浓度间的线性关系,并根据所述线性关系建立PLS校正模型;
S3:利用所述PLS校正模型对待测样品进行预测,获取所述待测样品的预测值;
S4:同时,定期对所述待测样品进行化验,并采集与化验的待测样品对应的待测样品的样品数据;
S5:采用Gamma Test对样本光谱数据和化验后采集到的待测样品的样品数据进行噪声统计值的计算,获取系统噪声的精确噪声方差值;
S6:根据所述精确噪声方差值、所述待测样品的样品数据以及当前时刻所述待测样品的预测值,通过KF算法修正当前时刻所述PLS校正模型的主因子系数。
此外,优选的方案是,在步骤S2中,在采用PLS法对所述建模样品建立近红外光谱数据与浓度间的线性关系的过程中,
设Xm×n为m个样品在n个波长上的光谱参数矩阵,Ym×p为m个样品p种成分含量构成的浓度矩阵,将Xm×n与Ym×p分解为如下形式:
X=TPt+E
Y=UQt+F
其中,矩阵T和矩阵U分别表示去掉大部分噪声后的光谱信息和浓度信息;E和F表示误差;
由于Xm×n与Ym×p存在线性关系Y=P'X,在分解时,矩阵T和U之间的线性关系为:U=TB;通过交换迭代矢量而使两个分解过程合二为一。
此外,优选的方案是,在步骤S5中,在采用Gamma Test对样本光谱数据和化验后采集到的待测样品的样品数据进行噪声统计值的计算,获取系统噪声的精确噪声方差值的过程中,获取第i个样本点对应的系统噪声方差值包括:
S51:假定数据之间的关系:Y=h(X)+r,
其中,h(X)表示光滑函数;r表示噪声变量;
S52:使用kd-tree算法在输入空间对各输入样本点Xi(1≤i≤M)进行计算,得到输入样本Xi(1≤i≤M)的第K(1≤K≤P)近邻域点XN[i,K](1≤i≤M),
S53:计算所有Xi(1≤i≤M)的第P近邻域点的最小均方距离δ(K)以及输出空间相应的最小均方距离γ(K),
S54:对(δ(K),γ(K))K(1≤K≤P)、(δ(K),γ(K))K(1≤K≤P),按公式γ=Aδ+R进行一次线性回归,所得一次线性函数的截距,即系统噪声方差值R;
S55:当新增加一个标准样品时,重复步骤S51至步骤S54,,得到每个样品对应噪声方差。
此外,优选的方案是,在步骤S6中,根据所述精确噪声方差值、所述待测样品的样品数据以及当前时刻所述待测样品的预测值,通过KF算法修正当前时刻所述PLS校正模型的主因子系数过程中,
设PLS初始模型主因子数为l,主因子系数为:
w1,t1,v1,p1;w2,t2,v2,p2;……;wi,ti,vi,pi(i=1,2,3…,l);
其中:
vi=(tTy)/(tTt)=[vi1 vi2 ... vip]
将所述PLS初始模型中的所有系数值组成状态向量:
W=[w1Tt1Tv1p1T...wiTtiTvipiT]T(i=1,2,3…,l)
系统的状态方程和观测方程表示为:
其中,Yek为标样浓度;Wk为第k个标样修正时刻的主因子系数;Xk为第k个样品光谱矢量;Yrk为预测浓度;
Vk为观测噪声,其统计特性为:
令
则观测方程为:Yek=HkWk+Dk+Vk;
其中,H表示状态变量Wk对测量变量Yk的增益;
Wk表示k时刻状态变量,即:用第k个标样修正时的PLS主因子系数;
Dk为中间变量。
本发明还提供一种基于GT-KF-PLS近红外光谱自适应模型校正系统,包括:建模样品选取单元,用于利用K/S算法从标准样品中选择有代表性的建模样品;
PLS校正模型建立单元,用于采用PLS法对所述建模样品建立近红外光谱数据与浓度间的线性关系,并根据所述线性关系建立PLS校正模型;
预测值获取单元,用于利用所述PLS校正模型对待测样品进行预测,获取所述待测样品的预测值;
样品数据获取单元,用于定期对所述待测样品进行化验,并采集与化验的待测样品对应的待测样品的样品数据;
精确噪声方差值获取单元,用于采用Gamma Test对样本光谱数据和化验后采集到的待测样品的样品数据进行噪声统计值的计算,获取系统噪声的精确噪声方差值;
PLS校正模型的主因子系数修正单元,用于根据所述精确噪声方差值、所述待测样品的样品数据以及当前时刻所述待测样品的预测值,通过KF算法修正当前时刻所述PLS校正模型的主因子系数。
此外,优选的方案是,所述PLS校正模型建立单元在采用PLS法对所述建模样品建立近红外光谱数据与浓度间的线性关系的过程中,
设Xm×n为m个样品在n个波长上的光谱参数矩阵,Ym×p为m个样品p种成分含量构成的浓度矩阵,将Xm×n与Ym×p分解为如下形式:
X=TPt+E
Y=UQt+F
其中,矩阵T和矩阵U分别表示去掉大部分噪声后的光谱信息和浓度信息;E和F表示误差;
由于Xm×n与Ym×p存在线性关系Y=P'X,在分解时,矩阵T和U之间的线性关系为U=TB;通过交换迭代矢量而使两个分解过程合二为一。
此外,优选的方案是,所述精确噪声方差值获取单元在采用Gamma Test对样本光谱数据和化验后采集到的待测样品的样品数据进行噪声统计值的计算,获取系统噪声的精确噪声方差值的过程中,获取第i个样本点对应的系统噪声方差值包括:
S51:假定数据之间的关系:Y=h(X)+r,
其中,h(X)表示光滑函数;r表示噪声变量;
S52:使用kd-tree算法在输入空间对各输入样本点Xi(1≤i≤M)进行计算,得到输入样本Xi(1≤i≤M)的第K(1≤K≤P)近邻域点XN[i,K](1≤i≤M),
S53:计算所有Xi(1≤i≤M)的第P近邻域点的最小均方距离δ(K)以及输出空间相应的最小均方距离γ(K),
S54:对(δ(K),γ(K))K(1≤K≤P)、(δ(K),γ(K))K(1≤K≤P),按公式γ=Aδ+R进行一次线性回归,所得一次线性函数的截距,即系统噪声方差值R;
S55:当新增加一个标准样品时,重复步骤S51至步骤S54,,得到每个样品对应噪声方差。
此外,优选的方案是,所述PLS校正模型的主因子系数修正单元根据所述精确噪声方差值、所述待测样品的样品数据以及当前时刻所述待测样品的预测值,通过KF算法修正当前时刻所述PLS校正模型的主因子系数过程中,
设PLS初始模型主因子数为l,主因子系数为:
w1,t1,v1,p1;w2,t2,v2,p2;……;wi,ti,vi,pi(i=1,2,3…,l);
其中:
vi=(tTy)/(tTt)=[vi1 vi2 ... vip]
将所述PLS初始模型中的所有系数值组成状态向量:
W=[w1Tt1Tv1p1T...wiTtiTvipiT]T(i=1,2,3…,l)
系统的状态方程和观测方程表示为:
其中,Yek为标样浓度;Wk为第k个标样修正时刻的主因子系数;Xk为第k个样品光谱矢量;Yrk为预测浓度;
Vk为观测噪声,其统计特性为:
令
则观测方程为:Yek=HkWk+Dk+Vk;
其中,H表示状态变量Wk对测量变量Yk的增益;
Wk表示k时刻状态变量,即:用第k个标样修正时的PLS主因子系数;
Dk为中间变量。
从上面的技术方案可知,本发明提供的基于GT-KF-PLS近红外光谱自适应模型校正方法及系统,有效利用观测输入(样本近红外光谱)输出数据(样本化验值),提出样本有效噪声方差(Gamma test,GT)改进的KF-PLS模型校正方法;采用衰减记忆的GT对输入输出数据进行实时方差估计,得到准确的观测噪声方差值,再利用KF-PLS实现精确模型校正,能够保证近红外光谱自适应校正模型的稳定性,最终实现基于近红外光谱技术的在线分析。
附图说明
通过参考以下结合附图的说明及权利要求书的内容,并且随着对本发明的更全面理解,本发明的其它目的及结果将更加明白及易于理解。在附图中:
图1为根据本发明实施例的基于GT-KF-PLS近红外光谱自适应模型校正方法流程示意图;
图2为根据本发明实施例的基于GT-KF-PLS近红外光谱自适应模型校正系统结构框图。
具体实施方式
在下面的描述中,出于说明的目的,为了提供对一个或多个实施例的全面理解,阐述了许多具体细节。然而,很明显,也可以在没有这些具体细节的情况下实现这些实施例。
以下将结合附图对本发明的具体实施例进行详细描述。
为了说明本发明提供的基于GT-KF-PLS近红外光谱自适应模型校正方法,图1示出了根据本发明实施例的基于GT-KF-PLS近红外光谱自适应模型校正方法流程。
如图1所示,本发明提供的基于GT-KF-PLS近红外光谱自适应模型校正方法包括:
S1:利用K/S算法从标准样品中选择有代表性的建模样品;
S2:采用PLS法对所述建模样品建立近红外光谱数据与浓度间的线性关系,并根据所述线性关系建立PLS校正模型;
S3:利用所述PLS校正模型对待测样品进行预测,获取所述待测样品的预测值;
S4:同时,定期对所述待测样品进行化验,并采集与化验的待测样品对应的待测样品的样品数据;
S5:采用Gamma Test对样本光谱数据和化验后采集到的待测样品的样品数据进行噪声统计值的计算,获取系统噪声的精确噪声方差值;
S6:根据所述精确噪声方差值、所述待测样品的样品数据以及当前时刻所述待测样品的预测值,通过KF算法修正当前时刻所述PLS校正模型的主因子系数。
在进行步骤S1之前,首先对原始光谱数据进行矢量归一化预处理。
在步骤S1中,采用K/S(Kennard-Stone)算法从标准样品中选择有代表性的校正集样品,并按留一交叉验证的方法确定最终的主因子数。
在步骤S2中,采用PLS法建立近红外光谱数据与浓度间的线性关系。
设Xm×n为m个样品在n个波长上的光谱参数矩阵,Ym×p为m个样品p种成分含量构成的浓度矩阵,PLS法不直接建立每种成份与光谱参数向量的关系方程,而是考虑Xm×n与Ym×p的外部关系和联系二者的内部关系,将Xm×n与Ym×p分解为如下形式:
X=TPt+E
Y=UQt+F
矩阵T和矩阵U分别表示去掉大部分噪声后的光谱信息和浓度信息;E和F表示误差;
由于Xm×n与Ym×p存在线性关系Y=P'X,在分解时,还考虑矩阵T和U之间的线性关系为:U=TB;通过交换迭代矢量而使两个分解过程合二为一。
在步骤S5中,采用Gamma Test对已知样本光谱和样本化验值数据进行噪声统计值的计算,得到系统噪声的精确信息,其中第i个样本点对应的系统噪声方差计算过程如下:
S51:假定样品光谱与标准浓度之间的关系如下Y=h(X)+r,式中h(X)表示光滑函数;r表示噪声变量。
S52:首先使用kd-tree算法在输入空间对各输入样本点Xi(1≤i≤M)进行计算,得到输入样本Xi(1≤i≤M)的第K(1≤K≤P)近邻域点XN[i,K](1≤i≤M),
S53:计算所有Xi(1≤i≤M)的第P近邻域点的最小均方距离δ(K)以及输出空间相应的最小均方距离γ(K),
S54:最后,对(δ(K),γ(K))K(1≤K≤P)、(δ(K),γ(K))K(1≤K≤P),按公式γ=Aδ+R进行一次线性回归,所得一次线性函数的截距,即为gamma统计值,也即系统噪声方差值R;
S55:当新增加一个标准样品时,重复S51至S54的操作。可得到每个样品对应噪声方差。
在步骤S6中,基于KF-PLS近红外光谱与化验值的模型校正。
需要说明的是,KF因其有很好的自适应动态实时滤波能力,对PLS校正模型进行调整建模,可以有效地对动态系统模型进行子空间演化逼近,获得动态演化校正模型。利用KF对PLS主因子系数进行估计,将PLS模型主因子系数看作系统状态变量,标准样品待测化学值看作系统观测变量,这样就把问题转化为状态参数的估计问题;即可获得随环境噪声、设备老化和测量对象变化的动态演化校正模型,这种子空间逼近模型可准确反映系统的动态时变特性。
设PLS初始模型主因子数为l,主因子系数为:
w1,t1,v1,p1;w2,t2,v2,p2;……;wi,ti,vi,pi(i=1,2,3…,l);
其中:
vi=(tTy)/(tTt)=[vi1 vi2 ... vip]
为使上述问题中主因子系数的计算转化为滤波递推估计形式,将模型中的所有系数值组成状态向量:
W=[w1Tt1Tv1p1T...wiTtiTvipiT]T(i=1,2,3…,l)
系统的状态方程和观测方程表示为:
其中,Yek为标样浓度,Wk为第k个标样修正时刻的主因子系数,Xk为第k个样品光谱矢量,Yrk为预测浓度。Vk为观测噪声,其统计特性为:
令
则观测方程为:Yek=HkWk+Dk+Vk。
其中,H表示状态变量Wk对测量变量Yk的增益。
Wk:表示k时刻状态变量,在这里也即是用第k个标样修正时的PLS主因子系数。Dk没有特别的含义,只是个中间变量而已。这里是令所以则
在系统的精确噪声方差值的基础上,可以采用基于KF-PLS进行精确建模。在仪器使用过程中,新增一个标样,修正一次模型,并引入遗忘因子,逐渐遗忘陈旧样品的作用,使校正模型具有适应设备陈旧化、环境变化和模型界外样品的自适应性。
与上述方法相对应,本发明还提供一种基于GT-KF-PLS近红外光谱自适应模型校正系统,图2示出了根据本发明实施例的基于GT-KF-PLS近红外光谱自适应模型校正系统逻辑结构。
如图2所示,本发明提供的基于GT-KF-PLS近红外光谱自适应模型校正系统200,包括:建模样品选取单元210、PLS校正模型建立单元220、预测值获取单元230、样品数据获取单元240、精确噪声方差值获取单元250和PLS校正模型的主因子系数修正单元260。
具体地,建模样品选取单元210,用于利用K/S算法从标准样品中选择有代表性的建模样品;
PLS校正模型建立单元220,用于采用PLS法对所述建模样品建立近红外光谱数据与浓度间的线性关系,并根据所述线性关系建立PLS校正模型;
预测值获取单元230,用于利用所述PLS校正模型对待测样品进行预测,获取所述待测样品的预测值;
样品数据获取单元240,用于定期对所述待测样品进行化验,并采集与化验的待测样品对应的待测样品的样品数据;
精确噪声方差值获取单元250,用于采用Gamma Test对样本光谱数据和化验后采集到的待测样品的样品数据进行噪声统计值的计算,获取系统噪声的精确噪声方差值;
PLS校正模型的主因子系数修正单元260,用于根据所述精确噪声方差值、所述待测样品的样品数据以及当前时刻所述待测样品的预测值,通过KF算法修正当前时刻所述PLS校正模型的主因子系数。
其中,PLS校正模型建立单元220在采用PLS法对所述建模样品建立近红外光谱数据与浓度间的线性关系的过程中,
设Xm×n为m个样品在n个波长上的光谱参数矩阵,Ym×p为m个样品p种成分含量构成的浓度矩阵,将Xm×n与Ym×p分解为如下形式:
X=TPt+E
Y=UQt+F
其中,矩阵T和矩阵U分别表示去掉大部分噪声后的光谱信息和浓度信息;E和F表示误差;
由于Xm×n与Ym×p存在线性关系Y=P'X,在分解时,矩阵T和U之间的线性关系为:U=TB;通过交换迭代矢量而使两个分解过程合二为一。
其中,所述精确噪声方差值获取单元250,在采用Gamma Test对样本光谱数据和化验后采集到的待测样品的样品数据进行噪声统计值的计算,获取系统噪声的精确噪声方差值的过程中,
获取第i个样本点对应的系统噪声方差值包括:
S51:假定数据之间的关系:Y=h(X)+r,
其中,h(X)表示光滑函数;r表示噪声变量;
S52:使用kd-tree算法在输入空间对各输入样本点Xi(1≤i≤M)进行计算,得到输入样本Xi(1≤i≤M)的第K(1≤K≤P)近邻域点XN[i,K](1≤i≤M),
S53:计算所有Xi(1≤i≤M)的第P近邻域点的最小均方距离δ(K)以及输出空间相应的最小均方距离γ(K),,对(δ(K),γ(K))K(1≤K≤P)
S54:对(δ(K),γ(K))K(1≤K≤P),按公式γ=Aδ+R进行一次线性回归,所得一次线性函数的截距,即系统噪声方差值R;
S55:当新增加一个标准样品时,重复步骤S51至步骤S54,,得到每个样品对应噪声方差。
其中,所述PLS校正模型的主因子系数修正单元260根据所述精确噪声方差值、所述待测样品的样品数据以及当前时刻所述待测样品的预测值,通过KF算法修正当前时刻所述PLS校正模型的主因子系数过程中,
设PLS初始模型主因子数为l,主因子系数为:
w1,t1,v1,p1;w2,t2,v2,p2;……;wi,ti,vi,pi(i=1,2,3…,l);
其中:
vi=(tTy)/(tTt)=[vi1 vi2 ... vip]
将所述PLS初始模型中的所有系数值组成状态向量:
W=[w1Tt1Tv1p1T...wiTtiTvipiT]T(i=1,2,3…,l)
系统的状态方程和观测方程表示为:
其中,Yek为标样浓度,Wk为第k个标样修正时刻的主因子系数,Xk为第k个样品光谱矢量,Yrk为预测浓度。Vk为观测噪声,其统计特性为:
令
则观测方程为:Yek=HkWk+Dk+Vk。
其中,H表示状态变量Wk对测量变量Yk的增益。
Wk表示k时刻状态变量,在这里也即是用第k个标样修正时的PLS主因子系数。Dk没有特别的含义,只是个中间变量而已。这里是令所以则
通过上述实施方式可以看出,本发明提供的基于GT-KF-PLS近红外光谱自适应模型校正方法及系统,有效利用观测输入(样本近红外光谱)输出数据(样本化验值),提出样本有效噪声方差(Gamma test,GT)改进的KF-PLS模型校正方法;采用衰减记忆的GT对输入输出数据进行实时方差,得到准确的观测噪声方差值,再利用KF-PLS实现精确校正模型,能够保证近红外光谱自适应校正模型的稳定性,最终实现基于近红外光谱技术的在线分析。
如上参照附图以示例的方式描述了根据本发明提出的基于GT-KF-PLS近红外光谱自适应模型校正方法及系统。但是,本领域技术人员应当理解,对于上述本发明所提出的基于GT-KF-PLS近红外光谱自适应模型校正方法及系统,还可以在不脱离本发明内容的基础上做出各种改进。因此,本发明的保护范围应当由所附的权利要求书的内容确定。