张海君，吴奕枢，王 飞，谢昆仑

(1.鄂尔多斯市国源矿业开发有限责任公司，内蒙古 鄂尔多斯 017000；2.郑州大学 力学与安全工程学院，河南 郑州 450001)

0 引 言

井下煤炭资源采出后，煤层顶板悬露，顶板原岩应力遭到扰动，顶板岩石在重力作用下将出现移动变形。随着工作面的继续推进，顶板出现破坏、冒落，导致地表沉陷现象的发生，将会对矿区及矿区周边的建筑物、生产设施和人员的安全造成威胁。随着中国煤炭资源的持续勘探和开发，煤炭开采所面对的矿区环境更加复杂，如何有效预测和防范开采造成的地表沉降危害，已成为矿区整体规划设计的一个重点。

井下开采所致的工作面顶板岩移及地表沉陷现象是一个涉及复杂岩体力学的变化过程。以往国内外学者提出了许多地表沉陷模型和理论[1-4]，典型的如基于观测数据的正态理论、二等分线理论、自然斜面理论和开采沉陷“拱理论”[5-6]、工作面上覆岩层“三带”分布理论[7]等。随着科学技术水平的不断提升，开采沉陷理论得到进一步发展。Salamon将影响函数预计法与连续介质理论相结合，开创了用于探究开采沉陷规律的边界元法(BEM)[8]；DONNELLY等基于SWIFT技术对矿山地表沉陷问题进行预测与分析[9]，BARYAKH等基于线性粘弹性理论，利用蠕变函数开发了一种动态地表沉陷预计程序[10]；崔希民等应用流变模型对地表沉降规律进行分析[11]；朱广轶等基于概率积分法提出了一种动态地表沉陷时间函数，并通过实例分析证明其具有比Knote函数更优的精准度[12]；余学义等根据观测数据研究单一工作面开采和多工作面开采后地表的沉陷规律，并分别评估了2种开采场景下地表的损害情况[13]；刘玉成对Knote时间函数加以改良，建立了一种更加精准的地表沉陷预计模型[14]；石晓宇等探讨了2种灰色模型在地表沉陷监测与预计方面的应用，并指出二者适用的开采情况[15]；黄明江等指出地表下沉曲线符合分形增长规律，并借助分形插值法实现了有限观测数据下对地表沉陷连续变化情况的拟合[16]。此外，一些学者还采用GIS，InSAR，ArcObjects，RTK，无人机激光雷达等新型观测技术和地理信息系统对矿山地表沉陷情况展开监测与分析[17-20]。

中国许多大型矿区属于厚煤层，煤炭储量相当丰富。近年来，综放开采技术被运用于厚煤层的开采过程中，因其具有产量大、效率高、成本低等特点，而备受业内认可。关于厚煤层开采引起的覆岩和地表移动规律等方面目前已有不少研究成果，如郭文兵等对“三软”特厚煤层地表沉陷情况进行观测，并基于概率积分法运用Matlab拟合了地表下沉曲线[21]；胡青峰等建立了塔山煤矿特厚煤层倾向主断面相似模拟实验，并研究了重复采动下覆岩变形和地表沉陷规律[22]；高超等针对东坡煤矿浅埋深、特厚煤层综采等特点，结合地表实测数据，对地表动态移动变形特征进行分析[23]；赵兵朝等通过对郭家河煤矿特厚煤层综放开采实地观测和分析求取了地表移动变形有关参数，并综合运用反演模拟和边坡稳定性分析对所得参数进行修正[24]，但所研究的对象较少涉及20 m以上的大采放比特厚煤层综放开采，故仍需对此类作业引起的覆岩与地表移动规律进行系统地研究。

文中将根据准格尔煤田龙王沟煤矿特厚煤层一次采全厚这一开采特征，建立基于概率积分法的地表沉陷预计模型，运用该模型分析煤层开采引起的地表移动规律与变形特征，同时讨论该模型存在的局限性，并采用FLAC3D数值模拟方法加以补充，从而获得较为准确、全面的特厚煤层综放开采的地表沉陷规律及其特征，为合理规划煤矿开采作业方案，保护矿区人员生产生活的安全提供指导方向和理论支持。

