0.   引　　言

为了贯彻我国生态文明建设和高质量发展政策，露天矿作为能源消费结构的重要组成部分，必须在发展过程中实现绿色、安全及高效。在露天矿区，外排土场是必不可少的场所，用于排放露天矿开采过程中产生的大量土石方。为了减少土地资源占用，通常采用增高扩容方式最大限度地利用外排场地。然而，这种做法往往会导致排土场的稳定性下降，周边地表变形加剧。尤其当排土场周边存在高速铁路、公路等刚性建构筑物时，变形引发的安全问题就变得更加突出，必须引起足够的重视。

国内外许多专家学者在解决土体变形问题方面采用了层状弹性力学理论，并做出了许多创新性的理论推导。TERASAWA等^[1]、IX^[2]、ACUM等^[3]、XPYCTAЛEB等^[4]、DE等^[5]不同的理论方法对双层和三层弹性体系在荷载作用下的应力和位移变形等问题进行了深入研究和分析，并提出了相应的理论公式和计算程序。国内学者朱照宏等^[6]将Love位移函数引入到高阶矩阵代数法中进行分析，求出基于单个圆形垂直均布荷载作用下双层、三层弹性体系解析解。吴晋伟^[7]对多层弹性体系进行了应力与位移分析，得出存在于多层体系中的应力与位移的计算公式，并编制了多层路面中的应力与位移计算程序。梁锡三^[8]以层间的等效代换为基础，提出了圆形多层体系在竖向均布荷载作用下的一些常用坐标点的竖向位移、竖向应力的解析方法。王凯^[9-13]采用递推回代法，完成了多层弹性体系应力与位移分析。同时，姜朋明^[14]、支喜兰等^[15]还采用了传递矩阵法、黏弹性层状地基轴对称问题等不同的方法进行了研究。这些研究为计算地表变形方法提供了解析基础。也有许多学者采用理论推导与数值模拟相结合的方法研究基底土体力学特性。吕惠卿等^[16]结合宏、细观力学对弹性层状体系理论进行综述。杨国良等^[17]利用层状弹性体系理论和冲击荷载沿路面结构下传扩散的特点建立了路表变形响应与模量反算模型。艾智勇等^[18]通过推导单层黏弹性介质的解析层元解，得到多层黏弹性体系总刚度矩阵并分析了黏弹性系数对路基沉降的影响。阳恩慧等^[19]推导了直角坐标系下的三维多层弹性层状体系静力学数值解法，并借助ABAQUS模拟验证了其合理性。刘光伟等^[20]通过蠕变模型对排土场填方体细粒土质砂蠕变参数进行反演，发现在稳态蠕变阶段中数值计算与现场监测结果相吻合。毕仲辉等^[21]基于层状地基模型，通过地震动输入方法，发现层状边坡会产生更大的塑性应变，且动力安全系数更小。曹博等^[22]通过建立排土场沉降蠕变理论模型，得到了填方土蠕变沉降方程和下沉速度系数，和Burgers黏弹性蠕变模型模拟结果具有较好的一致性。秦哲等^[23]基于弹黏塑性理论建立了一种新型非线性蠕变损伤模型，该模型既能反映整个蠕变过程中损伤演化规律，还能准确描述加速蠕变阶段非线性变形规律。周瑞鹤等^[24]基于分数阶导数引入黏塑性蠕变启动元件，建立的非定常分数阶黏弹塑性卸荷蠕变模型，可较好表征该类岩石的蠕变特征。

在现有的研究中，大多数学者都集中于针对均布载荷作用下的地基弹性变形解法进行研究，虽然也有一些学者采用了黏弹性解法进行求解，但对于像排土场这类三角或梯形载荷作用下的土体黏弹性变形，尚未建立可推广的成型解法。为此，笔者首先基于层状弹性力学理论，推导了层状弹性排土场周边地表变形预测方法，然后引入广义Kelvin三轴蠕变模型，推导了具有流变特性的软弱基底排土场周边地表变形解析预测方法，为排土场周边土体变形预测、排土场设计以及周边设施的安全评价提供了科学手段。

