CN115373022B - 一种基于振幅相位校正的弹性波场Helmholtz分解方法 - Google Patents

一种基于振幅相位校正的弹性波场Helmholtz分解方法 Download PDF

Info

Publication number
CN115373022B
CN115373022B CN202210024605.9A CN202210024605A CN115373022B CN 115373022 B CN115373022 B CN 115373022B CN 202210024605 A CN202210024605 A CN 202210024605A CN 115373022 B CN115373022 B CN 115373022B
Authority
CN
China
Prior art keywords
wave
longitudinal
vector
wave field
transverse
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
CN202210024605.9A
Other languages
English (en)
Other versions
CN115373022A (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.)
Yangtze University
Original Assignee
Yangtze University
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 Yangtze University filed Critical Yangtze University
Priority to CN202210024605.9A priority Critical patent/CN115373022B/zh
Publication of CN115373022A publication Critical patent/CN115373022A/zh
Application granted granted Critical
Publication of CN115373022B publication Critical patent/CN115373022B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • 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
    • 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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/53—Statics correction, e.g. weathering layer or transformation to a datum
    • 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/62—Physical property of subsurface

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

本发明公开了属于地震波信息处理技术领域，具体公开了一种基于振幅相位校正的弹性波场Helmholtz分解方法，包括：获取总矢量地震波场；对总矢量地震波场进行Helmholtz分解以得到标量纵波和矢量横波根据标量纵波和矢量横波与校正后的纵波P及横波S的一阶导数关系求解纵波P及横波S。本方法的相位校正完全在时间空间域内完成，克服了现有技术中波场分离方法在相位校正时操作较复杂，计算量大的技术问题。

Description

