CN112116708B - 一种获取三维地质实体模型的方法及系统 - Google Patents

一种获取三维地质实体模型的方法及系统 Download PDF

Info

Publication number
CN112116708B
CN112116708B CN202010954732.XA CN202010954732A CN112116708B CN 112116708 B CN112116708 B CN 112116708B CN 202010954732 A CN202010954732 A CN 202010954732A CN 112116708 B CN112116708 B CN 112116708B
Authority
CN
China
Prior art keywords
stratum
sampling point
preset
attitude
point
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010954732.XA
Other languages
English (en)
Other versions
CN112116708A (zh
Inventor
张宝一
刘肖莉
陈桥
佟勇强
黄岸烁
王丽芳
邓浩
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202010954732.XA priority Critical patent/CN112116708B/zh
Publication of CN112116708A publication Critical patent/CN112116708A/zh
Application granted granted Critical
Publication of CN112116708B publication Critical patent/CN112116708B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/29—Geographical information databases
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Geometry (AREA)
  • Algebra (AREA)
  • Computer Graphics (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明涉及一种获取三维地质实体模型的方法，包括：S1、针对与预先设定的空间相应的地层界面采样点数据和地层产状采样点数据，获取所述预先设定的空间的地层位势场函数；所述地层界面采样点数据包括：地层界面采样点的位置、地层界面采样点的场值；所述地层产状采样点数据包括：地层产状采样点的位置、地层产状采样点的梯度矢量值；S2、基于所述预先设定空间的地层位势场函数，建立所述预先设定的空间的三维地质实体模型。本发明能够根据不同的条件对三维实体模型进行分割或者合并，同时也在数据更新的情况下能够动态地获得新的三维实体模型。

Description

一种获取三维地质实体模型的方法及系统
技术领域
本发明涉及三维地质建模技术领域，尤其涉及一种获取三维地质实体模型的方法及系统。
背景技术
三维地层位势场是对地下地层结构的定量化表征，可以根据不同的条件对地层界面进行分割和合并，具有动态划分的能力，在数据更新的情况下可以自动对模型进行同步动态更新，是地质构造定量分析的重要手段。三维地下地层结构建模可分为显式建模和隐式建模两种方式。显式建模需要大量的人机交互来连接地质边界线以围成地质体的三维模型，受地质采样数据密度和地质要素的复杂性的影响较大，建模过程繁琐、建模效率不高，且模型难以更新。隐式建模采用插值的方法构建隐式的数学曲面来表达地层的边界，与显式建模相比，其在建模速度、可重复性、拓扑一致性上具有巨大的优势。隐式建模方法来构建三维地层位势场，是用插值方法通过已知地层界面采样点和地层产状采样点来联合构建位势场，插值方案有连续的协同克里格(Cokriging)方法、径向基函数插值(RBFs，RadiusBasic Functions)方法和离散化插值方法。
Cokriging和RBFs都使用地层界面的点数据和地层产状点数据插值建立地层界面的场函数，并提取场函数的等值面作为地层界面。RBFs插值算法是一系列用于表面重建和三维对象表示的插值函数，RBFs隐式建模方法对每个地层界面都要单独插值建立标量场函数。普通的RBFs隐式建模方法，通常只使用地层界面点数据，而没有采用地层产状点数据进行插值。目前的Hermite-Birkhoff径向基函数(HRBFs，Hermite-Birkhoff Radius BasicFunctions)隐式曲面建模方法，重构地层界面时不但要求在地层边界线上标量场的函数值相等(通常为0值)，还要求在边界线上的梯度也等于测量值，然而要求地层界线点的每个位置上都有精确的梯度测量值在地质调查是难以实现的，只能用最近的产状测量值或线的法线方向来近似计算待插值表面的法向量，即梯度方向。
目前的RBFs或HRBFs隐式曲面建模方法通过插值建立单个地层界面的标量场函数，不同地层界面之间的拓扑一致性难以维护，地层内部的场属性和产状等要素难以表达；而且要求地层界线点的每个位置上都有梯度测量值，只能用最近的产状或线的法线来近似计算得到，也会再造成地层界面的产状与实际情况不一致。
发明内容
(一)要解决的技术问题
鉴于现有技术的上述缺点和不足，本发明提供一种获取三维地质实体模型的方法。
(二)技术方案
为了达到上述目的，本发明提供一种获取三维地质实体模型的方法，包括：
S1、针对与预先设定的空间相应的地层界面采样点数据和地层产状采样点数据，获取所述预先设定的空间的地层位势场函数；
所述地层界面采样点数据包括：地层界面采样点的位置、地层界面采样点的场值；
所述地层产状采样点数据包括：地层产状采样点的位置、地层产状采样点的梯度矢量值；
S2、基于所述预先设定空间的地层位势场函数，建立所述预先设定的空间的三维地质实体模型。
优选的，所述步骤S1之前还包括：
S·0、针对预先设定的三维空间和与所述三维空间所对应的平面地质图、剖面地质图和地层柱状图，将所述三维空间划分为多个子空间；
S0、基于与所述三维空间对应的平面地质图和剖面地质图，获取子空间对应的地层界面采样点数据和地层产状采样点数据；
相应的，所述S1为：
基于与所述子空间相应的地层界面采样点数据和地层产状采样点数据，获取所述子空间的地层位势场函数。
优选的，所述步骤S0包括：
S01、基于与所述三维空间对应的平面地质图和剖面地质图，获取地层界面点数据和地层产状数据，然后将地层界面点数据和地层产状数据从二维到三维映射，获取三维空间中的点集数据；
所述点集数据包括地层界面点的位置和地层产状点的位置以及地层产状点的走向、倾向和倾角；
S02、基于与所述三维空间对应的地层柱状图，获取地层界面点的场值；所述地层界面点的场值为所述地层柱状图中的地层相对埋深值；
S03、基于地层产状点的走向、倾向和倾角，获取地层产状点的梯度矢量值；
S04、基于三维空间中地层界面点的位置、地层界面点的场值、地层产状点的位置以及地层产状点的梯度矢量值，获取地层界面采样点数据和地层产状采样点数据；
所述地层界面采样点数据为与所述三维空间中任一子空间所对应的地层界面点的位置和场值；
所述地层产状采样点数据为与所述三维空间中任一子空间所对应的地层产状点的位置和梯度矢量值。
优选的，所述步骤S1包括：
基于与所述子空间相应的地层界面采样点数据和地层产状采样点数据，根据预先设定的条件，采用公式(1)、公式(2)确定所述子空间的地层位势场函数；
其中，公式(1)为：
Figure BDA0002678214340000041
其中，公式(2)为：
Figure BDA0002678214340000042
其中，||p-pi||为子空间中任意点p到地层界面采样点pi点的欧式距离；||p-qj||为子空间中任意点p到地层产状采样点qj点的欧式距离；
f(p)为地层位势场函数；
Figure BDA0002678214340000043
地层位势场的梯度场函数；
N为地层界面采样点个数；
M为地层产状采样点个数；
Figure BDA0002678214340000044
为预先设定的径向基函数；
Figure BDA0002678214340000045
为Hamilton算子，且
Figure BDA0002678214340000046
Figure BDA0002678214340000047
为求偏导的运算符；
Figure BDA0002678214340000048
为对x方向求偏导的运算符；其中，x方向为纬度的方向；
Figure BDA0002678214340000049
为对y方向求偏导的运算符；其中，y方向为经度的方向；
Figure BDA00026782143400000410
为对z方向求偏导的运算符；其中，z方向为预先设定的垂直于x和y的方向；
H为Hessian算子，且
Figure BDA0002678214340000051
Figure BDA0002678214340000052
为求二次偏导的运算符；
Figure BDA0002678214340000053
为对x方向求二次偏导的运算符；
Figure BDA0002678214340000054
为对y方向求二次偏导的运算符；
Figure BDA0002678214340000055
为对z方向求二次偏导的运算符；
<，>为两个矢量的内积运算符；
αi为地层界面采样点的系数；
βj为地层产状采样点的矢量系数；
C(p)＝c1+c2px+c3py+c4pz为预先设定的一次多项式；
px为任意点p在以预先设定的xyz-o坐标系中x轴上的坐标；
py为任意点p在以预先设定的xyz-o坐标系中y轴上的坐标；
pz为任意点p在以预先设定的xyz-o坐标系中z轴上的坐标；
其中，预先设定的xyz-o坐标系以预先设定的点为原点，以纬度的方向为x轴的方向，以经度的方向为y轴的方向，以预先设定的垂直于x和y的方向为z轴方向；
c1为预先设定的截距系数；
c2为预先设定的px的系数；
c3为预先设定的py的系数；
c4为预先设定的pz的系数；
所述预先设定的条件为f(p)的二阶导数的函数值最小。
优选的，所述步骤S1具体包括：
S11、基于所述地层界面采样点数据和地层产状采样点数据，确定公式(1)和公式(2)中的参数系数αi、βj及c1、c2、c3、c4的具体值；
S12、基于所述参数系数αi、βj及c1、c2、c3、c4的具体值，确定所述预先设定子空间的地层位势场函数。
优选的，所述步骤S11包括：
S111、所述地层界面采样点数据和地层产状采样点数据作为任意点p及其对应的地层位势场值及其梯度场值分别代入式(1)和公式(2)，可得：
Figure BDA0002678214340000061
Figure BDA0002678214340000062
pk为地层界面采样点的位置；fk为地层界面采样点的场值；与所有地层界面采样点pi和所有地层产状采样点qj构成方程(3)；
qk为地层产状采样点的位置；gk为地层界面采样点的场值；与所有地层界面采样点pi和所有地层产状采样点qj构成方程(4)；
S112、根据预先设定的条件，确定正交条件：
Figure BDA0002678214340000063
Figure BDA0002678214340000064
Figure BDA0002678214340000065
为第j个地层产状采样点的矢量系数βj在x方向上的分量；
Figure BDA0002678214340000066
为第j个地层产状采样点的矢量系数βj在y方向上的分量；
Figure BDA0002678214340000067
为第j个地层产状采样点的矢量系数βj在z方向上的分量；
S113、基于所述公式(3)、(4)、(5)、(6)，获取公式(7)；
Figure BDA0002678214340000071
其中，Φ为N×N的矩阵，且Φ的元素为
Figure BDA0002678214340000072
且
Figure BDA0002678214340000073
Figure BDA0002678214340000074
为N×3M的矩阵，且
Figure BDA0002678214340000075
的元素为
Figure BDA0002678214340000076
且
Figure BDA0002678214340000077
HΦ为3M×3M的矩阵，且HΦ的元素为
Figure BDA0002678214340000078
且
Figure BDA0002678214340000079
Figure BDA00026782143400000710
其中元素
Figure BDA00026782143400000711
元素
Figure BDA00026782143400000712
元素
Figure BDA00026782143400000713
地层界面采样点的场值数据f＝[f1 f2…fN]T，地层产状采样点的场值梯度矢量数据g＝[g1 g2…gM]T；
S114、基于公式(7)，确定系数αi、βj及c1、c2、c3、c4的具体值。
优选的，所述步骤S12具体包括：
将所述确定的系数αi、βj及c1、c2、c3、c4的具体值，代入到公式(1)和(2)中，确定所述预先设定子空间的地层位势场函数；
所述预先设定子空间的地层位势场函数为：
Figure BDA00026782143400000714
Figure BDA0002678214340000081
其中，A为系数αi的具体值；B为系数βj的具体值；D(p)＝d1+d2px+d3py+d4pz；d1为c1的具体值；d2为c2的具体值；d3为c3的具体值；d4为c4的具体值。
优选的，所述步骤S2包括：
S21、根据所述子空间的地层位势场函数，得到三维网格点的地层位势场的场值及其梯度矢量值；
所述三维网格点为在预先设定的xyz-o坐标系中沿x、y、z轴方向，在三维空间中按预先设定的分辨率间隔Δx、Δy、Δz规则采样获得的位置点；
S22、基于所述三维网格点的地层位势场的场值及其梯度矢量值，采用预先设定的数字高程模型DEM和预先确定的子空间边界获取预先设定的地层位势场值的等势面；
S23、采用所述数字高程模型DEM将所述地层位势场值的等势面、预先确定的子空间边界面联合围成体，获取三维地质实体模型。
另一方面，本实施例提供一种获取三维地质实体模型的的系统，包括：
至少一个处理器；以及
与所述处理器通信连接的至少一个存储器，其中，所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行上述任一获取三维地质实体模型的方法。
(三)有益效果
本发明的有益效果是：本发明的一种获取三维地质实体模型的方法，由于将地层界面的相对埋深作为位势场的场值，地层产状数据作为位势场的梯度约束，通过地层界面点数据和地层产状数据构建地层位势场函数，并基于所述预先设定空间的地层位势场函数，建立所述预先设定的空间的三维地质实体模型，相对于现有技术而言，其可以更准确地表示现实中的三维地质情况。
附图说明
图1为本发明的一种获取三维地质实体模型的方法流程图；
图2为本发明的一种获取三维地质实体模型的方法示意图；
图3为本发明实施例中的平面地质图；
图4为本发明实施例中的地层柱状图；
图5为本发明实施例中的剖面地质图；
图6为本发明实施例中地层产状采样点的走向、倾向、倾角与梯度矢量的关系图；
图7为本发明实施例中地层位势场；
图8为本发明实施例中地层面三维实体模型；
图9为本发明实施例中地层三维块体模型。
具体实施方式
为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。
本发明实施例提出的一种获取三维地质实体模型的方法，将地层界面的相对埋深作为位势场的场值，地层产状数据作为位势场的梯度约束，通过地层界面点数据和地层产状数据构建地层位势场函数，并基于所述预先设定空间的地层位势场函数，建立所述预先设定的空间的三维地质实体模型，本实施例通过追踪位势场的不同等势面可以获得不同细节程度的三维地质实体模型，能够根据不同的条件对三维实体模型进行分割或者合并，同时也在数据更新的情况下能够动态地获得新的三维实体模型。本实施例能够根据不同的条件对三维实体模型进行分割或者合并，同时也在数据更新的情况下能够动态地获得新的三维实体模型。
为了更好的理解上述技术方案，下面将参照附图更详细地描述本发明的示例性实施例。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更清楚、透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。
参见图1，本实施例中的一种获取三维地质实体模型的方法，包括：
S1、针对与预先设定的空间相应的地层界面采样点数据和地层产状采样点数据，获取所述预先设定的空间的地层位势场函数；
所述地层界面采样点数据包括：地层界面采样点的位置、地层界面采样点的场值；
所述地层产状采样点数据包括：地层产状采样点的位置、地层产状采样点的梯度矢量值；
S2、基于所述预先设定空间的地层位势场函数，建立所述预先设定的空间的三维地质实体模型。
参见图2，本实施例中的一种获取三维地质实体模型的方法，能同时重构若干整合(及平行不整合)地层界面(参见图2中的地层界面1、地层界面2、地层界面3)，同时在地层内部的地层产状点的标量场梯度也符合测量的数据，真实还原三维地质空间的地层结构，克服现有RBFs和HRBFs隐式建模方法难以保证地层拓扑一致性的问题，又能根据新的输入数据动态更新三维地层位势场。
本实施例中优选的，所述步骤S1之前还包括：
S·0、针对预先设定的三维空间和与所述三维空间所对应的平面地质图、剖面地质图、地层柱状图(参见图3、图4、图5)，将所述三维空间划分为多个子空间；
S0、基于与所述三维空间对应的平面地质图、剖面地质图，获取子空间对应的地层界面采样点数据和地层产状采样点数据；
相应的，所述S1为：
基于与所述子空间相应的地层界面采样点数据和地层产状采样点数据，获取所述子空间的地层位势场函数。
本实施例中优选的，所述步骤S0包括：
S01、基于与所述三维空间对应的平面地质图、剖面地质图，获取地层界面点数据和地层产状数据，然后将地层界面点数据和地层产状数据从二维到三维映射，获取三维空间中的点集数据；
所述点集数据包括地层界面点的位置和地层产状点的位置以及地层产状点的走向、倾向、倾角；
S02、基于与所述三维空间对应的的地层柱状图，获取地层界面点的场值；所述地层界面的场值即为地层界面的相对埋深值。
S03、基于地层产状采样点的走向、倾向、倾角，获取地层产状点的梯度矢量值；
本实施例中的剖面地质图中的地层产状数据(包括走向、倾向和倾角)对地层的形态和分布起着重要的控制作用，地质人员利用一些位置观察到的地层界线数据，结合其它位置上的地层产状测量数据来重构三维地层结构。将三维空间定义为一个标量函数f(p),f是三维空间中任意点p位置上地层的相对埋深值，而把要模拟的一系列地层界面表示为fk(i＝1，...，K)，相当于一系列特定的等深面，即满足位势场f(p)＝fk条件的曲面，地层则占据了其底面fk和顶面fk+1之间的空间。地层内部的每个点都有一个相对于第四系顶面的埋深值，在每个地层内部也存在着无数个互不相交的等势(深)面，地层内部的场值自底向顶是渐变的。
本实施例中，地层产状点的梯度矢量值由地层产状采样点的走向、倾向、倾角转换得到的。在地层位势场构建中，不但f(p)在相同地层边界线上的函数值相等，还要求在地层内部测点上的梯度也等于测量值gj，即场函数满足两个约束条件：f(pi)＝fi，i＝1，2，...，N和
Figure BDA0002678214340000121
因此，需要在另外一些已知梯度数据
Figure BDA0002678214340000122
的控制点qj来约束标量场的方向，这些梯度数据可由地层的产状测量数据转换得到。参见图6，梯度矢量g即地质界面的法向量n方向，指向标量场值增长最快的方向(即指向老的地层)，它与地层走向矢量s和倾角矢量d两两正交；走向θ1是地质界面与水平面的交线的延伸方向，它与走向矢量平行，一般用与正北方向的夹角来表示；倾向θ2为倾角矢量在水平面上的投影，也用与正北方向的夹角来表示，并且走向与倾向相互垂直；倾角θ3是倾角矢量与倾向之间的夹角，地质界面的走向、倾向和倾角产状三要素可以通过测量获取。
本实施例中梯度是一个有模长、有方向的矢量，而梯度的模长在地层产状测量中是难以获取的。本实施例中采用相对埋深作为标量场的属性值，它的确切定义就是相对特定界面沿梯度方向上变化的距离值，因此在地层均匀变化情况下可以假定
Figure BDA0002678214340000123
显然，如果采用地质年代或序号作为标量场的属性值，即场量和距离是不同含义的变量，这些情况下||gj||≠1。
S04、基于三维空间中地层界面点的位置、地层界面点的场值、地层产状点的位置以及地层产状点的梯度矢量值，获取地层界面采样点数据和地层产状采样点数据；
所述地层界面采样点数据为与所述三维空间中任一子空间对应的地层界面点的位置和场值；
所述地层产状采样点数据为与所述三维空间中任一子空间对应的地层产状点的位置和梯度矢量值。
本实施例中优选的，所述步骤S1包括：
基于与所述子空间相应的地层界面采样点数据和地层产状采样点数据，根据公式(1)、公式(2)以及预先设定的条件，确定所述子空间的位势场函数；
其中，公式(1)为：
Figure BDA0002678214340000131
其中，公式(2)为：
Figure BDA0002678214340000132
其中，||p-pi||为任意点p到地层界面采样点pi点的欧式距离；
||p-qj||为任意点p到地层产状采样点qj点的欧式距离；
f(p)为地层位势场函数；
Figure BDA0002678214340000135
地层位势场的梯度场函数；
N为地层界面采样点个数；
N为地层界面采样点个数；
M为地层产状采样点个数；
Figure BDA0002678214340000134
为预先设定的径向基函数；
其中，本实施例中用的径向基函数为Cubic函数：
Figure BDA0002678214340000133
Figure BDA0002678214340000141
r代表两点之间的欧式距离，exp()为以自然常数e为底的指数函数，β为实数常数，k和d为整数常数。
Figure BDA0002678214340000144
为Hamilton算子，且
Figure BDA0002678214340000142
Figure BDA0002678214340000145
为求偏导的运算符；
Figure BDA0002678214340000146
为对x方向求偏导的运算符；其中，x方向为纬度的方向；
Figure BDA0002678214340000147
为对y方向求偏导的运算符；其中，y方向为经度的方向；
Figure BDA0002678214340000148
为对z方向求偏导的运算符；其中，z方向为预先设定的垂直于x和y的方向；
H为Hessian算子，且
Figure BDA0002678214340000143
Figure BDA0002678214340000149
为求二次偏导的运算符；
Figure BDA00026782143400001410
为对x方向求二次偏导的运算符；
Figure BDA00026782143400001411
为对y方向求二次偏导的运算符；
Figure BDA00026782143400001412
为对z方向求二次偏导的运算符；
<，>为两个矢量的内积运算符；
αi为地层界面采样点的系数；
βj为地层产状采样点的矢量系数；
C(p)＝c1+c2px+c3py+c4pz为预先设定的一次多项式；
本实施例中添加的预先设定的一次多项式可以用来保证等值曲面光滑性和连续性。
px为任意点p在以预先设定的xyz-o坐标系中x轴上的坐标；
py为任意点p在以预先设定的xyz-o坐标系中y轴上的坐标；
pz为任意点p在以预先设定的xyz-o坐标系中z轴上的坐标；
其中，预先设定的xyz-o坐标系以预先设定的点为原点，以纬度的方向为x轴的方向，以经度的方向为y轴的方向，以预先设定的垂直于x和y的方向为z轴方向；
c1为预先设定的截距系数；
c2为预先设定的px的系数；
c3为预先设定的py的系数；
c4为预先设定的pz的系数；
所述预先设定的条件为f(p)的二阶导数的函数值最小。
本实施例中能量函数(a)是f(p)二阶导数的函数，反映了函数的凹凸性和光滑程度，为使隐式函数f(p)尽可能光滑，就要使能量函数最小。
Figure BDA0002678214340000151
其中
Figure BDA0002678214340000152
是隐式函数f(p)的二阶偏导数，可利用变分技术在约束条件下求解能量最小的问题。
本实施例中优选的，所述步骤S1具体包括：
S11、基于所述地层界面采样点数据和地层产状采样点数据，确定公式(1)和公式(2)中的参数系数αi、βj及c1、c2、c3、c4的具体值；
S12、基于所述参数系数αi、βj及c1、c2、c3、c4的具体值，确定所述预先设定区域的地层位势场函数。
本实施例中优选的，所述步骤S11包括：
S111、所述地层界面采样点数据和地层产状采样点数据代入式(1)和公式(2)，可得：
Figure BDA0002678214340000161
Figure BDA0002678214340000162
pk为地层界面采样点的位置；fk为地层界面采样点的场值；与所有地层界面采样点pi和所有地层产状采样点qj构成方程(3)。
qk为地层产状采样点的位置；gk为地层界面采样点的场值；与所有地层界面采样点pi和所有地层产状采样点qj构成方程(4)。
S112、根据预先设定的条件也就是使公式(a)能量函数最小，确定正交条件：
Figure BDA0002678214340000163
Figure BDA0002678214340000164
Figure BDA0002678214340000166
为第j个地层产状采样点的矢量系数βj在x方向上的分量；
Figure BDA0002678214340000167
为第j个地层产状采样点的矢量系数βj在y方向上的分量；
Figure BDA0002678214340000168
为第j个地层产状采样点的矢量系数βj在z方向上的分量；
S113、基于所述公式(3)、(4)、(5)、(6)，获取公式(7)；
Figure BDA0002678214340000165
其中，Φ为N×N的矩阵，且Φ的元素为
Figure BDA0002678214340000169
且
Figure BDA0002678214340000171
Figure BDA0002678214340000175
为N×3M的矩阵，且
Figure BDA0002678214340000176
的元素为
Figure BDA0002678214340000177
且
Figure BDA0002678214340000178
HΦ为3M×3M的矩阵，且HΦ的元素为
Figure BDA0002678214340000179
且
Figure BDA00026782143400001710
Figure BDA0002678214340000172
其中元素
Figure BDA00026782143400001711
元素
Figure BDA00026782143400001712
元素
Figure BDA00026782143400001713
地层界面采样点的场值数据f＝[f1 f2…fN]T，地层界面采样点的场值数据g＝[g1g2 … gM]T；
S114、基于公式(7)，确定系数αi、βj及c1、c2、c3、c4的具体值。
本实施例中优选的，所述步骤S12具体包括：
将所述确定的系数αi、βj及c1、c2、c3、c4的具体值，代入到公式(1)和(2)中，确定所述预先设定子空间的地层位势场函数；
所述预先设定子空间的位势场函数为：
Figure BDA0002678214340000173
Figure BDA0002678214340000174
其中，A为系数αi的具体值；B为系数βj的具体值；D(p)＝d1+d2px+d3py+d4pz；d1为c1的具体值；d2为c2的具体值；d3为c3的具体值；d4为c4的具体值。
本实施例中优选的，所述步骤S2包括：
S21、根据所述子空间的地层位势场函数，得到三维网格点的地层位势场的场值及其梯度矢量值；
所述三维网格点为在预先设定的xyz-o坐标系中沿x、y、z轴方向，在三维空间中按预先设定的分辨率间隔Δx、Δy、Δz规则采样获得的位置点；
S22、基于所述三维网格点的地层位势场的场值及其梯度矢量值，采用预先设定的数字高程模型DEM、和预先确定的子空间边界获取预先设定的场值的等势面；
S23、参见图8、采用所述数字高程模型DEM将所述地层位势场值的等势面、预先确定的子空间边界面联合围成体，获取三维地质实体模型。
参见图9，本实施例中还可以将三维地质实体模型分割成规则排列的立方体的集合，每个立方体有对应的地质体类型和地层位势场值，得到三维块体模型。
本实施例中的一种获取三维地质实体模型的方法，由于将地层界面的相对埋深作为位势场的场值，地层产状数据作为位势场的梯度约束，通过地层界面点数据和地层产状数据构建地层位势场函数，并基于所述预先设定空间的地层位势场函数，建立所述预先设定的空间的三维地质实体模型，相对于现有技术而言，其可以更准确的表示现实中的三维地质。
举例说明，采用本实施例中的一种获取三维地质实体模型的方法获取位于广西壮族自治区德保县凌念-那查地区(如图3所示)的三维地质实体模型。参见图3，其中地区内主要地层包含泥盆系(那高岭组D1n、郁江阶D1y、东岗岭阶D2d和上泥盆统D3)、石炭系(下石炭统C1、中石炭统C2和上石炭统C3)、二叠系(栖霞组P1q和茅口组P1m)、三叠系(马脚岭组T1m、北泗组T1b、百逢组下段T2b1和百逢组上段T2b2)。根据广西壮族自治区德保县凌念-那查地区的地层柱状图(如图4所示)中各地层的厚度范围，确定了各地质界面相的相对埋深值，该相对埋深被作为隐式函数插值所使用的场值。从平面地质图和剖面地质图(如图5所示)中提取的各地层界面和断层面的地层界面采样点数据和地层产状采样点数据。然后使用HRBFs方法分别构建断层两侧地层位势场(参见图7)、地层面三维实体模型(参见图8)和地层三维块体模型(参见图9)。
另一方面，本实施例提供一种获取三维地质实体模型的的系统，包括：
至少一个处理器；以及
与所述处理器通信连接的至少一个存储器，其中，所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行上述任一获取三维地质实体模型的方法。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例，或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
应当注意的是，在权利要求中，不应将位于括号之间的任何附图标记理解成对权利要求的限制。词语“包含”不排除存在未列在权利要求中的部件或步骤。位于部件之前的词语“一”或“一个”不排除存在多个这样的部件。本发明可以借助于包括有若干不同部件的硬件以及借助于适当编程的计算机来实现。在列举了若干装置的权利要求中，这些装置中的若干个可以是通过同一个硬件来具体体现。词语第一、第二、第三等的使用，仅是为了表述方便，而不表示任何顺序。可将这些词语理解为部件名称的一部分。
此外，需要说明的是，在本说明书的描述中，术语“一个实施例”、“一些实施例”、“实施例”、“示例”、“具体示例”或“一些示例”等的描述，是指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
尽管已描述了本发明的优选实施例，但本领域的技术人员在得知了基本创造性概念后，则可对这些实施例作出另外的变更和修改。所以，权利要求应该解释为包括优选实施例以及落入本发明范围的所有变更和修改。
显然，本领域的技术人员可以对本发明进行各种修改和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也应该包含这些修改和变型在内。

Claims (6)

1.一种获取三维地质实体模型的方法，其特征在于，包括：
S1、针对与预先设定的空间相应的地层界面采样点数据和地层产状采样点数据，获取所述预先设定的空间的地层位势场函数；
所述地层界面采样点数据包括：地层界面采样点的位置、地层界面采样点的场值；
所述地层产状采样点数据包括：地层产状采样点的位置、地层产状采样点的梯度矢量值；
S2、基于所述预先设定空间的地层位势场函数，建立所述预先设定的空间的三维地质实体模型；
所述步骤S1之前还包括：
S001、针对预先设定的三维空间和与所述三维空间所对应的平面地质图、剖面地质图和地层柱状图，将所述三维空间划分为多个子空间；
S002、基于与所述三维空间对应的平面地质图和剖面地质图，获取子空间对应的地层界面采样点数据和地层产状采样点数据；
相应的，所述S1为：
基于与所述子空间相应的地层界面采样点数据和地层产状采样点数据，获取所述子空间的地层位势场函数；
所述步骤S002包括：
S0021、基于与所述三维空间对应的平面地质图和剖面地质图，获取地层界面点数据和地层产状数据，然后将地层界面点数据和地层产状数据从二维到三维映射，获取三维空间中的点集数据；
所述点集数据包括地层界面点的位置和地层产状点的位置以及地层产状点的走向、倾向和倾角；
S0022、基于与所述三维空间对应的地层柱状图，获取地层界面点的场值；所述地层界面点的场值为所述地层柱状图中的地层相对埋深值；
S0023、基于地层产状点的走向、倾向和倾角，获取地层产状点的梯度矢量值；
S0024、基于三维空间中地层界面点的位置、地层界面点的场值、地层产状点的位置以及地层产状点的梯度矢量值，获取地层界面采样点数据和地层产状采样点数据；
所述地层界面采样点数据为与所述三维空间中任一子空间所对应的地层界面点的位置和场值；
所述地层产状采样点数据为与所述三维空间中任一子空间所对应的地层产状点的位置和梯度矢量值；
所述步骤S1包括：
基于与所述子空间相应的地层界面采样点数据和地层产状采样点数据，根据预先设定的条件，采用公式(1)、公式(2)确定所述子空间的地层位势场函数；
其中，公式(1)为：
Figure FDA0003644306240000021
其中，公式(2)为：
Figure FDA0003644306240000022
其中，||p-pi||为子空间中任意点p到地层界面采样点pi点的欧式距离；||p-qj||为子空间中任意点p到地层产状采样点qj点的欧式距离；
f(p)为地层位势场函数；
Figure FDA0003644306240000023
地层位势场的梯度场函数；
N为地层界面采样点个数；
M为地层产状采样点个数；
Figure FDA0003644306240000024
为预先设定的径向基函数；
Figure FDA0003644306240000031
为Hamilton算子，且
Figure FDA0003644306240000032
Figure FDA0003644306240000033
为求偏导的运算符；
Figure FDA0003644306240000034
为对x方向求偏导的运算符；其中，x方向为纬度的方向；
Figure FDA0003644306240000035
为对y方向求偏导的运算符；其中，y方向为经度的方向；
Figure FDA0003644306240000036
为对z方向求偏导的运算符；其中，z方向为预先设定的垂直于x和y的方向；
H为Hessian算子，且
Figure FDA0003644306240000037
Figure FDA0003644306240000038
为求二次偏导的运算符；
Figure FDA0003644306240000039
为对x方向求二次偏导的运算符；
Figure FDA00036443062400000310
为对y方向求二次偏导的运算符；
Figure FDA00036443062400000311
为对z方向求二次偏导的运算符；
<，>为两个矢量的内积运算符；
αi为地层界面采样点的系数；
βj为地层产状采样点的矢量系数；
C(p)＝c1+c2px+c3py+c4pz为预先设定的一次多项式；
px为任意点p在以预先设定的xyz-o坐标系中x轴上的坐标；
py为任意点p在以预先设定的xyz-o坐标系中y轴上的坐标；
pz为任意点p在以预先设定的xyz-o坐标系中z轴上的坐标；
其中，预先设定的xyz-o坐标系以预先设定的点为原点，以纬度的方向为x轴的方向，以经度的方向为y轴的方向，以预先设定的垂直于x和y的方向为z轴方向；
c1为预先设定的截距系数；
c2为预先设定的px的系数；
c3为预先设定的py的系数；
c4为预先设定的pz的系数；
所述预先设定的条件为f(p)的二阶导数的函数值最小。
2.根据权利要求1所述的方法，其特征在于，所述步骤S1具体包括：
S11、基于所述地层界面采样点数据和地层产状采样点数据，确定公式(1)和公式(2)中的参数系数αi、βj及c1、c2、c3、c4的具体值；
S12、基于所述参数系数αi、βj及c1、c2、c3、c4的具体值，确定所述预先设定子空间的地层位势场函数。
3.根据权利要求2所述的方法，其特征在于，所述步骤S11包括：
S111、所述地层界面采样点数据和地层产状采样点数据作为任意点p及其对应的地层位势场值及其梯度场值分别代入式(1)和公式(2)，可得：
Figure FDA0003644306240000041
Figure FDA0003644306240000042
pk为地层界面采样点的位置；fk为地层界面采样点的场值；与所有地层界面采样点pi和所有地层产状采样点qj构成公式(3)；
qk为地层产状采样点的位置；gk为地层界面采样点的场值；与所有地层界面采样点pi和所有地层产状采样点qj构成公式(4)；
S112、根据预先设定的条件，确定正交条件：
Figure FDA0003644306240000051
Figure FDA0003644306240000052
Figure FDA0003644306240000053
为第j个地层产状采样点的矢量系数βj在x方向上的分量；
Figure FDA0003644306240000054
为第j个地层产状采样点的矢量系数βj在y方向上的分量；
Figure FDA0003644306240000055
为第j个地层产状采样点的矢量系数βj在z方向上的分量；
S113、基于所述公式(3)、公式(4)、公式(5)、公式(6)，获取公式(7)；
Figure FDA0003644306240000056
其中，Φ为N×N的矩阵，且Φ的元素为
Figure FDA0003644306240000057
且
Figure FDA0003644306240000058
Figure FDA0003644306240000059
为N×3M的矩阵，且
Figure FDA00036443062400000510
的元素为
Figure FDA00036443062400000511
且
Figure FDA00036443062400000512
HΦ为3M×3M的矩阵，且HΦ的元素为
Figure FDA00036443062400000513
且
Figure FDA00036443062400000514
Figure FDA00036443062400000515
其中元素
Figure FDA00036443062400000516
元素
Figure FDA00036443062400000517
元素
Figure FDA00036443062400000518
地层界面采样点的场值数据f＝[f1 f2 … fN]T，地层产状采样点的场值梯度矢量数据g＝[g1 g2 … gM]T；
S114、基于公式(7)，确定系数αi、βj及c1、c2、c3、c4的具体值。
4.根据权利要求3所述的方法，其特征在于，所述步骤S12具体包括：
将所述确定的系数αi、βj及c1、c2、c3、c4的具体值，代入到公式(1)和(2)中，确定所述预先设定子空间的地层位势场函数；
所述预先设定子空间的地层位势场函数为：
Figure FDA0003644306240000061
Figure FDA0003644306240000062
其中，A为系数αi的具体值；B为系数βj的具体值；D(p)＝d1+d2px+d3py+d4pz；d1为c1的具体值；d2为c2的具体值；d3为c3的具体值；d4为c4的具体值。
5.根据权利要求4所述的方法，其特征在于，所述步骤S2包括：
S21、根据所述子空间的地层位势场函数，得到三维网格点的地层位势场的场值及其梯度矢量值；
所述三维网格点为在预先设定的xyz-o坐标系中沿x、y、z轴方向，在三维空间中按预先设定的分辨率间隔Δx、Δy、Δz规则采样获得的位置点；
S22、基于所述三维网格点的地层位势场的场值及其梯度矢量值，采用预先设定的数字高程模型DEM和预先确定的子空间边界获取预先设定的地层位势场值的等势面；
S23、采用所述数字高程模型DEM将所述地层位势场值的等势面、预先确定的子空间边界面联合围成体，获取三维地质实体模型。
6.一种获取三维地质实体模型的系统，其特征在于，包括：
至少一个处理器；以及
与所述处理器通信连接的至少一个存储器，其中，所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行如权利要求1至5任一所述的获取三维地质实体模型的方法。
CN202010954732.XA 2020-09-11 2020-09-11 一种获取三维地质实体模型的方法及系统 Active CN112116708B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010954732.XA CN112116708B (zh) 2020-09-11 2020-09-11 一种获取三维地质实体模型的方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010954732.XA CN112116708B (zh) 2020-09-11 2020-09-11 一种获取三维地质实体模型的方法及系统

Publications (2)

Publication Number Publication Date
CN112116708A CN112116708A (zh) 2020-12-22
CN112116708B true CN112116708B (zh) 2022-06-21

Family

ID=73802458

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010954732.XA Active CN112116708B (zh) 2020-09-11 2020-09-11 一种获取三维地质实体模型的方法及系统

Country Status (1)

Country Link
CN (1) CN112116708B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113159321B (zh) * 2021-04-07 2022-05-20 中南大学 重力约束下断裂面形态的贝叶斯推断方法
CN115619983B (zh) * 2022-12-02 2023-04-07 中南大学 一种基于径向基函数的平均曲率可控隐式曲面生成方法
CN116486022A (zh) * 2023-03-23 2023-07-25 北京冽泉环保科技有限公司 一种三维地质模型的构建方法
CN116645484B (zh) * 2023-07-26 2023-10-03 航天宏图信息技术股份有限公司 地质曲面模型的构建方法、装置、电子设备及存储介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102254349A (zh) * 2011-06-30 2011-11-23 华东师范大学 一种使用钻孔数据构建沉积地层系统三维实体模型的方法
CN109389515A (zh) * 2018-10-11 2019-02-26 中石化石油工程技术服务有限公司 一种根据实钻地层界面埋深计算地层产状的方法及系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR3027944A1 (fr) * 2014-10-29 2016-05-06 Services Petroliers Schlumberger Generation d'elements structurels pour formation souterraine utilisant une fonction implicite stratigraphique

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102254349A (zh) * 2011-06-30 2011-11-23 华东师范大学 一种使用钻孔数据构建沉积地层系统三维实体模型的方法
CN109389515A (zh) * 2018-10-11 2019-02-26 中石化石油工程技术服务有限公司 一种根据实钻地层界面埋深计算地层产状的方法及系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
随钻电磁波测井地层产状因素响应特征三维数值模拟;邵才瑞 等;《地球物理学进展》;20170228;第32卷(第1期);第236-242页 *

Also Published As

Publication number Publication date
CN112116708A (zh) 2020-12-22

Similar Documents

Publication Publication Date Title
CN112116708B (zh) 一种获取三维地质实体模型的方法及系统
CN112381937B (zh) 一种基于钻孔和复杂地质剖面的多源地质数据耦合建模方法
Philip et al. A precise method for determining contoured surfaces
RU2306607C2 (ru) Способ, устройство и программный продукт для трехмерного моделирования геологического объема при помощи выбора трехмерных параметров геологической области
Natali et al. Modeling Terrains and Subsurface Geology.
Jones Data structures for three-dimensional spatial information systems in geology
Yilmaz The effect of interpolation methods in surface definition: an experimental study
CN110400371B (zh) 一种水平构造地貌实体的三维模型构建方法
EP1917542A1 (en) Method for tomographic inversion by matrix transformation
Molnár et al. Can the First Military Survey maps of the Habsburg Empire (1763–1790) be georeferenced by an accuracy of 200 meters
NO340138B1 (no) Fremgangsmåte, program og datamaskinsystem for modellering av paleoelvekanaler
CN109636912A (zh) 应用于三维声呐图像重构的四面体剖分有限元插值方法
CN114943178A (zh) 一种三维地质模型建模方法、装置及计算机设备
CN107221028B (zh) 一种基于地震解释数据的地质体闭合曲面三维重建方法
NO335056B1 (no) Metode for beregning av 3D gittermodeller av et reservoar
Dhont et al. 3-D modeling of geologic maps from surface data
CN106875484A (zh) 一种基于三维地形的地质堆积体快速拟合建模方法
CN108828669B (zh) 一种二维相交测线静校正处理方法、装置及系统
CN112346139B (zh) 一种重力数据多层等效源延拓与数据转换方法
CN104240301A (zh) 地质曲面重构方法及设备
CN108121853A (zh) 一种基于AutoCAD计算挖填方工程量的系统与方法
Holzrichter Processing and interpretation of satellite and ground based gravity data at different lithospheric scales
CN107730586B (zh) 一种地层建模的方法和系统
CN111951394A (zh) 基于地质图的断层构造单元三维模型构建方法及装置
Rózsa et al. Prediction of vertical gravity gradients using gravity and elevation data

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant