基于因果關系的故障傳播路徑辨識方法研究

呂佳朋, 史賢俊, 秦 亮, 趙超輪

(海軍航空大學岸防兵學院, 山東 煙臺 264001)

0 引 言

隨著科技的進步,裝備系統正向著大型化、復雜化、高耦合化的方向發展,這使得系統中各個變量之間的影響變得更加復雜。當系統發生故障時,故障往往會沿著系統結構,使得系統中多個檢測變量偏離正常狀態而發生報警,此時如果不能正確辨識出故障的源頭,正確隔離出故障發生的部件,那么便無法及時消除報警,使系統恢復正常狀態。因此,如何識別系統中故障傳播的途徑,從多個故障變量中確定故障發生的根源,對于系統裝備的警報消除及故障排除具有重要的意義[1-2]。

故障傳播路徑辨識,部分文獻稱其為故障溯源[3-4]或故障根源診斷[5-7],其可以找到受故障影響的變量,識別故障傳播路徑,定位故障源,為現場保障人員提供必要的操作參考。目前開展故障路徑識別的步驟主要包括如下3個[8]:① 故障相關變量辨識,即當故障發生時,辨識出的能夠受到故障影響的變量集合;② 相關變量因果辨識,即對于受故障影響的相關變量集合,研究變量之間的因果關系;③ 故障傳播路徑繪制,即繪制因果圖,展現變量之間的影響關系。目前,關于故障傳播路徑辨識的相關研究也主要從上述3個方面進行展開。

在故障相關變量辨識方面,目前生產實踐中最常用的方法是基于貢獻圖或重構貢獻圖的方法[9]。該方法的主要思想是設計出合理的變量貢獻指標,并認為貢獻指標越大,變量受到的故障影響也就越大,為故障相關變量。該方面的研究主要集中于對變量貢獻指標的研究上,如變量貢獻度指標[10]、基于貝葉斯的魯邦高斯混合貢獻指標、平方預測誤差(squared prediction error, SPE)、T2統計量[1]等指標。但是,這些方法存在的問題是容易受到涂抹效應的影響,該效應使得部分未發生故障的變量的貢獻值也超過閾值,從而出現錯誤的診斷結論。

在相關變量因果辨識方面,最具代表性的方法包括格蘭杰因果分析[11-13]、傳遞熵[14-15]等。何飛等[16]提出對比格蘭杰因果分析方法確定故障信息引起的異常因果關系。文獻[17]提出一種基于傳遞熵和修正條件互信息的方法,并且使用傳遞熵方法構建因果圖來實現故障定位。文獻[18]解決了傳遞熵在進行因果辨識時,由于變量中存在中間變量或者共因變量所導致的冗余變量問題。文獻[19]證明,對于高斯變量而言,傳遞熵及格蘭杰因果是等價的。該類方法存在的問題主要為:① 格蘭杰因果僅能適用于兩變量、平穩、線性的時間序列,同時格蘭杰因果更多展現的是格蘭杰原因事件對于格蘭杰結果事件的預測作用,而非日常所說的因果關系,從這個角度來看格蘭杰因果作為因果判據并不充分;② 作為一種信息論領域揭示因果關系的算法,當所需變量的個數過多時,傳遞熵巨大的計算量和計算成本會影響因果關系挖掘的效率[20]。文獻[21]指出,傳遞熵可能導致虛假因果的問題,同時使用傳遞熵衡量因果關系存在不一致性問題,即因果關系的方向可能會隨著系統狀態的變化而發生變化[22]。

在因果圖繪制方面,目前大部分文獻是將含有因果關系的變量直接用圖論中邊的概念表示出來,這種表示方法過于機械,往往會導致圖中存在過多的冗余連接而使圖的結構復雜,不易于清晰展現變量之間的關系。

