基于聯(lián)合典型變量矩陣的多階段發(fā)酵過程質(zhì)量相關(guān)故障監(jiān)測

高學(xué)金，何紫鶴，高慧慧，齊詠生

（1 北京工業(yè)大學(xué)信息學(xué)部，北京 100124； 2 計(jì)算智能與智能系統(tǒng)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100124；3 內(nèi)蒙古工業(yè)大學(xué)電力學(xué)院，內(nèi)蒙 古呼和浩特 010051）

引 言

間歇生產(chǎn)方式以其多品種、小批量、高附加值的優(yōu)勢成為保證產(chǎn)品高質(zhì)量的重要手段[1-4]。在間歇過程中，操作條件的異常容易引起過程故障，如果影響到產(chǎn)品質(zhì)量，則認(rèn)為是發(fā)生質(zhì)量相關(guān)的故障，可能造成重大經(jīng)濟(jì)損失。例如，發(fā)酵過程作為間歇過程的一種典型過程，由于生產(chǎn)進(jìn)程是不可人為逆轉(zhuǎn)的，并且發(fā)現(xiàn)故障往往有一定的延遲時(shí)間，一旦沒有及時(shí)檢測并排除故障，極易導(dǎo)致不必要的成本浪費(fèi)。所以，建立一個(gè)合理、有效的監(jiān)控體系，對(duì)生產(chǎn)過程進(jìn)行實(shí)時(shí)、精確的故障監(jiān)測，最大限度地抑制質(zhì)量下降和重大經(jīng)濟(jì)損失，就顯得至關(guān)重要[5-6]。合理的質(zhì)量相關(guān)的故障監(jiān)測方法能夠保證間歇過程的安全性、產(chǎn)品質(zhì)量的穩(wěn)定性和生產(chǎn)的高效性，主流方法包括典型相關(guān)分析(CCA)[7-8]以及偏最小二乘(partial least squares,PLS)[9-12]等。

多階段特性是發(fā)酵過程中的一個(gè)典型特性[13-14]，傳統(tǒng)的單一模型應(yīng)用于連續(xù)過程時(shí)能取得良好的效果，但是并不適用于間歇過程。因此，許多國內(nèi)外專家針對(duì)階段劃分問題進(jìn)行了探索。Yu 等[15]提出基于高斯混合模型(Gaussian mixture model, GMM)的方法，達(dá)到分階段建模的目標(biāo)。Lu 等[16]提出基于K-means 聚類的階段劃分方法，對(duì)時(shí)間片歸類并建立監(jiān)控模型。以上方法均將過渡階段硬性歸屬到每個(gè)穩(wěn)定階段中，不能很好地反映過渡階段的特性，相鄰兩個(gè)階段切換期間的過渡，是一種漸變的動(dòng)態(tài)趨勢。針對(duì)階段劃分中的硬分類問題，Zhao等[17]提出一種階段軟劃分方法，在K 均值方法的基礎(chǔ)上利用負(fù)載矩陣區(qū)分階段和階段間過渡的變化，揭示了過程潛在特性的行為規(guī)律。Ng 等[18]提出基于重疊PCA 的分段監(jiān)測方法，通過模糊C均值得到的隸屬度信息劃分穩(wěn)定階段與過渡階段。上述方法解決了硬劃分造成的監(jiān)控性能低下的問題，但均建立在過程數(shù)據(jù)服從靜態(tài)獨(dú)立分布的假設(shè)下，忽略了實(shí)際過程變量的測量點(diǎn)前后相互關(guān)聯(lián)的特性，因此對(duì)過程動(dòng)態(tài)特征變化的捕捉不夠靈敏。針對(duì)劃分階段中的動(dòng)態(tài)性問題，李征等[19]提出基于信息增量的階段劃分方法，利用信息增量結(jié)合滑動(dòng)窗技術(shù)捕獲系統(tǒng)的動(dòng)態(tài)特性。Zhao 等[20-21]提出一系列步進(jìn)有序的階段劃分方法，以過程特征對(duì)故障監(jiān)測性能的影響來捕捉過程的動(dòng)態(tài)特征變化。但是以上方法僅對(duì)時(shí)序相關(guān)性進(jìn)行解析，缺乏對(duì)過程數(shù)據(jù)變化的充分解析和判斷，無法區(qū)分過程強(qiáng)動(dòng)態(tài)變化與正常階段切換，往往得到粗糙甚至錯(cuò)誤的劃分結(jié)果，導(dǎo)致監(jiān)測結(jié)果跳變點(diǎn)控制限異常的情況。以上大多數(shù)階段劃分方法僅依據(jù)過程變量信息劃分階段，忽略了質(zhì)量變量信息，對(duì)于階段內(nèi)變量信息描述不完全，進(jìn)一步導(dǎo)致劃分結(jié)果不理想，在質(zhì)量相關(guān)故障監(jiān)測的過程中可能給出錯(cuò)誤的結(jié)論。在同一生產(chǎn)階段內(nèi)的不同時(shí)刻，過程變量與質(zhì)量變量間的相關(guān)關(guān)系具有較高的相似性，不同時(shí)段間則是有差異的[22]。因此，質(zhì)量變量是無法忽略的關(guān)鍵因素，也是客觀反映間歇過程多階段特性的重要指標(biāo)。

