1. 引言
青藏高原位于印度板块和亚欧板块碰撞的交界区域，现今深部构造运动仍然十分活跃，是陆陆碰撞研究的天然实验场。青藏高原从南向北大致可分为五大地体：拉萨地体、羌塘地体、松潘–甘孜地体、昆仑–柴达木地体和祁连山地体，这些地体被若干缝合带划分开来 [1] [2] 。鉴于青藏高原地理位置的关键性、深部活动和构造的复杂性及其现今构造运动的活跃程度，长期以来被国内外学者广泛关注。迄今为止，已开展了大量有价值的地学研究工作并取得了一系列研究进展。在地震学领域主要采用了主动源探测(包括深地震测深和深地震反射剖面) [1] [3] 和天然地震成像方法(远震体波接收函数偏移成像方法、接收函数反演方法等) [4] [5] 。
莫霍面作为全球普遍存在的地球内部主要间断面，给出了地壳–地幔分界深度，可以作为分析研究块体深部接触变形关系的重要依据，并能用于推断板块运动活跃区域的深部动力学过程。对于青藏高原莫霍面深度的探测，国内外学者主要采用了上述主动源测深方法和接收函数偏移成像方法进行了大量研究工作。综合分析这些研究成果，已经可以给出青藏高原地区地壳厚度分布的整体轮廓。高锐等对这方面的研究工作进行了系统分析总结 [1] 。这些研究工作认为，青藏高原的莫霍面深度横向分布特征为：西部较深，东部较浅；以羌塘地体为中心，南部较深，北部较浅；最深和最浅的地壳厚度差可达40 km [1] 。但是，由于各种成像技术存在各自的优势和局限性，多种方法的研究结果并不能很好的达成一致，对于青藏高原内部几个主要地体莫霍面深度的认识存在几公里甚至十公里以上的差别。因此，通过结合多种地震学数据综合分析，联合解释，获取对青藏高原地壳厚度的全面、深入的认识仍然是十分必要的。
岩石圈–软流圈边界(LAB)不属于全球普遍存在的间断面，这一特征与莫霍面完全不同。地球科学的各个分支对这一边界的定义均有所不同，识别方法也就相应的存在差别。远震体波接收函数方法观测到的来自LAB的转换震相振幅强度远比来自莫霍面的转换震相弱的多，使得识别该震像和进一步确认其埋藏深度难度加大。目前，在地震学领域识别LAB主要依靠对地球内部速度间断面敏感的接收函数方法，单一的接收函数成像技术缺少对绝对速度的约束，成像结果受初始模型影响较大。因此，迄今为止对青藏高原地区的岩石圈厚度研究工作十分有限，已有的结果也比较粗略。
本文采用近期发展的基于贝叶斯理论的P波接收函数和S波接收函数联合成像方法 [6] 和INDEPTH IV部分台站的数据尝试性的做了青藏高原中部的联合反演成像研究。同时利用远震体波中的纵、横波速度信息，高频和低频资料相互约束，有利于获取对青藏地区地壳和岩石圈厚度的新认识。
2. 数据和方法
用INDEPTH IV在青藏高原中部布设的26个宽频带流动地震台站于2007年5月至2009年5月记录到的5.5级以上远震事件进行接收函数分析和后续的联合反演。这些台站自南向北依次分布在拉萨地体、羌塘地体、松潘–甘孜地体和柴达木地体四个青藏高原中部的主要构造单元。为分析在青藏中部各个块体之间的深部接触变形关系和地壳、岩石圈厚度穿越不同块体的横向变化特征，我们做了两条南–北向剖面和两条东–西向剖面对联合反演结果进行分析(图1)。
在对SAC格式事件波形数据进行重采样、滤波等预处理工作后，我们选择了高信噪比的远震记录进行下一步接收函数分析。其中，对于P波接收函数(PRF)选用了震中距30˚~90˚的全部信噪比良好的远震数据；对于S波接收函数(SRF)，为获取岩石圈–软流圈边界(LAB)转换震相信息，选择震中距55˚~80˚的远震事件 [7] 。P波接收函数的估计采用时间域迭代反褶积算法 [8] ；S波接收函数估计采用接收函数最大或然性反褶积方法 [9] 。
基于贝叶斯理论的P波接收函数和S波接收函数联合反演方法是一种比较新的岩石圈地幔地震学成像方法。贝叶斯反演理论建立在概率论的基础上，该联合反演方法的本质是用P波接收函数的高频信息和S波接收函数的低频信息相互约束，从而获取更为可靠的岩石圈地幔(约300 km深度范围)的速度结构模型 [6] 。相比单一的P波接收函数反演增加了约束信息，同时利用S波接收函数对于LAB的识别能力获得对岩石圈底部边界的深度估计。关于该联合反演方法的详细解析读者可参阅文献 [6] 。
Figure 1. Topographic map with station locations (QD: Qaidam Terrane, QT: Qiangtang Terrane, SG: Songpan-Ganzi Terrane, LT: Lhasa Terrane)
图1. 地形和台站分布(QD：柴达木地体；QT：羌塘地体；SG：松潘–甘孜地体；LT：拉萨地体)
3. 结果
图2给出了本文采用的贝叶斯联合反演的两个实例。其中台站D01 (图2(a))位于羌塘地块和松潘–甘孜地块的交界部位，台站D24 (图2(b))位于羌塘地块内部。左图为初始P波速度模型和联合反演得到的P波速度模型，中图为初始S波速度模型和结果S波速度模型，右图为时间域的接收函数观测波形和合成理论波形。从反演速度模型中可以观察到，两个台站下方的中上地壳均存在约10 km厚的壳内低速层，这一特征可与5s左右P波接收函数中的负震相对应；两个台站的地壳厚度均为70 km。相比台站D01 (LAB深约100 km)，位于羌塘内部的台站D24下方岩石圈(约为150 km)更厚，这与S波接收函数波形相对应。需要说明的是，整个联合反演过程在频率域中进行，既提高了计算效率也节省了存储空间。我们为给出直观的分析结果，通过傅里叶逆变换将频率域波形转换到时间域由图2中的右图给出。图2表明，该联合反演方法能够在青藏高原这一深部结构复杂的区域得到比较可靠的速度结构模型，从而估计出地壳和岩石圈厚度。
为了分析位于青藏中部的研究区内各主要构造块体之间的深部结构关系，我们划分了四条剖面包括两条南北向剖面S1-N1，S2-N2和两条东西向剖面W1-E1，W2-E2 (剖面位置见图1)。其中，S1-N1剖面由6个台站组成，自南向北依次穿过拉萨地体、松潘–甘孜地体和柴达木地体(图3)。沿剖面走向在P波接收函数和S波接收函数剖面图(图3(a))中，均可追踪到清晰的来自莫霍面的转换震相Pms和Smp。该剖
面下方的地壳厚度自南向北逐渐变浅，这一分布特征还可以从速度结构剖面图(图3(b)，图3(d))中清楚的观察到。剖面S1-N1南部的拉萨地体(LS)下方莫霍面厚度超过70 km，前人研究结果认为拉萨地体下方平均莫霍面深度约为70 km [10] 。在松潘甘孜地体(SG)下方，莫霍面深度有明显的起伏变化(60~70 km)。由于松潘–甘孜地体的地理位置东西向跨度较大，高锐等分析认为，该地体下方的地壳厚度有自西向东逐渐变浅的趋势。本文的S1-N1剖面主要穿越了松潘–甘孜地体的西缘。剖面S1-N1最北端的两个台站逐渐进入柴达木地体(QD)南缘，从图3中可以看到明显的近地表沉积层，地壳厚度约为60 km。在图3(a)的S波接收函数剖面中，无法清晰追踪到来自LAB的转换震相SLp，因此未在图中标出。我们的联合反演方法同时拟合了时间窗内的转换震相及其多次波，因此对于整个反演深度范围内的结构信息均有约束。根据联合反演的结果，在速度结构剖面图中，标出了LAB的估计深度。沿着剖面S1-N1自南向北岩石圈厚度呈现浅–深–浅的横向变化特征，其中岩石圈最厚达约170 km，最浅处仅为100 km。此外，沿剖面的波速比横向变化(图3(c))与P波速度变化特征比较吻合。
Figure 3. Receiver functions and velocity structure distribution along profile S1-N1
图3. 沿剖面S1-N1的接收函数和速度结构横向变化
南北向另一条剖面S2-N2由12个台站组成，自南向北依次穿过羌塘地体、松潘–甘孜地体和柴达木地体(图4)。该剖面近地表构造复杂，相应的P波接收函数壳内波形比较复杂，Pms转换震相不易简单的追踪到。相反地，S波接收函数波形中Smp转换震相比较清晰(图4(a))。我们认为，这可能是由于该剖面下方的莫霍面不是简单的一阶速度间断面而是渐变的、梯度带类型的间断面，导致频率较高的P波接收函数无法有效识别，而相对低频的S波接收函数更适合识别这类间断面。从P波速度结构剖面(图4(b))可以观察到，壳内低速体的分布自南向北影响的深度范围逐渐变浅。羌塘地体和松潘–甘孜地体下方的壳内低速体几乎涉及整个地壳深度；剖面S2-N2最北端的柴达木地体下方的壳内低速体只存在于中上地壳。
该剖面下方地壳厚度起伏变化剧烈，位于松潘–甘孜地体的台站D01下方莫霍面达到84 km，位于柴达木地体的剖面最北端台站A01下方地壳厚度仅为56 km。剖面S2-N2下方的LAB深度的横向起伏变化更为剧烈，变化范围为100~160 km。分析认为，该剖面下方的岩石圈地幔结构横向变化如此剧烈可能是因为剖面恰好位于羌塘地体和松潘–甘孜地体的交界部位，这一区域地表构造复杂，深部运动可能仍然很活跃。
东西向剖面W1-E1自西向东依次穿过羌塘地体和松潘–甘孜地体，共由8个台站组成(图5)，位于羌塘地体的2个台站莫霍面较为一致，约为64 km；位于松潘–甘孜地体的6个台站地壳厚度横向变化明显，主要体现为西部较深，东部较浅，变化范围在64~84 km。台站D19位于两个块体的交界部位，壳内低速结构较为特殊。该剖面下方岩石圈厚度的横向分布特征与莫霍面分布特征相对应，在羌塘地体下方变化平缓，在松潘–甘孜地体下方起伏变化剧烈。剖面W1-E1下方的LAB深度变化范围是110~160 km。
Figure 4. Receiver functions and velocity structure distribution along profile S2-N2
图4. 沿剖面S2-N2的接收函数和速度结构横向变化
Figure 5. Receiver functions and velocity structure distribution along profile W1-E1
图5. 沿剖面W1-E1的接收函数和速度结构横向变化
上述结论通过对比分析剖面W1-E1的接收函数和速度结构结果得出。
东西向另一条剖面W2-E2位于研究区最南端，共7个台站(图6)。该剖面主要位于拉萨地体，最东端的2个台站位于羌塘地体。观察接收函数波形(图6(a))可以发现，在拉萨地体下方沿着该剖面高程起伏较大，但莫霍面深度横向变化平缓，变化范围为72~74 km。前人采用接收函数方法和深地震测深方法得出了相似的结果 [1] [10] ，他们认为拉萨地体下方的地壳厚度为70~80 km。剖面向东进入羌塘地体后莫霍
Figure 6. Receiver functions and velocity structure distribution along profile W2-E2
图6. 沿剖面W2-E2的接收函数和速度结构横向变化
面变深为76 km。在岩石圈地幔深度，LAB的横向变化特征与莫霍面类似，即拉萨地体台站下方岩石圈底边界平缓，在羌塘下方LAB深度变化大。
4. 结论和讨论
本文利用近期发展的基于贝叶斯理论的P波接收函数和S波接收函数联合反演方法和INDEPTH IV布设在青藏高原中部的部分台站记录到的远震事件，尝试性的研究了青藏高原中部岩石圈地幔深度的速度结构，得到了研究区地壳和岩石圈厚度在各块体下方的横向变化特征。鉴于青藏高原深部结构研究的重要性研究价值，众多国内外学者已在青藏地区做了大量观测和研究工作。对于莫霍面和LAB深度的地震学研究成果主要采用了基于主动源探测的深地震测深和深地震反射剖面方法，以及基于宽频带流动地震台阵观测的远震体波接收函数H-K扫描叠加方法和接收函数偏移成像方法。本文采用的联合反演方法充分利用了远震事件中体波所携带的转换震相和多次波信息，将高频和低频信息相结合，得到了青藏高原中部下方岩石圈地幔深度的绝对速度结构模型；在此基础上给出了研究区莫霍面和LAB深度变化分布特征。结果表明，青藏高原中部四个主要构造块体拉萨地体、羌塘地体、松潘–甘孜地体和柴达木地体在岩石圈地幔深度具有不同的深部结构特征；各块体下方的地壳和岩石圈厚度差别明显；块体交界部位的莫霍面和LAB深度变化尤为剧烈。其中，拉萨地体下方的地壳厚度约为70~76 km，岩石圈厚度约为110~160 km；羌塘地体下方64~76 km，岩石圈厚度在110~170 km变化；位于松潘–甘孜地体的台站D01下方莫霍面达到84 km，位于松潘–甘孜地块北端与柴达木地体交界部位的台站D10莫霍面深度为64 km。松潘甘孜地块的岩石圈厚度约为100~160 km。研究区中位于柴达木地体的台站较少，这些台站的莫霍面深度平均值为56 km，LAB深度平均值为120 km。
需要说明的是，本文采用的联合反演方法虽然综合利用了高低频信息和多种震相信息相互约束，但仍然缺少足够的绝对速度信息约束。加入面波频散数据的约束可能有利于获得更为稳定可靠的反演结果。
致谢
作者感谢IRIS DMC提供远震波形数据。
基金项目
国家自然科学基金(41574045，41774055)和地震动力学国家重点实验室自主基金(led2013a06)共同资助。