基于矩陣的故障樹分析方法

郭永晉，孫麗萍

（哈爾濱工程大學船舶工程學院，黑龍江哈爾濱150001）

基于矩陣的故障樹分析方法

郭永晉，孫麗萍

（哈爾濱工程大學船舶工程學院，黑龍江哈爾濱150001）

為開發大型通用故障樹分析程序、優化程序算法、降低NP困難問題，將矩陣引入到故障樹分析過程中，基于矩陣對故障樹進行結構編碼和參數轉化。闡述了應用矩陣求解故障樹最小割集、最小路集、不交化最小割集、頂事件發生概率和底事件重要度的方法和步驟。使用MATLAB軟件開發了相應的故障樹分析程序，并將其應用于風機齒輪箱失效的研究中，程序運行速度快、計算結果準確，表明基于矩陣的故障樹分析方法是有效可行的。

矩陣；故障樹分析法；定性分析；定量分析；MATLAB；NP困難問題

故障樹分析法 （fault tree analysis，FTA）是安全系統工程的重要方法之一，由美國貝爾實驗室在20世紀60年代提出，最初用來預測民兵導彈發射的失效概率。而后，FTA得到了極大的發展，如今已經廣泛應用于航天工業、機械制造、核電站、海洋工程等領域的系統可靠性分析中［1］。

在進行故障樹分析時，對故障樹的建模和計算是非常繁瑣費時的，當分析對象趨于大型化、復雜化，必須借助故障樹分析軟件來實現。在故障樹分析計算機化的過程中主要的難點是對于NP（nondeterministic polynomial）困難問題的處理，也就是FTA的計算量隨著故障樹規模的加大而呈指數增長這一問題。對于一些大型故障樹，其基本事件可能達到幾千個，其中還存在著大量的重復事件，如果不能解決組合爆炸問題，故障樹分析將難以進行［2］。針對這一難點，國內外研究者作了大量的工作［3］，其中一些學者在編程過程中應用矩陣來解決NP困難問題，但處理方法相對簡單，沒有完整介紹基于矩陣進行故障樹分析的全過程［4-6］。本文采用矩陣的方法對故障樹分析過程進行優化，闡述了優化算法及步驟。在此基礎上，采用MATLAB對故障樹分析算法進行編程，實現了故障樹最小割集、最小路集、頂事件發生概率及底事件重要度的求解。

1　故障樹編程預處理

1.1故障樹的規范化

故障樹由邏輯門與事件構成，常用的邏輯門是邏輯“與門”和邏輯“或門”，其他邏輯門在某種程度上都可以簡化為這兩種門類型［2，7］。在本文中故障樹經過規范化處理只含有“與門”、“或門”、底事件、中間事件和頂事件。圖1為一個經過規范化的故障樹，該故障樹中共有6個底事件，4個中間事件和5個邏輯門。

圖1　故障樹示意圖Fig.1 A fault tree

1.2故障樹結構編碼

為進行程序開發，需要采用一定的規則對故障樹進行數學描述。仔細考察故障樹可以發現，故障樹用門的輸入和輸出關系來描述基本事件和頂事件之間的關系，邏輯門在故障樹中起連接作用:對上，它代表一級原因事件；對下，它又代表一級結果事件。抓住了故障樹邏輯門的輸入輸出關系，整個故障樹也就清晰了。

本文采用矩陣對故障樹進行結構編碼，編碼規則如下:

1）底事件序號從1～2000依次編號，即故障樹底事件最大可能數目為2 000個。

2）頂事件和中間事件序號從2001開始編號。

3）用數字1表示邏輯“與門”，數字2表示邏輯“或門”。

4）各中間事件輸入事件個數不相等時，在輸入事件后使用數字0占位。

根據以上規則進行處理后，可將故障樹以矩陣形式表示。如圖1中的故障樹可表示為矩陣A。其中第1列是頂事件和中間事件的序號，第2列是連接輸出事件邏輯門的類型，第3、4、5列是相應邏輯門輸入事件的序號:

2　故障樹的定性分析

對故障樹進行定性分析的主要目的是找出它的所有最小割集和最小路集。

2.1全體割集的求法

用1×n（底事件總個數）的矩陣代表底事件，其中底事件序號位置元素1，其余位置元素為0。例如共有6個底事件，則序號為3的底事件可轉化為矩陣［0 0 1 0 0 0］。

