基于地质单元体的地震正演方法及装置

摘要

本发明提供了一种基于地质单元体的地震正演方法及装置，包括：将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体；根据相邻地质单元体之间的拓扑关系及公共边界，建立能够描述各个地质单元体之间拓扑关系的图形数据结构；根据地表观测系统中的炮点位置、检波点位置以及所研究的目标地质体遍历图形数据结构，计算可行的地震波射线路径及其传播时间；根据可行的射线路径及其传播时间，合成地震正演模拟记录；计算任务都在遍历的射线路径中进行，且参与计算的总网格数量极少，减小了正演计算量；采用最小二乘法能够快速准确的计算各个初始射线路径对应的可行射线路径，减少了计算量且适用于复杂地质地层模拟。

权利要求

1.一种基于地质单元体的地震正演方法，其特征在于，所述方法包括：将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体；

根据各个地质单元体的边界及相邻地质单元体之间的拓扑关系，建立各个地质单元体之间的图形数据结构；其中，各个所述地质单元体均包括多个边界；

根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，计算地震波的可行射线路径及所述可行射线路径的传播时间；其中，所述炮点位置和所述检波点位置所在的地质单元体分别对应所述可行射线路径中的下行路径的起始节点和上行路径的终止节点；

根据所述可行射线路径及所述可行射线路径的传播时间，计算地震正演模拟结果；

所述根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，计算地震波的可行射线路径及所述可行射线路径的传播时间，包括：根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，得到初始射线路径；其中，所述初始射线路径至少经过所述目标地质体包括的地质单元体；

在所述图形数据结构中查询所述初始射线路径经过的各个地质单元体的所有边界；

采用最小二乘法对所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界进行计算，得到各个所述初始射线路径对应的可行射线路径。

2.根据权利要求1所述的基于地质单元体的地震正演方法，其特征在于，将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体，包括：根据所述模拟工区对应的地层模型，获取模拟工区的弹性参数；所述弹性参数包括：岩石密度和地震波传播速度；

将岩石密度和地震波传播速度分布规律均相同的相邻地层或相邻地质体划分为地质单元体，得到多个地质单元体；

利用光滑曲线对各个地质单元体的边界进行拟合处理，得到拟合曲线；

根据所述拟合曲线，计算各个地质单元体边界上的入射点位置或出射点位置。

3.根据权利要求1所述的基于地质单元体的地震正演方法，其特征在于，所述采用最小二乘法对所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界进行计算，得到各个所述初始射线路径对应的可行射线路径，包括：根据所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界，构建射线追踪方程；

采用最小二乘法对所述射线追踪方程进行计算，得到各个所述初始射线路径对应的可行射线路径；其中，所述可行射线路径包括：各个地质单元体边界界面上入射点位置、透射/反射点位置和出射点位置。

4.根据权利要求1所述的基于地质单元体的地震正演方法，其特征在于，所述根据所述可行射线路径及所述可行射线路径的传播时间，计算地震正演模拟结果，包括：根据所述可行射线路径的传播时间，计算所述可行射线路径经过的各个地质单元体的各个边界界面上的反射系数和/或透射系数；

根据所述反射系数和/或透射系数，计算所述可行射线路径上检波点处的出射能量系数；

根据预设入射地震波能量以及所述出射能量系数，计算地震波在各个地质单元体的边界界面上的反射能量，得到地震模拟记录。

5.根据权利要求4所述的基于地质单元体的地震正演方法，其特征在于，所述根据预设入射地震波能量以及所述出射能量系数，计算地震波在各个地质单元体的边界界面上的反射能量，得到地震模拟记录，包括：将所述出射能量系数与预设地震子波进行褶积运算，得到对应射线的地震子波；

将所有所述可行射线路径上的地震子波均进行合成处理，得到地震模拟记录。

6.一种基于地质单元体的地震正演装置，其特征在于，所述装置包括：地质单元体划分模块，用于将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体；

图形数据结构建立模块，用于根据各个地质单元体的边界及相邻地质单元体之间的拓扑关系，建立各个地质单元体之间的图形数据结构；其中，各个所述地质单元体均包括多个边界；

第一计算模块，用于根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，计算地震波的可行射线路径及所述可行射线路径的传播时间；其中，所述炮点位置和所述检波点位置所在的地质单元体分别对应所述可行射线路径中的下行路径的起始节点和上行路径的终止节点；

第二计算模块，用于根据所述可行射线路径及所述可行射线路径的传播时间，计算地震正演模拟结果；

所述第一计算模块，包括：

遍历单元，用于根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，得到初始射线路径；其中，所述初始射线路径至少经过所述目标地质体包括的地质单元体；

查询单元，用于在所述图形数据结构中查询所述初始射线路径经过的各个地质单元体的所有边界；

第二计算单元，用于采用最小二乘法对所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界进行计算，得到各个所述初始射线路径对应的可行射线路径。

7.根据权利要求6所述的基于地质单元体的地震正演装置，其特征在于，所述地质单元体划分模块，包括：获取单元，用于根据所述模拟工区对应的地层模型，获取模拟工区的弹性参数；所述弹性参数包括：岩石密度和地震波传播速度；

