随着当前对油气藏的持续勘探开发，海洋油气勘探开发领域正逐步地由常规构造油气藏向非常规构造油气藏扩展，其研究目标日趋复杂，隐蔽性不断增强。只有不断提高地震勘探的采集处理水平才能应对更为复杂的地质条件对油气开发带来的挑战，为了实现勘探目标的高精度、高分辨率成像，特别关键的是要获取信噪比高、宽频带的地震数据。近年来，倾斜拖缆采集(也称斜缆采集)技术[1,2,3]得到了推广应用，其检波器的沉放深度是随偏移距增加而变化的，即从最小偏移距开始检波器沉放深度随着偏移距逐渐增大而增大或减小，使得获取的地震数据中含有陷波频带多样化的地震鬼波。经过鬼波衰减后，斜缆采集数据频谱的低频端与高频端明显宽于常规水平拖缆采集地震数据的频带，已成为海上宽频地震采集中最先进且相对经济可行的技术。

由于斜缆地震勘探主要用一次反射波成像，故所采集的地震数据中的直达波、浅层折射波以及多次波等都应包含在无效信息的界定范围内予以消除。海底(或软海底与基岩的界面)通常是强波阻抗界面，而海洋表面可以看成自由界面，因此在海水层中亦形成阶数多、能量强的全程多次波[4,5,6]。此外，从地下反射界面返回的地震波，在海面发生下行反射后，会在海水层中多次振荡传播形成微屈多次波。由于均含有海底子反射，这两种多次波具有相近的性质，如振幅均较强，且相关谱中的多次波周期基本一致，因此可统一归属到海底多次波(或水底多次波)的范畴。总的来讲，海洋地震数据中往往含有较多能量较强的海底多次波，如果不对其进行压制，则会使得不真实、不准确的成像结果严重影响后续的地震资料解释环节。

随着技术的发展，地震数据处理的首选方法已经转变为基于波动理论的多次波压制方法，其中波场外推法[7]与反馈环法[8,9,10,11]在实际生产中有着最为普遍的应用。这类方法通常包含多次波的预测和衰减两个步骤，即首先基于波动理论预测出地震记录中的多次波成分，然后通过自适应衰减方法去除原始地震记录中的多次波。然而，目前常用的波场外推等算法尚无法直接用来进行三维斜缆地震数据的多次波预测，其主要原因为：首先，基于波场外推等算法的多次波预测其本质上是通过两个子反射的空间褶积来实现多次波预测[12],如果缺少包含某子反射的地震道，将导致无法预测到基于该子反射的多次波信号，而海上三维拖缆地震数据中，联络测线方向的震源间距(即航行线间距)较为稀疏，且受缆数的限制该方向的检波点数极少、炮检距范围也很有限，直接抽取多次波预测必需的地震道进行三维自由界面多次波预测是不现实的[13,14];其次，基于波场外推等算法进行多次波预测时，若包含某子反射的地震道存在旅行时误差，将导致预测到多次波信号与原始记录中的存在显著差异[15,16],而斜缆采集的数据，由于检波器随着炮检距变化有不同的深度(最大深度通常可达50 m),使得地震道中记录的子反射(反射信号到达检波器)与多次波信号中的实际子反射(反射信号到达海面)具有较大时差，从而导致预测的多次波记录存在明显误差。

由上述分析可知，常规的基于波动方程理论的多次波预测算法难以实现三维斜缆地震数据的高精度多次波预测，因此本文提出并实现了基于射线理论的三维海底多次波预测方法，其首先基于叠加速度(或偏移速度)的射线追踪方法，根据炮点深度、检波点深度模拟出海底多次波的时间延迟量。若预测的海底多次波延迟时存在误差，则基于该初始值与给定的时间间隔计算系列海底多次波延迟时，并基于原始数据的相关谱分析确定出与海底多次波周期拟合最佳的延迟时，在此基础上，通过时空域(偏移距域)的褶积过程实现海底多次波预测。基于射线驱动的三维海底多次波预测方法，不需要空间褶积计算，避免了三维地震数据联络测线方向稀疏的影响；同时基于射线追踪与相关谱分析技术可获得高精度的斜缆地震数据其多次波的延迟时，从而实现精确的多次波预测。此外，该方法也同时避免了波场延拓计算，不需对输入数据与预测的多次波模型进行旅行时校正，还具有计算效率高的优势。理论模型试验和对野外数据的处理试验结果证明了该方法能够有效地预测海底多次波，其能在未来广泛应用于斜缆地震数据的多次波压制处理中。

