CN105980887A - 减轻拟声波波动传播算子中的不稳定性的系统和方法 - Google Patents

Abstract

描述了一种方法，该方法包括接收与多个源和接收器位置对应的一个或多个地震测量，并且提供用于具有非均匀倾斜的对称轴的地质介质的地球模型。该地球模型包括对称轴方向上的非零剪切速度。所述一个或多个地震测量是根据所述地球模型和能量守恒拟声波方程组在多个时间步长上传播的。所述能量守恒拟声波方程组描述一个或多个地震波场并且是通过对能量非守恒拟声波方程组中的一个或多个地震波场的一个或多个导数项进行近似而从所述能量非守恒拟声波方程组导出的。

Description

减轻拟声波波动传播算子中的不稳定性的系统和方法

相关申请

本申请涉及2013年12月30日提交的题为“SYSTEM ANDMETHOD OF MITIGATING INSTABILITIES IN APSEUDOACOUSTIC WAVE PROPAGATOR”的美国专利申请No.14/143571(代理人案卷号no.T-9578)，通过引用将其全部内容并入本文。

技术领域

所公开的实施例一般涉及用于地球物理勘探期间的地震数据采集、处理和解释的技术，并且更具体地涉及用于倾斜的横向各向同性介质的拟声波波动传播算子。

背景技术

地震勘探涉及为了烃类沉积物而勘测地下地质介质。一些勘测被称为“海洋”勘测，因为它们是在海洋环境中进行的。但是，“海洋”勘测不仅可以在咸水环境中进行，而且还可以在淡水和微咸水中进行。在一种类型的海洋勘测中，被称为“拖曳阵列”勘测，包含地震传感器的拖缆和源的阵列被拖曳在勘测船后面。

勘测通常涉及将(一个或多个)地震源和(一个或多个)地震传感器部署在预定位置处。源生成地震波，其传播进地质介质，沿其路径造成压力改变和振动。地质介质的物理性质的变化改变地震波，诸如它们的传播方向和其它性质。部分地震波到达地震传感器。一些地震传感器对压力改变敏感(水听器)，其它对粒子运动敏感(例如，地震检波器)，并且工业勘测可以部署仅一种类型的传感器或部署两者。响应于检测到的地震波，传感器生成对应的电信号并且将它们记录在存储介质中作为地震数据。

地震数据的分析可以被执行以将地震数据处理成地质介质的图像。逆时偏移(RTM)在时间上向前将源位置处的波场传播进地质介质并且在时间上向后将接收器位置处的记录的波场传播进地质介质，并且接着关联两种类型的波场以形成地质介质的图像。为了减少RTM的计算成本，同时仍然允许各向异性，构造与使用完全弹性系统相比不那么计算昂贵的微分方程拟声波系统。一些传统方式通过将S波速度设定为零，并且消除微分方程拟声波系统中的一个或多个对应项，来修改针对完全弹性系统或拟声波近似(例如，2×2二阶微分方程拟声波系统)的波散(dispersion)关系。不幸的是，因为在这种近似下可能不满足物理守恒定律，所以这样做可以导致引入不稳定性。这些不稳定性由随时间增长的稳态噪声引起，从而使得难以对地质介质中的较深层或复杂地质结构进行成像。此外，已经发现即使当S波速度被人为地设定为大于零时，这些不稳定性也可以保持。

发明内容

因此，需要减轻或消除这种不稳定性的波传播算子。根据一些实施例，在具有一个或多个处理器和存储器的电子设备处执行方法。该方法包括提供用于地质介质的地球模型。该地质介质具有非均匀倾斜的对称轴，并且针对地质介质内的至少一个位置子集该地球模型包括对称轴方向上的非零剪切速度。该电子设备使用根据地球模型应用的二阶拟声波方程组处理一个或多个地震测量。该拟声波方程组包括描述一个或多个地震波场的第一方程和描述一个或多个地震波场的第二方程。该处理包括根据拟声波方程组在多个时间步长上传播所述一个或多个地震波场，确定多个时间步长中的相应时间步长是否满足预定准则，并且根据相应时间步长满足预定准则的确定，应用与该地球模型和该拟声波方程组不同的一组约束来调整所述一个或多个地震波场。

根据一些实施例，在具有一个或多个处理器和存储器的电子设备处执行另一方法。该方法包括提供用于地质介质的地球模型。该地质介质具有非均匀倾斜的对称轴，并且该地球模型包括对称轴方向上的非零剪切速度。该设备使用根据地球模型应用的二阶拟声波方程组处理一个或多个地震测量。该拟声波方程组包括描述一个或多个地震波场的第一方程和描述一个或多个地震波场的第二方程。该处理包括根据拟声波方程组在多个时间步长上传播所述一个或多个地震波场。传播所述一个或多个地震波场包括计算地质介质内的多个位置处的所述一个或多个地震波场。针对地质介质内的多个位置的子集，在完成多个时间步长中的每个时间步长时，该设备应用与地球模型和拟声波方程组不同的一组约束以调整所述一个或多个地震波场。

根据一些实施例，在具有一个或多个处理器和存储器的电子设备处执行另一方法。该方法包括接收与多个源和接收器位置对应的一个或多个地震测量。该方法还包括提供用于地质介质的地球模型。该地质介质具有非均匀倾斜的对称轴并且针对地质介质内的至少一个位置子集该地球模型包括对称轴方向上的非零剪切速度。该方法还包括根据地球模型和能量守恒拟声波方程组在多个时间步长上传播所述一个或多个地震测量。该能量守恒拟声波方程组包括描述一个或多个地震波场的第一方程和描述一个或多个地震波场的第二方程。能量守恒拟声波方程组是通过对能量非守恒拟声波方程组中的一个或多个地震波场的一个或多个导数项进行近似从能量非守恒拟声波方程组导出的。

在本发明的另一方面中，为了解决上述问题，一些实施例提供了存储一个或多个程序的非临时性计算机可读存储介质。所述一个或多个程序包括指令，所述指令当被具有一个或多个处理器和存储器的电子设备执行时使电子设备执行本文中提供的方法中的任一个。

在本发明的又一方面中，为了解决上述问题，一些实施例提供了电子设备。该电子设备包括一个或多个处理器、存储器和一个或多个程序。所述一个或多个程序被存储在存储器中并且被配置为被一个或多个处理器执行。所述一个或多个程序包括操作系统和指令，其当被所述一个或多个处理器执行时使电子设备执行本文中提供的方法中的任一个。

附图说明

图1是根据一些实施例的海洋地震数据采集环境的示意图。

图2是根据一些实施例的说明地震建模系统的框图。

图3A-3D是根据一些实施例的减轻拟声波模型中的不稳定性的方法的示意性流程图。

图4是根据一些实施例的减轻拟声波模型中不稳定性的另一方法的示意性流程图。

图5说明了根据一些实施例计算的多个模拟波场快照。

图6说明了根据一些实施例计算的多个模拟波场快照。

图7说明了根据一些实施例计算的多个模拟叠前图像。

图8是根据一些实施例的减轻拟声波模型中的不稳定性的另一方法的示意性流程图。

相似的参考标号在整个附图中指对应的部分。

具体实施方式

下面详细描述的是对地质介质中的波场进行建模的方法，该方法减轻和/或消除与对完全弹性波动方程的拟声波近似相关联的不稳定性。根据一些实施例，提供了一种用于地质介质的地球模型(有时被称为“速度模型”)，其指定地质介质的各种特性，诸如剪切(S)波速度、压缩(P)波速度、刚度参数(例如，刚度张量)和/或各向异性参数(例如，倾斜角或Thomsen参数δ)、以及用于地质介质内的多个位置(例如，由计算栅格指定的位置)的密度。针对地质介质内的至少一个位置子集(例如，在地质上有利的地方)该地球模型包括对称轴方向上的非零剪切速度。基于该地球模型，使用至完全弹性波动方程的拟声波近似，在多个时间步长上传播(例如，时间步进)一个或多个地震波场。在一些实施例中，每个波场表示地质介质中的每个位置(例如，计算栅格上的每个位置)处的相应状态变量的值。为简单起见，术语“波场”和“状态变量”以及“地震波场”在本文中在一些实施例中被可交换地使用。

