CN102736108A - 基于样条拟合的真三维地震数据噪声压制方法 - Google Patents

Abstract

一种基于样条拟合的真三维地震数据噪声压制方法，所述方法包括：(1)对三维叠加数据体以窗口处理方式进行相位时间多项式的拟合处理；(2)利用正交多项式拟合纵横测线方向上各道的均方根振幅；(3)把窗口内纵横测线方向上的数据沿同相轴走向进行叠加，叠加后的结果再归一化；(4)将期望波形乘以各道拟合的均方根振幅；(5)把窗口按纵横测线方向和时间方向同时移动半个窗口步长，并重复步骤(1)至步骤(4)；(6)利用能量加权法对窗口的重叠部分的地震数据进行处理；(7)确定是否完成似合，如果没有完成拟合，则重复步骤(5)至步骤(6)的操作，如果完成拟合，则进行步骤(8)；(8)对拟合后的叠加数据进行混波处理。

Description

基于样条拟合的真三维地震数据噪声压制方法

技术领域

本发明涉及石油地震勘探，属于地震勘探资料处理与解释领域，具体地说，涉及一种针对复杂山地低信噪比资料中的噪声进行压制的基于样条拟合的真三维地震数据噪声压制方法。

背景技术

复杂山地地震资料噪声来源十分复杂，噪声对有效波有较强的干扰，尤其是无处不在的随机干扰，导致地震资料的信噪比低，严重影响了地震数据的分析、处理。而在地震资料处理中，信噪比和分辨率问题又一直是关注的焦点，这两者具有相辅相成、互为依赖的特点，提高地震资料的信噪比是提高地震勘探数据处理质量的一个关键步骤，而提高分辨率在很大程度上又受到地震数据信噪比的限制，所以人们追求一种在不损害信号分辨率的前提下提高资料的信噪比的新技术。

随着地震勘探技术的发展，三维地震资料采集在地震勘探中所占的比例逐渐增大。但是，目前对于三维数据的拟合方法大多是在二维上进行处理的，对于三维叠后数据体来说，二维拟合的方法只能在纵测线方向或者是横测线方向单独进行处理，而并不能对地震资料数据实现整体的拟合。这样往往会产生一些问题：当在地震资料的一个方向上进行拟合时，可以使该方向的同相轴更加平滑，但是可能会引起垂直于该方向的叠加剖面产生突跳现象，也就是某些地震道的能量会突然上移或者下移，形成一个垂直的断点，这样就影响了剖面的效果。

因此，基于上面提到的二维拟合产生的问题，以及三维地震资料在地震勘探中的普遍使用，对三维地震资料进行真三维处理，是有效利用三维资料采集的信息来获得更好地质成果的根本途径，同时也是三维采集对数据处理的要求。

发明内容

本发明的一方面在于提供了针对三维地震进行真三维处理的样条拟合技术，该技术不但能够提高同相轴的连续性，而且有效地解决了传统二维拟合后叠加剖面可能产生突跳现象的缺陷，提高了叠加剖面的质量。

本发明的另一方面在于提供了一种新的三维数据窗口的扫描方法，从而有效地确定三维数据体中同相轴的拟合多项式的系数。

本发明的另一方面在于提供了一种基于能量加权法的技术来处理窗口重叠部队发的三维数据的方法。

根据本发明的多方面，提供了一种基于样条拟合的真三维地震数据噪声压制方法，所述方法包括：(1)对三维叠加数据体以窗口处理方式进行相位时间多项式的拟合处理；(2)在窗口内，利用正交多项式拟合纵横测线方向上各道的均方根振幅；(3)把窗口内纵横测线方向上的数据沿同相轴走向进行叠加，叠加后的结果再归一化得到拟合信号的期望波形；(4)将期望波形乘以各道拟合的均方根振幅，并放到计算出的相位时间位置上，从而得到了一个窗口内模拟的地震剖面；(5)把窗口按纵横测线方向和时间方向同时移动半个窗口步长，并重复步骤(1)至步骤(4)；(6)利用能量加权法对窗口的重叠部分的地震数据进行处理；(7)确定是否完成似合，如果没有完成拟合，则重复步骤(5)至步骤(6)的操作，如果完成拟合，则进行步骤(8)；(8)对拟合后的叠加数据进行混波处理。

根据本发明的基于样条拟合的真三维地震数据噪声压制方法能够有效地压制实际地震资料中的噪声，并使有效信号更加突出，很好地保持实际地震记录所反映的地层信息和地质特征，提供了信噪比高的叠加剖面。根据本发明的基于样条拟合的真三维地震数据噪声压制方法不但能够提高同相轴的连续性，还有效地解决了传统二维拟合后叠加剖面产生突跳现象的缺陷，提高了叠加剖面的质量。尤其是对山地资料中信噪比较低的地震数据，根据本发明的基于样条拟合的真三维地震数据噪声压制方法的去噪效果较好，保真度更高。

附图说明

从以下结合附图对本发明的示例性实施例进行的描述中，本发明的这些和/或其它方面及优点将会变得清楚，并且更易于理解，其中：

图1是示出根据本发明的基于样条拟合的真三维地震数据噪声压制方法的流程图；

图2是在图1所示的方法中多道互相关系数扫描范围选取示图；

图3A是根据实施例的某地区三维叠后原始资料中横测线方向(CrossLine线)上的叠加剖面；

图3B是根据实施例的在某地区三维叠后原始资料的纵测线方向(InLine线)上做了常规二维拟合后横测线方向上的叠加剖面；

图3C是根据实施例的对某地区三维叠后原始资料做了根据本发明的真三维拟合后横测线方向上的叠加剖面；

图4A是根据另一实施例的某地区三维叠后原始资料中纵测线方向上的叠加剖面；

图4B是根据另一实施例的对某地区三维叠后原始资料做了根据本发明的真三维拟合后纵测线方向上的叠加剖面；

图4C是根据另一实施例的在某地区三维叠后原始资料中横测线方向上的叠加剖面；

图4D是根据另一实施例的对某地区三维叠后原始资料做了根据本发明的真三维拟合后横测线方向上的叠加剖面。

具体实施方式

以下，参照附图来详细描述示例性实施例以使本领域的普通技术人员更易于理解。目前声明于此的发明的示例性实施例实际上可包含各种形式，并不局限于显示和描述于此的示例。为了清晰，对公知的结构和功能的描述可被省略，并且在整个描述中，相同的标号表示相同的元件。

多项式拟合技术在提高信噪比的同时基本上不降低地震资料的分辨率，即信号的高频成分尽可能不受损失，也能保持各道的相对振幅，所以此方法在二维地震资料的处理上得到广泛的应用。

本发明提出一种基于样条拟合的真三维地震数据噪声压制方法，该方法基于地震道数据具有横向相干性的原理，用波形相似性准则假设地震记录相位时间横向变化可用一个二次多项式表示，沿着相位时间变化的各道振幅变化也可用一个待定系数多项式表示，通过多项式拟合求出地震信号相位时间，标准波形和振幅加权系数，然后将它们组合成拟合地震道，拟合地震道与原始道混波后得到最终信噪比、分辨率及保真度较高的地震数据。

图1是示出根据本发明的基于样条拟合的真三维地震数据噪声压制方法的流程图。

参照图1，在步骤101，对三维叠加数据体以窗口处理方式进行相位时间多项式的拟合处理。

在步骤102，在窗口内，利用正交多项式拟合纵横测线方向上各道的均方根振幅。

在步骤103，把窗口内纵横测线方向上的数据沿同相轴走向进行叠加，叠加后的结果再归一化得到拟合信号的期望波形。

在步骤104，将期望波形乘以各道拟合的均方根振幅，并放到计算出的相位时间位置上，从而得到了一个窗口内模拟的地震剖面。

在步骤105，把窗口按纵横测线方向和时间方向同时移动半个窗口步长，并重复步骤101至104。

在步骤106，利用能量加权法对窗口的重叠部分的地震数据进行处理。

在步骤107，确定是否完成似合。如果没有完成拟合，则重复步骤105至步骤106的操作，如果完成拟合，则可得到噪声压制后的三维拟合叠后数据，并进入步骤108。

在步骤108，对拟合后的叠加数据进行混波处理，最终得到信噪比、分辨率及保真度较高的地震数据。

以下将详细描述根据本发明的基于样条拟合的真三维地震数据噪声压制方法的各步骤。

现在对图1的步骤101(即，对三维叠加数据体以窗口处理方式进行相位时间多项式的拟合处理)进行详细描述。

在三维叠后地震资料上任取一个区域，对该区域按照时间轴划分一系列的窗口。设某一窗口W内横测线方向有2N+1个地震道，纵测线方向有2M+1个地震道，对该窗口在xy方向即纵横测线方向建立离散坐标系如下等式(1)：

D＝{(x，y)|x∈[-M，M]，y∈[-N，N]；x，y∈Z}.                (1)

在窗口W内，可用二元三次多项式来表示地震信号的到达时间(或窗口的中点时间)，如下等式(2)：

T(x，y)＝a₀₀+a₁₀x+a₁₁y+a₂₀x²+a₂₁xy+a₂₂y²+a₃₀x³+a₃₁x²y+a₃₂xy²+a₃₃y³    (2)

其中，x、y为某一道位置的相对道序号，a_ij为时间多项式系数，i为x项的乘方次数与y项的乘方次数之和，j为y项的乘方次数。

设窗口W时间长度为2L，采样时间间隔为Δt，则相关样点所在的空间为如下等式(3)：

$H=\left\{(x,y,t)\right|x\∈[-M,M],y\∈[-N,N],t\∈[T(x,y)-L,T(x,y)+L];x,y,\frac{t}{\mathrm{\Δt}}\∈Z\}.---\left(3\right)$

其中，xy为纵横测线方向，t为时间方向。

分别取[-M，M]与[-N，N]上的二次正交多项式，如以下等式(4)：

p₀(x)＝1，p₁(x)＝x， $p_2\left(x\right)=x^2-\frac{1}{3}M(M+1),$ $p_3\left(x\right)=x^3+\frac{1}{5}({3M}^2+3M-1)x,$ x∈[-M，M]

q₀(y)＝1，q₁(y)＝y， $q_2\left(y\right)=y^2-\frac{1}{3}N(N+1),$ $q_3\left(y\right)=y^3+\frac{1}{5}({3N}^2+3N-1)y,$ y∈[-N，N](4)

容易证明对于任意0≤i，j≤2，有如下等式(5)：

$\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}{p}_i\left(x\right)p_j\left(x\right)=\left\{\begin{array}{cc}0,& i\≠j\\ >0,& i=j\end{array}\right.,$

$\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{q}_i\left(x\right)q_j\left(x\right)=\left\{\begin{array}{cc}0,& i\≠j\\ >0,& i=j\end{array}\right..---\left(5\right)$

令F＝{F_ij|F_ij(x，y)＝p_i(x)q_j(y)，0≤i，j≤2}，则F为D上的正交多项式集，令G＝{G_i|0≤i≤8}＝F，易知G为D上最高次为四次的正交多项式集，将G中元素如下等式(6)：

G₀＝1，G₁＝x，G₂＝y，G₃＝xy，

$G_4=x^2-\frac{1}{3}M(M+1),$ $G_5=y^2-\frac{1}{3}N(N+1),$

$G_6=x^{2}y-\frac{1}{3}M(M+1)y,$ $G_7={\mathrm{xy}}^2-\frac{1}{3}N(N+1)x,$

$G_8=x^2{y}^2-\frac{1}{3}N(N+1)x^2-\frac{1}{3}M(M+1)y^2+\frac{1}{9}\mathrm{MN}(M+1)(N+1).---\left(6\right)$

则窗口W上的窗口中点时间拟合多项式如下等式(7)：

$T(x,y)=\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{c}_i{G}_i=c_0+\underset{i=1}{\overset{8}{\mathrm{\Σ}}}{c}_i{G}_i,---\left(7\right)$

其中，c_i是拟合系数。在等式(7)中，由于多项式集G中的G₄，G₅，G₈含有常数项，所以在拟合过程中拟合中心随多项式系数变化而变化，当地震数据的时空变化很剧烈的时候误差将会比较大。因此可对原算法进行改进：

取正交多项式集G＝{G₀，G₁，G₂，G₃，G₆，G₇}，则时间拟合多项式(7)变为如下等式(8)：

$T^{\′}(x,y)=\underset{i=0}{\overset{4}{\mathrm{\Σ}}}{c}_i{G}_i^{\′}=c_0+\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{c}_i{G}_i^{\′}---\left(8\right)$

由于G₁，G₂，G₃，G₆，G₇均不含常数项，故c₀作为拟合中心是不变的，可以有效地降低G作为拟合多项式拟合中心不定而产生的误差。

根据多道互相关性最强的原则确定最佳拟合系数。对窗口内地震数据S(x，y，t)可以计算归一化多道互相关系数，如下等式(9)：

$R(c_0,...,c_8)=\frac{\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}\{{\left[\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}S(x,y,T(x,y)+t)\right]}^2-\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{S}^2(x,y,T(x,y)+t)\}}{[(2M+1)\×(2N+1)-1]\×\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{S}^2(x,y,T(x,y)+t)},---\left(\text{9}\right)$

如果沿同相轴方向对同一时间层的各道数据进行相关性计算，则等式(9)变为如下等式(10)：

$R(c_0,...,c_8)=\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}\frac{{\left[\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}S(x,y,T(x,y)+t)\right]}^2-\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{S}^2(x,y,T(x,y)+t)}{[(2M+1)\×(2N+1)-1]\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{S}^2(x,y,T(x,y)+t)},---\left(10\right)$

如果任意两道按时间方向进行相关性计算，然后对所有可能的道组合累加，则等式(9)变为如下等式(11)：

$R(c_0,...,c_8)=\underset{x,x^{\′}=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y,y^{\′}=-N}{\overset{N}{\mathrm{\Σ}}}\frac{\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}s(x,y,T(x,y)+t)\×s(x^{\′},y^{\′},T(x^{\′},y^{\′})+t)}{{[\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}{s}^2(x,y,T(x,y)+t)\×s^2(x^{\′},y^{\′},T(x^{\′},y^{\′})+t)]}^{1/2}},---\left(11\right)$

就计算量而言，等式(9)或(10)的计算效率较高，而等式(11)的计算量较大。

系数c₀，...，c₈扫描的具体实现过程是：固定c₀，对c₁，c₂，......，c₈扫描。先扫描c₁，此时c₂，......，c_k＝0，比较多道互相关值确定c₁；保持已确定的多道互相关系数c₁，再扫描c₂；依此类推。窗口的形状随多项式系数的变化而变化，求系数的过程就是搜索信号的过程，由扫描确定的窗口就是期望信号的窗口。在扫描过程中，可以仅选择空间窗口中任意一对时间方向的对角边进行扫描(如图2所示)，这样，在保证有相同拟合效果的同时能够减少运算量。

按照一定步长(即，半个窗口步长)滑动窗口，可求出各个窗口内的相位时间拟合多项式。

现在对图1的步骤102(即，利用正交多项式拟合纵横测线方向上各道的均方根振幅)进行详细描述。

当时间拟合多项式系数确定c₀，...，c₈以后，便可对其所确定的窗口内地震信号的振幅来进行多项式拟合，可用下列多项式(12)来表示振幅拟合多项式：

A(x，y)＝b₀₀+b₁₀x+b₁₁y+b₂₀x²+b₂₁xy+b₂₂y²+b₃₀x³+b₃₁x²y+b₃₂xy²+b₃₃y³      (12)

在窗口内，各点的输出振幅通过对原始道振幅分别沿纵横测线的同相轴方向按距离加权叠加得到，具体公式如下等式(13)：

$A_{\mathrm{xy}}=\underset{i=-M}{\overset{M}{\mathrm{\Σ}}}\underset{j=-N}{\overset{N}{\mathrm{\Σ}}}{W}_{\mathrm{xyij}}{B}_{\mathrm{ij}},---\left(13\right)$

其中，A_xy表示坐标为(x，y)点的地震道的输出振幅，W_xyij为点(i，j)对点(x，y)的加权系数，一般为两点距离d((x，y)，(i，j))的函数，比如倒数函数；B_ij表示点(i，j)的原始振幅。

当得到A_xy之后，便可根据最小二乘法拟合出振幅多项式(12)的系数b₀₀，b₁₀，b₁₁，......，当多项式次数较高时，为了避免法方程系病态，也可将振幅多项式(12)改写成如下等式(14)：

A(x，y)＝d₀₀p₀(x)q₀(y)+d₁₀p₁(x)+d₁₁q₁(y)+d₂₀p₂(x)+d₂₁p₁(x)q₁(y)+d₂₂q₂(y)+......(14)

其中，p₀，q₀，p₁，q₁......为等式(4)所描述的正交多项式，采用最小二乘法确定系数d₀₀，d₁₀，d₁₁......。

现在对图1的步骤103(即，把窗口内纵横测线方向上的数据沿同相轴走向进行叠加，叠加后的结果再归一化得到拟合信号的期望波形)进行详细描述。

在窗口内，沿纵横测线上同相轴方向对信号振幅值进行累加，然后对累加后得到的单道的振幅值进行均方根平均，就能得到期望波形。这时，由于进行累加和均方根平均，期望波形中各道的振幅值比较均匀，各道的振幅值不是实际真实的振幅值。利用最小二乘方式拟合上述各道的振幅值，将该道拟合均方根振幅乘以标准波形，即可得到该道的实际波形。

现在对图1的步骤106(即，利用能量加权法对窗口的重叠部分的地震数据进行处理)进行详细描述。

窗口在时间方向上移动过程中，相邻窗口之间有半个窗口的重叠，为了不使剖面形态产生突变，对不同窗口重叠部分采用斜坡加权平均的方法来处理这个问题，即：对于窗口沿时间方向移动的情况，把窗口分成上下两个窗口，将上窗口中拟合的波形乘以权系数，其中，上窗口的权系数从上到下由1按样点数线性地减小到0，将下窗口中拟合的波形也乘以与上窗口的权系数的大小的方向相反的权系数，这样对于重叠部分任何采样点上，上窗口的权系数与下窗口的权系数之和恒为1，最后将得到的两个新波形相加即得到重叠部分波形。

