高阶交错网格和PML吸收边界在横向各向同性介质地震波场模拟中的应用

引用本文

陈洁, 朱守彪. 高阶交错网格和PML吸收边界在横向各向同性介质地震波场模拟中的应用[J]. 地震地质, 2020,42(3): 654-669.
CHEN Jie, ZHU Shou-biao. THE APPLICATION OF PML BOUNDARY CODITIONS IN THE SIMULATION OF SEISMIC WAVE BY THE HIGH-ORDER STAGGERED-GRID FINITE DIFFERENCES AT THE TRANSVERSELY ISOTROPIC MEDIA[J]. SEISMOLOGY AND GEOLOGY, 2020,42(3): 654-669. 复制到剪切板

DOI:10.3969/j.issn.0253-4967.2020.03.008
Permissions

《地震地质》编辑部所有

高阶交错网格和PML吸收边界在横向各向同性介质地震波场模拟中的应用

陈洁, 朱守彪

中国地震局地壳应力研究所, 地壳动力学重点实验室, 北京 100085

作者简介：陈洁, 女, 1995年生, 2018年于东华理工大学获地球物理学专业学士学位, 现为中国地震局地壳应力研究所固体地球物理学专业在读硕士研究生, 主要从事地壳动力学数值模拟方面的研究, E-mail: celeste21ccc@gmail.com。

收稿日期: 2019-09-09

基金项目: 国家自然科学基金(41874060)、国家重点研发计划项目(2017YFC1500104)和中国地震局地壳应力研究所基本科研业务专项(ZDJ2020-15)共同资助

摘要

完全匹配层吸收边界条件通常可以很好地吸收模型边界的地震反射波, 但对于横向各向同性介质的模拟效果欠佳, 且尚在发展阶段。为此, 文中推导了横向各向同性介质中弹性动力学波动方程, 给出了施加完全匹配层(PML)吸收边界条件下时间域二阶、空间域十阶精度的高阶交错网格的有限差分形式, 并分别建立了均匀的垂直向对称轴的横向各向同性介质(VTI介质)和倾斜向对称轴的横向各向同性介质(TTI介质)模型。计算结果表明, 对于对称轴为任意角度的横向各向同性介质, 当PML边界层厚度达到一定的数值时, 可以很好地抑制人工边界所产生的地震波反射效应, 且PML的吸收效果不会被入射角与入射波频率影响。

关键词: 完全匹配层; 吸收边界条件; 横向各向同性介质; 高阶交错网格; 有限差分

中图分类号:P315.3⁺1 文献标志码:A 文章编号:0253-4967(2020)03-0654-16

THE APPLICATION OF PML BOUNDARY CODITIONS IN THE SIMULATION OF SEISMIC WAVE BY THE HIGH-ORDER STAGGERED-GRID FINITE DIFFERENCES AT THE TRANSVERSELY ISOTROPIC MEDIA

CHEN Jie, ZHU Shou-biao

Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China

Abstract

In the realm of the numerical simulation, finite difference method and finite element method are more intuitive and effective than other simulation methods. In the process of simulating seismic wave propagation, the finite differences method is widely used because of its high computational efficiency and the advantage of the algorithm is more efficient. With the demand of precision, more and more researchers have proposed more effective methods of finite differences, such as the high-order staggered-grid finite differences method, which can restore the actual process of wave propagation on the premise of ensuring accuracy and improving the efficiency of operation. In the past numerical simulation of seismic wave field, different models of isotropic medium are mostly used, but it is difficult to reflect the true layer situation. With the research demand of natural seismology and seismic exploration, the research on anisotropic media is more and more extensive. Transversely isotropic(TI)media can well simulate the seismic wave propagation in the formation medium, such as gas-bearing sandstone, mudstone, shale et al., the character of TI media is reflected by introducing the Thomsen parameters to reflect its weak anisotropy of vertical direction by using Thomson parameter. Therefore, studying the process of seismic wave propagation in TI media can restore the true information of the formation to the greatest extent, and provide a more reliable simulation basis for the numerical simulation of seismic wave propagation. In the geodynamic simulation and the numerical simulation of the seismic wave field, under the limited influence of the calculation area, if no boundary conditions are added, a strong artificial boundary reflection will be generated, which greatly reduces the validity of the simulation. In order to minimize the influence of model boundaries on the reflection of seismic waves, it is often necessary to introduce absorbing boundary conditions. At present, there are three types of absorption boundary conditions: one-way wave absorption boundary, attenuation absorption boundary, and perfectly matched layer(PML)absorption boundary. In terms of numerical simulation of seismic waves, the boundary absorption effect of PML is stronger than the first two, which is currently the most commonly used method, and it also represents the cutting-edge development direction of absorption boundary technology. The perfectly matched layer absorbing boundary is effectively applied to eliminating the reflective waves from model boundaries, but for transversely isotropic medium, the effect of the absorbing is not very well. For this reason, the elastic dynamic wave equations in transversely isotropic media are derived, and we describe a second-order accurate time, tenth-order accurate space, formulation of the Madariaga-Virieux staggered-grid finite difference methods with the perfectly matched layer(PML)are given. In addition, we have established vertical transversely isotropic(VTI)media and arbitrary inclined tilted transversely isotropic(TTI)media models, using a uniform half-space velocity model and a two-layer velocity model, respectively. By combining the actual geoscience background, we set the corresponding parameters and simulation conditions in order to make our model more research-oriented. When setting model parameters, different PML thickness, incident angle, source frequency and velocity layer models were transformed to verify the inhibition of boundary reflection effect by PML absorption boundary layer. The implementations of this simulation show that the formula is correct and for the transversely isotropic(TI)media of any angular symmetry axis, when the thickness of the PML layer reaches a certain value, the seismic wave reflection effect generated by the artificial boundary can be well suppressed, and the absorption effect of PML is not subject to changes in incident angle and wave frequency. Therefore, the results of our study indicate that our research method can be used to simulate the propagation process of seismic waves in the transversely isotropic(TI)media without being affected by the reflected waves at the model boundary to restore the actual formation information and more valuable geological research.

Keyword: perfectly matched layer; boundary conditions; transversely isotropic media; high-order staggered grid; finite difference

Show Figures

0 引言

在地球动力学及地震波场的数值模拟中, 由于计算区域有限, 若不加入一定的边界条件, 将会产生很强的人工边界反射波, 这将大幅降低模拟的有效性(张衡等, 2017; 朱守彪等, 2017, 2018)。为尽量减少模型边界对地震波反射所造成的影响, 往往需要引入吸收边界条件。

目前, 常用的吸收边界条件主要有3种类型(崔丽苹等, 2018), 具体包括:

(1)单程波吸收边界, 包括旁轴近似法(Clayton et al., 1977)、多向透射边界等。这类方法计算简单、运算效率高, 但吸收效果会随地震波入射角的变化而相应发生变化。当地震波的入射角较小时边界的吸收效果较好, 但随着入射角增大, 吸收效果将越来越差。

(2)衰减型吸收边界。这类方法基于扩展模型计算区域的思路, 在扩充区内设置阻尼因子以压制地震波振幅, 进而达到吸收地震波的目的(Cerjan et al., 1985)。但由于阻尼因子难以确定, 往往导致吸收效果不佳, 同时扩展计算区域意味着计算量的增加, 故其计算效率较低。

(3)完全匹配层(Perfectly Matched Layer, 简称PML)吸收边界(Berenger, 1994)。该方法通过在计算区域外增加匹配层, 在匹配层内使用具有衰减效应的波动方程以达到吸收边界条件的效果。该方法最早被广泛应用于电磁学领域, 之后, Collino等(2001)在推广中得到各向异性介质二维弹性波一阶速度-应力方程的PML吸收边界条件。裴正林(2004a, b)又将其推广到三维弹性波动方程交错网格高阶有限差分数值模拟中。此后, 该方法在诸多学者(单启铜等, 2007; 赵海波等, 2007; 王维红等, 2013)的研究下取得了重大突破。 PML边界的吸收效果在地震波数值模拟方面强于前2种, 是目前最常用的方法。同时, 该边界条件也代表了吸收边界技术的前沿发展方向, 诸多学者基于该方法提出了新的不同类型的吸收衰减边界条件(Komatitsch et al., 2003; 孙林洁等, 2011)。由此可见, PML吸收边界条件在数值模拟应用中具有重要作用。

随着天然地震和地震勘探领域研究的不断推进, 对于各向异性介质的研究也越来越广泛。以往的地震波场数值模拟大多采用各向同性介质模型, 但这种模型很难体现真实的地层情况。在实际地层中, 例如含气的砂岩、泥岩、页岩等都具有横向各向同性(Transversely Isotropic, TI)介质的特征(金振民等, 1994; 邓继新等, 2004), 因此研究TI介质中地震波传播的过程可最大限度还原地层的真实信息, 为地震波传播的数值模拟提供较为可靠的模拟基础(侯安宁等, 1995; 李景叶等, 2006a, b)。

截至目前, 关于地震波在横向各向同性介质内传播的数值模拟中, 针对PML吸收边界效果的研究还不够深入。本文将PML吸收边界条件用于模拟地震波在横向各向同性介质中的传播过程, 利用高阶交错网格有限差分法模拟地震波场(Virieux, 1986; Levander, 1988; Saenger et al., 2000), 采用时间域二阶、空间域十阶精度, 通过改变吸收层的厚度并构建以不同角度入射VTI介质和TTI介质的情况, 考察PML吸收边界条件在TI介质中对于不同入射角度的地震波产生的反射效应, 以检验PML吸收边界的抑制效果。

1 波动方程及方法原理

1.1 弹性波动方程及弹性参数矩阵

横向各向同性介质是一种关于垂直于横向各向同性平面的轴对称的介质, 且横向平面是具有六边形对称性的无限平面, 在垂直向通过不同物理参数体现弱各向异性。其中, 对称轴为垂直向的横向各向同性介质称为VTI介质, 对称轴倾斜的横向各向同性介质称为TTI介质(图 1)。

	Figure Option View Download New Window
	图 1 横向各向同性介质 a 垂直向对称的横向各向同性介质; b 斜向对称的横向各向同性介质Fig. 1 Transversely isotropic media.

一般地, 在二维横向各向同性介质中, 一阶速度-应力弹性波方程可表达为(Alkhalifah, 2000; Hestholm, 2009)

$\{\begin{array}{l} \frac{\partial V_{x}}{\partial t} = \frac{1}{ρ} \frac{\partial σ_{x}}{\partial x} + \frac{1}{ρ} \frac{\partial σ_{xz}}{\partial z} \\ \frac{\partial V_{z}}{\partial t} = \frac{1}{ρ} \frac{\partial σ_{xz}}{\partial x} + \frac{1}{ρ} \frac{\partial σ_{z}}{\partial z} \\ \frac{\partial σ_{x}}{\partial t} = C_{11} \frac{\partial V_{x}}{\partial x} + C_{13} \frac{\partial V_{z}}{\partial z} \\ \frac{\partial σ_{z}}{\partial t} = C_{13} \frac{\partial V_{x}}{\partial x} + C_{33} \frac{\partial V_{z}}{\partial z} \\ \frac{\partial σ_{xz}}{\partial t} = C_{44} \frac{\partial V_{x}}{\partial x} + C_{44} \frac{\partial V_{z}}{\partial z} \end{array} (1)$

在垂直向对称的横向各向同性介质VTI介质中, 弹性参数矩阵可表达为

$C_{VTI} = [\begin{array}{l} C_{11}^{0} & C_{12}^{0} & C_{13}^{0} & 0 & 0 & 0 \\ C_{21}^{0} & C_{22}^{0} & C_{23}^{0} & 0 & 0 & 0 \\ C_{31}^{0} & C_{32}^{0} & C_{33}^{0} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44}^{0} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{55}^{0} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66}^{0} \end{array}] (2)$

式中, $C_{ij}^{0}$ 为VTI介质的弹性参数。

讨论垂直向对称的横向各向同性介质的物理属性时, 通常引入Thomsen参数来描述介质的弱各向异性特征, 各参数表达式为(Thomsen, 1986; 马德堂等, 2006)

