CN109901224A - 一种地震资料低频信号保护压制噪声方法 - Google Patents

Abstract

本发明公开了一种地震资料低频信号保护压制噪声方法，包括如下步骤：1)数据重排：根据信号的特点，重新选排道集；2)分频带：根据噪声所在的频带范围，通过MP谱分解成高、中、低频频带；3)在低频带内进行信噪分离；4)匹配相减，数据合成，输出去噪后的地震记录。本发明采用匹配追踪(MP，matching pursuit)算法，将地震记录进行分频带，对低频带的噪声能量进行预测压制，从而保护低频带内的有效信号，而不是采用上述带通滤波技术，在去除噪声的同时，有效的弱信号也消除掉了，采用本发明可有效消除低频端的噪声，而有效信号却不损失。该技术算法保真度高，参数易于选择。

Description

一种地震资料低频信号保护压制噪声方法

技术领域

本发明涉及保护低频信号技术领域，尤其是涉及一种地震资料低频信号保护压制噪声方法。

背景技术

低频信号(0-5Hz)对提高地震资料分辨率和反演稳定性具有非常重要的作用。但低频端的信息信噪比低，在地震资料处理过程中，往往采用带通滤波处理，在去掉低频噪声的过程时，微弱的信号也去掉了。保护低频信号处理是现代地震资料处理的重要研究内容之一，是当前研究上午热点和难点。

带通滤波技术是一种可以压制低频端噪声的常用技术。带通滤波能通过某一频率范围内的频率分量、将其他范围的频率分量衰减到极低水平的滤波器，因此，如设计0-5～130-150Hz的带通滤波器，可以压制0-5Hz的以下的频率分量，但这样，压制0-5Hz噪声的同时，0-5Hz的有用信号也消除了。

本申请是通过高精度谱分解技术进行分频带，在噪声所在的低频带区间内统计能量，利用均值加权思想来预测噪声，确定噪声所在频率区间和在相邻道间的能量和时间差异，并予以消除。

公开于该背景技术部分的信息仅仅旨在加深对本申请的总体背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。

发明内容

本发明的目的在于提供一种地震资料低频信号保护压制噪声方法，以解决现有技术中存在的技术问题。

为了实现上述目的，本发明采用以下技术方案：

本发明提供一种地震资料低频信号保护压制噪声方法，包括如下步骤：

1)数据重排：根据信号的特点，重新选排道集；

2)分频带：根据噪声所在的频带范围，通过MP谱分解成高、中、低频频带；

3)在低频带内进行信噪分离，具体如下：

①对于给定的地震道，给定其空间窗内的道数，估计该地震道在该空间窗的能量，估计方式可以采用中值、均值或者中值和均值的算法组合；

②计算原始地震道与对应时间窗的能量比例系数；

③对高出这个比例系数的值提取，作为噪声模型道；

4)匹配相减，数据合成，输出去噪后的地震记录。

作为一种进一步的技术方案，步骤2)具体为：

Morlet小波定义为：

式中ω_m为平均角频率，σ为尺度，φ为相移，u为时间延迟，对于Morlet小波，利用四个参数进行刻画，γ＝{u,σ,ω,φ}；

匹配追踪算法迭代地完成，每步迭代相应地提取最优子波在进行N步迭代后，地震道f(t)展开为如下形式：

其中，a_n是第n个子波的振幅，R^(N)f是残差，R⁽⁰⁾f＝f，在第n步分解迭代中，采用三个步骤实现：(1)有效地近似估计γ_n＝{u_n,σ_n,ω_n,φ_n}四个参数；(2)更新四个参数得到最优子波(3)估计振幅a_n，对于所有迭代三个步骤重复进行，具体说来:

第一步：利用复地震道属性估计{u_n,ω_n,φ_n}三个参数，设置复数道最大包络的时间为时间延迟u_n，瞬时频率为中心频率ω_n，瞬时相位为φ_n，再搜索第四个参数σ_n，利用下面方程：

其中，D＝{g_γ(t)}_γ∈Γ是构成子波的冗余字典，<f,h>为内积，为归一化子波利用方程(其中，D是地震记录，W为地震子波，R为地层反射系数)对固定u_n,ω_n,φ_n的均匀分布的σ值来搜寻最优参数σ_n；

第二步，利用方程在子字典中搜索这四个参数，搜索范围ξ为[ξ-△ξ,ξ+△ξ]；

第三步，估计最优子波的振幅，利用如下方程：

对于谱分解的当前迭代，其中正交于R⁽ⁿ⁺¹⁾f，有

对于方程的正交投影，我们选择使得最大，因此残差||R⁽ⁿ⁺¹⁾||f²最小，因此方程(4)满足最小残差条件：

在N步迭代后，最终的最小化残差R^(N)f^(t)认为是数据噪声；

将信号f(t)分解为一系列的子波后，在时频面信号f(t)的时频空间的振幅谱：

通过对信号进行匹配追踪谱分解，可以得到信号不同频率的频率切片。

作为一种进一步的技术方案，步骤3)中在低频带内进行信噪分离的具体计算方法如下：

记地震信号为xr(t)，利用希尔伯特变换，地震信号可以由瞬时相位和瞬时包络重构：

x_r(t)＝cosθ(t)·a(t) (8)

其中x_h(t)为信号的希尔伯特变换；因此，利用瞬时包络来统计全局平均能量Eglobal和给定时窗内的局部平均能量Elocal：

根据局部平均和全局平均能量，可以计算噪声检测系数：计算时窗沿着时间方向和空间方向滑动，既可计算出噪声检测函数k(i,j)，对函数k(i,j)进行空间平滑滤波，滤波器采用多点中值或平均函数f_n，得到检测标准函数g(k)：

g(k)＝f_n(k) (11)

对于信号，某一道的某一时刻，噪声检测函数k与检测标准函数g的比值在1附近变化，而记录中某些道上存在较强能量的噪声时，这个比值会显著变大，这样通过定义不同的门槛值C，可以设计不同能量的噪声压制函数：

α＝k(i,j)/g(i,j)

作为一种进一步的技术方案，步骤4)具体为：

将噪声压制函数作用于二维信号x(i,j),即得到去噪后记录：

x'(i,j)＝x(i,j)/s(i,j) (13)

对低频带内的记录进行噪声检测及衰减，然后再重构信号；对于分频处理，常用多个频带同时处理，然后再统一重构的方式：

这种方式会造成每个频带间因为截断而产生吉布斯效应，即使用窗函数也会造成重构信号的失真；为了避免分频重构时的滤波效应，采取多次迭代的方式，首先通过带通滤波得到要处理的某一频带的记录xf(i,j),相对于这一频带，可以得到不处理部分的记录xsave(i,j)：

x_save(i,j)＝x(i,j)-x_f(i,j) (15)

对低频带的记录进行噪声压制，然后和不需要处理的部分重构，得到新的记录：

通过上式，便得到了噪声部分，然后和原始数据相减，便得到低频带内最终去噪后的结果。

采用上述技术方案，本发明具有如下有益效果：

本发明采用匹配追踪(MP，matching pursuit)算法，将地震记录进行分频带，对低频带的噪声能量进行预测压制，从而保护低频带内的有效信号，而不是采用上述带通滤波技术，在去除噪声的同时，有效的弱信号也消除掉了，采用本发明可有效消除低频端的噪声，而有效信号却不损失。该技术算法保真度高，参数易于选择。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单的介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为某地区的实际地震原始记录波形图；

图2为通过常规带通滤波技术对图1中的波形图进行处理后的波形图；

图3为通过本发明的方法对图1中的波形图进行处理后的波形图。

具体实施方式

下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。

在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

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

本实施例提供一种地震资料低频信号保护压制噪声方法，包括如下步骤：

1)数据重排：根据信号的特点，重新选排道集；

2)分频带：根据噪声所在的频带范围，通过MP谱分解成高、中、低频频带；

3)在低频带内进行信噪分离，具体如下：

①对于给定的地震道，给定其空间窗内的道数，估计该地震道在该空间窗的能量，估计方式可以采用中值、均值或者中值和均值的算法组合；

②计算原始地震道与对应时间窗的能量比例系数；

③对高出这个比例系数的值提取，作为噪声模型道；

4)匹配相减，数据合成，输出去噪后的地震记录。

作为一种进一步的技术方案，步骤2)具体为：

Morlet小波定义为：

式中ω_m为平均角频率，σ为尺度，φ为相移，u为时间延迟，对于Morlet小波，利用四个参数进行刻画，γ＝{u,σ,ω,φ}；

匹配追踪算法迭代地完成，每步迭代相应地提取最优子波在进行N步迭代后，地震道f(t)展开为如下形式：

其中，a_n是第n个子波的振幅，R^(N)f是残差，R⁽⁰⁾f＝f，在第n步分解迭代中，采用三个步骤实现：(1)有效地近似估计γ_n＝{u_n,σ_n,ω_n,φ_n}四个参数；(2)更新四个参数得到最优子波(3)估计振幅a_n，对于所有迭代三个步骤重复进行，具体说来:

第一步：利用复地震道属性估计{un,ωn,φ_n}三个参数，设置复数道最大包络的时间为时间延迟u_n，瞬时频率为中心频率ω_n，瞬时相位为φ_n，再搜索第四个参数σ_n，利用下面方程：

其中，D＝{g_γ(t)}_γ∈Γ是构成子波的冗余字典，<f,h>为内积，为归一化子波利用方程(其中，D是地震记录，W为地震子波，R为地层反射系数)对固定u_n,ω_n,φ_n的均匀分布的σ值来搜寻最优参数σ_n；

第二步，利用方程在子字典中搜索这四个参数，搜索范围ξ为[ξ-△ξ,ξ+△ξ]；

第三步，估计最优子波的振幅，利用如下方程：

对于谱分解的当前迭代，其中正交于R⁽ⁿ⁺¹⁾f，有

对于方程的正交投影，我们选择使得最大，因此残差||R⁽ⁿ⁺¹⁾f||²最小，因此方程(4)满足最小残差条件：

在N步迭代后，最终的最小化残差R^(N)f(t)认为是数据噪声；

将信号f(t)分解为一系列的子波后，在时频面信号f(t)的时频空间的振幅谱：

通过对信号进行匹配追踪谱分解，可以得到信号不同频率的频率切片。

作为一种进一步的技术方案，步骤3)中在低频带内进行信噪分离的具体计算方法如下：

记地震信号为xr(t)，利用希尔伯特变换，地震信号可以由瞬时相位和瞬时包络重构：

x_r(t)＝cosθ(t)·a(t) (8)

其中x_h(t)为信号的希尔伯特变换；因此，利用瞬时包络来统计全局平均能量Eglobal和给定时窗内的局部平均能量Elocal：

根据局部平均和全局平均能量，可以计算噪声检测系数：计算时窗沿着时间方向和空间方向滑动，既可计算出噪声检测函数k(i,j)，对函数k(i,j)进行空间平滑滤波，滤波器采用多点中值或平均函数f_n，得到检测标准函数g(k)：

g(k)＝f_n(k) (11)

对于信号，某一道的某一时刻，噪声检测函数k与检测标准函数g的比值在1附近变化，而记录中某些道上存在较强能量的噪声时，这个比值会显著变大，这样通过定义不同的门槛值C，可以设计不同能量的噪声压制函数：

α＝k(i,j)/g(i,j)

作为一种进一步的技术方案，步骤4)具体为：

将噪声压制函数作用于二维信号x(i,j),即得到去噪后记录：

x'(i,j)＝x(i,j)/s(i,j) (13)

对低频带内的记录进行噪声检测及衰减，然后再重构信号；对于分频处理，常用多个频带同时处理，然后再统一重构的方式：

这种方式会造成每个频带间因为截断而产生吉布斯效应，即使用窗函数也会造成重构信号的失真；为了避免分频重构时的滤波效应，采取多次迭代的方式，首先通过带通滤波得到要处理的某一频带的记录xf(i,j),相对于这一频带，可以得到不处理部分的记录xsave(i,j)：

x_save(i,j)＝x(i,j)-x_f(i,j) (15)

对低频带的记录进行噪声压制，然后和不需要处理的部分重构，得到新的记录：

通过上式，便得到了噪声部分，然后和原始数据相减，便得到低频带内最终去噪后的结果。

本发明是基于地震信号在具有多次覆盖数据集中是相干的，振幅变化是光滑的，地震信号包络展现信号的振幅变化，信号包络是光滑的；且噪声具有随机性，不连续，因此采用统计均值能量做为标准，来预测噪声，并设计系数衰减噪声。

结合图1-3所示，图1是某地区的实际地震原始记录，有很强的低频噪音干扰；应用本发明提供的处理技术后，有效压制了低频噪音，地质信息更加丰富(明显优于图2中的处理结果)，说明该技术具有较好的振幅保真性，明显地增强了地震资料的信噪比和保真度。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (4)

1.一种地震资料低频信号保护压制噪声方法，其特征在于，包括如下步骤：

1)数据重排：根据信号的特点，重新选排道集；

2)分频带：根据噪声所在的频带范围，通过MP谱分解成高、中、低频频带；

3)在低频带内进行信噪分离，具体如下：

①对于给定的地震道，给定其空间窗内的道数，估计该地震道在该空间窗的能量，估计方式可以采用中值、均值或者中值和均值的算法组合；

