地殼分層各向異性介質(zhì)接收函數(shù)及其粒子群反演

齊少華， 劉啟元， 陳九輝， 郭飚

中國地震局地質(zhì)研究所地震動力學國家重點實驗室， 北京 100029

?

地殼分層各向異性介質(zhì)接收函數(shù)及其粒子群反演

齊少華， 劉啟元*， 陳九輝， 郭飚

中國地震局地質(zhì)研究所地震動力學國家重點實驗室， 北京 100029

地殼不同深度介質(zhì)的地震各向異性是研究地殼不同深度范圍變形方式的重要依據(jù).鑒于地殼介質(zhì)的復雜性，如何從遠震體波接收函數(shù)中提取不同深度的各向異性參數(shù)仍是一個有待深入研究的課題.在已有研究的基礎(chǔ)上，本文利用廣義反射-透射系數(shù)矩陣方法計算的合成地震圖，研究了復雜地殼分層各向異性介質(zhì)的接收函數(shù)隨反方位角(back azimuth)變化及不同層位各向異性參數(shù)對接收函數(shù)波場的影響，為各向異性介質(zhì)接收函數(shù)的解釋提供了新的理論依據(jù).通過引入粒子群優(yōu)化理論，發(fā)展了分層各向異性介質(zhì)接收函數(shù)全局反演算法.數(shù)值及觀測數(shù)據(jù)的驗證結(jié)果表明，在各向同性速度模型確定的前提下，我們的方法能夠可靠地提取地殼分層各向異性參數(shù)；在反演中引入曲波變換去噪技術(shù)，對于正確解析不同層位的各向異性參數(shù)具有重要價值.

地殼結(jié)構(gòu)； 分層各向異性； 接收函數(shù)； 粒子群反演

1 引言

地震各向異性與地球介質(zhì)的礦物組成、溫壓條件、應力狀態(tài)以及變形歷史密切相關(guān)，它已成為推測地球深部介質(zhì)變形的重要依據(jù)(Babuska and Cara,1991;Rabbel and Mooney, 1996).早期研究認為，地球內(nèi)部的地震各向異性主要表現(xiàn)在上地?；蛘呱系蒯ｍ敳?，它源于地幔物質(zhì)在高溫高壓環(huán)境下的塑性流動所形成的橄欖石晶格優(yōu)勢排列(Nicolas and Christensen,1987;Savage, 1999).遠震SKS分裂方法憑借其簡潔和直觀等優(yōu)點，被廣泛用于上地慢的方位各向異性研究(Vinnik et al., 1989;Silver and Chan,1991).

近年來，日本Chugoku地區(qū)的遠震P波接收函數(shù)的研究發(fā)現(xiàn)，莫霍面Ps轉(zhuǎn)換震相的快波與慢波到時差可達到0.2～0.7 s(Nagaya et al., 2008)，表明地殼各向異性對SKS的影響不可忽略.理論地震圖研究也證實，當?shù)貧て骄飨虍愋詮姸瘸^1%時，透射SKS波就可能發(fā)生分裂(姚陳等,2016).因此，地殼各向異性研究不僅對于理解地殼變形機制等動力學問題至關(guān)重要，而且對于上地幔各向異性的正確解釋也是極其必要的.

地殼各向異性的成因機理在不同深度上可能并不相同.上地殼的地震各向異性可能主要源于裂隙和微裂隙在應力驅(qū)動下的優(yōu)勢排列(Crampin,1987)，而中下地殼各向異性則可能是由于韌性變形或塑性流動導致的礦物優(yōu)勢排列(Weiss et al., 1999;Savage,1999).這意味著地殼各向異性具有分層的屬性，以至于基于單層模型只能得到地殼各向異性較為粗略的估計.接收函數(shù)對地殼的速度間斷面十分敏感，它在研究地殼各向異性及其分層結(jié)構(gòu)方面具有難以替代的優(yōu)勢，現(xiàn)今已成為人們關(guān)注和研究的新熱點(房立華和吳建平, 2009).

早在20世紀80年代末，姚陳(1989)就曾根據(jù)短周期遠震體波記錄上的莫霍面Ps轉(zhuǎn)換波的分裂特征研究了地殼裂隙介質(zhì)的各向異性.McNamara和Owens(1993)則利用接收函數(shù)(本文均指遠震P波接收函數(shù))莫霍面Ps轉(zhuǎn)換波的分裂特征研究了美國盆嶺地區(qū)的地殼各向異性.通過合成理論地震圖的研究，Levin和Park(1997)證明，不僅是傾斜界面，地殼各向異性也可以導致接收函數(shù)切向分量的能量.利用有限差分算法計算的合成地震圖，Okaya和McEvilly(2003)研究了各向異性對稱軸傾斜對地震波場的影響.Savage(1998)、Frederiksen和Bostock(2000)、徐震等(2006)則分別研究了傾斜界面與各向異性的接收函數(shù)波場差異.

在任意取向橫向各向同性(ATI)模型假定下，人們研究了地殼分層各向異性接收函數(shù)隨反方位角的周期性變化特征(Frederiksen and Bostock,2000;田寶峰等,2008).Nagaya等(2008)計算了雙層地殼各向異性情況下，各向異性對稱軸方位角不同的接收函數(shù).根據(jù)觀測數(shù)據(jù)的結(jié)果，齊少華等(2009)證明，水平對稱軸模型可能并不適合解釋變形復雜地區(qū)的地殼各向異性，在地殼范圍內(nèi)考慮ATI模型是極為必要的.但是，某些復雜的各向異性波場特征(如多層各向異性的相互影響、各向異性的沉積層和層間薄層界面的影響)，則尚未給予足夠的關(guān)注.

盡管如此，各向異性接收函數(shù)的理論研究已為利用接收函數(shù)波形反演方法獲取地殼不同層位的地震各向異性參數(shù)奠定了基礎(chǔ).各向異性的接收函數(shù)反演方法已在若干地區(qū)的地殼各向異性研究中取得了有價值的結(jié)果(Ozacar and Zandt,2004;Sherrington et al.,2004;田寶峰等, 2008).相較于剪切波分裂方法，基于波形反演的方法更適用于較為復雜的分層ATI介質(zhì).近年來，回避接收函數(shù)反演的各向異性參數(shù)提取方法也得到了發(fā)展(Liu and Niu,2012;Rümpker et al., 2014;Schulte-Pelkum and Mahan,2014)，并已引起了人們的興趣.與接收函數(shù)反演方法一樣，這些提取地殼各向異性參數(shù)的新方法主要依賴于單個臺站接收函數(shù)(徑向及切向分量)隨反方位角的變化特征.

實際觀測的接收函數(shù)存在介質(zhì)橫向非均勻散射的干擾，特別是接收函數(shù)的切向分量，它對殼內(nèi)橫向散射十分敏感(齊少華等,2016).另外，地表沉積層的混響效應也會導致后續(xù)界面轉(zhuǎn)換震相的信噪比降低.所有這些都增加了上述方法實際應用的難度.需要指出的是，現(xiàn)有的接收函數(shù)波場分析方法有賴于接收函數(shù)隨反方位角變化的完整性.但是，由于臺站(特別是流動臺站)所處的地理位置，采集的接收函數(shù)在方位分布上通常是間斷的，甚至可能存在大范圍的缺失，這給本來就存在的界面傾斜與各向異性的折衷(trade-off)帶來了更多的不確定性(Girardin and Farra,1998;Liu and Niu,2012;Rümpker et al., 2014).

本文的目的在于：(1)針對已有工作的不足，進一步研究較為復雜的地殼分層各向異性，包括層間薄層界面及沉積蓋層的接收函數(shù)；(2)我們將引入曲波變換和壓縮感知(齊少華等,2016)，壓制橫向散射噪聲，重構(gòu)接收函數(shù)波場；(3)在此基礎(chǔ)上，通過數(shù)值和實際數(shù)據(jù)的檢驗，研究接收函數(shù)反演方法提取復雜分層結(jié)構(gòu)不同層位各向異性參數(shù)的能力.這對于進一步理解各向異性接收函數(shù)波場和研究地殼分層各向異性都具有重要意義.

2 分層各向異性介質(zhì)的接收函數(shù)

2.1 地殼各向異性的類型

圖1 傾斜對稱軸各向異性模型(據(jù)Leidig和Zandt(2003)修改)(a) 慢軸各向異性模型； (b) 快軸各向異性模型.淺灰色細線橢圓代表各向同性面.對稱軸與豎直方向的夾角為各向異性對稱軸的傾角.Fig.1 Anisotropic models with a tilted axis of symmetry (modified from Leidig and Zandt, 2003)(a) Slow-axis anisotropy model; (b) Fast-axis anisotropy model. The light-grey thin ellipse represents the isotropic plane. The angle of the symmetry axis with the vertical direction is the plunge of the symmetry axis.

如圖1所示，快軸與慢軸各向異性分別代表兩種不同類型的各向異性模型.根據(jù)Levin和Park(1997)的描述，快軸各向異性的相速度曲面為長球面，形態(tài)類似“西瓜”；而慢軸各向異性其相速度曲面為扁球面，形態(tài)類似“南瓜”.礦物晶體的優(yōu)勢排列導致快軸各向異性，而上地殼裂隙的優(yōu)勢排列和云母的面理組構(gòu)(foliated mica fabric)會導致慢軸各向異性(Leidig and Zandt,2003).巖石學證據(jù)業(yè)已證明，云母面理組構(gòu)可能是中、下地殼各向異性的主要成因(Weiss et al., 1999;Lloyd et al., 2009;Ward et al., 2012).以往的研究通常根據(jù)實際情況，對各向異性的類型做出選擇(Sherrington et al.,2004;田寶峰等, 2008).但是，越來越多的證據(jù)表明，慢軸各向異性模型通?？梢垣@得更小的數(shù)據(jù)擬合殘差，似乎更為適合解釋地殼內(nèi)部的地震各向異性(Leidig and Zandt, 2003; Ozacar and Zandt, 2004; Zandt et al., 2004; Porter et al., 2011; Audet, 2015).因此，本文將統(tǒng)一采用慢軸各向異性模型來描述地殼各個層位的地震各向異性.

2.2 地殼分層各向異性接收函數(shù)的計算方法

為了研究地殼分層各向異性，本文將依據(jù)廣義反射透射系數(shù)矩陣方法(Chen,1993)進行理論接收函數(shù)的計算，并采用作者自行開發(fā)的程序代碼.這有助于考慮地殼復雜分層模型的完整波場.本文算法與Levin和Park(1997)沒有實質(zhì)性區(qū)別，但在以下方面采取了不同的處理：在構(gòu)建各向異性彈性參數(shù)時，我們采用了姚陳和蔡明剛(2009)推導的ATI彈性張量解析公式，在一定程度上提高了計算的精度.

根據(jù)弱各向異性假定(Backus,1965;Park and Yu,1992)，各向異性彈性參數(shù)可以近似地分為各向同性的背景值(由P波速度、S波速度和密度來描述)和擾動值(由各向異性對稱軸的方位角、傾角、P波和S波各向異性強度表征)，這種參數(shù)化方法便于正演模型的構(gòu)建，減少反演參數(shù)(Levin and Park,1997).

2.3 快軸與慢軸各向異性接收函數(shù)的比較

圖2給出了采用本文方法計算的單層地殼模型的快軸與慢軸各向異性的理論接收函數(shù).所用的水平慢度為0.0618 s·km-1，這大體相當于震中距為60°的情況.表1給出了相應的地殼模型參數(shù).其中，模型M1和M2分別為快軸各向異性(強度為正)和慢軸各向異性(強度為負).除了快軸與慢軸的方位角相差180°以外，其他參數(shù)沒有區(qū)別.

比較圖2a和2b可知，模型底部界面的透射Ps轉(zhuǎn)換震相及多次波的振幅及相位隨反方位角變化的周期性是一樣的：均以0°和180°反方位角為軸，徑向分量呈現(xiàn)對稱，切向分量呈現(xiàn)反對稱.但是，它們的振幅及相位隨反方位角變化的具體特征是不同的，特別是，切向分量的多次波呈現(xiàn)了相反的極性.這些差異或許可以作為我們判別各向異性類型的依據(jù).需要說明的是，快軸與慢軸方位角相差180°與傾角互補是等價的(Levin and Park,1997;Erickson, 2002).

表1 快軸與慢軸各向異性模型

注：H為層厚，vP和vS分別為P波和S波速度，ρ為密度，B和E分別為P波和S波各向異性強度(快軸為正，慢軸為負).θ和φ分別為各向異性對稱軸的傾角和方位角.所有符號下文所有各表中均相同.

2.4 地殼分層各向異性接收函數(shù)

接收函數(shù)是觀測臺站接收區(qū)底部入射波的脈沖響應.這意味著，對于分層各向異性的地殼，各界面的透射與反射波均會涉及所有層位的各向異性.因此，對于較為復雜的分層結(jié)構(gòu)，能否從接收函數(shù)正確地解析各個層位的各向異性參數(shù)無疑是一個值得關(guān)注的問題.另外，對于構(gòu)造活動區(qū)，地殼通常都存在不同程度的運動變形，并伴隨著殼內(nèi)的低速層(如青藏高原及我國的華北地區(qū)).在應力作用和應變調(diào)節(jié)下，它們有可能成為地殼物質(zhì)流動的通道(劉啟元等,2007;Beaumont et al.,2001;Liu et al., 2014)，意味著低速層可能具有較強的各向異性.劉啟元和邵學鐘(1985)曾證明，薄層的混響效應會導致相應的轉(zhuǎn)換震相能量的增強.如果它同時顯示了明顯的各向異性特征，那么是否會誤導分層各向異性參數(shù)的估計，也是一個有待研究的問題.因此，更為復雜的地殼分層各向異性將是本文關(guān)注的重點.我們將通過計算分層各向異性模型的接收函數(shù)隨反方位角的變化，來考察不同層位各向異性參數(shù)對接收函數(shù)波場的影響.

