CN106249285A - 一种有限频率s变换方法及地震信号的时频变换方法 - Google Patents

Abstract

本发明涉及石油地球物理勘探技术领域，公开了一种有限频率S变换方法、一种地震信号的时频变换方法和系统以及一种地震信号的AVAF分析方法，所述有限频率S变换方法包括：基于被测信号确定一个有限频率域；以及基于广义傅里叶变换的α域的表达形式，在所述有限频率域内对被测信号进行S变换。所述地震信号的时频变换方法包括：取地震信号为被测信号；采用所述有限频率S变换方法对地震信号进行时频变换。所述地震信号的AVAF分析方法包括：采用上述时频变换方法对地震信号进行时频变换，得到对应的时频谱；以及基于时频谱，进行地震信号的AVAF分析。本发明既保存了在频率域进行S变换高分辨率的特征，又极大减小了S变换的计算量，为AVAF反演提供了基础。

Description

一种有限频率S变换方法及地震信号的时频变换方法

技术领域

本发明涉及信号时频分析技术领域，特别是涉及应用于石油地球物理勘探技术领域的一种有限频率S变换方法、一种地震信号的时频变换方法和系统以及一种地震信号的AVAF分析方法。

背景技术

在电力系统谐波分析及地震波探测等多个方面，均涉及有对信号进行时频分析。特别是在石油地球物理勘探领域，通常需要针对地震数据进行AVAF反演分析，而进行AVAF反演分析之前，需要采用傅里叶变换(FT)、小波变换(WT)、S变换等时频变换方法对地震信号进行时频分析。

常规的傅里叶变换假设信号或者图像是稳定的，也就是说信号和图像在所有时间、所有位置包含的频率成分是一致的，但这种假设不能适应于分析非平稳信号。短时傅里叶变换(STFT)是其中一种加窗的时频变换方法，然而窗的选择需要在时间域与频率域之间做一个折中，窗越窄信号的时间分辨率越高，窗越宽信号的频率分辨率越高。

小波变换引入了多尺度的概念，使得在对信号进行小波分解时高频具有好的时间分辨率，低频具有比较好的频率分辨率。然而小波变换存变换域不再是频率域，而只是尺度的概念。另外，由于小波变换基的特性导致了相位的概念在小波域中是局部的。这两点均与AVAF分析的需要相矛盾。

S变换(Stockwell Transform)本质上来说是一种广义的傅里叶变换，S变换时窗的选择是与频率有关的，低频部分时窗大，高频部分时窗小，这也是为什么S变换具有多尺度的分辨率的关键。S变换具有非常适合AVF分析的特性，但是S变换的计算量特别大。现有技术提出在频率域的非冗余的S变换的快速算法，它的计算复杂度下降到和快速傅里叶变化相同的级别，但是它在用于AVF分析时存在着阶梯效应。现有技术中还提出了一种用脉冲函数的S变换响应来校正信号的S变换产生的阶梯效应的方法，该方法AVF分析的结果有明显的改善，但是依旧不能满足反演的数据要求。

本领域技术人员所公知的是，AVAF分析需要分频率的地震道集，地震道集的时频变换的精度和效率是反演中必须关注的问题，而快速、精确的时频变换方法严重影响AVAF分析在地震资料中的应用。因此，针对常规的傅里叶变换、小波变换、S变换在地震信号的时频分析中存在的缺点，本发明提出了一种改进的S变换，以实现快速、精确的时频变换。

发明内容

本发明的目的是提供一种有限频率S变换方法、一种地震信号的时频变换方法和系统以及一种地震信号的AVAF分析方法，用于解决目前的时频变换方法精度和效率无法满足反演要求的问题。

为了实现上述目的，本发明的技术方案提供了一种有限频率S变换方法，包括：基于被测信号确定一个有限频率域；以及基于广义傅里叶变换的α域的表达形式，在所述有限频率域内对被测信号进行S变换。

优选地，采用以下公式来表达有限频率域f：

