CN111045077B - 一种陆地地震数据的全波形反演方法 - Google Patents

一种陆地地震数据的全波形反演方法 Download PDF

Info

Publication number
CN111045077B
CN111045077B CN201911327139.6A CN201911327139A CN111045077B CN 111045077 B CN111045077 B CN 111045077B CN 201911327139 A CN201911327139 A CN 201911327139A CN 111045077 B CN111045077 B CN 111045077B
Authority
CN
China
Prior art keywords
seismic
data
land
seismic data
parameter model
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
CN201911327139.6A
Other languages
English (en)
Other versions
CN111045077A (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.)
Beijing Research Institute of Uranium Geology
Original Assignee
Beijing Research Institute of Uranium Geology
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 Beijing Research Institute of Uranium Geology filed Critical Beijing Research Institute of Uranium Geology
Priority to CN201911327139.6A priority Critical patent/CN111045077B/zh
Publication of CN111045077A publication Critical patent/CN111045077A/zh
Application granted granted Critical
Publication of CN111045077B publication Critical patent/CN111045077B/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/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
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/624—Reservoir parameters

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

本发明属于一种陆地勘探地震数据反演领域，具体公开一种陆地地震数据的全波形反演方法，建立光滑初始化各向异性参数模型；从单炮地震记录中获取分频带地震观测记录；根据参数模型和地震观测系统，获取地震正演模拟合成数据集；计算合成地震数据与实际观测地震数据的残差；计算基于KR‑OT目标函数的伴随场震源；计算高斯加权函数；计算KR‑OT伴随场震源；设定基于伴随场震源的目标函数，并利代替步骤三中参数模型；重复上述步骤三至步骤八，直至数据残差满足精度要求，即为全波形反演获得的模型。本发明的方法利用KR‑OT目标函数，通过拟合体波和面波的地震数据，实现陆地地震勘探数据浅部和深部地质结构的精细反演。

Description

