基于卷积完全匹配层的旋转交错网格高阶差分法模拟弹性波传播

doi:10.11720/wtyht.2018.1320

基于卷积完全匹配层的旋转交错网格高阶差分法模拟弹性波传播

冯德山¹^,², 王向宇¹^,²

1. 中南大学地球科学与信息物理学院,湖南长沙 410083

2. 有色资源与地质灾害探查湖南省重点实验室,湖南长沙 410083

Elastic wave propagation simulation in anisotropic media and random media using high-order difference method of rotation staggered grids based on convolutional perfectly matched layer

FENG De-Shan¹^,², WANG Xiang-Yu¹^,²

1. School of Geosciences and Info-Physics,Central South University,Changsha 410083,China

2. Key Laboratory of Non-ferrous Resources and Geological Detection of Hunan Province,Changsha 410083,China

责任编辑: 叶佩

收稿日期: 2017-07-17 修回日期: 2018-03-22 网络出版日期: 2018-08-05

基金资助:

国家自然科学基金资助项目.  41574116
中南大学创新驱动项目.  2015CX008
中南大学教师研究基金资助项目.  2014JSJJ001
中南大学升华育英人才计划.  2012
湖湘青年创新创业平台培养对象共同资助项目.  2013

Received: 2017-07-17 Revised: 2018-03-22 Online: 2018-08-05

Fund supported:

. 41574116

. 2015CX008

. 2014JSJJ001

. 2012

. 2013

作者简介 About authors

冯德山(1978-),博士,教授、博导,主要从事地球物理正反演与数据处理研究工作。Email:fengdeshan@126.com 。

摘要

从一阶速度—应力弹性波方程出发,基于旋转交错网格,推导了时间二阶精度空间2M阶精度的有限差分离散格式。阐述了递归卷积复频移完全匹配层(CPML)边界条件的原理,建立了一阶速度—应力弹性波高阶差分CPML边界条件的递推公式。开展了CPML边界中关键参数m、κ和α的选取实验,通过分析反射误差分布图,选取了CPML边界条件中最优参数。全局反射误差与波场快照都说明,CPML较PML对隐失波具有更优的吸收性能。基于Matlab平台,编写了基于CPML边界的旋转交错网格弹性波正演模拟程序,应用该程序对各向异性介质及随机介质进行了模拟,得到了弹性波正演剖面记录及波场快照,通过对正演剖面记录及波场快照的分析,可以更清楚地了解弹性波在各向异性介质及随机介质的传播特性,指导非均匀介质中地震勘探资料解释。

关键词： 旋转交错网格 ; 高阶差分法 ; 卷积完全匹配层 ; 各向异性介质 ; 随机介质 ; 正演

Abstract

On the basis of the first-order velocity-stress elastic wave equation,a finite-difference discretization scheme with 2M-order accuracy of second-order space accuracy is deduced based on rotational staggered grid.The principle of recursive convolutional complex frequency shift perfectly matched layer (CPML) boundary condition is expounded,and the recurrence formula of first order velocity-stress elastic wave high order difference CPML boundary condition is established.The selection of the key parameter m,and sum of CPML boundary is carried out.By analyzing the distribution of reflection error,the optimal parameters of CPML boundary condition are selected.Both global reflection errors and wavefield snapshots indicate that CPML has better absorption properties than evanescent waves.On the basis of the Matlab platform,a program of elastic staggered forward simulation of rotating staggered grid based on CPML boundary is developed.The program is used to simulate the anisotropic medium and random medium.The seismic wave forward record and wavefield snapshot are obtained.The forward propagation profile and the wavefield snapshots can be used to better understand the propagation characteristics of seismic waves in anisotropic media and random media and to guide the interpretation of seismic exploration data in heterogeneous media.

Keywords： rotation staggered grids ; high-order difference method ; convolutional perfectly matched layer ; anisotropic media ; random media ; forward

