地球物理学报  2020, Vol. 63 Issue (5): 2043-2055   PDF    
基于自适应变步长波场延拓的可控层分阶层间多次波模拟
匡伟康1, 胡天跃1, 段文胜2, 李劲松3, 李艳东3     
1. 北京大学地球与空间科学学院, 北京 100871;
2. 中国石油塔里木油田公司, 新疆库尔勒 841000;
3. 中国石油勘探开发研究院, 北京 100083
摘要:地下低速夹层的存在导致地震数据中包含较强能量的层间多次波,有效识别和预测深部储层上覆地层产生的层间多次波是提高深部储层解释精度的重要环节,而准确模拟层间多次波是辅助识别地震数据中层间多次波的一种非常有效的方法.本文提出了一种基于自适应变步长波场延拓的可控地层分阶层间多次波模拟方法,该方法基于自适应变步长波场延拓,以递归循环的方式实现分阶层间多次波的模拟.通过对模型添加双重层位约束,可以模拟指定地层产生的各阶层间多次波.利用二维反周期延拓方法压制波场延拓的边界反射优于传统方法,例如吸收边界法.提出自适应变步长波场延拓技术,大大提升了波场模拟的效率.理论和数值例子表明,本文方法模拟的一次波和各阶层间多次波与常用的有限差分方法模拟结果具有很好的一致性,且克服了有限差分方法无法分阶模拟波场的不足,显著提升了层间多次波识别的效率.
关键词: 自适应变步长波场延拓      可控层      分阶层间多次波      波场模拟     
Modeling inter-layer multiples based on adaptive step-length-varying wavefield extrapolation
KUANG WeiKang1, HU TianYue1, DUAN WenSheng2, LI JinSong3, LI YanDong3     
1. School of Earth and Space Science, Peking University, Beijing 100871, China;
2. Tarim Oilfield Company, PetroChina, Korla Xinjiang 841000, China;
3. Research Institute of Petroleum Exploration & Development, PetroChina, Beijing 100083, China
Abstract: Seismic data may contain some strong internal multiples resulted from low-velocity layers above the target reservoir. Identification and prediction of such internal multiples in seismic data are important procedures to achieve reliable interpretation. Internal multiples modeling is an effective way to distinguish the internal multiples in seismic data. We propose a method to simulate internal multiples of given layers based on adaptive step-length-varying wavefield extrapolation. Each order of internal multiple is simulated by the recursive loops. With the given layer constraints, internal multiples of varied orders can be simulated. A 2D antiperiodic extension method is implemented to suppress the boundary reflections, which is better than traditional ways, such as the absorbing boundary condition. A novel method of adaptive step-length-varying wavefield extrapolation is proposed to make the simulation more efficiently. Theoretical study and numerical examples show that the results of this method are consistent with the traditional finite difference method. Furthermore, this method can simulate internal multiples of different orders in given layers thus remarkably enhancing the efficiency of multiples identification, which the finite difference method cannot realize.
Keywords: Adaptive step length wavefield extrapolation    Given layers    Different order internal multiples    Wavefield modeling    
0 引言

我国油气勘探深度内普遍存在低速层,例如四川盆地的二叠系和奥陶系的低速泥岩页岩层,塔里木盆地的侏罗系煤层,伊利盆地的侏罗系煤层等.这些低速层的存在导致地震资料中层间多次波发育,严重影响对深部储层的处理和解释.有效识别和处理深部储层上覆地层产生的层间多次波有助于提高深部储层解释精度.

目前对层间多次波的识别方法主要可以分为两类,一类是数据驱动的层间多次波预测方法,另一类是模型驱动的波场模拟方法.数据驱动类的层间多次波预测方法在过去20年取得了很大的进展.逆散射级数方法(Weglein et al., 1997Ma and Weglein, 2015Innanen,2017Sun and Innanen, 2018金德刚等,2008李翔和胡天跃,2009),基于地震散射理论,可以一次性预测所有层的层间多次波主要能量,但该类算法普遍受限于庞大的计算量而难以大范围应用于实际数据中;虚同相轴方法(Ikelle,2006; Ikelle et al., 2009Liu et al., 2018吴静等,2013刘嘉辉等,2018),由拾取的一次波同相轴出发,预测与产生该一次波同相轴的层位相关的层间多次波,虚同相轴方法可以预测比较准确的层间多次波到时且具有较高的计算效率,但该类方法不能预测准确的层间多次波振幅;Marchenko方法(Broggini et al., 2012Wapenaar et al., 2014Meles et al., 2015da Costa Filho et al., 2016da Costa Filho and Curtis, 2018),是一种基于Marchenko自聚焦的层间多次波预测方法,在2015年Meles等首先将Marchenko方法引入到声波层间多次波的预测中,2016年Filho等又将该方法推广至弹性波情况,并在2018年将Marchenko方法应用于北海的实际数据中,此方法首先利用Marchenko自聚焦恢复地下虚源点的上行格林函数波场和下行格林函数波场,然后结合干涉法将其用于构建层间多次波,该方法可以预测某一目标层相关的层间多次波,但这种方法计算量也很大.上述数据驱动类层间多次波预测方法存在一个共同的问题,即无法对层间多次波进行分阶预测,虽然利用聚焦变换可以对表面多次波进行分阶预测(Berkhout and Verschuur, 2003宋鹏等,2015刘伊克等,2018),但是该方法还难以应用于层间多次波的分阶预测.

