0.   引言

中国地域辽阔,山脉众多,矿产资源十分丰富,矿石经过冶炼获得有用成分之后,会产生大量无法利用的尾矿,尾矿中富含重金属及有毒有害物质,为防止其任意排放造成大面积环境污染,一般需建立尾矿库集中对其进行存放。随着社会经济的快速发展,大量矿山被开采利用,尾矿库数量和规模也随之增加,目前我国尾矿库数量已有上万座,其中95%的尾矿库采用上游法筑坝。近年来,受到地震、降雨等的影响,尾矿库溃坝灾害频繁发生,危害十分严重,尤其是头顶库,对下游人民生命财产具有重大的威胁。

为保证尾矿坝渗流、变形及稳定性分析的准确性,合理确定尾矿材料的取值尤为重要。目前尾矿材料参数取值通常是通过三轴试验、共振柱试验和动单剪试验等室内试验获得。虽然通过试验测得的尾矿坝材料参数相对准确,能够较好反映尾矿材料的力学特性,但是费时费力,容易受到试验人员素质及试验数据处理方法的影响,导致参数取值与客观实际之间存在误差^[1]。相比之下,将现场实测数据与反演分析有机结合为确定尾矿坝材料参数取值提供了一条有效的途径^[2-7]。如林旭等^[2]运用BP神经网络模型基于实测数据对其尾矿坝内部初始地应力进行反演分析。徐维生等^[3]通过实测水头数据对尾矿坝渗透系数进行反演分析,为尾矿坝的安全运行和加高分析提供了数据来源。任赞松^[4]对洛阳大沟河尾矿库材料参数进行反演分析和溃坝模拟,为特大型尾矿库安全运营管理提供了第一手基础数据。李彪^[5]依据现场勘查数据对尹庄尾矿坝进行渗流场反演分析。Bharat等^[6]通过反演分析方法确定了尾矿材料压缩模量。金佳旭^[7]通过Fortran语言和VC++技术开发了尾矿砂孔压增长模型参数反演软件,对尾矿砂孔压增长应变模型进行参数计算,为研究尾矿库动力特性提供了基础资料。上述参数反演研究提高了尾矿材料取值精度,为尾矿坝安全运营提供了数据来源,但是大多将尾矿材料考虑为单一土层,忽略了尾矿材料参数不确定性的客观存在,均属于确定性反演分析范畴。

尾矿坝由于长期经受多循环水力充填及固结沉降作用,其物质组成和内部结构会发生一定的变化,尾矿材料存在较大的离散性和时空变异性。同时现场实测数据又非常有限,在获取尾矿材料参数时存在一定的测量误差和统计误差^[8-9]。本文结合贝叶斯更新理论和有限元分析提出了考虑不确定性的尾矿材料参数随机反演方法,以实际尾矿坝工程为例,基于有限的位移监测数据对多层尾矿材料参数进行随机反演分析,推断尾矿坝材料参数的概率分布,并识别不同材料参数（弹性模量和泊松比等）对尾矿坝变形的影响程度。

1.   尾矿坝材料参数随机反演方法

尾矿坝材料参数普遍存在一定的时空变异性、测量误差和统计误差^[8-9]。为获得精确的尾矿坝稳定性评估结果,需要从特定的场地收集现场试验数据及监测数据等进行参数概率反演分析来确定尾矿坝材料参数取值,这一过程可以通过贝叶斯分析估计参数的后验概率分布函数${f^{″}}_X(x)$。根据贝叶斯理论,${f^{″}}_X(x)$的计算表达式为^[10-11]

式中 ${f^{\prime }}_X(x)$表示尾矿材料参数的先验联合概率密度函数；a为比例常数,以保证在X的整个区域上对${f^{″}}_X(x)$的积分为1.0；L(x)为似然函数,表示在已知材料参数x=X前提下场地信息事件发生的概率。然而,一般情况下式（1）没有解析解,需要数值求解,本文采用Straub等^[10]于2015年提出的基于结构可靠度分析的贝叶斯更新方法,有效求解计算${f^{″}}_X(x)$数值解。该方法是将似然函数定义为一个新的失效区域来建立结构可靠度与贝叶斯更新之间的关系,将复杂的贝叶斯更新问题转化为一个等效的结构可靠度问题。然后采用子集模拟^[12]求解该等效的结构可靠度问题。首先,需要基于似然函数与似然函数有关的常数v,定义一个新的场地信息失效区域${\Omega }_X$：

