概率积分法在煤矿采空区地表变形评价中的应用

【摘 要】煤矿开采沉陷是由于煤矿地下煤层的开采导致的一种主要的地质灾害类型，对土地的破坏非常严重。本文通过对宁东某煤矿地质采矿条件的研究，利用概率积分法计算评价开采引起的地表变形及稳定性，为矿山企业安全生产及地质灾害防治提供科学参考依据。

【关键词】开采沉陷；地质灾害；概率积分法；地表变形；MathCAD软件

0 引言

煤炭从地下采出后，其上方覆盖岩层失去支撑，岩体内部的原有应力平衡状态受到破坏，引起岩体内应力的重新分布，使采空区周围的岩体产生位移、变形直至破坏。当采空区面积达到一定的范围后，引起采空区上方地表产生移动与变形，从而引发一种主要地质灾害类型，对矿区地质环境的破坏非常严重[1-2]，其主要表现为：地面塌陷及地裂缝使得地表耕地破坏；干旱山区地表水系受到破坏；平原地区大片地表移动形成盆地面积水；受采动影响的山地产生滑坡，危及工业及民用建筑物和生活生产设施的安全；道路塌陷、房屋变形破坏等[3-4]。

开采沉陷的预计方法主要有概率积分法、负指数函数法、典型曲线法和数值计算法（有限单元法、边界单元法和离散单元法等）。概率积分法源于上世纪50年代波兰学者李特威尼申引入岩层移动领域的随机介质理论，后由刘宝琛、廖国华等学者引入国内，经过几十年的研究，目前已成为我国较成熟的、应用最为广泛的预计方法之一[5]。

本文针对宁东某煤矿开发利用方案，通过研究煤矿采空塌陷与煤层厚度、倾角、采深、采厚、上覆岩层性质及开采方式等，建立计算模型并采用概率积分法计算走向及倾斜主断面上地表移动与变形，为该矿的安全生产和地质灾害防治提供可靠的科学依据。

1 地表移动变形预计[5-8]

1.1 半无限开时主断面上地表移动与变形的预计

地表移动观测的实践表明：当开采宽度达到相当大以后，再继续增加开采宽度，对地表下沉的最大值和停采线上方的地表变形几乎没有影响。将这种开采情况数学抽象，即认为停采线一侧煤层未被采动，而采区另一侧的煤层已全部采空，这就是所谓的半无限开采。

1.1.1 走向主断面

W（x）=■[erf（■）+1] T（x）=T■e■

K（x）=2？仔■（-■）e■ U（x）=bW■e■

ε（x）=2？仔b■（-■）e■ （1）

式中：W（x）、T（x）、K（x）、U（x）、ε（x）分别为主断面上下沉、倾斜、曲率、水平移动和水平变形值；W■、T■、K■、U■、ε■分别为最大下沉、最大倾斜、最大曲率、最大水平移动和最大水平变形，可按如下公式计算： W■=qmcosα、T■=W■/r、K■=1.52W■/r2、U■=bW■、ε■=1.52bW■/r；x――走向主断面上地表点至下沉曲线拐点的距离；e――自然对数的底；erf（x）――误差积分，erf（x）=■■e■dt；r-――主要影响半径；q――下沉系数；m――采厚；α――煤层倾角；b――水平移动系数，下同。

1.1.2 倾斜主断面

倾斜主断面上地表移动与变形的预计公式与式（1）相同，但需以y代替各式中的x。当预计下边界半个盆地移动与变形时，以下山方向的主要影响半径r1代替各式中的r；预计上边界半个盆地移动与变形时，用上山方向的主要影响半径r2代替各式中的r。其中，r1、r2按下式计算：

r1=H1/tanβ r2=H2/tanβ

式中，H1――工作面下边界的开采深度；H2――工作面上边界的开采深度。

1.2 有限开采时主断面上地表移动与变形的预计

有限开采是指回采面积未能使地表出现充分移动的开采。

1.2.1 走向主断面

W（x）=Cy[W3（x）-W4（x-l）]

T（x）=Cy[T3（x）-T4（x-l）]

K（x）=Cy[K3（x）-K4（x-l）]

U（x）=Cy[U3（x）-U4（x-l）]

ε（x）=Cy[ε3（x）-ε4（x-l）] （2）

式中：W3（x）、T3（x）、K3（x）、U3（x）、ε3（x）分别为工作面从左边界向右侧半无限开采的下沉、倾斜、曲率、水平移动和水平变形；W4（x-l）、T4（x-l）、K4（x-l）、U4（x-l）、ε4（x-l）分别为工作面从右边界向左侧半无限开采的下沉、倾斜、曲率、水平移动和水平变形；Cy――倾斜方向采动程度系数，l――工作面走向方向的计算长度，L――工作面倾斜方向的计算长度，分别按下式计算：

Cy =■erf■+1-■erf■+1

l=D■-2S■ L=（D■-S■-S■）■

其中，D■、D■――工作面倾斜和走向的长度；S■、S■、S2――走向、下山和上山方向的拐点偏移距；θ――开采影响传播角，下同。

1.2.2 倾斜主断面

W（y）=Cx[W1（y）-W2（y-L）]

T（y）=Cx[T1（y）-T2（y-L）]

K（y）=Cx[K1（y）-K2（y-L）]

U（y）=Cx[U1（y）-U2（y-L）]

ε（y）=Cx[ε1（y）-ε2（y-L）] （3）

式中：W1（y）、T1（y）、K1（y）、U1（y）、ε1（y）分别为工作面从下边界向上山方向半无限开采的下沉、倾斜、曲率、水平移动和水平变形；W2（y-L）、T2（y-L）、K2（y-L）、U2（y-L）、ε2（y-L）分别为工作面从右边界向左侧半无限开采的下沉、倾斜、曲率、水平移动和水平变形；Cx――走向方向采动程度系数，按下式计算：

Cx =■erf■+1-■erf■+1

1.3 地表移动盆地内任意点、任意方向移动与变形的预计

W（x，y）=■W（x）W（y）

T（x，y）=■[T（x）W（y）cosφ+T（y）W（x）sinφ]

K（x，y）=■[K（x）W（y）cos2φ+K（y）W（x）sin2φ+T（x）T（y）sin2φ]

U（x，y）=■[U（x）U（y）cosφ+U（y）U（x）sinφ]

ε（x，y）=■[ε（x）ε（y）cos2φ+ε（y）W（x）sin2φ+■（T（y）U（X）+T（x）U（y））sin2φ] （4）

式中：φ――从横坐标轴x轴方向逆时针旋转到给定的计算方向的角度（φ的区间为0-360°），其它符号同上。

2 概率积分法预计关键技术问题处理

2.1 多工作面开采影响的预计

首先计算出单个工作面开采引起地表预计点的下沉，然后根据叠加原理进行叠加计算获得多工作面开采的联合影响。

2.2 任意形状工作面开采影响的预计

将任意形状工作沿煤层走向或倾向划分成若干个狭长条带，近似看成矩形工作面，用公式（4）求出每个小矩形开采对地表任意点沿指定方向φ的移动和变形值，再将所有小矩形开采计算结果进行相应的叠加。据童立元等人[7]研究认为，当将采空区分割成0.1hav（hav为工作面平均采深，m）矩形单元时，精确解和近似解平均值的相对误差小于5%，按此分割叠加求解完全能满足工程实际的精度要求。

2.3 多层重复采动

重复采动是指岩层和地表已经受过一次开采的影响而产生移动、变形和破坏之后，再一次（或多次）经受开采（开采下部煤层，或同一煤层的下一个工作面）的影响，使得岩层和地表又一次受到采动。重复采动采空区地表和岩层移动过程加剧，通过下沉系数的增大进行处理，可按以下公式计算：

q复1=（1+a）q初 q复2=（1+a）q复1 （5）

式中：a―下沉“活化”系数，参考文献6附表5-5；q初、q复1、q复2分别为初采、第一次和第二次重复开采的下沉系数。

重复采动岩层和地表移动的预计可以用叠加法进行，步骤如下：①首先预计第一次煤开采的岩层与地表移动变形值；②按相同的方法和公式，用相应的重复采动时的预计参数预计各次采动移动与变形值；③将相同点沿相同方向初次开采和重复开采所得的同名移动与变形相加（代数和），即得各次开采总的移动与变形值。

2.4 预计参数的选取

影响地表沉陷的主要因素有采深H、采厚m、煤层倾角α、开采范围及地质构造等。这些因素反映到岩层移动和地表下沉的形式为岩移参数。主要岩移参数有下沉系数q、水平移动系数b、主要影响角移动角正切tanβ、下沉盆地的拐点偏移距S和开采影响角θ。原则上讲，这些参数应根据矿区实测资料确定，以提高预计的精度。如无实测资料时，可根据经验公式计算或统计数据确定[6]。

2.5 停采后地表变形总延续时间的估算

在无实测资料或实测表明采空区地表沉降尚未终止时，可按文献6中公式估算变形终止时间，即地表总移动期限计算公式：

tm=0.25hav （6）

计算的起算时间从矿区停采开始。

3 工程实例

3.1 工程概况

拟建矿山工程位于宁东矿区灵盐台地低缓的黄土丘陵地带，海拔高程+1213～+1230m，高差17m。矿区含煤地层为二叠系山西组，该煤层的上覆岩层分别为古近系泥岩、白垩系灰色砾岩夹薄层粉砂岩及二叠系系中上统灰紫、紫色粉砂岩、细砂岩。从整个岩层的岩性分析来看，该地区煤系地层的上覆岩层为中等硬度岩石。地表有一层2～10m表土层。

该矿主采煤层是一、三、五煤，开采方法为走向长壁全部陷落法，相关参数见表1、2。

表1 煤层地质开采参数

表2 地表移动主要参数

3.2 地表变形计算评价

由于采煤沉陷预计计算量大而繁琐，特别是误差函数erf（x）在概率积分法中有着最为基础的作用。因此，对其准确而又快速地求解对地表移动与变形值的计算非常重要。但e■没有原函数，故不能采用牛顿-莱布尼兹公式求解[9]，只能做近似计算。经过对常见计算工具Matlab、MathCAD等数学计算软件的对比，笔者选择了MathCAD14，主要是该软件集文本编辑、数学计算、图形绘制于一体，能够以手写格式直观地进行各种数值计算、矩阵计算、方程（含微分方程）或方程组的求解，并且内嵌数百个常用数学函数，其中就包括误差函数erf（x）。根据计算可得，回采工作面倾斜方向为非充分采动（见图1），走向方向为充分采动（见图2）。

（a）沿倾斜主断面下沉值（m） （b）沿倾斜主断面倾斜值（103mm/m）

（c）沿倾斜主断面曲率值（1/m） （d）沿倾斜主断面水平移动值（m）

（e）沿倾斜主断面水平变形值（103mm/m）

图1 沿倾斜主断面移动变形（非充分采动）

（a）沿走向主断面下沉值（m） （b）沿走向主断面倾斜值（103mm/m）

（c）沿走向主断面曲率值（1/m） （d）沿走向主断面水平移动值（m）

（e）沿走向主断面水平变形值（103mm/m）

图2 沿走向主断面移动变形（充分采动）

3.3 停采后地表变形总延续时间的估算

回采工作面平均采深hav=300m，运用公式（6）计算的tm=750天≈2.1年，即停采后约2.1年地表变形基本稳定。

4 结论

1）利用概率积分法对地表变形动态评价迅速准确，给定量评价提供了科学的依据。特别是，通过MathCAD数学软件的应用更是可以简化大量的计算、绘制各类变形曲线工作，使得技术人员能专心研究实际工程问题。

2）通过应用概率积分法对宁东某煤矿采空区地表变形计算可知，在缓倾角（小于25°）多煤层重复开采情况下，回采工作面倾斜方向为非充分采动，走向方向为充分采动。

3）地表移动变形计算表明，采空区上方最大下沉值W为4180 mm～5740mm；地表最大倾斜值（T）为30.3～40.6mm/m，最大曲率值K为（0.461～0.850）×10-3/m，最大水平移动值（U）为1850～2060mm，最大水平变形值（ε）为22.8～108.8mm/m。

4）停采后地表变形总延续时间估算结果表明，在评价区地质条件及开采方式下，采煤后约2.1年地表变形可达到稳定状态，经工程处理后土地可作为耕地或建设用地使用。

【参考文献】

[1]谢东海，冯涛，袁坚，等.建筑物下采煤的地表沉陷预测[J].西部探矿工程，2008，4：103-106.

[2]杨梅忠，任秀芳，于远祥.概率积分法在煤矿采空区地表变形动态评价中的应用[J].西安科技大学学报，2007，27（1）：39-42.

[3]陈祥恩，李德海，勾攀峰.巨厚松散层下开采及地表移动[M].徐州：中国矿业大学出版社，2001.

[4]中国科学院，编.矿山地质环境保护规定与矿山自然生态环境保护及治理标准实施手册[S].中国地质出版社，2009.

[5]孙文华，主编.三下采煤新技术应用与煤柱留设及压煤开采规程实用手册[S].北京：中国煤炭出版社，2005，8.

[6]建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程[M].北京：煤炭工业出版社，2000，6.

[7]童立元，刘松玉，邱钰，编.高速公路下伏采空区危害性评价与处治技术[M].东南大学出版社，2006，8.

[8]何国清，杨伦，等，编.矿山开采沉陷学[M].中国矿业大学出版社，1994.

[9]同济大学应用数学系.高等数学[M].北京：高等教育出版社，2002.