基于无人机建模的尾矿库溃坝事故推演*

曾 航，周虹延，胡少华，李墨潇，吴 浩，王世杰

(1.武汉理工大学 安全科学与应急管理学院，湖北 武汉 430070；2.华中师范大学 城市与环境科学学院，湖北 武汉 430079)

0 引言

尾矿库是1种人造高势能泥石流危险源，一旦发生溃坝，将造成尾矿库下游人员伤亡、财产损失以及环境污染[1-2]。2010年9月21日，广东信宜银岩尾矿库溃坝事故造成22人死亡[3]；2017年3月12日，黄石大冶铜绿山铜矿西北坝段溃坝事故造成直接经济损失4 518.28万元[4]；2019年1月25日，巴西Brumadinho铁矿尾矿库溃坝事故造成大量人员伤亡，对周边生态造成严重污染[5]。因此，研究尾矿库溃坝水砂演进过程，对预测尾矿库溃坝事故致灾后果及事故应急处置具有重要意义。

国内外学者主要通过模型试验法、理论分析法和数值模拟法对溃坝水砂演进进行研究[6-8]。模型试验法基于相似准则原理，仿照原体实物搭建模型进行试验研究：张力霆等[9]采用小比尺、大模型设计开展尾矿坝溃坝模型试验，通过对溃坝泥石流进行观测及反演分析，运用数值模拟验证试验结果正确性。理论分析法主要对溃坝过程中流速、冲击力、总泄沙量、淹没范围等规律进行研究：Shugan等[10]基于Benney浅水方程建立溃坝水流理论模型，描述含涡流溃坝运动过程，并通过试验验证不同初始水位条件下溃坝水流模型理论预测值;陈生水等[11]考虑渗流场、应力场、溃口及底床变化对溃坝水砂流量过程影响，建立连续降雨条件下尾矿库溃坝失稳模型，模拟尾矿库溃坝溃口变化过程;郑欣等[12]运用FLUENT软件模拟尾矿库溃坝后水砂演进过程，并基于流速、淹没范围等参数估算溃坝生命损失；陈俊[13]运用Flow-3D软件对马家沟尾矿库溃坝水砂在下游河道演进过程，并分析降渠高程和扩渠宽度对下游地区受灾程度影响。

但模型试验法造价成本高，需要巨大人力物力，试验结果缺乏可重复性及普适性；理论分析法与实际尾矿库水砂演进存在一定差别，无法解决复杂三维地形条件下水砂运移问题。随计算机技术发展，数值模拟方法可用于模拟泥沙在水流中运动状态，且计算模型可调试，能够较好地描述尾矿库溃坝水砂演进过程。但大部分数值模拟研究均未考虑尾矿库区地形因素对水砂演进过程影响。

近年来，无人机航测技术以精度高、灵活性强、成本低等优点得到广泛运用，无人机航测技术能快速获取高精度三维空间数据。Salvini[14]利用无人机航测技术对意大利某闭库大理岩矿区工程进行危险性分析；马国超等[15]结合激光扫描和无人机倾斜摄影技术，实现露天采场安全生产数字化监控；杨超等[16]基于无人机摄影测量对尾矿坝边坡表面变形进行整体监测。

本文采用无人机航测技术获取尾矿库三维空间数据，利用尾矿库数字高程模型进行溃坝事故模拟，得到水砂演进过程、各监测点水位高程及尾砂淤积厚度随时间变化规律，分析事故影响范围及影响程度。

1 基于无人机航测的三维地形建模

利用携带高性能相机或视频摄像机等低空数字摄影测量系统的无人驾驶飞行器，对地面进行有规划地重叠式拍摄，获取特定区域影像照片，对影像照片进行数据处理得到三维重建高分辨率地形。无人机航摄三维地形建模主要流程包括影像数据获取、数据预处理、空三加密及地理信息测绘产品生成4个部分，如图1所示。

图1 尾矿库三维地形建模

1.1 尾矿库影像数据获取

鑫荣尾矿库地处十堰市竹山县西沟谷地，为山谷型尾矿库，如图2所示。库区地质结构相对较完整，未见熔岩及断层分布，属区域性稳定地块内部地段，汇水面积1.36 km2，总库容304.2×104 m3，属三等库坝。尾矿库下游共有2处居民聚集点，其中位于上游沟谷的张家院距尾矿坝坝约500 m，有临时厂房分布；位于下游沟谷的银洞村距尾矿坝约1 300 m。

图2 鑫荣尾矿库基本特征

鑫荣尾矿库影像数据获取流程主要包括实地勘测、航迹规划、飞行作业3个部分。经实地勘测，基于下游居民聚集点及重要设施建筑，划定飞行区域为尾矿库库区至下游沟谷区域范围，覆盖面积达1.348 km2。选用哈瓦无人机(MEGA V8III)进行航摄作业，分辨率6 000像素×4 000像素，镜头焦距35 mm，具备高精度成图，支持PPK/RTK及其融合作业模式，智能、高效且环境适应能力强。

飞行作业前需确定摄影高度、航向重叠和旁向重叠。航向重叠指沿航线方向飞行的两相邻航拍照片对所摄地面重叠部分；旁向重叠指两相邻航带航拍照片对所摄地面重叠部分，如图3所示。像片重叠度如式(1)所示：

图3 航向重叠和旁向重叠

(1)

式中：lx、ly为像幅边长，m；Px、Py为航向和旁向重叠影像部分长度，m。

航迹规划需考虑作业区域环境要素如地形、地势对飞行计划影响。当航拍区域地势平坦，航向和旁向重叠率可适当减小，保证影像处理质量；当航拍区域地势崎岖，需对航摄参数进行调整以保证照片质量。尾矿库坝坡陡峭，取航向重叠度80%，旁向重叠度70%；下游沟谷地区及地势较为平坦处，取航向重叠度70%，旁向重叠度60%，经计算设定无人机摄影高度100 m。将整个尾矿库区划分为4个航拍区域，分别为尾矿坝库区(a)、上游沟谷区域(b)、上山主干道区域(c)和下游沟谷区域(d)，各飞行区域航摄参数见表1。

表1 各区域航摄参数

尾矿坝库区无人机航线规划及航拍影像如图4～5所示。

图4 尾矿库航线规划

图5 尾矿库航拍影像

1.2 空三加密

无人机飞行作业完成后，对航拍影像进行畸变差校正和空三加密处理。畸变校正目的是消除系统性因素引起的几何畸变，如可能存在于CCD阵列中的排列误差和因摄影镜头导致的非线性畸变[17-18]。运用数据处理软件中影像畸变差校正模块对原始航摄影像处理，消除影像畸变差和像主点偏移量。空三加密主要根据影像的像点测量坐标和少量控制点大地坐标，求解输出影像定向点大地坐标和影像外方位元素，主要包括相对定向、绝对定向解算和光束法区域网平差[19-20]。鑫荣尾矿库三维空间数据如图6所示。

图6 鑫荣尾矿库三维空间数据

1.3 漫顶溃坝模型

基于数字高程模型(DEM)数据，导入摄影测量三维重建软件，建立尾矿库重建三维模型，转化为stl格式导入三维流体计算软件，得到漫顶溃坝初始模型，如图7所示。

图7 鑫荣尾矿库漫顶溃坝模型

2 尾矿库漫顶溃坝模拟

2.1 尾矿库溃坝数学模型

采用流体计算软件模拟鑫荣尾矿库漫顶溃坝，软件基于重整化群(Re-Normalization Group，RNG)k-ε湍流模型与泥沙冲刷模型，采用优化Tru-VOF(Volume Of Fluid，VOF)函数捕获自由液面及其位置，结合面积/体积分数比(Fractional Area/Volume Ratios，FAVOR)技术划分网格，实现对泥沙沉积、侵蚀和移动的准确计算。湍流模型控制方程中连续方程如式(2)所示：

(2)

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

(3)

式中：ρ为流体密度，kg/m3；P为作用在流体微元压力，N/m2；Ax，Ay，Az分别为x，y，z方向可流动面积分数；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 模拟工况及边界条件

对模拟计算区域进行网格划分，网格划分4 m×4 m×4 m，溃坝模型网格划分如图8所示。

图8 溃坝模型网格划分

划分网格后需考虑边界条件，其中，Xmin为对称(Symmetry)边界条件；Xmax为自由出流(Outflow)边界条件；Ymin，Ymax为自由出流(Outflow)边界条件；Zmin为壁面(Wall)边界条件；Zmax为压力(Pressure)边界条件。为避免溃口出现的随机性，同时基于安全裕度考虑，模拟最危险工况发生，即当溃坝水流全部冲击坝体下游区域，才会造成居民财产最大损失，在坝顶中部位置预制1个宽30 m，深2 m的溃口，引导水砂沿预制溃口流出，保证计算结果准确性。泥沙模型添加尾砂特性见表2，其中临界希尔兹数设置为软件自动计算。

表2 尾砂特性参数

2.3 溃坝模拟结果

溃坝模拟结果表明，尾矿库漫顶溃坝后水砂积蓄的重力势能转化为动能，水砂沿下游方向以一定速度移动。不同时刻挟沙水流演进范围及流速分布如图9所示。漫顶60 s，水砂流到达下游坝脚处，对张家院村居民以及建筑物造成影响，此时溃坝挟沙水流表面流速较大，最快达10 m/s；溃决水流波前速度因下游地形表面粗糙度影响略有减小，由于下乡主干道存在一定地形高差，水流波前速度约8 m/s；漫顶100 s，挟沙水流部分向北冲击张家院村南侧，另一部分向南经由下乡主干道方向演进；漫顶300 s，水流到达下游1 000 m处；漫顶480 s，水流开始对2#沟谷中的银洞村村民聚居地造成影响；溃坝发生780 s，库内水位基本见底，下游2处村庄完全被水砂淹没；漫顶960 s，库内水位干涸，张家院村水流基本退去，但银洞村全部被水砂流淹没。

3 分析与讨论

为监测分析尾矿库漫顶溃坝后水砂运动状态及对下游村庄影响程度，在三维流体计算软件模拟水砂运动路径布置6个监测点，如图2所示。其中1#监测点(3-1)布设于张家院村北侧；2#监测点(3-2)布设于张家院村南侧，即尾矿库坝脚区域；3#监测点(3-3)布设于2#沟谷上游；4#监测点(3-4)布设于2条沟谷交汇处南；5#监测点(3-5)布设于银洞村西侧建筑区域位置；6#检测点(3-6)布设于银洞村东侧房屋聚居点。整理数据可知，点3-3未监测到数据变化，表明2#沟谷上游区域不会受水砂流影响；1、2监测点地理位置相近，4、5、6监测点均处下游沟谷，因此，将1、2监测点和4、5、6监测点监测数据分开整理。5个测点位置水位高程随时间变化及尾砂淤积厚度随时间变化如图10～11所示。

由图10可知，在尾矿库坝脚区域2#监测点最大水位高程达17 m。结合实际地形分析可知，2#监测点靠近坝址，地势地形较平缓，水砂流至2#监测点时，在地势低洼处汇集并迅速填满，水位迅速抬高至12 m。溃坝180 s，2#监测点水位下降，水流除正常下流且向北流动演进；溃坝240 s时，位于沟谷交叉处4#监测点与位于张家院北侧1#监测点几乎同时监测到水位变化，1#监测点水位在630 s内达到最大值8 m，随后立即下降，随水流消去而减小；4#监测点水位受影响后迅速上升至6 m，随后一直稳定不变；同时2#监测点水位抬高至17 m，后随水流演进不断在12 m上下波动；5#、6#监测点位于尾矿坝下游平缓区域，监测点水位总体呈上升趋势，5#监测点最大水位高程8 m，6#监测点最大水位高程10 m。

由图11可知，尾砂淤积主要集中在1#、2#监测点，2#监测点尾砂沉积厚度峰值达27 m，1#监测点尾砂沉积厚度最大达10.5 m，4#、5#、6#监测点尾砂沉积厚度均在2 m左右。监测点尾砂沉积厚度与尾矿坝距离成反比，因为尾砂随水流演进随地形地势沉降，同时随库水位下降，水流流速下降，水流对尾砂携带作用进一步减弱，2#监测点距尾矿坝最近，尾砂淤积厚度在660 s时达到最高，随后呈下降趋势；1#监测点位于尾矿坝上游，且受地形地势影响，尾砂沉积厚度达峰值后处于稳定状态；4#、5#、6#监测点因距尾矿坝较远，尾砂沉积处于较低水平。

图11 尾砂淤积厚度随时间变化

根据5个监测点水位高程与尾砂数据可知，水位高程及尾砂淤积厚度总体变化规律为先迅速增加后趋于稳定值，后缓慢下降。同时监测点尾砂沉积变化相对水位高程变化较慢，与尾砂流速小于水流速度一致。各监测点水位高程数据和尾砂淤积数据见表3。1#、2#监测点地处上游沟谷，水位高程平均可达5.6 m，尾砂淤积高度平均可达9.6 m，对该范围处临时厂房及居民聚集点采取溃坝预防措施；4#、5#、6#监测点地处下游沟谷，水位高程平均可达2.7 m，尾砂淤积高度平均可达0.96 m，采取预防措施可有效减轻下游居民聚集点受灾程度。

表3 各监测点水位高程和尾砂淤积厚度

4 结论

1)无人机倾斜摄影测量技术可快速获取高精度尾矿库三维空间数据及实景三维模型，实景三维模型可真实还原实地原貌，精确复现尾矿库溃坝灾变演化过程。

2)由溃坝模拟结果可知，溃坝水砂流速最大可达10 m/s，漫顶300 s时，水砂到达下游1 000 m处；位于上游沟谷的张家院南侧2#监测点水位高程最高处达17 m，尾砂淤积厚度最高达27 m。说明该尾矿库一旦溃坝，将会直接淹没上游沟谷临时厂房和下游沟谷房屋村落，并淹没和侵蚀乡道公路。各监测点水位高程及尾砂淤积厚度总体变化规律为先迅速增加后趋于某稳定值，后缓慢下降。因流速原因，尾砂淤积厚度变化相对于水位高程变化存在滞后现象。

3)研究结果基于无人机建模的尾矿库三维数字高程模型，模拟水砂淹没范围及影响区域更准确，对尾矿库溃坝事故应急措施更有针对性，同时结合水流演进和尾砂淤积范围可为尾矿库建设拦挡工程提供指导建议。