2016, Vol. 59
地震各向异性与地球介质的矿物组成、温压条件、应力状态以及变形历史密切相关，它已成为推测地球深部介质变形的重要依据(Babuska and Cara, 1991; Rabbel and Mooney, 1996).早期研究认为，地球内部的地震各向异性主要表现在上地幔或者上地幔顶部，它源于地幔物质在高温高压环境下的塑性流动所形成的橄榄石晶格优势排列(Nicolas and Christensen, 1987; Savage, 1999).远震SKS分裂方法凭借其简洁和直观等优点，被广泛用于上地慢的方位各向异性研究(Vinnik et al., 1989; Silver and Chan, 1991).
近年来，日本Chugoku地区的远震P波接收函数的研究发现，莫霍面Ps转换震相的快波与慢波到时差可达到0.2~0.7 s (Nagaya et al., 2008)，表明地壳各向异性对SKS的影响不可忽略.理论地震图研究也证实，当地壳平均各向异性强度超过1%时，透射SKS波就可能发生分裂(姚陈等, 2016).因此，地壳各向异性研究不仅对于理解地壳变形机制等动力学问题至关重要，而且对于上地幔各向异性的正确解释也是极其必要的.
地壳各向异性的成因机理在不同深度上可能并不相同.上地壳的地震各向异性可能主要源于裂隙和微裂隙在应力驱动下的优势排列(Crampin, 1987)，而中下地壳各向异性则可能是由于韧性变形或塑性流动导致的矿物优势排列(Weiss et al., 1999; Savage, 1999).这意味着地壳各向异性具有分层的属性，以至于基于单层模型只能得到地壳各向异性较为粗略的估计.接收函数对地壳的速度间断面十分敏感，它在研究地壳各向异性及其分层结构方面具有难以替代的优势，现今已成为人们关注和研究的新热点(房立华和吴建平, 2009).
早在20世纪80年代末，姚陈(1989)就曾根据短周期远震体波记录上的莫霍面Ps转换波的分裂特征研究了地壳裂隙介质的各向异性.McNamara和Owens (1993)则利用接收函数(本文均指远震P波接收函数)莫霍面Ps转换波的分裂特征研究了美国盆岭地区的地壳各向异性.通过合成理论地震图的研究，Levin和Park (1997)证明，不仅是倾斜界面，地壳各向异性也可以导致接收函数切向分量的能量.利用有限差分算法计算的合成地震图，Okaya和McEvilly (2003)研究了各向异性对称轴倾斜对地震波场的影响.Savage (1998)、Frederiksen和Bostock (2000)、徐震等(2006)则分别研究了倾斜界面与各向异性的接收函数波场差异.
在任意取向横向各向同性(ATI)模型假定下，人们研究了地壳分层各向异性接收函数随反方位角的周期性变化特征(Frederiksen and Bostock, 2000; 田宝峰等, 2008).Nagaya等(2008)计算了双层地壳各向异性情况下，各向异性对称轴方位角不同的接收函数.根据观测数据的结果，齐少华等(2009)证明，水平对称轴模型可能并不适合解释变形复杂地区的地壳各向异性，在地壳范围内考虑ATI模型是极为必要的.但是，某些复杂的各向异性波场特征(如多层各向异性的相互影响、各向异性的沉积层和层间薄层界面的影响)，则尚未给予足够的关注.
尽管如此，各向异性接收函数的理论研究已为利用接收函数波形反演方法获取地壳不同层位的地震各向异性参数奠定了基础.各向异性的接收函数反演方法已在若干地区的地壳各向异性研究中取得了有价值的结果(Ozacar and Zandt, 2004; Sherrington et al., 2004; 田宝峰等, 2008).相较于剪切波分裂方法，基于波形反演的方法更适用于较为复杂的分层ATI介质.近年来，回避接收函数反演的各向异性参数提取方法也得到了发展(Liu and Niu, 2012; Rümpker et al., 2014; Schulte-Pelkum and Mahan, 2014)，并已引起了人们的兴趣.与接收函数反演方法一样，这些提取地壳各向异性参数的新方法主要依赖于单个台站接收函数(径向及切向分量)随反方位角的变化特征.
实际观测的接收函数存在介质横向非均匀散射的干扰，特别是接收函数的切向分量，它对壳内横向散射十分敏感(齐少华等, 2016).另外，地表沉积层的混响效应也会导致后续界面转换震相的信噪比降低.所有这些都增加了上述方法实际应用的难度.需要指出的是，现有的接收函数波场分析方法有赖于接收函数随反方位角变化的完整性.但是，由于台站(特别是流动台站)所处的地理位置，采集的接收函数在方位分布上通常是间断的，甚至可能存在大范围的缺失，这给本来就存在的界面倾斜与各向异性的折衷(trade-off)带来了更多的不确定性(Girardin and Farra, 1998; Liu and Niu, 2012; Rümpker et al., 2014).
本文的目的在于：(1) 针对已有工作的不足，进一步研究较为复杂的地壳分层各向异性，包括层间薄层界面及沉积盖层的接收函数；(2) 我们将引入曲波变换和压缩感知(齐少华等, 2016)，压制横向散射噪声，重构接收函数波场；(3) 在此基础上，通过数值和实际数据的检验，研究接收函数反演方法提取复杂分层结构不同层位各向异性参数的能力.这对于进一步理解各向异性接收函数波场和研究地壳分层各向异性都具有重要意义.
2 分层各向异性介质的接收函数 2.1 地壳各向异性的类型如图 1所示，快轴与慢轴各向异性分别代表两种不同类型的各向异性模型.根据Levin和Park (1997)的描述，快轴各向异性的相速度曲面为长球面，形态类似“西瓜”；而慢轴各向异性其相速度曲面为扁球面，形态类似“南瓜”.矿物晶体的优势排列导致快轴各向异性，而上地壳裂隙的优势排列和云母的面理组构(foliated mica fabric)会导致慢轴各向异性(Leidig and Zandt, 2003).岩石学证据业已证明，云母面理组构可能是中、下地壳各向异性的主要成因(Weiss et al., 1999; Lloyd et al., 2009; Ward et al., 2012).以往的研究通常根据实际情况，对各向异性的类型做出选择(Sherrington et al., 2004; 田宝峰等, 2008).但是，越来越多的证据表明，慢轴各向异性模型通常可以获得更小的数据拟合残差，似乎更为适合解释地壳内部的地震各向异性(Leidig and Zandt, 2003; Ozacar and Zandt, 2004; Zandt et al., 2004; Porter et al., 2011; Audet, 2015).因此，本文将统一采用慢轴各向异性模型来描述地壳各个层位的地震各向异性.
|
图 1 倾斜对称轴各向异性模型(据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. |
为了研究地壳分层各向异性，本文将依据广义反射透射系数矩阵方法(Chen, 1993)进行理论接收函数的计算，并采用作者自行开发的程序代码.这有助于考虑地壳复杂分层模型的完整波场.本文算法与Levin和Park (1997)没有实质性区别，但在以下方面采取了不同的处理：在构建各向异性弹性参数时，我们采用了姚陈和蔡明刚(2009)推导的ATI弹性张量解析公式，在一定程度上提高了计算的精度.
根据弱各向异性假定(Backus, 1965; Park and Yu, 1992)，各向异性弹性参数可以近似地分为各向同性的背景值(由P波速度、S波速度和密度来描述)和扰动值(由各向异性对称轴的方位角、倾角、P波和S波各向异性强度表征)，这种参数化方法便于正演模型的构建，减少反演参数(Levin and Park, 1997).
2.3 快轴与慢轴各向异性接收函数的比较图 2给出了采用本文方法计算的单层地壳模型的快轴与慢轴各向异性的理论接收函数.所用的水平慢度为0.0618 s·km-1，这大体相当于震中距为60°的情况.表 1给出了相应的地壳模型参数.其中，模型M1和M2分别为快轴各向异性(强度为正)和慢轴各向异性(强度为负).除了快轴与慢轴的方位角相差180°以外，其他参数没有区别.
|
图 2 快轴与慢轴各向异性模型的接收函数波场 (a) 快轴：径向分量；(b) 快轴：切向分量；(c) 慢轴：径向分量；(d) 慢轴：切向分量；红色表示极性为正，蓝色表示极性为负.所有符号及颜色的表示在下文各图中均相同. Fig. 2 Receiver functions from the fast-axis and slow-axis anisotropic model (a) Fast-axis: R-component; (b) Fast-axis:T-component; (c) Slow-axis: R-component; (d) Slow-axis:T-component; Red color denotes positive amplitude, blue color negative amplitude. All symbols and colors have the same meaning in the following figures. |
|
表 1 快轴与慢轴各向异性模型 Table 1 Fast-axis and slow-axis anisotropic model |
比较图 2a和2b可知，模型底部界面的透射Ps转换震相及多次波的振幅及相位随反方位角变化的周期性是一样的：均以0°和180°反方位角为轴，径向分量呈现对称，切向分量呈现反对称.但是，它们的振幅及相位随反方位角变化的具体特征是不同的，特别是，切向分量的多次波呈现了相反的极性.这些差异或许可以作为我们判别各向异性类型的依据.需要说明的是，快轴与慢轴方位角相差180°与倾角互补是等价的(Levin and Park, 1997; Erickson, 2002).
2.4 地壳分层各向异性接收函数接收函数是观测台站接收区底部入射波的脉冲响应.这意味着，对于分层各向异性的地壳，各界面的透射与反射波均会涉及所有层位的各向异性.因此，对于较为复杂的分层结构，能否从接收函数正确地解析各个层位的各向异性参数无疑是一个值得关注的问题.另外，对于构造活动区，地壳通常都存在不同程度的运动变形，并伴随着壳内的低速层(如青藏高原及我国的华北地区).在应力作用和应变调节下，它们有可能成为地壳物质流动的通道(刘启元等, 2007; Beaumont et al., 2001; Liu et al., 2014)，意味着低速层可能具有较强的各向异性.刘启元和邵学钟(1985)曾证明，薄层的混响效应会导致相应的转换震相能量的增强.如果它同时显示了明显的各向异性特征，那么是否会误导分层各向异性参数的估计，也是一个有待研究的问题.因此，更为复杂的地壳分层各向异性将是本文关注的重点.我们将通过计算分层各向异性模型的接收函数随反方位角的变化，来考察不同层位各向异性参数对接收函数波场的影响.
我们首先考虑双层各向异性的情况，表 2给出了相应的地壳模型.首先，我们考虑上下两层各向异性参数相同的情况(模型M3).图 3a和3b分别给出了计算得到的接收函数径向和切向分量随反方位角的变化.比较图 2c-2d(单层地壳)和图 3a-3b(双层地壳)，可有如下观察：(1) 由于层间界面转换波(Ps1)及多次波(PpPs1和PsPs1+PpSs1)的出现，并不难做出双层地壳结构的判断；(2) 切向分量的初至P波投影的极性相同；(3) 层间界面Ps1震相的径向分量与其切向分量均显示了近似各向同性界面的特征；(4) 在径向和切向分量上，层间界面Ps1震相的多次波的各向异性特征明显；(5) 两者底部界面的转换波和多次波径向及切向分量没有明显差别.上述观察表明，当双层各向异性参数相同且不考虑多次波时，其与单层模型的各向异性特征基本相同.
|
|
表 2 双层各向异性模型 Table 2 Two-layered anisotropic model |
|
图 3 双层各向异性模型的接收函数波场：上下层各向异性参数相同 (a)径向分量；(b)切向分量. Fig. 3 RF-wavefield from two-layered anisotropic model: Same anisotropic parameters (a) R-component; (b) T-component. |
下面，我们考虑各向异性参数不同的双层地壳模型.我们先固定上层，仅改变下层的各向异性参数(模型M4-M6).图 4给出了相应的接收函数随反方位角的变化.
图 4a和4b分别给出了仅对称轴倾角不同时(M4)的接收函数径向与切向分量.与图 3a-3b相比，下层各向异性对称轴倾角的改变使接收函数径向及切向分量的各向异性特征(特别是层间界面Ps1震相)在保持周期性不变的情况下得到增强，意味着层间界面震相各向异性特征的强弱主要取决于各向异性参数的层间差异.
|
图 4 双层各向异性模型的接收函数波场：上层参数固定，下层参数改变 (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. |
图 4c和4d分别给出了仅对称轴方位角正交(M5)的接收函数径向与切向分量.与图 3a-3b和图 4a-4b对比，不难发现：(1) 在切向分量上，零时刻P波投影没有因为下层各向异性参数的改变而改变(如图 4d两个橙色箭头所示)；(2) 无论改变下层各向异性对称轴的倾角，还是方位角，层间界面的多次波(PpPs1和PsPs1+PpSs1)的径向分量均未有视觉上的变化；(3) 下层各向异性对称轴方位角的改变所引起的底部界面转换波(Ps2)的变化明显大于倾角改变所引起的变化，因为对称轴方位角的改变导致其对称位置平移了90°(如图 4d中的绿色箭头所示)；(4) 底部界面的多次波(PpPs2和PsPs2+PpSs2)的径向及切向分量的各向异性特征异常削弱.上述观察表明，接收函数波场对各向异性对称轴的方位角更为敏感，它可以引起接收函数随反方位角变化的对称位置的改变.特别是，当上下两层的各向异性对称轴的方位角正交时，下层界面转换波的各向异性特征会被严重削弱.
图 4e和4f分别给出了对称轴方位角正交且倾角不同时(M6)的接收函数径向与切向分量.与模型M5的情况相比，其波场特征未有明显变化，正如模型M3与M4的接收函数波场也没有明显差异的情况类似.这意味着接收函数波场的各向异性特征对倾角变化确实并不敏感.
表 2中模型M7-M9给出了固定下层各向异性参数，仅改变上层参数时的双层地壳模型，其接收函数波场如图 5所示.对比图 4和图 5，不难发现：对于双层地壳各向异性，(1) 仅改变下层的各向异性参数与仅改变上层各向异性参数的接收函数波场具有十分类似的现象，两者均表明，接收函数各向异性特征的变化对介质各向异性对称轴的方位角比对其倾角更为敏感；(2) 尽管图 4和图 5的波场特征明显不同，但两者呈现了某种程度的“互补性”，即与改变下层各向异性参数相反，上层各向异性参数的改变导致该层位相应震相的各向异性特征更为明显地改变.例如，与图 4不同，图 5橙色箭头指示了上层各向异性轴方位的改变，导致初至P波投影随反方位角变化的对称位置移动了90°，而图 5中绿色箭头指示的对称位置与图 3并没有什么不同，尽管下层各向异性参数将同时影响其上下两个界面的透射转换震相(图 4中的Ps1和Ps2)；(3) 无论改变上层，还是下层的各向异性参数，接收函数的切向分量较其径向分量对各向异性参数的改变更为敏感.
|
图 5 双层各向异性模型的接收函数波场：下层参数固定，上层参数改变 (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为各向异性低速薄层的上下均为各向同性介质的地壳模型，相应接收函数的径向和切向分量分别由图 6a和6b给出.观察可知，薄层界面的转换震相Ps1在径向和切向分量上均显示了清晰的各向异性特征，但该薄层界面上的反射波(PpPs1和PsPs1+PpSs1)的各向异性效应非常微弱，薄层对来自底部界面的震相Ps2并未形成可以察觉的干扰.这意味着，对于各向同性的双层地壳来说，各向异性薄层界面对接收函数波场的影响十分有限.值得关注的是，与各向同性薄层不同(刘启元和邵学钟, 1985)，在径向分量反方位角180°附近，各向异性导致了层间薄层界面转换震相能量的离散.
|
|
表 3 多层各向异性模型 Table 3 Multi-layered anisotropic model |
|
图 6 包含低速薄层的各向异性模型的接收函数波场 (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. |
模型M11是一个各向同性低速薄层与其上下层均为各向异性层的地壳模型，图 6c和6d分别给出了相应的接收函数径向和切向分量.模型M12与模型M11的不同是其薄层界面为各向异性低速薄层，图 6e和6f给出了M12的接收函数径向及切向分量.对比模型M11和M12的接收函数可知，两者的差别不大，但薄层界面的转换波随反方位角的变化特征是有差别的，这意味着接收函数波场的反方位完整性对于判别薄层是否存在各向异性十分重要.这也提示我们，根据有限反方位的接收函数波场只能对地壳分层各向异性做出粗略的估计.
模型M13是一个底部为过渡层(厚度8 km)，其上层为各向同性层的地壳模型.图 7a和7b分别给出了相应的接收函数径向及切向分量.由图 7a和7b可见，由于各向异性的存在，该过渡层上的转换震相(黑色箭头所示的震相)在反方位角180°附近出现了离散，而由于上行(下行)的次数增加，其相应多次波在过渡层顶部与底部彻底分开，且没有显现明显的各向异性效应.综合低速薄层(图 6a-6b)和过渡层(图 7a-7b)的接收函数波场可知，其各向异性对波场的影响局限于透射转换波，而多次波对其各向异性并不敏感，这与厚层介质的情况完全不同.图 7a-7b和图 2c-2d差异提供的证据表明，厚层和薄层各向异性对接收函数波场的影响是相互独立的.但是，在提取各向异性参数时，考虑多次反射与采用完整反方位的接收函数对于避免两者混淆均具有重要作用.
|
图 7 包含沉积层和过渡层的各向异性模型的接收函数波场 (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. |
模型M14是在模型M10的基础上增加了地表 2 km沉积薄层后的地壳模型.图 7c和7d分别给出了相应接收函数的径向和切向分量.图 7c和图 7d表明，沉积薄层对接收函数的径向和切向分量均产生了极为显著的影响.尽管如此，对比图 6a-6b可知，沉积薄层与低速薄层的各向异性特征(包括其对称性和周期性)依然可在切向分量上加以区分.
3 各向异性介质接收函数的粒子群反演本节将进一步研究从接收函数提取地壳分层各向异性参数的能力.如前所述，接收函数反演是提取地壳分层各向异性参数的常用方法.前人广泛采用了遗传、近邻和模拟退火等全局性算法，以避免陷于局部极小(Ozacar and Zandt, 2004; Sherrington et al., 2004; 田宝峰等, 2008).与各向同性介质的接收函数反演不同，地壳各向异性参数的提取需要对不同反方位角的接收函数径向与切向波形进行联合反演.因而，反演算法的效率，特别当针对大量台站的研究时，成为一个不可回避的问题.
近年来，一种称为“粒子群优化”(particle swarm optimization，简称PSO)的全局算法得到了发展和完善(Eberhart and Kennedy, 1995; Shi and Eberhart, 1998; Clerc and Kennedy, 2002).该算法基于鸟群和鱼群觅食过程中表现出的群体智能(swarm intelligence)，具有实现简单、易于并行、收敛速度快以及需要调整参数少等优势(Eberhart and Shi, 1998; Hassan et al., 2005).该方法在国内外正逐步得到日益广泛的关注和应用(Shaw and Srivastava, 2007; Martínez et al., 2010; 师学明等, 2009; 岳碧波等, 2009; 朱童等, 2011).我们将把该算法引入到地壳分层各向异性的接收函数反演.
在反演算法中，正演计算采用前文(2.2节)所述的方法，粒子的位置即为模型空间里的一个样点(Eberhart and Kennedy, 1995)，反演的目标函数则基于时间域波形互相关(Frederiksen et al., 2003).对于PSO算法，单个粒子依据其自身和其他粒子的飞行经验，调整自身的飞行速度，并以迭代方式向全局最优解的位置逼近.在迭代过程中，粒子既不会消亡，也不会被丢弃.粒子自身和其他粒子的飞行经验在算法执行的过程中被完整地保留，并指导着粒子的后续飞行.正是这种个体与自身、个体与群体之间信息分享的拓扑结构，使得PSO算法同时兼备全局与局部搜索的特性.
为了权衡PSO全局与局部搜索的能力，本文采用Shi和Eberhart (1998)发展的粒子速度更新方程.同时，为避免迭代中后期出现所谓的“粒子爆炸”现象(即粒子由于速度过大而逃离搜索空间)，引入“收缩因子”来避免“粒子爆炸”，同时使算法收敛速度进一步加快(Clerc and Kennedy, 2002).
理论上，PSO的收敛判据应该是粒子群中所有粒子均重叠在模型空间的同一个点上.然而，由于粒子存在惯性，且不同粒子的运行轨迹有很大的差异，所有粒子最终静止在同一点上可能需要很长的时间.为此，人们把所有粒子均进入一个合理的误差范围作为反演终止的判据，或者直接规定迭代的最大次数.
表 4给出了数值检验所使用的“真实”接收区模型.计算所用的水平慢度设定为0.0618 s·km-1.由表 4可知，我们在模型中考虑了低速层、高速层以及薄低速盖层.需要说明的是，在数值检验中，(1) 我们假定各向同性的速度模型已知，且各向异性参数的分层与速度参数保持一致，即我们仅把各向异性参数作为独立变量；(2) 初始各向异性参数均设为零；(3) 对反演参数作如下约束：P波与S波各向异性强度保持相同，取值范围限定在[-20%, 0]，各向异性对称轴的倾角限定在[0°, 90°]，方位角限定在[0°, 360°]；(4) 粒子数量设定为50，最大迭代次数为100次.
|
|
表 4 合成数据检验所用的各向异性模型 Table 4 Anisotropic model parameters used for synthetic data test |
图 8给出了数值检验结果.比较图 8a和8b可知，反演得到的接收函数波场与“真实”数据没有视觉上的差别.图 8f中红色圆圈描述的收敛曲线表明，大约经过50次左右的迭代后，目标函数的残差即已基本收敛.
|
图 8 各向异性接收函数反演的数值检验 (a) “真实”接收函数；(b) 预测接收函数；(c) 各向异性对称轴的方位角；(d) 各向异性对称轴的倾角；(e) 各向同性速度结构与各向异性强度；(f) “真实”与预测接收函数之间的残差. 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. |
比较表 4和图 8(c, d, e)可知，反演结果与“真实”的各向异性模型基本一致.这表明在各向同性速度模型已知，且不含噪声的条件下，利用本文方法可以从接收函数波场可靠地提取各个层位的各向异性参数.
与上述数值检验不同，实际观测数据难以避免横向非均匀构造的散射及各种噪声.这往往成为从观测接收函数中提取地震各向异性参数的难点(Savage, 1998).与理论的各向异性介质接收函数波场所展示的随反方位角的可追踪性不同，横向非均匀散射及各种噪声在空间取向上可能具有随机性.齐少华等(2016)给出的结果表明，利用曲波变换对接收函数波场进行二维去噪不仅可以明显地改善震相的可追踪性，而且可同时通过压缩感知波场重建技术获得全方位的接收函数，无疑有助于分层各向异性介质接收函数反演研究.
齐少华等(2016)曾利用曲波去噪方法处理过位于澳大利亚东北地区的全球数字地震台网(IRIS)台站CTAO (-20.0882°N，146.2545°E，台网编号IU)的接收函数.该台站所处的地理位置和构造背景如图 9所示.在此基础上，我们将进一步利用本文的分层各向异性介质接收函数反演技术，提取该台下方的地壳各向异性参数，其目的在于检验各向异性接收函数反演解释实际观测数据的能力.
|
图 9 台站CTAO的地理位置及其构造背景(修改自Ford等, 2010) 黑色实线勾画出了澳大利亚的主要构造区域；西部斜线区域代表太古宙克拉通；北部和南部横线区域代表元古宙克拉通；东部点阵区域显示了显生宙基底；黑色三角表示IRIS台网的台站，红色三角表示本文使用的台站. 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)已经描述了利用时间域迭代反褶积方法(Ligorría and Ammon, 1999)得到的151个远震事件(1999年2月-2015年7月，震中距30°≤Δ≤90°，Mb > 5.5)在CTAO台的接收函数.为此，本文不再展示未经曲波去噪处理的接收函数，仅简单综述基本的数据参数如下：(1) 利用时间域迭代反褶积方法提取接收函数时，采用的高斯系数为2.5，这大体上相当于0~1.2 Hz的带通滤波；(2) 波场重建的反方位角间隔为1°.图 10a给出了用于本文反演研究的接收函数，它是以10°反方位角间隔进行进一步分组叠加后的结果.
|
图 10 台站CTAO的地壳各向异性反演 (a) 观测接收函数；(b) 反演得到的接收函数；(c) 分层各向异性对称轴方位角；(d) 分层各向异性对称轴倾角；(e) 各向同性速度结构与各向异性强度；(f) 观测与预测接收函数之间的残差, (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. |
本文的各向异性介质接收函数反演分为以下步骤：(1) 首先对全方位叠加波形进行各向同性接收函数反演，以便确定台站下方地壳上地幔速度结构；(2) 对获得的地壳速度结构进行必要的简化，形成各向异性参数的分层结构，且在各向异性反演中固定不变.需要说明的是，速度模型中莫霍界面以下的部分不参与各向异性反演.
图 10给出了利用本文方法得到的反演结果.表 5给出了台站CTAO下方地壳速度结构及反演得到的分层各向异性参数.由表 5可知，该台下方的地壳厚度为36 km.Chevrot和van der Hilst (2000)利用接收函数H-κ扫描得到的地壳厚度为34 km，与我们的结果基本一致.
|
|
表 5 台站CTAO下方地壳各层的各向异性慢轴倾角和方位角 Table 5 The plunge and trend of the slow-axis crustal anisotropy beneath Station CTAO |
由图 10b可知，与观测的接收函数波场相比，前5 s窗口内，反演的接收函数径向与切向分量随反方位角的变化特征与观测数据吻合较好.需要说明的是，我们所用的目标函数以波形互相关为基础，这意味着主要强调相位拟合.观测与预测接收函数波场存在一定差异，一方面是由于水平分层地壳模型无法解释三维非均匀介质的波场；另一方面，则可能因为在反演中我们没有考虑莫霍界面以下上地幔的各向异性.台站下方各层的各向异性强度如图 10e所示.结果表明，台站CTAO下方地壳浅部各向异性强度较弱，而深部各向异性强度较强，其强度的大小与前人在其他地区得到的结果具有较好的可比性(Sherrington et al., 2004; 田宝峰等, 2008)，我们认为这是合理的.
图 10c和10d分别给出了台站CTAO下方地壳各层位各向异性慢轴的方位角及倾角，具体数值列于表 5.由于没有关于台站CTAO下方地壳各层位各向异性参数的研究报道，我们无法提供能够证明本文结果的直接可比证据.尽管如此，我们可以根据相关的间接证据来说明本文给出的该台下方各层位慢轴方位角及倾角的合理性.为此，我们给出如下讨论：
(1) 台站CTAO地壳慢轴各向异性倾角由浅到深逐渐趋于垂直方向，意味着中下地壳的各向异性组构(anisotropic fabric)可能趋于水平，且水平方向的波速快于垂直方向，这与Debayle和Kennett (2000)得到的岩石圈200 km以上SH波速明显快于SV的面波径向各向异性结果相一致.
(2) 台站CTAO的下地壳(第4和5层)慢轴各向异性方位角为近EW方向，意味着与其正交的快波方向为近NS方向，与根据瑞利波方位各向异性所得到的快波方向(80~100 km以上)基本一致(Debayle and Kennett, 2000; Simons and van der Hilst, 2003).
4 结论与讨论本文利用广义反射透射系数矩阵方法计算的合成地震图，研究了复杂地壳分层各向异性接收函数随反方位角的变化和不同层位各向异性参数对接收函数波场的影响，为各向异性介质接收函数的解释提供了新的理论依据；通过引入粒子群优化理论，发展了地壳分层各向异性接收函数反演方法，并在此基础上，通过数值和观测数据的检验，论证了分层各向异性接收函数反演提取地壳不同层位各向异性参数的能力.根据我们的结果，可以得到以下结论：
(1) 分层各向异性地壳接收函数的理论波场研究表明，快轴与慢轴模型的接收函数随反方位角变化的特征是不同的；接收函数波场的各向异性特征与各向异性对称轴方位角的关系较其倾角更为密切；接收函数切向分量较其径向分量对各向异性参数更为敏感；层间界面Ps震相随反方位角的变化特征是界面上下两层共同作用的结果，层间界面Ps震相各向异性特征的强弱主要取决于各向异性参数的层间差异；沉积盖层的混响效应可以对接收函数(径向和切向分量)造成干扰；各向异性薄层界面的转换震相可以导致震相的离散，但其多次波的各向异性效应微弱.本文的结果表明，偏振分析方法对于分层各向异性研究存在一定局限性.
(2) 数值检验的结果表明，在地壳各向同性速度结构确定的前提下，分层各向异性介质的粒子群接收函数反演可以正确提取各层位的各向异性参数；把多次反射纳入反演有利于正确解析不同层位的各向异性参数；由于接收函数对地壳介质的绝对速度并不敏感，在各向异性接收函数反演前，独立确定接收区各向同性速度结构，有助于分层各向异性接收函数反演的稳定性.
(3) 观测数据的验证结果表明，曲波变换去噪降低了横向非均匀散射和其他噪声的干扰，有效地完善了接收函数方位的完整性，这对正确解析不同层位的各向异性参数和避免接收函数反演中不同层位各向异性参数的多解性都具有重要价值.
(4) 虽然曲波变换去噪有效降低了横向非均匀散射和其他噪声的干扰，但并不意味着完全去除了三维地壳结构对接收函数波场的影响.理论上，把多次反射纳入接收函数反演有利于正确解析不同层位各向异性参数，但由于其路径放大效应，对于基于水平分层介质的接收函数反演来说，多次反射的各向异性特征解析仍然十分困难.
致谢作者感谢评阅人对本文提出的宝贵建议，这使得本文得到进一步完善.本文所用数据来源于IRIS全球地震台网.感谢Brecht Donckels博士提供粒子群优化的核心算法.
| Audet P. 2015. Layered crustal anisotropy around the San Andreas Fault near Parkfield, California. J. Geophys. Res., 120(5): 3527-3543. DOI:10.1002/2014JB011821 | |
| 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. DOI:10.1029/JZ070i014p03429 | |
| 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. DOI:10.1038/414738a | |
| Chen X F. 1993. A systematic and efficient method of computing normal modes for multilayered half-space. Geophys.J. Int., 115(2): 391-409. DOI:10.1111/gji.1993.115.issue-2 | |
| Chevrot S, van der Hilst R D. 2000. The Poisson ratio of the Australian crust: geological and geophysical implications. Earth Planet.Sci. Lett., 183(1-2): 121-132. DOI:10.1016/S0012-821X(00)00264-8 | |
| Clerc M, Kennedy J. 2002. The particle swarm-explosion, stability, and convergence in a multidimensional complex space. IEEE Transactions on Evolutionary Computation, 6(1): 58-73. DOI:10.1109/4235.985692 | |
| Crampin S. 1987. Geological and industrial implications of extensive-dilatancy anisotropy. Nature, 328(6130): 491-496. DOI:10.1038/328491a0 | |
| Debayle E, Kennett B L N. 2000. Anisotropy in the Australasian upper mantle from Love and Rayleigh waveform inversion. Earth Planet.Sci. Lett., 184(1): 339-351. DOI:10.1016/S0012-821X(00)00314-9 | |
| 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. Progress in Geophys., 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. Earth Planet.Sci. Lett., 300(3-4): 299-310. DOI:10.1016/j.epsl.2010.10.007 | |
| Frederiksen A W, Bostock M G. 2000. Modelling teleseismic waves in dipping anisotropic structures. Geophys.J. Int., 141(2): 401-412. DOI:10.1046/j.1365-246x.2000.00090.x | |
| 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. DOI:10.1046/j.1365-246X.2003.02043.x | |
| 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. DOI:10.1046/j.1365-246X.1998.00525.x | |
| 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. DOI:10.1111/gji.1997.131.issue-2 | |
| 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. DOI:10.1111/gji.2012.188.issue-1 | |
| Liu Q Y, Shao X Z. 1985. Study on the dynamic characteristics of PS converted waves. Acta Geophysica Sinica, 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. Earth Science Frontiers, 14(6): 205-213. DOI:10.1016/S1872-5791(08)60012-3 | |
| 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. Nature Geoscience, 7(5): 361-365. DOI:10.1038/ngeo2130 | |
| Lloyd G E, Butler R W H, Casey M, et al. 2009. Mica, deformation fabrics and the seismic properties of the continental crust. Earth Planet.Sci. Lett., 288(1-2): 320-328. DOI:10.1016/j.epsl.2009.09.035 | |
| 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. Journal of Applied Geophysics, 71(1): 13-25. DOI:10.1016/j.jappgeo.2010.02.001 | |
| 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. DOI:10.1029/93JB00711 | |
| 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. DOI:10.1785/0120080130 | |
| 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. DOI:10.1046/j.1365-246X.2003.01896.x | |
| 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. DOI:10.1111/gji.1992.110.issue-3 | |
| 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. DOI:10.1130/L126.1 | |
| Qi S H, Liu Q Y, Chen J H, et al. 2009. Wenchuan earthquake Ms8.0: Preliminary study of crustal anisotropy on both sides of the Longmenshan faults. Seismology and Geology, 31(3): 377-388. | |
| Qi S H, Liu Q Y, Chen J H, et al. 2016. Attenuation of noise in receiver functions using curvelettransform. Chinese J Geophys., 59(3): 884-896. DOI:10.6038/cjg20160311 | |
| Rabbel W, Mooney W D. 1996. Seismic anisotropy of the crystalline crust: What does it tell us?. Terra Nova, 8(1): 16-21. DOI:10.1111/ter.1996.8.issue-1 | |
| 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. DOI:10.1093/gji/ggu154 | |
| 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. DOI:10.1029/98JB00795 | |
| Savage M K. 1999. Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting. Rev. Geophys., 37(1): 65-106. DOI:10.1029/98RG02075 | |
| Schulte-Pelkum V, Mahan K H. 2014. A method for mapping crustal deformation and anisotropy with receiver functions and first results from USArray. Earth Planet.Sci. Lett., 402: 221-233. DOI:10.1016/j.epsl.2014.01.050 | |
| Shaw R, Srivastava S. 2007. Particle swarm optimization: A new tool to invert geophysical data. Geophysics, 72(2): F75-F83. DOI:10.1190/1.2432481 | |
| 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. Chinese J. Geophys., 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. DOI:10.1029/91JB00899 | |
| Simons F J, van der Hilst R D. 2003. Seismic and mechanical anisotropy and the past and present deformation of the Australian lithosphere. Earth Planet.Sci. Lett., 211(3-4): 271-286. DOI:10.1016/S0012-821X(03)00198-5 | |
| Tian B F, Li J, Wang W M, et al. 2008. Crust anisotropy of Taihangshan mountain range in north China inferred from receiver functions. Chinese J. Geophys., 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. DOI:10.1111/gji.2012.190.issue-2 | |
| 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. Chinese J. Geophys., 49(2): 438-448. | |
| Yao C. 1989. The splitting of the teleseismic PS waves propagating through cracked solids. Earthquake Research in China, 5(1): 38-47. | |
| Yao C, Cai M G. 2009. Analytical expression of TI elastic tensor with arbitrary orientation. Chinese J. Geophys., 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. Chinese J. Geophys., 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. Chinese J. Geophys., 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. DOI:10.1038/nature02847 | |
| Zhu T, Li X F, Li Y Q, et al. 2011. Seismic scalar wave equation inversion based on an improved particle swarm optimization algorithm. Chinese J. Geophys., 54(11): 2951-2959. DOI:10.3969/j.issn.0001-5733.2011.11.025 | |
| 房立华, 吴建平. 2009. 倾斜界面和各向异性介质对接收函数的影响. 地球物理学进展, 24(1): 42–50. | |
| 刘启元, 邵学钟. 1985. 天然地震PS转换波动力学特征的初步研究. 地球物理学报, 28(3): 291–302. | |
| 刘启元, 王峻, 陈九辉, 等. 2007. 1976年唐山大地震的孕震环境:密集地震台阵观测得到的结果. 地学前缘, 14(6): 205–213. | |
| 齐少华, 刘启元, 陈九辉, 等. 2009. 汶川MS8.0地震:龙门山断裂两侧地壳各向异性的初步研究. 地震地质, 31(3): 377–388. | |
| 齐少华, 刘启元, 陈九辉, 等. 2016. 接收函数的曲波变换去噪. 地球物理学报, 59(3): 884–896. DOI:10.6038/cjg20160311 | |
| 师学明, 肖敏, 范建柯, 等. 2009. 大地电磁阻尼粒子群优化反演法研究. 地球物理学报, 52(4): 1114–1120. DOI:10.3969/j.issn.0001-5733.2009.04.029 | |
| 田宝峰, 李娟, 王为民, 等. 2008. 华北太行山区地壳各向异性的接收函数证据. 地球物理学报, 51(5): 1459–1467. | |
| 徐震, 徐鸣洁, 王良书, 等. 2006. 用接收函数Ps转换波研究地壳各向异性──以哀牢山-红河断裂带为例. 地球物理学报, 49(2): 438–448. | |
| 姚陈. 1989. 穿透裂隙介质远震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. 基于粒子群优化算法的叠前角道集子波反演. 地球物理学报, 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 | |