CN103726481A - 一种重复开采地表变形预测方法 - Google Patents

Abstract

本发明一种重复开采地表变形预测方法，属于“三下”开采和数字矿山技术领域，本发明通过对地表倾斜值的叠加、曲率值的叠加、水平移动值的叠加和水平变形值的叠加，有效实现多工作面开采后的地表变形的预测；此外，本发明还可对条带开采地表变形、对同一层位多个工作面开采的地表变形和对不同层位开采后有影响重叠区域变形进行预测。

Description

一种重复开采地表变形预测方法

技术领域

本发明属于“三下”开采和数字矿山技术领域，具体涉及一种重复开采地表变形预测方法。

背景技术

随着煤炭资源的日渐贫乏，“三下”压煤问题越来越突出，煤炭资源开采出必然对地表产生影响，准确掌握开采后影响程度是提前做好地表保护和方案设计的一个重要依据，现有的地表开采变形预测中，得到现场验证和普遍认可的算法只有单工作面开采后的地表变形预测算法，不能对多工作面开采后的地表变形进行预测。

单工作面开采的地表移动变形过程相对重复开采地表变形过程简单很多，开采过程中不涉及到抵消和加剧过程，单工作开采的算法没有考虑重复采动对地表变形的影响过程，重复开采地表的破坏变形叠加除下沉值外其它的变形值不能应用算数和的办法叠加，在开采过程中倾斜值、曲率值、水平变形值和水平移动值都是在开采过程中存在抵消的过程，这样就需要一个合理可靠的算法来实现倾斜值的叠加、曲率值的叠加、水平移动值的叠加和水平变形值的叠加。

开发可以应用于重复开采后的地表变形预测的算法和软件势在必行，因为多工作面开采后地表变形移动过程相对于单工作面开采后的地表移动变形过程复杂得多，而且人工计算量非常大，要实现多工作面开采后的地表变形预测必须开发出一个行之有效的算法和计算机系统。

发明内容

针对现有技术存在的不足，本发明一种重复开采地表变形预测方法，以达到实现对重复开采地表变形值、条带开采地表变形值、同一层位多个工作面开采的地表变形值和开采后有影响重叠区域进行变形值进行预测的目的。

一种重复开采地表变形预测方法，包括以下步骤：

步骤1、确定预测地表的开采工作面的顺序，并在每个工作面中设定采样点的间距，确定每个工作面的参数，包括下沉系数、水平移动系数、拐点偏距、主要影响正切值、开采影响传播角、最大下沉角、煤层倾角、工作面尺寸、采深和采厚10个参数，并根据上述10个参数获得每个采样点的走向各变形值和倾斜方向各变形值，所述的变形值包括下沉值、倾斜值、曲率值、水平移动值和水平变形值五个值；

步骤2、根据每个采样点的走向各变形值和倾斜方向各变形值，采用概率积分法获得每个工作面的沿工作面走向各变形值和沿工作面倾斜方向各变形值，具体如下：

沿工作面走向各变形值的预测公式如下：

