面向地表沉陷动态预计的数值模型岩体力学参数加权反演方法与流程

本发明涉及一种煤矿开采数值模型岩土力学参数反演方法，尤其适用于地表沉陷动态预计的煤矿开采数值模型岩体力学参数加权反演方法。

技术背景

地下煤矿开采引起的地表沉陷在世界各矿区普遍存在。地表沉陷可能会导致地面建筑物损坏。为了保护地面建筑物，必须准确预测地面沉陷，使建筑物的损害程度小于其抗变形能力，从而合理设计地下开采。

考虑到采矿活动是一个时空过程，地表动态移动与变形对地面建筑物的影响可能更为显著。由于采动影响下地表动态移动与变形和工作面地质采矿条件、开采方式、工作面推进速度及开采时间有关，准确预测地表动态移动与变形是一件极其困难的事。部分学者围绕地表动态移动与变形预测开展了探索性研究，提出了多种时间影响函数模型。knothe提出了用于地表变形动态预计的knothe函数，这也是目前应用最广泛的地表动态移动与变形预测方法。俄罗斯a.a.baryakh等建立了考虑时间因素的动态地表沉陷预计方法；刘玉成利用数学分析方法改进了knothe时间函数，可求出沉陷速度最大的时间；徐洪钟针对knothe时间函数描述的不足，提出了基于logistic增长模型的时间函数；但目前大多数地表沉陷动态预测模型处在理论研究阶段，没有较多或较好的工程实例支撑。

随着计算机技术和现代力学理论的不断发展，数值模拟方法因能较为全面考虑矿区地质情况，越来越广泛的被用于地表沉陷预计。例如，xu利用a3-dfinitedifferencemethod预测wutongmine开采地表沉陷；vyazmensky利用等效连续体和混合方法实现了对地表沉陷的预测；guo采用数值模拟方法研究了wongawillistrippillar的地表沉陷规律；ma采用数值模拟方法研究了不同充填率下充填开采地表沉陷规律。但上述研究主要是关于开采稳定后的地表沉陷进行预测，而对于开采过程中地表动态移动与变形预测则鲜有研究。主要原因是由于选取的岩体力学参数难以满足地表沉陷动态预计的需求。因此，目前缺少面向地表沉陷动态预计的数值模型岩体力学参数反演方法。

技术实现要素：

针对上述技术问题，提供了一种解决了现有数值模型用于地表沉陷动态预计精度较差的难题。对于数值模型的建立与校准、地表沉陷动态预计、地面建筑物保护等均具有重要的实际意义和应用价值的面向地表沉陷动态预计的数值模型岩体力学参数加权反演方法。

为实现上述技术目的，本发明的面向地表沉陷动态预计的数值模型岩体力学参数加权反演方法，其特征在于步骤如下：

步骤1：建立研究区域初始数值模型，包括区域地质采矿条件、工作面分布情况、地应力监测结果及室内岩体力学参数，利用flac3d建立要研究区域的初始模型；

步骤2：利用遗传算法和初始数值模型，反演得到适用于非充分采动情况下地表沉陷预计的静态岩体力学参数；

步骤3：利用不同采动程度下的地表移动与变形实测数据，并根据不同岩体力学参数设置反演过程中各力学参数的相应权重，从而得到面向地表沉陷动态预计的数值模型岩体力学参数反演方法；

步骤4：结合不同采动程度下岩体力学参数的权重，利用加权反演方法，获取面向地表沉陷动态预计的数值模型岩体力学参数，根据得到的不同采动程度下岩体力学参数值及相应的权重，利用上述计算公式可以计算得到面向地表沉陷动态预计的全局数值模型岩体力学参数。

建立研究区域初始数值模型具体步骤为：

面向要进行地表沉陷动态预计的实验区域，收集研究区域地质采矿条件、工作面分布情况、地应力监测结果及室内岩体力学参数实验结果，根据获取的研究区域所有资料，利用flac3d建立要研究区域的初始模型，同时为了减小模型边界效应的影响，在模型边界预留1.4h的宽度，h为采深。

反演非充分采动和充分采动情况下数值模型静态岩体力学参数：

