一種分布式CPS異常檢測的無監(jiān)督圖模型

張 錦，程良倫

廣東工業(yè)大學(xué) 自動化學(xué)院，廣州 510006

1 引言

得益于泛在感知、先進計算以及超強連接，現(xiàn)代分布式信息物理融合系統(tǒng)（CPS）如交通網(wǎng)絡(luò)[1]、電廠、電網(wǎng)[2]均表現(xiàn)有強大的高效率、魯棒性以及靈活性的潛力。這一事實在世界各地以及一般工業(yè)[3]的研究當(dāng)中都有提及。然而要實現(xiàn)這樣的潛力，針對CPS的有效建模及分析方法必須具備這樣的性質(zhì)：可擴展、魯棒性、靈活性以及適應(yīng)性。目前大多數(shù)分析方法是在很大程度上依賴于領(lǐng)域知識的需要細(xì)致校準(zhǔn)和驗證的基于規(guī)則的模型[4]。

從性能監(jiān)控和分布式CPS診斷的角度來看，技術(shù)上的挑戰(zhàn)來自于分布式CPS有著巨大數(shù)量的強偶合子系統(tǒng)[5]，同時這些子系統(tǒng)存在不同的操作模式。要對各個子系統(tǒng)分別建模來捕獲其復(fù)雜的交互作用非常的棘手，但是基于數(shù)據(jù)驅(qū)動的建模方式可以有效地緩解這樣的問題[6]。然而，大多數(shù)基于數(shù)據(jù)驅(qū)動的建模方式需要大量的正樣本（無故障）及所有可能的負(fù)樣本（如，物理故障或網(wǎng)絡(luò)異常），負(fù)樣本在現(xiàn)實系統(tǒng)的生命周期中是很難獲取到的。因此，異常檢測方法應(yīng)該有具備如下性質(zhì)：（1）在沒有真正出現(xiàn)故障時識別大多數(shù)操作模式的潛力；（2）具備無監(jiān)督學(xué)習(xí)能力從標(biāo)稱模式中識別異常模式。此外，物理空間產(chǎn)生的主要是從傳感器和執(zhí)行器而來的連續(xù)時間信息，而網(wǎng)絡(luò)空間主要是處理物理信息而產(chǎn)生離散事件驅(qū)動的數(shù)據(jù)。這種基本屬性的差別和信息本質(zhì)性質(zhì)的差別，使得目前大多數(shù)的方法在處理網(wǎng)絡(luò)空間和物理空間時，分別進行建模及分析（具體細(xì)節(jié)可以參考文獻[7]）。

在這種背景下，本文提出一種用于分布式CPS在整個系統(tǒng)上異常檢測的框架，該框架采用表征CPS各個子系統(tǒng)間相互作用的基于符號動力學(xué)的時空特征提取方案。符號動態(tài)過濾（SDF）是一種復(fù)雜系統(tǒng)的數(shù)據(jù)驅(qū)動建模方式，其采用一種稱為數(shù)據(jù)抽象的統(tǒng)一表示描述不同類型的數(shù)據(jù)，在描述不同類型的數(shù)據(jù)方面具有優(yōu)勢。數(shù)據(jù)抽象預(yù)處理及對相關(guān)變量（如，傳感器的時間序列）進行數(shù)據(jù)空間劃分有助于對網(wǎng)絡(luò)和物理子系統(tǒng)進行統(tǒng)一建模。由SDF捕獲的特征用于構(gòu)造時空模式網(wǎng)絡(luò)（STPN）[8]——最近提出的一種因果關(guān)系圖模型概念。通過訓(xùn)練這樣的模型，本文為檢測低概率事件或異常模式提出了一個推理方案。

2 背景及基本理論

2.1 時空模式網(wǎng)絡(luò)

