王永辉, 冯 园, 李兰兰

(1.河南省地质调查院，郑州 450001；2.河南省城市地质工程技术研究中心，郑州 450001)

松散层指第三纪、第四纪以来尚未固结硬化成岩的疏散沉积层[1]，厚度>100 m可称巨厚松散层[2]。华中、华东、西北及东北大部分矿区都被巨厚松散层覆盖着[3]。巨厚松散层下的开采沉陷问题比一般地质采矿条件更加严重，出现地面下沉系数接近甚至大于1[4]、地面沉陷范围扩展、地表移动剧烈、对土地资源和建(构)筑物危害大等独特现象[3]，对沉陷影响区内的居民生命和财产安全造成了严重威胁。目前，对这一问题的系统性研究仍不多见。因此，对巨厚松散层这种特殊地质采矿条件下的开采沉陷岩土体移动变形特征的分析研究具有重要的现实意义。

1 采区的地质采矿条件

研究矿区的煤层缓倾、开采厚度大、上覆松散层超过100 m，是石炭系—二叠系煤层的典型代表。534E工作面(图1)位于该煤矿一采区，走向长壁后退式采煤，炮采工艺，垮落法顶板管理。

研究区属华北型含煤岩系，地层自上而下分别为：第四系，覆盖全区，上部棕黄色黏土、砂质黏土与黏土质砂、砂等相间成层；中部由黄色黏土、砂质黏土及黏土质砂组成，以黏土为主；下部由灰白色黏土质中、粗砂及砂砾石互层。侏罗系，上部灰绿色粉砂岩为主，夹细砂岩；下部红色粗、中、细粒砂岩。二叠系，上部为下石盒子组(P1x)泥岩；下部为山西组(P1s)灰白色砂岩，为该矿区主要煤系地层，3#煤层厚度在8.8 m左右，稳定可采。石炭系，上部为太原组(C3t)灰色砂岩、细砂岩，矿区主要煤系地层之一；下部本溪组(C2b)石灰岩夹铝质泥岩及铝土岩。奥陶系，煤系地层基底，上部马家沟组(O1m)褐灰色厚层状灰岩、白云质灰岩；下部冶里组(O1y)灰白色白云质灰岩(图2)。

2 地表沉陷变形规律

2.1 概率积分法

概率积分法是一种介于经验方法和理论方法之间的一种预测方法[1]。评价地表岩土体移动变形状态主要有5个常用指标：下沉、曲率、倾斜、水平变形和水平移动等[5]。大量实践应用表明，该方法对地表沉陷变形进行预测分析时，计算结果较为接近实际值[6]，是目前煤矿开采沉陷地表变形计算中应用最多的方法之一[7]。

2.2 预测模型

运用叠加原理，建立采空区走向、倾向主断面及平面任意点的数学模型。

a.采区走向主断面模型(图3)如下

(1)

b.采区倾向主断面模型(图4)如下

(2)

式中：W0(x)、i0(x)、K0(x)、U0(x)、ζ0(x)分别表示走向主断面x处沉陷、倾斜、曲率、水平移动、水平变形值；W0(y)、i0(y)、K0(y)、U0(y)、ζ0(y)分别表示倾向主断面上y处沉陷、倾斜、曲率、水平移动、水平变形值；erf为高斯误差函数；α为煤层倾角；W0为充分采动时的地表最大沉陷值，W0=Mqcosα；M为煤层可采厚度；q为下沉系数；r为煤层主要影响半径；b为水平移动系数；θ0为开采影响传播角；Cxm、Cym分别为走向、倾向采动程度系数；t1、t2分别为下山、上山边界相应参数r1、b1和r2、b2；l为走向计算长度，l=D3-S3-S4；D3为工作面走向长度；S3为走向左拐点偏距；S4为走向右拐点偏距；L为倾向计算长度；D1为工作面倾向长度；S1、S2分别为走向的下山、上山边界拐点偏距。

c.沉陷范围内地面任意点(x,y)移动变形预测数学模型如下

(3)

式中：W0′(x)为倾向充分采动时走向主断面的下沉值；W0′(y)为走向充分采动时倾向主断面的下沉值；i0′、K0′、U0′、ξ0′分别为倾向(走向)达到充分采动时，走向(倾向)主断面上的倾斜、曲率、水平移动、水平变形值；φ为从横坐标x方向逆时针旋转到待求方向的角度；C(x)、C(y)分别为主断面下沉分布系数。

2.3 概率积分法模型修正

巨厚松散层相比基岩属于典型的软弱层，二者物理力学性质相差巨大，所以，将534E工作面的工程地质模型修正为双层模型，即上部松散层，下部基岩层(图5)。在以往的大多数研究中，所采用的是综合角量参数，未对基岩和松散层各自的角量参数进行区分，其预测结果与实际有一定的偏差[8]。本文对概率积分法的数学模型进行了修正，基岩层和松散层具有各自的角量参数(图6)。