式中,w为[0,1]均匀分布随机变量的模拟值,v为似然函数乘子。为了保证获得的样本服从目标概率分布,对所有的x,都要满足vL(x)≤1.0。

接着,采用子集模拟求解式（2）所示的等效的结构可靠度问题,其中驱动变量Z与似然函数乘子v有关,子集模拟计算之前需要提前确定v值,v值的确定对参数随机反演计算精度和效率有重要的影响。为了保证计算精度和效率同时达到最佳,文献[11, 13]将场地信息失效区域${\Omega }_X$调整为

相应的子集模拟驱动变量Z调整为

驱动变量Z采用自然对数的形式是为了加快计算速度,则式（2）所示的目标失效区域转换为${\Omega }_X$ _{= {Z > b},对应的阈值b定义为}

改进的贝叶斯更新方法虽然无需提前确定v值,但是合理确定阈值b的值尤为重要。为了推断获得参数后验概率分布,b必须要大于b_min^{[11, 13]},b_min的计算表达式为

虽然b_min事先也不得知,但是这并不影响后续子集模拟计算和参数随机反演分析,可将P(Z > b)表达为一系列较大的中间事件条件概率的乘积：

式中 $P({\Omega }_{X_1})$为${\Omega }_{X_1}$的发生概率；$P({\Omega }_{X_i}|{\Omega }_{X_{i-1}})$为在${\Omega }_{X_{i-1}}$发生的条件下${\Omega }_{X_i}$发生的概率,i =1, 2, …, m,m为子集模拟达到失效区域${\Omega }_X$所需模拟的层数。中间事件${\Omega }_{X_i}=\{Z>b_i\}$,b_i为临界阈值。随着子集模拟层数的增加,当第m层子集模拟阈值b_i>b_min时,表示达到目标失效区域可终止子集模拟计算,然后从子集模拟最后一层中提取失效样本。

另外,为保证从第i层子集模拟中获得的失效样本服从目标概率分布,必须同时满足vL(x)不大于1.0和b_i >b_min这两个条件,即

若$B_i=\left\{e^{-b_i}L(x)>1.0\right\}$趋于空集,相应的概率$P(B_i)$接近于0,则式（8）便无条件满足,可终止子集模拟计算,从最后一层提取失效样本,推断尾矿坝材料参数后验概率分布并估计其后验统计特征。

式（1）中的似然函数表示尾矿坝材料参数与观测样本之间的概率关系以及实测数据与后验分布之间的联系,用于刻画场地信息。似然函数表示当尾矿材料参数取某一特定值时,与现场实测数据发生的概率成正比。对于不同类型的观测样本,需要建立多个与之对应的似然函数。当已知某点位移监测数据时,表征尾矿材料参数与位移监测样本之间概率关系的似然函数可构建为^[14]

式中 d为某监测点的位移监测值；u(x)为通过数值模拟软件或代理模型计算的某监测点的位移值；$\epsilon$为测量误差,通常模拟为均值为0标准差为某一常数${\sigma }_{\epsilon }$的正态分布^[14]；ϕ(.)为标准正态概率密度函数。当有多个点的位移监测数据时,应建立包含多个点监测数据相关性的联合似然函数为^[11]

式中,d为实测得到的某点位移值,q为实测样本量大小,协方差矩阵${\sum }_d$构成为

其中,${\rho }_{{\epsilon }_i,{\epsilon }_j}$表示不同点处监测数据测量误差之间的相关系数,为简化计算,取${\rho }_{{\epsilon }_i,{\epsilon }_j}$=0（i ≠ j）来推断参数后验分布。当然也可以将测量误差模拟为随机场,采用自相关函数及不同监测点的位置来计算测量误差之间的相关系数。

2.   工程案例

2.1   尾矿坝基本概况

某实际尾矿坝^[15]位于吉林省,该尾矿库与矿场的南侧紧密相连,最低海拔标高360 m左右,地处山区丘陵。库区位于东北,属于温带大陆性气候,年平均温度较低为5℃,年均降雨量649.6 mm。该尾矿库按照上游法建造,属于山谷型尾矿库,设计等别为三等尾矿库,初期坝采用堆石体砌筑,初级坝顶高程410 m,宽4 m,内坡与外坡比均为1∶2,坝基处高程为403 m。尾矿材料逐级往上堆积,堆积到最大坝高时的库容为5800×10⁴ m³,该尾矿库仍在运行当中,目前坝顶高程为470 m。计算模型长度为690 m,最低高程为367 m,具体剖面见图1。受管道水力冲击影响,一般越靠近输送口位置的尾矿材料越粗糙,越远离管口越细。在堆积某一层尾矿材料时,管口位置的移动也是随机的,排放一次就会形成小山包,向四处滑落,从而形成了尾矿材料之间错综复杂的关系。将尾矿材料概化为粗砂、细砂、粉砂、粉土4种土层,尾矿材料名称及分布见图2。

图  1  尾矿坝剖面图

Figure  1.  Profile of tailings dam

下载: 全尺寸图片 幻灯片

图  2  尾矿坝有限元模型及监测点布置

Figure  2.  Finite element model for tailings dam and layout of monitoring points

下载: 全尺寸图片 幻灯片

2.2   计算模型及参数统计特征

采用岩土工程有限元应力应变分析软件GEOSTUDIO中的SIGMA/W模块^[16]对黑山尾矿坝建立有限元模型,基岩下方是坚固岩石,可将基岩底部简化成水平线。有限元模型共离散了53425个单元和52723个节点。采用三角形和四边形网格,网格尺寸为1 m。先使用线弹性材料模型进行自重作用下尾矿坝原位分析,再将尾矿材料视为弹塑性模型进行荷载变形分析。位移监测点分布在430和445子坝上,各布置1个监测点,具体位置如图2所示。整个尾矿坝模型上部不施加任何约束,水平方向设置x方向的位移约束,基岩底部设置x,y方向的位移约束。

根据大黑山尾矿库岩土工程勘察报告,充分考虑放矿形式、时间效应和尾矿特性,并且参考借鉴前人学者对尾矿坝材料参数的试验研究成果,获得了不同尾矿材料物理力学参数取值,如表1所示。

表  1  尾矿材料物理力学参数

Table  1.  Physicomechanical parameters of tailings materials

土层	弹性模量/MPa	泊松比	黏聚力/kPa	内摩擦角/(°)	重度/(kN·m-³)
堆石	120	0.22	42.0	35	21.5
粗砂	50	0.30	7.9	32	22.3
细砂	40	0.30	8.1	28	21.5
粉砂	40	0.30	9.0	25	21.8
粉土	30	0.35	9.7	29	20.2
基岩	7700	0.28	30000	51	24.0

下载: 导出CSV 

| 显示表格

尾矿坝堆积过程中处于无其他荷载,只受重力作用的影响,尾矿坝材料主要是粗砂、粉砂等,发生变形所涉及到的主要参数应为刚度参数,故选取不同尾矿坝材料的弹性模量与泊松比来进行随机反演分析。基岩基本不发生变形,故在参数随机反演中不考虑基岩的作用,只对粗砂、细砂、粉砂、粉土、堆石体的刚度参数进行随机反演分析。尾矿坝弹性模量与泊松比的先验均值和变异系数借鉴相关土体材料参数资料及相关文献^[8-9]取值,不同尾矿材料参数（弹性模量、泊松比）均假设服从对数正态分布,先验均值和变异系数如表2所示。

表  2  不同尾矿材料参数先验和后验统计特征的比较

Table  2.  Comparison of prior and posterior statistics of parameters of tailings materials

土层	弹性模量先验信息		弹性模量后验信息		泊松比先验信息		泊松比后验信息
	均值/MPa	变异系数	均值/MPa	变异系数	均值	变异系数	均值	变异系数
堆石	120	0.15	124.40	0.144	0.22	0.1	0.221	0.100
粗砂	50	0.15	51.779	0.069	0.30	0.1	0.299	0.087
细砂	40	0.15	37.330	0.129	0.30	0.1	0.285	0.091
粉砂	40	0.15	39.871	0.150	0.30	0.1	0.302	0.099
粉土	30	0.15	30.221	0.149	0.35	0.1	0.350	0.100

