样条变换集成罚函数偏最小二乘方法用于光谱数据重构和定量分析

【摘要】 针对高维小样本光谱数据所显现的函数型数据（Functional data）特性、与性质参数的非线性关系及变量间存有的严重共线性，采用了样条变换集成罚函数偏最小二乘回归新技术。它首先以三次B基样条变换实现非线性光谱数据的线性化重构，随后将重构的新光谱矩阵交由罚函数偏最小二乘法（Penalized PLS）构建其与性质参变量间的校正模型，其中罚函数中的光滑因子由交叉验证优化确定以调控模型的拟合精度。最后，通过小麦样品水分含量的近红外光谱定量分析，结果显示该技术光谱数据重构稳健，去噪明显，并有效解决高维小样本的过拟合和变量间的共线性，而预测集的均方根误差（RMSEP）为0.1808%，方法的非线性校正模型预测能力得到了明显提高。

【关键词】 样条函数， 偏最小二乘， 粗糙惩罚， 近红外光谱， 定量分析， 小麦

1 引言 现代光谱以其分析速度快、重现性好、成本低、不消耗样品、易于实现在线分析等特点而得到广泛应用。而光谱化学计量学是近代红外光谱分析技术的重要组成部分，它通过多变量校正技术来进行数据（样本光谱和其性质参数）处理，以获得准确的分析结果[1，2]。考虑到近红外光谱数据通常呈多变量、强相关性，并与样品性质参变量间的非线性关系，适宜选用非线性偏最小二乘法（Nonlinear PLS，NLPLS）。目前，NLPLS实现方式有3种：一是基于样本矩阵的非线性变换，即在建模自变量中引入某些原始变量的非线性项，如二次项、交叉项等[3]；二是将建模变量投影到低维的曲线或曲面上得到非线性特征向量，再建立输入输出特征向量间的非线性关系[4]，但该方法计算复杂，建模受初值影响大；三是保留PLS的线性外部模型，而内部模型采用多项式、样条函数、模糊规则、神经网络、支持向量机等非线性形式[5～9]，该方法缺乏对建模物理变量的直观解释能力。 鉴于光谱变量与性质参变量间的具体非线性依存关系不明确，及样本个体光谱数据显现为波长变量的函数型数据（Functional data）特性[10]，本研究采用样条（Spline）变换集成罚函数偏最小二乘（Penalized PLS）回归新技术，记为SplinePPLS方法。首先利用样条基函数将光谱自变量与性质因变量之间的未知非线性关系按照各维自变量与因变量的拟线性关系相加展开[11]。由于样条函数分段拟合、可按需要裁剪以适应任意曲线连续变化的特点，使光谱的重构函数曲线适应光谱数据局部敏感特性的同时保持了函数的光滑性和连续性，从而可削减原始数据中的噪声。随后，考虑到光谱矩阵经样条变换后变量维数显著增加，将重构的新光谱矩阵交由罚函数偏最小二乘法构建其与性质参变量间的定量线性校正模型，其中基于转换权向量二阶导数的罚函数用以调控模型的拟合精度。为考察SplinePPLS方法的有效性及性能，对小麦近红外光谱数据进行了研究。

2 SplinePPLS方法的构建

2.1 B基样条曲线 设变量λ与x满足如下随机模型：x=s(λ)＋ε，Ｅ（ε）＝０，Ｖar(ε)＝σ２（１）若λ在区间［a，b］上的一个Ｍ段划分π∶a=ξ０＜ξ１＜…＜ξＭ＝b， 则式(1)中s(λ)的三次B基样条逼近曲线方程[11]为：s(λ)=∑M+2l=0clΩ3λ－ξl－1h， a≤λ≤b（2）式中Ω3λ－ξl－1h=13！h3∑4k=0(－１)k4

k(λ－ξl+k－３)３＋，是以ξl+k－３（k=0，1，2，3，4）为内控节点、步长为h的三次B基样条函数，它与x呈线性关系。对于分点ξl－１及其内控节点ξl+k－３位于划分的两侧，本研究取 ξ－３＝ξ－２＝ξ－１＝ξ０和ξＭ＝ξＭ＋１＝ξＭ＋２＝ξＭ＋３［１２］。 取观测位置λ１， λ2， λ…， λp相应数据点x1， x2， …， xp与它们在样条曲线上插值映射点的距离平方和最小为目标函数，优化求取各基函数的线性加权系数cl(l=0， 1， …， M＋２)，即共有Ｋ＝Ｍ＋３个基函数，从而完成该序列数据形如式（2）的B基样条逼近曲线构造。

