CN102073065B - 一种消除地震数据单频干扰的方法 - Google Patents

Abstract

本发明涉及油田的勘探、开发、开采技术，是消除地震数据单频干扰的方法。由快速频率扫描算法确定不等间隔频率采样频率，计算不等间隔频率采样线性调频谱频率样本、线性调频谱样本，将线性调频谱频率样本中最大值所对应的频率确定为单频干扰的频率，根据单频干扰余弦函数和正弦函数的振幅、地震记录的时间采样间隔，计算消除单频干扰后的地震有效信号，本发明速度快，不损害有效波，提高了频率分量的信噪比，本发明也可以消除地震数据中由周期性震动(如发电机等)产生的单频干扰。

Description

一种消除地震数据单频干扰的方法

技术领域

本发明涉及油田的勘探、开发、开采技术，具体是为反映地下地层层位、油藏描述提供高分辨率的地震图形和数据的一种消除地震数据单频干扰的方法。

背景技术

地震勘探的过程，就是在地面上的一系列点上，利用人工激发地震波，地震波向地下传播，当遇到波阻抗(地震波在地层介质中向地下传播的速度与介质密度的乘积)界面(即上下地层波阻抗不相等面)时，在波阻抗界面上地震波产生反射现象，地震波传播方向发生改变，地震波开始向上传播，在地面上的一系列接收位置上安置着接收器，接收向上传播的地震波数据，完成野外勘探。在野外地震数据采集过程中，如果在地面接收器附近存在高压输电线或周期性震动(如发电机等)，这样在地面接收器接收到的地震数据中就会存在很强的单频干扰，其频率在整个接收长度上是固定不变的，它与地下地震地质条件无关，与激发的地震信号无关，与地表地震地质条件无关。因此在地震勘探和地震数据处理中，这种波被看作为干扰，必须加以剔除。

在地震数据野外采集过程中，如果地震测线从高压输电线下面或者旁边通过，由于高压输电线的电流会产生很强的电磁场，这个电磁场也会引起地震检波器周期性振荡，在地震数据记录中记下这个周期性振荡，即单频干扰。单频干扰是地震数据中的干扰，它的存在，污染了地震反射信号，有时甚至完全掩盖了地震反射信号。野外采集时，地震观测系统已经经过认真仔细的设计，不可随意改动。这样在高压输电线通过的地区进行地震勘探时，单频干扰是不可避免的，并且野外采集过程中是无法克服的，只有在室内地震数据处理过程中，作为干扰加以消除。因此在地震数据处理中作为一种干扰，单频干扰必须加以消除。在地震记录中存在单频干扰时，常规的压制方法是在频率域内进行压制。频率域处理虽然简单、方便，但是存在以下问题，在浅层，当有效波与干扰的能量水平非常接近，或者有效波能量比干扰的能量强，则干扰不易识别；如果有效波的能量比干扰的能量弱，此时干扰容易识别。在深层，干扰易识别。同时频率域处理对于干扰仅仅在振幅上进行压制处理，压制量不易掌握，压制不足会在记录上存在残余的单频干扰，而压制过量会伤害有效信号。频率域压制还往往损害该频率附近有效波频率成分；为了减少对有效信号频率的损害，就要选取很窄的压制频带，这样对应的时间域算子很长，会产生严重的边界效应。同时由于单频干扰的频率受到周波不稳的影响，往往不是纯粹的50hz，同时还受到计算时窗选取的影响，使得快速傅里叶变换存在一些难以克服的问题。这些问题都使得在频率域内有效地压制单频干扰难以实现。

采用时间域单频干扰波压制方法(《石油地球物理勘探》，高少武，2001)中公开了在时间域消除地震数据中的单频干扰方法。该方法将单频干扰表示为振幅、频率和时延的余弦函数。可以有效地消除地震记录上的单频干扰。但由于单频干扰的频率和时延两个参数都采用扫描的方法求取，因此运算的速度比较缓慢，特别当数据量很大时非常费时。

快速傅立叶(FFT)的谱分辨率是1/(NΔ)，其中Δ是序列的时间采样间隔，N是序列长度，且N＝2**k，k称为变换的阶数。对于2ms采样的典型记录，0.001Hz的精度要求，意味着N必须不小于524288(即k不小于19)！！！而实际上所需要的只是以被估计频率为中心的一个有限的窄带频谱，在常规FFT方法中大量的计算和存储空间是无效的。因此，一个适应于窄带频谱的高分辨谱分析方法是必须的。

发明内容

本发明目的在于提供一种计算简单、效果明显的直接在时间域内识别和消除地震勘探单频干扰的方法。

本发明采用如下技术方案，包括以下步骤：

1)激发地震震源，采集地震数据并做预处理；

步骤1)所述的预处理是对地震数据进行置标签、定义观测系统的处理。

2)确定单频干扰的初始频率；

步骤2)所述的确定单频干扰的初始频率是指根据地震数据中原始波形数据x_i和它的振幅谱，由原始波形数据的震荡周期和它的振幅谱的最大位置所对应的频率，确定原始数据中单频干扰的初始频率f₀。

3)由快速频率扫描算法确定不等间隔频率采样频率f_k；

4)计算不等间隔频率采样线性调频谱频率样本X_k；

步骤4)所述的计算频率样本X_k采用如下计算方法：

X_k＝(g^Tc(k)-h^Ts(k))-j(h^Tc(k)+g^Ts(k))      (1)

其中

g^T＝(g₀，g₁，...，g_N-1)                   (2)

h^T＝(h₀，h₁，...，h_N-1)                   (3)

c^T(k)＝(c₀(k)，c₁(k)，...，c_N-1(k))       (4)

s^T(k)＝(s₀(k)，s₁(k)，...，s_N-1(k))       (5)

和

g_n＝x_n cos ω₀ n                          (6)

h_n＝x_n sin ω₀ n       (7)

c_n(k)＝cos ω^k n       (8)

s_n(k)＝sin ω^k n       (9)

其中

x_n-----表示地震数据样本，由输入地震数据确定；

ω₀-----表示单频干扰初始圆频率，由步骤2)确定；

ω^k-----表示第k个频率样本的圆频率，由频率算法确定；

s_n(k)-----表示圆频率为ω^k的正弦函数第n个样本，由公式(9)计算确定；

c_n(k)-----表示圆频率为ω^k的余弦函数第n个样本，由公式(8)计算确定；

h_n-----表示地震数据与初始圆频率正弦函数乘积的第n个样本，由公式(7)计算确定；

g_n-----表示地震数据与初始圆频率余弦函数乘积的第n个样本，由公式(6)计算确定；

s(k)-----表示圆频率为ω^k的正弦函数向量，由公式(5)构成；

c(k)-----表示圆频率为ω^k的余弦函数向量，由公式(4)构成；

h-----表示地震数据与初始圆频率正弦函数乘积向量，由公式(3)构成；

g-----表示地震数据与初始圆频率余弦函数乘积向量，由公式(2)构成；

N-----表示输入地震数据样本个数，由输入地震数据确定；

M-----表示线性调频谱频率样本个数，由频率算法确定；

5)计算线性调频谱样本Y_k；

步骤5)所述的计算线性调频谱样本Y_k公式为：

Y_k＝|X_k|，k＝1，2，3，...，M      (10)

6)将线性调频谱频率样本中最大值所对应的频率确定为单频干扰的频率f；

步骤6)所述的单频干扰的频率为：

$Y_{\max }=\underset{m\∈\[1,M\]}{\max }\left\{Y_k\right\},$ k＝1，2，3，...，M        (11)

f＝f_m                 (12)

其中

f-----表示单频干扰频率，由公式(12)确定；

f_m-----表示最大频率谱样本对应的频率，也是第m个频率样本的频率，由公式(11)确定；

Y_k-----表示线性调频谱第k个频率谱样本，由公式(10)计算确定；

X_k-----表示线性调频谱第k个频率样本，由公式(1)计算确定；

Y_max-----表示最大频率谱样本，也是第m个频率谱样本，由公式(11)确定；f就是所要确定的单频干扰的频率。

7)计算单频干扰余弦函数和正弦函数的振幅A和B；

步骤7)所述的计算单频干扰余弦函数和正弦函数的振幅A和B计算公式为：

$A=\frac{\left(x^{T}c\right)\left(s^{T}s\right)-\left(x^{T}s\right)\left(s^{T}c\right)}{\left(s^{T}s\right)\left(c^{T}c\right)-{\left(s^{T}c\right)}^2}---\left(13\right)$

$B=\frac{\left(x^{T}s\right)\left(c^{T}c\right)-\left(x^{T}c\right)\left(s^{T}c\right)}{\left(s^{T}s\right)\left(c^{T}c\right)-{\left(s^{T}c\right)}^2}---\left(14\right)$

其中

x^T＝(x₁，x₂，…，x_N)                          (15)

c^T＝(cos2πfΔt，cos4πfΔt，…，cos2NπfΔt) (16)

s^T＝(sin2πfΔt，sin4πfΔt，…，sin2NπfΔt) (17)

式中：

A-----表示单频干扰余弦函数的振幅；

B-----表示单频干扰正弦函数的振幅；

x-----表示地震记录的时间序列向量，野外地震数据采集得到；

c-----表示单频干扰余弦函数向量，根据频率、时间采样率计算得到；

s-----表示单频干扰正弦函数向量，根据频率、时间采样率计算得到；

f-----表示单频干扰频率，线性调频谱分析法确定；