針對(duì)以上問題，本文提出一種基于聯(lián)合典型變量矩陣的兩步階段劃分方法，并在子階段建立CCA故障監(jiān)測模型。該方法首先將沿批次展開的二維數(shù)據(jù)矩陣劃分為不同的時(shí)間片，利用對(duì)每個(gè)時(shí)間片矩陣進(jìn)行典型相關(guān)分析得到的原始典型變量矩陣構(gòu)建聯(lián)合典型變量矩陣，之后對(duì)聯(lián)合典型變量矩陣進(jìn)行K 均值聚類，實(shí)現(xiàn)基于靜態(tài)特征的第1 步劃分；然后采用慢特征分析(SFA)算法提取表征過程動(dòng)態(tài)性的慢特征，并用K 均值算法對(duì)其進(jìn)行第2 步劃分，綜合分析兩步劃分結(jié)果確定最終的穩(wěn)定階段和過渡階段，最后建立分階段CCA 監(jiān)測模型。

1 基本理論

1.1 典型相關(guān)分析

典型相關(guān)分析由Hotelling[23]提出，CCA 是研究兩組變量X與Y相關(guān)性的一種統(tǒng)計(jì)方法。給定一組從生產(chǎn)過程操作中收集到的數(shù)據(jù)，其中矩陣X是過程變量，矩陣Y是質(zhì)量變量。

在建立CCA 模型時(shí)，為了研究X與Y這兩組變量的相關(guān)性，需找到投影方向A和B，通過研究投影ATX與BTY的相關(guān)關(guān)系來確定兩組原始變量間的相關(guān)關(guān)系。

CCA 方法的核心思想就是尋找最優(yōu)的正交投影方向A和B，使得原始變量X和Y在此方向的投影的Person 系數(shù)最大化，其目標(biāo)最大化準(zhǔn)則函數(shù)為：

式中，E為與BTY弱相關(guān)的噪聲信號(hào)，但是與ATX不相關(guān)。

CCA 將X和Y看成整體，通過構(gòu)造新的典型變量矩陣U=ATX，V=BTY，分析兩個(gè)變量組之間的聯(lián)系。在同一生產(chǎn)階段內(nèi)的不同時(shí)刻，過程變量與質(zhì)量變量間的相關(guān)關(guān)系具有較高的相似性。典型變量矩陣U和V在分別包含過程變量與質(zhì)量變量信息的同時(shí)可以確保特征之間的關(guān)聯(lián)最大化[24]，從而突出和增強(qiáng)過程變量與質(zhì)量變量之間關(guān)系的變化對(duì)劃分階段的重要影響，因此典型變量矩陣的相似性可以表征生產(chǎn)過程的階段變化，從而提高階段劃分精度。

1.2 慢特征分析

慢特征分析是一種提取動(dòng)態(tài)信息的方法[25]。SFA通過特征映射可以將輸入的時(shí)間序列轉(zhuǎn)變?yōu)榫徛兓妮敵鲂盘?hào)，即輸入信號(hào)所隱含的慢特征(slow feature, SF)。給定一個(gè)m維輸入信號(hào)x(t)=[x1(t),x2(t),…,xm(t)]T，SFA 的目標(biāo)是找到一個(gè)特征映射g(x)=[g1(x),g2(x),… ,gQ(x)]T，使得慢特征s(t)=[s1(t),s2(t),… ,sQ(t)]T，sq(t)=gq(x(t))變化最慢，SFA 的優(yōu)化問題表達(dá)式如下：

發(fā)酵過程數(shù)據(jù)具有非線性特征，為了從非線性的過程數(shù)據(jù)中提取慢特征信息，眾多學(xué)者分別提出了多種基于SFA的非線性算法。其中比較典型的是二項(xiàng)式的非線性擴(kuò)展方法[25]和核SFA[26]。核SFA 雖然在一定程度上解決了非線性問題，但其主要缺點(diǎn)在于算法容易退化導(dǎo)致過擬合，且難以選擇合適的核函數(shù)進(jìn)行有效的非線性變換。因此，為豐富特征空間、兼顧算法可行性。本文采用二項(xiàng)式的非線性擴(kuò)展方法。給定輸入為x(t)={x1,x2,…,xm}，二次擴(kuò)展函數(shù)為v?(x)=[x1,x2,…,xm,x1x1,x1x2,…,xmxm]，然后中心化擴(kuò)展后的數(shù)據(jù)v?(x)使其均值為0 得到v(x)，其維度為m+m(m+1)/2，這樣就將非線性情況轉(zhuǎn)化成了擴(kuò)展空間里的線性情況。

對(duì)式(9)求取廣義特征值時(shí)，為了方便計(jì)算，利用奇異值分解進(jìn)行白化操作以滿足約束式(7)，SFA的二維幾何解釋如圖1 所示。輸入x(t)及其一階導(dǎo)數(shù)x?(t)的協(xié)方差矩陣用實(shí)線橢圓和虛線橢圓標(biāo)示于圖1(a)中，白化矩陣z(t)及其一階導(dǎo)數(shù)z?(t)的各向同性協(xié)方差用實(shí)線圓圈和虛線橢圓標(biāo)示于圖1(b)中，經(jīng)過白化的慢特征s滿足各向同性高斯分布，s?滿足各向異性高斯分布。因此，慢特征體現(xiàn)了發(fā)酵過程的動(dòng)態(tài)變化[27]。

圖1 SFA的二維幾何解釋示意圖Fig.1 The geometrical interpretation of SFA in two dimensions

發(fā)酵過程是一個(gè)緩慢時(shí)變的生化反應(yīng)，其運(yùn)行軌跡在同一間歇操作周期內(nèi)部呈現(xiàn)出局部動(dòng)態(tài)性，由于間歇過程受相同操作重復(fù)執(zhí)行的影響，初始培養(yǎng)基容量逐漸增長使得批次運(yùn)行軌跡在批次間逐漸變化，觀測數(shù)據(jù)產(chǎn)生慢時(shí)變特性，導(dǎo)致各個(gè)批次間亦呈現(xiàn)明顯的動(dòng)態(tài)變化。圖2 為運(yùn)行軌跡的示意圖。

圖2 慢時(shí)變發(fā)酵過程軌跡示意圖Fig.2 Illustration of slow-varying fermentation process trajectory

