面波, 又称为地滚波。在地震勘探领域, 面波主要是瑞利波, 由 P 波和 SV 波耦合而成, 在自由表面传播, 振幅随深度呈指数衰减。瑞利波是地震干扰波的一种, 具有高振幅、低频率、低视速度及强频散等特点, 在地震记录中呈扇形分布, 与有效波信号重叠, 降低地震数据的信噪比, 增加了地震资料处理的难度。利用其特点, 已经有数种方法实现对瑞利波的压制, 如高通滤波^[1]、f-k滤波^[2]、Radon 变换^[3]等方法, 后又引入小波变换^[4-5]、Rid-gelet变换^[6]、Curvelet变换^[7]等方法。以上6种方法都用不同的变换方法把原数据变换到另一个域, 利用瑞利波与有效波在新的域分离的特性实现滤波。高通滤波是频率域滤波, 去除的是低频部分, 可能会损伤部分低频有效波信号。f-k 滤波将数据转换到二维 f-k 域(频率‒波数域), 但需要规则的空间采样间隔, 同时也存在分离不完全的问题。Radon变换分为线性 Radon 变换、抛物线 Radon 变换和双曲线 Radon 变换等。与 f-k 滤波类似, 对于地质构造复杂的区域, 瑞利波趋向于在转换域中分散, 与有效波混叠在一起, 小波变换则转换到小波域, 具有良好的时域和频域的局部化特性, 与传统的傅里叶变换相比有优势, 但在高维线奇异的处理上存在缺陷。Ridgelet 变换在小波变换的基础上做了改进, 在对线和超平面的奇异性表示上优于小波变换。Curvelet 变换又对 Ridgelet 变换做了改进, 基于曲线分割的思想, 在无限小尺度内将每一段曲线近似为直线, 能更好地处理高维信号的线奇异特征。

本文引入超虚干涉法, 实现瑞利波的压制。干涉法的基础是用互相关的方法获得两个接收点之间的介质响应, 其基本形式为: 对同一扩散场的两个不同接收点记录做互相关, 可以得到当其中一个接收点作为虚源时, 另一个接收点观察到的格林函数。Wapenaar 等^[8]详细推导了声波和弹性波介质中用于地震干涉的格林函数表示定理, 为地震干涉法提供了理论基础。Dong 等^[9]基于互相关的方法首次实现瑞利波的提取, 并用最小平方滤波方法有效地去除瑞利波。Halliday 等^[10-15]对瑞利波干涉法做了大量的理论分析和数据验证, 并运用驻点相分析方法, 指出在自由表面存在的驻点是表面源干涉法提取瑞利波的关键, 且在强散射区域, 干涉法的处理效果比f-k方法更好。

超虚干涉法是 Bharadwaj 等^[16-17]和 Mallison 等^[18]处理折射波干涉法时首先提出的, 基于互相关得到的虚折射记录存在于虚源至检波点之间, 与实际地震资料炮点至检波点的记录不符。针对这一问题, 超虚干涉法在互相关之后增加一褶积步骤, 即对不同炮点的两个检波点记录做互相关并叠加, 将叠加记录与原始记录做褶积并二次叠加, 以获得炮点至检波点的预测折射波数据。基于互相关的瑞利波干涉法同样只能预测检波点间的瑞利波。本文引入超虚干涉法, 以得到炮点至检波点的预测瑞利波, 并利用基于多道L₁范数自适应相减方法压制瑞利波, 在合成数据的测试中取得一定的效果。

1 方法原理

1.1 超虚干涉法基本原理

如图 1(a)所示, s₁ 为震源, g 和 g′ 均为检波点, 瑞利波在自由表面传播。在g检波点的瑞利波地震记录在频率域可以近似地表示为

$u({{s}_{1}},\ g)=A({{s}_{1}},\ g){{\text{e}}^{\text{i}\omega {{t}_{{{s}_{_{1}}}g}}}},$ (1)

其中, A(s₁, g)是振幅项, ω为角频率, ${{t}_{s}}_{_{\text{1}}g}$为从震源s₁至检波点g的走时。

显示原图| 下载原图ZIP| 生成PPT

图1   瑞利波互相关示意图

