成連華，郭阿娟，劉 黎，郭慧敏

(1.西安科技大學(xué) 安全科學(xué)與工程學(xué)院，陜西 西安 710054；2.西安科技大學(xué) 西部礦井瓦斯智能抽采工程研究中心，陜西 西安 710054)

0 引 言

近年來(lái)，中國(guó)煤炭行業(yè)安全生產(chǎn)狀況持續(xù)好轉(zhuǎn)，但隨著煤炭開采向深部轉(zhuǎn)移，煤層瓦斯含量急劇升高，瓦斯爆炸風(fēng)險(xiǎn)依然存在。國(guó)家礦山安全監(jiān)察局統(tǒng)計(jì)數(shù)據(jù)顯示，2010—2020年煤礦事故死亡10 203人，因瓦斯爆炸事故死亡2 762人，約占煤礦事故死亡人數(shù)的27.10%，瓦斯爆炸風(fēng)險(xiǎn)嚴(yán)重威脅工人的生命安全。因此，對(duì)瓦斯爆炸風(fēng)險(xiǎn)演化過(guò)程定量分析，對(duì)制定預(yù)防措施具有重要現(xiàn)實(shí)意義。

目前國(guó)內(nèi)外學(xué)者圍繞煤礦瓦斯爆炸風(fēng)險(xiǎn)評(píng)價(jià)展開了大量研究。LI等基于改進(jìn)層次分析法、貝葉斯網(wǎng)絡(luò)等方法，評(píng)估了煤礦瓦斯爆炸風(fēng)險(xiǎn)[1-2]；TONG等運(yùn)用概率風(fēng)險(xiǎn)評(píng)估模型、蒙特卡羅法分析煤礦瓦斯爆炸事故中礦工不安全行為所造成的事故風(fēng)險(xiǎn)[3-4]；田水承等運(yùn)用文獻(xiàn)分析法和熵權(quán)法確定礦工不安全狀態(tài)評(píng)價(jià)指標(biāo)綜合權(quán)重，并進(jìn)行實(shí)例分析[5]；王婉青和李紅霞分別基于IAHP方法、模糊綜合評(píng)價(jià)法分析了瓦斯爆炸事故風(fēng)險(xiǎn)中的關(guān)鍵風(fēng)險(xiǎn)，同時(shí)對(duì)事故后果提出預(yù)防措施[6-7]。但大多數(shù)風(fēng)險(xiǎn)評(píng)價(jià)方法集中于對(duì)事故結(jié)果進(jìn)行分析，風(fēng)險(xiǎn)演化過(guò)程的定量分析有待進(jìn)一步研究。拓?fù)渚W(wǎng)絡(luò)是一種用來(lái)分析復(fù)雜系統(tǒng)相互作用及演化的工具，由于其契合復(fù)雜系統(tǒng)的特點(diǎn)，在城市交通[8-9]、電力系統(tǒng)[10]、房建事故[11]等諸多領(lǐng)域得到應(yīng)用。如ZHOU等構(gòu)建地鐵施工事故網(wǎng)絡(luò)(SCAN)探究26種事故之間的聯(lián)系[12]；丁明等運(yùn)用動(dòng)態(tài)故障樹理論，提出停電事故鏈網(wǎng)絡(luò)模型，并評(píng)估支路的關(guān)鍵風(fēng)險(xiǎn)程度[13]；周新宇等構(gòu)建房建事故行為致險(xiǎn)網(wǎng)絡(luò)模型，解析了人的不安全行為與事故間的作用機(jī)理[14]。然而當(dāng)前拓?fù)渚W(wǎng)絡(luò)應(yīng)用于煤礦瓦斯爆炸風(fēng)險(xiǎn)分析研究主要集中于無(wú)向無(wú)權(quán)網(wǎng)絡(luò)，對(duì)于有向加權(quán)網(wǎng)絡(luò)研究還處于探索階段。

鑒于此，運(yùn)用拓?fù)渚W(wǎng)絡(luò)，將影響瓦斯爆炸事故的風(fēng)險(xiǎn)因素作為節(jié)點(diǎn)，因素之間的演化關(guān)系作為連接邊，通過(guò)計(jì)算風(fēng)險(xiǎn)演化過(guò)程中直接和間接傳播概率，構(gòu)建煤礦瓦斯爆炸風(fēng)險(xiǎn)拓?fù)渚W(wǎng)絡(luò)度量模型。同時(shí)對(duì)節(jié)點(diǎn)出入度、聚類系數(shù)及最短路徑等拓?fù)鋮?shù)進(jìn)行分析，確定事故演化最短路徑和關(guān)鍵風(fēng)險(xiǎn)因素，以期為煤礦瓦斯爆炸風(fēng)險(xiǎn)量化分析提供參考。

1 瓦斯爆炸風(fēng)險(xiǎn)演化分析

1.1 瓦斯爆炸事故案例收集

以國(guó)家礦山安全監(jiān)察局官網(wǎng)中統(tǒng)計(jì)的煤礦瓦斯爆炸事故調(diào)查報(bào)告作為數(shù)據(jù)來(lái)源，通過(guò)搜集2010—2020年的瓦斯爆炸事故調(diào)查報(bào)告，梳理了152起瓦斯爆炸事故案例，主要包含事故基本情況、事故經(jīng)過(guò)、事故原因等基本信息。

1.2 瓦斯爆炸風(fēng)險(xiǎn)因素分析

基于層次全息建模(hierarchical holographic modeling，HHM)方法[15]從多維度、多層次進(jìn)行事故過(guò)程識(shí)別。結(jié)合安全系統(tǒng)工程理論從人、物、管理、環(huán)境、技術(shù)方法5個(gè)維度進(jìn)行分類，歸納分析出67個(gè)風(fēng)險(xiǎn)因素、2個(gè)風(fēng)險(xiǎn)事件。

統(tǒng)計(jì)瓦斯爆炸風(fēng)險(xiǎn)出現(xiàn)頻次，按公式(1)計(jì)算發(fā)生概率Pi。瓦斯爆炸事故Ck中如果出現(xiàn)風(fēng)險(xiǎn)因素類別i，則Ck=1，否則Ck=0。其中k為某個(gè)瓦斯爆炸事故案例；k={1,2,…,n}，n=152。

(1)

由于在一次瓦斯爆炸事故中可能存在多個(gè)風(fēng)險(xiǎn)，因此概率和大于1。統(tǒng)計(jì)分析瓦斯爆炸風(fēng)險(xiǎn)結(jié)果見表1。

表1 風(fēng)險(xiǎn)因素分析與編碼Table 1 Analysis and coding of risk factor

1.3 瓦斯爆炸風(fēng)險(xiǎn)演化路徑構(gòu)建

任何一起事故過(guò)程都可用一個(gè)或幾個(gè)事故鏈描述。在瓦斯爆炸風(fēng)險(xiǎn)演化過(guò)程中，當(dāng)風(fēng)險(xiǎn)節(jié)點(diǎn)狀態(tài)偏離時(shí)，風(fēng)險(xiǎn)累積達(dá)到風(fēng)險(xiǎn)閥值，觸發(fā)整個(gè)系統(tǒng)發(fā)生連鎖反應(yīng)，進(jìn)而導(dǎo)致事故發(fā)生。

