CN106546819B - 一种量测滞后下的发酵过程生物量在线检测方法 - Google Patents

Abstract

一种量测滞后下的发酵过程生物量在线检测方法，通过容积卡尔曼滤波对发酵过程生物量进行状态估计；利用离线分析结果在主要量测采样时刻进行容积卡尔曼滤波量测更新和量测修正；经过递推计算将估计结果从采样时刻递推到主要时刻；求取主要量测包含的特有信息；将主要量测含有的特有信息和次要量测估计结果进行协方差融合，获得既包含次要量测信息又包含量测滞后信息的状态估计结果，实现量测滞后下的发酵过程生物量在线检测。本发明在基于容积卡尔曼滤波器的发酵过程软测量技术基础上，采用协方差融合算法将量测滞后生物量信息引入发酵过程状态估计过程，提高了发酵过程生物量在线检测精度。

Description

一种量测滞后下的发酵过程生物量在线检测方法

技术领域

本发明涉及一种发酵过程生物量在线检测方法，属于信息数据处理及智能检测技术领域，尤其涉及一种基于容积卡尔曼滤波和协方差融合算法的量测滞后下发酵过程生物量在线检测方法。

背景技术

生物发酵过程是生物工程和现代生物技术及其产业化的基础，其机理十分复杂，具有高度的非线性、不确定性，不易实现关键生物量(如菌体浓度、代谢产物浓度和菌体比生长速率等)的在线检测，严重影响着发酵过程在线优化控制技术的实施，制约着先进控制算法和策略在发酵过程的实际应用。

目前，发酵过程中关键生物量的检测主要分为离线检测和在线检测。离线检测是经过采样进行离线分析，常用的离线分析方法有干重法、直接染色法和细胞计数法等，离线检测结果具有较高的精度，但采样间隔时间较长，检测结果存在滞后，难以满足发酵过程实时监测和控制要求。在线检测方法主要有硬传感方法和软测量方法。基于硬件传感器的生物量在线检测方法主要有光浊度法、荧光法、介电常数法、超声波法等，这类在线检测方法的仪器成本和维护费用高昂，而且有特定的适用范围，在实际的发酵过程应用中受到限制。软测量技术是解决发酵过程生物参量在线检测的一种非常实用的方法，利用系统辨识、神经元网络或支持向量机等方法建立发酵过程模型，进而通过可在线检测的过程参量(温度、pH值、溶氧浓度等)对生物量进行在线估计。但由于发酵过程的复杂性和模型不确定性，导致软测量精度较低，基于软测量方法的发酵过程生物量在线检测精度有待进一步提高。因此，在基于状态估计的软测量技术基础上，将离线分析获得的滞后生物量信息加以利用，能够进一步提高发酵过程生物量在线检测精度。

本发明在基于容积卡尔曼滤波器的发酵过程状态估计基础上，采用协方差融合算法将量测滞后生物量信息引入发酵过程状态估计过程，提出一种量测滞后下的发酵过程生物量在线检测方法，进而提高发酵过程生物量在线检测精度，能够为发酵过程监测与过程优化控制技术及其应用提供支持，在轻工、医药、食品、农业等领域具有很好的应用前景。

发明内容

本发明以提高量测滞后下的发酵过程生物量在线检测精度为目的。

发酵过程的连续时间非线性状态空间模型一般由一组微分方程组表示，对发酵过程的系统状态方程组和量测方程组进行离散化处理，得到的非线性系统就可以采用非线性滤波算法进行状态估计。考虑离散化得到的发酵过程非线性系统方程如下：

x_k+1＝f(x_k，u_k)+w_k

其中x_k是发酵过程的系统状态向量，包括n个状态量，u_k是输入变量，f(·)是表征发酵过程模型的非线性函数，w_k是协方差为Q均值为0的系统噪声。

发酵过程中能够规律地、实时地获得的过程参量称为次要量测，规律地、实时地获得的过程参量为温度、pH值等；需要人工离线分析或计算机自动控制取样分析的过程参量称为主要量测，需要人工离线分析或计算机自动控制取样分析的过程参量为菌体浓度、产物浓度和菌体比生长速率等。

次要量测的方程如下：

y_k＝h⁰(x_k)+v_k，y_k为次要量测函数的应变量。

其中h⁰(·)是次要量测函数，v_k是均值为0方差为R₀的量测噪声。次要量测以固定的频率在每个时刻k进行采样，采样值能够实时到达数据处理单元。

量测滞后的主要量测方程为：

z_g＝h¹(x_s)+r_s，z_g为主要量测函数的应变量。

