地球物理学报  2018, Vol. 61 Issue (9): 3689-3700   PDF    
一种改进的适应于倾斜Moho面的H-κ-θ叠加方法及应用
谭萍1,2,3, 陈赟1,2, 孙维昭4, 李玮1,2, 唐国彬5, 崔田1,2,3     
1. 中国科学院地质与地球物理研究所, 岩石圈演化国家重点实验室, 北京 100029;
2. 中国科学院地球科学研究院, 北京 100029;
3. 中国科学院大学, 北京 100049;
4. 中国石油集团东方地球物理勘探有限责任公司研究院, 河北涿州 072750;
5. 桂林理工大学地球科学学院, 广西桂林 541006
摘要:地壳厚度、波速比或泊松比是研究地壳结构和性质的基本地震学参数,对于研究地壳组分特征及构造演化具有重要意义.经典H-κ叠加方法是利用远震接收函数资料求取台站下方地壳厚度和波速比最为简便高效的方法.但该方法隐含着Moho面是水平界面的假设条件,意味着Ps转换波及后续多次波相对P波的走时主要取决于地壳厚度和纵横波速度,而忽略了界面产状的影响.理论模拟表明,如果不考虑Moho面的产状,特别是在Moho面倾角较大的情况下,利用经典H-κ叠加方法得到的地壳厚度和波速比会偏离实际模型,尤其会造成波速比的过高估计,从而影响到对地壳结构和性质的认知.为了解决Moho面倾斜条件下的地壳厚度和波速比求取问题,本文推导了界面倾斜条件下Ps转换波与后续多次波相对于直达P波的理论到时公式;基于经典H-κ叠加方法的思想,提出了一种可以同步求取地壳厚度-波速比-Moho面倾角的H-κ-θ叠加方法.通过理论模型测试,验证了该方法具备较高的稳定性和可靠性,并将此方法应用于青藏高原南部Hi-CLIMB台阵资料,显示出较好的应用潜力.
关键词: 地壳厚度      波速比      P波接收函数      H-κ叠加方法      倾斜Moho面     
An improved H-κ-θ stacking method to determine the crustal thickness and bulk vP/vS ratios in the case of a slant Moho interface
TAN Ping1,2,3, CHEN Yun1,2, SUN WeiZhao4, LI Wei1,2, TANG GuoBin5, CUI Tian1,2,3     
1. State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
2. Institutions of Earth Science, Chinese Academy of Sciences, Beijing 100029, China;
3. University of Chinese Academy of Sciences, Beijing 100049, China;
4. Geophysical Research Institute of BGP Inc., CNPC, Zhuozhou Hebei 072750, China;
5. College of Earth Sciences, Guilin University of Technology, Guilin Guangxi 541006, China
Abstract: Crustal thickness and vP/vS ratio (or Poisson's ratio) are the key seismological parameters for studying the crustal composition and properties. The H-κ stacking method is widely used as one of the most powerful techniques to determine the two parameters beneath a seismic station using teleseismic P-wave receiver functions. However, the H-κ stacking method assumes the Moho interface is flat, and thus the arrival times of converted Ps waves and subsequent multiples relative to the direct P-wave depend mainly on the crustal thickness and P-and S-wave velocities. The synthetic test shows that the crustal thickness and vP/vS ratios obtained by the H-κ stacking method would greatly deviate from the real model when the Moho interface is slant, and this deviation increases with the increasing dip angle of the Moho interface. In particular, this deviation will lead to the overestimation of vP/vS ratios, and eventually mislead the interpretation of the crustal structure and properties. In this study, the formulas of the theoretical arrival times of converted Ps waves and multiples between the slant Moho interface and the ground surface are derived, and then an improved H-κ-θ stacking method is developed to determine the crustal thickness, bulk vP/vS ratios and dip angle of the Moho interface (H-κ-θ) simultaneously based on the idea of the H-κ stacking method. The synthetic tests validate the good stability and reliability of the improved H-κ-θ stacking method. Applying the improved method to the Hi-CLIMB array data in southern Tibet further shows that the H-κ-θ stacking method has the good practicability for real seismic data.
Keywords: Crustal thickness    vP/vS ratios    P-wave receiver functions    H-κ stacking method    Slant Moho interface    
0 引言

地壳厚度、波速比(vP/vS)或泊松比是研究地壳结构和性质的基本地震学参数,对于研究地壳组分特征及构造演化具有重要意义(Christensen, 1996; 嵇少丞等, 2009).经典H-κ叠加方法是利用远震接收函数资料求取台站下方地壳厚度和波速比最为简便高效的方法(Zhu and Kanamori, 2000).该方法通过利用远震P波接收函数径向分量的Moho面转换横波(Ps)及后续多次波(Pps和Pss)相对远震P波的走时和振幅信息,联合求取地壳厚度和波速比,较好地解决了地壳厚度与波速比之间的折衷(trade-off)问题,具有数据自驱动程度高、实现简便、准确高效等优点,至今已作为接收函数常规处理的通用方法而得到普遍应用(例如,Chen et al., 2013, 2015).但H-κ叠加方法隐含着Moho面是水平界面的假设条件,即意味着Ps及后续多次波相对P波的走时主要取决于地壳厚度和纵横波速度,而忽略了界面产状的影响(Zhu and Kanamori, 2000).理论模拟表明,如果不考虑Moho面的产状,特别是在Moho面存在较大倾角的情况下,利用经典H-κ叠加方法得到的地壳厚度和速度比会在较大范围偏离实际模型,尤其会造成波速比的过高估计,从而影响对地壳结构和性质的解释(Lombardi et al., 2008; 查小惠等, 2013).

为了解决Moho面倾斜条件下的地壳厚度和波速比求取问题,前人通过系统分析认为,后方位角位于倾斜Moho面上倾方向的Ps波和多次波的走时与Moho面水平情况下的走时相对接近,因此提出只利用后方位角位于Moho面上倾方向的接收函数进行H-κ叠加,从而在一定程度上减小Moho面倾斜对准确求取地壳厚度和波速比的影响(Wang et al., 2010; 查小惠等, 2013).但是,这种方法只适应于Moho面倾角较小的情况,且实际资料处理中能够满足后方位角位于Moho面上倾方向的接收函数在数量上也往往非常有限.张洪双等(2009)在利用经典H-κ叠加方法思想的基础上,首先借助于研究区已有资料获得Moho面产状(倾向和倾角)信息,进而利用后方位角覆盖范围相对较全、信噪比较高的径向和切向接收函数资料来求取已知Moho面产状条件下的地壳波速比.然而,一方面,Moho面产状信息本身一般难以准确获取;另一方面,由于远震P波接收函数切向分量的信噪比往往较低,且受地震带分布本身的规律性以及观测时间的有限性等制约,实际资料中可以利用的远震事件后方位角覆盖范围往往是非常有限的,从而最终制约这一方法在实际资料处理中的应用和推广.因此,发展一种适应Moho面倾斜条件下的,可以同时估计地壳厚度-波速比-Moho面倾角(H-κ-θ)的实用方法,对于利用接收函数研究地壳组分和演化仍具有重要的现实意义.

为此,本文推导了Moho面倾斜条件下Ps转换波与后续多次波相对于直达P波的走时公式;基于经典H-κ叠加方法的思想,提出了一种可以同时估计地壳厚度-波速比-Moho面倾角的H-κ-θ叠加方法.通过理论模型测试,验证了该方法具备较高的稳定性和可靠性.最后,基于Hi-CLIMB台阵资料,在青藏高原南部Moho面倾斜地区进行了实际应用,显示出该方法具备较好的应用潜力.

1 基本原理 1.1 Moho面倾斜对经典H-κ叠加方法结果的影响

在Moho面倾斜的情况下,远震P波在Moho面发生转换,Ps转换波以及地表多次波的射线路径如图 1所示.

图 1 单个倾斜界面地震波射线路径图 Fig. 1 Seismic ray path on a single slant interface

为了研究Moho面倾斜条件下,远震P波接收函数的波形特征,以及Moho面倾角大小对经典H-κ叠加求取地壳厚度和波速比的影响,本文设计了Moho面倾斜的地壳模型,台站下方的Moho面深度为50 km,且向南倾斜(即界面倾向α=180°),地壳波速比为1.73,具体模型参数见表 1.

表 1 用于测试的地壳模型 Table 1 Crustal model for synthetic test

本文计算理论接收函数的过程,具体为给定Moho面倾角的情况下,基于表 1所示模型,利用RAYSUM程序(Frederiksen and Bostock, 2000)合成射线参数为0.06 s·km-1,后方位角为0°~360°的三分量地震记录;再利用时间域迭代反褶积算法(Ligorría and Ammon, 1999)计算得到给定Moho面倾角的理论接收函数.

对于Moho面倾斜条件下远震P波接收函数的波形特征,前人曾做过较为深入的研究(Langston, 1977; Owens et al., 1988; Cassidy, 1992; Savage, 1998; 陈九辉和刘启元,1999房立华和吴建平,2009查小惠等, 2013).本文基于表 1所示模型,分别计算了倾角为0°(水平界面)和10°(倾斜界面)情况下的远震P波理论接收函数以及直达P波振幅随后方位角的变化,如图 2所示.

图 2 远震P波理论接收函数以及直达P波振幅随后方位角的变化 (a) Moho面为水平界面的情况;(b) Moho面倾角为10°的情况.左侧为相应的地球模型;中间为径向和切向接收函数;右侧为直达P波振幅.红色波形所示对应振幅极性为正,蓝色波形所示对应振幅极性为负. Fig. 2 Synthetic P-wave receiver functions and the direct P-wave amplitude variations with back azimuths (a) Moho interface is flat; (b) Moho interface is slant with a dip angle of 10°. Left panels: earth models; Middle panels: synthetic radial- and tangential-component receiver functions; Right panels: amplitudes of the direct P-waves. The red-filled waveforms indicate the positive amplitudes, and the blue-filled waveforms indicate the negative amplitudes.

图 2可见,相对于Moho面水平的情况,当Moho面存在倾斜时,远震P波接收函数波形会出现以下特征:

(1) 切向分量出现了明显的能量;

(2) 直达P波的振幅随后方位角呈周期性变化.其中,径向分量振幅保持正极性不变,但振幅大小关于界面倾向对称;切向分量振幅极性发生反转,且振幅大小与极性均关于界面倾向反对称;

(3) 转换波Ps和多次波Pps、Pss的到时和振幅均随后方位角呈周期性变化.其中,径向分量振幅大小与极性均关于界面倾向对称;切向分量振幅极性发生反转,且振幅大小与极性均关于界面倾向反对称.

因此,在Moho面存在倾斜的情况下,远震P波接收函数中直达P波、转换波和后续多次波在到时、振幅乃至极性等方面,会出现随后方位角周期性系统变化的特征.总体上,在Moho面存在倾斜的情况下,后续多次波相对到时较水平界面到时提前,因此会出现界面深度偏小、波速比估计过高的现象,从而导致经典H-κ叠加方法不再适用于Moho面存在较大倾角的情况.

为了定量评估Moho面倾角对经典H-κ叠加方法求取地壳厚度和波速比的影响,我们基于表 1所示地壳模型,保持Moho面倾向(α=180°)不变,但在0°~15°之间由小到大、以1°间隔逐一增大Moho面倾角;利用经典H-κ叠加方法,求取相应的地壳厚度和波速比(或泊松比),同时采用bootstrap方法(Efron and Tibshirani, 1986)获得相应的统计误差.结果如图 3所示.

图 3 Moho面倾角为0°~15°情况下利用经典H-κ叠加方法获得的结果 (a)地壳厚度;(b)地壳波速比和泊松比.虚线所示为模型给定的真实值. Fig. 3 Results obtained by the H-κ stacking method for a slant Moho interface with dip angles of 0°~15° (a) Crustal thicknesses; (b) Bulk vP/vS ratios or Poisson′s ratios. The dashed line indicates the true value given by the model.

图 3可见,当Moho面存在倾斜时,利用经典H-κ叠加方法获得的结果与模型真实值之间存在一定程度偏差,且偏差随倾角的增大而增大.由于地壳厚度和波速比之间存在折衷,求得的地壳厚度相对真实值偏小,而波速比或泊松比则相对真实值偏大.泊松比σ是反映地壳组分特征的重要参数之一,其与波速比vP/vS之间存在换算关系:σ=0.5(1-((vP/vS)2-1)-1).通常认为,当σ≤0.26时,地壳总体组分偏酸性;当0.26<σ≤0.28时,地壳总体组分偏中性;当0.28<σ≤0.3时,地壳总体组分偏基性;σ>0.3时,壳内存在熔融或富含水等流体(嵇少丞等,2009).以此为例,本文测试模型给定地壳的平均波速比为1.73,对应泊松比为0.25,可认为相应的地壳总体组分偏酸性.但是,当Moho面倾角接近10°时,利用经典H-κ叠加方法获得的波速比已偏离真实值约0.04,对应泊松比大于0.26,相应的地壳总体组分却可解释为中性;当Moho面倾角接近15°时,估计的地壳厚度较真实值小15 km左右,波速比偏离真实值约0.18,对应的泊松比已大于0.3,则可能被解释为壳内存在熔融或者富集水等流体.由此可知,在Moho面倾角较大的情况下,利用经典H-κ叠加方法获得的地壳厚度和泊松比,均会大幅偏离真实值,会对后续解释带来严重误导.

1.2 Moho面倾斜条件下转换波和多次波理论到时的推导

考虑Moho面存在倾斜的情形,如图 4a所示,平面波从倾斜界面沿射线路径Γ到达地面接收点的走时t可以表示为:

(1)

图 4 (a) 单个倾斜界面地震波射线路径; (b) NEZ坐标系; (c)界面坐标系 Fig. 4 (a) Ray path of the seismic waves for a single slant interface; (b) NEZ coordinate system; (c) Interface coordinate system

其中,ξ为传播震相垂直慢度的绝对值.

于是,Ps转换波相对直达P波的走时可表示为:

(2)

其中,分别为在地壳中向上传播的P波和S波垂直慢度的绝对值.

PI为入射远震P波的慢度矢量,射线参数为p,后方位角为γ,倾斜界面下方P波速度为vP2,则在NEZ坐标系下,PI可以表示为:

(3)

PI从NEZ坐标系旋转到界面坐标系,坐标方向x1y1位于倾斜界面内,z1垂直于界面,旋转后可以得到:

(4)

其中α为倾斜界面倾向,θ为界面倾角.

令倾斜界面上方P波速度为vP1,应用Snell定律,则界面上方介质的P波慢度矢量可以表示为:

(5)

求该慢度矢量的垂直分量的绝对值,即为

(6)

倾斜界面的单位法向量n在NEZ坐标系下可表示为:

(7)

将(6)式化简,得到:

(8)

同理,可以得到

(9)

从而得到Ps转换波相对直达P波的走时:

(10)

由前文已知界面上方介质的P波慢度矢量P2I,将其从界面坐标系旋转到NEZ坐标系:

(11)

根据Snell定律,可以得到地表反射的P波与S波的慢度矢量PRSR分别为:

(12)

(13)

两者垂直分量的绝对值,即为在地壳中向下传播的P波垂直慢度的绝对值和S波垂直慢度的绝对值.

因此,多次波Pps和Pss相对直达P波的到时化简后可表示为:

(14)

(15)

(10)(14)(15)式即为适应倾斜Moho面的转换横波和地表多次波相对直达P波的理论到时公式.当界面水平时,即界面倾角θ=0°,上述三个公式即退化为经典H-κ叠加方法所利用的到时计算公式:

(16)

(17)

(18)

1.3 H-κ-θ叠加方法及模型测试

Zhu和Kanamori(2000)根据公式(16)—(18)提出了目前广泛用于求取地壳厚度和波速比的经典H-κ叠加方法,即

(19)

其中,r(t)为径向接收函数,t1t2t3分别对应指定地壳厚度h和波速比κ利用公式(16)—(18)计算得到的tPstPpstPssω1ω2ω3为加权系数,且Σωi=1.

在一定的变化范围内,按一定步长逐一给定hκ的取值,并利用公式(19)计算S(h, κ);当hκ的取值逼近真实值时,S(h, κ)的取值达到最大.反过来,即通过扫描搜索S(h, κ)取值最大时对应的h-κ,即可得到台站下方的地壳厚度和波速比,这就是经典H-κ叠加方法的基本原理.

比较公式(10)(14)(15)和(16)—(18)可知,在Moho面存在倾斜的情况下,转换波和多次波相对到时除主要受地壳厚度h和波速比κ影响外,还新增了Moho面的单位法向量n(取决于倾向α和倾角θ)、Moho面下方(即上地幔顶部)P波速度vP2等参数.为了定量测试新增参数αθvP2的影响,基于经典H-κ叠加方法的思想,构建新的目标函数如下:

(20)

这里t1t2t3分别对应利用公式(10)(14)(15)计算得到的tPstPpstPssx代表测试参数(即新增参数αθvP2中的任意一个,测试其中某个参数时,其他两个参数取真实值).与经典H-κ叠加方法的原理相同,通过扫描x, h, κ的取值,根据公式(20)求取x不同取值对应的S(h, κ, x)的最大值S(x);选取S(x)取值最大者对应的h-κ,即为求取的地壳厚度和波速比.

本文基于图 2b所示Moho面存在倾斜的情形,在0°≤α≤360°,0°≤θ≤20°,7.8 km·s-1vP2≤ 8.2 km·s-1范围内,分别测试不同参数取值对求取h, κ的影响,结果如图 5所示.

图 5 Moho面倾向、倾角和上地幔顶部P波速度对地壳厚度和波速比求取的影响 (a), (b)对应测试参数为Moho面倾向;(c), (d)对应测试参数为Moho面倾角;(e), (f)对应测试参数为上地幔顶部P波速度. Fig. 5 Impacts on estimation of the crustal thickness and bulk vP/vS ratios posed by dip directions, dip angles of the Moho interface and the P-wave velocities of the uppermost mantle (a) and (b) Test for the parameter of the dip directions of the slant Moho interface; (c) and (d) Test for the parameter of the dip angle of the slant Moho interface; (e) and (f) Test for the parameter of the P-wave velocities of the uppermost mantle.

图 5可知,当Moho面倾向α的误差≤30°时,对地壳厚度和速度比求取结果基本没有影响;当Moho面倾角θ的误差≤3°时,结果与真实值差别很小;上地幔顶部纵波速度vP2的误差对结果影响很小,完全可以忽略.因此,新增参数αθvP2之中,在Moho面倾向α的误差不大于30°的情况下,只有倾角θ对求取hκ影响较大.于是,在Moho面倾向α的误差不大于30°的情况下,公式(20)可进一步简化为:

(21)

也就是说,在Moho面倾向α的误差不大于30°的情况下,可以在一定取值范围内,通过扫描h, κ, θ的取值,根据公式(21)求取相应的S(h, κ, θ)的最大值S(θ);选取S(θ)取值最大者对应的h-κ-θ,即为台站下方的地壳厚度-波速比-Moho面倾角.这就是本文所提出的适应Moho面倾斜情况的H-κ-θ叠加方法的基本原理.

基于上述H-κ-θ叠加原理,我们计算了图 2b所示Moho面存在倾斜情况下的h-κ-θ.图 6给出了不同倾角θ(扫描范围为0°~20°,步长为1°)对应的S(h, κ, θ)的最大值S(θ).由图 6可见,当倾角θ的取值逼近真实值(~10°)时,由公式(21)计算得到的S(θ)也趋向于全局最大,意味着本文所提出的H-κ-θ叠加方法具有较高的稳定性,具备数据自驱动方面的适用性.当倾角θ的取值逼近真实值(~10°)时,对应的S(h, κ, θ)如图 7所示.由图 7可知,利用本文提出的H-κ-θ叠加方法所获得的能量团集中,误差小,可以准确求取h-κ-θ,而经典h-κ叠加方法却严重偏离真实值,进一步说明了本文提出的h-κ-θ方法的可靠性.

图 6 H-κ-θ叠加振幅最大值随倾角的变化 黑色圆点对应测试模型Moho面倾角的真实值10°,注意到此时叠加振幅达到全局最大. Fig. 6 Maximum of the stacking amplitude variations with the dip angles obtained by the H-κ-θ stacking method The black dot corresponds to the real dip angle (10°) of the Moho interface. Note that the stacking amplitude arrives at the global maximum when the testing dip angle closes to the real value.
图 7 本文提出的H-κ-θ叠加和经典H-κ叠加方法结果的比较 (a)不考虑Moho面倾斜,使用经典H-κ叠加获得的结果;(b)考虑Moho面倾斜,使用本文提出的H-κ-θ叠加获得的结果.虚线对应测试模型真实值(H=50 km,κ=1.73). Fig. 7 Comparison of the results from the improved H-κ-θ and the classical H-κ stacking methods (a) Classical H-κ stacking method without consideration of slant Moho; (b) H-κ-θ stacking method of this work considering slant Moho. The dashed line corresponds to the true values of the test model (H=50 km, κ=1.73).
2 实际资料处理

为了检验h-κ-θ叠加方法的实际应用效果,我们选择青藏高原Hi-CLIMB台阵资料进行了实际资料处理,剖面位置如图 8所示.现有的接收函数偏移成像结果显示(Nábělek et al., 2009; Wittlinger et al., 2009),沿Hi-CLIMB台阵剖面,喜马拉雅地区的Moho面存在大角度北倾特征.因此,该地区密集地震台阵资料为我们开展Moho面倾斜地区地壳厚度-波速比求取研究工作提供了较好条件.我们采用时间域迭代反褶积算法(Ligorría and Ammon, 1999)计算接收函数,所利用的有效事件震级大于6.0,震中距范围为30°~90°.

图 8 本文用到的Hi-CLIMB台站分布图 黑色三角形为宽频带地震台站,白色三角形为台站H0620,白色直线对应图 11所示剖面的投影位置. Fig. 8 Map showing distribution of Hi-CLIMB stations used in this study The black triangles indicate the broadband seismic stations, and white triangle indicates the station H0620. The white line is the projecting position of the profile shown in Fig. 11.

图 9图 10给出了台站H0620数据处理的实例.该台站(28.786°N,85.296°E;如图 8白色三角形所示)位于喜马拉雅块体内部.可利用的远震接收函数,以及直达P波振幅随后方位角的变化如图 9a所示.为了便于获得相对稳定可靠的、直达P波最大振幅随后方位角的变化规律,这里将直达P波最大振幅按后方位角每间隔20°进行算术平均.由图 9a可知,台站H0620的远震接收函数切向分量出现了明显的能量,且切向分量直达P波振幅在后方位角180°附近出现了极性反转,振幅大小和极性关于180°呈近似反对称;径向分量振幅保持极性不变,且振幅大小也关于180°呈近似对称,但180°附近对应出现振幅最大值.结合理论模拟测试结果所展示的规律性(如图 2所示),说明此处Moho界面的倾向应指向180°的反方向(即0°左右),大致为正北方向.因此,在利用本文发展的H-κ-θ叠加方法求取该台站下方地壳波速比的过程中,将Moho面倾向设为0°(可容许存在±30°误差),倾角的搜索范围为0°~20°,地壳平均P波速度设为6.4 km·s-1,上地幔顶部纵波速度设为8.3 km·s-1(Monsalve et al., 2008),结果如图 10b所示.

图 9 (a) H0620台站的接收函数和直达P波振幅随后方位角变化;(b)理论接收函数和直达P波振幅随后方位角变化 Fig. 9 (a) Receiver functions of station H0620 and the direct P-wave amplitude variations with the back azimuths; (b) Synthetic receiver functions and the direct P-wave amplitude variations with the back azimuths
图 10 H0620台站资料经典H-κ叠加和本文提出的H-κ-θ叠加方法获得结果的比较 (a)经典H-κ叠加方法得到的地壳厚度和平均波速比; (b)本文提出的H-κ-θ叠加方法获得结果. Fig. 10 Comparison of the results from the classical H-κ and the improved H-κ-θ stacking methods for the station H0620 (a) Crustal thickness and bulk vP/vS obtained by the classical H-κ stacking method; (b) Result obtained by the improved H-κ-θ stacking method.

H-κ-θ叠加方法获得的地壳厚度-波速比-Moho面倾角,以及给定的地壳平均P波速度,上地幔顶部纵波速度作为模型参数,计算理论三分量地震图并求取相应的理论接收函数,结果如图 9b所示.为便于对比,将理论接收函数的走时和振幅投射到相应的实际接收函数上,如图 9红色小圆点所示.通过比较可以看到,实际接收函数直达P波和后续多次波的走时和振幅与理论接收函数之间拟合程度较高,进一步说明利用本文发展的H-κ-θ叠加方法所获得的结果是可靠的.如果不考虑Moho面倾斜而使用经典H-κ叠加方法,H-κ域叠加能量团多而分散,最大能量对应的H-κ分别为43.1 km和1.99,远远偏离该地区地壳厚度和波速比,进一步说明经典H-κ叠加方法,已不再适用于该地区,即Moho面倾斜的影响已不可忽略.

基于以上流程,利用本文提出的H-κ-θ叠加方法,逐一处理可利用的每一个台站的接收函数,结果如图 11所示.通过与经典H-κ叠加方法获得的结果(引自Xu et al., 2013, 相应的vP/vS来自与此文通讯作者的私人通信)以及偏移成像获得的Moho面深度(Nábělek et al., 2009)之间的对比,可以看出利用H-κ-θ叠加方法得到的Moho面深度与接收函数偏移成像获得的Moho面形态之间吻合得较好,变化也相对稳定;而经典H-κ叠加方法获得的结果,变化剧烈,部分区域求得的波速比过高.因此利用本文提出的H-κ-θ叠加方法,能够得到较为合理的结果,显示出该方法具备较好的应用潜力.

图 11 不同方法获得的Hi-CLIMB剖面下方Moho面深度和波速比分布 为方便直接对比起见,此图中Moho面深度均参照Nábělek等(2009)统一校正为相对海平面的深度;剖面位置如图 8白色直线所示. Fig. 11 Moho depths and bulk vP/vS ratios beneath Hi-CLIMB profile obtained by different methods Note that for the direct comparison, Moho depths in this map are relative to sea level (Nábělek et al., 2009). The profile position is shown by the white line in Fig. 8.
3 结论

经典H-κ叠加方法隐含着Moho面是水平界面的假设条件,忽略了Moho面产状的影响.在Moho面存在较大倾角的情况下,利用经典H-κ叠加方法得到的地壳厚度和速度比会在较大范围偏离实际模型.本文推导了Moho面倾斜条件下Ps转换波与后续多次波相对于直达P波的理论到时公式,基于经典H-κ叠加方法的思想,提出了可以同步求取地壳厚度-波速比-Moho面倾角的H-κ-θ叠加方法.通过理论模型检验,验证了该方法具备较高的稳定性和可靠性,并将此方法应用于青藏高原南部Hi-CLIMB台阵资料,显示出较好的应用潜力.

需要说明的是,本文仅侧重于对当前流行的H-κ叠加方法在适应Moho面倾斜情况方面的改进.实际上,已有大量研究表明,除Moho面倾角外,下地壳各向异性的存在也会对H-κ叠加方法的结果产生重要影响(查小惠等, 2013; 任骏声和沈旭章, 2015; Shen et al., 2015).因此,为了准确求取台站下方的地壳结构和属性参数,需全面分析远震接收函数的径向和切向分量波形随后方位角变化的规律性,以便进一步区分和评估Moho面倾斜和地壳各向异性的影响,从而最终给出具有针对性的解决方案.

致谢  感谢Hi-CLIMB研究团队和IRIS网站提供了地震波形数据,Frederiksen等提供了合成理论地震图程序,以及宋晓东老师提供了Hi-CLIMB测线各流动台站利用经典H-κ叠加方法获得的地壳厚度和波速比结果.
References
Cassidy J F. 1992. Numerical experiments in broadband receiver function analysis. Bulletin of the Seismological Society of America, 82(3): 1453-1474.
Chen J H, Liu Q Y. 1999. Lithospheric receiver function in 3-D laterally inhomogeneous media using Maslov asymptotic theory. Chinese Journal of Geophysics (in Chinese), 42(1): 84-93.
Chen Y, Zhang Z J, Sun C Q, et al. 2013. Crustal anisotropy from Moho converted Ps wave splitting analysis and geodynamic implications beneath the eastern margin of Tibet and surrounding regions. Gondwana Research, 24(3-4): 946-957. DOI:10.1016/j.gr.2012.04.003
Chen Y, Xu Y G, Xu T, et al. 2015. Magmatic underplating and crustal growth in the Emeishan Large Igneous Province, SW China, revealed by a passive seismic experiment. Earth and Planetary Science Letters, 432: 103-114. DOI:10.1016/j.epsl.2015.09.048
Christensen N I. 1996. Poisson's ratio and crustal seismology. Journal of Geophysical Research:Solid Earth, 101(B2): 3139-3156. DOI:10.1029/95jb03446
Efron B, Tibshirani R. 1986. Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy. Statistical Science, 1(1): 54-75. DOI:10.1214/ss/1177013815
Fang L H, Wu J P. 2009. Effects of dipping boundaries and anisotropic media on receiver functions. Progress in Geophysics (in Chinese), 24(1): 42-50.
Frederiksen A W, Bostock M G. 2000. Modelling teleseismic waves in dipping anisotropic structures. Geophysical Journal International, 141(2): 401-412. DOI:10.1046/j.1365-246x.2000.00090.x
Ji S C, Wang Q, Yang W C. 2009. Correlation between crustal thickness and Poisson's ratio in the North China Craton and its implication for lithospheric thinning. Acta Geologica Sinica (in Chinese), 83(3): 324-330.
Langston C A. 1977. The effect of planar dipping structure on source and receiver responses for constant ray parameter. Bulletin of the Seismological Society of America, 67(4): 1029-1050.
Ligorría J P, Ammon C J. 1999. Iterative deconvolution and receiver-function estimation. Bulletin of the Seismological Society of America, 89(5): 1395-1400.
Lombardi D, Braunmiller J, Kissling E, et al. 2008. Moho depth and Poisson's ratio in the Western-Central Alps from receiver functions. Geophysical Journal International, 173(1): 249-264. DOI:10.1111/j.1365-246X.2007.03706.x
Monsalve G, Sheehan A, Rowe C, et al. 2008. Seismic structure of the crust and the upper mantle beneath the Himalayas:Evidence for eclogitization of lower crustal rocks in the Indian Plate. Journal of Geophysical Research:Solid Earth, 113(B8): B08315. DOI:10.1029/2007jb005424
Nábělek J, Hetényi G, Vergne J, et al. 2009. Underplating in the Himalaya-Tibet collision zone revealed by the Hi-CLIMB experiment. Science, 325(5946): 1371-1374. DOI:10.1126/science.1167719
Owens T J, Crosson R S, Hendrickson M A. 1988. Constraints on the subduction geometry beneath western Washington from broadband teleseismic waveform modeling. Bulletin of the Seismological Society of America, 78(3): 1319-1334.
Ren J S, Shen X Z. 2015. The neighborhood inversion of the crust anisotropic structures with receiver functions. Progress in Geophysics (in Chinese), 30(5): 2043-2054. DOI:10.6038/pg20150507
Savage M K. 1998. Lower crustal anisotropy or dipping boundaries? Effects on receiver functions and a case study in New Zealand. Journal of Geophysical Research:Solid Earth, 103(B7): 15069-15087. DOI:10.1029/98jb00795
Shen X Z, Yuan X H, Ren J S. 2015. Anisotropic low-velocity lower crust beneath the northeastern margin of Tibetan Plateau:Evidence for crustal channel flow. Geochemistry, Geophysics, Geosystems, 16(12): 4223-4236. DOI:10.1002/2015GC005952
Wang P, Wang L S, Mi N, et al. 2010. Crustal thickness and average Vp/Vs ratio variations in southwest Yunnan, China, from teleseismic receiver functions. Journal of Geophysical Research:Solid Earth, 115: B11308. DOI:10.1029/2009JB006651
Wittlinger G, Farra V, Hetényi G, et al. 2009. Seismic velocities in Southern Tibet lower crust:A receiver function approach for eclogite detection. Geophysical Journal International, 177(3): 1037-1049. DOI:10.1111/j.1365-246X.2008.04084.x
Xu Z J, Song X D, Zhu L P. 2013. Crustal and uppermost mantle S velocity structure under Hi-CLIMB seismic array in central Tibetan Plateau from joint inversion of surface wave dispersion and receiver function data. Tectonophysics, 584: 209-220. DOI:10.1016/j.tecto.2012.08.024
Zha X H, Sun C Q, Li C. 2013. The effects of dipping interface and anisotropic layer on the result of H-κ method. Progress in Geophysics (in Chinese), 28(1): 121-131. DOI:10.6038/pg20130113
Zhang H S, Tian X B, Teng J W. 2009. Estimation of crustal Vp/Vs with dipping Moho from receiver functions. Chinese Journal of Geophysics (in Chinese), 52(5): 1243-1252. DOI:10.3969/j.issn.0001-5733.2009.05.013
Zhu L P, Kanamori H. 2000. Moho depth variation in southern California from teleseismic receiver functions. Journal of Geophysical Research:Solid Earth, 105(B2): 2969-2980. DOI:10.1029/1999jb900322
陈九辉, 刘启元. 1999. 合成三维横向非均匀介质远震体波接收函数的Maslov方法. 地球物理学报, 42(1): 84-93.
房立华, 吴建平. 2009. 倾斜界面和各向异性介质对接收函数的影响. 地球物理学进展, 24(1): 42-50.
嵇少丞, 王茜, 杨文采. 2009. 华北克拉通泊松比与地壳厚度的关系及其大地构造意义. 地质学报, 83(3): 324-330.
任骏声, 沈旭章. 2015. 基于接收函数和邻域算法的地壳各向异性结构反演. 地球物理学进展, 30(5): 2043-2054. DOI:10.6038/pg20150507
查小惠, 孙长青, 李聪. 2013. 倾斜界面和各向异性地层对H-κ搜索结果的影响. 地球物理学进展, 28(1): 121-131. DOI:10.6038/pg20130113
张洪双, 田小波, 滕吉文. 2009. 接收函数方法估计Moho倾斜地区的地壳速度比. 地球物理学报, 52(5): 1243-1252. DOI:10.3969/j.issn.0001-5733.2009.05.013