郭小萍，李婷，李元

（沈陽化工大學(xué)信息工程學(xué)院，遼寧 沈陽 110142）

引 言

隨著現(xiàn)代社會(huì)對(duì)多品種、多規(guī)格和高質(zhì)量產(chǎn)品更迫切的市場(chǎng)需求，工業(yè)生產(chǎn)更加倚重于生產(chǎn)小批量、高附加值產(chǎn)品的間歇過程。在間歇生產(chǎn)過程中，原料組分、生產(chǎn)負(fù)荷、產(chǎn)品特性等因素的改變，過程數(shù)據(jù)呈現(xiàn)多工況特征。目前，多工況間歇過程的性能監(jiān)測(cè)與故障診斷技術(shù)正日益受到工業(yè)界和學(xué)術(shù)界的關(guān)注和重視[1-2]。

PCA方法[3-4]是目前廣泛應(yīng)用于過程監(jiān)測(cè)的數(shù)據(jù)驅(qū)動(dòng)方法，MacGregor等[5]提出的多向主元分析（multiway principal component analysis,MPCA）是PCA在間歇過程監(jiān)測(cè)中的應(yīng)用，表明了 MPCA能夠很好地處理間歇過程數(shù)據(jù)，達(dá)到過程檢測(cè)和診斷的目的。但MPCA在進(jìn)行故障檢測(cè)時(shí)假設(shè)過程數(shù)據(jù)服從正態(tài)分布且具有線性關(guān)系，不適用于多工況非線性過程，這在一定程度上限制了MPCA的應(yīng)用范圍[5]。核主元分析（kernel principal component analysis，KPCA)方法[6-7]是常用的非線性分析方法，基于核進(jìn)行非線性數(shù)據(jù)映射的分析方法，避免了對(duì)非線性映射具體形式的求取，只需定義特征空間中的內(nèi)積即可。但這些方法中核函數(shù)及其中參數(shù)的選取目前還沒有公認(rèn)有效的方法，限制了這些算法在更大范圍的應(yīng)用。He等[8-9]提出基于歐氏距離k近鄰的kNN算法，能夠比較有效地實(shí)現(xiàn)復(fù)雜批次過程故障檢測(cè)。該算法中局部近鄰的選取以最短歐氏距離為指標(biāo)，數(shù)據(jù)空間中兩點(diǎn)之間的歐氏距離是這兩點(diǎn)連線的直線距離，不能準(zhǔn)確描述非線性數(shù)據(jù)的復(fù)雜位置關(guān)系。測(cè)地線[10-11]是歐氏空間中的直線段在黎曼流形上的推廣，它是黎曼流形上的一條光滑曲線，在局部范圍內(nèi)測(cè)地線總是最短線，能夠很好地反映兩個(gè)樣本點(diǎn)在該流形上的幾何特性，因此在物理學(xué)、工程技術(shù)學(xué)、大地測(cè)量學(xué)和統(tǒng)計(jì)學(xué)等領(lǐng)域中越來越受到重視并得以廣泛研究和應(yīng)用[12-14]。在間歇過程故障監(jiān)測(cè)中，采用測(cè)地線計(jì)算批次樣本之間的最短距離，能夠準(zhǔn)確地體現(xiàn)具有非線性特點(diǎn)的多工況數(shù)據(jù)在空間上的復(fù)雜位置關(guān)系，為局部近鄰的選取提供準(zhǔn)確的距離指標(biāo)。

本文提出一種基于 GDS的多工況間歇過程監(jiān)測(cè)方法，通過構(gòu)造鄰接圖來表示多工況數(shù)據(jù)的局部幾何結(jié)構(gòu)，在一定假設(shè)條件下確定數(shù)據(jù)之間的相似性關(guān)系，構(gòu)建全局模型進(jìn)行故障監(jiān)測(cè)。首先，將按批次展開和標(biāo)準(zhǔn)化的多工況過程數(shù)據(jù)采用 PCA方法進(jìn)行降維；然后，提出采用 IDijkstra算法代替Dijkstra算法[15]，獲得賦權(quán)鄰接矩陣，計(jì)算各批次之間的最短測(cè)地線距離，通過構(gòu)建統(tǒng)計(jì)量Dα進(jìn)行故障監(jiān)測(cè)；最后，在人工三臂螺旋數(shù)值仿真和半導(dǎo)體工業(yè)實(shí)例中驗(yàn)證所提算法的有效性。

1 基于測(cè)地線距離統(tǒng)計(jì)量的間歇過程監(jiān)測(cè)

1.1 基于IDijkstra算法的測(cè)地線距離

測(cè)地線距離是黎曼流形上的一條光滑曲線，對(duì)于具有非線性特點(diǎn)的批次過程數(shù)據(jù)，歐氏距離把兩數(shù)據(jù)點(diǎn)之間的直線距離作為兩點(diǎn)之間的最短距離，然而它們沿著非線性流形的距離遠(yuǎn)遠(yuǎn)大于它們的輸入距離（圖1）[16]。對(duì)于具有強(qiáng)非線性流形結(jié)構(gòu)的swiss-roll數(shù)據(jù)集，從圖1（a）可以看出，高維空間中實(shí)線表示的兩點(diǎn)之間最短測(cè)地線距離遠(yuǎn)大于虛線表示的歐氏距離；從圖 1（c）可以看出，高維樣本在低維空間投影后，紅線表示的是最短測(cè)地線距離與藍(lán)線表示的歐氏距離相比，能夠更準(zhǔn)確地體現(xiàn)數(shù)據(jù)間的非線性位置關(guān)系。GDS的主要思想就是首先使用最近鄰圖上的最短路徑得到近似的測(cè)地線距離，代替不能表示數(shù)據(jù)內(nèi)部非線性流形結(jié)構(gòu)的歐氏距離。測(cè)地線距離的計(jì)算一般采用Dijkstra算法。

圖1 測(cè)地線距離Fig.1 Schematic diagram of geodesic distance

Dijkstra算法[17]在圖論中得到廣泛應(yīng)用，可以用來計(jì)算加權(quán)圖里的兩個(gè)指定頂點(diǎn)u0與v0之間的最短距離，也可計(jì)算u0到所有頂點(diǎn)之間的最短距離。已知圖G(V,E)及每條邊的權(quán)w(e)，對(duì)于任意指定的二點(diǎn)u0,v0∈V(G)，尋找路P(u0,v0)，使得

其中，Ω是從u0到v0所有路的集合，w(P)是路P上的各邊權(quán)之和。這樣的路P(u0,v0)稱為從u0到v0的最短路，w(P(u0,v0))稱為從u0到v0的最短距離。主要方法步驟如下。

① 置l(u0)=0；對(duì)v≠v0，置l(v)=∞；S0={u0}；i=0。

② 對(duì)每個(gè)v?Si，用min{l(v),l(ui)+w(ui,v)}代替l(v)。計(jì)算l(v)}，并把達(dá)到這個(gè)最小值的一個(gè)頂點(diǎn)記為ui+1，置

③ 若i=V? 1，則停止；若i＜V? 1，則用i+1代替i，并返回步驟②進(jìn)入下一次迭代。

本文對(duì)該方法進(jìn)行改進(jìn)提出了 IDijkstra算法，通過計(jì)算各個(gè)批次數(shù)據(jù)之間的鄰接關(guān)系權(quán)值構(gòu)造鄰接加權(quán)圖，鄰接加權(quán)圖上的每個(gè)頂點(diǎn)表示多工況數(shù)據(jù)的每個(gè)批次，得到任意兩批次之間的最短距離矩陣，使其更容易在計(jì)算機(jī)上進(jìn)行矩陣運(yùn)算。具體步驟為：將已知圖G(V,E)存儲(chǔ)在賦權(quán)鄰接矩陣W(wij)n×n里，其中，n為批次數(shù)，xi表示各批次對(duì)應(yīng)的數(shù)據(jù)點(diǎn)

可求出點(diǎn)xi(t=1,2,3,…,n)到其他各點(diǎn)之間的最短距離，其中d(xi,xj)為兩批次xi與xj之間的歐氏距離；h為鄰接加權(quán)圖的近鄰個(gè)數(shù)；將Dijkstra算法運(yùn)用在W矩陣的第t行進(jìn)而得到各批次數(shù)據(jù)點(diǎn)之間的測(cè)地線最短距離矩陣D。

1.2 多工況批次數(shù)據(jù)預(yù)處理

