基于辛-譜元方法的地球自由振蕩保彌散衰減數(shù)值模擬

李冰非， 李小凡， 李峰， 龔飛

1 中國地震災(zāi)害防御中心， 北京 100029 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 中國地質(zhì)大學(xué)地球物理與空間信息學(xué)院， 武漢 430074

0 引言

地球自由振蕩是由大地震在有限尺度的地球內(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ì)對波的衰減(Davies and Bunge, 2001; Alterman et al., 1974; Pollitz et al., 1987; Park et al., 2005)等性質(zhì).由于應(yīng)用駐波方法研究地球內(nèi)部結(jié)構(gòu)具有對采樣點(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ù)，得到具有對角化的質(zhì)量矩陣的顯式時間域微分方程組，且數(shù)值解更好地趨近于理論解，模擬精度和計算效率得以提高；而后又將譜元方法應(yīng)用于全球地震波傳播模擬研究(2002).Tromp等(2010)在多年研究的基礎(chǔ)上，將數(shù)值方法和高性能并行計算相結(jié)合，研制譜元法程序包SPECFEM3D_GLOBE，可以用來精細(xì)剖分地球模型、模擬三維地球模型或區(qū)域模型中地震波傳播.趙明等(2013)拓展了一種譜元-簡正振型耦合方法適用于各向異性介質(zhì)中地震波傳播數(shù)值模擬.近二十年來，研究者們將超大規(guī)模高性能并行計算數(shù)值模擬方法應(yīng)用于全球地震波傳播研究，將計算結(jié)果與實際觀測的超長波記錄對比研究，獲得了一系列較為令人滿意的結(jié)果(劉少林等, 2021).特別應(yīng)該提到的是，嚴(yán)珍珍等(2008a)應(yīng)用譜元法結(jié)合彈性波理論及高性能并行計算對由特大地震激發(fā)的地球自由振蕩進(jìn)行了成功的大規(guī)模數(shù)值模擬，從而證實了這類方法的可行性，為由大地震激發(fā)地球自由振蕩的數(shù)值模擬研究奠定了基礎(chǔ).

盡管理論上時間域數(shù)值方法均可處理諸如地球內(nèi)部橫向非均勻性及地震激發(fā)地球自由振蕩之全過程等問題，但現(xiàn)有地球自由振蕩時間域數(shù)值模擬技術(shù)均基于非保結(jié)構(gòu)計算，其長時程計算中的高積累誤差難以避免，很難精確處理地球自由振蕩模擬這一長時程追蹤問題.近年來已有聲波及彈性波保結(jié)構(gòu)模擬方法問世(羅明秋等, 2001; Li et al., 2012; 汪文帥等, 2012; 陳世仲等, 2016; Liu et al., 2017b,c; 李冰非等, 2019)；Li 等(2012)及Gao等(2016)進(jìn)一步將逆時離散變換與三階辛相結(jié)合，獲得了長時程的低頻散數(shù)值模擬，這類保結(jié)構(gòu)方法非常適合處理有關(guān)保守動力學(xué)體系的問題.若將地球近似為彈性球體，用上述保結(jié)構(gòu)方法模擬其自由振蕩問題將是很好的選擇；然而，地球本身是一具有衰減性質(zhì)的非完全彈性球體，真實地球介質(zhì)極其復(fù)雜, 除了各向異性以外,可能還有多相多孔性以及黏性(Kuznetsov, 1997; Kustowski et al., 2008a)，其自由振蕩(駐波)問題涉及非保守動力學(xué)系統(tǒng)問題.針對這一問題，Li等(2015)發(fā)展了阻尼聲波保結(jié)構(gòu)辛算法相關(guān)理論，數(shù)理邏輯的合理性和方法的可行性均得到驗證.Cai等(2017)基于同樣的模型方程提出了另一種耗散保結(jié)構(gòu)模擬方法，對于帶耗散的地震波運(yùn)動波動方程，可以分為一個守恒系統(tǒng)和一個耗散系統(tǒng)來求解.為更精確地模擬衰減模式下的地球自由振蕩，需要發(fā)展新的、高效駐波保結(jié)構(gòu)數(shù)值方法，研究適用于衰減介質(zhì)的保彌散數(shù)值模擬方法.

本文將保結(jié)構(gòu)計算延拓至非保守動力學(xué)體系，發(fā)展了一種更切合實際的、適用于處理地球自由振蕩衰減問題的高精度、高效率模擬計算方法，以2012年印度尼西亞蘇門答臘西北海域8.6級大地震觀測記錄為數(shù)據(jù)源，進(jìn)行了衰減模式下的全球自由振蕩數(shù)值實驗，研究地球深部(如地核)彌散衰減對地球自由振蕩的影響和各種衰減模式對地球自由振蕩的影響.

1 基于辛-譜元法的地球自由振蕩保彌散衰減數(shù)值模擬

1.1 地球自由振蕩波動方程

地球的一般彈性運(yùn)動方程為

(1)

其中ρ為密度，u為質(zhì)點(diǎn)位移，ω為地球角轉(zhuǎn)動矢量，f表示外力和重力的擾動，λ，μ為Lame常數(shù)(嚴(yán)珍珍等，2008b).

為便于表達(dá)，我們對地球的一般彈性運(yùn)動方程進(jìn)行簡化，忽略大氣變化因素，忽略自轉(zhuǎn)、重力、地球橢率以及橫向不均勻性，考慮簡化的地球彈性運(yùn)動方程為：

(2)

其中F表示震源.數(shù)值模擬地球自由振蕩，本質(zhì)上也是求滿足適當(dāng)初邊值條件的地球振動的常系數(shù)偏微分方程組的數(shù)值解.

1.2 基于辛-譜元法的地球自由振蕩保彌散衰減數(shù)值模擬

用譜元法將方程組(2)離散化后可寫為：

(3)

其中M為全局質(zhì)量矩陣，K為全局剛度矩陣，U表示全局位移向量.

對于本文研究的問題，彌散衰減介質(zhì)中彈性波方程可寫為：

(4)

其中u為位移矢量，σ為應(yīng)力張量，ε為應(yīng)變張量，c為四階彈性張量，D為彌散阻尼系數(shù)，f為外力.譜元法離散空間域后的微分方程組可寫為：

(5)

將式(5)應(yīng)用阻尼項保結(jié)構(gòu)辛格式(Li et al., 2015)進(jìn)行時間離散，可得如下基于譜元法的保彌散衰減波動方程的時間格式：

(6)

其中V為速度向量. Δt為時間步長，c2，d1，d2為辛系數(shù)，U(n)和V(n)的上標(biāo)表示第n時間步，U1和V1為中間變量.

本文將采用基于(6)式的辛格式-譜元法開展衰減模式下的地球自由振蕩數(shù)值模擬研究，并探究地核中不同衰減模式及其對地球自由振蕩的影響.

2 數(shù)值實驗

本節(jié)通過數(shù)值模擬實驗驗證本文提出的方法的有效性.采用2012年4月11日印度尼西亞蘇門答臘西北海域8.6級地震(2.327°N，93.063°E)作為震源，圖1所示為本次地震事件震源位置和Global Centroid-Moment-Tensor (CMT) Project提供的8.6級地震的震源球以及全球地震臺網(wǎng)(Global Seismographic Network, GSN)公布的152個臺站所在位置(https:∥www.iris.edu/hq/programs/gsn).數(shù)值模擬使用Computational Infrastructure for Geodynamics(CIG)提供的SPECFEM3D-GLOBE譜元法程序包(經(jīng)本文保彌散-辛格式升級).

圖1 地震事件震源和全球152個臺站位置示意圖(https:∥www.iris.edu/hq/programs/gsn)Fig.1 Epicenter, 152 stations on globe and earthquake mechanism

2.1 長時程數(shù)值模擬實驗

本節(jié)分別采用保彌散-辛格式和SPECFEM3D-GLOBE譜元法程序包中所采用的Newmark算法模擬同一模型，考慮地形、橢率、重力及自轉(zhuǎn)等因素，模擬地震波傳播600 min.采用S362ANI地球模型(Kustowski et al., 2008a,b)，并將數(shù)值模擬結(jié)果進(jìn)行同參數(shù)處理.圖2所示為NNA臺站Z分量地震記錄波形對比，從圖中波形放大圖可以看出，模擬時長約388 min時Newmark算法的結(jié)果開始顯現(xiàn)出高頻數(shù)值頻散(波形曲線出現(xiàn)高頻毛刺)，本文算法的波形記錄則比Newmark算法在長時間計算后更加穩(wěn)定，顯示出在模擬長時程地震波傳播時的優(yōu)越性.采用同樣的時間步長(dt=0.1425 s)，在同樣的計算節(jié)點(diǎn)條件下，本文保彌散-辛格式模擬耗時68 h 23 m 24 s，Newmark算法模擬耗時71 h 05 m 10 s，本文保彌散-辛格式模擬長時程地震波傳播時更加穩(wěn)定，同時不增加計算量甚至計算效率略有提高.本實驗僅簡要展示低階辛格式處理大尺度地震波傳播問題的有效性及抑制數(shù)值頻散的能力，同類辛格式的數(shù)值頻散前人已有詳盡的分析(Li et al., 2015).

2.2 內(nèi)外核彌散阻尼相等時地球模型自由振蕩數(shù)值實驗

本小節(jié)設(shè)計兩組實驗，其整體均采用S362ANI地球模型，并考慮地形、橢率、重力及自轉(zhuǎn)等因素.第一組實驗采用本文保彌散衰減方法模擬外核和內(nèi)核彌散衰減模型(設(shè)內(nèi)核和外核彌散阻尼D=-0.1 s-1),和采用SPECFEM3D-GLOBE譜元法程序包自帶的Newmark算法模擬無彌散衰減的模型，模擬2012年印度尼西亞北蘇門答臘8.6級地震發(fā)生后地震波傳播600 min (dt=0.1425 s，模擬252700步).將數(shù)值模擬結(jié)果與實際地震觀測數(shù)據(jù)進(jìn)行對比分析，本文采用的實際地震觀測數(shù)據(jù)來源為IRIS(Incorporated Research Institutions for Seismology).根據(jù)地球自由振蕩的甚低頻特性，我們將兩種數(shù)值模擬算法的結(jié)果和實際地震觀測結(jié)果(IRIS下載)采用0.0039 Hz的低通濾波，選擇Chirp-z變換算法求取頻譜進(jìn)行對比，并引入地球自由振蕩理論值(簡振模)作為參考.由于外核中不存在環(huán)型振蕩，本文僅研究地核彌散阻尼對球型振蕩的影響.且隨著簡正模球諧階l的增大，自由振蕩所發(fā)生的深度越來越淺(Jeans, 1923)，故針對地核的研究我們選擇頻率較低的低階簡正振型來分析.

圖3所示為本文數(shù)值模擬算法的結(jié)果與無衰減(彌散衰減或黏滯衰減)的數(shù)值模擬結(jié)果和實際地震觀測結(jié)果頻譜圖，對比所用的簡正模采用的是由一維地球模型得到的自由振蕩簡正模(圖中虛線)，可作為近似參考或約束.數(shù)值模擬結(jié)果頻譜圖顯示，相較于無衰減模型，本文實驗在地核加入了彌散阻尼后，低階的簡正振型0S2、0S3、…、0S9和3S1振型的振幅衰減明顯，模擬結(jié)果更加接近實測數(shù)據(jù).而0S10之后的振型未見振幅的衰減.這一結(jié)果與方明(1991)進(jìn)行的對球型振蕩的深度與球諧階之間的關(guān)系的數(shù)值實算結(jié)果一致.下面，我們在模擬研究地核中不同的彌散衰減對地球自由振蕩的影響時，選擇0S2、0S3、…、0S9和3S1振型的頻譜來進(jìn)行對比分析.

