基于三重震相波形非线性反演的俯冲带410-km间断面起伏研究

李嘉琪， 宁杰远， 蔡晨， 鲍铁钊

北京大学地球与空间科学学院， 北京 100871

0 引言

410-km间断面的精细结构，特别是俯冲板块内部及周围410-km间断面的起伏情况涉及到地幔对流模式、俯冲带滞留状态、深源地震形成机制等重要科学问题，受到地球科学家的广泛关注.

Helffrich等(1989)以及Kirby等(1991)指出俯冲板块内部橄榄石(olivine)到瓦兹利石(wadsleyite)的平衡态相变界面的深度可抬升至350 km.Thirot等(1998)用接收函数方法研究了410-km间断面在太平洋俯冲板块内部的起伏情况，发现其抬升到了350 km的深度，与理论预测的抬升情况一致.但是，Tonegawa等(2005)用接收函数方法研究了伊豆—小笠原俯冲带410-km间断面的起伏情况，发现只有约30 km的抬升.Niu等(2005)也用接收函数方法研究了伊豆—小笠原俯冲带附近410-km间断面的起伏情况，发现俯冲板块周围410-km间断面存在明显的抬升，但在俯冲带内部却没有清晰的410-km间断面信号.

Collier和Helffrich(1997)利用俯冲带内向下传播的地震体波的转换波资料，通过倾斜叠加的方法，发现伊豆—小笠原俯冲带内410-km间断面最高抬升到了350 km的深度.蒋志勇等(2003)利用类似的方法得到了同样的结论.但是，Revenaugh和Jordan(1991)利用ScS回折波研究了410-km间断面的起伏，认为在俯冲板块内部有小幅度抬升.

另外，Vidale和Benz(1992)利用台阵观测的反射波资料进行倾斜叠加，认为俯冲板块内部410-km间断面有不超过30 km的抬升，而Flanagan和Shearer(1998)的SS前驱波结果没有观测到西太平洋俯冲带地区的410-km间断面的抬升.

同时，Lidaka和Suetsugu(1992)利用体波走时残差分析，认为日本海俯冲带存在亚稳态橄榄石楔.Kawakatsu和Yoshioka(2011)利用接收函数叠加，同样发现日本海俯冲带下方存在的亚稳态橄榄石楔.Jiang和Zhao(2011)以及Jiang等(2015)利用双差地震走时分析，也认为在日本海俯冲带存在亚稳态橄榄石楔.这时，410-km间断面不升反降，与相变动力学的预测一致(见Sung and Burns，1976).

三重震相是组合震相，在高速间断面附近射线覆盖很密，非常适合研究地幔间断面的精细结构(如Brudzinski and Chen，2003；Gao et al.，2006；Zhang et al.，2008；Wang and Chen，2009； Wang and Niu，2010；Zhang et al.，2012; Chu et al.，2012).Wright和Kuo(2007)曾利用这一方法研究吕宋下方俯冲板块内部的410-km间断面起伏，认为抬升到了325 km的深度.但是，他们只用了走时进行研究，反演结果有较强的非唯一性.

本文利用NECESSArray密集宽频带地震台阵记录的发生在千岛俯冲带的地震所产生的三重震相波形资料，对俯冲板块内部及附近410-km间断面的结构进行非线性反演.密集的台阵资料，使得能够通过对比射线回折点通过或不通过俯冲板块内部时的情况.射线回折点在鞑靼海峡下方，此处射线路径方向与俯冲板块的走向大致一致，实现了最大限度地对受到俯冲板块影响的410-km间断面进行采样，克服小波长起伏不易识别的困难，给出了俯冲带内部及附近410-km间断面结构的稳定结果.

1 研究区域与资料

本研究的地震资料来源于在中国东北地区进行的持续两年的野外流动地震观测项目NECESSArray (NorthEast China Extended SeiSmic Array).该台阵共有流动台站127个，平均台间距约为80 km，东西跨度约1300 km，南北跨度约为600 km.图1中用黑色三角符号展示了本文所用的36个台站的位置.

我们根据国际地震中心(International seismological center，ISC) 提供的地震目录，选择了一个射线回折点附近的射线路径大致平行于俯冲板块走向的地震，该地震发震时刻为2009年10月10日21时24分(GMT时间)，震级为Mw5.9.所用资料中，该事件到NECESSArray的震中距范围为12°～20°，方位角范围为264°～280°.图1中用白色圆点展示了射线回折点(Turning Point，射线穿透最深处)的位置.由图1可以看到，这些回折点位于鞑靼海峡附近区域下方.图1中画出了射线在地表的投影，本文对北区(275°～280°，绿色曲线所示)、中区(269°～274°，红色曲线所示)、南区(264°～266°，蓝色曲线所示)三个方位角范围不同的地区进行对比研究.在具体反演过程中，我们先去除仪器响应，得到地震位移记录.然后进行一阶的巴特沃斯零相移滤波，其中的带通滤波范围为0.05～1 Hz.

为了得到更准确的地震发震深度，我们通过IRIS网站的breq_fast网页脚本，下载了该地震发生后1 h的连续记录来对发震深度进行精定位.为了保证数据质量，利用Crazyseismic(Yu et al.，2017)软件挑选出100多个信噪比高的记录.最终挑选出18个信噪比高、深度震相又清晰的记录.在读取了P、pP以及sP在各个远震台站的到时之后，固定经纬度不变，通过全局搜索确定震源深度.当震源深度是114 km时，残差达到最小值0.8 s.此外，作为对比，再采用GCMT(Dziewonski and Woodhouse，1983)给出的震源机制解(走向/倾角/滑动角分别为128°/11°/191°)，通过直达波以及深度震相的波形拟合，来反演地震深度(McCaffrey and Abers，1988).最终结果如图2所示.由波形反演得到的震源深度也是114 km，与利用到时搜索得到的结果一致.

图1 研究区域与所用地震、台站分布图中红色震源球为所用地震，黑色三角形为所用台站，绿、红、蓝的曲线代表北、中、南三区射线路径在地表投影，白色散点是射线穿透最深处的地理位置.黑色虚线为俯冲带等深线，黑色数字为等深线所对应的深度.Fig.1 Research region and the distribution of stations and eventsThe red beach ball represents the Mw5.9 Kuril subduction zone event, and black triangles show the stations used. Green, red and blue curves are the projections of the ray paths for northern, middle and southern regions, respectively. White dots represent the location of the turning points for the ray paths. Black dashed lines are the depth contours of the subduction zone, and the black numbers show the depth for the contours.

图2 震源机制解体波波形反演图图正中为震源球下半球的等面积投影，压缩象限为红色.其中三角形为台站所在位置示意，邻近区域内邻近台站用同个三角表示，18个台站共分为黑、黄、蓝、绿四个区域.蓝色实线波形为实际记录，红色虚线波形为正演结果，黑色竖杠代表残差计算时窗，波形左上方字母代表台站名，颜色与震源球中三角相对应.波形右上方数字代表残差值.左上角三角形为震源时间函数.波形拟合反演得到的震源深度为114 km.Fig.2 Focal depth inversion from depth phasesIn the equal-area projection of the lower hemisphere of the focal sphere, the compressional quadrant shaded in red and the 18 stations used are concentrated in 4 areas represented by the black, yellow, blue and green triangles. Also shown are the vertical-component records (solid blue lines) and synthetics (dashed red lines). The vertical black bars indicate the time windows for the residual calculation. For each waveform trace, the station name is given on the left with the same color as the triangle in the focal sphere, while the waveform residual is given on the right. The source time function is given in the top-left corner. The best fitting depth is ～114 km.

图3 410-km间断面对应的三重震相示意图.(a) 三重震相射线路径示意图，其中红色、绿色、蓝色实线分别为直达波、反射波、透射波，在410-km间断面附近形成了密集的射线覆盖； (b) IASP91模型(Kennett and Engdahl，1991)P波速度图，图中小框内展示的是小于120 km深度的浅部模型； (c) 三重震相理论波形图，其中黑色波形为用QSEIS (Wang，1999)正演的理论位移地震图，红色、绿色、蓝色实线为用Taup(Crotwell et al.，1999)计算得到的走时曲线图.AB，BC，CD定义与文中一致，黑色圆圈圈定的O点为AB支与CD支的交点位置.Fig.3 Overview of triplications from the 410-km discontinuity(a) Ray paths of P-wave triplications. The solid red, green and blue lines are direct waves, reflected waves and refracted waves, respectively; (b) P-wave velocity from IASP91 model (Kennett and Engdahl, 1991), the subfigure shows the structure shallower than 120 km; (c) Synthetic waveforms showing triplications of the 410-km discontinuity. The black waveforms are synthetics calculated by QSEIS (Wang, 1999), and solid red, green and blue lines are travel time curves calculated by Taup (Crotwell et al., 1999). AB, BC, and CD branches represent the direct waves, reflected waves and refracted waves, respectively. The O point in the black circle shows the crossover point of the AB and CD branch.

2 研究方法

2.1 三重震相方法简介

如前所述，本文使用的是三重震相方法.当地震体波(如P波)在地球内部传播，遇到地球内部的高速间断面(如图3b所示的410-km间断面)时，波的传播方式会发生改变，因此有三类传播路径不同的体波(如图3a所示):射线路径在间断面上方的直达波AB，间断面的反射波BC以及透过间断面下方传播的透射波CD.在一定震中距(如P波探测410-km间断面时是17°左右)台站的地震记录图上同时记录到多个震相并形成随震中距规律变化的三重震相(如图3c所示)，组成了能够约束间断面结构的可测量波形组合.此外，三重震相的三条射线在浅部经过基本相同的路径，浅部的异常结构会对三个震相有相似的影响；而在间断面附近三条路径分开，深部的结构会对三个震相有不同的影响.因此使用三重震相不仅可以减小浅部未知结构对深部结构反演的影响，还可以根据三重震相之间到时、振幅的相对变化，形成对深部间断面附近结构的有效约束.尤其是一维结构反演，因为反演参数少，反演结果的确定性非常好.

2.2 小生境遗传算法

早期的三重震相反演研究，为了减小台站稀少带来的限制，多选择单台多地震方案进行研究，但需要大量的高质量地震.2000年以后，随着宽频带密集地震台网的建设，仅一个地震就可以得到完整的走时曲线以及波形图，所以有条件选择信噪比高、受干扰少，震源函数简单的地震进行研究.

Johnson(1967)首先利用三重震相的走时信息反演上地幔结构，发现了410 km与660 km附近的两个高速间断面，并与Anderson(1967)的固-固相变理论计算相一致.后来，研究者通过比较理论地震图与实际记录的三重震相波形信息，分析波形异常，人工修改模型，获得更好的拟合(如Tajima and Grand，1995; Brudzinski and Chen，2003).但是要通过改变这些参数来拟合复杂的波形，需要丰富的经验及大量的时间.Gao等(2006)将共轭梯度法引入到三重震相的反演，利用计算机代替人去搜索模型.但是，试错法和共轭梯度法都面临着落入局部极值的风险，也存在最终结果依赖于初始模型的问题.

为了避免试错法反演对研究者经验的要求，以及共轭梯度法反演对初始模型的依赖与落入局部极值的风险，我们将一种基于搜索类反演的小生境遗传算法应用到了三重震相反演中(Koper et al.，1999；李少华等，2012).小生境遗传算法，继承了遗传算法的三个优点：其一是不依赖于初始模型，只需要给出模型的搜索范围，遗传算法将在此范围内随机地产生大量模型，并通过“杂交”和“筛选”进行之后的模型更新；其二是在模型迭代的过程中，允许“变异”的存在，模型中的参数在每一次迭代中，均有1%至5%的概率改变其数值，以此来避免落入局部极值；其三是由于搜索了大量的模型，遗传算法可以最终输出一系列可接受的模型集合，这些可接受模型的平均值和方差可以帮助估计最终模型的不确定度.此外，优于传统遗传算法的是，小生境遗传算法将解空间划分为十个 “群落”，每个“群落”之中模型参数相似，“群落”之间的模型存在差异，方便考察问题的多解性.为了克服搜索类算法计算代价大的问题，我们设计了CPU并行策略并编写了程序.CPU多核心在每一“代”内进行并行计算，极大地提高了反演效率.当使用100个CPU核心时，同时计算每一“代”内的100个模型，可以达到并行的最大效率.在本文反演中，我们使用了25个CPU核心，20,000次正演与反演搜索可以在4小时内完成.

本文设计了理论测试.首先，比较了使用IASP91(Kennett and Engdahl，1991)模型时，本文使用的QSEIS程序计算出的震相到时和Taup程序(Crotwell et al.，1999)计算出的到时是否一致,以及QSEIS计算的波形和DSM(Geller and Takeuchi，1995)方法计算出的波形是否一致.我们分别计算了震源在100 km至600 km，每隔100 km的情形.我们还改变了矩张量，测试了使用爆破源和双力偶源时各自的理论地震图.QSEIS在30度震中距之内的震相到时和波形分别与Taup及DSM的结果一致.

验证了QSEIS的正确性之后，我们进一步验证了小生境遗传算法的正确性和稳定性.我们根据IASP91模型计算理论位移波形图.在图4a蓝色虚线给出的搜索范围内，每隔30 km设置一个速度值的待反演点，其间的速度由相邻的两个反演点插值得到.需要注意的是，除了常规的速度值的待反演点外，我们还设置了反映间断面性质的三个待反演参数，分别是间断面的位置以及间断面两侧的速度值.利用小生境遗传算法进行20,000次搜索后，得到的最优模型与输入的IASP91模型十分接近(如图4a所示)，并且结果的收敛性很好.我们以最终模型集合的平均值作为正演的模型，计算理论地震图，与IASP91模型的理论波形对比.可以看到，其符合程度很高(如图4b所示).该方法收敛速度也很快，10代之后残差明显减小，50代之后再次减小，80代之后趋于稳定(如图4c所示).

2.3 台阵观测波形的整体归一化

为了避免台站条件差异对计算结果的影响，过去的三重震相反演计算中多对振幅进行单台归一化计算.本文对振幅采用对台阵中所有台站的记录波形进行整体归一化处理，可以保留台站之间的振幅相对变化信息，以形成更多的约束.为了避免个别台站资料质量不好的影响，可以采取降低权重或个别分析的方法加以处理.

下面用两个模型正演结果的比较来说明进行波形整体归一化的必要性.图5a展示的黑色模型是IASP91模型，红色模型是设计的仅在浅部存在低速的对比模型.如图5b所示，从波形整体归一化振幅中可以看到，浅部不同的结构除了带来整体到时移动外，主要对直达波振幅造成影响，后至波振幅基本未变，见图5b中黄色椭圆形虚线框所示区域.但是当进行单台归一化振幅处理时，由于在三重震相交点之前的震中距范围内直达波振幅总是最大的，归一化后，直达波振幅一直是“单位一”.而由于单台归一化，原来振幅未变的后至波，在归一化后振幅减小，如图5c中蓝色矩形虚线框所示.由于后至波对应的是410-km间断面处的反射波及其下的透射波，因此不管是手动试错法，还是自动搜索，都会错误地去调整深处的速度结构，见图5a中黄色与蓝色矩形虚线框所示意.因此可以看出，使用单台归一化的振幅，不仅损失了台站间的波形变化信息，还有可能使错误地理解波形中的异常来源，影响反演结果.因此，使用整体归一化振幅是必要的.

2.4 基于三重震相到时的初步分析及消除浅部影响的“先对齐、后反演”计算方案

我们手动提取了按方位角划分的三个区域中(见图1)各个台站所记录到的震相到时，对各个台站、各个震相分别减去Taup程序根据IASP91模型计算得到的理论到时后，得到其到时差.

如图6所示，全部三个区域的所有台站的震相到时均慢于IASP91模型，滞后时间在13°震中距附近达到7 s，这很可能是浅部存在大规模的低速异常所致(见图7的层析成像结果)，也可能有地震震源位置以及时间不准的贡献.

此外，在北区与中区，直达波的滞后时间以17°为界，在13°～17°之间随震中距增大而减小，大于17°后滞后时间基本不变；南区在17°之前似乎有同样趋势，但在17°后观测不到直达波.

13°～17°之间台站的直达波到时差最大有2 s的差异.这可能对应于三种典型的情形或其组合：第一种对应于200 km(13°的射线穿透深度)至320 km(17°的射线穿透深度)之间的区域存在高速异常；第二种对应于台站下方的浅部低速异常存在差异；第三种对应于震源区附近由于俯冲板块存在而出现的高速异常.17°之后，北区与中区直达波的滞后时间基本不变，反映320 km至410-km间断面之间不存在显著的异常或路径上异常的贡献被抵消了.而南区17°之后没有观测到直达波，可能反映在410 km深度上方存在高速间断面.

图4 小生境遗传算法反演模拟数据测试图(a) 反演模型图，黑色实线为IASP91模型，红色实线为最终可接受模型集合，蓝色虚线框为模型搜索范围； (b) 波形拟合图，黑色波形为IASP91模型的理论地震图，红色的为最终模型的平均模型对应的理论地震图； (c) 残差随遗传代数变化图.Fig.4 Synthetic tests for Niche Genetic Algorithm(a) Inverted model results. Black solid line is the IASP91 model, red lines are inverted acceptable model sets, and blue dashed lines show the model searching range; (b) Waveform fitting results. Black waveforms are synthetics for IASP91 model and red waveforms are synthetics for the inverted model; (c) Residual between data and synthetics with respect to generations.

图5 单台归一化振幅与全局归一化振幅的比较(a) 速度模型图，黑色实线为IASP91模型，红色实线为浅部低速的模型.黄色虚线框表示两模型存在速度梯度变化的位置，蓝色虚线框表示当采取了单台归一化后可能会错误修改模型的区域示意； (b) 全局归一化振幅示意图，黑色波形为IASP91模型的理论地震图，红色为浅部低速模型的理论地震图，黄色虚线框标识波形振幅差异的部分； (c) 单台归一化振幅示意图，蓝色虚线框标识波形振幅差异的部分，其他图示同(b).Fig.5 Comparison between trace normalization and array normalization(a) P-wave velocity profile. Black solid line is the IASP91 model, red solid line is a model with a low velocity zone in the shallower part. Yellow dashed box shows the region where the velocity gradient changes between two models. Blue dashed box roughly shows the region where we tend to modify when applying trace normalization; (b) Array normalized waveforms. Black waveforms are synthetics for the IASP91 model and red waveforms are synthetics for the red model in (a). Yellow dashed oval shows where the waveform amplitudes are different; (c) Trace normalized waveforms. Blue dashed oval shows where the waveform amplitudes are different. Other symbols are the same as (b).

体波层析结果及面波层析成像结果(如Fukao et al., 2001; Tao et al., 2018; Kang et al.，2016)都显示，在图7中震中距大于15°的台站下方，南区相比于中区速度更低，可能与长白山下方存在的热物质有关.但是，在图6中，中区与南区在15.7°附近的直达波到时却完全一致，这表示南区相比于中区，在该射线穿透深度附近存在着高速异常体，这与图7中南区的俯冲板块存在深度比中区更浅的相对应.

图6 三重震相走时拾取图横坐标为震中距，纵坐标为与IASP91模型预测到时之差.圆形代表北区的拾取结果、三角形代表中区的拾取结果、正方形代表南区拾取结果.其中实心标记代表410-km之上直达波的拾取、空心标记代表410-km之下透射波的拾取.Fig.6 Travel time analysis for triplication The x-axis is epicentral distance and the y-axis is the travel time differences compared with the IASP91 model′s predictions. Circles show the results for the northern region, triangles show the results for the middle region and squares show the results for the southern region. The solid marks represent the direct waves above the 410-km discontinuity and the hollow marks represent the refracted waves below the 410-km discontinuity.

全部三个区域透射波的滞后时间小于直达波的滞后时间，表示透射波与直达波射线之间的区域存在高于IASP91模型的高速异常，这与图7层析成像结果中震源区以及410-km间断面附近的高速俯冲板块相对应，也可能对应于410-km间断面的抬升.

比较北区与中区，直达波的滞后时间基本一致，但透射波与直达波的滞后时间之差，在北区约为2 s，在中区则约为3 s.这表示中区的直达波与透射波射线之间的部分相比于北区速度更高，这与图7层析成像结果里，射线穿过中区俯冲板块的部分比北区俯冲板块更多相对应，也可能对应于中区的410-km间断面比北区存在更大的抬升量.

北区与中区的透射波的滞后时间随震中距增大而逐渐减小，但是南区透射波的滞后时间基本保持不变.这可能是因为，在北区与中区随着震中距增大，射线逐渐更多地经过深部的俯冲板块.而南区可能由于俯冲板块更浅，射线已全部位于俯冲板块内部，因而不会呈现出随震中距增大，穿过俯冲板块更多的现象；也可能由于俯冲板块内部存在亚稳态橄榄石，该低速的亚稳态橄榄石区，部分抵消了俯冲板块因低温而带来的高速异常.

从走时资料中可以直接获取许多重要信息，但为了更好地确定结构，还需要通过波形拟合进行非线性反演.然而在三重震相一维反演研究中，缺乏对浅部及源区结构的约束.若直接对绝对到时进行反演，浅部的未知结构会对深部结构的反演造成污染(李嘉琪等, 2016).有的学者利用已有的浅部速度模型进行修正，可以一定程度上减少浅部结构的污染(Ye et al.，2011)，但当浅部结构不准确时，污染仍然存在.

先将理论波形与观测波形对齐，然后进行反演是一种可以选择的策略(Grand and Helmberger，1984；LeFevre and Helmberger，1989；Brudzinski and Chen，2003; Wang and Niu，2010; Chu et al.，2012；Zhang et al.，2012).这样忽略绝对到时的对齐操作可以很大程度上减小浅部结构的影响，凸显对三重震相之间相对到时更为敏感的深部结构.以下是我们对这种方法的合理性进行进一步的理论测试.

2.4.1 关于浅部未知一维结构影响的测试

理论测试中的参考模型是IASP91模型.对于待反演模型的小于50 km的浅部区域，我们固定其速度始终大于IASP91模型0.4 km·s-1，以代表不准确的浅部结构(图8c).由于反演中所使用的第一个台站的震中距为13°，其穿透深度约为200 km，因此我们只反演200 km以下的速度结构.

在计算残差时，我们一般选择P波波列之前5 s作为时窗的起始，P波波列之后2 s作为时窗的结束，如图10中红色理论波形持续时间所示.对于误差函数的计算，我们先对每一个台站的观测波形和理论波形通过波形互相关，获得到时差Δti.然后根据获得的到时差Δti将理论波形与实际记录进行对齐.随后再选取对齐后的时间域理论地震图与实际记录地震图的残差的L2范数作为误差函数L2:

(1)

其中，d(xi,t)为台站i记录到的观测位移记录，u(xi,t+Δti)为台站i延时Δti后的理论地震图，t1与t2分别为三重震相时窗的始末时刻，N为总台站数.

反演结果如图8b所示.可以看到，深部的结构可以得到很好恢复.当采取了对齐操作后，每道记录三重震相之间的波形信息被有效地凸显，同时台站间振幅的变化可以帮助约束深部结构.具体地，两种模型虽然浅部结构有所不同，但唯有在深部的结构相同时，才可以使得同一台站的各个波形之间以及各个台站之间波形变化的理论与“观测”结果一致.需要指出的是，忽略了绝对到时也会增加解的非唯一性.但这种非唯一性，主要是以模型的整体移动为主，但模型的相对形态并不会有很大改变(Gao et al.，2006).

2.4.2 关于台站下方二维结构影响的测试

在到时分析中，我们还观测到13°和17°的台站直达波到时差之间存在2s的差异.暗示浅部速度可能存在非均匀分布的情形.理论测试已经表明，三重震相由于其浅部射线路径相似，受浅部一维结构不准确的影响较弱.但能否将这种可忽略浅部影响的优势拓宽到浅部二维非均匀情形，必须进行进一步的理论测试.

我们在震中距小于13.5°的区域内设置厚为50 km的低速异常区，该区域P波速度低于IAPS91模型0.3 km·s-1，具体模型如图9a所示.从图9c中有限差分方法(Li et al., 2014)计算的黑色正演地震图可以看到，震中距小于13.5°的台站的到时相比于震中距大于13.5°的台站有一个明显的整体滞后，这是台站下方的低速异常带来的影响.

观察反演的波形，绿色虚线的波形是QSEIS计算的一维波形(红色实线)进行了互相关对齐后得到的.除去整体振幅略微的偏小，和二维模型产生的波形(黑色实线)符合很好.从反演模型可以看到，模型虽有些误差，但基本上收敛于IASP91模型，受浅部二维非均匀影响很小.这来源于两方面的共同约束：其一是台站间振幅变化带来的约束.虽然二维模型在浅部和一维模型差别很大，但在深部仍是IASP91模型.唯有深部模型一致，才可以拟合台站间的振幅变化.具体地，13.6°的台站相比于13.3°的台站，有一个明显的时间相对提前，但波形振幅没有很明显的变化，这暗示着这个时间提前是浅部的影响.若此变化来自深部，则主要来源于13.6°的射线和13.3°的射线穿透深度之间的区域，但这个局域的高速将极大地增大13.6°台站波形振幅，而实际波形中并未出现此现象；第二种约束是来源于每道记录三重震相本身，即直达波与反射波、透射波之间的相对到时.若异常来源于浅部，三重震相的三个震相将整体地移动；若异常来源于深部，则对三个震相的到时影响不同.具体地，若异常位于250～410 km范围内，对直达波的影响最大；若异常位于410 km之下，对透射波影响最大.不论哪种，都会改变震相之间的相对到时，无法同时拟合这三类震相.

图7 该地区层析成像结果图(a) 北区层析成像结果1(Fukao et al., 2001)，上方为地形图，下方为红蓝色表示由层析成像得到的速度异常，其中黑色三角形代表台站，红色五角星代表地震位置，黑色细实线为根据本文反演模型计算的射线路径，黑色虚线为IASP91模型中410-km与660-km间断面位置，桃红色虚线为波形反演得到的‘410-km’间断面位置； (b) 中区层析成像结果1，标识同(a); (c) 南区层析成像结果1，标识同(a)； (d) 北区层析成像结果2(Tao et al., 2018)，标识同(a)，由于Tao等(2018)研究区域未到达千岛俯冲带东北端，因而本文所用震源附近部分区域无层析成像2的结果； (e) 中区层析成像结果2，标识同(a); (f) 南区层析成像结果2，标识同(a).Fig.7 Tomography results in this region(a) Tomography result 1 (Fukao et al., 2001) for the northern region. The upper panel is the topography, the lower panel is the 2-D velocity variations. The red star and black triangles represent the locations of the earthquake in the Kurile subduction zone and stations in northeast China, respectively. Ray paths are calculated using one of the final acceptable models. Pink dashed line shows the depth of the 410-km discontinuity obtained in this study; (b) Tomography result 1 for the middle region; (c) Tomography result 1 for the southern region; (d), (e), and (f) are tomography result 2 (Tao et al., 2018) for the northern, middle and southern region, respectively.

图8 浅部一维不准确结构对反演的影响(a) 波形拟合图，黑色波形为IASP91模型的理论地震图，红色为最终平均模型的理论地震图，绿色虚线是红色波形互相关平移后的波形； (b) 反演模型图，黑色实线为IASP91模型，红色实线为最终模型集合； (c) 浅部速度模型图，黑色实线为IASP91模型，红色实线为反演模型被固定的浅部部分.Fig.8 Synthetic tests for inaccurate 1-D structure at shallow depth(a) Waveform comparison. Black and red lines are synthetics for model IASP91 (black line in b) and one of the acceptable models (red lines in b), respectively, whereas the green dotted lines are the red synthetics after aligning with the black ones by cross correlation. (b) 1-D inverted models. The solid black line is the IASP91 model, and the solid red lines are the final acceptable model sets from the inversion. (c) Shallow portions of the 1-D models. The solid black line is the IASP91 model, and the solid red line is the inverted model with the shallow part adjusted so that the velocity is 0.4 km·s-1 faster than the IASP91 model.

图9 浅部二维非均匀模型及反演结果图(a) 理论测试设计的2D模型图，黑色五角星为震源、黑色三角形为台站位置示意，黄色长方形区域为13.5°之前的低速异常区.红色、绿色、蓝色曲线分别为直达波、反射波、透射波的射线路径； (b) 反演模型图，黑色实线为二维实际模型150 km之下的部分(同IASP91模型)，红色实线为反演得到的模型集合； (c) 反演测试波形图，黑色波形为二维模型正演的理论地震图，红色为反演的平均模型对应的地震图，绿色虚线是红色波形互相关平移后的波形.Fig.9 Synthetic tests for inaccurate 2-D structure at shallow depth(a) The 2-D model used in the synthetic test. The background model is IASP91, and the velocity in the yellow region shows a -0.3 km·s-1 low-velocity anomaly. The black star is the earthquake, and the black triangles are the stations. The maximum epicentral distance influenced by the low-velocity anomaly is around 13.5°; (b) 1-D inverted models. The solid black line is the reference IASP91 model, and the solid red lines are the final acceptable model sets from the inversion; (c) Waveform comparison. The black lines are the synthetics for the 2-D reference model in (a) using a 2-D finite difference algorithm (Li et al., 2014) and the red lines are the synthetics for one of the acceptable models in (b). The green dashed lines are the red synthetics after aligning with the black ones by cross correlation.

图10 垂直分量实际观测波形与理论波形图(a) 北区理论波形与实际观测波形对比图，黑色波形为实际波形(已与理论波形通过互相关对齐)，红色波形为反演模型对应的波形(其长度代表反演所取时窗)，蓝色实线为Taup计算的走时曲线； (b) 中区理论波形与实际观测波形对比图，标识同(a)； (c) 南区理论波形与实际观测波形对比图，蓝色虚线是Taup计算的走时曲线延伸，其他标识同(a).Fig.10 Vertical-component displacement seismograms corresponding to ray paths in the (a) northern, (b) middle and (c) southern regionsThe horizontal axis is reduced time, and the vertical axis is epicentral distance. The red waveforms in (a), (b) and (c) are synthetics for one of the acceptable models from the inversion. The duration of the synthetics represents the misfit window in the inversion. The black waveforms are records after alignment with the synthetics by cross correlation. And blue lines show the travel-time curves of the direct, reflected and refracted waves.

图11 间断面抬升及间断面两侧速度差增大情形正演测试图(a) 间断面抬升30 km时，正演波形图，其中黑色虚线为走时曲线； (b) IASP91模型正演波形图，标识同(a)； (c) 间断面速度差增加60%时，正演波形图，标识同(a)； (d) 三种波形所对应模型图，其中黑色为IASP91模型，蓝色为间断面抬升30 km模型，红色为间断面速度差增大抬升60%的模型.蓝色、黑色、红色模型的波形图分别对应于(a)、(b)与(c).Fig.11 Modeling tests on how the model influences the cross-over point O(a) Waveforms for the model with a +30 km uplift of the discontinuity. Black dashed lines are traveltime curves for the model. (b) Waveforms for the IASP91 model. (c) Waveforms for the model with a +60% velocity jump across the discontinuity. (d) Corresponding models. Black line is the IASP91 model, blue line is for the model with a +30 km uplift of the discontinuity, and red line is for the model with a +60% velocity jump across the discontinuity. And the blue, black and red models correspond to (a), (b) and (c) respectively.

图12 波形反演结果图(a) 北区反演结果，黑色实线为IASP91模型，绿色实线为反演模型集合； (b) 中区反演结果，黑色实线为IASP91模型，红色实线为反演模型集合； (c) 南区反演结果，黑色实线为IASP91模型，蓝色实线为反演模型集合.Fig.12 P-wave velocity inversion results obtained in this study in the three cross sections (a) north, (b) middle and (c) southGreen, red and blue lines show the inverted acceptable model sets in this study, whereas the black line indicates the IASP91 model.

图13 模型trade-off反演测试图(a) 固定间断面位置370 km时，波形反演结果图，黑色实线为实际资料，红色实线为反演波形，黑色箭头指示拟合不好处； (b) 固定间断面位置380 km时，波形反演结果图，标识同(a)； (c) 固定间断面位置390 km时，波形反演结果图，标识同(a)； (d) 固定间断面位置400 km时，波形反演结果图，标识同(a)； (e) 固定四种深度时对应的反演模型图，其中黑色实线为IASP91模型，红色实线为固定四种深度时，反演得到的速度模型，其中红色粗线代表可接受模型，红色细线为拟合程度较差的模型.Fig.13 Synthetic tests for the trade-off of the model parameters(a) Waveform comparison when the discontinuity is fixed at 370 km. The black waveforms are the data and the red waveforms are the synthetics. Black arrow points at the waveform which has large mismatch; (b) Waveforms comparison when the discontinuity is fixed at 380 km; (c) Waveforms comparison when the discontinuity is fixed at 390 km; (d) Waveforms comparison when the discontinuity is fixed at 400 km; (e) The inverted models for these four cases with fixed discontinuity depth. The black line is the IASP91 model, the thick red lines are the acceptable models and the thin red lines are the models with worse waveform fitting.

以上的两个测试表明，我们开发的基于小生境遗传算法的三重震相一维反演方法可以正确并快速地得到反演结果，并且适用于某些浅部存在横向非均匀的情形.

3 鞑靼海峡下方410-km间断面的起伏

3.1 理论波形与实际观测波形的对比

图10展示了如前所述的北(方位角范围为275°～280°)、中(方位角范围为269°～274°)、南(方位角范围为264°～266°)三个区域台站的对齐后的实际观测波形(黑色曲线)和用反演得到的模型计算的理论波形(红色曲线)及其走时曲线.波形图是按震中距排列的，横轴为折合走时.需要指出的是，我们首先观察到了不同方位角的台站观测波形存在明显的差异.因为不同，只有分成三个区，才能对观测波形进行最好拟合.最重要的是，这种差异恰恰反映了射线路径与俯冲带有不同的几何关系，通过对比，正好用于研究俯冲板块内外间断面深度的差异.三个区的边界，由观测波形与正演波形的符合程度决定.具体地，我们先选取一个较窄的方位角范围进行反演(1°)，得到一个模型.对此模型进行更大方位角的正演并与实际资料相比较，当观测波形与正演波形有较大差异时，将此观测波形划分到下一个区域内.

可以看出，三个区域的观测波形随方位角变化而变化，即从北到南，三重震相交点位置O向更小震中距移动.由图11a、图11c的正演测试可知，O点出现的震中距减小可能对应着410-km间断面的抬升或间断面两侧速度差的增大.此外，OB支的振幅从北到南，逐渐减小.由图11b、图11c所示的正演测试可知，间断面两侧速度差增大，并不会改变OB支的振幅；而410-km间断面的抬升，会导致OB支振幅减小.因而从北到南，可以定性地判断间断面位置应存在逐渐抬升.具体的抬升量，以及间断面两侧的速度差的定量计算，需要由波形反演所确定.

通过与理论波形对比，发现北区和中区近震中距实际观测波形的第二个峰幅度偏小.我们认为，可能是由于间断面存在局部起伏所造成的.在南区，由于近震中距缺少台站且远震中距波形振幅较小，对模型的约束程度不如北区与中区.

3.2 鞑靼海峡下方410-km间断面附近的一维波速结构

图12为上面提及的三个区域(方位角分别为275°～280°、269°～274° 和264°～266°)的410-km间断面附近,反演得到的一维波速结构.可以看到三个区域的间断面均存在抬升，抬升量逐渐增大，从北区的10～20 km，到中区的20～30 km，至南区的60～70 km，与3.1节中定性的分析一致.并且，所有模型都显示，在220 km以下100 km左右的深度范围内存在小幅度高速异常；间断面两侧速度差(从北到南依次约为10%、10%、7%)均比IASP91模型的3.5%的速度差大；间断面下方(约500 km附近)，在北区与中区存在另一个高速异常的极大值，但南区未有此现象.

需要注意的是，反演结果存在一定的不确定性.以北区为例，间断面在400 km深度附近的模型，与间断面在390 km深度附近但其上方存在低速度梯度的模型会有类似的波形拟合效果.这是因为虽然抬升的间断面会使得OB支的振幅较小，但当间断面上方速度梯度降低时，OB支会延伸到更远的震中距，因而在相同震中距下振幅增大.

为了考察此种模型非唯一性的程度，我们以中区为例，分别固定间断面深度为370 km、380 km、390 km、400 km，考察在此固定位置下反演得到的对应的模型速度结构，以及波形拟合情况.从图13e中的反演模型结果可以看到，越浅的间断面位置会伴随着间断面上方更低的速度梯度值，与我们定性的分析一致.同时我们可以发现，这种不确定性的影响有限.具体地，间断面在390 km时波形拟合程度最好(图13c)；当间断面深度变化过大时，13度的后至波拟合程度变差(图13a、图13d).考虑到波形中误差的存在，我们认为中区内，间断面位于380 km与390 km的模型的拟合程度均可接受，因此最终的模型的深度不确定度为约10 km.其他区域情况类似，这对应于图12中每个区域内存在两类模型，其间断面深度大约有10 km的差异.

需要注意的是，在目前的台站分布、噪音水平以及P波主频2.5 s的情况下，我们无法区分这两种模型的优劣，因而我们无法确认410-km间断面上方是否存在低速区.尽管间断面深度存在约10 km的不确定度，但并不会对间断面抬升的结论产生大的影响.

4 结论与讨论

4.1 410-km间断面两侧的速度异常

反演结果中，间断面上方至220 km范围内的速度普遍大于IASP91模型的值.以射线覆盖较好的北区与中区为例，在波形拟合的意义上，此结果主要对应于近震中距台的较其他台更大的振幅.具体地，北区与中区的NEAD与NE9E台站(射线穿透深度约为200 km)中的AB支的振幅约为IASP91模型对应波形的1.5倍.更大的波形振幅，对应于模型中该深度更大速度梯度的存在.因而在反演结果中该处存在一个高速度梯度区.再往深，虽然速度值依旧增大，但速度梯度与IASP91模型类似.值得注意的是，由于我们使用互相关来消除绝对到时的影响，虽然模型的相对变化可以被约束，但整体上与真实模型会存在一个平移.因而只能说明在220 km附近存在高速度梯度区，但对其绝对速度无法约束.Fukao等(2001)与Tao等(2018)层析成像的结果(图7)一致地展示出，该区域浅部存在大规模的低速，且此低速终止于200 km附近.我们反演得到的高速度梯度带，对应于由较低波速变为正常波速的深度范围.

反演结果中，间断面两侧速度差(北区、中区约为10%，南区约为7%)普遍大于IASP91模型的值(3.5%).在到时分析中，直达波与透射波之间存在的大于IASP91模型约3 s的到时差，也支持间断面两侧更大的速度差.

间断面两侧的速度差是由三重震相的相对到时和波形确定的，且由于是一维反演，待反演参数较少，模型确定性较好.但同样由于是一维反演，此速度差的具体数值，会受到二维结构的影响.

具体地，利用三重震相进行一维反演时，一般把异常都归给了深部区域.这是因为三重震相的射线在浅部路径接近，在穿透深度处分开.并且，由于在转折点处行进路径最长，因而对此区域速度变化最为敏感.此假设在大多数情况下是成立的.但在震源区附近，因为俯冲带的存在，有很强的局部高速异常.并且此俯冲带对应的高速区异常大致平行于射线离开震源时的行进方向(如图7)，其累积效应不可忽略.因而，源区的高速异常可能会被折合到穿透深度附近，因而造成对该区域速度值的高估.一种可能的解决方案是利用李嘉琪等(2016)提出的借助于层析成像结果修正的反演方案来尽量消除此影响.在本文中，我们采取另一种较为简单与半定量的估计方式.具体地，我们参考其他研究者的层析成像结果(如Fukao et al.，2001；Tao et al.，2018)，扣除层析成像给出的源区异常值后，估计间断面下方的速度.

不同研究者的层析成像结果对于俯冲板块具体位置的估计存在差异(图7).具体地，Fukao等(2001)模型中俯冲板块位置较Tao等(2018)的结果更浅，但在中区二者的位置较一致.中区的间断面速度跃变为10%.尽管震源区的俯冲板块会对此异常值有贡献，但即使扣除了源区3%左右的高速贡献，仍剩余7%.中区显著大于IASP91模型的速度跃变可能来自于物质相变与俯冲板块内外温度差异的双重贡献.这也与图7b、图7e中三叉震相反演得到的间断面位置与层析成像结果中高速区域上表面接近相对应.进一步扣除了层析成像显示的目标区3%的速度异常的影响后，剩余约4%的波速跃变可能由橄榄石-瓦兹利石的相变所产生，与IASP91模型的速度跃变结果接近.

对于北区，俯冲板块更深，但我们仍然观察到了7%的剩余速度差.扣除了约4%的由橄榄石-瓦兹利石的相变所产生的速度跃变后，仍剩余3%的高速异常.这表明着此处的地幔转换带中，410-km间断面之下、俯冲板块上表面之上的区域存在高速异常体.图7d中Tao等(2018)的层析成像结果在射线转折点处(距震源8°～9°)确实存在局部高速结构，尽管幅度较弱.

对于南区，由于缺乏更多的尤其是近震中距的台站资料，对间断面的约束较弱，不利于讨论速度跃变值.

反演结果中，间断面下方存在一个速度梯度减小的区域.在波形拟合的意义上，此结果主要对应于CD支更小的振幅.具体地，震中距16°～18°的台站的初动震相波峰不再尖锐.较为平缓的透射波波形，对应于模型中其射线(CD支)穿透深度(约430 km)附近的低速度梯度区.Fukao等(2001)与Tao等(2018)层析成像的结果(图7)一致地展示出，在射线最敏感的转折点位置(距震源6°～8°)，400 km至更深的区域，存在低速异常区，此处的原本高速的俯冲板块也存在“间断”的现象，与本文反演得到的该处速度梯度的降低一致.但源区附近俯冲板块的二维结构，也会导致穿出其中进入正常地幔的射线振幅减小.因而，如要对此速度梯度减小量进行更为准确的估计，需要进行考虑了源区结构的二维或三维反演.

4.2 410-km间断面的抬升

反演结果表明410-km附近的间断面存在抬升，与矿物物理学的预测结果一致(Katsura andIto，1989; Bina and Helffrich，1994; Katsura et al.，2004).并且反演得到的间断面从北区的10～20 km，到中区的20～30 km，至南区的60～70 km，对应的平均位置分别为395 km、385 km、345 km.进一步，我们可以对抬升的间断面的温度进行具体估计.根据橄榄石平衡态相变界面的实验结果(Bina and Helffrich，1994)，未抬升的410-km间断面处温度为1400 K，当410-km间断面分别抬升15 km、25 km、65 km时，其温度相比于周围地幔物质降低约185 K、300 K、765 K.因此北区、中区、南区410-km间断面处温度约为1215 K，1100 K，635 K.

不同研究者的层析成像结果均显示俯冲板块从北到南逐渐变浅(图7).对于北区，俯冲板块上表面位于410-km间断面下方，温度效应最弱，410-km间断面抬升量最小；对于中区，间断面位于俯冲板块上表面附近，温度效应较明显，因此抬升量增大.

对于南区，在Fukao等(2001)的模型中，其俯冲板块上表面浅于本文反演得到的间断面.在此情形下，反演得到的365 km附近的间断面应是处于俯冲板块内部中上位置(图7c)的抬升的410-km间断面，推测俯冲板块的冷的温度效应最为明显，因而抬升量最大.但在Tao等(2018)的模型中，俯冲板块上表面位于本文反演得到的间断面附近(图7f).因此，南区反演得到的间断面也可能是俯冲板块的上表面.

总地来说，目前的研究结果表明俯冲带内部似乎不存在大量的亚稳态橄榄石.但是，进一步观察南区的波形，可以看到在17°到18°台站的透射波与直达波之间，似乎存在“中间震相”(图10c).但由于南区资料较少，波形振幅也较弱，因而我们仅拟合了初至震相与后置震相.今后将基于更多的资料，并对更多细节进行拟合，进一步研究上地幔结构细节.

致谢本文所用资料来自北京大学、美国、日本合作布设的NECESSArray观测数据.在论文完成过程中，使用了汪荣江教授的QSEIS程序，和钮凤林教授、王曙光博士、吴庆举研究员、陈望平教授、岳汉研究员、周莹教授、陈棋福教授，陈敏教授以及Peter Shearer, Steve Grand等教授进行了讨论.感谢Ross Maguire博士对英文摘要的语言方面帮助.本研究工作得到北京大学高性能计算校级公共平台支持.也感谢三位匿名评审人给出的细致又有建设性的意见.