CN102681014B - 基于多项式拟合的规则线性干扰压制方法 - Google Patents

Abstract

提供了一种基于多项式拟合的规则线性干扰压制方法，包括：(a)人工识别规则线性干扰的视速度范围；(b)利用小波变换对采集的地震数据进行分频处理为低频数据和高频数据；(c)根据噪声数据的视速度范围，确定噪声数据的视倾角范围，进行视倾角最佳中值扫描，确定当前视倾角下的规则线性干扰的最佳方向；(d)沿规则线性干扰的最佳方向对噪声数据进行多项式拟合；(e)重复步骤(c)和(d)，直到低频数据中的所有规则线性干扰被拟合出来；(f)对分频出来的原始低频数据减去拟合出来的规则干扰数据，得到压制规则干扰后的低频数据；(g)把压制规则干扰后的低频数据和分频后的原始高频数据，最终得到经过多项式拟合规则噪声压制后的地震数据。

Description

基于多项式拟合的规则线性干扰压制方法

技术领域

本发明涉及石油地震勘探,属于地震勘探资料处理与解释领域，更具体地讲，本发明涉及一种规则线性干扰压制方法。

背景技术

复杂山地地震资料噪声来源十分复杂，噪声对有效波有较强的干扰，尤其是地震资料中普遍存在着相当严重的规则线性干扰，如浅层强能量的地表多次折射波和规则线性干扰，中、深层的强声波干扰，以及频带很宽的地表直达波等。一般情况下，这些干扰破坏了有效反射信号，严重时覆盖了整个地震记录，完全淹没了有效信号，大大降低了地震资料的信噪比。对这些干扰，常规的处理方法是采用f-k滤波法，τ-p变换，中值滤波法，速度陷波法等。这些方法虽然对规则线性干扰的消除有一定的作用，但它们本身都存在一定的局限性。如f-k滤波技术要求均匀的空间采样，而陆地数据尤其3D数据常常难以满足这一苛刻条件；经典的τ-p变换也存在同样的问题；速度陷波滤波器则要求满足噪声为线状噪声的这一条件。这些方法都对道间振幅的变化相当敏感，即使细微的变化也会使滤波的效果变得极差。所以如果这些规则线性干扰不能得到很好的压制，最终会影响地震记录的叠加成像质量及地震剖面的横向分辨率。

因此，需要一种能够很好地压制这些规则线性干扰的方法。

发明内容

本发明提供了一种能够针对复杂山地低信噪比资料中的规则线性干扰压制方法，从而提高地震数据的处理质量。

为了实现上述目的，提供了一种基于多项式拟合的规则线性干扰压制方法，包括：(a)采集地震数据，并人工识别规则线性干扰的视速度范围和频率范围。(b)利用小波变换对采集的地震数据进行分频处理，从而地震数据被分为低频地震数据和高频地震数据；(c)根据噪声数据的视速度范围，确定噪声数据的视倾角范围，对分频后的低频地震数据在某个视倾角下进行视倾角最佳中值扫描，确定当前视倾角下的规则线性干扰的最佳方向；(d)沿规则线性干扰的最佳方向对噪声数据进行多项式拟合，拟合出此视倾角下不含有效信号的规则干扰；(e)确定低频地震数据中的所有规则线性干扰是否都被拟合出来，如果没有都被拟合出来，则返回步骤(c),如果全部都被拟合出来，则执行步骤(f)；(f)对分频出来的原始低频地震数据减去拟合出来的规则干扰数据，得到压制规则干扰后的低频地震数据；(g)把压制规则干扰后的低频地震数据和分频后的原始高频地震数据，通过小波变换进行重构，最终得到经过多项式拟合规则噪声压制后的地震数据。

优选地，在步骤(b)中，利用小波变换将采集的地震数据进行分频处理为低频低波数分量、低频高波数分量、高频低波数分量和高频高波数分量，并提取其中的低频低波数分量和低频高波数分量。

优选地，视倾角最佳中值扫描基于两个假设条件：(1)在地震记录中，沿着某一个视倾角方向取一组样点值组成一个序列{x_i,t}，其中i为道号,t为时间，如果{x_i,t}中只有随机噪声，当i足够大时,则{x_i,t}的中值为零；(2)如果在{x_i,t}中既有相干信号，又有随机噪声，则{x_i,t}的中值是该序列取值方向的相干信号。

优选地，根据下述等式(1)来确定出视速度下的视倾角范围：视倾角＝(1000×道间距)/视速度  (1)。

优选地，当视倾角范围确定后在指定的视倾角范围内进行扫描，求取一个中值序列{m_j}，其中j为扫描序号，通过下面的等式(2)来从中值序列{m_j}中找到最佳中值M₀：

c₀＝max{c_j}   (2)

其中，c_j是m_j与{x_i,t}的相关系数，由{m_j}与{x_i,t}可获得一个相关系数序列{c_j}，从而求取{c_j}的最大值c₀，c₀所对应的中值就是要求取的最佳中值M₀。

附图说明

图1是根据本发明的基于多项式拟合的规则线性干扰压制方法的流程图。

图2是描述根据本发明第一示例性实施例的基于多项式拟合的规则线性干扰压制方法运用于含一个规则线性干扰、两层有效同相轴的理论模型的示例的示图。