依照事故發(fā)展流程，確定風(fēng)險(xiǎn)節(jié)點(diǎn)的時(shí)序關(guān)系和邏輯順序，將風(fēng)險(xiǎn)節(jié)點(diǎn)依序相連形成事故鏈[16-17]，實(shí)現(xiàn)對(duì)瓦斯爆炸風(fēng)險(xiǎn)演化過(guò)程的描述。以山西平遙二畝溝煤業(yè)瓦斯爆炸事故為例進(jìn)行分析，對(duì)事故直接原因①和②進(jìn)行分析，可先后提取出D17與B22的2個(gè)風(fēng)險(xiǎn)因素。事故經(jīng)過(guò)表明：炮采組工人張某(無(wú)爆破工特種作業(yè)人員證件)未執(zhí)行“一炮三檢”和“三人連鎖爆破”制度違章爆破，爆破明火引爆9103工作面采空區(qū)涌入瓦斯，發(fā)生瓦斯爆炸，可先后提取出A16，A17和D21的3個(gè)風(fēng)險(xiǎn)因素和F11這1個(gè)風(fēng)險(xiǎn)事件。最終從該案例中提取出8個(gè)風(fēng)險(xiǎn)因素、1個(gè)風(fēng)險(xiǎn)事件和4條事故鏈，見表2。

表2 瓦斯爆炸事故信息描述及事故鏈構(gòu)建Table 2 Information description of gas explosion accident and construction of accident chain

分析152起煤礦瓦斯爆炸事故案例，將含有共同風(fēng)險(xiǎn)因素的事故鏈融合到一個(gè)全局網(wǎng)絡(luò)中[18]，運(yùn)用系統(tǒng)動(dòng)力學(xué)軟件構(gòu)建煤礦瓦斯爆炸風(fēng)險(xiǎn)演化路徑，如圖1所示。

圖1 煤礦瓦斯爆炸風(fēng)險(xiǎn)演化路徑Fig.1 Risk evolution path of gas explosion in coal mine

2 瓦斯爆炸風(fēng)險(xiǎn)度量模型構(gòu)建

2.1 風(fēng)險(xiǎn)網(wǎng)絡(luò)連邊計(jì)算

引入交互作用矩陣(factor interaction matrix，F(xiàn)IM)理論計(jì)算風(fēng)險(xiǎn)連邊權(quán)重，以分析系統(tǒng)影響因素的交互作用關(guān)系，如圖2所示。

圖2 因子交互作用矩陣

因子交互矩陣共有n個(gè)風(fēng)險(xiǎn)因素，主對(duì)角上各風(fēng)險(xiǎn)因素表示為Pi(i=1,2,…,n)，非對(duì)角上方區(qū)域表示對(duì)角線上其他風(fēng)險(xiǎn)因素對(duì)相應(yīng)Pi的影響xij，非對(duì)角下方區(qū)域表示對(duì)角線上Pi對(duì)其他各風(fēng)險(xiǎn)因素的影響xji。每一行對(duì)應(yīng)的數(shù)值相加表示為Ci，每一列的對(duì)應(yīng)數(shù)值相加表示為Ei，則風(fēng)險(xiǎn)網(wǎng)絡(luò)連邊權(quán)重wi的計(jì)算表達(dá)式為

(2)

由于文章篇幅所限，以安全檢查不到位→線頭外露→電路老化→電氣火花→引爆瓦斯(D19→B14→B11→B22→F11)風(fēng)險(xiǎn)演化路徑為例進(jìn)行計(jì)算，由FIM理論計(jì)算可得風(fēng)險(xiǎn)網(wǎng)絡(luò)連邊權(quán)重wi=(0.55，0.48，0.35，1.08)。

2.2 風(fēng)險(xiǎn)網(wǎng)絡(luò)迭代加權(quán)

設(shè)計(jì)結(jié)構(gòu)矩陣(design structure matrix，DSM)是一種用來(lái)表示系統(tǒng)中元素及相互作用的結(jié)構(gòu)化建模方法，反映風(fēng)險(xiǎn)傳播網(wǎng)絡(luò)中一個(gè)元素通過(guò)所有可能路徑影響另一個(gè)元素的集成風(fēng)險(xiǎn)傳播概率，可直觀對(duì)復(fù)雜系統(tǒng)進(jìn)行可視化分析[19-20]。筆者將迭代加權(quán)思想與DSM相結(jié)合，構(gòu)建集成風(fēng)險(xiǎn)傳播DSM矩陣。

