CN105974470A - 一种多分量地震资料最小二乘逆时偏移成像方法及系统 - Google Patents
一种多分量地震资料最小二乘逆时偏移成像方法及系统 Download PDFInfo
- Publication number
- CN105974470A CN105974470A CN201610520574.0A CN201610520574A CN105974470A CN 105974470 A CN105974470 A CN 105974470A CN 201610520574 A CN201610520574 A CN 201610520574A CN 105974470 A CN105974470 A CN 105974470A
- Authority
- CN
- China
- Prior art keywords
- big gun
- wave field
- record
- components
- 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
Links
- 238000013508 migration Methods 0.000 title claims abstract description 201
- 230000005012 migration Effects 0.000 title claims abstract description 201
- 238000003384 imaging method Methods 0.000 title claims abstract description 37
- 238000000034 method Methods 0.000 claims abstract description 86
- 230000008878 coupling Effects 0.000 claims description 24
- 238000010168 coupling process Methods 0.000 claims description 24
- 238000005859 coupling reaction Methods 0.000 claims description 24
- 238000005457 optimization Methods 0.000 claims description 24
- 238000000354 decomposition reaction Methods 0.000 claims description 16
- 230000015572 biosynthetic process Effects 0.000 claims description 12
- 238000005070 sampling Methods 0.000 claims description 12
- 238000004422 calculation algorithm Methods 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000001514 detection method Methods 0.000 claims description 3
- 238000001615 p wave Methods 0.000 claims description 3
- 238000004064 recycling Methods 0.000 claims description 3
- 238000004088 simulation Methods 0.000 claims description 3
- 230000001066 destructive effect Effects 0.000 abstract description 4
- 238000012986 modification Methods 0.000 abstract 1
- 230000004048 modification Effects 0.000 abstract 1
- 238000004321 preservation Methods 0.000 abstract 1
- 230000000694 effects Effects 0.000 description 5
- 150000003839 salts Chemical class 0.000 description 4
- 230000008569 process Effects 0.000 description 3
- 239000011800 void material Substances 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 238000004880 explosion Methods 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000003989 dielectric material Substances 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 230000005477 standard model Effects 0.000 description 1
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
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
- G01V2210/512—Pre-stack
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
技术领域
本发明属于地球物理勘探领域,涉及多分量地震资料叠前偏移成像处理,特别涉及一种多分量地震资料最小二乘逆时偏移成像方法及系统。
背景技术
在过去的几十年中,油气勘探中绝大多数以纵波勘探方法为主,该方法是基于声学介质假设的,它认为地下介质中只存在纵波。然而,随着油气勘探的深入发展,勘探的目标日趋复杂,这种基于声波方程的单分量地震数据处理方法越来越显得有些能力不足,尤其对于复杂构造成像。因为,地下介质中传播的波场不仅有纵波,而且还有横波以及转换波等波型,也就是说地震波场是弹性波场。多分量地震资料自然地包含了纵波,横波以及转换波等波型,与单一的纵波相比,纵横波之间的耦合可以更好的保持地震波场的运动学(走时、路径等)和动力学(波形、振幅、频率、相位、偏振特性等)特性。此外,由于横波的存在,使得多分量地震资料可以更有效地识别亮点反射,更好的进行储层流体识别、裂缝分布估计、各向异性分析等。因此针对多分量地震资料,发展基于弹性介质假设的偏移成像方法成为今后偏移成像的研究重点。
众所周知,叠前深度偏移已经成为当前油气勘探业界的主流方法,根据偏移过程中所用波场延拓算子的不同,目前已有的多分量地震资料叠前深度偏移成像方法可以分为三类。一是射线类偏移方法,如Kirchhoff偏移和束偏移等;二是单程波类偏移方法,如弹性屏法等;三是基于双程波动方程的逆时偏移方法。对于复杂构造来讲,逆时偏移是精度最高,算法最稳健的一种。对于多分量地震资料逆时偏移,其实现方式有多种,可以概括为两类:一是基于声波方程,如将多分量地震资料分离成纵波资料和横波资料,对每一种资料都采用传统的声波逆时偏移进行成像。二是基于弹性波方程,如直接对多分量资料进行笛卡尔分量成像,或利用波场的矢量势和标量势进行成像等。对基于声波方程的多分量地震资料偏移方法而言,该实现方式忽略了多分量数据的矢量特性,造成地下介质属性信息被忽略,这必然带来成像误差。此外,对多分量地震资料进行的纵横波分离往往不是完美的,这种不彻底的分离也必然会造成偏移剖面中出现串扰假象。对基于弹性波方程的多分量地震资料偏移方法而言,纵横波之间的串扰,横波极性反转等问题,在不同程度和方面影响着偏移的精度。除此之外,在逆时偏移中所采用的波场延拓算子存在误差,它并非波场正向延拓算子的精确逆,而是利用波场正向延拓算子的伴随算子来替代的。再加上,实际地震资料是不完整的,采集孔径有限,而且资料中存在假频、噪音等,综合这些因素,使得逆时偏移剖面存在采集脚印、分辨率低、振幅不均衡、假象严重等问题。因此,基于常规的弹性波逆时偏移方法无法很好的实现对多分量地震资料的高精度、高分辨率、搞信噪比、保幅成像。为此,必须建立一套新的基于弹性波方程的、以“高精度、高分辨率、高信噪比、保幅”为特点的多分量地震资料逆时偏移成像方法。
发明内容
本发明的目的在于针对现有技术存在的上述缺陷,提供一种基于弹性波方程的、以“高精度、高分辨率、高信噪比、保幅”为特点的多分量地震资料最小二乘逆时偏移成像方法及系统。
为了实现以上目的,本发明提供的一种多分量地震资料最小二乘逆时偏移成像方法,包括以下步骤:
(a):读取预设参数、给定的背景模型及观测的多炮多分量地震记录D=(Dx,Dy,Dz),确定进行偏移成像的所有多分量地震记录;
(b):针对每一炮,在该炮对应炮点位置设置震源子波,基于各向同性介质纵横波分解的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的每一个时刻的震源波场,该震源波场包含矢量纵波震源波场;从而获得每一炮对应的震源波场;
(c):针对每一炮,首先对与该炮对应的观测的多分量地震记录进行预处理,之后同样基于各向同性介质纵横波分解的弹性波方程,对该炮的经过预处理的多分量地震记录进行波场逆时延拓,获得该炮的每一个时刻的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于该炮的震源波场以及检波点波场应用成像条件,获得该炮的单炮偏移剖面;进而获得每一炮对应的单炮偏移剖面;将所有炮的单炮偏移剖面进行叠加,获得初始偏移剖面,也即第0次迭代的偏移剖面
所述的进行预处理具体为:
其中,t表示波的传播时间;
(d):针对每一炮,在该炮对应炮点位置设置震源子波,基于各向同性介质纵横波耦合的弹性波方程,利用(a)中给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用(c)中所述的初始偏移剖面I0、所述的该炮的背景波场以及(a)中给定的背景模型,构建该炮的虚拟震源;基于带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用给定的背景模型对该炮点进行波正向延拓,获得该炮的反偏移波场,对该炮的反偏移波场进行记录采样,获得该炮的多分量反偏移记录;进而获得每一炮对应的多分量反偏移记录,所有炮的多分量反偏移记录构成了初始多分量反偏移记录,也即第0次迭代的多分量反偏移记录此时,设置当前的迭代次数为i=0;
(e):设置当前的迭代次数i=i+1,基于第i-1次迭代所得到的所述偏移剖面和第i-1次迭代所得到的所述多分量反偏移记录进行反演迭代;
在第i次迭代中,针对每一炮,利用第i-1次迭代所得到的该炮的多分量反偏移记录di-1和与该炮对应的观测的多分量地震记录计算该炮的多分量残差记录,进而获得所有炮对应的多分量残差记录
所述的利用第i-1次迭代所得到的该炮的多分量反偏移记录di-1和与该炮对应的观测的多分量地震记录计算该炮的多分量残差记录具体为:
(f):针对每一炮,首先对该炮的多分量残差记录进行预处理,之后同样基于各向同性介质纵横波分解的弹性波方程,对经过预处理的多分量残差记录进行波场逆时延拓获得该炮的每一个时刻的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于该炮的检波点波场以及该炮的由步骤(b)得到的震源波场应用所述成像条件,获得该炮的单炮梯度剖面;进而获得每一炮对应的单炮梯度剖面;将所有炮的单炮梯度剖面进行叠加,获得本次迭代的梯度剖面,也即第i次迭代的梯度剖面
所述的对该炮的所述多分量残差记录进行预处理具体为:
(g):利用最优化反演算法,构建第i次迭代的下降方向剖面
(h):针对每一炮,在该炮对应炮点位置设置震源子波,同样基于各向同性介质纵横波耦合的弹性波方程,利用(a)中给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用步骤(g)中所获得的第i次迭代的下降方向剖面ri、该炮的背景波场以及给定的背景模型,构建该炮的虚拟震源;基于带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用给定的背景模型对该炮点进行波场正向延拓,获得该炮的反偏移波场,对该炮的反偏移波场进行记录采样,获得该炮的多分量扰动反偏移记录;进而获得每一炮对应的多分量扰动反偏移记录,所有炮的多分量扰动反偏移记录构成了第i次迭代的多分量扰动反偏移记录再利用步长公式计算本次迭代的优化步长αi;
所述的步长公式为
所述方程(4)中,求和变量k包含笛卡尔坐标的x,y和z三个方向;xs表示震源点位置,xr表示检波点位置;
(i):利用步骤(h)得到的优化步长αi及步骤(g)得到的下降方向剖面ri,更新第i次迭代的偏移剖面Ii=Ii-1+αiri;利用优化步长及步骤(h)得到的多分量扰动反偏移记录δdi,更新第i次迭代的多分量反偏移记录di=di-1+αiδdi;
所述更新第i次迭代的所述偏移剖面Ii=Ii-1+αiri具体为:
所述更新第i次迭代的所述多分量反偏移记录di=di-1+αiδdi具体为:
(j):第i次迭代完后,计算第i次迭代的目标泛函值fi,判断当前迭代是否满足收敛标准,如果满足则输出最新的偏移剖面为最终的偏移剖面I=(Ipp,Ips);否则重复步骤(e)至(i),直至获得最终偏移剖面;
所述计算第i次迭代的目标泛函值fi具体为:
所述方程(7)中,求和变量k包含笛卡尔坐标的x,y和z三个方向;xs表示震源点位置,xr表示检波点位置;
所述的收敛标准具体为:
所述方程(8)中,theshold表示迭代停止的阈值标准,通常选取0.00001为最佳。
本发明实施采用的技术方案还包括:在步骤(b),步骤(c)和步骤(f)中所述的各向同性介质纵横波分解的弹性波方程为:
所述方程(9)中,Vp=(Vxp,Vyp,Vzp)T表示矢量纵波速度波场;Vs=(Vxs,Vys,Vzs)T表示矢量横波速度波场;V=Vp+Vs=(Vx,Vy,Vz)T表示矢量速度波场;σ=(σxx,σyy,σzz,σxy,σyz,σzx)T表示矢量应力波场;λ和μ表示介质拉梅常数;ρ表示介质密度;t表示时间;x,y和z分别表示笛卡尔坐标的x,y和z三个方向。
本发明实施采用的技术方案还包括:在步骤(d)和步骤(h)中所述的各向同性介质纵横波耦合的弹性波方程为:
本发明实施采用的技术方案还包括:步骤(c)和步骤(f)中所述成像条件为:
所述方程(11)中,在步骤(c)中,在步骤(f)中,
所述方程(11)中,表示震源波场中的矢量纵波速度波场,表示检波点波场中的矢量纵波速度波场,表示检波点波场中的矢量横波速度波场。
本发明实施采用的技术方案还包括:在步骤(d)和步骤(h)中所述基于带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用给定的背景模型对该炮点进行波正向延拓,获得该炮的矢量反偏移波场,其特征在于此处的带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程可表示为:
所述方程(12)中,虚拟震源矢量F=(Fxx,Fyy,Fzz,Fxy,Fyz,Fzx)T可以表示为:
所述方程(12)中,δV=(δVx,δVy,δVz)T表示反偏移矢量速度波场;δσ=(δσxx,δσyy,δσzz,δσxy,δσyz,δσzx)T表示反偏移矢量应力波场;
所述方程(13)中,在步骤(d)中,在步骤(h)中,m1=ri pp,m2=ri ps。
本发明还提供一种多分量地震资料最小二乘逆时偏移成像系统,包括初始化模块、震源波场计算模块、初始偏移剖面计算模块、初始多分量反偏移记录计算模块和反演迭代模块;
所述初始化模块用于读取预设参数、给定的背景模型及观测的多分量地震记录,确定进行偏移成像的 所有多分量地震记录;
所述震源波场计算模块用于利用各向同性介质纵横波分解的弹性波方程进行波场正向延拓获得每一炮对应的每一个时刻的震源波场;
所述初始偏移剖面计算模拟用于计算所有炮的叠加偏移剖面,具体为针对每一炮,首先对所述该炮的观测的多分量地震记录进行预处理,之后基于所述各向同性介质纵横波分解的弹性波方程,对经过所述预处理的该炮的多分量地震记录进行波场逆时延拓,获得该炮的每一个时刻的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于所述该炮的震源波场以及所述该炮的检波点波场应用成像条件,获得该炮的单炮偏移剖面;对所有炮执行该操作,获得每一炮对应的所述单炮偏移剖面;将所有炮的所述单炮偏移剖面进行叠加,获得初始偏移剖面;
所述初始多分量反偏移记录计算模块用于计算所有炮的多分量反偏移记录,具体为:针对每一炮,在该炮对应炮点位置设置所述震源子波,基于各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用所述初始偏移剖面、所述的该炮的背景波场以及所述给定的背景模型,构建该炮的虚拟震源;基于带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波正向延拓,获得该炮的反偏移波场,对该所述炮的反偏移波场进行记录采样,获得该炮的多分量反偏移记录;对所有炮执行该操作,获得每一炮对应的所述多分量反偏移记录,所有炮的所述多分量反偏移记录构成了初始多分量反偏移记录;
所述反演迭代模块是用于利用最优化反演算法对偏移剖面进行更新,获得最终的偏移剖面;
本发明实施采用的技术方案还包括:所述反演迭代模块包括多分量残差记录计算单元、梯度剖面计算单元、下降方向剖面计算单元、扰动多分量反偏移记录计算单元、步长计算单元、偏移剖面及多分量反偏移记录更新计算单元、收敛条件判断单元:
所述分量残差记录计算单元用于利用多分量反偏移记录和观测的多分量地震记录计算多分量残差记录,具体为,设置当前迭代次数,针对每一炮,利用前一次迭代所得到的该炮的所述多分量反偏移记录和该炮的所述观测的多分量地震记录计算该炮的多分量残差记录,对所有炮执行该操作,获得每一炮对应的所述多分量残差记录;
所述梯度剖面计算单元用于计算本次迭代的梯度剖面,具体为,针对每一炮,首先对所述该炮的多分量残差记录进行预处理,之后基于所述各向同性介质纵横波分解的弹性波方程,对经过所述预处理的该炮的多分量残差记录进行波场逆时延拓获得该炮的每一个时刻的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于所述该炮的震源波场以及所述该炮的检波点波场应用所述成像条件,获得该炮的单炮梯度剖面;对所有炮执行该操作,获得每一炮对应的所述单炮梯度剖面;将所有炮的所述单炮梯度剖面进行叠加,获得本次迭代的梯度剖面;
所述下降方向剖面计算单元用于利用最优化反演算法,构建本次迭代的下降方向剖面;
所述扰动多分量反偏移记录计算单元用于计算每一炮对应的扰动多分量反偏移记录,具体为,针对每一炮,在该炮对应炮点位置设置所述震源子波,基于所述各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用获得的本次迭代的所述下降方向剖面、所述的该炮的背景波场以及所述给定的背景模型,构建该炮的虚拟震源;基于所述带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的反偏移波场,对该炮的所述反偏移波场进行记录采样,获得该炮的多分量扰动反偏移记录;对所有炮执行该操作,获得每一炮对应的所述多分量扰动反偏移记录,所有炮的所述多分量扰动反偏移记录构成了本次迭代的多分量扰动反偏移记录;
所述步长计算单元用于利用步长公式计算本次迭代的优化步长;
所述偏移剖面及多分量反偏移记录更新计算单元用于更新所述偏移剖面和所述多分量反偏移记录,具体为,利用所述优化步长、前一次迭代所得所述偏移剖面以及本次迭代所得的所述下降方向剖面更新所述偏移剖面;利用所述优化步长、前一次迭代所得所述多分量反偏移记录以及本次迭代所得的所述多分量扰动反偏移记录更新所述多分量反偏移记录;
利用所述优化步长及下降方向剖面,更新本次迭代的偏移剖面;利用所述优化步长及所述多分量扰动反偏移记录,更新本次迭代的多分量反偏移记录;
所述收敛条件判断单元用于判断本次迭代是否满足收敛停止标准,具体为,计算本次迭代的目标泛函值,判断本次迭代是否满足收敛标准,如果满足则输出最新的所述偏移剖面为最终的偏移剖面;否则重复此反演模块,直至获得所述最终偏移剖面。
本发明由于采取以上技术方案,其具有以下优点:1)本发明将反演的思想引入弹性波逆时偏移中,与常规弹性波逆时偏移相比,本发明可以获得高精度、高分辨率、高信噪比、振幅保真的偏移剖面;2)本发明所用成像条件自然地校正了横波极性,无需极性反转校正运算便可有效克服横波极性反转造成的偏移剖面相消干涉;3)本发明采用的成像方式,无需额外的纵横波分离运算,而且完整地保持了纵横波的矢量特性、振幅以及相位特性,有效地消除了纵横波之间串扰造成的偏移假象,极大地提高了成像的精度;4)本发明所得成像剖面物理意义非常明确,便于后期地质解译;5)本发明可以广泛用于油气勘探领域中,特别是对于复杂构造深部的成像效果更加明显。
附图说明
图1是本发明提供的多分量地震资料最小二乘逆时偏移成像方法流程示意图;
图2是本发明提供的二维陡倾角断层模型图;其中,图2(a)是拉梅常数λ模型;图2(b)是拉梅常数μ的模型;
图3是图2所示二维陡倾角断层模型的多炮叠加偏移剖面:其中,图3(a)是利用传统方法所得的水平分量剖面,图3(b)是利用传统方法所得的垂直分量剖面,图3(c)是利用传统方法所得的PP剖面,图3(d)是利用传统方法所得的PS剖面,图3(e)是利用本发明方法所得的初始PP剖面,图3(f)是利用本发明方法所得的初始PS剖面;
图4是图2所示二维陡倾角断层模型的多炮叠加偏移剖面:其中,图4(a)是利用本发明方法所得的PP剖面;图4(b)是利用本发明方法所得的PS剖面;
图5是本发明提供的切除直达波后的单炮地震记录:其中,图5(a)是观测记录水平分量;图5(b)是观测记录垂直分量;图5(c)是利用本发明方法所得的反偏移记录水平分量;图5(d)是利用本发明方法所得的反偏移记录垂直分量;
图6是本发明提供的Marmousi-ii模型:其中,图6(a)是拉梅常数λ模型;图6(b)是拉梅常数μ模型;
图7是图6所示Marmousi-ii模型的多炮叠加偏移剖面:其中,图7(a)是利用传统方法所得的水平分量剖面;图7(b)是利用传统方法所得的垂直分量剖面;图7(c)是利用传统方法所得的PP剖面;图7(d)是利用传统方法所得的PS剖面;图7(e)是利用本发明方法所得的PP剖面;图7(f)是利用本发明方法所得的PS剖面;
图8是本发明提供的SEG/EAGE Salt模型:其中,图8(a)是纵波速度模型;图8(b)是横波速度模型;
图9是图8所示SEG/EAGE Salt模型的多炮叠加偏移剖面:其中,图9(a)是利用传统方法所得的PP剖面;图9(b)是利用传统方法所得的PS剖面(经过极性反转校正);图9(c)是利用本发明方法所得的PP剖面;图9(d)是利用本发明方法所得的PS剖面;
图10是本发明提供的多分量地震资料最小二乘逆时偏移成像系统结构示意图。
具体实施方式
为了是本发明的目的,技术方案及优点更加清楚明白,一下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅用于解释本发明,并不用于限定本发明。
请参阅图1,是本发明提供的多分量地震资料最小二乘逆时偏移成像方法流程示意图。本发明实施例的多分量地震资料最小二乘逆时偏移成像方法包括以下步骤:
步骤S100:读取预设参数、给定的背景模型及观测的多炮多分量地震记录D=(Dx,Dy,Dz),确定进行偏移成像的所有多分量地震记录;
步骤S200:针对每一炮,在该炮对应炮点位置设置震源子波,基于各向同性介质纵横波分解的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的每一个时刻的震源波场,该震源 波场包含矢量纵波震源波场;对所有炮均执行该操作,获得每一炮对应的所述震源波场;
所述的各向同性介质纵横波分解的弹性波方程为:
所述方程(9)中,Vp=(Vxp,Vyp,Vzp)T表示矢量纵波速度波场;Vs=(Vxs,Vys,Vzs)T表示矢量横波速度波场;V=(Vx,Vy,Vz)T表示矢量速度波场;σ=(σxx,σyy,σzz,σxy,σyz,σzx)T表示矢量应力波场;λ和μ表示介质拉梅常数;ρ表示介质密度;t表示时间;x,y和z分别表示笛卡尔坐标的x,y和z三个方向;
步骤S300:针对每一炮,首先对该炮的所述观测的多分量地震记录进行预处理,之后基于所述各向同性介质纵横波分解的弹性波方程,对经过所述预处理的该炮的多分量地震记录进行波场逆时延拓,获得该炮的每一个时刻的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于该炮的所述震源波场以及该炮的所述检波点波场应用成像条件,获得该炮的单炮偏移剖面;对所有炮执行该操作,获得每一炮对应的所述单炮偏移剖面;将所有炮的所述单炮偏移剖 面进行叠加,获得初始偏移剖面,也即第0次迭代的偏移剖面
所述的对该炮的所述观测的多分量地震记录进行预处理具体为:
其中,t表示波的传播时间;
所述成像条件为:
所述方程(11)中,偏移剖面具体为:
所述方程(11)中,表示震源波场中的矢量纵波速度波场,表示检波点波场中的矢量纵波速度波场,表示检波点波场中的矢量横波速度波场;
步骤S400:针对每一炮,在该炮对应炮点位置设置所述震源子波,基于各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用步骤S100中所述的初始偏移剖面I0、所述的该炮的背景波场以及所述给定的背景模型,构建该炮的虚拟震源;基于带有虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波正向延拓,获得该炮的反偏移波场,对该炮的所述反偏移波场进行记录采样,获得该炮的多分量反偏移记录;对所有炮执行该操作,获得每一炮对应的所述多分量反偏移记录,所有炮的所述多分量反偏移记录构成了初始多分量反偏移记录,也即第0次迭代的多分量反偏移记录此时,设置当前的迭代次数为i=0;
在步骤S400中所述的各向同性介质纵横波耦合的弹性波方程为:
在步骤S400中所述带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程可表示为:
所述方程(12)中,虚拟震源矢量F=(Fxx,Fyy,Fzz,Fxy,Fyz,Fzx)T可以表示为:
所述方程(12)中,δV=(δVx,δVy,δVz)T表示反偏移矢量速度波场;δσ=(δσxx,δσyy,δσzz,δσxy,δσyz,δσzx)T表示反偏移矢量应力波场;
所述方程(13)中,
步骤S500:设置当前的迭代次数i=i+1,基于第i-1次迭代所得到的所述偏移剖面和第i-1次迭代所得到的所述多分量反偏移记录进行反演迭代;
在第i次迭代中,针对每一炮,利用第i-1次迭代所得到的该炮的所述多分量反偏移记录di-1和该炮的所述观测的多分量地震记录D计算该炮的多分量残差记录,对所有炮执行该操作,获得每一炮对应的所述多分量残差记录
所述的利用第i-1次迭代所得到的该炮的所述多分量反偏移记录di-1和该炮的所述观测的多分量地震记录D计算该炮的多分量残差记录具体为:
步骤S600:针对每一炮,首先对该炮的所述多分量残差记录进行预处理,之后基于所述各向同性介质纵横波分解的弹性波方程,对经过所述预处理的该炮的多分量残差记录进行波场逆时延拓获得该炮的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于该炮的所述震源波场以及该炮的所述检波点波场应用所述成像条件,获得该炮的单炮梯度剖面;对所有炮执行该操作,获得每一炮对应的所述单炮梯度剖面;将所有炮的所述单炮梯度剖面进行叠加,获得本次迭代的梯度剖面,也即第i次迭代的梯度剖面
所述的对该炮的所述多分量残差记录进行预处理具体为:
所述对于该炮的所述震源波场以及该炮的所述检波点波场应用所述成像条件,具体为基于所述方程(11),其中,
步骤S700:利用最优化反演算法,构建第i次迭代的下降方向剖面ri=(ri pp,ri ps);
步骤S800:针对每一炮,在该炮对应炮点位置设置所述震源子波,基于所述各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用步骤S700中所获得的第i次迭代的所述下降方向剖面ri、所述的该炮的背景波场以及所述给定的背景模型,构建该炮的虚拟震源;基于所述的各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型以及所述的该炮的虚拟震源对该炮点进行波场正向延拓,获得该炮的反偏移波场,对该炮的所述反偏移波场进行记录采样,获得该炮的多分量扰动反偏移记录;对所有炮执行该操作,获得每一炮对应的所述多分量扰动反偏移记录,所有炮的所述多分量扰动反偏移记录构成了第i次迭代的多分量扰动反偏移记录 再利用步长公式计算本次迭代的优化步长αi;
所述利用步骤S700中所获得的第i次迭代的所述下降方向剖面ri、所述的该炮的背景波场以及所述给定的背景模型,构建该炮的虚拟震源,具体为基于所述方程(13),其中,m1=ri pp,m2=ri ps;
所述的步长公式为
所述方程(4)中,求和变量k包含笛卡尔坐标的x,y和z三个方向;xs表示震源点位置,xr表示检波点位置;
步骤S900:利用步骤S800得到的所述优化步长αi及步骤S700得到的所述下降方向剖面ri,更新第i次迭代的所述偏移剖面Ii=Ii-1+αiri;利用所述优化步长及步骤S800得到的所述多分量扰动反偏移记录δdi,更新第i次迭代的所述多分量反偏移记录di=di-1+αiδdi;
所述更新第i次迭代的所述偏移剖面Ii=Ii-1+αiri具体为:
所述更新第i次迭代的所述多分量反偏移记录di=di-1+αiδdi具体为:
步骤S1000:第i次迭代完后,计算第i次迭代的目标泛函值fi,判断当前迭代是否满足收敛标准,如果满足则输出最新的所述偏移剖面为最终的偏移剖面I=(Ipp,Ips);否则重复步骤S500至S1000,直至获得所述最终偏移剖面;
所述计算第i次迭代的目标泛函值fi具体为:
所述方程(7)中,求和变量k包含笛卡尔坐标的x,y和z三个方向;xs表示震源点位置,xr表示检波点位置;
所述的收敛标准具体为:
所述方程(8)中,theshold表示迭代停止的阈值标准,通常选取0.00001为最佳。
本发明涉及一种多分量地震资料最小二乘逆时偏移成像方法及系统,本发明是在弹性波逆时偏移方法的基础上进行了改进,特征是直接以多分量地震资料为输入,在反演的框架下,通过在偏移的不同步骤采用不同的波场延拓算子以及采用新的成像条件,从而实现基于反演的弹性波逆时偏移。本发明将反演的思想引入弹性波逆时偏移中,与常规弹性波逆时偏移相比,本发明可以获得高精度、高分辨率、高信噪比、振幅保真的叠前深度偏移剖面。本发明有效克服了横波极性反转造成的偏移剖面相消干涉,而且在完整地保持纵横波的矢量特性、振幅以及相位特性的同时,有效地消除了纵横波之间串扰造成的偏移假象,极大地提高了成像的精度。本发明可以应用到各种复杂介质模型的多分量地震资料偏移中,它可以获得模型不同参数的各自构造,且成像剖面物理意义非常明确,便于后期地质解译。
相应的,本发明提供一种多分量地震资料最小二乘逆时偏移成像系统,如图10所示,该系统包括:
初始化模块、震源波场计算模块、初始偏移剖面计算模块、初始多分量反偏移记录计算模块和反演迭代模块;
所述初始化模块用于读取预设参数、给定的背景模型及观测的多炮多分量地震记录,确定进行偏移成像的所有多分量地震记录;
所述震源波场计算模块用于利用各向同性介质纵横波分解的弹性波方程进行波场正向延拓获得每一炮对应的每一个时刻的震源波场;
所述初始偏移剖面计算模拟用于计算所有炮的叠加偏移剖面,具体为针对每一炮,首先对所述该炮的观测的多分量地震记录进行预处理,之后基于所述各向同性介质纵横波分解的弹性波方程,对经过所述预 处理的该炮的多分量地震记录进行波场逆时延拓,获得该炮的每一个时刻的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于所述该炮的震源波场以及所述该炮的检波点波场应用成像条件,获得该炮的单炮偏移剖面;对所有炮执行该操作,获得每一炮对应的所述单炮偏移剖面;将所有炮的所述单炮偏移剖面进行叠加,获得初始偏移剖面;
所述初始多分量反偏移记录计算模块用于计算所有炮的多分量反偏移记录,具体为:针对每一炮,在该炮对应炮点位置设置所述震源子波,基于各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用所述初始偏移剖面、所述的该炮的背景波场以及所述给定的背景模型,构建该炮的虚拟震源;基于带有虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波正向延拓,获得该炮的反偏移波场,对该所述炮的反偏移波场进行记录采样,获得该炮的多分量反偏移记录;对所有炮执行该操作,获得每一炮对应的所述多分量反偏移记录,所有炮的所述多分量反偏移记录构成了初始多分量反偏移记录;
所述反演迭代模块是用于利用最优化反演算法对偏移剖面进行更新,获得最终的偏移剖面;
本发明实施采用的技术方案还包括:所述反演迭代模块包括多分量残差记录计算单元、梯度剖面计算单元、下降方向剖面计算单元、扰动多分量反偏移记录计算单元、步长计算单元、偏移剖面及多分量反偏移记录更新计算单元、收敛条件判断单元:
所述分量残差记录计算单元用于利用多分量反偏移记录和观测的多分量地震记录计算多分量残差记录,具体为,设置当前迭代次数,针对每一炮,利用前一次迭代所得到的该炮的所述多分量反偏移记录和该炮的所述观测的多分量地震记录计算该炮的多分量残差记录,对所有炮执行该操作,获得每一炮对应的所述多分量残差记录;
所述梯度剖面计算单元用于计算本次迭代的梯度剖面,具体为,针对每一炮,首先对所述该炮的多分量残差记录进行预处理,之后基于所述各向同性介质纵横波分解的弹性波方程,对经过所述预处理的该炮的多分量残差记录进行波场逆时延拓获得该炮的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于所述该炮的震源波场以及所述该炮的检波点波场应用所述成像条件,获得该炮的单炮梯度剖面;对所有炮执行该操作,获得每一炮对应的所述单炮梯度剖面;将所有炮的所述单炮梯度剖面进行叠加,获得本次迭代的梯度剖面;
所述下降方向剖面计算单元用于利用最优化反演算法,构建本次迭代的下降方向剖面;
所述扰动多分量反偏移记录计算单元用于计算每一炮对应的扰动多分量反偏移记录,具体为,针对每一炮,在该炮对应炮点位置设置所述震源子波,基于所述各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用获得的本次迭代的所述下降方向剖面、所述的该炮的背景波场以及所述给定的背景模型,构建该炮的虚拟震源;基于所述带有虚拟震源 的各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的反偏移波场,对该炮的所述反偏移波场进行记录采样,获得该炮的多分量扰动反偏移记录;对所有炮执行该操作,获得每一炮对应的所述多分量扰动反偏移记录,所有炮的所述多分量扰动反偏移记录构成了本次迭代的多分量扰动反偏移记录;
所述步长计算单元用于利用步长公式计算本次迭代的优化步长;
所述偏移剖面及多分量反偏移记录更新计算单元用于更新所述偏移剖面和所述多分量反偏移记录,具体为,利用所述优化步长、前一次迭代所得所述偏移剖面以及本次迭代所得的所述下降方向剖面更新所述偏移剖面;利用所述优化步长、前一次迭代所得所述多分量反偏移记录以及本次迭代所得的所述多分量扰动反偏移记录更新所述多分量反偏移记录;
利用所述优化步长及下降方向剖面,更新本次迭代的偏移剖面;利用所述优化步长及所述多分量扰动反偏移记录,更新本次迭代的多分量反偏移记录;
所述收敛条件判断单元用于判断本次迭代是否满足收敛停止标准,具体为,计算本次迭代的目标泛函值,判断本次迭代是否满足收敛标准,如果满足则输出最新的所述偏移剖面为最终的偏移剖面;否则重复此反演模块,直至获得所述最终偏移剖面。
为进一步说明本发明的可行性和有效性,下面举三个实例:
实例1:
图2是二维陡倾角断层模型图;其中,图2(a)是拉梅常数λ模型;图2(b)是拉梅常数μ的模型。在此模型上设置49个爆炸震源,震源子波设定为雷克子波,主频为15赫兹,起始震源点的位于(150m,100m)处,炮间隔为100m。采用中间放炮两边接收观测系统,单边最大偏移距2300m,最小偏移距为150m,道间距为10m。是图3是图2所示二维陡倾角断层模型的多炮叠加偏移剖面:其中,图3(a)是利用传统方法所得的水平分量剖面,图3(b)是利用传统方法所得的垂直分量剖面,图3(c)是利用传统方法所得的PP剖面,图3(d)是利用传统方法所得的PS剖面,图3(e)是利用本发明方法所得的初始PP剖面,图3(f)是利用本发明方法所得的初始PS剖面。从图3(a)和图3(b)中可以看出,两个剖面均存在比较明显的纵横波串扰噪声,且每个剖面均包含两个模型参数的构造,剖面的分辨率较低,振幅不均衡。从图3(c)和图3(d)中可以看出,尽管两个剖面纵横波串扰假象被压制,横波极性反转造成了PS剖面出现严重的相消干涉,同时两个剖面的分辨率和振幅均衡性较差,每个剖面均包含两个模型参数的构造。从图3(e)和图3(f)中可以看出,本发明方法所得的初始剖面较好地压制了纵横波之间的串扰假象,横波极性反转也被自动校正,PP剖面只包含拉梅常数λ的构造,PS剖面只包含拉梅常数μ的构造,但是两个剖面的分辨率比较低,振幅也不均衡,图4是图2所示二维陡倾角断层模型的多炮叠加偏移剖面:其中,图4(a)是利用本发明方法所得的PP剖面;图4(b)是利用本发明方法所得的PS剖面。从图4(a)和图4(b)中可以看出,本发明方法所得的最终剖面具有很高的精度、分辨率和信噪比,而且振幅是均衡性很好。对比图3和图4可以证明,本发明方法能够获得高质量的偏移剖面,同时也证明了本发明方法的可行性及有效性。图5是切除直达波后的单炮地震记录:其中,图5(a)是观测记录水平分量;图5(b)是观测记录垂直分量;图5(c)是利用本发明方法所得的反偏移记录水平分量;图5(d)是利用本发明方法所得的反偏移记录垂直分量。从图5中可以看出,本发明方法所得反偏移记录与观测记录匹配非常好,这间接地验证了本发明方法的有效性。
实例2:
图6是Marmousi-ii模型:其中,图6(a)是拉梅常数λ模型;图6(b)是拉梅常数μ模型。该模型是验证各种偏移方法成像效果的国际标准模型之一。在此模型上设置47个爆炸震源,震源子波设定为雷克子波,主频为15赫兹,起始震源点的位于(100m,100m)处,炮间隔为200m。采用中间放炮两边接收观测系统,单边最大偏移距3500m,最小偏移距为100m,道间距为10m。图7是图6所示Marmousi-ii模型的多炮叠加偏移剖面:其中,图7(a)是利用传统方法所得的水平分量剖面;图7(b)是利用传统方法所得的垂直分量剖面;图7(c)是利用传统方法所得的PP剖面;图7(d)是利用传统方法所得的PS剖面;图7(e)是利用本发明方法所得的PP剖面;图7(f)是利用本发明方法所得的PS剖面。从图7(a)和图7(b)中可以看出,传统的笛卡尔分量成像条件不能压制纵横波之间的串扰,造成成像结果存在严重的串扰假象。从图7(c)和图7(d)中可以看出,传统方法所得的PP剖面和PS剖面在压制纵横波串扰方面相比于传统的笛卡尔分量成像条件要好,但纵横波串扰假象并没有被完全压制,此外,横波极性反转造成了PS剖面存在严重的同相轴干涉相消,使得剖面质量大幅度降低,严重影响后续的解释工作。由图7也可以看出,两种传统方法的成像剖面振幅严重失衡,深部剖面相对振幅并不保真。与其相比,本发明方法所得成像剖面可见图7(e)和图7(f),从中可以看出,本发明方法成像效果更好,纵横波串扰几乎被完全消除,横波极性反转自动得到校正,偏移剖面的分辨率和精度更高,信噪比也更好,振幅均衡性更佳。
实例3:
图8是SEG/EAGE Salt模型:其中,图8(a)是纵波速度模型;图8(b)是横波速度模型。该模型是验证各种偏移方法成像效果的国际标准盐丘模型之一。在此模型上设置39个爆炸震源,震源子波设定为雷克子波,主频为12赫兹,起始震源点的位于(170m,100m)处,炮间隔为150m。采用中间放炮两边接收观测系统,单边最大偏移距4800m,最小偏移距为170m,道间距为10m。图9是图8所示SEG/EAGE Salt模型的多炮叠加偏移剖面:其中,图9(a)是利用传统方法所得的PP剖面;图9(b)是利用传统方法所得的PS剖面(经过极性反转校正);图9(c)是利用本发明方法所得的PP剖面;图9(d)是利用本发明方法所得的PS剖面。从图9中可以看出,本发明方法所得偏移剖面具有更高的分辨率、更好的信噪比、精度也更高,此外我们也可以清晰地看出本发明方法所得剖面具有很高的振幅保真性。
以上描述和实施例仅用于说明本发明,凡是在本发明技术方案的基础上进行的等同变化和改进,均应包含在本发明的保护范围之内。
Claims (7)
1.一种多分量地震资料最小二乘逆时偏移成像方法,其特征在于包括以下步骤:
(a):读取预设参数、给定的背景模型及观测的多炮多分量地震记录D=(Dx,Dy,Dz),确定进行偏移成像的多分量地震记录;
(b):针对每一炮,在该炮对应炮点位置设置震源子波,基于各向同性介质纵横波分解的弹性波方程,利用上述给定的背景模型对该炮点进行波场正向延拓,获得该炮的每一个时刻的震源波场,该震源波场包含矢量纵波震源波场,从而获得每一炮对应的震源波场;
(c):针对每一炮,首先对与该炮对应的观测的多分量地震记录进行预处理,之后同样基于各向同性介质纵横波分解的弹性波方程,对该炮的经过预处理的多分量地震记录进行波场逆时延拓,获得该炮的每一个时刻的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于该炮的震源波场以及检波点波场应用成像条件,获得该炮的单炮偏移剖面;进而获得每一炮对应的单炮偏移剖面;将所有炮的单炮偏移剖面进行叠加,获得初始偏移剖面,也即第0次迭代的偏移剖面
所述的预处理具体为:
其中,t表示波的传播时间;
(d):针对每一炮,在该炮对应炮点位置设置震源子波,基于各向同性介质纵横波耦合的弹性波方程。利用(a)中给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用(c)中所述的初始偏移剖面I0、所述的该炮的背景波场以及(a)中给定的背景模型,构建该炮的虚拟震源;基于带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用给定的背景模型对该炮点进行波正向延拓,获得该炮的反偏移波场,对该炮的反偏移波场进行记录采样,获得该炮的多分量反偏移记录;进而获得每一炮对应的多分量反偏移记录,所有炮的多分量反偏移记录构成了初始多分量反偏移记录,也即第0次迭代的多分量反偏移记录此时,设置当前的迭代次数为i=0;
(e):设置当前的迭代次数i=i+1,基于第i-1次迭代所得到的所述偏移剖面和第i-1次迭代所得到的所述多分量反偏移记录进行反演迭代;
在第i次迭代中,针对每一炮,利用第i-1次迭代所得到的该炮的多分量反偏移记录di-1和与该炮对应的观测的多分量地震记录计算该炮的多分量残差记录,进而获得所有炮的多分量残差记录
所述的利用第i-1次迭代所得到的该炮的多分量反偏移记录di-1和与该炮对应的观测的多分量地震记录计算该炮的多分量残差记录具体为:
(f):针对每一炮,首先对该炮的多分量残差记录进行预处理,之后同样基于各向同性介质纵横波分解的弹性波方程,对经过预处理的多分量残差记录进行波场逆时延拓获得该炮的每一个时刻的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于该炮的检波点波场以及该炮的由步骤(b)得到的震源波场应用所述成像条件,获得该炮的单炮梯度剖面;进而获得每一炮对应的单炮梯度剖面;将所有炮的单炮梯度剖面进行叠加,获得本次迭代的梯度剖面,也即第i次迭代的梯度剖面
所述的对该炮的多分量残差记录进行预处理具体为:
(g):利用最优化反演算法,构建第i次迭代的下降方向剖面
(h):针对每一炮,在该炮对应炮点位置设置震源子波,同样基于各向同性介质纵横波耦合的弹性波方程,利用(a)中给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用步骤(g)中所获得的第i次迭代的下降方向剖面ri、该炮的背景波场以及给定的背景模型,构建该炮的虚拟震源;基于带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用给定的背景模型对该炮点进行波场正向延拓,获得该炮的反偏移波场,对该炮的反偏移波场进行记录采样,获得该炮的多分量扰动反偏移记录;进而获得每一炮对应的多分量扰动反偏移记录,所有炮的多分量扰动反偏移记录构成了第i次迭代的多分量扰动反偏移记录再利用步长公式计算本次迭代的优化步长αi;
所述的步长公式为
所述方程(4)中,求和变量k包含笛卡尔坐标的x,y和z三个方向;xs表示震源点位置,xr表示检波点位置;
(i):利用步骤(h)得到的优化步长αi及步骤(g)得到的下降方向剖面ri,更新第i次迭代的偏移剖面Ii=Ii-1+αiri;利用优化步长及步骤(h)得到的多分量扰动反偏移记录δdi,更新第i次迭代的多分量反偏移记录di=di-1+αiδdi;
所述更新第i次迭代的偏移剖面Ii=Ii-1+αiri具体为:
所述更新第i次迭代的多分量反偏移记录di=di-1+αiδdi具体为:
(j):第i次迭代完后,计算第i次迭代的目标泛函值fi,判断当前迭代是否满足收敛标准,如果满足则输出最新的偏移剖面为最终的偏移剖面I=(Ipp,Ips);否则重复步骤(e)至(i),直至获得最终偏移剖面;
所述计算第i次迭代的目标泛函值fi具体为:
所述方程(7)中,求和变量k包含笛卡尔坐标的x,y和z三个方向;xs表示震源点位置,xr表示检波点位置;
所述的收敛标准具体为:
所述方程(8)中,theshold表示迭代停止的阈值标准,通常选取0.00001。
2.根据权利要求1所述的一种多分量地震资料最小二乘逆时偏移成像方法,其特征在于步骤(b),步骤(c)和步骤(f)中所述的各向同性介质纵横波分解的弹性波方程为:
所述方程(9)中,Vp=(Vxp,Vyp,Vzp)T表示矢量纵波速度波场;Vs=(Vxs,Vys,Vzs)T表示矢量横波速度波场;V=Vp+Vs=(Vx,Vy,Vz)T表示矢量速度波场;σ=(σxx,σyy,σzz,σxy,σyz,σzx)T表示矢量应力波场;λ和μ表示介质拉梅常数;ρ表示介质密度;t表示时间;x,y和z分别表示笛卡尔坐标的x,y和z三个方向。
3.根据权利要求1或2所述的一种多分量地震资料最小二乘逆时偏移成像方法,其特征在于步骤(d)和步骤(h)中所述的各向同性介质纵横波耦合的弹性波方程为:
4.根据权利要求1或2或3所述的一种多分量地震资料最小二乘逆时偏移成像方法,其特征在于步骤(c)和步骤(f)中所述成像条件为:
所述方程(11)中,在步骤(c)中,在步骤(f)中,
所述方程(11)中,表示震源波场中的矢量纵波速度波场,表示检波点波场中的矢量纵波速度波场,表示检波点波场中的矢量横波速度波场。
5.根据权利要求1或2或3或4所述的一种多分量地震资料最小二乘逆时偏移成像方法,其特征在于步骤(d)和步骤(h)中所述基于带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用给定的背景模型对该炮点进行波正向延拓,获得该炮的矢量反偏移波场,此处的带有所述虚拟震源的各向同 性介质纵横波耦合的弹性波方程表示为:
所述方程(12)中,虚拟震源矢量F=(Fxx,Fyy,Fzz,Fxy,Fyz,Fzx)T表示为:
所述方程(12)中,δV=(δVx,δVy,δVz)T表示反偏移矢量速度波场;δσ=(δσxx,δσyy,δσzz,δσxy,δσyz,δσzx)T表示反偏移矢量应力波场;
所述方程(13)中,在步骤(d)中,在步骤(h)中,m1=ri PP,m2=ri ps。
6.一种多分量地震资料最小二乘逆时偏移成像系统,其特征在于包括初始化模块、震源波场计算模块、初始偏移剖面计算模块、初始多分量反偏移记录计算模块和反演迭代模块;
所述初始化模块用于读取预设参数、给定的背景模型及观测的多炮多分量地震记录,确定进行偏移成像的所有多分量地震记录;
所述震源波场计算模块用于利用各向同性介质纵横波分解的弹性波方程进行波场正向延拓获得每一炮对应的每一个时刻的震源波场;
所述初始偏移剖面计算模拟用于计算所有炮的叠加偏移剖面,具体为针对每一炮,首先对所述该炮的观测的多分量地震记录进行预处理,之后基于所述各向同性介质纵横波分解的弹性波方程,对经过所述预处理的该炮的多分量地震记录进行波场逆时延拓,获得该炮的每一个时刻的检波点波场,该检波点波场包含矢量纵波检波点波场和矢量横波检波点波场;在相同的时刻,对于所述该炮的震源波场以及所述该炮的检波点波场应用成像条件,获得该炮的单炮偏移剖面;对所有炮执行该操作,获得每一炮对应的所述单炮偏移剖面;将所有炮的所述单炮偏移剖面进行叠加,获得初始偏移剖面;
所述初始多分量反偏移记录计算模块用于计算所有炮的多分量反偏移记录,具体为:针对每一炮,在该炮对应炮点位置设置所述震源子波,基于各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用所述初始偏移剖面、所述的该炮的背景波场以及所述给定的背景模型,构建该炮的虚拟震源;基于带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波正向延拓,获得该炮的反偏移波场,对该所述炮的反偏移波场进行记录采样,获得该炮的多分量反偏移记录;对所有炮执行该操作,获得每一炮对应的所述多分量反偏移记录,所有炮的所述多分量反偏移记录构成了初始多分量反偏移记录;
所述反演迭代模块是用于利用最优化反演算法对偏移剖面进行更新,获得最终的偏移剖面。
7.根据权利要求6所述的一种多分量地震资料最小二乘逆时偏移成像系统,其特征在于:所述反演迭代模块包括多分量残差记录计算单元、梯度剖面计算单元、下降方向剖面计算单元、扰动多分量反偏移记录计算单元、步长计算单元、偏移剖面及多分量反偏移记录更新计算单元、收敛条件判断单元:
所述分量残差记录计算单元用于利用多分量反偏移记录和观测的多分量地震记录计算多分量残差记录,具体为,设置当前迭代次数,针对每一炮,利用前一次迭代所得到的该炮的所述多分量反偏移记录和该炮的所述观测的多分量地震记录计算该炮的多分量残差记录,对所有炮执行该操作,获得每一炮对应的所述多分量残差记录;
所述梯度剖面计算单元用于计算本次迭代的梯度剖面,具体为,针对每一炮,首先对所述该炮的多分量残差记录进行预处理,之后基于所述各向同性介质纵横波分解的弹性波方程,对经过所述预处理的该炮的多分量残差记录进行波场逆时延拓获得该炮的检波点波场,该检波点波场包含矢量纵波检波点波场和矢 量横波检波点波场;在相同的时刻,对于所述该炮的震源波场以及所述该炮的检波点波场应用所述成像条件,获得该炮的单炮梯度剖面;对所有炮执行该操作,获得每一炮对应的所述单炮梯度剖面;将所有炮的所述单炮梯度剖面进行叠加,获得本次迭代的梯度剖面;
所述下降方向剖面计算单元用于利用最优化反演算法,构建本次迭代的下降方向剖面;
所述扰动多分量反偏移记录计算单元用于计算每一炮对应的扰动多分量反偏移记录,具体为,针对每一炮,在该炮对应炮点位置设置所述震源子波,基于所述各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的背景波场;利用获得的本次迭代的所述下降方向剖面、所述的该炮的背景波场以及所述给定的背景模型,构建该炮的虚拟震源;基于所述的带有所述虚拟震源的各向同性介质纵横波耦合的弹性波方程,利用所述给定的背景模型对该炮点进行波场正向延拓,获得该炮的反偏移波场,对该炮的所述反偏移波场进行记录采样,获得该炮的多分量扰动反偏移记录;对所有炮执行该操作,获得每一炮对应的所述多分量扰动反偏移记录,所有炮的所述多分量扰动反偏移记录构成了本次迭代的多分量扰动反偏移记录;
所述步长计算单元用于利用步长公式计算本次迭代的优化步长;
所述偏移剖面及多分量反偏移记录更新计算单元用于更新所述偏移剖面和所述多分量反偏移记录,具体为,利用所述优化步长、前一次迭代所得所述偏移剖面以及本次迭代所得的所述下降方向剖面更新所述偏移剖面;利用所述优化步长、前一次迭代所得所述多分量反偏移记录以及本次迭代所得的所述多分量扰动反偏移记录更新所述多分量反偏移记录;
利用所述优化步长及下降方向剖面,更新本次迭代的偏移剖面;利用所述优化步长及所述多分量扰动反偏移记录,更新本次迭代的多分量反偏移记录;
所述收敛条件判断单元用于判断本次迭代是否满足收敛停止标准,具体为,计算本次迭代的目标泛函值,判断本次迭代是否满足收敛标准,如果满足则输出最新的所述偏移剖面为最终的偏移剖面;否则重复此反演模块,直至获得所述最终偏移剖面。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610520574.0A CN105974470B (zh) | 2016-07-04 | 2016-07-04 | 一种多分量地震资料最小二乘逆时偏移成像方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610520574.0A CN105974470B (zh) | 2016-07-04 | 2016-07-04 | 一种多分量地震资料最小二乘逆时偏移成像方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105974470A true CN105974470A (zh) | 2016-09-28 |
CN105974470B CN105974470B (zh) | 2017-06-16 |
Family
ID=56953691
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610520574.0A Active CN105974470B (zh) | 2016-07-04 | 2016-07-04 | 一种多分量地震资料最小二乘逆时偏移成像方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105974470B (zh) |
Cited By (16)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106842300A (zh) * | 2016-12-21 | 2017-06-13 | 中国石油大学(华东) | 一种高效率多分量地震资料真振幅偏移成像方法 |
CN106970416A (zh) * | 2017-03-17 | 2017-07-21 | 中国地质科学院地球物理地球化学勘查研究所 | 基于波场分离的弹性波最小二乘逆时偏移系统及方法 |
CN107589443A (zh) * | 2017-08-16 | 2018-01-16 | 东北石油大学 | 基于弹性波最小二乘逆时偏移成像的方法及系统 |
CN108241173A (zh) * | 2017-12-28 | 2018-07-03 | 中国石油大学(华东) | 一种地震资料偏移成像方法及系统 |
CN108802813A (zh) * | 2018-06-13 | 2018-11-13 | 中国石油大学(华东) | 一种多分量地震资料偏移成像方法及系统 |
CN109188519A (zh) * | 2018-09-17 | 2019-01-11 | 中国石油大学(华东) | 一种极坐标下的弹性波纵横波速度反演系统及方法 |
CN109471171A (zh) * | 2018-09-21 | 2019-03-15 | 中国石油天然气集团有限公司 | 一种混叠地震数据分离的方法、装置及系统 |
WO2019242045A1 (zh) * | 2018-06-21 | 2019-12-26 | 成都启泰智联信息科技有限公司 | 一种虚拟震源二维波前构建地震波走时计算方法 |
CN111158049A (zh) * | 2019-12-27 | 2020-05-15 | 同济大学 | 一种基于散射积分法的地震逆时偏移成像方法 |
US10788597B2 (en) | 2017-12-11 | 2020-09-29 | Saudi Arabian Oil Company | Generating a reflectivity model of subsurface structures |
CN112462428A (zh) * | 2020-11-13 | 2021-03-09 | 中国地质科学院 | 一种多分量地震资料偏移成像方法及系统 |
CN112904426A (zh) * | 2021-03-27 | 2021-06-04 | 中国石油大学(华东) | 一种解耦弹性波逆时偏移方法、系统及应用 |
CN113406698A (zh) * | 2021-05-24 | 2021-09-17 | 中国石油大学(华东) | 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法 |
US11320557B2 (en) | 2020-03-30 | 2022-05-03 | Saudi Arabian Oil Company | Post-stack time domain image with broadened spectrum |
CN114814944A (zh) * | 2022-03-15 | 2022-07-29 | 长江大学 | 基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法 |
CN116413801A (zh) * | 2023-02-13 | 2023-07-11 | 中国石油大学(北京) | 一种各向异性介质弹性波高精度成像方法和系统 |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109557582B (zh) * | 2018-12-17 | 2019-11-29 | 中国石油大学(华东) | 一种二维多分量地震资料偏移成像方法及系统 |
Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20010004030A1 (en) * | 1999-12-13 | 2001-06-21 | Ryuji Yano | Seat belt control device |
US20040196738A1 (en) * | 2003-04-07 | 2004-10-07 | Hillel Tal-Ezer | Wave migration by a krylov space expansion of the square root exponent operator, for use in seismic imaging |
WO2012141805A2 (en) * | 2011-04-13 | 2012-10-18 | Chevron U.S.A. Inc. | Stable shot illumination compensation |
CN102879816A (zh) * | 2012-07-17 | 2013-01-16 | 中国科学院地质与地球物理研究所 | 一种地震多次波偏移方法 |
GB2503326A (en) * | 2012-04-19 | 2013-12-25 | Cggveritas Services Sa | Seismic data processing including predicting multiples using a reverse time demigration |
US20140301158A1 (en) * | 2013-04-03 | 2014-10-09 | Cgg Services Sa | Device and method for stable least-squares reverse time migration |
CN104360381A (zh) * | 2014-10-20 | 2015-02-18 | 李闯 | 一种地震资料的偏移成像处理方法 |
CN104391327A (zh) * | 2014-12-04 | 2015-03-04 | 中国海洋石油总公司 | 一种海上斜井井间地震叠前逆时深度偏移成像方法 |
CN105005076A (zh) * | 2015-06-02 | 2015-10-28 | 中国海洋石油总公司 | 基于最小二乘梯度更新速度模型的地震波全波形反演方法 |
-
2016
- 2016-07-04 CN CN201610520574.0A patent/CN105974470B/zh active Active
Patent Citations (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20010004030A1 (en) * | 1999-12-13 | 2001-06-21 | Ryuji Yano | Seat belt control device |
US20040196738A1 (en) * | 2003-04-07 | 2004-10-07 | Hillel Tal-Ezer | Wave migration by a krylov space expansion of the square root exponent operator, for use in seismic imaging |
US6819628B2 (en) * | 2003-04-07 | 2004-11-16 | Paradigm Geophysical (Luxembourg) S.A.R.L. | Wave migration by a krylov space expansion of the square root exponent operator, for use in seismic imaging |
WO2012141805A2 (en) * | 2011-04-13 | 2012-10-18 | Chevron U.S.A. Inc. | Stable shot illumination compensation |
GB2503326A (en) * | 2012-04-19 | 2013-12-25 | Cggveritas Services Sa | Seismic data processing including predicting multiples using a reverse time demigration |
CN102879816A (zh) * | 2012-07-17 | 2013-01-16 | 中国科学院地质与地球物理研究所 | 一种地震多次波偏移方法 |
US20140301158A1 (en) * | 2013-04-03 | 2014-10-09 | Cgg Services Sa | Device and method for stable least-squares reverse time migration |
CN104360381A (zh) * | 2014-10-20 | 2015-02-18 | 李闯 | 一种地震资料的偏移成像处理方法 |
CN104391327A (zh) * | 2014-12-04 | 2015-03-04 | 中国海洋石油总公司 | 一种海上斜井井间地震叠前逆时深度偏移成像方法 |
CN105005076A (zh) * | 2015-06-02 | 2015-10-28 | 中国海洋石油总公司 | 基于最小二乘梯度更新速度模型的地震波全波形反演方法 |
Non-Patent Citations (2)
Title |
---|
刘玉金等人: "扩展成像条件下的最小二乘逆时偏移", 《地球物理学报》 * |
郭书娟等人: "最小二乘逆时偏移成像方法的实现与应用研究", 《石油物探》 * |
Cited By (23)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106842300B (zh) * | 2016-12-21 | 2018-10-30 | 中国石油大学(华东) | 一种高效率多分量地震资料真振幅偏移成像方法 |
CN106842300A (zh) * | 2016-12-21 | 2017-06-13 | 中国石油大学(华东) | 一种高效率多分量地震资料真振幅偏移成像方法 |
CN106970416A (zh) * | 2017-03-17 | 2017-07-21 | 中国地质科学院地球物理地球化学勘查研究所 | 基于波场分离的弹性波最小二乘逆时偏移系统及方法 |
CN107589443A (zh) * | 2017-08-16 | 2018-01-16 | 东北石油大学 | 基于弹性波最小二乘逆时偏移成像的方法及系统 |
US10788597B2 (en) | 2017-12-11 | 2020-09-29 | Saudi Arabian Oil Company | Generating a reflectivity model of subsurface structures |
CN108241173A (zh) * | 2017-12-28 | 2018-07-03 | 中国石油大学(华东) | 一种地震资料偏移成像方法及系统 |
CN108802813A (zh) * | 2018-06-13 | 2018-11-13 | 中国石油大学(华东) | 一种多分量地震资料偏移成像方法及系统 |
CN108802813B (zh) * | 2018-06-13 | 2019-07-23 | 中国石油大学(华东) | 一种多分量地震资料偏移成像方法及系统 |
WO2019242045A1 (zh) * | 2018-06-21 | 2019-12-26 | 成都启泰智联信息科技有限公司 | 一种虚拟震源二维波前构建地震波走时计算方法 |
CN109188519A (zh) * | 2018-09-17 | 2019-01-11 | 中国石油大学(华东) | 一种极坐标下的弹性波纵横波速度反演系统及方法 |
CN109471171A (zh) * | 2018-09-21 | 2019-03-15 | 中国石油天然气集团有限公司 | 一种混叠地震数据分离的方法、装置及系统 |
CN111158049A (zh) * | 2019-12-27 | 2020-05-15 | 同济大学 | 一种基于散射积分法的地震逆时偏移成像方法 |
US11320557B2 (en) | 2020-03-30 | 2022-05-03 | Saudi Arabian Oil Company | Post-stack time domain image with broadened spectrum |
CN112462428A (zh) * | 2020-11-13 | 2021-03-09 | 中国地质科学院 | 一种多分量地震资料偏移成像方法及系统 |
CN112462428B (zh) * | 2020-11-13 | 2024-02-20 | 中国地质科学院 | 一种多分量地震资料偏移成像方法及系统 |
CN112904426A (zh) * | 2021-03-27 | 2021-06-04 | 中国石油大学(华东) | 一种解耦弹性波逆时偏移方法、系统及应用 |
CN112904426B (zh) * | 2021-03-27 | 2022-09-30 | 中国石油大学(华东) | 一种解耦弹性波逆时偏移方法、系统及应用 |
CN113406698A (zh) * | 2021-05-24 | 2021-09-17 | 中国石油大学(华东) | 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法 |
CN113406698B (zh) * | 2021-05-24 | 2023-04-21 | 中国石油大学(华东) | 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法 |
CN114814944A (zh) * | 2022-03-15 | 2022-07-29 | 长江大学 | 基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法 |
CN114814944B (zh) * | 2022-03-15 | 2022-10-21 | 长江大学 | 基于散度和旋度的弹性纵横波场分离和逆时偏移成像方法 |
CN116413801A (zh) * | 2023-02-13 | 2023-07-11 | 中国石油大学(北京) | 一种各向异性介质弹性波高精度成像方法和系统 |
CN116413801B (zh) * | 2023-02-13 | 2023-09-29 | 中国石油大学(北京) | 一种各向异性介质弹性波高精度成像方法和系统 |
Also Published As
Publication number | Publication date |
---|---|
CN105974470B (zh) | 2017-06-16 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105974470B (zh) | 一种多分量地震资料最小二乘逆时偏移成像方法及系统 | |
US8892410B2 (en) | Estimation of soil properties using waveforms of seismic surface waves | |
Zhang et al. | Robust source-independent elastic full-waveform inversion in the time domain | |
US6041018A (en) | Method for correcting amplitude and phase differences between time-lapse seismic surveys | |
EP2696217B1 (en) | Device and method for directional designature of seismic data | |
CN111221037B (zh) | 解耦弹性逆时偏移成像方法和装置 | |
Borisov et al. | Application of 2D full-waveform inversion on exploration land data | |
CN102156296A (zh) | 地震多分量联合弹性逆时偏移成像方法 | |
CN108241173B (zh) | 一种地震资料偏移成像方法及系统 | |
CN108139498A (zh) | 具有振幅保持的fwi模型域角度叠加 | |
CN109946741B (zh) | 一种TTI介质中纯qP波最小二乘逆时偏移成像方法 | |
Pecholcs et al. | A broadband full azimuth land seismic case study from Saudi Arabia using a 100,000 channel recording system at 6 terabytes per day: acquisition and processing lessons learned | |
CN106842300B (zh) | 一种高效率多分量地震资料真振幅偏移成像方法 | |
WO2010027682A2 (en) | Processing seismic data in common group-center gathers | |
Ibrahim | Separating simultaneous seismic sources using robust inversion of Radon and migration operators | |
US11635540B2 (en) | Methods and devices performing adaptive quadratic Wasserstein full-waveform inversion | |
Kamei et al. | Application of waveform tomography to a crooked-line 2D land seismic data set | |
Xin et al. | 3-D tomographic Q inversion for compensating frequency dependent attenuation and dispersion | |
CN110873893A (zh) | 深水自由表面多次波预测和压制方法及其系统 | |
Huang et al. | Elastic Gaussian beam migration method for rugged topography | |
Hou et al. | True amplitude imaging for incomplete seismic data | |
Brandsberg-Dahl et al. | Beam-wave imaging | |
Wu et al. | Attenuation compensation in multicomponent Gaussian beam prestack depth migration | |
Xie et al. | 3D prestack beam migration with compensation for frequency dependent absorption and dispersion | |
Hanafy et al. | Early arrival waveform inversion of shallow seismic land data |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant |