CN110135030A - 一种采空区地表沉降的预测方法 - Google Patents

Abstract

本发明涉及一种采空区地表沉降的预测方法。首先拟定采空区的采深、采深采厚比、采宽采深比、上覆岩性四因素以及各自四个水平，采用正交实验设计方法对四因素以及四水平按正交表进行设计，得到16种计算工况。采用有限元计算方法，获得16种工况对应的地表沉降数据。运用多种合适的预测模型，对地表沉降数据进行拟合，拟合前将所有原始数据做归一化处理，避免了数据不同量纲和数量级对模型的扰动。对不同函数拟合结果进行误差分析，最终优选二次函数模型为地表沉降预测模型，并给出预测模型中各参数的取值。本发明所述多因素采空区地表沉降预测模型可以较为迅速可靠的在采空区开挖早期即可初步预测出采空区地表沉降的一个具体数值。

Description

一种采空区地表沉降的预测方法

技术领域

本发明涉及属于地下工程技术领域，特别涉及一种采空区地表沉降的预测方法。

背景技术

国外很早就开始了采空区地表沉降的研究和观测工作，尤其是早期的前苏联、波兰、德国、加拿大和美国等采矿业较为发达的国家对开采沉陷理论和实践都进行了较为深入的研究。早在1867年德国工程师A.Schulz就提出了破裂角和保护地表必要支柱尺寸的观点，在此之后，各国均开始对采空区沉降进行测量、研究。1930年，英国开始了对采空区移动变形的观测并于1950年发现了观测和地表变形之间关系的重要性，建立了不同采动程度下的下沉系数修正体系。学者K.Wardell改进了沉降观测的方法，并对开采沉降的理论方法做出了巨大的贡献。除此之外，还有大家较为熟知的随机介质理论法就是由波兰学者李特维尼申于20世纪50年代创立，并将其引申到岩层移动变形的研究中。

我国对于采空区沉降预测的研究由于历史原因，起步较晚，从上世纪五十年代开始，中国矿业大学、中国科学技术研究院等所及企业，从理论和实践两个方面开始深入的研究采空区引起的地表变形规律，其研究的重点主要集中在残余沉降的检测、残余沉降的分析以及数值模拟预测采空区等这几个方面。后续的十年间，我国开始了相似材料模拟研究；同时，采空区引起的地表沉降规律可以借助各种数学描述语言绘制出地表下沉的曲线、地表水平位移的曲线，这样更形象的展示出变化的过程；在数值模拟方法的应用方面，有限元法、离散单元法、边界元法等计算方法运用的非常之广泛。其中极具代表性的人物就是中国工程院院士刘宝琛。刘宝琛院士长期从事岩石力学有关方面研究，形成了独树一帜的开采影响下地表移动级变形计算方法并开发了系列软件程序，创建发展了时空统一随机介质理论，并将其应用于分析研究地表移动基本规律，对隧道、建筑物下、河下及铁路下开采地表保护工程等方面。

现有的采空区地表表现的计算方法一般有实验法、模拟数值仿真法、概率积分法等。但现有的计算方法实施起来都较为复杂不便，需测量采空区有效参数较多，计算量较大且时间和经济成本较高并且考虑因素较为单一，并没有结合多个因素综合考虑对采空区地表沉降的影响。无法快速便捷的初步计算出一个采空区地表沉降的具体数值。

发明内容

本发明的目的是建立一个在采空区开挖前期，即可根据采空区的采深、采深采厚比、采宽采深比和上覆岩性四个因素计算得出采空区地表沉降的计算模型。并拥有快速便捷、结合多因素考虑采空区地表沉降、预测值较为准确、应用采空区工况较广泛、对所需采空区有效参数较少等一系列优点。较为迅速可靠的在采空区开挖早期即可初步预测出采空区地表沉降的一个具体数值。

本发明采用以下技术方案：

一种采空区地表沉降的预测方法，所述方法包括以下步骤:

步骤一：拟定采空区的采深、采深采厚比、采宽采深比、上覆岩性四因素以及各自四个水平；

步骤二：采用正交实验设计方法对步骤一四因素以及各自四水平按正交表进行设计，得出四因素四水平正交实验表及16种计算工况；

步骤三：运用ANSYS有限元仿真软件对上述工况进行仿真计算，得出各工况的地表沉降数据；

步骤四：采用二次函数拟合，得到采空区地表沉降预测模型；

步骤五：根据模型预测得出模型中各待定参数的取值。

优选地，所述步骤一采深的1、2、3、4水平分别对应数值为：50、100、150、200；采深采厚比的1、2、3、4水平分别对应数值为：20、30、40、50；采宽采深比的1、2、3、4水平分别对应数值为：0.8、1、1.2、1.4；上覆岩性的1、2、3、4水平分别对应数值为：1、0.95、0.9、0.85。

进一步优选地，所述步骤三仿真计算具体包括以下步骤：

3.1地质参数取值：

在ANSYS程序中，采空区地表变形仿真计算需要输入参数的选取如表1:

表1仿真计算岩土力学参数取值

3.2、边界条件的处理和计算范围选择：

(1)边界条件的处理：

位移边界为：底部边界的水平与垂直位移固定约束；左右边界的水平位移固定约束，垂直位移不约束；上边界为地面线，用自由边界；

(2)计算范围的选择：以采空区为中心区域，四周范围取2-4倍采宽；

3.3、仿真计算模型：利用ANSYS建立不同工况的上覆岩土层和采空区的三维地质模型；

3.4、利用3.3得到的地质模型进行地表变形计算，得到地表沉降数据。

更进一步优选地，所述步骤三各工况的地表沉降数据为表2：

表2采空区地表沉降数据

工况	1	2	3	4	5	6	7	8
									W(mm)	0.85	0.70	0.88	0.94	4.91	2.31	5.86	3.47
工况	9	10	11	12	13	14	15	16
									W(mm)	15.90	15.92	4.63	6.04	39.29	23.44	10.71	6.17

进一步优选地，所述步骤四二次函数拟合的拟合公式为

其中：x₁,x₂,x₃,x₄分别代表采深，采深采厚比，采宽采深比，上覆岩性。

更进一步优选地，所述步骤四二次函数拟合前将采深、采深采厚比数据按公式

做归一化处理，数据归化到区间[0.1,0.9]，式中，x_i’-归一化处理后得到的数据；x_max-采深、采深采厚比对应的最大数值；x_min-采深、采深采厚比对应的最小数值；x_i-采深、采深采厚比需要做归一化处理的数据。

更进一步优选地，所述步骤四沉降预测模型为

本发明具有以下有益效果：

1.本发明综合考虑了采深、采深采厚比、采宽采深比、上覆岩性多因素得出多因素采空区地表沉降预测模型较单因素来说更加准确、全面。

2.本发明所述多因素采空区地表沉降预测模型可以较为迅速可靠的在采空区开挖早期即可初步预测出采空区地表沉降的一个具体数值。

3.本发明所述多因素采空区地表沉降预测模型具有一定普遍性，可以较为普遍的运用于各类采空区工况。

4.经检验，本发明所述多因素采空区地表沉降预测模型误差较小，可以较为准确迅速的在开挖采空区前期计算出一个采空区地表沉降的具体值。

5.拟合前将所有原始数据做归一化处理，避免了数据不同量纲和数量级对模型的扰动。

附图说明

图1为采空区影响因素示意图；

图2地质结构三维有限元模型；

图3工况1采空区有限元模型及总位移云图；

图4拟合结果散点图。

具体实施方式

下面结合具体实施例，对本发明作进一步详细的阐述，但本发明的实施方式并不局限于实施例表示的范围。这些实施例仅用于说明本发明，而非用于限制本发明的范围。

实施例1

采空区地表变形ANSYS有限元计算原理

1、地表变形影响因素分析

采空区等特殊地质结构对地表稳定性的影响除了受采空区本身的几何尺寸，开采深度等因素影响以外，还受到复杂的上覆岩性、松散层厚度、地形、地质构造因素以及开采方法、采空区布置方式等诸多因素影响(图1)。本发明通过建立简单的ANSYS模型，对上述复杂影响因素分别进行模拟，分析其对地表稳定性的影响。

(1)开采厚度的影响

开采厚度对上覆岩层及地表的沉陷过程的性质有重要的影响。开采厚度越大，冒落带、导水裂缝带的高度越大，地表移动变形值也越大，移动过程表现得越剧烈，因此移动和变形值与开采厚度成正比。