开采沉陷领域通常用下沉、移动、倾斜、水平变形及曲率来评价地表移动与变形，参数反演时候可根据实际情况将地表下沉、移动、倾斜、水平变形及曲率作为评价指标，根据研究区域非充分采动情况下实测的地表下沉、水平移动、地表倾斜、水平变形及曲率等评价指标，利用遗传算法和上述flac3d建立的初始数值模型，反演得到适用于非充分采动情况下地表沉陷预计的静态岩体力学参数，岩体力学参数包括弹性模量、泊松比、抗拉强度、内聚力和内摩擦角，其中弹性模量和泊松比称为岩体力学刚度参数，抗拉强度、内聚力和内摩擦角称为岩体力学强度参数，设上述得到弹性模量e1、泊松比p1、抗拉强度s1、内聚力c1和内摩擦角f1，同时根据研究区域充分采动情况下实测的地表下沉、水平移动、地表倾斜、水平变形及曲率的评价指标，利用遗传算法和上述flacd建立的初始数值模型，反演得到适用于充分采动情况下地表沉陷动态预计的静态岩体力学参数，此时得到弹性模量e2、泊松比p2、抗拉强度s2、内聚力c2和内摩擦角f2。

确定不同采动程度下岩体力学参数对评价指标的权重具体步骤为：

开采沉陷是一个典型的时空过程，不同采动程度下模型岩体力学参数对地表下沉、水平移动、地表倾斜、水平变形及曲率评价指标的敏感性不同，可以充分利用不同采动程度下的地表移动与变形实测数据，并根据不同岩体力学参数设置反演过程中各力学参数的相应权重，从而得到面向地表沉陷动态预计的数值模型岩体力学参数反演方法，

为得到上述目的，基于上述flac3d建立的数值模型，就岩体力学参数对地表下沉、水平移动、倾斜、水平变形和曲率的影响进行单因素分析。所谓单因素分析就是改变弹性模量、泊松比、抗拉强度、内聚力、内摩擦角的某一参数时，其他参数取为flac3d建立的初始模型中的相应值。考虑到岩体力学参数之间的量纲差异较大，将flac3d建立的初始模型评价权重设为1，依据差分法定义相应权重，利用公式：式中，为指标ui对力学参数xj的权重；dui(x)为指标变化值；dxj为参数变化值；

根据上述建立的公式和单因素试验结果，可以得到不同采动程度下岩体力学参数对评价指标的权重，包括充分采动和非充分采动情况下岩体力学参数对评价指标的权重，得到非充分采动情况下弹性模量的权重为pe1、泊松比的权重为pp1、抗拉强度的权重为ps1、内聚力的权重为pc1、内摩擦角的权重为pf1；充分采动情况下弹性模量的权重为pe2、泊松比的权重为pp2、抗拉强度的权重为ps2、内聚力的权重为pc2、内摩擦角的权重为pf2。

获取面向地表沉陷动态预计的数值模型岩体力学参数具体步骤为：

反演得到了非充分采动情况下弹性模量e1、泊松比p1、抗拉强度s1、内聚力c1、内摩擦角f1以及充分采动情况下弹性模量e2、泊松比p2、抗拉强度s2、内聚力c2、内摩擦角f2，上述步骤3得到了非充分采动情况下弹性模量权重pe1、泊松比权重pp1、抗拉强度权重ps1、内聚力权重pc1、内摩擦角权重pf1以及充分采动情况下弹性模量权重pe2、泊松比权重pp2、抗拉强度权重ps2、内聚力权重pc2、内摩擦角pf2，结合不同采动程度下岩体力学参数的权重，利用加权反演方法，就可以获取面向地表沉陷动态预计的数值模型岩体力学参数，计算公式为：

x＝x1×px1+x2×px2

根据得到的不同采动程度下岩体力学参数值及相应的权重，利用上述计算公式可以计算得到面向地表沉陷动态预计的全局数值模型岩体力学参数，式中x1为非充分采动参数取值，px1为非充分采动参数权重，x2为充分采动参数取值，px2为充分采动参数权重。

有益效果：本发明兼顾了开采沉陷的时空过程，创造性的利用不同采动程度下岩体力学参数对评价指标的动态敏感性，基于加权的思想，提出了面向地表沉陷动态预计的数值模型岩体力学参数加权反演方法，解决了现在岩体力学参数反演方法应用于地表沉陷动态预计难以满足工程需求的难题，改变了数值模型只能应用于地表沉陷静态分析与预计的现状，提高了数值模型应用于地表沉陷动态预计的精度。本方法对于数值模型的建立与校准、地表沉陷动态预计、地面建筑物保护等均具有重要的实际意义和应用价值。

附图说明

图1为本发明实施的一种面向地表沉陷动态预计的数值模型岩体力学参数加权反演方法流程图；

图2为本发明实施例得到的岩体力学参数权重与塑性区占开采影响区百分比的关系图；

图3为本发明实施例得到的不同岩体力学参数反演方法用于地表沉陷动态预计结果对比图。

具体实施方式

下面将结合图和具体实施过程对本发明做进一步详细说明：

如图1所示，本发明的面向地表沉陷动态预计的数值模型岩体力学参数加权反演方法，包括如下步骤：

