CN113777652A - 海洋地震数据震源子波提取方法、装置和存储介质 - Google Patents

海洋地震数据震源子波提取方法、装置和存储介质 Download PDF

Info

Publication number
CN113777652A
CN113777652A CN202111083115.8A CN202111083115A CN113777652A CN 113777652 A CN113777652 A CN 113777652A CN 202111083115 A CN202111083115 A CN 202111083115A CN 113777652 A CN113777652 A CN 113777652A
Authority
CN
China
Prior art keywords
seismic
current
signal
channel
frequency domain
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
CN202111083115.8A
Other languages
English (en)
Other versions
CN113777652B (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.)
Sun Yat Sen University
Original Assignee
Sun Yat Sen University
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 Sun Yat Sen University filed Critical Sun Yat Sen University
Priority to CN202111083115.8A priority Critical patent/CN113777652B/zh
Publication of CN113777652A publication Critical patent/CN113777652A/zh
Application granted granted Critical
Publication of CN113777652B publication Critical patent/CN113777652B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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/32—Transforming one recording into another or one representation into another
    • G01V1/325—Transforming one representation into another
    • 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/20—Trace signal pre-filtering to select, remove or transform specific events or signal components, i.e. trace-in/trace-out
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/40—Transforming data representation

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

本发明公开了海洋地震数据震源子波提取方法、装置和存储介质，该方法包括：从观测海洋地震数据中截取包含鬼波的直达波窗口，以获得各地震道的叠加信号Wdg(t)，所述叠加信号Wdg(t)包含有当前地震道的直达波信号Wd(t)和震源鬼波信号Wg(t)；对任一地震道的叠加信号Wdg(t)进行频域转换，获得当前地震道的频域叠加信号
Figure DDA0003264710040000011
计算当前地震道上鬼波与直达波的到时差Δt；依据当前地震道上鬼波与直达波的到时差Δt，计算不同圆频率ω时当前地震道对应的1‑e‑iωΔt的数值；依据当前地震道的频域叠加信号
Figure DDA0003264710040000012
和不同圆频率ω时当前地震道对应的1‑e‑iωΔt的数值，计算当前地震道的频域直达波信号
Figure DDA0003264710040000013
对当前地震道的频域直达波信号
Figure DDA0003264710040000014
进行时域转换，获得当前地震道的地震震源子波Wd(t)；本发明能够从直达波和鬼波叠加信号中直接求取地震震源子波。

Description

