聚谷氨酸发酵过程中ATR-FTIR光谱信号的分数阶基线校正
何年, 单鹏*, 贺忠海, 王巧云, 李志刚, 吴缀
东北大学秦皇岛分校控制工程学院, 河北 秦皇岛 066000
*通讯作者 e-mail: peng.shan@neuq.edu.cn

作者简介: 单 鹏, 1985年生,东北大学秦皇岛分校控制工程学院讲师 e-mail: 2031816313@qq.com

摘要

采用衰减全反射傅里叶变换红外光谱法(ATR-FTIR), 结合多元校正模型对γ-聚谷氨酸(γ-PGA)发酵过程中两种主要底物葡萄糖和谷氨酸钠的浓度进行间接测量, 为优化发酵系统控制提供重要的反馈信息。 光谱测量中经常出现的基线漂移会严重影响后续多元校正模型的性能, 需要采用基线校正算法对光谱进行预处理。 现有流行的基线校正算法多数是基于Whittaker Smoother(WS)平滑算法, 这些算法均采用整数阶微分对拟合基线进行约束, 表达能力有限。 针对现有基线校正算法中的整数阶微分自适应性差的问题, 利用更加灵活的分数阶微分对基线进行约束, 提出了一种基于分数阶的基线校正算法, 实现对整数阶基线校正的扩展。 总共进行了5个批次的γ-PGA发酵实验, 并对不同批次和全部批次的ATR-FTIR光谱数据分别进行了分数阶基线校正, 模型的预测精度均得到不同程度的提升。 实验结果表明, 只有在批次2时, 基于整数阶的基线校正效果最好; 其他批次的基线校正效果最好时的阶次均为分数阶。 这也表明了分数阶微分(包含整数阶微分)对基线的约束更加合理。 同时发现全部批次的整体基线校正效果远远差于单一批次的效果, 原因可能是各批次发酵光谱的基线是不同的, 对不同的批次需要选用不同的阶次以获得最佳的基线校正。 此外, γ-PGA发酵样品的ATR-FTIR光谱测量是以蒸馏水为背景, 会在3 100~3 600 cm-1波数范围内出现负水峰, 形成有害的干扰信号; 分数阶基线校正后的光谱表明, 分数阶基线校正算法将负的水峰当作基线, 在一定程度上进行了消除。 综上分析, 分数阶基线校正算法不仅扩展了传统整数阶基线校正算法的应用范围, 也为消除ATR光谱中负的水峰提供了新的解决思路。

关键词: 分数阶微分; 光谱预处理; 基线校正
中图分类号:O657.33 文献标志码:A
Study on the Fractional Baseline Correction Method of ATR-FTIR Spectral Signal in the Fermentation Process of Sodium Glutamate
HE Nian, SHAN Peng*, HE Zhong-hai, WANG Qiao-yun, LI Zhi-gang, WU Zhui
School of Control Engineering, Northeastern University at Qinhuangdao, Qinhuangdao 066000, China
*Corresponding author
Abstract

In this paper, Attenuated Total Reflection Fourier Transformed Infrared Spectroscopy (ATR-FTIR) combined with the multivariate calibration model was used to realize the indirect measurement of the concentration of two main substrates (glucose and sodium glutamate) during the fermentation process of γ-polyglutamic acid (γ-PGA), which could provide feedback information for the fermentation process. The frequent baseline drift phenomenon in the spectrum measurement will seriously affect the performance of the subsequent multivariate calibration model, and it is necessary to use the baseline calibration algorithm to preprocess the spectrum. Most of the popular baseline correction algorithms are based on the Whittaker Smoother (WS) smoothing algorithm. And use integer-order differentials with limited expressive power to constrain the fitted baseline. Because of the poor adaptability of integer-order differential in the existing baseline correction algorithms, we use more flexible fractional-order differentials to constrain the baseline and then propose a baseline correction algorithm based on fractional-order, which realizes the extension of the integral order baseline correction. 5 batches of γ-PGA fermentation experiments were carried out, and the ATR-FTIR spectra of different batches and all batches were subjected to fractional baseline correction respectively; subsequently, the prediction accuracy of each model was improved to some extent. The experimental results show that only in batch 2 the baseline correction effect based on the integer-order is the best; the orders to obtain the best baseline correction effect for other batches were all fractional-order. Italso reflects that the constraint of the fractional-order derivative (including the integer-order derivative) on the baseline is reasonable. At the same time, it is found that the overall baseline correction effect of all batches is far worse than that of a single batch. The reason may be that the baseline of the spectra for each fermentation batch is different. Different orders need to be selected for different batches to achieve the best effect of baseline correction. In addition, the background spectrum was acquired with distilled water as the reference before measuring each γ-PGA fermentation sample. Anegative water peak thus inevitably appears in the wavenumber range of 3 100~3 600 cm-1 and forms harmful interference signals; the fractional baseline-corrected spectra show that the fractional-order baseline correction algorithm regards the negative water peak as the baseline and eliminates it to a certain extent. In summary, the fractional-order baseline correction algorithm expands the application range of the traditional integer-order baseline correction algorithm and provides a new solution to eliminate negative water peaks in the ATR spectra with water as the background spectrum.

