摘要尾矿库上游式筑坝法：

基于无人机航测生成的高精度数字高程模型(DEM)数据尾矿库上游式筑坝法，采用深度积分的Massflow数值分析方法，开展了新疆某山区尾矿库溃坝数值模拟研究。结果显示了该尾矿库溃坝后矿砂运动过程，分析数值模拟结果得出：(1)溃坝后矿砂最终堆积呈现两部分，一部分在尾矿库库区下游，平均厚度约13 m,另一部分则堆积于环保库和下游河道附近，平均厚度达25 m,其中堆积厚度最深处在初期坝坝底正对的河流处，约29 m。(2)尾矿库溃坝整个运动过程可分为失稳启动、加速下滑以及减速堆积三个阶段，该尾矿库溃坝从启动到最终静止整个过程约3 000 s,最大运动速度介于21～29 m/s之间，其中失稳启动和加速下滑阶段均在短时间内迅速完成，减速堆积阶段时间占整个过程的93.33%,表明尾矿库溃坝初期矿砂流速较快，具有较强的破坏性。研究成果可用于分析尾矿库矿砂运动路径、淹没区域及沿程堆积情况等，对降低尾矿库溃坝风险和应急响应有重要指导意义。

关键词：

尾矿库溃坝; Massflow; 数值模拟; 摄影测量; DEM;

作者简介：

林小莉(1995—),女，硕士研究生，主要从事遥感、GIS、地质灾害研究。E-mail:linxiaoli18@mails.ucas.ac.cn;*张清凌

(1974—),男，研究员，博士，主要从事遥感、GIS,地质灾害研究。E-mail:zhangqling@mail.sysu.edu.cn;

基金：科技部国家重点研发计划(YFB0504204);新疆维吾尔自治区重点实验室开放课题(2018D04027);

引用：

林小莉，刘英，张清凌，等． 基于无人机航空摄影的尾矿库溃坝数值模拟［J］. 水利水电技术( 中英文) ，2021，52( 5) : 71-81.

LIN Xiaoli，LIU Ying，ZHANG Qingling，et al． UAV aerial photography-based numerical simulation on failure of tailings reservoir dam［J］.Water Ｒesources and Hydropower Engineering，2021，52( 5) : 71-81.

尾矿库是由筑坝拦截谷口或围地构成，主要用于贮存金属或非金属矿山矿石选别后排出的尾矿和其他工业废渣，是矿山正常生产的必要设施。但尾矿库也是矿山的重大危险源之一，一旦发生溃坝将形成具有高势能的人造泥石流，将对其下游居民和设施安全造成严重的威胁，容易导致重大人员伤亡和财产损失。目前，我国共有尾矿库近8 000座，总量居世界第一 ,普遍具有“数量多且规模小、采用上游式筑坝工艺、安全度低、选址很难避免居民区、尾矿库技术力量薄弱 ”等特点。约90%的尾矿库采用上游式筑坝法，这种方法虽然工艺简单、基建投资少，但稳定性差，也最容易发生溃坝事故。例如山西省临汾市襄汾县新塔矿业有限公司980沟尾矿库发生了277人死亡、4人失踪的“9.8”溃坝事故；2019年巴西布鲁马迪市发生了165人死亡、160人失踪的“1.25”溃坝事故 。

尾矿库溃坝属于典型的灾害应急事件，具有尾矿砂流速快、突发性强等特点。应急管理部门在此类事件中，难以及时有效地实施应急处置措施。通过极端情况模拟，预测尾矿库溃坝模拟淹没范围及矿砂流动过程，可以事先规划居民和设施区位置，避开风险区，事后可以依据模拟结果并结合实际情况快速采取响应措施。目前针对尾矿库溃坝尾矿砂运动特征的主要研究方法有经验公式法、模型试验法和数值模拟法 。经验公式法主要是在水库溃坝后泥石流的演进计算公式的基础上进行修正和完善 ;模型试验法基于相似原理，利用地形、尾矿库设计资料等信息制作实验模型来模拟溃坝运动状况 ,例如敬小非等 利用物理模型试验方法模拟了尾矿库堆积过程，并且在此基础上针对不同洪水工况进行模拟，但该方法运算效率低且局限性强，难以开展多条件和多情景下的模拟。数值模拟法 则是将计算机与数理模型进行结合，例如李火坤等 利用Flow 3D软件模拟某尾矿库逐渐溃坝时溃口随时间的变化过程以及溃坝后尾矿砂的淹没范围，其特点是模拟过程不需要专门的试验仪器，且可以通过反复调试，开展多条件和多情景的模拟。但目前大多数研究都对尾矿库模型进行了不同程度简化，并没有考虑实际地形对尾砂运动的影响。

本文在现场勘察的基础上，利用无人机摄影测量三维重建的高分地形作为基础数据，采用深度积分的Massflow

数值分析方法，对新疆山区某尾矿库溃坝进行模拟分析，得出溃坝后不同时间内尾矿砂体的堆积范围和流速变化，可为山区尾矿库溃坝正演提供依据，进而为应急管理部门快速进行应急响应提出具有参考价值的危险范围，旨在减少溃坝造成的人员伤亡及财产损失，以及为该流域的生态环境安全提供数据支撑。

本文研究区为新疆某山区尾矿库，该尾矿库呈南北向条带分布，矿区距县城62 km, 地处西天山博罗科努山系科古琴山脊南坡的河间地块式山岭地带，地形切割角强烈，沟谷发育。该地区地貌属剥蚀中山地形，总体地势东北高西南低，海拔高程1 260～1 700 m, 相对高差350～420 m, 地形坡度大，达30°～35°,当地侵蚀基准面标高1 261 m,大部分矿体位于侵蚀基准面以上。

该尾矿库占地面积约17万m ,坝高45 m, 坝长218 m, 初期库容90万m ,最终库容360万m ,服务期20年，每年尾矿排放量约32万t, 采用上游式筑坝工艺，按照其库容分类属于四等库，按照其地理位置分类属于典型山谷型尾矿库(见图1)。

2.1 技术流程设计

技术流程主要分为三个部分(见图2):第一部分，通过无人机航摄并处理得到数字正射影像(DOM)和数字地形模型(DEM)这两个数据。第二部分，开展基于Massflow模型的溃坝模拟，即在搭建模型时，首先通过遥感影像、高精度DEM、现场调查以及尾矿库初始设计标准确定模拟尾矿库溃坝所需的初始物源量以及设计网格精度；然后主要通过现场试验以及历史事件反演确定模拟所需的基底摩擦模型以及相关系数；最后，进行模拟得到尾矿砂流深、流速数据。第三部分，利用ArcGIS软件，以DOM影像作为背景，实现尾矿砂流深、流速的可视化分析。

图1 尾矿库地理位置

图2 技术流程

2.2 无人机航拍测量数据获取

无人机常常以高时效、高分辨率和高机动性能为人们应对突发事件提供新的途径，其最突出的优点在于它可以用最短时间，提供最大限度接近灾情现场的地形数据，为监测突发事件提供重要依据。

尾矿库及其下游区域地形数据由QC-2型无人机，SONY RX1R M2数字航空摄影仪飞行获得。该地区地形高低起伏，周围地形比较复杂，考虑地物遮挡和变形关系的影响，在航测时为保证影像的处理质量，拍摄照片航向重叠率设为85%,旁向重叠率设为65%,同时为了保证照片的质量，采用交叉飞行的方式设计航线，共采集航摄影像1 625张，影像坐标系采用WGS-84。航摄参数如表1所列。为提高摄影测量精度，在飞行区布设了一些特征比较明显的地物点，通过地面控制点校正，获得准确三维影像。

表1 航空摄影基本参数

利用摄影三维建模软件PhotoScan, 对研究区的无人机航拍照片进行分析处理，并生成研究区的DEM数据(分辨率为0.1 m)。

根据无人机航摄所得数字正射影像和高精度数字高程模型对该矿区实施三维场景搭建，这种方法不仅能够获得研究区的高精度高程信息，准确地反映研究区地形地貌情况，也能够自动、快速且逼真地再现矿区的人造地物和自然地物，如图3所示为该矿区三维场景以及局部三维视图。

3.1 模型介绍

Massflow从流体力学的Navier-Stokes方程出发，基于深度积分的连续介质力学方程，利用改进的MacCormack-TVD有限差分方法，并结合MPICH和OpenMP并行处理技术，能够考虑复杂地形地貌和沟床侵蚀等特征，具有二阶精度和自适应求解功能 。基于Navier-Stokes方程的质量和动量守恒基本物理方程如下。

质量守恒方程

动量守恒方程

式中，p为流体密度；h为流体厚度；g为重力加速度；β为动量分布系数；ξ为几何校正系数；zb为基底边界；ka/p为侧向压力系数；E为侵蚀率；u,v分别为x,y方向的速度矢量；(τ ) ,(τ ) 分别为基底摩擦阻力模型。

