CN108594302A - 一种地震子波的提取方法及处理终端 - Google Patents

Abstract

本发明涉及一种地震子波的提取方法及处理终端，所述提取方法包括依次进行的以下步骤：步骤S1：获取地震道中直达波s_j(t)，并经过傅里叶变换获得s_j(ω)；步骤S2：计算得出反射系数R谱；步骤S3：获取检波器到震源和虚震源的距离，计算出中间变量B；步骤S4：将步骤S3中的中间变量B、步骤S1中的s_j(ω)和步骤S2中的反射系数R代入公式计算出地震子波w(ω)。本发明不依赖于地震褶积理论模型，不需要测井资料，不需要对地震数据和地层反射系数的分布做任何假设，也不需要对地震子波的相位特性做任何假设；同时，本发明适用于海上拖缆地震资料的地震子波提取，计算结果精度高，可以支持超长序列地震子波的提取，满足地震资料批量处理的要求。

Description

一种地震子波的提取方法及处理终端

技术领域

本发明涉及地震数据处理技术领域，具体是一种地震子波的提取方法及处理终端。

背景技术

地震子波是地震资料数据处理的重要参数。在偏移成像处理、正演模拟、叠前/叠后属性反演、波形反演、多次波压制、鬼波压制、宽频处理等处理中起到关键作用，更为一般地，地震子波是地震勘探系统的激励函数，也是地震波弹性波动方程中的源函数，因此获得地震子波是研究地震资料数据处理非常重要的一环。

目前，地震子波的获取主要有以下四种类型方法：

1)直接记录法，即在野外一定条件限制下，直接观测记录震源子波；

2)室内模拟法，利用震源建模软件理论模拟地震子波；

3)近场拟合法，即利用震源附近的近场检波器记录的近场子波按照一定的数学方法拟合远场震源子波；

4)地震资料提取法，即利用地震数据本身提取地震子波，又分为确定性和统计性两类方法。

但以上方法均存在各自的补正，地震数据采集受调查环境、装备的影响，野外直接观测震源子波十分困难；利用近场子波记录拟合时，由于近场记录中经常会有干扰、畸变，拟合结果往往存在较大误差，在很多情况下地震数据中并没有记录近场子波，该方法的使用也受到限制；室内软件模拟方法受软件所使用理论模型、震源类型以及实际环境参数等因素的限制，模拟子波和实际子波存在差异，如果再考虑到检波器特性、地震仪器因素的影响，从地震资料中直接提取地震子波更为合理、也更实用。

随着数字信号处理理论、非线性理论、优化理论的发展，从地震数据中直接提取地震子波的方法从理论到实践都取得巨大进展，此类方法目前研究的最多、应用也最广泛。地震资料提取法依据的理论是地震记录的褶积理论模型，即地震记录是地震子波和地层反射系数的褶积，如下公式所示：

s(t)＝r(t)*w(t)

式中，s(t)表示地震记录，r(t)表示地层反射系数，w(t)表示地震子波，

根据上面的公式，在获得地层反射系数以及地震记录后，地震记录也即地震资料，即可计算得到地震子波；地震资料提取法分为确定性和统计性两类方法，其中，确定性子波提取方法不需要对地层反射系数序列的分布作任何假设,利用测井资料计算反射系数，进而推算地震子波，但对测井数据的精度有较高的要求，这一条件在很多情况下是不满足的，在构造起伏变化大的地区，也局限于提取井旁地震资料的子波，在测井数据存在误差时提取的子波有振幅畸变和相位谱扭曲的现象，另外，地震子波和震源是相关的，理论上每个地震单炮记录就有一个地震子波，确定性方法不适用于每个单炮地震子波的提取；统计性子波提取方法不需要测井信息，但需要对地震资料和地层反射系数序列的分布进行某种假设，所得到子波精度与假设条件的满足程度有关，但要对地震资料和地层反射系数的分布做某种假设、有些方法还对地震子波的相位特性有一定的要求，在实际应用中还存在一些问题，如适用性、稳定性、提取子波精度、抗噪性、子波长度限制等。

