李吉俊 章智杰 李其操 董自健

摘? 要：文章使用自回歸因子分析模型（ARFA）對數(shù)據(jù)樣本進(jìn)行動態(tài)過程建模，分析了卡爾曼濾波和EM算法在估計ARFA模型中回歸矩陣參數(shù)A和載荷矩陣參數(shù)C的方法。在此基礎(chǔ)上，提出了一種使用變分貝葉斯EM（VBEM）故障檢測方法，對ARFA模型參數(shù)A和C進(jìn)行推斷和動態(tài)過程故障檢測。仿真實(shí)驗結(jié)果表明，在ARFA模型下，VBEM方法對下文所述的階躍信號、斜坡信號等四類故障的檢測效果要優(yōu)于EM方法對該類故障的檢測效果，并且降低了平均迭代次數(shù)。

關(guān)鍵詞：自回歸因子分析模型;EM算法;變分貝葉斯EM;故障檢測

中圖分類號：TP277;TP391.9? ? ? 文獻(xiàn)標(biāo)識碼：A 文章編號：2096-4706（2020）24-0001-06

Parameter Derivation and Fault Detection of ARFA Model Based on VBEM

LI Jijun，ZHANG Zhijie，LI Qicao，DONG Zijian

（School of Electronic Engineering，Jiangsu Ocean University，Lianyungang? 222005，China）

Abstract：In this paper，the autoregressive factor analysis（ARFA）model is used to model the dynamic process of data samples，and the methods of Kalman filter and EM algorithm to estimate the regression matrix parameter A and load matrix parameter C in ARFA model are analyzed. On this basis，a fault detection method based on variational Bayes EM（VBEM）is proposed to infer ARFA model parameters A and C and detect dynamic process faults. The simulation results show that under ARFA model，the detection effect of VBEM method is better than that of EM method for four kinds of faults，such as step signal and slope signal，and the average number of iterations is reduced.

Keywords：autoregressive factor analysis model;EM algorithm;variational Bayes EM;fault detection

0? 引? 言

隨著科學(xué)和技術(shù)的發(fā)展，在信息化、數(shù)據(jù)化的今天，工業(yè)系統(tǒng)已逐漸呈現(xiàn)出大規(guī)模、強(qiáng)復(fù)雜性、高集成化等特點(diǎn)，對生產(chǎn)過程的安全、平穩(wěn)、產(chǎn)品質(zhì)量等提出了更高的要求。一旦復(fù)雜工程系統(tǒng)發(fā)生故障而又未能得到有效處理，不僅可能導(dǎo)致財產(chǎn)損失或人員傷亡，還可能導(dǎo)致環(huán)境污染造成生態(tài)災(zāi)難[1]。因此，常常采用有效的工業(yè)過程故障檢測與診斷技術(shù)，來提高復(fù)雜工業(yè)過程控制系統(tǒng)的可靠性和安全性，及時有效地發(fā)現(xiàn)和檢測出系統(tǒng)中存在的故障，以此來避免由故障造成的人員傷亡和經(jīng)濟(jì)損失[2]。傳統(tǒng)的主元分析方法（Principal Com-ponent Analysis，PCA）是一種廣泛應(yīng)用于過程監(jiān)控的線性方法，通過將原始高維空間的數(shù)據(jù)投影到低維空間中來消除數(shù)據(jù)間的相關(guān)性。然而在實(shí)際工業(yè)過程短采樣中，還存在著數(shù)據(jù)的自相關(guān)性[3]。為解決這一問題，很多學(xué)者提出了多變量動態(tài)過程建模和故障檢測的方法。KU[4]等人在1995年提出動態(tài)PCA（DPCA）算法，該方法通過把當(dāng)前樣本的前L個樣本加入到數(shù)據(jù)矩陣中，構(gòu)成時間滯后的數(shù)據(jù)增廣矩陣，來提取時序相關(guān)的關(guān)系，從而消除數(shù)據(jù)的自相關(guān)性。LI[5]等提出了一種間接DPCA模型用于動態(tài)過程監(jiān)測，該方法可以有效提取數(shù)據(jù)的自相關(guān)性，卻忽略了數(shù)據(jù)的靜態(tài)特性。

