ABSTRACT

As a geophysical exploration method, VSP plays an increasingly important role in oil field exploration and development. Through the response of wave field on the vertical arrangement of wells, it can obtain geophysical parameters, velocity field information, horizon identification, well-side structural and lithologic characteristics of strata, etc. The ability to solve these inverse problems has been greatly improved. In order to solve the inverse problem, we must also deal with the positive problem-VSP forward modeling. In this paper, taking well XX as an example, two VSP forward modeling methods based on ray theory and wave equation are implemented by using seismic data and logging data. The differences and advantages and disadvantages of the two methods are compared with the field data.

Keywords:VSP Forward Modeling, Ray Theory, Wave Equation, Ray Tracing, Forward Track Record

两种VSP正演模拟方法应用对比研究

王雪峰¹，范增辉²，刘聪伟¹，郭鸿喜³，李飞¹，李娜¹，黄建华¹，张彦斌¹

¹中国石油东方地球物理公司新兴物探开发处，河北 保定

²中国石油东方地球物理公司西南物探研究院，四川 成都

³中国石油西南油气田分公司勘探事业部，四川 成都

收稿日期：2018年12月19日；录用日期：2019年1月2日；发布日期：2019年1月9日

摘 要

VSP作为一种地球物理勘探方法在油田勘探开发中扮演着越来越重要的角色。它通过井中垂直排列上的波场响应，求地球物理参数、速度场信息、层位识别、井旁构造形态和地层岩性特征等，解决这些反问题的能力得到很大提高。同时为了更好的解决反问题首先也要处理好正问题—VSP正演模拟。下文以XX井为例，综合应用地震资料，测井资料，分别进行了基于射线理论和基于波动方程的两种VSP正演模拟的实现过程，并结合野外实际采集资料对比了二者方法的差异和优劣势。

关键词 :VSP正演模拟，射线理论，波动方程，射线追踪，正演道集记录

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 前言

近年来VSP在“三高”地震勘探中扮演着重要角色，在油田勘探开发中发挥重要作用，解决反问题的能力得到很大提高，而应用这些的前提是解决好正问题。因此，VSP正演模拟受到越来越多的关注。正演模拟主要是基于射线理论和波动方程理论，两种方法都能实现VSP正演模拟，可以从理论上研究VSP的波场特征，指导VSP施工方案设计，提供资料采集质控依据，是VSP资料处理解释中必不可少的环节。迄今为止，国内外很多学者已经提出多种提高数值模拟精度的方法 [1] [2] [3] [4] ，Dablain [5] 应用不同差分近似阶数情况下的声波方程进行波场的对比研究，认为差分近似阶数越高，精度也越高；Alford [6] 指出网格大小和地震波传播方向是影响声波方程有限差分法频散的两个主要因素；裴正林、陈可洋等 [7] [8] 提出了声波、弹性波等的高阶交错网格有限差分方法，实现了地震波场高精度正演数值模拟；Fei T等 [9] 应用通量传输校正(FCT)技术有效压制了数值频散问题；董良国、吴国忱、宁刚、陈可洋等 [10] [11] [12] [13] 研究了正演模拟的频散问题、分析了误差来源及边界吸收问题，并提出了相应的解决办法。在前人研究的基础上，下文重点讨论基于射线理论方程的射线迭代法和基于波动方程的交错网格高阶有限差分法的原理，正演模拟实现过程，射线追踪，正演道集记录，进行了两者正演道集记录对比分析以及与实际野外采集道集记录的比较，同时总结了两者方法的优势和劣势。

2. 射线理论

2.1. 基本方程

各向同性介质波的运动方程如下式 [14] ：

$\rho \frac{\partial \text{2w}}{\partial 2t}=\left(\lambda +\mu \right)\nabla \left(\nabla \cdot \text{w}\right)+\mu \nabla 2\text{w}+\nabla \lambda \left(\nabla \cdot \text{w}\right)+\nabla \mu \times \left(\nabla \times \text{w}\right)+\text{2}\left(\nabla \mu \cdot \nabla \right)\text{w}$ (1)

式中： $\rho$ 是弹性介质密度；W为波函数； $\lambda$ 、 $\mu$ 是拉梅系数； $\nabla$ 是哈密顿算子。

$\nabla =\stackrel{\to }{i}\frac{\partial }{\partial x}+\stackrel{\to }{j}\frac{\partial }{\partial y}+\stackrel{\to }{k}\frac{\partial }{\partial z}$ (2)

