CN113126159B - 一种地震波场坡印廷矢量确定方法及系统 - Google Patents

一种地震波场坡印廷矢量确定方法及系统 Download PDF

Info

Publication number
CN113126159B
CN113126159B CN202110436362.5A CN202110436362A CN113126159B CN 113126159 B CN113126159 B CN 113126159B CN 202110436362 A CN202110436362 A CN 202110436362A CN 113126159 B CN113126159 B CN 113126159B
Authority
CN
China
Prior art keywords
poynting vector
time
seismic
vector
time derivative
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110436362.5A
Other languages
English (en)
Other versions
CN113126159A (zh
Inventor
李志远
梁光河
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Geology and Geophysics of CAS
Original Assignee
Institute of Geology and Geophysics of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Geology and Geophysics of CAS filed Critical Institute of Geology and Geophysics of CAS
Priority to CN202110436362.5A priority Critical patent/CN113126159B/zh
Publication of CN113126159A publication Critical patent/CN113126159A/zh
Priority to US17/725,971 priority patent/US20220350043A1/en
Application granted granted Critical
Publication of CN113126159B publication Critical patent/CN113126159B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • G01V1/366—Seismic filtering by correlation of seismic signals
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/63—Seismic attributes, e.g. amplitude, polarity, instant phase

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种地震波场坡印廷矢量确定方法及系统。所述方法包括：对地震波场的时间导数的幅值进行调整，计算调整后的地震波场的时间导数的坡印廷矢量，得到第一坡印廷矢量；第一坡印廷矢量的幅值与第二坡印廷矢量的幅值的差值处于设定范围内；第二坡印廷矢量为地震波场的坡印廷矢量；对第二坡印廷矢量和第一坡印廷矢量进行运算，得到地震波场的最终坡印廷矢量。本发明解决坡印廷矢量计算不稳定的问题。

Description

一种地震波场坡印廷矢量确定方法及系统
技术领域
本发明涉及地震勘探领域，特别是涉及一种地震波场坡印廷矢量确定方法及系统。
背景技术
波印廷矢量与地震波场的传播矢量方向一致，因为计算效率高，所以在地震勘探领域得到广泛的应用。
坡印廷矢量会有计算不稳定的问题，这是因为坡印廷矢量需要计算地震波场的空间和时间导数，而地震波场中地震子波波峰和波谷处的导数值为零，因此坡印廷矢量无法获得这些点处的传播方向，如图1所示的地震子波及其一阶时间导数，图1中圆圈标注的位置为计算坡印廷矢量不稳定的位置，这些位置会影响后续的数据处理，比如会降低角度域共成像点道集的提取质量。
发明内容
基于此，有必要提供一种地震波场坡印廷矢量确定方法及系统，以解决坡印廷矢量计算不稳定的问题。
为实现上述目的，本发明提供了如下方案：
一种地震波场坡印廷矢量确定方法，包括：
对地震波场的时间导数的幅值进行调整，计算调整后的地震波场的时间导数的坡印廷矢量，得到第一坡印廷矢量；所述第一坡印廷矢量的幅值与第二坡印廷矢量的幅值的差值处于设定范围内；所述第二坡印廷矢量为所述地震波场的坡印廷矢量；
对所述第二坡印廷矢量和所述第一坡印廷矢量进行运算，得到所述地震波场的最终坡印廷矢量。
可选的，所述对地震波场的时间导数的幅值进行调整，计算调整后的地震波场的时间导数的坡印廷矢量，得到第一坡印廷矢量，具体包括：
计算所述地震波场的一阶时间导数；
计算所述地震波场的一阶时间导数的幅值调整因子；
采用所述幅值调整因子调整所述地震波场的一阶时间导数的幅值，得到调整时间导数；
计算所述调整时间导数的坡印廷矢量，得到所述第一坡印廷矢量。
可选的，所述对所述第二坡印廷矢量和所述第一坡印廷矢量进行运算，得到所述地震波场的最终坡印廷矢量，具体包括：
将所述第二坡印廷矢量和所述第一坡印廷矢量相加，得到所述最终坡印廷矢量。
可选的，所述计算所述地震波场的一阶时间导数，具体包括：
采用有限差分法计算所述地震波场的一阶时间导数。
可选的，所述幅值调整因子的计算公式为：
Figure BDA0003033256130000021
其中，a为幅值调整因子；fd为地震波场的主频；Δt表示有限差分的时间间隔。
可选的，所述调整时间导数的计算公式为：
Figure BDA0003033256130000022
其中，
Figure BDA0003033256130000023
为t时刻的调整时间导数；a为幅值调整因子；t为时间；ut为t时刻的地震波场；
Figure BDA0003033256130000024
为对ut求时间偏导数得到的一阶时间导数。
可选的，所述第一坡印廷矢量的计算公式为：
Figure BDA0003033256130000025
其中，
Figure BDA0003033256130000026
为第一坡印廷矢量；
Figure BDA0003033256130000027
为t时刻的调整时间导数；t为时间；
Figure BDA0003033256130000028
为
Figure BDA0003033256130000029
的梯度。
可选的，所述第一坡印廷矢量的计算公式为：
Figure BDA0003033256130000031
其中，
Figure BDA0003033256130000032
为第一坡印廷矢量；ut为t时刻的调整时间导数；t为时间；
Figure BDA0003033256130000033
为对
Figure BDA0003033256130000034
求梯度；ut为t时刻的地震波场；fd为地震波场的主频；Δt表示有限差分的时间间隔。
可选的，所述地震波场的一阶时间导数的计算公式为：
Figure BDA0003033256130000035
其中，ut为t时刻的地震波场；t为时间；Δt表示有限差分的时间间隔；ut+Δt为t+Δt时刻的地震波场；ut-Δt为t-Δt时刻的地震波场。
本发明还提供了一种地震波场坡印廷矢量确定系统，包括：
幅值调整模块，用于对地震波场的时间导数的幅值进行调整，计算调整后的地震波场的时间导数的坡印廷矢量，得到第一坡印廷矢量；所述第一坡印廷矢量的幅值与第二坡印廷矢量的幅值的差值处于设定范围内；所述第二坡印廷矢量为所述地震波场的坡印廷矢量；
坡印廷矢量运算模块，用于对所述第二坡印廷矢量和所述第一坡印廷矢量进行运算，得到所述地震波场的最终坡印廷矢量。
与现有技术相比，本发明的有益效果是：
本发明提出了一种地震波场坡印廷矢量确定方法及系统，其对地震波场的时间导数的幅值进行调整，使得调整后的地震波场的时间导数的坡印廷矢量与地震波场的坡印廷矢量的幅值的差值处于设定范围内，之后对调整后的地震波场的时间导数的坡印廷矢量和地震波场的坡印廷矢量进行运算，得到最终的地震波场的坡印廷矢量。地震波场的坡印廷矢量和其时间导数的坡印廷矢量所表示的传播方向是一致的，幅值调整并不会改变坡印廷矢量的方向。而由于地震波场的坡印廷矢量及其时间导数的坡印廷矢量的幅值变化趋势互补，一个的极大值对应另一个的极小值，因此两者运算便可以抵消极小值带来的计算不稳定，提高了坡印廷矢量计算的稳定性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为地震子波及其一阶时间导数的时间-幅值示意图；
图2为本发明实施例提供的地震波场坡印廷矢量确定方法的流程图；
图3为本发明实施例提供的地震波场及其一阶时间导数的坡印廷矢量的时间-幅值示意图；
图4为均匀介质中的0.4s时刻的地震波场示意图；
图5为对图4所示的地震波场利用原始的坡印廷矢量方程、空间光滑法、光流法、时移法和本发明所计算出的坡印廷矢量示意图；其中，图5(a)为利用原始坡印廷矢量计算方程获得的坡印廷矢量的x分量示意图，图5(b)为利用原始坡印廷矢量计算方程获得的坡印廷矢量的z分量示意图，图5(c)为利用平滑法计算的坡印廷矢量的x分量示意图，图5(d)为利用平滑法计算的坡印廷矢量的z分量示意图，图5(e)为利用光流法迭代了5次的坡印廷矢量的x分量示意图，图5(f)为利用光流法迭代了5次的坡印廷矢量的z分量示意图，图5(g)为利用时移法计算的坡印廷矢量的x分量示意图，图5(h)为利用时移法计算的坡印廷矢量的z分量示意图，图5(i)为利用本发明计算的坡印廷矢量的x分量示意图，图5(j)为利用本发明计算的坡印廷矢量的z分量示意图；
图6为坡印廷矢量计算的传播角与真实传播角之间的误差示意图；其中，图6(a)为原始坡印廷矢量计算的传播角与真实传播角之间的误差示意图，图6(b)为平滑法计算的传播角与真实传播角之间的误差示意图，图6(c)为光流法计算的传播角与真实传播角之间的误差示意图，图6(d)为时移法计算的传播角与真实传播角之间的误差示意图，图6(e)为本发明计算的传播角与真实传播角之间的误差示意图；
图7为Marmousi模型示意图；
图8为Marmousi模型在1.0秒时刻的地震波场示意图；
图9为对图8所示的地震波场利用原始坡印廷方程计算得到的坡印廷矢量示意图；其中，图9(a)为对图8所示的地震波场利用原始坡印廷方程计算得到的坡印廷矢量的x分量的示意图，图9(b)为对图8所示的地震波场利用原始坡印廷方程计算得到的坡印廷矢量的z分量的示意图；
图10为对图8所示的地震波场利用本发明得到的坡印廷矢量示意图；其中，图10(a)为对图8所示的地震波场利用本发明得到的坡印廷矢量的x分量的示意图，图10(b)为对图8所示的地震波场利用本发明得到的坡印廷矢量的z分量的示意图；
图11为本发明实施例提供的地震波场坡印廷矢量确定系统的结构图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
地震波场的传播矢量在地震勘探领域有着非常重要的作用。它是生成角度域共成像点道集的关键参数，在逆时偏移中被广泛地用于压制背向散射导致的低频噪音、照明度补偿和转换波极性反转的校正等。传播矢量可以在空间-时间域求取，也可以在频率-波数域求取。频率-波数域方法需要进行大量的傅里叶变换，计算效率低。由于波印廷矢量的方向与传播矢量一致，同时它只需要计算地震波场的空间和时间导数，其计算效率远高于频率-波数域的方法，因此得到广泛的应用。然而，坡印廷矢量会有计算不稳定的问题。
为了解决上述问题，实施例提供了一种地震波场坡印廷矢量确定方法，图2为本发明实施例提供的地震波场坡印廷矢量确定方法的流程图。参见图2，该方法，具体包括：
步骤101：对地震波场的时间导数的幅值进行调整，计算调整后的地震波场的时间导数的坡印廷矢量，得到第一坡印廷矢量；所述第一坡印廷矢量的幅值与第二坡印廷矢量的幅值的差值处于设定范围内；所述第二坡印廷矢量为所述地震波场的坡印廷矢量。
步骤102：对所述第二坡印廷矢量和所述第一坡印廷矢量进行运算，得到所述地震波场的最终坡印廷矢量。
作为一种可选的实施方式，所述步骤101中的设定范围可以为一个数量级，即该实施方式可以对地震波场的时间导数的幅值进行调整，使得调整后的地震波场的时间导数的坡印廷矢量(第一坡印廷矢量)的幅值与第二坡印廷矢量的幅值处于同一个数量级内，此处将二者调整到同一个数量级的目的是为了提高计算精度。
作为一种可选的实施方式，所述步骤101，具体包括：
1)计算所述地震波场的一阶时间导数。具体的，可以采用有限差分法计算所述地震波场的一阶时间导数。所述一阶时间导数的计算公式为：
Figure BDA0003033256130000061
其中，
Figure BDA0003033256130000062
为对ut求时间偏导数得到的一阶时间导数；ut为t时刻的地震波场；t为时间；Δt表示有限差分的时间间隔；ut+Δt为t+Δt时刻的地震波场；ut-Δt为t-Δt时刻的地震波场。
2)计算地震波场的坡印廷矢量，得到第二坡印廷矢量。具体为：
Figure BDA0003033256130000063
其中，Pu为第二坡印廷矢量；
Figure BDA0003033256130000064
为ut的梯度。
3)计算所述地震波场的一阶时间导数的幅值调整因子。具体为：
ut+Δt和ut-Δt可以表示为
ut+Δt＝utexp(-iωΔt)；
ut-Δt＝utexp(iωΔt)；
时间差分为
Figure BDA0003033256130000071
iut和ut幅值一样，只是相位相差90度。相差90度的好处是，两者坡印廷矢量的幅值互补，一个的极大值，对应着另一个的极小值。由上式可得：
Figure BDA0003033256130000072
由于ω＝2πf，ω为地震波场的角频率，f为地震波场的频率，则有
Figure BDA0003033256130000073
为了简便，令f＝fd，则
Figure BDA0003033256130000074
从而得到幅值调整因子的计算公式为：
Figure BDA0003033256130000075
其中，a为幅值调整因子；fd为地震波场的主频。
4)采用所述幅值调整因子调整所述地震波场的一阶时间导数的幅值，得到调整时间导数。具体的，所述调整时间导数的计算公式为：
Figure BDA0003033256130000076
其中，
Figure BDA0003033256130000081
为t时刻的调整时间导数。
5)计算所述调整时间导数的坡印廷矢量，得到所述第一坡印廷矢量。第一坡印廷矢量可以采用以下两种方法中的任意一种计算。
方法一：所述第一坡印廷矢量的计算公式为：
Figure BDA0003033256130000082
其中，
Figure BDA0003033256130000083
为第一坡印廷矢量；
Figure BDA0003033256130000084
为
Figure BDA0003033256130000085
的梯度。
方法二：所述第一坡印廷矢量的计算公式为：
Figure BDA0003033256130000086
上述方法二的计算过程避免了求取
Figure BDA0003033256130000087
的时间偏导数，减少了计算量，计算效率更高。
作为一种可选的实施方式，所述步骤102，具体包括：
将所述第二坡印廷矢量和所述第一坡印廷矢量相加，得到所述最终坡印廷矢量。具体计算公式为：
Figure BDA0003033256130000088
其中，P为最终坡印廷矢量。
上述实施例的地震波场坡印廷矢量确定方法，由于地震波场的波印廷矢量所与其时间导数的坡印廷矢量方向一致，虽然地震波场的坡印廷矢量的幅值比其时间导数的坡印廷矢量小很多，但是两者的变化趋势是互补的，如图3所示，地震波场的坡印廷矢量的极小值对应着时间导数坡印廷矢量的极大值。因此，将两者的幅值调整到一个设定范围(例如同一个数量级)内，然后运算(例如求和)，便可以避免零值或极小值的不稳定问题。
为了说明本发明的优势，将本发明与其他为了解决坡印廷矢量计算不稳定而提出的方法进行对比。
现有的解决坡印廷矢量计算不稳定的问题的方法有以下几种：
(1)空间平滑法。在获得坡印廷矢量后，利用高斯函数对其进行平滑。
(2)多重一维平滑法。对于二维和三维数据，直接利用高斯平滑函数，计算量比较大。为节省计算量，采用多次一维平滑的方法。比如，对于坡印廷矢量的x分量，首先在x方向上进行一维平滑，然后再在z方向上进行一维平滑。
(3)时间平滑法。首先在一个时间段内，计算出空间每个位置点上的坡印廷矢量，然后在时间上利用高斯函数进行平滑。
(4)时间积分法。首先在一个时间段内，计算出空间每个位置点上的坡印廷矢量，然后对其进行时间积分。将该积分值作为计算点处的波印廷矢量。
(5)时移法。首先在一个时间段内，计算出空间每个位置点上的坡印廷矢量，然后搜索该时间段内坡印廷矢量绝对值的最大值。将该最大值处的坡印廷矢量作为计算点处的坡印廷矢量。
(6)光流法。该方法从计算机视觉中发展而来，通过迭代的方式，逐步获得精确的传播方向。其迭代公式为：
Figure BDA0003033256130000091
其中，px、py和pz分别表示的是坡印廷矢量P在笛卡尔坐标下，x、y和z三个方向上的分量；上角标n表示的是迭代次数；
Figure BDA0003033256130000092
Figure BDA0003033256130000093
α是加权系数。
在以上这些方法中，方法(1)-方法(4)，虽然增加了坡印廷矢量的稳定性，但是降低了其计算精度。方法(5)和方法(6)，虽然计算精度和稳定性得到了保证，但是需要的进行大量的搜索运算(时移法)或迭代运算(光流法)，其计算效率不高。而本发明的利用地震波场及其一阶时间导数计算坡印廷矢量的方法，在确保计算精度和稳定性的同时，还有着较高计算效率。
为了进一步说明本发明的优势，下面采用两个具体实例将本发明提出的方法与现有的方法进行了对比。
本实例采用简单模型。图4为均匀介质中的0.4s时刻的地震波场示意图，介质的速度为2.0km/s，主频为30Hz。图5为利用原始的坡印廷矢量方程、空间光滑法、光流法、时移法和本发明所计算出的坡印廷矢量示意图，其中，原始的坡印廷矢量方程为：
Figure BDA0003033256130000101
其中，P表示的是坡印廷矢量，u表示的是地震波场。
从图4和图5中可以看出，原始的坡印廷矢量中白色条带和黑色条带中有明显的条纹，这些条纹上的值接近或等于零值。等于零值的点为不稳定点。而其他方法无条纹或条纹不明显。图6为利用图5所示的坡印廷矢量计算出的地震波场的传播角与真实传播角的误差。图6a中的幅值比较高的点，是坡印廷矢量计算不稳定的点。从图6a-图6e中可以看到平滑法、光流法、时移法以及本发明都提高了坡印廷矢量稳定性。本发明和时移法的误差最小，计算精度最高。
相对于现有技术，本发明有着较高的精度和计算效率。对于平滑算法来说，平滑起到了滤除高频保留低频的作用，虽然避免了计算不稳定的点，但是精度提高的有限。本方法计算精度较高是因为地震波场的坡印廷矢量和其一阶时间导数的坡印廷矢量所表示的传播方向是一致的，幅值调整因子不改变坡印廷矢量的方向。其次，地震波场的坡印廷矢量及其一阶时间导数的坡印廷矢量的幅值变化趋势互补，一个的极大值对应另一个的极小值，因此两者相加便可以抵消极小值带来的计算不稳定。
与时移法和光流法相比，本发明的计算精度与时移法相仿，高于光流法，但是本发明的计算效率高，因为时移法需要大量的搜索计算，光流法需要大量的迭代运算，而本发明只是需要多计算一次时间导数的坡印廷矢量。
本实例采用复杂模型。图7为Marmousi模型，图8为1.0s时的地震波场，图9为利用原始的坡印廷矢量方程计算得的坡印廷矢量，图10为利用原始的坡印廷矢量方程计算得的坡印廷矢量。对比图10和图9可以看出，图9中黑色和白色条带中的条纹非常明显，不稳定点多，图10中条纹不明显，稳定性有了很大的改善。
采用本发明提供的方法得到的坡印廷矢量在地震勘探中有着非常重要的作用，可应用多种不同的场景，例如，在逆时偏移中压制背向散射导致的低频噪音、照明度补偿和转换波极性反转的校正等。
逆时偏移中压制背向散射导致的低频噪音，利用本实施例提供的方法分别求取出震源波场的坡印廷矢量Psrc和检波点波场的坡印廷矢量Prec，其中下角标src表示震源波场，rec表示检波点波场，利用下式得到反射角
Figure BDA0003033256130000111
然后利用下式成像便可以压制低频噪音
Figure BDA0003033256130000112
其中I表示逆时偏移成像剖面，x表示空间位置，Tmax表示最大成像时间，n为大于1的正整数，一般取值大于等于3，效果比较好。
本发明还提供了一种地震波场坡印廷矢量确定系统，图11为本发明实施例提供的地震波场坡印廷矢量确定系统的结构图。参见图11，该系统包括：
幅值调整模块201，用于地震波场的时间导数的幅值进行调整，计算调整后的地震波场的时间导数的坡印廷矢量，得到第一坡印廷矢量；所述第一坡印廷矢量的幅值与第二坡印廷矢量的幅值的差值处于设定范围内；所述第二坡印廷矢量为所述地震波场的坡印廷矢量。
坡印廷矢量运算模块202，用于对所述第二坡印廷矢量和所述第一坡印廷矢量进行运算，得到所述地震波场的最终坡印廷矢量。
本实施例提供的地震波场坡印廷矢量确定系统，由于地震波场的波印廷矢量所与其时间导数的坡印廷矢量方向一致，虽然地震波场的坡印廷矢量的幅值比其时间导数的坡印廷矢量小很多，但是两者的变化趋势是互补的，地震波场的坡印廷矢量的极小值对应着时间导数坡印廷矢量的极大值。因此，将两者的幅值调整到一个设定范围内，然后运算，便可以避免零值或极小值的不稳定问题。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (9)