SFA 可以將快速變化的成分摒棄，提取出時(shí)序樣本中變化最緩慢的本質(zhì)特征，這些特征包含著批次間緩慢波動(dòng)的動(dòng)態(tài)演化信息，這種動(dòng)態(tài)信息不僅能夠表征導(dǎo)致過程動(dòng)態(tài)時(shí)變的潛在驅(qū)動(dòng)力，而且能夠真實(shí)地反映過程變量隨時(shí)間變化的總體趨勢，因此SFA 能更加充分地提取發(fā)酵過程的動(dòng)態(tài)特征。

2 基于聯(lián)合典型變量矩陣的多階段發(fā)酵過程質(zhì)量監(jiān)測

2.1 基于聯(lián)合典型變量矩陣的兩步階段劃分

發(fā)酵過程的訓(xùn)練集包括歷史過程數(shù)據(jù)和歷史質(zhì)量數(shù)據(jù)，分別以三維矩陣X(I×Jx×K)和Y(I×Jy×K)表示。其中，I為批次數(shù)；Jx、Jy分別為過程、質(zhì)量變量個(gè)數(shù)，質(zhì)量變量包括菌體濃度和產(chǎn)物濃度等；K為每個(gè)批次的采樣時(shí)刻數(shù)。本文假設(shè)建模所需的各批次操作時(shí)間都是等長的。

為考慮質(zhì)量變量對(duì)生產(chǎn)過程階段劃分的影響，本文結(jié)合過程變量的典型變量矩陣U和質(zhì)量變量的典型變量矩陣V，構(gòu)建聯(lián)合典型變量矩陣。將每個(gè)時(shí)間片的聯(lián)合典型變量矩陣使用K-means 聚類算法進(jìn)行初步劃分，由于每個(gè)聯(lián)合典型變量矩陣是對(duì)每個(gè)獨(dú)立時(shí)間片采用CCA 算法得到的，同時(shí)傳統(tǒng)的CCA 算法是靜態(tài)的，所以得到的結(jié)果忽略了每個(gè)時(shí)間片的動(dòng)態(tài)特征，因此認(rèn)為使用K-means 聚類算法是對(duì)原始數(shù)據(jù)的靜態(tài)特征進(jìn)行的聚類，得到的聚類結(jié)果只能反映出原始數(shù)據(jù)大致的分布特征。發(fā)酵過程的動(dòng)態(tài)性往往與生產(chǎn)階段密切相關(guān)，因此有必要將能表征過程動(dòng)態(tài)性的動(dòng)態(tài)特征提取出來進(jìn)行再次劃分。

使用SFA 進(jìn)行動(dòng)態(tài)特征提取，最小化時(shí)間片矩陣在每個(gè)采樣時(shí)刻上不同批次間的變化，然后將所有批次得到的慢特征進(jìn)行組合得到慢特征時(shí)間片，使用K-means 聚類算法對(duì)所有時(shí)間片進(jìn)行再次劃分。為了精確地區(qū)分過程強(qiáng)動(dòng)態(tài)變化時(shí)刻與階段切換時(shí)刻，利用DPRF(differential phase recognition factor)指標(biāo)觀測動(dòng)態(tài)波動(dòng)情況，有效綜合利用靜動(dòng)態(tài)劃分結(jié)果，對(duì)階段切換點(diǎn)精準(zhǔn)定位，依據(jù)切換點(diǎn)對(duì)應(yīng)的采樣時(shí)刻k，實(shí)現(xiàn)間歇過程的階段劃分。圖3 為基于聯(lián)合典型變量矩陣的兩步階段劃分方法的流程圖，具體步驟如下。

圖3 基于聯(lián)合典型變量矩陣的兩步階段劃分Fig.3 Two-step phase partition based on joint canonical variable matrix

(1) 將三維過程測量數(shù)據(jù)X(I×Jx×K)沿批次方向展開為K個(gè)二維時(shí)間片矩陣Xk(I×Jx)，k=1,2,…,K。同理，將質(zhì)量測量數(shù)據(jù)Y(I×Jy×K)進(jìn)行相同的處理得到Y(jié)k(I×Jy)。

(2)將沿批次展開后得到的二維時(shí)間片矩陣Xk(I×Jx)進(jìn)一步進(jìn)行中心化和歸一化：

式中，j=1,2,…,Jx；k=1,2,…,K；xˉk,j為第k采樣時(shí)刻第j個(gè)過程變量的均值；xk,j為第k采樣時(shí)刻第j個(gè)過程變量；Sk,j為第k采樣時(shí)刻第j個(gè)過程變量的標(biāo)準(zhǔn)差。

(3)對(duì)每個(gè)標(biāo)準(zhǔn)化時(shí)間片數(shù)據(jù)進(jìn)行典型相關(guān)分析，得到相應(yīng)的Xk的典型變量矩陣Uk和Yk的典型變量矩陣Vk，Uk和Vk分別完整表達(dá)了Xk、Yk的每個(gè)特征之間的聯(lián)系，以及分別對(duì)Xk、Yk整體相關(guān)關(guān)系的影響，形式如式(13)和式(14)所示：

式中，Ak為第k時(shí)刻過程數(shù)據(jù)Xk(I×Jx)的投影矩陣；Bk為第k時(shí)刻質(zhì)量數(shù)據(jù)Yk(I×Jy)的投影矩陣。

為了最大程度地融合原始數(shù)據(jù)的基本信息，同時(shí)確保特征之間關(guān)聯(lián)最大化。針對(duì)每個(gè)采樣時(shí)刻k(k=1,2,…,K)，將質(zhì)量變量的典型變量矩陣Vk擴(kuò)展到過程變量的典型變量矩陣Uk中，由于Uk和Vk有最大的Person 相關(guān)系數(shù)，因此構(gòu)建相應(yīng)聯(lián)合典型變量矩陣Tk，得到式(15)：

