CN101852865B - 一种自适应消除地震勘探工业交流电干扰的方法 - Google Patents

Abstract

本发明是石油勘探地震数据采集的自适应识别和消除地震勘探工业交流电干扰的方法。步骤是确定工业交流电干扰的初始频率，计算工业交流电干扰余弦函数和正弦函数的振幅，根据地震数据向量x和工业交流电干扰初始频率f₀和工业交流电干扰频率f，根据地震记录的时间采样间隔Δt，计算工业交流电干扰y_i和地震有效信号。本发明流电干扰的缺点，提高了频率分量的信噪比，运算速度要快，不损害有效波的特点。

Description

一种自适应消除地震勘探工业交流电干扰的方法

技术领域

本发明涉及油田的勘探、开发、开采技术，具体是为反映地层层位、油藏描述提供高分辨率的地震图形，特别适用于野外地震数据采集测线上空有高压输电线通过时的自适应识别和消除地震勘探工业交流电干扰的方法。

背景技术

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

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

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

发明内容

本发明目的在于提供一种运算速度快、不损害有效波效、提高频率分量的信噪比果的自适应识别和消除地震勘探工业交流电干扰的方法。

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

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

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

2)确定工业交流电干扰的初始频率；

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

3)计算工业交流电干扰余弦函数和正弦函数的振幅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(1\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(2\right)$ 其中

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

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

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

式中：

A  -----表示工业交流电干扰余弦函数的振幅；

B  -----表示工业交流电干扰正弦函数的振幅；

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

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

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

f  -----表示工业交流电干扰频率，自适应频率计算算法确定；

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

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

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

4)根据地震数据向量x和工业交流电干扰初始频率f₀，采用以下频率修正公式确定工业交流电干扰频率f，

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

其中f₀是初始频率；γ是迭代频率自适应步长，通过实验确定频率自适应步长；Q、

和

是目标函数、目标函数对频率的一阶导数和目标函数对频率的二阶导数。目标函数Q的表达式为

Q＝(x-Ac-Bs)^T(x-Ac-Bs)

                        (7)

 ＝x^Tx-2Bx^Ts-2Ax^Tc+B²s^Ts+2ABs^Tc+A²c^Tc

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

为

${\▿}_{f}Q=4\mathrm{\π\Δt}(B^2-A^2)s^T\mathrm{Vc}-4\mathrm{\π\ΔtB}{x}^T\mathrm{Vc}$ (8)

$+4\mathrm{\π\ΔtA}{x}^T\mathrm{Vs}+4\mathrm{\π\ΔtAB}{(c+s)}^{T}V(c-s)$

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

为

${\▿}_f^{2}Q=8{\mathrm{\π}}^2\mathrm{\Δ}{t}^2(B^2-A^2){(c+s)}^T{V}^{T}V(c-s)+8{\mathrm{\π}}^2\mathrm{\Δ}{t}^{2}B{x}^T{V}^T\mathrm{Vs}$ (9)

$+8{\mathrm{\π}}^2\mathrm{\Δ}{t}^{2}A{x}^T{V}^T\mathrm{Vc}-32{\mathrm{\π}}^2\mathrm{\Δ}{t}^2\mathrm{AB}{s}^T{V}^T\mathrm{Vc}$

这里“T”表示向量或者矩阵转置，V是系数矩阵，

$V=\left[\begin{array}{cccc}1& 0& \·\·\·& 0\\ 0& 2& \·\·\·& 0\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ 0& 0& \·\·\·& N\end{array}\right]---\left(10\right)$

初始频率f₀由步骤2)确定；计算第一步中，f为f₀；

5)根据工业交流电干扰频率f、工业交流电干扰余弦函数和正弦函数的振幅A和B、地震记录的时间采样间隔Δt，按照以下公式计算工业交流电干扰y_i＝Acos 2πfiΔt+Bsin 2πfiΔt     (10)

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

6)按照以下公式计算地震有效信号：

S_i＝x_i-y_i                (11)

式中：原始地震数据x_i，由野外数据采集得到；参考的工业交流电干扰y_i，由工业交流电干扰计算公式计算得到；S_i是消除工业交流电干扰后的地震有效信号。

7)采用通常的方法根据得出消除工业交流电干扰的地震数据绘制消除工业交流电干扰后的地震数据剖面和存储消除工业交流电干扰后的地震数据。

