CN103376462A - 一种自动检测强能量噪声的方法 - Google Patents

Abstract

本发明是石油勘探地震数据处理中自动识别强能量噪声的方法，对目标道集的每道数据沿时间方向做傅氏变换，计算道集数据振幅谱，抽取任意频点序列，计算点值序列最大值、中值和概率密度等参数，完成所有频率样点值序列的强能量噪声检测，得到强能量噪声的频率域道集数据，沿时间方向做反傅氏变换，得到返回时间域的强能量噪声道集数据；本发明使用最大似然准则在变换域中通过概率统计自动搜索检测强能量噪声，对于在不同的道集、不同的时窗中强能量噪声振幅变化的情况，能够准确检测出强能量噪声。

Description

一种自动检测强能量噪声的方法

技术领域

本发明涉及石油勘探技术，是地震数据处理中一种自动检测强能量噪声的方法。

背景技术

在地震数据的采集过程中，由于各种因素的影响，使得原始记录中出现了不同种类的干扰波，如声波干扰、50Hz交流电干扰、碎发脉冲噪声、在记录深部出现的高频干扰及涌浪噪声等。这些干扰波具有相对较强的能量，占据一定的频带范围。常用的处理方法是滤波或道剔除，但无论是交互剔除废道还是人工编辑，都要花费大量的时间，而且淹没在能量干扰中的有效信号也损失掉了。这种处理方法已不能满足高分辨率勘探的精度要求。高分辨率勘探对资料处理的精度提出了更高的要求，特别是目前广泛使用的叠前多道处理技术(如地表一致性振幅补偿、统计子波反褶积等)要求输人数据有较高的质量，强能量噪声处理不当直接影响着这些方法的效果，限制了某些技术的应用。

目前工业上常规的强能量噪声衰减方法，是在一定的变换域(如时频域)中，按一定的方向在一定的时窗内得到样点振幅的中值、平均值或均方根值，通过试验得到一个门槛因子，门槛因子乘中值、平均值或均方根值得到门槛值，样点振幅值大于门槛值，则该样点值被视为强能量噪声，对其进行衰减处理；但是强能量噪声在不同的道集，在不同的时窗可能不同，此时应用小的门槛因子会伤害有效信号，应用大的门槛因子强能量噪声衰减不完全。

发明内容

本发明目的是提供一种能够很好将记录中的强能量噪声自动检测出来，将检测出的强能量噪声振幅值置零或旁道替换则达到衰减强能量噪声的自动检测强能量噪声的方法。

本发明通过以下技术方案实现，具体实施步骤是：

1)激发并采集地震数据；

2、根据权利要求1的方法，特征是步骤1)所述的地震数据是从采集的地震数据中任意取出一个地震道数据作为目标道集。

2)对目标道集的每道数据沿时间方向做傅氏变换，得到频率域道集数据；

3)计算频率域道集数据的振幅值，得到道集数据振幅谱；

4)抽取道集数据振幅谱的任意一个频率样点值序列；

5)计算频率样点值序列的最大值和中值；

6)按照以下公式计算频率样点值序列中值对应先验概率密度p₀(i)：

$p_0\left(i\right)=\frac{1}{{\mathrm{\λ}}_0}\exp (-\frac{a\left(i\right)}{{\mathrm{\λ}}_0}),$ 式中：

i是频率样点值序列的序号；

λ₀是步骤5)中得到的频率样点值序列的中值；

a(i)是步骤5)中得到的频率样点值序列。

7)按照以下公式计算频率样点值序列最大值对应先验概率密度p₁(i)：

$p_1\left(i\right)=\frac{1}{{\mathrm{\λ}}_1}\exp (-\frac{a\left(i\right)}{{\mathrm{\λ}}_1}),$ 式中：

i是频率样点值序列的序号；

λ₁是步骤5)中得到的频率样点值序列的最大值；

a(i)是步骤5)中得到的频率样点值序列；

8)按照以下公式计算后验概率U(i)：

$U\left(i\right)=\frac{{\mathrm{\ωp}}_1\left(i\right)}{{\mathrm{\ωp}}_1\left(i\right)+(1-\mathrm{\ω})p_0\left(i\right)},$ 式中：

