上地壳纵横波速度结构相关反演成像方法

王晓， 周小鹏， 张新彦， 白志明， 滕吉文

1 中国科学院地质与地球物理研究所，岩石圈演化国家重点实验室， 北京 100029 2 中国科学院大学， 北京 100049



上地壳纵横波速度结构相关反演成像方法

王晓1,2， 周小鹏1,2， 张新彦1,2， 白志明1*， 滕吉文1

1 中国科学院地质与地球物理研究所，岩石圈演化国家重点实验室， 北京 100029 2 中国科学院大学， 北京 100049

基于纵横波初至走时数据的层析成像方法越来越广泛地被应用于揭示不同构造域壳幔速度结构特征.我们从同一地质体的纵横波速度属性相关这一基本思想出发，提出一种相关反演成像的方法：纵横波速度反演交替进行，在迭代反演过程中每通过一次反演获得相应的纵波速度(或横波速度)结构后，更新相应的纵横波速度比模型以及相应的横波(或纵波)速度反演的初始模型，然后继续开展后续横波(或纵波)速度反演工作.在反演过程中依据纵横波速度的相关性信息和射线路径长度将走时残差以不同权重分配到射线路径经过的单元，依据网格节点周围平均的慢度扰动更新速度模型.正反演过程分别基于有限差分走时计算方法和反投影成像方法.两种典型模型试验表明，该技术应用于上地壳速度结构反演成像过程，可有效提高反演结果的可靠性，在很大程度上避免了常规单独反演纵波和横波速度过程容易带来的畸变和失真.该方法应用于重建青藏高原西部札达—泉水沟深地震测深(DSS)剖面下方的上地壳速度结构，揭示出与青藏高原西缘板块碰撞相关的上地壳速度结构特征.

上部地壳； 纵横波速度结构； 相关反演； 有限差分； 反投影

1 引言

自20世纪80年代以来，基于深地震测深纵横波初至走时数据的层析成像结果越来越广泛地被应用于揭示盆岭、裂谷、断陷等不同构造的深部形态及上地壳内高速岩体、岩浆囊和断层等结构的精细特征(Tarantola and Valette, 1982; Thurber, 1983; Vidale, 1988; Bregman et al., 1989; Saito, 1989; Hole and Zelt, 1995; 张先康，1996；Zhang and Toksozf, 1998; 丁志峰等，1999; Xu et al., 2006,2010,2014; Teng et al., 2013).Vidale(1988)最先提出了二维介质的有限差分正演模拟算法，Hole (1992)将其改进，推广应用于速度横向变化剧烈的三维介质，并提出了一种利用初至波走时数据Pg和/或Sg开展层析反演的方法.由于有限差分方法的通用性，该方法自提出以来不断得到应用和改进.经过Ammon和Vidale (1993)改进后，该技术已被广泛地应用于利用深地震测深(DSS)数据重建上地壳速度结构领域.王椿镛等(1997)应用该方法得到大别造山带人工地震测深剖面上部地壳横向不均匀结构图像，为大别造山带演化过程的研究提供了重要的深部地球物理证据；段永红等(1999)详细论述了二、三维层析成像的基本原理，并以唐山震区的实际反演结果验证了方法的有效性.之后，越来越多的学者将该方法用于研究地壳结构和构造特征，并取得了一系列的成果(Zelt, 1999; Bai et al., 2007; 王夫运等，2008；Thybo and Nielson, 2009; Zhang et al., 2009b,2011;Karplus et al., 2011).但是，该方法仍然存在一定的局限性.

一般情况下由于受到地震资料信噪比限制，有效的S波到时数据相比P波数据数据量少、质量差，单独反演所得横波速度结构与单独反演所得纵波速度结构经常存在显著区别，由此计算得到的VP/VS图像往往会产生畸变现象，这种不合理性就不可避免地给地质解释带来巨大挑战(Eberhart-Phillips, 1990).

针对以上缺陷，Thurber和Atre (1993)首先提出了一种基于P波到时数据及S-P时间差反演VP和VS的算法，这种算法假定P波和S波具有基本相同的射线路径.显然这一假定并不适合复杂三维速度结构：此时由于P波射线路径和S波射线路径差异较大，基于反演得到的P波速度模型和S波初始速度模型并不能可靠地反演获得S波速度结构(Wagner et al., 2005)，就像由S波速度模型和P波初始速度模型并不能可靠地反演获得P波速度结构一样(Eberhart-Phillips, 1990).

为此，Conder和Wiens (2006)提出采用奇异值分解法(SVD)联合反演VP和VP/VS，该方法不受P波路径和S波路径相同这一基本假设的限制.其射线追踪过程考虑了地球的曲率，采用梯形网格，用非线性方程表征S波走时和P波慢度、VP/VS比值之间的关系，并且将平滑约束加入偏导数矩阵中，几乎利用了所有的SVD解，提高了反演的稳定性，同时有效减小了速度模型的方差.虽然该方法给出了分辨率矩阵的定量描述，但是并不能给出相对分辨率.Zhang等(2009b)提出了另一种改进方法：一是在每次迭代后检查P波和S波射线路径，在后续反演计算VP/VS时排除射线路径差异较大的走时数据；二是使用了具有共同台站事件对的差分到时数据以确定震源区VP/VS模型.该算法避免了一些假象的产生，例如直接依据VP和VS得到地壳VP/VS模型等，但在实践中仍然存在一些不足，如需要人为确定射线路径差异的门槛值，以及某些情况下反演得到的VP与VS图像差异较大，地质解释存在困难等.为了有效解决P波路径与S波路径差异问题，Wang (2014)采用同一炮检对的P波和S波走时残差，直接将VP/VS异常作为反演的模型参数，计算P波射线路径上每一小段的相对P波慢度扰动q，将q用于迭代反演沿S波射线路径上每个格点的相对S波慢度扰动r，直到每个网格点的r达到误差允许范围内.相比以往常规反演中假设P波、S波射线路径相同而言，该方法反演得到的VP/VS比值误差更小.但是，该方法的VS模型由S波数据单独反演得到，其结果在很大程度上受到S波数据质量、反演参数及迭代次数等限制.

基于此，本文提出了一种综合利用深地震测深Pg和Sg走时数据反演上地壳纵横波速度结构的相关反演成像方法(图4).与以上方法相比，该方法从假定Pg和Sg数据的走时残差来自地下同一地质异常体、且该地质体的纵横波速度异常图像相关的这一基本思想出发，在交叉迭代反演纵横波速度过程中使用纵横波速度异常的相关性作为约束，将走时残差更准确可靠地反投影到目标地质异常体，从而有效保证了成像结果可靠性.

在资料处理和解释过程中，如果测线弯曲程度不大，一般按照直线方式来处理；如果测线弯曲程度较大，最好将测线进行分段直线划分后开展后续工作.本文提出的相关反演方法是针对三维介质进行反演的，在反演开始的数据输入过程中包含了接收点的x坐标(包括测点相对原点的坐标和测点相对炮点的坐标)、y坐标、z坐标(高程)等有效信息，有效解决了检波器的高程差异问题.几种典型模型的对比试验验证了该方法的有效性和优势，作为示例本文将该方法应用于青藏高原西部真实的DSS数据，揭示了上地壳内印度-欧亚板块陆陆碰撞的证据.

2 上地壳纵横波速度结构相关反演成像方法

2.1 构建一维速度模型

利用初至走时数据进行有限差分反演成像的方法具有震相拾取简单、成像速度快、容许对模型进行精细剖分、适应性强等优点(兰海强等，2012a,2012b；赵烽帆等，2014)，因此近二十年以来该方法在重建上部地壳结构，研究构造精细特征方面得到了广泛应用(Vidale, 1990; 王椿镛等，1997，2002；段永红等，1999；Zhao et al., 2004; 赵爱华和丁志峰，2005；Zhang et al., 2009b；刘一峰和兰海强，2012；Teng et al., 2013).Hole (1992)提出了利用初至走时数据反演地壳结构的反投影重建方法.该方法基本思想是，以小的矩形网格单元划分表征模型空间,近似认为地震波在网格单元上以平面波形式传播.在此基础上通过程函方程的有限差分算子计算出地震波从震源出发到所有网格点的初至走时场，根据计算得到的走时场沿最陡下降的方向由检波器到炮点快速追踪获得射线路径.反演过程采用反投影算法，即将走时残差平均分配到射线路径上，获得各网格单元内的慢度扰动值，进而获得更新的速度模型.

反演过程中对走时扰动与慢度扰动的关系作线性化处理，即

(1)

通过迭代方法以简单反投影近似求解该方程，最终给出慢度扰动模型：

(2)

其中，δtj和lj分别是第j条射线的走时残差和总长度，gj(r)是数据核，A是以射线为中心的横截面积.假设射线路径穿过K个节点，则每个节点上的慢度扰动为：

(3)

其中，δtk和lk分别是第k个网格节点的走时残差和射线长度.

这种反演方法的计算效率允许对模型采取密集采样，提供了空间分辨率较好的三维层析图像.鉴于该方法在反演过程使用了拟线性最小二乘反演，因此选择恰当的初始模型十分重要.Kissling等(1994,1995)在解决天然源地震的精定位和一维速度模型问题时提出确定最佳一维速度模型的方法, 最小一维速度模型方法计算的速度模型可使定位结果的走时残差均方根最小, 定位精度更高.理论上最佳一维模型应逼近真实的地球模型，然而在宽角反射/折射地震资料处理解释领域，作为利用Pg或Sg走时数据反演上地壳速度结构全过程的第一步，如何构建上地壳速度结构的最佳一维模型长期以来国内外文献一直鲜有论述且未见得到应用.一种变通的方法是尝试比较若干不同一维初始模型，及不同参数进行反演，然后比较不同反演结果对应的射线覆盖和走时拟合效果(徐涛等，2004；李飞等，2013)，从中选择较为合理的速度模型作为真实地球的近似模型，这在某种程度上的确有助于逼近真实地球模型，然而初始模型选择的主观随意性大大减少了结果的可信度.

为得到可靠的一维模型，基于Hole(1992)，Ammon和Vidale (1993)的反演代码，我们在初始反演过程中将反演参数的横向网格单元和节点数设定为最大，覆盖全剖面长度，横向光滑长度也设定为剖面整个长度，从而通过10次迭代反演巧妙地得到一维速度模型.为验证该方法的有效性，我们构建了两个理论速度模型，并计算相应的初至走时数据(图1，2)，目的是为了验证针对不同的异常体(简单规则异常和复杂不规则异常)，我们的反演方法是否均行之有效.剖面长度均为380 km，地表有9炮和380个检波器，间隔为1 km.需要注意的是，每张图中的(a)和(b)分别表示扰动的P波和S波速度模型，(c)和(d)是它们的理论折合走时曲线.这两种速度模型具有相同的初始一维P波和S波背景(图1a, 1b的插图).与此同时，我们简单地假设两个模型中VP/VS比值为恒定值1.80.

图3a和3b分别表示利用10次迭代后得到的VP、VS模型最终反演得到的一维速度-深度曲线.模型1(图1a，1b)的灰色曲线和模型2(图2a，2b)的虚线基本上接近原始理论曲线(黑色曲线)，这表明即使在初始一维模型与理论一维模型差别较大的情况下，反演得到的一维速度模型仍然能很好地逼近理论模型，因此我们的方法是稳定有效的，它可以广泛地用于合理构建与Pg和Sg数据对应的一维纵横波速度模型，从而稳定有效地开展后续二维反演工作.

2.2 相关反演成像方法

尽管获得了可靠的一维速度模型，但是通过反演构建可靠的二维VP、VS和VP/VS模型仍然是一个很大的挑战.图5和图6的右图分别显示了利用Ammon和Vidale(1993)的反演代码最终得到的模型1和模型2的VP、VS和VP/VS误差.这两个图表明尽管重建的VP与VS模型显示了与初始模型相似的特性，但是VS模型已经在很大程度上发生了扭曲，最终得到的VP/VS模型已经严重偏离了给定的初始值1.80.这些结果清楚地表明，分别单独反演VP和VS不能可靠地推断出最终的VP/VS模型，类似地，不能简单地依据VP和VP/VS的值推断VS(Eberhart-Phillips,1990；Wagner et al., 2005；Conder and Wiens,2006; Wang, 2014).

为构建可靠的二维模型，尝试在反演过程中添加紧约束：纵横波速度反演交叉进行，且以每次迭代后生成的VP或VS模型以及VP/VS比值作为下一步反演的VS或VP初始速度模型.反演方法见图4，“相关反演”的概念描述如下：

假定Pg和Sg数据的走时残差完全来自地下同一地质异常体，且该地质体的纵横波速度相对于围岩同步地表现出速度高或速度低的特点，则该地质体附近的VP和VS属性的图像特征存在明显的正相关或负相关.据此每次迭代后计算VP和VS模型的相关系数，分析每个网格点存在真实地质体的概率.三维空间围绕一个特定格点共有7个相邻网格点(包括其本身)，因而每次迭代后网格点的系数被定义为：

图1 模型1：具有矩形异常的上地壳速度模型((a) P波，(b)S波)以及对应的理论折合走时曲线((c)，(d))Fig.1 Model 1: upper crust velocity models with rectangular anomaly ((a) for P-wave and (b) for S-wave),and their corresponding theoretical reduced traveltime curve ((c) and (d))

图2 模型2：正弦速度模型((a) P波，(b)S波)以及对应的理论折合走时曲线((c)，(d))加入的异常延伸至上地壳底部，该模型与图1具有相同的一维速度背景.插图：(a)和(b)右下角的说明分别表示初始一维VP或VS模型.Fig.2 Model 2: sine-type velocity models ((a) for P-wave and (b) for S-wave) and their corresponding theoretical reduced traveltime curve ((c) and (d)) The added anomalies extend downward to the bottom of the upper crust, and this model has the same 1-D velocity background as that of Fig.1.Inset: illustrations in (a) and (b) at the bottom right corner respectively denotes the initial 1-D VP or VS model.