KIM[6]等提出了一種基于概率主成分分析（PPCA）的多變量過程監(jiān)控方法。并指出當(dāng)概率生成模型應(yīng)用于過程監(jiān)控時，大多數(shù)的統(tǒng)計量如單變量Shewhart圖、多變量Shewhart圖、T2和SPE統(tǒng)計量都可以被統(tǒng)一到概率模型的框架中。楊沛武[7]等提出了一種基于動態(tài)PPCA（Dynamic PPCA）的故障檢測技術(shù)，通過增廣時間序列變量來擴(kuò)大數(shù)據(jù)矩陣，進(jìn)而將數(shù)據(jù)的自相關(guān)性納入相應(yīng)的增廣矩陣中來進(jìn)行DPPCA建模。然而該方法無法全部獲取數(shù)據(jù)的互相關(guān)特征。LI[8]等提出了一種新的動態(tài)隱變量模型（DLV）。在DLV中，采用AR模型創(chuàng)建了一個新的殘差，并通過PCA分解來提取剩余部分的靜態(tài)變化。而? 和? 統(tǒng)計量則分別用來監(jiān)測動態(tài)和靜態(tài)變化。因此，DLV能夠同時提取數(shù)據(jù)的自相關(guān)性和互相關(guān)性。WEN[9]等提出了一種基于線性高斯?fàn)顟B(tài)空間模型（LGSSM）的動態(tài)過程監(jiān)測方法。LGSSM通過一階馬爾科夫特性來描述過程的動態(tài)關(guān)系，并對過程數(shù)據(jù)做了降維處理。雖然該方法也能夠同時提取數(shù)據(jù)的自相關(guān)性和互相關(guān)性，但對高階動態(tài)過程的特征提取能力較差。周樂[10]提出了一種自回歸因子分析模型（ARFA），它具有普適性的動態(tài)模型，能夠同時提取數(shù)據(jù)間的動態(tài)和靜態(tài)關(guān)系，有效提髙動態(tài)過程建模和故障檢測效果。

本文在ARFA的等價模型和卡爾曼濾波的基礎(chǔ)上，提出了一種新的變分貝葉斯EM的動態(tài)過程故障檢測方法。該方法推導(dǎo)模型回歸矩陣參數(shù)A和載荷矩陣參數(shù)C的期望和故障檢測的總體思路是：由參數(shù)和數(shù)據(jù)的聯(lián)合分布，根據(jù)變分理論求出A和C的邊緣分布，再求出參數(shù)A和C的后驗期望，并把這兩個參數(shù)期望代入到ARFA的標(biāo)準(zhǔn)化模型和卡爾曼濾波方程中進(jìn)行迭代，直到模型參數(shù)A和C收斂為止。把收斂后的模型參數(shù)代入T2和SPE統(tǒng)計量的相關(guān)公式中進(jìn)行計算。當(dāng)樣本的兩個統(tǒng)計量均在控制限以下時，則認(rèn)為沒有故障發(fā)生;均在控制限以上時，則認(rèn)為有故障發(fā)生[3]。

本文安排為，第一節(jié)介紹ARFA模型及其EM算法;第二節(jié)為利用變分貝葉斯EM方法對ARFA模型進(jìn)行推斷;第三節(jié)為仿真實(shí)驗;第四節(jié)為結(jié)論。

1? ARFA模型及其EM算法

1.1? ARFA模型及其等價系統(tǒng)

ARFA模型形式為：

fk=Azk-1+wk

jk=Cfk+vk? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（1）

其中jk∈RM為測量變量，fk∈RD為動態(tài)因子，zk-1= [fk-1T，fk-2T，…，fk-LT]T∈RDL包含了動態(tài)因子的前L個值。A∈RD×DL為回歸矩陣，C∈RM×D為載荷矩陣。系統(tǒng)的動態(tài)噪聲和測量噪聲分別為wk和vk，它們分別服從均值為零，方差為Q和R的高斯分布[10]。綜上所述，ARFA的模型有效的參數(shù)為{A，C，Q，R}。

