CN103901473A - 一种基于非高斯性最大化的双检信号上下行波场分离方法 - Google Patents

Abstract

本发明公开了一种基于非高斯性最大化的双检信号上下行波场分离方法，涉及地震勘探技术中地震信号处理领域，本发明的目的在于提供一种利用地震信号的超高斯分布特性，在时间-空间域中逐道实现双检信号的上下行波场分离的方法。本发明利用L1范数表征地震信号的非高斯性，采用非高斯性最大化准则实现上下行波场分离。克服了采用传统的频率-波数域方法需要假设检波器是均匀分布，且要求记录面是平面的前提，本方法不需要采用上述假设，更好的适应了实际数据的处理，满足了实际应用中的多种需要。

Description

一种基于非高斯性最大化的双检信号上下行波场分离方法

技术领域

本发明涉及地震勘探技术中地震信号处理领域，尤其涉及一种基于非高斯性最大化的双检信号上下行波场分离方法。

背景技术

在海上地震勘探中，由于海水和空气接触面是一个强反射面，导致采集到的地震数据中存在激发点和接收点鬼波。鬼波的存在，会抑制低频分量，并造成频率域的陷波点，从而降低地震记录的分辨率，使地震记录的同相轴发生混叠甚至产生假的同相轴。为了获得准确的地震记录剖面，就要尽量压制鬼波。

Van Melle和Weatherborn^[1]最早在物理本质上对鬼波进行了详细的阐述，并能够从地震记录中有效的识别出鬼波，但是没有确定出相应的鬼波参数来设计鬼波的逆滤波器。Lindsey^[2]等人认为通过单炮点和单缆采集到的地震数据很难设计出有效的滤波器来压制鬼波效应，对于陆地上的地震记录可以通过合并不同炮点深度采集到的地震记录来压制鬼波。Hamarbitan和Margrave^[3]又一次通过实验说明单炮单缆采集到的地震数据很难有效地消除鬼波，而通过合并不同炮点深度采集到的数据能够更好的压制鬼波。

自1953年Melle和Weatherborn提出鬼波的概念以来，鬼波的压制一直是一个永恒的主题。经过了多年的研究和试验，鬼波的压制技术取得了很大的进展，特别是改进采集方式，利用上下缆和双检技术，可以更有效地压制鬼波。

相对于传统的单检波器只采集压力波场，双检波器可以在同一个检波点上使用陆检和水检两种检波器，同时获取速度和压力数据，这两种数据都包含上行波和下行波^[4]（其中，经过地下反射的地震波向上传播被拖缆直接接收，称之为上行波；而对于那些从地下反射回来，继续传播到达海面，并经过海面反射后向下传播才被拖缆接收到的波，被称之为下行波，也就是鬼波）。由于下行波被陆检和水检这两种检波器接收时会产生极性相反的响应，因此通过对接收到的水检和陆检数据进行上下行波场分离，达到压制鬼波的目的。

在共炮集中，双检信号可以表示如下^[5]：

$u(t,x,y)=\frac{1}{2}(p(t,x,y)-g(x,y)v(t,x,y))---\left(1a\right)$

$d(t,x,y)=\frac{1}{2}(p(t,x,y)+g(x,y)v(t,x,y))---\left(1b\right)$

其中：t是时间，x,y是检波器坐标，p(t,x,y),v(t,x,y)分别是水检和陆检信号，u(t,x,y),d(t,x,y)是待求的上行和下行信号，g(x,y)是一个随角度变化的空变尺度函数。

现有技术中采用频率-波数域分离方法（即假设记录面是平的，且检波器是均匀分布的^[7]）实现双检信号的上下行波场分离，即将地震数据进行三维傅立叶变换到f-k_x-k_y域中，公式(1a)和(1b)改写为：

$U(\mathrm{\ω},k_x,k_y)=\frac{1}{2}(P(\mathrm{\ω},k_x,k_y)-G(\mathrm{\ω},k_x,k_y)V(\mathrm{\ω},k_x,k_y))---\left(2a\right)$

$D(\mathrm{\ω},k_x,k_y)=\frac{1}{2}(P(\mathrm{\ω},k_x,k_y)+G(\mathrm{\ω},k_x,k_y)V(\mathrm{\ω},k_x,k_y))---\left(2b\right)$

式中：G(ω,k_x,k_y)对应空变尺度函数g(x,y)^[6]，定义如下：

$G(\mathrm{\ω},k_x,k_y)=\frac{\mathrm{\ρ\ω}}{k_z}---\left(3\right)$

式中：

$k_z=\sqrt{{\left(\frac{\mathrm{\ω}}{v}\right)}^2-{k_x}^2-{k_y}^2}---\left(4\right)$

其中：k_x,k_y,k_z是波数向量，ω是频率，ρ,v分别是水的密度和声波传播速度。

发明人在实现本发明过程中，发现现有技术中至少存在以下缺点和不足：

现有的双检信号波场分离方法一般采用公式（2～4），在频率-波数域中实现双检信号的上下行波场分离。但实际海上地震数据采集时，上述假设一般难以满足，不适合对现场实际数据的处理，无法满足实际应用中的需要。

发明内容

本发明提供了一种基于非高斯性最大化的双检信号上下行波场分离方法，本发明通过利用地震信号的超高斯分布特性，在时间-空间域中实现上下行波场分离，详见下文描述：

一种基于非高斯性最大化的双检信号上下行波场分离方法，所述方法包括以下步骤：

（1）建立基于非高斯性最大化的目标函数；

（2）对每一地震道输入地震双检信号，获取水检和陆检信号，分别构造出向量X(t)和Z(t)；

（3）根据非高斯性最大化的目标函数和迭代加权最小二乘方法更新加权对角矩阵；

（4）判断迭代次数是否达到最大迭代次数N，如果否，返回步骤（3）；如果是，判断所有地震道是否处理完毕，如果没有处理完，返回步骤（2），处理下一地震道，直至所有地震道处理完毕，执行步骤（5）；

（5）输出分离结果。

所述建立基于非高斯性最大化的目标函数的步骤具体为：

目标函数为：

$\underset{g}{\arg \min {\left|\right|Y\left(t\right)\left|\right|}_1}$

其中，Y(t)＝X(t)-gZ(t)；

式中：令 $Y\left(t\right)=\left[\begin{array}{c}u\left(t\right)\\ d\left(t\right)\end{array}\right],X\left(t\right)=\left[\begin{array}{c}p\left(t\right)\\ p\left(t\right)\end{array}\right],Z\left(t\right)=\left[\begin{array}{c}-v\left(t\right)\\ v\left(t\right)\end{array}\right],$ g为一未知变量，t是时间，p(t)和v(t)分别是观测到的水检和陆检信号，u(t)和d(t)分别是待求的上行和下行信号。

所述根据非高斯性最大化的目标函数和迭代加权最小二乘方法更新加权对角矩阵W的步骤具体为：

1）每一步迭代中，根据公式 $g=\frac{{\left[X\left(t\right)\right]}^T{W}^T\mathrm{WZ}\left(t\right)}{{\left[Z\left(t\right)\right]}^T{W}^T\mathrm{WZ}\left(t\right)}$ 得到变量g；

2）获取上下行波场的估计；

3）根据公式更新加权对角矩阵，ε是阈值系数。

本发明提供的技术方案的有益效果是：传统的频率-波数域方法需要假设检波器是均匀分布，且要求记录面是平面的，本方法不需要采用上述假设，通过利用地震信号的超高斯分布特性，在时间-空间域中实现上下行波场分离，更适应实际数据的处理，满足了实际应用中的需要。

附图说明

图1为基于非高斯性最大化的双检信号上下行波场分离方法的流程图；

图2为本方法的计算机程序流程图；

图3为模拟出的上行波的示意图；

图4为模拟出的下行波的示意图；

图5为模拟出的水检信号的示意图；

图6为模拟出的陆检信号的示意图；

图7为本方法分离出的上行波的示意图；

图8为本方法分离出的下行波的示意图；

图9为空变尺度函数g(x)的示意图：其中，真实的用单实线表示；通过L2范数得到的初始估计用双实线表示；采用本方法迭代10次得到的最终估计用圆圈表示。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步地详细描述。