在一些实施例中，对于地质介质内的选定位置(例如，包括不支持剪切波的盐或水介质的位置)，一组约束(例如，与椭圆形各向异性一致的约束)被应用以“重置”地震波场并且从而防止地震波场的不稳定方面超出可接受水平。选择性地或此外，在一些实施例中，根据时间步长准则(例如，地质介质中的所有位置每1000个时间步长就被重置)，“重置”操作被应用到地质介质内的一些或所有位置。仍然可选择地或此外，拟声波近似包括通过对能量非守恒拟声波方程组中的地震波场的一个或多个导数项进行近似从能量非守恒拟声波方程(例如，耦合方程)组导出的能量守恒拟声波方程(例如，耦合方程)组。

现在将详细参照各种实施例，实施例的示例被在附图中说明。在下面的具体实施方式中，许多具体的细节被阐述以便提供对本公开内容以及本文中描述的实施例的全面理解。然而，在没有这些具体细节的情况下，本文中描述的实施例可以被实践。在其它情况下，众所周知的方法、过程、组件和机械装置没有被详细描述，以免不必要地模糊实施例的各方面。

图1是根据一些实施例的海洋地震数据采集环境100的示意图。在环境100中，勘测船102将一条或多条地震拖缆(图1中描绘的一个示例性拖缆104)拖曳在船102后面。地震拖缆104可以是几千米长并且可以包含可以用于支持沿着拖缆104的通信的各种支撑线缆(未示出)以及布线和/或电路(未示出)。通常，拖缆104包括记录地震信号的地震传感器106(例如，地震传感器106-a、106-b、106-c至地震传感器106-n)所安装到的主缆。

在一些实施例中，地震传感器106可以仅是压力传感器或可以是多分量地震传感器。对于多分量地震传感器的情况，每个传感器能够检测压力波场以及与邻近多分量地震传感器的声学信号相关联的粒子运动的至少一个分量。粒子运动的示例包括粒子位移的一个或多个分量(例如，主测线(inline)(x)、联络线(crossline)(y)和/或垂直(z)分量中的一个或多个(例如，如轴108中所示))、粒子速度的一个或多个分量、以及粒子加速度的一个或多个分量。

在一些实施例中，多分量地震传感器可以包括一个或多个水听器、地震检波器、粒子位移传感器、粒子速度传感器、加速计、压力梯度传感器、或它们的组合。

例如，在一些实施例中，特定的多分量地震传感器可以包括用于测量压力的地震检波器和测量地震传感器附近的粒子速度和/或加速度的三个对应正交分量的三个正交对准加速度计。注意的是多分量地震传感器可以被实现为单个设备或可以被实现为多个设备。特定的多分量地震传感器也可以包括一个或多个压力梯度传感器，其构成另一种类型的粒子运动传感器。每个压力梯度传感器测量相对于特定方向的特定点处的压力波场的改变。例如，压力梯度传感器中的一个可以获得指示特定点处压力波场相对于联络线方向的偏导数的地震数据，并且压力梯度传感器中的另一个可以获得特定点处指示相对于主测线方向的压力数据的地震数据。

海洋地震数据采集环境100包括一个或多个地震源阵列110。源阵列110又包括诸如空气枪的地震源(例如，地震源112-a、112-b、112-c至地震源112-n)的一个或多个串。在一些实施例中，地震源112可以被耦接到勘测船102或被勘测船102拖曳。选择性地，地震源112可以独立于勘测船102操作，因为源元件112可以被耦接到例如其它船或浮标。

当地震拖缆104被拖曳在勘测船102后面时，声学信号114(有时被称为“放炮(shots)”)由地震源112产生并且通过水柱116向下引入水底面120之下的岩层(strata)118(例如，岩层118-a、118-b、118-c、118-d和118-e，各表示地质介质的相应层)中。反射声学信号122被从各种地下地质特征(诸如示例性地层124(例如，盐丘))反射。

入射声学信号114产生对应的反射声学信号或压力波，其被地震传感器106感测。注意的是被地震传感器106接收和感测的压力波包括传播至传感器106而无反射的“上行”压力波以及由压力波从空气-水边界126反射而产生的“下行”压力波。

地震传感器126生成称为“道”的信号(例如，数字信号)，这种信号指示所获得的压力波场和粒子运动(如果传感器是粒子运动传感器)的测量。根据一些实施方式，道被记录并且可以至少部分地被部署在勘测船102上的信号处理单元128处理。

地震采集的目的是为了识别地下地质介质(例如示例性地质地层124(例如，盐丘))而构建勘测区域的图像。对表示的后续分析可以揭示地下地质介质中烃类沉积物的可能位置。在一些实施例中，对表示所进行的部分分析可以在地震勘测船102上执行，诸如由信号处理单元128执行。在一些实施例中，表示可以由例如可位于陆地上或者船102上的地震建模系统(诸如在图2中描绘的并在下面进一步描述的示例性地震建模系统200)处理。因此，许多变化是可能的，并且在所附权利要求的范围之内。

本领域技术人员将理解的是，以上描述的海洋地震数据采集环境100仅是可以使用的许多不同类型的地震数据采集环境中的一个的示例。例如，在一些实施例中，地震数据采集环境可以使用设置在海床上的固定传感器线缆。作为另一示例，在一些实施例中，地震数据采集环境可以是其中传感器线缆被埋在地球中的基于陆地的环境。因此，许多变化被构想，并且在所附权利要求的范围之内。

地震偏移是为了对地球中的反射性分布进行成像而典型使用的处理。在一些实施例中，地震偏移涉及在时间上向后传播(例如，时间步进)记录在地震传感器处的信号以计算在对应反射点处的一个或多个波场。在一些实施例中，地震偏移也涉及在时间上向前传播(时间步进)在地震源处产生的信号以计算对应反射点处的波场。地震偏移可以使地震各向异性(例如，具有陡峭或非均匀的对称轴方向的倾斜的横向各向异性)复杂化，地震各向异性可以在相应信号的向后和向前传播中引入不稳定性，从而破坏成像处理。下面更详细描述的是通过重置在某些传播时间的波场和/或重置地质介质内的位置处的波场来减轻或消除这些不稳定性的系统和方法(例如，方法400，图4和/或方法500，图5)。下面还详细描述的是通过使用能量守恒拟声波传播算子(例如，拟声波方程组)传播波场来减轻或消除这些不稳定性的系统和方法。

图2是根据一些实施例的说明地震建模系统200的框图。尽管某些具体特征被说明，但本领域技术人员将从本公开内容理解的是为了简明起见各种其它特征没有被说明，以免模糊本文中公开的实施例的更多相关方面。

为此，地震建模系统200包括一个或多个处理单元(CPU)202、一个或多个网络或其它通信接口208、存储器206和用于互连这些和各种其它组件的一个或多个通信总线204。地震建模系统200还可选地包括一个或多个地震传感器106(例如，地震检波器和/或水听器)、一个或多个地震源112(例如，空气枪)。通信总线204可以包括互连和控制系统组件之间的通信的电路(有时被称为芯片组)。存储器206包括高速随机存取存储器，诸如DRAM、SRAM、DDR RAM或其它随机存取固态存储器设备；并且可以包括非易失性存储器，诸如一个或多个磁盘存储设备、光盘存储设备、闪存设备或其它非易失性固态存储设备。存储器206可以可选地包括远离(一个或多个)CPU202的一个或多个存储设备。包括存储器206内的(一个或多个)非易失性和易失性存储器设备的存储器206包括非暂时性计算机可读存储介质。