實(shí)際上，風(fēng)險(xiǎn)傳播不僅存在直接傳播，也存在多條路徑的間接傳播關(guān)系。R_DSM(Ri,Rj)為Rj對(duì)Ri的直接傳播概率(傳播路徑是1階)。

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

(3)

Rj對(duì)Ri的2階路徑風(fēng)險(xiǎn)傳播(經(jīng)過(guò)中間風(fēng)險(xiǎn)因素Rp)概率為

(4)

Rj對(duì)Ri的3階路徑風(fēng)險(xiǎn)傳播(經(jīng)過(guò)中間風(fēng)險(xiǎn)因素Rp和Rq)概率為

(5)

風(fēng)險(xiǎn)演化過(guò)程中可能存在多條1階、2階或高階傳播路徑[21]。綜合考慮所有風(fēng)險(xiǎn)傳播路徑，得到Rj對(duì)Ri的集成風(fēng)險(xiǎn)傳播概率CP為

(6)

以上述風(fēng)險(xiǎn)網(wǎng)絡(luò)連邊計(jì)算的風(fēng)險(xiǎn)演化路徑(D19→B14→B11→B22→F11)為例進(jìn)行分析，RD19直接影響下一階RB14，通過(guò)RB14間接影響RB11(路徑長(zhǎng)度為2)，如圖3(a)所示。由于RD19導(dǎo)致RB14發(fā)生概率為0.55，而RB14導(dǎo)致下一階RB11發(fā)生概率為0.48，由式(4)計(jì)算得RD19與RB11風(fēng)險(xiǎn)間接(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 風(fēng)險(xiǎn)集成迭代加和概率Fig.3 Risk-integrated iterative summation probability

2.3 風(fēng)險(xiǎn)拓?fù)渚W(wǎng)絡(luò)形成

將風(fēng)險(xiǎn)迭代加和權(quán)重轉(zhuǎn)化成CSV文件，并運(yùn)用Gephi 0.9.2軟件構(gòu)建煤礦瓦斯爆炸風(fēng)險(xiǎn)加權(quán)拓?fù)渚W(wǎng)絡(luò)，如圖4所示。

圖4 風(fēng)險(xiǎn)加權(quán)拓?fù)渚W(wǎng)絡(luò)模型Fig.4 Risk-weighted topological network model

從圖4可以看出，煤礦瓦斯爆炸風(fēng)險(xiǎn)拓?fù)渚W(wǎng)絡(luò)模型是一個(gè)完整的連通網(wǎng)絡(luò)。圓圈代表67個(gè)風(fēng)險(xiǎn)節(jié)點(diǎn)，連線代表214條網(wǎng)絡(luò)連邊，連線的粗細(xì)程度代表節(jié)點(diǎn)間的集成傳播概率。如B25線條最粗，可直觀看出與該風(fēng)險(xiǎn)鏈接路徑最多，風(fēng)險(xiǎn)權(quán)值變化最大。通過(guò)定量描述瓦斯爆炸風(fēng)險(xiǎn)加權(quán)網(wǎng)絡(luò)中風(fēng)險(xiǎn)因素之間的復(fù)雜關(guān)系，為煤礦瓦斯爆炸風(fēng)險(xiǎn)度量提供了一個(gè)新的視角。

3 拓?fù)渚W(wǎng)絡(luò)風(fēng)險(xiǎn)度量分析

3.1 風(fēng)險(xiǎn)節(jié)點(diǎn)度分析

節(jié)點(diǎn)度分析是網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)最重要的指標(biāo)之一，反映了節(jié)點(diǎn)在網(wǎng)絡(luò)中的直接影響力和重要性[22]。節(jié)點(diǎn)的度值越大，表明節(jié)點(diǎn)越重要。節(jié)點(diǎn)的度可表示為

(7)