1、基于射线追踪驱动的斜缆地震资料海底多次波预测

对于斜缆采集数据，海底多次波与一次波的时差较小，但其随偏移距的变化量却比较显著，因此建立空变、时变的海底多次波“模型”至关重要[17]。本文提出了基于叠加速度(或偏移速度)的射线追踪方法，其可根据炮点深度、缆深模拟出海底多次波粗略时间延迟量，然后再通过原始数据的相关谱分析确定出与海底多次波周期拟合最佳的延迟时，具体步骤如图1所示。

图1 基于射线追踪驱动的海底多次波预测流程图

1.1 基于水平层状介质模型的多次波旅行时追踪

基于水平层状介质模型的多次波旅行时追踪其目的是基于速度场进行多次波旅行时射线追踪，获得多次波的旅行时信息。由于精确的速度模型难以构建，因此为提高多次波预测的适应性，本文仅采用叠加速度(或偏移速度)曲线建立近似的水平层状介质模型，然后基于该模型求取一个精度不是很高的多次波走时信息，为下一步基于自相关谱的多次波旅行时信息优化提供初始旅行时模型。

选取与炮点位置距离最近的速度曲线，其通常为包含N条速度-时间采样点的曲线，设各点的均方根速度与双程垂直旅行时分别为

输入利用公式建立的水平层状介质模型，根据严格的炮点、检波点坐标与深度进行二分法射线追踪，从而获得与模型中各反射界面相关的一次波与海底多次波传播路径。在射线追踪过程中，基于斯奈尔定律确定射线在介质中的传播方向。与界面i(1≤i≤N)有关的一次波与海底多次波传播路径在界面i上发生上行反射，而其在所经界面j(1≤j<i)上发生下行透射(或上行透射)。根据斯奈尔定律，透射角(或反射角)与入射角满足以下关系[19]:

式中：θj为射线在界面j处的入射角；θ′j为射线在界面j处的透射角(或反射角);vj与vj+1分别代表入射射线和透射射线所在介质的地震波速度。

随着地震波在介质中的传播，其波前曲率将不断变化，地震波的能量与其传播距离成反比，因此需要考虑几何扩散效应对地震波振幅衰减的影响。在水平界面情况下，几何扩散量L为

K为射线总段数，lj为第j段射线的长度，vj表示该射线段所经介质的地震波速度；θm为第m段射线在界面m(1≤m<j)上的入射角，θ′m为第m段射线在界面m上的透射角(或反射角)。

对于含有N层介质的水平层状模型，将分别获得来自于地下各反射界面的N-1个一次波与海底多次波同相轴；其中与反射界面i(1≤i<N)有关的一次波同相轴其传播路径包含K(K=2i)段射线，由于海底多次波又在海水层(第1层介质)中多往返传播了一次，因此其传播路径包含K+2段射线。根据各段射线的长度与所经地层的速度可以计算出总的旅行时。

式中：ti是与界面i有关的一次波旅行时；lj为一次波传播路径中第j段射线的长度；vj表示其第j段射线所经介质的速度；

i是与界面i有关的海底多次波旅行时；

j为海底多次波传播路径第j段射线的长度。

获得与界面i有关的一次波与海底多次波旅行时(见公式(4))之后，可以求出该海底多次波与其对应一次波的时间延迟量Δti,即

式中：i表示模型界面序号；ti是与界面i有关的一次波旅行时，

i是与界面i有关的海底多次波旅行时。

在计算一次波与海底多次波旅行时的同时，将其传播路径中各段射线的长度与所经地层的速度、射线的入射角与透射角(或反射角)代入公式(3)中，可以分别计算出一次波与海底多次波的几何扩散量Li与

式中：Li是与界面i有关的一次波几何扩散量；lj为一次波传播路径中第j段射线的长度；vj表示其第j段射线所经介质的速度；θm为第m段射线在界面m(1≤m<j)上的入射角；θ′m为第m段射线在界面m上的透射角(或反射角);

i是与界面i有关的海底多次波几何扩散量；

i为海底多次波传播路径中第j段射线的长度；vj表示其第j段射线所经介质的速度；

m为其第m段射线在界面m(1≤m<j)上的入射角；

m为其第m段射线在界面m上的透射角(或反射角)。

图1展示了基于射线追踪的一次波与多次波模拟过程，首先在图2(a)所示的偏移速度谱中拾取速度曲线，并利用Dix公式创建水平层状介质模型(见图2(b)),然后根据炮点、检波点的位置及深度进行射线追踪正演模拟，据此求取同一地震道中海底多次波与一次波的时差及振幅变化规律。

1.2 海底多次波延迟时校正与预测

基于射线追踪技术获得的旅行时，由于速度模型以及射线追踪方法的精度原因，导致获得的多次波旅行时信息精度较低。众所周知，对于具有周期性特征的多次波干扰，自相关谱分析是获取其准确延迟时的有效办法[20],因此，本文在基于射线追踪获得低精度多次波旅行时的基础上，再通过自相关谱分析进一步优化多次波旅行时信息的精度。

以追踪的各同相轴旅行时t

为中心，从原始记录d(n,t)中截取出长度为时窗长度l的各道记录，即

式中：ei表示截取的地震记录；n为输入数据的道号；t(0≤t≤l)表示旅行时。

令最大的相关延迟时为m,创建所截取记录

式中：ri表示截取记录的自相关谱；τ为相关延迟时。

图2 基于射线追踪的一次波与多次波模拟示意图

(a)偏移速度谱及拾取的速度曲线；(b)基于水平层状介质模型的一次波与多次波射线追踪；(c)合成的多次波模型。

由于海面的反射系数接近-1,在基于式所创建的自相关谱ri中，一次波与海底多次波因极性相反其相关值为负值，可将谱中的正值充零，此时该谱仅含有一次波与海底多次波的相关能量。以射线追踪所获取的周期(或延迟时)Δti为中心，根据一定时窗长度截取出包含一次波与海底多次波相关能量的区域。在适度平滑的情况下，给定等值线追踪的能量阈值参数，将在谱中追踪出包含了一次波与海底多次波相关能量的区域范围，然后根据该区域内极值出现的时间来确定与二者周期拟合得最佳的延迟时

为了求出目标地震道中各时刻的海底多次波时间延迟量，需要引入与各反射界面i有关的时间延迟量Δt′i,然后以旅行时t为准进行线性插值，即

式中：t表示旅行时，ti是与界面i有关的一次波旅行时；i表示水平层状介质模型中的界面序号，且有1<i≤N。

同理，目标地震道中各时刻的海底多次波与其对应一次波的振幅比为

式中：t表示旅行时，Ai是与界面i有关的海底多次波与一次波振幅比；i表示水平层状界面模型中的界面序号，且有1<i≤N。

由于海面与海底均为强波阻抗界面，从地下反射回来的地震波传播到海面时会发生下行反射，然后经过海底向上反射后被检波器所接收，从而形成了强振幅海底多次波。与一次反射波相比，海底多次波增加了海面位置下行与海底位置上行的传播过程，而且二者在海底及其下部地层中的传播路径也存在一定差异。设海面的反射系数为-1,因传播路径差异导致的海底多次波与一次波振幅比为A(t),则可利用原始地震数据来预测海底多次波数据，时间域的海底多次波预测方程可表示为

