CN105372704B - 一种获取地震波传播方向的方法及装置 - Google Patents

Abstract

本申请实施例公开了一种获取地震波传播方向的方法及装置。该方法包括利用所获取的地震波的单频波场，构造结构张量；计算所述结构张量的特征值；利用所述结构张量的特征值，确定所述地震波的传播方向。通过本申请实施例所提供的方法可以提高确定出地震波传播方向的准确性，还可以提高地震数据的处理效率。

Description

一种获取地震波传播方向的方法及装置

技术领域

本申请涉及石油勘探技术领域，特别涉及一种获取地震波传播方向的方法及装置。

背景技术

地震波的传播方向是地震波的一种重要属性，对地震资料的处理具有重大的意义。目前，现有技术中主要通过射线追踪技术和求解函数方程来实现地震波传播方向的获取。其中，射线追踪主要是指根据地震波的传播规律来确定地震波在实际地层中传播的射线路径，其主要包括打靶法和弯曲法。打靶法主要是通过给定的初始条件来确定震源发出的射线的起始点和终点之间的最短路径。而弯曲法主要是两点边值问题，首先在射线路径的起始点到终点之间给定一条路径，然后对给定的路径进行循环修正，直到修正量很小为止。

在实现本申请过程中，发明人发现现有技术中至少存在如下问题：

当介质比较复杂时，地震波射线在该介质中的传播速度可能会不断变化，从而导致通过所给定路径的射线可能会发生交叉或盲区等现象，这可能会导致无法准确确定出地震波的传播方向，因此可能需要一种新的方法来确定出地震波的传播方向。

发明内容

本申请实施例的目的是提供一种获取地震波传播方向的方法及装置，以提高确定出地震波传播方向的准确性。

为解决上述技术问题，本申请实施例提供一种获取地震波传播方向的方法及装置是这样实现的：

本申请实施例提供了一种获取地震波传播方向的方法，包括：

利用所获取的所述地震波的单频波场，构造结构张量；

计算所述结构张量的特征值；

利用所述结构张量的特征值，确定所述地震波的传播方向。

在一实施例中，所述单频波场通过以下方式获取：

获取目标区域中第一预设范围内的地震波场以及地震波速度数据；

利用所获取的地震波场以及所述地震波速度数据，获取预设频率所对应的单频波场。

在一实施例中，所述获取预设频率所对应的单频波场包括：

利用单程波延拓方法或双程波延拓方法，计算所述第一预设范围内地下不同深度处的单频波场。

在一实施例中，所述计算所述结构张量的特征值包括：

对所述结构张量进行平滑；

计算平滑后的结构张量的特征值。

在一实施例中，所述对所述结构张量进行平滑包括：

将所述结构张量的一个三维数据体中的一坐标点作为当前坐标点；

以所述当前坐标点为中心，计算所述当前坐标点在第二预设范围内的平均值；

将所得到的平均值作为所述当前坐标点的数值；

按照上述计算一坐标点的数值的方法，依次计算所述三维数据体中剩余N-1个坐标点的数值；

按照上述计算一个三维数据体中所有坐标点的数值的方法，依次计算所述结构张量中剩余M-1个三维数据体中所有坐标点的数值。

在一实施例中，所述利用所述结构张量的特征值，确定所述地震波的传播方向包括：

从所述结构张量的多个特征值中选取绝对值最大的特征值；

将所选取的特征值所对应的特征向量确定为所述地震波波前的法线方向。

本申请实施例还提供了一种获取地震波传播方向的装置，包括：

构造单元，用于利用所获取的所述地震波的单频波场，构造结构张量；

计算单元，用于计算所述结构张量的特征值；

确定单元，用于利用所述结构张量的特征值，确定所述地震波的传播方向。

在一实施例中，所述计算单元包括：

平滑子单元，用于对所述结构张量进行平滑；

计算子单元，用于计算平滑后的结构张量的特征值。

在一实施例中，所述平滑子单元具体用于：

将所述结构张量的一个三维数据体中的一坐标点作为当前坐标点；

