引言

复杂场地地震波场模拟一直是地震学领域中比较活跃的研究课题（孙吉泽等，2017）。目前较为成熟的数值模拟方法主要包括有限元法（finite element method，缩写为FEM）、有限差分法（finite difference mehod，缩写为FDM）和谱元法（spectral element method，缩写为SEM）等，其中谱元法自二十世纪末被引入，该方法集合伪谱法高精度和有限元法高灵活性的优势，目前已成为大尺度复杂场地地震波场模拟的主流方法（刘少林等，2021）。

随着模拟场地范围的扩大和模拟频率的提高，整体模型的网格数量急剧上升，庞大的计算量和存储需求给数值模拟造成了一定困难。借助高性能计算手段可以解决这一问题，高性能计算主要包括中央处理器（central processing unit，缩写为CPU）并行和图形处理器（graphics processing unit，缩写为GPU）并行两种方法，其中CPU并行主要通过提高CPU并行的数量达到提高计算速度的目的（刘春成等，2019）。国内外已有相关学者应用CPU并行方法进行大尺度复杂场地的地震动SEM数值模拟（Komatitsch et al，2004；Stupazzini et al，2009；胡元鑫等，2011；Chaljub et al，2015；于彦彦，2016；刘启方，2020）。但由于CPU并行规模受计算机性能限制，并行规模越大维护成本相对越高，且当CPU核数达到一定数量时，加速度比很难进一步提高，这些因素限制了CPU并行计算在大尺度复杂场地波场数值模拟上的应用（陈曦等，2016）。近年来，GPU在通用计算领域应用十分的广泛，高端GPU显卡的浮点运算最高约达每秒万亿次，在台式机或工作站内配置数块GPU计算卡即可实现大规模并行计算（熊琛等，2016）。这种并行计算方式相较于CPU并行方法量级更轻，加速效果更为明显，而且在成本和功耗上也不需要付出额外代价（Li et al，2012）。基于上述GPU计算能力优势，Komatitsch 等（2009）将GPU并行方法引入SEM数值模拟之中，并取得较好的加速效果。然而目前GPU并行计算在SEM数值模拟中的应用尚处于基础阶段，应用于大尺度复杂场地近断层地震动模拟的研究较少，因此，基于CPU-GPU异构并行架构进行复杂场地地震波场模拟，进一步研究GPU对于数值模拟的加速效果具有十分重要的意义。

基于上述分析，本文拟通过建立基本的动力学断层模型，使用开源谱元程序SPECFEM3D分别在CPU并行和CPU-GPU异构并行系统下进行计算，比较分析CPU-GPU异构并行计算的精度和效率；然后建立北京三河—平谷地区三维速度结构模型并根据已有文献确定近断层位置，构建含动力学断层的三维地形复杂模型，利用CPU-GPU异构并行模拟1679年三河—平谷M8.0地震的近断层地震动作用下北京地区的地震动响应，以期得到近断层地震动特性以及复杂地形对近断层地震动产生的影响，可为近场复杂场地地震动估计、工程抗震设计以及利用GPU加速技术进行复杂场地近断层地震动模拟提供参考。

1.   CPU-GPU异构并行方法

1.1   GPU通用计算基本原理

目前流行的GPU通用计算平台是NVIDIA公司2006年推出一种新的计算模式即统一计算设备架构（compute unified device architecture，缩写为CUDA），其具体内部架构如图1所示，可以看出，GPU相较于CPU控制单元数量少，指令流控制能力较差，但其内计算单元的数量远多于逻辑控制单元，因此在CUDA架构中，首先构造简单的数据计算模型，避免复杂的指令流控制；然后创建大规模的线程计算模型，在模型中由CPU来控制主程序，GPU作为协处理器，由GPU处理模型中大量的并行计算部分，以达到提升计算速度的目的（Shi et al，2012），并且由于GPU与主机通过卡槽相连，数据传输带宽高于传统的通信网络（Campa et al，2014）。CPU-GPU异构并行系统搭建方式如图2所示，通过CUDA并行计算架构将GPU的计算任务映射为大量的可以并行执行的线程，并由硬件执行这些线程。内核函数（kernel）以线程网格（grid）的形式组织，每个线程网格由多个线程块（block）组成，每个线程块由多个线程（thread）组成（Machidon et al，2020）。

图  1  CPU （a）和GPU （b）内部结构

Figure  1.  Internal structure of CPU （a） and GPU （b）

下载: 全尺寸图片 幻灯片

图  2  CPU-GPU异构并行架构

Figure  2.  CPU-GPU heterogeneous parallel architecture

