刘璐
, 安圣培
Lu LIU
中图分类号:
P315
通讯作者:
Corresponding authors:
收稿日期:
2015-09-20
修回日期:
2015-12-29
网络出版日期:
2016-04-28
版权声明:
2017 北京大学
基金资助:
展开
摘要 针对虚同相轴法在实际应用中存在多次波的高效匹配、其他噪声的干扰以及观测系统的影响等问题, 通过对模型数据的研究, 提出一种基于高阶高分辨率Radon变换的高效匹配方法, 在匹配时间和效率两方面有很大提高。对传统范数匹配方法中的滤波矩阵、时窗长度进行定量分析和优化选取, 对虚同相轴方法的抗噪能力以及观测系统的影响进行定量分析, 并讨论其适用范围。研究结果对该方法的实际应用有一定的参考意义。
关键词：
Abstract There are problems such as efficient match of multiples, influence of other noise and geometry when the virtual events method is used in real data. By studying on model data, an efficient matching method based on high-order and high-resolution Radon transform was presented, which had advantages in matching time and effect. There are also some quantitative analyses and parametric optimization on filter matrix and length of time window of the traditional matching methods. The authors made some quantitative analyses on the anti-noise capacity of virtual events method and discussed its range of application. The results have some referential meaning when the method is used in real data.
Keywords：
多次波的分离是常规地震资料处理过程中的一个重要环节[1], 无论是将其归为噪声, 还是用于成像[2], 多次波与一次波的分离都是必不可少的。传统的基于规则多次波的分离方法主要有两类[3]: 一类是基于有效波与多次波之间速度差异的方法, 如Radon 变换和 F-K 滤波; 另一类是基于波动方程的预测减去法, 如预测反褶积法[4]、波场外推法等。近年来, 由于需要较少的先验信息, 基于数据驱动的方法[5-6]更具优势, 主要有 SRME[7]、逆散射级数法[8]、稀疏反演法[9-11]、共聚焦点技术[12]和虚同相轴法[13-14]。
对于层间多次波而言, 目前比较先进的方法主要有两种。一种是 Weglein等[15]提出的基于点散射模型的逆散射级数法, 根据 Lippmann-Schwinger 方程分别构造对应于表面多次波和层间多次波的子级数来计算多次波。这是一种基于波动方程的数据驱动方法, 不需要先验的速度结构模型, 对于复杂构造有较好的适应能力, 但是计算量非常大。另一种是 Berkhout[16-17]提出的基于共聚焦点技术的层间多次波预测方法, 是在常规 SRME 方法的基础上发展而来, 其基本原理是利用共聚焦点算子, 将地震数据向下延拓至产生多次波的界面, 从而层间多次波的预测就可以用类似于 SRME 的褶积计算来实现, 用共聚焦点道集记录代替炮集记录。该方法的优势是容易适应复杂的地质条件, 实现起来比较简单, 但需要依赖于初始的速度模型, 以便产生较精确的聚焦算子。
Ikelle[18-19]基于克希霍夫积分表示定理, 提出通过构建虚同相轴来预测层间多次波的方法。虚同相轴是通过形成一个虚的向下散射点, 巧妙地将反射点从地下层界面移到表面, 利用两个形式上类似反射波的同相轴褶积来预测层间多次波。在分层情况比较好的条件下, 可以有效地实现层间多次波的预测。
虚同相轴方法能够较准确地预测层间多次波的到达时间, 但与实际数据中的层间多次波在相位和振幅上有差异, 直接相减会带来误差, 因而需要做匹配。Guitton 等[20]提出 L1 范数匹配法, 主要思想是求取原始地震道与匹配道的 L1 范数最小来计算滤波器, 通过迭代加权最小二乘法求取最小化目标函数, 有效地避免了 L2 范数中最小能量的假设, 同时利用了 L1 范数对于大异常值保持稳健性的特点。Guo 等[21]提出一种利用预测多次波特性来进行匹配的方法, 避免了自适应滤波器的求取, 并且充分利用了预测数据的信息。李学聪等[22]提出均衡多道 L1 范数匹配方法, 通过在空间方向上对地震记录进行均衡, 解决多次波与一次波非正交性问题。Huang 等[23]提出混合 L1/L2 范数匹配法, 通过一个自适应权重因子来调节目标函数中 L1 范数和L2范数的比例。权重因子通过原始数据和预测数据的能量对比和正交性来选取, 采用共轭梯度法求取非线性函数。该方法有效地解决了复杂模型 L1 范数匹配中多次波残留问题, 也能有效地解决有效波与多次波横向相交的问题, 但计算速度较慢。
本文在虚同相轴方法精确预测出层间多次波的情况下, 对 L1 范数和 L2 范数匹配法的参数选择和定义提出理论和实验依据。我们参考变换域滤波的理念, 提出一种基于高阶高分辨率 Radon 变换的快速匹配方法, 实际证明匹配效果很好。同时, 对虚同相轴方法的稳定性和抗噪性做一定的分析, 有助于今后该方法在实际数据中的应用。
虚同相轴方法原理如图 1 所示, 十字星表示炮点, 倒三角表示检波点, 箭头表示位于同一位置, *表示褶积。实线表示地震波从震源出发, 经过一定时间后到达检波点。虚线表示逆时间, 称为非因果同相轴, 可以通过两个一次波相关得到。SRME 方法说明通过褶积两个一次波可以得到表面多次波。同理, 将非因果同相轴与一次波褶积则可以得到层间多次波。
通过对两个一次波求互相关来构建虚同相轴, 相当于将一个一次波求逆时与另一个一次波进行褶积。将虚同相轴与一次波褶积, 就形成层间多次波。由于褶积运算在频率域运算比较简单, 所以对于二维情况, 虚同相轴和层间多次波在频率域的构建公式如下:
${{D}_{v}}({{x}_{r}},\omega ;{{x}_{s}})=\underset{{{x}_{k}}}{\mathop \sum }\,D_{u}^{*}({{x}_{r}},\omega ;{{x}_{k}}){{D}_{d}}({{x}_{k}},\omega ;{{x}_{s}})$, (1)
${{D}_{I}}({{x}_{r}},\omega ;{{x}_{s}})=\underset{{{x}_{k}}}{\mathop \sum }\,{{D}_{d}}({{x}_{r}},\omega ;{{x}_{k}}){{D}_{v}}({{x}_{k}},\omega ;{{x}_{s}})$, (2)
式(1)和(2)中, Dv(xr, ω; xs)表示所构建的虚同相轴, DI(xr, ω; xs)表示预测的层间多次波, Du*(xr, ω; xk)和Dd(xk, ω; xs)分别表示上层数据和下层数据, *表示对数据取共轭。
从式(1)和(2)可以看出,用虚同相轴法预测多次波,对原始数据进行了一次相关和褶积, 使得地震记录中的子波被拉伸, 降低了地层分辨率。因此在进行匹配前, 需要进行反褶积来恢复地震子波的原始形态。
对预测出的层间多次波进行反褶积处理后的信号与原始信号在振幅和相位上存在偏差, 直接相减会造成较大的残留, 因此需要进行匹配。目前比较常用的匹配方法是基于能量最小的 L2 范数匹配法和基于长度最小的 L1 范数匹配法, 前者需要满足多次波与一次波正交, 且需避免匹配时窗内出现大值, 后者可以有效地避免这两个假设条件。相对于单道匹配, 多道匹配效果更佳。
考虑到滤波器正负时移对褶积结果的贡献效果, 本文对文献[22]中的 M 矩阵做了改进, 针对褶积的定义重新排列 M 矩阵, 在多次波与一次波正交时, 能够做到完全匹配。设滤波器长度为l, 默认为奇数, 每个量表示为$a\left( -\frac{l-1}{2} \right)$, $a\left( -\frac{l-1}{2}+1 \right)$, …, a(0), a(1), …, $a\left( \frac{l-1}{2} \right)$, 改进后的 M 矩阵表示如下(原矩阵可参见文献[22]):
${{M}_{i}}=\left[ \begin{matrix} {{m}_{i}}\left( 1+\frac{l-1}{2} \right) & \cdots & {{m}_{i}}\left( 2 \right) & {{m}_{i}}\left( 1 \right) & 0 & 0 & 0 \\ {{m}_{i}}\left( 1+\frac{l-1}{2}+1 \right) & \cdots & {{m}_{i}}\left( 3 \right) & {{m}_{i}}\left( 2 \right) & {{m}_{i}}\left( 1 \right) & \ldots & 0 \\ \vdots & {} & \vdots & \vdots & \vdots & {} & \vdots \\ 0 & \ldots & {{m}_{i}}\left( m \right) & {{m}_{i}}\left( m-1 \right) & {{m}_{i}}\left( m-2 \right) & \cdots & {{m}_{i}}\left( m-\frac{l-1}{2}-1 \right) \\ 0 & 0 & 0 & {{m}_{i}}\left( m \right) & {{m}_{i}}\left( m-1 \right) & \ldots & {{m}_{i}}\left( m-\frac{l-1}{2} \right) \\\end{matrix} \right]$。 (3)
为了防止匹配时因分窗造成信号截断, 采用相邻时窗部分重合, 并且在重合部分做斜坡处理。在设计时窗时, 可根据预测多次波的形态设计双曲线形的时窗, 将预测多次波完整地放在一个匹配时窗中, 截断效应也会大大减弱。
匹配时窗长度的选取一般参考3.5倍子波长度, 具体的效果对比见2.2节。
总体来说, 范数匹配方法计算量大, 需要做精细调整的地方很多, 很难进行自动化。
本文提出一种基于变换域的快速匹配方法, 将原始数据和预测数据变换到另一域。因为一次波和多次波同相轴变换前后均可看成抛物线, 所以在抛物 Radon 空间均可映射成点。从预测数据的变换域中可以看出层间多次波存在的位置, 然后在原始数据变换域中相应位置做去除, 再反变换回时空域。由此, 可对层间多次波进行完全消除。为了防止变换到 Radon 域后, 因存在X状截断效应, 两者重叠在一起难以区分, 使用高阶高分辨率抛物Radon 变换法。高分辨率 Radon 变换法可以使数据收敛, 但在数据信噪比较差的情况下难以得到稳定解, 并使反变换后的数据振幅衰减或畸变。高阶变换(即正交多项式变换)用正交多项式表征同相轴AVO 特性[24], 可以很好地保存同相轴的振幅特性。将 Radon 变换的方向特性与正交多项式的振幅描述优势相结合, 可以得到高阶高分辨率抛物Radon变换[25]:
${{m}_{j}}(\tau ,\ q)=\underset{i=0}{\overset{N}{\mathop \sum }}\,d(t=\tau +qx_{i}^{2},{{x}_{i}}){{p}_{j}}({{x}_{i}}),$ (4)
式中, pj(x)是根据地层模型计算出的正交多项式, mj(τ, q)描述在时刻τ沿曲率q的抛物路径上地震数据振幅的第 j 阶正交多项式分解系数, 可以实现地震数据任意时刻沿不同曲率路径的正交多项式变换。其反变换形式为
$d(t,\ x)=\underset{q}{\mathop \sum }\,\underset{j}{\mathop \sum }\,{{m}_{j}}(\tau =t-q{{x}^{2}},\ q){{p}_{j}}(x),$ (5)
高阶高分辨率抛物Radon变换的实现是先构造代价函数:
$J(m)=||{{W}_{d}}(d-Lm)||_{2}^{2}+\lambda ||{{W}_{m}}m||_{2}^{2},$ (6)
将上式对mT求导, 可得
$m={{({{L}^{\text{T}}}W_{d}^{\text{T}}{{W}_{d}}L-\lambda W_{m}^{\text{T}}{{W}_{m}})}^{-1}}{{L}^{\text{T}}}W_{d}^{\text{T}}{{W}_{d}}d$, (7)
式中, Wd为数据加权对角矩阵, 反映数据的标准方差; Wm是模型加权对角矩阵。
用如图 2 所示的斜层模型进行声波正演, 密度均为 2 g/cm3, 第二层 3000 m/s 的高速斜层易产生层间多次波。模型大小是 900 m×900 m, 采用时间2阶、空间4阶的有限差分法, 表面设为吸收边界, 合成地震记录中不包含表面多次波。地震记录共200 炮, 每炮 200 道, 道间距和炮间距均为 3 m, 采样间隔为 0.00025 s, 记录长度为 1 s。单炮记录如图 3(a)所示, 箭头指向一次波。图 3(b)和(c)分别为预测多次波和匹配相减之后的结果。
图3 斜层模型单炮记录及预测匹配效果 Fig. 3 Synthetic record of oblique layer model and its predicting and matching result 图 4 表示模拟的南黄海地质模型, 密度均为 2 g/cm3。模型大小为 5000 m×2000 m, 声波正演, 采用时间 2 阶、空间 4 阶的有限差分法, 表面设为吸收边界。地震记录共 480 炮, 每炮 480 道, 道间距和炮间距均为 4 m, 采样间隔为 0.00005 s, 记录长度为3.5 s。
图4 模拟南黄海地质模型(据文献[16]修改) Fig. 4 Simulated South Yellow Sea geological model (after Ref. [16])
因层位较多, 若直接用虚同相轴方法, 会在底层产生很多假象, 所以需要计算易产生层间多次波的层位。数据的分割线是与一次波曲率相近的双曲线, 这样可以一次性地将该层的相关层间多次波全部预测出来, 然后进行匹配和消除, 接着再预测其他层位的相关多次波, 避免残留而造成混淆。图 5 (a)为单炮记录, 箭头指向一次波。图中显示为 AGC (automatic gain control)增益, 不可避免地将一次波的边界反射信号加强。图 5(b)为用本文方法进行多次波匹配消除的结果。
图 3 和
对预测的层间多次波进行反褶积, 从图 6 可以看出, 反褶积不仅压缩了地震子波的长度, 还将预测多次波记录中无效的微小旁瓣去除掉, 这对之后的匹配处理非常重要。
从图 7 可以看出, 在有效波与多次波正交情况下, 利用改进后的M矩阵, 可以做到精确匹配。
图7 M矩阵改进前后的L2范数单道匹配效果对比 Fig. 7 Single L2 norm matching result before and after improving M martrix
M矩阵改进后, 检验不同时窗长度对匹配效果的影响。首先应保证所匹配段只含多次波, 从而可以用前后能量差来判定匹配的效果。从图 8 可知, 匹配时窗不同, 匹配效果有较大差异。通过其他数据做相同实验, 得出结论: 最佳匹配时窗的长度和子波长度有关, 一般可取 3.5 倍的子波长度。图 8中用于测试的数据, 其子波长度为 321 个采样点, 最佳匹配时窗约为1200个采样点。
随着数据复杂程度增加, 层间多次波同相轴越来越多, 并与有效波同相轴位置近似, 上述匹配方法在实现过程中会出现耗时长、匹配不完全和损伤有效波等问题。下面讨论高阶高分辨率抛物 Radon变换法匹配的结果。
从图9可知, 经过高阶高分辨率抛物Radon变换后, 原始数据和预测数据在变换域中呈现为一个点。根据预测数据在变换域中的到达时间, 在原始数据中减去具有相同到达时间的数据, 再反变换回时空域, 就达到去除多次波的目的。相对于范数匹配法, 该方法省时很多。从图3和图5也可以看出, 匹配效果很好。
图9 原始数据(a)与预测数据(b)的高阶高分辨率抛物Radon变换结果 Fig. 9 Higher-order high resolution parabolic Radon transform of raw data (a) and predicted data (b)
对一个简单水平模型中的单炮记录分别加不同信噪比的高斯白噪、低频背景噪声和低频线性噪声, 检验虚同相轴方法的抗噪能力。
高斯白噪为随机背景噪声, 在实际勘探过程中不可避免, 从图 10 可以看出, 信噪比大于-10 时, 白噪对该方法基本上没有影响, 但信噪比过低, 会造成有效数据完全被掩盖。从图 11 可以看出, 当信噪比低于-15 时, 残差开始成倍增加, 虚同相轴的方法效果开始变差。所以用虚同相轴预测多次波前, 应保证背景噪声信噪比不低于-15。
图10 加不同信噪比的高斯白噪显示与预测效果 Fig. 10 Result of adding different Gauss white noise and its predicted consequence
图11 加不同信噪比白噪后的预测数据与原始预测数据残差对比 Fig. 11 Residual comparison of raw data and the predicted consequence of different white noise added data
当在沼泽、流沙等疏松介质中激发地震波时, 这些介质的固有振动会构成低频背景噪声, 从图12和 13 可知, 存在低频背景噪声时, 在信噪比低于-1时, 虚同相轴法的效果开始变差。
图12 加不同信噪比的低频背景噪声的预测效果 Fig. 12 Result of adding different low-frequence background noise and its predicted consequence
图13 加不同信噪比低频背景噪声后的预测数据与原始预测数据残差对比 Fig. 13 Residual comparison of raw data and the predicted consequence of different low-frequence background noise added data
当震源较浅时, 在地面与空气的分界面附近, 由震源激发可直接产生面波, 其特点是强度大、频率低, 在地震记录上呈线性。图 14 为弹性波正演的含面波及层间多次波记录, 图 15 为含不同信噪比的面波对虚同相轴方法的干扰。根据图 14 和 15可以判定, 当低频线性噪声信噪比低于 16时, 对虚同相轴方法的影响较大。
图14 正演数据含面波情况下的预测情况 Fig. 14 Forward data including surface wave and its predicted consequence
图15 加不同信噪比低频线性噪声后的预测数据与原始预测数据残差对比 Fig. 15 Residual comparison of raw data and the predicted consequence of different low-frequence linear noise added data
从以上分析可知, 各类噪声对虚同相轴的影响程度从小到大依次为: 高斯白噪、低频背景噪声、低频线性噪声。所以, 在对实际数据使用虚同相轴方法预测前, 应优先去除低频线性噪声。
以上是对虚同相轴方法在有噪声情况下的可靠性和稳定性分析。有这些数据做参考, 在用虚同相轴方法处理实际数据时, 就可以对该方法的适用范围有一个很好的把握。
用图 2 的斜层模型, 将放炮间隔分别设为 5道、10 道、20 道, 获得的虚同相轴预测结果和原始预测结果的比较如图 16 所示。可以看出, 当炮间距大于5倍道间距时, 虚同相轴方法的效果开始变差。从克希霍夫积分表示定理来看, 当炮间距与道间距的差距过大时, 用于构建虚同相轴波场的原始波场减少[26], 造成虚同相轴波场的不稳定, 从而导致预测层间多次波波场的不稳定。
图16 虚同相轴方法在不同观测系统下的预测效果 Fig. 16 The predicted consequence of different geometries
本文对构建虚同相轴预测层间多次波方法的实际应用做了几点探讨: 1) 对传统匹配方法中的滤波矩阵进行了改进, 定量分析了时窗长度的最优选取; 2) 提出一种基于变换域的高效匹配方法, 即高阶高分辨率Radon变换, 充分利用原始数据中同相轴的曲率和振幅信息, 避免了范数匹配中长度最小或能量最小的局限, 同时在计算速度上更快; 3) 讨论了虚同相轴方法的稳定性和抗噪性, 指出虚同相轴方法能否正确预测层间多次波, 在很大程度上取决于噪声(特别是面波)是否得到有效压制; 4) 虚同相轴法对观测系统的要求较高, 炮点必须在检波器的位置或其附近, 炮间距不能大于 5 倍道间距。本文分析结果对今后该方法在实际数据中的应用具有一定的参考意义。
The authors have declared that no competing interests exist.
多次波压制的研究现状与进展 地震偏移成像技术研究现状与发展趋势 应用反褶积技术压制多次波的实例分析 44(增刊1): 39‒43 Estimation of multiple scattering by iterative nversion, PartⅠ: theoretical considerations Estimation of multiple scattering by iterative inversion, Part Ⅱ: practical aspects and examples 3D surface-realted multiple elimination using parabolic sparse inversion The first test and evaluation of the inverse scattering series internal multiple attenua-tion algorithm for an attenuating medium // 84th Annual International Meeting, SEG, Expanded Abs-tracts 基于地表数据分离的层间多次波压制 基于稀疏反演的OBS数据多次波压制方法
Latest views of the sparse Radon transform CFP技术层间多次波预测及Curvelet域相减方法 基于构建虚同相轴压制地震层间多次波 基于构建虚同相轴压制地震层间多次波 Using the inverse scattering series to predict the wave-field at depth and the transmitted wave-field without an assumption about the phase of the measured reflection data or back propagation in the overburden Removeal of internal multiples with the common-focus-point (CFP) approach, Part 1: explanation of the theory Removeal of internal multiples with the common-focus-point (CFP) approach, Part 2: Application strategies and data examples A construct of internal multiples form surface data only: the concept of virtual seismic events Scattering diagrams in seismic imaging: More insights into the construction of virtual events and internal multiples Adaptive subtraction of multiples using the L1-norm A new multiple subtraction method using the attributes of predicted multiples // SEG Technical Program Expanded Abs-tracts 均衡多道1范数匹配多次波衰减的方法与应用研究 hybrid L1/L2-norm scheme with an adaptive weight factor for multiple subtraction // SEG Technical Program Expanded Abs-tracts Tracking the amplitude offset by using orthogonal polynomials 高阶高分辨率Radon变换地震数据重建方法
/
, Shengpei AN
1 方法原理
1.1 虚同相轴方法原理
1.2 高效匹配原理
2 数值实现
2.1 虚同相轴法预测多次波
2.2 层间多次波的匹配
2.3 虚同相轴法的抗噪能力
(a1)~(a4) 依次为加信噪比-100, -10, -1高斯白噪的数据及原始数据; (b1)~(b4) 依次为加信噪比-100, -10, -1高斯白噪数据的预测结果以及原始数据的预测结果
(a1)~(a4) 依次为加信噪比-3.77, -1.8, -0.67, 0.66低频背景干扰的显示结果; (b1)~(b4) 依次为加信噪比-3.77, -1.8, -0.67, 0.66低频背景干扰的预测结果
(a) 正演数据, 含面波和层间多次波; (b) 虚同相轴法预测的数据
2.4 观测系统对虚同相轴法的影响
3 结论
参考文献
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]
〈
〉