$\nabla \mu =\stackrel{\to }{i}\frac{\partial u}{\partial x}+\stackrel{\to }{j}\frac{\partial u}{\partial y}+\stackrel{\to }{k}\frac{\partial u}{\partial z}$ (3)

$\nabla \cdot \text{W}=\stackrel{\to }{i}\frac{\partial \text{Wx}}{\partial x}+\stackrel{\to }{j}\frac{\partial \text{Wy}}{\partial y}+\stackrel{\to }{k}\frac{\partial \text{Wz}}{\partial z}$ (4)

$\nabla \times \text{W}=\left|\begin{array}{ccc}\stackrel{\to }{i}& \stackrel{\to }{j}& \stackrel{\to }{k}\\ \frac{\partial }{x}& \frac{\partial }{y}& \frac{\partial }{z}\\ Wx& Wy& Wz\end{array}\right|$ (5)

假定方程(1)的时间简谐函数的解可用 $\omega$ 的负幂数表示为：

$\text{W}=\exp \left[j\omega \left(t-\tau \right)\right]{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{\left(j\omega \right)}^{-k}{W}_k}$ (6)

式中 $\tau$ 、 $W_k$ 与 $\omega$ 、t没有关系。(6)式为射线级数 [14] ， $W_k\left(k=0,1,2,\cdots \right)$ 为射线级数的振幅系数， $\tau$ 是相位函数，波前是运动着的等相位面 $t=\tau \left(x,y,z\right)$ ，与波前正交的轨迹线称为射线，射线在其与波前交点的方向与该点波前的梯度 $\nabla \tau$ 的方向相同。假定 $\nabla \tau \ne 0$ ，否则(6)式将不表示传播着的波。 $\tau$ 和 $W_k\left(k=0,1,2,\cdots \right)$ 都是坐标的未知函数，如果知道某些初始条件，将(6)式代入(1)式后，可确定 $\tau$ 和 $W_k\left(k=0,1,2,\cdots \right)$ 。

为了求相位函数 $\tau$ 和振幅系数 $W_k$ ，先给出射线理论基本方程组。

$\{\begin{array}{l}N\left(W_0\right)=0\\ N\left(W_1\right)-M\left(W_0\right)=0\\ N\left(W_k\right)-M\left(W_{k-1}\right)+L\left(W_{k-2}\right)=0\left(k\ge 2\right)\end{array}$ (7)

(7)式中

$\{\begin{array}{l}N\left(W_k\right)=-\rho W_k+\left(\lambda +\mu \right)\left(W_k\cdot \nabla \tau \right)\nabla \tau +\mu {\left(\Delta \tau \right)}^2{W}_k\\ M\left(W_k\right)=\left(\lambda +\mu \right)\left[\nabla \left(W_k\cdot \nabla \tau \right)+\Delta \tau \left(\Delta \cdot W_k\right)\right]+\mu \left[2\left(\nabla \tau \cdot \nabla \right)W_k+\left({\nabla }^2\tau \right)W_k\right]\\ +\nabla \lambda \left(W_k\cdot \Delta \tau \right)+\nabla \mu \times \left(\nabla \tau \times W_k\right)+2\left(\Delta \mu \cdot \Delta \tau \right)W_k\\ L\left(W_k\right)=\left(\lambda +\mu \right)\nabla \left(\nabla \cdot W_k\right)+\mu {\nabla }^2{W}_k+\nabla \lambda \left(\nabla \cdot W_k\right)+\nabla \mu \times \left(\nabla \times W_k\right)+z\left(\nabla \mu \cdot \nabla \right)W_k\end{array}$ (8)

且 $\lambda$ 、 $\mu$ 、 $\rho$ 及其相应的导数连续； $\tau$ 是解析函数，且 $\Delta \tau \ne 0,W_0\ne 0$ ； $\lambda$ 、 $\mu$ 、 $\rho$ 及其相应的导数在短距离内的变化不迅速。

方程(7)是射线理论的基本方程，由此方程组可以求 $\tau$ 及所有的 $W_k\left(k=0,1,2,\cdots \right)$ 。

而程函数方程是解决运动学问题的基本方程，能够确定射线、旅行时、波前等波的运动学特性。方程表达式如下：

$V_p={\left[\left(\lambda +2\mu \right)/\rho \right]}^{1/2}$ (9)

$V_s={\left[\mu /\rho \right]}^{1/2}$ (10)

方程(9)是描述P波传播的运动学特征， $V_p$ 是P波速度。

方程(10)是描述S波传播的运动学特征， $V_s$ 是S波速度。

2.2. 试射迭代方法

射线法主要包含两方面内容：1) 在已知射线起始位置和知道射线出射方位的前提下，确定射线路径；2) 给出起始出射方向，逐渐修改出射方向，直到射线足够精确地到达预定的接收点 [15] 。

