CN113376689B - 一种考虑层间多次波的弹性反射波走时反演方法 - Google Patents

一种考虑层间多次波的弹性反射波走时反演方法 Download PDF

Info

Publication number
CN113376689B
CN113376689B CN202110340994.1A CN202110340994A CN113376689B CN 113376689 B CN113376689 B CN 113376689B CN 202110340994 A CN202110340994 A CN 202110340994A CN 113376689 B CN113376689 B CN 113376689B
Authority
CN
China
Prior art keywords
wave
inversion
elastic
model
time
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
CN202110340994.1A
Other languages
English (en)
Other versions
CN113376689A (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.)
China Railway Design Corp
Original Assignee
China Railway Design Corp
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 China Railway Design Corp filed Critical China Railway Design Corp
Priority to CN202110340994.1A priority Critical patent/CN113376689B/zh
Publication of CN113376689A publication Critical patent/CN113376689A/zh
Application granted granted Critical
Publication of CN113376689B publication Critical patent/CN113376689B/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/303—Analysis for determining velocity profiles or travel times
    • 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/51—Migration
    • G01V2210/512—Pre-stack
    • 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
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6222—Velocity; travel time

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

本发明公开了一种考虑层间多次波的弹性反射波走时反演方法，包括以下步骤：输入反射地震记录以及建立的线性增加的线性初始速度模型、开展最小二乘真振幅偏移获取弹性波偏移剖面、进行震源侧非线性Born正演和检波点侧非线性Borm正演获取含多次波的背向扰动波场、建立时移互相关目标函数，并利用模型参数梯度计算公式构建反演参数梯度、将所得反演参数梯度对线性初始速度模型进行迭代更新等步骤。该方法给出了弹性多次反射波走时反演技术框架，并对结构复杂的Marmousi模型进行了测试，多次反射波反演结果能够大致刻画Marmousi模型轮廓，以此开展的常规FWI结果更加稳健。

Description

一种考虑层间多次波的弹性反射波走时反演方法
技术领域
本发明涉及地球物理勘探叠前速度参数反演技术领域，特别是涉及一种考虑层间多次波的弹性反射波走时反演方法。
背景技术
反射波反演(Reflection waveform inversion:RWI)的概念最初是由Xu Sheng在2012年EAGE会议上提出的，主要思路是采纳偏移/反偏移技术，利用次级震源在极其平滑的背景速度场中生成反射波能量，进而对速度场进行低波数迭代更新。经过近十年的发展，反射波反演技术逐渐从模型试算走向了实际生产，成为了目前最受欢迎的初始速度建模工具。在反射波反演的发展历程中出现了很多新的技术方法，譬如为避免真振幅偏移降低计算量，发展了基于走时互相关目标函数和不依赖于振幅目标函数的反射波反演方法；为减少梯度中的高波数偏移分量，在反射波反演中引入了上、下行波分离和角度域滤波等技术；为适用于多波多维度地震勘探，反射波反演也从最初的声介质拓展至弹性介质。然而，所有这些反演方法的基础均是一阶Born近似，能够有效利用一次反射波信息，但无法利用地震记录中蕴含的多次散射能量，因此传统RWI技术对多次波发育、地质结构复杂区域的反演显得无能为力。
作为目前地震勘探领域获取深度域速度场最为先进的技术手段之一，FWI长期受制于初始输入的准确性。随着初始速度建模方法的不断发展，反射波反演已成为低波数模型构建的主流方法。但是该类方法始终没有脱离一阶Born近似的局限，反演信息主要来源于一次反射波，换言之，传统反射波反演对模型低波数信息的恢复非常有限。实际情况下，地震记录中往往包含多次反射能量，能够提供更为丰富的反射路径及反射角度，可以有效拓宽反演波数覆盖范围。因此，如何在反射波反演中充分利用多次波信号是一个非常有意义的探索方向，这也是当前亟需解决的主要问题之一。
反射波反演技术利用反射波的传播路径进行低波数成分恢复，是目前勘探领域较为常用的方法。然而传统反射波反演方法受限于线性一阶Born近似，仅利用了一次反射波路径，在构造复杂的多次波发育地区很难得到较好的背景速度模型。因此，有必要发明一种适用于构造复杂区域的新型反射波反演技术。
发明内容
本发明的目的是针对现有反射波反演技术均未能利用层间多次反射信息，波路径覆盖范围比较有限，不适用于复杂地质情况的技术缺陷，而提供一种考虑层间多次波的弹性反射波走时反演方法。
为实现本发明的目的所采用的技术方案是：
一种考虑层间多次波的弹性反射波走时反演方法，包括以下步骤：
步骤1：输入反射地震记录以及建立的线性增加的线性初始速度模型m；
步骤2：开展最小二乘真振幅偏移，获取弹性波偏移剖面I；
步骤3：进行震源侧非线性Born正演和检波点侧非线性Borm正演，获取含多次波的背向扰动波场；
上述过程表示为：
Lδus＝(us+δus)I (1)
Lδur＝(ur+δur)I (2)
其中，L为弹性波动方程算子，δu为背向扰动波场，u为背景波场，下标s和r分别表示震源侧和检波点侧；
步骤4：建立时移互相关目标函数，并利用模型参数梯度计算公式构建反演参数梯度；
所述模型参数梯度计算公式为：
Figure BDA0002999579340000021
其中，us(x，t)和ur(x，t)分别表示空间位置点x处t时刻的震源波场与检波波场，τ∈[-T，T]为时移量，T为最大时移量；
上标·和··分别代表一阶时间导数和二阶时间导数；
步骤5：将所得反演参数梯度对线性初始速度模型进行迭代更新；
步骤6：判断迭代更新后的速度模型是否满足收敛条件；若满足收敛条件，进入步骤7；若不满足收敛条件，重复步骤2-6；
步骤7：输出低波数速度场。
在上述技术方案中，步骤4中，所述时移互相关目标函数定义为：
Figure BDA0002999579340000031
在上述技术方案中，步骤4中，所得反演参数梯度采取波数域滤波的方式对高波数成分进行剔除。
与现有技术相比，本发明的有益效果是：
1.本发明提供的考虑层间多次波的弹性反射波走时反演方法，考虑到地震波的走时属性和反射波路径均与模型参数的低波数成分密切相关，本发明联合走时与反射波恢复速度中的低波数成分，在初始模型远离真实模型的情况下，可最大限度提高反演的稳定性。
2.本发明提供的考虑层间多次波的弹性反射波走时反演方法，结合非线性Born正演算子，提出了弹性多次反射波走时反演方法，能够提高照明强度，扩大反演结果中的波数覆盖范围。
附图说明
图1所示为实施例1中的方法流程图；
图2：图2a为线性反偏移示意图；图2b为非线性反偏移示意图。
图3：图3a为部分Marmousi模型的真实纵波速度；图3b为Marmousi模型的真实横波速度；图3c为纵波速度的线性初始输入；图3d为横波速度的线性初始输入。
图4a为利用对比例1中的方法输出的低波数纵波速度场；图4b为利用对比例1中的方法输出的低波数横波速度场。图4c为利用实施例1中的方法输出的低波数纵波速度场；图4d为利用实施例1中的方法输出的低波数横波速度场；
图5：图5a为图4a进行全波形反演后所得的图像；图5b为图4b进行全波形反演后所得的图像；图5c为图4c进行全波形反演后所得的图像；图5d为图4d进行全波形反演后所得的图像。
具体实施方式
以下结合具体实施例对本发明作进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
实施例1
一种考虑层间多次波的弹性反射波走时反演方法，包括以下步骤：
步骤1：输入反射地震记录以及建立的线性增加的线性初始速度模型m；
步骤2：开展最小二乘真振幅偏移，获取弹性波偏移剖面I；
由于非线性反偏移的次级震源能量与弹性波偏移剖面I的幅值之间的正相关关系，当弹性波偏移剖面I仅能表示地下反射结构而无法反映反射系数绝对值时，背向扰动波场和背景波场会存在数量级上的差异，这将会导致后续过程中方程(1)(2)传播异常。因此，在考虑层间多次波的弹性反射波走时反演方法中采用最小二乘真振幅偏移来提供幅值准确的反射系数剖面是非常有必要的。
步骤3：进行震源侧非线性Born正演和检波点侧非线性Borm正演(即非线性反偏移)，获取含多次波的背向扰动波场δu；上述过程表示为：
Lδus＝(us+δus)I (1)
Lδur＝(ur+δur)I (2)
其中，L为弹性波动方程算子，δu为背向扰动波场，u为背景波场，下标s和r分别表示震源侧和检波点侧；
图2a所示为线性反偏移示意图，其中仅包含了一次反射能量。图2b所示为非线性反偏移示意图，其中包含一次反射和多次反射能量。非线性反偏移与线性反偏移的主要不同是一阶扰动波场的二次利用，从而产生更高阶的反射波，从而提高照明强度。
步骤4：建立时移互相关目标函数，并利用模型参数梯度计算公式构建反演参数梯度，即计算参数更新量；所述时移互相关目标函数定义为：
Figure BDA0002999579340000041
其中，us(x，t)和ur(x，t)分别表示空间位置点x处t时刻的震源波场与检波波场，τ∈[-T，T]为时移量，T为最大时移量；
所述模型参数梯度计算公式为：
Figure BDA0002999579340000042
其中上标·和··分别代表一阶时间导数和二阶时间导数；式(3)由三部分组成，其中前两个积分项分别表示地下成像点位置与震源点和检波点之间的波传播路径(包括一次反射波路径和多次反射波路径)；第三项表示地下空间位置相互之间的多次传播路径。除低波数层析分量之外，梯度中难免会存在小角度反射互相关导致的高波数偏移分量，可采取波数域滤波的方式对高波数成分进行剔除。
步骤5：将所得反演参数梯度对线性初始速度模型m进行迭代更新；
步骤6：判断迭代更新后的速度模型是否满足收敛条件；若满足收敛条件，进入步骤7；若不满足收敛条件，重复步骤2-6；
步骤7：输出低波数速度场。
对比例1
一种反射波反演方法，包括以下步骤：
步骤1：输入反射地震记录，建立线性增加的线性初始速度模型m；
步骤2：利用线性初始速度模型m进行偏移，获取弹性波偏移剖面I，此时I为源检波场的零延迟互相关；
步骤3：通过反偏移技术生成扰动波场，该过程表示为：
Lδus＝usI (5)
Lδur＝urI (6)
其中，公式(5)和(6)分别代表震源侧和检波点侧的反偏移过程，下标s和r分别表示震源侧和检波点侧，L为弹性波动方程算子，δu为扰动波场，u为背景波场；
上述反偏移过程是一个明显的线性过程，方程右侧次级震源uI产生的一阶扰动波场δu将直接传播到地表，不会与弹性波偏移剖面I中的任何构造相互作用进一步激发产生多次波，如图2a所示，由于仅包含一次反射波，其照明强度较低。
步骤4：根据伴随状态法以及隐式函数求导法则，求取线性初始速度模型m的层析更新量；所述层析更新量表示为：
Figure BDA0002999579340000051
其中上标·和··分别代表一阶时间导数和二阶时间导数；
步骤5：利用最优化算法对线性初始速度模型m进行迭代更新，并输出速度场模型。
上述公式(7)中的第一项和第二项分别代表成像点位置到震源和检波点的一次反射波路径，在模型比较简单的情况下对低波数成分的恢复效果比较明显，但不适用于复杂多变的地质结构。
对比实施例1中的步骤3和对比例中的步骤3，实施例1中的步骤3所用公式是让线性反偏移中生成的一阶扰动波场继续与弹性波偏移剖面I相互作用形成非线性的闭环。能够产生层间多次反射波，使反射波反演方法适用于复杂地下构造。
实施例2
本实施例是基于实施例1和对比例1的方法进行的实际应用。
图3所示为实施例1和对比例1中输入的反射地震记录以及建立的线性增加的线性初始速度模型m。其中，图3a为部分Marmousi模型的真实纵波速度；图3b为部分Marmousi模型的真实横波速度；依据真实纵波速度和真实横波速度获得反射地震记录。图3c为纵波速度的线性初始输入；图3d为横波速度的线性初始输入。纵波速度的线性初始输入和横波速度的线性初始输入组成线性初始速度模型m。
图4为实施例1和对比例1输出的低波数速度场模型。其中，图4a为利用对比例1中的方法输出的低波数纵波速度场；图4b为利用对比例1中的方法输出的低波数横波速度场。图4c为利用实施例1中的方法输出的低波数纵波速度场；图4d为利用实施例1中的方法输出的低波数横波速度场。
图4a、4b相比于图4c、4d，由于图4c、4d充分利用了层间多次波，因此其低波数成分相比于图4a、4b更加丰富。
为更加形象直观的验证图4a、4b与4c、4d之间的区别，分别对图4a、4b和4c、4d进行全波形反演，结果如图5所示。图5a、5b为图4a、4b进行全波形反演后所得的图像；图5c、5d为图4c、4d进行全波形反演后所得的图像。图5c相比于图5a以及图5d相比于图5b，其中下方更加清晰。由此可以验证，图4c、4d相比于图4a、4b其低波数更加丰富。由此可知，实施例1中的方法更加适合复杂构造地区低波数信息重建。
以上所述仅是本发明的优选实施方式，应当指出的是，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (3)

