修正三角模糊事故树算法在化工生产事故预测中的应用

张园园，张巨伟，赵爱彬，姚彦桃

（辽宁石油化工大学环境与安全工程学院，辽宁抚顺113001）

化工生产过程连续性较强，生产设备存在超温超压导致事故的风险，且一旦发生事故则可能产生极大的经济损失和人员伤亡，所以预测化工生产事故对安全管理和生产具有非常重要的指导意义。事故树分析法在我国化工事故预测中得到了广泛的应用，但是采用事故树分析法对实际问题进行分析时，事故树基本事件发生的概率往往难以确定，需要借助专家的力量进行抉择。因此，H.Tanaka等[1]首次把模糊数学引入到事故树分析中；G.S.Liang 等[2]在H.Tanaka 提出的理论基础上进行了进一步的分析。P.Gmytrasiewicz 等[3]和D.Singer[4]把模糊集合的理论引入到事故树分析中；T.Hideo 等[5]最早把三角模糊数学及梯形模糊数学应用于计算事故树顶上事件发生的概率；D.Singer[4]基于模糊截集的理论计算顶上事件发生的概率；P.V.Suresh等[6]基于不同的基本事件对顶上事件的影响，提出了一种新的进行基本事件模糊重要度排序的方法；J.Dunyak 等[7]提出了一种新的模糊事故树可靠性计算模型。在此基础上，模糊事故树理论在化工领域得到了广泛的应用。Y.H.Dong 等[8]应用模糊事故树理论分析了石化系统介质传输失效的可能性；王楷[9]采用模糊事故树分析氯气泄漏事故时引入维度的概念；S.M.Lavasani 等[10]应用模糊事故树评估化工厂中脱乙烷塔的失效概率，弥补了化工生产过程中由于现场数据不足导致管理措施不完善的不足；张园园等[11]把三角模糊数引入到化工生产事故分析中；赵莹莹等[12]把三角模糊数和事故树结合分析压力管道爆炸事故；蒲伟等[13]用三角模糊事故树理论评估石油泄漏的风险。上述研究结合模糊数学的理论和事故树进行事故分析或进行可靠性评估，将三角模糊数或梯形模糊数或抛物线型模糊数与专家的意见相结合，把事故树基本事件发生的概率模糊化，顶上事件发生的概率采用模糊算子进行基本运算。但是，这种方法在运算过程中存在不足之处，大部分对模糊数相乘的结果进行了线性的近似，而在实际分析过程中发现这种近似存在较大的误差。所以，本文针对这一点进行了修正分析，以期获得更加精确的模糊计算结果，降低事故概率预测的误差，减少化工事故的损失。

1 传统的三角模糊事故树预测模型

1.1 三角模糊数函数的定义

（1）模糊集的定义。若A 表示X 到[0，1]上的映射，X 为论域，用数学符号表示为A:X →[0,1],x ↦A，则称A 为模糊集，A(x)为隶属度函数。

（2）L-R 型模糊函数。A∈R，仅满足式（1）定义的函数，称为L-R 型函数。

式中，L(x) 为单调递增函数，0 ≤L(x) ＜1 且为单调递减函数，0 ≤R(x)＜1，且为函数左半部分隶属度1 所对应的概率；n 为函数右半部分隶属度1 所对应的概率。

（3）三角模糊函数。三角模糊函数是L-R 型函数的一种，三角模糊数函数曲线示意图如图1 所示。图1 中，m 表示A 核；l＋u 表示A 的盲度，记 为A=(m－l，m，m＋u)；Z 为x 轴上某一点；假设三角模糊数的隶属度函数曲线与x 轴围成图形的面积用S 表示，则S1、S2分别为被过Z 点的垂直于x 轴的线分解为两部分的图形的面积。

1.2 模糊事故树计算步骤

（1）构建事故树，确定最小割集或最小径集；

（2）用三角模糊数表示基本事件发生的概率。基本事件发生的概率可以根据三角模糊数的区间运算法则得到。若基本事件x1和x2发生的概率q1和q2的三角模糊概率表示为(m1-l1,m1,m1+u1)和(m2-l2,m2,m2+u2)，则三角模糊数的运算规则为：

（3）确定顶上事件发生的三角模糊概率。顶上事件的模糊发生概率可以通过模糊逻辑“或门”和模糊逻辑“与门”进行处理。模糊逻辑“与门”和模糊逻辑“或门”的计算公式为：

（4）运用中值法计算各个基本事件模糊重要度。图1 中，三角模糊函数的隶属度函数曲线与x轴围成图形的面积用S 表示，若x 轴上有一点Z，过该点作一条垂直于x 轴的垂线，把整个封闭面积S分成S1和S2两部分，若满足S1＝S2，则点Z 称为该三角模糊数的中位数。利用三角模糊数的中位数进行模糊重要度比较的规则如下：两个三角模糊数A1和A2的模糊中位数分别为Z1和Z2，如果Z1＜Z2，则A1＜A2；如果Z1＞Z2，则A1＞A2；如 果Z1＝Z2，则A1＝A2。对于有n 个基本事件组成的事故树，其结构函数为φ( x1,x2,…,xn)，三角模糊数的模糊中位数为Z，去掉第i 个基本事件的事故树的结构函数为φ( x1,x2,…,xi-1,0,xi+1,…,xn)，模糊中位数 为Zi，则第i 个基本事件的模糊重要度的表达式为：

2 传统预测模型误差分析

2.1 理论依据

给定论域U，对于任何模糊集合A，若λ ∈[0,1]，则称集合Aλ={x|A(x) ≥λ}为A 的水平截集。 若f:X = Xi→Y,Ai(i= 1,2,…,n)为 模糊集，那么，∀λ∈[0,1]，f (A1,A2,…,An)λ=f ((A1)λ,(A2)λ,…,(An)λ)，当 且 仅 当 ∀y ∈Y，使得f (A1,A2,…,An)(y)=

2.2 误差分析

在传统的三角模糊预测模型中，模糊事故树计算步骤（3）认为n 个三角模糊数相乘的结果是一个新的三角模糊数，该三角模糊数的隶属度函数曲线是由三角模糊数的左端点、中间点、右端点与其相应隶属度值组成的坐标用直线连接而成的，即用直线连接（三角模糊数左端点，0）、（三角模糊数右端点，0）、（三角模糊数中间点，1）3 点围成的三角形（见图1）。在模糊事故树计算步骤（4）中，根据截集定理可知，可以把事故树基本事件用模糊数表示，然后对每个模糊数取λ 水平截集，根据事故树的运算法则以及λ 水平截集的运算法则进行运算，得到用λ 水平截集表示的顶上事件发生的概率。n 个三角模糊数的截集相乘的结果是λ 的n 次方，当n＞1时，两端点与顶点之间的连线是曲线，只有在n=1时，两端点与顶点用直线连接才合适。所以，传统三角模糊事故树顶上事件模糊隶属度左端点与顶点直接用直线连接是存在误差的。因此，本文依据三角模糊数学的相关定理以及MATLAB 软件，对传统的三角模糊算法加以修正，绘制改进后三角模糊数相乘的隶属度函数图，确定每个基本事件较精确模糊重要度值。

3 修正的三角模糊事故树算法

3.1 理论依据

根据模糊数的表现定理，已知A=(ml,m,m+u)，则有：

3.2 修正策略

（1）若事故树有n 个基本事件，第i 个基本事件qi的概率用三角模糊数表示为Ai=(mili,mi,mi+ui)，且在(mi-li)、(mi+ui)点的隶属度为0，即Ai(mi-li)=0，Ai(mi+ui)=0，Ai在mi点的隶属度为1，即Ai(mi)=1。n 个三角模糊数的乘积依然是三角模糊数，在点和点的隶属度为0，在点的隶属度为1。在之间隶属度的求解策略是，对n 个三角模糊数取λ截集，第i 个三角模糊数的λ 截集区间为：[mili(1-λ),mi+ui(1-λ)]；区间左端点的乘积为：区间右端点的乘积为：