(2)开采深度的影响

同样条件下的采空区，当其拉伸变形超过岩石的允许抗拉强度时，直接顶板及其上部的部分岩层便与整体分开，破碎成大小不一的岩块，无规律地充填采空区。直接顶板岩层垮落并充填采空区后，由于破碎其体积增大，致使其上部的岩层移动逐渐减弱，埋深越深其对地表的影响越小。

(3)开采宽度的影响

采空区尺寸的大小可影响到地表的充分采动程度。当工作面宽度小于某一极限宽度时，受采动影响的最上部基岩控制层发生最大可能的弯曲，但不出现断裂，该控制层对上覆岩层起着支撑作用，冲积层基本上是随顶部控制层的弯曲而下沉，因而，此时将表现为地面下沉量较小，下沉盆地较缓。但当工作面宽度超过其极限宽度时，控制层因拉应力超过其抗拉强度而断裂，随着下覆岩层的冒落而下沉，地面下沉量将会明显增大。

(4)综合岩性的影响

当矿层采出以后，在采出空间附近形成了三个移动特征区，即充分采动区、最大弯曲区、岩石压缩区。无论是处在哪个区域内的岩层，沿层面方向均承受拉伸变形，最终在竖向方向形成垮落带、裂隙带(离层带)和弯曲下沉带达到新的平衡。由岩性的差异带来的是岩层物理力学参数的变化，由于坚硬岩层的抗拉伸能力和抗压缩能力强于软弱岩层，因此，坚硬岩性的区域比软弱岩性的区域完成三个移动特征区、最终在竖向形成“三带”的时间过程要长，也就是说，同样埋藏条件下和开采条件，综合岩性较硬相对于综合岩性较软的地表移动的时间要长，趋稳慢。另外，当上覆岩层为坚硬岩层时，在岩层的破坏和移动过程中，裂隙带内岩层除产生垂直层面的裂隙外，还产生大量顺层面离开的离层裂隙，而上覆岩体为软岩层时，在岩层移动过程中，上覆岩层不易产生离层裂隙，因此在其它条件相同的情况下，坚硬岩层的地表沉降值小于软弱岩层的地表沉降值。

基于探地雷达的探测数据，分析沿线区域地质构造情况。应用岩土力学方法建立力学模型，利用有限差分法对采动区应力分布进行数值模拟仿真，研究采空区上覆岩层及地表的移动规律及变形特征，分析开采深度、开采厚度、开采宽度、采空区上覆岩性等因素对地表变形的影响程度，并对不同的影响因素进行敏感度分析。

结合探地雷达图像，根据反射波组的波形与强度特征，通过同相轴的追踪，确定反射波组的地质含义，构筑地质——地球物理解释剖面。通过探地雷达剖面图谱的波形、振幅和频率特征，确定不同地质结构的位置、走向及空间分布。并基于不同地质结构的空间分布情况，借助岩体力学理论计算特殊地质结构的应力分布，分析地质结构稳定性，总结地下开采引起的地表移动变形时空规律。

2、基于有限元理论的地表变形数值模拟分析

(1)有限元理论的基本原理

有限元法最突出地优点就是适合处理非线性、非均质和复杂边界等问题，而在采空区围岩体地应力变形分析中就恰恰存在这些难题，因此很适合用有限元法。有限元法是用有限个单元体所构成的离散化结构，代替原来的连续体结构，来分析应力变形。它的求解过程一般有以下几个步骤：

1)研究区域的离散化

离散化就是将所研究的区域划分成有限个大小的单元体，并在单元体的指定点设置节点，把相邻单元体在节点处连接起来组成单元的集合体，以代替所研究问题的区域，并以所离散单元节点处的位移作为基本未知量。

2)选择位移模式

由于采用节点位移为基本未知量，因此需要用节点位移表示单元体的位移。为此，必须对单元中位移分布作出一定假定，一般假定位移是坐标的某种简单函数，这种函数称为位移模式或位移函数。根据所选定的位移模式，即可导出用节点位移表示单元内任意一点位移的关系式，其矩阵形式为

{f}＝[N]{U}e 2-1

式中，{f}为单元内任一点的位移阵列；{U}为单元节点的位移阵列；[N]为形函数矩阵，其元素是位置坐标函数。