符號動態(tài)過濾（SDF）已被證明從異常檢測和模式分類的時間序列數(shù)據(jù)中提取關(guān)鍵詞是非常有效的[9]。其核心思想是，從一個過程發(fā)出的符號序列（如，離散時間序列）可近似為D階馬爾可夫模型（D-Markov model），用于捕獲過程當(dāng)中潛在的關(guān)鍵行為。離散化或符號化的過程成為分區(qū)[9]?；诓煌哪繕?biāo)函數(shù)，一些研究當(dāng)中提出了不同的分區(qū)方法：如均勻分區(qū)（UP）、最大熵分區(qū)（MEP）、最大目標(biāo)離散化（MBD）等[10]。本文采用的是均勻分區(qū)方式。D階馬爾可夫模型本質(zhì)上是一個概率有限狀態(tài)自動機（PFSA），可以由狀態(tài)（表示數(shù)據(jù)空間的各個部分）以及狀態(tài)之間的轉(zhuǎn)換概率（可從數(shù)據(jù)中學(xué)習(xí)得到）描述。

時空模式網(wǎng)絡(luò)（STPN）可定義如下：

一個STPN可以由一個四元組表示W(wǎng)D=(Qa,Σb,Πab,Λab)，a,b表示STPN的節(jié)點。

（1）Qa={q1,q2,…,是與符號序列Sa相對應(yīng)的狀態(tài)集。

（2）是符號序列Sa的集合。

（3）Πab是一個| Qa|×| Σb|的符號生成矩陣，第ij個元素表示在符號序列中從狀態(tài)qi轉(zhuǎn)變時在符號序列中找到符號σi的概率。自身符號生成矩陣稱為原子模式（APs），當(dāng)a=b時；交叉符號生成矩陣稱為交互模式（RPs），當(dāng)a!=b時。

（4）Λab表示一個交互模式重要性指標(biāo)（或因果關(guān)系的程度），可以采用基于信息論的指標(biāo)，如傳遞熵或互信息來定義。

一個STPN可以用圖1描述。

圖1 子系統(tǒng)自身行為及子系統(tǒng)間交互行為的原子模式（APs）和交互模式（RPs）提取

2.2 受限玻爾茲曼機

受限玻爾茲曼機（RBM）近期在深度學(xué)習(xí)領(lǐng)域的無監(jiān)督特征提取中受到了很大的關(guān)注[11-12]。受限玻爾茲曼機的基本結(jié)構(gòu)如圖2的無監(jiān)督學(xué)習(xí)層（左上角）。作為基于能量的模型，通過學(xué)習(xí)權(quán)重和偏執(zhí)，使得系統(tǒng)在正常執(zhí)行期間觀察到的特征屬性獲得低能量（或高概率）?？紤]由一組可見變量v=(v1,v2,…,vD)和一組隱藏變量h=(h1,h2,…,hF)描述的系統(tǒng)狀態(tài)。變量的值根據(jù)實際需求可以是二值變量或?qū)崝?shù)變量。這些變量確定了系統(tǒng)的一個特定狀態(tài)和與其相關(guān)的能量值E(v,h)。能量值是變量之間連接權(quán)重的函數(shù)（對于RBM內(nèi)部連續(xù)僅考慮可見變量，不考慮隱藏變量），偏執(zhí)則是變量的函數(shù)。

圖2 一個基于數(shù)據(jù)驅(qū)動的分布式CPS系統(tǒng)異常檢測框架

通過以上描述可得到，狀態(tài) p(v,h)的概率僅依賴于連接(v,h)的能量，同時滿足玻爾茲曼分布：

在訓(xùn)練過程當(dāng)中，通常最大化訓(xùn)練數(shù)據(jù)的似然函數(shù)來得到權(quán)重和偏置。

3 提出模型

圖2描述了以上提出的基于數(shù)據(jù)驅(qū)動的系統(tǒng)異常檢測框架。在訓(xùn)練階段，學(xué)習(xí)STPN+RBM模型的步驟如下：

