物探与化探 2018 , 42 (1): 134-143 https://doi.org/10.11720/wtyht.2018.1.16

Orginal Article

TTI介质拟声波方程数值模拟

黄杰, 杨国权, 李振春, 谷丙洛

中国石油大学(华东) 地球科学与技术学院,山东青岛 266580

Numerical simulation of pseudo acoustic wave equation in TTI medium

HUANG Jie, YANG Guo-Quan, LI Zhen-Chun, GU Bing-Luo

School of Geosciences,China University of Petroleum (East China),Qingdao 266580,China

中图分类号: P631.4

文献标识码: A

文章编号: 1000-8918(2018)01-0134-10

责任编辑: HUANG JieYANG Guo-QuanLI Zhen-ChunGU Bing-Luo

收稿日期: 2016-11-18

修回日期: 2017-10-23

网络出版日期: 2018-01-20

基金资助: 国家重点基础研究发展计划(“973”计划)课题(2014CB239006),国家自然科学基金项目(41504100)

作者简介:

作者简介: 黄杰(1991-),男,硕士研究生,主要从事地震波正演模拟与逆时偏移工作。Email:18205420120@163.com

展开

摘要

TTI介质qP波数值模拟方法因为考虑了倾角因素,可以比VTI介质qP波数值模拟方法更加准确地描述各向异性介质中地震波场的传播规律。文中用拟声波方程对TTI介质中的地震波场进行了高阶有限差分数值模拟,在改进衰减函数分布方式后,通过坐标变换,利用改进的完全匹配层(perfectly matched layer,PML)边界控制方程对波场边界进行吸收处理,取得了良好的效果;然后分析了拟声波方程数值模拟中的稳定性问题,并对波场中的伪横波进行压制。通过对不同模型的数值模拟,验证了文中使用的TTI介质拟声波波动方程的稳定性以及所采用的PML边界控制方程的可靠性和适用性。

关键词： TTI介质 ; 拟声波方程 ; PML边界控制方程 ; 稳定性 ; 伪横波

Abstract

TTI qP wave numerical simulation method is more accurate than VTI qP wave numerical simulation method in describing the propagation law of seismic wavefield in anisotropic medium because of considering the dip angle factor.Using the pseudo acoustic wave equation,the authors carried out high order finite difference numerical simulation of the seismic wavefield in the TTI medium in this paper.After improving the distribution of the attenuation function,the authors used the improved perfectly matched layer (PML) boundary control equation to deal with the wavefield boundary and achieved good results.Then the pseudo shear wave in the numerical simulation of pseudo acoustic wave equation was suppressed and the stability problem was analyzed.The numerical simulation of different models prove that the TTI medium pseudo acoustic wave equation is stable and the PML boundary control equations have high reliability and applicability.

Keywords： TTI medium ; pseudo acoustic wave equation ; PML boundary control equation ; stability ; pseudo shear wave

PDF (2200KB) 元数据多维度评价相关文章收藏文章

本文引用格式导出 EndNote Ris Bibtex

黄杰, 杨国权, 李振春, 谷丙洛. TTI介质拟声波方程数值模拟[J]. 物探与化探, 2018, 42(1): 134-143 https://doi.org/10.11720/wtyht.2018.1.16

HUANG Jie, YANG Guo-Quan, LI Zhen-Chun, GU Bing-Luo. Numerical simulation of pseudo acoustic wave equation in TTI medium[J]. , 2018, 42(1): 134-143 https://doi.org/10.11720/wtyht.2018.1.16

0 引言

地震波场数值模拟是研究地震波在地下介质中传播规律的有效手段;有限差分法由于其计算速度快、精度高的优点,成为目前应用最为广泛的方法^[1]。

大量的研究表明,各向异性广泛存在于地下介质中。基于各向同性的假设无法准确地描述地震波在地下实际介质中的传播规律。前人对各向同性的研究,无论是理论上还是实际应用都发展得相当成熟了,考虑了各向异性的因素,对地震勘探来说有着十分重要的作用。若只是考虑各向同性或基于VTI介质的假设,在对盐丘侧翼、盐下等的复杂构造进行数值模拟时,就会产生运动学和动力学上的误差,并最终影响偏移成像的效果。在TTI介质中,其对称轴方向与垂直方向已不再一致;但从物理上来说,TTI介质与VTI介质并没有本质区别,只不过是观测坐标系与描述介质对称性的物理坐标系有一定夹角而已。