海洋地震数据震源子波提取方法、装置和存储介质
技术领域
本发明涉及地震数据处理技术领域，尤其涉及海洋地震数据震源子波提取方法、装置和存储介质。
背景技术
在海洋地震资料处理时，子波的提取具有重要的作用。地震波激发后在地层介质中传播时，由于大地滤波作用，尖锐的激发脉冲波形变成具有一定延续时间和有限能量的地震震源子波形态。为提高地震数据的分辨率，通常针对子波的形态进行子波整形等处理，此外，子波也是数据滤波、正演模拟和波形反演等海洋地震数据处理阶段的重要参数，但这些处理都是以子波的形态已知为前提的。因此，从观测的海洋地震数据中提取出地震震源子波一直是地震数据处理中的重点和难点。
如图1所示，地震数据中的直达波，是地震震源激发后由激发源直接传播到接收电缆(拖缆)被记录的波形。因未曾在地下界面发生反射、折射，故只携带最浅层的速度信息，不携带地下地层界面的信息。即直达波的波形直观表征了地震震源信号子波的形态，因此相对于反射波等其它地震波型来说，利用直达波提取地震震源子波具有一定的优势。此外，直达波传播路径最短，在多数地震道上通常最早出现，不受反射波、折射波等其它波形的干扰，相对容易识别。因此，利用直达波提取地震震源子波是地震数据处理中的一个主要研究方向。
2019年李福元等在石油地球物理勘探期刊发表了一篇文章“从海洋地震资料直达波提取地震震源子波”。该文通过推导直达波时距曲线，考虑了震源和接收系统组合效应，结合气泡振荡理论，得到直达波与震源信号之间的关系式，并给出了频率域利用直达波计算震源信号远场子波的解析解。当地震数据接收位置满足震源组合阵列远场条件时，可利用直达波计算震源信号远场子波。
利用直达波提取地震震源子波具有一定的优势，但在海洋地震勘探中通常存在震源鬼波(海洋地震勘探中，因激发源和接收电缆均沉放在水面以下一定深度，地震波激发后，除直接传播到接收电缆上水听器处的直达波外，还有激发后先向上传播到水面反射一次后再传播到接收电缆上水听器处的波，即伴生直达波存在的“震源鬼波”)，震源鬼波与直达波在到达时间上相差不大，二者叠加在一起，使从直达波中直接提取地震震源子波具有较大的困难。
图2示出了直达波和震源鬼波传播路径示意图，朝下射线为由震源直接传播到水听器的直达波，虚线射线为根据虚震源原理得到的由虚震源激发后直接传播到水听器的震源鬼波路径(与震源激发后在水面反射一次再到达水听器的路径是等同的)。激发源和接收电缆的沉放深度相对于地震波的总传播路径较小，因此震源鬼波信号的到达时间与直达波的到达时间非常接近。但因鬼波有在水面反射的过程，而水面处地震波的反射系数通常为-1。因此，震源鬼波在地震道上的特征通常为在直达波之后较小的时间间隔以符号相反的起跳方向出现。又因为子波通常具有一定的时间延续长度，因此直达波和紧随其后的震源鬼波叠加在一起，在时间域无法区分。
发明内容
本发明的目的是提供海洋地震数据震源子波提取方法、装置和存储介质，能够从直达波和震源鬼波叠加信号中直接求取地震震源子波，以供服务于后续的子波整形、正演模拟及波形反演等处理需求。
为了实现上有目的，本发明公开了一种海洋地震数据震源子波提取方法，其包括如下步骤：
S1、从观测海洋地震数据中截取包含震源鬼波的直达波窗口，以获得各地震道的叠加信号Wdg(t)，所述叠加信号Wdg(t)包含有当前地震道的直达波信号Wd(t)和震源鬼波信号Wg(t)；
S2、对任一地震道的叠加信号Wdg(t)进行频域转换，获得当前地震道的频域叠加信号
Figure BDA0003264710020000021
S3、计算当前地震道上震源鬼波与直达波的到时差Δt；
S4、依据当前地震道上震源鬼波与直达波的到时差Δt，计算不同圆频率ω时当前地震道对应的1-e-iωΔt的数值；
S5、依据当前地震道的频域叠加信号
Figure BDA0003264710020000031
和不同圆频率ω时当前地震道对应的1-e-iωΔt的数值，计算当前地震道的频域直达波信号
Figure BDA0003264710020000032
S6、对当前地震道的频域直达波信号
Figure BDA0003264710020000033
进行时域转换，获得当前地震道的地震震源子波Wd(t)。
与现有技术相比，本发明在包含有当前地震道的直达波信号Wd(t)和震源鬼波信号Wg(t)的叠加信号Wdg(t)的基础上，通过计算当前地震道上震源鬼波与直达波的到时差Δt和不同圆频率ω时当前地震道对应的1-e-iωΔt的数值，直接从直达波和震源鬼波叠加信号中直接求取地震震源子波Wd(t)，以供服务于后续的子波整形、正演模拟及波形反演等处理需求，其地震震源子波Wd(t)的提取方式简单、快捷，适于通过算法实现。
较佳地，所述步骤S1中，从观测海洋地震数据中截取包含震源鬼波的直达波窗口，以获得各地震道的叠加信号Wdg(t)，进一步包括：
S11、以任一地震道上的直达波起跳时间为标准时间值，将各地震道上的起跳时间拉平至所述标准时间值。
具体地，以第一个地震道上的直达波起跳时间作为所述标准时间值。
较佳地，所述步骤S2具体包括：
S21、利用傅里叶变换计算任一地震道的叠加信号Wdg(t)对应的频域，获得当前地震道的频域叠加信号
Figure BDA0003264710020000034
较佳地，所述S3具体包括：
S31、从观测系统中获取震源与水面的垂向距离zs、拖缆与水面的垂向距离zg、震源与水听器的水平距离xsg，设震源至水听器的距离为dsg，震源相对于水面的虚震源至水听器的距离为ds′g，地震波在水中的传播速度为Vw，当前地震道上直达波的到达时间为td，当前地震道上震源鬼波的到达时间为tg；
S32、依据公式
Figure BDA0003264710020000041
计算当前地震道上直达波的到达时间td，及依据公式
Figure BDA0003264710020000042
计算当前地震道上震源鬼波的到达时间tg；
S33、依据公式
Figure BDA0003264710020000043
计算当前地震道上震源鬼波与直达波的到时差Δt。
较佳地，所述步骤S5具体包括：
S51、依据公式
Figure BDA0003264710020000044
计算当前地震道的频域直达波信号
Figure BDA0003264710020000045
较佳地，所述步骤S6具体包括：
S61、利用傅里叶反变换计算当前地震道的频域直达波信号
Figure BDA0003264710020000046
的时域，获得当前地震道的地震震源子波Wd(t)。
相应地，本发明还公开了一种海洋地震数据震源子波提取装置，其包括：
获取模块，被配置为从观测海洋地震数据中截取包含震源鬼波的直达波窗口，以获得各地震道的叠加信号Wdg(t)，所述叠加信号Wdg(t)包含有当前地震道的直达波信号Wd(t)和震源鬼波信号Wg(t)；
第一转换模块，被配置为对任一地震道的叠加信号Wdg(t)进行频域转换，获得当前地震道的频域叠加信号
Figure BDA0003264710020000047
第一处理模块，被配置为计算当前地震道上震源鬼波与直达波的到时差Δt；
第二处理模块，被配置为依据当前地震道上震源鬼波与直达波的到时差Δt，计算不同圆频率ω时当前地震道对应的1-e-iωΔt的数值；
第三处理模块，被配置为依据当前地震道的频域叠加信号
Figure BDA0003264710020000048
和不同圆频率ω时当前地震道对应的1-e-iωΔt的数值，计算当前地震道的频域直达波信号
Figure BDA0003264710020000051
第二转换模块，被配置为对当前地震道的频域直达波信号
Figure BDA0003264710020000052
进行时域转换，获得当前地震道的地震震源子波Wd(t)。
较佳地，所述第一转换模块，进一步被配置为以任一地震道上的直达波起跳时间为标准时间值，将各地震道上的起跳时间拉平至所述标准时间值。
相应地，本发明还公开了一种存储介质，用于存储计算机程序，所述程序被处理器执行时实现如上所述的海洋地震数据震源子波提取方法。
附图说明
图1是直达波和反射波的传播路径及波形示意图；
图2是直达波和震源鬼波传播路径示意图；
图3是本发明的海洋地震数据子波提取方法的流程框图；
图4是本发明利用数值模拟获得的脉冲形式的地震剖面；
图5是本发明各地震道求取的震源鬼波与直达波的到时差；
图6是本发明利用数值模拟获得的子波形式的地震剖面；
图7是图6拉平后的地震剖面；
图8是图7拉平后地震剖面的振幅谱；
图9是图7拉平后地震剖面的解缠后相位谱；
图10是本发明处理效果时间域信号对比图；
图11是本发明处理效果振幅谱对比图；
图12是本发明处理效果相位谱对比图；
图13是图7中各地震道利用本发明得到的子波形态图；
图14是本发明各地震道利用本发明得到的子波的振幅谱；
图15是本发明各地震道利用本发明得到的子波的相位谱；
图16是本发明的海洋地震数据子波提取装置的结构框图。
具体实施方式
为详细说明本发明的技术内容、构造特征、所实现目的及效果，以下结合实施方式并配合附图详予说明。
请参阅图3所示，本实施例的海洋地震数据震源子波提取方法，是一种基于直达波的海洋地震数据子波提取方法，其适于从直达波和震源鬼波叠加信号中直接求取地震震源子波。
这里的地震震源子波具体为地震波激发后在地层介质中传播时，由于大地滤波作用，尖锐的激发脉冲波形变成具有一定延续时间和有限能量的形态，此时大地相当于一个低通滤波器，宽频带的尖锐脉冲经过大地滤波器后，高频成分相对低频损耗更多，使得信号的频带宽度变窄，在时间域的延续时间长度变大。
请参阅图1-图3所示，本实施例的海洋地震数据震源子波提取方法包括如下步骤：
S1、从观测海洋地震数据中截取包含震源鬼波的直达波窗口，以获得各地震道的叠加信号Wdg(t)，所述叠加信号Wdg(t)包含有当前地震道的直达波信号Wd(t)和震源鬼波信号Wg(t)。
实际上，由于不同地震道上的直达波起跳时间不同，为了统一起跳时间以简化后续计算，较佳地，所述步骤S1中，从观测海洋地震数据中截取包含震源鬼波的直达波窗口，以获得各地震道的叠加信号Wdg(t)，进一步包括：
S11、以任一地震道上的直达波起跳时间为标准时间值，将各地震道上的起跳时间拉平至所述标准时间值。
具体地，以第一个地震道上的直达波起跳时间作为所述标准时间值。这里的第一个地震道具体为所有地震道中直达波起跳时间最早的那一个地震道。
将各地震道上的起跳时间拉平至所述标准时间值的目的是便于后续处理和各地震道波形的比较，如不执行上述拉平步骤可能会影响后续叠加信号Wdg(t)的相位谱。
S2、对任一地震道的叠加信号Wdg(t)进行频域转换，获得当前地震道的频域叠加信号
Figure BDA0003264710020000071
较佳地，所述步骤S2具体包括：
S21、利用傅里叶变换计算任一地震道的叠加信号Wdg(t)对应的频域，获得当前地震道的频域叠加信号
Figure BDA0003264710020000072
将时域的叠加信号Wdg(t)通过傅里叶变换为频域叠加信号
Figure BDA0003264710020000073
以便于后续在频域上进行计算，简化计算过程。
S3、计算当前地震道上震源鬼波与直达波的到时差Δt。
较佳地，所述S3具体包括：
S31、从观测系统中获取震源与水面的垂向距离zs、拖缆(或水听器)与水面的垂向距离zg、震源与水听器的水平距离xsg，上述数据均能够从海洋地震勘探数据采集的观测系统中直接获取。
本实施例能够直接在已有数据的基础上通过计算获得当前地震道的地震震源子波Wd(t)。
S32、依据公式
Figure BDA0003264710020000074
计算当前地震道上直达波的到达时间td，及依据公式
Figure BDA0003264710020000075
计算当前地震道上震源鬼波的到达时间tg。
其中dsg为震源至水听器的距离，ds′g为震源相对于水面的虚震源至水听器的距离，Vw为地震波在水中的传播速度，td为当前地震道上直达波的到达时间，tg为当前地震道上震源鬼波的到达时间。
可以理解的是，因海洋地震勘探中，直达波只携带水层速度的信息，即各地震道上直达波的到达时间为震源至水听器的距离(该距离可根据震源和水听器的相对坐标直接计算)除以地震波在水层的传播速度(通常为1500m/s的常速度)，故各地震道上直达波和震源鬼波的到达时间可通过上述公式直接计算获得。
S33、依据公式
Figure BDA0003264710020000081
计算当前地震道上震源鬼波与直达波的到时差Δt。
由图2可知，对同一地震道而言，震源鬼波的到达时间总是滞后于直达波的到达时间，因此通过将当前地震道震源鬼波的到达时间tg和直达波的到达时间td进行作差处理，即可获得当前地震道上震源鬼波与直达波的到时差Δt。
S4、依据当前地震道上震源鬼波与直达波的到时差Δt，计算不同圆频率ω时当前地震道对应的1-e-iωΔt的数值。
S5、依据当前地震道的频域叠加信号
Figure BDA0003264710020000082
和不同圆频率ω时当前地震道对应的1-e-iωΔt的数值，计算当前地震道的频域直达波信号
Figure BDA0003264710020000083
较佳地，所述步骤S5具体包括：
S51、依据公式
Figure BDA0003264710020000084
计算当前地震道的频域直达波信号
Figure BDA0003264710020000085
S6、对当前地震道的频域直达波信号
Figure BDA0003264710020000086
进行时域转换，获得当前地震道的地震震源子波Wd(t)。
较佳地，所述步骤S6具体包括：
S61、利用傅里叶反变换计算当前地震道的频域直达波信号
Figure BDA0003264710020000087
的时域，获得当前地震道的地震震源子波Wd(t)。
对当前地震道而言，因震源鬼波在水面反射时的反射系数约为-1，且时间滞后直达波Δt，故震源鬼波信号与直达波信号有如下关系：
Wg(t)＝-Wd(t-Δt)。
此时，直达波与震源鬼波的时间域叠加信号可用公式表述为：
Wdg(t)＝Wd(t)+Wg(t)＝Wd(t)-Wd(t-Δt)，
时域信号经过傅里叶变换到频域后，根据傅里叶变换的性质，公式Wdg(t)＝Wd(t)+Wg(t)＝Wd(t)-Wd(t-Δt)
傅里叶变换后得到公式
Figure BDA0003264710020000091
其中
Figure BDA0003264710020000092
对应Wdg(t)的傅里叶变换，
Figure BDA0003264710020000093
对应Wd(t)的傅里叶变换，
Figure BDA0003264710020000094
对应Wg(t)的傅里叶变换，由此可推出
Figure BDA0003264710020000095
即直达波(无震源鬼波干扰时即为地震震源子波)的频谱为直达波与震源鬼波叠加波的频谱除以1-e-iωΔt，其中ω为圆频率。利用傅里叶反变换将直达波的频谱变换到时间域即可求得地震震源子波。
请参阅图1-图15所示，为了验证本实施例的效果，下面利用数值模拟的方式对本实施例进行验证：
图4所示为本实施例利用数值模拟方式得到的脉冲形式的地震剖面，设接收电缆上共有48个水听器(每个水听器接收到的随时间变化的振动图称为一个地震道)，震源与水面的垂向距离为12米，水听器与水面的垂向距离为18米，震源与第一个水听器的水平距离为150米，水听器间的水平距离为12.5米，地震波在水中的传播速度设为1500米/秒，模拟的时间域地震剖面每个地震道共有6000个采样点，时间采样间隔为0.0001秒。图4中每个地震道有一正一负两个脉冲，分别表示直达波和鬼波到达各地震道上的时间。从图4中可见，随着水听器与震源水平距离的增大，地震道上鬼波与直达波的到时差越来越小。图5表示各地震道求取的震源鬼波与直达波的到时差，其横轴表示各地震道对应的水听器与震源的水平距离。
图6为本实施例利用数值模拟获得的子波形式的地震剖面，是图4与主频30Hz的Ricker子波褶积后的结果，即主频30Hz的Ricker子波(地震震源子波的时间延续长度是该主频所对应的周期值)就是本实施例的数值模拟要最终求取的地震震源子波形态，后续本实施例处理得到的子波形态与其一致则表明本实施例的处理效果较好。图6示出了模拟直达波被震源鬼波影响的观测数据。从图6中可见，受震源鬼波影响，无法从直达波中直接看出主频30Hz的Ricker子波形态，而且随着水听器与震源水平距离的增大，直达波和震源鬼波的到时差减小，叠加后的波形幅值越小，越偏离Ricker子波形态。
对图6中各地震道中的波形采用拉平处理，即将所有地震道上直达波的起跳时间均校正到第1个地震道上的起跳时间，拉平后地震剖面如图7所示。拉平的目的是便于后续处理和各道波形的比较，如不拉平可能会影响信号的相位谱。图8和图9为拉平后地震剖面的振幅谱和相位谱，由图7利用傅里叶变换获得。从图8中可见，受震源鬼波的影响，各地震道上信号的振幅谱幅值逐渐减小、频带变窄。从图9可见，由于采用了拉平处理，各地震道上信号的相位谱基本一致。
图10为利用实施例处理效果时间域信号对比图。图中实线为数值模拟时采用的Ricker子波波形，也是本实施例要求取的时间域子波信号形态，二者一致表明本实施例的效果较好。虚线为利用图7中第10个地震道的数据(在图10中振幅较小的那条线)计算得到的本实施例处理得到的子波。与实线对比可见，本实施例得到的子波和期望得到的子波在时间域信号形态上具有高度的一致性。
图11为利用实施例处理效果振幅谱对比图，为图10利用傅里叶变换得到频谱后进一步求取振幅谱得到。图中实线为数值模拟时采用的主频为30Hz的Ricker子波的振幅谱。虚线为利用图7中第10个地震道的数据(在图10中振幅较小的那条线)计算得到的本实施例处理得到的子波的振幅谱。虚线与实线对比可见本实施例得到的子波和期望得到的子波在频域振幅谱上具有高度的一致性。
图12为利用实施例处理效果相位谱(相位谱均为解缠后结果)对比图，为图10利用傅里叶变换得到频谱后进一步求取相位谱得到。图中实线为数值模拟时采用的Ricker子波波形的相位谱。虚线为利用图7中第10个地震道的数据(在图10中振幅较小的那条线)计算得到的实施例处理得到的子波的相位谱。实线和虚线对比可见实施例得到的子波和期望得到的子波在频域相位谱上具有高度的一致性，与绿线所示数据的相位谱有一个常数的相移。
图10-图13仅为利用图7中第10个地震道采取实施例获得的结果，为比较不同的地震道获得的结果，将图7中每一个地震道均进行同样的处理，得到的地震震源子波波形如图13所示。从图中可见，利用图7中不同的地震道进行处理，尽管每一地震道上直达波和震源鬼波的到时差不同，叠加波的波形不同，但利用实施例得到的子波形态在时间域上具有高度一致性。
图14为图13利用傅里叶变换得到频谱后进一步计算振幅谱获得。从图中可见，尽管不同地震道上直达波和震源鬼波的到时差不同，叠加波的波形不同，但不同地震道利用实施例获得的子波在振幅谱上主频、幅值和带宽等方面具有良好的一致性。
图15为图13利用傅里叶变换得到频谱后进一步计算解缠后相位谱获得。从图中可见，尽管不同地震道上直达波和震源鬼波的到时差不同，叠加波的波形不同，但不同地震道利用实施例获得的子波在相位谱上具有良好的一致性。
综上，从图10-图15可知，利用数值模拟的方式证明了实施例的有效性。
请参阅图16所示，相应地，本发明还公开了一种海洋地震数据震源子波提取装置，其包括：
获取模块10，被配置为从观测海洋地震数据中截取包含震源鬼波的直达波窗口，以获得各地震道的叠加信号Wdg(t)，所述叠加信号Wdg(t)包含有当前地震道的直达波信号Wd(t)和震源鬼波信号Wg(t)；
第一转换模块20，被配置为对任一地震道的叠加信号Wdg(t)进行频域转换，获得当前地震道的频域叠加信号
Figure BDA0003264710020000121
第一处理模块30，被配置为计算当前地震道上震源鬼波与直达波的到时差Δt；
第二处理模块40，被配置为依据当前地震道上震源鬼波与直达波的到时差Δt，计算不同圆频率ω时当前地震道对应的1-e-iωΔt的数值；
第三处理模块50，被配置为依据当前地震道的频域叠加信号
Figure BDA0003264710020000122
和不同圆频率ω时当前地震道对应的1-e-iωΔt的数值，计算当前地震道的频域直达波信号
Figure BDA0003264710020000123
第二转换模块60，被配置为对当前地震道的频域直达波信号
Figure BDA0003264710020000124
进行时域转换，获得当前地震道的地震震源子波Wd(t)。
较佳地，所述第一转换模块20进一步被配置为以任一地震道上的直达波起跳时间为标准时间值，将各地震道上的起跳时间拉平至所述标准时间值。
相应地，本发明还公开了一种存储介质，用于存储计算机程序，所述程序被处理器执行时实现如上所述的海洋地震数据震源子波提取方法。
结合图1-图16，本发明在包含有当前地震道的直达波信号Wd(t)和震源鬼波信号Wg(t)的叠加信号Wdg(t)的基础上，通过计算当前地震道上震源鬼波与直达波的到时差Δt和不同圆频率ω时当前地震道对应的1-e-iωΔt的数值，直接从直达波和震源鬼波叠加信号中直接求取地震震源子波Wd(t)，以供服务于后续的子波整形、正演模拟及波形反演等处理需求，其地震震源子波Wd(t)的提取方式简单、快捷，适于通过算法实现。
以上所揭露的仅为本发明的优选实施例而已，当然不能以此来限定本发明之权利范围，因此依本发明申请专利范围所作的等同变化，仍属本发明所涵盖的范围。