图3是描述根据本发明第二示例性实施例的基于多项式拟合的规则线性干扰压制方法运用于二维叠前实际资料的示图。

图4是描述根据本发明第三示例性实施例的基于多项式拟合的规则线性干扰压制方法运用于三维叠前实际资料的示图。

具体实施方式

现在，详细描述本发明的实施例，其示例在附图中表示，其中，相同的标号始终表示相同的部件。以下通过参考附图描述实施例以解释本发明。

多项式拟合技术是假设信号在各道上出现时间符合一个待定的多项式，而不采用时间相同或线性变化的假设，同时假设各道的振幅变化也可用一个多项式来近似，这样相关性较好的规则线性干扰就可以通过一个多项式拟合出来，从而可以进行规则线性干扰压制处理，原始数据与拟合出的规则线性干扰数据之差就得到去噪后的数据，这样处理后的数据尽可能地保持了原始的分辨率，数据的信噪比有很大提高，且数据的高频成分不受损失，能保持原有信号的分辨率，同时也能保持原始各道的相对振幅，具有较高的保真度。

图1是根据本发明的基于多项式拟合的规则线性干扰压制方法的流程图，该方法包括以下步骤：

在步骤S101，采集地震数据，并人工识别规则线性干扰的视速度范围和频率范围。

在步骤S102，利用小波变换对采集的地震数据进行分频处理，从而地震数据被分为低频地震数据和高频地震数据；

在步骤S103，根据噪声数据的视速度范围，确定噪声数据的视倾角范围，对分频后的低频地震数据在某个视倾角下进行视倾角最佳中值扫描，确定当前视倾角下的规则线性干扰的最佳方向；

在步骤S104，沿规则线性干扰的最佳方向对噪声数据进行多项式拟合，拟合出此视倾角下不含有效信号的规则干扰；

在步骤S105，确定低频地震数据中的所有规则线性干扰是否都被拟合出来。如果没有都被拟合出来，则返回步骤S103。如果全部都被拟合出来，则执行步骤S106。

在步骤S106，对分频出来的原始低频地震数据减去拟合出来的规则干扰数据，得到压制规则干扰后的低频地震数据。

在步骤S107，把压制规则干扰后的低频地震数据和分频后的原始高频地震数据，通过小波变换进行重构，最终得到经过多项式拟合规则噪声压制后的地震数据。

下面将分别对上述各步骤进行详细描述。

对步骤S101(即采集地震数据，并人工识别规则线性干扰的视速度范围和频率范围)和步骤S102(即，利用小波变换对采集的地震数据进行分频处理，从而地震数据被分为低频地震数据和高频地震数据)进行详细描述。

叠前规则线性干扰如地表折射多次波、规则线性干扰和声波及直达波等，由于传播机理各不相同，使得它们相互之间以及与有效反射波之间存在诸多方面的差异，其显著特征之一就是频谱的衰减规律有较大不同，即不同的干扰波，其优势频段各不相同。因此，针对不同的规则线性干扰，在各自的优势频带范围内进行识别，不仅可以更准确地识别出规则线性干扰波，而且可以保证其它频段内的有效信号不受影响，使处理结果具有更高的保真度，而规则线性干扰一般集中于地震数据所有频率中的低频部分。

小波变换是一种时间窗和频率窗都可以改变的时频局部化分析方法。具有多分辨率的特点，即在低频部分具有较高的频率分辨率和较低的时间分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率。由于小波变换的频域和时域可调，可以解决傅立叶变换不能解决的问题，同时小波变换具备时域—频域的双重良好的局部性和随尺度变化的自动调焦功能，使得它在叠前分频处理方面具有独特的优势。

利用小波变换可以将x-t域的地震记录转换到二维小波变换域，即时间、频率、空间和波数四维域中，得到4个小波系数分量，即低频低波数分量、低频高波数分量、高频低波数分量和高频高波数分量。因为规则线性干扰主要集中在低频低波数分量和低频高波数分量中，因此只需对低频低波数分量和低频高波数分量做规则干扰压制处理，然后进行信号重构，即可消除信号中的规则线性干扰于扰。如果一次分解没有将规则线性干扰完全分离出来，可以继续进行分解。此时只需对第一次分解的低频低波数分量再进行二维小波分解，同样分离出4个小波系数分量，将低频低波数分量和低频高波数分量做规则干扰压制处理，然后与第一次分解得到的其余小波系数分量进行信号重构，即获得了消除规则线性干扰的记录。理论上，分离规则线性干扰的迭代过程可以无限进行下去，直到规则线性干扰被完全分离为止。

通过小波变换就可以找出数据中含规则线性干扰的低频数据，而高频数据部分是只含有效信号的数据，不进行噪声处理，从而有效的保护了高频部分的有效信号。

下面对步骤S103(即，根据噪声数据的视速度范围，确定噪声数据的视倾角范围，对分频后的低频地震数据在某个视倾角下进行视倾角最佳中值扫描，确定当前视倾角下的规则线性干扰的最佳方向)进行详细描述。

在地震记录中，有效反射波以及各种规则干扰波具有很强的相干性，都可看成是一种相干信号，在理想情况下，我们顺着某一组相干信号的同相轴取数，获得一个样点序列{x_i}，在有限区间内，{x_i}可以看成是一个常数序列。事实上，在地震记录中，不仅存在随机噪声，而且还存在许多不同视倾角的规则干扰信号。如果我们只希望预测某一个视倾角的规则干扰信号，那么其它规则干扰信号相对该视倾角规则干扰信号均可看成是随机噪声。

视倾角最佳中值扫描法基于两个假设条件：(1)在地震记录中，沿着某一个视倾角方向取一组样点值组成一个序列{x_i,t}(其中i为道号,t为时间)，如果{x_i,t}中只有随机噪声，当i足够大时,则{x_i,t}的中值为零。(2)如果在{x_i,t}中既有相干信号，又有随机噪声，则{x_i,t}的中值是该序列取值方向的相干信号。

在地震记录中，我们所希望预测的规则干扰信号的视倾角是千变万化的，让计算机自动地沿着我们所希望的视倾角方向取序列{x_i,t}是一个难以解决的技术问题。但是通过人工识别的视速度范围根据下述等式(1)来确定出视速度下的视倾角范围：

视倾角＝(1000×道间距)/视速度  (1)

当视倾角范围确定后让计算机在指定的视倾角范围内进行扫描，求取一个中值序列{m_j}(其中j为扫描序号)。并从该中值序列{m_j}中获得最佳中值，该最佳中值对应的倾角即为最佳视倾角。下面描述如何从中值序列{m_j}中找到最佳中值M₀。

由条件(2)(即，如果在{x_i,t}中既有相干信号，又有随机噪声，则{x_i,t}的中值是该序列取值方向的相干信号)可知：M₀就是我们所求取的规则干扰信号，即，M₀与{x_i,t}的相关系数达到最大。因此，我们可以用求取最大相关系数的方法来从{m_j}中找到M₀。

设c_j是m_j与{x_i,t}的相关系数，由{m_j}与{x_i,t}可获得一个相关系数序列{c_j}。求取{c_j}的最大值c₀，即：

c₀＝max{c_j}   (2)

c₀所对应的中值就是要求取的最佳中值M₀。

最佳中值M₀所在的规则线性干扰方向是当前视速度下规则线性干扰的最佳方向，即规则线性干扰所在的方向，但此规则线性干扰里面含有有效波，所以还必须进行后续处理。

下面对步骤S104(即，沿规则线性干扰的最佳方向对噪声数据进行多项式拟合，拟合出此视倾角下不含有效信号的规则干扰)进行详细描述。

在地震数据中进行拟合技术处理中常采用正交多项式拟合方法，这种拟合方法的拟合结果存在拟合不够准确的缺点，从而导致处理后地震数据的主频降低、短点变形，如把此方法运用与叠前规则线性干扰的拟合上，则不能准确、完整地拟合出规则干扰，从而影响压噪效果，究其原因，是正交多项式拟合过程中拟合中心点随正交多项式系数的改变而改变，从而导致拟合不够准确。

为了克服上述缺陷，这里采用二次多项式进行取代正交多项式进行规则干扰拟合，然后采用先线性拟合再微调的办法予以解决针对扫描过程中拟合系数不能直接采用最优化算法的问题，最后采取了先大步长等间距扫描确定最佳拟合系数的单极值区间，然后再用二分法进行最佳拟合系数扫描的措施,节省了计算时间针对拟合过程中出现断点模糊这一问题。

在某一种数据集上,可用下面的正交多项式来描述地震波的到达相位时间的表达式：

T(x)＝a₀+a₁p₁(x)+a₂p₂(x)

其中， $p_1\left(x\right)=x,p_2=x^2-M(M+1)/3,\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}{p}_2\left(x\right)p_1\left(x\right)=0---\left(3\right)$

在上述等式(3)中：T为拟合时间函数；a₁,a₂分别为一次和二次多项式的系数，a₀为初始拟合时窗的中心位置；x为观测点横坐标，且x＝-M,-M+1,......,-1,0,1,......M,M取拟合道数的1/2,p₁(x),p₂(x)为彼此正交的正交多项式。

由上述等式(3)可见,二次多项式含有常数项，拟合时窗的中心随a₂的变化而变化,因此只能假定正交二次多项式表示的是某一时空域中地震数据的整个趋势。为了避免每个时窗分界处由于采用不同的拟合系数而产生突变,实际处理中采取了滑动时窗的办法，但这种处理方式由于地震数据本身具有的时变空变性而显得较为粗糙。

如果二次多项式p₂(x)中不含常数项,便可以解决拟合中心点随二次项系数变化问题。这里采用二次拟合多项式表达相位时间：

t(x)＝t₀+a₁x+a₂x²

其中，x＝-M,-M+1,......,-1,0,1,......M,(4)

上述等式(4)中t₀为拟合中心点的时间，x为相对拟合中心道的道号；M为拟合道数的一半。由于等式(4)中的x,x²不正交，a₁,a₂不能通过正交多项式的特点独立扫描,从而使扫描的计算量增加。但是，我们要拟合的规则线性干扰是呈线性或拟线性的,而且我们同样可以通过视倾角最佳中值扫描法得到一次项系数,然后在此基础上再应用二次项系数扫描得到某点的最佳二次拟合多项式。具体实现方法是首先令a₂＝0，先扫描求取a₁，然后代人式(4),再扫描求取a₂,此时对a₂的扫描只是在a₁确定情况下的一个微调。这样二次多项式拟合出来的2m+1个时间值是以t₀和x＝0为拟合点的拟合结果,通过它可以实现对某点的准确拟合,而且计算量与正交二次多项式相同。

根据通过视倾角最佳中值扫描法可以确定拟合的初始扫描范围，但还需要通过下式(5)确定最佳拟合系数。最佳拟合系数的确定采用相似性最大原则，假设一次项系数a₁的扫描范围为[N₁,N₂]，扫描步长为Δa₁，那么a₁(i)＝N₁+i×Δa₁(i＝0,1,......(N₂-N₁)/Δa₁)为扫描范围内的具体扫描值。由相似准则有：

$R\left(i\right)=\underset{k,j=-m}{\overset{m}{\mathrm{\Σ}}}\frac{\underset{t=-l/2}{\overset{l/2}{\mathrm{\Σ}}}{s}_k(t_0+a_1\left(i\right)\×k+t)\×s_j(t_0+a_1\left(i\right)\×j+t)}{{[\underset{t=-l/2}{\overset{l/2}{\mathrm{\Σ}}}{s}_k^2(t_0+a_1\left(i\right)\×k+t)\×s_j^2(t_0+a_1\left(i\right)\×j+t)]}^{1/2}}---\left(5\right)$

在等式(5)中：R(i)为拟合道互相关之和；s_k,s_j为拟合道的样点值；l为相关时窗长；t₀为拟合点的时间；k和j为相对拟合道号；2m+1为拟合道数。

对于不同的a₁(i),使R(i)为最大的a₁(i),即为一次项系数的最佳拟合系数。用该最佳拟合系数进行地震数据拟合显然十分粗糙，所以还需确定扫描步长使得最佳拟合系数能够做到精细扫描。

由于在初始范围内计算相关函数R(i)时可能会产生串相位问题，因此不能确定所有参与扫描的a₁(i)使R(i)一定为单极值函数。但是,在初始范围内的某一局部区域不会发生串相位现象,即R(i)为单极值函数。也就是说在使R(i)为最大的以a₁(i)为中心的一个局部区域为单极值函数。然后以(a₁(i)-1,a₁(i)+1)为新的扫描区间，采用二分法进行新一轮的a₁(i)扫描。扫描步长为以下的等式(6)：

在等式(6)中，Δa₁为本次扫描步长；为前一次扫描步长。显然在扫描区间共有5个扫描值，且3个扫描值是上一轮扫描所得，这一轮仅对2个a₁(i)扫描即可。将这一轮扫描中使得R(i)为最大的a₁(i)定为新的扫描中心，再由[a₁(i)-Δa₁,a₁(i)+Δa₁]为新的扫描范围，采用等式(6)式确定新的扫描步长，进行下一轮扫描。

当时间多项式确定以后，便可对其所确定的时窗内地震信号的振幅来进行多项式拟合。可用下列振幅多项式(7)来表示：

A(x)＝b₀+b₁x+b₂x²+b₃x²   (7)

要求取上式的系数b₀,b₁,......，首先求对应时窗内各道记录的均方根振幅A(x_n)，A(x_n)的表达式为下述等式(8)：

$A\left(x_n\right)=\sqrt{\frac{\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}{S}^2(x_n,t_n+t)}{2L+1}}---\left(8\right)$

在等式(8)中,L为时间方向上窗口移动的大小，S(x_n,t)为第n道上窗口内地震道的样点值，当得到A(x_n)之后，便可根据最小二乘法拟合出振幅多项式(7)的系数b₀,b₁,......。

当振幅多项式(7)也确定以后，接下来的工作是形成期望规则干扰波形。其具体的方法是将时间多项式(3)确定的同一时窗内2N+1的道记录波形沿拟合出的同相轴方向相加，并对相加结果缩放，使其均方根振幅归一，波形相加计算等式(9)为：

$A_i=\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}S(n,t_n+i-L-1)$

i＝1,2,......2L+1(9)

在等式(9)中，t_n为第n道的时窗中点时间，由时间多项式确定。其归一化等式(10)为：

${\stackrel{\‾}{A}}_i=A_i/\sqrt{\frac{\underset{i=1}{\overset{2L+1}{\mathrm{\Σ}}}{A}_i^2}{2L+1}}---\left(10\right)$

当时间多项式、振幅多项式都已经确定，且期望波形形成之后，就可以形成期望剖面，即最终的拟合结果。具体实现如下：

首先将按振幅多项式求出的各道振幅值，乘上当前窗口内计算出的标准波形，即得窗口内各道的拟合波形，然后将得到的拟合波形放到时间多项式计算出的位置上便形成了当前视倾角最佳方向上的规则干扰。

在拟合出当前视倾角最佳方向上的规则干扰之后，在步骤S105确定低频地震数据中的所有规则线性干扰是否都被拟合出来。如果没有，则继续执行步骤S103和S104，换句话说，需要把低频数据中的视倾角范围内的所有规则线性干扰都拟合出来。

之后，在步骤S106，对分频出来的原始低频地震数据减去拟合出来的规则干扰数据，得到压制规则干扰后的低频地震数据。在步骤S107，把压制规则干扰后的低频地震数据和分频后的原始高频地震数据，通过小波变换进行重构，最终得到经过多项式拟合规则噪声压制后的地震数据。

本发明可以压制二维地震数据的规则线性干扰，而且可以压制三维地震数据的规则线性干扰。于三维叠前数据来说，由于三维野外采集时炮点不在检波线上，道间距是变化的，在小偏移距上规则线性干扰和直达波不呈线性变化，这时可以根据偏移距来计算规则线性干扰的视倾角，才能得到正确的结果。

下面参照图2来描述根据本发明第一示例性实施例的基于多项式拟合的规则线性干扰压制方法，其中该方法运用于含一个规则线性干扰、两层有效同相轴的理论模型。从图1中的(a)可以看到规则线性干扰与第一层有效同相轴耦合在一起，严重降低了实际资料的信噪比。

根据图1中的(a)，先分析地震数据的道头信息，考虑观测系统对规则噪声分布形态的影响，确定规则线性干扰具有的特点，判断出规则线性干扰所在的大致速度范围(v₁,v₂)，由于只是单一的规则线性干扰和有效信号，所以可不需要做分频处理，直接进行视倾角最佳中值扫描。具体实施如下：

由于我们不知道所要需找的规则线性干扰的具体速度，只能人为的判断规则线性干扰所在的一个速度范围(v₁,v₂)，我们利用等式(1)计算其视倾角范围(α₁,α₂)。

视倾角＝(1000×道间距)/视速度  (1)

这里我们的角度增量按Δα，对每一个视倾角进行递增扫描，而对于我们要扫描的每一个视倾角必须要找到他的最佳方向才是这个规则线性干扰所在方向的速度。这里采用最佳中值扫描法来确定视速度方向：在指定的视倾角范围内(α₁,α₂)进行扫描，先在初始扫描角度α₁所在方向求取一个中值序列{m_j}(其中j为扫描序号)。那么在{m_j}中总会有一个中值M₀所对应的倾角是我们所希望的视倾角，我们称M₀为最佳中值。为了找到这个最佳中值M₀，我们可以用求取最大相关系数的方法来从{m_j}中找到M₀，如等式(2)。

设c_j是m_j与{x_i,t}的相关系数，由{m_j}与{x_i,t}可获得一个相关系数序列{c_j}。求取{c_j}的最大值c₀，即：

c₀＝max{c_j}   (2)

c₀所对应的中值就是我们要求取的最佳中值M₀。

由于理论模型中只有一个线性干扰，如果在此角度下没有找到那个最佳中值，则说明此角度下无线性干扰，那么通过倾角增量，进行下一个倾角扫描，直到找到理论模型线性干扰所对应的真正倾角a_m,[a_m∈(a₁,a₂)]，那么通过等式(2)也必能找到此倾角a_m下的最佳中值M₀，即规则干扰下的最佳方向。此时就可扫描出当前视倾角a_m下的规则线性干扰，但此规则线性干扰里面含有有效波，所以必须进行后续处理。

然后沿规则线性干扰的最佳方向上对噪声数据进行多项式拟合处理，拟合出此视倾角下不含有效信号的规则干扰。具体实施如下：

以时窗扫描的形式，利用式(4)对线性干扰进行相位时间的二次多项式扫描，

t(x)＝t₀+a₁x+a₂x²   (4)

上式中的系数扫描通过式(5)相似性最大原则进行判断以找到准确的系数值a₁,a₂。

$R\left(i\right)=\underset{k,j=-m}{\overset{m}{\mathrm{\Σ}}}\frac{\underset{t=-l/2}{\overset{l/2}{\mathrm{\Σ}}}{s}_k(t_0+a_1\left(i\right)\×k+t)\×s_j(t_0+a_1\left(i\right)\×j+t)}{{[\underset{t=-l/2}{\overset{l/2}{\mathrm{\Σ}}}{s}_k^2(t_0+a_1\left(i\right)\×k+t)\×s_j^2(t_0+a_1\left(i\right)\×j+t)]}^{1/2}}---\left(5\right)$

为使最佳拟合系数进行地震数据拟合时做到精细扫描，还需要通过式(6)确定扫描步长。

当时间多项式(4)式确定以后，便可用式(7)对其所确定的时窗内地震信号的振幅来进行多项式拟合。

A(x)＝b₀+b₁x+b₂x²+b₃x²   (7)

要求取上式的系数b₀,b₁,......，需利用式(8)求对应时窗内各道记录的均方根振幅A(x_n)：

$A\left(x_n\right)=\sqrt{\frac{\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}{S}^2(x_n,t_n+t)}{2L+1}}---\left(8\right)$

当得到A(x_n)之后，便可根据最小二乘法拟合出振幅多项式(7)的系数b₀,b₁,......。

当振幅多项式(7)也确定以后，接下来的工作是形成期望规则干扰波形。其具体的方法是将时间多项式(3)确定的同一时窗内2N+1的道记录波形沿拟合出的同相轴方向按式

(9)相加：

$A_i=\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}S(n,t_n+i-L-1)$

i＝1,2,......2L+1   (9)

并对相加结果缩放，按式(10)使其均方根振幅归一：

${\stackrel{\‾}{A}}_i=A_i/\sqrt{\frac{\underset{i=1}{\overset{2L+1}{\mathrm{\Σ}}}{A}_i^2}{2L+1}}---\left(10\right)$

当时间多项式、振幅多项式已经确定，且期望波形形成之后，就可以形成期望剖面，即最终的拟合结果。具体实施如下：