Δt-----表示地震记录时间采样间隔，野外地震数据采集得到；

T-----表示向量转置；

N-----表示地震记录时间采样长度，野外地震数据采集得到；

8)根据单频干扰频率f、单频干扰余弦函数和正弦函数的振幅A和B、地震记录的时间采样间隔Δt，按照以下公式计算单频干扰y_i；

步骤8)所述的单频干扰y_i按照以下公式计算：

y_i＝Acos2πfiΔt+Bsin2πfiΔt        (18)

式中频率f由自适应频率计算算法确定，单频干扰余弦函数和正弦函数的振幅A和B由振幅计算公式计算确定；

9)计算消除单频干扰后的地震有效信号：

步骤9)所述的地震有效信号按照以下公式计算：

S_i＝x_i-y_i                            (19)

式中：原始地震数据x_i，由野外数据采集得到；估算的单频干扰y_i，由单频干扰计算公式计算得到；S_i是消除单频干扰后的地震有效信号。

10)根据得出消除单频干扰的地震数据绘制消除单频干扰后的地震数据剖面和存储消除单频干扰后的地震数据。

本发明在去除单频干扰频率分量时，仅仅对单频干扰进行处理，而不伤害数据中的有效信号，因此是消除单频干扰有效的方法。

本发明使用深层地震数据或者初至到达时间之前的地震数据来估计单频干扰，最有效地估算单频干扰的能量，可以达到最大限度地压制单频干扰频率成分，而使该频率分量上的有效波受到的损害最小，提高了该频率分量的信噪比。本发明既可以消除地震数据中由高压输电线产生的单频干扰，也可以消除地震数据中由周期性震动(如发电机等)产生的单频干扰。

附图说明

图1为单位圆上频率谱样本对比

a傅立叶变换

b等间隔线性调频变换

c不等间隔线性调频变换。

图2为线性调频谱频率样本分析计算算法框图

图3为理论数据的谱样本对比

a理论单频样本

b DFT谱样本

c等间隔线性调频谱样本

d不等间隔线性调频谱样本。

图4为理论数据单频干扰压制效果对比

a合成单频干扰；

b实际信号；

c理论数据；

d陷频法；

e本方法；

f本方法检测的单频干扰。

图5为理论数据DFT谱分析对比

a合成单频干扰；

b实际信号；

c理论数据；

d陷频法；

e本方法；

f本方法检测的单频干扰。

图6为理论数据等间隔线性调频谱对比

a合成单频干扰；

b实际信号；

c理论数据；

d陷频法；

e本方法；

f本方法检测的单频干扰。

图7为理论数据不等间隔线性调频谱对比

a合成单频干扰；

b实际信号；

c理论数据；

d陷频法；

e本方法；

f本方法检测的单频干扰。

图8为实际数据单频干扰波压制效果对比

a原始数据；

b时间域单频干扰压制法；

c线性调频谱分析法

d本方法检测的单频干扰。

图9实际数据单频干扰波压制效果DFT谱对比

a原始数据；

b时间域单频干扰压制法；

c线性调频谱分析法；

d本方法检测的单频干扰。

图10为实际数据单频干扰波压制效果等间隔线性调频谱对比

a原始数据；

b时间域单频干扰压制法；

c线性调频谱分析法；

d本方法检测的单频干扰。

图11为实际数据单频干扰波压制效果不等间隔线性调频谱对比

a原始数据；

b时间域单频干扰压制法；

c线性调频谱分析法；

d本发明检测的单频干扰。

具体实施方式

本发明把单频干扰表示为同频率不同振幅的余弦函数和正弦函数之和，并从记录中减去的方法来消除单频干扰。其余弦函数和正弦函数的频率采用自适应方法进行估算，余弦函数和正弦函数的振幅采用直接计算方法。

本发明包括以下步骤：

1)用地震震源激发和采集地震数据并做预处理，所述的预处理是指对地震数据置标签、定义观测系统。

2)确定单频干扰的初始频率。采用频谱分析方法，分析地震数据中原始波形数据x_i的频谱，根据地震数据中原始波形数据x_i和它的振幅谱，由原始波形数据的震荡周期和它的振幅谱的最大位置所对应的频率，确定原始数据中单频干扰的初始频率f₀。它并不是原始数据中单频干扰的实际频率。实际频率并不知道，是所要求取的。

3)地震记录为地震有效信号和单频干扰的和，即

x_i＝S_i+y_i                        (1)

式中，x_i表示地震记录；S_i表示地震有效信号；y_i表示单频干扰；i表示地震记录的时间采样序号；

4)单频干扰为同频率不同振幅余弦函数和正弦函数之和。

单频干扰的频率、振幅和时延在整个地震记录道内是稳定不变的，且为常数。本发明使用余弦函数和同频率不同振幅的正弦函数之和来表示单频干扰。其表达式是：

y_i＝Acos2πfiΔt+Bsin2πfiΔt    (2)

式中，y_i表示单频干扰；A表示单频干扰余弦函数的振幅；B表示单频干扰正弦函数的振幅；f表示单频干扰频率；Δt表示地震记录的时间采样间隔；i表示地震记录的时间采样序号；显然参数A和B与y_i是线性关系，而f与y_i是非线性关系，因此A和B的计算要比f容易和快速得多。

5)不等间隔频率采样频率f_k的确定。不等间隔频率采样线性调频谱分析方法的频率由快速频率扫描算法确定。

采用快速频率扫描所使用的频率，即对于初始频率f₀，第一组频率取值为

f_k＝f₀+0.1k，ω^k＝(2π)0.1k，-10≤k≤10   (3)

通过线性调频谱样本计算，可以确定第一组最佳频率f₁。将f₁作为第二组频率取值的初始频率，第二组频率取值为

f_k＝f₁+0.01k，ω^k＝(2π)0.01k，-5≤k≤5   (4)

通过线性调频谱样本计算，可以确定第二组最佳频率f₂。将f₂作为第二组频率取值的初始频率，第三组频率取值为

f_k＝f₂+0.001k，ω^k＝(2π)0.001k，-5≤k≤5 (5)

通过线性调频谱样本计算，可以确定第三组最佳频率f₃。将f₃作为第二组频率取值的初始频率，第四组频率取值为

f_k＝f₃+0.0001k，ω^k＝(2π)0.0001k，-5≤k≤5 (6)

通过线性调频谱样本计算，可以确定第四组最佳频率f₄。依次类推，可以确定满足频率精度要求的单频干扰频率f。这样总共频率扫描次数为54次。这样可以在频率范围[f₀-1，f₀+1]之内通过54次频率扫描可以确定单频干扰的最佳频率f，且其扫描精度为0.0001。这样为了达到频率计算精度Δf＝0.001，则需要计算M＝43个频率样本。为了达到频率计算精度Δf＝0.0001，则需要计算M＝54个频率样本。为了计算确定一个单频干扰频率，需要计算M个频率样本，共需要计算NM+1次正弦函数计算，NM+1次余弦函数计算，8NM次实数乘法，2NM次实数加法。显然运算次数要远远小于等频率采用线性调频谱分析方法频率样本运算次数。对于几十万甚至几千万道的叠前3D地震数据来说，运算效率是相当可观的。

6)线性调频谱分析方法确定单频干扰的频率f。采用线性调频谱分析方法计算确定单频干扰的频率。

线性调频变换分析方法是把离散傅立叶变换(DFT)表示成为褶积为基础的运算方法。这种方法在减少计算复杂性方面比傅立叶变换具有明显的优点。线性调频变换分析方法比快速离散傅立叶变换(DFT)具有更广泛的适应性，它可以用于计算在单位圆上傅立叶变换的任意一组等间隔样本值。令x^T＝(x₁，x₂，…，x_N)表示一个N点序列向量，X(e^jω)表示其傅立叶变换，X^T＝(X₁，X₂，…，X_N)表示其离散傅里叶变换DFT样本向量。则有

$X\left(e^{\mathrm{j\ω}}\right)=\underset{n=1}{\overset{N-1}{\mathrm{\Σ}}}{x}_n{e}^{-\mathrm{j\ωn}}---\left(7\right)$

和

$X_k=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_n{e}^{-j2\mathrm{\πkn}/N}---\left(8\right)$

这里，k＝0，1，…，N-1表示DFT样本序号。研究计算X(e^jω)的M个样本，这些样本在单位圆上以等角度间隔排列，即在如下频率处：

ω_k＝ω₀+kΔω，k＝0，1，2，...，M-1   (9)

这里，ω₀是初始频率，决定频率样本的起始位置，Δω是频率增量，决定频率样本采样的步长，它们可以任意选择。令ω₀＝0，Δω＝2π/N和M＝N，就成为通常的DFT。如图1所示，(a)是DFT变换频率样本，(b)是等间隔线性调频变换频率样本。与这个频率样本集对应的DFT是