在一些实施例中，存储器206或存储器206的非暂时性计算机可读存储介质存储下面的程序、模块和数据结构或包括操作系统216、网络通信模块218、地震建模模块220的它们的子集。

操作系统216包括用于处理各种基本系统服务和用于执行硬件相关任务的过程。

网络通信模块218经由通信网络接口208(有线或无线)和一个或多个通信网络(诸如互联网、其它广域网、局域网、城域网等)促进与其它设备的通信(例如，促进与地震源112和/或地震传感器106(如果不被包括在系统200中的话)的通信，或促进与其它基于陆地的组件的通信)(例如，在一些实施例中，地震建模系统200远离地震源112和/或地震传感器106)。

在一些实施例中，地震建模模块220被配置为接收用于地质介质的地球模型。地震建模模块220使用根据地球模型应用的二阶拟声波方程组处理(例如，从地震传感器106接收的)一个或多个地震测量以计算一个或多个地震波场。地震建模模块220根据拟声波方程组在多个时间步长上传播地震波场。根据(下面参照方法300，图3，和方法400，图4描述的)预定准则，地震建模模块220应用一组约束(例如，椭圆形约束或各向同性约束)来调整(或“重置”)地震波场。

为此，地震建模模块220可选地包括一个或多个子模块，每个子模块包括一组指令并且可选地包括元数据和参数。例如，在一些实施例中，地震建模模块220使用传播子模块224(其包括一组指令224-1和元数据和参数224-2，其中元数据和参数可以包括时间步长持续时间)传播地震波场，使用约束子模块222(其包括一组指令222-1和元数据和参数222-2)将该组约束应用到地震波场，并且将数据存储在数据子模块226(其包括地球模型226-1和地震数据226-2)中。

图3A-3D是根据一些实施例的减轻拟声波模型中的不稳定性的方法300的示意性流程图。该方法可选地由指令支配，所述指令存储在计算机存储器或非暂时性计算机可读存储介质(例如，图2中的存储器212)中并且由一个或多个计算机系统的一个或多个处理器(例如，(一个或多个)处理器202)执行，所述一个或多个计算机系统包括但不限于信号处理单元128(图1)和/或系统200(图2)。计算机可读存储介质可以包括磁盘或光盘存储设备、诸如闪存的固态存储设备、或一个或多个其它非易失性存储器设备。存储在计算机可读存储介质上的计算机可读指令可以包括以下中的一个或多个：源代码、汇编语言代码、目标码、或由一个或多个处理器解译的其它指令格式。在各种实施例中，每种方法中的一些操作可以被组合和/或一些操作的次序可以被从图中示出的次序改变。此外，在一些实施例中，单独图中示出的和/或与单独方法(例如，方法400，图4和/或方法800，图8A-8C)相关联讨论的操作可以被组合以形成其它方法，并且同一图中示出的和/或与同一方法相关联讨论的操作可以被分离成不同的方法。此外，在一些实施例中，方法中的一个或多个操作被图2中示出的系统200的模块(例如，(一个或多个)处理器202、地震建模模块220、存储器212、网络接口210和/或它们的任意子模块)执行。

地震建模系统提供(302)用于地质介质的地球模型。该地质介质具有非均匀倾斜的对称轴，并且针对地质介质内的至少一个位置子集该地球模型包括对称轴方向上的非零剪切速度。该地球模型可以包括压力波和剪切波速度、各向异性参数(例如，对称轴的方向，有时被称为“倾斜角”)和/或地质介质中的多个位置处(例如，计算栅格中的每个位置处)的刚度张量C的一个或多个值的任意组合。在一些实施例中，C是与Thomsen参数δ和∈相关的对角6×6刚度张量，其中：

$C=\left[\begin{array}{cccccc}{C}_{11}& C_{11}-2C_{66}& C_{13}& 0& 0& 0\\ C_{11}-2C_{66}& C_{11}& C_{13}& 0& 0& 0\\ C_{13}& C_{13}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right]---\left(1\right)$

其中：

$C_{33}={\mathrm{\ρV}}_{pz}^2---\left(2\right)$

C₁₁＝(1+2∈)C₃₃ (3)

$C_{13}=(\sqrt{f\left(f+2\mathrm{\δ}\right)}-(1-f\left)\right)C_{33}---\left(4\right)$

C₄₄＝(1-f)C₃₃ (5)

包括C₆₆的刚度张量C的所有其他元素等于零。并且，f是由以下等式给出的所谓的“f因子”：

$f=1-V_{SZ}^2/V_{PZ}^2---\left(6\right)$

f因子与垂直方向上的剪切波速度V_SZ和垂直方向上的压力波速度V_PZ有关。如本文中使用的，垂直和水平是相对于各向异性定义的：水平方向意味着任何垂直于倾斜的横向各向同性(TTI)介质的对称轴的方向并且垂直方向意味着沿着TTI介质的对称轴的方向(例如，垂直方向和水平方向被例如在计算栅格上的每个位置处定义)。在一些情况中，可以通过在空间中的每个点处(例如，依赖于对称轴的非均匀方向在空间中的每个点处不同地)旋转垂直横向的各向同性的刚度张量来概念地获取TTI介质。

地震建模系统使用根据地球模型应用的二阶拟声波方程组处理(304)一个或多个地震测量。在各种实施例中，所述一个或多个地震测量是现场记录的(例如，通过地震传感器106测量的，图1)、部分处理的和/或完全合成的(例如，对应于地震源112的模拟源脉冲(图1)，或不考虑任何特定源和/或传感器模拟的选择性模拟源和/或传感器测量)、或它们的组合。该拟声波方程组包括描述一个或多个地震波场的第一等式以及描述地震波场的第二等式。

作为示例，下面的2×2二阶拟声波微分方程(例如，耦合方程，有时被称为“波传播算子”)系统提供至计算上更昂贵的完全弹性TTI方程的近似：

${\∂}_t^{2}P=C_{11}{\mathrm{\Δ}}_{h}P+(C_{13}+C_{44}){\mathrm{\Δ}}_{v}Q+C_{44}{\mathrm{\Δ}}_{v}P---\left(7\right)$

${\∂}_t^{2}Q=(C_{13}+C_{44}){\mathrm{\Δ}}_{h}P+C_{33}{\mathrm{\Δ}}_{v}Q+C_{44}{\mathrm{\Δ}}_{h}Q---\left(8\right)$

其中，P和Q是描述波场的状态变量。无论如何，△_v是由介质的密度ρ归一化的垂直方向上的拉普拉斯算子(因为△_v是相对于单方向定义的，所以它也被称为和等于相对于垂直方向的“第二导数”)；并且△_h是由介质的密度ρ归一化的水平方向上的拉普拉斯算子(即，二维拉普拉斯算子)。更具体地，△_v和△_h由下面的等式给出，该等式也定义了相应的拉普拉斯算子被介质的密度ρ归一化的方式：

${\mathrm{\Δ}}_{h}u={\widetilde{\∂}}_{\hat{x}}\left(\frac{1}{\mathrm{\ρ}}{\∂}_{\hat{x}}u\right)+{\widetilde{\∂}}_{\hat{y}}\left(\frac{1}{\mathrm{\ρ}}{\∂}_{\hat{y}}u\right)---\left(9\right)$

${\mathrm{\Δ}}_{v}u={\widetilde{\∂}}_{\hat{z}}\left(\frac{1}{\mathrm{\ρ}}{\∂}_{\hat{z}}u\right)---\left(10\right)$

其中，是TTI介质的对称轴的笛卡尔(Cartesian)方向，并且和是垂直于TTI介质的对称轴的正交笛卡尔方向(u是用于任何适当的操作数(operand)的占位符)。微分算子和被在下面定义。拉普拉斯算子Δ_v和Δ_h通过下面的等式相关到“三维”拉普拉斯算子Δ_3D：

$\begin{array}{c}{\mathrm{\Δ}}_{3D}u={\mathrm{\Δ}}_{v}u+{\mathrm{\Δ}}_{h}u={\widetilde{\∂}}_{\hat{x}}\left(\frac{1}{\mathrm{\ρ}}{\∂}_{\hat{x}}u\right)+{\widetilde{\∂}}_{\hat{y}}\left(\frac{1}{\mathrm{\ρ}}{\∂}_{\hat{y}}u\right)+{\widetilde{\∂}}_{\hat{z}}\left(\frac{1}{\mathrm{\ρ}}{\∂}_{\hat{z}}u\right)\\ ={\∂}_x\left(\frac{1}{\mathrm{\ρ}}{\∂}_{x}u\right)+{\∂}_y\left(\frac{1}{\mathrm{\ρ}}{\∂}_{y}u\right)+{\∂}_z\left(\frac{1}{\mathrm{\ρ}}{\∂}_{z}u\right)\end{array}---\left(11\right)$

在以上的等式中，(x，y，z)指的是“非旋转”坐标(例如，由轴108指定的，图1)并且和是沿着非旋转坐标的“正常”偏导数。此外，在以上的等式中，微分算子被如下定义(例如，下面的算子实现了相对于非旋转坐标的旋转)：

${\∂}_{\hat{x}}u=cos\mathrm{\θ}cos\mathrm{\φ}{\∂}_{x}u+cos\mathrm{\θ}\sin \mathrm{\φ}{\∂}_{y}u-sin\mathrm{\θ}{\∂}_{z}u---\left(12\right)$

${\∂}_{\hat{y}}u=-sin\mathrm{\φ}{\∂}_{x}u+cos\mathrm{\φ}{\∂}_{y}u---\left(13\right)$

${\∂}_{\hat{z}}u=sin\mathrm{\θ}cos\mathrm{\φ}{\∂}_{x}u+sin\mathrm{\θ}sin\mathrm{\φ}{\∂}_{y}u+cos\mathrm{\θ}{\∂}_{z}u---\left(14\right)$

${\widetilde{\∂}}_{\hat{x}}u={\∂}_x\left(cos\mathrm{\θ}cos\mathrm{\φ}u\right)+{\∂}_y\left(cos\mathrm{\θ}sin\mathrm{\φ}u\right)-{\∂}_z\left(sin\mathrm{\θ}u\right)---\left(15\right)$

${\widetilde{\∂}}_{\hat{y}}u=-{\∂}_x\left(sin\mathrm{\φ}u\right)+{\∂}_y\left(cos\mathrm{\φ}u\right)---\left(16\right)$

${\widetilde{\∂}}_{\hat{z}}u={\∂}_x\left(sin\mathrm{\θ}cos\mathrm{\φ}u\right)+\∂y\left(sin\mathrm{\θ}sin\mathrm{\φ}u\right)+{\∂}_z\left(cos\mathrm{\θ}u\right)---\left(17\right)$

其中□是方位角并且θ是倾斜角。

应该理解的是，由于在如上等式中指示的适当位置处包含由1/ρ给出的浮力项，所以算子Δ_3D、Δ_v和Δ_h是广义拉普拉斯算子。为了便于解释，这些算子被简单称为“拉普拉斯算子”。此外，这些实施例不是旨在限制随后的权利要求；当在权利要求中使用时，术语“拉普拉斯”应当被解释为包括标准拉普拉斯算子以及诸如以上的更一般的形式，除非权利要求另外明确指示。

无论如何，继续方法300的描述，处理测量包括根据拟声波方程组在多个时间步长上传播(306)地震波场。例如，在一些实施例中，所述一个或多个测量被以前向时间步进方向或逆向时间步进方向传播。选择性地，在一些实施例中，一个或多个测量的第一子集被以前向时间步进方向传播并且一个或多个测量的第二子集被以逆向时间步进方向传播(例如，第一子集包括合成源“测量”并且第二子集包括测量的接收器测量)。

在一些实施例中，传播地震波场包括计算(308)介质内的多个位置处的地震波场(例如，由P和Q描述的地震波场)。在一些实施例中，在根据稳定性条件选择时间步长的同时，多个位置(例如，有时被称为“计算栅格”)被选择以避免空间波散。

例如，在一些实施例中，传播是使用应用到等式(7)和(8)的高阶有限差分算法(例如，时间上二阶并且空间上八阶的有限差分算法)实现的。在一些实施例中，在地质介质的边缘(例如，模型的边缘)处的边界条件被包括在有限差分算法中以完全指定该模型。例如，在一些情况(例如，诸如其中地质介质的顶部是水的海上勘探)下，边界条件包括地质介质顶部上的自由表面边界条件(例如，水表面)和模型的相应侧和底部上的吸收边界条件。

在一些实施例中，脉冲被注入等式(7)和(8)中并且随后被以前向时间方向(例如，当脉冲对应于已知的或估计的源信号时)或逆向时间方向(例如，当脉冲对应于测量的传感器信号时)在时间上(例如，时间步进地)传播。如本文中使用的，术语脉冲意味着用于拟声波方程组的适当的强制函数。例如，在一些实施例中，状态变量P对应于压力波场，并且此外可以被分解成响应函数(其中并且t是时间)以及强制函数使得在一些这种情况中，强制函数可以被写为：

$f(\stackrel{\‾}{r},t)=\mathrm{\γ}\mathrm{\δ}(\stackrel{\‾}{r}-{\stackrel{\‾}{r}}_0)f_s\left(t\right)---\left(18\right)$

其中γ是比例因子，是相应地震源或相应地震传感器的位置，δ是狄拉克△函数，以及f_s(t)是与相应的地震源或相应的地震传感器对应的源或传感器信号。

在一些情况中，前向建模包括当脉冲被以前向时间方向或逆向时间方向传播达当脉冲被地质特征反射时的相应时间时的偏移。在一些实施例中，在已知为逆时偏移(RTM)的过程中，源信号被以前向时间方向传播达相应时间并且传感器信号被以逆向时间方向传播达相应时间。前向时间传播的结果与逆向时间传播的结果之间的相关性(例如，结果的交叉相关性或结果的卷积)被用于产生地质介质的图像。

处理测量还包括确定(314)多个时间步长中的相应时间步长是否满足预定准则。在一些实施例中，确定多个时间步长中的相应时间步长是否满足预定准则包括确定(316)以时间步长的第一预定义数目为模的相应时间步长的数目是否等于预定义值。例如，在一些实施例中，每第1000个时间步长满足预定准则。这可以被数学表示为：

A＝t_YES(mod1,000) (19)

其中A是预定义值(例如，0、1、2、3等)，并且t_YES指示时间步长满足预定义准则。如下描述的，这容许地震波场每1000个时间步长就被重置，其限制了不稳定性的增长。然而，应当理解，第一预定义数目可以采取任意值(例如，500、2000等)。

在一些实施例中，确定多个时间步长中的相应时间步长是否满足预定义准则还包括确定(318)相应时间步长的数目是否大于时间步长的第二预定义数目(例如，t_YES>B，其中B是第二预定义数目)。

在一些实施例中，确定多个时间步长中的相应时间步长是否满足预定义准则包括(320)：计算(322)与相应时间步长相关联的误差度量并且比较(324)与相应时间步长相关联的该误差度量与预定义误差阈值。例如，地震建模系统可以计算与不稳定性增长对应的误差度量并且比较该误差度量与预定义误差阈值以便确定一个或多个地震波场是否正在变得不稳定而超出可接受水平。

