CN113296146A - 一种基于梯度道集相关加权的全波形反演梯度预处理方法 - Google Patents

一种基于梯度道集相关加权的全波形反演梯度预处理方法 Download PDF

Info

Publication number
CN113296146A
CN113296146A CN202110560874.2A CN202110560874A CN113296146A CN 113296146 A CN113296146 A CN 113296146A CN 202110560874 A CN202110560874 A CN 202110560874A CN 113296146 A CN113296146 A CN 113296146A
Authority
CN
China
Prior art keywords
gradient
formula
gather
waveform inversion
profile
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.)
Granted
Application number
CN202110560874.2A
Other languages
English (en)
Other versions
CN113296146B (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.)
Ocean University of China
Original Assignee
Ocean University of China
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 Ocean University of China filed Critical Ocean University of China
Priority to CN202110560874.2A priority Critical patent/CN113296146B/zh
Publication of CN113296146A publication Critical patent/CN113296146A/zh
Application granted granted Critical
Publication of CN113296146B publication Critical patent/CN113296146B/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

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

Abstract

本发明涉及一种基于梯度道集相关加权的全波形反演梯度预处理方法，属于地震偏移成像领域，具体包括以下步骤：首先在全波形反演中得到每一炮的梯度剖面，在各炮的梯度剖面中将同一位置处的梯度道分别抽取出来，并按照道与炮之间的距离由小到大进行排列，以此得到各个位置处的梯度道集，然后在每个梯度道集中，对远道的噪音部分进行初步切除，并将梯度道集中的所有梯度道与参考道逐一进行相关，选取相关值较大的N道进行叠加，从而形成一个新的梯度道，对各个梯度道集进行上述处理，并将各个新的梯度道放于对应位置处，得到最终的梯度剖面。本发明方法可提高梯度预处理的精度，减少迭代次数，节约计算时间，从而显著提升全波形反演的效果。

Description

