2013, Vol. 56
2. 国土资源部青岛海洋地质研究所, 青岛 266071
2. Qingdao Institute of Marine Geology, Ministry of Land Resources, Qingdao 266071, China
多次波作为一种主要的叠前地震噪音, 其存在难以识别, 影响地震成像的真实性和可靠性, 使目的层反射波的振幅、频率和相位发生畸变, 误导地震资料的解释, 断裂的描述和反演的研究, 增加了目标层位圈定和构造解释的不确定性[1-4].目前的偏移成像算法大多是假定地震记录中仅存在一次反射, 因而多次波的压制显得尤为重要.多次波分为表面多次波和层间多次波两种, 由于前者的能量比较强, 研究人员将大部分的努力投入到它的识别和衰减上.但是, 随着勘探程度的增大, 层间多次波问题逐渐突出, 由于传统方法的局限性, 其压制更加困难.
通常人们把压制多次波的方法分为两大类:一类是基于信号分析处理的滤波方法, 如τ-p变换, f-k变换和预测反褶积等; 另一类是基于波动方程的预测减去法[5], 如波场外推法, 反馈法和逆散射级数方法等.第一类各方法都存在一定的限制性:多次波和一次波有相似的动校正速度减弱了τ-p变换的有效性, 多次波和一次波具有相似的频率和振幅, 减弱了预测滤波的有效性; 多次波和一次波具有相似的斜率, 减弱了f-k滤波的有效性[6].第二类方法主要是利用多次波的产生机理来进行预测和压制多次波, 它综合考虑了多次波传播的运动学和动力学特征, 能适应复杂的地下结构, 利用本身的数据作为多次波预测算子, 需要较少或不需要关于地下结构的假定, 在不损坏一次反射信号的情况下衰减多次波.
对于层间多次波, 目前有效的是预测减去方法, 其中数据驱动型的方法由于需要较少的先验信息, 更具优势[7].El-Emam(2005)[8]给出了一维的利用偏移后数据进行层间多次波建模的方法, 取得一定的压制效果.Berkhout(2005)[9]提出基于共聚焦点技术的层间多次波预测方法, 它有效地将表面多次波压制技术转移到了地下散射点, 从而实现层间多次波的预测, 但是该方法依赖于速度模型, 并且需要复杂的基准面重建.Weglein(1997, 2003)[10-11]提出基于点散射模型的逆散射级数层间多次波预测方法, 能够一次预测出给定阶次的所有层间多次波, 同时较为准确地提供振幅信息.但是它要求有常速度背景, 以保证方法的收敛性, 并且要求满足垂直走时单调性假设条件, 计算量非常大.
针对上述方法在实现时需要满足各项假设的限制以及复杂的计算工作量, Ikelle(2006, 2009)[12-13]提出了通过构建虚地震同相轴估计层间多次波的方法.本文基于该思想准确地预测了层间多次波, 并在此基础上通过L1范数匹配算法实现多次波的压制.该方法摆脱了其他方法对特定界面或初始速度的依赖, 大大减少了计算量, 提高了层间多次波压制的效率.
2 方法原理 2.1 虚同相轴的概念在地震波成像理论中有两种类型的波:延迟波和超前波, 它们分别在正时间和负时间传播.延迟波随时间向前传播, 它与地震波接收的路径一致, 即在离开震源一定时间后到达检波点.超前波相反, 它随时间向后传播, 也就是说在离开震源之前已经到达了检波点.
在散射图中(如图 1), 地震波的传播从左侧开始, 在右侧结束.实线表示随时间正向传播, 虚线表示随时间反向传播.延迟波从左侧开始, 右侧结束, 超前波反之.
|
图 1 表面多次波和层间多次波的构建示意图 (a)表面多次波；(b)虚地震同相轴；(c)层间多次波. Fig. 1 Scattering diagram of surface multiple and IM (internal multiple) construction (a) Surface multiple; (b) Virtual seismic event; (c) IM. |
箭头表示了传播的方向.两条线相交的点就是所谓的散射点.真同相轴代表真实的地震记录, 为正向传播; 而那些对应于向后传播的非因果信号被称作非因果同相轴.混合了上述两种信号的同相轴并不能被实际记录到, 称之为虚地震同相轴, 它只出现在构建真同相轴的中间阶段.
“虚同相轴”, 顾名思义, 它是不能在地震记录中被直接记录到的, 但是它的存在却可以使我们成功地利用表面的散射点来构建层间多次波, 实现思路类似于基于波动原理构建表面多次波(图 1a)的方法, 即将表面的一个散射点作为向下反射点, 利用以该点作为接收点和以该点作为震源点的两个一次波的褶积来实现多次波的预测.
2.2 层间多次波的预测预测层间多次波分两步:(1)构建虚同相轴(图 1b), 通过对两个一次波求互相关实现; (2)预测层间多次波(图 1c), 通过虚同相轴和一次波褶积得到.由此我们可以清晰地看出利用虚同相轴预测层间多次波和其他基于波动理论预测层间多次波方法的区别:这里巧妙地把向下反射点从地下层界面移到了表面, 利用两个形式上类似于反射波褶积的工作来完成多次波的预测; 其他方法都是通过设法寻找地下的散射点以实现预测.然而由于地下的地质情况非常复杂, 很难准确定位向下反射的散射点, 影响了层间多次波预测的精确性, 同时使工作非常复杂.
3 数值实现首先, 在构建虚同相轴时, 我们要将数据分成上下两部分, 因为这一步是两组数据进行互相关, 如果数据中包含相同的同相轴, 则直接做相关后结果中会出现类似于直达波的记录, 进而导致第二步计算中预测出一次波, 使结果不准确.
通过散射图(如图 2)我们可以看到将数据分成上下两部分仍然能够预测出所有的层间多次波.图 2a等号左侧代表了一个四层的模型产生的四个一次反射波.蓝线代表分界面, 把数据分成两部分.将数据进行互相关(图 2a)得到虚同相轴, 然后将虚同相轴和下半部分的数据进行褶积可以预测出一系列层间多次波(图 2b), 这些层间多次波包含了分界面之上层界面(即第一个层界面)产生的所有一阶层间多次波.
|
图 2 层间多次波预测的散射示意图 (a)构建虚地震同相轴；（b)预测蓝色分界面以上层界面产生的层间多次波. Fig. 2 Scattering diagram of IM prediction with a four layered model (a) Construct virtual events; (b) Predict IMs generated by the layers above blue artificial interface. |
当然, 由于分界面以下还包含着各个阶次的层间多次波, 因而这个过程可以预测出所有阶次的上层界面产生的层间多次波.
对于二维情况, 虚同相轴和层间多次波在频率域的构建公式如下:
|
(1) |
|
(2) |
其中dV(xs, w, xr)表示构建的虚同相轴, dI(xs, w, xr)表示构建的层间多次波, d0(x, ω, xr)和d′0(x, ω, xr)分别表示原始数据的上下两部分, d0*(x, w, xr)表示d0的共轭, xs表示震源坐标, xr表示接收点坐标.
我们通过具体的数据来显示虚同相轴构建、层间多次波预测和压制的效果.如图 3所示, 将含四层反射界面的速度模型形成的24道地震记录d(x, t)分成两部分:d0(x, t)和d′0(x, t).通过计算得到了虚同相轴(如图 3d)和预测的层间多次波(如图 3e).这里我们在一二层之间选分界面.由图 3e可见所有与第一层相关的层间多次波都得到了准确的预测.经过匹配压制得到了减去层间多次波之后的结果(如图 3f), 可见效果比较理想.这里的匹配算法采用的是滑动窗口均衡多道L1范数方法, 后续会有详细的介绍.
|
图 3 层间多次波的预测和压制流程示意图 (a)合成记录；（b)第一次循环数据的上半部分；(c)第一次循环数据的下半部分；(d)第一次循环虚同相轴；(e)第一次循环预测的层间多次波；(f)第一次循环多次波压制之后的记录；(g)第二次循环虚同相轴；(h)第二次循环预测的层间多次波；(i)第二次循环多次波压制之后的记录. Fig. 3 The flow of IM prediction and subtraction (a) Synthetic record; (b) Upper part data in first loop; (c) Lower part data in first loop; (d) Virtual events in first loop; (e) Predicted IMs in first loop; (f) After IM subtraction in first loop; (g) Virtual events in second loop; (h) Predicted IMs in second loop; (i) After IM subtraction in second loop. |
为突显效果, 将图 3e对应的预测多次波放大, 提取1.4~2.4s的记录(如图 4), 将各个多次波所对应的散射示意图列在周边, 可见第一层产生的具有有效反射能量的层间多次波都得到了比较好的预测.
|
图 4 第一次循环预测出的层间多次波在1.4~2.4 s的记录与散射图对应情况 Fig. 4 The corresponding scattering diagrams to the predicted IMs in 1.4~2.4 s for the first loop |
同理, 我们继续选取新的分界面, 这次选在第二和第三个反射层之间, 即可消除第二层产生的所有层间多次波, 如图 3(g-i)所示, 得到了虚地震同相轴、预测的层间多次波和消除层间多次波之后的地震记录.由于第一层和第二层是两个主要的反射界面, 因而这里我们只针对它们进行计算.对比图 3a原始数据和图 3i多次波压制之后的结果可以看到效果已经非常明显.
4 滑动窗口多道L1范数多次波匹配算法从公式(2)可以看出, 预测出来的层间多次波相较于实际数据中的多次波在振幅和频率上存在偏差, 另外在时间上也存在或提前或延迟的一个短的移动, 因此需要通过一定的算法与实际数据进行匹配, 以实现最终的压制.
目前最常用且最易于实现的方法就是基于能量最小的最小二乘匹配算法和基于长度最小的L1范数匹配算法, 为了改善效果, 人们尝试单道和多道、时间域和频率域等多种方法.最小二乘算法根据多次波压制以后一次波的能量最小, 即L2范数最小实现匹配, 基于两个假设条件:一次波与多次波具有正交性, 多次波压制后的地震记录能量最小.尽管改善后的算法一定程度上减弱了正交性的要求, 但是仍然无法避免能量最小假设.基于L1范数最小的匹配方法首先由Guitton[14]提出, 它是在单道实现, 对相对大的异常保持稳健但是仍要求多次波和一次波有正交性.为此李学聪[15]提出了均衡多道L1范数自适应匹配滤波方法, 有效地避免了最小二乘方法的两个假设.
本文利用滑动窗口均衡多道L1范数匹配滤波方法, 通过迭代收敛算法[16]有效地压制了预测的层间多次波, 减去多次波后的一次波记录为:
|
(3) |
其中, d1(t), d2(t), d3(t), …dk(t)是k道地震记录, m1(t), m2(t), m3(t), …, mk(t)是与之对应的多次波模型记录, a(t)是自适应滤波器, d01(t), d02(t), d03(t), …, d0k(t)是滤波后的一次波输出.
将(3)式写成向量乘积形式为
|
(4) |
其中
|
(5) |
自适应滤波器a(t)的求取是通过最小化目标函数实现:
|
(6) |
由于目标函数为奇异函数, Bube[17]提出了一种L1/L2范数混合迭代的加权最小平方方法, 将公式(6)的L1范数最小化问题转化为加权的L2范数最小化问题:
|
(7) |
其中, Wi为加权矩阵:Wi=
在公式(7)两边对算子a求偏导并令其等于零, 即可得到公式(7)的最小化, 进而转化为线性方程的求解:
|
(8) |
我们通过一定的线性迭代算法即可求出所需的滤波算子.
通过滑动窗口, 即在时间方向选取一个窗口, 每次对相应的窗口进行计算; 同样, 在偏移距方向, 每次选择一定数量的地震道进行计算.相邻的时间窗口和道窗口都有重叠以保证结果的连续性和平滑性.
下面通过图 3对应数据的预测多次波和匹配层间多次波的对比来验证L1范数方法的效果.如图 5所示, (a)表示预测的层间多次波, 因为经过了多次的相乘和相加, 振幅能量相对强, 应用多道L1范数最小算法匹配后, 得到了真实的多次波(b).对比图 3a的原始数据和图 3f的多次波压制之后的记录亦可证实匹配的有效性.
|
图 5 L1范数匹配效果 (a)The predicted IM; (b)The multiple after match. Fig. 5 Match result by using L1 norm method (a) The predicted IM; (b) The multiple after match. |
我们采用五层的水平速度模型(如图 6), 表面设为吸收边界, 则合成地震记录中不包含表面多次波.地震记录参数包括120炮, 每炮120道, 炮间距20m, 道间距20m, 采用单边放炮, 添加随机噪音.
|
图 6 五层速度模型 Fig. 6 Five layered horizontal velocity model |
合成的单炮记录如图 7a所示, 层间多次波和一次有效波同时出现在地震记录中, 严重干扰了有效信号.依次在一二层之间和三四层之间设置分界面, 进行虚同相轴的构建和层间多次波的预测、匹配及压制, 所得结果如图 7b和7c所示, 多次波能量得到了有效的压制.为了进一步显示压制的效果, 这里依次给出了原始的和处理后的叠加剖面(如图 7(d-f)).通过对比可以发现多次波的能量除了部分层位由于匹配的原因有微小的残余之外, 大部分已经得到了很好的压制, 同时有效地保留了一次波的能量.
|
图 7 多次波压制前后效果对比 (a)原始地震记录；(b)第一次循环多次波压制结果；(c)第二次循环多次波压制结果；(d)原始叠加剖面；(e)第一次循环后叠加剖面(f)第二次循环后叠加剖面.图中红箭头：一次波；蓝箭头：第一层界面产生的层间多次波；绿箭头：第二层界面产生的层间多次波. Fig. 7 Comparison between the records before and after IM subtraction (a)Synthetic record; (b) Result after first loop; (c) Result after second loop; (d) Raw stack profile; (e) Stack profile after firtt loop; (f) Stack profile after second loop. Red arrow: primary; Blue arrow: IM generated by the first layer; Green arrow: IM generated by the second layer. |
南黄海地区地质条件复杂, 海水浅, 层间多次波非常发育.本文根据该地区的实际地质情况[18-19]建立了12层的水平模型(图 8), 具体参数见表 1.合成记录中仅包含一次波和层间多次波, 由于1, 6, 7, 10, 11界面上下两层的波阻抗差异大, 因而主要的层间多次波由这些界面产生.
|
图 8 南黄海地区水平层位地质模型 Fig. 8 Horizontal geological model of South Yellow Sea |
|
表 1 南黄海地质模型参数列表 Table 1 Parameter list of South Yellow Sea geological model |
利用有限差分法合成地震记录(图 9a), 共200炮, 每炮200道, 炮间距20m, 道间距20m, 采样间隔4 ms, 记录时间3500 ms.从剖面上可以看出1800~2200 ms区间的层间多次波非常发育; 在2200~3400ms层间多次波与一次波混合在一起, 远道可以明显地看出层间多次波的存在; 另外2700ms以下层间多次波和一次波速度近乎相同, 大大减弱了传统方法的可行性.经过本文的方法处理以后(图 9b), 效果明显改善, 大部分的多次波能量得到了很好的压制.这里还给出了1400~2400ms(如图 10)和2400~3400ms的放大图(如图 11), 可见多次波能量得到了很好的压制, 同时一次波能量得到了保留.
|
图 9 多次波压制前后效果对比 (a)原始合成记录；（b)多次波压制之后的记录.其中蓝箭头：一次波；红箭头和圆框：层间多次波. Fig. 9 Comparison between the records before and after IM subtraction (a)Synthetic record; (b)Result after IM subtraction. Blue arrow: primary; Red arrow and Circle: IM. |
|
图 10 1.4~2.4 s的图 9(a, b)放大显示结果 其中蓝箭头:一次波. Fig. 10 Enlargement display in 1.4~2.4 s of fig.9(a) and (b) The blue arrow: primary. |
|
图 11 2.4~3.4 s的图 9(a, b)放大显示结果 其中蓝箭头:一次波. Fig. 11 Enlargement display in 2.4~3.4 s of fig.9(a) and (b) The blue arrow: primary. |
为了进一步显示方法的有效性, 文中给出了多次波压制前后的叠加剖面(图 12(a, b)).对比处理前后叠加剖面, 一次信号得到了显著增强.
|
图 12 多次波压制前后叠加剖面 (a)原始数据叠加剖面；(b)多次波压制后的叠加剖面.其中蓝箭头：一次波；红箭头：层间多次波. Fig. 12 Stack profile before and after IM subtraction Blue arrow: primary; Red arrow: IM. |
Jakubowicz(1998)[20]介绍了一种基于波动方程的层间多次波预测方法, 与本文的方法在原理上有相似之处.该方法利用一个公式实现了多次波的预测, 同时需要保证该公式中共轭部分的数据对应的是一次波的成分, 在具体实现的过程中, 采用了类似于表面多次波压制反馈环方法的迭代原理, 对各层界面对应的层间多次波进行多次预测计算以得到所有阶次的层间多次波.而本文方法进行一次计算即可预测出所有阶次的层间多次波, 无需进行迭代计算, 提高了计算的效率.
本文采用通过构建虚地震同相轴的方法来增加一个虚的地表向下散射点, 从而实现层间多次波的预测.这种方法巧妙地将地下散射点移至地表, 减弱了常规方法进行层间多次波预测的复杂性, 大大减少了程序的运行时间.该方法需要增加一个人工的分界面, 将数据分成上下两部分, 同时上半部分是一次波, 这样每次循环能预测出分界面以上的层界面产生的所有阶次的层间多次波.我们每次向下挪动分界面, 最终实现所有层间多次波的预测.由于产生层间多次波的主要层位一般不多, 因而只要通过少量的循环就能有效地完成预测工作.
预测出的多次波需要采用有效的匹配方法与实际多次波在振幅、频率和相位上进行拟合, 这里采用了滑动窗口多道L1范数最小匹配方法, 它减弱了常规自适应方法对多次波和有效波正交性的要求, 同时便于控制各个参数的大小, 以使预测效果达到更好.对简单模型和南黄海地质模型的数据处理实例验证了本文研究方法的正确性和有效性.
| [1] | Yilmaz O. Seismic Data Analysis:Processing, Inversion, and Interpretation of seismic Data. SEG, 2001. |
| [2] | Sheriff R E, Geldart L P. Exploration Seismology. (2nd ed). Cambridge: Cambridge University Press, 1995 . |
| [3] | Berkhout A J. Seismic Migration:Imaging of Acoustic Energy by Wave Field Extrapolation. New York:Elsevier Science Pub1. Co., Inc., 1982. |
| [4] | 胡天跃. 地震资料叠前去噪技术的现状与未来. 地球物理学进展 , 2002, 17(2): 218–223. Hu T Y. The current situation and future of seismic data prestack noise attenuation techniques. Progress in Geophysics (in Chinese) , 2002, 17(2): 218-223. |
| [5] | Weglein A B. Multiple attenuation:an overview of recent advances and the road ahead (1999). The Leading Edge , 1999, 18(1): 40-44. DOI:10.1190/1.1438150 |
| [6] | Griffiths M, Hembd J, Prigent H. Applications of interbed multiple attenuation. The Leading Edge , 2011, 30(8): 906-912. DOI:10.1190/1.3626498 |
| [7] | Wu Z J, Sonika, Dragoset B, et al. Robust internal multiple prediction algorithm. 81st Annual International Meeting, SEG, Expanded Abstracts, 2011:3541-3545. |
| [8] | El-Emam A, Moore I, Shabrawi A. Interbed multiple prediction and attenuation:Case history from Kuwait. 75th Annual International Meeting, SEG, Expanded Abstracts, 2005:448-451. |
| [9] | Berkhout A J, Verschuur D J. Removal of internal multiples with the common-focus-point(CFP) approach:Part1-Explanation of the theory. Geophysics , 2005, 70(3): 45-60. |
| [10] | Weglein A B, Gasparotto F A, Carvalho P M, et al. An inverse-scattering series method for attenuating multiples in seismic reflection data. Geophysics , 1997, 62(6): 1975-1989. DOI:10.1190/1.1444298 |
| [11] | Weglein A B, Araujo F A, Carvalho P M, et al. Inverse scattering series and seismic exploration. Inverse Problem , 2003, 19(6): R27-R83. DOI:10.1088/0266-5611/19/6/R01 |
| [12] | Ikelle L T. A construct of internal multiples from surface data only:the concept of virtual seismic events. Geophys. J. Int. , 2006, 164(2): 383-393. DOI:10.1111/gji.2006.164.issue-2 |
| [13] | Ikelle L T, Erez I, Yang X J. Scattering diagrams in seismic imaging:More insights into the construction of virtual events and internal multiples. Journal of Applied Geophysics , 2009, 67(2): 150-170. DOI:10.1016/j.jappgeo.2008.10.009 |
| [14] | Guitton A, Verschuur D J. Adaptive subtraction of multiples using the L1-norm. Geophysical Prospecting , 2004, 52(1): 27-38. DOI:10.1046/j.1365-2478.2004.00401.x |
| [15] | 李学聪, 刘伊克, 常旭, 等. 均衡多道1范数匹配多次波衰减的方法与应用研究. 地球物理学报 , 2010, 53(4): 963–973. Li X C, Liu Y K, Chang X, et al. The adaptive subtraction of multiple using the equipoise multichannel L1 norm matching. Chinese J. Geophys. (in Chinese) , 2010, 53(4): 963-973. |
| [16] | 庞廷华.基于匹配滤波器的多次波自适应相减.北京:清华大学, 2009. Pang T H. Adaptive multiple subtraction based on matching filters (in Chinese). Beijing:Tsinghua University, 2009. |
| [17] | Bube K P, Langan R T. Hybrid l1/l2 minimization with applications to tomography. Geophysics , 1997, 62(4): 1183-1195. DOI:10.1190/1.1444219 |
| [18] | 张一波, 刘怀山, 吴志强. 南黄海多次波特征及其速度分析. 海洋地质动态 , 2008, 24(8): 8–13. Zhang Y B, Liu H S, Wu Z Q. Multiple characteristics and velocity analysis in south yellow sea. Marine Geology Letters (in Chinese) , 2008, 24(8): 8-13. |
| [19] | 童思友.南黄海地震资料多次波形成机理及压制技术研究[博士论文].青岛:中国海洋大学, 2010. Tong S Y. The formation mechanism and the suppression technology study of multiple wave of the South Yellow Sea seismic data (in Chinese)[Ph.D.thesis]. Qingdao:Ocean University of China, 2010. |
| [20] | Jakubowicz H. Wave equation prediction and removal of interbed multiples.//68th Annual International Meeting, SEG, Expanded Abstracts, 1998:1527-1530. |