CN113126151B - 一种基于纯波延拓方程的弹性反射波旅行时反演方法 - Google Patents

一种基于纯波延拓方程的弹性反射波旅行时反演方法 Download PDF

Info

Publication number
CN113126151B
CN113126151B CN202110260212.3A CN202110260212A CN113126151B CN 113126151 B CN113126151 B CN 113126151B CN 202110260212 A CN202110260212 A CN 202110260212A CN 113126151 B CN113126151 B CN 113126151B
Authority
CN
China
Prior art keywords
wave
pure
background
equation
velocity
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
CN202110260212.3A
Other languages
English (en)
Other versions
CN113126151A (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 Oceanographic Instrumentation Shandong Academy of Sciences
Original Assignee
Institute of Oceanographic Instrumentation Shandong Academy of Sciences
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 Oceanographic Instrumentation Shandong Academy of Sciences filed Critical Institute of Oceanographic Instrumentation Shandong Academy of Sciences
Priority to CN202110260212.3A priority Critical patent/CN113126151B/zh
Publication of CN113126151A publication Critical patent/CN113126151A/zh
Application granted granted Critical
Publication of CN113126151B publication Critical patent/CN113126151B/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/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/67—Wave propagation modeling
    • G01V2210/675—Wave equation; Green's functions

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)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种基于纯波延拓方程的弹性反射波旅行时反演方法，包括如下步骤：(1)将近似弹性波解耦延拓方程省略交叉项得到纯波延拓方程，基于纯波延拓方程计算背景波场；(2)利用得到的背景波场计算逆时偏移结果；(3)基于纯波延拓方程和逆时偏移结果求解扰动纯波延拓方程，再利用扰动纯波延拓方程计算扰动波场；(4)基于步骤(1)得到的背景波场和步骤(3)得到的扰动波场，分别计算纵横波梯度；(5)根据纵横波梯度计算步长大小，进行纵横波速度更新，当迭代值收敛到全局最小值时，输出更新的纵横波速度即为最终的反演结果。本发明所公开的方法可以提高反演计算效率和反演精度，为弹性波地震勘探提供技术支持。

Description

一种基于纯波延拓方程的弹性反射波旅行时反演方法
技术领域
本发明属于弹性波地震勘探领域，特别涉及一种基于纯波延拓方程的弹性反射波旅行时反演方法。
背景技术
弹性反射波旅行时反演以弹性波动理论为基础，匹配模拟地震反射数据与观测地震反射数据，沿着弹性波的传播路径恢复中低波数的速度分量，可以为弹性全波形反演、弹性逆时偏移等提供准确的初始速度模型，是目前弹性波地震勘探领域的前沿技术之一。但是，计算弹性反射波旅行时反演的梯度敏感核函数需要引入偏移/反偏移过程，在一次迭代过程中求解6次弹性波波动方程，巨大的计算量和海量的存储需求制约了其走向实际生产应用的步伐。
另外，由于弹性波场中多波模式的同时存在，弹性反射波旅行时反演存在严重的高波数串扰，这主要体现在横波梯度的敏感核函数中。弹性波解耦延拓可以有效的减弱弹性参数反演的耦合效应，但是反向延拓的横波应力波场解耦依然会引入高波数噪音，降低横波速度反演的精度。
发明内容
为解决上述技术问题，本发明提供了一种基于纯波延拓方程的弹性反射波旅行时反演方法，以达到提高反演计算效率和反演精度的目的。
为达到上述目的，本发明的技术方案如下：
一种基于纯波延拓方程的弹性反射波旅行时反演方法，包括如下步骤：
(1)将近似弹性波解耦延拓方程省略交叉项得到纯波延拓方程，基于纯波延拓方程计算背景波场；
(2)利用得到的背景波场计算逆时偏移结果；
(3)基于纯波延拓方程和逆时偏移结果求解扰动纯波延拓方程，再利用扰动纯波延拓方程计算扰动波场；
(4)基于步骤(1)得到的背景波场和步骤(3)得到的扰动波场，分别计算纵横波梯度；
(5)根据纵横波梯度计算步长大小，进行纵横波速度更新，当迭代值收敛到全局最小值时，输出更新的纵横波速度即为最终的反演结果。
上述方案中，步骤(1)具体如下：纯波延拓方程包括纯P波延拓方程和纯S波延拓方程，分别基于纯P波延拓方程计算正向传播的P波振动速度背景波场和反向传播的P波振动速度背景波场，基于纯S波延拓方程计算反向传播的S波振动速度背景波场，基于纯P波延拓方程计算反向传播的P波应力背景波场，基于纯S波延拓方程计算反向传播的S波应力背景波场，背景波场的数值模拟计算采用高阶交错网格有限差分实施。
进一步的技术方案中，纯P波延拓方程为：
Figure BDA0002969635470000021
其中，下标“0”表示背景，x和z分别表示空间水平方向和垂直方向，t表示时间，vPx0和vPz0分别表示P波振动速度背景波场的水平分量和垂直分量，τPxx0和τPzz0分别表示xx和zz方向的P波应力背景波场，ρ0和VP0分别表示背景密度和背景纵波速度；
纯S波延拓方程为：
Figure BDA0002969635470000022
其中，vSx0和vSz0分别表示S波振动速度背景波场的水平分量和垂直分量，τSxz0和τSzx0分别表示xz和zx方向的S波应力背景波场，VS0表示背景横波速度。
上述方案中，步骤(2)具体如下：分别用计算得到的正向传播的P波振动速度背景波场和反向传播的P波振动速度背景波场互相关计算PP波逆时偏移结果，用正向传播的P波振动速度背景波场和反向传播的S波振动速度背景波场计算PS波逆时偏移结果。
上述方案中，步骤(3)具体如下：以背景波场和速度扰动的相互作用为震源构建Born正演方程，即扰动的纯波延拓方程，扰动的纯波延拓方程包括扰动的纯P波延拓方程和扰动的纯S波延拓方程；其中，速度扰动用逆时偏移结果代替；
基于扰动的纯P波延拓方程计算正向传播的P波振动速度扰动波场和反向传播的P波应力扰动波场，基于扰动的纯S波延拓方程计算正向传播的S波振动速度扰动波场。
上述方案中，扰动的纯P波延拓方程为：
Figure BDA0002969635470000031
其中，符号“δ”表示扰动，ρ0表示背景密度，δVP表示纵波速度扰动，δVS表示横波速度扰动，vPx0和vPz0分别表示P波振动速度背景波场的水平分量和垂直分量，VP0表示背景纵波速度，δvPx和δvPz分别表示P波振动速度扰动波场的水平分量和垂直分量，δτPxx和δτPzz分别表示xx和zz方向的P波应力扰动波场；在实施过程中，方程(3)中的δVP和δVS用PP波逆时偏移结果表示；
扰动的纯S波延拓方程为：
Figure BDA0002969635470000032
在实施过程中，方程(4)中的δVS用PS波逆时偏移结果表示，δvSx和δvSz分别表示S波振动速度扰动波场的水平分量和垂直分量，δτSxz和δτSzx分别表示xz和zx方向的S波应力扰动波场，VS0表示背景横波速度。
上述方案中，步骤(4)具体方法如下：读取步骤(1)计算的背景波场和步骤(3)计算的扰动波场，分别求取纵波速度梯度
Figure BDA0002969635470000033
和横波速度梯度
Figure BDA0002969635470000034
整个计算过程中保持背景密度ρ0为常数；
Figure BDA0002969635470000041
其中，χP和χS分别表示匹配P波数据和匹配S波数据的目标函数，VP0表示背景纵波速度，vPx0和vPz0分别表示正向传播P波振动速度背景波场的水平分量和垂直分量，δvPx和δvPz分别表示正向传播P波振动速度扰动波场的水平分量和垂直分量，
Figure BDA0002969635470000042
和
Figure BDA0002969635470000043
分别表示反向传播P波应力扰动波场的xx分量和zz分量，
Figure BDA0002969635470000044
和
Figure BDA0002969635470000045
分别表示反向传播P波应力背景波场的xx分量和zz分量，VS0表示背景横波速度，δvSx和δvSz分别表示正向传播S波振动速度扰动波场的水平分量和垂直分量，
Figure BDA0002969635470000046
和
Figure BDA0002969635470000047
分别表示反向传播S波应力背景波场的xz分量和zx分量。
上述方案中，步骤(5)中，保证每次迭代的纵波速度更新量在20m/s-100m/s之间，横波速度更新量在10m/s-60m/s之间。
通过上述技术方案，本发明提供的基于纯波延拓方程的弹性反射波旅行时反演方法具有如下有益效果：
本发明利用近似弹性波解耦延拓方程得到相应的纯波延拓方程，在弹性反射波旅行时反演的实施流程中，基于纯波延拓方程计算所有的波场，可以有效的降低弹性反射波旅行时反演的计算量和存储量，提高3倍以上的计算效率。另外，弹性波解耦延拓计算的横波应力波场中会存在纵波分量，会给横波梯度计算带来高波数噪音，采用纯波延拓方程可以有效的消除横波应力波场中的纵波分量，提高反演精度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。
图1是基于纯波延拓方程的弹性反射波旅行时反演方法实施的流程示意图；
图2a-2d是单层速度模型，其中图2a为真实的纵波速度模型；图2b为真实的横波速度模型；图2c为初始的纵波速度模型；图2d为初始的横波速度模型
图3a-3f是基于弹性波解耦延拓方程和纯波延拓方程计算的反向横波应力波场的对比，其中图3a，3b和3c分别为基于弹性波解耦延拓方程计算的反向横波应力波场的xx，zz和xz分量；图3d，3e和3f分别为基于纯波延拓方程计算的反向横波应力波场的xx，zz和xz分量；
图4a-4b是基于弹性波解耦延拓方程和纯波延拓方程计算的横波梯度对比，其中图4a为基于弹性波解耦延拓方程计算的横波梯度；图4b为基于纯波延拓方程计算的横波梯度；
图5a-5d是Sigsbee2A速度模型，其中图5a为真实的纵波速度模型；图5b为真实的横波速度模型；图5c为初始的纵波速度模型；图5d为初始的横波速度模型；
图6a-6b是初始速度的成像结果，其中图6a为PP波成像结果；图6b为PS波成像结果；
图7a-7b是最终的速度反演结果，其中图7a为纵波速度；图7b为横波速度；
图8a-8f是真实速度、初始速度以及反演速度的伪井曲线对比，其中图8a，8b和8c分别为横向距离1200m，1500m和2200m处的纵波速度伪井曲线对比；图8d，8e和8f分别为横向距离1200m，1500m和2200m处的横波速度伪井曲线对比；
图9a-9b是反演速度的成像结果，其中图9a为PP波成像结果；图9b为PS波成像结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。
本发明提供了一种基于纯波延拓方程的弹性反射波旅行时反演方法，如图1所示，通过模型测试说明具体的技术方案：
第一步，输入数据
输入各炮点的地震子波数据；观测地震数据，包括纵横波矢量分离的P波地震数据、S波地震数据；纵横波初始速度模型，初始速度的速度值一般是从浅层至深层逐渐增大。
第二步，基于纯波延拓方程计算背景波场
利用近似弹性波解耦延拓方程，省略交叉项，得到纯波延拓方程，包括纯P波延拓方程和纯S波延拓方程。
其中纯P波延拓方程为：
Figure BDA0002969635470000051
其中，下标“0”表示背景，x和z分别表示空间水平方向和垂直方向，t表示时间，vPx0和vPz0分别表示P波振动速度背景波场的水平分量和垂直分量，τPxx0和τPzz0分别表示xx和zz方向的P波应力背景波场，ρ0和VP0分别表示背景密度和背景纵波速度；
纯S波延拓方程为：
Figure BDA0002969635470000061
其中，vSx0和vSz0分别表示S波振动速度背景波场的水平分量和垂直分量，τSxz0和τSzx0分别表示xz和zx方向的S波应力背景波场，VS0表示背景横波速度。
分别基于纯P波延拓方程(方程(1))计算正向传播的P波振动速度背景波场和反向传播的P波振动速度背景波场，基于纯S波延拓方程(方程(2))计算反向传播的S波振动速度背景波场，波场的数值模拟计算采用高阶交错网格有限差分实施。
建立单层速度模型(图2a至图2d)进行数值模拟测试。图3a-图3f展示是基于弹性波解耦延拓方程和纯波延拓方程计算的反向横波应力波场的对比，可以看出基于弹性波解耦延拓方程计算的反向横波应力波场中存在纵波成分的干扰，而基于纯波延拓方程计算的横波应力波场则可以消除纵波分量的干扰。从图4a可以看出，这种纵波成分的干扰会在横波梯度计算中产生高波数噪音(图4a中箭头处)，从图4b可以看出，基于纯波延拓方程计算的横波梯度则可以避免此类高波数噪音的影响。
第三步，计算逆时偏移结果
利用第二步得到的背景波场，分别用正向传播的P波振动速度波场和反向传播的P波振动速度波场互相关计算PP波逆时偏移结果，用正向传播的P波振动速度波场和反向传播的S波振动速度波场互相关计算PS波逆时偏移结果。
图5-图9以Sigsbee2A模型进行基于弹性反射波旅行时反演的数值测试。其中，图5a至图5d展示了Sigsbee2A模型的初始速度和真实速度。图6a和图6b展示了初始速度的逆时偏移结果，其中图6a为PP波逆时偏移结果，图6b为PS波逆时偏移结果，可以看出由于初始速度中的低波数速度分量不准确，PP波和PS波的逆时偏移结果没有得到有效的归位，特别是中深层绕射体的位置出现了“画弧”现象。
第四步，基于纯波延拓方程和逆时偏移结果计算扰动波场
以背景波场和速度扰动的相互作用为震源构建Born正演方程，即扰动的纯波延拓方程，其中速度扰动用第三步的逆时偏移结果代替。
扰动的纯P波延拓方程为：
Figure BDA0002969635470000071
其中，符号“δ”表示扰动，ρ0表示背景密度，δVP表示纵波速度扰动，δVS表示横波速度扰动，vPx0和vPz0分别表示P波振动速度背景波场的水平分量和垂直分量，VP0表示背景纵波速度，δvPx和δvPz分别表示P波振动速度扰动波场的水平分量和垂直分量，δτPxx和δτPzz分别表示xx和zz方向的P波应力扰动波场。在实施过程中，方程(3)中的δVP和δVS用PP波逆时偏移结果(图6a)表示。
扰动的纯S波延拓方程为：
Figure BDA0002969635470000072
其中，δvSx和δvSz分别表示S波振动速度扰动波场的水平分量和垂直分量，δτSxz和δτSzx分别表示xz和zx方向的S波应力扰动波场，VS0表示背景横波速度。在实施过程中，方程(4)中的δVS用PS波逆时偏移结果(图6b)表示。
分别基于纯P波延拓方程(方程(1))计算反向传播的P波应力背景波场，基于扰动的纯P波延拓方程(方程(3))计算正向传播的P波振动速度扰动波场和反向传播的P波应力扰动波场；基于纯S波延拓方程(方程(2))计算反向传播的S波应力背景波场，基于扰动的纯S波延拓方程(方程(4))计算正向传播的S波振动速度扰动波场。
表1和表2中对比了基于弹性波解耦延拓方程和纯波延拓方程计算的波场时间的对比，输入的速度为图2所示。
表1纵波梯度的计算时间对比
Figure BDA0002969635470000081
其中，uP0为正向传播的P波背景波场；δuP为正向传播的P波扰动波场；
Figure BDA00029696354700000812
为反向传播的P波背景波场；
Figure BDA00029696354700000813
为反向传播的P波扰动波场。
表2横波梯度的计算时间对比
Figure BDA0002969635470000082
其中，uP0为正向传播的P波背景波场；δuS为正向传播的S波扰动波场；
Figure BDA00029696354700000814
为反向传播的S波背景波场。
对比可以看出，基于纯波延拓方程计算波场的时间大约是基于弹性波解耦延拓方程计算波场时间的1/3，说明本发明可以有效的提弹性反射波旅行时反演的计算效率。
第五步，纵横波梯度计算
读取第二步和第四步计算的波场，根据方程(5)分别计算纵波速度梯度
Figure BDA0002969635470000083
和横波速度梯度
Figure BDA0002969635470000084
整个计算过程中保持背景密度ρ0为常数。
Figure BDA0002969635470000085
其中，χP和χS分别表示匹配P波数据和匹配S波数据的目标函数，VP0表示背景纵波速度，vPx0和vPz0分别表示正向传播P波振动速度背景波场的水平分量和垂直分量，δvPx和δvPz分别表示正向传播P波振动速度扰动波场的水平分量和垂直分量，
Figure BDA0002969635470000086
和
Figure BDA0002969635470000087
分别表示反向传播P波应力扰动波场的xx分量和zz分量，
Figure BDA0002969635470000088
和
Figure BDA0002969635470000089
分别表示反向传播P波应力背景波场的xx分量和zz分量，VS0表示背景横波速度，δvSx和δvSz分别表示正向传播S波振动速度扰动波场的水平分量和垂直分量，
Figure BDA00029696354700000810
和
Figure BDA00029696354700000811
分别表示反向传播S波应力背景波场的xz分量和zx分量。
其中，目标函数χP和χS具体如下：在地表位置处记录正向传播的振动速度扰动波场数据作为模拟地震数据，分别匹配模拟的P波数据与观测的P波地震数据建立目标函数见方程(6)，匹配模拟的S波数据与观测的S波地震数据建立目标函数见方程(8)。
Figure BDA0002969635470000091
其中，χP表示以P波数据建立的目标函数，||·||2表示L2范数，ΔτP表示模拟的P波数据与观测的P波地震数据之间的旅行时差。ΔτP具体计算采用方程(7)：
Figure BDA0002969635470000092
其中，
Figure BDA0002969635470000093
和
Figure BDA0002969635470000094
分别表示模拟的P波数据与观测的P波地震数据。当互相关值CP(τ)取最大时，τ与ΔτP相等。
Figure BDA0002969635470000095
其中，χS表示以S波数据建立的目标函数，ΔτS表示模拟的S波数据与观测的S波地震数据之间的旅行时差。ΔτS具体计算采用方程(9)：
Figure BDA0002969635470000096
其中，
Figure BDA0002969635470000097
和
Figure BDA0002969635470000098
分别表示模拟的S波数据与观测的S波地震数据。当互相关值CS(τ)取最大时，τ与ΔτS相等。
至此，通过第一步至第五步的计算得到了单炮的纵横波梯度，随后累加所有炮的梯度，完成纵横波梯度的计算。
第六步，纵横波速度更新
根据纵横波梯度计算步长大小，保证每次迭代的纵波速度更新量在20m/s-100m/s之间，横波速度更新量在10m/s-60m/s之间。
第七步，判断收敛性
计算本次迭代||ΔτP||2+||ΔτS||2的值，与上一次迭代的值进行比较，倘若数值减小，则继续第二步的计算。
第八步，得到纵横波速度的反演结果
当||ΔτP||2+||ΔτS||2的值连续5次处于非下降状态，则认为反演方法已经收敛到全局最小值，此时，输出纵横波速度数据即为最终的反演结果。
图7展示的是收敛后的纵横波速度的反演结果，可以看出基于纯波延拓方程的弹性反射波旅行时反演主要更新的是速度模型中的低波数分量，分别抽取不同横向距离处的真实速度、初始速度以及反演速度的伪井曲线对比于图8，经过更新后的速度伪井曲线与真实速度的更加贴合，说明本发明可以有效的恢复大尺度的背景速度分量。以反演速度进行逆时偏移成像结果展示于图9，对比图6中的结果可以看出，在速度中的低波数分量准确恢复之后，PP波和PS波的逆时偏移结果皆有较高程度的改善，成像结果更加聚焦，消除了中深层位置处绕射体的“画弧”现象，绕射体得到了很好的偏移归位。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (6)