模型驱动的波场模拟方法在地震勘探领域发挥了巨大的作用,如有限差分方法(Virieux,1984Liu and Sen, 2012Liu,2013王洋等,2015李振春等,2016杨鹏等,2017高静怀等,2018),以有限差分的方式对时空域波动方程进行求解,具有编程简单、计算效率高的特点,广泛应用于地震勘探领域;伪谱法(Kosloff and Baysal, 1982Carcione,2010Li et al., 2018),在傅里叶变换域完成对波动方程空间方向的求导,大大降低了空间差分方向的数值频散,但是相比于时空域有限差分方法,在一定程度上增加了计算量;有限元法(Smith,1975Liu et al., 2015Huang et al., 2017杨顶辉,2002刘少林等,2014张金波等,2018苏波等,2019)跳出有限差分的限制,不再采用传统的矩形网格划分模型,可以用更加灵活的方式实现网格的划分,如三角网格,不规则网格等,从而对复杂界面有更好的刻画,如地下断层,起伏地表等.

上述波场模拟方法均可模拟总波场信号,但是无法将一次波和层间多次波信号分离,更直接地实现对层间多次波的辅助识别.层间多次波直接模拟方法可以实现一次波和层间多次波的分离,且模拟波场的时间和振幅都具有较高的精度.Kennett(1979)提出了一种基于反射率法的层间多次波模拟方法,但该方法目前无法很好地适应介质的横向变化,且无法实现指定地层之间的层间多次波分阶模拟;Covey等(1989)提出了一种基于射线追踪的层间多次波模拟方法,该方法运算效率较高,能处理简单介质中的层间多次波模拟问题,但无法应对复杂变化的模型.近年来,Berkhout等提出了一种基于波场延拓的全波场模拟方法(Berkhout and Verschuur, 2011Berkhout,2014),弥补了有限差分方法只能一次性模拟总波场的缺陷,可以分阶次模拟波场,且能够处理相对复杂的地质模型,Verschuur和Davydenko等将其应用于模拟数据和实际数据全波场成像取得了一定的效果(Verschuur et al., 2016Davydenko and Verschuur, 20172018),但是,该方法目前只能一次性模拟所有地层的各阶层间多次波,不能直接模拟给定地层之间各阶层间多次波,且计算效率有待提升.

本文基于全波场模拟方法,在算法中添加双重层位约束,即指定产生下行反射的地层和产生上行反射的地层,使层间多次波的模拟更加灵活;引入了时间-空间域二维反周期延拓方法处理波场模拟过程中的人工边界反射;提出了自适应变步长波场延拓技术,用以提升层间多次波波场模拟的效率.数值例子表明,本文的层间多次波模拟方法得到的模拟波场与常用的有限差分方法得到的模拟波场具有很好的一致性,且克服了有限差分方法只能模拟总波场的缺点;添加双重层位约束后,可以直接模拟任意两个或多个地层之间的各阶层间多次波,使层间多次波的模拟更加灵活,提升了层间多次波识别的精度;反周期延拓可以得到比传统衰减吸收边界更佳的压制边界反射效果.本文方法模拟得到的指定地层之间的各阶层间多次波可以更好地辅助实际数据中的层间多次波识别,有助于指导制定更加精细的层间多次波处理策略,从而提高深部储层解释的准确度.

1 基本理论 1.1 分阶层间多次波模拟

本文中以层位n表示模型纵向第n层网格.将层间多次波模拟过程中层位n处的波场分为四类,即层位上方的下行波场和上行波场,层位下方的下行波场和上行波场,在频率域分别用AdnAunBdnBun单频列向量表示,向量中的元素对应层位n处不同的水平坐标,A表示层位上方(Above)波场,B表示层位下方(Below)波场,上标n表示层位位置,下标u表示上行波,下标d表示下行波.波场在传播过程中的关系如图 1所示,图 1a表示下行延拓过程中的波场关系,图 1b表示上行延拓过程中的波场关系.

图 1 波场延拓过程中的波场关系 (a)下行延拓过程中的波场关系;(b)上行延拓过程中的波场关系. Fig. 1 Relationship of wavefield in wavefield extrapolation (a) Relationship of wavefield in down-going wavefield extrapolation; (b) Relationship of wavefield in up-going wavefield extrapolation.

在下行延拓的过程中,层位n-1下方的下行波场Bdn-1经过下行延拓至层位n上方,得到Adn

(1)

层位n上方的下行波场Adn经过层位n的透射波场加上层位n下方的上行波场Bun经层位n的反射波场得到层位n下方的下行波场Bdn

(2)

在上行延拓的过程中,层位n+1上方的上行波场Aun+1经过上行延拓至层位n下方,得到Bun

(3)

层位n下方的上行波场Bun经过层位n的透射波场加上层位n上方的下行波场Adn经层位n的反射波场得到层位n上方的上行波场Aun

(4)

式(1)中Γd(Zn-1, Zn)表示下行延拓矩阵算子,将Zn-1处的波场延拓至Zn处;式(3)中Γu(Zn+1, Zn)表示上行延拓矩阵算子,将Zn+1处的波场延拓至Zn处;式(2)中Tdn表示层位n处的下行透射系数矩阵,Run表示层位n处的上行反射系数矩阵;式(4)中Tun表示层位n处的上行透射系数矩阵,Rdn表示层位n处的下行反射系数矩阵.

若层位m处的下行震源向量用Sdm表示,在下行循环的过程中,由(1)式和(2)式循环迭代,可以得到层位n上方的下行波场表达式如下:

(5)

式中Sdm表示下行一次源,在震源位于地表的情况下Sdm=Sd0RdmAdm+RumBum可理解为在模型中激发的二次源.

若层位m处的上行震源向量用Sum表示,在上行循环的过程中,由(3)式和(4)式循环迭代,可以得到层位n下方的上行波场表达式如下:

(6)

