基于ATR-FTIR光谱测量结晶过程溶液浓度的变量稳定加权混合收缩方法
徐啟蕾, 郭鲁钰, 杜康, 单宝明, 张方坤*
青岛科技大学自动化与电子工程学院, 山东 青岛 266061
*通讯作者 e-mail: f.k.zhang@hotmail.com

作者简介: 徐啟蕾, 女, 1980年生, 青岛科技大学自动化与电子工程学院副教授 e-mail: 1255707511@qq.com

摘要

针对衰减全反射-傅里叶变换红外(ATR-FTIR)光谱仪用于测量结晶过程溶液浓度时, 因光谱谱线维度高、 无关变量多, 导致的标定模型预测精度低、 可解释性差等问题, 提出了一种变量稳定加权混合收缩的新方法。 首先提出对光谱变量进行随机二进采样, 将建立的优秀子模型中变量被选频率与所有子模型中变量回归系数的稳定性指标进行加权评价的稳定加权变量种群分析法(SWVCPA)。 通过对变量的重要性进行排序, 采用指数递减函数在迭代过程中逐渐强制滤除重要性低的变量, 实现了对光谱变量空间的初步收缩, 并大幅提高了收缩的稳定性。 然后在收缩后的子空间继续使用一种新的动态麻雀算法(DSSA), 以最小化训练预测均方根误差(RMSEC)为适应度函数进一步优化变量组合。 这种混合优化方式融合了两类变量选择算法的优点, 通过子模型竞争的方法确保了前期变量收缩的稳定性, 防止算法陷入局部最优; 通过智能优化算法避免了对剩余变量组合的遍历寻优, 允许保留更多的变量进行精准选择。 为了验证新方法的性能, 使用L-谷氨酸溶液冷却结晶过程中6种不同浓度下采集到的ATR-FTIR光谱数据进行测试。 结果表明, 新方法将光谱变量数从613个减少到46个, 与原始光谱相比, 使用选择后变量建立的偏最小二乘法(PLSR)模型其预测均方根误差(RMSEP)为从1.727 9降低到0.165 4, 预测决定系数(R2)从0.973 7提高到0.999 7。 另外相比于特征谱段、 遗传算法(GA)以及变量种群组合分析法(VCPA)选择变量建立的模型, 使用新方法建立的溶液浓度预测模型具有更高的准确性和稳定性, 说明该方法对提高使用ATR-FTIR光谱法测量冷却结晶过程溶液浓度准确性和可靠性具有一定的实际应用价值。

关键词: 变量选择; 溶液浓度测量; 结晶; ATR-FTIR光谱; 智能优化算法
中图分类号:O657.3 文献标志码:A
A Hybrid Shrinkage Strategy Based on Variable Stable Weighted for Solution Concentration Measurement in Crystallization Via ATR-FTIR Spectroscopy
XU Qi-lei, GUO Lu-yu, DU Kang, SHAN Bao-ming, ZHANG Fang-kun*
College of Automation and Electronic Engineering, Qingdao University of Science and Technology, Qingdao 266061, China
*Corresponding author
Abstract

In this paper, a new method for stable weighted mixture contraction of variables is proposed to address the problems of low prediction accuracy and poor interpretability of calibration models due to high spectral line dimensionality and many irrelevant variables when using attenuated total reflection-Fourier transform infrared (ATR-FTIR) spectrometers for measuring solution concentrations in crystallization processes. The method first proposes a stable weighted variable population analysis (SWVCPA) with a random binary sampling of the spectral variables and a weighted evaluation of the selected frequencies of the variables in the established superior sub-models and the stability indicators of the regression coefficients of the variables in all sub-models. By ranking the importance of variables and using an exponentially decreasing function to gradually force the filtering out of variables of low importance during the iterative process, an initial contraction of the spectral variable space is achieved, and the stability of the contraction is substantially improved. Then a new Dynamic Sparrow Search Algorithm (DSSA) is continued on the shrunken subspace to optimize the combination of variables further using the minimization of the root mean square error of training prediction (RMSEC) as the fitness function. This hybrid optimization approach combines the advantages of both types of variable selection algorithms, ensuring the stability of the prior variable contraction through a sub-model competition approach, preventing the algorithm from falling into a local optimum, and avoiding the traversal search for the remaining variable combinations through an intelligent optimization algorithm, allowing more variables to be retained for accurate selection. ATR-FTIR spectral data collected at six different concentrations during the cooling and crystallization of L-glutamic acid solutions were tested. The results showed that the new method reduced the number of spectral variables from 613 to 46 and that the root mean square error of prediction (RMSEP) was reduced from 1.727 9 to 0.165 4, and the coefficient of determination of prediction (R2) improved from 0.973 7 to 0.999 7 for the partial least squares (PLSR) model built using the selected variables compared to the original spectra. Genetic algorithm (GA) and variable population combination analysis (VCPA) for selecting variables, the solution concentration prediction model developed using the new method has higher accuracy and stability, indicating the practical application of the method to improve the accuracy and reliability of measuring solution concentration in cooling crystallization processes using ATR-FTIR spectroscopy.

