CN102073066A - 一种消除地震数据谐波干扰的方法 - Google Patents

Abstract

本发明涉及油田的勘探中消除地震数据谐波干扰的方法，根据地震数据中原始波形数据，由原始波形震荡周期和振幅谱的最大位置，确定原始数据中谐波干扰的初始频率，计算谐波干扰余弦函数的振幅，根据地震数据向量、谐波干扰余弦函数的振幅、谐波干扰余弦函数向量，目标函数目导数，谐波干扰频率修正量和时延修正量，采用频率修正公式确定谐波干扰频率，计算地震有效信号。本发明克服了频率域压制谐波干扰的缺点，比时间域谐波干扰波压制方法速度快，提高信噪比，既可以消除地震数据中由高压输电线产生的谐波干扰，也可以消除地震数据中由周期性震动产生的谐波干扰。

Description

一种消除地震数据谐波干扰的方法

技术领域

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

背景技术

地震勘探的过程，就是在地面上的一系列点上，利用人工激发地震波，地震波向地下传播，当遇到波阻抗(地震波在地层介质中向地下传播的速度与介质密度的乘积)界面(即上下地层波阻抗不相等面)时，在波阻抗界面上地震波产生反射现象，地震波传播方向发生改变，地震波开始向上传播，在地面上的一系列接收位置上安置着接收器，接收向上传播的地震波数据，完成野外勘探。

在地震数据野外采集过程中，如果地震测线从高压输电线下面或者旁边通过，由于高压输电线的电流会产生很强的电磁场，这个电磁场也会引起地震检波器周期性振荡，在地震数据记录中记下这个周期性振荡，即工业交流电谐波。谐波是地震数据中的干扰，它的存在，污染了地震反射信号，有时甚至完全掩盖了地震反射信号。野外采集时，地震观测系统已经经过认真仔细的设计，不可随意改动。这样在高压输电线通过的地区进行地震勘探时，谐波干扰是不可避免的，并且野外采集过程中是无法克服的，只有在室内地震数据处理过程中，作为干扰加以消除。在野外地震数据采集过程中，如果在地面接收器附近存在高压输电线或周期性震动(如发电机等)，这样在地面接收器接收到的地震数据中就会存在很强的谐波干扰，其频率在整个接收长度上是固定不变的，它与地下地震地质条件无关，与激发的地震信号无关，与地表地震地质条件无关。因此在地震勘探和地震数据处理中，这种波被看作为干扰，必须加以剔除。

在地震记录中存在谐波干扰时，常规的压制方法是在频率域内进行压制。频率域处理虽然简单、方便，但是存在以下问题，在浅层，当有效波与干扰的能量水平非常接近，或者有效波能量比干扰的能量强，则干扰不易识别；如果有效波的能量比干扰的能量弱，此时干扰容易识别。在深层，干扰易识别。同时频率域处理对于干扰仅仅在振幅上进行压制处理，压制量不易掌握，压制不足会在记录上存在残余的谐波干扰，而压制过量会伤害有效信号。频率域压制还往往损害该频率附近有效波频率成分；为了减少对有效信号频率的损害，就要选取很窄的压制频带，这样对应的时间域算子很长，会产生严重的边界效应。同时由于谐波干扰的频率受到周波不稳的影响，往往不是纯粹的50hz，同时还受到计算时窗选取的影响，使得快速傅里叶变换存在一些难以克服的问题。这些问题都使得在频率域内有效地压制谐波干扰难以实现。

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

发明内容

本发明目的在于提供一种计算简单、效果显著的直接在时间域内消除地震数据谐波干扰的方法。

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

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

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

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

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

3)计算谐波干扰余弦函数的振幅A；

步骤3)所述的计算谐波干扰余弦函数的振幅A，其计算公式为：

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

其中

x^T＝(x₁，x₂，Λ，x_N)(2)

c^T＝(cos2πf(1+τ)Δt，cos2πf(2+τ)Δt，Λ，cos2πf(N+τ)Δt)(3)

式中：

A -----表示谐波干扰余弦函数的振幅；

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

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

f -----表示谐波干扰频率，由自适应频率计算算法确定；

τ-----表示谐波干扰时延，由自适应时延计算算法确定；

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

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

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

4)根据地震数据向量x、谐波干扰余弦函数的振幅A、谐波干扰余弦函数向量c，计算目标函数目Q，标函数Q的计算公式为

Q＝(x-Ac)^T(x-Ac)＝x^Tx-2Ax^Tc+A²c^Tc    (4)

5)根据地震数据向量x、谐波干扰余弦函数的振幅A、谐波干扰余弦函数向量c、谐波干扰正弦函数向量s、时延系数矩阵W、地震记录时间采样间隔Δt，计算目标函数目一阶导数

和二阶导数和

其中：Q是目标函数，

和

是目标函数对频率的一阶导数和目标函数对频率的二阶导数；和

是目标函数、目标函数对时延的一阶导数和目标函数对时延的二阶导数；

目标函数Q对频率f的一阶导数

为

${\▿}_{f}Q=4\mathrm{\πA\Δt}{x}^T\mathrm{Ws}-4\mathrm{\π}{A}^2\mathrm{\Δt}{c}^T\mathrm{Ws}---\left(5\right)$

目标函数Q对频率f的二阶导数

为

${\▿}_f^{2}Q=8{\mathrm{\π}}^{2}A{\mathrm{\Δt}}^2{x}^T{W}^T\mathrm{Wc}-8{\mathrm{\π}}^2{A}^2{\mathrm{\Δt}}^2{(c+s)}^T{W}^{T}W(c-s)---\left(6\right)$

目标函数Q对时延τ的一阶导数

为

${\▿}_{\mathrm{\τ}}{Q}_k=4\mathrm{\πAf\Δt}{x}^{T}s-4\mathrm{\π}{A}^2\mathrm{f\Δt}{c}^{T}s---\left(7\right)$

目标函数Q对时延τ的二阶导数

为

${\▿}_{\mathrm{\τ}}^{2}Q=8{\mathrm{\π}}^2{A}^2{f}^2{\mathrm{\Δt}}^2{x}^{T}c-8{\mathrm{\π}}^2{A}^2{f}^2{\mathrm{\Δt}}^2{(c+s)}^T(c-s)---\left(8\right)$

这里

s^T＝(sin 2πf(1+τ)Δt，sin 2πf(2+τ)Δt，Λ，sin2πf(N+τ)Δt)(9)

$W=\left[\begin{array}{cccc}1+\mathrm{\τ}& 0& \mathrm{\Λ}& 0\\ 0& 2+\mathrm{\τ}& \mathrm{\Λ}& 0\\ M& M& \mathrm{\Λ}& M\\ 0& 0& \mathrm{\Λ}& N+\mathrm{\τ}\end{array}\right]---\left(10\right)$

式中：

A  -----表示谐波干扰余弦函数的振幅；

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

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

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

f  -----表示谐波干扰频率，由自适应频率计算算法确定；

τ-----表示谐波干扰时延，由自适应时延计算算法确定；

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

T  -----表示向量转置或者矩阵转置；

W  -----表示系数矩阵，由谐波干扰时延确定；

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

6)根据目标函数一阶导数

和

二阶导数

和计算谐波干扰频率修正量Δf和时延修正量Δτ，其计算公式为：

采用以下公式计算确定谐波干扰频率修正量Δf，

$\mathrm{\Δf}=\frac{1}{{\▿}_f^{2}Q}{\▿}_{f}Q---\left(11\right)$

采用以下公式计算确定谐波干扰时延修正量Δτ，

$\mathrm{\Δ\τ}=\frac{1}{{\▿}_{\mathrm{\τ}}^{2}Q}{\▿}_{\mathrm{\τ}}Q---\left(12\right)$

7)根据谐波干扰初始频率f₀、初始时延τ₀、频率修正量Δf和时延修正量Δτ，计算确定谐波干扰频率f和时延τ，采用以下频率修正公式确定谐波干扰频率f，

f＝f₀-γΔf    (13)

采用以下时延修正公式确定谐波干扰时延τ，

τ＝τ₀-μΔτ(14)