（3）隶属度用AR(xR)表示，在（0，1）等间距取若干值，求解所有的Rλ；由于在的λ 截集为单调递减函数，根据公式（9）可知，给定任意一个xR，大于等于xR的最小Rλ对应的λ 值即为对应的隶属度值AR(xR)；找到尽可能多的(xR,AR(xR))值，连接这些点以及点得到的曲线即为的隶属度曲线。

（4）在第（3）步绘制精确隶属度曲线的基础上，在图形中找到模糊中位数Z，根据式（7）即可求出各个基本事件的精确模糊重要度值。

3.3 修正策略MATLAB 程序

其中，mm1、a1、b1 为原始数据，是基本事件的三角模糊概率的中间值及左右区间值。

ma1=mm1－a1；mb1=mm1＋b1

la=prod(ma1)%，基本事件相乘三角模糊数左端点值；

ma=prod(mm1)%，基本事件相乘三角模糊数中间点值；

ra=prod(mb1)%，基本事件相乘三角模糊数右端点值；

q1=0.0001；

for k=1:9999

lan1=k*q1；

m1(k)=prod(mm1+b1.*(1-lan1))；

end

%在(0,1)等间距取若干λ 值，求解所有的Rλ，放在数组m1 中。

xR1=[]；b01=[]；c01=[]；result=[]

for i=ma:0.001:ra

result1=m1(find(m1＞=i))

b1=min(result1)

D1=[b01,b1]

b01=D1

xR1=[xR1,i]

end

%求数组中大于等于xR的最小Rλ值。

ARxR1=[]；n1=length(b01)

for j=1:n1

id1=find(m1=b01(j))

E1=[ARxR1,id1]

ARxR1=E1

end

ARxR1=ARxR1*q1

%确定xR处的隶属度值AR(xR)。

4 修正前后预测结果对比分析

4.1 油库静电积累三角模糊事故树分析

静电积累是导致油库发生火灾爆炸的重要原因，因此分析油库静电具有重要意义。采用三角模糊事故树分析法查找导致油库静电积累的原因事件，分析常见的原因事件有8 个：进油时油速过快（X1）、油罐内壁凹凸不平（X2）、流动的油液与金属容器碰撞（X3）、进油时人员操作失误（X4）、液体油滴与气体摩擦（X5）、油流中混有金属漂浮物（X6）、测量的工具不符合要求（X7）、油流静置的时间过短（X8）。8个基本事件发生的概率及不发生的概率用三角模糊数表示（见表1）。

表1 基本事件的三角模糊数

由表1 可知，每个基本事件都会导致顶上事件静电积累事故的发生，所以这些事件应该用“或门”连接。

采用修正算法，根据表1 中数据，采用MATLAB 软件，绘制所有事件同时发生时顶上事件三角模糊隶属度曲线，以及任何一个基本事件（Xi、i=1～8）不发生但其他基本事件同时发生时顶上事件的三角模糊隶属度曲线，得到每个事件的精确模糊重要度，再与通过传统三角模糊重要度计算方法得到的每个事件的近似模糊重要度进行对比并进行误差分析。

4.2 基于修正算法的基本事件精确模糊重要度的求解

4.2.1 所有事件同时发生 所有事件同时发生时的原始数据如下。

全部基本事件的三角模糊数中间值数据矩阵：

全部基本事件的三角模糊数左半部分区间值数据矩阵：

全部基本事件的三角模糊数右半部分区间值数据矩阵：

利用MATLAB 软件，通过修正算法，绘制精确隶属度图形，结果如图2 所示，以图2 为基础确定模糊中位数Z，依据模糊重要度的计算原理，点Z 把图2 中由三角模糊数的隶属度曲线与x 轴围成的图形分 为S1和S2两部分，S1=S2=0.105 9，S=S1＋S2=0.211 8，Z=0.838 5。

4.2.2 X1不发生 X1不发生时的原始数据如下。

除X1的其他基本事件的三角模糊数中间值数据矩阵：

除X1的其他基本事件的三角模糊数左半部分区间值数据矩阵：

除X1的其他基本事件的三角模糊数右半部分区间值数据矩阵：