i是频率样点值序列的序号；

p₀(i)步骤6)是中得到的频率样点值序列中值对应先验概率密度；

p₁(i)步骤7)是中得到的频率样点值序列最大值对应先验概率密度；

ω是调试参数，表示强能量噪声样点数在总体样点数中的百分比；

步骤8)所述调试参数初始值等于0.1；

9)按照以下公式计算调试参数更新的强能量样点数在总体样点数中的百分比ω′、有效信号样点和总体样点振幅值方差λ′₀和强能量噪声样点和总体样点振幅值方差λ′₁：

${\mathrm{\ω}}^{\′}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}U\left(i\right);$

${\mathrm{\λ}}_1^{\′}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(U\left(i\right)a\left(i\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}U\left(i\right));$

${\mathrm{\λ}}_1^{\′}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(1-U\left(i\right))a\left(i\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(1-U\left(i\right));$

式中：

i是频率样点值序列的序号；

N是频率样点值序列的样点总数；

U(i)是步骤8)中得到的后验概率；

10)如果ω′和ω之差小于0.0001，λ′₁和λ₁之差小于0.001，同时λ′₀和λ₀之差小于0.001，则进行下一步，

否则将：

强能量样点数在总体样点数中的百分比ω′的值赋给步骤8)表示强能量噪声样点数在总体样点数中的百分比ω；

有效信号样点和总体样点振幅值方差λ′₀的值赋给步骤5)中得到的频率样点值序列的中值λ₀；

强能量噪声样点和总体样点振幅值方差λ′₁的值赋给步骤5)中得到的频率样点值序列的最大值λ₁；

11)按照以下公式计算相关系数R(i)：

$R\left(i\right)=|\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}{\left(\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}A(k,l)\right)}^2-\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}{A}^2(k,l)|/2*\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}{A}^2(k,l),$ 式中：

i是频率样点值序列的序号，j是空间样点值序列的序号；

A是步骤3)中得到的道集数据振幅谱；

k是频率方向变量，1是空间方向变量；

12)计算相关加权后验概率

$\stackrel{\‾}{U}\left(i\right)=\frac{1}{R\left(i\right)+1}U\left(i\right),$ 式中：

i是频率样点值序列的序号；

R(i)是步骤11)中得到的相关系数；

U(i)是步骤8)中得到的后验概率；

13)相关加权后验概率大于0.5的频率样点为强能量噪声样点，对步骤2)中相关加权后验概率小于或等于0.5的频率样点的频率值置零，得到强能量噪声的频率域道集数据；

14)按上述4)至13)步完成所有频率样点值序列的强能量噪声检测，得到强能量噪声的频率域道集数据；

15)对步骤14)强能量噪声的频率域道集数据每道数据沿时间方向做反傅氏变换，得到返回时间域的强能量噪声道集数据。

本发明不需要试验门槛因子，在时频域中通过概率统计自动检测强能量的门槛值，使用最大似然准则自动统计搜索强能量门槛值，能够适应强能量门槛值变化的情况，能够很好检测强能量噪声，在保护有效信号的同时衰减强能量噪声。

附图说明

图1为一个炮集海洋地震数据，由408道组成，横坐标为道数(道)，纵坐标为时间(秒)；

图2为一个时窗数据及其振幅谱，a)是图1中炮集数据抽取的一个时窗数据，b)是a)对应数据多道振幅谱，c)是a)对应数据平均振幅谱；其中a)的横坐标为道数(道)，纵坐标为时间(秒)；b)的横坐标为道数(道)，纵坐标为频率(赫兹)；c)的横坐标为频率(赫兹)，纵坐标为归一化振幅值；

图3为6Hz频率样点序列振幅谱及其相关加权概率，a)是6Hz频率样点序列振幅谱，b)是6Hz频率样点序列相关加权概率；a)的横坐标为道数(道)，纵坐标为振幅值；b)的横坐标为道数(道)，纵坐标为相关加权概率；

图4为20Hz频率样点序列振幅谱及其相关加权概率，a)是20Hz频率样点序列振幅谱，b)是20Hz频率样点序列相关加权概率；a)的横坐标为道数(道)，纵坐标为振幅值；b)的横坐标为道数(道)，纵坐标为相关加权概率；；

图5a)是图2a)时窗数据衰减涌浪噪声后数据，横坐标为道数(道)，纵坐标为时间(秒)；b)是衰减掉的涌浪噪声，横坐标为道数(道)，纵坐标为时间(秒)；c)是a)数据的多道振幅谱横坐标为道数(道)，纵坐标为频率(赫兹)；d)是a)数据的平均振幅谱，横坐标为频率(赫兹)，纵坐标为归一化振幅值；

图6a)是图1炮集数据衰减涌浪噪声后数据，横坐标为道数(道)，纵坐标为时间(秒)；b)是图1中检测并衰减掉的涌浪噪声，横坐标为道数(道)，纵坐标为时间(秒)；

图7a)是陆地单炮数据，横坐标为道数(道)，纵坐标为时间(秒)；b)b)是从a)数据中检测到的强能量噪声数据，横坐标为道数(道)，纵坐标为时间(秒)；c)是将检测到的强能量噪声数据从a)数据减去后的数据，横坐标为道数(道)，纵坐标为时间(秒)。

具体实施方式

以下结合附图和通过实例详细说明本发明。

本发明采用数据是一个炮集的海洋地震数据，其中有较强的涌浪噪声(涌浪噪声是一类低频强能量噪声)，408道，每道4000个采样点，采样间隔为2ms，如图1所示。通过以下步骤完成该炮集数据强能量噪声自动检测：

1)按时间方向和空间方向对炮集数据进行分块；

2)抽取任意一个分块数据(如图2-a所示)，对分块数据每道数据沿时间方向做傅氏变换，得到频率域分块数据；

3)计算频率域分块数据的振幅值，得到分块数据振幅谱；

4)抽取分块数据振幅谱的任意一个频率样点值序列a(i)，i是频率样点值序列的序号，i＝1，2，…，N(如图3-a和图4-a所示)；

5)计算频率样点值序列的最大值λ₁和中值λ₀；

6)计算频率样点值序列中值对应先验概率密度

7)计算频率样点值序列最大值对应先验概率密度

8)取ω等于0.1，计算后验概率

9)计算 ${\mathrm{\ω}}^{\′}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}U\left(i\right),$ ${\mathrm{\λ}}_1^{\′}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(U\left(i\right)a\left(i\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}U\left(i\right))$

和 ${\mathrm{\λ}}_1^{\′}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(1-U\left(i\right))a\left(i\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(1-U\left(i\right));$

10)如果ω′和ω之差小于0.0001，λ′₁和λ₁之差小于0.001，同时λ′₀和λ₀之差小于0.001，则进入到第11)步，否则将ω′赋给ω，λ′₁武给λ₁，λ′₀赋给λ₀，返回第6)步；

11)计算相关系数 $R\left(i\right)=|\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}{\left(\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}A(k,l)\right)}^2-\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}{A}^2(k,l)|/2*\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}{A}^2(k,l);$

12)计算相关加权后验概率

(如图3-b和图4-b所示)；

13)相关加权后验概率大于0.5的频率样点为强能量噪声样点，对步骤2)中相关加权后验概率小于或等于0.5的频率样点的频率值置零，得到强能量噪声的频率域道集数据；

14)按上述4)至13)步完成所有频率样点值序列的强能量噪声检测，得到强能量噪声的频率域分块数据；

15)对步骤14)强能量噪声的频率域道集数据每道数据沿时间方向做反傅氏变换，得到返回时间强能量噪声的分块数据(如图5-b所示)；

按上述2)至15)步骤完成所有分块数据的强能量噪声检测，得到整个炮集数据检测到的强能量噪声数据(如图6-b所示)，将检测到的强能量噪声从初始炮集(如图1所示)减去后的数据(如图6-a所示)可以看到涌浪噪声得到了很好衰减，说明本发明方法对于检测涌浪噪声的有效性。

为了检验自动检测涌浪噪声的效果，对图2-a中抽取时窗数据进行多道振幅谱(如图2-b所示)和平均振幅谱(如图2-c所示)计算，从多道振幅谱或平均振幅谱都可以看到抽取时窗数据低频段存在较强涌浪噪声；对衰减涌浪噪声后的时窗数据进行多道振幅谱(如图5-c所示)和平均振幅谱(如图5-d所示)计算，从多道振幅谱或平均振幅谱都可以看到抽取时窗数据低频段的较强涌浪噪声得到有效衰减，从衰减掉涌浪噪声后的数据(如图5-a所示)可以看到有效信号得到保护。

本发明采用另一个数据是一个陆上炮集地震数据(如图7-a所示)，其中有较强50Hz交流电干扰、高频噪声等强能量噪声；数据400道，采样点数为2500，采样间隔为2ms。通过以下步骤完成该炮集数据强能量噪声自动检测：

1)按空间方向对炮集数据进行分块；

2)抽取任意一个分块数据，对分块数据每道数据沿时间方向做傅氏变换，得到频率域分块数据；

3)计算频率域分块数据的振幅值，得到分块数据振幅谱；

4)抽取分块数据振幅谱的任意一个频率样点值序列a(i)，i是频率样点值序列的序号，i＝1，2，…，N；

5)计算频率样点值序列的最大值λ₁和中值λ₀；

6)计算频率样点值序列中值对应先验概率密度

7)计算频率样点值序列最大值对应先验概率密度

8)取ω等于0.1，计算后验概率

9)计算 ${\mathrm{\ω}}^{\′}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}U\left(i\right),$ ${\mathrm{\λ}}_1^{\′}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(U\left(i\right)a\left(i\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}U\left(i\right))$

和 ${\mathrm{\λ}}_1^{\′}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(1-U\left(i\right))a\left(i\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(1-U\left(i\right));$

10)如果ω′和ω之差小于0.0001，λ′₁和λ₁之差小于0.001，同时λ′₀和λ₀之差小于0.001，则进入到第11)步，否则将ω′赋给ω，λ′₁武给λ₁，λ′₀赋给λ₀，返回第6)步；

11)计算相关系数 $R\left(i\right)=|\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}{\left(\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}A(k,l)\right)}^2-\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}{A}^2(k,l)|/2*\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}{A}^2(k,l);$

12)计算相关加权后验概率

13)相关加权后验概率

大于0.5的频率样点为强能量噪声样点，对步骤2)中相关加权后验概率小于或等于0.5的频率样点的频率值置零，得到强能量噪声的频率域道集数据；

14)按上述4)至13)步完成所有频率样点值序列的强能量噪声检测，得到强能量噪声的频率域分块数据；

15)对步骤14)检测到的强能量噪声的频率域道集数据每道数据沿时间方向做反傅氏变换，得到返回时间强能量噪声的分块数据；

按上述2)至15)步骤完成所有分块数据的强能量噪声检测，得到整个炮集数据检测到的强能量噪声数据(如图7-b所示)，将检测到的强能量噪声从初始炮集(如图7-a所示)减去后的数据(如图7-c所示)可以看到强能量噪声得到了很好衰减，说明本发明方法对于检测陆地数据中强能量噪声的有效性。

本发明不需要试验门槛因子，在时频域中通过概率统计自动检测强能量的门槛值，使用最大似然准则自动统计搜索强能量门槛值，对于在不同的道集、不同的时窗中强能量噪声振幅变化的情况，也能够很好检测出强能量噪声，将检测出的强能量噪声从记录作衰减处理，达到压制强能量噪声的目的。

Claims (6)

1.一种自动检测强能量噪声的方法，特征是通过以下具体步骤实现：

1)激发并采集地震数据；

2)对目标道集的每道数据沿时间方向做傅氏变换，得到频率域道集数据；

3)计算频率域道集数据的振幅值，得到道集数据振幅谱；

4)抽取道集数据振幅谱的任意一个频率样点值序列；

5)计算频率样点值序列的最大值和中值；

6)按照以下公式计算频率样点值序列中值对应先验概率密度p₀(i)：

$p_0\left(i\right)=\frac{1}{{\mathrm{\λ}}_0}\exp (-\frac{a\left(i\right)}{{\mathrm{\λ}}_0}),$ 式中：

i是频率样点值序列的序号；

λ₀是步骤5)中得到的频率样点值序列的中值；

a(i)是步骤5)中得到的频率样点值序列；