1 龙王沟煤矿61601工作面概述

龙王沟井田位于内蒙古自治区鄂尔多斯市准格尔煤田中北部，面积41.35 km2，属典型黄土高原地貌。依据地表出露及钻孔揭露情况，该井田地层自下而上分别为：奥陶系中统马家沟组(O2m)，石炭系上统太原组(C2t)，二叠系下统山西组(P1s)、下石盒子组(P1x)，二叠系上统上石盒子组(P2s)、石千峰组(P2sh)，三叠系下统刘家沟组(T1l)，第三系上新统(N2)和第四系(Q)等。煤田地层结构形态为一近似于南北走向、呈波状起伏、向西侧倾斜的单斜构造，地层倾角较小，低于10°。61601工作面布置于龙王沟煤田太原组上部6煤处，为首采面，平均采深399.9 m，工作面倾向长度254.6 m，走向可采长度624 m，煤均厚23.15 m，煤层稳定。工作面设计采高5.1 m，顶煤厚18.05 m，采放比1∶3.54。煤层倾角0～2°，接近水平煤层。煤层顶底板各岩层如图1所示。61601工作面开采选用后退式走向长壁采煤方法，综采放顶煤一次采全高开采工艺，完全垮落法管理顶板。

为获得特厚煤层开采条件下地表沉陷的实测数据，分析地表沉陷规律，根据《煤矿测量规程》在61601工作面外设置地表沉陷观测站，根据《煤矿测量规程》规定，当工作面走向长度大于0.9H0(H0为煤层埋深，399.9 m)时，可只设半条走向观测线，故本工作面设置半条走向观测线。观测点按便于观测塌陷区的下沉和能够准确反映工作面下沉情况进行布置，并以观测站为起点布置观测点，首个观测点距开切眼上方地表水平距离288 m。之后每隔25 m设置一个观测点，直至工作面走向推进长度中点之后25 m，以确保观测数据能够监测到下沉量的极值点。观测线全长625 m，各点标记为T1，T2，T3，…，T26。各测点的布置如图2所示，图中矩形部分为61601工作面示意图，三角形中心为观测点位置。走向测线下沉测量工作起始于2019年1月17日，测量周期1年5个月，各测点平均每月测量一次，直至测量完毕。

2 基于概率积分法的地表沉陷预计模型

概率积分法是基于随机介质理论发展而来的岩体移动与变形计算模型，经典的概率积分法表述的地表任意点下沉位移计算公式见式(1)。

(1)

公式(1)表示地表任意点(x，y)在开采区域为(x1-x0)*(y1-y0)的煤层开采后产生的下沉量。Wmax为地表最大下沉量；x0，x1，y0，y1分别是煤层沿倾向开采和走向开采的起止点；r为主要影响半径，可由公式(2)得到。

(2)

式中H为采深；tanβ为主要影响角正切值。

文中使用Python Integrate模块下的Dblquad函数编写基于概率积分法的地表任意点移动与变形模拟程序，设地表为一倾向长945 m，走向长1 392 m的水平平面，并取左下角为坐标原点(0，0)。平面采用矩形网格进行划分，依照61601工作面实测数据确定开采范围和岩移参数，煤层厚23.15 m，下沉系数0.53，主要影响角正切值2.9。地表下倾向开采区域位于345～600 m之间，走向开采区域位于384～1 008 m之间。

模拟结果将保存为记录有地表各节点坐标及下沉量的数据文件。将该文件导入Surfer中，绘制地表开采下沉量等值线图，如图3所示，同时比较计算值与实测值，绘制相应的地表下沉曲线图，如图4所示(部分观测点受采动影响发生破坏，这些点未在图中绘出)。

