DOI: 10.11817/j.issn.1672-7207.2018.02.023

采空区复杂顶板的三维分形特性及曲面模拟

付建新^{1, 2}，宋卫东^{1, 2}，谭玉叶^{1, 2}

(1. 北京科技大学 土木与资源工程学院，北京，100083;

2. 金属矿山高效开采与安全教育部重点实验室，北京，100083)

摘要：以石人沟铁矿采空区探测为工程背景，首先基于分形理论，对采空区复杂顶板的三维分形特性进行定量研究，得到对应的分维值。针对采空区顶板具有自相似性的特点，应用参数方程分形插值迭代函数系统，对采空区顶板进行模拟分析。在现场有限的已知点的基础上，对采空区顶板曲面进行模拟重构。研究结果表明：采空区顶板具有明显的分形特性，顶板的三维分维值综合反映了边界线的复杂程度、暴露面积及垂直方向上起伏的剧烈程度，定量地表征了顶板的非线性程度；最终模拟重构的精度与参数c_i密切相关，c_i越小模拟精度越大，越接近实际情况。采空区的分形特性从本质上反映了采空区边界的非线性特征。

关键词：采空区；复杂顶板；分形理论；分形插值；曲面模拟

中图分类号：TD803             文献标志码：A         文章编号：1672-7207(2018)02-0431-08

Three dimensional fractal characteristics of complex roof of gob and its surface modelling

FU Jianxin^{1, 2}, SONG Weidong^{1, 2}, TAN Yuye^{1, 2}

(1. School of Civil and Resource Engineering,University of Science and Technology Beijing, Beijing 100083, China;

2. State Key Laboratory of High-Efficient Mining and Safety of Metal Mines,Ministry of Education, Beijing 100083, China)

Abstract: Taking gob detection of Shirengou iron mine as the engineering background, based on the fractal theory, the fractal characteristics of the complex roof of gob were quantitatively studied, and the corresponding fractal dimension values were obtained. According to the characteristics of self-similarity of the roof of gob, the parameter equation and iterated function system were applied to the simulation and analysis of the roof. On the basis of the finite known points of the field, the surface of the gob roof was simulated and reconstructed. The results show that the roof of the gob has obvious fractal characteristics. The three-dimensional fractal dimension of the roof reflects the complexity of the boundary line, the exposed area and the degree of fluctuation in the vertical direction, which quantitatively represents the nonlinear degree of the roof. The final precision of the simulation is closely related to the parameters of c_i, The smaller the c_i, the greater the simulation accuracy, and the closer the simulation surface is to the actual situation. The fractal characteristics of the gob reflect the nonlinear characteristics of the gob boundary in essence.

Key words: gob; complex roof; fractal theory; fractal interpolation; surface modeling

采空区稳定性分析研究一直是一个热点问题^[1-4]，也是一个极其复杂的问题。工程实际表明，采空区的稳定性与复杂程度有着密切的关系，目前，对于采空区复杂程度没有有效的定量指标进行表述，大多通过定性描述进行复杂程度的分级^[5-8]。地质活动的剧烈和开采过程中岩体爆破断裂过程的不确定性和非线性，造成了采空区本质上的复杂性。激光精细探测技术的应用，实现了采空区的三维空间可视化^[9]，并可进行较精确的定性描述和简单的定量表征，如采空区体积、暴露面积等^[10-11]，但由于采空区本质上的复杂性，这些常规的研究方法存在局限性，采用非线性的理论方法深入研究采空区空间特征及边界特性，是未来发展的必然。实际的顶板则是典型的空间曲面，复杂多样，无法用常规的几何理论进行定量描述。虽然，通过精细探测和三维模型重构可实现采空区顶板的可视化和定性描述，但无法揭示顶板的内在规律。由于地下环境的复杂性，大量采空区人员与设备无法进入，因此，无法采用精细探测，只能采用物探或者钻探，而由于场地的限制，只能获得很少的数据，需要采用合适的方法对未知点进行模拟预测。本文作者首先基于分形理论，对采空区三维顶板的分形特性进行研究，然后，基于分形插值理论对采空区边界进行模拟及预测，通过有限已知点，实现采空区顶板的插值模拟，为采空区稳定性分析及控制提供参考。

1  分形理论概述

