CN110456417B - 一种地震数据多次波压制方法 - Google Patents

一种地震数据多次波压制方法 Download PDF

Info

Publication number
CN110456417B
CN110456417B CN201910784112.3A CN201910784112A CN110456417B CN 110456417 B CN110456417 B CN 110456417B CN 201910784112 A CN201910784112 A CN 201910784112A CN 110456417 B CN110456417 B CN 110456417B
Authority
CN
China
Prior art keywords
dimensional
data
seismic data
formula
frequency
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
CN201910784112.3A
Other languages
English (en)
Other versions
CN110456417A (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 National Offshore Oil Corp CNOOC
CNOOC China Ltd Zhanjiang Branch
Original Assignee
China National Offshore Oil Corp CNOOC
CNOOC China Ltd Zhanjiang Branch
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 National Offshore Oil Corp CNOOC, CNOOC China Ltd Zhanjiang Branch filed Critical China National Offshore Oil Corp CNOOC
Priority to CN201910784112.3A priority Critical patent/CN110456417B/zh
Publication of CN110456417A publication Critical patent/CN110456417A/zh
Application granted granted Critical
Publication of CN110456417B publication Critical patent/CN110456417B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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/36Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy

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,输入三维数据d(offx,t),依据三维空间坐标信息,道集抽取为三维面元道集dbin(y,x,t);S2,计算非线性趋势拟合系数p(x),p(y),写出变换算子Lx,Ly,Lxy,计算Radon域数据M;S3,Radon域数据切除处理,得到多次波模型mra,并计算变换得到时空域多次波模型
Figure DDA0002177471070000011
S4,相减得到多次波压制后的三维面元道集p;S5,将道集按标量偏移距大小依次排列,得到多次波压制后三维数据,通过低频结果约束高频的计算,低频结果在运算中不会产生假频,在采样稀疏的情况下仍能取得比较好的高分辨率效果;同时避免了常规迭代方法提高分辨率的过程,大大提高了计算效率。

Description

一种地震数据多次波压制方法
技术领域
本发明涉及勘探地震信号处理领域,具体涉及一种地震数据多次波压制方法。
背景技术
海洋地震数据采集时,由于受到海水表面强反射系数的影响,海洋地震数据中会存在强烈的多次波;多次波的存在严重影响着地震数据速度分析和成像;如何有效去除多次波是地震资料处理中的一个重要研究内容。
Radon变换是工业界多次波压制最为常用的有效手段之一,同时也被广泛用于地震数据重建,波场分离,速度分析等。Radon变换的二维算法没有考虑地震波传播的三维效应,常规三维Radon变换算法没有考虑地层性质的各向异性,分辨率也较低,而且常规Radon变换只是沿着特定路径的叠加求和,并没有考虑横向的振幅变化,因此也存在不保幅的问题。
发明内容
本发明的目的是在于提供一种地震数据多次波压制方法,利用非线性趋势预测理论,对地震数据横向振幅变化进行拟合并存储;利用计算时低频的计算结果对高频的计算进行约束提高分辨率,实现地震数据三维高分辨率保幅Radon变换,得到更精确的多次波模型,进而达到更好的多次波压制效果;利用椭圆模型应对地层参数各向异性所带来的地震数据的方位效应,能够更好的进行复杂数据多次波的压制。
为解决上述发明的目的,本发明提供技术方案如下:
一种地震数据多次波压制方法,其特征在于,包括以下步骤:
S1,输入三维数据d(offx,t),依据三维空间坐标信息,道集抽取为三维面元道集dbin(y,x,t):输入三维数据d(offx,t),对于给定的以标量偏移距offx排列的三维共中心点道集地震数据d(offx,t),依据该道集每个标量偏移距offx所对应的纵测线和非纵测线空间坐标,将其转变抽取为三维面元道集dbin(y,x,t);其中,y为地震数据非纵测线方向的偏移距,x为地震数据纵测线方向的偏移距,t为纵坐标时间;
S2,计算非线性趋势拟合系数p(x),p(y),写出变换算子Lx,Ly,Lxy,计算Radon域数据M;计算Radon变换算子Lx,Ly,Lxy,同时将地震数据dbin(y,x,t)变换到频率域,并将如上结果代入公式中,求取得到频率域Radon域数据M;并将其做反变换至时间域,得到时间域的Radon域地震数据dra,依据一次波和多次波在Radon域的空间位置关系,对其进行切除操作,得到Radon域多次波模型mra
S3,Radon域数据切除处理,得到多次波模型mra,并依据公式变换得到时空域多次波模型
Figure BDA0002177471050000024
依据一次波和多次波在Radon域的空间位置关系,对其进行切除操作,得到Radon域多次波模型mra;将其变换到频率域Mra,结合Radon变换算子,进行反变换,得到频率域多次波模型Mra-freq,将Mra-freq变换回时间域,得到时空域的多次波模型
Figure BDA0002177471050000021
S4,相减得到多次波压制后的三维面元道集p;将得到的多次波模型
Figure BDA0002177471050000022
和三维面元道集dbin做减法,得到多次波压制后的三维面元道集p:
Figure BDA0002177471050000023
S5,将道集按标量偏移距大小依次排列,得到多次波压制后三维数据;将p按照纵测线和非纵测线偏移距大小,计算标量偏移距大小,并按照标量偏移距大小一次排列地震道,得到标量偏移距offx排列的、多次波压制后的三维共中心点道集地震数据。
进一步地,所述步骤S1中,三维面元道集dbin(y,x,t)基于椭圆模型,其表达为:
d(y,x,t)=∫∫∫m(qxy,qy,qx,τ=t-qxx2-qyy2-qxyxy,)dqxdqydqxy (1)
其中,qxy为考虑了纵测线和非纵测线方向椭圆效应的综合曲率参数;
qy为y方向曲率参数;qx为x方向曲率参数;τ为Radon域的时间;
对公式(1)离散并变换到频率域可得到:
Figure BDA0002177471050000031
可将上式写成矩阵算子相乘的形式:
D=LxMLyLxy (3)
进一步地,所述步骤S1中,对应于公式(1),考虑到同相轴振幅变化、离散化的反变换公式为:
Figure BDA0002177471050000032
进一步地,步骤S2中,计算非线性趋势拟合系数p(x),p(y),其公式为:
Figure BDA0002177471050000033
所述计算Radon域数据M,其公式为:
Figure BDA0002177471050000034
式中,Mn,Dn分别是第n个频率的多次波数据和地震数据;
Lx,Lxy,Ly分别为Radon变换计算时,包含了非线性趋势拟合系数的纵测线和非纵测线算子,可用如下公式表示:
Figure BDA0002177471050000035
进一步地,所述步骤S5中,其判断方法为:
对公式(6)中的Wn,Un,Vn由上一个频率的计算结果得出:
Wiin)=||Min-1)|| (8)
其中,Mi为在上一频率运算得到计算结果,W若为纵测线方向的加权矩阵,则i的范围为1~Nqx;
若为纵测线和非纵测线两个方向的加权矩阵,则i的范围为1~N qxy,若为横测线方向的加权矩阵V,则i的范围为1~Nqy。
本发明相对于现有技术的有益效果是:
1、三维方法考虑了地震波场的三维传播效应,能够更加精确的估计多次波的空间传播旅行时;
2、椭圆模型能够更好的应对地层性质存在各向异性时,采集得到的存在方位效应的地震数据,基于椭圆模型的Radon变换算子,在处理存在方位效应的地震数据能实现更好的Radon域能量集中;
3、非线性趋势预测能够较好的保存地震数据空间的振幅变化,从而能够更加精确的估计多次波的振幅;
4、利用低频结果约束高频的计算,更看重低频结果对整体结果的约束,低频结果在运算中不会产生假频,在采样稀疏的情况下仍能取得比较好的高分辨率效果;
5、利用低频计算结果对高频进行约束,避免了常规迭代方法提高分辨率的过程,可以大大提高计算效率。
附图说明
图1是本发明实施例中地震数据多次波压制方法的流程图;
图2是本发明实施例中标量偏移距为横坐标的三维道集地震数据图;
图3是本发明实施例中以非纵测线偏移距大小依次排列的三维小面元道集地震数据图;
图4是传统最小二乘三维抛物线Radon变换压制的多次波面元道集图;
图5是本发明实施例中最小二乘三维抛物线Radon变换多次波压制结果面元道集图;
图6是本发明实施例中标量偏移距为横坐标的传统最小二乘三维抛物线Radon变换多次波压制结果图;
图7是本发明实施例中压制的多次波面元道集图;
图8是本发明实施例中多次波压制结果小面元道集图;
图9是本发明实施例中以标量偏移距为横坐标的多次波压制结果图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚,下面将结合本发明具体实施例及相应的附图对本发明技术方案进行清楚、完整地描述。显然,所描述的实施例仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
结合图1所示,在均匀介质中,地震波前面是一个球面,有效波的波前面遇到地下反射面发生反射的时候,反射波前面也是标准的球面,当地震反射波被接收之后,对于相同的反射时间,接收到地震波的检波器表现为一个圆形的组成,因此将三维Radon变换中地震波旅行时和非纵测线、纵测线方向偏移距的关系可以表示为t=τ+qxx2+qyy2。但是当地下介质不满足密度均匀的情况下,这个模型也会变现出一定的误差。对于三维地震数据,如果地下介质在不同方向上密度变化较为明显,地震波旅行时相同时,三维的检波器分布不再是标准的圆,而是在圆的基础上发生了平滑的调整。当地震数据存在方位效应时,即不同方向上地层参数变现为各向异性,我们需要用新的模型来适用该种地层接收到的地震数据,当三维地震数据存在方位效应时,地震波的形态不再是标准的圆,而表现为椭圆t=τ+qxx2+qyy2+qxyxy,可以通过椭圆模型得到基于椭圆模型的Radon变换算子,在处理存在方位效应的地震数据实现更好的Radon域能量集中。
基于椭圆模型,可以将三维Radon变换表示为:
d(y,x,t)=∫∫∫m(qxy,qy,qx,τ=t-qxx2-qyy2-qxyxy,)dqxdqydqxy (1)
其中,qxy为考虑了纵测线和非纵测线方向椭圆效应的综合曲率参数,
qy为y方向曲率参数,qx为x方向曲率参数,τ为Radon域的时间。
相比于常规三维Radon变换,引入了qxy来拟合地震波的真实传播路径,这对存在方向各向异性的地震资料收敛地震波能量是有帮助的,将公式(1)离散并变换到频率域可得到:
Figure BDA0002177471050000061
可将上式写成矩阵算子相乘的形式:
D=LxMLyLxy (3)
在进行Radon变换求和时,通过在求和过程中加入拟合系数p(x)、p(y),可以考虑地震数据沿x和y方向空间振幅变化,从而估计出更加准确的多次波模型,得到更好的多次波压制效果。因此,对三维共中心点道集地震数据dbin(y,x,t),可进行保幅三维Radon变换,得到Radon域地震数据m(j,i,qxy,qy,qx,τ),其中,j,i分别为非纵测线和纵测线方向的非线性趋势拟合阶数,对应于公式(1),考虑了同相轴振幅变化了的、离散化的反变换公式为:
Figure BDA0002177471050000062
非线性趋势拟合系数p(x),p(y)可以依据下式进行计算:
Figure BDA0002177471050000063
进一步地,为提高m的计算精度,可在频率域依据如下公式计算高分辨率频率域Radon域地震数据M:
Figure BDA0002177471050000064
其中,Mn,Dn分别是第n个频率的多次波数据和地震数据,Lx,Lxy,Ly分别为Radon变换计算时,包含了非线性趋势拟合系数的纵测线和非纵测线算子,可用如下公式表示:
Figure BDA0002177471050000071
公式(6)中的Wn,Un,Vn由上一个频率的计算结果得到,可由下式得出:
Wiin)=||Min-1)|| (8)
其中,Mi为在上一频率运算得到计算结果,W若为纵测线方向的加权矩阵,则i的范围为1~Nqx,若为纵测线和非纵测线两个方向的加权矩阵,则i的范围为1~Nqxy,若为横测线方向的加权矩阵V,则i的范围为1~Nqy
利用低频运算结果约束高频计算的方法,避免了常规迭代方法提高分辨率的过程,可以大大提高计算效率。
结合图2-图9所示,图2为理论数据是标量偏移距为横坐标的三维道集地震数据,该道集中,一次波和多次波受到方位各向异性的影响,同相轴有一定的离散。
图3为以非纵测线偏移距大小依次排列的三维小面元道集地震数据。该道集中,沿着横测线方向较平的是一次波对应的同相轴,一次波同相轴沿y方向有较小幅度的弯曲,是由横测线间距引起的;沿y方向弯曲,且同相轴本身具有一定曲率的是多次波所对应的同相轴。该道集中,对同相轴沿纵测线方向的振幅加以变化,以模拟实际数据沿偏移距方向的振幅变化。
图4为传统最小二乘三维抛物线Radon变换压制的多次波面元道集。由于最小二乘抛物线Radon变换方法分辨率较低,不能较好的将一次波和多次波分离开,导致变换回的多次波中有一次波的信息,从而使得最终的多次波压制结果不保幅;图中可以看到有明显的一次波能量。
图5为最小二乘三维抛物线Radon变换多次波压制结果面元道集。该道集可以看出,由于传统最小二乘三维抛物线Radon变换方法未考虑同相轴振幅沿空间方向的变化,使得多次波模型振幅和原数据多次波振幅不匹配,最终导致多次波能量有残余。
图6为标量偏移距为横坐标的传统最小二乘三维抛物线Radon变换多次波压制结果。该道集中可以看到有明显的多次波能量残余。
图7为本发明方法压制的多次波面元道集。可以看出,多次波中一次波能量与图4相比明显减少,显示了本文方法分辨率高的优势。
图8为本发明方法多次波压制结果小面元道集。可以看出,由于本方法利用非线性趋势拟合,充分考虑了振幅随偏移距的变化情况,使得多次波压制后的结果中,多次波残余能量明显降低。
图9为以标量偏移距为横坐标的本发明方法多次波压制结果。可以看出,多次波压制彻底,同时一次波能量得到了很好的保护。
以上所述仅为本发明的实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的权利要求范围之内。

