甘智慧，尚 慧，杜荣军，占惠珠
(西安科技大学 地质与环境学院，陕西 西安 710054)
近年来，随着地下煤炭资源的采掘活动日益剧增，从而引发了一系列采空区地面沉陷问题，严重影响了社会经济发展和人类正常的生产生活。
众多学者借助现场实测[1]、理论分析[2-3]、数值模拟[4-5]和物理相似模拟试验[6-7]等方法对采煤引起的地面沉陷问题进行了大量研究。许家林等[8]采用物理和数值模拟方法，就覆岩主关键层对地表下沉动态过程的影响进行了研究；尹光志等[9]通过相似模型试验与数值分析相结合的方法，对大倾角煤层深部开采的采岩移动基本规律、矿山压力分布规律和地表移动变形规律进行了研究；孙学阳等[10-11]采用数值模拟与相似模拟实验相结合的方法研究双煤层开采对覆岩的破坏影响，结果显示，工作面留设煤柱宽度越大，煤层开采对其覆岩的影响越小；王永国等[12-13]采用理论分析、现场实测、数值模拟、水文监测等相结合的研究方法，深入研究覆岩运移特征及其主控影响因素；黄庆享等[14-15]采用物理模拟、数值计算和现场实测相结合的方法，研究浅埋煤层群开采覆岩及地表裂隙演化规律及机理。
综上可知，煤层开采引起的地面沉陷问题多采用2 种或2 种以上研究方法，但利用多时相DEM数据进行采煤沉陷的时空演化规律分析鲜见报道。笔者以宁夏石嘴山矿区为对象，建立研究区三维地质模型，模拟实际开采情况，分析不同开采期覆岩移动变形、应力分布及塑性区变化规律，并将地表位移模拟结果与多时相DEM 数据叠加分析结果进行相互验证，以期为缓倾斜煤层开采导致地面沉陷问题的研究提供新思路，为提高矿区开采沉陷监测效率提供新方法。
宁夏石嘴山矿区位于黄河中游上段，宁夏回族自治区石嘴山市惠农区境内，东临黄河，西靠贺兰山，北与内蒙古自治区阿拉善盟、乌海市伊克昭盟接壤，属典型的大陆半干旱气候，常年干燥多风，日照充足、雨量稀少，蒸发强烈。
石嘴山矿区地貌属于石嘴山侵蚀台地，地势西高东低；构造位置位于东鄂博梁背斜以东，总体呈北东向较宽缓的向斜构造，向斜轴自北东向南西蜿蜒伸展，略呈“S”型，两翼地层不对称，东南翼稍缓，地层倾角18°～32°，西北翼略陡，地层倾角13°～45°；出露地层自上而下依次为：第四系、新近系、三叠系、二叠系、石炭系、寒武系、蓟县系及长城系，主要含煤地层为石炭系上统太原组和二叠系下统山西组[16]，含煤地层总厚为185.36 m；水文地质条件中等，地下水类型属坚硬基岩裂隙水。
宁夏石嘴山矿区主要包括石嘴山一矿和石嘴山二矿(图1)，为地下开采[17]，发育煤层9 层，煤层总厚为32.09 m，平均采深285 m。石嘴山矿区已有50余年开采历史，地表变形严重，总体上形成7 个较大的塌陷坑，塌陷总面积为9.1 km2，最大塌陷深度为24.39 m[16]。采空区塌陷造成矿区附近房屋产生裂缝、有效耕地面积减少、井水河流污染等一系列地质环境问题[3]，因此，对该矿区采煤引发地面沉陷问题进行研究已迫在眉睫。
本文选取石嘴山一矿Ⅱ—Ⅱ′剖面作为研究对象，该剖面沿煤层倾向横跨3 号、4 号塌陷坑，具有一定代表性(图1)。石嘴山一矿走向长约4.35 km，倾向宽约 1.19 km，面积 5.18 km2。开采高程为+1 055～+600 m，煤层平均倾角20°，为缓倾斜煤层，其中二、三、五、六、七、九号煤属于稳定可开采煤层，煤层覆岩主要为砂岩和砂质页岩(图2)。根据井田内煤层赋存条件、多煤层可分组特点，将可采煤层分为上下两组，其中二、三煤层为上组煤，五、六、七、九煤层为下组煤。煤层开采采用斜井、分别集中运输的联合开拓方式，分为3 个水平。一水平+1 050～+974 m；二水平+974～+830 m；三水平+830～+600 m。一水平已于1968 年回采完毕，二水平已于1989 年2 月回采完毕[16]，目前正在开采三水平，亚阶段+680～+600 m。
图1 石嘴 山矿区位置Fig.1 Location of Shizuishan Mining Area
图2 石嘴山矿区Ⅱ—Ⅱ′工程地质剖面Fig.2 Engineering geological section at row of Ⅱ-Ⅱ′ in Shizuishan Mining Area
FLAC(Fast Lagrangian Analysis of Continua)是由美国Itasca 公司开发的(连续介质力学分析软件)显式有限差分程序[18]，采用混合离散法和动态松弛法，能较好地模拟地质材料渐近破坏和失稳，适用于模拟大变形、非线性问题；在模拟材料的屈服过程中，采用混合离散化方法模拟塑性破坏与塑性流动，比有限元有效[19]。
为研究石嘴山矿区地面沉陷及覆岩移动变形规律，以Ⅱ—Ⅱ′剖面为对象，运用FLAC3D软件进行数值模拟分析。依据Ⅱ-Ⅱ′剖面上钻孔所揭露的地层情况(图2)、工作面范围及采空区特点建立三维数值模型(图3)。该区域地层结构复杂，在数值模型建立过程中，对物理力学性质相近的岩土体进行适当简化，将煤层按上下两组简化为两层煤，平均厚度均为10 m，两组煤同时开采，煤层开采遵循由上向下分3 水平开采的顺序，最终确定模型尺寸长(X)×宽(Y)×高(Z)设置为2 200 m×600 m×850 m，包括79 560 个节点和86 121 个单元格。结合实际开采条件，石嘴山矿区煤层开采模拟时计算步骤如下：①单元网格化分利用FLAC3D5.0 Extrusion 建模；② 初始应力场求解；③初始应力位移清零；④ 根据矿体实际开挖顺序分3 水平进行开采，并在各开采水平间预留煤柱15 m。
图3 三维数值模型Fig.3 Three dimensional numerical model
模型基本假定为各向同性连续均匀介质，力学模型采用莫尔–库伦弹塑性模型。模型顶部地表为自由面，底部采用固定约束，四周采用水平位移约束。矿区初始应力以自重应力为主，设置重力加速度为10 m/s2，方向垂直向下，岩层物理力学参数见表1。
表1 岩层物理力学参数[16]Table 1 Physical and mechanical parameters of rock strata[16]
2.3.1 应力场分析
采矿工程破坏原岩应力，周围岩体应力重新平衡，维持采空区稳定的部分围岩会产生局部应力集中现象。随着开采深度的增加，应力逐渐增加，当超过其强度极限时就会发生损伤、破坏甚至采空区整体失稳。本模型分3 水平开挖，各水平开挖后采场周边岩层的煤层走向方向和倾向方向的应力分布情况如图4 所示。在沿矿体走向方向中心处(Y=300 m)设置平行矿体倾向方向剖面，该剖面各水平开采后应力场变化如下：
图4 各开采水平倾向及走向方向应力云图Fig.4 Stress cloud maps of inclination and strike direction of each mining level
①第一水平开采应力场分布呈现从上往下迭代递增的变化趋势，最大主应力(压应力)均小于22.1 MPa。煤柱出现明显的应力集中，竖向压应力约为4 MPa。煤层顶板出现拉应力，最大拉应力出现在顶板中间偏上的位置(图4a)；
② 随着工作面推进，第二水平开采的采空区围岩压应力扰动效应较为明显，采空区两侧煤柱出现高压应力集中，随着与煤柱距离的增大，压应力逐渐变小，但变化梯度远小于围岩应力变化梯度。采空区顶板出现拉应力集中，受拉破坏并产生断裂破碎、垮落(图4b)；
③第三水平开采后，采空区围岩最大拉应力达到3.2 MPa。顶板拉应力也不断增大，顶板角隅处出现拉应力集中，是影响采空区稳定性的主要因素。煤柱应力集中较为明显，随着开采深度增加，压应力增加明显，最大压应力约5 MPa。拉应力和压应力共同作用形成明显的剪切破坏，以压应力为主，其最大值达到19.28 MPa。上下两层煤层距离较近，采空区围岩出现应力叠加，表现出“群效应”，应力叠加区域的岩体更容易发生破坏(图4c)。
针对沿矿体走向方向剖面(X=600 m)在Y轴方向上的应力场进行分析。随着开挖进程不断推进，采空区体积不断增大，围岩应力释放明显。走向应力形态基本呈现对称分布，且由上到下逐渐增大。
开采第一水平，右围岩主要表现为压应力，最大值约为17.9 MPa，在采空区中央形成应力等值线轨迹拱，角隅处出现应力集中，采空区直接顶板也表现为应力，其值为0.065 MPa(图4d)；开采第二水平后上组煤顶板压应力由 0.065 MPa 减小为0.042 MPa，采空区顶板处于卸压状态(图4e)；开采第三水平主要表现为压应力，最大值为17.6 MPa。采空区顶板应力最小值都超过了顶板岩体的极限抗拉强度，顶板岩层出现拉张破坏，引起采空区顶板垮落(图4f)。
2.3.2 位移场分析
矿体开采必然引起围岩的变形和位移，随着采空区面积和体积的增大，围岩变形也持续增大，最终出现顶板垮落。各开采水平倾向及走向垂直位移数值模拟结果如图5 所示。
从倾向方向的位移云图(图 5a—图 5c)可以得出：煤层开挖后，煤层顶板岩层出现悬空，顶板下方没有支撑体，其顶板岩层出现垮落、坍塌，远离采空区的覆岩受矿层开采影响较小。随着采空区推进，移动现象延伸至地表并形成沉陷盆地，其在地表的影响范围远大于采空区实际尺寸。沉陷盆地不对称，上部岩层不仅向采空区移动，而且因为受重力作用由上山方向向下山方向位移。
图5 各开采水平倾向及走向方向位移云图Fig.5 Displacement cloud maps of inclination and strike direction of each mining level
一水平煤层开采时，采空区顶板垂直位移最大值–1.3 m，影响范围在开采部位正上方；地表沉降最大值约0.6 m。上组煤底板出现隆起现象，其最大位移量为0.17 m，这是由于煤层回采后导致采空区应力释放，底板岩层受到向上的力，导致出现上拱现象(图5a)。
二水平开采完成后，一水平最大变形位置的垂直位移由–1.3 m 增加到–5 m；随着开采深度加大，最大位移出现在上组煤第二水平顶板中心偏上山方向，约–7.4 m。采空区底板隆起达到最大值，位移量为0.73 m，采空区中间煤柱受到的上、下围岩挤压作用最为严重。地表沉降明显增大，地表沉降最大值从第一水平的 0.6 m 增加到第二水平的3.5 m(图5b)。
开采第三水平后最大垂直位移位置与二水平一致，从二水平–7.4 m 增加到–19.1 m，其正上方地表沉降值为12 m，对地表沉降影响显著；而三水平采空区顶板垂直位移–2～–8 m，对地表沉降影响较小，这是因为煤层三水平采深为 320～550 m，采高10～15 m，采高比均大于30(图5c)。
对平行煤层走向剖面(X＝600 m)在Y轴方向上的位移场进行分析。第一水平开采导致顶板沉降，位移最大值约1.0 m，地表沉降约为0.8 m，煤层底板稍有隆起(图5d)；随着工作面继续推进，开采第二水平时，底板出现隆起，最大底鼓量为0.09 m，顶板位移量比底板大，最大位移出现在上组煤采空区顶板中间偏上山方向，约为5.18 m。采空区周围覆岩位移最大，往外距离采空区越远位移越小，且围岩位移方向都指向采空区(图5e)；开采第三水平后最大沉降量为10 m，由于煤层开采引起采空区底板由三向受力状态变为水平二向受力状态，垂直方向应力释放引起底板回弹隆起，最大达到0.19 m。整体位移分布具有左右对称性规律，采空区围岩位移变化峰值均处在顶板及底板中间位置(图5f)。
通过 FLAC3D数值模拟得到地表沉降等值线(图6)。煤层回采后在地表形成2 个巨大的沉陷盆地，其中较小的沉陷区面积约为36 520 m2，最大沉降为7 m；另一个较大的地面沉陷区域呈“类W 型”，最大沉降量为12 m，影响范围约259 380 m2，地表垂直位移略小于采空区最大沉降量，但地面沉陷范围远大于采空区范围。地表垂直位移等值线图近似呈现椭圆形，采空区地表影响区分别由2 个沉降中心向四周扩散，等值线逐渐稀疏。由于研究区为倾斜煤层，因此，地面沉陷盆地不对称，2 个沉降中心均发生在沉陷盆地中部且偏下山方向，下山方向比上山方向影响范围更大。
图6 最终地表位移下沉等值线Fig.6 The final surface displacement subsidence contour
2.3.3 塑性区
塑性区形态可作为开挖扰动区的参考依据，因此，可基于数值模拟结果对围岩开挖塑性区扩展进行分析[20]。基于FLAC3D模拟得到塑性区分布状况(图7)，模型单元在塑性状态下将直接呈现出拉张或剪切破坏，可直观反映围岩稳定性。通过塑性区分析可知，塑性区集中在煤层顶板、底板及煤柱。由于煤层倾斜，煤柱及邻近采空区围岩基本呈剪切破坏，地表及底板局部受到拉张破坏，表明剪切破坏对采空区的稳定性影响较大。开采第一水平顶板出现贯通剪切破坏，顶板垂直位移较大，应力集中明显(图 7a)；随着煤层开采深度增加，塑性区向深度扩展，开采第二水平围岩剪切破坏扰动范围扩大，而拉张破坏影响范围较小，主要发生在地表及煤层底板角隅处(图 7b)。巷道底部煤层受到剪切破坏，导致开采层下部卸压煤体瓦斯流入开采层采空区，应对巷道下部的卸压煤体区域重点进行瓦斯抽采工作；随着掘进的深入塑性区不断向围岩深部扩展，开采第三水平塑性区影响范围扩大到采空区地表，地面沉陷量达到最大值(图7c)。
图7 各开采水平塑性区分布Fig.7 The horizontal plastic zone distribution of each mining level
数字高程模型(Digital Elevation Model，DEM)作为地形表面形态的一种离散数字表达，能够直观、全面地反映地表位移变化情况。为实现对采空区塌陷坑深度及范围监测，分别选取矿区20 世纪70 年代(1︰50 000)和2003 年(1︰10 000)等高线生成DEM(图8)，且将两时相DEM 进行叠加分析，得到20 世纪70 年代至2003 年期间地表垂直位移变化图(图9)。小部分高程增加(即图9 中红色区域)，但多数是由于城镇化建设造成。其他区域高程均有所降低，沉降幅度较大处均为塌陷坑，其中Ⅱ—Ⅱ′剖面穿过的3 号、4 号塌陷坑塌陷形态与数值模拟地表位移下沉范围和形状(图6)基本一致。
图8 石嘴山矿区两时相DEMFig.8 Bi-temporal images DEM in Shizuishan Mining Area
图9 石嘴山矿区20 世纪70 年代—2003 年沉降变化幅度Fig.9 The subsidence variation range chart of Shizuishan Mining Area from 1970s to 2003
为更进一步了解采空区地面塌陷垂直位移的变化情况，基于DEM 自动提取两时相地形进行叠加对比，得到Ⅱ—Ⅱ′剖面地形线变化对比图(图10)。靠近惠农区城区的0～1 000 m 处地形明显上升，大多是由于在城镇化建设中新增建(构)筑物造成；1 000～3 000 m 处地形出现明显不均匀沉降，分别形成3 号和4 号塌陷坑。其中3 号塌陷坑范围为167 m，竖向位移最大值为–4.2 m；4 号塌陷坑范围为933 m，竖向位移最大值为–12.2 m。
将Ⅱ—Ⅱ′剖面基于两时相DEM 差值运算结果与FLAC3D数模模拟得到的地表沉降曲线进行对比分析。由于两沉降曲线不完全重合，选取FLAC3D模型原点为基准点，即Ⅱ—Ⅱ′剖面地形线变化对比图1 100 m 位置处(图10)。数值模拟与DEM 分析结果对比如图11 所示，其中模拟结果的3 号塌陷中心水平位置在240 m 处且垂直位移约为–7 m，而在30～70 m 位置引起地表隆起现象，其最大值为1.24 m，这是由于模型边界效应引起。两时相DEM叠加统计得到3 号塌陷坑塌陷中心在60 m 处，垂直位移为–4.2 m。对比发现2 种方法所得3 号塌陷坑位置及沉降量相差较大，其主要原因是一水平煤层开采深度较浅，切眼接近地表，因此，存在一定的误差。4 号塌陷坑处2 种方法确定的塌陷中心水平位置基本相同，沉降趋势基本吻合，且两者变形量也基本一致，相差均不超过3.08 m，大部分点位相差在1 m 以下。数值模拟得到的地面沉陷值整体略大于两时相 DEM 数据垂直变化量，主要原因是DEM 数据来源分别是20 世纪70 年代和2003 年，而FLAC3D数值模拟的时间范围是从1956 年规划开采到矿山开采初期结束，FLAC3D的时间数据范围远大于DEM，其开采深度及开采范围也更大。综合以上分析可知，数值模拟计算的沉降量与两时相DEM叠加统计分析的变化量结果及趋势基本一致，计算结果合理。
图10 Ⅱ—Ⅱ′剖面地形线变化对比图(20 世纪70 年代—2003 年)Fig.10 Contrast of section Ⅱ-Ⅱ′ topographic line change contrast diagram(1970s-2003)
图11 两时相DEM 叠加分析与FLAC3D 数值模拟地表位移对比Fig.11 Comparison of surface displacement by bi-temporal images DEM superposition analysis and FLAC3D numerical Simulation
a.通过FLAC3D模拟采煤过程，煤层开采后岩层主应力呈现自上而下逐渐增大的变化趋势，采空区上山方向出现拉应力，最大值约3.2 MPa，煤柱出现明显的应力集中，最大应力值达到5 MPa。
b.地面沉陷盆地不对称，两个塌陷坑沉降中心均发生在沉陷盆地中部且偏下山方向，下山方向比上山方向影响范围大，2 个沉陷盆地最大沉降量分别为7 m 和12 m。
c.采空区围岩以剪切破坏为主，塑性区范围随煤层开采深度增加逐渐扩展，主要发生在地表及采空区围岩角隅处。
d.通过矿区20 世纪70 年代与2003 年DEM数据进行叠加分析得到：3 号塌陷坑最大垂直位移为–4.2 m，4 号塌陷坑最大垂直位移为–12.2 m。而数值模拟计算得到的两塌陷坑最大沉降量分别为–7 m 和–12 m，两者计算结果对比发现，3 号塌陷坑塌陷中心水平位置差异较大，4 号塌陷坑塌陷中心水平位置与沉降趋势基本吻合。DEM 数据可与数值模拟相互验证，对提高矿区沉陷监测效率有一定参考价值。