基于局部相似系数的三维地震数据多次波压制方法

摘要

本发明属于地球物理勘探信号处理技术领域，涉及一种基于局部相似系数的三维地震数据多次波压制方法，是在常规的三维表面相关多次波压制技术的基础上，计算联络测线多次波贡献道集与Radon变换结果局部相似系数谱，根据相似系数差异构建加权阈值算子，利用原始数据与加权阈值算子的克罗内克积，确定稳相叠加空间，采用菲涅尔空间的叠加代替多次波联络测线的直接求和过程，实现三维海洋地震数据的高分辨率多次波预测方法，常规方法产生的空间假频得到有效压制，使得多次波预测的更加准确，提高了多次波压制的精度与准确度，为地震数据处理的后续处理与解释带来了诸多便利。

权利要求

1.一种基于局部相似系数的三维地震数据多次波压制方法，其特征在于，包括以下步骤：a、根据地震数据道头信息，提取三维海洋地震勘探采集的炮道相同的共炮点道集d(xk,yk,t,xs,ys)与共检波点道集r(xr,yr,t；xk,yk)，利用快速傅里叶变换将其转换到频率域即为共炮点道集D(xr,yr,w,xs,ys)与共检波点道集R(xr,yr,w；xs,ys)，其中，(xk,yk)表示波场中任一点的位置，(xs,ys)表示震源位置，(xr,yr)表示检波点位置，t表示时间，w表示频率；

b、利用波场传播的运动学特征，共炮点道集D(xr,yr,w,xs,ys)与共检波点道集R(xr,yr,w；xs,ys)在频率域乘积得到不同位置的纵测线多次波贡献道集：Mxy(xr,yr,xs,ys,xk,yk,w)＝r0R(xr,yr,w；xk,yk)D(xk,yk,w；xs,ys)其中，Mxy(xr,yr,xs,ys,xk,yk,w)表示频率域的纵测线多次波贡献道集；

c、不同位置的纵测线多次波贡献道集沿纵测线方向求和可得到联络测线多次波贡献道集，根据费马最小时间定理，多次波贡献道集的菲涅尔空间最能准确体现地下结构的信息，多次波贡献道集的求和过程转换为菲涅尔空间的求和过程；利用F‑K滤波方法即能够得到其菲涅尔空间，再对其求和得到联络测线多次波贡献道集，并通过傅里叶逆变换变换回时间域my(xr,yr,xs,ys,yk,t)；

Mxyf(xr,yr,xs,ys,xk,yk,w)＝F(Mxy(xr,yr,xs,ys,xk,yk,w))其中，F为倾角滤波算子，Mxyf(xr,yr,xs,ys,yk,w)为纵测线菲涅尔空间，My(xr,yr,xs,ys,yk,w)为联络测线多次波贡献道集；

d、联络测线多次波贡献道集沿联络测线方向求和即可得到最终预测的多次波；对稀疏的联络测线多次波贡献道集进行抛物Radon变换，计算联络测线多次波贡献道集a与Radon变换结果b的局部相似系数c；

其中c1，c2是上式的最小二乘解，Ai和Bi分别是由ai和bi中元素组成的对角线算子，i表示道数，S是圆滑算子，λ1，λ2是在迭代反演过程中实现快速收敛和控制物理维度的参数，其T T

中λ1＝||AA||2，λ2＝||BB||2；

e、对局部相似系数c设置高通低切阈值参数s1、s2，s1和s2实质为分段函数，当c大于s1时，取值为1，当c小于s2时，取值为0，当c为两者之间按距离加权；构建加权阈值算子W(th1,th2)，利用加权阈值算子与原始数据克罗内克积确定菲涅尔空间；

W(th1,th2)＝Y(c,s1,s2)Y为高通低切阈值函数，myf(xr,yr,xs,ys,yk,t)为联络测线多次波贡献道集的菲涅尔空间；

f、对联络测线多次波贡献道集的菲涅耳空间沿联络测线方向求和得到预测的多次波；

说明书

基于局部相似系数的三维地震数据多次波压制方法

技术领域

[0001] 本发明属于地球物理勘探信号处理技术领域，具体涉及一种三维海洋地震资料多次波压制方法，特别涉及一种基于局部相似系数的三维地震数据多次波压制方法。

背景技术

[0002] 多次波是地震波在地表或地下波阻抗差异大的界面经多次反射被地表检波器接收的地下响应，与一次波形态相似，干扰一次波的有效识别，影响地震勘探数据处理与解释
的真实性与可靠性，针对海洋地震资料，多次波的压制尤为重要。SRME方法是工业界主流的
多次波压制方法，同其它多次波压制方法相比，具有适应于复杂地形，精度高的特点。其核
心思想为利用多次波传播的运动学特征，通过共炮点道集与共检波点道集的褶积再求和的
过程实现多次波的预测，通过匹配减去达到压制多次波的目的。目前，工业界常用的多次波
压制方法局限于二维，然而多次波的三维传播与反射效应严重影响多次波的压制效果，处
理对象逐渐转向三维。三维SRME是在二维的基础上，加入联络测线多次波贡献道集，通过联
络测线多次波贡献道集求和得到预测多次波，但由于实际三维海上观测系统在联络测线方
向数据采集较为稀疏，使得在预测多次波过程中产生了空间假频，降低了地震数据的信噪
比与精确度，对后续的地震数据处理造成了一定的困难。

[0003] Verschuur(Van Dedem E J,Verschuur D J.3D surface‑related multiple prediction:Asparse inversion approach[J].Seg Technical Program Expanded 
Abstracts,2005,18(1):2061.)提出了基于稀疏反演的多次波压制方法，对稀疏的联络测
线多次波贡献道集利用柯西约束反演在模型空间重建，再通过振幅和相位进行校正叠加得
到预测多次波。Ketil Hokstad等人(Hokstad K,Sollie R.3D surface‑related multiple 
elimination using parabolic sparse inversion[J].GEOPHYSICS,2006,71(6):V145‑
V152.)在Verschuur方法的基础上，根据地震波场的抛物线特征进行稀疏反演。Anatoly 
Baumstein等人(Baumstein A.3D SRME:Data reconstruction and application to 
field data[J].SEG Technical Program Expanded Abstracts,1999,23(1).)提出联络测
线方向数据重建得到3D SRME需要的密集采样的数据。方云峰，石颖等人(方云峰.深海地震
资料全三维表面多次波预测技术研究[D].2015.)提出利用求和孔径优化后的多次波贡献
道集进行多次波压制。但是，由于上述这些方法采用的稀疏反演方法，数据重构或者优化孔
径的不稳定性以及主观控制因素过多，使得对于不同的复杂地形影响较大。

[0004] 局部相似系数可测量两信号间各部分的相似性(Fomel S.Local seismic attributes[J].Geophysics,2007,72(3):A29.)，其实现方式是在全局相似系数的基础上，
结合最小二乘反演与局部互相关算法，目前主要应用于图像处理、图像元素校正等领域。

发明内容

[0005] 本发明的目在于提供一种基于局部相似系数的三维地震数据多次波压制方法，以克服常规的三维表面相关多次波压制技术由于联络测线方向采样稀疏，导致多次波贡献道
集求和结果存在大量空间假频的缺点，实现三维海洋地震数据的高分辨率多次波预测，使
常规方法产生的空间假频得到有效压制。

[0006] 本发明的目的是通过以下技术方案实现的：

[0007] 首先，了解三维海洋地震勘探数据采集观测系统所采集的地震数据特征，从中提取出共炮点道集与共检波点道集，共炮点道集与共检波点道集褶积得到不同位置的纵测线
多次波贡献道集；通过F‑K滤波方法仅保留纵测线多次波贡献道集的菲涅尔空间，对菲涅尔
空间沿纵测线方向求和得到联络测线多次波贡献道集；接下来，对稀疏的联络测线多次波
贡献道集进行抛物Radon变换，计算Radon变换结果与联络测线多次波贡献道集的局部相似
系数。对局部相似系数设置高通低切阈值参数，构建加权阈值算子，利用加权阈值算子与原
始数据克罗内克积确定菲涅尔空间，用菲涅尔空间的求和代替常规的联络测线多次波贡献
道集直接求和，得到最终预测的多次波。

[0008] 一种基于局部相似系数的三维地震数据多次波压制方法，包括以下步骤：

[0009] a、根据地震数据道头信息，提取三维海洋地震勘探采集的炮道相同的共炮点道集d(xk,yk,t,xs,ys)与共检波点道集r(xr,yr,t；xk,yk)，利用快速傅里叶变换将其转换到频率
域即为共炮点道集D(xr,yr,w,xs,ys)与共检波点道集R(xr,yr,w；xs,ys)。其中(xk,yk)表示波
场中任一点的位置，(xs,ys)表示震源位置，(xr,yr)表示检波点位置，t表示时间，w表示频
率。

[0010] b、利用波场传播的运动学特征，共炮点道集D(xr,yr,w,xs,ys)与共检波点道集R(xr,yr,w；xs,ys)在频率域乘积得到不同位置的纵测线多次波贡献道集：

[0011] Mxy(xr,yr,xs,ys,xk,yk,w)＝r0R(xr,yr,w；xk,yk)D(xk,yk,w；xs,ys)

[0012] 其中，Mxy(xr,yr,xs,ys,xk,yk,w)表示频率域的纵测线多次波贡献道集。

[0013] c、不同位置的纵测线多次波贡献道集沿纵测线方向求和可得到联络测线多次波贡献道集，根据费马最小时间定理，多次波贡献道集的菲涅尔空间最能准确体现地下结构
的信息，多次波贡献道集的求和过程转换为菲涅尔空间的求和过程。利用F‑K滤波方法即能
够得到其菲涅尔空间，再对其求和得到联络测线多次波贡献道集，并通过傅里叶逆变换变
换回时间域my(xr,yr,xs,ys,yk,t)。

[0014] Mxyf(xr,yr,xs,ys,xk,yk,w)＝F(Mxy(xr,yr,xs,ys,xk,yk,w))

[0015]

[0016] 其中，F为倾角滤波算子，Mxyf(xr,yr,xs,ys,yk,w)为纵测线菲涅尔空间，My(xr,yr,xs,ys,yk,w)为联络测线多次波贡献道集。

[0017] d、联络测线多次波贡献道集沿联络测线方向求和即可得到最终预测的多次波。然而由于联络测线分布过于稀疏，对其直接求和或采用上述所说的F‑K滤波方法得到菲涅尔
空间信息会带来严重的空间假频。这里对稀疏的联络测线多次波贡献道集进行抛物Radon
变换，计算联络测线多次波贡献道集a与Radon变换结果b的局部相似系数c。

[0018]

[0019]

[0020]

[0021] 其中c1，c2是上式的最小二乘解，Ai和Bi分别是由ai和bi中元素组成的对角线算子，i表示道数，S是圆滑算子，λ1，λ2是在迭代反演过程中实现快速收敛和控制物理维度的参数，
T T
其中λ1＝||AA||2，λ2＝||BB||2。

[0022] e、对局部相似系数c设置高通低切阈值参数s1、s2，实质为分段函数，大于s1等于1，小于s2等于0，两者之间按距离加权。构建加权阈值算子W(th1,th2)，利用加权阈值算子与原
始数据克罗内克积确定菲涅尔空间。

[0023] W(th1,th2)＝Y(c,s1,s2)

[0024]

[0025] Y为高通低切阈值函数，myf(xr,yr,xs,ys,yk,t)为联络测线多次波贡献道集的菲涅尔空间。

[0026] f、对联络测线多次波贡献道集的菲涅耳空间沿联络测线方向求和得到预测的多次波。

[0027]

[0028] 与现有技术相比，本发明的有益效果：