為了利用卡爾曼濾波來估計系統(tǒng)狀態(tài)，ARFA模型可以等價為下述線性系統(tǒng)[10]：

（2）

系統(tǒng)（2）的定義與ARFA相似，因此，我們可以得到兩者之間的關(guān)系為和? 。相應(yīng)的動態(tài)噪聲，其中，

1.2? 利用EM算法來估計ARFA模型參數(shù)A和C

K個輸出序列J=（j1，j2，…，jk）T和它們相應(yīng)的動態(tài)因子F=（f1，f2，…，fk）T的聯(lián)合對數(shù)似然函數(shù)為：

（3）

在M步中，利用對數(shù)極大似然函數(shù)分別對A和C求一階偏導(dǎo)數(shù)，得到的模型參數(shù)更新值為：

（4）

（5）

其中，矩陣、和? 分別由式（6）、式（7）、式（8）給出：

（6）

（7）

（8）

1.3? 基于ARFA模型的動態(tài)過程故障檢測

我們構(gòu)建T2和SPE兩種統(tǒng)計量來檢測動態(tài)過程是否發(fā)生故障。與傳統(tǒng)的PCA構(gòu)建T2統(tǒng)計量的方法不同，ARFA中動態(tài)因子fk由于包含了數(shù)據(jù)的動態(tài)關(guān)系而不再相互獨(dú)立，因此，fk已經(jīng)不再適用于直接構(gòu)建統(tǒng)計量以檢測隱變量空間的變化[10]。為此，我們利用動態(tài)因子的估計殘差來構(gòu)建T2統(tǒng)計量為：

（9）

（10）

其中， 由式（10）給出：

（11）

此外，根據(jù)模型的預(yù)測誤差，我們可以構(gòu)造SPE統(tǒng)計量為[10]：

（12）

SPE=ejkTejk? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （13）

根據(jù)正常樣本的統(tǒng)計量，我們可以估計故障檢測的統(tǒng)計限。T2統(tǒng)計量的控制限可以由χ2分布估計為[11]：

T2=T2lim=χ2（D）? ? ? ? ? ? ? ? ? ? ? ? （14）

其中，α為顯著性水平，D為動態(tài)因子的維度。

SPE統(tǒng)計量的控制限同樣由近似的χ2分布估計為[12]：

SPE=SPElim～g·χ2（h）? ? ? ? ? ? ? ? ?（15）

g和h滿足以下方程，其中，α為顯著性水平，g為常數(shù)，h為自由度。

gh=mean（SPEnormal）? ? ? ? ? ? ? ? ? ?（16）

2g2h=var（SPEnormal）? ? ? ? ? ? ? ? ? ?（17）

2? 利用變分貝葉斯EM方法對ARFA模型進(jìn)行推斷

參數(shù)和數(shù)據(jù)的聯(lián)合分布為：

P（A，C，fk，jk）=P（jk|A，C，fk）·P（A，C，fk）（18）

變分近似推導(dǎo)的目的就是為了尋找一個可以逼近待求后驗分布、可以處理的函數(shù)。根據(jù)變分理論[13，14]有：

P（A，C，fk）≈q（A，C，fk）=q（A）q（C）q（fk）（19）

A的邊緣分布為[15]：

（20）

其中，a（j）為A矩陣的第j列， 為協(xié)方差矩陣，SA，（j）為數(shù)量矩陣SA的第j列。

C的邊緣分布為[15]：

（21）

其中，c（j）為C矩陣的第j列， 為協(xié)方差矩陣，SC，（j）為數(shù)量矩陣SC的第j列。

假定A和C矩陣中的超參數(shù)α和γ均如服從以下形式的伽馬分布：

（22）

（23）

相關(guān)數(shù)據(jù)統(tǒng)計矩陣的表達(dá)式由式（24）、式（25）、式（26）、式（27）給出：

