CN109387874B - 一种混合相位子波提取方法 - Google Patents

Abstract

本发明提出了混合相位子波提取方法，包括如下步骤：S1、获取先验子波的振幅谱；S2、获取先验子波的相位谱；S3、融合先验子波的振幅谱和先验子波的相位谱，获取先验子波；S4、在最小二乘反演中引入先验子波约束项获取混合相位子波；S5、优化混合相位子波。本申请通过结合二维相移‑时移扫描技术获取先验子波的相位的信息，利用先验子波作为约束从而有效提高提取的混合相位子波的稳定性。

Description

一种混合相位子波提取方法

技术领域

本发明涉及油气勘探中地震资料信号处理技术领域，尤其是涉及一种混合相位子波提取方法。

背景技术

地震子波为具有确定的起始时间、能量有限且具有一定连续长度的信号，它是地震记录中的基本单元，因此子波提取在地震勘探中占据着重要的位置。在正演问题中，提取的子波配合波动方程或者褶积模型来形成正演地震数据。而在地震反演\反褶积问题中，子波提取的精度决定了最终反演\反褶积的精度。

地震子波的提取包含确定性子波提取和统计性子波提取两大类。确定性子波提取方法首先需要从测井数据中获取反射系数，然后结合井旁地震道提取出子波。统计性子波提取方法需要通过地震数据本身来估计子波，因此又可以分为二阶统计量方法和高阶统计量方法，其中二阶统计量方法无法确定子波的相位，而高阶统计量方法对噪声的敏感性较强。

在已知测井资料和深时关系的情况下，通过给出子波的起始时间和子波的长度，子波的相位可以唯一确定，然而实际地震资料中由于高频噪声的影响导致地震数据的相位谱在高频段异常不稳定，因此借助于最小二乘算法提取的子波并不稳定，常常表现为子波旁瓣起伏较大。

发明内容

针对现有技术中所存在的上述技术问题，本发明提出了混合相位子波提取方法。该方法通过结合二维相移-时移扫描技术获取先验子波的相位的信息，利用先验子波作为约束从而有效提高提取的混合相位子波的稳定性。

本申请提出了混合相位子波提取方法，包括如下步骤：

确定井旁道地震数据的输入时窗范围，并在该时窗范围内获取其自相关，从而获取先验子波的振幅谱AS_wav；

根据提取子波的时窗范围确定测井数据的输入时窗范围，在深度域对测井曲线进行滤波后，根据时深关系转换至时间域，得到反射系数；根据时移-相移二维扫描技术获取先验子波的相位谱PS_wav；

融合先验子波的振幅谱AS_wav和相位谱PS_wav，并进行反傅立叶变换，得到先验子波W_prior；

结合褶积模型，通过在最小二乘反演中引入先验子波约束项以提取混合相位的子波W_ext。

通过子波旁瓣压制、去除直流分量和高频修饰以优化混合相位子波W_ext，从而得到最终的子波。

本申请针对常规最小二乘法提取的子波相位不稳定问题，通过结合二维时移-相移扫描技术获取先验子波的相位信息，利用先验子波作为约束以有效提高提取的混合相位子波的稳定性。

作为对本申请的进一步改进，所述确定井旁道地震数据的时窗范围，并在该时窗范围内获取其自相关，从而获取先验子波的振幅谱AS_wav，具体包括：

根据井口坐标确定井旁道地震数据；

确定待提取的子波的起始时间为T_w0，子波的延续时间长度为T_{w_length}；

确定井旁地震道的时间采样率为dt，子波提取的时窗范围为T_top～T_base；

则根据褶积模型的定义，确定参与提取子波的井旁道数据的输入时窗范围为T_{seis_top}～T_{seis_base}。

考虑到子波的起始时间为T_w0，因此参与获取的井旁道的起始时间为：

T_{seis_top}＝T_top+T_w0 (3)

根据褶积模型的定义可知，褶积后新信号的采样点个数等于两个输入信号采样个数之和减去1，因此参与获取的井旁道的时间总长度为：

T_{seis_length}＝T_base-T_top+T_{w_length}-dt (4)

因此：

T_{seis_base}＝T_{seis_length}+T_{seis_top}＝T_base+T_w0+T_{w_length}-dt (5)

在所述井旁道地震数据的输入时窗范围T_{seis_top}～T_{seis_base}内获取其自相关，并对自相关结果做傅立叶变换，提取自相关的能量谱，借助于子波的自相关等于井旁道的自相关的假设，则自相关的能量谱等于先验子波的能量谱：

ES_wav＝fft[Auto_cor(seis)] (6)

式中ES_wav表示先验子波的能量谱，fft[*]表示对*进行快速傅立叶变换，Seis表示井旁道，Auto_cor(*)表示求*的自相关，其中自相关的计算公式为：

其中y(t)表示自相关后序列，x(τ)表示长度为N的原始序列。

根据振幅谱和能量谱的关系，对能量谱进行开平方计算，可以得到先验子波的振幅谱AS_wav：

作为对本申请的进一步改进，根据提取子波的时窗范围确定测井数据的输入时窗范围，在深度域对测井曲线进行滤波后，根据时深关系转换至时间域，并获取反射系数，具体包括：

根据叠后或角度叠加数据，结合测井时间曲线，确定提取子波的输入时窗范围；

输入测井曲线，当井旁道为叠后地震数据时，测井数据需要读入声波时差AC、密度ρ和时间曲线time；

当井旁道为角度叠加数据时，测井数据需要读入声波时差AC、密度ρ、时间曲线time和横波时差曲线SAC，还需给定入射角α₁。

在深度域对测井曲线进行Backus(巴克斯)滤波以消除地震数据和测井数据的尺度差异，得到滤波步长；其中Backus滤波的频率是根据地震数据的奈奎斯特采样频率确定的，即：

f_backus＝f_nyq＝0.5/dt (9)

则可以得到滤波步长：

step(t)＝V(t)/f_backus＝1/[AC(t)·f_backus] (10)

确定滤波步长后，则经过Backus滤波后的声波时差、横波时差以及密度为：

AC_Backus(t)＝〈AC〉_step(t) (11)

SAC_Backus(t)＝〈SAC〉_step(t) (12)

ρ_Backus(t)＝〈ρ〉_step(t) (13)

式中AC_Backus、SAC_Backus和ρ_Backus分别表示经过Backus滤波后的声波时差、横波时差和密度，而AC、SAC和ρ分别表示原始的声波时差、横波时差和密度。符号〈*〉_step(t)表示对*进行步长为step(t)的算术平均。

根据输入的时间曲线确立井的深时关系，并根据深时关系，对经过Backus滤波后的测井曲线进行深时转换，具体包括：

利用一维线性插值算法对经过Backus滤波后的深度域的测井曲线进行重新采样到时间域：

AC_time＝interp(AC_Backus) (14)

SAC_time＝interp(SAC_Backus) (15)

ρ_time＝interp(ρ_Backus) (16)

式中，AC_Backus、SAC_Backus和ρ_Backus分别表示经过Backus滤波后的声波时差、横波时差和密度，AC_time、SAC_time和ρ_time分别表示经过时深转换后的时间域上的声波时差、横波时差和密度。

根据时间域的测井曲线计算反射系数，具体包括：

当井旁道为叠后地震数据时，则计算公式为：

I_mp＝ρ_time/AC_time (18)

当井旁道为角度叠加数据时，利用Zoeppritz方程计算PP波的反射系数R_pp，其中Zoeppritz方程的计算公式为：

其中，V_p1＝1/AC_time(t)，V_s1＝1/SAC_time(t)。

作为对本申请的进一步改进，所述根据时移-相移二维扫描技术获取先验子波的相位谱PS_wav，具体包括：

确定先验子波时移和相移的扫描范围，并针对每一对时移和相移值获取一个对应的相位谱PS_scan。具体的：

①确定相移的扫描范围为

扫描步长为

则共需要扫描

个相移值；

②确定时移的扫描范围[t_min,t_max]，扫描步长t_step，则共需要扫描

个时移值；

③通过二维循环，则共需要扫描

