郭永晋，孙丽萍

（哈尔滨工程大学船舶工程学院，黑龙江哈尔滨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