引言

起伏界面波场复杂，在不规则界面上既产生波的反射、散射和绕射，又产生面波，更有受界面起伏特征影响的散射波等（孙建国，2007）。一般情况下，为处理复杂边界，需增加很多甚至是数倍于规则边界的计算量。为精确模拟起伏界面条件下的地震波传播特征，需解决两个关键问题——计算方法选择与模型网格精确剖分。近年来，提出的数值模拟计算方法包括有限差分法（FEM）（Levander，1988；Olsen等，1995，1997；Graves，1996；Olsen，2000；Moczo等，2001）、伪谱法（PSM）（Patera，1984；Furumura等，1998）、有限元法（FDM）（Smith，1975；Koketsu等，2004；）、谱元法（Patera，1984；Seriani，1998；Komatitsch等，1998，1999，2000，2002）。网格剖分方式可分为网格法和非网格法两大类（孙建国，2007）。网格法分为规则网格（Moczo，1989；Jastram等，1994）、非结构网格（Tessmer等，1992；Zahradník，1993）和贴体网格（垂向变换）（Tessmer等，1992；Hestholm等，1994）等。近年来，国内外对于网格法开展了系列研究工作。规则网格差分法最早由Alterman等（1968）引入地震波数值模拟计算中，并利用规则网格差分法完成了层状介质中二阶弹性波波动方程的模拟，随后Alford等（1974）分析了网格选取与模拟精度之间的关系。Madariaga（1976）提出了交错网格差分格法，并很好地应用于正演模拟中，但该方法使用低阶差分，精度较低，随后Dablain（1986）提出了任意阶差分精度算法，Crase等（1990）实现了高阶差分算法，因此有限差分法数值模拟精度被大大提高。但采用规则网格剖分时，如果界面起伏较大，会导致模型界面阶梯状离散，进而产生数值散射（Xue，2014）。六面体（四边形）单元的谱元法已成功应用于模拟波的传播(Priolo等，1994；Komatitsch等，2002，2005)，在处理不规则界面时，显示均匀介质内差异过大的网格会激发强的网格反射波（Zhou等，2010），由于其均基于规则网格实现，无法从根本上消除阶梯状网格引起的虚假绕射波，为此学者们发展了不规则网格有限差分法。

不规则网格有限差分法在有限元法中的非结构网格适用性较好，易生成网格，但对精度有一定影响，三角形单元增加了存储量，降低了效率和准确性（Komatitsch等，2000）。贴体网格（垂向变换）思路相近，依照实际地表起伏剖分曲线网格，使用较方便，适用于同位网格和交错网格（Hestholm等，2002；Zhang等，2006），但对网格变形范围具有一定限制（Tessmer等，1992）。Pitarka（1999）提出了各向同性介质中矩形不规则网格有限差分法。Hestholm等（2002）和Tarrass等（2011）提出了以贴体网格的方式对起伏地表进行剖分。Zhang等（2006）和Lan等（2011）利用同位网格的差分方式对曲坐标系下的波动方程进行差分求解，实现了正演模拟，且取得了较高精度。Zhang等（2019）以不同网格分辨率的模拟，强调了复杂情况下大规模、高分辨率地震波模拟的重要性。

1.   谱元法

本文通过谱元法分析六面体单元剖分起伏界面的不同方式对地震波场的影响。对于地下介质中传播的地震波，可由波动方程、初始条件和边界条件进行描述，即在$ \mathrm{\Omega } $中地震波位移满足以下方程：

式中，$ {\boldsymbol{u}} $表示位移向量，$ {\boldsymbol{\sigma}} $表示应力张量，$ \boldsymbol{C} $表示刚度张量，$ \rho $表示密度，$ f $表示外力，“：”表示张量相乘。

在谱元法中，方程（1）点乘任意试验权函数$ w $，并进行分步积分，可得到：

式中，$ \mathrm{\Omega } $表示物理计算域，$ \Gamma $为其边界，包括自由地表边界$ {\Gamma }_{{\rm{f}}} $和人工吸收边界$ {\Gamma }_{{\rm{abs}}} $。

由于地表应力为零，$ {\Gamma }_{\mathrm{f}} $在部分过程中消失，因此式（2）自然满足地表条件。

在参考单元Λ中，引入N阶多项式的基函数，其在单元$ {\mathrm{\Omega }}_{\mathrm{e}} $上对应1组节点，$ {u}_{N}^{{\rm{e}}} $和$ {w}_{N}^{{\rm{e}}} $分别为$ u $和$ w $在这些节点上的Legrange插值。这些节点在[−1，1]区间内取为Gauss-Lobatto-Legendre点，是下式的根：

式中，${P}_{N}^{\mathrm{{'}}}\left(\xi \right)$为Legendre N阶多项式的偏导数。

在参考单元Λ中，$ {u}_{N}^{e} $表示为：

式中，$ h\left({\xi }_{p}\right) $表示$ p $项一维Legrange插值。

引入多项式近似，方程（2）的积分可利用Gauss-Lobatto-Legendre积分规则得到：

式中，$ {w}_{i} $为积分系数，$ {\boldsymbol{J}}_{i} $为Jacobian矩阵，由影射$ {\mathrm{\Omega }}_{\mathrm{e}} $到Λ的关系$ {\Gamma }_{e} $决定。空间离散化后，进行时间离散化，用Ｕ表示介质位移未知量，得到：

式中，$ \mathit{M} $为质量矩阵，$ \mathit{K} $为刚度矩阵，Ｆ为源项。

Legendre谱元的不重叠单元的空间离散，以Gauss-Lobatto-Legendre节点插值多项式，结合Gauss-Lobatto-Legendre积分规则，使质量矩阵$ \mathit{M} $为对角线矩阵，$ \mathit{K} $为带宽矩阵。$ \mathit{M} $的对角性使谱元法计算效率较有限元法有了较大提高。

2.   3种剖分方式的建立

2.1   问题的提出

常规的网格剖分方式中多采用规则的阶梯状网格离散地质模型（Dablain，1986；Levander，1988；董良国等，2000，2004），如图1所示。在进行数值模拟计算时，当遇到有起伏界面的地质模型时，直接使用矩形网格离散界面会导致阶梯状的空间离散。如果网格剖分过于稀疏，如图2（a）所示，单元尺寸差异很大，在计算中会导致发散；如果网格剖分过于密集，如图2（b）所示，在加密网格后采用台阶状网格近似描绘起伏地表（Robertsson，1996）时，会大大增加计算量，占用计算机内存，且无法彻底消除阶梯的散射问题。为解决上述问题，须寻找使网格尽可能保持均匀、剖分尽可能贴合起伏界面的剖分方式。

图 1  矩形网格剖分起伏界面

Figure 1.  Rectangular mesh generation of undulating interface

下载: 全尺寸图片 幻灯片

图 2  矩形网格剖分示意图

Figure 2.  Schematic diagram of rectangular mesh generation

下载: 全尺寸图片 幻灯片

（1）问题一：深度方向（z轴方向）不同单元个数之间的过渡

取起伏界面处的1个台阶状网格进行研究（图3），为使网格贴合起伏界面，并保证均为正交或略倾斜的六面体单元，需考虑深度方向（z轴方向）不同单元个数之间的过渡。由图3可知，二维矩形网格剖分产生了阶梯状界面。

图 3  台阶状网格截面示意图

Figure 3.  Schematic diagram of the stepped mesh section

下载: 全尺寸图片 幻灯片

三维模型在地形等高线图上所处位置如图4（a）所示，其中x轴为南北向，y轴为东西向，z轴为深度方向，模型东西长6 000 m，南北长5 000 m，深度为1 100 m。起伏界面位于深度600 m处，最深处为700 m，起伏界面将模型分为上下两部分，图4（a）中红色外框对应上半部分模型深度为600 m的范围，红色内框对应上半部分模型深度为700 m的范围。图4（b）中z轴方向上红色外框范围的单元较红色内框范围单元少1个，产生了阶梯状网格。考虑到三维模型在2个方向（xz平面与yz平面）上均有阶梯状网格（图5），因此不同网格个数之间的过渡在2个方向上均需进行。当2个方向上单元过渡后，将拐角处（图5中点A、B、C、D）的单元进行了2次叠加二合一处理，导致单元发生扭曲畸变。

图 4  三维模型在地形中所处位置示意图

Figure 4.  Location of the 3D calculation model in terrain

下载: 全尺寸图片 幻灯片

图 5  起伏界面上与下的阶梯状网格完成单元二合一处理后变为斜面网格示意图

Figure 5.  The stepped mesh above and below the undulating interface become inclined grids using two-in-one unit method

下载: 全尺寸图片 幻灯片

（2）问题二：扭曲单元的衔接

当2个方向（xz平面与yz平面）均完成问题一中深度方向上的单元个数过渡后，需考虑2个方向拐角处因进行了2次过渡叠加而发生扭曲的单元衔接问题。

2.2   解决方案

本文引入TrueGrid网格前处理软件，针对问题一，对含起伏界面的研究区域进行六面体单元网格剖分，在固有台阶状网格处使用六面体单元二合一的方式，如图6（a）、6（b）所示，将1个台阶状单元网格通过扯动单元的方式，通过1个变形的六面体单元进行过渡衔接，将深度方向（z轴方向）上的不同单元个数进行过渡衔接。在上下台阶网格处使用该方法，将图3所示台阶状网格变为图6（c）所示的连续网格。

图 6  台阶状网格处使用单元二合一的方式处理过程示意图

Figure 6.  process of using two-in-one unit method on stepped mesh

下载: 全尺寸图片 幻灯片

完成单元二合一处理后的网格模型如图7所示，由图7可知，模型消除了台阶状网格，将台阶状网格变为了斜面网格，可将斜面投影至起伏界面，使网格完全贴合界面。但在4个拐角处，由于叠加了2次二合一处理，使拐角处单元发生扭曲畸变，因此，需单独处理4个拐角处的单元，即解决问题二。

图 7  起伏界面上与下的阶梯状网格完成单元二合一处理后变为斜面网格的模型示意图

Figure 7.  The stepped mesh above and below the undulating interface become inclined grids using two-in-one unit method

下载: 全尺寸图片 幻灯片

本文采用3种剖分方式，以对比不同方式的影响。由于4个拐角处理方式一致，本文仅以1个拐角为例。

（1）剖分方式一：在拐角处删除四纵列单元

在网格剖分时单元个数变化策略中，将16个单元合并为12个是网格变化策略中常用的方式（蒙茂洲，2010），如图8所示，该方法可使网格处理较规整。拐角处单元发生的扭曲畸变如图9（a）所示，删除四纵列扭曲单元如图9（b）所示，采用将16个单元合并为12个单元的方法调整其余单元位置，补充至删除的四纵列单元的位置上。为满足计算要求，略调整了12个单元的中心位置，如图9（c）所示，至此完成了剖分。

图 8  16个单元合并为12个单元的过程示意图

Figure 8.  The process of merging 16 units into 12 units

下载: 全尺寸图片 幻灯片

图 9  剖分方式一的过程示意图

Figure 9.  The process of mesh generation of the first method

下载: 全尺寸图片 幻灯片

（2）剖分方式二：在拐角处删除一纵列单元

将16个单元合并为12个单元的网格数量调整方法删除网格数量较多，且调整后倾斜的单元数量较多，采用单元个数变化策略中的将4个单元合并为3个单元的方式，如图10所示，局部调整方式可减少倾斜的单元个数（罗宏伟等，2013）。拐角处的单元发生扭曲畸变，如图11（a）所示，删除一纵列扭曲单元后如图11（b）所示。采用将4个单元合并为3个单元的方法，调整其余单元的位置，补充至删除的一纵列单元的位置上，如图11（c）所示。

图 10  4个单元合并为3个单元的过程示意图

Figure 10.  The process of merging 4 units into 3 units.

下载: 全尺寸图片 幻灯片

图 11  剖分方式二的过程示意图

Figure 11.  The process of mesh generation of the second method

下载: 全尺寸图片 幻灯片

（3）剖分方式三：将图7中的斜面部分单独剖分为1个过渡环

过渡环内、外块体相差1个单元数，即1个台阶。过渡环内、外侧网格充分贴合完成剖分，如图12所示，可知过渡环的设置使拐角处单元变形较小。

图 12  剖分方式三的过程示意图

Figure 12.  The process of mesh generation of the third method

下载: 全尺寸图片 幻灯片

由TrueGrid软件可进行3种剖分方式生成的整体网格质量评价，主要指标包括单元扭曲度及正交性等，评价结果如图13所示。由图13可知，网格质量良好，基本为正交的六面体单元。扭曲度是检查相对的2个角之间的差值，该值越小，单元扭曲度越小，由图13（a）可知，大部分单元相对2个角之间的差值为0。正交性是检查每个单元节点处的角度与90°的差值，该值越小，单元正交性越好，由图13（b）可知，大部分单元节点处的角度与90°的差值为0。

图 13  整体网格质量评价结果

Figure 13.  Evaluation results of overall mesh quality

下载: 全尺寸图片 幻灯片

3.   数值模拟结果对比

3.1   模型参数确定

本文选择的研究区域（图14红色框）位于京津冀地区的第四系地层局部起伏界面，附近地震活动频繁，是强震灾害多发区。根据物探及钻孔资料（金安蜀等，1980；国家地震局，1986；孙若昧等，1996；赵金仁等，1999；王夫运等，2005），在既有北京盆地模型（高孟潭等，2002）的基础上叠加了地形条件，建立了北京盆地含起伏地形的地下介质三维速度结构模型（付长华等，2012）。x轴为南、北方向，y轴为东、西方向，z轴为深度方向，整体深度为1100 m，盆地最深处为700 m，深度分布如图15所示，剖分网格的单元长度为100 m，表1所示为介质和震源模型参数，震源时间函数通过主频为4 Hz的Ricker子波表示，模拟的整个地震波传播时长为4 s。在三维地震波场模拟中计算量仍较大，且对储存能力的要求较高，为此，使用天河二号超级计算机进行并行计算，计算区域被剖分为9个进程。整个区域除自由地表外，其余表面均设置为吸收边界，以消除地震波在区域边界的反射。

图 14  研究区域第四系沉积地层地形

Figure 14.  Topographic map of quaternary sedimentary strata in the study area

下载: 全尺寸图片 幻灯片

图 15  研究区域深度分布

Figure 15.  Depth distribution map of the study area

下载: 全尺寸图片 幻灯片

表 1  模型基本参数

Table 1.  Basic parameters of the model

模型尺寸/m	总步长/步	格点间距/m	时间步长/s	震源时间函数主频/Hz	震源时间函数
6 000×5 000×1 100（长×宽×深）	4 000	100	0.001	4	主频为4 Hz的Ricker子波

下载: 导出CSV 

| 显示表格

3.2   模拟结果

观测点布置位置如下：①测点位置一如图16中红色点所示，3种剖分方式最大的区别在于拐角处理，在3种方式网格模型中六边形单元畸变最大的位置，即在拐角处布置4组测点如图17所示，在3种模型中取同一位置，以反映在拐角连接处不同剖分方式在拐角处的单元模拟结果的差异；②测点位置二如图16中蓝色点所示，在建立的三维网格模型中水平地表布置测点，在3种模型中取同一位置，以反映在水平地表处不同剖分方式对距拐角不同距离处的单元模拟结果的差异。

图 16  模型中测点位置的俯视图

Figure 16.  Top view of test points position in the model

下载: 全尺寸图片 幻灯片

图 17  模型中测点位置一的截面图

Figure 17.  Section view of test point position one in the model

下载: 全尺寸图片 幻灯片

3种剖分方式在测点位置一处的模拟结果对比如图18、19所示。根据Ricker子波理论波形可知，起伏界面位于模型600 m深度附近，当波传至起伏界面时，剖分方式二、三波形较规则，剖分方式一在0.2、2 s左右波形出现了不规则散射，剖分方式二、三时程曲线较吻合。当模型取均匀介质，即介质波速比为1∶1时，波的持时为2.7 s左右；当模型取非均匀介质，即介质波速比为0.8∶1时，波的持时为3.3 s左右，这符合软弱介质使持时增长的原理。

图 18  均匀介质中3种剖分方式在测点位置一处四纵列单元的模拟时程曲线结果对比

Figure 18.  Comparison of the simulation time history curves of three mesh generation methods in a uniform medium with four rows of elements at the test point one

下载: 全尺寸图片 幻灯片

图 19  非均匀介质中3种剖分方式在测点位置一处四纵列单元的模拟时程曲线结果对比

Figure 19.  Comparison of the simulation time history curves of three mesh generation methods in a four-row unit at a test point in heterogeneous medium

下载: 全尺寸图片 幻灯片

由图18可知，在深度700 m台阶网格上下附近，波形出现不规则的散射。在网格中取相同位置的单元直接使用阶梯状网格模拟的时程曲线如图20所示，对比图18、19可知，剖分方式二、三改善了台阶处网格在数值模拟时产生的散射现象。

图 20  非均匀介质中直接使用台阶状网格在测点位置一处四纵列单元的模拟时程曲线结果对比

Figure 20.  Comparison of simulation time history curve results of four rows of units at a test point directly using stepmesh in non-uniform medium

下载: 全尺寸图片 幻灯片

测点位置二如图16中蓝色点所示，分布于模型水平面的8个点位，每个点位间隔100 m，且从1点至8点离拐角单元距离越来越近。3种剖分方式在测点位置二处的模拟结果对比如图21所示，由图21可知，对于水平地表单元而言，当距拐角900 m左右时，3种剖分方式时程曲线基本相同，且距拐角距离越远，时程曲线越吻合；剖分方式二、三时程曲线较一致，而当距拐角距离＜300 m时，剖分方式一存在较大差异。

图 21  3种剖分方式在测点位置二的时程曲线对比

Figure 21.  Comparison of the time history curves of the three mesh generation methods at the second test point

下载: 全尺寸图片 幻灯片

4.   结论与展望

采用谱元法进行数值模拟的优势在于计算精度高、可较好地适应起伏界面、可克服有限元计算量大的缺点，因此被广泛使用。在以往的研究中，通常缺少灵活的六面体单元贴合起伏界面的剖分方式。利用谱元法规则六面体单元进行网格剖分时，界面起伏较大处会导致阶梯状离散，进而产生数值散射。为消除阶梯网格对起伏界面地震动模拟计算结果的影响，本文采用的剖分方式使网格尽可能保持均匀，且剖分尽可能贴合起伏界面。本文取得的研究进展主要有：

（1）本文基于TureGrid软件提出了模拟计算使用的起伏界面处地形六面体单元网格剖分方式，该剖分方式在固有台阶状网格处使用六面体单元二合一的方式将深度方向上的不同单元个数进行过渡衔接。采用Fortran语言编写了应用程序，通过该应用程序可快速建立起伏界面处的地形六面体网格模型。考虑到三维模型在2个方向（xz平面与yz平面）上均有台阶状网格，本文采取了删除四纵列拐角单元、删除一纵列拐角单元及构造单元过渡环3种剖分方式，解决2个方向上单元二合一后拐角处产生扭曲单元的问题，做到了六面体单元尽可能正交，网格尽可能均匀，剖分尽可能贴合起伏界面。

（2）本文通过3种剖分方式建立整体六面体网格模型，通过TrueGrid软件进行整体网格质量评价，单元扭曲度及正交性评价结果显示网格质量良好，基本为正交的规则六面体单元，符合网格剖分要求，满足谱元法计算需求。

（3）3种模型通过谱元法进行数值模拟计算，验证了,剖分方式的正确性。对于拐角处的单元而言，剖分方式一会引起波形散射，而剖分方式二、三均为规则波形。对于远离拐角处的单元而言，当距拐角距离＞900 m时，3种剖分方式模拟出的结果较一致；当距拐角距离为300～900 m时，剖分方式一与剖分方式二、三存在较小差异；当距拐角距离＜300 m时，剖分方式一与剖分方式二、三存在较大差异，剖分方式二、三差异始终较小。这表明删除拐角处一纵列单元方式与设置过渡环的方式均可使用，而删除四纵列拐角单元方式不推荐使用。

在起伏界面处，本文涉及的网格剖分后的剖分处理方式优势包括网格优势和程序优势。网格优势体现在：①剖分使用六面体网格，六面体网格在计算精度及存储空间等方面具有显著优势（王瑞等，2020）；②改善了规则网格剖分时利用台阶状网格近似描绘起伏界面的问题，消除了数值模拟时阶梯状边界引起的绕射和散射等问题；③剖分的网格尺寸尽可能保持均匀，避免因网格尺寸差异导致计算过程中出现发散；④剖分贴合起伏界面，提高了数值模拟的精度。程序优势体现在：基于TrueGrid软件编写了应用程序，提出了模拟计算用的起伏界面处地形六面体单元网格剖分方式，通过该程序可快速建立起伏界面处的六面体网格模型。

在今后的工作中，可将剖分方式二或三改为循环程序，通过输入起伏界面不同深度的边界坐标可得到该界面内三维精细化六面体网格剖分结果，可将程序运用到大尺度、大范围地形效应的模拟中。通过精准的网格剖分投影贴合起伏界面，以达到较高的计算精度，从而解决用六面体网格离散起伏界面模型中产生的阶梯状界面，提高谱元法处理起伏界面问题的灵活性。

致谢 国家超级计算天津中心提供了超算计算平台，在此表示感谢。