Fig. 1   Scheme of correlation for Rayleigh waves

对g 和 g′ 处的地震记录做互相关, 得到两点之间的介质响应, 即以 g 为虚源, 在 g′ 接收的地震记录在频率域的表达式如下:

$\varphi (g,\ {g}')=u({{s}_{1}},\ {g}')u{{({{s}_{1}},\ g)}^{*}}=A({{s}_{1}},\ {g}')A({{s}_{1}},g){{\text{e}}^{\text{i}\omega ({{t}_{{{s}_{_{1}}}{g}'}}-{{t}_{{{s}_{_{1}}}g}})}}$

$\ \ \ \ \ \ \ \ \ \ \ \ =A({{s}_{1}},\ {g}')A({{s}_{1}},\ g){{\text{e}}^{\text{i}\omega {{t}_{g{g}'}}}}$ (2)

其中, 上角标*表示复共轭。我们注意到 ${{t}_{s}}_{_{\text{1}}g\prime }-{{t}_{s}}_{_{\text{1}}g}$ =${{t}_{g}}_{{{g}'}},$这是因为瑞利波在自由表面传播, 二者具 有一段相同的射线路径, 在互相关的过程中相互抵消。如图 1(b)和(c)所示, 炮点 s₂ 或 s₃ 在检波点 g和 g′ 的瑞利波记录做互相关得到的介质响应与炮点s₁的相同, 因而可以通过多炮点的互相关记录进行叠加予以增强, 频率域表达式如下:

$\varphi (g,\ {g}')=\sum\nolimits_{j=1}^{N}{{{\varphi }_{j}}}(g,\ {g}')={{\text{e}}^{\text{i}\omega {{t}_{g{g}'}}}}\sum\nolimits_{j=1}^{N}{A}({{s}_{j}},\ {g}')A({{s}_{j}},\ g)\ $(3)

瑞利波通过叠加的方式得到增强, 而其他地震信号(如反射波、随机噪声等)由于不同炮点做互相关后走时差不一样, 在叠加过程中并不能一致地增强。因此, 在提取瑞利波的过程中, 瑞利波得到增强, 其他信号则被抑制。

Dong 等^[9]通过以上步骤提取出从一个检波点至另一个检波点的瑞利波信号, 并利用最小平方滤波法, 从原始数据中减去预测的瑞利波信号, 取得不错的结果。但是, 这一方法将其中一个检波点转换成虚源, 导致为了自适应相减时更准确, 必须要求该检波点附近存在真实炮点, 这在实际勘探地震中不一定能够满足。

超虚干涉法在互相关的基础上增加一次褶积步骤。在这种情况下, 虚拟炮点仅作为中间连接点出现, 对其位置没有严格要求。基于该虚拟炮点位置的不确定性, 可以根据不同的虚拟炮点预测出的瑞利波数据做第二次叠加, 进一步增强瑞利波信号。如图2所示, 从炮点 s′ 至检波点 g′, 中间存在可以作为虚拟炮点的检波点g₁, g₂, g₃等, 任意选择一个检波点g_k (k = 1, 2, 3, …), 将从炮点 s′ 至检波点g′的原始地震记录与检波点g_k至检波点 g′ 的叠加数据做褶积, 可得炮点 s′ 至检波点 g ′的合成瑞利波数据, 对不同中间检波点可做二次叠加。频率域的叠加公式如下:

$\begin{align} & \ \ \ \varphi ({s}',\ {g}') \\ & =\sum\nolimits_{k=1}^{M}{u}({s}',\ {{g}_{k}})\varphi ({{g}_{k}},\ {g}') \\ & =\sum\nolimits_{k=1}^{M}{A}({s}',\ {{g}_{k}}){{\text{e}}^{\text{i}\omega {{t}_{{s}'{{g}_{k}}}}}}{{\text{e}}^{\text{i}\omega {{t}_{{{g}_{k}}{g}'}}}}\sum\nolimits_{j=1}^{N}{A}({{s}_{j}},\ {g}')A({{s}_{j}},\ {{g}_{k}}) \\ \end{align}$

$={{\text{e}}^{\text{i}\omega {{t}_{{s}'{g}'}}}}\sum\nolimits_{k=1}^{M}{\sum\nolimits_{j=1}^{N}{A}}({s}',\ {{g}_{k}})A({{s}_{j}},\ {g}')A({{s}_{j}},\ {{g}_{k}})$。 (4)

显示原图| 下载原图ZIP| 生成PPT

图2   瑞利波褶积示意图

Fig. 2   Scheme of convolution for Rayleigh waves

1.2 L₁多道自适应相减

超虚干涉法预测的瑞利波信号与原始地震记录中的瑞利波通常不匹配, 二者在振幅、相位和到时上不完全一致。因此, 需要设计一个滤波器, 用滤波的方法从原始数据中减去瑞利波信号。这一步骤通常称为自适应相减方法。

最常见的是 L₂范数自适应相减法^[19](即维纳滤波法或最小平方滤波法)以及 L₁ 范数自适应相减 法^[20-21]。在多次波压制方面, 自适应相减法已得到广泛应用, 并且从单道滤波扩展到多道滤波。自适应相减的公式表示如下:

d₀(t) = d(t) - s(t)*a(t) , (5)

其中, d₀(t) (t = 0, 1, 2, …, n-1)为压制瑞利波后的地震数据, d(t) (t = 0, 1, 2, …, n-1)为原始地震数据, s(t) (t = 0, 1, 2, …, n-1)为预测出的瑞利波数据, a(t) (t=0, 1, 2, …, l -1)为滤波器, *表示褶积, n为截取的时窗长度, l为滤波器的长度。

为了方便, 用矩阵乘以向量的形式表示褶积:

d₀ = d - Sa, (6)

其中,

$S\ =\ \left[ \begin{matrix} s\left( 0 \right) & 0 & 0 & \cdots & 0 \\ s\left( 1 \right) & s\left( 0 \right) & 0 & \cdots & 0 \\ s\left( 2 \right) & s\left( 1 \right) & s\left( 0 \right) & \cdots & 0 \\ \vdots & \vdots & \vdots & {} & \vdots \\ s\left( n-1 \right) & s\left( n-2 \right) & s\left( n-3 \right) & \cdots & s\left( n-l \right) \\\end{matrix} \right]$。 (7)

Guitton 等^[20]用一个简单例子证明: 当有效波能量强于干扰噪声(如多次波)时, 用L₂范数自适应相减法会削弱高振幅有效波, 得到的解是能量最小化解, 却不是最优解。Guitton在此基础上提出L₁范数自适应相减法, 将两种信号差的绝对值之和作为目标函数, 对高振幅异常仍能保持良好的稳定性。目标函数表示如下:

e₁(a) = |d - Sa|₁, (8)

其中, |…|₁为对向量求1范数。目标函数e₁(a)为奇异函数, 非连续可导, 常规的方法无法求得使e₁(a)为极小值时的滤波器a。一种改进的方法是Bube等^[22]提出的L₁范数与L₂范数混合迭代的加权最小平方方法, 将式(8)的极值问题转化为加权的L₂范数极值问题。将式(8)转变为如下形式:

${{e}_{1}}(a)=||W(d-Sa)||_{2}^{2}$, (9)

其中, $||\cdot \cdot \cdot ||_{2}^{2}$为对向量求2范数; 加权矩阵

$W=\text{diag}\left( \frac{1}{{{\left( 1+{{\left( {{r}_{i}}/\varepsilon \right)}^{2}} \right)}^{1/4}}} \right)$

(r_i = (d - Sa)_i (i = 0, 1, 2, …, n - 1))

是第 i+1 个采样点进行瑞利波压制后的残余; ε 为先验常数, 一般取$\varepsilon =\frac{\max (d)}{100}$, 事实上ε 的取值将影响L₁范数和L₂范数在计算中的权重。

将式(9)用求和形式展开, 得

${{e}_{1}}(a)=\sum\nolimits_{i=0}^{n-1}{q}({{r}_{i}})=\sum\nolimits_{i=0}^{n-1}{\frac{r_{i}^{2}}{\sqrt{1+{{({{r}_{i}}/\varepsilon )}^{2}}}}}$。 (10)

对任意残余值r_i, 易得

$q({{r}_{i}})\approx \left\{ \begin{array}{*{35}{l}} r_{i}^{2},\ \ \ \ \ 当\left| {{r}_{i}} \right|/\varepsilon \ 较小时, \\ \varepsilon \left| {{r}_{i}} \right|,\ \ 当\left| {{r}_{i}} \right|/\varepsilon \ 较大时。 \\\end{array} \right.$ (11)

因此, 当残余较大时, 滤波方法近似为 L₁范数; 当残余较小时, 近似为 L₂范数。由于残余较大时对应于 L₂ 范数不能很好处理的高振幅异常, 因此 L₁ 范数与 L₂ 范数混合迭代的方法能有效地回避L₂范数的局限性。

由于加权矩阵 W 是滤波因子 a 的函数, 直接求导很困难。可以用迭代的方法求解, 每一步迭代都假设 W 是一个常量。此时, 对式(9)求偏导并令导数为零, 得

S^TW^TWSa = S^TW^TWSd。 (12)

对于多道L₁自适应相减法, 则需要对多道进行求和运算, 式(12)改为

$\sum\nolimits_{i=1}^{k}{S_{i}^{\text{T}}}W_{i}^{\text{T}}{{W}_{i}}{{S}_{i}}a=\sum\nolimits_{i=1}^{k}{S_{i}^{\text{T}}}W_{i}^{\text{T}}{{W}_{i}}{{d}_{i}}$。 (13)

设加权矩阵 W 初始值为单位矩阵, 代入式(13)得L₂范数滤波因子 a, 将 a 代入 r=d-Sa 得到滤波残留 r, 由 r 又可得到一个新的加权矩阵 W。实际操作中, 可以设定最大迭代次数和阈值, 当前后两次滤波残留 r 的 L₁ 范数差值小于阀值时, 可得到收敛的r, 即为L₁范数自适应相减结果。

2 实例分析

2.1 简单层状模型

首先对简单层状模型进行处理, 以验证方法的有效性。模型参数见表 1, 地表类型为起伏地表。

图 3(a)为某一炮的共炮点道集, 可以看到两条两边延伸的能量远大于反射波的瑞利波, 瑞利波与有效波有部分区域重叠。图 3(b)为该道集用超虚干涉法瑞利波压制后的结果, 可以看到瑞利波能量大部分被消除了, 在近偏移距处还能看到少许瑞利波残留, 远偏移距压制效果更好。在瑞利波与有效波交叉的地方, 有效波信号得到较好的恢复。实际上, 偏移距最小的 10 道未做瑞利波压制处理, 因为近偏移距瑞利波信号会与直达波耦合在一起, 造成瑞利波预测不准确。对于这一段, 由于没有反射波的存在, 后续处理时直接切除便可。

显示原图| 下载原图ZIP| 生成PPT

图3   水平层状模型瑞利波压制结果

Fig. 3   Suppression of Rayleigh waves for the horizental layered model

为了进一步展示瑞利波压制效果, 经过基本的地震数据处理步骤, 包括球面扩散补偿、直达波切除、动校正及叠加, 得到时间叠加剖面。在叠加剖面上, 超虚干涉法也基本上消除了瑞利波的影响。图 4(a)显示瑞利波压制前的剖面图, 图 4(b)是超虚干涉法瑞利波压制后的剖面图(图 4 中CMP 为共中心点道集)。可以看到, 瑞利波基本上被消除, 层位较清晰。

显示原图| 下载原图ZIP| 生成PPT

图4   瑞利波压制前后剖面对比

Fig. 4   Comparison of stacked sections for surface wave suppression

2.2 复杂模型

为了检测地下结构更复杂情况下的瑞利波压制效果, 设置S波速度模型(图5)。图5中, 地表类型为起伏地表, 浅层存在 S 波速度横向变化, P 波无变化, 在地下加入倾斜地层。

显示原图| 下载原图ZIP| 生成PPT

图5   存在表层S波速度横向变化的复杂模型

Fig. 5   S-wave model with lateral velocity variation

由于浅层 S 波速度存在横向变化, 瑞利波速度也随之存在横向变化。图 6(a)显示原始共炮点道集, 未做增益时反射波几乎不可见。图 6(b)为干涉法瑞利波压制后的结果, 可以看到瑞利波能量大部分被消除, 但残余的瑞利波能量仍然很大, 其量级与反射波能量相近, 反射波同相轴与瑞利波交叠的部分仍然不能清晰地分辨。

显示原图| 下载原图ZIP| 生成PPT

图6   复杂模型的瑞利波压制结果

Fig. 6   Suppression of Rayleigh waves for the complex model

为了进一步对比干涉法与其他一些主流压制瑞利波方法, 图 7 分别显示未做瑞利波压制、局域滤波(在瑞利波出现的区域内对数据的低频部分进行处理, 达到压制瑞利波的目的)、f-k 滤波和超虚干涉法压制瑞利波的叠加剖面。图 7(a)为未做瑞利波压制的原始叠加剖面, 瑞利波清晰可见。图 7(b)用局域滤波方法压制瑞利波, 可以发现瑞利波压制效果并不好, 主要原因是该合成数据中瑞利波与有效波频率较为接近, 仅根据频率的差别很难有效区分瑞利波与有效波, 压制瑞利波会损伤有效波信号(图 7(b)中也反映出这一点)。图 7(c)是用f-k方法进行瑞利波压制, 可以发现压制瑞利波较为彻底, 但有效波也受到一定的损伤。从图 7(c)中可以发现多处有效波同相轴形态发生变化, 尤其是在瑞利波与有效波交叠的部分, 其原因也是瑞利波与有效波分离不彻底, 对瑞利波滤波压制的同时也对有效波的部分信号做了同样处理, 导致部分有效波失真。图7(d)是用超虚干涉法做的瑞利波压制, 瑞利波残留以及对有效波的损伤都比局域滤波和 f-k 滤波法小, 且不会产生有效波同相轴的形态变化。

显示原图| 下载原图ZIP| 生成PPT

图7   不同方法压制瑞利波的叠加剖面对比

Fig. 7   Comparison of stacked sections by different Rayleigh wave suppression methods

进一步, 我们通过信噪比分析, 定量对比这 3种方法的效果。由于已知真实层位的位置, 可以确定叠加剖面中有效波的位置, 因此可以大致求得有效波的平均能量。信噪比公式:

$\text{SNR}=10\cdot {{\log }_{10}}\left( \frac{\text{PS}}{\text{PN}} \right)$, (14)

其中, PS为有效波的平均能量, PN为瑞利波残留的平均能量。由于深部的剖面包含多次波噪声, 为了只体现瑞利波对信噪比的影响, 仅选取图 7 中时间小于60 ms的剖面进行信噪比求取。从表 2 可以看到, 超虚干涉法瑞利波压制后信噪比最高, 残留的瑞利波能量最小。

3 结论

超虚干涉法压制瑞利波主要包括瑞利波预测和自适应相减两个步骤。本文在互相关预测瑞利波的基础上增加褶积步骤, 直接提取出炮点至检波点的瑞利波。从合成数据处理可以发现, 瑞利波能较好地从原始数据里预测出来。通过多道 L₁ 范数自适应相减办法, 可以较为有效地从原始数据里减去瑞利波。由于超虚干涉法具有完全数据驱动的特点, 不依赖于瑞利波与有效波频率和视速度的差异, 对于复杂地表有较好的适应性, 与常规的 f-k 方法相比, 在某些情况下瑞利波压制效果更好。

超虚干涉法压制瑞利波仍然有很多不完善的地方, 技术上还不成熟。基于自适应相减的方法会有压制不彻底的可能性。对于信噪比非常低的区域, 过高的瑞利波压制残留会干扰反射波信号, 使得地震资料的解释达不到预期值。可以从以下两方面加以改进: 1) 改进瑞利波的预测方法, 使得预测出的瑞利波与原始数据的瑞利波有更高的相似度, 降低自适应相减的难度; 2) 改进自适应相减算法, 使得在尽可能压制瑞利波的同时, 尽量不损伤有效波。进一步的工作有可能联合其他瑞利波压制方法共同处理, 并尝试该方法在实际数据上的应用。

The authors have declared that no competing interests exist.

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子