1.一种考虑层间多次波的弹性反射波走时反演方法，其特征在于：包括以下步骤：
步骤1：输入反射地震记录以及建立的线性增加的线性初始速度模型m；
步骤2：开展最小二乘真振幅偏移，获取弹性波偏移剖面I；
步骤3：进行震源侧非线性Born正演和检波点侧非线性Borm正演，获取含多次波的背向扰动波场；
上述过程表示为：
Lδus＝(us+δus)I (1)
Lδur＝(ur+δur)I (2)
其中，L为弹性波动方程算子，δu为背向扰动波场，u为背景波场，下标s和r分别表示震源侧和检波点侧；
步骤4：建立时移互相关目标函数，并利用模型参数梯度计算公式构建反演参数梯度；
所述模型参数梯度计算公式为：
Figure FDA0002999579330000011
其中，us(x，t)和ur(x，t)分别表示空间位置点x处t时刻的震源波场与检波场，τ∈[-T，T]为时移量，T为最大时移量；
上标·和··分别代表一阶时间导数和二阶时间导数；
步骤5：将所得反演参数梯度对线性初始速度模型进行迭代更新；
步骤6：判断迭代更新后的速度模型是否满足收敛条件；若满足收敛条件，进入步骤7；若不满足收敛条件，重复步骤2-6；
步骤7：输出低波数速度场。
2.如权利要求1所述的弹性反射波走时反演方法，其特征在于：步骤4中，所述时移互相关目标函数定义为：
Figure FDA0002999579330000012
。
3.如权利要求1所述的弹性反射波走时反演方法，其特征在于：步骤4中，所得反演参数梯度采取波数域滤波的方式对高波数成分进行剔除。
CN202110340994.1A 2021-03-30 2021-03-30 一种考虑层间多次波的弹性反射波走时反演方法 Active CN113376689B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110340994.1A CN113376689B (zh) 2021-03-30 2021-03-30 一种考虑层间多次波的弹性反射波走时反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110340994.1A CN113376689B (zh) 2021-03-30 2021-03-30 一种考虑层间多次波的弹性反射波走时反演方法

