CN104483702B - 一种适用于非均匀运动水体的地震正演模拟方法 - Google Patents

Abstract

本发明提供了一种适用于非均匀运动水体的地震正演模拟方法，该方法首先获取海水和地层中的相关数据、观测系统和震源参数；接着从质量守恒方程、欧拉方程和绝热状态方程出发，利用渐进方法和基于海水实际情况的适当假设，推导适用于非均匀运动水体的一阶应力‑速度方程；然后，给出合理的边界条件；最后采用高精度有限差分将控制方程和边界条件离散化，模拟地震波在海水和海底弹性介质中的传播过程。本方法综合考虑了海水中密度、声速和流速对地震波传播的影响，可以真实地反映地震波在复杂海洋环境中的传播情况，在海洋油气勘探中，尤其是海洋环境比较复杂的深水油气勘探中具有广泛的应用价值。

Description

一种适用于非均匀运动水体的地震正演模拟方法

技术领域

本发明属于海洋地震勘探领域，具体地为深水油气勘探领域，涉及一种针对复杂海洋环境的地震波模拟方法。

背景技术

提高深水区地震成像精度以得到勘探目标真实的深度构造图对深海油气田开发具有重要意义。但是，深水区地震成像目前仍存在很多问题，其中深海海洋环境条件极其复杂就是主要因素之一。海水中存在着中尺度涡、内波和海洋锋等中尺度现象。观测资料发现，当地震波穿过深水海域时中尺度现象会时时地反映在地震数据中(Biescas等，2008；Song等，2008；Ruddick，2003；Pinheiro等，2010)，同时这些非均匀水体所引起的多次超强反射还会覆盖中深层数据。因此我们有必要研究复杂海洋环境(即非均匀运动水体)下地震波的传播特征。

目前的地震勘探中，人们仍主要采用常规的声波方程，将非均匀运动的海水视为均匀介质。Kormann等(Kormann等，2008；2009)采用简化的波动方程和二阶完全匹配层吸收边界条件，模拟了海水层产生的地震反射波。但是在该研究中假设海水的流速为0并且海水中的密度和声速满足(其中K为流体压缩模量)，以上假设有一定的局限性，因为海水是在不停运动的，并且海水中声速和密度的关系比较复杂。姬莉莉等(姬莉莉等，2013)推导了含流速的声波方程，研究了中尺度涡对声波传播的影响，但是其假设海底地层只存在纵波。在海洋油气勘探中，气枪产生的信号在海底会产生纵波和横波，同时由油气层反射的横波信息在海底会重新转换为一些纵波传入水中。因此，为了全面研究复杂海洋环境对地震勘探的影响，海底地层不能简单的处理为声波层。若将海底视为弹性基底，同时考虑纵波和横波的传播，则需要对现有的控制方程和边界条件做出改进，使其在适用于非均匀运动水体的同时，又与海底弹性基底相匹配。

发明内容

针对现有技术的不足及问题，本发明提供了一种适用于非均匀运动水体的地震波正演模拟方法，可以综合全面地考虑海水中密度，声速、流速等因素对地震波传播的影响，又可以模拟海底弹性波的矢量波场。

依据本发明的技术方案，提供一种适用于非均匀运动水体的地震波正演模拟方法，具体包括以下步骤：

1)获取海水和地层中的相关数据、观测系统和震源参数；

2)利用渐进方法和基于海水实际情况的适当假设，推导出适用于非均匀运动水体的一阶应力-速度方程；

3)给出合理的边界条件；

4)采用高精度有限差分将控制方程和边界条件离散化；

5)计算非均匀运动水体影响下接收的地震波场和地震记录。

其中所述步骤1)获取海水和地层中的相关数据、观测系统和震源参数，包括：

获取海洋中的声速、密度和流速数据，海底地层中的密度、纵波速度和横波速度，并把这些数据按网格间距网格化，形成网格化数据文件；

获取观测系统参数包括道间距、最大、最小偏移距和时间采样间隔；

获取震源参数包括震源坐标和震源频率。

此外步骤2)从质量守恒方程、欧拉方程和绝热状态方程出发，利用渐进方法和基于海水实际情况的适当假设，推导出适用于非均匀运动水体的一阶应力-速度方程：

其中，为水中速度分量，是水中密度，是水中的压力，水中的声速，为水中的流速。

进一步地，步骤3)中在模型四周引入PML边界层，同时由于水中的波阻抗较小，为了能够看到水中的反射波，边界层引起的反射必须比水中的反射小一个数量级。基于以上要求，水中PML边界层中控制方程为：

其中，为在x和z方向上的分量， η^T-η₀是PML边界层厚度。

地层中PML边界条件采用与水中类似的方法，对各个变量进行分裂：

得到地层中PML边界层中控制方程为：

其中，为弹性介质中速度分量，是弹性介质中密度， 和是弹性介质中的应力，λ和μ为介质的弹性常数。

优选地，步骤4)中采用交错网格有限差分方法将控制方程和边界条件离散化。

步骤5)用于计算非均匀运动水体影响下接收的地震波场和地震记录，获得非均匀运动水体对接收到地震记录的影响，以提高地震数据质量。

与现有技术相比，本发明的有益效果是：给出了适用于非均匀运动水体的一阶应力-速度方程，从而可以综合考虑水体中密度、流速和声速等因素对地震波传播的影响；采用交错网格有限差分格式，以保证计算的精度；给出的PML边界条件，保证了边界反射比水体中的反射波小一个数量级，这样可以使我们更清楚地观察海水对地震波的反射作用，并精确计算非均匀运动水体对接收的地震数据的影响，以提高地震数据质量。

附图说明

图1是依据本发明的地球物理模型示意图；

图2(a)-图2(b)分别是含有中尺度涡的水中声速(a)和流速(b)分布示意图；

图3是计算模型示意图；

图4(a)-图4(b)分别是水中含中尺度涡时1.67s时刻计算的水平分量(a)和垂直分量(b)的波场快照图；

图5(a)-图5(b)分别是水中不含中尺度涡时1.67s时刻计算的水平分量(a)和垂直分量(b)的波场快照图；

图6(a)-图6(b)分别是有涡和没有涡时在海底接收到的地震波水平分量(a)和垂直分量(b)的对比图；

图7是中尺度涡影响范围内的最小偏移距剖面。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。另外地，不应当将本发明的保护范围仅仅限制至下述具体结构或部件或具体参数。

本发明首先获取海水和地层中的相关数据、观测系统和震源参数；接着从质量守恒方程、欧拉方程和绝热状态方程出发，利用渐进方法和基于海水实际情况的适当假设，推导适用于非均匀运动水体的一阶应力-速度方程；然后给出合理的边界条件；最后用交错网格有限差分将控制方程和边界条件离散化，计算地震波场和地震记录。

本发明提供的适用于非均匀运动水体的地震波正演模拟方法，具体内容包括：

步骤1：获取海水和地层中的相关数据、观测系统和震源参数。

获取海洋中的声速、密度和流速数据，海底地层中的密度、纵波速度和横波速度，并把这些数据按网格间距网格化，形成网格化数据文件；获取观测系统参数包括道间距、最大、最小偏移距和时间采样间隔；获取震源参数包括震源坐标和震源频率。

步骤2：利用渐进方法和基于海水实际情况的适当假设，推导适用于非均匀运动水体的一阶应力-速度方程：

其中，为水中速度分量，是水中密度，是水中的压力，水中的声速，为水中的流速。与目前正演模拟中常用的一阶应力-速度方程相比，上式多了含流速项。之所以写成一阶应力-速度方程的形式，是为了与地层中的弹性波方程相呼应，方便于整体求解。

地层中的控制方程仍采用常规的一阶应力-速度方程，为：

其中，为弹性介质中速度分量，是弹性介质中密度， 和是弹性介质中的应力，λ和μ为介质的弹性常数。

步骤3：给出合理的边界条件。

海底(液固交界面)的边界条件为：法向位移和法向应力连续。

在物理区域四周(上下左右四个面)加入完全匹配边界层(PML)。

步骤4：采用交错网格有限差分法将控制方程和边界条件离散化