式中：D(t)为原始地震数据中的一道记录，t表示旅行时；

为预测的多次波记录，t′表示海底多次波旅行时；Δt为海底多次波的时间延迟量，A(t)表示海底多次波的振幅系数。

2、模型数据的多次波预测与衰减实验

2.1 模型建立与理论地震记录的模拟

为了模拟含有复杂多次波的模型数据，设置一个三维层状介质模型，其x, y, z三个方向长度分别为5 000、5 000与2 000 m, 图3(a)显示了该模型各反射界面起伏形态。图3(b)与3(c)显示了分别沿主测线方向、联络测线方向提取的速度模型切片。在1 500～3 000 m/s的速度范围内将模型设为含有7套匀速地层的介质，其中第5层存在明显的速度翻转；沿主测线方向，模型左右两侧分别发育穹窿和凹陷构造，且穹窿两侧均发育断层构造，其走向均垂直于主测线方向；模型顶部为海水层，其深度范围为100～200 m。

图3 三维理论模型示例

(a)理论模型的反射界面形态；(b)主测线方向的模型切片(y=2 500 m);(c)联络测线方向的模型切片(x=2 500 m)。

将上述模型的海面设置为自由界面，利用主频为30Hz的雷克子波，采用奇偶炮激发的方式并基于声波波动方程的有限差分法实现上述模型的正演模拟，从而获取一组含有多次波干扰的炮集记录(详细观测系统参数见表1)。令总数为44条的奇、偶号航迹线的方向分别为正向与负向，在每条航迹线中设置119炮，则一共获得5 246个炮集记录。图4显示了第1401炮的正演模拟结果，可以看出炮集记录中含有大量强能量、速度与一次波有效信号较为接近的多次波干扰，如箭头指向的多次波同相轴，其对有效信号的提取分析造成了严重的干扰，从而降低了后续地震成像的真实性与可靠性。另外，炮集记录中的中、远偏移距道位置存在振幅较强的直达波与折射波等干扰，这会影响多次波预测与衰减结果，因此利用FK视速度滤波获得了信噪比更高的炮集记录(见图5)。

表1 三维海上地震观测系统参数表

图4 原始炮集记录示例

图5 消除中、远偏移距道中折射波干扰的炮集记录示例

2.2 海底多次波预测分析

输入图5所示消除了折射波干扰的炮集记录，创建叠前时间偏移速度谱(见图6(a)),据此拾取偏移速度曲线，然后利用Dix公式创建近似水平层状介质模型(见图6(b))。由于速度模型海底较为平坦，采用该模型针对所有炮进行海底多次波预测，即基于严格的炮点、检波点位置进行射线追踪，获得与模型中各反射界面相关的一次波与海底多次波传播路径，据此计算海底多次波相对于一次反射的时间延迟量与振幅变化量。

图6 近似水平层状介质模型的建立

输入预处理后的炮集记录(示例见图5),利用海底多次波的延迟时及其与一次波的振幅比值，基于公式(11)实现海底多次波的预测，获得了仅包含多次波成分的地震记录(见图7(b)),同时为便于对比，应用常规的波场延拓技术[21]预测的多次波见图7(a)。

图7 预测的多次波记录

(a)传统方法；(b)本文方法。炮号为1401。

通过对比可知，相比于常规的波场延拓法，应用本文方法预测的多次波其旅行时以及波形信息与输入记录中的多次波更为一致(见箭头指向的位置),这证明了本文中提出的多次波预测方法的精度和可靠性。

2.3 基于曲波域扩展滤波的多次波匹配衰减

针对图4所示原始炮集记录，输入图7所示的多次波预测记录，利用曲波域扩展滤波法进行多次波的匹配衰减[22,23]。图8与9分别给出了衰减多次波后的炮集记录与减去的多次波干扰。通过对比可知，图8所示记录中的多次波得到了明显压制，而减去的多次波干扰中并不存在明显的一次波同相轴(见图9),从而进一步证明了多次波预测结果的有效性。