应当理解的是(例如，相对于操作316-324)以上描述的准则可以被单独使用或可以(例如，使用逻辑算子(诸如“和”、“或”等))逻辑地组合。

根据相应时间步长满足预定义准则的确定，地震建模系统应用(326)与地球模型和拟声波方程组不同的一组约束来调整地震波场。该操作有时被称作“重置”操作，因为它限制了不稳定的增长。

例如，在一些情况中，地球模型包括(328)倾斜的横向各向同性模型并且应用该组约束包括应用与椭圆形各向异性模型一致的约束。在一些实施例中，地震波场包括(330)第一地震波场和第二地震波场，并且应用该组约束包括强制第一地震波场与第二地震波场之间的比例：

Q＝αP (20)

其中α是与椭圆形各向异性相关的比例常数。在一些实施例中，α通过以下关系式相关至Thomsen参数∈：

$\mathrm{\α}=\frac{1}{\sqrt{1+2\∈}}---\left(21\right)$

在一些实施例中，α等于1(unity)。这种情况有时被称作“各向同性约束”。

在一些实施例中，存在(336)强制第一地震波场与第二地震波场之间的比例的至少两种方式，其可以被单独使用或组合使用。例如，可以在多个时间步长的子集强制(337)跨整个地球模型的第一地震波场与第二地震波场之间的比例，并且在子集中的任意两个相邻时间步长之间存在多时间步长的预定义时间间隙。但是在可以被看作椭圆形各向异性的地球模型的子区内，可以在多个时间步长内的连续时间步长的子集强制(338)第一地震波场与第二地震波场之间的比例。

在一些实施例中，处理一个或多个地震测量包括对调整的一个或多个地震波场执行(331)逆时偏移以生成地质介质的地震图像。

图4是根据一些实施例的减轻拟声波模型中的不稳定性的方法400的示意性流程图。该方法可选地由指令支配，所述指令被存储在计算机存储器或非暂时性计算机可读存储介质(例如，图2中的存储器212)中并且由一个或多个计算机系统的一个或多个处理器(例如，(一个或多个)处理器202)执行，所述一个或多个计算机系统包括但不限于信号处理单元128(图1)和/或系统200(图2)。计算机可读存储介质可以包括磁盘或光盘存储设备、诸如闪存的固态存储设备或一个或多个其它非易失性存储器设备。存储在计算机可读存储介质上的计算机可读指令可以包括以下中的一个或多个：源代码、汇编语言代码、目标码、或由一个或多个处理器解译的其它指令格式。在各种实施例中，每种方法中的一些操作可以被组合和/或一些操作的次序可以被从图中示出的次序改变。此外，在一些实施例中，单独图中示出的和/或与单独方法(例如，方法300，图3和/或方法800，图8A-8C)相关联讨论的操作可以被组合以形成其它方法，并且同一图中示出的和/或与同一方法相关联讨论的操作可以被分离成不同的方法。此外，在一些实施例中，方法中的一个或多个操作被图2中示出的系统200的模块(包括，例如，(一个或多个)处理器202、地震建模模块220、存储器212、网络接口210和/或它们的任意子模块)执行。

地震建模系统提供(402)用于地质介质的地球模型。该地质介质具有非均匀倾斜的对称轴，并且针对地质介质内的至少一个位置子集该地球模型包括对称轴方向上的非零剪切速度。该地球模型可以包括压力波和剪切波速度、各向异性参数(例如，对称轴的方向)和/或以上参照等式(1)描述的刚度张量C的一个或多个值的任意组合。

地震建模系统使用根据地球模型应用的二阶拟声波方程组处理(404)一个或多个地震测量。在各种实施例中，所述一个或多个地震测量是现场记录的(例如，通过地震传感器106测量的，图1)、部分处理的和/或完全合成的(例如，对应于地震源112的模拟源脉冲(图1)，或不考虑任何特定源和/或传感器模拟的选择性模拟源和/或传感器测量)、或它们的组合。该拟声波方程组包括描述一个或多个地震波场的第一等式(例如，等式(7))以及描述地震波场的第二等式(例如，等式(8))。

处理地震测量包括根据拟声波方程组在多个时间步长上传播(406)地震波场。传播地震波场包括计算地质介质内的多个位置处的地震波场。例如，在一些实施例中，一个或多个测量被以前向时间步进方向或逆向时间步进方向传播。选择性地，在一些实施例中，一个或多个测量的第一子集被以前向时间步进方向传播并且一个或多个测量的第二子集被以逆向时间步进方向传播(例如，第一子集包括合成源“测量”并且第二子集包括测量的接收器测量)。

对于地质介质内的多个位置的子集，在完成多个时间步长中的每个时间步长时，地震建模系统应用(408)与地球模型和拟声波方程组不同的一组约束来调整地震波场(例如，如上参照方法300描述的约束中的任一约束)。在一些实施例中，地质介质内的多个位置的子集包括地质介质中的位置，该地质介质包括水和/或盐中的至少一个(例如，计算栅格中的多个位置的子集是不支持剪切波和/或倾斜的横向各向同性的地质介质的区)。在一些实施例中，在多个位置的子集处，地震波场在每个时间步长后就被“重置”。

图5说明了根据一些实施例计算的多个模拟波场快照。特别地，面板500-1和面板500-2均是通过将子波注入等式(7)和(8)并且在注入子波后将对应的波场传播达时间T＝4秒(例如，子波的注入模拟来自地震源的脉冲)而计算的波场快照。面板500-1与面板500-2之间的差异在于面板500-2中的波场是通过每1000个时间步长(在该示例中，时间步长等于0.7毫秒)就应用椭圆形约束(例如，应用等式(14)，被称作“重置”操作)而计算的，而没有重置操作用于面板500-1中的波场。面板500-3说明了面板500-1与面板500-2中的波场之间的计算差(为了清楚起见，在面板500-3中差乘以因子10)。

图6说明了根据一些实施例计算的多个模拟波场快照。特别地，面板600-1和面板600-2说明了与面板500-1和面板500-2(图5)相同的情形，然而，在面板600-1和600-2中，相应波场在注入子波之后已经被时间步进达时间T＝12秒(与面板500-1和500-2中的T＝4秒相比)。面板600-3说明了面板600-1和面板600-2中的波场之间的计算差(为了清楚起见，在面板600-3中差乘以因子10)。

图7说明了根据一些实施例计算的多个模拟叠前图像(示出具有用于球面散度的校正)。特别地，面板700-1说明了其中每个垂直列表示将被定位在对应水平位置处的地震传感器接收的道的叠前图像。面板700-2也说明了其中每个垂直列表示将被定位在对应水平位置处的地震传感器接收的道的叠前图像。面板700-1与面板700-2之间的差异在于面板700-2中的叠前图像是通过每1000个时间步长(在该示例中，时间步长等于0.7毫秒)就应用椭圆形约束(例如，应用等式(13)，被称作“重置”操作)而计算的，而没有重置操作用于面板700-1中的叠前图像。面板700-3说明了面板700-1和面板700-2中的叠前图像之间的计算差(为了清楚起见，在面板700-3中差乘以因子10)。从图像中清楚的是在没有重置操作的情况下不稳定性随时间增长，并且该不稳定性致使在大约8秒的时间后面板700-1中的叠前图像无效。通过比较，面板700-2中的图像在8秒后很长时间仍保持清楚，其又蕴含着地质介质的更大深度处的有用信息。

图8A-8C是根据一些实施例的减轻拟声波模型中的不稳定性的方法300的示意性流程图。该方法可选地由指令支配，所述指令被存储在计算机存储器或非暂时性计算机可读存储介质(例如，图2中的存储器212)中并且由一个或多个计算机系统的一个或多个处理器(例如，(一个或多个)处理器202)执行，所述一个或多个计算机系统包括但不限于信号处理单元128(图1)和/或系统200(图2)。计算机可读存储介质可以包括磁盘或光盘存储设备、诸如闪存的固态存储设备或一个或多个其它非易失性存储器设备。存储在计算机可读存储介质上的计算机可读指令可以包括以下中的一个或多个：源代码、汇编语言代码、目标码、或由一个或多个处理器解译的其它指令格式。在各种实施例中，每种方法中的一些操作可以被组合和/或一些操作的次序可以被从图中示出的次序改变。此外，在一些实施例中，单独图中示出的和/或与单独方法(例如，方法300，图3和/或方法400，图4)相关联讨论的操作可以被组合以形成其它方法，并且同一图中示出的和/或与同一方法相关联讨论的操作可以被分离成不同的方法。此外，在一些实施例中，方法中的一个或多个操作被图2中示出的系统200的模块(包括，例如，(一个或多个)处理器202、地震建模模块220、存储器212、网络接口210和/或它们的任意子模块)执行。

地震建模系统接收(802)与多个源和接收器位置(例如，地震传感器106和/或地震源112，图1)对应的一个或多个地震测量。在一些实施例中，地震建模系统接收与多个接收器位置对应的一个或多个地震测量并且也接收与多个源位置对应的一个或多个模拟源脉冲。模拟源脉冲是真实脉冲(例如，从地震源112释放的脉冲)的模拟表示，真实脉冲传播通过地质介质(例如，包括水)并且其对应的波场被地震传感器106接收和测量，从而导致地震测量。为了便于解释，术语“地震测量”应该理解为包括来自地震传感器的测量和来自地震源的模拟脉冲两者。这是因为在一些环境中，这两者是以类似的方式被操作的(例如，被处理和/或传播通过地质介质)。在各种实施例中，所述一个或多个地震测量是现场记录的(例如，通过地震传感器106测量的，图1)、部分处理的和/或完全合成的(例如，对应于地震源112的模拟源脉冲(图1)，或不考虑任何特定源和/或传感器模拟的选择性模拟源和/或传感器测量)、或它们的组合。

地震建模系统提供(804)用于地质介质(例如，包括岩层118和/或示例性地层124的地质介质，图1)的地球模型。地质介质具有非均匀倾斜的对称轴，并且该地球模型包括对称轴方向上的非零剪切速度。在一些实施例中，地球模型被存储在地震建模系统的非暂时性计算机可读存储介质(例如，图2中的存储器212)中并且被地震建模系统提供给一个或多个处理器(例如，(一个或多个)处理器202，图2)用于处理。在一些实施例中，方法800被迭代地执行并且地球模型被(例如，根据误差度量)在迭代之间更新。在一些实施例中，地震建模系统接收来自远程系统(例如，服务器)的地球模型，远程系统可以是基于陆地的(例如，在地震建模系统在航海船上的情况下)。在一些实施例中，地球模型包括(806)倾斜的横向各向同性模型。该地球模型可以包括压力波和剪切波速度、各向异性参数(例如，对称轴的方向)和/或以上参照等式(1)描述的刚度张量C的一个或多个值的任意组合。

地震建模系统根据地球模型和能量守恒拟声波方程组在多个时间步长上传播(808)地震测量。例如，在一些实施例中，所述一个或多个测量被以前向时间步进方向或逆向时间步进方向传播。选择性地，在一些实施例中，一个或多个测量的第一子集被以前向时间步进方向传播并且一个或多个测量的第二子集被以逆向时间步进方向传播(例如，第一子集包括合成源“测量”并且第二子集包括测量的接收器测量)。

能量守恒拟声波方程组包括描述一个或多个地震波场的第一等式(例如，等式(37))和描述地震波场的第二等式(例如，等式(38))。此外，能量守恒拟声波方程组是通过对能量非守恒拟声波方程组中的地震波场的一个或多个导数项进行近似从能量非守恒拟声波方程组中导出的(810)。

考虑下面的非限制性示例，其中描述一个或多个地震波场的能量守恒拟声波方程组是从能量非守恒拟声波方程组导出的。等式(7)和(8)构成能量非守恒拟声波方程组并且可以以下面的形式使用Thomsen参数和f值重写：

${\∂}_t^{2}P=C_{33}\{(f+2\∈){\mathrm{\Δ}}_{h}P+\sqrt{f(f+2\mathrm{\δ})}{\mathrm{\Δ}}_{v}Q+(1-f){\mathrm{\Δ}}_{3D}P\}---\left(22\right)$

${\∂}_t^{2}Q=C_{33}\{\sqrt{f(f+2\mathrm{\δ})}{\mathrm{\Δ}}_{h}P+{\mathrm{f\Δ}}_{v}Q+(1-f){\mathrm{\Δ}}_{3D}Q\}---\left(23\right)$

在一些实施例中，地震波场包括(812)第一地震波场和第二地震波场。第一地震波场和第二地震波场中的一个或多个包括压力波场和剪切波场的线性组合。包括波场P和波场Q的线性组合的新波场M被引入以进一步简化等式(22)和(23)：

$M=\frac{(f+2\∈)Q-\sqrt{f(f+2\mathrm{\δ})}P}{\sqrt{2f(\∈-\mathrm{\δ})}}---\left(24\right)$

相对于时间对等式(24)微分两次得到：

${\∂}_t^{2}M=\frac{(f+2\∈){\∂}_t^{2}Q-\sqrt{f(f+2\mathrm{\δ})}{\∂}_t^{2}P}{\sqrt{2f(\∈-\mathrm{\δ})}}---\left(25\right)$

等式(24)可以(即，通过代数地求解Q)被等价地写为：

$Q=\frac{\sqrt{f(f+2\mathrm{\δ})}P+\sqrt{2f(\∈-\mathrm{\δ})}M}{(f+2\∈)}---\left(26\right)$

将拉普拉斯导数Δ(例如，由Δ_3D、Δ_h、和/或Δ_v组成的组中的任一)应用到等式(26)得到：

$\mathrm{\Δ}Q=\mathrm{\Δ}\left(\frac{\sqrt{f(f+2\mathrm{\δ})}P+\sqrt{2f(\∈-\mathrm{\δ})}M}{(f+2\∈)}\right)---\left(27\right)$

如上所述，能量守恒拟声波方程组是通过对地震波场的一个或多个导数项进行近似导出的。在一些实施例中，近似的导数项包括(814)一个或多个拉普拉斯导数项。在一些实施例中，近似的导数项沿着相应方向包括(816)至少一个导数项，该导数项通过相对于相应方向将各向异性参数视为常数而被近似。例如，下面的拉普拉斯导数项可以通过相对于任一或所有方向(例如，1、2或3个不同的笛卡尔空间方向)将各向异性参数∈和δ视为常数而被近似：

$\mathrm{\Δ}\left(\frac{\sqrt{f(f+2\mathrm{\δ})}}{(f+2\∈)}P\right)\≈\frac{\sqrt{f(f+2\mathrm{\δ})}}{(f+2\∈)}\mathrm{\Δ}P---\left(28\right)$

$\mathrm{\Δ}\left(\frac{\sqrt{2f(\∈-\mathrm{\δ})}}{(f+2\∈)}M\right)\≈\frac{\sqrt{2f(\∈-\mathrm{\δ})}}{(f+2\∈)}\mathrm{\Δ}M---\left(29\right)$

因此，等式(27)通过如下等式被近似：

$\mathrm{\Δ}Q\≈\frac{\sqrt{f(f+2\mathrm{\δ})}\mathrm{\Δ}P+\sqrt{2f(\∈-\mathrm{\δ})}\mathrm{\Δ}M}{(f+2\∈)}---\left(30\right)$

使用以上近似，由等式(22)和(23)组成的方程组可以被以以下近似的形式表示，尽管仍然作为能量非守恒拟声波方程组：