（1）從多元訓(xùn)練符號序列中學(xué)習(xí)APs和RPs（單節(jié)點行為和節(jié)點對的交互行為）。

（2）考慮來自訓(xùn)練樣本的短符號子序列，同時對其進行評估求得Λij，i和 j對應(yīng)每一個子序列。

（3）對于每個子序列，基于在Λij中用戶設(shè)定的閾值，為每個AP和RP分配狀態(tài)0或1；從而每個子序列轉(zhuǎn)換成一個長為L的二值向量，L=#AP+#RP。

（4）然后使用RBM對APs和RPs對應(yīng)的可見層的節(jié)點進行系統(tǒng)行為建模。

（5）將從子序列生成的二值向量作為特征進行訓(xùn)練RBM模型。

3.1 訓(xùn)練STPN+RBM模型

考慮訓(xùn)練一個多元時間序列（上文提到的標(biāo)稱操作數(shù)據(jù)），X={XA(t)},t∈N,A=1,2,…,f,f是時間序列變量的個數(shù)或維數(shù)。首先，執(zhí)行符號化及學(xué)習(xí)PFSA來抽取原子模式及交互模式，將其用于對應(yīng)的 f個頂點和 f2條邊的時空模式網(wǎng)絡(luò)（STPN）。在這種情形下，定義符號序列集為S={SA}。然后，定義一個短序列X={XA(t)},t∈N*,A=1,2,…,f，其中 N*是 N 的子集。本質(zhì)上講，從整體訓(xùn)練數(shù)據(jù)中按不同的時窗（由N*表示）抽出的序列均可視為短序列。類似之前的定義，通過不同時間窗抽取的一組符號子序列可以定義為S={SA}。

下一步是為從整個時間序列提取的每個短子序列計算Λij。雖然基于信息論的度量是一個非常好的選擇，但這種度量方式需要大量的樣本點，因此在異常檢測這種缺乏負(fù)樣本的場景當(dāng)中不太適用。本文中采用文獻[13-14]中提到的統(tǒng)計推理策略度量方式，該度量的計算步驟分兩個階段：建模過程、預(yù)測過程。

建模過程：對整個訓(xùn)練數(shù)據(jù)的一個STPN可以描述為WD=(Qa,Σb,Πab,Λab)。整個符號序列集合定義為S；整個序列的狀態(tài)表示產(chǎn)為Q={Qa,a=1,2,…,f}；一個模式的狀態(tài)Πab取決于符號序列Qa及序列的狀態(tài)Sb。在該步驟，Πab的每一行都初始化為一個隨機向量。對于第m行，隨機向量的先驗概率密度函數(shù)| {Qa,Sb}由聯(lián)合狀態(tài)符號序列{Qa,Sb}決定，并且服從狄利克雷分布。概率密度函數(shù)如下：

其中是隨機向量的一個實現(xiàn)：

其中

符號序列的時間長度；是在狀態(tài)之后估計的。

Qa(k)是狀態(tài)序列Qa的第k個狀態(tài)；Sb(k +1)是符號序列Sb的第k+1個符號。

從而公式（3）可以改寫為：

其中T(n)=(n-1)!。

根據(jù)馬爾可夫的屬性可知的行向量之間統(tǒng)計意義上相互獨立。通過公式（2）、（5）可以求得概率轉(zhuǎn)換矩陣關(guān)于狀態(tài)符號序列{Qa,Sb}的條件聯(lián)合概率密函數(shù)：

其中

預(yù)測過程：在對整個訓(xùn)練序列集建完模之后，預(yù)測的主要目的是計算給定的子序列（由Q和S描述）的度量Λab(Q)。度量的值大小表示模式Πab的重要度或短子序列a→b的因果關(guān)系。根據(jù)這一觀點，可以用公式表示如下：

當(dāng)轉(zhuǎn)換矩陣是已知的情況下，聯(lián)合狀態(tài)符號子序列的概率是獨立多項式分布的乘積：