一种基于振幅相位校正的弹性波场Helmholtz分解方法
技术领域
本发明涉及地震波信息处理技术领域，具体涉及一种基于振幅相位校正的弹性波场Helmholtz分解方法。
背景技术
当今油气勘探正在由构造油气藏向岩性和隐蔽油气藏、由常规油气开发向页岩气和致密油等非常规油气开发转变，在此趋势下，多波地震勘探的优势意义重大。多波地震勘探技术能够为复杂储层油气藏勘探提供有利的解决方案，其优势在于有效的转换波成像。相比于常规纵波勘探可以获得更多的地下地质信息。多波勘探可以有效减小地球物理反演的多解性，不仅可以改善构造成像，在岩性描述、流体预测、裂缝检测、各向异性检测等油气藏勘探开发关键问题上，也能够提供更佳的解决方案。其中，弹性逆时偏移(elasticreverse-time migration，ERTM)是专门为多分量地震数据而发展的偏移成像技术。弹性逆时偏移直接以对多分量数据进行偏移处理，使用弹性波动方程进行波场延拓，能够保持波场的弹性特征和矢量特征，更为真实地模拟地震波场在地下介质中的传播过程；弹性逆时偏移可同时得到PP、PS、SP和SS等多种成像结果，为地震解释提供更可靠的依据。而充分体现这些特点和优势的关键在于提供有效的转换波成像。
波场分离是进行弹性逆时偏移的前提条件，最为传统和常用的波场分离方法为Helmholtz分解，该方法利用各向同性介质中“纵波无旋、横波无散”的性质，通过对弹性波场求取散度(▽·)得到标量纵波求取旋度(▽×)得到矢量横波/>但是这种方法分离出来的纵波和横波的相位和振幅都发生了改变，相位发生了π/2的改变，横纵波振幅比从c变成了cα/β，其中α和β分别为纵波和横波速度。现有技术中，Sun(2001,2011)给出了相位和纵横波振幅比校正的方法。对于振幅校正只需要在分离的横波上乘以纵横波速度比，但是对于相位校正，需要沿时间轴做Hilbert变换或者在频率域进行相位校正，操作较复杂，计算量大。
发明内容
本发明目的在于提供一种计算简单快速的基于时空域振幅相位校正的弹性波场Helmholtz分离方法，以克服现有技术中波场分离方法在相位校正时操作较复杂，计算量大的技术问题。
本发明所提供的基于振幅相位校正的弹性波场Helmholtz分解方法，包括：
获取总矢量地震波场；
对总矢量地震波场进行Helmholtz分解以得到标量纵波和矢量横波/>
根据以下关系求解经振幅相位校正以后的纵波P及横波S：
其中，α和β分别为纵波和横波速度。
进一步的，利用有限差分算法计算离散化的Helmholtz分解得到的标量纵波和矢量横波/>其中i、j、n分别是沿水平、垂直方向以及时间轴方向离散化后的位置点数；
通过以下的时间迭代计算获得离散化的经振幅相位校正以后的纵波和横波
其中，初始值和/>均设为0，Δt为所述有限差分算法在时间上的采样间隔。
进一步的，利用有限差分算法计算离散化的总矢量地震波场u、w分别为x和z方向的质点位移；
通过以下计算求得离散化的标量纵波和矢量横波/>
其中，Δx为所述有限差分算法在空间上的采样间隔。
进一步的，离散化的矢量地震波场依照如下的时间迭代计算得到:
本发明的原理和有益效果在于，已知在波数域中，Helmholtz分解的纵波和横波波场和/>的表达式为：
而校正以后的分离的纵波和横波波场P和S的表达式为：
式中，U＝(u,w)为总矢量波场，∧表示对应在波数域中的变量，k＝(kx,ky,kz)为波数向量，与传播方向一致，同时也代表着P波的振动方向，I＝(kx/|k|,ky/|k|,kz/|k|)为单位传播向量。所以/>和/>与P和S存在以下关系：
将P波频散关系|k|＝ω/α和S波频散关系|k|＝ω/β分别代入上式，并转到时间空间域，则可得到Helmholtz分解的纵波和横波波场和/>与校正以后的分离的纵波和横波波场P和S存在以下关系：
可以看出，根据本发明，分解后的振幅校正即是将Helmholtz分解的纵横波波场和/>分别乘以各自的纵波和横波速，而对于相位校正，由于校正后波场的一阶时间导数等于未经相位校正的波场，因此相位校正过程可以通过时间迭代，利用前时刻的波场值计算得到。于是，振幅和相位校正完全可以在时间空间域完成，计算简单方便，另外，利用有限差分算法，还可随着总地震矢量波场的正演模拟计算的时间迭代过程同步的进行校正，与现有技术中的分解方法及相应的振幅和相位校正方法相比，本发明计算过程简单方便，成本更低。
附图说明
图1为本发明实施例中的基于波场模型的模拟实验的具体实施流程图。
图2为本发明实施例的模拟实验中1.2s时刻的波场快照，其中子图a和b分别是原始矢量波场水平(u)和垂直方向位移(w)，子图c和d分别是传统Helmholtz分解得到的纵波波场和横波波场/>子图e和f分别是振幅和相位校正后的纵波波场(P)和横波波场(S)。
图3为本发明实施例的模拟实验中，(534m,400m)处检波点接收到的波形图(子图a)和能量振幅图(子图b)。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本实施例中，利用基于波场模型的模拟实验示例性的说明本发明的技术方案及技术效果。
本实施例中的基于振幅相位校正的弹性波场Helmholtz分解方法包括，
获取总地震矢量波场；
对总地震矢量波场进行Helmholtz分解以得到标量纵波和矢量横波/>
根据以下关系求解经振幅相位校正以后的纵波P及横波S：
其中，α和β分别为纵波和横波速度。
于是，本实施例中首先基于二维弹性波动方程在时间空间域进行有限差分正演模拟，计算出矢量地震波场。
以二阶位移方程为例，有：
式中，u、w分别为x和z方向的质点位移，α、β分别为地下介质P波和S波速度。
于是，根据有限差分数值计算原理进行离散化，得到u、w分别的时间迭代公式：
式中，Δx和Δt分别是有限差分算法在空间和时间上的采样间隔，i、j、n分别是沿水平、垂直方向以及时间轴方向离散化后的位置点数。在模拟计算中，各采样点的初始位移可设定为0，震源位置(nsx,nxz)则加载震源时间函数fs，于是根据上式，矢量地震波场中各采样点在各采样时刻的波场值皆可计算得到。
而在实际工作中，震源位置和震源时间函数可以通过对地震的观测或通过观测数据推算得到，进而模拟出总地震矢量波场，但不排除直接测量或推算总地震矢量波场的情况，无论何种方式获得的总地震矢量波场(包括其离散表达)，皆不影响本发明中的分解方法的实施，依旧属于本发明的保护范围。
然后，基于Helmholtz分解定理进行散度和旋度运算，分别计算出纵波和横波波场和/>此波场是未经校正的，与原纵横波场存在振幅和相位上的差异，无明确的物理意义。对矢量波场进行散度和旋度运算如下：
离散化以后的表达式为：
最后，进行振幅和相位校正，我们已知在波数域中，Helmholtz分解的纵波和横波波场和/>的表达式为：
而校正以后的分离的纵波和横波波场P和S的表达式为：
式中，U＝(u,w)为总矢量波场，∧表示对应在波数域中的变量，k＝(kx,ky,kz)为波数向量，与传播方向一致，同时也代表着P波的振动方向，I＝(kx/|k|,ky/|k|,kz/|k|)为单位传播向量。所以/>和/>与P和S存在以下关系：
将P波频散关系|k|＝ω/α和S波频散关系|k|＝ω/β分别代入上式，并转到时间空间域，则可得到Helmholtz分解的纵波和横波波场和/>与校正以后的分离的纵波和横波波场P和S存在以下关系：
可以看出，振幅校正即是将Helmholtz分解的纵横波波场和/>分别乘以各自的纵波和横波速。对于相位校正，由于校正后波场的一阶时间导数等于未经相位校正的波场，因此相位校正过程需要进行时间迭代，利用前时刻波场值计算得到，可随着正演模拟计算总波场的时间迭代过程进行校正。离散化上式以后的表达式为：
其中，初始值和/>均设为0，P和S即为最后计算得到的经振幅和相位校正以后的分离的纵波和横波波场。
本实施例中，模拟实验采用以下相关参数设定：
1.模型参数：模型大小为1600×1600m2，均匀各向同性介质，其中P波和S波速度分别是α＝3000m/s、β＝1500m/s。
2.震源参数：震源类型为P波和S波混合震源，震源时间函数fs为25Hz的Ricker子波，震源位于模型中心(800m,800m)处。
3.有限差分数值算法参数：空间采样间隔Δx＝2m，时间采样间隔Δt＝0.5ms，网格点数800×800，总模拟时长为3s，时间迭代次数为nt＝6000。边界条件采用PML吸收边界，厚度为50层。
模拟实验的具体实施流程基本如图1所示：
s1.设定有限差分算法相关参数，包括空间采样间隔Δx、时间采样间隔Δt、水平和垂直网格点数nx和nz，模拟时间长度，读取速度模型包括模型大小、纵波速度α和横波速度β，并将模型参数网格离散化，读取震源参数包括震源位置、震源类型、震源时间函数，设定输出参数包括检波点个数和位置、波场快照时间等。
s2.进入时间循环，直至达到有限差分时间循环总次数nt。
s3.加载震源，根据震源类型在震源位置(nsx,nxz)处加载震源项fs，本实施例中，设定震源为爆炸源，有：
unsx,nxz＝unsx,nxz+fs
wnsx,nxz＝wnsx,nxz+fs。
s4.波场迭代：根据波场有限差分时间迭代公式，利用前一时刻波场值当前时刻波场值/>纵波速度α和横波速度β以及有限差分空间和时间采样间隔Δx、Δt，计算得到下一时刻的波场值/>同时为防止模型四周出现强反射，需要加载PML吸收边界。
s5.利用散度和旋度算子计算得到Helmholtz分解波场，对总矢量波场求取散度(▽·)得到纵波波场/>求取旋度(▽×)得到横波波场/>
s6.计算得到经振幅和相位校正的Helmholtz分解波场。利用当前时刻的Helmholtz分解的纵波波场和横波波场/>前时刻经振幅和相位校正的Helmholtz分解波场/> 纵波速度α和横波速度β和有限差分时间采样间隔Δt，计算得到下一时刻经振幅和相位校正的Helmholtz分解波场/>
s7.存储波场快照和地震记录。根据设定的时间存储相应的波场值，包括总矢量波场Helmholtz分解波场/>和横波波场/>经振幅和相位校正的Helmholtz分解波场/>另外根据检波点位置储存每一时刻的波场值。
s8.判断是否达到有限差分时间循环总次数nt，未达到则返回第2步，否则时间循环终止。
s9.输出波场快照图和地震记录图。
模拟结果如图2和图3所示，其中，图2是1.2s时刻的波场快照，a和b子图分别是原始矢量波场水平(u)和垂直方向位移(w)，c和d分别是传统Helmholtz分解得到的纵波波场和横波波场/>e和f分别是振幅和相位校正后的纵波波场(P)和横波波场(S)。图3是(534m,400m)处检波点接收到的波形图(子图a)和能量振幅图(子图b)。可以看出，传统Helmholtz分解得到的波场/>和/>(图2c和2d，图3中以圆点标记的实线和以星号标记的虚线)振幅较小，纵横波振幅比与原始波场(图2a和2b，图3中实线和虚线)相比发生变化，相位也改变了π/2，总之，振幅相位的变化导致传统Helmholtz分解得到的波场无实际物理意义。经振幅相位校正的波场P和S(图2e和2f，图3中以方块标记的实线和以菱形标记的虚线)与原始波场相比，相位相同、保持了纵横波的振幅比，且如图3的b子图中所示，能量也与原始波场一致。
本发明中未涉及部分均与现有技术相同或可采用现有技术加以实现。尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (3)

