基于無人機建模的尾礦庫潰壩事故推演*

曾 航，周虹延，胡少華，李墨瀟，吳 浩，王世杰

(1.武漢理工大學 安全科學與應急管理學院，湖北 武漢 430070；2.華中師范大學 城市與環(huán)境科學學院，湖北 武漢 430079)

0 引言

尾礦庫是1種人造高勢能泥石流危險源，一旦發(fā)生潰壩，將造成尾礦庫下游人員傷亡、財產損失以及環(huán)境污染[1-2]。2010年9月21日，廣東信宜銀巖尾礦庫潰壩事故造成22人死亡[3]；2017年3月12日，黃石大冶銅綠山銅礦西北壩段潰壩事故造成直接經濟損失4 518.28萬元[4]；2019年1月25日，巴西Brumadinho鐵礦尾礦庫潰壩事故造成大量人員傷亡，對周邊生態(tài)造成嚴重污染[5]。因此，研究尾礦庫潰壩水砂演進過程，對預測尾礦庫潰壩事故致災后果及事故應急處置具有重要意義。

國內外學者主要通過模型試驗法、理論分析法和數(shù)值模擬法對潰壩水砂演進進行研究[6-8]。模型試驗法基于相似準則原理，仿照原體實物搭建模型進行試驗研究：張力霆等[9]采用小比尺、大模型設計開展尾礦壩潰壩模型試驗，通過對潰壩泥石流進行觀測及反演分析，運用數(shù)值模擬驗證試驗結果正確性。理論分析法主要對潰壩過程中流速、沖擊力、總泄沙量、淹沒范圍等規(guī)律進行研究：Shugan等[10]基于Benney淺水方程建立潰壩水流理論模型，描述含渦流潰壩運動過程，并通過試驗驗證不同初始水位條件下潰壩水流模型理論預測值;陳生水等[11]考慮滲流場、應力場、潰口及底床變化對潰壩水砂流量過程影響，建立連續(xù)降雨條件下尾礦庫潰壩失穩(wěn)模型，模擬尾礦庫潰壩潰口變化過程;鄭欣等[12]運用FLUENT軟件模擬尾礦庫潰壩后水砂演進過程，并基于流速、淹沒范圍等參數(shù)估算潰壩生命損失；陳俊[13]運用Flow-3D軟件對馬家溝尾礦庫潰壩水砂在下游河道演進過程，并分析降渠高程和擴渠寬度對下游地區(qū)受災程度影響。

但模型試驗法造價成本高，需要巨大人力物力，試驗結果缺乏可重復性及普適性；理論分析法與實際尾礦庫水砂演進存在一定差別，無法解決復雜三維地形條件下水砂運移問題。隨計算機技術發(fā)展，數(shù)值模擬方法可用于模擬泥沙在水流中運動狀態(tài)，且計算模型可調試，能夠較好地描述尾礦庫潰壩水砂演進過程。但大部分數(shù)值模擬研究均未考慮尾礦庫區(qū)地形因素對水砂演進過程影響。

近年來，無人機航測技術以精度高、靈活性強、成本低等優(yōu)點得到廣泛運用，無人機航測技術能快速獲取高精度三維空間數(shù)據(jù)。Salvini[14]利用無人機航測技術對意大利某閉庫大理巖礦區(qū)工程進行危險性分析；馬國超等[15]結合激光掃描和無人機傾斜攝影技術，實現(xiàn)露天采場安全生產數(shù)字化監(jiān)控；楊超等[16]基于無人機攝影測量對尾礦壩邊坡表面變形進行整體監(jiān)測。

本文采用無人機航測技術獲取尾礦庫三維空間數(shù)據(jù)，利用尾礦庫數(shù)字高程模型進行潰壩事故模擬，得到水砂演進過程、各監(jiān)測點水位高程及尾砂淤積厚度隨時間變化規(guī)律，分析事故影響范圍及影響程度。

1 基于無人機航測的三維地形建模

利用攜帶高性能相機或視頻攝像機等低空數(shù)字攝影測量系統(tǒng)的無人駕駛飛行器，對地面進行有規(guī)劃地重疊式拍攝，獲取特定區(qū)域影像照片，對影像照片進行數(shù)據(jù)處理得到三維重建高分辨率地形。無人機航攝三維地形建模主要流程包括影像數(shù)據(jù)獲取、數(shù)據(jù)預處理、空三加密及地理信息測繪產品生成4個部分，如圖1所示。