综上，目前的地震子波获取方法均存在各种不足，特别现有地震子波获取方法大多采用地震资料提取法，而地震资料提取法的核心是消除地震资料中地层反射信息的影响，但在实际应用中很难完全消除地层反射信息的影响，且基于褶积理论模型，而实际地震资料往往不能满足褶积理论模型的要求，因此急需一种能提取出高精度的地震子波的方法。

关于地震子波的概念及获取地震子波的方法的相关文献资料如下：

[1]百度百科：地震勘探、地震子波

[2]Gary F.Margrave，《Methods of Seismic Data Processing-Course LectureNotes》，The Department of Geology and Geophysics，2006

[3](苏)列.布列霍夫斯基赫，《海洋声学基础》,海洋出版社，1985.05

[4]刘伯胜，《水声学原理》,哈尔滨工业大学出版社，2010.03

[5]陈健等，《基于高阶统计量的地震子波提取方法评价》，石油地球物理勘探，2013，48(3)

[6]杨培杰，印兴耀，《地震子波提取方法综述》，石油地球物理勘探，2008，43(1)

[7]梁光河，《地震子波提取方法研究》，石油物探，1998，37(1)

[8]郑晓东，《转换波和非转换波与射线参数的关系》，地球物理学报，1991，34(6)

[9]陈浩林倪成洲等，《气枪阵列远场子波模拟与应用》，石油地球物理勘探，2008，43(6)

[10]A.Ziolkowski，《Why don’t we measure seismic signatures？》，Geophysics，1991，56(2)。

发明内容

针对现有技术的不足，本发明的目的之一提供一种地震子波的提取方法，其能够解决不依赖于褶积理论模型而能直接从地震资料中提取出地震子波；其思路是这样的：

目前对地震勘探中的地震子波的定义更多是以描述的方式进行，大致可以定义为“地震子波是由震源激发的尖脉冲经过大地滤波器滤波作用后形成的具有一定延续时间的稳定波形，是组成地震记录的基本单位”(参见：云美厚，丁伟，《石油物探》，2005年第44卷第6期)，基于对以上的认识，本发明提出地震子波仅与震源有关、不依赖于地层信息，基于这样的概念，如果可以从地震资料中分离出满足如下条件的波场信息：

1)该波场传播的地层介质信息已知或可以计算求得；

2)该波场传播与未知的地层信息无关；

3)可以建立地震子波和该波场之间的关系方程；

在满足上述条件后，就可以避免未知地层的反射信息对地震子波的影响，利用分离出的波场和地震子波之间的关系方程，计算出地震子波。

通过我们的研究和分析得出，直达波是满足上述调节的一类波场，本发明的方案也正是基于直达波来获取地震子波。

本发明的目的之二提供一种处理终端，能够解决不依赖于褶积理论模型而能直接从地震资料中提取出地震子波。

实现本发明的目的之一的技术方案为：一种地震子波的提取方法，其包括依次进行的以下步骤：

步骤S1：获取地震道中直达波s_j(t)，并经过傅里叶变换获得s_j(ω)；

步骤S2：按公式①计算得出反射系数R：

式中，R₀表示水平海面反射系数，为常数，射线参数θ为入射波的入射角，c表示水中的声波速度，为常数，ω表示入射波的圆周频率，σ表示海面波浪谱；

步骤S3：获取检波器到震源和虚震源的距离，并代入公式②，计算出中间变量B：

式中，l表示每个地震道检波器的组合个数，n表示地震道记录的个数，i为虚数单位，ω表示入射波的圆周频率，r_jk表示震源到第j个地震道中第k个检波器的直达波传播距离，r′_jk表示虚震源到第j个地震道中第k个检波器的直达波传播距离；

步骤S4：将步骤S3中的中间变量B、步骤S1中的s_j(ω)和步骤S2中的反射系数R代入公式③，计算出地震子波w(ω)：

经过公式③即计算出地震子波w(ω)。

实现本发明的目的之二的技术方案为：一种处理终端，其包括，

存储器，用于存储程序指令；

处理器，用于运行所述程序指令，以执行如下步骤：

步骤S1：获取地震道中直达波s_j(t)，并经过傅里叶变换获得s_j(ω)；

步骤S2：按公式①计算得出反射系数R：

式中，R₀表示水平海面反射系数，为常数，射线参数θ为入射波的入射角，c表示水中的声波速度，为常数，ω表示入射波的圆周频率，σ表示海面波浪谱；

步骤S3：获取检波器到震源和虚震源的距离，并代入公式②，计算出中间变量B：

式中，l表示每个地震道检波器的组合个数，n表示地震道记录的个数，i为虚数单位，ω表示入射波的圆周频率，r_jk表示震源到第j个地震道中第k个检波器的直达波传播距离，r′_jk表示虚震源到第j个地震道中第k个检波器的直达波传播距离；

步骤S4：将步骤S3中的中间变量B、步骤S1中的s_j(ω)和步骤S2中的反射系数R代入公式③，计算出地震子波w(ω)：

经过公式③即计算出地震子波w(ω)。

本发明的有益效果为：本发明从地震资料中分离出与地层信息无关的直达波，利用直达波进行计算出地震子波，使得本发明不依赖于地震褶积理论模型，不需要测井资料，不需要对地震数据和地层反射系数的分布做任何假设，也不需要对地震子波的相位特性做任何假设；同时，本发明适用于海上拖缆地震资料的地震子波提取，其计算模型简单，计算方法稳定有效，计算结果精度高、唯一性强，可以支持超长序列地震子波的提取，也可以为每一个单炮地震记录提取一个地震子波，满足地震资料批量处理的要求。

附图说明

图1为震源直达波和虚震源直达波的传播示意图；

图2为合成数据测试结果示意图；

图3为某实际测得的地震数据；

图4为用于提取震源信号远场地震子波的直达波数据；

图5为对应图4从直达波提取出的震源信号远场地震子波；

图6为采用本发明提取的地震子波合成直达波与实际记录的直达波对比图；

图7本发明较佳实施例的流程图；

图8本发明一种处理终端的结构示意图；

其中，图2至图6中的Time(s)表示时间(秒)，图6中的Preesure(mbars)表示压力(毫巴)，图2至5图中的Offset(m)表示震源与接收点间距(米)，图2和3中的Trace number表示接收点道号，图4和图5中的Shot number表示炮号。

具体实施方式

下面，结合附图以及具体实施方式，对本发明做进一步描述：

如图1至图8所示，本发明的方案是基于直达波来获取地震子波，直达波是指从震源沿介质直接传播到接收点的地震波，如果传播速度是匀速的，其传播路径是直线，在海洋地震勘探中，由于声波在水层中的传播速度可以被认为是匀速的，其直达波沿直线传播。

如图1所示，震源能量在水面的反射波被称为震源鬼波，可以看作是从虚震源(即震源关于水面的镜像点)发出的直达波，因此直达波是震源和虚震源直达波的叠合；图1中的S为震源，S’为S的镜像点，S’表示虚震源，g为检波点，d_s为震源S距离海平面的深度，d_g为检波点g距离海平面的的深度，x表示震源S与检波点g之间的水平距离，θ为海面反射入射角，r表示震源到检波点的直达波传播距离，r’为虚震源到到检波点的直达波传播距离。

从图1得出如下的关系方程式：

则震源直达波的时距曲线方程为：

式中，t表示地震波传播时间，c表示声波在水中的传播速度，为常数，通常取值为1500m/s-1540m/s。

虚震源直达波的时距曲线方程为：

式中，t′表示虚震源的地震波传播时间，c表示声波在水中的传播速度，为常数，通常取值为1500m/s-1540m/s。

根据Keller和Kolodner在1956年提出的气泡振荡理论，以及Ziolkowski et al.(1982)和Parkes et al.(1984)研究的气枪阵列特性，可知，调谐组合地震子波仅是其单个地震子波振幅的简单叠加，震源的组合地震子波信号有公式(5)：

式中，m表示震源个数，n_u(t)表示震源阵列中第u个震源的震源子波信号，表示n_u(t)沿着t轴左移个单位，t表示时间变量，r_u表示第u个检波点到震源的距离，c表示声波在水中的传播速度，p(t)表示检波点处相干叠加后的组合地震子波信号，

同样的，虚震源的组合地震子波信号有公式(6)：

式中，R表示反射系数，r′_u表示第u个检波点到虚震源的距离，p′(t)表示检波点处虚振源相干叠加后的组合地震子波信号。

而根据散射理论，可知反射系数R和入射波的圆周频率ω、海面波浪谱σ、射线参数p有如公式(7)的关系：

式中，R₀表示水平海面反射系数，为常数，通常R₀的取值为-1至-0.9之间，射线参数θ为入射波的入射角，ω为入射波的圆周频率，σ为海面波浪谱，可实际测得。

设震源阵列由m个互不相干的单震源组合而成，检波点由l个检波器并联组合而成，每炮地震数据有n道的地震道记录，u表示震源阵列中的第u个震源，j表示第j个地震道记录，k表示第k个检波器，s_j(t)表示第j个地震道的直达波，h_jk(t)表示第j个地震道中第k个检波器记录的直达波，n_i(t)表示震源阵列中第i个单震源的震源信号，因此得到公式(8)：

h_jk(t)＝p(t)+p′(t) (8)

将公式(5)和(6)代入公式(8)，得到单检波器记录的直达波和单震源信号间的关系方程为公式(9)：

式中，r_ujk表示第u个震源到第j个地震道中的第k个检波点的距离，r′_ujk表示第u个虚震源到第j个地震道中的第k个检波点的距离，

根据检波器组合理论，地震道中的直达波s_j(t)有公式(10)：

s_j(t)＝∑_kh_jk(t) (10)

将公式(9)代入公式(10),得到公式(11)：

将公式(11)进行傅里叶变换，得到频率域的表达公式(12)：

式中，ω表示入射波的圆周频率，i为虚数单位，s_j(ω)表示s_j(t)的频率域表达形式，

为便于理解公式(12)，对于n个地震道记录，可将公式(12)转换成矩阵方程(13)：

S＝AX (13)

其中，

而海上拖缆检波点的接收位置一般满足远场条件，在远场条件下，震源阵列可以被视为位于震源阵列中心的点震源，震源信号相位谱与炮检距不相关，因此方程(13)中的X向量退化为描述远场子波的单变量，记为w′(ω)，则方程(12)退化为(14)：

由公式(14)得到公式(15)：

公式(15)即为直达波计算震源远场子波的解析公式，而在实际地震数据处理中，可以使用叠加方法增强远场子波的抗干扰能力，也即是公式(16)：

公式(16)即为本发明的提取地震子波的计算公式，其中，w(ω)为本发明需要求取的地震子波，r_jk表示震源到第j个地震道中第k个检波器的直达波传播距离，r′_jk表示虚震源到第j个地震道中第k个检波器的直达波传播距离，s_j(ω)表示第j个地震道中的直达波的傅里叶变换，是s_j(t)的频率域表达形式，可从海上拖缆地震数据中直接分离得到。

综上，可以由公式(16)提取地震子波，更具体，提取地震子波的方法包括以下步骤：

步骤S1：获取地震道中直达波s_j(t)，并经过傅里叶变换获得s_j(ω)；

