技术领域
[0001] 本
发明涉及互联网技术领域，特别设计一种基于克里金法矿体储量估算设计方法及系统。
背景技术
[0002] 我国地下矿产资源非常丰富，对底下矿物质开采过程主要包括矿物资源开发、矿物加工、安全工程。地质人员需要对矿山的开发进行规划和设计，所以如何对地质矿产资源进行合理的开发和科学的开采变得非常重要。地质勘探需要开采过程既满足开采需求又可以避免对周边地质环境的破坏，所以开采前期工作必须要做充足的准备，地下矿体储量的估算就是其中重要的一个环节。对地下
矿石分布构造进行分析需要地统计学相关理论进行研究。
[0003] 地统计学已经形成了一套完整的理论体系，克里金法是其研究的一种插值方法，估值方法主要以矿体储量品位的估计为目的，对矿物质空间地质结构分析。该理论以区域化变量为
基础理论，利用变异函数对区域化变量进行空间
位置分析，在估计方差最小的情况下进行品位值加权来估计出总的矿体储量。
[0004] 随着科技不断的发展，地质开采人员不仅仅依赖于人工或器械的手段对地质空间进行分析，如今大量地质开发会采用计算机和信息技术作为开采工具的技术
支撑。通过计算机
软件进行地质开采的管理与
可视化，使开采人员更加方便的操作和管理数据，也可以直观清晰的通过统计图像对数据统计进行分析。可视化也使地质工作者在分析过程中可以通过三维图像对矿床进行整体编辑与操作，三维的可视化改变了传统建模中平面图和剖面图的展示方式。
[0005] 本系统采用B/S结构进行开发，运用地质统计学方法实现空间矿体估值与应用，通过B/S 模式统一客户端使用，用
服务器集中实现系统关键功能的研发。大大方便了系统的开发、维护和使用。使传统的手工操作与C/S结构开发的软件变得更加便捷。目的在于减少地质工作人员的工作量，提高估值
质量，增加地质人员对估值过程的三维立体的分析。
发明内容
[0006] 为克服
现有技术的不足，本发明旨在最大限度的利用勘探过程所提供的
采样数据计算出矿体储量，使地质工作人员更加的方便的对矿山进行设计与开采。为此，本发明采取的技术方案是，基于克里金法矿体储量估算设计方法，步骤如下:
[0007] 导入钻孔数据，矿体数据；
[0008] 对数据进行特高品处理，并生成累计直方图与品位处理前后对比表格，特高品处理是将超过特高品味的数据统一设置为
指定值；
[0009] 对数据进行样长组合，并生成样长直方图与样长表格，将处理好的样长数据保存到
数据库中，样长是数据样段通过地质统计学方法分割为指定值长度，样长组合是将数据样段通过地质统计学方法分割为长度为指定值的样段作为以后操作所用的数据源；
[0010] 根据样长数据生成概率QQ图与PP图，QQ图可以用来在分布的位置-尺度范畴上可视化的评估参数，PP图是根据变量的累积概率对应于所指定的理论分布累积概率绘制的散点图，以便查看数据分布是否符合进行克里金估值的条件；
[0011] 根据样长数据计算得到实验变异函数，拟合实验变异函数生成球状模型与搜索椭球
体模型，其中，其中搜索椭球体模型用于在插值过程中确定数据的影响范围，球状模型为矿体储量计算提供参数；
[0013] 对矿体划分后得到的每个小块进行克里金插值估值。累加全部小块的估值结果得到最终结果。
[0014] 在三维展示中加载样长与椭球体模型，进行展示；对矿体进行块体划分具体采用八叉树对矿体进行块体划分。
[0015] 其中，克里金插值
算法是通过空间上所有已知点的采样数据加权求和获得未知点数据的值，权重由计算得到，式2表示克里金插值法，其区域化变量需满足二阶平稳假设：
[0016] (其中λ是权重，Z是已知点样品品位，n是样品点的个数)
[0017] 1)无偏性
[0018] 克里金在保证估值 是无偏估计，见式(3)：
[0019]
[0020] 当 时满足式(4)，其中m为区域化变量的期望：
[0021]
[0022] 其中m为区域化变量的期望；
[0023] 2)最优性
[0024] 同时满足估计方差最小的前提下求n个权重系数λi：
[0025]
[0026] 这时的估计方差为克里格方差，记为
[0027]
[0028] 其中x为待估点位置，为满足估计方差最小，通过使用拉格朗日数乘原理对式(6)求条件极值，然后与式(4)联立得到克里格方程：
[0029]
[0030] 由变异函数与协方差之间方程γ(h)＝c(0)-c(h)，则克里金方程用变异函数表示如下：
[0031]
[0032] 其中μ为估值结果，将式(8)展开可以用矩阵方式(9)表示：
[0033] Kλ＝M (9)
[0034] 其中
[0035]
[0036] 根据式(9)求得λ矩阵，λ矩阵的每一个元素即为相应已知点的采样数据的权重，累加过后即可求得估值结果。
[0037] 在块体模型构建成功后，对该块体模型进行遍历用克里金法求出每个块的矿体储量，然后汇总得到整个矿体的储量，在对单个块体进行插值的过程中需要确定由哪些样本点需要参与计算，对被估块体有影响的样本点都需要参与插值计算，影响范围即计算出的搜索椭球体的范围，将搜索椭球体中心转换到该块体中心进行搜索，得到插值的样本点集。
[0038] 并行计算的克里金插值法的构建流程如下：
[0039] 1)优化系数矩阵计算过程
[0040] 对K矩阵采用三
角分解LU分解法进行分解，LU分解将一个矩阵分解为两个单位三角矩阵的乘积，方便用来求逆矩阵。在程序设计中可以对三角矩阵进行压缩存储来减少运行时占用的内存，LU分解如式(10)：
[0041] K＝LU (10)
[0042] 其中L为下三角矩阵，U为上三角矩阵，将式(10)代入克里金等式中有(LU)λ＝M，则分解等式成两个等式(11)，先求解y再求解λ：
[0043]
[0044] 由此将系数矩阵λ代入克里金方程即可得出插值结果；
[0045] 2)并行计算克里金插值。
[0046] 其中，根据样长数据计算得到实验变异函数，拟合实验变异函数生成球状模型与搜索椭球体模型，具体步骤如下：
[0047] 首先找到与矿体走向一致的平面，需要根据用户设置的变异函数个数对平面进行划分，然后针对每个方向根据实验变异函数理论公式式1计算出实验变异函数：
[0048] 将样长数据间隔相同距离的点找出，将其样品值带入公式中求得相应的变异函数值，以变异函数值为y轴，间隔h为x轴绘制该变异函数折线图供用户选择主方向，将该方向设为搜索椭球体的
主轴方向。
[0049]
[0050] 其中，h指两样品之间的距离，Z(i)指i位置的样品品位，根据主轴方向找到其垂平面进行二次划分，并绘制次轴平面的实验变异函数折线图供用户选择，将用户选择的方向设为搜索椭球体的次轴方向；
[0051] 根据主轴方向和次轴方向生成与其垂直的短轴方向，并计算出该短轴方向的实验变异函数，并将该方向设为搜索椭球体的主轴方向；
[0052] 然后，将三个轴的实验变异函数拟合为对应的三个轴的球状模型，即对之前的折线图找到一条合适的圆滑曲线或直线，并运用球状模型函数来描述。
[0053] 最后，根据球状模型的变程得到搜索椭球体三轴的半径长度，再结合实验变异函数所选择的方向生成搜索椭球体模型。用该椭球体表示空间区域变量的变异性，反映不同方向的自相关距离。
[0054] 本发明的特点及有益效果是：
[0055] 本发明一种基于克里金法矿体储量估算设计方法及系统，可以在浏览器直接实现矿体估值及三维展示，最大限度的利用勘探过程所提供的采样数据计算出矿体储量与生成矿体模型，使地质工作人员更加的方便的对矿山进行设计与开采；本发明经过误差统计，误差不超过5％，符合
地质学要求；本发明提供的系统界面简单，对操作人员没有技术要求。
附图说明：
[0056] 为了更清楚地说明本发明
实施例的技术方案，下面将对实施例描述所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0057] 图1是本发明一种基于克里金法矿体储量估算系统的实验变异函数空间计算方向图(a) 主轴平面搜索(b)次轴平面搜索；
[0058] 图2是本发明一种基于克里金法矿体储量估算系统的样段可视化图(a)钻孔样段信息整体效果(b)样段局部信息效果；
[0059] 图3是本发明一种基于克里金法矿体储量估算系统的椭球体模型可视化图；
[0060] 图4是本发明一种基于克里金法矿体储量估算系统的克里金插值
流程图。
[0061] 图5是本发明一种基于克里金法矿体储量估算系统的技术
框架示意图。
[0062] 图6是本发明一种基于克里金法矿体储量估算系统的系统逻辑架构图。
[0063] 图7是本发明一种基于克里金法矿体储量估算系统的并行克里金插值流程图。
具体实施方式
[0064] 本发明一种基于克里金法矿体储量估算方法及系统,该方法包括:
[0065] 导入钻孔数据、矿体数据；
[0066] 对数据进行处理与统计分析；
[0067] 利用变异函数的计算对地质空间进行深入分析，拟合实验变异函数并生成球状模型与搜索椭球体模型，其中搜索椭球体模型用于在插值过程中确定数据的影响范围，球状模型为矿体储量计算提供参数。；
[0068] 通过八叉树思想对待估值矿体进行块体划分；
[0069] 采用并行计算克里金插值算法快速估算出矿体的储量；
[0070] 在三维展示模块加载样长和搜索椭球体模型；
[0071] 本发明提供一种基于克里金法矿体储量估算系统，该系统执行模块依先后顺序如下：数据导入模块，数据分析处理模块，拟合理论变异函数模块，克里金插值估算模块，三维展示模块。所述数据导入模块包括钻孔信息导入，矿体信息导入；所述数据分析处理模块包括特高品处理功能，样长组合功能，绘制累积直方图、直方图、QQ图、PP图功能；所述拟合实验变异函数模块包括拟合球状模型功能，生成搜索椭球体功能；所述克里金插值估算模块包括矿体块体划分功能、矿体估值功能；所述三维展示模块包括加载样长功能、加载搜索椭球体模型功能。
[0072] 进一步地，所述特高品处理功能是将超过特高品味的数据统一设置为指定值或者不进行特高品味处理；
[0073] 进一步地，所述样长组合功能是将数据样段通过地质统计学方法分割为长度为指定值的样段作为以后操作所用的数据源；
[0074] 进一步地，所述绘制累积直方图、直方图、QQ图、PP图功能是绘制相应的统计图；
[0075] 统计学里QQ图(Q代表分位数)是一个概率图，用图形的方式比较两个概率分布，把他们的两个分位数放在一起比较。首先选好分位数间隔。图上的点(x,y)反映出其中一个第二个分布(y坐标)的分位数和与之对应的第一分布(x坐标)的相同分位数。因此，这条线是一条以分位数间隔为参数的曲线。如果两个分布相似，则该QQ图趋近于落在y＝x线上。如果两分布线性相关，则点在QQ图上趋近于落在一条直线上，但不一定在y＝x线上。QQ图可以用来在分布的位置-尺度范畴上可视化的评估参数。
[0076] PP图是根据变量的累积概率对应于所指定的理论分布累积概率绘制的散点图，用于直观地检测样本数据是否符合某一概率分布。如果被检验的数据符合所指定的分布，则代表样本数据的点应当基本在代表理论分布的对角线上。
[0077] 进一步地，所述拟合球状模型功能和生成搜索椭球体功能是对数据进行计算变异函数操作，并生球状模型和搜索椭球体模型的过程。
[0078] 首先找到与矿体走向一致的平面，需要根据用户设置的变异函数个数对平面进行划分，如图1(a)所示。然后针对每个方向根据实验变异函数理论公式(式1)计算出实验变异函数：将样长数据间隔相同距离的点找出，将其样品值带入公式中求得相应的变异函数值，以变异函数值为y轴，间隔h为x轴绘制该变异函数折线图供用户选择主方向，将该方向设为搜索椭球体的主轴方向。
[0079] (其中h指两样品之间的距离，Z(i)指i位置的样品品位)
[0080] 根据主轴方向找到其垂平面进行二次划分(如图1(b)所示)，并绘制次轴平面的实验变异函数折线图供用户选择，将用户选择的方向设为搜索椭球体的次轴方向。
[0081] 根据主轴方向和次轴方向生成与其垂直的短轴方向，并计算出该短轴方向的实验变异函数，并将该方向设为搜索椭球体的主轴方向。
[0082] 然后，将三个轴的实验变异函数拟合为对应的三个轴的球状模型，即对之前的折线图找到一条合适的圆滑曲线或直线，并运用球状模型函数来描述。
[0083] 最后，根据球状模型的变程得到搜索椭球体三轴的半径长度，再结合实验变异函数所选择的方向生成搜索椭球体模型。用该椭球体表示空间区域变量的变异性，反映不同方向的自相关距离。其中，椭球体模型用于在后面插值过程中确定数据的影响范围。
[0084] 进一步地，所述矿体块体划分功能是通过八叉树思想对待估值矿体进行块体划分；
[0085] 块体模型在结合实体模型约束的条件下进行构建，它是一种面向空间的体模型。构建该块体模型的目的是对矿体进行分解，将其网格离散化生成三维栅格模型。构建矿体时应尽量对目标矿体进行
覆盖。
[0086] 进一步地，所述矿体估值功能是对上一步中块体划分得到的每一个小块进行克里金插值，最后统计所有小块估值结果得到最终估值结果，在下文会具体介绍；
[0087] 进一步地，所述加载样长功能是系统根据样长组合后数据库获取到样段起始坐标s与结束坐标e，使用three.js(Three.js是一款运行在浏览器中的3D引擎)绘制空间中的线段，并标记处样段中心点(如图2(b))。构建样段步骤如下：
[0088] 遍历从数据库获取到的样段集合。将坐标s与e转换为系统相对坐标。并计算出每个样段的中点坐标mid。
[0089] 绘制线段。使用three.js的Geometry(three.js中创建点的方法)构建线段
顶点，用 LineBasicMaterial(three.js中一种材质)设置线段材质。然后通过THREE.Line(three.js 中创建线段的方)构建线段模型。
[0090] 绘制中心点。用SphereGeometry(three.js中创建球体的方法)建立半径为0.1球，通过MeshLambertMaterial(three.js中的一种材质)方法设置球的材质，最后Mesh(将物体与材质结合创建模型的方法)方法建立球模型。
[0091] 进一步地，所述加载搜索椭球体模型功能是读取数据中的所述变异函数信息，然后在所述三维展示模块中绘制；
[0092] 椭球体模型在three.js中没有相应的绘制方法。本系统对球状模型SphereGeometry方法进行更改。自定义方法为NewESphereGeometry，该方法可以绘制椭球体几何形状。具体构建步骤如下：
[0093] 获取椭球体参数。得到椭球体三轴方向 与三轴长度a、b、c。
[0094] 通过event.which(three.js中获取点击事件的方法)属性获取右键点击的坐标点位置 P，将P设为椭球体中心点坐标，以P为中心进行绘制。
[0095] 由方向 绘制椭球体三轴。采用ArrowHelper(three.js中绘制箭头的方法)方法绘制箭头几何体模型，并设置其材质。
[0096] 建立椭球几何体。用NewESphereGeometry方法实现椭球在三维空间的
定位定向，并生成
几何模型。
[0097] 使用MeshBasicMaterial(three.js中的一种材质)设置椭球体材质。并生成椭球体模型(如图3所示)。
[0098] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
[0099] 如图4所示，为本发明一种基于克里金法矿体储量估算系统的流程图,该方法步骤包括:
[0100] 1、导入钻孔数据、矿体数据；
[0101] 2、对数据进行处理与统计分析；
[0102] 3、拟合实验变异函数生成球状模型与搜索椭球体模型；
[0103] 4、通过八叉树思想对待估值矿体进行块体划分；
[0104] 5、用克里金法对步骤4中得到的每个块体进行插值并累加求值，其中，搜索椭球体模型用于在插值过程中确定数据的影响范围，球状模型为矿体储量计算提供参数；
[0105] 6、在三维展示中加载样长、椭球体模型。
[0106] 如图5所示，是本发明一种基于克里金法矿体储量估算系统的技术框架示意图，系统开发过程中涉及两端，客户端和服务端，客户端通过JavaScript(JavaScript一种直译式脚本语言，是一种动态类型、弱类型、基于
原型的语言，内置支持类型)编写脚本，利用WebGL (是一种3D绘图协议)进行3D
渲染，使用html5(HTML5是Web中核心语言HTML的规范) 和css(层叠样式表是一种用来表现HTML或XML等文件样式的计算机语言)搭建前端页面，并使用Canvas(Canvas是HTML5中用于图形的绘制元素)进行二维设计。服务器端的数据存储到MongoDB数据库中，通过Beego框架(Beego是一个快速开发Go应用的HTTP框架) 建立服务，通过c#建立WebService(WebService是一种跨编程语言和跨
操作系统平台的远程调用技术)实现Dimine(Dimine是提供一系列矿体处理
接口的第三方平台)接口调用，客户端和服务端通过WebSocket(WebSocket是一种在单个TCP连接上进行全双工通信的协议) 实现通信连接。
[0107] 如图6所示，是本发明一种基于克里金法矿体储量估算系统的系统逻辑架构图，本系统采用基于MVC架构的Beego框架进行开发。后端开发语言主要为Go语言(Go是Google开发的一种静态强类型、编译型语言)，以Python语言(Python是一种跨平台的
计算机程序设计语言)辅助开发数据分析部分的功能。前端页面展示部分采用Html5、JavaScript、Ajax (是指一种创建交互式网页应用的网页开发技术)等前端技术，三维可视化部分采用基于 WebGL的three.js进行构建。使用MongoDB数据库(MongoDB是一个基于分布式文件存储的数据库)对数据进行存储。
[0108] 下面介绍克里金插值过程。如图7所示，是克里金插值流程图。
[0109] 克里金插值算法是通过空间上所有已知点的采样数据加权求和获得未知点数据的值，权重由计算得到。克里金插值法(式2)的区域化变量需满足二阶平稳假设。
[0110]
[0111] 3)无偏性
[0112] 克里金在保证估值 是无偏估计(式(3))。
[0113]
[0114] 当 (其中m为区域化变量的期望)时满足式(4)
[0115]
[0116] 4)最优性
[0117] 同时满足估计方差最小(式(5))的前提下求n个权重系数λi。
[0118]
[0119] 这时的估计方差为克里格方差，记为
[0120]
[0121] 为满足估计方差最小，通过使用拉格朗日数乘原理对式(6)求条件极值，然后与式(4)联立得到克里格方程。
[0122]
[0123] 由变异函数与协方差之间方程γ(h)＝c(0)-c(h)，则克里金方程可用变异函数表示如下：
[0124]
[0125] 将式(8)展开可以用矩阵方式(9)表示：
[0126] Kλ＝M (9)
[0127] 其中
[0128]
[0129] 根据式(9)可求得λ矩阵。λ矩阵的每一个元素即为相应已知点的采样数据的权重。累加过后即可求得估值结果。
[0130] 在块体模型构建成功后，对该块体模型进行遍历用克里金法求出每个块的矿体储量，然后汇总得到整个矿体的储量。在对单个块体进行插值的过程中需要确定由哪些样本点需要参与计算，一般情况下，对被估块体有影响的样本点都需要参与插值计算。影响范围即我们计算出的搜索椭球体的范围，将搜索椭球体中心转换到该块体中心进行搜索，得到插值的样本点集。
[0131] 影响范围对插值结果的准确度起了相当重要的作用，所选取的影响范围不同，矿体含量计算结果会有很大的误差。所以确定合理的搜索椭球体非常重要，需要做好前期工作，若搜索到的样本点集少于插值允许的最少点数量，则需要对搜索范围进行扩大，扩大倍数由用户定义。扩大次数本系统定位三次，若该搜索椭球体扩大三次都没有满足点数要求则放弃对该块体进行估值。
[0132] 下面介绍克里金插值的并行实现过程：
[0133] 串行克里金法通过遍历划分好的块体模型，然后对每个块体进行一一赋值。当前一个块体进行插值结束后才能对下一个块体进行插值，这样消耗了大量计算时间。所以在并行算法中应尽量提高各个块体估值的并行性能。
[0134] 设有N个待估块体，则需要N个协程并发对该矿体进行估值，并调用CPU多核并行处理这N个协程，设定CPU核数来提高并行效率，CPU会根据
进程合理分配计算。由于每个待估块体估值过程都是独立的互不影响，只有在初始阶段传入相同的样品点集和搜索约束，结束时将插值结果写入该块体即可。
[0135] 为了减少并行时整体插值的时间，可以从提高单个块体插值速度上考虑。对一个块体插值的关键一步是计算克里金系数矩阵，克里金等式Kλ＝M。当计算数据量较大时，矩阵的阶数随之增大，计算矩阵求逆则会消耗巨大的内存，同时计算的时间也会大幅度提高。从而要提高速度，优化大矩阵求逆是需要考虑的核心问题。
[0136] 基于并行计算的克里金插值法的构建流程如下：
[0137] 1.优化系数矩阵计算过程
[0138] 对K矩阵采用三角分解法(LU分解)进行分解，LU分解将一个矩阵分解为两个单位三角矩阵的乘积，方便用来求逆矩阵。在程序设计中可以对三角矩阵进行压缩存储来减少运行时占用的内存。LU分解如式(10)
[0139] K＝LU (10)
[0140] 其中L为下三角矩阵，U为上三角矩阵。
[0141] 将式(10)代入克里金等式中有(LU)λ＝M，则分解等式成两个等式(11)，先求解y再求解λ。
[0142]
[0143] 由此将系数矩阵λ代入克里金方程即可得出插值结果。其中LU分解算法如算法1所示。
[0144]
[0145] 2.并行计算克里金插值
[0146] 本算法阐述利用Go语言进行并发编程和调用多核CPU并行计算克里格的过程。Go的并发机制最小单位是协程，它可以让某个函数区别与其他函数独立运行。Go通过GOMAXPROCS 方法设置逻辑处理器个数来实现并行，但并不需要维护协程数量与逻辑处理器的对应关系， Go会通过调度器把协程平等分配到每个逻辑处理器上。具体执行伪代码如算法2所示。
[0147]
[0148]
[0149] 并行计算伪代码过程如下：
[0150] 先判断当前块体是否在矿体内部，如果在的话继续划分，直到这个块体不能再划分为止，对当前块体进行克里格插值，并将结果累加得到最终结果。使用WaitGroup(go语言中协程同步的方法)来阻塞主协程的执行，通过add(),wait()函数(go语言中协程同步的方法)管理协程直到所有的协程执行完毕后再进行主协程的执行。在计算过程中待估总值需要所有的协程都来
访问并
修改，当两个协程同时修改该变量时结果就会出错。这就需要对该变量进行加
锁，使用go语言同步包(sync)中的Mutex(go中用于对代码块加锁的方法)提供的加锁机制对valuation_results变量进行加锁。加锁过的变量不会同时被多个协程访问，只有一个协程修改完毕并写回内存后才会解锁，下一个协程才可以继续修改。
