李宇斌，于 濤

(北京化工大學(xué)信息科學(xué)與技術(shù)學(xué)院，北京 100029)

0 引 言

間歇過程經(jīng)常用于生產(chǎn)高附加值的產(chǎn)品，在生物制藥、食品加工、聚合物反應(yīng)及半導(dǎo)體生產(chǎn)等諸多領(lǐng)域具有廣泛應(yīng)用[1-3]。間歇過程復(fù)雜多變，容易受各種因素的影響，一些微小故障的存在就可能導(dǎo)致產(chǎn)品質(zhì)量降低，甚至產(chǎn)生安全隱患。對間歇過程進(jìn)行有效的故障監(jiān)測，對確保生產(chǎn)安全，保證經(jīng)濟效益具有重要意義。多時段特性是間歇過程的本質(zhì)之一，間歇過程不同時段內(nèi)的變量特征不盡相同，需要對各個時段分別建立模型來描述過程特性。因此，正確有效地對多時段間歇過程進(jìn)行時段劃分，對過程監(jiān)測和階段建模具有重要意義[4-7]。

目前，針對間歇過程的時段劃分方法，基于過程變量數(shù)據(jù)多元統(tǒng)計分析的時段劃分方法應(yīng)用最為廣泛[8]。多向主元分析(Multi-way Principal Component Analysis， MPCA)[9]和多向偏最小二乘(Multi-way Partial Least Squares， MPLS)[10]是基于過程變量多元統(tǒng)計分析的時段劃分方法的基礎(chǔ)。Lu等[11]對間歇過程時間片矩陣進(jìn)行主元分析(Principal Component Analysis， PCA)得到反映過程相關(guān)性的負(fù)載矩陣，利用K-means算法對負(fù)載矩陣聚類，劃分時段；Lu等[12]利用K-means聚類對間歇過程時間片矩陣的偏最小二乘(Partial Least Squares， PLS)回歸參數(shù)矩陣進(jìn)行聚類分析，劃分時段，能夠確定過程變量在特定階段對產(chǎn)品質(zhì)量的影響；于濤等[13]提出一種基于滑動時間窗加權(quán)MPCA的時段劃分方法，對間歇過程數(shù)據(jù)通過滑動窗和加權(quán)處理得到局部特征，利用K-means聚類劃分時段，降低了非平穩(wěn)過程對時段劃分的影響；高學(xué)金等[14]對間歇過程時間片矩陣進(jìn)行核熵成分分析，得到主元，完成初次階段劃分，利用模糊C均值聚類(Fuzzy C-means， FCM)對添加時間片矩陣的擴展核熵負(fù)載矩陣的相似度聚類進(jìn)行二次階段劃分，有效克服了時段劃分過程跳變點錯分的現(xiàn)象。以上時段劃分方法是基于過程變量多元統(tǒng)計分析實現(xiàn)的，對過程變量數(shù)據(jù)突變較為敏感。

頻率響應(yīng)函數(shù)(Frequency Response Function， FRF)反映了系統(tǒng)的動態(tài)特性，具有廣泛的應(yīng)用[15]。Zenzen等[16]將FRF和蝙蝠算法相結(jié)合，實現(xiàn)了對桁架結(jié)構(gòu)損傷位置和嚴(yán)重程度的檢測；單衛(wèi)東等[17]利用FRF識別直升機尾傳動軸系非線性模態(tài)參數(shù)；李志農(nóng)等[18]將非線性輸出FRF用于轉(zhuǎn)子不對中-碰摩耦合故障診斷，能夠識別不對中的嚴(yán)重程度；Lin等[19]建立了高階FRF和裂紋物理參數(shù)之間的聯(lián)系，使用二階FRF識別呼吸裂紋的物理參數(shù)。

間歇過程數(shù)據(jù)具有高維度、非線性的特點。核主元分析(Kernel PCA， KPCA)經(jīng)常用于非線性高維數(shù)據(jù)的降維和特征提取[20]。梁京章等[21]將KPCA與融合密度思想的K-means算法相結(jié)合，實現(xiàn)非線性電力負(fù)荷曲線的有效降維和精確聚類；王玲等[22]利用KPCA提取多元時間序列特征，結(jié)合Gath-Geva聚類實現(xiàn)多元時間序列的模糊分段；Ding等[23]利用KPCA對網(wǎng)絡(luò)數(shù)據(jù)降維，縮短了反向傳播神經(jīng)網(wǎng)絡(luò)的訓(xùn)練時間，提高了網(wǎng)絡(luò)入侵檢測的速度。

傳統(tǒng)的時段劃分方法針對間歇過程的過程變量數(shù)據(jù)劃分時段，對輸入輸出數(shù)據(jù)突變較為敏感。本文基于小波變換估計間歇過程的瞬時頻率響應(yīng)函數(shù)(Instantaneous FRF， IFRF)，利用基于間歇過程的瞬時動態(tài)特性進(jìn)行時段劃分，提出基于IFRF的劃分方法。該方法通過小波變換估計系統(tǒng)的IFRF并進(jìn)行KPCA降維，利用FCM進(jìn)行聚類分析，實現(xiàn)基于小波域IFRF的間歇過程時段劃分。

1 基于非線性MPCA的間歇過程時段劃分

1.1 間歇過程樣本數(shù)據(jù)展開

間歇過程作為典型的批處理過程，其過程數(shù)據(jù)可以表示成一個三維張量F(I×J×K)，其中I為批次個數(shù)，J為變量個數(shù)，K為單批次的采樣點個數(shù)。在對間歇過程數(shù)據(jù)進(jìn)行處理時，通常需要將張量數(shù)據(jù)展開成二維矩陣，按批次展開、按變量展開和按時間展開是最為常用的3種展開方式。按批次展開為：Fi(J×K)，每個切片矩陣表示為第i個批次內(nèi)全部J個變量在全部K個采樣點的數(shù)據(jù)，其中i=1,…,I；按變量展開為：Fj(I×K)，每個切片矩陣表示為第j個變量在全部I個批次全部K個采樣點的數(shù)據(jù)，其中j=1,…,J；按時間展開為：Fk(I×J)，每個切片矩陣表示為第k個采樣點變量全部I個批次全部J個的數(shù)據(jù)，其中k=1,…,K。間歇過程的時段劃分具有一定的時間規(guī)律性，因此本文采用按時間展開的方式對間歇過程數(shù)據(jù)進(jìn)行展開。

1.2 基于KPCA的數(shù)據(jù)降維

PCA方法是一種有效的線性數(shù)據(jù)降維方法，它能夠用一組互不相關(guān)的主元代替原始數(shù)據(jù)，這些主元是原始數(shù)據(jù)的線性組合，包含原始數(shù)據(jù)的主要信息。

為了消除數(shù)據(jù)中量綱的影響，需要對數(shù)據(jù)集X進(jìn)行標(biāo)準(zhǔn)化處理。假設(shè)對原始數(shù)據(jù)標(biāo)準(zhǔn)化處理后得到矩陣X=[x1，x2，…，xN]，xj∈RM，(j=1,2,…,N)，N為樣本數(shù)量，M為變量個數(shù)。協(xié)方差矩陣為：

(1)

對協(xié)方差矩陣進(jìn)行特征值分解，λ為特征值，p為特征向量。按累計方差貢獻(xiàn)率來確定主元個數(shù)v，則PCA降維后的數(shù)據(jù)矩陣為：

T=XPv

(2)

PCA降維主要用于線性數(shù)據(jù)的降維，處理非線性數(shù)據(jù)的能力較差。KPCA方法在PCA方法的基礎(chǔ)上引入核函數(shù)，是對PCA方法的非線性擴展。KPCA方法通過非線性映射將原始數(shù)據(jù)投影到高維空間，將線性不可分?jǐn)?shù)據(jù)轉(zhuǎn)化為線性可分?jǐn)?shù)據(jù)，之后再對高維數(shù)據(jù)進(jìn)行線性PCA降維處理。

計算協(xié)方差矩陣為：

(3)

λp=Cp

(4)

(5)

其中，α1，α2，…，αN為常數(shù)。引入N×N維核函數(shù)矩陣：

Kij=K(xi,xj)=Φ(xi)TΦ(xj)

(6)

則有：

Nλα=Kα

(7)

其中，α是核矩陣K的特征向量。對特征向量p進(jìn)行歸一化后，可得到原始數(shù)據(jù)經(jīng)過非線性映射后在第m個特征向量pm上的投影為：

(8)

應(yīng)用KPCA，通常需要對核矩陣K進(jìn)行處理實現(xiàn)高維空間中心均值化，處理方式如下:

(9)

其中，IN為系數(shù)為1/N的N階單位矩陣。

1.3 基于FCM的時段劃分

FCM聚類是一種無監(jiān)督聚類算法，其引入隸屬度函數(shù)，通過隸屬度的大小對數(shù)據(jù)點進(jìn)行分類，定義如下：

(10)

聚類中心為：

(11)

第j個數(shù)據(jù)xj對第i類的隸屬度函數(shù)為：

(12)

基于非線性MPCA的間歇過程時段劃分通過使用FCM對KPCA降維后的數(shù)據(jù)矩陣進(jìn)行聚類分析來實現(xiàn)，具體步驟如下：