在進(jìn)行 GDS過程監(jiān)測(cè)之前，需要對(duì)多工況間歇過程數(shù)據(jù)進(jìn)行預(yù)處理，如圖2所示。首先，將每個(gè)工況的三維數(shù)據(jù)(I×J×K)按照批次方向展開成二維數(shù)據(jù)XC(I×JK)，其中I表示批次，J表示變量，K表示采樣時(shí)刻。然后，將展開后的各工況數(shù)據(jù)按照X= [X1;X2;… ;XC]的方式進(jìn)行組合，形成全局的建模數(shù)據(jù)X(cNc×JK)。用 Z-score方法對(duì)X進(jìn)行標(biāo)準(zhǔn)化，形成零均值方差為1的整體建模數(shù)據(jù)。

1.3 測(cè)地線距離統(tǒng)計(jì)量的監(jiān)測(cè)方法

基于測(cè)地線距離統(tǒng)計(jì)量GDS（geodesic distance statistic）的故障監(jiān)測(cè)方法由兩部分構(gòu)成：模型建立和故障監(jiān)測(cè)，如圖3所示。

（1）模型建立主要步驟

圖2 多工況批次過程數(shù)據(jù)預(yù)處理Fig.2 Pretreatment of multimode batch process data

① 多工況批次數(shù)據(jù)預(yù)處理，產(chǎn)生全局建模數(shù)據(jù)，并利用 PCA方法進(jìn)行降維，降維后的矩陣為X(m×l)，其中，m為批次數(shù)，l為降維后的維數(shù)。

圖3 基于GDS的過程監(jiān)測(cè)方法Fig.3 Process monitoring based on GDS

② 計(jì)算X各數(shù)據(jù)點(diǎn)之間的歐氏距離d(xi,xj)，利用式（2）計(jì)算各數(shù)據(jù)點(diǎn)之間的權(quán)重，求得賦權(quán)鄰接矩陣W(wij)。

③ 利用 IDijkstra算法計(jì)算各批次之間的最短測(cè)地線距離矩陣D。根據(jù)D選擇各批次數(shù)據(jù)點(diǎn)的k近鄰批次。

（2）故障監(jiān)測(cè)主要步驟

① 對(duì)于新來的一批待檢測(cè)樣本，數(shù)據(jù)預(yù)處理后按批次展開，利用建模數(shù)據(jù)的均值和方差進(jìn)行標(biāo)準(zhǔn)化，然后用建模的負(fù)載矩陣進(jìn)行PCA降維，降維后的數(shù)據(jù)用xnew(1×l)表示。

② 計(jì)算xnew與每一批建模數(shù)據(jù)的歐氏距離dnew(xnew,xi)，組成新的距離矩陣Dnew是m+1維的方陣。利用式（2）計(jì)算m+1個(gè)點(diǎn)之間的權(quán)重，求得賦權(quán)鄰接矩陣CW。

③ 利用IDijkstra算法計(jì)算m+1個(gè)點(diǎn)之間的測(cè)地線最短距離矩陣CD。根據(jù)最短距離CD選擇各數(shù)據(jù)點(diǎn)的k近鄰。

2 仿真實(shí)例

2.1 數(shù)值仿真例子

人工三臂螺旋的每一個(gè)臂作為一個(gè)工況的數(shù)據(jù)，共產(chǎn)生3類工況的建模和檢測(cè)數(shù)據(jù)。螺旋一個(gè)臂表示為

正常條件下采集螺旋每個(gè)臂的351個(gè)樣本作為一個(gè)工況的數(shù)據(jù)，3個(gè)工況共產(chǎn)生1053個(gè)采樣作為建模數(shù)據(jù)，進(jìn)行全局建模。通過調(diào)節(jié)噪聲幅度a=12來產(chǎn)生故障數(shù)據(jù)，每個(gè)臂產(chǎn)生21個(gè)故障數(shù)據(jù)，3個(gè)工況共產(chǎn)生63個(gè)采樣作為故障監(jiān)測(cè)數(shù)據(jù)。數(shù)據(jù)分布如圖4所示。

圖4 人工三臂螺旋的三維空間圖Fig.4 Three-dimensional map of artificial three-arm spiral