其中g＝s+d，h¹(·)是主要量测函数，r_s是均值为0方差为R₁的量测噪声。主要量测的目标是离线分析获得的生物量，s时刻是主要量测的采样时刻，分析结果延迟一段时间d并在主要时刻g到达数据处理单元，主要量测采样时刻s和延迟时间d在主要时刻g才已知。

针对上述系统方程和量测方程描述的量测滞后下发酵过程，一种量测滞后下的发酵过程生物量在线检测方法包括以下步骤：

步骤一 设定系统初始参数

进行初始参数设定，包括初始时刻的估计值初始时刻系统协方差协方差阵P_0|0。

步骤二 进行容积卡尔曼滤波时间更新并保存结果。

根据上一时刻的协方差阵P_k-1|k-1、状态估计值计算k时刻的容积点：

S_k1|k-1＝(chol(P_k-1|k-1))^T

其中S_k-1|k-1是P_k-1|k-1的平方根，chol表示Cholesky分解，i＝1，2，3，...，m，m＝2n。[1]_i表示[1]的第i列。[1]表示一组由n维单位列向量全排列和变换符号构成的全对称点集；若则[1]代表：

计算经系统方程非线性传播后的容积点状态预测值及误差协方差P_k|k-1：

其中X_i，k-1|k-1是容积点，u_k-1是k-1时刻的输入变量，ω是各个容积点的权值，Q是系统噪声协方差。

为了在主要量测到达后利用主要量测信息在主要量测采样时刻进行量测更新，而主要量测采样时刻在主要时刻g才已知，因此需要保存上述仅利用次要量测进行容积卡尔曼滤波得到的每个时刻时间更新获得的状态预测值误差协方差P_k|k-1。

步骤三 利用次要量测进行容积卡尔曼滤波量测更新

计算容积点：

S_k|k-1＝(chol(P_k|k-1))^T

计算经次要量测方程传播后的容积点：

Y_i，k|k-1＝h⁰(X_i，k|k-1)

计算量测预测值：

计算自协方差阵和互协方差阵：

其中是容积点的权值，R₀是次要量测噪声的方差。

步骤四利用次要量测信息进行量测修正并保存结果

W_y，k、P_y，k|k分别为利用次要量测信息计算得到的卡尔曼增益阵、状态预测值和误差协方差阵。

为了完成后续步骤求取主要量测含有的特有信息，而主要量测采样时刻在主要时刻才已知，需要保存每个时刻量测修正得到的P_y，k|k。

步骤五进行下一时刻容积卡尔曼滤波或者利用主要量测进行量测更新和量测修正

如果主要量测信息还未到达，量测信息仅有次要量测，则返回步骤二，进行下一时刻的容积卡尔曼滤波。

如果主要量测信息到达，即主要量测信息已经可用，那么当前时刻k是主要时刻，且知道了主要量测采样时刻s，则从步骤二中已保存数据中调取s时刻的时间更新结果，并利用主要量测值，进行主要量测采样时刻的量测更新和量测修正。

计算容积点：

S_s|s-1＝(chol(P_s|s-1))^T

计算经主要量测方程传播后的容积点：

Z_i，s|s-1＝h¹(X_i，s|s-1)

计算量测预测值：

计算自协方差和互协方差：

其中是容积点的权值，R₁是主要量测噪声的方差。

利用主要量测进行量测修正：

W_z，s＝P_xz，s|s-1(P_zz，s|s-1)^-1

其中W_z，s|s、P_z，s|s分别为利用主要量测信息计算得到的主要量测采样时刻s的卡尔曼增益阵、状态预测值和误差协方差阵。

步骤六将主要量测采样时刻s对应的结果递推到主要时刻g。

将s时刻的状态估计值和误差协方差P_z，s|s从s时刻进行时间更新：

S_z，s|s＝(chol(P_z，s|s))T

对状态估计值误差协方差P_z，s+1|s再次进行时间更新，共进行延迟时间d次时间更新，最终递推得到主要时刻g的状态估计值误差协方差P_z，g|s。

至此，在主要时刻有两个状态估计值，分别是利用主要量测并递推得到的状态估计值和仅利用次要量测得到的状态估计值这两个状态值包含了部分共有信息，这部分共有信息是将s时刻仅利用次要量测的P_y，s|s经时间更新递推到主要时刻g的状态值误差协方差P_y，g|s。

步骤七在主要时刻g进行协方差融合。

计算主要量测所含有的特有信息，其协方差阵和状态估计值如下：

把主要时刻g仅利用次要量测的容积卡尔曼滤波结果和主要量测含有的特有信息进行协方差融合：

