CN109946741B - 一种TTI介质中纯qP波最小二乘逆时偏移成像方法 - Google Patents

一种TTI介质中纯qP波最小二乘逆时偏移成像方法 Download PDF

Info

Publication number
CN109946741B
CN109946741B CN201910250669.9A CN201910250669A CN109946741B CN 109946741 B CN109946741 B CN 109946741B CN 201910250669 A CN201910250669 A CN 201910250669A CN 109946741 B CN109946741 B CN 109946741B
Authority
CN
China
Prior art keywords
field
imaging
reverse time
reflection coefficient
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.)
Active
Application number
CN201910250669.9A
Other languages
English (en)
Other versions
CN109946741A (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 CN201910250669.9A priority Critical patent/CN109946741B/zh
Publication of CN109946741A publication Critical patent/CN109946741A/zh
Application granted granted Critical
Publication of CN109946741B publication Critical patent/CN109946741B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Geophysics And Detection Of Objects (AREA)

Abstract

本说明书实施例公开了一种TTI介质中纯qP波最小二乘逆时偏移成像方法。本发明能够考虑具有倾斜对称轴的复杂各向异性介质实际波场传播情况,通过提供一种稳定的TTI介质最小二乘偏移算子和逆时反偏移算子,既能克服传统成像方法在处理存在复杂TTI介质时的不足,又通过最小二乘方式迭代方式压制了偏移成像噪声,提高了成像分辨率,得到保幅的成像剖面,开发基于TTI介质的最小二乘逆时偏移成像技术,为存在复杂各向异性探区的地震数据处理提供了高精度成像保障,提高了后续解释工作的质量。

Description

一种TTI介质中纯qP波最小二乘逆时偏移成像方法
技术领域
本说明书涉及勘探地球物理学领域,尤其涉及一种TTI介质中纯qP波最小二乘逆时偏移成像方法。
背景技术
地下介质广泛发育各向异性特征,有倾斜对称轴的横向各向同性(tiltedtransverse isotropy,TTI)介质分布广泛。忽略各向异性引起的地震波传播旅行时和相位差异将严重影响成像结果的准确性。因此为得到高精度的成像结果,需校正各向异性对波场传播的影响。常规逆时偏移成像在处理各向异性介质时,存在振幅不均衡、分辨率低、偏移噪声严重等问题。
基于此,需要一种适合于TTI介质中的最小二乘逆时偏移成像方法。
发明内容
本发明的目的在于,提供一种准确、稳定的TTI介质中纯qP波最小二乘逆时偏移成像方法。
为解决上述技术问题,本发明采用如下技术方案:
输入偏移参数场vp0,各向异性Thomsen参数模型ε和δ,各向异性倾角θ和观测数据dobs
采用预设的反传波场传播算子LT,对所述观测数据进行逆时偏移,得到初始反射系数场m(1),其中,所述反传波场传播算子LT与偏移参数相关;
采用如下方式对反射系数场m进行迭代,直至数据残差小于阈值:
根据预设的线性正演传播算子L和第k次迭代得到的反射系数场m(k)计算逆时反偏移记录Lm(k)
确定数据残差为Lm(k)-dobs,当数据残差大于阈值时,确定更新最速下降法下降方向g(k+1)和步长β(k+1),并根据最速下降法下降方向和步长计算共轭梯度法下降方向α(k+1)和步长y(k+1),更新反射系数场m(k+1)=m(k)(k+1)y(k+1)
当迭代停止时时,输出当前的反射系数场m为成像结果;
其中,反传波场传播算子LT的表达式为:
Figure BDA0002012314370000021
其中,p*(x,t)和q*(x,t)分别是背景波场p0(x,t)和背景辅助波场q0(x,t)的伴随波场,△precv(x,t)表示观测数据和正演数据的残差,x和z分别表示横坐标和纵向坐标,t表示波场传播时间;
线性正演传播算子L的表达式为:
Figure BDA0002012314370000022
其中,
Figure BDA0002012314370000031
为反射系数,vp0为背景速度,vps为速度扰动,ps(x,t)为应力场扰动波场,qs(x,t)为辅助波场的扰动波场。
本说明书实施例采用的上述至少一个技术方案能够达到以下有益效果:
与现有技术相比,本发明能够考虑具有倾斜对称轴的复杂各向异性介质实际波场传播情况,通过提供一种稳定的TTI介质最小二乘偏移算子和逆时反偏移算子,既能克服传统成像方法在处理存在复杂TTI介质时的不足,又通过最小二乘方式迭代方式压制了偏移成像噪声,提高了成像分辨率,得到保幅的成像剖面,开发基于TTI介质的最小二乘逆时偏移成像技术,为存在复杂各向异性探区的地震数据处理提供了高精度成像保障,提高了后续解释工作的质量。
附图说明
图1为本说明书实施例所提供的流程示意图;
图2为真实速度vp模型;
图3为各向异性Thomsen参数ε模型;
图4为各向异性Thomsen参数δ模型;
图5为各向异性倾角θ模型;
图6为偏移速度模型;
图7为相对慢度扰动模型;
图8为TTI介质逆时偏移成像结果;
图9为各向同性最小二乘逆时偏移成像结果;
图10为VTI介质中的最小二乘逆时偏移成像结果;
图11为本说明书实施例所提供的成像结果;
图12为本说明书实施例与其他方法的迭代收敛速度对比示意图。
具体实施方式
为使本申请的目的、技术方案和优点更加清楚,下面将结合本申请具体实施例及相应的附图对本申请技术方案进行清楚、完整地描述。显然,所描述的实施例仅是本申请一部分实施例,而不是全部的实施例。基于本说明书中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本申请保护的范围。
首先对本说明书实施例中所使用的反传波场传播算子LT和线性正演传播算子L进行解释说明。
在TTI介质的qP波正演模拟中一般采用如下的波动方程:
Figure BDA0002012314370000041
其中,vp表示qP波传播速度,p(x,t)表示应力场值,q(x,t)表示辅助应力场值,θ表示各向异性倾角,ε和δ表示介质Thomsen参数值,x和z分别表示纵横向坐标,f(xs,t)表示震源项,xs表示震源位置,t表示波场传播时间。
将辅助方程
Figure BDA0002012314370000042
在波数域进行求解,求解方式如下:
Figure BDA0002012314370000043
其中,FFT表示快速傅里叶变换,FFT-1表示快速傅里叶变换逆变换,kx和kz分别表示纵横向波数,pn表示离散的应力场,qn表示离散的辅助应力场。
时间差分离散形式:
Figure BDA0002012314370000051
其中,△t表示时间采样间隔,pn表示此时刻应力场值,pn-1表示上一时刻应力场值,pn+1表示下一时刻应力场值。
公式(1)的空间差分离散形式为:
Figure BDA0002012314370000052
其中,
Figure BDA0002012314370000053
分别表示不同偏导方式的空间四阶混合偏导数。i和j为x和z方向的空间坐标。
带入波动方程(1),得到TTI介质正演模拟差分递推公式:
Figure BDA0002012314370000054
求取正演算子的伴随算子LT(即反传波场传播算子),也即反传波场满足的波动方程如下:
Figure BDA0002012314370000061
其中,p*(x,t)和q*(x,t)分别是背景波场p0(x,t)和背景辅助波场q0(x,t)的伴随波场,△precv(x,t)表示观测数据和正演数据的残差。
Born假设下的线性正演炮记录如下,即线性正演传播算子L的表达式如下:
Figure BDA0002012314370000062
其中,
Figure BDA0002012314370000063
为定义的反射系数,vp0为背景速度,vps为速度扰动,ps(x,t)为应力场扰动波场,qs(x,t)为辅助波场的扰动波场,t和x分别表示空间坐标和时间
可以看到,本说明书实施例所提供的算子基于各向异性参数,换言之,在本说明书实施例所提供的方案考虑了各向异性影响。
上述部分对于本说明书实施例中所采用的算子的原理和形式进行了具体的说明。基于前述内容,本说明书的实施例提供的一种TTI介质中纯qP波最小二乘逆时偏移成像方法,如图1所示,该过程具体包括以下步骤:
S101,输入偏移参数场vp0,各向异性Thomsen参数模型ε和δ,各向异性倾角θ和观测数据dobs
需要说明的是,在本说明书实施例中,vp0,各向异性Thomsen参数模型ε和δ,各向异性倾角θ均和位置有关,即,与横坐标x和纵坐标z有关。
S103,采用预设的反传波场传播算子LT,对所述观测数据进行逆时偏移,得到初始反射系数场m(1),其中,所述反传波场传播算子LT与偏移参数相关。
反传波场传播算子LT的原理和表达形式已经在前文进行了说明,此处不再赘述。
S105,采用如下方式对反射系数场m进行迭代,直至数据残差小于阈值:根据预先得到的线性正演传播算子L和第k次迭代得到的反射系数场m(k)计算逆时反偏移记录Lm(k)
确定数据残差Lm(k)-dobs,当数据残差大于阈值时,确定更新最速下降法下降方向g(k+1)(最速下降法下降方向即梯度)和步长β(k+1),并根据最速下降法下降方向和步长计算共轭梯度法下降方向α(k+1)和步长y(k+1),更新反射系数场m(k+1)=m(k)(k+1)y(k+1)
具体而言,本发明建立在观测记录与反偏移数据的L2模的拟合基础上,其目标泛函如下:
Figure BDA0002012314370000071
在求取下降方向时,可以采用多种方式实现。在不同的方式下,下降方向及其步长的形式不同。
以使用共轭梯度法实现TTI介质最小二乘逆时偏移为例,基本迭代公式为:
m(k+1)=m(k)(k+1)y(k+1) (9)
其中,m表示所求的反射系数;α(k+1)和y(k+1)分别表示第k+1次迭代的共轭梯度法下降方向和步长,共轭梯度法是在最速下降法基础上发展而来,其具有比最速下降法收敛速度更快的优点。因此,共轭梯度法下降方向和步长可通过如下四个方程求解:
g(k+1)=LT[Lm(k)-dobs] (10)
Figure BDA0002012314370000081
y(k+1)=g(k+1)(k)y(k) (12)
Figure BDA0002012314370000082
其中,g(k+1)和β(k+1)表示第k+1次迭代的最速下降法下降方向(即梯度)和步长。
其中,梯度的计算方式为
Figure BDA0002012314370000083
即,梯度与L和LT相关。
S107,当迭代停止时时,输出当前的反射系数场为成像结果。
与现有技术相比,本发明能够考虑具有倾斜对称轴的复杂各向异性介质实际波场传播情况,通过提供一种稳定的TTI介质最小二乘偏移算子和逆时反偏移算子,既能克服传统成像方法在处理存在复杂TTI介质时的不足,又通过最小二乘方式迭代方式压制了偏移成像噪声,提高了成像分辨率,得到保幅的成像剖面,开发基于TTI介质的最小二乘逆时偏移成像技术,为存在复杂各向异性探区的地震数据处理提供了高精度成像保障,提高了后续解释工作的质量。
以下给出本说明是实施例在模型中的实际效果说明。
将本发明所提供的方法应用于国际标准Marmousi模型成像,取得了较理想的成像效果。真实速度模型(如图2所示)、各向异性Thomsen参数Epsilon模型(如图3所示)、各向异性Thomsen参数Delta模型(如图4示)、各向异性倾角模型(如图5所示)、偏移速度模型(如图6所示)、慢度扰动模型(如图7所示);建立全接收的观测系统,输入慢度扰动,偏移速度场和各向异性参数,线性正演模拟得到观测炮记录;输入偏移速度场和各向异性参数进行波场正传,并反传观测炮记录,采用互相关成像条件将正传和反传波场互相关,得到常规TTI介质逆时偏移成像结果(如图8所示)。
采用本说明书实施例所提供的方案,对常规逆时偏移成像结果应用TTI介质逆时反偏移算子进行反偏移,得到Born近似下的线性正演炮记录,将模拟炮记录与实际观测炮记录相减,判断数据残差是否小于给定阈值,若大于给定阈值则求取梯度更新方向和步长,然后更新成像剖面,再次进行逆时反偏移并计算模拟记录与观测记录数据残差,直至数据残差小于给定阈值,如果残差小于给定阈值则输出最终的偏移成像剖面,如图11所示。
相比各向同性声波最小二乘偏移结果(如图9所示)和VTI介质最小二乘偏移结果(如图10所示),本说明书实施例的方案考虑了所有各向异性因素的影响,模拟数据与观测数据匹配较好,成像结果随着迭代次数增加变得更加理想。反射同相轴成像位置准确,波形畸变引起的同相轴成像紊乱被消除。相比TTI介质逆时偏移成像结果(如图8所示),基本压制了浅层偏移噪音和震源采集脚印,分辨率得到明显提高,成像剖面中深层振幅更加均衡。在前述示意图中,distance对应于横坐标x,depth对应于纵坐标z。
图12为用不同成像方法得到的数据残差收敛曲线,在该示意图中,横坐标为迭代次数,纵坐标为第n次迭代的数据残差与初始数据残差的百分比。可以看出来,本发明的成像方法在迭代过程中的数据残差收敛曲线相比于其他成像方法收敛速度更快,最终稳定的残差值更小
为本说明书实施例所提供的方案,为复杂地质构造解释工作提供较准确的成像基础,为老油田二次勘探开发提供强有力的技术支持。
对应的,本申请实施例还提供一种计算机设备,所述设备包括包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,其中,所述处理器执行所述程序时实现如前所述的TTI介质中纯qP波最小二乘逆时偏移成像方法所述的方法。
本说明书中的各个实施例均采用递进的方式描述,各个实施例之间相同相似的部分互相参见即可,每个实施例重点说明的都是与其他实施例的不同之处。尤其,对于装置、设备和介质类实施例而言,由于其基本相似于方法实施例,所以描述的比较简单,相关之处参见方法实施例的部分说明即可,这里就不再一一赘述。
上述对本说明书特定实施例进行了描述。其它实施例在所附权利要求书的范围内。在一些情况下,在权利要求书中记载的动作或步骤或模块可以按照不同于实施例中的顺序来执行并且仍然可以实现期望的结果。另外,在附图中描绘的过程不一定要求示出的特定顺序或者连续顺序才能实现期望的结果。在某些实施方式中,多任务处理和并行处理也是可以的或者可能是有利的。

Claims (4)

1.一种TTI介质中纯qP波最小二乘逆时偏移成像方法,包括:
输入偏移参数场vp0,各向异性Thomsen参数模型ε和δ,各向异性倾角θ和观测数据dobs
采用预设的反传波场传播算子LT,对所述观测数据进行逆时偏移,得到初始反射系数场m(1),其中,所述反传波场传播算子LT与偏移参数相关;
采用如下方式对反射系数场m进行迭代,直至数据残差小于阈值:
根据预设的线性正演传播算子L和第k次迭代得到的反射系数场m(k)计算逆时反偏移记录Lm(k)
确定数据残差为Lm(k)-dobs,当数据残差大于阈值时,确定更新最速下降法下降方向g(k +1)和步长β(k+1),并根据最速下降法下降方向和步长计算共轭梯度法下降方向α(k+1)和步长y(k+1),更新反射系数场m(k+1)=m(k)(k+1)y(k+1)
当迭代停止时,输出当前的反射系数场m为成像结果;
其中,反传波场传播算子LT的表达式为:
Figure FDA0002528698850000011
其中,p*(x,t)和q*(x,t)分别是背景波场p0(x,t)和背景辅助波场q0(x,t)的伴随波场,Δprecv(x,t)表示观测数据和正演数据的残差,x和z分别表示横坐标和纵向坐标,t表示波场传播时间;
线性正演传播算子L的表达式为:
Figure FDA0002528698850000021
其中,
Figure FDA0002528698850000022
为反射系数,vps为速度扰动,ps(x,t)为应力场扰动波场,qs(x,t)为辅助波场的扰动波场。
2.如权利要求1所述的方法,所述最速下降法下降方向g由如下方式确定:
Figure FDA0002528698850000023
其矩阵形式为g(k+1)=LT[Lm(k)-dobs];所述最速下降法的下降步长为,
Figure FDA0002528698850000024
β(k)为第k次迭代时的步长。
3.如权利要求2所述的方法,所述共轭梯度法中的下降方向为
Figure FDA0002528698850000025
共轭梯度法中的下降步长为y(k+1)=g(k+1)(k)y(k)
4.一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,其中,所述处理器执行所述程序时实现如权利要求1至3任一所述的方法。
CN201910250669.9A 2019-03-29 2019-03-29 一种TTI介质中纯qP波最小二乘逆时偏移成像方法 Active CN109946741B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910250669.9A CN109946741B (zh) 2019-03-29 2019-03-29 一种TTI介质中纯qP波最小二乘逆时偏移成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910250669.9A CN109946741B (zh) 2019-03-29 2019-03-29 一种TTI介质中纯qP波最小二乘逆时偏移成像方法

Publications (2)

Publication Number Publication Date
CN109946741A CN109946741A (zh) 2019-06-28
CN109946741B true CN109946741B (zh) 2020-09-11

Family

ID=67012958

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910250669.9A Active CN109946741B (zh) 2019-03-29 2019-03-29 一种TTI介质中纯qP波最小二乘逆时偏移成像方法

Country Status (1)

Country Link
CN (1) CN109946741B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112285778B (zh) * 2020-10-29 2022-05-27 中国石油大学(华东) 一种粘声TTI介质中纯qP波的逆时偏移成像方法
US11733413B2 (en) 2021-04-30 2023-08-22 Saudi Arabian Oil Company Method and system for super resolution least-squares reverse time migration
CN114924313B (zh) * 2022-05-16 2022-12-13 中国海洋大学 基于行波分离的声波最小二乘逆时偏移梯度求取方法
CN115951401B (zh) * 2022-07-19 2023-09-15 中山大学 成像条件驱动的最小二乘逆时偏移成像方法、设备及存储介质
CN115201913B (zh) * 2022-07-27 2023-05-12 中山大学 基于无网格有限差分法的最小二乘逆时偏移成像方法、系统及存储介质

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20160306058A1 (en) * 2015-04-16 2016-10-20 The Board Of Trustees Of The Leland Stanford Junior Univerisity Reverse time migration based on geometric mean for imaging seismic sources
CN105717539B (zh) * 2016-01-28 2018-01-30 中国地质大学(北京) 一种基于多gpu计算的三维tti介质逆时偏移成像方法
CN108333628B (zh) * 2018-01-17 2019-09-03 中国石油大学(华东) 基于正则化约束的弹性波最小二乘逆时偏移方法

Also Published As

Publication number Publication date
CN109946741A (zh) 2019-06-28

Similar Documents

Publication Publication Date Title
CN109946741B (zh) 一种TTI介质中纯qP波最小二乘逆时偏移成像方法
CN108873066B (zh) 弹性介质波动方程反射波旅行时反演方法
US8892410B2 (en) Estimation of soil properties using waveforms of seismic surface waves
CN111221037B (zh) 解耦弹性逆时偏移成像方法和装置
US9625593B2 (en) Seismic data processing
CN105974470A (zh) 一种多分量地震资料最小二乘逆时偏移成像方法及系统
US20120275267A1 (en) Seismic Data Processing
CN105652321A (zh) 一种粘声各向异性最小二乘逆时偏移成像方法
CN110531410B (zh) 一种基于直达波场的最小二乘逆时偏移梯度预条件方法
CN112327358B (zh) 一种粘滞性介质中声波地震数据正演模拟方法
CN113064203B (zh) 共轭梯度归一化lsrtm方法、系统、存储介质及应用
CN112965102A (zh) 一种基于阻抗敏感核函数的最小二乘逆时偏移成像方法
CN111665556B (zh) 地层声波传播速度模型构建方法
CN106842300A (zh) 一种高效率多分量地震资料真振幅偏移成像方法
Lu et al. Broadband least-squares wave-equation migration
CN112285778B (zh) 一种粘声TTI介质中纯qP波的逆时偏移成像方法
CN113866823A (zh) 一种粘声各向异性介质中的正演成像方法
CN111175822B (zh) 改进直接包络反演与扰动分解的强散射介质反演方法
CN110888158B (zh) 一种基于rtm约束的全波形反演方法
CN114609671A (zh) 鬼波衰减方法、装置、计算机设备及可读存储介质
CN114325829B (zh) 一种基于双差思想的全波形反演方法
CN111665550A (zh) 地下介质密度信息反演方法
CN111665546A (zh) 用于可燃冰探测的声学参数获取方法
CN111665549A (zh) 地层声波衰减因子反演方法
CN113687415B (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