以所述当前坐标点为中心，计算所述当前坐标点在第二预设范围内的平均值；

将所得到的平均值作为所述当前坐标点的数值；

按照上述计算一坐标点的数值的方法，依次计算所述三维数据体中剩余N-1个坐标点的数值；

按照上述计算一个三维数据体中所有坐标点的数值的方法，依次计算所述结构张量中剩余M-1个三维数据体中所有坐标点的数值。

在一实施例中，所述确定单元包括：

选取子单元，用于从所述结构张量的多个特征值中选取绝对值最大的特征值；

确定子单元，用于将所选取的特征值所对应的特征向量确定为所述地震波波前的法线方向。

由以上本申请实施例提供的技术方案可见，本申请实施例通过利用所获取的所述地震波的单频波场，构造结构张量；计算所述结构张量的特征值；利用所述结构张量的特征值可以得到地震波波前法线方向，从而实现了提高确定出地震波传播方向的准确性的目的。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是本申请中一种获取地震波传播方向的方法的一个实施例的方法流程图。

图2是利用单程波延拓方法所得到的15Hz单频波场的示意图。

图3是上行波的波场延拓示意图。

图4是利用双程波延拓方法所得到的15Hz单频波场的示意图。

图5是单频波场的倾角扫描结果示意图。

图6是图5中所示的倾角扫描结果与单频波场的叠加效果图。

图7是本申请中一种获取地震波传播方向的装置的一个实施例的结构示意图。

具体实施方式

本申请实施例提供一种获取地震波传播方向的方法及装置。

为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。

下面结合附图对本申请所述的获取地震波传播方向的方法进行详细的说明。虽然本申请提供了如下述实施例或流程图所述的方法操作步骤，但基于常规或者无需创造性的劳动在所述方法中可以包括更多或者更少的操作步骤。在逻辑性上不存在必要因果关系的步骤中，这些步骤的执行顺序不限于本申请实施例提供的执行顺序。所述的方法的在实际中的装置或终端产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行。

图1是本申请中一种获取地震波传播方向的方法的一个实施例的方法流程图。该方法包括：

S110：利用所获取的所述地震波的单频波场，构造结构张量。

所述单频波场可以是指单个频率所对应的地震波场。地震波场可以是指地震波传播的空间。在这个空间的每一点上，一定时刻都有一定的波前通过。在这个空间中波的能量也按一定的规律传播，所有这些规律可以是由震源的特点以及在此空间内介质的物理性质(主要是弹性)和几何结构来决定。所述地震波场可以是在对地震数据进行偏移计算的过程中获取的，通常可以用P(t,x,y,z)来确定。其中，t表示地震波的传播时间，x表示x轴方向上地震波的传播距离；y表示y轴方向上地震波的传播距离，z表示z轴方向上地震波的传播距离。在本申请实施例中，z可以表示地下深度。

所述单频波场可以通过以下方式来获取：获取目标区域中第一预设范围内的地震波场以及地震波速度数据；利用所获取的地震波场以及所述地震波速度数据，获取预设频率所对应的单频波场。

所述地震波场可以在对地震数据进行偏移计算的过程中来获取。

所述地震波速度数据可以是指在目标区域中第一预设范围内每一个地下位置处(即不同深度处)的速度值。每一地下位置处的速度可以表示为v(x,y,z)。所述第一预设范围可以是以炮点为中心来设定，也可以以检波点为中心来设定。所述地震波速度数据可以从针对目标区域所建立的速度模型中来选取。

在对地震数据进行偏移计算的过程中，确定出目标区域的地震波场后，可以利用波场延拓方法来计算地下每个深度的单频波场值。所述波场延拓方法可以为单程波延拓方法或双程波延拓方法。

