CN116840916B

CN116840916B - 一种地震速度信号和加速度信号联合子波提取方法

Info

Publication number: CN116840916B
Application number: CN202310807076.4A
Authority: CN
Inventors: 游李; 杨宇勇; 罗川; 汪政洋
Original assignee: Chengdu Univeristy of Technology
Current assignee: Chengdu Univeristy of Technology
Priority date: 2023-07-04
Filing date: 2023-07-04
Publication date: 2024-03-26
Anticipated expiration: 2043-07-04
Also published as: CN116840916A

Abstract

本发明公开了一种地震速度信号和加速度信号联合子波提取方法，包括以下步骤：获取地震的加速度信号和速度信号，并对其进行标准化处理；进行常相位子波和最小相位子波的提取或者进行混合相位子波的提取；常相位子波和最小相位子波的提取具体包括以下子步骤：S1：基于速度信号提取获得最小相位子波一；S2：基于加速度信号提取获得常相位子波，并对其进行转换，获得最小相位子波二；S3：计算两个最小相位子波之间的残差；若残差大于残差阈值，则返回步骤S2，重新进行常相位子波提取；反之则输出最小相位子波一和常相位子波。本发明通过速度信号与加速度信号的一致性联合提取出不同信号类型的子波，实现了地震速度信号和加速度信号的子波联合提取。

Description

一种地震速度信号和加速度信号联合子波提取方法

技术领域

本发明涉及石油地球物理勘探技术领域，特别涉及一种地震速度信号和加速度信号联合子波提取方法。

背景技术

石油和天然气是国家重要的战略能源，但近年来，随着石油对外依存度的不断提高，石油和天然气的供应和安全问题日益突出，加强技术研究已成为当务之急。反射式地震勘探是利用反射波资料来获得地层结构与岩性的重要手段，而提取子波的准确性严重制约着认识地层结构与岩性的准确度。

目前，进行子波提取时，通常仅采用一种信号进行子波提取，如此对选取信号以及相位校正因子的依赖性大，因此，亟需一种能够提高子波提取准确性的子波提取方法。

发明内容

针对上述问题，本发明旨在提供一种地震速度信号和加速度信号联合子波提取方法。

本发明的技术方案如下：

一种地震速度信号和加速度信号联合子波提取方法，包括以下步骤：

获取地震的加速度信号和速度信号；

对所述加速度信号和所述速度信号进行标准化处理；

基于子波的相位类型需求，设定参数，进行常相位子波和最小相位子波的提取或者进行混合相位子波的提取。

作为优选，进行常相位子波和最小相位子波的提取具体包括以下子步骤：

S1：基于所述速度信号进行最小相位子波提取，获得最小相位子波一；

S2：基于所述加速度信号进行常相位子波提取，获得常相位子波，并对所述常相位子波进行转换，获得最小相位子波二；

S3：计算所述最小相位子波一与所述最小相位子波二之间的残差；

若所述残差大于残差阈值，则返回步骤S2，重新进行常相位子波提取；

若所述残差小于等于残差阈值，则输出所述最小相位子波一和所述常相位子波。

作为优选，步骤S1中，基于所述速度信号进行最小相位子波提取时，采用下式进行提取：

a(n)＝u(n)a_o(n) (1)

a_o(n)＝IFFT(A_o(ω)) (3)

A_o(ω)＝ln|X(ω)| (4)

式中：a(n)为最小相位子波；u(n)为系数；a_o(n)为a(n)的奇部序列；n为时间序列；IFFT表示求取傅里叶逆变换；A_o(ω))为a_o(n)的傅里叶变换；|X(ω)|为地震记录的振幅谱。

作为优选，步骤S2中，基于所述加速度信号进行常相位子波提取时，采用下式进行提取：

g(a，n)＝b(n)cos(a)-h_b(n)sin(a) (5)

b(n)＝IFFT(|X(ω)|) (6)

式中：g(a，n)为常相位子波；b(n)为零相位子波；a为相位校正因子；h_b(n)为b(n)的Hilbert变换。

作为优选，所述相位校正因子为：

式中：k为正整数。

作为优选，进行混合相位子波的提取时，采用下式进行提取：

w＝IFFT(W) (8)

W＝sin(Φ′)|W|+jcos(Φ′)|W| (9)

式中：w为混合相位子波；W为子波的频谱；Φ′为子波的相位谱；|W|为子波的振幅谱；j为虚数。

作为优选，子波的相位谱Φ′通过下式进行计算：

式中：A_pha为(N²/4)×N的矩阵；T表示矩阵的转置；-1表示矩阵的逆；为(N²/4)×1的列向量；N为高阶累积量矩阵的大小。

作为优选，子波的振幅谱|W|通过下式进行计算：

式中：|W(n)|为子波振幅谱；e为自然底数；W′(n)为子波振幅谱的自然对数；|W(N-n)|为子波振幅谱大于二分之一倍子波长度的部分；W′为(N/2)×1的列向量；A_amp为(N²/16)×(N/2)的矩阵；X′为(N²/16)×1的列向量。

本发明的有益效果是：

本发明通过联合地震速度信号和加速度信号提取子波，其能够提高地震信号的子波提取准确性，解决现有技术仅用一种信号进行子波提取对选取信号以及相位校正因子依赖性大的问题。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明地震速度信号和加速度信号联合子波提取方法的流程示意图；

图2为实施例1中的地震信号示意图；其中，图2(a)为速度信号，图2(b)为加速度信号；

图3为实施例1中本发明联合提取出的速度子波和加速度子波结果示意图；其中，图3(a)为速度子波，图3(b)为加速度子波；

图4为实施例1中本发明联合提取出的混合相位子波结果示意图；

图5为实施例2中的地震信号示意图；其中，图5(a)为速度信号，图5(b)为加速度信号；

图6为实施例2中本发明联合提取出的速度子波和加速度子波结果示意图；其中，图6(a)为速度子波，图6(b)为加速度子波；

图7为实施例2中现有技术单独提取出的速度子波和加速度子波结果示意图；其中，图7(a)为速度子波，图7(b)为加速度子波；

图8为实施例2中本发明联合提取出的混合相位子波结果示意图；

图9为实施例2中现有技术提取出的混合相位子波结果示意图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的技术特征可以相互结合。需要指出的是，除非另有指明，本申请使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。本发明公开使用的“包括”或者“包含”等类似的词语意指出现该词前面的元件或者物件涵盖出现在该词后面列举的元件或者物件及其等同，而不排除其他元件或者物件。

如图1所示，本发明提供一种地震速度信号和加速度信号联合子波提取方法，包括以下步骤：获取地震的加速度信号和速度信号；对所述加速度信号和所述速度信号进行标准化处理；基于子波的相位类型需求，设定参数，进行常相位子波和最小相位子波的提取或者进行混合相位子波的提取。

加速度信号对比速度信号，在频带宽度、倍频程等方面有巨大优势；速度信号对比加速度信号在低频部分具有有效信息更丰富等特点。本发明联合速度信号与加速度信号进行子波提取，不仅能保留频带宽度等因素，也能保留低频有效信号，提高子波提取结果的准确性。

需要说明的是，对信号进行标准化处理为现有技术，具体步骤在此不再赘述。设定参数时，包括设定处理道数、时窗等，其也为现有技术，具体在此不再赘述。

在一个具体的实施例中，进行常相位子波和最小相位子波的提取具体包括以下子步骤：

S1：基于所述速度信号进行最小相位子波提取，获得最小相位子波一。

在一个具体的实施例中，基于所述速度信号进行最小相位子波提取时，采用下式进行提取：

a(n)＝u(n)a_o(n) (1)

a_o(n)＝IFFT(A_o(ω)) (3)

A_o(ω)＝ln|X(ω)| (4)

上述实施例中，最小相位子波提取公式通过以下步骤推导得到：

假设地震子波为b(n)，子波的对数谱A(ω)可写为：

A(ω)＝lnB(ω) (17)

式中：B(ω)为地震子波的频谱；ln表示求以e为底的对数。

假设a(n)为地震子波b(n)的对数谱时间序列(即最小相位子波)。理论上证明，子波若为最小相位子波，那么其对数谱的时间序列是实的因果序列。任何实序列都可以写为其奇部序列与偶部序列之和，故a(n)可写为：

a(n)＝a_o(n)+a_e(n) (18)

式中：a_e(n)为对数谱时间序列a(n)的偶部序列。

由于a(n)具有因果性，a(n)又可写为式(1)所示的表达式，式(1)中的u(n)如式(2)所示，因此，只需得到奇部序列a_o(n)，通过公式(1)就可获得最小相位子波。

对式(31)左右两边取对数，可得：

ln|B(ω)|＝ln|X(ω)| (19)

同时，子波的频谱可以表示为：

式中：为子波的相位。

那么，将式(19)代入式(17)，可得：

同时，有：

A_o(ω)＝A_R(ω) (22)

A_e(ω)＝A_I(ω) (23)

A(ω)＝A_o(ω)+iA_e(ω) (24)

式中：A_o(ω)和A_e(ω)分别a_o(n)和a_e(n)的傅里叶变换；A_R(ω)和A_I(ω)分别为子波对数谱A(ω)的实部和虚部，i为虚数单位。

联立式(21)-(24)，可得：

A_o(ω)＝ln|B(ω)| (25)

即：

A_o(ω)＝ln|B(ω)|＝ln|X(ω)| (26)

从而得到式(4)所示的A_o(ω)的计算公式，并由此得到式(3)所示的a_o(n)的计算公式。最后，根据n的取值确定u(n)，并通过式(3)计算出奇部序列a_o(n)就可通过式(1)提取出基于速度信号的最小相位子波。

S2：基于所述加速度信号进行常相位子波提取，获得常相位子波，并对所述常相位子波进行转换，获得最小相位子波二。

在一个具体的实施例中，基于所述加速度信号进行常相位子波提取时，采用下式进行提取：

g(a,n)＝b(n)cos(a)-h_b(n)sin(a) (5)

b(n)＝IFFT(|X(ω)|) (6)

可选地，所述相位校正因子为：

式中：k为正整数。

上述实施例中，常相位子波提取公式通过以下步骤推导得到：

假设地震记录为x(n)，反射系数为r(n)，地震子波为b(n)(当子波相位为零时所述地震子波即为零相位子波)。那么，地震褶积模型可表达为：

x(n)＝b(n)*r(n) (27)

首先，对式(27)两边均进行傅里叶变换，即可得到褶积模型在频域上的表示：

X(ω)＝B(ω)*R(ω) (28)

式中：X(ω)为x(n)的傅里叶变换，B(ω)为b(n)的傅里叶变换，R(ω)为r(n)的傅里叶变换。

假设反射系数为白噪序列，则其振幅谱为：

|R(ω)|＝1 (29)

那么，地震记录x(n)的振幅谱可写为：

|X(ω)|＝|B(ω)||R(ω)|＝|B(ω)| (30)

即：

|B(ω)|＝|X(ω)| (31)

式(31)表明，可以通过地震记录的振幅谱|X(ω)|估计出子波的振幅谱|B(ω)|。

同时，子波的频谱可以表示为：

式中：为子波的相位。

其次，当子波相位为零时，即式(32)中的可得：

b(n)＝IFFT(B(ω))＝IFFT(|B(ω)|)＝IFFT(|X(ω)|) (33)

由此可以得到式(6)所示的零相位子波计算公式，通过地震记录的振幅谱|X(ω)|进行傅里叶逆变换后即可获得零相位子波b(n)。

再次，对零相位子波进行相位校正，可获得常相位子波在频域的表达式：

式中：G(a，ω)为常相位子波的频谱。

接下来，假设Hilbert变换因子的频率响应为H(ω)，则：

将式(35)代入式(34)，可得：

G(a，ω)＝B(ω)cos(a)-H(ω)sin(a) (36)

由此可以得到式(5)所示的常相位子波提取公式，利用该公式即可提取出加速度信号的常相位子波。

S3：计算所述最小相位子波一与所述最小相位子波二之间的残差；若所述残差大于残差阈值，则返回步骤S2，重新进行常相位子波提取；若所述残差小于等于残差阈值，则输出所述最小相位子波一和所述常相位子波。

需要说明的是，当所述残差大于残差阈值返回步骤S2时，通过重新设置新的相位校正因子，从而提取获得新的常相位因子；重新设置新的相位校正因子时，可采用式(7)所示的相位校正因子计算公式进行计算获得，每更新一次相位校正因子，k的取值往后取一位。

需要说明的是，所述残差阈值为人为设定阈值，所述残差阈值越小，获得的结果精度越高。本发明利用地震加速度信号等于地震速度信号的导数，以速度信号提取获得的最小相位子波一为基准，通过对加速度信号提取获得的常相位子波进行转换，利用其转换得到的最小相位子波二与所述最小相位子波一进行比较，从而确定常相位因子，避免人为设定常相位因子带来的不确定性，最终提高常相位子波的提取精度。

在一个具体的实施例中，进行混合相位子波的提取时，采用下式进行提取：

w＝IFFT(W) (8)

w＝sin(Φ′)|W|+jcos(Φ′)|W| (9)

可选地，子波的相位谱Φ′通过下式进行计算：

可选地，子波的振幅谱|W|通过下式进行计算：

在上述实施例中，利用加速度信号和速度信号的高阶累积量提取混合相位子波，其能够满足子波经过衰减后的特征，提取结果更符合实际。

上述实施例中，混合相位子波提取公式通过以下步骤推导得到：

假设存在一个零均值平稳随机过程y(n)，其中n＝0，±1，±2，...。

则y(n)三阶累计量可表示为：

cum_3y(τ₁，τ₂)＝E[y(n)y(n+τ₁)y(n+τ₂)] (37)

式中：E[.]表示计算数学期望。

上式可改写为：

k阶累积量的k-1次傅里叶变换即为其的k阶谱，也称为k-1谱。那么，上述三阶累计量的三阶谱(双谱)可表示为：

将式(38)代入式(39)，可得：

即：

式中：Y为y(n)的傅里叶变换。

那么，三阶谱的振幅谱和相位谱可分别表示为：

|cum_3y(ω₁，ω₂)|＝|Y(ω₁)||Y(ω₂)||Y(ω₁+ω₂)| (42)

假设地震记录为x(t)，则褶积模型可以表示为：

x(t)＝w(t)*r(t)+v(t) (44)

式中：w(t)为地震子波序列；r(t)为反射系数序列；v(t)为高斯噪音序列；*表示进行褶积。

将式(44)进行傅里叶变换，可得：

X(ω)＝W(ω)*R(ω)+V(ω) (45)

式中：X(ω)，W(ω)，R(ω)，V(ω)分别为x(t)，w(t)，r(t)，v(t)的傅里叶变换。

那么，x(t)的双谱为：

C_3x(ω₁，ω₂)＝C_3w(ω₁，ω₂)C_3r(ω₁，ω₂)+C_3v(ω₁，ω₂) (46)

因为高斯随机变量的高阶累积量为零，则：

C_3v(ω₁，ω₂)＝0 (47)

因此，将公式(47)代入公式(46)，可得：

C_3x(ω₁，ω₂)＝C_3w(ω₁，ω₂)C_3r(ω₁，ω₂) (48)

通常，假设反射系数的高阶累积量为多维冲击函数，因此，可将式(48)写为：

C_3x(ω₁，ω₂)＝C_3w(ω₁，ω₂)/σ² (49)

式中：σ为反射系数序列的方差。

公式(49)表明子波的三阶谱和地震记录的三阶谱之间只差了一个比例系数。因此，可利用地震记录的三阶谱确定子波的频谱。

对式(42)等式左右两边取对数，可得：

ln|cum_3y(ω₁，ω₂)|＝ln|Y(ω₁)|+ln|Y(ω₂)|+ln|Y(ω₁+ω₂)| (50)

将式(50)应用到地震数据中，可得：

ln|cum_3x(ω₁，ω₂)|＝ln|W(ω₁)|+ln|W(ω₂)|+ln|W(ω₁+ω₂)| (51)

令ln|cum_3x(ω₁，ω₂)|＝X′(k，l)，ln|W|＝W′，式(51)可写为：

X′(k，l)＝W′(k)+W′(l)+W′(k+l) (52)

式(52)表明X′是W′的线性组合，因此可将式(52)写为矩阵形式：

X′＝A_ampW′ (53)

其中，X′和A_amp的具体形式分别如式(16)和式(15)所示。

A_amp为列满秩的矩阵，因此，可用广义逆矩阵进行求解公式(53)，解得式(14)所示的W′的表达式，然后利用式(13)即可求解出子波的振幅谱。

将式(13)应用到地震数据中，可得：

式(54)表明是Φ_w的线性组合，因此可将式(54)写为矩阵形式：

式中：Φ′为N×1的列向量。

式(45)中，和A_pha的具体形式分别如式(12)和式(11)所示。

A_pha为列满秩的矩阵，因此，可用广义逆矩阵进行求解式(55)，解得式(10)所示的子波的相位谱的表达式，然后利用式(10)即可求解出子波的相位谱。

利用上述相位恢复的方法利用所有可能得到的双谱相位值，是非递推算法，因此在计算过程中没有累计误差。

根据式(10)和式(13)的计算结果，利用式(9)即可计算获得子波的频谱。然后对式(9)进行傅里叶逆变换，即可得到式(8)所示的混合相位子波提取公式，通过式(8)即可求解获得地震记录的混合相位子波。

在一个具体的实施例中，采用本发明所述地震速度信号和加速度信号联合子波提取方法对目标地震信号进行子波提取，具体包括以下步骤：

(1)获取地震的加速度信号和速度信号，并对所述速度信号和所述加速度信号进行标准化处理；

(2)设定参数(处理道数、时窗等)，对时窗范围内的数据进行平均化，获得平均地震记录及其振幅谱；

(3)对加速度信号获得的振幅谱进行傅里叶逆变换，获得零相位子波；

(4)对所述零相位子波进行希尔伯特变换并设置常相位因子，提取出加速度子波(即所述常相位子波)；

(5)对速度信号获得的振幅谱进行对数操作，获得对数谱序列；

(6)将对数谱进行傅里叶逆变换，获得频谱，然后提取出速度子波(即所述最小相位子波一)；

(7)将加速度子波转换为速度子波(即所述最小相位子波二)，计算转换的速度子波和原始速度子波之间的残差；若残差满足精度需求，则停止计算，输出联合提取的速度子波和加速度子波；若不满足精度要求，则重新选择相位校正因子，重复上述操作，直至满足精度需求；

(8)在提取出常相位子波和最小相位子波之后，将预处理(标准化、设定参数、平均化)之后的数据转换到加速度域，并与预处理之后的加速度信号进行拼接，获得融合后的数据；

(9)计算融合后的数据的三阶累积量及其三阶谱(双谱)，利用parzen窗对三阶累积量进行滤波，之后利用双谱计算子波的振幅谱和相位谱，最终，重构成混合相位子波。

在实施例1中，选取了某工区如图2所示的地震信号，时间范围为0-2s，道数范围为1-853道，将残差阈值设为0.001，最终输出联合提取出的速度子波和加速度子波如图3所示，利用高阶累积量的双谱重构出的混合相位子波如图4所示。

在实施例2中，选取了某工区如图5所示marmousi2模型的地震信号，时间范围为0-1.2s，道数范围为1-3500道，将残差阈值设为0.001，最终输出联合提取出的速度子波和加速度子波如图6所示，利用高阶累积量的双谱重构出的混合相位子波如图8所示。

在对上述Marmousi2模型进行子波提取时，不采用本发明的联合子波提取方法，在其他参数(处理道数、时窗范围)均一致的情况下，采用现有技术单独提取速度子波和加速度子波，结果如图7所示，采用现有技术(不加滤波器平滑三阶累积量)提取出的混合相位子波结果如图9所示。

对比图6和图7可知，现有技术单独提取的速度子波与本发明联合提取的速度子波基本一致，但是单独提取的加速度子波在振幅强弱等方面上均比本发明联合提取的加速度子波弱。

对比图8和图9可知，本发明联合提取得到的混合相位子波比现有技术单独提取得到的混合相位子波更接近理论子波，精确度更高。

综上所述，本发明一方面能够提供常相位子波和最小相位子波这种理论子波的提取，也能提取满足子波衰减后特征的混合相位子波，且两种提取方法均联合了地震速度信号和加速度信号，其能够提高提取结果的准确性。与现有技术相比，本发明具有显著的进步。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims

1.一种地震速度信号和加速度信号联合子波提取方法，其特征在于，包括以下步骤：

获取地震的加速度信号和速度信号；

对所述加速度信号和所述速度信号进行标准化处理；

基于子波的相位类型需求，设定参数，进行常相位子波和最小相位子波的提取或者进行混合相位子波的提取；

进行常相位子波和最小相位子波的提取具体包括以下子步骤：

若所述残差小于等于残差阈值，则输出所述最小相位子波一和所述常相位子波；

在提取出常相位子波和最小相位子波之后，将预处理之后的速度信号转换到加速度域，并与预处理之后的加速度信号进行拼接，获得融合后的数据；

计算融合后的数据的三阶累积量及其三阶谱，利用parzen窗对三阶累积量进行滤波，之后利用双谱计算子波的振幅谱和相位谱，最终，重构成混合相位子波；

进行混合相位子波的提取时，采用下式进行提取：

w＝IFFT(W) (8)

W＝sir(Φ′)|W|+jcos(Φ′)|W| (9)

式中：w为混合相位子波；IFFT表示求取傅里叶逆变换；W为子波的频谱；Φ′为子波的相位谱；|W|为子波的振幅谱；j为虚数。

2.根据权利要求1所述的地震速度信号和加速度信号联合子波提取方法，其特征在于，步骤S1中，基于所述速度信号进行最小相位子波提取时，采用下式进行提取：

a(n)＝u(n)a_o(n) (1)

a_o(n)＝IFFT(A_o(ω)) (3)

A_o(ω)＝ln|X(ω)| (4)

3.根据权利要求1所述的地震速度信号和加速度信号联合子波提取方法，其特征在于，步骤S2中，基于所述加速度信号进行常相位子波提取时，采用下式进行提取：

g(a，n)＝b(n)cos(a)-h_b(n)sin(a) (5)

b(n)＝IFFT(|X(ω)|) (6)

式中：g(a，n)为常相位子波；b(n)为零相位子波；a为相位校正因子；h_b(n)为b(n)的Hilbert变换；IFFT表示求取傅里叶逆变换；|X(ω)|为地震记录的振幅谱。

4.根据权利要求3所述的地震速度信号和加速度信号联合子波提取方法，其特征在于，所述相位校正因子为：

式中：k为正整数。

5.根据权利要求1所述的地震速度信号和加速度信号联合子波提取方法，其特征在于，子波的相位谱Φ′通过下式进行计算：

6.根据权利要求1所述的地震速度信号和加速度信号联合子波提取方法，其特征在于，子波的振幅谱|W|通过下式进行计算：

式中：|W(n)|为子波振幅谱；e为自然底数；W′(n)为子波振幅谱的自然对数；n为时间序列；N为高阶累积量矩阵的大小；|W(N-n)|为子波振幅谱大于二分之一倍子波长度的部分；W′为(N/2)×1的列向量；A_amp为(N²/16)×(N/2)的矩阵；T表示矩阵的转置；-1表示矩阵的逆；X′为(N²/16)×1的列向量。