PDF (7848KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

冯德山, 王向宇. 基于卷积完全匹配层的旋转交错网格高阶差分法模拟弹性波传播. 物探与化探[J], 2018, 42(4): 766-776 doi:10.11720/wtyht.2018.1320

FENG De-Shan, WANG Xiang-Yu. Elastic wave propagation simulation in anisotropic media and random media using high-order difference method of rotation staggered grids based on convolutional perfectly matched layer. Geophysical and Geochemical Exploration[J], 2018, 42(4): 766-776 doi:10.11720/wtyht.2018.1320

0 引言

弹性波正演是人们认识地下复杂介质的重要手段,在弹性波有限差分正演模拟中,除了二阶微分方程外,还常常采用一阶速度—应力方程,它的主要优点是不需对弹性常数进行空间微分^[1,2],正演计算能得到完整的波场响应。除了规则网格有限差分法外,Madariaga^[3]提出一种较为先进的交错网格,Virieux^[4,5]将其应用到各向同性介质中的SH波和P-SV波模拟中,它在不增加计算工作量和存储空间的前提下,和常规网格相比局部精度提高4倍,收敛速度也很快。裴正林^[6,7]、孙卫涛和杨慧珠^[8]将交错网格高阶有限差分法应用到双相各向异性介质弹性波传播中;陈祖银等^[9]、Chen^[10]将高阶交错网格有限差分法应用到弹性波场分离的数值模拟中,以便能更精确地研究纵波、横波在介质的传播规律;奚先和姚姚^[11,12,13]将交错网格有限差分法应用于二维弹性随机介质的正演及波场特征研究;但是,与标准交错网格有限差分方法相比,旋转交错网格有限差分的好处在于不同的物理量只位于两个不同的位置:应力和应变(或质点速度和位移)位于离散单元的中心,而质点速度或位移(或应力和应变)位于单元的顶点。经过这样的处理后,弹性常数就不再需要进行平均(模拟非均匀介质时)或内插(模拟各向异性介质时),因为此时所有的弹性模量都位于相同的位置,且与应力或应变的位置相对应^[14]。

应用有限差分法开展弹性波模型时,需要在截断人工边界处引入吸收边界条件。目前吸收效果较好的边界条件为完全匹配层(PML)边界^[15,16,17],它可以使弹性波在匹配层中无反射地通过并按指数规律衰减,适合于吸收较大范围的入射角。Hastings等^[18]将PML方法应用于弹性波的有限差分正演中;Drosseart^[19,20]提出了一种基于递归积分复频移PML算法;Zhang和Shen^[21]、Martin等^[22] 及张显文等^[23]将基于辅助微分方程的复频移PML应用到弹性波数值模拟中,该方法易于扩展到高阶差分形式;Roden和Gedney^[24]应用卷积技术实现了完全匹配层加载,大大提高了运算效率;张鲁新等^[25]将不分裂的卷积完全匹配层(CPML)结合旋转交错网格有限差分技术应用到孔隙弹性介质模拟中。使用常规的PML算法处理截断边界时仍会留下一些反射波和大量的隐失波,而采用CPML算法处理截断边界则可以减少反射波的存在,并且对隐失波的吸收具有很好的效果^[23]。

为了提高差分精度,减小网格频散,笔者将速度(应力)对时间的奇数阶高阶导数转化为应力(速度)对空间的导数,运用时间二阶精度、空间十六阶精度高阶差分法,由于随机介质中存在大量的波阻抗差较大的界面,并且需要尽可能减小界面反射,消除隐失波对波剖面的影响,因此采用标准交错网格结合常规PML显然不适合模拟随机介质,更为主要的是在截断界面处的波反射会干扰波剖面的分析^[13],因此需要引入旋转交错网格结合递归卷积复频移PML技术,开展各向异性介质与随机介质模型的弹性波正演。

1 旋转交错网格有限差分法

一阶速度—应力弹性波方程不需对弹性常数进行空间微分,且方程中仅含有速度和应力变量关于时间和空间的一阶偏导数,容易提高差分近似的阶数。为此,将二维各向异性介质及随机介质的正演控制方程表示为^[1]:

(1)

\{\begin{array}{l} ρ \frac{\partial V_{x}}{\partial t} = \frac{\partial τ_{xx}}{\partial x} + \frac{\partial τ_{xz}}{\partial z}, \\ ρ \frac{\partial V_{z}}{\partial t} = \frac{\partial τ_{zx}}{\partial x} + \frac{\partial τ_{zz}}{\partial z}, \\ \frac{\partial τ_{xx}}{\partial t} = C_{11} \frac{\partial V_{x}}{\partial x} + C_{13} \frac{\partial V_{z}}{\partial z} + C_{15} (\frac{\partial V_{z}}{\partial x} + \frac{\partial V_{x}}{\partial z}), \\ \frac{\partial τ_{zz}}{\partial t} = C_{31} \frac{\partial V_{x}}{\partial x} + C_{33} \frac{\partial V_{z}}{\partial z} + C_{35} (\frac{\partial V_{z}}{\partial x} + \frac{\partial V_{x}}{\partial z}), \\ \frac{\partial τ_{zx}}{\partial t} = C_{51} \frac{\partial V_{x}}{\partial x} + C_{53} \frac{\partial V_{z}}{\partial z} + C_{55} (\frac{\partial V_{z}}{\partial x} + \frac{\partial V_{x}}{\partial z}) 。 \end{array}

式(1)中V_x与V_z为别为x方向和z方向的质点速度分量,τ_xx与τ_zz为正应力,τ_xz与τ_zx为切应力,ρ为密度,C_ij(i,j=1,3,5)为介质的弹性参数。目前普遍采用图1a所示的标准交错网格有限差分^[3]对式(1)进行数值模拟。但是,当遇到波阻抗较大的界面时,在自由界面处,边界条件必须在差分算法中显式表示出来,标准交错网格有限差分变得不再稳定。Saenger等人^[26]为了改善标准交错网格的性能,提出了图1b所示的旋转交错网格,图中实线和虚线分别表示空间上的整网格和半网格。旋转交错网格在标准网格交错网格的基础上增加了网格旋转,将同一物理量的不同分量定义在同一网格点上,其中速度分量和密度分量定义在整网格点,应力分量和弹性常量定义在半网格点上。因此,在模拟非均匀性比较强或各向异性介质时,不需要对介质参数进行插值处理,在遇到自由界面时也不需要进行单独的显式处理。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 普通交错网格与旋转交错网格单元示意

a—标准交错网格;b—旋转交错网格

在旋转交错网格计算过程中,需要变换求导方向,先计算沿对角线的物理量的差分,然后利用对角线差分值的线性组合计算沿坐标轴的差分。采用两套坐标系,旧坐标系(x,z)的方向与整网格平行,新坐标系( $\overset{̅}{x}$ , $\overset{̅}{z}$ )与半网格平行,两坐标系的变换关系如下:

(2)

\{\begin{array}{l} \overset{̅}{z} = (Δ x / Δ r) \cdot x + (Δ z / Δ r) \cdot y, \\ \overset{̅}{x} = (Δ x / Δ r) \cdot x - (Δ z / Δ r) \cdot y 。 \end{array}

式(2)中Δx、Δz为网格间距,Δr= $\sqrt[]{Δ x^{2} + Δ z^{2}}$ ,两坐标系下的求导算子有如下关系:

(3)

\{\begin{array}{l} \partial / \partial z = (\partial / \partial \overset{̅}{z} - \partial / \partial \overset{̅}{x}) \cdot Δ r / 2 Δ z, \\ \partial / \partial x = (\partial / \partial \overset{̅}{z} + \partial / \partial \overset{̅}{x}) \cdot Δ r / 2 Δ x 。 \end{array}

式中 $\overset{̅}{x}$ 和 $\overset{̅}{z}$ 代表对角线方向,Δx和Δz是沿坐标轴方向的空间步长,Δr是对角线的长度,其中Δr= $\sqrt[]{Δ x^{2} + Δ z^{2}}$ 。∂/∂ $\overset{̅}{x}$ 和∂/∂ $\overset{̅}{z}$ 的离散差分算子的计算类似于标准交错网格中∂/∂x和∂/∂z的离散差分算子。应用旋转交错网格FDTD对式(1)进行离散,取x=iΔx,z=jΔz,t=kΔt,i、j、k分别表示空间和时间网格点,U、W分别代表速度分量V_x、V_z的离散值,P、Q、S分别代表应力τ_xx、τ_zz、τ_xz的离散值,得到时间二阶精度空间2M阶(本文取16阶)精度有限差分的离散格式见附录A。

根据中的旋转交错网格高阶差分格式,可以开展各向异性、随机介质中弹性波波动方程的正演。

下面分别给出标准交错网格和旋转交错网格的纽曼稳定性标准条件,对于二维一阶速度—应力弹性波方程时间2M阶精度、空间2N阶精度的标准交错网格有限差分方程的稳定性条件为:

(4)

\frac{Δ t v_{p}}{Δ h} \leq \frac{1}{\sqrt[]{2} \overset{N}{\sum_{n = 1}} | K_{n}^{(N)} |} 。

对于二维一阶速度—应力弹性波方程时间2M阶精度、空间2N阶精度的旋转交错网格有限差分方程的稳定性条件为:

(5)

\frac{Δ t v_{p}}{Δ h} \leq \frac{1}{\overset{N}{\sum_{n = 1}} | K_{n}^{(N)} |},

式中: $K_{n}^{(N)}$ 为差分系数,v_p为最大纵波速度,Δh为空间步长,Δt为时间步长,通过对比标准交错网格和旋转交错网格的稳定性条件可以发现,标准交错网格的稳定性条件较旋转交错网格的稳定性条件更为严格,即标准交错网格的数值稳定性较差。

2 递归卷积复频移完全匹配层边界条件

根据Kuzuoglu和Mittra^[27]的研究,定义对应坐标的微分算子和复频移拉伸函数分别为:

(6)

\partial_{\overset{̅}{p}} = \partial_{p} / S_{p},

(7)

S_{p} = κ_{p} + σ_{p} / (α_{p} + i ω), p = x, y, z

式中:κ_p为辅助衰减系数,它的引入是为了改善PML对表面波的吸收特性,σ_p为PML层内p方向衰减因子,α_p值则是为了改善PML对低频分量的吸收特性,σ_p和α_p为大于零的正实数,κ_p≥1的正实数。将式(7)代入式(6),并做傅立叶逆变换,得

(8)

\partial_{\overset{̅}{p}} = {\overset{̅}{S}}_{p} * \partial_{p},

式中, ${\overset{̅}{S}}_{p}$ 为1/S_p的傅立叶逆变换。

(9)

{\overset{̅}{S}}_{p} (t) = δ (t) / κ_{p} - [σ_{p} (t) / κ_{p}^{2}] \cdot H (t) e^{- [σ_{p} / (κ_{p} + α_{p})] t}

式(9)中δ(t)为冲激函数,H(t)为单位阶跃函数。令

(10)

ζ_{p} (t) = - [σ_{p} (t) / κ_{p}^{2}] \cdot H (t) e^{- [σ_{p} / (κ_{p} + α_{p})] t},

则

(11)

\partial_{\overset{̅}{p}} = δ (t) / κ_{p} + ζ_{p} * \partial_{p} 。

式(10)可分为两部分,第一部分只要计算空间偏导数除以系数κ_p就可以了,容易计算。而第二部分的计算需要进行卷积。在离散情况下,将nΔt时刻的卷积写为:

(12)

\begin{array}{r} ψ_{p}^{n Δ t} & = (ζ_{p} * \partial_{p})^{n Δ t} \\ = \int_{0}^{n} \partial_{p})^{n Δ t - τ} ζ_{p} (τ) d τ \\ = \overset{n - 1}{\sum_{m = 0}} \int_{m}^{Δ} \partial_{p})^{n Δ t - τ} ζ_{p} (τ) d τ, \end{array}

由于采用交错网格,则

(13)

\begin{array}{r} ψ_{p}^{n Δ t} & = \overset{n - 1}{\sum_{m = 0}} (\partial_{p})^{[n - (m + 1) / 2] Δ t} \int_{m}^{Δ} ζ_{p} (τ) d τ \\ = \overset{n - 1}{\sum_{m = 0}} Z_{p} (m) (\partial_{p})^{[n - (m + 1) / 2] Δ t}, \end{array}

其中:

(14)

\begin{array}{r} Z_{p} (m) & = \int_{m}^{Δ} ζ_{p} (τ) d τ \\ = - \frac{σ_{p} (t)}{κ_{p}^{2}} \int_{m}^{Δ} e^{- (σ_{p} / κ_{p} + α_{p}) t} d τ \\ = a_{p} e^{- (σ_{p} / κ_{p} + α_{p}) m Δ t}, \end{array}

式(13)中,

(15)

\begin{array}{l} a_{p} = \frac{σ_{p}}{κ_{p} (σ_{p} + κ_{p} α_{p})} (b_{p} - 1), \\ b_{p} = e^{- (σ_{p} / κ_{p} + α_{p}) Δ t} 。 \end{array}

式(14)为简单的指数形式,因此可写为递推形式:

$ψ_{p}^{n}$ =b_p $ψ_{p}^{n - 1}$ +a_p(α_p)ⁿ^-1^/²。

对于一阶速度—应力弹性波方程,在CPML吸收边界下的波动方程见附录B。

根据CPML边载方式可将整个模拟区域划分为内部计算区域和CPML区域。为了使加载的CPML吸收效果最好,需要进一步探讨吸收边界参数的选取。根据文献[26]对σ_p的设置方法为:

(16)

σ_{p} = \frac{(m + 1) V_{p}}{2 d_{PML}} \ln (\frac{1}{R}) {(\frac{x}{d_{PML}})}^{m},

式中:d_PML为吸收层厚度,x是CPML内计算点到CPML内边界的距离,R为理论反射系数,常取R=10^-6,m为多项式的阶数。对κ_p与α_p的设置方法为:

(17)

\{\begin{array}{l} κ_{p} = 1 + (κ_{\max} - 1) (x / d_{PML})^{m}, \\ α_{p} = π α_{\max} (1 - x / d_{PML}) 。 \end{array}

为定量比较CPML与PML在不同入射角情况下的吸收效果,采用均匀介质模型,比较不同入射角处检波器记录到的波形。图2所示模型C以B₁为左边界,大小为1 300 m×2 300 m。将模型C的左端延伸构造参考模型D,模型D以B₂为左边界,大小为2 300 m×2 300 m。设置纵波速度v_p=3 200 m/s,横波速度v_s=2 600 m/s,密度ρ=2 000 kg/m³。模拟使用的网格单元大小为Δx=Δz=5 m,时间步长为0.5 ms,记录长度为400 ms,震源位于模型D的中心,子波主频为30 Hz的零相位雷克子波。接收排列与吸收边界B₁平行,间距为5 m,距离吸收边界B₁右端50 m,最上面的接收点距离上边界为450 m,最下面的接收点与下边界距离450 m,接收排列共280个接收点。分别对模型C和模型D进行正演,震源和接收点位置相同,对于模型C,震源激发的波以不同入射角穿过接收排列,进入吸收边界B₁,如果仍不能被完全吸收,将有部分波反射到C被接收排列记录;而对模型D,由于吸收边界B₁位置没有吸收边界,因此弹性波在此处不会发生反射。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 观测系统和模型设置

取CPML参数R=10^-6,κ_max=1,α_max=0,即为常规的PML边界条件时,得到参数m和反射系数R的反射误差分布图3,通过分析接收排列关于参数m反射误差分布图可知,综合各个角度入射的反射误差分布结果,当m=3时反射误差达到最小。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 接收排列与参数m关系的反射误差分布

当m=3,R=10^-6,α_max=0时,讨论参数κ_max对反射误差的影响,模型与接收排列设置同上,得到接收排列关于参数κ_max的反射误差分布图4。通过分析接收排列关于参数κ_max反射误差分布图可知,综合各个角度入射的反射误差分布结果,当κ_max=3时反射误差达到最小。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 接收排列与参数κ_max关系的反射误差分布

当m=3,R=10^-6,κ_max=3时,讨论参数α_max对反射误差的影响,得到接收排列关于参数α_max的反射误差分布如图5所示。通过分析接收排列关于参数α_max反射误差分布图,综合各个角度入射的反射误差分布结果可知,当α_max=90时反射误差达到最小即α_max=3f₀,其中f₀为雷克子波的主频。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 接收排列与参数α_max关系的反射误差分布

为了比较CPML与PML对于大角度入射的掠角波的吸收效果,取第280个接收点的检波记录,PML吸收边界与CPML吸收边界反射误差对比如图6。图中可见,PML的最大误差在-57.9377 dB的数量级,而CPML的最大误差则得到了很好地改善,在-68.1107 dB的数量级这比传统的PML提高了大约10 dB,由此证明CPML吸收边界条件的反射误差较PML要更小,说明CPML效果更优。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 CPML与PML反射误差对比

设置窄长模型,进一步对比CPML与PML对大角度入射隐失波的吸收效果,模型大小为500 m×2 000 m,模型中纵波速度v_p=3 200 m/s,横波速度v_s=1 847.5 m/s,密度ρ=2 500 kg/m³,网格单元大小为Δx=Δz=5 m,时间步长为0.5 ms,记录长度为500 ms,震源为30 Hz的零相位雷克子波,位于模型上边界附近,加载CPML与PML两种边界条件得到的波场快照如图7所示。图中可见,CPML吸收边界算法对入射到边界层内波的能量进行了充分的吸收,而不像PML在上边界附近有规律出现的隐失波不能得到充分地吸收。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 加载CPML与PML两种边界条件的波场快照对比

a—PML吸收边界200 ms波场快照;b—PML吸收边界300 ms波场快照;c—PML吸收边界400 ms波场快照;d—CPML吸收边界200 ms波场快照;e—CPML吸收边界300 ms波场快照;f—CPML吸收边界400 ms波场快照

3 各向异性及随机介质模型正演实例

3.1 各向异性层状模型

建立图8所示大小为2 500 m×1 000 m各向异性VTI介质层状模型,模型参数如表1,介质第一、二层均为各向异性VTI介质。应用旋转交错网格高阶差分法结合CPML边界条件对该模型进行了正演,计算时加载主频为30 Hz的雷克子波,坐标为(x=1 000 m,z=0 m),接收排列位于地面z=0 m,空间步长为Δx=Δz=2.5 m,时间步长为Δt=0.2 ms,记录长度2 s,吸收边界参数选择m=3,R=10^-6,κ_max=3,α_max=120,此时正演的x、z分量剖面如图9所示。而使用PML吸收边界处理的波场快照(图10),在各向异性VTI介质中弹性波存在三叉区现象,加之转换波的存在,使得剖面记录变得复杂,并且由于各向异性的存在,使反射波和绕射波的形状发生了很大改变,容易造成地震数据的误解释。

表1 各向异性断层介质模型中的弹性参数(C_ij/GPa)

介质层	C₁₁	C₁₃	C₃₃	C₄₄	ε	δ	ρ/(kg/m³)
Ⅰ	18.1	7.11	9.68	1.7	0.43	0.09	2250
Ⅱ	25.5	5.0	18.4	5.6	0.19	-0.11	2440

新窗口打开| 下载CSV

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 各向异性VTI介质层状模型

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 各向异性VTI介质层状CPML边界条件下的剖面

a—x分量剖面;b—z分量剖面;1—反射qP波;2—反射转换qSV波;5—反射qSV波;6—反射转换qP波

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 各向异性VTI介质层状PML边界条件下的剖面

a—x分量剖面;b—z分量剖面;1—反射qP波;2—反射转换qSV波;5—反射qSV波;6—反射转换qP波

分析图11波场快照可知,随着时间的增加,qP波波场的变化平缓简单,qSV波波场变化剧烈,qP波的波前面呈现椭圆形状,qSV波出现三叉区,x分量关于过震源点的z方向直线相位相反,z分量关于过震源点的z方向直线相位相同,当qP波遇到波阻抗分界面时会发生反射和透射现象,除了反射qP波和透射qP波,还会产生反射转换qSV波和透射转换qSV波;同样地,qSV波遇到波阻抗分界面时也会产生反射qSV波、透射qSV波、反射转换qP波和透射转换qP波。再加上qSV波的波前面存在三叉区现象,使得波场变得十分复杂。

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 各向异性介质断层模型的正演波场快照

a—各向异性VTI层状介质x分量300 ms波场快照;b—各向异性VTI层状介质z分量300 ms波场快照;c—各向异性VTI层状介质x分量500 ms波场快照;d—各向异性VTI层状介质z分量500 ms波场快照;1—反射qP波;2—反射转换qSV波;3—透射转换qSV波;4—透射qP波;5—反射qSV波;6—反射转换qP波;7—透射qSV波;8—透射转换qP波

从图9可知,单炮记录十分清晰,基本没有数值频散且没有人工边界反射,在单炮记录上可以看到直达波、反射波、转换波及qSV波的三叉区,并且可以观测到x分量的相位相反现象。

3.2 随机介质分层模型

为了更真实地模拟实际地下介质,采用指数型椭圆自相关函数来构造随机介质,其函数表达式为:

(18)

φ (x, z) = \exp [- \sqrt[]{(x^{2} / a^{2} + z^{2} / b^{2})}],

式中,a、b分别表示介质在x和z方向上的自相关长度。构造步骤如下:

1)选择指数型椭圆自相关函数式,选择自相关长度a、b,再求得自相关函数φ(x,z)的二维傅里叶变换ϕ(k_x,k_z),即随机扰动σ(x,z)的功率谱;

2)由功率谱函数ϕ(k_x,k_z),根据随机过程方法产生对应的随机谱函数;

3)应用随机过程的谱展开公式,产生由给定的自相关函数φ(x,z)描述的随机扰动σ(x,z);

4)将随机扰动σ(x,z)规范化,从而得到以 φ(x,z)为自相关函数、具有指定均值及方差的随机介质模型。

图12、13、14分别为均匀层状介质模型、上层为参数a=b=1随机介质模型、上层为参数a=b=10随机介质模型,模型大小均为400 m×300 m,震源主频分别为20 Hz和40 Hz的雷克子波,炮点坐标(x=200 m,z=0 m),接收排列位于地面z=0 m,网格大小为Δx=Δz=1 m,计算时间步长为Δt=0.5 ms,记录长度1 s,密度为常数。

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12 均匀层状介质模型

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13 随机介质模型(a=b=1)

图14

新窗口打开| 下载原图ZIP| 生成PPT

图14 随机介质模型(a=b=10)

图15为震源主频为20 Hz的雷克子波z分量单炮剖面记录,其中图15a为主频为20 Hz的雷克子波在均匀介质第200道共炮点剖面记录,图15b为主频为20 Hz的雷克子波在随机介质(a=b=1)第200道共炮点剖面记录,图15c为主频为20 Hz的雷克子波在随机介质(a=b=10)第200道共炮点剖面记录。

图15

新窗口打开| 下载原图ZIP| 生成PPT

图15 主频为20 Hz的z分量第200道共炮点剖面记录

a—均匀层状介质剖面记录;b—随机层状介质(a=b=1)剖面记录;c—随机层状介质(a=b=10)剖面记录;1—反射P波;2—反射转换SV波;3—反射转换P波;4—反射SV波

图16为震源主频为40 Hz的雷克子波z分量单炮剖面记录,其中图16a为主频为40 Hz的雷克子波在均匀介质第200道共炮点剖面记录,图16b为主频为40 Hz的雷克子波在随机介质(a=b=1)第200道共炮点剖面记录,图16c为主频为40 Hz的雷克子波在随机介质(a=b=10)第200道共炮点剖面记录。

图16

新窗口打开| 下载原图ZIP| 生成PPT

图16 主频为40 Hz的z分量第200道共炮点剖面记录

a—均匀层状介质剖面记录;b—随机层状介质(a=b=1)剖面记录;c—随机层状介质(a=b=10)剖面记录;1—反射P波;2—反射转换SV波;3—反射转换P波;4—反射SV波

图17为随机模型z分量正演模拟波场快照,其中图17a为在500 ms时主频为20 Hz的雷克子波在均匀介质的震源正演模拟波场,图17b为在500 ms时主频为20 Hz的雷克子波在随机介质(a=b=1)的震源正演模拟波场,图17c为在500 ms时主频为20 Hz的雷克子波在随机介质(a=b=10)的震源正演模拟波场。

图17

新窗口打开| 下载原图ZIP| 生成PPT

图17 随机介质模型正演波场快照(20 Hz)

a—均匀层状介质500 ms波场快照;b—随机层状介质(a=b=1)500 ms波场快照;c—随机层状介质(a=b=10)500 ms波场快照;1—投射转换P波;2—透射SV波;3—反射SV波;4—反射转换P波;5—反射转换SV波;6—直达SV波

图18为随机模型z分量正演模拟波场快照,其中图18a为在500 ms时主频为40 Hz的雷克子波在均匀介质的震源正演模拟波场,图18b为在500 ms时主频为40 Hz的雷克子波在随机介质(a=b=1)的震源正演模拟波场,图18c为在500 ms时主频为40 Hz的雷克子波在随机介质(a=b=10)的震源正演模拟波场。

图18

新窗口打开| 下载原图ZIP| 生成PPT

图18 随机介质模型正演波场快照(40 Hz)

通过纽曼稳定性条件(式(4)、(5))分别算出,当模型介质纵波速度v_p=1 200 m/s,空间步长Δh=1 m时稳定的标准交错网格的时间步长为Δt≤0.43 ms,稳定的旋转交错网格的时间步长为Δt≤0.60 ms,下面实验证明了这个结论。

使用标准交错网格算法,取上层为参数a=b=1随机介质模型,模型大小均为400 m×300 m(图13模型),震源主频为40 Hz的雷克子波,炮点坐标 (x=200 m,z=0 m),接收排列位于地面z=0 m,网格大小为Δx=Δz=1 m,计算时间步长为Δt=0.5 ms,记录长度1 s,密度为常数。图19为标准交错网格和旋转交错网格当时间步长取Δt=0.5 ms时,在 88 ms 时的正演模拟波场对比。

图19

新窗口打开| 下载原图ZIP| 生成PPT

图19 标准交错网格与旋转交错网格正演波场快照对比

a—标准交错网格88 ms波场快照;b—旋转交错网格88 ms波场快照

由于取时间步长为Δt=0.5 ms,超出了标准交错网格的时间步长取值范围,因此,在随机介质与均匀介质的分界面处导致算法不稳定,出现了数值发散的现象,相比于旋转交错网格而言,时间步长的取值较为宽松,因此算法稳定,不会出现数值发散的现象。

通过以上两个例子验证了使用传统的标准交错网格进行正演模拟时,当界面存在波阻抗突变时,会使差分算法不稳定,导致数值发散的现象,而解决这一问题的另一个方法是使用更小的时间步长,从而降低了算法的效率,因此采用旋转交错网格,在增加算法稳定性的同时提高算法的效率,采用常规的PML边界条件则会导致吸收效果不佳,不能对大角度入射的隐失波充分吸收,采用CPML边界条件不仅充分吸收了隐失波,而且进一步提高了整体的吸收效果。

分析以上结果可知,雷克子波主频越高则分辨率也越高,当自相关长度很小(a=b=1) 时,反射波同相轴连续性较好,随着自相关长度的增加(a=b=10),反射波同相轴的连续性逐渐变差,甚至错断;同时,当自相关长度很小(a=b=1) 时,随机介质的非匀质性对直达波波前面几乎没有影响,随着自相关长度的增大(a=b=10),直达波波前面开始不再呈球面扩散;然而,当自相关长度很小(a=b=1) 时,波的散射现象严重,波形杂乱且反射波形受到散射波的严重覆盖,随着自相关长度的增加(a=b=10),弹性波具有较轻微的散射现象,反射波形受散射波的影响较小。基于以上原因,均匀介质模型的剖面记录中的直达波和反射波十分清晰,而随机介质模型的自相关长度较小(a=b=1)时,弹性波记录的同相轴连续性较好,但是由于波的散射现象严重,波形更为杂乱无章,反射波形受散射波影响严重,反射波形记录很难辨认出来,随着自相关长度的增大(a=b=10),弹性波记录的同相轴连续性变差,但是由于弹性波具有较轻微的散射现象,波形较小尺度随机介质更为清晰,反射波形受散射波影响较小,部分反射波形记录能够很容易地辨认出来。随机介质模型的剖面记录较均匀介质模型更为复杂且更难解释,然而随机介质显然更符合地下介质的实际情况,以随机介质的波形模拟结果能更好地指导地震资料的解释。

4 结论

1) 旋转交错网格在标准网格交错网格的基础上增加了网格旋转,将同一物理量的不同分量定义在同一网格点上。因此,在模拟非均匀性比较强或各向异性介质时,不需要对介质参数进行插值处理,在遇到自由界面时也不需要进行单独的显式处理,从而提高了算法的稳定性。

2) 通过分析接收排列关于参数m反射误差分布图,参数m=3,κ=3吸收边界效果最优,随着α_max的增大吸收效果会逐渐变好,当α_max=4f₀时已经可以取得很好的吸收效果。而且CPML为基于复频移和卷积计算技术的完全匹配层,对隐失波,特别是低频隐失波,较PML有着更出色的吸收效果。

3) 旋转交错网格高阶差分法能够有效模拟弹性波在各向异性介质与随机介质中的波场传播特征,可以作为复杂介质波传规律研究的有效数值工具。相比各向同性介质,各向异性介质中的反射波和绕射波的形状发生了很大畸变,可清晰观察到纵横波分裂与方位特征差异;而通过增加随机介质中自相关长度,地震记录将变得更加杂乱。由于实际的地震资料中,绝大部分介质是随机、各向异性的,通过对随机介质和各向异性介质的正演模拟,通过地震记录和波场快照获得了弹性波在这两种介质中的波形特征,分析了为地震剖面解释带来的影响,对实际地震资料的解释做了初步的指导。

附录

A)时间二阶精度空间2M阶(本文取16阶)精度有限差分的离散格式

(1)

\{\begin{array}{l} U_{i, j}^{k + 1 / 2} = U_{i, j}^{k - 1 / 2} + [{\overset{̅}{L}}_{x} (P_{i, j}^{k}) + {\overset{̅}{L}}_{z} (S_{i, j}^{k})] Δ t / ρ_{i, j}, \\ W_{i, j}^{k + 1 / 2} = W_{i, j}^{k - 1 / 2} + [{\overset{̅}{L}}_{x} (S_{i, j}^{k}) + {\overset{̅}{L}}_{z} (Q_{i, j}^{k})] Δ t / ρ_{i, j}, \\ P_{i, j}^{k + 1} = P_{i, j}^{k} + Δ t [C_{11 i, j} {\overset{̅}{L}}_{x} (U_{i, j}^{k + 1 / 2}) + C_{15 i, j} {\overset{̅}{L}}_{x} (W_{i, j}^{k + 1 / 2})] + Δ t [C_{13 i, j} {\overset{̅}{L}}_{z} (W_{i, j}^{k + 1 / 2}) + C_{15 i, j} {\overset{̅}{L}}_{z} (U_{i, j}^{k + 1 / 2})], \\ Q_{i, j}^{k + 1} = Q_{i, j}^{k} + Δ t [C_{13 i, j} {\overset{̅}{L}}_{x} (U_{i, j}^{k + 1 / 2}) + C_{35 i, j} {\overset{̅}{L}}_{x} (W_{i, j}^{k + 1 / 2})] + Δ t [C_{33 i, j} {\overset{̅}{L}}_{z} (W_{i, j}^{k + 1 / 2}) + C_{35 i, j} {\overset{̅}{L}}_{z} (U_{i, j}^{k + 1 / 2})], \\ S_{i, j}^{k + 1} = S_{i, j}^{k} + Δ t [C_{15 i, j} {\overset{̅}{L}}_{x} (U_{i, j}^{k + 1 / 2}) + C_{55 i, j} {\overset{̅}{L}}_{x} (W_{i, j}^{k + 1 / 2})] + Δ t [C_{35 i, j} {\overset{̅}{L}}_{z} (W_{i, j}^{k + 1 / 2}) + C_{55 i, j} {\overset{̅}{L}}_{z} (U_{i, j}^{k + 1 / 2})] 。 \end{array}

式(1)中:

(2)

{\overset{̅}{L}}_{x} [u {(x)}_{i, j}] = {\overset{N}{\sum_{n = 1}} K}_{n}^{(N)} [u {(x)}_{i + \frac{2 n - 1}{2}, j + \frac{2 n - 1}{2}} - u {(x)}_{i - \frac{2 n - 1}{2}, j - \frac{2 n - 1}{2}}] / 2 Δ x + \overset{N}{\sum_{n = 1}} K_{n}^{(N)} [u {(x)}_{i + \frac{2 n - 1}{2}, j - \frac{2 n - 1}{2}} - u {(x)}_{i - \frac{2 n - 1}{2}, j + \frac{2 n - 1}{2}}] / 2 Δ x,

(3)

{\overset{̅}{L}}_{z} [u {(x)}_{i, j}] = \overset{N}{\sum_{n = 1}} K_{n}^{(N)} [u {(z)}_{i + \frac{2 n - 1}{2}, j + \frac{2 n - 1}{2}} - u {(z)}_{i - \frac{2 n - 1}{2}, j - \frac{2 n - 1}{2}}] / 2 Δ z - \overset{N}{\sum_{n = 1}} K_{n}^{(N)} [u {(z)}_{i + \frac{2 n - 1}{2}, j - \frac{2 n - 1}{2}} - u {(z)}_{i - \frac{2 n - 1}{2}, j + \frac{2 n - 1}{2}}] / 2 Δ z 。

式(2)、(3)中差分系数 $K_{n}^{(N)}$ 与标准交错网格的空间差分系数求取方法一致^[1,2]。

B)CPML吸收边界下的弹性波波动方程:

(4)

\{\begin{array}{l} ρ \frac{\partial V_{x}}{\partial t} = \frac{\partial τ_{xx}}{\partial x} + ψ_{xx} + \frac{\partial τ_{xz}}{\partial z} + ψ_{xz}, \\ ρ \frac{\partial V_{z}}{\partial t} = \frac{\partial τ_{zx}}{\partial x} + ψ_{zx} + \frac{\partial τ_{zz}}{\partial z} + ψ_{zz}, \\ \frac{\partial τ_{xx}}{\partial t} = C_{11} (\frac{\partial V_{x}}{\partial x} + Ω_{xx}) + C_{13} (\frac{\partial V_{z}}{\partial z} + Ω_{zz}) + C_{15} (\frac{\partial V_{z}}{\partial x} + Ω_{zx} + \frac{\partial V_{x}}{\partial z} + Ω_{xz}), \\ \frac{\partial τ_{zz}}{\partial t} = C_{31} (\frac{\partial V_{x}}{\partial x} + Ω_{xx}) + C_{33} (\frac{\partial V_{z}}{\partial z} + Ω_{zz}) + C_{35} (\frac{\partial V_{z}}{\partial x} + Ω_{zx} + \frac{\partial V_{x}}{\partial z} + Ω_{xz}), \\ \frac{\partial τ_{zx}}{\partial t} = C_{51} (\frac{\partial V_{x}}{\partial x} + Ω_{xx}) + C_{53} (\frac{\partial V_{z}}{\partial z} + Ω_{zz}) + C_{55} (\frac{\partial V_{z}}{\partial x} + Ω_{zx} + \frac{\partial V_{x}}{\partial z} + Ω_{xz}) 。 \end{array}

式(4)中:

(5)

ψ_{ij} = b_{j} ψ_{ij} + a_{j} \cdot \partial τ_{ij} / \partial j, Ω_{ij} = b_{j} Ω_{ji} + a_{j} \cdot \partial V_{i} / \partial j 。

The authors have declared that no competing interests exist.

作者已声明无竞争性利益关系。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

董良国, 马在田, 曹景忠 , 等.

一阶弹性波方程交错网格高阶差分解法

[J]. 地球物理学报, 2000,43(3):411-419.

[本文引用: 3]

[2]

董良国, 马在田, 曹景忠 .

一阶弹性波方程交错网格高阶差分解法稳定性研究

[J]. 地球物理学报, 2000,43(6):856-864.

[本文引用: 2]

[3]

Madariaga

Dynamics of an expanding circular fault

[J]. Bulletin of the Seismological Society of America, 1976,66(3):163-182.

[本文引用: 2]

[4]

Virieux

P-SV wave propagation in heterogeneous media:Velocity stress finite-difference method

[J]. Geophysics, 1986,51(4):889-901.

DOI:10.1190/1.1442147 URL [本文引用: 1]

[5]

Virieux

SH-wave propagation in heterogeneous media:Velocity stress finite-difference method

[J]. Geophysics, 1984,51(4):1933-1942.

[本文引用: 1]

[6]

裴正林

双相各向异性介质弹性波传播交错网格高阶有限整分法模拟

[J]. 石油地球物理勘探, 2006,41(2):137-143.

DOI:10.3321/j.issn:1000-7210.2006.02.005 URL [本文引用: 1]

含流体孔隙各向异性介质中的弹性波传播问题是目前石油勘探和地震学研究的难点之一。基于Biot理论,本文提出了二维双相任意倾斜各向异性介质三分量弹性波方程交错网格任意偶阶精度有限差分解法,并对均匀及两层双相VTI介质和TTI介质中的弹性波场进行了模拟。波场快照和合成记录表明:当横波源的偏振方向与TI介质的对称轴存在一定夹角、且耗散系数较小时,在固相与流相中存在快横波、慢横波、快纵波和慢纵波;双相TI介质中横波分裂现象及波前面尖角和三分叉现象与TI介质中的基本相同,只是黏滞相界中弹性波速度有所减小,且具有频散和衰减特性;从直达波转换为反射波后,该波型具有反射波的性质。

[7]

裴正林

三维双相各向异性介质弹性波方程交错网格高阶有限差分法模拟

[J]. 中国石油大学学报:自然科学版, 2006,30(2):16-20.

DOI:10.3321/j.issn:1000-5870.2006.02.004 URL [本文引用: 1]

基于Biot理论,给出了三维双相各向异性介质应力-速度弹性波方程交错网格任意偶阶精度有限差分解法,并对三维双相各向异性(TIH)介质中弹性波场进行了模拟.结果表明,弹性波在三维双相各向异性介质中传播时存在快纵波qP1、快横波qS1、慢横波qS2和慢纵波qP2,并清楚地观测到了横波分裂、横波分裂盲点、波面三分叉等特殊现象.快纵波在固相和流相中的相位相同;而慢纵波在固相和流相中的相位相反,即慢纵波在流相中振幅大,而在固相中的振幅较小.另外,三维双相各向异性介质中弹性波场的快纵波和快横波的耦合关系、波的类型、能量分布和相位等都是在三维空间中变化的.

[8]

孙卫涛, 杨慧珠 .

双相各向异性介质弹性波场有限差分正演模拟

[J]. 固体力学学报, 2004,25(1):21-28.

DOI:10.3969/j.issn.0254-7805.2004.01.005 URL [本文引用: 1]

从双相各向异性介质模型出发,以Boit理论为基础,推导了斜方晶系各向异性介质一阶弹性波动方程.引入固、流体密度比和孔隙几何参数,将Biot方程系数简化为测量简单、物理意义明确的物理量.采用交错网格技术建立了各向异性孔隙介质波动方程的高精度差分格式,并首次对这类差分格式的频散特性和稳定性作了详细分析讨论,解决了计算稳定性和边界反射问题.与解析解的对比以及理论模型的数值模拟都表明,该方法不仅大大降低了计算量,提高了正演速度,并且具有良好的稳定性和精确性.

[9]

陈祖银, 王用军, 彭达 .

弹性波波场分离数值模拟

[J]. 中国石油勘探, 2014,19(6):47-53.

[本文引用: 1]

[10]

Chen K

Finite-difference simulation of elastic wave with separation in pure P-and S-modes

[J]. Journal of Computation methods in Physics, 2014,( 3):1-14.

DOI:10.1155/2014/108713 URL [本文引用: 1]

Elastic wave equation simulation offers a way to study the wave propagation when creating seismic data. We implement an equivalent dual elastic wave separation equation to simulate the velocity, pressure, divergence, and curl fields in pure P- and S-modes, and apply it in full elastic wave numerical simulation. We give the complete derivations of explicit high-order staggered-grid finite-difference operators, stability condition, dispersion relation, and perfectly matched layer (PML) absorbing boundary condition, and present the resulting discretized formulas for the proposed elastic wave equation. The final numerical results of pure P- and S-modes are completely separated. Storage and computing time requirements are strongly reduced compared to the previous works. Numerical testing is used further to demonstrate the performance of the presented method.

[11]

奚先, 姚姚 .

二维弹性随机介质中的波场特征

[J].地球物理学进展, 2005, 20(1): 147-154 .

[本文引用: 1]

[12]

奚先, 姚姚 .

随机介质模型的模拟与混合型随机介质

[J]. 中国地质大学学报:地球科学, 2002,27(1):67-71.

DOI:10.3321/j.issn:1000-2383.2002.01.014 URL [本文引用: 1]

讨论了随机介质模型的基本概念及指数型和高斯型椭圆自相关函数所描述的随机介质模型的特点,并提出了混合型随机介质模型的概念,该随机介质模型能更加灵活、准确地描述实际介质.通过选择在水平方向和垂直方向上的自相关长度a、b以及粗糙度r,可以产生出各种不同形式的混合型随机介质模型.模拟结果显示,混合型随机介质模型能更加灵活地描述实际介质,具有适应性强,使用方便、灵活,能有效地模拟油气藏细节的优点.

[13]

奚先, 姚姚 .

二维随机介质及波动方程正演模拟

[J]. 石油地球物理勘探, 2001,31(5):546-552.

DOI:10.3321/j.issn:1000-7210.2001.05.005 URL [本文引用: 2]

本文介绍了随机介质模型的基本概念,提出了一种根据随机过程的谱分解理论构造随机介质模型的新方法;并以二维指数型椭圆自相关函数为例进行了若干二维随机介质模型的正演模拟.结果显示,这些随机介质模型能灵活地描述实际介质.本文还利用有限差分法正演模拟了随机介质模型中弹性波的传播及其自激自收时间记录.正演计算结果表明,随机介质模型理论具有适应性强,使用方便、灵活,能有效地模拟油气藏细节的优点,是油藏(或油储)地球物理中值得进一步深入研究的课题,也是重要的发展方向.

[14]

陈浩, 王秀明, 赵海波 .

旋转交错网格有限差分及其完全匹配层吸收边界条件

[J]. 科学通报, 2006,51(17):1950-1994.

DOI:10.3321/j.issn:0023-074X.2006.17.002 URL [本文引用: 1]

将完全匹配层吸收边界条件引入到旋转网格有限差分中，用以解决在非均匀弹性和孔隙弹性介质情况下数值模拟中的吸收边界问题，同时，将旋转交错网格有限差分法用于数值求解等效弹性介质、孔隙弹性介质和各向异性弹性介质的波动方程．与普通的交错网格有限差分方法相比，旋转交错网格有限差分的好处在于不同的物理量只位于两个不同的位置：应力和应变（或质点速度和位移）位于离散单元的中心，而质点速度或位移（或应力和应变）位于单元的顶点．经过这样的处理后，弹性常数就不再需要进行平均（模拟非均匀介质时）或内插（模拟各向异性介质时），因为此时所有的弹性模量都位于相同的位置，且与应力或应变的位置相对应．为了验证新算法，对相同的模型采用了不同的算法计算和比较．研究结果表明，旋转交错网格有限差分算法和普通交错网格有限差分算法的结果吻合很好．不仅如此，新算法能很方便地处理声阻抗差别较大的非均匀介质，特别是在模拟充有液体或气体的裂缝介质时更具优势．另一方面，应用完全匹配层吸收边界条件可以大大减少人工界面产生的反射波．当边界层的厚度超过半个波长时，在人工界面几乎无反射波产生．理论和数值模拟结果表明，旋转网格中的完全匹配层吸收条件与普通交错网格中的吸收效果和处理方法几乎相同．此外，建立了在等效弹性、孔隙弹性和各向异性弹性介质中的完全匹配层吸收条件的速度-应力差分方程系统．

[15]

Berenger J

A perfectly matched layer for the absorption of electromagnetic waves

[J]. Journal of Computational Physics, 1994,114(2):185-200.

DOI:10.1006/jcph.1994.1159 URL [本文引用: 1]

[16]

Berenger J

Three-dimensional perfectly matched layer for the absorption of electromagnetic waves

[J]. Journal of Computational Physics, 1996,127(2):363-379.

DOI:10.1006/jcph.1996.0181 URL [本文引用: 1]

[17]

, Zeng Z

, Huang

, et al.

GPR simulation based on complex frequency shifted recursive integration PML boundary of 3D high order FDTD

[J]. Computers & Geosciences, 2012,49(4):121-130.

[本文引用: 1]

[18]

Hastings F

, Schneider J

, Broschat S

Application of the perfectly matched layer (PML) absorbing boundary condition toelastic wave propagation

[J]. Journal of the Acoustical Society of America, 1996,100(5):3061-3069.

DOI:10.1121/1.417118 URL [本文引用: 1]

[19]

Drossaert F

, Giannopoulos

A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves

[J]. Geophysics, 2007,72(2):9-17.

[本文引用: 1]

[20]

Drossaert F

, Giannopoulos

Complex frequency shifted convolution PML for FDTD modeling of elastic waves

[J]. Wave Motion, 2007,44(7-8):593-604.

DOI:10.1016/j.wavemoti.2007.03.003 URL [本文引用: 1]

The perfectly matched layer (PML) is nowadays considered as the best optimum absorbing boundary condition available. However, the PML with the classical stretching tensor has certain limitations. Strangely, these limitations have rarely been addressed in elastic wave modelling. For example, substantial reflections occur when strong evanescent waves are propagating parallel to the interface. To circumvent problems like this, the complex frequency shifted stretching tensor has been introduced in electromagnetic modelling. In this paper we show that the convolution PML with this stretching tensor as used in electromagnetic modelling can be adapted for elastic wave modelling. Numerical results of a model where the presence of evanescent waves is predominant show that the PML based on the complex frequency shifted stretching tensor can improve the performance of the absorbing boundary layer considerably.

[21]

Zhang

, Shen

Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling

[J]. Geophysics, 2010,75(4):141-154.

DOI:10.1190/1.3463431 URL [本文引用: 1]

[22]

Martin

, Komatitsch

, Gedney S

, et al.

A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML)

