采样误差分布分析法建立氯丙醇含量的稳健近红外光谱分析模型

陈万超 陶鑫 范长春 张飞宇 杜一平

摘 要 利用近红外光谱技术结合采样误差分布分析（SEPA）方法建立了二氯丙醇产品生产过程的氯化液杂质3-氯-1，2-丙二醇浓度的分析模型。对样本数据进行1000次随机划分，建立1000个子模型，获得多个潜变量数下的交互检验误差，进行统计分析。绘制了误差分布图，计算其中位数、标准偏差、偏斜度和分布峰度等统计指标，通过这些指标的综合分析对近红外光谱分析模型进行条件优化、建模和模型评价等。4种光谱处理方法显示出比较理想的模型性能，作为候选与不同波长区域的选择相结合，继续运用SEPA运算，进一步优化模型。最终优化的建模条件为：一阶导数结合标准正态变换; 6931～6017 vcm1波数区间;使用5个偏最小二乘潜变量。校正、交互检验和独立验证误差分别为0.881%、1.282%和1.167%。所选择的波长具有可解释性，模型的各项统计参数合理、可信。研究结果表明，SEPA能全面、合理地考察多项统计指标，可以建立实用、稳健的近红外光谱分析模型。

关键词 采样误差分布分析; 近红外光谱; 建模; 稳健

1 引 言

近红外光谱分析技术具有简单、快速、无损、准确等优点，近年来应用广泛。应用近红外光谱法时，通常需要借助化学计量学建立光谱多元分析模型，实现定性和定量分析，在实际应用时，建模常是难点，也容易出现错误。近红外光谱分析建模一般包括条件优化（光谱预处理方法选择和波长选择等）、建模和模型评价等部分，这些工作均离不开数据集合的划分，如将数据分为校正集、交互检验集和独立验证集等，不同的集合划分常获得不同的建模结果。近年来，集合划分方法及其相关研究备受关注[1～5]。划分方法大体可分为人为划分和随机划分两大类，前者根据指标性质数值或光谱之间距离[6]，或者从指标和光谱相关关系考虑样品之间的差异[2]划分样本集;后者通过多次随机采样，达到选择具有代表性样本数据集的目的[1]。对于复杂的近红外光谱数据，经常包含特殊数据，如指标性质数值相差不大，或光谱却有较大差异。这种数据采用人为划分方法难以得到理想的结果，而随机方法，只要样本采样次数足够大，各种情况的数据都能被考虑到，因此随机划分方法极具竞争力。Xu等[1]提出的蒙特卡罗交互检验法（MCCV）很好地利用了这种随机采样的特性，多次随机划分数据集，建立多个子模型，对这些子模型的误差取平均，用于交互检验。Li等[7]提出的模型集群分析方法（MPA）采用类似的策略建立多个子模型，并考虑了集群的子模型多种统计性质用于建模。本研究组在MPA算法的基础上，近期提出了采样误差分布分析（SEPA）算法[8]，利用多种统计指标综合分析随机采样误差，如误差分布的中位数和标准偏差、分布偏斜度和分布峰度等，用于近红外光谱分析中的异常点筛查、光谱预处理方法和波长选择、交互检验、模型评价，还用于模型转移[9]和波长选择[10]等。多次随机采样结合多种统计指标的综合考虑，能更全面地考察数据的各种状况，在此基础之上的条件优化和模型建立及评价可以更加合理，所建立的模型更加稳健。

二氯丙醇化工生产中的关键副产物一氯丙二醇对产品质量控制至关重要。由于主产物和被测物质结构非常相似，同时样品中还含有含量较小的其它化合物，所测的近红外光谱中难以获得被测组分清晰的特征谱峰，光谱干扰严重。本研究以此副产物为研究目标，采用近红外光谱进行分析研究， 利用SEPA方法对近红外光谱数据进行了细致研究，建立了合理、稳健的分析模型。

2 实验方法与理论

2.1 仪器与数据

样本取自江苏省某化工企业二氯丙醇生产工段的氯化液，样品为液体。用Antaris II近红外光谱仪（美国赛默飞公司）测定样品的透射光谱，比色皿光程为8 mm，波数范围3799～11999 cm1，分辨率4 cm1，去除吸光度饱和的数据，实际使用的光谱范围为5361～11999 cm1。扫描次数为16，取平均光谱。3-氯-1，2-丙二醇的浓度采用GC7820气相色譜仪（美国安捷伦公司）进行测定，方法参考国标法GB/T 13097-2007[11]。