式中Sum表示上行一次源,一般情况下,其取值为零向量.公式(5)和公式(6)的推导过程中,包含了如下条件:

(7)

(8)

(9)

(10)

式(7)和式(8)表明反射系数矩阵和透射系数矩阵之间的关系,式中I表示单位对角矩阵,式(9)和式(10)表明延拓算子的累乘传递性质.

层间多次波模拟循环过程如图 2所示,通过一次下行循环和一次上行循环,模拟得到一次波场,分别如图 2a2b所示,图中的圆点表示延拓过程中会激发的二次源;每增加一次下行循环和上行循环可模拟更高一阶的层间多次波,图 2c2d为一阶层间多次波的模拟过程,层间多次波的产生如图中的红色路径所示.若要模拟更高阶次的层间多次波,则需要更多的循环迭代过程.

图 2 层间多次波模拟过程中的递归循环示意图 (a)第一次下行循环;(b)第一次上行循环;(c)第二次下行循环;(d)第二次上行循环. Fig. 2 Recursive loops in internal multiples modeling (a) The first recursive loop; (b) The second recursive loop; (c) The third recursive loop; (d) The fourth recursive loop.
1.2 双重层位约束

1.1节所述方法可以模拟所有地层的一次波及各阶层间多次波,但是所得到的每一阶的层间多次波都包含所有地层产生的层间多次波的总波场,我们无法从得到的层间多次波波场判断产生层间多次波的地层.为了实现更精细的层间多次波模拟识别,本文在分阶层间多次波模拟方法中加入双重层位约束,即指定产生下行反射的地层和产生上行反射的地层,达到模拟任意两个或多个层位之间的任意阶次的层间多次波的目的.

在公式(2)和(4)中加入双重层位约束矩阵MdnMun如下:

(11)

(12)

(13)

(14)

Mdn表示下行反射层位约束矩阵,指定产生下行反射的地层,Mun表示上行反射层位约束矩阵,指定产生上行反射的地层,gdngun分别为下行反射约束模型Gd和上行反射约束模型Gu的第n行组成的向量.GdGu的计算方式为:若Gd的第n行第m列位于地质模型中拾取的下行反射地层上,则(Gd)n, m=1,否则(Gd)n, m=0;若Gu的第n行第m列位于地质模型中拾取的上行反射地层上,则(Gu)n, m=1,否则(Gu)n, m=0.

经过层位约束修正后,同理可以得到层位n上方的下行波场和层位n下方的上行波场表达式如下:

(15)

(16)

1.3 二维反周期延拓

波场延拓算子由波动方程出发,将波场分解为上行波场和下行波场,延拓过程需要将波场变换到频率-波数域完成相应的延拓操作,常见的相移法、广义屏法、傅里叶有限差分法等都涉及到频率-波数域的计算,由于傅里叶变换假设数据是周期性的,在将数据正变换到频率-波数域,在频率-波数域乘上延拓算子,然后再反变换到原始域的过程中会产生人工边界反射,一般采用衰减边界的方法压制这种人工边界反射,但是衰减边界方法只能压制部分的边界反射,在多次延拓的过程中效果更差.本文引入反周期延拓方法(Furumura and Takenaka, 1995),并将其拓展至时间-空间域二维反周期延拓,可以得到比衰减边界更好的压制效果,实际计算中可以根据需要采取不同阶次的反周期延拓,阶次越高压制效果越好,但是相应的计算量也越大.

时间-空间域二维K阶反周期延拓如下式所示:

(17)

上式中s表示原始时间域震源矩阵,s表示反周期延拓后的时间域震源矩阵,k1k2满足0≤k1, k2K,通常情况下取K=2便可得到比较理想的边界压制效果,I表示水平方向空间采样点数,J表示时间采样点数,上式中所有自变量均为整数.将时间域震源矩阵s和反周期延拓震源矩阵s变换到频率域,取其单频向量便得到输入震源向量Sd0和反周期延拓震源向量Sd0.用原始震源Sd0模拟得到原始波场,用反周期延拓的震源Sd0替换原始震源Sd0模拟得到反周期延拓波场,截取两次模拟波场的重叠部分取平均,即得到压制人工边界反射后的模拟波场.

1.4 自适应变步长波场延拓

刘洪等(2006)的研究表明,采用大步长指数变换实现相移加校正的波场延拓方法可以在保证稳定性的前提下提升波场延拓的效率和精度,其核心在于波场延拓算子采用近似展开得到,必然会存在一定的近似误差,大步长延拓可以减少多步小步长延拓过程中的误差迭代,从而在整体上提高波场延拓的精度.相移加校正类波场延拓的一般表达式为:

(18)

式中Ww表示空间域波场和波数域波场,eiф(x, y, z, kx, ky, ω)表示波场延拓因子.在非均匀介质中,将延拓因子分解为相移项和校正项来处理介质的横向变化:

(19)

式中eiα(z, kx, ky, ω)处理背景介质中的延拓,称为相移算子,可在波数域实现,eiβ(x, y, z, kx, ky, ω)为处理介质速度横向变化的校正项,可在空间域实现.

由于大步长波场延拓的优势,在不降低模拟精度的前提下,可以采取自适应变步长波场延拓的方式提高波场模拟的效率,根据地下介质的复杂程度,在构造简单的地方采用大步长,在构造复杂的地方采用小步长.

(20)

式中Cn为模型第n个层位的变步长因子,VnVn-1分别为模型第n个层位和第n-1个层位的速度向量.当Cn=0时,将第n个层位与第n-1个层位合并为一个延拓步长,当Cn≠0时,将第n个层位单独作为一个延拓步长.

