地球自由振蕩的保結(jié)構(gòu)模擬

陳世仲 ， 李小凡， 汪文帥

1 中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院地球與行星物理重點(diǎn)實(shí)驗(yàn)室， 北京　100029 2 中國科學(xué)院大學(xué)， 北京　100049 3 華北水利水電大學(xué)資源與環(huán)境學(xué)院， 鄭州　450045 4 寧夏大學(xué)數(shù)學(xué)計(jì)算機(jī)學(xué)院， 銀川　750021

地球自由振蕩的保結(jié)構(gòu)模擬

陳世仲1,2,3， 李小凡1， 汪文帥4

1 中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院地球與行星物理重點(diǎn)實(shí)驗(yàn)室， 北京100029 2 中國科學(xué)院大學(xué)， 北京100049 3 華北水利水電大學(xué)資源與環(huán)境學(xué)院， 鄭州450045 4 寧夏大學(xué)數(shù)學(xué)計(jì)算機(jī)學(xué)院， 銀川750021

摘要大地震等諸多激勵(lì)均能激發(fā)全球自由振蕩現(xiàn)象，通常表現(xiàn)為駐波形式的全球整體振蕩.現(xiàn)有的地震波數(shù)值模擬方法多為非保結(jié)構(gòu)方法，無法壓制長時(shí)程計(jì)算中的積累誤差.本文采用優(yōu)化的三階辛格式譜元法，對(duì)地球自由振蕩及全球尺度的地震波傳播進(jìn)行了長時(shí)程模擬.通過與傳統(tǒng)的基于Newmark算法的譜元方法結(jié)果對(duì)比分析，明確驗(yàn)證了本文所得優(yōu)化的三階辛格式譜元法在模擬地球自由振蕩等大規(guī)模長時(shí)程問題上的優(yōu)越性和準(zhǔn)確性.上述進(jìn)展在方法論層面為今后探測、刻畫全球尺度地球非均勻結(jié)構(gòu)的駐波數(shù)值方法奠定了部分基礎(chǔ)，并為相關(guān)研究領(lǐng)域提供了新的選擇.

關(guān)鍵詞地球自由振蕩； 保結(jié)構(gòu)模擬； 長時(shí)程追蹤； 三階辛格式譜元法

1引言

大地震、火山爆發(fā)、大氣流動(dòng)等地球內(nèi)部運(yùn)動(dòng)均能激發(fā)地球作為類球體的本征振蕩.地球自由振蕩現(xiàn)象從基本特征來看，是在整個(gè)地球尺度上擾動(dòng)的駐波(Bolt, 2012).地球自由振蕩信號(hào)包含了地球內(nèi)部物質(zhì)的物理特性，包括地球內(nèi)部非均勻結(jié)構(gòu)、地球自身物理狀態(tài)等諸多重要的地球物理信息.因此，精確模擬、分析地球自由振蕩是解讀地球內(nèi)部物質(zhì)結(jié)構(gòu)和特征的關(guān)鍵手段之一.

1829年，Poisson最早研究了完全彈性固體球的振動(dòng)問題(Shearer, 2009)對(duì)地球自由振蕩的認(rèn)識(shí)是從解析方法研究開始的，比如球諧展開.1882年，Lamb詳細(xì)討論了均勻球體模型的自由振蕩規(guī)律(Lapwood and Usami, 1981).1911年，Love在考慮重力作用的情況下，探討了可壓縮球體自由振蕩現(xiàn)象(Bolt, 2012).

地球自由振蕩頻譜特征受多種地球結(jié)構(gòu)因素及狀態(tài)的影響，包括：地球介質(zhì)的橫向非均勻性及徑向非均勻性、各向異性、介質(zhì)能量衰減效應(yīng)、地球的自轉(zhuǎn)、地球的橢率、地形以及斷裂機(jī)制等因素.因此研究地球的本征振蕩對(duì)研究地球內(nèi)部結(jié)構(gòu)及物理性質(zhì)具有重要意義.

鑒于目前相關(guān)研究的匱乏，本文從波動(dòng)理論出發(fā)，探尋求解地球自由振蕩問題的方法，并依托大規(guī)模并行計(jì)算集群付諸實(shí)踐.通過哈密爾頓體系運(yùn)動(dòng)方程的推導(dǎo)，結(jié)合譜元法，得出一個(gè)新的求解波動(dòng)方程的三級(jí)三階辛格式譜元算法(汪文帥等, 2012b)，并將其應(yīng)用于地球自由振蕩及全球尺度地震波的長時(shí)程模擬.數(shù)值模擬結(jié)果與實(shí)際觀測資料的對(duì)比分析表明，相對(duì)于采用Newmark算法的傳統(tǒng)譜元方法，本文提出的算法具有良好的保結(jié)構(gòu)特性，非常適合地球自由振蕩及全球尺度地震波的長時(shí)程模擬這類有高精度長時(shí)程求解變系數(shù)偏微分方程需求的物理問題(陳世仲等, 2012).本文的實(shí)驗(yàn)及分析還進(jìn)一步證實(shí)，地球內(nèi)部大尺度非均勻結(jié)構(gòu)對(duì)地球自由振蕩有著明顯的影響(陳世仲等, 2013).

2地球自由振蕩求解理論

2.1地球自由振蕩波動(dòng)方程建立