一种陆地地震数据的全波形反演方法
技术领域
本发明属于陆地勘探地震数据反演领域，具体涉及一种全波形目标函数的设计及反演方法。
背景技术
面波对于陆地数据全波形反演带来极大挑战。面波振幅强，但是穿透深度浅。而深部成像依靠的是穿透深度较深的体波(包括压缩波P波和剪切波S波)，面波的出现会掩盖体波(尤其是反射波部分)对于反演的贡献。目前流行的方法包括两种。第一，将观测数据中的面波滤除，使用声波方程来拟合观测到的体波部分。这种折衷方案减少了计算代价，但是只能在低频带适用。当频率增高之后，弹性效应明显，体波部分受到剪切波速度参数的影响也相应增大；此时用声波模拟体波并不合适。第二，放弃声波反演，应用弹性全波形反演，同时保留陆地数据中的面波。由于面波的频散特性，传统的2范数目标函数经常会出现跳周问题；目前流行的解决方法是使用其它目标函数拟合面波，比如包络线目标函数。虽然面波会被部分拟合，但面波残差依然是比体波部分大很多，所以体波携带的模型深部信息无法被提取。
发明内容
本发明的目的在于提供一种陆地地震数据的全波形反演方法，该方法利用KR-OT目标函数，通过同时拟合体波和面波的地震数据，实现陆地地震勘探数据浅部和深部地质结构的精细反演。
实现本发明目的的技术方案：
一种陆地地震数据的全波形反演方法，该方法具体包括以下步骤：
步骤一，建立陆地地震的光滑的初始化各向异性参数模型；
步骤二，从单炮地震记录中获取分频带的地震观测记录；
步骤三，根据参数模型和实际的地震观测系统，获取地震正演模拟合成数据集；
步骤四，计算合成地震数据与实际观测地震数据的残差δds,r；
步骤五，对于每一地震道，以体波第一个波包的中心作为高斯函数的中心点，计算高斯加权函数Ws,r(t)；
步骤六，利用上述步骤五求取的高斯加权函数Ws,r(t)，计算KR-OT伴随场震源Dd；
步骤七，设定基于KR-OT伴随场震源的目标函数
Figure BDA0002328663500000021
计算目标函数参数模型更新的方向和大小，并利用目标函数参数模型代替步骤三中使用的参数模型；
步骤八，重复上述步骤三至步骤七，直至最终的目标函数满足精度要求，此时的参数模型即为全波形反演在当前频带所获得的陆地地震数据的各项异性参数模型。
所述的步骤二的具体步骤如下：并从低频到高频，对观测数据进行带通滤波，获取分频带的地震观测记录。
所述的步骤三的具体步骤如下：从最低频带开始，采用交错网格，开展复杂地质情形下的地震波正演模拟计算，空间精度8阶，时间精度2阶，获取地震正演模拟合成数据集。
所述的步骤四中的合成地震数据与实际观测地震数据的残差
Figure BDA0002328663500000022
所述的步骤五中的高斯加权函数
Figure BDA0002328663500000031
所述的步骤六中的伴随场震源
Figure BDA0002328663500000032
所述的步骤七中的伴随场震源的目标函数
Figure BDA0002328663500000033
本发明的有益技术效果在于：本发明的陆地地震数据的全波形反演方法充分利用新目标函数的优点，同时增加高斯时间窗抑制远偏移距面波引起的跳周问题，最终充分发挥体波和面波对于反演模型的贡献，实现深部和浅部速度结构的同时重构。
对于面波问题，本发明的陆地地震数据的全波形反演方法利用面波，设计新的目标函数，构建新的反演流程，最终同时反演出模型深部和浅部的结构信息。新目标函数KR-OT(Kantorovich-Rubinstein Optimal Transport based objective function)是基于最优化传输的，具有两个优点。第一，比传统2范数更凸，对于全波形反演的初始模型精度的要求可以放低；第二，能自动平衡体波和面波的振幅。这两个特征使得KR-OT区别于其它目标函数。KR-OT目标函数的原理在于，它是一个拟合模式的函数，或者说是匹配观测数据与模拟数据相似性的函数。KR-OT可以沿着时间轴和空间轴匹配事件。由于地震事件在被比较的时候可以沿着时间轴和空间轴移动，所以无论是面波部分还是体波部分，其对伴随场的贡献都很平衡。在反演过程中，由于新目标函数自动平衡体波和面波的贡献，所以无论是浅部还是深部，地下速度结构都被重构。对于地震图拟合，虽然有部分数据因为高斯窗的使用而没有参与反演，但是最终的速度模型可以解释整张地震图上的各个震相，同时实现反射波到面波的拟合。
附图说明
图1为本发明所提供的模拟测试中的真实参数模型和初始化参数模型的对比示意图；
其中，图1(a-1)是真实的Vn速度模型的示意图；图1(a-2)是光滑的初始Vn速度模型的示意图；图1(b-1)是真实的Vs速度模型的示意图；图1(b-2)是光滑的初始Vs模型的示意图；图1(c-1)是真实的密度模型的示意图；图1(c-2)光滑的初始密度模型的示意图，根据初始Vn计算出来；图1(d)是真实各向异性参数η的示意图；图1(e)是真实各向异性参数δ的示意图，初始的两个各向异性参数为零；
图2为本发明所使用的震源子波的示意图；
图3为本发明所提供的地震数据与实际观测地震数据的残差的示意图；
其中，图3(a)为输入的数据残差的示意图，即2范数目标函数下的伴随场震源；图3(b)为高斯时间窗函数的示意图；图3(c)为基于高斯窗函数的KR-OT伴随场震源的示意图；
图4为本发明所提供的全波形反演的最终参数模型结果的示意图；
其中，图4(a)为纵波速度场Vn的示意图；图4(b)为各项异性参数η的示意图；图4(c)为各向异性参数δ的示意图；图4(d)为横波速度场Vs的示意图；图4(e)为介质密度的示意图；
图5为本发明所提供的位于不同空间位置处的两道检波器数据对比示意图。
具体实施方式
下面结合附图和实施例对本发明作进一步详细说明。
如图1至图5所示，本发明所提供的一种陆地地震数据的全波形反演方法，该方法具体包括以下步骤：
步骤一，建立陆地地震的光滑的初始化各向异性参数模型。
图1为模拟测试中的真实参数模型和初始化参数模型的对比示意图。真实模型是各项异性的VTI模型，初始模型是光滑的各项同性模型。图1(a-1)是真实的Vn速度模型，图1(a-2)是光滑的初始Vn速度模型；图1(b-1)是真实的Vs速度模型，图1(b-2)是光滑的初始Vs模型；图1(c-1)是真实的密度模型，图1(c-2)是光滑的初始密度模型，根据初始Vn计算出来；图1(d)和图1(e)分别是真实的各项异性参数eta和delta。
根据已有的地质、钻井、测井和地震资料，建立陆地地震的光滑的初始化各向异性参数模型。上述光滑的初始化各向异性参数模型为光滑的初始Vn速度模型、光滑的初始Vs模型、光滑的初始密度模型。
步骤二，从单炮地震记录中获取分频带的地震观测记录，具体步骤如下：
从低频到高频，对观测数据进行带通滤波，获取分频带的地震观测记录。图2为地震记录中的震源子波-Ricker子波示意图，其中心频率为3赫兹。
步骤三，根据参数模型和实际的地震观测系统，获取地震正演模拟合成数据集，具体步骤如下：
从最低频带开始，根据参数模型和实际的地震观测系统，采用交错网格，开展复杂地质情形下的地震波正演模拟计算，空间精度8阶，时间精度2阶，获取地震正演模拟合成数据集。
步骤四，计算合成地震数据与实际观测地震数据的残差δds,r。其计算公式如下：
Figure BDA0002328663500000061
其中，δds,r代表地震观测数据集与合成数据集的残差，
Figure BDA0002328663500000062
代表地震观测数据集，
Figure BDA0002328663500000063
代表地震合成数据集，下标s代表震源角标，下标r代表检波器角标。图3(a)所示即为单炮地震记录的数据残差。
步骤五，对于每一地震道，以体波第一个波包的中心作为高斯函数的中心点，计算高斯加权函数Ws,r(t)，其具体公式为：
Figure BDA0002328663500000064
其中，t代表时间变量，
Figure BDA0002328663500000065
代表体波第一个波包的中心位置，σ代表方差，Ws,r(t)代表求取的高斯加权函数。
图3(b)所示为单个单炮记录求取的高斯加权函数。
步骤六，基于上述步骤四求取的数据残差δds,r，计算基于KR-OT目标函数的优化函数
Figure BDA0002328663500000066
其中，xr代表检波器的空间位置坐标，t代表时间变量，
Figure BDA0002328663500000067
为Blip1，1-Lipschitz函数，其公式如下：
Figure BDA0002328663500000068
Figure BDA0002328663500000069
Figure BDA00023286635000000610
其中，δxr代表空间方向上的任意增量，δt代表时间方向上的任意增量。采用SDMM(Simultaneous Direction method of Multipliers)的方法，求解出最优的
Figure BDA00023286635000000611
进而获取基于KR-OT目标函数的伴随场震源
Figure BDA00023286635000000612
图3(c)所示为远偏移距面波得到压制后的KR-OT伴随场震源。
步骤七，设定基于优化的KR-OT伴随场震源的目标函数
Figure BDA00023286635000000613
计算目标函数参数模型更新的方向和大小，并利用新计算的目标函数参数模型代替步骤三中使用的参数模型。
基于优化的KR-OT伴随场震源的目标函数
Figure BDA0002328663500000071
的具体公式为：
Figure BDA0002328663500000072
其中，max代表求取其中最大值，Δxr代表检波器间距。
基于以上目标函数，采用准牛顿迭代法，计算参数模型更新的方向和大小，并利用新计算的参数模型代替步骤三中使用的参数模型。
步骤八，重复上述步骤三至步骤八，直至最终的目标函数满足精度要求，即最终的目标函数是初始目标函数的百分之十以内精度，此时的参数模型即为全波形反演在当前频带所获得的陆地地震数据的各项异性参数模型。
步骤九，移动到下个地震频带，并把上一个地震频带的最终速度模型当成下个地震频带的初始模型，重复步骤八，得到最终陆地地震数据的全波形反演模型。
图4所示即为全波形反演的最终结果。其中图4(a)为纵波速度场Vn；图4(b)为各项异性参数η；图4(c)为各向异性参数δ；图4(d)为横波速度场Vs；图4(e)为介质密度。可见无论浅层还是深层，模型参数都得到很好恢复。
图5所示为不同空间位置处的两道检波器的数据对比。从中可见，在初始数据存在跳周情况下，本发明的全波形反演方法，仍然能够很好地拟合观测数据。
上面结合附图和实施例对本发明作了详细说明，但是本发明并不限于上述实施例，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。本发明中未作详细描述的内容均可以采用现有技术。

