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

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

1 中國地震災害防御中心， 北京 100029 2 中國科學院地質與地球物理研究所， 北京 100029 3 中國地質大學地球物理與空間信息學院， 武漢 430074

0 引言

地球自由振蕩是由大地震在有限尺度的地球內部激發的三維駐波，即整體地球的固有振蕩(Benioff et al., 1961).地球自由振蕩的固有頻率和頻譜形態與地球內部結構密切相關(Lapwood and Usami, 1981)，包括地球內部橫向非均勻性及徑向非均勻性(Tsuboi et al., 1985; Tsuboi,1995)、地球內部各向異性(Woodhouse et al., 1986; Tromp, 1993,1995)、介質非完全彈性(Masters and Widmer, 1995)、以及介質對波的衰減(Davies and Bunge, 2001; Alterman et al., 1974; Pollitz et al., 1987; Park et al., 2005)等性質.由于應用駐波方法研究地球內部結構具有對采樣點要求不高的優勢，地球自由振蕩數值模擬研究可以方便、經濟地重建全球多尺度三維地球內部結構及獲得地球內部某些物理場的三維特征參數分布.

譜元法是近年來大尺度地震波及地球自由振蕩數值模擬的重要方法之一(Komatitsch et al., 2002, 2004; Komatitsch and Tromp, 2002; Tsuboi et al., 2004; Yan et al., 2014).譜元法將有限元算法和譜方法相結合，兼顧了有限元方法的靈活性和譜方法的高精度(Komatitsch and Vilotte, 1998; Liu et al., 2017a).譜元法是彈性波方程在空間近似的高階變分方法，Komatitsch和Tromp(1999)在譜元單元上采用Gauss-Lobatto-Legendre(GLL)積分并結合高階精度的拉格朗日插值基函數，得到具有對角化的質量矩陣的顯式時間域微分方程組，且數值解更好地趨近于理論解，模擬精度和計算效率得以提高；而后又將譜元方法應用于全球地震波傳播模擬研究(2002).Tromp等(2010)在多年研究的基礎上，將數值方法和高性能并行計算相結合，研制譜元法程序包SPECFEM3D_GLOBE，可以用來精細剖分地球模型、模擬三維地球模型或區域模型中地震波傳播.趙明等(2013)拓展了一種譜元-簡正振型耦合方法適用于各向異性介質中地震波傳播數值模擬.近二十年來，研究者們將超大規模高性能并行計算數值模擬方法應用于全球地震波傳播研究，將計算結果與實際觀測的超長波記錄對比研究，獲得了一系列較為令人滿意的結果(劉少林等, 2021).特別應該提到的是，嚴珍珍等(2008a)應用譜元法結合彈性波理論及高性能并行計算對由特大地震激發的地球自由振蕩進行了成功的大規模數值模擬，從而證實了這類方法的可行性，為由大地震激發地球自由振蕩的數值模擬研究奠定了基礎.

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

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

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

1.1 地球自由振蕩波動方程

地球的一般彈性運動方程為

(1)

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

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

(2)

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

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

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

(3)

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

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

(4)

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

(5)

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

(6)

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

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

2 數值實驗

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

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

2.1 長時程數值模擬實驗

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

2.2 內外核彌散阻尼相等時地球模型自由振蕩數值實驗

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

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

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

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

圖4 地核彌散阻尼D=-0.25 s-1模型(藍)、D=-0.1 s-1模型(黑)和無衰減模型以及實測數據頻譜圖對比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

在此基礎上我們試驗了僅增大內核彌散阻尼，實驗結果表明，當內核彌散阻尼系數增高至大于外核彌散阻尼系數時，部分振型頻率和振幅偏離實際觀測數據.

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

本小節進行僅增大外核彌散阻尼時自由振蕩的數值實驗.給出兩組模型實驗結果，分別為將外核彌散阻尼設為DOC=-0.9 s-1、內核彌散阻尼設為DIC=-0.3 s-1，以及將外核彌散阻尼設為DOC=-1.1 s-1、內核彌散阻尼設為DIC=-0.1 s-1.由圖5可見，外核彌散阻尼較大時模型模擬結果最接近實際地震觀測數據， 0S5、3S1、0S6、0S7、0S8、0S9振型的模擬結果與觀測數據在本征頻率對應性、幅值關系等方面擬合較好，說明外核彌散衰減對地球自由振蕩的影響大于內核彌散衰減.圖6為地核敏感振型3S1、3S2頻譜分析結果，從圖中可以看出，外核彌散阻尼DOC=-1.1 s-1時3S1、3S2振型頻譜幅值衰減，更加接近實測數據.表1和表2分別給出外核彌散阻尼DOC=-0.9 s-1模型和DOC=-1.1 s-1模型的數值模擬結果和實測數據的頻譜分析結果，兩個數值模型識別出的球型振蕩0S3到0S10振型本征頻率值相同，與實測頻率的誤差如表1所示；表2可以看出外核彌散阻尼DOC=-1.1 s-1模型比DOC=-0.9 s-1模型幅值均減小，更加接近實測數據幅值.在頻譜的甚低頻段，頻率的可識別度對地球自由振蕩的模擬研究至關重要.因此，本文采用的地核彌散阻尼等效模型可作為一種模擬甚低頻地球自由振蕩的有效模型.

圖5 不同地核彌散阻尼DIC=-0.3 s-1& DOC=-0.9 s-1模型(藍)、DIC=-0.1 s-1& DOC=-1.1 s-1模型(紅)和無衰減模型以及實測數據頻譜圖對比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模型(藍)模擬結果地核敏感振型(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 數值模擬結果與實測數據頻譜分析——頻率Table 1 Spectrum analysis of simulation and observation data—frequency

表2 數值模擬結果與實測數據頻譜分析——幅值Table 2 Spectrum analysis of simulation and observation data—amplitude

繼續增大外核彌散阻尼設為DOC=-1.5 s-1，內核彌散阻尼設為DIC=-0.1 s-1, 模擬結果如圖7所示，部分振型的頻率難以分辨，說明此模型外核彌散阻尼過大，與實際地球內部結構不相符.

圖7 不同地核彌散阻尼DIC=-0.1 s-1& DOC=-1.5 s-1模型(藍)、DIC=-0.3 s-1& DOC=-0.9 s-1模型(紅)和無衰減模型以及實測數據頻譜圖對比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

綜上，采用內核彌散阻尼為0.1 s-1、外核彌散阻尼為1.1 s-1的地核衰減模型時，地球自由振蕩數值模擬結果最接近實際觀測記錄.以上數值模擬實驗結果表明，對甚低頻段的地球自由振蕩而言，彌散衰減為地核自由振蕩衰減的重要成分.對于低頻段，主波長遠遠大于等效非均勻尺度, 引起較強的Rayleigh散射，使得等效散射衰減(彌散衰減)較黏滯衰減更為明顯；對于高頻段，由于主波長遠遠小于等效非均勻尺度，其散射主要為前向散射(沿傳播方向)，其等效散射衰減則遠弱于高頻情況下的黏滯衰減.前期數值實驗結果表明，當λ/a?1(Rayleigh散射模態，其中λ為波長，a為等效平均自由程)或λ/a～1(Mie散射模態)時，彌散衰減為衰減之主要成分；反之黏滯衰減占主導地位(李冰非，2019).就地球自由振蕩與地球內部結構而言，大多為λ/a?1情形，因而彌散衰減為衰減的主要成分.

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

3 總結

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

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

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