針對上述分析中存在的問題,本文將因果關系的概念[23-24]引入到對故障傳播路徑的辨識當中,從系統的歷史數據出發,挖掘系統中存在的因果關系,提出一種基于因果關系的故障傳播路徑辨識方法。該方法首先利用系統中故障發生的因果性,確定故障發生時受影響的變量,構建故障相關變量集合;其次,通過因果關系指標判斷不同變量之間的因果關系,構建了因果關系矩陣;最后,基于圖論理論,提出了?？蛇_性的賦權有向圖最小生成樹算法,該算法可以圖形化表示因果關系矩陣的內容,并去除其中的冗余連接,從而展現故障傳播路徑,有利于故障發生源頭的判斷。

本文結構如下:第1節分析了因果關系三要素及因果關系的數據表示方法,并從因果關系的視角分析了故障的發生和傳播;第2節詳細描述了本文所提基于因果關系的故障傳播路徑辨識方法的具體流程;第3節通過一個電路的實例驗證了本文方法的合理性和優越性。

1 因果關系視角下的系統故障分析

本節首先介紹因果關系中的相關概念,然后從因果關系的角度描述系統中故障發生的過程。

1.1 因果關系

因果關系是自然界中普遍存在的一種自然關系,一直以來,學者們對于因果關系的概念爭論不休,因此給出因果關系嚴格的定義是非常困難的,但是一般有定義如下。

定義 1(因果關系) 因果關系是“因”事件和“果”事件之間客觀存在的關系,其中“因”事件是導致“果”事件發生的原因。

傳統意義上對于因果關系的確定主要通過基于實證的方法,即通過實施隨機對照實驗來確定變量之間的因果關系,但這種方法容易受到倫理限制、個體不依從、費用時間花銷等因素的影響,因而在某些情況下具有不可操作性[25]。隨著大數據及計算機算力的發展,越來越多的學者開始嘗試從已知的大量經驗(實驗或觀察數據)中發現事物之間的因果關系[26-28]。

為了能夠從數據的角度處理因果關系的問題,需要一種能夠公式化表達因果關系的方法。Pearl[23-24]提出了結構因果模型(structural causal model, SCM)來描述數據之間的因果機制。

定義 2(函數因果模型) 在一般化的形式中,一個SCM由以下一組方程構成:

xi=fi(pai,ui),i=1,2,…,n

(1)

式中:pai表示變量xi的父變量的集合;ui表示由遺漏因子所導致的誤差。

一個完整的因果關系一般涉及3個要素:原因變量(式(1)中的pai),結果變量(式(1)中的xi)及因果機制(式(1)中的fi)。而式(1)說明,從數據的角度看因果關系,因果關系實質上是一種將原因變量映射為結果變量的映射。

1.2 故障發生及其中的因果性分析

一個簡易的信號與模型的關系如圖1所示,該圖表明了系統中模塊和信號之間存在的關系,而一個系統無論多么復雜(見圖2),其內部均是由許多大大小小類似的模塊構成。

圖2 系統結構示意圖Fig.2 Schematic diagram of system structure

信號在系統中進行傳遞,經過模塊后即可變為另一種信號。該過程使用頻域進行表示:

Y(S)=G(S)·X(S)

(2)

式中:S表示對時域函數進行拉普拉斯變換后得到的新函數的自變量。

式(2)可以看作是經過G(S)的映射作用將X(S)映射為Y(S)的過程,通過對比可以認為輸入信號經過模塊的映射作用映射為輸出信號,該過程類似于因果機制將原因變量映射為結果變量的過程。其中,輸入信號為因變量,輸出變量為果變量,模塊為因果映射機制。

信號經過模塊傳遞后其因果關系的示意圖如圖3所示。

圖3 模塊信號和因果關系對應示意圖Fig.3 Schematic diagram of model signal and causality

對于一個模塊或者系統而言,其輸入信號和輸出信號之間是存在因果關系的。如圖1所示的信號傳遞過程,該過程采用因果關系的SCM可表示為

Y=G(X)

式中:X表示模塊的輸入信號;Y表示模塊的輸出信號;G表示因果映射機制。

如前所述,在信號傳遞的過程中,系統模塊從因果關系的角度來看實質上相當于因果映射機制G,而系統模塊并不隨著外界信號的變化而發生改變,因此G僅僅與系統的結構或結構相關參數有關。故障是一種系統內部由器件老化或損壞等原因導致的系統不能執行規定功能的狀態。故障的發生會導致系統的結構或功能參數發生變化,這種改變使得G發生變化,最終使得輸出信號發生畸變。故障發生時,系統的傳遞過程變為

