0 引言

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

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

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

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

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

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

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

在下行延拓的过程中，层位n－1下方的下行波场B_d^n－1经过下行延拓至层位n上方，得到A_dⁿ：

(1)

层位n上方的下行波场A_dⁿ经过层位n的透射波场加上层位n下方的上行波场B_uⁿ经层位n的反射波场得到层位n下方的下行波场B_dⁿ：

(2)

在上行延拓的过程中，层位n+1上方的上行波场A_uⁿ⁺¹经过上行延拓至层位n下方，得到B_uⁿ：

(3)

层位n下方的上行波场B_uⁿ经过层位n的透射波场加上层位n上方的下行波场A_dⁿ经层位n的反射波场得到层位n上方的上行波场A_uⁿ：

(4)

式(1)中Γ_d(Z_n－1, Z_n)表示下行延拓矩阵算子，将Z_n－1处的波场延拓至Z_n处；式(3)中Γ_u(Z_n+1, Z_n)表示上行延拓矩阵算子，将Z_n+1处的波场延拓至Z_n处；式(2)中T_dⁿ表示层位n处的下行透射系数矩阵，R_uⁿ表示层位n处的上行反射系数矩阵；式(4)中T_uⁿ表示层位n处的上行透射系数矩阵，R_dⁿ表示层位n处的下行反射系数矩阵.

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

(5)

式中S_d^m表示下行一次源，在震源位于地表的情况下S_d^m=S_d⁰，R_d^mA_d^m+R_u^mB_u^m可理解为在模型中激发的二次源.

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

(6)

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

(7)

(8)

(9)

(10)

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

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

1.2 双重层位约束

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

在公式(2)和(4)中加入双重层位约束矩阵M_dⁿ和M_uⁿ如下：

(11)

(12)

(13)

(14)

M_dⁿ表示下行反射层位约束矩阵，指定产生下行反射的地层，M_uⁿ表示上行反射层位约束矩阵，指定产生上行反射的地层，g_dⁿ和g_uⁿ分别为下行反射约束模型G_d和上行反射约束模型G_u的第n行组成的向量.G_d和G_u的计算方式为：若G_d的第n行第m列位于地质模型中拾取的下行反射地层上，则(G_d)_{n, m}=1，否则(G_d)_{n, m}=0；若G_u的第n行第m列位于地质模型中拾取的上行反射地层上，则(G_u)_{n, m}=1，否则(G_u)_{n, m}=0.

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

(15)

(16)

1.3 二维反周期延拓

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

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

(17)

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

1.4 自适应变步长波场延拓

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

(18)

式中W和w表示空间域波场和波数域波场，e^{iф(x, y, z, k_x, k_y, ω)}表示波场延拓因子.在非均匀介质中，将延拓因子分解为相移项和校正项来处理介质的横向变化：

(19)

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

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

(20)

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

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

2 数值验证

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

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

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

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

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

图 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三组地层之间产生层间多次波的三种不同路径，路径中的反射均表示地层组中的综合反射效应.另外，由于一阶层间多次波占据层间多次波的主要能量，二阶层间多次波的能量与有效能量相比可以忽略不计，本节中的模拟数据总波场均指一次波场和一阶层间多次波场的叠加波场.

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

模拟数据最大偏移距为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道集识别出来，如图 10d和10f中的红色箭头和蓝色箭头所指；3.2~3.4 s的层间多次波动校时差明显大于一次波动校时差，相比之下更易于从原始CMP道集识别出来，如图 10d和10f中的绿色箭头所指.

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

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

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

4 结论

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