②计算原始地震道与对应时间窗的能量比例系数；

③对高出这个比例系数的值提取，作为噪声模型道；

4)匹配相减，数据合成，输出去噪后的地震记录。

2.根据权利要求1所述的地震资料低频信号保护压制噪声方法，其特征在于，步骤2)具体为：

Morlet小波定义为：

式中ω_m为平均角频率，σ为尺度，φ为相移，u为时间延迟，对于Morlet小波，利用四个参数进行刻画，γ＝{u,σ,ω,φ}；

匹配追踪算法迭代地完成，每步迭代相应地提取最优子波在进行N步迭代后，地震道f(t)展开为如下形式：

其中，a_n是第n个子波的振幅，R^(N)f是残差，R⁽⁰⁾f＝f，在第n步分解迭代中，采用三个步骤实现：(1)有效地近似估计γ_n＝{u_n,σ_n,ω_n,φ_n}四个参数；(2)更新四个参数得到最优子波(3)估计振幅a_n，对于所有迭代三个步骤重复进行，具体说来:

第一步：利用复地震道属性估计{u_n,ω_n,φ_n}三个参数，设置复数道最大包络的时间为时间延迟u_n，瞬时频率为中心频率ω_n，瞬时相位为φ_n，再搜索第四个参数σ_n，利用下面方程：

其中，D＝{g_γ(t)}_γ∈Γ是构成子波的冗余字典，<f,h>为内积，为归一化子波利用方程(其中，D是地震记录，W为地震子波，R为地层反射系数)对固定u_n,ω_n,φ_n的均匀分布的σ值来搜寻最优参数σ_n；

第二步，利用方程在子字典中搜索这四个参数，搜索范围ξ为[ξ-△ξ,ξ+△ξ]；

第三步，估计最优子波的振幅，利用如下方程：

对于谱分解的当前迭代，其中正交于R⁽ⁿ⁺¹⁾f，有

对于方程的正交投影，我们选择使得最大，因此残差||R⁽ⁿ⁺¹⁾f||²最小，因此方程(4)满足最小残差条件：

在N步迭代后，最终的最小化残差R^(N)f(t)认为是数据噪声；

将信号f(t)分解为一系列的子波后，在时频面信号f(t)的时频空间的振幅谱：

通过对信号进行匹配追踪谱分解，可以得到信号不同频率的频率切片。

3.根据权利要求1所述的地震资料低频信号保护压制噪声方法，其特征在于，步骤3)中在低频带内进行信噪分离的具体计算方法如下：

记地震信号为xr(t)，利用希尔伯特变换，地震信号可以由瞬时相位和瞬时包络重构：

x_r(t)＝cosθ(t)·a(t) (8)

其中x_h(t)为信号的希尔伯特变换；因此，利用瞬时包络来统计全局平均能量Eglobal和给定时窗内的局部平均能量Elocal：

根据局部平均和全局平均能量，可以计算噪声检测系数：计算时窗沿着时间方向和空间方向滑动，既可计算出噪声检测函数k(i,j)，对函数k(i,j)进行空间平滑滤波，滤波器采用多点中值或平均函数f_n，得到检测标准函数g(k)：

g(k)＝f_n(k) (11)

对于信号，某一道的某一时刻，噪声检测函数k与检测标准函数g的比值在1附近变化，而记录中某些道上存在较强能量的噪声时，这个比值会显著变大，这样通过定义不同的门槛值C，可以设计不同能量的噪声压制函数：

α＝k(i,j)/g(i,j)

4.根据权利要求1所述的地震资料低频信号保护压制噪声方法，其特征在于，步骤4)具体为：

将噪声压制函数作用于二维信号x(i,j),即得到去噪后记录：

x'(i,j)＝x(i,j)/s(i,j) (13)

对低频带内的记录进行噪声检测及衰减，然后再重构信号；对于分频处理，常用多个频带同时处理，然后再统一重构的方式：

这种方式会造成每个频带间因为截断而产生吉布斯效应，即使用窗函数也会造成重构信号的失真；为了避免分频重构时的滤波效应，采取多次迭代的方式，首先通过带通滤波得到要处理的某一频带的记录xf(i,j),相对于这一频带，可以得到不处理部分的记录xsave(i,j)：

x_save(i,j)＝x(i,j)-x_f(i,j) (15)

对低频带的记录进行噪声压制，然后和不需要处理的部分重构，得到新的记录：

通过上式，便得到了噪声部分，然后和原始数据相减，便得到低频带内最终去噪后的结果。