下载: 全尺寸图片 幻灯片

1.2   SEM在CPU-GPU异构并行系统中的实现

在SEM数值模拟中，时间递进采用串行计算的方式，通常基于显式的二阶有限差分时间方法。在串行时间递进中进行大量迭代时所用的时间决定整个模拟过程的总时间，而在SEM模拟中其它部分包括预处理和后处理阶段，其时间成本可以忽略不计，故本文将利用GPU可并行计算大规模线程的优点，尽量缩短每个时间步内迭代的时间，有效地提高计算效率。并且由于数值计算中选用了静态网格划分方法和显示计算方法，导致时间递进中的每次迭代在内存和计算时间上均保持相同的成本，更有利于实现GPU并行系统的加速效果。

在CPU-GPU异构并行系统内，首先根据CPU核心数将三维模型分解为相应个数的子域，各个子域通过消息传递接口（message passing interface，缩写为MPI）在每个时间步后进行各个子域之间的数据通讯，在每个时间步内采用CPU-GPU异构并行计算模式，由CPU负责逻辑运算和数据传输，GPU负责执行大量线程的并行计算，计算完成后将数据传输至CPU端。在单个时间步内GPU执行的计算可以概括为以下三个步骤，其对应于图3中的三个CUDA内核函数。

图  3  SEM中CPU-GPU异构并行流程

Figure  3.  CPU-GPU heterogeneous parallel process in SEM

下载: 全尺寸图片 幻灯片

1） 基于前一时间步计算的加速度a，采用纽马克（Newmark）时间方法对整体模型中所有网格点的全局位移u和速度v进行计算更新。

2） 首先将全局的位移结果映射到局部单元上，在每一个网格单元上的5³个高斯-洛巴托-勒让德（Gauss-Lobatto-Legendre，缩写为GLL）积分节点中进行数组计算，得到GLL节点的应力应变值，此时的单元之间应是相互独立的；然后将计算得到的值进行装配，相邻的单元之间共享的网格节点进行叠加；最后将局部映射到全局编号，再次计算加速度a。

3） 基于第二步计算的加速度a对全局坐标下所有GLL节点的速度v进行更新。由于第二步计算中加速度结果发生了改变，因此无法进行合并只能在第三步中选择计算更新。

2.   基本模型试算

本节通过南加州地震中心（Southern California Earthquake Center，缩写为SECE）和美国地质调查局（SUnited States Geological Survey，缩写为USGS）自发破裂程序验证项目（spontaneous rupture code verification project）提供的基准断层模型TPV15，分别在CPU并行与CPU-GPU异构并行系统下进行计算，通过对比同一观测点的位移和速度结果，测试在异构并行系统下计算结果的正确性和有效性，并在此基础上通过改变TPV15模型的网格数量，对比工作站内单块CPU共36核并行、两块CPU共72核并行以及CPU-GPU异构并行系统下的计算时间，进行并行计算效率分析。工作站内CPU和GPU型号说明列于表1。

表  1  CPU和GPU硬件参数

Table  1.  CPU and GPU hardware parameters

CPU型号 Xeon Gold 6240	主频/GHz	内存容量/GB	核数
	2.6	256	36
GPU型号 GeForce RTX2080 TI	显存容量/GB	显存带宽/（GB·s−1）	计算能力	流处理器单元
	11	616	7.5	4 352

下载: 导出CSV 

| 显示表格

TPV15断层模型的在三维空间的示意如图4所示。图中断层分为主断层和分支断层，分支断层与主断层的夹角为30°，其中主断层可破裂区域为28 km×15 km，分支断层可破裂的区域为15 km×12 km，其余部分为不可破裂区域，成核区位于主断层面上，大小为3 km×3 km。对整体模型进行网格划分，其中断层周围网格大小为300 m×300 m，模型最外部网格放大为900 m×900 m，网格总数量为1 748 400，时间步长设置为0.001 s，模拟总时长12 s。TPV15模型采用均匀介质，参数为ρ＝2 670 kg/m³，v_P＝6 000 m/s，v_S＝6 000 m/s，选取基于滑动弱化摩擦准则的动力学断层破裂方式，动摩擦系数μ_s＝0.677，静摩擦系数μ_d＝0.525， 滑动弱化距离D_c＝0.4 m，断层面上的初始正应力σⁱⁿⁱ＝120.0 MPa，主断层面的剪应力$\tau _{\rm{m}}^{{\rm{ini}}} $＝−70.0 MPa；分支断层面的剪应力$\tau _{\rm{b}}^{{\rm{ini}}} $＝−78.0 MPa，成核区的剪应力$\tau _{\rm{n}}^{{\rm{ini}}} $＝−81.6 MPa。