${\∂}_t^{2}P=C_{33}\[(f+2\∈){\mathrm{\Δ}}_{h}P+\sqrt{f(f+2\mathrm{\δ})}{\mathrm{\Δ}}_{v}Q+(1-f){\mathrm{\Δ}}_{3D}P\]---\left(31\right)$

${\∂}_t^{2}M=C_{33}\[\sqrt{2f(\∈-\mathrm{\δ})}{\mathrm{\Δ}}_{v}Q+(1-f){\mathrm{\Δ}}_{3D}M\]---\left(32\right)$

继续近似至地震波场的一个或多个导数项，在一些实施例中，近似的导数项中的相应一个是(818)地震波场中的相应一个的拉普拉斯导数并且近似的导数项中的相应一个是通过利用包括涉及一个或多个各向异性参数的拉普拉斯导数(例如，拉普拉斯导数对包括一个或多个各向异性参数的操作数进行运算)的项对地震波场中的相应一个的拉普拉斯导数进行近似而被近似的。例如，下面的用于u(例如，其中u是用于P或M的占位符)的三维拉普拉斯导数项可以被如下近似：

$(1-f){\mathrm{\Δ}}_{3D}u\≈(f+2\∈)\sqrt{\frac{1-f}{f+2\∈}}{\mathrm{\Δ}}_{3D}\left(\sqrt{\frac{1-f}{f+2\∈}}u\right)---\left(33\right)$

这导致下面的能量守恒拟声波方程组：

$\begin{array}{c}{\∂}_t^{2}P=C_{33}\[(f+2\∈){\mathrm{\Δ}}_{h}P+\sqrt{f(f+2\mathrm{\δ})}{\mathrm{\Δ}}_{v}Q+(f\\ +2\∈)\sqrt{\frac{1-f}{f+2\∈}}{\mathrm{\Δ}}_{3D}\left(\sqrt{\frac{1-f}{f+2\∈}}P\right)\]\end{array}---\left(34\right)$

${\∂}_t^{2}M=C_{33}\[\sqrt{2f(\∈-\mathrm{\δ})}{\mathrm{\Δ}}_{v}Q+(f+2\∈)\sqrt{\frac{1-f}{f+2\∈}}{\mathrm{\Δ}}_{3D}\left(\sqrt{\frac{1-f}{f+2\∈}}M\right)\]---\left(35\right)$

$Q=\frac{\sqrt{f(f+2\mathrm{\δ})}P+\sqrt{2f(\∈-\mathrm{\δ})}M}{(f+2\∈)}---\left(36\right)$

如以上可以看到的，在一些实施例中，传播地震测量包括(820)计算四个二阶导数(例如，恰好四个)，包括：

·相对于各向异性平面的第一地震波场的二维拉普拉斯导数，各向异性平面局部垂直于对称轴(例如，Δ_hP，等式(34))；

·对应于第一地震波场的项的三维拉普拉斯导数(例如，)，等式(34))；

·对应于第二地震波场的项的三维拉普拉斯导数(例如，)，等式(35))；以及

·相对于对称轴方向的第三地震波场的二阶导数(例如，Δ_vQ，等式(35))。

第三地震波场是第一地震波场和第二地震波场的线性组合(例如，如等式(36)指示的)。

应当注意的是，根据一些实施例等价的能量守恒拟声波方程组(例如，等价于等式(37)和(38))可以被如下表示：

$\begin{array}{c}{\∂}_t^{2}P=C_{33}\[(f+2\∈){\mathrm{\Δ}}_{h}P\\ +\sqrt{f(f+2\mathrm{\δ})}{\mathrm{\Δ}}_v\left(\frac{\sqrt{f(f+2\mathrm{\δ})}P+\sqrt{2f(\∈-\mathrm{\δ})}M}{(f+2\∈)}\right)\\ +(f+2\∈)\sqrt{\frac{1-f}{f+2\∈}}{\mathrm{\Δ}}_{3D}\left(\sqrt{\frac{1-f}{f+2\∈}}P\right)\]\end{array}---\left(37\right)$

$\begin{array}{c}{\∂}_t^{2}M=C_{33}\[\sqrt{2f(\∈-\mathrm{\δ})}{\mathrm{\Δ}}_v\left(\frac{\sqrt{f(f+2\mathrm{\δ})}P+\sqrt{2f(\∈-\mathrm{\δ})}M}{(f+2\∈)}\right)+(f\\ +2\∈)\sqrt{\frac{1-f}{f+2\∈}}{\mathrm{\Δ}}_{3D}\left(\sqrt{\frac{1-f}{f+2\∈}}M\right)\]\end{array}---\left(38\right)$

在等式(37)和(38)的情况中，传播地震测量包括计算五个二阶导数(例如，恰好五个)，包括：

·相对于各向异性平面的第一地震波场的二维拉普拉斯导数，各向异性平面局部垂直于对称轴(例如，Δ_hP，等式(37))；

·对应于第一地震波场的项的三维拉普拉斯导数(例如，)，等式(37))；

·对应于第二地震波场的项的三维拉普拉斯导数(例如，)，等式(38))；以及

·相对于对称轴方向的对应于第一地震波场的项的二阶导数(例如，)，等式(37)和(38))。

·相对于对称轴方向的对应于第二地震波场的项的二阶导数(例如，)，等式(37)和(38))。

在一些实施例中，传播地震波场包括(822)在多个时间步长中的每个时间步长处计算地质介质内的多个位置处的地震波场(例如，使用等式(34)-(36)计算三个地震波场和/或使用等式(37)和(38)计算两个地震波场)。在一些实施例中，传播地震测量还包括(824)对地震波场执行逆时偏移。

备选的能量守恒拟声波方程组也可以被导出为：

$\begin{array}{c}{\∂}_t^{2}P=C_{33}\[(f+2\mathrm{\ϵ}){\mathrm{\Δ}}_{h}P\\ +\sqrt{f(f+2\mathrm{\δ})}{\mathrm{\Δ}}_{v}Q+(f\\ +2\mathrm{\ϵ})\underset{j=1}{\overset{3}{\Σ}}{\∂}_j\left(\frac{1}{\mathrm{\ρ}}\frac{1-f}{f+2\mathrm{\ϵ}}{\∂}_j\right(P\left)\right)\]\end{array}---\left(39\right)$

$\begin{array}{c}{\∂}_t^{2}M=C_{33}\[\sqrt{2f(\mathrm{\ϵ}-\mathrm{\δ})}{\mathrm{\Δ}}_{v}Q+(f\\ +2\mathrm{\ϵ})\underset{j=1}{\overset{3}{\Σ}}{\∂}_j\left(\frac{1}{\mathrm{\ρ}}\frac{1-f}{f+2\mathrm{\ϵ}}{\∂}_j\right(M\left)\right)\]\end{array}---\left(40\right)$

等式(39)和(40)的系统也完全能量守恒并且因此是稳定的。也可以示出的是该等式系统在一些情况中在保持原始系统幅度方面对于地球模型参数变化是更稳健(robust)的。该等式系统与等式(37)和(38)的不同在于其中作用在缩放的波场上的拉普拉斯算子通过在两个一阶导数之间应用项(诸如)而被取代，从而导致针对更宽范围的各向异性参数的不同的更稳健的幅度成缩放的性质。

在一些情况中，地震波场中的每个在多个时间步长中的每个处通过存储在其中的相应能量被表征(826)。对于多个时间步长中的每个时间步长，存储在地震波场中的相应能量的和是常数。以这种方式，拟声波方程组是能量守恒拟声波方程组。

