柠檬酸被广泛应用于食品、药品、化工等领域，我国是世界上最大的柠檬酸生产国与出口国[1]。我国现有的柠檬酸生产方式主要是发酵法，将玉米、木薯等粉碎、液化后进行带渣发酵。然而带渣的混液含有大量粗蛋白，容易导致菌体疯长，糖酸转化率低[2];同时带渣玉米液化液较黏稠、溶氧差，生产能耗大。而利用去除残渣后的液化清液进行发酵，可以有效克服以上不足，清液黏度低利于溶氧，降低能耗，提高转化率[3]。监控柠檬酸发酵液化清液的生产过程对整体柠檬酸发酵过程至关重要，而传统方法只能通过对产品中各种理化值进行事后检验并由相关工艺人员判断清液生产过程是否正常。
近红外光谱从分子振动层面收集了大量过程信息，因其非破坏性、分析快、效率高而引起人们的广泛关注，现已在农业、石油、医药、食品和环境等领域得到应用[4-8]。但是现有的利用近红外光谱的监控方法是将测量得到的过程变量信息转化为理化值并与经验得到的理化值阈值对比进行监控[9-13]，这样忽略了近红外光谱本身的统计特性，导致漏报等很多异常状况出现。本文提出分析近红外光谱的统计特性，通过判断光谱的分布信息来判断光谱表示的生产过程是否异常，这样能够充分利用近红外不同波长上的信息，并根据统计分布的变化实现事前预警。
本文将基于概率偏最小二乘(probability partial least squares，PPLS)的多元统计过程监控方法与近红外光谱技术相结合来对柠檬酸发酵液化清液生产过程进行监控。通过对柠檬酸发酵液化清液光谱进行联合区间概率偏最小二乘法(synergy interval probability partial least squares，siPPLS)波段优选、建立PPLS模型、计算光谱的分布特性、建立统计监控指标与置信限对比，得到测试数据是否存在异常。最后与工艺人员给出的故障数据进行比对，得到最后漏报率与错报率，并与传统利用近红外光谱进行监控的方法进行对比。
1 材料与方法
1.1 样本制备
实验所用样本是某柠檬酸工厂生产过程中产生的发酵液化清液，按不同批次，不同日期取样，每天取8个样本，共计得236组数据。将提取到的清液样品以500 r/min搅拌1 min，谱图稳定。液化清液中的总糖由费林法测得，总氮由凯氏定氮法测得。
1.2 光谱数据采集及数据处理
近红外光谱仪(MATRIX-F型傅里叶近红外光谱仪，德国Bruker公司)用于液化清液光谱的采集，光谱波长范围为3 996～11 988 cm-1，光谱扫描分辨率为16 cm-1，扫描次数为64次。光谱采集软件为OPUS，采集的每条光谱包含1 037个数据点。预处理软件：The Unscrambler X 10.3,挪威CAMO公司。波段选择及模型构建软件:Matlab 7.11,美国MathWorks公司。
主成分分析(principal component analysis，PCA)[14-15]、偏最小二乘法(partial least squares，PLS)[16]是常见的特征提取方法，PLS在PCA的基础上考虑到了输入输出的回归关系。LI[17]等将概率分布引入PLS模型提出PPLS方法，考虑到每个变量的概率分布情况，在主元和误差都服从高斯分布的条件下，通过求解极大似然函数得到主元、残差等模型参数。
1.3 PPLS模型
PPLS是一种高斯隐变量模型，假设经过归一化后数据的过程变量和输出变量分别为
分别为过程的过程变量和输出变量个数。PPLS模型如公式(1)、公式(2)所示：
xn=Ptn+μx+ξn
(1)
yn=Ctn+μy+εn
(2)
式中：P∈RDx×K， P∈RDx×K,K<Dx是因子个数，μx、μy分别为X、Y的均值，残差
主元tn～N(0,IK)。
可用贝叶斯公式求tn的后验分布，如公式(3)所示：
(3)
式中：
由PPLS模型可知，主元满足标准正态分布，测量数据由均值、主元组合及未知噪声生成。模型中需要被估计的未知参数记做
可由交叉验证方法得到。将主元tn当作隐变量处理，用最大期望算法(expectation maximization algorithm，EM)求解未知参数。EM及收敛证明见文献[18]。具体求解过程也可参考文献[17]。最后得到yn的近似值[yn|xn]如公式(4)所示：
yn≈E(yn|xn)=CE(tn|xn)+μy=
(4)
1.4 siPPLS及PPLS模型评价指标
模型评价是判断模型建立好坏的重要指标。评价化学计量学模型一般用均方根误差，均方根误差越小，模型性能越好。校正集均方根误差(root mean square error of calibration，RMSEC)、预测均方根误差(root mean square error of prediction，RMSEP)、交叉验证均方根误差(root mean square error of cross validation，RMSECV)计算如公式(5)～公式(7)所示：
(5)
(6)
(7)
式中：
表示第i个样品建立模型时的估计值，
表示第i个样品用化学方法测量的真实值，nc表示样品数。
2 结果与分析
2.1 光谱预处理
光谱采集过程会受到测量条件、外部环境等影响而产生噪声，因此在分析数据前要对近红外光谱进行预处理，通过预处理可以减少背景噪声、基线漂移等对光谱的干扰。一阶导数、二阶导数、多元散射校正(multiplicative scatter correction，MSC)、平滑是常用预处理方法。本文用一阶导数、MSC、9点SG(Savitzky-Golay)平滑分别对原始光谱进行预处理。原始光谱及预处理后光谱如图1所示。SG平滑可以有效消除噪声但是平滑图(图1-b)没有消除光谱基线漂移，导数可以有效消除基线漂移和背景干扰，但是会放大信噪比；一阶导数图(图1-d)虽然消除了光谱的基线漂移，但是会导致光谱吸光度变化太小，不利于建模；MSC图(图1-c)既消除了基线漂移，吸光度变化也更适合建模，同时实验所用清液是过滤所得，会有大小不同的颗粒，颗粒大小会影响光谱，而MSC可以用来减小颗粒大小及分布不均匀产生的散射对近红外光谱的影响，综合考虑，本文选用MSC预处理方法。
a-原始光谱;b-SG平滑处理光谱;c-MSC处理光谱；d-阶导数处理光谱
图1 原始光谱及预处理后的柠檬酸发酵液化清液光谱
Fig.1 Original spectra and preprocessed near infrared spectra of clear solution for citric acid fermentation
2.2 校正集样本的选择
通过上述方法一共获得2组数据，第一组171组正常数据，第二组24组正常数据与31组异常数据混合。根据约2∶1在第一组正常数据中取111组数据作为校正集，60组作为验证集，第二组正常数据及异常数据混合共55组作为过程监控方法的样本集。由表1、表2可以看出，验证集的分组数据在校正集内，分组合理。
表1 Kennard-Stone分组结果(总糖)
Table 1 Results of Kennard-Stone(total sugar)
校正集验证集样品数最大值最小值平均值标准差样品数最大值最小值平均值标准差11120.5515.5817.981.356020.3215.717.951.32
表2 Kennard-Stone分组结果(总氮)
Table 2 Results of Kennard-Stone(total nitrogen)
校正集验证集样品数最大值最小值平均值标准差样品数最大值最小值平均值标准差1110.0660.0360.049 90.008600.0630.0450.049 50.006
2.3 光谱波段选择
常用的波长选择方法有无信息变量消除方法(uniformative variable elimination,UVE)、区间偏最小二乘(interval partial least squares,iPLS)[19]与联合区间偏最小二乘(synergy interval partial least squares,siPLS)[20-21]等。本文在PPLS的基础上，参考由NORGAARD等[22]提出的iPLS方法及siPLS提出siPPLS来选择波段。
siPPLS是建立在PPLS基础上的波段选择方法，siPPLS具体算法步骤如下:(1) 将整个光谱区域分成n个等宽的子区间;(2) 将数量为 2、3、4的所有等宽子区间排列组合建立PPLS回归模型; (3) 将不同子区间组合的变量建立PPLS模型并计算每个模型的RMSECV;(4) 选出最小的RMSECV对应的子区间组合对应的光谱波段。
由于还不能从理论上确定参加联合子区间间隔的数目，在应用siPPLS时，尝试将整个光谱区域分别划分为 10、11、12、…、25 个子区间。分别计算总糖、总氮的各种区间划分及子区间最佳组合下建立模型的RMSECV。对于总糖，如表3所示，siPPLS波长选择算法在光谱划分为18个子区间并按[7,11,14,15]区间组合选择波长时能获得最小的RMSECV，其值为0.347 8。对于总氮，如表4所示，siPPLS波长选择算法在光谱划分为20个子区间并按[4,5,10,16]区间组合选择波长时能获得最小的RMSECV，其值为2.400 5×10-3。
由表3可知，经过siPPLS波段选择，建立光谱与含糖量相关的PPLS模型时，光谱波段不是全部光谱，而是将光谱分为18个区间时的第7、11、14、15个区间的组合对应的光谱，即在全部光谱的1 037个数据点中343～399、571～627、742～855的数据点。
表3 siPPLS方法划分的不同间隔数的特征
光谱区间筛选结果(总糖)
Table 3 Spectral interval screening results of different
intervals of siPPLS method (total sugar)
区间数主成分数最佳间隔组合RMSECV109[ 3 5 8]0.388 81110[1 2 5 9]0.367 71210[ 5 7 10]0.372 4138[5 6 7 11]0.366 3149[ 2 7 8 11]0.374 41510[ 2 8 12]0.369 81610[ 2 7 8 13]0.366 7179[ 5 10 13 14]0.361 41810[7 11 14 15]0.347 8199[1 3 8 15]0.351 4209[3 10 12 16]0.350 12110[3 11 16 17]0.352 72210[3 11 17 18]0.354 32310[8 12 18 19]0.350 92410[2 12 19 20]0.353 32510[3 12 13 15]0.349 9
表4 siPPLS方法划分的不同间隔数的特征
光谱区间筛选结果(总氮)
Table 4 Spectral interval screening results of different
intervals of siPPLS method (total nitrogen)
区间数主成分数最佳间隔组合RMSECV109[ 3 5 8]2.408 6×10-31110[ 3 6 9]2.405 5×10-31210[4 5 6 7]2.402 8×10-31310[1 6 7 10]2.411 0×10-3149[4 7 11]2.407 1×10-31510[4 8 11]2.404 2×10-31610[4 6 7 13]2.405 5×10-3179[7 9 10 13]2.403 1×10-31810[5 9 10 14]2.400 8×10-3199[5 9 10 15]2.402 6×10-32010[4 5 10 16]2.400 5×10-32110[5 15 16 17]2.402 1×10-32210[6 10 12 17]2.402 6×10-32310[6 11 12 18]2.404 2×10-32410[5 12 18 19]2.401 6×10-32510[6 16 19 20]2.402 5×10-3
由表4可知，经过siPPLS波段选择，建立光谱与含氮量相关的PPLS模型时，光谱波段不是全部光谱，而是将光谱分为20个区间时的第4、5、10、16个区间的组合对应的光谱，即在全部光谱的1 037个数据点中154～255、460～510、766～816的数据点。
2.4 PPLS建模
将经过预处理以及波长选择的光谱数据进行PPLS建模。通过RMSEC、RMSEP对模型的准确性进行评价。图2、图3分别为校正集的PPLS模型及RMSEC指标、验证集的PPLS模型及RMSEP评价指标。由图2、图3可知，校正集与验证集的均方根误差之比在0.8～1.2，说明建立的PPLS模型是有效的。
a-含糖量拟合值与真实值对比；b-含氮量拟合值与真实值对比
图2 PPLS模型的样品拟合值与真实值对比图(校正集)
Fig.2 Correlation of predicted and measured value
based on PPLS
a-含糖量拟合值与真实值对比；b-含氮量拟合值与真实值对比
图3 PPLS模型的样品拟合值与真实值对比图
Fig.3 Correlation of predicted and measured value
based on PPLS
2.5 基于近红外的PPLS监控指标
现有基于近红外的监控通常建立近红外光谱与某种理化值之间的回归关系，通过判断理化值是否在阈值内判定过程是否异常。但是这样忽略了近红外光谱本身的统计特性，会降低监控结果的准确度。
本文提出通过分析近红外光谱的统计特性，对柠檬酸发酵液化清液过程进行监控。统计监控理论认为在随机干扰下，过程变量满足某种正态分布，如果过程出现故障该分布会变化。光谱不同波长及其对应的吸光度都可以看做不同的过程变量，这些过程变量本身都满足不同的正态分布，即在每个波长点处的所有171个样本都符合一定的分布，如图1-a中样本点1和2的分布如图4、图5所示，图4是波长在7 200 cm-1处所有171个样本点的分布图，图5是波长在8 904 cm-1处所有171个样本点的分布图，这些分布对于判断过程运行是否正常至关重要，而只用质量变量不足以反映出过程的这些信息。本文考虑了这些过程变量本身的分布情况，将特征提取方法PPLS与统计监控方法结合起来。根据建模的数据得到正常运行状况下的主元应该满足的正态分布，通过监控指标计算测试样本的主元与主元正常分布之间的马氏距离。具体地，首先得到样品的光谱xn和对应的理化值含糖量和含氮量yn，对校正集中的样本进行PPLS建模，可以得到公式(3)及公式(4)所需的参数值，根据公式(4)可以计算出校正集及验证集的含糖量及含氮量的预测值，计算校正集的RMSEC与验证集的RMSEP，判断模型是否可靠，其次在建立可靠PPLS模型的基础上可以得到主元相关参数如公式(4)，由于主元tn的测量值在建模时得不到，故用tn的后验分布的期望值即公式(3)的期望值来代替tn。根据概率统计常识，基于马氏距离的tn满足卡方统计规律。因此，可以计算出测试集的监控指标GT2如公式(8)所示：
(8)
式中：监控指标GT2服从置信度为α、自由度为k的卡方统计规律，自由度即为PPLS模型中的主元个数。
图4 样本点1的分布图
Fig.4 Distribution of sample point 1
图5 样本点2的分布图
Fig.5 Distribution of sample point 2
2.6 基于近红外的PPLS模型监控结果
现有基于近红外的监控方法通常建立近红外光谱与质量变量之间的回归关系，将光谱信息转化为质量变量，然后通过判断质量变量是否在阈值内判断过程是否正常。基于此柠檬酸发酵液化清液的监控结果如图6所示，图中光谱信息首先转化为总糖含量，然后对总糖含量进行监控，得到的糖含量全部在阈值范围内，但是实际过程存在异常。在实际工业过程中，总氮并没有设定严格的阈值，故不需要将近红外光谱转化为总氮进行监控，因此将近红外光谱转化为质量变量的方法忽略了很多有用过程信息，不能有效监控生产过程。
本文直接计算近红外光谱的统计特性，统计监控理论认为在正常情况下，过程变量满足正态分布，如果过程出现故障则该分布发生变化。本文首先对近红外光谱进行PPLS特征提取，然后通过公式计算出各个测试样本关于总糖、总氮的监控指标GT2的值，结果如图7、图8所示，将PPLS模型应用于近红外光谱的监控中，结合监控指标发现第8～17、第20～31、第34～39、第42～47及第52～53的指标超过控制线，工艺人员给出的故障数据为第1、第6～17、第20～31及第42～47。PPLS方法得到的监控指标漏报第1、6、7三个数据，错报第34、35、36、37、38、39、52、53八个数据，漏报率为9.68%，错报率为25.81%。相比于将近红外光谱信息转化为质量变量的方法，漏报率与错报率都有了很大的改进。
图6 近红外光谱转化为总糖的监控结果
Fig.6 Monitoring result of conversion of near-infrared
spectroscopy to total sugar
图7 测试样本的监控指标(总糖)
Fig.7 Monitoring indices of test samples (total sugar)
图8 测试样本的监控指标(总氮)
Fig.8 Monitoring indices of test samples (total nitrogen)
3 结论
针对柠檬酸发酵液化清液的生产过程的监控，本文从近红外光谱的统计特性层面提出了基于光谱的PPLS监控方法。过程中采用siPPLS方法选波段所得RMSECV为0.347 8、2.400 5×10-3，在理想范围内，说明siPPLS优选波段的方法是有效的。对近红外光谱的PPLS统计监控比传统方法将近红外光谱转化为质量变量得到的监控结果更为有效，漏报率为9.68%，错报率为25.81%，相较于传统方法全部检测不出来有了很大的提升，说明基于近红外光谱的PPLS过程监控是有效的，为实际柠檬酸工业过程监控提供了一种可靠的参考方法。
[1] 王宝石,陈坚,孙福新,等.发酵法生产柠檬酸的研究进展[J].食品与发酵工业, 2016,42(9):251-256.
[2] 潘声龙.玉米清液发酵生产柠檬酸的工艺研究[J].安徽农业科学, 2010,38(29):16 700-16 701.
[3] 孙荣.柠檬酸新型清液发酵工艺研究[D].济南:山东轻工业学院, 2011.
[4] ALVE J C, POPPI R J.Biodiesel content determination in diesel fuel blends using near infrared (NIR) spectroscopy and support vector machines (SVM)[J].Talanta, 2013,104(2):155-161.
[5] SASKIA L, MANFRED S, ROSNER S.Hydraulic traits of Norway spruce sapwood estimated by Fourier transform near-infrared spectroscopy (FT-NIR)[J].Canadian Journal of Forest Research, 2015,45(6):625-631.
[6] 战皓,方婧,杨滨,等.近红外光谱法测定不同产地独活中蛇床子素和二氢欧山芹醇当归酸酯含量[J].光谱学与光谱分析, 2017,37(4):1 110-1 113.
[7] 狄天云,高晓娟,张霞,等.基于指标成分与近红外光谱对宁夏野生和栽培甘草的比较鉴别研究[J].中国中药杂志, 2017,42(3):536-541.
[8] 王元忠,赵艳丽,张霁,等.近红外光谱信息筛选在玛咖产地鉴别中的应用[J].光谱学与光谱分析, 2016,36(2):394-400.
[9] ESPINOSA F, ANTONIO J.Combining PLS regression with portable NIR spectroscopy to on-line monitor quality parameters in intact olives for determining optimal harvesting time[J].Talanta, 2016,148:216-228.
[10] MIRSCHEL G, HELMSTEDT U, SCHERZER T, et al.Monitoring of the degree of condensation in alkoxysiloxane layers by nir reflection spectroscopy[J].Industrial & Engineering Chemistry Research, 2016,53(43):16 813-16 819.
[11] WEI Luo,JIA Wu,WANG Xuekai,et al.Near infrared spectroscopy combination with PLS to monitor the parameters of naproxen tablet preparation process[J].Analytical Methods, 2013,5(5):1 337-1 345.
[12] KANG Qian,RU Qingguo,LIU Yan,et al.On-line monitoring the extract process of Fu-fang Shuanghua oral solution using near infrared spectroscopy and different PLS algorithms[J].Spectrochim Acta A, 2016,152:431-437.
[13] 王永香,郑伟然,米慧娟,等.热毒宁注射液青蒿金银花浓缩过程近红外快速定量检测方法的建立[J].中草药, 2017，48(1):102-108.
[14] REED J P, DEVLIN D, ESTEVES S R, et al.Integration of NIRS and PCA techniques for the process monitoring of a sewage sludge anaerobic digester[J].Bioresour Technol, 2013,133(4):398-404.
[15] ZHOU Bo,YE Hao,ZHANG Haifeng,et al.Process monitoring of iron-making process in a blast furnace with PCA-based methods[J].Control Engineering Practice, 2016,47:1-14.
[16] HUCHO F, OBERTHUR W, LOTTSPEICH F.PLS-based EWMA fault detection strategy for process monitoring[J].Journal of Loss Prevention in the Process Industries, 2015,36(1):108-119.
[17] LI S, GAO J, NYAGILO J O, et al. Probabilistic partial least square regression: A robust model for quantitative analysis of raman spectroscopy data[C]. San Diego: IEEE International Conference on Bioinformatics & Biomedicine. 2012.
[18] STOCIA P, XU L Z,LI J.A new type of parameter estimation algorithm for missing data problems[J].Statistics & Probability Letters, 2005,75(3):219-229.
[19] RADY A M, GUYER D E.Evaluation of sugar content in potatoes using NIR reflectance and wavelength selection techniques[J].Postharvest Biology & Technology, 2015,103:17-26.
[20] LI Yang,GUO Mingye,SHI Xinyuan,et al.Online near-infrared analysis coupled with MWPLS and SiPLS models for the multi-ingredient and multi-phase extraction of licorice (Gancao)[J].Chinese Medicine, 2015,10(1):1-10.
[21] 郭慧娴,朱思祁,黎远鹏,等.基于iPLS和SiPLS算法的人体血清胆红素含量的可见-近红外光谱建模[J].光电子·激光, 2016,27(10):1 136-1 144.
[22] NORGAARD L, SAUDLAND A, WAGNER J, et al.Interval partial least-squares regression(iPLS): A comparative chemometric study with an example from near-infrared spectroscopy[J].Applied Spectroscopy, 2000,54(3):413-419.