在該數(shù)值例子中，p=3,h=10,k=3,α=0.7。分別用KPCA、kNN、PC-kNN[18]以及所提的 GDS算法對(duì)數(shù)值例子進(jìn)行檢測(cè)，對(duì)比效果如圖5所示，各種方法檢測(cè)結(jié)果對(duì)比如表1所示，3個(gè)工況數(shù)據(jù)各方法檢測(cè)結(jié)果對(duì)比如表2～表4所示。

圖5 基于不同方法的過程監(jiān)測(cè)Fig.5 Process monitoring based on different methods

表1 數(shù)值例子各方法檢測(cè)對(duì)比Table 1 Detection contrast of numerical example with different methods

表2 工況1檢測(cè)結(jié)果對(duì)比Table 2 Test results contrastion of mode 1

從表1可以看出，對(duì)于具有多工況特點(diǎn)的三臂螺旋數(shù)值仿真，KPCA算法檢測(cè)效果不好，這與算法中核函數(shù)及其參數(shù)的選擇有很大關(guān)系；基于歐氏距離k近鄰的kNN和PC-kNN算法能夠檢測(cè)出大部分的故障數(shù)據(jù)，但三臂螺旋正常樣本邊緣的故障樣本檢測(cè)效果不好；而本文提出的 GDS算法能夠檢測(cè)出全部的63個(gè)故障樣本，數(shù)據(jù)點(diǎn)之間測(cè)地線距離的計(jì)算更好地體現(xiàn)了數(shù)據(jù)之間的位置關(guān)系，與D2統(tǒng)計(jì)量相結(jié)合增大了建模數(shù)據(jù)邊緣正常樣本與故障樣本的差異，檢測(cè)效果明顯提高。從表2～表4可以看出，在三臂螺旋數(shù)據(jù)的多工況數(shù)值仿真中，GDS算法對(duì)各個(gè)工況的故障敏感度都很高，檢測(cè)效果非常明顯。

表3 工況2檢測(cè)結(jié)果對(duì)比Table 3 Test results contrastion of mode 2

表4 工況3檢測(cè)結(jié)果對(duì)比Table 4 Test results contrastion of mode 3

2.2 半導(dǎo)體工業(yè)過程應(yīng)用

半導(dǎo)體工業(yè)數(shù)據(jù)來源于在Lam 9600上進(jìn)行的半導(dǎo)體鋁蝕反應(yīng)[19-20]，是包括多個(gè)工況的間歇生產(chǎn)過程。本文選取半導(dǎo)體3個(gè)工況的數(shù)據(jù)進(jìn)行多工況故障監(jiān)測(cè)，從40個(gè)測(cè)量變量中選擇17個(gè)變量作為監(jiān)測(cè)變量。采集107個(gè)批次的正常數(shù)據(jù)，95批次用來建模，其中1～30批次取自第1個(gè)工況，31～62批次取自第2個(gè)工況，63～95批次取自第3個(gè)工況；12批次用來驗(yàn)證建模的準(zhǔn)確性，其中1～4批取自第1個(gè)工況，5～8批取自第2個(gè)工況，9～12批取自第 3個(gè)工況。同時(shí)，通過人為改變 TCP功率、RF功率等參數(shù)分別采集來自3個(gè)工況20批次的故障數(shù)據(jù)，其中 1～6批次屬于工況一故障，7～13批次屬于工況二故障，14～20批次屬于工況三故障。通過95個(gè)3個(gè)工況混合批次的數(shù)據(jù)進(jìn)行全局建模，利用3個(gè)工況12批正常數(shù)據(jù)檢驗(yàn)建模準(zhǔn)確性，最終驗(yàn)證各工況故障批次能否及時(shí)準(zhǔn)確地檢測(cè)出來。

在本例中，PCA降維的維數(shù)l=3，鄰接加權(quán)圖的近鄰數(shù)h=10；Dα統(tǒng)計(jì)量的近鄰數(shù)k=3，偏離系數(shù)α=0.7。將所提GDS算法與KPCA、kNN以及PC-kNN算法對(duì)半導(dǎo)體工業(yè)過程進(jìn)行故障監(jiān)測(cè)，對(duì)比結(jié)果如圖6所示。各種方法檢測(cè)結(jié)果對(duì)比如表5所示。

表5 不同方法對(duì)半導(dǎo)體3個(gè)工況20種故障的檢測(cè)對(duì)比Table 5 Detection contrast of 20 semiconductor faults from 3 modes with different methods