本发明实施例利用地震信号的超高斯分布特性[8]，在时间-空间域中实现上下行波场分离。根据公式(1a)和(1b)，本发明实施例逐道进行波场分离，从而避免假设检波器是均匀分布，且记录面是平面，更符合实际数据的处理，参见图1和图2，详见下文描述：

101：建立基于非高斯性最大化的目标函数；

考虑到本发明实施例是在时间-空间域中逐道进行上下行波场分离，省略检波器坐标，将公式(1a)和(1b)改写成向量形式如下：

Y(t)＝X(t)-gZ(t)      （5）

式中：令 $Y\left(t\right)=\left[\begin{array}{c}u\left(t\right)\\ d\left(t\right)\end{array}\right],X\left(t\right)=\left[\begin{array}{c}p\left(t\right)\\ p\left(t\right)\end{array}\right],Z\left(t\right)=\left[\begin{array}{c}-v\left(t\right)\\ v\left(t\right)\end{array}\right],$ g为一未知变量，t是时间，p(t)和v(t)分别是观测到的水检和陆检信号，u(t)和d(t)分别是待求的上行和下行信号。

通过对不同超高斯性度量进行实验，在本发明实施例中，最后选择L1范数，构建如下优化目标函数：

$\underset{g}{\arg \min {\left|\right|Y\left(t\right)\left|\right|}_1---\left(6\right)}$

为了求解公式（6）给定的优化问题，本发明实施例采用了混合L1/L2范数最小化方法^[9]该方法是一个迭代加权最小二乘方法，即利用最小L2范数方法迭代逼近最小L1范数解。在迭代的每一步，求解如下最小L2范数问题：

$\underset{g}{\arg \min {\left|\right|W(X\left(t\right)-\mathrm{gZ}\left(t\right))\left|\right|}_2^2---\left(7\right)}$

$W=\mathrm{diag}\left(\frac{1}{{(1+Y{\left(t\right)}^2/{\mathrm{\ϵ}}^2)}^{1/4}}\right)---\left(8\right)$

其中，W是加权对角矩阵，ε是阈值系数。

在初次迭代中，W通常设置为单位矩阵W，即使得初次迭代求解的结果与L2范数方法得到的结果相同。

每一步迭代中，式（7）的最小二乘解为：

$g=\frac{{\left[X\left(t\right)\right]}^T{W}^T\mathrm{WZ}\left(t\right)}{{\left[Z\left(t\right)\right]}^T{W}^T\mathrm{WZ}\left(t\right)}---\left(9\right)$

102：对于每一地震道输入地震双检信号p(t,x,y)和v(t,x,y)，获取水检和陆检信号p(t)和v(t)，分别构造出向量X(t)和Z(t)；

其中，该步骤还输入阈值系数ε、最大迭代次数N和初始迭代次数为0，阈值系数ε和最大迭代次数N的具体取值根据实际应用中的需要确定。

103：根据非高斯性最大化的目标函数和迭代加权最小二乘方法更新加权对角矩阵；

其中，该步骤具体为：

1）每一步迭代中，根据公式（9）得到变量g；

2）根据公式（5）得到上下行波场的估计；

3）根据公式（8）更新加权对角矩阵。

104：判断迭代次数是否达到最大迭代次数N，如果否，返回步骤103；如果是，判断所有地震道是否处理完毕，如果没有处理完，返回步骤102，处理下一地震道，直至所有地震道处理完毕后执行步骤105；

105：输出分离结果。

下面以具体的试验来验证本方法的可行性，详见下文描述：

在仿真实验中，本发明实施例以二维数据实验为例。利用褶积模型模拟出二维上行波（图3）和下行波（图4），通过给定一个随角度变化的空变尺度函数g(x)，利用公式（1）生成水检数据（图5）和陆检数据（图6）。

图7和图8分别给出了通过本方法处理分离得到的上行波和下行波。图9给出了真实的空变尺度函数g(x)（单实线表示），利用L2范数得到的g(x)初始估计（双实线表示），以及采用本方法迭代10次得到的g(x)最终估计（圆圈表示，即混合L1/L2范数估计方法）。通过上述实验可以看出本方法能够有效地分离上下行波场，满足了实际应用中的多种需要。

综上所述，本发明实施例利用地震信号的非高斯特性建立合适的目标函数，在时间-空间域中逐道实现上下行波场分离，可以适应检波器不是均匀分布，且记录面不是平面的情况。

参考文献

[1]Van Melle F A,Weatherborn K R.,Ghost reflections caused by energy initially reflectedabove the level of the shot.Geophysics,1953,18(4):793-804

[2]Lindsey J P.Elimination of seismic ghost reflections by means of a linear filter.Geophysics,1960,25(1):130-140

[3]Hamarbatan N S,Margrave G F.The ghost in the spectrum.69th Annual InternationalMeeting,SEG,Expanded Abstracts,1999,617-620

[4]Widmaier M.T.,Day D.,Reiser C.,and Long A.,Enhanced seismic resolution from dualsensor towed streamer acquisition,EAGE/SEG research workshop,2009

[5]Widmaier M.T.,Day D.,Reiser C.,and Long A.,Enhanced seismic resolution from dualsensor towed streamer acquisition,EAGE/SEG research workshop,2009

[6]Amundsen,L.,1993,Wave-number-based filtering of marine point sourcedata:Geophysics,58,1335–1348

[7]Sollner,W.,A.Day,and H.Tabti,2008,Space-frequency domain processing of irregulardual-sensor towed streamer data:78th Annual International Meeting,SEG,Expanded Abstracts,1078–1082.

[8]Walden,A.,1985,Non-Gaussian reflectivity,entropy,and deconvolution:Geophysics,50,2862～2888

[9]Gersztenkorn，A.，Bednar，J.B.，and Lines，L.R.，1986，Robust iterative inversionfor the one-dimensional acoustic wave equation：Geophysics，51，357－368

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于非高斯性最大化的双检信号上下行波场分离方法，其特征在于，所述方法包括以下步骤：

（1）建立基于非高斯性最大化的目标函数；

（2）对每一地震道输入地震双检信号，获取水检和陆检信号，分别构造出向量X(t)和Z(t)；

（3）根据非高斯性最大化的目标函数和迭代加权最小二乘方法更新加权对角矩阵；

（4）判断迭代次数是否达到最大迭代次数N，如果否，返回步骤（3）；如果是，判断所有地震道是否处理完毕，如果没有处理完，返回步骤（2），处理下一地震道，直至所有地震道处理完毕，执行步骤（5）；

（5）输出分离结果。

2.根据权利要求1所述的一种基于非高斯性最大化的双检信号上下行波场分离方法，其特征在于，所述建立基于非高斯性最大化的目标函数的步骤具体为：

目标函数为：

$\underset{g}{\arg \min {\left|\right|Y\left(t\right)\left|\right|}_1}$

其中，Y(t)＝X(t)-gZ(t)；

式中：令 $Y\left(t\right)=\left[\begin{array}{c}u\left(t\right)\\ d\left(t\right)\end{array}\right],X\left(t\right)=\left[\begin{array}{c}p\left(t\right)\\ p\left(t\right)\end{array}\right],Z\left(t\right)=\left[\begin{array}{c}-v\left(t\right)\\ v\left(t\right)\end{array}\right],$ g为一未知变量，t是时间，p(t)和v(t)分别是观测到的水检和陆检信号，u(t)和d(t)分别是待求的上行和下行信号。

3.根据权利要求2所述的一种基于非高斯性最大化的双检信号上下行波场分离方法，其特征在于，所述根据非高斯性最大化的目标函数和迭代加权最小二乘方法更新加权对角矩阵W的步骤具体为：

1）每一步迭代中，根据公式 $g=\frac{{\left[X\left(t\right)\right]}^T{W}^T\mathrm{WZ}\left(t\right)}{{\left[Z\left(t\right)\right]}^T{W}^T\mathrm{WZ}\left(t\right)}$ 得到变量g；

2）获取上下行波场的估计；

3）根据公式更新加权对角矩阵，ε是阈值系数。

Images