复杂地壳接收函数H-κ叠加——以安纳托利亚板块为例

危自根，储日升*，陈凌，崇加军，李志伟

1 中国科学院测量与地球物理研究所大地测量与地球动力学国家重点实验室,武汉 4300772 中国科学院地质与地球物理研究所岩石圈演化国家重点实验室,北京 1000293 中国科学院青藏高原地球科学卓越创新中心,北京 100101



复杂地壳接收函数H-κ叠加
——以安纳托利亚板块为例

危自根1，储日升1*，陈凌2,3，崇加军1，李志伟1

1 中国科学院测量与地球物理研究所大地测量与地球动力学国家重点实验室,武汉 4300772 中国科学院地质与地球物理研究所岩石圈演化国家重点实验室,北京 1000293 中国科学院青藏高原地球科学卓越创新中心,北京 100101

本文理论分析了具有不同沉积层和壳幔过渡带结构的接收函数及其相关的H-κ叠加结果，然后采用接收函数H-κ叠加和波形反演方法获得了具有复杂构造演化历史的中北安纳托利亚板块的地壳厚度(H)、VP/VS(κ)和VS结构.理论分析表明：厚的沉积层或沉积层和厚的壳幔过渡带共存都会使H-κ叠加失效；渐变型壳幔过渡带导致H-κ叠加的H位于过渡带中间，且随着频率增大逐渐靠近过渡带上方；倒转型壳幔过渡带导致H-κ叠加具有多极值，其结果可能反应过渡带内最大波阻抗界面上的地壳结构；1 km·s-1的VP变化会导致H-κ叠加的H变化7 km，而κ变化较小.实际资料分析表明：中北安纳托利亚H，κ和VS具有强烈的横向不均匀性，大部分区域沉积层厚度<0.5 km，局部地区壳幔过渡带厚度>3 km；北安纳托利亚断层切穿地壳，在局部地区可能存在流体；研究区存在残留古老的小陆块体.本文研究表明，仔细分析接收函数波形和其随方位角的变化特征且用其他地震学方法进行约束，有助于采用H-κ叠加方法获取复杂地壳结构信息.

接收函数H-κ叠加；波形反演；沉积层；壳幔过渡带；安纳托利亚板块

1 引言

接收函数由于其对间断面的特殊敏感性被广泛用来获取地球内部结构信息.接收函数H-κ叠加方法(Zhu and Kanamori,2000，以下简称H-κ方法)通过联合利用Moho面转换波和多次反射波，能同时反演地壳厚度(H)和平均VP/VS(κ)，进而约束区域地壳成分和构造演化，获得了广泛的应用.然而，H-κ方法假设地壳为水平单层均匀结构模型，这与广泛存在倾斜Moho界面、各向异性地层和厚沉积层等复杂的实际地壳结构不相符.真实地壳的不均匀性特征往往造成不同震中距和方位角接收函数的转换波和多次波震相的到时存在差异(Peng and Humphreys,1997;Zheng et al.,2005;Yeck et al.,2013)，进而影响H-κ叠加的可靠性.因此，理论与实际分析复杂地壳结构对接收函数H-κ叠加结果的影响十分必要.

目前，已有不少学者通过理论合成复杂地壳结构的接收函数，分析其对H-κ叠加的影响.罗艳等(2008)分析了具有逐渐递增速度的沉积层地壳的理论接收函数和H-κ叠加结果，建立了首到波峰与直达P波到时差与沉积层厚度的定量关系，认为H-κ方法不能有效获得存在巨厚沉积层的地壳结构.Lombardi等(2008)和Wang等(2010)通过分析倾斜Moho界面的接收函数H-κ叠加，发现上倾方向的结果更接近于地壳真实值.房立华和吴建平(2009)通过研究不同地壳模型的接收函数认为，分析不同方位角的径向和切向接收函数直达P波和Ps转换波振幅与到时随方位角的变化规律，可以较好地区分倾斜界面和介质各向异性的影响.查小惠等(2013)通过对不同倾斜Moho界面和各向异性地层的接收函数分析，认为H-κ叠加结果的偏离程度和Moho倾斜角度成正比，水平轴各向异性地层使所有方位角的H估计值偏大，与快轴夹角较小方位角的κ值严重偏小.个别学者通过改进H-κ方法以提高其适用性和结果精度：比如，N次方根H-κ叠加搜索(Niu et al.,2007)，两层或三层地壳的H-κ方法(Tang et al.,2008；Yeck et al.,2013)等.上述研究为采用H-κ方法获取存在倾斜Moho、各向异性等地壳的H和κ提供了参考.然而，目前对地壳中广泛存在的沉积层，厚的壳幔过渡带结构对接收函数H-κ叠加的影响还缺乏系统的对比分析，因此有必要做进一步研究.

除了复杂地壳结构等客观因素外，接收函数挑选和反演参数设置等人为因素，也会影响H-κ叠加结果的可靠性.因此，系统分析人为因素对H-κ叠加的影响也十分必要.安纳托利亚板块是由阿拉伯板块和欧亚板块于新生代碰撞汇聚形成的中特提斯造山带的一部分，正进行着多块体碰撞拼贴造山过程(Berberian and King,1981；Pearce et al.,1990)，地壳结构十分复杂.为了研究北安纳托利亚断层附近地壳结构特征，科研人员在该区域布设了由39个地震台站组成的台阵(YL)，且已经采用H-κ方法获得了台阵下方的地壳厚度与VP/VS(Vanacore et al.,2013).可公开使用的地震波形数据(IRIS网站)和已发表的H-κ叠加结果，使北安纳托利亚地区成为研究复杂地壳客观和主观因素影响H-κ叠加的理想场所.

鉴于前人已经系统理论分析了倾斜Moho界面和各向异性层对接收函数H-κ叠加的影响(房立华和吴建平，2009；查小惠等，2013)，本文将主要理论分析不同沉积层和壳幔过渡带结构对H-κ叠加结果的影响.除此外，我们采用和Vanacore等(2013)相同参数的H-κ方法获得了YL台阵下方的H和κ，并分析了两者结果的异同及其差异的来源；采用不同频率H-κ方法分析了研究区沉积层和壳幔过渡带的结构特征；采用多频接收函数波形反演方法获得了台阵下方的S波速度结构.本文试图通过理论和实际资料的接收函数H-κ叠加分析，为采用该方法获取具有沉积层或者壳幔过渡带等复杂结构的地壳厚度和VP/VS提供参考，并探讨安纳托利亚板块地壳的结构特征和演化进程.

2 沉积层和壳幔过渡带对接收函数H-κ叠加响应

2.1 接收函数H-κ叠加方法介绍

接收函数主要由转换波和多次反射波组成，记录了台站下方速度间断面对地震波的响应.在远震接收函数波形记录中，一般初至P波之后30 s内主要包含Moho面的Ps转换波和PpPs及PpSs+PsPs地表多次反射波.联合这几个震相可以约束地壳厚度与平均P波和S波速度比信息.接收函数H-κ叠加方法假设在一维水平均匀地壳模型下，构造一个H-κ平面内的叠加函数S(H,κ)，

S(H,κ)=WPsr(tPs)+WPpPsr(tPpPs)

-WPpSs+PsPsr(tPpSs+PsPs),

(1)

其中，r(t)为径向接收函数，WPs，WPpPs，WPpSs+PsPs为满足三者之和为1的不同震相的加权系数.最佳的地壳厚度与波速比对应叠加函数S(H,κ)的最大能量值.上述方法即为接收函数的H-κ叠加方法(Zhu and Kanamori,2000;Wei et al.,2015;危自根和陈凌，2012)，其主要优点是不用人工挑选震相走时，从而可以避免挑选震相时人为因素的影响.

接收函数H-κ叠加过程中，需要给定地壳P波速度和不同震相的加权系数等信息.我们采用反射率法(Langston,1979;Levin and Park,1997)计算水平单层模型下的合成地震图，然后由合成的地震事件进一步计算理论接收函数(图1a，后面合成理论接收函数方法相同)，并分析不同参数和接收函数分布的H-κ叠加结果与理论模型的差异(图1b).结果表明，随着射线参数值增大，Ps转换震相到时逐渐增加，而PpPs和PpSs+PsPs震相到时逐渐减小；输入P波速度每增加0.1 km·s-1，H大约增加0.7 km，κ减小不超过0.004；不同震相加权值组合(0.6/0.3/0.1,0.5/0.5/0.0,0.4/0.3/0.3)得到的H一样，κ变化小于0.01；不同接收函数分布得到(0.04～0.05 s·km-1，0.07～0.08 s·km-1，0.04～0.08 s·km-1)的H变化小于0.5 km，κ变化小于0.02.上述分析表明，在单层水平均匀地壳模型下，接收函数H-κ叠加获得的VP/VS对输入参数和接收函数分布不敏感，而地壳厚度则对输入的P波速度敏感.

图1 H为40 km，VP为6.4 km·s-1，κ为1.78，高斯系数为2.0的不同射线参数的理论接收函数(a)和不同输入参数和接收函数分布的H-κ叠加结果(b).图a黑色实线标注了Moho面典型震相的理论到时Fig.1 Synthetic receiver functions for crustal model with H of 40 km,VP of 6.4 km·s-1,κ of 1.78 and a gauss value of 2.0 (a) and corresponding H-κ stacking results for different given parameters and distribution of receiver functions (b).The black lines in (a) denote the theoretical arrival times of typical phases from the Moho

2.2 沉积层对接收函数H-κ叠加影响

沉积层广泛分布在地球大陆和海洋板块之中.由于沉积层相对较低的速度和密度，在其底界面形成的转换波或多次反射波可能会干扰来自Moho面的转换震相，从而造成接收函数H-κ叠加结果与真实值存在偏差(Zelt and Ellis,1999).已有学者提出一些方法来减弱接收函数中沉积层带来的影响，比如波形拟合(Zheng et al.,2005；Wei et al.,2016)，利用沉积层先验信息进行波场延拓(Langston,2011)，对接收函数进行特定的带通或者带阻滤波等(Leahy et al.,2012).对于H-κ方法，Yeck等(2013)提出了连续两层H-κ叠加技术.该方法首先利用高频H-κ叠加获得沉积层结构信息，然后采用低频和改进的H-κ公式求取沉积层底界面到Moho面的厚度和VP/VS，在美国Powder River和Deven盆地Moho面Ps震相尚未完全被沉积层多次波震相干扰的台站中得到了有效的应用.Yu等(2015)提出利用接收函数自相关构建的共振消除滤波器来消除沉积层多次波的干扰，进而获得沉积层和其基底到Moho的厚度与VP/VS，在薄于5 km的沉积层情况下效果较好.上述对接收函数的预处理以及H-κ方法的改进有效地减弱了薄的沉积层对Moho面带来的影响，但是对厚的沉积层或者Moho面Ps转换震相受到沉积层多次波震相干扰时效果较差.

我们计算了具有0～8 km厚的不同类型沉积层地壳的理论接收函数，并分析其对H-κ叠加结果的影响(图2).结果表明，单层(图2a,2b)和渐变型(图2c,2d)沉积层存在时，接收函数波形变复杂，Moho面Ps和多次转换波与直达P波的到时差随着沉积层厚度的增厚逐渐增大，H-κ叠加结果与模型值存在差异.在渐变型沉积层情况下，当厚度在0～8 km变化时，H-κ叠加获得的H和模型值差异不超过1 km，κ和模型值差异不超过0.03(图2e)，表明其对H-κ叠加结果影响较小.这种小的差异与其Moho面Ps震相没受到沉积层震相干扰一致(图2d).在单层沉积层情况下，当厚度不超过2 km时，直达P波受到沉积层Ps震相干扰导致最大振幅偏离0时刻；当厚度不超过4 km时，H-κ叠加的H和模型值差异不超过1 km，κ和模型值差异不超过0.03(图2e)，这与其Moho面Ps震相没受到沉积层震相干扰一致；当沉积层为6 km或者8 km厚时，H-κ叠加的H和模型值差异超过1 km，κ和模型值差异超过0.06(图2e)，这与其Moho面Ps震相受到沉积层多次波震相干扰有关(图2b).上述观测表明，单层沉积层比渐变型沉积层对接收函数H-κ叠加影响更大，且巨厚的沉积层会使H-κ叠加结果严重偏离地壳的真实值.

2.3 壳幔过渡带对接收函数H-κ叠加影响

壳幔过渡带是壳幔物质与能量交换的动力边界，对地壳和地幔的形成与演化以及深层动力过程有着重要作用.已有研究表明，壳幔过渡带存在三种基本类型(Davydova et al.,1972；滕吉文，2006)：一级间断面型(尖锐型)，连续或不连续的速度梯度带型(渐变型)，和高、低速相间的薄层束型(倒转型).这三种不同的壳幔过渡带在世界大陆地区都存在，比如壳幔过渡带比较尖锐的中国燕山带(赵俊猛等，1999)，5～10 km厚渐变型壳幔过渡带的乌克兰地盾地区(Pavlenkova,1988)，倒转型壳幔过渡带的美国犹他山脉和东部盆地地区(Braile,1977).

复杂的壳幔过渡带会干扰接收函数来自Moho面震相的到时与波形信息，进而使H-κ叠加结果偏离真实值.已有学者通过分析具有渐变型壳幔过渡带的理论接收函数H-κ叠加发现，不同频率获得的地壳厚度存在差异(Rumpfhuber et al.,2009).基于上述因素，我们分析了不同性质壳幔过渡带在不同频率下的理论接收函数和相应的H-κ叠加结果(图3).理论接收函数表明，壳幔过渡带主要影响Moho面PpPs多次波波列宽度和振幅信息，且随着壳幔过渡带厚度的增加，多次波波列更加复杂(图3(b,c,d,e)).对具有不同类型壳幔过渡带的接收函数进行H-κ叠加发现，渐变型壳幔过渡带极值范围明显比倒转型范围大且得到的H几乎都位于过渡带内，倒转型壳幔过渡带接收函数H-κ叠加呈现多极值特征.当渐变型壳幔过渡带厚度为4 km时，H-κ叠加获得的H随着频率增高逐渐靠近过渡带顶层，κ和模型值差异从～0.6%(G=1.0)增加到～2.6%(G=4.0)；当厚度为8 km时，H-κ叠加获得的H随着频率增高逐渐靠近过渡带上方，κ和模型值差异不超过1%(图3f).倒转型壳幔过渡带不同频率H-κ叠加获得的H接近过渡带底层，与模型值最大差异～3.5%，κ和模型值差异不超过1%(图3f).通过对比地壳模型发现(图3a,3f)，倒转型壳幔过渡带H与κ可能对应着过渡带内具有最大波阻抗界面上的地壳结构特征.

图2 采用高斯系数为2.0，射线参数为0.065 s·km-1计算的具有不同沉积结构地壳模型(a,c)的理论接收函数(b,d)和其H-κ叠加结果与模型值的差异(e).H-κ叠加给定的初始P波速度和模型值相同.(b,d)中短棒标明Moho震相(Ps，PpPs，PpSs+PsPs)和沉积层震相(s.Ps，s.PpPs)的理论到时Fig.2 Synthetic receiver functions (b,d) for crustal models with different sediment structure (a,c) with a gauss value of 2.0 and ray parameter of 0.065 s·km-1 and the differences between the H-κ stacking results and the real values.The given VP in the H-κ stacking is the same to the value of crustal model.Short bars in (b,d) show the theoretical arrivals of Moho phases (Ps,PpPs,PpSs+PsPs) and sediment phases (s.Ps,s.PpPs)

图3 采用不同高斯系数(G)和射线参数为0.065 s/km计算的具有不同壳幔过渡带类型地壳模型(a)的理论接收函数(b,c,d,e)和其H-κ叠加结果与模型值的差异(f).H-κ叠加给定的初始P波速度和模型值相同.b,c,d,e中短棒标明尖锐型壳幔界面Moho震相的到达时间Fig.3 Synthetic receiver functions (b,c,d,e) for crustal models (a) with different crust-mantle transition with different gauss values and ray parameter of 0.065 s/km and the differences between the H-κ stacking results and the real values (f).The given VP in the H-κ stacking is the same to the value of crustal model.Short bars in (b,c,d,e) show the theoretical arrivals of Moho phases from sharp crust-mantle transition

图4 采用高斯系数为2.0和射线参数为0.065 s·km-1计算的沉积层和壳幔过渡带共存的地壳模型(a)的理论接收函数(b)和其H-κ叠加结果与模型值的差异(c).H-κ叠加给定的初始P波速度和模型值相同.b中短棒标明尖锐型壳幔界面Moho震相的到达时间Fig.4 Synthetic receiver functions (b) for crustal models (a) with coexistence of sediment layer and crust-mantle transition with a gauss value of 2.0 and ray parameter of 0.065 s·km-1 and the differences between the H-κ stacking results and the real values (c).The given VP in the H-κ stacking is the same to the value of crustal model.Short bars in (b) show the theoretical arrivals of Moho phases from sharp crust-mantle transition

图5 安纳托利亚板块流动台阵(蓝色正方形)和本研究所用地震事件分布(右上角插图黑点).红点表示板块内1905到2014年大于5级的地震，数据来自美国地质调查局.带白色三角形的黑线标明目前正活跃的俯冲消减带，带黑色三角形的黑线标明块体缝合边界(Okay and Tüysüz,1999).黑色箭头表示块体相对欧亚板块运动的方向Fig.5 Distribution of temporary array (blue squares) in the Anatolian plate and teleseismic events (solid dots in right inset) used in the study.Red dots show earthquakes with magnitudes >5 from 1905 to 2014 from USGS.Black lines with white triangles show the active subduction zones currently.Black lines with black triangles show the suture boundaries of different blocks (Okay and Tüysüz,1999).The black arrows show the plate movement directions relative to the Eurasian Plate

图6 YL台阵单台所有径向接收函数叠加波形和接收函数数量分布(a)以及典型台站的H-κ叠加示意图(b,c).图a中细线和粗线分别为本文和Vanacore等(2013)使用的单台接收函数数量.十字叉(b,c)表明最佳地壳厚度和波速比值Fig.6 Stacked receiver function waveforms and number of receiver functions for each station in YL array (a) and receiver functions sorted by back azimuths and H-κ stacking results for representative stations.Thin and thick lines (a) show numbers of selected receiver functions for each station from this study and Vanacore et al.(2013),respectively.Crosses (b,c) represent the best estimates of crustal thickness and VP/VS ratio

2.4 沉积层和壳幔过渡带共存对接收函数H-κ叠加影响

第2.2和2.3节分析表明，2 km和4 km厚的沉积层对接收函数H-κ叠加结果影响相对较小，渐变型壳幔过渡带H-κ叠加的地壳厚度位于壳幔过渡带内而VP/VS和真实值相差不大，倒转型壳幔过渡带的H-κ叠加结果与模型值差异也不大.在实际地壳结构中，沉积层和壳幔过渡带可能会同时存在，进而对接收函数H-κ叠加结果产生影响.基于上述分析，我们测试了沉积层(2 km和4 km)和渐变以及倒转型过渡带(6 km)同时存在的接收函数H-κ叠加响应(图4).理论接收函数表明，Moho面转换波和多次反射波都受到沉积层震相或者壳幔过渡带震相干扰，波形非常复杂.接收函数H-κ叠加结果表明，在渐变型壳幔过渡带和沉积层共存情况下，获得的H和模型值的差别高于10%，κ和模型值差异大于6%；在倒转型壳幔过渡带和沉积层共存情况下，获得的H和模型值的差别高于9%，κ和模型值差异大于4%，且随着沉积层厚度增加差别越大.上述结果表明，接收函数H-κ叠加方法不能有效获得厚的沉积层和壳幔过渡带共存下的地壳结构信息.

3 中北安纳托利亚板块YL台阵接收函数H-κ和波形反演研究

3.1 安纳托利亚板块构造背景与YL台阵介绍

安纳托利亚板块是由多个小陆块在新生代碰撞拼贴而成(Berberian and King,1981；Pearce et al.,1990;Taymaz et al.,2007)，构造和动力学演化特征十分复杂(图5).东安纳托利亚板块与青藏高原具有类似的构造特征和岩浆作用，目前以南北向缩短和东西向伸展为主；中安纳托利亚板块主要受到东西向缩短和南北向拉伸作用；西安纳托利亚板块主要受到南北向拉伸作用(常承法，1986;Taymaz et al.,2007).多期次碰撞俯冲使安纳托利亚板块具有多个缝合带和复杂的地壳结构，不仅形成了长达1500 km的近东西走向的右旋安纳托利亚走滑断层，还形成了多个大小不等的沉积盆地(Dewey et al.,1986;Taymaz et al.,2007).北安纳托利亚断层自形成后5个百万年来断层积累位移达到85 km，地震频发，在1973年到2000年发生过10个6级以上的地震(Taymaz et al.,2007).已有的地质地球物理研究表明安纳托利亚板块存在广泛的可能形成厚壳幔过渡带的基性-酸性的岩浆(Eyuboglu et al.,2012)，存在局部的各向异性(Peng and Ben-Zion,2004)和小尺度强烈的Moho面深度变化(Vanacore et al.,2013).为了研究北安纳托利亚断层附近地壳结构特征，科研人员在该区域布设了流动地震台阵YL(图5)，运行时间从2005到2008年.目前，已有学者用接收函数H-κ叠加方法获得了该台阵35个台站下方的地壳厚度和VP/VS(Vanacore et al.,2013)以及地震层和三维P波速度结构(Yolsal-Çevikbilen et al.,2012)，发现该区域地壳厚度与VP/VS横向变化强烈，地震层厚度约为15 km且大多数地震发生在支链上.

3.2 YL台阵接收函数H-κ叠加分析

YL台阵记录了大量高质量的远震信息，为研究北安纳托利亚板块的地壳结构提供了数据基础.为了和Vanacore等(2013)结果进行对比，与之相同的参数被用来计算接收函数和进行H-κ叠加.我们选取5.5～7.0震级和30°～90°震中距的远震事件，截取P波到时前20 s和后100 s波形，采用1 Hz的截止频率和0.2的水平因子，使用最大熵谱反褶积方法计算接收函数(吴庆举和曾融生，1998).为确保用高质量的接收函数进行H-κ叠加，我们首先删除信噪比小于3和Ps震相不清晰或受到明显干扰的接收函数，然后将单台单个接收函数与所有接收函数叠加后的波形做互相关，删除相关系数小于0.6的波形，最后我们再一次人工挑选接收函数.对于挑选后的接收函数，我们选取VP为6.2 km·s-1，Ps，PpPs和PpSs/PsPs震相加权值分别为0.5，0.25和0.25来进行接收函数H-κ叠加搜索.图6展示了本文最终用到的28个台站的径向接收函数叠加后的波形，以及本文和Vanacore等(2013)接收函数数量分布以及两个典型台站的H-κ叠加结果.从图6a可以看出，几乎所有台站Moho面转换震相都较为清晰，但多次波却不是非常明显，部分台站直达P波和Ps转换波之间有震相，个别台站(DOGL)直达P波受到干扰，表明研究区域地壳结构强烈的横向不均匀性.本文和Vanacore等(2013)使用的单台接收函数数量具有大致的变化趋势，在个别台站数量几乎一致(比如YESI)，而在某些台站数量差别较大(比如YIKI).

图7展示了YL台阵下方的地壳厚度与VP/VS以及本文与Vanacore等(2013)结果的差异.十字叉标明Vanacore等(2013)展示了H-κ叠加结果，而本文却没获得有效的地壳厚度与VP/VS台站位置.这些台站主要分布在具有5～8 km厚的Cank盆地(Taymaz et al.,2007)，接收函数Moho面Ps震相不清晰或者存在不可分辨的双极值.本文研究表明，研究区地壳厚度与VP/VS分别在30 km到42 km和1.65到2.0之间大尺度波动；北安纳托利亚断层两侧地壳结构变化剧烈，比如台站KGAC和PELI，地壳厚度和VP/VS变化分别达6 km和0.2；断层附近VP/VS跳跃大，最小值1.70，比如BEDI，最大值接近2.0，比如INCE.本文和Vanacore等(2013)大的地壳结构差异主要分布在北安纳托利亚断层附近，比如台站INCE，地壳厚度和VP/VS差分别达到5 km和超过0.1.研究区地壳结构的强烈横向变化与其是由多个小陆块和增生楔在新生代拼合而成有关.在断层局部区域，高的VP/VS可能与其下方存在破碎岩石、地幔物质底侵或者流体有关；低的VP/VS可能与其下地壳的减薄从而导致石英质岩石占主要部分有关.理论接收函数H-κ叠加分析表明，倾斜Moho(Lombardi et al.,2008;Wang et al.,2010)和各项异性(房立华和吴建平，2009；查小惠等，2013)等复杂的地壳结构会使H-κ叠加结果随着不同的方位角存在差异，Moho面倾角超过15°以上下倾方向的结果会严重偏离地壳真实值.我们分析了典型台站INCE的接收函数特征和H-κ叠加结果来探讨造成本文和Vanacore等(2013)结果差异的原因(图6).对于该台站，接收函数波形没有显示出倾斜Moho和各向异性具有的周期性和转换震相到时差的特征(Wang et al.,2010)，Vanacore等(2013)和本文用于进行H-κ叠加的接收函数数量分别为110和67个，两者之比为1.65.该台站接收函数Moho面PpPs多次波震相一致性较差，且随着反方位角变化明显.不同反方位角范围内(0～360°，0～120°，120～200°，200～360°)的H-κ叠加结果表明，地壳厚度和VP/VS的最大差别分别超过6 km和达到0.15，表明台站下方地壳结构小尺度的横向变化.上述分析表明，不同的接收函数数量和分布特征可能是造成本文和Vanacore等(2013)H-κ叠加结果差异大的原因.

Yeck等(2013)通过理论分析和实际数据观测表明，高频接收函数H-κ叠加可以获得沉积层的厚度与VP/VS信息.选取高斯系数为10和VP为2.2 km·s-1(参考全球地壳模型Crust 1.0，Laske et al.,2013)，我们对YL台阵进行了接收函数H-κ叠加.结果表明，除了DOGL台站获得的沉积层厚度为0.9 km之外(图7e)，其余台站得到的沉积层厚度都小于0.5 km(图7f).对于DOGL台站，0.9 km的沉积层厚度与全球地壳模型Crust 1.0是一致的，沉积层不同震相的理论到时大致对应了最大的能量值，且在G=2.0情况下的直达P震相(图6a)显示了与存在1 km厚沉积层的理论接收函数波形相似的延迟现象特征(图2b).而对于ALIC等台站，H-κ叠加结果的沉积层震相都集中在P波后，其获得的<0.5 km厚的沉积层有两种可能性.第一，这些台站沉积层厚度确实非常薄，与其具有高的海拔与短时间的沉积一致，是本文支持的观点.第二，这些沉积层H-κ叠加结果反映的是一种假象.台站ALIC在2.5 s左右有一个一致的震相且在2 km和4 km附近有一些弱的能量集中区域，代表真的沉积层结构信息.上述分析表明，在采用简单的高频H-κ叠加获取沉积层信息时，需要先验沉积层信息来约束H-κ搜索范围.理论分析表明(图3)，存在厚壳幔过渡带的地壳的接收函数H-κ叠加结果随着频率变化而变化.Wei等(2011)对比鄂尔多斯南侧流动剖面的H-κ叠加结果和波形反演(Zheng et al.,2009)得到的壳幔过渡带S波速度结构发现，H-κ叠加得到的厚度位于壳幔过渡带速度最大变化处.本文计算了高斯系数为1.0，2.0，3.0和4.0的接收函数H-κ叠加结果(图7f)，来试图约束研究区的壳幔过渡带结构.结果表明，在高斯系数为2.0,3.0和4.0情况下，除台站INSU和KGAC之外，其他台站地壳厚度差异几乎不超过1 km；绝大部分台站(ALIN，KAVA，KGAC，KUZA除外)VP/VS差都不超过0.05.除台站KAVA，ALIN，BEDI和KUZA外(～3 km)，高斯系数为1.0得到的地壳厚度与其他频率的差别普遍不超过2 km.除台站TEPE，KUZA和KAVA外，高斯系数为1.0得到的VP/VS与其他频率的差别普遍在0.1以内.不同频率的H-κ叠加获得的地壳厚度和VP/VS差异表明，研究区部分台站下方的壳幔过渡带不是尖锐型的，存在一定的厚度，这与前人获得的该地区7.9～8.0 km·s-1的Pn波速度(Kuleli et al.,2004)和大量的基性-酸性的岩浆分布一致(Eyuboglu et al.,2012)，表明其强烈的壳幔相互作用.

图7 YL台阵下方地壳厚度(a)和平均VP/VS(c)，本文与Vanacore等(2013)结果差(b,d)，高斯系数为10的典型台站沉积层的H-κ叠加(e,f)和高斯系数为1，2，3，4的H-κ叠加结果(g).灰线(a,b,c,d)为北安纳托利亚断层.红线(e,f)为H-κ叠加得到的沉积层不同震相的理论到时Fig.7 The distributions of crustal thickness (a) and average VP/VS ratio (c) and the differences of H-κ stacking results between this study and Vanacore et al.(2013) (b,d) beneath the YL array,the H-κ stacking for sediment lay beneath typical stations with a gauss factor of 10 (e,f) and the H-κ stacking resutls for different gauss factor of 1.0,2.0,3.0 and 4.0 (g),respectively.Grey lines (a,b,c,d) show the North Anatolian fault.The red lines (e,d) show the theoretical arrival times of different phases

3.3 YL台阵接收函数波形反演分析

中北安纳托利亚区域构造复杂，研究其地壳速度分布能为探讨其构造演化提供有用的信息.上述接收函数H-κ叠加方法获得的地壳厚度、平均VP/VS、沉积层和壳幔过渡带特征为采用波形反演方法反演地壳的速度结构提供了约束.作为对照，接收函数波形反演获得的地壳速度模型能用来验证H-κ叠加结果的可靠性.Li等(2010)发展了接收函数差异演化非线性全局优化反演方法，能较好反演4～5层地壳S波速度结构，被Wei等(2016)用来联合H-κ方法共同约束四川盆地的地壳结构.Li等(2016)进一步发展了该方法，通过同时拟合不同频率的接收函数来约束地壳不同尺度的S波速度结构.利用H-κ方法获得的地壳厚度、VP/VS和研究区P波层析成像结果(Vanacore et al.,2013)作为约束，本文通过同时拟合高斯系数为1.0,2.0,3.0和4.0的接收函数反演了YL台阵下方地壳的S波速度结构.图8展示了两个典型台站的S波速度反演结果.从图中可以看出，在经过20000次的运算之后，目标函数收敛较好(图8e,8k)，不同频率段的合成接收函数与实际资料拟合较好，尤其是Moho面的主要震相.值得注意的是，台站DOGL反演结果显示存在1.4 km厚、S波速度仅为1.71 m·s-1的沉积层，与高频沉积层H-κ叠加的结果一致.

图9展示了研究区28个台站下方的S波速度结构分布.图9b,9c,9d分别为4～8 km，14～18 km和24～28 km内平均S波速度(采用GMT软件中surface命令，网格化数据间隔为6′).图9e,9f为两条大致与北安纳托利亚断层主链平行和垂直的剖面下方地壳的S波速度.结果表明，研究区地壳不同深度都存在强烈的小尺度的S波速度变化.在中下地壳(图9c，9d)，研究区断层中部附近(台站TEPE,ARSL,ISKE,DUMA,KARA等)显示出相对高速的特征，速度主要在3.9～4.0 km·s-1范围内变化；在断层右下方区域(KUZA,BAGB,CAKM)则显示出相对较低的速度特征，速度主要分布范围为3.2～3.7 km·s-1.平行断层走向的剖面(图9e)显示，台站KGAC和KARA下方存在从20 km直到Moho面的高速体异常，而台站KUZA,BAGB,CAKM下方中下地壳速度明显偏低，高速体和低速体速度差异达到10%；垂直断层走向的剖面(图9f)下方大部分区域在20 km深度附近存在高速异常，最高速度达到3.8 km·s-1，这种异常在台站DUMA下方直接延伸到上地幔顶部，而其周边台站CALT和DERE中下地壳速度明显要低.地质学研究认为安纳托利亚板块是随着新特提斯洋闭合，由多个古老小陆块拼合而成的(Berberian and King,1981；Pearce et al.,1990;Taymaz et al.,2007).区域P波层析成像(Yolsal-Çevikbilen et al.,2012)发现研究区地壳存在局部的P波高速物质，并推测这些异常体可能代表古老陆块的残留.本文观测到的研究区局部的中下地壳高速体和Yolsal-Çevikbilen等(2012)P波高速体异常位置一致，支持我们观测到的S波高速异常区域可能是古老陆块残留体的推测.

3.4 复杂地壳结构接收函数H-κ叠加分析

上述理论计算和实际资料分析表明，在存在沉积层和厚的壳幔过渡带等复杂结构地区，采用接收函数H-κ叠加方法求取地壳厚度与VP/VS时，需要非常谨慎.在做H-κ叠加之前，需要仔细分析接收函数不同震相波形的特征，比如直达P波的延迟性，其与Moho面Ps震相之间的复杂性，多次波震相波列结构和频谱特征，并基于已有的地质和地球物理资料大致判断台站下方是否具有沉积层或者厚的过渡带等，删除Moho面Ps震相受到严重干扰且转换波和多次波一致性较差的接收函数.除此外，还需分析接收函数随方位角或者震中距的变化特征，大致判断台站下方Moho面的分布状态.在进行H-κ叠加时，利用已有的地质和地球物理结果作为先验资料，联合其他地震学方法，约束地壳厚度与平均VP/VS，删除具有不可分辨的多极值的台站，对接收函数随方位角变化剧烈的台站分不同的方位角范围进行叠加等是十分有必要的.

图8 典型台站多频接收函数波形反演示意图图(a,b,c,d,g,h,i,j)中红色为拟合的接收函数，灰色为实际数据.(e,k)为反演资料与实际数据的拟合程度.(f,l)为反演的速度模型，其中黑色实线为初始模型，虚线为最终反演结果.Fig.8 Schematic diagrams for multiple frequency waveform inversion of receiver functions for representative stationsGrey and red lines (a,b,c,d,g,h,i,j) denote the real and recovered receiver functions,respectively.Panels e and k show objective function decreasing with iteration increasing.Panels f and l show the recovered S-wave velocity models,where black lines show the initial models and dashed lines show the last models.

图9 YL台阵下方S波速度结构图(b,c,d)为不同深度范围的平均S波速度;(e,f)为a中两条典型剖面下方的S波速度分布.十字叉代表台站位置.Fig.9 S-wave velocity distribution beneath the YL array Panels b,c and d show the average S-wave velocity in different depth range in the study region. Panels e and f show the crustal S-velocity along two typical profiles in a.Crosses show the station location.

4 结论

本文理论计算了具有不同沉积层和壳幔过渡带结构地壳的接收函数，并进一步分析其对H-κ叠加的影响.结果表明：沉积层的存在会增加接收函数Moho面转换波和多次反射震相的到时且有可能干扰Moho面的Ps波震相；壳幔过渡带的存在主要影响Moho面多次波震相的到时和波形；>6 km厚的沉积层会导致接收函数H-κ叠加的地壳厚度与VP/VS误差都大于3%；渐变型壳幔过渡带会导致H-κ叠加的地壳厚度位于壳幔过渡带之间，并且随着频率增高逐渐靠近过渡带上方；倒转型壳幔过渡带会导致H-κ叠加结果出现多极值现象，其结果可能对应着过渡带内具有最大波阻抗界面上的地壳结构特征；沉积层和厚的壳幔过渡带共存时会导致H-κ叠加方法失效.

我们采用接收函数H-κ叠加和波形反演方法获得了中北安纳托利亚板块YL台阵下方的地壳厚度、VP/VS和S波速度结构信息.接收函数H-κ叠加不同参数理论分析表明1 km·s-1的VP变化会产生7 km的地壳厚度差异，而对VP/VS的影响较小，不同加权系数对H-κ叠加结果影响较小.本文与Vanacore等(2013)H-κ结果的差异主要来自不同的接收函数数量以及分布特征.中北安纳托利亚地区地壳厚度、VP/VS和S波速度具有强烈的横向不均匀性，大部分台站上方的沉积层厚度比较薄，少数台站壳幔过渡带相对较厚，与研究区年轻的构造演化和沉积历史以及强烈的壳幔相互作用一致.沿断层附近区域部分台站下方波速比超过1.9，暗示可能存在地幔物质底侵或者局部流体.断层两侧大的地壳厚度差，表明其可能是穿透地壳的大型断层.北安纳托利亚地区局部中下地壳的高速异常体可能代表着古老陆块的残留.

致谢 感谢IRIS数据中心提供YL台阵地震事件波形数据，感谢两位匿名审稿专家的建设性意见.

Berberian M,King G C P.1981.Towards apaleogeography and tectonic evolution of Iran.Can.J.Earth Sci.,18(2):210-265.

Braile L W.1977.Interpretation of crustal velocity gradients and Q structure using amplitude-corrected seismic refraction profiles.∥Heacock J G,Keller G V,Oliver J E,et al eds.The Earth′s Crust.

Washington D C:American Geophysical Union,20:427-439.

Chang C F.1986.Tectonic comparison between Tibet plateaupost-collision and Turkey-Iranian plateau (II).Translated Seismology Geology (in Chinese),8(3):1-8.

Davydova N I,Kosminskaya I P,Kapustian N K,et al.1972.Models of the Earth′s crust and M-boundary.Z.Geophys.,38(3):369-393.

Dewey J F,Hempton M R,Kidd W S F,et al.1986.Shortening of continental lithosphere:The neotectonics of Eastern Anatolia—a young collision zone.∥Coward M P,Rtes A C eds.Collision Tectonics.Geological Society,London,SpecialPublications,19(1):1-36.

Eyuboglu Y,Santosh M,Yi K,et al.2012.Discovery of Miocene adakitic dacite from the Eastern Pontides Belt (NE Turkey) and a revised geodynamic model for the late Cenozoic evolution of the Eastern Mediterranean region.Lithos,146-147:218-232.

Fang L H,Wu J P.2009.Effects of dipping boundaries and anisotropic media on receiver functions.Progress in Geophysics (in Chinese),21(4):42-50.

Kuleli S,Toksoz M,Gurbuz C,et al.2014.Crustal structure study in Turkey with controlled seismic sources.∥AGU Fall Meeting 2004.Abstracts,S13B-1058.

Langston C A.1979.Structure under Mount Rainier,Washington,inferred from teleseismic body waves.J.Geophys.Res.,84(B9):4749-4762.

Langston C A.2011.Wave-field continuation and decomposition for passive seismic imaging under deep unconsolidated sediments.Bull.Seism.Soc.Amer.,101(5):2176-2190.

Laske G,Masters G,Ma Z T,et al.2013.Update on CRUST1.0-A1-degree global model of Earth′s crust.∥EGU General Assembly 2013.Abstracts,EGU2013-2658.

Leahy G M,Saltzer R L,Schmedes J.2012.Imaging the shallow crust with teleseismic receiver functions.Geophys.J.Int.,191(2):627-636.

Levin V,Park J.1997.Crustal anisotropy in the Ural Mountains Foredeepfrom teleseismic receiver functions.Geophysical Research Letters,24(11):1283-1286.

Li X L,Li Z W,Hao T Y,et al.2016.Multi-frequency receiver function inversion with a global optimizing approach.Computers and Geosciences,under review.

Li Z W,Hao T Y,Xu Y,et al.2010.A global optimizing approach for waveform inversion of receiver functions.Computers &Geosciences,36(7):871-880.

Lombardi D,Braunmiller J,Kissling E,et al.2008.Moho depth and Poisson′s ratio in the Western-Central Alps from receiver functions.Geophys.J.Int.,173(1):249-264.

Luo Y,Chong J J,Ni S D,et al.2008.Moho depth and sedimentary thickness in Capital region.Chinese J.Geophys.(in Chinese),51(4):1135-1145,doi:10.3321/j.issn:0001-5733.2008.04.022.

Niu F L,Bravo T,Pavlis G,et al.2007.Receiver function study of the crustal structure of the southeastern Caribbean plate boundary and Venezuela.J.Geophys.Res.,112(B11):B11308.doi:10.1029/2006JB004802.

Okay A I,Tüysüz O.1999.Tethyan sutures of northern Turkey.∥Durand B,Jolivet L,Horvárt,et al eds.The Mediterranean Basins:Tertiary Extension within the Alpine Orogen.Geological Society,London,Special Publications,156(1):475-515.

Pavlenkova N I.1988.The nature of seismic boundaries in the continental lithosphere.Tectonophysics,154(3-4):221-225.

Pearce J A,Bender J F,De Long S E,1990.Genesis of collision volcanism in Eastern Anatolia,Turkey.Journal of Volcanology and Geothermal Research,44(1-2):189-229.

Peng X H,Humphreys E D.1997.Moho dip and crustal anisotropy in Northwestern Nevada from teleseismic receiver functions.Bull.Seism.Soc.Amer.,87(3):745-754.

Peng Z G,Ben-Zion Y.2004.Systematic analysis of crustal anisotropy along the Karadere-Düzce branch of the North Anatolian fault.Geophys.J.Int.,159(1):253-274.

Rumpfhuber E M,Keller G R,Sandvol E,et al.2009.Rocky Mountain evolution:Tying continental dynamics of the Rocky Mountains and Deep Probe seismic experiments with receiver functions.J.Geophys.Res.,114(B8):B08301,doi:10.1029/2008JB005726.

Tang C C,Chen C H,Teng T L.2008.Receiver functions for three-layer media.Pure Appl.Geophys.,165(7):1249-1262.

Taymaz T,Yllmaz Y,Dilek Y.2007.The geodynamics of the Aegean and Anatolia:Introduction.∥Taymaz T,Yllmaz Y,Dilek Yeds.The geodynamics of the Aegean and Anatolia.Geological Society,London,Special Publications,291(1):1-16.

Teng J W.2006.Research on layer-bundle fine structures and physical attributes of crust-mantle boundary in deep earth.Journal of Jilin University (Earth Science Edition) (in Chinese),36(1):1-23.

Vanacore E A,Taymaz T,Saygin E.2013.Moho structure of the Anatolian Plate from receiver function analysis.Geophys.J.Int.,193(1):329-337.

Wang P,Wang L S,Mi N,et al.2010.Crustal thickness and average VP/VSratio variations in southwest Yunnan,China,from teleseismic receiver functions.J.Geophys.Res.,115(B11):B11308,doi:10.1029/2009JB006651.

Wei Z G,Chen L,Xu W W.2011.Crustal thickness and VP/VSratio of the central and western North China Craton and its tectonic implications.Geophys.J.Int.,186(2):385-389.

Wei Z G,Chen L.2012.Regional differences in crustal structure beneath northeastern China and northern North China Craton:Constraints from crustal thickness and VP/VSratio.Chinese J.Geophys.(in Chinese),55 (11):3601-3614,doi:10.6038/j.issn.0001-5733.2012.11.009.

Wei Z G,Chu R S,Chen L.2015.Regional differences in crustal structure of the North China Craton from receiver functions.Science China Earth Sciences,58(12):2200-2210.

Wei Z G,Chen L,Li Z W,et al.2016.Regional variation in Moho depth and Poisson′s ratio beneath eastern China and its tectonic implications.Journal of Asian Earth Sciences,115:308-320.

