CN101551465A - 一种自适应识别和消除地震勘探单频干扰的方法 - Google Patents

Abstract

本发明涉及野外地震数据采集中自适应识别和消除地震勘探单频干扰的方法。包括以下步骤：确定单频干扰波的初始频率；计算单频干扰波余弦函数和正弦函数的振幅；根据单频干扰波初始频率确定单频干扰波频率；计算单频干扰波和地震有效信号：得出消除单频干扰波的地震数据绘制消除单频干扰波后的地震数据剖面和存储消除单频干扰波后的地震数据剖面。本发明克服了频率域压制单频波的缺点，而且比时间域单频干扰波压制方法运算速度快，能够有效的消除地震记录上的单频波，保留了时间域单频干扰波压制方法不损害有效波的特点。

Description

一种自适应识别和消除地震勘探单频干扰的方法

技术领域

本发明涉及油田的勘探、开发、开采技术，具体是一种在野外地震数据采集中自适应识别和消除地震勘探单频干扰的方法。

背景技术

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

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

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

发明内容

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

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

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

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

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

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

3)计算单频干扰波余弦函数和正弦函数的振幅A和B，其计算公式为：

$A=\frac{\mathrm{dc}-\mathrm{eb}}{\mathrm{ac}-b^2}$

$B=\frac{\mathrm{db}-\mathrm{ea}}{b^2-\mathrm{ac}}$

其中

$a=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\cos }^{2}2\mathrm{\πfi\Δt}$

$b=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\cos 2\mathrm{\πfi\Δ}t\sin 2\mathrm{\πfi\Δt}$

$c=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\sin }^{2}2\mathrm{\πfi\Δt}$

$d=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\cos 2\mathrm{\πfi\Δt}$

$e=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\sin 2\mathrm{\πfi\Δt}$

式中：

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

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

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

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

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

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

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

4)根据地震数据S_i和单频干扰波初始频率f₀，采用以下频率修正公式确定单频干扰波频率f，

f＝f₀+Δf

其中：

f₀是初始频率，Δf是自适应频率修正值，

$\mathrm{\Δf}=\frac{2(\mathrm{VA}-\mathrm{UB})}{N(A^2+B^2)}$

$U=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\cos 2\mathrm{\πfi\Δ}t\sin \frac{2\mathrm{\πi}}{N}$

$V=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\sin 2\mathrm{\πfi\Δ}t\sin \frac{2\mathrm{\πi}}{N}$

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

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

y_i＝Acos2πfiΔt+Bsin2πfiΔt；

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

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

x_i＝s_i-y_i

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

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

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

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

附图说明

图1是本发明理论数据试算对比图，每种数据显示十道，(a)一个实际地震道，(b)理论单频干扰波记录，(c)合成理论记录，(d)使用本发明去除单频干扰波后的地震道，(e)使用本发明检测出的单频干扰波。

图2是本发明理论数据频谱对比图，(a)实际地震道频谱，(b)理论单频干扰波记录频谱，(c)合成理论记录频谱，(d)使用发明去除单频干扰波后的地震道频谱，(e)使用本发明测出的单频干扰波频谱。

图3是本发明VSP数据三分量数据处理对比图。(a)原始VSP三分量数据，(b)VSP软件单频干扰波消除之后的VSP三分量数据，(c)本发明单频干扰波消除之后的VSP三分量数据。

图4是VSP数据三分量频谱对比图。(a)原始分量的频谱，(b)VSP软件单频干扰波消除方法处理后的频谱，(c)本发明单频干扰波消除方法处理后的频谱。

图5是地震炮集数据处理对比图。(a)地震数据，(b)时间域单频干扰波压制方法消除单频干扰波之后的地震数据，(c)本发明单频干扰波消除之后的地震数据。

图6是频谱对比图。(a)原始分量的频谱，(b)时间域单频干扰波压制方法消除单频干扰波之后的频谱，(c)本发明单频干扰波消除方法处理后的频谱。

具体实施方式

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

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

本发明包括以下步骤：

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

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

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

S_i＝x_i+y_i              (1)

式中，

S_i-----表示地震记录

x_i-----表示地震有效信号

y_i-----表示单频干扰波

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

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

本发明使用余弦函数和同频率不同振幅的正弦函数之和来表示单频干扰波。其表达式是：

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容易和快速得多。

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

对于第k次迭代，建立目标函数：

$Q_k=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[S_i-y_i]}^2$

$=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{[S_i-A\cos 2\mathrm{\π}{f}_k\mathrm{i\Δt}-B\sin 2\mathrm{\π}{f}_k\mathrm{i\Δt}]}^2\→\min ---\left(3\right)$

式中S_i是原始地震记录。对于第k次迭代给定的频率f_k，为了振幅A和B的确定，由

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

得出单频干扰波的振幅A和B分别为

$A=\frac{\mathrm{dc}-\mathrm{eb}}{\mathrm{ac}-b^2}$

$B=\frac{\mathrm{db}-\mathrm{ea}}{b^2-\mathrm{ac}}---\left(5\right)$

其中

$a=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\cos }^{2}2\mathrm{\π}{f}_k\mathrm{i\Δt}$

$b=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\cos 2\mathrm{\π}{f}_k\mathrm{i\Δ}t\sin 2\mathrm{\π}{f}_k\mathrm{i\Δt}$

$c=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\sin }^{2}2\mathrm{\π}{f}_k\mathrm{i\Δt}$

$d=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\cos 2\mathrm{\π}{f}_k\mathrm{i\Δt}$

$e=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\sin 2\mathrm{\π}{f}_k\mathrm{i\Δt}---\left(6\right)$

这样，对于一个给定频率f_k，通过方程(6)可以确定出系数a、b、c、d和e，然后由方程(5)可以求解出单频干扰波的振幅A和B。

5)自适应频率计算算法。我们采用自适应频率计算算法估算单频干扰波的频率。首先定义

$U=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\cos 2\mathrm{\π}{f}_k\mathrm{i\Δ}t\sin \frac{2\mathrm{\πi}}{N}$

$V=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\sin 2\mathrm{\π}{f}_k\mathrm{i\Δ}t\sin \frac{2\mathrm{\πi}}{N}---\left(7\right)$

则频率的增量估算公式为

$\mathrm{\Δf}=\frac{2(\mathrm{VA}-\mathrm{UB})}{N(A^2+B^2)}---\left(8\right)$

单频干扰波频率的修正公式为

f_k+1＝f_k+Δf    (9)

这里f₀是初始频率，k是迭代次数，Δf是频率修正值。

自适应频率计算算法为：

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

(2)分析确定扫描初始频率f₀；

(3)设置k＝0；

(4)方程(6)可以确定出系数a、b、c、d和e，然后由方程(5)可以求解出单频干扰波的振幅A和B；

(5)由方程(7)式计算数值U和V；

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

(7)由方程(9)式计算频率修正量f_k+1；

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

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

(10)频率f_k+1就是所要求的单频波频率f。

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

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

x_i＝s_i-y_i                   (10)

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

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

本发明实施例。

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

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

这里，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(12\right)$

和

A＝Dcos2πfτΔt

B＝Dsin2πfτΔt            (13)

合成余弦函数的振幅D为3.726241E+12，频率f为50.450Hz，延迟τ为3.150，时间采样间隔Δt为1ms。由方程式(13)可以计算出对应余弦函数和正弦函数的振幅A和B分别为2.017972E+12和-3.132517E+12。

图1-2是本方法合成理论数据试算对比图。图1理论数据试算对比图，每种数据显示十道，(a)是一个实际地震道，其数据的最大值是3.726241E+13，数据道长是5000ms，时间采样间隔是1ms，(b)是理论单频干扰波记录，由方程(11)制作的合成理论单频干扰波记录，(c)是合成理论记录，是实际地震道(a)与理论单频干扰波记录(b)之和，(d)是使用本发明去除单频干扰波后的地震道，(e)是使用本发明检测出的单频干扰波，其中频率是50.450Hz，正弦函数的振幅是-3.132525E+12，余弦函数的振幅是2.017978E+12，由方程式(12)可以计算出对应方程(11)的振幅D为3.726251E+12，延迟τ为3.150。图2理论数据频谱对比图，(a)是一个实际地震道频谱，(b)是理论单频干扰波记录频谱，(c)是合成理论记录频谱，(d)是使用本发明去除单频干扰波后的地震道频谱，(e)是使用本发明检测出的单频干扰波频谱。从数据显示和频谱图中可以看出，本发明有效地消除了单频干扰波。

图3-4是本方法与商用VSP处理软件包处理VSP数据对比图。VSP数据是某地实际测量的三分量数据，每个接收点上记录X，Y和Z三个分量，排列上共有170个接收器，记录长度5000ms，采样间隔1ms，这样我们得到X，Y和Z三个分量各170道的记录。仅显示Z分量700～4000ms记录。图3是VSP数据Z分量处理对比图。(a)是原始VSP数据Z分量，(b)是VSP软件单频干扰波消除之后的VSP数据Z分量，(c)是本发明单频干扰波消除之后的VSP数据Z分量。对于Z分量，本发明和VSP软件单频干扰波消除方法处理效果之间的差异在图上可以明显看出来，VSP软件单频干扰波消除方法处理的结果上明显还存在剩余的单频干扰波，而本方法完全消除了单频干扰波。图4是VSP数据Z分量频谱对比图。(a)是原始分量的频谱，(b)是VSP软件单频干扰波消除方法处理后的频谱，(c)是本方法单频干扰波消除方法处理后的频谱。从频谱图上可以看出，VSP软件单频干扰波消除方法处理的结果上明显还存在剩余的频率为50HZ单频干扰波，而本发明完全消除了频率为50HZ的单频干扰波，同时有效地消除了频率为150HZ和250HZ的单频干扰波。

图5-6是本方法与时间域单频干扰波压制方法处理对比图。地震数据是某地3D叠前炮集数据，每炮1680道20个接收排列，记录长度6000ms，采样间隔2ms，仅显示了第一炮的1511-1680道2个接收排列800～6000ms，它们包含着强弱不等的50HZ附近的干扰波，图5是地震炮集数据处理对比图。(a)是地震数据，(b)是时间域单频干扰波压制方法消除单频干扰波之后的地震数据，(c)是本发明单频干扰波消除之后的地震数据。图6是不同地震道频谱对比图。(a)是原始分量的频谱，(b)是时间域单频干扰波压制方法消除单频干扰波之后的频谱，(c)是本发明单频干扰波消除方法处理后的频谱。从频谱图上可以看出，本方法完全消除了频率为50HZ的单频干扰波，而其它频率成分的有效波不但没有受到影响，而且在干扰波存在的频率成分上有效波也保存完好，说明本方法消除干扰波的方法并不伤害有效波，因此提高了单频附近信号的信噪比。很明显，本方法与时间域单频干扰波压制方法一样，有效地消除了单频干扰波，且两种算法看不出差异。在DELL工作站上，该数据第一炮的1680道，本发明的运算时间是34s，而同样的机器同样的参数，时间域单频干扰波压制方法需要146s，明显优于时间域单频干扰波压制方法。

Claims (3)

1、一种自适应识别和消除地震勘探单频干扰的方法，其特征在于采用如下步骤实现：

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

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

3)计算单频干扰波余弦函数和正弦函数的振幅A和B，其计算公式为：

$A=\frac{\mathrm{dc}-\mathrm{eb}}{\mathrm{ac}-b^2}$

$B=\frac{\mathrm{db}-\mathrm{ea}}{b^2-\mathrm{ac}}$

其中

$a=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\cos }^{2}2\mathrm{\πfi\Δt}$

$b=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\cos 2\mathrm{\πfi\Δ}t\sin 2\mathrm{\πfi\Δt}$

$c=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\sin }^{2}2\mathrm{\πfi\Δt}$

$d=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\cos 2\mathrm{\πfi\Δt}$

$e=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\sin 2\mathrm{\πfi\Δt}$

式中：

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

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

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

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

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

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

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

4)根据地震数据S_i和单频干扰波初始频率f₀，采用以下频率修正公式确定单频干扰波频率f，

f＝f₀+Δf

其中：

f₀是初始频率，Δf是自适应频率修正值，

$\mathrm{\Δf}=\frac{2(\mathrm{VA}-\mathrm{UB})}{N(A^2+B^2)}$

$U=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\cos 2\mathrm{\πfi\Δ}t\sin \frac{2\mathrm{\πi}}{N}$

$V=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{S}_i\sin 2\mathrm{\πfi\Δ}t\sin \frac{2\mathrm{\πi}}{N}$

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

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

y_i＝Acos 2πfiΔt+Bsin 2πfiΔt；

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

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

x_i＝s_i-y_i

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

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

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

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

Images