图3 反演获得的P波和S波一维速度曲线黑色实线：初始理论一维深度-速度曲线；点划线：用于反演得到更加可靠的一维曲线的初始一维曲线，其适合于接下来的二维反演；灰色实线：模型1 10次迭代后的反演结果；虚线：模型2 10次迭代后的反演结果.Fig.3 Derived 1-D velocity curve of P-wave (a) and S-wave (b) using inversion methodBlack solid line: original theoretical 1-D depth-velocity curve; Dash dot line: initial 1-D curve for inversion to get a more reliable 1-D curve which is appropriate for following 2-D inversion; Gray solid line: inversion results after 10 times iteration for model 1; Dotted line: inversion results after 10 times iteration for model 2.

×(|ai|+|bi|)0.125,

(4)

其中，

ai=(VPnew-VPori)/VPori，

bi=(VSnew-VSori)/VSori.

事实上，一般固定C0为0.00001，VPnew和VSnew表示每次迭代后所考虑点的新速度值，VPori和VSori是初始一维速度模型的速度值.(4)式表示围绕特定格点平均的相关系数.用这个公式，经过一次迭代后，如果邻近点的VP和VS没有发生变化，即ai=bi=0.0，系数将被保持在最小值C0；但是比如ai=bi=0.2时，我们将得到一个较大的系数，其值等于0.8917+C0.这意味着，如果每一步反演后围绕一个网格点存在相似的速度变化，那么计算出的相关系数将是一个较大的值.另一方面可以使用上述相关系数来描述现有的地质异常，计算的相关系数越大，对应地质体存在的概率越高.基于网格点的系数Cor1，可以很容易地通过线性插值计算出一个单元格的相关系数Cor2.然后利用这样一种“相关反演”的想法对Hole (1992)的反演过程进行了如下修改：

假设单元格的慢度扰动与计算得到的相关系数Cor2成比例，即模型空间某一格点周围纵横波速度图像特征的相关性越大，相应地质异常存在的可能性越大，该潜在异常体越可能导致较大的慢度扰动.上述可以被描述为：

(5)

(5)式假定射线路径穿过n个单元格，单元格k的长度为lk，Δtj是第j条射线的走时残差，Δui和Cori分别是第i个单元格的慢度扰动和相关系数.(5)式表示针对第j条射线的走时残差Δtj，根据第i个单元格的相关系数Cori，更新第i个单元格的慢度扰动Δui.如此将相关系数成功引入到慢度扰动计算中，进而更新速度模型.如果相关系数Cori恒等于1.0，该种情况下，(5)式可简化为Hole (1992)的常规表达式，所计算的慢度扰动变为：

(6)

其中L是射线总长度.

由以上分析可知，Hole(1992)有限差分反演的基本思想是，射线路径上每个网格单元均可能引起走时残差，且它们的几率是相同的.VP与VS模型单独反演，彼此之间没有任何约束.与此相比，我们在反演过程中通过引入P波和S波的相关系数Cor，将P波速度与S波速度紧密结合起来：利用每次迭代后的VP反演结果来约束VS反演过程，类似地，用VS的反演结果来约束VP反演，然后利用VP和VS模型之间的相关性来指导反演过程，从而得到更好的反演结果.

相应地我们改进了Hole(1992)，Ammon和Vidale (1993)的反演代码，得到改进的慢度扰动.射线追踪和慢度扰动计算后，取每个网格点慢度平均值，速度模型和VP/VS比值模型将被更新.图5和图6的右图分别表示使用该方法最终反演得到模型1和2的VP、VS和VP/VS比值误差的模型.单元格或网格点的反演参数重采样尺寸和平滑尺寸保持最小.这两张图表明，与依据Hole (1992)等算法反演得到的结果相比，由我们改进方法反演得到的VP和VS模型有明显改善，显然更接近初始理论模型(图1和图2)，而依据Hole (1992)算法分别反演纵横波速度结构，进而计算得到的速度比模型却严重失真；依据我们方法反演得到的VP/VS比值误差最大值仅为0.02～0.03，相比之下依据传统方法算得的VP/VS比值误差最大值高达0.2.以上试验表明，采用不同反演方法得到的最终模型可能具有相同水准的走时拟合效果，但这并不能保证反演结果的可靠性.即使在速度模型横向变化较大情形下，我们的二维反演方法对于不同类型的合成走时数据集也取得满意效果(图7).

图4 利用Pg和Sg数据获取上地壳结构的改进二维反演流程图Fig.4 Flow-chart of improved 2-D inversion procedure for upper crust structure using Pg and Sg data

图5 利用合成Pg和Sg走时数据反演获得的速度结构(模型1)(a) 利用Ammon和Vidale(1993)反演代码得到的上地壳VP结构； (b) 利用我们的改进方法得到的上地壳VP结构； (c) 利用Ammon和Vidale(1993)反演代码得到的上地壳VS结构； (d) 利用我们的改进方法得到的上地壳VS结构； (e) 利用Ammon和Vidale(1993)反演代码得到的VP/VS误差； (f) 利用我们的改进方法得到的VP/VS误差.所有反演结果均经过10次迭代获得.Fig.5 Final 2-D inversion results for Model 1 using synthetic Pg and Sg traveltime data(a) Upper crust VP structure derived by the inversion code of Ammon and Vidale (1993); (b) Upper crust VP structure derived by our improved approach; (c) Upper crust VS structure derived by the inversion code of Ammon and Vidale (1993); (d) Upper crust VS structure derived by our improved approach; (e) VP/VS error derived by the inversion code of Ammon and Vidale (1993); (f) VP/VS error derived by our improved approach. All the inversion results were acquired after 10 times iteration.

图6 利用合成Pg和Sg走时数据反演获得的速度结构(模型2)(a) 利用Ammon和Vidale(1993)反演代码得到的上地壳VP结构； (b) 利用我们的改进方法得到的上地壳VP结构； (c) 利用Ammon和Vidale(1993)反演代码得到的上地壳VS结构； (d) 利用我们的改进方法得到的上地壳VS结构； (e) 利用Ammon和Vidale(1993)反演代码得到的VP/VS误差； (f) 利用我们的改进方法得到的VP/VS误差.所有反演结果均经过10次迭代获得.Fig.6 Final 2-D inversion results for Model 2 using synthetic Pg and Sg traveltime data(a) Upper crust VP structure derived by the inversion code of Ammon and Vidale (1993); (b) Upper crust VP structure derived by our improved approach; (c) Upper crust VS structure derived by the inversion code of Ammon and Vidale (1993); (d) Upper crust VS structure derived by our improved approach; (e) VP/VS error derived by the inversion code of Ammon and Vidale (1993); (f) VP/VS error derived by our improved approach. All the inversion results were acquired after 10 times iteration.

图7 二维反演结果走时残差均方根(a) 对应于图1和图2中的VP模型; (b) 对应于图1和图2中的VS模型.Fig.7 Root Mean Square (RMS) of time residual of the 2-D inversion results(a) For the VP models in Figs.1 and 2; (b) For the VS models in Figs.1 and 2.

2.3 模型的定量评估

(7)

对模型空间每一格点扫描计算速度扰动引起的目标函数(走时残差均方根)变化，得到的残差均方根变化值ΔF定量描述了模型空间某点速度变化的重要性，或观测系统对模型空间某一点速度变化的分辨能力：图8和图9同时也分别展示了模型1和2的射线覆盖(图8a和8b，图9a和9b)和以上扰动函数ΔF值.对比二图可见，速度的可靠性在很大程度上取决于射线覆盖，然而在某些情况下，密集射线覆盖并不总意味着反演得到的速度变化特征的可靠性，一般情形下采用扰动函数ΔF的方法可以更简单直观地描述模型的可靠性：某格点附近扰动函数ΔF值越大，该点附近反演得到的速度图像特征越可靠(图8c和8d，图9c和9d).

