CN1129801C - 地震数据噪声的过滤方法 - Google Patents

Abstract

一种用于在时间间隔([1，...，T])从离散带噪声的地震信号中过滤噪声的方法被提出，所述方法包括以下步骤：利用所述带噪声的地震信号对所述时间间隔确定至少一个参考通道(x(t))作为所述噪声的估计；确定M个时间局部滤波器(w(i，t))的系数，所述滤波器构成滤波器组，M为大于或等于二的数，通过最小化代表所述滤波器组输出的误差的量度的代价函数(J(t))；并将所述滤波器组用于所述的至少一个估计(x(t))以确定所述噪声的M个滤波估计。滤波估计与时间窗(hi(i，t))相乘。

Description

地震数据噪声的过滤方法

本发明涉及从数据信号中过滤噪声的方法，特别涉及到用来减轻附加噪声对在地震勘探中采集到的数据信号的影响的噪声过滤方法。

在陆地地震勘测中，地震源引起了在地表或近地表的地震波。这些波穿过地面传播，并从地球内部的不同地层被反射，反射波可被在地表的传感器或地震检波器检测到。地震源的振动作用到地表还会产生所谓表面波或地面(ground roll)波，其通过地球的浅层进行传播。在地震检波器上，低频低速地面波的发生时间与地震勘测中感兴趣的深层反射波典型一致。地面波与反射信号的同时出现使得充分利用地震数据很困难，因为地面波经常掩盖反射波。

衰减地面波干扰从而减轻它对感兴趣的地震信号的影响的方法已经知道一些了。典型地，地震检波器不是独立被使用，而是以分-阵列或组的形式连接在一起，其被硬件连接或并在一起。这是一种与数据无关的聚束(beamforming)形式。同样，应用自适应信号处理抑制地震勘测中的地面波的尝试也被采用。

美国专利4,556,962试图通过在靠近地震源处放置一传感器检测干扰噪声来衰减从表面地震源发出的地面波。干扰噪声被换算，延时并与从较远的地震检波器阵列发出的信号相加，然后与初始振动信号互相关。该专利还建议使用自适应滤波器修正延时信号，使其更紧密对应被更远处的地震检波器阵列检测的信号。然而，靠近震源检测到的地面波可以与被地震检波器阵列接收到的信号有很大不同，自适应滤波器可能不能处理这种情况。

美国专利4,890,264中，抑制由表面波传播，风，及机械产生的非均匀分布的噪声的方法被描述。用来检测表面波的水平地震检波器与传统的用来检测地震能量的垂直地震检波器共同使用。表面波检测器的输出结合自适应滤波器被使用，以消除表面波干扰的影响。抑制地面波的方法是常用的多分量方法。有些地震波能量也可通过水平地震检波器被检测到，这会导致信号的消除。

在英国专利申请GB-A-2273358中，提出使用线性约束自适应聚束和自适应干扰消除聚束来抑制地面波的方法。该方法过滤被地震检波器阵列检测到的信号并将它们以某种方式相加，以维护从优先方向发出的信号，同时抑制从其他方向发出的干扰。滤波通过使用连续地自适应方法连同对用来构成主要通道和参考通道的地震反射和地面波进行扩散差分实现。建议的应用是钻井，当使用一钻头作为地震源时，由于钻头的缓慢移动，地面波是有效平稳的，且每个源接收器位置产生大量数据。这保证了用于本方法的自适应滤波器的算法的随机变化率类型是能够收敛的。然而，在表面地震实验中，地面波的表现通常是不平稳且不均匀的，算法的随机变化率类型可能非常慢而不能在信号包络内收敛。

美国专利5,237,538提出从地震数据中消除附加噪声的方法。该方法首先指出噪声的扩散特征，定义并提取了包含噪声的空间-时间门，并消除扩散以平滑噪声系列。这样，幅度和时间的变化就从门被消除。附加噪声通过聚束控制(beamsteer)算子(在此情况下是常规的叠加)被估计或通过在傅利叶变换域的f-x滤波被估计。使用短的三到五个点的单独滤波器，滤波噪声估计被从包含信号-加-噪声的数据道中减去。相反的幅度尺度被用来消除早期幅度补偿的影响。信号于是被扩散存储到原始地震记录。该方法用于地面波衰减时有一些显著的缺点。首先，特别地对于较短的阵列，信号总是泄漏到地面波估计中。实际上，总有一个信号分量出现在干扰通道中，其在时间上与主要通道中的信号共位。在另一方面，当阵列被容许更长一些时，地面波的分散出现使得难以完成有效的聚束。

根据本发明，一种在离散时间间隔([1，2，...，T])中从离散的带噪声地震信号中过滤噪声的方法，所述方法包括以下步骤

-使用所述带噪声地震信号对所述时间间隔决定至少一个参考通道(x(t))作为对所述噪声的估计；

-决定M个时间局部滤波器(w(i，t))的系数，所述滤波器构成滤波器组，M为大于等于2的数，通过最小化表示所述滤波器组的输出误差的量度的代价函数(J(t))；及

-应用所述的滤波器组到至少一个估计(x(t))以决定所述噪声的M个滤波估计。

代价函数(J(t))最好是时间的全局量，代表所述滤波器组在时间间隔([1，2，...，T])的全部输出的误差量度。

进而代价函数(J(t))可被最小化，利用当在时间域所述窗函数重叠时，对相同的信号，M个滤波器(w(i，t))的相邻滤波器的输出，被窗函数加权后的和是相等的近似。

推荐的方法包括将M个滤波估计与时间窗函数(h(i，t))相乘的步骤。

时间窗函数的特征在于要求只有相邻窗重叠。

时间窗函数，及所得时间窗对混合分量的应用，保证了滤波器在时间上是局部的，且容许该方法参照一全局最优判据从地震信号中自适应地消去噪声，例如通过最小化滤波信号在时间上的均方值来解决滤波信号的最优化问题。

本发明可用于二维(2D)和三维(3D)地震测量，可被用于陆地地震，海洋地震包括海底地震，及过渡区地震。

该方法可在存储数据或刚被采集的原始地震数据上实现。原始地震数据可根据数据采集点处的方法被过滤。这保证了从数据采集点得到“干净的”信号，且可被以这种形式从站点直接下载。这减少了必须被发送以进行离站分析的数据量及减少了与积累足够多的带噪声数据用于离站分析有关的代价和存储问题。该方法可用于单传感器记录，即用于将两个或多个地震传感器的信号混合成的任意组形式之前的记录。

噪声或干扰信号可在经过自适应滤波器之前被预处理，通过将信号分解到不同的频段，例如通过使用正交镜象滤波器。其允许要被处理的数据点在数量上的减少，也允许在自适应滤波器中的系数在数量上的减少。

数据选择时间窗函数最好由两条要求决定，其一是在任何给定时间所有窗的和等于一个单位，其二是只有相邻窗重叠。这些要求保证了被过滤信号的全局优化可通过使用一个近似来解决，对所有时间和所有滤波器和所有相邻滤波器的和，相邻滤波器的误差函数被与滤波器自身相关的误差函数代替。

数据选择时间窗函数的应用将被用来解决过滤信号的优化问题的方程解耦合。所得的解耦合方程可被例如主分量方法或代替地衰减的最小平方方法求解。

当对方程解耦合使用主分量方法被解时，主分量的个数可以被调整以改变滤波的程度并得到所希望的滤波准确度。信噪比越小，用于完成滤波的主分量的个数越多。

推荐的自适应滤波通过使用滤波器，或由许多局部滤波器组成的滤波器组实现。每个局部滤波器可以是多通道滤波器。在本发明进一步推荐的实施方式中，另一信号被用于配置自适应滤波。

在基于陆地地震的测量中，参照本发明的方法中用到的地震数据信号可通过至少两种地震传感器得到，典型地，地震检波器，该传感器分布于网络上，其中每个传感器都是独立的，每个传感器传送数据到一处理装置以被参照本发明滤波。网络被定义为传感器或地震检波器的分布阵列，其中每个传感器独立于其他所有传感器。这不同于现有技术的分布，其中每个地震检波器分布于内在相连的阵列上，且地震检波器被在全部阵列平均以试图在离站处理前减少噪声如地面波的影响。

在本发明的另一实施方式中，滤波器的输出信号可以被以迭代方式再次处理以进一步过滤噪声，典型地滤波器输出信号被反馈到至少一个参考信号和滤波器中。

参考信号可由一系列技术产生，例如通过计算数据信号的扩散差分得到。

在另一技术中，参考信号可通过中值叠加(median stacking)得到。它抑制了与要被过滤的噪声信号不同的扩散的地震信号。这就减少了参考信号对感兴趣地震信号的损害。

当这两种技术混合时，叠后信号的扩散差分可跟在中值叠加法之后。

当从多分量传感器得到地震数据后，参考信号可通过对传感器的每个分量极化滤波后得到。该方法可被用来增强从相邻传感器得到的参考信号，其对三维勘探非常先进，即三个分量，也可与扩散差分混合。

通过下面的详细说明及附图，本发明的这些和其他特点，推荐的实施方式及其变形，可能的应用及先进性将会被熟练的技术人员理解。

图1示出了使用聚束方法进行自适应消除的现有技术方法的处理框图。

图2A示出了根据本发明的自适应干扰消除方法。

图2B示出了图2A方法的一个变形。

图3示出了在自适应干扰消除中使用的用来支持滤波器系数的数据信号采样。

图4示出了根据本发明对参考通道采用扩散(moveout)差分的数据自适应多通道滤波器组。

图5示出了在自适应干扰消除中使用的用来支持滤波器系数的数据信号采样的替换方式。

图6示出了根据本发明基于扩散差分和中值叠加方法的参考通道的产生。

图7示出了基于中值叠加和扩散差分方法的参考通道的产生。

图8示出了基于极化差分的参考通道的产生。

图9示出了基于极化和扩散差分的参考通道的产生

现在参看附图，图1示出了在例如GB-A-2273358中说明的线性抑制自适应多通道聚束阵列。相同输入10，12从两列传感器产生，输入10从主要源或传感器发出，随后的输入12′，12″…，12^N从参考源或辅助传感器发出。输入10包括带噪声的信息负载信号。

参考输入12′，12″…，12^N被包含局部自适应滤波器14′，14″…，14^N的自适应多通道滤波器处理，处理后的分量在加法器16上被混合以产生输出信号y(t)18。

该滤波使用时间连续的自适应方法实现，即最小均方法(LMS)以优化滤波器性能。然而，随机变化率类型算法如最小均方法可能会收敛太慢，以致不能产生无噪声信号的估计，当附加噪声为非平稳和非均匀的时，例如在表面地震时遇到的情况。地面波噪声通常是弥散的，这经常表现为非平稳特性。

根据本发明的自适应干扰消除方法，当计算自适应滤波器系数时，估计了相位和幅度扰动的影响及地震信号传播的影响。这使得自适应干扰消除方法对数据中出现的相位和幅度扰动具有耐韧性。

图2A示出了根据本发明从地震信号过滤噪声的方法实现的系统。该系统包括许多地震检波器Gn，n＝1，…，N24′，24″，…，24^N，其记录了一般形式的信号{g_nj(t)：n＝1，…，N；j＝1，2 3}，其中j指多分量记录，例如3C地震检波器记录。为避免与图1混淆，应注意此时字母N用来表示传感器的个数，而字母K(见下)值表示通道的最大数目。

地震检波器24为在一网络中放置在地表的分开的传感器，用来检测相应于源波在地球内的地层处通过反射产生的地震波。地震检波器彼此互不相连，其用来提供一系列输入到主要和参考通道产生器26。使用压力传感器如水听器的等效系统，可代替地震检波器被使用。

收到所述记录信号的信号产生器26指明了主要通道d(t)28，其带有被附加干扰n₀(t)损害的信息负载信号s(t)，即

[1]      d(t)＝s(t)+n0(t)

及干扰参考信号通道x_k(t)30： $\left[2\right]----x_k\left(t\right)=F\left\{g_{\mathrm{nj}}\left(t\right)\right\}$

其中k＝1，…，K表示参考通道的数目。运算符F可以是用于记录信号并产生一个或多个参考通道的任何运算，参考通道可被看作附加干扰的第一次估计。因此，有以下假设：相对于x_k(t)与n₀(t)的互相关，信号s(t)与主要通道中的干扰n₀(t)的互相关，及与参考通道x_k(t)的互相关可忽略不计。例如，运算符F可为从已知的许多方法中选出的一种。一些用来产生参考信号30的特别先进的方法在下面被讨论。

参考通道x_k(t)被作为自适应滤波器组34的一系列局部自适应滤波器32′，32″，…32^M、加法器36’，36”，…36^M、乘法器38’，38”，…38^M和混合器40处理，以产生输出信号 $\left[3\right]----y\left(t\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{\mathrm{\τ}=0}{\overset{L-1}{\mathrm{\Σ}}}{h}_i\left(t\right)w_{\mathrm{ik\τ}}{x}_k(t-\mathrm{\τ})$

其中W_ikτ(t)是自适应滤波器的可调整的系数，h_i(t)是在输出末端用到的窗函数，M是局部多通道自适应滤波器的个数(或输出窗的个数)，L是每个通道的的系数的个数。每个自适应滤波器可以是图1所示类型的多

式[3]可被改写为使用抽头-输入向量x(t)在时间t的标量积的(加窗)和，定义为：

      x(t)≡[x₁(t)，...，x₁(t-L+1)，[4]       x₂(t)，...，x₂(t-L+1)，

      x_K(t)，...，x_K(t-L+1)]^T

抽头-加权向量w_i定义为

      w_i(t)≡[w_i10(t)，...，w_i1(L-1)(t-L+1)，[5]       w_i20(t)，...，w_i2(L-1)(t-L+1)，

      w_iK0(t)，...，w_iK(L-1)(t-L+1)]^T

滤波器组输出被与主输入混合以产生一过滤后的输出信号e(t)44，

[6]   e(t)＝d(t)-y(t)

在图2A中主输入在每个滤波器32的输出端求和。一个决定e(t)的等价方法如图2B所示，其中窗函数h_i(t)与每个滤波器的输出相乘，在所有输出被混合后加上主要通道。

信号e(t)也可作为误差信号，代价函数J被定义为等于对全部采样值输入的e²(t)的和，即 $\left[7\right]----J=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{e}^2\left(t\right)$

其中T是采样值的总个数。代价函数J对应于全部的误差能量，通过最小化代价函数，也就是最小化全局判据，e(t)就变为s(t)即无噪声纯信号的最小平方估计。

最小化全局判据的一个方法是使用窗函数h_i(t)以保证滤波器在时间上是局部的，即阵列具有有限的时间跨度，且可处理非平稳的噪声。窗口可以是100个采样值长的矩形窗。

窗42被下式约束 $\left[8\right]----\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{h}_i\left(t\right)=1$

对t＝1，2，...，T，其中T为采样值的总个数，且

                [9]      h_i(t)h_j(t)＝0

对j≠i-1，i，i+1。第一个限制保证了滤波器组34等效于单滤波器情况，如果所有局部滤波器32是相同的。第二个限制保证了窗口具有紧凑支持性。利用这两个限制能够解决误差信号e(t)的全局最优化。

优化问题是 $\left[10\right]----\underset{w_1,...,w_M}{\min }\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{e}^2\left(t\right).$

总的误差能量可被扩展为 $\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{e}^2\left(t\right)=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}[\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{h}_i\left(t\right)f_i\left(t\right){]}^2=$ $\left[11\right]----=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}[\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{h_i}^2\left(t\right){f_i}^2\left(t\right)+\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=i-1}{\overset{i+1}{\mathrm{\Σ}}}{h}_i\left(t\right)h_j\left(t\right)f_i\left(t\right)f_j\left(t\right)]$

其中f_i(t)是与局部滤波器w_i相关的局部误差函数： $\left[12\right]----f_i\left(t\right)=d\left(t\right)-w_i^T\left(t\right)x\left(t\right).$

从式[11]可见最优化问题导致了关于滤波器的大量联立方程组。使用第二个条件，式[9]加在窗上及下面对方程[11]第二项的近似[13]，方程可被分解。 $\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=i-1}{\overset{i+1}{\mathrm{\Σ}}}{h}_i\left(t\right)h_j\left(t\right)f_i\left(t\right)f_j\left(t\right)\≈$ $\left[13\right]----\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=i-1}{\overset{i+1}{\mathrm{\Σ}}}{h}_i\left(t\right)h_j\left(t\right)f_i^2\left(t\right)$

首先这是一个非常适度的近似，式[11]中近似的第二项远小于第一项，其次该近似比w_i≈w_i+l或w_i≈w_i-1更适合。式[13]的近似要求，当用于在相邻窗重叠的时间域上相同的输入数据时，相邻的滤波器产生相似的结果，而不是要求相邻滤波器在点到点的基础上相似。因此，该近似等同于要求两个函数的积分是近似的，而不是要求函数本身。

虽然优化导致了大量关于滤波器的联立方程组的产生，式[13]的近似分解了优化因此 $\left[14\right]----\underset{w_1,...,w_M}{\min }\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{e}^2\left(t\right)\⇔\underset{w_i}{\min }\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{h}_i\left(t\right)f^2\left(t\right)$

其中I＝1，2，0...M。

求滤波器组分量的矩阵标准方程可被得到，例如通过设变化率为零，或通过正交定理，对 $\left[15\right]----\left[\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{h}_i\left(t\right)x\left(t\right)x^T\left(t\right)\right]w_i=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{h}_i\left(t\right)d\left(t\right)x\left(t\right).$

这可通过不同的方法被解，其中两种方法是主分量法和衰减最小平方法。

式[15]的右边可被解释为相关矩阵Φ与滤波器系数的乘积。Φ的奇异值分解是 $\left[15.1\right]----\mathrm{\Φ}=\underset{n=1}{\overset{\mathrm{KL}}{\mathrm{\Σ}}}{\mathrm{\σ}}_n{u}_n{u}_{N^{\′}}^T$

其中σ₁＞σ₂＞...＞σ_NL是奇异值。奇异值分解给出了好的结果，即使矩阵是病态条件的。应用奇异值分解的另一优点是它允许主分量方法的使用。

在系统和控制理论文献中，术语′主分量′是指在平方矩阵的奇异值分解中的优势项，例如参见B.C.Moore，″Principal component analysis inlinear systems：controllability，observability，and model reduction″，IEEETrans.Autom.Control，Vol.AC-26，pp.17-31，1981。

对应于矩阵Φ的主分量矩阵Φ₁是 $\left[15.2\right]----{\mathrm{\Φ}}_1=\underset{n=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\σ}}_n{u}_n{u}_n^T$

其中1≤p≤NL。主分量矩阵Φ₁提供给Φ最优低秩近似值。式[15]的最优化问题可通过确定逆Φ₁ ^-1解决。

主分量方法混合时间延迟以定义滤波器支持(support)的使用允许滤波器集中于波的分量(即地面波)，其具有一定的扩散并附加在扩展的数据集范围中。进一步，只使用给定数目的主分量导致只有最具优势的附加分量会影响滤波器系数的计算。主分量的数目是可调整的参数，其允许用户改变滤波器的长度。使用的主分量的数目越多，滤波越强。通常，信噪比越小，应该被应用的滤波越强。

解矩阵方程的另一替换方法是衰减最小平方方法，使用： ${\mathrm{\Φ}}_2=\mathrm{\Φ}+{\mathrm{\ϵ}}^{2}I=\underset{n=1}{\overset{\mathrm{NL}}{\mathrm{\Σ}}}({\mathrm{\σ}}_n+{\mathrm{\ϵ}}^2)u_n{u}_n^T$ $\left[15.3\right]----{\mathrm{\Φ}}_2^{-1}={(\mathrm{\Φ}+{\mathrm{\ϵ}}^{2}I)}^{-1}=\underset{n=1}{\overset{\mathrm{NL}}{\mathrm{\Σ}}}{({\mathrm{\σ}}_n+{\mathrm{\ϵ}}^2)}^{-1}{u}_n{u}_n^T$

其中I是NL×NL单位矩阵，ε²是一个小的参数。在衰减最小平方方法中，在具有小的相关奇异值的分量被忽略的意义上，病态条件被防止，主分量结果被近似。主分量的分布只有轻微的改变。该方法的主要优点在于，由于矩阵Φ+ε²I是对称的，它可被使用如Cholesky分解解决，该方法比奇异值分解快得多。

而在如上所述的本发明的变形中，依靠主要地面波速度，最小信号保护间隔和滤波器长度，多通道滤波器支持在应用的开始被固定。该建立对水平分量作为参考的应用不理想，由于垂直和水平分量间的有效扩散在三个通道间不固定。改变滤波器长度可能会不错，为了局部地调整多通道滤波器支持，以便跟随最高相关的扩散。

下面，产生根据本发明的方法使用的干扰参考通道30的几个方法将被通过举例进行描述。

方法1-基于扩散差分的多参考通道

在地震采集中，反射能量的扩散通常与附加噪声如地面波的扩散不同。该性能可被用到自适应干扰消除中以抑制地面波。图3示出了从地震测量中得到的地震轨迹(traces)的实例。地震反射46的扩散和地面波48的扩散有不同的变化率。可以看到，地面波具有比地震反射陡峭的斜度。在主要通道50，即中央轨迹上，地震信号在地面波的某个相位被叠加。在相邻轨迹上，地面波的相同相位没有与地震信号重叠。因此，通过对相邻地震道的适度采样，可得到用于自适应干扰消除的噪声参考。采样区52、54提供作用于参考通道的滤波器系数的支持(support)，且由主要轨迹50的当前采样值定义的距离得到。多通道滤波器分量的相应支持可通过适当地延迟主要通道和参考通道输入被调整。采用所描述的产生参考通道的方法，所得到的数据自适应多通道如图4所示，其中与图2所示的已描述过的部件等效的部件用相应的参考数字显示出来。在图3和4中，

[16 ]    Δτ₁＝τ₀-τ₁

并且

[17]     Δτ₂＝τ₂-τ₀-L

其中L为与第二种采样相关的自适应滤波器使用的系数的个数。

为了图示的目的，在本应用余下部分，在要被滤波的轨迹的两边最近的轨迹被作为噪声参考。

Δτ_i的选择主要依靠两个因素，即为地震反射波形的自相关函数，和地面波的扩散。地震反射波的自相关函数依赖于地震反射的功率谱。通常，在时间上被Δτ_i定义的窗被选择大于延迟，其地震反射波的自相关函数具有可忽略的值。其次，最近轨迹的采样被选择，所有地面波的主要扩散平分了采样区域。这减少了影响，由于地面波的扩散性和它的明显速度范围。

方法2-基于扩散和时空相关差分的多参考通道

目标信号和干扰可以有不同的扩散(视在速度)和时空相关。这经常发生在表面地震采集中。地面波不只在扩散上与反射信号不同，而且有比地震反射强的多的时空相关，特别是在时间方向。这使得用来定义方法1中列出的多通道滤波器分量的支持的准则可被确定。图5示出该准则是如何用在从地震测量中得到的地震道的。地震信号56的扩散和地面波58的扩散具有不同的变化率。滤波器支持的两个域被说明为主要通道60的每个最近的地震道的两边，至于图3中的采样，采样区或滤波器支持62、64、66、68留出由t₀定义的中央采样值周围的一个时间采样值的距离。这是为了保护地震信号。

地面波通常具有比地震信号大的多的在时间尺度上的时间相关，在这种情况下，中央轨迹(即主要通道)的适当部分70，72可被作为参考。

该准则在很多应用中都很有用。当选择两个参考通道，相对于方法1中描述的直接地面波的发出方向，一个通道是随机的，另一个是反随机的，输出将包含反射波，以及趋前波(一个领先的事件)，其大约为该事件低通的一半，以及一个相似的滞后波(一个延迟的事件)。如果在所有处理步骤中都要求最小相位响应，就不希望出现趋前波。如果参考通道是通过使用相对于主要通道更接近信源的参考轨迹确定的，或通过使用相对于输出点是随机的参考轨迹的部分，则只有延迟波会出现在滤波数据中。

当干扰参考只通过使用时间延迟以定义滤波器支持而产生，而没有使用任何预处理，参考通道的反射波的出现会在处理过的主要通道上产生小的趋前波和/或滞后波。

例如，当两个相邻轨迹被用作参考时，趋前波和滞后波被产生。该影响也能在频域被观察到：当参考通道集中于与主要通道中的地面波有相同极性的地面波相位时，每个滤波器的低频增益约为-1/K，其中K是参考通道的个数。这样，在零频率处的总响应为空。

衰减(roll-off)依赖于K；对K＝4，角频率约为5Hz，对K＝2，角频率约为11-12Hz。响应在中间频率也显示出相同的增益。当参考通道集中在与主要通道中的地面波相位相反极性的地面波相位时，每个滤波器的低频增益约为+1/K。于是就没有低频衰减，通常在中间频率具有衰减。

由于具有正和负极性的响应非常互补，对同样的数据分别用不同的极性实现两次自适应干扰消除是可能的，在每个频率从两个实验中选择响应幅度最高的响应，其在绝大部分频率上大于1个单位。这样，在绝大多数频率上反卷积增益小于1个单位，这就减少了噪声而更理想地维持了地震反射波形。

方法3-基于中值叠加的单个参考通道

对一系列数字的中值滤波是通过将一矩形窗传至该序列，并将该序列中的每个点用集中于当前点的窗内点的中值替代的方法实现的。该1-D中值滤波运算可应用于在一特定方向上的2-D数据集以得到在那个方向的发出波形分量估计。

参考通道可通过使用图6所示的中值叠加运算沿主要扩散增加地面波被产生。叠加运算采用了典型的3到5个相邻轨迹被完成。如图6所示，传感器，或地震检波器，72，74，76产生参考轨迹78、80和主要轨迹82。这些轨迹用于叠加运算。与该法或传统运算法不同，中值运算并不抹掉地震反射(其与地面波不具有相同的扩散)，而是抑制它们。因此，由中值叠加运算产生的参考通道90包含了对干扰地面波的更好估计，且可能比传统平均叠加法更少地损坏感兴趣的地震信号。该方法可在重叠窗上比实现以重新调整叠加扩散。

方法4-基于中值叠加和扩散差分的多参考通道

通过先应用中值叠加114再进行时移τ₁、τ₂以产生参考通道x₁(t)、x₂(t)可进一步发展先前的实施方式。图7证明了五个轨迹91、92、93、94、96是怎样被用来产生两个参考通道的。通过中值滤波被泄漏的信号能量被进一步抑制。本实施方式在过滤从少量地震检波器或听音器发出的少量轨迹时特别有用。

方法5-基于极化差分的单参考通道

地面波和地震反射经常有不同极性，地震反射是线性极化的而地面波经常表现为椭圆性极化的。当地震波被多分量传感器或地震检波器测量的时候，极化滤波通过固有的极性被用来区别地震事件。

从多分量传感器产生的分量的极化滤波得到参考通道的方法示于图8，多分量传感器116检测反射波的三个正交分量u_x118、u_y120、u_z122。极化滤波器124过滤椭圆型极化的地面波信号，且过滤后的信号被用作自适应干扰消除的噪声参考126。初始的未过滤轨迹128被作为主要通道。自适应滤波器的使用减少了残余能量因此改进了地震波衰减。该方法使用从单个三分量(3C)地震检波器发出的数据并可被用于每个独立分量以恢复地面波衰减多分量数据。

方法6-基于极化和扩散差分的多参考通道

如图9所示，该方法使用至少是3C的参考地震检波器130，132和主要地震检波器134。主要检波器产生主要通道136，如在时间t₀处定义。极化滤波器138和140如方法4中所用，用来过滤椭圆型极化的地面波信号和产生参考轨迹142、144。参考轨迹142、144被延迟t₁和t₂以便根据方法1中描述的扩散差分方法采样参考轨迹。延迟的参考轨迹被用来提供参考通道。采用两个参考检波器实现该方法是可能的，主要通道从过滤前的一个分量得到。

在得到参考通道的该方法中，两个完全独立的判据扩散和极化被同时用来在地震信号和地面波之间求差分。该方法比方法4具有更强的衰减地面波的能力。

方法7-使用迭代自适应干扰消除产生的多参考通道

如上所述，自适应干扰消除的概念依赖于在噪声参考通道中信号分量的消失。事实上，这是不可能的，总是有一些信号的泄漏。从刚开始运行的自适应干扰消除中估计的噪声可以被用作第二个运行的噪声参考，其改进了从原始数据得到的信号和噪声的分离。采用迭代方法，利用过滤后的输出信号持续地改进了输出信号，可被使用多次且认为是先进的。对具有多参考通道的方法，该方法仍有增加有效阵列口径且只在计算复杂性上有线性增加的效果。

方法8-参量波形域分解作为参考产生器

在参量波形域分解中，分解地震数据集到其分支波形域的问题被作为参量求逆，其中每个波形域被它的傅利叶分量和频率无关参数(如明显的慢速)建模。该方法通过最小化频域中模型与数据波形间的最小平方误差解得波形域参数。由于明显慢速的参数非线性地进入问题，基于变化率的搜索技术被用来决定最优参数值。通过在每一频率解一线性系统然后进行傅利叶反变换，波形域被从数据中重构。

参量波形域分解方法可用于最初将数据分为地面波和地震反射信号分量中。由参量波形域分解提供的地面波估计随后被自适应干扰消除法作为噪声参考使用。参量波形域分解通常可完成好的地面波与信号的分离；然而它受数据中幅度和时间波动的影响。混合使用参量波形域分解和自适应干扰消除则有改进的效果。

本实施方式可进一步发展为，首先使用参量波形域分解然后时间移动来产生参考通道。在参量波形域分解中可能泄漏的信号能量被进一步抑制。参量波形域分解的实例可见C.Esmersoy，″Inversion of P and SV wavesfrom multicomponent offset vertical seismic profiles，″Geophysics，Vol.55，pp.39-50，1990。

方法9-水平分量作为参考通道

在多分量记录中，绝大多数被压缩的体波能量在垂直分量。水平分量(横向和径向的)包含地面波和剪切波能量，因此可被作为自适应干扰消除中的干扰参考。与基于假设地面波为椭圆型极化的极化滤波方法不同，这种3-C自适应干扰消除假设局部地(在时间域)，垂直与水平地面波分量之间有某种未知的联系，并有效地估计了该联系。由于地震反射信号也会出现在水平分量中，参考通道通过对水平分量进行时移而产生。

采用水平传感检波器作为表面波检测器用于自适应过滤垂直分量的一般思想参见G.A.Crews和D.R.Martinez的专利(美国_A_4,890,264)的描述。

现在描述根据本发明实现三维表面测量，即基于陆地的测量的方法。

3-D勘测的实现

在3-D表面地震堪测中，目前的实践是利用检波器的2-D地震检波器的天线阵列且对源位置的变化的水平分布收集数据。这些阵列典型地包含24个地震检波器，输出信号通过硬连接检波器及对它们的输出求和得到。

本文描述的自适应干扰消除可用于3-D勘测中，其具有地震检波器阵列模式，或网络，例如带中央检波器的六边形分布。如果直接地面波的发出方向(或，等效地，源位置)没有在处理时刻被提供，它可通过互相关任意参考轨迹与其他轨迹被估计。每个轨迹被依次作为主要通道，并选择一系列相应的参考轨迹，这样在主要轨迹和参考轨迹之间就存在足够的地面波扩散。通过使用根据本发明的方法得到的所有过滤过的轨迹被以常规聚束(平均运算)或中值叠加方法叠加。每个过滤后的轨迹的延迟为不同的1和2值保证了信号的趋前和滞后特征通常是不被强调的，而地震反射被强调。

这些检波器网的表面维数远小于目前使用的连接天线阵列。另外，过滤后的轨迹的随后叠加比目前的实践更多地维持了信号带宽，特别对于大的源到接收器的偏置。

使用完备重构滤波器组的实现

在绝大多数地震勘探中，附加噪声只占可得时间带宽的小部分。例如，在本文使用的测试数据，奈奎斯特频率为250Hz，而绝大部分地面波能量在30Hz以下。为减轻计算负担，集中滤波作用到附加噪声所在频段是希望的。达到此目标的一种方法是增加QMF(正交镜象滤波器)完备重构滤波器组到使用自适应多道滤波器组的噪声抑制系统。因此，两个滤波器组被有效地包含在系统中。QMF滤波器组在使用自适应滤波前将轨迹分解到频段并抽取，然后用于重新合成。使用QMF滤波器组抽取减少了要处理的点数并允许自适应滤波器系数的减少，明显节省了计算机处理时间和对计算机内存的要求。

由于在QMF滤波器组设计中需要频率重叠，在频带间存在一些泄漏。该泄漏的影响可通过展宽处理频段被解决。对残余地面波的附加抑制可通过构成被视为简单K-域滤波器或聚束器的常规组被完成。根据本发明的本方法避免了如地面波之类的混叠噪声出现的数据处理问题。

当地震检波器(或源)的采样间隔大于奈奎斯特定理的要求时出现混叠的地面波。在有些采集几何中，这可通过设计减少在该域采用的地震检波器(与/或源)的个数来实现。

要求合适的采样的方法，如k或f-k域滤波器不能有效地分离信号和地面波，如频域(f或f-k)滤波器在变换域为信号和噪声分配明确的区域来将地震信号和地面波一步分离那样。混叠的出现阻止了这样做。根据本发明的方法不变换数据，所以它不被混叠影响。它集中于噪声并消除它；剩余的信号为过滤后的信号。

现在讨论根据本发明的方法用于真实数据集合时与传统方法的比较。所用数据集合从一般的接收器集中得到，其中连续射点的位置相距4米。为了能考察对叠前信号和噪声两者的处理效果，增加了两个合成反射。对处理的进一步评价在给数据集增加任意尖峰噪声后被实现。

传统的分组技术的性能被相邻道在滑动窗内的算术平均仿真。对一传统的六个检波器的组(20米组)，某些地面波抑制可被完成，但增加了地震事件的失真。

当在自适应干扰消除中使用已知的归一化算法分析数据时，归一化的LMS算法无法跟踪地面波。采用LMS算法自适应干扰消除抑制地面波的效果有限。LMS算法不能足够快地收敛以跟上地面波的瞬变和非平稳性质。

根据本发明的方法用于测试数据集，多通道滤波器组的使用很好地维护了地震反射，地面波被有效地抑制。在本例设k＝2，即对每个主要通道有两个参考通道，L＝30，即每个道有30个系数，输出窗是长度为100的矩形窗，Δτ₁＝Δτ₂＝3个采样值；p＝5是使用的主分量个数。

根据本发明的方法除用于地面波外还可用于不相干的噪声源，如海洋地震勘探或海底勘探中的轮船噪声和钻井噪声。

在海洋环境中的实现

在海洋地震勘探中，声音源产生从水中传到地球内的波。这些波随后被表面下的地层结构反射或折射，返回水中并通过常的水听器阵列被记录，水听器阵列被拖到地震船后水面的附近。水听器被固定在拖缆上，通常有2-12个拖缆被拖着，每个拖缆通常有几公里长。拖缆由100米长的部分组成；每部分由里面的水听器和浸油外壳组成。拉力绳和垫圈构成了拖缆的骨架。

当拖缆被拖到船后时，由于源的变化产生了自-噪声。船的摇晃，特别在浪大的海上，导致拉力绳的振荡，拉力绳与连接器和油浸外壳相互作用，产生沿着拖缆传播的凸起波(或呼吸波)。压力变化被水听器检测到，增加并损坏了检测到的地震信号。随着拖缆在水中的移动，边界层扰动导致在外层壁的压力波动，其也传到水听器中。

根据本发明的自适应干扰消除方法可直接用于凸起波。先前的针对地震勘测的自适应干扰消除不得不使用特殊的传感器(传感器或过载传感器)来形成噪声参考通道，参见美国4821241，美国5251183。根据本发明的方法只通过水听器数据来形成噪声参考通道。使用3分量数据或极化滤波的方法不能用于拖缆噪声问题，但其他方法可用。

在有些海洋地震勘测中，存在从交叉线方向发出到拖缆的干扰。出现这种干扰的一种情况为由油井的钻探活动产生。由于可以有至少两个拖缆在船后，在干扰方向构成的水听器阵列将没有足够的水听器阵列以使用常规的聚束，或数据无关多通道滤波，例如f-k滤波。另一方面，根据本发明的自适应干扰消除方法可直接用于此问题。沿着交叉线的阵列，干扰将有更大的扩散，和明显地慢于地震信号，在应用自适应干扰消除方法时这一特点被利用。

根据本发明的方法的使用允许用来采集地震数据的地震检波器或水听器在数量上有明显减少。

Claims (19)

1.在时间间隔([1，2，...，T])中从离散的带噪声地震信号中过滤噪声的方法，所述方法包括以下步骤

-使用所述带噪声地震信号对所述时间间隔决定至少一个参考通道(x(t))作为对所述噪声的估计；

-决定M个时间局部滤波器(w(i，t))的系数，所述滤波器构成滤波器组，M为大于或等于2的数，通过最小化表示所述滤波器组的输出误差的量度的代价函数(J(t))；及

-对所述至少一个估计(x(t))应用所述的滤波器组，以决定所述噪声的M个滤波估计。

2.根据权利要求1的方法，其中代价函数(J(t))是时间的全局函数，代表所述滤波器组在时间间隔([1，2，...，T])的总输出的误差量度。

3.根据权利要求1的方法，包括将M个滤波估计与时间窗函数(h(i，t))相乘的步骤。

4.根据权利要求3的方法，其中时间窗函数的特征在于要求只有相邻窗重叠。

5.权利要求3的方法，其中，利用当在所述窗函数重叠的时间域应用于相同信号时M个滤波器(w(i，t))的相邻滤波器的输出由窗函数加权后的和相等这一近似，将代价函数(J(t))最小化。

6.权利要求1的方法，其中作为输入的离散带噪声信号是来自在任何组形式技术前的独立地震传感器中的记录。

7.根据权利要求1到5中任一项的方法，其中过滤输出信号的全局最优化是通过主分量方法解决的。

8.根据权利要求1到4中任一项的方法，其中过滤输出信号的全局最优化是通过衰减最小平方方法解决的。

9.根据权利要求7的方法，其中主分量的个数被调整以改变滤波程度。

10.根据权利要求1到6之一的方法，其中通过使用包括一系列局部滤波器的滤波器装置进行自适应滤波。

11.根据权利要求10的方法，其中每局部滤波器是多通道滤波器。

12.根据权利要求1到6之一的方法，其中过滤输出信号被用迭代的方式再处理以进一步过滤噪声。

13.根据权利要求1到6之一的方法，其中参考信号是通过数据信号的扩散差分得到的。

14.根据权利要求1到6中任一项的方法，其中参考信号是通过中值叠加方法得到的。

15.根据权利要求14的方法，其中中值叠加后面是叠加信号的扩散差分。

16.根据权利要求1到6中任一项的方法，其中地震数据从多分量传感器得到，参考信号通过对传感器得到的每个分量进行极化滤波得到。

17.根据权利要求15的方法，其中极化滤波与扩散差分混合。

18.根据权利要求1的方法的迭代的使用，其中在该滤波方法的前一个应用中确定的滤波噪声估计被用作随后的滤波步骤中的噪声估计。

19.权利要求1的方法在陆地、海洋或过渡区的2D或3D地震勘测中的使用。

Images