Keyword: Fractional differentiation; Spectral preprocessing; Baseline correction
引言

发酵过程中对物料浓度的在线检测是保证产品质量的关键, 目前仍有大部分的工厂采用传统的分离分析方法, 费时繁琐, 无法实现发酵的实时控制。 而衰减全反射傅里叶变换红外光谱(attenuation reflectance Fourier transformation infrared spectrometry, ATR-FTIR)技术具有快速、 绿色、 无损等优点, 应用广泛[1, 2, 3], 十分适合发酵过程中参数的在线检测。 本文通过对聚谷氨酸发酵过程中采集的ATR-FTIR光谱与建立的多元校正模型间接的测量发酵过程中主要底物葡萄糖和谷氨酸钠的浓度。 因光谱中存在基线漂移, 影响定量分析精度, 所以在建立多元校正模型前需进行光谱预处理, 去除噪声和基线。

其中基于WS平滑算法的基线校正算法应用广泛, 已成功应用于各种光谱分析中。 Eilers等首先将WS算法应用于信号的平滑和插值[4], 随后通过添加非对称权重, 在峰值信号处实现较强的光滑, 将光滑曲线作为拟合基线, 提出非对称最小二乘基线校正算法[5](asymmetric least squares, AsLS); Zhang等通过改进非对称权重, 实现自适应迭代, 提出自适应迭代加权惩罚最小二乘[6](adaptive iteratively reweighted penalized least squares, airpls); Baek[7]等将权重按logistic函数分配, 实现权重信息柔化。 (asymmetrically reweighted penalized least squares smoothing, arpls); 姜安[8]等通过引入拟合残差的一阶导数作为新的惩罚项, 用直方图估计背景设定阈值, 加快迭代速度(jiang improved asymmetric least squares baseline correction algorithm, jasls); He[9]等通过改进jasls的阈值, 形成(improved asymmetric least squares method, iasls); 近年Xu[10]等给惩罚项加上新的权重, 形成重加权惩罚基线校正(doubly reweighted penalized least squares, drpls); Ye[11]等通过改进arpls的权重赋值函数, 实现进一步优化(improved asymmetrically reweighted penalized least squares, iarpls)。

这些算法多是以改变权重或自动迭代权重的方式进行改进, 但对基线的约束都采取相同策略, 选择阶次较低的整数阶微分; 此外还引入新的惩罚项来改进算法性能。 这种改进可以通过引入不同的惩罚项, 实现对基线不同的约束, 原模型中的惩罚项是对基线的粗糙度进行约束, 在jpls算法中加入了真实光谱和拟合基线残差的一阶微分作为新的惩罚项, 他们认为拟合出的基线不仅与原始数据之间的误差很小, 而且还要求它们的一阶导数很接近。 这使得对拟合基线的约束性更强, 能够适应不同的光谱。