对于窗口沿纵横测线方向移动的情况，也用同样的方法处理，把窗口分成左右两个窗口，加权方式采用能量加权平均方式，其某道的权系数由左右两个窗口在该道的能量值决定，左窗口某重叠道权系数为左窗口在该道的能量比前一个右窗口在该道的能量与左窗口在该道的能量之和；右窗口某重叠道权系数为右窗口在该道的能量比后一个左窗口在该道的能量与右窗口在该道的能量之和。这样对于重叠部分任何采样点上，左窗口的权系数与右窗口的权系数之和恒为1，最后将得到的两个新波形相加即得到重叠部分波形。

例如，窗口A、B在同一道的能量为Ea、Eb，则该道能量为 $\left(\frac{\mathrm{Ea}}{\mathrm{Ea}+\mathrm{Eb}}\right)\·\mathrm{Ea}+\left(\frac{\mathrm{Eb}}{\mathrm{Ea}+\mathrm{Eb}}\right)\·\mathrm{Eb}.$ 其中，

为左窗口权系数，

为右窗口权系数， $\frac{\mathrm{Ea}}{\mathrm{Ea}+\mathrm{Eb}}+\frac{\mathrm{Eb}}{\mathrm{Ea}+\mathrm{Eb}}\≡1.$

现在对图1的步骤108(即，对拟合后的叠加数据进行混波处理，最终得到信噪比、分辨率及保真度较高的地震数据)进行详细描述。

将拟合出来的叠加数据与原始输入数据按以下的等式(15)进行混波处理：

S′(x，y，t)＝S(x，y，t)×P+T(x，y，t)(1-P)                    (15)

其中，S′(x，y，t)为输出记录；S(x，y，t)为原始输入记录；T(x，y，t)为拟合记录；P为混波比。

以下，将详细描述根据本发明的基于样条拟合的真三维地震数据噪声压制方法应用到实际地震资料的效果。

图3A是根据实施例的某地区三维叠后原始资料中横测线方向(CrossLine线)上的叠加剖面；图3B是根据实施例的在某地区三维叠后原始资料的纵测线方向(InLine线)上做了常规二维拟合后横测线方向上的叠加剖面；图3C是根据实施例的对某地区三维叠后原始资料做了根据本发明的真三维拟合后横测线方向上的叠加剖面。以下将简单说明在实际应用中真三维拟合方法。

在三维叠后地震数据体上任取一个区域，对该区域按时间轴划分一系列窗口。如假设某一窗口W内横测线方向上2N+1个地震道，纵测线方向上有2M+1个地震道，对该窗口在xy方向沿纵横向方向上，利用正交多项式确定这个窗口内的相位时间同相轴，如等式(8)：

$T^{\′}(x,y)=\underset{i=0}{\overset{4}{\mathrm{\Σ}}}{c}_i{G}_i^{\′}=c_0+\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{c}_i{G}_i^{\′}---\left(8\right)$

根据多道互相关性最强的原则确定最佳拟合系数(即，等式(8)的系数c_i，对窗口内地震数据S(x，y，t)可以计算归一化多道互相关系数，如等式(10)：

$R(c_0,...,c_8)=\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}\frac{{\left[\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}S(x,y,T(x,y)+t)\right]}^2-\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{S}^2(x,y,T(x,y)+t)}{[(2M+1)\×(2N+1)-1]\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{S}^2(x,y,T(x,y)+t)},---\left(10\right)$

系数扫描的具体过程实现过程是：固定c₀，对c₁，c₂，......，c₈扫描。先扫描c₁，此时c₂，......，c_k＝0，比较多道互相关值确定c₁；保持已确定的多道互相关系数，再扫描c₂；依此类推。窗口的形状随多项式系数的变化而变化，求系数的过程就是搜索信号的过程，由扫描确定的窗口就是期望信号的窗口。在扫描过程中，可以仅选择空间窗口中任意一对时间方向的对角边进行扫描(如图2所示)，这样，在保证有相同拟合效果的同时能够减少运算量。

当窗口W中的时间多项式确定以后，利用正交多项式拟合纵横测线方向上各道的均方根振幅，即在窗口W内，各点的输出振幅通过对原始道振幅分别沿纵横测线的同相轴方向按距离加权叠加得到，具体等式(13)如下：

$A_{\mathrm{xy}}=\underset{i=-M}{\overset{M}{\mathrm{\Σ}}}\underset{j=-N}{\overset{N}{\mathrm{\Σ}}}{W}_{\mathrm{xyij}}{B}_{\mathrm{ij}},---\left(13\right)$

当振幅多项式确定之后，根据以下步骤得到拟合信号的期望波形：将时间多项式确定的同一窗口W内横测线方向上2N+1的道记录波形和纵测线方向上2M+1的道记录波形同时沿拟合出的同相轴方向相加，并对相加结果缩放，使其均方根振幅归一化得到拟合信号的期望波形。

对期望波形乘上各道拟合的均方根振幅并放到计算出的相位时间位置上，就得到了一个窗口内模拟的地震剖面。

把窗口按纵横测线方向和时间方向同时移动半个窗口步长，重复上述步骤的操作，就可得到最终拟合后的叠加数据。

这里需注意的是窗口在时间方向上移动过程中，相邻窗口之间有半个窗口的重叠，为了不使剖面形态产生突变，对不同窗口重叠部分采用斜坡加权平均的方法来处理这个问题，即：对于窗口沿时间方向移动的情况，把窗口分成上下两个窗口，将上窗口中拟合的波形乘以权系数，其中，上窗口的权系数从上到下由1按样点数线性地减小到0，将下窗口中拟合的波形也乘以与上窗口的权系数的大小的方向相反的权系数，这样对于重叠部分任何采样点上，上窗口的权系数与下窗口的权系数之和恒为1，最后将得到的两个新波形相加即得到重叠部分波形。

对于窗口沿纵横测线方向移动的情况，也用同样的方法处理，把窗口分成左右两个窗口，加权方式采用能量加权平均方式，其某道的权系数由左右两个窗口在该道的能量值决定，左窗口某重叠道权系数为左窗口在该道的能量比前一个右窗口在该道的能量与左窗口在该道的能量之和；右窗口某重叠道权系数为右窗口在该道的能量比后一个左窗口在该道的能量与右窗口在该道的能量之和。这样对于重叠部分任何采样点上，左窗口的权系数与右窗口的权系数之和恒为1，最后将得到的两个新波形相加即得到重叠部分波形。

对拟合后的地震数据进行如下等式(15)的混波处理，将拟合出来的叠加数据与原始输数据按下式进行混波处理：

S′(x，y，t)＝S(x，y，t)×P+T(x，y，t)(1-P)                (15)

其中，S′(x，y，t)为输出记录；S(x，y，t)为原始输入记录；T(x，y，t)为拟合记录；P为混波比。

通过上面的步骤最终完成基于多项式拟合的真三维地震数据噪声压制处理。其压制效果如图3A至3C所示，与图3A相比，图3B和图3C中同相轴的连续性都有比较明显的改善，信噪比得到明显的提高。但是图3B的同相轴有明显的突跳出现，这种突跳现象破坏了整张剖面的效果，降低了剖面的质量，这也是常规二维拟合方法经常出现的问题。图3C却没有这种突跳现象，这说明了根据本发明的真三维拟合方法能够有效地提高地震资料剖面的质量，拟合的地震资料为后续的精细解释提供信噪比更高、保真度更好的地震剖面。

图4A至图4D相应于地质情况更为复杂的实际地震数据。高陡复杂地区信噪比较低，尤其浅层有效信号被噪声淹没，虽然本身的复杂性导致噪声处理难度较大，但是我们可通过一个合适的拟合窗口，对地震数据进行最优拟合去噪处理。

图4A是根据另一实施例的某地区三维叠后原始资料中纵测线方向上的叠加剖面；图4B是根据另一实施例的对某地区三维叠后原始资料做了根据本发明的真三维拟合后纵测线方向上的叠加剖面；图4C是根据另一实施例的在某地区三维叠后原始资料中横测线方向上的叠加剖面；图4D是根据另一实施例的对某地区三维叠后原始资料做了根据本发明的真三维拟合后横测线方向上的叠加剖面。

从图4A至图4D多项式拟合去噪前后对比分析知道：真三维拟合噪声压制技术能使得噪声得到更好的压制，尤其是浅层的有效信号更加突出，同相轴的连续性也更好，能更好地保持实际地震记录所反映的地层信息和地质特征，为后续的处理、解释提供信噪比更高的剖面。

虽然已参照实施例显示和描述本发明，但是本领域的技术人员应该理解，在不脱离由权利要求及其等同物限定的本发明的精神和范围的情况下，可以对其形式和细节进行各种改变。

Claims (8)

1.一种基于样条拟合的真三维地震数据噪声压制方法，所述方法包括：

(1)对三维叠加数据体以窗口处理方式进行相位时间多项式的拟合处理；

(2)在窗口内，利用正交多项式拟合纵横测线方向上各道的均方根振幅；

(3)把窗口内纵横测线方向上的数据沿同相轴走向进行叠加，叠加后的结果再归一化得到拟合信号的期望波形；

(4)将期望波形乘以各道拟合的均方根振幅，并放到计算出的相位时间位置上，从而得到了一个窗口内模拟的地震剖面；

(5)把窗口按纵横测线方向和时间方向同时移动半个窗口步长，并重复步骤(1)至步骤(4)；

(6)利用能量加权法对窗口的重叠部分的地震数据进行处理；

(7)确定是否完成似合，如果没有完成拟合，则重复步骤(5)至步骤(6)的操作，如果完成拟合，则进行步骤(8)；

(8)对拟合后的叠加数据进行混波处理。

2.如权利要求1所述的方法，其中，步骤(1)包括以下步骤：

在三维叠后地震资料上任取一个区域，对该区域按照时间轴划分一系列的窗口，设某一窗口W内横测线方向有2N+1个地震道，纵测线方向有2M+1个地震道，对该窗口在xy方向即纵横测线方向建立离散坐标系如下等式(1)：

D＝{(x，y)|x∈[-M，M]，y∈[-N，N]；x，y∈Z}.                (1)

在窗口W内，用二元三次多项式来表示地震信号的到达时间如下等式(2)：

T(x，y)＝a₀₀+a₁₀x+a₁₁y+a₂₀x²+a₂₁xy+a₂₂y²+a₃₀x³+a₃₁x²y+a₃₂xy²+a₃₃y³      (2)

其中，x、y为某一道位置的相对道序号，a_ij为时间多项式系数，i为x项的乘方次数与y项的乘方次数之和，j为y项的乘方次数，

设窗口W时间长度为2L，采样时间间隔为Δt，则相关样点所在的空间为如下等式(3)：

$H=\left\{(x,y,t)\right|x\∈[-M,M],y\∈[-N,N],t\∈[T(x,y)-L,T(x,y)+L];x,y,\frac{t}{\mathrm{\Δt}}\∈Z\}.---\left(3\right)$

其中，xy为纵横测线方向，t为时间方向，

分别取[-M，M]与[-N，N]上的二次正交多项式，如以下等式(4)：

p₀(x)＝1，p₁(x)＝x， $p_2\left(x\right)=x^2-\frac{1}{3}M(M+1),$ $p_3\left(x\right)=x^3+\frac{1}{5}({3M}^2+3M-1)x,$ x∈[-M，M]

q₀(y)＝1，q₁(y)＝y， $q_2\left(y\right)=y^2-\frac{1}{3}N(N+1),$ $q_3\left(y\right)=y^3+\frac{1}{5}({3N}^2+3N-1)y,$ y∈[-N，N](4)

容易证明对于任意0≤i，j≤2，有如下等式(5)：

$\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}{p}_i\left(x\right)p_j\left(x\right)=\left\{\begin{array}{cc}0,& i\≠j\\ >0,& i=j\end{array}\right.,$

$\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{q}_i\left(x\right)q_j\left(x\right)=\left\{\begin{array}{cc}0,& i\≠j\\ >0,& i=j\end{array}\right..---\left(5\right)$

令F＝{F_ij|F_ij(x，y)＝p_i(x)q_j(y)，0≤i，j≤2}，则F为D上的正交多项式集，令G＝{G_i|0≤i≤8}＝F，G为D上最高次为四次的正交多项式集，将G中元素如下等式(6)：

G₀＝1，G₁＝x，G₂＝y，G₃＝xy，

$G_4=x^2-\frac{1}{3}M(M+1),$ $G_5=y^2-\frac{1}{3}N(N+1),$

$G_6=x^{2}y-\frac{1}{3}M(M+1)y,$ $G_7={\mathrm{xy}}^2-\frac{1}{3}N(N+1)x,$

$G_8=x^2{y}^2-\frac{1}{3}N(N+1)x^2-\frac{1}{3}M(M+1)y^2+\frac{1}{9}\mathrm{MN}(M+1)(N+1).---\left(6\right)$

则窗口W上的窗口中点时间拟合多项式如下等式(7)：

$T(x,y)=\underset{i=0}{\overset{8}{\mathrm{\Σ}}}{c}_i{G}_i=c_0+\underset{i=1}{\overset{8}{\mathrm{\Σ}}}{c}_i{G}_i,---\left(7\right)$

其中，c_i是拟合系数，取正交多项式集G′＝{G₀，G₁，G₂，G₃，G₆，G₇}，则将时间拟合多项式(7)改变为如下等式(8)：

$T^{\′}(x,y)=\underset{i=0}{\overset{4}{\mathrm{\Σ}}}{c}_i{G}_i^{\′}=c_0+\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{c}_i{G}_i^{\′}---\left(8\right)$

根据多道互相关性最强的原则确定最佳拟合系数，对窗口内地震数据S(x，y，t)计算归一化多道互相关系数，如下等式(9)：

$R(c_0,...,c_8)=\frac{\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}\{{\left[\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}S(x,y,T(x,y)+t)\right]}^2-\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{S}^2(x,y,T(x,y)+t)\}}{[(2M+1)\×(2N+1)-1]\×\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{S}^2(x,y,T(x,y)+t)},---\left(\text{9}\right)$

系数c₀，...，c₈扫描的具体实现过程是：固定c₀，对c₁，c₂，......，c₈扫描，先扫描c₁，此时c₂，......，c_k＝0，比较多道互相关值确定c₁，保持已确定的多道互相关系数c₁，再扫描c₂，依此类推。

3.如权利要求2所述的方法，其中，步骤(1)还包括：

如果沿同相轴方向对同一时间层的各道数据进行相关性计算，则等式(9)为如下等式(10)：

$R(c_0,...,c_8)=\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}\frac{{\left[\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}S(x,y,T(x,y)+t)\right]}^2-\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{S}^2(x,y,T(x,y)+t)}{[(2M+1)\×(2N+1)-1]\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-N}{\overset{N}{\mathrm{\Σ}}}{S}^2(x,y,T(x,y)+t)},---\left(10\right)$

如果任意两道按时间方向进行相关性计算，然后对所有可能的道组合累加，则等式(9)为如下等式(11)：

$R(c_0,...,c_8)=\underset{x,x^{\′}=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y,y^{\′}=-N}{\overset{N}{\mathrm{\Σ}}}\frac{\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}s(x,y,T(x,y)+t)\×s(x^{\′},y^{\′},T(x^{\′},y^{\′})+t)}{{[\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}{s}^2(x,y,T(x,y)+t)\×s^2(x^{\′},y^{\′},T(x^{\′},y^{\′})+t)]}^{1/2}},---\left(11\right).$

4.如权利要求1所述的方法，其中，步骤(2)包括以下步骤：

用以下多项式(12)来对所述窗口内地震信号的振幅进行多项式拟合：

A(x，y)＝b₀₀+b₁₀x+b₁₁y+b₂₀x²+b₂₁xy+b₂₂y²+b₃₀x³+b₃₁x²y+b₃₂xy²+b₃₃y³      (12)

在所述窗口内，各点的输出振幅通过对原始道振幅分别沿纵横测线的同相轴方向按距离加权叠加得到，如下：

$A_{\mathrm{xy}}=\underset{i=-M}{\overset{M}{\mathrm{\Σ}}}\underset{j=-N}{\overset{N}{\mathrm{\Σ}}}{W}_{\mathrm{xyij}}{B}_{\mathrm{ij}},---\left(13\right)$

其中，A_xy表示坐标为(x，y)点的地震道的输出振幅，W_xyij为点(i，j)对点(x，y)的加权系数，B_ij表示点(i，j)的原始振幅，

根据最小二乘法拟合出振幅多项式(12)的系数b₀₀，b₁₀，b₁₁，......。

5.如权利要求4所述的方法，其中，当多项式次数较高时，将振幅多项式(12)改写成：

A(x，y)＝d₀₀p₀(x)q₀(y)+d₁₀p₁(x)+d₁₁q₁(y)+d₂₀p₂(x)+d₂₁p₁(x)q₁(y)+d₂₂q₂(y)+......(14)

其中，p₀，q₀，p₁，q₁......为等式(4)所描述的正交多项式，采用最小二乘法确定系数d₀₀，d₁₀，d₁₁......。

6.如权利要求4所述的方法，其中，W_xyij为点(i，j)和点(x，y)的距离d((x，y)，(i，j))的倒数函数。

7.如权利要求1所述的方法，其中，步骤(6)包括以下步骤：

对于窗口沿时间方向移动的情况，把窗口分成上下两个窗口，将上窗口中拟合的波形乘以权系数，其中，上窗口的权系数从上到下由1按样点数线性地减小到0，将下窗口中拟合的波形也乘以与上窗口的权系数的大小的方向相反的权系数，将得到的两个新波形相加来得到重叠部分波形；

对于窗口沿纵横测线方向移动的情况，把窗口分成左右两个窗口，加权方式采用能量加权平均方式，其某道的权系数由左右两个窗口在该道的能量值决定，左窗口某重叠道权系数为左窗口在该道的能量比前一个右窗口在该道的能量与左窗口在该道的能量之和，右窗口某重叠道权系数为右窗口在该道的能量比后一个左窗口在该道的能量与右窗口在该道的能量之和，将得到的两个新波形相加即得到重叠部分波形。

8.如权利要求1所述的方法，其中，步骤(8)包括：

利用等式(15)将拟合出来的叠加数据与原始输入数据进行混波处理：

S′(x，y，t)＝S(x，y，t)×P+T(x，y，t)(1-P)               (15)

其中，S′(x，y，t)为输出记录，S(x，y，t)为原始输入记录，T(x，y，t)为拟合记录，P为混波比。

Images