Wu Q J,Zeng R S.1998.The crustal structure of Qinghai-Xizang plateau inferred from broadband teleseismic waveform.Chinese J.Geophys.(in Chinese),41(5):669-679.

Yeck W L,Sheehan A F,Schulte-Pelkum V.2013.Sequential H-κstacking to obtain accurate crustal thicknesses beneath sedimentary basins.Bull.Seism.Soc.Amer.,103(3):2142-2150.

Yolsal-Çevikbilen S,Biryol C B,Beck S,et al.2012.3-D crustal structure along the North Anatolian Fault Zone in north-central Anatolia revealed by local earthquake tomography.Geophys.J.Int.,188(3):819-849.

Yu Y Q,Song J G,Liu K H,et al.2015.Determining crustal structure beneath seismic stations overlying a low-velocity sedimentary layer using receiver functions.J.Geophys.Res.,120(5):3208-3218,doi:10.1002/2014JB011610.

Zelt B C,Ellis R M.1999.Receiver-function studies in the Trans-Hudson Orogen,Saskatchewan.Can.J.Earth Sci.,36(4):585-603.

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 Geophys.(in Chinese),28(1):121-131,doi:10.6038/pg20130113.

Zhao J M,Zhang X K,Zhao G Z,et al.1999.Structure of crust-mantle transitional zone in different tectonic environments.Earth Science Frontiers (in Chinese),6(3):165-172.

Zheng T Y,Zhao L,Chen L.2005.A detailed receiver function image of the sedimentary structure in the Bohai Bay Basin.Physics of the Earth and Planetary Interiors,152(3):129-143.

Zheng T Y,Zhao L,Zhu R X.2009.New evidence from seismic imaging for subduction during assembly of the North China Craton.Geology,37(5):395-398.

Zhu L P,Kanamori H.2000.Moho depth variation in southern California from teleseismic receiver functions.J.Geophys.Res.,105(B2):2969-2980.

附中文参考文献

常承法.1986.青藏高原碰撞期后构造与土耳其—伊朗高原的比较(II).地震地质译丛,8(3):1-8.

房立华,吴建平.2009.倾斜界面和各向异性介质对接收函数的影响.地球物理学进展,24(1):42-50.

罗艳,崇加军,倪四道等.2008.首都圈地区莫霍面起伏及沉积层厚度.地球物理学报,51(4):1135-1145,doi:10.3321/j.issn:0001-5733.2008.04.022.

滕吉文.2006.地球深部壳-幔边界的层束精细结构与物理属性研究.吉林大学学报(地球科学版),36(1):1-23.

危自根,陈凌.2012.东北地区至华北北缘地壳结构的区域差异:地壳厚度与波速比的联合约束.地球物理学报,55(11):3601-3614,doi:10.6038/j.issn.0001-5733.2012.11.009.

吴庆举,曾融生.1998.用宽频带远震接收函数研究青藏高原的地壳结构.地球物理学报,41(5):669-679.

查小惠,孙长青,李聪.2013.倾斜界面和各向异性地层对H-κ搜索结果的影响.地球物理学进展,28(1):121-131,doi:10.6038/pg20130113.

赵俊猛,张先康,赵国泽等.1999.不同构造环境下的壳-幔过渡带结构.地学前缘,6(3):165-172.

(本文编辑 胡素芳)

Analysis of H-κ stacking of receiver functions beneath crust with complex structure:taking the Anatolia Plate as an example

WEI Zi-Gen1,CHU Ri-Sheng1*,CHEN Ling2,3,CHONG Jia-Jun1,LI Zhi-Wei1

1 State Key Laboratory of Geodesy and Earth′s Dynamics,Institute of Geodesy and Geophysics,Wuhan 430077,China2 State Key Laboratory of Lithospheric Evolution,Institute of Geology and Geophysics, Chinese Academy of Sciences,Beijing 100029,China3 Chinese Academy of Sciences Center for Excellence in Tibetan Plateau Earth Sciences,Beijing 100101,China

We calculated the synthetic receiver functions and their H-κstacking results for crust with different sediments and crust-mantle transition zones,and then inverted the crustal thickness (H) and VP/VSratio (κ) and S-wave velocity beneath the Anatolia plate with complex tectonic evolution history by H-κstacking and waveform inversion of receiver functions.The synthetic tests show:H-κstacking method is invalid in cases of existence of thick sediment or coexistence of sediment and crust-mantle transition zone;H obtained by H-κstacking is located in the transition zone and close to the top of transition zone with the increase of frequency in gradual crust-mantle transition zone;multiple extreme values exist in the H-κstacking and the last result is probably related to the crustal structure above the biggest wave impedance interface for inverted crust-mantle transition zone;the change of 1 km·s-1of P velocity results in 7 km variation of H but has little impact onκ.The real data analysis show:prominent lateral inhomogeneity in H and k and S-wave velocity were observed beneath the north Anatolian plate,with sediment of <0.5 km thick in most of the region and >3 km thick crust-mantle transition zone in local area;the north Anatolia fault extends to the mantle with existence of fluid in local area;residual old blocks may exist in the crust of the study region.Based on this study,we propose that the careful analysis of the waveform and changing trend along the azimuth of the receiver functions and the combination with other seismic methods can help to obtain reliable structural information for crust with complex structure by the H-κstacking method.

H-κstacking of receiver functions;Waveform inversion;Sediment;Crust-mantle transition zone;Anatolia plate

危自根，储日升，陈凌等.2016.复杂地壳接收函数H-κ叠加——以安纳托利亚板块为例.地球物理学报，59(11):4048-4062,

10.6038/cjg20161110.

Wei Z G,Chu R S,Chen L,et al.2016.Analysis of H-κstacking of receiver functions beneath crust with complex structure:taking the Anatolia Plate as an example.Chinese J.Geophys.(in Chinese),59(11):4048-4062,doi:10.6038/cjg20161110.

国家重点基础研究发展计划项目(2013CB733203)，国家自然科学基金项目(41225016,41604056)和大地测量与地球动力学国家重点实验室自主课题(Y509169216)资助.

危自根，男，副研究员，主要从事壳幔结构研究.E-mail:weizigen@whigg.ac.cn

*通讯作者 储日升，男，研究员，主要从事地震学研究.Email:chur@asch.whigg.ac.cn

10.6038/cjg20161110

P315

2016-06-08，2016-08-29收修定稿