利用MATLAB 软件，通过修正算法，绘制精确隶属度图形，结果如图3 所示，以图3 为基础确定模糊中位数Z1。由模糊中位数定义可知，Z1把图3 中由三角模糊数的隶属度曲线与x 轴围成的图形分为中分 成S1和S2两部分，S1=S2=0.123 4，S=S1＋S2=0.246 8，Z1=0.784 6。由模糊重要度计算原理可知，X1的模糊重要度为：M1=Z－Z1=0.838 5－0.784 6＝0.053 9。

4.2.3 X2不发生 X2不发生时的原始数据如下。

除X2的其他基本事件的三角模糊数中间值数据矩阵：

除X2的其他基本事件的三角模糊数左半部分区间值数据矩阵：

除X2的其他基本事件的三角模糊数右半部分区间值数据矩阵：

利用MATLAB 软件，通过修正算法，绘制精确隶属度图形，结果如图4 所示，以图4 为基础确定模糊中位数Z2。由模糊中位数定义可知，Z2把图4 中由三角模糊数的隶属度曲线与x 轴围成的图形分为分 成S1和S2两部分，S1=S2=0.110 2，S=S1＋S2=0.220 4，Z2=0.813 1。由模糊重要度计算原理可知，X2的模糊重要度为：M2=Z－Z2=0.838 5－0.813 1=0.025 4。

4.2.4 X3不发生 X3不发生时的原始数据如下。

除X3的其他基本事件的三角模糊数中间值数据矩阵：

除X3的其他基本事件的三角模糊数左半部分区间值数据矩阵：

除X3的其他基本事件的三角模糊数右半部分区间值数据矩阵：

利用MATLAB 软件，通过修正算法，绘制精确隶属度图形，结果如图5 所示，以图5 为基础确定模糊中位数Z3。由模糊中位数定义可知，Z3把图5 中由三角模糊数的隶属度曲线与x 轴围成的图形分为S1和S2两部分，S1=S2=0.111 6，S=S1＋S2=0.223 2，Z3=0.809 9。由模糊重要度计算原理可知，X3的模糊重要度为：M3=Z－Z3=0.838 5－0.80 99=0.028 6。

4.2.5 X4不发生 X4不发生时的原始数据如下。

除X4的其他基本事件的三角模糊数中间值数据矩阵：

除X4的其他基本事件的三角模糊数左半部分区间值矩阵：

除X4的其他基本事件的三角模糊数右半部分区间值矩阵：

利用MATLAB 软件，通过修正算法，绘制精确隶属度图形，结果如图6 所示，以图6 为基础确定模糊中位数Z4。由模糊中位数定义可知，Z4把图6 中由三角模糊数的隶属度曲线与x 轴围成的图形分为S1和S2两部分，S1=S2=0.127 7，S=S1＋S2=0.255 4，Z4=0.772 0。由模糊重要度计算原理可知，X4的模糊重要度为：M4=Z－Z4=0.838 5－0.772 0=0.066 5。

4.2.6 X5不发生 X5不发生时的原始数据如下。

除X5的其他基本事件的三角模糊数中间值数据矩阵：

除X5的其他基本事件的三角模糊数左半部分区间值矩阵：

除X5的其他基本事件的三角模糊数右半部分区间值矩阵：

利用MATLAB 软件，通过修正算法，绘制精确隶属度图形，结果如图7 所示，以图7 为基础确定模糊中位数Z5。由模糊中位数定义可知，Z5把图7 中由三角模糊数的隶属度曲线与x 轴围成的图形分为S1和S2两部分，S1=S2=0.119 6，S=S1＋S2=0.239 2，Z5=0.791 8。由模糊重要度计算原理可知，X5的模糊重要度为：M5=Z－Z5=0.838 5－0.791 8=0.046 7。

4.2.7 X6不发生 X6不发生时的原始数据如下。

除X6的其他基本事件的三角模糊数左半部分区间值矩阵：

除X6的其他基本事件的三角模糊数右半部分区间值矩阵：

