CN1015835B - 估算及抑制地震数据中异常干扰幅值的方法 - Google Patents

Abstract

本发明揭示了一种识别地震资料中异常于扰幅值的方法。对构成所述地震资料的预定数目的轨迹，在预定的时间期间内导出一个实测幅值因子，采用高斯-赛得尔矩阵转换法也可计算出对应于相同的时间期间及轨迹数目的无干扰辐值因子。计算出实测幅值因子与无干扰幅值因子之比；从而确定出一个较正因子，将其与一个预先定义的临界值相比较，来确定地震记录中异常干扰幅值的存在。用户可确定是否需要抑制或校正受到异常干扰幅值影响的轨迹。

Description

本发明一般说来涉及地震数据处理，特别是涉及估算及抑制地震数据中异常干扰幅值的方法。

过去数十年中，对地震记录资料的幅值问题做了广泛的研究。

Gutenberg（1936）参考Knott（1899）和Zoepprits等人（1912）早期所作的研究，讨论了来自地表下地震反射的期望幅值。在引入了自动增益控制（AGC）手段和记录设备之后，地球物理学更直接对时间异常现象和反射特性进行描述。但直到发展了二进制增益和浮点记录之后，反射幅值的研究才得到了普及。这些发展导致了“明点”解释技术，即将幅值异常的研究引入到地震解释中去。“明点”技术已被普遍接受。

这一技术首先被引用到，主要在三角洲区域的海洋地震数据的研究中，科学家们很快认识到，反射的亮度或暗度对地震解释具有重要的意义。Hermont（1969），O′Doherty和Anstey（1971），以及Lamer等人（1973）证明了反射幅值的重大价值。

将此项技术应用于陆地条件时，需要对影响记录幅值的因子进行更深入细致地研究。Sheriff（1973）、Taner和Koehler（1974）对这些因子作了考察，而Houba等人（1973）、Disher和Raudle（1973）讨论了一种陆地数据的处理方法。人们普遍认识到，有几个因子对接收到的地震信号的幅值有作用。

Larner等人（1973）、Larner（1975）、Taner和Koehler（1981）以及其他一些人确定了影响所接收的地震信号幅值的一些地表和地表下的特征。本发明的受让人早在1977年10月就提出了一种从地震数据中校正或消除异常干扰幅值的方法。Taner    Koehler（1981）提出了一个类似的计算实例，试图判别由于外界干扰。如近地表下层（风化层）的干扰所导致的异常幅值，此干扰可能会改变或影响所接收的地震信号的特性。

以上所有对地震数据处理中的异常幅值问题进行研究的所有人中，还没有一个人考虑到记算/校正过程中的任何干扰问题。

本发明的一个目的是要提供一种新的方法，在传统的数据处理之前进行质量控制并抑制地震数据中的异常干扰幅值信息。

本发明的另一个目的是要提供一种上述的方法，其中，由于幅值所相关的问题，在校正或抑制某一特定的地震轨迹（或轨迹的一部分）之前要主要考虑地震数据中的干扰。

设计这一计算机辅助方法，其目的是在传统的数据处理之前，估算、校正或抑制地震数据中存在的异常幅值。此方法给地球物理学家提供了在最后阶段处理之前对一预定轨迹估算干扰量的机会，这是因为在最后阶段要确定和消除干扰的影响非常困难。

通常对于每一冲击点的每一条信号分离地震数据轨迹要计算作为时间函数的幅值因子。时间函数是时限的序列，即计算所确定的和用户定义的时间期限。对每一轨迹计算与引起幅值变化的几个可变参数有关的数学模型，以确定一个实测幅值因子。对此数学模型的两边分别取对数，形成方程式从而计算出无干扰幅值因子，由计算出的无干扰幅值因子与实测幅值因子产生一个校正比率，并与从预定曲线和计算机打印的计算结果中选出的定义值相比较，如果校正比率大于一个给定的临界值，那么整条轨迹，或显示异常幅值的轨迹的那部分，在传统的数据处理之前要被校正或消除。

对本发明优点的更好的了解可通过以下根据附图作更详细的描述来实现。其中：

图1是传统的地震反射技术的略图;

图2是用于典型的地震探测线中的源接收器几何形状的一般示意图;

图3是与图2的几何形状相关的探测冲击中所作的冲击记录的一般示意图;

图4是图3中所示的冲击记录的放大部分，表示通常在冲击记录中所发现的异常幅值;

图5是经过所描述的方法校正后的图3所示冲击记录的实例。

地球物理探测的一个主要目的是对地球上某一特定区域的地表下结构提供一个精确的描述。所采用的手段之一即为地震反射剖面图。图1是传统的际地地震反射技术的一般示意图，其中，通过一振动或爆炸源在地球（P）表面的一个“冲击点”产生一个地震信号（S），此信号从一扩展的球形波前向下传播，受到具有可变声阻抗的地表下层的反射，反射的地震信号到达地表，被沿地球表面一条线或格上的诸点（R）所接收。被探测器（R）所接收的信号被转换为另一种信号，并通常被安置在一辆车（truck）上的一个中心记录设备所记录。对于某一给定的冲击点，一旦所有的数据都收集起来后，声信号发生器就被移到新的冲击点重新进行以上过程。

图2是用于上述地震探测中的一个源接收器的几何形状的一般示意图。图中包括五条标记1至5的接收器线，或称为“接收器分布”，表示沿地球表面的一条线A-B的探测线。每一接收器分布中还包括“X”表示的一个点，并分别标记为S₁、S₂、S₃等等，代表每一分布1至5中的冲击点。为解释方便起见，分布1至5彼此不重合，并且是通常记录在一典型的探测中的许多探测线的一个实例。实际上，每一地震探测分布1-5沿线A-B按一定距离彼此交错叠加在上述的探测分布上。

在一分布中的冲击点（源）到一特定的接收器的距离称之为偏移距离，每一分布1-5还有一个与探测中其它分布冲击共同的冲击点，这是在与沿一分布的总的接收器数目相关的同一区域所出现的冲击点，因此在每一分布中从冲击点到一特定接收器的距离（偏移）基本上是相同的，这处型式的探测本身还会导致地震数据的公共中点（CMP）或公共深度点（CDP）的积累。以上为本领域的普通技术人员所熟悉，因此其细节将不在此作进一步的讨论。

通常认为要获得对地表下的理想描述，探测和记录的信号除了反射信号的结构外，最好不要受到外界的影响。现在科学还没能提供出一种接收这种纯数据的方法。所接收到的数据受到许多的外界影响，因而数据的质量大受影响。下面所描述的本发明技术提供了一种能更好地区别影响地震数据质量的各干扰的作用的方法。

一般说来地震数据的质量与两大类因子有关。其一是用于探测的仪器，这包括用于产生信号的源，用于接收信号的探测器以及记录装置;第二个因子是地球地表下的地质结构，包括信号传播所经过的材料，沿地表风化层的深度、地表下层的密度，地表下层的连续性，等等。如上面的介绍中所简单提及的，地球科学家长期以来一直对地震信号的幅值大感兴趣。所接收的信号的幅值不仅是所接收的信号的强度的一般指示，而且代表了与所产生的信号相关的震相（位）。地震区域相邻两条轨迹的幅值的组合可帮助科学家确定地表下反射层位的特性。

在采用本发明的方法估算异常干扰幅值时，从具有完整地震资料的区域获取数据。图形信息提供 了沿地球表面探测的位置信息以及每个接收器或接收器组所处位置的信息。

图3一般地示出了与图2的几何结构相关的一个冲击点所作的冲击记录。记录的每条轨迹根据其沿探测线的位置在记录的顶部标上标记，而且其与冲击点的相对距离也能很容易地计算出来。垂向刻度为二向运行时间。对每条线所作的记录，如图3所示，包括1024不同的轨迹，甚至更多，每条轨迹代表由一单一探测器或探测器组在地球表面的一个已知区域所接收到的信号。如图所示，大多数的轨迹都有相对连续的幅值，一个值得注意的例外由标有AA的部分表示，它表示包含有带有一异常幅值的轨迹的记录的一部分。

图4为图3中AA部分所示的冲击记录的放大部分，由本图可比由图3更为清楚地区别轨迹与轨迹幅值的不同之处。要消除这些异常干扰幅值，我们遵循下面所描述的方法。

在此方法中，对冲击记录的每条轨迹都按图4所示t₁、t₂…t_n的时限序列进行检查。每一时限可能具有任意所需的时间长度并由用户所选择。在此我们假定时限的时间长度200毫秒。一旦用户选定了时限，在整个操作/估算过程中此时限都保持不变。但在对下一部分数据进行估算时可改变时限。

一旦确定了时限以及要估算的轨迹，下一步就是对每一轨迹计算作为时限函数的幅值。可采用RMS（均方根值）或峰幅值来代表幅值。当干扰持续较长时间（20到300毫秒）时最好采用RMS，此干扰在工业上称之为“斑点”（blobs）。一个斑点是指由以下因子所引起的干扰，这些因子是：擦过地震检波器的杂物、附近公路上的交通以及在记录过程中所出现的其它外来事件，图4标记N_b的幅值就是由上述干扰物引起的。当要识别由于记录中由一位或多位数字码弄污引起数字记录的错误从而导致的单一样本（尖峰信号）干扰时，则最好采用峰幅值，图4中标记N_s的幅值就属于此种情况。

假设源附近和接收器附近的幅值效应是地表连续的，而地表下的幅值变化是地表下连续的，从而对引起被观测幅值变化的参数建立起一个数学模型。更具体地说，考虑源i和接收器j的一条轨迹，在不同的水平时限h对其幅值作了测量：幅值因子Aijn可表示为冲击和接收器地表连续项、地表下连续地质项（CDP）、偏移项，以及通道数项的乘积。根据此数模，实测幅值Aijn可由下式表述;

Aijn＝Si＊Rj＊Gkn＊Oln＊Cc＊N（1）

其中：Aijn是实测幅值;

Si是与源位置i相关的幅值变化分量;

Rj是与探测器位置j相关的幅值变化分量;

Gkh是在时限h内与公共深度点位置（也叫中点）k相关的幅值变化分量;

Olh是在时限h内与偏移l相关的幅值变化分量;

Cc是与记录通道C相关的幅值变化分量;

N是干扰量。

对表达式（1）的等号两边取对数，得到下面的表达式（2）。

logAijn＝logSi+logRj+logGkh+logolh+logCc+logN    （2）

采用矩阵转换叠加法如高斯-赛得尔（Gauss-siedel）方法求解（2）式。在本方法中，表达式（2）中的项数可能很多。在进行地震数据处理时，要处理这样多数目的线性方程式，方程的解值必须精确而有效。可以采用多种求解方法，但我们最好是使用前面所述的高斯-赛得尔（Gauss-Seidel）代数矩阵求解法，下面我们参考一个线性拟合实例求解表达式（2.1）中a和b对这一求解方法的过程进行详细描述。

（2.1）

将（2.1）写成常规表达式形式，可得到

4 a+10 b＝13 （2.2）

以及10 a+30 b＝35.4 （2.3）

上述表达式又可改写成下述形式：

a＝3.25-2.500 b（2.4）

以及 b＝1.18-0.333 a（2.5）

因为高斯-赛得尔代数矩阵求解法是一种迭加方法，因此需要初值。为了开始进行迭加，我们使a＝b＝0，将b＝0代入（2.4）式，可得到a＝3.25，将这一解值代入（2.5）式可得到b＝0.0977。在进行第二轮迭加时，将b＝0.0977代入表达式（2.4）得到a＝3.0075，将a＝3.0075代入表达式（2.5）得到b＝0.1785，如此继续迭加，如表Ⅰ所示

对（2.1）式求解a和b的高斯-赛得尔迭加法列入表中的数据进行了舍入处理（表Ⅰ见后）

由迭加过程所得到的解值缓慢收敛于真解a＝1.8、b＝0.58收敛性并非总能保证，但如果初值选择合适，那么采用上述的高斯-赛得尔方法是可以收敛的。高斯-赛得尔方法的优点是它能快速求解大量相似的方程式。当将此方法用于方程式（2）时，转换过程中的一个基本假设是式（1）变量具有初值零。最好首先求出对地质或CDP位置有影响的平均幅值，利用与地质结构相关的这第一幅值。采用同样的过程得到偏移的幅值，然后求出通道的幅度值。除了确定干扰值以外，对所有轨迹的每一个变量都要重复以上过程。在进行了所有的计算之后，对与干扰量有关的幅值就不预考虑了。因为采用上述的高斯-赛得尔计算方法没有计算出干扰变量，因此，干扰变量的值在本方法的后面将进行论述。

对每一变量平均值的小误差几乎不必考虑，本方法的这一步骤只与每一信号相对于本记录中其它信号的平均幅值有关。采用矩阵转换时，一旦对每一变量的值确定后，将其以表达式（1）的相同方式重新组合，可产生一计算的幅值因子Bijh，这是一个无干扰幅值。每一轨迹时限都有两种幅值：一个实测幅值和一近似无干扰幅值。这两个值由表达式（3）以逐个时限（Window    by    window）的形式定标。

其中上标n是由用户定义的（或者为1或者为2）。由表达式（3）得到一个校正比率，指示是否要对轨迹作异常干扰幅值校正，或是否要在最后的数据处理步骤中从其后的相叠序列中将轨迹消除。校正值与用户定义的临界值（通常为1.0）相比较，以确定轨迹中是否存在着异常噪声干扰。可接受的数值范围在传统处理之前通过对曲线的分析来选择。

图4是由上述方法进行校正后的地震记录实例。要注意的是由线2开始的数据所标记的轨迹对幅值不再产生影响，并与其它的轨迹基本相似。

虽然本发明描述为计算机辅助计算、校正或抑制异常干扰幅值的方法，但熟悉本领域的技术人员完全可以用硬件来实施本发明，从而不必在数据处理中心也可以实施本方法。实施本发明的硬件，可在地震领域数据探测单元的多路转换器和记录仪的前面偶合。由海洋环境所采集的地震数据可同样采用相同的设备进行处理。

表Ⅰ

迭加次数    a    b

1    3.25    0.0977

5    2.4918    0.3502

10    2.0712    0.4902

15    1.9031    0.5462

20    1.8358    0.5686

25    1.8090    0.5779

29    1.7997    0.5807

Claims (7)

1、一种用于识别地震资料中异常干扰幅值数据的方法，包括如下步骤：

(a)对地震资料中的每一轨迹，计算出一个作为时间函数的实测幅值因子；

(b)对地震资料中一条轨迹的一部分，计算出一个作为时间函数的无干扰幅值因子；

(c)通过无干扰幅值因子和实测幅值因子计算出一个校正比率；

(d)将校正比率与一个预定值相比较，从而确定一条轨迹中是否存在一异常干扰幅值，以及在采用传统的地震数据处理方法之前是否需要校正异常值或消去轨迹；

其中对地震资料的每一轨迹计算作为时间函数的一个实测幅值因子时采用下面的表达式：

Aijh=Si＊Rj＊Gkh＊Olh＊Cc＊N

其中：Aijh是实测幅值；

Si是与源位置i相关的幅值变化分量；

Rj是与探测器位置j相关的幅值变化分量；

Gkh是在时限h内与公共深度点位置k(也叫中点)相关的幅值变化分量；

Olh是在时限h内与偏移1相关的幅值变化分量；

Cc是与记录通道c相关的幅值变化分量；

N是干扰量；

在计算无干扰幅值因子的步骤中进一步包括根据下面的表达式叠加计算矩阵转换：

