CN113376698B - 一种低漏检的地震数据尖峰噪声检测及压制方法 - Google Patents

Abstract

本发明涉及一种低漏检的地震数据尖峰噪声检测及压制方法，首先应用高阈值对尖峰噪声进行检测，以保证有用信号的低误检；然后以已检测到的尖峰噪声为中心，取一定长度的窗口，应用低阈值对窗口内的数据进行二次检测，从而保证尖峰噪声的低漏检；最后整合检测到的全部尖峰噪声并进行尖峰噪声簇判断，进而对尖峰噪声进行压制。经试验，本方法能够实现尖峰噪声的检测及压制，特别是对尖峰噪声簇尤为适用。与现有的尖峰噪声检测及压制方法相比，本方法在保证有用信号低误检的同时，还实现了尖峰噪声的低漏检，使得压噪后的地震数据质量大大提高，特别适合于受人员行走和风吹树动等干扰的地震数据去噪。

Description

一种低漏检的地震数据尖峰噪声检测及压制方法

技术领域：

本发明涉及一种地震数据处理方法，尤其是低漏检的地震数据尖峰噪声检测及压制方法，特别适用于受尖峰噪声簇干扰严重的地震数据去噪。

背景技术：

在地震勘探、地震监测和使用地震数据对火车、飞机和船舶等进行检测时，地震记录中常常会伴有强烈的尖峰噪声干扰，大大降低了地震数据质量。尖峰噪声常表现为单尖峰或尖峰噪声簇的形式。尖峰噪声簇由多个幅值和持续时间各异的单尖峰组成，这些单尖峰在时间轴上无规则紧密排列。尖峰噪声多以尖峰噪声簇的形式出现，其主要来源于人员行走、风吹树动和物体的突然放电等。此外，尖峰噪声在时域上表现为持续时间短、幅值与有用信号相当或远高于有用信号等特征。

目前，中值滤波器和均值滤波器是被广泛应用的尖峰噪声抑制方法，二者原理相近，对单尖峰噪声压制效果明显。但是，由于它们对数据中所有采样点都采用无差别平滑滤波策略，故而在用长滤波窗口滤除尖峰噪声簇时，会严重损坏有用信号。CN110542927A提出了一种变窗口加权地震数据尖峰噪声压制方法，该方法对不同持续时间的单尖峰噪声进行压制，但无法压制尖峰噪声簇。CN110542925A提出了一种基于峰值包络线的地震数据尖峰噪声识别及压制方法，但该方法主要适用于尖峰噪声与有用信号不重叠的情况，当二者混叠时则无效。

CN110542926A提出了一种地震数据尖峰噪声簇的自主检测和压制方法，该方法首先应用一个高阈值对尖峰噪声进行检测，然后进行尖峰噪声簇判断并实现噪声压制。该方法中的阈值设置值较高，从而确保了有用信号的低误检，但这也意味着尖峰噪声的高漏检。高阈值可以确保有用信号的低误检，低阈值则可以保证噪声的低漏检，而一个阈值难以兼顾有用信号的低误检和尖峰噪声的低漏检，因此，需要一种双阈值的地震数据尖峰噪声检测及压制方法，使误检和漏检得以兼顾，以解决现有技术在进行尖峰噪声检测时存在的高漏检问题。

发明内容：

本发明的目的就在于针对上述现有技术的不足，提供一种低漏检的地震数据尖峰噪声检测及压制方法。

本发明的技术要点是，首先应用高阈值对尖峰噪声进行检测，以保证有用信号的低误检，称此高阈值为全阈值；然后以已检测到的尖峰噪声为中心，取一定长度的窗口，应用低阈值对窗口内的数据进行二次检测，从而保证尖峰噪声的低漏检，称此低阈值为局部阈值；最后整合检测到的全部尖峰噪声并进行尖峰噪声簇判断，进而对尖峰噪声进行压制。

本发明的目的是通过以下技术方案实现的：

一种低漏检的地震数据尖峰噪声检测及压制方法，包括以下步骤：

a、读取一道包含尖峰噪声的地震数据X(n)，n＝1，2，...，N，N为数据总采样点数，采样率为Fs，观察时域波形并确定尖峰噪声的最短持续时间t_p；

b、应用STA/LTA方法对X(n)进行处理，记振幅的短时变化量与长时变化量的比值为R(n)，其中，特征函数、振幅的短时变化量、振幅的长时变化量和短时窗口长度L_s的计算方法参考CN110542926A，则0.01Fs≤L_s≤9Fs·t_p，取长时窗口长度L₁＝8L_s；

c、参考CN110542926A，计算全阈值

H＝Γ(R(n))+2γ(R(n)) (1)

