山西沉积型地下铝土矿体开采方式通常为房柱法^[1].为了防止上覆硬质黏土层风化冒落，矿山通常采取预留一定厚度的铝土矿护顶层、选择合理空区结构、保留合适尺寸的矿柱等措施，形成协同支撑系统，共同维持采空区的稳定^[2].当系统内部含有结构弱面、空区尺寸选择不合理、积蓄能量超过承载极限阈值时，将导致矿柱发生滑移破坏，进而造成采空区突然坍塌.因此，开展沉积型铝土矿支撑系统滑移突变失稳机理研究，具有重要的理论和工程意义.

目前，国内外已有采空区支撑系统稳定性研究的相关文献报道.Qin等^[3]通过构建顶板-煤柱力学系统失稳演化模型，给出了采空区失稳判据及位移突跳计算式.姜立春等^[4]在构建矿柱-顶板受力模型的基础上，研究了某铜矿采空区的结构参数与稳定性之间的关系，提出了矿柱安全系数的计算方法.Wang等^[5]在构建由多个矿柱构成的矿柱体系受力模型的基础上，研究了单轴抗压强度对串行和并行矿柱体系稳定性的影响.Singh等^[6]在研究煤矿开采深度与垮落带内破裂岩体性质基础上，提出了矿柱极限强度的计算式.Wagner等^[7]运用数值模拟的方法，分析了采空区参数对矿柱支撑系统的荷载作用，研究了局部矿柱破坏对支撑系统稳定性的影响.综上，现有文献鲜有涉及铝土矿床采空区，少有考虑裂隙或优势节理面等结构弱面的影响，忽略弱面控制型矿柱失稳破坏问题; 仅考虑矿柱单独变形，忽视顶板与矿柱的协同变形作用.在实际工程中，矿柱、顶板为一个整体协同构成采空区支撑系统，共同维系采空区稳定.对沉积型铝土矿体而言，岩层内部普遍存在大量结构弱面，矿柱变形破坏几乎均从斜向节理裂隙扩展开始，最终发生弱面控制型破坏.

本文在山西沉积型铝土矿采空区调研的基础上，综合此类采空区破坏特征，通过对矿柱-顶板受力分析，构建采空区矿柱-顶板支撑系统滑移突变模型，研究含弱面支撑系统的破坏失稳机制，分析系统内控参数对采空区稳定性的影响，利用工程实录数据验证计算结果的合理性和可靠性.

1 含弱面矿柱-顶板支撑系统力学模型

以山西晨光地下铝土矿典型采空区为对象(图 1).空间结构示意图(图 1b)中，矿柱宽度为a，采空区跨度为b，矿柱单侧屈服区宽度为c，上覆岩土层厚度为H.

由于采空区的径向长度远大于采场跨度，其变形问题可按弹性平面问题处理.设定采空区底板不发生变形，矿柱底部不产生线向位移和角位移.将上覆岩层和顶板自重作用简化为均匀荷载q作用，空区顶板视为刚度为EI的弹性梁，矿柱-顶板发生协同变形时，采空区沿矿柱(主优势节理、裂隙等)结构弱面滑移破坏，顶板岩层不发生破坏仅参与释放积聚能量.

由于采空区为对称分布，可将矿柱及两侧空区中心区域视为单元结构体，建立含弱面矿柱-顶板滑移失稳的力学结构简化模型(图 2).图 2中，在均匀荷载q作用下，顶板将发生挠曲变形，内部积聚弹性形变能，同时将部分能量传递给矿柱.随着矿柱累积能量的不断增加，系统在失稳破坏之前，将沿着矿柱内部某一结构弱面(优势节理面、裂隙等)产生一定滑移变形.

当矿柱几何尺寸选择和空间布置合理时，支撑系统内部能量释放与其结构变形相互协调，系统将保持稳定; 当矿柱的几何尺寸选择或空间布置不合理时，系统内部累积势能将转化为动能瞬间释放，支撑系统将沿结构弱面突然滑移、失去平衡，进而引起采空区坍塌.

根据A.H.Wilson两区约束理论^[8]，在上覆岩层荷载q长期作用下，矿柱两侧区域的岩体应力已超过屈服点，将产生应力屈服区(图 2).此时，矿柱屈服区内岩体介质强度降低，表现为明显的损伤弱化特征，可采用负指数模型来表示其本构关系^[9]:

(1)

式中：τ_p为结构弱面介质在屈服区段的剪应力; G_p为初始剪切模量; ε为结构弱面介质的应变; ε₀为峰值剪应力对应的应变.

矿柱中心区域岩体相对变形较小，为弹性核区(图 2)，该区域岩体介质应力没有超过屈服点，其应力-应变关系为

(2)

式中:τ_e为结构弱面介质在弹性核区段的剪应力; G_e为剪切模量.

2 支撑系统滑移突变模型 2.1 顶板挠曲线方程

建立矿柱-顶板支撑系统变形的y-x坐标系(图 2).假设在上部均匀荷载q作用下，系统沿结构弱面向下产生滑移位移u.根据顶板变形边界条件，设顶板的挠度曲线方程为^[9]

(3)

其中：A为顶板的挠度系数; θ为结构弱面的倾角.

顶板与矿柱应满足协调变形条件：

(4)

将式(3)代入式(4)，求解可得

(5)

令，可得顶板的挠曲线方程为

(6)

2.2 支撑系统的势能函数Π

在含弱面矿柱-顶板支撑系统中，顶板因挠曲变形将累积弹性应变能U_q; 矿柱受到压应力作用，将在结构弱面附近形成微裂纹富集软弱带，软弱带内部的微裂纹扩展、连通，将耗散能量U_s(屈服段应变能U_p和弹性段的弹性势能U_e)，宏观上表现为矿柱产生滑移位移u.顶板弹性应变能和矿柱耗散能量均来源于荷载q做功所转化的能量.因此，支撑系统的总势能为

(7)

式中：U_q为顶板积聚的挠曲应变能; U_p为矿柱结构弱面在屈服段内的应变能; U_e为矿柱结构弱面在弹性核区段内的弹性势能; W为荷载q所做的功.

其中，

(8)

(9)

(10)

(11)

式中：EI为弹性梁的刚度; B为软弱带的厚度; u₀为峰值荷载对应的位移.

将式(8)~式(11)代入式(7)，可得系统的总势能：

(12)

2.3 支撑系统滑移突变模型参数

研究表明^[3]，常见的基本突变模型有7种，其中尖点突变模型在采矿工程中应用最广，其平衡曲面标准形式方程为

(13)

式中:x为系统的状态变量; m, n为平衡曲面的控制变量.

由式(13)绘制流形曲面图如图 3所示.平衡曲面共有3叶:下叶表示支撑系统滑移失稳的孕育过程，此时系统内部势能不断累积，为准稳定状态; 中叶表示系统不稳定状态; 上叶表示系统失稳能量释放后的新稳定状态.

随着控制参数m, n的连续变化，如果系统状态沿路径Ⅰ演化，状态变量x变化连续，支撑系统将不发生突变; 如果系统状态沿路径Ⅱ演化，当接近折叠翼边缘(图 3圆点处)时，控制参数m, n仅发生微小变化，支撑系统状态也将从折叠翼的下叶跃迁到折叠翼上叶，即支撑系统将发生突变.

针对沉积型铝土矿采空区含弱面支撑系统滑移突变失稳的特点，根据系统总势能函数表达式(12)，利用数学方法将其转化为式(13)的标准形式，即可确定m, n两个控制参数.

取矿柱的滑移位移u为状态变量，对总势能П求一阶导数，令其等于0，可得含弱面矿柱-顶板支撑系统稳定性分析的平衡曲面方程：

(14)

根据平衡曲面的光滑性质，在尖点处应有

(15)

求解式(15)，可得尖点处应有u=u^*=2u₀.

将平衡曲面方程(14)在尖点处按Taylor级数展开，截除三次幂以上高阶项，经简化可得

(16)

令状态变量x=(u-u^*)/u^*，结合式(13)可得

(17)

由式(17)可知，控制参数m主要与采空区跨度b、矿柱宽度a、结构弱面倾角θ、软弱带厚度B、屈服区宽度c等采空区内控参数密切相关; 控制参数n除与上述采空区内控参数相关外，还受外控参数采空区上覆荷载q的影响.

3 支撑系统滑移失稳机制 3.1 充分必要条件

从前面分析可知，支撑系统滑移失稳的过程实质表现为系统状态从下叶到上叶的跃迁过程，仅当系统状态跨越折叠翼边缘时，支撑系统才可能发生滑移失稳.平衡曲面的折叠翼在m-n平面上的投影(图 3)称为支撑系统的分叉集，其代数式为

(18)

当D>0时，系统处于稳定状态; D=0为系统失稳的临界状态; D < 0时，系统状态已跨过折叠翼边缘，处于失稳后的新稳定状态.

1) 充分条件.若D≤0成立，有：当n>0，即在分叉集的右支区域，状态变量x没有发生跃升，此时支撑系统仅在数学结构上发生突变; 当n < 0，即在分叉集的左支区域，支撑系统将发生滑移突变失稳.

因此，支撑系统滑移失稳的充分条件为D=0，且满足n < 0.

2) 必要条件.若D≤0成立，须满足m≤0，即为系统滑移失稳的必要条件.

结合前面分析可知，在采空区岩体组成既定条件下，支撑系统滑移失稳的必要条件仅与采空区内控参数密切相关，与外控参数无关.

3) 充分必要条件.综合1)，2)可知，支撑系统滑移突变失稳的充要条件为

(19)

综合分析可知，支撑系统发生滑移失稳破坏主要是由内控因素与外控因素协同作用的结果.内控因素包括采空区跨度b、矿柱宽度a、结构弱面倾角θ、软弱带厚度B、屈服区宽度c等，外控因素主要为采空区上覆荷载q作用.

3.2 临界位移与能量释放量

为预判采空区滑移失稳及后果，对系统发生滑移失稳的临界位移及能量释放量进行计算.当支撑系统满足滑移突变失稳的充分必要条件时，有

1) 若m=0，式(13)有三重零根;

2) 若m < 0，式(13)有3个实根，分别为x₁= .

支撑系统跨越分叉集时，状态变量x的移动轨迹如图 4所示.支撑系统跨越分叉集前的滑移变形突跳量为

(20)

此为支撑系统滑移失稳时的临界位移，可作为预判失稳时的位移判据.

系统失稳前后的变形突跳量为

(21)

图 4中M，N两点的势能之差，可以用来估算矿柱-顶板支撑体系在突发失稳过程中释放的能量.将势能函数方程式(12)在尖点u^*处Taylor展开，截除四次幂以上高阶项，引入无量纲参数式(13)，可得系统突发失稳过程中释放的能量为

(22)

式中:参数m, n的含义与式(17)相同.

3.3 控制因素影响分析

对某矿山具体采空区而言，采空区的上覆围岩高度既定，围岩的构成和岩性既定，因此，在分析控制因素对支撑系统的影响时，仅需要分析内控参数的影响.以山西晨光地下铝土矿山为例，进行各内控参数对采空区稳定性影响的计算分析.

1) 控制参数与内控因素的关系.由失稳充分必要条件式(19)计算可得控制参数n与各内控参数之间曲线关系，结果如图 5所示.

由图 5a可知，在其他外控因素不变情况下，控制参数n与采空区跨度b呈负相关，随b增大而减小，但减小的速度逐渐减缓(图 5a).这是因为b越大，矿柱所受荷载越大，支撑系统越容易发生滑移失稳.

由图 5b可知，控制参数n与矿柱宽度a呈正相关，随a增大而增大，且增大的速度逐渐增快(图 5b).这是a增大增强了矿柱承载能力的结果.

由图 5c可知，控制参数n与结构弱面倾角θ呈负相关，随θ增大而减小，减小的速度逐渐减缓(图 5c).究其原因，θ越大，上覆荷载沿结构弱面的分力越大，越易诱发支撑系统滑移失稳.

由图 5d可知，控制参数n与软弱带厚度B呈负相关，随B增大而减小(图 5d).随着B增大，矿柱沿结构弱面的抗剪能力减小，支撑系统自稳能力变弱.

由图 5e可知，控制参数n与屈服区宽度c呈负相关，随c增大而减小，但减小的速度逐渐减缓(图 5e).这是由于c增大，矿柱损伤弱化区域变大，矿柱强度变弱，导致系统失稳.

综上可见，随着采空区内控因素b增大、a减小、θ增大、B增大、c增大，将更易满足支撑系统滑移突变失稳的充分必要条件，引发采空区失稳.

另外，由图 5分析可知，随着b增大、a减小、c增大，n减小速度逐渐减缓，这是采空区空间尺寸效应导致的; 随着θ增大，n减小速度逐渐减缓，归因于倾角效应的作用.

同理，分析控制参数m与各内控参数之间关系，可得相似结论，这里不再展开论述.

2) 内控因素敏感度.由图 5还可以发现，各参数在可取区间内变化时，均与n呈现出似线性关系.b, a, θ, B及c等内控参数的线性拟合系数分别为R_a=0.995 6, R_b=0.995 9, R_θ=0.991 6, R_B=0.998 8, R_c=0.985 1.

参考文献[10]的敏感度计算方法，可得控制参数n与a, b, θ, B及c等5个内控参数的敏感度，大小分别为f_a=3.56, f_b=4.02, f_θ=2.60, f_B=2.92, f_c=0.50(图 6).据此可知，各内控因素的敏感度大小依次为f_b>f_a>f_B>f_θ>f_c.

同理，计算可得各内控因素对控制参数m的敏感度分别为f_a=2.79, f_b=3.83, f_θ=1.90, f_B=2.26, f_c=0.71.各内控因素的敏感度大小依次为：f_b>f_a>f_B>f_θ>f_c，结果与控制参数n的相同.

根据上述计算结果可知，软弱带厚度B和结构弱面倾角θ的敏感度仅略小于采空区跨度b、矿柱宽度a，因此，矿山在开采设计过程时，应密切关注结构弱面相关因素的影响.

4 工程验证

山西晨光铝土矿矿区为典型的黄土高原地貌.矿区岩土层自上至下依次为第四系中上更新统(Q₂₊₃)、上第三系上新统保德组(N₂)、石炭系上统太原组(C₃t)、石炭系中统本溪组(C₂b)、奥陶系中统(O₂f).其中石炭系中统本溪组为含铝土矿层，铝土矿体的分布与形态严格受奥陶系侵蚀面的控制，矿体形态呈层状、似层状、透镜状产出，倾角5°~12°，一般不大于10°(图 1b).矿山主要采用房柱法开采，经过多年累积开采，遗留下大量采空区群.空区直接顶板为预留铝土矿体护顶岩层，平均厚度小于0.5 m.

拟研究区域见图 7，采空区结构参数及岩土体物理力学参数如表 1所示^[11].

根据现场调研结果，取滑移面倾角θ=42°，软弱带厚度B=0.3 m.由支撑系统发生滑移突变失稳的必要条件m≤0，计算可得：m=-0.816 < 0，可以判断采空区内部结构参数符合失稳的必要条件.

由系统失稳的充分必要条件式(19)，计算可得n=-0.271 < 0，此时有D=-1.283 < 0，可以判定该采空区已经达到了失稳的充要条件，支撑系统处于不稳定状态.

现场调查表明，2018年7~9月间，该区域内多次发生采空区垮塌事故(图 8)，矿柱多处发生局部滑移崩塌现象(图 8a)，与理论分析结果相符，验证了理论计算结果的合理性.

5 结论

1) 针对沉积型铝土矿采空区破坏主要表现为矿柱沿剪切弱面滑移破坏的特征，借助能量耗散理论和突变理论，构建了矿柱-顶板支撑系统滑移突变失稳模型，研究受结构弱面影响的支撑系统失稳机制.

2) 采空区支撑系统的突发滑移失稳主要是由于上覆荷载作用等外控因素和采空区跨度b、矿柱宽度a、结构弱面倾角θ、软弱带厚度B、屈服区宽度c等内控因素协同作用的结果.在外部荷载因素q与岩性既定的条件下，各内控因素的敏感度大小依次为：f_b>f_a>f_B>f_θ>f_c.其中f_B, f_θ略小于f_b, f_a.在实际工程中，应密切关注结构弱面因素的影响.

3) 计算结果表明，山西晨光铝土矿采空区结构参数已经达到了失稳的充要条件，支撑系统处于不稳定状态，现场实录数据验证了理论计算结果的合理性.