CN108181652A - 一种海底节点地震资料上下行波场数值模拟方法 - Google Patents

Abstract

本发明公开了一种海底节点地震资料上下行波场数值模拟方法，包括以下步骤：设计海洋勘探的速度模型和海底节点的观测系统以及相关的地球物理参数；通过声波波动方程交错网格有限差分正演模拟；得到压力分量波场和速度分量波场v；对压力分量与速度分量做振幅、相位的匹配，对于正演数据，求取每一道最大振幅比a，将系数乘以速度分量，上行波记录和下行波记录。该方法适应了当今海底节点地震资料研究的大趋势，对海底节点的波场研究具有现实意义和实用价值。

Description

一种海底节点地震资料上下行波场数值模拟方法

技术领域

本发明属于海洋地球物理地震探勘技术领域，具体地说，涉及一种海底节点地震资料上下行波场数值模拟方法。

背景技术

随着海底节点技术的兴起，多波多分量的波场特征研究越来越受到地球物理学家的关注。现阶段的海底节点波场正演模拟大多关注于三分量波场的研究，包括P分量、Z分量和水平分量。然而现阶段上下行波场的研究已经成为海底节点地震资料研究的重点，现阶段已有的正演模拟软件中，大多数服务于单炮记录的整个波场模拟，在整个波场记录的基础上通过特定的数学手段实现上下行波波场分离，分离的效果与波场记录的复杂程度、应用的数学手段有关，但无法实现简单高效。

发明内容

有鉴于此，本发明针对已有的正演模拟软件无法简单高效地满足上下行波波场分离，为上下行波研究提供技术前提的问题，提供了一种海底节点地震资料上下行波场数值模拟方法，该方法适应了当今海底节点地震资料研究的大趋势，对海底节点的波场研究具有现实意义和实用价值。

为了解决上述技术问题，本发明公开了一种海底节点地震资料上下行波场数值模拟方法，包括以下步骤：

步骤1、设计海洋勘探的速度模型；

步骤2、通过声波波动方程交错网格有限差分正演模拟；

步骤3、得到压力分量波场和速度分量波场，压力分量波场P＝P_x+P_z；速度分量波场为v＝vx+vz；

步骤4、对压力分量与速度分量做振幅、相位的匹配，得到上行波记录和下行波记录。

进一步地，所述海洋勘探的速度模型如下：海底深度450m，自上而下依次设置有第一层速度、第二层速度和第三层速度；所述第一层海水速度1500m/s，第二层海水速度2000m/s，第三层海水速度2500m/s，检波器在海底，炮点在海面，网格间距3m，采样率0.5ms，记录长度2s。

进一步地，所述步骤2中的通过声波波动方程交错网格有限差分正演模拟具体为：

可压缩流体中的标量声波方程为：

其中，P代表压力，C代表声波的速度，t是旅行时，而对于二维速度-深度模型有：

其中，x表示横向坐标，z则表示纵向坐标，由牛顿第二定理得：

其中，ρ表示密度，v_x、v_z分别为x、z方向的速度分量；

将式(2)、(3)代入(1)中，因为求导结果与顺序无关，则可得到：

其中，v表示速度；

等式(4)两边同时对时间积分，得：故地震波传播规律能近似地用一阶速度-应力声波方程组表达：

对一阶速度-应力声波方程组应用交错网格技术，则在空间网格中，应力P位于网格的各个网格点上，如(i，j)；v_x位于i方向上为整数、j方向为半网格数的点上，如v_z位于i方向上为半网格数、j方向上为整数的点上，如

与空间网格类似，在时间网格中应力P位于整数网格点上，而v_x、v_z位于时间半网格点上。

进一步地，时间二阶精度交错网格有限差分格式具体如下：

设v_x在m阶时间可导，在t处利用泰勒展式对和展开：

其中，Δt表示时间网格步距，O(Δt^m+1)表示Δt的m+1阶无穷小量。

两式相减得：

同理得：

将P对t+Δt和t在处展开，得：

即得到了时间二阶差分格式，若用网格位置则可以表示为：

进一步地，对于空间任意偶数阶精度交错网格有限差分格式具体如下：

对于存在2L+1阶导数的任意函数U(x)，则可对i+1和i在点处展开：

其中，Δx表示空间横向网格步距，O(Δx^2L+2)表示Δx的2L+1阶无穷小量；

两式相减得：

以此类推，得到对任意的m＝0,1,2,3,…,L有：

故表示为：

式中，a_m表示差分系数；

根据线性方程组(17)求得a_m，其中m＝1,2,3,····,L

解得：

则对于P表示为：

进一步地，，时间二阶、空间任意偶数阶精度交错网格有限差分算法实现如下：

将式(6)、式(18)代入式(9)，得：

将式(7)、式(23)代入式(13)，得：

将式(5)、式(25)、式(26)代入式(15)得：

式中，K＝ρv²，是介质的物性参数；

式(27)、(28)、(29)也就是时间精度为二阶、空间为任意偶数阶的一阶速度-应力方程组交错网格有限差分格式；

对于一阶速度-应力方程组，PML处理是将速度场与应力场分别分解成平行及垂直于透射界面的两个分量，然后再分别对其做吸收衰减，其衰减函数为：

其中，R是理论反射系数，δ是PML的厚度，V_p是纵波速度；

由于PML把应力场分解成平行及垂直于透射界面的两个分量，故应力P写成P＝P_x+P_z，对于一阶速度-应力方程组，加PML之后的差分形式为：

进一步地，所述步骤4中的对压力分量与速度分量做振幅、相位的匹配，得到上行波记录和下行波记录具体为：

步骤4.1、切除直达波，得到切除直达波之后的模拟的海底节点共检波点道集数据；

步骤4.2、计算压力分量P、速度分量Z对应的最大振幅比a；

步骤4.3、将系数分别乘以速度分量的对应道，得到匹配之后的速度分量z1，上行波记录U＝(P+z1)/2；下行波记录D＝(P-z1)/2；

其中，P表示模拟得到压力分量，z1表示模拟得到的速度分量z经过上述方式匹配得到的速度分量，U表示最后模拟得到的上行波记录，D表示模拟得到的下行波记录。

与现有技术相比，本发明可以获得包括以下技术效果：

交错网格有限差分模拟海底节点记录的速度分量和压力分量；基于速度分量与压力分量响应特征，直接模拟出上行波和下行波记录，直接模拟的上下行波能够服务于上下行波场的研究，如镜像偏移、上下行波反褶积等。

当然，实施本发明的任一产品并不一定需要同时达到以上所述的所有技术效果。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明海底节点波场P分量正演结果图；

图2是本发明模拟所用的雷克子波图；

图3是本发明空间网格分布；

图4是本发明时间网格分布；

图5是本发明完全匹配层的添加方法；

图6是本发明P、Z分量极性响应示意图；

图7是本发明水、陆检在海底接收到的上下行波场；

图8是本发明海底节点波场P分量正演结果图；

图9是本发明海底节点波场Z分量正演结果图；

图10是本发明海底节点波场上行波正演结果图；

图11是本发明海底节点波场下行波正演结果图；

图12是本发明模拟P分量及Z分量；其中，a代表模拟P分量，b代表模拟Z分量；

图13是本发明某软件合成的上行波和下行波；其中，a代表某软件合成的上行波，b代表某软件合成的下行波；

图14是本发明本发明合成的上行波和下行波；其中，a代表本发明合成的上行波，b代表本发明合成的下行波。

具体实施方式

以下将配合实施例来详细说明本发明的实施方式，藉此对本发明如何应用技术手段来解决技术问题并达成技术功效的实现过程能充分理解并据以实施。

本发明公开了一种海底节点地震资料上下行波场数值模拟方法，包括以下步骤：

步骤1、设计海洋勘探的速度模型；

下面以设计的水平层状海底模型为例，本次模拟的海洋勘探的速度模型如附图1中所示，海底深度450m，海水速度1500m/s，第二层速度2000m/s，第三层速度2500m/s，边界用PML边界吸收，检波器在海底，炮点在海面，网格间距3m，采样率0.5ms，记录长度2s。本次模拟用的子波是地震模拟常用的雷克子波，具体形态如附图2所示。

步骤2、通过声波波动方程交错网格有限差分正演模拟；具体过程如下：

可压缩流体中的标量声波方程为：

上式，P代表压力，C代表声波的速度，t是旅行时。而对于二维速度-深度模型有：

式中，x表示横向坐标，z则表示纵向坐标。由牛顿第二定理可得：

其中，ρ表示密度，v_x、v_z分别为x、z方向的速度分量。

将式(2)、(3)代入(1)中，因为求导结果与顺序无关，则可得到：

其中，v表示速度。

等式(4)两边同时对时间积分，可得：故地震波传播规律能近似地用一阶速度-应力声波方程组表达：

对一阶速度-应力声波方程组应用交错网格技术，则在空间网格中，应力P位于网格的各个网格点上，如(i，j)等；v_x位于i方向上为整数、j方向为半网格数的点上，如等；v_z位于i方向上为半网格数、j方向上为整数的点上，如等，相应的示意图(见图3)。

与空间网格类似，在时间网格中应力P位于整数网格点上，而v_x、v_z位于时间半网格点上(见图4)。

对于时间二阶精度交错网格有限差分格式：

设v_x在m阶时间可导，在t处利用泰勒展式对和展开：

其中，Δt表示时间网格步距，O(Δt^m+1)表示Δt的m+1阶无穷小量。

两式相减得：

同理得：

将P对t+Δt和t在处展开，得：

即得到了时间二阶差分格式，若用网格位置则可以表示为：

对于空间任意偶数阶精度交错网格有限差分格式：

对于存在2L+1阶导数的任意函数U(x)，则可对i+1和i在点处展开：

其中，Δx表示空间横向网格步距，O(Δx^2L+2)表示Δx的2L+1阶无穷小量。

两式相减得：

以此类推，可以得到对任意的m＝0,1,2,3,…,L有：

故可以表示为：

式中，a_m表示差分系数。

根据线性方程组(17)可以求得a_m，其中m＝1,2,3,····,L

解得：

则对于P可表示为：

对于时间二阶、空间任意偶数阶精度交错网格有限差分算法实现：

将式(6)、式(18)代入式(9)，得：

将式(7)、式(23)代入式(13)，得：

将式(5)、式(25)、式(26)代入式(15)得：

式中，K＝ρv²，是介质的物性参数。

式(27)、(28)、(29)也就是时间精度为二阶、空间为任意偶数阶的一阶速度-应力方程组交错网格有限差分格式。

完全匹配层是一种理论上的各向异性介质。由于任何频率的波以任意角度穿过它与自由空间的界面后，其相速度和特征阻抗均不发生变化，所以理论上不会发生波的反射现象。而透射波的能量则在介质中被迅速地吸收衰减，因此用PML处理边界问题效果非常理想。对于一阶速度-应力方程组，PML实际上是将速度场与应力场分别分解成平行及垂直于透射界面的两个分量，然后再分别对其做吸收衰减。其衰减函数为：

其中，R是理论反射系数，δ是PML的厚度，V_p是纵波速度。

完全匹配层的添加方法见图5。图中灰色区域为模拟地震波传播的区域，其周围是完全匹配层。d(x)、d(z)各自代表了x及z方向上的阻尼因子，区域1中d(x)≠0，d(z)≠0，速度V取角点的速度。区域2中d(x)＝0，d(z)≠0，速度V取边界速度。区域3中d(x)≠0，d(z)＝0，速度V取边界速度。这样当地震波传播至完全匹配层与区域的边界时就不发生反射，并且其能量会在完全匹配层中随传播距离呈指数衰减，最终在到达完全匹配层边界时约等于零。

由于PML把应力场分解成平行及垂直于透射界面的两个分量，故应力P可写成P＝P_x+P_z，对于一阶速度-应力方程组，加PML之后的差分形式为：

步骤3、根据以上的交错网格差分格式，得到压力分量P＝P_x+P_z，速度分量v＝vx+vz；

为阐述上下行波的合成原理，首先需要说明的是速度分量与压力分量的响应原理。

针对压缩波场和膨胀波场，P分量记录和Z分量记录具有不同的响应，如图6。P分量响应的水中压力随着压缩和膨胀作用的改变而发生极性变化。当P分量检波器受到压缩时产生负脉冲,受到膨胀时产生正脉冲。Z分量检波器是响应质点运动方向的，它随着质点运动方向的变化而发生极性转换。当质点运动方向向上时，Z分量检波器表现出负极性，当质点运动方向向下时,Z分量检波器表现出正极性，如图7所示。

步骤4、对压力分量与速度分量做振幅、相位的匹配，对于正演数据，首先切除直达波，得到切除直达波之后的模拟的海底节点共检波点道集数据。然后，计算压力分量P、速度分量Z对应的最大振幅比a。然后，将系数分别乘以速度分量的对应道，得到匹配之后的速度分量z1，上行波记录U＝(P+z1)/2；下行波记录D＝(P-z1)/2；其中，P表示模拟得到压力分量，z1表示模拟得到的速度分量z经过上述方式匹配得到的速度分量，U表示最后模拟得到的上行波记录，D表示模拟得到的下行波记录。

本发明结合当今海底节点研究的趋势，鉴于目前缺少上下行波波场正演的软件，本发明基于波动方程正演理论，提出了海底节点上下行波场的正演方法，填补了这方面的空白。

本发明给出了海底节点地震资料上下行波场模拟的方法，P、Z分量的正演(图8和图9)对于研究双检机制、双检合成等有重要的作用；上下行波的正演(图10和图11)对于研究上下行波反褶积、OBS/OBC/OBN数据下行波镜像偏移等提供理论研究基础。

为说明本发明的实用性，以某市场软件的模拟结果与本发明的模拟结果做对比，对比结果如附图中，图12是模拟的原始速度分量和压力分量的原始数据，切除直达波。图13是某市场软件分离的上下行波，图14是本发明分离的上下行波，两者对比，明显看到本发明的分离效果更好，上行波记录中下行波极大被压制，下行波记录中上行波也极大得到压制，避免两者相互干扰，服务于上下行波的进一步深入研究。

上述说明示出并描述了发明的若干优选实施例，但如前所述，应当理解发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离发明的精神和范围，则都应在发明所附权利要求的保护范围内。

Claims (7)

1.一种海底节点地震资料上下行波场数值模拟方法，其特征在于，包括以下步骤：

步骤1、设计海洋勘探的速度模型；

步骤2、通过声波波动方程交错网格有限差分正演模拟；

步骤3、得到压力分量波场和速度分量波场，压力分量波场P＝P_x+P_z；速度分量波场为v＝vx+vz；

步骤4、对压力分量与速度分量做振幅、相位的匹配，得到上行波记录和下行波记录。

2.根据权利要求1所述的模拟方法，其特征在于，所述海洋勘探的速度模型如下：海底深度450m，自上而下依次设置有第一层速度、第二层速度和第三层速度；所述第一层海水速度1500m/s，第二层海水速度2000m/s，第三层海水速度2500m/s，检波器在海底，炮点在海面，网格间距3m，采样率0.5ms，记录长度2s。

3.根据权利要求1所述的模拟方法，其特征在于，所述步骤2中的通过声波波动方程交错网格有限差分正演模拟具体为：

可压缩流体中的标量声波方程为：

其中，P代表压力，C代表声波的速度，t是旅行时，而对于二维速度-深度模型有：

其中，x表示横向坐标，z则表示纵向坐标，由牛顿第二定理得：

其中，ρ表示密度，v_x、v_z分别为x、z方向的速度分量；

将式(2)、(3)代入(1)中，因为求导结果与顺序无关，则可得到：

其中，v表示速度；

等式(4)两边同时对时间积分，得：故地震波传播规律能近似地用一阶速度-应力声波方程组表达：

对一阶速度-应力声波方程组应用交错网格技术，则在空间网格中，应力P位于网格的各个网格点上，如(i，j)；v_x位于i方向上为整数、j方向为半网格数的点上，如v_z位于i方向上为半网格数、j方向上为整数的点上，如

与空间网格类似，在时间网格中应力P位于整数网格点上，而v_x、v_z位于时间半网格点上。

4.根据权利要求3所述的模拟方法，其特征在于，时间二阶精度交错网格有限差分格式具体如下：

设v_x在m阶时间可导，在t处利用泰勒展式对和展开：

其中，Δt表示时间网格步距，O(Δt^m+1)表示Δt的m+1阶无穷小量。

两式相减得：

同理得：

将P对t+Δt和t在处展开，得：

即得到了时间二阶差分格式，若用网格位置则可以表示为：

5.根据权利要求4所述的模拟方法，其特征在于，对于空间任意偶数阶精度交错网格有限差分格式具体如下：

对于存在2L+1阶导数的任意函数U(x)，则可对i+1和i在点处展开：

其中，Δx表示空间横向网格步距，O(Δx^2L+2)表示Δx的2L+1阶无穷小量；

两式相减得：

以此类推，得到对任意的m＝0,1,2,3,…,L有：

故表示为：

式中，a_m表示差分系数；

根据线性方程组(17)求得a_m，其中m＝1,2,3,....,L

解得：

则对于P表示为：

6.根据权利要求5所述的模拟方法，其特征在于，时间二阶、空间任意偶数阶精度交错网格有限差分算法实现如下：

将式(6)、式(18)代入式(9)，得：

将式(7)、式(23)代入式(13)，得：

将式(5)、式(25)、式(26)代入式(15)得：

式中，K＝ρv²，是介质的物性参数；

式(27)、(28)、(29)也就是时间精度为二阶、空间为任意偶数阶的一阶速度-应力方程组交错网格有限差分格式；

对于一阶速度-应力方程组，PML处理是将速度场与应力场分别分解成平行及垂直于透射界面的两个分量，然后再分别对其做吸收衰减，其衰减函数为：

其中，R是理论反射系数，δ是PML的厚度，V_p是纵波速度；

由于PML把应力场分解成平行及垂直于透射界面的两个分量，故应力P写成P＝P_x+P_z，对于一阶速度-应力方程组，加PML之后的差分形式为：

7.根据权利要求1所述的模拟方法，其特征在于，所述步骤4中的对压力分量与速度分量做振幅、相位的匹配，得到上行波记录和下行波记录具体为：

步骤4.1、切除直达波，得到切除直达波之后的模拟的海底节点共检波点道集数据；

步骤4.2、计算压力分量P、速度分量Z对应的最大振幅比a；

步骤4.3、将系数分别乘以速度分量的对应道，得到匹配之后的速度分量z1，上行波记录U＝(P+z1)/2；下行波记录D＝(P-z1)/2；

其中，P表示模拟得到压力分量，z1表示模拟得到的速度分量z经过上述方式匹配得到的速度分量，U表示最后模拟得到的上行波记录，D表示模拟得到的下行波记录。