地质单元体划分单元，用于将岩石密度和地震波传播速度分布规律均相同的相邻地层或相邻地质体划分为地质单元体，得到多个地质单元体；

拟合处理单元，用于利用光滑曲线对各个地质单元体的边界进行拟合处理，得到拟合曲线；

第一计算单元，用于根据所述拟合曲线，计算各个地质单元体边界上的入射点位置或出射点位置。

8.根据权利要求6所述的基于地质单元体的地震正演装置，其特征在于，所述第二计算单元，包括：射线追踪方程构建子单元，用于根据所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界，构建射线追踪方程；

计算子单元，用于采用最小二乘法对所述射线追踪方程进行计算，得到各个所述初始射线路径对应的可行射线路径；其中，所述可行射线路径包括：各个地质单元体边界界面上入射点位置、透射/反射点位置和出射点位置。

说明书

基于地质单元体的地震正演方法及装置

技术领域

[0001] 本发明涉及地震正演模拟技术领域，具体而言，涉及一种基于地质单元体的地震正演方法及装置。

背景技术

[0002] 地震反射正演模拟在地震勘探理论研究和数据处理方法中具有重要的作用和意义。目前主要有射线追踪和波动方程两大类方法：两者比较起来，射线追踪类方法由于进行高频近似，使得其波场特征可控，但该方法不适宜于复杂模型；而波动方程类方法是通过不同方式推演波在介质中的传播，其可用于复杂模型，但该方法波场特征复杂，计算难度较大。实际中，应用上述两类方法进行地震反射正演模拟时，均需要首先选定波场模拟区域，然后将波场模拟区域进行相应的网格剖分，如采用矩形网格剖分或者三角网格剖分。但是，上述两类方法在剖分方面会考量计算负荷和复杂地质边界的近似程度，剖分网格面积越小，复杂地质边界的拟合近似程度越高，但计算量越大，误差累计效应越严重；而剖分网格面积越大，计算量越小，但地质边界近似程度越差。

[0003] 另外，在确定的观测系统中，即确定的炮点和检波点位置下，上述两类计算方法中都会出现冗余计算：在射线追踪类方法中，需要计算若干条射线路径，然后从中选择出射点位置最接近检波器的一条路径，而其余其他路径的计算均属于冗余计算；在波动方程类方法中，需要计算各个网格点上的波场特征，然后从中选择检波器位置处的波场特征并进行记录，而其余点的波场特征计算均属于冗余计算。

[0004] 发明人在研究中发现，现有技术中的地震反射正演模拟的方法均不理想，针对这一问题，目前尚未提供较好的解决方式。

发明内容

[0005] 有鉴于此，本发明实施例的目的在于提供一种基于地质单元体的地震正演方法及装置，以快速、准确的计算地震正演模拟结果。

[0006] 第一方面，本发明实施例提供了一种基于地质单元体的地震正演方法，所述方法包括：

[0007] 将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体；

[0008] 根据各个地质单元体的边界及相邻地质单元体之间的拓扑关系，建立各个地质单元体之间的图形数据结构；其中，各个所述地质单元体均包括多个边界；

[0009] 根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，计算地震波的可行射线路径及所述可行射线路径的传播时间；其中，所述炮点位置和所述检波点位置所在的地质单元体分别对应所述可行射线路径中的下行路径的起始节点和上行路径的终止节点；

[0010] 根据所述可行射线路径及所述可行射线路径的传播时间，计算地震正演模拟结果。

[0011] 结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体，包括：

[0012] 根据所述模拟工区对应的地层模型，获取模拟工区的弹性参数；所述弹性参数包括：岩石密度和地震波传播速度；

[0013] 将岩石密度和地震波传播速度分布规律均相同的相邻地层或相邻地质体划分为地质单元体，得到多个地质单元体；

[0014] 利用光滑曲线对各个地质单元体的边界进行拟合处理，得到拟合曲线；

[0015] 根据所述拟合曲线，计算各个地质单元体边界上的入射点位置或出射点位置。

[0016] 结合第一方面的第一种可能的实施方式，本发明实施例提供了第一方面的第二种可能的实施方式，其中，所述根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，计算地震波的可行射线路径及所述可行射线路径的传播时间，包括：

[0017] 根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，得到初始射线路径；其中，所述初始射线路径至少经过所述目标地质体包括的地质单元体；

[0018] 在所述图形数据结构中查询所述初始射线路径经过的各个地质单元体的所有边界；

[0019] 采用最小二乘法对所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界进行计算，得到各个所述初始射线路径对应的可行射线路径。

[0020] 结合第一方面的第二种可能的实施方式，本发明实施例提供了第一方面的第三种可能的实施方式，其中，所述采用最小二乘法对所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界进行计算，得到各个所述初始射线路径对应的可行射线路径，包括：

[0021] 根据所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界，构建射线追踪方程；

[0022] 采用最小二乘法对所述射线追踪方程进行计算，得到各个所述初始射线路径对应的可行射线路径；其中，所述可行射线路径包括：各个地质单元体边界界面上入射点位置、透射/反射点位置和出射点位置。

[0023] 结合第一方面，本发明实施例提供了第一方面的第四种可能的实施方式，其中，所述根据所述可行射线路径及所述可行射线路径的传播时间，计算地震正演模拟结果，包括：

[0024] 根据所述可行射线路径的传播时间，计算所述可行射线路径经过的各个地质单元体的各个边界界面上的反射系数和/或透射系数；

[0025] 根据所述反射系数和/或透射系数，计算所述可行射线路径上检波点处的出射能量系数；

[0026] 根据预设入射地震波能量以及所述出射能量系数，计算地震波在各个地质单元体的边界界面上的反射能量，得到地震模拟记录。

[0027] 结合第一方面的第四种可能的实施方式，本发明实施例提供了第一方面的第五种可能的实施方式，其中，所述根据预设入射地震波能量以及所述出射能量系数，计算地震波在各个地质单元体的边界界面上的反射能量，得到地震模拟记录，包括：

[0028] 将所述出射能量系数与预设地震子波进行褶积运算，得到对应射线的地震子波；

[0029] 将所有所述可行射线路径上的地震子波均进行合成处理，得到地震模拟记录。

[0030] 第二方面，本发明实施例还提供了一种基于地质单元体的地震正演装置，所述装置包括：

[0031] 地质单元体划分模块，用于将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体；

[0032] 图形数据结构建立模块，用于根据各个地质单元体的边界及相邻地质单元体之间的拓扑关系，建立各个地质单元体之间的图形数据结构；其中，各个所述地质单元体均包括多个边界；

[0033] 第一计算模块，用于根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，计算地震波的可行射线路径及所述可行射线路径的传播时间；其中，所述炮点位置和所述检波点位置所在的地质单元体分别对应所述可行射线路径中的下行路径的起始节点和上行路径的终止节点；

[0034] 第二计算模块，用于根据所述可行射线路径及所述可行射线路径的传播时间，计算地震正演模拟结果。

[0035] 结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方式，其中，所述地质单元体划分模块，包括：

[0036] 获取单元，用于根据所述模拟工区对应的地层模型，获取模拟工区的弹性参数；所述弹性参数包括：岩石密度和地震波传播速度；

[0037] 地质单元体划分单元，用于将岩石密度和地震波传播速度分布规律均相同的相邻地层或相邻地质体划分为地质单元体，得到多个地质单元体；

[0038] 拟合处理单元，用于利用光滑曲线对各个地质单元体的边界进行拟合处理，得到拟合曲线；

[0039] 第一计算单元，用于根据所述拟合曲线，计算各个地质单元体边界上的入射点位置或出射点位置。

[0040] 结合第二方面的第一种可能的实施方式，本发明实施例提供了第二方面的第二种可能的实施方式，其中，所述第一计算模块，包括：

[0041] 遍历单元，用于根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，得到初始射线路径；其中，所述初始射线路径至少经过所述目标地质体包括的地质单元体；

[0042] 查询单元，用于在所述图形数据结构中查询所述初始射线路径经过的各个地质单元体的所有边界；

[0043] 第二计算单元，用于采用最小二乘法对所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界进行计算，得到各个所述初始射线路径对应的可行射线路径。

[0044] 结合第二方面的第二种可能的实施方式，本发明实施例提供了第二方面的第三种可能的实施方式，其中，所述第二计算单元，包括：

[0045] 射线追踪方程构建子单元，用于根据所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界，构建射线追踪方程；

[0046] 计算子单元，用于采用最小二乘法对所述射线追踪方程进行计算，得到各个所述初始射线路径对应的可行射线路径；其中，所述可行射线路径包括：各个地质单元体边界界面上入射点位置、透射/反射点位置和出射点位置。

[0047] 本发明实施例提供的一种基于地质单元体的地震正演方法及装置，采用地质单元体对模拟工区的地层模型进行网格剖分处理，将弹性参数相同的地层或地质体划分一个地质单元体，同时，采用图形数据结构表示根据各个地质单元体的边界及相邻地质单元体之间的拓扑关系，并通过遍历上述图形数据结构来计算地震波的可行射线路径及可行射线路径的传播时间；根据可行射线路径及可行射线路径的传播时间，计算地震正演模拟结果，与现有技术中的射线追踪和波动方程两大类方法均存在计算量越大且误差累计效应严重相比，其具有计算量小、计算精度高、且适用于复杂地质构造的优点，具体体现在如下方面：

[0048] (1)由于地质单元体的面积远比现有技术中的矩形网格和三角网格要大许多，因此，参与计算的总网格数量极少，减小了正演计算量。(2)采用图形数据结构管理网格节点，并通过图的遍历找出所有可能的射线路径。后续的计算任务都限制在遍历出来的射线路径中进行，而不是在所有地质单元体范围内进行，同样减少了计算量。

[0049] 进一步的，本发明实施例提供的一种基于地质单元体的地震正演方法及装置，采用最小二乘法能够快速准确的计算各个初始射线路径对应的可行射线路径，大大减少了计算量且适用于复杂地质地层模拟。

[0050] 为使本发明的上述目的、特征和优点能更明显易懂，下文特举一个实施例，并配合所附附图，作详细说明如下。

附图说明

[0051] 为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

[0052] 图1示出了本发明实施例所提供的一种基于地质单元体的地震正演方法的流程图；