比照实际观测成果可知，61601工作面地表下沉变形趋势总体上与概率积分法沉陷预计模型相吻合。从图3可看出，采动后地表出现了一个以工作面所在位置上方为中心并向外扩大的类椭圆形沉降区域。地表下沉值较大，达到-13.32 m，等值线在-2.5～-10 m之间密集分布，显示此区间内地表高程存在较大幅度的下移变化，形成了坡度较陡的地表沉降区域边界，下沉量大于-10 m的等值线则较为稀疏，表明地表形成了上宽下窄的“漏斗状”下沉盆地。从模型的岩移参数上看，相比普通煤层的开采[25-26]，特厚煤层综放开采的地表下沉系数较大，边界角值比较小，主要影响角正切值显著偏大，反映出大采高、快速推进的综采放顶煤开采地表下沉速度快，岩体动态变形幅度大、运动剧烈的特征。此外，对比发现实测值具有更为陡峭的地表下沉边界，更为集中的地表下沉区域，实测地表发生沉陷的范围更小，下沉曲线在盆地边界收敛更快，充分地体现出特厚煤层综放开采地表沉陷相比一般地表沉陷规律具有下沉幅度剧烈、变形集中等特点。

概率积分法具有公式简单、计算高效快捷的特点，其只需借助较少的岩移参数，就能计算地表沉陷情况，为研究地质条件简单、岩层达到稳定状态的工作面沉陷状况提供了快速预计和分析的方法。但是概率积分法是一种基于随机介质理论的理想化模型，其不能充分考虑开采区各岩层的空间构型、力学特征的差异对地表移动变形产生的影响；同时作为一种静态地表沉陷预计方法，概率积分法难以直接建立地表移动变化趋势与进尺长度、时间等因素之间的关系等。为提高地表沉陷预计模型的精准性，进一步解释实地观测和概率积分法模型预计结果的形成机理，从动态的角度探究开采沉陷随时间的变化规律，这里讨论采用FLAC3D数值模拟方法对地表沉陷预计模型进行补充。

3 基于FLAC3D的61601工作面地表沉陷规律模拟

3.1 FLAC3D数值模拟模型的建立

根据61601工作面相关地质资料建立几何模型，如图5所示。该模型沿X轴方向全长1 195.4 m，对应开采工作面的倾向方向；沿Y轴方向全长2 000 m，对应开采工作面的走向方向。模型依据各地层实地勘测到的岩土组成情况分别建立了28个地层，并根据地质结构添加部分节理接触面以提高精度，总深度417.8 m。其中开采工作面位于第24岩层(六煤)的下部，开采区域所在位置为：沿X轴方向0～254.6 m，沿Y轴方向0～624 m。模型选用矩形网格进行划分，根据矿区岩层倾角很小这一特征将模型岩层形态设置为水平岩层。

以往的研究表明，直接将岩土力学试验获得的力学参数应用于数值模拟当中往往容易发生错误，通常需取实际力学参数的1/3～1/5作为数值模拟的模拟值，并进行反演分析多次计算方可取得较为准确的模型力学参数[27-28]。根据模拟比例经验值运用反演计算确定岩层模型最优力学参数，并导入模型中，主要岩层模型的力学参数见表1。

表1 主要岩层模型岩性及力学参数

本次模拟基于摩尔-库伦弹塑性模型进行预测计算，分别记录工作面推进169 m(开挖第37天)、244 m(开挖第77天)、342 m(开挖第181天)、496 m(开挖第310天)和624 m(开挖第369天)时地表的走向下沉情况。开挖模型如图6所示，其中2个三角形的位置对应地表下沉观测点的起止范围(T1～T26)，模型下部的矩形代表61601工作面投影，模型开挖进度与61601工作面实际开挖进度保持一致。

3.2 结果及分析

选取工作面进尺169，244，342，496和624 m时的走向地表Z向位移拟合数据，与地表观测站实测下沉量进行比照(表2)，绘制下沉曲线瀑布图，如图7所示(部分观测点受采动影响发生破坏掩埋无法记录，这些点未在图中绘出)。

