基于程函方程与三维速度模型的中国东北地区地幔过渡带接收函数研究

张炎， 钮凤林， 宁杰远,4*

1 北京大学地球与空间科学学院， 北京 100871 2 美国莱斯大学地球、环境与行星科学系， 休斯顿 77005 3 中国石油大学石油资源与勘探国家重点实验室和非常规天然气研究所， 北京 102249 4 河北红山巨厚沉积与地震灾害国家野外科学观测研究站， 北京 100871

0 引言

中国东北地区位于欧亚板块东缘，中朝块体和西伯利亚地台之间，西北部与中亚造山带相连，西南部与华北克拉通相接，东临日本海.受构造运动及俯冲板块的影响，该地区地质历史复杂，经历了蒙古—鄂霍次克洋的张开与闭合，古亚洲洋的关闭，阿穆尔板块与西伯利亚板块的碰撞，太平洋板块俯冲等地质历史事件，以盆地-山脉相间的形式为主，形成了三个主要的组成部分：中部的松辽盆地，西部的大兴安岭以及东南部的长白山山脉(杜晓娟等，2009；Zhou and Wilde，2013).早白垩时期，东北地区经历了大规模的地壳拉张运动，开始形成大量的断陷盆地，并伴随火山活动.晚白垩纪以来，该地区发育了大量的火山活动(Liu et al.，2001；Chen et al.，2007).中国东北地区的新生代火山岩主要是碱性玄武岩，火山活动以零星偶发为主(Fan and Hooper，1991；Liu et al.，2001)，广泛分布于盆地边缘与山脉中，包括长白山火山群、五大连池火山群、阿尔山火山群和阿巴嘎火山群等(刘嘉麒，1999)(图1a).研究该区域的地幔结构，可为该地区上下地幔之间的动力学关系，俯冲板块的现状和火山起源等科学问题的研究提供重要线索.

对于该地区新生代火山活动，普遍认为与太平洋板块的俯冲有一定的关系，但对板内火山活动起源的详细机制，仍然存在着争议.通过不同的观测资料和研究方法，前人给出了该地区新生代火山形成原因的几种观点.全球和区域层析成像结果显示(如Fukao et al.，1992；Bijwaard et al.，1998；Zhao，2004；Huang and Zhao，2006)，中国东北地区下方，地幔过渡带中存在近水平延伸约1500 km的地震波高速异常体，被解释为在地幔过渡带中滞留的俯冲板块；很多研究者(如Lei and Zhao, 2005，2006; Huang and Zhao, 2006；Zhao et al.，2007，2009；Duan et al.，2009；Li and Van Der Hilst，2010；Wei et al.，2012，2015；Lei et al.，2013)认为长白山火山与滞留在地幔过渡带内的太平洋板块脱水造成的物质在“大地幔楔”中上涌有关；Tang等(2014)则认为太平洋滞留板块存在“空缺”， 来自下地幔的热物质沿着“空缺”上涌是长白山火山的深部动力学机制，该观点也受到一些研究结果的支持(Liu et al.，2015；Zhang et al.，2016)；另外，Zou等(2008)认为俯冲板块在地幔过渡带内的停滞堆积，可能驱动软流圈上涌至浅部导致减压熔融；Zhang等(2016)认为五大连池和阿尔山火山的活动可能是地幔岩石圈的拆沉造成软流圈物质上涌所引起的，同时指出在阿巴嘎火山群附近可能存在小型地幔柱.

地幔内部在全球尺度存在两个地震波速度和密度发生快速跃变的间断面，其全球平均深度约为410 km和660 km，因此被称为410-km和660-km间断面.两个间断面间的地幔被称为是上下地幔的地幔过渡层或地幔过渡带.通常认为410-km间断面是由地幔中矿物橄榄石到瓦兹利石，660-km间断面是由林伍德石到布里奇曼石与镁方铁石的高温高压相变所致(Katsura and Ito，1989；Helffrich and Wood，2001；周春银等，2010).由于与410-km和660-km两个地幔间断面相关的相变分别具有正的和负的克拉伯龙斜率(Ito and Takahashi，1989)，在温度低于周围地幔的区域，410-km间断面深度将会变浅而660-km间断面深度则会变深，形成增厚的地幔过渡带；在温度高于周围地幔的区域，情况则正好相反.由于地幔物质相变对温度的依赖性，可以通过测量这两个间断面深度起伏及地幔过渡带厚度(410-km和660-km间断面的深度差)的横向变化，约束地幔过渡带内的低温(如板块俯冲区域)与高温(如热物质上涌区域)结构的分布.

地震学探测是认识地球深部结构、物质组成和动力学过程的重要手段.对于地球内部间断面结构的研究，主要依靠间断面处体波的折射、反射以及震相转换等信息，其中，接收函数是获取间断面深度最有效的方法之一.前人在中国东北地区的接收函数研究工作中(如Ai et al.，2003, 2008；Li and Yuan，2003；Liu et al.，2015；Tian et al.，2016；Zhang et al.，2016；杨凡等，2021)，由于使用的方法和参考速度模型不同，导致各项研究获得的间断面深度与过渡带厚度结果存在一定差异.例如Li和Yuan(2003)与Ai 等(2003)发现在长白山火山下方660-km间断面有20～30 km的下陷，而在其西侧下陷程度明显减小；Liu等(2015)得到了相似的结果，同时在东北地区北部得到了大范围660-km间断面略微下陷区和长白山火山群以西660-km间断面局部抬升的成像结果；Zhang等(2016)同样发现了长白山火山下方和西侧660-km间断面下陷和局部抬升的结构特征，但在北部没有发现大范围的660-km间断面下陷迹象.在进行410-km和660-km间断面接收函数成像时，通常采用共转换点(Common Conversion Point，CCP)叠加方法(Dueker and Sheehan，1997；Zhu，2000；Gilbert et al.，2003).为了获得高质量的地震图像，则需要一个相对准确的三维速度参考模型，采用一种既高效又准确的方法计算地震波在三维模型中的走时.目前，大多数研究采用一维射线路径走时校正算法，即首先利用一个一维参考速度模型进行射线追踪，计算一维射线路径，然后计算三维速度模型在一维路径中走时扰动的总和来获取近似的三维走时(Ai et al.，2008； Liu et al., 2015； Zhang et al., 2016).在后文中，我们将简称该近似三维走时计算方法为射线追踪校正法.当参考模型中存在较强速度异常时，三维射线路径的影响就会变得十分明显，从而导致射线追踪校正法计算走时的准确度大为降低.同时，射线追踪校正法需要对每个接收函数进行射线追踪与走时校正计算，计算效率较低.

本文采用了一种利用快速行进法(Fast Marching Method，FM法)求解程函方程的方法，以提高计算三维模型中地震波传播时间的准确度和效率(Rawlinson and Sambridge，2004a，2004b；De Kool et al.，2006；Guan and Niu，2018).FM法是采用有限差分法求解程函方程，追踪走时变化(Sethian and Popovici，1999).Guan和Niu(2018)优化了利用FM法计算P-S转换波相对走时的参数设置，并通过对一维与三维合成模型的数值模拟计算，验证了该方法的高效性.该方法考虑了三维路径效应的影响，提高了走时计算的准确度；在计算过程中先求出台站下方目标区域内P波与S波走时场，再根据台站与转换点的地理位置进行插值运算，求取P-S转换波相对走时，从而提高了走时计算效率.Cheng等(2016)在Kirchhoff叠前偏移成像中，使用了FM法计算转换波相对走时，取得了较好的结果.该方法横向分辨率高，对复杂地下结构成像较好，但是受台站分布与计算效率较低的影响，不太适用于台站分布不均的较大尺度区域研究.相比之下，本文所用的基于FM法的CCP叠加方法兼顾了计算效率与准确性，更适合于大尺度区域地幔过渡带研究.

本文中，我们首先利用国际合作项目NECESSArray以及中国国家地震台网记录到的实际波形数据，对该方法在大尺度研究区域中计算的可行性进行了评估，然后利用Tao等(2018)全波形反演获得的东亚地区高分辨率三维速度模型作为参考模型，对中国东北地区地幔过渡带结构进行了接收函数成像，最后依据本文获取的成像结果，探讨了太平洋板块在地幔过渡带的存在状态，并对板内火山的起源进行了尝试性解释.

1 数据与方法

我们采用了Liu等(2015)的接收函数数据，该数据是通过国际合作在中国东北地区布设的临时地震仪台站组成的大型地震台阵(NECESSArray)和中国国家地震台网，共255个台站观测的远震记录计算得到(图1a)，远震地震事件数为788个(图1b)，挑选出接收函数总计45505条.台站较为均匀地分布在东经115°—135°，北纬40°—49°的研究区域内，观测时间为2009年9月至2011年8月.三维速度参考模型采用的是Tao等(2018)通过三维全波形反演方法获得的东亚地区P波和S波速度模型(FWEA18).

1.1 一维模型相对走时计算

CCP叠加成像方法将转换震相振幅按照射线路径反向投影到其真实的空间位置，利用叠加提高信噪比来获得间断面的空间分布特征.首先通过选取适当的参考速度模型计算不同深度的P-S转换波相对走时，再将观测到的时间序列中各个采样点转换到相应的转换深度，从而得到转换震相与转换点位置的对应关系，完成时深转换.接收函数上各个点对应的振幅可以视为在这个转换位置上生成的转换波振幅.最后将相同转换点附近的接收函数对应的转换震相的振幅进行叠加，就可以得到该转换点的转换波振幅.

参考速度模型的选取会对CCP叠加结果产生较大的影响.参考速度模型更接近真实的地幔，计算得到的转换波相对走时则更准确，成像结果也更可靠.对于一条接收函数，将从震源到台站的直达P波的走时记作tP0，P-S转换点的深度记作d，从震源到该转换点的P波和从该转换点到台站的S波的走时分别记作tP和tS，则该P-S转换波相对直达P波的相对走时ΔtPds可以表示如下：

ΔtPds=tP+tS-tP0.

(1)

本研究采用了De Kool等(2006)通过快速行进法计算走时的程序(程序包名称为fm3d)，该程序采用了高效且无条件稳定的快速行进算法数值求解球面程函方程(Sethian，1996；Popovici and Sethian，2002).以一维IASP91速度模型为例，对于任一地震事件E到任一台站S的ΔtPds计算过程如下.

首先，选择一个范围合适的计算区域，并选择合适的水平向与垂直向间距划分出走时场的三维网格点.对于地震事件E，运行fm3d代码，计算远震P波的走时场，并将三维网格中每个网格点对应的走时储存起来(图2a).直达P波走时tP0则可以由地表离台站S位置最近的4个网格点储存的走时，通过插值得到.对于转换深度d，通过一维射线追踪得到该转换点的地理位置，在P波走时场内找到对应转换点附近8个网格点储存的走时信息，插值得到tP(图2b).然后，通过相同的计算过程，可以求出S波的走时场，并将从震源到台站的S波走时记为tS0，从震源到转换深度d的S波走时记为tSd.通过(tS0-tSd)可以计算得到约等于公式(1)里的tS的值(图2b).

图2 FM法计算相对走时示意图(a)、(b)与(c)、(d)分别表示在一维与三维速度模型中计算转换波相对走时的过程. 红色五角星表示计算区域的地表中心点；绿色倒三角代表任一台站位置；红色与紫色圆点分别代表计算出的P波与S波走时计算区域的节点.(a)Δl与Δh分别表示走时场内水平方向与垂直方向网格间隔；黑色箭头表示走时场外沿一维IASP91模型入射到走时场边界的初始远震体波.(c)黑色箭头表示以任一台站为虚拟源产生的S波.(b)、(d)黑色实线代表地表与任一转换深度；蓝色射线代表直达P波；灰色方形代表任一转换点位置；红色和紫色射线分别代表直达转换点的P波和S波；紫色虚线代表转换S波；红色和紫色虚线方框分别代表离转换点最近的P波和S波走时场网格范围.Fig.2 Illustration of the FM method for calculating relative traveltimes(a), (b) and (c), (d) show the process of calculating the relative traveltimes of the converted waves in the 1-D and 3-D velocity models, respectively. The red pentagrams represent the central receiving points on the surface; the green inverted triangles represent the location of any station; the red and purple dots represent the calculated P- and S-wave traveltime field nodes, respectively. (a) Δl and Δh denote the horizontal and vertical intervals in the traveltime field, respectively; the black arrows indicate the initial teleseismic waves incident to the traveltime field boundary along the IASP91 model outside the field. (c) The black arrows represent the S-wave generated by any station as a virtual source. (b), (d) The black solid lines represent the surface and any conversion depth; the blue lines indicate the direct P-wave; the gray square means the location of any conversion point; the red and purple lines represent the direct P- and S-wave at the conversion point, respectively; the purple dashed line represents the converted S-wave; the red and purple dashed boxes represent the nearest P- and S-wave traveltime field grids to the conversion point, respectively.

最后，通过P波和S波走时场以及台站和转换点地理位置，由插值得到每个接收函数在各个转换深度P-S转换波相对直达P波的相对走时ΔtPds.对于每个接收函数，用射线追踪法(一维模型中射线追踪法不需要进行三维速度校正)得到的ΔtPds减去FM法得到的ΔtPds，获得两种方法转换波相对走时之间的差值ΔtErrors.在一维速度模型中，FM法与射线追踪法计算走时的结果本质上应该一致，因此可以通过ΔtErrors的大小来判断FM法结果的可靠性并确定计算参数.

对于走时场计算区域范围的参数设置，主要包括三维体外部大小(水平方向L与垂直方向H)和内部网格点间距(水平方向Δl与垂直方向Δh)两个部分(图3c和图2a).Guan 和Niu (2018)对网格点内部间距，用了6组不同的参数进行了测试，发现在水平方向上间隔0.05°(Δl)，深度间隔5 km(Δh)的情况下，求出的ΔtErrors相对较小.经过确认其正确性，本研究的计算中，对网格点内部间距沿用了其参数设置.

由于Guan和Niu(2018)的研究区域较小，没有对计算区域外部范围进行讨论，而是直接选取水平方向边长为10°，垂直方向800 km的走时场范围进行计算.然而，当研究区域水平范围较大时，如果选取一个可以包含整个研究区域在内的走时场进行计算，会非常耗时.因此，本研究中将整个研究区域分为若干子区域进行计算(图3a).

图3 子区域划分与走时场范围示意图蓝色三角形表示台站分布位置；红色五角星表示中心接收点的位置.(a)黑色长方形表示研究区域范围；红色虚线将该区域划分为15个子区域.(b)5个不同大小的走时场范围俯视图. 边长4°的正方形代表研究区域的任一子区域.(c) 任一走时场范围立体图. L和H分别表示走时场水平与垂直范围.Fig.3 Illustration of sub-area division and traveltime field rangeBlue triangles indicate the distribution of stations; red pentagrams denote the location of the central receiving point. (a) Black rectangle shows the study area; the red dashed lines divide the area into 15 sub-areas. (b) The top view of 5 different sizes of the traveltime field. The square with 4° represent any sub-area. (c) 3-D view of any traveltime field. L and H indicate the horizontal and vertical ranges of the traveltime field, respectively.

在选择走时场外部范围时，考虑到要将地幔过渡带包括在内，因此垂直方向的深度定为800 km(图3c中H).如果选择的计算区域水平范围太小，可能导致射线路径从侧面入射，在660-km之上传入；反之，如果水平范围太大，会因为三维速度模型空间范围不够而产生假象；同时，水平范围太大或太小也可能导致原本震中距在30°～90°的接收函数超出震中距阈值范围，受到其他震相的影响.本文选取了8°、10°、12°、14°和16°(图3c中L对应的值)共5个不同大小的走时场计算区域水平范围进行测试(图3b)，在每个计算区域内，通过FM法和射线追踪法分别计算了45505条接收函数P410s(410-km间断面P-S转换波)和P660s(660-km间断面P-S转换波)的ΔtErrors.

1.2 三维模型相对走时计算

本研究中三维速度参考模型选取的是Tao等(2018)全波形反演得到的东亚地区高分辨率P和S波三维速度模型，通过插值的方法，将三维模型中的地震波波速赋值到计算区内各个网格点上.在FM法的计算过程中，可以将三维速度模型看作特殊的只有一层的一维模型.因此，tP的计算与1.1节中介绍的方法相同，而tS的计算过程存在一定的差异.

首先选择一个范围合适的三维计算区域，从震源到该区域边界采用一维IASP91模型计算.在三维计算区域内，采用1.1节中介绍的方法，计算出tP0与tP(图2a、d).然后，对于任一台站，将该台站视为一个虚拟源，在整个走时场三维体内传播S波，计算出S波的走时场，并将三维网格中每个网格点的走时储存起来(图2c).对于转换深度d，在S波走时场内找到对应转换点附近8个网格点储存的走时信息，插值得到从台站到该转换点的S波走时，由互易原理可知该走时即为tS.最后，通过P波和S波走时场以及台站和转换点地理位置，由插值得到每个接收函数在各个转换深度，P-S转换波相对直达P波的相对走时ΔtPds(图2d).为了提高计算速度，接收函数的转换点选取为一维射线追踪得到的地理位置.

2 结果

2.1 基于一维模型的ΔtErrors和CCP叠加成像

利用45505条接收函数中的P410s和P660s转换震相，在3个不同水平范围的计算区域中(8°、12°和16°)，通过射线追踪法和FM法得到了相对走时差ΔtErrors随震中距的分布情况，并通过柱状图进行了统计(图4).

从图4a—c可以看出，当计算区域水平范围为8°时，4729个接收函数P410s相对走时差大于0.1 s或小于-0.1 s.当水平范围大于或等于12°时，几乎全部的相对走时差都在±0.1 s之间，仅有2个(12°)大于0.1 s.

图4d—f的结果显示，P660s相对走时差在各个计算区域水平范围内，大于0.1 s或小于-0.1 s的数量都占有一定的比例.其中，计算区域水平范围为8°时计算出的相对走时差，约有14.7%超出了±0.1 s的范围，最大值达到了约7.5 s.随着计算区域水平范围的增加，超出±0.1 s的比例逐渐降低，水平范围12°时，比例最低，只有约2.5%.而随着走时场水平范围继续增加，比例又开始有所增长.同时，随着水平范围的增加，计算每个走时场的运算时间以及储存走时场数据占用的空间也随之增长.

图4 一维模型中P410s和P660s相对走时差计算结果(a)、(b)、(c) 分别表示走时场水平范围为8°、12°和16°时，一维模型中基于射线追踪法和FM法得到的P410s相对走时之差分布的统计结果; (d)、(e)、(f) 分别表示走时场水平范围为8°、12°和16°时，一维模型中基于射线追踪法和FM法得到的P660s相对走时之差分布的统计结果.Fig.4 Relative traveltime differences of P410s and P660s in the 1-D model(a), (b), (c) represent the statistical results about the distribution of the relative traveltime difference of P410s obtained by ray tracing method and FM method in the 1-D model when the horizontal range of traveltime field is 8°, 12° and 16°, respectively, while (d), (e), (f) represent the corresponding statistical results of P660s.

因此，本研究在一维模型中，走时场计算区域最终选择水平范围12°，深度范围800 km，并使用与Liu等(2015)相同的CCP叠加方法进行成像.通过对比发现，基于FM法的CCP叠加成像结果与基于射线追踪法得到的结果在相同纬度基本一致.

2.2 基于三维模型的ΔtErrors和CCP叠加成像

基于一维模型的ΔtErrors结果，本文在基于三维模型的计算中，选择计算区域水平范围12°的走时场进行计算.数值检验结果表明，基于三维速度模型的相对走时差ΔtErrors的分布情况与一维速度模型时类似.P410s相对走时差在各个计算区域水平范围内，±0.5 s之间的数量占据的比例都超过了90%，其中使用12°的水平范围计算结果占据的比例最高，达到了约91.9%.同样，在12°的计算区域水平范围时，得出的P660s相对走时差在±0.5 s之间的比例最高，达到约89.3%.

图5展示了410-km和660-km间断面的P-S转换点的地理位置分布，同时通过色标表示出共享相同转换点附近参与叠加的接收函数数量.CCP叠加过程中划分的网格大小，深度间隔为1 km，水平方向0.1°×0.1°；图5所示转换点密度，是高度为1 km，半径1°的圆柱形叠加体内参与叠加的转换点数量(具体方法参见Liu et al.，2015).可以看出研究区域中绝大部分区域叠加的接收函数数量都超过了200条，靠近中部的大部分区域普遍在1000条以上，最多达到2000条左右，表明本研究最终获得的成像结果拥有较高的可信度.

图5 P410s和P660s转换点位置分布图(a)、(b)分别展示了P410s和P660s转换点位置的分布情况. 红色圆点表示转换点的地理位置，背景颜色的深浅代表转换点的密集程度，即高度为1 km，半径1°的圆柱形叠加体内参与叠加的转换点数量.Fig.5 Distribution of conversion point locations for P410s and P660s(a) and (b) show the distribution of conversion point locations for P410s and P660s, respectively. The red dots indicate the geographical location of the conversion points, and the shade of the background color represents the density of the conversion points, which is the number of conversion points involved in stacking within a cylinder with 1 km height and 1° radius.

最终的成像区域，即是图5展示的接收函数叠加数量大于200条的区域.基于FM法计算出的P-S转换波相对走时，将接收函数进行CCP叠加，获得了该区域410-km和660-km间断面的起伏形态，并通过两个间断面深度之差，得到了过渡带厚度的横向变化情况(图6a、b、c).同时，本文展示了使用相同接收函数和三维速度模型的情况下，基于射线追踪校正法的CCP叠加成像结果(图6d、e、f).东北角(图6a、d中A区域)，410-km间断面转换波波峰并不明显，接收函数数量也不是特别多，叠加结果不太稳定，推测该区域间断面深度的结果并不可靠，因此不在本文讨论范围内.从图6中可以看出，通过对比，使用相同的三维参考模型的情况下，两种方法得到的间断面起伏形态与过渡带厚度的成像结果基本一致，但是由于三维路径效应的影响，导致在范围、深度与厚度方面存在一定的差别.

FM法得到的410-km间断面深度在402～426 km之间变化(图6a)，660-km间断面深度在646～702 km之间变化(图6b)，过渡带厚度在232～282 km之间变化(图6c).660-km间断面最明显的下陷区位于东经128°—131°之间，靠近长白山火山群，呈近南北向展布，集中于北纬40°—46°之间(图6b中F区域)，该下陷区域的位置和走向与Wadati-Benioff带600 km深度的位置对应较好.在研究区域北部和西南部(图6b中G、H区域)，同样可以观察到一定范围的660-km间断面下陷区，其中G区域的660-km间断面下陷程度较深，而H区域的660-km间断面下陷较浅.下陷区F西侧约200 km处，660-km间断面存在显著的抬升，范围在东经122°—126°，北纬40°—45°之间(图6b中I区域).另外，在东部(图6b中J区域)与阿尔山火山附近区域，660-km间断面也有明显抬升.

图6 地幔过渡带成像结果左列和右列分别展示了基于FM法和射线追踪校正法的CCP叠加成像结果. 其中(a)、(d)与(b)、(e)分别表示410-km间断面和 6 60-km间断面深度的横向变化；(c)、(f)表示过渡带厚度的横向变化. 图中注释与图1a相同.Fig.6 Mantle transition zone imaging resultsThe left and right columns show the imaging results of CCP stacking based on FM method and ray tracing correction method, respectively. (a), (d) and (b), (e) show the lateral variation of the 410-km and 660-km discontinuity, respectively; (c) and (f) show the lateral variation of the transition zone thickness. The figures are annotated as in Fig.1a.

410-km间断面较为明显的下陷区位于松辽盆地(图6a中B区域)和长白山火山附近区域，其中松辽盆地410-km间断面下陷大致集中在北纬41°—46°和东经122°—126°之间，呈北东-南西走向，长约900 km，宽约200 km，下陷程度约为12 km；长白山火山西南和东北方向的下陷区，程度约10 km.大兴安岭西南与东北部 (即阿巴嘎火山和阿尔山火山附近)也可以观察到410-km间断面在小范围内的微弱下陷.在大兴安岭西南和中部，以及研究区域北部(图6a中C、D、E区域)，存在几处独立的410-km间断面微升区域.

射线追踪校正法得到的410-km间断面深度在402～428 km之间变化(图6d)，660-km间断面深度在650～704 km之间变化(图6e)，过渡带厚度在232～284 km之间变化(图6f).与FM法的结果相比，射线追踪校正法得到的660-km间断面起伏形态，在长白山火山下方的下陷程度更深一些；而在北部，下陷的范围与程度略小一些，同时整体上660-km间断面的抬升程度更小(图6e).410-km间断面的整体抬升相对于FM法略小，而松辽盆地410-km间断面的下陷更明显(图6d).这种差异主要是FM法计算转换波相对走时的时候，直接在三维速度模型中进行计算，减小了三维路径效应的影响，使得相对走时的计算结果更加接近真实情况.

与间断面深度的对比结果类似，从图中可以看出，两种方法得到的过渡带厚度的横向变化基本一致，而在增厚与减薄的程度上略有不同(图6c、f).由于410-km间断面深度整体起伏不大，过渡带厚度受660-km间断面深度影响较大，在长白山火山周围、研究区域北部、大兴安岭中部与西南部等区域，分别发现了增厚的过渡带；在松辽盆地周边、阿尔山火山周围和研究区域东部则出现了减薄的过渡带.

3 讨论

如前文所述，在使用FM法计算三维速度模型中转换波相对走时的时候，沿用了在一维速度模型中通过射线追踪法得到的转换点地理位置(以下称其为“一维转换点”).图7a绘制出了在三维模型中，射线追踪校正法所使用的一维射线路径(红色与蓝色实线)和FM法直接在三维速度模型中使用的射线路径(红色与蓝色虚线).由于路径上存在速度异常区域，虽然转换点的地理位置保持不变，但是FM法得到的直达P波与P-S转换波的射线路径与一维速度模型中得到的射线路径相比，存在一定的差异，差异的大小与射线路径经过的速度异常程度和范围有关.因此，总体而言，在相同的三维速度模型中，射线路径上速度异常区的存在，对最终得到的转换波相对走时存在明显的影响.

图7 (a)三维S波速度模型剖面中直达P波与转换P410s波射线路径. 红色和蓝色实线代表在一维速度模型中通过射线追踪分别求出的直达P波和转换P410s波射线路径；红色和蓝色虚线代表在三维速度模型中分别求出的直达P波和转换P410s波射线路径. 黑色虚线表示410-km和660-km深度位置. (b)展示了基于FM法，使用三维速度模型中的转换点地理位置求出的P410s相对走时，与使用一维速度模型中转换点地理位置求出的P410s相对走时之差的统计结果Fig.7 (a) Direct P-wave and converted P410s ray paths in the 3-D S-wave velocity model profile. The red and blue solid lines represent the direct P-wave and converted P410s ray paths derived by ray tracing in the 1-D velocity model; the red and blue dashed lines represent the direct P-wave and converted P410s ray paths in the 3-D velocity model. The black dashed lines indicate the positions of 410-km and 660-km. (b) shows the statistical results of the difference between the relative traveltime of P410s based on FM method by using the geographic location of the converted points in the 3-D and 1-D velocity models

本文对FM法使用一维转换点，在三维速度模型中计算转换波相对走时会对结果造成的影响进行了讨论.从前文所提到的45505条接收函数中随机抽取2000条，首先使用FM法，沿用一维转换点，在三维速度模型中计算P410s的相对走时.然后直接在三维速度模型中搜索转换点地理位置(以下称其为“三维转换点”)，并基于三维速度模型计算P410s的相对走时.由于实际当中，一维转换点与三维转换点距离不会太远，本文搜索三维转换点的方法是，在三维模型中，以一维转换点为中心，在边长3°的正方形范围内，水平间距每0.001°取一个点，假设其为三维转换点，计算P410s的相对走时，则这些假设的三维转换点中P410s相对走时最小的，即对应真实的三维转换点.最后将FM法使用三维转换点求出的P410s相对走时，与FM法使用一维转换点求出的相对走时相减，取其绝对值为两个转换点相对走时之差.由于FM法使用三维转换点计算出的相对走时更加接近实际情况，因此相对走时差越小，说明FM法使用一维转换点计算相对走时的结果越准确.

本文随机抽取的2000条接收函数，基于FM法计算相对走时，使用一维转换点与使用三维转换点得到结果的差别，大部分小于0.1 s(约86%)，最大值不超过0.6 s，平均差异值约为0.05 s(图7b)，即使用一维转换点得到的相对走时与真实情况十分接近，因此不需要再计算三维速度模型中转换点的位置，而是沿用一维转换点即可获得接近真实的结果.

图8展示了东经126.6°，北纬41.0°的CCP叠加成像结果，蓝色实线表示在计算相对走时的过程中，只使用一维速度模型而不进行三维校正的方法 (以下简称为1D法).从图中可以看出，基于FM法(红色实线)与射线追踪校正法(绿色实线)，P410s和P660s的叠加转换振幅接近，其中基于FM法的结果振幅略大一些，同时两种方法均比基于1D法的振幅大.基于FM法得到的振幅与1D法相比，在整个区域内的振幅平均值比例为AP410s_FM/AP410s_1D=1.368，AP660s_FM/AP660s_1D=1.127(AP410s_FM表示基于FM法的P410s振幅平均值，AP410s_1D表示基于1D法的P410s振幅平均值；AP660s_FM表示基于FM法的P660s振幅平均值，AP660s_1D表示基于1D法的P660s振幅平均值).由此可以看出，选择准确的参考速度模型并考虑三维路径效应的影响对于结果的可靠性至关重要.FWEA18模型通过全波形反演，结合多种体波震相和面波信息，同时约束了P波和S波速度结构，具有较高的分辨率.我们认为基于该模型的FM法，能够使得接收函数的CCP叠加成像结果更加可信.

图8 CCP叠加结果对比图Fig.8 Comparison of CCP stacking imaging

结果显示，在整个研究区域中，410-km和660-km间断面的平均深度与过渡带的平均厚度分别为414 km、667 km和253 km，与Liu 等(2015)和Zhang 等(2016)的研究结果相比，410-km间断面深度的结果基本一致，660-km间断面深度和过渡带厚度存在一定的差异，比前人结果减小了约10 km.图6a中显示的410-km间断面抬升与下陷的区域，与Tao 等(2018)的三维模型S波速度的异常区有较好的对应.其中C、D和E区域410-km间断面的抬升与Zhang 等(2016)的成像结果基本一致，由于这三个区域离俯冲板块的作用区域较远，同时结合新生代火山分布与全波形反演结果，本文支持Zhang等(2016)的观点，认为该区域410-km间断面的抬升，主要的原因可能是由于岩石圈拆沉，致使低温物质下沉，从而造成410-km附近温度降低，相变位置变浅.

结合Tao等(2018)的全波形反演结果，长白山—五大连池连线一带660-km间断面下陷的形态与范围和太平洋俯冲板块在过渡带底部的滞留有较好的对应关系，其中南部F区域对应于日本海俯冲带在地幔过渡带的滞留区；北部G区域对应于千岛俯冲带在地幔过渡带的滞留区；H区域的660-km间断面下陷则可能与岩石圈物质拆沉有关.660-km间断面最明显的下陷区域(图6b中F区域)，下陷幅度约为20～40 km，与Wadati-Benioff带的600 km轮廓线形态基本一致，位于深源地震带的延伸方向上.这一结果与Liu 等(2015)和Zhang 等(2016)的研究结果相似，但是在范围与深度方面有一定的差异.更深的660-km间断面相变位置通常归因于较低的过渡带温度，根据实验室测试结果，660-km间断面的克拉伯龙斜率变化范围在-1.3～-2.8 MPa/K之间(Ito and Takahashi，1989；Fei et al.，2004)，以-2.0 MPa/K为例计算，取图6b中F区域660-km间断面平均下陷深度30 km进行计算，则该区域660-km间断面的下陷对应大约-500 K的温度变化.Tao 等(2018)全波形反演结果显示在该区域存在P波高速异常区，与本文中下陷范围有较好的对应，异常值约为2%，考虑到地震波速对温度的敏感度约为-0.43%/100K(Cammarano et al.，2003)，则该高速异常对应的温度变化约为-465 K，与本文估算的结果大致一致.图6b中F区域的下陷范围向西延伸到约东经128°，横向范围约300 km，大致说明了俯冲的太平洋板块在地幔过渡带的滞留范围，与前人(如Zhao，2004；Huang and Zhao，2006)显示的滞留板块近水平向西延伸约1500 km的结果并不一致.

沿北纬42°由西向东，在深度剖面中(图9c)可以较为清晰地识别出P410s和P660s转换波，从而判断出两个间断面大致的深度分布.410-km间断面在东经116.4°附近有小幅度抬升(图9c中红色椭圆位置)，在松辽盆地与长白山火山下方区域明显下陷(图9c中蓝色椭圆位置)；660-km间断面在大兴安岭与长白山火山附近区域有明显下陷(图9c中绿色椭圆位置)，在松辽盆地下方区域显著抬升(图9c中黄色椭圆位置).深度剖面的结果与前文所描述的间断面深度异常情况有较好的对应关系，同时更加直观地显示出了部分异常所处的位置、范围以及间断面起伏的程度.

图9 (a)地表地形起伏示意图. 浅蓝色表示海平面以下；红色三角形表示长白山火山群. (b)红色与蓝色虚线分别表示对应地理位置410-km和660-km间断面参与叠加的接收函数数量；黑色虚线的位置表示叠加数量为200条.(c)沿北纬42°的CCP叠加深度剖面. 黑色虚线分别表示410-km和660-km深度位置；红色表示转换波振幅大于0的部分Fig.9 (a) Topographic relief map. The light blue color indicates below sea level; the red triangle indicates the Changbaishan Volcanic zones. (b) The red and blue dashed lines represent the number of receiver function involved at the corresponding geographic locations of 410-km and 660-km discontinuity, respectively; the position of the black dashed line indicates the number of receiver function involved in stacking is 200. (c) The depth profile along 42°N. The black dashed lines indicate the position of 410 km and 660 km, respectively; the red color represents the part of the conversion wave with amplitude greater than 0

研究区域内660-km间断面明显抬升的位置分别位于图6b中的I、J和阿尔山火山附近区域，与全波形层析成像结果S波的低速异常区有较好的对应，说明这三处位置660-km间断面附近可能有高温异常存在.其中最显著的抬升区I，位于俯冲板块前端西侧，其南北向狭长的形状与F区域类似，空间上近似平行，抬升幅度约为10 km，推测可能与太平洋俯冲板块的破坏或下沉有关，从而导致局部下地幔高温物质上涌.同时长白山火山的岩浆活动，也有可能是热物质沿俯冲板块前端破裂缺口上涌形成的.上涌的热物质受到浅层结构的影响，在410-km间断面附近形成不同的热物质上升区，造成了松辽盆地及长白山火山410-km间断面附近温度上升、深度加深，为该区域高温物质来源与长白山火山起源问题提供了约束.图6b中与F和I区域空间上近似平行的J区域660-km间断面抬升，则可能是由于板块后端破坏所造成的高温物质上涌，致使660-km间断面附近温度升高，相变位置变浅.阿尔山火山附近区域660-km间断面的抬升，也可能是受到深部热流作用而形成的.由于阿尔山火山与J区域均位于成像区边缘，其图像需要补充更多数据，结合更大范围的成像结果进行确认.

本文主要讨论了温度与地幔过渡带速度间断面深度之间的关系.410-km和660-km间断面的深度除了受温度影响以外，还可能受到水含量和化学成分等因素的影响，从而使得实际情况变得更加复杂.许多矿物物理学实验与地球动力学模型表明，古老冷却的海洋板块中含水矿物能够将水带入地幔过渡带内，停滞于过渡带内的海洋板块可以将水释放到周围的地幔中，从而对410-km和660-km间断面深度产生很大影响(如Thompson，1992；Helffrich，2000；Higo et al.，2001；Harte，2010).一些研究表明，水的存在对地幔过渡带内间断面深度的影响有着类似于低温的作用，尤其是过渡带底部，水的存在对660-km间断面下陷的影响甚至比温度降低还要大，俯冲板块中2.0 wt%的水含量可能会引起660-km间断面下陷约15 km (Litasov et al.，2005；Helffrich，2000；Pearson et al.，2014).也有实验证据表明化学成分对间断面深度的影响，比如，含铁量对410-km间断面深度的影响与含水量相似，这意味着在富含铁元素的区域，410-km间断面深度将变得更浅(Katsura et al.，2004).我们会在未来的研究工作中对这部分内容进行更加详细的分析与讨论.

4 结论

本文使用了中国东北地区NECESSArray记录到的远震事件接收函数，采用快速行进法计算P-S转换波相对走时，通过与射线追踪校正法得到的相对走时以及成像结果的对比，说明了FM法在大尺度区域研究中的可行性.基于东亚地区高精度三维速度模型，通过FM法得到了中国东北地区地幔过渡带的三维结构.成像结果显示，在长白山—五大连池连线一带，有明显的660-km间断面南北向下陷和过渡带增厚区，东西方向延伸约300 km.同时在该区域东西两侧，发现了在空间上与下陷区近似平行延伸的660-km间断面抬升与过渡带减薄区，其中660-km间断面最显著抬升区位于长白山西侧约200 km.表明日本海俯冲带在地幔过渡带中近水平滞留的延伸距离有限，而且可能造成地幔热物质沿其所造成的破裂缺口上涌，提供了长白山等新生代火山活动的岩浆源.而松辽盆地西侧660-km间断面的下陷，反映的更有可能是该地区岩石圈物质拆沉的结果.研究区域内410-km间断面显示出的所有下陷区，与地表新生代火山活动区/拉伸区有很好的对应关系，表明这些地表构造与深部热物质上涌有关，同时也表明地表构造可能对深部热物质上涌通道有塑造作用.

致谢感谢NECESSArray国际合作项目和中国地震局地球物理研究所“中国地震台网数据管理中心”提供的波形数据.感谢中国地震局地球物理研究所柳正博士提供的接收函数数据.感谢匿名审稿人提出的很好的建设性意见.本研究工作得到北京大学高性能计算校级公共平台支持.