Y′(S)=G′(S)·X(S)

(3)

因此,從因果關系的角度來看,可以認為模塊作為一種因果機制,其作用機理是由模塊內部的結構所決定的。當故障發生時,模塊內部結構發生變化,導致其因果機制發生變化,因此模塊的輸出信號發生畸變。從因果視角下看系統故障的發生,即因果機制發生了變化。

故障發生時,系統的SCM變為

Y=G′(X)

式中:G′表示故障發生后該模塊的映射機制。

2 基于因果關系推斷的故障傳播路徑辨識方法

根據引言所述,為了對故障傳播路徑進行辨識,首先,需要對故障相關變量進行辨識,找出受故障影響的變量集合;其次,是要對受故障影響的觀測量之間的因果關系進行挖掘,辨識出變量之間的因果關系;最后,以圖的方式對故障的傳播路徑進行展現。

為使行文流暢,現不加證明地給出本文所用的一些定理和假設。

假設 1(原因變量及因果機制獨立性) 如果在因果機制G的作用下,X是Y的原因變量,則因果映射機制G與原因變量X的分布是相互獨立的。

定理 1基于SCM的因果關系表達式為

Y=G(X)

該表達式可以改寫為

D(PY‖εY)=D(PX‖εX)+D(uG‖εY)

式中:PX和PY分別表示變量X和變量Y的概率密度函數;D(·‖·)表示兩概率密度分布的相對熵距離;u表示PX在εX上的投影;uG是u在因果映射機制G下的像;εX和εY表示X和Y在光滑指數族分布ε上的參考分布,即

定理1的相關推導過程可以參照文獻[29],此處不再敷述。定理1表明,因果機制G將X映射為Y,Y的不規則度不僅取決于映射G,還與X的不規則度有關,既Y的不規則度等于X的不規則度疊加映射G的不規則度[30]。定理1揭示了在因果關系中,原因變量與結果變量之間存在的不規則度上的不對稱性,這種不對稱性可以作為推斷因果關系方向的依據。

2.1 故障相關變量辨識方法

假設系統的狀態集合為F={F0,F1,…,Fm},其中m表示系統故障狀態的數目,F0表示正常狀態;系統所要測量得到的變量集合為S,S={Sin,S1,…,Sn},其中Sin表示整個系統的輸入信號,n表示除輸入信號系統所要測量的信號的個數。

如第1.2節所述,對于同一個模塊的輸入信號和輸出信號之間存在的因果關系,可以用SCM進行表達。但是根據假設,系統的故障模式集合F和系統的信號集合S是已知的,而系統構型等內部的結構信息對于本文而言是未知的,即并不能判斷信號集合S中的信號是否存在因果關系。圖2所示的結構中,測試點5的信號與測試點7的信號之間存在因果關系,但是測試點5的信號和測試點2的信號之間卻不存在因果關系。

需注意到,第1.2節所述的模塊不一定指的是一個實際存在的系統,模塊可以通過多個模塊進行復合而形成新的模塊,如在圖2中,模塊3、模塊4及模塊5可以進行模塊的復合。復合以后的模塊依舊滿足第1.2節中的論述,因此對于復合系統而言,其輸入信號和輸出信號也存在著因果關系,如在圖2中,模塊3、模塊4及模塊5復合所得的模塊保證了測試點2和測試點3之間存在因果關系。進一步地,在系統中,各個模塊都可以通過一定的方式進行復合,因此在互為因果的輸入信號和輸出信號之間形成一個虛擬的模塊。而在整個系統之中,均可以建立從輸入信號到可觀信號之間的通路。在集合S中,輸入信號Sin可以看作是量測信號的原因變量,因此在集合S中可以寫出SCM方程組

(4)

根據定理1,當j=0時,式(4)可改寫為

(5)

式中:PSin和PSi分別表示信號Sin和Si的概率密度;εSi和εSin分別表示Si和Sin在光滑指數族分布上的參考分布;u表示PSin在εSin上的投影。

