CN110058303B - 声波各向异性逆时偏移混合方法 - Google Patents

Abstract

本发明涉及一种声波各向异性逆时偏移混合方法，将逆时偏移成像过程分为：1、在未输入地震数据时设定观测系统以及地震震源位置来得到各向异性介质中的地震数据。采用混合伪谱法和有限差分法来求解各向异性声学波动方程。应用有限差分法时，采用高阶旋转交错网格来保证计算精度以及稳定性，且在每个时间步长只需要一次快速傅里叶变换。2、根据输入的地震数据来应用改进的激发振幅成像条件进行逆时偏移成像。其基于上述算法进行地震正传过程以及地震记录的反传过程，再进行逆时偏移成像。在大范围的计算区域中，减少对计算内存的需求。本发明提高了逆时偏移的计算效率，在大规模的复杂各向异性介质中模拟地震声学记录并进行逆时偏移成像。

Description

声波各向异性逆时偏移混合方法

技术领域

本发明属于地震勘探技术领域，具体涉及一种逆时偏移方法，特别涉及一种声波各向异性逆时偏移混合方法，适用于具有倾斜对称轴的横向各向同性介质中的声波逆时偏移。

背景技术

根据目前的国际形势与科技水平，油气能源仍是各个国家的重要战略资源并且具有重要作用。由于现阶段全球油气资源分配不均，中东地区的油田开发吸引了许多国际能源公司的关注，甚至影响了一些发达国家的政治决策。因此，全球油气能源的勘探开发迎来了前所未有的机遇与挑战。随着陆地勘探的不断深入与海洋勘探的发展，陆上可开发油气田与海底浅层油气田越来越少，所以深水油气变得愈发重要。在实际生产工作中，地震勘探是深部油气资源勘探的一种重要方法，其勘探深度大且精度高。首先对采集到的地震记录做预处理，然后经过速度分析、反演、偏移等步骤得到详细的深部构造信息。

在海洋油气田的地震勘探中，船拖着震源与拖缆中的检波器连续航行作业，地震记录由位于水面下一定深度的水听器接收。由于剪切波不能在水中传播，接收到的地震记录只含有纵波分量，类似声学信号。另一方面，海底以下的介质在复杂应力条件下通常表现出各向异性，因此针对海洋地震勘探生产过程，声学逆时偏移方法已经被引入到了各向异性介质中。

模拟声波在各向异性介质中的传播过程是声学逆时偏移的前提与基础。对各向异性数值模拟的研究在20世纪80年代以后得到飞速的发展并逐步走向完善。1986年，Thomsen提出用Thomsen系数来表征介质的各向异性。1996年，Tsvankin将沿着对称轴的剪切波速度近似为零，进而推导出了横向各向同性介质中的声波各向异性波动方程。研究发现，在横向各向同性介质中实现精确的声学数值模拟对逆时偏移(RTM)结果有非常重要的作用。

此后，许多研究都集中于横向各向同性介质中的声波传播过程。Zhou在2006年提出耦合方程，并将其从具有垂直对称轴的横向各向同性介质推广到具有倾斜对称轴的横向各向同性介质。但是，虽然沿对称轴的剪切波的速度被设置为零，但这些方法仍然将纵波与剪切波耦合在一起。为了进一步消除不需要的横波分量，2009年Etgen提出了具有垂直对称轴的横向各向同性介质中纯纵波传播的解耦方程，然后将包含复杂伪微分算子的解耦方程引入到了具有倾斜对称轴的横向各向同性介质中。然而，由于在每个时间点都需要进行多次快速傅里叶变换，计算成本很高。

声学各向异性逆时偏移方法面临的主要问题是计算效率较低。这个问题有两方面的原因，首先由于各向异性介质中的传播较为复杂，正传波场与反传波场的模拟产生了大量的计算量；第二个原因是计算成像值时，每一时刻的正传波场快照都需要保存在计算机内存中，较大区域的计算会导致内存溢出。2013年，Zhan推导出了一个包含波数项与空间导数的纯纵波波动方程并应用混合法求解。由于将波数项的数量减少为两个并且完全地去除了波场中的横波分量，这一理论在保证计算效率的基础上实现了各向异性介质中的纯纵波模拟。然而计算过程中需要三次快速傅里叶变换，计算量仍然是需要解决的问题。同时波动方程中的二阶差分项的计算精度较低，难以适应较大的时间与空间网格间隔。将这种方法应用到逆时偏移成像方法时，结合成像条件进行波场的正传与反传会显著放大这一计算效率的问题。综上所述，在具有倾斜对称轴的横向各向同性介质中的声波逆时偏移技术存在着计算效率低的缺陷。

发明内容

本发明的目的是针对上述现有技术的不足，提供一种高计算效率的声波各向异性逆时偏移混合方法。

本发明的目的是通过以下技术方案实现的：

本发明声波各向异性逆时偏移混合方法，将逆时偏移成像过程分为两个部分：第一部分是在没有输入地震数据的情况下设定观测系统以及地震震源位置来得到各向异性介质中的地震数据；第二部分是根据输入的地震数据来应用改进的激发振幅成像条件进行逆时偏移成像。

一种声波各向异性逆时偏移混合方法，具体包括以下步骤：

A、当已知多炮地震记录时，直接输入地震记录；当地震记录未知时，输入速度模型与各向异性参数模型，基于混合格式声波波动方程模拟地震数据，然后输入得到的地震记录；

B、切除输入的地震记录中的直达波部分；

C、输入圆滑后的速度模型与各向异性参数模型；

D、从输入的地震记录中抽取出单炮的地震记录进行逆时偏移。这一步骤分为两个部分，第一部分是基于混合格式声波波动方程计算出相同震源位置的正传波场。第二部分同样是基于混合格式声波波动方程，根据单炮地震记录计算由数据激发的反传波场；

E、记录正传波场中一个网格点振幅最大的时刻，将正传波场与反传波场在这一时刻的振幅值相乘，结果为该点的成像值，依次计算所有网格点的成像值。

F、提取每一单炮记录并重复步骤D-E，然后叠加每一炮的逆时偏移结果。最终输出炮集覆盖区域的逆时偏移结果。

与现有技术相比，本发明的有益效果在于：本发明声波各向异性逆时偏移混合方法，提高了计算效率；提出了采用伪谱法与旋转交错网格有限差分方法计算的声波各向异性波动方程，其中只含有一次快速傅里叶变换；在逆时偏移过程中采用OpenMP并行计算技术显著提高了计算效率；另外，采用激发振幅成像条件进行逆时偏移减少对计算内存的需求。本发明在具有倾斜对称轴的横向各向同性介质中是高效的逆时偏移成像方法。与此同时，本发明解决了以下问题：

1、在逆时偏移计算中，含有全部反射信息的地震记录与对精度要求较低的模糊的速度模型与各向异性参数模型是必要的输入数据；本发明包含了一个可以设置震源与观测系统的正演数值模拟功能，即在地震记录未知的情况下，根据输入的速度模型与各向异性模型计算多炮声波地震记录。

2、现有的涉及伪谱法的纯纵波波动方程中都包含有两个及两个以上的波数项，应用伪谱法进行数值模拟时，在每个时间节点都需要对波数项进行快速傅里叶变换。多次在时间空间域与频率波数域之间的转换产生大量的计算量，导致计算效率降低。本发明中的波动方程只含有一个需要进行快速傅里叶变换的波数项，达到伪谱法中最简单的格式，提高逆时偏移的计算效率。

3、对于求解波动方程中的差分项通常采用规则网格计算。由于常规网格方法的稳定性较低，它难以适应较大的空间与时间步长，也难以适应变化剧烈的速度与各向异性参数。本发明整理出含有一阶差分项的波动方程，但由于交叉导数项的存在难以直接应用交错网格来进行计算，所以采用了旋转交错网格。

