CN105277986A - 基于自适应匹配滤波算子的可控震源谐波压制方法 - Google Patents

Abstract

本发明提供一种基于自适应匹配滤波算子的可控震源谐波压制方法，该基于自适应匹配滤波算子的可控震源谐波压制方法包括：步骤1，用传统预测滤波法对谐波干扰进行预测；步骤2，用伪多道匹配方法对步骤1中预测的谐波进行自适应修正，结合传统的多道匹配和约束道均衡方法，得到扩展的伪多道匹配滤波方法；步骤3，从原始记录中减去步骤2中经过修正后的预测出的谐波，得到谐波压制后的记录。该基于自适应匹配滤波算子的可控震源谐波压制方法效果好，精度高，能够很好地达到可控震源数据处理的要求的，并且为后续地震资料成像和解释提供可靠保证。

Description

基于自适应匹配滤波算子的可控震源谐波压制方法

技术领域

本发明涉及地震资料处理领域，特别是涉及到一种基于自适应匹配滤波算子的可控震源谐波压制方法。

背景技术

随着可控震源地震勘探的出现，人们很早就认识了谐波，由于滑动扫描出现之前的扫描方式中，谐波干扰对数据品质的影响很小，几乎可以忽略，所以谐波压制技术发展缓慢。为了提高生产效率，出现了滑动扫描等高效采集技术，由于人们进一步追究采集速度，滑动时间越来越小，谐波干扰影响也就越来越大，因此压制谐波干扰的技术不断更新，特别是国外，近些年提出了一系列压制谐波干扰的方法。如：野外采集中普遍使用的变频变相位扫描技术，该技术是基于扫描信号和谐波畸变信号的极性特征，在采集过程中通过改变扫描信号的初始振动相位，之后采用先相关后叠加的方式可以对谐波干扰进行有效的压制，但是该方法针对谐波出现的具体位置来选择合适的参数压制谐波，实现起来比较困难，质量监控也很难进行。Li等人(1994)基于线性频率扫描方式，提出了纯相移滤波(PPSM)方法来压制谐波干扰，该方法简单高效直观稳定，特别适用于VSP资料，但不适用于滑动扫描记录。对此，黄建平等人(2012)对其进行改进，提出了一种可以适用于滑动扫描技术的相移滤波法，该方法能够有效的压制本炮干扰，但是对于邻炮干扰的压制效果不是很好。

在可控震源地震勘探过程中，可控震源在向地下传输扫描信号的同时，由于种种因素，谐波干扰不可避免。基于上述原因在可控震源进行地震勘探时，必然存在谐波干扰。而谐波的存在对可控震源地震数据后续的处理、成像和解释带来了巨大的问题。

发明内容

本发明的目的是提供一种效果好，精度高，能够很好地达到可控震源数据处理的要求的基于自适应匹配滤波算子的可控震源谐波压制方法。

本发明的目的是通过如下技术措施来实现的：基于自适应匹配滤波算子的可控震源谐波压制方法，该基于自适应匹配滤波算子的可控震源谐波压制方法包括：步骤1，用传统预测滤波法对谐波干扰进行预测；步骤2，用伪多道匹配方法对步骤1中预测的谐波进行自适应修正，结合传统的多道匹配和约束道均衡方法，得到扩展的伪多道匹配滤波方法；步骤3，从原始记录中减去步骤2中经过修正后的预测出的谐波，得到谐波压制后的记录。

本发明的目的还可通过如下技术措施来实现的：

步骤1包括：

A)采用相移法分离地面力信号求取ρ_m，令

$P=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1}$

其中P为预测滤波算子，M为谐波的阶次，H_m为m次谐波的相位谱，H₁为基波的相位谱，ρ_m为m次谐波振幅与基波振幅的比值；

B)对谐波干扰进行预测和压制:

$H=R{a}_1{H}_1\·P=R{a}_1{H}_1\·\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1}$

$D=R\·(a_1{H}_1+\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{a}_m{H}_m)$

$D-H=R{a}_1{H}_1\·(1+\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1})-R{a}_1{H}_1\·\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1}=R{a}_1{H}_1$

其中，H是谐波干扰，D为频率域的炮记录，H_m是m次谐波的相位谱，H₁为基波的相位谱，R为地震震动向大地施加的外力，a₁为基波的振幅。

步骤2包括以下步骤：

a)假设初始地震记录为p，有效波记录为p₀，真实的谐波干扰为m，则

p₀＝p-m

假设预测的谐波分量为m₀，用a表示滤波因子，则真实的谐波干扰可由下式求得：

m＝a*m₀

则有效波记录可表示为：

p₀＝p-a*m₀

计算滤波因子a：

e(a)＝||p-Ma||²

其中M为预测谐波噪音的矩阵表示形式，e(a)为a的最小化目标函数；

b)在最小二乘意义下通过计算基于L2模的谐波自适应相减的滤波因子a来最小化预测的谐波信号与实际谐波信号的差，在一个可滑动的窗口中匹配预测的谐波干扰，使用多道方程式，

$p_0\left(t\right)=p\left(t\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{a}_i\left(t\right)*m_i\left(t\right)$

其中，p(t)是原始地震记录，m_i(t)是预测的第i道谐波干扰，a_i(t)是多道匹配因子，N是匹配的道数，p₀(t)为预测滤波后的地震记录；

根据约束道均衡的方法，将单道扩展为四道，其表达式为：

$p_0\left(t\right)=p\left(t\right)-[w_{1}m\left(t\right)+w_2\stackrel{.}{m}\left(t\right)+w_3{m}^H\left(t\right)+w_4{\stackrel{.}{m}}^H\left(t\right)]$

其中：p₀(t)为预测滤波后的地震记录，m^H(t)表示谐波的希尔伯特变换道，和为m(t)和m^H(t)导数道，w_i表示每一道的滤波算子；

采用褶积算子来调整上式：

$\begin{array}{c}{p}_0\left(t\right)=p\left(t\right)-[f_1\left(t\right)m\left(t\right)+f_2\left(t\right)\stackrel{.}{m}\left(t\right)+f_3\left(t\right)m^H\left(t\right)\\ +f_4\left(t\right){\stackrel{.}{m}}^H\left(t\right)]\end{array}$

其中：p₀(t)为预测滤波后的地震记录，和为_m(_t)和m^H(t)导数道，f_n(t)为各道频率表达式，m^H(t)表示谐波的希尔伯特变换道；

结合传统的多道匹配和约束道均衡方法，得到扩展的伪多道匹配滤波方法：

$\begin{array}{c}{p}_0\left(t\right)=p\left(t\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[f_{1,i}\left(t\right)m\left(t\right)+f_{2,i}\left(t\right)\stackrel{.}{m}\left(t\right)+f_{3,i}\left(t\right)m^H\left(t\right)\\ +f_{4,i}\left(t\right){\stackrel{.}{m}}^H\left(t\right)]\end{array}$

其中：p₀(t)为预测滤波后的地震记录，和为m(t)和m^H(t)导数道，f_n,i(t)为伪多道道频率表达式，m^H(t)表示谐波的希尔伯特变换道。

本发明中的基于自适应匹配滤波算子的可控震源谐波压制方法，针对可控震源谐波干扰，提出了自适应匹配预测滤波谐波压制技术。采用该方法压制谐波的地震数据残余谐波能量很少，并且不会减少有效波的能量。该方法的效果好，精度高，能够很好地达到可控震源数据处理的要求，为后续反演成像研究提供良好的技术支持。

附图说明

图1为谐波畸变信号的模拟；

图2为本发明的基于自适应匹配滤波算子的可控震源谐波压制方法的一具体实施例的流程图；

图3为本发明模拟试算的某工区速度模型；

图4为可控震源滑动扫描正演模拟多炮记录相关前的原始记录；

图5为可控震源滑动扫描正演模拟多炮记录相关后的相关记录；

图6为采用本发明方法进行谐波压制后的炮记录；

图7为采用本发明方法压制掉的谐波；

图8为传统预测滤波算法压制谐波的结果；

图9为采用预测滤波法压制掉的谐波；

图10为采用本发明提出的谐波压制方法得到的第二炮炮谐波压制结果；

图11为采用本发明提出的谐波压制方法得到的第三炮炮谐波压制结果；

图12为采用本发明提出的谐波压制方法得到的第四炮炮谐波压制结果。

具体实施方式

为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合附图所示，作详细说明如下。

本文中假设使用的是正弦扫描信号，扫描信号的表达式如式(1)所示：

式中，a₁为振幅随时间的表达式，ψ₁为初始相位，T为扫描周期，为频率表达式，可由瞬时频率积分得到

我们给出k次谐波的表达式如式(2)所示：

谐波畸变信号(如图1所示)由基波信号和各次谐波信号组成，因此我们将谐波畸变信号由下式表示：

如图2所示，图2为本发明的基于自适应匹配滤波算子的可控震源谐波压制方法的流程图。

在步骤101，用传统预测滤波法对谐波干扰进行预测。谐波畸变信号由基波信号和各次谐波信号组成，因此我们将谐波畸变信号由下式表示：

式中，a_m为振幅随时间的表达式，ψ_m为初始相位，T为扫描周期，t为时间，为频率表达式，令式(3)可以简化为式(4)所示的基波信号和各次谐波信号的和的形式。

$s\left(t\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_m{h}_m\left(t\right)=a_1\left(t\right)h_1\left(t\right)+\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{a}_m{h}_m\left(t\right)---\left(4\right)$

将该谐波畸变信号与反射系数做褶积运算，就可以得到含有谐波干扰：

$d\left(t\right)=r\left(t\right)*(a_1\left(t\right)h_1\left(t\right)+\underset{m=2}{\mathrm{\Σ}}{d}_m\left(t\right)h_m\left(t\right))---\left(5\right)$

其中：d(t)为含谐波干扰的炮记录，a₁为振幅随时间的表达式，t为时间，d_m(t)为谐波分量产生的炮记录表达式，‘＊’为褶积运算符号，r(t)为反射系数，为地面力信号。将式(5)变换到频率域得：

$D=R\·(a_1{H}_1+\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{a}_m{H}_m)---\left(6\right)$

其中：D为频率域的炮记录，H_m是m次谐波的相位谱，R为地震震动向大地施加的外力，a_m为傅立叶变换后振幅谱的权系数，m为谐波的阶次。将式(6)变形得

$D=R{a}_1{H}_1\·(1+\underset{m=2}{\overset{M}{\mathrm{\Σ}}}\frac{a_m{H}_m}{a_1{H}_1})=R{a}_1{H}_1\·(1+\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1})---\left(7\right)$

其中：D为频率域的炮记录，R为地震震动向大地施加的外力，H_m是m次谐波的相位谱，a_m为傅立叶变换后振幅谱的权系数，m为谐波的阶次，ρ_m是m次谐波振幅谱与基波信号振幅谱的比值。本发明采用相移法分离地面力信号求取ρ_m，令

$P=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1}---\left(8\right)$

其中P为预测滤波算子，则谐波干扰的预测和压制由式(9)(10)表示:

$H=R{a}_1{H}_1\·P=R{a}_1{H}_1\·\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1}---\left(9\right)$

$D-H=R{a}_1{H}_1\·(1+\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1})-R{a}_1{H}_1\·\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1}=R{a}_1{H}_1---\left(10\right)$

在步骤102，用伪多道匹配方法对预测的谐波进行自适应修正，结合传统的多道匹配和约束道均衡方法，得到扩展的伪多道匹配滤波方法。假设初始地震记录为p，有效波记录为p₀，真实的谐波干扰为m，则

p₀＝p-m(11)

假设预测的谐波分量为m₀，用a表示滤波因子，则真实的谐波干扰可由下式求得

m＝a*m₀(12)

则有效波记录可表示为

p₀＝p-a*m₀(13)

通过采用最小二乘的思想最小化目标函数计算滤波因子a：

e(a)＝||p-Ma||²(14)

其中：e(a)是求解a的最优化问题的函数，其中M为预测谐波噪音的矩阵表示形式。

在最小二乘意义下通过计算基于L2模的谐波自适应相减的滤波因子a来最小化预测的谐波信号与实际谐波信号的差。这需要有效波能量与谐波能量满足正交性，如果不满足正交性，滤波因子的求取就会出现问题，导致谐波压制不完全或者有效能量的损失。

在通常情况下，我们在一个可滑动的窗口中匹配预测的谐波干扰。使用多道方程代替式(13)。

$p_0\left(t\right)=p\left(t\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{a}_i\left(t\right)*m_i\left(t\right)---\left(15\right)$

其中，p(t)是原始地震记录，m_i(t)是预测的第i道谐波干扰，a_i(t)是多道匹配因子，N是匹配的道数，p₀(t)为预测滤波后的地震记录。如果我们仅仅使用最小二乘原理而不加任何约束，匹配的结果将会不准确。

我们结合传统的多道匹配和约束道均衡方法，得到扩展的伪多道匹配滤波方法：

$\begin{array}{c}{p}_0\left(t\right)=p\left(t\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[f_{1,i}\left(t\right)m\left(t\right)+f_{2,i}\left(t\right)\stackrel{.}{m}\left(t\right)+f_{3,i}\left(t\right)m^H\left(t\right)]\\ +f_{4,i}\left(t\right){\stackrel{.}{m}}^H\left(t\right)]\end{array}---\left(16\right)$

其中：p₀(t)为预测滤波后的地震记录，和为m(t)和m^H(t)导数道，f_n,i(t)为伪多道道频率表达式，m^H(t)表示谐波的希尔伯特变换道。

在步骤103，从原始记录中减去预测出的谐波，得到谐波压制后的记录。我们对如图3所示的国内速度模型进行试算。该速度模型存在较多的起伏地层，同时在中深层存在一个椭圆形高速异常体。模型大小为5000m×2500m，采用的网格间距为5m×5m。对可控震源的滑动扫描方式进行模拟，正演模拟的参数如下：采用线性升频扫描方式，起始频率f1＝5hz，终止频率f2＝80hz，扫描周期T＝10s。滑动时间和记录时间都是3s，时间采样间隔为0.5ms。四个震源车分别分别在1km，2km，3km，4km处激发。正演模拟得到的可控震源滑动扫描正演模拟多炮记录如图4所示，图4为相关前的原始记录，图中5为图4所示的原始记录与扫描信号互相关得到的相关记录。从图5中，我们可以清楚的看到，谐波干扰产生在负时间轴上，与前面炮记录的有效信号混叠在一起，对前面炮记录的有效信号造成干扰，如图中黑色箭头所示。

为了还原较为准确的野外采集情况，我们对真实地面力信号的振幅、相位和频率随机加一定范围内的误差，并分别采用本文提出的自适应匹配预测滤波方法和传统的预测滤波算法进行试算，得到的第一炮的谐波干扰的压制结果如图6-9所示。其中图6为谐波压制后的炮记录，图7为压制的谐波。图6中结果中，谐波干扰被压制的非常干净，而从图7中也可以看出压制的谐波干扰中不含有有效波，换言之，通过本发明的谐波压制方法，基本没有造成有效能量的损失，图8为传统预测滤波算法压制谐波的结果，谐波干扰没有压制干净，而且造成了有效信号的损失。而图9所示的传统预测滤波算法压制谐波的结果中，谐波干扰没有压制干净，而且造成了有效信号的损失。图10-12为采用本文提出的自适应匹配预测滤波方法得到的其他炮的谐波压制结果，谐波干扰都得到了明显的消除，并没有造成有效能量的损失。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，应此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (3)

1.一种基于自适应匹配滤波算子的可控震源谐波压制方法，其特征在于，该基于自适应匹配滤波算子的可控震源谐波压制方法包括：

步骤1，用传统预测滤波法对谐波干扰进行预测；

步骤2，用伪多道匹配方法对步骤1中预测的谐波进行自适应修正，结合传统的多道匹配和约束道均衡方法，得到扩展的伪多道匹配滤波方法；

步骤3，从原始记录中减去步骤2中经过修正后的预测出的谐波，得到谐波压制后的记录。

2.根据权利要求1所述的基于自适应匹配滤波算子的可控震源谐波压制方法，其特征在于，步骤1包括：

A)采用相移法分离地面力信号求取ρ_m，令

$P=\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1}$

其中P为预测滤波算子，M为谐波的阶次，H_m为m次谐波的相位谱，H₁为基波的相位谱，ρ_m为m次谐波振幅与基波振幅的比值；

B)对谐波干扰进行预测和压制:

$H={\mathrm{Ra}}_1{H}_1\·P={\mathrm{Ra}}_1{H}_1\·\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1}$

$D=R\·(a_1{H}_1+\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{a}_m{H}_m)$

$D-H={\mathrm{Ra}}_1{H}_1\·(1+\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1})-{\mathrm{Ra}}_1{H}_1\·\underset{m=2}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ρ}}_m\frac{H_m}{H_1}={\mathrm{Ra}}_1{H}_1$

其中，H是谐波干扰，D为频率域的炮记录，H_m是m次谐波的相位谱，H₁为基波的相位谱，R为地震震动向大地施加的外力，a₁为基波的振幅。

3.根据权利要求1所述的基于自适应匹配滤波算子的可控震源谐波压制方法，其特征在于，步骤2包括以下步骤：

a)假设初始地震记录为p，有效波记录为p₀，真实的谐波干扰为m，则

p₀＝p-m

假设预测的谐波分量为m₀，用a表示滤波因子，则真实的谐波干扰可由下式求得：

m＝a*m₀

则有效波记录可表示为：

p₀＝p-a*m₀

计算滤波因子a：

e(a)＝||p-Ma||²

其中M为预测谐波噪音的矩阵表示形式，e(a)为a的最小化目标函数；

b)在最小二乘意义下通过计算基于L2模的谐波自适应相减的滤波因子a来最小化预测的谐波信号与实际谐波信号的差，在一个可滑动的窗口中匹配预测的谐波干扰，使用多道方程式，

$p_0\left(t\right)=p\left(t\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{a}_i\left(t\right)*m_i\left(t\right)$

其中，p(t)是原始地震记录，m_i(t)是预测的第i道谐波干扰，a_i(t)是多道匹配因子，N是匹配的道数，p₀(t)为预测滤波后的地震记录；

根据约束道均衡的方法，将单道扩展为四道，其表达式为：

$p_0\left(t\right)=p\left(t\right)-[w_{1}m\left(t\right)+w_2\stackrel{\·}{m}\left(t\right)+w_3{m}^H\left(t\right)+w_4{\stackrel{\·}{m}}^H\left(t\right)]$

其中：p₀(t)为预测滤波后的地震记录，m^H(t)表示谐波的希尔伯特变换道，和为m(t)和m^H(t)导数道，w_i表示每一道的滤波算子；

采用褶积算子来调整上式：

$p_0\left(t\right)=p\left(t\right)-[f_1\left(t\right)m\left(t\right)+f_2\left(t\right)\stackrel{\·}{m}\left(t\right)+f_3\left(t\right)m^H\left(t\right)$

$+f_4\left(t\right){\stackrel{\·}{m}}^H\left(t\right)]$

其中：p₀(t)为预测滤波后的地震记录，和为m(t)和m^H(t)导数道，f_n(t)为各道频率表达式，m^H(t)表示谐波的希尔伯特变换道；

结合传统的多道匹配和约束道均衡方法，得到扩展的伪多道匹配滤波方法：

$p_0\left(t\right)=p\left(t\right)-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[f_{1,i}\left(t\right)m\left(t\right)+f_{2,i}\left(t\right)\stackrel{\·}{m}\left(t\right)+f_{3,i}\left(t\right)m^H\left(t\right)$

$+f_{4,i}\left(t\right){\stackrel{\·}{m}}^H\left(t\right)]$

其中：p₀(t)为预测滤波后的地震记录，和为m(t)和m^H(t)导数道，f_n,i(t)为伪多道道频率表达式，m^H(t)表示谐波的希尔伯特变换道。