$X\left(e^{j{\mathrm{\ω}}_k}\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_n{e}^{-j{\mathrm{\ω}}_{k}n},$ k＝0，1，2，...，M-1    (10)

定义

W＝e^-jΔω          (11)

$A=e^{j{\mathrm{\ω}}_0}---\left(12\right)$

有，

$X\left(e^{j{\mathrm{\ω}}_k}\right)=W^{k^2/2}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{g}_n{W}^{-{(k-n)}^2/2},$ k＝0，1，2，...，M-1    (13)

这里，

$g_n=x_n{A}^{-n}{W}^{n^2/2}---\left(14\right)$

令

$h_n=W^{{-n}^2/2}---\left(15\right)$

$X_n=X\left(e^{j{\mathrm{\ω}}_n}\right),$ n＝0，1，2，...，M-1    (16)

则有

$X_n=W^{n^2/2}{g}_n*h_n---\left(17\right)$

这里，符号“*”表示褶积运算。这样，相当于序列g_n和序列h_n的褶积，然后乘以

序列

是频率线性增加的复指数序列，在信号分析中，这种信号称为线性调频信号，所以这种信号分析方法称为线性调频变换分析方法。线性调频频率样本分析计算过程如图2所示。

等频率采样线性调频变换分析方法具有三个显著优点：(1)可以计算单位圆上任意一段曲线上的等间隔频率样本；而且起始频率可以任意给定；(2)频率样本的点数M可以与输入数据的点数N不相同，即可以计算很少的点；(3)频率采样间隔Δω可以取值很小，以达到细化频率的目的。

对于单频干扰初始频率f₀，在频率范围[f₀-1，f₀+1]内计算线性调频变换样本，取频率样本精度(也是频率采样间隔)Δf＝0.001，则有

ω₀＝2π(f₀-1)，Δω＝2πΔf＝0.002π      (18)

$M=[\frac{2}{\mathrm{\Δf}}+1]=2001---\left(19\right)$

这里，符合“[·]”表示取整。这样为了达到频率计算精度Δf＝0.001，则需要计算2001个频率样本。为了达到频率计算精度Δf＝0.0001，则需要计算20001个频率样本。为了计算确定一个单频干扰频率，需要计算M个频率样本，共需要计算2N次正弦函数计算，2N次余弦函数计算，N+2M次复指数计算，4N+M次复数乘法，一个N点复数序列与一个M点复数序列的褶积计算。

在消除地震数据单频干扰时，希望计算少量DFT样本，最好采用不等间隔频率采样，以达到减少计算工作量的目的。在单位圆上以不等角度间隔排列，即在如下频率处：

ω_k＝ω₀+ω^k，ω^k＝2πf_k，k＝1，2，3，...，M  (20)

由方程(10)，有

$X\left(e^{j{\mathrm{\ω}}_k}\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_n\cos {\mathrm{\ω}}_{k}n-j\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_{n+1}\sin {\mathrm{\ω}}_{k}n---\left(21\right)$

将方程(20)代入方程(21)，令

g_n＝x_n cos ω₀ n，h_n＝x_n sin ω₀ n            (22)

c_n(k)＝cos ω^k n，s_n(k)＝sin ω^k n            (23)

用向量符号表示，有

g^T＝(g₀，g₁，...，g_N-1)，h^T＝(h₀，h₁，...，h_N-1)  (24)

c^T(k)＝(c₀(k)，c₁(k)，...，c_N-1(k))，s^T(k)＝(s₀(k)，s₁(k)，...，s_N-1(k))  (25)

则有

X_k＝(g^Tc(k)-h^Ts(k))-j(h^Tc(k)+g^Ts(k))         (26)

这样，方程(26)就是不等频率采样间隔线性调频变换分析公式。不等频率采样线性调频变换分析方法具有三个显著优点：(1)可以计算单位圆上任意点频率样本；而且起始频率可以任意给定；(2)频率样本的点数M可以与输入数据的点数N不相同，即可以计算很少的点；(3)频率采样间隔Δω可以取值很小，以达到细化频率的目的。这三个优点使得不等频率采样线性调频变换分析方法在消除地震数据单频干扰方面具有明显的优势。

采用线性调频谱分析方法估算单频干扰的频率。在给定初始频率f₀和频率计算精度Δf时，对于地震数据向量x，计算确定调频变换频率样本向量X。计算X的谱为

Y_k＝|X_k|，k＝1，2，3，...，M     (27)

这就是线性调频谱频率样本。取向量Y^T＝(Y₁，Y₂，…，Y_M)中最大值所对应的频率就是单频干扰的频率，即

$Y_{\max }=\underset{k\∈\[0,M-1\]}{\max }\left\{Y_k\right\}---\left(28\right)$

f＝f_k                            (29)

这样就可以确定单频干扰频率f。线性调频谱分析计算单频干扰频率算法如下：

(1)对应地震记录向量x；给定单频干扰初始频率f₀和频率样本精度Δf；

(2)由公式(18)计算初始圆频率ω₀，对第一组频率取值，由公式(3)确定频率f_k和圆频率ω^k，频率样本序号数k＝1，2，...，21；

(3)由公式(22)计算g_n和h_n，由(23)计算c_n(k)和s_n(k)；

(4)由公式(24)和(25)构成向量g、h和c(k)、s(k)；

(5)由公式(26)计算线性调频谱样本X_k；

(6)对第一组21个频率取值，重复(2)----(5)的计算；

(7)由公式(27)计算谱序列Y_k；

(8)由公式(28)和(29)确定f₁；

(9)对第二组、第三组、第四组频率取值，分别由公式(24)、(25)、(26)、确定频率f_k和圆频率ω^k，频率样本数k＝1，2，...，11

(10)重复(2)----(9)的计算，在步骤(8)中，对第二组、第三组、第四组频率取值对应的频率分别是f₂、f₃、f₄；

(11)根据不同精度要求，由公式(28)和(29)可以确定出最佳单频干扰频率f；

7)单频干扰的振幅估算。在时间剖面上，深层时间段高频有效波的能量比浅层高频有效波的能量弱得多。因此可以利用深层时间段来估算单频干扰的振幅和频率，把它们作为整道地震记录上的单频干扰。采用最小二乘法估算单频干扰的振幅。建立目标函数：

$Q=\underset{i-1}{\overset{N}{\mathrm{\Σ}}}{\[x_i-y_i\]}^2={(x-\mathrm{Ac}-\mathrm{Bs})}^T(x-\mathrm{Ac}-\mathrm{Bs})---\left(30\right)$

$=x^{T}x-2B{x}^{T}s-2A{x}^{T}c+B^2{s}^{T}s+2\mathrm{AB}{s}^{T}c+A^2{c}^{T}c$

其中

x^T＝(x₁，x₂，…，x_N)                          (31)

c^T＝(cos2πfΔt，cos4πfΔt，…，cos2NπfΔt) (32)

s^T＝(sin2πfΔt，sin4πfΔt，…，sin2NπfΔt) (33)

式中，T表示向量转置；x表示地震记录的时间序列向量，野外地震数据采集得到；c表示单频干扰余弦函数向量，根据频率、时间采样率计算得到；s表示单频干扰正弦函数向量，根据频率、时间采样率计算得到。

为了得到振幅A和B的确定，目标函数Q对A和B分别求导，且令

$\frac{\∂Q}{\∂A}=0,$ $\frac{\∂Q}{\∂B}=0,---\left(34\right)$

得出单频干扰正余弦函数的振幅A和B分别为

$A=\frac{\left(x^{T}c\right)\left(s^{T}s\right)-\left(x^{T}s\right)\left(s^{T}c\right)}{{\left(s^{T}s\right)\left(c^{T}c\right)-\left(s^{T}c\right)}^2}---\left(35\right)$

和

$B=\frac{\left(x^{T}s\right)\left(c^{T}c\right)-\left(x^{T}c\right)\left(s^{T}c\right)}{{\left(s^{T}s\right)\left(c^{T}c\right)-\left(s^{T}c\right)}^2}---\left(36\right)$

这样，对于一个给定频率f，通过方程(32)和(33)可以确定出单频干扰余弦函数向量c和正弦函数向量s，由地震数据向量x，根据方程(35)和(36)可以求解出单频干扰正余弦函数的振幅A和B。

8)计算单频干扰。根据单频干扰频率f、单频干扰余弦函数和正弦函数的振幅A和B、地震记录的时间采样间隔Δt，按照公式(2)计算单频干扰y_i。

9)计算地震有效信号。消除地震记录上的单频干扰就是在已知地震记录x_i的情况下，通过估算出单频干扰y_i以恢复地震有效信号S_i的处理。即

S_i＝x_i-y_i           (37)

这里S_i是消除单频干扰后的地震记录，即地震有效信号。

10)采用通常的方法根据得出消除单频干扰的地震数据绘制消除单频干扰后的地震数据剖面和存储消除单频干扰后的地震数据。

本发明克服了频率域压制单频干扰的缺点，而且比时间域单频干扰波压制方法运算速度要快得多，不但能够有效的消除地震记录上的单频干扰，而且还保留了时间域单频干扰波压制方法不损害有效波的特点。本发明提高了该频率分量的信噪比，为地震数据的后续处理提供了必要的输入地震数据。本发明既可以消除地震数据中由高压输电线产生的单频干扰，也可以消除地震数据中由周期性震动(如发电机等)产生的单频干扰。

本发明实例如下：

首先计算理论单频数据的DFT样本谱和线性调频谱。

理论单频数据采用方程(2)生成。使用的参数是：单频余弦函数和正弦函数的振幅A和B分别是2.5432和4.8166；单频频率f是50Hz；时间采样间隔Δt是1ms，样本个数N是200，计算的理论数据如图3(a)所示。在DFT中，取初始圆频率ω₀为0Hz，样本个数M与数据样本个数N相同，即M＝N＝200，频率增量Δω＝2π/M＝0.01π，那么频率样本精度(也是DFT频率采用间隔)Δf＝1/MΔt＝5Hz，计算的DFT谱频率样本如图3(b)所示。如果采用DFT实现频率样本精度(也是DFT频率采用间隔)Δf＝0.001Hz，那么DFT样本个数M＝1/ΔfΔt＝1000000，即要计算1000000个DFT样本，是数据样本个数N的5000倍。在等间隔线性调频谱中，如果取初始圆频率ω₀＝49(2π)，终止圆频率ω₀＝51(2π)，频率样本精度(也是线性调频谱频率采用间隔)Δf＝0.001Hz，那么线性调频谱样本个数M＝2001，计算的等间隔线性调频谱谱频率样本如图3(c)所示。而在不等间隔线性调频谱中，取相同的初始圆频率ω₀＝49(2π)，终止圆频率ω₀＝51(2π)，频率样本精度(也是线性调频谱频率采用间隔)Δf＝0.001Hz，那么不等间隔线性调频谱样本个数M＝43，计算的不等间隔线性调频谱谱频率样本如图3(d)所示。