步骤1：建立研究区域初始数值模型，包括区域地质采矿条件、工作面分布情况、地应力监测结果及室内岩体力学参数，利用flac3d建立要研究区域的初始模型；

步骤2：利用遗传算法和初始数值模型，反演得到适用于非充分采动情况下地表沉陷预计的静态岩体力学参数；

步骤3：利用不同采动程度下的地表移动与变形实测数据，并根据不同岩体力学参数设置反演过程中各力学参数的相应权重，从而得到面向地表沉陷动态预计的数值模型岩体力学参数反演方法；

步骤4：结合不同采动程度下岩体力学参数的权重，利用加权反演方法，获取面向地表沉陷动态预计的数值模型岩体力学参数，根据得到的不同采动程度下岩体力学参数值及相应的权重，利用上述计算公式可以计算得到面向地表沉陷动态预计的全局数值模型岩体力学参数。

具体的：

建立研究区域初始数值模型具体步骤为：

面向要进行地表沉陷动态预计的实验区域，收集研究区域地质采矿条件、工作面分布情况、地应力监测结果及室内岩体力学参数实验结果，根据获取的研究区域所有资料，利用flac3d建立要研究区域的初始模型，同时为了减小模型边界效应的影响，在模型边界预留1.4h的宽度，h为采深。

反演非充分采动和充分采动情况下数值模型静态岩体力学参数：

开采沉陷领域通常用下沉、移动、倾斜、水平变形及曲率来评价地表移动与变形，参数反演时候可根据实际情况将地表下沉、移动、倾斜、水平变形及曲率作为评价指标，根据研究区域非充分采动情况下实测的地表下沉、水平移动、地表倾斜、水平变形及曲率等评价指标，利用遗传算法和上述flac3d建立的初始数值模型，反演得到适用于非充分采动情况下地表沉陷预计的静态岩体力学参数，岩体力学参数包括弹性模量、泊松比、抗拉强度、内聚力和内摩擦角，其中弹性模量和泊松比称为岩体力学刚度参数，抗拉强度、内聚力和内摩擦角称为岩体力学强度参数，设上述得到弹性模量e1、泊松比p1、抗拉强度s1、内聚力c1和内摩擦角f1，同时根据研究区域充分采动情况下实测的地表下沉、水平移动、地表倾斜、水平变形及曲率的评价指标，利用遗传算法和上述flacd建立的初始数值模型，反演得到适用于充分采动情况下地表沉陷动态预计的静态岩体力学参数，此时得到弹性模量e2、泊松比p2、抗拉强度s2、内聚力c2和内摩擦角f2。

确定不同采动程度下岩体力学参数对评价指标的权重具体步骤为：

开采沉陷是一个典型的时空过程，不同采动程度下模型岩体力学参数对地表下沉、水平移动、地表倾斜、水平变形及曲率评价指标的敏感性不同，可以充分利用不同采动程度下的地表移动与变形实测数据，并根据不同岩体力学参数设置反演过程中各力学参数的相应权重，从而得到面向地表沉陷动态预计的数值模型岩体力学参数反演方法，

为得到上述目的，基于上述flac3d建立的数值模型，就岩体力学参数对地表下沉、水平移动、倾斜、水平变形和曲率的影响进行单因素分析。所谓单因素分析就是改变弹性模量、泊松比、抗拉强度、内聚力、内摩擦角的某一参数时，其他参数取为flac3d建立的初始模型中的相应值。考虑到岩体力学参数之间的量纲差异较大，将flac3d建立的初始模型评价权重设为1，依据差分法定义相应权重，利用公式：式中，为指标ui对力学参数xj的权重；dui(x)为指标变化值；dxj为参数变化值；

根据上述建立的公式和单因素试验结果，可以得到不同采动程度下岩体力学参数对评价指标的权重，包括充分采动和非充分采动情况下岩体力学参数对评价指标的权重，得到非充分采动情况下弹性模量的权重为pe1、泊松比的权重为pp1、抗拉强度的权重为ps1、内聚力的权重为pc1、内摩擦角的权重为pf1；充分采动情况下弹性模量的权重为pe2、泊松比的权重为pp2、抗拉强度的权重为ps2、内聚力的权重为pc2、内摩擦角的权重为pf2。

获取面向地表沉陷动态预计的数值模型岩体力学参数具体步骤为：