Claims (10)

1.一种海洋地震数据震源子波提取方法，其特征在于，包括如下步骤：
从观测海洋地震数据中截取包含鬼波的直达波窗口，以获得各地震道的叠加信号Wdg(t)，所述叠加信号Wdg(t)包含有当前地震道的直达波信号Wd(t)和震源鬼波信号Wg(t)；
对任一地震道的叠加信号Wdg(t)进行频域转换，获得当前地震道的频域叠加信号
Figure FDA0003264710010000011
计算当前地震道上鬼波与直达波的到时差Δt；
依据当前地震道上鬼波与直达波的到时差Δt，计算不同圆频率ω时当前地震道对应的1-e-iωΔt的数值；
依据当前地震道的频域叠加信号
Figure FDA0003264710010000012
和不同圆频率ω时当前地震道对应的1-e-iωΔt的数值，计算当前地震道的频域直达波信号
Figure FDA0003264710010000013
对当前地震道的频域直达波信号
Figure FDA0003264710010000014
进行时域转换，获得当前地震道的地震震源子波Wd(t)。
2.如权利要求1所述的海洋地震数据震源子波提取方法，其特征在于，所述从观测海洋地震数据中截取包含鬼波的直达波窗口，以获得各地震道的叠加信号Wdg(t)，进一步包括：
以任一地震道上的直达波起跳时间为标准时间值，将各地震道上的起跳时间拉平至所述标准时间值。
3.如权利要求2所述的海洋地震数据震源子波提取方法，其特征在于，以第一个地震道上的直达波起跳时间作为所述标准时间值。
4.如权利要求1所述的海洋地震数据震源子波提取方法，其特征在于，所述对任一地震道的叠加信号Wdg(t)进行频域转换，获得当前地震道的频域叠加信号
Figure FDA0003264710010000021
具体包括：
利用傅里叶变换计算任一地震道的叠加信号Wdg(t)对应的频域，获得当前地震道的频域叠加信号
Figure FDA0003264710010000022
5.如权利要求1所述的海洋地震数据震源子波提取方法，其特征在于，所述计算当前地震道上鬼波与直达波的到时差Δt，具体包括：
从观测系统中获取震源与水面的垂向距离zs、拖缆与水面的垂向距离zg、震源与水听器的水平距离xsg，设震源至水听器的距离为dsg，震源相对于水面的虚震源至水听器的距离为ds′g，地震波在水中的传播速度为Vw，当前地震道上直达波的到达时间为td，当前地震道上鬼波的到达时间为tg；
依据公式
Figure FDA0003264710010000023
计算当前地震道上直达波的到达时间td，及依据公式
Figure FDA0003264710010000024
计算当前地震道上鬼波的到达时间tg；
依据公式
Figure FDA0003264710010000025
计算当前地震道上鬼波与直达波的到时差Δt。
6.如权利要求1所述的海洋地震数据震源子波提取方法，其特征在于，所述依据当前地震道的频域叠加信号
Figure FDA0003264710010000026
和不同圆频率ω时当前地震道对应的1-e-iωΔt的数值，计算当前地震道的频域直达波信号
Figure FDA0003264710010000027
具体包括：
依据公式
Figure FDA0003264710010000028
计算当前地震道的频域直达波信号
Figure FDA0003264710010000029
7.如权利要求1所述的海洋地震数据震源子波提取方法，其特征在于，所述对当前地震道的频域直达波信号
Figure FDA0003264710010000031
进行时域转换，获得当前地震道的地震震源子波Wd(t)，具体包括：
利用傅里叶反变换计算当前地震道的频域直达波信号
Figure FDA0003264710010000032
的时域，获得当前地震道的地震震源子波Wd(t)。
8.一种海洋地震数据震源子波提取装置，其特征在于，包括：
获取模块，被配置为从观测海洋地震数据中截取包含鬼波的直达波窗口，以获得各地震道的叠加信号Wdg(t)，所述叠加信号Wdg(t)包含有当前地震道的直达波信号Wd(t)和震源鬼波信号Wg(t)；
第一转换模块，被配置为对任一地震道的叠加信号Wdg(t)进行频域转换，获得当前地震道的频域叠加信号
Figure FDA0003264710010000033
第一处理模块，被配置为计算当前地震道上鬼波与直达波的到时差Δt；
第二处理模块，被配置为依据当前地震道上鬼波与直达波的到时差Δt，计算不同圆频率ω时当前地震道对应的1-e-iωΔt的数值；
第三处理模块，被配置为依据当前地震道的频域叠加信号
Figure FDA0003264710010000034
和不同圆频率ω时当前地震道对应的1-e-iωΔt的数值，计算当前地震道的频域直达波信号
Figure FDA0003264710010000035
第二转换模块，被配置为对当前地震道的频域直达波信号
Figure FDA0003264710010000036
进行时域转换，获得当前地震道的地震震源子波Wd(t)。
9.如权利要求8所述的海洋地震数据震源子波提取装置，其特征在于，所述第一转换模块，进一步被配置为以任一地震道上的直达波起跳时间为标准时间值，将各地震道上的起跳时间拉平至所述标准时间值。
10.一种存储介质，用于存储计算机程序，其特征在于：所述程序被处理器执行时实现如权利要求1～7中任一项所述的海洋地震数据震源子波提取方法。
CN202111083115.8A 2021-09-15 2021-09-15 海洋地震数据震源子波提取方法、装置和存储介质 Active CN113777652B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111083115.8A CN113777652B (zh) 2021-09-15 2021-09-15 海洋地震数据震源子波提取方法、装置和存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111083115.8A CN113777652B (zh) 2021-09-15 2021-09-15 海洋地震数据震源子波提取方法、装置和存储介质