其中f₀是初始频率，γ是迭代频率自适应步长，通过实验确定频率自适应步长；τ₀是初始时延，在计算时初始时延值可以设置为0；μ是迭代时延自适应步长，通过实验确定频率自适应步长；初始频率f₀由步骤2)确定；计算第一步中，f为f₀，τ为τ₀；

8)根据谐波干扰频率f、谐波干扰时延τ、谐波干扰余弦函数的振幅A、地震记录的时间采样间隔Δt；

步骤8)所述的谐波干扰y_i按照以下公式计算；

y_i＝Acos2πf(i+τ)Δt    (15)

式中频率f由自适应频率计算算法确定，时延τ由自适应时延计算算法确定，谐波干扰余弦函数的振幅A由振幅计算公式计算确定；

9)计算地震有效信号；

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

S_i＝x_i-y_i    (16)

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

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

本发明克服了频率域压制谐波干扰的缺点，比时间域谐波干扰波压制方法速度快，能够有效的消除地震记录上的谐波干扰，还保留了时间域单频干扰波压制方法不损害有效波的特点。

本发明提高了该频率分量的信噪比，既可以消除地震数据中由高压输电线产生的谐波干扰，也可以消除地震数据中由周期性震动(如发电机等)产生的谐波干扰。

附图说明

图1合成理论数据对比

  (a)原始数据

  (b)合成谐波干扰

  (c)合成数据

  (d)消除谐波干扰后数据

  (e)本发明检测出谐波干扰

图2合成理论数据频谱对比

  (a)原始数据频谱

  (b)合成谐波干扰频谱

  (c)合成数据频谱

  (d)消除谐波干扰后数据频谱

  (e)本发明检测出谐波干扰频谱

图3实际炮集数据对比

  (a)原始炮集数据

  (b)陷波处理后炮集数据

  (c)时间域单频干扰压制方法消除谐波干扰处理后炮集数据

  (d)本发明谐波干扰消除处理后炮集数据

  (e)本发明检测出的谐波干扰炮集数据

图4实际数据频谱对比

  (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  -----表示地震记录的时间采样序号；

为了有效地消除地震记录上的谐波干扰，设谐波干扰的频率、振幅和时延在整个地震记录道内是稳定不变的，且为常数。其表达式是：

y_i＝AcoS2πf(i+τ)Δt    (2)

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

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

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

其中

x^T＝(x₁，x₂，Λ，x_N)(4)

c^T＝(cos2πf(1+τ)Δt，cos2πf(2+τ)Δt，Λ，cos2πf(N+τ)Δt)(5)

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

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

$\frac{\∂Q}{\∂A}=0---\left(6\right)$

得出谐波干扰余弦函数的振幅A为

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

这样，对于一个给定频率f和时延τ，通过方程(5)可以确定出谐波干扰余弦函数向量c，由地震数据向量x，根据方程(7)可以求解出谐波干扰余弦函数的振幅A。

5)自适应频率和时延计算算法。采用自适应频率计算算法估算谐波干扰的频率。首先计算目标函数Q对频率f的一阶导数

${\▿}_{f}Q=4\mathrm{\πA\Δt}{x}^T\mathrm{Ws}-4\mathrm{\π}{A}^2\mathrm{\Δt}{c}^T\mathrm{Ws}---\left(8\right)$

目标函数Q对频率f的二阶导数

${\▿}_f^{2}Q=8{\mathrm{\π}}^{2}A{\mathrm{\Δt}}^2{x}^T{W}^T\mathrm{Wc}-8{\mathrm{\π}}^2{A}^2{\mathrm{\Δt}}^2{(c+s)}^T{W}^{T}W(c-s)---\left(9\right)$

目标函数Q对时延τ的一阶导数为

${\▿}_{\mathrm{\τ}}{Q}_k=4\mathrm{\πAf\Δt}{x}^{T}s-4\mathrm{\π}{A}^2\mathrm{f\Δt}{c}^{T}s---\left(10\right)$

目标函数Q对时延τ的二阶导数

为

${\▿}_{\mathrm{\τ}}^{2}Q=8{\mathrm{\π}}^2{A}^2{f}^2{\mathrm{\Δt}}^2{x}^{T}c-8{\mathrm{\π}}^2{A}^2{f}^2{\mathrm{\Δt}}^2{(c+s)}^T(c-s)---\left(11\right)$

这里

s^T＝(sin 2πf(1+τ)Δt，sin2πf(2+τ)Δt，Λ，sin 2πf(N+τ)Δt)(12)

$W=\left[\begin{array}{cccc}1+\mathrm{\τ}& 0& \mathrm{\Λ}& 0\\ 0& 2+\mathrm{\τ}& \mathrm{\Λ}& 0\\ M& M& \mathrm{\Λ}& M\\ 0& 0& \mathrm{\Λ}& N+\mathrm{\τ}\end{array}\right]---\left(13\right)$

式中，A表示谐波干扰余弦函数的振幅；x表示地震记录的时间序列向量，野外地震数据采集得到；c表示谐波干扰余弦函数向量，根据频率、时间采样率计算得到；s表示谐波干扰正弦函数向量，根据频率、时间采样率计算得到；f表示谐波干扰频率，由自适应频率计算算法确定；τ表示谐波干扰时延，由自适应时延计算算法确定；Δt表示地震记录时间采样间隔，野外地震数据采集得到T表示向量转置或者矩阵转置；W表示系数矩阵，由谐波干扰时延确定；N表示地震记录时间采样长度，野外地震数据采集得到；初始频率f₀由步骤2)确定；计算第一步中，f为f₀，τ为τ₀。则频率和时延的增量估算公式分别为

$\mathrm{\Δf}=-\frac{1}{{\▿}_f^{2}Q}{\▿}_{f}Q---\left(14\right)$

$\mathrm{\Δ\τ}=-\frac{1}{{\▿}_{\mathrm{\τ}}^{2}Q}{\▿}_{\mathrm{\τ}}Q---\left(15\right)$

谐波干扰频率和时延的修正公式分别为

f＝f₀+γΔf    (16)

τ＝τ₀+μΔτ(17)

这里f₀是初始频率，Δf是频率修正值，γ是迭代频率自适应步长，其取值范围为0＜γ≤1，通过实验确定频率自适应步长；τ₀是初始时延，Δτ是时延修正值，μ是迭代时延自适应步长，其取值范围为0＜μ≤1，通过实验确定时延自适应步长。

自适应频率和时延计算为：

(1)给定地震记录的时间序列向量x、地震记录的时间采样间隔Δt，给定最大迭代次数或者最小误差；

(2)分析确定谐波干扰初始频率f₀，设置谐波干扰初始时延τ₀＝0；

(3)设置k＝0；

(4)设置f＝f₀和τ＝τ₀。通过方程(5)和(12)可以确定出谐波干扰余弦函数向量c和谐波干扰正弦函数向量s，由方程(7)计算出谐波干扰余弦函数的振幅A。

(5)由方程(8)、(9)和(10)、(11)式分别计算目标函数对频率一阶导数

二阶导数

和目标函数对时延一阶导数二阶导数

(6)由方程(14)和(15)式分别计算频率修正值Δf和时延修正值Δτ；

(7)由方程(16)和(17)式分别计算频率修正量f和时延修正量τ；

(8)由方程(3)式计算误差Q；

(9)设置f₀＝f，τ₀＝τ；

(10)判断迭代次数或者误差是否满足要求，如果不满足要求转向步骤(4)执行，满足要求则执行步骤(11)；

(11)频率f和时延τ就是所要求的谐波干扰频率f和时延τ。

6)计算谐波干扰。根据谐波干扰频率f、时延τ、谐波干扰余弦函数的振幅A、地震记录的时间采样间隔Δt，按照公式(2)计算谐波干扰y_i。

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

S_i＝x_i-y_i    (18)

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

8)根据得出消除谐波干扰的地震数据绘制消除谐波干扰后的地震数据剖面和存储消除谐波干扰后的地震数据。

本发明的一种消除地震数据谐波干扰方法，是一个完整的流程，能够正确有效地识别和消除地震数据中的谐波干扰，是识别和消除地震数据中谐波干扰方法中最有效的方法。

本发明所述的一种消除地震数据谐波干扰方法，其特征在于为识别和消除地震数据中谐波干扰，找到了一条有效的途经，并应用于实际地震数据处理之中。

本发明所述的一种消除地震数据谐波干扰方法，其特征在于把谐波干扰表示为同频率不同振幅的余弦函数和正弦函数之和，并从记录中减去的方法，来消除谐波干扰。

本发明所述的一种消除地震数据谐波干扰方法，其特征在于去除谐波干扰频率分量时，对于其它频率成分没有任何损害，因此是消除谐波干扰最有效的方法。

本发明所述的一种消除地震数据谐波干扰方法，其特征在于使用深层地震数据来估计谐波干扰，最有效地估算谐波干扰的能量，这样可以达到最大限度地消除谐波干扰频率成分，而使该频率分量上的有效波受到的损害最小，提高了该频率分量附近的信噪比。

本发明消除野外采集的震源信号的地震数据谐波干扰方法，其特征在于仅仅对谐波干扰进行处理，而不伤害数据中的有效信号。从而为地震数据的后续处理提供了必要的输入地震数据。

本发明既可以消除地震数据中由高压输电线产生的谐波干扰，其特征在于也可以消除地震数据中由周期性震动(如发电机等)产生的谐波干扰。

本发明实例如下：

首先使用理论合成数据进行试算。理论合成数据采用一个实际地震道与一个正余弦函数之和。正余弦函数使用公式：

y_i＝Ccos2πfiΔt+D sin2πfiΔt    (19)

这里，C和D分别为余弦函数和正弦函数的振幅，f为频率，τ为延迟，Δt为时间采样间隔。正余弦函数振幅C和D与对应余弦函数的振幅A和延迟τ的关系为

$A=\sqrt{C^2+D^2},$ $\mathrm{\τ}=\frac{1}{2\mathrm{\πf\Δt}}\mathrm{arctg}(-\frac{D}{C})---\left(20\right)$

C＝Acos2πfτΔt，D＝-Asin2πfτΔt(21)

合成余弦函数和正弦函数的振幅C和D分别为7.612063E+11和7.847490E+11，频率f为49.970Hz，时间采样间隔Δt为1ms。由方程式(20)可以计算出对应余弦函数A为1.093282E+12，延迟τ为-2.550。

图1-2是本方法合成理论数据试算对比图。图1理论数据试算对比图，每种数据显示五道，(a)是一个实际地震道，其数据的最大值是1.093282E+13，数据道长是5000ms，时间采样间隔是1ms，这里显示了前4000ms。(b)是理论谐波干扰记录，由方程(19)制作的合成理论谐波干扰记录，(c)是合成理论记录，是实际地震道(a)与理论谐波干扰记录(b)之和，(d)是使用本发明去除谐波干扰后的地震道，(e)是使用本发明检测出的谐波干扰，其中频率是49.970Hz，振幅是1.093277E+12，时延是-2.550，由方程式(21)可以计算出对应方程(19)余弦函数和正弦函数的振幅C和D分别为7.612026E+11和7.847452E+11。图2理论数据频谱对比图，(a)是一个实际地震道频谱，(b)是理论谐波干扰记录频谱，(c)是合成理论记录频谱，(d)是使用本发明去除谐波干扰后的地震道频谱，(e)是使用本发明检测出的谐波干扰频谱。从数据显示和频谱图中可以看出，本发明有效地消除了谐波干扰。