Publications (2)

Publication Number Publication Date
CN113376689A CN113376689A (zh) 2021-09-10
CN113376689B true CN113376689B (zh) 2022-04-12

Family

ID=77570620

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110340994.1A Active CN113376689B (zh) 2021-03-30 2021-03-30 一种考虑层间多次波的弹性反射波走时反演方法

Country Status (1)

Country Link
CN (1) CN113376689B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113791447B (zh) * 2021-10-12 2023-06-20 同济大学 一种反射结构导引的反射波层析反演方法
CN117630174A (zh) * 2024-01-25 2024-03-01 中国铁路设计集团有限公司 板式混凝土多通道-多自由度脉冲波无损检测方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104391323A (zh) * 2014-11-21 2015-03-04 中国石油大学(华东) 一种利用反射波信息反演速度场中低波数成分的方法
CN108241173A (zh) * 2017-12-28 2018-07-03 中国石油大学(华东) 一种地震资料偏移成像方法及系统
CN108873066A (zh) * 2018-06-26 2018-11-23 中国石油大学(华东) 弹性介质波动方程反射波旅行时反演方法
CN112099082A (zh) * 2019-06-17 2020-12-18 中国海洋大学 一种共面元共方位角道集的地震回折波走时反演方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10459096B2 (en) * 2016-07-13 2019-10-29 Exxonmobil Upstream Research Company Joint full wavefield inversion of P-wave velocity and attenuation using an efficient first order optimization
CN107203002B (zh) * 2017-06-12 2019-05-24 中国科学院地质与地球物理研究所 反演速度模型的建立方法和地下结构的像的获得方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104391323A (zh) * 2014-11-21 2015-03-04 中国石油大学(华东) 一种利用反射波信息反演速度场中低波数成分的方法
CN108241173A (zh) * 2017-12-28 2018-07-03 中国石油大学(华东) 一种地震资料偏移成像方法及系统
CN108873066A (zh) * 2018-06-26 2018-11-23 中国石油大学(华东) 弹性介质波动方程反射波旅行时反演方法
CN112099082A (zh) * 2019-06-17 2020-12-18 中国海洋大学 一种共面元共方位角道集的地震回折波走时反演方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
解耦纵横波反射波走时反演;宋建勇等;《地球物理学报》;20190831;第62卷(第8期);第3130-3139页 *

Also Published As

Publication number Publication date
CN113376689A (zh) 2021-09-10

Similar Documents

Publication Publication Date Title
AU2016331881B8 (en) Q-compensated full wavefield inversion
EP2715405B1 (en) Method of processing seismic data by providing surface offset common image gathers
CN107817526B (zh) 叠前地震道集分段式振幅能量补偿方法及系统
CN113376689B (zh) 一种考虑层间多次波的弹性反射波走时反演方法
Guitton et al. A parameterization study for elastic VTI full-waveform inversion of hydrophone components: Synthetic and North Sea field data examples
Yao et al. Reflection-waveform inversion regularized with structure-oriented smoothing shaping
Sun et al. Born modeling based adjustive reflection full waveform inversion
CN113740901A (zh) 基于复杂起伏地表的陆上地震数据全波形反演方法及装置
Beasley et al. Wave equation receiver deghosting: A provocative example
RU2570827C2 (ru) Гибридный способ для полноволновой инверсии с использованием способа одновременных и последовательных источников
Raknes et al. Seismic imaging of the carbon dioxide gas cloud at Sleipner using 3D elastic time-lapse full waveform inversion
Alam Near-surface characterization using traveltime and full-waveform inversion with vertical and horizontal component seismic data
Jiao et al. Elastic migration for improving salt and subsalt imaging and inversion
Yang et al. Mitigating velocity errors in least-squares imaging using angle-dependent forward and adjoint Gaussian beam operators
Zhu et al. Scattering effect on shallow gas-obscured zone imaging in Bohai PL19-3 area
Gonçalves et al. Flexible layer-based 2D refraction tomography method for statics corrections
Ni* et al. Preliminary practice of stereotomography
Guerra et al. Wave-equation tomography using image-space phase encoded data
Galuzzi Modelling and optimization techniques for acoustic Full Waveform Inversion in seismic exploration
Browett et al. Campeche 3D WAZ imaging challenges
Moussavi Alashloo et al. Least-Squares Reverse Time Migration Using Generalised Diffraction-Stack Imaging Condition
Kakhkhorov Seismic imaging with primaries and multiples
Arias et al. Imaging of thrust structures under Colombian foothills landscape through Marchenko approach
Tian et al. The effect of initial models on seismic impedance inversion accuracy based on model
Nejati et al. Migrated Exploding Reflectors in Evaluation of Finite Difference Solution for Inhomogeneous Seismic Models

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