CN104240301B - 地质曲面重构方法及设备 - Google Patents
地质曲面重构方法及设备 Download PDFInfo
- Publication number
- CN104240301B CN104240301B CN201410461559.4A CN201410461559A CN104240301B CN 104240301 B CN104240301 B CN 104240301B CN 201410461559 A CN201410461559 A CN 201410461559A CN 104240301 B CN104240301 B CN 104240301B
- Authority
- CN
- China
- Prior art keywords
- mesh point
- property value
- unknown
- weight
- mesh
- 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
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
提供一种地质曲面重构方法及设备。所述方法包括：a)将基于采集到的地质数据生成的三维的空间散点沿Z轴方向投影到二维平面，以获得各空间散点在二维平面上的投影点；b)在由所有投影点构成的区域内建立N重矩形网格拓扑；c)根据各投影点的属性值对N重矩形网格拓扑的至少一个网格点进行插值计算，以获得至少一个网格点的属性值；d)根据已知网格点的属性值计算未知网格点的初始属性值；e)根据已知网格点的属性值和未知网格点的初始属性值，通过有限差分求解格式对几何偏微分方程进行迭代求解，以获得未知网格点的属性值；f)获得由N重矩形网格拓扑的所有网格点及其属性值表示的地质曲面。
Description
技术领域
本发明涉及地球物理勘探技术领域，更具体地讲，涉及一种地质曲面重构方法及设备。
背景技术
曲面可视化重构是油气等矿产资源评价、地质构造和各种数据场分布特征研究等最为直接的研究手段。随着各种三维空间数据获取技术的不断发展，特别是地球物理探测技术的发展使包含被测物体更多细节的地质数据的获取成为可能，地质勘探结果大多用一些离散的、不规则分布的数据点来反映，基于这些地质数据的三维地质建模技术已成为综合地质解释工作的发展方向。
为了通过这些离散的地质数据建立起区域性连续的整体模型，需要利用插值、逼近等曲面构造方法。曲面插值(Surface Interpolation)方法是严格通过给定的数据点来构造曲面，并根据原始数据点的值来插补空白区的值而不改变原始数据点的值。而曲面逼近(Surface Approximation)方法则是利用相对简单的数学曲面来近似地构造复杂的地质曲面，根据一定的数学准则使所给出的数学曲面最大限度地逼近地质曲面。运用较多的插值和拟合方法包括：按近点距离加权平均法、按方位取点加权法、双线性插值法、移动曲面插值法、二元三点插值法、三角剖分方法、克里金(Kriging)插值法和三次样条函数拟合法、趋势面拟合法、加权最小二乘拟合法等方法。
但由于地质数据包含层位数据、断层数据等，数据量较大，且由于受到断层约束，边界条件较为复杂，目前所使用的曲面重构方式的重构过程较为复杂(例如，需要进行复杂的参数计算)、重构时间长，并且得到的地质曲面往往不能令人满意(例如，需要对得到的地质曲面进行平滑处理)。
发明内容
本发明的示例性实施例在于提供一种地质曲面重构方法及设备，其能够得到更真实、更高质量的光滑地质曲面，且重构过程简便、能够实现快速重构。
根据本发明的一方面，提供一种地质曲面重构方法，包括：a)将基于采集到的地质数据生成的三维的空间散点沿Z轴方向投影到二维平面，以获得各空间散点在二维平面上的投影点，其中，各空间散点的投影点的属性值为所述空间散点沿Z轴方向的坐标的值；b)在由所有投影点构成的区域内建立N重矩形网格拓扑，其中，N为大于0的整数，且所述区域越大，N越大；c)根据各投影点的属性值对N重矩形网格拓扑的至少一个网格点进行插值计算，以获得至少一个网格点的属性值，其中，所述至少一个网格点作为已知网格点，未获得属性值的网格点作为未知网格点；d)根据已知网格点的属性值计算未知网格点的初始属性值；e)根据已知网格点的属性值和未知网格点的初始属性值，通过有限差分求解格式对几何偏微分方程进行迭代求解，以获得未知网格点的属性值；f)获得由N重矩形网格拓扑的所有网格点及其属性值表示的地质曲面。
可选地，步骤c)包括：按照从一重网格点至N重网格点的顺序，依次对N重矩形网格拓扑的每重网格点进行插值计算，其中，当对i重网格点进行插值计算时，如果所述i重网格点所在的i重网格中存在投影点，则根据存在的投影点的属性值对所述i重网格点进行插值计算，其中，i为大于0小于等于N的整数。
可选地，在步骤c)中，通过加权距离反比插值法或克里金插值法根据存在的投影点的属性值对所述i重网格点进行插值计算。
可选地，步骤d)包括：D1)根据已知网格点的属性值对至少一个未知网格点进行插值计算，以获得所述至少一个未知网格点的初始属性值；D2)针对未通过步骤D1)进行插值计算的未知网格点，将距离所述未知网格点最近的网格点的属性值或初始属性值的平均值作为所述未知网格点的初始属性值。
可选地，步骤D1)包括：D11)根据断层多边形约束将所述由所有投影点构成的区域划分为多个子区域；D12)对至少一个子区域内的未知网格点进行插值计算，以获得所述至少一个子区域内的未知网格点的初始属性值。
可选地，在步骤D12)中：对所述子区域内具有3个以上已知网格点的行进行B样条曲线插值，以获得所述行上的未知网格点的初始属性值；对所述子区域内具有3个以上已知网格点的列进行B样条曲线插值，以获得所述列上的未知网格点的初始属性值，其中，如果对所述子区域内的未知网格点所在的行和列分别进行了B样条曲线插值，则将所述未知网格点的初始属性值更新为分别获得的所述未知网格点的初始属性值的平均值。
可选地，N＝2，所述几何偏微分方程为四阶几何偏微分方程，步骤e)包括：E1)根据一重已知网格点的属性值和一重未知网格点的初始属性值，在第一预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得一重未知网格点的属性值；E2)分别针对各一重网格中的二重未知网格点，根据所述一重网格中的二重已知网格点的属性值和二重未知网格点的初始属性值，在第二预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得所述一重网格中的二重未知网格点的属性值；E3)根据所有已知网格点的属性值和通过步骤E1)和E2)获得的未知网格点的属性值，在第三预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以更新所有未知网格点的属性值。
可选地，所述第一预设条件为迭代次数达到第一预设迭代次数，或本次迭代求解得到的各一重未知网格点的属性值与上次迭代求解得到的各一重未知网格点的属性值的平均误差小于等于第一预设误差阈值，所述第二预设条件为迭代次数达到第二预设迭代次数，或本次迭代求解得到的各二重未知网格点的属性值与上次迭代求解得到的各二重未知网格点的属性值的平均误差小于等于第二预设误差阈值，所述第三预设条件为迭代次数达到第三预设迭代次数，或本次迭代求解得到的各未知网格点的属性值与上次迭代求解得到的各未知网格点的属性值的平均误差小于等于第三预设误差阈值。
根据本发明的另一方面，提供一种地质曲面重构设备，包括：投影装置，将基于采集到的地质数据生成的三维的空间散点沿Z轴方向投影到二维平面，以获得各空间散点在二维平面上的投影点，其中，各空间散点的投影点的属性值为所述空间散点沿Z轴方向的坐标的值；网格拓扑建立装置，在由所有投影点构成的区域内建立N重矩形网格拓扑，其中，N为大于0的整数，且所述区域越大，N越大；网格点插值装置，根据各投影点的属性值对N重矩形网格拓扑的至少一个网格点进行插值计算，以获得至少一个网格点的属性值，其中，所述至少一个网格点作为已知网格点，未获得属性值的网格点作为未知网格点；初始属性值计算装置，根据已知网格点的属性值计算未知网格点的初始属性值；迭代装置，根据已知网格点的属性值和未知网格点的初始属性值，通过有限差分求解格式对几何偏微分方程进行迭代求解，以获得未知网格点的属性值；地质曲面表示装置，获得由N重矩形网格拓扑的所有网格点及其属性值表示的地质曲面。
可选地，N＝2，所述几何偏微分方程为四阶几何偏微分方程，迭代装置包括：一重网格点迭代单元，根据一重已知网格点的属性值和一重未知网格点的初始属性值，在第一预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得一重未知网格点的属性值；二重网格点迭代单元，分别针对各一重网格中的二重未知网格点，根据所述一重网格中的二重已知网格点的属性值和二重未知网格点的初始属性值，在第二预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得所述一重网格中的二重未知网格点的属性值；全局迭代单元，根据所有已知网格点的属性值和通过一重网格点迭代单元和二重网格点迭代单元获得的未知网格点的属性值，在第三预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以更新所有未知网格点的属性值。
根据本发明示例性实施例的地质曲面重构方法及设备，可以得到更真实、更高质量的光滑地质曲面，且重构过程简便、能够实现快速重构。
将在接下来的描述中部分阐述本发明总体构思另外的方面和/或优点，还有一部分通过描述将是清楚的，或者可以经过本发明总体构思的实施而得知。
附图说明
通过下面结合示例性地示出实施例的附图进行的描述，本发明示例性实施例的上述和其他目的和特点将会变得更加清楚，其中：
图1示出根据本发明示例性实施例的地质曲面重构方法的流程图；
图2示出根据本发明示例性实施例的二重矩形网格拓扑的示例；
图3示出根据本发明示例性实施例的根据已知网格点的属性值对至少一个未知网格点进行插值计算的方法的流程图；
图4示出根据本发明示例性实施例的通过有限差分求解格式对几何偏微分方程进行迭代求解来获得未知网格点的属性值的方法的流程图；
图5示出根据本发明示例性实施例的地质曲面重构设备的框图；
图6示出根据本发明示例性实施例的插值单元的框图；
图7示出根据本发明示例性实施例的迭代装置的框图。
具体实施方式
现将详细参照本发明的实施例，所述实施例的示例在附图中示出，其中，相同的标号始终指的是相同的部件。以下将通过参照附图来说明所述实施例，以便解释本发明。
图1示出根据本发明示例性实施例的地质曲面重构方法的流程图。
如图1所示，在步骤S100，将基于采集到的地质数据生成的三维的空间散点沿Z轴方向投影到二维平面，以获得各空间散点在二维平面上的投影点。各空间散点的投影点的属性值为该空间散点沿Z轴方向的坐标的值。
具体说来，将三维的各空间散点(x，y，z)沿Z轴方向投影到二维平面，以获得该空间散点的投影点(u，v)，且获得的投影点的属性值为该空间散点的z值，应该理解，u＝x，v＝y。
在步骤S200，在由所有投影点构成的区域内建立N重矩形网格拓扑。
N为大于0的整数，且所述区域越大，N越大。例如，如果所有投影点构成的区域较小，则可只建立一重或二重矩形网格拓扑，如果所有投影点构成的区域较大，则可建立二重以上的矩形网格拓扑。
以二重矩形网格拓扑为例，可先在所有投影点构成的区域内建立一重网格拓扑，以通过一重网格控制整体曲面的形态和趋势，然后在各一重网格内建立二重网格，以将各一重网格划分为更小的区域，通过分区域构造曲面，使构造的曲面更精细和准确。也因此，N越大，网格越密，构造的曲面越精细和准确。图2示出根据本发明示例性实施例的二重矩形网格拓扑的示例，参照图2，由实线围成的4个网格表示一重网格，实心网格点表示一重网格点；在各一重网格内部由虚线围成的网格、或虚线和实线围成的网格表示二重网格，空心网格点表示二重网格点，应该理解，一重网格点也是二重网格点。
关于确定所有投影点构成的区域，作为示例，可根据投影点的坐标来确定所有投影点在二维平面上构成的区域。例如，假设所有投影点的u坐标的最大值为umax，u坐标的最小值为umin，则沿u轴方向的边界为u＝umin-C，u＝umax+C，C为非负常数，同理，假设所有投影点的v坐标的最大值为vmax，v坐标的最小值为vmin，沿v轴方向的边界为v＝vmin-C、v＝vmax+C，通过上述边界即可确定所有投影点构成的区域。
在步骤S300，根据各投影点的属性值对N重矩形网格拓扑的至少一个网格点进行插值计算，以获得至少一个网格点的属性值。至少一个网格点作为已知网格点，未获得属性值的网格点作为未知网格点。
作为示例，可按照从一重网格点至N重网格点的顺序，依次对N重矩形网格拓扑的每重网格点进行插值计算。当对i重网格点进行插值计算时，如果该i重网格点所在的i重网格中存在投影点，则根据存在的投影点的属性值对该i重网格点进行插值计算。i为大于0小于等于N的整数。
如果该i重网格点所在的i重网格中不存在投影点，则不对该i重网格点进行插值计算。
作为示例，可通过加权距离反比插值法或克里金插值法根据存在的投影点的属性值对该i重网格点进行插值计算。也可通过其它适合的方法根据存在的投影点的属性值对该i重网格点进行插值计算。
作为示例，当N＝2，即，N重矩形网格拓扑为二重矩形网格拓扑时，先对一重网格点进行插值计算，然后对二重网格点进行插值计算。应该理解，一重网格点同时也是二重网格点，因此，作为优选方式，如果对某一重网格点进行插值计算得到了一个属性值，然后将其作为二重网格点进行插值计算得到了另一属性值，则由于其作为一重网格点时，对其进行插值计算时所利用的投影点更多，插值精度更高，因此，保留其作为一重网格点时的插值结果。
在步骤S400，根据已知网格点的属性值计算未知网格点的初始属性值。
可使用各种适合的方法根据已知网格点的属性值计算未知网格点的初始属性值。
作为示例，可首先根据已知网格点的属性值对至少一个未知网格点进行插值计算，以获得所述至少一个未知网格点的初始属性值。然后，针对未通过上述插值计算的未知网格点，将距离该未知网格点最近的网格点的属性值或初始属性值的平均值作为该未知网格点的初始属性值。应该理解，通过上述步骤所有未知网格点都具有了初始属性值。
在步骤S500，根据已知网格点的属性值和未知网格点的初始属性值，通过有限差分求解格式对几何偏微分方程进行迭代求解，以获得未知网格点的属性值。
通过有限差分求解格式对几何偏微分方程进行迭代求解来获得未知网格点的属性值，避免了复杂的参数计算，大大降低了对操作者的数学专业知识的要求；可适用于含垂直断层、正断层和逆断层等复杂断层约束的曲面重构；引进的参数少，不需要大量的运算，能够实现含多种断层约束的地质曲面的快速重构。
几何偏微分方程可以是二阶几何偏微分方程、四阶几何偏微分方程等。选择不同的几何偏微分方程，可以实现不同特性的曲面的重建。例如，选择二阶几何偏微分方程，生成的曲面可在连接处(边界处)达到C0连续，即，相邻两个曲面具有交点；选择四阶几何偏微分方程，生成的曲面可在连接处达到C1连续，即，相邻两个曲面在交点处具有相同的一阶导数。
作为示例，四阶几何偏微分方程表示如下：
其中，Δs为拉普拉斯算子，t为时间变量，X(u，v)指示坐标为(u，v)的网格点的属性值。
可按照如下方式得到式(1)的有限差分求解格式：
对式(1)进行时间方向的离散：
其中，指示网格点(ui，vj)在tk+1时刻的属性值，指示网格点(ui，vj)在tk时刻的属性值，τ为tk+1时刻与tk时刻的时间差。
对式(1)进行空间方向的离散：
其中，Xi+1，j指示(ui，vj)沿u轴方向的下一个网格点(ui+1，vj)在tk时刻的属性值，Xi，j+1指示(ui，vj)沿v轴方向的下一个网格点(ui，vj+1)在tk时刻的属性值，相应地，Xi+2，j、Xi-1，j、Xi-2，j、Xi，j+2、Xi，j-1、Xi，j-2、Xi+1，j+1、Xi+1，j-1、Xi，j-1、Xi-1，j+1、Xi-1，j-1指示类似的含义，hu为沿u轴方向的相邻两个网格点间的距离，hv为沿v轴方向的相邻两个网格点间的距离。应该理解，如果网格点(ui，vj)为一重网格点，则网格点(ui+1，vj)、(ui，vj+1)等相关网格点均为一重网格点，hu为沿u轴方向的相邻两个一重网格点间的距离，hv为沿v轴方向的相邻两个一重网格点间的距离。如果网格点(ui，vj)为二重网格点时，则网格点(ui+1，vj)、(ui，vj+1)等相关网格点均为二重网格点，hu为沿u轴方向的相邻两个二重网格点间的距离，hv为沿v轴方向的相邻两个二重网格点间的距离。
如上所述，式(2)～(5)为式(1)的有限差分求解格式，因此，网格点(ui，vj)在tk+1时刻的值可由式(2)～(5)联合求得。
可将已知网格点的属性值和未知网格点的初始属性值作为初始值，通过式(2)～(5)对几何偏微分方程进行迭代求解，以获得未知网格点的属性值。
通过有限差分求解格式对四阶几何偏微分方程进行迭代求解来获得未知网格点的属性值，计算过程简单、计算时间短，且基于获得的未知网格点的属性值得到的曲面更为光滑。
在步骤S600，获得由N重矩形网格拓扑的所有网格点及其属性值表示的地质曲面，即，重构的地质曲面。
通过上述方式，可以实现含垂直断层、正断层和逆断层等复杂断层约束的曲面快速重构，且生成的曲面具有自然光滑度，无需再进行平滑处理。
图3示出根据本发明示例性实施例的根据已知网格点的属性值对至少一个未知网格点进行插值计算的方法的流程图。
如图3所示，在步骤S401，根据断层多边形约束将所述由所有投影点构成的区域划分为多个子区域。应该理解，子区域的划分与断层多边形约束有关，与网格无关。
在步骤S402，对至少一个子区域内的未知网格点进行插值计算，以获得所述至少一个子区域内的未知网格点的初始属性值。
作为示例，可分别针对上述子区域，对该子区域内具有3个以上已知网格点的行进行B样条曲线插值，以获得该行上的未知网格点的初始属性值。
并对该子区域内具有3个以上已知网格点的列进行B样条曲线插值，以获得该列上的未知网格点的初始属性值。
如果对该子区域内的未知网格点所在的行和列分别进行了B样条曲线插值，则将该未知网格点的初始属性值更新为分别获得的该未知网格点的初始属性值的平均值。
应该理解，也可使用其它适合的方法对至少一个子区域内的未知网格点进行插值计算。
图4示出根据本发明示例性实施例的通过有限差分求解格式对几何偏微分方程进行迭代求解来获得未知网格点的属性值的方法的流程图。
这里，N＝2，即，N重矩形网格拓扑为二重矩形网格拓扑，几何偏微分方程为四阶几何偏微分方程。
如图4所述，在步骤S501，根据一重已知网格点的属性值和一重未知网格点的初始属性值，在第一预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得一重未知网格点的属性值。
一重已知网格点为一重网格点中的已知网格点，一重未知网格点为一重网格点中的未知网格点。
第一预设条件可为迭代次数达到第一预设迭代次数，或本次迭代求解得到的各一重未知网格点的属性值与上次迭代求解得到的各一重未知网格点的属性值的平均误差小于等于第一预设误差阈值。
具体说来，将一重已知网格点的属性值和一重未知网格点的初始属性值作为迭代的初始值，通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，当满足第一预设条件时，认为达到稳态，迭代停止。
在步骤S502，分别针对各一重网格中的二重未知网格点，根据该一重网格中的二重已知网格点的属性值和二重未知网格点的初始属性值，在第二预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得所述一重网格中的二重未知网格点的属性值。
二重已知网格点为二重网格点中的已知网格点，二重未知网格点为二重网格点中的未知网格点。
第二预设条件可为迭代次数达到第二预设迭代次数，或本次迭代求解得到的各二重未知网格点的属性值与上次迭代求解得到的各二重未知网格点的属性值的平均误差小于等于第二预设误差阈值。
具体说来，分别针对每个一重网格，将该一重网格中的二重已知网格点的属性值和二重未知网格点的初始属性值作为迭代的初始值，通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，当满足第二预设条件时，认为达到稳态，迭代停止。
应该理解，第二预设误差阈值可小于第一预设误差阈值。
在步骤S503，根据所有已知网格点的属性值和通过步骤S501和S502获得的未知网格点的属性值，在第三预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以更新所有未知网格点的属性值。
所有已知网格点包括一重已知网格点和二重已知网格点，所有未知网格点包括一重未知网格点和二重未知网格点。
第三预设条件可为迭代次数达到第三预设迭代次数，或本次迭代求解得到的各未知网格点的属性值与上次迭代求解得到的各未知网格点的属性值的平均误差小于等于第三预设误差阈值。
具体说来，将所有已知网格点的属性值和通过步骤S501和S502获得的未知网格点的属性值作为迭代的初始值，通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，当满足第三预设条件时，认为达到稳态，迭代停止。
分别迭代求解一重未知网格点和二重未知网格点的属性值，可以加快几何偏微分方程的收敛速度，实现曲面的快速重构。由于迭代求解各一重网格中的二重未知网格点时，一重网格作为边界不能跨越，可能使得一重网格两侧的二重网格点的属性值相差过大而导致最终生成的曲面不够平滑，因此通过步骤S503进行全局迭代，可实现整个曲面在一重网格上实现光滑拼接，使重构的曲面达到期望的光滑度。
图5示出根据本发明示例性实施例的地质曲面重构设备的框图。
如图5所示，根据本发明示例性实施例的地质曲面重构设备包括：投影装置100、网格拓扑建立装置200、网格点插值装置300、初始属性值计算装置400、迭代装置500和地质曲面表示装置600。这些装置可由专门的传感器等器件来实现，或可由数字信号处理器、现场可编程门阵列等通用硬件处理器来实现，也可通过专用芯片等专用硬件处理器来实现，还可完全通过计算机程序来以软件方式实现。
具体说来，投影装置100用于将基于采集到的地质数据生成的三维的空间散点沿Z轴方向投影到二维平面，以获得各空间散点在二维平面上的投影点。各空间散点的投影点的属性值为该空间散点沿Z轴方向的坐标的值。
具体说来，投影装置100将三维的各空间散点(x，y，z)沿Z轴方向投影到二维平面，以获得该空间散点的投影点(u，v)，且获得的投影点的属性值为该空间散点的z值，应该理解，u＝x，v＝y。
网格拓扑建立装置200用于在由所有投影点构成的区域内建立N重矩形网格拓扑。
N为大于0的整数，且所述区域越大，N越大。例如，如果所有投影点构成的区域较小，则网格拓扑建立装置200可只建立一重或二重矩形网格拓扑，如果所有投影点构成的区域较大，则网格拓扑建立装置200可建立二重以上的矩形网格拓扑。
以二重矩形网格拓扑为例，网格拓扑建立装置200可先在所有投影点构成的区域内建立一重网格拓扑，以通过一重网格控制整体曲面的形态和趋势，然后在各一重网格内建立二重网格，以将各一重网格划分为更小的区域，通过分区域构造曲面，使构造的曲面更精细和准确。也因此，N越大，网格越密，构造的曲面越精细和准确。
关于确定所有投影点构成的区域，作为示例，网格拓扑建立装置200可根据投影点的坐标来确定所有投影点在二维平面上构成的区域。例如，假设所有投影点的u坐标的最大值为umax，u坐标的最小值为umin，则沿u轴方向的边界为u＝umin-C，u＝umax+C，C为非负常数，同理，假设所有投影点的v坐标的最大值为vmax，v坐标的最小值为vmin，沿v轴方向的边界为v＝vmin-C、v＝vmax+C，通过上述边界即可确定所有投影点构成的区域。
网格点插值装置300用于根据各投影点的属性值对N重矩形网格拓扑的至少一个网格点进行插值计算，以获得至少一个网格点的属性值。至少一个网格点作为已知网格点，未获得属性值的网格点作为未知网格点。
作为示例，网格点插值装置300可按照从一重网格点至N重网格点的顺序，依次对N重矩形网格拓扑的每重网格点进行插值计算。当对i重网格点进行插值计算时，如果该i重网格点所在的i重网格中存在投影点，则根据存在的投影点的属性值对该i重网格点进行插值计算。i为大于0小于等于N的整数。
如果该i重网格点所在的i重网格中不存在投影点，则不对该i重网格点进行插值计算。
作为示例，网格点插值装置300可通过加权距离反比插值法或克里金插值法根据存在的投影点的属性值对该i重网格点进行插值计算。也可通过其它适合的方法根据存在的投影点的属性值对该i重网格点进行插值计算。
初始属性值计算装置400用于根据已知网格点的属性值计算未知网格点的初始属性值。
初始属性值计算装置400可使用各种适合的方法根据已知网格点的属性值计算未知网格点的初始属性值。
作为示例，初始属性值计算装置400可包括：插值单元(未示出)和计算单元(未示出)。
插值单元用于根据已知网格点的属性值对至少一个未知网格点进行插值计算，以获得所述至少一个未知网格点的初始属性值。
计算单元用于针对未通过插值单元进行插值计算的未知网格点，将距离该未知网格点最近的网格点的属性值或初始属性值的平均值作为该未知网格点的初始属性值。
迭代装置500用于根据已知网格点的属性值和未知网格点的初始属性值，通过有限差分求解格式对几何偏微分方程进行迭代求解，以获得未知网格点的属性值。
通过有限差分求解格式对几何偏微分方程进行迭代求解来获得未知网格点的属性值，避免了复杂的参数计算，大大降低了对操作者的数学专业知识的要求；可适用于含垂直断层、正断层和逆断层等复杂断层约束的曲面重构；引进的参数少，不需要大量的运算，能够实现含多种断层约束的地质曲面的快速重构。
几何偏微分方程可以是二阶几何偏微分方程、四阶几何偏微分方程等。选择不同的几何偏微分方程，可以实现不同特性的曲面的重建。例如，选择二阶几何偏微分方程，生成的曲面可在连接处(边界处)达到C0连续，即，相邻两个曲面具有交点；选择四阶几何偏微分方程，生成的曲面可在连接处达到C1连续，即，相邻两个曲面在交点处具有相同的一阶导数。
迭代装置500可将已知网格点的属性值和未知网格点的初始属性值作为初始值，通过式(2)～(5)对几何偏微分方程进行迭代求解，以获得未知网格点的属性值。
通过有限差分求解格式对四阶几何偏微分方程进行迭代求解来获得未知网格点的属性值，计算过程简单、计算时间短，且基于获得的未知网格点的属性值得到的曲面更为光滑。
地质曲面表示装置600获得由N重矩形网格拓扑的所有网格点及其属性值表示的地质曲面，即，重构的地质曲面。
图6示出根据本发明示例性实施例的插值单元的框图。
如图6所示，根据本发明示例性实施例的插值单元可包括：区域划分单元401和区域插值单元402。
具体说来，区域划分单元401用来根据断层多边形约束将所述由所有投影点构成的区域划分为多个子区域。
区域插值单元402用来对至少一个子区域内的未知网格点进行插值计算，以获得所述至少一个子区域内的未知网格点的初始属性值。
作为示例，区域插值单元402可分别针对上述子区域，对该子区域内具有3个以上已知网格点的行进行B样条曲线插值，以获得该行上的未知网格点的初始属性值。并对该子区域内具有3个以上已知网格点的列进行B样条曲线插值，以获得该列上的未知网格点的初始属性值。如果对该子区域内的未知网格点所在的行和列分别进行了B样条曲线插值，则将该未知网格点的初始属性值更新为分别获得的该未知网格点的初始属性值的平均值。
应该理解，区域插值单元402也可使用其它适合的方法对至少一个子区域内的未知网格点进行插值计算。
图7示出根据本发明示例性实施例的迭代装置的框图。
这里，N＝2，即，N重矩形网格拓扑为二重矩形网格拓扑，几何偏微分方程为四阶几何偏微分方程。
如图7所示，根据本发明示例性实施例的迭代装置500可包括：一重网格点迭代单元501、二重网格点迭代单元502和全局迭代单元503。
具体说来，一重网格点迭代单元501根据一重已知网格点的属性值和一重未知网格点的初始属性值，在第一预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得一重未知网格点的属性值。
一重已知网格点为一重网格点中的已知网格点，一重未知网格点为一重网格点中的未知网格点。
第一预设条件可为迭代次数达到第一预设迭代次数，或本次迭代求解得到的各一重未知网格点的属性值与上次迭代求解得到的各一重未知网格点的属性值的平均误差小于等于第一预设误差阈值。
具体说来，一重网格点迭代单元501将一重已知网格点的属性值和一重未知网格点的初始属性值作为迭代的初始值，通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，当满足第一预设条件时，认为达到稳态，迭代停止。
二重网格点迭代单元502分别针对各一重网格中的二重未知网格点，根据该一重网格中的二重已知网格点的属性值和二重未知网格点的初始属性值，在第二预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得所述一重网格中的二重未知网格点的属性值。
二重已知网格点为二重网格点中的已知网格点，二重未知网格点为二重网格点中的未知网格点。
第二预设条件可为迭代次数达到第二预设迭代次数，或本次迭代求解得到的各二重未知网格点的属性值与上次迭代求解得到的各二重未知网格点的属性值的平均误差小于等于第二预设误差阈值。
具体说来，二重网格点迭代单元502分别针对每个一重网格，将该一重网格中的二重已知网格点的属性值和二重未知网格点的初始属性值作为迭代的初始值，通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，当满足第二预设条件时，认为达到稳态，迭代停止。
应该理解，第二预设误差阈值可小于第一预设误差阈值。
全局迭代单元503根据所有已知网格点的属性值和通过一重网格点迭代单元501和二重网格点迭代单元502获得的未知网格点的属性值，在第三预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以更新所有未知网格点的属性值。
所有已知网格点包括一重已知网格点和二重已知网格点，所有未知网格点包括一重未知网格点和二重未知网格点。
第三预设条件可为迭代次数达到第三预设迭代次数，或本次迭代求解得到的各未知网格点的属性值与上次迭代求解得到的各未知网格点的属性值的平均误差小于等于第三预设误差阈值。
具体说来，全局迭代单元503将所有已知网格点的属性值和通过一重网格点迭代单元501和二重网格点迭代单元502获得的未知网格点的属性值作为迭代的初始值，通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，当满足第三预设条件时，认为达到稳态，迭代停止。
通过一重网格点迭代单元501和二重网格点迭代单元502分别迭代求解一重未知网格点和二重未知网格点的属性值，可以加快几何偏微分方程的收敛速度，实现曲面的快速重构。由于迭代求解各一重网格中的二重未知网格点时，一重网格作为边界不能跨越，可能使得一重网格两侧的二重网格点的属性值相差过大而导致最终生成的曲面不够平滑，因此通过全局迭代单元503进行全局迭代，可实现整个曲面在一重网格上实现光滑拼接，使重构的曲面达到期望的光滑度。
根据本发明示例性实施例的地质曲面重构方法及设备，可以实现含垂直断层、正断层和逆断层等复杂断层约束的曲面快速重构，且生成的曲面具有自然光滑度，无需再进行平滑处理。
虽然已表示和描述了本发明的一些示例性实施例，但本领域技术人员应该理解，在不脱离由权利要求及其等同物限定其范围的本发明的原理和精神的情况下，可以对这些实施例进行修改。
Claims (9)
1.一种地质曲面重构方法，包括：
a)将基于采集到的地质数据生成的三维的空间散点沿Z轴方向投影到二维平面，以获得各空间散点在二维平面上的投影点，其中，各空间散点的投影点的属性值为所述空间散点沿Z轴方向的坐标的值；
b)在由所有投影点构成的区域内建立N重矩形网格拓扑，其中，N为大于0的整数，且所述区域越大，N越大；
c)根据各投影点的属性值对N重矩形网格拓扑的至少一个网格点进行插值计算，以获得所述至少一个网格点的属性值，其中，所述至少一个网格点作为已知网格点，未获得属性值的网格点作为未知网格点；
d)根据已知网格点的属性值计算未知网格点的初始属性值；
e)根据已知网格点的属性值和未知网格点的初始属性值，通过有限差分求解格式对几何偏微分方程进行迭代求解，以获得未知网格点的属性值；
f)获得由N重矩形网格拓扑的所有网格点及其属性值表示的地质曲面。
2.如权利要求1所述的方法，其中，步骤c)包括：
按照从一重网格点至N重网格点的顺序，依次对N重矩形网格拓扑的每重网格点进行插值计算，其中，当对i重网格点进行插值计算时，如果所述i重网格点所在的i重网格中存在投影点，则根据存在的投影点的属性值对所述i重网格点进行插值计算，其中，i为大于0小于等于N的整数。
3.如权利要求2所述的方法，其中，在步骤c)中，通过加权距离反比插值法或克里金插值法根据存在的投影点的属性值对所述i重网格点进行插值计算。
4.如权利要求1所述的方法，其中，步骤d)包括：
D1)根据已知网格点的属性值对至少一个未知网格点进行插值计算，以获得所述至少一个未知网格点的初始属性值；
D2)针对未通过步骤D1)进行插值计算的未知网格点，将距离所述未知网格点最近的网格点的属性值或初始属性值的平均值作为所述未知网格点的初始属性值。
5.如权利要求4所述的方法，其中，步骤D1)包括：
D11)根据断层多边形约束将所述由所有投影点构成的区域划分为多个子区域；
D12)对至少一个子区域内的未知网格点进行插值计算，以获得所述至少一个子区域内的未知网格点的初始属性值。
6.如权利要求5所述的方法，其中，在步骤D12)中：
对所述子区域内具有3个以上已知网格点的行进行B样条曲线插值，以获得所述行上的未知网格点的初始属性值；
对所述子区域内具有3个以上已知网格点的列进行B样条曲线插值，以获得所述列上的未知网格点的初始属性值，
其中，如果对所述子区域内的未知网格点所在的行和列分别进行了B样条曲线插值，则将所述未知网格点的初始属性值更新为分别获得的所述未知网格点的初始属性值的平均值。
7.如权利要求1所述的方法，其中，N＝2，所述几何偏微分方程为四阶几何偏微分方程，步骤e)包括：
E1)根据一重已知网格点的属性值和一重未知网格点的初始属性值，在第一预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得一重未知网格点的属性值，其中，所述第一预设条件为迭代次数达到第一预设迭代次数，或本次迭代求解得到的各一重未知网格点的属性值与上次迭代求解得到的各一重未知网格点的属性值的平均误差小于等于第一预设误差阈值；
E2)分别针对各一重网格中的二重未知网格点，根据所述一重网格中的二重已知网格点的属性值和二重未知网格点的初始属性值，在第二预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得所述一重网格中的二重未知网格点的属性值，其中，所述第二预设条件为迭代次数达到第二预设迭代次数，或本次迭代求解得到的各二重未知网格点的属性值与上次迭代求解得到的各二重未知网格点的属性值的平均误差小于等于第二预设误差阈值；
E3)根据所有已知网格点的属性值和通过步骤E1)和E2)获得的未知网格点的属性值，在第三预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以更新所有未知网格点的属性值，其中，所述第三预设条件为迭代次数达到第三预设迭代次数，或本次迭代求解得到的各未知网格点的属性值与上次迭代求解得到的各未知网格点的属性值的平均误差小于等于第三预设误差阈值。
8.一种地质曲面重构设备，包括：
投影装置，将基于采集到的地质数据生成的三维的空间散点沿Z轴方向投影到二维平面，以获得各空间散点在二维平面上的投影点，其中，各空间散点的投影点的属性值为所述空间散点沿Z轴方向的坐标的值；
网格拓扑建立装置，在由所有投影点构成的区域内建立N重矩形网格拓扑，其中，N为大于0的整数，且所述区域越大，N越大；
网格点插值装置，根据各投影点的属性值对N重矩形网格拓扑的至少一个网格点进行插值计算，以获得所述至少一个网格点的属性值，其中，所述至少一个网格点作为已知网格点，未获得属性值的网格点作为未知网格点；
初始属性值计算装置，根据已知网格点的属性值计算未知网格点的初始属性值；
迭代装置，根据已知网格点的属性值和未知网格点的初始属性值，通过有限差分求解格式对几何偏微分方程进行迭代求解，以获得未知网格点的属性值；
地质曲面表示装置，获得由N重矩形网格拓扑的所有网格点及其属性值表示的地质曲面。
9.如权利要求8所述的设备，其中，N＝2，所述几何偏微分方程为四阶几何偏微分方程，迭代装置包括：
一重网格点迭代单元，根据一重已知网格点的属性值和一重未知网格点的初始属性值，在第一预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得一重未知网格点的属性值，其中，所述第一预设条件为迭代次数达到第一预设迭代次数，或本次迭代求解得到的各一重未知网格点的属性值与上次迭代求解得到的各一重未知网格点的属性值的平均误差小于等于第一预设误差阈值；
二重网格点迭代单元，分别针对各一重网格中的二重未知网格点，根据所述一重网格中的二重已知网格点的属性值和二重未知网格点的初始属性值，在第二预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以获得所述一重网格中的二重未知网格点的属性值，其中，所述第二预设条件为迭代次数达到第二预设迭代次数，或本次迭代求解得到的各二重未知网格点的属性值与上次迭代求解得到的各二重未知网格点的属性值的平均误差小于等于第二预设误差阈值；
全局迭代单元，根据所有已知网格点的属性值和通过一重网格点迭代单元和二重网格点迭代单元获得的未知网格点的属性值，在第三预设条件下通过有限差分求解格式对四阶几何偏微分方程进行迭代求解，以更新所有未知网格点的属性值，其中，所述第三预设条件为迭代次数达到第三预设迭代次数，或本次迭代求解得到的各未知网格点的属性值与上次迭代求解得到的各未知网格点的属性值的平均误差小于等于第三预设误差阈值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410461559.4A CN104240301B (zh) | 2014-09-11 | 2014-09-11 | 地质曲面重构方法及设备 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410461559.4A CN104240301B (zh) | 2014-09-11 | 2014-09-11 | 地质曲面重构方法及设备 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104240301A CN104240301A (zh) | 2014-12-24 |
| CN104240301B true CN104240301B (zh) | 2017-03-15 |
Family
ID=52228302
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410461559.4A Active CN104240301B (zh) | 2014-09-11 | 2014-09-11 | 地质曲面重构方法及设备 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104240301B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111721914A (zh) * | 2020-05-25 | 2020-09-29 | 成都理工大学 | 一种元素迁移能力的度量方法、系统、装置和存储介质 |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105118091B (zh) * | 2015-08-26 | 2019-01-18 | 中国电建集团北京勘测设计研究院有限公司 | 一种构建多精度非均匀地质网格曲面模型的方法和系统 |
| CN107590855B (zh) * | 2017-10-25 | 2020-12-08 | 中国石油集团东方地球物理勘探有限责任公司 | 曲面重构的多目标优化模型、插值方法及曲面重构方法 |
| CN111179401B (zh) * | 2019-12-31 | 2023-07-25 | 北京真景科技有限公司 | 基于3d曲面数据的拓扑分组方法 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7924278B2 (en) * | 2006-07-28 | 2011-04-12 | Microsoft Corporation | Real-time GPU rendering of piecewise algebraic surfaces |
| CN101582173B (zh) * | 2009-06-24 | 2012-07-11 | 中国石油天然气集团公司 | 复杂地质构造块状模型构建方法 |
| CN102074052A (zh) * | 2011-01-20 | 2011-05-25 | 山东理工大学 | 基于样点拓扑近邻的散乱点云曲面拓扑重建方法 |
| CN102222365B (zh) * | 2011-07-29 | 2012-08-29 | 电子科技大学 | 复杂空间的曲面重构方法 |
| CN102867330B (zh) * | 2012-08-29 | 2014-10-01 | 电子科技大学 | 基于区域划分的空间复杂层位重构方法 |
-
2014
- 2014-09-11 CN CN201410461559.4A patent/CN104240301B/zh active Active
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111721914A (zh) * | 2020-05-25 | 2020-09-29 | 成都理工大学 | 一种元素迁移能力的度量方法、系统、装置和存储介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104240301A (zh) | 2014-12-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104268934B (zh) | 一种由点云直接重建三维曲面的方法 | |
| EP3293552B1 (en) | System and method for editing geological models by switching between volume-based models and surface-based structural models augmented with stratigraphic fiber bundles | |
| CN104240301B (zh) | 地质曲面重构方法及设备 | |
| CN105354873B (zh) | 用于多孔介质三维重构的模式密度函数模拟方法 | |
| Yilmaz | The effect of interpolation methods in surface definition: an experimental study | |
| CN109584357A (zh) | 基于多轮廓线的三维建模方法、装置、系统及存储介质 | |
| CN107248142B (zh) | 一种文物碎片自动拼接方法 | |
| CN112116708B (zh) | 一种获取三维地质实体模型的方法及系统 | |
| US20220244424A1 (en) | Geological Grid Analysis | |
| CN107221028A (zh) | 一种基于地震解释数据的地质体闭合曲面三维重建方法 | |
| CN106875487B (zh) | 一种基于邻域作用力的地质六面体网格平滑方法 | |
| EP3828822A1 (en) | Civil engineering | |
| CN105844710B (zh) | 一种地质体网格化过程中的数据检测方法 | |
| CN106875484A (zh) | 一种基于三维地形的地质堆积体快速拟合建模方法 | |
| Kao et al. | Using particle swarm optimization to establish a local geometric geoid model | |
| Gosciewski | Reduction of deformations of the digital terrain model by merging interpolation algorithms | |
| CN105931297A (zh) | 三维地质表面模型中的数据处理方法 | |
| US11300706B2 (en) | Designing a geological simulation grid | |
| CN105069844B (zh) | 基于逻辑邻域的地质曲面拟合方法 | |
| CN112053438A (zh) | 基于水平集的成矿构造深部推断建模方法 | |
| US20210165935A1 (en) | Polyline contributor in civil engineering | |
| CN108876916A (zh) | 辫状河训练图像生成方法及设备 | |
| Vasta et al. | Modified radial basis functions approximation respecting data local features | |
| CN105869209A (zh) | 三维地质表面模型中的畸形三角形数据处理方法 | |
| CN105785444B (zh) | 速度场构建方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20180207 Address after: 072751 Zhuozhou, Baoding, Fan Yang Road West, No. 189 Patentee after: Dongfang Geophysical Exploration Co., Ltd., China Petrochemical Corp. Address before: Shuangliu County Huayang Huayang Road in Chengdu city of Sichuan Province in 610213 section of No. 216, Igawa geophysical exploration company of the Ministry of science and technology Patentee before: China National Petroleum Corporation Chuanqing Drilling Engineering Geophysical Exploration Company Ltd. |
|
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20200921 Address after: 100007 Beijing, Dongzhimen, North Street, No. 9, No. Patentee after: CHINA NATIONAL PETROLEUM Corp. Patentee after: BGP Inc., China National Petroleum Corp. Address before: 072751 Zhuozhou, Baoding, Fan Yang Road West, No. 189 Patentee before: BGP Inc., China National Petroleum Corp. |