步骤S2：按公式①计算得出反射系数R：

式中，R₀表示水平海面反射系数，为常数，射线参数θ为入射波的入射角，c表示水中的声波速度，为常数，ω表示入射波的圆周频率，σ表示海面波浪谱；

步骤S3：获取检波器到震源和虚震源的距离，并代入公式②，计算出中间变量B：

式中，l表示每个地震道检波器的组合个数，n表示地震道记录的个数，i为虚数单位，ω表示入射波的圆周频率，r_jk表示震源到第j个地震道中第k个检波器的直达波传播距离，r′_jk表示虚震源到第j个地震道中第k个检波器的直达波传播距离；

步骤S4：将步骤S3中的中间变量B、步骤S1中的s_j(ω)和步骤S2中的反射系数R代入公式③，计算出地震子波w(ω)：

经过公式③即计算出地震子波w(ω)。

对以上方法进行验证：

a)合成数据验证：

如图2所示，采用合成了共有36道记录的直达波作为输入数据进行测试，合成记录使用了一个有52ms延迟、主频为75Hz的Ricker子波作为震源信号，接收道偏移距范围为225m至662.5m，道间距12.5m，每个接收道由16个检波器组合而成，组合间距11.719m，检波器等距组合、道内组合间距0.781m，水速1540m/s，海水面反射系数R₀使用常数-1，震源深度7m，电缆深度15m，合成后的数据如图2，图2中的图(a)表示震源信号，图(b)表示由图(a)震源信号合成36道的直达波，图(c)表示图(b)加入10％随机干扰的直达波，图(d)表示图(b)加入10％相干扰的直达波，图(e)表示由图(b)按上述方法得当的地震子波，图(f)表示由图(c)按上述方法得当的地震子波,图(g)表示由图(d)按上述方法得当的地震子波。从图中可以看出，无干扰情况下计算的结果和输入的震源信号一致；由于使用了叠加方法，加入随机干扰后得到的计算结果也有较高的信噪比，加入相干干扰的输入数据对解的影响也很小。以上说明本实施例是完全正确的。

b)实际地震数据验证：

如图3所示，为了测试利用地震数据中的直达波提取地震子波方法的有效性和实用性，我们选取了2015年5月采集于中国南海北部陆坡的某二维地震测线的部分单炮数据资料进行实际计算。野外采集使用总容量6400in3的40枪组合阵列作为震源，枪阵由四排子阵组成，每排子阵阵列长度为15m。枪沉放深度8.5m，电缆沉放深度16m。接收电缆是Sercel公司生产的Sentinel型24位数字固体电缆，每个地震道由8个检波器组成，检波器组的灵敏度为19.7μV/μbar(20℃)。

如图3和图4所示，地震数据记录长度8s，采样率2ms，480道接收，接收道偏移距范围为225m至6212.5m，道间距12.5m。图3中左上角白色矩形框圈出的是用于测试计算的数据范围。为了测试的针对性，我们回避了直达波和反射波混合的情况，选取的数据海底深度大约在1600m左右，海底反射波出现在2s以后。由枪阵组合尺寸、采样率的Nyquist频率，及1500m/s的声波速度，计算出震源信号的远场距离是37.5m，远小于地震数据的最小偏移距。选取了单炮近偏移距36道数据作为测试的输入数据，数据的偏移距范围大约在225m至663m之间，截取1.5s的数据长度可以确保计算的子波有1s的时长。

如图5所示，采用本实施例的方法，逐炮逐道计算，对于每一道数据，还要根据检波器组合参数计算出每个检波器的位置，每一炮的子波是其36道数据计算的36个子波的叠加结果，计算出的结果即如图5所示。图5是从直达波提取的震源信号远场的地震子波，10地震个子波分别对应于图4中的10个炮记录。

