基于扩展滤波的二维小波域多次波匹配衰减方法与流程

[0001]
本发明涉及地震资料处理与分析领域，是一种应用于地震数据多次波压制的基于扩展滤波的二维小波域匹配衰减方法。


背景技术：

[0002]
在海洋地震勘探中，因为海面可近似为自由界面，而海底一般为强波阻抗界面，所以相应地震资料中通常包含较大范围的海水鸣震、交混回响等多次波。多次波的存在，降低了地震资料的信噪比，干扰了人们对有效波的识别，同时严重影响了地震资料的成像质量。近年来，随着油气勘探向崎岖海底或含有复杂构造的海域发展，相应数据的多次波压制又成为地震数据处理的关键问题。
[0003]
在崎岖复杂海域的地震数据处理中，时空域扩展滤波固然能够在一定程度上克服多次波预测不准的问题，显著提升维纳滤波类方法剔除多次波的能力，但在一次波与多次波交叉位置有效波损伤也愈加严重。


技术实现要素：

[0004]
本发明要解决的技术问题是提供基于扩展滤波的二维小波域多次波匹配衰减方法，所述方法将扩展滤波引入至二维小波域中，形成了基于小波域扩展滤波的多次波匹配衰减方法，该方法借助了二维小波变换能够实现不同频率、方向信号分离的优势，降低了一次波与多次波的耦合程度，据此来改善常规时空域多次波匹配衰减的效果。
[0005]
本发明是通过如下技术方案来实现的：
[0006]
基于扩展滤波的二维小波域多次波匹配衰减方法，所述方法具体步骤如下：
[0007]
1)创建多次波预测记录的变换记录输入原始地震记录d(x,t)，通过基于波动理论的多次波预测方法获得仅包含多次波的记录m(x,t)，其中x、t分别表示偏移距与旅行时；若预测的多次波与原始地震记录中的多次波存在振幅缩放、常相位旋转和较小的时间延迟的差异，那么原始地震记录中的多次波信号可通过预测多次波信号、预测多次波信号的希尔伯特变换以及二者时间导数之和来表示，因此计算多次波记录m(x,t)的希尔伯特变换m
h
(x,t)、时间导数m'(x,t)以及希尔伯特变换结果的时间导数[m
h
(x,t)]'；
[0008]
2)将原始地震记录、多次波记录以及多次波记录的变换记录变换到二维小波域对原始地震记录d(x,t)作二维小波变换，可得
[0009][0010]
式中，表示小波基函数，w
d
(a,τ
x
,τ
t
)代表小波系数，a为伸缩因子，x为偏移距，t为旅行时，τ
x
和τ
t
为x、t维度上的平移因子；
[0011]
将多次波记录m(x,t)及其变换记录m
h
(x,t)、m'(x,t)、[m
h
(x,t)]'变换到二维小波域，可得
[0012][0013]
由于原始地震记录中的多次波可通过预测多次波、预测多次波的希尔伯特变换以及二者时间导数之和来表示，即
[0014]
d
m
(x,t)＝l1m(x,t)+l2m
′
(x,t)+l3m
h
(x,t)+l4[m
h
(x,t)]
′ꢀꢀ
(3)
[0015]
式中，x、t分别为偏移距与旅行时；d
m
(x,t)表示原始地震记录d(x,t)中的多次波；m(x,t)为预测的多次波记录，m
h
(x,t)为m(x,t)的希尔伯特变换，m'(x,t)、[m
h
(x,t)]'分别为m(x,t)、m
h
(x,t)的时间导数；l1、l2、l3和l4分别表示m(x,t)、m
h
(x,t)、m'(x,t)与[m
h
(x,t)]'对应的权重因子；
[0016]
将式(1)代入公式(3)中可得
[0017][0018]
根据积分运算的性质，式(3)右端积分项可展开为
[0019][0020]
综合公式(2)和公式(5)可得
[0021][0022]
式中，w
d
表示原始地震记录中多次波信号的小波系数；和分别为m(x,t)、m
h
(x,t)、m'(x,t)与[m
h
(x,t)]'的小波系数；l1、l2、l3和l4表示将w
d
展开为和之和各项对应的权重因子；a为伸缩因子，τ
x
和τ
t
为x、t维度上的平移因子；
[0023]
3)二维小波域的多次波自适应相减为了保证滤波过程的稳定性，对二维小波系数矩阵中位置为k的样点，以其为中心设置矩形窗口，截取相应数据构建求取滤波因子的表达式
[0024][0025]
式中，f1、f2、f3和f4为待求取的滤波因子；w
dn
、和分别为
所截取记录块中的样点，n为截取记录块中的样点序号，其中n≥4；
[0026]
基于最小平方准则求取滤波因子，则公式(7)转化为下式描述的最小二乘问题
[0027][0028]
式中，f为滤波因子向量[f1,f2,f3,f4]
t
，m表示截取记录块和所构成的矩阵，w
d
表示记录块中各样点w
dn
组成的向量；
[0029]
公式(8)所描述的方程组为超定方程组，因此采用阻尼最小二乘方法进行求解；求取滤波因子f的表达式表示为
[0030]
f＝(m
t
m+εi)-1
(m
t
w
d
)
ꢀꢀ
(9)
[0031]
式中，m
t
为m的转置矩阵，m表示截取记录块和所构成的矩阵，i表示单位矩阵，ε为阻尼系数，ε∈(0,1)，w
d
表示记录块中各样点w
dn
组成的向量；求出滤波因子f之后，即能够实现针对原始数据小波系数目标样点的多次波匹配相减处理；
[0032]
4)对小波系数中的每个样点重复步骤3)，并将最终的结果反变换回时空域。
[0033]
本发明与现有技术相比的有益效果：
[0034]
本方法借助了二维小波变换能够实现不同频率、方向信号分离的优势，降低了一次波与多次波的耦合程度，据此改善了常规时空域多次波匹配衰减容易造成有效波损伤的问题；与曲波域多次波匹配衰减方法相比，小波域扩展滤波所耗机时约为前者的十分之一，显著地提高了计算效率。
附图说明
[0035]
图1为基于扩展滤波的二维小波域多次波匹配衰减方法的流程图；
[0036]
图2为sigsbee2b模型的速度剖面；
[0037]
图3为炮点与接收排列均位于盐丘上方的炮集记录示例；
[0038]
图4为与图3所示炮集相对应的多次波预测记录；
[0039]
图5为消除多次波的炮集记录示例；
[0040]
图6为去除的多次波干扰示例。
具体实施方式
[0041]
本发明提出的基于扩展滤波的二维小波域多次波匹配衰减方法，具体的实施过程(详细流程见图1)主要分为以下四步：1)针对预测的多次波记录创建其希尔伯特变换记录，然后分别计算出二者的时间导数记录；2)将原始地震记录、预测的多次波记录及其变换记录全部转换至小波域；3)对于小波域中的每个样点，以其为中心截取矩形数据块，通过最小平方滤波获得该样点的多次波匹配衰减结果；4)对小波系数中的每个样点重复步骤3)，并将最终的结果反变换回时空域。
[0042]
实施例1
[0043]
sigsbee2b模型是由smaart(subsalt multiples attenuation and reduction team)组织于2002年设计的标准模型(见图2)，主要用于检验地震成像与多次波压制效果。模型中的海底由浅至深变化，盐丘的顶、底均为起伏不平的复杂界面，导致原始地震记录中
富含海底及盐丘上界面有关的鸣震多次波。smaart组织提供了sigsbee2b模型的原始地震记录，相应的观测参数为：采用单边放炮的激发方式，震源位于拖缆右侧，最小偏移距为零；共激发了496炮，单炮数据的最大道数为348道，而炮间距与道间距均分别为45.72米与22.86米。现针对该模型的原始地震记录进行自由界面多次波的预测，然后应用小波域的扩展滤波方法进行多次波的匹配衰减。
[0044]
以下结合附图详细阐述本实施例的具体实施过程：
[0045]
1)创建多次波预测记录的变换记录输入图3所示的原始地震记录d(x,t)，通过自由界面多次波预测方法获得仅包含多次波的记录m(x,t)(见图4)；预测的多次波与原始地震记录中的多次波存在振幅缩放、常相位旋转和较小的时间延迟的差异，那么原始地震记录中的多次波信号可通过预测多次波信号、预测多次波信号的希尔伯特变换以及二者时间导数之和来表示，因此计算多次波记录m(x,t)的希尔伯特变换m
h
(x,t)、时间导数m'(x,t)以及希尔伯特变换结果的时间导数[m
h
(x,t)]'；
[0046]
2)将原始地震记录、多次波记录以及多次波记录的变换记录变换到二维小波域对原始地震记录d(x,t)作二维小波变换，可得
[0047][0048]
式中，表示小波基函数，w
d
(a,τ
x
,τ
t
)代表小波系数，a为伸缩因子，x为偏移距，t表示旅行时，τ
x
和τ
t
为x、t维度上的平移因子；
[0049]
将多次波记录m(x,t)及其变换记录m
h
(x,t)、m'(x,t)、[m
h
(x,t)]'变换到二维小波域，可得
[0050][0051]
由于原始地震记录中的多次波可通过预测多次波、预测多次波的希尔伯特变换以及二者时间导数之和来表示，即
[0052]
d
m
(x,t)＝l1m(x,t)+l2m
′
(x,t)+l3m
h
(x,t)+l4[m
h
(x,t)]
′ꢀꢀ
(3)
[0053]
式中，x、t分别为偏移距与旅行时；d
m
(x,t)表示原始地震记录d(x,t)中的多次波；m(x,t)为预测的多次波记录，m
h
(x,t)为m(x,t)的希尔伯特变换，m'(x,t)、[m
h
(x,t)]'分别为m(x,t)、m
h
(x,t)的时间导数；l1、l2、l3和l4分别表示m(x,t)、m
h
(x,t)、m'(x,t)与[m
h
(x,t)]'对应的权重因子；
[0054]
将式(1)代入公式(3)中可得
[0055][0056]
根据积分运算的性质，式(3)右端积分项可展开为
[0057][0058]
综合公式(2)和公式(5)可得
[0059][0060]
式中，w
d
表示原始地震记录中多次波信号的小波系数；和分别为m(x,t)、m
h
(x,t)、m'(x,t)与[m
h
(x,t)]'的小波系数；l1、l2、l3和l4表示将w
d
展开为和之和各项对应的权重因子；a为伸缩因子，τ
x
和τ
t
为x、t维度上的平移因子；
[0061]
3)二维小波域的多次波自适应相减为了保证滤波过程的稳定性，对二维小波系数矩阵中位置为k的样点，以其为中心设置矩形窗口，截取相应数据构建求取滤波因子的表达式
[0062][0063]
式中，f1、f2、f3和f4为待求取的滤波因子；w
dn
、和分别为所截取记录块中的样点，n为截取记录块中的样点序号，其中n≥4；
[0064]
基于最小平方准则求取滤波因子，则公式(7)转化为下式描述的最小二乘问题
[0065][0066]
式中，f为滤波因子向量[f1,f2,f3,f4]
t
，m表示截取记录块和所构成的矩阵，w
d
表示记录块中各样点w
dn
组成的向量；
[0067]
公式(8)所描述的方程组为超定方程组，因此采用阻尼最小二乘方法进行求解；求取滤波因子f的表达式表示为
[0068]
f＝(m
t
m+εi)-1
(m
t
w
d
)
ꢀꢀ
(9)
[0069]
式中，m
t
为m的转置矩阵，m表示截取记录块和所构成的矩阵，i表示单位矩阵，ε为阻尼系数，ε∈(0,1)，w
d
表示记录块中各样点w
dn
组成的向量；求出滤波因子f之后，即能够实现针对原始数据小波系数目标样点的多次波匹配相减处理；
[0070]
4)对小波系数中的每个样点重复步骤(3)，并将最终的结果反变换回时空域。
[0071]
通过上述步骤将预测的多次波记录(见图4)及其希尔伯特变换以及二者的时间导数均变换到小波域，通过式(9)给出的阻尼最小二乘法求取每个滤波窗的滤波因子，并采用与小波域常规滤波法相同的参数对原始地震记录(见图3)进行多次波自适应相减，所得结
果如图5所示，而去除的多次波干扰见图6。与原始地震记录(见图3)相比，箭头指向的强多次波同相轴几乎被完全消除，椭圆范围内的一次波几乎未受损伤，一次波同相轴的连续性明显变好。