（24）

（25）

（26）

（27）

定義和的逆矩陣由式（28）和式（29）給出[15]：

（28）

（29）

qθ（θ）關(guān)于參數(shù)向量θ=（A，C）的變分后驗的期望由式（30）和式（31）給出：

（30）

（31）

3? 仿真實(shí)驗

本文在處理器為CPU Intel i5-5200U 2.20 GHz的筆記本中，采用Matlab（2016b版本）編程語言來對一個二階動態(tài)系統(tǒng)的數(shù)值例子進(jìn)行仿真實(shí)驗。數(shù)值例子為：

mk=Azk-1+wk

nk=Cmk+vk

其中，nk∈R4為k時刻的觀測值，mk∈R2是其相應(yīng)的動態(tài)因子，zk-1=[mk-1T，mk-2T]T∈R4包含了動態(tài)因子的前兩個值。系統(tǒng)的動態(tài)噪聲和測量噪聲分別為wk和vk，它們均服從均值為0，方差為0.01的高斯分布?；貧w矩陣和載荷矩陣分別為：

我們從標(biāo)準(zhǔn)正態(tài)分布中隨機(jī)取出100個2維數(shù)據(jù)作為正常樣本。為了進(jìn)行故障檢測實(shí)驗，我們產(chǎn)生了四種測試集，每個測試集當(dāng)中仍然包括100個樣本點(diǎn)。這些測試集分別引入了四種不同的常見故障，它們分別是：故障1：在第51～100個樣本點(diǎn)上引入幅度為8的階躍信號故障;故障2：在第51～100個樣本點(diǎn)上引入幅度為1.6×（50-t）的斜坡信號故障;故障3：在第51～100個樣本點(diǎn)上引入幅度為0.48×（50-t）2的拋物信號故障;故障4：在第51～100個樣本點(diǎn)上引入均值為12，方差為4的高斯噪聲。

為了驗證本文所提方法的有效性，我們在每種方法下對各種故障分別重復(fù)3 000次實(shí)驗。并記錄每次實(shí)驗中T2統(tǒng)計量和SPE統(tǒng)計量的檢測正確率以及所用方法的迭代次數(shù)，最后再把少部分的異常數(shù)據(jù)剔除掉并求取剩余數(shù)據(jù)在兩種方法下的均值和方差，統(tǒng)計結(jié)果如表1所示。

由表1可以看出，EM方法中的T2統(tǒng)計量對故障1檢測正確率的均值為0.277 3，而在VBEM方法中的均值為0.284 6。兩種方法下的SPE統(tǒng)計量對故障1檢測正確率的均值幾乎相等。VBEM中的T2統(tǒng)計量和SPE統(tǒng)計量對故障2、3、4檢測正確率的均值均高于EM方法。另外，從最后一列的平均值中也可以看到VBEM中的T2統(tǒng)計量和SPE統(tǒng)計量對總的故障檢測正確率的均值也均高于EM方法。所用VBEM方法的平均迭代次數(shù)為4.061 8低于EM方法的平均迭代次數(shù)4.666 2。總的來說，VBEM方法對本文所述四類故障的檢測效果要優(yōu)于EM方法。

如圖1～8所示為在某次實(shí)驗中所得到的兩種方法對各故障的檢測結(jié)果。圖中水平虛線為正常數(shù)據(jù)下的統(tǒng)計量的控制限，實(shí)線曲線分別為數(shù)據(jù)的T2統(tǒng)計量和SPE統(tǒng)計量，當(dāng)統(tǒng)計量超過它們的控制限時，則表明檢測到系統(tǒng)有故障發(fā)生。

如表2所示，在EM方法中T2統(tǒng)計量的漏報率較為嚴(yán)重。在VBEM方法中圖2和圖8的T2統(tǒng)計量的漏報率也較為嚴(yán)重。在EM方法中SPE統(tǒng)計量的漏報率較小。在圖3中SPE統(tǒng)計量的漏報率僅為2%。在圖5中SPE統(tǒng)計量的漏報率僅為14%。在VBEM方法中SPE統(tǒng)計量的檢測正確率為100%，漏報率為0%。

4? 結(jié)? 論

本文針對ARFA模型下的動態(tài)過程故障檢測問題，構(gòu)建了相應(yīng)的T2和SPE兩種統(tǒng)計量，并提出了一種新的變分貝葉斯EM故障檢測方法。在對上文所述四類故障的正確檢測率上，本方法要優(yōu)于傳統(tǒng)的EM故障檢測方法。但本文未對高階動態(tài)系統(tǒng)和實(shí)際的工業(yè)數(shù)據(jù)進(jìn)行仿真實(shí)驗，因此這也是今后的一個研究方向。

參考文獻(xiàn)：

[1] 何瀟，郭亞琦，張召，等.動態(tài)系統(tǒng)的主動故障診斷技術(shù) [J].自動化學(xué)報，2020，46（8）：1557-1570.

[2] 王世林.基于多元統(tǒng)計分析的工業(yè)過程故障檢測研究 [D].北京：華北電力大學(xué)（北京），2018.

[3] 劉春燕.基于PCA的過程工業(yè)故障檢測 [D].綿陽：西南科技大學(xué)，2017.

[4] KU W，STORER R H，GEORGAKIS C. Disturbance detection and isolation by dynamic principal component analysis [J].Chemometrics and Intelligent Laboratory Systems，1995，30（1）：179-196.

[5] LI W H，QIN S J. Consistent dynamic PCA based on errors-in-variables subspace identification [J].Journal of Process Control，2001，11（6）：661-678.

[6] KIM D，LEE I B. Process monitoring based on probabilistic PCA [J].Chemometrics and Intelligent Laboratory Systems，2003，67（2）：109-123.

[7] 楊沛武.概率PCA多元統(tǒng)計方法在過程監(jiān)控中的應(yīng)用研究 [D].無錫：江南大學(xué)，2008.

[8] LI G，QIN S J，ZHOU D H. A New Method of Dynamic Latent-Variable Modeling for Process Monitoring [J].IEEE Transactions on Industrial Electronics，2014，61（11）：6438-6445.

[9] WEN Q J，GE Z Q，SONG Z H. Data-based linear Gaussian state-space model for dynamic process monitoring [J].AlChE Journal，2012，58（12）：3763-3776.

[10] 周樂.基于概率的工業(yè)過程數(shù)據(jù)建模與故障檢測 [D].杭州：浙江大學(xué)，2015.

[11] QIN J S. Statistical Process Monitoring：Basics and Beyond [J].Journal of Chemometrics，2003，17（8-9）：480-502.

[12] BOX G E P. Some Theorems on Quadratic Forms Applied in the Study of Analysis of Variance Problems，I. Effect of Inequality of Variance in the One-Way Clas-sification [J].The Annals of Mathematical Statistics，1954，25（2）：290-302.

[13] JAAKKOLA T S. Tutorial on variational approximation methods [M]//OPPER M，SAAD D. Advanced Mean Field Methods：Theory and Practice.Massachusetts：MIT Press，2001：129-159.

[14] JORDAN M I，GHAHRAMANI Z，JAAKKOLA T S，et al. An Introduction to Variational Methods for Graphical Models [J].Machine Learning，1999，37（2）：183-233.

[15] BEAL M J. Variational algorithms for approximate Bayesian inference [D].London：University College London，2003.

作者簡介：李吉?。?995—），男，漢族，河南鄧州人，碩士研究生在讀，研究方向：故障檢測、圖像處理;章智杰（1996—），男，漢族，江蘇連云港人，碩士研究生在讀，研究方向：圖像處

理;李其操（1997—），男，漢族，浙江諸暨人，碩士研究生在讀，研究方向：圖像處理;通訊作者：董自健（1973—），男，漢族，江蘇連云港人，教授，博士研究生，研究方向：檢測與控制、通信技術(shù)。