[0053] 图2中左图示出了地质单元体网格剖分模型示意图；图2右图示出了地质单元体网格剖分后得到图形数据结构的示意图；

[0054] 图3中示出了用于测试最小二乘法射线追踪计算量的模型示例的示意图(距离单位：米)；其中，图3中，星形代表炮点，位于原点；三角形代表检波点，间距100米；

[0055] 图4示出了将将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体的流程图；

[0056] 图5示出了根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，计算地震波的可行射线路径及所述可行射线路径的传播时间的流程图；

[0057] 图6示出了根据所述可行射线路径及所述可行射线路径的传播时间，计算地震正演模拟结果的流程图；

[0058] 图7示出了本发明实施例所提供的一种基于地质单元体的地震正演装置的结构示意图；

[0059] 图8示出了本发明实施例所提供的一种基于地质单元体的地震正演装置中地质单元体划分模块和第一计算模块的结构示意图；

[0060] 图9示出了本发明实施例所提供的一种基于地质单元体的地震正演装置中第二计算单元和第二计算模块的结构示意图；

[0061] 图10示出了本发明实施例所提供的一种基于地质单元体的地震正演装置中第五计算单元的结构示意图。

[0062] 主要标号说明：

[0063] 10、地质单元体划分模块；20、图形数据结构建立模块；30、第一计算模块；40、第二计算模块；101、获取单元；102、地质单元体划分单元；103、拟合处理单元；104、第一计算单元；301、遍历单元；302、查询单元；303、第二计算单元；3031、射线追踪方程构建子单元；3032、计算子单元；401、第三计算单元；402、第四计算单元；403、第五计算单元；4031、褶积运算子单元；4032、合成处理子单元。

具体实施方式

[0064] 为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0065] 现有技术中，通常是通过射线追踪和波动方程两大类方法计算地震反射正演模拟的方法，在应用上述两大类方法进行地震反射正演模拟时，均是采用矩形网格剖分或者三角网格剖分将预先选定的波场模拟区域进行相应的网格剖分。但是，上述两种网格剖分方法在剖分方面会考量计算负荷和复杂地质边界的近似程度，剖分网格面积越小，复杂地质边界的拟合近似程度越高，但计算量越大，误差累计效应越严重；而剖分网格面积越大，计算量越小，但地质边界近似程度越差。

[0066] 另外，在确定的观测系统中，即确定的炮点和检波点位置下，上述两类计算方法中都会出现冗余计算：射线追踪中需要计算若干条射线路径，但只从中选择出射点位置最接近检波器的那条路径，其余属于冗余计算；波动方程类方法中需要计算各个网格点上的波场特征，但只记录检波器位置处的波场，其余属于冗余计算。

[0067] 基于上述两类方法的缺陷，本发明实施例提供了一种基于地质单元体的地震正演方法及装置，该方法包括三项技术措施：(一)基于地质单元体进行网格剖分，地质单元体的边界利用连续光滑函数进行拟合。(二)利用图形数据结构来管理地质单元体的相邻关系，通过图的遍历算法可以得到所有可能的射线路径。(三)在选定的射线路径上，利用最小二乘法来快速计算较精确的射线路径。本发明技术具有计算速度快，计算误差小，可定量分析地震波能量在各个地质单元体内的传播情况。下面通过实施例进行描述。

[0068] 参考图1，本发明实施例提供了一种基于地质单元体的地震正演方法，所述方法包括：

[0069] S101、将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体。

[0070] 上述弹性参数可以为地震波传播速度、岩石密度。在本发明实施例的基于地质单元体的地震正演方法中，将地震波传播速度、岩石密度等属性相同的各个地质单元作为正演模拟计算的基本单元，即按照地质单元体来划分研究区内的地层模型。而待计算的地震射线在各个地质单元体内部只需要计算一次，地质单元体边界上的入射点和出射点准确位置直接按照光滑曲线拟合方程进行计算。

[0071] 其中，划分的上述地质单元体的边界可根据实际情况选择不同类型的曲线进行拟合，如果地质单元体边界复杂，可以选择分段光滑曲线拟合，该种拟合方式比采用矩形或三角形网格剖分更能精确刻画地质单元体边界。

[0072] S102、根据各个地质单元体的边界及相邻地质单元体之间的拓扑关系，建立各个地质单元体之间的图形数据结构；其中，各个所述地质单元体均包括多个边界。

[0073] 现有技术中对地层模型进行网格剖分采用的矩形网格剖分方法或者三角网格剖分方法，得到的各个网格的相邻网格个数是有规律的，而本发明实施例中按照地质单元体对地层模型进行网格剖分以后，各个网格的相邻网格个数不再是有规律的，其是与地质单元体的空间尺寸和地质复杂程度相关。为了能便捷的确定每次射线追踪的路径经过哪些地质单元体，本发明实施例采用图形数据结构管理地质单元体相互之间的拓扑关系，即根据各个地质单元体的边界及相邻地质单元体之间的拓扑关系，建立各个地质单元体之间的图形数据结构，在射线追踪过程中，通过图的遍历来构建下行路径和上行路径。

[0074] 本发明实施例中，下行路径的起始节点为地表下最浅的地质单元体，只要该地质单元体的边界与地表重合，则该地质单元体可以成为下行路径的起始节点。上行路径的终止节点同样是地表下最浅的地质单元体，该地质单元体的边界与地表重合。下行路径的终止节点在地下某个指定层或者是研究区的最深层，下行路径的终止节点是上行路径的起始节点。

[0075] 如图2所示，图2中左图展示的是一个简单的基于地质单元体剖分的地层模型示意图，从图2左图中可以清楚看到各个地质单元体之间的相互接壤关系，同时各个地质单元体的边界采用光滑的曲线进行了较好的拟合。图2右图展示的是该地层模型所对应的地质单元体相邻关系的图形数据结构，通过该图形数据结构可以快捷地查询到各个地质单元地之间的是否相邻。

[0076] 其中，建立的上述地质单元体的图形数据结构，表征各相邻地质单元体的拓扑关系。

[0077] 对于所有的地质单元体，顺次按照其各条边界找到相邻的地质单元体，在图形数据结构中将相邻的地质单元体连接起来，如图2中右图所示，可以快捷地查询到所有可能的射线路径。

[0078] 表1本发明方法在图1模型中的拓扑关系表

[0079]

[0080] 另外，建立所有边界及其对应的地质单元体的拓扑关系表，描述与其相邻的地质单元体之间的关系。按照一定的规则确定所有的边界序号和节点序号，图2左图所示的序号命名规则是按照地质单元体的先后顺序进行，同时边界的走向是围绕地质单元体顺时针方向。如此，容易确定各个地质单元体的拓扑关系，其中，除了命名的地质单元体以外，还增加描述了模型的顶边界(即地表)、右边界、底边界和左边界。表1是图2左图所示的模型中地质单元体对应的拓扑关系表，通过该表可以查询到所有射线路径上各条边界及其对应节点，便于后续的最小二乘法计算。

[0081] S103、根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，计算地震波的可行射线路径及所述可行射线路径的传播时间；其中，所述炮点位置和所述检波点位置所在的地质单元体分别对应所述可行射线路径中的下行路径的起始节点和上行路径的终止节点。

[0082] 在本发明实施例的基于地质单元体的地震正演计算方法中，分别选择炮点位置和检波点位置所在的地质单元体作为下行路径的起始节点和上行路径的终止节点，遍历构建的图形数据结构来计算所有可能的射线路径，然后只针对各条射线路径进行射线追踪正演模拟计算，可以较大的减少正演计算工作量。由于图形数据结构中包括相邻地质单元体之间的拓扑关系，因此通过遍历图形数据结构的方法，能够通过极少的计算量，替代了目前常规射线追踪正演中试射追踪的巨大计算量。

[0083] 另外，采用图形数据结构管理，将非常方便地针对地质单元体进行射线追踪；如在固定炮点和检波点位置后，只要沿下行追踪到该节点，再从该节点上行追踪到检波点所在的终止节点即可。

[0084] 同时，采用图遍历的方法，可以找出所有可能的射线路径，较好的解决目前射线追踪算法中多条路径情况下的路径遗漏问题。对于找到的所有可能的射线路径，本发明实施例中采用最小二乘法进行射线追踪计算，以计算每个可能的射线路径对应的可行射线路径及可行射线路的传播时间。

[0085] S104、根据所述可行射线路径及所述可行射线路径的传播时间，计算地震正演模拟结果。

