CN102116868B - 一种地震波分解方法 - Google Patents

Abstract

一种地震波分解的方法，该方法包括选择以Gabor小波或Morlet小波作为子波形式；对获取的地震信号进行多子波的分解；在每一次分解，通过地震信号的复数属性计算小波的三个参数：时移u、中心频率ω、相移φ；根据这三个参数计算时间尺度参数σ；估计小波的振幅a；通过迭代确定最终参数及分解出不同主频的子波。本发明提供的地震子波分解方法，对传统的匹配追踪方法进行了改进，提高了计算效率；同时σ是可自适应可调节的。

Description

一种地震波分解方法

技术领域

本发明涉及地球物理勘探领域，更特别地，涉及石油地球物理勘探领域的一种地震波分解方法。

背景技术

地震勘探技术是利用弹性振动研究地层中不同岩石的弹性差异来认识地层的分布情况和地质构造状况的一种地球物理勘探方法，主要应用于石油、天然气勘探开发、矿产资源普查、工程地质勘查等领域，也是石油勘探中寻找储油圈闭最普遍最有效的一种途径。由地震所提供的地下信息不仅可以有效地划分圈闭范围，而且还可以研究判断地层的岩石性质和流体性质，甚至还可以扩展到油藏特征描述等方面。目前，存在多种对地震信号进行处理的方法，其中地震信号多子波分解方法可以将采集到的地震信号（地震道数据）分解为各种不同的子波，这些子波具有不同主频、不同子波宽度以及不同时间位置。将这些子波叠加在一起即可精确恢复原有的地震信号以用于地震信号分析。

近年来匹配追踪技术应用到地震信号处理中，将一个地震波分解成一系列子波，这些子波在子波时频分解中被称为时频单位。在匹配追踪过程中，为了合理地表征连续小波，需要利用以下五个参数：振幅A、时移u、尺度（时间轴长度）σ、平均频率ω_m和相位φ。常规的匹配追踪算法是一个很麻烦的迭代过程，需要从大量丰富的小波字典中选取出合适参数的子波。

对于连续小波而言，重要参数之一是尺度σ，它控制了频域带宽和时域子波长度。但是基于窗式傅式变换或者Gabor变换时进行谱分解却忽略了尺度这个可调参数，甚至于在一些匹配追踪技术应用中也被忽略了，所以现有的子波长度是不变的。比如基于Liu&Marfurt思想的则将小波宽度当作时中心频率的函数，在这种情况下实际就是常规小波变换。因此现有的采用匹配追踪技术的地震波分解方法存在问题：一、由于大量翻阅小波字典，导致计算速度很慢；二、尺度σ为固定值，子波为固定长度。

发明内容

为解决上述问题，本发明提供一种地震波分解方法，该方法包括：

1）选择子波形式，该子波形式为由公式（1）定义的Gabor小波形式：

$g_r\left(t\right)=w\left(\frac{t-u}{\mathrm{\σ}}\right)\exp [i\left(\mathrm{\ω}(t-u)\right)+\mathrm{\φ}]$ 公式（1），

其中g_r(t)是Gabor小波函数、t是时间变量、

是高斯窗函数、u是时移、σ是时间轴长度、ω是中心频率、φ是相移、i表示虚数单位；

2）获取地震信号，并对获取的地震信号进行多次分解，直至满足分解完成条件，分解后的地震信号的形式为公式（2）所定义的形式：

$f\left(t\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{a}_n{g}_{r_n}\left(t\right)+R^{(N+1)}f$ 公式（2），

其中f(t)是地震信号的函数，n为分解次数，a_n是与第n次分解的小波

对应的振幅值，N为正整数，R^(N+1)f为第N次分解后的残差，且R⁽¹⁾f＝f，标记j为1至N的正整数；

根据所述地震信号的复数属性来计算第j次分解时的参数集γ_j＝{u_j,σ_j,ω_j,φ_j}；

在计算出所述参数集γ_j＝{u_j,σ_j,ω_j,φ_j}的基础上计算第j次分解的小波的振幅值a_j；

3）将所述计算出的γ_j＝{u_j,σ_j,ω_j,φ_j}代入到公式（1）以得到在第j次分解出的小波将计算出的振幅值a_j与该小波

相乘以得到第j次分解出的子波。

根据本发明的可替换方法，其中公式（1）可以由以下公式来替代：

$m\left(t\right)=\exp [-\left(\frac{\ln 2}{{\mathrm{\π}}^2}\right)\frac{{\mathrm{\ω}}_m^2{(t-u)}^2}{{\mathrm{\σ}}^2}]\exp \left[i({\mathrm{\ω}}_m(t-u)+\mathrm{\φ})\right]$    公式（6），

其中m(t)是Morlet小波函数、u是时移、σ是时间轴长度、ω_m是平均角频率、φ是相移。

本发明提供的地震波分解方法，对传统的匹配追踪方法进行了改进，其尺度参数σ是可自适应调节的，通过σ的变化可以控制子波的长度。在对分解后得到的多子波进行分析中，还可以通过限制σ的大小来抑制脉冲或者正弦信号等干扰信号。

附图说明

图1是根据本发明的一个实施方式的一种地震波分解的方法的流程图；

图2是根据本发明的另一个实施方式的一种地震波分解的方法的流程图。

具体实施方式

下面结合附图对本发明提供的地震波分解方法做进一步描述。

如图1所示，根据本发明的一个实施方式，提供一种地震波分解方法，包括：

1）选择子波形式，该子波形式为由公式（1）定义的Gabor小波形式：

$g_r\left(t\right)=w\left(\frac{t-u}{\mathrm{\σ}}\right)\exp [i\left(\mathrm{\ω}(t-u)\right)+\mathrm{\φ}]$    公式（1），

其中g_r(t)是Gabor小波函数，

表示高斯窗函数、u是时移、σ是时间轴长度、ω是中心频率、φ是相移、i表示虚数单位；选择子波形式可以是从本领域技术人员公知的Gabor小波字典、离散Dirac函数和基于复指数的离散傅里叶函数中选择匹配地震信号的成一系列小波或小波原子，优选地本发明提供的地震波分解方法选择由公式（1）定义的Gabor小波形式作为子波形式。Gabor小波是一种时域的复函数，波形类似于复正弦信号与高斯信号包络相乘的结果。w(t)可以例如为

2）获取地震信号，并对获取的地震信号进行多次分解，直至满足分解完成条件，分解后的地震信号的形式为公式（2）所定义的形式：

$f\left(t\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{a}_n{g}_{r_n}\left(t\right)+R^{(N+1)}f$    公式（2），

其中f(t)是地震信号的函数，n为分解次数，a_n是与第n次分解的小波

对应的振幅值，N为正整数，R^(N+1)f为第N次分解后的残差，且R⁽¹⁾f＝f，标记j为1-N的正整数；对于每一次分解，都有 $R^{\left(j\right)}f=a_j{g}_{r_j}+R^{(j+1)}f,$ 其中

与R^(j+1)f正交。

根据所述地震信号的复数属性来计算第j次分解时的参数集γ_j＝{u_j,σ_j,ω_j,φ_j}，其中计算参数集γ_j＝{u_j,σ_j,ω_j,φ_j}包括：对初始信号进行Hilbert变换以构建复数信号，将该复数信号的最大振幅包络处所对应的时间作为时移u_j，并将该最大振幅包络处对应的时刻处的瞬时频率和瞬时相位分别作为中心频率ω_j和相移

根据公式（3）估算参数σ_j：

$g_{r_j}\left(t\right)=\arg \underset{g_{r_j}\&Element;D}{\max }\frac{|<R^{\left(j\right)}f,g_{r_j}>|}{\left|\right|g_{r_j}\left|\right|}$    公式（3），

其中D是小波集合，表示R^(j)f和

的内积，为第j次分解的小波

的标准化；其中，如果j=1，则所述初始信号为所获取的地震信号；如果j>1，则所述初始信号为在第j-1次分解后的残差R^(j)f。

具体地，在根据Hilbert变换得到参数u_j,ω_j,φ_j后，给定参数σ_j的范围，在该给定范围中选取多个值作为σ_j分别与参数u_j,ω_j,φ_j一起代入到公式（1）中，由此得到一组小波集合，设该小波集合为D。也就是说小波集合D由相同值的参数u_j,ω_j,φ_j和一组不同值的参数σ_j确定的一组小波。然后根据公式（3）在小波集合D中寻找满足公式（3）的小波

该小波中的σ_j即为估算出的σ_j。这样就避免了传统方法中需要翻阅小波字典。所述给定范围可以是例如1-20。在给定范围中多个值的选取可以根据估算精度的需求而有所不同，例如，在给定范围中选择多个值，相邻两个值之间的间隔可以是0.5-2，优选为1。

在计算出参数集γ_j＝{u_j,σ_j,ω_j,φ_j}的基础上计算第j次分解时的小波的振幅值a_j，计算振幅值a_j包括根据公式（4）来估算a_j：

$a_j=\frac{|<R^{\left(j\right)}f,g_{r_j}>|}{{\left|\right|g_{r_j}\left|\right|}^2}$    公式（4），

其中 $R^{\left(j\right)}f=a_j{g}_{r_j}+R^{(j+1)}f,$

且公式（4）满足条件 $\{a_j,{\mathrm{\&gamma;}}_j\}=\arg \underset{a_j,{\mathrm{\&gamma;}}_j\&Element;\mathrm{\&Gamma;}}{\min }{\left|\right|R^{(j+1)}f\left(t\right)\left|\right|}^2$    （公式5）；

其中Γ表示复数域。

这里，将计算出的参数集γ_j＝{u_j,σ_j,ω_j,φ_j}代入到公式（1）得到

然后根据以上公式（4）可以估算出a_j。

3）将所述计算出的γ_j＝{u_j,σ_j,ω_j,φ_j}代入到公式（1）以得到在第j次分解出的小波

将计算出的振幅值a_j与该小波

相乘以得到第j次分解出的子波；

其中所述分解完成条件可以为所述残差R^(N+1)f的能量小于所述地震信号平均能量的2%，优选为小于1%。地震信号的能量的计算方法可以采用本领域技术人员公知的方法，这里不再赘述。

对于当前分解而言，

其中

与R^(j+1)f是正交函数，满足 ${\left|\right|R^{(j+1)}f\left|\right|}^2={|<R^{\left(j\right)}f,g_{r_j}>|}^2/{\left|\right|g_{r_j}\left|\right|}^2+{\left|\right|R^{(j+1)}f\left|\right|}^2.$ 对于正交投影（公式3），选取使

最大，即R^(j+1)f²最小时的

因此，公式（4）满足残差最小的条件：

$\{a_j,{\mathrm{\&gamma;}}_j\}=\arg \underset{a_j,{\mathrm{\&gamma;}}_j\&Element;\mathrm{\&Gamma;}}{\min }{\left|\right|R^{(j+1)}f\left(t\right)\left|\right|}^2$    公式（5）。

本发明提供的地震波分解方法是多次迭代分解的过程。对于第一次迭代分解，有

即

此时将f(t)作为初始信号并根据上述方法确定γ₁＝{u₁,σ₁,ω₁,φ₁}和a₁，由此得到第一次迭代分解时分解出的子波接下来，由于

则有

在第二次迭代分解，将R⁽²⁾f作为初始信号，根据上述方法确定γ₂＝{u₂,σ₂,ω₂,φ₂}和a₂，由此得到第二次迭代分解时分解出的子波

在形式上有

在第三次迭代分解，将R⁽³⁾f作为初始信号，并根据上述方法确定γ₃＝{u₃,σ₃,ω₃,φ₃}和a₃。依次类推，在经过N次迭代分解后就能得到N个子波和残差R^(N+1)f，即如公式（2）描述的形式。当残差满足一定的条件后，不再对残差进行分解，该地震波分解结束。残差满足的条件根据处理地震信号的类型、环境、分解精度的不同而有所不同，例如待分解的地震信号的类型、特性以及所需的分解精度等这些因素发生变化，残差满足的条件也可以发生变化。残差满足的条件可以是如上所述残差的能量小于所述地震信号平均能量的2%，优选地，可以是1%。在分解结束后，残差R^(N+1)f可以被当作噪音数据对其进行处理，例如将其过滤掉。

可替换地，可以采用Morlet小波作为时频分解的子波单位，因为Morlet小波能够较好地表征弹性波在孔隙介质中传播所产生的能量衰减和速度频散。因此可以选择Morlet小波形式作为子波形式，该Morlet小波由公式（6）定义：

$m\left(t\right)=\exp [-\left(\frac{\ln 2}{{\mathrm{\&pi;}}^2}\right)\frac{{\mathrm{\&omega;}}_m^2{(t-u)}^2}{{\mathrm{\&sigma;}}^2}]\exp \left[i({\mathrm{\&omega;}}_m(t-u)+\mathrm{\&phi;})\right]$    公式（6），

其中m(t)是Morlet小波函数、t是时间变量、u是时移、σ是时间轴长度、ω_m是平均角频率、φ是相移、i表示虚数单位。采用Morlet小波对地震信号进行分解的方法与采用Gabor小波的相同，根据上述的方法计算参数集γ＝{u,σ,ω_m,φ}和a_n，即计算参数集γ＝{u,σ,ω_m,φ}的方法与计算参数集γ＝{u,σ,ω,φ}的方法相同，具体为：

1）选择子波形式，该子波形式为由（6）定义的Morlet小波形式：

$m\left(t\right)=\exp [-\left(\frac{\ln 2}{{\mathrm{\&pi;}}^2}\right)\frac{{\mathrm{\&omega;}}_m^2{(t-u)}^2}{{\mathrm{\&sigma;}}^2}]\exp \left[i({\mathrm{\&omega;}}_m(t-u)+\mathrm{\&phi;})\right]$    公式（6），

其中m(t)是Morlet小波函数、t是时间变量、u是时移、σ是时间轴长度、ω_m是平均角频率、φ是相移、i表示虚数单位；

2）获取地震信号，并对获取的地震信号进行多次分解，直至满足分解完成条件，分解后的地震信号的形式为公式（7）所定义的形式：

$f\left(t\right)=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_n{m}_n\left(t\right)+R^{(N+1)}f$    公式（7），

其中f(t)是地震信号的函数，n为分解次数，a_n是与第n次分解的小波m_n对应的振幅值，N为正整数，R^(N+1)f为第N次分解后的残差，R^(j)f＝a_jm_j+R^(j+1)f且R⁽⁰⁾f＝f，标记j为1至N的正整数；

根据所述地震信号的复数属性来计算第j次分解时的参数集其中计算第j次分解时的参数集

包括：对初始信号进行Hilbert变换以构建复数信号，将该复数信号的最大振幅包络处所对应的时间作为时移u_j，并将该最大振幅包络处对应的时刻处的瞬时频率和瞬时相位分别作为平均角频率

和相移

根据公式（8）估算参数σ_j：

$m_j\left(t\right)=\arg \underset{m_j\&Element;D}{\max }\frac{|<R^{\left(j\right)}f,m_j>|}{\left|\right|m_j\left|\right|}$    公式（8），

其中D是小波集合，R^(j)f,m_j表示R^(j)f和m_j的内积，

为第j次分解的小波m_j的标准化；

其中，如果j=1，则所述初始信号为所获取的地震信号；如果j>1，则所述初始信号为在第j-1次分解后的残差R^(j)f。

与Gabor小波的情况类似，在根据Hilbert变换得到参数

后，给定参数σ_j的范围，在该给定范围中选取多个值作为σ_j分别与参数

一起代入到公式（6）中，由此得到一组小波集合，设该小波集合为D。也就是说小波集合D由相同值的参数

和一组不同值的参数σ_j确定的一组小波。然后根据公式（8）在小波集合D中寻找满足公式（8）的小波m_j，该小波m_j中的σ_j即为估算出的σ_j。所述给定范围可以例如是1-20。在给定范围中多个值的选取可以根据估算精度的需求而有所不同，例如，相邻两个值之间的间隔可以是0.5-2，优选为1。

在计算出所述参数集

的基础上计算第j次分解时的小波的振幅值a_j；计算第j次分解的小波的振幅值a_j包括根据公式（9）来计算该振幅值a_j：

$a_j=\frac{|<R^{\left(j\right)}f,m_j>|}{{\left|\right|m_j\left|\right|}^2}$    公式（9），

其中R⁽ⁱ⁾f＝a_im_i+R⁽ⁱ⁺¹⁾f，

且公式（9）满足条件 $\{a_j,{\mathrm{\&gamma;}}_j\}=\arg \underset{a_j,{\mathrm{\&gamma;}}_j\&Element;\mathrm{\&Gamma;}}{\min }{\left|\right|R^{(j+1)}f\left(t\right)\left|\right|}^2$    （公式10）；

这里，将计算出的参数集代入到公式（6）得到m_j，然后根据以上公式（9）可以估算出a_j。

3）将所述计算出的

代入到公式（6）以得到在第j次分解出的小波m_j，将计算出的振幅值a_j与该小波m_j相乘以得到第j次分解出的子波。

与采用Gabor小波的情形类似，如果所述残差R^(j)f满足一定条件，则分解完成。残差满足的条件可以是残差的能量小于所述地震信号平均能量的2%，优选为1%。

另外，采用Morlet小波将获得的原始地震信号（地震道数据）分解为多个子波，这样能够很明显地看出能量在时频域的分布。

在该方法中除了可输出子波重构的数据外还可以输出子波的时频谱，进而进行时频分析：

将地震信号f(t)分解为一系列小波

n＝0,1,2,...,N-1，这样能够很明显地看出能量在时频域的分布。在时频域信号f(t)的能量密度表达式为：

$\mathrm{Ef}(t,\mathrm{\&omega;})=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\frac{a_n^2}{{\left|\right|g_{r_n}\left|\right|}^2}{\mathrm{Wg}}_{r_n}(t,\mathrm{\&omega;})$    公式（11）

其中

是某一选定子波的Wigner分布。同样地，我们将时频域振幅谱定义为：

$\mathrm{Af}(t,\mathrm{\&omega;})=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\frac{a_n}{g_{r_n}}\sqrt{{\mathrm{Wg}}_{r_n}(t,\mathrm{\&omega;})}$    公式（12）

对于公式（6）所给出的Morlet小波m(t)，其时频域的Wigner分布可定义为：

${\mathrm{Wg}}_r(t,\mathrm{\&omega;})=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{+\&infin;}m(t+\frac{\mathrm{\&tau;}}{2})\stackrel{\&OverBar;}{m}(t-\frac{\mathrm{\&tau;}}{2})\exp [-\mathrm{i\&omega;\&tau;}]\mathrm{d\&tau;}$    公式（13）

其中为m(t)的复共轭。再次利用积分变换式 ${\&Integral;}_{-\&infin;}^{+\&infin;}\exp [-a{\mathrm{\&tau;}}^2-\mathrm{b\&tau;}-c]\mathrm{d\&tau;}=\sqrt{\frac{\mathrm{\&pi;}}{a}}\exp [\frac{b^2}{4a}-c],$ 可得：

${\mathrm{Wg}}_r(t,\mathrm{\&omega;})=\sqrt{\frac{\mathrm{\&pi;}}{2\ln 2}}\frac{\mathrm{\&sigma;}}{{\mathrm{\&omega;}}_m}\exp [-\left(\frac{{\mathrm{\&pi;}}^2}{2\ln 2}\right)\frac{{\mathrm{\&sigma;}}^2{(\mathrm{\&omega;}-{\mathrm{\&omega;}}_m)}^2}{{\mathrm{\&omega;}}_m^2}]\&times;\exp [-\left(\frac{2\ln 2}{{\mathrm{\&pi;}}^2}\right)\frac{{\mathrm{\&omega;}}_m^2{(t-u)}^2}{{\mathrm{\&sigma;}}^2}]$    公式（14）

因此其时频域振幅谱解析式可写为：

$\mathrm{Af}(t,\mathrm{\&omega;})=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\frac{a_n}{\left|\right|g_{r_n}\left|\right|}{\left(\sqrt{\frac{\mathrm{\&pi;}}{2\ln 2}}\frac{{\mathrm{\&sigma;}}_n}{{\mathrm{\&omega;}}_n}\right)}^{1/2}\&times;\exp [-\left(\frac{{\mathrm{\&pi;}}^2}{4\ln 2}\right)\frac{{\mathrm{\&sigma;}}_n^2{(\mathrm{\&omega;}-{\mathrm{\&omega;}}_n)}^2}{{\mathrm{\&omega;}}_n^2}]\&times;\exp [-\left(\frac{\ln 2}{{\mathrm{\&pi;}}^2}\right)\frac{{\mathrm{\&omega;}}_n^2{(t-u_n)}^2}{{\mathrm{\&sigma;}}_n^2}]$

                      公式（15）

其中

为第n个小波的平均频率，

为归一化因数。

为了提高对地震信号进行分解时的分解速度，一个方面是要提高模值

或

的计算速度，因为在上述计算参数σ_j和a_j的公式（3）（公式（8））和公式（4）（公式（9））中都需要对或（||m_j||²或||m_j||）进行计算。如果能提高

或

（||m_j||²或||m_j||）的计算速度，那么对地震波分解的速度会有很大提高。因为地震信号是实值函数，因此公式（1）可以改写成或替换成公式（16）：

$g_r\left(t\right)=w\left(\frac{t-u}{\mathrm{\&sigma;}}\right)\cos [\mathrm{\&omega;}(t-u)+\mathrm{\&phi;}]$    公式（16）

其中，g_r(t)是Gabor小波函数，是高斯窗函数、u是时移、σ是时间轴长度、ω是中心频率、φ是相移；

同样地公式（6）可以改写成或替换成公式（17）：

$m\left(t\right)=\exp [-\left(\frac{\ln 2}{{\mathrm{\&pi;}}^2}\right)\frac{{\mathrm{\&omega;}}_m^2{(t-u)}^2}{{\mathrm{\&sigma;}}^2}]\cos [{\mathrm{\&omega;}}_m(t-u)+\mathrm{\&phi;}]$    公式（17）

其中m(t)是Morlet小波函数、u是时移、σ是时间轴长度、ω_m是平均角频率、φ是相移。

通过下面的公式（18）和公式（19）

$D_0={\&Integral;}_{-\&infin;}^{+\&infin;}{\left[w\left(\frac{t-u}{\mathrm{\&sigma;}}\right)\right]}^2\mathrm{dt}$    公式（18）

$D_1\left(\mathrm{\&omega;}\right)={\&Integral;}_{-\&infin;}^{+\&infin;}{\left[w\left(\frac{t-u}{\mathrm{\&pi;}}\right)\right]}^2\exp \left[i2(\mathrm{\&omega;}(t-u)+\mathrm{\&phi;})\right]\mathrm{dt}$    公式（19）

可以计算：

${\left|\right|g_{r_j}\left|\right|}^2=\frac{1}{2}(D_0+\mathrm{Re}[D_1\left(\mathrm{\&omega;}\right)\left]\right)$    公式（20）。

根据公式（18）和公式（19）计算：

$D_0=\frac{\mathrm{\&sigma;\&pi;}}{{\mathrm{\&omega;}}_j}\sqrt{\frac{\mathrm{\&pi;}}{2\ln 2}};$

$D_1\left({\mathrm{\&omega;}}_j\right)=\frac{\mathrm{\&sigma;\&pi;}}{{\mathrm{\&omega;}}_j}\sqrt{\frac{\mathrm{\&pi;}}{2\ln 2}}\exp [-\frac{{\mathrm{\&pi;}}^2{\mathrm{\&sigma;}}^2}{2\ln 2}+i2\mathrm{\&phi;}],$ 由此可以得出

${\left|\right|g_{r_j}\left|\right|}^2=\frac{\mathrm{\&pi;}}{2}\sqrt{\frac{\mathrm{\&pi;}}{2\ln 2}}\frac{{\mathrm{\&sigma;}}_j}{{\mathrm{\&omega;}}_j}(1+\exp [-\frac{{\mathrm{\&pi;}}^2{{\mathrm{\&sigma;}}_j}^2}{2\ln 2}\left]\cos 2{\mathrm{\&phi;}}_j\right)$    公式（21）

以及类似地

${\left|\right|m_j\left|\right|}^2=\frac{\mathrm{\&pi;}}{2}\sqrt{\frac{\mathrm{\&pi;}}{2\ln 2}}\frac{{\mathrm{\&sigma;}}_j}{{\mathrm{\&omega;}}_{m_j}}(1+\exp [-\frac{{\mathrm{\&pi;}}^2{{\mathrm{\&sigma;}}_j}^2}{2\ln 2}\left]\cos 2{\mathrm{\&phi;}}_j\right)$    公式（22）

在Gabor小波的情况中，将上述计算得到的参数集γ_j＝{u_j,σ_j,ω_j,φ_j}代入到公式（1）和公式（21）；将根据公式（1）得到的

（计算的振幅能量）和根据公式（21）得到的利用公式（3）以重新计算σ_j；将根据公式（1）得到的

和根据公式（21）得到的

利用公式（4）以重新计算a_j，这一步就是匹配；用重新计算得到的σ_j和a_j分别替换之前估算的σ_j和a_j。根据Hilbert变换得到的参数u_j,ω_j,φ_j以及这里重新计算得到的参数σ_j和a_j即为进行匹配后的参数，利用参数u_j,ω_j,φ_j以及这里重新计算得到的参数σ_j和a_j进行更精确的追踪方法以重新计算γ_j＝{u_j,σ_j,ω_j,φ_j}和a_j。所述追踪方法可以是本领域技术人员公知的匹配追踪算法中的追踪方法。例如，如果对γ_j＝{u_j,σ_j,ω_j,φ_j}中的u_j进行追踪，则可以先保持其它三个参数σ_j,ω_j,φ_j的值不变，在给定范围（根据半时长表达式计算宽度，其中l_t为半时长，ε为最小振幅值，单位为dB，可以设定ε＝-120dB，该宽度即为给定范围）内进行扫描，扫描的间隔例如可以是1，在扫描范围中计算地震信号与gabor小波的内积，寻找其能量与小波模平方的比值最大的点，此时对应的u_j就是最优的u_j。追踪其它三个参数σ_j,ω_j,φ_j的方法与追踪u_j的类似。在进行追踪得到γ_j＝{u_j,σ_j,ω_j,φ_j}后，可以根据该γ_j＝{u_j,σ_j,ω_j,φ_j}和公式（4）得到追踪后的a_j。最后利用追踪后的γ_j＝{u_j,σ_j,ω_j,φ_j}和a_j可以得到分解出的子波。

与Gabor子波的情况类似，在Morlet小波的情况中，将上述计算得到的参数集

代入到公式（6）和公式（22）；将根据公式（6）得到的m_j（计算的振幅能量）和根据公式（22）得到的||m_j||利用公式（8）以重新计算σ_j；将根据公式（6）得到的m_j和根据公式（22）得到的||m_j||²利用公式（9）以重新计算a_j，这一步就是匹配；用重新计算得到的σ_j和a_j分别替换之前估算的σ_j和a_j。根据Hilbert变换得到的参数以及这里重新计算得到的参数σ_j和a_j即为进行匹配后的参数，再利用根据Hilbert变换得到的参数

以及这里重新计算得到的参数σ_j和a_j进行更精确的追踪方法以重新计算

和a_j。这里的追踪方法可以是本领域技术人员公知的匹配追踪算法中的追踪方法。例如，如果对

中的u_j进行追踪，则可以先保持其它三个参数

的值不变，在给定范围（例如，根据半时长表达式计算宽度，其中l_t为半时长，ε为最小振幅值，单位为dB，可以设定ε＝-120dB，该宽度即为给定范围）内进行扫描，扫描的间隔例如可以是1，在扫描范围中计算地震信号与morlet小波m_j的内积，寻找其能量与小波模平方||m_j||²的比值最大的点，此时对应的u_j就是最优的u_j。追踪其它三个参数的方法与追踪u_j的类似，也是先保持参数集

中的三个参数的值固定，对剩下的一个参数在给定范围内进行扫描，找到满足公式（8）的最优的该参数值。在进行追踪方法得到

后，可以根据该

和公式（9）得到追踪后的a_j。最后利用追踪后的

和a_j可以得到分解出的子波。

在追踪过程中也应用了公式（21）或（22）进行模值的计算。

这样分别采用公式（21）和（22）在匹配和追踪过程中都可以提高公式（3）和公式（4）以及公式（8）和公式（9）的计算速度。

在传统方法中由于没有小波模的解析表达式，需要大量查阅小波字典，通过多次的匹配来完成，所以运算效率较低，而且尺度参数σ为一固定值使得子波长度是固定的。使用本发明提供的地震波分解方法，应用小波模及时频谱的解析表达式进行计算，不需要大量查阅小波字典，从而提高分解速度；而且尺度参数σ可以调节，从而子波的长度是可以变化的；不仅能减少分解后的残差，还能丰富小波字典。

在对地震信号进行分解之后，可以对分解出的各个子波进行处理，例如滤波、去噪等。对这些子波进行叠加就能重构地震信号。但是由于地震波分解时的频率间隔不可能无限的小，这就会使得一些信号缺失，由此使得重构的地震信号反射同相轴不连续。为此，本发明提供的一种地震波分解的方法还包括对多个子波进行重构，在地震波重构时进行反褶积处理以保证同相轴的连续性。所采用的反褶积处理可以是本领域技术人员公知的反褶积方法，优选为多道预测反褶积方法。

另外，当尺度参数σ的值极小时，对应的子波为脉冲信号；当σ值很大时，对应的子波为正弦信号。在重构时可以通过限定尺度参数σ的范围来剔除尖脉冲、单频谐振等噪声。例如可以将σ＜0.4以及σ＞10所对应的子波过滤掉，利用σ滤波，可以使重构后的地震信号的时频谱更为清晰，能更有效地进行分析。

本发明提供的方法可以通过硬件、软件或者硬件与软件的结合来实现，这些硬件为所属领域技术人员公知的器件，例如检波器、滤波器、处理器等；所述软件也为所属领域技术人员公知的载有采用计算机语言编程的程序的介质，例如软盘、硬盘、光盘以及移动硬盘等；采用的计算机语言为所属领域公知的语言。

Claims (7)

1.一种地震波分解方法，该方法包括：

1）选择子波形式，该子波形式为由公式（6）定义的Morlet小波形式：

$m\left(t\right)=\exp [-\left(\frac{\ln 2}{{\mathrm{\&pi;}}^2}\right)\frac{{\mathrm{\&omega;}}_m^2{(t-u)}^2}{{\mathrm{\&sigma;}}^2}]\exp \left[i({\mathrm{\&omega;}}_m(t-u)+\mathrm{\&phi;})\right]$ 公式（6），

其中m(t)是Morlet小波函数、t是时间变量、u是时移、σ是时间轴长度、ω_m是平均角频率、φ是相移、i表示虚数单位；

2）获取地震信号，并对获取的地震信号进行多次分解，直至满足分解完成条件，分解后的地震信号的形式为公式（7）所定义的形式：

$f\left(t\right)=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{a}_n{m}_n\left(t\right)+R^{(N+1)}f$ 公式（7），

其中f(t)是地震信号的函数，n为分解次数，a_n是与第n次分解的小波m_n对应的振幅值，N为正整数，R^(N+1)f为第N次分解后的残差，R^(j)f＝a_jm_j+R^(j+1)f且R⁽⁰⁾f＝f，标记j为1至N的正整数，m_j是第j次分解的小波；

根据所述地震信号的复数属性来计算第j次分解时的参数集 ${\mathrm{\&gamma;}}_j=\{u_j,{\mathrm{\&sigma;}}_j,{\mathrm{\&omega;}}_{m_j},{\mathrm{\&phi;}}_j\};$

在计算出所述参数集

的基础上计算第j次分解时的小波的振幅值a_j；

3）将所述计算出的

代入到公式（6）以得到在第j次分解出的小波m_j，将计算出的振幅值a_j与该小波m_j相乘以得到第j次分解出的子波；

其中计算第j次分解时的参数集

包括：

对初始信号进行Hilbert变换以构建复数信号，将该复数信号的最大振幅包络处所对应的时间作为时移u_j，并将该最大振幅包络处对应的时刻处的瞬时频率和瞬时相位分别作为平均角频率

和相移根据公式（8）估算参数σ_j：

$m_j\left(t\right)=\arg \underset{m_j\&Element;D}{\max }\frac{|<R^{\left(j\right)}f,m_j>|}{\left|\right|m_j\left|\right|}$ 公式（8），

其中D是小波集合，<R^(j)f,m_j>表示R^(j)f和m_j的内积，

为第j次分解的小波m_j的标准化；

其中，如果j=1，则所述初始信号为所获取的地震信号；如果j>1，则所述初始信号为在第j-1次分解后的残差R^(j)f；

其中根据公式（8）估算参数σ_j包括：

给定参数σ_j的范围，在该范围中选取多个值；将该多个值分别与所述参数u_j、

和

一起代入到公式（6）以得到多个小波，该多个小波即为小波集合D；根据公式（8）在该小波集合D中确定满足公式（8）的一个小波m_j，该小波m_j对应的σ_j为所估算的σ_j；

其中计算第j次分解的小波的振幅值a_j包括：根据公式（9）来计算该振幅值a_j：

$a_j=\frac{|<R^{\left(j\right)}f,m_j>|}{{\left|\right|m_j\left|\right|}^2}$ 公式（9），

其中公式（9）满足条件 $\{a_j,{\mathrm{\&gamma;}}_j\}=\arg \underset{a_j,{\mathrm{\&gamma;}}_j\&Element;\mathrm{\&Gamma;}}{\min }{\left|\right|R^{(j+1)}f\left(t\right)\left|\right|}^2$ 公式（10），

Γ表示复数域；

其中该方法还包括根据公式（22）来计算||m_j||²或||m_j||：

${\left|\right|m_j\left|\right|}^2=\frac{\mathrm{\&pi;}}{2}\sqrt{\frac{\mathrm{\&pi;}}{2\ln 2}}\frac{{\mathrm{\&sigma;}}_j}{{\mathrm{\&omega;}}_{m_j}}(1+\exp [-\frac{{\mathrm{\&pi;}}^2{{\mathrm{\&sigma;}}_j}^2}{2\ln 2}\left]\cos 2{\mathrm{\&phi;}}_j\right)$ 公式（22）；

其中根据公式（22）来计算||m_j||²或||m_j||包括：

将对初始信号进行Hilbert变换以构建复数信号后获得的参数

以及根据公式（8）估算的参数σ_j代入到公式（22）以得到||m_j||²；

计算||m_j||²的平方根以得到||m_j||。

2.根据权利要求1所述的地震波分解方法，其中所述参数σ_j的范围为1-20，所述多个值中相邻两个值之间的间隔为1。

3.根据权利要求1所述的地震波分解方法，其中所述分解完成条件为所述残差R^(N+1)f的能量小于所述地震信号平均能量的2%。

4.根据权利要求1所述的地震波分解方法，其中估算参数σ_j还包括：将根据公式（22）计算得到的||m_j||代入到公式（8），并根据公式（8）重新估算参数σ_j；用该重新估算的参数σ_j替换之前估算的参数σ_j。

5.根据权利要求4所述的地震波分解方法，其中计算振幅值a_j还包括：将根据公式（22）计算得到的||m_j||²代入到公式（9），根据公式（9）重新计算振幅值a_j；用该重新计算的振幅值a_j替换之前计算的振幅值a_j。

6.根据权利要求5所述的地震波分解方法，其中计算第j次分解时的参数集

还包括：根据公式（8）使用追踪方法对参数集

进行追踪以得到追踪后的参数集

以及

计算第j次分解时的小波的振幅值a_j还包括：使用追踪后的参数集

根据公式（9）计算追踪后的振幅值a_j。

7.根据权利要求6所述的地震波分解方法，其中在对参数集

进行追踪的过程中利用公式（22）计算公式（8）中的||m_j||；以及

在计算追踪后的振幅值a_j的过程中利用公式（22）计算公式（9）中的||m_j||²。

Images