f∈[f_min,f_max]

其中f_min和f_max分别是所述被测信号的最低频率和最高频率。

优选地，在确定的有限频率域中，每个频率的范围固定在所述被测信号的主频段。

优选地，设α域的表达形式为α(v′,v)，则基于广义傅里叶变换的α域的表达形式为：

α(v′,v)＝G(v′+v)·W(v′,σ)

式中，v是频率，v′是局部时间τ的傅里叶变换，G(v′+v)、W(v′,σ)分别是被测信号g(t)和窗函数w(τ-t,σ)的傅里叶谱，τ是局部时间。

优选地，所述在有限频率域内对被测信号进行S变换，具体包括：

从α域对有限频率域范围内的α(v′,v)进行反傅里叶变换，得到S变换的公式：

$S(\mathrm{\τ},v)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\mathrm{\α}(v^{\′},v^{\′})\·e^{i2{\mathrm{\πv}}^{\′}\mathrm{\τ}}{\mathrm{dv}}^{\′}$

式中，S(τ,v)表示被测信号的S变换，v表示频率，v′是局部时间τ的傅里叶变换，τ是局部时间。

本发明的技术方案还提供了一种地震信号的时频变换方法，包括：取地震信号为被测信号；采用上述的有限频率S变换方法对地震信号进行时频变换。

优选地，在确定的有限频率域中，每个频率的范围固定在所述地震信号的主频段。

优选地，设N是地震信号的长度，M是进行时频变换需要分频的个数，则将有限频率域中每个频率的范围固定在所述地震信号的主频段后，采用有限频率S变换方法对地震信号进行时频变换的速度提高N/M倍。

本发明的技术方案还提供了一种地震信号的时频变换系统，包括：有限频率域确定模块，用于基于地震信号确定一个有限频率域；以及时频变换模块，用于采用所述的有限频率S变换方法对地震信号进行时频变换。

本发明的技术方案还提供了一种地震信号的AVAF分析方法，包括：采用上述时频变换方法对地震信号进行时频变换，得到对应的时频谱；以及基于所述时频谱，进行地震信号的AVAF分析。

本发明的有益效果是：本发明相比已有的时频分析的技术，在有限频率域进行S变换，既可以保存在频率域进行S变换高分辨率的特征，又可以极大减小S变换的计算量，为下一步的AVAF反演分析提供了坚实的基础。

本发明的其它特征和优点将在随后的具体实施方式部分予以详细说明。

附图说明

附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明，但并不构成对本发明的限制。在附图中：

图1是本发明的实施一中有限频率S变换方法的流程示意图；

图2是本发明的实施二中地震信号的时频变换方法的流程示意图；

图3是本发明的实施三中地震信号的时频变换系统的结构示意图；

图4是本发明的实施四中地震信号的AVAF分析方法的流程示意图；

图5(a)是本发明的实施方式中用于示例的25Hz雷克子波示意图；

图5(b)是图5(a)的雷克子波的傅里叶谱示意图；

图5(c)是图5(a)的雷克子波的5Hz-100Hz的傅里叶谱示意图；

图5(d)是图5(a)的雷克子波的5Hz-100Hz的窗函数谱示意图；

图5(e)是图5(a)的雷克子波的5Hz-100Hz的有限频率域的示意图；

图5(f)是图5(a)的雷克子波的5Hz-100Hz的S变换的谱的示意图；

图6是本发明的实施方式中有限频率域S变换与非冗余S变换的精度对比图。

具体实施方式

以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。

实施例一

本实施例以信号g(t)为例，其进行S变换正变换的公式如下：

$S(\mathrm{\τ},v)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}g\left(t\right)\frac{\left|v\right|}{\sqrt{2\mathrm{\π}}}{e}^{\frac{{(\mathrm{\τ}-t)}^2{v}^2}{2}}{e}^{-i2\mathrm{\π}vt}dt---\left(1\right)$