3)单元分析

将位移模式代入几何方程，可导出用节点位移表示的单元应变公式

{ε}＝[B]{U}e 2-2

式中，[B]称为应变矩阵。

利用物理方程，由以上应变表达式导出用节点位移表示的单元应力计算公式

{σ}＝[G]{U} 2-3

式中，[G]为应力矩阵。

利用虚功原理建立作用于单元上的节点力和就节点位移之间的关系，即单元刚度方程

{F}＝[K]{U} 2-4

式中，[K]为单元刚度矩阵

4)计算等效节点荷载

将研究区域离散化后，即假定力是通过节点从一个单元传到另一个单元的，但作为实际的连续体区域，力是从单元的公共边界上进行传递的。因此，这种作用在单元边界上的表面力以及作用于单元上的体积力、集中力等都需要等效的移置到节点上，也就是用等效的节点荷载来代替作用在单元上的力。这种力的移置必须遵循静力等效或虚功等效原则业。

5)集合所有单元的刚度矩阵，建立整个结构的平衡方程

此过程包括两方面的内容：一是各个单元的刚度矩阵集合成描述结构平衡条件的整体刚度矩阵：二是将作用于各单元的等效节点列阵集合成总的荷载列阵，最常用的集合刚度矩阵方法是直接刚度法。一般来说，集合所依据的理由是要求所有相邻的单元在公共节点处的位移相等，于是得到以形函数矩阵[N]总体刚度矩阵[K]、荷载列阵{R}以及节点位移列阵{U}表示的整个结构的平衡方程

N[K]{U}＝{R} 2-5

6)引入位移边界条件，修正总体平衡方程

在已形成的总体平衡方程中，由于总体刚度矩阵[K]为一奇异矩阵，即其不存在逆矩阵[K]^-1，因而需考虑所研究区域的位移边界条件(或约束条件),对总体平衡方程组进行修改，以消除总体刚度矩阵的奇异性(从力学意义讲，是消除结构的刚体运动)，这样，才能由总体平衡方程组解出所有未知节点的位移。

7)解方程，求未知节点位移及单元应力

整体平衡方程为一线性迭代方程组，求解此线性代数方程组得到节点的全部位移值{U}。根据求得的位移，利用单元分析结构，可计算出各单元体的应变及应力等。

由于采空区围岩体应力应变关系为非线性的，反映到式(4-5)中，矩阵刚度矩阵[K]是节点位移的非线性函数，因而在结构的总体平衡方程[K]{U}＝{R}中，有限元的总体平衡方程变为

[K(U)]{U}＝{R} 2-6

它步骤在是节点位移的线性方程组，而是非线性方程组。

采空区地表变形ANSYS仿真计算内容

1、地表变形影响因素选择

基于上述地表变形影响因素综合分析，选取采深、采深采厚比、釆宽采深比和上覆岩性为主要影响因素，将这4种主要影响因素和4种水平值进行汇总，如表3所示。生产中实际需要的地表变形评价判据,也就是表中所示某些因素不同水平的组合。

表3釆空区地表变形影响因素和水平

其中上覆岩性的四个水平用如下线性函数进行折减。

G_k＝h·G

C_k＝h·C

E_k＝h·E

υ_k＝υ+(0.5-υ)·(1-h)

式中：h为折减系数，G_k为等效重度，C_k为等效粘聚力，E_k为等效弹性模量，ν_k为等效泊松比，为等效内摩擦角。

上覆岩性I级水平对应材料折减系数取1.0，上覆岩性II级水平对应材料折减系数取0.95，上覆岩性III级水平对应材料折减系数取0.9，上覆岩性IV级水平对应材料折减系数取0.85.

2、地表变形影响因素工况组合

根据表2釆空区地表变形影响因素和水平，将影响因素包括采深、采深采厚比、釆宽采深比和回采率等因素各4个水平组合。研究地表沉降随采深变化规律。保持采深变化，采厚、采宽和上覆岩性不变的工况组合如表4。

表4采空区计算模型正交模型表L₁₆(4⁴)

3、ANSYS数值模型的建立

(1)地质参数取值