图  4  TPV15断层模型

Figure  4.  TPV15 fault model

下载: 全尺寸图片 幻灯片

选取与SCEC/USGS自发破裂程序验证项目TPV15算例中一个相同观测点，位置如图5所示（沿走向方向2 km，沿垂直方向−0.6 km，沿深度方向0 km），基于CPU-GPU异构并行系统下计算得出该观测点的位移和速度三分量结果，与提供的结果进行对比，如图6所示，可见两曲线在波形和幅值上基本吻合，验证在CPU-GPU异构并行系统下模拟动力学断层震源具有较高的计算精度。

图  5  观测点位置示意图

Figure  5.  The diagram of observation point position

下载: 全尺寸图片 幻灯片

图  6  CPU-GPU异构并行系统计算的观测点位移、速度三分量与SECE/USGS提供结果对比

Figure  6.  The results of displacement and velocity in three directions at the observation point compared between CPU-GPU heterogeneous parallelism and SECE/USGS

下载: 全尺寸图片 幻灯片

同时为了测试CPU-GPU异构并行系统的计算性能和计算效率，通过改变TPV15模型的网格尺寸来改变网格的单元数量并分别应用CPU36核、72核和CPU-GPU异构并行进行计算，对比三者完成计算的时间（图7），发现应用GPU参与并行计算取得的加速效果十分明显。相较于CPU36核的计算时间，GPU的加速比最高为其3.04倍；与CPU72核的计算时间相比加速比最高可达2.16倍。综上表明CPU-GPU异构并行系统在进行大尺度复杂区域地震波场模拟计算方面具有明显的优势。

图  7  不同网格数量计算时间对比

Figure  7.  Calculation time between different grid numbers

下载: 全尺寸图片 幻灯片

3.   示意性算例

1679年发生的三河—平谷地震是北京周边发生的有资料记载的最大地震，推断震级为M8.0，震中距北京市中心约50 km （刘培玄等，2019）。据史料记载，此次地震由于震级较大、波及范围广、地震余震的持续时间长，造成了大量的人员伤亡和财产损失（朱耿尚，2014）。然而由于很多历史大震特别是古地震缺乏地震记录，难以进行震源破裂过程的反演，只能采用实际地质考察和数值模拟结合的方法，大致复现出古地震发生时周围区域的地震动响应。在考察分析相关资料的基础上，研究人员对1679年三河—平谷地震的震源参数进行了假定，并基于不同的数值方法对三河—平谷地震中不同的物理问题进行较为详细的研究，对深入理解震源模型的设置以及强地面运动模拟都起到一定的指导意义（高孟潭等，2002；刘博研等，2007；潘波等，2009；付长华，2012；张文强，2020）。本文将参考上述研究成果，在CPU-GPU异构并行系统下应用动力学断层震源模拟在三河—平谷地区发生M8.0地震时北京地区强地面运动，研究近断层地震动特征及实际地形对近断层地震动产生的影响。本次模拟计算区域如图8所示，其中模型经纬度范围为（116.0°—117.5°E，39.6°—40.5°N），根据目标区域得到计算模型沿东西长约为130 km，南北宽约为120 km，纵向高度约为40 km。借鉴付长华等（2015）对三河—平谷地震介质模型、断层面大小以及位置，其中断层面应力以及摩擦参数参考TPV15模型进行设置。为保证计算结果的精确可靠，要求SEM中最短波长中至少包含5个网格点，本文设定能模拟的最大频率为1 Hz，网格总数约为270万，GLL节点个数高达1.94亿，时间步距取为0.005 s。模拟90 s内的地震波传播，分别基于CPU36核、CPU72核、CPU-GPU异构并行系统进行计算，计算时间分别为25.3 h，14.6 h和9.5 h，表明CPU-GPU异构系统计算加速效果显著。三河—平谷地震发生时地表在不同时刻的地震波波场三分量快照如图9所示。

图  8  北京地区断层投影位置及整体三维物理模型

Figure  8.  Fault projection location and 3-D physical model in Beijing

下载: 全尺寸图片 幻灯片

图  9  不同周期下三维速度波场快照

（a） 东西方向；（b） 南北方向；（c） 竖直方向

Figure  9.  The snapshots of the three-dimensional velocity wave field under different cycles

（a） EW direction；（b） NS direction；（c） Vertical direction

下载: 全尺寸图片 幻灯片

为研究近断层地震动对周围重点区域的影响，本文在三河—平谷断层周边选取8个观测点，分别给出水平（EW）方向和竖直方向的加速度、速度和位移结果（图10），并将各个观测点对应的震中距和峰值地震动列于表2，以进一步研究近断层地震动特性，同时观察近断层破裂时对周边区域造成的影响。

图  10  各观测点水平 （左） 与竖直 （右）方向的加速度（a）、速度（b）和位移（c）时程

Figure  10.  Acceleration （a）， velocity （b） and displacement （c） time histories corresponding to each observation points in horizontal （left） and vertical （right） directions

下载: 全尺寸图片 幻灯片

表  2  不同观测点对应的震中距以及地震动峰值

Table  2.  Different observation points of epicenter distance and the peak value of ground motion

观测点	震中距/km	PGA/g	PGV/m·s−1	PGD/m
平谷	8.149	0.286	1.535	−1.216
三河	16.566	0.233	1.121	−0.945
通县	19.201	0.151	0.814	−0.497
兴隆	34.051	0.137	1.015	−0.753
北京	60.783	0.098	0.573	−0.413
大兴	70.293	0.043	0.307	0.196
怀柔	71.105	0.031	0.201	0.083
昌平	74.179	0.015	0.151	−0.061

下载: 导出CSV 

| 显示表格

由图10可见，此次模拟很好地反映了近断层效应：加速度峰值主要集中在沿发震断层两侧一个狭窄的范围内，其中平谷地区的峰值加速度最大达到0.286g，并且随着断层距的增加，加速度峰值迅速衰减，昌平地区的峰值加速度只达到0.015g，这一近断层特征被总结为近断层的集中效应，该特征在一些大地震观测记录中也得到了实证（如1970年通海地震、1995年日本阪神地震、1999年集集地震、2008年汶川地震和2018年精河地震等）（刘启方等，2006）；近断层速度大脉冲现象显著，对比图10b可以看出，靠近断层位置的三河、平谷等地区出现了速度峰值大、持时短的速度脉冲现象，平谷地区水平方向的峰值速度为1.535 m/s，三河地区水平方向的峰值速度为1.121 m/s，同时观察到由于断层设置为右旋的垂直走滑断层，沿水平方向的速度时程结果均高于竖直方向的结果，体现了近断层的方向性效应；近断层效应中比较重要的特性为地面的永久变形，该现象在靠近断层的位置相对明显，远离近断层的区域则永久位移现象并不明显，以平谷地区为例，模拟中水平最大位移达到1.216 m，永久位移达到0.63 m。

结果中除近断层地震动效应比较明显之外，地形效应对近断层地震动影响同样显著。以兴隆观测点为例，其地震动响应明显高于平原地区的位置水平峰值加速度达到0.137g，水平峰值速度达到1.015 m/s，同时三河—平谷地区地震动持时为25 s，通州地区地震动持时为30 s，对比可以看出兴隆地区地震动持时明显延长，地震动持续时间接近45 s。出现上述现象的主要是由于兴隆地区位于燕山山脉中段，地震波传播至山体地形时出现了明显的能量聚集效应（郝明辉，张郁山，2019），导致其地震动峰值增大，持时明显延长。

4.   讨论与结论

本文为拓展单台工作站模拟大尺度复杂场地地震动的能力，有效提高谱元法模拟的计算效率，基于CUDA编程平台实现了工作站级的CPU-GPU异构并行方法，并利用该异构并行方法分别进行了SECE/USGS提供的自发破裂模型TPV15和1679年三河—平谷M8.0强地面运动模拟，得到主要结论如下：

1） 根据SECE/USGS提供的自发破裂模型TPV15以及计算结果，对比同一观测点的速度、位移三分量结果，验证了CPU-GPU异构并行系统下谱元法模拟动力学断层的精度。改变网格数量对比CPU36核、CPU72核以及CPU-GPU异构并行的计算时间，测试结果显示加速比最高可达到CPU 36核和72核并行的3.04倍和2.16倍，表明CPU-GPU异构并行方法不仅满足计算精度要求，同时大大降低了计算时间，体现了工作站级CPU-GPU异构并行方法的计算优势。

2） 本文模拟了1679年三河—平谷M8.0地震中北京地区强地面运动的取得良好的加速效果，模拟结果体现了包括集中性、方向性、地面永久位移等近断层地震动的基本特征，以及地形效应对近断层地震动产生的显著影响，验证了所提出的异构并行方法对真实设定地震动模拟的适用性。在未来可进一步增加GPU并行数量，在工作站内基于CPU-GPU异构并行系统拓展到更精细化和范围更广的复杂场地地震动模拟，为地震区划、抗震设防等工作提供参考依据。