其中τ和v分别是S变换之后的时间坐标和频率坐标，这种离散S变换(DST)的计算量非常大，因为它不是一个正交变换，同时该S变换还会得到很多的冗余信息，用公式(1)在时间域内计算长度为N的时间信号的离散S变换需要作N²次积分，它的计算复杂度是O(N³)，需要O(N²)的存储空间。

因此，可知常规的S变换计算都是使用全频带的数据，复杂度过高，会影响进行时频分析的效率的精度。

针对该问题，本实施例提出了一种有限频率S变换方法，如图1所示，包括：基于被测信号确定一个有限频率域；以及基于广义傅里叶变换的α域的表达形式，在所述有限频率域内对被测信号进行S变换。

其中，采用以下公式来表达有限频率域f：

f∈[f_min,f_max]

式中，f_min和f_max分别是所述被测信号的最低频率和最高频率，即表示了本实施例中可根据被测信号的最低频率和最高频率来确定需要的有限频率域。

本实施例中，设α域的表达形式为α(v′,v)，则基于广义傅里叶变换(GFT)的α域的表达形式为：

α(v′,v)＝G(v′+v)·W(v′,σ) (2)

式中，v表示频率，v′是局部时间τ的傅里叶变换，G(v′+v)、W(v′,σ)分别是被测信号g(t)和窗函数w(τ-t,σ)的傅里叶谱，τ是局部时间。其中，所述窗函数w(τ-t,σ)优选为高斯窗函数。

α(v′,v)＝G(v′+v)·W(v′,σ)

式中，v是频率，v′是局部时间τ的傅里叶变换，G(v′+v)、W(v′,σ)分别是被测信号g(t)和窗函数w(τ-t,σ)的傅里叶谱，τ是局部时间。

因此，α域是基于公式(2)来表达的，它是信号的傅里叶谱在平面上展开平移，同时乘以一个高斯窗函数来得到的。

上述公式(2)是结合公式(1)获得的，即是将常规S变换与有限频率域的理论相结合，以确定适用的有限频率域。从公式(1)演化得到公式(2)的具体过程如下的公式(3)所示：

$\begin{array}{c}\mathrm{\α}(v^{\′},v)=\underset{+\mathrm{\∞}}{\overset{-\mathrm{\∞}}{\∫}}S(\mathrm{\τ},v)\·e^{-2{\mathrm{\πv}}^{\′}\mathrm{\τ}}d\mathrm{\τ}=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\{\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}g\left(t\right)\·w(\mathrm{\τ}-t,\mathrm{\σ})\·e^{-2\mathrm{\π}itv}\·dt\}\·e^{-2{\mathrm{\πi\τv}}^{\′}}\·d\mathrm{\τ}\\ =\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}g\left(t\right)\·e^{-2\mathrm{\π}itv}\·dt\·\[W(v^{\′},\mathrm{\σ})\·e^{-2{\mathrm{\πitv}}^{\′}}\]=W(v^{\′},\mathrm{\σ})\·\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}g\left(t\right)\·e^{-2\mathrm{\π}it(v+v^{\′})}\·dt\\ =G(v^{\′}+v)\·W(v^{\′},\mathrm{\σ})\end{array}---\left(3\right)$

公式(3)中α域的表达式实际是S变换的结果在局部时间τ的方向上再进行一次正向的傅里叶变换，这样它就变成另一种二维谱，它的值是信号的傅里叶谱在平面上展开平移，同时乘以一个高斯窗函数的结果。

得到公式(3)后，在确定的有限频率域f∈[f_min,f_max]内对被测信号进行S变换，得到S变换的公式如下：

$S(\mathrm{\τ},v)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\mathrm{\α}(v^{\′},v^{\′})\·e^{i2{\mathrm{\πv}}^{\′}\mathrm{\τ}}{\mathrm{dv}}^{\′}---\left(4\right)$

式中，S(τ,v)表示被测信号的S变换，v是频率，v′是局部时间τ的傅里叶变换，τ是局部时间。