[0029] 本发明应用local similarity优化3D SRME，对联络测线多次波贡献道集作抛物Radon变换，计算联络测线多次波贡献道集与Radon变换结果局部相似系数谱(local 
similarity)，根据相似系数差异构建加权阈值算子，利用原始数据与加权阈值算子克罗内
克积，确定菲涅尔空间，再对菲涅尔空间沿联络测线方向求和的过程，解决了由于联络测线
方向采集数据较为稀疏产生的空间假频问题，使得多次波预测的更加准确，提高了多次波
压制的精度与准确度，为地震数据处理的后续处理与解释带来了诸多便利。

[0030] 本发明有以下优点：

[0031] 1、基于局部相似系数确定菲涅尔空间不依赖于地下的假设和信息，属于数据驱动，仅需要获取震源子波的信息以及震源和检波器的空间分布就可以实现菲涅尔空间优
化，使预测结果不受空间假频影响。

[0032] 2、该方法同常规表面多次波压制方法相比，克服了联络测线采集较为稀疏的限制，具有更高的精度及准确度，稳定性高，且计算效率高，极大的降低了地震数据处理与解
释的成本。

附图说明

[0033] 图1基于局部相似系数的三维地震数据多次波压制方法流程图；

[0034] 图2三维海上地震数据采集观测系统示意图；

[0035] 图3三维多次波传播路径示意图；

[0036] 图4菲涅尔空间说明示意图；图4a密集采样的联络测线多次波贡献道集，图4b稀疏采样的联络测线多次波贡献道集，图4c密集采样的联络测线多次波贡献道集，密集采样的
菲涅尔空间，稀疏采样的菲涅尔空间的求和结果对比；

[0037] 图5纵测线多次波贡献道集示意图；

[0038] 图6纵测线多次波贡献道集F‑K滤波结果示意图；

[0039] 图7联络测线多次波贡献道集示意图；

[0040] 图8联络测线多次波贡献道集抛物Radon变换示意图；

[0041] 图9local similarity示意图；

[0042] 图10三维倾斜模型数值算例；图10a三维倾斜模型，图10b常规三维SRME的多次波预测结果，图10c基于局部相似系数拾取菲涅尔空间的表面多次波预测结果。

具体实施方式

[0043] 下面结合附图和实例对本发明进行进一步的详细说明。

[0044] 本发明是在常规的三维表面相关多次波压制技术(3D surface‑related multiples elimination,3D SRME)的基础上，计算联络测线多次波贡献道集与Radon变换
结果局部相似系数谱(local similarity)，根据相似系数差异构建加权阈值算子，利用原
始数据与加权阈值算子的克罗内克积，确定稳相叠加空间(菲涅尔空间)，采用菲涅尔空间
的叠加代替多次波联络测线的直接求和过程，实现三维海洋地震数据的高分辨率多次波预
测方法，常规方法产生的空间假频得到有效压制。基于局部相似系数的三维地震数据多次
波压制方法是通过Seismic Unix和MATLAB平台实现的。

[0045] 本发明基于局部相似系数的三维地震数据多次波压制方法，包括以下步骤：

[0046] a、三维海洋地震勘探述采集观测系统为拖缆船拖拽震源与布满检波器的多条测线前进采集，如图2所示，采集到的数据需要根据互易定理，按照地震数据道头信息，补全并
提取三维海洋地震勘探采集的炮道相同的共炮点道集d(xk,yk,t,xs,ys)与共检波点道集r
(xr,yr,t；xk,yk)，利用快速傅里叶变换将其转换到频率域即为共炮点道集D(xr,yr,w,xs,ys)
与共检波点道集R(xr,yr,w；xs,ys)。其中(xk,yk)表示波场中任一点的位置，(xs,ys)表示震源
位置，(xr,yr)表示检波点位置，t表示时间，w表示频率。

[0047] b、利用波场传播的运动学特征，如图3所示，共炮点道集D(xr,yr,w,xs,ys)与共检波点道集R(xr,yr,w；xs,ys)在频率域乘积得到不同位置的纵测线多次波贡献道集，如图5所示：

