贾庆仁, 车德福, 李佳徐, 陈凯
东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2017-05-16
基金项目：国家自然科学基金资助项目(51474050，U1602232)；中央高校基本科研业务费专项资金资助项目(N160106002)。
作者简介：贾庆仁(1990－)，男，山东日照人，东北大学博士研究生;
车德福(1970-)，男，山东海阳人，东北大学教授。

摘要：现有的煤层建模方法大都从多源数据中提取点或边界来构建实体模型, 对开采过程中煤层动态变化的过程缺乏足够关注.提出一种逐步细化的煤层精细建模方法, 首先, 对开采边界内进行网格划分后将其转换为三角网, 通过全局数据中提取的采样点和普通克里金法对三角网节点的高程及厚度属性进行插值, 之后将其拓展为三棱柱作为初始煤层模型; 其次, 确认开采数据的影响范围, 并对初始模型进行动态修正.在唐山市钱家营煤矿进行应用的结果表明, 该方法充分运用了开采数据, 使煤层的真实变化能够及时反映在煤层三维模型上, 提高模型的精度.
关键词：多源数据煤层建模普通克里金法断层模拟变异函数
Three-Dimensional Modeling Method of Coal Seam with Gradual Refinement
JIA Qing-ren, CHE De-fu, LI Jia-xu, CHEN Kai
School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: CHE De-fu, professor, E-mail: chedefu@mail.neu.edu.cn
Abstract: The existing coal seam modeling methods using discrete points or geological boundaries from multi-source data lack enough attention to the dynamic changes of coal seam in the mining process. A new 3D coal seam fine modeling method is proposed. Firstly, creating regional grid of the mining area and converted to triangulated irregular network (TIN), inter polating attributes of TIN nodes (elevation and thickness)using ordinary Kriging method, and then a generalized tri-prism(GTP) model of coal seam will be built as the initial model. Secondly, the initial model is refined gradually by the dynamic data within the influence domain of these data. A prototype system is created and implemented in Qianjiaying coal mine, Tangshan, China. The result shows that this method can accurately reflect the change of 3D coal seam model from the dynamic data, thus the accuracy of the 3D coal seam model can be improved gradually.
Key Words: multi-source datacoal seam modelingordinary krigingfault simulationvariogram
传统的图纸在地下地质构造发育、地下拓扑关系表达上受限于其平面化表达方式, 往往相互割裂甚至冲突.而三维模型由于其直观的空间关系、可视化的操作分析及数值模拟等优点, 可显著提升现有数据的表达和价值.这使得整合矿山生产不同阶段的各类数据、建立三维地质模型、表达乃至再现地下地层及构造发育情况成为一项重要任务^[1-3].
近年来随着三维技术的发展和数字矿山概念的提出, 国内外对三维地质建模及应用进行了大量研究, 对建模所需要的三维数据结构进行了界定, 将其划分为面元模型、体元模型和混合模型^[4].提出了包括基于序列剖面自动、半自动连接的地质体建模方法^[5]; 基于钻孔数据构建的广义三棱柱建模方法等.上述研究充分运用区域数据, 部分解决了稀疏数据下地质体及地质构造三维重建的问题并得到广泛应用.
与地质体建模不同, 煤层体三维建模的空间形态及关系相对简单, 但对精度要求更高, 因为精度的提升对煤层开采环境构建、井下工程建设、储量计算等具有重要意义.相较于地质体建模, 专门针对煤层的地质建模研究较少, 且大多沿用三维地质建模的方法.如基于地统计学、结合钻孔及地震剖面数据的煤层建模方法^[6]; 针对边界约束及钻孔数据的线框、块段建模的方法^[7]等.上述研究取得了一定的成果, 但忽略了对于局部具有重要意义的动态开采数据^[8].文献[9]尝试加入开采数据, 使用数学曲面造型技术在煤层模型的基础上构建细部曲面, 并进行拼接操作, 但由于与初始模型的构造方法不同, 割裂了区域整体数据和局部数据之间的联系, 未能实现方法上的统一.
针对上述问题, 本文提出了一种逐步细化的煤层精细建模方法, 该方法加入开采数据, 将建模过程分为两步:第一步, 通过OK(ordinary Kriging)计算建模区域划分的格网点的高程和厚度属性数据, 然后建立TIN网, 并最终拓展为GTP(generalized tri-prism)煤层实体模型^[10]; 第二步, 使用开采数据对已有模型局部进行动态更新, 提高模型精度.
1 动态精化的煤层三维建模方法煤矿所有生产活动都是围绕煤层展开, 由此产生的数据形式众多.本文将涉及的煤层数据根据全局规律性和局部异质性的特点分为两类.
1) 静态数据.静态数据包括在煤层开采之前及初期, 由地质勘探及巷道掘进过程中获取到的钻孔数据、DEM数据、实测数据.从静态数据中提取到的离散点, 可作为建立三维煤层模型的基础数据.但由于数据分布较为稀疏, 根据其建立的煤层模型在数据稀疏区域形态及空间分布上具有较大的不确定性, 对于小范围的分析计算难以适用, 可用作部分生产规划及计算.
2) 动态数据.动态数据是指在煤层工作面回采过程中获得的煤层相关数据, 包括石门见煤点、井巷数据、井下钻孔数据、工作面素描地质数据.由于动态数据是在采掘过程中产生的真实、局部的煤层空间分布及属性数据, 具有范围小、密度高的特点, 可以作为静态数据的补充对模型进行修正, 因此如何更好地利用动态数据决定了煤层模型精细程度.
根据两类数据特点提出的煤层精细模型的建立方法主要包含两部分, 即煤层初始建模(2节)和煤层模型动态精化(3节).技术路线图见图 1.
图 1(Fig. 1)

