CN110456417A - 一种地震数据多次波压制方法 - Google Patents

Abstract

本发明公开了一种地震数据多次波压制方法，包括以下步骤:S1，输入三维数据d(offx,t)，依据三维空间坐标信息，道集抽取为三维面元道集d_bin(y,x,t)；S2，计算非线性趋势拟合系数p(x)，p(y)，写出变换算子L_x，L_y，L_xy，计算Radon域数据M；S3，Radon域数据切除处理，得到多次波模型m_ra，并计算变换得到时空域多次波模型S4，相减得到多次波压制后的三维面元道集p；S5，将道集按标量偏移距大小依次排列，得到多次波压制后三维数据，通过低频结果约束高频的计算，低频结果在运算中不会产生假频，在采样稀疏的情况下仍能取得比较好的高分辨率效果；同时避免了常规迭代方法提高分辨率的过程，大大提高了计算效率。

Description

一种地震数据多次波压制方法

技术领域

本发明涉及勘探地震信号处理领域，具体涉及一种地震数据多次波压制方法。

背景技术

海洋地震数据采集时，由于受到海水表面强反射系数的影响，海洋地震数据中会存在强烈的多次波；多次波的存在严重影响着地震数据速度分析和成像；如何有效去除多次波是地震资料处理中的一个重要研究内容。

Radon变换是工业界多次波压制最为常用的有效手段之一，同时也被广泛用于地震数据重建，波场分离，速度分析等。Radon变换的二维算法没有考虑地震波传播的三维效应，常规三维Radon变换算法没有考虑地层性质的各向异性，分辨率也较低，而且常规Radon变换只是沿着特定路径的叠加求和，并没有考虑横向的振幅变化，因此也存在不保幅的问题。

发明内容

本发明的目的是在于提供一种地震数据多次波压制方法，利用非线性趋势预测理论，对地震数据横向振幅变化进行拟合并存储；利用计算时低频的计算结果对高频的计算进行约束提高分辨率，实现地震数据三维高分辨率保幅Radon变换，得到更精确的多次波模型，进而达到更好的多次波压制效果；利用椭圆模型应对地层参数各向异性所带来的地震数据的方位效应，能够更好的进行复杂数据多次波的压制。

为解决上述发明的目的，本发明提供技术方案如下：

一种地震数据多次波压制方法，其特征在于，包括以下步骤:

S1，输入三维数据d(offx,t)，依据三维空间坐标信息，道集抽取为三维面元道集d_bin(y,x,t)：输入三维数据d(offx,t)，对于给定的以标量偏移距offx排列的三维共中心点道集地震数据d(offx,t)，依据该道集每个标量偏移距offx所对应的纵测线和非纵测线空间坐标，将其转变抽取为三维面元道集d_bin(y,x,t)；其中，y为地震数据非纵测线方向的偏移距，x为地震数据纵测线方向的偏移距，t为纵坐标时间；

S2，计算非线性趋势拟合系数p(x)，p(y)，写出变换算子L_x，L_y，L_xy，计算Radon域数据M；计算Radon变换算子L_x,L_y，L_xy，同时将地震数据d_bin(y,x,t)变换到频率域，并将如上结果代入公式中，求取得到频率域Radon域数据M；并将其做反变换至时间域，得到时间域的Radon域地震数据d_ra，依据一次波和多次波在Radon域的空间位置关系，对其进行切除操作，得到Radon域多次波模型m_ra；

S3，Radon域数据切除处理，得到多次波模型m_ra，并依据公式变换得到时空域多次波模型依据一次波和多次波在Radon域的空间位置关系，对其进行切除操作，得到Radon域多次波模型m_ra；将其变换到频率域M_ra，结合Radon变换算子，进行反变换，得到频率域多次波模型M_ra-freq，将M_ra-freq变换回时间域，得到时空域的多次波模型

S4，相减得到多次波压制后的三维面元道集p；将得到的多次波模型和三维面元道集d_bin做减法，得到多次波压制后的三维面元道集p：