其中，Γ()为计算均值的函数，γ()为计算标准差的函数；拾取所有满足R(n)＞H条件的采样点序号，将这些采样点序号构成的集合记为W；

d、拾取W中的所有采样点序号，对其进行从小到大排序，构成矩阵A，A为1行J列的矩阵，A＝[a₁，a₂，...，a_J]，J为W中采样点序号的总个数；令矩阵B＝[a₂，a₃，...，a_J，0]，计算B与A的差值

C＝B-A (2)

将C中的任意整数记为c_j，j＝1，2，...，J；拾取所有满足c_j≠1条件的j，并进行从小到大排序构成一个序列，记为(b₁，b₂，...，b_M)，M为满足c_j≠1条件的j的个数，将该序列中的任意项记为b_m，m＝1，2，...，M；将A分为M块，每个子块记为A_m，

其中，b₀＝0；将此过程称为拾取W中的采样点序号所对应的尖峰噪声的时间段；

e、令m＝1；

f、计算A_m中采样点序号的个数

确定用于计算A_m局部阈值的采样点序号的个数

L_m′＝10L_m (4)

g、拾取小于且紧邻

的采样点序号，这些采样点序号构成的集合为

拾取大于且紧邻

的采样点序号，这些采样点序号构成的集合为

计算P_m′和P_m″的并集

P_m＝P_m′∪P_m″ (7)

h、令T_m＝P_m-W，将T_m中每个采样点序号记为t_m，g，g＝1，2，...，G，G为T_m中所有采样点序号的总个数；拾取T_m中采样点序号的最小值和最大值，分别记为

和

i、拾取所有t_m，g对应的R(n)值，构成集合R_m′；

j、计算A_m局部阈值

H_m′＝Γ(R_m′)+1.5γ(R_m′) (8)

拾取所有满足R(v)＞H_m′条件的采样点序号，

将这些采样点序号构成的集合记为Z_m；

k、令m＝m+1；

l、重复执行步骤f至步骤k，直到m＝M+1时停止，得到M个集合Z₁，Z₂，...，Z_M；计算这M个集合的并集

W′＝Z₁∪Z₂∪...∪Z_M (9)

m、参考步骤c，拾取W′中的采样点序号所对应的尖峰噪声的时间段，记为矩阵D_q，q＝1，2，...，Q，Q为尖峰噪声的个数，D_q＝[d_q，d_q+1，d_q+2，...，d_q′]，d_q和d_q′分别为D_q的起始采样点序号和终止采样点序号，将由所有D_q构成的集合记为D，D＝{D₁，D₂，...，D_Q}；

n、对D中的所有D_q进行尖峰噪声簇判断，然后进行噪声压制，方法参考CN110542926A的步骤e至步骤h，记压制尖峰噪声后的数据为X(n)′。

有益效果：

经试验，本发明公开的一种低漏检的地震数据尖峰噪声检测及压制方法，能够实现尖峰噪声的检测及压制，特别是对尖峰噪声簇尤为适用。在利用可控震源或冲击震源进行地震勘探、进行地震监测或对火车、飞机和船舶进行检测时，所采集到的地震数据通常会受到人员行走、风吹树动和物体的突然放电等引起的尖峰噪声簇干扰。现有的尖峰噪声簇检测及压制方法只考虑到避免有用信号误检，而未考虑尖峰噪声的漏检问题，而本方法在保证有用信号低误检的同时，还可以有效解决漏检问题，实现了尖峰噪声的低漏检，使得压噪后的地震数据质量大大提高。

附图说明：

图1地震数据尖峰噪声压制结果比较

(a)不含噪声的理想可控震源输出信号s(n)

(b)s(n)中加入随机噪声和尖峰噪声后的地震数据X(n)

(c)中值滤波压制尖峰噪声的结果

(d)CN110542926A方法压制尖峰噪声的结果

(e)本方法压制尖峰噪声的结果

具体实施方式：

下面结合附图和实施例对本发明做进一步的详细说明。

低漏检的地震数据尖峰噪声检测及压制方法，包括以下步骤：

a、给定不含噪声的理想可控震源输出信号s(n)，n＝1，2，...，14001，将其作为理想信号，如图1(a)所示。对s(n)叠加随机噪声和尖峰噪声，其中尖峰噪声包括短持续时间的单尖峰噪声r₁(n)、尖峰噪声簇r₂(n)和长持续时间的单尖峰噪声r₃(n)，r₁(n)、r₂(n)和r₃(n)的持续时间分别为0.002s、0.241s和0.03s，如图1(b)中三个虚线框所示，得到含噪地震数据X(n)。采样率Fs＝1000。尖峰噪声的最短持续时间t_p＝0.002。此时X(n)信噪比SNR₀＝0.2108；

