川滇地区地处中国西南，地质背景复杂。西北部是构造运动活跃的青藏高原，东部是稳定的扬子地块。作为两者的连接处，川滇地区断层密布，地质构造复杂，地震频发^[1]。川滇地区壳幔结构，是研究青藏高原隆升和扩张机制、地震活动性等相关问题的重要依据。

川滇地区的壳幔结构研究早在20世纪80年代就已展开，早期研究多采用人工震源测深的方法^[2-3]。随着技术的发展和数据的增加，各种被动源方法，比如接收函数^[4-6]、面波频散^[7-8]和联合反演^[9-12]等手段都在川滇地区得到应用。

虚拟地震测深方法（Virtual Deep Seismic Sounding，VDSS）是近些年新发展的S波数据测深技术。VDSS利用天然远震的全反射SsPmp波来研究结构，具有信号能量大、信噪比高、无需叠加、不容易受到近地表沉积和非均匀性影响的优势。Tseng等^[13]利用波形拟合测量到时差，给出了青藏高原一条约500 km长的测线下方反射点的地壳厚度估计，与接收函数偏移结果有很好的一致性；Yu等^[14]利用偏振分析的方法分离震相，提取了入射Ss波，扩展了VDSS方法的适用性；Yu等^[15]提出了根据直达Ss到时异常校正到时差的方法，降低了地壳厚度的测量误差；Kang等^[16]利用线性拟合，同时给出地壳P波速度和厚度的估计；Parker等^[17]提出利用包络线测量到时差的方法；Liu等^[18]提出了同时测量相位和到时的方法，丰富了到时差的测量方法。VDSS方法也在鄂尔多斯^[19-21]、川滇^[22-23]、美国西部^[15]、美国东南部^[17]等地区的地壳结构探测得到应用。

本文利用VDSS方法研究川滇地区的地壳厚度分布。首先测量了研究区域内台站上的SsPmp震相与Ss震相的到时差，并根据直达Ss的到时异常对到时差进行校正，最后偏移到深度域，得到了川滇地区地壳厚度分布，为川滇地区的壳幔结构研究提供了新的约束。

研究数据来自中国地震局在中国西南布设的116个宽频带台站^[24]（图1）。这些台站较好地覆盖了本次研究区域（98 °～108 °E，20 °～34 °N）。VDSS方法在地震筛选时通常有以下3个要求。

图  1  研究区域内台站分布图

1）震中距在30 °~55 °之间。震中距小于30 °，可能会记录到Ss三叉震相，增加复杂性；而震中距大于55 °，SsPmp不再是全反射波，振幅下降。

2）地震震级5级以上，深震更佳，这两点可以使得震源时间函数相对简单。

3）具有清晰的Ss和SsPmp震相，方便提取和分析。

根据以上要求，本文选取一个最佳地震进行研究。本文选取的地震，发震时刻为2020年07月06日22时54分48.8秒，震源位置为（5.69 °S，110.55 °E），深度535 km，体波震级6.7（参数来自GCMT catalog^[25-26]）。共有100个台站记录到该事件，记录的典型波形如图2b中垂向分量V和径向分量R所示。

图  2  KMI台站的波形图

VDSS方法是一种使用远震全反射SsPmp震相研究地壳结构的有效方法。SsPmp震相的射线路径如图3所示：远震Ss波到达地表D点后，转换为P波向下传播，在莫霍面C点全反射，在地表B点被仪器记录（图3）。

图  3  SsPmp射线路径示意图

对于单层均匀水平模型，有

式中：T为B点SsPmp波与Ss波到时差；H为地壳厚度；$ v $为地壳P波速度；$ p $为射线参数^[13]。VDSS的数据处理流程，一般是先进行震相分离，提取入射Ss波，再测量到时差，最后把到时差偏移到深度。

震相分离是将垂直方向记录V（t）和径向记录R（t），按照Ss波和SsPmp波偏振的方向，分解成伪S分量Si（t）和伪P分量Pi（t）的过程。Yu et al^[14]提出，根据偏振分析进行分解的方法：将V（t）和R（t）合成为矢量记录，遍历Ss波和SsPmp波的潜在偏振方向$ {\theta }_{\mathrm{S}} $和$ {\theta }_{\mathrm{P}} $，在每一种偏振方向组合下计算每个矢量记录点到2个偏振方向的距离${d}_{\mathrm{P},\mathrm{S}}\left(t\right)=V\left(t\right) $$ \mathrm{sin}{\theta }_{\mathrm{P},\mathrm{S}}-R\left(t\right)\mathrm{cos}{\theta }_{\mathrm{P},\mathrm{S}}$（下角标P、S表示伪P分量和伪S分量的潜在偏振方向），计算使得每个点到更近的直线的距离平方和最小的偏振方向$[{\theta }_{\mathrm{P}},{\theta }_{\mathrm{S}}]=\underset{{\theta }_{\mathrm{P}},{\theta }_{\mathrm{S}}}{\mathrm{m}\mathrm{i}\mathrm{n}}{\displaystyle\sum }_{t}{\mathrm{m}\mathrm{i}\mathrm{n}}^{2} $$ ({{d}}_{\mathrm{P}}\left(t\right),{{d}}_{\mathrm{S}}\left(t\right))$。一个偏振分析的例子如图2a所示，分解结果如图2b中的Si和Pi所示。

因为SsPmp震相是全反射震相，相对于Ss产生相移，所以两个震相的到时差并不能直接通过简单的峰值时刻相减或互相关等方法得到。常用的测量到时差的方法有两种，一种是波形拟合方法^{[14-16, 19]}，另外一种是同时测量相移和到时的方法^{[18, 20-21]}。

波形拟合方法，通常假定地壳为单层介质，通过给定地壳和上地幔的波速和密度，变化地壳厚度计算一系列介质响应，卷积上之前分解得到的入射Ss波，得到模拟记录；寻找与观测波形残差最小的模拟记录，将对应的理论到时差作为测量结果。

同时测量相移和到时的方法，是对分解后的伪S波分量进行相移和到时差遍历，然后和伪P波分量进行互相关，取互相关系数最大时的相移和到时差为测量结果。

拟合波形的方法，因给定了速度模型，到时测量可能出现系统误差；而同时测量相移和到时的方法中，相位测量的随机误差较大。本文结合两者，给定地壳P波速度6.0～6.2 km/s、上地幔P波速度7.9～8.1 km/s、波速比固定为1.732。对于每个台站的记录，先根据射线参数计算出理论相移范围，然后在此相移范围内同时测量到时和相移。

地壳和地幔中的异常体，会导致真实模型与平层假设不一致，直接使用公式进行到时差深度转换会存在误差，通常根据直达Ss波异常对公式进行修正^[15,21,23]。

式中：t为到时；第一个下角标S，P（后文中）分别表示Ss和SsPmp震相，第二个下角标D，B表示SsPmp震相的地表转换点和地表接收点（图3），第三个下角标O，T，A分别表示观测到时、理论到时和到时异常。

其中：$ p $表示理论射线参数，给定震源，台站位置和一维参考模型后，可以通过TauP^[27-28]方法计算得到；$ l $表示两个下标对应点间的距离。

式中：$ v $为地壳平均P波速度。令

联立上式，有

式中：$ T $为接收点B的SsPmp与Ss的到时差，与式（1）中的$ T $一致；dt为B点和D点的Ss到时异常差，也就是待求的修正量。

在数据处理中，需要先根据每个台站的直达Ss波到时异常，通过线性插值给出研究区域的Ss到时异常分布图。对于每个台站B的记录，先测量到时差$ T $，然后需要找到对应的地表反射点D的Ss波到时异常值$ {T}_{\mathrm{S}DA} $，与台站自身的到时异常值${T}_{{\rm S}BA}$做差，得到修正量$ dt $，两者相加得到修正后的到时差$ {T}' $，再进行深度偏移。因为地表反射点的位置与$ {T}' $耦合，所以通常采用迭代的思路，利用$ {T}'=T $作为初值，计算地表虚拟震源位置、莫霍面反射点位置和对应深度，然后根据目前的虚拟震源位置计算修正量dt，再令$ {T}'=T+dt $，重新计算虚拟震源位置和对应深度，循环计算直至收敛。

选取地震后，对每个台站的波形数据去均值、去线性趋势和去尖灭，做2～50 s带通滤波。根据地震和台站的位置，计算反方位角、理论到时和射线参数，并截取理论到时前20 s至后30 s的数据，将北、东分量旋转到径向、切向分量。对径向、垂直分量进行偏振分析，分解得到伪S、伪P分量。对不同台站的伪S分量进行两两互相关，计算到时差和互相关系数，选取互相关系数大于0.5、0.7和0.8的结果分别进行多道互相关，得到每个台站的到时异常。除去没有记录的台站和3个结果中最大值与最小值差异大于1 s的3个台站，以剩下97个互相关系数大于0.8台站的结果进行线性插值，得到研究区域内的直达Ss波到时异常分布图（图4）。