MANDELBROT在研究海岸线长度时提出了分形的概念^[12]，目前最常用的是盒维数^[13]。

盒维数又称计盒维数，设C是属于Rⁿ的非空子集，尺度为r且能完全覆盖C的正方形或立方体，所需的最少数量记为N_r(C)，则计盒维数D_B(C)为

            (1)

当式(1)的极限不收敛时，应分别计算上极限和下极限，即顶盒维数和底盒维数，根据维数基本定义可知，顶、底盒维数分别为

           (2)

若与C相交的边长为1/2ⁿ的正方形盒子数目为N_n(C)，则盒维数D_B为

            (3)

对于单位长度的正方体，若用边长r=1/k(k=1，2，3，…，n)的小立方体进行覆盖，则需要的数目N(r)为

            (4)

对于正方体，当用r=1/k(k=1，2，3，…，n)的立方体盒子进行覆盖时，数目为

            (5)

继续进行扩展到D维物体，盒子的数目N(r)与边长r的关系为

               (6)

式中：D为采空区顶板三维分形的计盒维数。

对式(6)取对数变换，可得

               (7)

该方法采用正方形、立方体、圆形或球体等具有可定量化尺度的图形去覆盖或近似分形对象，并记录对应尺度的数量N(r)。将图形尺度按照一定比例缩小，得到对应的N(r)。绘制N(r)与r双对数散点图并进行拟合，若具有明显的线性特征，即无标度特性，则说明具有分形特征。由分形定义可知

               (8)

对式(8)取对数变换，得

            (9)

线性拟合的直线斜率即为分维值。

2  采空区自相似性分析

自分形理论诞生以来，大量学者证明了矿床具有很好的分形特性^[14]，自相似性不仅存在于宏观空间分布中，而且存在于微观元素富集及品位分布中，甚至矿床形成的时空演化也具有自相似性。杨成杰等^[15-16]详细地论述了矿床中存在的自相似现象，体现在以下方面：1) 成矿作用的时间演变；2) 区域矿产分布；3) 矿化的空间分布；4) 矿石矿物中的元素分布；5) 矿床构造形迹。

采空区是开采矿体后遗留的结构，因此，与矿体有着相同的构造，相似的空间形态。另外，在开采过程中，频繁的爆破采动应力使采空区周围岩体充满了裂隙。研究表明，金属矿山采空区在空间结构域构造上具有以下共性：

1) 隐伏性强，空间分布规律性差；

2) 为最大限度地回收矿体资源，并保持较小的贫化率，总是尽可能地使开采边界接近矿体边界，并避免废石的混入。因此，采空区的空间形态与矿体有很高的相似度；

3) 在开采过程中，往往经过了大规模的爆破震动，岩体充满原生及次生裂隙，采空区边界是典型的爆破断裂边界。

综上所述，采空区在空间形态上与矿体具有高度的相似性，因此，从理论上讲采空区具有分形特性。

3  采空区三维模型及空间信息的获取

3.1  采空区空间形态的数据化

精确的采空区模型是准确进行采空区分维计算的基础，目前常用的三维探测系统为CMS空区探测系统。测量得到的是边界各点的三维坐标数据，可以直接被多种三维建模软件(如SURPAC、3DMine等)调用进行后处理，从而快速对采空区进行三维重构，从中获取点、线、面、体等空间信息。

采空区精细探测及建模技术在石人沟铁矿中得到了应用，得到了采空区的表面轮廓数据点云集，进行三维建模，得到三维实体图如图1所示。

图1  斜井采区39号采空区实体图

Fig. 1  Solid model of No.39 gob in Xiejing area

借助SURPAC等三维建模软件，可得该采空区的体积、顶板暴露面积、长、宽、高等几何参数，对于该采空区，体积为7 701 m³，顶板暴露面积为275 m²，长×宽×高为27.5 m×10.0 m×28.0 m。

虽然通过后期建模软件从采空区精细模型中得到采空区顶板若干定量指标，但这些指标并不能全面地反映采空区的复杂程度。

3.2  采空区顶板三维要素提取及模型重构

通过激光探测得到了采空区整个表面的数据点三维数据，需要去掉不需要的点集，保留顶板数据，进行三维重构。仍以石人沟铁矿斜井采区39号采空区顶板进行对比分析。图2所示为斜井采区39号采空区顶板三维模型图。

39号采空区顶板长度L_x基本和L_y相等，但z方向上则波动较大，具有明显的高度差。通过软件查询，可得到顶板上任意点的三维坐标。

图2  采空区顶板三维模型图

Fig. 2  Three-dimensional models of roof of gob

4  采空区复杂顶板三维分形特性计算

通常意义上的盒维数可较好地求解二维平面分形问题，但对于三维不规则曲面，往往需要进行切剖面或投影等处理，从而导致某一维度分维值的压缩或缺失，如将顶板投影到xOy平面，进行盒维数计算，得到xOy平面的分维数，基于欧氏几何理论体系设z方向维数为1，进行相加，得到顶板的三维盒维数。由于顶板在z方向上并不是均匀分布，因此该方法不能反映z方向上的复杂程度，得到的结果不能准确揭示顶板内在规律。

三维盒维数计算的关键问题是如何准确得到完全覆盖z方向上的立方体个数。针对采空区顶板三维形态的特征，提出了xOy-z步骤计算法。算法的基本过程如下。

1) 将顶板三维模型向xOy平面投影，确保顶板上的点全部都可以落到平面之上，不出现重叠与覆盖现象。采用边长为r的网格进行覆盖，确定包含顶板的非空网格r_i,j。图3所示为xOy平面投影网格覆盖图。图3中阴影为非空网格。

2) 确定非空盒子r_i,j内顶板上的高度坐标z_i,j，其中1≤i，j≤n-1。则r_i,j内完全覆盖z方向上的边长为r的立方体个数为

      (10)

式中：Int(x)为取整函数，则覆盖整个顶板的立方体的个数为

             (11)

图4所示为三维顶板覆盖立方体个数计算示意图。由图4可见：xOz平面网格覆盖图显示了第j行上，各列非空网格在z方向上完全覆盖顶板的立方体数量。图4(b)中散点填充域为顶板。

图3  xOy平面投影网格覆盖图

Fig. 3  Grid coverage map of xOy plane projection

图4  三维顶板覆盖立方体个数计算示意图

Fig. 4  Schematic diagram of covered cube of three-dimensional roof

以1/2倍率缩小立方体的尺度，重复上述算法，得到一系列N_n(r)~r数据对，绘制双对数曲线，进行线性拟合，得到函数关系式，即

           (12)

采用该方法进行求解时，每个步骤都是精确的数据，因此，得到的分维数接近真实的分维值。

进行三维计盒维数计算前，首先要确定顶板上任一点的坐标，尤其是z坐标数据。基于三维激光探测数据，借助于SURPAC软件，进行数据的获取，流程如图5所示。

图5  顶板z坐标数据提取流程

Fig. 5  Data extraction process of roof in z direction

利用得到的每个非空正方形网格角点对应的z坐标数据，基于MATLAB编制计算程序，计算任意尺度立方体覆盖整个顶板所需要的数量N_n(r)，进而求解三维分维数D。

利用上述算法，对石人沟铁矿斜井采区39号采空区顶板三维分型特性进行了研究。

首先得到不同尺度网格覆盖顶板xOy平面投影的非空网格分布及数量。根据表中数据，初始覆盖网格尺度为10 m，以1/2的倍率分别进行缩小并重新覆盖，得到对应的非空网格分布及数量。

表1所示为2个采空区顶板被不同尺度立方体完全覆盖时计算得到的数量。两者的双对数散点拟合曲线如图6所示。

斜井采区39号采空区顶板在xOy平面投影的计盒维数为1.845 2，而三维空间中分维值则为2.264 0，两者相差0.418 8。

由于考虑了z方向上维度的影响，顶板三维分维值出现了增大，斜井采区39号采空区顶板高差大，且起伏程度剧烈，三维分维值增加了0.418 8。顶板三维分维值可以较全面地反映采空区顶板的边界复杂程度、平整起伏程度及暴露面积大小，定量地表征了顶板的非线性程度。xOy平面分维值与三维分维值之间的差值代表了z方向(即高度)的维度，是一个小于1的分维数，与顶板的起伏程度有密切关系。

表1  立方体尺度与对应立方体数量列表

Table 1  Cube scale and corresponding cube number list  个

图6  顶板lnN_n(r)~lnr线性拟合图

Fig.6 Linear fitting diagram of lnN_n(r)-lnr of roof

5  采空区顶板边界分形插值模拟

5.1  采空区顶板二维剖面分形插值数学模型

由于采空区，通过试验或监测往往只能得到部分数据，为了更加精确地得到事物内在的本质规律，需要通过插值来进行补充。插值实际上是一种近似的过程，因此，插值函数的精度决定了最终结果的准确程度^[17]。对于高度非线性数据的差值计算，分形插值理论具有明显的优势，通过特殊的迭代函数IFS实现最终的迭代^[18-19]。

由前述章节分析可知，采空区具有明显的分形特性，边界线高度非线性，采用分形插值得到的结果更加准确。

一般的构造迭代函数IFS如下：

          (13)

式中：ω_i为变换矩阵，有

    (14)

式中：M为伸缩变换矩阵；N为平移变换矩阵。

设给定集合R，2个端点坐标分别为(x₀，y₀)，(x_L，y_L)，其中2个相邻的插值点为

           (15)

式中：i=1，2，…，P；P≤L。则有

，     (16)

一般地，设d为常数，定义为变换矩阵ω_i纵向压缩因子，取0≤d＜1。纵向压缩因子的大小决定了分形插值曲线模拟的准确程度，d越大，准确程度越大。

对于单因变量函数曲线，x_i不随y_i变化而改变，因此，b=0，进而求得其余系数：

     (17)

对于上述IFS迭代函数，有且只有一个吸引子，插值曲线会通过所有插值点并不断趋向于某个连续函数图形^[20]。

将已知点的数据代入迭代函数中，进行运算，直到收敛至吸引子。

采空区边界在平面上表现为封闭的曲线，可采用参数方程进行表征，因此，进行采空区边界分形插值模拟实际上就是对参数方程的迭代运算过程。

若已知点集，其中，且，设向量函数为，且。

若上式满足，即

               (18)

其中：i=1，2，…，n。

式(18)为采空区边界的参数方程，且可进行分形插值，其中，为插值区域。记，则。

根据分形插值迭代函数IFS构建原理，存在映射变换矩阵ω_i，使

(19)

式中：a_i为ω_i中t的压缩因子，且0＜a_i＜1。设，则称P_i为ω_i中x和y的压缩矩阵，且其范数满足0＜||P_i||₂＜1，上述映射满足

，       (20)

由式(20)可得

        (21)

经过计算可求得

(22)

将边界上已得到的点的坐标数据代入式(21)进行迭代运算，直到最终收敛至吸引子，完成最终的分形插值过程。

要进行采空区边界线的迭代计算，还需满足以下条件：

1) 对于边界线上的x和y两者互相独立，因此，d_i=0，g_i=0。

2) c_i与h_i分别表示t对x和y的影响程度，可设h_i=k_ic_i，k_i由下式求得：

                 (23)

为保证为压缩映射，应使0＜||P_i||₂＜1，即

           (24)

3) a_i为参数t的压缩因子，其取法与普通函数分形插值相同。

根据参数形式分形插值IFS迭代原理及上述条件，将已知点代入进行迭代求解，完成对边界线的插值运算，实现对边界线的模拟。

基于分形插值理论的采空区边界模拟及预测步骤如下：

1) 为保证数据的唯一性，首先要进行坐标变换和归一化处理，变换后坐标原点(x₀，y₀)为

，      (25)

则边界上的数据点(x_i，y_i)归一化后为

，        (26)

2) 确定t的压缩因子，为c_i设定初始值，并求取k_i。

3) 求取区间内的仿射变换迭代函数系统IFS，并得到每个参数的数值。

4) 进行迭代计算，实现采空区边界模拟及预测。

5.2  实例分析

对于部分人员、设备无法进入的采空区，物探和钻探是主要的探测手段，但由于场地或条件的限制，物探和钻探只能通过定点的方式获取少量的数据点，探测精度较低，因此，需要根据已获得的点对未知的点进行插值模拟。

进行物探之后，往往需要进行钻探定位，钻探定位时通常将探测区域划分为若干剖面，每个剖面进行选择性的钻探定位。

以石人沟铁矿为实例，在探测过程中，将某采空区区域划分为17个剖面，每个剖面选择7~12个位置进行定位。图7所示为某一剖面得到的8个点，坐标分别为：

图7  某剖面现场钻探点示意图

Fig. 7  Schematic diagram of site drill of some cross-sections

采用本文的分形插值方法，设c_i=0.2，代入已知点坐标，进行迭代运算，完成插值运算，最终的差值结果如图8所示，图8中粗线为分形插值后得到的顶板边界线。

图8  分形插值后得到的顶板线

Fig. 8  Roof line obtained by fractal interpolation

由图8可以看出：分形插值得到的采空区边界线具有一定的波动性，分形插值曲线的波动程度与已知点所在位置和c_i有密切关系。

采用同样的方法对其他剖面进行插值运算，得到最终的顶板曲面如图9所示。

由图9可知：通过分形插值运算之后，对顶板进行了有效的模拟，大大提高了采空区的探测精度。

图9  分形插值后的顶板三维曲面

Fig. 9  Three-dimensional surface of roof by fractal interpolation

6  结论

1) 采空区分形特性研究表明，采空区三维顶板都具有明显的自相似性和无标度特性，可采用分形几何理论体系进行研究。

2) 分维值采空区的非线性定量表征，具有很好的物理意义。任意剖面边界线分维值表征了边界线的复杂程度及变化规律，剖面二维平面区域分维值表征了区域复杂程度，而三维顶板分维值则表征了采空区顶板空间形态的复杂程度。三维顶板分维值与顶板xOy平面投影分维值插值为小于1的分数，表示高度方向起伏的复杂程度。

3) 根据分形插值理论，可运用少量数据，插值得到顶板的分形曲面，并使插值得到的顶板曲面通过已知的插值数值点，这对复杂的采空区顶板的模拟研究和直观显示具有重要的意义。

4) 采空区分形特性的研究，具有重要的理论和实际意义，可从本质上揭示采空区的非线性规律及特征，并可进行进一步理论研究，应用前景广阔。

参考文献：

[1] SAMUEL A L, JRGEN F B, GREGORY E B, et al. Computational fluid dynamics simulation on the longwall gob breathing [J] International Journal of Mining Science and Technology, 2017, 27(2): 185-189.

[2] ZHU Defu, TU Shihao. Mechanisms of support failure induced by repeated mining under gobs created by two-seam room mining and prevention measures[J]. Engineering Failure Analysis, 2017, 82: 161-178.

[3] TANG Mingyun, JIANG Bingyou, ZHANG Ruiqing, et al. Numerical analysis on the influence of gas extraction on air leakage in the gob[J]. Journal of Natural Gas Science & Engineering, 2016, 33: 278-286.

[4] SAKAMOTO A, YAMADA N, IWAKI K, et al. Applicability of recycling materials to cavity filling materials[J]. Journal of the Society of Materials Science, 2005, 54(11): 1123-1128.

[5] WHITTLES D N, LOWNDES I S, KINGQMAN S W, et al. The stability of methane capture boreholes around a long wall coal panel[J]. International Journal of Coal Geology, 2007, 71(2/3): 313-328.

[6] SU Hetao, ZHOU Fubao, SONG Xiaolin, et al. Risk analysis of coal self-ignition in longwall gob: a modeling study on three-dimensional hazard zones[J]. Fire Safety Journal, 2016, 83: 54-65.

[7] 来兴平, 张立杰, 蔡美峰. 神经网络在大尺度采空区细观损伤统计与预测中应用[J]. 北京科技大学学报, 2003, 25(4): 300-303.

LAI Xingping, ZHANG Lijie, CAI Meifeng. Application of neural network to the statistics and prediction of dynamical damage and evolvement in the large scale mine-out area supported by rock-based composite materials[J]. Journal of University of Science and Technology Beijing, 2003, 25(4): 300-303.

[8] 刘科伟, 李夕兵, 刘希灵, 等. 复杂空区群露天开采境界三维可视化及其应用[J]. 中南大学学报(自然科学版), 2011, 42(10): 3118-3124.

LIU Kewei, LI Xibing, LIU Xiling. 3D visualization of complicated cavity group under open-pit limit and its application[J]. Journal of Central South University (Science and Technology), 2011, 42(10): 3118-3124.

[9] 李夕兵, 李地元, 赵国彦, 等. 金属矿地下采空区探测、处理与安全评判[J]. 采矿与安全工程学报, 2006, 23(1): 24-29.

LI Xibing, LI Diyuan, ZHAO Guoyan, et al. Detecting disposal and safety evaluation of the underground gob in metal mines[J]. Journal of Mining & Safety Engineering, 2006, 23(1): 24-29.

[10] 宋卫东, 付建新, 杜建华, 等. 基于精密探测的金属矿山采空区群稳定性分析[J]. 岩土力学, 2012, 33(12): 3781-3787.

SONG Weidong, FU Jianxin, DU Jianhua, et al. Analysis of stability of gob group in metal mines based on precision detection[J]. Rock and Soil Mechanics, 2012, 33(12): 3781-3787.

[11] 熊立新. 采空区信息网络管理系统研究[D]. 长沙: 中南大学资源与安全工程学院, 2010: 100-105.

XIONG Lixin. The research on network management system of gob information[D]. Changsha: Central South University. School of Resources and Safety Engineering, 2010: 100-105.

[12] BERRY M V, LEWIS Z V. On the Weierstrass-Mandelbrol fractal function[J]. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1980, 370(1743): 459-484.

[13] 王国艳, 于广明, 于永江, 等. 采动岩体裂隙分维演化规律分析[J]. 采矿与安全工程学报, 2012, 29(6): 859-863.

WANG Guoyan, YU Guangming, YU Yongjiang, et al. Study on cracks fractal evolution laws of mining rock mass[J]. Journal of Mining & Safety Engineering, 2012, 29(6): 859-863.

[14] 姜仁义. 矿体的相似原理与标准采矿方法设计[J]. 现代矿业, 2014, 30(1): 17-21.

JIANG Renyi. Similarity principle of the orebody and the standard mining method designing[J]. Modern Mining, 2014, 30(1): 17-21.

[15] 杨成杰, 吴冲龙, 张夏林. 基于实体与块体混合模型的三维矿体可视化建模技术[J]. 煤炭学报, 2012, 37(4): 553-558.

YANG Chengjie, WU Chonglong, ZHANG Xialin. Modeling technology of 3D visualization orebody based on solid and block model[J]. Journal of China Coal Society, 2012, 37(4): 553-558.

[16] 张哲儒, 毛华海. 分形理论与成矿作用[J]. 地学前沿, 2000, 7(1): 195-204.

ZHANG Zheru, MAO Huahai. Fractal theory and mineralization[J]. Earth Science Frontier, 2000, 7(1): 195-204.

[17] 王威, 苏经宇, 马东辉, 等. 岩体质量评价的分形插值模型[J]. 中南大学学报(自然科学版), 2012, 43(12): 4827-4833.

WANG Wei, SU Jingyu, MA Donghui, et al. Evaluation of rock mass quality based on fractal interpolation model[J]. Journal of Central South University (Science and Technology), 2012, 43(12): 4827-4833.

[18] WANG Yidong, DAN Wenjiao, XU Yongfu, et al. Fractal and morphological characteristics of single marble particle crushing in uniaxial compression tests[J]. Advances in Materials Science and Engineering, 2015, 2015(1): 537692-1-9.

[19] 桂蕾, 殷坤龙. 基于多重分形理论的滑坡地表监测位移分析[J]. 中南大学学报(自然科学版), 2014, 45(11): 3908-3914.

GUI Lei, YIN Kunlong. Analysis of landslide surface monitoring displacement based on multi-fractal theory[J]. Journal of Central South University (Science and Technology), 2014, 45(11): 3908-3914.

[20] FENG Zhigang, PAN Xuezai, DAI Guoxing, et al. Research on rock fracture surface morphology characterization under brazilian test[J]. Abstract and Applied Analysis, 2014, 2014(2): 434898-1-10.

(编辑  刘锦伟)

收稿日期：2017-02-17；修回日期：2017-04-11

基金项目(Foundation item)：国家自然科学基金资助项目(51274023)；中央高校基本科研业务费专项资金资助项目(FRF-TD-17-025A2) (Project(51274023) supported by the National Natural Science Foundation of China; Project(FRF-TD-17-025A2) supported by the Fundamental Research Funds for the Central Universities)

通信作者：宋卫东，博士，教授，从事矿山开采安全控制及岩石力学研究；E-mail：songwd@ustb.edu.cn