基于權值時變模型的礦井突水最優逃生路徑的動態選擇

于丹， 顏偉,2*， 李劭昱

(1.山東科技大學能源與礦業工程學院， 青島 266590； 2.山東科技大學礦山災害預防控制省部共建國家重點實驗室培育基地， 青島 266590)

在煤礦建設和采礦生產中，礦井突水是僅次于瓦斯爆炸的主要突發性災害之一。盡管礦井突水監測及預警技術研究及應用已進行多年，但隨著煤礦開采深度的逐步增加，采煤環境的水文地質條件將更加復雜多變，礦井突水事故仍時有發生[1]。礦井突水一般來勢兇猛，不僅會在短時間內淹沒煤礦巷道和生產設備，破壞生態環境，還會造成人員重大傷亡。因此，為了讓井下礦工在發生突水之后，快速、準確地選擇出最優的逃生路徑(即從逃生地到安全地所用時間最短、危險性最低的安全撤離路線)[2]是應對礦井突水事故的關鍵。

當前，關于礦井突水逃生路徑選擇方面的相關研究主要分為兩個方面：一是在研究影響井下礦工逃生因素方面。除考慮突水災害水位高度對逃生影響外，還應考慮巷道類型、風速風向、坡度、局部障礙物數量[3]以及是否有交通工具[4]等因素，這些影響因素使得巷道節點之間的距離不再是絕對距離而需要用當量長度來表示；二是在逃生路徑選擇模型構建及算法方面。主要是對傳統算法的改進以及朝著建立更符合現實情況的動態逃生路徑選擇模型這兩個方向進行研究。蔡明杰等[5-6]通過對交叉、變異算子進行重構，在適應度函數中加入安全通過概率來改進傳統的遺傳算法；通過對松弛點入隊位置調整來優化SPFA(shortest path faster algorithm)算法，從而提高算法搜索突水救援路線的效率。劉夢杰等[7]提出了一種新型路徑搜索算法——雙向A*算法，根據評估函數來確定搜索方向，使得選擇的逃生路徑節點更少，大大提高了算法的搜索效率。Yan等[8]對蟻群算法進行了改進：一是利用巷道分區特性影響螞蟻的行為，提高螞蟻的搜索效率；二是修改了蟻群會合策略和蟻群死亡規則，以提高蟻群算法的性能。王鵬等[9]將改進螢火蟲算法應用到突水環境的路徑規劃中，并與A*算法進行比較，結果表明：改進后的螢火蟲算法選擇的路徑更優，具有大范圍搜索優化能力。周越等[10]為研究實時水位高度變化對逃生路徑選擇影響，提出時間當量長度模型，并通過改進Dijkstra算法進行求解。Zhao等[11]建立了三維礦井巷道網絡，結合實時變化的水位及巷道的可通行性，運用Dijkstra算法實現了基于3D網絡模型的最佳逃生路徑搜索。Du等[12]、Wu等[13]考慮實時水位高度變化對逃生時間的影響，通過改進Dijkstra算法，提出了并行時變最早到達算法；實現人員逃生智能動態路徑搜索，生成安全迅速的井下人員最優逃生路徑。

綜上可知，雖然學者們研究了基于實時災害影響的逃生路徑選擇模型及算法，但在考慮突水災害因素方面只考慮水位高度，而實際情況是不同時刻巷道內的水流速同樣對逃生路徑選擇影響很大。為此，提出一種同時考慮水位高度和水流速的礦井突水最優路徑選擇模型，將所研究的時間進行分段，分別估算各時間段內各巷道的水位高度和水流速，用于判斷不同時段內各巷道的危險性，再將其加權到巷道初始當量長度中形成權值時變模型，并通過改進的Dijkstra算法來對權值時變模型進行求解，從而動態選擇出安全可靠的最優逃生路徑。

1 最優逃生路徑選擇原則分析

根據實際礦井突水應急避災的情況，結合煤礦安全管理，提出礦井突水災害逃生路徑選擇的兩大原則：安全性原則、時間最短原則。