个相移-时移对，假设当前扫描到第m个相移值，第n个时移值记为

则对应的相位为：

其中j表示虚部，ω表示圆频率，jωt_n表示由于时移引起的相位变化，是通过傅立叶变换的时移性质推导而来的。则当前相移-时移对对应的相位谱为：

④融合所述相位谱PS_scan和先验子波的振幅谱AS_scan，得到子波的频谱FS_scan；

FS_scan＝AS_wav·PS_scan (22)

⑤将子波的频谱FS_scan进行反傅立叶变换得到线性相位的子波：

W_scan＝ifft(FS_scan) (23)

⑥将该线性相位子波和反射系数R_pp褶积获取合成道：

Syn_scan＝W_scan*R_pp (24)

⑦计算合成道与井旁道的相关系数，其中相关系数的计算公式为：

式中

和

分别表示x和y的平均值。

⑧重复第③至第⑦步，扫描完所有的时移值和相移值，统计出相关系数最大值对应的时移和相移值，从而获取先验子波的相位谱PS_wav。

获取先验子波的相位谱PS_wav后，融合先验子波的振幅谱AS_wav和相位谱PS_wav，并进行反傅立叶变换，得到先验子波W_prior；

作为对本申请的进一步改进，所述通过在最小二乘反演中引入先验子波约束项以提高子波提取的稳定性，具体包括：

给定先验子波约束项的权重因子，结合褶积模型，将权重因子加入提取子波的目标函数中，并令所述目标函数的一阶偏导数为零，从而提取混合相位的子波W_ext。

结合褶积模型，利用最小二乘算法从井旁道地震道中提取子波。已知反射系数R_pp和井旁道Seis，根据褶积模型有如下表达式：

R_pp*W_ext＝Seis (26)

其中*表示褶积过程，褶积过程可以写成矩阵相乘的形式，为了将(26)式转换成矩阵乘积形式，将向量R_pp表达为反射系数褶积矩阵G_R：

式中

表示向量R_pp的第j个元素，M表示子波的采样个数。因此有：

G_RW_ext＝Seis (28)

通过最小二乘算法可以直接求出提起的子波W_ext，但是由于地震数据噪声的存在，导致求出的混合相位子波并不稳定，主要表现在子波的旁瓣起伏较大，因此为了稳定子波的相位，通过引入先验子波约束项，构建提取子波的目标函数为：

f(W_ext)＝(G_RW_ext-Seis)^T(G_RW_ext-Seis)+λ(W_ext-W_prior)^T(W_ext-W_prior) (29)

对公式(1)求关于W_ext的一阶偏导数，并令其为0，则：

W_ext＝(G_R ^TG_R+λI)^-1(G_R ^TSeis+λW_prior) (30)

式中，T为转置，λ表示先验子波约束项权重，根据目标函数得到先验子波约束项的权重，权重越大表示提取的子波的形态与先验子波越接近，提取的子波越稳定；权重越小，表示井旁道和合成道的残差越小，子波的旁瓣起伏越大，提取的子波不稳定。

作为对本申请的进一步改进，通过对子波旁瓣压制以优化混合相位的子波W_ext，具体包括：

对混合相位的子波W_ext进行傅立叶变换得到其频谱FS_ext，并分别计算其振幅谱AS_ext和相位谱PS_ext；

对所述振幅谱AS_ext进行反傅立叶变换获取混合相位子波W_ext对应的零相位子波W₀，用Papoulis窗函数同W₀进行点积，压制W₀的子波旁瓣：

W_{0_H}＝W₀·H_p (31)

式中W_{0_H}为W₀经过旁瓣压制后的子波，H_p为Papoulis窗函数组成的向量，其中Papoulis窗函数的数学表达式为：

将经过子波旁瓣压制后的零相位子波W_{0_H}进行傅里叶变换，得到其振幅谱AS_{ext_H}，然后将其和相位谱PS_ext融合得到经过子波旁瓣压制后的混合相位子波的频谱。

FS_{ext_H}＝AS_{ext_H}·exp(PS_ext) (33)

作为对本申请的进一步改进，通过去除直流分量和高频修饰以优化混合相位子波W_ext，具体包括：

通过令混合相位子波的频谱FS_{ext_H}中的零频率成分为零，以去除直流分量；

通过对混合相位子波的频谱FS_{ext_H}进行高切频滤波，以去除高频噪音；

将经过去除直流分量和高频噪音的混合相位子波FS_{ext_H}的频谱进行反傅立叶变换，得到最终的子波。

与现有技术相比，本申请解决了借助最小二乘算法提取的子波不稳定的问题，通过结合二维相移-时移扫描技术获取先验子波的相位的信息，利用先验子波作为约束从而有效提高提取的混合相位子波的稳定性。

附图说明

下面将结合附图来对本发明的优选实施例进行详细地描述。在图中：

图1显示了根据本发明的实施例的混合相位子波提取流程图。

图2显示了本发明实施例选用的测井曲线示意图。

图3显示了本发明实施例选用的井旁道地震数据示意图。

图4显示了采用本发明的实施例所述的提取方法对图2进行Backus滤波以及深时转换后的测井曲线示意图。

图5显示了根据本发明的实施例所述的提取方法获取的井旁道自相关示意图。

图6显示了根据本发明的实施例的二维时移-相移扫描对应的合成道与井旁道的互相关值示意图。

图7显示了根据本发明的实施例所述的提取方法获取的先验子波的示意图。

图8显示了根据本发明的实施例的不同先验子波权重对应的混合相位子波示意图。

图9显示了根据本发明的实施例的经过子波旁瓣压制、去除直流分量和高频噪音前后的混合相位子波。

图10显示了根据本发明的实施例所述的提取方法获取的井震标定界面。

图11显示了根据本发明的实施例所述的提取方法获取的井震标定界面的井旁道以及合成地震记录的放大图。

图12显示了根据本发明的实施例的PetroV软件中子波提取模块的界面示意图。

图13显示了根据本发明的实施例的提取的混合相位子波以及对应的振幅谱和相位谱。

在附图中，相同的部件使用相同的附图标记。附图并未按照实际的比例绘制。

具体实施方式

下面将结合附图对本发明做进一步说明。

混合相位子波提取方法，参考图1，给出了本申请的总体流程，包括如下步骤：

确定井旁道地震数据的输入时窗范围，并在该时窗范围内计算其自相关，从而获取先验子波的振幅谱AS_wav；

根据提取子波的时窗范围确定测井数据的输入时窗范围，在深度域对测井曲线进行滤波后，根据时深关系转换至时间域，得到反射系数；根据时移-相移二维扫描技术获取先验子波的相位谱PS_wav；

融合先验子波的振幅谱AS_wav和相位谱PS_wav，并进行反傅立叶变换，得到先验子波W_prior；

结合褶积模型，通过在最小二乘反演中引入先验子波约束项以提取混合相位的子波W_ext；

通过子波旁瓣压制、去除直流分量和高频修饰以优化混合相位子波W_ext，从而得到最终的子波。

本申请针对常规最小二乘法提取的子波相位不稳定问题，通过结合二维时移-相移扫描技术获取先验子波的相位信息，利用先验子波作为约束以有效提高提取的混合相位子波的稳定性。

在本申请的实施例中，所述确定井旁道地震数据的时窗范围，并在该时窗范围内计算其自相关，从而获取先验子波的振幅谱AS_wav，具体包括：

根据井口坐标确定井旁道地震数据；

确定待提取的子波的起始时间为T_w0，子波的延续时间长度为T_{w_length}；

确定井旁地震道的时间采样率为dt，子波提取的时窗范围为T_top～T_base；

则根据褶积模型的定义，确定参与提取子波的井旁道数据的输入时窗范围为T_{seis_top}～T_{seis_base}。

考虑到子波的起始时间为T_w0，因此参与计算的井旁道的起始时间为：

T_{seis_top}＝T_top+T_w0 (3)

根据褶积模型的定义可知，褶积后新信号的采样点个数等于两个输入信号采样个数之和减去1，因此参与计算的井旁道的时间总长度为：

T_{seis_length}＝T_base-T_top+T_{w_length}-dt (4)

因此：

T_{seis_base}＝T_{seis_length}+T_{seis_top}＝T_base+T_w0+T_{w_length}-dt (5)

在所述井旁道地震数据的输入时窗范围T_{seis_top}～T_{seis_base}内计算其自相关，并对自相关结果做傅立叶变换，提取自相关的能量谱，借助于子波的自相关等于井旁道的自相关的假设，则自相关的能量谱等于先验子波的能量谱：

ES_wav＝fft[Auto_cor(seis)] (6)

式中ES_wav表示先验子波的能量谱，fft[*]表示对*进行快速傅立叶变换，Seis表示井旁道，Auto_cor(*)表示求*的自相关，其中自相关的计算公式为：

其中y(t)表示自相关后序列，x(τ)表示长度为N的原始序列。

根据振幅谱和能量谱的关系，对能量谱进行开平方计算，可以得到先验子波的振幅谱AS_wav：

在本申请的实施例中，根据提取子波的时窗范围确定测井数据的输入时窗范围，在深度域对测井曲线进行滤波后，根据时深关系转换至时间域，并计算反射系数，具体包括：

根据叠后或角度叠加数据，结合测井时间曲线，确定提取子波的输入时窗范围；

输入测井曲线，当井旁道为叠后地震数据时，测井数据需要读入声波时差AC、密度ρ和时间曲线time；

当井旁道为角度叠加数据时，测井数据需要读入声波时差AC、密度ρ、时间曲线time和横波时差曲线SAC，还需给定入射角α₁。

图2为验证本申请所涉及的测井资料，从左到右依次是声波时差、横波时差、密度以及深时关系，纵轴为TVD，深度间隔为0.125m。

图3为验证本申请所涉及的角度叠加数据，其中中心角为18度，提取子波的时窗范围为1300ms～1900ms，子波的起始时间为-50ms，子波的时间长度为100ms，因此输入的角度叠加地震数据的时间范围为1250ms～1848ms，井旁道为第6CDP处。

由于地震数据和测井数据在垂向上的尺度差异，在将测井曲线进行深时转换前，需要先对测井曲线在深度域进行Backus滤波，计算测井曲线在地震尺度的等效弹性性质，图4中左边三附图中蓝色曲线1、蓝色曲线3和蓝色曲线5依次为经过Backus滤波后的声波时差、横波时差和密度，红色虚线2、红色虚线4和红色虚线6分别对应经过深时转换后的声波时差、横波时差和密度，通过在时间域叠置，可以看出经过Backus滤波后，测井曲线在深度域与时间域较为接近，尺度差异能够较为明显的消除，而图2中的原始曲线从视觉效果上来看，要比时间域的测井曲线稠密，存在较为明显的尺度差异。在得到时间域的测井曲线后，结合Zoeppritz方程，计算中心角为18度的角度反射系数，其反射系数在图4中第4道所示。

在深度域对测井曲线进行Backus滤波以消除地震数据和测井数据的尺度差异，得到滤波步长；其中Backus滤波的频率是根据地震数据的奈奎斯特采样频率确定的，即：

f_backus＝f_nyq＝0.5/dt (9)

则可以得到滤波步长：

step(t)＝V(t)/f_backus＝1/[AC(t)·f_backus] (10)

确定滤波步长后，则经过Backus滤波后的声波时差、横波时差以及密度为：

AC_Backus(t)＝<AC>_step(t) (11)

SAC_Backus(t)＝＜SAC>_step(t) (12)

ρ_Backus(t)＝<ρ>_step(t) (13)

式中AC_Backus、SAC_Backus和ρ_Backus分别表示经过Backus滤波后的声波时差、横波时差和密度，而AC、SAC和ρ分别表示原始的声波时差、横波时差和密度。符号<*>_step(t)表示对*进行步长为step(t)的算术平均。

根据输入的时间曲线确立井的深时关系，并根据深时关系，对经过Backus滤波后的测井曲线进行深时转换，具体包括：

利用一维线性插值算法对经过Backus滤波后的深度域的测井曲线进行重新采样到时间域：

AC_time＝interp(AC_Backus) (14)

SAC_time＝interp(SAC_Backus) (15)

ρ_time＝interp(ρ_Backus) (16)

式中，AC_Backus、SAC_Backus和ρ_Backus分别表示经过Backus滤波后的声波时差、横波时差和密度，AC_time、SAC_time和ρ_time分别表示经过时深转换后的时间域上的声波时差、横波时差和密度。

根据时间域的测井曲线计算反射系数，具体包括：

当井旁道为叠后地震数据时，则计算公式为：

I_mp＝ρ_time/AC_time (18)

当井旁道为角度叠加数据时，利用Zoeppritz方程计算PP波的反射系数R_pp，其中Zoeppritz方程的计算公式为：

其中，V_p1＝1/AC_time(t)，V_s1＝1/SAC_time(t)。

在本申请的实施例中，所述根据时移-相移二维扫描技术获取先验子波的相位谱PS_wav，具体包括：

确定先验子波时移和相移的扫描范围，并针对每一对时移和相移值获取一个对应的相位谱PS_scan。具体的：

①确定相移的扫描范围为

扫描步长为

则共需要扫描

个相移值；

②确定时移的扫描范围[t_min,t_max]，扫描步长t_step，则共需要扫描

个时移值；

③通过二维循环，则共需要扫描

个相移-时移对，假设当前扫描到第m个相移值，第n个时移值记为

则对应的相位为：

其中j表示虚部，ω表示圆频率，jωt_n表示由于时移引起的相位变化，是通过傅立叶变换的时移性质推导而来的。则当前相移-时移对对应的相位谱为：

④融合所述相位谱PS_scan和先验子波的振幅谱AS_scan，得到子波的频谱FS_scan；

FS_scan＝AS_wav·PS_scan (22)

⑤将子波的频谱FS_scan进行反傅立叶变换得到线性相位的子波：

W_scan＝ifft(FS_scan) (23)

⑥将该线性相位子波和反射系数R_pp褶积获取合成道：

Syn_scan＝W_scan*R_pp (24)

⑦计算合成道与井旁道的相关系数，其中相关系数的计算公式为：

式中

和

分别表示x和y的平均值。

⑧重复第③至第⑦步，扫描完所有的时移值和相移值，统计出相关系数最大值对应的时移和相移值，从而获取先验子波的相位谱PS_wav。

获取先验子波的相位谱PS_wav后，融合先验子波的振幅谱AS_wav和相位谱PS_wav，并进行反傅立叶变换，得到先验子波W_prior；

将井旁道时窗范围1250ms～1848ms内的地震数据输入，计算其自相关，并取其-50ms到50ms参与计算先验子波的振幅谱，将先验子波的振幅谱进行反傅里叶变换便可以得到先验子波对应的零相位子波，如图5所示。

优选地，确定相移的扫描范围[-90°,90°]，扫描步长为1°，时移的扫描范围[-20,20]，扫描步长为2ms，通过扫描一系列的线性相位子波与反射系数褶积得到的合成道与井旁道的互相关系数，可以得到互相关系数最大值对应的相移和时移值，从而可以确定出先验子波的线性相位，结合先验子波的振幅谱便可以得到先验子波的形态，图6给出了不同相移、时移对应的线性相位子波与反射系数褶积的合成道与井旁道的互相关系数分布图，从分布图中可以确定出互相关值最大为相移42度，时移2ms处，因此可以得到先验子波的形态如图7所示。

在得到先验子波后，给定先验子波约束项权重系数，通过最小二乘反演的思想可以从井旁道中提取出混合相位子波，图8给出了基于不同先验子波权重系数提取的混合相位子波。其中点线的先验子波权重为0，实线的先验子波权重为0.1，虚线的先验子波权重为2。可以看出随着权重系数的增加，提取的混合相位子波的旁瓣抖动越小，形状也更加接近先验子波。

