1 前言

化工园区是目前石油和化学工业发展的新模式和重要载体，然而产业联合、集约化、一体化的化工园区发展模式带来生态效益和经济效益的同时，也增加了一定的安全风险^[1]：化工园区内危险源数量多、密度大；生产工艺复杂、工序繁多；周边环境复杂、影响面大。化工园区局部的事故能够沿着多条路径传播扩散，少数具有高危传播状态的路径甚至会引起整个化工园区灾难性事故的发生^[2-5]。面对庞大复杂的化工园区事故链网络，如何高效、有针对性地进行重点监测和预防，提升化工园区整体安全性这一问题亟待解决。

目前，国内外学者主要对化工园区事故风险评价^[6]、事故多米诺效应^[7]、事故类型后果评价及模拟^[8]等方面进行研究，而针对化工园区事故链模型建立并识别其重要节点和关键事故传播路径的相关研究则较少，KHAKZAD等^[9]在考虑多种事故场景下结合图论利用图中心性度量识别化工厂中最脆弱的装置；陈国华等^[10]提出基于多灾种耦合关系模型研究自然灾害诱发的化工园区事故，通过对比不同节点对系统的扰动值进行重要性排序，但除自然灾害外其他影响因素并未考虑其中。

在应用复杂网络对节点重要性分析方面，国内外学者做了大量的研究：王纪洋等^[11]构建了AHP-TOPSIS模型来评价危险化工工艺风险等级，该方法对各指标进行主观赋权，解决了指标的量纲问题和信息易丢失问题，提高了结果准确性；LÜ等^[12]提出的LeaderRank算法在不考虑节点的结构特性和节点间影响的情况下，排序结果唯一且对网络噪音有更好的容忍性；姜英等^[13]采用客观赋权的TOPSIS方法初步确定关键子系统，然后利用LeaderRank算法识别子系统中的关键节点，降低了建模和计算的复杂性；LI等^[14]提出的加权LeaderRank算法可以使网络中其他节点从Ground节点中获取到部分权值，该方法对错误数据有很强的抗干扰能力，但没有考虑到LR值均分这一问题。

在应用网络传播行为研究传播路径方面，LIU等^[15]提出基于小世界聚类的故障传播模型并利用Dijkstra算法对故障传播路径进行了分析，该方法理论完善且运算稳定性强，但其应用对多源有向加权的化工园区事故链网络具有局限性；白翠粉等^[16]通过蚁群算法求取电力变压器网络的最高风险故障传播路径并获得了与实际故障统计相一致的结果，该方法具有正反馈机制和较强的鲁棒性，但是其将所有节点的故障传播概率设为固定值忽略了不同节点影响的差异性。

针对以上问题，本文首先综合考虑人、物、环境和管理等因素对化工园区事故的影响，构建了化工园区事故链网络模型；然后基于主客观组合赋权的TOPSIS法确定无向网络的重要节点，利用引进链入链出因子的加权LeaderRank算法解决LR值均分问题并确定有向网络的重要节点，结合这2种方法全面评价化工园区事故链网络的节点重要性；最后提出基于事故传播强度的方法识别关键事故传播路径，为化工园区的安全生产和监测控制提供理论支持。

2 化工园区事故链的复杂网络模型建立 2.1 确定化工园区事故链网络模型结构

基于复杂网络理论^[13]及前人对事故链模型^{[5, 17-18]}的相关研究，按照致灾因子→承灾体→事故类型→事故后果的事故链路线，结合事故发生的因果关系及事故演化机理建立了一般性的化工园区事故链网络模型，依据化工事故数据^[19-20]以及实际调研进一步完善了该模型。其中，事故链路线的各组成部分可继续递阶层次划分，如致灾因子中人的因素包括误操作、疲劳作业等，事故类型中火灾包括池火、喷射火等，具体影响因素间的相关性体现在事故链网络模型结构图中，如图 1所示。

2.2 构建化工园区事故链的复杂网络模型

利用复杂网络理论将化工园区事故链中的各影响因素抽象成网络的节点，影响因素之间的关系抽象成网络的边，定义邻接矩阵见式(1)。

$ {\rm{A}} = \left[ {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}&{...}&{{a_{1n}}}\\ {{a_{21}}}&{{a_{22}}}&{...}&{{a_{2n}}}\\ {...}&{...}&{...}&{...}\\ {{a_{n1}}}&{{a_{n2}}}&{...}&{{a_{nn}}} \end{array}} \right] $

(1)

在构建有向网络模型中，当节点i指向节点j时，a_ij=1；否则a_ij=0。将邻接矩阵输入到Ucinet6.0软件中，可得到化工园区事故链网络拓扑结构。

3 化工园区事故链网络的节点重要性评价 3.1 基于改进的TOPSIS方法的节点重要性评价

为了更加全面、准确量化节点在整个网络的重要程度，本文选取了点度中心性(DC)、接近中心性(CC)、中间中心性(BC)和特征向量中心性(EC)等指标对节点进行综合评价^[21]。

3.1.1 权重的确定

综合考虑了影响因素原始数据间的关系及专家根据实际数据做出的判断，利用组合权重的方法来确定权重。

3.1.1.1 熵权法客观指标权重确定

熵权法是根据影响因子的变异程度，通过信息熵得到各影响因子的熵权，进而获取各影响因子权重。其计算步骤如下：

(1) 对各评价指标数值进行标准化处理

$ {p_{ij}} = \frac{{{x_{ij}}/{x_{j\max }}}}{{\sum\limits_{i = 1}^n {{x_{ij}}/{x_{j\max }}} }} $

(2)

当x_ij/x_jmax=0，p_ij=0时，公式无意义。因此需要对p_ij进行修正：

$ {p_{ij}} = \frac{{{x_{ij}}/{x_{j\max }} + {{10}^{ - 4}}}}{{\sum\limits_{i = 1}^n {\left( {{x_{ij}}/{x_{j\max }}} \right) + {{10}^{ - 4}}} }} $

(3)

(2) 确定各影响因子熵权及权重

$ {e_j} = - \frac{1}{{\ln n}}\sum\limits_{i = 1}^n {{p_{ij}}} \cdot \ln {p_{ij}} $

(4)

$ {w_j} = \frac{{1 - {e_j}}}{{\sum\limits_{j = 1}^m {1 - {e_j}} }} $

(5)

3.1.1.2 层次分析法(AHP)主观指标权重确定

利用AHP进行主观赋权计算时首先根据三标度法规则建立比较矩阵，再构建判断矩阵，最后进行一致性检验，若一致性检验指标C.I.≤0.10，则接受判断矩阵计算得到的权重数据^[21]。

3.1.1.3 组合权重的确定

引入拉格朗日函数和欧式距离函数使两个权重值与偏好系数保持一致性，将AHP和熵权法组合计算以得到更加合理的指标权重^[22]。其计算步骤如下：

(1) 构建优化决策模型