所述单程波延拓方法也可以称为基于单程波方程的偏移方法，其可以是指将声波方程沿深度方向分解为上、下行独立传播的方程，然后可以对上行波或下行波沿深度方向进行递推计算。具体的，可以是指将地下震源波场从炮点位置延拓到反射界面；然后将反射点作为新的震源，依据波场扩散规律，将反射波场上行延拓到地面，最后将同一时间到达检波点的波场进行叠加。所述单程波延拓方法还可以是指将地下震源波场从炮点延拓到反射界面；然后将反射点作为新的震源，将反射波场下行延拓，最后将同一时间到达预设点的波场进行叠加。目前，常用的单程波延拓方法主要包括单程波有限差分法，***步相移法、傅里叶有限差分法以及广义屏法等。利用单程波延拓方法进行递推计算，可以节省计算量和存储资源。

所述双程波延拓方法可以是指利用基于双程波方程的有限差分偏移方法通过在时间域上进行递推计算来求解全声波方程。

在一实施例中，可以利用单程波延拓方法来计算地下每个深度处的单频波场可以包括以下子步骤：

S111：利用所获取的地震波场，构建三维单程波方程。具体的，

在确定出目标区域的地震波场后，利用所述地震波场，构建各向同性完全弹性介质中的标量三维声波方程，所构建的三维声波方程可以表示如下：

上式中，P(t,x,y,z)为波场，v(x,y,z)为声波传播速度。

将上述三维声波方程分别对x、y、z和t进行傅里叶变换。在假设v(x,y,z)为定值时，所得到的传播方程如下：

上式中，k_x、k_y和k_z分别为对应于x、y和z的波数域变量；ω＝2πf。其中，正、负号分别表示上行波和下行波。

对于上行波，对式(2)中的z变量进行傅里叶反变换，可以得到下述三维单程波方程：

对于下行波，所构建的三维单程波方程如下：

S112：利用所获取的地震波速度数据以及所构建的三维单程波方程，计算预设频率所对应的单频波场。

在构建上行波或下行波的三维单声波方程后，可以利用所获取的地震波速度数据，计算预设频率所对应的单频波场。具体的，

在一实施例中，首先可以利用所获取的地震波速度数据，计算出预设范围内每个深度位置处的单程波场；然后根据所得到的每个深度位置处的单程波场，计算出预设频率所对应的单频波场，即将预设频率(例如15Hz)代入到所得到的单程波场中，即可得到预设频率所对应的单频波场。图2示出了利用单程波延拓方法所得到的15Hz单频波场。

所述计算出预设范围内每个深度位置处的单程波场可以包括：首先可以对式(3)进行积分；然后对x和y变量进行傅里叶反变换，得到地面单程波场的表达式；再将所获取的地震波速度数据中地震波在地面上的传播速度代入所得到的表达式中，从而可以计算出地面单程波场，即深度为0时的单程波场，可以用P(x,y,0,t)来表示；最后将所得到的地面单程波场沿z轴反方向进行递推计算，并利用所获取的地震波速度数据中每个深度位置处的速度，计算出预设范围内每个深度位置处的单程波场。图3示出了上行波的波场延拓示意图，图中v1、v2和v3分别表示在深度Δz、2Δz和3Δz处地震波的传播速度。

所述计算出预设范围内每个深度位置处的单程波场也可以包括：首先可以对式(4)进行积分；然后对x和y变量进行傅里叶反变换，得到地面单程波场的表达式；再将所获取的地震波速度数据中地震波在炮点处的传播速度代入所得到的表达式中，从而可以计算出炮点处的单程波场，例如深度为h时的单程波场，可以用P(x,y,h,t)来表示；最后将所得到的单程波场沿z轴正方向进行递推计算，并利用所获取的地震波速度数据中每个深度位置处的速度，计算出预设范围内每个深度位置处的单程波场。

可以利用现有技术中的方法来对地震波场沿z轴反方向或正方向进行递推计算，在此不再赘叙。

在另一实施例中，可以利用双程波延拓方法来计算地下每个深度处的单频波场可以包括以下子步骤：

S111’：利用所获取的地震波场，构建三维双程波方程。

S112’：利用所获取的地震波速度数据以及所构建的三维双程波方程，计算预设频率所对应的单频波场。

首先可以利用所获取的地震波速度数据，计算出预设范围内每个深度位置处的双程波场；然后根据所得到的每个深度位置处的双程波场，计算出预设频率所对应的单频波场。具体的，可以通过将时间域的双程波场乘以指数函数e^-iωt来将双程波波场转换成频率域的双程波场，再将预设频率(例如15Hz)代入到所得到的双程波场中，获取该双程波场的实部数值，即可得到预设频率所对应的单频波场。图4示出了利用双程波延拓方法所得到的15Hz的单频波场。

上述两个步骤的具体执行过程可以参考现有技术中的方法，在此不再赘叙。

利用单程波延拓方法或双程波延拓方法所得到的预设频率所对应的单频波场构成了三维波场数据体。所述三维波场数据体可以理解为一个三维波场，该波场中每个点具有波场值。

在得到第一预设范围内预设频率所对应的单频波场后，可以利用所得到的单频波场构造结构张量。

首先，可以计算所述单频波场中每个位置(即每个坐标点)处的梯度，所得到的梯度可以表示如下：

上式中，为单频波场I的梯度。当x,y和z取不同的数值时，即可得到单频波场在各个点处的梯度。所述单频波场中每个坐标点处的梯度为1x3(1行3列)维度的矩阵。矩阵中的每个元素分别是单频波场在x,y或z方向上的偏导数，每个元素均可以是一个三维数据体，即每个元素均可以是以三维坐标(x,y,z)为位置函数的数据集合体。

然后，可以利用所得到的单频波场的梯度，构造结构张量。具体的可以通过对所述单频波场的梯度进行转置运算来构造结构张量，即将所述梯度的转置与所述梯度进行相乘，得到结构张量。

所构造的结构张量T可以表示如下：

上述结构张量T是一个3x3维度的矩阵，矩阵中的每个元素也均是三维数据体，该三维数据体中每一坐标点的数值为单频波场在x,y和z方向的偏导数两两之间的乘积。

所述结构张量可以表示目标区域的变化方向和沿变化方向的变化量大小。其中，所述结构张量的特征向量可以反映目标区域中局部区域变化的方向，其特征值可以反映所述变化量的大小。对于三维数据，所述结构张量可以为半正定矩阵，该半正定矩阵的主特征向量可以表示正交于地层的梯度方向，另外两个特征向量可以表示平行于地层方向的平面方向。

S120：计算所述结构张量的特征值。

在构造出结构张量后，可以首先对所述结构张量进行平滑，然后可以计算平滑后的结构张量的特征值。

所述对所述结构张量进行平滑可以是指对所述结构张量中每个元素的数值进行平滑，具体的可以包括：将所述矩阵的一个三维数据体中的一坐标点作为当前坐标点；然后以所述当前坐标点为中心，计算所述当前坐标点在第二预设范围内的平均值，例如计算该当前坐标点的数值与其周围相邻坐标点的数值之间的平均值；再将所得到的平均值作为所述当前坐标点的数值；接着按照上述计算一坐标点的数值的方法，依次计算所述三维数据体中剩余N-1个坐标点的数值；最后按照上述计算一个三维数据体中所有坐标点的数值的方法，依次计算所述结构张量中剩余M-1个三维数据体中所有坐标点的数值。其中，N为一个三维数据体中坐标点的个数；M为三维数据体的个数，在一实施例中M为9。

对所述结构张量进行平滑可以用公式表示如下：

上式中，ST为平滑后的结构张量；求和符号∑代表平滑计算。

所述结构张量为3x3的矩阵，因而该结构张量有三个特征值。可以利用现有的数学方法来计算结构张量的特征值，在此不再赘叙。

S130：确定所述地震波的传播方向。

在计算出所述结构张量的特征值后，可以将所有的特征值进行对比，从所有的特征值中选取出绝对值最大的特征值。将所选取的绝对值最大的特征值所对应的特征向量为主特征向量。该主特征向量是地震波波前的法线方向，可以用两个xdip、ydip来描述该方向，其中xdip是地震波波前法线在xz平面的投影与z轴的夹角，ydip是地震波波前法线在yz平面的投影与z轴的夹角。

由此可见，在得到所述结构张量的特征值后，可以确定出地震波波前法线在xz平面的投影与z轴的夹角以及该法线在yz平面的投影与z轴的夹角，即确定出所述地震波的传播方向。