上述提及的各种改进算法都未曾涉及对粗糙度描述方法进行改进, 都沿用整数阶微分或直接固定阶次, 使得对基线的约束不够灵活。 低阶整数阶微分通常只取1, 2和3阶, 可选择性极差; 同时考虑到整数阶微分不能很好的描述基线的特点, 且在实际信号中整数阶阶次的信号很少见。 故引入分数阶微分的概念, 提高算法灵活性, 扩展对粗糙度的描述, 从而进一步研究微分阶次对基线校正效果的影响。 我们提出的分数阶基线校正算法涵盖了原来的整数阶算法, 理论上认为分数阶基线校正效果不会差于整数阶基线校正效果; 这一推断在γ -PGA发酵光谱分析中进行了仔细的检验。

1 原理
1.1 AsLS基线校正算法原理

AsLS基线校正是在WS平滑算法的基础上改进而来, 首先对WS平滑算法进行简单介绍[4]。 WS优化目标如式(1)

argminz[i(yi-zi)2+λi(Δzi)2](1)

式(1)中, yi为原始信号的第i个点, zi为平滑序列的第i个点, λ为正则化系数, Δ为微分算子, 一阶微分可表示Δzizizi-1。 Eilers在文中将最小化目标函数的第一项称为保真度, 第二项称为粗糙度。 通过调节λ得到合适的平滑序列。 当λ越大, 对粗糙度的惩罚越大, 就要求序列越平滑。 对于一阶微分, 当λ越大时, 平滑序列z就越趋向于一条水平的直线。 同理二阶微分, 在λ越大时, 滑序列z就越趋向于一条倾斜的直线。

为简化代数运算用矩阵表示目标函数

Q=y-z|2+λDdz|2(2)

最小化式(2)可得

z=(E+λDTdDd)-1y

式(2)中, y为原始信号列向量, z为平滑序列列向量, E为单位矩阵, Dd为微分矩阵, 只能表示整数阶。 当信号长度为5, 取整数一阶、 二阶微分时, 表示如下

D1=-110000-110000-110000-11D2=1-210001-210001-21

AsLS在WS算法基础上引入非对称权重作用于保真度[5]

Q=W|y-z|2+λDdz|2(3)

最小化式(3)可得

z=(W+λDTdDd)-1Wy

式(3)中, W为权重对角矩阵, 式中权重系数Wi根据非对称的方式选择

Wi=pyi> zi1-pyizi

一般p取值范围为0.001~0.1。 λ取值范围为102~109, 固定迭代次数, W一般迭代10次。

AsLS也可以看作是平滑算法, 通过对含有峰的光谱信号进行平滑, 得到一条光滑的不含有峰的曲线作为基线, 这种基线校正方法不需要任何的先验信息, 只需要通过调节反对称权重和正则化系数就能得到一条适合的基线。

1.2 基于分数阶的基线校正算法原理

AsLS算法中的Dd微分算子只适用于整数阶, 在Dd的基础上扩展到分数阶, 较为简便的实现分数阶基线校正。 分数阶微分定义有不同的形式, 为了更好的包含原有整数阶, 选用Grumwald-Letnikov(GL)分数阶微分定义, 表示如式(4)[12, 13, 14]

dαf(x)=limh01hαm=0t-ah(-1)mΓ(α+1)m!Γ(α-m+1)f(x-mh)(4)

式(4)中: α 为阶数; h为微分步长; ta分别为微分的上、 下限; Γ (x)为Gamma函数。 当函数f(x)定义域为x∈ [a, t]且h=1时, 由式(4)可得出f(x)的分数阶微分表达式

dαf(x)dxαf(x)+(-α)f(x-1)+(-α)(-α+1)2f(x-2)++Γ(-α+n)n!Γ(-α)f(x-n)(5)

式(5)中: α 为阶数, f(x)的0阶微分为f(x)本身。 同样的将微分差值运算构造成矩阵的形式, 记为Dα

Dα=d1000d2d10d3d2d100dkd2d10dn+1dkd3d2d1

其中dn+1=Γ(-α+n)n!Γ(-α), n=t-a; 通常真实光谱波数数目在1 000以上, 但当n较大时,dn+1很小, 基本可以忽略不计, 同时为了加快计算速度, 节省存储空间, 设置一个长度k, 将k以后的系数全部置为零。 矩阵结构如下

Dα=d1000d2d10dkd2d1000dkd2d1000dkd2d1