反演得到了非充分采动情况下弹性模量e1、泊松比p1、抗拉强度s1、内聚力c1、内摩擦角f1以及充分采动情况下弹性模量e2、泊松比p2、抗拉强度s2、内聚力c2、内摩擦角f2，上述步骤3得到了非充分采动情况下弹性模量权重pe1、泊松比权重pp1、抗拉强度权重ps1、内聚力权重pc1、内摩擦角权重pf1以及充分采动情况下弹性模量权重pe2、泊松比权重pp2、抗拉强度权重ps2、内聚力权重pc2、内摩擦角pf2，结合不同采动程度下岩体力学参数的权重，利用加权反演方法，就可以获取面向地表沉陷动态预计的数值模型岩体力学参数，所述参数共有5个，每个参数均包含非充分采动和充分采动的取值，计算公式为：

x＝x1×px1+x2×px2

根据得到的不同采动程度下岩体力学参数值及相应的权重，利用上述计算公式可以计算得到面向地表沉陷动态预计的全局数值模型岩体力学参数，式中x1为非充分采动参数取值，px1为非充分采动参数权重，x2为充分采动参数取值，px2为充分采动参数权重。

具体实施例：

步骤1：建立研究区域初始数值模型

选择兖州矿区某工作面为例进行本发明方法的详细说明。该工作面煤层为近水平煤层，采厚是5m。岩层从下到上依次为粉砂岩、煤、细砂岩、粉砂岩、砾岩、细砂岩、粉砂岩、细砂岩、中砂岩及表土层。根据该工作面的地质采矿条件、地应力监测结果、工作面分布及室内岩体力学参数实验结果，利用flac3d软件建立了一个长1200m，宽1200m，高300m的一个数值模型。模型网格平面尺寸为30m×30m，高度随岩层高度略有不同，共划分12600单元，23064个节点。模型底部边界取u＝v＝0全约束边界，顶部为自由边界，左右边界为水平位移约束边界。

步骤2：反演非充分采动和充分采动情况下数值模型静态岩体力学参数

选择地表下沉和地表水平移动为评价指标。根据该工作面第6期非充分采动情况下的实测地表下沉和地表水平移动值，利用遗传算法和初始模型，反演得到该工作面非充分采动情况下数值模型弹性模量、泊松比、内聚力、抗拉强度及内摩擦角等岩体力学静态参数；同时根据该工作面第13期充分采动情况下的实测地表下沉和地表水平移动值，利用遗传算法和初始模型，反演得到该工作面充分采动情况下数值模型弹性模量、泊松比、内聚力、抗拉强度及内摩擦角等岩体力学静态参数。

步骤3：确定不同采动程度下岩体力学参数对评价指标的权重

根据上述flac3d建立的初始模型，利用单因素实验方法研究不同采动程度下岩体力学参数对地表下沉和水平移动的权重。考虑到岩体力学参数之间的量纲差异较大，将flac3d建立的初始模型评价权重设为1，依据差分法定义相应权重，计算方式如下：

式中，为指标ui对力学参数xj的权重；dui(x)为指标变化值；dxj为参数变化值

如图2所示，根据上述建立的公式和单因素试验结果，可以得到不同采动程度下岩体力学参数对评价指标的权重，包括充分采动和非充分采动情况下岩体力学参数对评价指标的权重。根据图2可以得到非充分采动和充分采动情况下岩体力学刚度参数和强度参数的权重。岩体力学刚度参数包括弹性模型和泊松比，岩体力学强度参数包括抗拉强度、内聚力和内摩擦角。

步骤4：获取面向地表沉陷动态预计的数值模型岩体力学参数

步骤2中已通过遗传算法反演得到了适用于非充分采动和充分采动情况下地表沉陷预计的静态岩体力学参数，再结合步骤3得到的实验区不同采动程度下岩体力学参数对评价指标的权重，利用加权的思想，就可以得到面向地表下沉动态预计的数值模型岩体力学参数。

本发明方法应用于了兖州矿区某工作面。图3为本发明实施例得到的不同岩体力学参数反演方法用于地表沉陷动态预计结果对比图。图中实测值为该工作面某点实测的地表下沉值；加权反演是指利用本发明中的方法得到的岩体力学参数进行的地表沉陷动态预计结果；常规反演是指基于第13期实测数据利用遗传算法得到的岩体力学参数进行的地表沉陷动态预计结果。6期时，常规反演和实测值相差120mm，加权反演和实测值仅相差17mm；9期时，常规反演和实测值相差184mm，加权反演和实测值相差52mm；11期时，常规反演和实测值相差89mm，加权反演和实测值相差40mm；13期时，加权反演和实测值相差65mm；加权反演和实测值相差31mm；14期时，常规反演和实测值相差265mm，加权反演和实测值相差58mm。因此，本发明能显著提高地表沉陷动态预计精度，兼顾了煤矿开采的时空特点，克服了常规反演得到的数值模型岩体力学参数难以满足地表沉陷动态预计的需求。