尽管如此，即使采用上述评估方法，我们仍然不知道是否能够正确揭示潜在的地质体.幸运的是，由于我们采用了相关反演的办法，反演所得纵横波速度图像特征的相关性在很大程度上指示了目标地质体的存在与否：图8和9清楚显示，采用相关反演所得纵横波速度结构的相关性图像能够准确揭示出速度异常体在剖面中的位置和大致形态，这是其他评估办法无法做到的，由此也充分体现了我们采用反演方法的优越性.因此，计算VP和VS图像特征相关性的方法可以帮助我们揭示VP与VS模型可能的同步速度异常，从而进一步验证反演结果的可靠性.

3 札达—泉水沟深地震测深剖面上地壳速度结构

为了验证我们的反演成像方法及解评价方法的适用性和优越性，我们尝试基于以上方法利用实际深地震测深(DSS)数据反演青藏高原西部札达—泉水沟剖面下方的上地壳速度结构.

印度板块向欧亚大陆的俯冲碰撞对区域地形和环境带来了显著影响，如一千多公里的地壳缩短，西藏地壳的东向挤压以及中亚和喜马拉雅周边许多著名的活动高山.青藏高原西缘浓缩了青藏高原古生代和新生代的造山历史，作为研究陆陆碰撞的理想窗口正日益成为国内外地学领域关注的焦点(Murphy et al., 2000; Lacassin et al., 2004; Wittlinger et al., 2004; Valli et al., 2007; 李海兵等，2007,2008; Caldwell et al., 2009;许志琴等，2011).研究区位于西构造结东侧，处于喜马拉雅、拉萨、羌塘、塔里木、拉达克、甜水海及西昆仑等不同地体叠合衔接部位，且发育巨大的喀喇昆仑走滑断层.为探讨研究区新生代以来的深部造山过程及其对环境的影响，中国科学院地质与地球物理研究所从2011年9月至11月在该区域开展了380 km长的深地震测深(DSS)项目(图10).南北走向的深地震测深(DSS)剖面始于札达盆地，然后穿越阿伊拉日居山、喀喇昆仑断裂(KF)、拉萨地块、班公湖—怒江缝合带(BNS)、羌塘地块、郭扎—龙木错断裂，并最终延伸到北端的甜水海地块.本文基于前文所述方法利用获得的DSS数据揭示青藏高原西部与俯冲碰撞在剖面上地壳留下的构造遗迹.

图8 模型1二维反演结果评估(a) 所得VP模型的射线覆盖； (b) 所得VS模型的射线覆盖； (c) VP模型每一格点速度扰动引起的目标函数(走时残差均方根)变化值； (d) VS模型每一格点速度扰动引起的目标函数(走时残差均方根)变化值； (e) 计算得到的VP与VS模型之间的相关系数.Fig.8 Assessment of the 2-D inversion results as for Model 1 using our tomographic system(a) Ray-coverage of derived VP model; (b) Ray-coverage of derived VS model; (c) RMS variation on VP model that corresponds to the velocity perturbation at each point; (d) RMS variation on VS model that corresponds to the velocity perturbation at each point; (e) Computed correlation coefficient between VP and VS model.

图9 模型2二维反演结果的评估(a) 所得VP模型的射线覆盖； (b) 所得VS模型的射线覆盖； (c) VP模型每一格点速度扰动引起的目标函数(走时残差均方根)变化； (d) VS模型每一格点速度扰动引起的目标函数(走时残差均方根)变化； (e) 计算得到的VP与VS模型之间的相关系数.Fig.9 Assessment of the 2-D inversion results as for Model 2 using our tomographic system(a) Ray-coverage of derived VP model; (b) Ray-coverage of derived VS model; (c) RMS variation on VP model that corresponds to the velocity perturbation at each point; (d) RMS variation on VS model that corresponds to the velocity perturbation at each point; (e) Computed correlation coefficient between VP and VS model.

图10 青藏高原西部札达—泉水沟深地震测深剖面及其构造背景主要断裂：MBT，主边界冲断裂；MCT，主中央冲断裂；STDS，藏南拆离系；IYS，印度—雅鲁藏布江缝合带；KKF，喀喇昆仑断裂；SF，狮泉河断裂；BNS，班公—怒江缝合带；DT，多玛冲断裂；LCF，龙木错断裂. 红色五角星表示炮点，蓝色三角形表示接收点.Fig.10 Tectonic setting and the deep seismic sounding profile in Western Tibetan Plateau which was carried out in September and October, 2011Main Faults: MBT, Main Boundary Thrust; MCT, Main Central/Panjal Thrust; STDS: Southern Tibetan Detachment System; IYS, Indus-Yarlung Zangbo River Suture Belt; KKF, Karakorum Fault; SF, Shiquanhe Fault; BNS, Bangong-Nujiang Suture Belt; DT: Domar Thrust; LCF, Longmu Co Fault. The red star denotes the shot point, and the blue triangle shows the receiver.

该DSS测线海拔高程在4.0～5.8 km，大部分检波器被部署在海拔4.6～5.1 km的高度.在不同构造单元(图10)总共触发了七炮，工作区每一炮点深度～40 m，TNT炸药量3000 kg.共有200台便携式三分量数字地震仪沿剖面以2.0 km左右的间距接收地震数据.这些地震仪记录的每炮信号长度为5 min，采样间隔5 ms.对P波和S波进行带通滤波保留1～10 Hz的频率成分.图11和图12分别展示了折合速度为6.0 km·s-1和3.5 km·s-1的P波和S波的所有折合地震记录.可以利用走时互换原理从这些记录中拾取出清晰的P波和S波初至(Pg和Sg)，一般来说Pg的信噪比(SNR)要比Sg的高得多.

Pg震相一般可以在偏移距0～100 km观测到，某些炮可观测到最大偏移距达120 km的地方，如第一炮(SP1)的北支和第五炮(SP5)的南支(图11a，11e).在SP2的南支(图11b)偏移距-50～30 km，SP3的北支(图11c)10～50 km，SP6的北支(图11f)5～40 km，Pg走时发生明显延迟.这些特征也可以由拾取的Sg数据(图12)观察到.S波记录显示在以下偏移距范围内Sg的到时数据也发生了明显延迟：SP2的南支(图12b)-30～5 km，SP3(图12c)的-90～-60 km和10～50 km，SP5的南支(图12e)-80～-60 km，SP6的北支(图12f)5～40 km.这些延迟的Pg和Sg走时数据彼此相关性明显，炮点SP3和SP6等记录剖面上纵横波初至走时同步滞后表明这两炮点附近存在对应的沉积盆地.

从记录剖面上共拾取了398个Pg和382个Sg走时数据用于反演.正演过程将模型划分为1 km×1 km的矩形单元.首先，利用前文改进算法分别反演纵横波一维速度结构，10次迭代后得到用于二维反演的一维初始速度模型(图13).然后，使用前文所述相关反演成像方法交叉反演，10次迭代后得到稳定的纵横波速度模型和VP/VS比值模型(图15).Pg数据走时残差均方根从初次迭代的0.336 s下降到10次迭代后的0.119 s，Sg数据走时残差均方根从初次迭代的0.442 s下降到10次迭代后的0.169 s.与最终速度模型相应的计算走时数据也被标注在地震记录剖面上(图11，图12)，可以看出我们的合成走时数据与观测数据有良好的匹配效果.

图11 札达—泉水沟宽角地震剖面P波地震记录记录剖面折合速度为6.0 km·s-1，十字表示拾取的Pg数据，曲线表示利用我们的层析成像系统计算的对应最终VP模型的折合走时数据.Fig.11 P-wave seismic record of Zhada-Quanshuigou wide-angle seismic profileThe record sections were plotted with a reduced velocity of 6.0 km·s-1, and the cross denotes the picked Pg data, while the curve represents the computed reduced traveltime data correspond to our final VP model using our tomographic system.

图12 札达—泉水沟宽角地震剖面S波地震记录记录剖面折合速度为3.5 km·s-1，十字表示拾取的Sg数据，曲线表示利用我们的层析成像系统计算的对应最终VS模型的折合走时数据.Fig.12 S-wave seismic record of Zhada-Quanshuigou wide-angle seismic profileThe record sections were plotted with a reduced velocity of 3.5 km·s-1, and the cross denotes the picked Sg data, while the curve represents the computed reduced traveltime data correspond to our final VS model using our tomographic system.

图13 用于札达—泉水沟宽角地震剖面Pg和Sg数据反演的一维速度模型Fig.13 The 1-D velocity models for inversion as for the Pg and Sg data of Zhada-Quanshuigou wide-angle seismic profile

图14和图15分别展示了使用Hole(1992)，Ammon和Vidale(1993)的反演代码以及我们的相关反演成像方法所获得的VP、VS和VP/VS比值模型.可以看出，基于Hole(1992)，Ammon和Vidale(1993)的反演方法分别反演纵横波速度结构的过程，由于缺少相互关联约束，反演得到的VP和VS模型之间图像特征存在较大区别(图14a，14b)，互相关性较差，这无疑会使得VP/VS比值产生偏差或扭曲.与之对比用我们的相关反演成像方法得到的VP和VS模型(图15a，15b)相互关联的非常好，图像特征较相似，证明我们采用新反演方法得到的上地壳结构更为合理.图16显示的射线路径和走时残差均方根扰动表明我们模型大多数地方的速度变化是可靠的.图16e表示计算得到的VP与VS模型之间的相关性，其表明在剖面的某些位置存在明显的正异常，如横向坐标50 km，深度10 km的位置，以及横向坐标120～140 km的位置，存在着明显的地质异常.

图15a、15b显示，该剖面最南端札达盆地下方，及狮泉河盆地下方，均存在明显的低VP、低VS及低VP/VS异常，且异常向下延伸的深度达到8～10 km.这可能分别与印度洋板块的北向俯冲及古班公—怒江洋盆的闭合有关，其北向俯冲向下延伸的几何形态表明它可能是不同地质历史时期板块阶段俯冲碰撞所遗留的构造痕迹或“化石”.推测其基本原理为，伴随俯冲碰撞前陆盆地堆积大量海洋沉积覆盖.在山脉隆起地形抬升过程中，水体被逐渐排泄出高原，从而使得岩土较为干燥.岩石在干燥条件下横波速度伴随压力下降减小的幅度较慢，而纵波速度伴随压力下降减小的幅度较快，这一岩石物理基础可以很好地解释上述纵横波速度比偏低的现象.

Kapp等(2003)的研究表明，狮泉河逆冲断裂以南的狮泉河盆地覆盖有大量中生代沉积，且狮泉河北向俯冲向下延深到10 km以下，这与我们的以上分析是一致的.

阿伊拉日居山的隆升剥蚀、岩浆活动与喀喇昆仑断裂的活动密切相关.对喀喇昆仑东南段阿伊拉日居山地区同构造片麻岩-花岗岩的年代学研究认为喀喇昆仑断裂形成时代为23～25 Ma,持续变形到约12 Ma,之后阿伊拉日居山快速隆升(李海兵等，2006,2007,2008).王世锋等(2008)对札达盆地内新生代沉积给出了全面而系统的磁性地层学的时代约束,并探讨了盆地的控盆构造问题, 揭示盆地内新生代地层的古地磁年龄为9.5～2.6 Ma，且控盆构造为喀喇昆仑断裂.我们的成像结果显示，与札达盆地及狮泉河盆地8～10 km的厚度相比较，噶尔盆地沉积厚度仅3～4 km，这可能与盆地形成机制不同有关：噶尔盆地为走滑断陷盆地，而前两者为板块俯冲碰撞的前陆盆地.新生代以来伴随喀喇昆仑断裂在不同阶段的走滑拉张运动，可能同步产生了阿伊拉日居山逐步抬升和噶尔盆地逐步沉降、断层切割越来越深等地质现象.

日松背斜及炮点SP6附近存在明显高速异常体.炮点SP4南侧的异常体具有较高的VP、VS和VP/VS比值，上述性质表明近地表岩石主要由基性或超基性材料构成，这也由部分暴露于地表的超基性铁镁质岩石得到验证.

我们的VP与VS模型还显示，从SP6到SP7的剖面下方存在另一个浅盆地，其南缘可能被一个相对较深的断裂(深度大约6 km)所切割.除此之外，在著名的班公湖—怒江缝合带北侧还有可能存在着最大深度达8 km的浅断裂.与此同时，尽管龙木错断裂被认为是羌塘地体和甜水海地体之间的主要边界断裂，我们的VP与VS模型并没有显示出与之相应的显著速度变化，这可能意味着它新生代的活动并不活跃.

图14 利用Ammon和Vidale (1993)反演代码得到的札达—泉水沟剖面下的上地壳VP、VS和VP/VS模型对比(a) VP速度模型； (b) VS速度模型； (c) VP/VS模型.Fig.14 The upper crustal VP, VS and VP/VS model beneath the Zhada-Quanshuigou for comparison which were derived using the inversion codes of Ammon and Vidale (1993)(a) The derived VP model; (b) The derived VS model; (c) The resulted VP/VS model.

图15 利用我们反演方法重建的札达—泉水沟宽角地震剖面下方的上地壳结构(a) VP速度模型； (b) VS速度模型； (c) VP/VS模型.图顶说明表示沿剖面的地势与构造背景.Fig.15 The upper crustal structure beneath the Zhada-Quanshuigou wide-angle seismic profile, which were inverted using our tomographic system(a) The derived VP model; (b) The derived VS model; (c) The resulted VP/VS model.The illustration on the top in this figure shows the relief and the tectonic setting along the profile.

图16 札达—泉水沟剖面下方上地壳速度模型二维反演结果评估(a) VP模型对应的射线覆盖； (b) VS模型对应的射线覆盖； (c) VP模型格点速度扰动引起的走时残差均方根变化值； (d) VS模型格点速度扰动引起的走时残差均方根变化值； (e) 计算得到VP与VS模型之间的相关系数.Fig.16 Assessment of the 2-D inversion results as for our upper crustal velocity model beneath Zhada-Quanshuigou profile(a) Ray-coverage of derived VP model; (b) Ray-coverage of derived VS model; (c) RMS variation on VP model that corresponds to the velocity perturbation at each point; (d) RMS variation on VS model that corresponds to the velocity perturbation at each point；(e) The computed correlation coefficient between VP and VS model.

