一种基于大规模地震数据的三维地质曲面重建方法

一种基于大规模地震数据的三维地质曲面重建方法
【专利摘要】本发明提供了一种基于大规模地震数据的三维地质曲面重建方法。所述方法包括以下步骤：首先进行大规模地震数据的地形分块，具体包括将层位数据和地质约束多边形数据划分入不同地形块内并保存至外存中；然后以每个地形块为基本单位进行曲面重建，所述曲面重建步骤包括：从对应路径下的外存中读取地形块内部的层位数据和地质约束多边形数据，然后在地形块内构建初始网格，在每个初始网格内进行约束Delaunay三角剖分并进行三维空间三角网恢复，得到三维拓扑，最后对三维拓扑上的点进行曲面插值，得到三维地质曲面。本发明针对数据量很大的情况提出一种新的曲面重建方法，有效地利用有限的内存完成海量数据的处理。
【专利说明】
一种基于大规模地震数据的三维地质曲面重建方法
技术领域
[0001] 本发明涉及地质勘探技术领域，更具体地讲，涉及一种大规模工区产生的海量勘 探数据的曲面重建方法。
【背景技术】
[0002] 在三维地质建模的领域，对三维大规模海量地震数据进行曲面重建、三维可视化 乃至进行人机交互的操作，都已经成为常规而且必要的工业需求。然而，在大规模地震数据 进行曲面重建有一定的难度，原因之一在于原始数据量过大，数据范围过广。在现有的条件 下，一台普通PC机的硬件条件有限，尤其是内存方面，一般来说只有几GB到十几GB，大规模 地震勘探数据，原始数据一般有GB级别，在曲面重建过程中还会出现大量中间数据，因而难 以将所有大规模的地震勘探数据同时放入电脑内存中进行曲面重建。

【发明内容】

[0003] 针对现有技术中存在的不足，本发明的目的之一在于解决上述现有技术中存在的 一个或多个问题。例如，本发明的目的之一在于针对原始地震数据量很大的情况下提出一 种新的曲面重建方法，以有效地利用有限的内存完成海量数据的处理。
[0004] 为了实现上述目的，本发明提供了一种基于大规模地震数据的三维地质曲面重建 方法。所述方法包括以下步骤:首先进行大规模地震数据的地形分块，具体包括将层位数据 和地质约束多边形数据划分入不同地形块内并保存至外存中；然后以每个地形块为基本单 位进行曲面重建，所述曲面重建步骤包括:从对应路径下的外存中读取地形块内部的层位 数据和地质约束多边形数据，然后在地形块内构建初始网格，在每个初始网格内进行约束 Delaunay三角剖分并进行三维空间三角网恢复，得到三维拓扑，最后对三维拓扑上的点进 行曲面插值，得到三维地质曲面。
[0005] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，在对地 质约束多边形数据进行分块处理时，当一个地形块内只有一个地质约束多边形点时，求出 两个相邻地质约束多边形点连成的线段与地形块边界的交点，并加入对应的地形块内的地 质约束多边形数据中。
[0006] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，所述方 法还可以包括保留一份整体的封闭多边形数据，并且，在对每个地形块进行曲面重建操作 时，都要调用封闭多边形数据进行属性判断;具体地，定义封闭多边形数据集合P:
[0007] P:^^，···^}等式（1)
[0008] 其中
[0009] Pi= {upi，lpi，idi}等式（2)
[0010] 在等式⑴和等式(2)中，P表示封闭多边形的集合，Ρ,表示封闭多边形集合P中的 第i个封闭多边形，为封闭多边形的总个数，u Pi表示第i个封闭多边形的上盘多 边形点集合，Ip:表示第i个封闭多边形的下盘多边形点集合，icU表示第i个封闭多边形的根 多边形id。
[0011] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，所述方 法还可以包括引入地质约束多边形点的索引属性，以标识地质约束多边形点在整个地质约 束多边形中的位置，并用于判断地质约束多边形点之间是否相邻，只有索引相邻的两个地 质约束多边形点才能连接形成约束边，具体地，定义地质约束多边形点集合：
[0012] 等式(3)
[0013] 其中
[0014] Vj={Xj，Y」，IDj，Index」}等式（4)
[0015] 在等式(3)和等式(4)中，VP表示一个地质约束多边形上所有地质约束多边形点的 集合，Vj表示地质约束多边形上第j个地质约束多边形点，l<j<n，n为一个地质约束多边 形上地质约束多边形点的总个数，Xj和Yj表示地质约束多边形点的位置，IDj表示第j个地质 约束多边形点所在的地质约束多边形在地质约束多边形集合中的ID，Indexj表示第j个地 质约束多边形点在地质约束多边形中的位置索引。
[0016] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，所述方 法还可以包括在内存开辟缓存区，在曲面重建操作中，执行外存与缓存区间的交互以及内 存与缓存区间的交互，具体地，在进行缓存区地形块调度时，首先判断当前处理的地形块及 其邻近地形块（即与当前处理的地形块邻近的地形块)是否读入缓存区，如果在缓存区中找 到当前处理地形块及其临近地形块，则直接从缓存区中选取对应地形块;如果在缓存区中 没有找到当前处理地形块及其临近地形块，则判断缓存区是否已满，若缓存区未满且缓存 区剩余空间足够，从外存中将需要的数据块添加到缓存区中；若缓存区已满或者缓存区剩 余空间小于处理当前块所需要的空间，则采用地形块置换策略，将距离当前处理的地形块 最远的地形块移出缓存区，将当前处理的地形块的临近块移入缓存区。
[0017] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，所述插 值步骤包括对地形块的非边界点和边界点进行深度插值，其中，对非边界点进行深度插值 的步骤为对网格点插值，再对断层多边形点进行插值，最后进行补偿插值;对地形块的边界 点进行深度插值的过程中，所有与当前边界点拓扑连接的点都要作为种子点参与当前边界 点的深度插值计算;并且，对边界点添加插值与否的属性，将插值后的边界点的属性设置为 已插值，对已插值的点不进行二次插值。
[0018] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，所述方 法读取大规模地震数据中的预定行，并对所述预定行的地震数据进行地形分块处理，然后 重复进行所述读取和地形分块操作，直至大规模地震数据被读取和处理完。
[0019] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，所述对 初始给定的地形块进行扩充时保证整个区域呈矩形。
[0020] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，在所述 保存步骤中，为每个地形块设置一个id，并采用二维数组的方式进行地形块的id索引，确定 一个二维数组a[M][N]，其中所述二维数组的第一维表示地形块的X方向，二维数组的第二 维表示地形块的Y方向，二维数组中的值为地形块的id。
[0021] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，所述方 法还可以包括在所述地形块内构建初始网格步骤之前，对地形块中的层位数据和地质约束 多边形数据进行预处理。
[0022] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，所述方 法对不同的地形块，根据地形的分布确定不同的分辨率，具体地，采用能量函数法计算地形 块初始的分辨率，然后对计算出来的分辨率值结合地形块内的断层分布状况进行微调，如 果一个地形块内存在的断层约束多边形越多，分辨率越高。
[0023] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，所述方 法采用离散L0D模型进行曲面重建。
[0024] 根据本发明基于大规模地震数据的三维地质曲面重建方法的一个实施例，所述方 法采用Constrained Delaunay约束三角剖分，例如，循环换边法。
[0025] 与现有技术相比，本发明的有益效果包括:有效地解决了单台PC机内存难以实现 对整个大规模地质工区的勘探数据进行整体曲面重建的问题。
【附图说明】
[0026] 通过下面结合附图进行的描述，本发明的上述和其他目的和特点将会变得更加清 楚，其中：
[0027] 图1是对原始点集进行分割归并法进行三角剖分示意图。
[0028] 图2是Constrained Delaunay三角剖分流程示意图。
[0029] 图3基于地形分块的大规模地震数据曲面重建整体流程。
[0030] 图4是根据本发明示例性实施例的基于大规模地震数据的三维地质曲面重建方法 对大规模地震原始数据进行地形分块示意图。
[0031] 图5是根据本发明示例性实施例的基于大规模地震数据的三维地质曲面重建方法 动态地形分块的示意图。
[0032] 图6是地形块数据在硬盘中分布示意图。
[0033] 图7是根据本发明示例性实施例的基于大规模地震数据的三维地质曲面重建方法 地质多边形数据划分至地形块示意图。
[0034] 图8是地形块中封闭多边形示意图。
[0035] 图9是地质多边形重复跨越同一地形块示意图。
[0036] 图10是地根据本发明示例性实施例的基于大规模地震数据的三维地质曲面重建 方法形块与缓存区域示意图。
[0037] 图11是根据本发明示例性实施例的基于大规模地震数据的三维地质曲面重建方 法缓存区置换地形块示意图。
[0038] 图12是根据本发明示例性实施例的基于大规模地震数据的三维地质曲面重建方 法采用缓存区对地形块进行数据调度的整体流程示意图。
[0039] 图13是根据本发明示例性实施例的基于大规模地震数据的三维地质曲面重建方 法地形块边界点插值示意图。
【具体实施方式】
[0040] 在下文中，将结合附图和示例性实施例详细地描述根据本发明的基于大规模地震 数据的三维地质曲面重建方法。
[0041] 根据本发明示例性实施例的基于大规模地震数据的三维地质曲面重建方法首先 进行大规模地震数据的地形分块，然后以每个地形块(简称为块)为基本单位进行曲面重 建。主要创新点在于采用三维空间三角网恢复的方法解决大规模地震数据的曲面重建难 题。
[0042] 本发明主要解决的问题是：首先要完成大规模地震数据的地形分块，其中涉及到 层位数据的分块以及约束多边形数据的分块。其次，为了进一步简化地形数据，引入多分辨 的概念，不同地形块的数据分辨率与具体块内的地形相关，相邻的地形块可能由于分辨率 的不同，而产生裂缝现象，应该采取一定的措施加以避免。再次，将整个大规模地质工区划 分为地形块，虽然解决了大规模地质数据下曲面重建内存不足的问题，但是地形块数量假 如过多，势必会导致内存与硬盘间的频繁交互，因此，确定一个合适的内存调度算法，以实 现相对高效的曲面重建过程也是很必要的。最后，要解决多个地形块之间边界的插值问题。 [0043]为便于理解，
【申请人】就本申请所涉及的理论进行简单介绍。
[0044] 1、Delaunay三角剖分技术
[0045] Delaunay三角剖分的定义如下：假设平面上点集合为V，边集合为E，对于E中的任 意一条边e，如果存在圆C，只包含e的两个端点，而不包含其他任意点，则边e为一条 Delaunay边。假设T为平面散点集合的一个三角剖分，如果T的所有边均为Delaunay边，则T 为Delaunay三角剖分。
[0046] Delaunay三角剖分是一种优化的三角剖分方法，产生的三角网有诸多良好的特 性，比如外接圆的空圆性、最小内角最大性、唯一性、局部性、凸包性等。
[0047]下面简要介绍一下Delaunay三角剖分的几个常用方法：
[0048] (1)逐点插入法
[0049]逐点插入法的流程是:首先，确定能够涵盖所有离散点数据的一个凸包作为边界 多边形。其次，依次将原始点数据插入边界多边形中，每次插入点时，对新插入点的位置加 以判断，连接新插入点与其周围的三角形顶点。再次，对连接好的新边，依照外接圆准则拉 进行L0P优化。最终，所有点都被插入，完成三角网的构建。现在常用的一些De launay三角剖 分，生成三角网格的方法都可归类为逐点生成法，只是在一些数据结构还有初始的边界多 边形构建方面有所区分。
[0050] (2)分割归并算法(分治法）
[0051] 将复杂的问题先划分为比较简单的子问题进行求解，再将所有子问题以一定的规 则合并，解决原问题的思想称为分治思想。采用分割归并法进行Delaunay三角剖分主要有 以下几个步骤:首先，对离散的点集合数据进行排序，一般的顺序是X方向由左到右，Y方向 由上到下，先对X方向排完序以后再对Y方向排序。其次，对散点集合数据进行递归划分，直 到所有划分出来的子集合内的点数不大于3。再次，将所有子集合中的点连接为三角形。最 终，合并子集三角网，每次合并时依照外接圆准则进行局部优化，使合并之后仍为Delaunay 三角网，最终得到整体的Delaunay三角网。如图1所示为对6个原始散点采用分割归并算法 进行三角剖分的过程示意图，其中（a)显示了经过排序后的原始数据，（b)表示分别对子集 进行划分，（c)表示对子集进行合并。
[0052] (3)三角网增长法
[0053]三角网增长法的主要流程如下：首先，在所有离散点集合V中选取一点A，并找到欧 氏距离上距离A最近的一点B，得到边AB。其次，以AB边为基础，寻找符合Delaunay三角剖分 外接圆判定准则的第三点C，得到AABC。再次，以AABC的三条边为基础，依次寻找符合条件 的三角形，最终得到Delaunay三角剖分产生的三角网。
[0054] 2 ·约束Delaunay三角网技术
[0055] 在实际的地质建模过程中，地质勘探数据往往不是独立的，而是彼此之间存在一 定的约束关系，比较常见的约束关系有地层的断裂线、逆断层的断层上下盘走向线、地质勘 探人员人为绘制的地势走向约束线等等。由于这些约束的存在，采用普通的Delaunay剖分 方法产生的三角网虽然能够得到一定的数字地面模型，但是无法呈现出地表实际的复杂约 束关系，也就对实际地质建模的需求上贡献较小。于是，在Delaunay三角剖分的过程中加入 一定的约束关系，成为非常关键的技术及研究方向。实现中常用的是Constrained Delaunay约束三角剖分。
[0056] Constrained Delaunay约束三角剖分是将约束边进行整条嵌入的，并不进行加点 操作，也不进行分割，所以对原始拓扑的影响小于Conforming Delaunay方法。也是本发明 在复杂约束下进行Delaunay三角剖分使用的方法。循环换边法作为一种比较有代表性的 Constrained Delaunay约束三角剖分方法，其步骤如下：
[0057] (1)将约束边两个端点插入原始的非约束三角网。
[0058] (2)计算三角网中与约束边相交的非端点相交边，得到相交边集合E={eo，ei, ...， θη} 〇
[0059] (3)遍历找到的相交边集合Ε，对每条边&，找到以&为对角线的两个三角形，判断 两个三角形所组成的四边形的凹凸性。如果四边形为严格凸多边形，则交换四边形对角线， 否则跳过中的其他边进行判断。如确定要交换四边形对角线，则将边 eiWE中去除，判 断交换对角线后产生的新边是否与约束边非端点相交。对相交边集合E进行更新。
[0060] (4) 一轮遍历结束后，如果不存在非端点相交的边，说明约束边已经被完整嵌入。 如果仍存在，则重复步骤(3)，直至约束边与每条三角网格中的边均无非端点相交。
[0061 ] Constrained Delaunay约束三角剖分的流程图如图2所示，流程顺序按照从(a)到 ⑴。
[0062] 3、大规模地震数据下曲面重建的理论基础
[0063]大规模的工区地质勘探数据，特点为数据区域范围大，数据点多。在某些超大的工 区下，地质勘探的原始数据甚至达到了G或者T级别，而现有的普通PC机的内存范围却只有 4G-32G不等。由于受到客观硬件条件的制约，为了解决大规模地震数据下进行三维地质建 模的问题，很多专家学者进行了研究和探索，主要方法是对原始地震数据进行处理，减少需 要处理的总数据量和一次读入内存的数据量。地形分块技术和层次细节技术也是针对于大 规模数据进行地质建模的过程中数据组织技术，接下来对其加以介绍。
[0064] 3.1地形分块技术
[0065]地形分块技术，顾名思义，就是将整个地质工区按照所在的地理坐标，划分为面积 一样的若干块。在单机模式下，通过合理的分块方案，可以解决任意大小的大规模地震数据 的曲面重建与绘制问题。对存放在外存的大规模地震数据进行基于空间分布的分块操作， 使每个产生的地形块数据均可以一次性装入内存。地形分块步骤有以下几步：
[0066] (1)确定每一块的大小，常用的为大小相同的正方形块。
[0067] (2)将地震数据放入到块中，若某一块不存在数据，则为了整体的完整性添加空 白。
[0068] (3)如果块内存在网格，则尽量保证每块的网格点数为(2k+l)X(2 k+l)。
[0069] (4)将生成好的地形块数据文件保存在硬盘中，方便实时调用。
[0070] 在实际的地形分块过程中，有几个问题是需要考虑的：首先就是分块的大小，若地 形分块过大，可能内存无法存放单个块文件，影响内存效率。如果分块过小，虽然内存方面 不存在问题，但是可能会导致过多的内外存I/O交互，影响整体效率，同时，块文件过多，也 会增加额外的外存存储管理开销。其次，每个块边界上的数据点，不止受当前块的影响，周 围块也会对其有影响，所以进入内存处理时，往往是同时使几个相邻块进入内存，这时就涉 及到缓存设计的问题，即如何建立缓存机制来处理内存与外存间的频繁交互问题。
[0071] 3.2层次细节技术
[0072] 层次细节L0D(Level of Detail)技术常用于实时显不，在实际曲面重建以及显不 的过程中，有的区域地质勘探点比较密集，且存在断层或者褶皱等地质规律约束，在这些区 域往往要应用较多的原始点，构建比较密集的网格数据。而在另一些区域，地形比较平坦， 主要以层位为主，不存在约束，用少部分原始点构建简单网格即可完成区域重建。这样，在 不同层次和不同范围，采用不同分辨率的网格对原始数据进行重建，既保证了运行的效率， 又保证了重建的曲面效果。
[0073]目前常用的L0D模型有以下两种：
[0074] (1)离散1/?模型
[0075]离散L0D模型指的是事先通过计算，得到整个工区的多个副本，每个副本对应着一 重分辨率，这样得到了一个基于不同分辨率的金字塔状数据模型。在曲面重建与显示时，根 据实际要求选取对应分辨率。该方法实现较简单，但是在不同分辨率的区域切换时，会产生 比较明显的跳变感。而且会存储冗余的分辨率网格点，造成存储的浪费。
[0076] (2)连续的 LOD(CLOD)模型
[0077] 因为离散的L0D模型在很多地方有所不足，连续L0D模型的概念被提出。连续L0D模 型的概念和离散L0D类似，区别在于连续L0D模型并不事先存储好不同分辨率的数据，而是 实时计算具体的分辨率。
[0078] 离散L0D模型实现较简单，计算量不大，一次生成多次使用，对CPU的负担也较小。 缺点在于要存储多个分辨率下的网格数据，还有跳变现象。连续L0D模型在显示上效果要优 于离散L0D模型，但计算较多，复杂且耗时。由于本发明并不涉及分辨率的动态变化，在曲面 重建的过程中，数据只需保存一次，所以综合考虑采用离散L0D模型的思想。
[0079] 在一个示例性实施例中，采用基于地形分块的方法进行大规模地震数据下的曲面 重建，具体包括以下操作：
[0080] A、大规模地震数据的地形分块整体流程
[0081] 首先，将原始的大规模地震数据和约束多边形数据按照地形划分入不同地形块 内，并保存至外存(例如，硬盘)中。其次，以每个地形块为基本单位进行曲面重建，先对地形 块中的层位数据和断层数据进行预处理，然后对处理后的数据进行二维网格构建，根据网 格数据构建二维Delaunay三角剖分并进行三维空间三角网恢复，最后进行曲面插值并保存 重建后的曲面拓扑数据，在显示时将所有地形块内的曲面拓扑一起显示，从而实现整个大 规模工区的曲面重建和显示。大规模地质工区地震数据曲面重建方法的整体流程如图3所 示。这里，本申请中出现的约束多边形和地质多边形是地质约束多边形的简称，本发明说的 地质约束多边形特指断层多边形。
[0082]由于单台PC机内存难以实现对整个大规模地质工区的勘探数据进行整体曲面重 建，所以使用地形分块技术，将整体区域划分为若干地形块，如图4所示。由于在整体数据全 部读入完之前，无法确定整体的最大值最小值，也无法确定地形块的个数，所以要使用地形 块动态添加的方法，动态对地形块的地形和属性进行判定。
[0083] a、大规模地震数据下地形分炔基本流程
[0084]大规模地震数据的地形分块的步骤主要有以下几步：
[0085] (1)给定初始地形块属性
[0086]主要有以下几项:地形块的长度、X方向地形块个数、Y方向地形块个数。这几项属 性由工区的实际大小决定，可以由外部进行指定，X方向与Y方向的地形块的最小值为1，即 最起码对于大规模地震工区数据存在一个地形块，此时该块的范围为整个地质工区。
[0087] (2)读取原始数据
[0088] 读取地质勘探原始数据，包括层位数据和断层多边形数据。在大规模工区数据下 原始地质勘探数据的数据量较大，如果在所有数据均从文件中读取之后再进行地形分块操 作，数据操作量较大，花费时间较长。而如果每读入一个数据就进行分块操作，操作会相对 频繁，影响效率。因此使用类似分块的思想进行数据读入，即每次读取地震数据文件中的若 干行，对这若干行数据进行地形分块操作，操作完以后再继续进行读入操作，直至地震数据 文件被读取和处理完。
[0089] 在一个实施例中，地形分块操作具体为:初始化对全部原始地震数据进行读入，获 取地震数据的最大最小值，根据地震数据的具体范围以及制定的单个块长度确定地形块在 X方向和Y方向的个数。但是，在大规模工区地震数据的曲面重建过程中，数据量较大，如果 按照上述操作，预先对原始数据进行全部读入，会消耗一定的时间，为了提高处理效率，本 发明优选使用一种动态添加地形块的方法，即初始给定一定范围的地形块，如果读入的地 震数据在初始给定的地形块内，则将地震数据分配到地形块中。如果读入的地震数据在初 始给定的地形块之外，则对给定的地形块进行扩充，使之包含读入地震数据。进一步优选 地，为了方便曲面重建操作和曲面显示，对地形块扩充时要保证整个区域呈矩形状。
[0090] 如图5所示，初始设置默认X方向和Y方向各有两个地形块，对数据进行读入时，h、 V2、V3、V 4等点经判断在已经存在的地形块内，就直接添加进相应的地形块内部。点V5不在已 经存在的地形块内部，就对现有地形块进行动态添加，添加新的地形块以包含点V 5，同时为 了保证整体区域的完整性，添加了空白的地形块。
[0091] (4)将数据分布到地形块中
[0092] 原始地质勘探数据主要是层位数据和断层数据，按照地震数据的二维坐标，寻找 对应的地形块，将地震数据放入对应的地形块即可。
[0093] (5)保存地形块内数据
[0094]将地形块内的地震数据保存至硬盘中，为了方便确定保存位置，对每个地形块给 定一个块id，本发明使用二维数组的方式进行块的id索引，确定一个二维数组a[M][N]，其 中数组的第一维表示地形块的X方向，数组的第二维表示地形块的Y方向，数组中的值为地 形块的id。例如，a [ Ο ] [ Ο ]的值为Ο，说明X方向的第0个、Y方向的第0个地形块中块id为0。随 着地形块的动态添加，二维数组a[M][N]也会动态变化，时刻改变对应位置的地形块id。同 一地形块内的文件最终保存在同一个地形块id对应的路径下。
[0095] 经过以上的操作后，完成对所有原始地震数据的读入，原始的大规模地震数据转 化为地形块数据并保存在硬盘中，为下一步操作做准备，地形块数据在硬盘中的分布如图6 所示，后缀名为block的文件为地形块文件，地形块文件以"地形块idj也形块分辨率"的规 则命名。
[0096] b、大规模地震数据下复杂约束分块处理。
[0097] 在大规模地震数据的处理中，除原始的层位数据之外，还存在多种复杂约束信息， 以地质约束多边形的形式存在。其中，断层多边形数据为地质勘探人员人工绘制的多边形， 在本发明中，地质约束多边形指断层多边形，另外，封闭多边形也特指断层多边形，也就是 说本发明所有的多边形都是指断层多边形。在大规模地震数据下进行地形分块处理时，地 质约束多边形数据的处理有以下几个难点：
[0098]首先要考虑的是地质约束多边形点的分布问题。地质多边形在二维平面上均显示 为折线段的形式，地质约束多边形数据相邻的两个点往往间距较远，在实际地质数据分布 中，很有可能出现两个相邻的多边形点跨越地形块的情况，此时可能在一个地形块内只有 一个地质多边形点。由于本发明使用地形块作为曲面重建的基本单位，在每个地形块内只 有一个地质多边形点的情况下，无法构建约束边实现后续的Delaunay约束三角剖分运算， 所以要做出相应处理。在一个示例性实施例中，采取的处理方法是求出这两个相邻地质多 边形点连成的线段与地形块边界的交点，并加入对应的地形块内的地质多边形数据中。如 图7所示，为地质多边形上的三个连接点，首先找到地质约束多边形与地形块的交 点V 4、V5，将V4、V5点分别放入对应的地形块，此时要注意的是由于是边界点，所以在相邻的 地形块中都要放入。进行操作以后，地形块1内的地质多边形点为V 3、V4,块2内的多边形点为 ，块3内的地质多边形点为V5、V 2。
[0099]大规模地震数据下的复杂约束信息划分至地形块内的另一个复杂之处在于，地质 约束多边形中存在不封闭的垂直断层多边形，也存在逆断层多边形、正断层多边形、倒转多 边形等封闭的多边形，如图8所示。进行地质曲面重建的过程中会应用多边形的边界数据对 点的重数和属性进行判定。而封闭多边形的分块操作可能会出现与实际不符的情况，如图8 中的块1内部，点A为地质约束多边形点，点B、C为地质约束多边形与地形块的交点，块1中的 约束多边形被判断为点A、B、C连接而成的封闭多边形，明显与实际不符。所以一方面要将地 质约束多边形的点数据划分到各个地形块内，另一方面要保留一份整体的封闭地质多边形 数据(简称封闭多边形）。
[0100]定义封闭多边形集合P:
[0101] 卩二沾七，…'}等式(1)
[0102] 多边形数据如下，包含上盘多边形点集合，下盘多边形点集合与根多边形id
[0103] Pi= {upi，lpi，idi}等式（2)
[0104] 也即，在等式(2)中，uPi表示上盘多边形点集合，lPi表示下盘多边形点集合，idi表 示根多边形id，在二维平面上，对所有的点都给了一个id，这里根多边形icU指的是当前第i 个多边形中该点在二维点集中的i d。
[0105] 也就是说，在等式（1)和等式(2)中，P表示所有封闭多边形数据的集合，每一fPi 都是一个封闭多边形，集Ρ:表示封闭多边形集合P中的第i个封闭多边形，为封闭 多边形的总个数，等式(2)是针对等式（1)某一个多边形Pi的数据结构，每一个多边形结构 都是一样，u Pl表示第i个封闭多边形的上盘多边形点集合，lPl表示第i个封闭多边形的下盘 多边形点集合，icU表示第i个封闭多边形的根多边形id。
[0106] 封闭多边形数据是作为一个整体存在的，对每个地形块进行曲面重建操作时，都 要调用封闭多边形数据进行属性判断。
[0107] 还有一点需要考虑的，就是地质多边形线段重复经过同一个地形块的情况。如图9 所示，地质多边形重复进入地形块1和地形块2，此时，如果按照上文提到的划分方法，则块1 内部的多边形点依次为:乂:^^^^^由于地形块是进行曲面重建的基本单位名块 内进行约束Delaunay三角剖分时是以多边形的线段组合作为约束的，即V 2、V5组成的线段也 会成为约束，而实际上这两个点并不相连，也不存在约束。只有索引相邻的两点才能连接形 成约束边。为了解决这个问题，引入一个属性:地质多边形点的索引。即当前的一个地质多 边形点在整个多边形中处于怎样的一个位置。在此定义多边形点集合：
[0108] VelVhVf.Jn}等式（3)
[0109] 其中Vi为第i个点，Vi中的数据包括多边形点的位置、点对应的多边形i d、点在多边 形中的位置索引。
[0110] Vi={Xi,Yi JDiJndexi}等式（4)
[0111] 也即，在等式(4)中，Xi和Yi表示地质约束多边形点的位置，IDi表示第i个点对应的 地质约束多边形ID，Indexi表示第i个点在地质约束多边形中的位置索引；l$i彡n，n为地 质约束多边形点的总个数。
[0112]需要补充说明的是，等式(3)中的VP指的是某一个地质约束多边形上所有点（即地 质约束多边形点)集合，等式(4)中Vi为该地质约束多边形上第i个点，Xi和Yi表示地质约束 多边形点的位置。由于每个多边形都有一个ID，所以IDi表示的是该点所在的地质多边形它 在地质多边形集合的多边形ID，Inde Xl才表示该点在该地质约束多边形上的位置索引，是 该点在该地质约束多边形上的唯一标识。
[0113]以图9为例，引入多边形点的索引属性后，经判断可知…、…在原始地质多边形上的 位置值并不相邻，不会连接成为约束边，保证了结果的正确性。
[0114] B、大规模地震数据的多分辨率设置
[0115] 为了进一步对大规模的地震数据进行简化，提高曲面重建的效率，本发明示例性 实施例使用多分辨的层次细节方法，对不同的地形块，根据地形的分布确定不同的分辨率。
[0116] 在一个示例性实施例中，采用的是能量函数法，设z = H(x，y)为点（x，y)处的深度 值，对于一块η个数据点的地形区域，定义f。为该区域所有点的能量函数。
[0117] 1 等式(5)
[0118] 其中Pi为该区域的相对深度分布：
等式(6)
[0120] 在等式(6)中，Zo为该区域中最小深度值，ε是为了防止分母为零添加的很小的安 全的数。
[0121] 采用能量函数的方法，分辨率能够很好的体现出地形的起伏状况，由于所有点都 参与了运算，防止了出现局部点过密的情况导致分辨率偏差。本发明采用能量函数法计算 初始的分辨率，然后对计算出来的分辨率值结合块内的断层分布状况进行微调，如果一个 地形块内存在的断层或者倒转约束较多，分辨率就会相对较高。
[0122] 在本发明示例性实施例中，实际曲面拓扑重建的方法是在地形块内构建初始四边 形网格，在四边形网格内进行约束De 1 aunay三角剖分得到二维拓扑，最后将二维拓扑恢复 为三维拓扑。每个地形块内需要构建一定过个数的四边形网格，地形块内的分辨率与网格 密度直接相关，假设在一个地形块内，分辨率为P，则该地形块内构建的初始四边形网格的 个数为2 PX2P。
[0123] C、大规模地震数据下的地形块调度
[0124] 由于数据规模和指定块长度的问题，有可能进行地形分块操作之后产生多个地形 块，此时如果逐一对地形块进行曲面重建操作，可能会涉及到多次内存与硬盘的交互，影响 曲面重建效率。考虑到以上的问题，本发明示例性实施例进一步开辟一部分内存空间作为 缓存区，在曲面重建操作中，执行硬盘-缓存区间的交互以及内存-缓存区间的交互，这样一 方面避免了地形块个数过多时产生的频繁交互，另一方面可以保证地形块内拓扑点的深度 值的正确性。
[0125] 对于每个地形块而言，在曲面三维拓扑重建过程中，地形块边界上的网格点的深 度值不能仅由当前地形块内部的点进行插值操作获取，还要将临近块内的点的深度值考虑 进来，涉及到了地形块的边界缝合问题，这一点就对缓存区大小的设置以及地形块在缓存 区内的调度提出了要求，大规模地震数据下的地形块调度主要有以下几个方面：
[0126] (1)缓存区设置
[0127] 缓存区的大小与地形块的个数和每个地形块的大小相关，如果单个地形块数据量 较大，整个内存都无法放入多个地形块，自然会限制缓存区内的地形块个数。如果地形块个 数比较少且数据量不大，完全可以将所有地形块放入缓存区。对于一般的大规模地震数据， 地形块较多，整体数据较大，对每个地形块操作时不仅要用到当前地形块内的数据，还要用 到其周围8个地形块的数据(边界区域地形块除外），因此缓存区设置的最小值至少为9个地 形块。
[0128] (2)缓存区内添加地形块
[0129] 进行缓存区地形块调度时，首先在缓存区中寻找当前处理的地形块和其临近地形 块。在上文地形分块的过程中，已经建立了表示块索引关系的二维数组a[M][N]，通过该二 维数组，对某一个地形块进行处理时，能够快速的得到其相邻地形块的id。即二维数组中的 相邻关系就对应了地形块间实际的相邻关系，如当前地形块的位置对应到二维数组中为a [m] [η]，则其相邻的地形块为a[m-l][n-l]、a[m-l] [n]、a[m-l][n+l]、a[m][n+l]、a[m+l][n + 1]、a[m+l] [n]、a[m+l] [n-1]、a[m] [n-1]所对应的地形块，进而得到地形块id以及对应的 块数据文件。如果在缓存区中找到当前处理地形块及其临近地形块，直接从缓存区中选取 对应地形块进行操作。否则判断缓存区是否已满，如果缓存区未满且缓存区剩余空间足够， 从硬盘文件中将需要的数据块文件添加到缓存区中。如图10所示，红色区域为当前处理的 地形块1，块2、3、4为其临近地形块，在对块1进行处理时，在缓存区内寻找块1、2、3、4，未找 到对应地形块且缓存区未满，则从硬盘中添加块1、2、3、4的数据到缓存区域内。
[0130] (3)缓存区内进行地形块置换
[0131]如果缓存区已满或者缓存区剩余空间小于处理当前块所需要的空间，要采用地形 块置换策略，用新的地形块数据替换将缓存区内的一部分地形块。常用的缓存区调度方法 主要根据地形块进入缓存区的时间以及使用次数决定。对于本发明而言，要考虑地形块之 间的临近关系，保证对某个地形块进行处理时，该地形块的临近地形块均在缓存区内。在此 基础上本发明使用一种基于距离的地形块置换方法，对当前地形块进行处理时，如果缓存 区空间不足，将距离当前地形块最远的地形块移出缓存区，将当前地形块的临近块移入缓 存区。如图11所示，初始状态处理的地形块为块12,块12及其临近块都在缓存区内。处理完 块12后对块13进行处理，需要将块13的临近块移入缓存区，此时缓存区已满，经过位置判 断，块13的临近块中块7、块14、块17不在缓存区内，将缓存区中距离块13最远的三块:块10、 块11、块20移出缓存区，将块7、块14、块17移入缓存区。
[0132] 采用缓存区对地形块进行数据调度的整体流程示意图如图12所示。
[0133] D、大规模地震数据下的地形块插值
[0134] 地形块内执行的操作，是在一个地形块内，先从对应路径下的硬盘文件中读取该 地形块内部的层位数据文件和约束多边形数据文件，再根据本地形块的分辨率建立一定步 长的初始网格，在每个网格内部进行约束Delaunay三角剖分构建，得到三角形网格拓扑，再 对构建好的二维拓扑进行拼接操作，得到三维拓扑，最终对拓扑上的点进行深度插值，得到 真实的三维地质曲面。地形块内点深度插值主要有以下两步：
[0135] (1)非边界点插值
[0136] 对于非地形块边界点，深度插值步骤是先对网格点插值，再对断层多边形点进行 插值，最后进行补偿插值。
[0137] 首先对网格点进行插值，连接好的曲面拓扑，网格点占据点的绝大多数，所以先对 网格点进行插值，在曲面拓扑的建立过程中，已经把原始的层位点作为种子点分配到对应 的网格中，网格点的值由对应网格的种子点进行插值而得。其次对多边形上的点进行插值。 多边形上的点相对于网格点来说数量较小，而且一般均与网格点连成拓扑边。完成网格点 的插值过程后，得到了网格点的实际深度值，可以使用网格点数据对与网格点相连接的多 边形点进行插值。
[0138] 最后进行补偿插值，经过网格插值和多边形点插值两步，只有少部分点可能依然 没有实际深度值，此时要进行补偿插值操作，使用经过上面两步所得到的拓扑中已经有实 际深度值的点，对仍然没有实际深度值的点进行插值。补偿插值不规定次数，直至所有拓扑 连接上的点均具有实际深度值。
[0139] (2)边界点插值
[0140] 对地形块的边界点进行深度插值的过程中，不仅要考虑当前地形块内点的影响， 还要考虑相邻块内部的点对其进行的影响。如图13所示，在左上角地形块的边界点A，对点A 进行深度值插值时，利用的是与点A连接拓扑的其他点，在本地形块内部的拓扑中点A和点 B、C、F、G相连接，而在其临近地形块内点A和点C、D、E、F相连接，以上所有与点A拓扑连接的 点都要作为种子点参与点A的深度插值计算。
[0141] 另一个需要注意的地形块边界点插值问题的是不同地形块之间的重复插值问题， 即多个地形块均对同一个点进行插值的情况。以图13中的点A为例，在对左下地形块进行点 插值时，已经寻找了各个地形块中与点A相连接的点对点A进行插值，在对左上地形块中进 行点插值操作时还会对点A进行一次插值，对相同点的重复插值操作势必会对整体算法效 率产生影响。为了避免重复插值，对边界点添加插值与否的属性，每次对边界点插值后该点 的属性设置为已插值，对已插值的点不进行二次插值。
[0142] 综上所述，本发明结合大规模地质工区地震数据的实际情况，提出了大规模地震 数据三维地质曲面重建的方法，将整体区域划分为若干地形块，有效地解决了单台PC机内 存难以实现对整个大规模地质工区的勘探数据进行整体曲面重建的问题。本发明对海量数 据的处理有着非常重要的意义，在物探领域应用前景广阔。
[0143] 尽管上面已经通过结合示例性实施例描述了本发明，但是本领域技术人员应该清 楚，在不脱离权利要求所限定的精神和范围的情况下，可对本发明的示例性实施例进行各 种修改和改变。
【主权项】
1. 一种基于大规模地震数据的三维地质曲面重建方法，其特征在于，所述方法包括以 下步骤： 首先进行大规模地震数据的地形分块，具体包括将层位数据和地质约束多边形数据划 分入不同地形块内并保存至外存中； 然后以每个地形块为基本单位进行曲面重建，所述曲面重建步骤包括:从对应路径下 的外存中读取地形块内部的层位数据和地质约束多边形数据，然后在地形块内构建初始网 格，在每个初始网格内进行约束Delaunay三角剖分并进行三维空间三角网恢复，得到三维 拓扑，最后对三维拓扑上的点进行曲面插值，得到三维地质曲面。2. 根据权利要求1所述的基于大规模地震数据的三维地质曲面重建方法，其特征在于， 在对地质约束多边形数据进行分块处理时，当一个地形块内只有一个地质约束多边形点 时，求出两个相邻地质约束多边形点连成的线段与地形块边界的交点，并加入对应的地形 块内的地质约束多边形数据中。3. 根据权利要求1所述的基于大规模地震数据的三维地质曲面重建方法，其特征在于， 所述方法还包括保留一份整体的封闭多边形数据，并且，在对每个地形块进行曲面重建操 作时，都要调用封闭多边形数据进行属性判断，具体地， 定义封闭多边形数据集合P: Ρ={Ρι，Ρ2，···，Ρη} 等式（1) 其中， Pi={upi，lpi，idi} 等式（2) 在等式（1)和等式(2)中，P表示封闭多边形的集合，Pi表示封闭多边形集合P中的第i个 封闭多边形，1 <m，m为封闭多边形的总个数，uPi表示第i个封闭多边形的上盘多边形点 集合，lPl表示第i个封闭多边形的下盘多边形点集合，icU表示第i个封闭多边形的根多边形 id〇4. 根据权利要求1所述的基于大规模地震数据的三维地质曲面重建方法，其特征在于， 所述方法还包括引入地质约束多边形点的索引属性，以标识地质约束多边形点在整个地质 约束多边形中的位置，并用于判断地质约束多边形点之间是否相邻，只有索引相邻的两个 地质约束多边形点才能连接形成约束边，具体地， 定义地质约束多边形点集合： VP={Vi，V2，，，_，Vn} 等式(3) 其中， Vj = {Xj，Yj，IDj，Indexj} 等式(4) 在等式(3)和等式(4)中，VP表示一个地质约束多边形上所有地质约束多边形点的集合， Vj表示地质约束多边形上第j个地质约束多边形点，l<j<n，n为一个地质约束多边形上地 质约束多边形点的总个数，Xj和Yj表示地质约束多边形点的位置，IDj表示第j个地质约束多 边形点所在的地质约束多边形在地质约束多边形集合中的ID，Indexj表示第j个地质约束 多边形点在地质约束多边形中的位置索引。5. 根据权利要求1所述的基于大规模地震数据的三维地质曲面重建方法，其特征在于， 所述方法采用离散L0D模型进行曲面重建。6. 根据权利要求1所述的基于大规模地震数据的三维地质曲面重建方法，其特征在于， 所述插值步骤包括对地形块的非边界点和边界点进行深度插值，其中， 对非边界点进行深度插值的步骤为对网格点插值，再对断层多边形点进行插值，最后 进行补偿插值； 对地形块的边界点进行深度插值的过程中，所有与当前边界点拓扑连接的点都要作为 种子点参与当前边界点的深度插值计算;并且，对边界点添加插值与否的属性，将插值后的 边界点的属性设置为已插值，对已插值的点不进行二次插值。7. 根据权利要求1所述的基于大规模地震数据的三维地质曲面重建方法，其特征在于， 所述方法读取大规模地震数据中的预定行，并对所述预定行的地震数据进行地形分块处 理，然后重复进行所述读取和地形分块操作，直至大规模地震数据被读取和处理完。8. 根据权利要求1所述的基于大规模地震数据的三维地质曲面重建方法，其特征在于， 所述对初始给定的地形块进行扩充时保证整个区域呈矩形。9. 根据权利要求1所述的基于大规模地震数据的三维地质曲面重建方法，其特征在于， 在所述保存步骤中，为每个地形块设置一个id，并采用二维数组的方式进行地形块的id索 弓丨，确定一个二维数组a[M][N]，其中所述二维数组的第一维表示地形块的X方向，二维数组 的第二维表示地形块的Y方向，二维数组中的值为地形块的i d。10. 根据权利要求1所述的基于大规模地震数据的三维地质曲面重建方法，其特征在 于，所述方法采用能量函数法计算地形块初始的分辨率，然后对计算出来的分辨率值结合 地形块内的断层分布状况进行微调，如果一个地形块内存在的断层约束多边形越多，分辨 率越高。
【文档编号】G06T17/05GK106097446SQ201610397592
【公开日】2016年11月9日
【申请日】2016年6月7日
【发明人】王颀, 赵洲, 孙健, 陆林超, 张亨, 罗文 , 唐浩, 刘璞, 范晓, 刘开元, 吴秋波, 邹文
【申请人】中国石油集团川庆钻探工程有限公司地球物理勘探公司