如图6所示，图中细虚线是实际地震记录的直达波，粗实线是采用本方法提取出的地震子波进行合成出的直达波，细虚线和粗实线只在波峰和波谷的地方未完全重合，从图6中也可以看出细虚线和粗实线几乎重叠而无法区分出，只在波峰和波谷处，细虚线和粗实线能够区分出来，说明本实施例所提出方法是可行的、利用本方法提取的地震子波是正确的。

本发明还涉及一种实现以上方法的实体装置的处理终端100，其包括，

存储器101，用于存储程序指令；

处理器102，用于运行所述程序指令，以执行如下步骤：

步骤S1：获取地震道中直达波s_j(t)，并经过傅里叶变换获得s_j(ω)；

步骤S2：按公式①计算得出反射系数R：

式中，R₀表示水平海面反射系数，为常数，射线参数θ为入射波的入射角，ω为圆周频率，σ为海面波浪谱，可实际测得；

步骤S3：获取检波器到震源和虚震源的距离，并代入公式②，计算出中间变量B：

式中，l表示每个地震道检波器的组合个数，n表示地震道记录的个数，i为虚数单位，ω表示入射波的圆周频率，r_jk表示震源到第j个地震道中第k个检波器的直达波传播距离，r′_jk表示虚震源到第j个地震道中第k个检波器的直达波传播距离；

步骤S4：将步骤S3中的中间变量B、步骤S1中的s_j(ω)和步骤S2中的反射系数R代入公式③，计算出地震子波w(ω)：

经过公式③即计算出地震子波w(ω)。

对于本领域的技术人员来说，可根据以上描述的技术方案以及构思，做出其它各种相应的改变以及变形，而所有的这些改变以及变形都应该属于本发明权利要求的保护范围之内。

Claims (2)

1.一种地震子波的提取方法，其特征在于：其包括依次进行的以下步骤：

步骤S1：获取地震道中直达波s_j(t)，并经过傅里叶变换获得s_j(ω)；

步骤S2：按公式①计算得出反射系数R：

式中，R₀表示水平海面反射系数，为常数，射线参数θ为入射波的入射角，c表示水中的声波速度，为常数，ω表示入射波的圆周频率，σ表示海面波浪谱；

步骤S3：获取检波器到震源和虚震源的距离，并代入公式②，计算出中间变量B：

式中，l表示每个地震道检波器的组合个数，n表示地震道记录的个数，i为虚数单位，ω表示入射波的圆周频率，r_jk表示震源到第j个地震道中第k个检波器的直达波传播距离，r′_jk表示虚震源到第j个地震道中第k个检波器的直达波传播距离；

步骤S4：将步骤S3中的中间变量B、步骤S1中的s_j(ω)和步骤S2中的反射系数R代入公式③，计算出地震子波w(ω)：

经过公式③即计算出地震子波w(ω)。

2.一种处理终端，其特征在于：其包括，

存储器，用于存储程序指令；

处理器，用于运行所述程序指令，以执行如下步骤：

步骤S1：获取地震道中直达波s_j(t)，并经过傅里叶变换获得s_j(ω)；

步骤S2：按公式①计算得出反射系数R：

式中，R₀表示水平海面反射系数，为常数，射线参数θ为入射波的入射角，c表示水中的声波速度，为常数，ω表示入射波的圆周频率，σ表示海面波浪谱；

步骤S3：获取检波器到震源和虚震源的距离，并代入公式②，计算出中间变量B：

式中，l表示每个地震道检波器的组合个数，n表示地震道记录的个数，i为虚数单位，ω表示入射波的圆周频率，r_jk表示震源到第j个地震道中第k个检波器的直达波传播距离，r′_jk表示虚震源到第j个地震道中第k个检波器的直达波传播距离；

步骤S4：将步骤S3中的中间变量B、步骤S1中的s_j(ω)和步骤S2中的反射系数R代入公式③，计算出地震子波w(ω)：

经过公式③即计算出地震子波w(ω)。