的定義與上文對子序列的定義一致。

做與文獻[13]中類似的推導(dǎo)，可以得到度量矩陣：

其中K是一個常量系數(shù)。

根據(jù)公式（9）可以對相應(yīng)短子序列獲得 APs和RPs重要性指標(biāo)。為了訓(xùn)練整個系統(tǒng)的RBM模型，可進一步將APs和RPs的度量指標(biāo)標(biāo)準(zhǔn)化為二值狀態(tài)（低值置0，高值置1）。注意到，每個子序列所有的APs和RPs一起形成長度為L=f2的二值向量(L=#AP+#RP，其中：#AP=f,#RP=f*(f-1))。一個這樣的二值向量就是一條用于訓(xùn)練系統(tǒng)級RBM的樣本（有f2個可見單元），從不同子序列抽取的許多這樣的樣本則構(gòu)成了整個訓(xùn)練樣本集。然后如圖2所示，使用最大似然估計來訓(xùn)練RBM。本文采用的方式是將度量矩陣強制轉(zhuǎn)換成二值向量來簡化RBM的訓(xùn)練，對于訓(xùn)練來說，這一過程是非必須的。

說明：在STPN學(xué)習(xí)過程中，加入一個滯后時間，STPN+RBM可以處理可變時滯系統(tǒng)，同時會大降低整個學(xué)習(xí)過程的復(fù)雜性。這與遞歸神經(jīng)網(wǎng)絡(luò)（RNNs）類似。

3.2 異常檢測過程

異常檢測過程采用RBM——一個基于能量的概率圖模型的自由能概念。RBM的能量函數(shù)可以定義如下：

W是隱藏單元的權(quán)重，b和c分別是可見單元和隱藏單元的偏置。

利用RBM的權(quán)重和偏置，可以獲得自由能的表示：

自由能的另一種表達可以是[8]：

在訓(xùn)練階段，最小化能量函數(shù)則可以得到模型的權(quán)重和偏置。通常異常模式應(yīng)該表現(xiàn)為低概率（高能量）的狀態(tài)。因此，在測試階段，可以根據(jù)模型對樣本的概率估計來檢測異常模式。在做異常檢測的過程中，測試子序列按照訓(xùn)練數(shù)據(jù)的轉(zhuǎn)換方式，將其轉(zhuǎn)換為 f2維二值向量。對多個測試子序列進行估計，則可以得到自收能的分布情況。對于標(biāo)稱模式，自由能的分布應(yīng)該與訓(xùn)練數(shù)據(jù)的分布比較接近，對于異常模式，自由能的分布將不同于訓(xùn)練數(shù)據(jù)的分布。

對比訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)的分布，可以采用KL散度這一指標(biāo)[15]。因為KL散度是一個非對稱的描述兩個分布P和Q的指標(biāo)，是一種非對稱的信息度量。因此這里可以采用一種對稱KL距離的度量，即：

4 實驗與討論

4.1 數(shù)據(jù)生成

向量自回歸（VAR）是處理多個相關(guān)指標(biāo)分析與預(yù)測最易操作的模型之一，在多元時間序列數(shù)據(jù)中非常的靈活與簡單，在經(jīng)濟學(xué)等學(xué)科中得到了廣泛的應(yīng)用[16]。向量自回歸的基本模型(Y(t)=yi,j,i=(1,2,…,f),t∈N)可以用如下公式表示：其中 p是時間滯后階數(shù)，Ai,j是第 j個時間序列對第i個時間序列的影響系數(shù)，μt是均值E()μt=0；協(xié)方差矩陣E( )