图  4  Ss波到时异常图

将全部台站的直达Ss波按照到时异常校正后，叠加得到平均入射Ss波M（t）（图2b）。根据射线参数计算理论相移范围，利用M（t）在每个台站P分量上同时测量相移和到时差（图5）。同时，利用包络线的方法测量了到时差^[17]，当两个到时测量结果差距大于0.5 s时（图1中灰色三角形）就舍去此数据，以控制数据质量。假设$ v $=6.1 km/s，利用台站的到时差和研究区域的直达波到时异常图，迭代求解地壳厚度（图6）。

图  5  台站VDSS到时差分布图

图  6  地壳厚度分布图

全部台站中，直达Ss波到时异常最小值为−3.94 s，最大值为3.52 s，异常的分布表现出规律性，与大的地质块体分区较为一致（图4）。东北部四川盆地表现为一致的负异常，即实际到时早于理论到时。这与面波成像^[8-9]显示四川盆地深部具有高速异常特征的结果相符。南部云南地区表现为一致正异常，实际到时晚于理论到时，这也与面波成像^[8-9]所显示的云南地区深部具有低速异常特征的结果相符，可能与缅甸俯冲带弧后拉张有关。在峨眉山大火成岩省内环区，到时异常显著小于南边相邻的云南地区，可能与相对应的高速异常有关。西北部青藏高原区域到时异常横向波动较大，规律性不显著，可能反映了其深部结构的复杂性。

全部台站中，SsPmp与Ss到时差最小值为3.19 s，最大值为15.07 s（图5）。到时差的巨大差异，主要反映的是研究区域巨大的岩石圈厚度差异和复杂的结构。到时差最小的区域出现在南边的云南地区和四川盆地，这两个区域内部不同台站的到时差差异较小，一致性较好；到时差最大的区域，出现在云南地区与青藏高原的交界处（26°N，103°E附近），地壳开始变厚的区域。而青藏高原内部，到时差小于交界处，横向波动较大，规律性不显著，说明结构比较复杂。

进行直达Ss到时异常校正后，四川盆地地壳厚度为32～42 km，云南南部的地壳厚度为32～37 km，云南地区与青藏高原交界处往西北延伸，地壳厚度从40 km递增到60 km。校正后地壳厚度与前人结果^{[8-9, 11]}更为一致。在云南地区与青藏高原的交界处（26°N，103°E）附近，若直接利用式（1）进行到时差深度转换，交界处地壳厚度在60～68 km之间；而进行直达Ss到时异常校正之后，地壳厚度降低到49～60 km，相同台站对应的莫霍面反射点的深度差异可达13 km，这说明了直达Ss到时异常校正的必要性和有效性。同时，校正后相邻反射点的深度值更为接近，也说明了校正的有效性。

校正后，本文研究结果（图6b）与前人结果^[4-6]在部分点位上也存在差异。差异较大的点主要集中在两个区域。一是部分研究区域边缘的反射点，台站稀疏是直达波到时异常与相邻反射点差异较大的原因，校正量误差较大导致深度误差较大；另外一个区域是地壳结构剧烈变化的地区，比如龙门山断层，其他断层交汇的区域附近。本文计算直达Ss波到时异常分布图时，是利用每个台站的异常值，采用的线性插值的方式得到的。当地壳结构剧烈变化时，简单的线性插值不能反映真实的到时异常。同时VDSS方法反映的是地表转换点到台站之间（水平距离大约120 km）的一个平均结果，与反射点处确切的地壳厚度存在差异。亢豆等^[22]的研究发现断层附近VDSS方法得到的地壳厚度偏小，与本文结果类似。另外，本文的结果也显示了龙门山断裂带两侧从东南往西北地壳厚度突然增加的现象。

利用VDSS方法研究川滇地区地壳厚度分布显示，在地质结构复杂的地区，直达波Ss到时异常校正是必要的。在地壳厚度变化剧烈的地区，VDSS方法的结果存在较大的误差。校正后四川盆地地壳厚度为32～42 km，云南地区南部地壳厚度为32～37 km，云南地区与青藏高原交界处，地壳厚度从40 km递增到60 km，地壳厚度结果与前人一致。本文用S波列相关震相确定地壳厚度，为川滇地区的壳幔结构研究提供了新的约束。

致谢　研究工作得到北京大学高性能计算校级公共平台支持，中国地震局地震预测研究所和地球物理研究所为本研究提供了地震波形数据，在此一并感谢。