1 引言

波动方程法是地震波场数值模拟的重要方法，其本质是求解地震波动方程.有限差分法^[1-2]是迄今最常用而且较成熟的波动方程正演模拟方法，其算法简单，易于实现，而且精度较高.但是有限差分方法只能使用规则网格剖分，当模拟起伏地表或地层等复杂构造的地质体时，具有很大的局限性，其精度不够理想.随着计算机计算速度的不断提高，内存存储的不断扩大，地震波动方程有限元模拟方法^[3-5]也逐渐发展起来.有限元方法可以使用不规则网格剖分，克服有限差分的上述缺陷，保证了模拟的逼真性.但是该方法在每个时间步上都需要求解大型的线性常微分方程组，造成其计算和存储量巨大、效率相对较低.这些缺点对于模型数据巨大、效率要求较高的地震波场模拟问题来说，尤其不能接受.如果能有效地改进有限元方法的效率，则有限元方法将大有可为.谱元法将有限元方法与谱方法相结合，具有高并行性和高精度的特点，但是目前该方法对网格要求较严格，通用性不好.基于以上优缺点，近年来也出现了一些将有限元和有限差分方法相结合的尝试^[6-7].

间断有限元方法(Discontinuous Galerkin finite element methods，简写为DG)是将有限体积方法和有限元方法相结合形成的一类新型高精度方法，其既具有有限元方法可以模拟任何复杂介质模型的韧性，又具有极高的并行性，能够克服有限元方法计算效率低的缺点.另外，由于结合了高精度算法的斜率限制等技术^[8]，间断有限元方法在处理断裂、高陡等变化剧烈的地质构造时具有明显的优势.该类方法自Reed和Hill^[9]提出之后，在近年来得到了迅速的发展，多种形式的间断有限元方法被提出，其中包括龙格-库塔间断有限元(简写为RKDG)方法^[10-11]、局部间断有限元(简写为LDG)方法^[12-13]、IP方法等，并已经成功应用于计算流体动力学、气象学、海洋学、电磁学等领域.近年来，间断有限元方法在地震勘探领域的应用也受到越来越多的关注，国内鲜有成果见于文献，但国外研究早已蓬勃开展.在弹性波动方程模拟中，Käser等人提出了ADER-DG方法^[14-18]，该方法是将间断有限元方法与任意高阶时间积分方法相结合形成的，可以使得时间离散达到与空间离散相同的精度.在声波方程模拟中，IP方法^[19]和LDG方法^[20]是主要的选用方法.相比于IP方法，LDG方法具有全局和局部能量守恒特性，有助于控制波传播过程中耗散作用.

本文选取了局部间断有限元方法对二维地震声波波动方程进行数值模拟.该方法的主要思想是将高阶微分方程分解成一组耦合的一阶微分方程系统，然后使用RKDG方法求解此系统.在空间上，采用LDG方法进行离散，进而得到一个常微分方程系统，然后选取了蛙跳格式(Leap-Frog scheme)求解此系统，从而得到了LDG数值模拟算法.最后设计了计算实例进行模拟.在算例一中，分别用本文算法和连续有限元方法模拟了lamb问题，并进行了比较，验证了本文算法在计算的规模、精度和效率等方面对连续有限元方法的明显改进效果.算例二验证了该算法对起伏高陡构造的适应性，证实了算法的优越性.

2 方法

考虑二维标量声波波动方程(含震源)：

(1)

其中p=p(x，z，t)表示声波波场，υ=υ(x，z)是声波速度，f=f(x，z，t)表示震源.Ω表示空间计算区域，I=[0，T_s]表示时间计算区间.边界Γ=Γ_N∪Γ_abc，Γ_N和Γ_abc分别表示其自由边界和吸收边界.我们使用了旁轴近似理论吸收边界条件，即用入射波方程代替全波动方程.n=(n_x，n_z)^T为边界的外法向量，上标T表示向量的转置.

LDG方法的主要思想是先对高阶微分项进行“降阶”处理，将空间二阶微分方程分解转化成耦合一阶微分方程组，再利用间断有限元方法来求解.为此，引入新变量

(2)

则方程(1)转化为关于未知量(p，q)^T的微分方程：

(3)