Claims (6)

1.一种陆地地震数据的全波形反演方法，其特征在于，该方法具体包括以下步骤：
步骤一，建立陆地地震的光滑的初始化各向异性参数模型；
步骤二，从单炮地震记录中获取分频带的地震观测记录；
步骤三，根据参数模型和实际的地震观测系统，获取地震正演模拟合成数据集；
步骤四，计算合成地震数据与实际观测地震数据的残差δds,r；
步骤五，对于每一地震道，以体波第一个波包的中心作为高斯函数的中心点，计算高斯加权函数Ws,r(t)；其中，t代表时间变量；
步骤六，基于上述步骤四求取的数据残差δds,r，计算基于KR-OT目标函数的优化函数
Figure FDA0003339774010000015
其中，xr代表检波器的空间位置坐标；
步骤七，利用上述步骤五求取的高斯加权函数Ws,r(t)，计算优化后的KR-OT伴随场震源Dd；
步骤八，设定基于优化的KR-OT伴随场震源的目标函数
Figure FDA0003339774010000011
计算目标函数参数模型更新的方向和大小，并利用目标函数参数模型代替步骤三中使用的参数模型；所述的步骤八中的基于优化的KR-OT伴随场震源的目标函数
Figure FDA0003339774010000012
其中，max代表求取其中最大值，Δxr代表检波器间距，
Figure FDA0003339774010000013
代表地震观测数据集，
Figure FDA0003339774010000014
代表地震合成数据集，下标s代表震源角标，下标r代表检波器角标；
步骤九，重复上述步骤三至步骤八，直至最终的数据残差满足精度要求，此时的参数模型即为全波形反演在当前频带所获得的陆地地震数据的各项异性参数模型。
2.根据权利要求1所述的一种陆地地震数据的全波形反演方法，其特征在于：所述的步骤二的具体步骤如下：从单炮地震记录中提取每炮记录的震源子波，并从低频到高频，对观测数据进行带通滤波，获取分频带的地震观测记录。
3.根据权利要求2所述的一种陆地地震数据的全波形反演方法，其特征在于：所述的步骤三的具体步骤如下：从最低频带开始，采用交错网格，开展复杂地质情形下的地震波正演模拟计算，空间精度8阶，时间精度2阶，获取地震正演模拟合成数据集。
4.根据权利要求3所述的一种陆地地震数据的全波形反演方法，其特征在于：所述的步骤四中的合成地震数据与实际观测地震数据的残差
Figure FDA0003339774010000021
5.根据权利要求4所述的一种陆地地震数据的全波形反演方法，其特征在于：所述的步骤五中的高斯加权函数
Figure FDA0003339774010000022
其中，
Figure FDA0003339774010000023
代表体波第一个波包的中心位置，σ代表方差。
6.根据权利要求5所述的一种陆地地震数据的全波形反演方法，其特征在于：所述的步骤七中的伴随场震源
Figure FDA0003339774010000024
其中，
Figure FDA0003339774010000025
为Blip1，1-Lipschitz函数。
CN201911327139.6A 2019-12-20 2019-12-20 一种陆地地震数据的全波形反演方法 Active CN111045077B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911327139.6A CN111045077B (zh) 2019-12-20 2019-12-20 一种陆地地震数据的全波形反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911327139.6A CN111045077B (zh) 2019-12-20 2019-12-20 一种陆地地震数据的全波形反演方法

Publications (2)

Publication Number Publication Date
CN111045077A CN111045077A (zh) 2020-04-21
CN111045077B true CN111045077B (zh) 2022-03-18

Family

ID=70238076

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911327139.6A Active CN111045077B (zh) 2019-12-20 2019-12-20 一种陆地地震数据的全波形反演方法

Country Status (1)

Country Link
CN (1) CN111045077B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114442169B (zh) * 2020-11-03 2024-05-28 中国石油天然气股份有限公司 地震资料中近源信号压制方法及装置
CN117687083A (zh) * 2022-09-09 2024-03-12 中国石油化工股份有限公司 全波形反演方法、设备及存储介质
CN116626751B (zh) * 2023-05-26 2024-03-19 北京中矿大地地球探测工程技术有限公司 基于多目标函数的黏弹性参数同步反演方法、装置和设备
CN116660981B (zh) * 2023-07-25 2023-10-24 北京中矿大地地球探测工程技术有限公司 基于包络的各向异性参数反演方法、装置及存储介质

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012141799A2 (en) * 2011-02-25 2012-10-18 University Of Florida Research Foundation, Inc. Detection of sinkholes or anomalies
CN103913772A (zh) * 2014-04-02 2014-07-09 西南石油大学 基于储层地质力学参数的微地震事件正演模拟方法
CN105467444A (zh) * 2015-12-10 2016-04-06 中国石油天然气集团公司 一种弹性波全波形反演方法及装置
CN106230028A (zh) * 2016-09-08 2016-12-14 安徽电气工程职业技术学院 一种风电—抽水蓄能联合系统的多目标优化方法
CN107422379A (zh) * 2017-07-27 2017-12-01 中国海洋石油总公司 基于局部自适应凸化方法的多尺度地震全波形反演方法
CN109975869A (zh) * 2019-03-27 2019-07-05 中国石油大学(北京) 一种沿地层走向光滑约束的反射波波形反演方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150362623A1 (en) * 2014-06-12 2015-12-17 Westerngeco, Llc Joint inversion of attributes
US11360223B2 (en) * 2017-09-21 2022-06-14 Chevron U.S.A. Inc. System and method for improved full waveform inversion

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012141799A2 (en) * 2011-02-25 2012-10-18 University Of Florida Research Foundation, Inc. Detection of sinkholes or anomalies
CN103913772A (zh) * 2014-04-02 2014-07-09 西南石油大学 基于储层地质力学参数的微地震事件正演模拟方法
CN105467444A (zh) * 2015-12-10 2016-04-06 中国石油天然气集团公司 一种弹性波全波形反演方法及装置
CN106230028A (zh) * 2016-09-08 2016-12-14 安徽电气工程职业技术学院 一种风电—抽水蓄能联合系统的多目标优化方法
CN107422379A (zh) * 2017-07-27 2017-12-01 中国海洋石油总公司 基于局部自适应凸化方法的多尺度地震全波形反演方法
CN109975869A (zh) * 2019-03-27 2019-07-05 中国石油大学(北京) 一种沿地层走向光滑约束的反射波波形反演方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"互相关与最小二乘加权目标函数全波形反演";梁煌 等;《世界地质》;20170630;第36卷(第2期);第588-594页 *