修正后的影响半径的计算

(4)

(5)

式中：δ0为基岩走向边界角；φ为松散层移动角。

2.4 概率积分法预测参数

概率积分法中参与计算地表沉陷变形的主要有6个参数：下沉系数q、水平移动系数b、主要影响角正切tanβ、开采影响传播角θ0、采动程度系数C和拐点偏距S等[7]。

a.下沉系数q的大小受地质因素和采矿因素的综合影响。从宏观上，岩体硬度越大，q值越小；但二者之间的定量关系及其他因素的影响，还处于研究阶段[9]。

b.水平移动系数b的取值与矿区具体的地质采矿条件相关，包括：松散层厚度、采深、岩层倾角等因素[10]。松散层本身具有流变特性，使得水平移动系数增大。

c.主要影响半径r与主要影响角正切tanβ。

在修正模型中，tanβ是分配到各层中的开采深度与模型分层中的主要影响半径r的比值。tanβ的值主要受岩性的影响，对采深的变化反应不明显[1]。求取rN时，岩性越硬，tanβ越小，r就越大。

d.开采影响传播角θ0主要与岩性、煤层倾角有关。与倾角的关系为[6]

θ0=90°-kα

(6)

式中：k值与岩性有关，覆岩越坚硬k值越大，软弱岩层取0.5～0.6，经计算可得开采影响传播角θ0=87°。

e.拐点偏移距S产生的原因是采空区边界附近的覆岩受煤柱的支撑，起到了“悬臂梁”的作用，使得计算边界和开采边界之间产生了一定相隔距离。拐点偏移距的大小主要受开采深度、采动程度、顶板控制方法及岩性的影响[11]。

f.采动程度系数Cxm和Cym最大下沉值点坐标可用(7)式来计算[1]

(7)

式中D1=129 m，得ym=42.5 m，计算可得W0(ym)=3 323 mm，地表最大沉陷值在充分采动时W0=8 184 mm，倾向采动系数Cym=0.41。走向若为有限开采，最大沉陷点坐标在l/2，走向主断面模型中计算W0(l/2)=8 154 mm，则Cxm≈1.0。

此次预测使用的参数是对鲍店1308观测站的数据基于最小二乘原理进行拟合得到的(表1)。

表1 概率积分法预测模型实测参数Table 1 Measured parameters by probabilistic integration method prediction model

2.5 地表沉陷分布规律预测

2.5.1 建立坐标系

根据前述各项参数建立沉陷区地表的xoy直角坐标系如图7所示。

2.5.2 主断面移动变形分析

基于修正的预测模型(1)式、(2)式得到534E工作面开采后走向及倾向主断面上各点移动变形值。下面就下沉、地表移动和地表变形等主要指标进行分析。

a.下沉分析

在走向上，采动系数达到1，在M点附近出现一段平底，移动盆地略显盆形(图8)；在倾向上，由于远未达到充分采动，移动盆地呈尖底(图9)。

b.水平移动分析

在水平移动预测曲线上(图10)，水平移动值在下山拐点O1(1 791 mm)和上山拐点O2(-1 467 mm)达到了最大值，所有点的水平移动均指向沉陷中心点。

如图11所示，由于岩层倾角的影响，双层介质的运动形式略有差异。基岩的移动Sj指向煤层法向，其水平移动全部指向上山方向；松散层的水平移动指向沉陷中心，在下山方向，基岩和松散层水平移动的方向一致，地表点的水平移动值为二者叠加；在上山方向，二者水平移动方向相反，地表点的水平移动值为二者之差，所以上山的水平移动小于下山。这种影响会随着煤层倾角的增大而增大、松散层厚度的增大而减小。

如图12所示，在移动盆地的边缘，下沉值(实线)小于水平移动值(虚线)，使得水平移动的范围大于沉陷的范围，这也造成了移动盆地的扩展。

c.水平变形分析

图13中，下山点A至拐点O1水平变形大于0，地表呈拉伸状态，受双向拉应力影响，在地面容易产生拉裂缝；从拐点至最大下沉点水平变形小于0，地表呈压缩状态，受双向压应力影响，先前出现的拉张裂缝会不同程度闭合。

2.5.3 地表移动变形分析

通过修正的模型(3)式，得到沉陷区内地表任一点的下沉值W(x,y)。以10 mm的下沉值作为地表移动盆地起始点，描绘地表下沉等值线图(图14)。

从图14可以看出沉陷盆地以采空区为轴具有良好的对称性。巨厚松散层的边界角和移动角较小，影响半径较大；同时，因松散层具有流变性，在下沉的同时又以水平流动的形式填充下部的沉陷空间，水平移动系数较大，使得地表移动盆地边缘收敛缓慢，范围远超采空区的范围。松散层较软，下沉时不出现层离裂缝，呈整体性下沉，在压应力、导水裂隙及抽水等因素的综合影响下，黏土层失水压缩固结，使下沉系数和下沉量都较大。

