一种地震勘探中的能流密度矢量计算方法

1.本发明涉及地震波油气勘探领域，特别是涉及一种地震勘探中的能流密度矢量计算方法。


背景技术：

2.叠前深度偏移成像是用于对地下地质结构进行成像、帮助油气勘探和生产的地震数据处理技术，虽然偏移成像技术是目前工业界应用最为广泛的技术之一，但由于该技术受限于输入速度场精度的影响，因此目前行业内的重要研究方向之一就是地下速度场的构建。随着地震勘探技术的不断发展，速度建模方法逐渐从常规交互式nmo速度分析方法转向偏移速度分析方法，在众多的速度建模方法中，基于角度域共成像点道集(角道集)的层析方法得到广泛的应用。其中，基于方向矢量方法的角道集计算方法由于理论简单、易于实现而逐渐在工业界得到应用。
3.目前已有多种方向矢量方法可以计算地震波的传播角度进而计算角道集，其中，基于能流密度矢量(poynting矢量)的计算方法(yoon k.和marfurt k.，2006，reverse-time migration using the poynting vector.exploration geophysics，37(1)：102-107)理论明确而逐渐得到应用。但是本技术的发明人发现：由于地震波在地下传播时较为复杂，易出现在同一空间成像点不同时刻波场混叠的问题，因此造成能流密度计算不够准确，进而影响最终的速度建模精度。


技术实现要素：

4.本发明的目的在于提供一种地震勘探中的能流密度矢量计算方法，通过扫描地下每一空间点的波场振幅值，然后获得准确度更高的能流密度矢量，以此为基础可获得高精度的角度域共成像点道集。
5.为实现上述目的，本发明提供了一种地震勘探中的能流密度矢量计算方法，其包括以下步骤：
6.步骤一，基于波动方程对波场进行正向延拓，获得震源波场的位移场；
7.步骤二，将步骤一保存的震源波场的位移场进行梯度运算获得震源波场的应力场，将步骤一保存的震源波场的位移场进行传播时间的导数运算获得震源波场的速度场；
8.步骤三，在传播时间上扫描所述步骤二的震源波场的速度场，并保存每一个空间点在速度场振幅模值最大时的应力场和速度场；
9.步骤四，根据能流密度的计算方法获得波场的能流密度矢量。
10.作为本发明优选的方案，还包括以下验证步骤，
11.步骤五，利用地层倾角和步骤四中获得的震源波场的能流密度矢量，根据波场传播方向计算公式计算波场的入射角；
12.步骤六，基于波动方程对波场进行反向延拓，获得检波波场的位移场；
13.步骤七，利用所述震源波场和检波波场位移场、所述波场的入射角以及角度域共
成像点道集计算公式计算角度域共成像点道集，对比理论入射角与角度域共成像点道集中的入射角，以验证所述步骤四计算的能流密度矢量的计算准确性。
14.作为本发明优选的方案，所述步骤一中，纵波介质速度场基于二阶声波波动方程，以获得不同正向延拓时间的位移场并保存每一个空间点不同传播时间的震源波场的位移场ps(t,x)，
15.采用如下二阶声波波动方程计算震源波场的位移场ps(t,x)：
[0016][0017]
式中，x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；t为波场的传播时间；p为波场的位移场，上标s为震源波场；v
p
(x)为位于空间点x处的纵波速度。
[0018]
作为本发明优选的方案，所述步骤二中，
[0019]
采用如下震源波场的位移场梯度运算公式计算震源波场的应力场：
[0020][0021]
式中式中为空间点x处的位移场在x方向上的导数，为空间点x处的位移场在z方向上的导数；为空间点x处的震源波场应力场；x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；p为波场的位移场，上标s为震源波场；t为波场的传播时间；为梯度算子。
[0022]
采用如下公式计算震源波场的位移场关于波场的传播时间t的导数：
[0023][0024]
式中，是位移场模值最大时关于震源波场的传播时间波场的传播时间的导数，p为波场的位移场，上标s为震源波场；t为波场的传播时间，δt为波场的传播时间的采样间隔；
[0025]
作为本发明优选的方案，所述步骤三中，
[0026]
采用如下位移场最大振幅模值公式计算每一个地下空间点在所有传播时间上的最大振幅模值同时保存出现最大振幅模值的传播时间：
[0027][0028]
式中，为位于空间点x处的位移场模的最大值；x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；| |表示求模运算，max表示求最大值运算，t
max
为最大传播时间，δt为波场的传播时间的采样间隔。
[0029]
作为本发明优选的方案，所述步骤四中，
[0030]
采用如下公式计算地震波场中的能流密度矢量：
[0031]
[0032]
式中，pv表示地震波场中的能流密度矢量，v表示波场传播的质点振动速度，τ表示应力；
[0033]
采用如下公式计算震源波场中的能流密度矢量：
[0034][0035]
式中，为震源波场在空间点x处的能流密度矢量沿x方向的分量，为震源波场在空间点x处的能流密度矢量沿z方向的分量；为空间点x处的位移场模的最大值在x方向上的导数；为空间点x处的位移场模的最大值在z方向上的导数；为空间点x处的位移场模的最大值关于传播时间t的导数；x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；p为波场的位移场，上标s为震源波场；t为波场的传播时间。
[0036]
作为本发明优选的方案，所述步骤五中，
[0037]
采用如下波场传播方向计算公式计算波场到达成像点处的入射角：
[0038][0039]
式中，α(x)表示波场到达成像点x处的入射角；为震源波场在空间点x处的能流密度矢量沿x方向的分量，为震源波场在空间点x处的能流密度矢量沿z方向的分量；θ(x)为地层倾角。
[0040]
作为本发明优选的方案，所述步骤六中，
[0041]
将地面记录的地震数据作为边值条件，纵波介质速度场基于二阶声波波动方程进行反向时间延拓，并保存每一个空间点不同传播时间的检波波场的位移场pr(t,x)，
[0042]
采用如下二阶声波波动方程确定检波波场的位移场pr(t,x)：
[0043][0044]
式中，x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；t为波场的传播时间；p为波场的位移场，上标r为检波波场；v
p
(x)为位于空间点x处的纵波速度。
[0045]
作为本发明优选的方案，所述步骤七中，
[0046]
采用如下角度域共成像点道集计算公式计算角度域共成像点道集：
[0047][0048]
式中，i
adcigs
表示角度域共成像点道集，x表示笛卡尔直角坐标系中一个空间点的坐标，t为波场的传播时间，t
max
为最大传播时间；θ(x)为波场到达成像点为x处的入射角，θk为角度域共成像点道集中的入射角，δ表示δ函数；p为波场的位移场，上标s表示震源波场，上标r表示检波波场，*表示卷积运算符。
[0049]
作为本发明优选的方案，所述步骤七中，
[0050]
选取水平位置1.75km处的角度域共成像点道集，通过snell’s定律和三角关系计算位于水平位置1.75km处的理论入射角，并通过角度域共成像点道集计算公式提取水平位置1.75km处的角度域共成像点道集，对比理论入射角与通过角度域共成像点道集计算公式中的得出的角度域共成像点道集中的入射角，以验证所述步骤四中的能流密度矢量计算的准确性。
[0051]
本发明实施例一种地震勘探中的能流密度矢量计算方法与现有技术相比，其有益效果在于：
[0052]
本发明通过扫描地下每一空间点的波场振幅值，然后获得准确度更高的能流密度矢量，以高准确度的能流密度矢量为基础能够获得高精度的角度域共成像点道集。
附图说明
[0053]
为了更清楚地说明本发明实施例的技术方案，下面将对实施例的附图作简单地介绍。
[0054]
图1是本发明实施例提供的步骤流程图；
[0055]
图2是本发明实施例提供根据水平和倾斜界面构建的纵波介质速度场模型的示意图；
[0056]
图3是扫描图1中的纵波介质速度场获得的反射界面地层倾角的示意图；
[0057]
图4是正向延拓至750ms时的震源波场位移场的示意图；
[0058]
图5是正向延拓至750ms时的震源波场振幅值最大的速度场的示意图；
[0059]
图6a是正向延拓至750ms时的能流密度矢量的沿x方向的分量的示意图；
[0060]
图6b是正向延拓至750ms时的能流密度矢量的沿z方向的分量的示意图；
[0061]
图7是根据所求能流密度矢量以及地层倾角获得的每一空间点的成像角度的示意图；
[0062]
图8是逆向延拓至750ms时的检波波场位移场的示意图；
[0063]
图9a是标有炮点、检波点、射线路径以及理论计算的波场入射角的纵波介质速度场的示意图；
[0064]
图9b是标有偏移计算的波场入射角以及理论计算的波场入射角的角道集的示意图。
具体实施方式
[0065]
下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
[0066]
在本发明的描述中，需要理解的是，术语“上”、“下”、“左”、“右”、“顶”、“底”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。应当理解的是，本发明中采用术语“第一”、“第二”等来描述各种信息，但这些信息不应限于这些术语，这些术语仅用来将同一类型的信息彼此区分开。例如，在不脱离本发明范围的情况下，“第一”信息也可以被称为“第二”信息，类似的，“第二”信息也可以被称为“第一”信息。
[0067]
如图1至图9b所示，本发明实施例优选实施例的一种地震勘探中的能流密度矢量计算方法，包括以下步骤：
[0068]
步骤一，基于波动方程对波场进行正向延拓，获得震源波场的位移场(s1)；在本实施例中，具体以图2所示的纵波介质速度场模型为例，基于二阶声波波动方程，以获得不同正向延拓时间的位移场并保存每一个空间点不同传播时间的震源波场的位移场ps(t,x)，
[0069]
采用如下二阶声波波动方程计算震源波场的位移场ps(t,x)：
[0070][0071]
式中，x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；t为波场的传播时间；p为波场的位移场，上标s为震源波场；v
p
(x)为位于空间点x处的纵波速度。这样能够获得正向延拓至750ms时刻基于二阶声波波动方程得到如图4所示的震源波场的位移场；
[0072]
步骤二，将步骤一保存的震源波场的位移场进行梯度运算获得震源波场的应力场，将步骤一保存的震源波场的位移场进行传播时间的导数运算获得震源波场的速度场(s2)；
[0073]
具体的，采用如下震源波场的位移场梯度运算公式计算震源波场的应力场：
[0074][0075]
式中，式中，为空间点x处的位移场在x方向上的导数，为空间点x处的位移场在z方向上的导数；为空间点x处的震源波场应力场；x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；p为波场的位移场，上标s为震源波场；t为波场的传播时间；为梯度算子；
[0076]
采用如下公式计算震源波场的位移场关于波场的传播时间t的导数：
[0077][0078]
式中，是位移场模值最大时关于波场的传播时间的导数，p为波场的位移场，上标s为震源波场；t为波场的传播时间。
[0079]
步骤三，在传播时间上扫描所述步骤二的震源波场的速度场，并保存每一个空间点在速度场振幅模值最大时的应力场和速度场(s3)；
[0080]
具体的，采用如下位移场最大振幅模值公式计算每一个地下空间点在所有传播时间上的最大振幅模值同时保存出现最大振幅模值的传播时间t
max
：
[0081][0082]
式中，为位于空间点x处的位移场模的最大值；x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；| |表示求模运算，max表示求最大值运算，t
max
为最大传播时间，t为波场的传播时间，δt为波场的传播时间的采样间隔。在本实施例中，这样能够获得在在正向延拓至750ms时出现振幅模最大值时的地下
空间点(如图5所示)。
[0083]
步骤四，根据能流密度的计算方法获得波场的能流密度矢量；
[0084]
具体的，采用如下公式计算地震波场中的能流密度矢量：
[0085][0086]
式中，pv表示地震波场中的能流密度矢量，v表示波场传播的质点振动速度，τ表示应力；
[0087]
采用如下公式计算震源波场中的能流密度矢量：
[0088][0089]
式中，为震源波场在空间点x处的能流密度矢量沿x方向的分量，为震源波场在空间点x处的能流密度矢量沿z方向的分量；为空间点x处的位移场模的最大值在x方向上的导数；为空间点x处的位移场模的最大值在z方向上的导数；为空间点x处的位移场模的最大值关于传播时间t的导数；x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z为垂直方向的坐标；p为波场的位移场，上标s为震源波场；t为波场的传播时间。在本实施例中，这样能够获得在正向延拓至750ms时出现模最大值时的空间点处的能流密度矢量沿x方向的分量(如图6a所示)，以及，在正向延拓至750ms时出现模最大值时的空间点处的能流密度矢量沿z方向的分量(如图6b所示)。
[0090]
示例性的，所述地震勘探中的能流密度矢量计算还包括以下验证步骤，以验证波场的能流密度矢量计算的准确性，
[0091]
步骤五，利用已知的地层倾角(如图3所示)和步骤四中获得的震源波场的能流密度矢量，根据波场传播方向计算公式计算波场的入射角(s4)；具体的，所述步骤五中，
[0092]
采用如下波场传播方向计算公式计算波场到达成像点处的入射角：
[0093][0094]
式中，α(x)表示波场到达成像点x处的入射角；为震源波场在空间点x处的能流密度矢量沿x方向的分量，为震源波场在空间点x处的能流密度矢量沿z方向的分量；θ(x)为地层倾角。这样能够获得整个空间上各成像点x上的入射角(如图7所示)。
[0095]
步骤六，基于波动方程对波场进行反向延拓，获得检波波场的位移场(s5)；在本实施例中，具体以图2所示的纵波介质速度场模型为例，将地面记录的地震数据作为边值条件，基于二阶声波波动方程进行反向时间延拓，并保存每一个空间点不同传播时间的检波波场的位移场pr(t,x)，
[0096]
采用如下二阶声波波动方程确定检波波场的位移场pr(t,x)：
[0097][0098]
式中，x＝x,z表示笛卡尔直角坐标系中一个空间点的坐标，x为水平方向的坐标，z
为垂直方向的坐标；t为波场的传播时间；p为波场的位移场，上标r为检波波场；v
p
(x)为位于空间点x处的纵波速度。这样能够获得在反向延拓至750ms时刻基于二阶声波波动方程得到如图8所示的检波波场的位移场。
[0099]
步骤七，利用所述震源波场和检波波场位移场、所述波场的入射角以及角度域共成像点道集计算公式计算角度域共成像点道集，对比理论入射角与角度域共成像点道集中的入射角，以验证所述步骤四计算的能流密度矢量的计算准确性(s6)。具体的，采用如下角度域共成像点道集计算公式计算角度域共成像点道集：
[0100][0101]
式中，i
adcigs
表示角度域共成像点道集，x表示笛卡尔直角坐标系中一个空间点的坐标，t为波场的传播时间，t
max
为最大传播时间；θ(x)为波场到达成像点为x处的入射角，θk为角度域共成像点道集中的入射角，δ表示δ函数；p为波场的位移场，上标s表示震源波场，上标r表示检波波场，*表示卷积运算符。
[0102]
示例性的，在所述步骤七中，选取水平位置1.75km处的角度域共成像点道集，通过snell’s定律和三角关系计算位于水平位置1.75km处的理论入射角(如图9a所示)，并通过角度域共成像点道集计算公式提取水平位置1.75km处的角度域共成像点道集，对比理论入射角与通过角度域共成像点道集计算公式中的得出的角度域共成像点道集中的入射角，以验证所述步骤四中的能流密度矢量计算的准确性。当理论入射角与计算得出的角度域共成像点道集中的入射角相同，则上述能流密度矢量计算准确。其中，在所述步骤七中，可根据实际情况选取水平位置任一处的角度域共成像点道集，再进行后续测试。
[0103]
本发明通过扫描地下每一空间点的波场振幅值，然后获得准确度更高的能流密度矢量，以高准确度的能流密度矢量为基础能够获得高精度的角度域共成像点道集。
[0104]
并且，本发明还提供了地震勘探中的能流密度矢量计算装置，震源波场位移场获取模块，用于对波场基于波动方程进行正向延拓，获得震源波场的位移场；
[0105]
应力场以及速度场获取模块，用于对与保存的位移场进行梯度运算和关于传播时间的导数运算，分别获震源波场的应力场和速度场；
[0106]
空间点最大振幅模值应力场以及速度场获取模块，用于在传播时间上扫描所述震源波场速度场，并保存每一个空间点速度场振幅模值最大时的应力场和速度场；
[0107]
能流密度矢量获取模块，用于根据能流密度的计算方法获得波场的能流密度矢量；
[0108]
波场入射角获取模块，用于利用地层倾角和能流密度矢量，根据波场传播方向计算公式计算波场的入射角；
[0109]
检波波场位移场获取模块，用于对波场基于波动方程进行逆向延拓，获得检波波场的位移场；
[0110]
角度域共成像点道集获取模块，用于利用震源波场的位移场和检波波场的位移场、波场的入射角以及角度域共成像点道集计算公式计算角度域共成像点道集，并对比理论入射角与角度域共成像点道集中的入射角。
[0111]
此外，本发明还提供了地震勘探中的能流密度矢量计算设备，其包括处理器、存储器，所述存储器用于存放至少一可执行指令，所述可执行指令使所述处理器执行所述的能
流密度矢量的计算方法对应的操作，具体的操作步骤为：
[0112]
基于波动方程对波场进行正向延拓，获得震源波场的位移场；
[0113]
对与所述保存的位移场进行梯度运算和关于传播时间的导数运算，分别获震源波场的应力场和速度场；
[0114]
在传播时间上扫描所述震源波场速度场，并保存每一个空间点速度场振幅模值最大时的应力场和速度场；
[0115]
根据能流密度的计算方法获得波场的能流密度矢量；
[0116]
利用地层倾角和所述能流密度矢量，根据波场传播方向计算公式计算波场的入射角；
[0117]
基于波动方程对波场进行逆向延拓，获得检波波场的位移场；
[0118]
利用所述震源波场和检波波场位移场、所述波场的入射角以及角度域共成像点道集计算公式计算角度域共成像点道集，并与理论值入射角进行比较。
[0119]
在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
[0120]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和替换，这些改进和替换也应视为本发明的保护范围。