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

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

引 言

間歇生產(chǎn)方式以其多品種、小批量、高附加值的優(yōu)勢成為保證產(chǎn)品高質(zhì)量的重要手段[1-4]。在間歇過程中，操作條件的異常容易引起過程故障，如果影響到產(chǎn)品質(zhì)量，則認為是發(fā)生質(zhì)量相關(guān)的故障，可能造成重大經(jīng)濟損失。例如，發(fā)酵過程作為間歇過程的一種典型過程，由于生產(chǎn)進程是不可人為逆轉(zhuǎn)的，并且發(fā)現(xiàn)故障往往有一定的延遲時間，一旦沒有及時檢測并排除故障，極易導(dǎo)致不必要的成本浪費。所以，建立一個合理、有效的監(jiān)控體系，對生產(chǎn)過程進行實時、精確的故障監(jiān)測，最大限度地抑制質(zhì)量下降和重大經(jīng)濟損失，就顯得至關(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ā)酵過程中的一個典型特性[13-14]，傳統(tǒng)的單一模型應(yīng)用于連續(xù)過程時能取得良好的效果，但是并不適用于間歇過程。因此，許多國內(nèi)外專家針對階段劃分問題進行了探索。Yu 等[15]提出基于高斯混合模型(Gaussian mixture model, GMM)的方法，達到分階段建模的目標。Lu 等[16]提出基于K-means 聚類的階段劃分方法，對時間片歸類并建立監(jiān)控模型。以上方法均將過渡階段硬性歸屬到每個穩(wěn)定階段中，不能很好地反映過渡階段的特性，相鄰兩個階段切換期間的過渡，是一種漸變的動態(tài)趨勢。針對階段劃分中的硬分類問題，Zhao等[17]提出一種階段軟劃分方法，在K 均值方法的基礎(chǔ)上利用負載矩陣區(qū)分階段和階段間過渡的變化，揭示了過程潛在特性的行為規(guī)律。Ng 等[18]提出基于重疊PCA 的分段監(jiān)測方法，通過模糊C均值得到的隸屬度信息劃分穩(wěn)定階段與過渡階段。上述方法解決了硬劃分造成的監(jiān)控性能低下的問題，但均建立在過程數(shù)據(jù)服從靜態(tài)獨立分布的假設(shè)下，忽略了實際過程變量的測量點前后相互關(guān)聯(lián)的特性，因此對過程動態(tài)特征變化的捕捉不夠靈敏。針對劃分階段中的動態(tài)性問題，李征等[19]提出基于信息增量的階段劃分方法，利用信息增量結(jié)合滑動窗技術(shù)捕獲系統(tǒng)的動態(tài)特性。Zhao 等[20-21]提出一系列步進有序的階段劃分方法，以過程特征對故障監(jiān)測性能的影響來捕捉過程的動態(tài)特征變化。但是以上方法僅對時序相關(guān)性進行解析，缺乏對過程數(shù)據(jù)變化的充分解析和判斷，無法區(qū)分過程強動態(tài)變化與正常階段切換，往往得到粗糙甚至錯誤的劃分結(jié)果，導(dǎo)致監(jiān)測結(jié)果跳變點控制限異常的情況。以上大多數(shù)階段劃分方法僅依據(jù)過程變量信息劃分階段，忽略了質(zhì)量變量信息，對于階段內(nèi)變量信息描述不完全，進一步導(dǎo)致劃分結(jié)果不理想，在質(zhì)量相關(guān)故障監(jiān)測的過程中可能給出錯誤的結(jié)論。在同一生產(chǎn)階段內(nèi)的不同時刻，過程變量與質(zhì)量變量間的相關(guān)關(guān)系具有較高的相似性，不同時段間則是有差異的[22]。因此，質(zhì)量變量是無法忽略的關(guān)鍵因素，也是客觀反映間歇過程多階段特性的重要指標。

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

1 基本理論

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

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

在建立CCA 模型時，為了研究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ù)最大化，其目標最大化準則函數(shù)為：

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

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

1.2 慢特征分析

慢特征分析是一種提取動態(tài)信息的方法[25]。SFA通過特征映射可以將輸入的時間序列轉(zhuǎn)變?yōu)榫徛兓妮敵鲂盘?，即輸入信號所隱含的慢特征(slow feature, SF)。給定一個m維輸入信號x(t)=[x1(t),x2(t),…,xm(t)]T，SFA 的目標是找到一個特征映射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)化問題表達式如下：

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

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

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

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

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

SFA 可以將快速變化的成分摒棄，提取出時序樣本中變化最緩慢的本質(zhì)特征，這些特征包含著批次間緩慢波動的動態(tài)演化信息，這種動態(tài)信息不僅能夠表征導(dǎo)致過程動態(tài)時變的潛在驅(qū)動力，而且能夠真實地反映過程變量隨時間變化的總體趨勢，因此SFA 能更加充分地提取發(fā)酵過程的動態(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ì)量變量個數(shù)，質(zhì)量變量包括菌體濃度和產(chǎn)物濃度等；K為每個批次的采樣時刻數(shù)。本文假設(shè)建模所需的各批次操作時間都是等長的。

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

使用SFA 進行動態(tài)特征提取，最小化時間片矩陣在每個采樣時刻上不同批次間的變化，然后將所有批次得到的慢特征進行組合得到慢特征時間片，使用K-means 聚類算法對所有時間片進行再次劃分。為了精確地區(qū)分過程強動態(tài)變化時刻與階段切換時刻，利用DPRF(differential phase recognition factor)指標觀測動態(tài)波動情況，有效綜合利用靜動態(tài)劃分結(jié)果，對階段切換點精準定位，依據(jù)切換點對應(yīng)的采樣時刻k，實現(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個二維時間片矩陣Xk(I×Jx)，k=1,2,…,K。同理，將質(zhì)量測量數(shù)據(jù)Y(I×Jy×K)進行相同的處理得到Y(jié)k(I×Jy)。

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

式中，j=1,2,…,Jx；k=1,2,…,K；xˉk,j為第k采樣時刻第j個過程變量的均值；xk,j為第k采樣時刻第j個過程變量；Sk,j為第k采樣時刻第j個過程變量的標準差。

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

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

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

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

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

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

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

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

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

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

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

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

2.2 離線建模

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

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

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

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

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

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

2.3 在線監(jiān)測

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

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

若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ā)生故障且會導(dǎo)致質(zhì)量故障。

3 實驗結(jié)果與分析

3.1 仿真實驗

青霉素發(fā)酵過程是典型的間歇過程，整個發(fā)酵過程，一般最為關(guān)注的質(zhì)量變量為青霉素濃度，過程監(jiān)測的焦點集中于影響最終產(chǎn)品質(zhì)量的故障。本文基于青霉素仿真平臺Pensim[31-32]進行實驗仿真，表1為初始值和參數(shù)值的取值范圍，具體設(shè)定步驟為：設(shè)定每批次發(fā)酵時間為400 h，采樣時間間隔為1 h，按照表1 對初始條件和參數(shù)值隨機取值，10 個過程變量和2 個質(zhì)量變量用于監(jiān)測，如表2 所示。為了模擬青霉素實際發(fā)酵過程，對變量均進行白噪聲處理。本文選取40批正常數(shù)據(jù)組成樣本集，其中包括過程測量數(shù)據(jù)X(40×10×400)和質(zhì)量數(shù)據(jù)Y(40×2×400)，另外2 個故障批次用于測試，故障設(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ù)時，根據(jù)95%的累計方差貢獻率[33]選取主元數(shù)目為2。按式(16)計算相似度指標Sk(k=1,2,…,400)，作為K-means 算法的輸入樣本。觀察劃分結(jié)果可明顯看出整個生產(chǎn)過程被明顯地初步劃分為4 個階段：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)確定最佳的慢特征個數(shù)，最終選取前5 個慢特征。然后按式(19)計算相似度指標S′k(k=1,2,…,400)，輸入Kmeans 算法進行再次劃分。觀察結(jié)果可知，過程被再次劃分為4 個階段：1～23 h，24～174 h，175～238 h，239～400 h。上述兩步劃分結(jié)果不一定嚴格與菌體生長階段一致，但能夠獲得分階段建立故障監(jiān)測模型時注重的數(shù)據(jù)特征。為了區(qū)分強動態(tài)變化點與階段切換點，利用DPRF 指標觀測動態(tài)波動情況。為清晰地呈現(xiàn)DPRF 指標于兩步階段劃分的當前結(jié)果圖中，此處將所有DPRF 值垂直向上平移了2 個單位。DPRF 絕對值越大表示波動水平越高，最大值對應(yīng)的采樣時刻是最可能的切換點。為了合理判定波動水平高低，通過多次實驗將采樣點DPRF 指標絕對值的平均值作為臨界值，在該仿真實驗中將臨界值確定為1。如果DPRF 的絕對值比1 大則認為波動水平很高，并且隸屬于兩次劃分結(jié)果的非重疊區(qū)域，則可以定位階段切換點，圖5 中使用虛線橢圓標出了6 個階段間的切換點，分別是32、58、180、191、262、275 h 時刻，最后綜合分析劃分1～31 h，59～179 h，192～261 h，276～400 h時刻為穩(wěn)定階段，32～58 h，180～191 h，262～275 h時刻為過渡階段。

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

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

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

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

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

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

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

3.2 算法實際驗證

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

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

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

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

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

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

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

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

圖9 是三種方法對故障批次監(jiān)測結(jié)果。圖9(a)在故障起始階段的漏報較多，圖9(b)故障發(fā)生階段存在較多漏報，對比之下，圖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)計量的故障檢測率顯著提高。說明本文所提階段劃分方法和監(jiān)測模型的有效性。表7為三種方法在故障批次下誤報率的對比結(jié)果。

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

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

4 結(jié) 論

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

符 號 說 明

F——F分布

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

H——保留的慢特征集合

s——慢特征

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

x——輸入信號

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

λ——特征值

ρ——Person系數(shù)

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

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

·t——時間平均

下角標

c——階段

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

new——當前運行批次

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

α——顯著性水平

κ——主元個數(shù)