上述利用单频波场来构造结构张量，并对结构张量进行平滑以及求取特征值的过程可以简称为利用结构张量倾角扫描法对单频波场进行扫描。

图5-图6分别示出了单频波场的倾角扫描结果的示意图以及所述倾角扫描结果与单频波场的叠合效果示意图。从这两幅图中可以看出利用结构张量对单频波场处理后所得到的地震波的走势与所述单频波场中地震波的传播方向相吻合，由此可见利用本申请实施例提供的方法可以准确的确定出地震波的传播方向。

通过上述描述可以看出，本申请实施例通过利用所获取的地震波场以及地震波速度数据，获取预设频率所对应的单频波场；利用所述单频波场，构造结构张量；获取所述结构张量的特征值，利用所述特征值可以得到地震波波前法线在xz平面的投影与z轴的夹角以及该法线在yz平面的投影与z轴的夹角，从而实现了提高确定出地震波传播方向的准确性的目的。而且，本申请实施例中在获取预设频率所对应的单频波场时，在进行递推计算的过程中所利用的地震波传播速度可以是变化的，这正好可以与实际地下介质中速度不断变化的情形吻合，因而可以进一步保证所确定出的地震波传播方向的准确性。此外，本申请实施例中利用的是叠前偏移计算过程中所得到的地震波场来计算地震波的方向，而不是在进行叠前偏移计算后，又重新求解波动方程，因此利用本申请实施例所提供的方法，可以提高地震数据的处理效率，还有利于后续入射角道集的提取。

本申请实施例还提供了一种获取地震波传播方向的装置，如图7所示。该装置可以包括：构造单元710、计算单元720以及确定单元730。其中，构造单元710可以用于利用所获取的所述地震波的单频波场，构造结构张量。计算单元720可以用于计算所述结构张量的特征值；确定单元730可以用于利用所述结构张量的特征值，确定所述地震波的传播方向。

在一实施例中，计算单元720可以包括(图中未示出)：

平滑子单元，用于对所述结构张量进行平滑。该平滑单元具体的可以用于将所述结构张量的一个三维数据体中的一坐标点作为当前坐标点；以所述当前坐标点为中心，计算所述当前坐标点在第二预设范围内的平均值；将所得到的平均值作为所述当前坐标点的数值；按照上述计算一坐标点的数值的方法，依次计算所述三维数据体中剩余N-1个坐标点的数值；按照上述计算一个三维数据体中所有坐标点的数值的方法，依次计算所述结构张量中剩余M-1个三维数据体中所有坐标点的数值。

计算子单元，用于计算平滑后的结构张量的特征值。

在一实施例中，确定单元730可以包括(图中未示出)：

选取子单元，用于从所述结构张量的多个特征值中选取绝对值最大的特征值；

确定子单元，用于将所选取的特征值所对应的特征向量确定为所述地震波波前的法线方向。

上述实施例阐明的***、装置或单元，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。

为了描述的方便，描述以上装置时以功能分为各种单元分别描述。当然，在实施本申请时可以把各单元的功能在同一个或多个软件和/或硬件中实现。

通过以上的实施方式的描述可知，本领域技术人员还可以了解到本发明实施例列出的各种说明性逻辑块、单元和步骤可以通过硬件、软件或两者的结合来实现。至于是通过硬件还是软件来实现取决于特定的应用和整个***的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。

本发明实施例中所描述的各种说明性的逻辑块或单元都可以通过通用处理器，数字信号处理器，专用集成电路(AS工C)，现场可编程门阵列或其它可编程逻辑装置，离散门或晶体管逻辑，离散硬件部件，或上述任何组合的设计来实现或操作所描述的功能。通用处理器可以为微处理器，可选地，该通用处理器也可以为任何传统的处理器、控制器、微控制器或状态机。处理器也可以通过计算装置的组合来实现，例如数字信号处理器和微处理器，多个微处理器，一个或多个微处理器联合一个数字信号处理器核，或任何其它类似的配置来实现。

本发明实施例中所描述的方法或算法的步骤可以直接嵌入硬件、处理器执行的软件模块、或者这两者的结合。软件模块可以存储于RAM存储器、闪存、ROM存储器、EPROM存储器、EEPROM存储器、寄存器、硬盘、可移动磁盘、CD-ROM或本领域中其它任意形式的存储媒介中。示例性地，存储媒介可以与处理器连接，以使得处理器可以从存储媒介中读取信息，并可以向存储媒介存写信息。可选地，存储媒介还可以集成到处理器中。处理器和存储媒介可以设置于ASIC中，ASIC可以设置于用户终端中。可选地，处理器和存储媒介也可以设置于用户终端中的不同的部件中。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于***实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。

虽然通过实施例描绘了本申请，本领域普通技术人员知道，本申请有许多变形和变化而不脱离本申请的精神，希望所附的权利要求包括这些变形和变化而不脱离本申请的精神。

Claims (6)

1.一种获取地震波传播方向的方法，其特征在于，包括：

利用所获取的地震波的单频波场，构造结构张量，所述结构张量中的每个元素均是三维数据体；

计算所述结构张量的特征值；

利用所述结构张量的特征值，确定所述地震波的传播方向，

其中，计算所述结构张量的特征值包括：

对所述结构张量进行平滑；

计算平滑后的结构张量的特征值；

其中，对所述结构张量进行平滑包括：

将所述结构张量的一个三维数据体中的一坐标点作为当前坐标点；

以所述当前坐标点为中心，计算所述当前坐标点在第二预设范围内的平均值；

将所得到的平均值作为所述当前坐标点的数值；

按照上述计算一坐标点的数值的方法，依次计算所述三维数据体中剩余N-1个坐标点的数值，N为一个三维数据体中坐标点的个数；

按照上述计算一个三维数据体中所有坐标点的数值的方法，依次计算所述结构张量中剩余M-1个三维数据体中所有坐标点的数值，M为三维数据体的个数。

2.根据权利要求1所述的方法，其特征在于，所述单频波场通过以下方式获取：

获取目标区域中第一预设范围内的地震波场以及地震波速度数据；

利用所获取的地震波场以及所述地震波速度数据，获取预设频率所对应的单频波场。

3.根据权利要求2所述的方法，其特征在于，所述获取预设频率所对应的单频波场包括：

利用单程波延拓方法或双程波延拓方法，计算所述第一预设范围内地下不同深度处的单频波场。

4.根据权利要求1所述的方法，其特征在于，所述利用所述结构张量的特征值，确定所述地震波的传播方向包括：

从所述结构张量的多个特征值中选取绝对值最大的特征值；

将所选取的特征值所对应的特征向量确定为所述地震波波前的法线方向。

5.一种获取地震波传播方向的装置，其特征在于，包括：

构造单元，用于利用所获取的所述地震波的单频波场，构造结构张量，所述结构张量中的每个元素均是三维数据体；

计算单元，用于计算所述结构张量的特征值；

确定单元，用于利用所述结构张量的特征值，确定所述地震波的传播方向，

其中，所述计算单元包括：

平滑子单元，用于对所述结构张量进行平滑；

计算子单元，用于计算平滑后的结构张量的特征值；

其中，所述平滑子单元具体用于：

将所述结构张量的一个三维数据体中的一坐标点作为当前坐标点；

以所述当前坐标点为中心，计算所述当前坐标点在第二预设范围内的平均值；

将所得到的平均值作为所述当前坐标点的数值；

按照上述计算一坐标点的数值的方法，依次计算所述三维数据体中剩余N-1个坐标点的数值；

按照上述计算一个三维数据体中所有坐标点的数值的方法，依次计算所述结构张量中剩余M-1个三维数据体中所有坐标点的数值。

6.根据权利要求5所述的装置，其特征在于，所述确定单元包括：

选取子单元，用于从所述结构张量的多个特征值中选取绝对值最大的特征值；

确定子单元，用于将所选取的特征值所对应的特征向量确定为所述地震波波前的法线方向。