CN107356972A - 一种各向异性介质的成像方法 - Google Patents

一种各向异性介质的成像方法 Download PDF

Info

Publication number
CN107356972A
CN107356972A CN201710505548.5A CN201710505548A CN107356972A CN 107356972 A CN107356972 A CN 107356972A CN 201710505548 A CN201710505548 A CN 201710505548A CN 107356972 A CN107356972 A CN 107356972A
Authority
CN
China
Prior art keywords
data
continuation
residual error
field
wave
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
CN201710505548.5A
Other languages
English (en)
Other versions
CN107356972B (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 University of Petroleum East China
Original Assignee
China University of Petroleum East 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 China University of Petroleum East China filed Critical China University of Petroleum East China
Priority to CN201710505548.5A priority Critical patent/CN107356972B/zh
Publication of CN107356972A publication Critical patent/CN107356972A/zh
Application granted granted Critical
Publication of CN107356972B publication Critical patent/CN107356972B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/34Displaying seismic recordings or visualisation of seismic data or attributes
    • G01V1/345Visualisation of seismic data or attributes, e.g. in 3D cubes
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/50Corrections or adjustments related to wave propagation
    • G01V2210/58Media-related
    • G01V2210/586Anisotropic media

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

本发明公开了一种各向异性介质的成像方法,包括如下步骤:S1,获得观测数据,根据观测数据获得速度场和各相异性参数;S2,对观测数据进行偏移得到反射系数;S3,根据速度场和各相异性参数,对反射系数进行反偏移得到模拟数据;S4,计算观测数据和模拟数据的数据残差;S5,根据数据残差获得最终成像结果。通过采用本发明提供的方法,充分考虑了地下介质的各向异性这一实际情况,能够更加精确的模拟地震波在真实地下介质中的传播规律,进一步提高偏移成像结果的分辨率和保真度,使得偏移成像结果中不同深度的同相轴振幅更加均衡,成像结果更为准确。

Description

一种各向异性介质的成像方法
技术领域
本发明涉及一种各向异性介质的成像方法,属于地球物理勘探领域。
背景技术
各向异性介质即物理性质具有方向特性的地球物理介质。通常,地下沉积地层大多成层 分布,呈现横向各相同性(VTI,vertical transverse isotropic)介质特征。随着高精度地震 勘探技术的发展,能够更加真实且全面的反映地震波在地下传播规律的各向异性介质成像方 法越来越引人关注。当前的地下介质反演过程中,通常使用迭代算法,以获得成像结果。
但是,现有技术中的缺陷在于,通常的反演过程中,大多是基地下介质的各向同性假设, 忽略介质的各向异性,导致模拟波场和真实波场在传播过程中出现偏差,降低成像时的准确 性和精度。
发明内容
本发明的目的在于,提供一种各向异性介质的成像方法,它可以解决当前技术中存在的 问题,根据地下介质的各向异性这一实际情况,获得更准确的反演成像结果。
为解决上述技术问题,本发明采用如下的技术方案:一种各向异性介质的成像方法,包 括以下步骤:
S1,获得观测数据,根据观测数据获得速度场和各相异性参数;
S2,对观测数据进行偏移得到反射系数;
S3,根据速度场和各相异性参数,对反射系数进行的反偏移得到模拟数据;
S4,计算观测数据和模拟数据的数据残差;
S5,根据数据残差获得最终成像结果。
本发明提供的方法,充分考虑了地下介质的各向异性这一实际情况,能够更加精确的模 拟地震波在真实地下介质中的传播规律,基于线性反演思想通过反复迭代进一步提高偏移成 像结果的分辨率和保真度,使得偏移成像结果中不同深度的同相轴振幅更加均衡,成像结果 更为准确。此外,本发明继承了传统单程波类偏移成像方法高精度、高效率的优势,同时占 用较小的存储空间,偏移成像噪音少,对速度误差的敏感性更弱。
附图说明
图1为本发明的一个实施例流程示意图。
图2为本发明本发明的一个实施例中的速度场和各相异性参数图,其中
(a)为速度场,(b)为描述纵波水平速度和垂向速度关系的函数图像ε,(c)为描述纵 波各向异性强度δ图像;
图3a为对一个模型现有技术的成像结果,图3b为对前述模型本发明的成像结果。
下面结合附图和具体实施方式对本发明作进一步的说明。
具体实施方式
实施例1
一种各向异性介质的成像方法,如图1所示,包括以下步骤:
S1,获得观测数据,根据观测数据获得速度场和各相异性参数。
观测数据dobs根据各向异性Hess模型获得,模型主要由左侧的高速体岩丘,中间的强 各向异性异常体以及右侧的高陡断层组成。模型大小为901-376,纵横向网格间距均为10m。
观测数据共151炮,每炮401道接收,采用中间放炮两边接收的滚动观测系统,炮间隔 60m,道间隔10m。记录时间3s,时间采样间隔3ms。根据模型得到的速度和各向异性参数正演使用雷克子波,主频20Hz。
所述的各向异性参数包括描述纵波水平速度和垂向速度关系的函数ε,和,描述纵波各 向异性强度的参数δ,如图2所示,
S2,对观测数据进行偏移得到反射系数。
例如可以采用高阶傅里叶有限差分算子(FFD)偏移得到偏移成像结果,这里可以采用 的方法很多,属于常规技术,就不再赘述。
S3,根据速度场和各相异性参数,对反射系数进行反偏移得到模拟数据。
具体的说,根据得到的反射系数mk,根据速度场和各向异性参数进行正演模拟(也叫反 偏移)得到模拟数据Lmk。其中k表示迭代次数,L表示Born线性正演算子。
这里首先需要说明的是,地下波在反射时通常存在多次反射波,一次反射波是指地震波 向下传播时遇到地层界面形成的反射波,一次反射波是相对多次波而言的,但是除了一次波 以外的其他反射波通常能量较弱被当作噪音进行压制。因此,实际处理中所述的反射波通常 指一次反射波。
频率域中的反射波Us(xr|xs,ω)可以表示成以下公式
其中f(ω)是频率域震源子波,G(x'|xs,ω)是VTI介质中由震源位置xs传播到散射点x'的 Green函数,G(xr|x',ω)是由散射点x'传播到接收点xr的Green函数。m(x')是x'点的类反 射系数,即最小二乘偏移成像的目标。
需要说明的是,在各向异性介质中m(x')不仅仅与速度有关还受到各向异性参数ε,δ的影 响。U0(x'|xs,ω)是由震源传播到散射点x'的入射波场,Us(xr|xs,ω)是在检波点xr接受到的一 次反射波。结合公式(1),具体的反偏移过程如下:
S31,根据速度场和各相异性参数,从地表开始,对震源进行逐层延拓,在地下每一点 记录并存储入射波场,直到延拓到最深层。
具体的说,根据速度场和各相异性参数,以及震源(震源可由测数据得到),从地表开 始,利用高阶傅里叶有限差分(即FFD)算子震源进行逐层延拓得到入射波场,在地下每一 点记录并存储入射波场,直到延拓到最深层。
在数值模拟的过程中,会对地下介质进行网格剖分,所述的每一点在模拟过程中是指网 格的交点,其物理含义为地下介质的每一点。而所述的延拓即指的是模拟地震波在地下的传 播情况。所述的高阶FFD波场延拓包括频率波数域的相移,频率空间域时移和频率空间域有 限差分补偿,具体方式如下:
I,频率波数域的相移:
U1(kx,zj+1;ω)=U(kx,zj;ω)exp(ikz0Δz) (2)
其中U(kx,zj;ω)是深度zj处入射的总波场,j表示地层深度序号,kx表示水平波数,ω表示 圆频率,相移算子exp(ikz0Δz)中kz0表示背景垂直波数,i表示虚数单位,Δz表示延拓深度 间隔。经过相移处理延拓得到的出射波场U1(kx,zj+1;ω)为背景波场。
式(3)中c(z),ε0(z),δ0(z)是背景速度和各向异性参数。他们只随深度z发生变化。
II,频率空间域时移:
U2(x,zj+1;ω)=U1(x,zj+1;ω)exp(iωΔsΔz), (4)
其中U1(x,zj+1;ω)是由U1(kx,zj+1;ω)经过有关x的傅里叶反变换得到的。是介质的慢度 扰动。时移校正主要是对小角度入射的波场由于背景速度和真实速度之间的偏差引起的波场 传播误差进行补偿。
III,频率空间域有限差分补偿:
上式在计算时按照级联算法分两步进行具体计算
上计算得到的波场U(x,zj+1;ω)即为zj+1处最终的出射波场。这里的有限差分补偿主要是对由于 各向异性参数偏差引起的高角度传播的波场误差进行校正。其中i是虚数单位,l1,l2,l3是 各向异性参数的函数。具体转换公式如下:
ξ=(1+2ε),η=2(ε-δ),ξ0=(1+2ε0),η0=2(ε00) (7)
换言之,由公式(3)、(4)、(6)联合得到出射波场。需要说明的是,波场从上一层出射, 然后入射到下一层,上一层的出射波场就是下一层的入射波场,
S32,根据每一点的入射波场与反射系数,获得二次震源。
即会获得每一点上各自的二次震源,二次震源的作用是形成反射波,因此加载反射波的 过程已经包含了二次震源的加载。
S33,在最深层将反射波初始化为0,向地表逐层延拓,在延拓过程中,反射波场叠加来 自每一层的反射波,直至地表;
具体的说,在最深层将反射波初始化为0,并利用高阶FFD算子逐层延拓,在延拓过程 中,反射波场叠加上来自每一层的反射波,直至地表;这里在采用高阶FFD算子逐层延拓时 也包括如步骤S31中所述的频率波数域的相移,频率空间域时移和频率空间域有限差分补偿, 这里就不再重复说明。
S34,根据延拓到地表时的反射波场通过傅里叶反变换得到模拟数据。
S4,计算观测数据和模拟数据的数据残差,计算方法为模拟数据Lmk减去观测数据dobs即可。
S5,根据数据残差获得最终成像结果,具体的方法如下:
S51,判断数据残差是否小于给定阈值,如果是,则输出反射系数作为成像结果;否者 进入步骤S52。
S52,根据数据残差得到梯度,具体的说,包括如下步骤:
根据速度场v和各向异性参数,加载震源和当前数据残差,gk+1=LT(Lmk-dobs),其中LT表示偏移算子;
利用傅里叶有限差分正向延拓算子逐层延拓入射波场,利用反向延拓算子逐层延拓数据 残差;
对延拓到每一层的入射波场取共轭,并与延拓到相应位置的观测记录按照以下公式所示 的互相关成像条件提取成像值:
其中*表示取共轭;
最后根据成像值获得梯度值g。
利用傅里叶有限差分进行计算综合了效率和精度两方面的优势,傅里叶有限差分算法比 现有技术相比精度和效率更高。
在这个过程中的利用FFD算子逐层延拓也包含如步骤S31和S33中所述的频率波数域 的相移,频率空间域时移和频率空间域有限差分补偿,这里就不再重复说明。
S53,根据梯度计算共轭梯度方向的修正参数。
具体的说,根据梯度计算共轭梯度方向的修正参数其中P-1表示照明 预条件算子,T表示转置,共轭梯度是由梯度计算得到的,共轭梯度具有更好的全局收敛性, 梯度仅具有较好的局部收敛性。
S54,根据修正参数计算共轭梯度方向和共轭梯度更新步长。
即,计算共轭梯度方向zk+1=P-1gk+1+βzk,和共轭梯度更新步长
S55,利用共轭梯度方向和步长更新反射系数,并返回步骤S3。
即,利用共轭梯度方向和步长更新反射系数:mk+1=mk-αzk+1,并返回步骤S3,重新开 始迭代,直至满足S51中的终止条件数据残差小于给定阈值。最终得到的成像结果如图3b 所示。
传统的不考虑各向异性的方法所得到的结果图3a所示,对模型中间部分强各向异性区 域,偏移结果存在明显的绕射波不收敛现象,强各向异性体的上边界几乎没有对应同相轴。 而对比本发明的结果图3b来看,无论是对高陡的岩丘左侧边界,高陡断层,还是岩丘下深 部层状构造均达到了比较理想的成像效果。同时子波被压缩,剖面整体的同相轴更加细致清 晰,分辨率有明显提高。
通过采用本发明提供的方法,充分考虑了地下介质的各向异性这一实际情况,能够更加 精确的模拟地震波在真实地下介质中的传播规律,基于线性反演思想通过反复迭代,本发明 进一步提高偏移成像结果的分辨率和保真度,使得偏移成像结果中不同深度的同相轴振幅更 加均衡,成像结果更为准确。此外,本发明继承了传统单程波类偏移成像方法的优点,占用 较小的存储空间,偏移成像噪音少,对速度误差的敏感性更弱。
当然,上述实施例说明并非是对本发明的限制,本发明不仅限于上述举例,本技术领域 的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换,也应属于本发明的保 护范围。

Claims (6)

1.一种各向异性介质的成像方法,其特征在于,包括如下步骤:
S1,获得观测数据,根据观测数据获得速度场和各相异性参数;
S2,对观测数据进行偏移得到反射系数;
S3,根据速度场和各相异性参数,对反射系数进行反偏移得到模拟数据;
S4,计算观测数据和模拟数据的数据残差;
S5,根据数据残差获得最终成像结果。
2.如权利要求1所述的一种各向异性介质的成像方法,其特征在于,所述步骤S1中的各相异性参数包括:纵波水平速度和垂向速度关系的函数,和,纵波各向异性强度,所述纵波是根据所述观测数据得到的。
3.如权利要求2所述的一种各向异性介质的成像方法,其特征在于,所述步骤S3,根据速度场和各相异性参数,对反射系数进行反偏移得到模拟数据,具体包括如下步骤:
S31,根据速度场和各相异性参数,从地表开始,对震源进行逐层延拓,在地下每一点记录并存储入射波场,直到延拓到最深层;
S32,根据每一点的入射波场与反射系数,获得二次震源;
S33,在最深层将反射波初始化为0,向地表逐层延拓,在延拓过程中,反射波场叠加来自每一层的反射波,直至地表;
S34,根据延拓到地表时的反射波场得到模拟数据。
4.根据权利要求3所述的一种各向异性介质的成像方法,其特征在于,所述步骤S31中,对震源进行逐层延拓,具体包括:对延拓波场进行频率波数域的相移,频率空间域时移和频率空间域有限差分补偿。
5.如权利要求1所述的一种各向异性介质的成像方法,其特征在于,所述步骤S5中,根据数据残差获得最终成像结果,具体包括:
S51,判断数据残差是否小于给定阈值,如果是,则输出反射系数作为成像结果;否者进入步骤S52;
S52,根据数据残差得到梯度;
S53,根据梯度计算共轭梯度方向的修正参数
S54,根据修正参数计算共轭梯度方向和共轭梯度更新步长;
S55,利用共轭梯度方向和步长更新反射系数,并返回步骤S3。
6.如权利要求5所述的一种各向异性介质的成像方法,其特征在于,所述步骤S52中,根据数据残差得到梯度,具体包括:
S61,根据速度场v和各向异性参数,加载震源和当前数据残差;
S62,利用傅里叶有限差分正向延拓算子逐层延拓入射波场,利用反向延拓算子逐层延拓数据残差
S63,对延拓到每一层的入射波场取共轭,并与延拓到相应位置的观测记录按照以下公式所示的互相关成像条件提取成像值:
S64,根据成像值获得梯度值。
CN201710505548.5A 2017-06-28 2017-06-28 一种各向异性介质的成像方法 Expired - Fee Related CN107356972B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710505548.5A CN107356972B (zh) 2017-06-28 2017-06-28 一种各向异性介质的成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710505548.5A CN107356972B (zh) 2017-06-28 2017-06-28 一种各向异性介质的成像方法

Publications (2)

Publication Number Publication Date
CN107356972A true CN107356972A (zh) 2017-11-17
CN107356972B CN107356972B (zh) 2019-09-06

Family

ID=60273470

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710505548.5A Expired - Fee Related CN107356972B (zh) 2017-06-28 2017-06-28 一种各向异性介质的成像方法

Country Status (1)

Country Link
CN (1) CN107356972B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108562937A (zh) * 2018-03-15 2018-09-21 东北石油大学 一种地震成像方法
CN108845355A (zh) * 2018-09-26 2018-11-20 中国矿业大学(北京) 地震偏移成像方法及装置
CN110907993A (zh) * 2018-09-18 2020-03-24 中国石油化工股份有限公司 一种最小二乘偏移成像方法及系统
CN112630824A (zh) * 2019-10-09 2021-04-09 中国石油化工股份有限公司 一种地震成像中的离散点扩散函数生成方法及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2093591A1 (en) * 2008-02-22 2009-08-26 PGS Geophysical AS Method for Three Dimensional Seismic Travel Time Tomography in Transversely Isotropic Media
US20120099396A1 (en) * 2010-10-22 2012-04-26 Chevron U.S.A. Inc. System and method for characterization with non-unique solutions of anisotropic velocities
CN102508291A (zh) * 2011-10-17 2012-06-20 中国石油天然气股份有限公司 横向变速小尺度体反射系数公式应用方法
CN105319590A (zh) * 2014-07-30 2016-02-10 中国石油化工股份有限公司 一种基于hti介质的各向异性单参数反演方法
CN105487112A (zh) * 2014-09-18 2016-04-13 中国石油化工股份有限公司 一种构建地层反射系数的方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2093591A1 (en) * 2008-02-22 2009-08-26 PGS Geophysical AS Method for Three Dimensional Seismic Travel Time Tomography in Transversely Isotropic Media
US20120099396A1 (en) * 2010-10-22 2012-04-26 Chevron U.S.A. Inc. System and method for characterization with non-unique solutions of anisotropic velocities
CN102508291A (zh) * 2011-10-17 2012-06-20 中国石油天然气股份有限公司 横向变速小尺度体反射系数公式应用方法
CN105319590A (zh) * 2014-07-30 2016-02-10 中国石油化工股份有限公司 一种基于hti介质的各向异性单参数反演方法
CN105487112A (zh) * 2014-09-18 2016-04-13 中国石油化工股份有限公司 一种构建地层反射系数的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
李振春 等: ""基于平面波加速的VTI介质最小二乘逆时偏移"", 《地球物理学报》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108562937A (zh) * 2018-03-15 2018-09-21 东北石油大学 一种地震成像方法
CN110907993A (zh) * 2018-09-18 2020-03-24 中国石油化工股份有限公司 一种最小二乘偏移成像方法及系统
CN108845355A (zh) * 2018-09-26 2018-11-20 中国矿业大学(北京) 地震偏移成像方法及装置
CN112630824A (zh) * 2019-10-09 2021-04-09 中国石油化工股份有限公司 一种地震成像中的离散点扩散函数生成方法及系统
CN112630824B (zh) * 2019-10-09 2024-03-22 中国石油化工股份有限公司 一种地震成像中的离散点扩散函数生成方法及系统

Also Published As

Publication number Publication date
CN107356972B (zh) 2019-09-06

Similar Documents

Publication Publication Date Title
AU2012220584B2 (en) Sensitivity kernel-based migration velocity analysis in 3D anisotropic media
US8560242B2 (en) Pseudo-analytical method for the solution of wave equations
KR101797451B1 (ko) 상호상관 목적 함수를 통한 해양 스트리머 데이터에 대한 동시 소스 반전
US10234582B2 (en) Joint inversion of seismic data
CN107356972B (zh) 一种各向异性介质的成像方法
CN107462924B (zh) 一种不依赖于测井资料的绝对波阻抗反演方法
CN108241173B (zh) 一种地震资料偏移成像方法及系统
EP3067718B1 (en) Boundary layer tomography method and device
CN102116869A (zh) 高精度叠前域最小二乘偏移地震成像技术
CN109633749A (zh) 基于散射积分法的非线性菲涅尔体地震走时层析成像方法
CN109655890B (zh) 一种深度域浅中深层联合层析反演速度建模方法及系统
Guo et al. Topography-dependent eikonal tomography based on the fast-sweeping scheme and the adjoint-state technique
Li et al. Waveform inversion of seismic first arrivals acquired on irregular surface
CN104280774B (zh) 一种单频地震散射噪声的定量分析方法
CN108919351A (zh) 基于逆时聚焦原理进行观测系统双向聚焦性的评价方法
Wang et al. Monochromatic wave-equation traveltime inversion in vertical transverse isotropic media: Theory and application
Kazei et al. Acquisition and near-surface impacts on VSP mini-batch FWI and RTM imaging in desert environment
Bychkov et al. Near-surface correction on seismic and gravity data
Xing et al. Application of 3D stereotomography to the deep-sea data acquired in the South China Sea: a tomography inversion case
Ni* et al. Preliminary practice of stereotomography
Lee et al. Ongoing development of a 3D seismic velocity model for Canterbury, New Zealand for broadband ground motion simulation
Barnes et al. Diving wave tomography: A robust method for velocity estimation in a foothills geological context
Shi et al. True-amplitude Gaussian-beam Migration for Acoustic Transversely Isotropic Media with a Vertical Symmetry Axis
Taillandier et al. Refraction traveltime tomography based on adjoint state techniques
Yang et al. Least-squares extended reverse time migration in the presence of velocity errors

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20190906

Termination date: 20200628

CF01 Termination of patent right due to non-payment of annual fee