另外还有一种“降阶”方法，即再引入一个变量，则方程(1)转化为

这两种分解方法的LDG空间离散方法是类似的，区别在于后者得到的等价形式与应力-速度方程或者弹性波方程相似，更符合物理意义，而前者可以直接求得未知量p.为了便于与连续有限元方法比较，本文中考虑由(2)和(3)组成的微分方程系统.

记网格剖分为T_h，定义在其上的有限元空间记为S _h^m={φ|φ|_K∈P^m(K)，K∈Τ_h}，P^m(K)表示剖分单元K上的最大次数不超过m的多项式的集合.h表示剖分中所有单元边长的最大值.

由于时间上采用显式格式时，算法在每个单元上的计算都是相同而且相互独立的，因此我们只在某一个单元上推导算法即可.用(·，·)_k和〈·，·〉_e分别表示在单元K上的面积分和在单元边界e上的线积分，将方程(3)的两端同乘以检验函数φ∈ P^m(K)，然后在K上积分，

对于在Γ=Γ_abc∪Γ_N上的线积分，应用(3)式中的边界条件，即对于Γ_N上的积分，以0替代q·n，对于Γ_abc上的积分，以替代q·n，可得

以离散解p_h和q_h分别替换p和q，由于离散解在单元边界上是间断的，因此对于边e上的积分，如果e∈Γ，则直接用q·n在边界上的值替代即可.如果，就以数值通量Φ(q_h)来代替q，可知离散解满足方程

(4)

其中向量Φ(q_h)的定义为

C₁₂是以C₁₂为对角线元素的二阶矩阵，C₁₁和C₁₂为常数.平均{·}和跳跃[[·]]的定义如下：边e为单元K和K^*的公共边，n_k和n_k^* 分别为它们的外法向量，假设v和w分别为标量和向量函数，则在边e上有

方程(4)等价于如下矩阵形式的常微分方程系统

(5)

其中质量矩阵；；刚度矩阵，N为单元K上的未知量个数.最后一项中，向量v_h是由单元K和与K相邻的所有单元的所有解向量组成的，即若假设K₁，K₂，K₃表示与单元K₀=K相邻的三个单元，v_i=(p_h，q_h，x，q_h，z)^T表示单元K_i上的解向量，大小为3N×1，则v_h=(v ₀^T，v ₁^T，v ₂^T，v₃^T) ^T，大小为12N×1.相应的R_p为方程(4)中右端第三项生成的矩阵，大小为N×12N，其形式较为复杂，此处略去.

注意到常微分方程组(5)中同时含有一阶和二阶的时间导数，因此我们选用了蛙跳格式(Leap-Frog scheme)来进行时间离散.令Δt=T_s/N_t表示时间步长，N_t为时间步个数，则对于第n+1(≤N_t)步，有

(6)

类似地，在方程(2)两端同乘以ψ=(φ，φ)^T∈ P^m(K)×P^m(K)，并在K上积分，可得

(7)

其中数值通量为

方程(7)的等价矩阵形式为

(8)

至此，我们的算法可以描述如下：定义p _h⁰为方程(p _h⁰，φ)_k=(p(x，z，0)，φ)_k，的解.假设p _h^-1，q _h^-1均为零向量，则q _h⁰可通过方程(8)求出.假设第n个时间步已经计算结束，则对于第n+ 1步，

1) 并行地执行：在单元K上，计算方程组(6)求出p _hⁿ⁺¹；

2) 并行地执行：将p _hⁿ⁺¹代入方程(8)求出q _hⁿ⁺¹；循环以上步骤，直至计算所有时间步.

3 计算实例

在下面的算例中，均选用了定点上延迟的雷克子波，可描述为

其中(x_s，z_s)为源点坐标，设，

其中峰值频率取为f_p=30Hz.时间步长Δt=0.1ms，采样间隔取为0.5ms，道间距为10m.计算时，采用了最大边长度不大于10m的三角形网格，参数选取为C₁₁=C₂₂=0，C₁₂=1/h，计算设备选用了12核的工作站.

算例1

考虑lamb问题，计算区域为Ω={(x，z)|0≤x≤800，0≤z≤800}，波速υ=2000m/s，震源位置(x_s=400，z_s=400).四条边均取吸收边界条件.计算时长T=0.4s.我们采用连续有限元方法(记为CFEM)和局部间断有限元算法分别求解该问题，对于每种算法，实验了m=0，1，2，3次的多项式作为基函数的情形.为了便于对比，所有的实验都使用了同一个三角形网格剖分.

对于LDG算法，我们实验了串行(记为LDG-Serial)和简单并行(记为LDG-Parallel)两种策略，这里的简单并行指使用多线程并行求解各个剖分单元，使用OpenMP实现线程并行，线程数设置为30.这样的并行方式是基于LDG算法各个剖分单元的计算是独立的特点之上的，CFEM不具备这样的特点，此处不考虑CFEM的并行.事实上，CFEM的并行需要借助于区域分解方法^[21]，算法十分复杂，此方法同样可用于LDG-Parallel策略中，即进程并行的同时，每个进程中进行多线程并行，这种方式可进一步降低存储压力和提升效率，但本文未作考虑.下面从问题规模、计算效率和计算精度三个方面进行对比.

问题规模  理论上，LDG方法的未知量个数d_LDG总是剖分单元个数N_e的倍数，而CFEM方法的未知量个数d_CFEM一般依赖于剖分中单元顶点的个数.分析图 1可以发现，当基函数次数为0时，LDG方法与CFEM方法是等效的，d_LDG和d_CFEM都等于剖分单元的个数N_e.随着次数m的增加，d_LDG总是大于d_CFEM，但二者的比值d_LDG/d_CFEM在不断减小.这表明在三角形网格下，随着次数的增大，两种方法的计算和存储规模的差距在逐渐缩小.

计算效率   图 2描述了CFEM，LDG-Serial和LDG-Parallel的单步计算时间随基函数次数的变换情况，记CFEM，LDG-Serial和LDG-Parallel方法的计算时间分别为t_C，t_S和t_P.由于LDG方法比CFEM方法要多计算在边上的积分，所以即使在m=0时，t_S仍远大于t_C，但是并行之后，计算时间t_P大幅度缩短至与t_C相当.随着m的增大，各个时间都有不同程度的增大，其中t_C的增速最快，并在m=3时超过了t_S.t_P增速最缓慢，t_P与t_C的差距逐渐拉大.因此LDG方法可获得较高的计算效率.

计算精度 图 3中给出了接收点在震源正上方250m处时接收到的波形.LDG、CFEM的模拟结果以及lamb问题解析解分别用红色、蓝色和黑色线条标记，图 3a和3b分别为采用一次元和二次元的结果.从图中可以发现，相比于有限元方法，间断有限元方法在模拟精度和压制频散方面具有更好的效果.

算例2

考虑如图 4所示的起伏地层模型.模型长2000m，深1000m，分为3层，第2层上边界位于400m深处，设置了背斜和向斜构造，其中顶点1和2，3和4的距离为400m，弧高300m，属于高陡构造.第3层位于900m深.自上而下速度依次为2000m/s，3000m/s和5000m/s.震源点坐标为(1000，0).地表取自由边界条件，其余边界取吸收边界条件.为了压制“假频”，基函数采用了二次多项式.

图 5中给出了0.18s，0.25s，0.29s，0.38s，0.43s和0.71s时刻的波场快照，从图中可以发现，地震波在背斜顶部以及1~4四个顶点处产生的绕射波、深水平层的反射波、以及绕射波形成的多次反射波都能清晰地捕捉到.单炮记录如图 6所示，由于使用了旁轴近似理论吸收边界条件，所以边界反射现象未能全部消除.记录中，同相轴清晰，无“假频”干扰，即使深层的反射波也得到了较好的表现.可见，本算法对此类构造具有较好的模拟效果.

4 结论

本文应用局部间断有限元方法离散二维声波方程，在时间上选用了显式的蛙跳格式，形成了一种新型模拟算法.该算法结合了有限差分方法和有限元方法的优点，如可以使用三角形单元，因而可以很好地模拟起伏地形；具有天然的并行性，因而计算效率比有限元方法高出很多；在单个单元内的未知量个数比有限元方法和有限差分方法都多，因而计算精度更高等等.给出的数值算例证实了这些优点.基于以上的优越性，该算法可以进一步推广到弹性波数值模拟以及偏移成像等领域.