基于譜元－簡(jiǎn)正振型耦合方法的核幔邊界D″區(qū)地震波波形模擬方法研究

趙 明，趙 亮，Yann Capdeville

1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所 巖石圈演化國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029

2 中國(guó)科學(xué)院大學(xué)，北京 100049

3 Laboratoire de planétologie et de géodynamique de Nantes，CNRS，F(xiàn)rance

1 引 言

地球核幔邊界D″區(qū)是一個(gè)構(gòu)造復(fù)雜和強(qiáng)烈變化的區(qū)域，其中最重要的特征之一便是具有顯著的地震波各向異性．通過對(duì)地震各向異性的觀測(cè)研究，可以獲取地球內(nèi)部構(gòu)造變形、動(dòng)力學(xué)過程等信息．目前，地震學(xué)上主要采用橫波分裂觀測(cè)來研究D″區(qū)各向異性．比較常用的方法是測(cè)量地幔射線路徑相似的橫波震相如SKS－SKKS、S－ScS或Sdiff波的相對(duì)到時(shí)，這樣可以校正臺(tái)站端和震源端上地幔各向異性的影響，經(jīng)過校正后的剪切波分裂時(shí)差δt和快軸方向φ可以認(rèn)為代表了D″各向異性的大小和方向［2－3］．大多數(shù)相對(duì)走時(shí)分析的觀測(cè)結(jié)果通常與垂直對(duì)稱軸橫向各向同性（簡(jiǎn)稱VTI）介質(zhì)結(jié)構(gòu)相符［4－6］．

受地震和接收臺(tái)站分布的影響，全球尺度D″區(qū)各向異性的地震學(xué)觀測(cè)覆蓋率和分辨率目前還很有限，存在的問題主要集中在兩方面：有限的射線方位覆蓋和上地幔各向異性校正的誤差．因此一方面需要布設(shè)大規(guī)模密集寬頻帶數(shù)字地震臺(tái)陣，另一方面還需要繼續(xù)發(fā)展新的地震成像和計(jì)算方法提高觀測(cè)的分辨率．特別是其中一個(gè)重要的研究方法是將正演模擬與波形分析相結(jié)合，以便更好地約束各向異性的幾何形狀，減少分裂時(shí)間分析中模型的不確定性［7］．對(duì)于邊界形態(tài)、物質(zhì)屬性極為復(fù)雜的D″區(qū)而言，需要發(fā)展能夠處理三維非均勻、各向異性介質(zhì)復(fù)雜模型的高精度地震波波形模擬方法．

目前有著廣泛應(yīng)用的波動(dòng)方程數(shù)值解法主要有：有限差分法［8］、譜方法、有限元法以及譜元法等．這些方法各有其優(yōu)缺點(diǎn)：有限差分法以差分代替微分，算法簡(jiǎn)單，效率高，但其應(yīng)用受到難以處理自由邊界和復(fù)雜結(jié)構(gòu)模型的限制；偽譜法在空間域基于快速Fourier變換或Chebyshev變換，精度高，空間頻散小，收斂速度快，但是同樣也只能處理簡(jiǎn)單模型結(jié)構(gòu)和邊界；有限元法雖然可以模擬任意復(fù)雜介質(zhì)模型和適應(yīng)自由邊界，但其低階方法存在嚴(yán)重?cái)?shù)值頻散，高階方法存在偽波的問題，同時(shí)計(jì)算量很大；而譜元法是一種融合有限元法和偽譜法思想的方法，它保留了有限元處理邊界和復(fù)雜結(jié)構(gòu)的靈活性和偽譜法精度高并且快速收斂的優(yōu)點(diǎn)，因此對(duì)于全球地震波模擬而言是一種比較理想的數(shù)值方法．譜元法經(jīng)過近十幾年的發(fā)展，其對(duì)復(fù)雜三維介質(zhì)模型的模擬精度已經(jīng)得到很好的檢驗(yàn)，可以綜合考慮各向異性、彈性衰減、速度和密度非均勻性、地球自轉(zhuǎn)、重力、橢球率、地表起伏和海洋等因素［9－11］，被證明是一種在地震波全球和區(qū)域數(shù)值模擬中應(yīng)用前景非常廣闊的方法，其主要缺點(diǎn)在于計(jì)算量大，對(duì)計(jì)算資源的要求相當(dāng)高．

為了充分地利用譜元法靈活處理三維復(fù)雜模型區(qū)域的優(yōu)勢(shì)，同時(shí)在保證計(jì)算結(jié)果精度的前提下盡可能控制運(yùn)算量，Capdeville等提出一種將譜元法和經(jīng)典解析解方法簡(jiǎn)正振型法耦合的方法——譜元－簡(jiǎn)正振型法［12］，同時(shí)設(shè)計(jì)了專門針對(duì)D″區(qū)展開研究的“sandwich”模型架構(gòu)（圖1），即只對(duì)核幔邊界以上約300km 范圍使用譜元法，對(duì)于D″區(qū)上下兩個(gè)球?qū)ΨQ區(qū)域采用計(jì)算需求更低的簡(jiǎn)正振型法計(jì)算［13］．兩種方法采用DtN 算子耦合，其主要作用是：對(duì)給定的位移邊界條件（又稱Dirichlet邊界條件），通過“DtN”算子可以計(jì)算對(duì)應(yīng)邊界處的牽引力（也稱Neumann邊界條件）［14］．耦合方法的計(jì)算效率和精度在“地球模擬器”上得到了很好的證實(shí)：對(duì)于各向同性PREM 模型計(jì)算大于8s周期的體波波形，配備128 個(gè)CPU 核、160GB 內(nèi) 存，計(jì) 算2000s理論地震圖需約28～29h．經(jīng)過測(cè)試，對(duì)比全局的譜元方法，如SPECFEM3D對(duì)相同模型計(jì)算大于9s周期2000s長(zhǎng)理論地震圖，最低需要占用384個(gè)CPU核，538GB內(nèi)存，模擬時(shí)間約16h．并且，隨著模擬頻率的增高，全局譜元法的計(jì)算需求還會(huì)呈指數(shù)級(jí)增長(zhǎng)．盡管如此，目前Capdeville等的方法還只局限于針對(duì)非均勻各向同性介質(zhì)，因此將該方法發(fā)展為能夠模擬地震波在地球核幔邊界D″地區(qū)各向異性介質(zhì)中傳播的數(shù)值方法具有十分重要的科學(xué)和實(shí)用意義．

2 問題描述和方法介紹

2．1 地震波波動(dòng)方程

對(duì)初始處于靜力學(xué)平衡狀態(tài)的三維地球模型Ω，忽略地球自轉(zhuǎn)，其地震波運(yùn)動(dòng)方程可寫成如下形式：

ρ（r）為密度，隨深度變化；r表示地球半徑；u為位移場(chǎng)，其對(duì)時(shí)間的二階導(dǎo)數(shù)記為ü；f代表震源，分點(diǎn)力源和矩張量震源的等效體力兩種，后者可表示為：f（r，t）＝－M?rδ（r－r0）g（t－t0），其中M為矩張量，g（t－t0）為震源時(shí)間函數(shù)．

方程（1）中Κ定義為彈性－重力算符，包含彈性項(xiàng)和重力項(xiàng)，忽略地震波傳播引起的地球質(zhì)量重新分布效應(yīng)［15］（又稱柯林近似）后得到：

右端第一項(xiàng)為彈性項(xiàng)，T（u）表示彈性張量；右端第二項(xiàng)表示靜力學(xué)平衡狀態(tài)下只由重力引起的壓力；右端第三項(xiàng)表示密度擾動(dòng)引起的重力變化．

對(duì)于線彈性介質(zhì)，（2）式中彈性張量T（u），根據(jù)胡克定律，有

其中Cijkl（i，j，k，l＝1，2，3）為廣義彈性系數(shù)矩陣，該矩陣具有多種對(duì)稱性，Cijkl＝Cjikl＝Cijlk＝Cklij．對(duì)于最一般的彈性各向異性介質(zhì)，用21個(gè)獨(dú)立彈性系數(shù)即可描述．此外，采用張量變換符號(hào)，可將四階張量Cijkl變?yōu)槎ACij（i，j＝1，2…，6）．

對(duì)于VTI介質(zhì)，彈性系數(shù)矩陣為：

其中只有C11，C31，C33，C44，C66是獨(dú)立參數(shù)，且與地震波速存在如下關(guān)系［16］：

實(shí)際地球內(nèi)的黏彈性介質(zhì)會(huì)引起地震波衰減，Komatitsch等［17］用若干標(biāo)準(zhǔn)線彈性體的組合近似模擬地震波衰減效應(yīng)，譜元－簡(jiǎn)正振型法沿用這一方案．

假定模擬開始時(shí)（t＝0）介質(zhì)處于靜止?fàn)顟B(tài)，則初始條件為：

同時(shí)地表應(yīng)滿足自由邊界條件：

此外，在不同計(jì)算區(qū)域的邊界處需滿足應(yīng)力和位移連續(xù)條件（n表示外法向單位向量）．

固固邊界Γ1（Ω1下邊界，ΩD″上邊界）：

固液邊界Γ2（Ω2上邊界，ΩD″下邊界）：

2．2 空間域波動(dòng)方程“弱”形式及求解簡(jiǎn)述

對(duì)于核幔邊界D″區(qū)（記為ΩD″）采用譜元法計(jì)算，譜元法的求解基于方程（1）的“弱”形式：

w為任意測(cè)試向量；Γ－i（i＝1，2）代表ΩD″的上、下耦合邊界，（·，·）為內(nèi)積形式，例如

“弱”形式的初始條件變?yōu)椋?/p>

方程（12）的譜元法求解需要將模型區(qū)域ΩD″劃分成一系列互不重疊單元Ωe，e＝1，2，…，n，且ΩD″＝Ωe．由于譜元法求解的單元是規(guī)則四邊形或立方體，對(duì)于球殼狀ΩD″，需要先建立“立方－球”映射，使得球殼可以分解為六個(gè)相同的立方體，再對(duì)每個(gè)立方體進(jìn)一步劃分成三維六面體網(wǎng)格計(jì)算［18］．接著對(duì)每個(gè)單元函數(shù)采用高階Lagrange插值多項(xiàng)式近似，單元積分采用Gauss－Lobatto－Legendre（GLL）積分，積分點(diǎn)與插值點(diǎn)相同，從而可以形成對(duì)角質(zhì)量矩陣簡(jiǎn)化運(yùn)算［9－10］．

2．3 頻率域求解

對(duì)D″區(qū)以上和以下的球?qū)ΨQ區(qū)域Ω1（3850km＜r＜6371km）、Ω2（0＜r＜3480km），需要用簡(jiǎn)正振型法求解如下頻率域波動(dòng)方程：

根據(jù)Phinney等［19］，方程在球?qū)ΨQ模型和給定位移邊界條件下有解并能用廣義球諧基函數(shù)γl，m（θ，φ）（l為角階數(shù)，m為方位階數(shù)）展開．考慮固固邊界條件（8）的情形，其解的形式如下：

先不考慮方程（17）的右端項(xiàng)，將（18）代入（17），對(duì)每一個(gè)給定的角階數(shù)l和頻率ω，都能得到三個(gè)相互獨(dú)立且符合初始正則性條件的解，其中兩個(gè)為球振型，一個(gè)為扭振型，記為｛qdl（r，ω），q＝1，2，3｝．若方程（17）加上右端項(xiàng)，相應(yīng)地可得到一個(gè)的特解，記為（r，ω）．此時(shí)方程（17）的解可寫為：

系數(shù)｛qal，m（ω），m∈［－l，l］｝的值隨耦合邊界條件（8）而定．

固液邊界條件（10）、（11）下的求解過程與固固邊界情形基本類似，不同之處在于只有一個(gè)滿足條件的解dl（r，ω）．

2．4 邊界耦合

在譜元法中對(duì)（13）式，即ΩD″邊界的牽引力（i＝1，2）的計(jì)算需要用到“DtN”算子．“DtN”算子首先在頻率域中建立，以固固邊界耦合為例，對(duì)每一個(gè)角階數(shù)l和頻率ω，邊界位移為（rΓ，ω），利用（19）式及連續(xù)性條件（11）都可得到對(duì)應(yīng)的頻率域“DtN”算子lD＋i（ω）及邊界牽引力（rΓ，ω）值．

“DtN”算子在頻率－波數(shù)域和時(shí)間－空間域建立的詳細(xì)推導(dǎo)過程見Capdeville等［12］的文章．對(duì)位于頻率域Ω1的震源和臺(tái)站的影響，需要對(duì)“DtN”算子D＋1和地表位移引入額外的算子加以考慮，具體方法見Capdeville等［13］的附錄．

2．5 并行執(zhí)行

為得到具有足夠分辨率的經(jīng)過D″區(qū)的SKS、SKKS、S、ScS波形，觀察其隨模型改動(dòng)產(chǎn)生的變化，本研究模擬的最短周期為6s，理論地震圖計(jì)算時(shí)間為2000s，時(shí)間步長(zhǎng)0．08s．為此對(duì)ΩD″劃分9×9×6×6共2916 個(gè)六面體單元，單元采用6～8 階Lagrange多項(xiàng)式插值．計(jì)算在中科院地球深部結(jié)構(gòu)重點(diǎn)實(shí)驗(yàn)室基于OPENMPI并行技術(shù)的地球模擬器上調(diào)用128個(gè)CPU 核進(jìn)行，約需36h．相比于譜元法，譜元簡(jiǎn)正振型方法由于只對(duì)核幔邊界以上300km 左右區(qū)域使用譜元法，大大降低了數(shù)值模擬所需要?jiǎng)澐值膯卧獢?shù)量和計(jì)算資源使用量．

本研究所用譜元簡(jiǎn)正振型法為hsemm＿classic5．1［20］，并根據(jù)本研究需要拓展為能夠模擬VTI介質(zhì)．此外，對(duì)程序在mpi通信、網(wǎng)格生成速度方面分別做了優(yōu)化，使得程序在地球模擬器上運(yùn)行更為流暢．

3 方法驗(yàn)證

對(duì)表1所示的PREM 各向同性D″區(qū)模型（其中Cij值由（5）式推出），譜元簡(jiǎn)正振型耦合方法與簡(jiǎn)正振型法已經(jīng)做過對(duì)比驗(yàn)證，二者計(jì)算結(jié)果相當(dāng)一致［21］．

表1 PREM 模型D″區(qū)參數(shù)Table 1 Parameters of PREM model′s D″area

為了證明本研究對(duì)hsemm＿classic5．1各向異性拓展的有效性，在此基礎(chǔ)上設(shè)置了如表2 所示各向異性參數(shù)ξ＝≈1．3的VTI－PREM 模型，并使用拓展方法與簡(jiǎn)正振型法進(jìn)行對(duì)比驗(yàn)證．

表2 VTI－PREM 模型D″區(qū)參數(shù)Table 2 Parameters of model VTI－PREM model′s D″area

由于對(duì)小于8s周期的球振型理論地震圖簡(jiǎn)正振型法計(jì)算結(jié)果尚不夠理想［22］，本文只對(duì)比兩種方法在8s時(shí)的計(jì)算結(jié)果．對(duì)比所用震源時(shí)間函數(shù)如圖1a所示，其中拐角頻率ωc≈0．13Hz（圖1b）決定了模擬的最短周期為8s．計(jì)算結(jié)果如圖2所示，在震中距為65°時(shí)，三分量上二者波形幾乎沒有差別，當(dāng)震中距為125°時(shí)，噪音會(huì)相應(yīng)增加，但不影響震相的匹配程度．值得一提的是從算法角度說，本文采用的PREM 和VTI－PREM 模型可以看做徑向分層、橫向均勻的三維模型，其在hsemm＿classic中的代入處理上與三維模型一致，因此驗(yàn)證該方法對(duì)一維模型的計(jì)算精度實(shí)際上也驗(yàn)證了該方法對(duì)三維模型的計(jì)算精度［23］．

圖1 對(duì)比所用的Heavis震源時(shí)間函數(shù)（a）及其振幅譜（b），拐角頻率0．13Hz（約8s）．Fig．1 （a）The Heavis source－time function used for benchmark between CSEM（coupled spectral element method with mode）and NMS（normal mode summation）method；（b）The corresponding source spectrum to the left．The central frequency is 0．07Hz and the corner frequency is 0．13Hz（～8s）．

圖2 簡(jiǎn)正振型法（紅）和耦合CSEM 方法（黑）分別對(duì)VTI－PREM 模型計(jì)算得到的三分量理論地震圖對(duì)比．（a）（b）（c）臺(tái)站震中距為65°，三分量震相到時(shí)和振幅上符合很好；（d）（e）（f）臺(tái)站震中距為125°，隨著震中距增大，CSEM 方法的噪音增大，但主要震相符合得很好Fig．2 Comparison between NMS（red）synthetics and CSEM synthetics（black）in three components（Z，R，T）．（a）（b）（c）Epicenter distance 65°，the waveforms match perfectly in both arrival time and amplitude；（d）（e）（f）Epicenter distance 125°，the noise of CSEM result increases，but the amplitude and arrival time still match very well

4 一維各向異性模型中的S、ScS、SKS波形模擬

本研究應(yīng)用VTI－PREM 模型，同時(shí)選取震源深度分別為312．6km（模擬S－ScS）和612．4km（模擬SKS－SKKS）的兩個(gè)地震事件（圖3），并采用更簡(jiǎn)單的Gauss震源時(shí)間函數(shù)（圖4），探討D″區(qū)在VTI各向異性情形下S，ScS，SKS、SKKS 波的傳播和變化．

圖3 （a）事件1和事件2，震源深度分別為612．4km 和312．6km，臺(tái)站震中距范圍分別為60°～130°和65°～90°，背景為CMB附近S波層析成像圖［24］；（b）S，ScS，SKS二維走時(shí)路徑圖［25］（圖中數(shù)字為震中距，單位°）Fig．3 （a）A map of two events used for waveform modelling．The depths of the two events are 612．4km and 312．6km，seperately，and the epicenter distances range from 60°～130°and 65°～90°，seperately；（b）The ray paths of S（red），ScS（purple），SKS（green），SKKS（blue）．

圖4 各向異性模擬所用的高斯震源時(shí)間函數(shù)（a）及震幅譜（b），拐角頻率為0．125Hz（8s）Fig．4 （a）Gauss source－time function used for the anisotropic modelling．The source－time function is similar to a triangle with half－duration of 8s；（b）Spectral amplitude of the source，the corner frequency is 0．125Hz．

4．1 S，ScS震相

在VTI－PREM 模型中，預(yù)設(shè)波在VTI各向異性介質(zhì)中平均SV 速度（6．265km·s－1），在核幔邊界以上150km 范圍比SH 速度（7．265km·s－1）慢1km·s－1，ξ≈1．3．模擬結(jié)果印證了這一點(diǎn)：如圖5，S震相的徑向和橫向分量在震中距范圍65°～90°到時(shí)和波形基本一致，這主要是因?yàn)镾震相沒有通過各向異性D″區(qū)；而ScS的橫向分量比徑向分量在震中距65°～74°及82°～95°都更早出現(xiàn)，這反映了ScS震相通過各向異性D″區(qū)引起的S波分裂效應(yīng)．

圖5 VTI－PREM 模型計(jì)算得到的理論地震圖SV（紅虛）與SH（黑）位移分量對(duì)比，按S參考到時(shí)排列并對(duì)波形自歸一，截取S震相以后100s．震中距范圍65°～90°，共18臺(tái)站Fig．5 Synthetic displacement seismograms of SH（black）and SV （red）components for the VTI－PREM model，aligned by S arrival predicted by PREM and each seismogram self－normalized by its maximum amplitude．The epicenter distances of the 18stations range from 65°～90°．

4．2 SKS和SKKS震相

從圖6（a、c）看，VTI－PREM 模型計(jì)算的理論地震圖SKS、SKKS徑向分量相比各向同性情形的SKS（ref）、SKKS（ref）都發(fā)生了延遲，而橫向分量上（圖6（b、d））則沒有明顯的SKS、SKKS震相出現(xiàn)．這主要有兩種可能的原因：一是由于SKS、SKKS且快速近垂直穿過D″區(qū)，這么短的路徑和時(shí)間里很難發(fā)生顯著的分裂現(xiàn)象；二是SKS、SKKS均為外核內(nèi)的P波在核幔邊界轉(zhuǎn)換成的SV波，而在VTI介質(zhì)情形下，徑向偏振的S 波幾乎不會(huì)發(fā)生分裂．如果要討論SKS 和SKKS的橫向分量及變化，還需要考慮介質(zhì)的方位角各向異性進(jìn)行模擬，這將在今后下一步的工作中完成．

圖6 VTI－PREM 模型計(jì)算得到的理論地震圖徑向分量和橫向分量：（a）（b）按SKS參考到時(shí)排列并對(duì)波形自歸一，截取SKS震相以后100s，震中距范圍90°～120°，共10臺(tái)站；（c）（d）按SKKS參考到時(shí)排列并對(duì)波形自歸一，截取SKKS震相以后100s，震中距范圍90°～130°，共14臺(tái)站Fig．6 Synthetic displacement seismograms of radial components and transverse components for the VTI－PREM model：（a）（b）aligned by SKS arrival predicted by PREM and each seismogram self－normalized by its maximum amplitude．The epicenter distances of the 10stations range from 90°～120°；（c）（d）aligned by SKKS arrival predicted by PREM and each seismogram self－normalized by its max amplitude．The epicenter distances of the 14stations range from 90°～130°

5 討 論

本研究擴(kuò)展了Capdeville等提出的并行譜元－簡(jiǎn)正振型方法，使其具備了模擬垂直對(duì)稱軸各向異性介質(zhì)的功能，與經(jīng)典解析解方法簡(jiǎn)正振型法的對(duì)比證明了該擴(kuò)展的有效性．應(yīng)用新程序，模擬了VTI－PREM 模型下的SV 和SH 波形變化，所得到的分裂結(jié)果正確反映了介質(zhì)結(jié)構(gòu)的變化．其中對(duì)SKS波的模擬顯示，VTI介質(zhì)情形下的D″區(qū)會(huì)引起SKS波在徑向分量上的快慢變化，而不會(huì)引起橫向分量上的變化，如果要觀測(cè)到SKS和SKKS的橫向分量變化，需要進(jìn)一步模擬方位角各向異性介質(zhì)結(jié)構(gòu)下（如TTI）的S波分裂現(xiàn)象．

目前D″區(qū)實(shí)際觀測(cè)認(rèn)為該區(qū)在核幔邊界以上約300km 范圍S波各向異性強(qiáng)度總體而言dlnξ＜3%［26］，從數(shù)值模擬角度而言這屬于較“弱”的各向異性，需要計(jì)算的地震波最短周期為1～2s，是一項(xiàng)頗具挑戰(zhàn)性的工作．其最大的困難主要來自兩方面：一是簡(jiǎn)正振型方法在球振型計(jì)算上的局限，其低于7～8s周期的結(jié)果會(huì)出現(xiàn)很多噪音，不過近幾年簡(jiǎn)正振型法方法上和理論上都有很大改進(jìn)，為譜元簡(jiǎn)正振型耦合方法對(duì)更高頻率的應(yīng)用提供了新的發(fā)展空間［27－28］；二是計(jì)算量相當(dāng)大，全球模擬需要用到大規(guī)模并行計(jì)算，對(duì)計(jì)算機(jī)軟硬件投入很大，目前尚不具備普遍開展此項(xiàng)研究的條件．單獨(dú)從方法的角度說，譜元法對(duì)于模擬的最短周期并沒有任何限制，因此隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，譜元法必將在核幔邊界各向異性模擬中發(fā)揮更大的作用．

6 結(jié) 論

本文將一種用于模擬核幔邊界D″區(qū)非均勻性模擬的數(shù)值模擬方法——譜元簡(jiǎn)正振型法擴(kuò)展為各向異性，并通過與簡(jiǎn)正振型方法的對(duì)比驗(yàn)證，證明譜元簡(jiǎn)正振型方法對(duì)于＞8s頻率的理論地震圖計(jì)算模擬結(jié)果是可靠的．同時(shí)對(duì)VTI介質(zhì)結(jié)構(gòu)模型的計(jì)算模擬，表明譜元－簡(jiǎn)振振型法可以模擬各向異性介質(zhì)中的體波波形變化．下一步的工作將比較方位角各向異性介質(zhì)中SV 和SH 波形的不同，并將模擬提升到更高頻率，以應(yīng)對(duì)精確模擬S波分裂現(xiàn)象的需求．

致 謝 感謝何玉梅博士提供的畫圖程序以及與羅玉來進(jìn)行的有益探討；感謝中國(guó)科學(xué)院地質(zhì)與地球物理研究所深部結(jié)構(gòu)重點(diǎn)實(shí)驗(yàn)室提供Earth 模擬器；感謝Earth模擬器負(fù)責(zé)人張志剛博士在并行計(jì)算中提供的幫助．

（References）

［1］ Dziewonski A M，Anderson D L．Preliminary reference Earth model．Phys．EarthPlanet．Inter．，1981，25（4）：297－356．

［2］ Thomas C，Kendall J M．The lowermost mantle beneath northern Asia－II．Evidence for lower－mantle anisotropy．Geophys．J．Int．，2002，151（1）：296－308．

［3］ Wookey J，Kendall J M，Rümpker G．Lowermost mantle anisotropy beneath the north Pacific from differential S－ScS splitting．GeophysicalJournalInternational，2005，161（3）：829－838．

［4］ Nowacki A，Wookey J，Kendall J M．New advances in using seismic anisotropy，mineral physics and geodynamics to understand deformation in the lowermost mantle．J．Geod．，2011，52（3－4）：205－228．

［5］ 楊鳳琴，劉斌，倪四道等．西伯利亞下地幔底部的剪切波各向異性．地震學(xué)報(bào)，2008，30（2）：209－213．

Yang F Q，Liu B，Ni S D，et al．Shear velocity anisotropy of the lowermost mantle beneath the Siberia．ActaSeismologica Sinica（in Chinese），2008，30（2）：209－213．

［6］ 戴志陽(yáng)，劉斌，王宵翔等．西太平洋下D″區(qū)的剪切波各向異性．地震學(xué)報(bào)，2007，29（5）：459－466．

Dai Z Y，Liu B，Wang X X，et al．Shear wave anisotropy in D"region beneath the western Pacific．ActaSeismologica Sinica（in Chinese），2007，29（5）：459－466．

［7］ 戴志陽(yáng)．西太平洋下地幔D"層的地震波速度各向異性研究［博士論文］．合肥：中國(guó)科學(xué)技術(shù)大學(xué)，2008．

Dai Z Y．Seismic velocity anisotropy in D′layer of mantle beneath the Western Pacific（in Chinese）［Ph．D．Thesis］．Hefei：University of Science and Technology of China．

［8］ Zhao L，Wen L X，Chen L，et al．A two－dimensional hybrid method for modeling seismic wave propagation in anisotropic media．J．Geophys．Res．，2008，113：B12307，doi：10．1029/2008JB005733．

［9］ Komatitsch D，Tromp J．Spectral－element simulations of global seismic wave propagation－I．Validation．Geophys．J．Int．，2002，149（2）：390－412．

［10］ Komatitsch D，Tromp J．Spectral－element simulations of global seismic wave propagation－II．3－D models，oceans，rotation and self－gravitation．Geophys．J．Int．，2002，150（1）：303－318．

［11］ Liu Q，Polet J，Komatitsch D，et al．Spectral－element moment tensor inversions for earthquakes in southern California．Bull．Seismol．Soc．Am．，2004，94（5）：1748－1761．

［12］ Capdeville Y，Chaljub E，Vilotte J P，et al．Coupling the spectral element method with a modal solution for elastic wave propagation in global earth models．Geophys．J．Int．，2003，152（1）：34－66．

［13］ Capdeville Y，To A，Romanowicz B．Coupling spectral elements and modes in a spherical Earth：an extension to the＇sandwich＇case．Geophys．J．Int．，2003，154（1）：44－57．

［14］ Givoli D，Keller J B．Non－reflecting boundary conditions for elastic waves．WaveMotion，1990，12（3）：261－279．

［15］ Chaljub E，Komatitsch D，Vilotte J P，et al．Spectralelement analysis in seismology．∥ Wu R S，Maupin V．Advances in Geophysics，Advances in Wave Propagation in Heterogeneous Earth．San Diego：Elsevier Academic Press Inc．，2007，48：365－419．

［16］ Love A E H．A Treatise on the Mathematical Theory of Elasticity．New York：Dover Publications，1944．

［17］ Komatitsch D，Vilotte J P，Vai R，et al．The spectral element method for elastic wave equations：application to 2D and 3D seismic problems．InternationalJournalfor NumericalMethodsinEngineering，1999，45：1139－1164．

［18］ Ronchi C，Ianoco R，Paolucci P S．The Cubed Sphere：a new method for the solution of partial differential equations in spherical geometry．J．Comput．Phys．，1996，124（1）：93－114．

［19］ Phinney R A，Burridge R．Representation of the elasticgravitational excitation of a spherical earth model by generalized spherical harmonics．GeophysicalJournal International，1973，34（4）：451－487．

［20］ To A，Romanowicz B，Capdeville Y，et al．3D effects of sharp boundaries at the borders of the African and Pacific superplumes：Observation and modeling．EarthPlanet．Sci．Lett．，2005，233（1－2）：137－153．

［21］ Qin Y L，Capdeville Y，Montagner J P，et al．Reliability of mantle tomography models assessed by spectral element simulation．GeophysicalJournalInternational，2009，177（1）：125－144．

［22］ Dahlen F A，Tromp J．Theoretical Global Seismology．Princeton：Princeton University Press，1998．

［23］ Komatitsch D，Vinnik L P，Chevrot S．SHdiff－SVdiff splitting in an isotropic Earth．J．Geophys．Res．，2010，115（B7）：B07312．

［24］ Grand S P．Mantle shear－wave tomography and the fate of subducted slabs．Philosophical Transactions of the Royal Society of London．Series A：Mathematical．Physicaland EngineeringSciences，2002，360（1800）：2475－2491．

［25］ He Y M，Wen L X，Zheng T Y．Geographic boundary and shear wave velocity structure of the"Pacific anomaly"near the core－mantle boundary beneath western Pacific．EarthandPlanetaryScienceLetters，2006，244（1－2）：302－314．

［26］ Panning M，Romanowicz B．A three－dimensional radially anisotropic model of shear velocity in the whole mantle．GeophysicalJournalInternational，2006，167（1）：361－379．

［27］ Yang H Y Y，Zhao L，Hung S H．Synthetic seismograms by normal－mode summation：a new derivation and numerical examples．GeophysicalJournalInternational，2010，183（3）：1613－1632．

［28］ Al－Attar D，Woodhouse J H，Deuss A．Calculation of normal mode spectra in laterally heterogeneous earth models using an iterative direct solution method．GeophysicalJournal International，2012，189（2）：1038－1046．