S5，将道集按标量偏移距大小依次排列，得到多次波压制后三维数据；将p按照纵测线和非纵测线偏移距大小，计算标量偏移距大小，并按照标量偏移距大小一次排列地震道，得到标量偏移距offx排列的、多次波压制后的三维共中心点道集地震数据。

进一步地，所述步骤S1中，三维面元道集d_bin(y,x,t)基于椭圆模型，其表达为：

d(y,x,t)＝∫∫∫m(q_xy,q_y,q_x,τ＝t-q_xx²-q_yy²-q_xyxy,)dq_xdq_ydq_xy (1)

其中，q_xy为考虑了纵测线和非纵测线方向椭圆效应的综合曲率参数；

q_y为y方向曲率参数；q_x为x方向曲率参数；τ为Radon域的时间；

对公式(1)离散并变换到频率域可得到：

可将上式写成矩阵算子相乘的形式：

D＝L_xML_yL_xy (3)

进一步地，所述步骤S1中，对应于公式(1)，考虑到同相轴振幅变化、离散化的反变换公式为：

进一步地，步骤S2中，计算非线性趋势拟合系数p(x)，p(y)，其公式为：

所述计算Radon域数据M，其公式为：

式中，M_n,D_n分别是第n个频率的多次波数据和地震数据；

L_x，L_xy，L_y分别为Radon变换计算时，包含了非线性趋势拟合系数的纵测线和非纵测线算子，可用如下公式表示：

进一步地，所述步骤S5中，其判断方法为：

对公式(6)中的Wn,Un,Vn由上一个频率的计算结果得出：

W_ii(ω_n)＝||M_i(ω_n-1)|| (8)

其中，Mi为在上一频率运算得到计算结果，W若为纵测线方向的加权矩阵，则i的范围为1～Nqx；

若为纵测线和非纵测线两个方向的加权矩阵，则i的范围为1～N qxy,若为横测线方向的加权矩阵V，则i的范围为1～Nqy。

本发明相对于现有技术的有益效果是：

1、三维方法考虑了地震波场的三维传播效应，能够更加精确的估计多次波的空间传播旅行时；

2、椭圆模型能够更好的应对地层性质存在各向异性时，采集得到的存在方位效应的地震数据，基于椭圆模型的Radon变换算子，在处理存在方位效应的地震数据能实现更好的Radon域能量集中；

3、非线性趋势预测能够较好的保存地震数据空间的振幅变化，从而能够更加精确的估计多次波的振幅；

4、利用低频结果约束高频的计算，更看重低频结果对整体结果的约束，低频结果在运算中不会产生假频，在采样稀疏的情况下仍能取得比较好的高分辨率效果；

5、利用低频计算结果对高频进行约束，避免了常规迭代方法提高分辨率的过程，可以大大提高计算效率。

附图说明

图1是本发明实施例中地震数据多次波压制方法的流程图；

图2是本发明实施例中标量偏移距为横坐标的三维道集地震数据图；

图3是本发明实施例中以非纵测线偏移距大小依次排列的三维小面元道集地震数据图；

图4是传统最小二乘三维抛物线Radon变换压制的多次波面元道集图；

图5是本发明实施例中最小二乘三维抛物线Radon变换多次波压制结果面元道集图；

图6是本发明实施例中标量偏移距为横坐标的传统最小二乘三维抛物线Radon变换多次波压制结果图；

图7是本发明实施例中压制的多次波面元道集图；

图8是本发明实施例中多次波压制结果小面元道集图；

图9是本发明实施例中以标量偏移距为横坐标的多次波压制结果图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明具体实施例及相应的附图对本发明技术方案进行清楚、完整地描述。显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