2.2 采样误差分布分析（SEPA）

利用SEPA[8]对近红外光谱数据进行光谱预处理、波长选择、建模和模型评价，均以交互检验的形式进行工作。因此，在此以交互检验过程为例说明SEPA。用X表示光谱，维度nm，y表示指标浓度，维度n1，n和m分别表示样本数和波长数，随机采样次数为N，最高潜变量数为Ncv，建模方法为偏最小二乘法（PLS）。首先通过随机采样将样本集合划分为校正集和交互验证集，采样N次，获得N个不同的样本集合划分;然后分别对每个样本集合用校正集建立校正子模型，并针对各个潜变量数（nLVs）计算交互验证预测均方根误差（RMSECV）;对所有得到的N个子模型的RMSECV 误差数据（共计NNcv）进行统计分析，在每个潜变量数下，绘制误差分布图，并计算其误差中位数、标准偏差（STD）、分布偏斜度和分布峰度;最终综合考虑上述统计结果， 进行条件选择和模型评价。

采用中位数而非平均值估计RMSECV，是因为前者更加稳健，标准偏差STD能反映误差分布宽度，其值越小， 说明采样偏差越小，对建模越有利;而分布偏斜度和分布峰度则是判断误差分布的合理性，前者反映分布的对称性，理想值为0，后者表示分布的峰形态（简单理解就是峰太瘦或太胖的程度），其理想值为3。由于采用随机选择样本集划分，并多指标综合考虑误差分布，因此所选择的参数和建立的模型更加合理和稳健。

2.3 数据分析

样本总数为83，被测指标为3-氯-1，2-丙二醇的含量，光谱有1722个波长点。采用Kennard -Stone（K-S）方法[7]将83个样本划分为校正集和验证集，前者用于条件优化和建模，后者用来对模型进行独立的检验，以此设置校正集68个样本，验证集15个样本。条件优化过程均采用交互检验的方法，在建模集68个样本中，随机选取20个样本组成交互检验集，其它样本进行建模（称为校正集）。在进行SEPA运算时，随机采样次数为1000。表1列出了被测指标3-氯-1，2-丙二醇含量的相关数据。所有计算均在Matlab（R2010a）环境下自编程序完成。

3 结果与讨论

3.1 近红外光谱

在所测得的近红外光谱中，低波数的吸光度较高，但噪声也高，而高波数位置，除了在8600 cm1附近有一个小峰外，基本无吸收。明显的吸收峰包括8600、6900、6300、6000和5800 cm1（这些谱峰在建模中的作用参见3.5节）。

3.2 采用原始数据初步建模

SEPA的核心是多次随机采集校正集和交互检验集，建立多个子模型，通过对子模型的综合分析进行条件优化和建模。为了说明SEPA的核心工作过程和数据分析方法，用建模集的68个样本的原始光谱数据，采用SEPA方法初步建模。交互检验集样本数为20，校正集样本数为48，随机采集这些样本次数800，即可建立800个子模型，每个子模型都产生各个潜变量数（nLVs=1～15）的一组交互检验误差，对每个潜变量数的800个误差进行统计分析，绘制误差分布图，在每个潜变量数下，取其中位数作为交互检验误差，对nLVs作图。

在潜变量数为8时交互检验误差最小，其标准偏差也最小。此时的误差分布均匀，相对比潜变量数少于8时，更能体现出正态分布的轮廓，作为比较给出4个潜变量的误差分布，其分布不如8个潜变量的误差分布。但8个潜变量的误差分布左右略不对称，计算发现其偏斜度为1.479，峰度为12.32，距的理想值0和3相差较远。这些不足可以通过光谱预处理和波长选择进一步得到改善（参见3.4和3.5节）。

3.3 随机采样次数的选择

数据集划分通过随机采样完成，随机采样次数分别选择600、800、1000、1200，按3.2节的方法进行交互检验，重复计算5次，考察结果的稳定性。结果表明，采样次数为600时，潜变量数有时选择8，有时选择9;采样次数为800时，潜变量数多选择8，偶尔选择9;当采样次数为1000或更高时，潜变量数均选择8。而重复计算误差的相对标准偏差RSD（%）分别为0.7578、0.8845、0.1361和0.2366。因此，选择随机采样次数为1000，计算时间约20 s。

3.4 光谱预处理方法的选择

