, 芦鹏程, 刘晓东, 齐天舒 东北大学秦皇岛分校 计算机与通信工程学院, 河北 秦皇岛 066000
收稿日期：2020-09-09
基金项目：国家自然科学基金青年基金资助项目(61601104)。
作者简介：赵煜辉(1971-)，男，河北秦皇岛人，东北大学秦皇岛分校教授。
摘要：针对近红外光谱数据的维度高、特征之间存在严重的多重共线性的特点, 提出了无迁移标准的通过校正分布差异的标定迁移方法(calibration transfer via correcting distributions difference, CT-CDD).CT-CDD首先建立主仪器的偏最小二乘模型, 然后通过偏最小二乘模型提取主仪器和从仪器的潜变量，并且分别对主仪器和从仪器的潜变量进行聚类.该方法基于这样的假设: 聚类后的主仪器和从仪器的每一部分特征光谱均服从单高斯分布.最后, 找到2个仪器的最接近的子分布, 通过校正均值和方差来校正数据分布的差异.实验结果表明CT-CDD通常更加鲁棒并且还可以实现最低的均方根预测误差.
关键词：近红外光谱迁移标准标定迁移偏最小二乘模型分布差异
Research on Calibration Transfer via Correcting Distributions Difference
ZHAO Yu-hui
, LU Peng-cheng, LIU Xiao-dong, QI Tian-shu School of Computer and Communication Engineering, Northeastern University at Qinhuangdao, Qinhuangdao 066000, China
Corresponding author: ZHAO Yu-hui, E-mail: 1000272@neuq.edu.cn.
Abstract: Aiming at the characteristics of high dimensionality of near-infrared spectroscopy data and serious multi-collinearity between features, a method of calibration transfer via correcting distributions difference(CT-CDD) was proposed without transfer standards. CT-CDD firstly establishes PLS(partial least square) model of the master instrument, and then latent variables of both the master instrument and the slave instrument are extracted by the PLS model. Next, the latent variables of the two instruments are clustered. The method is based on the assumption that the characteristic spectra of each part of the master and slave instruments after clustering is a single Gaussian distribution. Finally, the nearest sub-distribution of the two instruments is found, and the differences in data distribution is corrected by correcting mean and variance. The experimental results show that CT-CDD is generally more robust and can also achieve the lowest RMSEP (root mean squared error on prediction).
Key words: near-infrared spectroscopytransfer standardcalibration transferpartial least square (PLS) modeldistributions difference
近红外光谱分析技术具备仪器操作简单、分析数据速度快、成本较低、不污染样品等优势, 已在各领域得到了普遍应用, 如农产品生产、化工产品生产、食品生产以及环境监测领域[1-4].近红外光谱分析技术在定性分析和快速物质成分定量以及实现在线检测方面具有独特优势[5], 这主要由于化学计量学方法的不断发展和商业化高精度光谱仪器的日益成熟.
近红外光谱分析技术的重要内容是多元校正模型的建立[6-7]，即通过使用一定的数学分析方法, 对近红外光谱数据进行分析建模, 从而达到对一些指标进行预测的目的.但是, 在实际的工业生产中, 测量仪器和测量条件往往并不稳定或是前后不一致.通过对原有近红外光谱数据进行分析建立的模型往往对新的数据集并不适用, 从而导致原有模型失效.如果在测量环境变化后, 便因为原有模型对新的测量环境测得的数据不适用, 而抛弃原有的多元校正模型, 选择重建新模型, 往往会耗费巨大的时间和空间的资源以及造成对原有学习成果的浪费.在此种情况下, 标定迁移具有重大的意义.标定迁移是指在不同测量仪器或测量状态下的多元标定模型迁移方法, 通过将从光谱数据迁移至主光谱数据空间, 进而实现主光谱数据模型对从光谱数据模型的预测, 避免重复建模[8].
目前应用较为广泛的标定迁移方法有斜率截距校正(slope bias correction, SBC)算法、直接标准化(direct standard, DS)方法[9]和分段直接标准化(paragraphing direct standard, PDS)方法、典型相关分析标定迁移(canonical correlation analysis calibration transfer, CCACT)、多元散射校正(multiple scattering correction, MSC)[10]以及迁移成分回归(transfer component regression, TCR).其中, DS和PDS方法通过一组标准样本来纠正主仪器和从仪器之间光谱的差异.在DS中, 主仪器的每个波长与从仪器的所有波长相关.在PDS中, 主仪器的每个波长与从仪器的波长窗口相关, 并且最终基于每个窗口的回归系数形成带状迁移矩阵.实验结果与假设是一致的, 即在各种迁移方法中, 主仪器和从仪器之间的频谱相关性被限制在较小的区域.PDS的关键是窗口大小的选择和标准样本数目的确定.CCACT方法指对主从仪器的标样提取特征.MSC方法对光谱数据进行预处理, 会简化数学模型, 让模型具备更强的解释性.然而在实际的工业生产中, MSC预处理方法并不能显著提高模型的预测能力[11]; SBC算法为一种单变量方法, 因此在测量仪器和测量条件变化引起系统化的光谱差异的情况下, 才能取得较好的效果.现实生活中, 光谱差异往往比较复杂, 此时它的预测能力是不确定的; TCR是一种无标准的迁移方法, 它结合了迁移成分分析TCA和普通最小二乘法(OLS).TCA的基本思想是在再生希尔伯特空间中投影2个仪器的数据, 在这个空间中, 主仪器和从仪器的数据分布尽可能地接近, 同时保留原始数据的关键属性.TCR是一个具有良好泛化能力的稳健模型, 但无法实现更准确的预测.
对此, 本文提出了一种新颖的没有迁移标准的标定迁移方法: 通过校正PLS[12]子空间的数据分布差异(CT-CDD)的方法进行标定迁移.该方法首先建立了主仪器的PLS模型, 同时将主仪器和从仪器的光谱投影到PLS子空间.其次, 将不同光谱的潜变量分别进行聚类分析, 同时利用普通最小二乘法建立主仪器潜变量与浓度信息之间的回归模型.此外, 找到2台仪器之间的对应的数据分布最接近的特征光谱, 并分别计算转换函数.最后, 预测结果可以通过各自的转换函数得到校正.通过2个NIR数据集, 将CT-CDD的性能与SBC, CCACT, PDS, TCR和MSC的预测性能进行比较.
1 理论知识1.1 偏最小二乘法在化学计量学中, 偏最小二乘(PLS)算法是一种很有效的多元标定方法.PLS算法结合了多元线性回归、主成分分析、典型相关分析的优点, 通过分数向量建立输入空间和响应空间之间的关系.PLS模型的目的是确保最佳的潜变量数量.潜变量是原始变量的线性组合，它包含了关于X和y之间关系的最大相关信息.在数学上, 由以下目标函数表示:
 | (1) |
假设一个PLS模型被建立在光谱矩阵X∈Rn×p和浓度矢量y∈Rn×1之间, 该模型称为PLS1(n表示样本数, p表示原始光谱的维度).在这个算法中, 第一个权重向量必须是矩阵XTyyTX的主要的特征向量.从第二个潜变量开始, 它要求随后的潜变量与前面的潜变量正交(不相关).因此, 此后的权重向量也是矩阵的主要特征向量, 重复这一系列步骤直到收敛.PLS1模型通过如下等式构建:
 |
最后, 模型的回归系数β表示如下:
 | (2) |
1.2 校正分布差异的标定迁移方法一个域由两个主要部分组成: 输入空间R, 其对应的边际概率分布为P(R).给定主仪器的有标签光谱{Xm, ym}和从仪器的无标签光谱{Xs}.本文的目标是预测从仪器的测试数据输出{ys}.不同仪器测量的光谱不同, 因此两种仪器之间的数据分布是不同的.
相对熵或KL散度可用于表示2个数据分布的距离.KL散度的公式写成如下形式:
 | (3) |
 | (4) |
 | (5) |
 | (6) |
降维后, 2个仪器的数据分布之间的距离为
 | (7) |
然后需要找到一个转换函数, 使主光谱和从光谱之间的距离最小.当主光谱和从光谱之间的距离为零时, 确定最佳转换函数.由于每种仪器的数据分布是混合高斯分布, 因此需要对数据进行聚类.聚类之后的数据近似为单高斯分布.此外, 当2个仪器数据分布之间的距离最近时, 执行转换函数, 进行迁移.目标函数表述为如下形式:
 | (8) |
