李冰非, 李小凡, 李峰, 龔飛
1 中國地震災(zāi)害防御中心, 北京 100029 2 中國科學(xué)院地質(zhì)與地球物理研究所, 北京 100029 3 中國地質(zhì)大學(xué)地球物理與空間信息學(xué)院, 武漢 430074
地球自由振蕩是由大地震在有限尺度的地球內(nèi)部激發(fā)的三維駐波,即整體地球的固有振蕩(Benioff et al., 1961).地球自由振蕩的固有頻率和頻譜形態(tài)與地球內(nèi)部結(jié)構(gòu)密切相關(guān)(Lapwood and Usami, 1981),包括地球內(nèi)部橫向非均勻性及徑向非均勻性(Tsuboi et al., 1985; Tsuboi,1995)、地球內(nèi)部各向異性(Woodhouse et al., 1986; Tromp, 1993,1995)、介質(zhì)非完全彈性(Masters and Widmer, 1995)、以及介質(zhì)對(duì)波的衰減(Davies and Bunge, 2001; Alterman et al., 1974; Pollitz et al., 1987; Park et al., 2005)等性質(zhì).由于應(yīng)用駐波方法研究地球內(nèi)部結(jié)構(gòu)具有對(duì)采樣點(diǎn)要求不高的優(yōu)勢,地球自由振蕩數(shù)值模擬研究可以方便、經(jīng)濟(jì)地重建全球多尺度三維地球內(nèi)部結(jié)構(gòu)及獲得地球內(nèi)部某些物理場的三維特征參數(shù)分布.
譜元法是近年來大尺度地震波及地球自由振蕩數(shù)值模擬的重要方法之一(Komatitsch et al., 2002, 2004; Komatitsch and Tromp, 2002; Tsuboi et al., 2004; Yan et al., 2014).譜元法將有限元算法和譜方法相結(jié)合,兼顧了有限元方法的靈活性和譜方法的高精度(Komatitsch and Vilotte, 1998; Liu et al., 2017a).譜元法是彈性波方程在空間近似的高階變分方法,Komatitsch和Tromp(1999)在譜元單元上采用Gauss-Lobatto-Legendre(GLL)積分并結(jié)合高階精度的拉格朗日插值基函數(shù),得到具有對(duì)角化的質(zhì)量矩陣的顯式時(shí)間域微分方程組,且數(shù)值解更好地趨近于理論解,模擬精度和計(jì)算效率得以提高;而后又將譜元方法應(yīng)用于全球地震波傳播模擬研究(2002).Tromp等(2010)在多年研究的基礎(chǔ)上,將數(shù)值方法和高性能并行計(jì)算相結(jié)合,研制譜元法程序包SPECFEM3D_GLOBE,可以用來精細(xì)剖分地球模型、模擬三維地球模型或區(qū)域模型中地震波傳播.趙明等(2013)拓展了一種譜元-簡正振型耦合方法適用于各向異性介質(zhì)中地震波傳播數(shù)值模擬.近二十年來,研究者們將超大規(guī)模高性能并行計(jì)算數(shù)值模擬方法應(yīng)用于全球地震波傳播研究,將計(jì)算結(jié)果與實(shí)際觀測的超長波記錄對(duì)比研究,獲得了一系列較為令人滿意的結(jié)果(劉少林等, 2021).特別應(yīng)該提到的是,嚴(yán)珍珍等(2008a)應(yīng)用譜元法結(jié)合彈性波理論及高性能并行計(jì)算對(duì)由特大地震激發(fā)的地球自由振蕩進(jìn)行了成功的大規(guī)模數(shù)值模擬,從而證實(shí)了這類方法的可行性,為由大地震激發(fā)地球自由振蕩的數(shù)值模擬研究奠定了基礎(chǔ).
盡管理論上時(shí)間域數(shù)值方法均可處理諸如地球內(nèi)部橫向非均勻性及地震激發(fā)地球自由振蕩之全過程等問題,但現(xiàn)有地球自由振蕩時(shí)間域數(shù)值模擬技術(shù)均基于非保結(jié)構(gòu)計(jì)算,其長時(shí)程計(jì)算中的高積累誤差難以避免,很難精確處理地球自由振蕩模擬這一長時(shí)程追蹤問題.近年來已有聲波及彈性波保結(jié)構(gòu)模擬方法問世(羅明秋等, 2001; Li et al., 2012; 汪文帥等, 2012; 陳世仲等, 2016; Liu et al., 2017b,c; 李冰非等, 2019);Li 等(2012)及Gao等(2016)進(jìn)一步將逆時(shí)離散變換與三階辛相結(jié)合,獲得了長時(shí)程的低頻散數(shù)值模擬,這類保結(jié)構(gòu)方法非常適合處理有關(guān)保守動(dòng)力學(xué)體系的問題.若將地球近似為彈性球體,用上述保結(jié)構(gòu)方法模擬其自由振蕩問題將是很好的選擇;然而,地球本身是一具有衰減性質(zhì)的非完全彈性球體,真實(shí)地球介質(zhì)極其復(fù)雜, 除了各向異性以外,可能還有多相多孔性以及黏性(Kuznetsov, 1997; Kustowski et al., 2008a),其自由振蕩(駐波)問題涉及非保守動(dòng)力學(xué)系統(tǒng)問題.針對(duì)這一問題,Li等(2015)發(fā)展了阻尼聲波保結(jié)構(gòu)辛算法相關(guān)理論,數(shù)理邏輯的合理性和方法的可行性均得到驗(yàn)證.Cai等(2017)基于同樣的模型方程提出了另一種耗散保結(jié)構(gòu)模擬方法,對(duì)于帶耗散的地震波運(yùn)動(dòng)波動(dòng)方程,可以分為一個(gè)守恒系統(tǒng)和一個(gè)耗散系統(tǒng)來求解.為更精確地模擬衰減模式下的地球自由振蕩,需要發(fā)展新的、高效駐波保結(jié)構(gòu)數(shù)值方法,研究適用于衰減介質(zhì)的保彌散數(shù)值模擬方法.
本文將保結(jié)構(gòu)計(jì)算延拓至非保守動(dòng)力學(xué)體系,發(fā)展了一種更切合實(shí)際的、適用于處理地球自由振蕩衰減問題的高精度、高效率模擬計(jì)算方法,以2012年印度尼西亞蘇門答臘西北海域8.6級(jí)大地震觀測記錄為數(shù)據(jù)源,進(jìn)行了衰減模式下的全球自由振蕩數(shù)值實(shí)驗(yàn),研究地球深部(如地核)彌散衰減對(duì)地球自由振蕩的影響和各種衰減模式對(duì)地球自由振蕩的影響.
地球的一般彈性運(yùn)動(dòng)方程為
(1)
其中ρ為密度,u為質(zhì)點(diǎn)位移,ω為地球角轉(zhuǎn)動(dòng)矢量,f表示外力和重力的擾動(dòng),λ,μ為Lame常數(shù)(嚴(yán)珍珍等,2008b).
為便于表達(dá),我們對(duì)地球的一般彈性運(yùn)動(dòng)方程進(jìn)行簡化,忽略大氣變化因素,忽略自轉(zhuǎn)、重力、地球橢率以及橫向不均勻性,考慮簡化的地球彈性運(yùn)動(dòng)方程為:
(2)
其中F表示震源.數(shù)值模擬地球自由振蕩,本質(zhì)上也是求滿足適當(dāng)初邊值條件的地球振動(dòng)的常系數(shù)偏微分方程組的數(shù)值解.
用譜元法將方程組(2)離散化后可寫為:
(3)
其中M為全局質(zhì)量矩陣,K為全局剛度矩陣,U表示全局位移向量.
對(duì)于本文研究的問題,彌散衰減介質(zhì)中彈性波方程可寫為:
(4)
其中u為位移矢量,σ為應(yīng)力張量,ε為應(yīng)變張量,c為四階彈性張量,D為彌散阻尼系數(shù),f為外力.譜元法離散空間域后的微分方程組可寫為:
(5)
將式(5)應(yīng)用阻尼項(xiàng)保結(jié)構(gòu)辛格式(Li et al., 2015)進(jìn)行時(shí)間離散,可得如下基于譜元法的保彌散衰減波動(dòng)方程的時(shí)間格式:
(6)
其中V為速度向量. Δt為時(shí)間步長,c2,d1,d2為辛系數(shù),U(n)和V(n)的上標(biāo)表示第n時(shí)間步,U1和V1為中間變量.
本文將采用基于(6)式的辛格式-譜元法開展衰減模式下的地球自由振蕩數(shù)值模擬研究,并探究地核中不同衰減模式及其對(duì)地球自由振蕩的影響.
本節(jié)通過數(shù)值模擬實(shí)驗(yàn)驗(yàn)證本文提出的方法的有效性.采用2012年4月11日印度尼西亞蘇門答臘西北海域8.6級(jí)地震(2.327°N,93.063°E)作為震源,圖1所示為本次地震事件震源位置和Global Centroid-Moment-Tensor (CMT) Project提供的8.6級(jí)地震的震源球以及全球地震臺(tái)網(wǎng)(Global Seismographic Network, GSN)公布的152個(gè)臺(tái)站所在位置(https:∥www.iris.edu/hq/programs/gsn).數(shù)值模擬使用Computational Infrastructure for Geodynamics(CIG)提供的SPECFEM3D-GLOBE譜元法程序包(經(jīng)本文保彌散-辛格式升級(jí)).
圖1 地震事件震源和全球152個(gè)臺(tái)站位置示意圖(https:∥www.iris.edu/hq/programs/gsn)Fig.1 Epicenter, 152 stations on globe and earthquake mechanism
本節(jié)分別采用保彌散-辛格式和SPECFEM3D-GLOBE譜元法程序包中所采用的Newmark算法模擬同一模型,考慮地形、橢率、重力及自轉(zhuǎn)等因素,模擬地震波傳播600 min.采用S362ANI地球模型(Kustowski et al., 2008a,b),并將數(shù)值模擬結(jié)果進(jìn)行同參數(shù)處理.圖2所示為NNA臺(tái)站Z分量地震記錄波形對(duì)比,從圖中波形放大圖可以看出,模擬時(shí)長約388 min時(shí)Newmark算法的結(jié)果開始顯現(xiàn)出高頻數(shù)值頻散(波形曲線出現(xiàn)高頻毛刺),本文算法的波形記錄則比Newmark算法在長時(shí)間計(jì)算后更加穩(wěn)定,顯示出在模擬長時(shí)程地震波傳播時(shí)的優(yōu)越性.采用同樣的時(shí)間步長(dt=0.1425 s),在同樣的計(jì)算節(jié)點(diǎn)條件下,本文保彌散-辛格式模擬耗時(shí)68 h 23 m 24 s,Newmark算法模擬耗時(shí)71 h 05 m 10 s,本文保彌散-辛格式模擬長時(shí)程地震波傳播時(shí)更加穩(wěn)定,同時(shí)不增加計(jì)算量甚至計(jì)算效率略有提高.本實(shí)驗(yàn)僅簡要展示低階辛格式處理大尺度地震波傳播問題的有效性及抑制數(shù)值頻散的能力,同類辛格式的數(shù)值頻散前人已有詳盡的分析(Li et al., 2015).
本小節(jié)設(shè)計(jì)兩組實(shí)驗(yàn),其整體均采用S362ANI地球模型,并考慮地形、橢率、重力及自轉(zhuǎn)等因素.第一組實(shí)驗(yàn)采用本文保彌散衰減方法模擬外核和內(nèi)核彌散衰減模型(設(shè)內(nèi)核和外核彌散阻尼D=-0.1 s-1),和采用SPECFEM3D-GLOBE譜元法程序包自帶的Newmark算法模擬無彌散衰減的模型,模擬2012年印度尼西亞北蘇門答臘8.6級(jí)地震發(fā)生后地震波傳播600 min (dt=0.1425 s,模擬252700步).將數(shù)值模擬結(jié)果與實(shí)際地震觀測數(shù)據(jù)進(jìn)行對(duì)比分析,本文采用的實(shí)際地震觀測數(shù)據(jù)來源為IRIS(Incorporated Research Institutions for Seismology).根據(jù)地球自由振蕩的甚低頻特性,我們將兩種數(shù)值模擬算法的結(jié)果和實(shí)際地震觀測結(jié)果(IRIS下載)采用0.0039 Hz的低通濾波,選擇Chirp-z變換算法求取頻譜進(jìn)行對(duì)比,并引入地球自由振蕩理論值(簡振模)作為參考.由于外核中不存在環(huán)型振蕩,本文僅研究地核彌散阻尼對(duì)球型振蕩的影響.且隨著簡正模球諧階l的增大,自由振蕩所發(fā)生的深度越來越淺(Jeans, 1923),故針對(duì)地核的研究我們選擇頻率較低的低階簡正振型來分析.
圖3所示為本文數(shù)值模擬算法的結(jié)果與無衰減(彌散衰減或黏滯衰減)的數(shù)值模擬結(jié)果和實(shí)際地震觀測結(jié)果頻譜圖,對(duì)比所用的簡正模采用的是由一維地球模型得到的自由振蕩簡正模(圖中虛線),可作為近似參考或約束.數(shù)值模擬結(jié)果頻譜圖顯示,相較于無衰減模型,本文實(shí)驗(yàn)在地核加入了彌散阻尼后,低階的簡正振型0S2、0S3、…、0S9和3S1振型的振幅衰減明顯,模擬結(jié)果更加接近實(shí)測數(shù)據(jù).而0S10之后的振型未見振幅的衰減.這一結(jié)果與方明(1991)進(jìn)行的對(duì)球型振蕩的深度與球諧階之間的關(guān)系的數(shù)值實(shí)算結(jié)果一致.下面,我們在模擬研究地核中不同的彌散衰減對(duì)地球自由振蕩的影響時(shí),選擇0S2、0S3、…、0S9和3S1振型的頻譜來進(jìn)行對(duì)比分析.
圖3 引入地核彌散阻尼模型和無衰減模型以及實(shí)測數(shù)據(jù)(IRIS)頻譜圖對(duì)比Fig.3 Spectrum comparison of models with and without the core dispersive attenuation and IRIS data
第二組實(shí)驗(yàn)將內(nèi)核和外核彌散阻尼均設(shè)為D=-0.25 s-1,其他參數(shù)設(shè)置與上一組相同.圖4所示為不同彌散阻尼模型數(shù)值模擬實(shí)驗(yàn)、無衰減(彌散衰減或黏滯衰減)數(shù)值模擬結(jié)果和實(shí)際地震觀測結(jié)果頻譜圖對(duì)比,可看出地核彌散阻尼加大后,不影響地球自由振蕩的頻率,即本征頻率對(duì)應(yīng)性良好,只影響振幅衰減.由于地球衰減和儀器響應(yīng)等因素,實(shí)際地震臺(tái)陣記錄到的長時(shí)間傳播的地震波信號(hào)能量會(huì)大大衰減,超長周期信號(hào)較難識(shí)別,故實(shí)測數(shù)據(jù)的甚低頻段能量非常低,可以識(shí)別出的自由振蕩振型僅為0S5、0S6、0S7、0S8、0S9振型,在這些振型處本文模擬的地核自由振蕩振幅更加接近實(shí)測數(shù)據(jù).由于計(jì)算所用服務(wù)器和計(jì)算穩(wěn)定性的約束,數(shù)值模擬無法達(dá)到臺(tái)網(wǎng)地震儀的采樣頻率.實(shí)測數(shù)據(jù)來源為IRIS,地震儀器的采樣率為50 sps (sample per second),采樣間隔為0.05 s.本文方法采用的時(shí)間步長為0.1425 s.因此數(shù)值模擬結(jié)果的有效信息遠(yuǎn)少于實(shí)測數(shù)據(jù)的有效信息,信息量約相當(dāng)于實(shí)測數(shù)據(jù)的三分之一.但即使是在信息量相對(duì)較少的情況下,數(shù)值模擬的結(jié)果依然是可靠的,說明了本文算法具有高效求解全球尺度地震波長時(shí)程傳播問題及地球自由振蕩問題的優(yōu)勢.
圖4 地核彌散阻尼D=-0.25 s-1模型(藍(lán))、D=-0.1 s-1模型(黑)和無衰減模型以及實(shí)測數(shù)據(jù)頻譜圖對(duì)比Fig.4 Spectrum comparison of models with core dispersive attenuation D=-0.25 s-1 (blue),D=-0.1 s-1 (black) and model without the core dispersive attenuation and IRIS data
在此基礎(chǔ)上我們試驗(yàn)了僅增大內(nèi)核彌散阻尼,實(shí)驗(yàn)結(jié)果表明,當(dāng)內(nèi)核彌散阻尼系數(shù)增高至大于外核彌散阻尼系數(shù)時(shí),部分振型頻率和振幅偏離實(shí)際觀測數(shù)據(jù).
本小節(jié)進(jìn)行僅增大外核彌散阻尼時(shí)自由振蕩的數(shù)值實(shí)驗(yàn).給出兩組模型實(shí)驗(yàn)結(jié)果,分別為將外核彌散阻尼設(shè)為DOC=-0.9 s-1、內(nèi)核彌散阻尼設(shè)為DIC=-0.3 s-1,以及將外核彌散阻尼設(shè)為DOC=-1.1 s-1、內(nèi)核彌散阻尼設(shè)為DIC=-0.1 s-1.由圖5可見,外核彌散阻尼較大時(shí)模型模擬結(jié)果最接近實(shí)際地震觀測數(shù)據(jù), 0S5、3S1、0S6、0S7、0S8、0S9振型的模擬結(jié)果與觀測數(shù)據(jù)在本征頻率對(duì)應(yīng)性、幅值關(guān)系等方面擬合較好,說明外核彌散衰減對(duì)地球自由振蕩的影響大于內(nèi)核彌散衰減.圖6為地核敏感振型3S1、3S2頻譜分析結(jié)果,從圖中可以看出,外核彌散阻尼DOC=-1.1 s-1時(shí)3S1、3S2振型頻譜幅值衰減,更加接近實(shí)測數(shù)據(jù).表1和表2分別給出外核彌散阻尼DOC=-0.9 s-1模型和DOC=-1.1 s-1模型的數(shù)值模擬結(jié)果和實(shí)測數(shù)據(jù)的頻譜分析結(jié)果,兩個(gè)數(shù)值模型識(shí)別出的球型振蕩0S3到0S10振型本征頻率值相同,與實(shí)測頻率的誤差如表1所示;表2可以看出外核彌散阻尼DOC=-1.1 s-1模型比DOC=-0.9 s-1模型幅值均減小,更加接近實(shí)測數(shù)據(jù)幅值.在頻譜的甚低頻段,頻率的可識(shí)別度對(duì)地球自由振蕩的模擬研究至關(guān)重要.因此,本文采用的地核彌散阻尼等效模型可作為一種模擬甚低頻地球自由振蕩的有效模型.
圖5 不同地核彌散阻尼DIC=-0.3 s-1& DOC=-0.9 s-1模型(藍(lán))、DIC=-0.1 s-1& DOC=-1.1 s-1模型(紅)和無衰減模型以及實(shí)測數(shù)據(jù)頻譜圖對(duì)比Fig.5 Spectrum comparison of models with core dispersive attenuation DIC=-0.3 s-1& DOC=-0.9 s-1 (blue),DIC=-0.1 s-1& DOC=-1.1 s-1 (red) and model without the core dispersive attenuation and IRIS data
圖6 不同地核彌散阻尼DIC=-0.3 s-1& DOC=-0.9 s-1模型(紅)、DIC=-0.1 s-1& DOC=-1.1 s-1模型(藍(lán))模擬結(jié)果地核敏感振型(3S1、3S2)頻譜圖對(duì)比Fig.6 Core sensitive normal modes (3S1,3S2) spectrum comparison of models with core dispersive attenuation DIC=-0.3 s-1& DOC=-0.9 s-1 (red) and DIC=-0.1 s-1& DOC=-1.1 s-1 (blue)
表1 數(shù)值模擬結(jié)果與實(shí)測數(shù)據(jù)頻譜分析——頻率Table 1 Spectrum analysis of simulation and observation data—frequency
表2 數(shù)值模擬結(jié)果與實(shí)測數(shù)據(jù)頻譜分析——幅值Table 2 Spectrum analysis of simulation and observation data—amplitude
繼續(xù)增大外核彌散阻尼設(shè)為DOC=-1.5 s-1,內(nèi)核彌散阻尼設(shè)為DIC=-0.1 s-1, 模擬結(jié)果如圖7所示,部分振型的頻率難以分辨,說明此模型外核彌散阻尼過大,與實(shí)際地球內(nèi)部結(jié)構(gòu)不相符.
圖7 不同地核彌散阻尼DIC=-0.1 s-1& DOC=-1.5 s-1模型(藍(lán))、DIC=-0.3 s-1& DOC=-0.9 s-1模型(紅)和無衰減模型以及實(shí)測數(shù)據(jù)頻譜圖對(duì)比Fig.7 Spectrum comparison of models with core dispersive attenuation DIC=-0.1 s-1& DOC=-1.5 s-1 (blue),DIC=-0.3 s-1& DOC=-0.9 s-1 (red) and model without the core dispersive attenuation and IRIS data
綜上,采用內(nèi)核彌散阻尼為0.1 s-1、外核彌散阻尼為1.1 s-1的地核衰減模型時(shí),地球自由振蕩數(shù)值模擬結(jié)果最接近實(shí)際觀測記錄.以上數(shù)值模擬實(shí)驗(yàn)結(jié)果表明,對(duì)甚低頻段的地球自由振蕩而言,彌散衰減為地核自由振蕩衰減的重要成分.對(duì)于低頻段,主波長遠(yuǎn)遠(yuǎn)大于等效非均勻尺度, 引起較強(qiáng)的Rayleigh散射,使得等效散射衰減(彌散衰減)較黏滯衰減更為明顯;對(duì)于高頻段,由于主波長遠(yuǎn)遠(yuǎn)小于等效非均勻尺度,其散射主要為前向散射(沿傳播方向),其等效散射衰減則遠(yuǎn)弱于高頻情況下的黏滯衰減.前期數(shù)值實(shí)驗(yàn)結(jié)果表明,當(dāng)λ/a?1(Rayleigh散射模態(tài),其中λ為波長,a為等效平均自由程)或λ/a~1(Mie散射模態(tài))時(shí),彌散衰減為衰減之主要成分;反之黏滯衰減占主導(dǎo)地位(李冰非,2019).就地球自由振蕩與地球內(nèi)部結(jié)構(gòu)而言,大多為λ/a?1情形,因而彌散衰減為衰減的主要成分.
數(shù)值實(shí)驗(yàn)結(jié)果還顯示對(duì)于地球自由振蕩超低頻段外核彌散衰減的影響大于內(nèi)核彌散衰減,這種現(xiàn)象可能是由外核和內(nèi)核非均勻尺度的差異以及外核和內(nèi)核體積的差異造成.
本文基于地震臺(tái)網(wǎng)實(shí)際觀測記錄,得到了更加適用于真實(shí)地球自由振蕩數(shù)值模擬的地核等效彌散衰減系數(shù),但鑒于計(jì)算資源的限制,本文選擇的模型單元剖分不夠精細(xì),未來預(yù)計(jì)進(jìn)一步增加剖分精細(xì)程度;本文通過數(shù)值實(shí)驗(yàn)?zāi)M給出地核等效彌散阻尼,但限于模型單元剖分精細(xì)程度不夠及未考慮介質(zhì)各向異性等因素,其精度評(píng)估亦有難度.未來計(jì)劃對(duì)本文等效地核彌散衰減開展更深入的研究.
對(duì)于存在衰減效應(yīng)的地球自由振蕩問題,本文發(fā)展了具有保彌散衰減性質(zhì)、適用于低頻問題的數(shù)值模擬方法,將其應(yīng)用于全球尺度地球自由振蕩衰減效應(yīng)數(shù)值模擬,發(fā)展了基于辛格式-譜元法的保彌散數(shù)值模擬方法,提供了一種更切合實(shí)際的、適用于處理地球自由振蕩衰減問題的高精度、高效率模擬計(jì)算方法.數(shù)值模擬實(shí)驗(yàn)結(jié)果表明,在地核區(qū)域,對(duì)甚低頻段的地球自由振蕩而言,外核中的自由振蕩衰減主要源自彌散衰減.本文的研究方法提供了一種可通過自由振蕩這樣的甚低頻振蕩研究探索地核非均勻性的有效途徑.
致謝非常感謝中國科學(xué)院地質(zhì)與地球物理研究所計(jì)算模擬實(shí)驗(yàn)室為本文研究提供高性能計(jì)算集群,非常感謝IRIS提供地震數(shù)據(jù)以及CIG提供SPECFEM3D-GLOBE程序.