$\left\{\begin{array}{c}{W}^0\left(x\right)=C_{\mathrm{ym}}[W(x;t_1)-W(x-l;t_2)]\\ i^0\left(x\right)=C_{\mathrm{ym}}[i(x;t_1)-i(x-l;t_2)]\\ K^0\left(x\right)=C_{\mathrm{ym}}[K(x;t_1)-K(x-l;t_2)]\\ U^0\left(x\right)=C_{\mathrm{ym}}[U(x;t_1)-U(x-l;t_2)]\\ {\mathrm{\ϵ}}^0\left(x\right)=C_{\mathrm{ym}}[\mathrm{\ϵ}(x;t_1)-\mathrm{\ϵ}(x-l;t_2)]\\ C_{\mathrm{ym}}=\frac{W_{\mathrm{my}}^0}{W^0}\end{array}\right.---\left(1\right)$

其中，l为工作面走向长度；x为走向坐标，W(x;t₁)表示位于工作面走向方向距离选定原点的x处的下沉值；W(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处下沉值；i(x;t₁)表示位于工作面走向方向距离选定原点的x处倾斜值；i(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处倾斜值；K(x;t₁)表示位于工作面走向方向距离选定原点的x处曲率值；K(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处曲率值；U(x;t₁)表示位于工作面走向方向距离选定原点的x处水平移动值；U(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处水平移动值；ε(x;t₁)表示位于工作面走向方向距离选定原点的x处水平变形值；ε(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处水平变形值；t₁表示下山方向开采边界的10个参数、t₂表示上山方向开采边界的10个参数；W⁰(x)表示走向下沉值，i⁰(x)表示走向倾斜值，K⁰(x)表示走向曲率值，U⁰(x)表示走向水平移动值，ε⁰(x)表示走向水平变形值：

表示走向充分采动倾向主断面的最大下沉值；W⁰表示走向和倾向均为充分采动时的地表最大下沉值，C_ym表示倾向采动程度系数；

沿工作面倾斜方向各变形值的预测公式为：

$\left\{\begin{array}{c}{W}^0\left(y\right)=C_{\mathrm{xm}}[W(y;t_3)-W(y-L;t_4)]\\ i^0\left(y\right)=C_{\mathrm{xm}}[i(y;t_3)-i(y-L;t_4)]\\ K^0\left(y\right)=C_{\mathrm{xm}}[K(y;t_3)-K(y-L;t_4)]\\ U^0\left(y\right)=C_{\mathrm{xm}}[U(y;t_3)-U(y-L;t_4)]\\ {\mathrm{\ϵ}}^0\left(y\right)=C_{\mathrm{xm}}[\mathrm{\ϵ}(y;t_3)-\mathrm{\ϵ}(y-L;t_4)]\\ C_{\mathrm{xm}}=\frac{W_{\mathrm{mx}}^0}{W^0}\end{array}\right.---\left(2\right)$

其中，L表示工作面倾斜长度；

表示倾向充分采动时走向主断面的最大下沉值；C_xm表示走向采动程度系数；y表示倾斜方向坐标；W⁰(y)表示倾斜方向下沉值，i⁰(y)表示倾斜方向倾斜值，K⁰(y)表示倾斜方向曲率值，U⁰(y)表示倾斜方向水平移动值，ε⁰(y)表示倾斜方向水平变形值；t₃表示走向左侧开采边界的10个参数；t₄表示走向右侧开采边界的10个参数；W(y;t₃)表示位于工作面倾斜方向距离选定原点的y处的下沉值；W(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处的下沉值；i(y;t₃)表示位于工作面倾斜方向距离选定原点的y处倾斜值；i(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处倾斜值；K(y;t₃)表示位于工作面倾斜方向距离选定原点的y处曲率值；K(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处曲率值；U(y;t₃)表示位于工作面倾斜方向距离选定原点的y处水平移动值；U(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处水平移动值；ε(y;t₃)表示位于工作面倾斜方向距离选定原点的y处水平变形值；ε(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处水平变形值；

步骤3、采用距离幂次反比法将每个工作面的沿工作面走向各变形值进行叠加，将每个工作面的沿工作面倾斜方向各变形值进行叠加，具体如下：

步骤3-1、按照地表开采工作面的顺序，确定第一次开采的工作面和第二次开采的工作面之间的重叠区域，将重叠区域中第一次开采的工作面上某一采样点作为圆心，以第一次开采的工作面中采样点的间距作为半径，确定搜索范围；

步骤3-2、在所确定的搜索范围内对应搜索第二次开采的工作面中采样点的数据，包括沿工作面走向各变形值和沿工作面倾斜方向各变形值，并采用距离幂次反比法对搜索到的采样点数据进行叠加计算，公式如下：

$Z_d=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{Z_i}{d_i}}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{d_i}}---\left(3\right)$

其中，Z_d表示搜索到的点经过距离幂次反比法得到的修正变形值；Z_i表示搜索到的采样点变形值，包括下沉值、倾斜值、曲率值、水平移动值和水平变形值；d_i表示搜索到第二次开采工作面的采样点距离搜索范围圆心的距离；n表示搜索到的点的个数；

确定叠加后的采样点变形值：

Z=Z_x+Z_d    (4)

其中，Z表示叠加后的搜索范围圆心的变形值；Z_x表示搜索范围圆心点的原始变形值；

步骤3-3、依次将重叠区域中第一次开采的工作面上其他采样点作为圆心，以第一次开采的工作面中采样点的间距作为半径，确定新的搜索范围，并重复执行步骤3-2，直至重叠区域中第一次开采工作面的所有采样点均搜索完成；

步骤3-4、将经过叠加合成后的第一次开采的工作面和第二次开采的工作面作为一个整体，重复步骤3-1至步骤3-3与下一次开采的工作面进行叠加合成，直至所有工作面采样点的沿工作面走向各变形值和沿工作面倾斜方向各变形值均完成叠加合成；

步骤4、将每个采样点的变形值进行合成，获得采样点实际预测变形值，具体为：

步骤4-1、对采样点的沿工作面走向下沉值和沿工作面倾斜方向下沉值进行合成，获得实际预测下沉值，公式如下：

$W(x,y)=\frac{1}{W^0}{W}^0\left(x\right)W^0\left(y\right)---\left(5\right)$

其中，W(x,y)表示坐标为(x,y)处的采样点实际下沉值；

步骤4-2、对采样点的沿工作面走向倾斜值和沿工作面倾斜方向倾斜值进行合成，获得实际预测倾斜值，公式如下：

$i(x,y,\mathrm{\α})=\frac{1}{W^0}[i^0\left(x\right)W^0\left(y\right)\cos \mathrm{\α}+W^0\left(x\right)i^0\left(y\right)\sin \mathrm{\α}]---\left(6\right)$

其中，i(x,y,α)表示点(x,y)沿α方向上的实际倾斜值；

根据实际变形值的条件，获得如下公式；

$\frac{\∂i(x,y,\mathrm{\α})}{\∂\mathrm{\α}}{|}_{\mathrm{\α}={\mathrm{\α}}_i}=0---\left(7\right)$

此时获得

$\mathrm{\α}={\mathrm{\α}}_i=\arctan \frac{W^0\left(x\right)i^0\left(y\right)}{i^0\left(x\right)W^0\left(y\right)}---\left(8\right)$

其中，α_i表示倾斜变形值实际值方向；

将获得的方向α_i代入公式(6)中，即获得实际倾斜值i(x,y,α_i)；

步骤4-3、对采样点的沿工作面走向曲率值和沿工作面倾斜方向曲率值进行合成，获得实际预测曲率值，公式如下：

$K(x,y,\mathrm{\α})=\frac{1}{W^0}[K^0\left(x\right)W^0\left(y\right){\cos }^2\mathrm{\α}+K^0\left(y\right)W^0\left(x\right){\sin }^2\mathrm{\α}+i^0\left(x\right)i^0\left(y\right)\sin 2\mathrm{\α}]---\left(9\right)$

其中，K(x,y,α)表示点(x,y)沿α方向上的实际曲率值；

根据实际变形值的条件，获得如下公式；

$\frac{\∂K(x,y,\mathrm{\α})}{\∂k}{|}_{\mathrm{\α}={\mathrm{\α}}_k}=0---\left(10\right)$

此时获得

$\mathrm{\α}={\mathrm{\α}}_K=\frac{1}{2}\arctan \frac{2i^0\left(x\right)i^0\left(y\right)}{K^0\left(x\right)W^0\left(y\right)-K^0\left(y\right)W^0\left(x\right)}---\left(11\right)$

其中，α_K表示曲率值实际值方向；

将α_K值代入式(9)获得实际曲率值K(x,y,α_K)；

步骤4-4、对采样点的沿工作面走向水平位移值和沿工作面倾斜方向水平位移值进行合成，获得实际预测水平位移值，公式如下：

$U(x,y.\mathrm{\α})=\frac{1}{W^0}[U^0\left(x\right)W^0\left(y\right)\cos \mathrm{\α}+U^0\left(y\right)W^0\left(x\right)\sin \mathrm{\α}]---\left(12\right)$

其中，U(x,y,α)表示点(x,y)沿α方向上的实际水平位移值；此时，α=α_i；

步骤4-5、对采样点的沿工作面走向水平变形值和沿工作面倾斜方向水平变形值进行合成，获得实际预测水平变形值，公式如下：

$\begin{array}{c}\mathrm{\ϵ}(x,y,\mathrm{\α})=\frac{1}{W^0}\left\{\right[{\mathrm{\ϵ}}^0\left(x\right)W^0\left(y\right){\cos }^2\mathrm{\α}+{\mathrm{\ϵ}}^0\left(y\right)W^0\left(x\right){\sin }^2\mathrm{\α}]\\ +[U^0\left(x\right)i^0\left(y\right)+U^0\left(y\right)i^0\left(x\right)]\sin \mathrm{\α}\cos \mathrm{\α}\}\end{array}---\left(13\right)$

其中，ε(x,y,α)表示点(x,y)沿α方向上的实际水平变形值；此时，α=α_K；

步骤5、判断确定获得的采样点实际预测变形值所属范围，采取措施，具体为：

当水平变形值在0～2范围内，或曲率值在0～0.2范围内，或倾斜值在0～3范围内时，则建筑物为轻微损坏，建筑物可以不修或者简单维修；所述的简单维修；

当水平变形值在2～4范围内，或曲率值在0.2～0.4范围内，或倾斜值在3～际范围内时，则建筑物为轻度损坏，需要加强建筑物的监测，并及时进行简单维修和加固；

当水平变形值在4～6范围内，或曲率值在0.4～0.6范围内，或倾斜值在6～10范围内时，则建筑物为中度损坏，需要实时检测、及时维修和加强防护；或者采取降低采高或采取充填开采措施，减少地表的沉陷，控制地表变形在保护标准以内；

当水平变形值大于6，或曲率值大于0.6，或倾斜值大于10时，则建筑物为严重损坏，需要人员撤离，并进行大规模维修或拆建。

步骤1所述的根据10个参数获得每个采样点的走向各变形值和倾斜方向各变形值，具体如下：

(1)计算每个采样点的走向各变形值如下：

计算下沉值公式：

$W\left(x\right)=\frac{W^0}{2}[\mathrm{erf}\left(\frac{\sqrt{\mathrm{\π}}}{r}x\right)+1]---\left(14\right)$

其中，W(x)表示x处的下沉值；W⁰表示走向和倾向均为充分采动时的地表最大下沉值，W⁰=mqcosα，其中，m表示煤层开采厚度；δ表示煤层倾角；q表示下沉系数；

表示概率积分函数，其中，r表示主要影响半径，

字表示采深，tanβ表示主要影响正切值；

计算倾斜值公式：

$i\left(x\right)=\frac{W^0}{r}{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(15\right)$

其中，i(x)表示x处的倾斜值；

计算曲率值公式：

$K\left(x\right)=-\frac{2\mathrm{\π}{W}^0}{r^3}x{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(16\right)$

其中，K(x)表示x处的曲率值；

计算水平移动值公式：

$U\left(x\right)=\mathrm{bri}\left(x\right)=b{W}^0{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(17\right)$

其中，U(x)表示x处的水平移动值；b表示水平移动系数；

计算水平变形值公式：

$\mathrm{\ϵ}\left(x\right)=\mathrm{brK}\left(x\right)=-\frac{2\mathrm{\πb}{W}^0}{r^2}x{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(18\right)$

其中，ε(x)表示x处的水平变形值；

Ⅲ2Ⅳ计算每个采样点的倾斜方向各变形值如下：

计算下沉值公式：

$W\left(y\right)=\frac{W^0}{2}[\mathrm{erf}\left(\frac{\sqrt{p}}{r}y\right)+1]---\left(19\right)$

其中，W(y)表示y处的下沉值；

计算倾斜值公式：

$i\left(y\right)=\frac{W^0}{r}{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(20\right)$

其中，i(y)表示y处的倾斜值；

计算曲率值公式：

$K\left(y\right)=-\frac{2\mathrm{\π}{W}^0}{r^3}y{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(21\right)$

其中，K(y)表示y处的曲率值；

计算水平移动值公式：

$U\left(y\right)=\mathrm{bri}\left(y\right)=b{W}^0{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(22\right)$

其中，U(y)表示y处的水平移动值；

计算水平变形值公式：

$\mathrm{\ϵ}\left(y\right)=\mathrm{brK}\left(y\right)=-\frac{2\mathrm{\πb}{W}^0}{r^2}x{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(23\right)$

其中，ε(y)表示y处的水平变形值。

本发明优点：

本发明一种重复开采地表变形预测方法，具有以下方面的优点：

1、可以实现重复开采地表五个变形值的预测；

2、可以对条带开采地表变形进行预测；

3、可以对同一层位多个工作面开采的地表变形进行预测；

4、可以对不同层位，开采后有影响重叠区域进行变形预测。

附图说明

图1为本发明一种实施例的重复采动地表变形过程图；

图2为本发明一种实施例的重复开采地表变形预测方法流程图；

图3为本发明一种实施例的最大下沉角示意图，其中，图(a)为非充分采动下沉角示意图，图(b)为充分采动下沉角示意图

图4为本发明一种实施例的开采影响传播角示意图；

图5为本发明一种实施例的主要影响正切值示意图；

图6为本发明一种实施例的拐点偏距示意图；

图7为本发明一种实施例的走向倾向倾斜值示意图，其中，图(a)为沿走向方向的倾斜值等值线示意图，图(b)为沿倾斜方向的倾斜值等值线示意图；

图8为本发明一种实施例的走向倾向水平移动值示意图，其中，图(a)为沿走向方向的水平移动值等值线示意图，图(b)为沿倾斜方向的水平移动值等值线示意图；

图9为本发明一种实施例的走向倾向曲率值示意图，其中，图(a)为沿走向方向的曲率值等值线示意图，图(b)为沿倾斜方向的曲率值等值线示意图；

图10为本发明一种实施例的走向倾向水平变形值示意图，其中，图(a)为沿走向方向的水平变形值等值线示意图，图(b)为沿倾斜方向的水平变形值等值线示意图；

图11为本发明一种实施例的工作面位置示意图；

图12为本发明一种实施例的下沉值比较图，其中，图(a)为分采下沉值示意图，图(b)为单采下沉值示意图；

图13为本发明一种实施例的倾斜值比较图，其中，图(a)为分采倾斜值示意图，图(b)为单采倾斜值示意图；

图14为本发明一种实施例的曲率值对比图，其中，图(a)为分采曲率值示意图，图(b)为全采曲率值示意图；

图15为本发明一种实施例的水平变形值对比图，其中，图为(a)分采水平变形值示意图，图(b)为全采水平变形值示意图；

图16为本发明一种实施例的水平移动值对比图，其中，图(a)为分采水平移动值示意图，(b)为单采水平移动值示意图。

具体实施方式

下面结合附图对本发明一种实施例做进一步说明。

多工作面开采的地表移动变形过程与单工作面开采地表移动变形过程不同：

(1)单工作面开采地表移动变形过程

地表移动盆地是在工作面推进过程中逐渐形成的。当回采工作面自开切眼开始向前推进的过程中随着工作面推进尺度的增加和时间的增长，上覆岩层逐渐下沉，开采影响逐渐波及到地表，引起地表下沉，继而形成倾斜变形、曲率变形、水平移动和水平变形等。随着工作面继续向前推进和时间的推移，采空区面积不断的增大，地表的影响范围也在逐步的扩大，地表形成的盆地也逐渐扩大，而地表的下沉首先是逐渐在变大，如果工作面达到充分采动，以后下沉值会保持在一个定值不会再变大，随着时间的推移在地表形成一个比采空区面积大得多的下沉盆地。

(2)多工作面开采地表移动变形形成过程

多工作面的依次开采地表沉陷盆地的形成过程如图1所示，图1展示了地表移动盆地随开采的工作面个数逐渐的增加而形成地表沉陷盆地的过程。当第一工作面开采结束时地表沉陷盆地如图1中W₁盆地所示，在地表形成一个小盆地W₁。当开采第二个工作面时沉陷范围扩大为W₂，当开采第三个工作面时沉陷范围扩大为W₃，当开采第四个工作面时沉陷范围扩大为W₄，依次类推，如果中间过程中达到充分采动则下沉的最大值改变，而下沉范围在扩大。当工作开采到第四个工作面不在向前开采时，地表仍继续下沉，大约经过3～5年的时间地表沉陷逐渐稳定，如图1中的W'₀₄所示。并且在地表移动向前发展的过程中，其后方的地表点仍在继续沉陷，但其沉陷的速度逐渐在减慢，一直到稳定。通过日常的实践积累一般是先进入移动的点，同样先稳定下来。

一种重复开采地表变形预测方法，方法流程图如图2所示，包括以下步骤：

步骤1、确定预测地表的开采工作面的顺序，并在每个工作面中设定采样点的间距，确定每个工作面的参数，包括下沉系数、水平移动系数、拐点偏距、主要影响正切值、开采影响传播角、最大下沉角、煤层倾角、工作面尺寸、采深和采厚10个参数，并根据上述10个参数获得每个采样点的走向各变形值和倾斜方向各变形值，所述的变形值包括下沉值、倾斜值、曲率值、水平移动值和水平变形值五个值；

1、下沉系数及水平移动系数

下沉系数是指在达到充分采动时，地表最大下沉量W_max与采高M之间的比值，下沉系数与煤层的埋藏深度、地质条件、采矿技术条件、煤岩的物理力学性质、工作面的开采尺寸和重复采动情况有关。不同的覆岩性质区分地表下沉系数的取值也不同，覆岩坚硬时为0.4～0.65，覆岩中硬时为0.658～0.85，覆岩软弱时为0.8～1.0。

水平移动系数是指在达到充分采动时，地表最大水平移动值U_max与最大下沉值W_max之间的比值，在一般的地质采矿条件下水平移动系数在0.1～0.4之间。

2、最大下沉角及煤层倾角

在移动盆地倾斜主断面上，地表移动盆地下山方向的影响范围扩大，最大下沉点不在采空区中央的正上方，而是向下山方向偏移，其位置用最大下沉角θ确定。非充分采动和充分采动时，在移动盆地倾斜主断面上(图3中图(a))，实测地表下沉曲线的最大下沉点在地表水平线上的投影点O与采空区中点的连线与水平线之间在煤层下山方向一侧的夹角θ，称为最大下沉角。超充分采动时(图3中图(b))，可根据下山充分采动角和上山充分采动角作两直线，其交点P至采空区中点的连线与水平线在煤层下山方向一侧的夹角θ，为此时的最大下沉角。

最大下沉角除与岩性有关外，还与煤层倾角δ有关。在缓倾斜和倾斜煤层条件下，θ值随煤层倾角的增大而减小，一般用θ=90°-kδ表示，式中的k为与岩性相关的系数。覆岩坚硬时，θ=90°-(0.7～0.8)δ；覆岩中硬时，θ=90°-(0.6～0.7)δ；覆岩软弱时，θ=90°-(0.5～0.6)δ。

3、开采影响传播角

开采影响传播角是倾向主断面特有的参数。在图4中，设A，B为实际开采边界，由于顶板的悬臂梁影响所确定的计算边界分别为C、D，其下山和上山方向的拐点偏距分别为S₁，S₂。由于煤层的倾斜，点C到上山无穷远处的点G之间煤层的半无限开采引起的下沉曲线的拐点，不位于计算边界点C的正上方地表，而是向下山方向偏移，位于O处。CO线与水平线的夹角θ₀称为开采影响传播角。半无限开采DG引起的地表下沉曲线的拐点出现在O₁点，DO₁线与水平线的夹角为θ₀。开采影响传播角θ₀一般认为与煤层倾角有关，即：θ₀=90°-kδ。式中的k为小于1的常数，一般在0.5到0.8之间。覆岩坚硬时，θ₀=90°-(0.7～0.8)δ；覆岩中硬时，θ₀=90°-(0.6～0.7)δ；覆岩软弱时，θ₀=90°-(0.5～0.6)δ。

4、主要影响正切值

充分采动时地表移动盆地平底最大下沉点至盆地边缘的距离称为主要影响半径r(如图5所示)。地表移动和破坏主要发生在-r到r之间的范围内。将x=±r的地表点与煤壁相连，其连线与水平线之间的夹角为主要影响角，其正切值tanβ称为主要影响角正切。tanβ=H/r，其中H是采深。不同性质覆岩情况下，主要影响角正切经验值分别为：覆岩坚硬时，tanβ=1.2～1.6；覆岩中硬时，tanβ=1.4～2.2；覆岩软弱时，tanβ=1.8～2.6。

5、拐点偏距

下沉曲线由凸变凹的地表分界点D称为下沉曲线的拐点(如图6所示)。由于工作面开采过程中在采空区靠近煤壁处的顶板存在着悬臂梁，使得在煤壁的O₁的正上方位置不是下沉曲线的拐点，而是向右偏移了一段距离S₀，S₀称为拐点偏距。拐点在采空区侧，则拐点偏距取正值，在煤壁侧取负值。拐点偏距与采深、覆岩岩性和煤层硬度有关。采深越大拐点偏距越大，覆岩岩性和煤层越坚硬拐点偏距越大，反之则越小。不同性质覆岩情况下，拐点偏距的取值分别为：坚硬覆岩时，S₀=(0.15～0.2)×H；中硬覆岩时，S₀=(0.1～0.15)×H；软弱覆岩时，S₀=(0.05～0.1)×H。

获得每个采样点的走向各变形值和倾斜方向各变形值，具体如下：

(1)计算每个采样点的走向各变形值如下：

计算下沉值公式：

$W\left(x\right)=\frac{W^0}{2}[\mathrm{erf}\left(\frac{\sqrt{\mathrm{\π}}}{r}x\right)+1]---\left(14\right)$

其中，W(x)表示x处的下沉值；W⁰表示走向和倾向均为充分采动时的地表最大下沉值，W⁰=mqcosδ，其中，m表示煤层开采厚度；δ表示煤层倾角；q表示下沉系数；

表示概率积分函数，其中，r表示主要影响半径，

字表示采深，tanβ表示主要影响正切值；

计算倾斜值公式：

$i\left(x\right)=\frac{W^0}{r}{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(15\right)$

其中，i(x)表示x处的倾斜值；

计算曲率值公式：

$K\left(x\right)=-\frac{2\mathrm{\π}{W}^0}{r^3}x{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(16\right)$

其中，K(x)表示x处的曲率值；

计算水平移动值公式：

$U\left(x\right)=\mathrm{bri}\left(x\right)=b{W}^0{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(17\right)$

其中，U(x)表示x处的水平移动值；b表示水平移动系数；

计算水平变形值公式：

$\mathrm{\ϵ}\left(x\right)=\mathrm{brK}\left(x\right)=-\frac{2\mathrm{\πb}{W}^0}{r^2}x{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(18\right)$

其中，ε(x)表示x处的水平变形值；

Ⅲ2Ⅳ计算每个采样点的倾斜方向各变形值如下：

计算下沉值公式：

$W\left(y\right)=\frac{W^0}{2}[\mathrm{erf}\left(\frac{\sqrt{p}}{r}y\right)+1]---\left(19\right)$

其中，W(y)表示y处的下沉值；

计算倾斜值公式：

$i\left(y\right)=\frac{W^0}{r}{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(20\right)$

其中，i(y)表示y处的倾斜值；

计算曲率值公式：

$K\left(y\right)=-\frac{2\mathrm{\π}{W}^0}{r^3}y{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(21\right)$

其中，K(y)表示y处的曲率值；

计算水平移动值公式：

$U\left(y\right)=\mathrm{bri}\left(y\right)=b{W}^0{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(22\right)$

其中，U(y)表示y处的水平移动值；

计算水平变形值公式：

$\mathrm{\ϵ}\left(y\right)=\mathrm{brK}\left(y\right)=-\frac{2\mathrm{\πb}{W}^0}{r^2}x{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(23\right)$

其中，ε(y)表示y处的水平变形值。

步骤2、根据每个采样点的走向各变形值和倾斜方向各变形值，采用概率积分法获得每个工作面的沿工作面走向各变形值和沿工作面倾斜方向各变形值，具体如下：

沿工作面走向各变形值的预测公式如下：

$\left\{\begin{array}{c}{W}^0\left(x\right)=C_{\mathrm{ym}}[W(x;t_1)-W(x-l;t_2)]\\ i^0\left(x\right)=C_{\mathrm{ym}}[i(x;t_1)-i(x-l;t_2)]\\ K^0\left(x\right)=C_{\mathrm{ym}}[K(x;t_1)-K(x-l;t_2)]\\ U^0\left(x\right)=C_{\mathrm{ym}}[U(x;t_1)-U(x-l;t_2)]\\ {\mathrm{\ϵ}}^0\left(x\right)=C_{\mathrm{ym}}[\mathrm{\ϵ}(x;t_1)-\mathrm{\ϵ}(x-l;t_2)]\\ C_{\mathrm{ym}}=\frac{W_{\mathrm{my}}^0}{W^0}\end{array}\right.---\left(1\right)$

其中，l为工作面走向长度；x为走向坐标，W(x;t₁)表示位于工作面走向方向距离选定原点的x处的下沉值；W(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处下沉值；i(x;t₁)表示位于工作面走向方向距离选定原点的x处倾斜值；i(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处倾斜值；K(x;t₁)表示位于工作面走向方向距离选定原点的x处曲率值；K(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处曲率值；U(x;t₁)表示位于工作面走向方向距离选定原点的x处水平移动值；U(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处水平移动值；ε(x;t₁)表示位于工作面走向方向距离选定原点的x处水平变形值；ε(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处水平变形值；t₁表示下山方向开采边界的10个参数、t₂表示上山方向开采边界的10个参数；W⁰(x)表示走向下沉值，i⁰(x)表示走向倾斜值，K⁰(x)表示走向曲率值，U⁰(x)表示走向水平移动值，ε⁰(x)表示走向水平变形值；表示走向充分采动倾向主断面的最大下沉值；W⁰表示走向和倾向均为充分采动时的地表最大下沉值，C_ym表示倾向采动程度系数；

沿工作面倾斜方向各变形值的预测公式为：

$\left\{\begin{array}{c}{W}^0\left(y\right)=C_{\mathrm{xm}}[W(y;t_3)-W(y-L;t_4)]\\ i^0\left(y\right)=C_{\mathrm{xm}}[i(y;t_3)-i(y-L;t_4)]\\ K^0\left(y\right)=C_{\mathrm{xm}}[K(y;t_3)-K(y-L;t_4)]\\ U^0\left(y\right)=C_{\mathrm{xm}}[U(y;t_3)-U(y-L;t_4)]\\ {\mathrm{\ϵ}}^0\left(y\right)=C_{\mathrm{xm}}[\mathrm{\ϵ}(y;t_3)-\mathrm{\ϵ}(y-L;t_4)]\\ C_{\mathrm{xm}}=\frac{W_{\mathrm{mx}}^0}{W^0}\end{array}\right.---\left(2\right)$

其中，L表示工作面倾斜长度；表示倾向充分采动时走向主断面的最大下沉值；C_xm表示走向采动程度系数；y表示倾斜方向坐标；W⁰(y)表示倾斜方向下沉值，i⁰(y)表示倾斜方向倾斜值，K⁰(y)表示倾斜方向曲率值，U⁰(y)表示倾斜方向水平移动值，ε⁰(y)表示倾斜方向水平变形值；t₃表示走向左侧开采边界的10个参数；t₄表示走向右侧开采边界的10个参数；W(y;t₃)表示位于工作面倾斜方向距离选定原点的y处的下沉值；W(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处的下沉值；i(y;t₃)表示位于工作面倾斜方向距离选定原点的y处倾斜值；i(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处倾斜值；K(y;t₃)表示位于工作面倾斜方向距离选定原点的_y处曲率值；K(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处曲率值；U(y;t₃)表示位于工作面倾斜方向距离选定原点的y处水平移动值；U(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处水平移动值；ε(y;t₃)表示位于工作面倾斜方向距离选定原点的y处水平变形值；ε(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处水平变形值；

本发明实施例中，采用概率积分法算出每个工作面开采后沿着走向和倾向方向上的各个变形值，而且在概率积分法中定义了，除下沉值外每个工作面开采形成盆地特定位置变形值的正与负。

本发明实施例中，采用矩形网格模型的等值线绘制方法对地表变形情况进行解释说明：地表下沉盆地内各点的倾斜和水平移动，除了与点的位置有关外，还与给定方向有关。设定一个开采工作平面图，煤层倾斜方向指向图下方，原点在最大下沉点外，亦可在开采边界交点上。y轴指向煤层上山方向，由y轴正方向顺时针旋转90°为x轴正方向。如若给定方向与x轴正方向平行，则φ=0°；如若给定方向和y轴正方向平行，则φ=90°。沿走向(φ=0°)和沿倾斜(φ=90°)方向的倾斜值和水平移动等值线全面积分布均为两组椭圆(如图7中图(a)图(b)、图8中图(a)图(b))，实线的一组表示正值，虚线表示负值。

曲率值和水平变形值等值线的全面积分布均为四组椭圆形(如图9中图(a)图(b)、图10中图(a)图(b))。工作面中部上方为负值，用虚线表示，煤柱上方为正值，用实线表示。在多工作面开采后地表移动变形的预测过程中可以将单个工作面开采后沿着走向和倾向的各个变形值进行代数和的叠加，在叠加的过程中充分的体现了开采过程中的抵消和加剧的实际变化过程。

步骤3、采用距离幂次反比法将每个工作面的沿工作面走向各变形值进行叠加，将每个工作面的沿工作面倾斜方向各变形值进行叠加，具体如下：

步骤3-1、按照地表开采工作面的顺序，确定第一次开采的工作面和第二次开采的工作面之间的重叠区域，将重叠区域中第一次开采的工作面中某一采样点作为圆心，以第一次开采的工作面中采样点的间距作为半径，确定搜索范围；

本发明实施例中，将A与B工作面的相重叠区域的数据进行叠加，以A工作面的数据作为母本，在B工作的数据中进行搜索，在搜过程中以A工作面的计算步长为搜索半径，即X步长和Y步长(即采样点X方向和Y方向的间距)。搜索到的所有数据利用距离幂次反比法进行相加确定Z_d；

步骤3-2、在所确定的搜索范围内对应搜索第二次开采的工作面中采样点的数据，包括沿工作面走向各变形值和沿工作面倾斜方向各变形值，并采用距离幂次反比法对搜索到的采样点数据进行叠加计算，公式如下：

$Z_d=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{Z_i}{d_i}}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{d_i}}---\left(3\right)$

其中，Z_d表示搜索到的点经过距离幂次反比法得到的修正变形值；Z_i表示搜索到的采样点变形值，包括下沉值、倾斜值、曲率值、水平移动值和水平变形值；d_i表示搜索到第二次开采工作面的采样点距离搜索范围圆心的距离；n表示搜索到的点的个数；

确定叠加后的采样点变形值：

Z=Z_x+Z_d     (4)

其中，Z表示叠加后的搜索范围圆心的变形值；Z_x表示搜索范围圆心点的原始变形值；

步骤3-3、依次将重叠区域中第一次开采的工作面中其他采样点作为圆心，以第一次开采的工作面中采样点的间距作为半径，确定新的搜索范围，并重复执行步骤3-2，直至重叠区域中第一次开采工作面的所有采样点均搜索完成；

步骤3-4、将经过叠加合成后的第一次开采的工作面和第二次开采的工作面作为一个整体，重复步骤3-1至步骤3-3与下一次开采的工作面进行叠加合成，直至所有工作面采样点的沿工作面走向各变形值和沿工作面倾斜方向各变形值均完成叠加合成；

本发明实施例中，在提取数据之前将重叠部分的B工作面中的点删除；当工作面个数大于两个工作面时，首先确定开采顺序，然后将第一次和第二次开采的工作面利用上述方法叠加后的数据作为一个整体工作面开采后的数据再与第三个工作面开采后的数据进行叠加，以此类推多个工作面开采后的沿着走向和倾向的变形值。

步骤4、将每个采样点的变形值进行合成，获得采样点实际预测变形值，具体为：

步骤4-1、对采样点的沿工作面走向下沉值和沿工作面倾斜方向下沉值进行合成，获得实际预测下沉值；

在概率积分法中，采空区地表的下沉量为W₀=mq cosδ，采空区的范围是一个矩形，矩形两个边长分别为L和D。开采影响范围为以采空区为中心向周围的发散的区域，发散区域的大小是一个影响半径的范围。则在整个影响范围内的下沉值用积分函数表示为；

$W(x,y)=W^0{\∫}_0^D{\∫}_0^L\frac{1}{r^2}{e}^{-\mathrm{\π}\frac{{(x-s)}^2+{(y-t)}^2}{r^2}}\mathrm{dtds}---\left(24\right)$

其中，t和s分别表示预测点在矩形两个边上的位置，则

$W(x,y)=W^0{\∫}_0^D{\∫}_0^L\frac{1}{r}{e}^{-\mathrm{\π}\frac{{(x-s)}^2}{r^2}}\frac{1}{r}{e}^{-\mathrm{\π}\frac{{(y-t)}^2}{r^2}}\·\mathrm{dtds}---\left(25\right)$

将上式进行变换，上式可写成；

$\begin{array}{c}W(x,y)=W^0{\∫}_0^D\frac{1}{r}{e}^{-\mathrm{\π}\frac{{(x-s)}^2}{r^2}}\mathrm{ds}{\∫}_0^L\frac{1}{r}{e}^{-\mathrm{\π}\frac{{(y-t)}^2}{r^2}}\mathrm{dt}\\ =\frac{1}{W^0}{W}^0[{\∫}_0^{\∞}\frac{1}{r}{e}^{-\mathrm{\π}\frac{{(x-s)}^2}{r^2}}\mathrm{ds}-{\∫}_D^{\∞}\frac{1}{r}{e}^{-\mathrm{\π}\frac{{(x-s)}^2}{r^2}}\mathrm{ds}]\·W^0[{\∫}_0^{\∞}\frac{1}{r}{e}^{-\mathrm{\π}\frac{{(y-t)}^2}{r^2}}\mathrm{dt}-{\∫}_L^{\∞}\frac{1}{r}{e}^{-\mathrm{\π}\frac{{(y-t)}^2}{r^2}}\mathrm{dt}]\end{array}---\left(26\right)$

上式简写为：

$W(x,y)=\frac{1}{W^0}{W}^0\left(x\right)W^0\left(y\right)---\left(5\right)$

其中，W(x,y)表示坐标为(x,y)处的采样点实际下沉值；

步骤4-2、对采样点的沿工作面走向倾斜值和沿工作面倾斜方向倾斜值进行合成，获得实际预测倾斜值；

本发明实施例中，在开采沉陷过程中水平移动值、倾斜值、曲率值和水平变形值都是因为下沉产生的。针对任意一点A(x,y)点处沿指定的方向上的四个变形值。其中α为从x轴的正向按照逆时针旋转时计算得到的角值。A(x,y)点沿α方向上的倾斜值i(x,y,α)为该点的下沉W(x,y)在α方向上的变化率，在数学计算中即为α方向上的方向导数，则有：

$i(x,y,\mathrm{\α})=\frac{\∂W(x,y)}{\∂\mathrm{\α}}=\frac{\∂W(x,y)}{\∂x}\cos \mathrm{\α}+\frac{\∂W(x,y)}{\∂y}\sin \mathrm{\α}---\left(27\right)$

上式化简后为：

$i(x,y,\mathrm{\α})=\frac{1}{W^0}{i}^0\left(x\right)W^0\left(y\right)\cos \mathrm{\α}+\frac{1}{W^0}{W}^0\left(x\right)i^0\left(y\right)\sin \mathrm{\α}---\left(28\right)$

即为：

$i(x,y,\mathrm{\α})=\frac{1}{W^0}[i^0\left(x\right)W^0\left(y\right)\cos \mathrm{\α}+W^0\left(x\right)i^0\left(y\right)\sin \mathrm{\α}]---\left(6\right)$

其中，i(x,y,α)表示点(x,y)沿α方向上的实际倾斜值；

对地表移动盆地内任意点A(x,y)，方向α取值不同其移动和变形值(除该点下沉外)也不同。α的取值为地表的移动和变形为最大时，若在A点修筑建筑物，应将建筑物的长边尽量避开最大变形值的方向，而使建筑物的短边尽量沿最大变形值的方向。

最大倾斜和水平移动方向的确定，过A(x,y)点实际倾斜出现在α=α_i处，据实际变形值的条件，应有；

$\frac{\∂i(x,y,\mathrm{\α})}{\∂\mathrm{\α}}{|}_{\mathrm{\α}={\mathrm{\α}}_i}=0---\left(7\right)$

此时获得

$\mathrm{\α}={\mathrm{\α}}_i=\arctan \frac{W^0\left(x\right)i^0\left(y\right)}{i^0\left(x\right)W^0\left(y\right)}---\left(8\right)$

其中，α_i表示倾斜变形值实际值方向；

将获得的方向α_i代入公式(6)中，即获得实际倾斜值i(x,y,α_i)；

步骤4-3、对采样点的沿工作面走向曲率值和沿工作面倾斜方向曲率值进行合成，获得实际预测曲率值；

任意的A(x,y)点沿α方向的曲率K(x,y,α)为该点的倾斜i(x,y,α)在α方向上的变化率，在数学计算中即为i(x,y,α)在α方向上的方向导数，则有；

$K(x,y,\mathrm{\α})=\frac{\∂i(x,y,\mathrm{\α})}{\∂\mathrm{\α}}=\frac{\∂i(x,y,\mathrm{\α})}{\∂x}\cos \mathrm{\α}+\frac{\∂i(x,y,\mathrm{\α})}{\∂y}\sin \mathrm{\α}---\left(29\right)$

上式化简后为：

$K(x,y,\mathrm{\α})=\frac{1}{W^0}[K^0\left(x\right)W^0\left(y\right){\cos }^2\mathrm{\α}+K^0\left(y\right)W^0\left(x\right){\sin }^2\mathrm{\α}+i^0\left(x\right)i^0\left(y\right)\sin 2\mathrm{\α}]---\left(9\right)$

其中，K(x,y,α)表示点(x,y)沿α方向上的实际曲率值；

根据实际变形值的条件，获得如下公式；

$\frac{\∂K(x,y,\mathrm{\α})}{\∂k}{|}_{\mathrm{\α}={\mathrm{\α}}_k}=0---\left(10\right)$

此时获得

$\mathrm{\α}={\mathrm{\α}}_K=\frac{1}{2}\arctan \frac{2i^0\left(x\right)i^0\left(y\right)}{K^0\left(x\right)W^0\left(y\right)-K^0\left(y\right)W^0\left(x\right)}---\left(11\right)$

其中，α_K表示曲率值实际值方向；

将α_K值代入式(9)获得实际曲率值K(x,y,α_K)；

步骤4-4、对采样点的沿工作面走向水平位移值和沿工作面倾斜方向水平位移值进行合成，获得实际预测水平位移值；

概率积分法适合于倾角小于25°的煤层，所以讨论的煤层是小倾角的，各向同性的，根据式U⁰(x)=bri⁰(x)，则水平位移与倾斜成正比的关系，获得

$U(x,y,\mathrm{\α})=\mathrm{b\γi}(x,y,\mathrm{\α})=\frac{1}{W^0}[{\mathrm{b\γi}}^0\left(x\right)W^0\left(y\right)\cos \mathrm{\α}+W^0\left(x\right)\mathrm{b\γ}{i}^0\left(y\right)\sin \mathrm{\α}]---\left(30\right)$

即：

$U(x,y.\mathrm{\α})=\frac{1}{W^0}[U^0\left(x\right)W^0\left(y\right)\cos \mathrm{\α}+U^0\left(y\right)W^0\left(x\right)\sin \mathrm{\α}]---\left(12\right)$

其中，U(x,y,α)表示点(x,y)沿α方向上的实际水平位移值；由于水平移动和倾斜成正比，所以倾斜实际值的方向α_i也是水平位移实际值出现的方向，此时，α=α_i；

步骤4-5、对采样点的沿工作面走向水平变形值和沿工作面倾斜方向水平变形值进行合成，获得实际预测水平变形值；

根据式ε⁰(x)=brK⁰(x)得；

$\begin{array}{c}\mathrm{\ϵ}(x,y,\mathrm{\α})=\frac{1}{W^0}\left\{\right[{\mathrm{\ϵ}}^0\left(x\right)W^0\left(y\right){\cos }^2\mathrm{\α}+{\mathrm{\ϵ}}^0\left(y\right)W^0\left(x\right){\sin }^2\mathrm{\α}]\\ +[U^0\left(x\right)i^0\left(y\right)+U^0\left(y\right)i^0\left(x\right)]\sin \mathrm{\α}\cos \mathrm{\α}\}\end{array}---\left(13\right)$

其中，ε(x,y,α)表示点(x,y)沿α方向上的实际水平变形值；实际水平变形也出现在α_K方向，此时，α=α_K；

步骤5、判断确定获得的采样点实际预测变形值所属范围，采取措施；

本发明实施例中，采用矩形网格模型的等值线绘制将获得的采样点实际变形值进行可视化，具体如下：

步骤5-1、根据所有采样点的经纬坐标，确定采样点在矩形网格模型中的位置；

步骤5-2、将采样点沿着一方向无折返寻找等值点并进行连线；

步骤5-3、获得每条连线上的等值点，并判断该位置的变形值是否超过表1中的要求；

表1砖混结构建筑物损坏等级

即：

当水平变形值在0～2范围内，或曲率值在0～0.2范围内，或倾斜值在0～3范围内时，则建筑物为轻微损坏，建筑物可以不修或者简单维修；所述的简单维修，例如地表用木质材料支撑等。

当水平变形值在2～4范围内，或曲率值在0.2～0.4范围内，或倾斜值在3～际范围内时，则建筑物为轻度损坏，需要加强建筑物的监测，并及时进行简单维修和加固；

当水平变形值在4～6范围内，或曲率值在0.4～0.6范围内，或倾斜值在6～10范围内时，则建筑物为中度损坏，需要实时检测、及时维修和加强防护；或者采取降低采高或采取充填开采措施，减少地表的沉陷，控制地表变形在保护标准以内；

当水平变形值大于6，或曲率值大于0.6，或倾斜值大于10时，则建筑物为严重损坏，需要人员撤离，并进行大规模维修或拆建，所述的大规模维修例如局部破坏严重，将其翻新重建或整体拆除。

相邻工作面无煤柱开采与整体单工作面开采的地表变形结果相同，以此原则检验算法的可靠性。

如图11所示，是三个走向长度和倾向长度都相等多工作面开采。要验证的结果是三个工作面依次全采后的地表移动变形情况与一次全采后的地表移动变形进行对比。三个工作面的名称分别为1号、2号和3号，并将三个工作面总体的理想大工作面起名为0号工作面。

分别用单工作面开采0号工作面和多工作开采1号、2号和3号工作面叠加后的预测结果，绘制下沉等值线图(图12中图(a)图(b))、倾斜值等值线图(图13中图(a)图(b))、曲率值等着线图(图14中图(a)图(b))、水平变形等值线图(图15中图(a)图(b))和水平移动等值线图(图16中图(a)图(b))进行比较分析。

从图12中可以看出多工作面开采下沉值叠加结果与单工作面开采下沉值结果的最大值和影响范围相同。从图13中可以看出多工作面开采倾斜值叠加结果与单工作面开采倾斜值结果的最大值和影响范围相同。从图14中可以看出多工作面开采曲率值叠加结果与单工作面开采曲率值结果的最大值和影响范围相同。从图15中可以看出多工作面开采水平变形值叠加结果与单工作面开采水平变形值结果的最大值和影响范围相同。从图16中可以看出多工作面开采水平移动值叠加结果与单工作面开采水平移动值结果的最大值和影响范围相同。

根据以上的结果分析得知，重复开采地表移动变形算法是正确可靠的，可以用在重复开采地表移动变形的预测。

Claims (2)

1.一种重复开采地表变形预测方法，其特征在于，包括以下步骤：

步骤1、确定预测地表的开采工作面的顺序，并在每个工作面中设定采样点的间距，确定每个工作面的参数，包括下沉系数、水平移动系数、拐点偏距、主要影响正切值、开采影响传播角、最大下沉角、煤层倾角、工作面尺寸、采深和采厚10个参数，并根据上述10个参数获得每个采样点的走向各变形值和倾斜方向各变形值，所述的变形值包括下沉值、倾斜值、曲率值、水平移动值和水平变形值五个值；

步骤2、根据每个采样点的走向各变形值和倾斜方向各变形值，采用概率积分法获得每个工作面的沿工作面走向各变形值和沿工作面倾斜方向各变形值，具体如下：

沿工作面走向各变形值的预测公式如下：

$\left\{\begin{array}{c}{W}^0\left(x\right)=C_{\mathrm{ym}}[W(x;t_1)-W(x-l;t_2)]\\ i^0\left(x\right)=C_{\mathrm{ym}}[i(x;t_1)-i(x-l;t_2)]\\ K^0\left(x\right)=C_{\mathrm{ym}}[K(x;t_1)-K(x-l;t_2)]\\ U^0\left(x\right)=C_{\mathrm{ym}}[U(x;t_1)-U(x-l;t_2)]\\ {\mathrm{\ϵ}}^0\left(x\right)=C_{\mathrm{ym}}[\mathrm{\ϵ}(x;t_1)-\mathrm{\ϵ}(x-l;t_2)]\\ C_{\mathrm{ym}}=\frac{W_{\mathrm{my}}^0}{W^0}\end{array}\right.---\left(1\right)$

其中，l为工作面走向长度；x为走向坐标，W(x;t₁)表示位于工作面走向方向距离选定原点的x处的下沉值；W(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处下沉值；i(x;t₁)表示位于工作面走向方向距离选定原点的x处倾斜值；i(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处倾斜值；K(x;t₁)表示位于工作面走向方向距离选定原点的x处曲率值；K(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处曲率值；U(x;t₁)表示位于工作面走向方向距离选定原点的x处水平移动值；U(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处水平移动值；ε(x;t₁)表示位于工作面走向方向距离选定原点的x处水平变形值；ε(x-l;t₂)表示位于工作面走向方向距离选定原点的x-l处水平变形值；t₁表示下山方向开采边界的10个参数、t₂表示上山方向开采边界的10个参数；W⁰(x)表示走向下沉值，i⁰(x)表示走向倾斜值，K⁰(x)表示走向曲率值，U⁰(x)表示走向水平移动值，ε⁰(x)表示走向水平变形值；

表示走向充分采动倾向主断面的最大下沉值；W⁰表示走向和倾向均为充分采动时的地表最大下沉值，C_ym表示倾向采动程度系数；

沿工作面倾斜方向各变形值的预测公式为：

$\left\{\begin{array}{c}{W}^0\left(y\right)=C_{\mathrm{xm}}[W(y;t_3)-W(y-L;t_4)]\\ i^0\left(y\right)=C_{\mathrm{xm}}[i(y;t_3)-i(y-L;t_4)]\\ K^0\left(y\right)=C_{\mathrm{xm}}[K(y;t_3)-K(y-L;t_4)]\\ U^0\left(y\right)=C_{\mathrm{xm}}[U(y;t_3)-U(y-L;t_4)]\\ {\mathrm{\ϵ}}^0\left(y\right)=C_{\mathrm{xm}}[\mathrm{\ϵ}(y;t_3)-\mathrm{\ϵ}(y-L;t_4)]\\ C_{\mathrm{xm}}=\frac{W_{\mathrm{mx}}^0}{W^0}\end{array}\right.---\left(2\right)$

其中，L表示工作面倾斜长度；表示倾向充分采动时走向主断面的最大下沉值；C_xm表示走向采动程度系数；y表示倾斜方向坐标；W⁰(y)表示倾斜方向下沉值，i⁰(y)表示倾斜方向倾斜值，K⁰(y)表示倾斜方向曲率值，U⁰(y)表示倾斜方向水平移动值，ε⁰(y)表示倾斜方向水平变形值；t₃表示走向左侧开采边界的10个参数；t₄表示走向右侧开采边界的10个参数；W(y;t₃)表示位于工作面倾斜方向距离选定原点的y处的下沉值；W(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处的下沉值；i(y;t₃)表示位于工作面倾斜方向距离选定原点的y处倾斜值；i(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处倾斜值；K(y;t₃)表示位于工作面倾斜方向距离选定原点的y处曲率值；K(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处曲率值；U(y;t₃)表示位于工作面倾斜方向距离选定原点的y处水平移动值；U(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处水平移动值；ε(y;t₃)表示位于工作面倾斜方向距离选定原点的y处水平变形值；ε(y-L;t₄)表示位于工作面倾斜方向距离选定原点的y-L处水平变形值；

步骤3、采用距离幂次反比法将每个工作面的沿工作面走向各变形值进行叠加，将每个工作面的沿工作面倾斜方向各变形值进行叠加，具体如下：

步骤3-1、按照地表开采工作面的顺序，确定第一次开采的工作面和第二次开采的工作面之间的重叠区域，将重叠区域中第一次开采的工作面上某一采样点作为圆心，以第一次开采的工作面中采样点的间距作为半径，确定搜索范围；

步骤3-2、在所确定的搜索范围内对应搜索第二次开采的工作面中采样点的数据，包括沿工作面走向各变形值和沿工作面倾斜方向各变形值，并采用距离幂次反比法对搜索到的采样点数据进行叠加计算，公式如下：

$Z_d=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{Z_i}{d_i}}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{d_i}}---\left(3\right)$

其中，Z_d表示搜索到的点经过距离幂次反比法得到的修正变形值；Z_i表示搜索到的采样点变形值，包括下沉值、倾斜值、曲率值、水平移动值和水平变形值；d_i表示搜索到第二次开采工作面的采样点距离搜索范围圆心的距离；n表示搜索到的点的个数；

确定叠加后的采样点变形值：

Z=Z_x+Z_d      (4)

其中，Z表示叠加后的搜索范围圆心的变形值；Z_x表示搜索范围圆心点的原始变形值；

步骤3-3、依次将重叠区域中第一次开采的工作面上其他采样点作为圆心，以第一次开采的工作面中采样点的间距作为半径，确定新的搜索范围，并重复执行步骤3-2，直至重叠区域中第一次开采工作面的所有采样点均搜索完成；

步骤3-4、将经过叠加合成后的第一次开采的工作面和第二次开采的工作面作为一个整体，重复步骤3-1至步骤3-3与下一次开采的工作面进行叠加合成，直至所有工作面采样点的沿工作面走向各变形值和沿工作面倾斜方向各变形值均完成叠加合成；

步骤4、将每个采样点的变形值进行合成，获得采样点实际预测变形值，具体为：

步骤4-1、对采样点的沿工作面走向下沉值和沿工作面倾斜方向下沉值进行合成，获得实际预测下沉值，公式如下：

$W(x,y)=\frac{1}{W^0}{W}^0\left(x\right)W^0\left(y\right)---\left(5\right)$

其中，W(x,y)表示坐标为(x,y)处的采样点实际下沉值；

步骤4-2、对采样点的沿工作面走向倾斜值和沿工作面倾斜方向倾斜值进行合成，获得实际预测倾斜值，公式如下：

$i(x,y,\mathrm{\α})=\frac{1}{W^0}[i^0\left(x\right)W^0\left(y\right)\cos \mathrm{\α}+W^0\left(x\right)i^0\left(y\right)\sin \mathrm{\α}]---\left(6\right)$

其中，i(x,y,α)表示点(x,y)沿α方向上的实际倾斜值；

根据实际变形值的条件，获得如下公式；

$\frac{\∂i(x,y,\mathrm{\α})}{\∂\mathrm{\α}}{|}_{\mathrm{\α}={\mathrm{\α}}_i}=0---\left(7\right)$

此时获得

$\mathrm{\α}={\mathrm{\α}}_i=\arctan \frac{W^0\left(x\right)i^0\left(y\right)}{i^0\left(x\right)W^0\left(y\right)}---\left(8\right)$

其中，α_i表示倾斜变形值实际值方向；

将获得的方向α_i代入公式(6)中，即获得实际倾斜值i(x,y,α_i)；

步骤4-3、对采样点的沿工作面走向曲率值和沿工作面倾斜方向曲率值进行合成，获得实际预测曲率值，公式如下：

$K(x,y,\mathrm{\α})=\frac{1}{W^0}[K^0\left(x\right)W^0\left(y\right){\cos }^2\mathrm{\α}+K^0\left(y\right)W^0\left(x\right){\sin }^2\mathrm{\α}+i^0\left(x\right)i^0\left(y\right)\sin 2\mathrm{\α}]---\left(9\right)$

其中，K(x,y,α)表示点(x,y)沿α方向上的实际曲率值；

根据实际变形值的条件，获得如下公式；

$\frac{\∂K(x,y,\mathrm{\α})}{\∂k}{|}_{\mathrm{\α}={\mathrm{\α}}_k}=0---\left(10\right)$

此时获得

$\mathrm{\α}={\mathrm{\α}}_K=\frac{1}{2}\arctan \frac{2i_0\left(x\right)i_0\left(y\right)}{K_0\left(x\right)W_0\left(y\right)-K_0\left(y\right)W_0\left(x\right)}---\left(11\right)$

其中，α_K表示曲率值实际值方向；

将α_K值代入式(9)获得实际曲率值K(x,y,α_K)；

步骤4-4、对采样点的沿工作面走向水平位移值和沿工作面倾斜方向水平位移值进行合成，获得实际预测水平位移值，公式如下：

$U(x,y.\mathrm{\α})=\frac{1}{W^0}[U^0\left(x\right)W^0\left(y\right)\cos \mathrm{\α}+U^0\left(y\right)W^0\left(x\right)\sin \mathrm{\α}]---\left(12\right)$

其中，U(x,y,α)表示点(x,y)沿α方向上的实际水平位移值；此时，α=α_i；

步骤4-5、对采样点的沿工作面走向水平变形值和沿工作面倾斜方向水平变形值进行合成，获得实际预测水平变形值，公式如下：

$\begin{array}{c}\mathrm{\ϵ}(x,y,\mathrm{\α})=\frac{1}{W^0}\left\{\right[{\mathrm{\ϵ}}^0\left(x\right)W^0\left(y\right){\cos }^2\mathrm{\α}+{\mathrm{\ϵ}}^0\left(y\right)W^0\left(x\right){\sin }^2\mathrm{\α}]\\ +[U^0\left(x\right)i^0\left(y\right)+U^0\left(y\right)i^0\left(x\right)]\sin \mathrm{\α}\cos \mathrm{\α}\}\end{array}---\left(13\right)$

其中，ε(x,y,α)表示点(x,y)沿α方向上的实际水平变形值；此时，α=α_K；

步骤5、判断确定获得的采样点实际预测变形值所属范围，采取措施，具体为：

当水平变形值在0～2范围内，或曲率值在0～0.2范围内，或倾斜值在0～3范围内时，则建筑物为轻微损坏，需要加强建筑物的监测，并及时进行加固；

当水平变形值在2～4范围内，或曲率值在0.2～0.4范围内，或倾斜值在3～6范围内时，则建筑物为轻度损坏，需要降低采高或采取充填开采措施，减少地表的沉陷，控制地表变形在保护标准以内；

当水平变形值在4～6范围内，或曲率值在0.4～0.6范围内，或倾斜值在6～10范围内时，则建筑物为中度损坏，需要降低采高或采取充填开采措施，减少地表的沉陷，控制地表变形在保护标准以内；

当水平变形值大于6，或曲率值大于0.6，或倾斜值大于10时，则建筑物为严重损坏，需要人员撤离，并进行维修或拆建。

2.根据权利要求1所述的重复开采地表变形预测方法，其特征在于，步骤1所述的根据10个参数获得每个采样点的走向各变形值和倾斜方向各变形值，具体如下：

(1)计算每个采样点的走向各变形值如下：

计算下沉值公式：

$W\left(x\right)=\frac{W^0}{2}[\mathrm{erf}\left(\frac{\sqrt{\mathrm{\π}}}{r}x\right)+1]---\left(14\right)$

其中，W(x)表示x处的下沉值；W⁰表示走向和倾向均为充分采动时的地表最大下沉值，W⁰=mqcosα，其中，m表示煤层开采厚度；α表示煤层倾角；q表示下沉系数；表示概率积分函数，其中，r表示主要影响半径，H表示采深，tanβ表示主要影响正切值；

计算倾斜值公式：

$i\left(x\right)=\frac{W^0}{r}{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(15\right)$

其中，i(x)表示x处的倾斜值；

计算曲率值公式：

$K\left(x\right)=-\frac{2\mathrm{\π}{W}^0}{r^3}x{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(16\right)$

其中，K(x)表示x处的曲率值；

计算水平移动值公式：

$U\left(x\right)=\mathrm{bri}\left(x\right)=b{W}^0{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(17\right)$

其中，U(x)表示x处的水平移动值；b表示水平移动系数；

计算水平变形值公式：

$\mathrm{\ϵ}\left(x\right)=\mathrm{brK}\left(x\right)=-\frac{2\mathrm{\πb}{W}^0}{r^2}x{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(18\right)$

其中，ε(x)表示x处的水平变形值；

(2)计算每个采样点的倾斜方向各变形值如下：

计算下沉值公式：

$W\left(y\right)=\frac{W^0}{2}[\mathrm{erf}\left(\frac{\sqrt{p}}{r}y\right)+1]---\left(19\right)$

其中，W(y)表示y处的下沉值；

计算倾斜值公式：

$i\left(y\right)=\frac{W^0}{r}{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(20\right)$

其中，i(y)表示y处的倾斜值；

计算曲率值公式：

$K\left(y\right)=-\frac{2\mathrm{\π}{W}^0}{r^3}y{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(21\right)$

其中，K(y)表示y处的曲率值；

计算水平移动值公式：

$U\left(y\right)=\mathrm{bri}\left(y\right)=b{W}^0{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(22\right)$

其中，U(y)表示y处的水平移动值；

计算水平变形值公式：

$\mathrm{\ϵ}\left(y\right)=\mathrm{brK}\left(y\right)=-\frac{2\mathrm{\πb}{W}^0}{r^2}x{e}^{-\mathrm{\π}\frac{y^2}{r^2}}---\left(23\right)$

其中，ε(y)表示y处的水平变形值。

Images