解决问题的关键是求解目标函数的转换函数.由于聚类之后的数据分布近似为单高斯分布, 主光谱和从光谱的数据分布分别为tim~N(μim, Σim)和tis~N(μis, Σis).高斯分布由参数均值和协方差确定，然后分别对均值和协方差进行校正.均值的校正十分容易实现，协方差的校正需要在协方差纠正之前, 对数据进行中心化.2个仪器的数据分布分别为ti_cenm~N(0, Σim)和ti_cens~N(0, Σis).对于协方差, 通过矩阵M执行转换函数φ(ts)=Mts, M∈RA×A(即协方差变换), 使得校正的从光谱与主光谱之间的距离最小化.目标函数最终表达如下:
 | (9) |
1) 由于聚类后的每组数据近似正态分布, 因此数据的均值为0, 概率密度函数Pi(ti_cenm)和Qi(ti_cens)由式(10)给出:
 | (10) |
 | (11) |
 | (12) |
 | (13) |
 | (14) |
2 算法建模CT-CDD算法的实现过程如下:
输入: 主仪器的光谱和对应的样本浓度信息(Xm, ym), 从仪器无标签的光谱数据Xs和从仪器的测试集(Xtests, ytests).
输出: 从仪器测试集的预测值?tests.
步骤1??在主仪器的标定集(Xcalm, ycalm)上建立PLS模型, 计算出PLS模型的回归系数βm, 权重矩阵Wm和载荷矩阵Pm.
1) 使用KS算法划分数据集:
 |
 |
 |
1) 计算出主仪器的特征光谱和预测值(Tm, ?m)以及从仪器的伪特征光谱
 |
 |
3) 使用k-means将主、从仪器的光谱分别划分为K份:
 |
4) 对主仪器聚类后的子光谱特征, 分别建立最小二乘模型.
for i in range(K):
 |
 |
步骤3??预测阶段: 将从仪器的待测样本转换至主仪器的模型, 使用主仪器模型对其进行预测.
1) 从仪器测试集的伪特征和伪预测值(
 |
 |
 |
 |
3.1 数据集介绍玉米数据集包含在3台NIR光谱仪(M5, MP5和MP6)上测量的80个样品.其中, 这3台不同的红外光谱仪因其工作原理不同, 所以得到的近红外光谱略有差异, 但对绝大多数谷物而言, 仪器的工作原理不同所产生的误差并不会影响实验结果, 所以本文采用这3台仪器测量的80个玉米样本的近红外光谱作分析.玉米数据集中每个样品含有4种成分: 水分、油、蛋白质和淀粉.波长范围为1 100~2 498 nm(700点), 间隔为2 nm.该数据集可以从http://www.eigenvector.com/Data/Corn/下载.本文使用3台近红外光谱仪和水分含量进行研究讨论.不同仪器之间的光谱差异见图 1.由图 1可知，不同仪器测量的同一份谷物样本光谱是有差异的.
图 1(Fig. 1)
 | 图 1 不同仪器之间的光谱差异Fig.1 Spectral differences between different instruments (a)-M5与MP5；(b)-M5与MP6；(c)-MP5与MP6. |
3.2 数据处理对于玉米数据集, 通过Kennard-Stone(KS)算法将80个玉米样本分成两组: 其中第一组是64个, 被用作标定集的样本; 第二组是16个, 被用作测试集的样本.主仪器的PLS模型建立在标定集上.对于有迁移标准的迁移方法, 使用Kennard-Stone方法在标定样本上选择若干个标准样品.
3.3 性能评估在该实验中, 均方根误差RMSE被用作参数选择和模型评估的指标.此外, RMSEC表示标定集的训练误差, RMSECV表示交叉验证误差, RMSEP表示测试集的预测误差. RMSE计算方法为
 | (15) |
偏差(bias)和标准误差(SE)也被用作模型评估的参考指标.bias和SE的表达式如下所示:
 | (16) |
为了更加直观地比较CT-CDD和其他标定迁移方法在预测性能上的差异程度, 另一个重要参数改善率(h)被引用来进行比较, 定义如下:
 | (17) |
3.4 实验结果对比及分析本文选用玉米的光谱数据集检验模型的性能.使用SBC, PDS, CCACT, TCR和MSC作为比较方法, 进行对比实验.对于SBC, PDS, CCACT, MSC和TCR算法均采用PLS算法作为主体算法, 使用主仪器的光谱数据建立多元标定模型作为参考模型, 用于对从仪器的待测样本进行预测.实验结果主要包含以下几个部分:
1) 主仪器PLS模型的潜变量数目的选取过程以及RMSEC, RMSEP, RMSECV的结果展示;
2) CT-CDD在迁移前后预测结果的拟合能力示意图;
3) CT-CDD和其他五种迁移方法的RMSEP比较;
4) CT-CDD和其他五种迁移方法的预测结果的拟合能力示意图.
玉米数据集实验结果分析如下:
在玉米数据集上, M5, MP5和MP6三台仪器的PLS模型的RMSECVmin, RMSEC, RMSEP和LV如表 1所示.仪器M5的PLS模型的RMSECVmin, RMSEC, RMSEP分别为0.010 66, 0.005 99和0.007 64.可以看出三个均方根误差相差不大, 因此PLS模型比较稳定, 不存在过拟合和欠拟合的现象.仪器MP5的PLS模型的RMSECVmin, RMSEC, RMSEP分别为0.130 35, 0.094 58和0.124 45.和M5的PLS模型类似, 没有欠拟合和过拟合的现象.在MP6的PLS模型上再一次得到了相同的结论.
表 1(Table 1)
 表 1 PLS模型相关结果的总结Table 1 Summary of relevant results of the PLS models
| 表 1 PLS模型相关结果的总结 Table 1 Summary of relevant results of the PLS models | ||||||||||||||||||||||||||||||
参数的选择通过10折交叉验证进行选取, 并且基于最低的RMSECV准则来确定最佳的潜变量数量.PLS模型中最佳潜变量的数量分别为14, 15和10.主仪器能够建立一个预测性能更优的模型是十分重要的, 因此, 选择具有良好的预测性能的仪器作为主仪器.从表 1可知, 仪器MP6的预测误差大于仪器MP5, 仪器MP5的预测误差大于仪器M5.因此, 选择两种组合(M5*-MP5和M5*-MP6)进行模型检查，上标*表示主仪器, 另一个表示从仪器.
为了证明所提出的标定迁移方法(CT-CDD)的预测性能, 将CT-CDD与SBC, PDS, CCACT, TCR和MSC五种标定迁移方法进行比较.在CT-CDD中, 聚类的数目通过10折交叉验证确定.由于玉米数据集含有80个样本, 因此聚类后的子模型的最大数目设置为3.否则, 所计算的迁移矩阵是欠秩的, 导致最后的预测结果将无穷大.此外由于样本数目的限制, 当聚类数目较大时, 聚类后的特征光谱也没有足够的样本建立一个稳定的模型.在图 2中, 可以看到仪器M5的聚类数目为2时, 对应的RMSECV是最低的.因此, 选择聚类数目为2, 做k-means的初始参数.在PDS中, 选择PLS作为迁移函数.PDS的窗口大小通过5折交叉验证进行选取, 范围在3到16之间, 间隔为2，PDS的最佳窗口大小如表 2所示.可以看到M5作为主仪器, MP5作为从仪器时, PDS的最佳窗口大小均为15.当M5作为主仪器, MP6作为从仪器, 标准样本数目为5时, 最佳窗口大小为3.其余的PDS的最佳窗口大小均为15.在TCR中, 搜索最佳子空间的维数从1到24，子空间的最佳维数分别为6和4.
图 2(Fig. 2)
 | 图 2 由M5确定的不同子模型数量的交叉验证误差Fig.2 Cross validation errors for the number of different sub-models determined by M5 |
