CN103048684B - 一种多分量地震资料面波压制方法 - Google Patents

Abstract

本发明公开了一种多分量地震资料面波压制方法，该方法包括：获取在多个时间采样点采样得到的单个多分量矢量地震道信号S_k(t)；对S_k(t)进行时频双域分解，以得到时频谱信号W_gS_k(t，a)；对时频谱信号W_gS_k(t，a)进行瞬时极化分析，其中该瞬时极化分析包括计算出针对各个时间采样点的瞬时极化椭圆率P(t，a)；基于视速度大小，将W_gS_k(t，a)分为非面波区信号和面波区信号，并对面波区信号进行面波信号压制处理，以得到压制面波信号后的时频谱信号W′_gS_k(t，a)；对W′_gS_k(t，a)进行反变换，并重构以得到压制面波后的矢量地震道信号S′_k(t)；其中，面波信号压制处理包括：统计出在各个频率下的非面波区信号的均方根能量E_nonGR(a)；针对每个频率下的面波区信号，以该频率下的均方根能量E_nonGR(a)为能量阈值将面波区信号分为强能量信号和弱能量信号后进行分类约束极化滤波。

Description

一种多分量地震资料面波压制方法

技术领域

本发明涉及地震资料数字信号处理领域，具体涉及多分量地震信号时频域矢量面波压制方法。

背景技术

在多分量地震资料处理领域中，通常需要将强面波噪声与有效体波信号进行分离。其有效信号中弱转换波信号通常被强面波所淹没，而且转换波频带与面波频带有很大部分重叠，这使得压制强面波的同时保护弱转换波的问题变得非常棘手。利用小波变换(或者S变换)对多分量地震信号进行多尺度分解，可以在时频域将重叠信号分开；将分开后的信号根据面波强能量、大椭圆极化率的矢量极化特性进行面波识别，并从分解后的信号中减去识别出的面波信号，然后进行重构便可以得到面波压制后的信号。

面波压制所涉及的主要数字信号处理技术有：小波变换(S变换)、希尔伯特变换、瞬时极化分析。小波变换是现代谱分析的强有力工具，它既能考察信号局部时域的频变特性，又能考察局部频域的时变特性。对于频率域重叠的非平稳信号采用小波变换进行多尺度分解分析同样适用于地震资料信号处理。通过希尔伯特变换，将原信号变换为解析信号后，可以很方便地应用瞬时极化分析算法进行瞬时极化率的计算，从而可以有效地将多分量矢量处理技术结合到该项面波压制技术中。据资料记载，采用双域极化滤波进行面波压制最为典型的代表为Diallo等人(2003年，2005年)，Kulesh(2006年，2007年)提出的小波域自适应瞬时极化分析面波压制技术。但是该方法在极大压制面波的同时难以保护与该面波频带重叠的有效信号不被压制，从而在面波压制和保护有效信号这两个方面难以兼顾。

发明内容

为解决上述问题，本发明的目的是提供一种能够压制多波地震资料中强面波的同时又能保护与其频带重叠的有效信号的面波压制方法。

为实现上述目的，本发明提供了一种多分量地震资料面波压制方法，该方法包括：

获取在多个时间采样点采样得到的单个多分量矢量地震道信号S_k(t)，k为正整数；

对所述S_k(t)进行时频双域分解，以得到时频谱信号W_gS_k(t，a)，t表示时间，a表示频率；

对所述时频谱信号W_gS_k(t，a)进行瞬时极化分析，其中该瞬时极化分析包括计算出针对各个时间采样点的瞬时极化椭圆率P(t，a)；

基于视速度大小，将所述W_gS_k(t，a)分为非面波区信号和面波区信号，并对所述面波区信号进行面波信号压制处理，以得到压制面波信号后的时频谱信号W′_gS_k(t，a)；

对所述W′_gS_k(t，a)进行反变换，并重构以得到压制面波后的矢量地震道信号S′_k(t)；

其中，所述面波信号压制处理包括：

统计出在各个频率下的所述非面波区信号的均方根能量E_nonGR(a)；

针对每个频率下的所述面波区信号，以该频率下的所述均方根能量E_nonGR(a)为能量阈值将所述面波区信号分为强能量信号和弱能量信号后进行分类约束极化滤波。

通过上述技术方案，可以有效地压制多分量地震资料中的强面波，同时保护低频弱转换信号不被压制，从而极大提高多波多分量地震资料的信噪比。

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

附图说明

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

图1是根据本发明的一个实施方式的多分量地震资料面波压制方法的流程图；

图2示出了通过使用本发明的实施方式的方法的单个矢量地震道面波压制效果图；

图3示出了通过使用本发明的实施方式的方法的四川新场地区多分量资料面波压制效果图；以及

图4示出了使用本发明的实施方式提供的方法与常用商业软件压制效果对比图。

具体实施方式

图1是根据本发明的一个实施方式的多分量地震资料面波压制方法的流程图。如图1所示，根据本发明的一个实施方式，提供了一种多分量地震资料面波压制方法，该方法可以包括：

获取在多个时间采样点采样得到的单个多分量矢量地震道信号S_k(t)，k为正整数；

对所述S_k(t)进行时频双域分解，以得到时频谱信号W_gS_k(t，a)，t表示时间，a表示频率(即，尺度因子)；

对所述时频谱信号W_gS_k(t，a)进行瞬时极化分析，其中该瞬时极化分析包括计算出针对各个时间采样点的瞬时极化椭圆率P(t，a)；

基于视速度大小，将所述W_gS_k(t，a)分为非面波区信号和面波区信号，并对所述面波区信号进行面波信号压制处理，以得到压制面波信号后的时频谱信号W′_gS_k(t，a)；

对所述W′_gS_k(t，a)进行反变换，并重构以得到压制面波后的矢量地震道信号S′_k(t)；

其中，所述面波信号压制处理包括：

统计出在各个频率下的所述非面波区信号的均方根能量E_nonGR(a)；

针对每个频率下的所述面波区信号，以该频率下的所述均方根能量E_nonGR(a)为能量阈值将所述面波区信号分为强能量信号和弱能量信号后进行分类约束极化滤波。

其中，所述强能量信号可以采用由等式(1)描述的滤波器进行滤波：

等式(1)，其中F_strong(t，a)表示

滤波后的强能量信号，P_strong表示强能量(面波)信号所在的极化椭圆率分布范围；以及

所述弱能力信号可以采用由等式(2)描述的滤波器进行滤波：

等式(2)，其中F_weak(t，a)表示滤波后的弱能量信号，P_weak表示弱能量(面波)信号所在的极化椭圆率分布范围。

在本发明的优选实施方式中，所述P_stonrg的取值可以为P_strong∈(0，0.25]，所述P_waek的取值可以为P_weak∈[0.45，1)。

在本发明的一个实施方式中，在频率a下的所述均方根能量E_nonGR(a)可以根据等式(3)来计算：

$E_{\mathrm{nonGR}}\left(a\right)=\sqrt{1/M\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\left|W_g{S}_k(\mathrm{i\Δt},a)\right|}^2}$ 等式(3)

其中M为所述非面波区信号的时间采样点个数，Δt为采样时间间隔。

对S_k(t)可以采用本领域技术人员公知的小波变换或S变换来分解得到W_gS_k(t，a)。以小波变换为例，可以利用复数母小波函数g(t)的连续小波变换根据等式(4)来分解所述S_k(t)：

$W_g{S}_k(t,a)=<g_{t,a},S_k>=\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}\frac{1}{a}{g}^{*}\left(\frac{\mathrm{\&tau;}-t}{a}\right)S_k\left(\mathrm{\&tau;}\right)\mathrm{d\&tau;},$ a∈R，t∈R    等式(4)

其中，R表示实数域，(·)^*表示复共轭，所采用的复数母小波(例如复数Morlet母小波)函数

这里，W_gS_k(t，a)的虚部是实部的希尔伯特变换，所以该信号W_gS_k(t，a)实际就是一个解析信号，可以直接进行瞬时极化分析。但是如果采用实数母小波进行小波变换或者S变换得到W_gS_k(t，a)，则随后还需采用希尔伯特变换将时频谱信号构造成解析时频谱信号。

在本发明中，计算瞬时极化椭圆率的方法可以包括例如协方差矩阵极化分析法、单点瞬时极化分析法等。例如，所述进行瞬时极化分析可以包括：

根据等式(5)采用希尔伯特变换构造解析时频谱

$W_g{S}_k^C(t,a)=W_g{S}_k(t,a)+\mathrm{iH}\left(W_g{S}_k(t,a)\right)$ 等式(5)，H表示希尔伯特变换；这里，如果是采用复小波变换得到所述W_gS_k(t，a)，则该步骤可以省略。

将所述改写为 表示复数信号的模，并求取使得值最大时的相位函数φ₀(t，a)， ${\mathrm{\&phi;}}_0(t,a)=\frac{1}{2}\arg [B(t,a)+\mathrm{\&epsiv;C}(t,a)]+\mathrm{\&pi;n},$ n∈N    N为正整数，其中 $B(t,a)=\frac{1}{2}\underset{k}{\mathrm{\&Sigma;}}{\left(W_g{S}_k^C(t,a)\right)}^2,$ $C(t,a)=\frac{1}{2}{\left(\underset{k}{\mathrm{\&Sigma;}}{W}_g{S}_k^C(t,a)\right)}^2,$ ε＜＜1，n为周期信号的周期数。