Keyword: Variable Selection; Solution concentration monitoring; Crystallization; ATR-FTIR spectra; Intelligent optimization algorithm
引言

准确监测溶液浓度对优化结晶过程控制具有重要意义。 作为一种原位测量手段, 红外光谱具有快速、 简便、 无损和多组分检测等特点, 已广泛应用于石油、 食品、 制药等领域[1]。 近红外光谱包含从样品属性以及环境和仪器中获得的大量信息, 导致变量数远大于样本数, 这种高维数据带来了很多传统统计方法无法应对的“ 维数灾难” [2], 使得预测模型容易过拟合或不准确。 由于光谱变量间存在严重的共线性, 使得建立的预测模型不稳定、 可解释性差[3]。 为解决上述问题, 提出了变量选择技术。 变量选择是特征提取过程, 即从全近红外光谱中选择特定波长来建模。 近几十年来, 变量选择作为近红外光谱定量分析中的一个重要统计工具, 已经开发了大量新方法; 包括基于蒙特卡罗无信息变量消除(monte carlo-uninformative variable elimination, MC-UVE)方法[4]、 遗传算法(genetic lgorithm, GA)[5]、 移动窗口法(moving window, MW)[6]等。 近年来, 随着计算机算力的提升, 基于模型总体分析(model population analysis, MPA)[7]一类方法成为变量选择算法的研究重点。 Yun等[8]结合达尔文进化论所基于的简单但有效的“ 适者生存” 原则, 开发了一种用于选择多分量光谱数据关键波长的最佳组合的新策略, 称为竞争自适应重加权采样(competitive adaptive reweighted sampling, CARS)。 于雷等[9]基于迭代保留信息变量(iteratively retaining informative variables, IRIV)的方法, 对高光谱变量进行选择, 用以估测大豆叶片SPAD的值, 这是一种通过随机组合考虑变量之间可能交互作用的策略, 每一轮迭代中只保留强信息和弱信息变量, 直到不存在无信息和干扰变量。 Deng等[10]提出了自举软收缩(bootstrapping soft shrinkage, BOSS)方法, 结合了MPA、 软收缩、 加权抽样(weighted bootstrap sampling, WBS)策略, 通过将众多偏最小二乘(partial least squares, PLS)子模型的回归系数转化为可变权重, 并基于加权采样的步更新, 实现变量空间的软收缩。 Yun等[11]提出了变量组合种群分析(variable combination population analysis, VCPA)方法, 采用指数递减函数(exponentially decreasing function, EDF)剔除贡献不大的变量, 从而缩小变量空间。 龙燕等[12]使用变量迭代空间收缩法(variable iterative space shrinkage approach, VISSA), 提取了与干旱胁迫高度相关的叶绿素荧光参数, 准确检测苗期番茄干旱胁迫状态。 然而上述方法均使用单一指标如加权、 回归系数、 频数评价变量重要性, 其稳定性有待进一步提高[3]