Publications (2)

Publication Number Publication Date
CN113777652A true CN113777652A (zh) 2021-12-10
CN113777652B CN113777652B (zh) 2022-05-27

Family

ID=78844303

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111083115.8A Active CN113777652B (zh) 2021-09-15 2021-09-15 海洋地震数据震源子波提取方法、装置和存储介质

Country Status (1)

Country Link
CN (1) CN113777652B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120092961A1 (en) * 2010-10-13 2012-04-19 The Petroleum Institute Analysis and filtering of surface waves in shallow water environment using s and t-f-k transform
CN104536045A (zh) * 2015-01-16 2015-04-22 中国海洋石油总公司 一种基于子波处理的鬼波压制方法
US20160109601A1 (en) * 2014-10-20 2016-04-21 Pgs Geophysical As Method of acquiring ghost-free signatures for broadband source calibration
CN108594302A (zh) * 2018-07-26 2018-09-28 广州海洋地质调查局 一种地震子波的提取方法及处理终端
CN110850474A (zh) * 2018-08-20 2020-02-28 中国石油化工股份有限公司 一种海洋地震资料震源鬼波压制方法和系统
CN112817047A (zh) * 2020-12-31 2021-05-18 北京东方联创地球物理技术有限公司 海洋地震自适应去鬼波方法、装置、电子设备及介质

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120092961A1 (en) * 2010-10-13 2012-04-19 The Petroleum Institute Analysis and filtering of surface waves in shallow water environment using s and t-f-k transform
US20160109601A1 (en) * 2014-10-20 2016-04-21 Pgs Geophysical As Method of acquiring ghost-free signatures for broadband source calibration
CN104536045A (zh) * 2015-01-16 2015-04-22 中国海洋石油总公司 一种基于子波处理的鬼波压制方法
CN108594302A (zh) * 2018-07-26 2018-09-28 广州海洋地质调查局 一种地震子波的提取方法及处理终端
CN110850474A (zh) * 2018-08-20 2020-02-28 中国石油化工股份有限公司 一种海洋地震资料震源鬼波压制方法和系统
CN112817047A (zh) * 2020-12-31 2021-05-18 北京东方联创地球物理技术有限公司 海洋地震自适应去鬼波方法、装置、电子设备及介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
成谷 等: "Morlet小波匹配追踪方法的时频原子参数分析", 《地球物理学进展》 *
范兴利 等: "基于Morlet小波尺度参数寻优的匹配追踪时频分析", 《中山大学学报（自然科学版）》 *