对于同一个模型,基于波场延拓的层间多次波模拟的计算效率在一定程度上取决于计算所需的延拓步数,每一步波场延拓的计算量与步长无关.自适应变步长波场延拓技术的使用可以有效减少波场模拟所需要的延拓步数,从而提高整个波场模拟的效率,且局部大步长的使用可以在一定程度上减小小步长波场延拓过程中的误差累积,提高数值模拟的精度.

2 数值验证

为了验证本文方法的正确性,本文首先在一个简单的四层模型上对比本文方法所采用的反周期延拓法和常规的衰减边界压制人工边界反射的效果,并比较本文方法模拟得到的波场和有限差分方法得到的波场的一致性.该四层模型的速度结构如图 3所示,模型水平长度为3000 m,深度为1750 m.

图 3 四层介质速度模型 Fig. 3 Velocity model of four layers

在模拟过程中分别采用常规的衰减边界和二维反周期延拓方法压制人工边界反射,模拟结果如图 4所示.

图 4 压制人工边界反射效果对比 (a)原始未压制边界反射单炮道集;(b)常规衰减边界压制边界反射后道集;(c)二维反周期延拓压制边界反射后道集. Fig. 4 Comparison of suppressing boundary reflections result using different methods (a) Original source gather without boundary reflections suppressing; (b) Source gather after boundary reflections suppressing using traditional absorbing boundary condition; (c) Source gather after boundary reflections suppressing using 2D antiperiodic extension method.

由于波场延拓过程中正反傅里叶变换的周期性假设,导致波场在模型边界发生边界反射并传导到时间轴,使得边界反射严重干扰反射波场,如图 4a中箭头所示,使得模拟结果让人无法接受;常规的衰减边界是在波场模拟的过程中对衰减边界内的波场乘以一个衰减项,以达到衰减边界反射的目的,本文添加衰减边界压制边界反射的效果如图 4b所示,与图 4a相比,衰减边界压制了大部分的边界反射,但是仍然有部分的边界反射残留,如图 4b中箭头所示;而对比图 4b4c可以看到,本文的二维反周期延拓方法比常规的衰减边界方法有更佳的压制边界反射的效果.

为了进一步验证本文方法的正确性,将本文方法模拟得到的波场和有限差分模拟的结果进行比较如图 5所示.

图 5 本文方法模拟结果与有限差分模拟结果的比较 (a)有限差分模拟的总波场;(b)本文方法模拟的总波场;(c)本文方法模拟的一阶层间多次波波场. Fig. 5 Comparison of the results of our method and finite difference method (a) Source gather of finite difference method; (b) Source gather of our method; (c) First-order internal multiples source gather of our method.

对比图 5a5b可以看到,本文方法模拟得到的波场与有限差分方法模拟的波场具有很好的一致性,同时本文方法还克服了有限差分方法无法分阶模拟波场的缺陷,如图 5c所示,本文方法还可以单独模拟某一阶的层间多次波.为了更清晰地比较本文方法和有限差分方法的结果,我们在零偏移距单道上对比本文方法的结果和有限差分方法的结果,比较结果如图 6所示.

图 6 本文方法模拟结果和有限差分模拟结果的比较 (a)本文方法模拟的一次波和有限差分模拟结果的比较;(b)本文方法模拟的一阶层间多次波和有限差分模拟结果的比较;(c)本文方法模拟的二阶层间多次波和有限差分模拟结果的比较.蓝色实线表示有限差分模拟结果,红色实线为本文方法的模拟结果. Fig. 6 Comparison of the results of our method and finite difference method (a) Comparison of the primaries; (b) Comparison of the first-order internal multiples; (c) Comparison of the second-order internal multiples.

图 6中P表示一次波,P后面的数字表示产生一次波的层位,M表示层间多次波,M后面的数字表示发生下行反射的层位,一个数字表示仅发生一次下行反射,即一阶层间多次波,两个数字表示分别在两个数字所指的层位发生下行反射,即二阶层间多次波.从图 6的对比结果可以看到,本文的方法模拟得到的一次波和各阶层间多次波均与有限差分模拟结果很好的一致性,从而验证了本文方法的正确性.

3 塔里木实际资料层间多次波识别

将本文的方法应用到塔里木地区的二维模型上,该模型通过对塔里木地区的实际地震数据和测井数据建模得到的.该区目标层位于深度约5.5~6.5 km的奥陶系碳酸盐岩储层,目标区域上覆石炭系灰岩地层,目标层下伏寒武系基底,石炭系灰岩地层之上发育广泛的侏罗系低速层,侏罗系低速层与石炭系灰岩层之间存在显著的波阻抗差,产生能量较强的层间多次波,对奥陶系碳酸盐岩储层的反射波信号形成严重的干扰,正确识别处理干扰奥陶系碳酸盐储层反射信号的层间多次波能量有助于提升目标层的解释准确度.该区二维速度模型如图 7所示,图 8为塔里木地区2.5 km到6 km深度速度模型放大图,图 8中字母A、B、C标注三组地层,如图 8中虚线所示,数字1、2、3所标示的传播路径表示A、B、C三组地层之间产生层间多次波的三种不同路径,路径中的反射均表示地层组中的综合反射效应.另外,由于一阶层间多次波占据层间多次波的主要能量,二阶层间多次波的能量与有效能量相比可以忽略不计,本节中的模拟数据总波场均指一次波场和一阶层间多次波场的叠加波场.

图 7 塔里木地区二维速度模型 Fig. 7 2D velocity model of Tarim region
图 8 塔里木地区2.5~6 km深度速度模型放大图 Fig. 8 Enlarged view of Tarim region from 2.5 km to 6 km

图 9是对该区实际数据第940个CMP道集进行速度分析后得到的速度谱,从图 9中可以清楚的看到,实际数据中存在明显的层间多次波能量,如图 9中的黑色椭圆所示,准确识别实际数据道集中的这些层间多次波能量,并分析其具体来源对下一步的处理和解释具有重要的实际意义.

图 9 实际数据第940个CMP道集速度谱 Fig. 9 Stack velocity spectrum of CMP 940 of real data

模拟数据最大偏移距为6175 m,道间距50 m,双边采集,模拟炮数为900炮,起始炮水平坐标位置为6175 m,炮间距为50 m,时间记录长度为7 s,时间采样率为0.004 s.实际数据第940个CMP道集和模拟数据相应CMP道集如图 10所示,图 10a—c分别为动校正前实际观测波场的CMP道集、模拟总波场数据CMP道集、模拟一阶层间多次波CMP道集,图 10d—f分别为动校正后实际观测波场的CMP道集、模拟总波场CMP道集、模拟一阶层间多次波CMP道集,从图 10c可以看到,层间多次波能量主要集中在2.8~3.0 s和3.2~3.5 s,如图 10c中的箭头所指;2.8~3.0 s及3.4~3.5 s的层间多次波动校时差与一次波相近,难以从原始CMP道集识别出来,如图 10d10f中的红色箭头和蓝色箭头所指;3.2~3.4 s的层间多次波动校时差明显大于一次波动校时差,相比之下更易于从原始CMP道集识别出来,如图 10d10f中的绿色箭头所指.

图 10 动校正前后实际观测数据和模拟数据第940 CMP道集 (a)动校正前实际观测CMP道集;(b)动校正前模拟总波场CMP道集;(c)动校正前模拟一阶层间多次波CMP道集;(d)动校正后实际观测CMP道集;(e)动校正后模拟总波场CMP道集;(f)动校正后模拟一阶层间多次波CMP道集. Fig. 10 CMP 940 gather of real data and simulated data before and after NMO (a) The CMP gather of real data before NMO; (b) The CMP gather of simulated data before NMO; (c) The CMP gather of simulated first-order internal multiples before NMO; (d) The CMP gather of real data after NMO; (e) The CMP gather of simulated data after NMO; (f) The CMP gather of simulated first-order internal multiples after NMO.

为了分析影响储层反射信号的主要层间多次波能量的来源,我们添加双重层位约束控制波场模拟过程中产生上行反射和下行反射的层位,其结果如图 11所示.由于单一层间多次波同相轴在振幅上显著弱于一次波同相轴,可见的层间多次波往往由多组时间相近的层间多次波相互干涉的形态显现,本文不研究单一地层之间产生的层间多次波,而是研究地层组之间多次反射产生层间多次波,这样更接近于实际情况.图 11a为未添加层位约束时模拟的总波场动校正后的CMP道集,可以看到数据中存在明显的层间多次波能量,如图 11a中的箭头所示;图 11b为添加约束不允许图 8中地层组B发生上行反射和下行反射模拟得到的总波场动校正后的CMP道集,对比图 11a11b可以看到,添加约束后,层间多次波的主要能量基本压制掉了,由此可以判断,数据中的层间多次波主要能量均在地层组B发生了反射,图 11c图 11a11b中CMP道集数据的差异道集.更进一步,我们仅允许地层组B发生上行反射和下行反射得到了图 11d所示的动校正后的一阶层间多次波CMP道集;仅允许地层组A发生下行反射和地层组B发生上行反射,得到了图 11e所示的动校正后的一阶层间多次波CMP道集;仅允许地层组B发生上行反射下行反射,且允许地层组C发生上行反射,得到了图 11f所示的动校正后的一阶层间多次波CMP道集,图 11def中的层间多次波的射线路径分别如图 8中数字1、2、3所标注的路径,在图 11c中它们分别如红色箭头、绿色箭头和蓝色箭头所指.图 11d所示的层间多次波由地层组B内的多次反射产生,图 11f所示的层间多次波由地层组B内的多次反射波经地层组C向上反射后产生,这两组层间多次波可以归类为短程多次波,与一次波动校时差相近,在数据中往往难以识别和处理.图 11e所示的层间多次波由地层组A与地层组B之间的多次反射产生,这组层间多次波可以归类为长程多次波,与一次波动校时差差异明显,在数据中更加容易识别和处理.

图 11 添加层位约束前后模拟数据第940 CMP道集对比 (a)添加层位前模拟的总波场CMP道集;(b)不允许层位组B发生上行反射和下行反射后模拟总波场CMP道集;(c) a和b的CMP差异道集;(d)允许层位组B发生上行反射和下行反射模拟得到的一阶层间多次波CMP道集;(e)允许层位组A发生下行反射,且允许层位组B发生上行反射模拟得到的一阶层间多次波CMP道集;(f)允许层位组B发生上行反射和下行反射,且允许层位组C发生上行反射模拟得到的一阶层间多次波CMP道集. Fig. 11 CMP 940 gather of simulated data with and without reflected-layer constrains (a) Total wavefield of the CMP gather without reflected-layer constrains; (b) Total wavefield of the CMP gather simulated without reflection from layers group B; (c) Difference between (a) and (b); (d) First-order internal multiples simulated with up-going reflections and down-going reflections among layers group B; (e) First-order internal multiples simulated with down-going reflections from layers group A and up-going reflections from layers group B; (f) First-order internal multiples simulated with down-going and up- going reflections from layers group B, and up-going reflections from layers group C.