本发明克服了频率域压制工业交流电干扰的缺点，提高了频率分量的信噪比，而且比时间域工业交流电干扰波压制方法运算速度要快得多，不但能够有效的消除地震记录上的工业交流电干扰，而且还保留了时间域单频干扰波压制方法不损害有效波的特点。

本发明还可以消除消除地震数据中由周期性震动(如发电机等)产生的工业交流电干扰。

附图说明

图1合成理论数据对比；

(a)原始数据

(b)合成工业交流电干扰

(c)合成数据

(d)消除工业交流电干扰后数据

(e)本发明检测出工业交流电干扰

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

(a)原始数据频谱

(b)合成工业交流电干扰频谱

(c)合成数据频谱

(d)消除工业交流电干扰后数据频谱

(e)本发明检测出工业交流电干扰频谱

图3实际炮集数据对比；

(a)原始炮集数据

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

(c)时间域单频干扰压制方法消除工业交流电干扰处理后炮集数据

(d)本发明工业交流电干扰消除处理后炮集数据

(e)本发明检测出的工业交流电干扰炮集数据

图4实际数据第51炮第81道频谱对比。

(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＝Acos 2πfiΔt+Bsin 2πfiΔt    (2)式中，

y_i -----表示工业交流电干扰；

A  -----表示工业交流电干扰余弦函数的振幅；

B  -----表示工业交流电干扰正弦函数的振幅；

f  -----表示工业交流电干扰频率；

Δt-----表示地震记录的时间采样间隔；

i  -----表示地震记录的时间采样序号；

显然参数A和B与y_i是线性关系，而f与y_i是非线性关系，因此A和B的计算要比f容易和快速得多。

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

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

$=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)                             (4)

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

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

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

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

$\frac{\∂Q}{\∂A}=0,$ $\frac{\∂Q}{\∂B}=0,---\left(7\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(8\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(9\right)$

这样，对于一个给定频率f，通过方程(5)和(6)可以确定出工业交流电干扰余弦函数向量c和正弦函数向量s，由地震数据向量x，根据方程(8)和(9)可以求解出工业交流电干扰正余弦函数的振幅A和B。

5)自适应频率计算算法。采用自适应频率计算算法估算参考工业交流电干扰的频率。首先计算目标函数Q对频率f的一阶导数

${\▿}_{f}Q=4\mathrm{\π\Δt}(B^2-A^2)s^T\mathrm{Vc}-4\mathrm{\π\ΔtB}{x}^T\mathrm{Vc}$ (10)

$+4\mathrm{\π\ΔtA}{x}^T\mathrm{Vs}+4\mathrm{\π\ΔtAB}{(c+s)}^{T}V(c-s)$

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

${\▿}_f^{2}Q=8{\mathrm{\π}}^2\mathrm{\Δ}{t}^2(B^2-A^2){(c+s)}^T{V}^{T}V(c-s)+8{\mathrm{\π}}^2\mathrm{\Δ}{t}^{2}B{x}^T{V}^T\mathrm{Vs}$ (11)

$+8{\mathrm{\π}}^2\mathrm{\Δ}{t}^{2}A{x}^T{V}^T\mathrm{Vc}-32{\mathrm{\π}}^2\mathrm{\Δ}{t}^2\mathrm{AB}{s}^T{V}^T\mathrm{Vc}$

这里“T”表示向量或者矩阵转置，V是系数矩阵，

$V=\left[\begin{array}{cccc}1& 0& \·\·\·& 0\\ 0& 2& \·\·\·& 0\\ \·& \·& & \·\\ \·& \·& \·\·\·& \·\\ \·& \·& & \·\\ 0& 0& \·\·\·& N\end{array}\right]---\left(12\right)$

则频率的增量估算公式为

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

工业交流电干扰频率的修正公式为

f＝f₀+γΔf                    (14)

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

自适应频率计算算法为：

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

(2)分析确定工业交流电干扰初始频率f₀；

(3)设置k＝0；

(4)设置f＝f₀。通过方程(5)和(6)可以确定出工业交流电干扰余弦函数向量c和正弦函数向量s，由方程(8)和(9)计算出工业交流电干扰正余弦函数的振幅A和B。

(5)由方程(10)和(11)式计算目标函数一阶导数

和二阶导数

(6)由方程(13)式计算频率修正值Δf；

(7)由方程(14)式计算频率修正量f；

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

(9)设置f₀＝f；

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

(11)频率f就是所要求的工业交流电干扰频率f。

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

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

S_i＝x_i-y_i                (15)

这里S_i是消除工业交流电干扰后的地震记录，即地震有效信号。

8)采用通常的方法根据得出消除工业交流电干扰的地震数据绘制消除工业交流电干扰后的地震数据剖面和存储消除工业交流电干扰后的地震数据。