2)對矩陣Fv進(jìn)行FCM聚類，得到第k個采樣點數(shù)據(jù)對第l類的隸屬度函數(shù)ulk，根據(jù)隸屬度函數(shù)劃分類別，實現(xiàn)時段劃分。

2 基于IFRF的間歇過程時段劃分

多時段特性是間歇過程的本質(zhì)特性之一。頻率響應(yīng)函數(shù)反映了系統(tǒng)的動態(tài)特性，基于IFRF的間歇過程時段劃分針對系統(tǒng)的瞬時動態(tài)特性進(jìn)行時段劃分。

2.1 基于WT的IFRF估計

頻率響應(yīng)函數(shù)是系統(tǒng)中輸出信號和輸入信號的頻域形式之比，描述了系統(tǒng)的動態(tài)特性。對于瞬態(tài)激勵系統(tǒng)，通常采用輸出信號和輸入信號的傅里葉變換(FT)之比來估計FRF。

利用FT估計系統(tǒng)的FRF，具體形式如下：

(13)

FT將時域信號轉(zhuǎn)換為頻域信號，只保留了頻域信息，而丟失了時域信息。因此基于FT估計FRF無法反映系統(tǒng)的FRF隨時間的變化。用小波變換(WT)代替傅里葉變換能夠解決這個問題。

對于單輸入單輸出(Single Input Single Output， SISO)系統(tǒng)，使用小波變換估計其FRF為：

(14)

WT(k,ω)[y(k)]和WT(k,ω)[x(k)]分別為輸出信號和輸入信號的小波變換時頻表示。

在輸入輸出信號中經(jīng)常存在高斯白噪聲，設(shè)w(k)～N(0,σ2)，此時有：

(15)

(16)

進(jìn)行WT，可得到：

(17)

(18)

式中，WT(k,ω)[w(k)]為高斯白噪聲的小波變換時頻表示，且:

E|WT(k,ω)[w(k)]|2=σ2‖ψ‖2

(19)

其中，ψ表示為小波基函數(shù)。

得到基于WT的FRF如下：

(20)

在數(shù)值實現(xiàn)時，經(jīng)常采用互功率譜法估計FRF。輸入信號的自功率譜和輸入輸出信號的互功率譜如下：

(21)

(22)

其中，*表示復(fù)共軛。

基于WT自相關(guān)和互相關(guān)功率譜估計FRF為：

(23)

則H(k,ω)中，第k時刻對應(yīng)的所有數(shù)據(jù)即為當(dāng)前時刻的IFRF。

對于多輸入多輸出(Multiple Input Multiple Output， MIMO)系統(tǒng)，假定有R個輸入、S個輸出，每個輸入對每個輸出都存在對應(yīng)的FRF，定義第r個輸入對第s個輸出的FRF為Hrs(k,ω)，其中1≤r≤R，1≤s≤S。Hrs(k,ω)中，第k時刻對應(yīng)的所有數(shù)據(jù)即為當(dāng)前時刻第r個輸入對第s個輸出的IFRF。

MIMO系統(tǒng)基于WT的FRF表示為：

(24)

進(jìn)一步處理有：

(25)

2.2 基于IFRF的間歇過程時段劃分

基于IFRF的間歇過程時段劃分采用IFRF替代間歇過程數(shù)據(jù)進(jìn)行時段劃分，并結(jié)合KPCA算法降維和FCM算法聚類劃分間歇過程時段。

基于IFRF的間歇過程時段劃分步驟如下：

(26)

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

本文設(shè)計具有2個時段的數(shù)值仿真過程驗證基于IFRF的時段劃分方法；并將該方法應(yīng)用于青霉素發(fā)酵過程仿真實驗數(shù)據(jù)的時段劃分，驗證基于IFRF的時段劃分方法的魯棒性。

3.1 數(shù)值仿真實驗

本文設(shè)計具有2個時段的數(shù)值仿真過程，如表1所示。

表1中，v1為輸入變量，v2、v3、v4、v5為輸出變量，構(gòu)成一個單輸入四輸出系統(tǒng)。u為邊界值為[2,3.5]的均勻分布，ζ～(0,0.2)為高斯白噪聲，i為批次，k為采樣時刻。

采集多批次具有120采樣點的數(shù)據(jù)，前50采樣點為時段1，后70采樣點為時段2，采用本文方法進(jìn)行時段劃分，結(jié)果如圖1所示。

圖1 時段劃分結(jié)果

從圖1可以看出，此過程劃分為2個時段，分段點為50，與真值一致，驗證了基于IFRF的時段劃分方法的有效性。

3.2 青霉素發(fā)酵過程實驗