本实施例中，在用公式(4)进行有限频率S变换时，将有限频率域内每个频率的范围固定在被测信号主频段，有利于提高S变换的效率。

本实施例的这种在有限频率域的S变换利用了快速傅里叶变换的高效率，从而使得效率得到大幅提高，在计算长度为N的时间信号时计算复杂度为O(MNlogN)，相比于公式(1)所表示的常规S变换正变换，复杂度显著降低。

实施例二

实施例一中的有限频率S变换方法适用于需要进行时频分析的各个领域，特别是适用于对地震信号进行时频变换。据此，如图2所示，本实施例还给出了一种地震信号的时频变换方法，该方法主要步骤为：取地震信号为被测信号；采用实施例一的有限频率S变换方法对地震信号进行时频变换。

如针对一个地震时间序列x(t)，首先对其进行傅里叶变换得到它的傅里叶谱X(ν)，然后在起傅里叶谱的基础上将它扩展到二维的形式X(v′,ν)，然后再在二维谱的基础上乘以一个高斯窗函数W(v′,σ)，σ表示方差，据此就可以得到地震时间序列α域的表达形式α(v′,v)。最后对α域的数据在v′方向上进行反傅里叶变换，就得到了某一个频率的S变换结果。

在该地震信号的时频变换方法中，在有限频率域内进行S变换时，将频率的范围固定在地震信号的主频段。设N是地震信号的长度，M是进行时频变换需要分频的个数，则将有限频率域中每个频率的范围固定在所述地震信号的主频段后，采用有限频率S变换方法对地震信号进行时频变换的速度提高N/M倍。举例说明，设地震信号的长度是2000，进行时频变换需要分频的个数是100个，则将频率的范围固定在地震信号的主频段后，则S变换的速度可以提高2000/100＝20倍。

实施例三

另外，对应于实施例二的时频变换方法，如图3所示，本实施例还公开了一种地震信号的时频变换系统，包括：有限频率域确定模块，用于基于地震信号确定一个有限频率域；以及时频变换模块，用于采用实施例一中的有限频率S变换方法对地震信号进行时频变换。

该时频变换系统与上述的时频变换方法的具体实施步骤一致，这里不再累述。

实施例四

快速、精确的时频变换会影响AVAF分析在地震资料中的应用，因此本实施例考虑将实施例二中的地震信号的时频变换方法或实施例三的时频变换系统应用于AVAF分析，提出了一种地震信号的AVAF分析方法，如图4所示，包括：采用实施例二的时频变换方法对地震信号进行时频变换，得到对应的时频谱；以及基于所述时频谱，进行地震信号的AVAF分析。当数据量巨大时，在有限频率域的效率还能提高，这取决于反演的频率个数的选择，如果选择在10个以内，按照地震信号的长度为2000个采样点来计算，那么效率将提高200倍。同时在有限频率范围内进行S变换的方法能得到非常平滑的结果，这对反演的数据准备来说是至关重要的，如果利用非冗余S变换得到的谱来进行反演将是不可行的。

需说明的是，在石油地球物理勘探领域，如何根据得到的时频谱进行AVAF分析是已公开的成熟技术，这里不再详细描述。

地震资料的频带通常是几赫兹到一百多赫兹，通常用于AVAF反演只需要几个频率即可，这样既可以保证算法的效率，又能得到令人满意的结果。假如地震信号采样率是2ms，时长为2s，反演频率为10个，那么进行AVAF分析的效率可以提高100倍。