b、应用STA/LTA方法对X(n)进行处理，记振幅的短时变化量与长时变化量的比值为R(n)，其中，特征函数、振幅的短时变化量、振幅的长时变化量和短时窗口长度L_s的计算方法参考CN110542926A，则10≤L_s≤18，取长时窗口长度L₁＝8L_s；

c、参考CN110542926A，计算全阈值

H＝Γ(R(n))+2γ(R(n)) (1)

其中，Γ()为计算均值的函数，γ()为计算标准差的函数，本例中H＝2.284；拾取所有满足R(n)＞H条件的采样点序号，将这些采样点序号构成的集合记为W；

d、拾取W中的所有采样点序号，对其进行从小到大排序，构成矩阵A，A为1行98列的矩阵，A＝[a₁，a₂，...，a₉₈]；令矩阵B＝[a₂，a₃，...，a₉₈，0]，计算B与A的差值

C＝B-A (2)

将C中的任意整数记为c_j，j＝1，2，...，98；拾取所有满足c_j≠1条件的j，并进行从小到大排序构成一个序列，记为(b₁，b₂，...，b₆)，共有6个满足c_j≠1条件的j，将该序列中的任意项记为b_m，m＝1，2，...，6；将A分为6块，每个子块记为A_m，

其中，b₀＝0；将此过程称为拾取W中的采样点序号所对应的尖峰噪声的时间段；

e、令m＝1；

f、计算A_m中采样点序号的个数

确定用于计算A_m局部阈值的采样点序号的个数

L_m′＝10L_m (4)

g、拾取小于且紧邻

的采样点序号，这些采样点序号构成的集合为

拾取大于且紧邻

的采样点序号，这些采样点序号构成的集合为

计算P_m′和P_m″的并集

P_m＝P_m′∪P_m″ (7)

h、令T_m＝P_m-W，将T_m中每个采样点序号记为t_m，g，g＝1，2，...，G，G为T_m中所有采样点序号的总个数；拾取T_m中采样点序号的最小值和最大值，分别记为

和

i、拾取所有t_m，g对应的R(n)值，构成集合R_m′；

j、计算A_m局部阈值

H_m′＝Γ(R_m′)+1.5γ(R_m′) (8)

拾取所有满足R(v)＞H_m′条件的采样点序号，

将这些采样点序号构成的集合记为Z_m；

k、令m＝m+1；

l、重复执行步骤f至步骤k，直到m＝7时停止，得到6个集合Z₁，Z₂，...，Z₆；计算这6个集合的并集

W′＝Z₁∪Z₂∪...∪Z₆ (9)

m、参考步骤c，拾取W′中的采样点序号所对应的尖峰噪声的时间段，记为矩阵D_q，q＝1，2，...，9，共9个尖峰噪声，D_q＝[d_q，d_q+1，d_q+2，...，d_q′]，d_q和d_q′分别为D_q的起始采样点序号和终止采样点序号，将由所有D_q构成的集合记为D，D＝{D₁，D₂，...，D₉}；

n、对D中的所有D_q进行尖峰噪声簇判断，然后进行噪声压制，方法参考CN110542926A的步骤e至步骤h，记压制尖峰噪声后的数据为X(n)′，如图1(e)所示。

此时，经本发明提出的方法处理后的地震数据信噪比SNR′＝6.4708，尖峰噪声能量压制百分比Per′＝98.59％，数据质量大大提高，此外，可以看到r₁(n)和r₃(n)完全被压制，r₂(n)处理后振幅仅为原来的1/4，如图1(e)虚线框。图1(c)为应用中值滤波器处理后的结果，其中r₁(n)完全被压制，r₂(n)的振幅约为原来的2/5，r₃(n)基本没有改变，此时信噪比SNR″＝-0.3112，尖峰噪声能量压制百分比Per″＝81.50％，虽然大部分尖峰噪声的能量得到压制，但由于中值滤波器对有用信号的损坏严重，使得数据质量差于去噪前的数据。图1(d)为应用CN110542926A方法压制尖峰噪声后的结果，其中r₁(n)完全被压制，r₂(n)稍有削减，r₃(n)的振幅约为原来的3/5，此时信噪比SNR′″＝0.3446，尖峰噪声能量压制百分比Per′″＝35.08％，尖峰噪声漏检严重，从而使得噪声压制效果差，数据质量改善不明显。而本专利提出的方法可以有效检测尖峰噪声，实现尖峰噪声的低漏检，对于长持续时间有用信号中混叠的尖峰噪声的检测和压制效果明显，尤其适用于尖峰噪声簇的检测和压制。