利用MATLAB 软件，通过修正算法，绘制精确隶属度图形，结果如图8 所示，以图8 基础确定模糊中位数Z6。由模糊中位数定义可知，Z6把图8 中由三角模糊数的隶属度曲线与x 轴围成的图形分为S1和S2两部分，S1=S2=0.110 9，S=S1＋S2=0.221 8，Z6=0.810 0。由模糊重要度计算原理可知，X6的模糊重要度为：M6=Z－Z6=0.838 5－0.810 0=0.028 5。

4.2.8 X7不发生 X7不发生时的原始数据如下。

除X7的其他基本事件的三角模糊数中间值数据矩阵：

除X7的其他基本事件的三角模糊数左半部分区间值矩阵：

除X7的其他基本事件的三角模糊数右半部分区间值矩阵：

利用MATLAB 软件，通过修正算法，绘制精确隶属度图形，结果如图9 所示，以图9 为基础确定模糊中位数Z7。由模糊中位数定义可知，Z7把图9 中由三角模糊数的隶属度曲线与x 轴围成的图形分为S1和S2两部分，S1=S2=0.106 9，S=S1＋S2=0.213 8，Z7=0.818 0。由模糊重要度计算原理可知，X7的模糊重要度为：M7=Z－Z7=0.838 5－0.818 0=0.020 5。

4.2.9 X8不发生 X8不发生时的原始数据如下。

除X8的其他基本事件的三角模糊数中间值数据矩阵：

除X8的其他基本事件的三角模糊数左半部分区间值矩阵：

除X8的其他基本事件的三角模糊数右半部分区间值矩阵：

利用MATLAB 软件，通过修正算法，绘制精确隶属度图形，如图10 所示，以图10 为基础确定模糊中位数Z8。由模糊中位数定义可知，Z8把图10 中由三角模糊数的隶属度曲线与x 轴围成的图形分为S1和S2两部分，S1=S2=0.126 1，S=S1＋S2=0.252 2，Z8=0.774 9。由模糊重要度计算原理可知，X8的模糊重要度为：M8=Z－Z8=0.838 5－0.774 9=0.063 6。

4.3 修正前后的误差分析

修正前的模糊中位数分别为：Z=0.882 5，Z1=0.835 6，Z2=0.861 1，Z3=0.857 9，Z4=0.823 0，Z5=0.840 8，Z6=0.857 0，Z7=0.864 0，Z8=0.825 9。修正前X1、X2、X3、X4、X5、X6、X7、X8的模糊重要度分别为：M1=0.046 9，M2=0.021 4，M3=0.024 6，M4=0.059 5，M5=0.041 7，M6=0.025 5，M7=0.018 5，M8=0.056 6。

修正后的模糊中位数分别为：Z=0.838 5，Z1=0.784 6，Z2=0.813 1，Z3=0.809 9，Z4=0.772 0，Z5=0.791 8，Z6=0.810 0，Z7=0.818 0，Z8=0.774 9。修正 后X1、X2、X3、X4、X5、X6、X7、X8的模糊重要度为：M1=0.053 9，M2=0.025 4，M3=0.028 6，M4=0.066 5，M5=0.046 7，M6=0.028 5，M7=0.020 5，M8=0.063 6。

通过修正三角模糊事故树预测方法得知，X1、X2、X3、X4、X5、X6、X7、X8的模糊重要度误差分别减小了 12.99%、15.75%、13.99%、10.53%、10.71%、10.53%、9.76%、11.00%，误差平均减小11.91%。

5 结 论

（1）基于截集定理进行分析得知，使用三角模糊事故进行预测，当基本事件的个数为2 个或2 个以上时，顶上事件三角模糊隶属度曲线是非线性的，所以通过传统三角模糊事故树进行线性预测存在较大误差。

（2）基于表现定理得到三角模糊事故树顶上事件非线性隶属度曲线的算法，以及基本事件模糊重要度的较精确计算方法。

（3）采用修正三角模糊事故树预测算法，以化工生产过程中因静电积累导致油库火灾爆炸事故为例，计算了8 个用“或门”连接的基本事件的精确模糊重要度。与通过传统三角模糊事故树预测基本事件的模糊重要度进行了对比，结果表明修正算法的模糊重要度误差明显减小，平均减小11.91%。