通過(guò)Gephi 0.9.2軟件分析得到節(jié)點(diǎn)度分布如圖5所示，橫坐標(biāo)為度貢獻(xiàn)值，縱坐標(biāo)為度分布數(shù)量。平均度為2.315，表明每個(gè)風(fēng)險(xiǎn)因素至少與2個(gè)不安全因素存在鏈接關(guān)系。

圖5 節(jié)點(diǎn)度分布Fig.5 Distribution of node degree

選取風(fēng)險(xiǎn)節(jié)點(diǎn)中出、入度數(shù)排前的節(jié)點(diǎn)共20個(gè)，度分布如圖6所示。從圖6可以看出，B20，B34的入度和出度最高，表明誘發(fā)該風(fēng)險(xiǎn)的途徑較多，因此對(duì)該風(fēng)險(xiǎn)進(jìn)行有效控制，抑制瓦斯爆炸發(fā)生的效果最佳。同理，A12，B17，B24等度分布較高的風(fēng)險(xiǎn)節(jié)點(diǎn)也是風(fēng)險(xiǎn)防控中需重點(diǎn)關(guān)注對(duì)象。

圖6 風(fēng)險(xiǎn)節(jié)點(diǎn)出、入度統(tǒng)計(jì)Fig.6 Outbound and inbound statistics of risk nodes

3.2 風(fēng)險(xiǎn)聚類分析

聚類系數(shù)是網(wǎng)絡(luò)節(jié)點(diǎn)中鄰接節(jié)點(diǎn)之間相互連接的比例，反映了網(wǎng)絡(luò)節(jié)點(diǎn)集團(tuán)化程度。對(duì)于與k個(gè)節(jié)點(diǎn)相連的節(jié)點(diǎn)i而言，它最多可與這些鄰居節(jié)點(diǎn)構(gòu)成k(k-1)/2條連邊，則實(shí)際存在連邊Ei與最多可能連邊個(gè)數(shù)之比是節(jié)點(diǎn)i的聚類系數(shù)Ci。

Ci=2Ei/k(k-1)

(8)

通過(guò)Gephi 0.9.2軟件分析得到風(fēng)險(xiǎn)網(wǎng)絡(luò)模型的聚類系數(shù)為0.212，表明在整個(gè)網(wǎng)絡(luò)演化中，導(dǎo)致瓦斯爆炸的成因復(fù)雜，但大部分風(fēng)險(xiǎn)因素之間的相互影響并不明顯。

篩選節(jié)點(diǎn)加權(quán)聚類系數(shù)排前的節(jié)點(diǎn)共20個(gè)，如圖7所示。從圖7可以看出，B20，B33聚類系數(shù)最高，表明該節(jié)點(diǎn)集團(tuán)化嚴(yán)重，若對(duì)該節(jié)點(diǎn)及鄰接節(jié)點(diǎn)構(gòu)成的子系統(tǒng)采取斷鏈控制措施，可有效減緩事故發(fā)生。

圖7 節(jié)點(diǎn)加權(quán)聚類系數(shù)Fig.7 Node-weighted clustering coefficients

3.3 最短路徑分析

節(jié)點(diǎn)度分布及聚類系數(shù)分別分析了風(fēng)險(xiǎn)節(jié)點(diǎn)的重要程度及聚類程度，但未考慮上一節(jié)點(diǎn)引發(fā)后續(xù)節(jié)點(diǎn)的風(fēng)險(xiǎn)傳遞。通過(guò)Gephi 0.9.2軟件分析得到平均最短路徑為2.506，表明瓦斯爆炸風(fēng)險(xiǎn)拓?fù)渚W(wǎng)絡(luò)中一個(gè)節(jié)點(diǎn)狀態(tài)發(fā)生變化平均僅通過(guò)2.506步就能夠引起鄰接其它節(jié)點(diǎn)發(fā)生變化。

Dijkstra算法是荷蘭著名計(jì)算機(jī)科學(xué)家E.W.Dijkstra于1959年提出關(guān)于圖論中求最短路徑常用算法[23]。借鑒Dijkstra算法，編制計(jì)算程序，求解瓦斯爆炸風(fēng)險(xiǎn)最短路徑，得到風(fēng)險(xiǎn)節(jié)點(diǎn)的最短路徑及其最短路徑值(風(fēng)險(xiǎn)權(quán)值)，見表3。