青霉素發(fā)酵過程是經(jīng)典的間歇過程，本文利用Pensim仿真平臺生成青霉素發(fā)酵過程仿真數(shù)據(jù)，設(shè)置每批次時長為400 h，采樣時間間隔為1 h，生成正常批次數(shù)據(jù)和底物流加速率在200 h～300 h時存在不同幅度階躍突變的批次數(shù)據(jù)。選擇底物流加速率作為輸入變量，底物濃度、菌體濃度、青霉素濃度、發(fā)酵液體積、二氧化碳濃度作為輸出變量，組成單輸入多輸出系統(tǒng)。

采用基于非線性MPCA的時段劃分方法對以上數(shù)據(jù)進(jìn)行時段劃分，正常批次過程數(shù)據(jù)展開向量強度圖如圖2(a)所示，降維后的特征向量強度圖如圖2(b)所示。

圖2 過程數(shù)據(jù)向量強度圖

從圖2的過程數(shù)據(jù)強度圖可以初步看出，間歇過程在100 h前存在一個較為明顯的分段點。結(jié)合間歇過程先驗知識，將間歇過程劃分為3個時段，各個突變幅度數(shù)據(jù)的時段劃分結(jié)果如圖3所示，時段劃分的分段點如表1所示。

圖3 基于非線性MPCA的青霉素發(fā)酵過程時段劃分結(jié)果

從圖3和表2可以看出，在底物流加速率存在不同幅度的階躍突變時，基于非線性MPCA的時段劃分方法的劃分結(jié)果中，分段點1和分段點2都存在波動，且分段點2變化較為明顯。

表2 基于非線性MPCA的時段劃分分段點

使用基于IFRF的時段劃分方法對以上數(shù)據(jù)進(jìn)行時段劃分，正常批次數(shù)據(jù)基于小波變換估計的IFRF展開向量強度圖如圖4(a)所示，KPCA降維后的特征向量強度圖如圖4(b)所示。

圖4 IFRF向量強度圖

從圖4的IFRF向量強度圖同樣可以看出，間歇過程在100 h前存在一個明顯的分段點，與基于非線性MPCA的時段劃分方法具有相似的結(jié)果。將間歇過程劃分為3個時段，各個突變幅度數(shù)據(jù)的時段劃分結(jié)果如圖5所示，時段劃分的分段點如表3所示。

圖5 基于IFRF的青霉素發(fā)酵過程時段劃分結(jié)果

表3 基于IFRF的時段劃分分段點

對表3中分段點進(jìn)行統(tǒng)計分析，計算分段點1均值為μ1=44，標(biāo)準(zhǔn)差為σ1=0；分段點2均值為μ2=160.4，標(biāo)準(zhǔn)差為σ2=5.86。結(jié)合圖5可以得出，在青霉素發(fā)酵過程的200 h～300 h底物流加速率存在不同幅度的階躍突變時，采用基于IFRF的時段劃分方法，分段點1基本不變，符合間歇過程實際情況，分段點2存在波動，集中在正常數(shù)據(jù)時段劃分的分段點2周圍。

計算基于非線性MPCA的時段劃分方法的不同幅度突變數(shù)據(jù)分段點1和分段點2的標(biāo)準(zhǔn)差并與本文所提方法進(jìn)行對比，如表4所示。

表4 分段點均值與標(biāo)準(zhǔn)差

從表2～表4可以看出，基于非線性MPCA的時段劃分方法的2個分段點受輸入數(shù)據(jù)突變的影響較大；本文所提基于IFRF的時段劃分方法在輸入數(shù)據(jù)中存在階躍突變時，分段點1不發(fā)生變化，分段點2變化較小，具有比基于非線性MPCA的時段劃分方法更小的標(biāo)準(zhǔn)差。本文所提基于IFRF的時段劃分方法能夠?qū)崿F(xiàn)間歇過程的時段劃分，且更加穩(wěn)定，魯棒性更強。

4 結(jié)束語

本文提出了一種基于IFRF的間歇過程時段劃分方法。利用系統(tǒng)的動態(tài)特性不會跟隨輸入數(shù)據(jù)的變化而改變的特點，用描述系統(tǒng)瞬時動態(tài)特性的瞬時頻率響應(yīng)函數(shù)代替非線性MPCA方法中的過程變量數(shù)據(jù)進(jìn)行時段劃分。利用數(shù)值仿真過程數(shù)據(jù)和青霉素發(fā)酵過程仿真實驗數(shù)據(jù)表明本文所提方法的有效性并將本文所提方法與基于非線性MPCA的時段劃分進(jìn)行對比實驗。結(jié)果表明，所提基于IFRF的時段劃分方法能夠減少輸入數(shù)據(jù)突變對時段劃分結(jié)果的影響，具有較高的魯棒性。