内容1)是积分问题，已知速度分布，应用程函数方程(9)和(10)，对微分方程作积分求解射线路径问题。内容2)是射线法的核心和难点，它涉及到解两个非线性方程的方程组 [15] 。

$\{\begin{array}{l}h\left(i_0,j_0\right)=H\\ g\left(i_0,j_0\right)=G\end{array}$ (11)

式中： $i_0$ 为起始入射角； $j_0$ 为起始方位角；

h和g为射线到达点的坐标，它是 $i_0$ 和 $j_0$ 的函数；

H和G是期望的射线到达点的坐标。

通常不可能精确求解该线性方程组，使用迭代的方法能近似求解。在每一迭代歩中先计算偏导数：

$\frac{\partial h}{\partial i_0},\frac{\partial h}{\partial j_0},\frac{\partial g}{\partial i_0},\frac{\partial g}{\partial j_0}$ (12)

而后解线性方程组：

$\left[\begin{array}{c}\frac{\partial h}{\partial i_0}\frac{\partial h}{\partial j_0}\\ \frac{\partial g}{\partial i_0}\frac{\partial g}{\partial j_0}\end{array}\right]\left[\begin{array}{c}{i}_0^{\left(n+1\right)}-i_0^{\left(n\right)}\\ j_0^{\left(n+1\right)}-j_0^{\left(n\right)}\end{array}\right]=\left[\begin{array}{c}H-h\left(i_0^{\left(n\right)},j_0^{\left(n\right)}\right)\\ G-g\left(i_0^{\left(n\right)},j_0^{\left(n\right)}\right)\end{array}\right]$ (13)

求估计值 $i_0$ 和 $j_0$ ，反复迭代，直到( $H-h$ )和( $G-g$ )的值在期望的误差范围之内。

3. 波动方程理论

基于波动方程的算法有很多，下文着重介绍基于波动方程的交错网格高阶有限差分法。

3.1. 方法原理

设二维声波方程在直角坐标系中表示为：

$U_{xx}+U_{zz}=V^{-2}\left(x,z\right)U_{xz}$ (14)

式中：U是声波响应；V是波速

下标表示对于相应变量的偏导数。

采用显示二阶差分格式，写出其有限差分解 [16] 。

$U_{k,j}^i=2\left(1-2A^2\right)U_{k,j}^{i-1}-U_{k,j}^{i-1}+A^2\left[U_{k+1,j}^{i-1}+U_{k-1,j}^{i-1}+U_{k,j+1}^{i-1}+U_{k,j-1}^{i-1}\right]$ (15)

$A=V\left(x_k,z_j\right)\Delta t/h$ (16)

$\begin{array}{c}U\left(x_k,z_j,t_i\right)=2\left(1-2A^2\right)U\left(x_k,z_j,t_{i-1}\right)-U\left(x_k,z_j,t_{i-2}\right)\\ +A^2\left[U\left(x_{k+1},z_j,t_{i-1}\right)+U\left(x_{k-1},z_j,t_{i-1}\right)+U\left(x_k,z_{j+1},t_{i-1}\right)+U\left(x_k,z_{j-1},t_{i-1}\right)\right]\end{array}$ (17)

式中： $x_k=k\Delta x$ ； $z_j=j\Delta z$ ； $t_i=i\Delta t$ ； $h=\Delta x=\Delta z$ ； $A=\left[V\left(x_k,z_j\right)\Delta t\right]/h$ 。

递推求解时，先要给出边界条件和初始条件。在已知x-z面内起始两个时间步的波场分布后，即可递推求出x-z面内各个时间步的波场分布。

4. 应用实例

无论是射线迭代法还是交错网格高阶有限差分法正演模拟，在资料准备、地质模型建立、观测系统布设上是一致的，这样才具有可比性 [17] [18] [19] 。

4.1. 资料准备

正演要求的基本数据是：正确的深度域地面地震剖面，地质层位底界深度，地质层位的纵横波速度、密度等地球物理参数和观测系统数据。

4.2. 地质模型建立

模型建立的一般流程是：1) 建立初始模型；2) 建立构造模型；3) 定义模型的介质属性；4) 检查模型合理性，修改模型。