[0048] Mxy(xr,yr,xs,ys,xk,yk,w)＝r0R(xr,yr,w；xk,yk)D(xk,yk,w；xs,ys)

[0049] 其中，Mxy(xr,yr,xs,ys,xk,yk,w)表示频率域的纵测线多次波贡献道集。

[0050] c、不同位置的纵测线多次波贡献道集沿纵测线方向求和可得到联络测线多次波贡献道集，根据费马最小时间定理，多次波贡献道集的菲涅尔空间最能准确体现地下结构
的信息，多次波贡献道集的求和过程转换为菲涅尔空间的求和过程。利用F‑K滤波方法即能
够得到其菲涅尔空间，如图6所示，再对其求和得到联络测线多次波贡献道集，并通过傅里
叶逆变换变换回时间域my(xr,yr,xs,ys,yk,t)。

[0051] Mxyf(xr,yr,xs,ys,xk,yk,w)＝F(Mxy(xr,yr,xs,ys,xk,yk,w))

[0052]

[0053] 其中，F为倾角滤波算子，Mxyf(xr,yr,xs,ys,yk,w)为纵测线菲涅尔空间，My(xr,yr,xs,ys,yk,w)为联络测线多次波贡献道集。

[0054] d、联络测线多次波贡献道集沿联络测线方向求和即可得到最终预测的多次波。然而由于联络测线分布过于稀疏，如图4所示，对其直接求和或采用上述所说的F‑K滤波方法
得到菲涅尔空间信息会带来严重的空间假频。这里对稀疏的联络测线多次波贡献道集进行
抛物Radon变换，计算联络测线多次波贡献道集a(图7)与Radon变换结果b(图8)的局部相似
系数c(图9)。

[0055]

[0056]

[0057]

[0058] 其中c1，c2是上式的最小二乘解，Ai和Bi分别是由ai和bi中元素组成的对角线算子，i表示道数，S是圆滑算子，λ1，λ2是在迭代反演过程中实现快速收敛和控制物理维度的参数，
T T
其中λ1＝||AA||2，λ2＝||BB||2。

[0059] e、对局部相似系数c设置高通低切阈值参数s1、s2，实质为分段函数，大于s1等于1，小于s2等于0，两者之间按距离加权。构建加权阈值算子W(th1,th2)，利用加权阈值算子与原
始数据克罗内克积确定菲涅尔空间。

[0060] W(th1,th2)＝Y(c,s1,s2)

[0061]

[0062] Y为高通低切阈值函数，myf(xr,yr,xs,ys,yk,t)为联络测线多次波贡献道集的菲涅尔空间。

[0063] f、对联络测线多次波贡献道集的菲涅耳空间沿联络测线方向求和得到预测的多次波。该预测的多次波与常规SRME方法预测的多次波相比，精度更高，压制空间假频效果较
为明显。

[0064]

[0065] 实施例1

[0066] a、地层模型为2000m*1000m的倾斜地层模型，如图10a所示，地层模型有三层界面，即四层介质。第一层界面在(0，0)处的深度为800m，第一层介质的速度为1200m/s，密度为
3
1000kg/m ，沿联络测线方向倾斜5°；第二层界面在(0，0)处的深度为1300m，第二层介质的
3
速度为1500m/s，密度为2000kg/m ，沿纵测线方向倾斜10°；第三层界面在(0，0)处的深度为
3
2200m，第三层介质的速度为2000m/s，密度为2500kg/m ，为水平层；第四层介质的速度为
3
2500m/s，密度为3000kg/m 。观测系统设置为11条测线长2000m，彼此间隔100m(‑500m—
500m)均匀分布，测线上间隔25m均匀分布着检波器接收震源激发的地震波，震源在彼此间
隔25m的位置处依次激发，得到炮检相同的地震勘探数据。本实施实例的目的是通过模拟数
据来测试本方法对多次波压制的准确性；

[0067] b、在上述炮检相同的地震数据中，根据数据道头信息提取出共炮点道集d(xk,yk,t,xs,ys)与共检波点道集r(xr,yr,t；xk,yk)，利用快速傅里叶变换将其转换到频率域即为共
炮点道集D(xr,yr,w,xs,ys)与共检波点道集R(xr,yr,w；xs,ys)。其中(xk,yk)表示波场中任一
点的位置，(xs,ys)表示震源位置，(xr,yr)表示检波点位置，t表示时间，w表示频率。

[0068] c、利用波场传播的运动学特征，共炮点道集D(xr,yr,w,xs,ys)与共检波点道集R(xr,yr,w；xs,ys)在频率域乘积得到不同位置的纵测线多次波贡献道集：

[0069] Mxy(xr,yr,xs,ys,xk,yk,w)＝r0R(xr,yr,w；xk,yk)D(xk,yk,w；xs,ys)

[0070] 其中，Mxy(xr,yr,xs,ys,xk,yk,w)表示频率域的纵测线多次波贡献道集,r0等于‑1。

[0071] d、不同位置的纵测线多次波贡献道集沿纵测线方向求和可得到联络测线多次波贡献道集，根据费马最小时间定理，多次波贡献道集的菲涅尔空间最能准确体现地下结构
的信息，多次波贡献道集的求和过程转换为菲涅尔空间的求和过程。利用F‑K滤波方法即能
够得到其菲涅尔空间，再对其求和得到联络测线多次波贡献道集，并通过傅里叶逆变换变
换回时间域my(xr,yr,xs,ys,yk,t)。其中，F‑K滤波参数选择为slopes＝‑4000,‑50,50,4000，
amps＝0,1,1,0其中slopes表示波数k，amps表示振幅。

[0072] Mxyf(xr,yr,xs,ys,xk,yk,w)＝F(Mxy(xr,yr,xs,ys,xk,yk,w))

[0073]

[0074] 其中，F为倾角滤波算子，Mxyf(xr,yr,xs,ys,yk,w)为纵测线菲涅尔空间，My(xr,yr,xs,ys,yk,w)为81条联络测线多次波贡献道集。

[0075] e、联络测线多次波贡献道集沿联络测线方向求和即可得到最终预测的多次波。然而由于联络测线分布过于稀疏，对其直接求和或采用上述所说的F‑K滤波方法得到菲涅尔
空间信息会带来严重的空间假频。这里对稀疏的联络测线多次波贡献道集进行抛物Radon
变换，Radon变换参数设置为参考最大偏移距为2000m，最小时差为200ms，最大时差为
1000ms。

[0076] 计算联络测线多次波贡献道集a与Radon变换结果b的局部相似系数c。

[0077]

[0078]

[0079]

[0080] 其中，c1，c2是上式的最小二乘解，Ai和Bi分别是由ai和bi中元素组成的对角线算子，i表示道数(1—11)，S是圆滑算子，沿时间方向为11，沿偏移距方向为5；λ1，λ2是在迭代反
演过程中实现快速收敛和控制物理维度的参数，λ1为0.0433乘以a的特征值，λ2为0.0076乘
以b的特征值。

[0081] f、对局部相似系数c设置高通低切阈值参数s1、s2，实质为分段函数，大于s1等于1，小于s2等于0，两者之间按距离加权。构建加权阈值算子W(th1,th2)，利用加权阈值算子与原
始数据克罗内克积确定菲涅尔空间。

[0082] W(th1,th2)＝Y(c,s1,s2)

[0083]

[0084] Y为高通低切阈值函数，myf(xr,yr,xs,ys,yk,t)为联络测线多次波贡献道集的菲涅尔空间。

[0085] 其中，s1＝1.5*10‑5,s2＝0.01。

[0086] g、对联络测线多次波贡献道集的菲涅耳空间沿联络测线方向求和得到预测的多次波。该预测的多次波，如图10c与常规SRME方法预测的多次波如图10b相比，精度更高，压
制空间假频效果较为明显。

[0087]