基于小波变换动态时间规整的近红外光谱模型传递方法

王其滨，杨辉华,，潘细朋，李灵巧,*

(1.桂林电子科技大学 电子工程与自动化学院，广西 桂林 541004；2.北京邮电大学 自动化学院，北京 100876)

近红外光谱分析技术因具有实时、高效、低成本以及不破坏样本等优点，被广泛应用于食品、药品、化工等领域[1-4]。该技术通过建立分析模型，将训练样本的近红外光谱与物理化学性质相联系，可计算出未知样本对应的性质参数[5]。

模型传递方法是近红外光谱分析技术中一个常见并具有重要意义的技术问题[6]。通常在实际的应用中，因为测量环境的改变或光谱仪器在制造过程的机械差异等原因，源机(Master)上建立的定量分析模型应用于目标机(Slave)测量的光谱数据时会产生较大的偏差[7]，不能得到正确的预测结果。如果在目标机上再次测量多组样本数据建立定量模型则费时费力[8]，此问题严重阻碍了近红外光谱仪在实际生产中的应用，而模型传递技术很好地解决了这一问题。模型传递是建立源机与目标机上测量光谱之间的数学模型[9]，该模型通过对目标机上的近红外光谱进行转换，可将其用源机上建立的模型进行参数分析。模型传递方法主要分为有标样和无标样模型传递。目前，有标样模型传递算法应用更为广泛，如直接校正算法(DS)[10]、分段直接校正算法(PDS)[11]以及Shenk's算法[12]等。DS算法通过选取少量样本建立标样集，将目标机上测量的全光谱数据转换成源机光谱数据来构建校正模型[13]。PDS算法在目标机光谱中增加窗口限制[14]，选取窗口内光谱波长点转换成源机光谱，代替DS算法中全光谱进行校正。Shenk's算法利用某波长点周围的光谱数据来计算该波长的校正参数。

光谱仪器采集的光谱数据集普遍存在基线漂移、随机噪声、漫反射等现象，如果直接对光谱进行建模分析，会出现模型传递效果差、不稳定等问题，故本文尝试对光谱数据进行小波变换预处理。小波变换作为一种源于傅里叶变换的信号处理工具[15]，为一种时间和频率的局域变换，因而能有效的从信号中提取信息[16]。小波可以沿时间轴前后平移，也可以按比例伸展和压缩以获取低频和高频小波，从而对函数或信号进行多尺度细化分析[17]。通过一维小波连续变换对近红外光谱进行预处理，能够实现对光谱数据平滑、降噪以及消除基线漂移的处理。动态时间规整算法是一个优化问题[18]，算法通过计算两个序列点之间最佳的对应关系，获得代价矩阵最小时的规整路径。在近红外光谱领域，由于设备老化、机械误差以及温度湿度等条件的影响，源机与目标机测得的光谱之间不但会产生基线漂移，而且波长点也会产生偏移[19]。通过寻找源机光谱与目标机光谱的最佳对应关系，可以建立更精确的校正模型。

1 WDTW算法

WDTW模型传递方法，首先将光谱进行小波变换预处理，然后计算源机光谱与目标机光谱上各个波长点的相关距离，通过DTW算法找到源机光谱与目标机光谱各个波长点的关联关系，使得整体的相关距离最小，最后再根据此关联关系构建回归模型。通过在近红外药品光谱数据集和汽油光谱数据集之间进行传递模型实验，实验结果相较于传统方法光谱差异和预测标准偏差均有所降低。

1.1 波长点相关距离计算

设两个序列分别表示为X=(x1,x2,x3,…,xN)和(y1,y2,y3,…,yM)。DTW算法可以将序列X和Y中的点进行最佳匹配，使两序列点之间的距离和最小。序列和匹配关系示意图如图1所示。

图1 序列的匹配关系Fig.1 The matching relation of sequence

(1)

式中，Cov表示光谱序列之间的协方差，Var表示光谱序列的方差。源机光谱和目标机光谱序列之间的相关距离可根据相关系数求出：

(2)

DTW算法通过构造一个代价矩阵C∈RN×M表示序列X与序列Y中每个元素之间的距离。

C∈RN×M∶ci,j=Di,j

(3)

通常在语音识别领域，DTW算法是通过ci,j=‖xi-yi‖计算元素之间距离，而光谱数据更关注于两个光谱序列波长点的相关性，因此用相关距离Di,j。

图2 最佳规整路径Fig.2 The optimal warping path

1.2 规整路径的计算

根据代价矩阵，DTW算法能够找到元素之间相关距离最小的路径，该路径代表了元素之间的最佳对应关系。图2表示了两个序列X和Y的最佳规整路径示意图。

规整路径的计算必须满足3个条件：