针对上述问题, 提出一种基于被选频率与回归系数稳定性结合的变量收缩算法(stabilized weights variable combination population analysis, SWVCPA), 并在变量收缩到设定个数后使用改进动态麻雀算法(dynamic sparrow search algorithm, DSSA)进一步寻找最优变量组合, 既可以防止单纯依靠一种评价指标导致的模型不稳定问题, 又可以克服直接在高维光谱中使用智能优化算法易陷入局部极小的问题。 避免了VCPA对剩余变量所有组合的遍历寻优过程, 减少了建模次数并允许设置更大的输出变量空间来寻找更优的变量组合。 通过使用L-谷氨酸溶液冷却结晶过程中采集到的ATR-FTIR光谱数据进行建模测试, 与使用全谱段、 特征谱段、 GA以及原始VCPA选择变量建立的偏最小二乘模型进行比较。 结果表明, 本方法能够有效过滤光谱数据中的无关变量, 且保留变量的重要性和稳健性更高, 可以有效提高冷却结晶过程溶液浓度监测模型的准确度。

1 算法原理
1.1 SWVCPA原理

为了改变VCPA采用单一变量重要性评价规则带来的问题, 首先提出一种被选频率与回归系数稳定性结合的SWVCPA算法。 与VCPA类似, 其也是基于MPA思想的建模方法。 实现步骤由二进制矩阵抽样(binary matrix sampling, BMS)、 变量重要性评价和指数递减函数(exponentially decreasing function, EDF)强制去除冗余变量构成。

首先生成k个随机二进制矩阵, 矩阵中“ 1” 表示该位置处的变量被选择用于建立PLSR模型, 0” 表示变量未被选择。 通过计算所有子集模型的交叉验证均方根误差(root mean squared error of cross-validation, RMSECV), 选择一定数量的较优模型, 统计其中变量被选的频数; 同时引入回归系数稳定性指标与频数加权综合评价变量的重要性。 回归系数稳定性指标如式(1)所示

cj=βj¯s(βj)(1)

式(1)中, cj表示第j个变量的回归系数稳定性指标。 βj¯表示第j个变量M次二进制采样中的回归系数的平均值, s(β j)表示第j个变量在M次采样中的标准差。 cj反映了第j个变量在模型中的重要性与稳定性。

将变量频率与回归系数稳定性进行加权, 作为SWVCPA的变量重要性评价指标, 如式(2)所示

wj=a×fjbest+(1-a)×cj(2)

式(2)中, fjbest为第j个变量在较优模型中的被选频率, cj为第j个变量在M次采样中的回归系数稳定性指标, a为加权系数, 范围为0到1, 默认可以设置为0.5。

相比VCPA只考虑较优模型中变量被选频率, SWVCPA采用了较优模型中变量频率和所有模型偏最小二乘回归系数稳定性加权的方式判断变量重要性。 这种方式使得评价指标受二进矩阵采样中的随机性影响更小, 稳定性更高, 并允许设置更小的BMS采样次数以加快算法运算速度。 同时在迭代后期对于无关变量的滤除速度更快, 最终选择的变量更加稳定。

使用如式(3)的指数衰减函数计算当前迭代次数i下应该保留的变量比率ri, 按较优模型中重要性排序高低进行保留。

ri=e-θi(3)

通过θ 控制EDF曲线, 表示为式(4)

θ=ln(P/L)K(4)

式(4)中, P为光谱变量数, K为EDF的最大迭代次数, L为迭代完成后保留的变量个数。

当SWVCPA完成对变量空间的初步压缩后, 使用改进的动态麻雀算法对剩余变量进一步优化。

1.2 动态麻雀优化算法

受麻雀在自然中觅食策略的启发, Yan等[13]提出了麻雀搜索算法(sparrow search algorithm, SSA)。 麻雀算法是一种种群优化算法, 其将种群划分为发现者和加入者, 并从中随机挑选麻雀作为警戒者。 其中发现者按一定规则更新位置, 负责为种群寻找适应度更好的区域。 加入者负责监视发现者, 并向最优的发现者移动。 警戒者负责随机选择麻雀加速收敛到最优位置或让最优位置的麻雀重新跳出该位置。 麻雀算法的具体信息可以参见文献[13]。 与其他智能优化算法相比, SSA算法具有稳定性较高、 收敛速度快, 更适用于光谱变量选择算法的需求。 但是SSA算法收敛速度较快导致容易陷入局部最优解。

考虑在原始麻雀优化算法中, 发现者和警戒者的比例是固定的, 如果可以根据迭代次数将种群中不同角色的数量进行动态变化, 将会获得更好的效果。 因此提出一种改进动态麻雀搜索算法, 其发现者和警戒者的比例按式(5)更新。