[J]. CMES, 2010,56(1):17-40.

[本文引用: 1]

[23]

张显文, 韩立国, 黄玲 , 等.

基于递归积分的复频移PML弹性波方程交错网格高阶差分法

[J]. 地球物理学报, 2009,52(7):1800-1807.

DOI:10.3969/j.issn.0001-5733.2009.07.014 URL [本文引用: 2]

常规完全匹配吸收边界（PML）对以近掠射角入射到界面上的波以及低频波、损耗波都会产生虚假边界反射．基于递归积分的不分裂复频移PML算法，利用复频移拉伸函数，极大地改善了PML边界条件的性能，我们进一步推导出基于递归积分的不分裂复频移PML弹性波方程交错网格高阶差分法，对长条形介质模型进行数值模拟，与常规PML算法进行比较说明该算法对以掠射角入射到PML界面的波以及PML层内损耗波的吸收效果．

[24]

Roden

, Gedney

S D

Convolution(PMLCPML):an effcient FDTD implementation of the CFS-PML for arbitrary media

[J]. Microwave & Optical Technology Letters, 2000,27(5):334-339.

[本文引用: 1]

[25]

张鲁新, 符力耘, 裴正林 .

不分裂卷积完全匹配层与旋转交错网格有限差分在孔隙弹性介质模拟中的应用

[J]. 地球物理学报, 2010,53(10):2470-2483.

DOI:10.3969/j.issn.0001-5733.2010.10.021 URL [本文引用: 1]

完全匹配层吸收边界在地震波模拟中已广泛使用，但常用的场分裂格式完全匹配层吸收边界（SPML）和传统的不分裂完全匹配层吸收边界（NPML）X4极低频入射波或大角度入射波的边界吸收效果不好．一种无需分裂和显式卷积计算的完全匹配层吸收边界（CPML）不仅能够解决常规PML吸收边界的不足，而且具有存储量小、计算效率高、易于编程实现的特点．本文将这种完全匹配层（CPML）吸收边界引入到孔隙弹性介质速度一应力格式的旋转交错网格有限差分算法中，对完全匹配层吸收边界参数进行数值分析，得到一组优化的参数．孔隙弹性介质数值模拟结果表明这种不分裂卷积完全匹配层的吸收效果优于常规完全匹配层．

[26]

Saenger E

, Gold

, Shapiro S

Modeling the propagation of elastic waves using a modified finite-difference grid

[J]. Wave Motion, 2000,31(1):77-92.

DOI:10.1016/S0165-2125(99)00023-2 URL [本文引用: 2]

The modeling of elastic waves with an explicit finite difference (FD) scheme on a staggered grid causes instability problems when the medium possesses high contrast discontinuities (strong heterogeneities). In this paper we have derived a new rotated staggered grid where all medium parameters are defined at appropriate positions within an elementary cell for the essential operations. Using this modified grid it is possible to simulate the propagation of elastic waves in a medium containing cracks, pores or free surfaces without applying boundary conditions. We compare the von Neumann stability criterion and the dispersion error for the new rotated staggered grid with the results of the standard staggered grid. Additionally, we show two synthetic examples and a comparison with a laboratory experiment to demonstrate advantages of the new rotated staggered grid in 2D and 3D.

[27]

Kuzuoglu

, Mittra

Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers

[J]. IEEE Microwave and Guided Wave Letters, 1996,6(12):447-449.

DOI:10.1109/75.544545 URL [本文引用: 1]

一阶弹性波方程交错网格高阶差分解法

2000

... 弹性波正演是人们认识地下复杂介质的重要手段,在弹性波有限差分正演模拟中,除了二阶微分方程外,还常常采用一阶速度—应力方程,它的主要优点是不需对弹性常数进行空间微分^[1,2],正演计算能得到完整的波场响应.除了规则网格有限差分法外,Madariaga^[3]提出一种较为先进的交错网格,Virieux^[4,5]将其应用到各向同性介质中的SH波和P-SV波模拟中,它在不增加计算工作量和存储空间的前提下,和常规网格相比局部精度提高4倍,收敛速度也很快.裴正林^[6,7]、孙卫涛和杨慧珠^[8]将交错网格高阶有限差分法应用到双相各向异性介质弹性波传播中;陈祖银等^[9]、Chen^[10]将高阶交错网格有限差分法应用到弹性波场分离的数值模拟中,以便能更精确地研究纵波、横波在介质的传播规律;奚先和姚姚^[11,12,13]将交错网格有限差分法应用于二维弹性随机介质的正演及波场特征研究;但是,与标准交错网格有限差分方法相比,旋转交错网格有限差分的好处在于不同的物理量只位于两个不同的位置:应力和应变(或质点速度和位移)位于离散单元的中心,而质点速度或位移(或应力和应变)位于单元的顶点.经过这样的处理后,弹性常数就不再需要进行平均(模拟非均匀介质时)或内插(模拟各向异性介质时),因为此时所有的弹性模量都位于相同的位置,且与应力或应变的位置相对应^[14]. ...

... 一阶速度—应力弹性波方程不需对弹性常数进行空间微分,且方程中仅含有速度和应力变量关于时间和空间的一阶偏导数,容易提高差分近似的阶数.为此,将二维各向异性介质及随机介质的正演控制方程表示为^[1]: ...

... 式(2)、(3)中差分系数

K_{n}^{(N)}

与标准交错网格的空间差分系数求取方法一致^[1,2]. ...

一阶弹性波方程交错网格高阶差分解法稳定性研究

2000

... 式(2)、(3)中差分系数

K_{n}^{(N)}

与标准交错网格的空间差分系数求取方法一致^[1,2]. ...

Dynamics of an expanding circular fault

1976

... 式(1)中V_x与V_z为别为x方向和z方向的质点速度分量,τ_xx与τ_zz为正应力,τ_xz与τ_zx为切应力,ρ为密度,C_ij(i,j=1,3,5)为介质的弹性参数.目前普遍采用图1a所示的标准交错网格有限差分^[3]对式(1)进行数值模拟.但是,当遇到波阻抗较大的界面时,在自由界面处,边界条件必须在差分算法中显式表示出来,标准交错网格有限差分变得不再稳定.Saenger等人^[26]为了改善标准交错网格的性能,提出了图1b所示的旋转交错网格,图中实线和虚线分别表示空间上的整网格和半网格.旋转交错网格在标准网格交错网格的基础上增加了网格旋转,将同一物理量的不同分量定义在同一网格点上,其中速度分量和密度分量定义在整网格点,应力分量和弹性常量定义在半网格点上.因此,在模拟非均匀性比较强或各向异性介质时,不需要对介质参数进行插值处理,在遇到自由界面时也不需要进行单独的显式处理. ...

P-SV wave propagation in heterogeneous media:Velocity stress finite-difference method

1986

SH-wave propagation in heterogeneous media:Velocity stress finite-difference method

1984

双相各向异性介质弹性波传播交错网格高阶有限整分法模拟

2006

三维双相各向异性介质弹性波方程交错网格高阶有限差分法模拟

2006

双相各向异性介质弹性波场有限差分正演模拟

2004

弹性波波场分离数值模拟

2014

Finite-difference simulation of elastic wave with separation in pure P-and S-modes

2014

二维弹性随机介质中的波场特征

2005

随机介质模型的模拟与混合型随机介质

2002

二维随机介质及波动方程正演模拟

2001

... 为了提高差分精度,减小网格频散,笔者将速度(应力)对时间的奇数阶高阶导数转化为应力(速度)对空间的导数,运用时间二阶精度、空间十六阶精度高阶差分法,由于随机介质中存在大量的波阻抗差较大的界面,并且需要尽可能减小界面反射,消除隐失波对波剖面的影响,因此采用标准交错网格结合常规PML显然不适合模拟随机介质,更为主要的是在截断界面处的波反射会干扰波剖面的分析^[13],因此需要引入旋转交错网格结合递归卷积复频移PML技术,开展各向异性介质与随机介质模型的弹性波正演. ...

旋转交错网格有限差分及其完全匹配层吸收边界条件

2006

A perfectly matched layer for the absorption of electromagnetic waves

1994

... 应用有限差分法开展弹性波模型时,需要在截断人工边界处引入吸收边界条件.目前吸收效果较好的边界条件为完全匹配层(PML)边界^[15,16,17],它可以使弹性波在匹配层中无反射地通过并按指数规律衰减,适合于吸收较大范围的入射角.Hastings等^[18]将PML方法应用于弹性波的有限差分正演中;Drosseart^[19,20]提出了一种基于递归积分复频移PML算法;Zhang和Shen^[21]、Martin等^[22] 及张显文等^[23]将基于辅助微分方程的复频移PML应用到弹性波数值模拟中,该方法易于扩展到高阶差分形式;Roden和Gedney^[24]应用卷积技术实现了完全匹配层加载,大大提高了运算效率;张鲁新等^[25]将不分裂的卷积完全匹配层(CPML)结合旋转交错网格有限差分技术应用到孔隙弹性介质模拟中.使用常规的PML算法处理截断边界时仍会留下一些反射波和大量的隐失波,而采用CPML算法处理截断边界则可以减少反射波的存在,并且对隐失波的吸收具有很好的效果^[23]. ...

Three-dimensional perfectly matched layer for the absorption of electromagnetic waves

1996

GPR simulation based on complex frequency shifted recursive integration PML boundary of 3D high order FDTD

2012

Application of the perfectly matched layer (PML) absorbing boundary condition toelastic wave propagation

1996

A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves

2007

Complex frequency shifted convolution PML for FDTD modeling of elastic waves

2007

Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling

2010

A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML)

2010

基于递归积分的复频移PML弹性波方程交错网格高阶差分法

2009

... [23]. ...

Convolution(PMLCPML):an effcient FDTD implementation of the CFS-PML for arbitrary media

2000

不分裂卷积完全匹配层与旋转交错网格有限差分在孔隙弹性介质模拟中的应用

2010

Modeling the propagation of elastic waves using a modified finite-difference grid

2000

... 根据CPML边载方式可将整个模拟区域划分为内部计算区域和CPML区域.为了使加载的CPML吸收效果最好,需要进一步探讨吸收边界参数的选取.根据文献[26]对σ_p的设置方法为: ...

Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers

1996

... 根据Kuzuoglu和Mittra^[27]的研究,定义对应坐标的微分算子和复频移拉伸函数分别为: ...

〈

〉