用列數為n的矩陣代表故障樹的頂事件和中間事件，它們由其輸入事件矩陣間的運算生成。定義“與”運算為矩陣的相加，即將多個矩陣每行分別依次相加生成新矩陣，則新矩陣的行數為所有矩陣行數的積，矩陣經過“與”運算后如果出現大于1的元素，將此元素賦值為1。定義“或”運算為矩陣的合并，即將多個矩陣豎向合并生成新矩陣，則新矩陣的行數為所有矩陣行數的和。按此規則對故障樹進行處理可得到割集矩陣B，矩陣中每一行代表一個割集，數字為1元素的列數代表底事件序號。

圖1故障樹的割集矩陣為式（2），共有4行，代表該故障樹共有4個割集，分別為｛1，2｝、｛1，4，5，6｝、｛2，3｝、｛2，3，4，5，6｝。

2.2割集的最小化

找出故障樹的全體割集后需要將割集簡化、吸收得到全部最小割集。在故障樹定性分析程序中，通常對全部割集進行兩兩比較，假如割集數量為p個，需要比較p（p-1）/2次。在一些大型故障樹中，特別是與或門交替出現時，割集數量會變得非常大，此時需要比較的次數相對也會變得很大。針對這一特點，本文提出采用分組對比的方法，流程如下:

1）提取割集矩陣B的每一行，計算非零元素個數q（即每個割集包含的底事件的個數）。

2）對所有割集進行分組，將t相等的割集歸在同一組。

3）將所有分組按q的數值由小到大依次排列，將每組中的割集與比其q值大的組中的割集依次比較，刪除可以被吸收的割集。

4）將每組中沒有被吸收的割集提取，生成最小割集矩陣C。

由于此方法避免了對q值相同的割集進行比較，并且對分析過程中可被吸收的割集直接刪除，所以計算量可在一定程度上減小，當分組較少時，減小更為明顯。此方法在下面的故障樹定量計算過程中也有應用。

圖1故障樹的最小割集矩陣為式（3），可以得出此故障樹共包含有3個最小割集，分別為｛1，2｝、｛1，4，5，6｝、｛2，3｝。

2.3最小路集的求法

最小路集就是原故障樹對偶樹的最小割集，只需將原故障樹的與或門類型進行互換并求解即可。用此方法可以求得圖1故障樹的最小路集分別為｛1，2｝、｛1，3｝、｛2，4｝、｛2，5｝、｛2，6｝。

3　故障樹的定量分析

故障樹定量分析的目的是求出頂事件發生概率及底事件重要度。一般情況下，最小割集彼此相交，根據相容事件的概率計算公式求解會產生組合爆炸問題，所以必須對最小割集進行不交化處理，把最小割集的相交和通過不交化方法變成不交和，再求頂事件發生概率［8］。

MATLAB的數據對象為數組，是進行運算的基本單元，其提供的數組類型主要有數值數組、字符串數組、元胞數組等。其中元胞數組是MATLAB的一種特殊數據類型，允許在一個數組里存放各種不同類型的數據。本文編寫的程序大量應用了數組。

3.1最小割集的不交化

應用遞歸法進行不交化運算，設故障樹有k個最小割集Kj（1≤j≤k），則該故障樹結構函數的不交型表達式為

應用矩陣進行最小割集不交化的過程主要依靠不交型積之和定理［2］，為便于描述不交化實現方法，首先定義兩個運算:

1）X運算:1×n割集矩陣C1與1×n割集矩陣C2進行比較，如果兩矩陣中相同位置上的元素都為1，則將矩陣C1此位置元素更改為0。現以割集｛1，2｝與｛1，4，5，6｝為例說明X運算:

2）Y運算:首先計算1×n矩陣D中數值為1的元素個數N，將矩陣D豎向羅列N次生成矩陣E，保留矩陣E中第i（1≤i≤N）行前i個元素1，將其他所有元素賦值為0。最后，將矩陣中每行最后一個數值為1的元素賦值為-1。例如對矩陣D1進行Y運算的計算過程如下

編寫定量分析程序時需要分別提取不交型表達式的每一項，生成每一項的不交化矩陣。對于不交項，首先將代表割集Kj（1≤j≤k-1）的矩陣分別存入元胞數組G，將G中每個元素分別與Kk進行X運算，得到Kbj，將此結果更新到元胞數組G中。隨后采用計算最小割集使用的分組對比方法，對G中的元素進行吸收化簡，將可以被吸收的Kbj刪除。最后分別提取G中每個元素進行Y運算，將結果存入元胞數組H中，求解不交化矩陣Jk，求解流程如圖2所示。

圖2　不交化矩陣求解流程圖Fig.2 Solution flow chart of non-cross matrix