3.2 模型搭建

3.2.1 物源数据

通过无人机航测、GPS校正得到尾矿库的高精度地形数据(DEM),分辨率为0.1 m, 采用遥感分类识别结果、现场调查以及尾矿库初始设计标准，确定尾矿库以及下泄尾矿砂范围，模拟区域长约1 353 m, 宽约2 469 m。由于本文考虑尾矿库溃坝后危害的极限化，以尾矿库坝体全部溃决为极端情景。此情景下，尾矿砂一次性全部泻出，所以启动物源方量以及启动物源的高度分布数据可以通过现有高精度DEM数据减去库区建成以后还未堆积矿砂时的DEM,表示为

式中，H1为现有DEM地形高程；H2为尾矿库建成后还未堆积尾矿砂时的DEM地形高程

由于库区建成时DEM数据的精度与目前所拥有DEM的精度不一致，通过ArcGIS的重采样方式将网格分辨率统一设为1 m。

3.2.2 模拟区域网格划分

考虑到Massflow模型采用网格式计算方式，计算网格的大小直接影响模拟结果的精度。理论上网格划分越精细，模拟结果越接近真实值，但是随着计算网格的增加，所需要的计算时间也会呈几何级数的增加。针对以上问题，本文结合该尾矿库溃坝模拟区域的实际情况、计算精度要求以及计算机计算能力，确定溃坝模拟区域选用1 m×1 m的均匀格网(DEM精度为1 m),共计334 056个计算网格，基于Courant准则的变时刻步进行计算。

3.3 基底摩擦模型选取

基于Massflow模型进行溃坝模拟时，需要确定基底摩擦模型，譬如Coulomb模型、Manning模型 、Voellmy模型 ,不同的基底摩擦模型所需的参数不同，尽量选择所需参数较少且易于获取的模型。本研究通过对已有文献的分析总结以及考虑研究区实际情况，在保证参数获取准确的前提下，选择适用于滑坡、岩崩等泥石流灾害，并且基于实验模型的Coulomb摩擦模型，其表达式如下

式中，τb为基底摩阻力；ρ为物体密度；c为黏聚力；φ为内摩擦角；λ为超孔隙水压力系数。

图3 矿区三维地形

图4 研究区尾矿砂颗粒级配曲线

3.4 基底模型参数选取

计算参数的选取对数值模拟的准确性有重要影响。本研究，选用的基底摩擦模型Coulomb模型的相关参数包括矿砂密度、黏聚力、内摩擦角、超孔隙水压力系数。

常规情况下，参数的确定可以根据现场调查、物理模拟试验、历史事件反演以及专家经验等方法 获取。根据调查分析和野外实验结果综合取值，该尾矿库矿砂的平均密度为1 490 kg/m 。

在极端条件下尾矿砂处于极度饱和的状态，黏聚力接近于0。该尾矿库的尾矿砂颗粒粒径主要集中在0.01～0.1mm, 属于极细粒尾砂(见图4),根据相关参考文献记录 ,极细粒尾砂的内摩擦角大致在20°～32°之间，同时类比相关历史事件 ,尾矿砂内摩擦角大致在20°～32°之间(见表2),本研究考虑阻力下限情况，故内摩擦角设置为20°。对于超孔隙水压力系数，由于尾矿库溃坝属于突发型地质灾害，具有很强的流态运动特征，在实际运动过程中，土体局部孔隙水压力是变化的，并且会对岩土体产生刮擦，但复杂的力学计算方程将使数值模拟难以进行，因此本文将整个运动过程中的孔隙水压力系数简化为定值，且不考虑溃坝滑坡运动过程中的砂流密度演化以及沟道侵蚀作用，并忽略动量分布系数(即β=1),根据已有参考文献记录 ,超孔隙水压力系数一般范围在0.60～0.80,本文综合考虑研究区实际情况，设为0.60(见表3)。

表2 相关案例取值

表3 尾矿库溃坝模拟参数

4.1 尾矿砂堆积情况分析

Massflow模型能够通过命令输出DEM格式的不同时刻堆积厚度文件，图5为不同时刻下尾矿砂的堆积情况图，为明细尾矿库溃坝模拟过程中各个区域的矿砂变化情况，分别统计了4个有代表性断面(见图5),在0 s(溃坝开始时间)和3000 s(溃坝结束时间)时尾矿砂堆积厚度。从堆积厚度图(见图6)和尾矿砂流深断面图(见图7)可以看出：

(1)整个过程中，尾矿库中的砂体会在溃坝发生后的极短时间内迅速冲向下游，并沿着库区主体朝向溃决漫流，最终砂体集中堆积于库区中下部、环保库周围、河道及其附近区域。

图5 尾矿砂堆积剖面

(2)溃坝后尾矿砂会依次在不同位置进行快速堆积，10 s左右砂流到达环保库；20 s时淹没环保库，辐射附近居民点并且顺势流向下方河道进行堆积；50～ 2 500 s顺着河流方向堆积，并向两侧扩张，直至3 000 s过程结束。

(3)由于库区内尾矿砂是沿着原有地形进行分布，在库区下部地形较低，所以当尾矿库溃坝时，仍然会有一部分尾矿砂堆积在库区内；加之该尾矿库夹在两山的山谷中，所以溃坝后尾矿砂将受到河道旁山体的阻挡，顺势堆积于河道区域。

(4)尾矿砂在剖面B-B′的50～200 m内堆积厚度变化最大，最大厚度约为20 m; 在剖面C-C′的65～150 m内堆积厚度变化最大，最大厚度约22 m; 在剖面D-D′的200～750 m范围内堆积厚度变化最大，最大厚度约29 m。

(5)最终溃坝后尾矿砂体堆积呈两部分：一部分仍留在尾矿库中，堆积长度约863.5 m, 最大宽度约300 m, 平均厚度约13 m; 另一部分堆积在尾矿库下游区域，辐射面积达234 971.60 m ,平均厚度达25 m, 其中在环保库堆积厚度约21 m, 居民点堆积厚度约11 m, 堆积在河道中的尾矿砂长约1 100 m, 平均向河道左侧扩张85 m, 堆积厚度最大的位置处于尾矿库初期坝坝底正对的河流处，最大厚度约29 m。

4.2 尾矿砂运动速度分析

基于Massflow模型计算的尾矿砂不同时刻下X、Y方向的运动速度，利用ArcGIS栅格计算器算出尾矿砂在溃坝后的合速度如下

式中，v为合速度；x为东西方向的速度；y为南北方向的速度。

图8为不同时刻尾矿砂的运动速度，可以看出：

(1)从溃坝开始，尾矿砂在10 s到达环保库上方，运动速度可达到28 m/s左右，并且这种高速将在20 s时一直持续到尾矿砂到达下方河道。从50 s左右开始，除去尾矿库坝顶至初期坝底，其余部分尾矿砂的运动速度开始出现减少的趋势，直到500 s之前，除去初期坝部分，尾矿砂的最大运动速度约为8 m/s, 在尾矿砂堆积的边缘区域。

图6 尾矿库溃坝关键时间点堆积范围

(2)在整个溃坝过程中，尾矿库坝顶至初期坝底这一部分的运动速度都较高，猜测其主要原因可能是因为地形起伏明显，高程变化较大。从10 s开始，直到100 s左右结束，该区域的最大运动速度虽然都维持在26～29 m/s, 但是覆盖的范围在逐渐减少，主要是因为尾矿砂所聚集的能量在受到地面摩擦的影响下，开始逐渐消散。从500 s开始，尾矿坝部分的最大运动速度以及高值速度分布区域开始减小，此区域的运动速度维持在13～23 m/s, 并且在环保库附近以及下方河道处开始以低速度进行堆积。

(3)整个溃坝过程根据速度变化情况，同滑坡一样可以分为失稳启动、加速下滑以及减速堆积三部分：0～10 s为失稳启动阶段，当开始溃坝时，尾矿砂拥有一定动能而向下滑动，速度也呈现逐渐增加趋势，其中溃坝龙头速度最大，超过25 m/s; 10～20 s为加速下滑阶段，整个尾矿砂运动速度在此阶段不断增加，经过附近居民点的速度也超过了15 m/s, 这种速度足以摧毁房屋，从而造成人员和财产的损失；20～3 000 s为减速堆积阶段，在50 s时尾矿砂到达河流时被对面山体所阻拦，就顺势堆积在河道以及附近区域，平均速度不超过10 m/s。

4.3 不同点位尾矿砂厚度和速度分析