4 结论

当前基于深地震测深纵横波初至走时数据的层析成像结果越来越广泛地被应用于揭示造山带、盆地等不同类型构造域下方的上地壳结构精细特征.本文基于纵横波初至走时数据异常来源于同一地质体，因而该地质体附近纵横波速度图像特征密切相关这一思想，在Hole(1992)及Ammon和Vidale(1993)工作基础上提出了反演上地壳纵横波速度结构的相关反演成像方法，这一方法可描述为：

(1)纵横波速度反演交替进行，在迭代反演过程中每通过一次反演获得相应的纵波速度(或横波速度)结构后，更新相应的纵横波速度比模型以及相应的横波(或纵波)速度反演的初始模型，然后开展下一步横波(或纵波)速度反演工作.

(2)在反演过程中依据纵横波速度结构的相关性信息和射线路径长度将走时残差以不同权重分配到射线路径经过的单元，依据网格节点周围平均的慢度扰动更新速度模型.

(3)正反演过程分别基于Vidale (1988)的有限差分走时计算方法和Hole (1992)的反投影成像方法.两种典型模型试验表明，该方法应用于上地壳速度结构反演成像过程，可有效提高反演结果的可靠性，在很大程度上避免了常规单独反演纵波和横波速度过程容易带来的畸变和失真.

此外我们发展了获取初始一维速度模型及直接定量描述每个网格点速度值可靠性的解评价方法：通过计算每个网格点参数单独变化引起的目标函数——走时残差均方根的改变量来评估观测系统对特定格点速度变化的分辨能力.这种方法相比射线覆盖数评价方法更为简洁直观.相关反演最终结果所显示纵横波速度图像特征之间的相关性，则可以直接用于估计地质异常体存在与否及可能的位置.

以上相关反演成像方法被成功地应用于重建青藏高原西部札达—泉水沟深地震测深(DSS)剖面下方的上地壳速度结构，揭示出与青藏高原西缘板块碰撞相关的上地壳速度结构特征：札达盆地及狮泉河盆地下方存在明显的VP和VS低速异常及低VP/VS异常，可能是印度板块北向俯冲及班公湖—怒江特提斯洋闭合在上地壳遗留的地球物理“化石”和证据.噶尔盆地沉积厚度较小，其形成与新生代以来喀喇昆仑断裂的走滑拉张运动有关.

致谢 感谢J.A.Hole教授、C.J.Ammon教授和J.E.Vidale教授提供用于该项研究的反演软件，并衷心感谢来自田小波研究员、徐涛副研究员和梁晓峰副研究员的宝贵意见和建议.谨以此文纪念中国科学院地质与地球物理研究所张忠杰研究员(1964—2013).

Ammon C J, Vidale J E. 1993. Tomography without rays.Bull.Seism.Soc.Am., 83(2): 509-528.

Bai Z M, Zhang Z J, Wang Y H. 2007. Crustal structure across the Dabie-Sulu orogenic belt revealed by seismic velocity profiles.J.Geophys.Eng., 4(4): 436-442.

Bregman N D, Bailey R C, Chapman C H. 1989. Crosshole seismic tomography.Geophysics, 54(2): 200-215, doi: 10.1190/1.1442644.

Caldwell W B, Klemperer S L, Rai S S, et al. 2009. Partial melt in the upper-middle crust of the northwest Himalaya revealed by Rayleigh wave dispersion.Tectonophysics, 477(1-2): 58-65.

Conder J A, Wiens D A. 2006. Seismic structure beneath the Tonga arc and Lau back-arc basin determined from jointVP,VP/VStomography.Geochem.Geophys.Geosyst., 7(3): Q03018, doi: 10.1029/2005GC001113.

Ding Z F, He Z Q, Sun W G, et al. 1999. 3-D crust and upper mantle velocity structure in eastern Tibetan plateau and its surrounding areas.ChineseJ.Geophys. (in Chinese), 42(2): 197-205.Duan Y H, Lai X L, Zhang X K, et al. 1999. Finite-difference 2-D and 3-D seismic traveltime velocity tomography.NorthChinaEarthquakeSciences(in Chinese), 17(4): 53-60.

Eberhart-Phillips D. 1990. Three-dimensional P and S velocity structure in the Coalinga Region, California.J.Geophys.Res.:SolidEarth, 95(B10): 15343-15363.

Hole J A. 1992. Nonlinear high-resolution three-dimensional seismic travel time tomography.J.Geophys.Res., 97(B5): 6553-6562.

Hole J A, Zelt B C.1995. 3-D finite-difference reflection traveltimes.Geophys.J.Int., 121(2): 427-434.

Kapp P, Murphy M A, Yin A, et al. 2003. Mesozoic and Cenozoic tectonic evolution of the Shiquanhe area of western Tibet.Tectonics, 22(4), doi: 10.1029/2001TC001332.

Karplus M S, Zhao W, Klemperer S L, et al. 2011. Injection of Tibetan crust beneath the south Qaidam Basin: Evidence from INDEPTH IV wide-angle seismic data.J.Geophys.Res., 116(B7): B07301, doi: 10.1029/2010JB007911.

Kissling E, Ellsworth W L, Eberhart-Phillips D E, et al. 1994. Initial reference models in local earthquake tomography.J.Geophys.Res.:SolidEarth, 99(B10): 19635-19646.

Kissling E, Solarino S, Cattaneo M. 1995. Improved seismic velocity reference model from local earthquake data in Northwestern Italy.TerraNova, 7(5): 528-534.

Lacassin R, Valli F, Arnaud N, et al. 2004. Large-scale geometry, offset and kinematic evolution of the Karakorum fault, Tibet.EarthPlanet.Sci.Lett., 219(3-4): 255-269.

Lan H Q, Zhang Z, Xu T, et al. 2012a. A comparative study on the fast marching and fast sweeping methods in the calculation of first-arrival traveltime field.ProgressinGeophys. (in Chinese), 27(5): 1863-1870, doi: 10.6038/j.issn.1004-2903.2012.02.005.Lan H Q, Zhang Z, Xu T, et al. 2012b. Effects due to the anisotropic stretching of the surface-fitting grid on the traveltime computation for irregular surface by the coordinate transforming method.ChineseJ.Geophys. (in Chinese), 55(10): 3355-3369, doi: 10.6038/j.issn.0001-5733.2012.10.018.

Li F, Xu T, Wu Z B, et al. 2013. Segmentally iterative ray tracing in 3-D heterogeneous geological models.ChineseJ.Geophys. (in Chinese), 56(10): 3514-3522, doi: 10.6038/cjg20131026.Li H B, Valli F, Xu Z Q, et al. 2006. Deformation and tectonic evolution of the Karakorum fault, western Tibet.GeologyinChina(in Chinese), 33(2): 239-255.

Li H B, Valli F, Liu D Y, et al. 2007. The formation age of the Karakorum Fault in western Tibet: Constraints from SHRIMP U-Pb dating of zircons.ChineseScienceBulletin, 52(4):438-447.Li H B, Valli F, Arnaud N, et al. 2008. Rapid uplifting in the process of strike-slip along the Karakorum fault zone in western Tibet: Evidence from40Ar/39Ar thermochronology.ActaPetrologicaSinica(in Chinese), 24(7): 1552-1556.