图 1 煤层精细建模方法技术路线图Fig.1 Fine modeling method of coal seam based on time series technology roadmap

2 基于克里金法的煤层初始建模2.1 克里金法基本原理相较于确定性方法(趋势面法、多项式法、数学曲面方法等), 克里金法对采样点进行插值的过程加入了空间结构化变异特征, 同时由于其最优无偏估计，能够给出计算误差等优点, 在土壤、大气等地学研究领域得到了广泛应用.其基本原理为：设区域化变量{Z(x)∈A}, 其中x表示空间位置.已知采样点x_i(i=1, 2, …, n)处的属性值为Z(x_i), 则根据OK原理, x₀处的属性值Z(x₀)的估计值是x_i属性值的加权和.

(1)

具体计算过程本文不再赘述.
2.2 煤层初始建模煤层是空间中的一种层状地质实体, 其概念空间形态与分布相对于非层状实体较为简单, 有较强的全局规律性, 在煤层地质体精细建模研究中, 由多源数据获取煤层表面面元模型, 再拓展到体元模型的建模方式能够在表面模拟效果和实体内部表达属性两个方面达到较好的融合, 而且维护模型所需要的数据量相对较小.
本文选择格网-TIN混合面元构建静态模型, 以获得更好的表面模拟效果和易于更新的几何结构.具体步骤如下：
1) 获取井田范围内所有煤层底空间离散点, 主要为钻孔点(图 2a);
图 2(Fig. 2)

图 2 煤层初始模型构建Fig.2 Construction of initial model (a)—空间散点分布；(b)—格网化并构建TIN网; (c)—高程、厚度插值；(d)—建立初始模型.

2) 根据井田边界约束生成规则格网(图 2b);
3) 通过OK法分别对格网点的高程、厚度属性进行插值(图 2c);
4) 由格网点根据Delaunay法则生成TIN网, 并根据各点厚度属性, 拓展连接煤层顶、底面对应三角形生成GTP模型, 最终构建起表面与内部拓扑一致的煤层实体模型(图 2d).
3 煤层模型动态精化更新时为保证格网变化的连续性, 首先需要确定数据的影响范围, 即选择哪些数据对格网重新插值.在插值过程中, 距估计点近的点会对远的点产生隔离作用, 可采用八分象限法等减少参加计算点的数量, 损失少量精度以换取更高的时间效率.获得某工作面一期新数据后, 将所有新数据加入点集P, {P(x_i)∈A, (i=1, 2, …, m)}, 针对格网更新范围的选取, 本文以象限法为基础, 加入距离因素, 提出了一种“八分象限距离加权”确定更新范围的方法.该方法分别计算P内采样点x_i对点集Q{Q(y_i)∈A, (i=1, 2, …, n)}(Q=A－P)内n个点的象限与距离加权和w₀至w_n后, 对结果进行排序, 取权值最大的8点加入数组R, {R(x)∈A}, 针对点x_i有

(2)

其中:d_j为x_i到y_j距离; d为x_i到Q的最小距离; n_j是对某一象限中根据距离排序Q在该象限的所有点后, 距离d_j所在的索引位置.
获取到R之后, 求取其凸包范围(图 3), 然后对该范围域内所有点的高程(厚度)值重新进行OK插值计算, 完成动态精化过程.
图 3(Fig. 3)

图 3 选取数据更新范围Fig.3 Select the data update range

4 断层处理本文采取“整体法”解决断层数据建模问题, 该方法不仅能解决多条断层相交难题, 也可以减少数据误差, 最大限度地保障煤层空间位置的准确性, 具体过程如下.
1) 计算断层面.断层由一个或几个断层点组成, 每点都包括走向、倾角、落差, 只记录一个断层点时还应该包含延伸数据.由断层初始点属性可计算断层三角面(图 4a).
图 4(Fig. 4)

图 4 煤层断层Fig.4 Coal sean fault (a)—二维图；(b)—三维渲染.

2) 计算断煤交线.断层面与煤层TIN网求交得到的交线即为断煤交线.由断层初始点落差及分配比率, 计算断煤交线上对应点落差.
3) 分割煤层TIN网.根据断层正、逆属性将地层分为断层上、下盘区域, 插入交线, 并将范围内的煤层三角形分割为上、下盘三角形.
4) 计算煤层错动.由断层面分别计算上、下盘三角形各角点于XOY平面错动距离和Z方向高程变化量.
5) 影响域计算.由断煤交线落差值计算断层影响域, 默认落差每米修正半径为10 m.计算在断层影响域内地层TIN网三角形点角点Z坐标变化量(图 4b).
6) 计算影响域内煤层错动.计算在正断层(逆断层忽略此步骤)影响域内地层TIN网角点XOY平面位置变化量.
5 实例研究通过VC++及OpenGL编程实现上述煤层精细建模方法, 并开发煤层三维建模系统.
5.1 基础煤层建模钱家营煤矿位于唐山市, 井田构造以褶曲为主、断裂为辅.井田内可采和局部可采煤层共8层, 即5, 6_-1/2, 7, 8, 9, 11和12_-1煤, 均属复杂结构的中厚-厚煤层, 井田范围88 km², 使用GeoMS3D对钱家营煤矿5, 7, 8, 9, 11, 12a煤层进行了三维建模.图 5a为整个井田边界内钻孔分布情况三维视图.钱家营煤矿共有勘探钻孔314个, 其中参与7煤层初始建模的钻孔有265个.
图 5(Fig. 5)

图 5 三维模型Fig.5 3D models (a)—钻孔分布情况；(b)—井田边界约束下的TIN网.

为建立钱家营煤矿煤层模型, 在井田边界范围内应用本文所述静态建模方法完成初始煤层模型的构建.由于井田范围过大, 因此建模过程中使用100 m×100 m的边长进行初始格网构建(图 5b), 对于需要进行加密的工作面或其他重点区域, 采用预先设定的20 m边长约束对格网进行加密.
使用R语言中的geoR开发包选择变异函数模型.首先建立5种变异函数模型(球状模型、指数模型、高斯模型、圆模型、三次曲线模型)对区域化变量的结构特征进行描述.采用加权最小二乘法来拟合参数, 通过交叉验证选取模型.
以265个钻孔资料为验证数据源, 以平均克里金方差(MKV)、平均误差(ME)、平均z分数(MZS)和平均累积概率(MCP)为统计学指标.其判断标准为:MKV和MZS尽可能小; ME尽可能接近于零; CP尽可能接近于1.表 1给出了5种模型统计数据的比较，图 6为5种模型的拟合曲线.
表 1(Table 1)

表 1 5种理论变异函数模型的精确度评价Table 1 Precision of 5 theoreticalvariogram 模型 MKV ME MZS MCP 高斯 53 953.04 -3.4×10-15 -1.1×10-16 0.503 261 指数 10 499.85 -2.163 53 -0.008 74 0.496 978 球状 6 854 -1.659 53 -0.008 27 0.497 146 圆形模型 6 317.06 -1.699 12 -0.009 25 0.496 754 三次曲线 567.38 -0.329 13 0.030 167 0.496 745

表 1 5种理论变异函数模型的精确度评价 Table 1 Precision of 5 theoreticalvariogram

图 6(Fig. 6)

图 6 5种理论变异函数模型的拟合曲线Fig.6 Fitting curves of 5 theoretical variational models

由表 1及图 6知, 球状模型、圆模型及三次曲线模型拟合效果较好.各模型MCP指标比较接近, 三次曲线模型MKV, ME指标最优且优势较大, 但是其MZS较大, 原因为MKV值小于其他任何模型的10%, 因此计算MZS时分母偏小.综合来看, 最终选择拟合效果最好的三次曲线模型作为变异函数.公式为
其中:块金值C₀=0；变程α=5 582.7；基台值s=61 760.3；h′为仅考虑平面各向异性, 经过各向异性变换后的各向同性距离：
其中, .
根据变异函数估计的井田范围内格网点的高程值, 建立的初始模型, 并按照同样的方法插值煤层厚度.
5.2 动态修正模型动态数据由现场人员进行整理并使用开发的功能模块录入数据库.如图 7a所示, 以7煤层2074工作面为例, 从2014年7月到2015年3月, 工作面范围内共新增点状数据(井下钻孔、石门见煤点)19个, 线状数据(素描煤层线)82条, 断层21条(其中可见断层3条), 使用这些数据对此范围内的煤层进行修正.其中断层数据由于来源不同, 其精确度也不尽相同.在开采过程中实测获得的、或由实测数据进行推断获得的断层精确度较高, 需要进行断层建模.
图 7(Fig. 7)

