三维波动方程时空域混合网格有限差分数值模拟方法

2021-08-03 11:11胡自多刘威雍学善王小卫韩令贺田彦灿

地球物理学报 2021年8期

胡自多，刘威，雍学善，王小卫，韩令贺，田彦灿

1 中国石油勘探开发研究院西北分院，兰州 730020 2 中国石油天然气集团有限公司油藏描述重点实验室，兰州 730020

0 引言

波动方程数值模拟是逆时偏移(Baysal et al., 1983; Virieux et al., 2011; 陈生昌和周华敏, 2018)和全波形反演(Tarantola, 1984; Pratt et al., 1998; 董良国等, 2015)的重要基础.有限元法(Marfurt, 1984; Moczo et al., 2010, 2011)、虚谱法(Reshef et al., 1988; 黄建平等, 2016)和有限差分法(Alterman and Karal, 1968; Liu and Sen, 2011a; 梁文全等, 2013)是求解波动方程的三大主流数值算法(皮红梅等，2009；胡恒山，2018).相比其他两类方法，有限差分法占用内存更小，计算效率更高，因而成为应用最普遍的波动方程数值模拟算法(蒋韬等，2008；张志禹等，2017；姜占东等，2021).但是，有限差分法采用差分算子近似波动方程中的时间和空间偏微分算子，导致时间和空间数值频散，降低了模拟精度(Alford et al., 1974; Dablain, 1986)；同时，逆时偏移和全波形反演巨大的计算量，发展高精度和高效率的有限差分模拟算法具有重要意义.

Dablain(1986)研究指出，时间和空间高阶差分格式能够有效压制数值频散、提高模拟精度.然而，时间高阶差分格式占用内存过大，稳定性降低(Chen, 2007, 2011).针对标量波动方程，常规高阶差分方法采用时间2阶和空间高阶(2M)差分格式，并利用空间域频散关系和泰勒级数展开求解差分系数，本文称为常规空间域高阶有限差分法(简称为CSD-FDM).CSD-FDM的差分系数算法仅衡量了空间差分算子(Laplace差分算子)的精度，而没有考虑差分离散波动方程的精度，空间差分算子虽然能达到2M阶精度，但差分离散波动方程仅具有2阶精度.实际上，有限差分法通过迭代求解差分离散波动方程实现波动方程数值模拟，差分离散波动方程的差分精度才能更准确地描述有限差分法的模拟精度.Liu和Sen (2009)保持CSD-FDM的差分格式不变，提出利用时空域频散关系改进差分系数算法，本文称为时空域高阶有限差分法(简称为TSD-FDM).TSD-FDM比CSD-FDM具有更高的模拟精度和更强的稳定性.TSD-FDM在二维和三维标量波动方程模拟中，差分离散波动方程分别沿8个和48个传播方向具有2M阶差分精度，其他传播方向只能达到2阶差分精度，存在数值各向异性.基于相同思路，Liu和Sen(2011b)提出了面向应力-速度声波方程的时空域高阶交错网格差分法，同样有效提高了模拟精度和稳定性.严红勇等进一步推广应用于声波、VTI介质、黏声介质逆时偏移中，有效压制了数值频散造成的成像假象(严红勇和刘洋, 2013；Yan and Liu，2013a，b).基于泰勒级数展开的差分系数算法计算效率高，通常低波数成分具有较高的模拟精度，但高波数模拟精度迅速降低.为了提高高波数成分的模拟精度，Zhang和Yao (2012,2013)、Liu(2013,2014)通过最小二乘优化频散误差，改进差分系数算法提高模拟精度.

除了改进差分系数算法，设计更合理的差分格式是提高模拟精度的另一重要途径.针对二维和三维应力-速度声波方程，Tan和Huang(2014a)利用坐标轴网格点和非坐标轴网格点构建空间差分算子近似一阶空间微分算子，提出了混合交错网格差分方法，并利用时空域频散关系和泰勒展开求解差分系数，差分离散波动方程沿任意传播方向可达到4阶、6阶差分精度，并具有更强的稳定性.Tan和Huang(2014b)利用最小二乘法优化差分系数，模拟精度进一步提高.Ren和Li(2017)针对二维弹性波方程，构建了更具一般性的混合交错网格差分格式，弹性波差分离散方程沿任意传播方向可达到2N(N<5)阶差分精度.混合交错网格差分格式，应力和速度的空间一阶偏导算子，各坐标轴相互独立计算，因此2D可较易推广到3D.

针对二维标量波动方程，Liu和Sen(2013)、张保庆等(2016)、胡自多等(2016)和Wang等(2016)先后提出将二维Laplace微分算子近似为常规笛卡尔坐标系中坐标轴网格点构建的M个Laplace差分算子和旋转笛卡尔坐标系中坐标轴网格点构建的N个Laplace差分算子的加权平均，构建了不同的混合网格差分格式，利用时空域频散关系和泰勒展开求解差分系数，差分离散波动方程沿任意传播方向理论上可达到任意偶数阶差分精度，明显提高了二维标量波动方程有限差分法的模拟精度，稳定性进一步增强，还能采用更大的时间采样间隔以获得更高的计算效率.不同于混合交错网格差分格式，标量波动方程压力对空间的2阶偏导(Laplace算子)作为一个整体计算，2D混合网格差分格式不能直接推广到3D.而地震勘探的对象通常是三维介质，发展三维混合网格差分格式更具实际意义.

如何利用三维笛卡尔坐标系中非坐标轴网格点构建Laplace差分算子是建立三维混合网格差分格式需要解决的关键问题.本文将三维笛卡尔坐标系中的非坐标轴网格点分为两类：坐标平面内的非坐标轴网格点和坐标平面外的非坐标轴网格点，并系统推导出了这两类非坐标轴网格点构建三维Laplace差分算子的方法.在此基础上，将三维Laplace微分算子近似为坐标轴网格点构建的M个Laplace差分算子和非坐标轴网格点构建的N个Laplace差分算子的加权平均，提出了一种面向三维标量波动方程的混合网格M+N型差分方法(简称为MG-FDM)，并根据时空域频散关系和泰勒级数展开建立了差分系数求解方程组，推导出差分系数的通解.增大N的取值，三维MG-FDM给出的差分离散波动方程沿任意传播方向能够达到4阶、6阶、甚至任意偶数阶差分精度.频散分析表明，相比三维CSD-FDM和TSD-FDM，计算效率基本相同时，MG-FDM的数值频散压制效果更好，模拟精度更高；模拟精度基本相当时，MG-FDM能采用更大的时间采样间隔，计算效率更高.稳定性分析表明，MG-FDM稳定性更强，为其采用更大的时间采样间隔以提高计算效率奠定了基础.三维数值模拟实验进一步证实了MG-FDM在提高模拟精度和计算效率方面的优越性.

1 三维混合网格差分格式构建方法

1.1 三维常规空间域高阶(CSD-FDM)和时空域高阶(TSD-FDM)差分格式

常密度介质中，三维标量波动方程为：

(1)

有限差分法数值求解波动方程普遍采用2阶时间差分方案，式(1)右边的时间偏微分可差分表示为：

(2)

三维常规空间域高阶差分(CSD-FDM)和三维时空域高阶差分(TSD-FDM)，均采用图1所示的差分格式，波动方程(1)中的三维Laplace微分算子可差分近似为：

(3)

图1 三维CSD-FDM和TSD-FDM的差分格式Fig.1 FD stencil of 3D CSD-FDM and TSD-FDM

式(3)、(2)代入式(1)：

(4)

1.2 三维混合网格M+N型差分格式构建

构建三维混合网格M+N型差分格式(MG-FDM)的基本思想是联合利用坐标轴网格点和非坐标轴网格点一起差分近似三维Laplace算子.如何利用非坐标轴网格点构建Laplace差分算子是建立三维MG-FDM需要解决的关键问题.

图2a—c给出了三维笛卡尔坐标系中与差分中心点距离相等的三组非坐标轴网格点，这三组网格点与差分中心点的距离依次增大.由于无法将任意一组非坐标轴网格点置于坐标原点位于差分中心点的三维旋转笛卡尔坐标系的坐标轴上，导致不能直接利用非坐标轴网格点构建三维Laplace差分算子.为了利用非坐标轴网格点构建三维Laplace差分算子，本文将非坐标轴网格点分成两类：坐标平面内的非坐标轴网格点(如图2a、c)和坐标平面外的非坐标轴网格点(如图2b).

(5)

(6)

将三维坐标系中坐标轴网格点差分格式(图1)与非坐标轴网格点差分格式(图2)进行组合，构建出三维MG-FDM的差分格式(图3).图1与图2a中的差分格式组合，可构建出三维MG-FDM(N=1)的差分格式(图3a).式(3)与式(5)加权求和：

(7)

其中cm(m=1,2,…,M)和c1,1,0为权系数.

图2 三维坐标系中非坐标轴网格点构建的差分格式(a) 坐标平面内的非坐标轴网格点构建的差分格式(网格点与中心点的距离是坐标平面外的非坐标轴网格点构建的差分格式(网格点与中心点的距离是坐标平面内的非坐标轴网格点构建的差分格式(网格点与中心点的距离是Fig.2 FD stencils constructed by off-axial grid nodes in the 3D Cartesian coordinate system(a) FD stencil constructed by the off-axial grid nodes in the coordinate plane (distance between grid node and center grid node is (b) FD stencil constructed by the off-axial grid nodes outside the coordinate plane (distance between grid node and center grid node is (c) FD stencil constructed by the off-axial grid nodes in the coordinate plane (distance between grid node and center grid node is

图3 三维MG-FDM的差分格式(a) MG-FDM(N=1)； (b) MG-FDM(N=2)； (c) MG-FDM(N=3).Fig.3 FD stencils of 3D MG-FDM

式(7)表明，三维MG-FDM(N=1)是将三维Laplace微分算子近似为坐标轴网格点构建的M个Laplace差分算子和非坐标轴网格点构建的1个Laplace差分算子的加权平均.式(7)和(2)代入(1)：

(8)

式(8)是三维MG-FDM(N=1)的差分离散波动方程，附录B给出了三维MG-FDM(N=2,3)对应的差分离散波动方程，同理可导出N取任意值时三维MG-FDM的差分离散波动方程.MG-FDM实质上是将三维Laplace微分算子近似为坐标轴网格点构建的M个Laplace差分算子和非坐标轴网格点构建的N个Laplace差分算子的加权平均.

2 差分系数计算方法和差分精度分析

三维CSD-FDM和TSD-FDM的差分系数计算，分别利用空间域频散关系和时空域频散关系.为了阐述三维MG-FDM差分系数计算方法和差分精度对比，我们首先列出三维CSD-FDM和TSD-FDM的差分系数算法，进而推导出三维TSD-FDM和MG-FDM差分系数的通解.

2.1 三维CSD-FDM的差分系数算法

在均匀介质中，三维标量波动方程具有如下离散形式的平面波解：

(9)

其中kx=ksinφcosθ,ky=ksinφsinθ,kz=kcosφ，ω为角频率，k为波数，φ和θ为平面波的传播方向角.

将式(9)代入式(3)中，得到CSD-FDM 的Laplace差分算子频散关系，即空间域频散关系：

+2cos(mkzh)-6].

(10)

对余弦函数进行泰勒展开，可导出三维CSD-FDM的差分系数通解为：

(m=1,2,…,M).

(11)

三维CSD-FDM差分系数计算，利用空间域频散关系，仅考虑空间差分算子(Laplace差分算子)的差分精度.

2.2 三维TSD-FDM的差分系数算法

三维TSD-FDM差分系数计算，考虑了差分离散波动方程的差分精度.将平面波解式(9)代入差分离散方程式(4)，得到差分离散波动方程的频散关系，即时空域频散关系：

(12)

(13)

其中r=υτ/h为Courant条件数，表示单位时间采样间隔(τ)内波传播的距离与空间采样间隔(h)之比.

(m=1,2,…,M)，

(14)

其中f(j,φ,θ)=sin2jφ(cos2jθ+sin2jθ)+cos2jφ.可以看出，三维TSD-FDM的差分系数am(m=0,1,2,…,M)是地震波传播方向角φ和θ的函数.Liu和Sen(2009)计算差分系数时，取φ=π/2,θ=π/8，分析出三维TSD-FDM的差分离散波动方程沿48个传播方向具有2M阶差分精度.

取φ=0,θ=0，式(14)可简化得到TSD-FDM差分系数通解的一个特殊形式：

(m=1,2,…,M)，

(15)

取φ=0,θ=0，三维TSD-FDM的差分离散波动方程仅沿6个传播方向(3个坐标轴的正负方向)具有2M阶差分精度.因此，取φ=π/2,θ=π/8计算的差分系数，TSD-FDM的差分离散波动方程能达到更高的模拟精度.

2.3 三维MG-FDM的差分系数算法

借鉴时空域频散关系和泰勒展开求解差分系数的方法，本文推导出了MG-FDM的差分系数通解.

将平面波解式(9)代入三维MG-FDM(N=1)的差分离散波动方程式(8)，得到：

+cos(kxh+kyh)+cos(kyh-kzh)+cos(kyh+kzh)+cos(kzh-kxh)+cos(kzh+kxh)-6]，

(16)

式(16)描述了三维MG-FDM(N=1)的差分离散波动方程的时空域频散关系，泰勒展开余弦函数可得：

+(ky-kz)2j+(ky+kz)2j+(kz-kx)2j+(kz+kx)2j]h2j-2.

(17)

(j=1,2,…,M)，

(18)

(19)

2.4 差分精度分析和对比

差分精度是一种定量描述有限差分法模拟精度的传统方法，通常将误差函数关于时间采样间隔τ或者空间采样间隔h的最小幂指数称为差分精度阶数.根据频散关系式(10)，定义三维CSD-FDM的Laplace差分算子的误差函数εCSD-FDM为：

(20)

根据频散关系式(12)，定义三维CSD-FDM的差分离散波动方程的误差函数ECSD-FDM为：

(21)

根据三维TSD-FDM差分离散波动方程的频散关系式(12)，结合式(13)，TSD-FDM的差分离散波动方程的误差函数：

(22)

计算差分系数时，取φ=π/2,θ=π/8，三维TSD-FDM的差分离散波动方程，沿特定的48个传播方向可达到2M阶差分精度，但沿其他传播方向仅具有2阶差分精度.整体来讲，三维TSD-FDM的差分离散波动方程仅具有2阶差分精度.

根据三维MG-FDM(N=1)的差分离散波动方程的频散关系式(16)，定义MG-FDM(N=1)的差分离散波动方程的误差函数EMG-FDM(N=1)为：

(23)

结合式(17)、(18)，误差函数EMG-FDM(N=1)可表示为：

(24)

同理，可分析出三维MG-FDM(N=2,3)的差分精度.表1给出了三维CSD-FDM、TSD-FDM和MG-FDM(N=1,2,3)的差分精度，可以看出，相比三维CSD-FDM和TSD-FDM，MG-FDM能够有效提高差分离散波动方程的差分精度.增大N的取值，MG-FDM的差分离散波动方程理论上可达到任意偶数阶差分精度.

表1 三维CSD-FDM、TSD-FDM和MG-FDM(N=1,2,3)的差分离散波动方程的差分精度Table 1 Difference accuracy of the discrete difference wave equation for 3D CSD-FDM, TSD-FDM and MG-FDM (N=1,2,3)

3 数值频散分析和稳定性分析

3.1 数值频散分析

数值频散使得相速度υph与地震波的真实传播速度不相等，并且随地震波的传播方向变化，呈现数值各向异性特征.本文采用归一化相速度δ(φ,θ)=υph/υ描述相速度的数值频散特性，根据相速度定义υph=ω/k和三维MG-FDM(N=1)的差分离散波动方程的频散关系式(16)，可得出MG-FDM(N=1)的归一化相速度δ(φ,θ)的表达式：

(25)

(26)

其中，G=λ/h，λ为波长，G为波长与空间采样间隔之比，则G值越大，1/G越小，表示空间采样越密，反之，空间采样越稀疏.

归一化相速度δ(φ,θ)的值越接近1，相速度数值频散越小；δ(φ,θ)>1，相速度偏大，称为时间数值频散；δ(φ,θ)<1，相速度偏小，称为空间数值频散.另外,δ(φ,θ)随地震波传播方向角φ和θ变化越小，数值各向异性越弱，反之，数值各向异性越强.相速度数值频散越小，数值各向异性越弱，则数值模拟精度越高；相反，相速度数值频散越大，数值各向异性越强，则数值模拟精度越低.

同样地，可导出三维CSD-FDM、TSD-FDM和MG-FDM(N=2,3)的归一化相速度δ(φ,θ)的表达式.δ(φ,θ)是传播方向角φ和θ的函数，分析相速度频散特性时给出了8个传播角度(φ=π/2;θ=0,π/16,2π/16,3π/16,4π/16)和(φ=π/12,2π/12,3π/12;θ=2π/16)的相速度频散曲线.

图4给出了三维CSD-FDM(M=2,6,10)、TSD-FDM(M=2,6,10)和MG-FDM(M=2,6,10;N=1)的相速度频散曲线，Courant条件数r的取值均为r=0.3.分析对比相速度频散曲线特征，可以得出：

(1)三维CSD-FDM(M=2)、TSD-FDM(M=2)和MG-FDM(M=2;N=1)的相速度频散特征相似，均具有明显的空间数值频散，模拟精度都很低.

(2)三维CSD-FDM(M=6,10)具有明显的时间数值频散；TSD-FDM(M=6,10)的相速度频散曲线较发散，同时存在空间和时间数值频散，表现出明显的数值各向异性(数值频散特征随地震波的传播方向变化)；MG-FDM(M=6,10;N=1)的相速度频散曲线收敛，数值各向异性明显减弱，数值频散幅值也明显减小.

图4 三维CSD-FDM、TSD-FDM和MG-FDM(N=1)的相速度频散曲线(r=0.3)(a)(b)(c) CSD-FDM(M=2,6,10); (d)(e)(f) TSD-FDM(M=2,6,10); (g)(h)(i) MG-FDM(M=2,6,10;N=1)Fig.4 Phase velocity dispersion curves of 3D CSD-FDM, TSD-FDM and MG-FDM (N=1) with r=0.3

综合(1)和(2)可以看出：相比三维CSD-FDM和TSD-FDM，M取值较小时(M=2左右)，MG-FDM在压制相速度数值频散方面无明显优势；M取值较大时(M≥6)，MG-FDM的相速度数值频散和数值各向异性均明显减小，因而具有更高的模拟精度.

图5给出了三维MG-FDM(M=8,15;N=1,2,3)的相速度频散曲线(r=0.3).需要注意，图5a、b、c和图5d、e、f两组频散曲线具有不同的纵轴刻度.

图5 三维MG-FDM的相速度频散曲线(r=0.3)(a)(b)(c) MG-FDM(M=8;N=1,2,3)； (d)(e)(f) MG-FDM(M=15;N=1,2,3).Fig.5 Phase velocity dispersion curves of 3D MG-FDM with r=0.3

(1)三维MG-FDM(M=8;N=1,2,3)的相速度数值频散特征基本相同；相比MG-FDM(M=15;N=1,2)，MG-FDM(M=15;N=3)的相速度频散曲线更收敛，数值频散幅值和数值各向异性均明显减小.因此，三维标量波动方程数值模拟对精度要求较高时，建议采用MG-FDM(N=1)，且M取值为8左右；对模拟精度要求极其苛刻时，可采用MG-FDM(N=3)，且M取值为15左右；N取值大于4，同时M取值更大的三维MG-FDM因为数值计算量巨大，计算效率极低而很少采用.三维MG-FDM可根据模拟精度要求选择适当的M和N值.

(2)三维MG-FDM(M=8;N=1,2)的相速度数值频散特征基本相同，同样地，MG-FDM(M=15;N=1,2)的相速度数值频散特征也基本相同，和三维MG-FDM(N=1,2)的差分离散波动方程均具有4阶差分精度的结论一致(见表1).

图6给出了r取值分别为r=0.1,0.2,0.4时，三维CSD-FDM(M=10)，TSD-FDM(M=10)和MG-FDM(M=8;N=1)的相速度频散曲线.根据r=υτ/h，空间采样间隔h保持不变时，r取值从0.1增大至0.2和0.4，相应的时间采样间隔τ将分别增大至2倍和4倍.

(1)r取值为0.1，三维CSD-FDM表现出轻微的时间频散，具有较高的模拟精度；r取值增大至0.2和0.4，时间频散显著增强，模拟精度低.

(2)r取值为0.1和0.2时，三维TSD-FDM的相速度频散曲线较收敛，表现出轻微的空间和时间频散，具有较高的模拟精度；r取值增大至0.4时，频散曲线严重发散，数值频散的幅值明显增大，数值各向异性显著增强，模拟精度低.

(3)当r取值为0.1、0.2和0.4时，三维MG-FDM的相速度频散曲线收敛性均较好，数值频散微弱，模拟精度高.

图6 三维CSD-FDM(M=10)，TSD-FDM(M=10)和MG-FDM(M=8;N=1)的相速度频散曲线(a)(b)(c) 分别为r=0.1,0.2,0.4时CSD-FDM(M=10)的相速度频散曲线； (d)(e)(f) 分别为r=0.1,0.2,0.4时TSD-FDM(M=10)的相速度频散曲线； (g)(h)(i) 分别为r=0.1,0.2,0.4时MG-FDM(M=8;N=1)的相速度频散曲线.Fig.6 Phase velocity dispersion curves of 3D CSD-FDM (M=10), TSD-FDM (M=10) and MG-FDM (M=8; N=1)(a),(b) and (c) CSD-FDM (M=10) with r=0.1,0.2,0.4; (d),(e) and (f) TSD-FDM (M=10) with r=0.1,0.2,0.4; (g),(h) and (i) MG-FDM (M=8;N=1) with r=0.1,0.2,0.4.

综合可知，三维MG-FDM(M=8;N=1)可以比CSD-FDM(M=10)和TSD-FDM(M=10)取更大的r值(速度模型υ和空间采样间隔h相同时，r取值越大等价于时间采样间隔τ越大)，同时数值频散大小基本相当.因此，MG-FDM(M=8;N=1)能采用较大的时间采样间隔以提高计算效率，同时保持较高的模拟精度.

3.2 稳定性分析

有限差分法通过迭代求解差分离散波动方程模拟地震波的传播过程，必须确保迭代过程稳定.根据式(16)可得出：

-cos(kxh+kyh)-cos(kyh-kzh)-cos(kyh+kzh)-cos(kzh-kxh)-cos(kzh+kxh)]，

(27)

取空间波数kx=ky=kz=π/h，并且0≤1-cos(ωτ)≤2，则可得到：

(28)

图7给出了三维CSD-FDM、TSD-FDM和MG-FDM(N=1,2,3)的稳定性曲线，描述了稳定性条件约束下的最大r取值随M的变化关系.稳定性条件约束前提下，r取值越小，稳定性越弱；相反，r取值越大，稳定性越强.分析图中的稳定性曲线可以看出：

(1)M的取值增大，三维CSD-FDM、TSD-FDM和MG-FDM(N=1,2,3)的稳定性均下降.

(2)N的取值增大，MG-FDM的稳定性增强，但差分离散波动方程的差分精度相同时，稳定性基本相同，MG-FDM(N=1,2)的差分离散波动方程均为4阶差分精度，稳定性基本一致.

(3)三维CSD-FDM的稳定性最弱，TSD-FDM的稳定性明显增强，MG-FDM(N=1,2,3)的稳定性进一步增强.MG-FDM(N=1,2,3)的强稳定性为数值模拟采用更大的时间采样进而提高计算效率奠定了基础.

图7 三维CSD-FDM，TSD-FDM和MG-FDM(N=1,2,3)的稳定性曲线Fig.7 Stability curves of 3D CSD-FDM, TSD-FDM and MG-FDM (N=1,2,3)

4 计算效率分析

计算效率是衡量有限差分格式优劣的另一项重要指标.图6中的相速度频散曲线表明，三维CSD-FDM和TSD-FDM均可采用较小的时间采样间隔来减小数值频散，以获得较高的模拟精度，但计算量会显著增加.MG-FDM的相速度频散曲线收敛性好，模拟精度高，稳定性强，可采用更大的时间采样间隔以提高计算效率，同时保持较高的模拟精度.

下面以三维CSD-FDM(M=10)、TSD-FDM(M=10)和MG-FDM(M=8;N=1)为例分析三种差分格式具有基本相同的模拟精度条件下的计算效率.这三种差分格式均包含61个网格点，单次时间迭代的浮点运算量基本相同，模拟时间长度相等时，计算效率之比约等于时间采样间隔之比.

分析图6中的相速度频散曲线，可以认为CSD-FDM(M=10)取r=0.1，TSD-FDM(M=10)取r=0.2和MG-FDM(M=8;N=1)取r=0.4基本能保持同样高的模拟精度.根据r=vτ/h,在速度v和空间采样间隔h保持不变的条件下，时间采样间隔τ之比等价于r之比.此分析表明，模拟精度基本相当的前提下，TSD-FDM可采用2倍于CSD-FDM的时间采样间隔，MG-FDM可采用4倍于CSD-FDM的时间采样间隔.因此，TSD-FDM的计算效率是CSD-FDM的2倍；MG-FDM的计算效率是CSD-FDM的4倍，是TSD-FDM的2倍.表2列出了模拟精度基本相当时，三种差分格式采用的时间采样间隔和理论加速比.

表2 三维CSD-FDM、TSD-FDM和MG-FDM三种差分格式的计算效率分析Table 2 Analysis of computational efficiency for 3D CSD-FDM, TSD-FDM and MG-FDM

5 数值模拟实例

5.1 层状介质模型

图8a给出了一个8 km×3 km×6 km的三层模型：v1=2000 m·s-1,h1=3.0 km;v2=2450 m·s-1,h2=1.2 km;v3=2680 m·s-1,h3=1.8 km.空间采样间隔Δx=Δy=Δz=h=10 m，模型网格数为nx×ny×nz=801×301×601，主频25 Hz的Ricker子波作为震源，位于网格点(51,21,3).三维CSD-FDM(M=10)、TSD-FDM(M=10)和MG-FDM(M=8;N=1)分别采用不同的时间采样间隔进行数值模拟.图8b给出了MG-FDM(M=8;N=1)采用时间采样间隔τ=1.5 ms进行数值模拟得到的炮集，为了便于对比分析，图9给出了三维CSD-FDM(M=10)采用时间采样间隔τ=0.5 ms(图9a)和τ=1.0 ms(图9b)，TSD-FDM(M=10)采用τ=1.0 ms(图9c)和τ=1.5 ms(图9d)，MG-FDM(M=8;N=1)采用τ=1.0 ms(图9e)和τ=1.5 ms(图9f)进行数值模拟得到的局部炮集，局部区域范围由坐标标出.

图8 三维层状速度模型及数值模拟炮集(a) 层状速度模型； (b) MG-FDM(M=8;N=1)数值模拟炮集，时间采样间隔τ=1.5 ms.Fig.8 3D layered velocity model and numerical modeling seismogram(a) 3D layered velocity model； (b) Numerical modeling seismogram using MG-FDM(M=8;N=1) with τ=1.5 ms.

图9 三维层状介质模型数值模拟的局部炮集(a)(b) CSD-FDM(M=10), 时间采样间隔τ=0.5 ms和τ=1.0 ms; (c)(d) TSD-FDM(M=10), 时间采样间隔τ=1.0 ms和τ=1.5 ms; (e)(f) MG-FDM(M=8;N=1), 时间采样间隔τ=1.0 ms和τ=1.5 ms.Fig.9 Local parts of the numerical modeling seismograms on 3D layered model(a) and (b) CSD-FDM(M=10) with τ=0.5 ms and τ=1.0 ms; (c) and (d) TSD-FDM(M=10) with τ=1.0 ms and τ=1.5 ms; (e) and (f) MG-FDM(M=8;N=1) with τ=1.0 ms and τ=1.5 ms.

图10给出了三维CSD-FDM(M=10)采用τ=0.5 ms,0.75 ms,1.0 ms，TSD-FDM(M=10)采用τ=0.5 ms,1.0 ms,1.25 ms,1.5 ms和MG-FDM(M=8;N=1)采用τ=1.0 ms,1.25 ms,1.5 ms进行数值模拟得到的单道波形图，检波器位于网格点(650,3,3).

图10 三维层状介质模型数值模拟单道波形对比图检波器位于(650,3,3)，蓝色为参考波形,CSD-FDM(M=10)采用极小时间采样间隔τ=0.1 ms得到；红色为不同差分格式模拟波形. ①②③分别为CSD-FDM(M=10), 时间采样间隔τ=0.5 ms,0.75 ms,1.0 ms；④⑤⑥⑦分别为TSD-FDM(M=10)，时间采样间隔τ=0.5 ms,1.0 ms,1.25 ms,1.5 ms; ⑧⑨⑩分别为M2M+N-FD(M=8;N=1), 时间采样间隔τ=1.0 ms,1.25 ms,1.5 ms.Fig.10 Comparison of numerical modeling single trace waveforms on 3D layered model The geophone is located at (650,3,3), and the blue line is the reference waveform which is obtained by CSD-FDM (M=10) with a very small time sampling interval τ=0.1 ms; The red lines are the modeling waveform with different FD schemes. ①②③CSD-FDM(M=10) with τ=0.5 ms,0.75 ms,1.0 ms；④⑤⑥⑦TSD-FDM(M=10) with τ=0.5 ms,1.0 ms,1.25 ms,1.5 ms； ⑧⑨⑩ M2M+N-FD(M=8; N=1) with τ=1.0 ms,1.25 ms,1.5 ms.

从图9和图10可以看出：CSD-FDM采用时间采样间隔τ=0.5 ms，未出现明显数值频散，采用τ=0.75 ms,1.0 ms，出现明显的时间频散；TSD-FDM采用τ=0.5 ms,1.0 ms，未出现明显的数值频散，采用τ=1.25 ms,1.5 ms，出现较明显的数值频散；MG-FDM采用τ=1.0 ms和τ=1.5 ms，均未出现明显的数值频散.表3给出了层状介质模型单炮模拟时，三种差分格式采用不同时间采样间隔的计算机耗时和加速比，数值模拟均采用Inter Xeon CPU E5-2670处理器.

层状模型数值模拟结果表明：确保不出现明显数值频散，实现高精度数值模拟，MG-FDM(M=8;N=1)可达到CSD-FDM(M=10)的计算效率的3.31倍，可达到TSD-FDM(M=10)的计算效率的1.5倍.

表3 三维层状介质模型CSD-FDM、TSD-FDM和MG-FDM数值模拟计算效率对比Table 3 Comparison of computational efficiency for numerical modeling on 3D layered model with CSD-FDM, TSD-FDM and MG-FDM

5.2 盐丘模型

图11b给出了MG-FDM(M=8;N=1)采用时间采样间隔τ=1.5 ms进行数值模拟得到的炮集，同样为了分析对比方便，图12仅给出三种差分格式数值模拟的两个局部炮集.图12a、b给出了CSD-FDM采用时间采样间隔τ=1.0 ms，图12c、d给出TSD-FDM采用τ=1.25 ms，图12e、f给出MG-FDM采用τ=1.5 ms进行数值模拟得到的两个局部炮集.图13给出了CSD-FDM采用τ=0.5 ms,1.0 ms，TSD-FDM(M=10)采用τ=1.0 ms,1.25 ms和MG-FDM采用τ=1.25 ms,1.5 ms的两个单道波形记录，两个检波器分别位于网格点(3,650,3)和(560,650,3).

图11 修改的三维SEG/EAGE盐丘速度模型及数值模拟炮集(a) 三维盐丘速度模型； (b) MG-FDM(M=8;N=1)数值模拟炮集，时间采样间隔τ=1.5 ms，顶部显示了3.0 s时刻的切片.Fig.11 Revised 3D SEG/EAGE salt model and numerical modeling seismogram(a) Revised 3D SEG/EAGE salt model；(b) Numerical modeling seismogram using MG-FDM(M=8;N=1) with τ=1.5 ms. The top shows the slice at 3.0 s.

图12 修改的三维SEG/EAGE盐丘模型数值模拟的两个局部炮集(a)(b) CSD-FDM(M=10), 时间采样间隔τ=1.0 ms; (c)(d) TSD-FDM(M=10), 时间采样间隔τ=1.25 ms; (e)(f) MG-FDM(M=8;N=1), 时间采样间隔τ=1.5 ms.Fig.12 Two local parts of the numerical modeling seismogram on the revised 3D SEG/EAGE salt model(a) and (b) CSD-FDM(M=10) with τ=1.0 ms； (c) and (d) TSD-FDM(M=10) with τ=1.25 ms； (e) and (f) MG-FDM(M=8;N=1) with τ=1.5 ms.

从图12和图13可以看出：CSD-FDM采用τ=1.0 ms，表现出明显的时间数值频散；TSD-FDM采用τ=1.25 ms，表现出较明显的数值频散；MG-FDM采用τ=1.5 ms，未出现明显的数值频散.表4给出了盐丘模型单炮模拟时，三种差分格式采用不同时间采样间隔的计算机耗时和加速比.

图13 修改的三维SEG/EAGE盐丘模型数值模拟单道波形对比图(a)(b)检波器分别位于(3,650,3)和(560,650,3)，蓝色为参考波形,CSD-FDM(M=10)采用极小时间采样间隔τ=0.1 ms得到；红色为不同差分格式模拟波形. ①②分别为CSD-FDM(M=10), 时间采样间隔τ=0.5 ms,1.0 ms；③④分别为TSD-FDM(M=10)，时间采样间隔τ=1.0 ms,1.25 ms;⑤⑥分别为M2M+N-FD(M=8;N=1), 时间采样间隔τ=1.25 ms,1.5 ms.Fig.13 Numerical modeling single trace waveform on the revised 3D SEG/EAGE salt model(a) and (b) Geophone at (3,650,3) and (560,650,3). ① and ② CSD-FDM(M=10) with τ=0.5 ms,1.0 ms；③ and ④ TSD-FDM(M=10) with τ=1.0 ms,1.25 ms； ⑤and ⑥ M2M+N-FD(M=8;N=1) with τ=1.25 ms,1.5 ms.

表4 修改的三维SEG/EAGE盐丘模型CSD-FDM、TSD-FDM和MG-FDM数值模拟计算效率对比Table 4 Comparison of computational efficiency for numerical modeling on revised 3D SEG/EAGE Salt model with CSD-FDM, TSD-FDM and MG-FDM

盐丘模型数值模拟结果表明，三维MG-FDM(M=8;N=1)能采用更大的时间采样间隔以提高计算效率，同时能更有效地减小数值频散，获得更高的模拟精度.

6 结论与讨论

本文提出了一种面向三维标量波动方程数值模拟的混合网格有限差分方法(MG-FDM)，将三维Laplace微分算子近似表示为坐标轴网格点构建的M个Laplace差分算子和非坐标轴网格点构建的N个Laplace差分算子的加权平均，充分利用了与差分中心点距离更近的非坐标轴网点，有效减小数值频散和数值各向异性，获得更高的模拟精度.通过对比分析，得出以下结论：

(1)三维CSD-FDM的空间差分算子(Laplace 差分算子)虽然能达到2M阶差分精度，但其差分离散波动方程仅具有2阶差分精度；TSD-FDM的差分离散波动方程沿特定的48个传播方向可达到2M阶差分精度，但沿其他传播方向仅具有2阶差分精度；MG-FDM的差分离散波动方程沿任意传播方向可达到4阶、6阶、甚至任意偶数阶差分精度.

(2)相比三维CSD-FDM和TSD-FDM，计算效率基本相同时，MG-FDM的数值频散压制效果更好，模拟精度更高；模拟精度基本相同时，MG-FDM可以采用更大的时间采样间隔以提高计算效率；稳定性分析表明，MG-FDM稳定性更强，为采用更大的时间采样间隔以提高计算效率奠定了基础.

(3)数值模拟实例进一步证实了三维MG-FDM在模拟精度和计算效率方面的优越性，也验证了其普遍适用性.

但是，本文提出的三维MG-FDM也具有一定的局限性，推导过程隐含了三个方向的空间采样间隔相等，即MG-FDM要求采用立方体网格.适用于长方体网格的更具一般性的MG-FDM，有待进一步研究.

附录A 三维坐标系中非坐标轴网格点构建Laplace差分算子