當故障Fj發生時,即j≠0時,式(4)可改寫為

(6)

將式(5)與式(6)兩式作差,可以得出

(7)

(8)

算法 1 故障相關變量篩選算法輸入:Q j=[x1,x2,…,xm]。輸出:故障相關變量集合S i={S i1,S i2,…,S ip},p表示故障相關變量的個數。初始化:初始化集合C1和C2,C1=?,C2=?。步驟 1 從Q j中任意選取兩個點,并將其作為均值點c1和c2。步驟 2 C1={c1},C2={c2}。步驟 3 對于Q j中剩下的任意點xi,分別求其與均值點c1和c2之間的距離|xi-cj|,j=1,2。步驟 4 根據最近距離,確定xi的歸屬集合,λi=argmin dij。步驟 5 Cλi=Cλi∪{xi}。步驟 6 更新均值點ci=1|ci|∑x∈cix,i=1,2。步驟 7 對于Q j中的任意點xi,分別求其與均值點c1和c2之間的距離|xi-cj|,j=1,2。步驟 8 根據最近距離,確定xi的歸屬集合,λi=argmin dij。步驟 9 Cλi=Cλi∪{xi}。步驟 10 重復步驟6～步驟9,直到兩個集合中的均值不再發生變化。步驟 11 取集合C1和C2中元素平均數值大的集合,其即為故障相關變量S i。

2.2 相關變量因果辨識方法

第2.1節確定了故障相關變量,本節對相關變量之間的因果辨識方法進行了研究。

根據定理1的描述,原因變量和結果變量之間的不規則度存在著不對等關系,利用這種不對等關系,可以用來指示因果變量之間的方向。構建因果關系指示指標來進行因果關系分析:

(9)

通過式(9),可以獲得因果矩陣,來表示變量之間的因果關系:

(10)

2.3 概率密度的微分熵估計方法

在對式(7)及式(9)進行計算時,需要知道信號的概率密度,但實際上概率分布是很難精確得知的,只能利用有限的樣本進行估算。因此,本節對式(7)及式(9)的估算進行討論。

不失一般性地,式(7)及式(9)均可以歸納成如下問題:設PX和PY分別表示變量X和Y的概率密度函數,D(·‖·)表示兩概率密度分布的相對熵,εX和εY表示變量X和Y在光滑指數族分布ε上的參考分布,求D(PY‖εY)-D(PX‖εX)。

引理 1[29]設u和v分別是PX和PY在εX和εY上的投影,則有

D(PX‖εX)=D(PX‖u)=-S(PX)+S(u)

式中:S(·)表示隨機變量的微分熵。

引理1的相關推導過程可以參照文獻[29],此處不再敷述。

根據引理1,D(PY‖εY)-D(PX‖εX)可以變形為

D(PY‖εY)-D(PX‖εX)= (-S(PX)+S(u))-(-S(PY)+S(v))= (S(u)-S(v))-(S(PX)-S(PY))

式中:S(u)-S(v)對于最終結果不產生影響[29];S(PX)-S(PY)使用熵估計器進行估計:

2.4 變量因果圖構建方法

根據式(10),進一步構建故障相關變量因果傳遞圖。

在一般的因果圖構建的過程中,均沒有考慮變量之間的因果變量是之間因果還是間接因果,只是機械地將因果矩陣表達的因果關系使用圖的方式表達了出來,這種做法會使得因果圖變得復雜繁瑣,特別是在涉及得變量很多時,不容易揭示變量之間的因果關系。

2.4.1 問題描述

不失一般性地,將因果圖構建的問題闡述如下。

由式(10)構建有向圖G(V,E),其中V為節點集合,代表故障相關變量,E為有向邊集合,代表變量之間的因果傳遞關系,且每條邊上有一個數作為權。在本文中,由于式(10)中的數字大小并不表示因果關系的強弱,僅僅定性地表明因果關系,故設有向邊的權值為1?，F需要求圖G的子圖T,T滿足如下條件:①T需要包含G中所有的頂點;②T需要具有最小的邊數;③T需要保持式(10)中的因果結構,即圖G(V,E)中任意兩個頂點之間的可達性是不變的。

2.4.2 基于?？蛇_性的賦權有向圖最小生成樹算法的因果圖構建方法

分析上述問題,從T的條件中可以看出,需要在賦權有向圖G中求得其最小生成樹,并且滿足頂點的可達性不變。本文根據文獻[31]提出一種保可達性的賦權有向圖最小生成樹算法,如算法2所示。

算法 2 ?？蛇_性的賦權有向圖最小生成樹算法輸入 因果矩陣C。輸出 有向圖T。步驟 1 根據因果矩陣C構建有向圖G(V,E),V為節點集合,E為有向邊集合。步驟 2 設子圖T(U,D)為所求最小生成樹,U為生成樹的頂點集合,D為生成樹的有向邊集合,初始化U和D為空集合。步驟 3 在圖G的集合V中,任意選取一個頂點v0放入U中,此時U={v0},V=V-{v0}。步驟 4 在V中,選擇頂點v1,要求:(1) v0與v1之間存在有向邊連接;(2) 在圖G中,v0與v1僅具有一條通路。將頂點v1加入到集合U中,將有向邊加入到集合D中,同時更新集合V。步驟 5 在V中,選擇頂點v2,要求:(1) v2與集合U中的點存在有向邊連接;(2) 該有向邊連接的兩個頂點在原圖G中僅存在一條通路;(3) 該有向邊的弧頭在子圖T中的入度為0(若弧頭不屬于子圖T的頂點集合,則默認弧頭的入度為0)。將頂點v2加入到集合U中,將有向邊加入到集合D中,同時更新集合V。步驟 6 重復步驟5,直到U=V為止。

2.5 算法流程

本文所提的基于因果關系的故障傳播路徑辨識算法示意圖如圖4所示。

圖4 基于因果關系的故障傳播路徑辨識算法示意圖Fig.4 Schematic diagram of fault propagation path identification algorithm based on causality

3 實例分析

3.1 實驗設置

本文選擇雙帶通濾波器作為研究對象對本文的方法進行驗證。雙帶通濾波器的電路結構圖如圖5所示,其中各元件的參數參照文獻[32],電路故障參數設置如表1所示。

表1 電路故障參數

圖5 雙帶通濾波器電路圖Fig.5 Circuit diagram of double bandpass filter

在進行實驗驗證時,按照步驟分別進行。首先進行故障相關變量辨識,采用本文方法及貢獻圖的方法進行對比;其次對相關變量因果關系進行辨識,采用本文方法及基于傳遞熵的方法進行對比;最后進行因果圖的繪制。

3.2 實驗結果

3.2.1 故障相關變量辨識結果

(1) 本文方法辨識結果

根據第2.1節所述的基于因果關系確定故障相關變量的方法,獲得系統不同故障狀態下的QFj,如表2所示。

表2 系統因果機制變化衡量結果

根據表2的結果進一步可以獲得故障相關變量辨識的結果,如表3所示。

表3 故障相關變量辨識結果

(2) 貢獻圖方法辨識結果

根據基于貢獻圖的方法,計算出各狀態下各測試點的貢獻圖,如圖6所示。

圖6 不同故障狀態的貢獻圖計算結果Fig.6 Contribution graph calculation results of different fault status

根據圖6結果確定的故障相關變量辨識的結果如表3所示。

3.2.2 相關變量因果關系辨識結果

根據表3所示的結果,進一步判斷各個相關變量之間的因果關系,由于故障狀態3及4僅有一個故障相關測試點,故選擇故障狀態1和2進行因果關系辨識。

(1) 本文方法辨識結果

通過第2.2節所述,對相關變量的因果關系進行識別,構建的因果矩陣分別如表4及表5所示。

表4 因果關系度量(故障狀態1,本文方法)

表5 因果關系度量(故障狀態2,本文方法)

(2) 基于傳遞熵的辨識結果

使用傳遞熵(在計算傳遞熵中使用分布直方圖的方式進行概率密度估算)構建的因果矩陣結果如表6及表7所示。