$\{\begin{array}{l} ε = \frac{C_{11}^{0} - C_{33}^{0}}{2 C_{33}^{0}} \\ δ = \frac{(C_{13}^{0} + C_{44}^{0})^{2} - (C_{33}^{0} - C_{44}^{0})^{2}}{2 C_{33}^{0} (C_{33}^{0} - C_{44}^{0})} \\ γ = \frac{C_{66}^{0} - C_{44}^{0}}{2 C_{44}^{0}} \end{array} (3)$

3个参数分别体现了P波和SV波各向异性的强弱。因此, 在VTI介质中, 各弹性参数可用Thomsen参数表达为(徐文才等, 2016)

$\{\begin{array}{l} C_{11}^{0} = C_{22}^{0} = ρ (1 + 2 ε) V_{P 0}^{2} \\ C_{33}^{0} = ρ V_{P 0}^{2} \\ C_{44}^{0} = C_{55}^{0} = ρ V_{S 0}^{2} \\ f = 1 - \frac{C_{44}^{0}}{C_{33}^{0}} \\ C_{13}^{0} = C_{23}^{0} = ρ V_{P 0}^{2} \sqrt{f (f + 2 γ)} - ρ V_{S 0}^{2} \\ C_{12}^{0} = ρ V_{P 0}^{2} (1 + 2 ε - 2 (1 - f) (1 + 2 δ)) \\ C_{66}^{0} = ρ V_{S 0}^{2} (1 + δ) \end{array} (4)$

式中, V_P₀和V_S₀分别为沿VTI介质的对称轴传播的准P波和准S波的相速度。

倾斜向的横向各向同性介质TTI介质中, 弹性参数矩阵可表达为

$C_{TTI} = [\begin{array}{l} C_{11} & C_{12} & C_{13} & 0 & C_{15} & 0 \\ C_{21} & C_{22} & C_{23} & 0 & C_{25} & 0 \\ C_{31} & C_{32} & C_{33} & 0 & C_{35} & 0 \\ 0 & 0 & 0 & C_{44} & 0 & C_{46} \\ C_{51} & C_{52} & C_{53} & 0 & C_{55} & 0 \\ 0 & 0 & 0 & C_{64} & 0 & C_{66} \end{array}] (5)$

其中, 各弹性参数与斜向旋转的角度有关, 其关系可表达为

$\{\begin{array}{l} C_{11} = si n^{2} θ (si n^{2} θ C_{33}^{0} + co s^{2} θ C_{13}^{0}) + co s^{2} θ (si n^{2} θ C_{13}^{0} + co s^{2} θ C_{11}^{0}) + si n^{2} 2 θ C_{44}^{0} \\ C_{12} = C_{21} = si n^{2} θ C_{13}^{0} + co s^{2} θ C_{12}^{0} \\ C_{13} = C_{31} = si n^{2} θ (si n^{2} θ C_{13}^{0} + co s^{2} θ C_{11}^{0}) + co s^{2} θ (si n^{2} θ C_{33}^{0} + co s^{2} θ C_{13}^{0}) - si n^{2} 2 θ C_{44}^{0} \\ C_{14} = C_{41} = 0 \\ C_{15} = C_{51} = \frac{1}{2} si n^{2} θ (si n^{2} θ C_{13}^{0} + co s^{2} θ C_{11}^{0}) - \frac{1}{2} si n^{2} θ (si n^{2} θ C_{33}^{0} + co s^{2} θ C_{13}^{0}) - si n^{2} 2 θ C_{44}^{0} (co s^{2} θ - si n^{2} θ) \\ C_{16} = C_{61} = 0 \\ C_{22} = C_{11}^{0} \\ C_{23} = C_{32} = si n^{2} 2 θ C_{12}^{0} + co s^{2} θ C_{13}^{0} \\ C_{24} = C_{42} = 0 \\ C_{25} = C_{52} = \frac{1}{2} (si n^{2} θ C_{12}^{0} - si n^{2} θ C_{13}^{0}) \\ C_{26} = C_{62} = 0 \\ C_{33} = si n^{2} θ (si n^{2} θ C_{11}^{0} + co s^{2} θ C_{13}^{0}) + co s^{2} θ (si n^{2} θ C_{13}^{0} + co s^{2} θ C_{33}^{0}) + si n^{2} 2 θ C_{44}^{0} \\ C_{34} = C_{43} = 0 \\ C_{35} = C_{53} = \frac{1}{2} \sin 2 θ (si n^{2} θ C_{11}^{0} + co s^{2} θ C_{13}^{0}) - \frac{1}{2} \sin 2 θ (si n^{2} θ C_{13}^{0} + co s^{2} θ C_{33}^{0}) + si n^{2} 2 θ C_{44}^{0} (co s^{2} θ - si n^{2} θ) \\ C_{36} = C_{63} = 0 \\ C_{44} = si n^{2} θ C_{66}^{0} + si n^{2} θ C_{44}^{0} \\ C_{45} = C_{54} = 0 \\ C_{46} = C_{64} = - \sin θ \cos θ C_{44}^{0} + \sin θ \cos θ C_{66}^{0} \\ C_{55} = \frac{1}{4} \sin 2 θ (\sin 2 θ C_{11}^{0} - \sin 2 θ C_{13}^{0}) - \frac{1}{4} \sin 2 θ (\sin 2 θ C_{13}^{0} - \sin 2 θ C_{33}^{0}) + C_{44}^{0} (co s^{2} 2 θ - si n^{2} 2 θ) \\ C_{56} = C_{65} = 0 \\ C_{66} = C_{44}^{0} si n^{2} θ + co s^{2} θ C_{66}^{0} \end{array} (6)$

θ 为偏离垂直向对称轴的旋转角度。可见, 当θ =0时即为VTI介质, 故VTI介质是TTI介质的一种特殊情况。

在横向各向同性介质中, 常用以下条件约束Thomsen参数(李磊等, 2011):

$\{\begin{array}{l} \frac{1}{4} < f < 1 \\ ε > - \frac{f}{2} \\ - \frac{1}{2} < γ < \frac{1 + 2 ε}{4 (1 - f)} - \frac{1}{2} \\ \frac{1}{2 f} - 1 < δ < 2 (\frac{1}{f} - 1) \end{array} (7)$

1.2 高阶交错网格有限差分格式及PML吸收边界条件

交错网格有限差分的原理是将速度和应力分量定义在整网格和半程网格上进行计算, 并分别对时间导数和空间导数进行泰勒展开, 得到精度为O(Δ t^M+Δ x^N)的差分格式(董良国等, 2000a, b)。

	Figure Option View Download New Window
	图 2 有限差分形式图Fig. 2 The finite-differences.

在数值模拟过程中, 由于计算机的储存有限, 计算区域会被截断, 产生人工边界反射层, 从而影响数值模拟结果的准确性。为了消除边界反射效应, 本文在数值模拟过程中采用Collino等(2001)改进后的PML吸收边界条件, 其原理是在模型计算区域的各个界面设置一个厚度可以拟定的吸收层, 对弹性波场的速度和应力分量在吸收边界设定的区域内进行分裂, 并分别对各个分裂后的波场进行不同程度的吸收和耗损(图 3)。

	Figure Option View Download New Window
	图 3 PML吸收边界条件示意图在计算区域内设置有限厚度的吸收层, 分为上、下、左、右和4个顶点范围的吸收边界Fig. 3 PML(perfectly matched layer) boundary conditions.

对不同问题进行PML吸收边界模拟计算时, 需要根据其各自的特点推导相应的PML方程(Turkel et al., 1988)。在弹性问题中, 使用PML吸收边界条件时, 往往将速度和应力分为水平和垂直2个方向(Hastings et al., 1996)。因此, 速度-应力方程可写为

$\{\begin{array}{l} v = v^{⊥} + v^{∥} \\ \frac{\partial v^{∥}}{\partial t} = D^{∥} \frac{1}{ρ} \frac{\partial σ}{\partial y} \\ \frac{\partial v^{⊥}}{\partial t} + d (x) v^{⊥} = D^{⊥} \frac{1}{ρ} \frac{\partial σ}{\partial x} \\ σ = σ^{⊥} + σ^{∥} \\ A \frac{\partial σ^{∥}}{\partial t} = E^{∥} \frac{\partial v}{\partial y} \\ A \frac{\partial σ^{⊥}}{\partial t} + d (x) σ^{⊥} = E^{⊥} \frac{\partial σ}{\partial x} \end{array} (8)$

其中,

$\begin{array}{l} D^{∥} = [\begin{array}{l} 0 & 0 & 1 \\ 0 & 1 & 0 \end{array}], D^{⊥} = [\begin{array}{l} 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}] \\ E^{∥} = [\begin{array}{l} 0 & 0 \\ 0 & 1 \\ \frac{1}{2} & 0 \end{array}], E^{⊥} = [\begin{array}{l} 1 & 0 \\ 0 & 0 \\ 0 & \frac{1}{2} \end{array}] (9) \end{array}$

A为速度的m× m矩阵; D^∥、 D^⊥分别代表平行和垂直于应力方向的基矩阵; E^∥、 E^⊥分别代表平行和垂直于速度方向的基矩阵; d(x)为阻尼因子, d(x)=d₀ ${(\frac{x}{δ_{L}})}^{2}$ , 其中x为计算过程中震源位置的坐标, d₀=log $(\frac{1}{R}) \frac{3 V_{P}}{2 δ_{L}}$ , R为反射系数, δ _L为PML的吸收层厚度, V_P为P波的最大速度。

在数值模拟过程中, 对PML的计算区域进行网格划分, 将v_x在节点(i, j)上展开, v_z在节点 $(i + \frac{1}{2}, j + \frac{1}{2})$ 展开, σ _xx和σ _zz在节点 $(i + \frac{1}{2}, j)$ 展开, σ _xz在节点 $(i, j + \frac{1}{2})$ 展开, 当PML边界层在x方向上计算时, 取d^z=0; 在z方向计算时, 取d^x=0。 PML边界条件下的高阶交错网格有限差分格式为式(10)-(14), 精度为时间域二阶、空间域2N阶, 其中 $L_{n}^{N}$ 为空间导数的阶数。

$\{\begin{array}{l} (v_{x})_{i, j}^{n} = (v_{x}^{x})_{i, j}^{n} + (v_{x}^{z})_{i, j}^{n} \\ \frac{(v_{x}^{x})_{i, j}^{n + 1} - (v_{x}^{x})_{i, j}^{n}}{Δ t} + d_{i}^{x} \frac{(v_{x}^{x})_{i, j}^{n + 1} + (v_{x}^{x})_{i, j}^{n}}{2} = \frac{1}{ρ Δ x} \overset{N}{\sum_{n = 1}} L_{n}^{N} ((σ_{xx})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} - (σ_{xx})_{i - \frac{1}{2}, j}^{n + \frac{1}{2}}) \\ \frac{(v_{x}^{z})_{i, j}^{n + 1} - (v_{x}^{z})_{i, j}^{n}}{Δ t} + d_{j}^{z} \frac{(v_{x}^{z})_{i, j}^{n + 1} + (v_{x}^{z})_{i, j}^{n}}{2} = \frac{1}{ρ Δ z} \overset{N}{\sum_{n = 1}} L_{n}^{N} ((σ_{xz})_{i, j + \frac{1}{2}}^{n + \frac{1}{2}} - (σ_{xz})_{i, j - \frac{1}{2}}^{n + \frac{1}{2}}), \end{array} (10)$

$\{\begin{array}{l} (v_{z})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n} = (v_{z}^{x})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n} + (v_{z}^{z})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n} \\ \frac{(v_{z}^{x})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n + 1} - (v_{z}^{x})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n}}{Δ t} + d_{i + \frac{1}{2}}^{x} \frac{(v_{z}^{x})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n + 1} + (v_{z}^{x})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n}}{2} = \frac{1}{ρ Δ x} \overset{N}{\sum_{n = 1}} L_{n}^{N} ((σ_{xz})_{i + 1, j + \frac{1}{2}}^{n + \frac{1}{2}} - (σ_{xz})_{i, j + \frac{1}{2}}^{n + \frac{1}{2}}) \\ \frac{(v_{z}^{z})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n + 1} - (v_{z}^{z})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n}}{Δ t} + d_{j + \frac{1}{2}}^{z} \frac{(v_{z}^{z})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n + 1} + (v_{z}^{z})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n}}{2} = \frac{1}{ρ Δ z} \overset{N}{\sum_{n = 1}} L_{n}^{N} ((σ_{zz})_{i + \frac{1}{2}, j + 1}^{n + \frac{1}{2}} - (σ_{zz})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}}) \end{array} (11)$

$\{\begin{array}{l} (σ_{x})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} = (σ_{x}^{x})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} + (σ_{x}^{z})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} \\ \frac{(σ_{x}^{x})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} - (σ_{x}^{x})_{i + \frac{1}{2}, j}^{n - \frac{1}{2}}}{Δ t} + d_{i + \frac{1}{2}}^{x} \frac{(σ_{x}^{x})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} + (σ_{x}^{x})_{i + \frac{1}{2}, j}^{n - \frac{1}{2}}}{2} = C_{11} \frac{1}{Δ x} \overset{N}{\sum_{n = 1}} L_{n}^{N} ((v_{x})_{i + 1, j}^{n} - (v_{x})_{i, j}^{n}) \\ \frac{(σ_{x}^{z})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} - (σ_{x}^{z})_{i + \frac{1}{2}, j}^{n - \frac{1}{2}}}{Δ t} + d_{j}^{z} \frac{(σ_{x}^{z})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} + (σ_{x}^{z})_{i + \frac{1}{2}, j}^{n - \frac{1}{2}}}{2} = C_{13} \frac{1}{Δ z} \overset{N}{\sum_{n = 1}} L_{n}^{N} ((v_{z})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n} - (v_{z})_{i + \frac{1}{2}, j - \frac{1}{2}}^{n}) \end{array} (12)$

$\{\begin{array}{l} (σ_{z})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} = (σ_{z}^{x})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} + (σ_{z}^{z})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} \\ \frac{(σ_{z}^{x})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} - (σ_{z}^{x})_{i + \frac{1}{2}, j}^{n - \frac{1}{2}}}{Δ t} + d_{i + \frac{1}{2}}^{x} \frac{(σ_{z}^{x})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} + (σ_{z}^{x})_{i + \frac{1}{2}, j}^{n - \frac{1}{2}}}{2} = C_{13} \frac{1}{Δ x} \overset{N}{\sum_{n = 1}} L_{n}^{N} ((v_{x})_{i + 1, j}^{n} - (v_{x})_{i, j}^{n}) \\ \frac{(σ_{z}^{z})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} - (σ_{z}^{z})_{i + \frac{1}{2}, j}^{n - \frac{1}{2}}}{Δ t} + d_{j}^{z} \frac{(σ_{z}^{z})_{i + \frac{1}{2}, j}^{n + \frac{1}{2}} + (σ_{z}^{z})_{i + \frac{1}{2}, j}^{n - \frac{1}{2}}}{2} = C_{33} \frac{1}{Δ z} \overset{N}{\sum_{n = 1}} L_{n}^{N} ((v_{z})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n} - (v_{z})_{i + \frac{1}{2}, j - \frac{1}{2}}^{n}) \end{array} (13)$

$\{\begin{array}{l} (σ_{xz})_{i, j + \frac{1}{2}}^{n + \frac{1}{2}} = (σ_{xz}^{x})_{i, j + \frac{1}{2}}^{n + \frac{1}{2}} + (σ_{xz}^{z})_{i, j + \frac{1}{2}}^{n + \frac{1}{2}} \\ \frac{(σ_{xz}^{x})_{i, j + \frac{1}{2}}^{n + \frac{1}{2}} - (σ_{xz}^{x})_{i, j + \frac{1}{2}}^{n - \frac{1}{2}}}{Δ t} + d_{i}^{x} \frac{(σ_{xz}^{x})_{i, j + \frac{1}{2}}^{n + \frac{1}{2}} + (σ_{xz}^{x})_{i, j + \frac{1}{2}}^{n - \frac{1}{2}}}{2} = C_{44} \frac{1}{Δ z} \overset{N}{\sum_{n = 1}} L_{n}^{N} ((v_{z})_{i + \frac{1}{2}, j + \frac{1}{2}}^{n} - (v_{z})_{i - \frac{1}{2}, j - \frac{1}{2}}^{n}) \\ \frac{(σ_{xz}^{z})_{i, j + \frac{1}{2}}^{n + \frac{1}{2}} - (σ_{xz}^{z})_{i, j + \frac{1}{2}}^{n - \frac{1}{2}}}{Δ t} + d_{j + \frac{1}{2}}^{z} \frac{(σ_{xz}^{z})_{i, j + \frac{1}{2}}^{n + \frac{1}{2}} + (σ_{xz}^{z})_{i, j + \frac{1}{2}}^{n - \frac{1}{2}}}{2} = C_{44} \frac{1}{Δ x} \overset{N}{\sum_{n = 1}} L_{n}^{N} ((v_{x})_{i, j + 1}^{n} - (v_{x})_{i, j}^{n}) \end{array} (14)$

2 模型实例模拟及分析

2.1 垂直向对称的横向各向同性介质(VTI介质)模型

设置模型大小为200m× 200m, 采用时间域二阶、空间域十阶精度的有限差分交错网格方法进行数值模拟。设置空间步长Δ x和Δ z均为10m, 时间步长Δ t为10ms, 选取短时间脉冲式的雷克子波为震源, 表达式为f(t)=[1-2(π f₀t)²](Ricker, 1944), 主频选为25Hz, 震源位置位于中心坐标(100m, 100m)处, 纵波速度为2.5km/s, 横波速度为1.8km/s, 密度为2i000kg/m³。由于模拟的为天然地震波, 故只考虑1个各向异性参数, 设Thomsen参数为ε =0.25、 δ =0.25、 γ =0.2(李磊等, 2011)。图 4 为VTI介质未加任何边界条件时在575ms的波场快照。

	Figure Option View Download New Window
	图 4 在575ms时VTI介质未加边界条件的波场快照图Fig. 4 Snapshots of the VTI media without boundary conditions at 575ms moment.

图 5 VTI介质加入PML边界条件的模拟结果, 展示了PML边界条件在模拟VTI介质时对边界反射效应的影响, 浅色代表经过衰减后较弱的能量, 深色则代表较强的能量。由模拟结果可见, 在VTI介质中准P波(qP)比准SV波(qSV)传播速度更快, 且能量先被吸收, 因此在快照图中的外圈为横向各向同性介质中的qP波, 而qSV波在传播时因纵向弱各向异性将被压缩成椭圆状。图5a-d为δ _L=10m时的模拟结果, 可见在加入PML边界条件后, 边界反射效应被削弱, 但上、下边界的吸收效果欠佳, 可明显见到少许反射效应; 图5e-h为δ _L=30m时的模拟结果, 与δ _L=10m时相比, 在增大吸收层厚度后, 475ms处的边界反射效应基本被完全消除。

	Figure Option View Download New Window
	图 5 均匀半空间VTI介质加入PML边界条件, 不同时刻速度波场的能量快照图Fig. 5 Snapshots of the energy of the velocity wavefields at different moments in the homogeneous acoustic VTI media within the PML boundary conditions. a-e δ _L=10m; f-j δ _L=30m

2.2 倾斜向对称的横向各向同性介质(TTI介质)模型

2.2.1 均匀半空间横向各向同性介质(TTI介质)模型

在对TTI介质模型进行模拟时, 模型大小、时间步长、空间步长保持不变。为验证PML边界条件的吸收效果是否受不同入射角度和频率的影响, 在入射角和主频分别为30° 和25Hz、 45° 和35Hz、 60° 和45Hz的情况下进行模拟。设Thomsen参数为ε =0.2、 δ =0.1、 γ =0.23, P波和S波的速度分别为2.5km/s和1.8km/s, 介质密度为2i000kg/m³, 取δ _L=30m。图 6 为TTI介质加入PML边界条件的模拟结果。

	Figure Option View Download New Window
	图 6 均匀半空间TTI介质加入PML边界条件, 不同时刻的波场快照Fig. 6 Snapshots of the homogeneous acoustic TTI media within the PML boundary conditions at different moments. a-e θ =30° , 主频为25Hz; f-j θ =45° , 主频为35Hz; k-m θ =60° , 主频为45Hz

模拟时, 给定的Thomsen参数条件中有ε ≠ δ 。 TTI介质斜向的倾斜角为上述3种角度时, 都出现了符合实际情况的qSV波前斜向三叉现象。图6a-e为入射角θ =30° 、主频为25Hz时的波场快照图, 图6f-j为入射角θ =45° 、主频为35Hz时的波场快照图, 图6k-o为入射角θ =60° 、主频为45Hz时的波场快照。对比3种入射角度和主频的结果可发现, 在模拟TTI介质中的弹性波传播时, PML吸收边界都可以很好地吸收以任意频率、角度入射的波所产生的反射波。

2.2.2 双层速度结构横向各向同性(TTI)介质模型

此外, 本研究还将PML吸收边界条件应用于横向各向同性介质(TTI介质)中的含气砂岩, 以模拟地震波的传播情况。为了提高精度并更好地还原地层信息, 模型大小、时间步长、空间步长保持不变, 设置主频为25Hz, θ =30° , 设上边界为自由边界, 左、右和底界面分别加入厚度为30m的PML吸收边界层。在Z=20m处设置一分界面, 放置一近地表震源, 震源坐标为(100m, 13m)。考虑该模型上覆沉积层砂岩、下伏含气砂岩层, 设第一层为沉积层砂岩, 其Thomsen参数为ε =0.3、 δ =0.08、 γ =0.23, 第二层为含气砂岩, ε =δ =0、 γ =0.23(Kim et al., 1993)。第一层内P波和S波的速度分别为2.8km/s和2.1km/s, 介质密度为2i100kg/m³; 第二层内P波和S波的速度分别为3.8km/s和2.8km/s, 介质密度为2i300kg/m³。双层速度结构的TTI介质模型加入PML吸收边界条件进行模拟的波场快照图见图 7。

	Figure Option View Download New Window
	图 7 双层速度结构TTI介质加入PML边界条件, 不同时刻波场快照Fig. 7 Snapshots of the TTI medium a near-surface acoustic isotropic layer of 0.02km depth within the PML boundary conditions at different moments.

由模拟结果可见, 随着时间的增加, 波前逐渐向外传播, 且因各向异性导致qP波出现近椭圆状的波前面, qSV波也出现了三叉现象; 由于模型顶部设置为自由界面, 故在模拟时产生了反射波(图7a), 当入射波传播到Z=20m分界面时, 发生了明显的反射和透射现象; 图7e、 f展示的结果说明当地震波传播到左、右边界时, 由于PML吸收边界层的作用, 并未产生明显的反射效应。可见, PML吸收边界层在双层速度结构模型中同样适用。

3 结论

本文通过对TI介质中地震波动力学方程的推导及数值计算, 得到如下初步结论:

(1)推导了横向各向同性介质中的弹性波动方程, 给出了施加完全匹配层(PML)吸收边界条件下的时间域二阶、空间域十阶精度的高阶交错网格的有限差分形式, 并通过数值计算验证所得公式是正确的。

(2)在VTI介质和TTI介质模型中模拟了不同厚度的PML吸收边界对地震波反射效应的影响。结果表明, 当吸收层厚度达到一定数值时, 该类型边界的吸收效果显著。

(3)通过建立VTI和TTI介质模型, 验证了PML吸收边界条件在任何倾斜角度的横向各向同性介质中都可以很好地抑制边界反射效应, 说明PML吸收边界条件的效果不受入射角度的影响。

(4)采用主频为25Hz、 35Hz、 45Hz的雷克子波进行TTI介质模拟。由模拟结果可见, 对于不同频率的地震波, PML吸收边界条件都可以很好地吸收边界效应产生的反射波。

(5)无论对于均匀半空间模型还是双层速度结构模型, PML吸收边界都能很好地削弱地震波在横向各向同性介质中传播时产生的边界效应, 也验证了其对于模拟不同满足地学条件的横向各向同性介质模型的适用性。

参考文献

View Option

[1]	崔丽苹, 张会星. 2018. 弹性波方程正演混合吸收边界的改进[J]. 中国海洋大学学报(自然科学版), 48(12): 87-92. CUI Li-ping, ZHANG Hui-xing. 2018. Forward modeling of elastic wave equation with the improved hybrid absorbing boundary condition[J]. Periodical of Ocean University of China, 48(12): 87-92(in Chinese). [本文引用:1]
[2]	邓继新, 史謌, 刘瑞珣, 等. 2004. 泥岩、页岩声速各向异性及其影响因素分析[J]. 地球物理学报, 47(5): 862-868. DENG Ji-xin, SHI Ge, LIU Rui-xun, et al. 2004. Analysis of the velocity anisotropy and its affection factors in shale and mudstone[J]. Chinese Journal of Geophysics, 47(5): 862-868(in Chinese). [本文引用:1]
[3]	董良国, 马在田, 曹景忠, 等. 2000 a. 一阶弹性波方程交错网格高阶差分解法[J]. 地球物理学报, 43(3): 411-419. DONG Liang-guo, MA Zai-tian, CAO Jing-zhong, et al. 2000 a. A staggered-grid high-order difference method of one-order elastic wave equation[J]. Chinese Journal of Geophysics, 43(3): 411-419(in Chinese). [本文引用:1]
[4]	董良国, 马在田, 曹景忠. 2000 b. 一阶弹性波方程交错网格高阶差分解法稳定性研究[J]. 地球物理学报, 43(6): 856-864. DONG Liang-guo, MA Zai-tian, CAO Jing-zhong. 2000 b. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation[J]. Chinese Journal of Geophysics, 43(6): 856-864(inChinese). [本文引用:1]
[5]	侯安宁, 何樵登. 1995. 各向异性介质中弹性波动高阶差分法及其稳定性的研究[J]. 地球物理学报, 38(2): 243-251. HOU An-ning, HE Qiao-deng. 1995. Study of an elastic wave high-order difference method and its stability in anisotropic media[J]. Chinese Journal of Geophysics, 38(2): 243-251(in Chinese). [本文引用:1]
[6]	金振民, Ji S, 金淑燕. 1994. 橄榄石晶格优选方位和上地幔地震波速各向异性[J]. 地球物理学报, 37(4): 469-477. JIN Zhen-min, Ji S, JIN Shu-yan. 1994. Lattice preferred orientation of olivines and seismic anisotropy in the upper mantle[J]. Chinese Journal of Geophysics, 37(4): 469-477(in Chinese). [本文引用:1]
[7]	李景叶, 陈小宏. 2006 a. TI介质地震波场数值模拟边界条件处理[J]. 西安石油大学学报(自然科学版), 21(4): 20-23. LI Jing-ye, CHEN Xiao-hong. 2006 a. Treatment of the boundary conditions in the numerical simulation of the seismic wave field in transversely isotropic(TI)medium[J]. Journal of Xi'an Shiyou University(Natural Science Edition), 21(4): 20-23(in Chinese). [本文引用:1]
[8]	李景叶, 陈小宏. 2006 b. 横向各向同性介质地震波场数值模拟研究[J]. 地球物理学进展, 21(3): 700-705. LI Jing-ye, CHEN Xiao-hong. 2006 b. Study on seismic wave field numerical simulation in transverse isotropic medium[J]. Progress in Geophysics, 21(3): 700-705(in Chinese). [本文引用:1]
[9]	李磊, 郝重涛. 2011. 横向各向同性介质和斜方介质各向异性参数的约束条件[J]. 地球物理学报, 54(11): 2819-2830. LI Lei, HAO Chong-tao. 2011. Constraints on anisotropic parameters in transversely isotropic media and the extensions to orthorhombic media[J]. Chinese Journal of Geophysics, 54(11): 2819-2830(in Chinese). [本文引用:1]
[10]	马德堂, 朱光明. 2006. 关于横向各向同性介质中的Thomsen参数取值的讨论[J]. 石油地球物理勘探, 41(4): 431-438. MA De-tang, ZHU Guang-ming. 2006. Discussion on taking values of Thomsen parameters in transverse isotropic media[J]. Oil Geophysical Prospecting, 41(4): 431-438(in Chinese). [本文引用:1]
[11]	裴正林. 2004 a. 三维各向异性介质中弹性波方程交错网格高阶有限差分法数值模拟[J]. 石油大学学报(自然科学版), 28(5): 23-29. PEI Zheng-lin. 2004 a. Three-dimensional numerical simulation of elastic wave propagation in 3-D anisotropic media with staggered-grid high-order difference method[J]. Journal of China University of Petroleum(Natural Science Edition), 28(5): 23-29(in Chinese). [本文引用:1]
[12]	裴正林. 2004 b. 任意起伏地表弹性波方程交错网格高阶有限差分法数值模拟[J]. 石油地球物理勘探, 39(6): 629-634. PEI Zheng-lin. 2004 b. Staggered-grid high-order finite difference numerical simulation of elastic wave equation with arbitrary undulations on the surface[J]. Oil Geophysical Prospecting, 39(6): 629-634. [本文引用:1]
[13]	单启铜, 乐友喜. 2007. PML边界条件下二维黏弹性介质波场模拟[J]. 石油物探, 46(2): 126-130. SHAN Qi-tong, YUE You-xi. 2007. Wavefield simulation of 2-D viscoelastic medium in perfectly matched layer boundary[J]. Geophysical Prospecting for Petroleum, 46(2): 126-130(in Chinese). [本文引用:1]
[14]	孙林洁, 印兴耀. 2011. 基于PML边界条件的高倍可变网格有限差分数值模拟方法[J]. 地球物理学报, 54(6): 1614-1623. SUN Lin-jie, YIN Xing-yao. 2011. A finite-difference scheme based on PML boundary condition with high power grid step variation[J]. Chinese Journal of Geophysics, 54(6): 1614-1623(in Chinese). [本文引用:1]
[15]	王维红, 柯璇, 裴江云. 2013. 完全匹配层吸收边界条件应用研究[J]. 地球物理学进展, 28(5): 2508-2514. WANG Wei-hong, KE Xuan, PEI Jiang-yun. 2013. Application investigation of perfectly matched layer absorbing boundary condition[J]. Progress in Geophysics, 28(5): 2508-2514(in Chinese). [本文引用:1]
[16]	徐文才, 杨国权, 李振春, 等. 2016. 横向各向同性介质拟声波一阶速度-应力方程[J]. 石油地球物理勘探, 51(1): 87-96. XU Wen-cai, YANG Guo-quan, LI Zhen-chun, et al. 2016. First order velocity-stress equation in TI media[J]. Oil Geophysical Prospecting, 51(1): 87-96(in Chinese). [本文引用:1]
[17]	张衡, 刘洪, 李博, 等. 2017. TTI介质声波方程分裂式PML吸收边界条件研究[J]. 石油物探, 56(3): 349-361. ZHANG Heng, LIU Hong, LI Bo, et al. 2017. The research on split PML absorbing boundary conditions of acoustic equation for TTI media[J]. Geophysical prospecting for Petroleum, 56(3): 349-361(in Chinese). [本文引用:2]
[18]	赵海波, 王秀明, 王东, 等. 2007. 完全匹配层吸收边界在孔隙介质弹性波模拟中的应用[J]. 地球物理学报, 50(2): 581-591. ZHAO Hai-bo, WANG Xiu-ming, WANG Dong, et al. 2007. Applications of the boundary absorption using a perfectly matched layer for elastic wave simulation in poroelastic media[J]. Chinese Journal of Geophysics, 50(2): 581-591(in Chinese). [本文引用:1]
[19]	朱守彪, 袁杰, 缪淼. 2017. 青海玉树地震( M_S=7. 1)产生超剪切破裂过程的动力学机制研究[J]. 地球物理学报, 60(10): 3832-3843. ZHU Shou-biao, YUAN Jie, MIAO Miao. 2017. Dynamic mechanisms for supershear rupture processes of the Yushu earthquake( M_S=7. 1)[J]. Chinese Journal of Geophysics, 60(10): 3832-3843(in Chinese). [本文引用:1]
[20]	朱守彪, 袁杰. 2018. 2008年汶川大地震中北川地区极重震害的物理机制研究[J]. 地球物理学报, 61(5): 1863-1873. ZHU Shou-biao, YUAN Jie. 2018. Physical mechanism for extremely serious seismic damage in the Beichuan area caused by the great 2008 Wenchuan earthquake[J]. Chinese Journal of Geophysics, 61(5): 1863-1873(in Chinese). [本文引用:1]
[21]	Alkhalifah T. 2000. An acoustic wave equation for anisotropic media[J]. Geophysics, 65(4): 1239-1250. [本文引用:1]
[22]	Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of computational physics, 114(2): 185-200. [本文引用:1]
[23]	Cerjan C, Kosloff D, Kosloff R, et al. 1985. A nonreflecting boundary condition for discrete acoustic and elastic wave equations[J]. Geophysics, 50(4): 705-708. [本文引用:1]
[24]	Clayton R W, Engquist B. 1977. Absorbing boundary conditions for acoustic and elastic wave equations[J]. Bulletin of the Seismological Society of America, 67(6): 1529-1540. [本文引用:1]
[25]	Collino F, Tsogka C. 2001. Application of the perfectly matched absorbing layer model to the linear electrodynamic problem in anisotropic heterogeneous media[J]. Geophysics, 66(1): 294-307. [本文引用:2]
[26]	Hastings F D, Schneider J B, Broschat S L. 1996. Application of the perfectly matched layer(PML)absorbing boundary condition to elastic wave propagation[J]. The Journal of the Acoustical Society of America, 100(5): 3061-3069. [本文引用:1]
[27]	Hestholm S. 2009. Acoustic VTI modeling using high-order finite differences[J]. Geophysics, 74(5): 67-73. [本文引用:1]
[28]	Kim K Y, Wrolstad K H, Aminzadeh F. 1993. Effects of transverse isotropy on P-wave AVO for gas sand s[J]. Geophysics, 58(6): 883-888. [本文引用:1]
[29]	Komatitsch D, Tromp J. 2003. A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation[J]. Geophysical Journal International, 154(1): 146-153. [本文引用:1]
[30]	Levand er A R. 1988. Fourth-order finite-difference P-SV seismograms[J]. Geophysics, 53(11): 1425-1436. [本文引用:1]
[31]	Ricker N. 1944. Wavelet functions and their polynomials[J]. Geophysics, 9(3): 314-323. [本文引用:1]
[32]	Saenger E H, Gold N, Shapiro S A. 2000. Modeling the propagation of elastic waves using a modified finite-difference grid[J]. Wave Motion, 31(1): 77-92. [本文引用:1]
[33]	Thomsen I. 1986. Weak elastic anisotropy[J]. Geophysics, 51(10): 1954-1966. [本文引用:1]
[34]	Turkel E, Yefet A. 1998. Absorbing PML boundary layers for wave-like equation[J]. Applied Numerical Mathematics, 27(4): 533-557. [本文引用:1]
[35]	Virieux J. 1986. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method[J]. Geophysics, 51(4): 889-901. [本文引用:]

2018

0.0

崔丽苹, 张会星. 2018. 弹性波方程正演混合吸收边界的改进[J]. 中国海洋大学学报(自然科学版), 48(12): 87-92.

CUI

Li-ping

, ZHANG

Hui-xing

. 2018. Forward modeling of elastic wave equation with the improved hybrid absorbing boundary condition[J]. Periodical of Ocean University of China, 48(12): 87-92(in Chinese).

... 目前, 常用的吸收边界条件主要有3种类型(崔丽苹等, 2018), 具体包括: ...

2004

0.0

邓继新, 史謌, 刘瑞珣, 等. 2004. 泥岩、页岩声速各向异性及其影响因素分析[J]. 地球物理学报, 47(5): 862-868.

DENG

Ji-xin

, SHI

, LIU

Rui-xun

, et al. 2004. Analysis of the velocity anisotropy and its affection factors in shale and mudstone[J]. Chinese Journal of Geophysics, 47(5): 862-868(in Chinese).

2000

0.0

董良国, 马在田, 曹景忠, 等. 2000 a. 一阶弹性波方程交错网格高阶差分解法[J]. 地球物理学报, 43(3): 411-419.

DONG

Liang-guo

, MA

Zai-tian

, CAO

Jing-zhong

, et al. 2000 a. A staggered-grid high-order difference method of one-order elastic wave equation[J]. Chinese Journal of Geophysics, 43(3): 411-419(in Chinese).

All methods of seismic wave-field simulation are trying to improve their computational accuracy and efficiency. After transforming the odd higher-order time derivatives of velocity （stress） to spatial derivatives of stress （velocity） in the high-order finite difference method, we can use any order of the temporal and spatial difference accuracy and staggered-grid technique for numerical calculation of the one-order elastic wave equations expressed with velocity and stress. The snapshots and simulated results of an actual model show that this method is more accurate and can decrease the grid dispersion in conventional difference method. Meanwhile, bigger grid space can be used in order to raise the computational efficiency.

Laboratory of Marine Geology of NDE, Tongji University, Shanghai 200092, China

提高计算精度和运算效率是所有波场正演方法所追求的目标，本文通过将速度（应力）对时间的奇数阶高阶寻数转化为应力（速度）对空间的导数，运用时间和空间差分精度均可达任意阶的高阶差分法，通过交错网格技术，对一阶速度－应力弹性波方程进行了数值求解．波场快照以及实际模型的正演结果表明，这种求解一阶弹性波方程的高阶差分解法，和常规的差分法相比网格频散显著减小，精度明显提高，而且可以取较大的空间步长，提高计算效率。

... 邓继新等, 2004), 因此研究TI介质中地震波传播的过程可最大限度还原地层的真实信息, 为地震波传播的数值模拟提供较为可靠的模拟基础(侯安宁等, 1995 ...

2000

0.0

董良国, 马在田, 曹景忠. 2000 b. 一阶弹性波方程交错网格高阶差分解法稳定性研究[J]. 地球物理学报, 43(6): 856-864.

DONG

Liang-guo

, MA

Zai-tian

, CAO

Jing-zhong

. 2000 b. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation[J]. Chinese Journal of Geophysics, 43(6): 856-864(inChinese).

1995

0.0

侯安宁, 何樵登. 1995. 各向异性介质中弹性波动高阶差分法及其稳定性的研究[J]. 地球物理学报, 38(2): 243-251.

HOU

An-ning

, HE

Qiao-deng

. 1995. Study of an elastic wave high-order difference method and its stability in anisotropic media[J]. Chinese Journal of Geophysics, 38(2): 243-251(in Chinese).

1994

0.0

金振民, Ji

, 金淑燕. 1994. 橄榄石晶格优选方位和上地幔地震波速各向异性[J]. 地球物理学报, 37(4): 469-477.

JIN

Zhen-min

, Ji

, JIN

Shu-yan

. 1994. Lattice preferred orientation of olivines and seismic anisotropy in the upper mantle[J]. Chinese Journal of Geophysics, 37(4): 469-477(in Chinese).

Olivine-petrofabric data in combination with crystalelastic constants have beenused to calculate velocity anisotropy from mantle-derived xenoliths（1herzolite） inMingxi, Fujian province.It is suggested that seismic anisotropy originates from lattice preferred orientation（LPO） of olivine induced by plastic flow of highly anisotropic olivines. The Patterns of olivine fabric and distribution of Vp, VSV, VSH and△VS with respect to the tectonic setting implies the direction of mantle-matter flowfrom NWW to SEE, which is coincided with the orientation of a-[100] axies（maximum VP direction） of olivine beneath East-South coast region in China.

1. China University of Geosciences Wuhan 430074; 2. Department of Geology, University of Montreal, Quebec, Canada

根据福建省明溪幔源包体（二辉橄榄岩）中橄榄石晶格优选方位（ＬＰＯ）及其晶体弹性刚度系数，计算了地震波速度及其各向异性．研究结果表明，该区地震波各向异性是由橄榄石塑性流动产生晶格优选方位而引起的．与构造背景有关的ＶＰ，Ｖｓ１，Ｖｓ２和△Ｖｓ分布特征表明，中国东南沿海地区上地幔物质流动方向（由ＮＷＷ向ＳＥＥ）与橄榄石［１００］定向排列方向（ａ轴）和ＶＰ最大方向有一致的趋势．

2006

0.0

李景叶, 陈小宏. 2006 a. TI介质地震波场数值模拟边界条件处理[J]. 西安石油大学学报(自然科学版), 21(4): 20-23.

Jing-ye

, CHEN

Xiao-hong

. 2006 a. Treatment of the boundary conditions in the numerical simulation of the seismic wave field in transversely isotropic(TI)medium[J]. Journal of Xi'an Shiyou University(Natural Science Edition), 21(4): 20-23(in Chinese).

... 李景叶等, 2006a, b) ...

