基于趋势分析的发酵过程监控

李艳坡， 陈进东， 潘 丰

(江南大学轻工过程先进控制教育部重点实验室，江苏无锡214122)

在补料分批发酵过程中，由于初始条件(包括菌体质量浓度、菌体的初始活性等)和操作条件(例如补料速率)的不同，各批次之间的过程变量的大小往往会有不同;另外，即使在同批次发酵过程的各个阶段(生长阶段、产酸阶段等)之间，过程变量的自身特性和变量之间的相互关系也会变化［1-2］，这些都给发酵状态的识别和诊断带来很大困难。以往的研究中，如数学模型法、基于数据的主元分析法，由于是以发酵的整个过程为基础建模，很难准确反映每个阶段的状态，甚至会出现同样的参数设置在其中一个阶段运行正常，而在另一个阶段则出现误判或漏报的情况［3］。

为此，很多学者以发酵阶段识别为基础对发酵的监控策略做了改进。Soumen K Maiti以主元分析法得到的主元之间的相关系数为基础辨识出发酵过程的各个阶段，但该方法未考虑随机波动的影响［4］。文献［5］以谷氨酸发酵过程中CO2的生产速率和—NH3消耗量在各个发酵阶段的特性，及时间范围为基础实现对各个阶段进行识别。然而，当菌体活性异常时将直接影响上述判断。文献［6-7］通过当前批次当前时段过程变量的变化趋势与之前预定义的状态模式比较，判断出当前的发酵生理状态，并以此实现发酵阶段的识别和相应控制策略，但作者没有找到设定当前时段长度的通用方法。

通过上述分析，文中提出以当前时段变量的变化趋势为基础的发酵过程监控方法，并实现对非正常发酵状态的识别和初步诊断。

1 趋势分析原理

根据文献［8］的分析，发酵的阶段检测可以归结为是对能充分反映发酵状态信息的在线变量的拐点的检测。根据文献［9-10］对代谢流向的分析和菌体对氧气利用量的关系，文中选择尾气中的O2含量来反映当前发酵的状态信息。用数理统计中的线性相关分析［11］可以实现对平滑曲线拐点的检测。

线性相关分析是对X，Y两组数据的线性关系程度分析的一种有效手段，其相关性大小由相关系数表示:

在发酵过程中，在线测量的参数通常会因各种干扰而造成参数的实时曲线伴随不规则的波动，如果不加处理，会因此干扰对拐点的判断。由于发酵过程的大滞后性，采用普通滤波(yout=a×yk－1+(1－a)×yk，a为滤波系数)时，如果滤波系数太小，则滤波效果不明显，如果太大，则会带来明显的滞后，而基于小波分析滤波方法在时间窗长度和核函数的选择上存在不确定性［1］。其实，普通滤波只是考虑到了当前值和上一时刻值之间的连续性，而发酵的当前值是和前一时段不仅仅是和前一时刻有关联，即要考虑当前值和前一时段之间的惯性，以此滤除期间随机的波动。如果采用和小波分析类似的移动时间窗法，则可以将对当前时刻数值的分析拓展到对当前时段的分析。因此，文中以时间窗的时间跨度为单位，通过对该时段的参数(O2)组成的向量和时间序列组成的向量做线性拟合，获取参数当前的变化趋势。线性拟合方程如下:

其中，X为固定时间长度为n的时间序列向量;Y为O2在当前时间窗时段内组成的向量;b为截距，K为线性拟合系数，ε为线性拟合误差。K即是该参数的当前变化趋势。由于发酵过程的连续性，相对于O2的曲线，K的曲线会平滑很多，可以通过相关系数识别其中的拐点。拐点识别的整个过程流程如图1所示。

图1 阶段识别过程流程Fig.1 Process of phase identification

由于发酵过程伴有不确定的干扰，当前时段设定的越长，越有利于分析出参数的真正变化趋势，而不至于被干扰所误导。但如果当前时段设定的太长，会因惯性太大而引发滞后，综合滤波和实时性的考虑，且谷氨酸发酵的每个阶段的时间长度都大于350 min。文中多次比较后在优先保证滤波效果的基础上，将当前时段长度n设定为180(从图2b可以看出当n设定为180时当前时段的变化趋势K已经很平滑，太长则会影响对每个发酵阶段的及时识别)。则 Y= ［yi－179yi－178…yi－1yi］180。图 1 中向量 Kt长度的设定也是在保证对K曲线中拐点的识别不产生明显滞后的前提下，尽可能减小Kt的长度。

