基于数值积分法的地下开采地表变形预测算法

孙士国，文小艳，胡业浩

(1. 华北科技学院，河北 三河065201； 2. 临沂大学机械工程学院，山东 临沂276005)

1 引言

随着国民经济的高速发展，各种能源逐渐被消耗，其中，矿产资源消耗量较高，合理开采与利用矿产资源，对社会与人类未来发展具有重大意义。由于大量的矿山开采活动，开采区域范围内与矿区内部的稳定性遭到破坏，地表沉陷问题日渐发生，岩层存在移动、变形等情况，随着开采范围逐渐扩大，岩层移动变形情况加剧，最终出现地表塌陷、崩塌等严重问题。

所谓的地下开采地表变形，即为地下资源被开采结束后，开采作业面与采空区边缘的原始应力平衡均被打破，因此，岩体或地表均会出现移动变形等问题，最终反映至地表，即构成了地表变形。在19世纪时，就出现过地表变形与沉陷最终使建筑物受到损伤的现象，这种现象由地下开采引起，为此，各开采单位开始重视地下开采造成的负面影响。由于地表变形的影响十分复杂，且不同地质与地下水分布均会导致变形情况的不同，所以，选取适当的预测方法，才能使地表变形预测更加精准。许多学者对变形预测方法进行研究，例如童凯等，研究基于盐岩腔体收敛过程的地表变形动态预测算法，但该算法预测不同开采深度时的地表水平位移误差较大。同时还有钟国强等，研究基于SFLA-GRNN模型的基坑地表最大沉降预测算法，但该算法针对相应开采宽度情况下地表竖向位移预测精度不高，因此，本文研究基于数值积分法的地下开采地表变形预测算法。

数值积分法是现阶段矿产开采过程中应用较为广泛的预测算法。该算法的应用能够使地表沉陷变形得到有效预测，并改善地下开采对地表变形的影响。本文依据数值积分法的基本原理，构建地表变形预测函数，使地表下沉值与地表相关变形值得到计算，通过该算法旋转坐标系与处理拐点偏距，实现地表变形预测，并通过预测程序，输出最终的预测结果。

2 基于数值积分法的地下开采地表变形预测算法

2.1 数值积分法基本原理

在随机介质原理的基础上，构建新型的开采地表变形预测函数，即为数值积分法。该算法是较为成熟的变形预测算法，通过式(1)计算地表下沉值

(1)

式(1)中，坐标(，)的下沉量由(，)表示；地表最大下沉值由表示；地下开采范围由表示；开采影响传播角由表示，表示下沉量，、分别表示地下开采宽度与深度，表示某段开采岩石。

通过式(2)计算地表某一点沿方向与方向的倾斜

(2)

式(2)中，(，)、(，)分别表示地表点(，)随着，方向的倾斜，∂表示所沿方向。

通过式(3)计算地表变形曲率

(3)

式(3)中，(，)、(，)依次表示地表点(，)随着，方向的曲率，表示曲率值。

通过式(4)计算地表变形水平位移

(4)

式(4)中，(，)、(，)依次表示地表点(，)随着，方向的水平位移，表示水平位移值。

通过式(5)计算地表变形的水平变形与扭曲值

(5)

式(5)中，(，)、(，)依次表示地表点(，)随着，方向的水平变形与扭曲值，表示水平变形与扭曲值。

求解地表变形的剪切值，该值可通过下沉与水平位移的求导获得，通过式(6)表示

(6)

式(6)中，(，)、(，)依次为地表点(，)的扭曲与剪切，表示扭曲值，表示剪切值。同时，该公式依据直角减小的剪切值为正进行设定。通过以上公式，即能够实现地下开采地表变形各参数的计算。

2.2 基于数值积分法的算法实现

221 坐标系旋转调整

通过以上数值积分法在计算地表变形参数时，通常采用轴平行于走向主断面、轴平行于倾向主断面的(，)坐标系，但现实情况下的坐标数据会出现不同的坐标形式，因此，需在预测地表变形前转换坐标位置，使作业面中的角点坐标调整为相对坐标，再进行地表变形预测。通过式(7)表示坐标旋转变换

=′cos+′sin

=′sin-′cos

(7)

式(7)中，新坐标点由、表示，而原始坐标点则由′、′表示；原坐标系′轴逆时针旋转至新坐标系轴的角度值则由表示，同时还表示走向方向方位角。通过Rotate函数实现作业面走向方向坐标系的转换，完成坐标系旋转的调整，为后续地表变差预测提供全面的坐标系数值。

222 拐点偏距处理

当悬臂作用力下引起的拐点出现偏移，偏移的长度即为拐点偏距，但地下开采地表变形情况下，通常存在任意多边形变化，因此拐点偏距会出现较大影响。所以，根据数值积分法计算得到作业面数据，能够得出角点沿作业面走向、倾角方向中的拐点偏距、。

通过加减的方式即能够计算走向方向内的拐点偏距，但倾斜方向存在资源层倾斜影响，因此需将角点的坐标调整为式(8)。

=-cos-(+sin)·cot

(8)

式(8)中，调整后的坐标由表示；调整前的坐标由表示；该角点区域的采深由表示；煤层倾角由表示。

223 地表变形预测

通过将原始作业面调整为计算作业面，实现作业面的分割预计。本文选取剖面线填充算法分割作业面，该算法能够快速有效实现作业面分割，并将分割后的作业面转化为一系列的小矩形以供后续计算，之后将小矩形叠加形成地表变形预计值。通过图1表示预测的步骤。

图1 地表变形预测流程

2.3 地下开采地表变形预测程序设计

依据数值积分法的计算原理，通过VC6.0设计地下开采地表变形计算分析程序，通过图2表示该程序设计结果。

图2 地表变形预测程序图

在图2中，前处理单元负责编写TXT格式的计算文件，输入地表变形预测所需条件等参数，之后启动程序，通过数值积分法计算地表变形参数并开始变形预测。预测结束后进入后处理单元，依据所测地表变形数值，保存成为Matlab图形文件。同时还可依据预测变形结果绘制等值线图，并通过AutoCAD软件内的DXF文件形式输出，还可以采用Matlab三维图形文件输出地表变形结果。

3 实例分析

通过仿真形式验证本文算法的预测性能，并选取文献[5]基于盐岩腔体收敛过程的地表变形动态预测算法、文献[6]基于SFLA-GRNN模型的基坑地表最大沉降预测算法与本文算法进行对比。

选取10个地表测试点，分析应用本文算法后地表变形数据的实测值与预测值误差，分析结果如表1所示。

表1 预测误差分析

根据表1可知，经本文算法预测后的下沉值中，仅有位置4、位置9存在较小误差，且最高误差仅为2mm；而经预测的曲率值仅在位置4、位置5与位置10存在误差；在倾斜值的预测中，出现误差的情况依然相对较小，最高仅为0.002m，由此说明本文算法的预测结果具有较小的误差。

分析不同开采深度情况下，地表竖向位移曲线，分析应用本文算法后的预测效果，实验结果如图3所示。

图3 不同开采深度竖向位移情况

根据图3可知，开采深度为20m时，竖向位移的变化要高于40m与60m的变化，而开采深度为60m时，竖向位移的变化最低，说明当开采深度越大，竖向位移的变化越低。同时，通过本文算法预测地表竖向位移变形过程中，本文预测值与竖向位移实测值十分接近，未出现较大变化，因此，本文算法对竖向位移变形预测效果较高。

分析不同开采深度情况下，地表水平位移曲线，分析应用本文算法后的预测效果，实验结果如图4所示。

图4 不同开采深度水平位移情况

根据图4可知，开采深度越高，地表的水平位移变形波动越小，当开采深度为20m时，水平位移约在0.2m～-0.2m之内波动，当开采深度达到60m时，水平位移约在0.05m～-0.05m之内波动，同时，经本文算法预测后与实际开采情况十分符合，依然未出现较大误差，因此，本文算法预测地表水平位移时误差较低。

分析不同开采宽度情况下，本文算法对地表竖向位移的预测情况，分析结果如图5所示。

图5 不同开采宽度竖向位移预测情况

根据图5可知，当开采宽度处于10m时，地表竖向位移变形最大，而开采宽度为30m时，地表竖向位移变形最小。采用本文算法对地表竖向位移进行预测，可知预测结果与实测结果均保持一致，所以，本文算法能够实现不同开采宽度的竖向位移变形值预测。

分析开采宽度为20m时，三种算法对地表竖向位移的预测情况，分析结果如图6所示。

图6 20m开采宽度时不同方法预测结果

根据图6可知，在处于20m的开采宽度时，实际竖向位移变形约在-0.2m～-0.8m之间波动，当三种算法开始预测后，文献[6]算法的预测结果与实测值存在较大偏离，该算法预测后，地表竖向位移约在-0.2m～-0.7m之间波动，而文献[5]算法地表竖向位移变形波动约在-0.2m～-0.65m以内，预测结果明显与实际结果不符，经本文算法预测后，地表的竖向位移未出现较大差别，不同开采区域的变形波动预测值均与实测值相符，因此，本文算法预测结果较为精准。

分析开采宽度为20m情况下，三种算法对地表水平位移的预测情况，分析结果如图7所示。

图7 20m宽度下不同算法水平位移预测结果

根据图7可知，开采深度为20m时，实际地表水平位移变化约在0.2m～-0.15m之间，而文献[6]算法的地表水平位移预测值约在-0.02m～0.12m之间波动，而文献[5]算法的预测结果也与实际结果存在较大误差，经本文算法预测后，得到的地表水平位移结果与实测结果并未存在较大误差，因此，本文预测结果较为精准。

4 结论

本文研究基于数值积分法的地下开采地表变形预测算法，通过数值积分法的基本原理计算地表变形的下沉值以及地表下沉后的变化值，并通过该算法调整坐标系位置，实现旋转坐标系，并完成拐点偏距处理，以调整后的坐标系与拐点偏距以及计算后的地表变形值为基础，完成地下开采地表变形的预测，并构建地下开采地表变形预测程序，实现地表变形预测结果的输出。经实验验证该算法的预测性能可知，在不同开采深度、开采宽度情况下，本文算法均能够实现精准的地表竖向、水平位移预测，预测结果要高于其它算法。在未来研究阶段，可依据现有研究，对地表变形预测算法继续优化，实现精度更高的变形预测。