Alterman等(1959)將地球自由振蕩滿足的二階偏微分方程本征值問題歸化為一階常微分方程組，奠定了現(xiàn)代計(jì)算機(jī)數(shù)值求解的基礎(chǔ).Dahlen(1968)采用“二階擾動(dòng)”理論研究了自轉(zhuǎn)、橢球地球的簡正模.方俊分別采用1066A模型和Harold Jeffreys的修正地球模型進(jìn)行了環(huán)型和球型自由振蕩簡正模的計(jì)算，給出了1066A模型的環(huán)型和球型自由振蕩本征周期較完整的理論計(jì)算值，并討論了地球內(nèi)部密度參數(shù)和彈性參數(shù)的變化對(duì)本征周期的影響(方俊, 1990a; 1990b).Tsuboi等(1985)結(jié)合擾動(dòng)技術(shù)和反演方法，利用牛頓方法及變分原理求解非線性方程，理論上分析了橫向不均勻性地球自由振蕩特征值的偏微分方程模型問題，并得出比一階退化擾動(dòng)理論更精確的特征頻率和特征函數(shù).隨后，Tsuboi(1995)又根據(jù)變分原理計(jì)算了橫向不均勻性和非彈性地球模型的基本模式，研究了地球自由振蕩譜峰與橫向不均勻性以及非彈性等之間的可能關(guān)系，即譜峰的中心頻率隨時(shí)間而改變.嚴(yán)珍珍等采用譜元法與并行計(jì)算數(shù)值模擬特大地震激發(fā)的彈性波在地球內(nèi)部傳播過程(Yan et al., 2014; 嚴(yán)珍珍等, 2008; 嚴(yán)珍珍等, 2012).

地震波的波動(dòng)方程數(shù)值模擬實(shí)質(zhì)上是求彈性地震波波動(dòng)方程的數(shù)值解，模擬的地震波場包含了地震波產(chǎn)生和傳播過程的所有運(yùn)動(dòng)學(xué)與動(dòng)力學(xué)信息，所以波動(dòng)方程數(shù)值模擬方法一直在地震模擬中占有重要地位.

類似地震波的波動(dòng)方程數(shù)值模擬，數(shù)值模擬地球自由振蕩，本質(zhì)上也是求滿足適當(dāng)初邊值條件的地球振動(dòng)的常系數(shù)偏微分方程組的數(shù)值解.地球的一般彈性運(yùn)動(dòng)方程為

(1)

其中，f表示為外力和重力的擾動(dòng)，一般情況大氣的變化因素不予考慮，則只考慮重力的影響.

據(jù)此，可對(duì)地球的一般彈性運(yùn)動(dòng)方程進(jìn)行簡化，以便與數(shù)值模擬求解.通常采取的簡化包括： (1) 忽略大氣因素； (2) 忽略科里奧利力； (3) 忽略重力的影響； (4) 忽略地球橢率，簡化為球體； (5) 忽略橫向不均勻性，簡化為圈層結(jié)構(gòu)；等.

可得到，無震源激發(fā)情況下，簡化的地球彈性運(yùn)動(dòng)方程：

(2)

其中，λ,μ只與r有關(guān).

欲求解式(2)，還需另加近似條件和邊界條件：

(2) 在地心r=0處正則，即u,v,w=0；

(3) 在圈層的固-固邊界上，應(yīng)力和位移都連續(xù)；

2.2波動(dòng)方程理論求解

由震源項(xiàng)F激發(fā)的波動(dòng)，在不考慮衰減、重力和科里奧利力的前提下，可將地球的彈性波動(dòng)方程簡化為

(3)

運(yùn)用亥姆霍茲定理求解，可將F和u表示為

(4)

(5)

將(4)式、(5)式代入方程(3)，化簡得到

(6)

求解(6)，得出

(7)

其中，

(8)

針對(duì)不同的震源機(jī)制F，求解積分，通過球諧函數(shù)展開，最終解出u.Lapwood和Usami(1981)曾討論過f(t)=exp(iωt)和f(t)為單位階躍函數(shù)時(shí)的解的表達(dá)式.亦可使用傅里葉變換，把方程從時(shí)間域變換到頻域中求解.

2.3三階辛格式譜元數(shù)值求解波動(dòng)方程

目前常用的波動(dòng)方程數(shù)值模擬方法多為非保結(jié)構(gòu)方法,如：有限差分法、偽譜法、時(shí)域非辛格式有限元法及時(shí)域非辛格式譜元法等.盡管這些方法簡潔、易于實(shí)現(xiàn)、執(zhí)行效率高，但存在固有的、數(shù)理邏輯上的局限性.尤其是其時(shí)間域離散化近似的非保結(jié)構(gòu)性質(zhì)，在一定程度和一定范圍內(nèi)對(duì)波動(dòng)方程及解的結(jié)構(gòu)破壞是無法避免的，其數(shù)值解對(duì)于理論解的明顯偏離也是不可避免的.這一先天不足對(duì)于短時(shí)程計(jì)算問題并不明顯，但對(duì)于地球自由振蕩這類長時(shí)程數(shù)值模擬則極為明顯.而辛算法則可有效解決此類問題，理論研究表明辛算法具有極佳的保結(jié)構(gòu)特性.在長時(shí)間模擬的情況下，辛算法依然能夠非常有效地控制累積誤差的增長.因此，辛算法非常適合需要進(jìn)行長時(shí)程模擬并對(duì)精度有較高要求的物理學(xué)數(shù)值模擬問題(Komatitsch and Tromp, 1999).

對(duì)于波動(dòng)方程空間域離散，譜元法的優(yōu)越性毋庸置疑.該方法兼具有限元處理邊界和結(jié)構(gòu)的靈活性和偽譜法的高精度和快速收斂特性等優(yōu)勢，在單元上采用 Gauss-Lobatto-Legendre(GLL)積分并結(jié)合基于Legendre正交多項(xiàng)式的高階拉格朗日型插值，使得質(zhì)量矩陣對(duì)角化，避免了單元內(nèi)高階插值所帶來的Runge現(xiàn)象.由于可以直接采用顯格式算法，避開了大規(guī)模線性方程組的并行求解，大大簡化了數(shù)值計(jì)算的計(jì)算量和實(shí)現(xiàn)方法，同時(shí)提高了模擬精度和算法的穩(wěn)定性.

根據(jù)譜元法的原理，將波動(dòng)方程在空間域離散，得到的微分方程組可寫為

(9)

將波動(dòng)方程空間域離散后，在時(shí)間域可使用多種算法來進(jìn)行離散，常用的非辛算法有差分算法、龍格庫塔算法等.這些非辛算法普遍存在精度較低，長時(shí)程計(jì)算后數(shù)值頻散明顯.目前，在譜元法數(shù)值模擬研究領(lǐng)域，較為認(rèn)可的時(shí)間域離散算法是Newmark算法，其從本質(zhì)上來說，屬于二階保動(dòng)量算法.理論研究表明，該方法在短時(shí)程的數(shù)值模擬中，體現(xiàn)出尚佳的性能，但該算法對(duì)積累誤差和相位誤差的控制一般，隨著模擬時(shí)長的增加，二者均有比較明顯的增加(Li et al., 2012; Zampieri and Pavarino, 2006).

針對(duì)長時(shí)程高精度模擬的物理問題的需求，為更高效地獲得更高精度的模擬結(jié)果，本文優(yōu)化了Li等(2012)提出的新三階辛算法(New Three-stage Third-order Symplectical Algorithm, NTSTO)，采用其進(jìn)行時(shí)間域離散，并將其結(jié)合進(jìn)行空間域離散的譜元法，進(jìn)而得到一種新的具有時(shí)-空保結(jié)構(gòu)特性的三階辛格式譜元法.

對(duì)于微分方程數(shù)值求解，可進(jìn)行如下三級(jí)三階保辛算法：

(10)

將式(10)代入式(9)，同時(shí)考慮到譜元法所得的質(zhì)量矩陣M為對(duì)角矩陣，非常便于線性方程求解，提高計(jì)算效率，從而可得到波動(dòng)方程的三階辛格式求解格式：

(11)

其中，Δt為時(shí)間步長，而c(i)、d(i)為辛條件系數(shù)(馮康和秦孟兆, 2003)，i=1,2,3.

由于NTSTO算法具有良好的辛算法保結(jié)構(gòu)特性，因此在進(jìn)行時(shí)間域離散時(shí)，可采取相對(duì)大的時(shí)間步長計(jì)算.如此，即可保證模擬結(jié)果的準(zhǔn)確性，同時(shí)又能彌補(bǔ)三階算法比較耗時(shí)的不足.

而關(guān)于辛條件系數(shù)的選取，可根據(jù)不同的最優(yōu)化原則進(jìn)行.汪文帥等(2012a)采用相位最小誤差原理，求取另一組辛條件系數(shù)，并通過多組數(shù)值試驗(yàn)，驗(yàn)證了NTSTO算法無論在內(nèi)存消耗、穩(wěn)定性及計(jì)算耗時(shí)，還是長時(shí)程跟蹤能力都相當(dāng)優(yōu)越，而且可以有效地處理復(fù)雜幾何形狀和復(fù)雜介質(zhì).但其探討的模型是物性參數(shù)和結(jié)構(gòu)特點(diǎn)相對(duì)簡單的二維模型，且模擬地震波時(shí)長較短(與地球自由振蕩模擬相比而言)，其選擇的辛條件系數(shù)足以達(dá)到研究要求.

然而，地球自由振蕩研究所涉及的是具有復(fù)雜物性變化及橫向非均勻結(jié)構(gòu)的三維地球模型，且模擬地震波時(shí)程極長.如果僅考慮相位最小誤差，數(shù)理邏輯角度將很難滿足研究的需要.因此，本文根據(jù)最小均方誤差原理，求取得到一組優(yōu)化的辛條件系數(shù)：

(12)

代入此優(yōu)化的條件系數(shù)，即得到優(yōu)化的三階辛譜元方法(NewThree-stageThird-orderSymplecticalSpectralElementMethod,NTSTO-SEM)，并將其應(yīng)用于后續(xù)的地球自由振蕩深入研究.

3地球自由振蕩數(shù)值模擬及分析

3.1震源及模型的選取

據(jù)地球自由振蕩的特性可知，只有大地震才能產(chǎn)生比較明顯的地球自由振蕩現(xiàn)象.因此在進(jìn)行自由振蕩研究時(shí)，選取8級(jí)以上的地震作為震源.本文選取2012年4月11日16時(shí)38分(北京時(shí)間)在印度尼西亞蘇門答臘西北海域(2.35°N、92.82°E)發(fā)生的8.6級(jí)強(qiáng)震(Magnitude 8.6-OFF THE WEST COAST OF NORTHERN SUMATRA. USGS于11 April 2012發(fā)布).

CMT公布的震源機(jī)制解為

(13)

圖1所示為本次地震震中和全球168個(gè)臺(tái)站所在位置以及CMT提供的本次地震的震源球.

文中模擬采用Preliminary Reference Earth Model(PREM)地球模型.該模型是Dziewonski和Anderson(1981)提出的一維參考模型，引入黏滯衰減和各向異性的平均值，因此其上地幔部分是隨頻率變化的橫向各向同性模型.

3.2數(shù)值模擬方案

本次數(shù)值模擬使用了294個(gè)進(jìn)程進(jìn)行并行計(jì)算.參考全球地震臺(tái)網(wǎng)(Global Seismographic Network, GSN)和中國地震臺(tái)網(wǎng)(China Seismographic Network, CSN)公布的臺(tái)站信息，在全球168個(gè)臺(tái)站所處經(jīng)緯位置分別設(shè)置地震波接收點(diǎn).使用的譜元程序?yàn)镃omputational Infrastructure for Geodynamics(CIG)提供的SPECFEM3D_GLOBE(Komatitsch and Tromp, 1999; Komatitsch et al., 1999).

基于PREM模型，考慮了地形和橢率，分別使用傳統(tǒng)譜元法和NTSTO-SEM算法各模擬了地震波在全球表面?zhèn)鞑?00 min的整個(gè)過程.其中，在采用NTSTO-SEM算法計(jì)算時(shí)，為降低NTSTO算法的耗時(shí)，經(jīng)過穩(wěn)定性分析得出，NTSTO算法可取的最大時(shí)間步長為Newmark時(shí)間步長的2.8倍.詳細(xì)地球模型參數(shù)、數(shù)值模擬參數(shù)及網(wǎng)格化參數(shù)如表1所示.

表1　數(shù)值模擬相關(guān)參數(shù)

注：*為GLL點(diǎn)距和譜單元大小均為地表平均值; ☆為NTSTO算法時(shí)間步長采用Newmark算法時(shí)間步長的2倍.

從理論角度分析，采用相同時(shí)間步長，三階算法耗時(shí)約為二階算法的1.5倍.本文在數(shù)值模擬時(shí)，NTSTO算法采用大時(shí)間步長.在相同震源、相同模型的情況下，采用Newmark離散的譜元法和NTSTO-SEM算法(下文簡稱Newmark算法和NTSTO算法)，分別模擬地震波在全球傳播500 min，所消耗的時(shí)間分別為463h41min5s、486h32min34s，用時(shí)比約為1∶1.04.該比例遠(yuǎn)低于1∶1.5的理論值，表明本文推導(dǎo)的三階辛格式譜元法在采用大時(shí)間步長運(yùn)算時(shí)，可非常有效地彌補(bǔ)三階算法耗時(shí)的不足.

3.3數(shù)值實(shí)驗(yàn)結(jié)果分析

根據(jù)本次地震以水平錯(cuò)動(dòng)為主要表現(xiàn)的特點(diǎn)，同時(shí)考慮研究橫向非均勻的需要，選取了分布于地球表面不同位置的168個(gè)臺(tái)站作為分析對(duì)象.文中分析的實(shí)際地震觀測數(shù)據(jù)，均是按照模擬采取的震源和傳播時(shí)長，從Incorporated Research Institutions for Seismology(IRIS)申請所得.限于篇幅，只列舉其中2個(gè)特征臺(tái)站CHTO.IU(98.94°E, 18.81°N)和PALK.II(80.70°E, 7.27°N)，將按照時(shí)間域波形和頻率域頻譜兩類進(jìn)行分析.

(1) 時(shí)間域波形曲線分析

圖2所示為CHTO.IU臺(tái)站模擬地震波Z分量波形對(duì)比圖，圖3所示為PALK.II臺(tái)站模擬地震波Z分量波形對(duì)比圖.圖2、圖3中，紅色曲線為臺(tái)網(wǎng)觀測數(shù)據(jù)，藍(lán)色曲線為數(shù)值模擬的結(jié)果.辛算法的優(yōu)勢在于長時(shí)程模擬控制誤差，為突出說明此特點(diǎn)，圖4和圖5分別選取圖2和圖3的時(shí)長200～500 min時(shí)段進(jìn)行展示.

由圖2和圖3可見，在0～20 min的初始階段，兩種算法所得的波形曲線與實(shí)測數(shù)據(jù)的波形契合度均較好，而NTSTO算法的結(jié)果明顯好于Newmark算法.但隨著模擬時(shí)間的推移，與實(shí)測數(shù)據(jù)相比，兩種算法所得地震波曲線均出現(xiàn)較明顯的偏差.這是因?yàn)閿?shù)值模擬采用的模型是PREM模型，不考慮衰減和橫向非均勻等影響因素，模擬的結(jié)果與實(shí)際的觀測數(shù)據(jù)必然存在差異.盡管如此，隨著模擬時(shí)長的增長，如200 min之后，NTSTO算法所得的波形曲線與觀測數(shù)據(jù)仍具有較好的對(duì)應(yīng)性，并且保持穩(wěn)定.從圖4和圖5中即可看出.

圖1　震中和全球168個(gè)臺(tái)站位置及震源球Fig.1　Epicenter, 168 stations on globe and earthquake mechanism

圖2　CHTO.IU臺(tái)站地震波Z分量波形對(duì)比(a) 譜元法; (b) 三階辛方法.Fig.2　Z component of seismic wave of station CHTO.IU(a) SEM; (b) NTSTO-SEM.

另外，對(duì)波形的頻譜分析表明，模擬結(jié)果與觀測數(shù)據(jù)存在相位差.在進(jìn)行時(shí)域誤差分析時(shí)，模擬結(jié)果與觀測數(shù)據(jù)不能直接做定量分析，特別是幅值的差異.因?yàn)檫@類對(duì)比所出現(xiàn)的“誤差”很大程度上來自模型與實(shí)際地球間的巨大差異，而目前階段這一“誤差”很難定量地從總誤差中予以區(qū)分.因此，今后的研究應(yīng)不斷完善地球模型，使其更貼近實(shí)際地球介質(zhì)構(gòu)造.

(2) 頻率域頻譜分析

因?yàn)樽杂烧袷幨潜菊黝l率問題，研究的主要是頻率域問題.圖6所示為CHTO.IU臺(tái)站模擬地震波Z分量頻譜對(duì)比圖，圖7所示為PALK.II臺(tái)站模擬地震波Z分量頻譜對(duì)比圖.圖6和圖7中，紅色曲線為臺(tái)網(wǎng)觀測數(shù)據(jù)頻譜分析結(jié)果，藍(lán)色曲線為數(shù)值模擬的頻譜分析結(jié)果.

圖3　PALK.II臺(tái)站地震波Z分量波形對(duì)比(a) 譜元法; (b) 三階辛方法.Fig.3　Z component of seismic wave of station PALK.II(a) SEM; (b) NTSTO-SEM.

圖4　CHTO.IU臺(tái)站地震波Z分量波形長時(shí)程對(duì)比(a) 譜元法; (b) 三階辛方法.Fig.4　Z component of seismic wave of station CHTO.IU(a) SEM; (b) NTSTO-SEM.

不同振型的地球自由振蕩的理論頻率值通常被當(dāng)作參照依據(jù)，故本文亦引用之.圖6和圖7中，上方橫坐標(biāo)軸分別為不同階的球型振蕩(從0S2到0S16)，連接上下橫坐標(biāo)軸的黑色虛線為各階球型振蕩對(duì)應(yīng)的理論頻率.理論研究表明，將無數(shù)個(gè)不同自由振蕩振型進(jìn)行疊加，即可得出地球振動(dòng)波形(Bolt, 2012; Lapwood and Usami, 1981).但這些理論值是通過一維簡化地球模型求解所得，此模型同地球?qū)嶋H的三維非均勻復(fù)雜結(jié)構(gòu)(包括橫向非均勻結(jié)構(gòu))有極大偏差，所以即使將這些理論值用于簡正模求和，也無法得到真實(shí)可信的振動(dòng)波形.故這些理論值對(duì)實(shí)際研究而言，僅可作為粗略的一級(jí)近似參考或約束.從圖6及圖7中也可發(fā)現(xiàn)，這些理論值大多與觀測數(shù)據(jù)相去甚遠(yuǎn).

圖5　PALK.II臺(tái)站地震波Z分量波形長時(shí)程對(duì)比(a) 譜元法; (b) 三階辛方法.Fig.5　Z component of seismic wave of station PALK.II(a) SEM; (b) NTSTO-SEM.

圖6　CHTO.IU臺(tái)站地震波Z分量頻譜對(duì)比(a) 譜元法; (b) 三階辛方法.Fig.6　Z component of spectrum of station CHTO.IU(a) SEM; (b) NTSTO-SEM.

圖7　PALK.II臺(tái)站地震波Z分量頻譜對(duì)比(a) 譜元法; (b) 三階辛方法.Fig.7　Z component of spectrum of station PALK.II(a) SEM; (b) NTSTO-SEM.

表2所示為CHTO.IU臺(tái)站模擬結(jié)果和臺(tái)網(wǎng)數(shù)據(jù)的頻譜分析結(jié)果，表3所示為PALK.II臺(tái)站模擬結(jié)果和臺(tái)網(wǎng)數(shù)據(jù)的頻譜分析結(jié)果.兩表中，第一列為臺(tái)網(wǎng)數(shù)據(jù)識(shí)別出的譜峰頻率值，第二列為Newmark算法模擬識(shí)別出的譜峰頻率值，第三列為Newmark算法與觀測數(shù)據(jù)的絕對(duì)誤差；第四列為NTSTO算法模擬識(shí)別出的譜峰頻率值，第五列為NTSTO算法與觀測數(shù)據(jù)的絕對(duì)誤差.

本文采用的是PREM模型，所以兩個(gè)臺(tái)站識(shí)別的振型頻率值應(yīng)基本一致.而實(shí)際地球是復(fù)雜的橫向非均勻衰減結(jié)構(gòu)，會(huì)造成兩個(gè)臺(tái)站觀測數(shù)據(jù)識(shí)別的振型頻率值略有差異或者部分缺失.因此，綜合兩個(gè)臺(tái)站的觀測數(shù)據(jù)，可得到最全的振型頻率參考值.表2和表3第一列的加括號(hào)數(shù)據(jù)即兩個(gè)臺(tái)站觀測數(shù)據(jù)的互補(bǔ).

表2和表3共38組對(duì)比數(shù)據(jù)，只有1組數(shù)據(jù)，NTSTO算法的結(jié)果稍遜于Newmark算法.這還建立在前者的數(shù)據(jù)量明顯少于后者以及觀測數(shù)據(jù)的基礎(chǔ)上，因?yàn)榍罢卟捎昧舜髸r(shí)間步長.如果計(jì)算條件允許，采用同樣小步長，NTSTO算法的精度將遠(yuǎn)高于Newmark算法.

對(duì)比圖6和表2、圖7和表3，可明確看出，NTSTO算法所得結(jié)果在頻譜的峰值對(duì)應(yīng)性、幅值關(guān)系、能量分布等多方面都明顯優(yōu)于Newmark算法所得結(jié)果.尤其是在甚低頻段，頻率可識(shí)別度明顯高于Newmark算法的結(jié)果，這點(diǎn)對(duì)于自由振蕩的研究至關(guān)重要，見表2和表3的前兩組數(shù)據(jù).自由振蕩的最長周期可達(dá)54 min，頻率約為0.29 mHz.鑒于介質(zhì)衰減等因素，長時(shí)間傳播的情況下，信號(hào)的能量會(huì)大幅降低，超長周期信號(hào)很難有效識(shí)別.因此在圖中，實(shí)測數(shù)據(jù)的甚低頻段的頻譜能量都非常低.

表2　CHTO.IU臺(tái)站模擬結(jié)果和臺(tái)網(wǎng)數(shù)據(jù)的

表3　PALK.II臺(tái)站模擬結(jié)果和臺(tái)網(wǎng)數(shù)據(jù)的

縱觀時(shí)間域和頻率域的結(jié)果，模擬結(jié)果與實(shí)際觀測資料均未能完全契合，主要原因是實(shí)際地球的內(nèi)部結(jié)構(gòu)及其介質(zhì)的物理性質(zhì)遠(yuǎn)比PREM模型復(fù)雜得多，表明地球內(nèi)部大尺度非均勻結(jié)構(gòu)對(duì)地球自由振蕩有著明顯的影響.另外，信息量的缺失也是不可忽視的原因.由于計(jì)算平臺(tái)限制及計(jì)算穩(wěn)定性的約束，數(shù)值模擬無法達(dá)到臺(tái)網(wǎng)地震儀的采樣頻率.從表1可知，Newmark算法采取的最小時(shí)間步長為0.07315 s，NTSTO算法采用前者的2.8倍的大步長，最小時(shí)間步長僅為0.20482 s，而此次下載實(shí)測數(shù)據(jù)的地震儀器的采樣率為50 sps，即采樣間隔為0.05 s.因此，實(shí)測數(shù)據(jù)的有效信息遠(yuǎn)多于模擬結(jié)果的有效信息.其中，NTSTO算法的信息量僅相當(dāng)于Newmark算法的約1/3，實(shí)測數(shù)據(jù)的約1/4.

盡管如此，NTSTO算法的結(jié)果依然與觀測數(shù)據(jù)有良好的擬合關(guān)系，能夠更準(zhǔn)確識(shí)別自由振蕩振型頻率值.模擬的結(jié)果表明本文推導(dǎo)的三階辛格式譜元法在采用大時(shí)間步長運(yùn)算時(shí)，精度依然明顯高于采用Newmark算法的傳統(tǒng)譜元法.如果計(jì)算平臺(tái)滿足，采用與Newmark算法相同的小時(shí)間步長，模擬結(jié)果的準(zhǔn)確性將極大提升，該算法的優(yōu)越性將更加明顯.同時(shí)，將其結(jié)合觀測數(shù)據(jù)也可大幅增高觀測記錄中振型頻率的可識(shí)別率，并且逐步完善對(duì)地球內(nèi)部結(jié)構(gòu)的認(rèn)知.

4結(jié)論

前人在數(shù)值模擬自由振蕩方面的研究多集中于同自由振蕩不同振型的理論值進(jìn)行對(duì)比.前述已知，這些理論值所基于的模型與實(shí)際地球結(jié)構(gòu)出入甚大.數(shù)值模擬結(jié)果與理論參考值對(duì)應(yīng)好壞實(shí)際意義一般，應(yīng)更關(guān)注同實(shí)際觀測資料的對(duì)比分析.

本文提供了一種保結(jié)構(gòu)模擬地球自由振蕩的新方法，并應(yīng)用此方法對(duì)2012年4月11日印度尼西亞蘇門答臘8.6級(jí)地震激發(fā)的地球自由振蕩進(jìn)行了數(shù)值模擬.數(shù)值實(shí)驗(yàn)結(jié)果顯示，無論是時(shí)間域理論地震記錄，還是頻率域波譜曲線分布，均與實(shí)際觀測數(shù)據(jù)具有良好的對(duì)應(yīng).研究分析及數(shù)值結(jié)果證實(shí)，本文推出的方法具有保結(jié)構(gòu)性質(zhì)，可有效控制波動(dòng)方程長時(shí)程模擬中的積累誤差及相位誤差.盡管數(shù)值結(jié)果與觀測數(shù)據(jù)尚有差別，但這些差別主要源于波動(dòng)方程空間離散時(shí)單元剖分較粗糙(限于計(jì)算平臺(tái)的條件)及所采用地球模型與實(shí)際地球結(jié)構(gòu)間的差異.

本文所涉及的工作將保結(jié)構(gòu)計(jì)算(模擬)拓展至地球自由振蕩研究領(lǐng)域，其方法適合有高精度長時(shí)程求解變系數(shù)偏微分方程需求的物理問題；在方法論層面上為探測、刻畫全球尺度地球非均勻(特別是橫向非均勻)結(jié)構(gòu)的駐波數(shù)值方法奠定了部分基礎(chǔ)，并為相關(guān)研究領(lǐng)域提供了新的選擇.

致謝非常感謝“中國科學(xué)院地質(zhì)與地球物理研究所計(jì)算模擬實(shí)驗(yàn)室(Computer Simulation Lab，IGGCAS)”為本文研究提供高性能計(jì)算集群.同時(shí)，非常感謝IRIS提供地震數(shù)據(jù)以及CIG提供SPECFEM3D_GLOBE程序.

