CN109188510A - 一种高精度自动识别和压制地震数据单频干扰的方法 - Google Patents

Abstract

本发明提供一种高精度自动识别和压制地震数据单频干扰的方法，包括输入原始地震数据，对深部地震信号进行求导运算得到地震信号导数的频谱，并对其进行叠加处理；根据定义的单频干扰判别准则在频谱上识别单频干扰，并将其频率作为初始频率；重构单频干扰信号，对构建的目标函数求导得到单频干扰振幅；计算由目标函数转换的归一化互相关函数值，将互相关函数的极大值对应的频率确定为单频干扰的最佳频率；将单频干扰的频率及其振幅代入正余弦函数公式，得到重构的单频干扰信号,从原始地震数据中减去重构的单频干扰信号，完成单频干扰的压制工作。本发明提供的方法能够自动识别和压制地震数据中的单频干扰，并能取得保真去噪的效果。

Description

一种高精度自动识别和压制地震数据单频干扰的方法

技术领域

本发明涉及一种高精度自动识别和压制地震数据单频干扰波的方法，尤其适用于消除地震数据中多个不同频率的单频干扰信号。

背景技术

单频干扰(又称为谐波干扰、工频干扰)是地震勘探中一种典型的干扰波，具有振幅固定、频率单一的特点。单频干扰产生于野外地震数据采集阶段，其产生机制：(1)当检波器附近存在高压输电线时，检波器会受到工业交流电产生的50Hz电磁场干扰，以及100Hz或150Hz等频率的伴随电磁场干扰；(2)当检波器附近存在频率固定的震源时(如停滞态汽车、匀速转动的马达、工作的发电机等)，检波器会受到单一频率的机械震动干扰；(3)当连接检波器的电缆漏电时，检波器也会记录下50Hz的电流信号干扰。与其它干扰波相比，单频干扰分布范围广，一般贯穿于整个地震道记录，而且深部数据中单频干扰的振幅明显大于有效反射信号的振幅，有效反射信号往往被单频干扰淹没，从而严重降低地震数据的信噪比。因此，消除单频干扰波已成为地震数据处理中的重要环节。

目前，国内外地球物理学者已研究出多种单频干扰波识别和压制方法，如陷频滤波法、带通滤波均方根能量单频识别法、归一化目标函数互相关系数识别法，余弦函数重构单频干扰压制等方法。当地震数据中只包含整50Hz频率的单频干扰时，现有方法能够对其压制。但是，由于实际地震数据中的单频干扰频率不是整数，而且数据中也不只包含一个50Hz的单频干扰，同时还存在多个50Hz的伴随单频干扰，如99.6Hz、149.8Hz等。因此，传统方法不能消除地震数据中的所有单频干扰。目前，大多数方法在单频干扰的识别问题上，一般是通过人为预先给定单频干扰的初始频率，然后采用相关方法逼近实际数据中的单频干扰频率，因而还不具有自动获得初始频率的功能。由于实际地震数据中除50Hz频率的单频干扰外，其它伴随单频干扰为随机分布，事先很难给定。因此人为给定的方法存在一定的局限性，在数据处理中经常出现由于初始频率不准或没有囊括所有单频干扰，导致去噪后的数据中仍有残留的单频干扰。此外，现有方法还存在计算量大、计算效率低的问题。因此，迫切需要研究出一种既快又准的单频干扰波识别和压制方法，以适应当今海量地震数据的去噪处理工作。

(1)陷波器压制单频干扰方法

早期对单频干扰的压制是从野外采集环节入手，通过地震仪器中的陷波器直接过滤这种50Hz的电磁场干扰(黄绪德等，1977；宋祈真等，1979；孙文森等，1990)。而陷波器一旦设置好后，只能针对固定频率的干扰信号进行压制。由于实际地震道中一般存在多个频率的单频干扰波，而且单频干扰频率不是整50Hz，因此这种过滤固定频率信号的压制方法，应对问题不够灵活，压制频率不够准确，难以获得理想的降噪效果。

(2)频率滤波法压制单频干扰

在1980至1990年代，随着计算机技术和数字信号处理技术的深入发展，单频干扰的压制由室外转移到室内。在此期间，一般采用数字信号处理技术中频率域带通滤波法压制单频干扰，如基于相位滤波理论，通过设计50Hz的滤波器，滤除地震信号中50Hz的单频干扰(郭其贵等，1988)；在频率域中采用中值滤波来过滤单频干扰(凌云等，1992)。上述方法均借用了数字信号处理领域相对成熟的技术，在实践中虽然取得了一定的噪声压制效果，但仍然存在较为严重的残存干扰。而且，这些方法在压制单频干扰的同时，也伤及了单频干扰附近频带的有效信号。

(3)重构单频干扰模型压制方法

近十年来，单频干扰压制方法的发展更为多样化，如：①通过构造频率整50Hz的单频干扰模型，计算构造信号和原始地震信号的互相关来压制50Hz单频干扰，显然这种方法只能针对50Hz的单频干扰进行压制，应用对象单一(李文艳等，2001)；②采用共轭梯度法估算单频干扰的频率、相位、振幅，并由该三参数重构单频干扰(刘洋等，2003；段云卿等，2006)；③余弦函数重构单频干扰的压制方法(吴小培等，2003；胡伟等，2005；高少武等，2008；陈可洋等，2014；陈可洋等，2017)。余弦函数构造单频干扰方法从频率、相位、振幅三个参数逐渐逼近单频干扰，在频率逼近能达到0.01Hz的精度，与其它方法相比，该方法计算精度高。但由于该方法需要对频率、相位参数进行扫描，因此计算量大，计算效率低，目前的余弦函数重构单频干扰方法还不能应对海量地震数据的处理。因此，改进现有技术方法，提高计算效率，以应对海量地震数据的降噪处理。

(4)单频干扰识别方法

关于单频干扰识别的方法并不多，目前主要有：带通滤波均方根能量单频干扰识别方法(李慧等，2012)，基于归一化互相关函数识别方法(陈可洋等，2014)、线性调频谱法识别单频干扰(高少武等，2010)。前两种方法是在时间域利用单频干扰的振幅特性来识别单频干扰，当地震数据中无其它干扰波或其它干扰波能量较弱时，能较准确识别出单频干扰，但当地震数据中存在其它干扰波且能量较强时，如面波、异常振幅等干扰波，采用这两种方法识别单频干扰就比较困难。而线性调频谱法识别单频干扰是在频率域通过频率扫描查找最大谱对应的频率来识别单频干扰，该方法可以扫描得到单频干扰的频率，但是用于频率扫描的初始值则需要人为给定，而且为了达到一定精度，需要进行大量的计算，因而该方法计算量大、效率低。为了应对海量地震数据的处理，发明扫描速度快、计算精度高的单频干扰自动识别和压制方法是适应现代地震勘探技术的需要。

发明内容

本发明的目的是在于提供一种高精度自动识别和压制地震数据单频干扰的方法，以消除地震数据中单频干扰对有效反射信号的影响，为后续地震数据处理工作奠定基础。该方法既不误伤有效反射信号，也不漏杀单频干扰，能够自动地进行保真去噪处理。

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

一种高精度自动识别和压制地震数据单频干扰的方法，包括如下步骤：

S1，输入原始地震数据，选取2s以下数据进行n阶求导，利用傅氏变换将求导信号转换到频率域，得到原来单频干扰振幅4倍以上的频谱；由于在直接利用原始地震数据进行傅氏变换获得的频谱中，低频段单频干扰的能量强于有效信号，而高频段单频干扰能量则不明显。因此，为了突出单频干扰能量，利于单频干扰在频谱上的识别，采用了一种首先计算所选地震信号的n阶导数，然后再对导数信号进行傅氏变换的方法。该方法可以获得单频干扰与有效信号具有较大能量差异的频谱。因此，在频谱上根据可能量差异定义判别准则，识别单频干扰。

S2，对得到的频谱进行叠加处理，消除锯齿状振幅起伏现象对单频干扰识别的影响；

S3，根据单频干扰的凸起频率振幅特征定义单频干扰判别准则，根据判别准则在频谱上一次性识别出所有单频干扰，并将单频干扰对应的频率作为后续计算的初始频率；

S4，利用正余弦函数逼近法重构单频干扰信号，将原始地震信号与重构单频干扰信号差的平方求和定义为目标函数，对目标函数求导得到单频干扰的振幅；

S5，将定义的目标函数转换为原始地震信号与重构单频干扰的归一化互相关函数，采用频率折半逼近算法进行频率扫描，得到多个互相关函数的值，将其极大值对应的频率视为单频干扰的最佳频率；

S6，将单频干扰频率及用其计算的振幅代入正余弦函数公式，得到重构的单频干扰信号；

S7，从原始地震数据中减去重构的单频干扰信号，获得无单频干扰的地震数据，完成消除单频干扰处理工作。

优选地，根据下述公式计算地震信号的n阶导数，并进行傅氏变换：

已知地震信号是离散采样序列，其求导公式如下：

公式(2-1)中X^T表示原始地震信号，(X^T)′表示一次求导信号，i表示采样点号，N表示采样点总数。同理，把一阶求导信号作为原始信号再次求导，可得到二阶求导函数。依此类推，可以得到n阶求导函数，记为：(X^T)⁽ⁿ⁾。

计算地震信号一阶导数频谱的傅氏变换公式如下：

F[(X^T)′]＝iωF(ω) (2-2)