1.1 安全性原則

由于礦井突水發生后，大量的地下水會迅速涌入巷道。水流作為最主要的災害因素，其高度及流速嚴重影響工作人員的逃生效率，因此，為了井下礦工的人身安全和提高逃生效率，在選擇路徑時，盡量避開有水災的路徑或者選擇水位較低的路徑。

1.2 時間最短原則

在礦井突水發生的情況下，受災礦工生還的可能性與逃生時間成反比，逃生時間越長，生還的可能性越低。因此，需要找出花費時間最短的路徑，而最短距離提供的路徑，其用時不一定最少，所以在研究這方面內容時，逃生的時效性一般以巷道的“當量長度”衡量,即找到距離安全出口當量長度最短的路徑[14]。

2 最優逃生路徑模型構建

2.1 模型整體設計

礦井突水發生后，由于受到水流因素的影響，井下礦工在逃生時不能隨意選擇逃生路徑，而應該結合巷道內的環境來選擇。依據水流漫延的方向，將巷道網絡圖由無向圖變成有向圖，有利于確定正確的疏散方向，進而排除無效疏散巷道；一些巷道被水淹沒或者危險系數過高而不能通行的巷道，也應該在搜索路徑之前排除；在計算出巷道的危險系數的基礎上，再考慮巷道固有環境對礦工通行的影響。據此，本文綜合考慮了突水水流和巷道固有環境兩類因素，為井下礦工設計了一種礦井突水最優逃生路徑模型，如圖1所示。

2.2 突水漫延方向

發生突水時，水流會在水力梯度的作用下迅速向周圍連通的巷道漫延，若突水點高于相鄰巷道高程時，會發生下向漫延；當突水點以下的巷道充滿水之后，水流就會沿著突水點之上的位置開始上向升漲。根據水流漫延的方向，將巷道網絡圖由無向圖變成有向圖，有利于確定正確的逃生方向。由于在水流下向漫延過程中，水流受到重力作用影響自動找到最低的巷道位置，所以從巷道突水點到這樣一段的巷道不應該作為選擇的逃生路徑。

圖1 礦井突水最優逃生路徑模型Fig.1 Optimal escape path model of mine water inrush

2.3 危險系數的計算

水流是影響逃生路徑選擇的最重要的動態因素，這里將水位高度和水的流速對人體穩定性的影響作為求解危險系數的影響因子。

2.3.1 計算巷道內水位高度

當已知突水流量和時間時，利用插值計算得出對應時間段內的水位高度(不考慮圍巖滲入量和設施排水量)，計算公式為[15]

(1)

(2)

式中：lx為流入巷道的水流長度，m；Q為突水流量(恒定)，m3/s；t為突水時間，s；S為巷道橫截面積，m2；hx為水位高度，m；h0為起點節點標高，m；h1為終點節點標高，m；l為涌入水流的巷道總長度，m。

基于上述所求的水位高度，來估算每段時間不同巷道內的水位上漲速度，假設當前時刻為0，t為突水時間(以秒為單位)，將突水時間進行時段劃分為T0

2.3.2 計算巷道內水流速

由于風速和障礙物對巷道內水流速的影響比較抽象，難以用數學方法進行定量計算，因此只考慮巷道的種類對水流速的影響。在不同類型的巷道中，由于巷道斷面、井巷坡度、巷道壁粗糙度等因素有所差異，水的流速是不同的且隨著時間變化。

根據文獻[16]分別求出水平、豎直、傾斜這3種巷道類型中的水流速，具體計算公式如下。

(1)水平巷道水流速。

(3)

式(3)中：g為重力加速度；λ為沿程阻力系數；hf為沿程水頭損失；R為水力半徑，R=A/χ，其中，A為過水斷面面積，χ為濕周，χ=b+2h，b為巷道底部寬度，h為水位高度；C為謝才系數，C=R1/6/n，其中n為粗糙系數，n取0.02；J為水力坡度。

(2)豎直巷道水流速。

(4)

式(4)中：m為水流單位質量；hvertical為豎井高度。

(3)傾斜巷道水流速。若斜井傾斜角度小于45°時，由水平巷道水流速度計算公式結合傾斜角度求得；若斜井傾斜角度大于等于45°時，由豎直巷道水流速度計算公式結合傾斜角度求得。

2.3.3 巷道危險性判斷

首先根據Xia等[17]通過研究洪水的水位高度及水流速對人體發生跌倒失穩的影響得到臨界流速Uc計算公式為

(5)

式(5)中：d為水位高度；ρf為水的密度，取1.0×103kg/m3；hp、mp分別為人的身高及體重(取礦工們的平均身高和平均體重分別為1.70 m、70 kg)；由中國成年人身體結構特征率定出經驗參數a1=0.633、b1=0.367、a2=1.015×10-3m3/kg、b2=-4.927×10-3m3；參數α和β與人體形狀有關，對于其率定，使用Karvonen等[18]提出的真實人體跌倒失穩的試驗方法得到α=7.867 m0.5/s，β=0.462。

然后根據巷道內水的平均流速與臨界流速的比值可以計算人體在每條巷道的危險系數Sij，計算公式為

Sij=Uij/Ucij

(6)

式(6)中：Uij為某時間段巷道Eij內的水的平均流速，m/s。Ucij為某時間段人體發生跌倒失穩時的巷道Eij內臨界流速，m/s；Sij為人在有水流的巷道Eij中逃生的危險系數。

最后根據Sij值將人體在巷道中通行的危險程度劃分為4個等級(參考文獻[19]中洪水對人的危險等級進行劃分)如表1所示。

表1 危險等級Table 1 Risk grade

2.4 巷道初始當量長度的計算

井下環境具有特殊性，含有多種阻礙礦工通行的因素，如巷道類型、風速風向、坡度、障礙物、照明情況、是否有交通工具等都會對逃生產生影響，相關影響因素可以轉化為巷道初始當量長度來表示。

2.4.1 巷道影響因子的確定

結合礦井的巷道環境以及現場模擬分析，選取巷道類型、巷道坡度、風速風向、局部障礙物的數量，作為影響井下礦工逃生的主要靜態影響因素。將巷道類型、巷道坡度、風速風向、局部障礙物的數量對巷道初始當量長度的影響系數分別記為ηt、ηs、ηv、ηr。影響系數的計算公式為[20]

(7)

式(7)中：η(Eij)為某因素對巷道Eij初始當量長度的影響系數；T(Eij)為受到某因素影響時通過巷道Eij的通行時間；t(Eij)為不受某因素影響時通過巷道Eij的通行時間。

根據式(7)所求得各因素對巷道初始當量長度的影響系數的計算結果如表2～表5所示。

表2 巷道類型對巷道初始當量長度的影響系數ηtTable 2 Influence coefficient ηt of roadway type on initial equivalent length of roadway

表4 風速風向對巷道初始當量長度的影響系數ηvTable 4 Influence coefficient ηv of wind speed and direction on initial equivalent length of roadway

表5 局部障礙物對巷道初始當量長度的影響系數ηrTable 5 Influence coefficient ηr of local obstacles on initial equivalent length of roadway

2.4.2 巷道初始當量長度的確定

假設節點Vi和Vj間巷道Eij的實際長度為l(Eij)，將巷道的各靜態因素影響系數加權到其對應的實際長度中，得到該條巷道初始當量長度值L(Eij)為

(8)

式(8)中：n為巷道Eij上障礙物數量；ηrm為第m個局部通行障礙物影響巷道正常通行的系數。

2.5 目標函數

作為應急逃生方案的決策者，在對逃生路徑進行選擇時，一般按照當量長度最小的路徑前進，所用時間最短，成功逃生的概率就越大。所以當煤礦發生突水事故時，礦工選擇的最優逃生路徑就是從井下突水點到安全出口的所有的安全、可能路徑中，找出當量長度最小的逃生路徑。假設逃生人員從受災點到安全出口的所有路徑集合為P，P中當量長度最小的路徑即所求的最優逃生路徑由n條巷道組成，Sij為巷道內不同時段下的通行危險系數，建立礦井突水災害逃生路徑選擇的目標優化數學模型，可表示為

(9)

約束條件：

(10)

(11)

(12)

xij={0,1},i=1,2,…n;j=1,2,…,n

(13)

式中：fL(P)為路徑P的當量長度；Lij為巷道(Vi,Vj)的初始當量長度，Lij=0,∞分別為Vi=Vj、Vi與Vj不相鄰，Lij>0即為式(8)所求得的值；式(11)表示xij的取值構成從源節點V1到目的節點Vn的一條可行逃生路徑；式(12)表示逃生路徑中不含回路；xij為決策變量，xij=0,1分別表示巷道(Vi,Vj)不在、在選定的路線P上。

3 改進的Dijkstra算法

由于水位高度和水的流速隨著時間變化，那么在特定時段內用靜態路徑選擇算法所求的最優路徑，無法證明在全時域保持最優；為了動態描述各巷道在不同時間下的危險性情況，根據文獻[21]，將所研究的時間分為K個時段，基于上述數學模型求解出每個時段內各巷道的危險系數，加權到巷道初始當量長度上以此作為逃生路徑的危險性權重，而巷道權值隨不同時段的危險性權重大小變化而變化。傳統的Dijkstra算法不能解決最優路徑的動態選擇問題，所以對該算法進行改進，步驟如下。

節點與節點之間的連線上的數值為巷道權值，如V16與V17之間的17(19)，其中17為當前時段下的巷道權值，括號里的19為更新后的巷道權值圖2 巷道網絡示意圖Fig.2 Roadway network diagram

步驟1初始化。確定突水點為Vi，安全節點為Vj，將突水點Vi到節點Vx的最小累計權值記為Wix,Pix為到節點Vx所經過的路徑節點集合。利用集合Q表示最小累計權值已知的節點，用集合U表示最小累計權值未知的點。初始時，Q只包含突水點Vi，U中包含除點Vi外的其他節點；突水點Vi相鄰節點的權值Wi和路徑Pi可以從巷道網絡圖中得到，不相鄰的節點的權值為∞。在圖2中，設突水點為V8，安全節點為V4，巷道網絡權值時變周期為T，時段編號為Tnum，TTnum為巷道權值更新時刻。初始時刻Q={V8}，U={V1,V2,V3,V4,V5,V6,V7,V9,V10,V11,V12,V13,V14,V15,V16,V17}。節點V9、V15、V16與節點V8相鄰，所以W8,9、W8,15、W8,16可以獲知，其余節點累計權值均為∞。

步驟2更新遍歷節點集和邊集。根據模型的篩選指標，從集合E中刪除無效邊，更新集合V和E。由于圖2中E9,7為不可通行路徑，在集合E中去掉E9,7，集合V={V1,V2,V3,V4,V5,V6,V10,V11,V12,V13,V14,V15,V16,V17}。

步驟3對未知節點集進行遍歷。從U中選出距源點Vi權值最小的節點Vk，將該點從U移到Q中，同時，將該點Vk的權值和對應的路徑分別存放到Wik和Pix中。在步驟1的基礎上，設集合U中節點V15距源點V8具有最小累計權值W8,16，故Q={V8,V16}，U={V1,V2,V3,V4,V5,V6,V10,V11,V12,V13,V14,V15,V17}。

步驟4比較更新鄰近節點權值。設節點Vk到鄰近子節點Vg的權值為dkg，判斷Wik+dkg與Wig的大小，更新i～g的最小權值Wig對應路徑Pig。在步驟3的基礎上，節點V16的鄰近節點包括V17和V8(為突水點排除)，由于一開始節點V17的累積權值為無窮大，故此時應當更新W17。

步驟5判斷巷道權值是否更新。若到節點k的最小權值Wk>TTnum，表明到達節點k的巷道將經歷權值更新。在前4步基礎上，當前，Q={V8,V16,V17}，U={V1,V2,V3,V4,V5,V6,V10,V11,V12,V13,V14,V15}。集合U中V14具有最小累計權值W14，且W14>TTnum，說明在E17,14巷道權值發生更新。

步驟6計算路徑實際累計權值。經過步驟5判斷之后，計算實際權值W14=W17+d17,14=32+25=57。

步驟7循環步驟5、6，直至找到突水點Vi到安全節點Vj的最小累計權值Wij結束，即得到最優逃生路徑。

一般情況下，突水事故點Vi固定，逃生出口一般至少有3個(主井、副井、風井)，因此逃生終點Vj在選取Vj1、Vj2、Vj3的情況下，按照上述算法步驟可分別計算出Vi～Vj的總共至少3條路徑。根據具體的礦井突水事故案例，可對其按照路徑當量長度進行排序，選取最優的逃生路線。

4 實例分析

4.1 基于Python仿真實驗

以W煤礦礦井的1035工作面實際巷道布置圖為基礎，簡化其節點，選出50個關鍵節點形成巷道網絡圖(圖3)進行實驗驗證，其中V48、V49、V50是該礦井的井口，將這3個節點設為安全節點。根據該礦井水文地質資料、采掘工程資料得知V1易發生突水，所以將節點V1設為突水點。在未發生突水之前，根據式(8)計算出每條巷道的初始當量長度，如表6所示(限于篇幅只列舉部分數據)。

圖3 1035工作面巷道網絡圖Fig.3 Roadway network diagram of 1035 working face

假設V1點發生突水時突水流量恒定為5.86 m3/s，根據文獻[22]求得突水從一個節點到另一個節點的漫延時間，再經過無向圖寬度優先搜索算法生成突水漫延路線網絡圖(圖4)。

水流下向漫延時間計算公式為

(14)

式(14)中：Lij為巷道初始當量長度，m；Vij為不同巷道內的水流速度，m/s。

水位上向升漲時間計算公式為

表6 巷道結構信息及初始當量長度計算結果(部分)Table 6 Roadway structure information and initial equivalent length calculation results (part )

(15)

式(15)中：V為根據巷道斷面與淹沒巷道的路徑長度計算得到的巷道空間體積，m3；Q為突水流量(恒定)，m3/s；Q1為單位設施排水能力(這里暫不考慮)。

表7 巷道水位上升情況(部分)Table 7 Water level rise inroadway (part)

圖4 突水漫延路線網絡圖Fig.4 Network diagram of water inrush spread route

表8 不同巷道的危險性(部分)Table 8 Risk of different roadways (part)

根據巷道的當量長度及突水漫延方向構建巷道加權有向網絡圖為：設加權有向圖G=[V,E,B,C,l(p,tspread)]，其中，V為巷道節點的集合，E為節點邊的集合，B為巷道節點的標高信息集，C為(Vi,Vj)兩點之間更新后的當量長度，l(p,tspread)為水漫延路徑集，其中，p為水流經過的節點，tspread為水漫延的時間；使用鄰接矩陣來進行存儲數據和描述。然后根據以上數據信息，在Python平臺上編寫程序，模擬在各段時間內，水流可能漫延的巷道，調取巷道數據信息庫，修改各段巷道的平均水位高度和平均水流速，重新計算每條巷道的危險性系數及巷道的當量長度。最終用改進Dijkstra算法計算出突水點到3個安全節點的最短路徑，然后對這3條路徑按照路徑當量長度從小到大排序，結果如表9所示。

由于只考慮在同一個工作面工作的礦工人員逃生，人數相對來說不是很多，所以從這3條中選取1條最優的路徑作為逃生路徑即可，則從突水點到安全節點的最優逃生路徑為：V1→V2→V6→V10→V11→V26→V34→V35→V37→V41→V44→V47→V48(圖5紅色虛線)，路徑當量長度最短為2 239.58 m；并且本結果也適用于當逃生人員比較多時，同時從三條路徑中進行逃生，可有效避免擁堵，提高井下人員的逃生效率。

表9 改進Dijkstra算法的仿真結果Table 9 Simulation results of improved Dijkstra algorithm

圖5 1035工作面的最優逃生路徑Fig.5 Optimal escape path of 1035 working face

4.2 結果比較分析

傳統Dijkstra算法只適合求解路徑權值不變的情況，而無法根據巷道中水位高度及水流速的變化來選擇逃生路徑，也就是說在本實例中只能求解出特定某T時段下的最優路徑，而無法根據整個突水漫延過程來動態選擇出最優路徑。

假設選取T1(k=1)這個時段下的路徑權值，用傳統Dijkstra算法來分別求解突水點V1到安全點V48、V49、V50的最優逃生路徑結果如表10所示。

表10 傳統Dijkstra算法的仿真結果Table 10 Simulation results of traditional Dijkstra algorithm

從表9、表10可知，表10中從突水點V1到安全點V48、V49、V50的最優逃生路徑的路徑當量長度都比表9中求得的結果小，這是因為在T1(k=1)時段下，突水漫延范圍相對較小，某些巷道的危險系數也較小，在巷道靜態因素不變的情況下，路徑當量長度也較小，所以選擇的最優逃生路徑與表9也不相同。

為了體現該優化算法的動態選擇性，根據表9得到的最優結果，以到達安全節點V48的逃生路徑動態選擇為例，分別選取T12(k=12)、T13(k=13)、T14(k=14)3個時段下運用改進的Dijkstra算法仿真結果進行對比，結果如表11所示。

表11 T12、T13、T14時段下改進Dijkstra算法的前2條最優逃生路徑仿真結果Table 11 Simulation results of the former two optimal escape paths of improved-Dijkstra algorithm at T12,T13,T14

比較表9～表11可知，傳統Dijkstra算法僅僅只能反映當前時間段下的最優路徑，不能根據突水災害實際情況動態選擇出合理的逃生路徑，而表9得到的最優逃生路徑是根據表11不同時段下巷道危險性程度進行動態選擇的，得到的路徑更加安全、可靠，改進的Dijkstra算法恰好彌補了傳統算法的不足。因此，本文模型以及改進的算法比傳統的模型和傳統的Dijkstra算法在礦井突水動態環境中求解最優逃生路徑更具優勢。

5 結論

(1)為了更符合礦井突水災害逃生路徑選擇的實際情況，同時考慮水位高度和水流速對人體穩定性的影響，并據此提出一種危險系數的計算方法，優化傳統最優逃生路徑數學模型。

(2)通過劃分時間段來估算隨時間變化的水位高度與水流速，從而計算各時間段內的預測危險系數，進而得到巷道當量長度，以此動態的選擇最優逃生路徑。

(3)改進后的Dijkstra算法，根據劃分的時間段求得對應的危險性權重，從而更新每條巷道當量長度的鄰接矩陣，突破了傳統的Dijkstra算法不能用于動態路徑搜索的不足，使得在礦井突水最優路徑的選擇中更具有精確性與實用性。通過仿真實驗分別求得了從突水點V1到3個安全地點V48、V49、V50的最優逃生路徑：①V1→V2→V6→V10→V11→V26→V34→V35→V37→V41→V44→V47→V48；②V1→V3→V5→V9→V13→V14→V17→V18→V26→V34→V40→V49；③V1→V2→V6→V10→V11→V26→V34→V35→V36→V39→V42→V45→V46→V50；然后對這3條路徑按照路徑當量長度從小到大排序，選取其中一條最優的路徑，即V1→V2→V6→V10→V11→V26→V34→V35→V37→V41→V44→V47→V48，路徑當量長度為2 239.58 m作為該工作面的突水逃生路徑。

(4)此項研究不僅適用于礦井突水動態最優逃生路徑的選擇，還可以應用于地鐵隧道透水事故等逃生路徑動態選擇中，具有普遍適用的意義。