图3-4是实际地震数据处理效果对比和频谱对比图。图3是某地区实际野外炮集数据处理效果对比，这里显示了一个炮集。(a)是原始炮集数据，炮集数据每炮321道，中间放炮两边接收，记录长度5000ms，采样间隔1ms。显然地震数据中存在着很强的谐波干扰，必须进行谐波干扰压制处理。(b)是陷波处理后炮集数据，陷波器算子长度为500ms，陷波器频带宽度为8Hz，即NH(46，50，54)，NH(146，150，154)和NH(246，250，254)，陷波处理后数据浅层去除了谐波干扰，而深层还存在很强的残留谐波干扰。而且在地震数据的起始和结尾处存在有很强的边界效应，这种边界效应是不可避免的，任何频率域处理方法都会存在。(c)是时间域单频干扰压制方法消除谐波干扰处理后炮集数据，显然有效地消除了地震数据上谐波干扰。(d)是本发明谐波干扰消除处理后炮集数据，显然有效地消除了地震数据上谐波干扰。(e)本发明检测出的谐波干扰炮集数据。图4是实际数据频谱对比，作频谱分析数据的时间范围是(1000-4500ms)，显示了三道地震数据的频谱，左边是第131道数据的频谱，其谐波干扰能量与有效波相比要弱一点，中间是第31道数据的频谱，其谐波干扰能量与有效波相比差不多，右边是第231道数据的频谱，其谐波干扰能量明显要强于有效波。(a)是原始数据振幅谱，(b)是陷波处理后数据振幅谱，(c)是时间域单频干扰压制方法消除谐波干扰处理后振幅谱，(d)是本发明处理后数据振幅谱。对于浅层地震数据，原始数据浅层谐波干扰和有效信号能量相当，在陷波处理后数据频谱上，谐波干扰得到了消除，但是谐波干扰附近有效信号也被消除，因此伤害了有效信号，而本方法一点也没有伤害有效信号。对于深层地震数据，原始数据谐波干扰能量比有效信号能量要强得多，在陷波处理后数据频谱上，谐波干扰得到了一定的压制，但是还存在很强的残余谐波干扰，而本方法有效地消除了谐波干扰，一点没有残余谐波干扰。总之，陷波法处理后浅层谐波有一定的压制，但伤害了谐波干扰附近有效信号，因此降低了谐波干扰附近有效信号的信噪比；而本方法有效地消除了谐波干扰，而并不伤害干扰频率附近有效信号，因此提高了谐波干扰频率附近有效信号的信噪比。同时很明显，本方法与时间域单频干扰压制方法一样，有效地消除了谐波干扰，且两种算法看不出差异。在DELL工作站上，该数据一炮321道，本发明的运算时间是45s，而同样的机器同样的参数，时间域单频干扰压制方法需要178s，本发明明显优于时间域单频干扰压制方法消除谐波干扰。因此本发明是消除地震数据谐波干扰的一种十分有效方法。

Claims (6)

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

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

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

3)计算谐波干扰余弦函数的振幅A；

4)根据地震数据向量x、谐波干扰余弦函数的振幅A、谐波干扰余弦函数向量c，计算目标函数目Q，标函数Q的计算公式为：

Q＝(x-Ac)^T(x-Ac)＝x^Tx-2Ax^Tc+A²c^Tc    (4)

5)根据地震数据向量x、谐波干扰余弦函数的振幅A、谐波干扰余弦函数向量c、谐波干扰正弦函数向量s、时延系数矩阵W、地震记录时间采样间隔Δt，计算目标函数目一阶导数和

二阶导数

和

 其中：Q是目标函数，

和是目标函数对频率的一阶导数和目标函数对频率的二阶导数；

和

是目标函数、目标函数对时延的一阶导数和目标函数对时延的二阶导数；

目标函数Q对频率f的一阶导数

为

${\▿}_{f}Q=4\mathrm{\πA\Δt}{x}^T\mathrm{Ws}-4\mathrm{\π}{A}^2\mathrm{\Δt}{c}^T\mathrm{Ws}---\left(5\right)$

目标函数Q对频率f的二阶导数为

${\▿}_f^{2}Q=8{\mathrm{\π}}^{2}A{\mathrm{\Δt}}^2{x}^T{W}^T\mathrm{Wc}-8{\mathrm{\π}}^2{A}^2{\mathrm{\Δt}}^2{(c+s)}^T{W}^{T}W(c-s)---\left(6\right)$

目标函数Q对时延τ的一阶导数

为

${\▿}_{\mathrm{\τ}}{Q}_k=4\mathrm{\πAf\Δt}{x}^{T}s-4\mathrm{\π}{A}^2\mathrm{f\Δt}{c}^{T}s---\left(7\right)$

目标函数Q对时延τ的二阶导数

为