下面以图5所示的正演数据为例，具体介绍将上述各实施例相配合的应用。图5(a)是主频为25HZ的雷克子波(地震子波的一种)，雷克子波在地震数据处理中经常出现，一般认为与地震信号的载波类似，它的时间采样率是1ms。图5(b)是图5(a)的频谱(振幅谱)，Nyquist频率为500hz，该频谱包括了正负频率的子波频谱，且它们是共轭的，其中，求雷克子波的频谱为本领域的成熟技术。图5c是图5(a)在地震波信号在5HZ-100HZ的傅里叶谱G(v′+v)，该傅里叶谱通过将图5(b)中的一维频谱进行频扩展到二维得到。图5(d)是频率域的高斯窗函数谱，示意了不同高斯窗函数，呈现放射状。图5(e)是该图5(a)的地震波信号在α域中的表达形式，是由图5(c)的频移结果G(v′+v)与图5(d)中的高斯窗函数相乘得到的α域的结果。图5(f)是该图5(a)的地震波信号沿着v′的方向进行傅里叶反变换得到的信号的S谱。

如图6所示，通过比较图5(a)至图5(f)有限频率域快速S变换与非冗余S变换的精度发现，可知本实施方式的有限频率S变换得到的谱具有更好的连续性，消除了台阶效应，为下一步的AVAF分析打下坚实的基础。

以上结合附图详细描述了本发明的优选实施例，但是，本发明并不限于上述实施例中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种简单变型，这些简单变型均属于本发明的保护范围。

另外需要说明的是，在上述具体实施例中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合，为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

此外，本发明的各种不同的实施例之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。

Claims (10)

1.一种有限频率S变换方法，其特征在于，包括：

基于被测信号确定一个有限频率域；以及

基于广义傅里叶变换的α域的表达形式，在所述有限频率域内对被测信号进行S变换。

2.根据权利要求1所述的有限频率S变换方法，其特征在于，采用以下公式来表达有限频率域f：

f∈[f_min,f_max]

其中f_min和f_max分别是所述被测信号的最低频率和最高频率。

3.根据权利要求2所述的有限频率S变换方法，其特征在于，在确定的有限频率域中，每个频率的范围固定在所述被测信号的主频段。

4.根据权利要求1至3中任一项所述的有限频率S变换方法，其特征在于，设α域的表达形式为α(v′,v)，则基于广义傅里叶变换的α域的表达形式为：

α(v′,v)＝G(v′+v)·W(v′,σ)

式中，v是频率，v′是局部时间τ的傅里叶变换，G(v′+v)、W(v′,σ)分别是被测信号g(t)和窗函数w(τ-t,σ)的傅里叶谱，τ是局部时间。

5.根据权利要求4所述的有限频率S变换方法，其特征在于，所述在有限频率域内对被测信号进行S变换，具体包括：

从α域对有限频率域范围内的α(v′,v)进行反傅里叶变换，得到S变换的公式：

$S(\mathrm{\τ},v)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\mathrm{\α}(v^{\′},v)\·e^{i2{\mathrm{\πv}}^{\′}\mathrm{\τ}}{\mathrm{dv}}^{\′}$

式中，S(τ,v)表示被测信号的S变换，v表示频率，v′是局部时间τ的傅里叶变换，τ是局部时间。

6.一种地震信号的时频变换方法，其特征在于，包括：

取地震信号为被测信号；

采用权利要求1至5中任一项所述的有限频率S变换方法对地震信号进行时频变换。

7.根据权利要求6所述的时频变换方法，其特征在于，在确定的有限频率域中，每个频率的范围固定在所述地震信号的主频段。

8.根据权利要求7所述的时频变换方法，其特征在于，设N是地震信号的长度，M是进行时频变换需要分频的个数，则将有限频率域中每个频率的范围固定在所述地震信号的主频段后，采用有限频率S变换方法对地震信号进行时频变换的速度提高N/M倍。

9.一种地震信号的时频变换系统，其特征在于，包括：

有限频率域确定模块，用于基于地震信号确定一个有限频率域；以及

时频变换模块，用于采用权利要求1至5中任一项所述的有限频率S变换方法对地震信号进行时频变换。

10.一种地震信号的AVAF分析方法，其特征在于，包括：

采用权利要求6至8中任一项所述时频变换方法对地震信号进行时频变换，得到对应的时频谱；以及

基于所述时频谱，进行地震信号的AVAF分析。