2006

0.0

李景叶, 陈小宏. 2006 b. 横向各向同性介质地震波场数值模拟研究[J]. 地球物理学进展, 21(3): 700-705.

Jing-ye

, CHEN

Xiao-hong

. 2006 b. Study on seismic wave field numerical simulation in transverse isotropic medium[J]. Progress in Geophysics, 21(3): 700-705(in Chinese).

Seismic wave field numerical simulation is one of effective methods for understanding seismic wave characterization during spreading in underground medium, and for helping to interpret the recording data. Improving computation precision and efficiency is the main objective of wave field numerical simulation studies. As one of the most efficient methods, the finite difference technique is widely used in wave equation solution. But the shortcoming of the traditional numerical frequency dispersion during wave equation finite difference solvingt affects the calculation efficiency and precision. In the paper, one-order velocity-stress elastic wave equations in transverse isotropic medium are solved by a staggered-grid high-order difference flux-corrected transport method. The snapshots computed by different difference methods show that the method study in the paper work more favorable in decreasing grid dispersion and improving computation precision compared with two-order and staggered-grid four-order difference methods. Therefore the longer spatial step length can be utilized to improve computation efficiency.

地震波场数值模拟是理解地震波在地下介质中的传播特点，帮助解释观测数据的有效手段，而提高计算精度和运算效率是所有波场数值模拟方法研究所追求的目标.有限差分技术是求解波动方程计算效率最高、应用最为广泛的方法之一.但传统的有限差分技术计算过程中的数值频散问题影响了该技术的计算精度与计算效率.本文通过交错网格高阶有限差分技术与通量校正传输方法（Flux|corrected transport method，FCT）相结合, 对横向各向同性介质（Transverse isotropic medium，TI）一阶速度|应力弹性波动方程组进行了数值求解研究.波场快照数值模拟结果表明，本文研究的数值模拟方法与波动方程二阶有限差分方法、交错网格四阶有限差分方法相比，在压制网格数值频散方面有明显的优势，计算精度提高，而且可以利用较大的空间步长，提高计算效率.