图 7 工作面修正过程Fig.7 Working face refining process (a)—工作面煤层相关数据(高程渲染)；(b)—修正前(厚度渲染)；(c)—修正后(厚度渲染).

图 7b、图 7c为2074工作面修正前、后的煤层模型, 按照时间顺序加入的动态数据改变了该区域的平均高程及厚度(表 2).在进行开采设计或工作统计的时候, 可以获得更好的效果.图 8为煤层及煤层相关工程建模效果.
表 2(Table 2)

表 2 煤层属性随时间推进变化情况Table 2 Change of coal seam properties with time m 变量 静态模型 (2014-05) 一期 (2014-07) 二期 (2014-09) 三期 (2014-12) 平均高程 -634.5 -627.1 -626.2 -623.9 平均厚度 3.5 3.6 3.9 4.1

表 2 煤层属性随时间推进变化情况 Table 2 Change of coal seam properties with time

图 8(Fig. 8)

图 8 所有煤层及相关工程建模效果图Fig.8 Coal seam and underground engineering models

6 结论1) 增加数据来源, 使用开采数据将煤层建模分为两个阶段, 并将其发展为持续过程, 通过增加可用数据提高了模型的精度.
2) 结合格网, TIN, GTP等二、三维数据结构, 建立了易于维护、可局部更新的煤层实体模型.
3) 使用“八分象限距离加权”方法确定格网更新范围, 确保了新数据与原始格网的连续性.
参考文献

[1]	Thore P, Shtuka A, Lecour M, et al. Structural uncertainties: determination, management, and applications[J].Geophysics, 2002, 67(3): 840–852.
[2]	Caumon G, Collon-Drouaillet P, Veslud C L C D, et al. Surface-based 3D modeling of geological structures[J].Mathematical Geosciences, 2009, 41(8): 927–945.
[3]	吴立新, 史文中, ChristopherG. 3D GIS与3D GMS中的空间构模技术[J].地理与地理信息科学, 2003, 19(1): 5–11. ( Wu Li-xin, Shi Wen-zhong, Christopher G. Spatial modeling technologies for 3D GIS and 3D GMS[J].Geography and Geo-Information Science, 2003, 19(1): 5–11.)
[4]	陈国良, 刘修国, 盛谦, 等. 一种基于交叉剖面的地质模型构建方法[J].岩土力学, 2011, 32(8): 2409–2415. ( Chen Guo-liang, Liu Xiu-guo, Sheng Qian, et al. A modeling method based on intersected geological sections[J].Rock and Soil Mechanics, 2011, 32(8): 2409–2415.)
[5]	Li X, Li P, Zhu H. Coal seam surface modeling and updating with multi-source data integration using Bayesian geostatistics[J].Engineering Geology, 2013, 164(4): 208–221.
[6]	Zhu L F, Li M J, Li C L, et al. Coupled modeling between geological structure fields and property parameter fields in 3D engineering geological space[J].Engineering Geology, 2013, 167(24): 105–116.
[7]	Tercan A, Eünver B, Hindistan M A, et al. Seam modeling and resource estimation in the coalfields of western Anatolia[J].International Journal of Coal Geology, 2013, 112(2): 94–106.
[8]	孙振明, 毛善君, 祁和刚, 等. 煤矿三维地质模型动态修正关键技术[J].煤炭学报, 2014, 39(5): 918–924. ( Sun Zhen-ming, Mao Shan-jun, Qi He-gang, et al. Dynamic correction of coal mine three-dimensional geological model[J].Journal of China Coal Society, 2014, 39(5): 918–924.)
[9]	车德福, 吴立新, 陈学习, 等. 基于GTP修正的R3DGM建模与可视化方法[J].煤炭学报, 2006, 31(5): 576–580. ( Che De-fu, Wu Li-xin, Chen Xue-xi, et al. Modeling and visualizing methods for real 3D geosciences model based on amended generalized tri-prism(GTP)[J].Journal of China Coal Society, 2006, 31(5): 576–580.)
[10]	Fasani G B, Bozzano F, Cardarelli E, et al. Underground cavity investigation within the city of Rome (Italy): a multi-disciplinary approach combining geological and geophysical data[J].Engineering Geology, 2013, 152(1): 109–121.