4、逆时偏移往往需要多炮叠加形成完整成像结果。实际地震勘探生产中的勘探区域通常很大，炮数很多，全部单炮逆时偏移结果的叠加对计算机的要求很高。本发明应用激发振幅成像条件来减小对计算内存的压力，并且在加入OpenMP并行计算后，具有更快的计算速度。

附图说明

图1声波各向异性逆时偏移混合方法流程图；

图2a Dog Creek页岩的简化相速度公式的误差分析；

图2b Wills Point页岩的简化相速度公式的误差分析；

图2c Niobrara shale页岩的简化相速度公式的误差分析；

图2d Pierre shale页岩的简化相速度公式的误差分析；

图3旋转交错网格有限差分方法示意图；

图4a纵波速度标准BP2007具有倾斜对称轴的横向各向同性模型示意图；

图4b介质倾角标准BP2007具有倾斜对称轴的横向各向同性模型示意图；

图4c Thomsen参数ε标准BP2007具有倾斜对称轴的横向各向同性模型示意图；

图4d Thomsen参数δ标准BP2007具有倾斜对称轴的横向各向同性模型示意图；

图5a纵波速度标准BP2007具有倾斜对称轴的横向各向同性模型右侧断层区域速度及正演波场示意图；

图5b 2s时刻纵波波场标准BP2007具有倾斜对称轴的横向各向同性模型右侧断层区域速度及正演波场示意图；

图5c 3s时刻纵波波场标准BP2007具有倾斜对称轴的横向各向同性模型右侧断层区域速度及正演波场示意图；

图5d 4s时刻纵波波场标准BP2007具有倾斜对称轴的横向各向同性模型右侧断层区域速度及正演波场示意图；

图6标准BP2007具有倾斜对称轴的横向各向同性模型1641炮叠加逆时偏移结果。

具体实施方式

下面结合附图与实例对本发明进行进一步的详细说明

本发明声波各向异性逆时偏移混合方法能在Linux系统下的GCC平台以及Windows系统下的VC平台实现。实现过程分为两个部分，第一部分是在地震数据以及圆滑的速度模型与各向异性参数模型已知情况下计算逆时偏移结果；第二部分是在速度模型与各向异性参数模型已知、地震记录未知情况下计算地震数据；

第一部分是逆时偏移过程，也是本发明的核心。首先需要输入多炮地震记录以及模糊的速度与各向异性参数模型，在具有垂直对称轴的横向各向同性介质中需要将对称轴的倾角φ设置为90°。根据成像条件将正传波场与反传波场中的振幅值相乘得到单炮记录的逆时偏移成像结果，多炮记录的结果叠加后得到全区域的逆时偏移成像结果；在第二部分中，根据输入的速度模型与各向异性参数模型设置多个炮点以及观测系统，然后采用伪谱法与旋转交错网格有限差分方法求解各向异性声波波动方程来计算多炮地震记录；本发明采用了高效的波动方程与稳定的计算策略，对于两个部分的计算均采用并行技术，实现了具有倾斜对称轴的横向各向同性介质中的逆时偏移成像过程。

在Linux系统中需要搭建OpenMP编译库；在Windows系统中需要开启VC平台中的OpenMP语言支持。

一种声波各向异性逆时偏移混合方法，具体包括以下步骤：

A、当已知多炮地震记录时，直接输入地震记录。当地震记录未知时，输入速度模型与各向异性参数模型。设定正演数值模拟的参数，包括模型的网格大小nz×nx，空间网格间距dz、dx，时间采样间距dt，时间采样点数nt。然后，确定炮点的位置并且设置观测系统。基于混合格式声波波动方程模拟地震数据，然后为逆时偏移输入得到的地震记录。

B、在输入地震数据后，需要切除数据中的直达波部分。

C、输入的速度模型及各向异性参数模型需要进行圆滑处理。

D、从输入的地震记录中抽取出单炮的地震记录进行逆时偏移。这一步骤分为两个部分，第一部分是基于混合格式声波波动方程计算出相同震源位置的正传波场。第二部分同样是基于混合格式声波波动方程，根据单炮地震记录计算由数据激发的反传波场；