上述P_g|g、即为主要时刻新的协方差阵和状态估计值，将此结果带回到步骤二即可进行下一时刻的状态估计。

以上就是量测滞后下的发酵过程生物量在线检测方法步骤。

本发明的有益效果在于：针对发酵过程这种高维复杂非线性系统中的生物量在线检测问题，将容积卡尔曼滤波算法和协方差融合算法结合，提出一种量测滞后下的发酵过程生物量在线检测方法。本发明将发酵过程中滞后的主要量测信息加以利用，在主要时刻采用协方差融合算法将主要量测的特有信息和仅利用次要量测的容积卡尔曼滤波结果进行融合，结果既包含了次要量测信息又包含了滞后的主要量测信息，从而有效地提高了量测滞后下发酵过程生物量在线检测的精度。

附图说明

图1是本发明所述的发酵过程量测滞后问题示意图。上部等间隔的矩形块表示次要量测y，其采样率固定且实时到达数据处理单元。非等间隔的矩形块表示主要量测z，其采样率不固定，s时刻是主要量测的采样时刻，分析结果随机延迟一段时间d，于主要时刻g到达数据处理单元。

图2是本发明所述的量测滞后下的发酵过程生物量在线检测算法流程图；

图3、4是本发明具体实施方式所述仿真实验的生物量在线检测结果图；

图5、6是本发明具体实施方式所述仿真实验的RMSE(均方根误差)结果图。

具体实施方式

下面结合实例及附图对本发明作进一步的描述，需要说明的是，实施例并不限定本发明要求保护的范围。

具体实施方式：发酵过程中的量测滞后问题描述如图1所示。选用量测滞后下的连续发酵罐模型进行仿真实验，仿真环境是Matlab 2015a，仿真长度200小时，步长0.1小时，仿真总长2000步。

该发酵过程的机理模型为：

系统维数n＝3，状态量x＝[X，S，P]，变量X、S、P分别表示菌体浓度、底物浓度、产物浓度，单位均为g/L。

输入量u＝[D，S_f]，D表示稀释率(h^-1)，输入量S_f表示流加底物浓度(g/L)，D和S_f分别等于0.2h^-1和20g/L。μ是P和S的函数，表示单位生长率。

Y_X/S、α、β是产量参数，μ_m是最大单位生长率，P_m、K_m、K_i是反应常数，上述参量的取值如表1所示。

表1 模型参数取值

采用欧拉法得到离散的系统方程，离散化步长h＝0.1小时。

X_k+1＝X_k+h(-DX_k+μ_kX_k)+w_1，k

S_k+1＝S_k+h(DS_f-S_k)-μ_kX_k/Y_x/s)+w_2，k

P_k+1＝P_k+h(-DP_k+(αμ_k+β)X_k)+w_3，k

其中w_k是系统噪声，其均值为0，方差为Q＝diag(0.001，0.001，0.1)，单位是g/L。

系统的次要量测方程如下：

y_k＝h⁰(x_k)+v_k＝[0 0 1]x_k+v_k＝P_k+v_k

P_k是k时刻状态变量中的产物浓度，采样间隔为0.1小时，每个时刻采样一次且采样值实时到达数据处理单元，v_k是均值为0方差为R₀＝[5]的量测噪声，单位是g/L。

主要量测方程为：

z_g＝h¹(x_k)+r_s＝[1 0 0]x_s+r_s＝X_s+r_s

X_s是主要量测采样时刻状态变量中的菌体浓度，需要经过实验室离线分析获得，采样频率较低且采样间隔不固定，分析结果滞后一段时间才能获得。其中g＝s+d，表示某一次的主要量测于s时刻完成采样，延迟一段时间d，于g时刻到达数据处理单元。r_s是均值为0方差为R₁＝[0.1]的量测噪声，单位是g/L。根据系统方程和量测方程，生成仿真数据，包括系统真值、次要量测值、主要量测值。生成仿真时间序列，包括主要量测采样时刻序列、延迟时间序列、主要时刻序列，其中主要量测每一小时左右(10个步长加上一个1～3的随机步数)采样一次，延迟时间d是一组在0～9之间的随机整数倍步长，主要时刻是主要量测采样时刻加上延迟时间。

采用量测滞后下的发酵过程生物量在线检测方法对上述连续发酵过程模型进行生物量在线检测，具体实现如下：

1、设置初始滤波参数

设置初始状态初值x(：，1)＝[1；5；1]，估计初始值协方差阵P₀＝diag(1，1，1)。

2、进行容积卡尔曼滤波时间更新

根据上一时刻的协方差阵求取容积点及其权值(i＝1，2，3，4，5，6)：

S_k-1|k-1＝chol(P_k-1|k-1)^T

其中此连续发酵过程模型的系统维数n＝3，式中[1]代表所以

计算经系统非线性函数传播后的容积点一步状态预测值预测误差协方差阵P_k|k-1：

保存每个时刻的上述两项结果。

3、仅利用次要量测进行容积卡尔曼滤波量测更新

计算容积点和经次要量测方程传播后的容积点(i＝1，2，3，4，5，6)：

S_k|k-1＝chol(P_k|k-1)^T

Y_i，k|k-1＝h⁰(X_i，k|k-1)＝X_i，k|k-1(3，：)

计算量测预测值、自协方差阵和互协方差阵：

4、进行量测修正

计算卡尔曼增益阵：

计算仅利用次要量测进行容积卡尔曼滤波得到的状态预测值和误差协方差阵P_y，k|k：

保存每个时刻的上述两项结果。

5、如果当前时刻不是主要时刻，主要量测信息不可用，则返回步骤2继续进行下一时刻的容积卡尔曼滤波。如果当前时刻是主要时刻，则进行协方差融合处理。

在主要时刻，主要量测信息已到达数据处理单元且主要量测采样时刻s和延迟时间d已知，从上述步骤2保存的结果中读取主要量测采样时刻的状态预测值预测误差协方差阵P_s|s-1，根据这两项信息和主要量测值进行量测更新以及量测修正：

S_s|s-1＝chol(P_s|s-1)^T

Z_i，s|s-1＝h¹(X_i，s|s-1)＝X_i，s|s-1(1，：)

W_z，s＝P_xz，s|s-1(P_zz，s|s-1)^-1

得到新的采样时刻的状态估计值和协方差阵P_z，s|s。

6、将采样时刻的结果递推到主要时刻：

将利用次要量测获得的s时刻的状态估计值、误差协方差从s时刻进行时间更新递推到g时刻，得到和P_y，g|s。

将利用主要量测获得的s时刻的状态估计值、误差协方差从s时刻进行时间更新递推到g时刻，得到和P_z，g|s。

7、求取主要量测包含的特有信息并进行协方差融合：

将此结果带回步骤2进行下一时刻的量测滞后下容积卡尔曼滤波，直至k等于仿真步长。

上述步骤即可完成该量测滞后下的连续发酵过程生物量在线检测，量菌体浓度、底物浓度在线检测曲线分别如图3、4所示，实线是仿真真实值，虚线是仅利用次要量测的容积卡尔曼滤波算法结果，点划线是本发明所述的量测滞后下容积卡尔曼滤波算法结果。进行50次蒙特卡洛仿真实验，对比仅利用次要量测的容积卡尔曼滤波算法与本发明所述算法的均方根误差(RMSE)，图5，6分别是菌体浓度、底物浓度的均方根误差对比图。实验结果表明，本发明提出的量测滞后下的发酵过程生物量在线检测算法能够有效利用量测滞后信息，有效地提高了发酵过程生物量软测量的检测精度。

Claims (1)

1.一种量测滞后下的发酵过程生物量在线检测方法，其特征在于：

发酵过程的连续时间非线性状态空间模型由一组微分方程组表示，对发酵过程的系统状态方程组和量测方程组进行离散化处理，得到的非线性系统采用非线性滤波算法进行状态估计；考虑离散化得到的发酵过程非线性系统方程如下：

x_k+1＝f(x_k,u_k)+w_k

其中x_k是发酵过程的系统状态向量，包括n个状态量，u_k是输入变量，f(·)是表征发酵过程模型的非线性函数，w_k是协方差为Q均值为0的系统噪声；

发酵过程中能够规律地、实时地获得的过程参量称为次要量测，规律地、实时地获得的过程参量为温度、pH值；需要人工离线分析或计算机自动控制取样分析的过程参量称为主要量测，需要人工离线分析或计算机自动控制取样分析的过程参量为菌体浓度、产物浓度和菌体比生长速率；

次要量测的方程如下：

y_k＝h⁰(x_k)+v_k，y_k为次要量测函数的应变量；

其中h⁰(·)是次要量测函数，v_k是均值为0方差为R₀的量测噪声；次要量测以固定的频率在每个时刻k进行采样，采样值能够实时到达数据处理单元；

量测滞后的主要量测方程为：

z_g＝h¹(x_s)+r_s，z_g为主要量测函数的应变量；