Claims (3)

1.一种地震数据多次波压制方法,其特征在于,包括以下步骤:
S1,输入三维数据d(offx,t),依据三维空间坐标信息,道集抽取为三维面元道集dbin(y,x,t);
S2,计算非线性趋势拟合系数p(x),p(y),写出变换算子Lx,Ly,Lxy,计算Radon域数据M;
S3,Radon域数据切除处理,得到多次波模型mra,并依据公式变换得到时空域多次波模型
Figure FDA0003038824360000011
S4,相减得到多次波压制后的三维面元道集p;
Figure FDA0003038824360000012
S5,将道集按标量偏移距大小依次排列,得到多次波压制后三维数据;
所述步骤S1中,三维面元道集dbin(y,x,t)基于椭圆模型,其表达为
din(y,x,t)=∫∫∫m(qxy,qy,qx,τ=t-qxx2-qyy2-qxyxy,)dqxdqydqxy (1)
其中,qxy为考虑了纵测线和非纵测线方向椭圆效应的综合曲率参数;
qy为y方向曲率参数;qx为x方向曲率参数;τ为Radon域的时间;
对公式(1)离散并变换到频率域可得到:
Figure FDA0003038824360000013
可将上式写成矩阵算子相乘的形式:
D=LxMLyLxy (3)
所述步骤S2中,计算非线性趋势拟合系数p(x),p(y),其公式为:
Figure FDA0003038824360000014
所述计算Radon域数据M,其公式为:
Figure FDA0003038824360000015
式中,Mn,Dn分别是第n个频率的多次波数据和地震数据;
Lx,Lxy,Ly可用如下公式表示:
Figure FDA0003038824360000021
2.如权利要求1所述的一种地震数据多次波压制方法,其特征在于,所述步骤S1中,对应于公式(1),考虑到同相轴振幅变化、离散化的反变换公式为:
Figure FDA0003038824360000022
3.如权利要求1所述的一种地震数据多次波压制方法,其特征在于,所述步骤S5中,其判断方法为:
对公式(5)中的Wn,Un,Vn由上一个频率的计算结果得出:
Wkkn)=||Mkn-1)|| (8)
其中,Mk为在上一频率运算得到计算结果,Wkk若为纵测线方向的加权矩阵,则k的范围为1~Nqx;若为纵测线和非纵测线两个方向的加权矩阵,则k的范围为1~Nqxy,若为横测线方向的加权矩阵,则k的范围为1~Nqy
CN201910784112.3A 2019-08-23 2019-08-23 一种地震数据多次波压制方法 Active CN110456417B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910784112.3A CN110456417B (zh) 2019-08-23 2019-08-23 一种地震数据多次波压制方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910784112.3A CN110456417B (zh) 2019-08-23 2019-08-23 一种地震数据多次波压制方法