Also Published As

Publication number Publication date
CN113777652B (zh) 2022-05-27

Similar Documents

Publication Publication Date Title
US10802167B2 (en) Seismic acquisition method and apparatus
CN107894613A (zh) 弹性波矢量成像方法、装置、存储介质及设备
CN102879817A (zh) 基于地面地震数据获取地下裂缝信息的控制方法
SG194307A1 (en) Seismic data processing including compensating for source and receiver ghost effects in reverse time migration
Marsset et al. Deep-towed high resolution seismic imaging II: Determination of P-wave velocity distribution
CN113777652B (zh) 海洋地震数据震源子波提取方法、装置和存储介质
Zhang et al. Deghosting towed streamer data in τ/p domain based on rough sea surface reflectivity
CN113514889B (zh) 一种提升海洋深反射地震数据中低频信号能量的处理方法
Jeng et al. A nonlinear method of removing harmonic noise in geophysical data
WO2017203482A1 (en) System and method for acquiring and inverting sparse-frequency data
Zhang et al. One-way wave propagation in the ray-centred coordinate system for vertical transversely isotropic media
Rabinovich et al. Modeling of a reservoir fracture zone formed by hydraulic fracturing
Brahmani et al. Estimation of power spectral density of seismic data using welch method
CN116125535B (zh) 三维vsp成像的方法及装置
Tan et al. Combined adaptive multiple subtraction based on event tracing and Wiener filtering
CN117555022A (zh) 一种地震子波的确定方法、装置、电子设备以及存储介质
Yablokov et al. AUTOMATIZATION OF SPECTRAL ANALYSIS OF SURFACE WAVE DATA
CN115774285A (zh) 一种基于衰减合成记录的提高地震剖面分辨率方法和系统
Sandmeier Non-destructive testing of concrete with electromagnetic and acoustic–elastic waves: data analysis
Tao et al. Extracting near borehole P and S reflections from acoustic logging measurements

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