一种基于波动方程边值反演的变深度缆鬼波压制方法与流程

本发明涉及一种能源开发勘探方法，具体是关于一种基于波动方程边值反演的变深度缆鬼波压制方法。

背景技术：

地震勘探是油气勘探的重要手段，特别是海上油气勘探，由于其钻井成本高、测井、地质等其它资料少，使地震勘探在海上油气勘探中具有更重要的作用。宽频带的地震数据可以解释为既能够提高分辨率又能够提高中深层信噪比的基础数据，提高构造解释和储层预测的精度和可靠性，是海上地震勘探追求的目标之一。然而由于种种条件的限制，这一目标并没有很好的实现。在进行海上地震勘探时，一般需要将震源和接收器分别置于海平面以下的一定深度，由于海水和空气的接触面是一个良好的反射面，因此接收器除了接收到来自地下界面的反射波外，还将接收到与海平面的反射作用相关的伴随波，这个伴随波就是鬼波。鬼波会对来自地下的反射信号产生陷波作用，使得地震数据的频带变窄，降低地震剖面的分辨率。当电缆的深度较大时，鬼波导致假的同相轴出现。在海上地震勘探中，拖缆是野外常用的采集作业方式，一般在同一次采集中，将接收缆放置在某一个固定的深度上，接收缆中所有的检波器的深度都相同。为得到频率更高的地震资料，应将拖缆放置在尽量浅的深度上，但是这时低频段的噪声会很强，将拖缆放置在较深的地方进行采集可以减少数据中的噪声，但是由于鬼波的陷波作用，使得地震数据频带变窄。过去，由于油气勘探的目标层埋藏较浅，大多数厚度较薄，需要高分辨率的地震资料才能识别出较薄的含油气储层，因此采集中通常将拖缆尽量放置在较浅的深度上，使鬼波的第一个陷波点位于地震成像所需要的频段之外。随着油气勘探的发展，在中浅层找到油气藏的几率越来越小，勘探的目标层逐渐转向中深层。但是由于上部地层的反射，地震信号到达中深层时能量大大减弱，地震资料的信噪比降低。为了提高中深层的信噪比，野外地震资料采集时，将电缆放置在较深的地方进行采集。这种方式在提高地震资料信噪比的同时，也使得鬼波的陷波影响大大增加。

为了压制鬼波的陷波影响，人们在地震数据采集和室内处理中采用各种办法，如上下缆采集方法就是其中一种。在20世纪50年代，Haggerty最早提出利用上下缆采集方法提高海洋拖缆资料品质的概念，上下缆采集方法就是采集船同时拖拽两条放置在不同深度上的电缆进行采集，较深的拖缆主要测量低频信息，较浅的拖缆填补深部拖缆的陷波点，然而这会增加采集的复杂程度和野外操作难度。双检采集方法是另一种鬼波压制的野外采集方法，双检采集方法在水检(压力检波器)之外增加了陆检(速度检波器)，两套检波器集成于电缆的同一位置，分别记录地震波在海水传播过程中引起的压力波场和速度波场变化，对两种检波器的输出信号进行最佳比例合并，即可压制鬼波。鬼波压制的第三个方法是变深度缆采集，在此种采集方法中，电缆放置在非零倾角的斜线上或者弧形线上，以使鬼波的陷波点分散化。变深度拖缆采集作为一种新的采集方法，它的目的是通过使电缆中不同检波器的陷波点分散，从而有利于压制鬼波，同时将拖缆上远离船的部分放置尽量深，使得低频段得到高信噪比的数据。和传统采集的固定深度不同，拖缆的深度随着偏移距而变化，拖缆中不同位置的检波器的放置深度不同，每个检波器的陷波频率不同，从而使最终综合了不同检波器的结果没有明显的陷波点。处理变深度采集数据关键步骤是检波点鬼波的压制。变深度缆采集的地震数据为我们进一步压制鬼波提供一个非常好的数据基础，但是这种数据要通过特定的处理方法才能达到压制鬼波和拓宽地震数据的频带的目的。2010年，Soubaras等人对变深度缆采集的地震数据的处理方法重新进行研究，提出通过联合反褶积方法来压制虚反射，该方法的基本思想是对变深度缆数据进行偏移和镜像偏移，在这两个偏移数据的基础上进行联合反褶积来压制鬼波，取得了较好的效果。2012年，Soubaras又将该方法推广到叠前数据，这个方法的优点是它把接收点的鬼波看作是一种特殊的波，但是该方法需要完成两次偏移，增加计算量，因此研究人员利用优化滤波算子来解决双偏移计算的问题，针对叠前数据实现鬼波压制。研究人员针对变深度缆地震数据提出一种频率域Radon变换方法，基于接收点鬼波的射线路径，给出一种Radon逆变换算子，该方法能够以较高的精度预测鬼波，然后从实际接收的地震数据中减去预测的鬼波，进而达到鬼波压制的目的。这些方法的应用大大提高了鬼波的压制效果，但是这些方法都是对波动方程进行简化后提出来的，适应于地下构造较简单的地下介质。随着油气勘探的发展，主要勘探的含油气地层越来越复杂，对于这种复杂的地下介质，上述方法压制鬼波的精度降低，很难完全鬼波压制。

技术实现要素：

针对上述问题，本发明的目的是提供一种适用于复杂地下介质且能够提高鬼波压制精度的基于波动方程边值反演的变深度缆鬼波压制方法。

为实现上述目的，本发明采取以下技术方案：一种基于波动方程边值反演的变深度缆鬼波压制方法，其特征在于，包括以下步骤：1)对海上野外变深度缆采集的地震数据抽取共炮点道集；2)对共炮点道集中某炮的地震数据进行一维傅里叶变换，得到含有鬼波的地震数据的频谱；3)计算频谱中某频率ω对应的含有鬼波的频率域地震数据U_AB；4)根据输入的检波点深度计算得到海水面下方基准面的GREEN函数G₁和基准面相对于海水面的对称线的GREEN函数G₂；5)根据得到的含有鬼波的频率域地震数据U_AB以及计算得到的G₁和G₂，计算得到频率ω对应的鬼波压制后的频率域地震数据F_CD；6)对共炮点道集中频谱上的所有频率ω均进行步骤3)～5)得到鬼波压制后的频率域地震数据；7)对共炮点道集中其他炮均进行步骤2)～6)得到鬼波压制后的频率域地震数据；8)对鬼波压制后的频率域地震数据F_CD进行傅里叶逆变换，得到鬼波压制后的共炮点道集。

进一步地，计算频谱中某频率ω对应的含有鬼波的频率域地震数据U_AB的具体过程为：定义基准面为CD，海水面为EF，变深度缆为AB，在CD、EC和FD段均进行观测，得到的地震记录分别为u_CD、u_EC和u_FD，在海水中变深度缆AB的地震记录为u_AB，根据波动方程理论，由u_CD到u_AB的正演问题写为散射形式的正问题：

u(x,z,t)|_(x，z)∈EC＝u_EC (8)

u(x,z,t)|_(x,z)∈FD＝u_FD (9)

取海水面EF的边界条件为自由边界条件，即：

在时间t＝0的时刻，波场尚未开始传播，初始条件取为：

其中，u(x,z,t)为位移波场，x,z分别为水平和垂直方向的坐标，v(x,z)为介质中(x,z)点的速度，f_CD为由u_CD积分得到的散射源函数，对于实际问题，u_EC和u_FD是无法求出的，因此在EC和FD上取吸收边界条件，由于正演问题是在海水中考虑的波的传播问题，因而速度可以取常值，即：

由于水表是负反射，反射系数为-1，因此地震波分成两部分，一部分是直接向上传播的波表示为u₁，另一部分是经过海水面EF反射后向下传播的波表示为u₂，总波场表示为u＝u₁+u₂，且u₁满足：

u₂满足：

其中，C’D’为基准面CD相对于海水面EF的对称线，f_C’D’为f_CD的镜像函数,且f_C’D’与f_CD大小相等，此时，公式(13)和(14)均为均匀无限介质中的方程，对公式(13)和(14)两边均关于时间取傅里叶变换得到：

将公式(15)和(16)通过GREEN函数表示并求解得到：

其中，G₁、G₂为GREEN函数，公式(17)和(18)在变深度缆AB上的接收点也成立，即：

将上述公式(19)和(20)写成矩阵形式得到：

U_AB1＝G₁F_CD (21)

U_AB2＝-G₂F_CD (22)

将公式(21)和(22)合并得到含有鬼波的频率域地震数据U_AB：

U_AB＝G₁F_CD-G₂F_CD＝(G₁-G₂)F_CD (23)

进一步地，根据输入的检波点深度计算得到海水面下方基准面的GREEN函数G₁和基准面相对于海水面的对称线的GREEN函数G₂具体过程为：在均介质中，GREEN函数的计算公式为：

其中，参数k＝ω/v；参数

根据输入的检波点深度及公式(24)和(25)分别计算得到G₁和G₂：

进一步地，根据含有鬼波的频率域地震数据U_AB以及计算得到的G₁和G₂，计算得到频率ω鬼波压制后的频率域地震数据F_CD具体过程为：公式(23)中，G₁-G₂是一个矩阵，根据公式(25)～(27)计算得到G₁-G₂，根据海上野外地震数据采集时变深度缆AB的沉放深度情况获得检波点深度F_CD是不含鬼波的频率域地震数据，为待求量；U_AB是野外地震采集时得到的含有鬼波的频率域地震数据，为已知量，通过公式(23)的求解得到不含有鬼波的频率域地震数据F_CD。

本发明由于采取以上技术方案，其具有以下优点：1、本发明利用波动方程来描述鬼波的产生过程，因而具有较高的精度，能够很好地描述地震波在地下介质中的传播及鬼波的产生过程。2、本发明基于波动方程边值反演，将鬼波的压制问题归结为波动方程的边值反演问题，利用波动方程反演方法解决鬼波压制的问题，鬼波压制效果好。3、本发明针对共炮点道集进行鬼波压制，可以较好地保持原始数据的叠前特征，且本发明在频率域进行，对每个频率可以独立进行处理，因而具有较好的并行性且并行效率较高，压制鬼波后可以拓宽地震信号频谱，提高地震资料的分辨率，本发明可以广泛用于海上变深度缆采集地震资料的处理过程中。

附图说明

图1是本发明的原理示意图，其中，是鬼波的传播路线，是有效波(不包含鬼波)的传播路线；

图2是本发明的实施流程图；

图3是本发明鬼波压制后的数值模拟实验结果与合成的不含鬼波的数值模拟实验结果对比图，其中，图3(a)是合成的无鬼波的单炮道集模拟地震记录，图3(b)是合成的含鬼波的单炮道集模拟地震记录，图3(c)是本发明鬼波压制后的单炮道集模拟地震记录；

图4是本发明鬼波压制后的数值模拟实验频谱与合成的不含鬼波的数值模拟实验频谱对比图，其中，图4(a)是合成的无鬼波的单炮道集振幅谱，图4(b)是合成的含鬼波的单炮道集振幅谱，图4(c)是本发明鬼波压制后的单炮道集振幅谱；

图5是本发明的实际变深度缆地震数据压制鬼波后的处理结果与鬼波压制前的结果对比图，其中，图5(a)是变深度缆采集的一个地震单炮道集数据，图5(b)是本发明压制鬼波后的地震单炮道集数据。

具体实施方式

以下结合附图来对本发明进行详细的描绘。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。

如图1所示，地震波从O点激发后向地下传播，遇到每一个地下反射界面后被反射向上传播，首先经过基准面CD，此时的地震波是不含有鬼波的，地震波继续向上传播，一部分直接传播到变深度缆AB，另一部分传播到海水面EF后被海水面全反射向下传播到变深度缆AB，这部分地震波就是鬼波。因此，变深度缆AB接收到的地震波中既包含了经过基准面CD后向上直接传播到变深度缆AB的地震波场，也包含了被海水面全反射向下传播的地震波场，即鬼波。所以从波场传播的角度看，鬼波压制就是已知变深度缆AB接收的地震记录，如何通过反演求取不含有鬼波的基准面CD的地震记录。

上述通过物理角度对鬼波压制进行说明，下面通过数学角度对鬼波压制进行说明，地震波在二维介质中的传播可用如下的波动方程描述：

其中u(x,z,t)为位移波场，x,z分别为水平和垂直方向的坐标，v(x,z)为介质中(x,z)点的速度，δ()为Diracδ-函数，g(t)为震源函数。假设在O点激发地震波，在CD、EC和FD段均进行观测，得到的地震记录分别为u_CD、u_EC和u_FD。通过地震记录在区域CDFE内考虑地震波的传播问题，主控方程为公式(1),将u_CD、u_EC和u_FD分别作为CD、EC、FD的边界条件，即：

u(x,z,t)|_(x,z)∈CD＝u_CD (2)

u(x,z,t)|_(x,z)∈EC＝u_EC (3)

u(x,z,t)|_(x,z)∈FD＝u_FD (4)

考虑到地震勘探的实际情况，取海水面EF的边界条件为自由边界条件，即：

初始条件取为：

通过波动方程公式(1)及其边界和初始条件公式(2)～(6)的数值求解可以得到海水中变深度缆AB上的地震记录u_AB，同时定义了一个由u_CD到u_AB的波场正向传播过程即正演问题。而鬼波的压制则是由u_AB求解u_CD，这是波场的反向传播过程即波动方程公式(1)的反演问题，在求解反演问题时，未知的是波动方程的边值。在数学领域，波动方程公式(1)是偏微分形式，可以称为偏微分方程或微分方程，求解该方程的反演问题，可称为反问题。上述鬼波压制是由u_AB求解u_CD的过程，数学上为求解微分方程的反问题，也称为边值反问题，地震勘探中称为波动方程边值反演。

如图2所示，基于上述物理角度和数学角度的分析，本发明的基于波动方程边值反演的变深度缆鬼波压制方法，具体包括以下步骤：

1)对海上野外变深度缆采集的地震数据抽取共炮点道集，该道集中包含鬼波数据。