Publications (2)

Publication Number Publication Date
CN110456417A CN110456417A (zh) 2019-11-15
CN110456417B true CN110456417B (zh) 2021-06-04

Family

ID=68488805

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910784112.3A Active CN110456417B (zh) 2019-08-23 2019-08-23 一种地震数据多次波压制方法

Country Status (1)

Country Link
CN (1) CN110456417B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111190222B (zh) * 2020-01-07 2021-06-25 中国海洋石油集团有限公司 一种基于滤波器形态检测的多次波自适应相减算法
CN111190226B (zh) * 2020-02-18 2022-03-25 中国石油大学(华东) 一种三维地震数据面波噪声压制方法
CN111239828B (zh) * 2020-03-09 2021-07-30 吉林大学 基于最优双曲线积分路径叠加的多次波压制方法
CN111239827B (zh) * 2020-03-09 2021-07-30 吉林大学 基于局部相似系数的三维地震数据多次波压制方法
CN111239830B (zh) * 2020-03-09 2021-07-16 吉林大学 基于局部相关加权函数的海洋地震数据自动速度分析方法
CN112327362B (zh) * 2020-10-30 2021-11-23 中国海洋大学 速度域的海底多次波预测与追踪衰减方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6636809B1 (en) * 1999-05-03 2003-10-21 Compagnie Generale De Geophysique High resolution radon transform seismic traces processing
US6574567B2 (en) * 2001-01-23 2003-06-03 Pgs Americas, Inc. Weighted slant stack for attenuating seismic noise
CN109188516B (zh) * 2018-10-31 2021-07-20 中国石油化工股份有限公司 Radon域能量扫描叠加的微地震事件定位方法
CN109633752B (zh) * 2019-01-04 2020-04-07 吉林大学 基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法