(4) 對(duì)聯(lián)合典型變量矩陣使用K-means 聚類算法進(jìn)行劃分，實(shí)現(xiàn)生產(chǎn)過程的初步劃分。K-means算法原理如下。

K-means 算法的輸入是數(shù)據(jù)點(diǎn)之間的相似度，考慮到第1 步劃分是基于數(shù)據(jù)的靜態(tài)特征，因此測量數(shù)據(jù)常常包含噪聲或奇異值，可能會(huì)導(dǎo)致把當(dāng)前階段的數(shù)據(jù)錯(cuò)分到其他階段中的現(xiàn)象，在線監(jiān)測時(shí)產(chǎn)生漏報(bào)和誤報(bào)，為此將采樣時(shí)間k擴(kuò)展到聯(lián)合典型變量矩陣Tk(k=1, 2, … ,K)中得到擴(kuò)展聯(lián)合典型變量矩陣T?k=[Tkk]。本文利用數(shù)據(jù)集的相似度指標(biāo)[28]定義每個(gè)擴(kuò)展聯(lián)合典型變量矩陣T?k與其他時(shí)刻擴(kuò)展聯(lián)合典型變量矩陣T?i(i=1,2,…,K，且k≠i)的相似度：

(5) 計(jì)算標(biāo)準(zhǔn)化時(shí)間片矩陣X?k的慢特征，首先進(jìn)行非線性擴(kuò)展，即對(duì)標(biāo)準(zhǔn)化數(shù)據(jù)中的輸入樣本進(jìn)行二階非線性擴(kuò)展并中心化擴(kuò)展后的數(shù)據(jù)得到v，用v替換式(8)中的x，即實(shí)現(xiàn)了非線性情況向線性情況的轉(zhuǎn)化，之后的步驟與線性求解相同，需要分別計(jì)算矩陣M和N，并根據(jù)式(9)做廣義特征值分解，求得向量矩陣W。采用式(11)求出所有的慢特征s={si,i=1,2,…,Jx+Jx(Jx+1)/2}，將不同時(shí)刻每個(gè)批次的慢特征組合成為慢特征時(shí)間片Sk(k=1,2,…,K)。通過剔除比所有輸入變量變化都快的慢特征信息[29]，從而確定所要提取的慢特征數(shù)目D，計(jì)算方法如下：

(6)間歇過程與另一種重要生產(chǎn)過程——連續(xù)過程相比，具有其特殊性，主要區(qū)別在于始終處于操作條件改變頻繁和非穩(wěn)定工作狀態(tài)下，即同一個(gè)生產(chǎn)批次中不存在穩(wěn)態(tài)工作點(diǎn)，產(chǎn)品和操作改變頻繁。如圖4 所示，討論了動(dòng)態(tài)波動(dòng)的兩種正常情況。一種是不同階段之間的正常切換，另一種是由于原料變化、外部操作環(huán)境等因素引起的某一階段呈現(xiàn)強(qiáng)動(dòng)態(tài)變化。這種動(dòng)態(tài)波動(dòng)是一種不同于階段切換的隨機(jī)行為。因此，為了提供更精確的劃分結(jié)果，需要對(duì)這兩種情況進(jìn)行區(qū)分。階段識(shí)別因子(phase recognition factor, PRF)[30]指標(biāo)度量不同批次運(yùn)行軌跡在相同采樣時(shí)刻上的變化快慢。

圖4 間歇過程動(dòng)態(tài)波動(dòng)的兩種正常情況Fig.4 Two normal cases for the dynamic deviations of batch process

利用慢特征分析得到的慢特征矩陣Sk,D計(jì)算不同批次在第k個(gè)采樣時(shí)刻上的緩慢變化程度：

PRF 值越大，不同的批次軌跡在相同采樣時(shí)刻上變化越快。相較于發(fā)酵過程的穩(wěn)定階段平穩(wěn)的運(yùn)行模式，過渡階段呈現(xiàn)出更強(qiáng)的過程波動(dòng)和更高的不確定性，因此不同的批次在穩(wěn)態(tài)階段的相似性更顯著，也就是穩(wěn)態(tài)階段的PRF 值通常遠(yuǎn)小于過渡階段的PRF 值。為了觀測動(dòng)態(tài)波動(dòng)情況，定義兩個(gè)相鄰時(shí)間切片之間的DPRF(differential PRF)指標(biāo)：

式中，k=2,3,…,K。基于DPRF 指標(biāo)，能夠觀測動(dòng)態(tài)波動(dòng)情況。當(dāng)DPRF 出現(xiàn)明顯變化時(shí)，其對(duì)應(yīng)的采樣時(shí)刻既可能是同一個(gè)階段內(nèi)的強(qiáng)動(dòng)態(tài)變化點(diǎn)，也可能是處于不同階段切換期間的階段切換點(diǎn)。

(7)為了精準(zhǔn)區(qū)分過程強(qiáng)動(dòng)態(tài)變化點(diǎn)與階段切換點(diǎn)，結(jié)合兩步劃分結(jié)果實(shí)現(xiàn)對(duì)跳變點(diǎn)的判斷：①兩步劃分結(jié)果重疊區(qū)域中間不會(huì)存在階段切換點(diǎn)，可能會(huì)出現(xiàn)強(qiáng)動(dòng)態(tài)變化點(diǎn)；②階段切換點(diǎn)均分布于兩步劃分結(jié)果非重疊區(qū)域。因此，忽略隸屬于兩次劃分結(jié)果的重疊區(qū)域的DPRF 跳變點(diǎn)，將隸屬于兩次劃分結(jié)果非重疊區(qū)域的DPRF 跳變點(diǎn)確定為階段切換點(diǎn)。依據(jù)切換點(diǎn)對(duì)應(yīng)的采樣時(shí)刻k，將生產(chǎn)過程分為穩(wěn)定階段和過渡階段。

