张小艳, 朱圣凯, 杨鑫磊

(西安科技大学计算机科学与技术学院，西安 710000)

传统的二维图表在地质构造规律及拓扑关系表达上受限于其平面化表达方式，存在着空间信息缺失、制图过程复杂的缺陷，越来越无法满足人们对地层认识和空间分析的需求。而三维模型由于其直观的空间关系、可视化的操作分析等优点，在表达矿山地质体的地层信息及空间构造特征方面具有巨大优势[1-4]。

近年来随着计算机技术的发展，以及数字矿山、智慧矿山等概念的提出，矿山地质体三维建模受到越来越多的关注[5]。中外学者对矿山地质三维建模的进行了大量的研究，吴志春等[6]利用PRB数字填图技术在GOCAD软件平台上构建浅部三维地质模型，所建立的地质界面与实际地质数据较为吻合，但所建三维模型难以在其他软件平台上展示，移植性较差；Yucel等[7]借助于Agisoft软件进行无人机图像的3D建模，实现了Etili Comakli露天煤矿三维地形模型的构建与可视化，高效地完成了三维地质建模过程中涉及的图像匹配、地理配准和数字高程建模过程；坚润堂等[8]基于DIMINE软件，采用距离幂次反比法对矿体进行品位推估并建立西藏白容—岗讲矿区矿床三维实体模型；魏永勇[9]采用C/C++及IDL混合编程，借助于二次曲面Shepard插值法构建了地层的数字高程模型，实现了三维地层建模，但所建模型中地层交界线缺乏必要的曲线平滑，预测精度需进一步提高；此外杨彪等[10]借助Surpac矿业软件结合VRML和HTML语言, 构建矿山地表、地层、断层等地质模型，实现了矿山在Web环境下三维可视化与交互。

综上，已有矿山地质三维建模的研究虽各有特色，但也有其缺陷和不足，借助ARC/INFO、3DMine等C/S结构建模软件实现的矿体三维模型难以跨平台使用，无法与矿山企业现有的Web管理系统相结合；而利用VRML、Flash等Web 3D[11]技术实现的矿体三维模型则需要额外加载浏览器插件或者虚拟机才能完成矿山地质体模型的渲染，限制了矿体三维模型在Web端的可视化；此外建模过程中采用单一空间插值预测算法难以处理复杂的地质构造，无法实现切合实际地质情况的三维地质模型。基于以上分析，利用一种可在变异过程中修正变异方向的自适应差分进化算法来优化Kriging空间插值算法(updating and modifying adaptive differential evolution Kriging algorithm, UMDE-Kriging)，结合规则格网法建立精确度更高的数字高程模型(digital elevation model, DEM)，并在此基础上研究基于Web GL标准、以浏览器为载体的三维数字化建模方式，利用Three.js和Catmull-Rom样条曲线实现无插件、可移植、跨平台、支持多浏览器运行的采煤工作面煤层三维可视化模型，为煤矿采煤工作面的开采掘进提供重要参考。

1 关键技术

1.1 Web GL

Web GL(Web graphics library)是一种跨平台的3D绘图规范，通过在OpenGL ES 2.0中添加JavaScript绑定，从而为HTML5中Canvas提供硬件3D加速渲染[12]。与早期的Web 3D实现相比，Web GL能够在大多数平台上运行，兼容性好，同时还支持CPU硬件加速，实现三维图形渲染。

1.2 Three.js介绍

Three.js依据Web GL规范，对底层Web GL代码进行简单封装，可以将二维数据以三维方式呈现出来，是一款开源的Web GL图形开发框架。简单、直观的封装Web GL常用的三维对象，通过掩饰细节来加快开发效率，减轻开发者的开发负担[13]。由于Three.js采用了JavaScript语言，因此能够很容易的与其他浏览器组件进行交互。

Three.js与其他的Web GL图形框架一样，在使用时都需先定义三维世界中所包含的基本元素，主要流程如下：

(1)设置三维场景(scene)：场景是用来存放一切所创建对象的容器，用THREE.Scene来表示。

(2)相机设置(camera)：相机决定了显示场景中哪一个角度的景色。

(3)设置灯光源(light)：Three.js中提供一系列的光源，包括环境光、方向光、点光源等。

(4)设置物体(object)：物体就是场景中需要处理的三维对象。Three.js通过加载3D建模工具导出的模型文件或者借助THREE.mesh三维网格动态创建三维场景中的模型。

(5)添加渲染器(renderer)：渲染器决定了3D空间中的物体对象映射到二维平面上的方式，以及渲染的结果应该绘制到场景中什么元素上。

