CN103729562A - 一种基于重构识别分析的青霉素发酵过程故障监测方法 - Google Patents

Abstract

本发明一种基于重构识别分析的青霉素发酵过程故障监测方法，属于故障监测与诊断技术领域，本发明通过分析正常数据和故障数据的相关性，获得对正常状况影响最大的故障方向，根据故障贡献方向，将数据中的故障消除，即将工业状况模型恢复到正常状态，并识别此故障类型，并通过实际故障解决措施将实际工况恢复正常；本发明通过分析故障与正常状况的相关性，使得监测的准确性有了很大提高，为复杂工业过程的安全加强了安全保障，减少了损失，提高了产品的质量。

Description

一种基于重构识别分析的青霉素发酵过程故障监测方法

技术领域

本发明属于故障监测与诊断技术领域，具体涉及一种基于重构识别分析的青霉素发酵过程故障监测方法。 

背景技术

工业生产过程中存在潜在危险因素。这些危险因素在一定条件下，使得生产过程失败，甚至造成事故的发生、生产的扰乱、生活的不安全、环境的严重污染和经济的大量损失。这就对生产过程的监控和诊断提出了更高的要求。由于实际的工业生产过程中的不确定性及复杂的过程，一般包括非线性、时变性、变量耦合、时间相关性、多模态、多周期、大规模和间歇等特性，使它难以建立精确的过程模型。为了避免上述问题，基于数据驱动和知识驱动的复杂工业过程的故障诊断一直被视为一个热门话题，广泛用于工业过程。 

在过去十年中，多元统计分析技术已被广泛用于过程分析、监测和故障诊断。从海量数据中提取过程信息在各个领域中被发展，如主元分析法，独立元分析法和偏最小二乘法已广泛应用于工业生产过程的监测和诊断。当过程监测的统计量超出所需的控制限制，可以认为生产过程中出现异常的变化。当故障被检测到时，它需要人们可以迅速地诊断出异常的情况，并得到一个确定的原因，使它可以迅速被恢复正常。主元分析法是一种流行的，已被开发和广泛使用在许多工业生产过程中。主元分析法就是将高维数据投影到低维正交子空间的过程中，并保留关键过程的数据信息的方法。主元分析法空间被划分为主成分的子空间和剩余成分的子空间，结合T²（Hotelling’s T2，即霍特林T2统计量）和SPE（squared predicted error，平方预测误差）统计监测图和贡献图来监测正常工业状况是否发生异常。主元分析法虽然已经成功地使用在许多工业过程中，并对过程监控和诊断提供了一个很大的推动作用，但在生产过程中的应用中仍然还存在一些问题。虽然主元分析法用于故障检测已相当完善，但仍非常困难的识别那种影响多个传感器而产生的故障。因此，贡献图给出的是故障发生的结果，但不是故障发生的原因。 

发明内容

针对现有技术存在的不足，本发明提出一种基于重构识别分析的青霉素发酵过程故障监测方法，以达到不仅能够监测到故障的发生，还能识别出故障的类型，加强复杂工业过程的安全安全保障、减少损失和提高产品质量的目的。 

一种基于重构识别分析的青霉素发酵过程故障监测方法，包括以下步骤： 

步骤1、采集青霉素发酵过程中的多组历史数据，每组数据中包括：通风率、搅拌器功 率、基质进给速率、基质进给温度、生成的热量、溶解氧浓度、pH值、培养容积和二氧化碳浓度九个变量； 

步骤2、确定采集的多组历史数据中的正常数据组和故障数据组，并通过求平均值和标准差的方法对数据进行预处理，具体为： 

步骤2-1、将正常数据组和故障数据组存储至矩阵中，形成历史正常数据矩阵和历史故障数据矩阵； 

步骤2-2、计算历史正常数据组中每个变量的平均值和标准差； 

步骤2-3、将历史正常数据矩阵和历史故障数据矩阵的每个变量与计算获得的平均值做差，并与标准差相除，将获得的结果替换矩阵中原始位置的元素，进而完成对历史正常数据矩阵和历史故障数据矩阵的预处理； 

步骤3、采用主元分析方法分别对预处理后获得的历史正常数据矩阵和历史故障数据矩阵进行分解，获得历史正常数据得分矩阵、历史正常数据负载矩阵、历史故障数据得分矩阵和历史故障数据负载矩阵； 

步骤4、建立历史正常数据矩阵与历史故障数据矩阵之间的相关性，具体如下： 

步骤4-1、根据预处理后的历史正常数据矩阵、历史故障数据矩阵和采集历史数据的组数，建立历史正常数据矩阵和历史故障数据矩阵之间的协方差矩阵； 

步骤4-2、采用奇异值分解法对获得的协方差矩阵进行分解，获得具有历史正常数据矩阵和历史故障数据矩阵之间相关性的数据相关负载矩阵； 

步骤4-3、将历史故障数据矩阵与获得的数据相关负载矩阵相乘，计算获得数据相关得分矩阵； 

步骤4-4、计算所获数据相关得分矩阵的方差和历史故障数据得分矩阵的方差，根据上述两个方差的比值，获得数据相关得分矩阵与历史故障数据得分矩阵之间的比率向量； 

步骤4-5、将获得的比率向量中的每个比率值与设定的阈值相比较，提取其中大于阈值的比率值，并根据所提取的比率值所在比率向量的列，提取历史正常数据负载矩阵中对应列的数据，获得历史正常数据负载矩阵的一个子矩阵，即获得故障相关负载矩阵，并将该子矩阵作为数据重构的故障方向； 

步骤5、将获得的故障相关负载矩阵分别投影到主元子空间和残差空间，获得主元子空间内的故障方向，残差空间的故障方向； 

步骤6、实时采集被测青霉素发酵过程中的多组实时数据，并根据步骤2中的历史正常数据组中每个变量的平均值和标准差，对采集的数据进行预处理，构成被测数据矩阵； 

步骤7、根据被测数据矩阵、历史正常数据矩阵的转换矩阵和数据重构的故障方向，获 得在第一目标函数值最小时，SPE统计量重构中的故障量值；并根据历史正常数据负载矩阵、被测数据矩阵和数据重构的故障方向获得在第二目标函数值最小时，T²统计量重构中的故障量值； 

步骤8、根据获得的SPE统计量重构中的故障量值、数据重构的故障方向和被测数据矩阵获得SPE统计量重构中被测数据中的正常数据的部分；再根据历史正常数据矩阵的转换矩阵，获得SPE统计量重构中被测数据矩阵的SPE统计值，并通过χ²分布获得SPE统计量重构中的置信限； 

步骤9、根据T²统计量重构中的故障量值、数据重构的故障方向和被测数据矩阵，获得T²统计量重构中被测数据中的正常数据的部分；再根据历史正常数据负载矩阵和该矩阵对应的特征值，获得在T²统计量重构中被测数据矩阵的T²统计值；并通过F分布获得T²统计量重构中的置信限； 

步骤10、判断在T²统计量重构中被测数据矩阵的T²统计值是否小于其置信限，在SPE统计量重构中被测数据矩阵的SPE统计值是否小于其置信限，若同时小于，则该被测青霉素发酵过程中的故障类型为历史数据中存储的故障类型，并执行步骤12；否则该被测青霉素发酵过程中的故障类型不是历史数据中存储的故障类型，并执行步骤11； 

步骤11、返回执行步骤1，重新采集历史数据获取其他故障类型； 

步骤12、根据判断所得被测青霉素发酵过程的故障类型，采取解决措施，具体为： 

当通风率超出8.60L/h±0.05范围时，则调节通风口阀门的开度，控制通风大小，使通风率控制在设定范围内； 

当搅拌器功率超出30.0W±0.5范围时，则修改搅拌器的功率值，使搅拌器功率控制在设定范围内； 

当基质进给速率超出0.040L/h～0.045L/h范围时，则修改基质进给的量值，控制基质进给量，使基质进给速率控制在设定范围内； 

当基质进给温度超出296.0K±0.1K范围时，则修改进给基质的温度值，使基质进给温度控制在设定范围内； 

当生成的热量超出70kCal～80kCal范围时，则调节冷水和热水阀门的开度，控制冷水和热水流量大小，使生成的热量控制在设定范围内； 

当溶解氧浓度超出1.11g/L～1.13g/L范围时，则调节搅拌器的速率，使溶解氧的浓度控制在设定范围内； 

当pH值超出5.0±0.5范围时，则调节酸液和碱液阀门的开度，控制酸液和碱液流量大 小，使pH值控制在设定范围内； 

当培养容积超出用户设定范围时，则改变容积的大小，使其控制在设定范围内； 

当二氧化碳浓度超出2.0g/L～2.5g/L范围时，则调节调节搅拌器的速率，使二氧化碳的浓度控制在设定范围内。 

步骤3所述的采用主元分析方法分别对预处理后获得的历史正常数据矩阵和历史故障数据矩阵进行分解，公式如下： 