2.2 離線建模

在完成基于聯(lián)合典型變量矩陣的階段劃分結(jié)束后，需要分別對(duì)每個(gè)階段數(shù)據(jù)建立監(jiān)測模型，計(jì)算出相應(yīng)統(tǒng)計(jì)量的控制限，用于后續(xù)的過程在線監(jiān)測。

(1)首先以批次展開的方式將每個(gè)階段的歷史過程數(shù)據(jù)Xc(I×Jx×K)轉(zhuǎn)換成Xc(I×KcJx)，其中Kc是階段c對(duì)應(yīng)的采樣時(shí)刻；其次對(duì)Xc(I×KcJx)按式(12)做標(biāo)準(zhǔn)化處理，再沿變量方向展開得到Xc(Jx×KcI)。對(duì)質(zhì)量數(shù)據(jù)Yc(I×Jy×K)做相同處理得到質(zhì)量數(shù)據(jù)矩陣Yc(Jy×KcI)。

(2)對(duì)以變量展開后每個(gè)階段的二維過程數(shù)據(jù)和質(zhì)量數(shù)據(jù)，建立CCA監(jiān)測模型，模型表達(dá)式如下：

式中，Ac、Bc為投影方向；Λκ,c為對(duì)角陣；Ec為殘差矩陣。

(3)分別計(jì)算正常操作條件下的質(zhì)量無關(guān)故障的監(jiān)測統(tǒng)計(jì)量Tx2統(tǒng)計(jì)量、質(zhì)量相關(guān)故障的監(jiān)測統(tǒng)計(jì)量Ty

2統(tǒng)計(jì)量及樣本控制限。

2.3 在線監(jiān)測

(1)在新批次的采樣時(shí)刻k，采集過程變量數(shù)據(jù)xnew,k(1×Jx)和質(zhì)量變量數(shù)據(jù)ynew,k(1×Jy)，采用離線建模的均值和方差對(duì)其標(biāo)準(zhǔn)化。

(2)根據(jù)采樣時(shí)刻k對(duì)應(yīng)的階段，選擇相應(yīng)的監(jiān)控模型，計(jì)算此時(shí)刻的統(tǒng)計(jì)量控制限：

若Tx2未超出控制限，Ty2超出控制限，則說明生產(chǎn)過程發(fā)生質(zhì)量有關(guān)的故障，過程數(shù)據(jù)未發(fā)生故障，質(zhì)量受到除了過程數(shù)據(jù)以外的因素影響造成故障；

若Tx2、Ty2均超出控制限，則說明生產(chǎn)過程發(fā)生質(zhì)量相關(guān)的故障，過程數(shù)據(jù)發(fā)生故障且會(huì)導(dǎo)致質(zhì)量故障。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 仿真實(shí)驗(yàn)

青霉素發(fā)酵過程是典型的間歇過程，整個(gè)發(fā)酵過程，一般最為關(guān)注的質(zhì)量變量為青霉素濃度，過程監(jiān)測的焦點(diǎn)集中于影響最終產(chǎn)品質(zhì)量的故障。本文基于青霉素仿真平臺(tái)Pensim[31-32]進(jìn)行實(shí)驗(yàn)仿真，表1為初始值和參數(shù)值的取值范圍，具體設(shè)定步驟為：設(shè)定每批次發(fā)酵時(shí)間為400 h，采樣時(shí)間間隔為1 h，按照表1 對(duì)初始條件和參數(shù)值隨機(jī)取值，10 個(gè)過程變量和2 個(gè)質(zhì)量變量用于監(jiān)測，如表2 所示。為了模擬青霉素實(shí)際發(fā)酵過程，對(duì)變量均進(jìn)行白噪聲處理。本文選取40批正常數(shù)據(jù)組成樣本集，其中包括過程測量數(shù)據(jù)X(40×10×400)和質(zhì)量數(shù)據(jù)Y(40×2×400)，另外2 個(gè)故障批次用于測試，故障設(shè)置如表3所示。

表1 初始條件及控制參數(shù)設(shè)置Table 1 Initial condition and control parameter settings

表2 建模所用變量Table 2 Variables used for modeling

表3 故障設(shè)置信息Table 3 Fault settings

圖5 為采用本文方法得到的劃分結(jié)果，本文用CCA 處理數(shù)據(jù)時(shí)，根據(jù)95%的累計(jì)方差貢獻(xiàn)率[33]選取主元數(shù)目為2。按式(16)計(jì)算相似度指標(biāo)Sk(k=1,2,…,400)，作為K-means 算法的輸入樣本。觀察劃分結(jié)果可明顯看出整個(gè)生產(chǎn)過程被明顯地初步劃分為4 個(gè)階段：1～86 h，87～200 h，201～286 h，287～400 h。

圖5 青霉素發(fā)酵過程階段劃分結(jié)果Fig.5 Phase partition result in penicillin fermentation process