1.   水平层状基底排土场周边地表变形黏弹性解法

在工程实践中，通常排土场基底土体在受到荷载作用时，初期会表现出近似线性弹性形变，这意味着土体会迅速恢复原状，没有持久的变形。然而，随着荷载的持续作用，逐渐观察到基底土体开始出现明显的持久性变形，这表现为土体的流变特性。弹性解法虽然能够很好地描述土体的瞬时响应，但它们无法准确捕捉基底土体的流变行为和时效性变形。

因此，在考虑排土场周边地表变形时，笔者首先考虑基底层状弹性解法，再引入黏弹性本构模型，以准确地描述基底土体的流变特性，更全面地理解和预测排土场基底土体的复杂行为。通过推导排土场周边地表的黏弹性解法，能够更好地理解和分析变形的时效性，这不仅在理论上具有重要意义，还对实际工程应用具有重要价值。

1.1   水平层状基底排土场周边地表变形弹性解法

排土场基底由多层不同岩性和厚度的土体组成，多可将其视为层状弹性体系，其中各土层的厚度、弹性模量和泊松比分别为h_i、E_i、μ_i，（i=1，2，…，N−1）。当基底为半无限弹性空间体时，其为N=1的特例。排土场荷载空间形态多为梯形台拟柱体^[25]，与弹性基底都对称通过排土场中心的某一轴。因此，轴对称弹性体内仅产生对称于该轴的应力和位移分量，可视为空间轴对称问题。通过建立柱坐标系rθz，可以方便地描述应力和位移，排土场与层状基底力学模型如图1所示，并用Lame-Navier方程式^[26]求解空间轴对称弹性力学变形问题。

图  1  排土场与层状基底力学模型

Figure  1.  Mechanical model of dump and layered basement

下载: 全尺寸图片 幻灯片

式中：$ {\nabla ^2} = \dfrac{{{\partial ^2}}}{{\partial {r^2}}} + \dfrac{1}{r}\dfrac{\partial }{{\partial r}} + \dfrac{1}{{{r^2}}}\dfrac{{{\partial ^2}}}{{\partial {\theta ^2}}} + \dfrac{{{\partial ^2}}}{{\partial {z^2}}} $为Laplace算子；e为总应变；$e = {\varepsilon _r} + {\varepsilon _z} + {\varepsilon _\theta } = \left( {\dfrac{\partial }{{\partial r}} + \dfrac{1}{r}} \right)u + \dfrac{1}{r}\dfrac{{\partial v}}{{\partial \theta }} + \dfrac{{\partial w}}{{\partial z}}$，m；u,v,w分别为r,θ, z方向的位移；ε_r为r方向的正应变；ε_θ为θ方向的正应变；ε_z为z方向的正应变。

处于轴对称情况下，设调和函数Φ=Φ(r,z)，∂Φ/∂θ=0，当Φ满足重调和方程：

则方程式(1)的解，z方向位移分量的表达式^[26]为

式中：E为弹性模量，MPa；μ为泊松比。

排土场基底成层状分布，对于其中某一i土层，运用汉克尔积分变换，将式(2)两边各乘rJ₀(${\xi r} $)并由0～∞对r求积分，可得

式中：J₀(${\xi r} $)为0阶第一类贝塞尔。

重调和方程式(2)转化为常微分方程，解得：

式中： $ A_i^{^{(1)}} $，$ B_i^{^{(1)}} $ ，$ C_i^{^{(1)}} $，$ D_i^{^{(1)}} $分别为积分常数。

根据汉克尔积分变换得逆变换式，由式(5)可得：

将式(5)中调和函数Φ的常微分解$\bar \varphi $(r,z)对z的各阶导数代入式(6)，并令式中$ {A^{(0)}} = {\xi ^{ - 2}}{A^{(1)}} $，$ {B^{(0)}} = {\xi ^{ - 3}}{B^{(1)}} $，$ {C^{(0)}} = {\xi ^{ - 2}}{C^{(1)}} $，$ {D^{(0)}} = {\xi ^{ - 3}}{D^{(1)}} $，即得到任意排土场荷载作用下N层基底内第i层土层应力和位移分量表达式：