從表5可以看出，對(duì)于具有3個(gè)工況的半導(dǎo)體間歇過程數(shù)據(jù)，基于核的KPCA算法故障監(jiān)測(cè)效果不理想，T2統(tǒng)計(jì)量只能檢測(cè)出7個(gè)故障批次，出現(xiàn)這種情況的原因是KPCA算法中核函數(shù)參數(shù)的選擇沒有確定的標(biāo)準(zhǔn)，在多次試驗(yàn)選優(yōu)的情況下檢測(cè)效果仍不理想?；跉W氏距離k近鄰的 kNN算法和PC-kNN算法能檢測(cè)出16個(gè)故障批次，對(duì)于與邊緣建模樣本距離很近的故障數(shù)據(jù)檢測(cè)效果不佳。而本文提出的GDS算法檢測(cè)出全部20個(gè)故障，檢測(cè)率達(dá)到 100%。利用測(cè)地線距離準(zhǔn)確描述了數(shù)據(jù)之間的非線性位置關(guān)系，構(gòu)造Dα統(tǒng)計(jì)量減小了邊緣建模數(shù)據(jù)的偏離程度，增大了故障數(shù)據(jù)與邊緣建模數(shù)據(jù)的差異，提高了控制閾值的準(zhǔn)確度，使得其他 3種方法無法檢測(cè)出來的第3和第9個(gè)微小故障也能成功檢測(cè)出來，故障監(jiān)測(cè)效果明顯提高。

3 結(jié) 論

本文提出了一種基于測(cè)地線距離統(tǒng)計(jì)量的多工況間歇過程監(jiān)測(cè) GDS算法。該算法通過改進(jìn)的IDijkstra算法獲得各數(shù)據(jù)點(diǎn)之間的最短測(cè)地線距離，準(zhǔn)確描述了多工況數(shù)據(jù)的非線性位置關(guān)系；通過構(gòu)造故障監(jiān)測(cè)統(tǒng)計(jì)量Dα減小了邊緣訓(xùn)練數(shù)據(jù)的偏離程度，提高了控制閾值的準(zhǔn)確度。利用人工三臂螺旋數(shù)值仿真和半導(dǎo)體工業(yè)過程應(yīng)用，驗(yàn)證了所提 GDS算法在故障檢測(cè)中的突出優(yōu)越性。但本文偏離系數(shù)的選擇對(duì)實(shí)驗(yàn)結(jié)果有一定影響，選擇經(jīng)驗(yàn)值不具有普遍性，因此偏離系數(shù)的選取原則仍待繼續(xù)研究。

圖6 基于不同方法的過程監(jiān)測(cè)Fig.6 Process monitoring based on different methods

符 號(hào) 說 明

a——決定噪聲幅度的參數(shù)

D,Dα——分別為算法模型的測(cè)地線距離矩陣，GDS統(tǒng)計(jì)量

Dα_lim——GDS統(tǒng)計(jì)量的控制閾值

Dα_new——待檢測(cè)批次數(shù)據(jù)的統(tǒng)計(jì)量

d(xi,xj)——i和j數(shù)據(jù)點(diǎn)之間的歐氏距離

h——鄰接加權(quán)圖的近鄰數(shù)

k——Dα統(tǒng)計(jì)量的近鄰數(shù)

m,l——分別為批次數(shù)，降維后數(shù)據(jù)維數(shù)

P(u0,v0)——鄰接圖中任意兩點(diǎn)之間的路徑

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

u0,v0——鄰接圖中任意兩個(gè)頂點(diǎn)號(hào)

W(wij)n×n——賦權(quán)鄰接矩陣

w(e)——鄰接圖中各頂點(diǎn)之間的權(quán)值

xi,xj——第i和j個(gè)數(shù)據(jù)點(diǎn)

α——距離的偏離系數(shù)

εe——服從N( 0 ,1)的噪聲

ξ,τ,ζ——原始數(shù)據(jù)空間3個(gè)方向的維度坐標(biāo)

[1]Wang Shu. A study of batch process fault diagnosis and prediction methods based on data [D]. Shenyang: Northeastern University, 2010

[2]Chen Yong(陳勇). Fault diagnosis of the production process based on multiway principal components analysis[D]. Hangzhou: Zhejiang University,2003

[3]Zhou Donghua(周東華), Li Gang(李鋼), Li Yuan(李元). Industrial Process Fault Diagnosis of Data-Driven(數(shù)據(jù)驅(qū)動(dòng)的工業(yè)過程故障診斷技術(shù))[M].Beijing: Science Press, 2011

[4]Wang Haiqing（王海清）, Song Zhihuan（宋執(zhí)環(huán)）, Wang Hui（王慧）. Fault detection behavior analysis based on PCA process monitoring [J].Journal of Chemical Industry and Engineering(China)(化工學(xué)報(bào)), 2002(3): 297-301

[5]Nomikos P, MacGregor J F. Monitoring batch processes using multiway principal components analysis [J].AIChE Journal, 1994,40(8): 1361-1375

[6]Scholk opf B, Smola A, Muller K R. Nolinear component analysis as a kernel eigenvalue problem [J].Neural Computation, 1998, 10(5):1299-1319

[7]Guo Xiaoping(郭小萍), Yuan Jie(袁杰), Li Yuan(李元). Nearest neighbors index of kernel features based batch processes monitoring and simulation research [J].Journal of System Simulation(系統(tǒng)仿真學(xué)報(bào)), 2014, 26(7): 1424-1429

[8]He Q P, Wang J. Fault detection using the k-nearest neighbor rule for semiconductor manufacturing processes [J].IEEE Transactions on Semiconductor Manufacturing, 2007, 20(4): 345-354

[9]Guo Xiaoping(郭小萍), Yuan Jie(袁杰), Li Yuan(李元). Bach process monitoring based on k-nearest neighbor of spatial [J].Journal of Automation(自動(dòng)化學(xué)報(bào)),2014,40(1): 135-142

[10]Bai Zhengguo(白正國), Shen Yibing(沈一兵), Shui Naixiang(水乃翔),et al. Riemannian Introduction(黎曼幾何引論)[M]. Beijing:Higher Education Press, 2004

[11]Jurgen Jost. Riemannian Geometry and Geometric Analysis[M].Berlin: Springer-Verlag, 2002

[12]Quan Yong, Yang Jie. Geodesic distance for support vector machines[J].Acta Automatic Sinica,2005, 31(2): 202-208

[13]Yang Qing(楊慶), Chen Guiming(陳桂明), Liu Jingjie(劉鯖潔), He Qingfei(何慶飛). Kernel principal components analysis based on geodesic distance and its application in gear pump fault recognition[J].Journal of Shanghai Jiao Tong University(上海交通大學(xué)學(xué)報(bào)),2011, 45(11):1632-1637

[14]Wu Xiaoting(吳曉婷). Data dimensionality reduction algorithm research based on manifold learning[D]. Dalian: Liaoning Normal University, 2010

[15]Zhang Fuhao(張福浩), Liu Jiping(劉記平), Li Qingyuan(李青元).Shortest path optimization algorithm based on Dijkstra algorithm [J].Remote Sensing Information(遙感信息), 2004, 2: 38-41

[16]Li Chenbo(李臣波), Liu Rentao(劉潤濤). A shortest path algorithm based on Dijkstra [J].Journal of Harbin University of Science and Technology(哈爾濱理工大學(xué)學(xué)報(bào)), 2008, 13(1): 35-37

[17]Su Baoli(蘇寶麗), Li Ning(李寧). The optimization algorithm of Dijkstra and its application for the best path in GIS system [J].Remote Sensing Technology and Application(遙感技術(shù)與應(yīng)用),2013, 28 (5):866-870

[18]He Q P, Wang J. Principal component based k-nearest- neighbor rule for semiconductor process fault detection//Proceedings of the 2008 American Control Conference[C]. Seattle, WA: IEEE, 2008:1606-1611

[19]Yu J. Multiway discrete hidden Markov model-based approach for dynamic batch process monitoring and fault classi fi cation [J].American Institute of Chemical Engineers Journal, 2012, 58(9):2714-2725

[20]Lee S P, Chao A K, Tsung F, Wong D S H, Tseng S T, Jang S S.Monitoring batch processes with multiple on-off steps in semiconductor manufacturing [J].Journal of Quality Technology,2011, 43(2): 142-157