CN116840916B - 一种地震速度信号和加速度信号联合子波提取方法 - Google Patents
一种地震速度信号和加速度信号联合子波提取方法 Download PDFInfo
- Publication number
- CN116840916B CN116840916B CN202310807076.4A CN202310807076A CN116840916B CN 116840916 B CN116840916 B CN 116840916B CN 202310807076 A CN202310807076 A CN 202310807076A CN 116840916 B CN116840916 B CN 116840916B
- Authority
- CN
- China
- Prior art keywords
- wavelet
- phase
- extracting
- signal
- acceleration
- 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
Links
- 230000001133 acceleration Effects 0.000 title claims abstract description 64
- 238000000034 method Methods 0.000 title claims abstract description 19
- 238000000605 extraction Methods 0.000 claims abstract description 49
- 238000012545 processing Methods 0.000 claims abstract description 8
- 238000001228 spectrum Methods 0.000 claims description 71
- 239000011159 matrix material Substances 0.000 claims description 26
- 238000012937 correction Methods 0.000 claims description 15
- 238000001914 filtration Methods 0.000 claims description 2
- 239000000284 extract Substances 0.000 abstract description 3
- 238000004364 calculation method Methods 0.000 description 6
- 230000001186 cumulative effect Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 230000001364 causal effect Effects 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 125000003275 alpha amino acid group Chemical group 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
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. analysis, for interpretation, for correction
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- 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. analysis, for interpretation, for correction
- 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
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/40—Transforming data representation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/624—Reservoir parameters
Abstract
本发明公开了一种地震速度信号和加速度信号联合子波提取方法，包括以下步骤：获取地震的加速度信号和速度信号，并对其进行标准化处理；进行常相位子波和最小相位子波的提取或者进行混合相位子波的提取；常相位子波和最小相位子波的提取具体包括以下子步骤：S1：基于速度信号提取获得最小相位子波一；S2：基于加速度信号提取获得常相位子波，并对其进行转换，获得最小相位子波二；S3：计算两个最小相位子波之间的残差；若残差大于残差阈值，则返回步骤S2，重新进行常相位子波提取；反之则输出最小相位子波一和常相位子波。本发明通过速度信号与加速度信号的一致性联合提取出不同信号类型的子波，实现了地震速度信号和加速度信号的子波联合提取。
Description
技术领域
本发明涉及石油地球物理勘探技术领域，特别涉及一种地震速度信号和加速度信号联合子波提取方法。
背景技术
石油和天然气是国家重要的战略能源，但近年来，随着石油对外依存度的不断提高，石油和天然气的供应和安全问题日益突出，加强技术研究已成为当务之急。反射式地震勘探是利用反射波资料来获得地层结构与岩性的重要手段，而提取子波的准确性严重制约着认识地层结构与岩性的准确度。
目前，进行子波提取时，通常仅采用一种信号进行子波提取，如此对选取信号以及相位校正因子的依赖性大，因此，亟需一种能够提高子波提取准确性的子波提取方法。
发明内容
针对上述问题，本发明旨在提供一种地震速度信号和加速度信号联合子波提取方法。
本发明的技术方案如下：
一种地震速度信号和加速度信号联合子波提取方法，包括以下步骤：
获取地震的加速度信号和速度信号；
对所述加速度信号和所述速度信号进行标准化处理；
基于子波的相位类型需求，设定参数，进行常相位子波和最小相位子波的提取或者进行混合相位子波的提取。
作为优选，进行常相位子波和最小相位子波的提取具体包括以下子步骤：
S1：基于所述速度信号进行最小相位子波提取，获得最小相位子波一；
S2：基于所述加速度信号进行常相位子波提取，获得常相位子波，并对所述常相位子波进行转换，获得最小相位子波二；
S3：计算所述最小相位子波一与所述最小相位子波二之间的残差；
若所述残差大于残差阈值，则返回步骤S2，重新进行常相位子波提取；
若所述残差小于等于残差阈值，则输出所述最小相位子波一和所述常相位子波。
作为优选，步骤S1中，基于所述速度信号进行最小相位子波提取时，采用下式进行提取：
a(n)＝u(n)ao(n) (1)
ao(n)＝IFFT(Ao(ω)) (3)
Ao(ω)＝ln|X(ω)| (4)
式中：a(n)为最小相位子波；u(n)为系数；ao(n)为a(n)的奇部序列；n为时间序列；IFFT表示求取傅里叶逆变换；Ao(ω))为ao(n)的傅里叶变换；|X(ω)|为地震记录的振幅谱。
作为优选，步骤S2中，基于所述加速度信号进行常相位子波提取时，采用下式进行提取：
g(a，n)＝b(n)cos(a)-hb(n)sin(a) (5)
b(n)＝IFFT(|X(ω)|) (6)
式中：g(a，n)为常相位子波；b(n)为零相位子波；a为相位校正因子；hb(n)为b(n)的Hilbert变换。
作为优选，所述相位校正因子为：
式中：k为正整数。
作为优选，进行混合相位子波的提取时，采用下式进行提取：
w＝IFFT(W) (8)
W＝sin(Φ′)|W|+jcos(Φ′)|W| (9)
式中：w为混合相位子波；W为子波的频谱；Φ′为子波的相位谱；|W|为子波的振幅谱；j为虚数。
作为优选，子波的相位谱Φ′通过下式进行计算：
式中：Apha为(N2/4)×N的矩阵；T表示矩阵的转置；-1表示矩阵的逆；为(N2/4)×1的列向量；N为高阶累积量矩阵的大小。
作为优选，子波的振幅谱|W|通过下式进行计算：
式中：|W(n)|为子波振幅谱；e为自然底数；W′(n)为子波振幅谱的自然对数；|W(N-n)|为子波振幅谱大于二分之一倍子波长度的部分；W′为(N/2)×1的列向量；Aamp为(N2/16)×(N/2)的矩阵；X′为(N2/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)ao(n) (1)
ao(n)＝IFFT(Ao(ω)) (3)
Ao(ω)＝ln|X(ω)| (4)
式中：a(n)为最小相位子波；u(n)为系数；ao(n)为a(n)的奇部序列；n为时间序列；IFFT表示求取傅里叶逆变换；Ao(ω))为ao(n)的傅里叶变换；|X(ω)|为地震记录的振幅谱。
上述实施例中，最小相位子波提取公式通过以下步骤推导得到：
假设地震子波为b(n)，子波的对数谱A(ω)可写为：
A(ω)＝lnB(ω) (17)
式中：B(ω)为地震子波的频谱；ln表示求以e为底的对数。
假设a(n)为地震子波b(n)的对数谱时间序列(即最小相位子波)。理论上证明，子波若为最小相位子波，那么其对数谱的时间序列是实的因果序列。任何实序列都可以写为其奇部序列与偶部序列之和，故a(n)可写为：
a(n)＝ao(n)+ae(n) (18)
式中：ae(n)为对数谱时间序列a(n)的偶部序列。
由于a(n)具有因果性，a(n)又可写为式(1)所示的表达式，式(1)中的u(n)如式(2)所示，因此，只需得到奇部序列ao(n)，通过公式(1)就可获得最小相位子波。
对式(31)左右两边取对数，可得：
ln|B(ω)|＝ln|X(ω)| (19)
同时，子波的频谱可以表示为：
式中：为子波的相位。
那么，将式(19)代入式(17)，可得：
同时，有：
Ao(ω)＝AR(ω) (22)
Ae(ω)＝AI(ω) (23)
A(ω)＝Ao(ω)+iAe(ω) (24)
式中：Ao(ω)和Ae(ω)分别ao(n)和ae(n)的傅里叶变换；AR(ω)和AI(ω)分别为子波对数谱A(ω)的实部和虚部，i为虚数单位。
联立式(21)-(24)，可得：
Ao(ω)＝ln|B(ω)| (25)
即：
Ao(ω)＝ln|B(ω)|＝ln|X(ω)| (26)
从而得到式(4)所示的Ao(ω)的计算公式，并由此得到式(3)所示的ao(n)的计算公式。最后，根据n的取值确定u(n)，并通过式(3)计算出奇部序列ao(n)就可通过式(1)提取出基于速度信号的最小相位子波。
S2：基于所述加速度信号进行常相位子波提取，获得常相位子波，并对所述常相位子波进行转换，获得最小相位子波二。
在一个具体的实施例中，基于所述加速度信号进行常相位子波提取时，采用下式进行提取：
g(a,n)＝b(n)cos(a)-hb(n)sin(a) (5)
b(n)＝IFFT(|X(ω)|) (6)
式中：g(a，n)为常相位子波；b(n)为零相位子波；a为相位校正因子；hb(n)为b(n)的Hilbert变换。
可选地，所述相位校正因子为：
式中：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为混合相位子波；W为子波的频谱；Φ′为子波的相位谱；|W|为子波的振幅谱；j为虚数。
可选地，子波的相位谱Φ′通过下式进行计算：
式中：Apha为(N2/4)×N的矩阵；T表示矩阵的转置；-1表示矩阵的逆；为(N2/4)×1的列向量；N为高阶累积量矩阵的大小。
可选地，子波的振幅谱|W|通过下式进行计算：
式中：|W(n)|为子波振幅谱；e为自然底数；W′(n)为子波振幅谱的自然对数；|W(N-n)|为子波振幅谱大于二分之一倍子波长度的部分；W′为(N/2)×1的列向量；Aamp为(N2/16)×(N/2)的矩阵；X′为(N2/16)×1的列向量。
在上述实施例中，利用加速度信号和速度信号的高阶累积量提取混合相位子波，其能够满足子波经过衰减后的特征，提取结果更符合实际。
上述实施例中，混合相位子波提取公式通过以下步骤推导得到：
假设存在一个零均值平稳随机过程y(n)，其中n＝0，±1，±2，...。
则y(n)三阶累计量可表示为：
cum3y(τ1，τ2)＝E[y(n)y(n+τ1)y(n+τ2)] (37)
式中：E[.]表示计算数学期望。
上式可改写为：
k阶累积量的k-1次傅里叶变换即为其的k阶谱，也称为k-1谱。那么，上述三阶累计量的三阶谱(双谱)可表示为：
将式(38)代入式(39)，可得：
即：
式中：Y为y(n)的傅里叶变换。
那么，三阶谱的振幅谱和相位谱可分别表示为：
|cum3y(ω1，ω2)|＝|Y(ω1)||Y(ω2)||Y(ω1+ω2)| (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)的双谱为：
C3x(ω1，ω2)＝C3w(ω1，ω2)C3r(ω1，ω2)+C3v(ω1，ω2) (46)
因为高斯随机变量的高阶累积量为零，则：
C3v(ω1，ω2)＝0 (47)
因此，将公式(47)代入公式(46)，可得：
C3x(ω1，ω2)＝C3w(ω1，ω2)C3r(ω1，ω2) (48)
通常，假设反射系数的高阶累积量为多维冲击函数，因此，可将式(48)写为：
C3x(ω1，ω2)＝C3w(ω1，ω2)/σ2 (49)
式中：σ为反射系数序列的方差。
公式(49)表明子波的三阶谱和地震记录的三阶谱之间只差了一个比例系数。因此，可利用地震记录的三阶谱确定子波的频谱。
对式(42)等式左右两边取对数，可得：
ln|cum3y(ω1，ω2)|＝ln|Y(ω1)|+ln|Y(ω2)|+ln|Y(ω1+ω2)| (50)
将式(50)应用到地震数据中，可得：
ln|cum3x(ω1，ω2)|＝ln|W(ω1)|+ln|W(ω2)|+ln|W(ω1+ω2)| (51)
令ln|cum3x(ω1，ω2)|＝X′(k，l)，ln|W|＝W′，式(51)可写为：
X′(k，l)＝W′(k)+W′(l)+W′(k+l) (52)
式(52)表明X′是W′的线性组合，因此可将式(52)写为矩阵形式：
X′＝AampW′ (53)
其中，X′和Aamp的具体形式分别如式(16)和式(15)所示。
Aamp为列满秩的矩阵，因此，可用广义逆矩阵进行求解公式(53)，解得式(14)所示的W′的表达式，然后利用式(13)即可求解出子波的振幅谱。
将式(13)应用到地震数据中，可得：
式(54)表明是Φw的线性组合，因此可将式(54)写为矩阵形式：
式中：Φ′为N×1的列向量。
式(45)中，和Apha的具体形式分别如式(12)和式(11)所示。
Apha为列满秩的矩阵，因此，可用广义逆矩阵进行求解式(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 (6)
1.一种地震速度信号和加速度信号联合子波提取方法，其特征在于，包括以下步骤：
获取地震的加速度信号和速度信号；
对所述加速度信号和所述速度信号进行标准化处理；
基于子波的相位类型需求，设定参数，进行常相位子波和最小相位子波的提取或者进行混合相位子波的提取；
进行常相位子波和最小相位子波的提取具体包括以下子步骤：
S1：基于所述速度信号进行最小相位子波提取，获得最小相位子波一；
S2：基于所述加速度信号进行常相位子波提取，获得常相位子波，并对所述常相位子波进行转换，获得最小相位子波二；
S3：计算所述最小相位子波一与所述最小相位子波二之间的残差；
若所述残差大于残差阈值，则返回步骤S2，重新进行常相位子波提取；
若所述残差小于等于残差阈值，则输出所述最小相位子波一和所述常相位子波；
在提取出常相位子波和最小相位子波之后，将预处理之后的速度信号转换到加速度域，并与预处理之后的加速度信号进行拼接，获得融合后的数据；
计算融合后的数据的三阶累积量及其三阶谱，利用parzen窗对三阶累积量进行滤波，之后利用双谱计算子波的振幅谱和相位谱，最终，重构成混合相位子波；
进行混合相位子波的提取时，采用下式进行提取：
w＝IFFT(W) (8)
W＝sir(Φ′)|W|+jcos(Φ′)|W| (9)
式中：w为混合相位子波；IFFT表示求取傅里叶逆变换；W为子波的频谱；Φ′为子波的相位谱；|W|为子波的振幅谱；j为虚数。
2.根据权利要求1所述的地震速度信号和加速度信号联合子波提取方法，其特征在于，步骤S1中，基于所述速度信号进行最小相位子波提取时，采用下式进行提取：
a(n)＝u(n)ao(n) (1)
ao(n)＝IFFT(Ao(ω)) (3)
Ao(ω)＝ln|X(ω)| (4)
式中：a(n)为最小相位子波；u(n)为系数；ao(n)为a(n)的奇部序列；n为时间序列；IFFT表示求取傅里叶逆变换；Ao(ω))为ao(n)的傅里叶变换；|X(ω)|为地震记录的振幅谱。
3.根据权利要求1所述的地震速度信号和加速度信号联合子波提取方法，其特征在于，步骤S2中，基于所述加速度信号进行常相位子波提取时，采用下式进行提取：
g(a，n)＝b(n)cos(a)-hb(n)sin(a) (5)
b(n)＝IFFT(|X(ω)|) (6)
式中：g(a，n)为常相位子波；b(n)为零相位子波；a为相位校正因子；hb(n)为b(n)的Hilbert变换；IFFT表示求取傅里叶逆变换；|X(ω)|为地震记录的振幅谱。
4.根据权利要求3所述的地震速度信号和加速度信号联合子波提取方法，其特征在于，所述相位校正因子为：
式中：k为正整数。
5.根据权利要求1所述的地震速度信号和加速度信号联合子波提取方法，其特征在于，子波的相位谱Φ′通过下式进行计算：
式中：Apha为(N2/4)×N的矩阵；T表示矩阵的转置；-1表示矩阵的逆；为(N2/4)×1的列向量；N为高阶累积量矩阵的大小。
6.根据权利要求1所述的地震速度信号和加速度信号联合子波提取方法，其特征在于，子波的振幅谱|W|通过下式进行计算：
式中：|W(n)|为子波振幅谱；e为自然底数；W′(n)为子波振幅谱的自然对数；n为时间序列；N为高阶累积量矩阵的大小；|W(N-n)|为子波振幅谱大于二分之一倍子波长度的部分；W′为(N/2)×1的列向量；Aamp为(N2/16)×(N/2)的矩阵；T表示矩阵的转置；-1表示矩阵的逆；X′为(N2/16)×1的列向量。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310807076.4A CN116840916B (zh) | 2023-07-04 | 2023-07-04 | 一种地震速度信号和加速度信号联合子波提取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310807076.4A CN116840916B (zh) | 2023-07-04 | 2023-07-04 | 一种地震速度信号和加速度信号联合子波提取方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116840916A CN116840916A (zh) | 2023-10-03 |
| CN116840916B true CN116840916B (zh) | 2024-03-26 |
Family
ID=88170211
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310807076.4A Active CN116840916B (zh) | 2023-07-04 | 2023-07-04 | 一种地震速度信号和加速度信号联合子波提取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116840916B (zh) |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102707314A (zh) * | 2012-05-28 | 2012-10-03 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 一种多路径双谱域混合相位子波反褶积方法 |
| CN103018775A (zh) * | 2012-11-15 | 2013-04-03 | 中国石油天然气股份有限公司 | 基于相位分解的混合相位子波反褶积方法 |
| WO2013152468A1 (zh) * | 2012-04-13 | 2013-10-17 | 中国石油天然气集团公司 | 一种地层品质因子反演方法 |
| CN105510968A (zh) * | 2015-12-31 | 2016-04-20 | 中国海洋大学 | 一种基于地震海洋学的海水物性测量方法 |
| CN111708083A (zh) * | 2020-06-05 | 2020-09-25 | 成都理工大学 | 一种基于模型的深度域地震子波提取方法 |
| CN113156493A (zh) * | 2021-05-06 | 2021-07-23 | 中国矿业大学 | 一种使用归一化震源的时频域全波形反演方法及装置 |
| US11500117B1 (en) * | 2021-09-13 | 2022-11-15 | Institute Of Geology And Geophysics, Chinese Academy Of Sciences | Method and system for evaluating filling characteristics of deep paleokarst reservoir through well-to-seismic integration |
| CN115453629A (zh) * | 2022-08-10 | 2022-12-09 | 中国石油大学(北京) | 一种基于时频域自适应滤波的波形校正方法及装置 |
| CN115980845A (zh) * | 2023-02-21 | 2023-04-18 | 中海石油(中国)有限公司 | 基于激发相位先验约束的子波提取方法、系统、设备 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110988986B (zh) * | 2019-12-25 | 2021-01-01 | 成都理工大学 | 改善深层碳酸盐岩储层刻画精度的地震资料低频增强方法 |
| CN111505708B (zh) * | 2020-04-28 | 2021-04-20 | 西安交通大学 | 一种基于深度学习的强反射层剥离方法 |
-
2023
- 2023-07-04 CN CN202310807076.4A patent/CN116840916B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2013152468A1 (zh) * | 2012-04-13 | 2013-10-17 | 中国石油天然气集团公司 | 一种地层品质因子反演方法 |
| CN102707314A (zh) * | 2012-05-28 | 2012-10-03 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 一种多路径双谱域混合相位子波反褶积方法 |
| CN103018775A (zh) * | 2012-11-15 | 2013-04-03 | 中国石油天然气股份有限公司 | 基于相位分解的混合相位子波反褶积方法 |
| CN105510968A (zh) * | 2015-12-31 | 2016-04-20 | 中国海洋大学 | 一种基于地震海洋学的海水物性测量方法 |
| CN111708083A (zh) * | 2020-06-05 | 2020-09-25 | 成都理工大学 | 一种基于模型的深度域地震子波提取方法 |
| CN113156493A (zh) * | 2021-05-06 | 2021-07-23 | 中国矿业大学 | 一种使用归一化震源的时频域全波形反演方法及装置 |
| US11500117B1 (en) * | 2021-09-13 | 2022-11-15 | Institute Of Geology And Geophysics, Chinese Academy Of Sciences | Method and system for evaluating filling characteristics of deep paleokarst reservoir through well-to-seismic integration |
| CN115453629A (zh) * | 2022-08-10 | 2022-12-09 | 中国石油大学(北京) | 一种基于时频域自适应滤波的波形校正方法及装置 |
| CN115980845A (zh) * | 2023-02-21 | 2023-04-18 | 中海石油(中国)有限公司 | 基于激发相位先验约束的子波提取方法、系统、设备 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116840916A (zh) | 2023-10-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Hale et al. | A fixed-point continuation method for l1-regularized minimization with applications to compressed sensing | |
| CN106405645B (zh) | 一种基于资料品质分析的信噪比可控的地震拓频处理方法 | |
| CN107179550B (zh) | 一种数据驱动的地震信号零相位反褶积方法 | |
| CN105427264A (zh) | 一种基于群稀疏系数估计的图像重构方法 | |
| CN110646841B (zh) | 时变稀疏反褶积方法及系统 | |
| Blu et al. | Self-similarity: Part II—Optimal estimation of fractal processes | |
| CN114200525B (zh) | 一种自适应的多道奇异谱分析地震数据去噪方法 | |
| CN112835103A (zh) | 自适应去鬼波与宽频准零相位反褶积联合处理方法及系统 | |
| CN104730576A (zh) | 基于Curvelet变换的地震信号去噪方法 | |
| Li | A principal component analysis approach to noise removal for speech denoising | |
| CN116840916B (zh) | 一种地震速度信号和加速度信号联合子波提取方法 | |
| CN112526599A (zh) | 基于加权l1范数稀疏准则的子波相位估计方法及系统 | |
| Shatilo | Seismic Phase Unwrapping: Methods, Results, PROBLEMS1 | |
| CN104422956A (zh) | 一种基于稀疏脉冲反演的高精度地震谱分解方法 | |
| CN113109873B (zh) | 一种基于秩残差约束的沙漠地区地震信号噪声抑制方法 | |
| CN113310684B (zh) | 基于尺度空间和改进稀疏表示的齿轮箱故障特征提取方法 | |
| Hao et al. | Denoising method based on spectral subtraction in time-frequency domain | |
| CN109991657B (zh) | 基于逆二分递推奇异值分解的地震资料高分辨率处理方法 | |
| Huang et al. | Laser gyro signal filtering by combining CEEMDAN and principal component analysis | |
| CN109085649B (zh) | 一种基于小波变换优化的地震资料去噪方法 | |
| Wang et al. | A self-supervised method using noise2noise strategy for denoising crp gathers | |
| CN117289343A (zh) | 基于小波倒谱的叠前道集提频方法、设备及介质 | |
| CN116719085A (zh) | 一种地震记录高分辨率处理方法、装置、设备及存储介质 | |
| CN113419275B (zh) | 一种基于稀疏字典学习的高分辨率地震处理方法 | |
| CN113341463B (zh) | 一种叠前地震资料非平稳盲反褶积方法及相关组件 |
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 |