... 李景叶等, 2006a, b) ...

2011

0.0

李磊, 郝重涛. 2011. 横向各向同性介质和斜方介质各向异性参数的约束条件[J]. 地球物理学报, 54(11): 2819-2830.

Lei

, HAO

Chong-tao

. 2011. Constraints on anisotropic parameters in transversely isotropic media and the extensions to orthorhombic media[J]. Chinese Journal of Geophysics, 54(11): 2819-2830(in Chinese).

The presence of anisotropy in earth materials has been widely accepted. Due to its high symmetry and relative simplicity in mathematical formulae, transversely isotropic (TI) medium has become one of the most important anisotropic media in seismology. In theoretical and practical studies, anisotropic parameters expressed in terms of elastic constants proposed by Thomsen in 1986 are widely adopted. Values of elastic constants and Thomsen's anisotropic parameters cannot be chosen freely. There are some constraints among those parameters. Unfortunately, such constraints were ignored by some authors. Sometimes improper parameters were selected and consequently physically meaningless and misleading results were generated. In this study, we proposed a set of simple inequalities for constraining the values of Thomsen's anisotropic parameters. The inequalities, given by 1/4 - f/2, 1/2f-1 = max(0, delta) and gamma >= 0. For some directions called longitudinal directions, at which phase velocity of P waves reaches an extremum, the direction of P polarization is coincident with wavenormal and ray direction. The sufficient and necessary condition for the presence of longitudinal direction is given by epsilon/delta(L) (1- f)/2 is the condition for basal triplication, and delta - epsilon >(1 - f)/2 for axial triplication, (epsilon-delta)/(1+delta)> 2(1 - f)/3 for oblique triplication. The inequality for oblique triplication is simplified from the inequality proposed by Thomsen and Dellinger (2003). We also extend the constraints on elastic constants and anisotropic parameters of TI media to those of orthorhombic media.

[Li Lei; Hao Chong-Tao] China Earthquake Adm, Inst Geol, State Key Lab Earthquake Dynam, Beijing 100029, Peoples R China.

