0 引言

传统的地震波场数值模拟模式将地下介质视为完全弹性介质.这一假设忽略了地震波传播过程中的能量衰减等问题, 造成理论研究结果与实际采集到的地震波记录在振幅、相位等方面有一定差异, 因而简单地使用完全弹性介质模型不能有效解释许多实际地震波传播现象.实际上，地下介质是非完全弹性、各向异性和物性参数连续变化的非均匀介质.相比完全弹性介质, 使用黏弹性介质模型能更好地模拟地震波传播.

黏弹性介质是一种介于固相完全弹性介质和流相黏性介质之间的介质.相应的, 黏弹性方程兼具了波动方程和扩散方程的特点.从而, 不论在理论分析方面, 还是数值模拟方面, 黏弹性方程均比单纯的波动方程或者单纯的扩散方程更复杂.在对黏弹性介质研究过程中, 地球科学研究人员采用的模型有Kelvin介质模型、Maxwell介质模型、标准线性模型等.其中, Kelvin介质模型作为地震学中的一种重要模型, 一直受到广大地球科学研究人员的重视.孙成禹等(2010)基于Kelvin模型对黏弹性波动方程的稳定性进行研究, 得到了Q_p较大时的稳定性条件.贺同江等(2010)基于褶积微分算子进行了二维黏弹性介质地震波数值模拟.严红勇和刘洋(2012)基于数值模拟结果分析了时间域和频率域黏弹性介质对地震波的吸收和衰减规律.

模拟地震波在黏弹性介质中传播的数值方法有多种, 常用的有:几何射线法(Červeny and Firbas, 1984)、反射率法(Booth and Crampin, 1983)、有限元法(Chen, 1984;冯德山和王珣, 2017)、有限差分法(Kelly et al., 1976; Dablain, 1986; Virieux, 1984, 1986;董良国等, 2000a;辛维等, 2015;胡自多等, 2016)、谱元法或伪谱法(Komatitsch and Vilotte, 1998; Komatitsch and Tromp, 1999)等.不同的方法各有优劣, 在近些年的研究中许多学者(裴正林和牟永光, 2004;王雪秋等, 2005; Lu and Yang, 2005;李信富等, 2007)都进行了详细的分析.其中, 有限差分法因具有编程简单、计算速度快、存储量小等优势受到地球科学研究人员的青睐.

自1968年Alterman和Karal(1968)第一次采用有限差分方法进行地震波数值模拟到目前为止, 有限差分方法有了很大的发展.从网格划分角度讲，有限差分方法实现了规则网格—交错网格(Madariaga, 1976)—可变网格(Jastram and Behle, 1992)—不规则网格(Moczo, 1989;宋国杰等, 2010;韩颜颜等, 2015;豆辉和张剑锋, 2016)的发展.其中SG方法因在不增加计算量的前提下提高了模拟精度而受到学者们的广泛关注.1976年, Madariaga(1976)首次提出SG方法, 并完成了弹性介质地震波模拟, 精度为O(Δt²+h²).之后, Virieux (1984)提出了代替二阶位移方程的一阶速度-应力方程, 并采用SG方法进行数值模拟.Levander(1988)在Virieux(1984, 1986)的基础上提出了精度为O(Δt²+h⁴)的SG差分格式.随后, 董良国等(2000b)给出了一阶速度-应力方程更高精度的SG差分格式, 精度达到O(Δt⁴+h^2N).杨庆节等(2014)在董良国等(2000b)的基础上完善了时间4阶、空间2N阶精度的差分格式.但是, 计算精度的提高并不意味着SG方法能在较大空间步长条件下有效压制数值频散.因此, 寻找新的方法以压制较大空间步长条件下的数值频散具有重要的理论价值和现实意义.

2003年, Yang等(2003)提出了近似解析离散化方法(NADM), 该方法与传统方法的不同之处在于对空间高阶偏导数逼近的时候, 采用位移场及梯度场共同逼近, 其中梯度场反映波走势性质, 在计算过程中可以保留更多的波场信息, 从而使得该方法突破了传统差分方法的限制, 可以在较大空间步长条件下有效压制数值频散.之后在NADM基础上发展而来的ONADM(Yang et al., 2004, 2007a)、INADM(Yang et al., 2007b)、SNADM(宋国杰, 2008)、RNADM(Yang et al., 2010)、WNADM(宋国杰等, 2010)、NAIDM(Tong et al., 2011)、NSPRK(Ma et al., 2011)、MNAD(Chen et al., 2014)等方法均保留了NADM能在较大空间步长条件下压制数值频散的特点.汪勇等(2017)基于多种近似解析离散化方法——NADM、RNADM、INADM、WNADM研究了声波在黏弹性介质中的传播规律.值得注意的是NAD-type方法单个时间层内计算比较复杂, 尤其是三维情况下, 单个时间层计算时间较大.

综上所述, SG方法能够在单个时间层计算时间消耗较低的条件下实现高精度模拟, 但是在较大空间步长条件下无法有效压制数值频散; ONADM方法能保证较大空间步长条件下有效压制数值频散, 但单个时间层内计算时间消耗较大.本文集两种方法之长, 提出了一种能在较大空间步长条件下有效压制数值频散且时间消耗较少的黏滞声波波场模拟混合方法(SG-ONADM).在本文中, 我们首先给出了黏滞声波方程的SG方法、ONADM方法和SG-ONADM方法的离散格式, 分析了SG-ONADM方法的算法准确性; 然后在各向同性均匀介质模型中, 比较了SG、ONADM和SG-ONADM三种方法的计算效率; 最后使用SG-ONADM方法研究了黏滞声波在复杂介质中的传播.

1 方法原理 1.1 黏滞声波波动方程

黏滞声波波动方程三维分量形式为

(1)

其中, i=1、2、3, 且1、2、3分别表示x、y、z方向; u_i表示相应方向的位移分量, σ_ii表示相应方向的应力分量, f_i表示相应方向的体力; l_t是关于时间的一阶偏导数算子, 即分别为弹性介质和黏弹性介质的拉梅系数.根据地震勘探原理相关知识知, 其中ρ为介质密度, v_p为纵波速度, ω为角频率, Q_p为纵波品质因子.

1.2 交错网格方法

一阶黏滞声波波动方程的具体表达式为

其中, v_x、v_y、v_z分别是x、y、z方向的速度, P是应力, 在(2d)式中, 引入加速度场, 可以得到:

(3)

根据SG方法理论可以得到v_x、v_y、v_z、P的离散格式, 其差分模板见图 1a, 其中将应力P定义在整网格(i, j, k)处, 速度分量v_x和加速度分量定义在半网格处, 速度分量v_y和加速度分量定义在半网格处, 速度分量v_z和加速度分量定义在半网格处.对于时间格点, P定义在整时间网格点t和t－1处, 速度和加速度分量v_x、v_y、v_z、均定义在半网格点和处.为了方便书写, 我们用U, V, W分别表示v_x、v_y、v_z, 同时用U1, V1, W1表示.这里以U, U1为例, 给出具体的离散格式:

(4)

其中, C_n^(N)为差分系数(严红勇和刘洋，2012).

依次类推, 可写出V、V1、W、W1、P的表达式, 从而得到SG方法的完整离散格式.

1.3 最优近似解析离散化方法

为使用ONADM方法进行黏滞声波波场模拟, 需对公式(1)进行变换.考虑到Kelvin介质同均匀弹性各向同性介质一样, 体积应变系数(θ)和位移矢量(u)满足θ=∇·u的关系, 带入式(1)中, 整理得到的二阶黏滞声波波动方程为

(5)

为了避免同时对(5)式右侧的时间和空间混合偏导数进行计算, 我们令, 可以得到方程组：

其中, ∇²为拉普拉斯算子.

为叙述方便, 将位移场、速度场及其相应的梯度场合写成向量, 加速度场定义为.对于(6b)式, 如果直接简单地使用中心差分方法, 就会导致ONADM方法精度降低, 这并不是我们希望的.因此, 我们通过三阶泰勒截断展开来获得速度场的迭代表达式(6c), 即:

(6c)

在(6a)式中仅含有关于时间、空间的二阶偏导数, 按照ONADM方法在弹性介质中的离散方式处理即可(差分模板见图 1b).为实现位移场和速度场的时间推移, 首先要计算p、θ对时间的高阶偏导数.在此, 我们采用Dablain (1986)的思想将其转换为关于空间的偏导数, 然后采用Yang等(2007a)给出的高阶空间偏导数逼近公式对所涉及到的高阶偏导数进行计算.至此, 我们可以采用ONADM方法对式(6a)和式(6b)进行离散计算, 其中涉及到的高阶空间偏导数逼近表达式见附录.

1.4 SG-ONADM混合方法

在SG-ONADM方法中, 对于v_x、v_y、v_z随时间的推移采用SG方法进行计算, 具体参考1.2小节; P随时间的推移同时采用SG方法和ONADM方法.需要注意的是ONADM是一种基于整网格的有限差分方法, 为了保证高阶偏导数在计算时格点位置位于整网格点, 我们在对SG方法的格点分布定义过程中, 将P的空间、时间位置定义在整网格处, 具体做法如下:

根据式(2a)—(2c)给出的速度、应力关系, 将式(2d)中的l_t算子转换为应力P关于空间的一阶偏导数, 可以得到:

(7)

(7) 式中关于速度场的一阶空间偏导数我们仍采用SG方法进行差分.

接下来我们给出P的具体表达式为

(8)

(8) 式中关于P的高阶空间偏导数采用Yang等(2007a)一文给出的表达式，即：

(9)

其中 同x方向相类似.在(9)式中含有空间一阶偏导数, 这里我们采用SG方法的一阶空间偏导数离散格式进行计算, 从而实现了两类方法的混合.

1.5 算法准确性

众多学者(Madariaga, 1976;董良国等, 2000a, b ;严红勇和刘洋, 2012)的研究结果表明采用SG方法能够准确模拟地震波在地下介质中的传播过程.我们可以把精细网格下SG方法的模拟结果当作“准解析解”来判断本文所提出的SG-ONADM方法的准确性.因此, 我们设计了二维黏弹性介质模型, 模型大小为3 km×3 km, 纵波速度为v=3 km·s^-1, 品质因子Q_p=100, 密度为ρ=2.3 g·cm^-3.

数值模拟中, SG方法参数为: Δt=1 ms, h=Δx=Δz=5 m, SG-ONADM方法的模拟参数为: Δt=2 ms, h=Δx=Δz=25 m.震源函数为

(10)

其中震源主频f₀=30 Hz.

在一个波长里, SG方法的采样点数N_SG=v/(h_SGf₀)=3/(0.005×30)=20, 在这样的采样率条件下, SG方法的计算结果已经足够精确, 可以把它的计算结果视为“准解析解”, 而SG-ONADM方法的采样点数仅为N_SG-ONADM=v/(h_SG-ONADMf₀)=3/(0.025×30)=4.

图 2给出了分别采用SG和SG-ONADM方法模拟得到的时长0.4 s, 坐标位置为(1.8 km, 1.8 km)处的波形记录.从图 2可以看出, 采用SG-ONADM方法得到的波形与准解析解(每个波长内采样点数为20的SG网格方法)的波形基本一致.这说明SG-ONADM方法可以在较大空间步长条件下进行黏滞声波波场模拟, 并得到与采用精细网格下SG方法一致的波场模拟结果.

2 数值实验

为了说明本文所提出方法的有效性, 我们设计了以下2个模型.

2.1 模型1

我们选择均匀介质模型作为本文的第一个实验模型.模拟区域为3 km×3 km×3 km, 速度为v=4 km·s^-1, 介质密度ρ=2.5 g·cm^-3, 品质因子Q_p=100.模拟过程中, 分别采用ONADM方法, SG方法(时间2阶、空间4阶精度和10阶精度)和本文提出的SG-ONADM进行模拟, 并比较它们的模拟结果.计算参数均采用空间步长Δx=Δy=Δz=30 m, 时间步长Δt=3 ms, 震源位置(1.5 km, 1.5 km, 1.5 km), 震源函数如(10)式所示, 其中震源主频f₀=30 Hz.

图 3是使用三种方法在0.3 s时刻的模拟结果.由于本模型是一个均匀各向同性介质模型, 地震波传播过程中波场快照的xoy、xoz、yoz截面具有一致性, 因此我们仅选择xoz截面进行观察.直观地看到, SG-ONADM方法(图 3a)和ONADM方法(图 3b)均得到了清晰、无频散的波场快照.SG方法在4阶(图 3c)、10阶(图 3d)精度条件下均频散严重, 虽然10阶精度模拟结果较4阶精度模拟结果有了一定改善, 但是仍然无法呈现有效的波场快照.这说明提高模拟精度在一定程度上能够有效改善模拟结果, 但并不能使SG方法在较大空间步长条件下有效压制数值频散.

为了得到清晰的、无数值频散的波场快照, 现将4阶SG方法的计算参数进行适当的调整: Δx=Δy=Δz=5 m，Δt=1 ms, 得到的无频散波场快照如图 4所示.

上述实验, 从较大空间步长条件下的模拟效果角度对三种方法进行了比较, 说明SG-ONADM方法可以在较大空间步长条件下实现对数值频散的有效压制.接下来我们从计算效率的角度对三种方法进行比较.

理论上影响计算效率的主要因素包括内存需求和计算量两个方面.首先我们来分析三种方法对内存的需求.假设我们模拟M m×N m×L m的计算区域, SG方法在5 m的空间步长条件下, 需要存储11个大小的数组, ONADM在30 m的空间步长条件下, 需要存储24个大小为的数组, 而本文提出的SG-ONADM在30 m的空间步长条件下, 只需存储7个大小的数组.易知, SG方法的存储量是SG-ONADM方法的339.43倍, ONADM方法的存储量是SG-ONADM的3.43倍.计算量方面(我们忽略了单个加、减同乘、除的差异, 假定每一个基本运算的计算量是1, 幂的运算量是2), 在对M m×N m×L m计算区域进行等时间T的有效模拟条件下, SG方法的计算量是, ONADM的计算量是, SG-ONADM的计算量是, 即SG方法的计算量是SG-ONADM计算量的661.41倍, ONADM方法的计算量是SG-ONADM计算量的8.35倍.可以看出, SG-ONADM的计算量小于SG方法和ONADM方法.综合存储量和计算量两个方面的分析, 我们可以发现SG-ONADM的计算效率高于SG方法和ONADM方法.

为进一步说明SG-ONADM方法的计算效率较SG方法和ONADM方法高.基于模型1, 我们进行了等时间T的试算(SG方法采用调整之后的参数).试算结果表明, SG方法耗时36310.20 s; ONADM方法耗时305.60 s; SG-ONADM方法耗时26.32 s; 即SG-ONADM方法的计算效率是SG方法的1379.57倍, 是ONADM方法的11.61倍, 依然表明SG-ONADM的计算效率高于SG和ONADM方法.

注意到, SG和ONADM方法的实际计算效率要低于纯粹由计算量比较得到的理论计算效率, 这是因为SG和ONADM方法需要较高的内存来存储计算数组, 而内存存储量的增大又导致内存Cache命中率降低, 从而减慢了计算速度.

图 5给出了SG-ONADM方法结合CPML边界对模型1进行模拟得到的0.45 s时刻波场快照, 此时声波传播到模拟区域边界.可以观察到, 在边界处无任何可见的边界反射, 说明SG-ONADM方法结合CPML边界很好地完成了对边界虚假反射的吸收.

2.2 模型2

为了验证SG-ONADM方法在复杂介质的有效性, 我们选择了局部盐丘模型(完整盐丘模型大小为676×676×210, 本文所选区域大小为400×400×200)进行数值模拟实验.盐丘模型(图 6)中速度v的变化范围为1.5~4.482 km·s^-1, 密度ρ变化范围为2.08~2.58 g·cm^-3.震源放置在(4 km, 4 km, 0.2 km)处, 即图中☆所示.空间步长Δx=Δy=Δz=20 m, 时间步长Δt=2 ms, 震源主频、品质因子分别为f₀=15 Hz和Q_p=150.震源函数见(10)式, 边界采用CPML条件.

图 7分别是4 s时刻z=0 km, y=4 km和z=0 km, x=4 km处地震记录.其中存在直达波、反射波以及大量的多次波, 无可见数值频散.说明本文提出的SG-ONADM应用于复杂介质中的地震波场模拟是有效的, 同时也说明本方法和CPML可以很好地结合在一起进行复杂介质黏滞声波波场模拟.

3 结论

研究地震波在黏弹性介质中的传播规律具有非常重要的理论指导意义和应用价值.本文在SG方法和ONADM方法的基础上提出了一种新的黏声波波场数值模拟的混合方法——SG-ONADM.该方法中对控制方程中的一阶空间偏导数采用SG方法进行近似; 对高阶空间偏导数采用ONADM的高阶偏导数逼近公式进行计算.理论分析与数值算例均表明SG-ONADM方法计算效率高于SG方法和ONADM方法, 能在较大空间步长条件下有效压制数值频散, 可以用来进行大规模黏滞声波波场模拟.

下一步的研究中, 我们拟将该方法推广至黏弹性地震波场数值模拟.

附录 ONADM方法计算过程中涉及到的高阶偏导数逼近公式

为了叙述方便, 我们定义了单位算子符I, 二阶中心差分算子δ²和位移运算符E, 具体表达式为

(A1)

ONADM方法计算过程中涉及到的高阶空间导数逼近表达式为

(A2)

(A3)

(A4)

(A5)

(A6)

(A7)

(A8)

(A9)

(A10)

(A4)中, g, e=x, y, z且g≠e.定义；(A6)中的三种情况分别为g=x, e=y、g=x, e=z和g=y, e=z；在(A8)和(A9)中, g, e=x, y, z, g≠e, 定义

；(A10)对应三种情况分别为: g=x, e=y, d=z或g=y, e=z, d=x或g=x, e=z, d=y.定义