PD=0.9-exp-imaxiSD=0.1+0.3exp-imaxi(5)

式(5)中i为迭代次数, imax为最大迭代次数。 可以看出, 在优化前期发现者的比例(PD)高, 便于获得更广泛的初始值, 避免模型过早收敛。 随着迭代次数的增加, 发现者比例降低, 加入者比例相应增加, 并在最优值附近进行更精确的搜索。 这种更强的局部搜索能力非常适合光谱邻近变量间共线性强的特点, 便于找到更优的特征变量。 另外, 警戒者的比例(SD)也随着迭代次数的增加而增加, 不仅增加了模型的收敛速度, 也防止种群聚集在极小值, 使模型后期依然具备一定的搜索能力, 提高了模型跳出局部极小值的能力。

结合上述两种改进算法, 提出一种新的混合收缩算法SWVCPA-DSSA, 其流程如图1所示, 由SWVCPA和DSSA两部分组成。 通过输入原始ATR-FTIR光谱数据, 并初始化基本参数, 最终输出经过选择的最优变量子集, 用于建立监测模型。

图1 SWVCPA-DSSA算法流程图Fig.1 Flow chart of SWVCPA-DSSA algorithm

2 实验部分
2.1 数据采集

在实验室设备中进行冷却结晶实验并获取光谱数据。 首先使用精确度为万分之一的高分辨率分析天平(梅特勒-托利多)称量9 g· L-1的谷氨酸(C5H9NO4)样品, 样本由Sigma公司制造, 纯度为99%。 将样本置于2L夹套结晶器中, 并向其中加入1 L的蒸馏水, 配置成初始的谷氨酸样品溶液。 将溶液加热至75 ℃, 并保持90 min, 整个实验过程中一直匀速搅拌, 以确保谷氨酸完全溶解。 谷氨酸溶液吸光度光谱由用金刚石制成的ATR浸没探针作为内部反射原件获得, 探针连接到FTIR分光光度计, 再将数据传输到一台运行着ReactIR15软件(梅特勒-托利多)的监控计算机上记录。 冷却过程保持溶液温度以0.8 ℃· min-1的快速冷却速率将溶液冷却至15 ℃。 并在冷却过程中每分钟采样一次ATR-FTIR光谱和溶液温度数据, 直到图像监测系统检测到溶液内结晶核的产生。 检测到结晶核出现后, 停止数据采样并将溶液温度重新升高至75 ℃。 再加入称量好的谷氨酸晶体配置成更高浓度的溶液样本, 保持90 min确保谷氨酸完全溶解再开始新的冷却过程。 最终获得的光谱数据如图2所示, 采集光谱范围为3 000~650 cm-1

图2 L-谷氨酸水溶液红外光谱图Fig.2 Infrared spectrum of L-glutamic acid solution

2.2 数据预处理

光谱预处理方法可以分为基线校正、 散射校正、 平滑处理和尺度缩放四类[14]。 为了避免不同预处理对变量选择算法结果可能存在的干扰, 也为了验证在不进行特殊预处理情况下不同变量选择方法的预测性能, 数据均使用标准化预处理。

2.3 模型建立与评价

实验共获得6种(9, 15, 21, 27, 33, 39 g· L-1)不同浓度样品的光谱数据, 从中等间隔抽取240组数据, 并随机将其中160组作为训练集, 用来训练溶液浓度标定模型。 其余80组作为测试集, 用来测试模型的预测效果。 分别使用决定系数(R2)、 建模均方根误差(root mean square correction error, RMSEC)和预测均方根误差(root mean square prediction error, RMSEP)作为模型的评价指标。

3 结果与讨论

几种方法运行50次平均选择的变量如图3所示。 由文献[15]可知, L-谷氨酸溶液的特征吸收峰主要集中在1 800~1 200 cm-1波长范围内, 单独选择该特征谱段的变量用PLS建模。 GA-PLS选择了几个主要特征区域, 缩减了变量空间, 但可以看出其选择的变量仍然较为分散, 且包含了许多冗余变量。 VCPA和SWVCPA选择的主要变量基本相同, 主要集中在1 450~1 400 cm-1波段, 这是CH2分子链形变和对称羧酸离子伸缩振动峰的特征区域。 另外还选择了1 200 cm-1波段的变量, 其代表C— O羧基伸缩的特征。 与VPCA相比, SWVPCA选择了更少的无关谱段(如2 300和850 cm-1), 表明通过引入回归系数稳定性指标使得SWVCPA的选择更加稳定, 更加集中于重要性高的特征变量。