图8 去除多次波后的炮集记录

(a)基于传统方法预测多次波后；(b)基于本文方法预测多次波后。炮号为1401。

图9 去除的多次波

(a)基于传统方法预测的多次波；(b)基于本文方法预测的多次波。炮号为1401。

对预处理炮集记录(见图5)与本文方法压制多次波后的炮集记录(见图8(b))分别进行基于一次波偏移速度场的叠前时间偏移处理，图10(a)与图10(b)分别显示了多次波剔除前后的炮集记录沿主测线方向的偏移剖面。通过比较可知，多次波剔除前的炮集记录其偏移剖面中含有明显的多次波同相轴(如图10(a)中箭头所示),应用本文方法剔除多次波后的剖面中箭头指向位置的多次波能量被显著衰减，说明通过文中方法的多次波预测与衰减，较好地消除了数据中的多次波干扰，数据的信噪比显著提升。

图10 叠前时间偏移剖面对比

((a)衰减多次波前的偏移剖面；(b)衰减多次波后的偏移剖面。y=2 500 m.)

3、东海浅水区域的多次波预测与衰减分析

3.1 东海某三维工区与原始数据简介

基于东海某三维区块的地震资料进行多次波预测与衰减实验，该工区水深在80～100 m间变化。野外地震勘探采用双船四源的激发方式，通过深度接近线性变化的倾斜缆进行地震信息接收，其检波器深度随着炮检距增大由7 m变化为40 m(见图11),详细的观测参数见表2。受野外勘探方式与地质条件的影响，采集的地震记录中存在较强的混叠干扰、鬼波与多次波，从而显著降低了原始数据的信噪比[24,25,26]。完成混叠干扰衰减、振幅平衡、野值压制、气泡噪音衰减以及鬼波剔除等环节的处理后，针对信噪比较高的地震数据运用本文方法进行多次波预测与剔除试验。

图11 倾斜缆地震采集的检波点深度变化图

为了分析地震数据中多次波的发育特征，抽取单一航迹线中同一震源炮集记录的共偏移距道，基于炮号递增的顺序组成共偏移距剖面(见图12),其中箭头所指位置即为2～4阶海底全程多次波同相轴，其可根据多次波同相轴的周期性与极性特征识别，这些多次波因其具有较强的能量故而将同时段内振幅较弱的一次波信号掩盖。创建图13所示剖面的自相关谱，所得结果见图13,其中箭头指向的同相轴为海底多次波与一次波的相关能量(周期介于120～140 ms),该同相轴振幅强且连续性好，这从另一个角度佐证了原始数据中海底多次波发育强的特征。

表2 三维海上地震观测系统参数表

图12 原始地震记录某航迹线的最小偏移距剖面

图13 如图11所示最小偏移距剖面的自相关谱

3.2 多次波预测分析

输入图14所示预处理后的炮集记录，利用创建的近似水平层状介质模型进行海底多次波的预测。首先基于严格的炮点、检波点位置进行射线追踪，求出海底多次波相对于一次反射的时间延迟量与振幅变化量后，利用公式(11)进行海底多次波预测，获得了如图15(b)所示的多次波记录，基于传统的波场延拓法预测的多次波记录见图15(a)。同样应用曲波域扩展滤波方法，分别基于传统方法预测的多次波记录和本文方法预测的多次波记录进行了多次波处理，剔除结果分别见图16(a)和图16(b)。

图14 原始炮集记录

对比图14与15可知，基于本文方法预测的多次波其多次波同相轴位置和波形都与原始记录中多次波更为一致(多次波位置如图14和15中箭头所示),这证明了本文多次波预测方法的准确性。

图15 预测的多次波记录

(a)传统方法；(b)本文方法。

基于传统波场延拓方法和基于本文方法的多次波预测方法分别预测了7 957炮多次波记录，计算时间如表3所示。