$ \left\{ \begin{array}{l} {W_j} = \alpha {w_{{\rm{A}}j}} + \beta {w_{{\rm{B}}j}}\\ \alpha + \beta = 1 \end{array} \right. $

(6)

(2) 建立权重关联方程

$ \left\{ \begin{array}{l} D\left( {{w_{{\rm{A}}j}}, {w_{{\rm{B}}j}}} \right) = \sqrt {\sum\limits_{j = 1}^m {{{\left( {{w_{{\rm{A}}j}} - {w_{{\rm{B}}j}}} \right)}^2}} } \\ D{\left( {{w_{{\rm{A}}j}}, {w_{{\rm{B}}j}}} \right)^2} = {\left( {\alpha - \beta } \right)^2} \end{array} \right. $

(7)

联立式(6)和式(7)即可得到综合权重W_j。

3.1.2 基于TOPSIS法计算节点重要性

TOPSIS法可以定性的描述节点在整个网络结构中的信息和重要程度。其具体计算步骤见文献[21]，经计算得到所有节点的综合评价指标C_i^*，C_i^*的值越大，说明节点的评价越好，节点越重要。

3.2 基于改进的加权LeaderRank算法的节点重要性评价

为解决LR值均分问题，本文引入链入链出因子对加权LeaderRank算法进行改进，主要思想是根据节点的重要程度分配LR值，进一步提高了排序的准确性。其计算步骤如下：

(1) 计算不同节点间边的权值

$ \omega _{j, i}^{{\rm{in}}} = \frac{{{I_i}}}{{\sum\limits_{p = B(j)} {{I_P}} }} $

(8)

$ \omega _{j, i}^{{\rm{out}}} = \frac{{{O_i}}}{{\sum\limits_{p = B(j)} {{O_P}} }} $

(9)

$ {\omega _{j, i}} = \mu \omega _{j, i}^{{\rm{out}}} + (1 - \mu )\omega _{j, i}^{{\rm{in}}} $

(10)

(2) 计算各节点的LR值

$ {\rm{L}}{{\rm{R}}_i}(t{\rm{ + 1}}) = \sum\limits_{j = 1}^{N + 1} {\frac{{{\omega _{ji}}}}{{\sum\limits_{l = 1}^{N + 1} {{\omega _{jl}}} }}} {\rm{L}}{{\rm{R}}_j}(t) + \frac{{{\omega _{gi}}}}{{\sum\limits_{l = 1}^{N + 1} {{\omega _{gl}}} }}{\rm{L}}{{\rm{R}}_g}(t) $

(11)

结合2种方法得到的C*值和LR值对节点重要性进行排序，确定整个化工园区事故链网络的重要节点，以提高结果的合理性及准确性。

4 基于事故传播强度的关键事故传播路径识别

本文提出基于节点度数的负载分配法来计算节点事故概率，将其与有向边事故影响度融合表征有向边的事故传播属性，结合有向边事故传播属性和模型结构特征边介数负荷属性，定义事故传播强度以寻找关键事故传播路径。

4.1 有向边事故传播属性评估 4.1.1 基于节点度数负载分配方法的节点事故概率

在现实网络中，过载传播路径事故的概率通常与该事故传播路径的负载容量及初始负载有着密切关系^[23]，故本文提出基于节点度数的负载分配方法来计算节点事故概率，其计算步骤如下：

(1) 确定节点i的初始负载N_i

$ {N_i} = {({k_i}(1 + (\sum\limits_{m \in {\tau _i}} {{k_m}} ) - {k_i}))^\delta } $

(12)

(2) 邻点所接受的负载分配ΔN_j

$ \Delta {N_j} = \frac{{{N_i}{{({k_j}(1 + (\sum\limits_{n \in {\lambda _j}} {{k_n}} ) - {k_j}))}^\delta }}}{{\sum\limits_{m \in {\lambda _j}} {({k_m}{{(1 + (\sum\limits_{h \in {\lambda _j}} {{k_h}} ) - {k_m})}^\delta })} }} $

(13)

(3) 节点事故概率P_i

负载分配后节点的负载变为N_i'=N_i+ΔN，网络中的节点在发生过载之后会有一定概率出现事故，即节点事故概率为：

$ {P_i} = \frac{{{N_i}^\prime - {N_i}}}{{{N_i}^\prime }} $

(14)

4.1.2 基于LR值的有向边事故影响度

若事故影响度值高的节点发生事故并向邻点传递，则两节点间有向边的事故影响度值也较高。本文选择3.2节中计算得到的LR值来表示节点的事故影响度，将有向边的事故影响度LR(e_ij)定义为：

$ {\rm{LR}}({e_{ij}}) = \sqrt {{\rm{L}}{{\rm{R}}_i} \cdot {{\rm{L}}{{\rm{R}}_j}}} $

(15)

4.1.3 有向边事故传播属性

节点的事故传播属性受节点自身事故概率及邻近节点的影响，而有向边事故影响度反映了节点发生事故时对相邻节点的影响度，因此将二者结合来表征有向边的事故传播属性Q(e_ij)，即：

$ Q({e_{ij}}) = {P_i} \cdot {\rm{LR}}({e_{ij}}) $

(16)

4.2 边介数负荷属性评估

边介数反映了事故在整个网络中的流经频率，边介数越大即负荷量越大，该有向边对其它节点的联络控制作用越强，容易导致事故的迅速传播。边介数L(e_ij)表示如下：

$ L({e_{ij}}) = \sum\limits_{s = 1}^N {{p_{{{({e_{ij}})}_s}}}} $

(17)

4.3 有向边事故传播强度评估

有向边事故传播强度是发生事故时通过某有向边进行传播的强值，事故传播强度越大，表示事故通过此支路的传播代价越小且后果越严重^[24]。依据化工园区事故链网络模型，定义有向边事故传播强度I(e_ij)为：

$ I({e_{ij}}) = \frac{{Q({e_{ij}})}}{{\sum {Q({e_{ij}})} }} + \frac{{L({e_{ij}})}}{{\sum {L({e_{ij}})} }} $

(18)

将事故影响度与事故概率乘积值最大的初始节点确定为事故源，然后根据有向边事故传播强度值来识别网络模型中的关键事故传播路径。对关键事故传播路径及其中涉及的节点进行重点保护和防范，以降低事故灾害的发生。

5 案例分析

选取广东省惠州市某一典型化工园区为研究对象，该化工园区目前主要依托大炼油、大乙烯项目，生产和发展炼油、芳烃下游、C2下游、C3下游、C4下游和炼化副产品等^[18]。通过资料收集和调研分析发现，该化工园区所在地的自然灾害主要有地震、飓风、暴雨、洪水等，化工园区企业的危险品储存种类复杂、数量繁多，生产工艺复杂，发生化工事故易引发连锁效应。

5.1 复杂网络模型构建

基于第2节相关理论以及该化工园区实际情况，得到如图 2所示的化工园区事故链复杂网络拓扑结构图，其变量定义如表 1所示。

5.2 节点重要性评价结果

利用式(2)~(7)计算该化工园区事故链网络的综合权重，即：DCw_ij为0.086 7，BCw_ij为0.276 6，CCw_ij为0.146 3，ECw_ij为0.490 4。将综合权重结果代入TOPSIS法，经计算可得到综合各指标数值的总体评价C*；利用改进的加权LeaderRank算法对该有向网络进行节点重要性评价得到LR值，结果如表 2所示。

结合2种算法结果得出该化工园区事故链网络的重要节点为3、8、22、25、26、35、36、37、57、60，即误操作易造成生产工艺条件失常，以及其他致灾因素导致的反应器失效、管道断裂、容器储罐破裂等进而引发火灾、爆炸、泄漏事故的发生，化工园区事故的传播和升级最终会造成人员伤亡、经济损失等严重后果。该结果与化工企业事故案例统计^[19-20]中高频率出现的事故原因、事故类型及事故后果等相符，验证了本文提出节点重要性评价方法的准确性。

5.3 基于事故传播强度识别关键事故传播路径

化工园区事故链网络模型中有向边的事故传播强度由各有向边的事故传播属性和边介数属性决定，根据式(12)~(18)计算得到该化工园区事故链复杂网络拓扑结构图中各有向边的事故传播强度如表 3所示。

通过比较初始节点事故影响度与节点事故概率乘积值确定节点3为事故源，根据表 3中各有向边事故传播强度值的大小得到该化工园区事故的关键传播路径为：3→8→37→36→42→35→41→51→52→60，即误操作→生产工艺条件失常→爆炸→泄漏→燃物泄漏→火灾→池火→火焰吞噬→安全通道与出口缺陷→人员伤亡。化工事故案例统计^[19-20]中出现的主要事故影响因素和事故传播链与本结论相吻合(如四川阳春工业园区的恒达科技有限公司“7•12”重大爆炸着火事故)，验证了基于事故传播强度识别关键事故传播路径的合理性。

6 结论

通过分析化工园区事故影响因素，利用复杂网络理论建立了化工园区事故链网络模型，并从其事故链的网络拓扑结构性质出发，以层次分析法、熵权法组合赋权的TOPSIS法和引入链入链出因子的加权LeaderRank算法分别从无向和有向两个角度综合评价了化工园区事故链网络中的节点重要性，最后基于事故传播强度识别化工园区事故链网络的关键事故传播路径。经案例分析验证了该方法的可行性和有效性。

该研究方法一方面丰富了基于复杂网络研究化工园区事故传播的相关理论基础，另一方面为化工园区智慧安全建设和选择重点监测变量提供了决策支持和依据。

符号说明: