巴振寧 趙靖軒 吳孟桃 梁建文1，

1） 中國天津 300354 中國地震局地震工程綜合模擬與城鄉(xiāng)抗震韌性重點(diǎn)實(shí)驗(yàn)室

2） 中國天津 300354 天津大學(xué)土木工程系

引言

復(fù)雜場地地震波場模擬一直是地震學(xué)領(lǐng)域中比較活躍的研究課題（孫吉澤等，2017）.目前較為成熟的數(shù)值模擬方法主要包括有限元法（finite element method，縮寫為FEM）、有限差分法（finite difference mehod，縮寫為FDM）和譜元法（spectral element method，縮寫為SEM）等，其中譜元法自二十世紀(jì)末被引入，該方法集合偽譜法高精度和有限元法高靈活性的優(yōu)勢，目前已成為大尺度復(fù)雜場地地震波場模擬的主流方法（劉少林等，2021）.

隨著模擬場地范圍的擴(kuò)大和模擬頻率的提高，整體模型的網(wǎng)格數(shù)量急劇上升，龐大的計算量和存儲需求給數(shù)值模擬造成了一定困難.借助高性能計算手段可以解決這一問題，高性能計算主要包括中央處理器（central processing unit，縮寫為CPU）并行和圖形處理器（graphics processing unit，縮寫為GPU）并行兩種方法，其中CPU 并行主要通過提高CPU 并行的數(shù)量達(dá)到提高計算速度的目的（劉春成等，2019）.國內(nèi)外已有相關(guān)學(xué)者應(yīng)用CPU 并行方法進(jìn)行大尺度復(fù)雜場地的地震動SEM 數(shù)值模擬（Komatitschet al，2004；Stupazziniet al，2009；胡元鑫等，2011；Chaljubet al，2015；于彥彥，2016；劉啟方，2020）.但由于CPU 并行規(guī)模受計算機(jī)性能限制，并行規(guī)模越大維護(hù)成本相對越高，且當(dāng)CPU 核數(shù)達(dá)到一定數(shù)量時，加速度比很難進(jìn)一步提高，這些因素限制了CPU 并行計算在大尺度復(fù)雜場地波場數(shù)值模擬上的應(yīng)用（陳曦等，2016）.近年來，GPU 在通用計算領(lǐng)域應(yīng)用十分的廣泛，高端GPU 顯卡的浮點(diǎn)運(yùn)算最高約達(dá)每秒萬億次，在臺式機(jī)或工作站內(nèi)配置數(shù)塊GPU 計算卡即可實(shí)現(xiàn)大規(guī)模并行計算（熊琛等，2016）.這種并行計算方式相較于CPU 并行方法量級更輕，加速效果更為明顯，而且在成本和功耗上也不需要付出額外代價（Liet al，2012）.基于上述GPU 計算能力優(yōu)勢，Komatitsch 等（2009）將GPU 并行方法引入SEM 數(shù)值模擬之中，并取得較好的加速效果.然而目前GPU 并行計算在SEM 數(shù)值模擬中的應(yīng)用尚處于基礎(chǔ)階段，應(yīng)用于大尺度復(fù)雜場地近斷層地震動模擬的研究較少，因此，基于CPU-GPU 異構(gòu)并行架構(gòu)進(jìn)行復(fù)雜場地地震波場模擬，進(jìn)一步研究GPU 對于數(shù)值模擬的加速效果具有十分重要的意義.

基于上述分析，本文擬通過建立基本的動力學(xué)斷層模型，使用開源譜元程序SPECFEM3D分別在CPU 并行和CPU-GPU 異構(gòu)并行系統(tǒng)下進(jìn)行計算，比較分析CPU-GPU 異構(gòu)并行計算的精度和效率；然后建立北京三河—平谷地區(qū)三維速度結(jié)構(gòu)模型并根據(jù)已有文獻(xiàn)確定近斷層位置，構(gòu)建含動力學(xué)斷層的三維地形復(fù)雜模型，利用CPU-GPU 異構(gòu)并行模擬1679 年三河—平谷M8.0 地震的近斷層地震動作用下北京地區(qū)的地震動響應(yīng)，以期得到近斷層地震動特性以及復(fù)雜地形對近斷層地震動產(chǎn)生的影響，可為近場復(fù)雜場地地震動估計、工程抗震設(shè)計以及利用GPU 加速技術(shù)進(jìn)行復(fù)雜場地近斷層地震動模擬提供參考.

1 CPU-GPU 異構(gòu)并行方法

1.1 GPU 通用計算基本原理

目前流行的GPU 通用計算平臺是NVIDIA 公司2006 年推出一種新的計算模式即統(tǒng)一計算設(shè)備架構(gòu)（compute unified device architecture，縮寫為CUDA），其具體內(nèi)部架構(gòu)如圖1 所示，可以看出，GPU 相較于CPU 控制單元數(shù)量少，指令流控制能力較差，但其內(nèi)計算單元的數(shù)量遠(yuǎn)多于邏輯控制單元，因此在CUDA 架構(gòu)中，首先構(gòu)造簡單的數(shù)據(jù)計算模型，避免復(fù)雜的指令流控制；然后創(chuàng)建大規(guī)模的線程計算模型，在模型中由CPU 來控制主程序，GPU 作為協(xié)處理器，由GPU 處理模型中大量的并行計算部分，以達(dá)到提升計算速度的目的（Shiet al，2012），并且由于GPU 與主機(jī)通過卡槽相連，數(shù)據(jù)傳輸帶寬高于傳統(tǒng)的通信網(wǎng)絡(luò)（Campaet al，2014）.CPU-GPU 異構(gòu)并行系統(tǒng)搭建方式如圖2 所示，通過CUDA 并行計算架構(gòu)將GPU 的計算任務(wù)映射為大量的可以并行執(zhí)行的線程，并由硬件執(zhí)行這些線程.內(nèi)核函數(shù)（kernel）以線程網(wǎng)格（grid）的形式組織，每個線程網(wǎng)格由多個線程塊（block）組成，每個線程塊由多個線程（thread）組成（Machidonet al，2020）.

圖1 CPU （a）和GPU （b）內(nèi)部結(jié)構(gòu)Fig. 1 Internal structure of CPU （a） and GPU （b）

圖2 CPU-GPU 異構(gòu)并行架構(gòu)Fig. 2 CPU-GPU heterogeneous parallel architecture

1.2 SEM 在CPU-GPU 異構(gòu)并行系統(tǒng)中的實(shí)現(xiàn)

在SEM 數(shù)值模擬中，時間遞進(jìn)采用串行計算的方式，通?；陲@式的二階有限差分時間方法.在串行時間遞進(jìn)中進(jìn)行大量迭代時所用的時間決定整個模擬過程的總時間，而在SEM 模擬中其它部分包括預(yù)處理和后處理階段，其時間成本可以忽略不計，故本文將利用GPU 可并行計算大規(guī)模線程的優(yōu)點(diǎn)，盡量縮短每個時間步內(nèi)迭代的時間，有效地提高計算效率.并且由于數(shù)值計算中選用了靜態(tài)網(wǎng)格劃分方法和顯示計算方法，導(dǎo)致時間遞進(jìn)中的每次迭代在內(nèi)存和計算時間上均保持相同的成本，更有利于實(shí)現(xiàn)GPU 并行系統(tǒng)的加速效果.

在CPU-GPU 異構(gòu)并行系統(tǒng)內(nèi)，首先根據(jù)CPU 核心數(shù)將三維模型分解為相應(yīng)個數(shù)的子域，各個子域通過消息傳遞接口（message passing interface，縮寫為MPI）在每個時間步后進(jìn)行各個子域之間的數(shù)據(jù)通訊，在每個時間步內(nèi)采用CPU-GPU 異構(gòu)并行計算模式，由CPU 負(fù)責(zé)邏輯運(yùn)算和數(shù)據(jù)傳輸，GPU 負(fù)責(zé)執(zhí)行大量線程的并行計算，計算完成后將數(shù)據(jù)傳輸至CPU 端.在單個時間步內(nèi)GPU 執(zhí)行的計算可以概括為以下三個步驟，其對應(yīng)于圖3 中的三個CUDA 內(nèi)核函數(shù).

圖3 SEM 中CPU-GPU 異構(gòu)并行流程Fig. 3 CPU-GPU heterogeneous parallel process in SEM

1） 基于前一時間步計算的加速度a，采用紐馬克（Newmark）時間方法對整體模型中所有網(wǎng)格點(diǎn)的全局位移u和速度v進(jìn)行計算更新.

2） 首先將全局的位移結(jié)果映射到局部單元上，在每一個網(wǎng)格單元上的53個高斯-洛巴托-勒讓德（Gauss-Lobatto-Legendre，縮寫為GLL）積分節(jié)點(diǎn)中進(jìn)行數(shù)組計算，得到GLL 節(jié)點(diǎn)的應(yīng)力應(yīng)變值，此時的單元之間應(yīng)是相互獨(dú)立的；然后將計算得到的值進(jìn)行裝配，相鄰的單元之間共享的網(wǎng)格節(jié)點(diǎn)進(jìn)行疊加；最后將局部映射到全局編號，再次計算加速度a.

3） 基于第二步計算的加速度a對全局坐標(biāo)下所有GLL 節(jié)點(diǎn)的速度v進(jìn)行更新.由于第二步計算中加速度結(jié)果發(fā)生了改變，因此無法進(jìn)行合并只能在第三步中選擇計算更新.

2 基本模型試算

本節(jié)通過南加州地震中心（Southern California Earthquake Center，縮寫為SECE）和美國地質(zhì)調(diào)查局（SUnited States Geological Survey，縮寫為USGS）自發(fā)破裂程序驗(yàn)證項(xiàng)目（spontaneous rupture code verification project）提供的基準(zhǔn)斷層模型TPV15，分別在CPU 并行與CPU-GPU 異構(gòu)并行系統(tǒng)下進(jìn)行計算，通過對比同一觀測點(diǎn)的位移和速度結(jié)果，測試在異構(gòu)并行系統(tǒng)下計算結(jié)果的正確性和有效性，并在此基礎(chǔ)上通過改變TPV15 模型的網(wǎng)格數(shù)量，對比工作站內(nèi)單塊CPU 共36 核并行、兩塊CPU 共72 核并行以及CPU-GPU 異構(gòu)并行系統(tǒng)下的計算時間，進(jìn)行并行計算效率分析.工作站內(nèi)CPU 和GPU 型號說明列于表1.

表1 CPU 和GPU 硬件參數(shù)Table 1 CPU and GPU hardware parameters

圖4 TPV15 斷層模型Fig. 4 TPV15 fault model

選取與SCEC/USGS 自發(fā)破裂程序驗(yàn)證項(xiàng)目TPV15 算例中一個相同觀測點(diǎn)，位置如圖5 所示（沿走向方向2 km，沿垂直方向?0.6 km，沿深度方向0 km），基于CPUGPU 異構(gòu)并行系統(tǒng)下計算得出該觀測點(diǎn)的位移和速度三分量結(jié)果，與提供的結(jié)果進(jìn)行對比，如圖6 所示，可見兩曲線在波形和幅值上基本吻合，驗(yàn)證在CPU-GPU異構(gòu)并行系統(tǒng)下模擬動力學(xué)斷層震源具有較高的計算精度.

圖5 觀測點(diǎn)位置示意圖Fig. 5 The diagram of observation point position

圖6 CPU-GPU 異構(gòu)并行系統(tǒng)計算的觀測點(diǎn)位移、速度三分量與SECE/USGS 提供結(jié)果對比Fig. 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 異構(gòu)并行系統(tǒng)的計算性能和計算效率，通過改變TPV15模型的網(wǎng)格尺寸來改變網(wǎng)格的單元數(shù)量并分別應(yīng)用CPU36 核、72 核和CPU-GPU 異構(gòu)并行進(jìn)行計算，對比三者完成計算的時間（圖7），發(fā)現(xiàn)應(yīng)用GPU 參與并行計算取得的加速效果十分明顯.相較于CPU36 核的計算時間，GPU 的加速比最高為其3.04 倍；與CPU72 核的計算時間相比加速比最高可達(dá)2.16 倍.綜上表明CPU-GPU 異構(gòu)并行系統(tǒng)在進(jìn)行大尺度復(fù)雜區(qū)域地震波場模擬計算方面具有明顯的優(yōu)勢.

圖7 不同網(wǎng)格數(shù)量計算時間對比Fig. 7 Calculation time between different grid numbers

3 示意性算例