图3 不同算法在L-谷氨酸ATR-FTIR光谱数据中选择的变量
(a): 原始光谱; (b): 不同算法在样品选择的变量
Fig.3 Variables selected by different algorithms on L-glutamate ATR-FTIR spectral data
(a): Original spectra; (b): Variables selected by different algorithms on L-glutamat ATR-FTIR spectral date

这一优势还可以从另一方面得到验证, 当BMS设置为1 000时, VCPA可优选出关键变量, 当减少BMS次数为100时VCPA无法选择出合适的变量, 如图4(a, b)所示。 而BMS设置为100时, SWVCPA与BMS为1 000时选择的变量基本相同, 如图4(c)所示。 当减少随机二进制矩阵抽样(BMS)次数为100时, VCPA无法选择出合适的变量, 而SWVCPA与采样1 000次的VCPA选择基本相同。 这种良好的稳定性使得SWVCPA可以在保证效果的情况下使用更少BMS来大幅提高程序的运行速度。

图4 (a) BMS设置为1 000的VCPA变量选择过程; (b) BMS设置为100时VCPA变量选择过程; (c) BMS设置为100时SWVCPA变量选择过程Fig.4 (a) VCPA variable selection process with BMS set to 1 000; (b) VCPA variable selection process with BMS set to 100; (c) SWVCPA variable selection process with BMS set to 100

另外由如图4中可以看出, 一些较弱的特征信息在最后几次迭代过程中被滤除。 比如1 730 cm-1附近的变量, 其代表C=O羧酸盐振动, 在酸性条件下的特征表现较弱。 然而由于计算机内存的限制, VCPA允许保留的最大变量个数为14[16], 而采用混合策略的SWVPCA-DSSA则可以保留更多的变量。 SWVCPA保留不同变量个数对最终模型的预测效果如图5所示。 可以看出, 模型的RMSEC呈现出先平稳后大幅增加的趋势, 说明如果最终保留变量数过少可能会导致一些特征变量的丢失, 使模型预测性能下降; 另外, 变量从100到45变化的平稳趋势, 证明了DSSA可以在缩减后的空间中得到较为稳定的模型, 体现了采用混合策略的必要性和优势。 与GA选择相比, 改进麻雀算法所选择变量建立模型准确度更高、 标准差更小、 模型更稳定。

图5 SWVCPA保留变量数对最终模型的影响Fig.5 Effect of the number of SWVCPA retention variables on the final model

综上所述, 可以得到不同选择策略在变量空间下的表现, 如图6所示。 假设每次滤除的变量都是使模型建模效果提升最大的变量, 画出图6中绿色曲线所示的最优曲线图。 其在达到最小值点后, 继续删除变量将导致模型效果变差。 红色区域与虚线为遗传算法的理论表现, 由于直接在全谱中使用GA其极易陷入局部极小, 选择的变量不稳定使得模型效果不可信。 VPCA的选择过程如黄色曲线所示, 其在达到最大允许保留变量前, 使用EDF进行收缩, 通过对最后保留的L个变量的所有组合进行建模找到最好的建模序列。 这种方式不仅效率低下, 且使得参数L不能设置的过大, 如果最优模型的变量数Lbest大于Lmax, 则无法得到最优模型。 混合搜索策略的优化过程如蓝线所示, 其吸收了两种算法的优点, 通过引入智能优化算法使得SWVCPA可以保留更多的变量Lh, 并在接近最优变量数时使用DSSA进行剩余优化, 使结果尽可能接近最优模型。

图6 不同选择策略在变量空间下的表现示意图Fig.6 Schematic representation of different selection strategies in the variable space

