基于測地線距離統計量的多工況間歇過程監測

郭小萍，李婷，李元

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

引 言

隨著現代社會對多品種、多規格和高質量產品更迫切的市場需求，工業生產更加倚重于生產小批量、高附加值產品的間歇過程。在間歇生產過程中，原料組分、生產負荷、產品特性等因素的改變，過程數據呈現多工況特征。目前，多工況間歇過程的性能監測與故障診斷技術正日益受到工業界和學術界的關注和重視[1-2]。

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

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

1 基于測地線距離統計量的間歇過程監測

1.1 基于IDijkstra算法的測地線距離

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

圖1 測地線距離Fig.1 Schematic diagram of geodesic distance

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

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

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

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

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

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

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

1.2 多工況批次數據預處理

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

1.3 測地線距離統計量的監測方法

基于測地線距離統計量GDS（geodesic distance statistic）的故障監測方法由兩部分構成：模型建立和故障監測，如圖3所示。

（1）模型建立主要步驟

圖2 多工況批次過程數據預處理Fig.2 Pretreatment of multimode batch process data

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

圖3 基于GDS的過程監測方法Fig.3 Process monitoring based on GDS

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

③ 利用 IDijkstra算法計算各批次之間的最短測地線距離矩陣D。根據D選擇各批次數據點的k近鄰批次。

（2）故障監測主要步驟

① 對于新來的一批待檢測樣本，數據預處理后按批次展開，利用建模數據的均值和方差進行標準化，然后用建模的負載矩陣進行PCA降維，降維后的數據用xnew(1×l)表示。

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

③ 利用IDijkstra算法計算m+1個點之間的測地線最短距離矩陣CD。根據最短距離CD選擇各數據點的k近鄰。

2 仿真實例

2.1 數值仿真例子

人工三臂螺旋的每一個臂作為一個工況的數據，共產生3類工況的建模和檢測數據。螺旋一個臂表示為

正常條件下采集螺旋每個臂的351個樣本作為一個工況的數據，3個工況共產生1053個采樣作為建模數據，進行全局建模。通過調節噪聲幅度a=12來產生故障數據，每個臂產生21個故障數據，3個工況共產生63個采樣作為故障監測數據。數據分布如圖4所示。

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

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

圖5 基于不同方法的過程監測Fig.5 Process monitoring based on different methods

表1 數值例子各方法檢測對比Table 1 Detection contrast of numerical example with different methods

表2 工況1檢測結果對比Table 2 Test results contrastion of mode 1

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

表3 工況2檢測結果對比Table 3 Test results contrastion of mode 2

表4 工況3檢測結果對比Table 4 Test results contrastion of mode 3

2.2 半導體工業過程應用

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

在本例中，PCA降維的維數l=3，鄰接加權圖的近鄰數h=10；Dα統計量的近鄰數k=3，偏離系數α=0.7。將所提GDS算法與KPCA、kNN以及PC-kNN算法對半導體工業過程進行故障監測，對比結果如圖6所示。各種方法檢測結果對比如表5所示。

表5 不同方法對半導體3個工況20種故障的檢測對比Table 5 Detection contrast of 20 semiconductor faults from 3 modes with different methods

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

3 結 論

本文提出了一種基于測地線距離統計量的多工況間歇過程監測 GDS算法。該算法通過改進的IDijkstra算法獲得各數據點之間的最短測地線距離，準確描述了多工況數據的非線性位置關系；通過構造故障監測統計量Dα減小了邊緣訓練數據的偏離程度，提高了控制閾值的準確度。利用人工三臂螺旋數值仿真和半導體工業過程應用，驗證了所提 GDS算法在故障檢測中的突出優越性。但本文偏離系數的選擇對實驗結果有一定影響，選擇經驗值不具有普遍性，因此偏離系數的選取原則仍待繼續研究。

圖6 基于不同方法的過程監測Fig.6 Process monitoring based on different methods

符 號 說 明

a——決定噪聲幅度的參數

D,Dα——分別為算法模型的測地線距離矩陣，GDS統計量

Dα_lim——GDS統計量的控制閾值

Dα_new——待檢測批次數據的統計量

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

h——鄰接加權圖的近鄰數

k——Dα統計量的近鄰數

m,l——分別為批次數，降維后數據維數

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

p——主元個數

u0,v0——鄰接圖中任意兩個頂點號

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

w(e)——鄰接圖中各頂點之間的權值

xi,xj——第i和j個數據點

α——距離的偏離系數

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

ξ,τ,ζ——原始數據空間3個方向的維度坐標

[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(數據驅動的工業過程故障診斷技術)[M].Beijing: Science Press, 2011

[4]Wang Haiqing（王海清）, Song Zhihuan（宋執環）, Wang Hui（王慧）. Fault detection behavior analysis based on PCA process monitoring [J].Journal of Chemical Industry and Engineering(China)(化工學報), 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(系統仿真學報), 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(自動化學報),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(上海交通大學學報),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(哈爾濱理工大學學報), 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(遙感技術與應用),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