一种基于梯度道集相关加权的全波形反演梯度预处理方法
技术领域
本发明属于地震速度建模领域，具体涉及一种基于梯度道集相关加权的全波形反演梯度预处理方法。
背景技术
全波形反演作为一种高精度的地震速度建模算法，其每次迭代过程中，首先基于初始速度模型，采用双程波方程有限差分数值模拟技术得到合成地震记录，然后对合成地震记录与实测地震记录的差值建立目标泛函，基于最小平方框架下的优化思想计算出梯度方向，从而得到模型修改量完成一次迭代。整个全波形反演通常需要几十次甚至上百次的迭代反演。
全波形反演不仅考虑了地震波传播的旅行时等运动学信息，而且充分利用了振幅、相位等动力学信息，理想条件下其反演精度可达到波长数量级，因此其在未来的油气勘探中被寄予厚望。
常规全波形反演算法求取梯度剖面时，总的梯度剖面是由每一炮的梯度剖面叠加而成，实际上，当梯度道和炮点距离很远时，当前炮对此道的梯度几乎无贡献，相反，会带来不规则干扰，影响梯度剖面的质量，进而增加迭代次数和计算时间。当前并无很好地选取有效梯度炮从而提升反演收敛效率的方法。
发明内容
本发明要解决的技术问题在于提供一种基于梯度道集相关加权的全波形反演梯度预处理方法。首先在全波形反演中得到每一炮的梯度剖面，在各炮的梯度剖面中将同一位置处的梯度道分别抽取出来，并按照道与炮之间的距离由小到大进行排列，以此得到各个位置处的梯度道集，然后在每个梯度道集中，对远道的噪音部分进行初步切除，并将梯度道集中的所有梯度道与参考道(即每个处理后的梯度道集所有道叠加的结果)逐一进行相关，选取相关值较大的N道(N为给定的有效叠加道数，一般取该条地震测线满覆盖次数的1/2)进行叠加，从而形成一个新的梯度道，对各个梯度道集进行上述处理，并将各个新的梯度道放于对应位置处，得到最终的梯度剖面。模型实验结果表明，本方法可提高梯度预处理的精度，减少迭代次数，节约计算时间，从而显著提升全波形反演的效果。
本发明采取以下技术方案：
一种基于梯度道集相关加权的全波形反演梯度预处理方法，其特征在于所述方法具体包括以下步骤：
(1)给定一个初始速度模型V(x,z)，其中x、z表示空间位置坐标(x＝1,2,3,…,Nx,z＝1,2,3,…,Nz)，Nx、Nz分别代表模型横向和纵向网格点总数；给定地震子波Wi(t)，i表示炮的序号(i＝1,2,3,…,S)，S表示总炮数，t表示时间，基于初始速度模型，应用声波方程通过有限差分进行正演模拟，得到各炮的正时波场Ui(x,z,t)和合成地震记录Cali(x,z,t),应用声波方程进行逆时延拓得到各炮的反传波场Uni(x,z,t)；
(2)给定实际地震记录Obsi(x,z,t)，将其与合成地震记录的差值作为逆时扰动进行逆时延拓得到各炮的逆时波场Ri(x,z,t)；
(3)利用公式(1)得到各炮的常规预处理梯度剖面Gi(x,z)，其中
Figure BDA0003073780080000021
表示对Ui(x,z,t)的时间微分，
Figure BDA0003073780080000022
表示对Ri(x,z,t)的时间微分，这样总共可以得到S个梯度剖面；所述的公式(1)为：
Figure BDA0003073780080000031
(4)对于这S个梯度剖面，将每一个剖面中的同一个x位置处的道提取出来，形成Nx个梯度道集，记为IGx(i,z)，其中x＝1,2,3,…,Nx,i＝1,2,3,…,S，z＝1,2,3,…,Nz，每个梯度道集中，按照炮点位置与当前x的距离由小到大排列；
(5)在各个梯度道集IGx(i,z)中，对远道的噪音部分进行初步切除，然后利用公式(2)将所有道进行叠加，形成参考道Bx(z)；所述的公式(2)为：
Figure BDA0003073780080000032
(6)利用公式(3)将各个梯度道集IGx(i,z)中的各道分别与该梯度道集所对应的参考道Bx(z)做相关，得到各道相对于参考道的相关值Cx(i)；所述的公式(3)为：
Figure BDA0003073780080000033
(7)在各个梯度道集IGx(i,z)中取相关值最大的前N道，利用公式(4)进行叠加，形成新的梯度剖面Gnew(x,z)，其中N为给定的有效叠加道数，所述公式(4)为：
Figure BDA0003073780080000034
全波形反演每次迭代中的梯度预处理均需要重复以上7个步骤。
进一步，所述步骤(7)中的N取该条地震测线满覆盖次数的1/2。
本发明与现有技术相比的有益效果：
本发明提出的一种基于梯度道集相关加权的全波形反演梯度预处理方法。首先在全波形反演中得到每一炮的梯度剖面，在各炮的梯度剖面中将同一位置处的梯度道分别抽取出来，并按照道与炮之间的距离由小到大进行排列，以此得到各个位置处的梯度道集，然后在每个梯度道集中，对远道的噪音部分进行初步切除，并将梯度道集中的所有梯度道与参考道(即每个处理后的梯度道集所有道叠加的结果)逐一进行相关，选取相关值较大的N道(N为给定的有效叠加道数，一般取该条地震测线满覆盖次数的1/2)进行叠加，从而形成一个新的梯度道，对各个梯度道集进行上述处理，并将各个新的梯度道放于对应位置处，得到最终的梯度剖面。模型实验结果表明，本方法可提高梯度预处理的精度，减少迭代次数，节约计算时间，从而显著提升全波形反演的效果。
附图说明
图1为全波形反演常规梯度预处理方法流程图；
图2为基于梯度道集相关加权的全波形反演梯度预处理方法流程图；
图3为Marmousi速度模型；
图4为第250道的梯度道集图；
图5为对第250道的梯度道集切除后的图片；
图6为第250道梯度道集的参考道；
图7为常规梯度预处理后的梯度剖面；
图8为基于道集相关加权的梯度预处理后的梯度剖面；
图9常规梯度预处理的第15次反演结果；
图10基于梯度道集相关加权的梯度预处理的第10次反演结果。
具体实施方式
本发明以Marmousi模型为例(如图3所示)进行阐述，Marmousi模型横向长度为9220m，纵向长度为3460m,模型x、z方向的网格步长均为20m。
基于该模型，一共进行461次放炮，炮间隔为20m，设置461个检波点，检波点以20m为间隔均匀的分布在模型的正上方，炮点深度和检波点深度均为0m。
下面详细阐述本发明具体实施方式，流程如图2所示，具体步骤如下：
(1)将横向长度为9220m，纵向长度为3460m的Marmousi模型网格化，模型x、z方向的网格步长均为20m，得到网格大小为461×173的模型，其中x、z表示网格位置坐标(x＝1,2,3,…,461,z＝1,2,3,…,173)；给定地震子波Wi(t)，进行461次放炮，i表示炮的序号(i＝1,2,3,…,461)，t表示时间，基于Marmousi速度模型，应用声波方程通过有限差分进行正演模拟，得到各炮的正时波场Ui(x,z,t)和合成地震记录Cali(x,z,t)，应用声波方程进行逆时延拓得到各炮的反传波场Uni(x,z,t)；
(2)给定实际地震记录Obsi(x,z,t)，将其与合成地震记录的差值作为逆时扰动进行逆时延拓得到各炮的逆时波场Ri(x,z,t)；
(3)利用公式(1)得到各炮的常规预处理梯度剖面Gi(x,z)，其中
Figure BDA0003073780080000051
表示对Ui(x,z,t)的时间微分，
Figure BDA0003073780080000052
表示对Ri(x,z,t)的时间微分，这样总共可以得到461个梯度剖面；所述的公式(1)为
Figure BDA0003073780080000053
对于这461个梯度剖面，将每一个剖面中的同一个x位置处的道提取出来，形成461个梯度道集，记为IGx(i,z)，其中x＝1,2,3,…,461,i＝1,2,3,…,461，z＝1,2,3,…,173，每个梯度道集中，按照炮点位置与当前x的距离由小到大排列；图4为第250道的梯度道集图；
(4)在各个梯度道集IGx(i,z)中，对远道的噪音部分进行初步切除，图5为对第250道的梯度道集切除后的图片；然后利用公式(2)将所有道进行叠加，形成参考道Bx(z)；所述的公式(2)为：
Figure BDA0003073780080000054
图6为第250道梯度道集的参考道；
(5)利用公式(3)将各个梯度道集IGx(i,z)中的各道分别与该梯度道集所对应的参考道Bx(z)做相关，得到各道相对于参考道的相关值Cx(i)；所述的公式(3)为：
Figure BDA0003073780080000061
(6)在各个梯度道集IGx(i,z)中取相关值最大的前250道利用公式(4)进行叠加，形成新的梯度剖面Gnew(x,z)，所述公式(4)为：
Figure BDA0003073780080000062
为了证明本发明所述方法的有益效果，这里与常规梯度预处理(流程如图1所示)后的梯度剖面进行对比。如图7为常规梯度预处理后的梯度剖面，图8为基于梯度道集相关加权的梯度预处理梯度剖面。对比两图可以发现，图8顶部低频噪音得到消除，其两侧的同相轴连续性得到增强，同相轴均衡性变好。图9为常规梯度预处理的第15次反演结果，图10为基于梯度道集相关加权的梯度预处理的第10次反演结果，对比发现，两图精度相当，则说明基于道集相关加权的梯度预处理方法可以减少全波形反演的迭代次数，从而减少计算时间。

Claims (2)

1.一种基于梯度道集相关加权的全波形反演梯度预处理方法，其特征在于所述方法具体包括以下步骤：
(1)给定一个初始速度模型V(x,z)，其中x、z表示空间位置坐标，x＝1,2,3,…,Nx,z＝1,2,3,…,Nz，Nx、Nz分别代表模型横向和纵向网格点总数；给定地震子波Wi(t)，i表示炮的序号，i＝1,2,3,…,S，S表示总炮数，t表示时间，基于初始速度模型，应用声波方程通过有限差分进行正演模拟，得到各炮的正时波场Ui(x,z,t)和合成地震记录Cali(x,z,t),应用声波方程进行逆时延拓得到各炮的反传波场Uni(x,z,t)；
(2)给定实际地震记录Obsi(x,z,t)，将其与合成地震记录的差值作为逆时扰动进行逆时延拓得到各炮的逆时波场Ri(x,z,t)；
(3)利用公式(1)得到各炮的常规预处理梯度剖面Gi(x,z)，其中
Figure FDA0003073780070000011
表示对Ui(x,z,t)的时间微分，
Figure FDA0003073780070000012
表示对Ri(x,z,t)的时间微分，这样总共可以得到S个梯度剖面；所述的公式(1)为：
Figure FDA0003073780070000013
(4)对于这S个梯度剖面，将每一个剖面中的同一个x位置处的道提取出来，形成Nx个梯度道集，记为IGx(i,z)，其中x＝1,2,3,…,Nx,i＝1,2,3,…,S，每个梯度道集中，按照炮点位置与当前x的距离由小到大排列；
(5)在各个梯度道集IGx(i,z)中，对远道的噪音部分进行初步切除，然后利用公式(2)将所有道进行叠加，形成参考道Bx(z)；所述的公式(2)为：
Figure FDA0003073780070000014
(6)利用公式(3)将各个梯度道集IGx(i,z)中的各道分别与该梯度道集所对应的参考道Bx(z)做相关，得到各道相对于参考道的相关值Cx(i)；所述的公式(3)为：
Figure FDA0003073780070000021
(7)在各个梯度道集IGx(i,z)中取相关值最大的前N道利用公式(4)进行叠加，形成新的梯度剖面Gnew(x,z)，其中N为给定的有效叠加道数，所述公式(4)为：
Figure FDA0003073780070000022
全波形反演每次迭代中的梯度预处理均需要重复以上7个步骤。
2.根据权利要求1所述的一种基于梯度道集相关加权的全波形反演梯度预处理方法，其特征在于N取该条地震测线满覆盖次数的1/2。
CN202110560874.2A 2021-05-19 2021-05-19 一种基于梯度道集相关加权的全波形反演梯度预处理方法 Active CN113296146B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110560874.2A CN113296146B (zh) 2021-05-19 2021-05-19 一种基于梯度道集相关加权的全波形反演梯度预处理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110560874.2A CN113296146B (zh) 2021-05-19 2021-05-19 一种基于梯度道集相关加权的全波形反演梯度预处理方法

Publications (2)

Publication Number Publication Date
CN113296146A true CN113296146A (zh) 2021-08-24
CN113296146B CN113296146B (zh) 2022-02-22

Family

ID=77323919

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110560874.2A Active CN113296146B (zh) 2021-05-19 2021-05-19 一种基于梯度道集相关加权的全波形反演梯度预处理方法

Country Status (1)

Country Link
CN (1) CN113296146B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113449440A (zh) * 2021-08-30 2021-09-28 中国海洋大学 一种最小二乘逆时偏移的梯度道集相关加权预处理方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106526674A (zh) * 2016-11-14 2017-03-22 中国石油化工股份有限公司 一种三维全波形反演能量加权梯度预处理方法
CN107894618A (zh) * 2017-11-10 2018-04-10 中国海洋大学 一种基于模型平滑算法的全波形反演梯度预处理方法
US20180156933A1 (en) * 2016-12-02 2018-06-07 Bp Corporation North America Inc. Seismic acquisition geometry full-waveform inversion
CN110058302A (zh) * 2019-05-05 2019-07-26 四川省地质工程勘察院 一种基于预条件共轭梯度加速算法的全波形反演方法
WO2020009752A1 (en) * 2018-07-02 2020-01-09 Exxonmobil Upstream Research Company Full wavefield inversion with an image-gather-flatness constraint
CN111766628A (zh) * 2020-07-29 2020-10-13 浪潮云信息技术股份公司 一种预条件的时间域弹性介质多参数全波形反演方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106526674A (zh) * 2016-11-14 2017-03-22 中国石油化工股份有限公司 一种三维全波形反演能量加权梯度预处理方法
US20180156933A1 (en) * 2016-12-02 2018-06-07 Bp Corporation North America Inc. Seismic acquisition geometry full-waveform inversion
CN107894618A (zh) * 2017-11-10 2018-04-10 中国海洋大学 一种基于模型平滑算法的全波形反演梯度预处理方法
WO2020009752A1 (en) * 2018-07-02 2020-01-09 Exxonmobil Upstream Research Company Full wavefield inversion with an image-gather-flatness constraint
CN110058302A (zh) * 2019-05-05 2019-07-26 四川省地质工程勘察院 一种基于预条件共轭梯度加速算法的全波形反演方法
CN111766628A (zh) * 2020-07-29 2020-10-13 浪潮云信息技术股份公司 一种预条件的时间域弹性介质多参数全波形反演方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
王庆等: "时间域地震全波形反演方法进展", 《地球物理学进展》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113449440A (zh) * 2021-08-30 2021-09-28 中国海洋大学 一种最小二乘逆时偏移的梯度道集相关加权预处理方法

Also Published As

Publication number Publication date
CN113296146B (zh) 2022-02-22

Similar Documents

Publication Publication Date Title
CN109669212B (zh) 地震数据处理方法、地层品质因子估算方法与装置
CN105388520B (zh) 一种地震资料叠前逆时偏移成像方法
CN105093301B (zh) 共成像点反射角角道集的生成方法及装置
CN109343118B (zh) 一种异常初至时间修正方法
CN103499836B (zh) 空变-多时窗融合高精度剩余静校正方法
CN110187382B (zh) 一种回折波和反射波波动方程旅行时反演方法
CN111158049A (zh) 一种基于散射积分法的地震逆时偏移成像方法
CN111045077B (zh) 一种陆地地震数据的全波形反演方法
CN110389377B (zh) 基于波形互相关系数相乘的微震偏移成像定位方法
CN113296146B (zh) 一种基于梯度道集相关加权的全波形反演梯度预处理方法
CN113031063B (zh) 一种基于成像道集相关加权的逆时偏移成像方法
CN107894618A (zh) 一种基于模型平滑算法的全波形反演梯度预处理方法
EP2321671A2 (en) Processing seismic data in common group-center gathers
CN102590857A (zh) 真地表起伏叠前深度域双程波成像方法
CN102162858A (zh) 利用非对称走时进行动校正速度分析的方法
CN106257309A (zh) 叠后地震数据体处理方法及装置
CN111999764A (zh) 基于时频域目标函数的盐下构造最小二乘逆时偏移方法
CN112630830A (zh) 一种基于高斯加权的反射波全波形反演方法及系统
CN112946742B (zh) 一种拾取精确叠加速度谱的方法
CN113449440A (zh) 一种最小二乘逆时偏移的梯度道集相关加权预处理方法
CN110888158B (zh) 一种基于rtm约束的全波形反演方法
CN113325467A (zh) 一种基于槽波频散特征的微地震震源定位方法
CN112327356A (zh) 基于同相轴迭代追踪提取的混叠记录分离方法
CN1215424C (zh) 消除地震勘探记录中强相关干扰波的处理方法
CN114924313B (zh) 基于行波分离的声波最小二乘逆时偏移梯度求取方法

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