劉少林, 楊頂輝, 徐錫偉, 李小凡, 申文豪, 劉有山

1 應(yīng)急管理部國(guó)家自然災(zāi)害防治研究院， 北京 100085 2 清華大學(xué)數(shù)學(xué)科學(xué)系， 北京 1000843 中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院， 武漢 430074 4 中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 巖石圈演化國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 100029

0 引言

基于地震波運(yùn)動(dòng)方程的地震成像方法越來(lái)越受到研究者的關(guān)注(Tape et al.,2009; Tan and Huang,2014; Liu et al.,2018; Lei et al.,2020). 在此類地震成像原理中,無(wú)論是伴隨成像方法(Tromp et al.,2005),還是逆時(shí)偏移方法(Dai et al.,2012),由于都需要反復(fù)數(shù)值求解地震波運(yùn)動(dòng)方程,因此數(shù)值算法的計(jì)算效率直接關(guān)系到地震成像效率. 經(jīng)過(guò)四十余年的發(fā)展,研究者已開(kāi)發(fā)出多種數(shù)值方法用于地震波運(yùn)動(dòng)方程的直接求解,如：有限差分法(FDM)(Madariaga,1976; Graves,1996; Zhang and Yao,2013),偽譜法(PSM)(Kosloff and Baysal,1982; Furumura et al.,1998; Wang and Takenaka,2011),有限元法(FEM)(Marfurt,1984; Padovani et al.,1994; Liu et al.,2014a),譜元法(SEM)(Seriani and Priolo,1994; Komatitsch and Vilotte,1998; Liu et al.,2017a)等.

在以上提到的數(shù)值算法中,FDM使用最為廣泛(Yang et al.,2013).通過(guò)對(duì)地震波方程直接離散,FDM的優(yōu)點(diǎn)在于算法較為直觀、程序易于實(shí)現(xiàn)和并行(Graves,1996).但FDM的不足也較為明顯,其缺點(diǎn)主要表現(xiàn)在三個(gè)方面.首先,傳統(tǒng)FDM在使用低階算子時(shí)數(shù)值頻散明顯(Zhang and Yao,2013),其較強(qiáng)的數(shù)值頻散必然降低算法精度.雖然高階插值能減弱FDM的數(shù)值頻散,但是高階插值不僅會(huì)增加計(jì)算量,同時(shí)差分算子半徑過(guò)長(zhǎng)造成并行計(jì)算時(shí)通信量的增加,降低了FDM并行計(jì)算效率(宋國(guó)杰等,2012).其次,FDM網(wǎng)格不靈活,難以運(yùn)用于復(fù)雜模型中的地震波場(chǎng)模擬(Shragge,2014).雖然可通過(guò)坐標(biāo)變換的方式將偽正交網(wǎng)格變換成正交網(wǎng)格(Rao and Wang,2013; Zhang et al.,2014),然后在正交網(wǎng)格上采用常規(guī)差分格式求解地震波運(yùn)動(dòng)方程,實(shí)現(xiàn)復(fù)雜不規(guī)則模型中的地震波場(chǎng)模擬.但當(dāng)模型復(fù)雜程度高時(shí)會(huì)造成網(wǎng)格嚴(yán)重畸變,使得曲線網(wǎng)格FDM穩(wěn)定性變差,以至于影響地震波模擬效率.第三,FDM較難處理自由地表邊界條件(Ruud and Hestholm,2001).雖然通過(guò)設(shè)置虛擬層的方式可方便滿足自由地表邊界條件,但該處理方式對(duì)自由地表邊界條件近似程度較低(蘇波等,2019).針對(duì)FDM存在的缺陷或不足,研究者做出了許多有益研究(Yang et al.,2003; Zhang and Yao,2013; 劉少林等,2013；Liu,2014),例如：Yang等(2003)引入質(zhì)點(diǎn)位移的梯度信息,聯(lián)合質(zhì)點(diǎn)位移用于空間微分算子的數(shù)值近似,該方法有效解決了FDM數(shù)值頻散大的弊端,但FDM網(wǎng)格的不靈活性以及自由邊界不易處理等問(wèn)題依然存在.

FEM可采用非結(jié)構(gòu)化網(wǎng)格離散復(fù)雜模型,具有很強(qiáng)的網(wǎng)格靈活性(Cohen,2002),同時(shí),將FEM形成的地震波運(yùn)動(dòng)方程弱形式所包含的邊界積分項(xiàng)直接令為零,即可方便實(shí)現(xiàn)自由地表邊界條件.雖然FEM具有一定的優(yōu)越性,但是FEM計(jì)算量和存儲(chǔ)量較大(Bielak et al.,2005),在大尺度地震波場(chǎng)模擬時(shí)對(duì)計(jì)算機(jī)的計(jì)算性能要求較高.此外,FEM在單元上插值逼近時(shí),往往只能采用低階插值,因此精度較低,在采用高階插值時(shí)出現(xiàn)Runge現(xiàn)象,造成精度的損失(Fichtner,2011).盡管研究者在提升FEM計(jì)算效率方面做出了許多改進(jìn),例如：通過(guò)構(gòu)造集中質(zhì)量矩陣的方式減少FEM的計(jì)算量(Padovani et al.,1994),采用存儲(chǔ)核矩陣的方式減少FEM的存儲(chǔ)量(Liu S L et al., 2014; Meng and Fu,2017),但FEM的計(jì)算與存儲(chǔ)量都明顯高于FDM.

在簡(jiǎn)單模型中,傅里葉PSM具有較高的計(jì)算精度,PSM算法的核心在于通過(guò)快速Fourier變換對(duì)空間微分算子逼近(Kosloff and Baysal,1982),其理論精度可達(dá)到FDM無(wú)窮階時(shí)的精度(龍桂華等,2009).在相同空間網(wǎng)格采樣的條件下,PSM的數(shù)值頻散小于FDM.但PSM由于使用全局長(zhǎng)算子,一方面較難適用于復(fù)雜分層、固液耦合模型(Wang et al.,2001); 另一方面長(zhǎng)算子用于空間微分近似,有悖于微分算子局部性的特點(diǎn)(Li et al.,2011).此外,與FDM類似,PSM難以具備網(wǎng)格的靈活性也較難處理自由地表邊界條件.

相對(duì)而言,SEM是較為高效的地震波模擬方法,該方法不僅具有譜方法的高精度性也具有FEM網(wǎng)格的靈活性.在20世紀(jì)80年代,SEM最早在流體力學(xué)領(lǐng)域提出并得到應(yīng)用(Patera,1984).隨后,Seriani和Priolo(1994)將SEM引入到地震學(xué)領(lǐng)域,并模擬了二維非均勻介質(zhì)中聲波傳播.Seriani和Priolo(1994)的數(shù)值模擬結(jié)果表明,在相同空間網(wǎng)格采樣條件下,譜元法的精度要高于FEM,例如：在空間插值階數(shù)為8階、地震波的最短波長(zhǎng)采樣點(diǎn)數(shù)為4.5時(shí)SEM即可以獲得高精度模擬結(jié)果(Seriani and Priolo,1994).早期的SEM在正交多項(xiàng)式的選取上往往采用Chebyshev正交多項(xiàng)式(Seriani and Priolo,1994; Seriani,1997,1998),雖然基于Chebyshev正交多項(xiàng)式的SEM取得了較好的數(shù)值模擬結(jié)果,但Chebyshev多項(xiàng)式帶權(quán)正交(Wang et al.,2007),以至于該類SEM無(wú)法形成對(duì)角質(zhì)量矩陣,在求解線性方程組時(shí)占用較大計(jì)算資源(Seriani,1998).雖然該類SEM可采用類似于FEM質(zhì)量集中的方式形成對(duì)角質(zhì)量矩陣,但這必然造成SEM精度損失(Seriani and Oliveira,2007).

為了保證SEM的計(jì)算精度,同時(shí)避免求解大型線性方程組,基于Legendre正交多項(xiàng)式的SEM被發(fā)展起來(lái),并被運(yùn)用于地震波場(chǎng)的數(shù)值模擬(Komatitsch and Vilotte,1998).該類SEM的優(yōu)勢(shì)在于積分點(diǎn)與插值點(diǎn)重合,使得質(zhì)量矩陣為對(duì)角矩陣,從而避免求解大型線性方程組.基于Legendre正交多項(xiàng)式的譜元法已被廣泛應(yīng)用于區(qū)域尺度和全球尺度的地震波模擬與成像之中(Komatitsch et al.,2002; Dong et al.,2019; Lei et al.,2020).

剛度矩陣與解向量的乘積為譜元法計(jì)算量的主要來(lái)源.在并行計(jì)算時(shí),由于模型的復(fù)雜性,不規(guī)則網(wǎng)格剖分模型會(huì)導(dǎo)致分配到每個(gè)計(jì)算核心的單元數(shù)不同,使得每個(gè)計(jì)算核心所分配的剛度矩陣數(shù)量不同,最終導(dǎo)致計(jì)算核心負(fù)載不均衡,影響了并行計(jì)算效率.為了提升SEM的并行計(jì)算效率,Liu等(2017a)提出逐元并行譜元法,其思想是將譜元單元平均分配到每個(gè)計(jì)算核心,在單元上進(jìn)行剛度矩陣與解向量相乘.數(shù)值實(shí)驗(yàn)表明,逐元并行策略能有效提升SEM的并行計(jì)算效率.逐元并行SEM除了提升并行效率以外,還將剛度矩陣的存儲(chǔ)轉(zhuǎn)化為64個(gè)子矩陣的存儲(chǔ),大幅降低了SEM的存儲(chǔ)量.

需要注意的是,本文作者前期對(duì)逐元并行SEM的研究中計(jì)算核心之間通信采用的是主從(Master-Slave)通信模型(Liu et al.,2017a),雖然主從通信模式程序上容易實(shí)現(xiàn)(Komatitsch and Tromp,1999),但計(jì)算核心通信等待時(shí)間較長(zhǎng).為了減少通信等待時(shí)間,本文采用緩存通信模式,在計(jì)算核心對(duì)之間通信,大幅減少SEM的通信時(shí)間.此外,本文不再采用存儲(chǔ)子矩陣的方式用于重建單元?jiǎng)偠染仃?而是存儲(chǔ)單元雅克比矩陣的逆及其行列式的值,再結(jié)合插值多項(xiàng)式直接對(duì)空間微分算子數(shù)值近似,因此避免了重建稠密單元?jiǎng)偠染仃?減少了計(jì)算量.為了消除邊界截?cái)喽鸬奶摷俜瓷?本文在模型邊界上利用逐元并行SEM求解二階位移形式的PML吸收邊界條件.

1 逐元譜元法

地震波在震源處激發(fā),向三維空間輻射能量,地震波的傳播滿足以下運(yùn)動(dòng)方程:

(1a)

(1b)

(2)

其中n為邊界外法向,Ω為研究區(qū)域.考慮自由地表處法向應(yīng)力為0,即T·n=0,認(rèn)為地震波在無(wú)限半空間中傳播,即無(wú)窮遠(yuǎn)處邊界上應(yīng)力為0,即T=0,公式(2)可簡(jiǎn)化為:

(3)

公式(3)自動(dòng)滿足自由地表邊界條件.將研究區(qū)域Ω離散成N個(gè)非重疊六面體單元,假設(shè)每個(gè)單元的控制點(diǎn)數(shù)為n,物理空間中的每個(gè)六面體單元都可變換成標(biāo)準(zhǔn)參考單元,同樣每個(gè)參考單元可變換回至六面體單元,物理域與參考域之間的坐標(biāo)滿足以下關(guān)系(Komatitsch and Tromp,1999):

(4)

其中xa為控制點(diǎn)在物理域中的坐標(biāo),Na(ξ,η,λ)(-1≤ξ≤1,-1≤η≤1,-1≤λ≤1)為單元形函數(shù).在參考單元上求解如下方程確定每個(gè)方向上插值節(jié)點(diǎn)的坐標(biāo):

(1-ξ2)Y′m(ξ)=0，

(5)

其中Ym為m階Legendre多項(xiàng)式,其上標(biāo)表示一階導(dǎo)數(shù).公式(5)可確定m+1個(gè)根ξi(i=0,…,m),即在每個(gè)方向上插值節(jié)點(diǎn)數(shù)為m+1,三維情況下一共有(m+1)3個(gè)插值節(jié)點(diǎn).在每個(gè)單元上利用Lagrange多項(xiàng)式近似函數(shù)值,對(duì)u和v的近似可表示為:

(6)

其中L為L(zhǎng)agrange插值多項(xiàng)式,α,β和γ為單元e上三個(gè)方向GLL(Gauss-Lobatto-Legendre)節(jié)點(diǎn)索引.由公式(6)可知,單元上的插值多項(xiàng)式由三個(gè)方向上的Lagrange插值多項(xiàng)式組合而成.利用公式(6),公式(2)等號(hào)左邊積分可以表示為:

(7)

(8)

其中∪為矩陣組裝算子.定義Γe算子將總體質(zhì)量矩陣投影至單元上:

Me′=Γe(M).

(9)

在同一個(gè)單元上Me′和Me不相等,原因在于相鄰單元共享GLL節(jié)點(diǎn),與共享節(jié)點(diǎn)對(duì)應(yīng)的質(zhì)量矩陣元素需累加,使得Me′中共享節(jié)點(diǎn)對(duì)應(yīng)的質(zhì)量矩陣元素值要大于Me中的元素值.

與公式(7)類似,公式(3)中等號(hào)右邊第一式也可寫(xiě)成離散形式:

(10)

(11)

將(6)和(11)式代入(10)式并利用GLL數(shù)值求積,(10)式可離散為:

(12)

(12)式可進(jìn)一步寫(xiě)成矩陣的形式:

(13)

(14)

(15)

(16)

(17)

(18)

假設(shè)地震震源為點(diǎn)源,且可以寫(xiě)成關(guān)于矩張量的形式,公式(1)中的體力f可表示為:

(19)

其中M為地震矩張量(Aki and Richards,1980),S(t)為震源時(shí)間函數(shù),xs為震源坐標(biāo).將(19)式代入(3)式,公式(3)等式右邊第二項(xiàng)可表示為:

(20)

+Lα(ξs)L′β(ηs)Lλ(λs)Gi2

+Lα(ξs)Lβ(ηs)L′λ(λs)Gi3]，

