CN101907726B - 一种自动识别和消除地震勘探工业电干扰的方法 - Google Patents

Abstract

本发明涉及油田的勘探、开发技术中自动识别和消除地震勘探工业电干扰的方法，确定工业电干扰的初始频率，计算工业电干扰余弦函数和正弦函数构成余弦函数矩阵和正弦函数矩阵，计算确定工业电干扰振幅向量，计算余弦函数自适应减和正余弦函数自适应减确定工业电干扰，得到地震有效信号。本发明不需要确定工业电干扰频率，有效提高计算效率。本发明既可以消除地震数据中由高压输电线产生的工业电干扰，也可以消除地震数据中由周期性震动(如发电机等)产生的干扰。

Description

一种自动识别和消除地震勘探工业电干扰的方法

技术领域

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

背景技术

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

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

自相关和褶积是信号处理中两种最基本、最常用的运算。对于工业电波干扰，我们提出了基于信号分析理论的自相关褶积分析识别和消除工业电干扰的方法。其基本原理是：通过自相关运算和褶积运算工业电干扰余弦函数和正弦函数，通过自适应减方法估算工业电干扰，然后从地震记录中将其减去，达到消除工业电干扰目的。本发明一种自动识别和消除地震勘探工业电干扰的方法，并不直接估算工业电干扰频率参数，因此计算速度快，工业电干扰估算效率高。

发明内容

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

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

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

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

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

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

3)自相关分析确定工业电干扰余弦函数；

步骤3)所述的确定工业电干扰余弦函数，就是计算地震数据自相关函数，进而由地震数据自相关函数确定工业电干扰余弦函数。其计算公式为

$\cos 2\mathrm{\πft}=\frac{R_{\mathrm{xx}}\left(t\right)}{R_{\mathrm{xx}}\left(0\right)}---\left(1\right)$

式中，

f----工业电干扰频率；

t----工业电干扰时间采样，也是地震数据自相关函数时间采样；

R_xx(t)----地震数据自相关函数；

R_xx(0)----地震数据自相关函数零延迟值。

4)褶积自相关分析确定工业电干扰正弦函数；

步骤4)所述的确定工业电干扰正弦函数，就是计算地震数据褶积函数，进而由地震数据自相关函数和褶积函数确定工业电干扰正弦函数。其计算公式为

$\sin 2\mathrm{\πft}=\frac{P_{\mathrm{xx}}\left(0\right)R_{\mathrm{xx}}\left(t\right)-R_{\mathrm{xx}}\left(0\right)P_{\mathrm{xx}}\left(t\right)}{R_{\mathrm{xx}}\left(0\right)\sqrt{R_{\mathrm{xx}}^2\left(0\right)-P_{\mathrm{xx}}^2\left(0\right)}}---\left(2\right)$

式中，

P_xx(t)----地震数据褶积函数；

P_xx(0)----地震数据褶积函数零延迟值。

5)根据工业电干扰余弦函数和正弦函数构成余弦函数矩阵和正弦函数矩阵；

步骤5)所述的根据工业电干扰余弦函数构成余弦函数矩阵，就是由步骤3)自相关分析所确定工业电干扰余弦函数，按照下列方程构成余弦函数矩阵：

$c=\left[\begin{array}{cccc}\cos 2\mathrm{\πf}(1-L)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2-L)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N-L)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}(1-L+1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2-L+1)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N-L+1)\mathrm{\Δt}\\ M& M& M& M\\ \cos 2\mathrm{\πf}(1-1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2-1)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N-1)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}\left(1\right)\mathrm{\Δt}& \cos 2\mathrm{\πf}\left(2\right)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}\left(N\right)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}(1+1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2+1)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N+1)\mathrm{\Δt}\\ M& M& M& M\\ \cos 2\mathrm{\πf}(1+L-1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2+L-1)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N+L-1)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}(1+L)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2+L)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N+L)\mathrm{\Δt}\end{array}\right]---\left(3\right)$

步骤5)所述的根据工业电干扰正弦函数构成正弦函数矩阵，就是由步骤4)自相关褶积分析所确定工业电干扰正弦函数，按照下列方程构成正弦函数矩阵：

$s=\left[\begin{array}{cccc}\sin 2\mathrm{\πf}(1-L)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2-L)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N-L)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}(1-L+1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2-L+1)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N-L+1)\mathrm{\Δt}\\ M& M& M& M\\ \sin 2\mathrm{\πf}(1-1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2-1)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N-1)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}\left(1\right)\mathrm{\Δt}& \sin 2\mathrm{\πf}\left(2\right)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}\left(N\right)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}(1+1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2+1)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N+1)\mathrm{\Δt}\\ M& M& M& M\\ \sin 2\mathrm{\πf}(1+L-1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2+L-1)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N+L-1)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}(1+L)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2+L)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N+L)\mathrm{\Δt}\end{array}\right]---\left(4\right)$

6)计算确定工业电干扰振幅向量

步骤6)所述的计算确定工业电干扰振幅向量，包括工业电余弦函数振幅向量和工业电正余弦函数振幅向量。就是根据余弦函数矩阵和正弦函数矩阵，计算工业电干扰振幅向量。工业电干扰振幅向量计算公式如下：

a＝(CC^T)^-1Cx                (5)

对于工业电余弦函数振幅向量，C就是工业电干扰余弦函数矩阵c，a表示工业电干扰余弦函数系数向量；而对于工业电正余弦函数振幅向量，C就是由工业电干扰余弦函数矩阵c和工业电干扰余弦函数矩阵s按照下式构成，

$C=\left(\begin{array}{c}c\\ s\end{array}\right)---\left(6\right)$

a就是由工业电干扰正弦函数系数向量A和余弦函数系数向量B构成，

$a=\left(\begin{array}{c}A\\ B\end{array}\right)---\left(7\right)$

7)计算确定工业电干扰

步骤7)所述的确定工业电干扰，包括余弦函数自适应减确定工业电干扰和正余弦函数自适应减确定工业电干扰。

余弦函数自适应减确定工业电干扰就是根据工业电干扰余弦函数，采用自适应减方法，计算确定地震数据工业电干扰。余弦函数自适应减工业电干扰可以表示为

y＝C^Ta                        (8)

式中，

y----工业电干扰向量；

C----工业电干扰余弦函数矩阵；

a----工业电干扰余弦函数系数向量；

T----向量或者矩阵转置。

且

y^T＝(y₁，y₂，Λ，y_N)

a^T＝(A_-L，A_-L+1，A_-L+2，Λ，A_-1，A₀，A₁，Λ，A_L-2，A_L-1，A_L)

正余弦函数自适应减确定工业电干扰就是根据工业电干扰正余弦函数，采用自适应减方法，计算确定地震数据工业电干扰。余弦函数自适应减工业电干扰可以表示为

y＝S^TA+c^TB                (9)

式中，

y-----工业电干扰向量；

s----工业电干扰正弦函数矩阵

c----工业电干扰余弦函数矩阵

A----工业电干扰正弦函数系数向量；

B----工业电干扰余弦函数系数向量。

y^T＝(y₁，y₂，Λ，y_N)

A^T＝(A_-L，A_-L+1，A_-L+2，Λ，A_-1，A0，A₁，Λ，A_L-2，A_L-1，A_L)

B^T＝(B_-L，B_-L+1，B_-L+2，Λ，B_-1，B₀，B₁，Λ，B_L-2，B_L-1，B_L)

8)计算确定地震有效信号；

计算确定出工业电干扰之后，由地震数据减去工业电干扰，得到地震有效信号。

S_i＝x_i-y_i                      (10)

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

9)绘制消除工业电干扰后的地震数据剖面和存储消除工业电干扰后的地震数据。

本发明不需要确定工业电干扰频率，可有效提高计算效率。

本发明使用深层地震数据或者初至到达时间之前的地震数据来估计工业电干扰，最有效地估算工业电干扰的能量，可以达到最大限度地压制工业电干扰频率成分，而使该频率分量上的有效波受到的损害最小，提高了该频率分量的信噪比。

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

附图说明

图1理论数据对比

(a)理论工业电数据；

(b)自相关函数；

(c)工业电余弦函数；

(d)工业电振幅向量；

(e)计算工业电；

(f)计算工业电与理论工业电误差

图2合成数据工业电干扰压制效果对比

a合成工业电干扰；

b实际信号；

c合成数据；

d陷频滤波法；

e计算余弦函数；

f计算工业电干扰；

g恢复信号

图3合成数据工业电干扰压制效果频谱对比

a合成工业电干扰；

b实际信号；

c合成数据；

d陷频滤波法；

e计算余弦函数；

f计算工业电干扰；

g恢复信号

图4实际数据工业电干扰压制效果对比

a原始数据；

b时间域工业电干扰压制法；

c自相关分析法；

d自相关分析法检测的工业电干扰

图5实际数据工业电干扰压制效果频谱对比

a原始数据；

b时间域工业电干扰压制法；

c自相关分析法；

d自相关分析法检测的工业电干扰

图6理论数据对比

(a)理论工业电数据；

(b)自相关函数；

(c)褶积函数；

(d)工业电余弦函数；

(e)工业电正弦函数；

(f)工业电振幅向量；

(g)计算工业电；

(h)计算工业电与理论工业电误差

图7合成数据工业电干扰压制效果对比

a合成工业电干扰；

b实际信号；

c合成数据；

d计算自相关函数；

e计算褶积函数；

f计算余弦函数；

g计算正弦函数；

h计算工业电干扰；

i陷频滤波法；

j恢复信号

图8合成数据工业电干扰压制效果频谱对比

a合成工业电干扰；

b实际信号；

c合成数据；

d计算自相关函数；

e计算褶积函数；

f计算余弦函数；

g计算正弦函数；

h计算工业电干扰；

i陷频滤波法；

j恢复信号

图9实际数据工业电干扰压制效果对比

(a)原始数据；

(b)陷频滤波法；

(c)自相关褶积分析法

(d)自相关褶积分析法检测的工业电干扰

图10实际数据第1345道工业电干扰压制效果频谱对比

(a)原始数据；

(b)陷频滤波法；

(c)自相关褶积分析法；

(d)自相关褶积分析法检测的工业电干扰

图11实际数据第1485道工业电干扰压制效果频谱对比

(a)原始数据；

(b)陷频滤波法；

(c)自相关褶积分析法；

(d)自相关褶积分析法检测的工业电干扰

图12实际数据第1645道工业电干扰压制效果频谱对比

(a)原始数据；

(b)陷频滤波法；

(c)自相关褶积分析法；

(d)自相关褶积分析法检测的工业电干扰

具体实施方式

高压输电线或周期性震动(如发电机等)会在地面附近产生周期性工业电干扰。在地震数据采集过程中，如果在地面接收器附近存在高压输电线或周期性震动(如发电机等)，那么地面接收器接收到的地震记录就是地震有效信号和工业电干扰的叠加。本发明的一种自适应识别和消除地震勘探工业电干扰方法就是识别并消除地震记录中的工业电干扰。

本发明首先计算地震数据的褶积函数和自相关函数，然后计算工业电干扰正弦函数和余弦函数；再利用自适应减方法直接计算工业电干扰；最后从地震记录中减去工业电干扰的方法来消除工业电干扰。

基于自相关的工业电干扰余弦函数的自适应减方法消除工业电干扰。本发明包括以下步骤：

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

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

(3)地震记录为地震有效信号和工业电干扰的和，即

x(t)＝S(t)+y(t)                            (1)

式中，

x(t)----地震记录；

S(t)----地震有效信号；

y(t)----工业电干扰；

t----工业电干扰的时间；

(4)确定工业电干扰余弦函数；

根据地震数据自相关函数，计算确定工业电干扰余弦函数。地震数据的自相关函数定义为

$R_{\mathrm{xx}}\left(\mathrm{\τ}\right)={\∫}_0^{T}x\left(t\right)x(t+\mathrm{\τ})\mathrm{dt}---\left(2\right)$

将方程(1)代入(2)，有

R_xx(τ)＝R_SS(τ)+R_Sy(τ)+R_yS(τ)+R_yy(τ)            (3)

其中，

R_xx(τ)----地震记录的自相关；

R_SS(τ)----地震有效信号的自相关；

R_Sy(τ)----地震有效信号和工业电干扰的互相关；

R_yS(τ)----工业电干扰和地震有效信号的互相关；

R_yy(τ)----工业电干扰的自相关；

且

$R_{\mathrm{SS}}\left(\mathrm{\τ}\right)={\∫}_0^{T}S\left(t\right)S(t+\mathrm{\τ})\mathrm{dt}$

$R_{\mathrm{Sy}}\left(\mathrm{\τ}\right)={\∫}_0^{T}S\left(t\right)y(t+\mathrm{\τ})\mathrm{dt}$

$R_{\mathrm{yS}}\left(\mathrm{\τ}\right)={\∫}_0^{T}y\left(t\right)S(t+\mathrm{\τ})\mathrm{dt}$

$R_{\mathrm{yy}}\left(\mathrm{\τ}\right)={\∫}_0^{T}y\left(t\right)y(t+\mathrm{\τ})\mathrm{dt}---\left(4\right)$

假设地震数据与工业电干扰是不相关的，则它们之间的互相关为零，即

R_Sy(τ)＝R_yS(τ)＝0               (5)

将方程(5)代入(3)，有

R_xx(τ)＝R_SS(τ)+R_yy(τ)          (6)

即，一个地震记录之间的自相关就是地震有效信号之间的自相关与工业电干扰之间的自相关之和。从方程(6)可以看出，如果地震有效信号之间的自相关为零，则工业电干扰之间的自相关就是地震记录之间的自相关。即

R_SS(τ)＝0                         (7)

则

R_xx(τ)＝R_yy(τ)                   (8)

由于在地震数据深层，地震有效信号能量与工业电干扰能量相比，要小的多，因此利用深层资料估算地震数据之间的自相关，方程(7)和(8)会近似满足。对于地震数据初至到达时间之前，由于还没有地震有效信号到达，方程(7)理论上绝对满足，这样方程(8)也完全满足。因此为了估算工业电干扰之间的自相关，可以使用地震数据初至到达时间之前的数据来估算，如果初至时间比较小，则可以使用深层地震数据来估算。

为了消除地震记录上的工业电干扰影响，假设工业电干扰的频率、振幅和相位在整个地震记录道内是稳定不变的，且为常数，则可以使用余弦函数来表示工业电干扰。其表达式是：

y(t)＝A cos(2πft+φ)                   (9)

式中，

A----工业电干扰的振幅

f----工业电干扰的频率

φ----工业电干扰的相位

t----工业电干扰的时间

将方程(9)代入方程(4)中，并经过简单运算，有

$R_{\mathrm{yy}}\left(\mathrm{\τ}\right)=\frac{A^{2}T}{2}\cos 2\mathrm{\πf\τ}+\frac{A^2}{4\mathrm{\ω}}[\sin (2\mathrm{\ωT}+2\mathrm{\φ}+\mathrm{\ω\τ})-\sin (2\mathrm{\φ}+\mathrm{\ω\τ})]---\left(10\right)$

如果T是工业电干扰的周期(或者周期的整数倍)，则有

$R_{\mathrm{yy}}\left(\mathrm{\τ}\right)=\frac{A^{2}T}{2}\cos 2\mathrm{\πf\τ}---\left(11\right)$

在方程(11)中，令τ＝0，则有

$R_{\mathrm{yy}}\left(0\right)=\frac{A^{2}T}{2}---\left(12\right)$

将方程(12)代入方程(11)，并把τ换成t，有

$\cos 2\mathrm{\πft}=\frac{R_{\mathrm{yy}}\left(t\right)}{R_{\mathrm{yy}}\left(0\right)}---\left(13\right)$

将方程(13)代入方程(8)中，有

$\cos 2\mathrm{\πft}=\frac{R_{\mathrm{xx}}\left(t\right)}{R_{\mathrm{xx}}\left(0\right)}---\left(14\right)$

方程(14)就是利用地震数据自相关函数计算工业电干扰余弦函数公式。这样对给定的地震数据，仅仅需要计算一次相关，就可以确定工业电干扰的余弦函数。显然运算次数要远远小于各种工业电干扰估算运算次数。对于几十万甚至几千万道的叠前3D地震数据来说，运算效率是相当可观的。

(5)自适应减方法确定工业电干扰函数；

估算出工业电干扰余弦函数之后，采样基于工业电干扰余弦函数的自适应减方法估算工业电干扰，这样可以省去估算工业电干扰的频率、振幅和相位三个参数。基于工业电干扰余弦函数的自适应减工业电干扰可以表示为：

$y_i=\underset{k=-L}{\overset{L}{\mathrm{\Σ}}}{A}_k\cos \mathrm{\ω}(i+k)\mathrm{\Δt}---\left(15\right)$

式中，

A_k----工业电干扰余弦函数的第k个系数；

Δt----地震数据时间采样间隔；

i----工业电干扰的时间样点序号；

k----工业电干扰的系数顺序号；

2L+1----工业电干扰系数总数。

把方程(15)用向量可以表示为

y＝C^Ta                   (16)

式中，

y----工业电干扰向量；

C----工业电干扰余弦函数矩阵；

a----工业电干扰余弦函数系数向量；

T----向量或者矩阵转置。

且

y^T＝(y₁，y₂，Λ，y_N)

a^T＝(A_-L，A_-L+1，A_-L+2，Λ，A_-1，A₀，A₁，Λ，A_L-2，A_L-1，A_L)

$C=\left[\begin{array}{cccc}\cos \mathrm{\ω}(1-L)\mathrm{\Δt}& \cos \mathrm{\ω}(2-L)\mathrm{\Δt}& \mathrm{\Λ}& \cos \mathrm{\ω}(N-L)\mathrm{\Δt}\\ \cos \mathrm{\ω}(1-L+1)\mathrm{\Δt}& \cos \mathrm{\ω}(2-L+1)\mathrm{\Δt}& \mathrm{\Λ}& \cos \mathrm{\ω}(N-L+1)\mathrm{\Δt}\\ M& M& M& M\\ \cos \mathrm{\ω}(1-1)\mathrm{\Δt}& \cos \mathrm{\ω}(2-1)\mathrm{\Δt}& \mathrm{\Λ}& \cos \mathrm{\ω}(N-1)\mathrm{\Δt}\\ \cos \mathrm{\ω}\left(1\right)\mathrm{\Δt}& \cos \mathrm{\ω}\left(2\right)\mathrm{\Δt}& \mathrm{\Λ}& \cos \mathrm{\ω}\left(N\right)\mathrm{\Δt}\\ \cos \mathrm{\ω}(1+1)\mathrm{\Δt}& \cos \mathrm{\ω}(2+1)\mathrm{\Δt}& \mathrm{\Λ}& \cos \mathrm{\ω}(N+1)\mathrm{\Δt}\\ M& M& M& M\\ \cos \mathrm{\ω}(1+L-1)\mathrm{\Δt}& \cos \mathrm{\ω}(2+L-1)\mathrm{\Δt}& \mathrm{\Λ}& \cos \mathrm{\ω}(N+L-1)\mathrm{\Δt}\\ \cos \mathrm{\ω}(1+L)\mathrm{\Δt}& \cos \mathrm{\ω}(2+L)\mathrm{\Δt}& \mathrm{\Λ}& \cos \mathrm{\ω}(N+L)\mathrm{\Δt}\end{array}\right]---\left(17\right)$

为了确定工业电干扰余弦函数系数向量a，建立如下目标函数

Q＝(x-y)^T(x-y)

＝x^Tx-2x^TC^Ta+a^TCC^Ta          (18)

在方程(18)中，令

$\frac{\∂Q}{\∂a}=-{2x}^T{C}^T+2a^T{\mathrm{CC}}^T=0---\left(19\right)$

有

由方程(19)，得到计算工业电干扰振幅向量的公式

a＝(CC^T)^-1Cx                 (20)

这样，对于工业电干扰余弦函数，首先由方程(20)求解出工业电干扰的系数向量a；然后由方程(16)可以计算出工业电干扰；再从原始地震道中减去估算的工业电干扰，就得到了去除工业电干扰的地震记录。

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

S_i＝x_i-y_i                     (21)

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

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

首先在理论工业电数据上进行计算和对比。理论工业电数据采用方程(22)生成。

y(t)＝A cos 2πft+Bsin 2πft              (22)

这里使用的参数是：工业电余弦函数和正弦函数的振幅A和B分别是2.543和4.816；工业电频率f是50.135Hz；时间采样间隔Δt是1ms，样本个数N是500，计算的理论数据如图1(a)所示。图1(b)是由公式(4)计算的自相关函数，其中T取值为80ms，图1(c)是由公式(14)计算的工业电余弦函数。对于L＝10，由公式(20)计算的工业电振幅向量如图1(d)所示。图1(e)是由公式(17)计算的工业电。我们计算了图1(a)理论工业电数据与图1(e)计算工业电之间的互相关，其互相关值为0.9988，即理论工业电数据与计算工业电数据之间完全一致。图1(f)是图1(a)理论工业电数据与图1(e)计算工业电之差，从图中可以看出，仅仅在数据的起始20ms和终止20ms部分，由于计算的边界效应产生了数据之间的较大误差之外，在数据中间部分误差很小，几乎为零。因此表明方法在理论上是有效的。

然后使用理论数据来说明方法的有效性。理论工业电数据采用方程(22)生成。这里使用的参数是：工业电余弦函数和正弦函数的振幅A和B分别是6.541和4.811；工业电频率f是50.254Hz；时间采样间隔Δt是1ms，样本个数N是5000，生成工业电干扰数据，一道的数据十道显示如图2(a)所示。信号采用一段实际地震数据，其最大值为工业电5倍，如图2(b)所示，生成的理论数据如图2(c)所示，陷频滤波法处理后的信号如图2(d)所示(陷频滤波参数为：陷频滤波器算子长度为500ms，频带宽度8Hz，即NH(46，50，54))。通过本方法计算出的工业电干扰余弦函数如图2(e)所示，计算生成的工业电干扰数据如图2(f)所示，本方法处理后恢复的信号如图2(g)所示。它们对应的频谱分别如图3(a)、(b)、(c)、(d)、(e)、(f)和(g)所示。从数据上看，陷频滤波法可以消除工业电干扰，但是在920ms处出现一些边界效应，从频谱上看，陷频滤波法在消除工业电干扰的同时，也严重损害工业电附近的有效信号。而本方法有效的识别并消除了工业电干扰，且没有伤害工业电干扰频率附近信号的频率成分，因此有效的提高了工业电干扰附近信号的信噪比。

实际数据是一个野外炮集数据，有180道，数据时间采样间隔是2ms，数据记录长度是6000ms。图4是炮集数据对比，图5是炮集数据第81道频谱对比：a是原始数据，b是时间域工业电干扰压制法，c是自相关分析法；d是自相关分析法检测的工业电干扰。显然原始数据中包含着很强的工业电干扰，从频谱中也可以清楚地看到工业电干扰。从数据以及频谱中可以看到时间域工业电干扰压制法和自相关分析法都非常有效地消除了地震数据上的工业电干扰，而且两种方法效果差异不大。本方法大大提高了运算效率，节省运算时间，更加适用于海量地震数据处理需要。因此是消除工业电干扰的最有效方法。

基于褶积和自相关的工业电干扰余弦函数的自适应减方法消除工业电干扰。本发明包括以下步骤：

(1)采用如实施例1的步骤(1)----步骤(4)处理。

(2)确定工业电干扰正弦函数；

根据地震数据自相关函数和褶积函数，计算确定工业电干扰正弦函数。地震数据的褶积定义为

$P_{\mathrm{xx}}\left(\mathrm{\τ}\right)={\∫}_0^{T}x\left(t\right)x(t-\mathrm{\τ})\mathrm{dt}---\left(23\right)$

将方程(1)代入(23)，有

P_xx(τ)＝P_SS(τ)+P_Sy(τ)+P_yS(τ)+P_yy(τ)            (24)

其中，

P_xx(τ)----地震记录的褶积；

P_SS(τ)----地震有效信号的褶积；

P_Sy(τ)----地震有效信号和工业电干扰的褶积；

P_yS(τ)----工业电干扰和地震有效信号的褶积；

P_yy(τ)----工业电干扰的褶积；

且

$P_{\mathrm{SS}}\left(\mathrm{\τ}\right)={\∫}_0^{T}S\left(t\right)S(t-\mathrm{\τ})\mathrm{dt}$

$P_{\mathrm{Sy}}\left(\mathrm{\τ}\right)={\∫}_0^{T}S\left(t\right)y(t-\mathrm{\τ})\mathrm{dt}$

$P_{\mathrm{yS}}\left(\mathrm{\τ}\right)={\∫}_0^{T}y\left(t\right)S(t-\mathrm{\τ})\mathrm{dt}$

$P_{\mathrm{yy}}\left(\mathrm{\τ}\right)={\∫}_0^{T}y\left(t\right)y(t-\mathrm{\τ})\mathrm{dt}---\left(25\right)$

同样假设地震数据与工业电干扰是不相关的，则它们之间的互相关为零，即

P_Sy(τ)＝P_yS(τ)＝0                     (26)

将方程(26)代入(24)，有

P_xx(τ)＝P_SS(τ)+P_yy(τ)                (27)

即，一个地震记录之间的褶积就是地震有效信号之间的褶积与工业电干扰之间的褶积之和，从方程(27)可以看出，如果地震有效信号之间的褶积为零，则工业电干扰之间的褶积就是地震记录之间的褶积。即

P_SS(τ)＝0                              (28)

则

P_xx(τ)＝P_yy(τ)                        (29)

同样由于在地震数据深层，地震有效信号能量与工业电干扰能量相比，要小的多，因此利用深层资料估算地震数据之间的褶积，方程(26)、(28)和(29)会近似满足。对于地震数据初至到达时间之前，由于还没有地震有效信号到达，方程(28)理论上绝对满足，这样方程(29)也完全满足。因此为了估算工业电干扰之间的褶积，可以使用地震数据初至到达时间之前的数据来估算，如果初至时间比较小，则可以使用深层地震数据来估算。

将方程(9)代入方程(24)中，并经过简单运算，有

$P_{\mathrm{yy}}\left(\mathrm{\τ}\right)=\frac{A^{2}T}{2}\cos (2\mathrm{\πf\τ}+2\mathrm{\φ})+\frac{A^2}{4\mathrm{\ω}}[\sin (2\mathrm{\ωT}-\mathrm{\ω\τ})+\sin \mathrm{\ωt}]---\left(29\right)$

如果T是工业电干扰的周期(或者周期的整数倍)，则有

$P_{\mathrm{yy}}\left(\mathrm{\τ}\right)=\frac{A^{2}T}{2}\cos (2\mathrm{\πf\τ}+2\mathrm{\φ})---\left(30\right)$

在方程(30)中，令τ＝0，则有

P_yy(0)＝R_yy(0)cos 2φ                    (31)

将方程(31)、(11)、(12)和(13)分别代入方程(30)，并把τ换成t，有

$\sin 2\mathrm{\πft}=\frac{P_{\mathrm{yy}}\left(0\right)R_{\mathrm{yy}}\left(t\right)-R_{\mathrm{yy}}\left(0\right)P_{\mathrm{yy}}\left(t\right)}{R_{\mathrm{yy}}\left(0\right)\sqrt{R_{\mathrm{yy}}^2\left(0\right)-P_{\mathrm{yy}}^2\left(0\right)}}---\left(32\right)$

将方程(8)和(28)代入方程(32)中，有

$\sin 2\mathrm{\πft}=\frac{P_{\mathrm{xx}}\left(0\right)R_{\mathrm{xx}}\left(t\right)-R_{\mathrm{xx}}\left(0\right)P_{\mathrm{xx}}\left(t\right)}{R_{\mathrm{xx}}\left(0\right)\sqrt{R_{\mathrm{xx}}^2\left(0\right)-P_{\mathrm{xx}}^2\left(0\right)}}---\left(33\right)$

方程(33)就是由地震数据自相关函数和褶积函数计算工业电干扰正弦函数公式。这样对给定的地震数据，仅仅需要计算一次褶积和一次相关，就可以确定工业电干扰的正弦函数。显然运算次数要远远小于各种工业电干扰估算运算次数。对于几十万甚至几千万道的叠前3D地震数据来说，运算效率是相当可观的。

(3)自适应减方法确定工业电干扰函数；

基于工业电干扰正余弦函数的自适应减方法。估算出工业电干扰余弦函数和正弦函数之后，采样基于工业电干扰正余弦函数的自适应减方法估算工业电干扰，这样可以省去估算工业电干扰的频率、振幅和相位三个参数。基于工业电干扰正余弦函数的自适应减工业电干扰可以表示为：

$y_i=\underset{k=-L}{\overset{L}{\mathrm{\Σ}}}{A}_k\sin 2\mathrm{\πf}(i+k)\mathrm{\Δt}+\underset{k=-L}{\overset{L}{\mathrm{\Σ}}}{B}_k\cos 2\mathrm{\πf}(i+k)\mathrm{\Δt}---\left(34\right)$

式中，

A_k----工业电干扰正弦函数的第k个系数；

B_k----工业电干扰余弦函数的第k个系数

f----工业电干扰的频率

Δt----地震数据时间采样间隔；

i----工业电干扰的时间样点序号

k----工业电干扰的系数顺序号

2L+1----工业电干扰系数总数。

把方程(34)用向量可以表示为

y＝s^TA+c^TB                    (35)

式中，

y-----工业电干扰向量；

s----工业电干扰正弦函数矩阵

c----工业电干扰余弦函数矩阵

A----工业电干扰正弦函数系数向量；

B----工业电干扰余弦函数系数向量。

y^T＝(y₁，y₂，Λ，y_N)

A^T＝(A_-L，A_-L+1，A_-L+2，Λ，A_-1，A₀，A₁，Λ，A_L-2，A_L-1，A_L)

B^T＝(B_-L，B_-L+1，B_-L+2，Λ，B_-1，B₀，B₁，Λ，B_L-2，B_L-1，B_L)

$s=\left[\begin{array}{cccc}\sin 2\mathrm{\πf}(1-L)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2-L)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N-L)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}(1-L+1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2-L+1)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N-L+1)\mathrm{\Δt}\\ M& M& M& M\\ \sin 2\mathrm{\πf}(1-1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2-1)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N-1)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}\left(1\right)\mathrm{\Δt}& \sin 2\mathrm{\πf}\left(2\right)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}\left(N\right)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}(1+1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2+1)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N+1)\mathrm{\Δt}\\ M& M& M& M\\ \sin 2\mathrm{\πf}(1+L-1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2+L-1)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N+L-1)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}(1+L)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2+L)\mathrm{\Δt}& \mathrm{\Λ}& \sin 2\mathrm{\πf}(N+L)\mathrm{\Δt}\end{array}\right]$

$c=\left[\begin{array}{cccc}\cos 2\mathrm{\πf}(1-L)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2-L)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N-L)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}(1-L+1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2-L+1)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N-L+1)\mathrm{\Δt}\\ M& M& M& M\\ \cos 2\mathrm{\πf}(1-1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2-1)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N-1)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}\left(1\right)\mathrm{\Δt}& \cos 2\mathrm{\πf}\left(2\right)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}\left(N\right)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}(1+1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2+1)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N+1)\mathrm{\Δt}\\ M& M& M& M\\ \cos 2\mathrm{\πf}(1+L-1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2+L-1)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N+L-1)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}(1+L)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2+L)\mathrm{\Δt}& \mathrm{\Λ}& \cos 2\mathrm{\πf}(N+L)\mathrm{\Δt}\end{array}\right]---\left(36\right)$

这里符号“T”表示向量或者矩阵转置。为了确定工业电干扰余弦函数系数向量a，建立如下目标函数

Q＝(x-y)^T(x-y)

＝x^Tx-2x^Ts^TA-2x^Tc^TB

+A^Tss^TA+2A^Tsc^TB+B^Tcc^TB    (37)

在方程(37)中，令

$\frac{\∂Q}{\∂A}=-2x^T{s}^T+2A^T{\mathrm{ss}}^T+2B^T{\mathrm{cs}}^T=0---\left(38\right)$

和

$\frac{\∂Q}{\∂B}=-2x^T{c}^T+2A^T{\mathrm{sc}}^T+2B^T{\mathrm{cc}}^T=0---\left(39\right)$

由方程(38)和(39)，得到计算工业电干扰振幅向量的公式

$\left(\begin{array}{c}A\\ B\end{array}\right){\left[\left(\begin{array}{c}c\\ s\end{array}\right){\left(\begin{array}{cc}c& s\end{array}\right)}^T\right]}^{-1}\left(\begin{array}{c}c\\ s\end{array}\right)x---\left(40\right)$

令

$a=\left(\begin{array}{c}A\\ B\end{array}\right),C=\left(\begin{array}{c}c\\ s\end{array}\right)---\left(41\right)$

则有

a＝(CC^T)^-1Cx             (42)

这与方程(20)完全一致。这样，对于工业电干扰正余弦函数，首先由方程(42)求解出工业电干扰的系数向量A和B；然后由方程(35)可以计算出工业电干扰；再从原始地震道中减去估算的工业电干扰，就得到了去除工业电干扰的地震记录。

(4)采用如实施例1的步骤(6)----步骤(7)处理。

首先在理论工业电数据上进行计算和对比。理论工业电数据采用方程(22)生成。这里使用的参数是：工业电余弦函数和正弦函数的振幅A和B分别是2.543和4.816；工业电频率f是50.135Hz；时间采样间隔Δt是1ms，样本个数N是500，计算的理论数据如图6(a)所示；图6(b)是由公式(4)计算的自相关函数，其中T取值为80ms；图6(c)是由公式(23)计算的褶积函数，其中T取值为80ms；图6(d)是由公式(14)计算的工业电余弦函数；图6(e)是由公式(33)计算的工业电正弦函数。对于L＝20，由公式(42)计算的工业电振幅向量A和B如图6(f)所示。图6(g)是由公式(35)计算的工业电。我们计算了图6(a)理论工业电数据与图6(g)计算工业电之间的互相关，其互相关值为0.9998，即理论工业电数据与计算工业电数据之间几乎完全一致。图6(h)是图1(a)理论工业电数据与图6(g)计算工业电之差，从图中可以看出，仅仅在数据的起始5ms部分，由于计算的边界效应产生了数据之间的较大误差之外，在数据中间部分误差很小，几乎为零。因此表明方法在理论上是有效的。

然后使用理论数据来说明方法的有效性。理论工业电数据采用方程(22)生成。这里使用的参数是：工业电余弦函数和正弦函数的振幅A和B分别是6.541和4.811；工业电频率f是50.254Hz；时间采样间隔Δt是1ms，样本个数N是5000，生成工业电干扰数据，一道的数据十道显示如图7(a)所示。信号采用一段实际地震数据，其最大值为工业电5倍，如图7(b)所示，生成的理论数据如图7(c)所示，陷频滤波法处理后的信号如图7(d)所示(陷频滤波参数为：陷频滤波器算子长度为500ms，频带宽度8Hz，即NH(46，50，54))。通过本方法计算出的工业电干扰余弦函数如图7(e)所示，计算生成的工业电干扰数据如图7(f)所示，本方法处理后恢复的信号如图7(g)所示。它们对应的频谱分别如图8(a)、(b)、(c)、(d)、(e)、(f)和(g)所示。从数据上看，陷频滤波法可以消除工业电干扰，但是在920ms处出现一些边界效应，从频谱上看，陷频滤波法在消除工业电干扰的同时，也严重损害工业电附近的有效信号。而本方法有效的识别并消除了工业电干扰，且没有伤害工业电干扰频率附近信号的频率成分，因此有效的提高了工业电干扰附近信号的信噪比。

实际数据是一个野外炮集数据，有180道，数据时间采样间隔是2ms，数据记录长度是6000ms。我们显示了900-4100ms。图9是炮集数据对比，图10是炮集数据第81道频谱对比：a是原始数据，b是时间域工业电干扰压制法，c是自相关分析法；d是自相关分析法检测的工业电干扰。显然原始数据中包含着很强的工业电干扰，从频谱中也可以清楚地看到工业电干扰。从数据以及频谱中可以看到时间域工业电干扰压制法和自相关分析法都非常有效地消除了地震数据上的工业电干扰，而且两种方法效果差异不大。本方法大大提高了运算效率，节省运算时间，更加适用于海量地震数据处理需要。因此是消除工业电干扰的最有效方法。

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

Claims (3)

1.一种自动识别和消除地震勘探工业电干扰的方法，其特征是包括以下步骤：

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

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

3)采用下式计算工业电干扰余弦函数和正弦函数：

$\cos 2\mathrm{\πft}=\frac{R_{\mathrm{xx}}\left(t\right)}{R_{\mathrm{xx}}\left(0\right)}---\left(1\right)$

式中，

f     ----工业电干扰频率；

t     ----工业电干扰时间采样，也是地震数据自相关函数时间采样；

R_xx(t)----地震数据自相关函数；

R_xx(0)----地震数据自相关函数零延迟值；

4)采用下式计算工业电干扰正弦函数：

$\sin 2\mathrm{\πft}=\frac{P_{\mathrm{xx}}\left(0\right)R_{\mathrm{xx}}\left(t\right)-R_{\mathrm{xx}}\left(0\right)R_{\mathrm{xx}}\left(t\right)}{R_{\mathrm{xx}}\left(0\right)\sqrt{R_{\mathrm{xx}}^2\left(0\right)-P_{\mathrm{xx}}^2\left(0\right)}}---\left(2\right)$

式中，

P_xx(t)----地震数据褶积函数；

P_xx(0)----地震数据褶积函数零延迟值；

5)根据工业电干扰余弦函数和正弦函数构成余弦函数矩阵和正弦函数矩阵；

所述的余弦函数矩阵按照下列方程构成：

$c=\left[\begin{array}{cccc}\cos 2\mathrm{\πf}(1-L)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2-L)\mathrm{\Δt}& \·\·\·& \cos 2\mathrm{\πf}(N-L)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}(1-L+1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2-L+1)\mathrm{\Δt}& \·\·\·& \cos 2\mathrm{\πf}(N-L+1)\mathrm{\Δt}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \cos 2\mathrm{\πf}(1-1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2-1)\mathrm{\Δt}& \·\·\·& \cos 2\mathrm{\πf}(N-1)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}\left(1\right)\mathrm{\Δt}& \cos 2\mathrm{\πf}\left(2\right)\mathrm{\Δt}& \·\·\·& \cos 2\mathrm{\πf}\left(N\right)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}(1+1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2+1)\mathrm{\Δt}& \·\·\·& \cos 2\mathrm{\πf}(N+1)\mathrm{\Δt}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \cos 2\mathrm{\πf}(1+L-1)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2+L-1)\mathrm{\Δt}& \·\·\·& \cos 2\mathrm{\πf}(N+L-1)\mathrm{\Δt}\\ \cos 2\mathrm{\πf}(1+L)\mathrm{\Δt}& \cos 2\mathrm{\πf}(2+L)\mathrm{\Δt}& \·\·\·& \cos 2\mathrm{\πf}(N+L)\mathrm{\Δt}\end{array}\right]---\left(3\right)$

正弦函数矩阵按照下列方程构成：

$s=\left[\begin{array}{cccc}\sin 2\mathrm{\πf}(1-L)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2-L)\mathrm{\Δt}& \·\·\·& \sin 2\mathrm{\πf}(N-L)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}(1-L+1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2-L+1)\mathrm{\Δt}& \·\·\·& \sin 2\mathrm{\πf}(N-L+1)\mathrm{\Δt}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \sin 2\mathrm{\πf}(1-1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2-1)\mathrm{\Δt}& \·\·\·& \sin 2\mathrm{\πf}(N-1)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}\left(1\right)\mathrm{\Δt}& \sin 2\mathrm{\πf}\left(2\right)\mathrm{\Δt}& \·\·\·& \sin 2\mathrm{\πf}\left(N\right)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}(1+1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2+1)\mathrm{\Δt}& \·\·\·& \sin 2\mathrm{\πf}(N+1)\mathrm{\Δt}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \sin 2\mathrm{\πf}(1+L-1)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2+L-1)\mathrm{\Δt}& \·\·\·& \sin 2\mathrm{\πf}(N+L-1)\mathrm{\Δt}\\ \sin 2\mathrm{\πf}(1+L)\mathrm{\Δt}& \sin 2\mathrm{\πf}(2+L)\mathrm{\Δt}& \·\·\·& \sin 2\mathrm{\πf}(N+L)\mathrm{\Δt}\end{array}\right]$

6)按照如下公式计算确定工业电干扰振幅向量：

a＝(CC^T)^-1Cx       (5)

对于工业电余弦函数振幅向量，C是工业电干扰余弦函数矩阵c，a表示工业电干扰余弦函数系数向量；

对于工业电正余弦函数振幅向量，C是由工业电干扰余弦函数矩阵c和工业电干扰余弦函数矩阵s按照下式构成，

$C=\left(\begin{array}{c}c\\ s\end{array}\right)---\left(6\right)$

a是由工业电干扰正弦函数系数向量A和余弦函数系数向量B构成，

$a=\left(\begin{array}{c}A\\ B\end{array}\right);---\left(7\right)$

7)采用下式计算余弦函数自适应减和正余弦函数自适应减确定工业电干扰：

余弦函数计算：

y＝C^Ta             (8)

式中，

y  ----工业电干扰向量；

C  ----工业电干扰余弦函数矩阵；

a  ----工业电干扰余弦函数系数向量；

T  ----向量或者矩阵转置。

且

y^T＝(y₁，y₂，…，y_N)

a^T＝(A_-L，A_-L+1，A_-L+2，…，A_-1，A₀，A₁，…，A_L-2，A_L-1，A_L)

正余弦函数计算：

y＝s^TA+c^TB                (9)

式中，

y  ----- 工业电干扰向量；

s  ----  工业电干扰正弦函数矩阵

c  ----  工业电干扰余弦函数矩阵

A  ----  工业电干扰正弦函数系数向量；

B  ----  工业电干扰余弦函数系数向量。

y^T＝(y₁，y₂，…，y_N)

A^T＝(A_-L，A_-L+1，A_-L+2，…，A_-1，A₀，A₁，…，A_L-2，A_L-1，A_L)

B^T＝(B_-L，B_-L+1，B_-L+2，…，B_-1，B₀，B₁，…，B_L-2，B_L-1，B_L)；

8)采用下述公式计算地震有效信号：

S_i＝x_i-y_i                 (10)

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

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

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

Images