实验表明, 当k较小时, 基线校正效果较差, 一般取20以上的值。 但计算速度会随k的增大而变慢。 用Dα 替换Dd即可实现分数阶基线校正(fractional differential asymmetric least squares, FdAsLS)。 当α 取整数阶时, 比原来的整数阶微分矩阵多了几项, 这对于一些信号的校正是不利的, 因为起始部分点的微分变化较大, 容易造成基线的突变, 但对于原始信号起始部分基本为零的光谱信号, 并不会产生影响。 在取分数阶时, 该现象尤为明显, 因为前k行的微分表达式都是不相同的, 从矩阵的前k行, 可以明显的观察到。 GL定义下的整数阶矩阵相比于原来整数阶矩阵, 略有不同。 当信号长度为5, GL定义下的整数一阶、 二阶微分时, 表示如下

D1α=10000-110000-110000-110000-11D2α=10000-210001-210001-210001-21

2 实验部分
2.1 数据采集

对于γ -PGA发酵实验, 选用的菌种为枯草芽孢杆菌亚种, 从中国工业微生物菌种保藏管理中心(China Center of Industrial Culture Collection, CICC)购买, 菌种编号为20643。 将以冻干粉的形式存储的菌种先进行活化培养, 然后在培养好的固体菌落中, 用接种环挑选一株生长状态良好的菌体, 接种于种子培养基(500 mL三角瓶装液量100 mL), 然后在37 ℃和180 r· min-1的恒温振荡培养箱中(THZ-92A, 跃进医疗器械有限公司, 中国上海)中培养10~16 h。 所用种子培养基为: 葡萄糖(10 g· L-1), 蛋白胨(10 g· L-1), 牛肉膏(5 g· L-1), 氯化钠(5 g· L-1)。 发酵培养基由葡萄糖(40 g· L-1), 酵母膏(5 g· L-1), 谷氨酸钠(35 g· L-1), 氯化铵(2 g· L-1), 磷酸氢二钾(5 g· L-1)和硫酸镁组成(0.5 g· L-1)组成。 种子培养基和发酵培养基均在121 ℃下灭菌20 min。 将经过种子培养的菌株接种到接种量为2%的发酵培养基中, 并将3 L的发酵培养液放入工作容积为5 L的发酵罐(GRJB-5D, 绿色生物工程有限公司, 中国镇江)中, 在37 ℃恒温和300 r· min-1搅拌速度的条件下进行发酵。

用配备有水平铂金钻石ATR采样附件(ZnSe, 单反射)的布鲁克Alpha型傅里叶变换红外光谱仪(德国, 埃特林根)上收集光谱数据。 在35 ℃下, 以8 cm-1的分辨率在4 000~600 cm-1的波数范围内进行64次扫描。 每个样品测量之前, 用蒸馏水作为参考获取背景光谱。 对于每个样品, 重复测量两次。 所得的平均光谱用于进一步分析。 作为主要底物, 葡萄糖(g· L-1)和谷氨酸钠(g· L-1)是用于监测γ -PGA发酵的参数。 总共进行了48 h的5次发酵实验, 获得151个样品的光谱及其发酵参数的标准值。 各批次样本数分别为14, 27, 40, 40和30, 它们的光谱如图1所示。

图1 各批次光谱Fig.1 Spectra of each batch

2.2 模型建立

各批次样本根据Kennard-Stone(KS)算法, 按3:1比例划分样本, 3份作为校正集, 建立模型。 其余1份构成测试集, 验证模型。 多元校正模型选择偏最小二乘回归(partial least squares regression, PLS)[15], PLS采用5折交叉验证, 从1~15中选出最佳潜变量个数。 评价指标选择PLS模型的校正集均方根误差(root mean square error of calibration, RMSEC)和测试集均方根误差(root mean square error of prediction, RMSEP)。 对5个批次数据分别建立原始光谱、 原始光谱+FdAsLS预处理、 原始光谱+AsLS预处理的PLS模型。 最后合并所有批次样本, 重新划分校正集和测试集, 建模过程不变。

3 结果与讨论