1.一种基于纯波延拓方程的弹性反射波旅行时反演方法，其特征在于，包括如下步骤：
(1)将近似弹性波解耦延拓方程省略交叉项得到纯波延拓方程，基于纯波延拓方程计算背景波场；
(2)利用得到的背景波场计算逆时偏移结果；
(3)基于纯波延拓方程和逆时偏移结果求解扰动纯波延拓方程，再利用扰动纯波延拓方程计算扰动波场；
(4)基于步骤(1)得到的背景波场和步骤(3)得到的扰动波场，分别计算纵横波梯度；
(5)根据纵横波梯度计算步长大小，进行纵横波速度更新，当迭代值收敛到全局最小值时，输出更新的纵横波速度即为最终的反演结果；
步骤(1)具体如下：纯波延拓方程包括纯P波延拓方程和纯S波延拓方程，分别基于纯P波延拓方程计算正向传播的P波振动速度背景波场和反向传播的P波振动速度背景波场，基于纯S波延拓方程计算反向传播的S波振动速度背景波场，基于纯P波延拓方程计算反向传播的P波应力背景波场，基于纯S波延拓方程计算反向传播的S波应力背景波场，背景波场的数值模拟计算采用高阶交错网格有限差分实施；
纯P波延拓方程为：
Figure FDA0003551905310000011
其中，下标“0”表示背景，x和z分别表示空间水平方向和垂直方向，t表示时间，vPx0和vPz0分别表示P波振动速度背景波场的水平分量和垂直分量，τPxx0和τPzz0分别表示xx和zz方向的P波应力背景波场，ρ0和VP0分别表示背景密度和背景纵波速度；
纯S波延拓方程为：
Figure FDA0003551905310000012
其中，vSx0和vSz0分别表示S波振动速度背景波场的水平分量和垂直分量，τSxz0和τSzx0分别表示xz和zx方向的S波应力背景波场，VS0表示背景横波速度。
2.根据权利要求1所述的一种基于纯波延拓方程的弹性反射波旅行时反演方法，其特征在于，步骤(2)具体如下：分别用计算得到的正向传播的P波振动速度背景波场和反向传播的P波振动速度背景波场互相关计算PP波逆时偏移结果，用正向传播的P波振动速度背景波场和反向传播的S波振动速度背景波场计算PS波逆时偏移结果。
3.根据权利要求2所述的一种基于纯波延拓方程的弹性反射波旅行时反演方法，其特征在于，步骤(3)具体如下：以背景波场和速度扰动的相互作用为震源构建Born正演方程，即扰动的纯波延拓方程，扰动的纯波延拓方程包括扰动的纯P波延拓方程和扰动的纯S波延拓方程；其中，速度扰动用逆时偏移结果代替；
基于扰动的纯P波延拓方程计算正向传播的P波振动速度扰动波场和反向传播的P波应力扰动波场，基于扰动的纯S波延拓方程计算正向传播的S波振动速度扰动波场。
4.根据权利要求3所述的一种基于纯波延拓方程的弹性反射波旅行时反演方法，其特征在于，扰动的纯P波延拓方程为：
Figure FDA0003551905310000021
其中，符号“δ”表示扰动，ρ0表示背景密度，δVP表示纵波速度扰动，δVS表示横波速度扰动，vPx0和vPz0分别表示P波振动速度背景波场的水平分量和垂直分量，VP0表示背景纵波速度，δvPx和δvPz分别表示P波振动速度扰动波场的水平分量和垂直分量，δτPxx和δτPzz分别表示xx和zz方向的P波应力扰动波场；在实施过程中，方程(3)中的δVP和δVS用PP波逆时偏移结果表示；
扰动的纯S波延拓方程为：
Figure FDA0003551905310000031
在实施过程中，方程(4)中的δVS用PS波逆时偏移结果表示，δvSx和δvSz分别表示S波振动速度扰动波场的水平分量和垂直分量，δτSxz和δτSzx分别表示xz和zx方向的S波应力扰动波场，VS0表示背景横波速度。
5.根据权利要求1所述的一种基于纯波延拓方程的弹性反射波旅行时反演方法，其特征在于，步骤(4)具体方法如下：读取步骤(1)计算的背景波场和步骤(3)计算的扰动波场，分别求取纵波速度梯度
Figure FDA0003551905310000032
和横波速度梯度
Figure FDA0003551905310000033
整个计算过程中保持背景密度ρ0为常数；
Figure FDA0003551905310000034
其中，χP和χS分别表示匹配P波数据和匹配S波数据的目标函数，VP0表示背景纵波速度，vPx0和vPz0分别表示正向传播P波振动速度背景波场的水平分量和垂直分量，δvPx和δvPz分别表示正向传播P波振动速度扰动波场的水平分量和垂直分量，
Figure FDA0003551905310000035
和
Figure FDA0003551905310000036
分别表示反向传播P波应力扰动波场的xx分量和zz分量，
Figure FDA0003551905310000037
和
Figure FDA0003551905310000038
分别表示反向传播P波应力背景波场的xx分量和zz分量，VS0表示背景横波速度，δvSx和δvSz分别表示正向传播S波振动速度扰动波场的水平分量和垂直分量，
Figure FDA0003551905310000039
和
Figure FDA00035519053100000310
分别表示反向传播S波应力背景波场的xz分量和zx分量。
6.根据权利要求1所述的一种基于纯波延拓方程的弹性反射波旅行时反演方法，其特征在于，步骤(5)中，保证每次迭代的纵波速度更新量在20m/s-100m/s之间，横波速度更新量在10m/s-60m/s之间。
CN202110260212.3A 2021-03-10 2021-03-10 一种基于纯波延拓方程的弹性反射波旅行时反演方法 Active CN113126151B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110260212.3A CN113126151B (zh) 2021-03-10 2021-03-10 一种基于纯波延拓方程的弹性反射波旅行时反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110260212.3A CN113126151B (zh) 2021-03-10 2021-03-10 一种基于纯波延拓方程的弹性反射波旅行时反演方法

Publications (2)

Publication Number Publication Date
CN113126151A CN113126151A (zh) 2021-07-16
CN113126151B true CN113126151B (zh) 2022-06-07

Family

ID=76772859

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110260212.3A Active CN113126151B (zh) 2021-03-10 2021-03-10 一种基于纯波延拓方程的弹性反射波旅行时反演方法

Country Status (1)

Country Link
CN (1) CN113126151B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101358867A (zh) * 2008-09-03 2009-02-04 山东省科学院海洋仪器仪表研究所 一种海洋水位实时监测***
CN103293553A (zh) * 2013-04-17 2013-09-11 中国海洋石油总公司 一种复杂海底上下缆地震采集数据边界元延拓校正方法
CN108873066A (zh) * 2018-06-26 2018-11-23 中国石油大学(华东) 弹性介质波动方程反射波旅行时反演方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090257308A1 (en) * 2008-04-11 2009-10-15 Dimitri Bevc Migration velocity analysis methods
CN111983703B (zh) * 2020-07-24 2023-07-25 中国石油天然气集团有限公司 井间电磁测量流体成像方法、***及装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101358867A (zh) * 2008-09-03 2009-02-04 山东省科学院海洋仪器仪表研究所 一种海洋水位实时监测***
CN103293553A (zh) * 2013-04-17 2013-09-11 中国海洋石油总公司 一种复杂海底上下缆地震采集数据边界元延拓校正方法
CN108873066A (zh) * 2018-06-26 2018-11-23 中国石油大学(华东) 弹性介质波动方程反射波旅行时反演方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于能量密度的自解耦互相关成像条件;张晓语，等;《地球物理学报》;20181231;第61卷(第12期);4965-4975 *

Also Published As

Publication number Publication date
CN113126151A (zh) 2021-07-16

Similar Documents

Publication Publication Date Title
CN108873066B (zh) 弹性介质波动方程反射波旅行时反演方法
KR102020759B1 (ko) Q-보상된 전 파동장 반전
Etgen et al. Computational methods for large-scale 3D acoustic finite-difference modeling: A tutorial
US8406081B2 (en) Seismic imaging systems and methods employing tomographic migration-velocity analysis using common angle image gathers
US20050174886A1 (en) Method for processing borehole seismic data
AU2015383134B2 (en) Multistage full wavefield inversion process that generates a multiple free data set
CN110007340B (zh) 基于角度域直接包络反演的盐丘速度密度估计方法
WO2017162731A1 (en) Method of operating a data-processing system for the simulation of the acoustic wave propagation in the transversely isotropic media comprising an hydrocarbon reservoir
CN110187382A (zh) 一种回折波和反射波波动方程旅行时反演方法
CN109946742A (zh) 一种TTI介质中纯qP波地震数据模拟方法
Datta et al. Full-waveform inversion of salt models using shape optimization and simulated annealing
Hu An improved immersed boundary finite-difference method for seismic wave propagation modeling with arbitrary surface topography
CN114460646B (zh) 一种基于波场激发近似的反射波旅行时反演方法
Bekar et al. Solving the eikonal equation for compressional and shear waves in anisotropic media using peridynamic differential operator
US5986973A (en) Energy minimization in surface multiple attenuation
CN111665556A (zh) 地层声波传播速度模型构建方法
Jia et al. Superwide-angle one-way wave propagator and its application in imaging steep salt flanks
CN113126151B (zh) 一种基于纯波延拓方程的弹性反射波旅行时反演方法
CN107340537A (zh) 一种p-sv转换波叠前逆时深度偏移的方法
Jiang et al. Full waveform inversion based on inversion network reparameterized velocity
Jin et al. 2D multiscale non‐linear velocity inversion
WO2017189180A1 (en) Fwi with areal and point sources
Moura et al. Progressive matching optimisation method for FWI
Saraswat et al. Simultaneous stochastic inversion of prestack seismic data using hybrid evolutionary algorithm
Bao et al. Elastic wave pre-stack reverse-time migration based on the second-order P-and S-wave decoupling equation

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