將所有項的不交化矩陣豎向合并可以生成矩陣J，該矩陣即為整個故障樹的不交化矩陣。下式為圖1故障樹的不交化矩陣:

3.2頂事件發生概率計算

在不交化矩陣J中，若各底事件的發生概率為fi（1≤i≤n），則矩陣中數字為1表示該位置發生概率為fi；數字-1表示該位置發生概率為 （1-fi）；數字0表示該位置不需計算。將對應概率值賦到矩陣J中，可生成不交化計算矩陣M。將M每行的概率值相乘，再將所有結果相加即是頂事件發生概率。

若圖1所示故障樹各底事件的發生概率分別為f1=f2=0.1，f3=f4=0.2，f5=f6=0.3，那么其頂事件的發生概率為

應用容斥定理概率計算公式可以求得頂事件發生概率的解析值為

可見兩種方法的計算結果相同，由此驗證了方法和程序的準確性。

3.3底事件重要度計算

故障樹中，用底事件重要度描述其對系統故障的影響大小，主要包括概率重要度和關鍵重要度的計算。若頂事件發生概率為F，則第i個底事件的概率重要度IPR和關鍵重要度ICR分別定義為

概率重要度就是故障樹失效概率計算公式分別對每個底事件概率fi進行求導得來的。所以可以得出兩個結論:1）若計算式中某一項沒有某個底事件概率的影響，則對該底事件概率求導過后該項為0。2）若計算式中某一項某個底事件概率的影響是底事件對立事件的發生概率，即對應概率為（1-fi），則求導后該底事件概率對該項的影響是使該項添加了一個負號。基于這兩個結論對不交化計算矩陣M做相應的變化即可求得底事件重要度。求第i個底事件的概率重要度實現過程為:

1）將矩陣M第i列刪去，得到矩陣Mi，將Mi每行概率值分別相乘，將得到的結果依次存入到1×n矩陣N中。

2）遍歷矩陣J中第i列元素，若值為0，則將矩陣N中該行概率相乘值賦為0；若值為-1，則將矩陣N中該行概率相乘值賦值為其相反數；若值為1，則不作改變。

3）矩陣N中所有元素的和即為該底事件的概率重要度。

底事件關鍵重要度將概率重要度乘以相應的概率比即可求得。經過計算圖1中故障樹6個底事件的概率重要度分別為0.096 2，0.278 2，0.090 0，0.008 1，0.005 4，0.005 4，關鍵重要度分別為0.324 8，0.939 2，0.607 7，0.054 7，0.054 7，0.054 7。

4　應用實例

目前國內安裝使用的風力發電機組大多為非直驅式機組。其中，齒輪箱是整套機組中非常重要的中間部件，其故障直接威脅風電機組的安全可靠運行，所以分析研究齒輪箱故障原因及預防措施具有重要意義［9］。圖3所示為以“齒輪箱失效”為頂事件的故障樹，該故障樹共有13個底事件，表1為所有事件的描述，各底事件發生概率如表2所示。

圖3　風機齒輪箱故障樹Fig.3 The fault tree of wind turbine gearbox

將故障樹初始數據輸入程序進行計算，結果顯示該故障樹共有9個最小割集，分別為｛1，2｝，｛1，3｝，｛1，4｝，｛1，5｝，｛1，6｝，｛1，7｝，｛1，13｝，｛1，11，12｝，｛1，8，9，10｝。頂事件發生概率為3.002 1× 10-3。各底事件重要度如表3所示。分析計算結果可知定期對齒輪箱進行維護非常重要，從底事件重要度的角度考慮，如果齒輪箱工作異常，可考慮首先檢查底事件2、4、5、6、7、13，然后再對齒輪箱其他故障進行檢測。

表1　事件描述Table 1 Description of events

表2　底事件發生概率Table 2 Probability of failure of basic events

表3　底事件重要度Table 3 Importance of basic events

5　結論

故障樹分析法廣泛應用于可靠性工程中，本文將矩陣引入到故障樹分析中，基于MATLAB軟件開發出了相應的分析程序，得到如下結論:1）使用矩陣可以方便地對故障樹進行結構編碼，表達出故障樹的邏輯關系。結構矩陣的產生可以為進一步的故障樹分析提供初始參數。2）在故障樹分析過程中應用分組對比的方法和不交型積之和定理可以有效地降低NP困難問題，基于矩陣的變化和矩陣間的運算可以開發出相應分析過程的優化算法，且易于計算機編程實現。3）由于矩陣和優化算法的應用，本文編寫的故障樹分析程序具有運算速度快，計算準確的優點。基于矩陣的故障樹分析方法為計算機輔助故障樹分析提供了一個新的思路。