为验证分数阶基线校正效果, 进行六组实验。 每组实验以PLS模型的RMSEP为评价指标。 固定反对称权重p=0.001和分数阶微分长度k=20。 微分阶次从0.5到4.5, 间隔0.1, 共41个阶次; λ从100到109, 幂指数间隔0.5, 共19个数, 通过网格搜索法筛选出最佳参数组合。 为了比较分数阶和整数阶基线校正算法性能的优劣, 我们将所有整数阶参数组合(其中0阶是用未经过基线校正的原始数据建模)和最佳分数阶参数组合(对应最小的RMSEP)在六组实验上的结果进行了汇总(见表1)。 考虑到两种整数阶微分表示不同, 所有实验, 在取整数阶时, 选择AsLS中的原始定义。 其中只有批次2的预测模型误差减小最多所对应的基线校正微分阶次为整数阶, 其余模型最佳阶次都为分数阶。 可以反映出分数阶微分基线校正有着不低于整数阶的校正效果, 大多数情况下都超过了整数阶。 同时分数阶微分阶次取值任意, 运用灵活。 通过网格搜索法得到的最佳阶次中, 取值基本没有重复, 足以说明分数阶微分具有自适应性, 能够灵活地提取真实基线, 因此相比于从前单一的整数阶微分, 分数阶微分能够充分发挥基线校正的优势, 能最大限度的提高多元校正模型的预测精度。 在批次1中, 底物的预测误差大幅减小, 葡萄糖预测模型的RMSEP从2.098降到0.857, 预测精度提高较多。 同时其他批次预测精度都大幅提升, 表明了基线校正有助于后续光谱定量分析。 将5个批次的样本合并后建立模型(未进行预处理的原始数据), 可以明显地观察到, 全局模型的预测结果不如各批次单独建模(批次3除外)的局部模型; 且各批次合并后的样本即使经过基线校正后, 模型精度提升也远差于各批次单独进行基线校正后的预测精度, 即校正效果不明显。 这一现象的产生可以归因于不同发酵批次的基线是不相同的, 通过固定一组参数对所有批次光谱进行相同的基线校正是不合理的。 从实际情况来考虑, 在不同的发酵批次过程中, 很多因素(如: 谷草芽孢杆菌的活性、 测量条件、 仪器性能的差异等)都会造成基线的变化, 因此应该对不同批次光谱单独进行基线校正, 才能充分发挥基线校正的效果。

表1 不同阶次的基线校正效果比较 Table 1 Comparison of baseline correction effects of different orders

以批次3为例, 对原始光谱建立PLS模型, 发现预测结果很差, 但经过基线校正以后, 预测均方根误差减小, 图2可以直观的反映出这种变化, 其余各批次经过基线校正以后, 结果如图3和图4所示, 其中图3是预测葡萄糖浓度RMSEP最小时所对应的基线校正处理后的结果, 图4对应于谷氨酸钠。

图2 偏最小二乘回归分析
(a): 原始光谱葡萄糖预测结果; (b): 基线校正后葡萄糖预测结果; (c): 原始光谱谷氨酸钠预测结果; (d): 基线校正后谷氨酸钠预测结果
Fig.2 Partial least squares regression analysis
(a): Predicted result of glucose concentration by the original spectrum; (b): Predicted result of glucose concentration after baseline correction; (c): Predicted result of sodium glutamate concentration by the original spectrum; (d): Predicted result of sodium glutamate conentration after baseline correction

图3 葡萄糖RMSEP最小时对应各批次基线校正后的光谱Fig.3 Baseline corrected spectra of glucose for each batch with the minimum RMSEP

图4 谷氨酸钠RMSEP最小时对应各批次基线校正后的光谱Fig.4 Baseline corrected spectra of sodium glutamate for each batch with the minimum RMSEP

经过基线校正后的光谱相较于原始光谱, 除批次4用以预测谷氨酸钠浓度的校正光谱还保留负的水峰, 其余校正光谱基本不含或含有少量负的水峰。 因此认为峰信号全部是有用的, 通常AsLS被用于只含有全部为正峰或全部为负峰的信号的校正, 但是对于用水溶液测得的ATR光谱, 负的水峰对于后续的分析是无用的甚至是有害的, 所以可以通过AsLS基线校正将水峰扣除。 改进后的FdAsLS同样具有该效果, 且对负峰的扣除效果更好。

