巴振寧， 趙靖軒， 張郁山， 梁建文， 張玉潔

1 天津大學(xué)中國地震局地震工程綜合模擬與城鄉(xiāng)抗震韌性重點實驗室， 天津 300350 2 天津大學(xué)土木工程系， 天津 300350 3 中國地震災(zāi)害防御中心， 北京 100029

0 引言

近年來，近場強地面運動模擬逐漸成為工程地震學(xué)和理論地震學(xué)的一個熱點課題，相關(guān)研究可為強震記錄不足地區(qū)的工程結(jié)構(gòu)抗震設(shè)防以及城市防震減災(zāi)提供依據(jù)(王海云和謝禮立, 2008).隨著城市現(xiàn)代化的發(fā)展，超高層建筑、長大橋梁、超高拱壩等一些重大工程結(jié)構(gòu)不斷涌現(xiàn)，對這些具有較長自振周期的工程結(jié)構(gòu)進(jìn)行抗震分析時，需要考慮低頻地震動(<1 Hz)對結(jié)構(gòu)非線性響應(yīng)產(chǎn)生的重要影響，尤其是近斷層區(qū)可能存在的長周期速度大脈沖；而高頻地震動(>1 Hz)控制著地震動的幅值，對結(jié)構(gòu)地震反應(yīng)同樣也不可忽視(許可, 2020).因此，開展寬頻帶地震動模擬研究以為工程結(jié)構(gòu)抗震設(shè)防提供合理輸入有著重要的科學(xué)意義和工程應(yīng)用價值.

目前寬頻地震模擬的主流方法是將確定性方法模擬的低頻地震動和隨機方法模擬的高頻地震動相結(jié)合的混合方法(Graves and Pitarka, 2010; Frankel, 2009; Mai et al., 2010; 孫曉丹和陶夏新, 2012).但該方法無法反映三維地形對高頻地震動產(chǎn)生的影響(Lin et al., 2011; Baltay et al., 2017; Sahakian et al., 2019).同時，隨機方法得到的結(jié)果僅是一條方向未知的水平向加速度時程，而確定性方法能夠獲得三分量地震動，兩者的疊加具有不協(xié)調(diào)之處(曹澤林, 2020).為避免上述問題，綜合考慮震源-傳播-場地整個物理機制對地震波場特性影響，采用完全確定性方法進(jìn)行的寬頻帶地震動模擬逐漸成為一種趨勢.

確定性的寬頻地震動模擬依賴于高效的解析或數(shù)值方法以及合理的震源模型.目前，確定性的寬頻地震動模擬可采用頻率波數(shù)域(FK)和各種數(shù)值方法.其中，F(xiàn)K方法作為一種半解析方法，適用于水平成層的一維波速結(jié)構(gòu)模型內(nèi)源引起的響應(yīng)問題，由于無需劃分網(wǎng)格且自動滿足無窮遠(yuǎn)波輻射條件，模擬精度高、模擬帶寬范圍較廣(Cao et al., 2019; 梁建文等, 2022).確定性方法中有限元法(FEM)、有限差分法(FDM)、譜元法(SEM)等數(shù)值方法能夠計算三維復(fù)雜場地的地震響應(yīng)，但由于對區(qū)域地殼內(nèi)小尺度不均勻體和震源的高頻特征缺乏足夠認(rèn)識，導(dǎo)致了上述方法在過去多限于低頻模擬.但由于地震和背景噪聲記錄的不斷豐富，相關(guān)學(xué)者通過反演等手段對地下波速結(jié)構(gòu)有了更加清晰的認(rèn)識.同時，已有學(xué)者提出了能夠激發(fā)寬頻地震波的合理運動學(xué)或動力學(xué)震源模型.針對運動學(xué)震源模型，目前主要包括確定性震源、隨機震源和混合震源模型(王海云等, 2008)：確定性震源模型主要強調(diào)了斷層面上凹凸體對震源破裂的控制作用(Irikura and Miyake, 2011)；隨機震源模型主要強調(diào)震源產(chǎn)生的隨機散射和介質(zhì)不均勻?qū)е碌母哳l擾動作用；混合震源模型由于同時考慮了破裂面上凹凸體對滑動量不均勻分布的主導(dǎo)作用和滑動量的隨機性，是激發(fā)寬頻地震波的合理震源模型.

目前，隨著相關(guān)學(xué)者針對地下波速結(jié)構(gòu)和震源破裂的認(rèn)識逐漸深刻，以及超級計算的出現(xiàn)突破了計算能力和計算耗時的瓶頸，推動了各數(shù)值方法在區(qū)域尺度寬頻帶地震動模擬方面的發(fā)展.基于數(shù)值方法的模擬區(qū)域范圍從幾十公里到幾百公里，網(wǎng)格單元的總數(shù)從數(shù)百億到數(shù)千億；模擬的頻率從2～5 Hz、10 Hz或更高，已基本可以覆蓋工程結(jié)構(gòu)敏感頻段.事實上，美國勞倫斯國家重點實驗室Rodgers等(2019, 2020)最新證實了0～5和0～10 Hz的寬頻模擬可采用完全確定性物理的技術(shù)實現(xiàn).Fu等(2017) 借助神威太湖之光超級計算機，針對1976年唐山大地震實現(xiàn)了整個華北地區(qū)的寬頻(0～18 Hz)地震動模擬.上述成果表明，基于確定性的寬頻帶地震動模擬是現(xiàn)代地震工程地震動模擬的重要發(fā)展方向.

因此，本文為實現(xiàn)基于譜元法的確定性寬頻地震動模擬，主要開展了兩項工作.首先，將確定性凹凸體震源模型和Graves和Pitarka(2010, 2015)提出的GP14.3隨機震源模型結(jié)合，進(jìn)而給出了一種能激發(fā)寬頻地震波的混合震源模型；將上述運動學(xué)混合震源模型開發(fā)到了SPECFEM 3D開源代碼中，實現(xiàn)了基于GP14.3混合運動學(xué)震源模型和SEM的全過程寬頻帶地震動模擬.

1 譜元法及混合運動學(xué)震源模型

1.1 SEM理論

SEM屬于廣義有限元法，建立在波動方程的弱形式理論基礎(chǔ)上，是求解偏微分方程的一種有效的數(shù)值計算方法.該方法集整了偽譜法和FEM的優(yōu)勢，兼?zhèn)淞藗巫V法的高精度、快速收斂和FEM高靈活性的特點.并且SEM劃分的單元更稀疏，可避免高階有限元可能存在的Runge現(xiàn)象.此外，由于單元質(zhì)量陣為對角陣，因此大大減少了計算存儲，易于在多節(jié)點上實現(xiàn)并行計算(李孝波等, 2014).近年來，隨著Komatitsch和Tromp(1999)將SEM推廣到地震波場的模擬并推出適用于二維和三維問題的譜元程序SPECFEM 2D及SPECFEM 3D，SEM作為一種高精度的數(shù)值方法得到越來越廣泛的應(yīng)用.在本節(jié)中簡要介紹SEM(基于Legendre展開)的求解過程.

震源激發(fā)的地震波滿足線彈性動力學(xué)控制方程：

(1)

σ=c∶ε,

(2)

式中，c為介質(zhì)剛度張量，“∶”為張量乘法，ε為應(yīng)變張量，其表達(dá)式為

(3)

震源項f為一雙力偶點源，一般可以表示為

(4)

式中，M為對稱矩張量，其各元素表示描述震源錯動中的經(jīng)分解后的不同位態(tài).S(t)為描述震源破裂過程的震源時間函數(shù)，x和xS為任一場點和震源位置處.由于在寬頻地震動受震源破裂過程的影響，使用單個點源描述發(fā)震機制不再符合實際情況，故需要將破裂過程等效為多個子源疊加，包括各個子源的矩張量和時間函數(shù)，描述混合震源模型破裂過程的時空分布.

將式(2)—(4)代入式(1)中，結(jié)合自由表面零應(yīng)力條件可得到計算域Ω的弱形式：

(5)

式中，w為引入的任意試函數(shù)，以此模型表面零應(yīng)力條件已自動滿足，這也是采用方程弱形式區(qū)別于強形式的最明顯特征.另一方面，采用nl階六面體對三維計算域進(jìn)行離散，進(jìn)而某特定譜單元Ωe上的函數(shù)可用其上的GLL插值節(jié)點表征：

(6)

式中，fα β γ為參考點處的插值形函數(shù)權(quán)值，lk(k=α,β,γ)為k階Lagrange多項式，(ξ,η,ζ)為離散六面體單元對應(yīng)的標(biāo)準(zhǔn)參考單元上的插值節(jié)點，它與六面體上的插值點x=(x,y,z)建立了一一映射的關(guān)系，且有

dxdydz=Jedξdηdζ,

(7)

式中，Je表示Ωe與其對應(yīng)參考單元的三階Jacobian矩陣.因此式(5)中函數(shù)在譜單元上的積分可利用GLL積分法則簡化為

(8)

式中，ωk(k=α,β,γ)為GLL點上的積分權(quán)值.進(jìn)一步通過方程集整，疊加鄰接譜單元對共享節(jié)點的貢獻(xiàn)，可建立起關(guān)于整體位移向量U的常微分方程：

(9)

式中，M、C和K分別為對角質(zhì)量矩陣、吸收邊界矩陣和剛度矩陣.上述方程可采用顯式Newmark-β法高效求解：

(10)

SEM模擬波動傳播時涉及相對周期誤差(物理頻散)和數(shù)值阻尼誤差，前者是對波傳播過程中由于網(wǎng)格尺寸、積分方法和單元類型的改變導(dǎo)致對波速的錯誤估計造成.采用CFL(Courant-Friedrichs-Lewy)條件對顯式格式中的步長施加約束：Δt≤γΔxmin/vmax，其中γ為嚴(yán)格小于1的常數(shù)，Δxmin和vmax為模型GLL節(jié)點間最小距離和最大波速，同時網(wǎng)格劃分尺寸應(yīng)滿足條件：s≤λminN/5，其中λmin為介質(zhì)中傳播最短波長，N為Lagrange多項式的階數(shù)，取4≤N≤8，一般取N=4.

1.2 基于GP14.3的混合震源在SEM中的開發(fā)

混合震源模型主要是將低波數(shù)確定性的凹凸體震源模型中引入高頻隨機成分，以確保合成寬頻地面運動的帶寬有效性.首先根據(jù)定標(biāo)率建立包含凹凸體震源的有限斷層模型，確定性的凹凸體震源建立過程不再過多闡述，具體細(xì)節(jié)可以參考王海云(2004).進(jìn)而在凹凸體震源分布的基礎(chǔ)上考慮在斷層面滑動分布中結(jié)合Graves和Pitarka(2010, 2015)提出的GP14.3隨機震源模型，該隨機震源模型采用半隨機方法，即調(diào)整隨機錯動分布使其具有k-2波數(shù)譜，是一種經(jīng)過多位學(xué)者驗證過，并且成功應(yīng)用到舊金山地區(qū)海沃德斷層的寬頻帶強地面模擬中較為成熟的隨機震源模型.本節(jié)將主要介紹GP14.3隨機震源公式的推導(dǎo)以及將混合震源開發(fā)到譜元程序SPECFEM 3D中的流程.

首先將已有的凹凸體斷層面上滑動量的空間分布經(jīng)二維傅里葉變換至波數(shù)域，得到斷層面上確定性的滑動波數(shù)譜Dslip(ks,kd):

Dslip(ks,kd)=?Uslip(x′,y′)e-iksxe-ikdydksdkd,

(11)

式中，ks、kd分別表示沿斷層面走向、傾向的波數(shù)，Uslip(x′,y′)表示斷層面上滑動量的空間分布，坐標(biāo)(x′,y′)位于斷層面上的局部坐標(biāo)系中.

其次采用波數(shù)衰減滿足von Karman自相關(guān)函數(shù)的波數(shù)譜并引入隨機數(shù)φ表達(dá)震源中的隨機成分(Mai and Beroza, 2002)：

(12)

式中，φ是(-π, π)區(qū)間上服從均勻分布的隨機數(shù)，H是Hurst指數(shù)，這里根據(jù)Graves和Pitarka(2015)設(shè)置為0.75，as和ad設(shè)置為與震級MW相關(guān)的兩個經(jīng)驗系數(shù)

(13)

然后在波數(shù)域中結(jié)合確定性的低波數(shù)譜和隨機高波數(shù)譜，并利用二維逆傅里葉變換至空間域中：

×Dslip(0,0)(1-W)]eiksxeikdydksdkd,

(14)

式中，Dslip(0,0)表示走向、傾向為0時的滑動波數(shù)譜，W表示波數(shù)結(jié)合函數(shù)，用以選擇確定性部分的低波數(shù)成分和隨機部分的高波數(shù)成分，其表達(dá)式為

(15)

式中，N控制結(jié)合的銳度，設(shè)置為4；kcs和kcd表示沿斷層走向和傾向的拐角波數(shù)，用以確定震源譜中低波數(shù)和高波數(shù)的結(jié)合界限，設(shè)置kcs=1/斷層沿走向長度，kcd=1/斷層沿傾向?qū)挾?由此，斷層面上的滑動量分布得以確定.

最后基于斷層面上的滑動量可以得到每個斷層上的M0，進(jìn)而根據(jù)斷層的走向、傾向和滑動角得到每個子源對應(yīng)的地震矩張量.

關(guān)于破裂起始時刻的分布(即斷層面上的破裂傳播過程)，首先需要確定破裂速度的初始分布.考慮淺地殼層區(qū)域的破裂能力相較于深層弱(Kagawa et al., 2004)，將地表以下5 km深度范圍內(nèi)的平均破裂速度設(shè)置為剪切波速的56%，深度大于8 km處破裂速度設(shè)定為剪切波速的80%，5～8 km的區(qū)域設(shè)置為線性過渡區(qū)域，因此破裂速度沿深度范圍的初始分布可以表示為

(16)

進(jìn)而根據(jù)子源中心至破裂起始點的直線距離和初始破裂速度分布即可確定斷層面上各子源的破裂起始時刻Tij0，而后考慮各子源滑動量對各自破裂起始時刻的時間擾動并引入隨機數(shù)，得到破裂起始時刻的空間分布：

(17)

式中，M0表示總地震矩大小，sij表示各子源的滑動量，sA表示斷層面上的平均滑動量，sM表示所有子源中最大滑動量，ε表示服從標(biāo)準(zhǔn)正態(tài)分布的隨機數(shù)，σT是對數(shù)正態(tài)標(biāo)準(zhǔn)差，此處設(shè)置為0.2.

考慮淺層地殼(z<5 km)的破裂持時是較深處的(z>8 km)的2倍，其間為線性過度，各子源的上升時間可以表示為

(18)

最后根據(jù)斷層面上的平均上升時間τA(公式(19))按比例調(diào)整各子源的上升時間，即可得其在斷層面上的空間分布.

(19)

式中，c1為計算常量，取為1.45×10-9，αT為縮放系數(shù)，與斷層的滑動角(λ)和傾角(δ)相關(guān):

αT=[1+FDFRcα]-1,

(20)

(21)

(22)

FD和FR分別為傾角和滑動角因子，cα為針對逆斷層的高頻輻射因子.

此外，同樣為了考慮各子源滑動量對各自上升時間的時間擾動通過引入隨機數(shù)，得到各個子源上升時間的空間分布：

τi=τ0iexp(εσR)，

(23)

式中，τ0ij表示各個子源的初始上升時間，ε表示服從標(biāo)準(zhǔn)正態(tài)分布的隨機數(shù)，σR是對數(shù)正態(tài)標(biāo)準(zhǔn)差，此處設(shè)置為0.5.

對于一個設(shè)定地震，基于SPECFEM 3D建立GP混合運動學(xué)震源模型的整個流程(圖1)可以歸納為如下7個步驟：

圖1 基于SPECFEM 3D混合運動學(xué)震源建模流程圖Fig.1 Flow chart of hybrid kinematic source based on SPECFEM 3D

(1)借助全局參數(shù)定標(biāo)律，估計破裂面的尺寸和平均滑動量.進(jìn)一步將整個破裂面劃分為NL×NM個矩形小網(wǎng)格.

(2)借助全局參數(shù)與局部參數(shù)間的定標(biāo)律，估計最大凹凸體的尺寸和平均錯動量，并將其定位在破裂面上，然后再根據(jù)整個破裂面上的平均錯動量估計凹凸體外斷層的錯動量.

(3)借助二維傅里葉變換將上述確定性的錯動分布轉(zhuǎn)換得到確定性波數(shù)譜，借助式(12)引入的隨機波數(shù)譜，將確定性的低波數(shù)譜和隨機高波數(shù)譜結(jié)合，最后二維逆傅里葉變換至空間域中，得到各個網(wǎng)格的錯動量.

(4)根據(jù)(3)計算得到的各個子源的錯動量，計算得到每個子源的矩張量，并確定各個子源的位置，生成譜元程序SPECFEM 3D中的要震源參數(shù)文件CMTSOLUTION.

(5)根據(jù)斷層埋深確定破裂速度和上升時間的初始分布，而后考慮各子源滑動量對各子源的破裂速度和上升時間的擾動并引入隨機數(shù)因子，得到破裂速度和上升時間的時間分布.

(6)根據(jù)子源對應(yīng)的破裂時間和上升時間，選定的震源函數(shù)形式以及時間步長的大小，得到每個子源對應(yīng)的震源時間.

(7)由于每個子源時間函數(shù)不同，故需要從主參數(shù)文件Par_file中設(shè)置EXTERNAL_SOURCE_FILE為true，并且在震源參數(shù)文件CMTSOLUTION中每個子源中額外加入對應(yīng)的時間函數(shù)文件的路徑.

2 精度驗證：一維波速結(jié)構(gòu)地震動模擬

為驗證SEM模擬寬頻帶地震動的精度，本節(jié)與基于修正動力剛度矩陣的FK方法(Ba et al., 2021, 2022)進(jìn)行對照.前文提及的FK方法作為一種半解析方法，模擬精度高、模擬帶寬范圍較廣，適合作為對比驗證的參照結(jié)果.故本節(jié)分別基于SEM和FK方法在一維波速空間下針對同一混合震源模型進(jìn)行模擬，比較兩種方法模擬得到0～10 Hz的地震動時程和頻譜曲線.

本文基于上述兩種方法分別建立了均勻和多層的一維波速結(jié)構(gòu)模型，兩種模型大小均為60 km×60 km×40 km，速度結(jié)構(gòu)參數(shù)分別見圖2(a,b).本節(jié)與FK方法對比的模擬結(jié)果的頻帶范圍為0～10 Hz，SEM中為保證計算結(jié)果的精確可靠，要求最短波長中至少包含5個GLL節(jié)點，需要根據(jù)模擬頻率范圍確定網(wǎng)格大小.模型采用Trelis軟件進(jìn)行劃分，圖2(a，b)分別為均勻和多層模型譜單元網(wǎng)格圖，其中單層模型網(wǎng)格總數(shù)約為270萬，GLL節(jié)點數(shù)量為1.94億個；多層模型網(wǎng)格總數(shù)約為340萬，GLL節(jié)點數(shù)量為2.54億個.

圖2 均勻(a)和多層(b)一維波速結(jié)構(gòu)模型圖中VP為P波波速，VS為S波波速，ρ為密度，Q為品質(zhì)因子.Fig.2 Homogeneous (a) and multi-layer (b) 1D wave velocity structure modelsVP is P-wave velocity, VS is S-wave velocity, ρ is density, Q is quality factor.

根據(jù)上述兩種一維波速結(jié)構(gòu)模型分別建立對應(yīng)的混合震源模型.本文選擇建立MW6.0走滑斷層的混合震源模型，首先根據(jù)定標(biāo)律確定斷層的全局震源參數(shù)，見表1所示.

表1 斷層全局震源參數(shù)Table 1 Fault global source parameters

考慮到本次地震震級小于6.5級，斷層面上采用單凹凸體，根據(jù)姜偉等(2017)中走滑斷層的定標(biāo)率公式，計算得到震源的局部參數(shù)見表2所示.由于斷層在均勻半空間和多層一維波速結(jié)構(gòu)情況下的埋深不同，斷層位置周圍對應(yīng)的波速不一致，根據(jù)1.2節(jié)中式(11)—(15)計算斷層面滑動量分布，式(16)—(23)計算得到斷層各個子源的破裂時間和上升時間，最終得到的均勻和多層一維波速結(jié)構(gòu)模型的震源模型如圖3所示.

圖3 均勻(a)和多層(b)一維波速結(jié)構(gòu)模型對應(yīng)的混合運動學(xué)震源每個子圖右上角的數(shù)字三值組分別表示給定分布的最小值、平均值和最大值.Fig.3 Homogeneous (a) and muti-layer (b) hybrid kinematic source model Triplet of numbers at top right of each panel indicates the minimum, mean and maximum values of the given distribution, respectively.

表2 斷層局部震源參數(shù)Table 2 Fault local source parameters

限于篇幅，僅給出了由表面處沿著斷層方向和垂直斷層方向震中距為10 km位置處(見圖4)加速度三分量時程和對應(yīng)的頻譜結(jié)果.其中，圖5為均勻一維波速結(jié)構(gòu)模型結(jié)果，圖6為多層一維波速結(jié)構(gòu)模型結(jié)果.結(jié)果表明：基于SEM模擬的加速度時程波形、持時和幅值以及頻譜結(jié)果在0～10 Hz的頻率范圍內(nèi)與FK方法的結(jié)果一致，驗證了基于GP混合運動學(xué)震源模型在SEM模擬的精度.

圖4 觀測點P1和P2位置Fig.4 Observation points P1 and P2

圖5 P1(a)和P2(b)均勻一維波速結(jié)構(gòu)模型計算結(jié)果Fig.5 Calculation results of P1 (a) and P2 (b) homogeneous 1D wave velocity structure models

圖6 P1(a)和P2(b)多層一維波速結(jié)構(gòu)模型計算結(jié)果Fig.6 Calculation results of P1 (a) and P2 (b) multilayer 1D wave velocity structure models

3 適用性檢驗：2021年漾濞6.4級地震動模擬

為檢驗SEM模擬寬頻地震動的適用性，本節(jié)選擇模擬2021年漾濞6.4級(0～5 Hz)地震，將模擬結(jié)果與臺站觀測記錄的時程和反應(yīng)譜的比較，及與NGA-West2地震動衰減方程的反應(yīng)譜曲線的比較，檢驗了方法的適用性.

2021年5月21日云南大理州漾濞縣6.4級地震造成數(shù)十萬人受災(zāi)，數(shù)萬間房屋倒塌，直接經(jīng)濟損失高達(dá)33.2億元.中國應(yīng)急管理部將此次地震列為2021年全國十大自然災(zāi)害.根據(jù)中國地震局工程力學(xué)研究所提供的強震動觀測數(shù)據(jù)，漾濞MS6.4地震共觸發(fā)28個自由場強震動觀測臺站，除053YPX臺站只獲取了兩個水平分量的記錄外，其余臺站均獲得了完整的三分量強震動加速度記錄.基于上述豐富地震記錄，許多學(xué)者和研究人員開展了近斷層地震動模擬特征的相關(guān)研究，何欣娟和潘華(2021)采用隨機有限斷層法模擬了目標(biāo)區(qū)域峰值加速度的分布；周紅等(2021)采用NNSIM隨機有限斷層法模擬了近斷層200 km范圍內(nèi)的峰值分布特征；強生銀等(2021)基于隨機有限斷層三維地震動模擬方法，給出了漾濞地震中28個強震動臺站的三分量加速度時程模擬結(jié)果.上述模擬方法大多數(shù)選用隨機方法進(jìn)行模擬，由于隨機方法僅能合成一維水平分量的高頻地震動并且無明確方位.本節(jié)采用GP14.3震源模型和SEM對該地震進(jìn)行0～5 Hz地震動模擬，模擬計算區(qū)域見圖7所示，其中模型經(jīng)度范圍99.4°E—100.4°E、緯度25.4°N—26.0°N，區(qū)域范圍內(nèi)包含053YBX、053DLY、053YPX和053BTH臺站.根據(jù)目標(biāo)區(qū)域得到計算模型沿東西長約100 km，南北寬約60 km，縱向約30 km.本文以“中國大陸淺層結(jié)構(gòu)模型”和“中國大陸巖石圈速度結(jié)構(gòu)模型”為主，參考CRUST1.0，并依據(jù)漾濞地區(qū)的地面高程數(shù)據(jù)，構(gòu)建了包含起伏地形的三維速度結(jié)構(gòu)模型，包含介質(zhì)的壓縮波速VP、介質(zhì)的剪切波速VS，質(zhì)量密度ρ、品質(zhì)因子Q以及層厚度等，具體參數(shù)見表3所示.

圖7 漾濞地震模擬區(qū)域Fig.7 Yangbi seismic simulation area

表3 漾濞地區(qū)速度結(jié)構(gòu)參數(shù)Table 3 Velocity structure parameters of Yangbi area

本文設(shè)定能模擬的最大頻率為5 Hz，網(wǎng)格總數(shù)約為1169×104，GLL節(jié)點數(shù)量高達(dá)7.65×108個，時間步距取為0.005 s，模擬90 s內(nèi)的地震波傳播，采用國家超算“天河一號”400核計算耗時7.2 h.

基于1.2節(jié)的GP14.3混合震源模型建立漾濞地震的混合運動學(xué)震源模型參數(shù)包括全局震源參數(shù)和局部震源參數(shù).本文的全局震源參數(shù)采用張斌等(2021)研究給出的包括斷層面的幾何信息、產(chǎn)狀(走向角、傾角及上邊緣埋深)及破裂方式(破裂形式和滑動角)等，具體見表4.將斷層面劃分為15×15個1.0 km×1.0 km的矩形子源.震源局部參數(shù)(主要是凹凸體模型參數(shù))基于局部參數(shù)定標(biāo)率確定，因本次地震震級小于6.5級，斷層面上僅采用單凹凸體，具體局部參數(shù)及定標(biāo)率見表5.根據(jù)1.2節(jié)公式計算斷層面各個子源的時空分布，得到漾濞地震混合運動學(xué)震源模型見圖8所示.

表4 漾濞地震斷層全局震源參數(shù)Table 4 Global source parameters of the Yangbi earthquake fault

表5 漾濞地震斷層局部震源參數(shù)Table 5 Local source parameters of the Yangbi earthquake fault

圖8 漾濞地震混合運動學(xué)震源模型每個子圖右上角的數(shù)字三值組分別表示給定分布的最小值、平均值和最大值.Fig.8 Hybrid kinematic source model of the Yangbi earthquake Triplet of numbers at top right of each panel indicates the minimum, mean and maximum values of the given distribution, respectively.

3.1 實際觀測記錄對比

圖9給出了模擬區(qū)域內(nèi)共4個臺站的模擬結(jié)果與強震記錄的加速度時程和反應(yīng)譜的結(jié)果對比.其中反應(yīng)譜結(jié)果左上角標(biāo)注了每個臺站的名稱，黑線代表NS方向，紅線代表EW方向，藍(lán)線代表UD方向.反應(yīng)譜結(jié)果下方為模擬結(jié)果與強震記錄對應(yīng)的加速度三分量時程結(jié)果，方向與反應(yīng)譜結(jié)果的方向保持一致，上方為觀測結(jié)果，下方為模擬結(jié)果.每條時程的持時均為50 s，PGA均標(biāo)注于曲線的右上方，單位為cm·s-2.值得注意的是，圖9中給出的強震記錄已經(jīng)帶通濾波至0.1～5 Hz.對比結(jié)果表明：本文模擬的加速度時程的波形、持時和幅值與強震記錄對比良好，特別是與震中距7.9 km的53YBX臺站對比結(jié)果較為一致，并且模擬結(jié)果在周期0.2～5 s的PGA反應(yīng)譜與強震記錄吻合的PGA反應(yīng)譜結(jié)果較為一致，總體上表明了本文提出的基于SEM進(jìn)行寬頻地震動模擬以及應(yīng)用GP14.3的混合震源模型的針對實際地震模擬的適用性.由于本文建立的含起伏地表的地殼層波速結(jié)構(gòu)模型，暫時未能考慮地表低波速土體對地震波的影響，故隨著震中距的增加，速度結(jié)構(gòu)對地震波傳播的影響更加明顯，導(dǎo)致震中距較遠(yuǎn)的臺站記錄與模擬結(jié)果差別相對明顯.

圖9 模擬與(YBXDLYYPXBTH)臺站記錄的加速度時程和反應(yīng)譜(Sa)對比Fig.9 The simulation is compared with acceleration time histories and reaction spectra (Sa) recorded at (YBXDLYYPXBTH) stations

3.2 GMPE衰減方程對比

為了進(jìn)一步對比各個周期下，特別是基于物理的寬頻地震動模擬得到高頻結(jié)果的適用性，本文將模擬結(jié)果與NGA-West2提供的地震動衰減關(guān)系曲線進(jìn)行對比，對比結(jié)果見圖10所示.其中，采用Dangkua等 (2018)基于地震動數(shù)據(jù)分析后，計算得到的適用于中國大陸地區(qū)的地震動衰減關(guān)系，即將NGA-West2中ASK14、BSSA14、CB14、CY14和IM14每個衰減關(guān)系均取0.2的權(quán)重計算得到地震動衰減關(guān)系.同時，為與上述衰減關(guān)系進(jìn)行對比，我們將模擬得到的兩個水平正交分量(NS和EW)分量的加速度時程在水平面進(jìn)行旋轉(zhuǎn)，計算0°～90°范圍內(nèi)每旋轉(zhuǎn)θ角度后的加速度時程，將不同周期下各水平分量的反應(yīng)譜值作幾何平均并求其中位數(shù)得到GMRotD50.將模擬結(jié)果得到不同周期下GMRotD50結(jié)果隨距破裂頂部邊緣的水平距離Rx(km)的變化與周期0.2～5 s的GMPE衰減方程進(jìn)行對比，圖中紅線代表反應(yīng)譜結(jié)果的平均值，上下兩條黑線分別代表±σ的誤差范圍.對比結(jié)果顯示，在0.2～5 s的周期內(nèi)高、低頻的模擬結(jié)果與衰減關(guān)系曲線均取得了良好的一致性，進(jìn)一步說明了上述模擬方法的可靠性.

圖10 不同周期下模擬的加速度反應(yīng)譜RotD50結(jié)果和GMPE曲線對比結(jié)果PSa是偽加速度反應(yīng)譜，σ是標(biāo)準(zhǔn)差，Rx是沿垂直于斷層走向方向、觀測點到斷層頂面的距離.Fig.10 The results of simulated acceleration response spectrum RotD50 are compared with the GMPE curves under different periods PSa is pseudo spectral acceleration, σ is standard deviation, Rx is horizontal distance from top of rupture measured perpendicular to fault strike.

根據(jù)上述地震動時程、峰值、反應(yīng)譜以及衰減關(guān)系的充分比較分析，驗證了本文提出方法的可靠性，總體上說明本文采用的SEM和混合震源模型適用于2021年漾濞6.4級地震的寬頻地震動模擬，同時，通過模擬給出了2021年漾濞地震下區(qū)域的峰值地震動空間分布特征，見3.3節(jié)所示.

3.3 地面運動峰值PGA和PGV分布

圖11(a，b)分別是目標(biāo)區(qū)域的水平方向峰值加速度(PGA)和峰值速度(PGV)分布圖，水平分量由NS分量和EW分量通過矢量合成方式獲得.觀測圖中峰值較大區(qū)域集中在斷層投影面附近，這與混合震源模型中凹凸體的位置相對應(yīng)，體現(xiàn)了近斷層地震動的集中性效應(yīng)；垂直于斷層走向的PGA和PGV比平行于斷層方向的衰減更快，這是由于斷層呈中心破裂，釋放的能量持續(xù)在破裂前端匯聚，體現(xiàn)了近斷層地震動的破裂方向性效應(yīng).

圖11 地震動PGA(a)和PGV(b)分布圖Fig.11 Distribution of ground motion PGA (a) and PGV (b)

研究區(qū)域PGA分布形態(tài)，震中點臨近漾濞縣城，最大值接近400 cm·s-2，在近斷層區(qū)域，整個斷層位置區(qū)域，PGA均高于300 cm·s-2，隨著震中距的增加，地震動出現(xiàn)快速衰減，區(qū)域內(nèi)從300～50 cm·s-2衰減迅速，在永平縣附近PGA低于50 cm·s-2.觀察區(qū)域內(nèi)PGV分布發(fā)現(xiàn)，極震區(qū)PGV接近45 cm·s-1，漾濞縣城范圍內(nèi)PGV達(dá)到40 cm·s-1，區(qū)域內(nèi)PGV比PGA的衰減相對緩慢，永平縣附近PGV達(dá)到20 cm·s-1.同時值得關(guān)注的是，由于地表起伏的影響，PGA和PGV的峰值分布出現(xiàn)明顯的不均勻性，特別是大理和洱海西側(cè)位置受局部凸起地形的影響導(dǎo)致地震動能量聚集，放大了地震動響應(yīng)，上述地區(qū)出現(xiàn)的高烈度異?，F(xiàn)象值得注意.

4 結(jié)論

本文將GP14.3有限斷層混合運動學(xué)震源開發(fā)到SPECFEM 3D開源代碼中，實現(xiàn)了基于確定性物理模型的全過程寬頻帶地震動模擬.文中首先根據(jù)定標(biāo)律建立了凹凸體震源模型，進(jìn)而通過在斷層面的滑動分布中結(jié)合Graves & Pitarka提出的GP14.3隨機震源模型得到混合運動學(xué)震源模型.最后將上述混合運動學(xué)震源模型計算得到每個子源的矩張量以及對應(yīng)的時間函數(shù)按照對應(yīng)格式輸入到SPECFEM 3D譜元程序中.

為檢驗該方法的精度和可靠性，文中分別建立一維(均勻和多層)和漾濞地區(qū)的波速結(jié)構(gòu)模型以及對應(yīng)的混合運動學(xué)震源模型，將得到的SEM模擬結(jié)果分別與FK方法的計算結(jié)果和2021年6.4級漾濞地震記錄進(jìn)行對比，并且論文最后還給出了漾濞地震下近場強地面運動的空間分布特征.主要結(jié)論如下：

(1)針對均勻和多層一維波速結(jié)構(gòu)模型，SEM模擬結(jié)果與作者建立的FK方法的計算結(jié)果在寬頻帶(0～10 Hz)吻合良好，驗證了提出方法的精度;針對云南漾濞6.4級地震，模擬結(jié)果與區(qū)域內(nèi)4個臺站觀測記錄的時程和反應(yīng)譜的結(jié)果接近，并且模擬得到地震動衰減規(guī)律與NGA-West2地震動預(yù)測方程在頻率0.1～5 Hz的反應(yīng)譜結(jié)果具有較好的一致性，檢驗了方法的適用性.

(2)2021年漾濞地震震中PGA接近400 cm·s-2，PGV達(dá)到45 cm·s-1，極震區(qū)烈度最高可達(dá)到Ⅸ度，并且受局部地形起伏影響，大理以及洱海西側(cè)位置出現(xiàn)了高烈度異常區(qū).

致謝感謝中國地震局工程力學(xué)研究所為本研究提供了2021年5月21日漾濞6.4級地震近場區(qū)域的地震觀測記錄．感謝國家超級計算天津中心“天河一號”為本文計算提供的幫助．感謝三位匿名審稿專家對本文提出的修改意見.感謝Dimitri Komatitsch和Jeroen Tromp等人提供的Specfem3D_Cartesian譜元法計算程序，程序的下載網(wǎng)址為(https:∥github.com/geodynamics/specfem3d.git).