圖3 引入地核彌散阻尼模型和無衰減模型以及實測數(shù)據(jù)(IRIS)頻譜圖對比Fig.3 Spectrum comparison of models with and without the core dispersive attenuation and IRIS data

第二組實驗將內(nèi)核和外核彌散阻尼均設(shè)為D=-0.25 s-1，其他參數(shù)設(shè)置與上一組相同.圖4所示為不同彌散阻尼模型數(shù)值模擬實驗、無衰減(彌散衰減或黏滯衰減)數(shù)值模擬結(jié)果和實際地震觀測結(jié)果頻譜圖對比，可看出地核彌散阻尼加大后，不影響地球自由振蕩的頻率，即本征頻率對應(yīng)性良好，只影響振幅衰減.由于地球衰減和儀器響應(yīng)等因素，實際地震臺陣記錄到的長時間傳播的地震波信號能量會大大衰減，超長周期信號較難識別，故實測數(shù)據(jù)的甚低頻段能量非常低，可以識別出的自由振蕩振型僅為0S5、0S6、0S7、0S8、0S9振型，在這些振型處本文模擬的地核自由振蕩振幅更加接近實測數(shù)據(jù).由于計算所用服務(wù)器和計算穩(wěn)定性的約束，數(shù)值模擬無法達(dá)到臺網(wǎng)地震儀的采樣頻率.實測數(shù)據(jù)來源為IRIS，地震儀器的采樣率為50 sps (sample per second)，采樣間隔為0.05 s.本文方法采用的時間步長為0.1425 s.因此數(shù)值模擬結(jié)果的有效信息遠(yuǎn)少于實測數(shù)據(jù)的有效信息，信息量約相當(dāng)于實測數(shù)據(jù)的三分之一.但即使是在信息量相對較少的情況下，數(shù)值模擬的結(jié)果依然是可靠的，說明了本文算法具有高效求解全球尺度地震波長時程傳播問題及地球自由振蕩問題的優(yōu)勢.

圖4 地核彌散阻尼D=-0.25 s-1模型(藍(lán))、D=-0.1 s-1模型(黑)和無衰減模型以及實測數(shù)據(jù)頻譜圖對比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ǔ)上我們試驗了僅增大內(nèi)核彌散阻尼，實驗結(jié)果表明，當(dāng)內(nèi)核彌散阻尼系數(shù)增高至大于外核彌散阻尼系數(shù)時，部分振型頻率和振幅偏離實際觀測數(shù)據(jù).

2.3 僅增大外核彌散阻尼時自由振蕩數(shù)值實驗

本小節(jié)進(jìn)行僅增大外核彌散阻尼時自由振蕩的數(shù)值實驗.給出兩組模型實驗結(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可見，外核彌散阻尼較大時模型模擬結(jié)果最接近實際地震觀測數(shù)據(jù)， 0S5、3S1、0S6、0S7、0S8、0S9振型的模擬結(jié)果與觀測數(shù)據(jù)在本征頻率對應(yīng)性、幅值關(guān)系等方面擬合較好，說明外核彌散衰減對地球自由振蕩的影響大于內(nèi)核彌散衰減.圖6為地核敏感振型3S1、3S2頻譜分析結(jié)果，從圖中可以看出，外核彌散阻尼DOC=-1.1 s-1時3S1、3S2振型頻譜幅值衰減，更加接近實測數(shù)據(jù).表1和表2分別給出外核彌散阻尼DOC=-0.9 s-1模型和DOC=-1.1 s-1模型的數(shù)值模擬結(jié)果和實測數(shù)據(jù)的頻譜分析結(jié)果，兩個數(shù)值模型識別出的球型振蕩0S3到0S10振型本征頻率值相同，與實測頻率的誤差如表1所示；表2可以看出外核彌散阻尼DOC=-1.1 s-1模型比DOC=-0.9 s-1模型幅值均減小，更加接近實測數(shù)據(jù)幅值.在頻譜的甚低頻段，頻率的可識別度對地球自由振蕩的模擬研究至關(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ù)據(jù)頻譜圖對比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)頻譜圖對比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ù)據(jù)頻譜分析——頻率Table 1 Spectrum analysis of simulation and observation data—frequency