理论单频数据采用方程(2)生成。使用的参数是：单频余弦函数和正弦函数的振幅A和B分别是6.541和4.811；单频频率f是50.254Hz；时间采样间隔Δt是1ms，样本个数N是5000，生成单频干扰数据，一道的数据十道显示如图4(a)所示。信号采用一段实际地震数据，其最大值为单频5倍，如图4(b)所示，生成的理论数据如图4(c)所示。使用陷频法对数据进行了处理，陷频器长度为512ms，频带宽度为8Hz，即NH(46，50，54)，陷波处理后数据浅层去除了单频干扰，而深层还存在很强的残留单频干扰，陷频法处理后的信号如图4(d)所示。通过本方法计算检测出的单频频率f是50.254Hz，与理论数据给定的频率完全相同。计算单频余弦函数和正弦函数的振幅A和B分别是6.542和4.815，与给定数据误差非常小，计算生成的单频干扰数据如图4(f)所示，本方法处理后恢复的信号如图4(e)所示。它们的DFT谱分别如图5(a)、(b)、(c)、(d)、(e)、(f)所示，它们的等间隔线性调频线性调频谱谱分别如图6(a)、(b)、(c)、(d)、(e)、(f)所示，它们的不等间隔线性调频线性调频谱谱分别如图7(a)、(b)、(c)、(d)、(e)、(f)所示。从数据上看，陷频法可以消除单频干扰，但是在920ms处出现一些边界效应，可是从DFT谱和线性调频谱谱上看，陷频法在消除单频干扰的同时，严重的消除了单频干扰以及单频附近的有效信号，因此降低了单频附近信号的信噪比。而本方法有效的识别并消除了单频干扰，且没有伤害单频干扰频率附近信号的频率成分，因此有效的提高了单频干扰附近信号的信噪比。

实际数据是一个野外炮集数据，每炮有180道，数据时间采样间隔是2ms，数据记录长度是6000ms。显示了900-4100ms。显示了其中的一炮。图8是炮集数据对比，图9是炮集数据第81道DFT谱对比，图10是炮集数据第81道等间隔线性调频谱谱对比，图10是炮集数据第81道等间隔线性调频谱谱对比：a是原始数据，b是时间域单频干扰压制法，c是线性调频谱分析法；d是本方法检测的单频干扰。显然原始数据中包含着很强的单频干扰，从DFT谱和线性调频谱谱中也可以清楚地看到单频干扰。从数据以及DFT谱和线性调频谱谱中可以看到时间域单频干扰压制法和线性调频谱分析法都非常有效地消除了地震数据上的单频干扰，而且两种方法效果差异不大。但是对一个1.6G共130559道的地震数据，不等间隔线性调频谱方法数据处理需要112s；等间隔线性调频谱方法数据处理需要557s，是不等间隔线性调频谱方法时间的5倍；而相同数据时间域单频干扰压制法需要208383s，是等间隔线性调频谱方法时间的374倍，是不等间隔线性调频谱方法时间的1860倍。不等间隔线性调频谱方法大大提高了运算效率，节省运算时间，更加适用于海量地震数据处理需要。因此是消除单频干扰的最有效方法。特别适用于野外地震数据采集过程中，地震测线上空或者附近有高压输电线通过时所采集的实际地震数据。

Claims (2)

1.一种消除地震勘探单频干扰的方法，其特征在于采用如下技术方案，包括以下步骤：

1)激发地震震源，采集地震数据并做预处理；

2)确定单频干扰的初始频率；

所述的确定单频干扰的初始频率是指根据地震数据中原始波形数据xi和它的振幅谱，由原始波形数据的震荡周期和它的振幅谱的最大位置所对应的频率，确定原始数据中单频干扰的初始频率f₀；

3)由快速频率扫描算法确定不等间隔频率采样频率f_k；

4)计算不等间隔频率采样线性调频谱频率样本X_k；采用如下计算方法：