［1］曾聲奎，趙廷弟，張建國，等.系統可靠性設計分析教程［M］.北京:北京航空航天大學出版社，2001:117-144. ZENG Shengkui，ZHAO Tingdi，ZHANG Jianguo，et al. System reliability design and analysis course［M］.Beijing: Beihang University Press，2001:117-144.

［2］金星，洪延姬，沈懷榮，等.工程系統可靠性數值分析方法［M］.北京:國防工業出版社，2002:142-215. JIN Xing，HONG Yanji，SHEN Huairong，et al.Numerical analysis methods of reliability for engineering systems［M］. Beijing:National Defense Industry Press，2002:142-215.

［3］RUIJTERS E，STOELINGA M.Fault tree analysis:a survey of the state-of-the-art in modeling，analysis and tools［J］. Computer science review，2014，15-16:29-62.

［4］胡云昌，陳金水.求系統失效樹最小割集的新方法［J］.中國造船，1989，104（1）:112-122. HU Yunchang，CHEN Jinshui.New method for searching the minimum cutsets of system fault tree［J］.Shipbuilding of China，1989，104（1）:112-122.

［5］陳金水，蔡佳昆，蔡惠明.利用矩陣運算實現故障樹參數間轉化［J］.天津大學學報，2004，37（8）:682-685. CHEN Jinshui，CAI Jiakun，CAI Huiming.Transformation of main parameters by matrixing method in FTA［J］.Journal of Tianjin university，2004，37（8）:682-685.

［6］曹利鋒，鄒樹梁，唐德文，等.基于VC++與MATLAB的故障樹分析系統［J］.計算機技術與發展，2014，24（1）: 77-80，84. CAO Lifeng，ZOU Shuliang，TANG Dewen，et al.Fault tree analysis system based on VC++and MATLAB［J］.Computer technology and development，2014，24（1）:77-80，84.

［7］STAPELBERG R F.Handbook of reliability，availability，maintainability and safety in engineering design［M］.London:Springer，2009:546-552.

［8］羅航，王厚軍，黃建國，等.故障樹定量分析及其交互方式的實現［J］.電子測量與儀器學報，2010，24（5）:473-480. LUO Hang，WANG Houjun，HUANG Jianguo，et al.Quantitative fault tree analysis and its realization of interaction mode［J］.Journal of electronic measurement and instrument，2010，24（5）:473-480.

［9］杜亮.風電機齒輪箱常見故障分析、診斷與預防［J］.內蒙古電力技術，2011，29（2）:15-18. DU Liang.Analysis，diagnose and prevention to common problems of wind mill generator gearbox［J］.Inner Mongolia electric power，2011，29（2）:15-18.

本文引用格式：

郭永晉，孫麗萍.基于矩陣的故障樹分析方法［J］.哈爾濱工程大學學報，2016，37（7）:896-900.

GUO Yongjin，SUN Liping.Fault tree analysis method based on a matrix［J］.Journal of Harbin Engineering University，2016，37（7）:896-900.

Fault tree analysis method based on a matrix

GUO Yongjin，SUN Liping
（College of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，China）

In order to develop a large general program for fault tree analysis（FTA），optimize the algorithm used，and reduce the NP-hard problem，we introduce a matrix into the FTA process，and conduct structural encoding and parameter conversion for the fault trees，based on this matrix.We illustrate the methods and procedure for using the matrix in calculating the minimal cut sets，minimal path sets，non-cross minimal cut sets，probability of occurrence of a top event，and the importance of bottom events.We develop a corresponding FTA program utilizing MATLAB software.In addition，we apply this program to the failure research for the gearbox of wind turbines.The program operating speed is fast and the calculation results are exact，which prove the effectiveness and feasibility of the matrix-based FTA methods.

matrix；fault tree analysis；qualitative analysis；quantitative analysis；MATLAB；NP-hard problem

10.11990/jheu.201505052

TB114.3

A

1006-7043（2016）07-896-05

2015-05-18.網絡出版日期：2016-06-13.

國家國際科技合作專項項目（2013DFE73060）.

郭永晉（1990-），男，碩士；孫麗萍（1962-），女，教授，博士生導師，博士.

孫麗萍，E-mail:sunliping@hrbeu.edu.cn.

網絡出版地址：http：//www.cnki.net/kcms/detail/23.1390.u.20160613.1436.002.html