表3 计算时间

由表3可知，本文方法的计算效率约为传统波场延拓法的36倍，因此本文方法相比于常规方法具有更高的计算效率，更适应于大规模实际数据的多次波预测。

3.3 基于曲波域扩展滤波的多次波匹配衰减分析

分别输入图15(a)和15(b)所示的多次波预测记录，利用曲波域扩展滤波法分别对其进行了多次波压制，压制结果分别如图16(a)与16(b)所示。通过对比可知，图16(b)所示记录中的多次波得到了明显压制，而图16(a)中上还有部分多次波残余(如图16(a)中箭头所示),这进一步验证了本文的多次波预测方法相比于传统的波场延拓法更具精度优势。

图16 去除多次波后的炮集记录

(a)基于传统方法预测的多次波；(b)基于本文方法预测的多次波。

图17(a)和17(b)分别显示了基于传统方法预测、压制多次波和基于本文方法预测、压制多次波后的共偏移距剖面，其对应的自相关谱分别如图18(a)和18(b)所示。从自相关谱上也可看出，基于本文方法预测的多次波，经过多次波剔除后其共偏移距剖面的残留多次波明显少于传统多次波预测方法得到的结果。

图17 衰减多次波后的共偏移距剖面

(a)基于传统方法预测的多次波；(b)基于本文方法预测的多次波。

图18 衰减多次波后的共偏移距剖面自相关谱

(a)基于传统方法预测的多次波；(b)基于本文方法预测的多次波。

4、结语

三维斜缆采集技术在宽频地震采集中得到了推广应用，然而其斜缆的特点导致其多次波的周期也发生明显变化，对多次波的准确预测提出更高挑战。针对该问题本文提出了射线追踪驱动的三维斜缆海底多次波预测方法，模型数据试验和对实际数据的处理试验结果均证实了本文提出的方法能够有效地预测斜缆数据的海底多次波，其预测精度明显高于传统的波场外推技术，且具有计算效率高的优势，其在实际斜缆地震数据的多次波压制处理中具有广阔的应用前景。

当然，本文方法也有一定的适应条件，当海底起伏剧烈，尤其是崎岖海底时难以实现多次波的高精度预测，发展新的算法提升崎岖海底地区的多次波预测精度是下一步的工作重点。

参考文献:

[3]张振波,李东方.斜缆宽频地震勘探技术在珠江口盆地的应用[J].石油地球物理勘探,2014,49(3):451-456.

[12]Verschuur D J.地震多次波去除技术的过去､现在和未来(陈浩林等译)[M].北京:石油工业出版社,2010:132-135.

[17]唐进,杨凯,顾汉明,等.海上变深度缆地震采集宽频机理分析[J].地球物理学进展,2015,30(5):2386-2392.

[18]马作彬.使Dix公式求得的层速度稳定的技术[J].国外油气勘探,1989,1(6):82-84.

[19]陆基孟.地震勘探原理[M].东营:中国石油大学出版社,2014:30-31.

[20]陈浩林.地震多次波去除技术的过去､现在和未来[M].北京:石油工业出版社,2010:55-70.

[9]董烈乾,李培明,张奎,等.利用曲波变换预测多次波模型[J].石油地球物理勘探,2015,50(6):1098-1104.

[22]王通,王德利,冯飞,等.三维稀疏反演多次波预测及曲波域匹配相减技术[J].吉林大学学报(地球科学版),2017,47(6):1865-1874.3

基金资助:国家重大科技专项(KJGG2022-0302); 山东省重大科技创新工程项目(2019JZZY010803); 崂山实验室科技创新项目(LSKJ202204803); 国家自然科学基金项目(42074138); 中海石油(中国)有限公司上海分公司技术服务课题(CCL2022SHP007RS)资助~~;

文章来源:李键,秦德文,尹文笋,等.射线追踪驱动的三维斜缆地震资料海底多次波预测[J].中国海洋大学学报(自然科学版),2024,54(06):133-144.