1679 年發(fā)生的三河—平谷地震是北京周邊發(fā)生的有資料記載的最大地震，推斷震級為M8.0，震中距北京市中心約50 km （劉培玄等，2019）.據(jù)史料記載，此次地震由于震級較大、波及范圍廣、地震余震的持續(xù)時間長，造成了大量的人員傷亡和財產(chǎn)損失（朱耿尚，2014）.然而由于很多歷史大震特別是古地震缺乏地震記錄，難以進(jìn)行震源破裂過程的反演，只能采用實(shí)際地質(zhì)考察和數(shù)值模擬結(jié)合的方法，大致復(fù)現(xiàn)出古地震發(fā)生時周圍區(qū)域的地震動響應(yīng).在考察分析相關(guān)資料的基礎(chǔ)上，研究人員對1679 年三河—平谷地震的震源參數(shù)進(jìn)行了假定，并基于不同的數(shù)值方法對三河—平谷地震中不同的物理問題進(jìn)行較為詳細(xì)的研究，對深入理解震源模型的設(shè)置以及強(qiáng)地面運(yùn)動模擬都起到一定的指導(dǎo)意義（高孟潭等，2002；劉博研等，2007；潘波等，2009；付長華，2012；張文強(qiáng)，2020）.本文將參考上述研究成果，在CPU-GPU 異構(gòu)并行系統(tǒng)下應(yīng)用動力學(xué)斷層震源模擬在三河—平谷地區(qū)發(fā)生M8.0 地震時北京地區(qū)強(qiáng)地面運(yùn)動，研究近斷層地震動特征及實(shí)際地形對近斷層地震動產(chǎn)生的影響.本次模擬計算區(qū)域如圖8 所示，其中模型經(jīng)緯度范圍為（116.0°—117.5°E，39.6°—40.5°N），根據(jù)目標(biāo)區(qū)域得到計算模型沿東西長約為130 km，南北寬約為120 km，縱向高度約為40 km.借鑒付長華等（2015）對三河—平谷地震介質(zhì)模型、斷層面大小以及位置，其中斷層面應(yīng)力以及摩擦參數(shù)參考TPV15 模型進(jìn)行設(shè)置.為保證計算結(jié)果的精確可靠，要求SEM 中最短波長中至少包含5 個網(wǎng)格點(diǎn)，本文設(shè)定能模擬的最大頻率為1 Hz，網(wǎng)格總數(shù)約為270 萬，GLL 節(jié)點(diǎn)個數(shù)高達(dá)1.94 億，時間步距取為0.005 s.模擬90 s 內(nèi)的地震波傳播，分別基于CPU36 核、CPU72 核、CPU-GPU 異構(gòu)并行系統(tǒng)進(jìn)行計算，計算時間分別為25.3 h，14.6 h 和9.5 h，表明CPU-GPU 異構(gòu)系統(tǒng)計算加速效果顯著.三河—平谷地震發(fā)生時地表在不同時刻的地震波波場三分量快照如圖9 所示.

圖8 北京地區(qū)斷層投影位置及整體三維物理模型Fig. 8 Fault projection location and 3-D physical model in Beijing

圖9 不同周期下三維速度波場快照（a） 東西方向；（b） 南北方向；（c） 豎直方向Fig. 9 The snapshots of the three-dimensional velocity wave field under different cycles（a） EW direction；（b） NS direction；（c） Vertical direction

為研究近斷層地震動對周圍重點(diǎn)區(qū)域的影響，本文在三河—平谷斷層周邊選取8 個觀測點(diǎn)，分別給出水平（EW）方向和豎直方向的加速度、速度和位移結(jié)果（圖10），并將各個觀測點(diǎn)對應(yīng)的震中距和峰值地震動列于表2，以進(jìn)一步研究近斷層地震動特性，同時觀察近斷層破裂時對周邊區(qū)域造成的影響.

表2 不同觀測點(diǎn)對應(yīng)的震中距以及地震動峰值Table 2 Different observation points of epicenter distance and the peak value of ground motion

圖10 各觀測點(diǎn)水平 （左） 與豎直 （右）方向的加速度（a）、速度（b）和位移（c）時程Fig. 10 Acceleration （a）， velocity （b） and displacement （c） time histories corresponding to each observation points in horizontal （left） and vertical （right） directions