本次发酵的识别效果如图2所示。在图a中细线为发酵尾气中O2含量的曲线，粗线为识别曲线，用I表示。在450 min左右的实线尖峰对应I=3+18(为了使识别曲线和O2曲线范围一致，加上18);图b为识别过程，其中曲线1(实线)为O2变化趋势K曲线，曲线2(虚线)为变化趋势向量Kt和时间序列T的相关系数，如图所示在斜率的拐点处相关系数出现了大幅跳跃。

图2 O2拐点识别过程示意Fig.2 Process of inflection point identification

2 谷氨酸发酵过程特点

2.1 发酵培养条件

菌株采用实验室保存的谷氨酸棒状杆菌S9114。发酵罐体积10 L，初始装液量为5 L，接种量为体积分数8%，初始葡萄糖质量浓度为120 g/L。通过流加氨水将pH控制在7.2 ～7.3。通风比固定为1.0。在长菌阶段，将温度设为32℃，溶氧(DO)体积分数控制在10% 左右;到产酸前期，将温度升至34℃，同时将DO体积分数逐步调整到15% ～25%;在产酸后期，将温度设为36℃，将DO体积分数逐步调整到10% ～20%。DO通过PID控制器调节搅拌转速自动控制。

2.2 发酵各个阶段的特点

谷氨酸补料发酵过程可以分为5个阶段，即长菌阶段、过渡阶段、产酸前期、产酸后期和接近放罐阶段。各个阶段的特点及特征如下:

1)长菌阶段。当菌体适应发酵罐环境之后，由于此时底物充足，菌体开始迅速生长，耗氧量也随之迅速增加。为了维持DO的稳定，需要不断提高搅拌速度。由于进气量恒定，整个阶段的转速接近直线上升，发酵液中溶解的氧气增多，对应逃逸到尾气中O2的含量会直线下降。

2)过渡阶段。由于菌体开始转向生产谷氨酸副产物需要暂时调整，菌体活性暂时降低，耗氧量随之下降，为了维持DO恒定，对应转速也会下降，空气溶解到发酵液中的量下降，故逃逸到尾气中的O2含量开始上升。因此，在过渡开始时，尾气中O2会出现一个凹拐点。在菌体调整完成后随着菌体活性的恢复，发酵液中需氧量的增加而使O2继续下降。

3)产酸前期。这期间菌体活性很高，产酸速度也快，菌体耗氧量增加，对应逃逸到尾气中的O2总体上仍会持续下降。

4)产酸后期。随着葡萄糖的限制，菌体含量开始稳定，菌体活性逐渐降低，转速随菌体需氧量的减少而降低，对应逃逸到尾气中的O2此时会缓慢上升。所以在该阶段的初期，尾气中O2也会出现一个凹拐点。

5)接近放罐阶段。此时菌体含量开始下跌，菌体活性明显降低，转速因维持DO的稳定而进一步降低，对应逃逸到尾气中的O2上升速度会比之前加快，所以在这一阶段的初期，O2也会出现一个拐点。整个发酵过程几个相关参数的变化趋势如图3所示。

2.3 各个阶段O2特征的提取

为了实现对谷氨酸发酵过程的在线监控，需要对当前的发酵生理状态做出实时判断，因此首先需要找出每种生理状态的正常变化范围。根据以上分析的各个阶段O2的基本特点以及多批次相同条件下发酵数据的统计，同时为了避免系统过于灵敏而造成误判，对每个阶段的O2正常状态下的特征大致归结如下:

第1阶段记为K1＜0;第2阶段记为K2，且K1/10＜K2＜0;第3阶段记为 K3，且 K1/3＜K3＜K2;第4阶段记为K4，且K4＞0;第5阶段记为K5，且 K5＞ K4。

对每个阶段整体线性拟合后可以得出每个阶段的特征如图3所示。图中对各个变量的范围做了拉伸处理，以使每条曲线都可以清楚显示。

图3 发酵过程中O2的变化特点Fig.3 Character of O2during fermentation process

3 谷氨酸异常发酵批次的识别

在实际生产中，常见的有染菌和因菌体代谢活性过低引起的低产两种不正常类型。如果能及时辨别出发酵过程中的异常状态，则可以避免发酵的失败或提早结束该批次，避免资源的浪费。

3.1 发酵罐染菌的识别

由染菌的特点可知，在发酵罐染菌后由于杂菌(一般都是需氧型)旺盛的代谢能力会使得尾气中O2含量持续下降(K'＜0)(这里假设杂菌尚未占绝对优势，即谷氨酸菌种在各个阶段的过渡过程表现出的拐点特征仍可以通过O2的变化反映)。并且其坡度(K')会小于正常值(人工操作时的判别依据)。K'的获取如图4所示。基于趋势分析识别染菌批次的过程和结果如图5所示。图6为3种不同批次的产酸量对比(采样间隔为2 h)。

从图5中可以看出，在第2阶段末(第2条竖直线，约14 h处)趋势法已识别出染菌症状。由于染菌后的主要特征为产酸量降低［12］，根据图6中染菌批次的离线数据分析，只有在18 h后才能识别出染菌，比趋势法晚约4 h。

3.2 低产批次的识别

对于低产批次，通常由于在产酸阶段菌体活性过低，对O2的需求降低，尾气中O2的坡度(K'值)相对于正常批次会增加至正常范围之外(人工操作的判别依据)。识别的过程和结果如图7所示。可以看出，在第2阶段末(第2条竖直线，约14 h处)趋势法已识别出低产特征，与图6中的离线分析结果(第12～14 h之间出现明显低产)比较接近。

图7 低产批次的识别结果Fig.7 Identification of batch of low product

4 结语

挖掘O2等可在线测量变量中蕴含的发酵状态信息，是研究发酵过程在线监控的重要内容。文中在总结现场人员的操作规律基础上通过线性拟合找出一种基于趋势的定性分析法，并成功用于谷氨酸发酵阶段的识别。以此为基础，根据谷氨酸发酵的特点对各个阶段的尾气中O2含量的特征做了总结，进而定义每个阶段的正常状况，最后通过与非正常发酵特征的对比实现了对非正常发酵过程的初步诊断。

文中的方法从人工操作的规律入手，属于知识工程范畴，它对于实现人工经验转换到程序式的自动控制做了进一步尝试。由于该方法只是从操作经验入手，针对不同的菌种发酵，只需要根据操作人员的经验做相应调整即可，因此具有比较好的推广性。

［1］Sebastien Regis，Andrei Doncescu，Jacky Desachy.Detection and characterization of physiological states in bioprocesses based on Holder exponent［J］.Knowledge-Based Systems，2008，21:70-79.

［2］WEI Su-na，YAN Men-ga，Ahmet Palazoglu，et al.A method for multiphase batch process monitoring based on auto phase identification［J］.Journal of Process Control，2011，21:627-638.

［3］Graells M，Venkatasubramanian V.Fault diagnosis of a benchmark fermentation process:a comparative study of feature extraction and classification techniques［J］.Bioprocess and Biosystem Engineering，2012，35:689-704.

［4］Soumen K Maiti，Rajesh K Srivastav，Mani Bhushan，et al.Real time phase detection based online monitoring of batch fermentation processes［J］.Process Biochemistry，2009，44:799-811.

［5］史仲平，潘丰.发酵过程解析控制与检测技术［M］.2版.北京:化学工业出版社，2010.

［6］Konstantin B Konstantinov，Toshiomi Yoshida.Real-time qualitative analysis of the temporal shapes of(Bio)process variables［J］.AICHE Journal，1992，38(11):1703-1705.

［7］Konstantin B Konstantinov，Toshiomi Yoshida.A knowledge-based pattern recognition approach for real-time diagnosis and control of fermentation processes as variable structure plants［J］.IEEE Transactions on Systems Man and Cybernetics，1991，21(4):908-914.

［8］Kashyap Muthuswamy，Rajagopalan Srinivasan.Phase-based supervisory control for fermentation process development［J］.Process Control，2003，13:367-382.

［9］肖杰.基于呼吸商在线监测的谷氨酸发酵过程控制研究［D］.无锡:江南大学，2008.

［10］储矩，李友荣.现代工业发酵调控学［M］.北京:化学工业出版社，2001.

［11］肖继先.概率论与数理统计［M］.北京:科学出版社，2012.

［12］付丙勇.谷氨酸生产染菌的原因及对策［J］.发酵科技通讯，2012，41(1):23-26.FU Bin-yong.Causes and countermeasures of bacterial contamination during glumute acid production［J］.Fermentation Technology Communication，2012，41(1):23-26.(in Chinese)