logAijh=logSi+logRj+logGkh+logOlh+logCc+logN

其中计算矩阵转换的方法通过高斯-赛得尔(Gauss-Seidel)矩阵转换计算法完成，

其中计算校正比率的步骤通过下面的表达式完成：

其中上标n和临界值是由用户定义的。

2、一种对具有许多冲击轨迹的地震资料识别异常干扰幅值数据的方法，包括如下步骤：

（a）对预定数目的冲击轨迹，测量出一个作为时间函数的实测幅值因子;

（b）对预定数目的冲击轨迹，计算出一个作为时间函数的无干扰幅值因子;

（c）将无干扰幅值因子和实测幅值因子相比，导出一个校正比率;

（d）将校正比率与一个预定值临界值相比较，确定一条轨迹中是否存在一异常干扰幅值，

其中对一预定的冲击轨迹测量一个作为时间函数的幅值因子时采用下面的表达式：

Aijh＝Si＊Rj＊Gkh＊Olh＊Cc＊N

其中：

Aijh是实测幅值;

Si是与源位置i相关的幅值变化分量;

Rj是与探测器位置j相关的幅值变化分量;

Gkh是在时限h内与公共深度点（也叫中点）位置k相关的幅值变化分量;

Olh是在时限h内与偏移l相关的幅值变化分量;

Cc是与记录通道c相关的幅值变化分量;

N是干扰量;

其中计算无干扰幅值因子的步骤中进一步包括根据下面的表达式叠加计算矩阵转换：

logAijg＝logSl+logRj+logGkh+logOlh+logCc+logN

其中计算矩阵转换的方法通过高斯-赛得尔（Gauss-Seidel）矩阵转换计算法完成。

其中由无干扰幅值因子和实测幅值因子导出校正比率的步骤根据下面的表达式完成：

其中上标n和临界值是由用户定义的。

3、一种用于识别地震资料中异常干扰幅值数据的方法，对构成资料的预定数目的轨迹，在预定的时间期间内，导出一实测幅值因子，同时计算出一个无干扰幅值因子。将这两个因子相比，从而确定一个校正因子，将校正因子与一个预定的临界值相比较来确定在地震资料中异常干扰幅值的存在，并抑制或校正受到异常干扰幅值影响的轨迹;

其中对地震资料的每一轨迹计算作为时间函数的一个实测幅值因子时采用下面的表达式：

Aijh＝Si＊Rj＊Gkh＊Olh＊Cc＊N

其中：

Aijh是实测幅值;

Si是与源位置i相关的幅值变化分量;

Rj是与探测器位置j相关的幅值变化分量;

Gkh是在时限h内与公共深度点位置k（也叫中点）相关的幅值变化分量;

Olh是在时限h内与偏移l相关的幅值变化分量;

Cc是与记录通道c相关的幅值变化分量;

N是干扰量;

其中计算无干扰幅值因子的步骤中进一步包括根据下面的表达式叠加计算矩阵转换：

logAijh＝logSi+logRj+logGkh+logOlh

+logCc+logN

上述的计算矩阵转换的方法通过高斯-赛得尔（Gauss-Seidel）矩阵转换计算法完成，

其中计算校正比率的步骤通过下面的表达式完成：

其中上标n和临界值是由用户定义的。

4、如权利要求1所述的方法，进一步包括在进行传统的地震数据处理之前，根据计算的校正比率校正受到异常干扰幅值影响的轨迹的步骤。

5、如权利要求1所述的方法，进一步包括在进行传统的地震数据处理之前，根据所计算的校正比率，抑制受到异常干扰幅值影响的轨迹的步骤。

6、如权利要求2所述的方法，进一步包括在进行传统的地震数据处理之前，根据计算的校正比率与预定的临界值的比较，对受到异常干扰幅值影响的轨迹进行校正步骤。

7、如权利要求2所述的方法，进一步包括在进行传统的地震数据处理之前，根据计算的校正比率与预定的临界值的比较，抑制受到异常干扰幅值影响的轨迹。

Images