Thomsen提出的横向各向同性(TI)介质各向异性参数( α 0 、 β 0 、 ε、δ、γ )是各向异性理论研究和实际资料处理中的常用参数,Thomsen参数的取值必须符合物理学定律和实际地学情况,随意的取值可能导致无意义乃至错误的结果.本文根据热力学定律和弹性常数的物理意义,结合大量的实测数据,提出常见TI介质的Thomsen参数需要满足以下约束条件:1/4 f =1- β 0 2 / α 0 2 ε >- f /2; 1/2 f -1 δ f -1); -1/2 γ ε )/4(1- f )-1/2.Thomsen参数 ε、δ、γ 的取值区间主要受P、S波参考速度比 β 0 / α 0 的约束,其值可正可负,实测参数中 ε和γ 正多负少; δ、γ、β 0 /α 0 的取值范围都有上下界,而 ε没有上限;γ与ε 正相关且取值上限受到 ε 的限制; ε、γ 与 δ 之间不存在明显的约束.鉴于椭圆各向异性介质和薄互层等效TI介质在实际应用中的普遍性,专门给出了这两种介质各向异性参数的附加约束条件.一些在实际资料中可观测到的P波特殊偏振方向和SV波三叉现象也能为各向异性参数提供额外的约束.之后,我们将TI介质弹性常数和各向异性参数的约束条件扩展到对正交各向异性介质弹性常数和各向异性参数的约束.本文提出的各向异性参数约束条件简单实用,为各向异性理论研究和数值模拟中参数的选择提供依据,为实际资料的各向异性参数反演提供约束,既能避免出现无物理意义的研究结果,又能加速反演的搜索过程,提高生产效率,具有理论指导意义和实际应用价值.

2006

0.0

马德堂, 朱光明. 2006. 关于横向各向同性介质中的Thomsen参数取值的讨论[J]. 石油地球物理勘探, 41(4): 431-438.

De-tang

, ZHU

Guang-ming

. 2006. Discussion on taking values of Thomsen parameters in transverse isotropic media[J]. Oil Geophysical Prospecting, 41(4): 431-438(in Chinese).

2004

0.0

裴正林. 2004 a. 三维各向异性介质中弹性波方程交错网格高阶有限差分法数值模拟[J]. 石油大学学报(自然科学版), 28(5): 23-29.

PEI

Zheng-lin

. 2004 a. Three-dimensional numerical simulation of elastic wave propagation in 3-D anisotropic media with staggered-grid high-order difference method[J]. Journal of China University of Petroleum(Natural Science Edition), 28(5): 23-29(in Chinese).

... 在实际地层中, 例如含气的砂岩、泥岩、页岩等都具有横向各向同性(Transversely Isotropic, TI)介质的特征(金振民等, 1994 ...

2004

0.0

裴正林. 2004 b. 任意起伏地表弹性波方程交错网格高阶有限差分法数值模拟[J]. 石油地球物理勘探, 39(6): 629-634.

PEI

Zheng-lin

. 2004 b. Staggered-grid high-order finite difference numerical simulation of elastic wave equation with arbitrary undulations on the surface[J]. Oil Geophysical Prospecting, 39(6): 629-634.

2007

0.0

单启铜, 乐友喜. 2007. PML边界条件下二维黏弹性介质波场模拟[J]. 石油物探, 46(2): 126-130.

SHAN

Qi-tong

, YUE

You-xi

. 2007. Wavefield simulation of 2-D viscoelastic medium in perfectly matched layer boundary[J]. Geophysical Prospecting for Petroleum, 46(2): 126-130(in Chinese).

2011

0.0

孙林洁, 印兴耀. 2011. 基于PML边界条件的高倍可变网格有限差分数值模拟方法[J]. 地球物理学报, 54(6): 1614-1623.

SUN

Lin-jie

, YIN

Xing-yao

. 2011. A finite-difference scheme based on PML boundary condition with high power grid step variation[J]. Chinese Journal of Geophysics, 54(6): 1614-1623(in Chinese).

Heterogeneous media is a kind of multi-scale media. The scale of the fracture reservoir could reach up to millimeter. Simulating a fractured carbonate model with finite-difference method is calling for the same scale of sample steps. A finite-difference method with high power of grid steps variation is presented here. Spatial and temporal sampling steps can be optimized according to material properties and model geometry. The ratio of grid steps could reach up to hundreds or more. We improve the method of numerical grid division in order to add more flexibility and efficiency to finite-difference method with grid variation. The numerical results support that the method can describe heterogeneous media in scale of millimeter accurately with efficiency, and the property made the method suitable for modeling cavity or fractured carbonate reservoirs.

[Sun Lin-Jie; Yin Xing-Yao] China Univ Petr, Sch Geosci, Qingdao 266555, Peoples R China.

非均匀介质一般为多尺度介质,对非均匀介质进行正演模拟需要多尺度的网格剖分.碳酸盐岩缝洞介质的尺度一般为厘米级甚至毫米级.使用有限差分方法对其进行精细模拟需要差分步长达到缝洞介质的尺度.为了提高有限差分数值模拟方法的精度和效率,使之可以应用于非均匀介质的正演模拟,本文推导了基于PML边界的空间和时间步长同时变化的高倍数可变网格差分格式,步长变化倍数可以达到百倍以上.并且在一般意义的变网格算法的基础上,改进了变网格算法的网格剖分方式,进一步减小了精细尺度模型数值模拟的内存消耗.数值试验表明,该方法可以精细描述毫米尺度的地质体,提高有限差分方法模拟精度,同时也节约了内存,提高了模拟效率.

2013

0.0

王维红, 柯璇, 裴江云. 2013. 完全匹配层吸收边界条件应用研究[J]. 地球物理学进展, 28(5): 2508-2514.

WANG

Wei-hong

, KE

Xuan

, PEI

Jiang-yun

. 2013. Application investigation of perfectly matched layer absorbing boundary condition[J]. Progress in Geophysics, 28(5): 2508-2514(in Chinese).

2016

0.0

徐文才, 杨国权, 李振春, 等. 2016. 横向各向同性介质拟声波一阶速度-应力方程[J]. 石油地球物理勘探, 51(1): 87-96.

Wen-cai

, YANG

Guo-quan

, LI

Zhen-chun

, et al. 2016. First order velocity-stress equation in TI media[J]. Oil Geophysical Prospecting, 51(1): 87-96(in Chinese).

2017

0.0

... 在横向各向同性介质中, 常用以下条件约束Thomsen参数(李磊等, 2011): ...

... 2(李磊等, 2011) ...

2007

0.0

赵海波, 王秀明, 王东, 等. 2007. 完全匹配层吸收边界在孔隙介质弹性波模拟中的应用[J]. 地球物理学报, 50(2): 581-591.

ZHAO

Hai-bo

, WANG

Xiu-ming

, WANG

Dong

, et al. 2007. Applications of the boundary absorption using a perfectly matched layer for elastic wave simulation in poroelastic media[J]. Chinese Journal of Geophysics, 50(2): 581-591(in Chinese).

2017

0.0

... 马德堂等, 2006) ...

2018

0.0

朱守彪, 袁杰. 2018. 2008年汶川大地震中北川地区极重震害的物理机制研究[J]. 地球物理学报, 61(5): 1863-1873.

ZHU

Shou-biao

, YUAN

Jie

. 2018. Physical mechanism for extremely serious seismic damage in the Beichuan area caused by the great 2008 Wenchuan earthquake[J]. Chinese Journal of Geophysics, 61(5): 1863-1873(in Chinese).

2000

0.0

1994

0.0

1985

0.0

1977

0.0

2001

0.0

... 此后, 该方法在诸多学者(单启铜等, 2007 ...

... 为了消除边界反射效应, 本文在数值模拟过程中采用Collino等(2001)改进后的PML吸收边界条件, 其原理是在模型计算区域的各个界面设置一个厚度可以拟定的吸收层, 对弹性波场的速度和应力分量在吸收边界设定的区域内进行分裂, 并分别对各个分裂后的波场进行不同程度的吸收和耗损(图 3) ...

1996

0.0

2009

0.0

... 孙林洁等, 2011) ...

1993

0.0

2003

0.0

... 王维红等, 2013)的研究下取得了重大突破 ...

1988

0.0

1944

0.0

... 因此, 在VTI介质中, 各弹性参数可用Thomsen参数表达为(徐文才等, 2016) ...

2000

0.0

1986

0.0

... 0 引言在地球动力学及地震波场的数值模拟中, 由于计算区域有限, 若不加入一定的边界条件, 将会产生很强的人工边界反射波, 这将大幅降低模拟的有效性(张衡等, 2017 ...

1998

0.0

... 对不同问题进行PML吸收边界模拟计算时, 需要根据其各自的特点推导相应的PML方程(Turkel et al ...

1986