${\▿}_{\mathrm{\τ}}^{2}Q=8{\mathrm{\π}}^2{A}^2{f}^2{\mathrm{\Δt}}^2{x}^{T}c-8{\mathrm{\π}}^2{A}^2{f}^2{\mathrm{\Δt}}^2{(c+s)}^T(c-s)---\left(8\right)$

这里

s^T＝(sin 2πf(1+τ)Δt，sin 2πf(2+τ)Δt，Λ，sin 2πf(N+τ)Δt)(9)

$W=\left[\begin{array}{cccc}1+\mathrm{\τ}& 0& \mathrm{\Λ}& 0\\ 0& 2+\mathrm{\τ}& \mathrm{\Λ}& 0\\ M& M& \mathrm{\Λ}& M\\ 0& 0& \mathrm{\Λ}& N+\mathrm{\τ}\end{array}\right]---\left(10\right)$

式中：

A  ------表示谐波干扰余弦函数的振幅；

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

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

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

f   -----表示谐波干扰频率，由自适应频率计算算法确定；

τ  -----表示谐波干扰时延，由自适应时延计算算法确定；

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

T  -----表示向量转置或者矩阵转置；

W  -----表示系数矩阵，由谐波干扰时延确定；

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

6)根据目标函数一阶导数

和

二阶导数

和

计算谐波干扰频率修正量Δf和时延修正量Δτ，其计算公式为：

采用以下公式计算确定谐波干扰频率修正量Δf，

$\mathrm{\Δf}=\frac{1}{{\▿}_f^{2}Q}{\▿}_{f}Q---\left(11\right)$

采用以下公式计算确定谐波干扰时延修正量Δτ，

$\mathrm{\Δ\τ}=\frac{1}{{\▿}_{\mathrm{\τ}}^{2}Q}{\▿}_{\mathrm{\τ}}Q---\left(12\right)$

7)根据谐波干扰初始频率f₀、初始时延τ₀、频率修正量Δf和时延修正量Δτ，计算确定谐波干扰频率f和时延τ，采用以下频率修正公式确定谐波干扰频率f，

f＝f₀-γΔf    (13)

采用以下时延修正公式确定谐波干扰时延τ，

τ＝τ₀-μΔτ(14)

其中f₀是初始频率，γ是迭代频率自适应步长，通过实验确定频率自适应步长；τ₀是初始时延，在计算时初始时延值可以设置为0；μ是迭代时延自适应步长，通过实验确定频率自适应步长；初始频率f₀由步骤2)确定；计算第一步中，f为f₀，τ为τ₀；

8)根据谐波干扰频率f、谐波干扰时延τ、谐波干扰余弦函数的振幅A、地震记录的时间采样间隔Δt；

9)计算地震有效信号；

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

2.根据权利要求1所述的方法，其特征是步骤1)所述的预处理是指对地震数据，进行置标签、定义观测系统的处理。

3.根据权利要求1所述的方法，其特征是步骤2)所述的确定谐波干扰的初始频率是指根据地震数据中原始波形数据x_i和它的振幅谱，由原始波形数据的震荡周期和它的振幅谱的最大位置所对应的频率，确定原始数据中谐波干扰的初始频率f₀。

4.根据权利要求1所述的方法，其特征是步骤3)所述的计算谐波干扰余弦函数的振幅A，其计算公式为：

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

其中

x^T＝(x₁，x₂，Λ，x_N)(2)

c^T＝(cos2πf(1+τ)Δt，cos2πf(2+τ)Δt，Λ，cos2πf(N+τ)Δt)(3)

式中：

A  -----表示谐波干扰余弦函数的振幅；

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

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

f  ------表示谐波干扰频率，由自适应频率计算算法确定；

τ-----表示谐波干扰时延，由自适应时延计算算法确定；

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

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

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

5.根据权利要求1所述的方法，其特征是步骤8)所述的谐波干扰yi按照以下公式计算；

y_i＝A cos2πf(i+τ)Δt(15)

式中频率f由自适应频率计算算法确定，时延τ由自适应时延计算算法确定，谐波干扰余弦函数的振幅A由振幅计算公式计算确定。

6.根据权利要求1所述的方法，其特征是步骤9)所述的地震有效信号按照以下公式计算：

S_i＝x_i-y_i    (16)

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

Images