其中，F[·]表示傅氏变换，ω为数字频率，导函数经过傅氏变换得到的各频率的幅度函数为F(ω)，它是原函数傅氏变换幅值的iω倍。

优选地，根据下述变换公式得到地震数据导数的频谱：

若函数X^T的n阶导数(X^T)⁽ⁿ⁾的傅氏变换存在，则二者关系式如下：

F[(X^T)⁽ⁿ⁾]＝(iω)ⁿF(ω) (2-3)

从式(2-3)可以看出，n阶导函数经过傅氏变换得到的幅度函数是原函数傅氏变换幅度的(iω)ⁿ倍。可见，单频干扰的能量得到成倍提高，利于单频干扰的识别。

优选地，对地震数据导数的频谱进行叠加处理，获得振幅平滑后的频谱，以消除锯齿状频谱现象对单频干扰识别的影响。频谱叠加公式如下：

对选取的原始地震数据求导后，得到采样点总数为K的地震信号，从该信号中分段连续截取N个样点，计算得到第i段的频谱：

F_i(n) 1≤n≤N且N≤K (2-4)

其中，i表示截取的起始样点，其取值范围为1≤i≤K+1-N，F₁(n)表示[1，N]范围内的样点计算得到的频谱，F₅(n)表示[5，4+N]范围内的样点计算得到频谱。

将所有分段频谱按照下述公式计算，可获得振幅平滑后的频谱：

优选地，根据单频干扰频谱特征，定义单频干扰判别准则，根据判别准则在频谱上识别单频干扰，其判别准则定义公式如下：

已知采样率为F，样点数为N的一段地震信号经傅氏变换后，将获得的N个复数点取模，即得到各复数点对应的幅值(振幅)，记为：

A(n) 1≤n≤N (2-6)

如果振幅A(n)满足下面的不等式：

A(n-1)＜A(n)，且A(n)＞A(n+1) 1＜n＜N (2-7)

则将其对应的频率定义为“凸起频率”，然后将每个“凸起频率”对应的振幅标记为P(i)，“凸起频率”总数记为H。

如果P(i)满足不等式(2-8)，则其频率对应的信号被判定为“疑似干扰”

P(i-1)＜P(i)，P(i)＞P(i+1) 1＜i＜H (2-8)

若第j个“疑似干扰”频率的振幅P(j)能使不等式(2-9)成立，则该“疑似干扰”即被判定为单频干扰，判别公式如下：

不等式(2-9)表示“疑似干扰”对应的振幅的4倍与周边4个“凸起频率”振幅和的比值，如果比值大于一定的阈值M，则“疑似干扰”被确定为单频干扰，并记录单频干扰的频率，为后续频率扫描提供初始参数。

优选地，在步骤S6中，利用正余弦函数扫描法重构单频干扰信号，将原始地震信号与重构单频干扰信号两者差的平方求和作为目标函数，对目标函数求导得到单频干扰的振幅，其计算公式如下：

y_i＝A cos 2πfiΔt+B sin 2πfiΔt (2-10)

其中，y_i表示重构的单频干扰，A表示余弦信号的振幅，B表示正弦信号的振幅，f表示单频干扰的频率，i为采样点号，Δt为地震数据采样时间间隔。

公式(2-10)的向量形式如下：

y＝Bs^T+Aa^T (2-11)

其中，s^T表示单频干扰的正弦序列，a^T表示单频干扰的余弦序列，两者展开式如下：

由原始采样序列估算出单频干扰序列，需定义的目标函数如下：

其中，x_i为原始地震信号，Q表示原始采样序列减去重构的单频干扰序列后的振幅能量值，X^T为原始采样序列。

利用公式(2-13)，对Q关于各变量求导，可获得Q的极小值，将能使Q为极小值的重构序列视为理想的单频干扰序列。具体计算如下：

假设频率已知，对Q关于正弦振幅B和余弦振幅A求导，并使导数为0，可以获得一个二元一次方程组：

求解方程组(2-14)，可以获得A和B的值，其计算公式：

如果确定了单频干扰的准确频率，则可由公式(2-15)计算得到重构单频干扰的振幅A和B，将频率、振幅代入公式(2-10)，即可得到重构的单频干扰。

优选地，将定义的目标函数(2-13)转换成向量形式，得到原始地震采样序列与重构单频干扰采样序列的归一化互相关函数，由归一化互相关函数的取值确定单频干扰的最佳频率。其公式如下：

将公式(2-13)改写成向量形式，可得到原始地震采样序列与重构单频干扰采样序列的归一化互相关函数，其公式为：

由公式(2-16)，不同频率f可以计算出不同的互相关函数值，即Corr(f)，将使得互相关函数值为极大值的频率视为单频干扰的最佳频率。

优选地，将公式(2-9)确定的单频干扰频率作为频率扫描的初始值，采用频率折半逼近算法计算公式(2-16)中的互相关函数值，函数取极大值的频率即为单频干扰的最佳频率。与传统频率扫描算法相比，频率折半逼近算法可以快速得到单频干扰的频率，因此计算效率高。

①线性频率扫描算法：

如给定扫描初始频率f₀，第一级频率扫描的样点如下：

f_k＝f₀+0.1k-10≤k≤10 (2-17)

在第一级扫描样点中找到最优频率f₁，将f₁作为第二级扫描的初始频率，则第二级频率扫描的样点如下：

f_k＝f₁+0.01k-10≤k≤10 (2-18)

以此类推，最终可以获得满足精度要求的单频干扰频率f。线性频率扫描算法采取对单频干扰的准确频率分级逼近的策略，每一级精度提升为上一级精度的0.1倍，例如：在49.5Hz到50.5Hz之间进行频率扫描，频率精度可以达到0.001Hz。但是频率精度每提高一级至少需要计算33个频率样点。因此，该方法计算量大，效率低。

②频率折半逼近算法

将识别出的单频干扰频率作为频率扫描的初始值f₀，初始频率扫描范围为Δf，则第一级频率扫描样点公式为：

f_k＝f₀+Δfk k＝[-1，0，1] (2-19)

其中，f_k为扫描出的频率样点，k为扫描的点数。

第一级频率扫描共计算了3个频率样点，其精度能达到Δf。

从第一级的三个频率样点中找出能使Corr(f)为最大时的频率，记为f₁。然后计算第二级的频率样点f₂，f₃，第二级频率扫描样点公式为：

第二级频率扫描实际计算了2个频率样点，该级频率精度达到了Δf/2，将f₁，f₂，f₃列为第二级频率样点。

从第二级的三个频率样点中找出Corr(f)为最大时的频率，同样记为f₁。然后计算第三级的频率样点f₂，f₃，第三级频率扫描样点公式为：

第三级频率扫描实际计算了2个频率样点，该级频率精度达到(1/2)²Δf，将f₁，f₂，f₃划归为第三级频率样点。

以此类推，将上一级Corr(f)为最大时的频率均记为f₁，然后计算第n级的频率样点f₂，f₃，第n级频率扫描样点公式为：

第n级频率扫描同样是计算了2个频率样点，而且频率精度达到了(1/2)^(n-1)ΔfHz，将f₁，f₂，f₃划归为第n级频率样点。当(1/2)^(n-1)Δf达到了预期的精度时，f₁即为逼近的单频干扰频率。

例如：频谱识别出的单频干扰为50Hz，将其作为初始频率f₀，单频干扰的初始扫描范围Δf为0.5Hz，只需计算3+9+9＝21个频率样点，单频干扰的精度便能达到0.00097656Hz(约0.001Hz)，。可见，频率折半逼近法较线性频率扫描方法减少了12个频率样点的计算量。因此，在相同的计算精度前提下，频率折半逼近法计算效率高。

优选地，将扫描得到的单频干扰频率和用其重新计算得到的振幅代入公式(2-10)后，可获得重构的单频干扰信号，从原始地震数据中减去重构的单频干扰信号，即得到无单频干扰的地震数据。其计算公式如下：

S＝X-Y (2-23)

其中，X为原始地震数据，Y为重构的单频干扰信号，S为去除单频干扰后的地震数据。

实际数据采集过程中，检波器是每隔一段时间记录一次振幅数据(如采样率为1ms)。因此，获得的地震数据是一系列的离散振幅值，故公式(2-23)中的三个参数可以展开成如下向量的形式：

S^T＝(s₁，s₂，…，s_i，…，s_k)

X^T＝(x₁，x₂，…，x_i，…，x_k) (2-24)

Y^T＝(y₁，y₂，…，y_i，…y_k)

公式(2-24)中X^T为原始地震信号采样序列，S^T为有效信号采样序列，Y^T为单频干扰采样序列，i为采样点编号，k为采样点总数，T表示向量的转置。用正余弦函数逼近法重构单频干扰就是在已知原始采样序列X^T的情况下，估算出干扰序列Y^T。

综上所述，本发明提供的方法能够自动识别和压制地震数据中的单频干扰，并能取得保真去噪的效果。

附图说明

图1是本发明提供的高精度自动识别和压制地震数据单频干扰方法的流程图；

图2(a)-图2(d)为同一时间段地震信号经过求导得到的频谱及进行叠加处理后的频谱；

图3为利用图2d频谱，按照单频干扰的频谱判别准则识别出的单频干扰；

图4(a)-图4(b)压制单频干扰后的频谱；

图5(a)-图5(c)是利用本发明提供的方法消除单频干扰前后的模型数据及残差图；

图6(a)-图6(c)是利用本发明提供的方法去除实际地震数据中单频干扰前后的结果及残差图.

具体实施方式

下面结合附图对本发明作进一步详述，以下实施例只是描述性的，不是限定性的，不能以此限定本发明的保护范围。

如图1所示，图1是本发明提供的高精度自动识别和压制地震数据单频干扰方法的流程图。

一种高精度自动识别和压制地震数据单频干扰的方法，包括如下步骤：

S1，输入原始地震数据，选取深层时间段地震信号，并进行n阶求导；利用傅氏变换将时间域的求导信号转换到频率域，获得单频干扰振幅有较大提高的频谱。

已知地震信号是离散采样序列，其求导公式如下：

公式(3-1)中X^T表示原始地震信号，(X^T)′表示一次求导信号，i表示采样点号，N表示采样点总数。同理，把一阶求导信号作为原始信号再次求导，可得到二阶求导函数。依此类推，可以得到n阶求导函数，记为：(X^T)⁽ⁿ⁾。其傅氏变换公式如下：

F[(X^T)⁽ⁿ⁾]＝(iω)ⁿF(ω) (3-2)

其中，F[·]傅氏变换符号，ω为数字频率，F(ω)为原始地震数据的频谱幅度函数。由式(3-2)中可知，对求导函数进行傅氏变换得到的幅度函数值是直接对原始数据进行傅氏变换得到的幅度函数值的(iω)ⁿ倍。因此，与原函数频谱相比，导函数频谱中单频干扰的能量会成倍提高。通过实验对地震信号经过1至3阶求导运算后，频谱上单频干扰的幅度会显著提升(如图2所示)。

图2展示了同一时间段地震信号分别经过一阶和二阶导数计算后，得到的频谱。图2a为原始地震信号的频谱，从该图上只能初步看出50Hz频率周围存在疑似单频干扰。图2b为原始地震信号一阶导数的频谱，与图2a相比，频率在50Hz、150Hz和250Hz的幅度较原来有了一定程度的提升。图2c为原始地震信号二阶导数的频谱。与图2b相比，疑似单频干扰频率的幅度又有很大提升，而且高频段疑似单频干扰的幅度较低频提升更大。由图2可见，随着对地震信号求导阶数增大，经傅氏变换获得的频谱中疑似单频干扰能量则成倍提升，而其有效信号的频谱幅度则变化不明显。因此，地震信号求导后的频谱更利于识别单频干扰。

S2，对得到的频谱进行叠加处理，消除“锯齿状”频谱现象对单频干扰识别的影响。

由地震信号生成的原始频谱，相邻频率样点的振幅大小参差不齐，“锯齿状”现象较为严重，在一定程度上会影响单频干扰的判别。因此，需要进行频谱平滑处理，以消除频谱的“锯齿状”现象(锯齿状指的是：不同频率的振幅大小不同，在频谱上会出现高低起伏现象，这种起伏现象即称为“锯齿状”现象)。

对选取的原始地震数据求导后得到采样点总数为K的地震信号，从该信号中分段连续截取N个样点，计算得到第i段的频谱：

F_i(n) 1≤n≤N且N≤K (3-3)

其中，i表示截取的起始样点，其取值范围为1≤i≤K+1-N，F₁(n)表示[1，N]范围内的样点计算得到的频谱，F₅(n)表示[5，4+N]范围内的样点计算得到频谱。

将所有分段频谱按照下述公式计算，可获得振幅平滑后的频谱：

具体实施如下：

由于地震数据中单频干扰的频率和振幅比较固定，因此选取任意时段的地震信号进行傅氏变换，其频谱特征基本不变，而单频干扰以外的其它信号的频谱则会有一定的变化。因此，通过对不同时间段信号的频谱进行叠加求和，再取其均值，从而得到一个新的频谱。在新得到的频谱上，单频干扰的幅度相对叠加处理前，基本保持不变，而其它频率成分的频谱则被平滑，“锯齿状”频谱现象会有明显改善。因此，频谱叠加可以降低其它频率成分被误判为单频干扰的风险。

例如：记录长度为6s的地震记录，选取2s-4s的地震信号进行频谱判别，对节选出的2s长记录按照1s时窗计算频谱，每计算一次频谱，时窗向后滑动1个采样点(也可以设置推进多个采样点)，待计算出所有时窗频谱后，对这些频谱求和再取均值，从而形成一个新的频谱。

图2c和图2d展示了对频谱进行叠加处理前后的结果。图2c为未进行叠加处理的频谱，从该图中可以看到，“锯齿状”频谱现象比较明显。图2d为叠加处理后的频谱，可以看到“锯齿状”现象基本消除，而单频干扰的频谱则基本未变。

S3，根据单频干扰“凸起频率”特征定义单频干扰判别准则，根据判别准则在频谱上一次性识别出所有单频干扰，并将单频干扰对应的频率作为后续计算的初始频率。具体实施如下：

采样率为F，样点数为N的一段地震信号经傅氏变换后，获得的N个复数点取模，可以得到各复数点对应的振幅，记为：

A(n) 1≤n≤N (3-5)

对获得的振幅参数按不等式(3-6)进行判断，如果振幅参数满足该不等式，则将其对应的频率定义为“凸起频率”。

A(n-1)＜A(n)，A(n)＞A(n+1) 1＜n＜N (3-6)

将每个“凸起频率”对应的振幅参数重新标记为P(i)，“凸起频率”总数记为H。如：P(1)表示第1个“凸起频率”振幅参，P(10)表示第10个“凸起频率”振幅。按照不等式(3-7)对P(i)进行判断，如果P(i)满足该判别式，则其对应的频率被判定为“疑似干扰”。

P(i-1)＜P(i)，P(i)＞P(i+1) 1＜i＜H (3-7)

如果第j个“凸起频率”被判定为“疑似干扰”，则继续进行排查。如果仍然满足不等式(3-8)，则该“疑似干扰”即为单频干扰。

不等式(3-8)表示“疑似干扰”对应的振幅与周边4个(也可以设定为6个)“凸起频率”的振幅的均值做比，如果比值大于一定的阈值M，则“疑似干扰”被判定为单频干扰(阈值M人为设定，取2或3即可)。

图3展示了利用频谱判别法识别出的单频干扰。根据本发明定义的单频干扰判别准则，经过“凸起频率”、“疑似干扰”的排查，最后确定的单频干扰的初始频率为50Hz、150Hz和250Hz，从而为后续计算提供了初始参数。

S4，利用正余弦函数逼近法重构单频干扰信号，将原始地震信号与重构单频干扰信号差的平方求和定义为目标函数，对目标函数求导得到单频干扰的振幅，具体公式如下：

当原始地震记录中含有单频干扰时，可以认为单频干扰是一个谐波信号，该信号在开始采样时具有一定的初始相位，在整个采样时间里频率和振幅保持不变，用正余弦函数表示如下：

y_i＝A cos 2πfiΔt+B sin 2πfiΔt (3-9)

其中，y_i表示重构的单频干扰，A表示余弦信号的振幅，B表示正弦信号的振幅，f表示单频干扰的频率，i为采样点号，Δt为地震数据采样时间间隔。

公式(3-9)的向量形式如下：

y＝Bs^T+Aa^T (3-10)

其中，s^T表示单频干扰的正弦序列，a^T表示单频干扰的余弦序列，两者展开式如下：

由原始采样序列估算出单频干扰序列，定义的目标函数如下：

其中，x_i为原始地震信号，Q表示原始采样序列减去重构的单频干扰序列后的振幅能量值，X^T为原始采样序列。

利用公式(3-12)，对Q关于各变量求导，可获得Q的极小值，将能使Q为极小值的重构序列视为理想的单频干扰序列。具体计算如下：

假设频率已知，对Q关于正弦振幅B和余弦振幅A求导，并使导数为0，可以获得一个二元一次方程组：

求解方程组(3-13)，可以获得A和B的值，其计算公式如下：

如果确定了单频干扰的最佳频率，则由公式(3-14)计算出重构单频干扰的振幅A和B，将频率、振幅代入公式(3-9)，即可得到重构的单频干扰。

S5，将目标函数转换为原始地震信号与重构单频干扰的互相关函数，采用频率折半逼近算法进行频率扫描，得到多个互相关函数值，将函数的极大值对应的频率确定为单频干扰的准确频率。具体公式如下：

将公式(3-12)改写成向量形式，可得到原始地震采样序列与重构单频干扰采样序列的归一化互相关函数，其公式为：

由公式(3-8)确定的单频干扰频率作为频率扫描的初始值，采用频率折半逼近算法由公式(3-15)计算互相关函数值Corr(f)，根据Corr(f)的极值确定单频干扰的频率，其计算公式如下：

给定单频干扰频率扫描的初始值f₀和初始频率扫描范围Δf，则第一级频率扫描样点公式如下：

f_k＝f₀+Δfk k＝[-1，0，1] (3-16)

公式(3-16)表明单频干扰频率存在于以f₀为中心，正负Δf范围之内。第一级频率扫描共计算了3个频率样点，其精度可达到Δf。

从第一级的三个频率样点中，找到Corr(f)为最大时的频率，记为f₁，按照公式(3-17)计算第二级的频率样点f₂，f₃。将f₁，f₂，f₃划归为第二级频率样点，第二级频率扫描实际计算了2个频率样点，该级频率精度达到Δf/2。

再从第二级的三个频率样点中，找到Corr(f)为最大时的频率，记为f₁。按照公式(3-18)计算第三级的频率样点f₂，f₃将f₁，f₂，f₃划归为第三级频率样点，第三级频率扫描实际计算了2个频率样点，该级频率精度达到(1/2)²Δf。

以此类推，将上一级Corr(f)中最大值对应的频率记为f₁，本级的频率样点f₂，f₃按照下面通式计算，

式中n表示扫描等级。将f₁，f₂，f₃划归为第n级频率样点，第n级频率扫描实际计算了2个频率样点，该级频率精度达到(1/2)^(n-1)ΔfHz，当频率精度达到了预期的精度时，则频率f₁被视为单频干扰的最佳频率。

S6，将得到的单频干扰频率及用其计算的振幅代入正余弦函数重构公式，得到重构的单频干扰信号。

S7，从原始地震数据中减去重构的单频干扰信号，获得无单频干扰的地震数据，完成消除单频干扰的处理工作。

当原始地震信号中含有单频干扰时，该地震信号可以看作有效信号和单频干扰信号的叠加，如下式所示：

X＝S+Y′ (3-20)

其中，X为原始地震信号，S为地震信号的有效成分，Y′为地震信号的单频干扰成分。在公式(3-20)中，X是已知的，S和Y′未知。由于单频干扰波具有周期性的特点，因此可通过正余弦函数方法重构出单频干扰Y，然后从原始地震记录中减去重构的单频干扰Y，即可获得无单频干扰的地震数据。若重构的单频干扰Y与实际地震信号中的单频干扰Y′极度逼近，则压制单频干扰后，可以得到比较理想的效果。具体实施如下：

在被工业交流电干扰污染的地震数据中，一般除存在一个频率大约50Hz的单频干扰之外，还会存在几个伴随频率的单频干扰。因此，要考虑多个单频干扰同时存在的压制策略。

利用原始地震信号先估算第一个频率的单频干扰，再估算第二个频率的单频干扰，依次类推，直到估算出所有频率的单频干扰后，将各个频率的单频干扰相加求和得到整体干扰，最后从原始地震信号减去整体干扰即完成压制全部单频干扰的处理工作。

Claims (11)

1.一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，包括如下步骤：

S1，输入原始地震数据，选取2s以下数据进行n阶求导，利用傅氏变换将求导信号转换到频率域，得到原来单频干扰振幅4倍以上的频谱；

S2，对得到的频谱进行叠加处理，消除锯齿状振幅起伏现象对单频干扰识别的影响；

S3，根据单频干扰的凸起频率振幅特征定义单频干扰判别准则，根据判别准则在频谱上一次性识别出所有单频干扰，并将单频干扰对应的频率作为后续计算的初始频率；

S4，利用正余弦函数逼近法重构单频干扰信号，将原始地震信号与重构单频干扰信号差的平方求和定义为目标函数，对目标函数求导得到单频干扰的振幅；

S5，将定义的目标函数转换为原始地震信号与重构单频干扰的归一化互相关函数，采用频率折半逼近算法进行频率扫描，得到多个互相关函数值，将其极大值对应的频率视为单频干扰的最佳频率；

S6，将单频干扰频率及用其计算的振幅代入正余弦函数公式，得到重构的单频干扰信号；

S7，从原始地震数据中减去重构的单频干扰信号，获得无单频干扰的地震数据，完成消除单频干扰处理工作。

2.根据权利要求1所述的一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，地震信号为离散时间序列，根据下面公式计算离散地震信号的n阶导数：

公式(1-1)中X^T表示原始地震信号，(X^T)′表示一阶求导信号，i表示采样点号，N表示采样点总数；

根据上述公式，把一阶求导信号作为原始信号再次求导，可得到二阶求导函数；依此类推，最终得到n阶求导函数(X^T)⁽ⁿ⁾。

3.根据权利要求2所述的一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，根据下述变换公式得到地震数据导数的频谱：

F[(X^T)⁽ⁿ⁾]＝(iω)ⁿF(ω) (1-2)

其中，(X^T)⁽ⁿ⁾为地震信号的n阶导数，F[·]表示傅氏变换符号，ω为数字频率，F(ω)为幅度函数。

4.根据权利要求3所述的一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，对地震数据导数的频谱进行叠加平滑处理，以消除锯齿状振幅起伏现象对单频干扰识别的影响。

5.根据权利要求4所述的一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，根据下述公式进行频谱叠加处理，获得振幅平滑后的频谱：

对所选取的原始地震数据求导后得到采样点总数为K的地震信号，从该信号中分段连续截取N个样点，通过傅氏变换得到的第i段频谱为：

F_i(n)1≤n≤N且N≤K (1-3)

其中，i表示截取的起始样点，其取值范围为1≤i≤K+1-N，F₁(n)表示[1，N]范围内的样点计算得到的频谱，F₅(n)表示[5，4+N]范围内的样点计算得到频谱。

将所有分段频谱按照下述公式进行叠加处理，可获得振幅平滑后的频谱：

6.根据权利要求1所述的一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，根据下面公式定义单频干扰判别准则：