[0086] 在射线追踪的过程中，对于可行射线路径上的每个边界界面依据佐普里兹方程(Zoeppritz's equafion)可以比较精确地研究地震波能量的变化规律，同时模拟比较保真的地震波场数据。

[0087] 在本发明实施例的整个正演计算过程中，误差较小，同时可以较精确的模拟各个地质单元体的边界，且能够计算出精确的入射角和反射角，因此，通过Zoepprittz方程可以计算出较精确的反射系数和透射系数，根据上述反射系数、透射系数以及预先设定的入射地震波能量，可以链式地计算地震波能量在各个地质单元体的边界界面上的反射能量和透射能量，以研究地震波场能量变化，最终在地表上能够获得较精确的地震波数据(该地震波数据实际中为地震波反射能量)。因此，采用本发明实施例的正演计算方法，不仅能够较精确地研究地震波场的运动学特征，而且也能够较精确地研究地震波场的动力学特征。

[0088] 在本发明实施例提供的基于地质单元体的地震正演方法中，由于采用地质单元体进行网格剖分，采用图形数据结构方式管理相邻之间的地质单元体并通过图的遍历得出所有可能的射线路径，然后采用最小二乘法进行射线追踪计算，其与现有技术相比，具有计算量小、计算精度高、能适用于复杂地质构造的优点：

[0089] 其一：计算量小。在本发明方法中采用了三项技术措施，极大地减少正演计算量。(1)采用地质单元体网格剖分方法。将地震波传播速度、岩石密度等属性相同的地质单元体作为一个计算网格，其面积远比目前的矩形网格和三角网格要大得许多，因此总的计算网格数量极少，则正演计算量就小。(2)采用图形数据结构管理网格节点，并通过图的遍历找出所有可能的射线路径。后续的计算任务都限制在遍历出来的射线路径中进行，而不是在所有地质单元体范围内进行，减少计算量。(3)采用最小二乘法快速计算准确的射线路径处理。目前受制于较小的试射入射角增量的限制，试射追踪次数也非常巨大，例如在图3所示的简单模型中，设定射线追踪距离误差为0.01米，最小二乘法计算大约4、5次，而按相同的第一层入射角随炮检距的变化率，如果扫描角度范围为30°，则常规的扫描法射线追踪需要试射300000次，采用二分法计算也要18次(详见表2)。

[0090] 表2本发明方法在图3模型参数下第三层底界面的射线追踪计算效率分析[0091]

[0092] 其二：计算精度高。(1)采用地质单元体网格剖分方法，计算网格数量少，误差累计效应小，计算精度高。(2)采用图形数据结构管理网格节点，可以遍历出所有可能的射线路径，解决了多路经中路径遗漏问题，计算精度高。(3)采用最小二乘法，解决了因角度扫描射线追踪中为了平衡计算量而选择较大角度增量所产生的计算误差问题。(4)复杂边界地质体的边界采用光滑曲线拟合表示，在射线追踪过程中能依据拟合曲线直接求解入射点和出射点精确位置，计算精度高。

[0093] 其三：能适用于复杂地质构造模型。(1)采用地质单元体网格剖分，使得复杂地质构造模型中网格剖分简单明了，同时能够直接赋予计算网格地质地震属性和意义，高精度的刻画地质构造单元及其相互关系，还能直观定量研究特殊地质构造中的地震响应。(2)采用地质单元网格剖分方法，对于复杂边界的地质体可以通过光滑曲线进行拟合，确保高精度地描述和模拟复杂地质体。(3)采用图形数据结构来管理计算网格，同时利用图的遍历方法可以获得全部可能的射线路径，因此在复杂地质构造中不会遗漏每一条射线路径。

[0094] 进一步的，参考图4，步骤101中将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体，首先需要分析所要模拟的工区的地层结构，然后按照相同弹性参数进行地质单元体的划分，即按照避免多解性、尽量单一解的原则将全工区划分成若干地质，上述划分地质单元体的具体方法包括：

[0095] S201、根据所述模拟工区对应的地层模型，获取模拟工区的弹性参数；所述弹性参数包括：岩石密度和地震波传播速度。

[0096] S202、将岩石密度和地震波传播速度分布规律均相同的相邻地层或相邻地质体划分为地质单元体，得到多个地质单元体。

[0097] S203、利用光滑曲线对各个地质单元体的边界进行拟合处理，得到拟合曲线。

[0098] S204、根据所述拟合曲线，计算各个地质单元体边界上的入射点位置或出射点位置。

[0099] 结合上述步骤201～步骤204，首先分析所要模拟的工区的地层结构，将全工区的地层结构按照分析后的相同弹性参数划分成若干个地质单元体，这样的划分方式避免了多解性、使划分的结果尽量满足单一解的原则。

[0100] 具体的，上述划分的方式如下，首先分析模拟工区的地层参数，将相同的岩石密度和地震波传播速度分布规律的地层划分为地质单元体，划分的单元体内部的密度和速度可以是均匀的，也可以随深度呈规律的线性变化或者指数变化。通常，可以将地质单元体表示为W(V,D,βv,βd)，其中，V和βv分别表示背景速度以及速度随深度变化的参数；D和βd分别表示背景密度以及密度随深度变化的参数，在均匀各向同性介质中，βv和βd都等于0。

[0101] 另外，依据避免多解性、尽量单一解的原则，需要对地质单元体的边界进行处理，在进行地质单元体的边界处理时，首先取凸多边形，同时还要进行分段拟合，保障后续的最小二乘法求解。具体的，在进行分段拟合时，可以根据地质单元体边界的复杂程度确定拟合方程的次数，一般表示为Z＝F(x,y)。

[0102] 本发明实施例中，根据上述避免多解性、尽量单一解的两条原则确定地质单元体剖分方法，并给各个地质单元体进行编号，同时确定各个地质单元体的弹性参数及其各条边界的拟合曲线。在本发明实施例的图2中，对所示的地层模型进行了网格剖分、边界划分和节点标识，其实施过程参考了上述地质单元体划分的原则。

[0103] 本发明实施例中，通过遍历图形数据结构的方法获取所有可能的射线路径，然后根据这些射线路径的边界及这些射线路径经过的地质单元体的弹性参数，计算这些射线路径中的可行射线路径以及可行射线路径匹配的时间，对应的，参考图5，上述步骤103的具体计算方法，包括：

[0104] S301、根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历所述图形数据结构，得到初始射线路径；其中，所述初始射线路径至少经过所述目标地质体包括的地质单元体。

[0105] 上述初始射线路径即通过遍历图形数据结构的方法获取所有可能的射线路径。

[0106] S302、在所述图形数据结构中查询所述初始射线路径经过的各个地质单元体的所有边界。

[0107] S303、采用最小二乘法对所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界进行计算，得到各个所述初始射线路径对应的可行射线路径。

[0108] 本步骤中，计算各个初始射线路径对应的可行射线路径的方法，具体包括：

[0109] 根据所述初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界，构建射线追踪方程；

[0110] 采用最小二乘法对所述射线追踪方程进行计算，得到各个所述初始射线路径对应的可行射线路径；其中，所述可行射线路径包括：各个地质单元体边界界面上入射点位置、透射/反射点位置和出射点位置。

[0111] 结合上述步骤301～步骤303，在获取的所有可能的射线路径中，采用最小二乘法，快速计算每个可能的射线路径中实际可行的射线路径并得出该实际可行的射线路径的射线传播时间。

[0112] 具体的，在地表所布置的观测系统中选定炮点位置和检波点位置，并确定研究的目标地层，然后从地质单元体的图形数据结构中查询到所有可能的射线路径，同时从边界对应的拓扑关系表(该拓扑关系表包括边界表)中查询到相应的边界，然后将各个可能的射线路径对应的地质单元体的弹性参数方程以及边界拟合方程代入到最小二乘法计算模块中，由最小二乘法计算模块进行计算，最终根据设定的误差选择较精确的射线路径，该较精确的射线路径即每个可能的射线路径对应的可行射线路径。

[0113] 实际中，根据炮点位置、检波点位置以及研究的目标地层，假设选定了其中的一条可能的射线路径，选择的该射线路径顺次通过了N个地质单元体，本发明实施例中允许射线多次穿过同一个地质单元体，该射线路径中的射线将旅行2N-1条次边界，在上述射线追踪计算中只有一个未知量，即射线在炮点位置的入射角。本发明实施例中，首先选择合适的入射角初值，在选择了合适的入射角初值后，对2N-1条次边界和N个次地质单元体所对应的参数构成的射线追踪方程进行最小二乘法计算，就能够得到合适的计算结果。

[0114] 当然，并不是所有的可能路径都能得到合适的计算结果，有的研究目标地质单元体处在射线的盲区。比如在图2左图所示的模型中，如果炮点和检波点都位于3号地质单元体的B10地表边界上，而研究目标是2号和5号地质单元体的交界B8，则从图2右图(即图形数据结构中)选定的经过3-2-5-2-3地质单元体的射线路径，通过最小二乘法的计算结果是不正确的，因为边界B8正好处在射线路径的盲区。如此类似的盲区问题，本发明实施例中可以通过最小二乘法对图形数据结构中的所有的可能的射线路径进行预分析处理，排出盲区，减少计算量。

[0115] 目前的射线追踪算法中，一般采用入射角按增量进行角度扫描，上述试射追踪如果是在简单的水平或单斜地层模型中，可以采用二分法的策略来减少入射角试射追踪次数。但是，在复杂地质构造模型中，二分法难以实施，只能采用入射角增量试射的办法，这会因为选择符合精度要求的较小的角度增量而造成巨大的计算量。而本发明实施例中采用最小二乘法的方法进行射线追踪，能够快速地通过计算获得精度较高的结果，极大的减少了计算量。

[0116] 本发明实施例中，参考图6，上述步骤104中，所述根据所述可行射线路径及所述可行射线路径的传播时间，计算地震正演模拟结果，具体包括：

[0117] S401、根据所述可行射线路径的传播时间，计算所述可行射线路径经过的各个地质单元体的各个边界界面上的反射系数和/或透射系数。

[0118] S402、根据所述反射系数和/或透射系数，计算所述可行射线路径上检波点处的出射能量系数。

[0119] 本发明实施例中，作为第一种具体的实施方式，根据所述反射系数，计算所述可行射线路径上检波点处的出射能量系数；作为第二种具体的实施方式，根据所述透射系数，计算所述可行射线路径上检波点处的出射能量系数；作为第三种具体的实施方式，根据所述反射系数和所述透射系数，计算所述可行射线路径上检波点处的出射能量系数；

[0120] S403、根据预设入射地震波能量以及所述出射能量系数，计算地震波在各个地质单元体的边界界面上的反射能量，得到地震模拟记录。

[0121] 结合上述步骤401～步骤403，本步骤是选择合适的射线追踪结果，得到正演模拟结果。

[0122] 如果是研究运动学问题，就选择一条旅行时间最短的射线路径，将其旅行时间作为最终的地震模拟结果；如果是研究动力学问题，可以在所有正确的射线路径上应用Zoepprittz方程继续计算各条界面上的透射系数或者反射系数，链式计算得到各条射线路径上检波点处的出射能量系数，再将出射能量系数与设定的地震子波褶积运算得到对应射线的地震子波，最后将所有正确路径上的地震子波合成得到最终的地震模拟记录。

[0123] 其中，上述链式计算即计算射线路径在当前界面的反射能量系数，同时计算该射线路径在当前界面的透射能量系数，然后，根据投射能量系数，在计算其在下一个界面的反射能量系数，根据当前射线路径的发射和透射次数进行上述反射能量系数和透射能量系数的计算，最后，得到所有到达地表的反射能量系数，并据此计算到达地表的地震波能量。

[0124] 本发明实施例还提供了一种基于地质单元体的地震正演装置，所述装置用于执行上述基于地质单元体的地震正演方法，参考图7，所述装置包括：

[0125] 地质单元体划分模块10，用于将模拟工区中弹性参数相同的相邻地层或相邻地质体划分成同一个地质单元体，得到多个地质单元体；

[0126] 图形数据结构建立模块20，用于根据各个地质单元体的边界及相邻地质单元体之间的拓扑关系，建立各个地质单元体之间的图形数据结构；其中，各个地质单元体均包括多个边界；

[0127] 第一计算模块30，用于根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历图形数据结构，计算地震波的可行射线路径及可行射线路径的传播时间；其中，炮点位置和检波点位置所在的地质单元体分别对应可行射线路径中的下行路径的起始节点和上行路径的终止节点；

[0128] 第二计算模块40，用于根据可行射线路径及可行射线路径的传播时间，计算地震正演模拟结果。

[0129] 进一步的，参考图8，地质单元体划分模块10，包括：

[0130] 获取单元101，用于根据所述模拟工区对应的地层模型，获取模拟工区的弹性参数；所述弹性参数包括：岩石密度和地震波传播速度；

[0131] 地质单元体划分单元102，用于将岩石密度和地震波传播速度分布规律均相同的相邻地层或相邻地质体划分为地质单元体，得到多个地质单元体；

[0132] 拟合处理单元103，用于利用光滑曲线对各个地质单元体的边界进行拟合处理，得到拟合曲线；

[0133] 第一计算单元104，用于根据所述拟合曲线，计算各个地质单元体边界上的入射点位置或出射点位置。

[0134] 进一步的，参考图8，第一计算模块30，包括：

[0135] 遍历单元301，用于根据地表观测系统中的炮点位置、检波点位置以及目标地质体遍历图形数据结构，得到初始射线路径；其中，初始射线路径至少经过目标地质体包括的地质单元体；

[0136] 查询单元302，用于在图形数据结构中查询初始射线路径经过的各个地质单元体的所有边界；

[0137] 第二计算单元303，用于采用最小二乘法对初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界进行计算，得到各个初始射线路径对应的可行射线路径。

[0138] 进一步的，参考图9，第二计算单元303，包括：

[0139] 射线追踪方程构建子单元3031，用于根据初始射线路径经过的各个地质单元体的弹性参数以及查询的所有边界，构建射线追踪方程；

[0140] 计算子单元3032，用于采用最小二乘法对射线追踪方程进行计算，得到各个初始射线路径对应的可行射线路径；其中，可行射线路径包括：各个地质单元体边界界面上入射点位置、透射/反射点位置和出射点位置

[0141] 进一步的，参考图9，第二计算模块40，包括：

[0142] 第三计算单元401，用于根据可行射线路径的传播时间，计算可行射线路径经过的各个地质单元体的各个边界界面上的反射系数和/或透射系数；

[0143] 第四计算单元402，用于根据反射系数和/或透射系数，计算可行射线路径上检波点处的出射能量系数；

[0144] 第五计算单元403，用于根据预设入射地震波能量以及出射能量系数，计算地震波在各个地质单元体的边界界面上的反射能量，得到地震模拟记录。

[0145] 进一步的，参考图10，第五计算单元403，包括：

[0146] 褶积运算子单元4031，用于将出射能量系数与预设地震子波进行褶积运算，得到对应射线的地震子波；

[0147] 合成处理子单元4032，用于将所有可行射线路径上的地震子波均进行合成处理，得到地震模拟记录。

[0148] 本发明实施例提供的一种基于地质单元体的地震正演装置，采用地质单元体对模拟工区的地层模型进行网关剖分处理，将弹性参数相同的地层或地质体划分一个地质单元体，同时，采用图形数据结构表示根据各个地质单元体的边界及相邻地质单元体之间的拓扑关系，并通过遍历上述图形数据结构来计算地震波的可行射线路径及可行射线路径的传播时间；根据可行射线路径及可行射线路径的传播时间，计算地震正演模拟结果，与现有技术中的射线追踪和波动方程两大类方法均存在计算量越大且误差累计效应严重相比，其具有计算量小、计算精度高、且适用于复杂地质构造的优点，具体体现在如下方面：

[0149] (1)由于地质单元体的面积远比现有技术中的矩形网格和三角网格要大许多，因此，参与计算的总网格数量极少，减小了正演计算量。(2)采用图形数据结构管理网格节点，并通过图的遍历找出所有可能的射线路径。后续的计算任务都限制在遍历出来的射线路径中进行，而不是在所有地质单元体范围内进行，同样减少了计算量。

[0150] 进一步的，本发明实施例提供的一种基于地质单元体的地震正演装置，采用最小二乘法能够快速准确的计算各个初始射线路径对应的可行射线路径，大大减少了计算量且适用于复杂地质地层模拟。

[0151] 本发明实施例所提供的基于地质单元体的地震正演装置可以为设备上的特定硬件或者安装于设备上的软件或固件等。本发明实施例所提供的装置，其实现原理及产生的技术效果和前述方法实施例相同，为简要描述，装置实施例部分未提及之处，可参考前述方法实施例中相应内容。所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，前述描述的系统、装置和单元的具体工作过程，均可以参考上述方法实施例中的对应过程，在此不再赘述。

[0152] 在本发明所提供的实施例中，应该理解到，所揭露装置和方法，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

[0153] 所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

[0154] 另外，在本发明提供的实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。

[0155] 所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

[0156] 应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释，此外，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。

[0157] 最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。