采用SFA 提取慢特征，通過式(18)確定最佳的慢特征個(gè)數(shù)，最終選取前5 個(gè)慢特征。然后按式(19)計(jì)算相似度指標(biāo)S′k(k=1,2,…,400)，輸入Kmeans 算法進(jìn)行再次劃分。觀察結(jié)果可知，過程被再次劃分為4 個(gè)階段：1～23 h，24～174 h，175～238 h，239～400 h。上述兩步劃分結(jié)果不一定嚴(yán)格與菌體生長階段一致，但能夠獲得分階段建立故障監(jiān)測模型時(shí)注重的數(shù)據(jù)特征。為了區(qū)分強(qiáng)動(dòng)態(tài)變化點(diǎn)與階段切換點(diǎn)，利用DPRF 指標(biāo)觀測動(dòng)態(tài)波動(dòng)情況。為清晰地呈現(xiàn)DPRF 指標(biāo)于兩步階段劃分的當(dāng)前結(jié)果圖中，此處將所有DPRF 值垂直向上平移了2 個(gè)單位。DPRF 絕對(duì)值越大表示波動(dòng)水平越高，最大值對(duì)應(yīng)的采樣時(shí)刻是最可能的切換點(diǎn)。為了合理判定波動(dòng)水平高低，通過多次實(shí)驗(yàn)將采樣點(diǎn)DPRF 指標(biāo)絕對(duì)值的平均值作為臨界值，在該仿真實(shí)驗(yàn)中將臨界值確定為1。如果DPRF 的絕對(duì)值比1 大則認(rèn)為波動(dòng)水平很高，并且隸屬于兩次劃分結(jié)果的非重疊區(qū)域，則可以定位階段切換點(diǎn)，圖5 中使用虛線橢圓標(biāo)出了6 個(gè)階段間的切換點(diǎn)，分別是32、58、180、191、262、275 h 時(shí)刻，最后綜合分析劃分1～31 h，59～179 h，192～261 h，276～400 h時(shí)刻為穩(wěn)定階段，32～58 h，180～191 h，262～275 h時(shí)刻為過渡階段。

本文將提出的算法與其他階段劃分方法進(jìn)行比較，包括步進(jìn)有序時(shí)段劃分(step-wise sequential phase partition, SSPP)方法[20]和基于信息增量矩陣-偏最小二乘(information increment matrix-partial least square, IIMPLS)方法[19]。對(duì)于SSPP 算法，當(dāng)容忍因子值α=4 時(shí)，整個(gè)過程被分為三個(gè)階段，即1～38 h，39～179 h，180～400 h。但是劃分結(jié)果并沒有反映出過渡階段的特性，并不符合實(shí)際的物理過程。圖6 為限制因子θ=3.3 和3.0 時(shí)，IIMPLS 的階段劃分結(jié)果。限制因子θ取值為3.0 時(shí)，46、67、149 和231 h 被確定為切換點(diǎn)。然而IIMPLS和SSPP 的階段數(shù)目均隨可調(diào)參數(shù)變化，需要根據(jù)先驗(yàn)知識(shí)選擇合適的參數(shù)取值才能保證合理的劃分結(jié)果。相比之下，本文提出的算法可以獲得更準(zhǔn)確、更可靠的結(jié)果，并且克服了可調(diào)參數(shù)的影響。

圖6 IIMPLS方法對(duì)青霉素發(fā)酵過程的階段劃分結(jié)果Fig.6 Phase partition result using the IIMPLS method in penicillin fermentation process

為了更好地表明本文方法的可行性和有效性，將監(jiān)控效果與不分階段的CCA 方法(CCA)、不考慮質(zhì)量變量只用過程變量的時(shí)間片矩陣求相似度進(jìn)行兩步階段劃分的CCA、用SSPP 進(jìn)行階段劃分的CCA 方法(SSPP-CCA)和用IIMPLS 進(jìn)行階段劃分的CCA 方法(IIMPLS-CCA)進(jìn)行了對(duì)比驗(yàn)證，其中，SSPP 算法的容忍因子α取值為4，IIMPLS 算法的限制因子θ取值為3.0。用故障檢測率和故障誤報(bào)率定量衡量本文方法的有效性。

首先，考慮的階躍故障1 是人為將攪拌速率變量降低10%，故障發(fā)生時(shí)間為200～400 h。這種故障是模擬青霉素進(jìn)入發(fā)酵后期可能發(fā)生的攪拌電機(jī)故障，攪拌速率的降低致使供氧情況改變，但是此階段需氧量減少，溶氧率超過了需氧率，這種故障引起溶氧下降的幅度對(duì)菌體濃度和產(chǎn)物濃度的影響甚小，因此這種故障是質(zhì)量不相關(guān)的。

其次，考慮的斜坡故障2 是人為將底物流加速率變量降低15%，故障發(fā)生時(shí)間為200～400 h。這種故障是模擬青霉素發(fā)酵后期可能發(fā)生的進(jìn)料泵泄漏故障，底物流加速率減小，直接導(dǎo)致青霉素培養(yǎng)基含量急劇下降，青霉素的產(chǎn)量明顯地受這種故障引起的基質(zhì)濃度下降幅度的影響，因此該故障為質(zhì)量相關(guān)的故障。

從圖7 展示的故障批次2 的監(jiān)測結(jié)果中可以發(fā)現(xiàn)，總地來說五種方法的Tx2統(tǒng)計(jì)量有相近的有效檢測率，但是在故障發(fā)生階段Ty2統(tǒng)計(jì)量較容易產(chǎn)生漏報(bào)，主要原因是與階躍式故障相比，斜坡式的故障發(fā)生比較緩慢，不是突變故障。圖7(a)、(c)和(d)質(zhì)量相關(guān)故障檢測率低下。相比較圖7(e)來看，圖7(b)Ty2統(tǒng)計(jì)量對(duì)200 h 加入的15%斜坡故障在286 h 才監(jiān)測到，存在很大的延遲，原因在于圖7(b)在階段劃分時(shí)未考慮質(zhì)量變量信息，因此在質(zhì)量相關(guān)的故障監(jiān)測時(shí)存在較多漏報(bào)。從圖7(e)的監(jiān)測結(jié)果可以發(fā)現(xiàn)，Ty2統(tǒng)計(jì)量故障檢測率明顯高于其他四種方法，本文所提方法更及時(shí)地顯示質(zhì)量相關(guān)故障的發(fā)生，減小了因傳統(tǒng)對(duì)整個(gè)批次建模而引起的偏差，從而保證了對(duì)質(zhì)量相關(guān)故障的及時(shí)檢出。表5 為五種方法在故障批次2 下誤報(bào)與檢測率的對(duì)比結(jié)果。

表5 故障批次2的誤報(bào)率和檢測率對(duì)比結(jié)果Table 5 False alarm rate and fault detection rate of false batch 2

3.2 算法實(shí)際驗(yàn)證

白介素-2 是治療腫瘤癌癥等疾病的藥物，是通過大腸桿菌發(fā)酵獲得的，其生產(chǎn)過程屬于多階段間歇過程。為了更多更好地獲取具有醫(yī)藥價(jià)值的白介素-2，對(duì)大腸桿菌發(fā)酵過程采取精細(xì)化質(zhì)量控制非常重要。本文選取北京某生物制藥廠的大腸桿菌發(fā)酵過程實(shí)際生產(chǎn)數(shù)據(jù)進(jìn)行算法驗(yàn)證，現(xiàn)場發(fā)酵過程的持續(xù)時(shí)間為6.5 h，每隔10 min 進(jìn)行一次采樣。與仿真實(shí)驗(yàn)相比，現(xiàn)場實(shí)驗(yàn)更易遭受外界環(huán)境的干擾，使得變量面臨的不確定性因素增加，表現(xiàn)出更為強(qiáng)烈的動(dòng)態(tài)特性，因此實(shí)際生產(chǎn)數(shù)據(jù)驗(yàn)證更突顯所提方法的實(shí)際意義和效果。實(shí)驗(yàn)中，選取8 個(gè)過程變量(酸堿度、溶解氧、濃度罐壓、溫度、攪拌速率、補(bǔ)碳、補(bǔ)氮、通風(fēng)速率)和1 個(gè)質(zhì)量變量(菌體濃度)。在實(shí)際發(fā)酵過程中，各批次實(shí)際操作的不同導(dǎo)致批次間數(shù)據(jù)長度不相等。對(duì)于這種情況，有兩種簡單的處理方法。第一種為最短長度法，即找到最短的一批數(shù)據(jù)，將其余批次的數(shù)據(jù)進(jìn)行截取使得它們也具有這個(gè)最短的長度[34]；第二種方法為利用最長批次的數(shù)據(jù)進(jìn)行建模，將較短批次所缺的那部分?jǐn)?shù)據(jù)當(dāng)作“缺損數(shù)據(jù)”來處理[35]。本實(shí)驗(yàn)采用最短長度法，即找到最短的一批數(shù)據(jù)，將其余批次的數(shù)據(jù)進(jìn)行截取使得它們也具有這個(gè)最短的長度，將數(shù)據(jù)處理成周期為390 min，即每個(gè)批次具有39 個(gè)采樣點(diǎn)。采集20 批實(shí)際生產(chǎn)正常數(shù)據(jù)組成歷史訓(xùn)練集，即過程測量數(shù)據(jù)X(20×8×39)和質(zhì)量數(shù)據(jù)Y(20×1×39)。在實(shí)際的工業(yè)發(fā)酵現(xiàn)場，大腸桿菌發(fā)酵過程有3 個(gè)階段：種子培養(yǎng)階段、快速生長階段和產(chǎn)物合成階段，其中第13 和21 個(gè)采樣點(diǎn)為階段轉(zhuǎn)換時(shí)刻，對(duì)應(yīng)菌體生長和產(chǎn)物合成的時(shí)刻[36]。另外由于測量條件的限制，只有在第10～39 個(gè)采樣點(diǎn)之間人為引入通風(fēng)速率的階躍故障批次數(shù)據(jù)，減少幅度為5%。

表4 故障批次1的誤報(bào)率和檢測率Table 4 False alarm rate and fault detection rate of false batch 1

本文所提方法的階段劃分結(jié)果如圖8 所示，首先對(duì)大腸桿菌數(shù)據(jù)的聯(lián)合典型變量矩陣采用式(16)計(jì)算相似度指標(biāo)Sk(k=1,2,…,39)，之后利用Kmeans 算法進(jìn)行第一步劃分。第二步劃分采用同青霉素發(fā)酵過程第二步劃分相同的方法即可。然后利用DPRF 指標(biāo)觀測動(dòng)態(tài)波動(dòng)情況。為清晰地呈現(xiàn)DPRF 指標(biāo)于兩步階段劃分的當(dāng)前結(jié)果圖中，此處將所有DPRF 值垂直向上平移了0.03 個(gè)單位。為了合理判定波動(dòng)水平高低，通過多次實(shí)驗(yàn)，將采樣點(diǎn)DPRF 指標(biāo)絕對(duì)值的平均值作為臨界值，在該實(shí)際驗(yàn)證中將臨界值確定為0.0351。忽略隸屬于兩次劃分結(jié)果的重疊區(qū)域的DPRF 值跳變點(diǎn)，圖8中使用虛線橢圓標(biāo)出了4 個(gè)階段間的切換點(diǎn)，分別是第12、14、21、23 個(gè)采樣點(diǎn)，結(jié)合兩步劃分的結(jié)果，得到最終的劃分結(jié)果為1～11，15～20，24～39 采樣點(diǎn)為穩(wěn)定階段，12～14,21～23 采樣點(diǎn)為過渡階段。

圖8 大腸桿菌發(fā)酵過程階段劃分結(jié)果Fig.8 Phase partition result of E.coli fermentation process