2)对共炮点道集中某炮的地震数据进行一维傅里叶变换，得到含有鬼波的地震数据的频谱。

3)计算频谱中某频率ω对应的含有鬼波的频率域地震数据U_AB。

根据波动方程理论，上述由u_CD到u_AB的正演问题可以写为散射形式的正问题：

u(x,z,t)|_(x,z)∈EC＝u_EC (8)

u(x,z,t)|_(x,z)∈FD＝u_FD (9)

取海水面EF的边界条件为自由边界条件，即：

在时间t＝0的时刻，波场尚未开始传播，初始条件取为：

其中，f_CD为由u_CD积分得到的散射源函数。对于实际问题，u_EC和u_FD是无法求出的，因此在EC和FD上取吸收边界条件。由于上述正演问题是在水中考虑的波的传播问题，因而速度可以取常值，公式(7)中v(x,z)＝v，即：

进一步，由于水表是负反射，通常反射系数为-1，因此地震波可以分成两部分，一部分是直接向上传播的波表示为u₁，另一部分是经过海水面EF反射后向下传播的波表示为u₂，总波场表示为u＝u₁+u₂，且u₁满足：

u₂满足：

其中，C’D’为基准面CD相对于海水面EF的对称线，f_C’D’为f_CD的镜像函数,且f_C’D’与f_CD大小相等。此时，公式(13)和(14)均为均匀无限介质中的方程，对公式(13)和(14)两边均关于时间取傅里叶变换得到：

将公式(15)和(16)通过GREEN函数表示并求解得到：

其中，G₁、G₂为GREEN函数，具体形式后面给出。公式(17)和(18)在变深度缆AB上的接收点也成立，即：

将上述公式(19)和(20)写成矩阵形式得到：

U_AB1＝G₁F_CD (21)

U_AB2＝-G₂F_CD (22)

将公式(21)和(22)合并得到含有鬼波的频率域地震数据U_AB：

U_AB＝G₁F_CD-G₂F_CD＝(G₁-G₂)F_CD (23)

4)根据输入的检波点深度计算得到海水面EF下方基准面CD的GREEN函数G₁和基准面CD相对于海水面EF的对称线C’D’的GREEN函数G₂。

在均介质中，GREEN函数的计算公式为：

其中，参数k＝ω/v；参数

根据输入的检波点深度及公式(24)和(25)分别计算得到G₁和G₂：

5)根据得到的含有鬼波的频率域地震数据U_AB以及计算得到的G₁和G₂，计算得到频率ω对应的鬼波压制后的频率域地震数据F_CD。

公式(23)中，G₁-G₂是一个矩阵，根据公式(25)～(27)可以计算得到G₁-G₂，根据海上野外地震数据采集时变深度缆AB的沉放深度情况，可以获得检波点深度F_CD是不含鬼波的频率域地震数据，为待求量；U_AB是海上野外地震采集时得到的含有鬼波的频率域地震数据，为已知量；通过公式(23)的求解可以得到不含有鬼波的频率域地震数据F_CD，从而达到压制鬼波的目的。

6)对共炮点道集中频谱上的所有频率ω均进行步骤3)～5)得到鬼波压制后的频率域地震数据。

7)对共炮点道集中其他炮均进行步骤2)～6)得到鬼波压制后的频率域地震数据。

8)对鬼波压制后的频率域地震数据F_CD进行傅里叶逆变换，得到鬼波压制后的共炮点道集。

如图3所示，本发明分析了鬼波压制前、后的数值模拟实验结果，图3(a)为合成的无鬼波的单炮道集模拟地震记录，图3(b)为合成的含鬼波的单炮道集模拟地震记录，由图3(b)的模拟地震记录利用本发明压制鬼波后得到图3(c),对比图3(c)和图3(a)可以看出，两幅图完全一致，说明图3(b)中的鬼波已经被压制，且有效信号得到很好的保留。

如图4所示，图4是图3中各个图所对应的频谱图，图4(a)为合成的无鬼波的单炮道集振幅谱，图4(b)为合成的含鬼波的单炮道集振幅谱，可以清楚的看到在图4(b)中鬼波的陷波作用，图4(c)是本发明鬼波压制后的单炮道集振幅谱，从图4(c)中可以看到陷波带的频谱已经恢复。

如图5所示，图5是应用本发明对实际变深度缆地震数据进行鬼波压制前、后的处理结果对比，图5(a)为变深度缆采集的一个地震单炮道集数据，在这幅图中可以清楚的看到鬼波，图5(b)为用本发明压制鬼波后的地震单炮道集数据，对比这两幅图可以清楚的看到，通过本发明处理后的原始地震单炮数据中的鬼波得到很好的压制，同时原始地震单炮数据中的有效波得到很好的保持。

上述各实施例仅用于说明本发明，其中方法的各步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。