Also Published As

Publication number Publication date
CN110456417A (zh) 2019-11-15

Similar Documents

Publication Publication Date Title
CN110456417B (zh) 一种地震数据多次波压制方法
Lin et al. Eikonal tomography: surface wave tomography by phase front tracking across a regional broad-band seismic array
Wu et al. Directional illumination analysis using beamlet decomposition and propagation
Theune et al. Analysis of prior models for a blocky inversion of seismic AVA data
US8103453B2 (en) Method of seismic data interpolation by projection on convex sets
Downton et al. Linearized amplitude variation with offset (AVO) inversion with supercritical angles
Delph et al. Constraining crustal properties using receiver functions and the autocorrelation of earthquake‐generated body waves
CN107894612B (zh) 一种q吸收衰减补偿的声波阻抗反演方法及系统
EP2360495B1 (en) DIP-based corrections for data reconstruction in three-dimensional surface-related multiple prediction
EP3749987B1 (en) Systems and methods to enhance 3-d prestack seismic data based on non-linear beamforming in the cross-spread domain
Eilon et al. An adaptive Bayesian inversion for upper-mantle structure using surface waves and scattered body waves
US9952341B2 (en) Systems and methods for aligning a monitor seismic survey with a baseline seismic survey
EA032186B1 (ru) Сейсмическая адаптивная фокусировка
US5515335A (en) Seismic trace overburden correction method
CN109375265A (zh) 一种基于变相位雷克子波匹配追踪的理想地震谱分解方法
Huang et al. Pre‐stack seismic inversion based on ℓ1− 2‐norm regularized logarithmic absolute misfit function
Jensen et al. Point‐spread function convolution to simulate prestack depth migrated images: A validation study
Chen et al. Joint data and model-driven simultaneous inversion of velocity and density
Guigné et al. Acoustic zoom high-resolution seismic beamforming for imaging specular and non-specular energy of deep oil and gas bearing geological formations
CN102478663B (zh) 一种三维地震观测系统偏移噪声获取方法及装置
Abedi et al. Three‐parameter Radon transform in layered transversely isotropic media
US20230184972A1 (en) Method and system for seismic processing using virtual trace bins based on offset attributes and azimuthal attributes
WO2022256666A1 (en) Method and system for reflection-based travel time inversion using segment dynamic image warping
Zhang et al. One-way wave propagation in the ray-centred coordinate system for vertical transversely isotropic media
Wang et al. Rock Fracture Monitoring Based on High‐Precision Microseismic Event Location Using 3D Multiscale Waveform Inversion

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