为了更清楚地了解不同点位在尾矿库溃坝所后面临的情况，本文从居民财产安全以及生态安全的角度，选取包括库区、坝体、居民点、河流在内的7个代表点(见图5),即P1位于尾矿砂堆积下游区、P2位于初期坝底、P3位于环保库处、P4位于居民点1、P5位于居民点2、P6位于河流处1、P7位于河流处2,分析溃坝后500 s内每5 s尾矿砂的堆积厚度和运动速度(见图9),可以得知：

(1)在0～500 s的时间内，除位于尾矿库内的P1监测点对应的堆积厚度在逐渐减少，其余监测点的堆积厚度都呈逐步增加的态势，这主要是由于溃坝后库区内的尾矿砂开始迅速往下运动，原始堆积厚度也随之逐渐减少，而相应的位于库区外围的监测点的堆积厚度在增加。

图7 不同时刻尾矿砂流淹没深度纵断面

(2)各个监测点的运动速度也基本上符合先突然增加后缓慢减小最后趋于平缓的规律，主要是因为尾矿库溃坝后尾矿砂拥有一定程度的动能，这种能量促使其以一定速度向下游冲击，随着时间的增加，其运动过程中因受到地面的摩擦，能量逐渐开始消散，从而其运动速度也逐渐减小。

(3)尾矿砂在流经P2(初期坝底)、P3(环保库)、P4(居民点1)以及P5(居民点2)的运动速度能达到15 m/s以上，其中流经居民点1的最大速度可以达到20 m/s以上，这种高速度足以对流经范围内的任何建筑设施进行毁灭性地破坏，若撤离不及时，将造成严重的人员和财产的损失，大约在300 s以后，所有监测点的运动速度都小于5 m/s; 对于P2(初期坝底),由于其地形高程变化较大，使得到达坝底的速度都维持在15 m/s左右。

本文利用Massflow模型对新疆某山区尾矿库进行了尾矿砂一次性全部泄出情景下的溃坝模拟，通过分析溃坝以后尾矿砂在不同时刻的堆积厚度以及运动速度，掌握了溃坝事件引起的下游河道堆积情况及对周边居民点的影响，得出如下结论：

(1)本文基于三维重建地形，将Massflow模型与GIS平台相结合，能够很好地实现尾矿砂流深、流速可视化，从而揭示溃坝后尾矿砂在高分辨率地形下的时空演化过程。

(2)通过Massflow模型分析，该尾矿库溃坝运动历时3000 s, 溃坝后尾矿砂堆积呈现两部分：一部分仍堆积在尾矿库的下游区域，平均厚度约13 m; 另一部分堆积于环保库周围以及河道中，影响范围约234 971.60 m ,平均厚度达25 m, 其中环保库堆积厚度约为21 m, 居民点堆积厚度约为11 m, 堆积在河道中的尾矿砂长约1 100 m, 平均向河道左侧扩张85 m, 堆积厚度最大的位置处于尾矿库初期坝坝底正对的河流处，最大厚度约29 m。

(3)整个溃坝过程可划分为失稳启动阶段、加速下滑阶段以及减速堆积阶段。其中失稳启动阶段和加速下滑阶段都在短时间内迅速完成，减速堆积阶段占整个过程的93.33%。整个运动过程速度呈现一个突然增加然后逐渐减小的趋势，表明尾矿库溃坝将瞬间造成巨大破坏。

图8 尾矿库溃坝关键时间点运动速度

图9 监测点流深、流速随时间的变化(P1-P7)

基于三维重建高分辨率地形模型的Massflow数值模拟得到的尾矿砂的堆积范围和运动速率，可为应急管理部门制定应急响应方案提供有效帮助，再结合下游的建筑物、耕地、人员等信息可对尾矿库溃坝进行风险评估研究，并针对不同的风险情况，制定相应的措施，对潜在危险范围内的建筑物、人员进行有效的保护，极大程度降低财产和人员损失。

水利水电技术

水利部《水利水电技术》杂志是中国水利水电行业的综合性技术期刊（月刊），为全国中文核心期刊，面向国内外公开发行。本刊以介绍我国水资源的开发、利用、治理、配置、节约和保护，以及水利水电工程的勘测、设计、施工、运行管理和科学研究等方面的技术经验为主，同时也报道国外的先进技术。期刊主要栏目有:水文水资源、水工建筑、工程施工、工程基础、水力学、机电技术、泥沙研究、水环境与水生态、运行管理、试验研究、工程地质、金属结构、水利经济、水利规划、防汛抗旱、建设管理、新能源、城市水利、农村水利、水土保持、水库移民、水利现代化、国际水利等。