下载: 导出CSV 

| 显示表格

2.3   应力变形分析

尾矿坝由尾矿材料逐级堆积形成的,堆积时先设置堆石体,防止后续堆积尾矿材料溢出。在堆积尾矿时,由堆积点向四周扩散,对周围有加载作用,进行尾矿坝变形分析时,以尾矿坝堆积高度470 m作为基准。通过计算分析发现x方向位移呈现上游小,下游较大的趋势,子坝以下部分x方向位移最大可达68 mm,整个库区呈非对称形态,符合尾矿坝上游有库区阻挡,向下游滑动的趋势。在自重作用下,尾矿坝高程越高,y方向沉降量越大,在坝顶位置y方向最大沉降值可达720 mm。故在自重作用下,尾矿坝总体呈现向下沉降,并向下游微小滑动的趋势。

2.4   参数随机反演分析

参数随机反演分析过程中通常需要进行大量的有限元计算尾矿坝代表性监测点的位移值构建似然函数,为提高反演分析计算效率,在进行尾矿材料参数随机反演分析之前,采用3阶Hermite随机多项式展开建立2个监测点处尾矿坝位移的代理模型^[17]。接着根据文献[18],取测量误差标准差${\sigma }_{\epsilon }$= 2.0 mm,基于获取430子坝D1监测点和445子坝D2监测点的位移值（70.6,220.35 mm）构建似然函数。进而采用提出方法推断不同尾矿材料弹性模量与泊松比的后验统计特征及概率分布。其中,贝叶斯更新中每层子集模拟样本数目为2000,条件概率取10%,重复独立计算10次子集模拟,取均值作为最终的计算结果^[11]。5种尾矿材料中,粗砂弹性模量的概率密度函数变化最大,细砂次之,而堆石体、粉砂、粉土弹性模量的概率密度函数变化很小,如图3所示。表明融合现场位移监测数据进行参数随机反演分析之后,堆石体、粉土和粉砂弹性模量的不确定性并没有得到有效缩减,这是因为粉土、粉砂距离监测点较远,因而对基于监测点位移的弹性模量随机反演的影响较小；而堆石体因为其弹性模量较大,自身变形量较小,所以也对这两监测点尾矿坝的变形影响较小。相比之下,不同尾矿材料的泊松比的概率密度函数也基本没有变化,说明泊松比对尾矿坝变形的影响程度较小。同时,表2进一步比较了不同尾矿材料弹性模量与泊松比先验和后验统计特征。由表2可见,粗砂、细砂弹性模量变异系数均有不同程度的缩减,其中粗砂弹性模量变异系数缩减的程度最大,进一步说明通过融合位移监测数据进行参数随机反演分析能够更加准确地表征尾矿坝材料参数的不确定性。

图  3  不同尾矿材料弹性模量先验和后验概率密度函数的比较

Figure  3.  Comparison of prior and posterior probability density functions of elastic modulus of tailings materials

下载: 全尺寸图片 幻灯片

3.   结论

本文在贝叶斯更新理论和有限元分析框架下提出了考虑不确定性的尾矿材料参数随机反演方法。基于有限的位移监测数据,对某实际尾矿坝5种材料参数进行随机反演分析。

（1）基于2个监测点的位移的数据进行尾矿坝材料参数随机反演,发现监测点的位置对参数随机反演有较大的影响,距离监测点较远的粉土、粉砂反演之后弹性模量与泊松比的概率分布基本不变,而离监测点位较近的粗砂、细砂经反演之后弹性模量与泊松比的不确定性有明显的缩减。另外,堆石体的泊松比和弹性模量的概率分布也几乎没有变化,这是因为堆石体本身的弹性模量较大,自身变形较小。

（2）提出方法可以有效基于位移监测值对尾矿坝5种材料进行参数随机反演分析,获得的粗砂、细砂弹性模量的变异系数的缩减程度明显大于泊松比,说明弹性模量对尾矿坝的变形影响程度更大。提出方法还可拓展到土石坝、堤防等水工岩土结构参数随机反演中,为确定有限数据条件下材料参数取值提供有效的工具。