为了分析层间多次波在叠加剖面上的特征,我们将本文方法模拟得到的层间多次波叠加剖面和实际数据叠加剖面进行对比如图 12所示,图 12a为实际数据叠加剖面,图 12b为本文方法模拟的一阶层间多次波叠加剖面.图中紫色箭头所指为地层组B内的多次反射产生的层间多次波同相轴,如图 8中的路径1所示,绿色箭头所指的层间多次波主要由地层组A和地层组B之间的多次反射产生同相轴,如图 8中的路径2所示,蓝色箭头所指为地层组B和地层组C之间多次反射产生层间多次波同相轴,如图 8中的路径3所示.层间多次波的存在导致叠加剖面的同相轴分辨率降低和连续性变差;图 12a中紫色箭头所指部位,一次波同相轴变宽并带有拖尾现象,分辨率明显下降;图 12a中绿色箭头和蓝色箭头所指同相轴难以追踪,连续性显著变差.可见,为了提高储层解释的准确度,避免假同相轴的干扰,必须进行有效的层间多次波识别和处理.

图 12 实际数据叠加剖面和本文方法模拟的一阶层间多次波数据叠加剖面对比 (a)实际数据叠加剖面;(b)本文方法模拟一阶层间多次波数据叠加剖面. Fig. 12 Comparison of stacked real data and stacked simulated first-order internal multiples (a) Stacked real data; (b) Stacked simulated first-order internal multiples.

另外,通过对比11e、11f的层间多次波与其对应的图 12b叠加剖面层间多次波同相轴可以发现,在道集上路径2所指的层间多次波振幅比路径3所指的层间多次波振幅略强,但是在叠加剖面上路径2所指的层间多次波同相轴能量却明显比路径3所指的层间多次波能量弱,这是由于路径2所指的层间多次波为长程多次波,其动校时差明显大于一次波,动校正无法将其拉平,导致其在叠加的过程中无法同相叠加,而路径3所指的层间多次波为短程多次波,其动校时差与一次波相近,动校正几乎将其拉平,导致其在叠加的过程中能实现同相叠加.由此可见,与长程多次波相比,短程多次波对实际数据的干扰大于长程多次波,且不易于识别和处理,在实际资料的处理中,利用本文的层间多次波模拟方法辅助识别这种短程多次波将有助于制定更加准确和精细的地震资料处理解释策略.

4 结论

地下低速夹层的存在导致地震数据中包含较强能量的层间多次波,有效识别和处理深部储层上覆地层产生的层间多次波有助于提高深部储层解释精度.理论和数值例子表明,本文基于自适应变步长波场延拓的可控层位分阶层间多次波模拟方法的模拟结果与常用的有限差分方法的模拟结果具有很好的一致性,且克服了有限差分方法只能模拟总场,无法对各阶波场分阶次模拟的不足;本文引入的二维反周期延拓方法在层间多次波模拟的过程中很好地压制了人工边界反射,效果优于传统的衰减边界方法;本文通过添加双重层位约束,大大提高了层间多次波模拟的精细度,实现了模拟指定层位之间的各阶层间多次波的目的,显著提升了层间多次波识别的准确度,自适应变步长波场延拓技术的使用,大大提升了层间多次波模拟的效率,使其更加适用于实际地质模型.本文方法在塔里木地区的应用取得了明显的效果,通过添加层位约束,分析了地震数据中主要层间多次波的来源,有助于后续制定更有针对性的层间多次波识别和处理策略,提升深部储层的解释精度.