References

Alterman Z, Jarosch H, Pekeris C L. 1959. Oscillations of the Earth.Proc.Roy.Soc.Lond.Ser.-A, 252(1268): 80-95.

Bolt B A. 2012. Seismology: Surface Waves and Earth Oscillations. Burlington, MA: Elsevier.

Bullen K E, Haddon R A. 1967. Derivation of an Earth model from free oscillation data.Proc.Natl.Acad.Sci.USA, 58(3): 846-852. Chen S Z, Li X F, Wang W S, et al. 2012. Simulation of the Earth′s free oscillation based on a lateral inhomogeneous earth model.∥ 28th Ann. Mtg., Chinese Geophysical Society (in Chinese). Beijing: Chinese Geophysical Society.

Chen S Z, Li X F, Wang W S, et al. 2013. The Earth′s Free Oscillation of Structure-preserving Simulation in Analysis of Lateral Heterogeneity of the Earth Structure on the Global Scale.∥ 29th Ann. Mtg., Chinese Geophysical Society (in Chinese). Beijing: Chinese Geophysical Society.

Dahlen F A. 1968. The normal modes of a rotating, elliptical earth.Geophys.J.Int., 16(4): 329-367.

Dziewonski A M, Anderson D L. 1981. Preliminary reference Earth model.Phys.EarthPlan.Int., 25(4): 297-356.

Fang J. 1990a. Parameter kenerls in linear inversion of the Earth lree oscillation (I).ScienceinChinaSeriesB(in Chinese), 1990, (2): 202-207.

Fang J. 1990b. Parameter kenerls in linear inversion of the earth free oscillation (II).ScienceinChinaSeriesB(in Chinese), (3): 329-336.

Feng K, Qin M Z. 2003. Symplectic Geometric Algorithms for Hamiltonian Systems (in Chinese). Hangzhou: Zhejiang Science and Technology Press.

Komatitsch D, Tromp J. 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation.GeophysicalJournalInternational, 139(3): 806-822. Komatitsch D, Vilotte J P, Vai R, et al. 1999. The spectral element method for elastic wave equations—Application to 2-D and 3-D seismic problems.Int.J.Numer.Meth.Eng., 45(9): 1139-1164. Komatitsch D, Tromp J. 2002. Spectral-element simulations of global seismic wave propagation—I. Validation.GeophysicalJournalInternational, 149(2): 390-412.

Lapwood E R, Usami T. 1981. Free Oscillations of the Earth. New York: Cambridge University Press.

Li X F, Wang W S, Lu M W, et al. 2012. Structure-preserving modelling of elastic waves: a symplectic discrete singular convolution differentiator method.GeophysicalJournalInternational, 188(3): 1382-1392.Lillie R J. 1999. Whole Earth Geophysics: An Introductory Textbook for Geologists and Geophysicists. New Jersey: Prentice Hall.Shearer P M. 2009. Introduction to Seismology. 2nd ed. New York: Cambridge University Press.Tsuboi S, Geller R J, Morris S P. 1985. Partial derivatives of the eigenfrequencies of a laterally heterogeneous earth model.Geophys.Res.Lett., 12(12): 817-820.

Tsuboi S. 1995. Free oscillations of a laterally heterogenous and anelastic earth.PureAppl.Geophys., 145(3-4): 445-457.

Wang W S, Li X F, Lu M W, et al. 2012a. Structure-preserving modeling for seismic wavefields based upon a multisymplectic spectral element method.ChineseJ.Geophys. (in Chinese), 55(10): 3427-3439, doi: 10.6038/j.issn.0001-5733.2012.10.026.

Wang W S, Li X F, Chen S Z, et al. 2012b. Spectral-element method for seismic wave equations in a 3-D spherical coordinate system.∥ 28th Ann. Mtg., Chinese Geophysical Society (in Chinese). Beijing: Chinese Geophysical Society.

Yan Z Z, Zhang H, Yang C C, et al. 2008. An initial study of the numerical simulation of the Earth′s free oscillations process excited by earthquake.AdvancesinEarthScience(in Chinese), 23(10): 1020-1026.

Yan Z Z, Zhang H, Fan X T, et al. 2012. Comparative analysis on the characteristics of low-frequency energy released by the Wenchuan earthquake and Kunlun Mountain earthquake.ChineseJ.Geophys. (in Chinese), 55(12): 4218-4230, doi: 10.6038/j.issn.0001-5733.2012.12.033.

Yan Z Z, Zhang H, Fan X T, et al. 2014. Comparative analysis of the characteristics of global seismic wave propagation excited by theMw9.0 Tohoku earthquake using numerical simulation.ScienceChinaEarthSciences, 57(7): 1626-1636.

Zampieri E, Pavarino L F. 2006. Approximation of acoustic waves by explicit Newmark′s schemes and spectral element methods.JournalofComputationalandAppliedMathematics, 185(2): 308-325.

附中文參考文獻(xiàn)

陳世仲, 李小凡, 汪文帥等. 2012. 基于橫向非均勻地球模型的地球自由振蕩模擬.∥中國地球物理學(xué)會(huì)第二十八屆年會(huì). 北京: 中國地球物理學(xué)會(huì).

陳世仲, 李小凡, 汪文帥等. 2013. 全球尺度地球內(nèi)部橫向非均勻結(jié)構(gòu)研究中的地球自由振蕩保結(jié)構(gòu)數(shù)值模擬.∥中國地球物理學(xué)會(huì)第二十九屆年會(huì). 北京: 中國地球物理學(xué)會(huì).

方俊. 1990a. 地球自由振蕩線性反演中的參數(shù)核(Ⅰ). 中國科學(xué)(B輯 化學(xué) 生命科學(xué) 地學(xué)), (2): 202-207.

方俊. 1990b. 地球自由振蕩線性反演的參數(shù)核(Ⅱ). 中國科學(xué)(B輯 化學(xué) 生命科學(xué) 地學(xué)), (3): 329-336.

馮康, 秦孟兆. 2003. 哈密爾頓系統(tǒng)的辛幾何算法. 杭州: 浙江科學(xué)技術(shù)出版社.

汪文帥, 李小凡, 魯明文等. 2012a. 基于多辛結(jié)構(gòu)譜元法的保結(jié)構(gòu)地震波場模擬. 地球物理學(xué)報(bào), 55(10): 3427-3439, doi: 10.6038/j.issn.0001-5733.2012.10.026.

汪文帥, 李小凡, 陳世仲等. 2012b. 三維球坐標(biāo)系下的地震波方程的譜元解法.∥中國地球物理學(xué)會(huì)第二十八屆年會(huì). 北京: 中國地球物理學(xué)會(huì).

嚴(yán)珍珍, 張懷, 楊長春等. 2008. 地震激發(fā)地球自由振蕩過程的數(shù)值模擬初步探索. 地球科學(xué)進(jìn)展, 23(10): 1020-1026.

嚴(yán)珍珍, 張懷, 范湘濤等. 2012. 汶川與昆侖山強(qiáng)烈地震激發(fā)的地球自由振蕩頻譜的對(duì)比分析. 地球物理學(xué)報(bào), 55(12): 4218-4230, doi: 10.6038/j.issn.0001-5733.2012.12.033.

(本文編輯汪海英)

Structure-preserving numerical simulation for Earth′s free oscillations

CHEN Shi-Zhong1,2,3, LI Xiao-Fan1, WANG Wen-Shuai4

1KeyLaboratoryofEarthandPlanetaryPhysics,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcademyofSciences,Beijing100049,China3SchoolofResourcesandEnvironment,NorthChinaUniversityofWaterResourcesandElectricPower,Zhengzhou450045,China4SchoolofMathematicsandComputerScience,NingxiaUniversity,Yinchuan750021,China

AbstractFree oscillations of the Earth are usually as standing waves of entire globe, which can be excited by many factors such as a great earthquake. The oscillation types of this physical phenomenon depend on Earth′s structure and physical parameters. Therefore, intrinsic characters of the Earth can be revealed by studying Earth′ free oscillations. Numerical simulation is an effective way to research seismic waves travelling in different earth models. In essence, the numerical simulation of Earth′s free oscillations is a long-term tracing of variable coefficient wave equation. The majority of numerical methods for seismic wave simulation are non-structure-preserving methods, which cannot suppress the accumulative error in long-time computation,e.g.,finite difference methods, pseudo-spectral methods, and finite element methods.In comparison, the symplectical methods based on the Hamilton system possess fine structure-preserving ability which makes them be suitable for the physical problem acquiring good precision in long-time simulation. The spectral element method combines the generality of the finite element method with the accuracy of spectral techniques. This method expands the solution in trigonometric series, of which a chief advantage is that the resulting method is of very high order. It also chooses high-degree piecewise polynomial basis functions, also achieving a very high order of accuracy. Symplectical methods introduced into spectral element methods may be an accurate and effective choice. In this research, a three-order symplectical spectral element method is modified for long-time simulation of Earth′s free oscillations and seismic waves propagating in the global range. The event simulated is a great earthquake (Mw8.6) off the west coast of northern Sumatra on 11 April 2012. The PREM earth model is used for simulation. Our method and traditional spectrum element method using Newmark time scheme are chosen for simulation respectively with the same parameters and computation scales. The time step adopted by our method is 2.8 times that adopted by the Newmark method. Even under this extreme condition, the feasibility and effectiveness of this method have been proved by analysis of the simulation results, especially, by comparing with real seismic data processed with the same procedures and parameters. Those data are provided by the Incorporated Research Institutions for Seismology. This improvement can be a new and valid alternative to numerical methods for investigation on heterogeneity (especially lateral heterogeneity) in the globe range. Meanwhile, Earth′s free oscillations influenced obviously by large inhomogeneous structure of Earth interiors has been confirmed by comparison of simulation results and real seismic data recorded by seismic networks.KeywordsEarth′s free oscillations; Structure-preserving simulation; Long-term tracing; Three-order symplectical spectral element method

基金項(xiàng)目國家自然科學(xué)基金項(xiàng)目(41174047,40874024,41204041)資助.

作者簡介陳世仲，男，1985年生．中國科學(xué)院地質(zhì)與地球物理研究所在讀博士研究生，主要從事地震波正反演及地球自由振蕩等相關(guān)研究．E-mail: csz042@mail.iggcas.ac.cn

doi:10.6038/cjg20160513 中圖分類號(hào)P315

收稿日期2015-04-20，2016-03-04收修定稿

陳世仲， 李小凡， 汪文帥. 2016. 地球自由振蕩的保結(jié)構(gòu)模擬.地球物理學(xué)報(bào)，59(5):1685-1695,doi:10.6038/cjg20160513.

Chen S Z, Li X F, Wang W S. 2016. Structure-preserving numerical simulation for Earth′s free oscillations.ChineseJ.Geophys. (in Chinese),59(5):1685-1695,doi:10.6038/cjg20160513.