CN105911587B - 一种利用单程波算子的双程波叠前深度偏移方法 - Google Patents
一种利用单程波算子的双程波叠前深度偏移方法 Download PDFInfo
- Publication number
- CN105911587B CN105911587B CN201610257053.0A CN201610257053A CN105911587B CN 105911587 B CN105911587 B CN 105911587B CN 201610257053 A CN201610257053 A CN 201610257053A CN 105911587 B CN105911587 B CN 105911587B
- Authority
- CN
- China
- Prior art keywords
- wave
- migration
- operator
- fourier
- way
- 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
Links
- 238000013508 migration Methods 0.000 title claims abstract description 99
- 230000005012 migration Effects 0.000 title claims abstract description 99
- 238000000034 method Methods 0.000 title claims abstract description 86
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 75
- 238000003384 imaging method Methods 0.000 claims abstract description 26
- 230000009466 transformation Effects 0.000 claims abstract description 10
- 238000004364 calculation method Methods 0.000 claims description 13
- 238000012545 processing Methods 0.000 description 19
- 238000010586 diagram Methods 0.000 description 13
- 238000005516 engineering process Methods 0.000 description 11
- 238000004458 analytical method Methods 0.000 description 4
- 238000010276 construction Methods 0.000 description 4
- 238000009776 industrial production Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 238000013213 extrapolation Methods 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 2
- 241001515806 Stictis Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000013480 data collection Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 239000007789 gas Substances 0.000 description 1
- 239000003345 natural gas Substances 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 238000004451 qualitative analysis Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000007787 solid Substances 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/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
- 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/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- 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
-
- 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/63—Seismic attributes, e.g. amplitude, polarity, instant phase
-
- 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/64—Geostructures, e.g. in 3D data cubes
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
本发明涉及一种利用单程波算子的双程波叠前深度偏移方法,包括如下步骤:1,采用上下双检波器地震数据采集系统,形成检波器波场;2,提取地震子波信息,形成震源波场;3,读取单炮地震数据、速度模型和震源波场;4,对震源波场和检波器波场在时间域作傅里叶变换得到频率‑空间域波场数据;5,对频率‑空间域波场数据进行深度延拓计算;6,利用互相关成像原理或反射系数成像原理进行偏移成像;7,重复步骤4~6延拓计算到最大深度,保存偏移结果,计算下一炮地震数据;8,重复步骤3~7至最后一炮,输出偏移剖面。本发明,利用单程波算子进行全波动方程深度偏移,实现了快速、准确的双程波深度偏移算法,更有工业利用价值。
Description
技术领域
本发明涉及地震勘探中的偏移成像方法,具体说是一种利用单程波算子的双程波叠前深度偏移方法。适用于油气勘探。单程波算子指单程波延拓算子。
背景技术
地震勘探方法是目前石油、天然气及其他能源勘探的重要方法,而地震偏移(指地震偏移成像技术)作为现代地震勘探数据处理的关键环节,为地震解释和岩性反演提供必要的技术支撑。因此,发展一种不仅能准确反映地下构造形态,还能反映地层岩性变化的真振幅地震偏移技术(真振幅地震偏移成像技术),一直是人们关注的问题,特别是为后偏移处理中的诸如地震岩性解释、属性分析和AVO(Amplitude Versus Offset,振幅随偏移距的变化)分析提供较为可靠的振幅信息。
现有技术一:常规地震偏移处理技术(即现代常规地震数据采集技术)
在地震数据采集方面只接收地表(或海面)的波场值,然后利用该地表的波场值成像地下构造。
与地震数据采集方案配套的地震偏移成像技术有:
(1)克希霍夫积分偏移法;
(2)有限差分偏移法;
(3)频率波数域偏移法;
(4)逆时偏移法。
这些方法的共同特点是:依据波动方程利用地表的波场值成像地下构造。
由于波动方程含有深度方向的二阶偏导数,数学上要使该波动方程能够求解,必须已知地表的波场值和地表波场在深度方向的偏导数。而现有技术没有接收地表波场在深度方向的偏导数,因此,在数学上严格意义的解是不存在的。
为此,现有技术对波动方程进行各种近似处理,以便只利用地表的波场值对地下构造进行成像。正是这种对波动方程的近似处理导致成像振幅失真。成像振幅不能真实的反映目的层的岩性变化。从而不利于岩性勘探。
现有技术二:上下双检波器地震数据采集系统的全波动方程真振幅叠前深度偏移[1]
该发明是利用上下双检波器记录的波场信息,通过差分公式计算地表波场对深度的导数,根据地表的波场值(上层检波器在地表处记录的波场值)及地表波场对深度的导数两个边界信息,在深度域上求解全波动方程。在该发明中,偏移算法没有使用任何形式的近似,所以理论上该发明可以偏移计算各种反射波及实现地震数据的真振幅偏移。
在该发明中,在深度域上进行偏移时,不可避免的要滤除耗损波。该发明采用了谱投影算法来消除深度延拓中产生的耗损波。在利用谱投影算法时,对每一层延拓深度都必须计算频带内所有频率成分的谱投影,对于宽频带地震数据,该步骤会制约地震深度偏移处理的效率,不利于工业生产实际操作。
现有技术三:采用单程波深度偏移算法的单程波波动方程偏移技术
该技术是将波动方程分解为上行波方程和下行波方程,结合地表的波场值实现上下行波的波场延拓成像。单程波波动方程偏移技术是解决地震偏移成像的有力技术,经过数十年的发展,其方法和技术取得了快速发展,使其在速度横向变化剧烈的介质中能较为准确的模拟波场传播。在单程波深度偏移算法中,根据使用的单程波算子的不同,单程波深度偏移算法包括:频率-空间域有限差分法[2]、裂步傅里叶法[3]、傅里叶有限差分法[4]及广义屏法[5]等。
该技术的假设前提是全波动方程可以准确的分解为上下行波方程。因为在实际应用中很难满足上述假设条件,因此单程波波动方程对复杂构造成像和陡倾角构造成像效果并不理想。基于同样的原因,该技术对于真振幅成像也无能为力。
本发明涉及以下定义:
1.网格剖分间距:Δz,Δx,Δy
Δz为z方向网格剖分间距,也是延拓步长;
Δx为x方向网格剖分间距;
Δy为y方向网格剖分间距;
2.角频率:ω
3.地震波场值:
u(x,y,z,t),将其在时间域上做傅里叶变换后为频率域波场值其在频率-波数域为u(kx,ky,kz,ω);
为在z=0深度处的波场值u(x,y,z=0,t)关于时间做傅里叶变换后的结果;
为在z=0深度处的波场值u(x,y,z=0,t)关于x,y,t做傅里叶变换后的结果;
为在z=0深度处的观测系统采集波场值d(x,y,z=0,t)关于x,y,t做傅里叶变换后的结果;
为在z=Δz深度处的波场值u(x,y,z=Δz,t)关于时间做傅里叶变换后的结果;
为在z=Δz深度处的波场值u(x,y,z=Δz,t)关于x,y,t做傅里叶变换后的结果;
为在z=Δz深度处的观测系统采集波场值d(x,y,z=Δz,t)关于x,y,t做傅里叶变换后的结果;
为在z=2Δz深度处的波场值u(x,y,z=2Δz,t)关于x,y,t做傅里叶变换后的结果;
4.模型速度:v(x,y,z)
5.频率域地震波场值对z方向的二阶偏导数:
频率波数域地震波场值对z方向的二阶偏导数:
6.z方向的波数:
7.为虚数单位。
发明内容
针对现有技术中存在的缺陷,本发明的目的在于提供一种利用单程波算子的双程波叠前深度偏移方法,以克服二阶偏导数波动方程在深度域上无法数值求解的问题。
为达到以上目的,本发明采取的技术方案是:
一种利用单程波算子的双程波叠前深度偏移方法,其特征在于,包括如下步骤:
步骤1,通过上下双检波器地震数据采集系统获取双检波器地震数据,形成检波器波场;
步骤2,从双检波器地震数据中提取地震子波信息,形成震源波场;
步骤3,读取单炮地震数据、速度模型和震源波场;
步骤4,对震源波场和检波器波场在时间域作傅里叶变换得到频率-空间域波场数据;
步骤5,对频率-空间域波场数据进行深度延拓计算;
步骤6,利用互相关成像原理或反射系数成像原理进行偏移成像;
步骤7,重复步骤4~6延拓计算到最大深度,保存偏移结果,计算下一炮地震数据;
步骤8,重复步骤3~7至最后一炮,输出偏移剖面。
在上述技术方案的基础上,步骤5具体包括:
输入z和z+Δz的波场值和根据下式计算z+2Δz的波场值
式(5)中,为单程波算子。
在上述技术方案的基础上,采用不同的近似算法来计算kz,则得到不同的双检波器叠前深度偏移算法。
在上述技术方案的基础上,具体计算单程波算子中的kz时,采用裂步傅里叶方法、傅里叶有限差分方法或高阶广义屏方法;
分别形成:双检波器裂步傅里叶叠前深度偏移、双检波器傅里叶有限差分叠前深度偏移和双检波器高阶广义屏叠前深度偏移。
在上述技术方案的基础上,以裂步傅里叶方法来近似单程波算子时,其中kz可以表示为
其中:vr为介质参考速度;kx和ky分别表示x方向和y方向波数;ASSF表示利用裂步傅里叶方法近似的单程波算子;表示利用裂步傅里叶方法近似计算的kz。
在上述技术方案的基础上,以傅里叶有限差分方法来近似单程波算子时,其中kz可以表示为
其中:
AFFD表示利用傅里叶有限差分方法近似的单程波算子;表示利用傅里叶有限差分方法近似计算的kz;
和分别表示对x和y的二阶偏导数。
在上述技术方案的基础上,以高阶广义屏方法来近似单程波算子时,其中kz可以表示为
其中:
AGSP表示利用高阶广义屏方法近似的单程波算子;表示利用高阶广义屏方法近似计算的kz。
本发明所述的利用单程波算子的双程波叠前深度偏移方法,为适应工业生产中对快速偏移算法的要求,同时,为了保证偏移振幅的可靠性而提出,基于上下双检波器地震数据采集系统,通过上下双检波器地震数据采集系统提供的两个边界条件,利用单程波算子(单程波延拓算子)进行全波动方程叠前深度偏移,建立双程波叠前深度偏移方案,实现了波动方程在深度域的完全求解。
本发明所述的利用单程波算子的双程波叠前深度偏移方法,基于上下双检波器地震数据采集系统,通过在双程波深度偏移中成功的引入单程波算子解决了上述谱投影算法的缺点,利用了单程波算子的计算形式,并进一步发展了单程波偏移方法,在构造成像和真振幅成像上取得了单程波算子无法比拟的效果,实现了快速、准确的双程波深度偏移算法,更有工业利用价值。
本发明所述的利用单程波算子的双程波叠前深度偏移方法,利用常规的单程波算子构建了新的双程波叠前深度偏移算法,不仅能快速成像各种复杂的地下构造形态,也能提供在一定程度上反映地下构造岩性变化的真振幅反射信息,这将为地震后期的AVO(AVA)反演、属性分析、全波形反演及岩性物理提供理论基础,具有较高的工业生产应用价值。
本发明所述的利用单程波算子的双程波叠前深度偏移方法,是基于上下双检波器地震数据采集系统提出来的一种利用单程波算子的双程波叠前深度偏移算法(双程波叠前深度延拓算法),达到了利用全波动方程求解地震波场信息的目的,并为后续的地质解释人员提供可靠的地下构造和岩性信息,增加了地质解释和地震反演分析的可信度。
本发明有以下有益效果:
1.计算效率高。本发明结合了单程波算法高效计算的优点,实现了双程波叠前深度偏移地震数据的快速处理,这将为工业生产节约成本。
2.本发明是基于全波动方程的深度偏移算法,因此,理论上本发明可以偏移处理各种波场信息,利于构造成像。
3.相较于传统的单程波叠前深度偏移算法,本发明的算法具有更好的岩性识别性能,为后偏移处理提供可靠的资料。
附图说明
本发明有如下附图:
图1.Marmousi速度模型;
图2.a单程波广义屏叠前深度偏移剖面;
图2.b双检波器高阶广义屏叠前深度偏移剖面;
图3.a叠前逆时偏移剖面;
图3.b双检波器裂步傅里叶叠前深度偏移剖面;
图3.c双检波器傅里叶有限差分叠前深度偏移剖面;
图3.d双检波器高阶广义屏叠前深度偏移剖面;
图4.a单程波广义屏叠前深度偏移算法在Marmousi模型的A位置计算的反射系数与相同位置处理论反射系数对比图;
图4.b叠前逆时偏移算法在Marmousi模型的A位置计算的反射系数与相同位置处理论反射系数对比图;
图4.c双检波器裂步傅里叶叠前深度偏移算法在Marmousi模型的A位置计算的反射系数与相同位置处理论反射系数对比图;
图4.d双检波器傅里叶有限差分叠前深度偏移算法在Marmousi模型的A位置计算的反射系数与相同位置处理论反射系数对比图;
图4.e双检波器高阶广义屏叠前深度偏移在Marmousi模型的A位置计算的反射系数与相同位置处理论反射系数对比图;
图5.a单程波广义屏叠前深度偏移算法在Marmousi模型的B位置计算的反射系数与相同位置处理论反射系数对比图;
图5.b叠前逆时偏移算法在Marmousi模型的B位置计算的反射系数与相同位置处理论反射系数对比图;
图5.c双检波器裂步傅里叶叠前深度偏移算法在Marmousi模型的B位置计算的反射系数与相同位置处理论反射系数对比图;
图5.d双检波器傅里叶有限差分叠前深度偏移算法在Marmousi模型的B位置计算的反射系数与相同位置处理论反射系数对比图;
图5.e双检波器高阶广义屏叠前深度偏移在Marmousi模型的B位置计算的反射系数与相同位置处理论反射系数对比图;
图6.a单程波广义屏叠前深度偏移算法在Marmousi模型的C位置计算的反射系数与相同位置处理论反射系数对比图;
图6.b叠前逆时偏移算法在Marmousi模型的C位置计算的反射系数与相同位置处理论反射系数对比图;
图6.c双检波器裂步傅里叶叠前深度偏移算法在Marmousi模型的C位置计算的反射系数与相同位置处理论反射系数对比图;
图6.d双检波器傅里叶有限差分叠前深度偏移算法在Marmousi模型的C位置计算的反射系数与相同位置处理论反射系数对比图;
图6.e双检波器高阶广义屏叠前深度偏移在Marmousi模型的C位置计算的反射系数与相同位置处理论反射系数对比图。
具体实施方式
以下结合附图对本发明作进一步详细说明。
本发明所述的利用单程波算子的双程波叠前深度偏移方法,包括如下步骤:
步骤1,通过上下双检波器地震数据采集系统获取双检波器地震数据,形成检波器波场;
步骤2,从双检波器地震数据中提取地震子波信息,形成震源波场;
步骤3,读取单炮地震数据、速度模型和震源波场;
步骤4,对震源波场和检波器波场在时间域作傅里叶变换得到频率-空间域波场数据;
步骤5,对频率-空间域波场数据进行深度延拓计算,具体包括:
输入z和z+Δz的波场值和根据下式计算z+2Δz的波场值
已知在表波场及地下一定深度的波场进行偏移成像的方法中,频率波数域的波动方程及其边界条件为:
常微分方程(1.a)的通解形式为:
式(2)中,上式右边第一项对应的是上行波延拓算子,右边第二项对应的是下行波延拓算子;
由于目前地震常规采集技术只能提供一个边界条件,只能求解一个待定系数,所以在单程波波场延拓时必须要舍弃其中一项,而本发明基于上下双检波器地震数据采集系统,所以在波场延拓时上述两项的待定系数均可准确求出,体现了本发明的双程波偏移特点;
式(2)中,C1和C2为待定系数,可有边值条件(1.b)、(1.c)为:
利用z=0和z=Δz处的波场值延拓计算z=2Δz处的波场值为:
根据欧拉公式,一般化的深度延拓公式可以进一步写为:
式(5)中,为单程波算子;
结合单程波偏移的研究成果,采用不同的近似算法来计算kz,可以得到不同的双检波器叠前深度偏移算法;
步骤6,利用互相关成像原理或反射系数成像原理进行偏移成像;
步骤7,重复步骤4~6延拓计算到最大深度,保存偏移结果,计算下一炮地震数据;
步骤8,重复步骤3~7至最后一炮,输出偏移剖面。
在上述技术方案的基础上,具体计算单程波算子中的kZ时,采用裂步傅里叶方法、傅里叶有限差分方法或高阶广义屏方法。
因此,形成了三种双检波器叠前深度偏移算法:双检波器裂步傅里叶叠前深度偏移、双检波器傅里叶有限差分叠前深度偏移和双检波器高阶广义屏叠前深度偏移。
除了上述列举的三种常规的单程波算子计算方法,还可以使用其他单程波算子计算方法来计算单程波算子中的kZ,以取得同样的甚至更好的效果。
在上述技术方案的基础上,以裂步傅里叶方法[3]来近似单程波算子时,其中kZ可以表示为
其中:vr为介质参考速度;kx和ky分别表示x方向和y方向波数;ASSF表示利用裂步傅里叶方法近似的单程波算子;表示利用裂步傅里叶方法近似计算的kZ。
在上述技术方案的基础上,以傅里叶有限差分方法[4]来近似单程波算子时,其中kZ可以表示为
其中:
AFFD表示利用傅里叶有限差分方法近似的单程波算子;表示利用傅里叶有限差分方法近似计算的kZ;
和分别表示对x和y的二阶偏导数。
在上述技术方案的基础上,以高阶广义屏方法[5]来近似单程波算子时,其中kZ可以表示为
其中:
AGSP表示利用高阶广义屏方法近似的单程波算子;表示利用高阶广义屏方法近似计算的kZ。
以下为具体实施例。
实例1:
本实例的目的是检验本发明对复杂构造的成像性能。
在本实例中使用了以下算法对Marmousi模型正演计算的炮集记录进行了偏移处理,速度模型见图1:
1、单程波高阶广义屏叠前深度偏移算法(算法1);
2、叠前逆时偏移算法(算法2);
3、双检波器裂步傅里叶叠前深度偏移算法(算法3);
4、双检波器傅里叶有限差分叠前深度偏移算法(算法4);
5、双检波器高阶广义屏叠前深度偏移算法(算法5)。
文献[6]指出采用高阶广义屏算子计算的波场值与真实波场值最为接近,因此在本发明中仅以高阶广义屏延拓算子为例说明单程波偏移算法对复杂构造的成像能力。
图2.a和图2.b分别是利用算法1和算法5计算的偏移剖面。对比上述偏移剖面(见图1中箭头所指区域),虽然本发明中的算法5利用了高阶广义屏偏移算法中相同的kZ计算公式,但本发明中的算法5突破了单程波偏移的缺点,可见本发明对Marmousi模型中复杂的构造走向和层位接触关系的成像要比常规的单程波偏移手段在相应构造位置的成像要清晰、明确,凸显了双检波器双程波深度偏移的优势。
图3.a~3.d都是基于全波动方程的偏移剖面,对比发现逆时偏移剖面在模型的浅部存在比本发明中的双程波偏移算法更加严重的低频噪音,而且本发明提出的双程波偏移方法比逆时偏移方法对模型中深部的背斜构造成像更加清晰。
对比图3.a~3.d的偏移剖面,虽然在构造成像上无法区分本发明三种偏移算法(算法3、4、5)的差异,但三者在计算构造保真性能上存在明显差异,整体上算法3的真振幅能力最好,算法4的真振幅能力次之,算法5的真振幅能力再次之。
为了对比五种偏移算法在真振幅成像上的差异,在本实例中速度模型的A,B和C位置提取了偏移算法计算的反射系数和理论反射系数,对比曲线见图4.a~6.e,并统计了两者的相关系数,见表1。
表1
在表1中,从计算反射系数与理论反射系数的相关系数可明显的发现本发明提出的双检波器叠前深度偏移算法具有更准确的真振幅估计性能。
在图4.a~4.e中,黑色箭头所指位置可以发现常规的偏移算法(算法1和算法2)对深部构造具备较差的成像能力,但本发明的双检波器偏移算法(算法3、算法4和算法5)都较好的成像了深部构造。
在图5.a~5.e中和图6.a~6.e中,黑色实线为上述五种偏移算法计算的反射系数,黑色虚线为理论反射系数。直观对比,发现利用本发明的双检波器偏移算法(算法3、算法4和算法5)计算的反射系数比常规的偏移方法(算法1和算法2)计算的反射系数与理论反射系数吻合的更好,表1中利用相关系数定性分析也说明了这一点。
本说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。
参考文献(如专利/论文/标准)
[1]刘学伟,尤加春.一种真振幅偏移成像方法.中国.201510388935.6[P].
[2]Clear J.F.Imaging of earth's interior.Blackwell ScientificPbulication.Inc.1985.
[3]Stoffa P.L.1990.Split-step Fourier migration.Geophysics,55,410-421.
[4]Ristow D.Rnlll T.1994.Fourier finite-differencemigration.Geophysics,59,1882-1893.
[5]de Hoop,M.V.,Le Rousseau J.H.,and Wu R.S.,2000,Generalization ofthe phase-screen approximation for the scattering of acoustic waves:WaveMotion,31,43–70.
[6]Le Rousseau,J.,and de Hoop,M.,2001,Modeling and imaging with thescalar generalized-screen algorithms in isotropic media.Geophysics,66,1551–1568.
Claims (2)
1.一种利用单程波算子的双程波叠前深度偏移方法,其特征在于,包括如下步骤:
步骤1,通过上下双检波器地震数据采集系统获取双检波器地震数据,形成检波器波场;
步骤2,从双检波器地震数据中提取地震子波信息,形成震源波场;
步骤3,读取单炮地震数据、速度模型和震源波场;
步骤4,对震源波场和检波器波场在时间域作傅里叶变换得到频率-空间域波场数据;
步骤5,对频率-空间域波场数据进行深度延拓计算;
步骤6,利用互相关成像原理或反射系数成像原理进行偏移成像;
步骤7,重复步骤4~6延拓计算到最大深度,保存偏移结果,计算下一炮地震数据;
步骤8,重复步骤3~7至最后一炮,输出偏移剖面;
步骤5具体包括:
输入z和z+Δz的波场值和根据下式计算z+2Δz的波场值
式(5)中,为单程波算子;
采用不同的近似算法来计算kz,则得到不同的双检波器叠前深度偏移算法;
具体计算单程波算子中的kz时,采用裂步傅里叶方法、傅里叶有限差分方法或高阶广义屏方法;
分别形成:双检波器裂步傅里叶叠前深度偏移、双检波器傅里叶有限差分叠前深度偏移和双检波器高阶广义屏叠前深度偏移;
以裂步傅里叶方法来近似单程波算子时,其中kz可以表示为:
其中:vr为介质参考速度;kx和ky分别表示x方向和y方向波数;ASSF表示利用裂步傅里叶方法近似的单程波算子;表示利用裂步傅里叶方法近似计算的kz;
以高阶广义屏方法来近似单程波算子时,其中kz可以表示为:
其中:
AGSP表示利用高阶广义屏方法近似的单程波算子;表示利用高阶广义屏方法近似计算的kz;
其中,ω为角频率;v为模型速度;Δz为z方向网格剖分间距,也是延拓步长;kZ为z方向的波数。
2.如权利要求1所述的利用单程波算子的双程波叠前深度偏移方法,其特征在于:以傅里叶有限差分方法来近似单程波算子时,其中kz可以表示为:
其中:
AFFD表示利用傅里叶有限差分方法近似的单程波算子;表示利用傅里叶有限差分方法近似计算的kz;
和分别表示对x和y的二阶偏导数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610257053.0A CN105911587B (zh) | 2016-04-22 | 2016-04-22 | 一种利用单程波算子的双程波叠前深度偏移方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610257053.0A CN105911587B (zh) | 2016-04-22 | 2016-04-22 | 一种利用单程波算子的双程波叠前深度偏移方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105911587A CN105911587A (zh) | 2016-08-31 |
CN105911587B true CN105911587B (zh) | 2019-04-09 |
Family
ID=56751736
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610257053.0A Active CN105911587B (zh) | 2016-04-22 | 2016-04-22 | 一种利用单程波算子的双程波叠前深度偏移方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105911587B (zh) |
Families Citing this family (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107179551B (zh) * | 2017-06-19 | 2018-04-06 | 吉林大学 | 一种利用微震记录对地下构造直接成像的方法 |
CN108181657B (zh) * | 2017-12-28 | 2019-02-12 | 中国石油大学(北京) | 全波形反演梯度计算中分离偏移和层析成像模式的方法 |
CN108259582A (zh) * | 2017-12-29 | 2018-07-06 | 上海威惠智能科技有限公司 | 定位系统及方法 |
CN111610559B (zh) * | 2019-02-22 | 2022-05-06 | 中国石油天然气股份有限公司 | 目的层叠前深度偏移成像方法及装置 |
CN111812708B (zh) * | 2019-04-11 | 2022-08-30 | 中国石油天然气股份有限公司 | 地震波成像方法及装置 |
CN110058304B (zh) * | 2019-05-06 | 2020-01-14 | 吉林大学 | 基于单程波传播算子的多次波可控照明定量分析方法 |
CN111077566B (zh) * | 2019-12-10 | 2021-02-09 | 成都理工大学 | 一种基于矩阵分解的双程波叠前深度偏移的方法 |
CN111077567B (zh) * | 2019-12-10 | 2021-02-19 | 成都理工大学 | 一种基于矩阵乘法的双程波叠前深度偏移的方法 |
CN111781647B (zh) * | 2020-07-13 | 2022-05-20 | 中油奥博(成都)科技有限公司 | 一种大斜井炮检移动vsp自由表面多次波成像方法和装置 |
CN112782761A (zh) * | 2020-10-16 | 2021-05-11 | 中国石油大学(华东) | 单程波正演模拟方法、装置、存储介质及处理器 |
CN112946737B (zh) * | 2021-01-20 | 2023-10-31 | 中国地质大学(北京) | 一种利用纵横波速度增量交汇图识别天然气水合物的方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101839999A (zh) * | 2009-03-20 | 2010-09-22 | 中国石油集团东方地球物理勘探有限责任公司 | 一种确定叠前时间偏移最佳速度剖面的方法 |
CN103293553A (zh) * | 2013-04-17 | 2013-09-11 | 中国海洋石油总公司 | 一种复杂海底上下缆地震采集数据边界元延拓校正方法 |
CN103901473A (zh) * | 2014-04-14 | 2014-07-02 | 中国海洋石油总公司 | 一种基于非高斯性最大化的双检信号上下行波场分离方法 |
WO2015143195A1 (en) * | 2014-03-20 | 2015-09-24 | Schlumberger Canada Limited | Reconstructing impulsive source seismic data from time distributed firing airgun array data |
CN104991268A (zh) * | 2015-07-03 | 2015-10-21 | 中国地质大学(北京) | 一种真振幅偏移成像方法 |
-
2016
- 2016-04-22 CN CN201610257053.0A patent/CN105911587B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101839999A (zh) * | 2009-03-20 | 2010-09-22 | 中国石油集团东方地球物理勘探有限责任公司 | 一种确定叠前时间偏移最佳速度剖面的方法 |
CN103293553A (zh) * | 2013-04-17 | 2013-09-11 | 中国海洋石油总公司 | 一种复杂海底上下缆地震采集数据边界元延拓校正方法 |
WO2015143195A1 (en) * | 2014-03-20 | 2015-09-24 | Schlumberger Canada Limited | Reconstructing impulsive source seismic data from time distributed firing airgun array data |
CN103901473A (zh) * | 2014-04-14 | 2014-07-02 | 中国海洋石油总公司 | 一种基于非高斯性最大化的双检信号上下行波场分离方法 |
CN104991268A (zh) * | 2015-07-03 | 2015-10-21 | 中国地质大学(北京) | 一种真振幅偏移成像方法 |
Non-Patent Citations (1)
Title |
---|
罗焕宏.基于单程波算子的地震波场模拟及叠前深度偏移.《中国优秀硕士学位论文全文数据库•基础科学辑》.2011,(第04期), |
Also Published As
Publication number | Publication date |
---|---|
CN105911587A (zh) | 2016-08-31 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105911587B (zh) | 一种利用单程波算子的双程波叠前深度偏移方法 | |
Song et al. | Moho imaging based on receiver function analysis with teleseismic wavefield reconstruction: Application to South China | |
Berkhout et al. | Imaging of multiple reflections | |
Mallick et al. | Amplitude-variation-with-offset and prestack-waveform inversion: A direct comparison using a real data example from the Rock Springs Uplift, Wyoming, USA | |
US10317548B2 (en) | Reflection seismic data Q tomography | |
Wang et al. | Inversion of seismic refraction and reflection data for building long-wavelength velocity models | |
Tauzin et al. | Receiver functions from seismic interferometry: a practical guide | |
EP2548052B1 (en) | System and method of 3d salt flank vsp imaging with transmitted waves | |
Wu et al. | Microseismic source locations with deconvolution migration | |
CN113805237B (zh) | 偏移陆地交叉排列地震的使用压缩感测模型的方法和系统 | |
Pafeng et al. | Prestack waveform inversion of three-dimensional seismic data—An example from the Rock Springs Uplift, Wyoming, USA | |
US11561312B2 (en) | Mapping near-surface heterogeneities in a subterranean formation | |
Li et al. | Periodic plane-wave least-squares reverse time migration for diffractions | |
Roy Chowdhury | Seismic data acquisition and processing | |
Li et al. | Waveform inversion of seismic first arrivals acquired on irregular surface | |
da Costa Filho et al. | Attenuating multiple-related imaging artifacts using combined imaging conditions | |
US10379245B2 (en) | Method and system for efficient extrapolation of a combined source-and-receiver wavefield | |
Nguyen et al. | Inversion of Scholte wave dispersion and waveform modeling for shallow structure of the Ninetyeast Ridge | |
Oren et al. | Image-domain DAS 3D VSP elastic transmission tomography | |
Chen et al. | Research on vertical cable seismic interferometry imaging | |
Jiang et al. | Fast least-squares prestack time migration via accelerating the explicit calculation of Hessian matrix with dip-angle Fresnel zone | |
Novotny | Trace interpolation by slant-stack migration | |
Sæther | Seismic Forward Modeling of Deltaic Sequences | |
Chen et al. | Research on improved cross-correlation seismic interferometry method based on marine vertical cable seismic data | |
Cabrera-Pérez | Multiscale Ambient Noise Tomography for volcano monitoring and geothermal exploration in volcanic areas |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
CB03 | Change of inventor or designer information |
Inventor after: Liu Xuewei Inventor after: You Jiachun Inventor after: Han Wengong Inventor after: Zhang Guangde Inventor after: Yang Dekuan Inventor before: Liu Xuewei Inventor before: You Jiachun |
|
CB03 | Change of inventor or designer information | ||
GR01 | Patent grant | ||
GR01 | Patent grant |