Liu Y F, Lan H Q. 2012. Study on the numerical solutions of the eikonal equation in curvilinear coordinate system.ChineseJ.Geophys. (in Chinese), 55(6): 2014-2026, doi: 10.6038/j.issn.0001-5733.2012.06.023.

Murphy M A, Yin A, Kapp P, et al. 2000. Southward propagation of the Karakoram fault system, southwest Tibet: Timing and magnitude of slip.Geology, 28(5): 451-454.

Saito H. 1989. Traveltimes and raypaths of first-arrival seismic waves: Computation method based on Huygens′ Principle. ∥ 59thAnn. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstract, 244-247.Tarantola A, Valette B. 1982. Generalized nonlinear inverse problems solved using the least squares criterion.Rev.Geophys.SpacePhys., 20(2): 219-232.Teng J W, Zhang Z J, Zhang X K, et al. 2013. Investigation of the Moho discontinuity beneath the Chinese mainland using deep seismic sounding profiles.Tectonophysics, 609: 202-216, doi: 10.1016/j.tecto.2012.11.024.

Thurber C H. 1983. Earthquake locations and three-dimensional crustal structure in the Coyote Lake area, central California.J.Geophys.Res., 88(B10): 8226-8236.

Thurber C H, Atre S R. 1993. Three-dimensionalVP/VSvariations along the Loma Prieta rupture zone.Bull.Seism.Soc.Am., 83(3): 717-736.

Thybo H, Nielsen C A. 2009. Magma-compensated crustal thinning in continental rift zones.Nature, 457(7231): 873-876, doi: 10.1038/nature07688.

Valli F, Arnaud N, Leloup P H, et al. 2007. Twenty million years of continuous deformation along the Karakorum fault, western Tibet: a thermochronological analysis.Tectonics, 26(4), doi: 10.1029/2005TC001913.

Vidale J E. 1988. Finite-difference traveltime calculation.Bull.Seism.Soc.Am., 78: 2062-2076.

Vidale J E. 1990. Finite-difference calculation of traveltimes in three dimensions.Geophysics, 55(5): 521-526.

Wagner L S, Beck S, Zandt G. 2005. Upper mantle structure in the south central Chilean subduction zone (30°to 36°S).J.Geophys.Res.:SolidEarth, 110(B1): B01308, doi: 10.1029/2004JB003238.

Wang C Y, Zhang X K, Ding Z F, et al. 1997. Finite-difference tomography of upper crustal structure in Dabieshan orogenic belt.ChineseJ.Geophys. (ActaGeophysicaSinica) (in Chinese), 40(4): 495-502.

Wang C Y, Mooney W D, Wang X L, et al. 2002. Study on 3-D velocity structure of crust and upper mantle in Sichuan-Yunnan region, China.ActaSeismologicaSinica(in Chinese), 24(1): 1-16, doi: 10.3321/j.issn:0253-3782.2002.01.001.Wang F Y, Duan Y H, Yang Z X, et al. 2008. Velocity structure and active fault of Yanyuan-Mabian seismic zone—The result of high-resolution seismic refraction experiment.ScienceinChinaSeriesD:EarthSciences, 51(9): 1284-1296.

Wang S F, Zhang W L, Fang X M, et al. 2008. Magnetic & stratigraphic characteristics and tectonic significance of Zhada Basin in Southwest Tibet.ScienceinChina(in Chinese), 53(6): 676-683.

Wang Z. 2014. Joint inversion of P-wave velocity andVP-VSratio: imaging the deep structure in NE Japan.AppliedGeophysics, 11(2): 119-127.Wittlinger G, Vergne J, Tapponnier P, et al. 2004. Teleseismic imaging of subducting lithosphere and Moho offsets beneath western Tibet.EarthPlanet.Sci.Lett., 221(1-4): 117-130.

Xu T, Xu G M, Gao E G, et al. 2004. Block modeling and shooting ray tracing in complex 3D media.ChineseJ.Geophys. (in Chinese), 47(6): 1118-1126, doi: 10.3321/j.issn.0001-5733. 2004.06.027.

Xu T, Xu G M, Gao E G, et al. 2006. Block modeling and segmentally iterative ray tracing in complex 3D media.Geophysics, 71(3): T41-T51.

Xu T, Zhang Z, Gao E, et al. 2010. Segmentally iterative ray tracing in complex 2D and 3D heterogeneous block models.Bull.Seismol.Soc.Am., 100(2): 841-850.

Xu T, Li F, Wu Z B, et al. 2014. A successive three-point perturbation method for fast ray tracing in complex 2D and 3D geological models.Tectonophysics, 627: 72-81, doi: 10.1016/j.tecto.2014.02.012.Xu Z Q, Li H B, Tang Z M, et al. 2011. The transformation of the terrain structures of the Tibet Plateau through large-scale strike-slip faults.ActaPetrologicaSinica(in Chinese), 27(11): 3157-3170.

Zelt C A. 1999. Modelling strategies and model assessment for wide-angle seismic traveltime data.Geophys.J.Int., 139(1): 183-204.Zhang H J, Thurber C, Bedrosian P. 2009a. Joint inversion forVP,VS, andVP/VSat SAFOD, Parkfield, California.Geochem.Geophys.Geosyst., 10(11), doi: 10.1029/2009GC002709.

Zhang J, Toksozf M N. 1998. Nonlinear refraction traveltime tomography.Geophysics, 63(5): 1726-1737.

Zhang X K, Wang C Y, Liu G D, et al. 1996. Fine crustal structure in Yanqing-Huailai region by deep seismic reflection profiling.ChineseJ.Geophys. (ActaGeophysicaSinica) (in Chinese), 39(3): 356-364.

Zhang Z J, Bai Z M, Mooney W, et al. 2009b. Crustal structure across the Three Gorges area of the Yangtze platform, central China, from seismic refraction/wide-angle reflection data.Tectonophysics, 475(3-4): 423-437, doi: 10.1016/j.tecto.2009.05.022.

Zhang Z J, Klemperer S, Bai Z M, et al. 2011. Crustal structure of the Paleozoic Kunlun orogeny from an active-source seismic profile between Moba and Guide in east Tibet, China.GondwanaResearch, 19(4): 994-1007.

Zhao A H, Zhang Z J, Teng J W. 2004. Minimum travel time tree algorithm for seismic ray tracing: Improvement in effciency.J.Geophys.Eng., 1(4): 245-251, doi:10.1088/1742-2132/1/4/001.

Zhao A H, Ding Z F. 2005. A double-grid algrithm for calculating traveltimes of wide-angle reflection waves.ChineseJ.Geophys. (in Chinese), 48(5) : 1141-1147.

Zhao F F, Ma T, Xu T. 2014. A review of the travel-time calculation methods of seismic first break.ProgressinGeophys. (in Chinese), 29(3): 1102-1113, doi: 10.6038/pg20140313.

附中文参考文献

丁志峰, 何正勤, 孙为国等. 1999. 青藏高原东部及其边缘地区的地壳上地幔三维速度结构. 地球物理学报, 42(2): 197-205.

段永红, 赖晓玲, 张先康等. 1999. 有限差分二、三维速度层析成象. 华北地震科学, 17(4): 53-60.

兰海强, 张智, 徐涛等. 2012a. 地震波走时场模拟的快速推进法和快速扫描法比较研究.地球物理学进展, 27(5): 1863-1870, doi: 10.6038/j.issn.1004-2903.2012.02.005.

兰海强, 张智, 徐涛等. 2012b. 贴体网格各向异性对坐标变换法求解起伏地表下地震初至波走时的影响.地球物理学报, 55(10):3355-3369, doi: 10.6038/j.issn.0001-5733.2012.10.018.

李飞, 徐涛, 武振波等. 2013. 三维非均匀地质模型中的逐段迭代射线追踪. 地球物理学报, 56(10): 3514-3522, doi: 10.6038/cjg20131026.

李海兵, Valli F, 许志琴等. 2006. 喀喇昆仑断裂的变形特征及构造演化. 中国地质, 33(2): 239-255.

李海兵, Valli F, 刘敦一等. 2007. 喀喇昆仑断裂的形成时代: 锆石SHRIMP U-Pb年龄的制约. 科学通报, 52(4): 438-447.

李海兵, Valli F, Arnaud N等. 2008. 喀喇昆仑断裂带走滑过程中伴随的快速隆升作用:热年代学证据. 岩石学报, 24(7): 1552-1556.

刘一峰, 兰海强. 2012. 曲线坐标系程函方程的求解方法研究．地球物理学报, 55(6): 2014-2026, doi: 10.6038/j.issn.0001-5733.2012.06.023.

王椿镛, 张先康, 丁志峰等. 1997. 大别造山带上部地壳结构的有限差分层析成像. 地球物理学报, 40(4): 495-502.

王椿镛, Mooney W D, 王溪莉等. 2002. 川滇地区地壳上地幔三维速度结构研究. 地震学报, 24(1): 1-16, doi: 10.3321/j.issn:0253-3782.2002.01.001.

王夫运, 段永红, 杨卓欣等. 2008. 川西盐源—马边地震带上地壳速度结构和活动断裂研究——高分辨率地震折射实验结果. 中国科学, 38(5): 611-621.

王世锋, 张伟林, 方小敏等. 2008. 藏西南札达盆地磁性地层学特征及其构造意义. 中国科学, 53(6): 676-683.

徐涛, 徐果明, 高尔根等. 2004. 三维复杂介质的块状建模和试射射线追踪. 地球物理学报, 47(6): 1118-1126, doi: 10.3321/j.issn:0001-5733.2004.06.027.

许志琴, 李海兵, 唐哲民等. 2011. 大型走滑断裂对青藏高原地体构架的改造. 岩石学报, 27(11): 3157-3170.

张先康, 王椿镛, 刘国栋等. 1996. 延庆—怀来地区地壳细结构——利用深地震反射剖面. 地球物理学报, 39(3): 356-364.

赵爱华, 丁志峰. 2005. 宽角反射地震波走时模拟的双重网格法. 地球物理学报, 48(5): 1141-1147.

赵烽帆, 马婷, 徐涛. 2014. 地震波初至走时的计算方法综述. 地球物理学进展, 29(3): 1102-1113, doi: 10.6038/pg20140313.

(本文编辑 何燕)

Tomographic imaging of velocity structure in upper crust based on correlated inversion ofVPandVS

WANG Xiao1,2, ZHOU Xiao-Peng1,2, ZHANG Xin-Yan1,2, BAI Zhi-Ming1*, TENG Ji-Wen1

1StateKeyLaboratoryofLithosphereEvolution,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2UniversityofChineseAcademyofSciences,Beijing100049,China

Presently seismic tomographic methods based on first arrivals of deep seismic sounding traveltime data have been widely used to image the upper crustal velocity structure. However, generally there are some difficulties in the inversion process using conventional methods: (1) How to construct reliable upper crust models using Pg and Sg data; (2) How to quantitatively estimate the reliability of reconstructed models; and (3) How to decrease the uncertainty and distortion validly.Based on the idea that theVPandVSanomalies from a same geological body are closely related each other, we develop a correlated inversion method:VPandVSinversion are carried out alternately, in which the correspondingVP/VSmodel and initialVS(orVP) model are updated using the obtainedVP(orVS) model after each inversion. The slowness perturbation resulted from the traveltime residuals is weighted and projected to the grid-points on the ray trace based on the updated correlation information, thus the velocity models are updated. Forward and inversion processes are respectively based on finite-difference calculations of travel times and the back-projection inversion method.Correlated inversion ofVPandVSmodels and sequential inversion technique are incorporated into the tomographic process to ensure the reliability of inversion results and avoid the distortion of the resultingVP/VSmodel. Comparison with the conventional inversion results using the code of Ammon and Vidale (1993), both of the derivedVPandVSmodels by our approach have been obviously improved, which are clearly closer to the original theoretical ones. While those constructed using the code of Ammon and Vidale (1993) have been seriously distorted. Moreover, the maximumVP/VSerror based on the conventional approach is up to 0.2, while the maximumVP/VSerror with our approach is only 0.02～0.03.Tests using two typical models show that the technology used in inversion tomography of the velocity structure within upper crust, which can effectively improve the reliability of the inversion results, can avoid the uncertainty and distortion resulted from the conventional singleVPandVSinversion process on a large scale. This method has been successfully applied to the western Tibetan plateau to reconstruct the velocity structure of upper crust beneath the Zhada-Quanshuigou deep seismic sounding profile. The derived models reveal the remained upper crustal velocity structure features resulted from plate collision in the western Tibetan plateau. Two similar obvious low-value anomalies ofVP,VSandVP/VSwith the maximum depth 8～10 km are present beneath the Zhada and Shiquanhe basin, which are probably the effect of the subduction of Indian plate beneath the Eurasian continent since Cenozoic or the closure of the Mesozoic Bangong-Nujiang oceanic basin. The disclosed Gar basin has a sedimentary thickness of 3～4 km, where the sedimentation and deep processes may be due to the gradual uplift of the Ayilariju Mountains and the strike-slip movement of the Karakoram fault.

Upper crust;VPandVSstructure; Correlated inversion; Finite-difference; Back-projection

10.6038/cjg20151011.

Wang X, Zhou X P, Zhang X Y, et al. Tomographic imaging of velocity structure in upper crust based on correlated inversion ofVPandVS.ChineseJ.Geophys. (in Chinese),58(10):3553-3570,doi:10.6038/cjg20151011.

中国地震局公益性行业科研专项(201408023)和国家自然科学基金(41374062,41174075，41274070，41474068)联合资助.

王晓，女，1989年生，博士研究生，主要从事地壳结构及地震资料偏移成像研究.E-mail: wangxiao@mail.iggcas.ac.cn

*通讯作者 白志明，男，1971年生，副研究员，主要研究方向为地壳结构成像及动力学研究.E-mail: bbzzmm@mail.iggcas.ac.cn

10.6038/cjg20151011

P315

2014-12-29，2015-09-15收修定稿

王晓， 周小鹏， 张新彦等. 2015. 上地壳纵横波速度结构相关反演成像方法.地球物理学报，58(10):3553-3570,