0 引言

准确而高效的声波模拟算法对增加我们认识波动方程在三维复杂介质中传播有着重要作用.以双程波方程为基础的全空间离散方法理论上可以处理任意复杂的介质，所以一直以来都是学术研究的热点.以均匀网格进行全空间离散的代表方法为有限差分法(Kelly et al., 1976)，由于网格的排布规则，这类方法的不需要网格的几何描述文件，从而具有易于实现、易于并行的特性，但是该方法不易于处理起伏地表模型.非规则网格进行全空间离散方法能够方便的应对起伏地表，代表方法有有限元法(Dupros et al., 2010; Komatitsch et al., 2010a)，谱元法(Komatitsch et al., 2005; Komatitsch and Tromp, 2002)和格子法等(Gao and Zhang, 2006; Zhang and Liu, 1999).通过非规则网格剖分和非规则网格算子，这类方法能够灵活的处理复杂界面问题，但是这类方法需要网格的几何描述文件记录网格格点位置，需提前剖分模型网格，较有限差分方法的数值实现较为复杂.

近年来，通用图像处理器(General Purpose Graphic Processing Unit, GPGPU)的发展提升了科学计算的高性能计算的计算能力，通过单指令多线程(Single Instruction Multiple Threads, SIMT)的计算框架实现计算加速.相比于中央处理器(Central Processing Unit, CPU)善于处理流控制和复杂的数据缓存，图形处理器(Graphic Processing Unit, GPU)则更擅长大规模的简单计算.正演模拟本质上是一个简单的、局部的正演算子作用于所有计算网格单元上，非常符合GPU的加速.近年来亦有很多针对GPU优化的波动方程正演方法提出，对原有的正演方法进行了加速，如有限差分法(Abdelkhalek et al., 2009; Micikevicius, 2009)、谱元法(Komatitsch et al., 2010b)等，并且加速效果非常显著.

格子法(Zhang and Liu, 1999)是一种基于非结构网格的弹性波正演模拟方法，它具有灵活处理起伏地表，自然满足自由表面边界条件，并且能够灵活处理裂缝介质(Zhang, 2005; Zhang and Gao, 2009)、高速间断面(Zhang, 2004a)，各向异性介质(Gao and Zhang, 2006; Zhang and Verschuur, 2002)，固液耦合介质(Pulkkinen et al., 2014; Zhang, 2004b).该方法发展了基于二维三角形(Zhang and Liu, 1999)、二维四边形(Zhang, 1997)、三维四面体(Gao and Zhang, 2006)网格的波动方程正演算子.双重网格(Yang et al., 2017)的提出极大的减少了需要网格剖分的个数，该方法结合了结构化网格和非结构化网格各自的优点，通过全局非结构网格局部结构化网格减少了几何描述文件的存储量，适用于GPU并行计算，从而提升了格子法的计算效率.二维双重网格给非规则网格正演模拟提供了新的思路，三维波动方程模拟的几何描述网格要比二维网格几何描述文件大一个数量级，能够更好的发挥双重网格带来的优势.

本文利用了六面体网格的拓扑特性，提出了基于非结构六面体双重网格格子法的声波方程模拟方法，该方法一方面继承了格子法灵活处理起伏地表、自由表面边界条件的特性.另一方面该方法提出了三维双重网格解决方案，使用GPU加速，提升了三维传统格子法的计算效率.本文的组成如下，首先我们将介绍三维六面体网格格子法的基本原理；接下来将介绍六面体双重网格的几何特性和GPU优化策略，然后通过三个数值算例验证本方法的有效性.最后，我们讨论该方法之后发展的方向和结论部分.

1 基本理论

格子法是一个基于非规则、非结构网格的波动方程数值模拟算法，该方法能够灵活处理起伏地表、高速间断面等，是一种灵活的波动方程正演模拟方法.

1.1 控制方程

忽略震源项，非均匀介质下的一阶声波声压-速度方程可以表示为

(1)

(2)

其中，κ(x)=ρv²为体涨模量，ρ为密度，v代表速度，b(x)=1/ρ(x)为浮力模量，v_i(x, t)分别为x_i方向分量的速度场，p(x, t)为声压场，在式中，i=1, 2, 3为爱因斯坦求和哑标.

(1) 式两边除以κ(x)，同时进行体积分，公式左边假设在积分区域Ω内为一个常量，等式右边利用高斯-格林公式，可以得到如下形式：

(3)

其中Ω为积分区域，为积分区域的外表面，α_i为外法向余弦在x_i方向的分量，V为积分区域的体积，S为积分包络面的面积，i=1, 2, 3为爱因斯坦哑标.格子法的空间离散形式如图 1所示，每个网格单元沿着网格单元中心，面中心和各边中心将将六面体网格单元分割为8个子部分，其中在网格内部有12个子平面.声压p定义在网格格点上(NODE A, B, …, H)，速度定义在每个子平面的中心处(point 1, 2, 3, …)，体涨模量和浮力模量在网格的中心处，并假设它们在网格单元内为一个常量.对于每个网格点的积分区域为包含该点所有单元的所有子部分区域的集合，以点A为例，点A在此网格的积分区域为黄色部分所围成的区域，速度场v在该区域的外表面.可以看到，该单元的每个网格点在其单元内对应的积分区域上内都定义有3个速度分量的离散点(点1，2，3).

按照上述方式对(3)式进行空间离散，可以得到离散形式的格子法控制方程，控制方程可以表示为

(4)

其中m为点A附近的网格单元数，l遍历包含格点A的所有网格单元，k遍历每个网格单元的3个离散点，即速度v的定义位置.

1.2 等参变换求导

定义于网格内部的速度v通过该网格单元的8个声压得到，这里面涉及到计算声压p的一阶空间导数.非规则四边形网格有限差分法(Zhang, 1997)通过等参变换计算四边形网格内部的空间一阶导数，该方法首先将笛卡尔坐标系下的非规则网格投影到等参坐标系下的规则网格中，然后在等参坐标系下进行求导，然后将该导数再变换回笛卡尔坐标系中，该方法能够使得计算参量在网格边缘处连续.我们将该方法推广至三维六面体网格中，六面体网格等参变换如下所示：

(5)

其中p为待求导的参量(声压)，x_i为笛卡尔坐标系坐标，ξ_i为等参坐标，i, j=1, 2, 3，等式右手边第一项为等参变换的雅可比矩阵，其中等参变换关系满足线性插值的关系，表示如下：

(6)

其中，p_i, x_i(i=1, 2, …, 8)分别为8个格点的声压值和坐标值，N_i为等参变换的形函数，形函数为ξ_i的函数，表示如下：

(7)

将式两侧乘以雅可比矩阵的逆矩阵，可得：

(8)

将(6)、(7)式代入式(8)，再将导数结果代入(2)式中即可求得速度.

1.3 时间离散

对于时间离散，格子法采用向前差分格式离散方程.格子法采用吝啬存储(Luo and Schuster, 1990)的方法，即只存储声压标量场而不存储速度矢量场，从而很大程度上节省内存开销.将格子法控制方程(4)、(2)，重写为如下形式：

(9)

(10)

其中Δt代表时间差分的时间采样间隔.格子法时间迭代过程如下：

1) 通过非规则网格一阶导数公式(8)计算得到；

2) 将代入(10)式计算更新速度场v_i(x, (n+0.5)Δt)；

3) 将v_i(x, (n+0.5)Δt)代入公式(9)更新声压p.

2 双重网格

双重网格最先被提出来用于减少二维三角网格格子法的三角网格的几何描述和GPU加速(Yang et al., 2017).双重网格是对非结构网格进行二次剖分，利用利用三角形的拓扑分割特性，生成一种结构和非结构网格合二为一的网格，即一种非结构粗网格和局部结构化的细网格，非结构粗网格用于描述地质构造，细网格用于满足频散关系，进行控制方程的计算.双重网格格子法具有以下三个优点：1)双重网格只需要存储粗网格的几何描述，细网格格点的坐标均可以通过拓扑关系计算得到，从而减少网格描述的IO；2)双重网格适用于粗粒度并行，可应用于GPU加速；3)局部结构化的细网格有利于合并访存，从而减少内存寻址的时间.

2.1 减少几何描述

双重网格三角形网格具有全等的特性，即三角粗网格下所有的三角形细网格均为全等.图 2展示了双重三角形网格.

然而在四边形网格和六面体网格双重网格中，粗网格中的四边形(六面体)细网格单元并不具有全等的特性.然而四边形(六面体)的细网格格点坐标可以通过等参变换建立关系，六面体网格细网格格点坐标拓扑关系为

(11)

其中，为等参变换形函数(7)式的ξ₁, ξ₂, ξ₃的值，x^lmn为向量表示序号为l, m, n的细网格格点的坐标，x_k为向量表示的粗网格8个格点的坐标，N为分割数.图 3展示了分割数为3的六面体双重网格.如此可以避免存储所有细网格的网格几何信息，从而减少了网格剖分的工作量，同时减少了几何描述文件的大小.

双重网格格子法的控制方程作用在所有细网格，因而需要所有细网格的雅克比矩阵，这里存在一个存储和计算的平衡问题.存储所有细网格的雅克比矩阵将耗费大量的存储量，而在每一步时间迭代临时计算雅克比矩阵将耗费大量的计算量，所以需要在存储量和计算量达到良好的平衡.我们利用双重网格的几何拓扑关系，可以通过粗网格和拓扑关系来计算得到细网格的雅克比矩阵.

将公式(11)代入公式(6)中可得：

(12)

其中x_ik(x_i^coarse, l, m, n)为粗网格内第(l, m, n)个细网格的第k个格点的i分量，x_i(ξ)为细网格的坐标的i分量，将(12)式代入雅克比矩阵中可以得到：

(13)

公式(8)的雅克比逆矩阵可以由伴随矩阵和行列式相除得到，表示为

(14)

其中J^*为伴随矩阵，|J|为雅克比矩阵的行列式.将(13)、(8)式，代入公式(14)，并将结果表示为三部分多项式相乘的形式，第一部分为粗网格的坐标信息，第二部分为细网格在粗网格下的位置信息，第三部分为细网格内等参坐标信息，表示为

(15)

其中i, j=1, 2, 3为雅克比矩阵的第ij项，ξ_a为离散点在细网格的等参坐标，P_ij(), Q()代表多项式.将公式(15)代入公式(14)中可得：

(16)

在数值实现中，我们将(16)式分子、分母的第一项作为粗网格的几何系数保存在内存中，经过化简同类项可以得到个23×9+65=272独立参数，通过这些参数就可以计算得到所有细网格的雅可比矩阵逆矩阵，同时这些参数与分割数无关.

2.2 双重网格GPU配置

双重网格同时适用于粗粒度并行，方便进行GPU加速.Nvidia的GPU使用的CUDA编程架构采用了分层资源管理机制，即grid，block和thread.每个block分配了共享内存(shared memory)，可以被该block下的thread共同访问.共享内存具有高速率和低延迟的特点，但是共享内存的容量有限.双重网格非常适用于这种分层资源配置，在双重网格配置中，grid对应于非规则网格，block对应处理每一个粗网格，thread对应于每一个细网格.上面提到细网格的几何描述可以通过粗网格拓扑关系得到(如式(16)中的参数)，这些粗网格的几何描述文件被存放在共享内存当中.在进行控制方程运算时，每个thread可先根据粗网格的几何文件计算得到细网格的几何描述，然后进行该细网格的控制方程时间迭代.该GPU配置可以提高并行的粒度，粗粒度的GPU程序可以获得更高的加速比，图 4展示了三维六面体双重网格格子法的GPU配置.

2.3 合并访存，减少内存寻址的时间

细网格的几何描述通过拓扑描述得到，每一个细网格格点需要有一个唯一编号来确定格点在内存当中的位置.根据细网格在粗网格的相对位置，双重网格每一个粗网格由角点、边点、面点和内点组成，双重六面体网格格点情况请见图 3.其中，除了角点有几何描述文件直接给出编号以外，其他三种格点(边点、面点和内点)的编号通过计算获得.

从图 3可以看出，网格点在内存的排布为，首先是角点，其次是边点，接下来是面点，最后存放内点.它们各自的编号通过表示为

(17)

其中，N为分割数，l, m, n为细网格在粗网格的相对位置值，N₀, N₁, N₂分别为角点，边点和面点的个数，f(l), f(m(N-1)+l)为边点和面点序号的模板.

在非结构网格当中，每个网格周围网格单元数不固定.在双重网格结构下，内点仅仅属于一个粗网格内部，所以他们的顺序可以仅由该粗网格的局部顺序决定.而边点和面点与其他粗网格的边点和面点相连，而它们在不同的粗网格局部顺序缺少一个结构化特征来保持一致.所以这两种格点(边点和面点)需要一个不依赖于网格单元的方式来确定，我们使用几何描述文件来唯一确定边点和面点的排列顺序，并通过模板将每个粗网格内的局部编号来转换为他们内存中的排列顺序.在双重网格格子法当中，每条边，每个面都有其唯一编号，并且它们分别存储了该边(该面)的点和边的信息.图 5展示了边点和面点的序号模板.

对于每一条边，每一个面，每个粗网格内部的内点，它们的内存地址都是连续的，从而实现了合并访存的效果，减少了GPU在内存当中的寻址时间和寻址次数，降低了延迟提升了计算效率.

3 数值算例

本节将通过3个数值算例来验证该方法的有效性.首先我们使用均匀模型进行理论验证，并测试不同分割数下双重网格的GPU计算时间.第二个模型为一个起伏地表的理论模型，通过该模型我们论证该方法能够适应起伏地表.第三个模型为SEG/EAGE C3盐丘模型，通过该模型我们将验证该方法能够处理速度间断面的特性.

3.1 均匀模型

第一个模型为均匀模型，模型大小为2 km×2 km×2 km，速度大小为3000 m·s^-1，密度为1200 kg·m^-3，细网格尺寸为10 m，时间步长为1.5 ms，震源为主频为15 Hz的雷克子波，模拟时间步长为500步(0.75 s)，使用的GPU为Nvidia GTX 980Ti.震源位置位于模型的中心处，我们分别记录了距离震源200 m，400 m和600 m位置处的波形，波形图如图 6所示.从图 6中可以看出，该方法在10 m，1.5 ms网格间距和时间步长的情况下没有出现失稳的情况，并且几乎没有数值频散，说明该方法的稳定性.

接下来比较不同分割数下的GPU模拟时间，我们使用的GPU为Nvidia GTX 980Ti.我们固定细网格大小和数量，通过调整分割数和粗网格的数量，模拟同一大小的均匀模型，使用Nvidia Visual Profiler分析每一个时间步所花费的时间(ms).模拟时间如图 7所示.

从图 7中可以看出在分割数为2的时候，每个CUDA block下的CUDA thread数量较小，并行度不够，计算时间为1351 ms.可以见到在分割数为3时，粗网格数量明显减少，并且计算时间从1351 ms降低到750 ms.在分割数为4时，即每个粗网格分割为64个细网格.根据CUDA架构，每次有32个运算单元(warp数)计算CUDA thread任务，分割数为4时刚好为32的2倍，即保证了GPU运算单元随时保持满载，进一步提升了计算效率.在分割数大于4的情况下，计算已经成为了该算法的主要瓶颈，计算性能不能进一步提升，并且warp数的非对齐会导致GPU的CUDA核不能保证时刻满载.在分割数为5和6时，计算时间分别为391.81和428 ms.当分割数大于6时，受到GPU寄存器的限制，程序无法正常运行.我们可以得到，一个分割数为4的双重网格配置可以达到计算效率的最大化.

3.2 起伏地表模型

第二个模型为起伏地表模型，我们在地表放置了一个半圆形柱的槽，并将震源放置在了地表的中心位置来验证该方法能够处理起伏地表.模型如图 8所示，模型大小为(1960 m, 1960 m, 960 m)，速度为3600 m·s^-1，密度为1242 kg·m^-3，震源位置位于(960 m, 960 m, 0 m)，震源类型为主频为20 Hz的雷克子波，速度为圆柱的圆半径为90.5 m，长度为320 m.我们抽取了4个平面，x=960 m, y=880 m, z=0 m, z=120 m，并展示了这4个切面在0.2 s时刻的波场快照，如图 9所示.

从图 9中可以看出，由于该方法通过自适应的非规则网格而非阶梯近似刻画起伏地表和速度间断面，所以不会出现阶梯近似带来的尾波效应.通过该算例证明了该方法能够高精度的模拟起伏地表问题.

3.3 SEG/EAGE C3模型

C3模型是SEG/EAGE于1996年提出的三维盐丘模型，我们使用该模型的部分模型来论证方法能够处理复杂三维介质.模型大小为3.3 km×1.6 km×1 km, 模型的速度如图 10所示，震源为主频20 Hz的雷克子波，时间步长为1 ms，并于在该模型的6个表面加入了40 m的吸收层，吸收边界使用的是非分裂卷积完美匹配层吸收边界条件(Komatitsch and Martin, 2007).

图 11展示了SEG C3模型0.45 s的波场快照，从波场快照中可以看到，波前面在经过盐丘的时候产生了明显的变化，该方法能够处理高速间断面，同时能够处理复杂地质构造.

六面体网格相比于四面体网格具有高稳定性，高精度的特点.对于相同精度的波动方程模拟，六面体网格允许比四面体网格更大的网格尺寸.对于剖分同一个地质体模型，六面体网格的数量要远小于四面体网格的数量.这些特点都表明了六面体良好的几何性质.但是非结构四面体网格的网格剖分理论部分已经比较成熟，使用维诺图(Voronoi Tessellation)可以实现非结构四面体网格对于任意复杂结构进行剖分(Du et al., 2010; Du and Wang, 2003).但是目前非结构六面体网格剖分仍存在很多问题没有解决，在数学界被称为‘神圣网格’.虽然目前已经有一些六面体网格剖分软件的出现，但是它们多用于较为规则的剖分而不是复杂地质体，这一部分的工作需要我们进一步的研究.

4 结论

我们实现了基于非规则六面体网格的三维格子法并且使用双重网格的网格配置实现了GPU加速.首先，使用该方法继承了格子法能够处理起伏地表、刻画速度间断面，满足自由表面边界条件的特点.第二，双重网格极大的提高了该方法的运行速度，减少了网格几何描述，实现了GPU粗粒度并行，通过合并访存降低了延迟提高了效率.通过数值算例和理论分析，我们得到了分割数为4的双重网格能够得到最优的加速比，同时使得几何描述文件大小是常规非规则网格几何描述文件的，极大的减少了程序的IO操作.通过三个数值算例，均匀模型、起伏地表模型和SEG/EAGE的三维盐丘模型，我们分别说明了该方法的精确性，起伏地表的特性和能够处理高速间断面的特性.我们讨论的非结构六面体网格剖分仍存在一些问题，该方法的实用性仍需要进一步的研究.