PML边界层中方程的离散化如下：

其中 η^T-η₀是PML边界层厚度， 是沿坐标轴进行分裂所得。这样给边界条件是因为水中波阻抗较小，为了观察水中产生的反射波，必须保证边界反射比水中反射小一个数量级。

同理可将地层中的控制方程以及地层周围的PML边界层离散化。

步骤5：计算接收的地震波场和地震记录，获得非均匀运动水体对接收到地震记录的影响。

下面结合附图，对本发明的适用于非均匀运动水体的地震波正演模拟方法进行具体步骤说明：

(1)获取海水和地层中的相关数据、观测系统和震源参数：

为了展示本专利技术的发明效果，给出我国南海深水区的具体实例进行说明。

图1为地球物理模型，上层为海水，水深为2200m。中尺度涡的有效影响半径和有效影响深度分别为150km和1600m。海底有两层弹性介质层。图2给出了含中尺度涡的海水内的密度、声速和流速分布图。表1给出了弹性介质层内的物理参数。

表1

震源的主频为30Hz，坐标深度为50m。

接收器的道间距为10m，最小偏移距为-1500m，最大偏移距为1500m。

(2)推导适用于非均匀运动水体的一阶应力-速度方程。过程具体如下：

理想流体中的质量守恒方程、欧拉方程和绝热状态方程为：

其中是梯度算子， 是流体密度，是流体的速度，是水中的压力，水中的声速。

海水的非均匀性可分为两类：静态非均匀(即把海水看成达到流体静力学平衡的分层流体，以下简称静态海水)和动态非均匀(即海水受到某种扰动，流体静力学平衡被破坏，如旋涡、内波、海洋锋等)，静态海水的非均匀性主要表现为密度和声速的非均匀。海水的运动形式虽然呈现多样化，但是从本质上说都可用密度场、声速场和流速场来描述，动态海水的非均匀性主要体现在这三个物理量是随时空变化的函数。这里假设方程组(5)中的量 由两部分组成(Robertson等，1985)：静态海水时的值；由海水运动(中尺度涡、内波和海洋锋等)所引起的扰动量。前者由下标0表示，后者用下标1表示。如下：

下面将方程(5)和(6)中各无物理量无量钢化，得：

其中和为特征长度和特征时间，在接下来的讨论中会确定出他们的值。接着，假设 和分别为参考密度、参考流速和参考声速，无量纲化下面的各个量

其中u₀，v₁，p₀，p₁，ρ₀和ρ₁为无量纲量，且大小为常数阶，i为沿x轴的单位向量，δ是代表旋涡扰动的无量纲小量。

将方程(6)-(8)代入到方程(5a)可得：

其中M表示马赫数

当没有声速扰动时(δ＝0)，因为只是z的函数，所以方程(9)右边等于0，这样方程左右两边平衡。取O(δ)的项，则得：

$\mathrm{\ζ}\frac{\∂{\mathrm{\ρ}}_1}{\∂t}+\▿\·({\mathrm{\ρ}}_0{v}_1+M{\mathrm{\ρ}}_1{v}_0)=0---\left(10\right)$

其中

类似地，可将方程(6)-(8)代入到方程(5b)，并去掉高阶项，可得：

$\mathrm{\ζ}\frac{{\∂v}_1}{\∂t}+M(v_0\·\▿)v_1+M(v_1\·\▿)v_0=-\frac{\mathrm{\ζ}}{{\mathrm{\ρ}}_0}\▿p_1---\left(12\right)$

最后，来化简方程(5c)。定义折射率并将方程(6)-(8)代入到方程(5c)，可得：

$\begin{array}{c}{n}^{-2}(\mathrm{\ζ\δ}\frac{{\∂\mathrm{\ρ}}_1}{\∂t}+\mathrm{\δ}{v}_1\·\▿{\mathrm{\ρ}}_0+\mathrm{\δ}{\mathrm{Mv}}_0\·{\▿\mathrm{\ρ}}_1+{\mathrm{\δ}}^2{v}_1\·\▿{\mathrm{\ρ}}_1)\\ =M^3{v}_0\·{\▿p}_0+\mathrm{\ζ}{\mathrm{M\δv}}_0\·{\▿p}_1+M^2\mathrm{\δ}{v}_1\·{\▿p}_0\\ +{\mathrm{\ζ\δ}}^2{v}_1\·{\▿p}_1+{\mathrm{\ζ}}^2\mathrm{\δ}\frac{{\∂p}_1}{\∂t}\end{array}---\left(13\right)$