2 自适应差分进化Kriging算法

2.1 Kriging算法

Kriging插值法又称空间自协方差最佳插值法，是依据协方差函数对随机过程/随机场进行空间建模和预测(插值)的回归算法[14]。根据采样点与预测点之间的空间位置和相关程度，确定一个影响待插点值的距离范围。对其中每一个样品赋予一定的权重，最后进行加权平均来估计待预测点属性值。

Kriging插值法的一般公式为

(1)

式(1)中：Z*(x)为待求插值点的估计值；Z(xi)为第i个采样点xi处的数值；λi是第i个采样点的插值权重系数，并且需满足无偏性和最优性两个条件。即：

(2)

σ2=E[Z*(x)-Z(x)]2=min

(3)

式中：σ为待求插值点的估计值Z*(x)和真实值Z(x)之间的标准差，在Kriging法中，用来衡量各个样本点之间空间相关程度的是变差函数。

(4)

式(4)中：h为两点之间距离；N是由h分开的成对样本点的数量；Z(xi)是点的属性值。当满足本征假设时，可用变差函数来表示Kriging方程组：

(5)

式(5)中：γ(xi，xj)表示两个采样点xi和xj的变差函数值；μ为拉格朗日系数；j=1，2,…，n。

2.2 差分进化算法

差分进化算法(DE)思想源于早期提出的遗传算法(GA)，模拟遗传学中的杂交、变异、复制来设计遗传算子，通过优化具体参数、值函数解决连续的非线性优化问题[15]。

使用DE算法中的变异、交叉和选择运算实现种群的进化，其中差分变异公式为

vi(g+1)=xr1(g)+F·[xr2(g)-xr3(g)]，

i=1,2,…,NP, 1≤r1≠r2≠r3≤NP

(6)

式(6)中：vi(g+1)为g+1代变异个体；xr1(g)、xr2(g)、xr3(g)是从t代种群中随机选取的三个个体；F为缩放比例因子；NP表示种群大小。

由g+1代变异个体vi(g+1)与上一代个体xi(g)交叉得到g+1代子个体ui(g+1)，交叉公式为

(7)

式(7)中，当rand(0,1)函数产生的随机数小于交叉概率RC时，ui(g+1)中对应参数来自变异个体vi(g+1)，否则该参数来自上一代个体xi(g)。

DE算法采用贪婪算法来选择进入下一代种群的个体，选择公式为

(8)

式(8)中，f为待优化求解的目标函数。

2.3 UMDE-Kriging算法实现

基于差分进化算法以及Kriging插值算法，结合煤矿采煤工作面实际地质情况以及矿井测量人员经验公式，设计合适的自适应放缩因子F，通过差分、变异、选择建立变差函数模型参数最优解，进而求解出满足精度的待预测点高程属性值。

UMDE-Kriging算法主要步骤如下：

2.3.1 确定变差函数参数

通过式(4)得出实验变差函数云图，选定合适的变差函数模型，并确定其参数c0、c、a，其中c0为块金值，c为偏基台值，a为变程。

2.3.2 选定合适目标函数

变差函数模型的拟合本质上是表征实验变差函数与理论变差函数的近似程度。因此，可以将各距离下的实验变差函数值与理论变差函数值之间的差值平方和作为目标函数。

(9)

式(9)中，γ*(h)为理论变差函数值；γ(h)为实验变差函数值；m为实验变差函数值的个数。

2.3.3 初始化种群

在解空间随机生成M个个体，因为解空间是模型参数c0、c、a的集合域，即每个个体由三维向量组成。

(10)

每个个体中的参数取值如下：

(11)

式(11)中，(C0_max-C0_min)、(Cmax-Cmin)、(Amax-Amin)分别为模型参数c0、c、a的取值区间，需要人为设定。

2.3.4 确定差分进化算法的有关运行参数

在差分进化算法的变异过程中通过增加自适应缩放因子F，修正变异方向。确定最大进化代数、种群大小和交叉概率。

(12)

式(12)中，F0为设置的缩放因子；fi(pre)为当前个体的上一代个体适应值；fmin(pre)为上一代最优个体适应值；fi(cur)为当代个体的适应值。

2.3.5 差分进化运算

在随机生成的初始种群基础上，进行变异[式(6)]、交叉[式(7)]和选择[式(8)]三个过程，多轮进化迭代直至得出最优解fmin，求得变差函数模型参数，拟合出相应的模型表达式。

2.3.6 预测值计算