各批次校正后的光谱各不相同, 间接反映出发酵过程的特殊性。 且校正后的光谱各批次存在较大的差异, 如批次4预测谷氨酸钠时, 基线校正效果作用不大, 校正后的光谱基本不变, 且经过基线校正后, RMSEC反而变差了。 理论上RMSEC应该小于RMSEP, 但是在该数据中, 各批次经过基线校正后部分批次出现RMSEP小于RMSEC的情况。 这可能是样本太小所造成的, 同时测试集样本中个别样本预测很差, 但经过基线校正后, 预测精度立刻提升, 图2(a, b)中红色标记的测试集样本就属于这种情况。 图2(a)中存在明显偏离斜线y=x(x, y分别为测量值和PLS模型预测值。 样本点越偏离斜线, 说明模型预测精度越差, 反之预测精度越高)的点, 但经过基线校正[图2(b)], 红色点十分贴近斜线且基本均匀分布在斜线的两侧。 一方面样本的均匀分布, 说明KS算法样本划分合理; 另一方面样本点经过基线校正后靠近斜线, 说明分数阶基线校正算法的有效性。 当样本量变大时, 即所有批次合并后, 重新划分校正集测试集, 所得实验结果与理论相符。

4 结论

基于WS平滑算法的各种基线校正算法(包括WS平滑算法)都可以通过微分算子扩展至任意阶次, 具有更好的灵活性和可选择性, 从而进一步提高多元校正模型的预测精度。 γ -PGA发酵实验的光谱数据分析结果表明, 分数阶微分基线校正效果优于整数阶。 同时发现AsLS和FdAsLS基线校正算法在去除基线的同时, 实现了对ATR光谱水峰的扣除, 消除水峰对后续光谱定量的影响, 扩展了该基线校正算法应用范围。 同时为用水溶液测得的ATR光谱消除水峰提供了新的思路。

参考文献
[1] Schorn-García D, Cavaglia J, Giussani B, et al. Microchemical Journal, 2021, 166: 106215. [本文引用:1]
[2] Pu Y Y, O'Donnell C, Tobin J T, et al. International Dairy Journal, 2020, 103: 104623. [本文引用:1]
[3] Cavaglia J, Schorn-García D, Giussani B, et al. Food Control, 2020, 109: 106947. [本文引用:1]
[4] Eilers P H C. Analytical Chemistry, 2003, 75: 3631. [本文引用:2]
[5] Eilers P, Boelens H. Baseline Correction with Asymmetric Least Squares Smoothing, Leiden University Medical Center Report 1, 2005. [本文引用:2]
[6] Zhang Zhimin, Chen Shan, Liang Yizeng. Analyst, 2010, 135(5): 1138. [本文引用:1]
[7] Baek S J, Park A, Ahn Y J, et al. Analyst, 2015, 140(1): 250. [本文引用:1]
[8] JIANG An, PENG Jiang-tao, XIE Qi-wei, et al(姜安, 彭江涛, 谢启伟, ). Computers and Applied Chemistry(计算机与应用化学), 2012, 29(5): 537. [本文引用:1]
[9] He Shixuan, Zhang Wei, Liu Lijuan, et al. Anal. Methods, 2014. [本文引用:1]
[10] Xu D G, Liu S, Cai Y Y, et al. Applied Optics, 2019, 58(14): 3913. [本文引用:1]
[11] Ye Jianfeng, Tian Ziyang, Wei Haoyun, et al. Applied Optics, 2020, 59(34): 10933. [本文引用:1]
[12] Loverro A. Fractional Calculus: History, Definitions and Applications for the Engineer. Report, 2004, Corpus ID: 18520731. [本文引用:1]
[13] Hong Yongsheng, Liu Yaolin, Chen Yiyun, et al. Geoderma, 2019, 337: 758. [本文引用:1]
[14] Xia Zhenzhen, Yang Jie, Wang Jing, et al. Applied Spectroscopy, 2020, 74(4): 417. [本文引用:1]
[15] Fourie E, Aleixand re-Tudo J L, Mihnea M, et al. Food Control, 2020, 115: 107303. [本文引用:1]