(21)

其中(ξs,ηs,λs)為震源坐標(biāo)xs所對(duì)應(yīng)的局部坐標(biāo),下標(biāo)es表示震源所在的單元.(21)式可寫(xiě)成向量相乘的形式:

(22)

其中單元上的載荷向量Fe具體形式為:

(23)

由(21)—(23)式不難看出,震源項(xiàng)離散以后所形成的載荷向量Fe只在一個(gè)單元上有值,在其他單元上均為零.如果震源恰好位于單元表面而被多個(gè)單元共享,此時(shí)Fe可在多個(gè)單元上有值.為了簡(jiǎn)化程序設(shè)計(jì),可在震源位置施加以小的擾動(dòng)(如單元尺寸的1/100，本文算例中選擇1/300)，將震源放置到單元內(nèi)部,保證Fe只在一個(gè)單元上有值.數(shù)值算例顯示該處理方式不會(huì)造成明顯誤差,關(guān)于精度的討論本文將在第4節(jié)數(shù)值算例部分具體給出.需要注意的是,本文公式(21)—(23)與Komatitsch和Tromp(1999)的表述不同,Komatitsch和Tromp(1999)的處理方式是將點(diǎn)源平滑至單元的每個(gè)GLL節(jié)點(diǎn)上,本文公式(21)—(23)仍將震源當(dāng)成點(diǎn)源.當(dāng)單元尺寸遠(yuǎn)小于地震波長(zhǎng)時(shí),Komatitsch和Tromp(1999)的震源處理方式不會(huì)產(chǎn)生較大誤差; 當(dāng)單元尺寸和地震波長(zhǎng)相當(dāng)時(shí),該處理方式可能產(chǎn)生明顯誤差.Komatitsch和Tromp(1999)震源處理存在的問(wèn)題已超出本文范圍,本文不再討論.

由(7)、(13)和(22)式可知,以下等式成立:

(24)

由于測(cè)試向量v的任意性,(24)式等價(jià)為:

(25)

(26)

從公式(26)容易看出,本文先在每個(gè)單元上做矩陣Ke和向量Pe相乘運(yùn)算,再將單元上的加速度向量組裝獲得總體加速度向量.因?yàn)楣?26)采用逐元(element-by-element)的方式計(jì)算地震波場(chǎng),所以本文稱之為逐元SEM.該計(jì)算方式與前人在合成總體剛度矩陣以后再做矩陣與向量相乘不同(Padovani et al.,1994; Liu Y S et al., 2014),(26)式所體現(xiàn)的優(yōu)勢(shì)主要表現(xiàn)在三方面.首先,極大簡(jiǎn)化了并行程序設(shè)計(jì).SEM的計(jì)算量主要來(lái)源于單元上的矩陣與向量乘積,將單元上的矩陣與向量乘積分配到計(jì)算核心就可以方便實(shí)現(xiàn)SEM的并行計(jì)算.其次,節(jié)省了存儲(chǔ)空間.由(7)、(15)—(17)式不難看出,逐元SEM不需要存儲(chǔ)稠密的單元?jiǎng)偠染仃?只用存儲(chǔ)單元GLL點(diǎn)上的雅克比矩陣(Jacobian matrix)行列式的值以及雅克比矩陣的逆(計(jì)算單元雅克比矩陣逆的過(guò)程見(jiàn)附錄A),也就是需要存儲(chǔ)元素的個(gè)數(shù)為GLL點(diǎn)數(shù)的10倍,因此逐元SEM的存儲(chǔ)量和FDM的存儲(chǔ)量在同一數(shù)量級(jí)(為FDM的5倍左右).最后,節(jié)省了計(jì)算量.由(12)式可知逐元SEM充分利用到GLL積分與插值重合的特點(diǎn),一個(gè)方向上的空間微分?jǐn)?shù)值近似,只保留對(duì)應(yīng)方向的非零元,因此極大減少了浮點(diǎn)運(yùn)算次數(shù).如果空間GLL點(diǎn)的個(gè)數(shù)為H,插值階數(shù)為m,那么逐元SEM的浮點(diǎn)運(yùn)算次數(shù)大約為O(18(m+1)H),其計(jì)算量與FDM的計(jì)算量大致相同.

2 逐元并行策略

在關(guān)于SEM的并行程序設(shè)計(jì)中,研究者往往采用主從(master-slave)通信模式(Komatitsch and Tromp,2002; Liu et al.,2017a)完成每個(gè)時(shí)間迭代步的數(shù)據(jù)通信.SEM的主從通信過(guò)程分為三步,首先,一個(gè)計(jì)算核心作為主計(jì)算核心收集從計(jì)算核心需要通信的信息(如加速度向量); 然后,主計(jì)算核心對(duì)通信信息處理(如不同計(jì)算核心共有的GLL節(jié)點(diǎn)對(duì)應(yīng)的加速度相加); 最后,主計(jì)算核心將處理好的信息分發(fā)給從計(jì)算核心.顯然這三步中主要通信和計(jì)算任務(wù)都由主計(jì)算核心完成,而從計(jì)算核心承擔(dān)的通信和計(jì)算任務(wù)遠(yuǎn)小于主計(jì)算核心,較多時(shí)間浪費(fèi)在通信等待上,以至于降低了SEM的并行計(jì)算效率.