尽管如上描述了特定实施例，但是应该明白的是其不意在将本发明限制于这些特定实施例。相反地，本发明包括在随附权利要求的精神和范围内的替代、修改和等价物。许多具体细节被阐述以便提供对本文中呈现的主题的全面理解。但是对本领域技术人员将明显的是在没有这些具体细节的情况下该主题也可以被实践。在其它情况中，已知的方法、过程、组件和电路没有被详细描述，以免不必要地模糊实施例的各方面。

本文中的本发明的说明书中使用的术语仅用于描述特定实施例的目的并且不意在限制本发明。如本发明的说明书和随附权利要求中所使用的，单数形式“一”、“一个”和“该”意在也包括复数形式，除非上下文另外明确指示。也将理解的是，如本文中使用的术语“和/或”指的是并且包含相关联列出的项中的一个或多个的任意和所有可能的组合。还将理解的是，术语“包括”、“包含”、“含有”和/或“包括有”当在说明书中使用时指定陈述的特征、操作、元件、和/或部件的存在，但不排除一个或多个其它特征、操作、元件、部件、和/或它们的组的存在或增加。

如本文中使用的，术语“如果”可以被解释为意味着“当……时”或“在……时”或“响应于确定……”或“根据……的确定”或“响应于检测到……”陈述的先决条件为真，这取决于上下文。类似地，短语“如果确定[陈述的先决条件为真]”或“如果[陈述的先决条件为真]”或“当[陈述的先决条件为真]时”可以被解释为意味着“在确定……时”或“响应于确定……”或“根据……的确定”或“在检测到……时”或“响应于检测……”陈述的先决条件为真，这取决于上下文。

尽管各种附图中的一些以特定次序说明了多个逻辑阶段，但不依赖于次序的阶段可以被重新排序并且其它阶段可以被组合或打破。虽然一些重新排序或其它分组被特别提到，但是对于本领域技术人员而言其它的将是明显的并且因此没有呈现替代的详尽清单。此外，应该认识到，阶段可以在硬件、固件、软件或它们的任意组合中实现。

出于解释的目的，已经参考具体实施例对前面的描述进行了描述。然而，上述说明性的讨论并非旨在穷举或将本发明限制于所公开的精确形式。鉴于上述教导许多修改和变化是可能的。实施例被选择并描述以便最好地解释本发明的原理及其实际应用，从而使得本领域技术人员能够以适合于预期的特定用途的各种修改来最好地利用本发明和各种实施例。

Claims (15)

1.一种计算机实现的方法，所述方法包括：

接收与多个源和接收器位置对应的一个或多个地震测量；

提供用于地质介质的地球模型，所述地质介质具有非均匀倾斜的对称轴，其中针对所述地质介质内的至少一个位置子集所述地球模型包括对称轴方向上的非零剪切速度；以及

根据所述地球模型和能量守恒拟声波方程组在多个时间步长上传播所述一个或多个地震测量，所述能量守恒拟声波方程组包括描述一个或多个地震波场的第一等式和描述一个或多个地震波场的第二等式；以及

其中，所述能量守恒拟声波方程组是通过对能量非守恒拟声波方程组中的所述一个或多个地震波场的一个或多个导数项进行近似而从所述能量非守恒拟声波方程组导出的。

2.根据权利要求1所述的方法，其中所述地球模型包括倾斜的横向各向同性模型。

3.根据权利要求1所述的方法，其中一个或多个近似的导数项沿着相应方向包括至少一个导数项，所述至少一个导数项通过相对于相应方向将各向异性参数视为常数而被近似。

4.根据权利要求1所述的方法，其中传播所述一个或多个地震测量包括计算四个二阶导数，包括：

相对于各向异性平面的第一地震波场的二维拉普拉斯导数，所述各向异性平面局部垂直于对称轴；

第一地震波场的三维拉普拉斯导数；

第二地震波场的三维拉普拉斯导数；以及

相对于对称轴方向的第三地震波场的二阶导数；

其中，第三地震波场是第一地震波场和第二地震波场的线性组合。

5.根据权利要求1所述的方法，其中传播所述一个或多个地震波场包括在所述多个时间步长中的每个时间步长计算所述地质介质内的多个位置处的所述一个或多个地震波场。

6.根据权利要求1所述的方法，其中传播所述一个或多个地震测量进一步包括对所述一个或多个地震波场执行逆时偏移。

7.根据权利要求1所述的方法，其中：

所述一个或多个地震波场中的每个在所述多个时间步长中的每个处通过存储在其中的相应能量被表征；以及

对于所述多个时间步长中的每个时间步长，存储在所述一个或多个地震波场中的相应能量的和是常数。

8.一种电子设备，所述电子设备包括：

一个或多个处理器；

存储器；以及

一个或多个程序，其中所述一个或多个程序被存储在存储器中并且被配置为由所述一个或多个处理器执行，所述一个或多个程序包括指令，所述指令当被所述一个或多个处理器执行时使设备：

接收与多个源和接收器位置对应的一个或多个地震测量；

提供用于地质介质的地球模型，所述地质介质具有非均匀倾斜的对称轴，其中针对所述地质介质内的至少一个位置子集所述地球模型包括对称轴方向上的非零剪切速度；以及

根据所述地球模型和能量守恒拟声波方程组在多个时间步长上传播所述一个或多个地震测量，所述能量守恒拟声波方程组包括描述一个或多个地震波场的第一等式和描述一个或多个地震波场的第二等式；以及

其中，所述能量守恒拟声波方程组是通过对能量非守恒拟声波方程组中的所述一个或多个地震波场的一个或多个导数项进行近似而从所述能量非守恒拟声波方程组导出的。

9.根据权利要求8所述的电子设备，其中所述地球模型包括倾斜的横向各向同性模型。

10.根据权利要求8所述的电子设备，其中一个或多个近似的导数项沿着相应方向包括至少一个导数项，所述至少一个导数项通过相对于相应方向将各向异性参数视为常数而被近似。

11.根据权利要求8所述的电子设备，其中传播所述一个或多个地震测量包括计算四个二阶导数，包括：

相对于各向异性平面的第一地震波场的二维拉普拉斯导数，所述各向异性平面局部垂直于对称轴；

第一地震波场的三维拉普拉斯导数；

第二地震波场的三维拉普拉斯导数；以及

相对于对称轴方向的第三地震波场的二阶导数；

其中，第三地震波场是第一地震波场和第二地震波场的线性组合。

12.根据权利要求8所述的电子设备，其中传播所述一个或多个地震波场包括在所述多个时间步长中的每个时间步长计算所述地质介质内的多个位置处的所述一个或多个地震波场。

13.根据权利要求8所述的电子设备，其中传播所述一个或多个地震测量进一步包括对所述一个或多个地震波场执行逆时偏移。

14.根据权利要求8所述的电子设备，其中：

所述一个或多个地震波场中的每个在所述多个时间步长中的每个处通过存储在其中的相应能量被表征；以及

对于所述多个时间步长中的每个时间步长，存储在所述一个或多个地震波场中的相应能量的和是常数。

15.一种存储一个或多个程序的非临时性计算机可读存储介质，所述一个或多个程序包括指令，所述指令当被具有一个或多个处理器和存储器的电子设备执行时使设备：

接收与多个源和接收器位置对应的一个或多个地震测量；

提供用于地质介质的地球模型，所述地质介质具有非均匀倾斜的对称轴，其中针对所述地质介质内的至少一个位置子集所述地球模型包括对称轴方向上的非零剪切速度；以及

根据所述地球模型和能量守恒拟声波方程组在多个时间步长上传播所述一个或多个地震测量，所述能量守恒拟声波方程组包括描述一个或多个地震波场的第一等式和描述一个或多个地震波场的第二等式；以及

其中，所述能量守恒拟声波方程组是通过对能量非守恒拟声波方程组中的所述一个或多个地震波场的一个或多个导数项进行近似而从所述能量非守恒拟声波方程组导出的。