Claims (1)

1.一种低漏检的地震数据尖峰噪声检测及压制方法，其特征在于，包括以下步骤：

a、读取一道包含尖峰噪声的地震数据X(n)，n＝1，2，...，N，N为数据总采样点数，采样率为Fs，观察时域波形并确定尖峰噪声的最短持续时间t_p；

b、应用STA/LTA方法对X(n)进行处理：

b1)、计算地震信号X(n)的特征函数

b2)、利用递归算法计算地震信号振幅的短时变化量

长时变化量

式中L_l和L_s为整数，且0.01Fs≤L_s≤9Fs·t_p，取长时窗口长度L_l＝8L_s；

b3)、记振幅的短时变化量与长时变化量的比值为R(n)，即

c、为了有效检测出尖峰噪声，采用阈值检测的方式；针对检测中出现的尖峰噪声漏检现象，需要进行二次阈值检测，为了区分两种阈值，将前者定义为全阈值，后者定义为局部阈值；

计算尖峰噪声全阈值

H＝Γ(R(n))+2γ(R(n)) (5)

其中，Γ()为计算均值的函数，γ()为计算标准差的函数；拾取所有满足R(n)＞H条件的采样点序号，将这些采样点序号构成的集合记为W；

d、拾取W中的所有采样点序号，对其进行从小到大排序，构成矩阵A，A为1行J列的矩阵，A＝[a₁，a₂，...，a_J]，J为W中采样点序号的总个数；令矩阵B＝[a₂，a₃，...，a_J，0]，计算B与A的差值

C＝B-A (6)

将C中的任意整数记为c_j，j＝1，2，...，J；拾取所有满足c_j≠1条件的j，并进行从小到大排序构成一个序列，记为(b₁，b₂，...，b_M)，M为满足c_j≠1条件的j的个数，将该序列中的任意项记为b_m，m＝1，2，...，M；将A分为M块，每个子块记为Am，

其中，b₀＝0；将此过程称为拾取W中的采样点序号所对应的尖峰噪声的时间段；

e、令m＝1；

f、计算A_m中采样点序号的个数

确定用于计算A_m局部阈值的采样点序号的个数

L_m′＝10L_m (8)

g、拾取小于且紧邻

的采样点序号，这些采样点序号构成的集合为

拾取大于且紧邻

的采样点序号，这些采样点序号构成的集合为

计算P_m′和P_m″的并集

P_m＝P_m′∪P_m″ (11)

h、令T_m＝P_m-W，将T_m中每个采样点序号记为t_m，g，g＝1，2，...，G，G为T_m中所有采样点序号的总个数；拾取T_m中采样点序号的最小值和最大值，分别记为

和

i、拾取所有t_m，g对应的R(n)值，构成集合R_m′；

j、计算A_m局部阈值

H_m′＝Γ(R_m′)+1.5γ(R_m′) (12)

拾取所有满足R(v)＞H_m′条件的采样点序号，

将这些采样点序号构成的集合记为Z_m；

k、令m＝m+1；

l、重复执行步骤f至步骤k，直到m＝M+1时停止，得到M个集合Z₁，Z₂，...，Z_M；计算这M个集合的并集

W′＝Z₁∪Z₂∪...∪Z_M (13)

m、参考步骤c，拾取W′中的采样点序号所对应的尖峰噪声的时间段，记为矩阵D_q，q＝1，2，...，Q，Q为尖峰噪声的个数，D_q＝[d_q，d_q+1，d_q+2，...，d_q′]，d_q和d_q′分别为D_q的起始采样点序号和终止采样点序号，将由所有D_q构成的集合记为D，D＝{D₁，D₂，...，D_Q}；

n、对D中的所有D_q进行尖峰噪声簇判断，然后进行噪声压制：

n1)、对每一个D_q判断其长度是否满足条件

d′_q-d_q+1≤F_S (14)

若满足条件，该区间被认为存在尖峰噪声，否则视为不存在尖峰噪声；

n2)、经过步骤n1)判断得出尖峰噪声的连续采样点区间D′＝{D′₁，D′₂，...，D′_K}，K为尖峰噪声个数，

为第k个尖峰噪声段，

和

为该尖峰噪声所在连续采样点区间的起始采样点序号和终止采样点序号，k∈[1，K]；定义第k个尖峰噪声的持续采样点数的一半为

为向上取整运算；

n3)、对于每一个D′_k的定义两个权值

式中

n4)、对于每一个D′_k的采样点，在其两侧各取L_k窗长的信号

令

式中Θ()为计算中值的函数；X′(n)为压制尖峰噪声后的信号。

Images