光谱预处理可以改善模型精度，而不同的方法常会导致预测误差和最佳潜变量数均发生明显的变化。常见的光谱预处理方法[12]包括Savitzky-Golay平滑（用S表示）、多元散射校正（MSC）、标准正态变换（SNV）、一阶导数（1D）、二阶导数（2D）以及它们的组合。考察了25种有效的单一方法及其组合，分别进行SEPA计算。不同光谱预处理方法的交互检验误差随潜变量数变化图，因为有些方法得到的变化图非常相似，图中只绘制了其中的一种。有4种方法（SNV、1D、1D+S和1D+SNV）的误差较低，而且潜变量数也较低。它们的交互检验误差图及标准偏差（误差棒），SNV的预测误差最低，为1.391%，但潜变量数较高（nLVs=7），另外3种方法的预测误差稍高，但可以选择更低的潜变量数，如果选择4，预测误差分别为1.568%、1.552%和1.495%。 本研究将这4种方法作为候选方法，通过波长选择进一步优化模型。

3.5 波长选择

近年来，发展了很多波长选择方法[13～20]，如移动窗口偏最小二乘法[13]、无信息变量剔除法[14]、竞争性自适应加权抽样法[15]和基于集群分析变量提取法[16]等，各有优缺点。为了简化，本研究采用最简单的相关系数法选择候选波长区域，即计算浓度与各波数处吸光度的相关系数R。各个波数吸光度与浓度的相关系数图，单一波数的相关系数R都不高，按R>0.4选择3个波数区间： R1=6931～6017 cm1、R2=5959～5897 cm1、R3=5758～5523 cm1。光谱中8600 cm1峰因为强度太低，在波長选择中没有体现。5800 cm1位置因为干扰物也同样具有强吸收，在波长选择中也没有体现。波长选择的R2和R3谱带为5800 cm1两边区域，避开了强干扰区域。R1主要反映6900、6300和6000 cm1 3个谱带。

被测组分3-氯-1，2-丙二醇在4769、4421和3995 cm1具有明显的吸收峰[21]，但样品中含有其它含氯丙醇，含量较高的是二氯丙醇，对被测组分构成比较严重的光谱干扰。3-氯-1，2-丙二醇在7000～5700 cm1范围有一个很宽的吸收带，而主要干扰物二氯丙醇除在7000和5900 cm1处和5700～5800 cm1范围有明显吸收峰，其它区域吸收较弱[21]。3个谱带正是反映了被测组分和干扰物之间的光谱差异。在R1（绿色）、R2（橙色）和R3（蓝色）波数范围，3-氯-1，2-丙二醇的吸收比较强，而二氯丙醇的吸收比较弱，尤其是R1带。

将3.4节中的4种光谱预处理方法与所选择的波长区间组合，利用SEPA综合考察模型性能。表2给出了模型的各项指标，以及在各种优化条件下，建模后采用验证集计算预测误差RMSEP。

与原始光谱相比，光谱预处理后模型得到明显改善，最佳方法为SNV方法，在7个潜变量下，模型的建模集和验证集误差RMSEC和RMSEP为0.985%和1.195%。进一步考虑波长选择，当使用3个波长区段建模时，各种误差并未降低，原因可能是在R2和R3段光谱干扰依然很严重。R1区域的光谱干扰相对比较弱，如果只选择R1进行建模，模型误差进一步降低。最佳结果是1D+SNV结合R1，在5个潜变量下，RMSEC和RMSEP降为0.881%和1.167%。1000次计算交互检验误差的分布图，4种方法的分布都比较均匀，分布的对称性也得到显著改善。最佳方法1D+SNV+R1的偏斜度为0.4726，峰度为3.9573，与原始光谱建模的1.4789和12.3246（见3.2节）相比，更接近理想值的0和3。另一外值得指出的是，标准偏差STD在选用不同的处理方法和波长时变化不大，表中STD均在0.186～0.270之间，因此STD指标在条件优化时作用不大。

3.6 最终模型

在最佳条件下，即光谱处理方法为一阶导数1D结合SNV，波长区间为R1（6931～6017 cm1），潜变量数为5，用建模集建立模型，对验证集进行预测。建模集和验证集样本的预测结果。各个样本基本上均匀分布在斜率为1的直线两侧，虽然误差并不太小，但可能是这套数据本身的特征，所建模型反映的正是数据的本质。 RMSEC和RMSEP为0.881%和1.167%，对应的相关系数分别为0.95和 0.88。

4 结 论

利用采样误差分布分析方法对二氯丙醇产品生产过程的氯化液进行近红外光谱分析，建立了杂质3-氯-1，2-丙二醇浓度的分析模型。SEPA方法通过对样本数据多次随机划分建立多子模型，进而计算中位数、标准偏差、偏斜度和分布峰度等统计指标，通过综合分析对近红外光谱分析模型进行条件优化、建模和模型评价等工作。结果表明，在光谱预处理方法选择和波长选择中，SEPA能选择更加合理的结果，波长选择的结果更具解释性。用独立验证集进行模型评价所建模型的稳健性。本研究为建立实用的近红外光谱分析模型提供了一种有效的方法。

Reference

1 Xu Q S， Liang Y Z. Chemom. Intell. Lab. Syst.， 2001， 56（1）： 1-11

2 Galvao R， Araujo M， Jose G， Pontes M， Silva E， Saldanha T. Talanta， 2005， 67（4）： 736-740

3 Deng B C， Yun Y H， Liang Y Z， Cao D S， Xu Q S， Yi L Z， Huang X. Anal. Chim. Acta， 2015， 880： 32-41

4 JIN Zhao-Xi， ZHANG Xiu-Juan， LUO Fu-Yi， AN Dong， ZHAO Sheng-Yi， RAN Hang， YAN Yan-Lu. Spectrosc. Spect. Anal.， 2016， 36（12）： 3920-3925

靳召晰， 张秀娟， 罗付义， 安 冬， 赵盛毅， 冉 航， 严衍禄. 光谱学与光谱分析， 2016， 36（12）： 3920-3925

5 YUN Yong-Huan， DENG Bai-Chuan， LIANG Yi-Zeng. Chinese J. Anal. Chem.， 2015， 43（11）： 1638-1647

云永欢， 邓百川， 梁逸曾. 分析化学， 2015， 43（11）： 1638-1647

6 Kennard R W， Stone L A. Technometrics， 1969， 11（1）： 137-148

7 Li H D， Zeng M M， Tan B B， Liang Y Z， Xu Q S， Cao D S. Metabolomics， 2010， 6（3）： 353-361

8 Chen W C， Du Y P， Zhang F Y， Zhang R Q， Ding B Y， Chen Z K， Xiong Q. J. Chemom.， 2017： e2933

9 Zhang F Y， Chen W C， Zhang R Q， Ding B Y， Yao H M， Ge J， Ju L， Yang W Y， Du Y P. Chemom. Intell. Lab. Syst.， 2017， 171： 234-240

10 Zhang R Q， Zhang F Y， Chen W C， Yao H M， Ge J， Wu S C， Wu T， Du Y P. Chemom. Intell. Lab. Syst.， 2018， 175： 47-54

11 GB/T 13097-2007， Determination of Phthalate Esters in Foods. National Standards of the People's Republic of China

工業用环氧氯丙烷. 中华人民共和国国家标准， GB/T 13097-2007

12 Burns D A.， Ciurczak E W. Handbook of Near-Infrared Analysis， 2008： CRC Press

13 Jiang J H， James B R， Siesler Heinz W， Yukihiro O. Anal. Chem.， 2002， 74（14）： 3555-3565

14 Cai W S， Li Y K， Shao X G. Chemom. Intell. Lab. Syst.， 2008， 90（2）： 188-194

15 Zheng K Y， Li Q Q， Wang J J， Geng J P， Cao P， Sui T， Wang X， Du Y P. Chemom. Intell. Lab. Syst.， 2012， 112（6）： 48-54

16 Yun Y H， Wang W T， Deng B C， Lai G B， Liu X B， Ren D B， Liang Y Z， Fan W， Xu Q S. Anal. Chim. Acta， 2015， 862： 14-23

17 Lin Y W， Xiao N， Wang L L， Li C Q， Xu Q S. Chemom. Intell. Lab. Syst.， 2017， 168： 62-71

18 Deng B C， Yun Y H， Liang Y Z， Yi L Z. Analyst， 2014， 139： 4836-4845

19 Deng B C， Yun Y H， Ma P， Lin C C， Ren D B， Liang Y Z. Analyst， 2015， 140（6）： 1876-1885

20 Wang L L， Lin Y W， Wang X F， Xiao N， Xu Y D， Li H D， XU Q S. Chemom. Intell. Lab. Syst.， 2018， 172： 229-240

21 Michael B， Hans P V. FT-NIR Atlas. 1993： VCH publisher. 76， 89