CN104422956B - 一种基于稀疏脉冲反演的高精度地震谱分解方法 - Google Patents

Abstract

本发明提供了一种基于稀疏脉冲反演的高精度地震谱分解方法，属于油气地球物理勘探领域。本发明方法包括：(1)提取地震道上的m个变子波，所述m为时间样点个数；(2)构建所述变子波的褶积模型的矩阵形式；(3)求取所述褶积模型矩阵形式中的反射系数的稀疏形式，即稀疏脉冲；(4)利用所有变子波与其对应的稀疏脉冲构成新的变子波序列；(5)计算所述新的变子波序列的时频谱。利用本发明方法得到的信号的时频谱是单个子波频谱的叠加，因此具有很高的时频聚焦性。

Description

一种基于稀疏脉冲反演的高精度地震谱分解方法

技术领域

本发明属于油气地球物理勘探领域，具体涉及一种基于稀疏脉冲反演的高精度地震谱分解方法。

背景技术

地震信号是典型的非平稳信号即频率随着时间变化而变化。地下介质对不同频率成分地震波的反射特征是不一样的，最常见的比如地震波在地下介质中的吸收衰减，高频成分衰减的快，低频衰减得慢，此时地震波能量补偿主要是高频成分的补偿，再比如近年来岩石物理实验室发现的含流体的岩体对低频成分有增强的作用等等。这些都是地震信号非平稳的原因，如何对地震信号时频分解呢？复数道地震信号是比较早的信号时频分析思想，随后短时傅立叶变换开启了地震信号时频分析的大门，目前地震信号分析中常用的时频分析方法如小波变换，S变换，广义S变换、以及衍生出的改进的小波变换，S变换等。随着油气勘探的深入高精度的地震信号时频分析方法呼之欲出，近年来以匹配追踪为代表的谱分解方法受到地球物理工作者的关注，其基本思想是将地震信号分解为一系列子波的叠加之和，分别对子波做谱分解然后叠加各个子波的时频谱。

现有的地震信号时频分析方法中，都是先给定一个母波，将母波进行伸缩变换成一些列的子波，将地震信号分解成这些子波的加权之和，每个子波的加权系数反映地震信号的频率成分的多少。上述的子波是人为给定的子波与实际的地震子波相差甚远。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种基于稀疏脉冲反演的高精度地震谱分解方法，提高时频聚焦性。

本发明是通过以下技术方案实现的：

一种基于稀疏脉冲反演的高精度地震谱分解方法，包括：

(1)提取地震道上的m个变子波，所述m为时间样点个数；

(2)构建所述变子波的褶积模型的矩阵形式；

(3)求取所述褶积模型矩阵形式中的反射系数的稀疏形式，即稀疏脉冲；

(4)利用所有变子波与其对应的稀疏脉冲构成新的变子波序列；

(5)计算所述新的变子波序列的时频谱。

所述步骤(1)是这样实现的：

设时窗内地震记录的平滑振幅谱为则该时窗内的变子波为：

其中ifft为反傅立叶变换算子。T_i为地震道上第i个时窗，时窗的长度是m个时间采样点的长度。

所述步骤(2)是这样实现的：

变子波的褶积模型的矩阵如公式(2)所示：

其中s(t)为一个地震道上的地震信号，t₁，t₂....t_m为一个地震道上地震信号的时间离散点，w(t′)为利用步骤(1)求得的变子波，下标对应地震道信号的离散时间样点，r(t)为反射系数，ε为噪音。

所述步骤(3)中是利用相关向量机方法反演得到所述稀疏脉冲的。

所述步骤(4)是这样实现的：

在地震信号的时间t_k处新的变子波为：

其中，为时间t_k处新的变子波，r(t_k)为通过步骤(3)得到的稀疏脉冲，为通过步骤(1)得到的变子波。

所述步骤(5)是这样实现的：

对通过步骤(4)得到的新的变子波做winger-ville分布得到该新的变子波的时频谱，然后对于下一个时间样点t_k+1重复上述步骤(4)和步骤(5)，直至所有时间样点的计算完毕；

然后叠加所有的时间样点处的新的变子波的时频谱即得到了新的变子波序列的时频谱(即该道地震道上的地震信号的时频谱)。

与现有技术相比，本发明的有益效果是：本发明方法是基于地震波的子波而非常规方法的人为给定一个小波，而且利用本发明方法得到的信号的时频谱是单个子波频谱的叠加，因此具有很高的时频聚焦性。

附图说明

图1a是地震信号。

图1b是地震信号的小波变换。

图1c是利用本发明方法得到的地震信号时频谱。

图2是本发明方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

本发明提取时变空变子波(以下都简称变子波)，这个变子波可以近似看成是地震的实际子波，将每一道地震信号分解为变子波与瞬时稀疏脉冲的褶积之和，对每个时变子波与瞬时稀疏脉冲卷积构成的新子波，并对新子波做谱分解构成信号的瞬时谱。

(1)变子波提取

地震子波经过地下介质的时候会发生变化，比如振幅、相位的变化、不同频率成分的吸收衰减等等，因此子波是随着深度是变化的，本发明首先在很小的时窗内(比如这个小时窗为15个时间采样点)取一段地震数据，此时数据的长度为15，由于数据样点太少本发明在数据的末尾充0，将截取的这段地震数据扩充为100。截取的这段地震数据的傅立叶谱可以看成是子波谱与反射系数谱的乘积，由于子波谱相对比较平滑，反射系数谱“杂乱”，因此对时窗内地震记录谱平滑后可以估计出子波的振幅谱。地震子波的振幅谱是很好获得的，难以准确获得的是相位谱，本发明中感兴趣的只是子波的振幅谱，子波的相位对最终时频分析结果没有影响，因此这里将子波的相位认为0，得到一个对称的子波。具体的做法是：设某时窗内地震记录的平滑振幅谱为则时窗内的子波为：

其中ifft为反傅立叶变换算子。T_i为地震道上第i个时窗。

(2)稀疏脉冲反演

设反射系数为r(t)，子波为w(t′)，地震记录s(t)＝w(t′)*r(t)，其中的*表示卷积运算。本发明推导出的瞬时子波褶积模型的矩阵形式为：

其中t₁，t₂....t_m为其道地震信号的时间离散点，r(t)为反射系数时间序列，ε为噪音。

w₁(t′)，w₂(t′).....w_m(t′)为瞬时子波，下标对应地震道信号的离散时间样点。这些瞬时子波是通过上面(1)中的方法得到的。

而经典褶积公式中子波w(t′)是不随时间变换的，也就是公式(2)中的w₁(t′)，w₂(t′)....w_m(t′)都为w(t′)，即不随时间变化。

方程(2)简记为：

S_k＝Y_k(w，r)+ε_k (3)

k＝1，2，.....m，为时间样点序列号，Y为子w波与r的褶积运算符(也称为卷积)，ε为随机噪音。假设ε服从于均值为0，方差为σ²的高斯分布，于是：

求解r和σ理论上可以通过(4)式的最大似然函数求得，但这样求得的反射系数大部分都是0，导致过学习，假设反射系数服从于均值为0，方差为1/h的超参数分布即：

p(r_k|h_k)＝N(r_k|0，1/h_k) (5)

通过相关向量机可以求得r的稀疏形式(请参考Michael E Tipping.TheRelevance Vector Machine.In Sara A Solla，Todd K Leen，and Klaus-Robert Muller，editors，Advances in Neural Information Processing Systems12.Cambridge，Mass：MIT Press，2000.，此文献给出了求解方程(3)中反射系数r的方法，本发明是应用这个文献中的方法来解方程的)。

(3)基于稀疏脉冲反演的时频分布

通过上述的发明内容(2)可以得到稀疏脉冲，设在地震信号的时间t_k处新的子波记为：

对新的子波做winger-ville分布，对于下一个时间样点t_k+1重复上述新子波重构以及winger-ville分布计算，所有时间样点计算完毕后，叠加所有的winger-ville分布就得到了一道地震信号时频谱。叠加是这样实现的：新子波都是有时间范围的，比如新子波1的时间范围为5ms-15ms，新子波2的时间范围为8ms-18ms，这样8ms-15ms是他们重合的时间范围，分别对新子波1，新子波2做winger-ville分布后，8ms-15ms范围内的时频谱数值就是新子波1新子波2时频谱(时频谱就是做了winger-ville分布)在8-15ms范围内数值之和。

(4)如图2所示，本发明方法实现的步骤如下：