在ANSYS程序中，采空区地表变形仿真计算需要输入4个值：天然密度、弹性模量、泊松比、粘聚力、内摩擦角及膨胀角。膨胀角用来控制体积膨胀的大小对压实的颗粒状材料，当材料受剪时，颗粒将会膨胀。如果膨胀角为零，则不会发生体积膨胀；如果膨胀角与内摩擦角相等，则会发生严重的体积膨胀。本发明对膨胀角取值为零。参数的选取如表5。

表5仿真计算岩土力学参数取值

3.1、边界条件的处理和计算范围选择

(1)边界条件的处理：

位移边界为：底部边界的水平与垂直位移固定约束；左右边界的水平位移固定约束(置零)，垂直位移不约束(自由)；上边界为地面线，用自由边界。

(2)计算范围的选择

计算模型范围必须大于拟建构筑物在地基中形成的应力范围和可能造成构筑物附加沉陷变形的老采空区范围。本次计算模型范围，以采空区为中心区域，四周范围取3倍采宽。

3.2、仿真计算模型

利用ANSYS建立的上覆岩土层和采空区的三维地质模型如图2所示。

4、工况1地表变形计算过程举例

依据表2采空区计算工况正交工况表L₁₆(4⁴)，选取工况1的采深为50m，采深采厚比为20，采宽采深比为0.8，上覆岩性为1，在建模中采用了分层建模并使用了包括DP-模型和单元生死等功能。考虑到了仿真的时间成本，计算量和计算机硬件成本等，划分网格则采用了精密度为8的自由网格划分方法。有限元的模型见图3。在有限元计算仿真中，由于划分网格是四面体，它们的顶点并不是规律的分布在地表面上，造成了地表面上的节点数据提取不是完全的成比例。但也表现出一定圆形的对称性。所以主要沿着水平某一方向取一系列离散点作为观测点，在工况1中。观测点包括(497，498，499，500，501)等节点。

表6工况1采空区地表观测点位移

W＝DY_max＝0.85(mm)

W_sp＝DX_max(DZ_max)＝0.32(mm)

其中：X，Y，Z-分别是观测点坐标，(mm)；DX，DY，DZ-分别是观测点沿X、Y、Z方向的位移量，(mm)；i-倾斜变形，(mm/m)；ε-水平变形，(mm/m)；W-地表沉降，(mm)；W_sp-水平位移，(mm)。

在地表所取的观测点中，我们发现节点501拥有较大的竖向变形值，而在节点499处拥有较大的水平位移值。根据公式计算出该模型下的最大沉降值，最大倾斜值，和最大水平变形值如上表6所示。

在采空区工况1中，图3的有限元模型位移云图中显示出地表沉降规律是微围绕着采空区呈现出一定的圆形塌陷，且越往中心沉降越严重。通过查阅相关资料及文献得到，地表沉降多表现为此类形式。

经仿真计算得到各工况的地表沉降数据为：

表7采空区地表沉降数据

工况	1	2	3	4	5	6	7	8
									W(mm)	0.85	0.70	0.88	0.94	4.91	2.31	5.86	3.47
工况	9	10	11	12	13	14	15	16
									W(mm)	15.90	15.92	4.63	6.04	39.29	23.44	10.71	6.17

采用二次函数拟合，得到采空区地表沉降预测模型

二次函数拟合的拟合形式为

其中：x₁,x₂,x₃,x₄分别代表采深，采深采厚比，采宽采深比，上覆岩性。

拟合前先将采深、采深采厚比数据按公式

做归一化处理，数据归化到区间[0.1,0.9]，式中，x_i’-归一化处理后得到的数据；x_max-采深、采深采厚比对应的最大数值；x_min-采深、采深采厚比对应的最小数值；x_i-采深、采深采厚比需要做归一化处理的数据。

拟合系数a₀,a₁,…,a₁₄分别取为(101.69,-79.42,-17.45,-200.52,17.47,84.69,10.36,73.15,-2.05,-33.56,39.54,-0.59,54.08,-17.2,0.96)。

得到的沉降预测模型为：

根据模型预测得出模型中各待定参数的取值：

运用ANSYS有限元软件进行数值模拟仿真，得出各个工况下采空区地表沉降数值与本发明所述多因素采空区地表沉降预测模型计算得出各个工况下采空区地表沉降数值进行对比，数值模拟仿真结果与预测模型计算结果列于下表8。