1.一种地震波场坡印廷矢量确定方法，其特征在于，包括：
对地震波场的时间导数的幅值进行调整，计算调整后的地震波场的时间导数的坡印廷矢量，得到第一坡印廷矢量；所述第一坡印廷矢量的幅值与第二坡印廷矢量的幅值的差值处于设定范围内；所述第二坡印廷矢量为所述地震波场的坡印廷矢量；
对所述第二坡印廷矢量和所述第一坡印廷矢量进行运算，得到所述地震波场的最终坡印廷矢量；
所述对地震波场的时间导数的幅值进行调整，计算调整后的地震波场的时间导数的坡印廷矢量，得到第一坡印廷矢量，具体包括：
计算所述地震波场的一阶时间导数；
计算所述地震波场的一阶时间导数的幅值调整因子；
采用所述幅值调整因子调整所述地震波场的一阶时间导数的幅值，得到调整时间导数；
计算所述调整时间导数的坡印廷矢量，得到所述第一坡印廷矢量。
2.根据权利要求1所述的一种地震波场坡印廷矢量确定方法，其特征在于，所述对所述第二坡印廷矢量和所述第一坡印廷矢量进行运算，得到所述地震波场的最终坡印廷矢量，具体包括：
将所述第二坡印廷矢量和所述第一坡印廷矢量相加，得到所述最终坡印廷矢量。
3.根据权利要求1所述的一种地震波场坡印廷矢量确定方法，其特征在于，所述计算所述地震波场的一阶时间导数，具体包括：
采用有限差分法计算所述地震波场的一阶时间导数。
4.根据权利要求3所述的一种地震波场坡印廷矢量确定方法，其特征在于，所述幅值调整因子的计算公式为：
Figure FDA0003549080100000021
其中，a为幅值调整因子；fd为地震波场的主频；Δt为有限差分的时间间隔。
5.根据权利要求1所述的一种地震波场坡印廷矢量确定方法，其特征在于，所述调整时间导数的计算公式为：
Figure FDA0003549080100000022
其中，
Figure FDA0003549080100000023
为t时刻的调整时间导数；a为幅值调整因子；t为时间；ut为t时刻的地震波场；
Figure FDA0003549080100000024
为对ut求时间偏导数得到的一阶时间导数。
6.根据权利要求1所述的一种地震波场坡印廷矢量确定方法，其特征在于，所述第一坡印廷矢量的计算公式为：
Figure FDA0003549080100000025
其中，
Figure FDA0003549080100000026
为第一坡印廷矢量；
Figure FDA0003549080100000027
为t时刻的调整时间导数；t为时间；
Figure FDA0003549080100000028
为
Figure FDA0003549080100000029
的梯度。
7.根据权利要求1所述的一种地震波场坡印廷矢量确定方法，其特征在于，所述第一坡印廷矢量的计算公式为：
Figure FDA0003549080100000031
其中，
Figure FDA0003549080100000032
为第一坡印廷矢量；
Figure FDA0003549080100000033
为t时刻的调整时间导数；t为时间；
Figure FDA0003549080100000034
为对
Figure FDA0003549080100000035
求梯度；ut为t时刻的地震波场；fd为地震波场的主频；Δt为有限差分的时间间隔。
8.根据权利要求3所述的一种地震波场坡印廷矢量确定方法，其特征在于，所述一阶时间导数的计算公式为：
Figure FDA0003549080100000036
其中，ut为t时刻的地震波场；t为时间；Δt为有限差分的时间间隔；ut+Δt为t+Δt时刻的地震波场；ut-Δt为t-Δt时刻的地震波场。
9.一种地震波场坡印廷矢量确定系统，其特征在于，包括：
幅值调整模块，用于对地震波场的时间导数的幅值进行调整，计算调整后的地震波场的时间导数的坡印廷矢量，得到第一坡印廷矢量；所述第一坡印廷矢量的幅值与第二坡印廷矢量的幅值的差值处于设定范围内；所述第二坡印廷矢量为所述地震波场的坡印廷矢量；
坡印廷矢量运算模块，用于对所述第二坡印廷矢量和所述第一坡印廷矢量进行运算，得到所述地震波场的最终坡印廷矢量；
所述对地震波场的时间导数的幅值进行调整，计算调整后的地震波场的时间导数的坡印廷矢量，得到第一坡印廷矢量，具体包括：
计算所述地震波场的一阶时间导数；
计算所述地震波场的一阶时间导数的幅值调整因子；
采用所述幅值调整因子调整所述地震波场的一阶时间导数的幅值，得到调整时间导数；
计算所述调整时间导数的坡印廷矢量，得到所述第一坡印廷矢量。
CN202110436362.5A 2021-04-22 2021-04-22 一种地震波场坡印廷矢量确定方法及系统 Active CN113126159B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202110436362.5A CN113126159B (zh) 2021-04-22 2021-04-22 一种地震波场坡印廷矢量确定方法及系统
US17/725,971 US20220350043A1 (en) 2021-04-22 2022-04-21 Method and system for stabilizing Poynting vector of seismic wavefield

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110436362.5A CN113126159B (zh) 2021-04-22 2021-04-22 一种地震波场坡印廷矢量确定方法及系统

Publications (2)

Publication Number Publication Date
CN113126159A CN113126159A (zh) 2021-07-16
CN113126159B true CN113126159B (zh) 2022-05-06

Family

ID=76779117

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110436362.5A Active CN113126159B (zh) 2021-04-22 2021-04-22 一种地震波场坡印廷矢量确定方法及系统

Country Status (2)

Country Link
US (1) US20220350043A1 (zh)
CN (1) CN113126159B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1664851A2 (en) * 2003-07-28 2006-06-07 Light Prescriptions Innovators LLC. Three-dimensional simultaneous multiple-surface method and free-form illumination-optics designed therefrom
CN104614765A (zh) * 2015-02-05 2015-05-13 西南石油大学 一种增强地震波激发照明的设计方法
CN106918840A (zh) * 2017-05-05 2017-07-04 中国石油化工股份有限公司 基于角度波场分步提取的逆时偏移角道集成像方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10514475B2 (en) * 2016-08-30 2019-12-24 Schlumberger Technology Corporation Post-critical reflection muting in seismic migration

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1664851A2 (en) * 2003-07-28 2006-06-07 Light Prescriptions Innovators LLC. Three-dimensional simultaneous multiple-surface method and free-form illumination-optics designed therefrom
CN104614765A (zh) * 2015-02-05 2015-05-13 西南石油大学 一种增强地震波激发照明的设计方法
CN106918840A (zh) * 2017-05-05 2017-07-04 中国石油化工股份有限公司 基于角度波场分步提取的逆时偏移角道集成像方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
高效逆时偏移角度道集生成方法研究;吴成梁 等;《石油物探》;20180531;第57卷(第3期);第404-418页 *