将得到的变差函数模型表达式代入式(5)中求解，得出采样点权重系数λi，进而通过式(1)得到待预测点属性值。

3 数据准备

选取某煤矿16402采煤工作面0～700 m深度范围内的54组巷道钻孔数据以及切眼地质采样数据作为样本数据集合(包括22组巷道钻孔数据和32组切眼采样数据)，如表1所示。

表1 采样点数据

所选工作面区域横向宽度为200 m，钻孔采样相对煤层底板的高程范围为0～2 m，样本点分布如图1所示。

图1 采样点分布

3.1 算法可行性验证

从样本数据集合中提取出同属于分一煤层下表面的控制点数据，选取Kriging插值法中的指数模型作为优化对象，利用修正变异方向的自适应差分进化算法(UMDE)进行优化。其中，差分进化算法的参数设置如表2所示，模型参数边界区间如表3所示。

表2 差分进化算法的参数设置

表3 模型参数的边界区间

首先根据54组初始样本点数据按照采样编号绘制出高程值分布折线A；之后分别根据UMDE-Kriging算法和普通Kriging算法计算出相应位置的高程理论值，绘制出相应折线B和折线C，如图2所示。

图2 采样点高程分布折线图

根据54组预测结果和实际高程值进行交叉验证，分别计算均方根误差(RMSE)和平均绝对百分误差(MAPE)。

(13)

(14)

表4 算法预测精度对比

由图2和表4可知，UMDE-Kriging算法插值的结果要比普通Kriging算法更加接近于与实际采样点高程值，误差更小，在高程预测求解过程中表现出更好的有效性和准确性。因此，将UMDE-Kriging插值算法应用于采煤工作面信息预测完全可行。

3.2 数据格网化处理

数字高程模型(DEM)是通过有限的地形高程数据实现地形表面的数字模拟或表面地形的数字表示[16]，其中规则格网DEM以其数据结构简单、运算速度快、效率高的优点大量应用于实际工程。

根据采煤工作面边界约束，利用规则格网法生成各地质分层的格网表面，并结合UMDE-Kriging空间插值算法进行任意格网点的属性值预测和采集，具体步骤如下：

3.2.1 地层划分

针对16402采煤工作面涉及的煤层地质数据进行综合分析，将相同的岩性并且垂直方向位置相似的看作同一层，根据地层分界控制点信息将其划分为分一煤层、夹矸层、分二煤层。

3.2.2 格网划分

分别计算采煤工作面各地层起始深度和终止深度的采样点水平坐标(x，y)的最大值和最小值，合理划分网格单元大小，计算每个格网点的x、y坐标。

3.2.3 格网点插值

利用UMDE-Kriging空间插值算法对各地层格网点的高程、煤质指标进行插值拟合处理，生成具有属性值的规则格网模型。

将16402采煤工作面0～700 m范围区域划分为11×11格网，以分一煤下表面为例，结合UMDE-Kriging算法建立规则格网模型，如图3所示。

图3 规则格网模型

4 煤层三维模型实现

煤层是空间中的一种层状地质实体，其空间形态与分布相对于非层状实体较为简单，具有较强的全局规律性[17]，针对采煤工作面的煤层地质情况，可将其分为开采煤分层和夹矸层。

以16402采煤工作面54组样本数据为基础，采用JetBrains WebStorm 2018作为开发工具，利用JavaScript与Three.js进行开发。借助于规则格网法和UMDE-Kriging插值法形成规则网格数据模型，根据点聚成线，线聚成面，面聚成体的原则逐步实现采煤分层模型与夹矸层模型；建立采煤工作面0～700 m范围内煤层三维体元模型，实现部分煤质信息的动态展示。其中，所建模型长度为700 m，宽度为200 m，高度为2 m。

4.1 采样点分布模型

将各煤层规则网格数据模型中每个数据点视为一个独立的球体模型，依据其三维坐标属性值(x，y，z)确定其空间位置，结合THREE.SphereGeometry()构造方法循环完成各煤分层和夹矸层所有数据点的绘制。其中点半径选取为30，球体维度方向和经度方向上的分段数均设置为10：var geometry=new THREE.SphereGeometry(30, 10, 10)，采样点分布模型如图4(a)所示。

4.2 煤层面平滑曲线模型

Catmull-Rom样条曲线是一种3次的分段插值样条，它利用前后的2个节点直接计算出各点处的切矢量，每相邻4个节点的信息就可以构造出一段经过中间两点的3次多项式函数，并且所构造的每一段样条能够保持连续[18]。Catmull-Rom曲线每个分段拥有四个控制点，假定这四个控制点分别为P0、P1、P2、P3。则Catmull-Rom样条函数的矩阵表达式为

(15)

式(15)中，浮点坐标t∈[0,1]，其中M为样条曲线插值矩阵：

(16)

式(16)中，T称为形状因子，当样条用于离散点连续平滑时，取T=0.5为最佳值。根据式(15)和式(16)可以在P1、P2之间构造出一段平滑的曲线，经过拼接每一段平滑曲线可构成经过所有控制点的完整连续曲线。

利用所建立的规则格网模型，分别存储其中同属一条线上的控制点信息，借助Catmull-Rom插值算法绘制Catmull-Rom样条曲线，建立煤层面光滑曲线模型，为煤层面曲面模型绘制打下基础。如图4(b)所示，以分一煤层下表面格网模型数据为基础，分别绘制巷道平行方向和垂直方向Catmull-Rom样条曲线。

4.3 多煤层面元模型

分层建模的基本思想是按照所划分的地层顺序，根据不同地质层面的控制点依次建立各地层面模型。基于所绘制的Catmull-Rom样条曲线，反算每一条曲线上的点的空间位置信息。依据点在曲线的相对位置，按照弧长比例返回其对应的三维坐标值，并将其作为所属地质分层曲面控制点；之后借助Three.js自带的参数化几何模型THREE.ParametricGeometry依次实现地层面曲面模型的绘制：var geometry=new THREE.ParametricGeometry(function(u,v), slices, stacks)，其中slices用于参数化函数的切片数量u；stacks用于参数化函数的堆栈数量v，循环调用function()返回所有曲面控制点的三维坐标信息。利用接收到控制点位置信息并绘制各煤层曲面模型，采煤工作面各地质分层(煤层与矸石层)曲面模型如图4(c)所示。

4.4 体模型绘制

通过获取采样点水平方向的空间分布信息和垂直方向的煤层分层信息，在生成多个描述不同煤层地质面模型的基础上，利用各煤层上下表面边界线，在这两条闭合轮廓线之间建立一系列的三角面片，设置当前层的颜色，构建出煤层侧面模型；结合各煤层顶板面模型、底板面模型共同组合成封闭立体空间，完成采煤工作面各地质分层模型绘制。如图4(d)所示，每个地质分层(煤层与矸石层)都是有上下两个表面模型和若个侧面模型“缝合”而成。

图4 采煤工作面模型

4.5 三维模型交互显示

4.5.1 人机交互

基于Three.js建立的工作面煤层三维模型与用户的交互采用THREE.OrbitControls控件实现：controls_camera=new THREE.OrbitControls(camera, renderer.domElement)。

在这种控制模式下，分别利用鼠标左键、中键、右键改变场景中相机的视点，实现三维模型的旋转、放缩、平移。从不同的角度观察采煤工作面煤层分布情况，如图5所示。

图5 人机交互展示

4.5.2 煤质信息展示

通过给采煤工作面三维模型绑定鼠标响应事件来完成各煤层煤质信息的展示：documenet.addEventListener(“mousemove”,函数,false)。

当鼠标光标移动至三维模型中某个平面时，利用THREE.Raycaster控件完成模型中平面对象的拾取，实现选中煤分层或矸石分层表面采样点三维坐标和煤质信息展示。如图6所示，在左侧信息显示区域展示选中曲面所有控制点的三维坐标以及部分煤质指标。

Td为真密度，Ad为灰分

5 结论

以某矿16402综采工作面的部分采样数据为基础，验证了UMDE-Kriging空间插值算法可行性，并结合规则格网法完成虚拟采样点高程值及煤质指标的插值拟合处理，建立了精确度更高的数字高程模型，并借助Three.js技术和Catmull-Rom样条曲线构建了能够反映实际地质特征的16402采煤工作面煤层三维可视化模型。

(1)所建模型克服了基于其他建模技术需要额外加载插件才能完成地质体模型渲染的缺陷，具有可移植、跨平台、支持多浏览器运行的优势。

(2)实现了三维模型的旋转、平移、放缩和分层显示，具有良好的人机交互性能。

(3)对工作面的地质空间特性及煤质指标进行了良好的展示，在表面模拟和内部属性表达两个方面达到较好的融合。

(4)能够完整、便捷的嵌入至煤炭企业Web端信息管理系统中，为矿山开采掘进提供技术支持。

但本模型的构建针对煤层中可能存在的断层、褶皱等特殊地质区域还不能很好地呈现，在后续研究中，将结合煤田地质学以及空间预测算法进行修正，进一步完善。