① 边界限制：规整路径必须从(1,1)点出发，在(N,M)点结束。

② 顺序限制：规整路径需要按照顺序匹配，不能交叉或跳过某点进行匹配。

③ 单调限制：路径匹配过程中，只能向右、向上或者右上方向移动，以避免无意义的循环。

如顺序限制和单调限制条件：设当前点为(i,j)，下一个将要匹配的点只能是(i+1,j)、(i,j+1)或者(i+1,j+1)。规整路径中相关联点代价的和即为DTW算法的损失函数：

(4)

DTW算法要求解的就是使关联代价最小化：

DTW(X,Y)=min{cp(X,Y),p∈PN×M}

(5)

通过应用动态规划来求解累加距离最小的路径，其计算公式如下：

D(i,j)=min{D(i-1,j-1),D(i-1,j),D(i,j-1)}+c(xi,yj),i∈[1,N],j∈[1,M]

(6)

式(6)中,c(xi,yj)为当前所在点的关联代价，累加距离最小即为其与前一匹配点相关距离的最小值之和。累加距离最小的路径即为规整路径。

1.3 校正模型的建立

(7)

(8)

通过PLS方法求得方程系数，将源机光谱点分别做如上运算，然后将求得的全部系数保存于矩阵F，F即为模型传递所求的传递矩阵。通过如下公式即求得传递后的光谱。

(9)

2 实验部分

2.1 实验数据