自从Alkhalifah^[2-3]在VTI介质研究中,从频散关系出发,提出了声学近似的思想(即假设沿对称轴方向的横波速度为零),并推导出VTI介质四阶拟声波偏微分方程,许多学者对VTI介质拟声波方程正演模拟及逆时偏移进行了更深入的研究,并逐渐向TTI介质发展。但由于四阶偏微分方程求解过程十分复杂,Zhou等^[4]、Du等^[5]和Hestholm^[6]推导了VTI介质2阶耦合的拟声波方程。Zhou等^[7]、Fletcher等^[8]基于声学近似的TTI介质频散关系,通过引入辅助波场函数,推导出不同的2阶耦合的拟声波方程,但辅助函数的物理意义并不明确。Duveneck^[9-10]从胡克定律和运动方程出发,推导出VTI介质拟声波波动方程,并赋予水平和垂直应力分量明确的物理意义,而后又推广到TTI介质。2010年,Fowler等^[11]给出了二阶拟声波偏微分方程的一般形式,并认为它们彼此等价,可以相互转化。但几种拟声波方程都存在低速、低振幅的qSV波干扰问题,这会对正演模拟和偏移成像的结果造成影响。针对这一问题,Liu^[12]、Pestana等^[13]、Chu等^[14]、Zhan等^[15]、Barrera等^[16]进行了较为深入的研究。

笔者基于二阶TTI介质拟声波方程,利用改进的PML边界控制方程对边界进行吸收处理,并进行了模型试算,数值模拟结果证明了拟声波方程的稳定性和改进的边界条件的有效性。

1 稳定的TTI介质拟声波方程

从Christoffel方程出发,可以得到精确的TTI介质qP-qSV波频散关系^[8]为:

$\begin{matrix} \begin{matrix} ω^{4} = [(v_{px}^{2} + v_{sz}^{2}) T_{2} + (v_{pz}^{2} + v_{sz}^{2}) T_{1}] ω^{2} + v_{px}^{2} v_{sz}^{2} T_{2}^{2} + \\ [v_{pz}^{2} (v_{pn}^{2} - v_{px}^{2}) + v_{sz}^{2} (v_{pn}^{2} + v_{pz}^{2})] T_{1} T_{2} + v_{pz}^{2} v_{sz}^{2} T_{1}^{2}, (1) \end{matrix} \end{matrix}$

其中:

$\begin{matrix} \{\begin{matrix} T_{1} = k_{x}^{2} si n^{2} θ + k_{z}^{2} co s^{2} θ + k_{x} k_{z} \sin 2 θ, \\ T_{2} = k_{x}^{2} co s^{2} θ + k_{z}^{2} si n^{2} θ - k_{x} k_{z} \sin 2 θ 。 \end{matrix} (2) \end{matrix}$

v_pn=v_pz $\begin{matrix} \sqrt[]{1 + 2 δ} \end{matrix}$ 、v_px=v_pz $\begin{matrix} \sqrt[]{1 + 2 ε} \end{matrix}$ 、v_pz分别为NMO速度、对称平面的切线速度、关于对称平面的法线速度,θ为对称轴的倾角,ε和δ为Thomsen各向异性参数^[17],ω为角频率,k_x、k_z为空间波数。

在 $\begin{matrix} v_{sz}^{2} \end{matrix}$ =0时,式(1)简化为:

$\begin{matrix} ω^{4} = (v_{px}^{2} T_{2} + v_{pz}^{2} T_{1}) ω^{2} + v_{pz}^{2} (v_{pn}^{2} - v_{px}^{2}) T_{1} T_{2} 。 (3) \end{matrix}$

基于以上关系,通过引入一个辅助波场函数 q(ω,k_x,k_z),可以得到TTI介质二阶耦合qP波波动方程^[10,18] :

$\begin{matrix} \{\begin{matrix} \partial^{2} \frac{p}{\partial} t^{2} = v_{px}^{2} H_{2} p + v_{pz} v_{pn} H_{1} q, \\ \partial^{2} \frac{q}{\partial} t^{2} = v_{pz} v_{pn} H_{2} p + v_{pz}^{2} H_{1} q 。 \end{matrix} (4) \end{matrix}$

偏微分算子为:

$\begin{matrix} \{\begin{matrix} H_{1} = si n^{2} θ \frac{\partial^{2}}{\partial x^{2}} + co s^{2} θ \frac{\partial^{2}}{\partial z^{2}} + \sin 2 θ \frac{\partial^{2}}{\partial x \partial z}, \\ H_{2} = co s^{2} θ \frac{\partial^{2}}{\partial x^{2}} + si n^{2} θ \frac{\partial^{2}}{\partial z^{2}} - \sin 2 θ \frac{\partial^{2}}{\partial x \partial z} 。 \end{matrix} (5) \end{matrix}$

在后面的公式中,为方便,均令v_pz=v_p,

那么,将式(4)整理得:

$\begin{matrix} \{\begin{matrix} \frac{\partial^{2} p}{\partial t^{2}} = v_{p}^{2} [(1 + 2 ε) (co s^{2} θ \frac{\partial^{2} p}{\partial x^{2}} + si n^{2} θ \frac{\partial^{2} p}{\partial z^{2}} - \sin 2 θ \frac{\partial^{2} p}{\partial x \partial z}) + \sqrt[]{1 + 2 δ} (si n^{2} θ \frac{\partial^{2} q}{\partial x^{2}} + co s^{2} θ \frac{\partial^{2} q}{\partial z^{2}} + \sin 2 θ \frac{\partial^{2} q}{\partial x \partial z})], \\ \frac{\partial^{2} q}{\partial t^{2}} = v_{p}^{2} [\sqrt[]{1 + 2 δ} (co s^{2} θ \frac{\partial^{2} p}{\partial x^{2}} + s i n^{2} θ \frac{\partial^{2} p}{\partial z^{2}} - \sin 2 θ \frac{\partial^{2} p}{\partial x \partial z}) + (si n^{2} θ \frac{\partial^{2} q}{\partial x^{2}} + co s^{2} θ \frac{\partial^{2} q}{\partial z^{2}} + \sin 2 θ \frac{\partial^{2} q}{\partial x \partial z})] 。 \end{matrix} (6) \end{matrix}$

2 TTI介质拟声波波动方程高阶有限差分

2.1 二阶偏导数高阶精度差分系数

对二阶偏导数项采用2N阶精度有限差分进行求解,以压制数值频散,提高数值模拟精度。差分系数的求取主要是通过泰勒级数展开法,通过求解以下系数方程

$\begin{matrix} (\begin{matrix} 1^{2} & 2^{2} & \dots & N^{2} \\ 1^{4} & 2^{4} & \dots & N^{4} \\ ︙ & ︙ & ⋱ & ︙ \\ 1^{2 N} & 2^{2 N} & \dots & N^{2 N} \end{matrix}) \times (\begin{matrix} c_{1} \\ c_{2} \\ ︙ \\ c_{N} \end{matrix}) = (\begin{matrix} 1 \\ 0 \\ ︙ \\ 0 \end{matrix}) (7) \end{matrix}$

可得差分系数^[19]。差分权系数为:

$\begin{matrix} \begin{matrix} c_{n} = \frac{{(- 1)}^{n + 1} \overset{N}{\prod_{i = 1, i \neq n}} i^{2}}{n^{2} \cdot \overset{n - 1}{\prod_{i = 1}} (n^{2} - i^{2}) \overset{N}{\prod_{i = n + 1}} (i^{2} - n^{2})} \\ (n = 1,2, \dots, N, N > 2) (8) \end{matrix} \end{matrix}$

其中: $\begin{matrix} c_{0} = - 2 \overset{N}{\sum_{i = 1}} c_{i} \end{matrix}$

则有:

$\begin{matrix} \begin{matrix} \frac{\partial^{2} u}{\partial x^{2}} = \frac{1}{{(Δx)}^{2}} {c_{0} u (x) + \\ \overset{N}{\sum_{n = 1}} c_{n} [u (x + nΔx) + u (x - nΔx)]} 。 (9) \end{matrix} \end{matrix}$

2.2 混合偏导数高阶精度差分系数

相比于各向同性介质和VTI介质波动方程,TTI介质波动方程中多了混合偏导数项,这在差分的过程中增加了计算复杂程度。空间混合偏导数可以先沿一个方向(如x)求取偏导数,再对其结果沿另一个方向(如z)求取偏导数得到。若函数u(x,z)的某阶混合偏导数连续,则该偏导数的结果与求导顺序无关。一阶偏导数的差分系数可以由方程

$\begin{matrix} (\begin{matrix} 1 & 2^{1} & \dots & N^{1} \\ 1 & 2^{3} & \dots & N^{3} \\ ︙ & ︙ & ⋱ & ︙ \\ 1 & 2^{2 N - 1} & \dots & N^{2 N - 1} \end{matrix}) (\begin{matrix} a_{1} \\ a_{2} \\ ︙ \\ a_{N} \end{matrix}) = (\begin{matrix} \frac{1}{2} \\ 0 \\ ︙ \\ 0 \end{matrix}) (10) \end{matrix}$

确定^[19]。解系数方程可得:

$\begin{matrix} a_{n} = \frac{{(- 1)}^{n + 1} \overset{N}{\prod_{i = 1, i \neq n}} i^{2}}{2 n \overset{n - 1}{\prod_{i = 1}} (n^{2} - i^{2}) \overset{N}{\prod_{i = n + 1}} (i^{2} - n^{2})}, (11) \end{matrix}$

该式所求即为一阶导数不同差分精度的差分权系数。则有:

$\begin{matrix} \frac{\partial u (x)}{\partial x} = \overset{N}{\sum_{n = 1}} \frac{a_{n} [u (x + nΔx) - u (x + nΔx)]}{Δx}, (12) \end{matrix}$

那么二阶混合偏导数可以写成:

$\begin{matrix} \begin{matrix} \frac{\partial^{2} u}{\partial x \partial z} = \frac{1}{ΔxΔz} \overset{N}{\sum_{\frac{n_{1} = - N}{n_{1} \neq 0}}} \overset{N}{\sum_{\frac{n_{2} = - N}{n_{2} \neq 0}}} \frac{n_{1}}{| n_{1} |} \cdot \\ \frac{n_{2}}{| n_{2} |} a_{| n_{1} |} a_{| n_{2} |} u (x + n_{1} Δx, z + n_{2} Δz), (13) \end{matrix} \end{matrix}$

这里, $\begin{matrix} a_{| n_{1} |} \end{matrix}$ 、 $\begin{matrix} a_{| n_{2} |} \end{matrix}$ 对应一阶导数的权系数值,且a₀=0。

根据以上推导可得二维 TTI 介质时间2阶、空间2N 阶差分精度的有限差分格式为:

$\begin{matrix} \{\begin{matrix} p_{i, j}^{k + 1} = 2 p_{i, j}^{k} - p_{i, j}^{k - 1} + v_{p}^{2} [(1 + 2 ε) (co s^{2} θ \cdot M_{1} + si n^{2} θ \cdot M_{2} - \sin 2 θ \cdot K_{1}) + \\ \sqrt[]{1 + 2 δ} (si n^{2} θ \cdot M_{3} + co s^{2} θ \cdot M_{4} + \sin 2 θ \cdot K_{2})], \\ q_{i, j}^{k + 1} = 2 q_{i, j}^{k} - q_{i, j}^{k - 1} + v_{p}^{2} [\sqrt[]{1 + 2 δ} (co s^{2} θ \cdot M_{1} + si n^{2} θ \cdot M_{2} - \sin 2 θ \cdot K_{1}) + \\ (si n^{2} θ \cdot M_{3} + co s^{2} θ \cdot M_{4} + \sin 2 θ \cdot K_{2})], \end{matrix} (14) \end{matrix}$

其中,

$\begin{matrix} \{\begin{matrix} M_{1} = \frac{1}{{(Δx)}^{2}} [c_{0} p_{i, j}^{k} + \overset{N}{\sum_{n = 1}} c_{n} (p_{i + n, j}^{k} + p_{i - n, j}^{k})] \\ M_{2} = \frac{1}{{(Δz)}^{2}} [c_{0} p_{i, j}^{k} + \overset{N}{\sum_{n = 1}} c_{n} (p_{i, j + n}^{k} + p_{i, j - n}^{k})] \\ M_{3} = \frac{1}{{(Δx)}^{2}} [c_{0} q_{i, j}^{k} + \overset{N}{\sum_{n = 1}} c_{n} (q_{i + n, j}^{k} + q_{i - n, j}^{k})] \\ M_{4} = \frac{1}{{(Δz)}^{2}} [c_{0} q_{i, j}^{k} + \overset{N}{\sum_{n = 1}} c_{n} (q_{i, j + n}^{k} + q_{i, j - n}^{k})] \\ K_{1} = \frac{1}{ΔxΔz} \overset{N}{\sum_{\frac{n_{1} = - N}{n_{1} \neq 0}}} \overset{N}{\sum_{\frac{n_{2} = - N}{n_{2} \neq 0}}} \frac{n_{1}}{| n_{1} |} \cdot \frac{n_{2}}{| n_{2} |} a_{| n_{1} |} a_{| n_{2} |} p_{i + n_{1}, j + n_{2}}^{k} \\ K_{2} = \frac{1}{ΔxΔz} \overset{N}{\sum_{\frac{n_{1} = - N}{n_{1} \neq 0}}} \overset{N}{\sum_{\frac{n_{2} = - N}{n_{2} \neq 0}}} \frac{n_{1}}{| n_{1} |} \cdot \frac{n_{2}}{| n_{2} |} a_{| n_{1} |} a_{| n_{2} |} q_{i + n_{1}, j + n_{2}}^{k} \end{matrix} \end{matrix}$

3 完全匹配层(PML)吸收边界条件

边界反射的处理是正演模拟中的一个重要问题。目前效果最好、应用最广泛的是1994年由Berenger^[20]从电磁领域引入的完全匹配层(PML)边界条件。理论上,该边界条件可以完全吸收以任意角度、任意频率入射的波动,很多学者在此基础上进行了研究^[21-22]。下面介绍了基于TTI介质的改进的PML边界控制方程,并给出了其相应的高阶有限差分格式。

首先,在复数空间进行坐标变换^[21-22],得到复数空间坐标系下的TTI介质二阶拟声波波动方程为(以方程(6)第一式的左边界为例,其它边界同理):

$\begin{matrix} \begin{matrix} {(iω)}^{2} P (x, z, ω) = (1 + 2 ε) \cdot v_{p}^{2} \cdot \{co s^{2} θ [{(\frac{iω}{iω + d_{x}})}^{2} \frac{\partial^{2} P (x, z, ω)}{\partial x^{2}} - \frac{{(iω)}^{2} d_{x}^{'}}{(iω + d_{x})^{3}} \frac{\partial P (x, z, ω)}{\partial x}] + \\ si n^{2} θ \frac{\partial^{2} P (x, z, ω)}{\partial z^{2}} - \sin 2 θ \frac{iω}{iω + d_{x}} \frac{\partial^{2} P (x, z, ω)}{\partial x \partial z}\} + \\ \sqrt[]{1 + 2 δ} \cdot v_{p}^{2} \cdot \{si n^{2} θ [{(\frac{iω}{iω + d_{x}})}^{2} \frac{\partial^{2} Q (x, z, ω)}{\partial x^{2}} - \frac{{(iω)}^{2} d_{x}^{'}}{(iω + d_{x})^{3}} \frac{\partial Q (x, z, ω)}{\partial x}] + \\ co s^{2} θ \frac{\partial^{2} Q (x, z, ω)}{\partial z^{2}} + \sin 2 θ \frac{iω}{iω + d_{x}} \frac{\partial^{2} Q (x, z, ω)}{\partial x \partial z}\}, (15) \end{matrix} \end{matrix}$

其中,P(x,z,ω),Q(x,z,ω)为p、q关于时间的傅里叶变换;d_x为衰减函数,x是到吸收层内边界的距离; $\begin{matrix} d_{x}^{'} \end{matrix}$ 为d_x'关于x的一阶导数。

对于TTI介质二阶拟声波方程(15),方程中等号右边共有6项二阶偏导数项,PML边界条件将波场分裂成6部分,即:

$\begin{matrix} P = P_{1} + P_{2} + P_{3} + P_{4} + P_{5} + P_{6} 。 (16) \end{matrix}$

结合方程(16),向方程(15)中引入中间变量G₁、G₂,得到PML吸收边界条件:

$\begin{matrix} \{\begin{matrix} (iω + d_{x})^{2} P_{1} (x, z, ω) = (1 + 2 ε) co s^{2} θ [v_{p}^{2} \frac{\partial^{2} P (x, z, ω)}{\partial x^{2}} + G_{1} (x, z, ω)], G_{1} (x, z, ω) = - v_{p}^{2} \frac{d_{x}^{'}}{iω + d_{x}} \frac{\partial P (x, z, ω)}{\partial x}, \\ {(iω)}^{2} P_{2} (x, z, ω) = v_{p}^{2} (1 + 2 ε) si n^{2} θ \frac{\partial^{2} P (x, z, ω)}{\partial z^{2}}, iω (iω + d_{x}) P_{3} (x, z, ω) = - v_{p}^{2} (1 + 2 ε) \sin 2 θ \frac{\partial^{2} P (x, z, ω)}{\partial x \partial z}, \\ (iω + d_{x})^{2} P_{4} (x, z, ω) = \sqrt[]{1 + 2 δ} si n^{2} θ [v_{p}^{2} \frac{\partial^{2} Q (x, z, ω)}{\partial x^{2}} + G_{2} (x, z, ω)], G_{2} (x, z, ω) = - v_{p}^{2} \frac{d_{x}^{'}}{iω + d_{x}} \frac{\partial Q (x, z, ω)}{\partial x}, \\ {(iω)}^{2} P_{5} (x, z, ω) = v_{p}^{2} \sqrt[]{1 + 2 δ} co s^{2} θ \frac{\partial^{2} Q (x, z, ω)}{\partial z^{2}}, iω (iω + d_{x}) P_{6} (x, z, ω) = v_{p}^{2} \sqrt[]{1 + 2 δ} \sin 2 θ \frac{\partial^{2} Q (x, z, ω)}{\partial x \partial z}, \\ P (x, z, ω) = P_{1} (x, z, ω) + P_{2} (x, z, ω) + P_{3} (x, z, ω) + P_{4} (x, z, ω) + P_{5} (x, z, ω) + P_{6} (x, z, ω), \end{matrix} (17) \end{matrix}$

对上式关于ω作傅里叶反变换,可以得到时间域的PML吸收边界条件:

$\begin{matrix} \{\begin{matrix} (\partial_{t} + d_{x})^{2} p_{1} = (1 + 2 ε) co s^{2} θ [v_{p}^{2} \cdot \frac{\partial^{2} p}{\partial x^{2}} + g_{1} (x, z, t)] \\ (\partial_{t} + d_{x}) g_{1} (x, z, t) = - v_{p}^{2} \cdot d_{x}^{'} \frac{\partial p}{\partial x} \\ \partial_{t}^{2} p_{2} = v_{p}^{2} (1 + 2 ε) si n^{2} θ \frac{\partial^{2} p}{\partial z^{2}} \\ (\partial_{t}^{2} + d_{x} \cdot \partial_{t}) p_{3} = - v_{p}^{2} (1 + 2 ε) \sin 2 θ \frac{\partial^{2} p}{\partial x \partial z} \\ (\partial_{t} + d_{x})^{2} p_{4} = \sqrt[]{1 + 2 δ} si n^{2} θ [v_{p}^{2} \cdot \frac{\partial^{2} q}{\partial x^{2}} + g_{2} (x, z, t)] \\ (\partial_{t} + d_{x}) g_{2} (x, z, t) = - v_{p}^{2} \cdot d_{x}^{'} \frac{\partial q}{\partial x} \\ \partial_{t}^{2} p_{5} = v_{p}^{2} \cdot \sqrt[]{1 + 2 δ} co s^{2} θ \frac{\partial^{2} q}{\partial z^{2}} \\ (\partial_{t}^{2} + d_{x} \cdot \partial_{t}) p_{6} = v_{p}^{2} \cdot \sqrt[]{1 + 2 δ} \sin 2 θ \frac{\partial^{2} q}{\partial x \partial z} \\ p = p_{1} + p_{2} + p_{3} + p_{4} + p_{5} + p_{6} \end{matrix} (18) \end{matrix}$

其中,g₁(x,z,t)、g₂(x,z,t)分别为G₁(x,z,ω)、G₂(x,z,ω)的傅里叶反变换。

那么,很容易得出方程(6)第一式左边界的PML边界高阶有限差分格式为:

$\begin{matrix} \{\begin{matrix} p_{1 i, j}^{k + 1} = \frac{1}{1 + Δt d_{x} + \frac{(Δt d_{x})^{2}}{2}} \{2 p_{1 i, j}^{k} + [Δt d_{x} - \frac{(Δt d_{x})^{2}}{2} - 1] p_{1 i, j}^{k - 1} + Δ t^{2} (1 + 2 ε) co s^{2} θ (v_{p}^{2} \cdot M_{1} + \frac{g_{1 i, j}^{k + 1} + g_{1 i, j}^{k}}{2})\} \\ g_{1 i, j}^{k + 1} = \frac{1}{1 + \frac{Δt \cdot d_{x}}{2}} [(1 - \frac{Δt \cdot d_{x}}{2}) g_{1 i, j}^{k} - v_{p}^{2} \cdot Δt \cdot d_{x}^{'} \cdot D_{1}] \\ p_{2 i, j}^{k + 1} = 2 p_{2 i, j}^{k} - p_{2 i, j}^{k - 1} + Δ t^{2} \cdot v_{p}^{2} (1 + 2 ε) si n^{2} θ \cdot M_{2} \\ p_{3 i, j}^{k + 1} = \frac{1}{1 + \frac{d_{x} \cdot Δt}{2}} [2 p_{3 i, j}^{k} + (\frac{d_{x} \cdot Δt}{2} - 1) \cdot p_{3 i, j}^{k - 1} - v_{p}^{2} Δ t^{2} (1 + 2 ε) \sin 2 θ \cdot K_{1}] \\ p_{4 i, j}^{k + 1} = \frac{1}{1 + Δt d_{x} + \frac{(Δt d_{x})^{2}}{2}} \{2 p_{4 i, j}^{k} + [Δt d_{x} - \frac{(Δt d_{x})^{2}}{2} - 1] p_{4 i, j}^{k - 1} + Δ t^{2} \sqrt[]{1 + 2 δ} si n^{2} θ (v_{p}^{2} \cdot M_{3} + \frac{g_{2 i, j}^{k + 1} + g_{2 i, j}^{k}}{2})\} \\ g_{2 i, j}^{k + 1} = \frac{1}{1 + \frac{Δt \cdot d_{x}}{2}} [(1 - \frac{Δt \cdot d_{x}}{2}) g_{2 i, j}^{k} - v_{p}^{2} \cdot Δt \cdot d_{x}^{'} \cdot D_{2}] \\ p_{5 i, j}^{k + 1} = 2 p_{5 i, j}^{k} - p_{5 i, j}^{k - 1} + Δ t^{2} \cdot v_{p}^{2} \sqrt[]{1 + 2 δ} co s^{2} θ \cdot M_{4} \\ p_{6 i, j}^{k + 1} = \frac{1}{1 + \frac{d_{x} \cdot Δt}{2}} [2 p_{6 i, j}^{k} + (\frac{d_{x} \cdot Δt}{2} - 1) \cdot p_{6 i, j}^{k - 1} + v_{p}^{2} Δ t^{2} \sqrt[]{1 + 2 δ} \sin 2 θ \cdot K_{2}] \\ p_{i, j}^{k + 1} = p_{1 i, j}^{k + 1} + p_{2 i, j}^{k + 1} + p_{3 i, j}^{k + 1} + p_{4 i, j}^{k + 1} + p_{5 i, j}^{k + 1} + p_{6 i, j}^{k + 1} \end{matrix} (19) \end{matrix}$

其中:

$\begin{matrix} \{\begin{matrix} M_{1} = \frac{1}{{(Δx)}^{2}} [c_{0} p_{i, j}^{k} + \overset{N}{\sum_{n = 1}} c_{n} (p_{i + n, j}^{k} + p_{i - n, j}^{k})] \\ M_{2} = \frac{1}{{(Δz)}^{2}} [c_{0} p_{i, j}^{k} + \overset{N}{\sum_{n = 1}} c_{n} (p_{i, j + n}^{k} + p_{i, j - n}^{k})] \\ M_{3} = \frac{1}{{(Δx)}^{2}} [c_{0} q_{i, j}^{k} + \overset{N}{\sum_{n = 1}} c_{n} (q_{i + n, j}^{k} + q_{i - n, j}^{k})] \\ M_{4} = \frac{1}{{(Δz)}^{2}} [c_{0} q_{i, j}^{k} + \overset{N}{\sum_{n = 1}} c_{n} (q_{i, j + n}^{k} + q_{i, j - n}^{k})] \\ D_{1} = \frac{1}{Δx} \overset{N}{\sum_{n = 1}} a_{n} [p_{i + n, j}^{k} - p_{i - n, j}^{k}] \\ D_{2} = \frac{1}{Δx} \overset{N}{\sum_{n = 1}} a_{n} [q_{i + n, j}^{k} - q_{i - n, j}^{k}] \\ K_{1} = \frac{1}{ΔxΔz} \overset{N}{\sum_{\frac{n_{1} = - N}{n_{1} \neq 0}}} \overset{N}{\sum_{\frac{n_{2} = - N}{n_{2} \neq 0}}} \frac{n_{1}}{| n_{1} |} \cdot \frac{n_{2}}{| n_{2} |} a_{| n_{1} |} a_{| n_{2} |} p_{i + n_{1}, j + n_{2}}^{k} \\ K_{2} = \frac{1}{ΔxΔz} \overset{N}{\sum_{\frac{n_{1} = - N}{n_{1} \neq 0}}} \overset{N}{\sum_{\frac{n_{2} = - N}{n_{2} \neq 0}}} \frac{n_{1}}{| n_{1} |} \cdot \frac{n_{2}}{| n_{2} |} a_{| n_{1} |} a_{| n_{2} |} q_{i + n_{1}, j + n_{2}}^{k} \end{matrix} (20) \end{matrix}$

同理可得方程(6)第二式左边界的PML边界控制方程及高阶有限差分格式。

在利用完全匹配层边界条件对边界进行处理时,衰减函数的选择至关重要。如果选择的衰减函数不合适,将会产生不同程度的边界反射,无法发挥出PML边界的优势。Hastings等^[23]、Collino等^[21]、Groby等^[24]都对衰减函数进行了发展,以下为文中所用的衰减函数:

$\begin{matrix} \{\begin{matrix} d_{i} (i) = d_{0} {(i / δ)}^{4} \\ d_{0} = \frac{3 v_{p}}{2 δ} \cdot \log (\frac{1}{R}) \end{matrix} (21) \end{matrix}$

其中,i为计算网格点到吸收层内边界的距离,分为x、z两个方向;δ为PML吸收层的厚度;R为理论反射系数,一般取0.01~1.0×10^-⁶;v_p为波速;d₀为反射系数函数。

根据PML边界条件的理论,分别在波场的x、z方向对波场边界进行衰减,那么,在计算区域外就会被划分为8个区域,如图1所示,其中箭头所指的区域为内部计算区域,剩下的区域为边界。文中采用与以往PML边界不同的划分方案,分别沿x、z方向把波场分成3个区域,两个衰减区域,一个内部计算区域。图2所示,为文中采用的衰减函数d_x(x)、d_z(z)的设置方式。可以看出,将图2所示的两个衰减函数衰减示意图重叠之后的未被衰减的区域与图1所示的内部计算区域是一致的。这样划分可以保证在4个角点区域波场沿x方向和z方向都被衰减,不需要单独对边界4个角点区域进行处理,同时又避免了以往PML边界8个衰减区域的繁琐,编程实现时更简便。PML边界处理示意如图2。

[1]	牟永光,裴正林.三维复杂介质地震数值模拟[M].北京:石油工业出版社,2005:1-13. [本文引用: 1]
[2]	Alkhalifah T. Acoustic approximations for processing in transversely isotropic media [J].Geophysics,1998,63(2):623-631. [本文引用: 1]
[3]	Alkhalifah T. An acoustic wave equation for anisotropic media [J].Geophysics,2000,65(4):1239-1250. [本文引用: 1]
[4]	Zhou H,Zhang G,Bloor R. An anisotropic acoustic wave equation for VTI media [C]//68^th EAGE Conference and Exhibition Incorporating SPE EUROPEC 2006.Vienna,Austria:EAGE,2006:H033. [本文引用: 1]
[5]	Du X,Fletcher R,Fowler P J. A new pseudo-acoustic wave equation for VTI media [C]//70^th EAGE Conference and Exhibition Incorporating SPE EUROPEC 2008.Rome:European Association of Geoscientists and Engineers,2008:1-5. [本文引用: 1]
[6]	Hestholm S. Acoustic VTI modeling using high-order finite differences [J].Geophysics,2009,74(5):T67-T73. [本文引用: 1]
[7]	Zhou H B,Zhang G Q,Bloor R. An anisotropic acoustic wave equation for modeling and migration in 2D TTI media [C]//Expanded Abstracts of 76^th SEG Annual International Meeting.Tulsa,OK:SEG,2006:194-198. [本文引用: 1]
[8]	Fletcher R,Du X,Fowler P J. A new pseudo-acoustic wave equation for TI media [C]//Expanded Abstracts of 78^th SEG Annual Internet Meeting.Tulsa,OK:SEG,2008:2082-3086. [本文引用: 2]
[9]	Duveneck E,Milcik P,Bakker P M,et al. Acoustic VTI wave equations and their application for anisotropic reverse-time migration [C]//Expanded Abstracts of 78^th SEG Annual International Meeting.Tulsa,OK:SEG,2008:2186-2190. [本文引用: 2]
[10]	Duveneck E,Bakker P M. Stable P-wave modeling for reverse-time migration in tilted TI media [J].Geophysics,2011,76(2):S65-S75. [本文引用: 2]
[11]	Fowler P J,Du X,Fletcher R P. Coupled equations for reverse time migration in transversely isotropic media [J].Geophysics,2010,75(1):S11-S22. [本文引用: 1]
[12]	Liu F Q,Morton S A,Jiang S S,et al. Decoupled wave equations for P and SV waves in an acoustic VTI media [C]//SEG Technical Program Expanded Abstracts 2009.Tulsa,OK:SEG,2009:2844-2848. [本文引用: 2]
[13]	Pestana R C,Ursin B,Stoffa P L. Separate P- and SV-wave equations for VTI media [C]//SEG Technical Program Expanded Abstracts 2011.Tulsa,OK:SEG,2011:163-167. [本文引用: 1]
[14]	Chu C L,Macy B K,Anno P D. An accurate and stable wave equation for pure acoustic TTI modeling [C]//SEG Technical Program Expanded Abstracts 2011.Tulsa,OK:SEG,2011:179-184. [本文引用: 3]
[15]	Zhan G,Pestana R C,Stoffa P L. An acoustic wave equation for pure P wave in 2D TTI media [C]//SEG Technical Program Expanded Abstracts 2011.Tulsa,OK:SEG,2011:168-173. [本文引用: 3]
[16]	Barrera Pacheco D F,Pestana R,Vivas F. New pseudo-acoustic wave equations for modeling and reverse time migration in TTI media [C]//75^th EAGE Conference & Exhibition Incorporating SPE EUROPEC 2013.London:SPE,2013. [本文引用: 1]
[17]	Thomsen L. Weak elastic anisotropy [J].Geophysics,1986,51(10):1954-1966. [本文引用: 1]
[18]	张岩,吴国忱. TTI介质qP波逆时偏移中伪横波噪声压制方法 [J].地球物理学报,2013,56(6):2065-2076. [本文引用: 2]
[19]	刘洋,李承楚,牟永光. 任意偶数阶精度有限差分法数值模拟 [J].石油地球物理勘探,1998,33(1):1-10. [本文引用: 2]
[20]	Berenger J P. A perfectly matched layer for the absorption of electromagnetic waves [J].Journal of Computational Physics,1994,114(2):185-200. [本文引用: 1]
[21]	Collino F,Tsogka C. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media [J].Geophysics,2001,66(1):294-307. [本文引用: 3]
[22]	Komatitsch D,Tromp J. A Perfectly Matched Layer absorbing boundary condition for the second-order seismic wave equation [J].Geophysical Journal International,2003,154(1):146-153. [本文引用: 2]
[23]	Hastings F D,Schneider J B,Broschat S L. Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation [J].The Journal of the Acoustical Society of America,1996,100(5):3061-3069. [本文引用: 1]
[24]	Groby J P,Tsogka C. A time domain method for modeling viscoacoustic wave propagation [J].Journal of Computational Acoustics,2004,14(2):201-236. [本文引用: 1]
[25]	Tsvankin I. Seismic signatures and analysis of reflection data in anisotropic media [M].London:Pergamon,2001. [本文引用: 1]
[26]	Grechka V,Zhang L B,Rector J W. Shear waves in acoustic anisotropic media [J].Geophysics,2004,69(2):576-582. [本文引用: 3]
[27]	Xu S,Zhou H B. Accurate simulations of pure quasi-P-waves in complex anisotropic media [J].Geophysics,2014,79(6):T341-T348. [本文引用: 1]

TTI介质拟声波方程数值模拟

Numerical simulation of pseudo acoustic wave equation in TTI medium

0 引言

1 稳定的TTI介质拟声波方程

2 TTI介质拟声波波动方程高阶有限差分

2.1 二阶偏导数高阶精度差分系数

2.2 混合偏导数高阶精度差分系数

3 完全匹配层(PML)吸收边界条件

4 稳定性分析

5 伪横波的压制方法

6 数值试验

6.1 均匀TTI介质模型

6.2 TTI介质BP_partI模型

6.3 TTI介质BP_partII模型

7 结论

参考文献 文献选项 原文顺序 文献年度倒序 文中引用次数倒序 被引期刊影响因子

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子