2.2 罚函数偏最小二乘方法（Penalized PLS） 数据点xj(j=1，2，…，p)在B基样条映射的线性空间中的插值映射点zj的分坐标定义为：zj，0=Ω３（λj－ξ－１h）， zj，l=Ω3(λj－ξ0h， …，，zj，M+2=Ω3(λj－ξＭ＋１h)（3）现将样本个体自变量x在p维变量空间中的取值，即x=(x1，x2，…，xp)T，计算其所有插值映射点zj各分坐标并加以组合，即可得到x的B基插值映射点矢量 z，即为z=(z1，0， z1，1， …， z1，Ｍ＋２， z２，0， z２，1，…， z２，Ｍ＋２， …， zp，0， zp，1， …，zp，M+2)T（4）

由于每一维变量xj有K个映射分坐标，故矢量z的空间维数将为p×Ｋ。 对于样本容量n的自变量矩阵Xn×p，欲构建其与性质矢量y间的非线性校正模型，则先实施X的每一样本个体xi(i=1，2，…，n)的B基样条变换（各样本个体选取基函数个数应相同，即K1=K2=…＝Ｋn=K），得到映射样本矩阵Zn×（p×K），而其与y已演变为拟线性关系。再实施Z与y的线性PLS算法。 鉴于Z较X的变量维数显著增加，变量间的相关性更为严重，本研究采用罚函数偏最小二乘（Penalized PLS，PPLS）方法[10]构建 Z与y间的校正模型。PPLS方法的目标函数为

arg maxwwTZTyyTZwwTw+P(w)(5)

上式分母中基于转换权向量w的罚函数P(w)=wT(ΔφK2)w， 其中K2=（DK－１ＤＫ）ＴＤＫ－１ＤＫ，而DK为(K－１)×Ｋ维的w一阶差分算子矩阵。另外， Δφ＝diag(φ１，φ２，…，φp)为各初始自变量的光滑因子对角阵。Pw实为二阶导数罚函数，意在增强w平滑其特征向量t=Zw的能力，以提高模型的稳健性。 由上述PPLS算法思想可见，它是以放弃拟合精度为代价寻求预测性能更优的一种改进PLS方法。同时，PPLS对病态数据的耐受性远强于普通PLS方法。现将基于初始样本阵X，y及各自变量光滑因子对角阵Δφ的B基样条变换的PPLS算法（SplinePPLS）步骤[10]归结如下：（1）实施Ｘ的n个B基样条逼近曲线的优化构造，并得其映射矩阵Z；（2）令h=1， Zh=Z，并计算P=ΔφK2及M=(Ip+P)－１，其中Ｉp为p维单位阵；（3）计算转换权向量wh=MZThy， 并规一化wh=wh/‖wh‖；（4）计算特征向量th=Zhwh，并规一化 th=th/‖th‖；（5）记Th=［t1，t2，…，th］，计算其正交投影矩阵Qh=Th(TThTh)+TTh， 式中“+”为矩阵广义逆；（6）计算剩余矩阵Zh+1=Zh－ＱhZh；（7）令h=h+1， 重复步骤（3）～（7），直至由交叉验证(Cross validation )法[13] 确定所需提取的最优成分数h后，将得到转换权矩阵W=[w1，w2，…，wh]， 进而可计算Ｚ与y间的线性回归系数β＝（Ｗ（ＷＴＺＴＺＷ）－１ＷＴＺＴy）zy。

3 SplinePPLS方法为小麦近红外光谱重构及定量分析

3.1 样本数据说明 小麦样品数据取自文献[14]，自变量取其近红外光谱在波长1100～2500 nm、扫描分辨率为2 nm的若干波长处的吸光度值log(1/R)， 即维数p=701，其中R为样本的反射率，样本容量n=100，谱图如图1所示。性质参变量为小麦水分质量百分含量，数值范围在12.45%～17.36%。从原始数据集中随机划出80个构成训练集用于光谱B基样条变换的优化确定及后继校正模型建立，其余20个组成独立测试集，用于检验B基样条的光谱插值重构能力及模型的预测性能。

Fig.1 NIR diffuse reflectance spectra of wheat samples

3.2 实验方式与性能评价指标 为检验SplinePPLS方法的性能，先将训练样本Xntrain×p以“变量‘留一’交叉验证”[13]选定K个样条变换的B基函数，其优化评定指标为式(6)中的RMSECVspline；再以“样本个体‘留一’交叉验证”选定各变量的光滑因子φ1，φ2，…，φp及校正模型所需PLS最优成分数h，它们的优化评定指标则为式(6)中的RMSECVppls。而光谱B基样条变换的插值重构性能及校正模型的预测性能，则交由测试样本Xntest×p计算，它们的评价指标分别为式(7)中的RMSEPspline和RMSEPppls。RMSECVspline=∑pj=1∑ntraini=1(xij－ij)2/(ntrain×p)， RMSECVppls=∑ntraini=1(yi－i)2/ntrain（6）

RMSEＰspline=∑pj=1∑ntesti=1(xij－ij)2/(ntest×p)， RMSEＰppls=∑ntesti=1(yi－i)2/ntest（7）式中xij和ij分别为第i样本个体、第j波长变量下吸光度的实验值和B基样条曲线的插值；yi和i则分别为第i样本个体性质参变量的实验测试值和模型预报值。

3.3 结果与分析 3.3.1 光谱数据的B基样条变换 现基于光谱阵Xntest×p，构造它们的３次B基样条逼近曲线，据此完成测试光谱阵Xntest×p的三次B基样条变换。即，先将光谱波长1100～2500 nm作预设的M段划分，再实施Xntrain×p的“变量‘留一’交叉验证”实验，即依次留用一个测量位置(变量)λj的观测数据矢量x·j=[x1j，x2j，…，xntrainj]T作内部验证，而剩余p－１列观测数据用于样条逼近曲线的最小二乘拟合，可得ntrain个样条曲线各K=M+3个基函数的线性加权系数ci，j(l=0，1，…，M+2，i=1，2，…，ntrain)， 并据此实现x·j的插值估计·j。最后，将p轮循环得到的所有·j(j=0，1，…，p)代入式(6)，即可计算该M取值下的RMSECVspline。改变M取值，并依据对插值精度RMSECVspline的要求，即可选定所需的基函数个数K。 图2 B基样条曲线的基函数个数优化确定及第1训练样本个体的光谱重构结果（略）

Fig.2 Selection of basic functions number and the rebuilding spectrum curve for the first sample data

a. Rootmean squared error at different numbers of basic functions; b. 28 Basic functions and its rebuilding spectrum curve.

图2a显示了RMSECVspline与K的相关关系，随着K的增多，RMSECVspline总体呈下降趋势，在K=28时，RMSECVspline已很小，而其后的RMSECVspline值下降有限。考虑到参数K取值越大，经B基样条变换后的拟线性变量个数将越多，且变量间将出现更多复共线性， B基样条逼近曲线易出现对训练数据的过拟合，而对包含于数据中的噪音削减不够；若K取值太小，B基样条逼近曲线对数据的插值能力将下降。因此，本研究选定K=28。图2b显示了这28个基函数及由它们所重构的第1训练样本个体光谱数据的三次B基样条逼近曲线。由图2b可见，该样条曲线实现了对光谱数据较高精度的插值拟合。另外，将这28个基函数用于Xntest×p的插值重构，其精度指标RMSEＰspline，表明B基样条具有极强的插值重构能力。

3.3.2 SplinePPLS方法中参数的选择 对于模型而言，预报性能最为重要。影响SplinePPLS模型预报精度的主要因素有：B基样条变换基函数个数K、PPLS模型中各自变量光滑因子φ1，φ2，…，φp及PLS最优成分数h等。其中，K的优选过程见3.3.1节，并将其取值为K=28。而表1则为PPLS方法施于样本阵｛Ｘ，y｝优化选择φ1， φ2， …， φp及h的过程结果。其中，为减少搜索空间维数，将各光滑因子简单取值相同φ=φ1=φ2=…，＝φp。

表1 PPLS方法中参数的优化选择（略）

Table 1 Parameter selection of the penalized PLS method

从表1可见，模型性能指标ＲＭＳＥＣＶppls在参数φ和h的二维搜索格子点φ＝１０００００， h=5位置达最小，由此选定它们为PPLS模型相应参数的最优取值。同时还发现，该位置ＲＭＳＥＣＶppls指标小于φ＝０， h=5的PLS方法的结果。由2.2节PPLS方法的目标函数式(5)可知，φ取值的大小将影响w平滑其特征向量t=Zw的程度，并进一步作用于PPLS方法的模型系数β＝（Ｗ（ＷＴＸＴＸＷ）－１ＷＴＴＴy）x，y。图3即为PLS和PPLS方法施于样本阵｛Ｘ，y｝的β结果比较。其中PPLS模型系数β曲线受到变量光滑因子φ的粗造惩罚而得到了平滑。这样，它既可最大限度保证分析信号不失真，又能进一步削除噪音。

图3 PLS和PPLS模型回归系数比较（略）

Fig.3 Comparison of model coefficients for the PLS and Penalized PLS method

3.3.3 模型精度比较和分析 表2为多元线性回归MLR（Multiple linear regression）， PLS， PPLS， SplinePLS及SplinePPLS的建模结果，各方法的参数最优值由交叉验证法[13]确定。另外，SplinePLS和SplinePPLS方法中Ｋ＝２８。先分析校正方法对模型精度的影响，SplinePLS和SplinePPLS属NLPLS校正方法，它们的RMSECV和RMSEP分别小于同属线性校正的PLS和PPLS方法，但它们提取的PLS最优成分数h均多于后两种线性方法。由此说明，B基样条变换在一定程度上实现了光谱数据与样品性质变量间的内在非线性关系，但需通过增加PLS成分将包含在拟线性化变量项中的非线性信息带入校正模型。而MLR方法未进行光谱数据的噪音削减，以及消除变量间的复共线性，RMSECV和MSSEP值均显著高于其它4种方法。PPLS和SplinePPLS方法的RMSECV和RMSEP值均分别小于未进行粗糙惩罚的PLS和SplinePLS方法的。由此说明，在用于PLS特征向量提取的目标函数中集成转换权向量二阶导数的罚函数，可平滑特征向量和进一步削减噪音，从而使模型的预测能力和稳定性得以提高。

表2 5种不同校正模型的性能比较（略）

Table 2 Comparison of model performance for five methods