本发明公开了一种光谱波数的选择方法，所述方法针对光谱的波数，多次随机抽样校正样本，建立偏最小二乘回归模型，计算每个波数的变量投影重要性系数，按降序排序，并得到与之相应排列顺序的波数集，对重排后的波数集进行逐步波数筛选，统计每次波数选择的结果，得到波数初选集；然后统计所有波数初选集内每个波数的偏最小二乘回归系数的绝对值并作相应的处理，并对处理后的偏最小二乘回归系数按降序排序，记录对应的波数排列次序，再采用反向剔除弱相关的波数的策略，从而获得最优特征波数集。本发明能充分挖掘波数中的有效信息，有效克服波数选择的主观性问题，最大限度地提取特征波数，消除弱相关因素的影响，显著提高了模型的稳健性与精度。
1.一种光谱波数的选择方法，其特征在于，所述光谱波数的选择方法针对光谱的波数，多次随机抽样校正样本，建立偏最小二乘回归模型，计算每个波数的变量投影重要性系数VIP，按降序排序，并得到与之相应排列顺序的波数集，对重排后的波数集进行逐步波数筛选，统计每次波数选择的结果，得到波数初选集；然后统计所有波数初选集内每个波数的偏最小二乘回归系数的绝对值并作相应的处理，并对处理后的偏最小二乘回归系数按降序排序，记录对应的波数排列次序，再采用反向剔除弱相关的波数的策略，从而获得最优特征波数集，具体步骤包括：步骤(1)，对于样本个数为m、波数个数为p的样本集，建立样本光谱矩阵X(m×p)，样本性质矩阵为Y(m×1)，样本原波数集V＝{x1,x2,…,xp}；设定最大统计次数T，设定建模过程中最大筛选次数KMAX，初始化统计次数t＝1；同时，将样本集划分为校正集及预测集；
步骤(2)，从步骤(1)所述的校正集内随机抽取S个样本作为校正样本集，对校正样本集建立偏最小二乘回归模型；
步骤(3)，基于步骤(2)中建立的偏最小二乘回归模型，通过变量投影重要性系数来判断每个波数对建立模型的重要程度，VIP值大的波数对模型的重要程度高，将VIP值按降序排列，并得到和VIP值排列顺序相对应的波数集Vt,1，然后，根据每个波数对模型的重要程度，逐步淘汰波数，并计算第k次淘汰波数时建模的交叉验证均方根误差RMSECVt,k，记录第k次淘汰波数时波数子集Vt,k，其中k为从1到KMAX的自然数；
步骤(4)，查找交叉验证均方根误差RMSECVt,k的最小值，将此次波数选择的波数子集Vt,k记为波数初选集selected_t，并记录波数子集selected_t对应的回归模型系数绝对值的矩阵w_t；
步骤(5)，重复步骤(2)到步骤(4)，直至达到最大统计次数T，获得波数子集selected_t，其中t＝1,2,…,T，同时得到每个波数子集对应的回归模型系数绝对值的矩阵集合{w_1,w_2,…,w_T}；
步骤(6)，统计出每个波数在步骤(5)中的集合{w_1,w_2,…,w_T}相应的回归模型系数值之和，并做归一化处理得到矩阵W1×p，将W1×p每个波数按降序排列，得到矩阵WS1×p；记与矩阵WS1×p相对应的原波数重排后的波数集为resort_wave＝{x′1,x′2,…,x′p}，其中x′1为原波数集内偏最小二乘回归系数最大的波数点，x′p为原波数集内偏最小二乘回归系数最小的波数点；
步骤(7)，针对步骤(6)所述的波数集采取反向剔除法，即从低位开始逐个剔除系数值小的波数；剔除个数设定阈值g，所述阈值g为自然数，所述阈值g大于等于2同时小于步骤(1)所述的波数个数p，在步骤(6)所述的波数集resort_wave＝{x′1,x′2,…,x′p}中，将波数子集{x′p-g+1,…,x′p-1,x′p}剔除，波数子集{x′1,…,x′p-g-1,x′p-g}保留并作为最优特征波数子集selected_wave。
2.根据权利要求1所述的一种光谱波数的选择方法，其特征在于，所述步骤(1)中，所述最大统计次数T的取值范围优选在100到1000之间。
3.根据权利要求1所述的一种光谱波数的选择方法，其特征在于，所述步骤(1)中，所述最大筛选次数KMAX的取值范围优选在20到200之间。
4.根据权利要求1所述的一种光谱波数的选择方法，其特征在于，所述步骤(7)中，所述阈值g的优选值设定方法为：基于所述步骤(6)划分的波数集，第1步首先从所述步骤(6)中所述的波数集resort_wave＝{x′1,x′2,…,x′p}中剔除偏最小二乘回归系数最小的波数点x′p，将剩余的波数子集{x′1,…,x′p-1}按所述步骤(1)中所述已划分好的校正集及预测集的样本建立PLS回归模型，计算波数集剔除数目为1时所建模型的综合评价指标CP1；第2步从波数子集{x′1,…,x′p-1}中剔除x′p-1，将剩余波数子集{x′1,…,x′p-2}按所述步骤(1)中所述已划分好的校正集及预测集的样本建立PLS回归模型，计算波数集剔除数目为2时所建模型的综合评价指标CP2；第k步从波数子集{x′1,…,x′p-k+1}中剔除x′p-k+1，将剩余的波数子集{x′1,…,x′p-k}按所述步骤(1)中所述已划分好的校正集及预测集的样本建立PLS回归模型，计算波数集剔除数目为k时模型的综合评价指标CPk；重复上述计算，直至第(p-2)步，计算波数集剔除数目为(p-2)时的模型综合评价指标CPp-2，此时仅剩偏最小二乘回归系数最大的两个波数x′1和x′2；所述CPi最小值所对应的波数剔除数目为阈值g的优选值，其中i＝
1,…,p-2；所述CPi的值为模型校正标准偏差RMSECV与模型预测标准偏差RMSEP的加和，其中i＝1,…,p-2。
一种光谱波数的选择方法
技术领域
[0001] 本发明涉及光谱分析领域，尤其是涉及一种光谱波数的选择方法。
背景技术
[0002] 光谱分析技术作为一种快速无损的定量分析方法，已成功应用于食品、农业、石油化工等领域。然而光谱检测中获得的谱图往往重叠严重，光谱信息冗余，特征吸收区域不明显。为了提高模型预测精度和简化模型，需要对波数进行优选，选择与待测样本信息最相关的特征波数用于模型的建立。
[0003] 目前，关于光谱分析中波长(波长是波数的倒数)选择方法的研究，主要有基于间隔的波长选择策略、无信息变量消除算法等方法。基于间隔的波长选择策略，波长区间大小的选择难以确定，且特征波长可能只是单独的点，而选择出的波长段中可能含有无关的波长，并不能最大程度地消除无关波长；无信息变量消除算法消除无关变量也是一种常用的波长选择方法，其采用对校正光谱矩阵人为产生一噪声矩阵，消除信息小于噪声变量的波长的方法进行特征波长提取，但这种方法也存在一定的主观性，并且噪声矩阵的选择影响着波长选择的结果，最终波长选择结果通常存在较多的冗余。上述现有方法的局限性，影响了模型的稳健性与精度。
发明内容
[0004] 本发明的目的在于，针对现有技术的不足，提供一种光谱波数的选择方法。
[0005] 为此，本发明采用如下解决方案：
[0006] 一种光谱波数的选择方法，所述光谱波数的选择方法针对光谱的波数，多次随机抽取校正样本，建立偏最小二乘回归模型，计算每个波数的变量投影重要性系数(Variable Importance in the Projection，VIP)，按降序排序，并得到与之相应排列顺序的波数集，对重排后的波数集进行逐步波数筛选，统计每次波数选择的结果，得到波数初选集；然后统计所有波数初选集内每个波数的偏最小二乘回归系数的绝对值并作相应的处理，对处理后的偏最小二乘回归系数按降序排序，记录对应的波数排列次序，再采用反向剔除弱相关的波数的策略，从而获得最优特征波数集。
[0007] 本发明的基于变量投影重要性系数与偏最小二乘回归系数协同的光谱特征波数选择方法，所述方法包括以下步骤：
[0008] 步骤(1)，对于样本个数为m、波数个数为p的样本集，建立样本光谱矩阵X(m×p)，样本性质矩阵为Y(m×1)，样本原波数集V＝{x1,x2,…,xp}；设定最大统计次数T，设定建模过程中最大筛选次数KMAX，初始化统计次数t＝1；同时，将样本集划分为校正集及预测集；
[0009] 步骤(2)，从步骤(1)所述的校正集内随机抽取S个样本作为校正样本集，对校正样本集建立偏最小二乘(Partial Least Squares,PLS)回归模型；
[0010] 步骤(3)，基于步骤(2)中建立的偏最小二乘回归模型，通过变量投影重要性系数(VIP)来判断每个波数对建立模型的重要程度，VIP值大的波数对模型的重要程度高，将VIP值按降序排列，并得到和VIP值排列顺序相对应的波数集Vt,1，然后，根据每个波数对模型的重要程度，逐步淘汰波数，并计算第k次淘汰波数时建模的交叉验证均方根误差RMSECVt,k(Root Mean Square Error of Cross Validation)，记录第k次淘汰波数时波数子集Vt,k，其中k为从1到KMAX的自然数；
[0011] 步骤(4)，查找交叉验证均方根误差RMSECVt,k(k＝1,2…,KMAX)的最小值，将此次波数选择的波数子集Vt,k记为波数初选集selected_t，并记录波数子集selected_t对应的回归模型系数矩阵绝对值w_t；
[0012] 步骤(5)，重复步骤(2)到步骤(4)，直至达到最大统计次数T，获得波数子集selected_t(t＝1,2,…,T)，同时得到每个波数子集对应的回归模型系数矩阵绝对值集合{w_1,w_2,…,w_T}；
[0013] 步骤(6)，统计出每个波数在步骤(5)中的集合{w_1,w_2,…,w_T}相应的回归模型系数值之和，并做归一化处理得到矩阵W1×p，并将W1×p每个波数按大小从高至低排列，得到矩阵WS1×p；记与矩阵WS1×p相对应的原波数重排后的波数集为resort_wave＝{x′1,x′2,…,x′p}，其中x′1为原波数集内偏最小二乘回归系数最大的波数点，x′p为原波数集内偏最小二乘回归系数最小的波数点；
[0014] 步骤(7)，针对步骤(6)所述的波数集采取反向剔除法，即从低位开始逐个剔除系数值小的波数；剔除个数设定阈值g，所述阈值g为自然数，所述阈值g大于等于2同时小于步骤(1)所述的波数个数p，在步骤(6)所述的波数集resort_wave＝{x′1,x′2,…,x′p}中，将波数子集{x′p-g+1,…,x′p-1,x′p}剔除，将波数子集{x′1,…,x′p-g-1,x′p-g}保留并作为最优特征波数子集selected_wave。
[0015] 优选地，所述步骤(1)中，所述最大统计次数T的取值范围优选在100到1000之间。
[0016] 优选地，所述最大筛选次数KMAX的取值范围为1到波数个数p，优选在20到200之间。
[0017] 优选地，所述步骤(7)中，所述阈值g的优选值设定方法为：基于所述步骤(6)划分的波数集，第1步首先从所述步骤(6)中所述的波数集resort_wave＝{x′1,x′2,…,x′p}中剔除偏最小二乘回归系数最小的波数点x′p，将剩余的波数子集{x′1,…,x′p-1}按所述步骤(1)中已划分好的校正集及预测集的样本建立PLS回归模型，计算波数集剔除数目为1时模型的综合评价指标CP1；第2步从波数子集{x′1,…,x′p-1}中剔除x′p-1，将剩余波数子集{x′1,…,x′p-2}按所述步骤(1)中已划分好的校正集及预测集的样本建立PLS回归模型，计算波数集剔除数目为2时模型的综合评价指标CP2；第k步从波数子集{x′1,…,x′p-k+1}中剔除x′p-k+1，将剩余的波数子集{x′1,…,x′p-k}按所述步骤(1)中已划分好的校正集及预测集的样本建立PLS回归模型，计算波数集剔除数目为k时模型的综合评价指标CPk；重复上述计算，直至第(p-2)步，计算波数集剔除数目为(p-2)时的综合评价指标CPp-2，此时仅剩偏最小二乘回归系数最大的两个波数x′1和x′2；所述CPi(i＝1,…,p-2)最小值所对应的波数剔除数目为阈值g的优选值；所述CPi(i＝1,…,p-2)的值，优选模型校正标准偏差RMSECV与模型预测标准偏差RMSEP的加和。
[0018] 本发明所提供的光谱波数的选择方法具有以下优点：
[0019] (1)采用随机采样的方式来确定校正样本，消除人为划分校正样本建模的主观因素；
[0020] (2)结合变量投影重要性系数与偏最小二乘回归系数共同进行波数筛选，最大化地甄别出有益波数；
[0021] (3)统计多次波数选择的结果，克服由于由所选择的校正样本对于波数重要性评价的影响，最大限度地提取特征波数；
[0022] (4)经过优选波数变量，能使用较少波数建立简化模型，且所建模型稳健性好、精度高。
附图说明
[0023] 图1为本发明基于变量投影重要性系数与偏最小二乘回归系数协同的光谱特征波数选择方法的流程图；
[0024] 图2为生物柴油调和油拉曼光谱的原始光谱图；
[0025] 图3为全部波数及其回归模型系数指标的统计图，圆圈标出的点是最终筛选出的特征波数；
[0026] 图4为波数子集数目与模型的综合评价指标CP的关系图；
[0027] 图5为筛选出的特征波数与拉曼光谱谱图的关系；
[0028] 图6为生物柴油调和油浓度的预测值与真实值的比较图。
具体实施方式
[0029] 本发明的具体实施方式，将结合附图和拉曼光谱分析的具体实施例进行详述。
[0030] 图1为本发明提出的基于变量投影重要性系数与偏最小二乘回归系数协同的光谱特征波数选择方法的流程图。
[0031] 本具体实施例，采用生物柴油调和油拉曼光谱数据样本来验证本发明的方法。该数据集包含62个拉曼光谱测量样本，生物柴油的质量含量范围为0％～100％(w/w)，考察拉曼光谱与调和油中生物柴油浓度之间的关系。拉曼光谱波数区间进行线性插值处理后得到共2033个波数，样本的原始光谱图如图2所示。
[0032] 针对生物柴油调和油拉曼光谱数据样本，实施本发明提出的基于变量投影重要性系数与偏最小二乘回归系数协同的光谱特征波数选择方法的步骤如下：
[0033] 步骤(1)，对于样本个数为62、波数个数为p＝2033、原波数集为V＝{x1,···,xp-1,xp}＝{1,···,2032,2033}的生物柴油调和油拉曼光谱样本集，建立样本光谱矩阵X(m×p)，表征调和油中生物柴油浓度的性质矩阵为Y(m×1)；由于建模样本的随机性，所述最大统计次数T其取值范围优选在100到1000之间，T值越大，波数选择所需的时间越长，T值太小又难以消除建模样本的随机性对结果的影响，本具体实施例中设定T＝500；所述最大筛选次数KMAX决定着每次波数初选过程的效率，其取值范围为1到波数个数p，本具体实施例中KMAX＝131；初始化统计次数t＝1。将原始62个生物柴油调和油拉曼光谱样本划分为包含44个校正集样本和18个预测集样本。
[0034] 步骤(2)，从步骤(1)所述的校正集内随机抽取S＝0.8×44(取整)，即35个样本作为校正样本集cal。对校正样本集按公式(1)建立潜变量个数为A的偏最小二乘(Partial Least Squares,PLS)回归模型；该模型的数学表达式如下：
[0035] Xcal＝TPT+E (1)
[0036] Ycal＝UQT+F (2)
[0037] T＝XW(PT W)-1 (3)
[0038] 其中校正光谱矩阵Xcal的得分矩阵是T，载荷矩阵是P，权重矩阵是W，残差光谱是E；校正样本浓度性质矩阵Ycal的得分矩阵是U，载荷矩阵是Q，残差光谱是F。
[0039] 步骤(3)，基于步骤(2)中建立的PLS回归模型，按照公式(2)计算每个波数对模型的重要程度指标：变量投影重要性系数VIP。第j(j＝1,···,p)个波数的VIP值的计算公式如下：
[0040]
[0041] 其中wa，ta和qa分别是矩阵W，T和Q的第a列。wja是第j个波数在第a个潜变量中的权重值。变量投影重要性系数VIPj越大，表示波数j对于建立回归模型的重要程度越高。
[0042] 然后，将变量投影重要性系数VIP按降序排序，采用半折法逐步淘汰波数：
[0043] 第1次淘汰波数时建模采用的波数个数为p，计算此时建模的交叉验证均方根误差RMSECVt,1，记录此时波数集Vt,1；第2次淘汰波数时被淘汰的波数是波数集Vt,1的后一半，采用剩余的波数进行建模，计算第2次淘汰波数时模型的交叉验证均方根误差RMSECVt,2，记录此时波数子集Vt,2；直到第k次淘汰波数时，判断上一次波数子集Vt,k-1中波数个数是否大于2且小于200，若大于2且小于200时淘汰波数集Vt,k-1中最后一个波数，否则淘汰波数子集Vt,k-1的后一半，剩余的波数进行建模，计算第k次淘汰波数时建模的交叉验证均方根误差RMSECVt,k，记录此时波数子集Vt,k；重复上述计算，直到第KMAX次淘汰波数时建模采用的波数只剩下两个重要程度最高的波数，计算第KMAX次淘汰波数时建模的交叉验证均方根误差RMSECVt,KMAX，记录此时波数子集Vt,KMAX。
[0044] 步骤(4)，查找交叉验证均方根误差RMSECVt,k(k＝1,2,···,KMAX)中的最小值，将此次波数选择的波数子集Vt,k记为波数初选集selected_t；记录波数初选集对应的回归系数矩阵w_t。矩阵w_t即公式(5)中矩阵B的绝对值：
[0045] YcalsX(cal,selected_t)B+e (5)
[0046] 其中X(cal,selected_t)是校正集光谱矩阵Xcal中对应的波数初选集selected_t列，校正模型的回归系数是B＝[b1,b2,···,bp]T，e是模型偏差。
[0047] 步骤(5)，重复步骤(2)到步骤(4)，直至达到最大统计次数T，获得波数回归系数绝对值矩阵集{w_1,w_2,···,w_T}，统计出达到最大统计次数T后每个波数对应的回归系数绝对值之和矩阵 本样本集的回归系数绝对值和矩阵为W1×2033＝[w1,w2,···,w2033]。
[0048] 步骤(6)，对步骤(5)所述的回归模型系数绝对值和矩阵W1×2033做归一化处理：得矩阵W′1×2033。回归模型系数均值与波数关系如图3所示；
然后将原波数按回归模型系数均值大小降序排序，得到波数集resort_wave＝
{x′1,···，x′p-1,x′p}，x′1为原波数集内系数最大的波数点，x′p为原波数集内系数最小的波数点；
[0049] 步骤(7)，为步骤(6)所述的波数子集个数设定阈值g，在步骤(6)所述波数集resort_wave＝{x′1,···，x′p-1,x′p}中，将波数子集{x′p-g+1,···，x′p-1,x′p}剔除，将波数子集{x′1,···，x′p-g-1,x′p-g}保留并作为最优特征波数子集selected_wave。本样本集的波数集为resort_wave＝{x′1,···，x′2032,x′2033}；
[0050] 步骤(7)中阈值g的优选值设定方法具体是：第1步首先从所述步骤(6)中所述的波数集{x′1,···，x′2032,x′2033}中剔除偏最小二乘回归系数最小的波数子集x′2033，将剩余的波数子集{x′1,···，x′2032}按所述步骤(1)中所述已划分好的校正集及预测集的样本建立PLS回归模型，计算波数集剔除数目为1时所建立模型的综合评价指标CP1；第2步从波数子集{x′1,···，x′2032}中剔除x′2032，将剩余的波数子集{x′1,···，x′2031}按所述步骤(1)中所述已划分好的校正集及预测集的样本建立PLS回归模型，计算波数子集剔除数目为2时所建立模型的综合评价指标CP2；第k步从波数子集{x′1,···，x′2034-k}中剔除所述步骤(6)中所述的x′2034-k，以剩余的波数子集{x′1,···，x′2033-k}对所述步骤(1)中所述样本集的已划分好校正集及预测集的样本建立PLS回归模型，计算波数子集剔除数目为k时所建立模型的综合评价指标CPk；重复上述计算，直至第2031步，计算波数子集剔除数目为2031时的综合评价指标CP2031，此时仅剩偏最小二乘回归系数最大的最后两个波数。
[0051] 图4为波数子集数目与模型的综合评价指标CP的关系图；查找上述指标CP所有2031个值中的CP最小值，模型综合指标CP为最小值时(图4中圆圈内的点)所对应的波数子集内波数数目为364；故阈值g的优选值为1669，将波数子集{x′365,···，x′2032,x′2033}删除，得到最优特征波数子集selected_wave＝{x′1,···，x′363,x′364}，其中x′1＝1964，x′p＝209；共364个特征波数，仅占原始波数的17.9％。
[0052] 图5是最优特征波数在生物柴油调和油样本浓度为2％的拉曼光谱中的分布情况。
[0053] 结果验证：采用上述方法所得的最优特征波数子集selected_wave建立校正模型，图6为采用上述选择的最优特征波数集selected_wave建立PLS回归模型得到的真实值与预测值关系。表1比较了本发明方法与全波数方法以及相关系数法(CC)、无信息变量消除算法(UVE)、基于间隔的波长选择策略(iPLS)等常用的波数选择方法，结果表明本发明方法采用了较少的波数个数，并获得了预测均方根误差RMSEP(Root Mean Square Error of Prediction)为最小的理想结果。
[0054] 表1的结果充分表明：本发明方法能高效地提取最优特征波数，可靠地精简模型，显著地提高模型的稳健性与精度。
[0055] 表1
[0056]
[0057] 本发明上述具体实施方式中所涉及的光谱为拉曼光谱，仅为优选实施例，具体实施时，也可以用于紫外光谱、可见光谱、近红外光谱、中红外光谱、远红外光谱、荧光光谱等光谱的特征波数选择。
[0058] 本发明涉及光谱特征波数选择方法，由于根据光学的基本原理，光的波数与光的波长、光的频率这三种参数之间可以直接相互换算，故本发明的方法也可以称为光谱特征波长选择方法或光谱特征频率选择方法。
[0059] 上述具体实施方式用来解释说明本发明，仅为本发明的优选实施例而已，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改、等同替换、改进等，都落入本发明的保护范围。