表8结果误差分析

上述预测模型的平方相关系数R²＝0.997，F值为243.24，p值为0.05，均方误差为0.48，说明预测模型得出的结果是高度精确的。

预测模型计算结果与数值模拟仿真计算结果配对散点图分布如图4：可以看出散点几乎分布在直线y＝x上，说明预测模型计算效果很好，验证成功。

Claims (7)

1.一种采空区地表沉降的预测方法，其特征在于,所述方法包括以下步骤:

步骤一：拟定采空区的采深、采深采厚比、采宽采深比、上覆岩性四因素以及各自四个水平；

步骤二：采用正交实验设计方法对步骤一四因素以及各自四水平按正交表进行设计；

步骤三：运用ANSYS有限元仿真软件对上述工况进行仿真计算，得出各工况的地表沉降数据；

步骤四：采用二次函数拟合，得到采空区地表沉降预测模型；

步骤五：根据模型预测得出模型中各待定参数的取值。

2.根据权利要求1所述的方法，其特征在于：所述步骤一采深的1、2、3、4水平分别对应数值为：50、100、150、200；采深采厚比的1、2、3、4水平分别对应数值为：20、30、40、50；采宽采深比的1、2、3、4水平分别对应数值为：0.8、1、1.2、1.4；上覆岩性的1、2、3、4水平分别对应数值为：1、0.95、0.9、0.85；

步骤二四因素四水平正交实验表及16种计算工况组合如下表1；

表1

工况 采深(m) 采厚(采深采厚比) 采宽(采宽采深比) 上覆岩性 1 50 2.5(20) 40(0.8) 1 2 50 1.67(30) 50(1) 2 3 50 1.25(40) 60(1.2) 3 4 50 1(50) 70(1.4) 4 5 100 5(20) 100(1) 3 6 100 3.33(30) 80(0.8) 4 7 100 2.5(40) 140(1.4) 1 8 100 2(50) 120(1.2) 2 9 150 7.5(20) 180(1.2) 4 10 150 5(30) 210(1.4) 3 11 150 3.75(40) 120(0.8) 2 12 150 3(50) 150(1) 1 13 200 10(20) 280(1.4) 2 14 200 6.67(30) 240(1.2) 1 15 200 5(40) 200(1) 4 16 200 4(50) 160(0.8) 3

3.根据权利要求2所述的方法，其特征在于：所述步骤三仿真计算具体包括以下步骤：

3.1地质参数取值：

在ANSYS程序中，采空区地表变形仿真计算需要输入参数的选取如表2:

表2仿真计算岩土力学参数取值

3.2、边界条件的处理和计算范围选择：

(1)边界条件的处理：

位移边界为：底部边界的水平与垂直位移固定约束；左右边界的水平位移固定约束，垂直位移不约束；上边界为地面线，用自由边界；

(2)计算范围的选择：以采空区为中心区域，四周范围取2-4倍采宽；

3.3、仿真计算模型：利用ANSYS建立不同工况的上覆岩土层和采空区的三维地质模型；

3.4、利用3.3得到的地质模型进行地表变形计算，得到地表沉降数据。

4.根据权利要求3所述的方法，其特征在于：所述步骤三各工况的地表沉降数据为表3：

表3采空区地表沉降数据

工况 1 2 3 4 5 6 7 8 W(mm) 0.85 0.70 0.88 0.94 4.91 2.31 5.86 3.47 工况 9 10 11 12 13 14 15 16 W(mm) 15.90 15.92 4.63 6.04 39.29 23.44 10.71 6.17

5.根据权利要求2所述的方法，其特征在于：所述步骤四二次函数拟合的拟合公式为

其中：x₁,x₂,x₃,x₄分别代表采深，采深采厚比，采宽采深比，上覆岩性。

6.根据权利要求5所述的方法，其特征在于：所述步骤四二次函数拟合前将采深、采深采厚比数据按公式

做归一化处理，数据归化到区间[0.1,0.9]，式中，x_i’-归一化处理后得到的数据；x_max-采深、采深采厚比对应的最大数值；x_min-采深、采深采厚比对应的最小数值；x_i-采深、采深采厚比需要做归一化处理的数据。

7.根据权利要求6所述的方法，其特征在于：所述步骤四沉降预测模型为：