3 地下岩土体移动变形数值模拟

地下岩土体变形移动是开采沉陷问题的核心，是地面沉陷的原因。

3.1 FLAC3D数值模型

以巨厚松散层为重点，对砂砾石、黏土、砂土和基岩进行性质区分，对薄且无控制作用的岩层进行有效合并，共划分了48 000个单元，52 173个节点。确定模型边界条件：①模型四周边界施加水平约束，水平位移为0，xdisp=ydisp=0；②模型底部边界固定，无垂直和水平位移，即zdisp=0；③顶部按照实际情况，无约束，即为自由边界。由于岩体和松散层属于弹塑性地质体，经历弹性变形达到屈服极限后会产生较大的塑性流动变形，故选择莫尔-库伦塑性模型。

3.2 岩体力学参数的选取

参数的选取至关重要，所需的参数包括[12]：密度、体积模量、剪切模量、内聚力、摩擦角和抗拉强度。第四系和煤层的各项物理-力学参数来自现场取样的测试结果；上覆岩层的各项参数中，根据实测的抗压强度，按照经验关系来求取抗拉强度[13]；参照该煤矿以西4 km的杨村矿区实测的岩层物理-力学参数，同时结合矿区地表移动变形的实测数据对参数进行一定的修正，确定的参数如表2所示。

表2 模型物理-力学参数Table 2 Physical-mechanical parameters used in model

3.3 模拟结果分析

应力-应变平衡之后地表最大沉陷量3.26 m，沉陷范围沿采空区呈椭圆形向外扩散，沉陷主要集中在采空区上方，这与概率积分法的计算结果基本吻合(图15)。最大水平移动量发生在采空区边界附近，即概率积分法中的计算边界(图16)，验证了此前关于地表移动变形规律的分析。

从图17、图18可以看出，岩土体移动变形从采空区附近的基岩至地表遵循先减小再增大、再减小的规律，这主要取决于岩土体的性质。根据关键层理论，煤层上覆的砂岩层在沉陷变形过程

中起控制作用，松散层相当于砂岩层上的荷载。由于砂岩层厚度较小，力学性质较差，在上覆巨厚松散层的垂向压力和自重应力下，发生垮落沉陷，过程中伴有一定的碎胀性和剪胀性，部分采空区被填充，从采空区至基岩上表面，最大下沉值会呈现减小的趋势。

砾石层中，下伏基岩移动变形造成的扰动使砾石颗粒之间的空隙迅速增大，沉降量大幅减小；黏土层中，下伏基岩和砾石在移动变形后具有良好的导水性，黏土层失水固结，且黏土层厚度较大，固结沉降量较大，完全抵消了基岩和砾石层沉降的减小的量；砂土层中，巨厚松散层的缓冲作用开始显现，砂土层移动过程中变得疏松，地表沉陷量略有减小。整体而言，基岩上表面至松散层表面最大下沉值呈增大的趋势，局部会因为剪胀作用有减小趋势。

松散层中某点的变形量包括基岩下沉引起的土体塌陷、应力集中造成的剪胀变形和剪缩变形、含水层失水造成的黏土层主固结沉降以及土体内附加应力产生的次固结沉降。固结沉降远大于剪胀变形，所以，地表的沉陷量大于基岩的沉陷量。

4 地表沉陷实测分析

在该煤矿1#采区闭坑后，对地表塌陷区进行了1∶1000的地形测绘和沉降监测工作，实测沉陷值和地裂缝分布如图19所示。同时，通过与20世纪60年代煤矿开采前的地形图对比，确定该区域内原始地面海拔高度为40.0 m。

沉陷区东部正处于刚刚放顶之后，地表处于持续沉陷过程中，受拉应力影响地表遍布地裂缝。在1#采区的中部和西部，属早期开采沉陷区域，地表已基本稳沉，由图可见，沉陷等值线分布也与前述预测结果基本吻合。地表最大沉陷值>8.0 m，<8.5 m，这与解析法分析充分采动情况下的倾向和走向最大沉陷值(W0=8 184 mm，W0(l/2)=8 154 mm)是基本一致的。由此可见，解析法建立的双层预测模型和选取的参数是正确的，预测结果是可信的。基于解析法物理模型建立的数值分析模型，在对松散地层进行了进一步细分的基础上，结合实测的物理、力学参数，对各类岩土体沉陷过程中移动变形规律的分析也是可信的。

5 结 论

巨厚松散层的沉陷特征明显有别于一般地质采矿条件。移动盆地具有良好的对称性。下沉系数、最大下沉值、水平移动系数、水平移动值和沉陷盆地的范围都比一般地质采矿条件要大，盆地边缘的水平移动值大于下沉值。而黏土层的压缩固结是造成沉降量增大的根本原因。在沉陷过程中，松散层同时扮演着荷载、固结沉降和缓冲系统的三重角色。这些特殊的移动变形规律取决于巨厚松散层颗粒疏散及富含地下水。