其中g＝s+d，h¹(·)是主要量测函数，r_s是均值为0方差为R₁的量测噪声；主要量测的目标是离线分析获得的生物量，s时刻是主要量测的采样时刻，分析结果延迟一段时间d并在主要时刻g到达数据处理单元，主要量测采样时刻s和延迟时间d在主要时刻g才已知；

该方法包括以下步骤：

步骤一 设定系统初始参数

进行初始参数设定，包括初始时刻的估计值初始时刻系统协方差协方差阵P_0|0；

步骤二 进行容积卡尔曼滤波时间更新并保存结果；

根据上一时刻的协方差阵P_k-1|k-1、状态估计值计算k时刻的容积点:

S_k-1|k-1＝(chol(P_k-1|k-1))^T

其中S_k-1|k-1是P_k-1|k-1的平方根，chol表示Cholesky分解，i＝1,2,3,…,m，m＝2n；[1]_i表示[1]的第i列；[1]表示一组由n维单位列向量全排列和变换符号构成的全对称点集；若则[1]代表：

计算经系统方程非线性传播后的容积点状态预测值及误差协方差P_k|k-1：

其中X_i,k-1|k-1是容积点，u_k-1是k-1时刻的输入变量，ω是各个容积点的权值，Q是系统噪声协方差；

为了在主要量测到达后利用主要量测信息在主要量测采样时刻进行量测更新，而主要量测采样时刻在主要时刻g才已知，因此需要保存仅利用次要量测进行容积卡尔曼滤波得到的每个时刻时间更新获得的状态预测值误差协方差P_k|k-1；

步骤三 利用次要量测进行容积卡尔曼滤波量测更新

计算容积点：

S_k|k-1＝(chol(P_k|k-1))^T

计算经次要量测方程传播后的容积点：

Y_i,k|k-1＝h⁰(X_i,k|k-1)

计算量测预测值：

计算自协方差阵和互协方差阵：

其中是容积点的权值，R₀是次要量测噪声的方差；

步骤四 利用次要量测信息进行量测修正并保存结果

W_y,k、P_y,k|k分别为利用次要量测信息计算得到的卡尔曼增益阵、状态预测值和误差协方差阵；

为了完成后续步骤求取主要量测含有的特有信息，而主要量测采样时刻在主要时刻才已知，需要保存每个时刻量测修正得到的P_y,k|k；

步骤五 进行下一时刻容积卡尔曼滤波或者利用主要量测进行量测更新和量测修正

如果主要量测信息还未到达，量测信息仅有次要量测，则返回步骤二，进行下一时刻的容积卡尔曼滤波；

如果主要量测信息到达，即主要量测信息已经可用，那么当前时刻k是主要时刻，且知道了主要量测采样时刻s，则从步骤二中已保存数据中调取s时刻的时间更新结果，并利用主要量测值，进行主要量测采样时刻的量测更新和量测修正；

计算容积点：

S_s|s-1＝(chol(P_s|s-1))^T

计算经主要量测方程传播后的容积点：

Z_i,s|s-1＝h¹(X_i,s|s-1)

计算量测预测值：

计算自协方差和互协方差：

其中是容积点的权值，R₁是主要量测噪声的方差；

利用主要量测进行量测修正：

W_z,s＝P_xz,s|s-1(P_zz,s|s-1)^-1

其中W_z,s|s、P_z,s|s分别为利用主要量测信息计算得到的主要量测采样时刻s的卡尔曼增益阵、状态预测值和误差协方差阵；

步骤六 将主要量测采样时刻s对应的结果递推到主要时刻g；

将s时刻的状态估计值和误差协方差P_z,s|s从s时刻进行时间更新:

S_z,s|s＝(chol(P_z,s|s))^T

对状态估计值误差协方差P_z,s+1|s再次进行时间更新，共进行延迟时间d次时间更新，最终递推得到主要时刻g的状态估计值误差协方差P_z,g|s；

至此，在主要时刻有两个状态估计值，分别是利用主要量测并递推得到的状态估计值和仅利用次要量测得到的状态估计值这两个状态值包含了部分共有信息，这部分共有信息是将s时刻仅利用次要量测的P_y,s|s经时间更新递推到主要时刻g的状态值误差协方差P_y,g|s；

步骤七 在主要时刻g进行协方差融合；

计算主要量测所含有的特有信息，其协方差阵和状态估计值如下：

把主要时刻g仅利用次要量测的容积卡尔曼滤波结果和主要量测含有的特有信息进行协方差融合：

P_g|g、即为主要时刻新的协方差阵和状态估计值，将此结果带回到步骤二即可进行下一时刻的状态估计。