表6 列出了采用IIMPLS(θ=2.1)和本文所提算法得到的劃分結(jié)果的比較，其中實(shí)際切換點(diǎn)之前的劃分點(diǎn)表示為誤差。從表6 可以看出，IIMPLS 算法劃分的第三階段即產(chǎn)物合成階段存在顯著延遲，相比之下，本文所提算法給出的劃分結(jié)果更合理，并且顯著減少延遲。

表6 大腸桿菌發(fā)酵過程階段劃分的對(duì)比結(jié)果Table 6 Comparison of partition results of E.coli fermentation process

為使實(shí)驗(yàn)結(jié)果更具說服力，將監(jiān)控效果與不分階段的CCA 方法(CCA)、不考慮質(zhì)量變量只用過程變量的時(shí)間片矩陣求相似度進(jìn)行兩步階段劃分的CCA方法進(jìn)行了對(duì)比驗(yàn)證。

故障批次是在發(fā)酵中期即產(chǎn)物開始大量產(chǎn)生的階段人為引入通風(fēng)速率故障，通風(fēng)速率變量從第10～39 采樣點(diǎn)發(fā)生階躍式降低，降低幅度為5%。通風(fēng)速率的降低致使供氧情況改變，產(chǎn)物合成階段需氧量增加，引起溶氧率下降，該故障會(huì)導(dǎo)致發(fā)酵罐內(nèi)需氧小于溶氧，導(dǎo)致溶氧率的急劇降低，成為合成大腸桿菌的限制因素，最終影響產(chǎn)物的數(shù)量與質(zhì)量。因此該故障為質(zhì)量相關(guān)的故障。

圖9 是三種方法對(duì)故障批次監(jiān)測結(jié)果。圖9(a)在故障起始階段的漏報(bào)較多，圖9(b)故障發(fā)生階段存在較多漏報(bào)，對(duì)比之下，圖9(c)中Tx

圖9 大腸桿菌發(fā)酵過程故障批次的在線監(jiān)測結(jié)果Fig.9 The online monitoring result of false batch of E.coli fermentation process

2統(tǒng)計(jì)量的故障檢測率顯著提高。說明本文所提階段劃分方法和監(jiān)測模型的有效性。表7為三種方法在故障批次下誤報(bào)率的對(duì)比結(jié)果。

表7 故障批次的誤報(bào)率和檢測率對(duì)比結(jié)果Table 7 False alarm rate and false detection rate of false batch

根據(jù)仿真數(shù)據(jù)和實(shí)際數(shù)據(jù)的階段劃分結(jié)果與在線監(jiān)控效果，一方面證實(shí)了本文所提方法可以將整個(gè)發(fā)酵生產(chǎn)過程準(zhǔn)確地劃分為不同的穩(wěn)定和過渡階段，另一方面證明了本文所提方法可以顯著提高質(zhì)量相關(guān)的故障檢測率，針對(duì)多階段發(fā)酵過程的質(zhì)量相關(guān)故障監(jiān)測提供了合理和有效的解決思路。

4 結(jié) 論

大多數(shù)階段劃分方法僅依據(jù)過程變量信息將整個(gè)生產(chǎn)過程劃分為多個(gè)子階段，忽略了質(zhì)量變量信息，且很少考慮過程的動(dòng)態(tài)性。本文提出一種基于聯(lián)合典型變量矩陣的階段劃分方法實(shí)現(xiàn)發(fā)酵過程質(zhì)量相關(guān)故障監(jiān)測，與前人工作相比，本文方法主要有以下優(yōu)點(diǎn)：(1)構(gòu)建的聯(lián)合典型變量矩陣能合理融合原始過程變量和質(zhì)量變量的特征信息，同時(shí)確保二者關(guān)聯(lián)最大化，解決了由于忽略質(zhì)量變量造成的劃分結(jié)果不理想的問題，提高了質(zhì)量相關(guān)故障檢測率;(2)在利用傳統(tǒng)CCA 提取靜態(tài)特征的基礎(chǔ)上，采用SFA 提取表征過程動(dòng)態(tài)性的慢特征，通過動(dòng)靜協(xié)同的思想，分析原始數(shù)據(jù)的靜態(tài)特征和動(dòng)態(tài)特征實(shí)現(xiàn)兩步劃分，克服了因?yàn)椴豢紤]數(shù)據(jù)動(dòng)態(tài)性而導(dǎo)致的劃分靈敏度下降的缺陷;(3)為了精確地區(qū)分過程強(qiáng)動(dòng)態(tài)變化與正常階段切換，引入相鄰階段識(shí)別系數(shù)(DPRF)指標(biāo)觀測動(dòng)態(tài)波動(dòng)情況，有效綜合利用靜動(dòng)態(tài)劃分結(jié)果，精準(zhǔn)定位階段切換點(diǎn)，從而將整個(gè)操作階段劃分為穩(wěn)定階段和過渡階段，使得劃分結(jié)果符合實(shí)際生產(chǎn)過程。青霉素仿真平臺(tái)與大腸桿菌實(shí)際生產(chǎn)數(shù)據(jù)證明該方法具有可行性和有效性，為發(fā)酵過程的質(zhì)量相關(guān)故障監(jiān)測提供了一種可行方案。

符 號(hào) 說 明

F——F分布

g(·)——慢特征分析的特征映射

H——保留的慢特征集合

s——慢特征

s?——慢特征的一階導(dǎo)數(shù)

x——輸入信號(hào)

x?——輸入信號(hào)的一階導(dǎo)數(shù)

λ——特征值

ρ——Person系數(shù)

v(x)——零均值的二次擴(kuò)展函數(shù)

v?(x)——二次擴(kuò)展函數(shù)

·t——時(shí)間平均

下角標(biāo)

c——階段

m——輸入信號(hào)的維數(shù)

new——當(dāng)前運(yùn)行批次

Q——慢特征分析特征映射的維數(shù)

α——顯著性水平

κ——主元個(gè)數(shù)