表 2(Table 2)
表 2 MSC, SBC, PDS, CCACT, TCR, 和CT-CDD六种迁移方法在玉米数据集下的RMSEPTable 2 RMSEP of corn datasets with SBC, PDS, CCACT, TCR, MSC and CT-CDD
| 表 2 MSC, SBC, PDS, CCACT, TCR, 和CT-CDD六种迁移方法在玉米数据集下的RMSEP Table 2 RMSEP of corn datasets with SBC, PDS, CCACT, TCR, MSC and CT-CDD | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
CT-CDD和其他5种标定迁移方法的预测误差见表 2.对于仪器MP5到仪器M5的光谱转移, 当标准样品的数量是35时, SBC达到最低RMSEP(0.288 72)；当标准样品数为5时, PDS达到最低的RMSEP(0.188 28)；当标准样品数为25时, CCACT达到最低RMSEP(0.186 99), CT-CDD的RMSEP小于PDS, SBC和CCACT这三种方法的预测结果的最小值.此外, CT-CDD的RMSEP也小于TCR(0.473 91)和MSC(2.078).对于从MP6到M5的光谱转移, CT-CDD的预测性能也是优于其他5种对比方法.从这两组实验可以看出CT-CDD模型在通常情况下能够取得最优的预测性能, 并具有更好的鲁棒性.
此外, 图 3和图 4显示了6种不同的标定迁移方法获得的预测值与测量值的关系图.预测浓度和测量浓度之间的零差异, 将会使得样本点在直线上.对于有标准样本的标定迁移方法, 在不同的标准样本下, 当预测性能最优时, 选择该组实验用于比较, 以便更加充分地体现出CT-CDD能够取得良好的预测性能.
图 3(Fig. 3)
 | 图 3 6种方法在仪器M5和仪器MP5之间预测结果的散点图Fig.3 Scatter plots for prediction between instruments M5 and MP5 using six methods |
图 4(Fig. 4)
 | 图 4 6种方法在仪器M5和仪器MP6之间预测结果的散点图Fig.4 Scatter plots for prediction between instruments M5 and MP6 using six methods |