X_k＝(g^Tc(k)-h^Ts(k))-j(h^Tc(k)+g^Ts(k))       (1)

其中

g^T＝(g₀，g₁，...，g_N-1)                    (2)

h^T＝(h₀，h₁，...，h_N-1)                    (3)

c^T(k)＝(c₀(k)，c₁(k)，...，c_N-1(k))        (4)

s^T(k)＝(s₀(k)，s₁(k)，...，s_N-1(k))        (5)

和

g_n＝x_n cos ω₀ n                           (6)

h_n＝x_n sin ω₀ n                           (7)

c_n(k)＝cos ω^k n                           (8)

s_n(k)＝sin ω^k n                           (9)

其中

x_n-----表示地震数据样本，由输入地震数据确定；

ω₀-----表示单频干扰初始圆频率，由步骤2)确定；

ω^k-----表示第k个频率样本的圆频率，由频率算法确定；

s_n(k)-----表示圆频率为ω^k的正弦函数第n个样本，由公式(9)计算确定；

c_n(k)-----表示圆频率为ω^k的余弦函数第n个样本，由公式(8)计算确定；

h_n-----表示地震数据与初始圆频率正弦函数乘积的第n个样本，由公式(7)计算确定；

g_n-----表示地震数据与初始圆频率余弦函数乘积的第n个样本，由公式(6)计算确定；

s(k)-----表示圆频率为ω^k的正弦函数向量，由公式(5)构成；

c(k)-----表示圆频率为ω^k的余弦函数向量，由公式(4)构成；

h-----表示地震数据与初始圆频率正弦函数乘积向量，由公式(3)构成；

g-----表示地震数据与初始圆频率余弦函数乘积向量，由公式(2)构成；

N-----表示输入地震数据样本个数，由输入地震数据确定；

M-----表示线性调频谱频率样本个数，由频率算法确定；

5)计算线性调频谱样本Y_k，公式为：

Y_k＝|X_k|，k＝1，2，3，...，M       (10)

6)将线性调频谱频率样本中最大值所对应的频率确定为单频干扰的频率f；

所述的单频干扰的频率为：

$Y_{\max }=\underset{m\∈\[1,M\]}{\max }\left\{Y_k\right\},$ k＝1，2，3，...，M       (11)

f＝f_m                (12)

其中

f-----表示单频干扰频率，由公式(12)确定；

f_m-----表示最大频率谱样本对应的频率，也是第m个频率样本的频率，由公式(11)确定；

Y_k-----表示线性调频谱第k个频率谱样本，由公式(10)计算确定；

X_k-----表示线性调频谱第k个频率样本，由公式(1)计算确定；

Y_max-----表示最大频率谱样本，也是第m个频率谱样本，由公式(11)确定；f就是所要确定的单频干扰的频率；

7)计算单频干扰余弦函数和正弦函数的振幅A和B；

所述的计算单频干扰余弦函数和正弦函数的振幅A和B计算公式为：

$A=\frac{\left(x^{T}c\right)\left(s^{T}s\right)-\left(x^{T}s\right)\left(s^{T}c\right)}{\left(s^{T}s\right)\left(c^{T}c\right)-{\left(s^{T}c\right)}^2}---\left(13\right)$

其中

x^T＝(x₁，x₂，…，x_N)                          (15)

c^T＝(cos2πfΔt，cos4πfΔt，…，cos2NπfΔt) (16)

s^T＝(sin2πfΔt，sin4πfΔt，…，sin2NπfΔt) (17)

式中：

A-----表示单频干扰余弦函数的振幅；

B-----表示单频干扰正弦函数的振幅；

x-----表示地震记录的时间序列向量，野外地震数据采集得到；

c-----表示单频干扰余弦函数向量，根据频率、时间采样率计算得到；

s-----表示单频干扰正弦函数向量，根据频率、时间采样率计算得到；

f-----表示单频干扰频率，线性调频谱分析法确定；

Δt-----表示地震记录时间采样间隔，野外地震数据采集得到；

T-----表示向量转置；

8)根据单频干扰频率f、单频干扰余弦函数和正弦函数的振幅A和B、地震记录的时间采样间隔Δt，按照以下公式计算单频干扰y_i；

y_i＝Acos2πfiΔt+Bsin2πfiΔt      (18)

式中频率f由自适应频率计算算法确定，单频干扰余弦函数和正弦函数的振幅A和B由振幅计算公式计算确定；

9)计算消除单频干扰后的地震有效信号：

10)根据得出消除单频干扰的地震数据绘制消除单频干扰后的地震数据剖面和存储消除单频干扰后的地震数据。

2.根据权利要求1所述的一种消除地震勘探单频干扰的方法，其特征在于步骤1)所述的预处理是对地震数据进行置标签、定义观测系统的处理。

Images