已知采样率为F，样点数为N的一段地震信号的导数进行快速傅氏变换，将获得的N个复数点取模，即得到各复数点对应的振幅，记为：

A(n)1≤n≤N (1-5)

如果振幅A(n)满足下面不等式：

A(n-1)＜A(n)，且A(n)＞A(n+1)1＜n＜N (1-6)

则将其对应的频率定义为凸起频率，并将每个凸起频率对应的振幅标记为P(i)，凸起频率总数记为H；

如果P(i)满足下面不等式：

P(i-1)＜P(i)，P(i)＞P(i+1)1＜i＜H (1-7)

则其对应频率的信号被判定为疑似干扰；

若第j个疑似干扰频率的振幅P(j)使下面不等式成立，则疑似干扰即被判定为单频干扰；

不等式(1-8)表示疑似干扰对应的振幅与周边4个凸起频率对应振幅的均值做比，如果比值大于设定的阈值M，则疑似干扰被确定为单频干扰；式(1-8)即为单频干扰判别准则。

7.根据权利要求1所述的一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，在步骤S6中，利用正余弦函数法重构单频干扰信号时，将原始地震采样序列与重构单频干扰序列之差的平方求和定义为目标函数，对目标函数求导可得到单频干扰的初始频率对应的振幅；正余弦函数重构单频干扰公式、目标函数及振幅计算公式如下：

y_i＝Acos2πfiΔt+Bsin2πfiΔt (1-9)

其中，y_i表示重构的单频干扰，A表示余弦信号的振幅，B表示正弦信号的振幅，f表示单频干扰的频率，i为采样点号，Δt为地震数据采样时间间隔。

公式(1-9)的向量表达式如下：

y＝Bs^T+Aa^T (1-10)

其中，s^T表示单频干扰的正弦序列，a^T表示单频干扰的余弦序列，两者展开式如下：

由原始采样序列估算出单频干扰序列，定义的目标函数如下：

其中，x_i为原始地震信号，Q表示原始采样序列减去重构的单频干扰序列后的能量差值，X^T为原始采样序列；

利用公式(1-12)，对Q关于各变量求导，可获得Q的极小值，将能使Q取极小值的重构序列视为实际的单频干扰序列；分别对Q关于正弦分量振幅B和余弦分量振幅A求导，并使导数为0，可得到一个二元一次方程组：

求解方程组(1-13)，可以获得A和B的值，计算公式如下：

8.根据权利要求7所述的一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，原始地震采样序列与重构单频干扰采样序列的归一化互相关函数式如下：

将公式(1-12)改写成向量形式，可得到原始地震采样序列与重构单频干扰采样序列的归一化互相关函数，其公式为：

其中，X，Y分别为原始地震信号和重构的单频干扰信号，Corr(f)为互相关函数值；

由公式(1-15)可知，不同的扫描频率f可以计算出不同的互相关函数值。

9.根据权利要求8所述的一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，利用公式(1-15)，将识别出的单频干扰频率作为扫描的初始频率，采用频率折半逼近算法进行频率扫描可得到多个互相关函数值，其极大值对应的频率即为单频干扰的最佳频率。

10.根据权利要求9所述的一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，频率折半逼近算法公式如下：

给定单频干扰频率扫描的初始值f₀和初始频率扫描范围Δf，则第一级频率扫描样点公式：

f_k＝f₀+Δfk k＝[-1，0，1] (1-16)

从第一级的三个频率样点中找到能使函数值Corr(f)为最大时的频率，并记为f₁，按照公式(1-17)计算第二级的频率样点f₂，f₃。将f₁，f₂，f₃列为第二级频率样点；

从第二级的三个频率样点中再次寻找Corr(f)为最大时的频率，记为f₁；按照公式(1-18)计算第三级的频率样点f₂，f₃；将f₁，f₂，f₃划归为第三级频率样点；

以此类推，将上一级使Corr(f)最大时的频率记为f₁，本级的频率样点f₂，f₃，按照公式如下：

式(1-19)中n表示扫描等级，将f₁，f₂，f₃划归为第n级频率样点，第n级频率扫描实际计算了2个频率样点，其精度达到(1/2)^(n-1)ΔfHz。当(1/2)^(n-1)Δf达到了预期的频率精度时，频率f₁即为逼近的单频干扰频率。

11.根据权利要求10所述的一种高精度自动识别和压制地震数据单频干扰的方法，其特征在于，将扫描得到的单频干扰的最佳频率f₁及按公式(1-14)重新计算出的振幅A和B代入公式(1-9)，得到重构的单频干扰信号；然后从原始地震数据中减去重构的单频干扰，即可得到无单频干扰的地震数据；

去除单频干扰计算公式如下：

S＝X-Y (1-20)

其中，X为原始地震数据，Y为重构的单频干扰，S为消除单频干扰后的地震数据。