在本申请的实施例中，所述通过在最小二乘反演中引入先验子波约束项以提高子波提取的稳定性，给定先验子波约束项的权重因子，结合褶积模型，将权重因子加入提取子波的目标函数中，并令所述目标函数的一阶偏导数为零，从而提取混合相位的子波W_ext。具体包括：

结合褶积模型，利用最小二乘算法从井旁道地震道中提取子波。已知反射系数R_pp和井旁道Seis，根据褶积模型有如下表达式：

R_pp*W_ext＝Seis (26)

其中*表示褶积过程，褶积过程可以写成矩阵相乘的形式，为了将(26)式转换成矩阵乘积形式，将向量R_pp表达为反射系数褶积矩阵G_R：

式中

表示向量R_pp的第j个元素，M表示子波的采样个数。因此有：

G_RW_ext＝Seis (28)

通过最小二乘算法可以直接求出提起的子波W_ext，但是由于地震数据噪声的存在，导致求出的混合相位子波并不稳定，主要表现在子波的旁瓣起伏较大，因此为了稳定子波的相位，通过引入先验子波约束项，构建提取子波的目标函数为：

f(W_ext)＝(G_RW_ext-Seis)^T(G_RW_ext-Seis)+λ(W_ext-W_prior)^T(W_ext-W_prior) (29)

对公式(1)求关于W_ext的一阶偏导数，并令其为0，则：

W_ext＝(G_R ^TG_R+λI)^-1(G_R ^TSeis+λW_prior) (30)

式中，T为转置，λ表示先验子波约束项权重，根据目标函数得到先验子波约束项的权重，权重越大表示提取的子波的形态与先验子波越接近，提取的子波越稳定；权重越小，表示井旁道和合成道的残差越小，子波的旁瓣起伏越大，提取的子波不稳定。

在本申请的实施例中，通过对子波旁瓣压制以优化混合相位的子波W_ext，具体包括：

对混合相位的子波W_ext进行傅立叶变换得到其频谱FS_ext，并分别计算其振幅谱AS_ext和相位谱PS_ext；

对所述振幅谱AS_ext进行反傅立叶变换获取混合相位子波W_ext对应的零相位子波W₀，用Papoulis窗函数同W₀进行点积，压制W₀的子波旁瓣：

W_{0_H}＝W₀·H_p (31)

式中W_{0_H}为W₀经过旁瓣压制后的子波，H_p为Papoulis窗函数组成的向量，其中Papoulis窗函数的数学表达式为：

将经过子波旁瓣压制后的零相位子波W_{0_H}进行傅里叶变换，得到其振幅谱AS_{ext_H}，然后将其和相位谱PS_ext融合得到经过子波旁瓣压制后的混合相位子波的频谱。

FS_{ext_H}＝AS_{ext_H}·exp(PS_ext) (33)

在本申请的实施例中，通过去除直流分量和高频修饰以优化混合相位子波W_ext，具体包括：

通过令混合相位子波的频谱FS_{ext_H}中的零频率成分为零，以去除直流分量；

通过对混合相位子波的频谱FS_{ext_H}进行高切频滤波，以去除高频噪音；

将经过去除直流分量和高频噪音的混合相位子波FS_{ext_H}的频谱进行反傅立叶变换，得到最终的子波。

图9给出了子波旁瓣压制、去除直流分量和高频噪音处理前后的混合相位子波(先验子波约束项权重为0.1)，通过对比可以看出经过子波旁瓣压制、去除直流分量和高频噪音处理后，提取的混合相位子波旁瓣抖动减小并且子波形态更加稳定光滑。

图10为集成了该子波提取算法的PetroV软件中的井震标定界面，该界面中从左到右依次为声波时差曲线、横波时差曲线、密度曲线、井旁道以及合成地震记录，图11为井旁道以及合成地震记录的放大图。设置子波的提取时窗为Inv_Hor_re1510到Inv_Hor_re2005之间，井旁道为小角度叠加数据体，其中心角度为6度，首先用主频为30Hz雷克子波产生合成道，通过井旁道和合成道波组关系对比来调整时深关系。在获取一个准确的时深关系的前提下，设置子波提取的起始时间为-50ms，子波的长度为100ms，高截频为80Hz，先验子波约束项权重因子为0.1，见图12，提取的混合相位子波如图13所示，图13中下方分别为混合相位子波的振幅谱和相位谱，其中相位谱值域在-8度到80度之间，可以看出相位较为稳定。