為了提升SEM的并行效率,本文采用緩存通信模式(Gropp et al.,1996),在需要通信的計(jì)算核心對(duì)之間兩兩通信,具體通信過(guò)程如圖1所示,為了計(jì)算流程更加完整,圖1也顯示了通信以外的其他過(guò)程.由圖1可知,在合成單元質(zhì)量矩陣以后,需要在計(jì)算核心之間通信完成之后獲得Γe(M).在每一個(gè)時(shí)間步,每個(gè)單元上需要計(jì)算KePe,然后將計(jì)算核心之間共享GLL點(diǎn)上的加速度在通信之后累加得到完整的質(zhì)點(diǎn)加速度.從圖1不難看出,與主從通信模式不同,本文的通信方式將總的通信量分配給計(jì)算核心對(duì),有效減少了通信等待時(shí)間.此外,在主從通信模式中主計(jì)算核心的數(shù)據(jù)處理任務(wù)(共享GLL點(diǎn)上加速度累加求和)被分配給各個(gè)計(jì)算核心,有效平衡了計(jì)算核心的數(shù)據(jù)計(jì)算加載.因此,本文通信模式有效提升了SEM的并行計(jì)算效率.有關(guān)并行計(jì)算效率,本文將在第4節(jié)進(jìn)一步討論.

需要指出的是,為了討論方便,圖1只顯示相鄰計(jì)算核心之間通信.事實(shí)上,任何計(jì)算核心可能與其他多個(gè)計(jì)算核心構(gòu)成通信對(duì)進(jìn)行通信.至于如何構(gòu)成通信對(duì),取決于模型復(fù)雜程度以及模型剖分以后生成的網(wǎng)格.

圖1 逐元并行SEM地震波場(chǎng)模擬流程圖Fig.1 Schematic showing the processes of seismic wave modeling using element-by-element parallel spectral-element method

3 PML吸收邊界條件

PML吸收邊界條件已被廣泛應(yīng)用于地震波數(shù)值模擬之中(羅玉欽和劉財(cái),2020),為了本文的完整性本節(jié)簡(jiǎn)要介紹本文所采用的二階位移形式PML吸收邊界條件.有限的計(jì)算資源難以模擬高頻(>2 Hz)地震波在整個(gè)地球中的傳播(Tong et al.,2014).實(shí)際地震波場(chǎng)模擬時(shí)往往將模型截?cái)?然后在小規(guī)模模型中模擬高頻地震波傳播.在使用SEM模擬地震波傳播時(shí),直接將模型截?cái)嗖蛔鋈魏翁幚?相當(dāng)于在截?cái)噙吔缣幰胱杂傻乇磉吔鐥l件,此時(shí)必然造成強(qiáng)烈的虛假反射而干擾模型內(nèi)部地震波場(chǎng)的可靠性.為了防止邊界截?cái)喽氲奶摷俜瓷?本文在邊界截?cái)嗵幰隤ML吸收邊界條件(Liu Y S et al.,2014; Liu et al.,2017a).

構(gòu)造PML吸收邊界條件的思想是在頻率域PML吸收層中做坐標(biāo)變換而引入衰減因子使得波場(chǎng)沿著垂直截?cái)噙吔绶较蛑笖?shù)衰減.PML吸收層中的坐標(biāo)變換公式如下(Komatitsch and Tromp,2003):

(27)

將(27)式代入(1)式,忽略震源項(xiàng),得到:

(28)

4 數(shù)值算例

本節(jié)通過(guò)兩個(gè)數(shù)值算例討論本文發(fā)展的逐元并行SEM在三維地震波場(chǎng)模擬時(shí)的有效性和高效性.數(shù)值算例采用二階中心差分用于時(shí)間離散(Liu et al.,2017b); PML吸收層厚度為5個(gè)譜單元.

圖2 均勻介質(zhì)模型用于測(cè)試逐元并行SEM的有效性坐標(biāo)原點(diǎn)位于模型上表面中心,五角心表示震源,坐標(biāo)為(0 km,0 km,0.5 km),三角形表示臺(tái)站,兩個(gè)臺(tái)站為R1和R2,分別位于(0 km,20 km,0 km)和(20 km,20 km,0 km)處.介質(zhì)參數(shù)P波速度為5.8 km·s-1,S波速度為3.46 km·s-1,密度為2.72 g·cm-3.Fig.2 3D homogeneous model for the benchmark test of EBE-SEMThe origin of the Cartesian coordinates is located at the center of the top surface of the model. The star denotes an explosive source, and is placed at (0 km,0 km,0.5 km). The triangles represent two stations R1 and R2,which are at (0 km,20 km,0 km) and (20 km,20 km,0 km),respectively. The media parameters are described by VP,VS,and ρ,which are 5.8 km·s-1,3.46 km·s-1,2.72 g·cm-3,respectively.

4.1 理論測(cè)試

為了測(cè)試逐元并行SEM的有效性,選擇如圖2所示的均勻各向同性介質(zhì)模型.模型水平尺度為70 km,深度為30 km.模型介質(zhì)參數(shù)由P波速度(5.8 km·s-1)、S波速度(3.46 km·s-1)和密度(2.72 g·cm-3)描述.笛卡爾坐標(biāo)原點(diǎn)位于模型上表面中心.如圖2所示x軸正向指向正北,y軸正向指向正東,z軸正向垂直向下.震源為爆炸源位于(0 km,0 km,0.5 km)處,其對(duì)應(yīng)的矩張量只在對(duì)角線上有值,對(duì)角線上的值都為2×1016N·m.兩個(gè)臺(tái)站R1和R2分別位于(0 km,20 km,0 km)和(20 km,20 km,0 km)處.震源時(shí)間函數(shù)選擇Ricker子波,其主頻為1 Hz,時(shí)間函數(shù)最高頻率約為2.5 Hz.采用正方體離散圖2中規(guī)則模型,正方體的邊長(zhǎng)為1 km.在每個(gè)單元上采用4階多項(xiàng)式插值.波形模擬結(jié)果如圖3所示.

圖3 逐元并行SEM合成波形與理論波形對(duì)比(a)、(c)和(e)分別為臺(tái)站R1處位移的x、y和z分量; (b)、(d)和(f)分別為臺(tái)站R2處位移的x、y和z分量.黑線由逐元并行SEM計(jì)算得到; 紅色線由解析解計(jì)算得到.Fig.3 Synthetic waveforms at stations R1 (a，c，e) and R2 (b，d，f) computed by EBE-SEM (black lines) and analytical method (red lines)(a) and (b) are for x component of displacement; (c) and (d) are for y component of displacement; (e) and (f) are for z component of displacement.

從圖3可以看出,逐元并行SEM得到的合成波形與解析解(由Cagniard-de Hoop(De Hoop,1959)方法得到)得到的波形幾乎完全重合,兩者之間并未出現(xiàn)明顯差別.無(wú)論是振幅較小的直達(dá)P波還是振幅較大的Rayleigh面波都得到了精確模擬,說(shuō)明逐元并行SEM用于地震波模擬的有效性.兩個(gè)臺(tái)站R1和R2靠近模型邊界,但從合成地震圖上并未發(fā)現(xiàn)來(lái)自模型邊界的虛假反射波,說(shuō)明地震波傳播至PML吸收層時(shí)已被有效吸收.

同樣采用圖2中的模型用于測(cè)試逐元并行SEM的并行計(jì)算效率.逐元并行SEM最少需要2個(gè)計(jì)算核心進(jìn)行并行計(jì)算,在使用2個(gè)計(jì)算核心計(jì)算時(shí)總耗時(shí)為t2,由于計(jì)算核心之間通信時(shí)間遠(yuǎn)小于用于矩陣與向量乘積的計(jì)算時(shí)間,此時(shí)可忽略通信耗時(shí),當(dāng)使用i個(gè)計(jì)算核心計(jì)算時(shí),理論計(jì)算時(shí)間為ti=2t2/i.由于計(jì)算核心增多時(shí),通信量隨之增加,因此實(shí)際計(jì)算時(shí)間要大于理論計(jì)算時(shí)間ti.如果并行計(jì)算時(shí)間越靠近ti說(shuō)明并行性越好.逐元并行SEM的并行計(jì)算效率如圖4所示.從圖4可看出,實(shí)際并行計(jì)算時(shí)間接近理論并行時(shí)間.當(dāng)計(jì)算核心數(shù)為素?cái)?shù)(如11、19)時(shí)每個(gè)計(jì)算核心所分配的空間網(wǎng)格較為復(fù)雜,使得通信量相對(duì)增加,以至于實(shí)際并行計(jì)算時(shí)間出現(xiàn)異常(圖4).即便出現(xiàn)并行計(jì)算時(shí)間異常(如計(jì)算核心數(shù)為11和19),也未出現(xiàn)在主從通信模式中計(jì)算核心增多而計(jì)算時(shí)間增加的情況(Liu et al.,2017a),說(shuō)明本文采用的計(jì)算核心對(duì)通信方式能有效平衡分配給計(jì)算核心的計(jì)算和通信任務(wù).事實(shí)上,將本文的通信模式換成主從通信模式,在使用20個(gè)計(jì)算核心并行計(jì)算時(shí)耗時(shí)增加約2倍,使用3至20個(gè)計(jì)算核心平均計(jì)算時(shí)間增加1倍左右.

圖4 逐元并行SEM并行效率曲線圖中黑點(diǎn)為實(shí)際計(jì)算時(shí)間,黑色曲線為理論計(jì)算時(shí)間.Fig.4 The parallel efficiency of the EBE-SEMThe black line with solid circles shows the real computing time,while the black line is the theoretical computing time.

應(yīng)該指出的是,逐元并行SEM除了并行效率提升以外,在節(jié)省內(nèi)存方面具有巨大優(yōu)勢(shì).與本文作者前期工作利用64個(gè)單元子矩陣重建單元?jiǎng)偠染仃?Liu et al.,2017a)不同,本文的逐元SEM不再需要重建單元?jiǎng)偠染仃?(12)式,(15)—(17)式),而只用存儲(chǔ)GLL點(diǎn)上雅克比矩陣的逆及行列式的值(一共10個(gè)元素),因此比Liu等(2017a)的方法節(jié)省約6倍的存儲(chǔ).對(duì)于圖2中的三維模型,單元個(gè)數(shù)為196875,每個(gè)單元上的GLL點(diǎn)數(shù)為125,存儲(chǔ)雅克比矩陣的逆及行列式的值一共消耗0.23 GB內(nèi)存,本次模擬一共消耗約2 GB內(nèi)存(額外需要存儲(chǔ)當(dāng)前時(shí)刻位移三分量、上一時(shí)刻位移三分量、加速度三分量、空間GLL點(diǎn)坐標(biāo)、空間GLL點(diǎn)編號(hào)等).

4.2 實(shí)際地震波場(chǎng)模擬

為了驗(yàn)證逐元并行SEM在實(shí)際地震模擬時(shí)的有效性,本節(jié)模擬2013年蘆山MW6.6地震波場(chǎng)傳播.蘆山地震震中區(qū)及周邊地形如圖5所示.從圖中可看出,模型存在較大的地形起伏.模型西部為青藏高原隆起區(qū),海拔較高; 東部主要為四川盆地,海拔較低.從全球地形網(wǎng)格文件(Tozer et al.,2019)中提取研究區(qū)域內(nèi)的高程,然后將高程文件導(dǎo)入到Hypermesh剖分軟件(僅限學(xué)術(shù)用途使用)中.除了地表地形以外,模型的其他界面(Conrad、Moho面等)由CRUST1.0模型(Laske et al.,2013)提供,模型深度到120 km.由Hypermesh剖分軟件生成的網(wǎng)格如圖6所示,離散模型采用六面體單元,其平均尺寸為4 km.本文采用張小艷等(2020)給出的震源機(jī)制解得到震源矩張量以及矩張量釋放率函數(shù)(圖7a).通過(guò)對(duì)矩張量釋放率函數(shù)做時(shí)間積分得到震源時(shí)間函數(shù)S(t)(圖7b).實(shí)際地震波場(chǎng)模擬結(jié)果如圖8和圖9所示.

圖5 2013年四川蘆山MW6.6地震震中及周邊地形圖五角心表示震中,其經(jīng)緯度坐標(biāo)為(102.888°E,30.308°N) (USGS提供,http:∥www.usgs.gov).小方塊表示震中附近國(guó)家測(cè)震臺(tái)網(wǎng)的4個(gè)固定臺(tái).Fig.5 Topography of the area around the epicenter of 2013 Lushan MW6.6 earthquakeThe star denotes the epicenter of the earthquake,the coordinates of which are (102.888°E,30.308°N) (http:∥www.usgs.gov). The squares are broadband seismic stations from China National Seismic Network.

圖6 2013年蘆山地震震源區(qū)網(wǎng)格化模型地表地形高程從全球地形網(wǎng)格中提取得到.模型內(nèi)部界面從CRUST1.0中提取得到.模型深度為120 km.Fig.6 Mesh for the model of the area around the 2013 Lushan EarthquakeTopography is derived from the global topography grid. The subsurface is extracted from the CRUST1.0. The depth of the model is 120 km.

圖7 (a) 矩張量釋放率隨時(shí)間變化; (b) 震源時(shí)間函數(shù)Fig.7 (a) Moment rate function; (b) Source time function

圖8中紅色線為實(shí)際地震儀記錄的波形經(jīng)過(guò)去儀器響應(yīng)、低通濾波(最高截止頻率為0.2 Hz)以后的波形,圖中黑色線為逐元并行SEM的合成波形.總體而言,數(shù)值方法能較好模擬直達(dá)P波和Rayleigh面波初至到時(shí),但相位有較大差別.需要注意的是,圖8a中數(shù)值解顯示Rayleigh面波主要能量到時(shí)要早于觀測(cè)到時(shí),主要原因在于模型速度不準(zhǔn).模型淺部(圖6,深度小于20 km)S波速度為3.55 km·s-1,而四川盆地淺部(深度小于20 km)S波平均速度約為2.5 km·s-1(Laske et al.,2013),如果利用S波速度近似面波速度,模擬的面波到時(shí)約比觀測(cè)早20 s,這與圖8a中的結(jié)果基本一致.另外,為了模擬得到更準(zhǔn)確的波形,在震源近似上需要采用有限斷層近似(公式(23)是對(duì)點(diǎn)源的離散而不是有限斷層的離散)以及采用更加精確的地震波速度結(jié)構(gòu)模型.圖9展示的是9個(gè)時(shí)刻的地面波場(chǎng)(垂直分量)快照.當(dāng)t=1 s時(shí)(圖9a),地震波主要能量還未傳播至地表,地表位移幾乎為0.當(dāng)t=10 s時(shí)(圖9b),震中處表現(xiàn)出較強(qiáng)的位移.除了震中處較強(qiáng)的位移以外,圖9也顯示地表傳播的直達(dá)P波和Rayleigh面波成分,其中Rayleigh面波振幅遠(yuǎn)大于P波.圖9f顯示地震波傳播到PML吸收層時(shí)被很好地吸收而未見(jiàn)來(lái)自模型邊界的反射波.

圖8 合成波形與實(shí)際觀測(cè)波形對(duì)比圖(a)—(d)中波形分別由臺(tái)站1—4記錄.紅色線為觀測(cè)波形; 黑線為合成波形.Fig.8 Waveform comparisonThe waveforms in plots (a)—(d) are recorded by stations 1—4,respectively. The red lines are the real observed seismograms,while the black lines are the synthetic seismograms.

圖9 2013年蘆山地震地表波場(chǎng)快照?qǐng)D(a)—(i)分別為t=1 s、10 s、20 s、30 s、40 s、50 s、60 s、70 s和80 s時(shí)的地表垂直位移.Fig.9 Wavefield snapshots of the 2013 Lushan Earthquake. The seismic wavefields of vertical component are on the free surfaceThe snapshots in plots (a)—(i) are taken at t=1 s,10 s,20 s,30 s,40 s,50 s,60 s,70 s,and 80 s,respectively.

5 結(jié)論和討論

本文發(fā)展了一種高效的逐元并行SEM用于地震波場(chǎng)模擬.并在邊界處構(gòu)造二階位移形式的PML吸收邊界條件避免邊界截?cái)喽a(chǎn)生的虛假反射.通過(guò)公式推導(dǎo)、數(shù)值實(shí)驗(yàn)、分析與對(duì)比,證實(shí)本文給出的逐元并行SEM主要表現(xiàn)出三方面的優(yōu)點(diǎn): (1)并行性強(qiáng); (2)計(jì)算量小; (3)內(nèi)存消耗少.此外,基于逐元并行SEM原理,我們開(kāi)發(fā)出一套三維地震波模擬軟件包,將其稱之為EBE-SEM3D.

通過(guò)逐元并行SEM求解地震波場(chǎng)和伴隨波場(chǎng)(Tromp et al.,2005),可在EBE-SEM3D的基礎(chǔ)上發(fā)展一套伴隨層析成像軟件包,通過(guò)構(gòu)造遠(yuǎn)震入射邊界條件(Liu et al,2017a)也可方便開(kāi)發(fā)出一套遠(yuǎn)震伴隨成像軟件包.關(guān)于區(qū)域地震和遠(yuǎn)震的伴隨成像軟件包本文作者將另文介紹和給出.

本文的EBE-SEM3D軟件包可無(wú)償用于學(xué)術(shù)用途,為了獲得軟件包請(qǐng)與本文第一作者郵件聯(lián)系.

致謝感謝兩位審稿專家提出的寶貴修改意見(jiàn).感謝Julien Diaz教授關(guān)于解析解的交流、討論以及為本研究提供解析解源代碼.感謝中國(guó)地震局地球物理研究所國(guó)家測(cè)震臺(tái)網(wǎng)數(shù)據(jù)備份中心(doi:10.11998/SeisDmc/SN)為本研究提供地震波形數(shù)據(jù).本文圖由開(kāi)源軟件GMT(Generic Mapping Tools) (Wessel et al.,2013)繪制而成.逐元并行SEM軟件包中并行通訊庫(kù)采用的是開(kāi)源庫(kù)MPICH (http:∥www.mpich.org).

附錄A

由公式(4)知,通過(guò)物理坐標(biāo)對(duì)局部坐標(biāo)求偏導(dǎo)可計(jì)算雅克比矩陣行列式的值,即:

(A1)

由(A1)式可獲得單元雅克比矩陣在每個(gè)GLL點(diǎn)處行列式的值.對(duì)雅克比矩陣求逆,可得到局部坐標(biāo)對(duì)物理坐標(biāo)偏導(dǎo)數(shù)的值,即:

(A2)