一种混合插值算法的露天矿产储量估算方法及其应用

王 璐,卢小平,李国利,李小雷,康爱峰

(1. 河南理工大学自然资源部矿山时空信息与生态修复重点实验室,河南 焦作454003;2. 鹤壁恒源矿业集团有限公司,河南 鹤壁458030)

传统二维GIS技术难以描述复杂环境地质体的三维空间信息,而三维GIS技术可以更好地解决这一问题[1-2]。近年来,随着智能化进程的加快和地下空间利用需求的激增,地质体三维建模成为三维GIS领域的研究热点。传统的地质体三维建模中,钻孔数据是最主要的数据来源[3]。但是钻孔数据往往成本较高、数量较少,用少量的钻孔数据表达复杂的地质体内部结构,往往需要对钻孔数据进行插值加密,生成虚拟钻孔,以提高建模的精度[4]。

地质体三维建模是矿产储量估算的基础,精确的地质体三维模型可以提高储量估算的精度,在矿区的生产规划和监测方面起到指导作用。地质体三维建模插值方法大体分为两类:一类是确定性插值方法,另一类是地统计插值方法。文献[5]利用局部多项式插值算法进行插值拟合,计算了矿山现开采范围内的土石方量值。文献[6—8]采用普通克里金和反距离加权插值法对矿体进行了储量估算。文献[9—11]利用克里金插值进行地质体三维建模,进而对矿产储量估算。文献[12]对比了克里金插值、样条插值及反距离加权插值,将矿体表面模型进行叠加,采用矿量计算公式,得到弓长岭独木矿区铁矿的总储量。但这些方法大都忽视了地质原理对插值方法的要求,对地质体采用单一的插值方法难以取得最优的插值结果。

针对以上问题,本文结合鹤壁市二道庄露天矿的实际情况,提出一种混合插值算法。一方面,对钻孔分层数据进行空间自相关分析,针对具有空间自相关性的地层采用克里金、规则样条、张力样条、反距离加权进行插值,对不具有空间相关性的地层,则利用规则样条、张力样条、反距离加权进行插值,通过对比分析选择出适用于每个地层的插值方法;另一方面,利用Dynamo进行可视化编程,建立基于钻孔数据的地质体三维模型,并基于地质体三维模型,利用SuperMap进行二次开发,建立矿产储量估算系统。

1 研究方法

本文研究方法的具体方法流程如图1所示。

图1 方法流程

1.1 插值算法

空间插值是将空间连续的点转化为连续的数据曲面。根据插值形式的不同,可分为空间外插和空间内插[13];根据插值方法的不同,可分为确定性方法和地质统计学方法[14]。确定性插值方法是基于样本点之间的相似程度创建一个插值拟合曲面,常见的有反距离加权法(IDW)、径向基函数法(RBF)等。地质统计学方法是根据样本点之间的统计规律与空间自相关性确定待测点的值,常见的有克里金(Kriging)插值法等。因此,可以根据已知点之间是否具有空间自相关性选择插值方法,从而达到更好的插值效果。

1.1.1 反距离加权法(IDW)

反距离加权法的表达式[15]为

(1)

式中,f(x,y)为待测点的属性值;n为已知点的个数;zi为第i个已知点的属性值,i=1, 2, …,n;di为(x,y)到(xi,yi)的欧式距离;p为幂指数,通常取值为1～3,但一般p为2时获得的效果更好。

1.1.2 径向基函数法

径向基函数较适合处理空间各向同性的问题。常见的径向基函数有薄板样条函数、张力样条函数、规则样条函数、高次曲面样条函数及反高次曲面样条函数[16],其表达式为

(2)

式中,ai为待定系数;di为欧式距离;φ为一种径向基函数。

1.1.3 Kriging法

Kriging算法包括普通克里金、简单克里金、泛克里金、协同克里金等[17]。半变异函数是克里金插值的一种模型,用于评估样本点之间的空间自相关性[18]。常用的半变异函数模型有球面模型、指数模型、高斯模型等。克里金插值算法的数学表达式为

(3)

式中,f(x,y)为待测点的属性值;zi为第i个已知点的属性值;λi为权重系数。能够满足点(xo,yo)处的插值f(xo,yo)与真实值zo差值最优系数为

(4)

(5)

同时,满足无偏估计条件为

(6)

(7)

1.1.4 空间自相关分析

在对地质体插值时,利用变异函数作为判断钻孔数据是否空间自相关的工具[19]。变异函数主要包括变程、基台值、块金值、偏基台值,如图2所示。通过绘制半变异函数云图表示自相关分析结果。半变异云图反映了采样点与相邻采样点之间的空间关系,距离越大,半变异值越大,空间相关性也越好。

图2 半变异函数结构

1.1.5 交叉验证与精度评价

为减少主观判断对差值结果客观性的影响,选择留一法进行交叉验证,并采用平均误差(ME)和均方根误差(RMSE)对插值结果进行精度评价[20]。

平均误差是预测值与实际值之间插值的平均值,计算公式为

(8)

均方根误差是用于描述预测值与实际值之间的偏离程度,计算公式为

(9)

由于自然界的地质体形状各异,任何插值方法都无法适应所有的地质体,因此需要对研究区进行具体分析,从而选出最合适的插值方法。本文首先对钻孔分层数据进行自相关分析,然后选择相应的插值方法。

1.2 地质体三维建模

建筑信息模型(BIM)用于解决复杂的三维地质体的可视化。文献[21—22]提出了基于BIM信息模型的三维滑坡地质灾害预测;文献[23]采用BIM软件,利用钻孔数据和地质剖面数据建立了地质体三维模型。

1.2.1 Revit和Dynamo可视化编程

近年来,越来越多的人基于BIM技术在地质体三维建模方面进行研究[24-27]。Revit软件提供了强大的族功能,Dynamo是Revit的一个可视化编程插件,具有强大的节点库,包含运算、分析、图像等功能,节点可以通过逻辑关系进行连接实现可视化,也可以采用Python等编程语言达到更多目的[28]。

1.2.2 建模流程

由于研究区范围小、地质构造相对简单,不涉及裂缝和倒置地层,因此本文基于钻孔数据使用Revit进行地质体三维建模。该方法首先通过对分层钻孔数据进行插值处理,选出每层数据最佳的插值方法;然后根据插值结果生成虚拟钻孔;最后基于Dynamo生成地质模型。具体建模流程如图3所示。

图3 建模流程

2 矿山储量估算系统设计

矿产储量估算系统采用C/S的结构模式,在Visual Studio 2013 环境下,运用C#语言,在SuperMap平台上开发,将矿山地质体三维模型与矿产储量估算系统进行融合。系统包含数据管理、矿山二三维显示、储量估算3个模块。数据管理模块包括矿区的DOM、DEM及钻孔数据的管理;矿山二三维显示模块包括矿体地质体三维模型、三维地形模型、矢量模型的显示,以及对模型的交互浏览、平移、缩放、旋转等操作。储量估算模块包括对开采量的计算及报表输出等。系统功能结构如图4所示。

图4 系统功能结构

3 试验与分析

3.1 研究区概况

研究区为鹤壁市二道庄矿区,面积约为5.2 km2, 海拔高度为150～440 m,如图5所示。

图5 研究区

3.2 空间数据采集与处理

利用无人机采集空间数据,经内业处理得到矿区的DEM与DOM。共收集40个工程地质钻孔,钻孔深度为40.89～121.56 m。原始钻孔柱状图中包含钻孔编号、孔口坐标、开孔日期、终孔日期、终孔深度、分层情况、岩性描述等信息。钻孔数据坐标采用CGCS2000。得到钻孔数据空间分布如图6所示。

图6 钻孔分布

3.3 三维地质体构建

3.3.1 钻孔插值及精度检验

钻孔数据按照岩性分为白云质灰岩、鲕状白云岩、白云岩、鲕状白云岩、灰岩共5层,分别用Z0、Z1、Z2、Z3、Z4表示,并导入数据库中统一管理。计算各地层数据的变异函数,并绘制变异函数云图如图7所示,据此分析各层数据的空间自相关性。

图7 各地层变异函数云图

由图7可知,Z0、Z2、Z3地层的云图点较散乱,变异函数值随已知点之间距离的增大呈先增大、后减小的趋势,不符合空间自相关的云图特点,即不具备空间自相关性;Z1、Z4的变异函数值随距离的增大有明显的增大趋势,符合空间自相关的特性。因此,对Z1、Z4采用克里金插值、径向基函数插值及反距离加权插值;对Z0、Z2、Z3则采用径向基函数插值和反距离加权插值。其中,克里金插值选取普通克里金插值,变异函数选择球面函数,径向基函数插值选取规则样条函数和张力样条函数两种插值方法。插值结果对比见表1、表2。

表1 Z1、Z4插值结果

表2 Z0、Z2、Z3插值结果

由表1可知,Z1、Z4地层较平坦,高低点起伏不大,整体呈上(北)高、下(南)低的情况。反距离加权插值得到的结果为,在高点和低点周围出现明显的凹包和凸包,“牛眼”现象较为明显,基本符合反距离加权插值的特点。克里金、规则样条、张力样条3种插值得到的结果大致相同,但克里金插结果较平滑,在高低点处的过渡较平缓,优于规则样条和张力样条插值;规则样条和张力样条插值结果,在Z1和Z4地层上部出现轻微的“牛眼”现象。综上分析可知,在符合空间自相关时,克里金插值优于其他3种插值方法。

由表2可知,Z0、Z2、Z3地层的地形起伏较大,其反距离加权插值都出现较明显的“牛眼”现象,特别是较高和较低区域的“牛眼”更加突显;规则样条和张力样条插值法在Z0和Z3地层的插值结果都较平滑,在Z2层两种方法都出现轻微的“牛眼”,这是由于该层局部起伏较大。

对插值结果进行对比后,对钻孔数据采用留一法验证各插值方法的精度,结果如图8所示,其中数据点到红色实线的距离表示各插值方法的误差值。

图8 留一法验证结果

由Z1和Z4地层的折线图可以看出,克里金插值结果最贴合红线,即最接近实测值;反距离加权插值结果误差最大;规则样条和张力样条结果的误差介于两者之间。结果表明,当地层符合空间自相关时,克里金插值方法优于其他插值方法。由Z0、Z2、Z3地层的折线图可以看出,由于3个地层的起伏较大,所有插值方法都与实测值有较大的误差。Z0和Z3中张力样条、规则样条、克里金插值折线大致呈相同趋势,几乎贴近于实测值,反距离加权法误差最为明显;Z2中张力样条、规则样条、克里金插值整体上优于反距离加权插值,但在个别点误差大于反距离加权插值。

对不同插值方法的插值精度用平均误差(ME)和均方根误差(RMSE)进行评价,结果见表3。张力样条插值法在Z0地层ME和RMSE最小,精度最好;Z1地层克里金插值法的ME和RMSE最小;Z2选择规则样条插值法的ME和RMSE最小;Z3和Z4交叉验证结果较特殊,Z3规则样条插值法的RMSE最小,张力样条插值法的ME最小,Z4克里金插值法的RMSE最小,规则样条插值法的ME最小。

表3 交叉验证结果 m

为了选择适用于各地层最优的插值方法,根据插值结果、交叉验证与精度评价可知:Z0选择张力样条插值法,Z1和Z4选择克里金插值法,Z2和Z3选择规则样条插值法。然后根据规则格网生成点模型,以便后边生成地质体模型。

3.3.2 地质体模型生成

将插值加密得到的虚拟钻孔点数据整理成Excel数据并导入Dynamo中,通过Topography.ByPoints节点生成虚拟钻孔点模型。生成的点模型如图9所示。

图9 地层点模型

将生成的地层点模型通过Topodraphy.PolySurface节点生成地层模型,如图10所示。通过Surface.PerimeterCurves节点生成地形曲面的边界,并将边界向上向下投影到一定的平面。

图10 地层面模型

将生成的地层面模型及两个投影面进行Loft放样,生成体模型,通过Geometry.SplitByTools节点与地形曲面进行裁剪,得到地质体模型,如图11所示。通过Springs.FamilyInstance.ByGeometry节点将生成的地质体导入Revit软件,可在Revit软件中更改各地层的属性,如图12所示。

图11 地质体模型

图12 Revit中三维地质模型

3.4 矿山储量估算系统实现与应用

将三维地质模型导入SuperMap中,与DOM和DEM生成的三维地形模型进行融合,存储至数据库中,基于C#进行二次开发,生成矿山储量估算系统。该系统从2022年初在矿区应用以来,满足了生产需求,提高了工作效率。全区储量报表如图13所示。

图13 全区储量报表

4 结 语

为满足露天矿山生产管理的实际需求,本文提出了一种基于混合插值算法的露天矿地质体三维建模方法,并据此进行矿山储量估算,为矿山企业实现生产科学管理提供了数据支撑。得出结论如下:①将钻孔数据分层,结合空间自相关分析,选出每个地层最适合的插值方法,不局限于一个区域的插值只选取单一的插值方法,提高了钻孔插值的精度,进而提高了地质体三维建模的精度;②通过Revit与Dynamo,采用“点—面—体”的建模思路构建三维地质体模型,基于Revit平台构建的模型可以与其他BIM及GIS软件进行交互。虽然本文方法及研究成果在露天矿山得到应用,但对于地质结构复杂情况的露天矿产储量估算仍需进一步研究。