首先将按振幅多项式求出的各道振幅值，乘上当前窗口内计算出的标准波形，即得窗口内各道的拟合波形，然后将得到的拟合波形放到时间多项式计算出的位置上便形成了当前视倾角a_m最佳方向上的规则干扰，如图2中的(b)所示。

最后，用原始数据减去拟合出来的规则干扰数据，得到压制规则干扰后的数据，如图2中的(c)所示，最终得到经过多项式拟合规则噪声压制后的数据，从而完成整个规则干扰压制操作过程。

压制效果如图2中的(a)-(c)所示，基于多项式拟合的规则线性干扰压制方法，可以完全拟合出规则线性干扰，而且得到的有效信号基本没有损伤，整个理论模型展示了此方法具有良好的去噪效果。

下面参照图3来描述根据本发明第二示例性实施例的基于多项式拟合的规则线性干扰压制方法，其中该方法运用于二维叠前实际资料，从图3中的(a)可以看到有很强的规则线性干扰，这些规则线性干扰与有效信号耦合在一起，严重降低了实际资料的信噪比。

根据图3中的(a)，先分析地震数据的道头信息，考虑观测系统对规则噪声分布形态的影响，确定规则线性干扰具有的特点，判断出规则线性干扰所在的大致速度范围(v₁,v₂)，由于实际资料中规则线性干扰主要分布在低频部分，所以必须先对数据做分频处理。由于小波变换有无限细分的优越性，这里选用小波变换进行分频处理,具体实施如下：

利用小波变换将图3中的(a)中时空域的地震记录转换到二维小波变换域，即：时间、频率、空间和波数四维域中，得到4个小波系数分量，即：低频低波数分量、低频高波数分量、高频低波数分量和高频高波数分量。因为规则线性干扰主要集中在低频低波数分量和低频高波数分量中，因此只需对低频低波数分量和低频高波数分量做后续规则干扰压制处理，然后把规则干扰消除后的低频分量与高频分量进行信号重构，即可得到消除了信号中的规则线性干扰的数据。如果一次分解没有将规则线性干扰完全分离出来，可以继续进行分解。此时只需对第一次分解的低频低波数分量再进行二维小波分解，同样分离出4个小波系数分量，将低频低波数分量和低频高波数分量做规则干扰压制处理，然后与第一次分解得到的其余小波系数分量进行信号重构，即获得了消除规则线性干扰的记录。理论上，分离规则线性干扰的迭代过程可以无限进行下去，直到规则线性干扰被完全分离为止。

当通过小波变换分离出含规则干扰的低频数据后，下面进行视倾角最佳中值扫描，以此找到每个视倾角下的规则线性干扰的最佳所在方向，具体实施如下：

由于我们不知道所要需找的规则线性干扰的具体速度，只能人为的判断规则线性干扰所在的一个速度范围(v₁,v₂)，但我们利用等式(1)计算其视倾角范围(α₁,α₂)。

视倾角＝(1000×道间距)/视速度  (1)

这里我们的角度增量按Δα，对每一个视倾角进行递增扫描，而对于我们要扫描的每一个视倾角必须要找到他的最佳方向才是这个规则线性干扰所在方向的速度。这里采用最佳中值扫描法来确定视速度方向：在指定的视倾角范围内(α₁,α₂)进行扫描，先在初始扫描角度α₁所在方向求取一个中值序列{m_j}(其中j为扫描序号)。那么在{m_j}中总会有一个中值M₀所对应的倾角是我们所希望的视倾角，我们称M₀为最佳中值。为了找到这个最佳中值M₀，我们可以用求取最大相关系数的方法来从{m_j}中找到M₀，如等式(2)。

设c_j是m_j与{x_i,t}的相关系数，由{m_j}与{x_i,t}可获得一个相关系数序列{c_j}。求取{c_j}的最大值c₀，即：

c₀＝max{c_j}   (2)

c₀所对应的中值就是我们要求取的最佳中值M₀。

通过上面的操作能扫描到a₁所在规则线性干扰的最佳方向，但此规则线性干扰里面含有有效波，所以必须进行后续拟合处理，才能拟合出不含有效信号的线性干扰。

我们这里的视倾角扫描是按一定的倾角扫描增量进行扫描，直到视倾角范围内的所有线性干扰的最佳方向都扫描到。而这里的多项式拟合处理是在每一个视倾角找到其线性干扰最佳方向后就进行处理，以此重复直到拟合出所有视倾角下的线性干扰。

然后沿规则线性干扰的最佳方向上对噪声数据进行多项式拟合处理，拟合出此视倾角下不含有效信号的规则干扰，具体实施如下：

以时窗扫描的形式，利用等式(4)对线性干扰进行相位时间的二次多项式扫描，

t(x)＝t₀+a₁x+a₂x²   (4)

上式中的系数扫描通过等式(5)相似性最大原则进行判断以找到准确的系数值a₁,a₂。

$R\left(i\right)=\underset{k,j=-m}{\overset{m}{\mathrm{\Σ}}}\frac{\underset{t=-l/2}{\overset{l/2}{\mathrm{\Σ}}}{s}_k(t_0+a_1\left(i\right)\×k+t)\×s_j(t_0+a_1\left(i\right)\×j+t)}{{[\underset{t=-l/2}{\overset{l/2}{\mathrm{\Σ}}}{s}_k^2(t_0+a_1\left(i\right)\×k+t)\×s_j^2(t_0+a_1\left(i\right)\×j+t)]}^{1/2}}---\left(5\right)$

为使最佳拟合系数进行地震数据拟合时做到精细扫描，还需要通过式(6)确定扫描步长。

当时间多项式(4)式确定以后，便可用等式(7)对其所确定的时窗内地震信号的振幅来进行多项式拟合。

A(x)＝b₀+b₁x+b₂x²+b₃x²   (7)

要求取上式的系数b₀,b₁,......，需利用等式(8)求对应时窗内各道记录的均方根振幅A(x_n)：

$A\left(x_n\right)=\sqrt{\frac{\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}{S}^2(x_n,t_n+t)}{2L+1}}---\left(8\right)$

当得到A(x_n)之后，便可根据最小二乘法拟合出振幅多项式(7)的系数b₀,b₁,......。

当振幅多项式(7)也确定以后，接下来的工作是形成期望规则干扰波形。其具体的方法是将时间多项式(3)确定的同一时窗内2N+1的道记录波形沿拟合出的同相轴方向按等式(9)相加：

$A_i=\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}S(n,t_n+i-L-1)---\left(9\right)$

i＝1,2,......2L+1

并对相加结果缩放，按等式(10)使其均方根振幅归一：

${\stackrel{\‾}{A}}_i=A_i/\sqrt{\frac{\underset{i=1}{\overset{2L+1}{\mathrm{\Σ}}}{A}_i^2}{2L+1}}---\left(10\right)$

当时间多项式、振幅多项式已经确定，且期望波形形成之后，就可以形成期望剖面，即最终的拟合结果。具体实现如下：

首先将按振幅多项式求出的各道振幅值乘上当前窗口内计算出的标准波形，即得窗口内各道的拟合波形，然后将得到的拟合波形放到时间多项式计算出的位置上便形成了当前视倾角a_m最佳方向上的规则干扰，如图3中的(b)所示。

最后，用原始数据减去拟合出来的规则干扰数据，得到压制规则干扰后的数据，如图3中的(c)所示，最终得到经过多项式拟合规则噪声压制后的数据，从而完成整个规则干扰压制操作过程。

压制效果如图3中的(a)-(c)所示，各组规则线性干扰得到了较好的滤除，原始炮集记录的信噪比得到明显的改善。

下面参照图4来描述根据本发明第三示例性实施例的基于多项式拟合的规则线性干扰压制方法，其中该方法运用于三维叠前实际资料。

对于三维实际地震数据，如图4中的(a)所示，可以看到有很强的规则线性干扰，这些规则线性干扰与有效信号耦合在一起，严重降低了实际资料的信噪比。对于三维资料的处理和二维资料的处理方式大致相同，不同的是由于三维野外采集时炮点不在检波线上，道间距是变化的，在小偏移距上面波和直达波不呈线性变化，这时需根据偏移距来计算规则线性干扰的视倾角。从压制效果图4中的(b)和图4中的(c)来看基于多项式拟合的规则线性干扰得到了预期的滤波效果，线性噪音得到更好的压制，同相轴的连续性也更好，对有效信号的保护也较好，能更好地保持实际地震记录所反映的地层信息和地质特征。

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

Claims (5)

1.一种基于多项式拟合的规则线性干扰压制方法，包括：

(a)采集地震数据，并人工识别规则线性干扰的视速度范围和频率范围；

(b)利用小波变换对采集的地震数据进行分频处理，从而地震数据被分为低频地震数据和高频地震数据；

(c)根据噪声数据的视速度范围，确定噪声数据的视倾角范围，对分频后的低频地震数据在某个视倾角下进行视倾角最佳中值扫描，确定当前视倾角下的规则线性干扰的最佳方向；

(d)沿规则线性干扰的最佳方向对噪声数据进行多项式拟合，拟合出此视倾角下不含有效信号的规则干扰；

(e)确定低频地震数据中的所有规则线性干扰是否都被拟合出来，如果没有都被拟合出来，则返回步骤(c),如果全部都被拟合出来，则执行步骤(f)；

(f)对分频出来的原始低频地震数据减去拟合出来的规则干扰数据，得到压制规则干扰后的低频地震数据；

(g)把压制规则干扰后的低频地震数据和分频后的原始高频地震数据，通过小波变换进行重构，最终得到经过多项式拟合规则噪声压制后的地震数据，

其中，步骤(d)包括如下步骤：

在某一种数据集上,用下面的正交多项式(3)来描述地震波的到达相位时间的表达式：

T(x)＝a₀+a₁p₁(x)+a₂p₂(x)

其中， $p_1\left(x\right)=x,p_2\left(x\right)=x^2-M(M+1)/3,\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}{p}_2\left(x\right)=0---\left(3\right)$

在上述等式(3)中：T为拟合时间函数；a₁,a₂分别为一次和二次多项式的系数，a₀为初始拟合时窗的中心位置；x为观测点横坐标，且x＝-M,-M+1,......,-1,0,1,......M,M取拟合道数的1/2,p₁(x),p₂(x)为彼此正交的正交多项式，

以时窗扫描的形式，利用式(4)对线性干扰进行相位时间的二次多项式扫描，

t(x)＝t₀+a₁x+a₂x²

其中，x＝-M,-M+1,......,-1,0,1,......M,           (4)

其中，t₀为拟合中心点的时间，x为相对拟合中心道的道号；M为拟合道数的一半，x,x²不正交，a₁,a₂分别为一次和二次多项式的系数，

其中，具体实现方法是首先令a₂＝0，先扫描求取a₁，然后代人式(4),再扫描求取a₂,此时对a₂的扫描只是在a₁确定情况下的一个微调，这样二次多项式拟合出来的2m+1个时间值是以t₀和x＝0为拟合点的拟合结果,通过它实现对某点的准确拟合,而且计算量与正交二次多项式相同，根据通过视倾角最佳中值扫描法确定拟合的初始扫描范围，但还需要通过下式(5)确定最佳拟合系数，最佳拟合系数的确定采用相似性最大原则，假设一次项系数a1的扫描范围为[N₁,N₂]，扫描步长为Δa₁，那么a₁(i)＝N₁+i×Δa₁为扫描范围内的具体扫描值，i＝0,1,......(N₂-N₁)/Δa₁；

等式(4)中的系数扫描通过等式(5)相似性最大原则进行判断以找到准确的系数值a₁,a₂：

$R\left(i\right)=\underset{k,j=-m}{\overset{m}{\mathrm{\Σ}}}\frac{\underset{t=-l/2}{\overset{l/2}{\mathrm{\Σ}}}{s}_k(t_0+a_1\left(i\right)\×k+t)\×s_j(t_0+a_1\left(i\right)\×j+t)}{{[\underset{t=-l/2}{\overset{l/2}{\mathrm{\Σ}}}{s}_k^2(t_0+a_1\left(i\right)\×k+t)\×s_j^2(t_0+a_1\left(i\right)\×j+t)]}^{1/2}}---\left(5\right)$

在等式(5)中：R(i)为拟合道互相关之和；s_k,s_j为拟合道的样点值；l为相关时窗长；t₀为拟合中心点的时间；k和j为相对拟合道号；2m+1为拟合道数，

为使最佳拟合系数进行地震数据拟合时做到精细扫描，需要通过等式(6)确定扫描步长：

在等式(6)中，Δa₁为本次扫描步长；为前一次扫描步长，

当等式(4)确定以后，通过等式(7)对其所确定的时窗内地震信号的振幅来进行多项式拟合：

A(x)＝b₀+b₁x+b₂x²+b₃x²           (7)

要求取等式(7)中的系数b₀,b₁,......，需利用等式(8)求对应时窗内各道记录的均方根振幅A(x_n)：

$A\left(x_n\right)=\sqrt{\frac{\underset{t=-L}{\overset{L}{\mathrm{\Σ}}}{S}^2(x_n,t_n+t)}{2L+1}}---\left(8\right)$

在等式(8)中,L为时间方向上窗口移动的大小，S(x_n,t_n+t)为第n道上窗口内地震道的样点值，

当得到A(x_n)之后，便根据最小二乘法拟合出振幅多项式(7)的系数b₀,b₁,......，

当振幅多项式(7)也确定以后，形成期望规则干扰波形，其中，形成期望规则干扰波形的步骤是将时间多项式(3)确定的同一时窗内2N+1的道记录波形沿拟合出的同相轴方向按等式(9)相加：

$A_i=\underset{n=-N}{\overset{N}{\mathrm{\Σ}}}S(n,t_n+i-L-1)$

i＝1,2,......2L+1     (9)

在等式(9)中，t_n为第n道的时窗中点时间，由时间多项式确定，

对相加结果缩放，按等式(10)使其均方根振幅归一：

${\stackrel{\‾}{A}}_i=A_i/\sqrt{\frac{\underset{i=1}{\overset{2L+1}{\mathrm{\Σ}}}{A}_i^2}{2L+1}}---\left(10\right).$

2.如权利要求1所述的规则线性干扰压制方法，其中，在步骤(b)中，利用小波变换将采集的地震数据进行分频处理为低频低波数分量、低频高波数分量、高频低波数分量和高频高波数分量，并提取其中的低频低波数分量和低频高波数分量。

3.如权利要求1所述的规则线性干扰压制方法，其中，视倾角最佳中值扫描基于两个假设条件：(1)在地震记录中，沿着某一个视倾角方向取一组样点值组成一个序列{x_i,t}，其中i为道号,t为时间，如果{x_i,t}中只有随机噪声，当i足够大时,则{x_i,t}的中值为零；(2)如果在{x_i,t}中既有相干信号，又有随机噪声，则{x_i,t}的中值是该序列取值方向的相干信号。

4.如权利要求3所述的规则线性干扰压制方法，其中，根据下述等式(1)来确定出视速度下的视倾角范围：

视倾角＝(1000×道间距)/视速度     (1)。

5.如权利要求4所述的规则线性干扰压制方法，其中，当视倾角范围确定后在指定的视倾角范围内进行扫描，求取一个中值序列{m_j}，其中j为扫描序号，通过下面的等式(2)来从中值序列{m_j}中找到最佳中值M₀：

c₀＝max{c_j}           (2)

其中，c_j是m_j与{x_i,t}的相关系数，由{m_j}与{x_i,t}可获得一个相关系数序列{c_j}，从而求取{c_j}的最大值c₀，c₀所对应的中值就是要求取的最佳中值M₀。