表6 因果關系度量(故障狀態1,基于傳遞熵的方法)

表7 因果關系度量(故障狀態2,基于傳遞熵的方法)

兩種方法在時間方面的開銷如圖7所示,圖中縱坐標軸使用的是對數坐標;橫坐標軸測試點對指的是在系統發生故障1時,測試點對1、2、4、5任意兩兩組合所組成的6對測試點對。

圖7 兩種算法的時間開銷Fig.7 Time cost of two algorithms

3.2.3 因果圖繪制結果

根據表4和表5的因果度量結果,繪制因果圖,展現故障傳播路徑,如圖8所示。

圖8 故障傳播路徑Fig.8 Propagation paths of fault

3.3 實驗分析

3.3.1 關于故障相關變量辨識結果的討論

從表3的結果來看,本文方法能對故障相關變量進行辨識,但是基于貢獻圖的方法對于部分故障未能給出辨識結果,同時基于貢獻圖的方法在部分變量辨識上,兩種統計指標的貢獻率相差較大。

下面通過對系統進行分析來說明辨識結果的正確性。

故障3及故障4均為電阻R18發生故障,由于故障發生在系統靠近末端的位置,所以僅測試點5能測出故障,而前置測試點1～4均未發生明顯變化。故障2為電阻R8發生故障,當故障發生時,測試點2的信號會發生變化,而與之有物理連接的測試點4及測試點5也會受到影響。故障的傳播及影響需要依托系統的拓撲結構,模塊之間的物理連接為故障影響的傳播提供了路徑。

但是,物理連接并不能保證故障的影響一定會傳遞到后續的模塊當中。故障1為電容C1發生故障,該故障的發生會導致前置的Sallen-Key電路的通帶發生變化,從而使得測試點1的輸出信號發生變化,測試點1在正常及故障狀態下的輸出結果如圖9所示。

圖9 測試點1在正常狀態及故障1下的輸出Fig.9 Output of measuring point 1 in normal status and fault 1

通過電路物理上的鏈接,測試點1的變化會影響測試點2、測試點4、測試點5,但是卻不能影響與之相連的測試點3。測試點2在正常狀態及故障1下的輸出結果如圖10所示。

圖10 測試點3在正常狀態及故障1下的輸出Fig.10 Output of measuring point 3 in normal status and fault 1

由圖10可以看出,測試點2信號在故障發生前后并沒有明顯的變化,這說明故障1的發生并沒有影響測試點2,故測試點2不應當被篩選為故障1變量影響集合,與本文方法結果一致。

3.3.2 關于因果關系的討論

從表6及表7的辨識結果來看,基于傳遞熵的辨識方法對于部分變量之間的因果方向指示并不正確。例如,如表6所示,當故障1發生時,測試點5的故障信號會傳播到測試點1、2、4;但從系統電路結構來看,測試點5位于整個電路系統的末端,且系統中并不存在跨模塊的反饋回路,因此當測試點5檢測出故障時,該故障并不能影響到測試點1、2、4的信號,這與表6的辨識結果不符。

表4及表5的辨識結果符合對于系統結構的分析,說明本文方法能夠指示測試點之間的因果關系。

從圖7可以看出,基于傳遞熵方法的時間成本遠遠超過本文所提方法的時間成本,這是因為傳遞熵在使用過程中涉及到概率密度估算問題,該問題目前并沒有較好的解決方案。而基于因果關系的方法,從根本上避免了概率密度估算的問題,因而時間成本被顯著壓縮。

4 結 論

當系統存在多個報警時,為了迅速定位根源故障并進行報警的清除,本文從因果關系的角度出發,揭示了故障發生的因果內涵,提出了一種基于因果關系的故障傳播路徑辨識方法。該方法通過數據之間的因果關系,首先進行故障相關變量的篩選,然后根據因果關系指標對變量之間的因果關系進行推斷,最后通過?？蛇_性的賦權有向圖最小生成樹算法,繪制因果圖,展現了故障在系統中的傳播路徑。雙帶通濾波器電路仿真實驗表明,本文算法對于故障傳播路徑的辨識具有一定的有效性和實用性。