本发明实施例如下：

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

y_i＝D cos 2πf(i+τ)Δt        (16)

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

$D=\sqrt{A^2+B^2}$

$\mathrm{\τ}=\frac{1}{2\mathrm{\πf\Δt}}\mathrm{arctg}(-\frac{A}{B})---\left(17\right)$

和

A＝D cos 2πfτΔt

B＝-D sin 2πfτΔt    (18)

合成余弦函数的振幅D为7.836130E+10，频率f为49.850Hz，延迟τ为2.750，时间采样间隔Δt为1ms。由方程式(18)可以计算出对应余弦函数和正弦函数的振幅A和B分别为5.104585E+10和-5.945430E+10。

图1-2是本方法合成理论数据试算对比图。

图1理论数据试算对比图，每种数据显示五道，(a)是一个实际地震道，其数据的最大值是7.836130E+11，数据道长是5000ms，时间采样间隔是1ms，这里显示了前4000ms。(b)是理论工业交流电干扰记录，由方程(16)制作的合成理论工业交流电干扰记录，(c)是合成理论记录，是(a)与(b)之和，(d)是使用本发明去除工业交流电干扰后的地震道，(e)是使用本发明检测出的工业交流电干扰，其中频率是49.850Hz，正弦函数的振幅是-5.945425E+10，余弦函数的振幅是5.104583E+10，由方程式(17)可以计算出对应方程(16)的振幅D为7.836124E+10，延迟τ为2.750。

图2理论数据频谱对比图，(a)是一个实际地震道频谱，(b)是理论工业交流电干扰记录频谱，(c)是合成理论记录频谱，(d)是使用本发明去除工业交流电干扰后的地震道频谱，(e)是使用本发明检测出的工业交流电干扰频谱。从数据显示和频谱图中可以看出，本发明有效地消除了工业交流电干扰。

图3-4是实际地震数据处理效果对比和频谱对比图。

图3是实际野外炮集数据处理效果对比，显示了两个炮集。(a)是原始炮集数据，炮集数据每炮106道，记录长度5000ms，采样间隔1ms。显然地震数据中存在着很强的工业交流电干扰，必须进行工业交流电干扰压制处理。(b)是陷波处理后炮集数据，陷波器算子长度为500ms，陷波器频带宽度为8Hz，即NH(46，50，54)，NH(146，150，154)和NH(246，250，254)，陷波处理后数据浅层去除了工业交流电干扰，而深层还存在很强的残留工业交流电干扰。而且在地震数据的起始和结尾处存在有很强的边界效应，这种边界效应是不可避免的，任何频率域处理方法都会存在。(c)是时间域单频干扰压制方法消除工业交流电干扰处理后炮集数据，显然有效地消除了地震数据上工业交流电干扰。(d)是本发明工业交流电干扰消除处理后炮集数据，显然有效地消除了地震数据上工业交流电干扰。(e)本发明检测出的工业交流电干扰炮集数据。

图4是实际数据第51炮第81道频谱对比，显示了第51炮第81道浅层(100-2500ms)和深层(2500-4900ms)频谱对比。(a)是原始数据振幅谱，(b)是陷波处理后数据振幅谱，(c)是时间域单频干扰压制方法消除工业交流电干扰处理后振幅谱，(d)是本发明处理后数据振幅谱。对于浅层地震数据，原始数据浅层工业交流电干扰和有效信号能量相当，在陷波处理后数据频谱上，工业交流电干扰得到了消除，但是工业交流电干扰附近有效信号也被消除，因此伤害了有效信号，而本方法一点也没有伤害有效信号。对于深层地震数据，原始数据工业交流电干扰能量比有效信号能量要强得多，在陷波处理后数据频谱上，工业交流电干扰得到了一定的压制，但是还存在很强的残余工业交流电干扰，而本方法有效地消除了工业交流电干扰，一点没有残余工业交流电干扰。总之，陷波法处理后浅层工业交流电有一定的压制，但伤害了工业交流电干扰附近有效信号，因此降低了工业交流电干扰附近有效信号的信噪比；而本方法有效地消除了工业交流电干扰，而并不伤害干扰频率附近有效信号，因此提高了工业交流电干扰频率附近有效信号的信噪比。同时很明显，本方法与时间域单频干扰压制方法一样，有效地消除了工业交流电干扰，且两种算法看不出差异。在DELL工作站上，该数据一炮106道，本发明的运算时间是14s，而同样的机器同样的参数，时间域单频干扰压制方法需要61s，明显优于时间域单频干扰压制方法消除工业交流电干扰。

Claims (2)

1.一种自适应识别和消除地震勘探工业交流电干扰的方法，其特征在于采用如下步骤实现：

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

所述的预处理是指对地震数据置标签、定义观测系统；

2)确定工业交流电干扰的初始频率；

3)计算工业交流电干扰余弦函数和正弦函数的振幅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(1\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(2\right)$

其中

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

c^T＝(cos2πfΔt，cos4πfΔt，Λ，cos2NπfΔt)    (4)

s^T＝(sin2πfΔt，sin4πfΔt，Λ，sin2NπfΔT)    (5)

式中：

A为工业交流电干扰余弦函数的振幅；

B为工业交流电干扰正弦函数的振幅；

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

c为工业交流电干扰余弦函数向量，根据频率、时间采样率计算得到；

s为工业交流电干扰正弦函数向量，根据频率、时间采样率计算得到；

f为工业交流电干扰频率，自适应频率计算算法确定；

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

T为向量转置；

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

4)根据地震数据向量x和工业交流电干扰初始频率f₀，采用以下频率修正公式确定工业交流电干扰频率f，

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

其中f₀是初始频率；γ是迭代频率自适应步长；Q、

和

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

目标函数Q的表达式为

Q＝(x-Ac-Bs)^T(x-Ac-Bs)

                                                (7)

 ＝x^Tx-2Bx^Ts-2Ax^Tc+B²s^Ts+2ABs^Tc+A²c^Tc

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

为

${\▿}_{f}Q=4\mathrm{\π\Δt}(B^2-A^2)s^T\mathrm{Vc}-4\mathrm{\π\Δt}{\mathrm{Bx}}^T\mathrm{Vc}$ (8)

$+4\mathrm{\π\Δt}{\mathrm{Ax}}^T\mathrm{Vs}+4\mathrm{\π\ΔtAB}{(c+s)}^{T}V(c-s)$

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

为

${\▿}_f^{2}Q=8{\mathrm{\π}}^2{\mathrm{\Δt}}^2(B^2-A^2){(c+s)}^T{V}^{T}V(c-s)+8{\mathrm{\π}}^2{\mathrm{\Δt}}^2{\mathrm{Bx}}^T{V}^T\mathrm{Vs}$ (9)

$+8{\mathrm{\π}}^2{\mathrm{\Δt}}^{2}A{x}^T{V}^T\mathrm{Vc}-32{\mathrm{\π}}^2\mathrm{\Δ}{t}^2\mathrm{AB}{s}^T{V}^T\mathrm{Vc}$

这里“T”表示向量或者矩阵转置，V是系数矩阵，

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

初始频率f₀由步骤2)确定；计算第一步中，f为f₀；

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

y_i＝Acos 2πfiΔt+Bsin 2πfiΔt    (10)

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

6)按照以下公式计算地震有效信号：

S_i＝x_i-y_i    (11)

式中：原始地震数据x_i，由野外数据采集得到；参考的工业交流电干扰y_i，由工业交流电干扰计算公式计算得到；S_i是消除工业交流电干扰后的地震有效信号。

2.根据权利要求1所述的自适应识别和消除地震勘探工业交流电干扰的方法，其特征在于步骤2)所述的确定工业交流电干扰的初始频率是指根据地震数据中原始波形数据x_i和它的振幅谱，由原始波形数据的震荡周期和它的振幅谱的最大位置所对应的频率，确定原始数据中工业交流电干扰的初始频率f₀。

Images