同样地，没有声速扰动时(δ＝0)，因为只是z的函数，所以方程(13)右边等于0，这样方程左右两边平衡。一般情况下，M＜＜1，在δ＝O(M²)的条件下，上面的方程可简化为：

$n^{-2}(\mathrm{\ζ}\frac{\∂{\mathrm{\ρ}}_1}{\∂t}+v_1\·\▿{\mathrm{\ρ}}_0+{\mathrm{Mv}}_0\·\▿{\mathrm{\ρ}}_1)=\mathrm{\ζ}{\mathrm{Mv}}_0\·\▿p_1+{\mathrm{\ζ}}^2\frac{{\∂p}_1}{\∂t}---\left(14\right)$

这里需要指出，δ＝O(M²)的条件是合理的。下面就来分析其合理性，由这个条件和方程(8)可得：

在海洋中，一般情况下M比10^-3小。这样就得到旋涡产生的流速扰动比10^-3ms^-1小。这是对流速的一个合理的估计。

至此，我们将方程(5)的三个方程化简为方程(10)、(12)和(14)。进一步地，可将方程(10)和(14)相加，并化简得：

$-{\mathrm{\ρ}}_0\▿\·v_1=\mathrm{\ζ}{n}^2{\mathrm{Mu}}_0\frac{{\∂p}_1}{\∂x}+{\mathrm{\ζ}}^2{n}^2\frac{{\∂p}_1}{\∂t}---\left(16\right)$

下面重点讨论方程(12)的化简：

${M(v_0\·\▿)v}_1+{M(v_1\·\▿)v}_0{=\mathrm{Mu}}_0\frac{{\∂v}_1}{\∂x}+{\mathrm{Mv}}_{1z}{u}_0^{\′}---\left(17\right)$

其中v_1z为声速v₁在z方向上的分量。

分析u₀′的量纲，以进一步化简方程(12)。其实分析u₀′的量纲是必要的，因为即使很小的流速也可能引起较大的流速梯度。将流速梯度写为如下形式：

u₀′＝ξG(z) (18)

其中G是常数阶量纲。这里定义ξ为|u₀′|在深度上的最大值。

为了得到ξ可能的最大值，我们考虑与最大剪切流速相对应的深度。在此深度附近，G≌1。为了便于分析，这里将方程(18)写成有量纲的形式：

其中为震源的频率，单位为赫兹。

我国南海观测数据表明，即使是在震源频率较低的情况下，也有例如：Wu(吴培木等，2001)观测的一个比较典型的中尺度海洋现象，有而且表面的最大流速为代入到方程(19)得：

当时，可得所以。再由Robertson曾证明当时，含u₀′的项可以忽略(Robertson等，1985)，则方程(12)可化简为如下形式：

$\mathrm{\ζ}\frac{{\∂v}_1}{\∂t}+{\mathrm{Mu}}_0\frac{{\∂v}_1}{\∂x}=-\frac{\mathrm{\ζ}}{{\mathrm{\ρ}}_0}\▿p_1---\left(21\right)$

将方程(16)和方程(21)写成有量纲的形式可得：

与常规的一阶压力-速度方程相比，上式多了含流速项，可以用来描述非均匀运动水体中声速的传播。

在海底弹性介质中，控制方程采用常规的

仍采用常规的一阶应力-速度方程，为：

其中各个变量的定义如前文所述。

(3)给出合理的边界条件。

海底(液固交界面)边界条件，满足法向位移和法向应力连续。

在物理区域四周(上下左右)加入完全匹配(PML)边界层。其中，水中PML边界层中控制方程为：

地层中PML边界条件采用与水中类似的方法，由于弹性介质中变量较多，所以方程比水中要复杂。首先对各个变量进行分裂：