E、记录正传波场中一个网格点振幅最大的时刻，将正传波场与反传波场在这一时刻的振幅值相乘，结果为该点的成像值。此成像条件为改进的激发振幅成像条件

这种计算方法不需要保存正传波场每一时间采样点的波场快照，避免了波场重建过程，减小了对计算内存的需求。计算所有网格点的成像值后得到单炮逆时偏移结果。

F、提取每一单炮记录并重复步骤D-E，然后叠加每一炮的逆时偏移结果。最终输出炮集覆盖区域的逆时偏移结果。

进一步地，步骤A的地震记录未知情况以及步骤D的两个部分所涉及的具有倾斜对称轴的横向各向同性介质中的声波波场模拟算法如下：

准确相速度表示为：

其中，

θ是相角，V_p与V_s分别是纵波与横波沿对称轴方向的速度,ε与δ是Thomsen各向异性参数。

设置公式(1)中的横波速度为零，得相速度与纵波速度的关系式为：

波数与相角的三角函数有如下关系：

其中，(k_x,k_z)是沿着x、z轴方向的波数，ω是角频率。

补偿算子

用来提高简化相速度的精度，得到只与纵波速度有关的简化相速度公式为：

根据四种页岩实例的速度及各项异性参数，以精确相速度为参照，对比分析了公式(4)与传统方法的误差(图2a-图2d)，结果显示公式(4)具有更高的精度。

将公式(4)转换到频率波数域，其中相角的三角函数由沿两个方向的波数代替：

此时对称轴倾角为0°，公式(5)只适用于具有垂直对称轴的横向各向同性介质。具有倾斜对称轴的横向各向同性介质中的纯纵波波动方程能够通过将波数旋转倾角φ得到，即当φ为0°时公式(8)与公式(5)具有相同的格式：

即原波数(k_x,k_z)可用旋转后的波数(k_x,k_z)表示为：

代入公式(5)可得旋转倾角φ后的相速度公式：

其中，

q是辅助变量。

在方程两边同时乘波场应力P，并做反傅里叶变换，可得具有倾斜对称轴的横向各向同性介质中的波动方程：

公式(9)中只有波数项Q需要进行傅里叶变换，将其转为一阶差分格式：

其中，P需要依次对x、z两个方向连续求一阶偏导数，所以采用旋转交错网格来计算差分算子。其中密度ρ与x、z方向的速度算子定义在整网格点上，其余项如应力、Thomsen参数、辅助变量Q等均定义在半网格点上。x与

为旋转后的方向，二者互相垂直(图3)。

在旋转坐标系下，原空间一阶差分项可表示为：

其中，Δr是旋转坐标系下的网格间距，等同于原坐标系下对角线方向，其长度可表示为

旋转交错网格下的M阶精度差分格式计算方法为：

本发明采用八阶旋转交错网格差分精度，其中差分系数采用最小二乘法计算。

实施例1

根据具体硬件环境，选择在Linux系统下将OpenMP编译库进行安装并完成GCC编译平台的搭建。

选取标准BP2007具有倾斜对称轴的横向各向同性模型以及标准地震记录来进行逆时偏移测试。模型大小为78.73km×11.26km,网格大小为12596×1801，网格间距为6.25m。地震记录中共包含1641炮，时间采样点为1151个，时间采样间隔为8毫秒。地震记录为单侧照明，炮间距为50m，测线上有800个检波器，道间距为12.5m，炮检距为50m。逆时偏移中计算正传波场与反传波场采用中心频率为30Hz的雷克子波作为震源。

图4a-图4d为标准BP2007具有倾斜对称轴的横向各向同性模型全区域的速度模型与各向异性参数模型。其中浅层为水层，具有最小的速度1492m/s。模型中左侧有两个盐丘构造，具有最大的速度4352m/s，模型右侧有复杂的断层构造。