式中：${A}_{i}^{(1)} $，${B}_{i}^{(1)} $，${C}_{i}^{(1)} $，${D}_{i}^{(1)} $为积分变量的函数。

令$ A_i^{(1)} = {A_i}\overline q \left( \xi \right) $，$ B_i^{(1)} = {B_i}\overline q \left( \xi \right) $，$ C_i^{(1)} = {C_i}\overline q \left( \xi \right) $，$ D_i^{(1)} = {D_i}\overline q \left( \xi \right) $，其中q为作用在弹性体系表面的圆形均布垂直荷载的荷载集度，荷载圆半径为δ，当受垂直均布荷载作用时，$ \overline q \left( \xi \right) = q\int_0^\delta {r{J_0}\left( {\xi r} \right)} {\mathrm{d}}r $，引入新积分变量x=ξδ。于是有：

排土场垂直荷载作用下，地表面应力边界条件通过汉克尔积分变换得：

在排土场荷载作用下，假设排土场基底内部土层界面完全连续，即$ {\sigma _{zi}} = {\sigma _{zi + 1}} $，$ {\tau _{zi}} = {\tau _{zi + 1}} $，$ {u_{ri}} = {u_{ri + 1}} $，$ {u_{zi}} = {u_{zi + 1}} $。设有N层基底土，第N层基底土视作半无限空间体，则上述式中i变为N−1，即：

式中：${z_i} = \displaystyle\sum_{j = 1}^i {h_j}$，$ {m_k} = \dfrac{{1 + {\mu _i}}}{{1 + {\mu _{i + 1}}}}\dfrac{{{E_{i + 1}}}}{{{E_i}}} $。N层弹性体系有N−1个层间界面，每个界面有4个层间结合条件，再加上2个表面应力边界条件(式9)和层间连续条件A_N=C_N=0，一共可建立4N个线性代数方程式，这样就可以求得系数A_i、B_i、C_i和D_i，(i=1，2，…，N−1)，排土场基底通常可以简化看成单层和多层情况，以下为具体求解过程。

基底为单层(N=1)时，即为弹性半空间体。由基底内部土层界面连续条件得C_N=C₁=0，D_N=D₁=0，由式(9)得A₁=−2μ₁，B₁=−1。圆形均布垂直总荷载为πδ²q，当圆心为0，作用点在z轴时，零阶汉克尔变换式为$ \overline q \left( \xi \right) = q\dfrac{\delta }{\xi }{J_1}\left( {\xi \delta } \right) $。设一垂直集中力Q与圆形均布垂直荷载相等，即Q=πδ²q，则当荷载圆半径δ→0并求极限，垂直集中力Q的零阶汉克尔变换式为$ \overline q \left( \xi \right) = \mathop {\lim }\limits_{\delta \to 0} \dfrac{{Q{J_1}\left( {\xi \delta } \right)}}{{\pi \xi \delta }} = \dfrac{Q}{{2\pi }} $。将排土场沿坡面竖向条分，各微分段dxdy均可视作作用于地表面的垂直集中力，以此积分叠加，所产生的地表位移u_z为

式中：h为边坡高度，m；α为边坡角，(°)；L为排土场坡顶平盘宽度，m；D为排土场走向长度，m；a为坡脚距地表A点距离，m；γ为排弃物料容重，kN/m³。

当基底为多层时，由式(10)层间接触面位移和应力连续条件以及式(9)地表边界条件构建的方程组解得地表弹性变形为

1.2   基底土体黏弹性本构模型选取

排土场黏土基底蠕变行为表现出时变收敛的特征，因此引入了描述土材料时变蠕变特性的广义开尔文模型。黏弹性模型通常由若干个胡克体H和黏滞体N组合而成。为准确描述排土场基底的蠕变特性，笔者选用如图2中2个Kelvin模型串联而成的广义Kelvin模型。

图  2  广义Kelvin模型元件

Figure  2.  Element diagram of generalized Kelvin model

下载: 全尺寸图片 幻灯片

Kelvin模型本构方程：

五单元广义Kelvin模型式：

蠕变柔量为

式中：Δσ=σ₁−σ₃，σ₁为轴向应力，MPa，σ₃为围压，MPa；I₁为直角坐标系下第一应力不变量，MPa；K₁，K₂为广义Kelvin模型体积模量，MPa；G^H为广义Kelvin模型中胡克体剪切模量，MPa；$G_1^{\mathrm{K}} $、$G_2^{\mathrm{K}} $为广义Kelvin模型中Kelvin体剪切模量，MPa；η₁、η₂为广义Kelvin模型中Kelvin体黏滞系数，MPa·h。

广义Kelvin模型的剪切模量G为蠕变柔量的倒数，G=1/J(t)。

1.3   均质基底排土场周边地表变形黏弹性解法

当排土场基底N=1时，可视为半无限均质基底。土体变形参数泊松比、弹性模量^[27]分别为

若在均质基底排土场上作用集中荷载Q，对柱坐标系下Boussinesq位移解做Laplace变换得：

式中：s为复变量，S⁻¹。

对广义Kelvin模型，应变偏量张量e_ij为

式中：G^K为广义开尔文模型开尔文体的剪切模量，MPa；G^H为广义开尔文模型中胡克体的剪切模量，MPa。

通过蠕变柔量分量做Laplace变换得到体积模量$\hat K$和剪切模量$\hat G$：

则集中力Q引起的地表变形：

对式(21)做Laplace逆变换，得到黏弹性均质基底表面作用一集中力Q的变形解答，将排土场荷载微分，叠加微分荷载对地表变形的作用，可求得黏弹性均质基底地表变形表达式。

1.4   层状基底排土场周边地表变形黏弹性解法

将基底的柔度系数V(t)代替弹性模量E表征黏弹性体的刚度特性，V(t)=1/E(t)。对于广义Kelvin模型的本构方程式(13)做Laplace变换和逆变换，获得柔度系数：

由弹黏性材料的应力与应变关系，设G(t−ζ)=(1+μ)V(t−ζ)^[27]，V(t−ζ)为t−ζ时刻的柔度系数，ζ为时间变量，得到任意时刻t的地表变形表达式：

由基底表面边界条件及层间连续条件可得

令G_i(t−ζ)=$ G_i^E $(t−ζ)+$ G_i^v $(t−ζ)，其中上标“E”表示弹性，上标“v”表示黏弹性。将含$ G_i^v $(t−ζ)的积分部分用梯形积分式表达。取两积分节点间足够小增量dT^[28]，dT=t/(t−1)，m为节点数，于是：

其中，

通过积分回代解出积分参量A_i、B_i、C_i、D_i，从而得到地表变形表达式(23)。

2.   工程应用

笔者以歪头山铁矿外排土场作为工程实例，将下盘排土场坡顶坡底视为矩形的梯形台，其东北侧临近高铁线，最小距离不足50 m，经测量下盘排土场与高铁走向的并行长度约1300 m，排土场在2020—2021年4月进行排土工程作业，边坡高度由128.31 m增高至164.67 m，边坡角提高4°～14°。在矿区地势近水平区域，与高铁并行范围内布设有1条GNSS监测线XP−3(监测点编号为X1～X3)，其中X1监测点位于排土场坡脚处(a=0 m)，X2监测点位于G304高速公路临近排土场一侧(a=20 m)，X3监测点位于高铁线墩台附近(a=50 m)，监测点于2021年4月开展沉降监测。以监测频率1次/周监测到2021年7月底，之后监测频率为1次/月。现场平面图以及监测线剖面图如图3所示。以下盘排土场基底域内中心为原点O，西北帮与高铁线并行方向为y轴，沉降值u_z、相邻墩台沉降差∆u_z为z轴，建立坐标系，工程示例如图4所示。

图  3  下盘排土场平面图及监测线剖面图

Figure  3.  Footwall dump plan and monitoring line profile

下载: 全尺寸图片 幻灯片

图  4  排土场与高铁线墩台相对位置示意

Figure  4.  Settlement diagram of high-speed railway pier under the action of dump

下载: 全尺寸图片 幻灯片

2.1   现状排土场周边地表变形解析

矿区排土场及基底的岩土体物理力学参数见表1。运用本文层状基底排土场周边地表变形计算方法对排土场现状形态进行数值计算，判断排土场荷载对周边高铁及墩台影响是否达到相关规范^[29-30]标准(规范中，高铁属于无砟轨道，其绝对沉降值不应超过20 mm，相邻墩台沉降差不应超过5 mm，高架桥相邻墩台间隔为25 m，即两墩台之间沉降差不应超过0.2 mm/m)。

表  1  地表变形分析的土体物理力学参数

Table  1.  Soil physical and mechanical parameters of surface deformation analysis

岩土	容重 γ/(kN·m−3)	弹性模量 E/MPa	泊松比 μ	K/ MPa	GH/ MPa	$G_1^{\mathrm{K}} $/ MPa	η1/ (MPa·h)	$G_2^{\mathrm{K}} $/ MPa	η2/ (MPa·h)
排弃物	18.6	500.0	0.200	—	—	—	—	—	—
基底黏土	16.5	46.3	0.258	1.02	3.12	0.03	2.70	0.06	0.47
花岗岩	26.1	644000.0	0.180	—	—	—	—	—	—

下载: 导出CSV 

| 显示表格

歪头山铁矿下盘排土场与高铁相距a=50 m，排土场基底为黏土层和混合花岗岩层组成的双层黏弹性体(N=2)，黏土层厚度H₁约5 m，混合花岗岩层可视作半无限弹性空间体。当N=2时，C_N=C₂=0，D_N=D₂=0，运用层状基底排土场周边地表变形黏弹性解析式(23)可求解出高铁及墩台受排土场影响的变形值。图5为高铁沉降u_z与时间t、距坡脚距离a的关系曲线，图6为高铁沉降值u_z、相邻墩台沉降差∆u_z与墩台位置关系曲线。

图  5  地表变形与时间和距坡脚距离的关系

Figure  5.  Relation between surface deformation and time and distance from slope foot

下载: 全尺寸图片 幻灯片

图  6  高铁沉降值、相邻墩台沉降差与墩台位置关系

Figure  6.  High railway settlement value, adjacent pier settlement difference and pier position relationship

下载: 全尺寸图片 幻灯片

可以看出排土场周边地表变形随着时间t演化逐渐增大后趋于稳定，同时距坡脚距离增大而减小；高铁线墩台沉降u_z和相邻墩台沉降差∆u_z会随着时间演化逐渐增大，后趋于平缓；高铁在−650～+650 m段沉降u_z=2.08 mm，高铁离开排土场并行范围(即≤−650 m和≥+650 m)，在±775 m时最大相邻墩台沉降差∆u_zmax=0.006 mm/m。因此，歪头山铁矿下盘现状排土场周边地表变形与高铁相邻墩台沉降差均低于规范阈值，处于安全稳定状态。

2.2   地表变形解析解与监测系统数据对比

由XP−3监测线上X1～X3监测点的沉降曲线(图7)，可以看出随着排土场建设完成后，3个监测点的变形值均不断增加，由于工程荷载非线性增加，变形增加无明显规律，直至2021年6月底，排土场周边地表沉降基本停止。运用本文地表变形解析方法，获得3个监测点位置地表沉降随时间的演化曲线，并加入监测数据如图8所示。可以看出，排土场周边地表沉降随时间演化而增长，在1～2个月内达到峰值后趋于稳定，对比监测数据，监测点位置的地表沉降预测值与监测数据基本吻合，误差低于1%。

图  7  X1～X3监测点监测结果

Figure  7.  Monitoring results of monitoring points X1−X3