分别采用药品和汽油2个数据集验证了WDTW模型传递方法的有效性。数据集1是由国际漫反射会议(IDRC)公布的药品近红外光谱数据集(http://www.eigenvector.com/data/tablets/index.html)，包括了药品的活性、重量和硬度3个性质参数。数据集2是由中国石油化工股份有限公司石油化工科学研究院提供的汽油近红外光谱数据集，其包括了汽油C7、 C8、 C9、C10 4个成分的性质参数。2个数据集样本数、样本维度、波长范围及波长间隔如表1所示。

表1 实验数据集的详细介绍Table 1 Detailed introduction of experimental data

2.2 数据预处理

光谱仪器采集的光谱数据集普遍存在基线漂移、随机噪声、漫反射等现象，故需要对光谱数据进行预处理。预处理过程首先采用数据平滑消除光谱数据中可能存在的随机噪声(平滑窗口为17，多项式阶数为3)，然后再使用小波变换(小波尺度Scale为4，母小波波长为32*Scale)进行数据处理。

2.3 建立定量模型及评价指标

采用K-S算法选取药品和汽油数据的标样集，通过计算近红外光谱之间的欧氏距离，选取最具代表性的光谱样本子集即为标样集。首先将实验数据成分进行排序，然后采用隔三取一的方式从汽油数据集中选取60个测试集，28个校正集；同样从药品数据集中选取460个测试集，155个校正集。最后采用10折交叉验证，据预测标准偏差(SEP)得到最优的PLS模型主因子数，用PLS建立源机定量模型。采用光谱平均差异(ARMS)、预测标准偏差(SEP)和光谱校正率(Prcorrected)作为WDTW模型传递方法的评价指标。其中，ARMS越小，说明两个不同的光谱仪器测得的光谱差异越小;SEP越小,表明所建模型的预测能力越强；Prcorrected值越小，代表模型的传递效果越好。ARMS和SEP计算公式如下：

(10)

(11)

(12)

3 结果与讨论

3.1 数据预处理结果

实验预处理过程先采用数据平滑，再进行小波变换处理。前者主要用于消除光谱数据中可能存在的随机噪声，后者用于消除数据可能存在的基线漂移。实验结果显示，通过将药品和汽油样本预处理前后的目标机仪器光谱图进行对比，可以发现预处理过程能够有效消除样本数据间普遍存在的基线漂移、随机噪声和漫反射等现象。图3为预处理前后药品样本目标机仪器光谱图。

图4 药品样本源机和目标机光谱间波长点匹配关系Fig.4 The match relation of the master and slave instrumental spectrograms wavenumber of the drug sample

3.2 WDTW算法波长点匹配

通过WDTW算法，可以得到源机光谱和目标机光谱波长点之间最优的对应关系。以药品近红外光谱数据为例，由于光谱波长点较多，只选取1 000～1 160 nm波长点显示，其中上方黑色线为目标机光谱(Slave)，下方蓝色线为源机光谱(Master)，中间红色线即连接了源机光谱和目标机光谱之间最优对应关系(见图4)。从图中红色连接线可看到，在波长点10～15 nm、20～25 nm和120～130 nm处，波长点之间发生明显偏移，而WDTW算法很好地找到了两者最优的对应关系。WDTW模型传递方法基于此对应关系建立了回归方程。

3.3 模型传递前后平均光谱差异对比

使用WDTW模型传递方法建立了近红外光谱校正模型，分别对比汽油、药品样本的原始平均光谱图和传递后的平均光谱图，可以很直观地看出WDTW模型的传递效果。图5左图为汽油样本测试集原始的平均光谱图，其中实线代表源机光谱，点线代表目标机光谱，可见源机光谱和目标机光谱之间有很大差异，如果直接将源机光谱建立的回归模型用于目标机光谱数据，必将产生很大误差。图5右图为汽油样本测试集模型传递后的平均光谱图，源机光谱和目标机光谱之间基本重合，可见WDTW模型传递方法减小了两仪器间光谱的差异。药品样本测试集校正前后的谱图效果与汽油样品相同，即采用WDTW模型传递方法于药品样本可有效减小源机和目标机光谱之间的差异。

表2记录了WDTW算法模型传递前后源机和目标机仪器间的光谱差异以及光谱校正率。以汽油C7成分为例，在目标机光谱取标样集为5，校正前光谱平均差异为0.018 8，校正后光谱的平均差异为0.002 1，通过校正率计算公式(12)可得其光谱校正率为98.75%，WDTW模型传递算法可以很好的消除源机光谱和目标机光谱之间的差异。

表2 WDTW算法传递前后测试集主从仪器间的光谱差异Table 2 Spectral difference between master and slave instruments before and after test set passing through WDTW algorithm

将WDTW算法与DS、PDS和DTW 3种算法模型传递后以及未校正的光谱平均差异进行对比，如表3所示，对比未校正的光谱，发现经过模型传递后的平均光谱差异均有减小，其中WDTW算法表现最好，模型传递后的平均光谱的差异整体最小，也说明此算法模型的传递效果最好。其中汽油样本经过WDTW算法模型传递后，C7、C8、C9和C10 4种成分的平均光谱差异均为0.002 1。在药品样本数据集中，WDTW算法同样取得了很好的效果，药品活性、重量和硬度在模型传递后，平均光谱差异分别为0.026 8、0.025 7、0.026 0。

表3 DS、PDS、DTW和WDTW传递后测试集主从仪器间的光谱差异Table 3 Spectral differences between master and slave devices of test set after DS,PDS,DTW and WDTW transmission

3.4 模型传递前后预测结果对比

为了进一步验证WDTW模型传递方法的有效性，实验采用源机光谱数据建立PLS回归模型，将经过模型传递后的光谱数据代入回归模型，并用预测标准偏差评价模型传递的效果。通过与DS、PDS、DTW算法比较(PDS算法的窗口大小用交叉验证确定),WDTW算法传递后，光谱预测结果的预测标准偏差整体最小，模型传递效果最好，结果如表4所示。

表4 DS、PDS、DTW和WDTW传递后测试集预测标准偏差(SEP)Table 4 Prediction standard deviation(SEP) of the test set after DS,PDS,DTW and WDTW pass

实验结果显示，在模型传递前，未校正的目标机光谱直接带入源机光谱数据建立的PLS回归模型进行预测有较大误差。通过模型传递之后，校正后的光谱数据预测标准偏差均有减小，其中WDTW算法整体的预测标准偏差最小，表明此时所建模型预测能力最好。DS算法通过全光谱建立校正模型，范围太大容易引起误差。PDS算法通过限定一个窗口，在窗口范围内建立校正模型，窗口大小将会直接影响结果的好坏。DTW算法很好的找到源机光谱和目标机光谱的最佳对应关系，其建立的模型相较DS和PDS算法在汽油数据集上有明显的优势，但在药品数据集中表现较差、原因可能是由于未经过数据的预处理，光谱中的一些噪声会引起DTW算法建立的模型不稳定。而WDTW算法在经过小波变换预处理后，有效消除了样本数据间存在的基线漂移、随机噪声等，之后通过找到波长点最佳对应关系，再根据此对应关系建立校正回归模型，可以更加精确的求得转移光谱，实现近红外光谱模型的传递。

4 结 论

本文针对近红外光谱模型传递问题，提出了基于小波变换动态时间规整的模型传递方法。该方法首先使用小波变换进行光谱预处理，然后利用动态时间规整算法求得传递光谱。通过一维小波连续变换对近红外光谱进行多尺度细化分析，以实现对光谱数据平滑、降噪以及消除基线漂移的处理。在经过小波变换预处理后，动态时间规整算法通过寻找不同光谱仪器之间波长点最佳的匹配关系并建立回归模型，求得传递矩阵实现模型传递。对汽油和药品两个数据集进行实验，证明了该方法能够消除不同仪器间光谱的平均差异，提高模型的预测能力和通用性。通过与其他模型传递方法DS、PDS和DTW算法进行对比后，表明基于WDTW算法所求的参数少，传递效果最好，整体性能最优。