从图 3中可以看出, CT-CDD的样本点更加接近直线, MSC，TCR和SBC在该组实验下, 拟合效果较差.同样, 在图 4中, 对于仪器MP6到仪器M5的光谱传输, 可以看到CT-CDD通常比其他五种方法更接近直线, 这再次确认了CT-CDD能够达到最佳预测性能.实验结果表明, CT-CDD能够在所有模型中实现最佳的预测性能, 同时具有更好的泛化能力.
4 结语本文研究了不同仪器上数据分布差异之间的关系, 并提出了一种无标准的通过校正PLS子空间数据分布差异(CT-CDD)的标定迁移方法.该方法试图找到一个转换函数, 使得当从仪器的数据投射到这个空间时, 主仪器和从仪器之间的数据分布距离可以大大减小.由于数据服从非单高斯分布, 因此需要对光谱进行聚类, 并通过各自的转换函数最小化两个仪器之间的每个子分布的距离.本方法中值得注意的是, 在相同的PLS子空间中能够保留2个仪器的重要属性并消除光谱的多重共线性.同时主仪器的特征和从仪器的伪特征之间的数据差异可以更精确地缩小.然后通过校正来自不同仪器的潜变量的每个部分的均值和方差来进一步校正数据分布上的差异.
使用SBC, PDS, CCACT, TCR和MSC作为对比实验, 使用2个NIR数据集检验CT-CDD方法的性能, 并且所提出的方法(CT-CDD)通常实现了最佳的RMSEP.结果清楚地表明, CT-CDD能够成功地用于校正在不同仪器上测量的光谱之间的差异.对于SBC, PDS和CCACT这三种迁移方法, 它们需要标准样品建立迁移模型.在TCR中, 从仪器样品还需要少量的参考值.这两个条件在实际应用中, 都会产生很昂贵的代价，甚至无法满足这一条件.因此, 当标准样品在实际应用中不可获得时, CT-CDD是一种有效的标定迁移方法.
参考文献
| [1] | Aryal G H, Hunter K W, Huang L M. A supramolecular red to near-infrared fluorescent probe for the detection of drugs in urine[J]. Organic & Biomolecular Chemistry, 2018, 16(40): 7425-7429. |
| [2] | Ali R, Ahmadi N H, Mohammad K. A case-study of NIRS application for infant cerebral hemodynamic monitoring: a report of data analysis for feature extraction and infant classification into healthy and unhealthy[J]. Informatics in Medicine Unlocked, 2018, 11: 44-50. DOI:10.1016/j.imu.2018.04.001 |
| [3] | Fukuda M. Optical topography as an auxiliary laboratory test for differential diagnosis of depressive state: clinical application of near-infrared spectroscopy (NIRS) as the first trial for approved laboratory tests in psychiatry[J]. Seishin Shinkeigaku Zasshi=Psychiatria et Neurologia Japonica, 2015, 117(2): 79-93. |
| [4] | Xue J T, Liu Y F, Ye L M, et al. Rapid and simultaneous analysis of five alkaloids in four parts of Coptidis Rhizoma by near-infrared spectroscopy[J]. Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy, 2018, 188(1): 611-618. |
| [5] | Harris C, Despa M, Kelly K. Design and fabrication of a cross flow micro heat exchanger[J]. Journal of Microelectromechanical Systems, 2002, 9(4): 502-508. |
| [6] | Giraudo A, Grassi S, Savorani F, et al. Determination of the geographical origin of green coffee beans using NIR spectroscopy and multivariate data analysis[J]. Food Control, 2019, 5: 137-145. |
| [7] | Shi Y Y, Li J Y, Chu X L. Progress and applications of multivariate calibration model transfer methods[J]. Chinese Journal of Analytical Chemistry, 2019, 47(4): 479-487. DOI:10.1016/S1872-2040(19)61152-7 |
| [8] | Li J, Yu X N, Ge W Z, et al. Qualitative analysis of maize haploid kernels based on calibration transfer by near-infrared spectroscopy[J]. Analytical Letters, 2019, 52(2): 249-267. DOI:10.1080/00032719.2018.1459656 |
| [9] | Zimmerman N, Presto A A, Kumar S P N, et al. A machine learning calibration model using random forests to improve sensor performance for lower-cost air quality monitoring[J]. Atmospheric Measurement Techniques, 2018, 11(1): 291-313. DOI:10.5194/amt-11-291-2018 |
| [10] | McGuire T, Harrison J A. Direct standardization[J]. International Journal of Urology, 1995, 5(4): 398-400. |
| [11] | Geladi P D, Macdougall D B, Martens H. Linearization and scatter-correction for near-infrared reflectance spectra of meat[J]. Applied Spectroscopy, 1985, 39(3): 491-500. DOI:10.1366/0003702854248656 |
| [12] | Li H D, Xu Q S, Liang Y Z. LibPLS: an integrated library for partial least squares regression and discriminant analysis[J]. Chemometrics and Intelligent Laboratory Systems, 2018, 176: 34-43. DOI:10.1016/j.chemolab.2018.03.003 |