由圖10 可見，此次模擬很好地反映了近斷層效應(yīng)：加速度峰值主要集中在沿發(fā)震斷層兩側(cè)一個狹窄的范圍內(nèi)，其中平谷地區(qū)的峰值加速度最大達(dá)到0.286g，并且隨著斷層距的增加，加速度峰值迅速衰減，昌平地區(qū)的峰值加速度只達(dá)到0.015g，這一近斷層特征被總結(jié)為近斷層的集中效應(yīng)，該特征在一些大地震觀測記錄中也得到了實(shí)證（如1970 年通海地震、1995 年日本阪神地震、1999 年集集地震、2008 年汶川地震和2018 年精河地震等）（劉啟方等，2006）；近斷層速度大脈沖現(xiàn)象顯著，對比圖10b 可以看出，靠近斷層位置的三河、平谷等地區(qū)出現(xiàn)了速度峰值大、持時短的速度脈沖現(xiàn)象，平谷地區(qū)水平方向的峰值速度為1.535 m/s，三河地區(qū)水平方向的峰值速度為1.121 m/s，同時觀察到由于斷層設(shè)置為右旋的垂直走滑斷層，沿水平方向的速度時程結(jié)果均高于豎直方向的結(jié)果，體現(xiàn)了近斷層的方向性效應(yīng)；近斷層效應(yīng)中比較重要的特性為地面的永久變形，該現(xiàn)象在靠近斷層的位置相對明顯，遠(yuǎn)離近斷層的區(qū)域則永久位移現(xiàn)象并不明顯，以平谷地區(qū)為例，模擬中水平最大位移達(dá)到1.216 m，永久位移達(dá)到0.63 m.

結(jié)果中除近斷層地震動效應(yīng)比較明顯之外，地形效應(yīng)對近斷層地震動影響同樣顯著.以興隆觀測點(diǎn)為例，其地震動響應(yīng)明顯高于平原地區(qū)的位置水平峰值加速度達(dá)到0.137g，水平峰值速度達(dá)到1.015 m/s，同時三河—平谷地區(qū)地震動持時為25 s，通州地區(qū)地震動持時為30 s，對比可以看出興隆地區(qū)地震動持時明顯延長，地震動持續(xù)時間接近45 s.出現(xiàn)上述現(xiàn)象的主要是由于興隆地區(qū)位于燕山山脈中段，地震波傳播至山體地形時出現(xiàn)了明顯的能量聚集效應(yīng)（郝明輝，張郁山，2019），導(dǎo)致其地震動峰值增大，持時明顯延長.

4 討論與結(jié)論

本文為拓展單臺工作站模擬大尺度復(fù)雜場地地震動的能力，有效提高譜元法模擬的計算效率，基于CUDA 編程平臺實(shí)現(xiàn)了工作站級的CPU-GPU 異構(gòu)并行方法，并利用該異構(gòu)并行方法分別進(jìn)行了SECE/USGS 提供的自發(fā)破裂模型TPV15 和1679 年三河—平谷M8.0 強(qiáng)地面運(yùn)動模擬，得到主要結(jié)論如下：

1） 根據(jù)SECE/USGS 提供的自發(fā)破裂模型TPV15 以及計算結(jié)果，對比同一觀測點(diǎn)的速度、位移三分量結(jié)果，驗(yàn)證了CPU-GPU 異構(gòu)并行系統(tǒng)下譜元法模擬動力學(xué)斷層的精度.改變網(wǎng)格數(shù)量對比CPU36 核、CPU72 核以及CPU-GPU 異構(gòu)并行的計算時間，測試結(jié)果顯示加速比最高可達(dá)到CPU 36 核和72 核并行的3.04 倍和2.16 倍，表明CPU-GPU 異構(gòu)并行方法不僅滿足計算精度要求，同時大大降低了計算時間，體現(xiàn)了工作站級CPU-GPU 異構(gòu)并行方法的計算優(yōu)勢.

2） 本文模擬了1679 年三河—平谷M8.0 地震中北京地區(qū)強(qiáng)地面運(yùn)動的取得良好的加速效果，模擬結(jié)果體現(xiàn)了包括集中性、方向性、地面永久位移等近斷層地震動的基本特征，以及地形效應(yīng)對近斷層地震動產(chǎn)生的顯著影響，驗(yàn)證了所提出的異構(gòu)并行方法對真實(shí)設(shè)定地震動模擬的適用性.在未來可進(jìn)一步增加GPU 并行數(shù)量，在工作站內(nèi)基于CPU-GPU 異構(gòu)并行系統(tǒng)拓展到更精細(xì)化和范圍更廣的復(fù)雜場地地震動模擬，為地震區(qū)劃、抗震設(shè)防等工作提供參考依據(jù).