基于瞬時(shí)頻率響應(yīng)函數(shù)的間歇過(guò)程時(shí)段劃分

李宇斌，于 濤

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

0 引 言

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

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

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

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

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

1 基于非線性MPCA的間歇過(guò)程時(shí)段劃分

1.1 間歇過(guò)程樣本數(shù)據(jù)展開(kāi)

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

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

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

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

(1)

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

T=XPv

(2)

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

計(jì)算協(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的特征向量。對(duì)特征向量p進(jìn)行歸一化后，可得到原始數(shù)據(jù)經(jīng)過(guò)非線性映射后在第m個(gè)特征向量pm上的投影為：

(8)

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

(9)

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

1.3 基于FCM的時(shí)段劃分

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

(10)

聚類中心為：

(11)

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

(12)

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

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

2 基于IFRF的間歇過(guò)程時(shí)段劃分

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

2.1 基于WT的IFRF估計(jì)

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

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

(13)

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

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

(14)

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

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

(15)

(16)

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

(17)

(18)

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

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

(19)

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

得到基于WT的FRF如下：

(20)

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

(21)

(22)

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

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

(23)

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

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

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

(24)

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

(25)

2.2 基于IFRF的間歇過(guò)程時(shí)段劃分

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

基于IFRF的間歇過(guò)程時(shí)段劃分步驟如下：

(26)

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

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

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

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

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

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

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

從圖1可以看出，此過(guò)程劃分為2個(gè)時(shí)段，分段點(diǎn)為50，與真值一致，驗(yàn)證了基于IFRF的時(shí)段劃分方法的有效性。

3.2 青霉素發(fā)酵過(guò)程實(shí)驗(yàn)

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

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

圖2 過(guò)程數(shù)據(jù)向量強(qiáng)度圖

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

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

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

表2 基于非線性MPCA的時(shí)段劃分分段點(diǎn)

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

圖4 IFRF向量強(qiáng)度圖

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

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

表3 基于IFRF的時(shí)段劃分分段點(diǎn)

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

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

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

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

4 結(jié)束語(yǔ)

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