首先，测试发明的正演数值模拟功能。截取速度模型与各向异性参数模型右侧的断层部分，图5a-图5d是这一部分的速度模型以及在2s,3s,4s时刻截取的波场快照。可以看出断层区域介质速度变化剧烈并且构造极其复杂，存在若干绕射点。在整个纵波的传播过程中波前面没有发生频散现象且观察到有效的反射波。本发明计算出复杂具有倾斜对称轴的横向各向同性介质中稳定的纵波波场，在地面设置观测系统可以记录地震数据。

然后，测试本发明的核心部分，进行整个模型的逆时偏移。首先，除去所有炮的地震数据中的直达波，输入此地震记录以及圆滑后的速度模型及各向异性参数模型。在地震记录中抽取出第一炮的记录，模拟出震源在相同位置的正传波场以及由地震记录激发的反传波场，根据改进的激发振幅成像条件提取出正传波场与反传波场对应的振幅值并进行计算，最终得到单炮逆时偏移结果。

图6为根据1641炮地震记录重复上述步骤叠加得到的逆时偏移结果。其中由于地面放炮在测线范围内难以接收到陡倾角区域的反射波，所以此部分垂向同相轴不连续，并且由于盐丘高速体上部分的遮挡，下部成像结果中有白噪现象。海底地表起伏情况与实际介质吻合，盐丘构造上表面的深度与精确模型保持一致，深部横向盐丘中没有界面假象。逆时偏移结果的整体连续性好、信噪比高，褶皱地层的走向以及复杂的断层构造信息很详细。

本发明中提供地震数据的功能可以用于复杂的实际构造，并且在速度剧烈变化的情况下仍然能够提供稳定的地震数据。在逆时偏移部分，本发明根据输入的低信噪比的实际数据、圆滑的速度模型及各向异性参数模型计算出逆时偏移结果。

Claims (1)

1.一种声波各向异性逆时偏移混合方法，其特征在于，包括以下步骤：

A、当已知多炮地震记录时，直接输入地震记录；当地震记录未知时，进行地震正演，设定正演数值模拟的参数，包括模型的网格大小、空间网格间距、时间采样间距和时间采样点数，输入速度模型与各向异性参数模型，设置震源和观测系统，基于混合格式声波波动方程模拟地震数据，然后输入得到的地震记录；

B、切除输入的地震记录中的直达波部分；

C、输入圆滑后的速度模型与各向异性参数模型；

D、从输入的地震记录中抽取出单炮的地震记录进行逆时偏移；

D1、基于混合格式声波波动方程进行相同震源位置的正传波场模拟；

D2、基于混合格式声波波动方程，根据单炮地震记录进行由数据激发的反传波场模拟；

E、记录正传波场中一个网格点振幅最大的时刻，将正传波场与反传波场在这一时刻的振幅值相乘，结果为该点的成像值，依次计算所有网格点的成像值；

F、提取每一单炮记录并重复步骤D-E，然后叠加每一炮的逆时偏移结果；最终，输出炮集覆盖区域的逆时偏移结果；

所述步骤A、步骤D1以及步骤D2中，所述混合格式声波波动方程同时含有差分项和一个波数项；

当横波速度为零时，在频率波数域中的相速度与纵波速度的关系式为：

上式中ω是角频率，(k_x,k_z)是沿着x、z轴方向的波数，V_p是纵波沿对称轴方向的速度,ε和δ是Thomsen各向异性参数，补偿算子

θ是相角；

倾角φ为横向各向同性介质的倾斜对称轴的倾角，将公式(1)中的波数旋转倾角φ后为：

上式中

k_x和k_z是旋转倾角φ后沿着x、z轴方向的波数，q是辅助变量；

在公式(2)两边同时乘波场的应力P，并做反傅里叶变换为：

上式中α、β和Q为辅助变量，v_x与v_z分别为x和z方向的速度算子，t为时间，ρ为密度；

公式(3)为所述的混合格式波动方程，步骤A的地震记录模拟以及步骤D的正传波场与反传波场的模拟是通过采用伪谱法和旋转交错网格有限差分方法求解公式(3)完成。

Images