μtμt=Σμ的噪聲。在本文中，將使用VAR模型生成多變量數(shù)據(jù)用于異常檢測模型的仿真。5個頂點的層次結(jié)構(gòu)表示各個頂點之間不同的相互作用。各種相互作用表示分布式CPS中的標(biāo)稱條件和異常條件。生成兩組數(shù)據(jù)用于兩種情況的分析工作：（1）定義6種模式，第1種屬于標(biāo)稱模式，其余5種為異常模式；（2）定義8種模式，前3種屬于標(biāo)稱模式，其余5種為異常模式。第2種情況用于模擬CPS中多個標(biāo)稱操作模式的情形。

然后從VAR生成的代表不同交互模式的原始數(shù)據(jù)中訓(xùn)練得到模型STPN。對于無監(jiān)督RBM的訓(xùn)練過程，僅僅從標(biāo)稱模式中學(xué)習(xí)權(quán)重和偏置。訓(xùn)練好的RBM模型將用來對所有操作模式的數(shù)據(jù)計算關(guān)于不同模式輸入的自由能。最終對多個服從高斯分布的測試數(shù)據(jù)進行自由能的評估，對比評估結(jié)果。

4.2 實驗一分析：單標(biāo)稱模式

圖3中表示了預(yù)定義的圖模型，異常發(fā)生的情況是一個頂點出錯并且與其連接的邊丟失。在實驗中，采用標(biāo)稱條件的原子模式和交互模式數(shù)據(jù)訓(xùn)練RBM模型。然后將訓(xùn)練好的RBM模型用于6個模式的所有測試數(shù)據(jù)，產(chǎn)生各個自由能的概率分布，如圖4。從圖中可以看出，標(biāo)稱模式的自由能分布與訓(xùn)練數(shù)據(jù)自由能分布非常接近，其他模式的自由能分布則非常的不同。這表明異常模式具有高自由能，同時發(fā)生的概率低。計算分布之間的相關(guān)性得到，訓(xùn)練/標(biāo)稱模式的自由能分布與所有其他測試模式的自由能分布之間的KLD指標(biāo)分別為：

圖3 合成數(shù)據(jù)的圖模型

圖4 合成數(shù)據(jù)的圖模型

4.3 實驗二分析：多標(biāo)稱模式

與先前的情況類似，預(yù)定義的圖模型如圖5，自由能的分布如圖6。與預(yù)期的一致前三種模式的自由能分布類似于訓(xùn)練數(shù)據(jù)的自由能分布。KLD分別為0.052、0.263、0.639。該圖還顯示，對于4至8，KLD值分別為5.06、5.46、6.24、4.91、9.97可以清楚地鑒定異常模式。總體上，結(jié)果清晰地表明，所提出的框架可以捕獲在一個模型內(nèi)的多個標(biāo)稱模式，同時還可以有效地檢測因果關(guān)系模式的輕微變化。

圖5 RBM模型對各種模式（單標(biāo)稱模式）的自由能分布

圖6 RBM模型對各種模式（多標(biāo)稱模式）的自由能分布

4.4 討論

在分布式CPS當(dāng)中，其異常機制、特點和持續(xù)時間各不相同，這使得異常檢測非常困難，特別是對所有可能的異常標(biāo)記數(shù)據(jù)的收集。本文提出的框架中只需要正例樣本，異常檢測則看作是以標(biāo)稱數(shù)據(jù)為條件的低概率事件。從實驗的結(jié)果可以看出，異常條件下的自由能分布不同于標(biāo)稱條件下自由能的分布，可以使用諸如KLD的變化來度量。此外，它有潛力監(jiān)控從小的物理退化到嚴(yán)重故障或網(wǎng)絡(luò)攻擊，其中小KLD表示因果模式的輕微變化，可預(yù)示早期退化或故障前兆。

5 結(jié)束語

本文提出了一個系統(tǒng)級的CPS異常檢測框架，該框架是一種基于數(shù)據(jù)驅(qū)動的方法，但不同于監(jiān)督學(xué)習(xí)方法，不需要標(biāo)記好的正負(fù)樣本。該框架包括了一個時空特征提取方法用于發(fā)現(xiàn)和表征CPS各子系統(tǒng)之間相互作用的因果關(guān)系，以及使用RBM的全系統(tǒng)模式的自由能估計。以上實驗表明，本文所提出的方法可以在單個概率圖模型中捕獲多個不同的標(biāo)稱模式，并通過識別低概率事件來檢測異常模式。同時準(zhǔn)確性和魯棒性方面也表現(xiàn)良好。

目前的工作重點主要集中在檢測異常情況以及量化假警報，進一步的工作將擴展到：（1）將圖模型應(yīng)用到各種異常情況的根本原因分析上；（2）采用疊加RBM的方法捕獲更復(fù)雜的標(biāo)稱模式；（3）檢測分布式CPS中同時出現(xiàn)的多個故障。

[1]Work D，Bayen A.Impacts of the mobile internet on transportation cyber-physical systems：Traffic monitoring using smartphones[C]//National Workshop for Research on High-Confidence Transportation Cyber-Physical Systems：Automotive，Aviation，and Rail，2008：18-29.

[2]趙慶周，李勇，田世明，等.基于智能配電網(wǎng)大數(shù)據(jù)分析的狀態(tài)監(jiān)測與故障處理方法[J].電網(wǎng)技術(shù)，2016，40（3）：774-780.

[3]Bradley J，Barbier J，Handler D.Embracing the Internet of everything to capture your share of$14.4 trillion[D].CISCO White Paper，2013.

[4]文成林，呂菲亞，包哲靜，等.基于數(shù)據(jù)驅(qū)動的微小故障診斷方法綜述[J].自動化學(xué)報，2016，42（9）：1285-1299.

[5]Sanislav T，Miclea L.Cyber-physical systems concept，challenges and research areas[J].Control Engineering&Applied Informatics，2012，14（2）：28-33.

[6]Choi A，Zheng L，Darwiche A，et al.A tutorial on Bayesian networks for system health management[J].Machine Learning and Knowledge Discovery for Engineering Systems Health Management，2011，10（1）：1-29.

[7]Krishnamurthy S，Sarkar S.Scalable anomaly detection and isolation in Cyber-physical systems using Bayesian networks[C]//Proceedings of ASME Dynamical Systems and Control Conference，2014.

[8]Jiang Z，Sarkar S.Understanding wind turbine interactions using spatiotemporal pattern network[C]//Proceedings of ASME Dynamical Systems and Control Conference，2015.

[9]Rao C，Sarkar S，Ray A，et al.Comparative evaluation of symbolic dynamic filtering for detection of anomaly patterns[J].Signal，Image and Video Processing，2009，3（2）：101-114.

[10]Sarkar S，Srivastav A，Shashanka M.Maximally bijective discretization for data-driven modeling of complex systems[C]//American Control Conference，2013：2674-2679.

[11]Hinton G，Salakhutdinov R.Supporting online material for reducing the dimensionality of data with neural networks[J].Science，2006，313：504-507.

[12]Roux N，Bengio Y.Representational power of restricted boltzmann machines and deep belief networks[J].Neural Computation，2008，20（6）：1631-1649.

[13]Akintayo A，Sarkar S.A symbolic dynamic filtering approach to unsupervised hierarchical feature extraction from time Series data[C]//American Control Conference，2015：5824-5829.

[14]Sarkar S，Mukherjee K，Sarkar S，et al.Symbolic dynamic analysis of transient time series for fault detection in gas turbine engines[J].Journal of Dynamic Systems Measurement&Control，2012，135（1）.

[15]Kullback S，Leibler R A.On information and sufficiency[J].The Annals of Mathematical Statistics，1951，22（22）：79-86.

[16]Goebel R，Roebroeck A，Kim D，et al.Investigating directed cortical interactions in time-resolved FMRI data using vector autoregressive modeling and Granger causality mapping[J].Magnetic Resonance Imaging，2003，21（10）：1251-1261.