然后，将上述变量代入到方程中，得：

首先选取旋涡中心3000m的范围进行计算，计算模型如图3所示，PML边界层的厚度为100m。

(4)采用交错网格有限差分法将控制方程和边界条件离散化

水中控制方程的离散化为：

水中PML边界层中方程的离散化如下：

地层中控制方程的离散化为：

地层中PML边界层中方程的离散化如下：

(5)计算地震波场和地震记录，获得中尺度涡地震波传播的影响。

计算了地震波场瞬时快照和地震记录。图4(a)和(b)为1.67s时的地震波场快照，可以看出：PML边界有效地吸收的了边界反射；交错网格有限差分算法没有产生数值频散，中尺度涡内产生的反射波清晰可见。图4和图5比较可以发现，有中尺度涡时水中会产生微弱的反射波；地震波穿过水层到达地层之后，有涡和没有涡的波阵面也有不同。图6给出了有涡和没有涡时在(3000m，2200m)接收到的地震波对比图。可以发现，两者的传播时间相差了15ms，有涡时振幅相对无涡时最大变化了10％。

计算整个中尺度涡影响范围的情况，图7给出了最小偏移距剖面。可以看出，水中上层的反射层比较密集。从第一层界面的反射波(2.9s到3.0s的放大图)可以看出，原本应该平直的界面发生弯曲。

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

Claims (3)

1.一种适用于非均匀运动水体的地震正演模拟方法，其特征在于包括以下步骤：

1)获取海水和地层中的相关数据、观测系统和震源参数；

2)利用渐进方法和基于海水实际情况的适当假设，推导出适用于非均匀运动水体的一阶应力-速度方程；

3)给出合理的边界条件；

4)采用高精度有限差分将控制方程和边界条件离散化；

5)计算非均匀运动水体影响下接收的地震波场和地震记录；

所述步骤1)获取海水和地层中的相关数据、观测系统和震源参数，包括：

获取海洋中的声速、密度和流速数据，海底地层中的密度、纵波速度和横波速度，并把这些数据按网格间距网格化，形成网格化数据文件；

获取观测系统参数包括道间距、最大、最小偏移距和时间采样间隔；

获取震源参数包括震源坐标和震源频率；

所述步骤2)从质量守恒方程、欧拉方程和绝热状态方程出发，利用渐进方法和基于海水实际情况的适当假设，推导出适用于非均匀运动水体的一阶应力-速度方程：

其中，为水中声波速度水平分量和垂直分量，是水中密度，是水中的压力，水中的声速，为水中的流速；为弹性波横向传播距离，为弹性波纵向传播深度，为弹性波传播时间；

所述步骤3)中在模型四周引入PML边界层，同时由于水中的波阻抗较小，为了能够看到水中的反射波，边界层引起的反射必须比水中的反射小一个数量级；基于以上要求，水中PML边界层中控制方程为：

其中，分别为在x和z方向上的分量，的定义如下：

上式中η代表和η^T表示完全匹配边界层(PML)外边界的位置，η₀表示完全匹配边界层(PML)内边界的位置(即边界层与模型交界处的位置)，η^T-η₀是完全匹配边界层(PML)厚度；

地层中PML边界条件采用与水中类似的方法，对各个变量进行分裂：

得到地层中PML边界层中控制方程为：

其中，为弹性介质中和方向上的速度分量，是弹性介质中密度，和是分别是弹性介质中和方向上的正应力，是弹性介质中的剪切应力；和分别是在和方向上的分量，和分别是在和方向上的分量，和分别是在和方向上的分量，λ和μ为介质的Lame第一参数和第二参数；和的定义和前面相同。

2.如权利1所述的适用于非均匀运动水体的地震正演模拟方法，其特征在于，所述步骤4)中采用交错网格有限差分方法将控制方程和边界条件离散化。

3.如权利1所述的适用于非均匀运动水体的地震正演模拟方法，其特征在于，所述步骤5)计算非均匀运动水体影响下接收的地震波场和地震记录，获得非均匀运动水体对接收到地震记录的影响，以提高地震数据质量。