①按照上述(1)中的方法提取m个变子波m为时间样点个数(时窗的长度就是m个时间采样点的长度)。

②按照上述公式(2)构建出变子波的矩阵褶积模型。

③依照上述(2)中的方法求解出稀疏的反射系数。

④按照上述(3)中的方法计算时频谱。

图1a为地震信号，图1b为小波变换，图1c为本发明的谱分解。对比后两个图可以看出本发明的谱分解时频聚焦性更好，不论是时间分辨率还是频率分辨率。

本发明直接从地震信号中提取时变的子波而非常规的人为给定一个子波，推导出变子波矩阵褶积模型，并利用相关向量机的方法计算出稀疏的地层反射系数，将时变子波与地震反射时间上的稀疏反射系数相乘构成新的子波，对新子波做高精度的winger-ville分布，然后叠加所有时间样点上的winger-ville谱。本发明的时频谱具有很高的时频聚焦性。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (4)

1.一种基于稀疏脉冲反演的高精度地震谱分解方法，其特征在于：所述方法包括：

(1)提取地震道上的m个变子波，所述m为时间样点个数；

(2)构建所述变子波的褶积模型的矩阵形式；

(3)求取所述褶积模型矩阵形式中的反射系数的稀疏形式，即稀疏脉冲；

(4)利用所有变子波与其对应的稀疏脉冲构成新的变子波序列；

(5)计算所述新的变子波序列的时频谱,

所述步骤(4)是这样实现的：

在地震信号的时间t_k处新的变子波为：

$w_{New,t_k}=w_{t_k}.r\left(t_k\right)---\left(6\right)$

其中，为时间t_k处新的变子波，r(t_k)为通过步骤(3)得到的稀疏脉冲，为通过步骤(1)得到的变子波,

所述步骤(5)是这样实现的：

对通过步骤(4)得到的新的变子波做winger—ville分布得到该新的变子波的时频谱，然后对于下一个时间样点t_k+1重复上述步骤(4)和步骤(5)，直至所有时间样点的计算完毕；

然后叠加所有的时间样点处的新的变子波的时频谱即得到了新的变子波序列的时频谱。

2.根据权利要求1所述的基于稀疏脉冲反演的高精度地震谱分解方法，其特征在于：所述步骤(1)是这样实现的：

设时窗内地震记录的平滑振幅谱为w_T(f)，则该时窗内的变子波为：

$w_{T_i}\left(t\right)=ifft\left(w_T\right(f\left)\right)---\left(1\right)$

其中ifft为反傅立叶变换算子，T_i为地震道上第i个时窗，时窗的长度是m个时间采样点的长度。

3.根据权利要求2所述的基于稀疏脉冲反演的高精度地震谱分解方法，其特征在于：所述步骤(2)是这样实现的：

变子波的褶积模型的矩阵如公式(2)所示：

$\left[\begin{array}{c}s\left(t_1\right)\\ s\left(t_2\right)\\ ...\\ s\left(t_m\right)\end{array}\right]=\left[\begin{array}{cccc}{w}_1\left({t^{\′}}_n\right)& & & \\ w_1\left({t^{\′}}_{n+1}\right)& w_2\left({t^{\′}}_n\right)& & \mathrm{....}\\ \mathrm{...}& w_2\left({t^{\′}}_{n+1}\right)& \mathrm{......}& w_m\left({t^{\′}}_{N-1}\right)\\ & \mathrm{...}& & w_m\left({t^{\′}}_N\right)\end{array}\right]\left[\begin{array}{c}r\left(t_1\right)\\ r\left(t_2\right)\\ \mathrm{....}\\ r\left(t_m\right)\end{array}\right]+\left[\begin{array}{c}{\mathrm{\ϵ}}_1\\ {\mathrm{\ϵ}}_2\\ \mathrm{....}\\ {\mathrm{\ϵ}}_m\end{array}\right]---\left(2\right)$

其中s(t)为一个地震道上的地震信号，t₁,t₂....t_m为一个地震道上地震信号的时间离散点，w(t')为利用步骤(1)求得的变子波，w(t')的下标对应地震道信号的时间离散点的序数，r(t)为反射系数，ε为噪音。

4.根据权利要求3所述的基于稀疏脉冲反演的高精度地震谱分解方法，其特征在于：所述步骤(3)中是利用相关向量机方法反演得到所述稀疏脉冲的。