通过图片拓绘方式，结合地面地震剖面(见图1)，地质层位深度，建立深度域模型，然后定义速度密度等常属性及梯度属性，得到XX井地质模型(见图2)。

图1. 过XX井Line190801地震剖面

图2. 过XX井Line190801地质模型

4.3. 观测系统布设

模拟观测系统布设与实际采集时一致，观测系统见表1。

表1. XX井观测系统表

4.4. 射线追踪

射线法正演是对分块变速块内常速的二维地质模型进行共炮点射线追踪的过程。SP1点和SP2点直达波射线追踪见图3和图4。

图3. SP1点直达波反射波射线追踪

图4. SP2点直达波反射波射线追踪

4.5. 正演道集记录

制作正演道集记录时可以使用的子波有Ricker子波、Gauss子波、Yushi子波，此次采用主频30 Hz的Ricker子波。SP1和SP2点野外原始X、Z分量见图5和图6。SP1点X、Z分量正演道集记录见图7和图8。SP2点X、Z分量正演道集记录见图9和图10。

图5. SP1点野外原始X、Z分量道集记录

图6. SP2点野外原始X、Z分量道集记录

图7. SP1点X、Z分量正演道集记录(射线理论)

图8. SP1点X、Z分量正演道集记录(波动方程)

图9. SP2点X、Z分量正演道集记录(射线理论)

图10. SP2点X、Z分量正演道集记录(波动方程)

基于射线理论正演道集记录与野外资料相比，符合不太好，存在较大差别，正演道集记录整体波组特征较单一，主要是一次下行波，一次上行波和一次上行转换波，上行波同相轴有些局部还存在缺失，波组频率和分辨率与实际野外资料差别很大。

基于波动方程正演道集记录与野外资料对比，波组特征一致，符合良好，但也有一些差别：1) 因为野外资料可控震源采集时存在时移，所以，实际资料同相轴向下移；2) 野外资料所用震源比计算VSP正演道集记录所用震源能生成较长的子波，并具有高的P波和S波能量比；3) 正演道集记录上的S波与直达波的耦合在野外记录上不太明显，主要是两者震源能量比不同，而且野外衰减较正演道集记录严重；4) 野外资料弹性波在三维空间球面扩散，正演道集记录弹性波在二维空间呈柱面扩散。

总的来说，基于射线理论的正演道集记录指导意义不是太大，尤其当模型复杂时而基于波动方程的有限差分对于弹性介质和复杂构造情况下计算VSP正演道集记录最为有效。

5. 结论

通过上述对两种方法的介绍，综合考虑模型精度和成本，每一种方法都有其应用条件和应用范围。

基于射线理论的试射迭代方法制作VSP射线追踪、射线路径、覆盖次数统计、角度统计、反射(透射)系数统计等方面有很大优势，在正演道集记录制作上存在限制。适用于构造相对简单模型的正演模拟，在界面变化剧烈、速度变化剧烈的模型上(薄互层、盐丘刺穿等)，将很难得到满意的结果。

基于波动理论的有限差分方法可以解标量波动方程，也可以解矢量波动方程。既适用于均匀地层，也适应用于不均匀地层。选用的坐标系可以是笛卡尔直角坐标系，也可以是其它坐标系，可以模拟不同的观测系统。对于研究广角入射、横向变化、散射及规模与波长相当的不均匀体等有明显的优越性。制作VSP正演道集记录可用于研究某些类型的波是否存在、波的振幅变化、波至时间、同相轴曲率、能量聚焦和分散、波的互相干涉等多种特征，模拟细致的构造和岩性特征特别有利。该法主要缺点在于计算工作量大，耗费机器时间长，当网络节点多时，要求计算机内存资源多，同时该方法突出的问题之一就是数值频散问题，这种频散现象是由于波场中的各种频率成分具有不同的传播速度，从而导致这些频率成分随着时间的推移在时域波形上逐渐分散开，但目前对各种频率成分如何分散开的问题研究相对较少 [20] ，实际应用时应该相互结合这两种或其他方法，实现优势互补，更好的指导野外勘探。

文章引用

王雪峰,范增辉,刘聪伟,郭鸿喜,李 飞,李 娜,黄建华,张彦斌. 两种VSP正演模拟方法应用对比研究
A Comparative Study of Two VSP Forward Modeling Methods[J]. 地球科学前沿, 2019, 09(01): 1-11. https://doi.org/10.12677/AG.2019.91001

参考文献