1 引言

自中生代开始，多个块体(如羌塘块体、拉萨块体等)向欧亚大陆南部拼贴，新特提斯洋关闭，印度与欧亚发生陆陆碰撞.在印度大陆持续向北推进的作用下，欧亚大陆南侧和拼贴块体的地表抬升，进而 形成青藏高原(Haines et al., 2003; Yin and Harrison， 2000). 欧亚大陆同印度板块间近南北向1400~2350 km巨大的缩短量(Harrison et al., 1992)，使得地壳厚度显著增加(70~80 km)，达正常大陆地壳厚度(Christensen and Mooney， 1995)的两倍，这一重要特征已被大量地球物理探测证实(滕吉文等， 1997，2012；吴庆举和曾融生，1998; 李秋生等，2004; Zhang et al., 2011;Nábělek et al., 2009; Shi et al., 2004; Tian et al., 2005; Tseng et al., 2009; Jiang et al., 2006; 高锐等2011).然而，众多研究表明：青藏高原的演化不是简单连续的，而是存在明显的时空脉动性变化，如21~18 Ma和8 Ma以来的两期高原大范围快速抬升(Harrison et al., 1992)，又如青藏高原火山岩随时间向北再向南的多期次分布(Chung et al., 2005; Ding et al., 2003).由于岩石圈不同的动力学过程如低角 度俯冲(underthrusting)、板片回退(roll-back)、板 片断离(break-off 或detachment)、 高角度俯冲(subduction)、对流移除(convective removement)，对应不同的地表形变模式和不同的火山时空及成分变化，基于此地质学家试图通过对这些事件分析研究，重建新生代青藏高原下方印度岩石圈、青藏高原岩石圈及欧亚大陆岩石圈的动力学过程(Chung et al., 2005; Ding et al., 2003; Molnar et al., 1993; Chen and Tseng， 2007)，所以探测现今地幔岩石圈的几何形态，有助于陆陆碰撞和青藏高原隆升的动力学过程重建和高原演化研究(吴福元等，2008；Guillot and Replumaz， 2013).

区域固定台网数据层析成像研究显示印度岩石圈低角度俯冲在青藏高原下(Huang and Zhao， 2006)，上地幔顶部回折波研究(McNamara et al., 1997)指示藏北上地幔顶部存在一个明显的Sn不发育的高衰减区(主要集中在羌塘和松潘—甘孜块体下方)，并被认为是缺乏地幔盖层，地壳与软流圈直接接触(Owens and Zandt， 1997)；体波远震层析成像也表明在藏北上地幔存在明显的低速异常(Liang et al., 2012; Tilmann et al., 2003; Wittlinger et al., 1996)，但整个高原下方地幔岩石圈在汇聚过程中如何形变以及目前的空间位置和几何形态的研究远未达成共识.

近20年来，中外科学家在青藏高原已开展多条近南北向的宽频带地震探测.在85°E以西的西部，接收函数界面成像(Zhao et al., 2010)和远震走 时地震波速度成像(Li et al., 2008，Zhao et al., 2014)显示印度岩石圈低角度俯冲在青藏高原下方；在东部(90°E以东)，接收函数界面成像(Kosarev et al., 1999; Zhao et al., 2011)显示印度岩石圈高角度俯冲至拉萨块体或班公—怒江缝合带下方；在85°E至90°E之间的中 部地区，接收函数界面成像研究(吴庆举等，2004; Zhao et al., 2010)和远震走时地震波速度成像(Zheng et al., 2007)显示印度俯冲板片小角度向北俯冲至羌塘块体下方; 而Tilmann等(2003)的远震P波走时地震波速度成像结果则显示汇聚增厚的岩石圈在班公—怒江缝合带南侧近垂直下沉(downwelling).这些不同研究方法获得的结果存在差异的原因可能来自多个方面，如岩石圈底界面深度较大导致其对应的地震信号信噪比较低，再如岩石圈结构在东西方向上可能存在明显的分带性而高原内部观测密度却是有限的(比如 Liang et al., 2012)，以及高原下方的岩石圈与软流圈之间地震波速度差异较小.

尽管INDEPTH-Ⅲ宽频带天然地震观测已完成十多年，并已发表多项地震学探测结果，但对该剖面下方的地幔岩石圈结构的研究仍存在较大争议.Tilmann等(2003)通过远震走时地震波速度成像获得了剖面下方P波速度扰动图像，结果显示汇聚增厚的岩石圈在班公—怒江缝合带南侧近垂直下沉；而在相邻的位置，Zheng等(2007)的成像结果却显示印度俯冲板片小角度向北俯冲至羌塘块体下方，显然该地区需要开展更多的地幔岩石圈探测研究工作.考虑到S波速度对温度变化比较敏感，本文利用INDEPTH-Ⅲ剖面S波波形记录，使用多道互相关准确提取相对走时，通过层析成像获取剖面下方上地幔S波速度扰动图像，为青藏高原动力学研究提供更多约束. 2 方法与原理

地学层析成像是将医学X射线CT理论应用于调查地下物性参数分布状况的地球物理探测技术(Aki and Lee， 1976; Aki et al., 1977; Culter et al., 1984)，其主要目的是确定地球内部的精细结构和局部非均匀性(杨文采，1993).地震层析成像通过对观察到的地震波各种震相的运动学(如走时、射线路径)和动力学(如波形、振幅)特性的资料进行分析，进而反演出地下介质的速度及其他物性参数结构等重要信息的一种地球物理方法(滕吉文等，2004).本文主要利用远震S波走时数据来反演上地幔介质的速度结构.

地震波走时是地震波能量沿射线经过不同速度区域所需时间的积分.一般，定义地震波的走时残差为地震波沿射线穿过实际地球介质和穿过参考模型之间走时的差异，它主要反映的是地震射线所经过区域的速度同参考模型之间的差异；其近似等同于沿射线所经过区域的慢度异常的积分值，同时受到地震震源定位误差和发震时刻误差等因素的影响.远震的传播路径在靠近震源的一侧基本相同，而在台站下方靠近台站部分存在明显的差别，因此不同台站所记录的远震震相波形相似，台站下方速度结构的横向非均匀性对台站间的相对走时残差起主要贡献，使用远震相对走时基本排除了研究区域外因素对反演研究区内速度异常的影响.所以，本研究中使用远震S波的相对走时残差作为约束，对台站下方的速度结构进行层析成像.

我们将研究区域沿剖面(28°N—36°N，80~400 km 之间)划分为反演所使用的规则四边形网格(0.25°×20 km)，共693个节点，并假设四边形网格节点上的慢度异常值为Δs_m，n.射线i经过第j个网格的慢度异常Δs_ij(中点(x_j，y_j))可以根据该网格周围四个节点(m，n)、(m+1，n)、(m，n+1)、(m+1，n+1)的慢度异常插值得到

射线i经过网格j的长度l_ij可以由相应节点坐标计算获得.每条射线的相对走时残差t_i可以表示为：

c_imn表示(m，n)节点对第i条射线的影响系数，是节点周围所有网格内射线长度l_ij及其在该网格内权重((1)式)的乘积之和，Δs_mn表示在(m，n)节点上的慢度扰动.

在方程组(3)中，节点上的慢度异常为待求解的未知变量，已知条件包括地震射线路径构建的系数矩阵和相对走时残差构建的数据向量.

这里我们使用阻尼最小二乘QR分解(least squares QR-factorization，LSQR)方法(Paige and Saunders， 1982)求解该线性方程组，进而得到INDEPTH-III剖面下方的S波速度扰动异常分布.

在计算初始射线路径时，本文使用ISAP91模型的S波速度模型作为初始模型，然后通过反演出来的速度扰动对初始模型进行更新，并对射线路径进行重新计算，直到慢度扰动结果最终收敛.

人工拾取震相到时的精度受人员熟练程度的影响很大，拾取误差也不容易估计.为了弥补这些缺点，本研究中采用多道互相关的方法计算台站间的 相对走时残差(Hales and Doyle， 1967，Vandecar and Crosson， 1990).该方法的另外一个优点是对信 号的信噪比要求较低(Vandecar and Crosson， 1990)，所以我们在筛选数据的时候对震级没有进行特别的限制.

在对接收到同一个地震的任意两个台站的地震波形进行互相关分析时，将截取的S波震相在±5 s进行扫描计算互相关系数，认为互相关系数最大(大于0.8)的时间差是两个台站(第i个和第j个台站)之间的相对走时残差Δt_ij，如果第i个台站相对于初始模型的相对走时残差记为t_i，那么

考虑到地震波形中总是有噪声的影响，互相关得到的相对走时也会受到噪声的影响，这时为了对上述方程组(4)进行求解，需要引入一个到时差之和为零的假设作为约束：

由方程组(4)和方程(5)组成的方程组可以方便地用最小二乘法求解，从而得到每条射线的相对走时残差(Vandecar and Crosson， 1990). 3 数据和结果

本文用到的数据是中美等国合作于1998年7 月至1999年6月间在青藏高原中部布设的INDEPTH-Ⅲ 剖面(西藏德庆—龙尾错剖面)记录的宽频带地震波形数据.整条剖面共60个宽频带地震台构成，从拉萨块体跨越班公—怒江缝合带进入羌塘块体内部(图 1a)，全长大约400 km.考虑到整条剖面基本沿着北北西走向在一条线上分布，而该区域内的主要构造带为近东西向延伸，我们忽略垂直测线方向的横向非均匀性.为了减小垂直于测线方向的射线可能引入的横向非均匀性的影响，我们仅使用反方位角与测线夹角在45°以内的地震事件(图 1b).在这些地震中，进一步挑选震中距在30°~150°之间，有清晰的S、SKS、PKS、ScS等S波震相的波形，并对数据进行0.05~0.25 Hz的带通滤波处理.为了保证数据的可靠性，要求每个地震事件至少被8个台站同时接收，并且剔除走时残差大于4 s的一些不合理的台站，确保排除仪器故障等特殊原因导致的严重异常.我们共筛选出49个宽频带地震台记录的30个远震事件，共计622条射线.

图1

Fig. 1

图 1 台站以及应用到的事件的分布情况 (a)中三角形表示台站的位置，(b)中空心圆表示INDEPTH-Ⅲ接收到的所有事件的位置，实心圆表示我们用到的事件位置. Fig.1 The distribution of stations and the event applied (a)The triangles shows the location of the stations.(b)The hollow circular shows the location of the events recorded by INDEPTH-Ⅲ，the solid circular shows the location of the events we used.

以1998年7月31日12时40分59秒发生在169.82°E，21.61°S震级M_w6.0的地震事件为例.通过对经过理论到时校正后的波形图和经过相对走时残差校正后的波形比较，可以明显发现经过相对走时残差校正后，所有台站的波形在到时和形态上都比较一致(图 2a)；通过多道互相关方法得到的相对走时残差表现为：在32°N以南相对走时残差为负值，以北为正，说明该事件的地震射线经过的32°N以南地区的速度可能大于32°N以北地区(图 2b).

图2

Fig. 2

图 2 多道互相关计算得出的相对走时残差 (a)中左图是将所有台站的理论到时进行校正后的波形，右图是进行相对走时残差校正后的波形. Fig. 2 Residual of travel times calculated by multichannel cross-correlation (a)shows the waveform before and after corrected by residual.

因为远震射线在接近地表的时候入射角比较小，射线交叉少，对地壳内的速度异常不能提供有效约束；并且，拉萨块体与羌塘块体之间地壳结构存在明显的差异(Tseng et al., 2009)，使用均匀地壳模型对远震走时残差进行校正显然是不合理的.所以我们根据Mechie等在2004利用INDEPTH Ⅲ宽角反射地震数据反演得到的研究区地壳S波速度(图 3)对多道互相关获得的相对走时残差进行校正，得出了排除地壳速度异常影响的相对走时残差(图 2).如图 3中所示，地壳速度结构在32°N班公—怒江缝合带附近非均匀性比较强.对地壳中的速度异常进行校正以后得到的80 km以下速度结构引起的相对走时残差(图 2)表明，32°N附近相对走时残差负值变大，深部的高速体相对于浅部可能更靠北.

图3

Fig. 3

图 3 INDEPTH Ⅲ宽角地震数据得到的该区域 地壳中的S波速度(Mechie et al., 2004) Fig. 3 Crustal S wave velocity reference model along the profile INDEPTH Ⅲ provided by Mechie et al.(2004)

由于反演问题本身的多解性，在反演过程中选择不同的反演参数可以获得多样的反演结果(张风雪等，2013)，这就决定了确定最终的反演所使用的参数具有一定的经验性.地震波走时层析成像的分辨率很大程度上依赖于研究区域的地震波射线覆盖和相互交叉情况(田有等，2009).本研究中实际使用的射线分布如图 4所示，在30°N—34°N之间射线密度最大，中间的大部分区域节点附近的射线总长度达到600 km以上，而射线覆盖范围随着深度从浅到深而增加，但是射线密度随深度增 加而减小，射线交叉比较好的区域大约在150~300 km之间.同时，我们对模型进行了检测板测试(Humphreys and Clayton， 1988)(如图 5)： 给定±10%速度异常相间的速度模型，并对正演计算的走时引入时间标准差0.1倍的随机噪声来模拟实际数据中存在的噪声干扰.两组大小不同的检测板实验表明，在射线密度比较大的范围内0.5°×80 km的检测板恢复情况比较清 楚，输入模型的特征基本能恢复(图 5a，5c)；而0.5°×40 km 的检测板在射线交叉少的模型两侧200~400 km深的位置对输入模型的恢复明显变差(图 5b，5d).

图4

Fig. 4

图 4 射线密度分布图 颜色深浅代表所有射线经过该节点附近射线的总长度. Fig. 4 Ray density distribution The darkness means the total length of all rays pass through the near nodes.

图5

Fig. 5

图 5 检测板实验 红色代表低速异常，蓝色代表高速异常.(a，c)的分辨率为0.5°×80 km，(b，d)的分辨率为0.5°×40 km. Fig. 5 Checkboard tests Red and blue colors denote low and high velocity abnormaly，respectively. The resolutions in(a，c)are horizontal 0.5° and 80 km in depth domain，while the resolution in depth domain decrease to 40 km in(b，d).

最小二乘反演中的阻尼系数起到调和模型参数 分辨率和误差之间关系的作用(Paige and Saunders， 1982; Nolet，1987). 不同阻尼系数的反演结果中，模型光滑程度与走时残差之间的关系如图 6a所示.经过多次试验，我们选取最佳的阻尼系数为10.经过四次迭代以后速度结构趋于稳定.每次迭代后的结果与上一次迭代结果之间的差异(其中第1次相对于均匀模型)如表 1所示.反演前后的相对走时残差的分布情况如图 6b所示，反演之后走时残差明显减小，同时残差的方差也明显减小.

图6

Fig. 6

图 6 不同阻尼系数的反演结果中，光滑程度与走时残 差之间的关系(a)和反演前后相对走时残差分布的对比(b) Fig. 6 The relationship between smooth degree and the residual travel time by different damping coefficient shows in(a). The contrast of the distribution between before and after inversion shows in(b)

表1(Table 1)

表 1 迭代次数与速度模型稳定性的变化趋势 Table 1 The relationship between the total times of iterations and the stability of the velocity model change trend 迭代次数 1 2 3 4 与上一次的误差 0.3753 0.0093 4.27×10-4 1.43×10-5

表 1 迭代次数与速度模型稳定性的变化趋势 Table 1 The relationship between the total times of iterations and the stability of the velocity model change trend

为避免反演图像中相邻网格点之间速度扰动的异常跳变，我们对反演图像进行了光滑处理，以获得我们所关心的大尺度结构图像.结果中每一个格点的异常值用它周围相邻的八个格点和它本身值的加 权平均代替(权系数取为：本身1/4，边上的临点1/8，角上的点1/16).结果如图 7所示：一个高速体自拉萨块体中部(约31°N)高角度向北并向深部延伸，倾斜角度约为65°，可分辨的深度可达300 km.

图7

Fig. 7

图 7 INDEPTH-Ⅲ剖面下S波慢度异常层析成像的结果 蓝色表示慢度小，代表高速异常，红色表示慢度大，代表低速异常. Fig. 7 S tomographic imaging results under INDEPTH-Ⅲ profile Blue color indicates small(negative)slowness and represents high seismic velocity，while red color indicates large(positive)slowness and represents low seismic velocity. The block arrows are the possible dynamic interpretations of downwelling and upwelling.

在进行地震波走时层析成像时，分布不合理的射线可能导致反演得到的异常形态同异常的实际形态之间存在明显差异，因此，有必要进行可恢复实验(Zhao et al., 1992)来验证我们所使用的数据和方法是否具有区分最终结果中得到高速异常为高角度异常或垂直异常的能力.为此，本文进行了三种不同形态高速异常模型的可恢复实验，如图 8所示，三个输入模型中分别包含竖直展布的高速异常体、向北倾斜65°的高速异常体和垂直向下的倒梯形高速异常体，输入S波速度异常的幅度设定为5%.使用与实际数据相同的震源分布和台站分布，通过正演计算获取模型所对应的走时异常，并加入走时异常标准差0.1倍的噪声，通过层析成像反演获得不同输入模型所对应的恢复结果.从恢复结果看(图 8)，在300 km深度以上，三种输入模型的基本特征都能够得到恢复并加以区分.

图8

Fig. 8

图 8 三个不同形态高速异常模型的可恢复实验 蓝色代表高速异常，红色代表低速异常.分别给出三种形态不同的高速异常(a)、(c)、(e)，恢复结果如图(b)、(d)、(f)所示. Fig. 8 Recoverable tests of anomaly model with three different shapes Three high velocity anomaly models(a)，(c)， and (e)are the input，while(b)，(d)， and (f)are output，respectively.

从图 8可以看出，在使用向北倾斜65°的高速异常体作为输入模型进行可恢复实验时，高速体北侧仅出现范围有限的微弱虚假低速异常，异常幅度和分布范围远不及实际数据反演结果(图 7)中班公—怒江缝合带以北羌塘块体上地幔的低速异常.为了检验该低速异常的可靠性，我们进行了另一种可恢复实验(图 9)，在倾斜高速异常体的北侧并排设置一个低速异常(正异常和负异常的幅度同样设定为 5%)，从恢复结果看，正负异常的基本特征在300 km 深度以上可以很好地恢复.两组可恢复实验表明，本文结果中高角度倾斜S波高速体及其上方和北侧的S波低速异常都是可信的异常体.

图9

Fig. 9

图 9 同时存在高速异常和低速异常倾斜模型的可恢复实验 Fig. 9 Recoverable test for model with a high and a low velocity anomaly coexisting

反演得到的S波速度成像结果(图 7)显示一个高速体自拉萨块体中部(约31°N)高角度向北并向 深部延伸，倾斜角度约为65°，可分辨的深度可达300 km.该S波高速体形态明显同Tilmann等(2003)获得的P波高速体存在差异(在32°N—33°N之间深度小于100 km的区域，我们的S波结果呈现低速异常，而Tilmann的结果表现为高速异常).这两者的差异可能主要来源于以下两个方面，首先是射线覆盖程度的不同造成的分辨率差异，因为P波通常有比S波更密集的射线分布；其次，S波和P波对温度异常和物性差异的敏感度不同，S波速度对温度异常更敏感，而P波速度对物质成分差异更敏感(Boyd et al., 2004).P波和S波层析成像分别在同样地区发现了高速异常，说明该异常的可靠性.Tilmann等(2003)认为他们反演的异常主要 是因为印度大陆岩石圈不断向北汇聚增厚，因为力 学不稳定性或者之前高密度的海洋岩石圈的拖曳，导致增厚的印度大陆岩石圈发生重力下沉. 同样的解释也可以用来帮助我们理解反演得到S波速度异常图像，这个异常也可能是岩石圈不断增厚进而发生拆沉后的残留部分.