与现有技术相比，本申请解决了借助最小二乘算法提取的子波不稳定的问题，通过结合二维相移-时移扫描技术获取先验子波的相位的信息，利用先验子波作为约束从而有效提高提取的混合相位子波的稳定性。

以上所述仅为本发明的优选实施方式，但本发明保护范围并不局限于此，任何本领域的技术人员在本发明公开的技术范围内，可容易地进行改变或变化，而这种改变或变化都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求书的保护范围为准。

Claims (7)

1.一种混合相位子波提取方法，其特征在于，包括如下步骤：

S1、获取先验子波的振幅谱；

S2、获取先验子波的相位谱；

S3、融合先验子波的振幅谱和先验子波的相位谱，获取先验子波；

S4、根据先验子波约束项获取混合相位子波；

所述步骤S2具体包括：

S21、获取测井数据的输入时窗范围，在深度域对测井曲线进行滤波后，根据时深关系转换至时间域，并获取反射系数；

S22、根据时移-相移二维扫描技术获取先验子波的相位谱；

所述步骤S22具体包括：

S221、确定先验子波时移和相移的扫描范围，并针对每一对时移和相移值获取一个对应的相位谱；

S222、将相位谱和先验子波的振幅谱相融合，得到子波的频谱；

S223、将子波的频谱进行傅立叶反变换得到线性相位的子波，并和反射系数褶积获取合成道，获取合成道与井旁道的相关系数；

S224、扫描完所有的时移值和相移值，统计出相关系数最大值对应的时移和相移值，从而获取先验子波的相位谱。

2.根据权利要求1所述的混合相位子波提取方法，其特征在于，所述步骤S1具体包括：

S11、根据井口坐标确定井旁道地震数据；

S12、根据子波的起始时间、子波的时间长度以及提取子波的时窗范围确定井旁道地震数据的输入时窗范围；

S13、在所述井旁道地震数据的输入时窗范围内获取其自相关，并对自相关结果做傅立叶变换，得到井旁道的自相关的振幅谱，并根据其振幅谱和能量谱的关系获取先验子波的振幅谱。

3.根据权利要求1所述的混合相位子波提取方法，其特征在于，所述步骤S21具体包括：

S211、根据叠后或角度叠加数据，结合测井时间曲线，确定提取子波的输入时窗范围；

S212、输入测井曲线，在深度域对测井曲线进行滤波以消除地震数据和测井数据的尺度差异，得到滤波步长；

S213、根据输入的时间曲线确立井的深时关系，并根据深时关系，对经过滤波后的测井曲线进行深时转换，得到时间域的测井曲线；

S214、根据时间域的测井曲线获取反射系数。

4.根据权利要求1所述的混合相位子波提取方法，其特征在于，所述步骤S4具体包括：

结合褶积模型，将先验子波约束项的权重因子引入提取子波的目标函数，以提取混合相位的子波。

5.根据权利要求4所述的混合相位子波提取方法，其特征在于，所述目标函数为：

f(W_ext)＝(G_RW_ext-Seis)^T(G_RW_ext-Seis)+λ(W_ext-W_prior)^T(W_ext-W_prior) (1)

令所述目标函数的一阶偏导数为零，从而得到混合相位的子波W_ext：

W_ext＝(G_R ^TG_R+λI)^-1(G_R ^TSeis+λW_prior) (2)

式中，λ为先验子波约束项权重因子，Seis为井旁道，G_R为反射系数褶积矩阵，T为转置，I为单位矩阵，W_prior为先验子波，W_ext为混合相位子波。

6.根据权利要求1所述的混合相位子波提取方法，其特征在于，所述提取方法还包括步骤S5：优化混合相位子波。

7.根据权利要求6所述的混合相位子波提取方法，其特征在于，所述步骤S5中，通过子波旁瓣压制、去除直流分量和高频修饰的方式优化混合相位子波。

Images