1.一种基于振幅相位校正的弹性波场Helmholtz分解方法，包括：
获取总矢量地震波场；
对总矢量地震波场进行Helmholtz分解以得到标量纵波和矢量横波/>
其特征在于，根据以下关系求解经振幅相位校正以后的纵波P及横波S：
其中，α和β分别为纵波和横波速度；
利用有限差分算法计算离散化的Helmholtz分解得到的标量纵波和矢量横波/>其中i、j、n分别是沿水平、垂直方向以及时间轴方向离散化后的位置点数；
通过以下的时间迭代计算获得离散化的经振幅相位校正以后的纵波和横波/>
其中，初始值和/>均设为0，Δt为所述有限差分算法在时间上的采样间隔。
2.根据权利要求1所述的方法，其特征在于，利用有限差分算法计算离散化的总矢量地震波场u、w分别为x和z方向的质点位移；
通过以下计算求得离散化的标量纵波和矢量横波/>
其中，Δx为所述有限差分算法在空间上的采样间隔。
3.根据权利要求2所述的方法，其特征在于，离散化的矢量地震波场依照如下的时间迭代计算得到:
CN202210024605.9A 2022-01-11 2022-01-11 一种基于振幅相位校正的弹性波场Helmholtz分解方法 Active CN115373022B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210024605.9A CN115373022B (zh) 2022-01-11 2022-01-11 一种基于振幅相位校正的弹性波场Helmholtz分解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210024605.9A CN115373022B (zh) 2022-01-11 2022-01-11 一种基于振幅相位校正的弹性波场Helmholtz分解方法

Publications (2)

Publication Number Publication Date
CN115373022A CN115373022A (zh) 2022-11-22
CN115373022B true CN115373022B (zh) 2024-04-19

Family

ID=84060816

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210024605.9A Active CN115373022B (zh) 2022-01-11 2022-01-11 一种基于振幅相位校正的弹性波场Helmholtz分解方法

Country Status (1)

Country Link
CN (1) CN115373022B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117233838B (zh) * 2023-09-20 2024-04-05 长江大学 一种二维vti介质中的弹性准纵横波场分离和逆时偏移成像方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4766574A (en) * 1987-03-31 1988-08-23 Amoco Corporation Method for depth imaging multicomponent seismic data
CN105242305A (zh) * 2015-09-06 2016-01-13 中国科学院地质与地球物理研究所 一种纵波和横波的分离方法及系统
CN112904426A (zh) * 2021-03-27 2021-06-04 中国石油大学(华东) 一种解耦弹性波逆时偏移方法、系统及应用

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4766574A (en) * 1987-03-31 1988-08-23 Amoco Corporation Method for depth imaging multicomponent seismic data
CN105242305A (zh) * 2015-09-06 2016-01-13 中国科学院地质与地球物理研究所 一种纵波和横波的分离方法及系统
CN112904426A (zh) * 2021-03-27 2021-06-04 中国石油大学(华东) 一种解耦弹性波逆时偏移方法、系统及应用

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Scalar imaging condition for elastic reverse time migration;Yuting Duan 等;GEOPHYSICS;20150831;第80卷(第4期);S127–S136 *
基于散度和旋度纵横波分离方法的改进;李志远 等;地球物理学报;20131231;第56卷(第6期);2012-2022 *

Also Published As

Publication number Publication date
CN115373022A (zh) 2022-11-22

Similar Documents

Publication Publication Date Title
RU2693495C1 (ru) Полная инверсия волнового поля с компенсацией показателя качества
EP3063562B1 (en) Methods of subsurface exploration, computer program product and computer-readable storage medium
CN111221037B (zh) 解耦弹性逆时偏移成像方法和装置
US11041971B2 (en) Full wavefield inversion with an image-gather-flatness constraint
CN113740901B (zh) 基于复杂起伏地表的陆上地震数据全波形反演方法及装置
CN111025387B (zh) 一种页岩储层的叠前地震多参数反演方法
Wang et al. 2D frequency-domain elastic full-waveform inversion using the block-diagonal pseudo-Hessian approximation
CN111766628A (zh) 一种预条件的时间域弹性介质多参数全波形反演方法
CN101201409A (zh) 一种地震数据变相位校正方法
Zhang et al. Elastic full waveform inversion with source-independent crosstalk-free source-encoding algorithm
Wang et al. High-resolution wave-equation AVA imaging: Algorithm and tests with a data set from the Western Canadian Sedimentary Basin
CN115373022B (zh) 一种基于振幅相位校正的弹性波场Helmholtz分解方法
Vigh et al. Elastic full-waveform inversion using OBN data acquisition
CN107229066A (zh) 基于地面地震构造约束的vsp数据全波形反演建模方法
CN114814944B (zh) 基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法
Zhong et al. Elastic reverse time migration method in vertical transversely isotropic media including surface topography
CN105319594B (zh) 一种基于最小二乘参数反演的傅里叶域地震数据重构方法
Turco et al. Geostatistical interpolation of non-stationary seismic data
Qu et al. An elastic full-waveform inversion based on wave-mode separation
CN113866823A (zh) 一种粘声各向异性介质中的正演成像方法
Cho et al. Accelerating 2D frequency-domain full-waveform inversion via fast wave modeling using a model reduction technique
Biondi et al. Preconditioned elastic full-waveform inversion with approximated Hessian
Roberts et al. Investigation into the use of 2D elastic waveform inversion from look‐ahead walk‐away VSP surveys
Plessix et al. Frequency-domain finite-difference migration with only few frequencies?
Jaimes-Osorio et al. Inversion comparison using an elastic local solver to recover elastic parameters

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