圖1 尾礦庫三維地形建模

1.1 尾礦庫影像數(shù)據(jù)獲取

鑫榮尾礦庫地處十堰市竹山縣西溝谷地，為山谷型尾礦庫，如圖2所示。庫區(qū)地質結構相對較完整，未見熔巖及斷層分布，屬區(qū)域性穩(wěn)定地塊內部地段，匯水面積1.36 km2，總庫容304.2×104 m3，屬三等庫壩。尾礦庫下游共有2處居民聚集點，其中位于上游溝谷的張家院距尾礦壩壩約500 m，有臨時廠房分布；位于下游溝谷的銀洞村距尾礦壩約1 300 m。

圖2 鑫榮尾礦庫基本特征

鑫榮尾礦庫影像數(shù)據(jù)獲取流程主要包括實地勘測、航跡規(guī)劃、飛行作業(yè)3個部分。經實地勘測，基于下游居民聚集點及重要設施建筑，劃定飛行區(qū)域為尾礦庫庫區(qū)至下游溝谷區(qū)域范圍，覆蓋面積達1.348 km2。選用哈瓦無人機(MEGA V8III)進行航攝作業(yè)，分辨率6 000像素×4 000像素，鏡頭焦距35 mm，具備高精度成圖，支持PPK/RTK及其融合作業(yè)模式，智能、高效且環(huán)境適應能力強。

飛行作業(yè)前需確定攝影高度、航向重疊和旁向重疊。航向重疊指沿航線方向飛行的兩相鄰航拍照片對所攝地面重疊部分；旁向重疊指兩相鄰航帶航拍照片對所攝地面重疊部分，如圖3所示。像片重疊度如式(1)所示：

圖3 航向重疊和旁向重疊

(1)

式中：lx、ly為像幅邊長，m；Px、Py為航向和旁向重疊影像部分長度，m。

航跡規(guī)劃需考慮作業(yè)區(qū)域環(huán)境要素如地形、地勢對飛行計劃影響。當航拍區(qū)域地勢平坦，航向和旁向重疊率可適當減小，保證影像處理質量；當航拍區(qū)域地勢崎嶇，需對航攝參數(shù)進行調整以保證照片質量。尾礦庫壩坡陡峭，取航向重疊度80%，旁向重疊度70%；下游溝谷地區(qū)及地勢較為平坦處，取航向重疊度70%，旁向重疊度60%，經計算設定無人機攝影高度100 m。將整個尾礦庫區(qū)劃分為4個航拍區(qū)域，分別為尾礦壩庫區(qū)(a)、上游溝谷區(qū)域(b)、上山主干道區(qū)域(c)和下游溝谷區(qū)域(d)，各飛行區(qū)域航攝參數(shù)見表1。

表1 各區(qū)域航攝參數(shù)

尾礦壩庫區(qū)無人機航線規(guī)劃及航拍影像如圖4～5所示。

圖4 尾礦庫航線規(guī)劃

圖5 尾礦庫航拍影像

1.2 空三加密

無人機飛行作業(yè)完成后，對航拍影像進行畸變差校正和空三加密處理。畸變校正目的是消除系統(tǒng)性因素引起的幾何畸變，如可能存在于CCD陣列中的排列誤差和因攝影鏡頭導致的非線性畸變[17-18]。運用數(shù)據(jù)處理軟件中影像畸變差校正模塊對原始航攝影像處理，消除影像畸變差和像主點偏移量。空三加密主要根據(jù)影像的像點測量坐標和少量控制點大地坐標，求解輸出影像定向點大地坐標和影像外方位元素，主要包括相對定向、絕對定向解算和光束法區(qū)域網平差[19-20]。鑫榮尾礦庫三維空間數(shù)據(jù)如圖6所示。

圖6 鑫榮尾礦庫三維空間數(shù)據(jù)

1.3 漫頂潰壩模型

基于數(shù)字高程模型(DEM)數(shù)據(jù)，導入攝影測量三維重建軟件，建立尾礦庫重建三維模型，轉化為stl格式導入三維流體計算軟件，得到漫頂潰壩初始模型，如圖7所示。

圖7 鑫榮尾礦庫漫頂潰壩模型

2 尾礦庫漫頂潰壩模擬

2.1 尾礦庫潰壩數(shù)學模型

采用流體計算軟件模擬鑫榮尾礦庫漫頂潰壩，軟件基于重整化群(Re-Normalization Group，RNG)k-ε湍流模型與泥沙沖刷模型，采用優(yōu)化Tru-VOF(Volume Of Fluid，VOF)函數(shù)捕獲自由液面及其位置，結合面積/體積分數(shù)比(Fractional Area/Volume Ratios，F(xiàn)AVOR)技術劃分網格，實現(xiàn)對泥沙沉積、侵蝕和移動的準確計算。湍流模型控制方程中連續(xù)方程如式(2)所示：

(2)

動量方程如式(3)所示：

(3)

式中：ρ為流體密度，kg/m3；P為作用在流體微元壓力，N/m2；Ax，Ay，Az分別為x，y，z方向可流動面積分數(shù)；u，v，w分別為x，y，z方向速度分量，m/s；Gx，Gy，Gz分別為重力加速度在x，y，z坐標軸方向上分量，m/s2；fx，fy，fz分別為黏滯力加速度在x，y，z軸方向上的分量，m/s2。

湍流控制方程如式(4)所示：

(4)

2.2 模擬工況及邊界條件

對模擬計算區(qū)域進行網格劃分，網格劃分4 m×4 m×4 m，潰壩模型網格劃分如圖8所示。

圖8 潰壩模型網格劃分

劃分網格后需考慮邊界條件，其中，Xmin為對稱(Symmetry)邊界條件；Xmax為自由出流(Outflow)邊界條件；Ymin，Ymax為自由出流(Outflow)邊界條件；Zmin為壁面(Wall)邊界條件；Zmax為壓力(Pressure)邊界條件。為避免潰口出現(xiàn)的隨機性，同時基于安全裕度考慮，模擬最危險工況發(fā)生，即當潰壩水流全部沖擊壩體下游區(qū)域，才會造成居民財產最大損失，在壩頂中部位置預制1個寬30 m，深2 m的潰口，引導水砂沿預制潰口流出，保證計算結果準確性。泥沙模型添加尾砂特性見表2，其中臨界希爾茲數(shù)設置為軟件自動計算。

表2 尾砂特性參數(shù)

2.3 潰壩模擬結果

潰壩模擬結果表明，尾礦庫漫頂潰壩后水砂積蓄的重力勢能轉化為動能，水砂沿下游方向以一定速度移動。不同時刻挾沙水流演進范圍及流速分布如圖9所示。漫頂60 s，水砂流到達下游壩腳處，對張家院村居民以及建筑物造成影響，此時潰壩挾沙水流表面流速較大，最快達10 m/s；潰決水流波前速度因下游地形表面粗糙度影響略有減小，由于下鄉(xiāng)主干道存在一定地形高差，水流波前速度約8 m/s；漫頂100 s，挾沙水流部分向北沖擊張家院村南側，另一部分向南經由下鄉(xiāng)主干道方向演進；漫頂300 s，水流到達下游1 000 m處；漫頂480 s，水流開始對2#溝谷中的銀洞村村民聚居地造成影響；潰壩發(fā)生780 s，庫內水位基本見底，下游2處村莊完全被水砂淹沒；漫頂960 s，庫內水位干涸，張家院村水流基本退去，但銀洞村全部被水砂流淹沒。

3 分析與討論

為監(jiān)測分析尾礦庫漫頂潰壩后水砂運動狀態(tài)及對下游村莊影響程度，在三維流體計算軟件模擬水砂運動路徑布置6個監(jiān)測點，如圖2所示。其中1#監(jiān)測點(3-1)布設于張家院村北側；2#監(jiān)測點(3-2)布設于張家院村南側，即尾礦庫壩腳區(qū)域；3#監(jiān)測點(3-3)布設于2#溝谷上游；4#監(jiān)測點(3-4)布設于2條溝谷交匯處南；5#監(jiān)測點(3-5)布設于銀洞村西側建筑區(qū)域位置；6#檢測點(3-6)布設于銀洞村東側房屋聚居點。整理數(shù)據(jù)可知，點3-3未監(jiān)測到數(shù)據(jù)變化，表明2#溝谷上游區(qū)域不會受水砂流影響；1、2監(jiān)測點地理位置相近，4、5、6監(jiān)測點均處下游溝谷，因此，將1、2監(jiān)測點和4、5、6監(jiān)測點監(jiān)測數(shù)據(jù)分開整理。5個測點位置水位高程隨時間變化及尾砂淤積厚度隨時間變化如圖10～11所示。

