CN114460646B - 一种基于波场激发近似的反射波旅行时反演方法 - Google Patents

一种基于波场激发近似的反射波旅行时反演方法 Download PDF

Info

Publication number
CN114460646B
CN114460646B CN202210381052.2A CN202210381052A CN114460646B CN 114460646 B CN114460646 B CN 114460646B CN 202210381052 A CN202210381052 A CN 202210381052A CN 114460646 B CN114460646 B CN 114460646B
Authority
CN
China
Prior art keywords
time
wave field
excitation
wave
background
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
CN202210381052.2A
Other languages
English (en)
Other versions
CN114460646A (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 CN202210381052.2A priority Critical patent/CN114460646B/zh
Publication of CN114460646A publication Critical patent/CN114460646A/zh
Application granted granted Critical
Publication of CN114460646B publication Critical patent/CN114460646B/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
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • 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
    • 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

本发明属于海洋地震勘探领域，具体公开了一种基于波场激发近似的反射波旅行时反演方法，包括如下步骤：以给定的地震子波和观测地震数据，结合初始速度模型计算背景波场；以观测地震数据作为震源，计算伴随背景波场，计算逆时偏移结果；恢复背景波场，计算扰动波场；基于扰动波场建立反射波旅行时反演的梯度方程，计算梯度，根据梯度确定步长大小，得到本次迭代的速度更新量；迭代更新速度参数，直到满足收敛条件，输出速度数据即为最终的反演结果。本发明所公开的方法只保存波场的激发振幅和激发时间，可以提高反射波旅行时反演的效率，大幅度降低存储消耗，为海洋地震勘探背景速度建模提供技术支持。

Description

一种基于波场激发近似的反射波旅行时反演方法
技术领域
本发明属于海洋地震勘探领域，特别涉及一种基于波场激发近似的反射波旅行时反演方法。
背景技术
海洋地震反演通过观测到的地震数据估算相应的地球物理参数，进而反推海底地下的结构形态及物质成分，可以有效地识别地质构造、预测自然灾害和勘探油气藏。地震波的传播速度不仅是处理、解释地震数据的主要依据，而且其本身就是反映底下介质构造和岩性的重要数据。因此，求取地震波在地下传播的精确速度显得尤为重要。
反射波旅行时反演以波动理论为基础，匹配模拟地震反射数据与观测地震反射数据的旅行时信息建立目标函数，利用局部最优化的方法迭代寻找目标函数的全局最小值，可以有效的恢复海底底质的背景速度，可以为全波形反演、逆时偏移等反演、成像技术提供准确的初始速度模型，是目前海洋地震勘探领域的前沿方向之一。但是，由于反射波旅行时反演算法需要引入偏移/反偏移过程来计算扰动波场，每次迭代过程需要同时存储背景波场和扰动波场的数据，海量的存储需求将产生巨大的时间消耗，制约了反射波旅行时反演算法走向实际生产应用的步伐。
发明内容
为解决上述技术问题，本发明提供了一种基于波场激发近似的反射波旅行时反演方法，以达到降低波场存储消耗，提高反演效率的目的。
为达到上述目的，本发明的技术方案如下：
一种基于波场激发近似的反射波旅行时反演方法，包括如下步骤：
（1）以给定的地震子波和观测地震数据，结合初始速度模型计算背景波场，并保存背景波场的激发振幅和激发时间；
（2）以观测地震数据作为震源，计算伴随背景波场，读取步骤（1）中的激发时间计算逆时偏移结果；
（3）读取步骤（1）中的激发振幅恢复背景波场，结合步骤（2）中的逆时偏移结果，计算扰动波场，保存上行方向扰动波场的激发振幅和激发时间；
（4）基于扰动波场建立反射波旅行时反演的梯度方程，计算梯度，根据梯度确定步长大小，得到本次迭代的速度更新量；
（5）迭代更新速度参数，直到满足收敛条件，输出速度数据即为最终的反演结果。
上述方案中，步骤（1）具体如下：输入已知的观测地震数据、初始速度模型和震源 子波，采用时间2阶、空间10阶的交错网格有限差分方法数值求解声波方程（1），得到背景波 场
Figure 50593DEST_PATH_IMAGE001
，并保存背景波场的激发振幅
Figure 597987DEST_PATH_IMAGE002
和激发时间
Figure 608669DEST_PATH_IMAGE003
；
Figure 451991DEST_PATH_IMAGE004
（1）
其中，
Figure 924954DEST_PATH_IMAGE005
表示声波传播算子，
Figure 636558DEST_PATH_IMAGE006
表示震源子波矢量。
上述方案中，步骤（2）具体如下：以观测地震数据替换声波方程（1）中的震源子波，采用时间2阶、空间10阶的交错网格有限差分方法数值求解声波方程（1），计算伴随背景波场，基于步骤（1）中的激发时间互相关背景波场和伴随背景波场，得到逆时偏移结果。
上述方案中，步骤（3）具体如下：
读取步骤（1）中的激发振幅
Figure 377112DEST_PATH_IMAGE007
，根据方程（2）恢复背景波场
Figure 515969DEST_PATH_IMAGE001
：
Figure 473299DEST_PATH_IMAGE008
（2）
其中，符号“*”表示褶积，
Figure 723014DEST_PATH_IMAGE006
表示震源子波矢量；
采用时间2阶、空间10阶的交错网格有限差分方法数值求解Born正演方程（3），计 算扰动波场
Figure 583654DEST_PATH_IMAGE009
，保存上行方向扰动波场的激发振幅和激发时间；
Figure 893413DEST_PATH_IMAGE010
（3）
其中，
Figure 964137DEST_PATH_IMAGE011
表示速度扰动，这里，速度扰动用逆时偏移结果代替；
扰动波场的计算：
Figure 394375DEST_PATH_IMAGE012
（4）
其中，
Figure 234155DEST_PATH_IMAGE005
表示传播算子，
Figure 714815DEST_PATH_IMAGE011
表示模型参数扰动算子，
Figure 882622DEST_PATH_IMAGE013
表示扰动波场的激发振幅。
上述方案中，步骤（4）具体如下：定义反射波旅行时反演的L2范数目标函数，依据目标函数得到基于波场激发近似的梯度方程，给予合适的步长，得到本次迭代的速度更新量。
进一步的技术方案中，定义反射波旅行时反演的L2范数目标函数：
Figure 739719DEST_PATH_IMAGE014
（5）
其中，
Figure 434006DEST_PATH_IMAGE015
表示目标函数，
Figure 85567DEST_PATH_IMAGE016
表示背景速度参数，
Figure 504785DEST_PATH_IMAGE017
表示观测地震数据与模拟地震数 据的旅行时差异，
Figure 899994DEST_PATH_IMAGE018
表示与
Figure 448787DEST_PATH_IMAGE017
有关的L2范数，
Figure 677774DEST_PATH_IMAGE017
由互相关辅助方程（6）计算：
Figure 944808DEST_PATH_IMAGE019
（6）
其中，
Figure 254960DEST_PATH_IMAGE020
和
Figure 658259DEST_PATH_IMAGE021
分别表示检波点位置处的观测地震数据与模拟地震数据，模拟地震 数据
Figure 917202DEST_PATH_IMAGE021
来自于扰动波场
Figure 546898DEST_PATH_IMAGE022
，
Figure 283910DEST_PATH_IMAGE023
表示时间；当方程（6）的互相关值达到最大时，对应的
Figure 807295DEST_PATH_IMAGE024
等于
Figure 611040DEST_PATH_IMAGE017
；
根据方程（5）和（6），基于伴随状态法推导反射波旅行时反演的梯度方程，并将方程（2）和（4）代入梯度方程得方程（7）：
Figure 587087DEST_PATH_IMAGE025
Figure 255353DEST_PATH_IMAGE026
Figure 898824DEST_PATH_IMAGE027
（7）
其中，
Figure 607892DEST_PATH_IMAGE028
表示梯度，
Figure 336813DEST_PATH_IMAGE029
表示两个矢量的内积运算，
Figure 556573DEST_PATH_IMAGE030
和
Figure 165802DEST_PATH_IMAGE031
分别表示背景波场 和扰动波场的二阶时间导数，
Figure 671870DEST_PATH_IMAGE032
和
Figure 888088DEST_PATH_IMAGE033
分别表示伴随背景波场和伴随扰动波场，
Figure 19861DEST_PATH_IMAGE034
表示震源 子波矢量的二阶时间导数，
Figure 247711DEST_PATH_IMAGE035
表示零延迟互相关；伴随背景波场和伴随扰动波场分别通过求 解方程（1）和（3）计算，其中，震源子波矢量
Figure 190259DEST_PATH_IMAGE006
替换为反射波旅行时反演方法的伴随震源
Figure 628194DEST_PATH_IMAGE036
， 见方程（8）；
Figure 691221DEST_PATH_IMAGE037
（8）
根据方程（7），给定合适的步长，得到本次迭代的速度更新量。
上述方案中，步骤（5）中，保证每次迭代的速度更新量在20m/s-100m/s之间，迭代更新速度参数，直到满足收敛条件。
进一步的技术方案中，收敛条件为：计算本次迭代模拟地震数据与观测地震数据旅行时误差的平方值，与上一次迭代的值进行比较，倘若数值减小，则继续进行步骤（1），当误差的平方值连续5次处于非下降状态，则认为反演方法已经收敛到全局最小值，满足收敛条件，此时，输出速度数据即为最终的反演结果。
通过上述技术方案，本发明提供的一种基于波场激发近似的反射波旅行时反演方法具有如下有益效果：
本发明对反射波旅行时反演方法中的背景波场和扰动波场进行激发近似，只保存激发振幅和激发时间，避免了全波场的存储与读写操作，可以有效降低反射波旅行时反演的存储消耗。另外，激发振幅与震源矢量的褶积可以有效的避免波场震源特征缺失的问题，保证了波形的完整性，而上行波激发振幅的保存解决了扰动波场的波路径问题，有利于提高基于波场激发近似的反射波旅行时反演方法的精度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍。
图1是本发明所公开的一种基于波场激发近似的反射波旅行时反演方法实施的流程示意图；
图2a-2b是Sigsbee2A速度模型，其中图2a为真实的速度模型；图2b为初始的速度模型；
图3a-3b是常规反射波旅行时反演梯度与本发明基于波场激发近似的反射波旅行时反演梯度的对比，其中图3a为常规反射波旅行时反演的梯度；图3b为本发明基于波场激发近似的反射波旅行时反演的梯度；
图4a-4b是常规反射波旅行时反演的结果，其中图4a为迭代40次后反演的纵波速度；图4b为反演纵波速度对应的RTM结果；
图5a-5b是本发明基于波场激发近似的反射波旅行时反演的结果，其中图5a为迭代40次后反演的纵波速度；图5b为反演纵波速度对应的RTM结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。
本发明提供了一种基于波场激发近似的反射波旅行时反演方法，如图1所示，通过模型测试说明具体的技术方案：
（1）以给定的地震子波和观测地震数据，结合初始速度模型计算背景波场，并保存背景波场的激发振幅和激发时间；
具体如下：
输入已知的观测地震数据、初始速度模型和震源子波，采用时间2阶、空间10阶的 交错网格有限差分方法数值求解声波方程（1），得到背景波场
Figure 773578DEST_PATH_IMAGE001
，并保存背景波场的激发振 幅
Figure 621448DEST_PATH_IMAGE002
和激发时间
Figure 546679DEST_PATH_IMAGE003
；
Figure 763861DEST_PATH_IMAGE004
（1）
其中，
Figure 90937DEST_PATH_IMAGE005
表示声波传播算子，
Figure 752119DEST_PATH_IMAGE006
表示震源子波矢量。
模型测试时的地震子波一般给定雷克子波，实际应用时的地震子波一般通过子波反演技术得到；观测地震数据，主要指的是地表检波器接收的垂直分量的地震数据；初始速度模型中初始速度的速度值一般是从浅层至深层逐渐增大。
（2）以观测地震数据作为震源，计算伴随背景波场，读取步骤（1）中的激发时间计算逆时偏移结果；
具体如下：以观测地震数据替换声波方程（1）中的震源子波，采用时间2阶、空间10阶的交错网格有限差分方法数值求解声波方程（1），计算伴随背景波场，基于步骤（1）中的激发时间互相关背景波场和伴随背景波场，得到逆时偏移结果。
（3）读取步骤（1）中的激发振幅恢复背景波场，结合步骤（2）中的逆时偏移结果，计算扰动波场，保存上行方向扰动波场的激发振幅和激发时间；
具体如下：
读取步骤（1）中的激发振幅
Figure 164645DEST_PATH_IMAGE007
，根据方程（2）恢复背景波场
Figure 176333DEST_PATH_IMAGE001
：
Figure 233282DEST_PATH_IMAGE008
（2）
其中，符号“*”表示褶积，
Figure 688534DEST_PATH_IMAGE006
表示震源子波矢量；激发近似后的波场失去了震源特 征，激发振幅褶积震源信号矢量用于恢复背景波场。
采用时间2阶、空间10阶的交错网格有限差分方法数值求解Born正演方程（3），计 算扰动波场
Figure 588357DEST_PATH_IMAGE009
，保存上行方向扰动波场的激发振幅和激发时间；
Figure 655932DEST_PATH_IMAGE010
（3）
其中，
Figure 68852DEST_PATH_IMAGE011
表示速度扰动，这里，速度扰动用逆时偏移结果代替；
扰动波场的计算：
Figure 570372DEST_PATH_IMAGE012
（4）
其中，
Figure 331392DEST_PATH_IMAGE005
表示传播算子，
Figure 435614DEST_PATH_IMAGE011
表示模型参数扰动算子，
Figure 467155DEST_PATH_IMAGE013
表示扰动波场的激发振幅。
扰动波场的激发振幅褶积震源信号矢量可以恢复扰动波场的震源特征，保存上行 扰动波场的激发振幅
Figure 264210DEST_PATH_IMAGE038
，用于解决扰动波场的多路径的问题。
（4）基于扰动波场建立反射波旅行时反演的梯度方程，计算梯度，根据梯度确定步长大小，得到本次迭代的速度更新量；
具体如下：定义反射波旅行时反演的L2范数目标函数，依据目标函数得到基于波场激发近似的梯度方程，给予合适的步长，得到本次迭代的速度更新量。
定义反射波旅行时反演的L2范数目标函数：
Figure 873046DEST_PATH_IMAGE014
（5）
其中，
Figure 298736DEST_PATH_IMAGE015
表示目标函数，
Figure 683318DEST_PATH_IMAGE016
表示背景速度参数，
Figure 651274DEST_PATH_IMAGE017
表示观测地震数据与模拟地震数 据的旅行时差异，
Figure 622773DEST_PATH_IMAGE018
表示与
Figure 68797DEST_PATH_IMAGE017
有关的L2范数，
Figure 199564DEST_PATH_IMAGE017
由互相关辅助方程（6）计算：
Figure 715253DEST_PATH_IMAGE019
（6）
其中，
Figure 298681DEST_PATH_IMAGE020
和
Figure 797664DEST_PATH_IMAGE021
分别表示检波点位置处的观测地震数据与模拟地震数据，模拟地震 数据
Figure 658304DEST_PATH_IMAGE021
来自于扰动波场
Figure 968063DEST_PATH_IMAGE022
，
Figure 415618DEST_PATH_IMAGE023
表示时间；当方程（6）的互相关值达到最大时，对应的
Figure 344391DEST_PATH_IMAGE024
等于
Figure 184171DEST_PATH_IMAGE017
；
根据方程（5）和（6），基于伴随状态法推导反射波旅行时反演的梯度方程，并将方程（2）和（4）代入梯度方程得方程（7）：
Figure 38732DEST_PATH_IMAGE039
Figure 331173DEST_PATH_IMAGE026
Figure 440468DEST_PATH_IMAGE027
（7）
其中，
Figure 10121DEST_PATH_IMAGE028
表示梯度，
Figure 661682DEST_PATH_IMAGE029
表示两个矢量的内积运算，
Figure 706998DEST_PATH_IMAGE030
和
Figure 476109DEST_PATH_IMAGE031
分别表示背景波场 和扰动波场的二阶时间导数，
Figure 900268DEST_PATH_IMAGE032
和
Figure 988310DEST_PATH_IMAGE033
分别表示伴随背景波场和伴随扰动波场，
Figure 255343DEST_PATH_IMAGE034
表示震源 子波矢量的二阶时间导数，
Figure 776799DEST_PATH_IMAGE035
表示零延迟互相关；伴随背景波场和伴随扰动波场分别通过求 解方程（1）和（3）计算，其中，震源子波矢量
Figure 445677DEST_PATH_IMAGE006
替换为反射波旅行时反演方法的伴随震源
Figure 78522DEST_PATH_IMAGE036
， 见方程（8）；
Figure 567272DEST_PATH_IMAGE040
（8）
根据方程（7），每次梯度计算只需要存储背景波场的激发振幅
Figure 697427DEST_PATH_IMAGE007
和上行方向扰 动波场的激发振幅
Figure 594713DEST_PATH_IMAGE038
以及相应的激发时间，避免了整个背景波场和扰动波场的存储。根 据方程（7），给定合适的步长，得到本次迭代的速度更新量。
图2-图5以Sigsbee2A模型进行基于波场激发近似的反射波旅行时反演的数值测试，其中图2a和图2b展示了Sigsbee2A模型的初始速度和真实速度。图3a和图3b分别展示了常规反射波旅行时反演和基于波场激发近似的反射波旅行时反演的梯度结果，二者在能量上存在略微的差异，这主要因为常规反射波旅行时反演的波场是多路径的。单炮梯度的计算需求见表1。
表 1. 单炮梯度的计算存储需求对比
方法 波场维度 存储量 时间
常规方法 351*184*2000 985.47MB 57.87s
激发近似 351*184 0.99MB 47.13s
常规反射波旅行时反演梯度需要存储四次波场
Figure 758978DEST_PATH_IMAGE041
，单个波场的维度是 351*184*2000，由于波场文件可以重复利用，最终只需要保存两个波场文件，常规反射波旅 行时反演的存储量是985.47MB。激发近似反射波旅行时反演算法需要保存激发振幅和激发 时间，同样由于文件可以重复利用，最终需要保存四个文件，其存储量是0.99MB，存储量降 低了近1000倍，这主要与时间采样点数有关。
（5）迭代更新速度参数，直到满足收敛条件，输出速度数据即为最终的反演结果。
根据梯度确定步长大小，保证每次迭代的速度更新量在20m/s-100m/s之间，迭代更新速度参数，直到满足收敛条件。
收敛条件为：计算本次迭代模拟地震数据与观测地震数据旅行时误差的平方值，与上一次迭代的值进行比较，倘若数值减小，则继续进行步骤（1），当误差的平方值连续5次处于非下降状态，则认为反演方法已经收敛到全局最小值，满足收敛条件，此时，输出速度数据即为最终的反演结果。
经过40次迭代后，常规反射波旅行时反演恢复的纵波速度见图4a，图4b是图4a对应的反射波旅行时反演结果。同样经过40次迭代后，激发近似反射波旅行时反演恢复的纵波速度见图5a，可以看出，激发近似反射波旅行时反演主要更新速度模型中的低波数分量，其精度与常规反射波旅行时反演几乎没有差异，以图5a所示的纵波速度作为初始模型，得到逆时偏移结果见图5b，成像结果聚焦，中深层位置处绕射体无“画弧”现象，绕射体得到了很好的偏移归位, 这说明了我们的激发近似反射波旅行时反演可以提供低波数分量准确的纵波速度模型。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (8)

1.一种基于波场激发近似的反射波旅行时反演方法，其特征在于，包括如下步骤：
（1）以给定的地震子波和观测地震数据，结合初始速度模型计算背景波场，并保存背景波场的激发振幅和激发时间；
（2）以观测地震数据作为震源，计算伴随背景波场，读取步骤（1）中的激发时间计算逆时偏移结果；
（3）读取步骤（1）中的激发振幅恢复背景波场，结合步骤（2）中的逆时偏移结果，计算扰动波场，保存上行方向扰动波场的激发振幅和激发时间；
（4）基于扰动波场建立反射波旅行时反演的梯度方程，计算梯度，根据梯度确定步长大小，得到本次迭代的速度更新量；
（5）迭代更新速度参数，直到满足收敛条件，输出速度数据即为最终的反演结果。
2.根据权利要求1所述的一种基于波场激发近似的反射波旅行时反演方法，其特征在 于，步骤（1）具体如下：输入已知的观测地震数据、初始速度模型和震源子波，采用时间2阶、 空间10阶的交错网格有限差分方法数值求解声波方程（1），得到背景波场
Figure 475233DEST_PATH_IMAGE001
，并保存背景波 场的激发振幅
Figure 310334DEST_PATH_IMAGE002
和激发时间
Figure 162621DEST_PATH_IMAGE003
；
Figure 274934DEST_PATH_IMAGE004
（1）
其中，
Figure 203576DEST_PATH_IMAGE005
表示声波传播算子，
Figure 401339DEST_PATH_IMAGE006
表示震源子波矢量。
3.根据权利要求2所述的一种基于波场激发近似的反射波旅行时反演方法，其特征在于，步骤（2）具体如下：以观测地震数据替换声波方程（1）中的震源子波，采用时间2阶、空间10阶的交错网格有限差分方法数值求解声波方程（1），计算伴随背景波场，基于步骤（1）中的激发时间互相关背景波场和伴随背景波场，得到逆时偏移结果。
4.根据权利要求2所述的一种基于波场激发近似的反射波旅行时反演方法，其特征在于，步骤（3）具体如下：
读取步骤（1）中的激发振幅
Figure 824361DEST_PATH_IMAGE007
，根据方程（2）恢复背景波场
Figure 56759DEST_PATH_IMAGE001
：
Figure 297248DEST_PATH_IMAGE008
（2）
其中，符号“*”表示褶积，
Figure 841362DEST_PATH_IMAGE009
表示震源子波矢量；
采用时间2阶、空间10阶的交错网格有限差分方法数值求解Born正演方程（3），计算扰 动波场
Figure 301031DEST_PATH_IMAGE010
，保存上行方向扰动波场的激发振幅和激发时间；
Figure 122356DEST_PATH_IMAGE011
（3）
其中，
Figure 658380DEST_PATH_IMAGE012
表示模型参数扰动算子，这里，速度扰动用逆时偏移结果代替；
扰动波场的计算：
Figure 565156DEST_PATH_IMAGE013
（4）
其中，
Figure 329981DEST_PATH_IMAGE014
表示扰动波场的激发振幅。
5.根据权利要求4所述的一种基于波场激发近似的反射波旅行时反演方法，其特征在于，步骤（4）具体如下：定义反射波旅行时反演的L2范数目标函数，依据目标函数得到基于波场激发近似的梯度方程，给予合适的步长，得到本次迭代的速度更新量。
6.根据权利要求5所述的一种基于波场激发近似的反射波旅行时反演方法，其特征在于，定义反射波旅行时反演的L2范数目标函数：
Figure 271392DEST_PATH_IMAGE015
（5）
其中，
Figure 119262DEST_PATH_IMAGE016
表示目标函数，
Figure 637968DEST_PATH_IMAGE017
表示背景速度参数，
Figure 331118DEST_PATH_IMAGE018
表示观测地震数据与模拟地震数据的 旅行时差异，
Figure 963919DEST_PATH_IMAGE019
表示与
Figure 841745DEST_PATH_IMAGE018
有关的L2范数，
Figure 723113DEST_PATH_IMAGE018
由互相关辅助方程（6）计算：
Figure 95320DEST_PATH_IMAGE020
（6）
其中，
Figure 745744DEST_PATH_IMAGE021
和
Figure 528892DEST_PATH_IMAGE022
分别表示检波点位置处的观测地震数据与模拟地震数据，模拟地震数据
Figure 897557DEST_PATH_IMAGE022
来自于扰动波场
Figure 306410DEST_PATH_IMAGE023
，
Figure 76920DEST_PATH_IMAGE024
表示时间；当方程（6）的互相关值达到最大时，对应的
Figure 765391DEST_PATH_IMAGE025
等于
Figure 621351DEST_PATH_IMAGE018
；
根据方程（5）和（6），基于伴随状态法推导反射波旅行时反演的梯度方程，并将方程（2）和（4）代入梯度方程得方程（7）：
Figure 600940DEST_PATH_IMAGE026
Figure 960377DEST_PATH_IMAGE027
Figure 350907DEST_PATH_IMAGE028
（7）
其中，
Figure 428584DEST_PATH_IMAGE029
表示梯度，
Figure 444819DEST_PATH_IMAGE030
表示两个矢量的内积运算，
Figure 924342DEST_PATH_IMAGE031
和
Figure 220194DEST_PATH_IMAGE032
分别表示背景波场和扰动 波场的二阶时间导数，
Figure 50747DEST_PATH_IMAGE033
和
Figure 106559DEST_PATH_IMAGE034
分别表示伴随背景波场和伴随扰动波场，
Figure 440588DEST_PATH_IMAGE035
表示震源子波矢 量的二阶时间导数，
Figure 172921DEST_PATH_IMAGE036
表示零延迟互相关，
Figure DEST_PATH_IMAGE037
表示上行扰动波场的激发振幅；伴随背景波 场和伴随扰动波场分别通过求解方程（1）和（3）计算，其中，震源子波矢量
Figure 67934DEST_PATH_IMAGE009
替换为反射波旅 行时反演方法的伴随震源
Figure 786491DEST_PATH_IMAGE038
，见方程（8）；
Figure 99661DEST_PATH_IMAGE039
（8）
根据方程（7），给定合适的步长，得到本次迭代的速度更新量。
7.根据权利要求1所述的一种基于波场激发近似的反射波旅行时反演方法，其特征在于，步骤（5）中，保证每次迭代的速度更新量在20m/s-100m/s之间，迭代更新速度参数，直到满足收敛条件。
8.根据权利要求7所述的一种基于波场激发近似的反射波旅行时反演方法，其特征在于，收敛条件为：计算本次迭代模拟地震数据与观测地震数据旅行时误差的平方值，与上一次迭代的值进行比较，倘若数值减小，则继续进行步骤（1），当误差的平方值连续5次处于非下降状态，则认为反演方法已经收敛到全局最小值，满足收敛条件，此时，输出速度数据即为最终的反演结果。
CN202210381052.2A 2022-04-13 2022-04-13 一种基于波场激发近似的反射波旅行时反演方法 Active CN114460646B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210381052.2A CN114460646B (zh) 2022-04-13 2022-04-13 一种基于波场激发近似的反射波旅行时反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210381052.2A CN114460646B (zh) 2022-04-13 2022-04-13 一种基于波场激发近似的反射波旅行时反演方法

Publications (2)

Publication Number Publication Date
CN114460646A CN114460646A (zh) 2022-05-10
CN114460646B true CN114460646B (zh) 2022-06-28

Family

ID=81418599

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210381052.2A Active CN114460646B (zh) 2022-04-13 2022-04-13 一种基于波场激发近似的反射波旅行时反演方法

Country Status (1)

Country Link
CN (1) CN114460646B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115755178B (zh) * 2023-01-06 2023-04-25 青岛欧谱赛斯海洋科技有限公司 一种基于积分地震子波的时间域全波形反演方法
CN116660981B (zh) * 2023-07-25 2023-10-24 北京中矿大地地球探测工程技术有限公司 基于包络的各向异性参数反演方法、装置及存储介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2605047A1 (en) * 2011-12-15 2013-06-19 CGGVeritas Services SA Wave-fields separation for seismic recorders distributed at non-flat recording surfaces
EP2863243A2 (en) * 2011-11-01 2015-04-22 Geco Technology B.V. Methods and devices for transformation of collected seismic data for improved visualization capability
CN108037526A (zh) * 2017-11-23 2018-05-15 中国石油大学(华东) 基于全波波场vsp/rvsp地震资料的逆时偏移方法
CN108873066A (zh) * 2018-06-26 2018-11-23 中国石油大学(华东) 弹性介质波动方程反射波旅行时反演方法
CN110187382A (zh) * 2019-03-05 2019-08-30 中国石油大学(华东) 一种回折波和反射波波动方程旅行时反演方法
CN111751881A (zh) * 2019-03-29 2020-10-09 中国石油天然气集团有限公司 一种海上采集地震数据旅行时的校正方法、装置及系统
CN113534259A (zh) * 2021-07-09 2021-10-22 中石化石油工程技术服务有限公司 一种可控震源高效采集实时叠前时间偏移成像方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2755059B1 (en) * 2013-01-15 2021-11-03 CGG Services SAS Seismic data processing including data-constrained surface-consistent correction
CN104391323B (zh) * 2014-11-21 2015-11-18 中国石油大学(华东) 一种利用反射波信息反演速度场中低波数成分的方法
GB2538807B (en) * 2015-05-29 2019-05-15 Sub Salt Solutions Ltd Method for improved geophysical investigation

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2863243A2 (en) * 2011-11-01 2015-04-22 Geco Technology B.V. Methods and devices for transformation of collected seismic data for improved visualization capability
EP2605047A1 (en) * 2011-12-15 2013-06-19 CGGVeritas Services SA Wave-fields separation for seismic recorders distributed at non-flat recording surfaces
CN108037526A (zh) * 2017-11-23 2018-05-15 中国石油大学(华东) 基于全波波场vsp/rvsp地震资料的逆时偏移方法
CN108873066A (zh) * 2018-06-26 2018-11-23 中国石油大学(华东) 弹性介质波动方程反射波旅行时反演方法
CN110187382A (zh) * 2019-03-05 2019-08-30 中国石油大学(华东) 一种回折波和反射波波动方程旅行时反演方法
CN111751881A (zh) * 2019-03-29 2020-10-09 中国石油天然气集团有限公司 一种海上采集地震数据旅行时的校正方法、装置及系统
CN113534259A (zh) * 2021-07-09 2021-10-22 中石化石油工程技术服务有限公司 一种可控震源高效采集实时叠前时间偏移成像方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于波动方程重建震源子波的三维全波形反演;梁展源等;《石油地球物理勘探》;20171231;第52卷(第6期);第1200-1207页 *
近地表结构调查及参数反演综述;沈鸿雁等;《石油物探》;20190731;第58卷(第4期);第471-485页 *

Also Published As

Publication number Publication date
CN114460646A (zh) 2022-05-10

Similar Documents

Publication Publication Date Title
RU2693495C1 (ru) Полная инверсия волнового поля с компенсацией показателя качества
CN114460646B (zh) 一种基于波场激发近似的反射波旅行时反演方法
Kim et al. 3-D traveltime computation using second-order ENO scheme
CN108873066B (zh) 弹性介质波动方程反射波旅行时反演方法
CN113221393B (zh) 一种基于非结构有限元法的三维大地电磁各向异性反演方法
US8406081B2 (en) Seismic imaging systems and methods employing tomographic migration-velocity analysis using common angle image gathers
US20100054082A1 (en) Reverse-time depth migration with reduced memory requirements
US5394325A (en) Robust, efficient three-dimensional finite-difference traveltime calculations
CN106526674B (zh) 一种三维全波形反演能量加权梯度预处理方法
CN103713315B (zh) 一种地震各向异性参数全波形反演方法及装置
US10310117B2 (en) Efficient seismic attribute gather generation with data synthesis and expectation method
WO2020009752A1 (en) Full wavefield inversion with an image-gather-flatness constraint
CN113064203B (zh) 共轭梯度归一化lsrtm方法、系统、存储介质及应用
WO2016193180A1 (en) Improved method for inversion modelling
WO2016093875A1 (en) Systems and methods for aligning a monitor seismic survey with a baseline seismic survey
Datta et al. Full-waveform inversion of salt models using shape optimization and simulated annealing
US6324478B1 (en) Second-and higher-order traveltimes for seismic imaging
Hole et al. Interface inversion using broadside seismic refraction data and three‐dimensional travel time calculations
NO20190489A1 (en) Seismic modeling
CN113219531A (zh) 致密砂岩气水分布的识别方法及装置
CN115755178A (zh) 一种基于积分地震子波的时间域全波形反演方法
CN113031067B (zh) 一种基于Rytov-WKBJ近似的叠前地震反演方法
CN111175822B (zh) 改进直接包络反演与扰动分解的强散射介质反演方法
CN114460632A (zh) 快速射线追踪方法、电子设备及存储介质
Khoury et al. DSR migration velocity analysis by differential semblance optimization

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