我們首先考慮雙層各向異性的情況，表2給出了相應的地殼模型.首先，我們考慮上下兩層各向異性參數(shù)相同的情況(模型M3).圖3a和3b分別給出了計算得到的接收函數(shù)徑向和切向分量隨反方位角的變化.比較圖2c—2d(單層地殼)和圖3a—3b(雙層地殼)，可有如下觀察：(1)由于層間界面轉(zhuǎn)換波(Ps1)及多次波(PpPs1和PsPs1+PpSs1)的出現(xiàn)，并不難做出雙層地殼結(jié)構(gòu)的判斷；(2)切向分量的初至P波投影的極性相同；(3)層間界面Ps1震相的徑向分量與其切向分量均顯示了近似各向同性界面的特征；(4)在徑向和切向分量上，層間界面Ps1震相的多次波的各向異性特征明顯；(5)兩者底部界面的轉(zhuǎn)換波和多次波徑向及切向分量沒有明顯差別.上述觀察表明，當雙層各向異性參數(shù)相同且不考慮多次波時，其與單層模型的各向異性特征基本相同.

圖3 雙層各向異性模型的接收函數(shù)波場：上下層各向異性參數(shù)相同(a) 徑向分量； (b) 切向分量.Fig.3 RF-wavefield from two-layered anisotropic model: Same anisotropic parameters(a) R-component; (b) T-component.

H(km)vP(km·s-1)vS(km·s-1)ρ(g·cm-3)BEθ(°)?(°)M3206.5003.8502.7142-5%-5%450157.2004.1143.0876-5%-5%450∞8.0404.4803.2976----M4206.5003.8502.7142-5%-5%450157.2004.1143.0876-5%-5%750∞8.0404.4803.2976----M5206.5003.8502.7142-5%-5%450157.2004.1143.0876-5%-5%4590∞8.0404.4803.2976----M6206.5003.8502.7142-5%-5%450157.2004.1143.0876-5%-5%7590∞8.0404.4803.2976----M7206.5003.8502.7142-5%-5%750157.2004.1143.0876-5%-5%450∞8.0404.4803.2976----M8206.5003.8502.7142-5%-5%4590157.2004.1143.0876-5%-5%450∞8.0404.4803.2976----M9206.5003.8502.7142-5%-5%7590157.2004.1143.0876-5%-5%450∞8.0404.4803.2976----

注：各向異性對稱軸的方位角按順時針方向旋轉(zhuǎn)并以正北方向為0°；傾角按順時針方向旋轉(zhuǎn)，豎直向上為0°.下文各表相同.

下面，我們考慮各向異性參數(shù)不同的雙層地殼模型.我們先固定上層，僅改變下層的各向異性參數(shù)(模型M4—M6).圖4給出了相應的接收函數(shù)隨反方位角的變化.

圖4a和4b分別給出了僅對稱軸傾角不同時(M4)的接收函數(shù)徑向與切向分量.與圖3a—3b相比，下層各向異性對稱軸傾角的改變使接收函數(shù)徑向及切向分量的各向異性特征(特別是層間界面Ps1震相)在保持周期性不變的情況下得到增強，意味著層間界面震相各向異性特征的強弱主要取決于各向異性參數(shù)的層間差異.

圖4c和4d分別給出了僅對稱軸方位角正交(M5)的接收函數(shù)徑向與切向分量.與圖3a—3b和圖4a—4b對比，不難發(fā)現(xiàn)：(1)在切向分量上，零時刻P波投影沒有因為下層各向異性參數(shù)的改變而改變(如圖4d兩個橙色箭頭所示)；(2)無論改變下層各向異性對稱軸的傾角，還是方位角，層間界面的多次波(PpPs1和PsPs1+PpSs1)的徑向分量均未有視覺上的變化；(3)下層各向異性對稱軸方位角的改變所引起的底部界面轉(zhuǎn)換波(Ps2)的變化明顯大于傾角改變所引起的變化，因為對稱軸方位角的改變導致其對稱位置平移了90°(如圖4d中的綠色箭頭所示)；(4)底部界面的多次波(PpPs2和PsPs2+PpSs2)的徑向及切向分量的各向異性特征異常削弱.上述觀察表明，接收函數(shù)波場對各向異性對稱軸的方位角更為敏感，它可以引起接收函數(shù)隨反方位角變化的對稱位置的改變.特別是，當上下兩層的各向異性對稱軸的方位角正交時，下層界面轉(zhuǎn)換波的各向異性特征會被嚴重削弱.

圖4 雙層各向異性模型的接收函數(shù)波場：上層參數(shù)固定，下層參數(shù)改變(a) M4，徑向； (b) M4，切向； (c) M5，徑向； (d) M5，切向； (e) M6，徑向； (f) M6，切向.Fig.4 RF wavefield from two-layered anisotropic model: Lower parameters modified(a) M4,R-component; (b) M4,T-component; (c) M5,R-component.(d) M5,T-component; (e) M6,R-component; (f) M6, T-component.

圖4e和4f分別給出了對稱軸方位角正交且傾角不同時(M6)的接收函數(shù)徑向與切向分量.與模型M5的情況相比，其波場特征未有明顯變化，正如模型M3與M4的接收函數(shù)波場也沒有明顯差異的情況類似.這意味著接收函數(shù)波場的各向異性特征對傾角變化確實并不敏感.

表2中模型M7—M9給出了固定下層各向異性參數(shù)，僅改變上層參數(shù)時的雙層地殼模型，其接收函數(shù)波場如圖5所示.對比圖4和圖5，不難發(fā)現(xiàn)：對于雙層地殼各向異性，(1)僅改變下層的各向異性參數(shù)與僅改變上層各向異性參數(shù)的接收函數(shù)波場具有十分類似的現(xiàn)象，兩者均表明，接收函數(shù)各向異性特征的變化對介質(zhì)各向異性對稱軸的方位角比對其傾角更為敏感；(2)盡管圖4和圖5的波場特征明顯不同，但兩者呈現(xiàn)了某種程度的“互補性”，即與改變下層各向異性參數(shù)相反，上層各向異性參數(shù)的改變導致該層位相應震相的各向異性特征更為明顯地改變.例如，與圖4不同，圖5橙色箭頭指示了上層各向異性軸方位的改變，導致初至P波投影隨反方位角變化的對稱位置移動了90°，而圖5中綠色箭頭指示的對稱位置與圖3并沒有什么不同，盡管下層各向異性參數(shù)將同時影響其上下兩個界面的透射轉(zhuǎn)換震相(圖4中的Ps1和Ps2)；(3)無論改變上層，還是下層的各向異性參數(shù)，接收函數(shù)的切向分量較其徑向分量對各向異性參數(shù)的改變更為敏感.

圖5 雙層各向異性模型的接收函數(shù)波場：下層參數(shù)固定，上層參數(shù)改變(a) M7，徑向； (b) M7，切向； (c) M8，徑向； (d) M8，切向； (e) M9，徑向； (f) M9，切向.Fig.5 RF wavefield from two-layered anisotropic model: upper parameters modified(a) M7,R-component; (b) M7,T-component; (c) M8,R-component; (d) M8,T-component; (e) M9,R-component; (f) M9,T-component.

表3給出了包含低速薄層、沉積層以及莫霍過渡層等更為復雜的分層各向異性模型.模型M10為各向異性低速薄層的上下均為各向同性介質(zhì)的地殼模型，相應接收函數(shù)的徑向和切向分量分別由圖6a和6b給出.觀察可知，薄層界面的轉(zhuǎn)換震相Ps1在徑向和切向分量上均顯示了清晰的各向異性特征，但該薄層界面上的反射波(PpPs1和PsPs1+PpSs1)的各向異性效應非常微弱，薄層對來自底部界面的震相Ps2并未形成可以察覺的干擾.這意味著，對于各向同性的雙層地殼來說，各向異性薄層界面對接收函數(shù)波場的影響十分有限.值得關(guān)注的是，與各向同性薄層不同(劉啟元和邵學鐘,1985)，在徑向分量反方位角180°附近，各向異性導致了層間薄層界面轉(zhuǎn)換震相能量的離散.

模型M11是一個各向同性低速薄層與其上下層均為各向異性層的地殼模型，圖6c和6d分別給出了相應的接收函數(shù)徑向和切向分量.模型M12與模型M11的不同是其薄層界面為各向異性低速薄層，圖6e和6f給出了M12的接收函數(shù)徑向及切向分量.對比模型M11和M12的接收函數(shù)可知，兩者的差別不大，但薄層界面的轉(zhuǎn)換波隨反方位角的變化特征是有差別的，這意味著接收函數(shù)波場的反方位完整性對于判別薄層是否存在各向異性十分重要.這也提示我們，根據(jù)有限反方位的接收函數(shù)波場只能對地殼分層各向異性做出粗略的估計.

表3 多層各向異性模型

模型M13是一個底部為過渡層(厚度8 km)，其上層為各向同性層的地殼模型.圖7a和7b分別給出了相應的接收函數(shù)徑向及切向分量.由圖7a和7b可見，由于各向異性的存在，該過渡層上的轉(zhuǎn)換震相(黑色箭頭所示的震相)在反方位角180°附近出現(xiàn)了離散，而由于上行(下行)的次數(shù)增加，其相應多次波在過渡層頂部與底部徹底分開，且沒有顯現(xiàn)明顯的各向異性效應.綜合低速薄層(圖6a—6b)和過渡層(圖7a—7b)的接收函數(shù)波場可知，其各向異性對波場的影響局限于透射轉(zhuǎn)換波，而多次波對其各向異性并不敏感，這與厚層介質(zhì)的情況完全不同.圖7a—7b和圖2c—2d差異提供的證據(jù)表明，厚層和薄層各向異性對接收函數(shù)波場的影響是相互獨立的.但是，在提取各向異性參數(shù)時，考慮多次反射與采用完整反方位的接收函數(shù)對于避免兩者混淆均具有重要作用.

模型M14是在模型M10的基礎(chǔ)上增加了地表2 km沉積薄層后的地殼模型.圖7c和7d分別給出了相應接收函數(shù)的徑向和切向分量.圖7c和圖7d表明，沉積薄層對接收函數(shù)的徑向和切向分量均產(chǎn)生了極為顯著的影響.盡管如此，對比圖6a—6b可知，沉積薄層與低速薄層的各向異性特征(包括其對稱性和周期性)依然可在切向分量上加以區(qū)分.

3 各向異性介質(zhì)接收函數(shù)的粒子群反演

本節(jié)將進一步研究從接收函數(shù)提取地殼分層各向異性參數(shù)的能力.如前所述，接收函數(shù)反演是提取地殼分層各向異性參數(shù)的常用方法.前人廣泛采用了遺傳、近鄰和模擬退火等全局性算法，以避免陷于局部極小(Ozacar and Zandt,2004;Sherrington et al., 2004;田寶峰等,2008).與各向同性介質(zhì)的接收函數(shù)反演不同，地殼各向異性參數(shù)的提取需要對不同反方位角的接收函數(shù)徑向與切向波形進行聯(lián)合反演.因而，反演算法的效率，特別當針對大量臺站的研究時，成為一個不可回避的問題.

圖6 包含低速薄層的各向異性模型的接收函數(shù)波場(a) M10，徑向； (b) M10，切向； (c) M11，徑向； (d) M11，切向； (e) M12，徑向； (f) M12，切向. Fig.6 RF wavefield from anisotropic model with low-velocity thin layer(a) M10,R-component; (b) M10,T-component; (c) M11,R-component; (d) M11,T-component; (e) M12,R-component; (f) M12,T-component.

圖7 包含沉積層和過渡層的各向異性模型的接收函數(shù)波場(a) M13，徑向； (b) M13，切向； (c)M14，徑向； (d) M14，切向.Fig.7 RF wavefield from anisotropic model with sedimentary and transition layer(a) M13,R-component; (b) M13,T-component; (c) M14,R-component; (d) M14,T-component.

圖8 各向異性接收函數(shù)反演的數(shù)值檢驗(a) “真實”接收函數(shù)； (b) 預測接收函數(shù)； (c) 各向異性對稱軸的方位角； (d) 各向異性對稱軸的傾角； (e) 各向同性速度結(jié)構(gòu)與各向異性強度； (f) “真實”與預測接收函數(shù)之間的殘差. sum表示疊加波形； Ani表示各向異性強度； 彩色圓圈代表每個模型的殘差； 紅色圓圈代表全局殘差，陰影部分為各向同性層.Fig.8 Numerical test on the inversion of seismic anisotropy from RF(a) ‘true’ RF; (b) Predicted RF; (c) Trend of anisotropy; (d) plunge of anisotropy; (e) Isotropic velocity model and magnitude of anisotropy; (f) Misfit between ‘true’ and predicted RF. ‘sum’ denotes the summation over all backazimuths;‘Ani’ denotes magnitude of anisotropy; colored solid circle denotes the misfit for each model;red solid circle denotes the global misfit; shaded area represents isotropic layer.

近年來，一種稱為“粒子群優(yōu)化”(particle swarm optimization，簡稱PSO)的全局算法得到了發(fā)展和完善(Eberhart and Kennedy,1995;Shi and Eberhart,1998;Clerc and Kennedy,2002).該算法基于鳥群和魚群覓食過程中表現(xiàn)出的群體智能(swarm intelligence)，具有實現(xiàn)簡單、易于并行、收斂速度快以及需要調(diào)整參數(shù)少等優(yōu)勢(Eberhart and Shi,1998;Hassan et al., 2005).該方法在國內(nèi)外正逐步得到日益廣泛的關(guān)注和應用(Shaw and Srivastava,2007;Martínez et al., 2010;師學明等,2009;岳碧波等,2009;朱童等,2011).我們將把該算法引入到地殼分層各向異性的接收函數(shù)反演.

在反演算法中，正演計算采用前文(2.2節(jié))所述的方法，粒子的位置即為模型空間里的一個樣點(Eberhart and Kennedy, 1995)，反演的目標函數(shù)則基于時間域波形互相關(guān)(Frederiksen et al., 2003).對于PSO算法，單個粒子依據(jù)其自身和其他粒子的飛行經(jīng)驗，調(diào)整自身的飛行速度，并以迭代方式向全局最優(yōu)解的位置逼近.在迭代過程中，粒子既不會消亡，也不會被丟棄.粒子自身和其他粒子的飛行經(jīng)驗在算法執(zhí)行的過程中被完整地保留，并指導著粒子的后續(xù)飛行.正是這種個體與自身、個體與群體之間信息分享的拓撲結(jié)構(gòu)，使得PSO算法同時兼?zhèn)淙峙c局部搜索的特性.

為了權(quán)衡PSO全局與局部搜索的能力，本文采用Shi和Eberhart(1998)發(fā)展的粒子速度更新方程.同時，為避免迭代中后期出現(xiàn)所謂的“粒子爆炸”現(xiàn)象(即粒子由于速度過大而逃離搜索空間)，引入“收縮因子”來避免“粒子爆炸”，同時使算法收斂速度進一步加快(Clerc and Kennedy,2002).

理論上，PSO的收斂判據(jù)應該是粒子群中所有粒子均重疊在模型空間的同一個點上.然而，由于粒子存在慣性，且不同粒子的運行軌跡有很大的差異，所有粒子最終靜止在同一點上可能需要很長的時間.為此，人們把所有粒子均進入一個合理的誤差范圍作為反演終止的判據(jù)，或者直接規(guī)定迭代的最大次數(shù).

表4給出了數(shù)值檢驗所使用的“真實”接收區(qū)模型.計算所用的水平慢度設(shè)定為0.0618 s·km-1.由表4可知，我們在模型中考慮了低速層、高速層以及薄低速蓋層.需要說明的是，在數(shù)值檢驗中，(1)我們假定各向同性的速度模型已知，且各向異性參數(shù)的分層與速度參數(shù)保持一致，即我們僅把各向異性參數(shù)作為獨立變量；(2)初始各向異性參數(shù)均設(shè)為零；(3)對反演參數(shù)作如下約束：P波與S波各向異性強度保持相同，取值范圍限定在[-20%,0]，各向異性對稱軸的傾角限定在[0°,90°]，方位角限定在[0°,360°]；(4)粒子數(shù)量設(shè)定為50，最大迭代次數(shù)為100次.

圖8給出了數(shù)值檢驗結(jié)果.比較圖8a和8b可知，反演得到的接收函數(shù)波場與“真實”數(shù)據(jù)沒有視覺上的差別.圖8f中紅色圓圈描述的收斂曲線表明，大約經(jīng)過50次左右的迭代后，目標函數(shù)的殘差即已基本收斂.

比較表4和圖8(c,d, e)可知，反演結(jié)果與“真實”的各向異性模型基本一致.這表明在各向同性速度模型已知，且不含噪聲的條件下，利用本文方法可以從接收函數(shù)波場可靠地提取各個層位的各向異性參數(shù).

與上述數(shù)值檢驗不同，實際觀測數(shù)據(jù)難以避免橫向非均勻構(gòu)造的散射及各種噪聲.這往往成為從觀測接收函數(shù)中提取地震各向異性參數(shù)的難點(Savage,1998).與理論的各向異性介質(zhì)接收函數(shù)波場所展示的隨反方位角的可追蹤性不同，橫向非均勻散射及各種噪聲在空間取向上可能具有隨機性.齊少華等(2016)給出的結(jié)果表明，利用曲波變換對接收函數(shù)波場進行二維去噪不僅可以明顯地改善震相的可追蹤性，而且可同時通過壓縮感知波場重建技術(shù)獲得全方位的接收函數(shù)，無疑有助于分層各向異性介質(zhì)接收函數(shù)反演研究.

表4 合成數(shù)據(jù)檢驗所用的各向異性模型

齊少華等(2016)曾利用曲波去噪方法處理過位于澳大利亞東北地區(qū)的全球數(shù)字地震臺網(wǎng)(IRIS)臺站CTAO(-20.0882°N，146.2545°E，臺網(wǎng)編號IU)的接收函數(shù).該臺站所處的地理位置和構(gòu)造背景如圖9所示.在此基礎(chǔ)上，我們將進一步利用本文的分層各向異性介質(zhì)接收函數(shù)反演技術(shù)，提取該臺下方的地殼各向異性參數(shù)，其目的在于檢驗各向異性接收函數(shù)反演解釋實際觀測數(shù)據(jù)的能力.

圖9 臺站CTAO的地理位置及其構(gòu)造背景(修改自Ford等,2010)黑色實線勾畫出了澳大利亞的主要構(gòu)造區(qū)域； 西部斜線區(qū)域代表太古宙克拉通； 北部和南部橫線區(qū)域代表元古宙克拉通； 東部點陣區(qū)域顯示了顯生宙基底； 黑色三角表示IRIS臺網(wǎng)的臺站， 紅色三角表示本文使用的臺站.Fig.9 Geographic location of Station CTAO and its tectonic background (modified from Ford et al.,2010).Black solid lines outline the major tectonic provinces in Australia;western regions with inclined lines indicate Archean cratons; northern and southern regions with horizontal lines indicate Proterozoic cratons; eastern regions with dots indicate Phanerozoic basements; black triangles denote IRIS stations,while red triangle represents station CTAO in this study.

齊少華等(2016)已經(jīng)描述了利用時間域迭代反褶積方法(Ligorría and Ammon,1999)得到的151個遠震事件(1999年2月—2015年7月，震中距30°≤Δ≤90°，Mb>5.5)在CTAO臺的接收函數(shù).為此，本文不再展示未經(jīng)曲波去噪處理的接收函數(shù)，僅簡單綜述基本的數(shù)據(jù)參數(shù)如下： (1) 利用時間域迭代反褶積方法提取接收函數(shù)時，采用的高斯系數(shù)為2.5，這大體上相當于0～1.2 Hz的帶通濾波； (2) 波場重建的反方位角間隔為1°.圖10a給出了用于本文反演研究的接收函數(shù)，它是以10°反方位角間隔進行進一步分組疊加后的結(jié)果.

本文的各向異性介質(zhì)接收函數(shù)反演分為以下步驟： (1) 首先對全方位疊加波形進行各向同性接收函數(shù)反演，以便確定臺站下方地殼上地幔速度結(jié)構(gòu)； (2) 對獲得的地殼速度結(jié)構(gòu)進行必要的簡化，形成各向異性參數(shù)的分層結(jié)構(gòu)，且在各向異性反演中固定不變.需要說明的是，速度模型中莫霍界面以下的部分不參與各向異性反演.

圖10給出了利用本文方法得到的反演結(jié)果.表5給出了臺站CTAO下方地殼速度結(jié)構(gòu)及反演得到的分層各向異性參數(shù).由表5可知，該臺下方的地殼厚度為36 km.Chevrot和 van der Hilst(2000)利用接收函數(shù)H-κ掃描得到的地殼厚度為34 km，與我們的結(jié)果基本一致.

由圖10b可知，與觀測的接收函數(shù)波場相比，前5 s窗口內(nèi)，反演的接收函數(shù)徑向與切向分量隨反方位角的變化特征與觀測數(shù)據(jù)吻合較好.需要說明的是，我們所用的目標函數(shù)以波形互相關(guān)為基礎(chǔ)，這意味著主要強調(diào)相位擬合.觀測與預測接收函數(shù)波場存在一定差異，一方面是由于水平分層地殼模型無法解釋三維非均勻介質(zhì)的波場；另一方面，則可能因為在反演中我們沒有考慮莫霍界面以下上地幔的各向異性.臺站下方各層的各向異性強度如圖10e所示.結(jié)果表明，臺站CTAO下方地殼淺部各向異性強度較弱，而深部各向異性強度較強，其強度的大小與前人在其他地區(qū)得到的結(jié)果具有較好的可比性(Sherrington et al., 2004;田寶峰等, 2008)，我們認為這是合理的.

圖10c和10d分別給出了臺站CTAO下方地殼各層位各向異性慢軸的方位角及傾角，具體數(shù)值列于表5.由于沒有關(guān)于臺站CTAO下方地殼各層位各向異性參數(shù)的研究報道，我們無法提供能夠證明本文結(jié)果的直接可比證據(jù).盡管如此，我們可以根據(jù)相關(guān)的間接證據(jù)來說明本文給出的該臺下方各層位慢軸方位角及傾角的合理性.為此，我們給出如下討論：

(1) 臺站CTAO地殼慢軸各向異性傾角由淺到深逐漸趨于垂直方向，意味著中下地殼的各向異性組構(gòu)(anisotropic fabric)可能趨于水平，且水平方向的波速快于垂直方向，這與Debayle和Kennett(2000)得到的巖石圈200 km以上SH波速明顯快于SV的面波徑向各向異性結(jié)果相一致.

(2) 臺站CTAO的下地殼(第4和5層)慢軸各向異性方位角為近EW方向，意味著與其正交的快波方向為近NS方向，與根據(jù)瑞利波方位各向異性所得到的快波方向(80～100 km以上)基本一致(Debayle and Kennett,2000;Simons and van der Hilst,2003).

圖10 臺站CTAO的地殼各向異性反演(a) 觀測接收函數(shù)； (b) 反演得到的接收函數(shù)； (c) 分層各向異性對稱軸方位角； (d) 分層各向異性對稱軸傾角； (e) 各向同性速度結(jié)構(gòu)與各向異性強度； (f) 觀測與預測接收函數(shù)之間的殘差, (e,f)的其他說明同圖8.Fig.10 Inversion of anisotropy from RF at Station CTAO(a) Observed RF; (b) Predicted RF; (c) Trend of anisotropy; (d) Plunge of anisotropy; (e) Isotropic velocity model and magnitude of anisotropy; (f) Misfit between observed and predicted RF. The Others are same as Fig.8.

LayerNo.H(km)vS(km·s-1)θ(°)?(°)14.03.45756.000.0124.03.52015.6926.4938.03.40210.3743.07412.03.5367.8495.3058.04.0613.0583.90

4 結(jié)論與討論

本文利用廣義反射透射系數(shù)矩陣方法計算的合成地震圖，研究了復雜地殼分層各向異性接收函數(shù)隨反方位角的變化和不同層位各向異性參數(shù)對接收函數(shù)波場的影響，為各向異性介質(zhì)接收函數(shù)的解釋提供了新的理論依據(jù)；通過引入粒子群優(yōu)化理論，發(fā)展了地殼分層各向異性接收函數(shù)反演方法，并在此基礎(chǔ)上，通過數(shù)值和觀測數(shù)據(jù)的檢驗，論證了分層各向異性接收函數(shù)反演提取地殼不同層位各向異性參數(shù)的能力.根據(jù)我們的結(jié)果，可以得到以下結(jié)論：

(1) 分層各向異性地殼接收函數(shù)的理論波場研究表明，快軸與慢軸模型的接收函數(shù)隨反方位角變化的特征是不同的；接收函數(shù)波場的各向異性特征與各向異性對稱軸方位角的關(guān)系較其傾角更為密切；接收函數(shù)切向分量較其徑向分量對各向異性參數(shù)更為敏感；層間界面Ps震相隨反方位角的變化特征是界面上下兩層共同作用的結(jié)果，層間界面Ps震相各向異性特征的強弱主要取決于各向異性參數(shù)的層間差異；沉積蓋層的混響效應可以對接收函數(shù)(徑向和切向分量)造成干擾；各向異性薄層界面的轉(zhuǎn)換震相可以導致震相的離散，但其多次波的各向異性效應微弱.本文的結(jié)果表明，偏振分析方法對于分層各向異性研究存在一定局限性.

(2) 數(shù)值檢驗的結(jié)果表明，在地殼各向同性速度結(jié)構(gòu)確定的前提下，分層各向異性介質(zhì)的粒子群接收函數(shù)反演可以正確提取各層位的各向異性參數(shù)；把多次反射納入反演有利于正確解析不同層位的各向異性參數(shù)；由于接收函數(shù)對地殼介質(zhì)的絕對速度并不敏感，在各向異性接收函數(shù)反演前，獨立確定接收區(qū)各向同性速度結(jié)構(gòu)，有助于分層各向異性接收函數(shù)反演的穩(wěn)定性.

(3) 觀測數(shù)據(jù)的驗證結(jié)果表明，曲波變換去噪降低了橫向非均勻散射和其他噪聲的干擾，有效地完善了接收函數(shù)方位的完整性，這對正確解析不同層位的各向異性參數(shù)和避免接收函數(shù)反演中不同層位各向異性參數(shù)的多解性都具有重要價值.

(4) 雖然曲波變換去噪有效降低了橫向非均勻散射和其他噪聲的干擾，但并不意味著完全去除了三維地殼結(jié)構(gòu)對接收函數(shù)波場的影響.理論上，把多次反射納入接收函數(shù)反演有利于正確解析不同層位各向異性參數(shù)，但由于其路徑放大效應，對于基于水平分層介質(zhì)的接收函數(shù)反演來說，多次反射的各向異性特征解析仍然十分困難.

致謝 作者感謝評閱人對本文提出的寶貴建議，這使得本文得到進一步完善.本文所用數(shù)據(jù)來源于IRIS全球地震臺網(wǎng).感謝Brecht Donckels博士提供粒子群優(yōu)化的核心算法.

Audet P. 2015.Layered crustal anisotropy around the San Andreas Fault near Parkfield, California.J.Geophys.Res.,120(5): 3527-3543.

Babuska V, Cara M. 1991. Seismic Anisotropy in the Earth.Netherlands: Kluwer Academic Publishers.

Backus G E. 1965. Possible forms of seismic anisotropy of the uppermost mantle under oceans.J.Geophys.Res.,70(14): 3429-3439.

Beaumont C, Jamieson R A, Nguyen M H, et al. 2001.Himalayantectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation.Nature, 414(6865): 738-742.

Chen X F.1993. A systematic and efficient method of computing normal modes for multilayered half-space.Geophys.J.Int.,115(2): 391-409.

Chevrot S, van der Hilst R D. 2000. The Poisson ratio of the Australian crust: geological and geophysical implications.EarthPlanet.Sci.Lett.,183(1-2): 121-132.

Clerc M, Kennedy J. 2002. The particle swarm-explosion, stability, and convergence in a multidimensional complex space.IEEETransactionsonEvolutionaryComputation, 6(1): 58-73.

Crampin S. 1987. Geological and industrial implications of extensive-dilatancy anisotropy.Nature, 328(6130): 491-496.

Debayle E, Kennett B L N. 2000. Anisotropy in the Australasian upper mantle from Love and Rayleigh waveform inversion.EarthPlanet.Sci.Lett.,184(1): 339-351.

Eberhart R C, Kennedy J. 1995. A new optimizer using particle swarm theory. ∥ Proceedings of the Sixth International Symposium on Micro Machine and Human Science.Nagoya: IEEE, 39-43.

Eberhart R C, Shi Y H. 1998. Comparison between genetic algorithms and particle swarm optimization.∥ Proceedings of the 7th International Conference on Evolutionary Programming. San Diego, California, USA:Springer, 611-616. Erickson J. 2002. Anisotropic crustal structure inversion using a niching genetic algorithm: A feasibility study[Master thesis].Tucson:Universityof Arizona.

Fang L H, Wu J P. 2009. Effects of dipping boundaries and anisotropic media on receiver functions.ProgressinGeophys.(in Chinese), 24(1): 42-50.

Ford H A, Fischer K M, Abt D L, et al. 2010. The lithosphere-asthenosphere boundary and cratonic lithospheric layering beneath Australia from Sp wave imaging.EarthPlanet.Sci.Lett.,300(3-4): 299-310.

Frederiksen A W, Bostock M G. 2000. Modelling teleseismic waves in dipping anisotropic structures.Geophys.J.Int.,141(2): 401-412.

Frederiksen A W, Folsom H, Zandt G. 2003. Neighbourhood inversion of teleseismic Ps conversions for anisotropy and layer dip.Geophys.J.Int.,155(1): 200-212.

Girardin N, Farra V. 1998. Azimuthal anisotropy in the upper mantle from observations of P-to-S converted phases: application to southeast Australia.Geophys.J.Int.,133(3): 615-629. Hassan R, Cohanim B, De Weck O, et al. 2005.A comparison of particle swarm optimization and the genetic algorithm.∥ Proceedings of the 46thAIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference. Austin, Texas: AIAA. Leidig M, Zandt G. 2003. Modeling of highly anisotropic crust and application to the Altiplano-Puna volcanic complex of the central Andes.J.Geophys.Res.,108(B1):ESE 5-1-ESE 5-15, doi: 10.1029/2001JB000649.

Levin V, Park J. 1997. P-SH conversions in a flat-layered medium with anisotropy of arbitrary orientation.Geophys.J.Int.,131(2): 253-266.

Ligorría J P, Ammon C J. 1999. Iterative deconvolution and receiver-function estimation.Bull.Seism.Soc.Am.,89(5): 1395-1400. Liu H F, Niu F L. 2012. Estimating crustal seismic anisotropy with a joint analysis of radial and transverse receiver function data.Geophys.J.Int.,188(1):144-164.

Liu Q Y, Shao X Z. 1985. Study on the dynamic characteristics of PS converted waves.ActaGeophysicaSinica(in Chinese),28(3): 291-302.

Liu Q Y, Wang J, Chen J H, et al. 2007.Seismogenic tectonic environment of 1976 great Tangshan earthquake: results given by dense seismic array observations.EarthScienceFrontiers(in Chinese), 14(6): 205-213.

Liu Q Y, van der Hilst R D, Li Y, et al. 2014. Eastward expansion of the Tibetan Plateau by crustal flow and strain partitioning across faults.NatureGeoscience, 7(5): 361-365.

Lloyd G E, Butler R W H, Casey M, et al. 2009. Mica, deformation fabrics and the seismic properties of the continental crust.EarthPlanet.Sci.Lett.,288(1-2): 320-328.

Martínez J L F, Gonzalo E G, álvarez J P F, et al. 2010. PSO: a powerful algorithm to solve geophysical inverse problems: application to a 1D-DC resistivity case.JournalofAppliedGeophysics, 71(1): 13-25.

McNamara D E, Owens T J. 1993. Azimuthal shear wave velocity anisotropy in the Basin and Range province using Moho Ps converted phases.J.Geophys.Res.,98(B7): 12003-12017.

Nagaya M, Oda H, Akazawa H, et al. 2008. Receiver functions of seismic waves in layered anisotropic media: Application to the estimate of seismic anisotropy.Bull.Seism.Soc.Am.,98(6): 2990-3006.

Nicolas A, Christensen N I. 1987.Formation of anisotropy in upper mantle peridotites—A review.∥ Froidevaux C, Fuchs K, eds.Composition, Structure and Dynamics of the Lithosphere-Asthenosphere System.Washington DC:AGU, 111-123.

Okaya D A, McEvilly T V. 2003. Elastic wave propagation in anisotropic crustal material possessing arbitrary internal tilt.Geophy.J.Int.,153(2): 344-358.

Ozacar A A, Zandt G. 2004. Crustal seismic anisotropy in central Tibet: Implications for deformational style and flow in the crust.Geophys.Res.Lett.,31(23): L23601. Park J, Yu Y. 1992. Anisotropy and coupled free oscillations: simplified models and surface wave observations.Geophys.J.Int.,110(3): 401-420. Porter R, Zandt G, McQuarrie N. 2011. Pervasive lower-crustal seismic anisotropy in Southern California: Evidence for under platedschists and active tectonics.Lithosphere, 3(3): 201-220. Qi S H, Liu Q Y, Chen J H, et al. 2009. Wenchuan earthquakeMs8.0: Preliminary study of crustal anisotropy on both sides of the Longmenshan faults.SeismologyandGeology(in Chinese), 31(3): 377-388. Qi S H, Liu Q Y, Chen J H, et al. 2016. Attenuation of noise in receiver functions using curvelettransform.ChineseJGeophys.(in Chinese), 59(3): 884-896, doi: 10.6038/cjg20160311.

Rabbel W, MooneyW D. 1996. Seismic anisotropy of the crystalline crust: What does it tell us?.TerraNova, 8(1): 16-21.

Rümpker G, Kaviani A, Latifi K. 2014. Ps-splitting analysis for multilayered anisotropic media by azimuthal stacking and layer stripping.Geophys.J.Int.,199(1): 146-163.

Savage M K. 1998.Lower crustal anisotropy or dipping boundaries?Effects on receiver functions and a case study in New Zealand.J.Geophys.Res.,103(B7): 15069-15087.

Savage M K. 1999.Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting.Rev.Geophys.,37(1): 65-106. Schulte-Pelkum V, Mahan K H. 2014.A method for mapping crustal deformation and anisotropy with receiver functions and first results from USArray.EarthPlanet.Sci.Lett.,402: 221-233.

Shaw R, Srivastava S. 2007. Particle swarm optimization: A new tool to invert geophysical data.Geophysics, 72(2): F75-F83.

Sherrington H F, Zandt G, Frederiksen A. 2004. Crustal fabric in the Tibetan Plateau based on waveform inversions for seismic anisotropy parameters.J.Geophys.Res.,109: B02312.

Shi X M, Xiao M, Fan J K, et al. 2009. The damped PSO algorithm and its application for magnetotelluric sounding data inversion.ChineseJ.Geophys.(in Chinese), 52(4): 1114-1120, doi: 10.3969/j.issn.0001-5733.2009.04.029.

Shi Y, Eberhart R. 1998. A modified particle swarm optimizer. ∥Proceedings of the 1998 IEEE International Conference on Evolutionary Computation Proceedings,1998. IEEE World Congress on Computational Intelligence.Anchorage, AK:IEEE, 69-73. Silver P G, Chan W W. 1991.Shear wave splitting and subcontinental mantle deformation.J.Geophys.Res.,96(B10): 16429-16454. Simons F J, van der Hilst R D. 2003.Seismic and mechanical anisotropy and the past and present deformation of the Australian lithosphere.EarthPlanet.Sci.Lett.,211(3-4): 271-286.

Tian B F, Li J, Wang W M, et al. 2008.Crust anisotropy of Taihangshan mountain range in north China inferred from receiver functions.ChineseJ.Geophys.(in Chinese), 51(5): 1459-1467.

Vinnik L P, Farra V, Romanowicz B. 1989. Azimuthal anisotropy in the earth from observations of SKS at GEOSCOPE and NARS broadband stations.Bull.Seism.Soc.Am.,79(5): 1542-1558. Ward D, Mahan K, Schulte-Pelkum V. 2012. Roles of quartz and mica in seismic anisotropy of mylonites.Geophys.J.Int.,190(2): 1123-1134.

Weiss T, Siegesmund S, Rabbel W, et al. 1999. Seismic velocities and anisotropy ofthe lower continental crust: A review. ∥ Seismic Exploration of the Deep Continental Crust. Basel:Birkh?user, 97-122.

Xu Z, Xu M J, Wang L S, et al. 2006. A study on crustal anisotropy using P to S converted phase of the receiver function: Application to Ailaoshan-Red River fault zone.ChineseJ.Geophys.(in Chinese), 49(2): 438-448.

Yao C. 1989.The splitting of the teleseismic PS waves propagating through cracked solids.EarthquakeResearchinChina(in Chinese), 5(1): 38-47.

Yao C, Cai M G. 2009. Analytical expression of TI elastic tensor with arbitrary orientation.ChineseJ.Geophys.(in Chinese), 52(9): 2345-2348, doi: 10.3969/j.issn.0001-5733.2009.09.019.

Yao C, Hao C T, Zhang G L. 2016.A study of synthetic seismograms for SKS-wave response to crustal fractured-induce anisotropy.ChineseJ.Geophys.(in Chinese), 59(7): 2498-2509, doi: 10.6038/cjg20160715.

Yue B B, Peng Z M, Hong Y G, et al. 2009.Wavelet inversion of pre-stack seismic angle-gather based on particle swarm optimization algorithm.ChineseJ.Geophys.(in Chinese), 52(12): 3116-3123, doi: 10.3969/j.issn.0001-5733.2009.12.021.

Zandt G, Gilbert H, Owens T J, et al. 2004. Active foundering of a continental arc root beneath the southern Sierra Nevada in California.Nature, 431(7004): 41-46.

Zhu T, Li X F, Li Y Q, et al. 2011. Seismic scalar wave equation inversion based on an improved particle swarm optimization algorithm.ChineseJ.Geophys.(in Chinese), 54(11): 2951-2959, doi: 10.3969/j.issn.0001-5733.2011.11.025.

附中文參考文獻

房立華, 吳建平. 2009. 傾斜界面和各向異性介質(zhì)對接收函數(shù)的影響. 地球物理學進展, 24(1): 42-50.

劉啟元, 邵學鐘. 1985. 天然地震PS轉(zhuǎn)換波動力學特征的初步研究. 地球物理學報, 28(3): 291-302.

劉啟元, 王峻, 陳九輝等. 2007. 1976年唐山大地震的孕震環(huán)境: 密集地震臺陣觀測得到的結(jié)果. 地學前緣, 14(6): 205-213.

齊少華, 劉啟元, 陳九輝等.2009. 汶川MS8.0地震: 龍門山斷裂兩側(cè)地殼各向異性的初步研究. 地震地質(zhì), 31(3): 377-388.

齊少華, 劉啟元, 陳九輝等. 2016. 接收函數(shù)的曲波變換去噪. 地球物理學報, 59(3): 884-896, doi: 10.6038/cjg20160311.

師學明, 肖敏, 范建柯等. 2009. 大地電磁阻尼粒子群優(yōu)化反演法研究. 地球物理學報, 52(4): 1114-1120, doi: 10.3969/j.issn.0001-5733.2009.04.029.

田寶峰, 李娟, 王為民等. 2008. 華北太行山區(qū)地殼各向異性的接收函數(shù)證據(jù). 地球物理學報, 51(5): 1459-1467.

徐震, 徐鳴潔, 王良書等. 2006. 用接收函數(shù)Ps轉(zhuǎn)換波研究地殼各向異性──以哀牢山―紅河斷裂帶為例. 地球物理學報, 49(2): 438-448.

姚陳.1989. 穿透裂隙介質(zhì)遠震PS波的分裂. 中國地震, 5(1): 38-47.姚陳, 蔡明剛. 2009. 任意空間取向TI彈性張量解析表述. 地球物理學報, 52(9): 2345-2348, doi: 10.3969/j.issn.0001-5733.2009.09.019.

姚陳, 郝重濤, 張廣利. 2016. SKS波對地殼裂隙各向異性的響應——理論地震圖研究. 地球物理學報, 59(7): 2498-2509, doi: 10.6038/cjg20160715.

岳碧波, 彭真明, 洪余剛等. 2009. 基于粒子群優(yōu)化算法的疊前角道集子波反演. 地球物理學報, 52(12): 3116-3123, doi: 10.3969/j.issn.0001-5733.2009.12.021.

朱童, 李小凡, 李一瓊等. 2011. 基于改進粒子群算法的地震標量波方程反演. 地球物理學報, 54(11): 2951-2959, doi: 10.3969/j.issn.0001-5733.2011.11.025.

(本文編輯 胡素芳)

Stratified crustal anisotropy from receiver function and its particle swarm inversion

QI Shao-Hua, LIU Qi-Yuan*, CHEN Jiu-Hui, GUO Biao

StateKeyLaboratoryofEarthquakeDynamics,InstituteofGeology,ChinaEarthquakeAdministration,Beijing100029,China

The crustal anisotropy at multiple depths is essential for studying the vertical variation of the crustal deformation. Due to the complexity of the crust, the extraction of the crustal anisotropic parameters at different depths from teleseismic receiver function (RF) is an ongoing subject that requires further studies. Based on the previous works, this paper goes further to discuss the RF wavefield patterns of the stratified crustal anisotropy in terms of the generalized reflection and transmission coefficients method, which provides us a new theoretical basis for understanding the RF data with stratified anisotropic crustal models. In addition, we develop a global inversion method for extracting the stratified anisotropic model from RF via the particle swarm algorithm. The test results from both numerical and real data show that our method can correctly recover the anisotropic parameters of the multi-layered crustal model and distinguish the different anisotropic layers at each depth of the crust, as long as the isotropic velocity model is well determined. Noise suppression in RFs via the curvelet transform is valuable for the correct analysis of the crustal anisotropy at multiple depths.

Crustal structure; Stratified anisotropy; Receiver function; Particle swarm inversion

10.6038/cjg20161217.

國家自然科學基金(41590862)資助.

齊少華，男，1983年生，博士研究生，主要從事地震各向異性和接收函數(shù)研究.E-mail：ricown@163.com

*通訊作者 劉啟元，男，1945年生，研究員，主要從事寬頻帶地震學和地球深部結(jié)構(gòu)研究. E-mail：qyliu@ies.ac.cn

10.6038/cjg20161217

P315

2016-06-13，2016-08-28收修定稿

齊少華， 劉啟元， 陳九輝等. 2016. 地殼分層各向異性介質(zhì)接收函數(shù)及其粒子群反演. 地球物理學報，59(12):4544-4559,

Qi S H, Liu Q Y, Chen J H, et al. 2016. Stratified crustal anisotropy from receiver function and its particle swarm inversion.ChineseJ.Geophys. (in Chinese),59(12):4544-4559,doi:10.6038/cjg20161217.