从图7(e)数据可看出，开采完毕后FLAC3D地表沉陷预计模型的最大下沉量为-13.39m，与实测值(-13.32 m)相差小于5%，可见FLAC3D数值模拟在最大下沉量方面显示出了很好的拟合度。实测的地表最大下沉点位于距观测线起点525 m的T22号观测点附近，而数值模拟的最大下沉点位于地表平面的中心点(T25点)，这可能是因为自然条件下矿区岩体仍存在一定倾角以及岩层间存在断层等不连续结构所致。此外，与模拟值相比，开采初期(图7(a)和(b))实测下沉曲线在拐点后的切线斜率迅速增大，形成更陡峭的下沉曲线，但随着工作面的推进，模拟下沉曲线与实测数据间的差距持续缩小(图7(c)和(d))，至停采线位置时(图7(e))，二者已基本趋于一致。总体而言，FLAC3D对特厚煤层综放开采地表沉陷的数模模拟效果相比经典概率积分法预计模型更为精确，可较为全面地拟合开采过程中地表下沉趋势和最大下沉情况。

从图7可看出，随着工作面的推进，地表岩层也随之发生移动变形，并逐渐形成一个规模远大于采空区的下沉区域。分析地表走向发生明显沉降(>1 m)区域的各测点实测下沉速度随工作面推进的变化情况(图8)可知，观测线一侧的地表在工作面推进169～244 m和496 m处形成2个沉陷速度峰值区域，其中T14和T15测点的最大下沉速度在推进距离169 m附近取得最大值，T15～T18测点最大下沉速度在推进距离244 m附近取得最大值，并在T17点处出现地表最大下沉速度12.72 cm/d，以上各点在244 m后下沉速度开始迅速降低，表明开采推进169～244 m时地表及工作面各上覆岩层处于剧烈的运动当中。采空区中部区域的地表各测点下沉速度则在推进距离496 m左右达到最大值，形成第2个较为明显的峰值区域，但低于首次出现峰值区域的下沉速度。地表走向下沉区域各点的沉降速度均呈现先上升后下降的周期性变化，且总体上距开切眼越近的点，其达到最大沉降速度的时间越早，开采中地表沉降速度最大值所在位置均滞后进尺线一定距离，平均最大下沉速度滞后距为132 m。当工作面推进到169 m处时，地表下沉较为平缓，随后地表开始急剧下沉，至推进长度244 m时最大下沉值达到-7.745 m，地表移动变形异常活跃，此后61601工作面长度中点(推进312 m)之前的下沉边界区各点下沉速度开始放缓并逐渐趋于稳定，中点之后的地表各点下沉速度开始加快，但此后也不断减慢，到开采作业结束前地表各点下沉速度均下降至趋近于零，地表总体移动活跃时间较短。造成这一特征的原因是由于快速推进，使覆岩由下而上传递变形的时间缩短，推进速度越快，上覆岩层越接近整体连续变形，动态变形过程相对缩短，压实效果愈好，从而使移动期缩短。此外，通过现场实测发现，工作面切眼位置两侧下沉盆地边缘出现动态裂缝，裂缝方向与采空区边界方向基本一致，宽度一般200～300 mm，且以工作面切眼位置最为发育，这表明下沉盆地边缘区域处于61601工作面裂缝带中。工作面走向边界角、移动角、裂隙角、超前影响角见表3。综合分析可知，特厚煤层综放开采下，地表沉陷呈现出盆地边界陡峭、变形分布集中，地表下沉速度大且变化幅度大等特征，这是由于工作面快速推进，煤层上覆岩层各层次的下沉速度普遍快速增加，相对悬露的时间减少，使得岩层发生形变、垮落并进而影响地表结构的范围趋于集中。另外，特厚煤层综放开采的采高大，岩体破碎、断裂区域扩大，形成的冒落带、裂缝带相应增高，弯曲带相对减小，也使地表移动变形分布相对集中。

表2 地表观测站部分测点实测下沉量

表3 工作面走向边界角、移动角、裂隙角、超前影响角

运用FLAC3D模拟倾向主断面处地表沉陷随推进长度的变化情况，并记录不同推进长度下地表最大下沉值(表4)。观察可发现，推进长度≤342 m时，地表下沉较不明显，下沉值均低于1 m。当工作面推进到395 m时，地表最大下沉值首次>1 m，开始进入下沉活跃阶段。绘制工作面推进长度在244,342,395,496和624 m处工作面及其周围(420～795 m)地表下沉曲线图(图9)，可知开采结束后倾向主断面地表沉陷同样具有下沉速度大、下沉幅度剧烈、变形分布集中等特厚煤层综放开采特征，且相比走向主断面在距采空区中点上方地表的任意观测点长度之间的倾斜程度更大，形成的下沉盆地边缘也更陡峭。开采初期，采空区范围较小且未到达倾向主断面正上方，地表受影响发生形变的范围及程度相对较小，倾向下沉现象直到工作面推进过半以后才较为明显，且随着工作面的不断推进并最终到达停采线，下沉曲线的整体对称性出现较为显著的下降。倾向主断面的下沉变化持续时间比走向主断面更短，这是因为在工作面推采过程中，地表倾向主断面上各点受采动影响较走向主断面更晚。

表4 倾向地表最大下沉值

采用FLAC3D和Tecplot 360绘制开采区域Z向位移云图和走向主断面Z向位移云图，如图10,11所示。与图7(a)和图7(b)对比可知，工作面推进196 m之前，煤层处于初始采动状态，走向下沉区域呈现上部平缓，下部深陷的“V”字型断面，且下沉点集中分布于最大下沉点两侧100 m的范围内，地表未形成明显的下沉盆地。此后，采空区范围不断扩展，工作面上覆岩层受应力影响发生变形断裂，并向上传递至地表，地表继而快速进入沉降活跃期，活跃期内的移动与变形剧烈且集中，走向地表下沉区域进一步扩大并形成底部开阔的“U”字型断面(图10，图11)，下沉盆地逐渐发育。随着工作面上覆岩层逐渐垮落并充填采空区，地表移动进入衰减阶段，开采至临界开采尺寸时，地表达到充分采动状态，此时开切眼附近的下沉量已基本保持稳定，地表下沉盆地近似椭圆形，底部区域保持平坦并逐渐向两侧扩展。

4 结 论

1)特厚煤层综放开采地表下沉情况基本符合概率积分法预计模型，但其发生沉陷的区域范围更小，沉陷区域边界坡度更大；FLAC3D数值模拟在下沉趋势和下沉量方面与实测值吻合度更高，且能从地表下沉速度、下沉盆地发育状况等动态变化量出发分析地表沉陷规律，对概率积分法地表沉陷预计模型具有较好地补充作用。

2)特厚煤层综放开采具有开采厚度大、推进速度快的特点，其地表沉陷情况相较一般开采呈现出下沉速度快、下沉幅度剧烈、变形分别集中的特征，形成的地表下沉盆地边界陡峭、中部区域扩大，地表最大下沉值达到-13.32 m，下沉系数0.53，主要影响角正切值2.9，主要影响半径137.90 m。

3)地表走向主断面各观测点的沉降速度均先后出现先上升后下降的变化，并在观测线一侧取得2个较为明显的峰值区域，最大下沉速度12.72 cm/d，最大下沉速度滞后距平均132 m。倾向主断面在工作面中点上方地表的任意观测点长度之间的倾斜程度比走向主断面更大，下沉盆地边缘更陡峭，下沉变化持续时间也更短。地表移动变形动态特征方面，工作面推进196 m之前走向下沉区域呈现“V”字型断面，且随着开采逐渐扩大，整体形成底部开阔的“U”字型断面并进一步发育为盆地，达到充分采动状态时，地表下沉盆地近似椭圆形，底部区域平坦，并逐渐向两侧扩展。