表3 瓦斯爆炸風(fēng)險(xiǎn)演化最短路徑分析Table 3 Analysis of the shortest path in risk evolution of gas explosion

通過(guò)分析表3可知，風(fēng)險(xiǎn)演化路徑越短，權(quán)重值越大，但在風(fēng)險(xiǎn)路徑傳遞過(guò)程中風(fēng)險(xiǎn)權(quán)值逐漸減少，說(shuō)明該過(guò)程存在風(fēng)險(xiǎn)損失；B31→F12風(fēng)險(xiǎn)路徑最短，風(fēng)險(xiǎn)權(quán)重變化顯著，表明B31(沖擊波蔓延)對(duì)瓦斯爆炸風(fēng)險(xiǎn)演化影響較大；雖然B11→B16→B25、A20→B20→F11和B17→B16→B25這3條最短演化路徑都是由3個(gè)風(fēng)險(xiǎn)因素構(gòu)成，但初始風(fēng)險(xiǎn)因素不同，風(fēng)險(xiǎn)演化路徑不同，所以風(fēng)險(xiǎn)權(quán)值也不同。此外，風(fēng)險(xiǎn)演化最短路徑也從側(cè)面反映了瓦斯爆炸事故發(fā)生的可能性。在實(shí)際瓦斯爆炸風(fēng)險(xiǎn)管控工作中，優(yōu)先對(duì)最短風(fēng)險(xiǎn)路徑上的風(fēng)險(xiǎn)因素加以控制，可以提高解決事故的效率。

4 結(jié) 論

1)通過(guò)統(tǒng)計(jì)分析2010—2020年共計(jì)152起瓦斯爆炸事故案例，采用HHM方法從人、物、管理、環(huán)境、技術(shù)方法5個(gè)維度進(jìn)行風(fēng)險(xiǎn)分析，得到67個(gè)風(fēng)險(xiǎn)因素，形成214條瓦斯爆炸風(fēng)險(xiǎn)演化路徑鏈。

2)以拓?fù)渚W(wǎng)絡(luò)算法為基礎(chǔ)，引入DSM對(duì)瓦斯爆炸風(fēng)險(xiǎn)進(jìn)行迭代加權(quán)計(jì)算，并運(yùn)用Gephi 0.9.2軟件構(gòu)建瓦斯爆炸風(fēng)險(xiǎn)拓?fù)渚W(wǎng)絡(luò)度量模型，實(shí)現(xiàn)瓦斯爆炸風(fēng)險(xiǎn)演化過(guò)程的量化分析，以揭示事故發(fā)生的潛在特征。

3)通過(guò)分析瓦斯爆炸風(fēng)險(xiǎn)有向帶權(quán)網(wǎng)絡(luò)參數(shù)，發(fā)現(xiàn)風(fēng)險(xiǎn)節(jié)點(diǎn)聚集程度較低，呈現(xiàn)明顯小世界網(wǎng)絡(luò)特征，表明在整個(gè)風(fēng)險(xiǎn)演化過(guò)程中大部分風(fēng)險(xiǎn)因素之間相互影響不明顯，存在顯著傳遞關(guān)系。

4)運(yùn)用Dijkstra算法搜尋風(fēng)險(xiǎn)演化最短路徑，得到7條最短風(fēng)險(xiǎn)演化路徑和3個(gè)關(guān)鍵風(fēng)險(xiǎn)因素，優(yōu)先對(duì)瓦斯爆炸最短風(fēng)險(xiǎn)路徑上關(guān)鍵風(fēng)險(xiǎn)加以管控，對(duì)預(yù)防瓦斯爆炸事故發(fā)生、保障企業(yè)安全運(yùn)行具有現(xiàn)實(shí)意義。