References
Berkhout A J, Verschuur D J. 2003. Transformation of multiples into primary reflections.//2003 SEG Annual Meeting. Dallas, Texas: SEG, 1925-1928. https://www.researchgate.net/publication/249858473_Transformation_of_multiples_into_primary_reflections
Berkhout A J, Verschuur D J. 2011. Full wavefield migration, utilizing surface and internal multiple scattering.//2011 SEG Annual Meeting. San Antonio, Texas: SEG, 3212-3216. https://www.researchgate.net/publication/268453662_Full_wavefield_migration_utilizing_surface_and_internal_multiple_scattering
Berkhout A J. 2014. Review paper:An outlook on the future of seismic imaging, Part Ⅰ:forward and reverse modelling. Geophysical Prospecting, 62(5): 911-930. DOI:10.1111/1365-2478.12161
Broggini F, Snieder R, Wapenaar K. 2012. Focusing the wavefield inside an unknown 1D medium:Beyond seismic interferometry. Geophysics, 77(5): A25-A28. DOI:10.1190/geo2012-0060.1
Carcione J M. 2010. A generalization of the Fourier pseudospectral method. Geophysics, 75(6): A53-A56. DOI:10.1190/1.3509472
Covey J D, Hron F, Peacock K L. 1989. On the role of partial ray expansion in the computation of ray synthetic seismograms for layered structures. Geophysical Prospecting, 37(8): 907-923. DOI:10.1111/j.1365-2478.1989.tb02240.x
da Costa Filho C A, Meles G A, Curtis A. 2016. Elastic internal multiple prediction using Marchenko and interferometric methods.//2016 SEG International Exposition and Annual Meeting. Dallas, Texas: SEG, 4545-4549.
da Costa Filho C A, Curtis A. 2018. Marchenko and interferometry based multiple attenuation of a North Sea field data set.//2018 SEG International Exposition and Annual Meeting. Anaheim, California, USA: SEG, 4518-4522. https://www.researchgate.net/publication/327536795_Marchenko_and_interferometry_based_multiple_attenuation_of_a_North_Sea_field_dataset
Davydenko M, Verschuur D J. 2017. Full-wavefield migration:using surface and internal multiples in imaging. Geophysical Prospecting, 65(1): 7-21. DOI:10.1111/1365-2478.12360
Davydenko M, Verschuur D J. 2018. Including and using internal multiples in closed-loop imaging-Field data examples. Geophysics, 83(4): R297-R305. DOI:10.1190/geo2017-0533.1
Furumura T, Takenaka H. 1995. A wraparound elimination technique for the pseudospectral wave synthesis using an antiperiodic extension of the wavefield. Geophysics, 60(1): 302-307. DOI:10.1190/1.1443760
Gao J H, Xu W H, Wu B Y, et al. 2018. Trapezoid grid finite difference seismic wavefield simulation with uniform depth sampling interval. Chinese J. Geophys. (in Chinese), 61(8): 3285-3296. DOI:10.6038/cjg2018M0313
Huang X X, Zhao J G, Di B R, et al. 2017. 2.5D frequency-domain finite-element modeling in viscoelastic media using unstructured mesh.//2017 SEG International Exposition and Annual Meeting. Houston, Texas: SEG, 4045-4049.
Ikelle L T. 2006. A construct of internal multiples form surface data only:the concept of virtual seismic events. Geophysical Journal International, 164(2): 383-393. DOI:10.1111/j.1365-246X.2006.02857.x
Ikelle L T, Erez I, Yang X J. 2009. Scattering diagrams in seismic imaging:More insights into the construction of virtual events and internal multiples. Journal of Applied Geophysics, 67(2): 150-170. DOI:10.1016/j.jappgeo.2008.10.009
Innanen K A. 2017. Time-and offset-domain internal multiple prediction with nonstationary parameters. Geophysics, 82(2): V105-V116.
Jin D G, Chang X, Liu Y K. 2008. Algorithm improvement and strategy of internal multiples prediction based on inverse scattering series method. Chinese J. Geophys. (in Chinese), 51(4): 1209-1217.
Kennett B L N. 1979. Theoretical reflection seismograms for elastic media. Geophysical Prospecting, 27(2): 301-321. DOI:10.1111/j.1365-2478.1979.tb00972.x
Kosloff D D, Baysal E. 1982. Forward modeling by a Fourier method. Geophysics, 47(10): 1402-1412. DOI:10.1190/1.1441288
Li J X, Innanen K A, Tao G. 2018. A 3D pseudo-spectral method for SH wave simulation in heterogeneous media.//2018 SEG International Exposition and Annual Meeting. Anaheim, California, USA: SEG, 4005-4009. https://www.researchgate.net/publication/327611850_A_3D_pseudospectral_method_for_SH-wave_simulation_in_heterogeneous_media
Li X, Hu T Y. 2009. Surface-related multiple removal with inverse scattering series method. Chinese J. Geophys. (in Chinese), 52(6): 1633-1640. DOI:10.3969/j.issn.0001-5733.2009.06.026
Li Z C, Yang F S, Wang X D. 2016. High-precision numerical simulation of first-order qP-waves in the transversely isotropic (TI) medium optimized by the LS-RSGFD method. Chinese J. Geopgys. (in Chinese), 59(4): 1477-1490. DOI:10.6038/cjg20160428
Liu H, Yuan J H, Chen J B, et al. 2006. Theory of large-step wavefield depth extrapolation. Chinese J. Geophys. (in Chinese), 49(6): 1779-1793.
Liu J H, Hu T Y, Peng G X, et al. 2018. Removal of internal multiples by iterative construction of virtual primaries. Geophysical Journal International, 215(1): 81-101. DOI:10.1093/gji/ggy270
Liu J H, Hu T Y, Peng G X. 2018. Suppressing seismic inter-bed multiples with the adaptive virtual events method. Chinese J. Geophys. (in Chinese), 61(3): 1196-1210. DOI:10.6038/cjg2018L0663
Liu S L, Li X F, Liu Y S, et al. 2014. Dispersion analysis of triangle-based finite element method for acoustic and elastic wave simulations. Chinese J. Geophys. (in Chinese), 57(8): 2620-2630. DOI:10.6038/cjg20140821
Liu S L, Li X F, Wang W S, et al. 2015. A modified symplectic scheme for seismic wave modeling. Journal of Applied Geophysics, 116: 110-120. DOI:10.1016/j.jappgeo.2015.03.007
Liu Y, Sen M K. 2012. A hybrid absorbing boundary condition for elastic staggered-grid modelling. Geophysical Prospecting, 60(6): 1114-1132. DOI:10.1111/j.1365-2478.2011.01051.x
Liu Y. 2013. Globally optimal finite-difference schemes based on least squares. Geophysics, 78(4): T113-T132. DOI:10.1190/geo2012-0480.1
Liu Y K, Liu X J, Zhang Y B. 2018. Migration of seismic multiple reflections. Chinese J. Geophys. (in Chinese), 61(3): 1025-1037. DOI:10.6038/cjg2018L0368
Ma C, Weglein A B. 2015. A new inverse scattering series (ISS) internal-multiple-attenuation algorithm that predicts the accurate time and approximate amplitude of the first-order internal multiples and addresses spurious events: analysis and tests in 2D.//2015 SEG Annual Meeting. New Orleans, Louisiana: SEG, 4402-4407.
Meles G A, Löer K, Ravasi M, et al. 2015. Internal multiple prediction and removal using Marchenko autofocusing and seismic interferometry. Geophysics, 80(1): A7-A11.
Smith W D. 1975. The application of finite element analysis to body wave propagation problems. Geophysical Journal of the Royal Astronomical Society, 42(2): 747-768.
Song P, Zhu B, Li J S, et al. 2015. Reverse time migration of divided-order multiples. Chinese J. Geophys. (in Chinese), 58(10): 3791-3803. DOI:10.6038/cjg20151029
Su B, Li H L, Liu S L, et al. 2019. Modified symplectic scheme with finite element method for seismic wavefield modeling. Chinese J. Geophys. (in Chinese), 62(4): 1440-1452. DOI:10.6038/cjg2019M0538
Sun J, Innanen K A. 2018. Multidimensional inverse-scattering series internal multiple prediction in the coupled plane-wave domain. Geophysics, 83(2): V73-V82.
Verschuur D J, Staal X R, Berkhout A J. 2016. Joint migration inversion:Simultaneous determination of velocity fields and depth images using all orders of scattering. The Leading Edge, 35(12): 1037-1046. DOI:10.1190/tle35121037.1
Virieux J. 1984. SH-wave propagation in heterogeneous media:Velocity-stress finite difference method. Geophysics, 49(11): 1933-1942. DOI:10.1190/1.1441605
Wang Y, Liu H, Zhang H, et al. 2015. A global optimized implicit staggered-grid finite-difference scheme for elastic wave modeling. Chinese J. Geophys. (in Chinese), 58(7): 2508-2524. DOI:10.6038/cjg20150726
Wapenaar K, Thorbecke J, van der Neut J, et al. 2014. Marchenko imaging. Geophysics, 79(3): WA39-WA57. DOI:10.1190/geo2013-0302.1
Weglein A B, Gasparotto F A, Carvalho P M, et al. 1997. An inverse-scattering series method for attenuating multiples in seismic reflection data. Geophysics, 62(6): 1975-1989. DOI:10.1190/1.1444298
Wu J, Wu Z Q, Hu T Y, et al. 2013. Seismic internal multiple attenuation based on constructing virtual events. Chinese J. Geophys. (in Chinese), 56(3): 985-994. DOI:10.6038/cjg20130326
Yang D H. 2002. Finite element method of the elastic wave equation and wavefield simulation in two-phase anisotropic media. Chinese J. Geophys. (in Chinese), 45(4): 575-583.
Yang P, Li Z C, Gu B L. 2017. Pure quasi-P wave forward modeling method in TI media and its application to RTM. Chinese J. Geophys. (in Chinese), 60(11): 4447-4467. DOI:10.6038/cjg20171130
Zhang J B, Yang D H, He X J, et al. 2018. Discontinuous Galerkin method for solving wave equations in two-phase and viscoelastic media. Chinese J. Geophys. (in Chinese), 61(3): 926-937. DOI:10.6038/cjg2018L0095
高静怀, 徐文豪, 吴帮玉, 等. 2018. 深度均匀采样梯形网格有限差分地震波场模拟方法. 地球物理学报, 61(8): 3285-3296. DOI:10.6038/cjg2018M0313
金德刚, 常旭, 刘伊克. 2008. 逆散射级数法预测层间多次波的算法改进及其策略. 地球物理学报, 51(4): 1209-1217. DOI:10.3321/j.issn:0001-5733.2008.04.032
李翔, 胡天跃. 2009. 逆散射级数法去除自由表面多次波. 地球物理学报, 52(6): 1633-1640. DOI:10.3969/j.issn.0001-5733.2009.06.026
李振春, 杨富森, 王小丹. 2016. 基于LS-RSGFD方法优化的横向各向同性(TI)介质一阶qP波高精度数值模拟. 地球物理学报, 59(4): 1477-1490. DOI:10.6038/cjg20160428
刘洪, 袁江华, 陈景波, 等. 2006. 大步长波场深度延拓的理论. 地球物理学报, 49(6): 1779-1793. DOI:10.3321/j.issn:0001-5733.2006.06.026
刘嘉辉, 胡天跃, 彭更新. 2018. 自适应虚同相轴方法压制地震层间多次波. 地球物理学报, 61(3): 1196-1210. DOI:10.6038/cjg2018L0663
刘少林, 李小凡, 刘有山, 等. 2014. 三角网格有限元法声波与弹性波模拟频散分析. 地球物理学报, 57(8): 2620-2630. DOI:10.6038/cjg20140821
刘伊克, 刘学建, 张延保. 2018. 地震多次波成像. 地球物理学报, 61(3): 1025-1037. DOI:10.6038/cjg2018L0368
宋鹏, 朱博, 李金山, 等. 2015. 多次波分阶逆时偏移成像. 地球物理学报, 58(10): 3791-3803. DOI:10.6038/cjg20151029
苏波, 李怀良, 刘少林, 等. 2019. 修正辛格式有限元法的地震波场模拟. 地球物理学报, 62(4): 1440-1452. DOI:10.6038/cjg2019M0538
王洋, 刘洪, 张衡, 等. 2015. 一种全局优化的隐式交错网格有限差分算法及其在弹性波数值模拟中的应用. 地球物理学报, 58(7): 2508-2524. DOI:10.6038/cjg20150726
吴静, 吴志强, 胡天跃, 等. 2013. 基于构建虚同相轴压制地震层间多次波. 地球物理学报, 56(3): 985-994. DOI:10.6038/cjg20130326
杨顶辉. 2002. 双相各向异性介质中弹性波方程的有限元解法及波场模拟. 地球物理学报, 45(4): 575-583. DOI:10.3321/j.issn:0001-5733.2002.04.015
杨鹏, 李振春, 谷丙洛. 2017. 一种TI介质纯qP波正演方法及其在逆时偏移中的应用. 地球物理学报, 60(11): 4447-4467. DOI:10.6038/cjg20171130
张金波, 杨顶辉, 贺茜君, 等. 2018. 求解双相和黏弹性介质波传播方程的间断有限元方法及其波场模拟. 地球物理学报, 61(3): 926-937. DOI:10.6038/cjg2018L0095