下载: 全尺寸图片 幻灯片

图  8  X1～X3监测点位置地表变形演化预测

Figure  8.  Prediction of surface deformation evolution at monitoring points X1−X3

下载: 全尺寸图片 幻灯片

随着排土场工程的发展，注意到基底土体在短期加载下会表现为近似瞬时的弹性变形，然后逐渐进入蠕变阶段，直至稳定在工程建设阶段性结束后的一段时间内。目前，由监测系统表明现有排土场周边地表土体变形满足相关规范的要求。

然而，考虑到实际现场条件的多样性，不仅要满足规范要求，还要对排土场进行增高和扩容的设计。因此，为了更好地指导排土场形态的设计，本文运用黏弹性解法，对排土场周边地表变形进行了预测和分析。这种方法更全面地了解了排土场周边地表变形的时空演化，以及边坡形态等重要参数对这些变形的影响。这不仅有助于设计更安全的排土场，还可以提供实际工程建设的重要指导。

2.3   排土场周边地表变形预测与时空演化规律分析

在选定下盘排土场的合理边坡高度为120～180 m和边坡角度为10°～20°后，针对不同的边坡高度和边坡角度，通过图9和图10的观察，发现地表变形与距离呈现出反比例关系。此外，随着边坡高度和角度的增加，地表变形也显著增加。

图  9  不同边坡高度排土场作用下周边地表变形规律

Figure  9.  Surrounding surface deformation law of dump with a different slope height

下载: 全尺寸图片 幻灯片

图  10  不同边坡角排土场作用下周边地表变形规律

Figure  10.  Surrounding surface deformation law of dump with different slope angles

下载: 全尺寸图片 幻灯片

另外，在不同的边坡高度和边坡角度下，排土场的作用导致黏弹性基底地表变形，如图11所示，可以观察到随着时间的推移，地表变形呈现出持续增加的趋势。在经过1～2个月后，土体变形趋于稳定。当排土场高度达到333 m，且边坡角度为20°时，土体变形值接近规范阈值，表明不宜继续增加排土场的高度。

图  11  黏弹性基底排土场周边地表变形随时间变化规律

Figure  11.  Surface deformation around dump with viscoelastic basement changes with time

下载: 全尺寸图片 幻灯片

从时空角度考察，随着排土场工程的不断进行，地表变形随时间的演化表现出多种趋势，并最终趋于稳定。此外，不同基底高度和不同边坡角度对地表变形产生不同的影响，这些影响在空间上表现出明显的特征，即距离坡脚越远，地表变形越小，而边坡高度和角度越大，地表变形越显著。因此，这些研究结果为实际工程提供了重要的指导，建议选择较小的边坡高度和角度以减小地表变形，同时需谨慎控制排土场的作用时间，以避免超过土体的承载能力，确保工程的安全和稳定。

3.   结　　论

1) 借助层状弹性理论，将微分的排土场载荷对周边土体的作用进行叠加，引入广义Kelvin三轴蠕变模型，推导了具有流变特性的黏弹性软弱基底排土场周边地表变形解析方法，为研究排土场载荷作用下周边地表变形时空演化规律奠定理论基础。

2) 将黏弹性基底排土场周边地表变形解析方法应用于歪头山铁矿下盘排土场工程实际，获得的周边地表变形预测解析结果与现场监测数据基本吻合，验证了方法的合理性；随着距坡脚距离增大，排土场周边地表变形逐渐减小；地表变形与相邻墩台沉降差均随着时间的推移不断增大，在短时间达到峰值，之后趋于稳定。

3) 在加载阶段，随着距离坡脚的增加，排土场周边地表变形的变化趋势可以用反比例曲线来描述，即变形随距离的增加而减小，但下降速率逐渐减缓，并且当边坡高度或边坡角增大时，地表变形会随之增大；当排土场压实稳固进入蠕变阶段，周边地表变形会在1～2个月内达到峰值，之后趋于稳定；运用本文预测方法，可指导水平层状基底排土场形态设计，保障周边设施安全。