0 引 言
山区地表采动损坏带来的最严重问题是引发了山体滑坡，许多学者开始研究对于开采沉陷引起的采动滑坡，揭示矿区采动滑坡的特性，分析采动滑坡的发生机理等。汤伏金等[1]探讨了地下开采诱发滑坡的基本规律，分析了采动滑坡特点，并对采动滑坡的成因进行分类，对采动滑坡的物理模拟进行实验研究。梁明等[2]将采动山体变形划分为3个阶段并提出相应的变形监测方案，讨论了采动山体滑坡监测网的设计原则及观测方法，提出了采动山体变形监测数据分析和采动引起滑坡的预报方法。汤伏全等[3]着重分析了滑动型失稳分裂面的形成特征，并导出其稳定性计算的解析式，并以陕西韩城象山采动滑坡为模型，通过相似材料模拟实验和有限元数值分析，对采动诱发山体顺层(面)滑坡的机制进行了研究。汤伏全等[4]通过对实际滑坡资料的分析研究，总结出采动滑坡的特点和规律，运用山区地表移动的理论，分析采动坡体时上覆岩层各部位应力、应变的变化规律，考虑岩石的不抗拉特性，用弹塑性有限元数值模拟对采动坡体的稳定性进行分析，并详细讨论建模时应解决的相关问题，并在此基础上提出采动滑坡贯通性滑动面的确定和采动坡体稳定性的计算方法。陆力等[5]就平朔煤炭工业公司安太堡不采区B900工作面开采地表沉陷实测研究特别是采用GPS动态观测研究的基础上，深入分析了随综放采场推进边坡地表沉陷发展的规律，包括地表采动边坡沉陷范围、最大沉陷值、随采场推进而趋于稳定的时间以及边坡破坏的主要形式，提出了应采取的控制开采方法和安全防范措施[6]，提出构建基于GIS的山区采动滑坡预报系统，从而实现矿区绿色开采及灾害及早防治[7]。
但是，对于采动条件下坡体的变形破坏问题以往只是侧重于对地表的移动和地面沉降的研究，而对坡体自身工程地质条件与采空区塌陷引发地表移动的相关性、协调性涉及相对较少。因而，准确的研究采动滑坡灾害的形成机理以及矿山开采系统是减少矿山滑坡灾害，合理利用矿产资源，改善矿山环境的关键，此研究具有重大的经济价值和社会意义。
基于此，笔者以山西省某采空区为例，对20世纪开采3煤形成的采空区与近10年开采8煤形成的采空区进行模拟，通过三维数值模拟[8]，对煤层开采条件下上覆岩层的位移、应变等进行分析，从定性、定量方面对天然采空工况下采动坡体的稳定性进行评价。
1 研究区工程地质条件
结合采空区地质资料及现场勘察，采动坡体位于太行山北段西侧，沟谷纵横，属中低山地貌，微地貌发育。区域地貌受地质构造、新构造运动及地层岩性的制约，按其形态及成因类型可以分为构造溶蚀侵蚀低中山区、构造剥蚀低中山区、构造侵蚀煤系地层低山丘陵和侵蚀堆积谷地。工作区西部及北部为构造剥蚀低中山区，标高为820～1 020 m，切割深度约200 m，坡脚较缓。根据野外地质调查并结合相关资料，区域内赋存有奥陶系、石炭系、二叠系及第四系地层。区内主要发育一组褶皱构造，形态以开阔褶皱为主，工作区处于向斜的东南翼。该向斜呈北东展布，轴长7 500 m，组成地层为二叠系，北西翼倾向南东，倾角7°～13°，东南翼倾向北西—南西，倾角3°～10°。区内断层不发育。
研究区内主要的含煤地层为石炭系上统太原组和二叠系下统山西组。太原组平均厚度125.50 m，含煤8层，煤层平均厚度16.88 m，含煤系数13.45%；山西组平均厚度56.25 m，含煤6层，煤层平均厚度4.18 m，含煤系数7.43%。其中，可采煤层主要有3煤、8煤、9煤、15煤，开采标高+550～+390 m。区内发育有5套含水层：奥陶系石灰岩岩溶裂隙含水层、石炭系上统太原组碎屑岩夹石灰岩岩溶裂隙含水层、二叠系下统山西组碎屑岩裂隙含水层、二叠系石盒子组碎屑岩裂隙含水层以及第四系全新统砂砾层孔隙含水层。各含水层之间有较多优良的隔水层段，使各含水层之间呈近平行的储水系统，起主要作用的隔水层段有：石炭系中统本溪组黏土泥岩、泥岩隔水层和石炭系上统太原组及二叠系下统山西组泥岩隔水层。区域构造不发育，含水层之间都存在着较稳定的相对隔水层，所以各含水层之间水利联系较差。
斜坡体地貌单元属侵蚀、剥蚀中低山，坡率平均约30°，上部坡体坡度为30～450°。因砂岩、泥谷长期的操露，风化利蚀严重，在排除人为因素的影响下，斜坡的安全系数也是较低的，是处于基本稳定状态。斜坡体邻近地下煤矿开采形成的采空区，位于该斜坡中上部的正下方。采空区上覆岩层在沉陷变形过程中当处于斜坡段时多伴随有水平或侧向位移，特别是在强风化破碎岩层和或坡积堆积体中尤为明显，这种侧向或水平变形产生的推力容易使下侧黄土形成挤压、鼓胀和隆起等不良形变。
2 计算原理
Midas GTS NX的计算原理是通过划分好的网格进行节点的坐标值进行计算来确定岩土体的应力应变值，相对初始状态进行模拟分析，对初始应力场进行模拟计算。得到结果以后通过设置一个安全系数，将岩土体的参数不断进行折减，临界折减系数就是安全系数，直到最终计算出来得到的结果不收敛。并根据计算出来的应力应变值绘制出应力应变云图,最终找到破坏点和破坏面。
岩土体参数的折减方法依据强度折减法，又称为强度储备法，该法认为材料由于外界环境因素的影响，其抗剪强度指标无法确定，而且波动比较大。因此，在稳定性分析中、在弹塑性有限元分析中，将岩土体抗剪强度参数黏聚力和内摩擦角同时除以折减系数 ，得到一组新的黏聚力和内摩擦角。计算公式如下：
cy=c′/SRF
(1)
tanφy=tanφ′/SRF
(2)
式中：SRF为有限元强度折减系数；c′、cy为折减前后的有效黏聚力；φ′、φy为折减前后的有效内摩擦角。
模型需要遵守摩尔-库仑准则，即岩土体的黏聚力以及内摩擦角是决定其剪切破坏大小的因素，当研究对象的剪应力超出岩土体的黏聚力，以及内摩擦力所控制的抗剪应力时，则对岩土体产生剪切破坏[9]。具体公式为
(3)
式中：Fs剪应力；σ1为最大主应力；σ3为最小主应力；C为岩土体的黏聚力；ø为岩土体内摩擦角。
3 网格划分及模型确定
根据前期的勘察资料，采空区所在山体的顶部以及中下部不同程度地发育有裂缝，因此，这座山体被确定为研究重点。依据山体地形特征，数值模拟计算范围确定为与山体两侧近于平行的长方形，面积约0.59 km2。网格划分的正确与否直接决定着后续计算分析的结果，在对模型进行网格划分时，不仅要考虑其地层、地形、网格规格、划分方式等因素，而且对于模型中采空区处的网格需要加密处理，从而让整个模型的地质力学特性，能够更详细的体现出来。在网格划分完以后，需要对整个模型的网格质量进行检查，本次模型的地层单元网格长度设置为4 m，模型共划分了27 005个节点，145 884个单元。
根据前期资料，模型岩体主要包括：砂泥岩、煤层、黄土层、第四纪冲积体。
图1 几何模型及网格划分示意
Fig.1 Geometric model and meshing diagram
4 模拟方案及物理力学参数
模型边界条件采用X轴方向侧面约束X方向位移，Y轴方向侧面约束Y方向位移，底面约束Z方向位移，顶部为自由边界。本次模拟选用摩尔-库仑模型(Mohr-Coulomb)和空模型分别模拟岩体和采空区，模拟工况为天然采空工况状态。
根据资料，山西某矿3煤为房柱式开采，分别于1982—1985年、1986—1992年、2004—2006年开采了与山体相邻的北部、中部和南部，目前3煤已经采完，矿柱已全部采出。8煤为壁式开采，目前8101、8102、8103工作面已采空。本次模拟简化开采过程，3煤和8煤均为3次开采，分北部、中部和南部3部分。黄土覆盖层的参数采用来自山西省交通科学院的测试结果。对于粗砂岩、细砂岩、粉砂岩、砂泥岩、煤层、第四纪冲积体、强风化带，物理力学参数的选取主要来自同类岩体的横向比较。具体数值见表1。
表1 天然采空工况参数
Table 1 Natural goaf parameters
岩类弹性模量/GPa泊松比黏聚力/MPa内摩擦角/(°)抗拉强度/kPa容重/(kg·m-3)强风化带3.120.310.80032.0780.02 290Q2+30.200.320.03222.41.01 420Q40.200.350.03020.02.01 740P2s3砂泥岩12.500.322.00036.51090.02 390P2s2、P2s1砂泥岩10.500.322.00030.0960.02 390P1x3、P1x2砂泥岩9.200.301.80037.1820.02 390P1s、P1x1砂泥岩9.200.301.80037.1820.02 450煤层0.500.340.30030.5510.01 378C3t砂泥岩4.260.311.40032.0780.02 541
5 稳定性数值模拟结果分析
在自然采空情况下主要模拟3煤、8煤的开采过程。模拟前首先将系统初始应力状态下的各个方向的位移和速度清零，确保之后的位移、速度完全是由煤层的开采引起。从开始开采3煤到8煤采完，系统都能调整应力达到稳定状态。
5.1 水平位移分析
3煤、8煤采空后的水平位移云图如图2所示。地表X方向上，在山坡的顶部产生了负向的水平位移，位移约为-0.022 6 m；在坡腰和坡脚地带，同样产生了负向的水平位移，其水平位移由坡顶向四周逐渐减小。表明地表坡体各部分的移动距离是受采空区的几何中心所控制，当采空区形成后，其几何中心处的顶板岩块首先开始下沉(将该岩块称之为关键块)，随后关键块的移动带动周边岩块与之下沉，在地表则呈现近似条带状的水平移动。由于采空区上方岩土体的重度分布不均匀，使得由采空而形成的沉陷，在水平方向上的位移量也极不规则，加上3煤工作面的推进距离明显大于其开采深度，所以造成移动盆地内的多个点均达到或者超过临界开采点，造成地表移动盆地的面积远大于采空区的面积。
图2 3煤、8煤采空后地表水平位移
Fig.2 Horizontal displacement cloud map of surface
after mining of No.3 and No.8 coal seam
在采空区两侧边缘附近，岩土体受采动影响，往采空区方向位移，呈受拉状态，易出现拉裂缝，与实际观测发现裂缝所在位置较吻合，且左侧为凸形坡顶，右侧为凹形坡脚，左侧拉裂程度较右侧强烈；而采空区中间上方地表，受到两侧岩土体的挤压，处于受压状态，X方向位移很小。8煤采后，负向最大位移增加3 mm，正向最大位移增加4 mm。因此水平位移主要由3煤的开采引起，8煤的开采对水平位移的影响很小。
沿着Y方向不同部位(50、100、150、200、300、400 m)切剖面如图3所示。不同剖面生成的X方向位移如图4所示。
图3 沿Y方向不同剖面位置
Fig.3 Location of different sections along the Y direction
图4 各部面水平位移分布
Fig.4 Horizontal displacement distribution cloud
由于煤层开挖形成采空区，在不同部位，产生的位移不同。在采空区南侧，无煤层开采部位，如Y=50 m剖面，X方向位移云图变化较为均匀，表明其位移变化较小。接近采空区的部位，如Y=100 m，Y=150 m剖面，由于采空区的存在产生了负向位移。在采空区中心部位，即煤层最大开采部位，如Y=200 m剖面，X方向正向位移最大。在采空区北侧，距离采空区中心位置越远，X方向位移越小。
5.2 垂直位移分析
煤层的采空必将导致围岩的变形与位移。采空区变形随其面积与体积的增大而不断加剧，甚至顶板冒落，还可能波及到周边采场、巷道的稳定。从开采的3煤、8煤采空后，垂直位移如图5所示，仅选取代表性剖面X=200 m进行叙述。
图5 3、8煤采空后X=200 m剖面垂直位移
Fig.5 Z-disp cloud diagram of X=200 m section after No.3 and No.8 coal mining seam
3煤开采后，地表最大沉陷位移为0.046 m，X=200 m剖面沉陷云图如图5a所示，此时顶板已沉陷到底板。开采8煤后，地表最大垂直位移0.067 1 m，X=200 m剖面垂直位移如图5b所示。
沿Y方向不同位置切剖面(图3)生成的不同剖面Y方向垂直位移如图6所示，顶板中间位移最大，顶板往上位移等值线呈弧形扩散，并越来越小。底板上有正向位移，说明底板有少量鼓起，是采煤后释放应力引起的底鼓现象。
图6 垂直位移
Fig.6 Vertical displacement cloud each prfile
由于3煤、8煤为缓倾斜赋存煤层，倾角约为4°。煤矿采出后，整体位移分布基本服从近对称特性规律，并略向煤层倾向方向倾斜。采场周边围岩位移最大，往外距离开挖边界越远，围岩位移就越小，围岩移动方向均指向采空区。在采空区顶板以上形成围岩移动位移等值拱，越往上发展拱径越大，近地表附近等值线则发散开来。从顶板往上，位移越来越小。采空区顶板中央沉陷量最大，向采空区两侧则急剧降低，采空区以外则沉陷量很小。在地表，沉陷位移等值线以采空区为中心呈弧状分布，同时沉陷量往外逐渐减小。
5.3 采动坡体稳定性定性分析
判断滑坡体的(潜在)滑动面(带)，可根据其剪应变增量的大小来判断：剪应变增量较大(绝对值)的部位，则为其(潜在) 滑动面(带)，变形破坏也多沿此处发生；剪应变增量较小或基本上没有发生变化的部位，一般不会有潜在滑动面的产生，因此，这些部位也不会发生较大的变形或破坏[10]。
3煤、8煤采空后不同剖面剪应变如图7所示。从图中可以看出，剪切带主要分布在采空区端部，表现为上覆岩体沿此剪切带的下切错落。近地表附近并没有有利于坡体滑动的剪切带出现。
图7 天然采空工况下剪应变增量云图
Fig.7 Shear strain incremental cloud map under natural mining conditions
根据3煤、8煤煤采空后水平位移云图可知，以采空区中部为界，上坡方向岩土体有向下坡方向的位移；下坡方向岩土体有向上坡方向的位移。这将导致采空区两边缘上覆岩土体出现拉裂，即在坡顶坡脚出现裂缝。但在采空区中部上方，岩土体受到两侧岩土的挤压，位移很小。总体上看，采空区两侧岩土体向采空区方向滑移。
根据变形体的基本特征和形成的主导因素，由于整个坡体坡度较大，坡面呈现出不平整的阶梯状，重心较高，且变形体上已经形成了众多的拉长裂缝，会更有利于降水的入渗，极易在其表层形成变形。
综上所述，采空区上坡方向岩土体的移动为开采引起的采动滑移。虽然在坡顶坡脚均有裂缝产生，但大多为拉张性裂缝，或由不均衡拉张引起的局部剪切裂缝。坡体也并非整体滑移，而是都向采空区中部移动，并且移动范围为以采空区为中心的圆弧状。在剖面上，并没有出现有利于坡体向下坡方向滑动的剪切带。因此，坡体在天然采空工况下会发生开采沉陷，但不会发生采动滑坡。
5.4 采动坡体稳定性定量分析
使用强度折减法对斜坡进行稳定性计算，为了定性评价变形体的稳定性，本次对变形体现场实测了3条剖面，选取其中1条典型剖面作为计算分析的剖面，使用之前所给的岩土力学参数，由于研究区抗震设防烈度为Ⅶ度，设计基本地震动峰值加速度为0.10 g(g为重力加速度)，根据相关规范计算时可不考虑地震对变形体的影响。由于变形体上方的斜坡下有采空区，根据前面的定性分析和数值模拟结果，变形体稳定性计算考虑以下2种工况：①工况一：自重；②工况二：自重+采空沉陷形成的侧向推力(斜坡岩体侧向变形产生的剪应力)。对所给的剖面进行多次计算，计算结果见表2。
表 2 变形体稳定性计算结果
Table 2 Calculation results of deformation stability
计算工况剖面号稳定系数稳定状态规范规定安全系数说明工况一1-11.738稳定1.300满足工况二1-11.154基本稳定1.250不满足
根据《建筑边坡工程技术规范》(GB 50330)、《工程地质勘察规范》(DBJ 50 T-043-2016)，其稳定性评价标准为：①当K<1时，不稳定；②1.00<K<1.05，欠稳定;③1.05<K<Kt (规定安全系数)，基本稳定;④K>Kt(规定安全系数)，稳定。
6 结 论
通过对山西某煤矿采空区形成过程以及采动坡体在新老采空区影响下的变形破坏的模拟，经过一系列分析研究，得出3条结论。
1)采空区及两侧边缘地区的水平位移具有一定的规律：采空区南侧部位水平位移变化较小；接近采空区的部位产生负向位移；在采空区中心部位，正向水平位移达到最大；在采空区北侧，距离采空区位置越远，水平位移越小。并且，水平位移主要由3煤的开采引起，8煤的开采对水平位移的影响很小。
2)在采空区顶板以上形成围岩移动位移等值拱，越往上发展拱径越大，近地表附近等值线则发散开来。在地表，沉陷位移等值线以采空区为中心呈弧状分布，同时沉陷量往外逐渐减小。
3)采动坡体稳定性定性分析的结果与定量分析的结果相吻合。在天然采空工况下，剪切带分布于采空区上方及两侧，完全采空后延伸到坡顶，没有出现有利于坡体滑动的剪切带；且通过定量分析，坡体的稳定系数为1.738，大于规范规定的安全系数1.300。因此，坡体是稳定的，不会出现滑坡。但在采空区岩体因沉陷变形产生的侧向推力为主的作用下，坡体的稳定系数为1.154，根据相关规范边坡稳定性评价标准，稳定系数大于1.050且小于规范规定的安全系数1.250，坡体处于基本稳定状态。
[1] 汤伏全,田家琦,于 根.地下开采诱发滑坡的机制分析[J].辽宁工程技术大学学报：自然科学版，1991,13(S1): 16-17.
TANG Fuquan, TIANJiaqi, YU Gen. Mechanism analysis of landslide induced by underground mining [J].Journal of Liaoning Technical University:Natural Science,1991,13(S1):16-17.
[2] 梁 明,汤伏全.山体下采煤的地面监测与数据分析[J].西安科技大学学报,1994,14(2):145-150.
LIANG Ming,TANG Fuquan.Surface monitoring and data analysis of coal mining under mountain [J].Journal of Xi′an University of Science and Technology,1994,14(2):145-150.
[3] 汤伏全,梁 明.地下开采影响下山体的稳定性分析与评价[J].矿山测量,1995(1):23-26.
TANG Fuquan,LIANG Ming.Stability analysis and evaluation of the downhill influenced by underground mining[J].Mine Surveying,1995(1):23-26.
[4] 汤伏全,梁 明.地下采矿诱发山体滑坡机制的研究[J].煤矿开采,1995,1(3):34-37.
TANG Fuquan,LIANG Ming.Research on the mechanism of landslide induced by underground mining[J].Coal Mining Technology,1995,1(3):34-37.
[5] 陆 力,魏衍栋,常夫伟,等.基于RS和GIS的矿山泥石流形成条件信息提取[J].山西煤炭,2014,34(11):1-3.
LU Li,WEI Yandong,CHANG Fuwei,et al. Extraction of conditions for formation of debris flow in mines based on RS and GIS[J].Shanxi Coal,2014,34(11):1-3.
[6] 何姣云.矿山采动灾害监测及控制技术研究[D]. 武汉：武汉理工大学,2007.
[7] 胡晋山,何宗宜,康建荣,等.山区地下开采坡体稳定性分析及滑坡预报系统的设计[J].中国矿业, 2009,18(11):112-115.
HU Jinshan,HE Zongyi,KANG Jianrong,et al. Stability analysis of landslides in mountainous areas and design of landslide prediction system[J].China Mining Magazine,2009,18(11):112-115.
[8] TANGSB,HUANG R,TANGCA,et al. The failure processes analysis of rock slope using numerical modelling techniques[J]. Engineering Failure Analysis,2017,79:999-1016.
[9] 毛新虎,周永昌,张国文.矿山泥石流地质环境的模糊综合评价[J].西南民族大学学报：自然科学版,2008,34(5):929-934.
MAO Xinhu,ZHOU Yongchang,ZHANG Guowen. Fuzzy comprehensive evaluation of geological environment of mine debris flow[J].Journal of Southwest University for Nationalities:Natural Science Edition,2008,34(5):929-934.
[10] 孟 伟,张立伟,曾德礼.不同因素作用下采空区对地表稳定性影响分析[J].路基工程,2009，37(1):84-86.
MENG Wei,ZHANG Liwei,ZENG Deli. Analysis of influence of goaf on surface stability under different factors[J].Subgrade Engineering，2009，37(1):84-86.
[11] 胡海峰,康建荣.基于有限元法的采动坡体稳定性计算[J].太原理工大学学报,2000，31(4):343-345,353.
HU Haifeng,KANG Jianrong. Calculation of stability of mining slope based on finite element method[J]. Journal of Taiyuan University of Technology,2000，31(4):343-345,353.
[12] 牛 威.煤矿采空塌陷导致土地破坏研究：以山西西山矿区为例[J].中国地质灾害与防治学报, 2006,17(4):163-164.
NIU Wei.Study on land destruction caused by coal mine caving collapse：taking Xishan Mining Area in Shanxi as an example[J]. The Chinese Journal of Geological Hazard and Control,2006,17(4):163-164.
[13] 丁勇春,王建华,徐 斌.基于FLAC3D的基坑开挖与支护三维数值分析[J].上海交通大学学报,2009，43(6):976-980.
DING Yongchun,WANG Jianhua,XU Bin. Three-dimensional numerical analysis of foundation pit excavation and support based on FLAC3D[J]. Journal of Shanghai Jiaotong University,2009，43(6): 976-980.
[14] 贾新果.采煤沉陷区地表残余沉陷时间函数模型研究[J].煤炭科学技术,2018,46(11)：157-162.
JIA Xinguo. Study on time function model of surface residual subsidence in coal mining subsidencearea[J].Coal Science and Technology,2018,46(11):157-162.
[15] 高富强,高新峰,康红普.动力扰动下深部巷道围岩力学响应FLAC分析[J].地下空间与工程学报, 2009，29(4):680-685.
GAO Fuqiang,GAO Xinfeng,KANG Hongpu.FLAC analysis of mechanical response of surrounding rock in deep roadway under dynamic disturbance[J]. Chinese Journal of Underground Space and Engineering,2009，29(4):680-685.
[16] YU Yang,WANG Enzhi,ZHONG Jianwen,et al. Stability analysis of abutment slopes based on long-term monitoring and numerical simulation[J].Engineering Geology,2014,183:159-169.
[17] 宋卫东,付建新,杜建华,等.基于精密探测的金属矿山采空区群稳定性分析[J].岩土力学,2012,33(12):3781-3787.
SONG Weidong,FU Jianxin,DU Jianhua,et al. Stability analysis of goafs in metal mines based on precision detection[J].Rock and Soil Mechanics,2012,33(12):3781-3787.
[18] 王海军.煤层顶板沉积环境对其稳定性影响研究[J].煤炭科学技术,2017,45(2):178-184.
WANG Haijun. Study on the influence of sedimentary environment of coal seam roof on its stability[J].Coal Science and Technology,2017,45(2):178-184.
[19] 郭 曼,钱素娟.基于FLAC3D的采空区群稳定性及处理方案模拟[J].矿业研究与开发,2016,36(7):93-96.
GUO Man,QIAN Sujuan. Simulation of goaf group stability and treatment scheme based on FLAC3D[J].Mining Research and Development,2016,36(7):93-96.
[20] 杨 鹏.采场上覆岩层采动裂隙演化规律相似模拟研究[J].煤炭科学技术,2014,42(8):121-124.
YANG Peng.Simulation of mining fracture evolution law of overlying strata in mining face[J].Coal Science and Technology,2014,42(8):121-124.