根据所述相位函数求取极化椭圆的主轴函数R(t，a)和次轴函数r(t，a)，其中， $R(t,a)=\mathrm{Real}[e^{-i{\mathrm{\&phi;}}_0(t,a)}\&CenterDot;W_g{S}_k^C(t,a)],$ $r(t,a)=\mathrm{Real}[e^{-i({\mathrm{\&phi;}}_0(t,a)+\mathrm{\&pi;}/2)}\&CenterDot;W_g{S}_k^C(t,a)],$ Real表示取复数信号的实部；

根据等式(6)计算所述极化椭圆率P(t，a)：等式(6)。

作为一个基本常识，用于划分非面波区信号和面波区信号的视速度一般小于2000m/s，视速度可以从1000m/s至2000m/s的范围中取值，取值原则以该视速度划分出的面波区能刚好全部囊括所有的面波信号为准。

上述方法可以通过软件、硬件或软件和硬件的结合来实现。例如可以通过计算机仿真技术来实现本发明的方法。所用到的编程语言可以是所属领域常用的计算机语言，例如C/C++、Fortran、Java等。

本发明的实施方式提供的方法在有效压制多分量地震资料中的强面波的同时保护低频弱转换波信号不被压制，与现有技术中的面波压制方法相比，具有以下几个优点：

1、算法基于单道处理设计，适用于任意排序的多分量地震道集的面波压制，适用性强；

2、算法基于单道处理设计，不会带进空间假频，优于传统的f-k法和τ-p面波压制方法；

3、压制面波的同时可以保护有效信号，克服了低信噪比多分量资料面波压制的难题；

4、面波压制效果好，信号保真度高，适用于工业化生产，为多波工业化处理和解释获得高信噪比资料提供了保障。

下面通过举例并结合图2-4来说明使用本发明的实施方式提供的方法在压制面波方面的效果。

图2示出了通过使用本发明的实施方式的方法的单个多分量矢量地震道面波压制效果图。以两分量矢量信号为例，根据本发明的实施方式的方法具体可以描述为：

第一步：输入单个矢量地震道信号，即图2中的图2a和图2e分别对应的原始R分量信号和Z分量信号(可以由同一个数字检波器接收地面震动产生)；

第二步：对每个分量进行时频分解，以小波变换为例，即为图2中的图2b和图2f分别对应的原始R分量和Z分量小波谱；

第三步：求取瞬时椭圆极化率，图2中没有具体对应的示意图。可以表述为：对R、Z小波谱分别抽取每个频率下的小波谱信号，即每个频率下存在一对小波谱信号序列(三分量则存在三个小波谱序列)，采用这一对(或者三个)序列联合计算出属于这一对序列(或者三个序列)的瞬时椭圆极化率(矢量地震道的极化属性)，具体计算公式参见上述步骤。

第四步：根据视速度大小划分非面波区和面波区，分别对应图2中的虚线矩形框和实线矩形框的分界点位置的确定。具体分解点位置的计算公式为：I_p＝abs(Offset)/(V_gRoll·dt)，其中abs(·)表示取绝对值，offset表示该地震道对应的偏移距，V_gRoll表示最大面波速度(可以由操作员人工给定)，dt为采样时间间隔。

第五步：非面波区统计信号均方根能量，面波区能量约束滤波，滤波后各分量的频谱可见图2中的图2c和图2g。

第六步：反变换到时间域得到滤波后的有效信号。如果采用的是小波变换，则反变换指小波反变换，如果采用的是S变换则反变换指S反变换。图2中对应的反变换结果为图2d和图2h。

图3示出了通过使用本发明的实施方式的方法的四川新场地区多分量资料面波压制效果图。使用本发明的方法对四川新场地区的多分量地震资料进行了工业化生产处理，得到了显著的效果。从图3中我们可以看出，本发明的方法确实能最大限度的压制面波信号，提高了多分量地震资料的信噪比。而且，从图3中的图3a和图3b的对比，图3c和图3d的对比可以看出，在压制面波的同时，低频信号被保护得很好，这说明本发明的方法确实能在不伤害有效信号的情况下对面波噪声进行很好的压制。

图4示出了使用本发明的实施方式提供的方法与常用商业软件压制效果对比图，其中从左到右依次表示使用CGG CANCL压制面波的效果，使用冰狐DENOISE压制面波的效果，使用本发明的方法压制面波的效果。资料均采用四川新场地区多分量地震资料。从图4的上三幅图，尤其是椭圆标记区域，我们可以看到，CGG CANCL和本发明的方法对有效信号的保护效果好，冰狐DENOISE保护差；从图4的下三幅图压制掉的噪声，尤其在矩形和椭圆标记区域，可以看出，本发明的方法的面波压制效果实际是CGG CANCL和冰狐DENOISE的一个联合压制效果。

以上实例，均有力地证明了本发明的有益效果：能极大压制低信噪比资料中的强面波，同时保护与其频带重叠的弱转换波信号不被压制，提高多分量资料信噪比。

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

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。

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

Claims (6)

1.一种多分量地震资料面波压制方法，该方法包括：

获取在多个时间采样点采样得到的单个多分量矢量地震道信号S_k(t)，k为正整数；

对所述S_k(t)进行时频双域分解，以得到时频谱信号W_gS_k(t,a)，t表示时间，a表示频率；

对所述时频谱信号W_gS_k(t,a)进行瞬时极化分析，其中该瞬时极化分析包括计算出针对各个时间采样点的瞬时极化椭圆率P(t,a)；

基于视速度大小，将所述W_gS_k(t,a)分为非面波区信号和面波区信号并对所述面波区信号进行面波信号压制处理，以得到压制面波信号后的时频谱信号W_g'S_k(t,a)；

对所述W_g'S_k(t,a)进行反变换并重构以得到压制面波后的矢量地震道信号S'_k(t)；

其中，所述面波信号压制处理包括：

统计出在各个频率下的所述非面波区信号的均方根能量E_nonGR(a)；

针对每个频率下的所述面波区信号，以该频率下的所述均方根能量E_nonGR(a)为能量阈值将所述面波区信号分为强能量信号和弱能量信号后进行分类约束极化滤波。

2.根据权利要求1所述的方法，其中，所述强能量信号采用由等式(1)描述的滤波器进行滤波：

等式(1)，其中F_strong(t,a)表示滤波后的强能量信号，P_strong表示强能量信号所在的极化椭圆率分布范围；

所述弱能量信号采用由等式(2)描述的滤波器进行滤波：

等式(2)，其中F_weak(t,a)表示滤波后的弱能量信号，P_weak表示弱能量信号所在的极化椭圆率分布范围。

3.根据权利要求2所述的方法，其中，所述P_stonrg的取值为P_strong∈(0,0.25]，所述P_waek的取值为P_weak∈[0.45,1)。

4.根据权利要求1所述的方法，其中，在频率a下的所述均方根能量E_nonGR(a)根据等式(3)来计算：

$E_{\mathrm{nonGR}}\left(a\right)=\sqrt{1/M\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{\left|W_g{S}_k(\mathrm{i\&Delta;t},a)\right|}^2}$        等式(3)

其中M为所述非面波区信号的时间采样点个数，Δt为采样时间间隔。

5.根据权利要求1所述的方法，其中，对所述S_k(t)进行时频双域分解包括：

利用复数母小波函数g(t)的连续小波变换根据等式(4)分解所述S_k(t)：

$W_g{S}_k(t,a)=<g_{t,a},S_k>=\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}\frac{1}{a}{g}^{*}\left(\frac{\mathrm{\&tau;}-t}{a}\right)S_k\left(\mathrm{\&tau;}\right)\mathrm{d\&tau;},a\&Element;R,t\&Element;R$    等式(4)

其中，R表示实数域，(·)^*表示复共轭，采用的复数母小波函数其中，ω₀表示复数母小波函数的中心频率。

6.根据权利要求1所述的方法，其中，所述进行瞬时极化分析包括：

根据等式(5)采用希尔伯特变换构造解析时频谱

$W_g{S}_k^C(t,a)=W_g{S}_k(t,a)+\mathrm{iH}\left(W_g{S}_k(t,a)\right)$ 等式(5)，H表示希尔伯特变换；

将所述改写为 $W_g{S}_k^C(t,a)=\stackrel{\&OverBar;}{A}(t,a)\&CenterDot;e^{\mathrm{i\&phi;}(t,a)},$ 表示复数信号的模，φ(t，a)表示t时刻尺度a下复数信号的幅角，并求取使得值最大时的相位函数φ₀(t,a)，n∈N N为正整数，其中

$B(t,a)=\frac{1}{2}\underset{k}{\mathrm{\&Sigma;}}{\left(W_g{S}_k^C(t,a)\right)}^2,C(t,a)=\frac{1}{2}{\left(\underset{k}{\mathrm{\&Sigma;}}{W}_g{S}_k^C(t,a)\right)}^2,\mathrm{\&epsiv;}<<1;$

根据所述相位函数求取极化椭圆的主轴函数R(t,a)和次轴函数r(t,a)，其中， $R(t,a)=\mathrm{Real}[e^{-i{\mathrm{\&phi;}}_0(t,a)}\&CenterDot;W_g{S}_k^C(t,a)],r(t,a)=\mathrm{Real}[e^{-i({\mathrm{\&phi;}}_0(t,a)+\mathrm{\&pi;}/2)}\&CenterDot;W_g{S}_k^C(t,a)],$ Real表示取复数信号的实部；以及

根据等式(6)计算所述瞬时极化椭圆率P(t,a)：等式(6)。