表2 數(shù)值模擬結(jié)果與實測數(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所示，部分振型的頻率難以分辨，說明此模型外核彌散阻尼過大，與實際地球內(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ù)據(jù)頻譜圖對比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ù)值模擬結(jié)果最接近實際觀測記錄.以上數(shù)值模擬實驗結(jié)果表明，對甚低頻段的地球自由振蕩而言，彌散衰減為地核自由振蕩衰減的重要成分.對于低頻段，主波長遠(yuǎn)遠(yuǎn)大于等效非均勻尺度, 引起較強(qiáng)的Rayleigh散射，使得等效散射衰減(彌散衰減)較黏滯衰減更為明顯；對于高頻段，由于主波長遠(yuǎn)遠(yuǎn)小于等效非均勻尺度，其散射主要為前向散射(沿傳播方向)，其等效散射衰減則遠(yuǎn)弱于高頻情況下的黏滯衰減.前期數(shù)值實驗結(jié)果表明，當(dāng)λ/a?1(Rayleigh散射模態(tài)，其中λ為波長，a為等效平均自由程)或λ/a～1(Mie散射模態(tài))時，彌散衰減為衰減之主要成分；反之黏滯衰減占主導(dǎo)地位(李冰非，2019).就地球自由振蕩與地球內(nèi)部結(jié)構(gòu)而言，大多為λ/a?1情形，因而彌散衰減為衰減的主要成分.

數(shù)值實驗結(jié)果還顯示對于地球自由振蕩超低頻段外核彌散衰減的影響大于內(nèi)核彌散衰減，這種現(xiàn)象可能是由外核和內(nèi)核非均勻尺度的差異以及外核和內(nèi)核體積的差異造成.

3 總結(jié)

本文基于地震臺網(wǎng)實際觀測記錄，得到了更加適用于真實地球自由振蕩數(shù)值模擬的地核等效彌散衰減系數(shù)，但鑒于計算資源的限制，本文選擇的模型單元剖分不夠精細(xì)，未來預(yù)計進(jìn)一步增加剖分精細(xì)程度；本文通過數(shù)值實驗?zāi)M給出地核等效彌散阻尼，但限于模型單元剖分精細(xì)程度不夠及未考慮介質(zhì)各向異性等因素，其精度評估亦有難度.未來計劃對本文等效地核彌散衰減開展更深入的研究.

對于存在衰減效應(yīng)的地球自由振蕩問題，本文發(fā)展了具有保彌散衰減性質(zhì)、適用于低頻問題的數(shù)值模擬方法，將其應(yīng)用于全球尺度地球自由振蕩衰減效應(yīng)數(shù)值模擬，發(fā)展了基于辛格式-譜元法的保彌散數(shù)值模擬方法，提供了一種更切合實際的、適用于處理地球自由振蕩衰減問題的高精度、高效率模擬計算方法.數(shù)值模擬實驗結(jié)果表明，在地核區(qū)域，對甚低頻段的地球自由振蕩而言，外核中的自由振蕩衰減主要源自彌散衰減.本文的研究方法提供了一種可通過自由振蕩這樣的甚低頻振蕩研究探索地核非均勻性的有效途徑.

致謝非常感謝中國科學(xué)院地質(zhì)與地球物理研究所計算模擬實驗室為本文研究提供高性能計算集群，非常感謝IRIS提供地震數(shù)據(jù)以及CIG提供SPECFEM3D-GLOBE程序．