通过将几种算法运行50次, 采用5折交叉验证训练模型, 比较表1中模型的结果验证了上述分析的正确性。 可以看出由于引入过多无关变量, 导致模型易受到数据中噪声的干扰, 使用全谱段建立的PLS模型预测效果最差。 通过选择特征谱段减少建模变量后建立的PLS模型的预测性能大幅提高, 证明了变量选择对于提高模型性能的必要性。 使用GA-PLS选择的变量数进一步减小, 但由于在高维变量中直接使用智能优化算法, 导致选择的变量数不稳定, 且使用了更多的因子数来建立PLS模型, 有陷入过拟合的风险。 VCPA选择了更少的变量, 模型的预测效果进一步提高, 但受到遍历组合所需时间的限制, 使得保留的变量数不能过多, 可能会丢失一部分有用信息。 相比之下, SWVCPA-DSSA结合了对比算法的优点, 允许设置更灵活的保留变量数, 并更稳定地选择特征变量, 因此建立了更优的模型。 结果显示, 使用SWVCPA-DSSA选择变量建立的PLS模型的预测均方根误差最小为0.165 4, 决定系数R2为0.999 7, 模型的预测效果最好, 且稳定性更高。

表1 不同方法结果对比表 Table 1 Comparison table of the results of different methods
4 结论

为了提高基于ATR-FTIR光谱建立的冷却结晶过程溶液浓度监测模型的准确性, 提出了一种SWVCPA-DSSA算法, 结合回归系数稳定性与被选频数进行变量选择, 并在缩减后的变量空间上使用动态麻雀算法进行进一步组合寻优。 以L-谷氨酸溶液冷却结晶过程为例, 通过与其他方法进行对比, 证明了新方法能够更加稳定、 准确地提取红外光谱中的有效信息, 提高预测模型的准确性, 具有一定的应用价值。

参考文献
[1] Gao Z G, Rohani S, Gong J B, et al. Engineering, 2017, 3(3): 343. [本文引用:1]
[2] Yun Y H, Li H D, Deng B C, et al. TrAC Trends in Analytical Chemistry, 2019, 113: 102. [本文引用:1]
[3] Yan H, Song X Z, Tian K D, et al. Spectrochimica Acta Part A: Molecular and Biomololecular Spectroscopy, 2019, 210: 362. [本文引用:2]
[4] Cai W S, Li Y K, Shao X G. Chemometrics and Intelligent Laboratory Systems, 2018, 90(2): 188. [本文引用:1]
[5] Yang Z F, Xiao H, Zhang L, et al. Spectrochimica Acta Part A: Molecular and Biomololecular Spectroscopy, 2019, 223: 117327. [本文引用:1]
[6] Li Q, Huang Y, Song X Z, et al. Spectrochimica Acta Part A: Molecular and Biomololecular Spectroscopy, 2019, 214: 129. [本文引用:1]
[7] Deng B C, Lu H M, Tan C Q, et al. Chemometrics and Intelligent Laboratory Systems, 2018, 172: 223. [本文引用:1]
[8] Yun Y H, Zhang J H, Chen H M, et al. Chemometrics and Intelligent Laboratory Systems, 2020, 197: 103920. [本文引用:1]
[9] YU Lei, ZHANG Tao, ZHU Ya-xing, et al(于雷, 章涛, 朱亚星, ). Transactions of the Chinese Society of Agricultural Engineering(农业工程学报), 2018, 34(16): 148. [本文引用:1]
[10] Deng B C, Yun Y H, Cao D S, et al. Analytica Chimica Acta, 2016, 908: 63. [本文引用:1]
[11] Yun Y H, Bin J, Liu D L, et al. Analytica Chimica Acta, 2019, 1058: 58. [本文引用:1]
[12] LONG Yan, MA Min-juan, WANG Ying-yun, et al(龙燕, 马敏娟, 王英允, ). Transactions of the Chinese Society of Agricultural Engineering(农业工程学报), 2021, 37(11): 172. [本文引用:1]
[13] Yan H, Song X Z, Tian K D, et al. Spectrochimica Acta Part A: Molecular and Biomololecular Spectroscopy, 2019, 210: 362. [本文引用:2]
[14] Xue J K, Shen B. Systems Science & Control Engineering, 2020, 8(1): 22. [本文引用:1]
[15] Jiao Y P, Li Z C, Chen X S, et al. Journal of Chemometrics, Vol. , 2020, 34(11): 3306. [本文引用:1]
[16] Zhang F K, Liu T, Wang X Z, et al. Journal of Crystal Growth, 2017, 459: 50. [本文引用:]