7)按照以下公式计算频率样点值序列最大值对应先验概率密度p₁(i)：

$p_1\left(i\right)=\frac{1}{{\mathrm{\λ}}_1}\exp (-\frac{a\left(i\right)}{{\mathrm{\λ}}_1}),$ 式中：

i是频率样点值序列的序号；

λ₁是步骤5)中得到的频率样点值序列的最大值；

a(i)是步骤5)中得到的频率样点值序列；

8)按照以下公式计算后验概率U(i)：

9)按照以下公式计算调试参数更新的强能量样点数在总体样点数中的百分比ω′、有效信号样点和总体样点振幅值方差λ′₀和强能量噪声样点和总

体样点振幅值方差λ′₁：

${\mathrm{\ω}}^{\′}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}U\left(i\right);$

${\mathrm{\λ}}_1^{\′}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(U\left(i\right)a\left(i\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}U\left(i\right));$

${\mathrm{\λ}}_1^{\′}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(1-U\left(i\right))a\left(i\right)/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(1-U\left(i\right));$

式中：

i是频率样点值序列的序号；

N是频率样点值序列的样点总数；

U(i)是步骤8)中得到的后验概率；

10)如果ω′和ω之差小于0.0001，λ′₁和λ₁之差小于0.001，同时λ′₀和λ₀之差小于0.001，则进行下一步，

ω是调试参数，表示强能量噪声样点数在总体样点数中的百分比；

否则将：

强能量样点数在总体样点数中的百分比ω′的值赋给步骤8)表示强能量噪声样点数在总体样点数中的百分比ω；

有效信号样点和总体样点振幅值方差λ′₀的值赋给步骤5)中得到的频率样点值序列的中值λ₀；

强能量噪声样点和总体样点振幅值方差λ′₁的值赋给步骤5)中得到的频率样点值序列的最大值λ₁；

11)计算相关系数R(i)：

12)计算相关加权后验概率

13)相关加权后验概率大于0.5的频率样点为强能量噪声样点，对步骤2)中相关加权后验概率小于或等于0.5的频率样点的频率值置零，得到强能量噪声的频率域道集数据；

14)按上述4)至13)步完成所有频率样点值序列的强能量噪声检测，得到强能量噪声的频率域道集数据；

15)对步骤14)强能量噪声的频率域道集数据每道数据沿时间方向做反傅氏变换，得到返回时间域的强能量噪声道集数据。

2.根据权利要求1的方法，特征是步骤1)所述的地震数据是从采集的地震数据中任意取出一个地震道数据作为目标道集。

3.根据权利要求1的方法，特征是步骤8)按照以下公式计算后验概率U(i)：

$U\left(i\right)=\frac{{\mathrm{\ωp}}_1\left(i\right)}{{\mathrm{\ωp}}_1\left(i\right)+(1-\mathrm{\ω})p_0\left(i\right)},$ 式中：

i是频率样点值序列的序号；

p₀(i)步骤6)是中得到的频率样点值序列中值对应先验概率密度；

p₁(i)步骤7)是中得到的频率样点值序列最大值对应先验概率密度；

ω是调试参数，表示强能量噪声样点数在总体样点数中的百分比。

4.根据权利要求3的方法，特征是所述调试参数ω初始值等于0.1。

5.根据权利要求1的方法，特征是步骤11)按照以下公式计算相关系数R(i)：

$R\left(i\right)=|\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}{\left(\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}A(k,l)\right)}^2-\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}{A}^2(k,l)|/2*\underset{k=i-1}{\overset{i+1}{\mathrm{\Σ}}}\underset{l=j-1}{\overset{j+1}{\mathrm{\Σ}}}{A}^2(k,l),$ 式中：

i是频率样点值序列的序号，j是空间样点值序列的序号；

A是步骤3)中得到的道集数据振幅谱；

k是频率方向变量，1是空间方向变量。

6.根据权利要求1的方法，特征是步骤12)计算相关加权后验概率

采用以下公式：

$\stackrel{\‾}{U}\left(i\right)=\frac{1}{R\left(i\right)+1}U\left(i\right),$ 式中：

i是频率样点值序列的序号；

R(i)是步骤11)中得到的相关系数；

U(i)是步骤8)中得到的后验概率。

Images