结合图1所示，在均匀介质中，地震波前面是一个球面，有效波的波前面遇到地下反射面发生反射的时候，反射波前面也是标准的球面，当地震反射波被接收之后，对于相同的反射时间，接收到地震波的检波器表现为一个圆形的组成，因此将三维Radon变换中地震波旅行时和非纵测线、纵测线方向偏移距的关系可以表示为t＝τ+q_xx²+q_yy²。但是当地下介质不满足密度均匀的情况下，这个模型也会变现出一定的误差。对于三维地震数据，如果地下介质在不同方向上密度变化较为明显，地震波旅行时相同时，三维的检波器分布不再是标准的圆，而是在圆的基础上发生了平滑的调整。当地震数据存在方位效应时，即不同方向上地层参数变现为各向异性，我们需要用新的模型来适用该种地层接收到的地震数据，当三维地震数据存在方位效应时，地震波的形态不再是标准的圆，而表现为椭圆t＝τ+q_xx²+q_yy²+q_xyxy，可以通过椭圆模型得到基于椭圆模型的Radon变换算子，在处理存在方位效应的地震数据实现更好的Radon域能量集中。

基于椭圆模型，可以将三维Radon变换表示为：

d(y,x,t)＝∫∫∫m(q_xy,q_y,q_x,τ＝t-q_xx²-q_yy²-q_xyxy,)dq_xdq_ydq_xy (1)

其中，qxy为考虑了纵测线和非纵测线方向椭圆效应的综合曲率参数，

qy为y方向曲率参数，qx为x方向曲率参数，τ为Radon域的时间。

相比于常规三维Radon变换，引入了qxy来拟合地震波的真实传播路径，这对存在方向各向异性的地震资料收敛地震波能量是有帮助的，将公式(1)离散并变换到频率域可得到：

可将上式写成矩阵算子相乘的形式：

D＝L_xML_yL_xy (3)

在进行Radon变换求和时，通过在求和过程中加入拟合系数p(x)、p(y)，可以考虑地震数据沿x和y方向空间振幅变化，从而估计出更加准确的多次波模型，得到更好的多次波压制效果。因此，对三维共中心点道集地震数据d_bin(y,x,t)，可进行保幅三维Radon变换，得到Radon域地震数据m(j,i,q_xy,q_y,q_x,τ)，其中，j,i分别为非纵测线和纵测线方向的非线性趋势拟合阶数，对应于公式(1)，考虑了同相轴振幅变化了的、离散化的反变换公式为：

非线性趋势拟合系数p(x),p(y)可以依据下式进行计算：

进一步地，为提高m的计算精度，可在频率域依据如下公式计算高分辨率频率域Radon域地震数据M：

其中，Mn,Dn分别是第n个频率的多次波数据和地震数据，L_x，L_xy，L_y分别为Radon变换计算时，包含了非线性趋势拟合系数的纵测线和非纵测线算子，可用如下公式表示：

公式(6)中的W_n,U_n,V_n由上一个频率的计算结果得到，可由下式得出：

W_ii(ω_n)＝||M_i(ω_n-1)|| (8)

其中，M_i为在上一频率运算得到计算结果，W若为纵测线方向的加权矩阵，则i的范围为1～N_qx，若为纵测线和非纵测线两个方向的加权矩阵，则i的范围为1～N_qxy,若为横测线方向的加权矩阵V，则i的范围为1～N_qy。

利用低频运算结果约束高频计算的方法，避免了常规迭代方法提高分辨率的过程，可以大大提高计算效率。

结合图2-图9所示，图2为理论数据是标量偏移距为横坐标的三维道集地震数据，该道集中，一次波和多次波受到方位各向异性的影响，同相轴有一定的离散。

图3为以非纵测线偏移距大小依次排列的三维小面元道集地震数据。该道集中，沿着横测线方向较平的是一次波对应的同相轴，一次波同相轴沿y方向有较小幅度的弯曲，是由横测线间距引起的；沿y方向弯曲，且同相轴本身具有一定曲率的是多次波所对应的同相轴。该道集中，对同相轴沿纵测线方向的振幅加以变化，以模拟实际数据沿偏移距方向的振幅变化。