$\left(\begin{array}{c}T=\mathrm{XP}\\ X={\mathrm{TP}}^T+E=\hat{X}+E\end{array}\right.---\left(1\right)$

其中，X表示采集的历史数据矩阵，包括历史正常数据

和历史工况故障数据

N表示数据的组数，M表示变量数目，取值为9，x₁₁,x₁₂,...,x_1N表示采集的第1组到第N组历史正常数据，x₂₁,x₂₂,...,x_2N表示采集的第1组到第N组历史故障数据，

T表示数据主元分解后的得分矩阵，T(N×R)=[t₁,t₂,…,t_N]，包括历史正常数据得分矩阵T₁和历史故障数据得分矩阵T₂；P表示数据主元分解后的负载矩阵，P(M×R)=[p₁,p₂,…,p_M]，包括历史正常数据负载矩阵P₁和历史故障数据负载矩阵P₂，其中，R表示数据中包含的主元个数，R的取值满足 

λ_i是数据主元分解后的负载矩阵P中p_i列向量对应的特征值；E表示数据分解后的残差矩阵，是一个N行M列矩阵。 

步骤4-1所述的建立历史正常数据矩阵和历史故障数据矩阵之间的协方差矩阵，公式如下： 

$S=\frac{1}{N-1}{X}_1^T{X}_2---\left(2\right)$

其中，S表示历史正常数据矩阵和历史故障数据矩阵之间的协方差矩阵，是一个M行M列的矩阵； 

步骤4-2所述的采用奇异值分解法对获得的协方差矩阵进行分解，获得具有历史正常数据矩阵和历史故障数据矩阵之间相关性的数据相关负载矩阵，公式如下： 

S=P₃Λ₃P₃ ^T             （3） 

其中，Λ₃表示一个对角阵，包含幅值递减的非负特征值λ₁≥λ₂≥…≥λ_M≥0，P₃表示历 史正常数据矩阵和历史故障数据矩阵之间相关性的数据相关负载矩阵； 

步骤4-4所述的获得数据相关得分矩阵与历史故障数据得分矩阵之间的比率向量，公式如下： 

${\mathrm{Ratio}}_i=\frac{\mathrm{var}\left(T_3\right)}{\mathrm{var}\left(T_2\right)}(i=\mathrm{1,2},...,R)---\left(4\right)$

其中，Ratio_i表示第i个比率值，var(·)表示对得分矩阵求方差，T₃表示数据相关得分矩阵，T₃=X₂×P₃； 

步骤4-5所述的阈值取值范围为0.8～1.2。 

步骤5所述的获得主元子空间内的故障方向，残差空间的故障方向，公式如下： 

$\left(\begin{array}{c}{P}^{*}=Z=\hat{Z}+\tilde{Z}\\ \tilde{Z}=\tilde{C}Z=(I-P_1{P}_1^T)Z=(I-P_1{P}_1^T)P^{*}\\ \hat{Z}=\hat{C}Z=P_1{P}_1^{T}Z=P_1{P}_1^T{P}^{*}\end{array}\right.---\left(5\right)$

其中，I表示一个单位矩阵；P^*表示故障相关负载矩阵；Ζ表示数据重构的故障方向；

表示主元子空间内的故障方向，

表示残差空间的故障方向，表示历史正常数据映射到残差空间的转换矩阵，

表示历史正常数据映射到主元子空间的转换矩阵，  $\hat{C}=P_1{P}_1^T.$

步骤7所述的获得在第一目标函数值最小时，SPE统计量重构中的故障量值，并获得在第二目标函数值最小时，T²统计量重构中的故障量值，具体如下： 

SPE统计量重构中的故障量值，公式如下： 

$\begin{array}{c}{f}_{\mathrm{SPE}}=\arg \min {\left|\right|\widetilde{C^{}}{x^{*}}_{\mathrm{SPE}}\left|\right|}^2=\arg \min {\left|\right|\tilde{C}(x-{\mathrm{Zf}}_{\mathrm{SPE}})\left|\right|}^2\\ =\arg \min {\left|\right|\tilde{x}-\tilde{Z}{f}_{\mathrm{SPE}}\left|\right|}^2={\left({\tilde{Z}}^T\tilde{Z}\right)}^{-1}{\tilde{Z}}^{T}x\end{array}---\left(6\right)$

其中，argmin||.||表示目标函数绝对值最小时，未知数的取值；x^* _SPE表示SPE统计量重构中被测数据中的正常数据的部分，x^* _SPE=x-Ζf_SPE，x表示被测数据矩阵；

表示被测数据矩阵在残差空间的投影矩阵； 

将公式（5）中的故障方向Ζ在残差子空间映射的矩阵

代入公式（6）中获得故障量值f_SPE： 

$\begin{array}{c}{f}_{\mathrm{SPE}}=\left\{{\left[(I-P_1{P}_1^T)P^{*}\right]}^T\right[(I-P_1{P}_1^T)P^{*}{\left]\right\}}^{-1}[(I-P_1{P}_1^T)P^{*}{]}^{T}x\\ ={\left[P^{*T}(I-P_1{P}_1^T)P^{*}\right]}^{-1}{P}^{*T}(I-P_1{P}_1^T)x\end{array}---\left(7\right)$

在T²统计量的重构中的故障量值，公式如下： 

$\begin{array}{c}{f}_{T^2}=\arg \min {\left|\right|{\mathrm{\Λ}}^{-1/2}{P_1}^T{x^{*}}_{T^2}\left|\right|}^2=\arg \min \left|\right|{\mathrm{\Λ}}^{-1/2}{P_1}^T{x^{*}}_{T^2}{\left|\right|}^2\\ =\arg \min {\left|\right|{P_1}^T(x-Z{f}_{T^2})\left|\right|}^2={\left({\hat{Z}}^T\hat{Z}\right)}^{-1}{\hat{Z}}^{T}x\end{array}---\left(8\right)$

其中，

表示T²统计量重构中被测数据中的正常数据的部分，

P₁表示历史正常数据负载矩阵，Λ表示历史正常数据负载矩阵对应的特征值构成的对角矩阵； 

将公式（5）中的故障方向Ζ在主元子空间映射的矩阵

代入公式（8）中，获得故障量值

的求解公式： 

$\begin{array}{c}{f}_{T^2}{=\left\{{\left[P_1{P}_1^T{P}^{*}\right]}^T\right[P_1{P}_1^T{P}^{*}\left]\right\}}^{-1}[P_1{P}_1^T{P}^{*}{]}^{T}x\\ =[P^{*T}{P}_1{P}_1^T{P}^{*}{]}^{-1}{P}^{*T}{P}_1{P}_1^{T}x\end{array}---\left(9\right).$

步骤9所述的获得在T²统计量重构中被测数据矩阵的T²统计值，公式如下： 

T²=t^TΛ^-1t               （10） 

其中，t=x^*P₁，Λ=diag{λ₁,λ₂,…,λ_R}是对角矩阵，λ₁,λ₂,…,λ_R是历史正常数据负载矩阵对应的特征值； 

通过F分布获得T²统计量重构中的置信限，公式如下： 

$T_{R,N,\mathrm{\α}}^2~\frac{R(N-1)}{N-R}{F}_{R,N-R,\mathrm{\α}}---\left(11\right)$

其中，

表示T²统计量重构中的置信限，α表示置信度，R表示主元的个数，F_R,N-R,α表示置信度为α、第一自由度为R、第二自由度为N-R的F分布； 

步骤8所述的获得SPE统计量重构中x^* _SPE被测数据矩阵的SPE统计值，公式如下： 

$\mathrm{SPE}={\left|\right|\tilde{C}{x^{*}}_{\mathrm{SPE}}\left|\right|}^2={\left(\tilde{C}{x^{*}}_{\mathrm{SPE}}\right)}^T\tilde{C}{x^{*}}_{\mathrm{SPE}}={x^{*T}}_{\mathrm{SPE}}\tilde{C}{x^{*}}_{\mathrm{SPE}}---\left(12\right)$

其中，SPE表示SPE统计量重构中被测数据矩阵的SPE统计值； 

通过χ²分布获得SPE统计量重构中的置信限，公式如下： 

${\mathrm{SPE}}_{\mathrm{\α}}~{\mathrm{g\χ}}_h^2---\left(13\right)$

其中，g=b/2a，h=2a²/b，其中，a表示历史正常数据的平均值；b表示历史正常数据的方差；h表示置信度，

表示置信度为h的χ²分布。 

本发明优点： 

本发明一种基于重构识别分析的青霉素发酵过程故障监测方法，通过此方法分析正常数据和故障数据的相关性，获得对正常状况影响最大的故障方向，根据故障贡献方向，将数据中的故障消除，即将工业状况模型恢复到正常状态，并识别此故障类型，并通过实际故障解决措施将实际工况恢复正常；本发明通过分析故障与正常状况的相关性，使得监测的准确性有了很大提高，为复杂工业过程的安全加强了安全保障，减少了损失，提高了产品的质量。 

附图说明

图1为本发明一种实施例的青霉素发酵过程示意图； 

图2为本发明一种实施例的基于重构识别分析的青霉素发酵过程故障监测方法流程图； 

图3为本发明一种实施例的采集各个变量的数据示意图，其中，图（a）为通风率；图（b）为搅拌器功率；图（c）为基质给进速率；图（d）为基质给进温度；图（e）为生成的热量；图（f）为溶解氧的浓度；图（g）为ph值；图（h）为二氧化碳浓度；图（i）为培养容积； 

图4为本发明一种实施例的青霉素发酵过程的通风率加入50%的阶跃跳变时T²和SPE统计量示意图，其中，图（a）为采用主元分析算法得出的青霉素发酵过程的T²统计量；图（b）为采用主元分析算法得出的青霉素发酵过程的SPE统计量；图（c）为采用故障相关重构算法得出的青霉素发酵过程的T²统计量；图（d）为采用故障相关重构算法得出的青霉素发酵过程的SPE统计量； 

图5为本发明一种实施例的青霉素发酵过程的搅拌器功率加入50%的阶跃跳变时T²和SPE统计量示意图，其中，图（a）为采用主元分析算法得出的青霉素发酵过程的T²统计量；图（b）为采用主元分析算法得出的青霉素发酵过程的SPE统计量；图（c）为采用故障相关重构算法得出的青霉素发酵过程的T²统计量；图（d）为采用故障相关重构算法得出的青霉素发酵过程的SPE统计量； 

图6为本发明一种实施例的青霉素发酵过程的基质进给速率加入50%的阶跃跳变时T²和SPE统计量示意图，其中，图（a）为采用主元分析算法得出的青霉素发酵过程的T²统计量；图（b）为采用主元分析算法得出的青霉素发酵过程的SPE统计量；图（c）为采用故障相关 重构算法得出的青霉素发酵过程的T²统计量；图（d）为采用故障相关重构算法得出的青霉素发酵过程的SPE统计量； 

图7为本发明一种实施例的青霉素发酵过程的通风率加入2%的斜坡跳变时T²和SPE统计量示意图，其中，图（a）为采用主元分析算法得出的青霉素发酵过程的T²统计量；图（b）为采用主元分析算法得出的青霉素发酵过程的SPE统计量；图（c）为采用故障相关重构算法得出的青霉素发酵过程的T²统计量；图（d）为采用故障相关重构算法得出的青霉素发酵过程的SPE统计量； 

图8为本发明一种实施例的青霉素发酵过程的基质进给速率加入2%的斜坡跳变时T²和SPE统计量示意图，其中，图（a）为采用主元分析算法得出的青霉素发酵过程的T²统计量；图（b）为采用主元分析算法得出的青霉素发酵过程的SPE统计量；图（c）为采用故障相关重构算法得出的青霉素发酵过程的T²统计量；图（d）为采用故障相关重构算法得出的青霉素发酵过程的SPE统计量； 

图9为本发明一种实施例的青霉素发酵过程的搅拌器功率加入2%的斜坡跳变时T²和SPE统计量示意图，其中，图（a）为采用主元分析算法得出的青霉素发酵过程的T²统计量；图（b）为采用主元分析算法得出的青霉素发酵过程的SPE统计量；图（c）为采用故障相关重构算法得出的青霉素发酵过程的T²统计量；图（d）为采用故障相关重构算法得出的青霉素发酵过程的SPE统计量。 

具体实施方式

下面结合附图对本发明一种实施例做进一步说明。 

本发明实施例中，青霉素发酵过程是青霉素产生菌在合适的培养基、pH值、温度、空气流量、搅拌等发酵条件下进行生长和合成抗生素的代谢活动。图1为青霉素生产发酵过程示意图，其中，被控变量包括发酵罐的pH值和温度，它们分别通过操纵变量：酸、碱流量和冷、热水流量将其控制在一定的值，主要采用控制器FC控制酸、碱流量和冷、热水的阀门的开度，来调节pH值和温度的。在青霉素发酵过程中温度和pH值采用闭环控制，而补料采用开环定值控制，青霉素发酵的每个批次持续时间是400小时，包含约45小时的预培养阶段和约355小时的间歇补料阶段。 

一种基于重构识别分析的青霉素发酵过程故障监测方法，方法流程图如图2所示，包括 以下步骤： 

步骤1、采集青霉素发酵过程中的多组历史数据，每组数据中包括：通风率、搅拌器功率、基质进给速率、基质进给温度、生成的热量、溶解氧浓度、pH值、培养容积和二氧化碳浓度九个变量； 

本发明实施例中，温度和pH值采用闭环控制，而补料采用开环定值控制；对于青霉素发酵过程采集两批数据，采用0.5h的采样间隔采集历史数据，作为训练数据用来建立青霉素发酵过程初始的模型，由400个观测数构成，其中在此数据建模和检测中取99%的置信限，每组样本数据包含9个变量，各个变量如图3中图（a）～图（i）的数据示意图，建立青霉素发酵过程初始的监测模型。本发明实施例中，随机选取了10组标准数据如表1所示，本发明实施例中以400组历史数据为例对以下内容进行说明。 

表1.青霉素发酵建模中的十组数据 

步骤2、确定采集的400组历史数据中的正常数据组和故障数据组，并通过求平均值和标准差的方法对数据进行预处理，具体为： 

步骤2-1、将正常数据组和故障数据组存储至矩阵中，形成历史正常数据矩阵和历史故障数据矩阵； 

步骤2-2、计算历史正常数据组中每个变量的平均值和标准差； 

步骤2-3、将历史正常数据矩阵和历史故障数据矩阵的每个变量与计算获得的平均值做差，并与标准差相除，将获得的结果替换矩阵中原始位置的元素，进而完成对历史正常数据 矩阵和历史故障数据矩阵的预处理； 

步骤3、采用主元分析方法分别对预处理后获得的历史正常数据矩阵  $X_1(400\×9)={[X_{11}^T,x_{12}^T,...,x_{1400}^T]}^T$ 和历史故障数据矩阵 $X_2(400\×9)={[x_{21}^T,x_{22}^T,...,x_{2400}^T]}^T$ 进行分解，获得历史正常数据得分矩阵T₁、历史正常数据负载矩阵P₁、历史故障数据得分矩阵T₂和历史故障数据负载矩阵P₂； 

本发明实施例中，将历史正常数据矩阵

和历史故障数据矩阵

代入公式（1）中进行分解，公式如下： 

$\left(\begin{array}{c}T=\mathrm{XP}\\ X={\mathrm{TP}}^T+E=\hat{X}+E\end{array}\right.---\left(1\right)$

其中，X表示采集的历史数据矩阵，包括历史正常数据和历史工况故障数据，x₁₁,x₁₂,...,x₁₄₀₀表示采集的第1组到第400组历史正常数据，x₂₁,x₂₂,...,x₂₄₀₀表示采集的第1组到第400组历史故障数据

T表示数据主元分解后的得分矩阵，T(400×R)=[t₁,t₂,…,t₄₀₀]，包括历史正常数据得分矩阵T₁和历史故障数据得分矩阵T₂；P表示数据主元分解后的负载矩阵，P(9×R)=[p₁,p₂,…,p₉]，包括历史正常数据负载矩阵P₁和历史故障数据负载矩阵P₂，其中，R表示数据中包含的主元个数，R的取值满足 

λ_i是数据主元分解后的负载矩阵P中p_i列向量对应的特征值；E表示数据分解后的残差矩阵，是一个400行9列矩阵。 

步骤4、建立历史正常数据矩阵X₁与历史故障数据矩阵X₂之间的相关性，具体如下： 

步骤4-1、根据预处理后的历史正常数据矩阵X₁、历史故障数据矩阵X₂和采集历史数据的组数400，建立历史正常数据矩阵和历史故障数据矩阵之间的协方差矩阵S； 

公式如下： 

$S=\frac{1}{N-1}{X}_1^T{X}_2---\left(2\right)$

其中，S表示历史正常数据矩阵和历史故障数据矩阵之间的协方差矩阵，是一个9行9列的矩阵； 

步骤4-2、采用奇异值分解法对获得的协方差矩阵S进行分解，获得具有历史正常数据矩 阵和历史故障数据矩阵之间相关性的数据相关负载矩阵； 

本发明实施例中，将获得的协方差矩阵S代入公式（3）中进行分解，公式如下： 

S=P₃Λ₃P₃ ^T             （3） 

其中，Λ₃表示一个对角阵，包含幅值递减的非负特征值λ₁≥λ₂≥…≥λ₉≥0（半正定矩阵），P₃表示历史正常数据矩阵X₁和历史故障数据矩阵X₂之间相关性的数据相关负载矩阵； 

步骤4-3、将历史故障数据矩阵X₂与获得的数据相关负载矩阵P₃相乘，计算获得数据相关得分矩阵T₃，即T₃=X₂×P₃； 

步骤4-4、计算所获数据相关得分矩阵T₃的方差和历史故障数据得分矩阵T₂的方差，根据上述两个方差的比值，获得数据相关得分矩阵T₃与历史故障数据得分矩阵T₂之间的比率向量； 

本发明实施例中，利用历史故障数据的得分矩阵T₂和数据相关得分矩阵T₃之间的比率Ratio_i(i=1,2,...,R)来截取到对正常状况影响较大的故障方向，得到故障相关负载矩阵P^*(M×R)。其中j∈[1,2,…,R]是主元故障情况i的子集。 

公式如下： 

${\mathrm{Ratio}}_i=\frac{\mathrm{var}\left(T_3\right)}{\mathrm{var}\left(T_2\right)}(i=\mathrm{1,2},...,R)---\left(4\right)$

其中，Ratio_i表示第i个比率值，var(·)表示对得分矩阵求方差，T₃表示数据相关得分矩阵； 

步骤4-5、将获得的比率向量中的每个比率值与设定的阈值（取值范围为0.8～1.2，本发明实施例中具体取值为1）相比较，，提取其中大于阈值的比率值，并根据所提取的比率值所在比率向量的列，提取历史正常数据负载矩阵中对应列的数据，获得历史正常数据负载矩阵的一个子矩阵，即获得故障相关负载矩阵P^*，并将该子矩阵作为数据重构的故障方向Ζ； 

本发明实施例中，如果比率值大于阙值，这就意味着故障对于正常状况的影响更显著，并能引起T²和SPE统计量的显著增加，偏离正常状况下的统计量控制线；如果比率值小于阙值，那故障对于正常状况的影响将被考虑很小或忽略不记。即当比较值满足Ratio_j>α时，可 以截取到对正常状况影响较大的故障方向，j∈[1,2,…,R]是主元故障情况i的子集。 

步骤5、将获得的故障相关负载矩阵P^*分别投影到主元子空间和残差空间，获得主元子空间内的故障方向残差空间的故障方向

公式如下： 

$\left(\begin{array}{c}{P}^{*}=Z=\hat{Z}+\tilde{Z}\\ \tilde{Z}=\tilde{C}Z=(I-P_1{P}_1^T)Z=(I-P_1{P}_1^T)P^{*}\\ \hat{Z}=\hat{C}Z=P_1{P}_1^{T}Z=P_1{P}_1^T{P}^{*}\end{array}\right.---\left(5\right)$

本发明实施例中，将故障相关负载矩阵P^*替代数据重构的故障方向Ζ，在公式（5）中进行计算，其中，I表示一个单位矩阵；P^*表示故障相关负载矩阵；Ζ表示数据重构的故障方向；

表示主元子空间内的故障方向，

表示残差空间的故障方向，

表示历史正常数据映射到残差空间的转换矩阵，

表示历史正常数据映射到主元子空间的转换矩阵， $\hat{C}=P_1{P}_1^T.$

步骤6、实时采集被测青霉素发酵过程中的多组实时数据，并根据步骤2中的历史正常数据组中每个变量的平均值和标准差，对采集的数据进行预处理，构成被测数据矩阵； 

本发明实施例中，对实时采集的被测数据随机抽取10组如表2所示： 

表2.青霉素发酵测试数据中的十组数据 

步骤7、根据被测数据矩阵x、历史正常数据矩阵的转换矩阵

和数据重构的故障方向Z，获得在第一目标函数值最小时，SPE统计量重构中的故障量值f_SPE；并根据历史正常数据负载矩阵P₁、被测数据矩阵x和数据重构的故障方向Ζ获得在第二目标函数值最小时，T²统计量重构中的故障量值

SPE统计量重构中的故障量值，公式如下： 

$\begin{array}{c}{f}_{\mathrm{SPE}}=\arg \min {\left|\right|\widetilde{C^{}}{x^{*}}_{\mathrm{SPE}}\left|\right|}^2=\arg \min {\left|\right|\tilde{C}(x-{\mathrm{Zf}}_{\mathrm{SPE}})\left|\right|}^2\\ =\arg \min {\left|\right|\tilde{x}-\tilde{Z}{f}_{\mathrm{SPE}}\left|\right|}^2={\left({\tilde{Z}}^T\tilde{Z}\right)}^{-1}{\tilde{Z}}^{T}x\end{array}---\left(6\right)$

本发明实施例中，根据故障相关方向对实时工况故障数据重构识别。在重构算法中，能够从工况数据中区分故障，就要得到精确地故障相关方向。所有可能的故障数据集合为F₁，而实际故障数据F₂则是F₁的子集，根据重构算法，获得x^* _SPE：x^* _SPE表示SPE统计量重构中被测数据中的正常数据的部分，x^* _SPE=x-Ζf_SPE，x表示被测数据矩阵；重构就是将正常部分x^*从扰动的数据x中沿着给定的故障方向Ζ中重构出来，

是一个正交矩阵； 

本发明实施例中，在SPE统计量的重构中，通过求得正常部分x^*的最好估计值来计算量值f_SPE，即用其到残差子空间(S_r)的最短距离这个约束条件求得。 

将公式（5）中的故障方向Ζ在残差子空间映射的矩阵

代入公式（6）中获得故障量值f_SPE： 

$\begin{array}{c}{f}_{\mathrm{SPE}}=\left\{{\left[(I-P_1{P}_1^T)P^{*}\right]}^T\right[(I-P_1{P}_1^T)P^{*}{\left]\right\}}^{-1}[(I-P_1{P}_1^T)P^{*}{]}^{T}x\\ ={\left[P^{*T}(I-P_1{P}_1^T)P^{*}\right]}^{-1}{P}^{*T}(I-P_1{P}_1^T)x\end{array}---\left(7\right)$

在T²统计量的重构中的故障量值

公式如下： 

$\begin{array}{c}{f}_{T^2}=\arg \min {\left|\right|{\mathrm{\Λ}}^{-1/2}{P_1}^T{x^{*}}_{T^2}\left|\right|}^2=\arg \min \left|\right|{\mathrm{\Λ}}^{-1/2}{P_1}^T{x^{*}}_{T^2}{\left|\right|}^2\\ =\arg \min {\left|\right|{P_1}^T(x-Z{f}_{T^2})\left|\right|}^2={\left({\hat{Z}}^T\hat{Z}\right)}^{-1}{\hat{Z}}^{T}x\end{array}---\left(8\right)$

其中，

表示T²统计量重构中被测数据中的正常数据的部分，

P₁表示历史正常数据负载矩阵，Λ表示历史正常数据负载矩阵对应的特征值构成的对角矩阵； 

将公式（5）中的故障方向Ζ在主元子空间映射的矩阵

代入公式（8）中，获得故障量 值

的求解公式： 

$\begin{array}{c}{f}_{T^2}{=\left\{{\left[P_1{P}_1^T{P}^{*}\right]}^T\right[P_1{P}_1^T{P}^{*}\left]\right\}}^{-1}[P_1{P}_1^T{P}^{*}{]}^{T}x\\ =[P^{*T}{P}_1{P}_1^T{P}^{*}{]}^{-1}{P}^{*T}{P}_1{P}_1^{T}x\end{array}---\left(9\right)$

步骤8、根据获得的SPE统计量重构中的故障量值f_SPE、数据重构的故障方向Ζ和被测数据矩阵x获得SPE统计量重构中被测数据中的正常数据的部分x^* _SPE；再根据历史正常数据矩阵的转换矩阵

获得SPE统计量重构中被测数据矩阵的SPE统计值，并通过χ²分布获得SPE统计量重构中的置信限； 

本发明实施例中，获得SPE统计值，公式如下： 

$\mathrm{SPE}={\left|\right|\tilde{C}{x^{*}}_{\mathrm{SPE}}\left|\right|}^2={\left(\tilde{C}{x^{*}}_{\mathrm{SPE}}\right)}^T\tilde{C}{x^{*}}_{\mathrm{SPE}}={x^{*T}}_{\mathrm{SPE}}\tilde{C}{x^{*}}_{\mathrm{SPE}}---\left(12\right)$

其中，SPE表示SPE统计量重构中被测数据矩阵的SPE统计值； 

通过χ²分布获得SPE统计量重构中的置信限,χ²分布是服从标准正态分布的相互独立的随机变量的平方和构成的一组新的随机变量的分布，公式如下： 

${\mathrm{SPE}}_{\mathrm{\α}}~{\mathrm{g\χ}}_h^2---\left(13\right)$

其中，g=b/2a，h=2a²/b，其中，a表示历史正常数据的平均值；b表示历史正常数据的方差；h表示置信度，

表示置信度为h的χ²分布。 

步骤9、根据T²统计量重构中的故障量值

数据重构的故障方向Ζ和被测数据矩阵x，获得T²统计量重构中被测数据中的正常数据的部分

再根据历史正常数据负载矩阵P₁和该矩阵对应的特征值λ₁,λ₂,…,λ_R，获得在T²统计量重构中被测数据矩阵的T²统计值；并通过F分布获得T²统计量重构中的置信限； 

获得在T²统计量重构中被测数据矩阵的T²统计值，公式如下： 

T²=t^TΛ^-1t                （10） 

其中，t=x^*P₁，Λ=diag{λ₁,λ₂,…,λ_R}是对角矩阵，λ₁,λ₂,…,λ_R是历史正常数据负载矩阵对应的特征值； 

通过F分布获得T²统计量重构中的置信限，F分布是两个独立的随机变量的卡方分布被各自的自由度除以后的比率这一统计量的分布，公式如下： 

$T_{R,N,\mathrm{\α}}^2~\frac{R(N-1)}{N-R}{F}_{R,N-R,\mathrm{\α}}---\left(11\right)$

其中，

表示T²统计量重构中的置信限，α表示置信度，R表示主元的个数，F_R,N-R,α表示置信度为α、第一自由度为R、第二自由度为N-R的F分布； 

步骤10、判断在T²统计量重构中被测数据矩阵的T²统计值是否小于其置信限，在SPE统计量重构中被测数据矩阵的SPE统计值是否小于其置信限，若同时小于，则该被测青霉素发酵过程中的故障类型为历史数据中存储的故障类型，并执行步骤12；否则该被测青霉素发酵过程中的故障类型不是历史数据中存储的故障类型，并执行步骤11； 

本发明实施例中，T²统计值和SPE统计值分别形成两条曲线，置信限为一条控制线，判断统计量的曲线是否超出控制线，当对某一类型实时工况故障数据用故障相关的重构算法处理后，观测数据的T²和SPE统计量如果与正常数据的统计量浮动相同，并且没有超出统计量规定的控制线时，则此算法识别出了此类型故障；反之如果观测数据的T²和SPE统计量与正常数据不相符，且超出统计量规定的控制线，则认为算法识别不出此类型故障。 

步骤11、返回执行步骤1，重新采集历史数据获取其他故障类型； 

步骤12、根据判断所得被测青霉素发酵过程的故障类型，采取解决措施，具体为： 

当通风率超出8.60L/h±0.05范围时，则调节通风口阀门的开度，控制通风大小，使通风率控制在设定范围内； 

当搅拌器功率超出30.0W±0.5范围时，则修改搅拌器的功率值，使搅拌器功率控制在设定范围内； 

当基质进给速率超出0.040L/h～0.045L/h范围时，则修改基质进给的量值，控制基质进给量，使基质进给速率控制在设定范围内； 

当基质进给温度超出296.0K±0.1K范围时，则修改进给基质的温度值，使基质进给温度控制在设定范围内； 

当生成的热量超出70kCal～80kCal范围时，则调节冷水和热水阀门的开度，控制冷水和热水流量大小，使生成的热量控制在设定范围内； 

当溶解氧浓度超出1.11g/L～1.13g/L范围时，则调节搅拌器的速率，使溶解氧的浓度控制在设定范围内； 

当pH值超出5.0±0.5范围时，则调节酸液和碱液阀门的开度，控制酸液和碱液流量大小，使pH值控制在设定范围内； 

当培养容积超出用户设定范围时，则改变容积的大小，使其控制在设定范围内； 

当二氧化碳浓度超出2.0g/L～2.5g/L范围时，则调节调节搅拌器的速率，使二氧化碳的浓度控制在设定范围内。 

本发明实施例中，在青霉素发酵过程中，青霉素发酵过程新采样x是工况时，采用重构分析方法更新步骤1素发酵过程的初始模型，并利用重构分析方法的求得T²统计和SPE统计中重构的正常部分

最后通过统计变量T²和SPE进行故障监测。本发明实施例中，将本方法与主元分析法进行比较。六个不正常的批次数据用于在线监测。前三个不正常的批次是分别在通风率、搅拌器功率、基质进给速率三个变量上增加了50%的阶跃跳变。如图4所示，故障从第100个样本开始持续到第400个样本。如图4中图（a）和图4中图（b）所示，主元分析法能够监测到故障并发出警报。如图4中图（c）和图4中图（d）所示，本发明实施例中所提出的基于重构识别分析方法不仅能够监测出故障，而且能够识别出故障类型，并将故障消除，恢复为正常数据。图5中图（a）、图（b）、图（c）、图（d）和图6中图（a）、图（b）、图（c）、图（d）所示的事例与图4类似。并在后三个不正常批次的通风率、搅拌器功率、基质进给速率三个变量上增加了2%的斜坡跳变。如图7所示，故障从第100个样本开始持续到第400个样本。如图7中图（a）和图7中图（b）所示，主元分析法能够监测到故障并发出警报。如图7中图（c）和图7中图（d）所示，本发明实施例中所提出的基于重构识别分析方法不仅能够监测出故障，而且能够识别出故障类型，并将故障消除，恢复为正常数据。图8中图（a）、图（b）、图（c）、图（d）和图9中图（a）、图（b）、图（c）、图（d）所示的事例与图4类似。 

由上述对比分析可得，本发明基于重构识别分析的青霉素发酵过程故障监测方法有效地找到对正常状况影响最大的故障方向，本发明所提出的方法比主元分析法更具优势。不仅能够监测到故障的发生，还能识别出故障的类型，为复杂工业过程的安全加强了安全保障，减少了损失，提高了产品的质量。 

Claims (7)

1.一种基于重构识别分析的青霉素发酵过程故障监测方法，其特征在于，包括以下步骤：

步骤1、采集青霉素发酵过程中的多组历史数据，每组数据中包括：通风率、搅拌器功率、基质进给速率、基质进给温度、生成的热量、溶解氧浓度、pH值、培养容积和二氧化碳浓度九个变量；

步骤2、确定采集的多组历史数据中的正常数据组和故障数据组，并通过求平均值和标准差的方法对数据进行预处理，具体为：

步骤2-1、将正常数据组和故障数据组存储至矩阵中，形成历史正常数据矩阵和历史故障数据矩阵；

步骤2-2、计算历史正常数据组中每个变量的平均值和标准差；

步骤2-3、将历史正常数据矩阵和历史故障数据矩阵的每个变量与计算获得的平均值做差，并与标准差相除，将获得的结果替换矩阵中原始位置的元素，进而完成对历史正常数据矩阵和历史故障数据矩阵的预处理；

步骤3、采用主元分析方法分别对预处理后获得的历史正常数据矩阵和历史故障数据矩阵进行分解，获得历史正常数据得分矩阵、历史正常数据负载矩阵、历史故障数据得分矩阵和历史故障数据负载矩阵；

步骤4、建立历史正常数据矩阵与历史故障数据矩阵之间的相关性，具体如下：

步骤4-1、根据预处理后的历史正常数据矩阵、历史故障数据矩阵和采集历史数据的组数，建立历史正常数据矩阵和历史故障数据矩阵之间的协方差矩阵；

步骤4-2、采用奇异值分解法对获得的协方差矩阵进行分解，获得具有历史正常数据矩阵和历史故障数据矩阵之间相关性的数据相关负载矩阵；

步骤4-3、将历史故障数据矩阵与获得的数据相关负载矩阵相乘，计算获得数据相关得分矩阵；

步骤4-4、计算所获数据相关得分矩阵的方差和历史故障数据得分矩阵的方差，根据上述两个方差的比值，获得数据相关得分矩阵与历史故障数据得分矩阵之间的比率向量；

步骤4-5、将获得的比率向量中的每个比率值与设定的阈值相比较，提取其中大于阈值的比率值，并根据所提取的比率值所在比率向量的列，提取历史正常数据负载矩阵中对应列的数据，获得历史正常数据负载矩阵的一个子矩阵，即获得故障相关负载矩阵，并将该子矩阵作为数据重构的故障方向；

步骤5、将获得的故障相关负载矩阵分别投影到主元子空间和残差空间，获得主元子空间内的故障方向，残差空间的故障方向；

步骤6、实时采集被测青霉素发酵过程中的多组实时数据，并根据步骤2中的历史正常数据组中每个变量的平均值和标准差，对采集的数据进行预处理，构成被测数据矩阵；

步骤7、根据被测数据矩阵、历史正常数据矩阵的转换矩阵和数据重构的故障方向，获得在第一目标函数值最小时，SPE统计量重构中的故障量值；并根据历史正常数据负载矩阵、被测数据矩阵和数据重构的故障方向获得在第二目标函数值最小时，T²统计量重构中的故障量值；

步骤8、根据获得的SPE统计量重构中的故障量值、数据重构的故障方向和被测数据矩阵获得SPE统计量重构中被测数据中的正常数据的部分；再根据历史正常数据矩阵的转换矩阵，获得SPE统计量重构中被测数据矩阵的SPE统计值，并通过χ²分布获得SPE统计量重构中的置信限；

步骤9、根据T²统计量重构中的故障量值、数据重构的故障方向和被测数据矩阵，获得T²统计量重构中被测数据中的正常数据的部分；再根据历史正常数据负载矩阵和该矩阵对应的特征值，获得在T²统计量重构中被测数据矩阵的T²统计值；并通过F分布获得T²统计量重构中的置信限；

步骤10、判断在T²统计量重构中被测数据矩阵的T²统计值是否小于其置信限，在SPE统计量重构中被测数据矩阵的SPE统计值是否小于其置信限，若同时小于，则该被测青霉素发酵过程中的故障类型为历史数据中存储的故障类型，并执行步骤12；否则该被测青霉素发酵过程中的故障类型不是历史数据中存储的故障类型，并执行步骤11；

步骤11、返回执行步骤1，重新采集历史数据获取其他故障类型；

步骤12、根据判断所得被测青霉素发酵过程的故障类型，采取解决措施，具体为：

当通风率超出8.60L/h±0.05范围时，则调节通风口阀门的开度，控制通风大小，使通风率控制在设定范围内；

当搅拌器功率超出30.0W±0.5范围时，则修改搅拌器的功率值，使搅拌器功率控制在设定范围内；

当基质进给速率超出0.040L/h～0.045L/h范围时，则修改基质进给的量值，控制基质进给量，使基质进给速率控制在设定范围内；

当基质进给温度超出296.0K±0.1K范围时，则修改进给基质的温度值，使基质进给温度控制在设定范围内；

当生成的热量超出70kCal～80kCal范围时，则调节冷水和热水阀门的开度，控制冷水和热水流量大小，使生成的热量控制在设定范围内；

当溶解氧浓度超出1.11g/L～1.13g/L范围时，则调节搅拌器的速率，使溶解氧的浓度控制在设定范围内；

当pH值超出5.0±0.5范围时，则调节酸液和碱液阀门的开度，控制酸液和碱液流量大小，使pH值控制在设定范围内；

当培养容积超出用户设定范围时，则改变容积的大小，使其控制在设定范围内；

当二氧化碳浓度超出2.0g/L～2.5g/L范围时，则调节调节搅拌器的速率，使二氧化碳的浓度控制在设定范围内。

2.根据权利要求1所述的基于重构识别分析的青霉素发酵过程故障监测方法，其特征在于，步骤3所述的采用主元分析方法分别对预处理后获得的历史正常数据矩阵和历史故障数据矩阵进行分解，公式如下：

$\left(\begin{array}{c}T=\mathrm{XP}\\ X={\mathrm{TP}}^T+E=\hat{X}+E\end{array}\right.---\left(1\right)$

其中，X表示采集的历史数据矩阵，包括历史正常数据

和历史工况故障数据

N表示数据的组数，M表示变量数目，取值为9，x₁₁,x₁₂,...,x_1N表示采集的第1组到第N组历史正常数据，x₂₁,x₂₂,...,x_2N表示采集的第1组到第N组历史故障数据，

T表示数据主元分解后的得分矩阵，T(N×R)=[t₁,t₂,…,t_N]，包括历史正常数据得分矩阵T₁和历史故障数据得分矩阵T₂；P表示数据主元分解后的负载矩阵，P(M×R)=[p₁,p₂,…,p_M]，包括历史正常数据负载矩阵P₁和历史故障数据负载矩阵P₂，其中，R表示数据中包含的主元个数，R的取值满足λ_i是数据主元分解后的负载矩阵P中p_i列向量对应的特征值；E表示数据分解后的残差矩阵，是一个N行M列矩阵。

3.根据权利要求1所述的基于重构识别分析的青霉素发酵过程故障监测方法，其特征在于，步骤4-1所述的建立历史正常数据矩阵和历史故障数据矩阵之间的协方差矩阵，公式如下：

$S=\frac{1}{N-1}{X}_1^T{X}_2---\left(2\right)$

其中，S表示历史正常数据矩阵和历史故障数据矩阵之间的协方差矩阵，是一个M行M列的矩阵；

步骤4-2所述的采用奇异值分解法对获得的协方差矩阵进行分解，获得具有历史正常数据矩阵和历史故障数据矩阵之间相关性的数据相关负载矩阵，公式如下：

S=P₃Λ₃P₃ ^T                （3）

其中，Λ₃表示一个对角阵，包含幅值递减的非负特征值λ₁≥λ₂≥…≥λ_M≥0，P₃表示历史正常数据矩阵和历史故障数据矩阵之间相关性的数据相关负载矩阵；

步骤4-4所述的获得数据相关得分矩阵与历史故障数据得分矩阵之间的比率向量，公式如下：

${\mathrm{Ratio}}_i=\frac{\mathrm{var}\left(T_3\right)}{\mathrm{var}\left(T_2\right)}(i=\mathrm{1,2},...,R)---\left(4\right)$

其中，Ratio_i表示第i个比率值，var(·)表示对得分矩阵求方差，T₃表示数据相关得分矩阵，T₃=X₂×P₃；

步骤4-5所述的阈值取值范围为0.8～1.2。

4.根据权利要求1所述的基于重构识别分析的青霉素发酵过程故障监测方法，其特征在于，步骤5所述的获得主元子空间内的故障方向，残差空间的故障方向，公式如下：

$\left(\begin{array}{c}{P}^{*}=Z=\hat{Z}+\tilde{Z}\\ \tilde{Z}=\tilde{C}Z=(I-P_1{P}_1^T)Z=(I-P_1{P}_1^T)P^{*}\\ \hat{Z}=\hat{C}Z=P_1{P}_1^{T}Z=P_1{P}_1^T{P}^{*}\end{array}\right.---\left(5\right)$

其中，I表示一个单位矩阵；P^*表示故障相关负载矩阵；Ζ表示数据重构的故障方向；

表示主元子空间内的故障方向，

表示残差空间的故障方向，

表示历史正常数据映射到残差空间的转换矩阵，

表示历史正常数据映射到主元子空间的转换矩阵， $\hat{C}=P_1{P}_1^T.$

5.根据权利要求1所述的基于重构识别分析的青霉素发酵过程故障监测方法，其特征在于，步骤7所述的获得在第一目标函数值最小时，SPE统计量重构中的故障量值，并获得在第二目标函数值最小时，T²统计量重构中的故障量值，具体如下：

SPE统计量重构中的故障量值，公式如下：

$\begin{array}{c}{f}_{\mathrm{SPE}}=\arg \min {\left|\right|\widetilde{C^{}}{x^{*}}_{\mathrm{SPE}}\left|\right|}^2=\arg \min {\left|\right|\tilde{C}(x-{\mathrm{Zf}}_{\mathrm{SPE}})\left|\right|}^2\\ =\arg \min {\left|\right|\tilde{x}-\tilde{Z}{f}_{\mathrm{SPE}}\left|\right|}^2={\left({\tilde{Z}}^T\tilde{Z}\right)}^{-1}{\tilde{Z}}^{T}x\end{array}---\left(6\right)$

其中，argmin||.||表示目标函数绝对值最小时，未知数的取值；x^* _SPE表示SPE统计量重构中被测数据中的正常数据的部分，x^* _SPE=x-Ζf_SPE，x表示被测数据矩阵；

表示被测数据矩阵在残差空间的投影矩阵；

将公式（5）中的故障方向Ζ在残差子空间映射的矩阵

代入公式（6）中获得故障量值f_SPE：

$\begin{array}{c}{f}_{\mathrm{SPE}}=\left\{{\left[(I-P_1{P}_1^T)P^{*}\right]}^T\right[(I-P_1{P}_1^T)P^{*}{\left]\right\}}^{-1}[(I-P_1{P}_1^T)P^{*}{]}^{T}x\\ ={\left[P^{*T}(I-P_1{P}_1^T)P^{*}\right]}^{-1}{P}^{*T}(I-P_1{P}_1^T)x\end{array}---\left(7\right)$

在T²统计量的重构中的故障量值，公式如下：

$\begin{array}{c}{f}_{T^2}=\arg \min {\left|\right|{\mathrm{\Λ}}^{-1/2}{P_1}^T{x^{*}}_{T^2}\left|\right|}^2=\arg \min \left|\right|{\mathrm{\Λ}}^{-1/2}{P_1}^T{x^{*}}_{T^2}{\left|\right|}^2\\ =\arg \min {\left|\right|{P_1}^T(x-Z{f}_{T^2})\left|\right|}^2={\left({\hat{Z}}^T\hat{Z}\right)}^{-1}{\hat{Z}}^{T}x\end{array}---\left(8\right)$

其中，

表示T²统计量重构中被测数据中的正常数据的部分，

P₁表示历史正常数据负载矩阵，Λ表示历史正常数据负载矩阵对应的特征值构成的对角矩阵；

将公式（5）中的故障方向Ζ在主元子空间映射的矩阵

代入公式（8）中，获得故障量值的求解公式：

$\begin{array}{c}{f}_{T^2}{=\left\{{\left[P_1{P}_1^T{P}^{*}\right]}^T\right[P_1{P}_1^T{P}^{*}\left]\right\}}^{-1}[P_1{P}_1^T{P}^{*}{]}^{T}x\\ =[P^{*T}{P}_1{P}_1^T{P}^{*}{]}^{-1}{P}^{*T}{P}_1{P}_1^{T}x\end{array}---\left(9\right).$

6.根据权利要求1所述的基于重构识别分析的青霉素发酵过程故障监测方法，其特征在于，步骤9所述的获得在T²统计量重构中被测数据矩阵的T²统计值，公式如下：

T²=t^TΛ^-1t                （10）

其中，t=x^*P₁，Λ=diag{λ₁,λ₂,…,λ_R}是对角矩阵，λ₁,λ₂,…,λ_R是历史正常数据负载矩阵对应的特征值；

通过F分布获得T²统计量重构中的置信限，公式如下：

$T_{R,N,\mathrm{\α}}^2~\frac{R(N-1)}{N-R}{F}_{R,N-R,\mathrm{\α}}---\left(11\right)$

其中，

表示T²统计量重构中的置信限，α表示置信度，R表示主元的个数，F_R,N-R,α表示置信度为α、第一自由度为R、第二自由度为N-R的F分布。

7.根据权利要求1所述的基于重构识别分析的青霉素发酵过程故障监测方法，其特征在于，步骤8所述的获得SPE统计量重构中x^* _SPE被测数据矩阵的SPE统计值，公式如下：

$\mathrm{SPE}={\left|\right|\tilde{C}{x^{*}}_{\mathrm{SPE}}\left|\right|}^2={\left(\tilde{C}{x^{*}}_{\mathrm{SPE}}\right)}^T\tilde{C}{x^{*}}_{\mathrm{SPE}}={x^{*T}}_{\mathrm{SPE}}\tilde{C}{x^{*}}_{\mathrm{SPE}}---\left(12\right)$

其中，SPE表示SPE统计量重构中被测数据矩阵的SPE统计值；

通过χ²分布获得SPE统计量重构中的置信限，公式如下：

${\mathrm{SPE}}_{\mathrm{\α}}~{\mathrm{g\χ}}_h^2---\left(13\right)$

其中，g=b/2a，h=2a²/b，其中，a表示历史正常数据的平均值；b表示历史正常数据的方差；h表示置信度

表示置信度为h的χ²分布。

Images