Also Published As

Publication number Publication date
US20220350043A1 (en) 2022-11-03
CN113126159A (zh) 2021-07-16

Similar Documents

Publication Publication Date Title
Hennenfent et al. Seismic denoising with nonuniformly sampled curvelets
RU2300123C2 (ru) Высокоразрешающее преобразование радона для обработки сейсмических данных
Yang et al. Optimal staggered-grid finite-difference schemes based on the minimax approximation method with the Remez algorithm
CN109100786B (zh) 深度域品质因子的确定方法和装置
Sun et al. Cross-correlation analysis and time delay estimation of a homologous micro-seismic signal based on the Hilbert–Huang transform
CN105277985A (zh) 一种基于图像处理的ovt域地震数据规则化方法
Liu et al. Seismic dip estimation based on the two-dimensional Hilbert transform and its application in random noise attenuation
Wang et al. Structure-oriented DTGV regularization for random noise attenuation in seismic data
Wason et al. Source separation via SVD-free rank minimization in the hierarchical semi-separable representation
CN111768349A (zh) 一种基于深度学习的espi图像降噪方法及系统
CN113126159B (zh) 一种地震波场坡印廷矢量确定方法及系统
Feng et al. Multigranularity feature fusion convolutional neural network for seismic data denoising
CN108828659B (zh) 基于傅里叶有限差分低秩分解的地震波场延拓方法及装置
Shao et al. Simultaneous inversion of Q and reflectivity using dictionary learning
Jahanjooy et al. Blind inversion of multidimensional seismic data using sequential Tikhonov and total variation regularizations
AlSalem et al. Embedded boundary methods for modeling 3D finite-difference Laplace-Fourier domain acoustic-wave equation with free-surface topography
CN109697704A (zh) 基于bm3d算法的自适应全变分espi图像降噪方法及系统
Champagnat et al. Deconvolution of sparse spike trains accounting for wavelet phase shifts and colored noise
CN111694056B (zh) 一种压制地震资料异常噪声的方法、存储介质及设备
Li Signal smoothing with time-space fractional order model
Kim et al. Least-squares reverse time migration using analytic-signal-based wavefield decomposition
CN107807391B (zh) 地震数据处理方法、装置、电子设备及计算机存储介质
Ma et al. Super-resolution stacking based on compressive sensing
Chen et al. Deblending based on the frequency-structure weighted threshold: The case of the slope-constrained f-k transform
CN110673207B (zh) 一种高频重建方法、装置及计算机存储介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant