2017, Vol. 38
2. 北京科技大学 土木与环境工程学院，北京 100083;
3. 渤海装备辽河重工有限公司，辽宁 盘锦 124010
2. School of Civil & Environmental Engineering, University of Science and Technology Beijing, Beijing 100083, China;
3. China Petroleum Liaohe Equipment Company, Panjin 124010, China
尾矿库是金属非金属矿山进行矿石选别后排出尾矿的场所，一般筑坝方式为拦截谷口或围地构成，是维持矿山生产的重要设施.同时尾矿库是一种人造的具有高势能的泥石流，威胁下游居民及设施的安全[1-3].尾矿库的建造已有较长历史，国内外学者对尾矿溃坝生命损失风险的研究很少，主要研究集中在坝体稳定性评价和环境影响.我国不同于其他国家，相当数量的尾矿坝下游就是密集的居民区，一旦溃坝可能引起大量生命伤亡，因此，应重视对尾矿溃坝生命损失风险控制的研究.溃坝生命损失风险控制研究包括坝体稳定性评价、溃坝事故死亡人数估算、确定尾矿溃坝生命损失可接受风险标准及最优安全投入几个方面的内容.
本文从上述几个方面对尾矿坝溃坝生命损失风险控制进行了系统全面的研究.首先选取蒙特卡罗方法对溃坝概率进行计算，建立尾矿坝坡失稳破坏的功能函数，确定模拟参数及模拟次数；其次对溃坝后尾矿砂流的下游演进过程进行模拟，估算生命损失；最后确定适合我国国情的溃坝生命损失风险可接受标准，并引入生活质量指数对溃坝生命价值进行估算，进而获得保证尾矿坝安全的最优安全投入.对尾矿坝溃坝生命损失风险控制进行研究，是减灾防灾工程和进行尾矿库风险管理的必要手段.21世纪以来世界各国越来越重视人的生命价值，生命损失已经成为公众和社会关注的焦点.如果能够提出适合我国的尾矿溃坝生命损失风险控制决策是非常有意义的.
1 基于蒙特卡罗理论的尾矿坝稳定性计算 1.1 基于蒙特卡罗的坝体稳定性分析步骤蒙特卡罗法是通过随机模拟和统计试验来求解结构可靠度的数值方法.具体分析步骤如下：
1) 研究尾矿坝坝坡失稳的破坏模式，选取合适的力学模型进行分析，确定功能函数，这是蒙特卡罗模拟的基础.
2) 找出影响坝坡稳定性的参数，确定这些参数的概率分布和分布参数.
3) 输入坝坡的几何形状，使之投绘在所分析的工程地质剖面上.
4) 在分析剖面上，根据已判定的破坏模式，确定最危险的潜在破坏面.
5) 根据步骤2) 所确定的各随机变量的分布类型和参数，抽样取得一组随机数X1, X2, X3, …, Xi，计算功能函数，Zi≤0，记L=1.
6) 重复步骤5)，直至达到满足预期精度要求的实验次数N, 坝坡失效概率Pf=L/N.
只要N足够大，就能保证所要求的尾矿坝坝坡失稳的溃坝概率与L/N之间的误差可以小于制定误差容许上限.
1.2 确定功能函数和实验次数尾矿坝坝坡是否失稳与筑坝材料的力学性质有着直接密切的关系，郑欣等[4]利用突变学理论研究发现尾矿砂的内摩擦角、黏聚力对坝体抗滑稳定性有显著影响.在可靠性计算时，由于土体的容重γ自身变异系数较小，所以将γ作为确定性变量不会带来过大误差[5-6].进行尾矿溃坝概率分析时将稳定性分析视为随机过程，把影响其稳定性的最主要因素黏聚力和内摩擦角作为随机变量进行蒙特卡罗模拟.简化毕肖普法模型同瑞典圆弧法、摩根斯坦-普利斯法等比较常见的刚体极限平衡法相比更加简单易于计算，而且与精确的严格条分法相比，误差在5%以下[7]，因此主要采用简化的毕肖普法构建功能函数.当不考虑地震的作用，功能函数简化为
|
(1) |
式中：φ′为土条i的有效内摩擦角；c′i为土条i的有效凝聚力；ui为作用于土条i底边上的孔隙水压力；bi为土条i的宽度；Wi为土条所受的重力；αi为土条i的角度.
以往研究结果表明，模拟次数在5 000~10 000次，就可以满足工程评价的精度要求.所以设定蒙特卡罗模拟次数为10 000次对尾矿坝溃坝概率进行计算.
按照确定的功能函数和随机参数，使用Matlab软件并自行编程对溃坝概率进行计算，由于所编程序篇幅过大，在这里就不给出详细的程序代码.
2 尾矿溃坝生命损失估算我国学者还没有对尾矿溃坝生命损失进行系统研究，国外学者在水库溃坝生命损失评估方面已经取得了一些研究成果[8-9]，所以目前对尾矿溃坝事故人员伤亡估算主要还是借鉴水库大坝溃坝生命损失评价的研究成果.李雷等[8]通过研究几种国外的溃坝人员伤亡估算方法, 总结出对生命损失估算影响比较大的几个重要参数.周克发等[10]调查了数座我国已溃水库大坝的现场, 并结合国内外溃坝死亡人数估算的研究成果, 初步改进了我国溃坝生命损失评价方法.我国重大危险源普查填表说明中给出了尾矿库事故可能造成的伤亡人数估算方法.
尾矿溃坝生命损失估算研究主要包括两个方面：溃坝砂流下游演进过程模拟，目的是得到比较符合实际情况的溃坝砂流影响范围 (尾矿坝下游不同位置的尾矿砂流的深度和流速)；溃坝死亡人数估算方法的研究.
对尾矿溃坝砂流下游演进过程的研究主要是借鉴泥石流的研究方法，采用数值方法求解溃坝砂流的控制方程.目前控制方程的数值求解方法主要有特征线法、有限元法、有限差分法、有限体积法等[11].有限体积法由于在国外研究的比较早，目前较为完善，一些大型流体计算软件如Fluent，Star-CD等也运用了有限体积法[12-13].本文主要采用Fluent软件对溃坝砂流演进进行模拟，由于尾矿砂粒径通常比较小，所以把溃坝后的尾矿和水的混合物看作是“液相”一相流，溃坝砂流在下游的演进过程当作气液两相流进行模拟，自由表面的追踪采用VOF法[14].
我国重大危险源普查填表说明中给出的尾矿坝溃坝事故可能造成的死亡人数简化公式为
|
(2) |
式中：S为尾矿库事故可能造成的死亡人数；Ni为第i个居民点的居民人数；i为尾矿坝下游10倍坝高范围内，n个居民点的顺序数；Mj为第j个居民点的居民人数；j为尾矿坝下游10倍坝高以外，40倍坝高以范围内，m个居民点的顺序数.
3 尾矿溃坝生命损失可接受风险标准 3.1 尾矿溃坝生命损失可接受风险国外很多国家已制定了国家和行业的可接受风险标准.目前我国很少学者创新性地建立了我国的可接受风险标准，主要是尝试着通过把国外的方法与我国的实际国情结合起来，确定我国的风险可接受标准，但更多的则是直接采用国外的风险标准.
不存在可接受风险，所谓的可接受风险是所承担的风险与所获收益权衡后的结果.可接受风险标准的确定过程是一个决策的过程, 所以必须是企业和公众共同来制定尾矿溃坝生命损失可接受风险标准.首先确定一个尾矿溃坝生命损失可接受风险的基准值，在很多会议中，超过半数以上赞成才算通过，由此可见1/2通常被作为保证质和量的界限，所以以基准值为中心，上下浮动1/2，得到一个尾矿坝溃坝生命损失风险可接受标准的协商区域，然后由企业和公众在这一区域进行讨论，最后确定可接受风险标准.通过这种方式得到的尾矿坝溃坝生命损失可接受风险标准比较客观，不会偏离实际风险.
通过分析我国从2003年到2012年的尾矿坝溃坝死亡人数, 确定我国尾矿溃坝生命损失可接受风险标准的基准值.具体统计情况见表 1.
|
表 1 我国尾矿溃坝死亡人数统计表 Table 1 Number of casualties of tailings dam break in China |
为了保证数据的平稳性，避免异常数据对研究结果的影响，尾矿溃坝生命损失可接受标准基准值的确定可采取中位数法，5×10-4人/(a·座) 为整组数据的中位数，但直接用过去10 a的统计数据得到的该值作为当前的可接受风险标准, 显然确定的标准值过低，因此需要参考国外现有的风险标准并适当提高该标准值，最终确定0.000 1人/(a·座) 为我国尾矿溃坝生命损失可接受风险标准的基准值.风险提供者和风险承受者可在如下的范围内协商确定最终的尾矿溃坝生命损失可接受风险标准：
|
Nathwani等[15]发展并使用了由Lind提出的生活质量指数 (LQI)，见式 (3).采取措施减小项目风险时，LQI会受到影响，见式 (4)[16]：
|
(3) |
式中：L为生活质量指数LQI；g为人均年国民生产总值；e为预期寿命；ω为工作在预期寿命中占的比例.
|
(4) |
式中：ΔL为生活质量指数的变化量；Δg为人均国民生产总值的变化量；Δe为预期寿命的变化量.
当LQI指数的变化为非负数时，见式 (5)[15], 增加安全投入意味着减少人均国民生产总值，式 (5) 可以转化为式 (6)[15], 当式 (6) 左右两端相等时可以得到安全投入的最优方案.
|
(5) |
|
(6) |
某尾矿堆积坝顶达196 m标高，筑坝总长度900 m, 平均坝坡为1:4.5，现滩面干坡段长度可达90 m.该尾矿坝下游村庄中的巨各庄 (约500人) 与尾矿坝的距离在10倍坝高以内，豆各庄 (约200人)、赵家庄 (约200人)、金山子村 (约200人) 与尾矿坝的距离为10倍坝高以外.表 2为该尾矿坝体材料力学参数.
|
|
表 2 尾矿砂力学参数 Table 2 Mechanical parameters of tailings sand |
首先确定最危险滑弧所在的位置，将土条分成30块, 土条宽度为5 m.选取黏聚力和内摩擦角作为模拟参数，考虑黏聚力与摩擦角的相关性, 从计算结果偏于安全考虑，黏聚力和内摩擦角采用正态分布的概率模型.模拟次数定为10 000次.使用Matlab软件，通过编制的程序对溃坝概率进行计算，经计算该尾矿坝溃坝概率为0.001 6.
本文的三维溃坝砂流演进数学模型以Navier-Stokes方程为基础，紊流方程为封闭模型，其基本微分方程包括连续性方程、动量方程、紊动能方程和紊动能耗散率方程；气液两相的自由表面的追踪采用VOF法；设定溃口宽度为238 m，下泄砂流总量 (假定尾矿坝内全部尾矿下泄) 为5×106m3，下泄泥砂密度为1.6 t/m3；采用Herschel-Bulkely黏度模型，该模型中的稠度系数为0.024，流性系数为0.7，屈服应力为0.923, 糙率系数为0.02, 时间步长设置为0.01 s；结果的收敛性很好，各项残差值均在0.1以下，模拟结果如图 1所示.模拟过程中没有考虑下游地形情况对模拟结果的影响，因此模拟得到的影响范围会比实际大.如果较少地考虑应急救援所需要的经济支出，更多地考虑救援准备的充足性，那么以此结果为依据进行应急救援准备会更合理.
|
图 1 砂流影响范围 Fig.1 Influence scope of discharged tailing flow |
尾矿坝溃决后50~100 s内尾矿砂流覆盖了巨各庄、豆各庄及赵家庄，溃决100~150 s内尾矿砂流到达了金山子村.巨各庄、豆各庄、赵家庄及金山子村最大泥深约为4 m，最大速度为20 m/s，最终泥深约为2 m.利用重大危险源普查中提供的溃坝后死亡人数估算公式计算该尾矿坝溃坝后可能的生命损失为325人.从溃坝砂流演进模拟可知一旦发生溃坝，即使及时得到警报，人员的逃生机会也是微乎其微，人员的死亡率几乎为100%.利用重大危险源普查中提供的方法进行估算的生命损失过于保守，因此对该尾矿坝溃坝生命损失的估算建议参考溃坝砂流下游演进的模拟结果，生命损失为1 100人.对该尾矿坝来说风险远远超过了风险可接受标准，风险不可接受.
把溃坝生命损失现实风险降低到可接受风险标准范围内的最优安全投入计算如下：
2014年人均国民生产总值g=46 600元，预期寿命e=71.4 a，ω=0.174 6(其中用于工作的时间按35 a计算)，由于死亡人员的年龄各异，建议个人预期寿命增加的平均年数为e/2.挽救一个人生命损失的安全成本C为社会挽救生命一年最优安全成本Δg与预期寿命增加Δe的乘积，即
|
(7) |
从计算结果可以看出，保证一个30岁左右的员工安全地生存到预期寿命的最优安全投入约为390万人民币.使用同样的方法计算出美国避免一个人因事故死亡的最优安全成本为600万美元，计算结果与美国实际的赔偿金额是基本吻合的.2008年美国环保署把命价从600万提高到910万美元，2014年美国纽约州皇后区法院更是做出了一项高达6 200万美元的工伤赔偿裁决，这些都体现出了世界各国都越来越重视人的生命价值.我国目前的工伤死亡赔偿金额仅为60万左右，远远低于本文计算出的390万的最优安全投入.尽管人的生命是无法用金钱来衡量的，但事故的发生往往是必然的，所以有必要对人的生命价值进行研究，得到较为合理的伤亡赔偿金额，这才是对生命真正的尊重.建议国家相关部门提高我国目前的工伤死亡赔偿金额.
用于挽救该尾矿坝下游1 100名居民生命的最优安全投入为3 932 279×1 100≈43亿元.显然这笔安全投入是巨大的，所以建议积极采取安全措施提高该尾矿坝的稳定性避免溃坝事故的发生，同时在尾矿坝下游栽种树木及构筑可以拦截溃坝砂流的防护体，保障下游居民的人身财产安全.
5 结论1) 根据建立的尾矿坝坡失稳破坏的功能函数，选择蒙特卡罗方法计算出尾矿坝的溃坝概率为0.001 6.
2) 对溃坝砂流演进过程使用fluent软件进行了两相流三维模拟，得到了溃坝后尾矿砂流在下游的影响范围及到达下游各居民点的时间、最大泥深及最大速度.根据模拟结果和重大危险源普查两种方式对溃坝后可能的生命损失进行估算，重大危险源普查中提供的方法估算的结果过于保守，最后根据fluent模拟结果对下游生命损失进行估算，该尾矿坝一旦溃坝，生命损失为1 100人.
3) 依据建立的我国尾矿溃坝个人生命可接受风险并采用基于LQI的方法对生命价值即保证员工安全的最优安全投入进行估算，计算得到的最优安全投入为390万人民币.并将美国和我国现行的工伤死亡赔偿标准进行比较研究，建议国家相关部门提高我国现行的工伤死亡赔偿标准.
4) 从对尾矿坝的风险评价研究中可以看出，有些尾矿坝的溃坝概率虽然很低，但我国大部分尾矿库的下游就是密集的居民居住区，一旦发生溃坝事故，生命损失是巨大的，现实风险远远超过可接受风险标准.因而对尾矿库进行风险管理，不仅要考虑如何采取措施提高坝体稳定性，还要弄清生命损失状况及存在规律，提出适合我国国情的溃坝生命损失风险控制决策.
| [1] | Nazli T O, Resat U, Nihat S I. A study on geotechnical characterization and stability of downstream slope of a tailings dam to improve its storage capacity (Turkey)[J]. Environmental Earth Sciences, 2013, 69: 1871–1890. DOI:10.1007/s12665-012-2016-1 |
| [2] | Strachan C. Tailings dam performance from USCOLD incident-survey data[J]. Mining Engineering, 2001, 53(3): 49–53. |
| [3] | Sharma R S, Al-Busaidii T S. Groundwater pollution due to a tailings dam[J]. Engineering Geology, 2001, 60(1): 235–244. |
| [4] |
郑欣, 秦华礼.
基于耗散结构理论的尾矿坝系统研究[J]. 金属矿山, 2009(01): 159–161.
( Zheng Xin, Qin Hua-li. Study of tailing dam system based on dissipative structure theory[J]. Metal Mine, 2009(01): 159–161. DOI:10.3321/j.issn:1001-1250.2009.01.045 ) |
| [5] | Carson M A, Kirk M J. Hill slope form and process[M]. Cambridge: Cambridge University Press, 1972. |
| [6] |
徐茂其, 张大泉. 斜坡稳定性综合评价探讨[C]//全国自然边坡稳定性分析暨华蓥山边坡变形研讨会. 北京, 1991: 112-117.
( Xu Mao-qi, Zhang Da-quan.Comprehensive evaluation of slope stability [C]// Analysis of China Natural Slop Stability—Huaying Mountain Slope Deformation Seminar.Beijing, 1991:112-117. ) |
| [7] |
任克峰. 边坡稳定性评估方法及应用[D]. 宁夏: 宁夏大学, 2008, 27-28.
( Ren Ke-feng.Slope stability evaluation method and its application [D].Ningxia:Ningxia University, 2008:27-28. ) |
| [8] |
李雷, 周克发.
大坝溃决导致的生命损失估算方法研究现状[J]. 水利水电科技进展, 2006, 26(2): 76–80.
( Li Lei, Zou Ke-fa. Methods for evaluation of life loss induced by dam failure[J]. Advances in Science and Technology of Water Resources, 2006, 26(2): 76–80. ) |
| [9] | Graham W J. A procedure for estimating loss of life caused by dam failure[J]. Sedimentation & River Hydraulics, 1999, 6(5): 1–43. |
| [10] |
周克发, 李雷.
我国溃坝生命损失评价模型初步研究[J]. 安全与环境学报, 2007, 7(3): 145–149.
( Zou Ke-fa, Li Lei. Evaluation model of loss of life due to dam breach in China[J]. Journal of Safety and Environment, 2007, 7(3): 145–149. ) |
| [11] |
谭维炎.
浅水动力学的回顾和当代前沿问题[J]. 水科学进展, 1999, 10(3): 296–303.
( Tan Wei-yan. A review of shallow-water hydrodynamics and some contemporary problems[J]. Advances in Water Science, 1999, 10(3): 296–303. ) |
| [12] | Levque J. Finite volume methods for hyperbolic problems[M]. Cambridge: Cambridge University Press, 2002. |
| [13] | Toro E F. Shock-capturing methods for free-surface shallow flows[M]. London: John Wily & Sons Ltd, 2001. |
| [14] | Hirt C W, Nichols B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]. Journal of Computational Physics, 1981, 39(1): 201–225. DOI:10.1016/0021-9991(81)90145-5 |
| [15] | Nathwani J S, Lind N C, Pandey M D. Affordable safety by choice:the life quality method[M]. Waterloo: University of Waterloo, 1997: 52-187. |