图4为传统最小二乘三维抛物线Radon变换压制的多次波面元道集。由于最小二乘抛物线Radon变换方法分辨率较低，不能较好的将一次波和多次波分离开，导致变换回的多次波中有一次波的信息，从而使得最终的多次波压制结果不保幅；图中可以看到有明显的一次波能量。

图5为最小二乘三维抛物线Radon变换多次波压制结果面元道集。该道集可以看出，由于传统最小二乘三维抛物线Radon变换方法未考虑同相轴振幅沿空间方向的变化，使得多次波模型振幅和原数据多次波振幅不匹配，最终导致多次波能量有残余。

图6为标量偏移距为横坐标的传统最小二乘三维抛物线Radon变换多次波压制结果。该道集中可以看到有明显的多次波能量残余。

图7为本发明方法压制的多次波面元道集。可以看出，多次波中一次波能量与图4相比明显减少，显示了本文方法分辨率高的优势。

图8为本发明方法多次波压制结果小面元道集。可以看出，由于本方法利用非线性趋势拟合，充分考虑了振幅随偏移距的变化情况，使得多次波压制后的结果中，多次波残余能量明显降低。

图9为以标量偏移距为横坐标的本发明方法多次波压制结果。可以看出，多次波压制彻底，同时一次波能量得到了很好的保护。

以上所述仅为本发明的实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (5)

1.一种地震数据多次波压制方法，其特征在于，包括以下步骤:

S1，输入三维数据d(offx,t)，依据三维空间坐标信息，道集抽取为三维面元道集d_bin(y,x,t)；

S2，计算非线性趋势拟合系数p(x)，p(y)，写出变换算子L_x，L_y，L_xy，计算Radon域数据M；

S3，Radon域数据切除处理，得到多次波模型m_ra，并依据公式变换得到时空域多次波模型

S4，相减得到多次波压制后的三维面元道集p；

S5，将道集按标量偏移距大小依次排列，得到多次波压制后三维数据。

2.如权利要求1所述的一种地震数据多次波压制方法，其特征于，所述步骤S1中，三维面元道集d_bin(y,x,t)基于椭圆模型，其表达为

d(y,x,t)＝∫∫∫m(q_xy,q_y,q_x,τ＝t-q_xx²-q_yy²-q_xyxy,)dq_xdq_ydq_xy (1)

其中，q_xy为考虑了纵测线和非纵测线方向椭圆效应的综合曲率参数；

q_y为y方向曲率参数；q_x为x方向曲率参数；τ为Radon域的时间；

对公式(1)离散并变换到频率域可得到：

可将上式写成矩阵算子相乘的形式：

D＝L_xML_yL_xy (3) 。

3.如权利要求2所述的一种地震数据多次波压制方法，其特征在于，所述步骤S1中，对应于公式(1)，考虑到同相轴振幅变化、离散化的反变换公式为：

。

4.如权利要求2所述的一种地震数据多次波压制方法，其特征在于，所述步骤S2中，计算非线性趋势拟合系数p(x)，p(y)，其公式为：

所述计算Radon域数据M，其公式为：

式中，M_n,D_n分别是第n个频率的多次波数据和地震数据；

L_x，L_xy，L_y分别为Radon变换计算时，包含了非线性趋势拟合系数的纵测线和非纵测线算子，可用如下公式表示：

。

5.如权利要求4所述的一种地震数据多次波压制方法，其特征在于，所述步骤S5中，其判断方法为：

对公式(6)中的Wn,Un,Vn由上一个频率的计算结果得出：

W_ii(ω_n)＝||M_i(ω_n-1)|| (8)

其中，Mi为在上一频率运算得到计算结果，W若为纵测线方向的加权矩阵，则i的范围为1～Nqx；

若为纵测线和非纵测线两个方向的加权矩阵，则i的范围为1～N qxy,若为横测线方向的加权矩阵V，则i的范围为1～Nqy。