模擬地震波傳播的三維逐元并行譜元法

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

1 應急管理部國家自然災害防治研究院， 北京 100085 2 清華大學數學科學系， 北京 1000843 中國地質大學地球物理與空間信息學院， 武漢 430074 4 中國科學院地質與地球物理研究所， 巖石圈演化國家重點實驗室， 北京 100029

0 引言

基于地震波運動方程的地震成像方法越來越受到研究者的關注(Tape et al.,2009; Tan and Huang,2014; Liu et al.,2018; Lei et al.,2020). 在此類地震成像原理中,無論是伴隨成像方法(Tromp et al.,2005),還是逆時偏移方法(Dai et al.,2012),由于都需要反復數值求解地震波運動方程,因此數值算法的計算效率直接關系到地震成像效率. 經過四十余年的發展,研究者已開發出多種數值方法用于地震波運動方程的直接求解,如：有限差分法(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)等.

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

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

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

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

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

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

需要注意的是,本文作者前期對逐元并行SEM的研究中計算核心之間通信采用的是主從(Master-Slave)通信模型(Liu et al.,2017a),雖然主從通信模式程序上容易實現(Komatitsch and Tromp,1999),但計算核心通信等待時間較長.為了減少通信等待時間,本文采用緩存通信模式,在計算核心對之間通信,大幅減少SEM的通信時間.此外,本文不再采用存儲子矩陣的方式用于重建單元剛度矩陣,而是存儲單元雅克比矩陣的逆及其行列式的值,再結合插值多項式直接對空間微分算子數值近似,因此避免了重建稠密單元剛度矩陣,減少了計算量.為了消除邊界截斷而引起的虛假反射,本文在模型邊界上利用逐元并行SEM求解二階位移形式的PML吸收邊界條件.

1 逐元譜元法

地震波在震源處激發,向三維空間輻射能量,地震波的傳播滿足以下運動方程:

(1a)

(1b)

(2)

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

(3)

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

(4)

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

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

(5)

其中Ym為m階Legendre多項式,其上標表示一階導數.公式(5)可確定m+1個根ξi(i=0,…,m),即在每個方向上插值節點數為m+1,三維情況下一共有(m+1)3個插值節點.在每個單元上利用Lagrange多項式近似函數值,對u和v的近似可表示為:

(6)

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

(7)

(8)

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

Me′=Γe(M).

(9)

在同一個單元上Me′和Me不相等,原因在于相鄰單元共享GLL節點,與共享節點對應的質量矩陣元素需累加,使得Me′中共享節點對應的質量矩陣元素值要大于Me中的元素值.

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

(10)

(11)

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

(12)

(12)式可進一步寫成矩陣的形式:

(13)

(14)

(15)

(16)

(17)

(18)

假設地震震源為點源,且可以寫成關于矩張量的形式,公式(1)中的體力f可表示為:

(19)

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

(20)

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

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

(21)

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

(22)

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

(23)

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

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

(24)

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

(25)

(26)

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

2 逐元并行策略

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

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

需要指出的是,為了討論方便,圖1只顯示相鄰計算核心之間通信.事實上,任何計算核心可能與其他多個計算核心構成通信對進行通信.至于如何構成通信對,取決于模型復雜程度以及模型剖分以后生成的網格.

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

3 PML吸收邊界條件

PML吸收邊界條件已被廣泛應用于地震波數值模擬之中(羅玉欽和劉財,2020),為了本文的完整性本節簡要介紹本文所采用的二階位移形式PML吸收邊界條件.有限的計算資源難以模擬高頻(>2 Hz)地震波在整個地球中的傳播(Tong et al.,2014).實際地震波場模擬時往往將模型截斷,然后在小規模模型中模擬高頻地震波傳播.在使用SEM模擬地震波傳播時,直接將模型截斷不做任何處理,相當于在截斷邊界處引入自由地表邊界條件,此時必然造成強烈的虛假反射而干擾模型內部地震波場的可靠性.為了防止邊界截斷而引入的虛假反射,本文在邊界截斷處引入PML吸收邊界條件(Liu Y S et al.,2014; Liu et al.,2017a).

構造PML吸收邊界條件的思想是在頻率域PML吸收層中做坐標變換而引入衰減因子使得波場沿著垂直截斷邊界方向指數衰減.PML吸收層中的坐標變換公式如下(Komatitsch and Tromp,2003):

(27)

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

(28)

4 數值算例

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

圖2 均勻介質模型用于測試逐元并行SEM的有效性坐標原點位于模型上表面中心,五角心表示震源,坐標為(0 km,0 km,0.5 km),三角形表示臺站,兩個臺站為R1和R2,分別位于(0 km,20 km,0 km)和(20 km,20 km,0 km)處.介質參數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 理論測試

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

圖3 逐元并行SEM合成波形與理論波形對比(a)、(c)和(e)分別為臺站R1處位移的x、y和z分量; (b)、(d)和(f)分別為臺站R2處位移的x、y和z分量.黑線由逐元并行SEM計算得到; 紅色線由解析解計算得到.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)方法得到)得到的波形幾乎完全重合,兩者之間并未出現明顯差別.無論是振幅較小的直達P波還是振幅較大的Rayleigh面波都得到了精確模擬,說明逐元并行SEM用于地震波模擬的有效性.兩個臺站R1和R2靠近模型邊界,但從合成地震圖上并未發現來自模型邊界的虛假反射波,說明地震波傳播至PML吸收層時已被有效吸收.

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

圖4 逐元并行SEM并行效率曲線圖中黑點為實際計算時間,黑色曲線為理論計算時間.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.

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

4.2 實際地震波場模擬

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

圖5 2013年四川蘆山MW6.6地震震中及周邊地形圖五角心表示震中,其經緯度坐標為(102.888°E,30.308°N) (USGS提供,http:∥www.usgs.gov).小方塊表示震中附近國家測震臺網的4個固定臺.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年蘆山地震震源區網格化模型地表地形高程從全球地形網格中提取得到.模型內部界面從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) 矩張量釋放率隨時間變化; (b) 震源時間函數Fig.7 (a) Moment rate function; (b) Source time function

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

圖8 合成波形與實際觀測波形對比圖(a)—(d)中波形分別由臺站1—4記錄.紅色線為觀測波形; 黑線為合成波形.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年蘆山地震地表波場快照圖(a)—(i)分別為t=1 s、10 s、20 s、30 s、40 s、50 s、60 s、70 s和80 s時的地表垂直位移.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 結論和討論

本文發展了一種高效的逐元并行SEM用于地震波場模擬.并在邊界處構造二階位移形式的PML吸收邊界條件避免邊界截斷而產生的虛假反射.通過公式推導、數值實驗、分析與對比,證實本文給出的逐元并行SEM主要表現出三方面的優點: (1)并行性強; (2)計算量小; (3)內存消耗少.此外,基于逐元并行SEM原理,我們開發出一套三維地震波模擬軟件包,將其稱之為EBE-SEM3D.

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

本文的EBE-SEM3D軟件包可無償用于學術用途,為了獲得軟件包請與本文第一作者郵件聯系.

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

附錄A

由公式(4)知,通過物理坐標對局部坐標求偏導可計算雅克比矩陣行列式的值,即:

(A1)

由(A1)式可獲得單元雅克比矩陣在每個GLL點處行列式的值.對雅克比矩陣求逆,可得到局部坐標對物理坐標偏導數的值,即:

(A2)