基于拓扑网络算法的煤矿瓦斯爆炸风险度量

成连华，郭阿娟，刘 黎，郭慧敏

(1.西安科技大学 安全科学与工程学院，陕西 西安 710054；2.西安科技大学 西部矿井瓦斯智能抽采工程研究中心，陕西 西安 710054)

0 引 言

近年来，中国煤炭行业安全生产状况持续好转，但随着煤炭开采向深部转移，煤层瓦斯含量急剧升高，瓦斯爆炸风险依然存在。国家矿山安全监察局统计数据显示，2010—2020年煤矿事故死亡10 203人，因瓦斯爆炸事故死亡2 762人，约占煤矿事故死亡人数的27.10%，瓦斯爆炸风险严重威胁工人的生命安全。因此，对瓦斯爆炸风险演化过程定量分析，对制定预防措施具有重要现实意义。

目前国内外学者围绕煤矿瓦斯爆炸风险评价展开了大量研究。LI等基于改进层次分析法、贝叶斯网络等方法，评估了煤矿瓦斯爆炸风险[1-2]；TONG等运用概率风险评估模型、蒙特卡罗法分析煤矿瓦斯爆炸事故中矿工不安全行为所造成的事故风险[3-4]；田水承等运用文献分析法和熵权法确定矿工不安全状态评价指标综合权重，并进行实例分析[5]；王婉青和李红霞分别基于IAHP方法、模糊综合评价法分析了瓦斯爆炸事故风险中的关键风险，同时对事故后果提出预防措施[6-7]。但大多数风险评价方法集中于对事故结果进行分析，风险演化过程的定量分析有待进一步研究。拓扑网络是一种用来分析复杂系统相互作用及演化的工具，由于其契合复杂系统的特点，在城市交通[8-9]、电力系统[10]、房建事故[11]等诸多领域得到应用。如ZHOU等构建地铁施工事故网络(SCAN)探究26种事故之间的联系[12]；丁明等运用动态故障树理论，提出停电事故链网络模型，并评估支路的关键风险程度[13]；周新宇等构建房建事故行为致险网络模型，解析了人的不安全行为与事故间的作用机理[14]。然而当前拓扑网络应用于煤矿瓦斯爆炸风险分析研究主要集中于无向无权网络，对于有向加权网络研究还处于探索阶段。

鉴于此，运用拓扑网络，将影响瓦斯爆炸事故的风险因素作为节点，因素之间的演化关系作为连接边，通过计算风险演化过程中直接和间接传播概率，构建煤矿瓦斯爆炸风险拓扑网络度量模型。同时对节点出入度、聚类系数及最短路径等拓扑参数进行分析，确定事故演化最短路径和关键风险因素，以期为煤矿瓦斯爆炸风险量化分析提供参考。

1 瓦斯爆炸风险演化分析

1.1 瓦斯爆炸事故案例收集

以国家矿山安全监察局官网中统计的煤矿瓦斯爆炸事故调查报告作为数据来源，通过搜集2010—2020年的瓦斯爆炸事故调查报告，梳理了152起瓦斯爆炸事故案例，主要包含事故基本情况、事故经过、事故原因等基本信息。

1.2 瓦斯爆炸风险因素分析

基于层次全息建模(hierarchical holographic modeling，HHM)方法[15]从多维度、多层次进行事故过程识别。结合安全系统工程理论从人、物、管理、环境、技术方法5个维度进行分类，归纳分析出67个风险因素、2个风险事件。

统计瓦斯爆炸风险出现频次，按公式(1)计算发生概率Pi。瓦斯爆炸事故Ck中如果出现风险因素类别i，则Ck=1，否则Ck=0。其中k为某个瓦斯爆炸事故案例；k={1,2,…,n}，n=152。

(1)

由于在一次瓦斯爆炸事故中可能存在多个风险，因此概率和大于1。统计分析瓦斯爆炸风险结果见表1。

表1 风险因素分析与编码Table 1 Analysis and coding of risk factor

1.3 瓦斯爆炸风险演化路径构建

任何一起事故过程都可用一个或几个事故链描述。在瓦斯爆炸风险演化过程中，当风险节点状态偏离时，风险累积达到风险阀值，触发整个系统发生连锁反应，进而导致事故发生。

依照事故发展流程，确定风险节点的时序关系和逻辑顺序，将风险节点依序相连形成事故链[16-17]，实现对瓦斯爆炸风险演化过程的描述。以山西平遥二亩沟煤业瓦斯爆炸事故为例进行分析，对事故直接原因①和②进行分析，可先后提取出D17与B22的2个风险因素。事故经过表明：炮采组工人张某(无爆破工特种作业人员证件)未执行“一炮三检”和“三人连锁爆破”制度违章爆破，爆破明火引爆9103工作面采空区涌入瓦斯，发生瓦斯爆炸，可先后提取出A16，A17和D21的3个风险因素和F11这1个风险事件。最终从该案例中提取出8个风险因素、1个风险事件和4条事故链，见表2。

表2 瓦斯爆炸事故信息描述及事故链构建Table 2 Information description of gas explosion accident and construction of accident chain

分析152起煤矿瓦斯爆炸事故案例，将含有共同风险因素的事故链融合到一个全局网络中[18]，运用系统动力学软件构建煤矿瓦斯爆炸风险演化路径，如图1所示。

图1 煤矿瓦斯爆炸风险演化路径Fig.1 Risk evolution path of gas explosion in coal mine

2 瓦斯爆炸风险度量模型构建

2.1 风险网络连边计算

引入交互作用矩阵(factor interaction matrix，FIM)理论计算风险连边权重，以分析系统影响因素的交互作用关系，如图2所示。

图2 因子交互作用矩阵

因子交互矩阵共有n个风险因素，主对角上各风险因素表示为Pi(i=1,2,…,n)，非对角上方区域表示对角线上其他风险因素对相应Pi的影响xij，非对角下方区域表示对角线上Pi对其他各风险因素的影响xji。每一行对应的数值相加表示为Ci，每一列的对应数值相加表示为Ei，则风险网络连边权重wi的计算表达式为

(2)

由于文章篇幅所限，以安全检查不到位→线头外露→电路老化→电气火花→引爆瓦斯(D19→B14→B11→B22→F11)风险演化路径为例进行计算，由FIM理论计算可得风险网络连边权重wi=(0.55，0.48，0.35，1.08)。

2.2 风险网络迭代加权

设计结构矩阵(design structure matrix，DSM)是一种用来表示系统中元素及相互作用的结构化建模方法，反映风险传播网络中一个元素通过所有可能路径影响另一个元素的集成风险传播概率，可直观对复杂系统进行可视化分析[19-20]。笔者将迭代加权思想与DSM相结合，构建集成风险传播DSM矩阵。

实际上，风险传播不仅存在直接传播，也存在多条路径的间接传播关系。R_DSM(Ri,Rj)为Rj对Ri的直接传播概率(传播路径是1阶)。

SP(1)(Ri,Rj)=R_DSM(Ri,Rj)

(3)

Rj对Ri的2阶路径风险传播(经过中间风险因素Rp)概率为

(4)

Rj对Ri的3阶路径风险传播(经过中间风险因素Rp和Rq)概率为

(5)

风险演化过程中可能存在多条1阶、2阶或高阶传播路径[21]。综合考虑所有风险传播路径，得到Rj对Ri的集成风险传播概率CP为

(6)

以上述风险网络连边计算的风险演化路径(D19→B14→B11→B22→F11)为例进行分析，RD19直接影响下一阶RB14，通过RB14间接影响RB11(路径长度为2)，如图3(a)所示。由于RD19导致RB14发生概率为0.55，而RB14导致下一阶RB11发生概率为0.48，由式(4)计算得RD19与RB11风险间接(2阶)传播概率为：R_DSM(RB11,RD19)=0.48×0.55=0.264。演化路径R1→R2→…→R5集成传播概率为：CP(R5,R1)=1-(1-0.26×0.35)(1-0.55×0.48×1.08)=0.35。如图3(b)所示。

图3 风险集成迭代加和概率Fig.3 Risk-integrated iterative summation probability

2.3 风险拓扑网络形成

将风险迭代加和权重转化成CSV文件，并运用Gephi 0.9.2软件构建煤矿瓦斯爆炸风险加权拓扑网络，如图4所示。

图4 风险加权拓扑网络模型Fig.4 Risk-weighted topological network model

从图4可以看出，煤矿瓦斯爆炸风险拓扑网络模型是一个完整的连通网络。圆圈代表67个风险节点，连线代表214条网络连边，连线的粗细程度代表节点间的集成传播概率。如B25线条最粗，可直观看出与该风险链接路径最多，风险权值变化最大。通过定量描述瓦斯爆炸风险加权网络中风险因素之间的复杂关系，为煤矿瓦斯爆炸风险度量提供了一个新的视角。

3 拓扑网络风险度量分析

3.1 风险节点度分析

节点度分析是网络拓扑结构最重要的指标之一，反映了节点在网络中的直接影响力和重要性[22]。节点的度值越大，表明节点越重要。节点的度可表示为

(7)

通过Gephi 0.9.2软件分析得到节点度分布如图5所示，横坐标为度贡献值，纵坐标为度分布数量。平均度为2.315，表明每个风险因素至少与2个不安全因素存在链接关系。

图5 节点度分布Fig.5 Distribution of node degree

选取风险节点中出、入度数排前的节点共20个，度分布如图6所示。从图6可以看出，B20，B34的入度和出度最高，表明诱发该风险的途径较多，因此对该风险进行有效控制，抑制瓦斯爆炸发生的效果最佳。同理，A12，B17，B24等度分布较高的风险节点也是风险防控中需重点关注对象。

图6 风险节点出、入度统计Fig.6 Outbound and inbound statistics of risk nodes

3.2 风险聚类分析

聚类系数是网络节点中邻接节点之间相互连接的比例，反映了网络节点集团化程度。对于与k个节点相连的节点i而言，它最多可与这些邻居节点构成k(k-1)/2条连边，则实际存在连边Ei与最多可能连边个数之比是节点i的聚类系数Ci。

Ci=2Ei/k(k-1)

(8)

通过Gephi 0.9.2软件分析得到风险网络模型的聚类系数为0.212，表明在整个网络演化中，导致瓦斯爆炸的成因复杂，但大部分风险因素之间的相互影响并不明显。

筛选节点加权聚类系数排前的节点共20个，如图7所示。从图7可以看出，B20，B33聚类系数最高，表明该节点集团化严重，若对该节点及邻接节点构成的子系统采取断链控制措施，可有效减缓事故发生。

图7 节点加权聚类系数Fig.7 Node-weighted clustering coefficients

3.3 最短路径分析

节点度分布及聚类系数分别分析了风险节点的重要程度及聚类程度，但未考虑上一节点引发后续节点的风险传递。通过Gephi 0.9.2软件分析得到平均最短路径为2.506，表明瓦斯爆炸风险拓扑网络中一个节点状态发生变化平均仅通过2.506步就能够引起邻接其它节点发生变化。

Dijkstra算法是荷兰著名计算机科学家E.W.Dijkstra于1959年提出关于图论中求最短路径常用算法[23]。借鉴Dijkstra算法，编制计算程序，求解瓦斯爆炸风险最短路径，得到风险节点的最短路径及其最短路径值(风险权值)，见表3。

表3 瓦斯爆炸风险演化最短路径分析Table 3 Analysis of the shortest path in risk evolution of gas explosion

通过分析表3可知，风险演化路径越短，权重值越大，但在风险路径传递过程中风险权值逐渐减少，说明该过程存在风险损失；B31→F12风险路径最短，风险权重变化显著，表明B31(冲击波蔓延)对瓦斯爆炸风险演化影响较大；虽然B11→B16→B25、A20→B20→F11和B17→B16→B25这3条最短演化路径都是由3个风险因素构成，但初始风险因素不同，风险演化路径不同，所以风险权值也不同。此外，风险演化最短路径也从侧面反映了瓦斯爆炸事故发生的可能性。在实际瓦斯爆炸风险管控工作中，优先对最短风险路径上的风险因素加以控制，可以提高解决事故的效率。

4 结 论

1)通过统计分析2010—2020年共计152起瓦斯爆炸事故案例，采用HHM方法从人、物、管理、环境、技术方法5个维度进行风险分析，得到67个风险因素，形成214条瓦斯爆炸风险演化路径链。

2)以拓扑网络算法为基础，引入DSM对瓦斯爆炸风险进行迭代加权计算，并运用Gephi 0.9.2软件构建瓦斯爆炸风险拓扑网络度量模型，实现瓦斯爆炸风险演化过程的量化分析，以揭示事故发生的潜在特征。

3)通过分析瓦斯爆炸风险有向带权网络参数，发现风险节点聚集程度较低，呈现明显小世界网络特征，表明在整个风险演化过程中大部分风险因素之间相互影响不明显，存在显著传递关系。

4)运用Dijkstra算法搜寻风险演化最短路径，得到7条最短风险演化路径和3个关键风险因素，优先对瓦斯爆炸最短风险路径上关键风险加以管控，对预防瓦斯爆炸事故发生、保障企业安全运行具有现实意义。