Also Published As

Publication number Publication date
CN111045077A (zh) 2020-04-21

Similar Documents

Publication Publication Date Title
CN111045077B (zh) 一种陆地地震数据的全波形反演方法
CN108345031B (zh) 弹性介质主动源和被动源混采地震数据全波形反演方法
CN109669212B (zh) 地震数据处理方法、地层品质因子估算方法与装置
KR101548976B1 (ko) 지진 표면파들의 파형들을 사용하는 토양 특성들의 추정
US10620331B2 (en) Reverse time migration in anisotropic media with stable attenuation compensation
CN103492910B (zh) 时域中的同步小波提取和反卷积
CN105467444B (zh) 一种弹性波全波形反演方法及装置
CN103713315B (zh) 一种地震各向异性参数全波形反演方法及装置
CN107505654A (zh) 基于地震记录积分的全波形反演方法
CN104237945B (zh) 一种地震资料自适应高分辨处理方法
CN111025387B (zh) 一种页岩储层的叠前地震多参数反演方法
CN110187382B (zh) 一种回折波和反射波波动方程旅行时反演方法
CN104570116A (zh) 基于地质标志层的时差分析校正方法
CN109655918B (zh) 地面浅井微地震监测观测台站位置确定方法及系统
CN110737018B (zh) Vsp地震资料各向异性建模方法
Sun et al. Joint minimization of the mean and information entropy of the matching filter distribution for a robust misfit function in full-waveform inversion
CN110780341B (zh) 一种各向异性地震成像方法
CN111158047B (zh) 一种三维弹性波场矢量分解法、装置及计算机存储介质
CN109143345B (zh) 基于模拟退火的品质因子q非线性反演方法及系统
CN110873895A (zh) 一种变网格微地震逆时干涉定位方法
CN111175822B (zh) 改进直接包络反演与扰动分解的强散射介质反演方法
CN108680957B (zh) 基于加权的局部互相关时频域相位反演方法
Kurzmann et al. Real data applications of seismic full waveform inversion
Li et al. Active-source Rayleigh wave dispersion by the Aki spectral formulation
Kazei et al. Acquisition and near-surface impacts on VSP mini-batch FWI and RTM imaging in desert environment

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