CN114332413A - 一种基于滑动克里金插值的地质体建模方法及装置 - Google Patents

Abstract

本发明公开了一种基于滑动克里金插值的地质体建模方法及装置，所述方法包括：获取地质数据，所述地质数据为结合属性数据与空间数据的四维数据；将每一组四维数据作为初始采样点，根据初始采样点的三维坐标计算深度值；对所有数值点坐标排序分组，获取变程dist；根据变程dist划分网格并确定每个采样点所属的网格；根据克里金方程获取待插点的估计值，并剔除深度小于0的插值点，得到插值点与采样点集合；将插值点与采样点集合的深度值转化为高程值；基于插值点和采样点构建三角网，得到具有属性数据的地质曲面模型。本发明提高了对数据的插值速度和精度，提升了地质体建模效果。

Description

一种基于滑动克里金插值的地质体建模方法及装置

技术领域

本发明属于计算机图形处理技术领域，具体涉及一种地质体建模方法及装置。

背景技术

地质体建模是依赖于地质勘探数据，通过计算机图形学和计算机理论的方法在计算机上展示出来，为地质勘探工程提供解决方案和策略的一种方法。

在地质勘探数据采集过程中，通常采用抽样采集和一定距离采集，达不到对每个空间位置点采集数据，采样获得的空间数据呈现离散分布。达不到地质建模所需要的数据量要求。因此需要采用空间插值方法对每个空间位置进行插值估计，以达到加密空间数据和规格化网格分布的效果。

克里金插值法是一种较为常用的空间插值方法，它以变差函数理论为基础，对限定区域内的变量做线性无偏的最优估计，参考输入离散三维空间数组的参数，以及待插点和离散的已知点位置关系，通过对已知点进行数据分析，对待插点的数据进行最优估计，在插值效果上较符合实际，精确度高。传统克里金插值是将所有己知点作为影响因子对待插点进行估算。克里金插值的关键在于求解克里金方程组，当选择全部已知点对任一待插点进行估算时，显然需要求解一个多维矩阵，多维矩阵的求解过程时间长、效率低，而且当采样规模很大时，计算机甚至无法求解。另外，地质统计学中的协方差是一个估计量，由于实际测量数据往往包含许多不确定性，可能导致在大范围内插值所产生的均方差比较小邻域内产生的均方差大很多。因此在进行待插值计算之前有必要首先确定一个邻域，以保证在划分邻域内的己知点对待插点的估计产生影响，而邻域外数据产生影响不大或者没有影响。

空间属性插值即通过已知采样点对任一坐标点(x,y,z)处的坐标值进行属性值估计得到(x,y,z,p)。然而，实践证明当三维空间数据加入属性数据即进行四维空间数据体插值时，该方法存在如下问题：其一，因为属性数据体规模通常很庞大，在克里金插值的最初进行变程值的计算时就存在困难，无法完成插值运算，建模时无法获得每个空间位置的估计值。其二，采样点往往呈现离散分布，分布不均的特点。通过原始克里金插值法插值出来的空间位置及属性值当采样点分布不均时，通过该方法得到的插值结果不精确，建模效果与实际建模需求误差较大。

发明内容

发明目的：针对现有技术的不足，本发明的目是提供一种基于滑动克里金插值的地质体建模方法及装置，利用更优秀的插值方法提高地质体建模的精度和效率。

技术方案：为了实现上述发明目的，本发明提出的一种基于滑动克里金插值的地质体建模方法，包括以下步骤：

获取地质数据，所述地质数据为结合属性数据与空间数据的四维数据；

将每一组四维数据作为初始采样点，根据初始采样点的三维坐标计算深度值；

对所有数值点坐标排序分组，获取变程dist；

根据变程dist划分网格并确定每个采样点所属的网格；

根据克里金方程获取待插点的估计值，并剔除深度小于0的待插点，得到插值点与采样点集合；

将插值点的深度值转化为高程值；

基于插值点和采样点构建三角网，得到具有属性数据的地质曲面模型。

其中，根据初始采样点的三维坐标计算深度值包括：将每个初始采样点三维坐标作为起点，方向向量Vec(0，0，1)作Z轴正方向的射线，求射线与上顶面的交点，交点坐标(x，y，z′)，空间节点坐标深度值D＝Z′-Z。

进一步地，对所有数值点坐标排序分组，获取变程dist包括：

根据x，y方向上的坐标计算所有采样点相互之间的距离h_ij和待插点与采样点的距离h_i0，h_ij和h_i0称为距离变量；

根据所有求得的距离变量，获取最大距离h_max，最小距离h_min，得到分组划分的变程

其中S为划分的组数。

进一步地，根据变程dist划分网格并确定每个采样点所属的网格包括：

根据变程dist对区域进行网格划分，得到大网格；

以指定的划分步长对每个大网格进行二次网格划分，根据已知采样点坐标，在二层网格划分中进行数据点分配，确定每个采样点在二层网格所属的网格单元；

使用漫水法对网格进行搜索，将搜索到的点填充大网格，直到大网格内的点数量满足预设的最小点数量。

进一步地，所述划分步长根据下式计算得到：

step＝sqrt(N*w*h/n)/k

式中sqrt表示平方根运算，N为采样点总数，w为待划分区域的宽度，h为待划分区域的长度，n为每个大网格内采样点数量，每个大网格又被划分为k²个小的网格。

进一步地，根据克里金方程获取待插点的估计值包括：

根据克里金方程组：

计算得到每个插值点的权重值

将计算出的权重值α_i带入公式

得到插值点t₀的估计值Z^*(t₀)；

其中，第i个空间点t_i的属性值Z(t_i)为相应的深度值D_i，则得到待插点t₀深度的估计值D₀；i∈[1，N]；γ_ij代表γ(h_ij)，h_ij表示第i个采样点和第j个采样点之间的距离，

E表示数学期望，μ为一辅助参数。

进一步地，将插值点的深度值转化为高程值包括：以当前插值点坐标(x₀，y₀，D₀)为起点，方向向量Vec(0，0，1)作Z轴正方向的射线，求该射线与上顶面的交点，得到交点坐标(x₀，y₀，z₀′)，则z₀＝z₀′-D₀，得到最终的插值点坐标为(x₀，y₀，z₀)。

本发明还提供一种基于滑动克里金插值的地质体建模装置，包括：

地质数据获取模块，用于获取地质数据，所述地质数据为结合属性数据与空间数据的四维数据；

深度值转换模块，用于将每一组四维数据作为初始采样点，根据初始采样点的三维坐标计算深度值；

变程确定模块，用于对所有数值点坐标排序分组，获取变程dist；

网格划分模块，用于根据变程dist划分网格并确定每个采样点所属的网格；

插值模块，用于根据克里金方程获取待插点的估计值，并剔除深度小于0的待插点，得到插值点与采样点集合；

高程转换模块，用于将插值点的深度值转化为高程值；

地质模型构建模块，用于基于插值点和采样点构建三角网，得到具有属性数据的地质曲面模型。

本发明还提供一种计算机设备，包括：

一个或多个处理器；

存储器；以及

一个或多个程序，其中所述一个或多个程序被存储在所述存储器中，并且被配置为由所述一个或多个处理器执行，所述程序被处理器执行时实现如上所述的基于滑动克里金插值的地质体建模方法的步骤。

本发明还提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现如上所述的基于滑动克里金插值的地质体建模方法的步骤。

有益效果：通过对耗时的分析可知克里金插值需要求解方程组，计算量大、速度慢。不管是三维高程插值还是四维属性数据体插值，本发明改进的克里金算法相比普通克里金和梅钢等人提出的滑动克里金在不过多损失精度的情况下具有明显的速度优势。本发明在四维属性数据体规模庞大时，能够更快地实现对高程数据和属性数据体的插值，从而提高地质体建模效率和效果。

附图说明

图1为根据本发明的基于滑动克里金插值的地质体建模方法流程图；

图2为普通克里金方法一层网格划分示意图；

图3为根据本发明的改进的滑动克里金领域划分示意图；

图4为根据本发明的克里金插值效果显示；

图5为根据本发明的基于深度的滑动克里金领域划分效果显示；

图6为根据本发明的最终地质体建模效果图。

具体实施方式

下面结合附图对本发明的技术方案作进一步说明。

参照图1，本发明的实施方式中，一种基于滑动克里金插值的地质体建模方法，包括以下步骤：

步骤S1，获取地质数据。

钻孔是地质体建模中的一个重要组成部分，是一种结合了属性数据与空间数据的四维数据结构体。地质钻孔是人类利用机器向地下挖的深井，空间数据包括钻孔的三维坐标，属性数据包括土层名称、土层颜色、土层厚度、深度等。在地质勘探过程中，通过钻孔，获取到相应区域的地质属性，再录入地质数据库中。本发明实施例中选取的钻孔数据为华北地台潘庄凸起构造单元研究区的钻孔数据库数据，作为实验的种子节点。

步骤S2，以一组空间四维点作为初始采样点，根据采样点三维坐标计算深度值。

将每个初始空间采样点三维坐标作为起点，方向向量Vec(0，0，1)作Z轴正方向的射线，求射线与上顶面的交点，交点坐标(x，y，z′)，空间节点坐标深度值D＝Z′-Z，更新节点坐标高程值Z为深度值D。其中上顶面为建立当前地层面时，位于上方的已存在的面。此时得到的位置空间数据是三维坐标和深度值构成的。

步骤S3，对所有数值点坐标排序分组，获取变程dist。

计算所有采样点相互之间的距离h_ij和待插点与采样点的距离h_i0，h_ij和h_i0称为距离变量。假设采样点i的坐标为(x_i，y_i，z_i)，采样点j的坐标为(x_j，y_j，z_j)，则采样点i和采样点j之间的距离为

采样间距是影响变差函数的唯一变量，而与各采样点所在的坐标位置没有直接关系。

根据求得的所有距离变量进行分组。设最大距离为h_max，最小距离为h_min，假设划分为S组，则划分的变程为

每个距离变量h_ij所在的分组可以用

表示。

步骤S4，划分网格并确定每个采样点所属的网格。

根据得到的变程dist，确定所有采样点所属的网格单元。通过网格划分区域在x，y方向上的最小边界值和划分步长可以确定每个网格单元所在的行和列。

将当前采样点按其位置加入到相应网格列表中，并在普通克里金划分方法进行一层网格划分中再进行更细致的网格划分，限定网格邻域内的最大点数和最小点数，假设种子节点(即采样点)数量为N，待划分区域的长和宽可以通过种子节点坐标确定范围，假设长为w，宽为h。如果每个大网格内的种子节点数量为n，每个大网格又进一步被划分为k²个小的网格，在均匀的采用分布的情况下，可以得到步长公式为step＝sqrt(N*w*h/n)/k。

其中网格划分方法为：用变程dist对区域进行网格划分，再以步长为step/k进行二次网格划分，根据已知采样点坐标，在二层网格划分中进行数据点分配，确定每个采样点在二层网格中所属的网格单元。再使用漫水法对网格进行搜索，将搜索到的点填充大网格，直到大网格内的点数量满足最小点数量。图2是一层网格划分，把采样点划分为各个区域；假设k＝3，则图3是第二次网格划分示意图。

步骤S5，获取待插点的估计值。设A为整个研究区域，t为研究区域内表示位置的随机变量，空间点t的属性值为Z(t)，Z(t)是定义在三维空间中的一个实函数，即区域化变量。区域化变量是地质统计学学术用语。以空间点t的三个直角坐标t_u，t_v，t_w为自变量的随机场，Z(t_u，t_v，t_w)＝z(t)称为一个区域化变量。E(t)表示t的数学期望。对研究区域A内的任一待插点t₀，如果该点区域化变量的真实值Z(t₀)是未知的，求得该点区域化变量的估计值z^*(t₀)并替换真实值。

(1)首先计算加权系数α_i。加权系数α_i即第i个位置处己知采样点的权重，表示己知采样点的值对待插点值的影响度或贡献度。α_i不仅取决于数据点之间的距离、计算位置，还取决于所有采样点的整体空间排列。

采用无偏估计，即对任一待插点t₀，其真实值Z(t₀)与估计值z^*(t₀)之间差的数学期望为零，公式表达式为E[z^*(t₀)-Z(t₀)]＝0。由平稳假设条件可知E[Z(t)]＝m为常数。根据公式

有：

从而得到关系式

采用最小估计方差，最小估计方差条件就是使所有待插点的估计值与真实值差的方差达到最小，即在权重系数一定时满足公式Var[Z(t₀)-Z^*(t₀)]＝min。该公式又可表示为：

结合无偏条件

并采用拉格朗日乘数法进行计算得到以下目标函数：

其中μ为一辅助参数，j和i表示第j和第i个采样点。对所有的α_i和μ依次进行一阶偏导数求解并令其为零，可得到N+1个方程：

因区域化变量满足二阶平稳假设时有E(Z(t_i)Z(t_j))＝C₀+m²的关系，又因为平稳条件下满足C(h)＝C(0)-γ(h)，其中h为距离差值，

为变差函数，是反应区域变化量Z(t)与Z(t+h)变异程度的量。C(h)称为协方差函数，C(0)为常数。在Z(t)的协方差存在且平稳时，有C(h)＝E[Z(t)Z(t+h)-m²，m＝E[Z(t)]为常数，C₀称为块金效应，是一个常数。可以得到克里金方程组：

其中γ_ij代表γ(h_ij)，h_ij表示第i和第j个采样点之间的距离差值。

计算方程得到每个插值点的权重值

(2)得到每个插值点的估计值。根据

在计算出每个采样点对当前插值点的权重值后，代入公式

中，可得到插值点t₀的估计值Z^*(t₀)。以此类推，得到每个插值点属性的估计值。这里使每个采样点的属性值Z(t_i)为相应的深度值D_i，其中i∈[1，N]，可以得到待插点t₀深度的估计值D₀，当前待插点的坐标更新为(x₀，y₀，D₀)。x₀、y₀、D₀分别表示当前待插点的x，y方向的坐标及深度，以此类推，可得到每个待插点的深度，并更新相应的坐标。同理可得到每个待插点的各项属性值数据。在本文中，待插点和插值点可互换地使用。

步骤S6，对深度小于0的插值点进行处理，得到插值点与采样点集合。

循环遍历每个插值点，过滤深度值小于0的插值点，假设此时每个插值点坐标为(x₀，y₀，D₀)，当D₀＜0时，D₀＝0。

步骤S7，利用上顶面将待插点的深度值转化为高程值。

从第一个待插点开始，以当前待插点坐标(x₀，y₀，D₀)为起点，Vec(0，0，1)为方向向量做射线，求该射线与上顶面的交点，得到交点坐标(x₀，y₀，z₀′)，则z₀＝z₀′-D₀，得到最终的待插点坐标为(x₀，y₀，z₀)。

步骤S8，构建曲面。使用所有种子节点和得到的插值点作为最终节点，每个节点与满足一定条件的两个邻居节点构建三角形，构建三角形的条件为最小内角最大化和最大外接圆最小化。得到delauney三角网，delauney三角网中每个三角形都满足上述两个条件。将每个面赋予相对应的属性值，最终得到一个具有高精度属性数据的地质曲面模型。

因四维属性数据体规模庞大，普通克里金和梅钢等人提出的滑动邻域克里金方法无法完成插值运算，本发明提供了一种基于深度插值的滑动克里金插值方法，通过在划分一层网格单元的基础上进行了更细小的二层网格划分，并且此处对邻域内规定了最大点数和最小点数要求，把基于高程插值改为基于深度插值，并筛选节点，保证了面的任意一个节点深度都小于位于上方的面的节点深度。经实验该条件下进行克里金插值既缩短了运算时间，又保证了插值效果以及最终建模的合理性。

实验中采用C++/CLR编程语言和vs2017以及Bentley平台进行仿真，其中，计算机具体配置为CPU-i7，DDR-8G，SATA1-TB，OperatingSystem-Win10。

图4是基于普通克里金插值效果示意图，图5是基于深度的滑动克里金领域划分效果在三维中显示。实验数据表明普通克里金方法无法完成四维属性数据体的插值运算，而本发明改进的基于深度插值的滑动克里金插值可实现四维属性数据体插值且插值速度相比反距离加权插值法和层次B样条插值插值法在可接受范围内，但精确度却高出很多，如表中所示。

图6是最后建模的效果展示。

基于方法实施例相同的技术构思，本发明的另一实施例还提供一种基于滑动克里金插值的地质体建模装置，包括：

地质数据获取模块，用于获取地质数据，所述地质数据为结合属性数据与空间数据的四维数据；

深度值转换模块，用于将每一组四维数据作为初始采样点，根据初始采样点的三维坐标计算深度值；

变程确定模块，用于对所有数值点坐标排序分组，获取变程dist；

网格划分模块，用于根据变程dist划分网格并确定每个采样点所属的网格；

插值模块，用于根据克里金方程获取待插点的估计值，并剔除深度小于0的插值点，得到插值点与采样点集合；

高程转换模块，用于将插值点与采样点集合的深度值转化为高程值；

地质模型构建模块，用于基于插值点和采样点构建三角网，得到具有属性数据的地质曲面模型。

应当理解，本发明实施例中的基于滑动克里金插值的地质体建模装置可以实现上述方法实施例中的全部技术方案，其各个功能模块的功能可以根据上述方法实施例中的方法具体实现，其具体实现过程可参照上述实施例中的相关描述，此处不再赘述。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。

Claims (10)

1.一种基于滑动克里金插值的地质体建模方法，其特征在于，包括以下步骤：

获取地质数据，所述地质数据为结合属性数据与空间数据的四维数据；

将每一组四维数据作为初始采样点，根据初始采样点的三维坐标计算深度值；

对所有数值点坐标排序分组，获取变程dist；

根据变程dist划分网格并确定每个采样点所属的网格；

根据克里金方程获取待插点的估计值，并剔除深度小于0的插值点，得到插值点与采样点集合；

将插值点与采样点集合的深度值转化为高程值；

基于插值点和采样点构建三角网，得到具有属性数据的地质曲面模型。

2.根据权利要求1所述的基于滑动克里金插值的地质体建模方法，其特征在于，根据初始采样点的三维坐标计算深度值包括：将每个初始采样点三维坐标作为起点，方向向量Vec(0,0,1)作Z轴正方向的射线，求射线与上顶面的交点，交点坐标(x,y,z’)，空间节点坐标深度值D＝Z’-Z。

3.根据权利要求1所述的基于滑动克里金插值的地质体建模方法，其特征在于，对所有数值点坐标排序分组，获取变程dist包括：

根据x，y方向上的坐标计算所有采样点相互之间的距离h_ij和待插点与采样点的距离h_i0，h_ij和h_i0称为距离变量；

根据所有求得的距离变量，获取最大距离h_max，最小距离h_min，得到分组划分的变程

其中S为划分的组数。

4.根据权利要求1所述的基于滑动克里金插值的地质体建模方法，其特征在于，根据变程dist划分网格并确定每个采样点所属的网格包括：

根据变程dist对区域进行网格划分，得到大网格；

以指定的划分步长对每个大网格进行二次网格划分，根据已知采样点坐标，在二层网格划分中进行数据点分配，确定每个采样点在二层网格所属的网格单元；

使用漫水法对网格进行搜索，将搜索到的点填充大网格，直到大网格内的点数量满足预设的最小点数量。

5.根据权利要求4所述的基于滑动克里金插值的地质体建模方法，其特征在于，所述划分步长根据下式计算得到：

step＝sqrt(N*w*h/n)/k

式中sqrt表示平方根运算，N为采样点总数，w为待划分区域的宽度，h为待划分区域的长度，n为每个大网格内采样点数量，每个大网格又被划分为k²个小的网格。

6.根据权利要求5所述的基于滑动克里金插值的地质体建模方法，其特征在于，根据克里金方程获取待插点的估计值包括：

根据克里金方程组：

计算得到每个插值点的权重值

将计算出的权重值α_i带入公式

得到插值点t₀的估计值Z^*(t₀)；

其中，第i个空间点t_i的属性值Z(t_i)为相应的深度值D_i，则得到待插点t₀深度的估计值D₀；i∈[1,N]；γ_ij代表γ(h_ij)，h_ij表示第i个采样点和第j个采样点之间的距离，

E表示数学期望，μ为一辅助参数。

7.根据权利要求1所述的基于滑动克里金插值的地质体建模方法，其特征在于，将插值点的深度值转化为高程值包括：以当前插值点坐标(x₀,y₀,D₀)为起点，vec(0,0,1)为向量做射线，求该射线与上顶面的交点，得到交点坐标(x₀,y₀,z₀′)，则z₀＝z₀′-D₀，得到最终的插值点坐标为(x₀,y₀,z₀)。

8.一种基于滑动克里金插值的地质体建模装置，其特征在于，包括：

地质数据获取模块，用于获取地质数据，所述地质数据为结合属性数据与空间数据的四维数据；

深度值转换模块，用于将每一组四维数据作为初始采样点，根据初始采样点的三维坐标计算深度值；

变程确定模块，用于对所有数值点坐标排序分组，获取变程dist；

网格划分模块，用于根据变程dist划分网格并确定每个采样点所属的网格；

插值模块，用于根据克里金方程获取待插点的估计值，并剔除深度小于0的待插点，得到插值点与采样点集合；

高程转换模块，用于将插值点的深度值转化为高程值；

地质模型构建模块，用于基于插值点和采样点构建三角网，得到具有属性数据的地质曲面模型。

9.一种计算机设备，其特征在于，包括：

一个或多个处理器；

存储器；以及

一个或多个程序，其中所述一个或多个程序被存储在所述存储器中，并且被配置为由所述一个或多个处理器执行，所述程序被处理器执行时实现如权利要求1-7中任一项所述的基于滑动克里金插值的地质体建模方法的步骤。

10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1-7中任一项所述的基于滑动克里金插值的地质体建模方法的步骤。

Images