由圖10可知，在尾礦庫壩腳區(qū)域2#監(jiān)測點最大水位高程達17 m。結合實際地形分析可知，2#監(jiān)測點靠近壩址，地勢地形較平緩，水砂流至2#監(jiān)測點時，在地勢低洼處匯集并迅速填滿，水位迅速抬高至12 m。潰壩180 s，2#監(jiān)測點水位下降，水流除正常下流且向北流動演進；潰壩240 s時，位于溝谷交叉處4#監(jiān)測點與位于張家院北側1#監(jiān)測點幾乎同時監(jiān)測到水位變化，1#監(jiān)測點水位在630 s內達到最大值8 m，隨后立即下降，隨水流消去而減小；4#監(jiān)測點水位受影響后迅速上升至6 m，隨后一直穩(wěn)定不變；同時2#監(jiān)測點水位抬高至17 m，后隨水流演進不斷在12 m上下波動；5#、6#監(jiān)測點位于尾礦壩下游平緩區(qū)域，監(jiān)測點水位總體呈上升趨勢，5#監(jiān)測點最大水位高程8 m，6#監(jiān)測點最大水位高程10 m。

由圖11可知，尾砂淤積主要集中在1#、2#監(jiān)測點，2#監(jiān)測點尾砂沉積厚度峰值達27 m，1#監(jiān)測點尾砂沉積厚度最大達10.5 m，4#、5#、6#監(jiān)測點尾砂沉積厚度均在2 m左右。監(jiān)測點尾砂沉積厚度與尾礦壩距離成反比，因為尾砂隨水流演進隨地形地勢沉降，同時隨庫水位下降，水流流速下降，水流對尾砂攜帶作用進一步減弱，2#監(jiān)測點距尾礦壩最近，尾砂淤積厚度在660 s時達到最高，隨后呈下降趨勢；1#監(jiān)測點位于尾礦壩上游，且受地形地勢影響，尾砂沉積厚度達峰值后處于穩(wěn)定狀態(tài)；4#、5#、6#監(jiān)測點因距尾礦壩較遠，尾砂沉積處于較低水平。

圖11 尾砂淤積厚度隨時間變化

根據(jù)5個監(jiān)測點水位高程與尾砂數(shù)據(jù)可知，水位高程及尾砂淤積厚度總體變化規(guī)律為先迅速增加后趨于穩(wěn)定值，后緩慢下降。同時監(jiān)測點尾砂沉積變化相對水位高程變化較慢，與尾砂流速小于水流速度一致。各監(jiān)測點水位高程數(shù)據(jù)和尾砂淤積數(shù)據(jù)見表3。1#、2#監(jiān)測點地處上游溝谷，水位高程平均可達5.6 m，尾砂淤積高度平均可達9.6 m，對該范圍處臨時廠房及居民聚集點采取潰壩預防措施；4#、5#、6#監(jiān)測點地處下游溝谷，水位高程平均可達2.7 m，尾砂淤積高度平均可達0.96 m，采取預防措施可有效減輕下游居民聚集點受災程度。

表3 各監(jiān)測點水位高程和尾砂淤積厚度

4 結論

1)無人機傾斜攝影測量技術可快速獲取高精度尾礦庫三維空間數(shù)據(jù)及實景三維模型，實景三維模型可真實還原實地原貌，精確復現(xiàn)尾礦庫潰壩災變演化過程。

2)由潰壩模擬結果可知，潰壩水砂流速最大可達10 m/s，漫頂300 s時，水砂到達下游1 000 m處；位于上游溝谷的張家院南側2#監(jiān)測點水位高程最高處達17 m，尾砂淤積厚度最高達27 m。說明該尾礦庫一旦潰壩，將會直接淹沒上游溝谷臨時廠房和下游溝谷房屋村落，并淹沒和侵蝕鄉(xiāng)道公路。各監(jiān)測點水位高程及尾砂淤積厚度總體變化規(guī)律為先迅速增加后趨于某穩(wěn)定值，后緩慢下降。因流速原因，尾砂淤積厚度變化相對于水位高程變化存在滯后現(xiàn)象。

3)研究結果基于無人機建模的尾礦庫三維數(shù)字高程模型，模擬水砂淹沒范圍及影響區(qū)域更準確，對尾礦庫潰壩事故應急措施更有針對性，同時結合水流演進和尾砂淤積范圍可為尾礦庫建設攔擋工程提供指導建議。