青藏高原中部INDEPTH-III剖面低速层研究

牛 潇 ， 郑洪伟 ， 贺日政 *， 李娱兰 ， 严江勇 ， 李宗旭

1)中国地质科学院地球深部探测中心， 北京 100037； 2)中国地质科学院地质研究所， 北京 100037；3)中国科学技术大学地球和空间科学学院， 安徽合肥 230026

由多个地体在不同时期拼合而成的青藏高原(常承法和郑锡澜， 1973)， 在新生代以来经历了印度大陆北向碰撞与俯冲(Yin and Harrison， 2000)， 形成了地球表面年代最新、规模最大、海拔最高的构造汇聚区并且碰撞过程持续至今， 是研究陆陆碰撞的天然实验室(高锐， 1997； Yin and Harrison， 2000；Haines et al.， 2003)。始于20世纪70年代多国科学家联合在高原内开展的深部探测研究(Gao et al.，2005)， 为探索青藏高原的隆升机制提供了初步框架认识， 如地壳挤压叠覆加厚(Dewey and Burke，1973； England and Houseman， 1986)、 地 壳俯 冲(Roger et al.， 2000； Tapponnier et al.， 2001)、底侵加厚(Molnar and Tapponnier， 1977)等模型。随着研究的深入， 青藏高原内部仍有很多深部结构、构造和性质的问题有待研究(高锐， 1997)， 如青藏高原低速层分布及其构造成因。本次研究选择探测研究程度较高的INDEPTH-III剖面所经过的班公湖—怒江缝合带中段(Gao et al.， 2005； 刘国成等， 2014)， 利用接收函数反演技术研究该区域下的低速层分布特征并探讨可能的构造成因。

位于藏北高原中部的班公湖—怒江缝合带(Yin and Harrison， 2000)是青藏高原重要的缝合带， 以此为界的拉萨地体和羌塘地体的壳幔属性具有明显的南北差异， 如羌塘下的壳内低速层分布特征(李永华等， 2006a； 严江勇等， 2019)， 青藏高原普遍存在低速层(李永华等， 2006b； 贺日政等， 2007； 滕吉文等， 2012)， 或低速层并不普遍存在(Owens and Zandt， 1997； 吴庆举和曾融生， 1998)。因此对班公湖—怒江缝合带两侧低速层分布特征的研究尤为重要。

利用接收函数一维波形反演方法(刘启元等，1996)获取跨越班公湖—怒江缝合带的INDEPTH-III剖面速度结构特征， 再收集最近20多年利用 INDEPTH-III项目数据发表的研究成果(Zhao et al.， 2001； Michie et al.， 2004； Ross et al.，2004； Tian et al.， 2005； He et al.， 2010)以及其它相关探测成果(邹长桥等， 2012； Gao et al.， 2013； 刘国成等， 2014； 吴蔚等， 2017； 严江勇等， 2019)基础上，综合分析班公湖—怒江缝合带两侧的壳内低速层分布特征，不仅对理解青藏高原内部物质变化、板块俯冲、地幔对流等深部物理化学过程具有非常重要的意义(滕吉文等， 2012)， 而且从侧面认识青藏高原的形成和演化机制。因此对藏北高原中部青藏高原不同地块内地壳结构以及低速层的研究尤为重要。

1 区域构造背景

青藏高原被普遍认为是由多个地体在不同地质时代在经历了多期次碰撞拼合而成的一个复杂的地质构造单元(图 1； Chang et al.， 1986)。在其内部的相邻地质构造单元间的标志为一系列近东西走向的缝合带， 自南向北依次为， 雅鲁藏布江缝合带、班公湖—怒江缝合带、龙木错—双湖缝合带、金沙江缝合带、阿尼玛卿—昆仑山缝合带， 祁连缝合带(图1)。位于青藏高原中部的班公湖—怒江缝合带， 作为中生代关闭消亡的中特提斯大洋的产物， 其西起班公湖， 向东经改则、东巧、丁青、嘉玉桥至八宿县的上林卡， 经左贡扎玉、梅里雪山西坡与昌宁、孟连带相通(潘桂棠等， 2004)， 东西绵延数千公里(Yin and Harrison， 2000)。

作为拉萨地体和羌塘地体间的重要构造边界，班公湖—怒江缝合带与其南北两侧， 不论是构造还是地表地质甚至是岩浆活动事件， 存在着显著差异(尹安， 2001； 潘桂棠等， 2004)， 如班公湖—怒江缝合带以北的新生代火山岩浆活动主要与亚洲地幔大面积拆离(罗照华等， 2006)或部分断离作用有关(He et al.， 2014)， 而其南部的冈底斯带大多与碰撞机制有关(罗照华等， 2006)。班公湖—怒江缝合带还是一个重要的多金属成矿带(唐菊兴， 2019)。因此， 探测班公怒江缝合带及其两侧的地壳上地幔深部结构特征并开展区域对比研究， 显得尤为重要。利用INDEPTH-III宽频带地震观测数据(Tilmann et al.，2003)开展接收函数分析， 研究班公湖—怒江缝合带两侧的地壳结构差异， 非常之必要。

2 接收函数方法及数据

接收函数是从三分量地震台站记录中分离出垂直分量与水平分量做反卷积后得到的地震时间序列，具有很好的台站下方结构分辩率， 避免了震源和传播路径因素的影响， 是研究地球深部重要的方法之一(Ammon et al.， 1990； Langston， 1979)。由于地下结构会对P-to-S转换波(简写为Ps)产生扰动， Ps传播过程中记录了地下介质结构信息， 不仅用于反演壳幔分层速度(Langston， 1979)， 地壳厚度和泊松比估算(Zhu and Kanamori， 2000)， 还用于壳幔结构成像(Tian et al.， 2005)。为了提取接收函数， 目前广泛采用方法主要是时间域迭代反褶积(Owens and Zandt，1997)和频率域反褶积(Langston， 1979)。在此基础上，利用不同的方法获取台站(或台阵)下方结构特征。

图1 研究区内的区域构造特征及本文所用的台站分布和地震事件分布图(缝合带行迹来自尹安, 2001)Fig. 1 The regional structural characteristics of the study area, the station distribution and seismic event distribution(suture belt after YIN, 2001)

接收函数反演分为线性反演和非线性反演。线性反演(Ammon et al.， 1990； Owens and Zandt， 1997)，具有运算速度快优点， 但需要先验结构模型， 而且要求台站下方介质结构特征光滑连续且不能存在断裂。显然， 该方法不能满足“地下介质属性及其结构特征具有明显的非连续性”。为克服线性反演方法的局限性， 一些非线性反演方法， 如遗传因子(Shibutani et al.， 1996)， 接收函数复谱比(刘启元等，1996)和模拟退火(高星等， 2005)等方法运用而生。本文采用了刘启元等(1996)提出的接收函数复谱比非线性反演方法。该方法不依赖初始模型而且还对初始模型没有过分严格要求。

接收函数复谱比反演按照 Tarantloa的构造公式所示的目标函数S(m)：

其中g(m)和d分别为反演中的理论和观测的数据矢量，g(m)和d的分量相应于不同频率的接收函数复谱比，m为模型参数矢量， 角标p表示其先验信息，CD和CM分别为数据和模型的协方差矩阵， *表示复共轭， 上脚标T表示转置。

之后对目标函数用共扼梯度法进行优化， 主要计算目标函数的梯度Nγ即：

g∂和mi分别是预测数据矢量和模型矢量的分量， I为单位矩阵， R表示模型空间中解的分辨率， 角标N表示第N次迭代的解。最后经过多次迭代， 选取最佳拟合波形， 得到反演结果。接收函数复谱比正演计算采用了反射率法。对于确定的慢度， 它仅依赖于地表附近介质绝对参数的量， 而数值检验结果表明接收函数复谱比反演的结果与初始模型的选择无关(刘启元等， 1996)。

本文所用的数据来源于INDEPTH-III计划布设的宽频带地震仪所记录的远震事件， 地震记录的采样率为50 Hz， 台间距大多在10 km左右。地震事件参数来自美国地质调查局(http：//www.usgs.gov)， 地震事件震级Mb≥5.5， 震中距在 30°～90°之间。

本次研究采用的数据处理流程如下： 首先对宽频带地震仪记录的三分量地震事件波形数据进行处理， 即(1)对原始数据以 20 Hz的采样率重采样， 截取P波初动前10 s和之后100 s的波形； (2)去均值、线性趋势及仪器响应等； (3)利用后方位角对原始的Z-N-E三分量地震波形数据旋转到 Z-R-T(垂向， 径向和切向)坐标系， 做带通滤波(带宽 0.05～2 Hz)；(4)应用时间域迭代反褶积对垂向和径向分量进行反褶积， 高斯系数选取 2.5； (5)利用 CrazySeismic 软件(Yu et al.， 2017)挑选高信噪比(不小于 2.0)的接收函数。而对于较低信噪比的接收函数， 采用基于S变换的相位滤波技术(曲中党等， 2015)获取较李永华等(2006a)更多台站下方的接收函数。(6)利用接收函数复谱比非线性反演(刘启元等， 1996)获取了台站下方的一维S波速度结构特征。(7)构建适合剖面所在区域的标准速度模型： 由于 IASPEI91标准模型的Moho深度较青藏高原实际Moho深度浅得多，依据前人 H-k 叠加结果(李永华等， 2006b； 刘国成等，2014)修改了IASPEI91标准速度模型中的Moho深度。(8)利用地震波层析成像思想， 基于改进的IASPEI速度模型， 计算改进后的模型与反演结果之间的扰动， 获取了剖面下的速度异常分布特征，突出壳内低速及莫霍面变化特征。

3 反演结果与分析

依据前述规则， 获得了INDEPTH-III剖面中信噪比较高的23个台站下方的接收函数， 其中ST02、ST03等18个台站(图1)下方的接收函数因信噪比很低而舍去。图2展示了高信噪比的23台站下方的接收函数的倾斜叠加(slant stack)结果， 清楚地展示 P波初至和P-s转换震相。之后利用接收函数复谱比非线性反演方法(刘启元等， 1996)获得了该23个台站反演结果(图 3)， 显示所有垂向分量的反演结果与实际波形相关系数都为 0.96， 而径向分量的相关系数大多都在0.95以上。这主要是由于加入了相位滤波技术(曲中党等， 2015)提高了拟合系数， 使得反演结果更加可靠。对每个台站下的速度结构进行中值滤波， 有效地抑制反演结果中的小尺度异常并保留了主要界面特征(图4)。

对于低速层的判定， 本文采取如下规则：(1)低速层的厚度不小于 5 km， (2)要比改进的IASPEI91速度模型的速度要低， 即同一深度的速度值低于改进后的IASPEI91速度模型。

本项研究结果与李永华等(2006a)的结果相比，二者的共同点是在拉萨块体北部和羌塘地体中下地壳有低速层分布； 二者的差异是， 本文获得的中上地壳内低速层分布特征比较精细， 部分低速区呈现速度反转即低速区比上下区域速度值都低， 更能体现低速特征(图 5)。这是由于使用了不同的反演算法， 李永华等(2006a)使用的是遗传算法的反演方法。尽管不论线性还是非线性一维波形反演都具有较大的不确定性(Ammon et al.， 1990)， 但相比之下， 本文所用的复谱比非线性反演方法(刘启元等，1996)不依赖于通过先验信息构建的初始速度模型(李永华等， 2006a)， 而取决于所挑选的接收函数的信噪比。

在不同台站下方的低速层分布特征如下： (1)在ST01、ST08、ST23、ST28、ST29、ST34、ST35、ST37台站下方低速层厚度较大； (2)在ST05、ST26、ST36台站下方的低速层较薄； (3)ST37、ST38台站下方的低速层与吴蔚等(2017)展示的临近台站(NQT10、NQT12)速度结构相似； (4)整体上拉萨地体比羌塘地体低速层埋深较浅。

跨班公怒江缝合带及两侧的INDEPTH-III深反射探测结果(Ross et al.， 2004)显示班公怒江缝合带顶部约有小于 2 km 厚度的沉积盖层， 以及上地壳(约 25 km以浅)表现为无反射， 本文 S波速度结构也在 32°N附近 ST15、ST16、ST19的台站下方近地表出现低速区。羌塘盆地内的S波速度结构(吴蔚等， 2017； 严江勇等， 2019)显示沿 88.5°E 测线南羌塘低速层埋深约 30 km， 北羌塘盆地内广泛存在的低速层明显变浅， 这与本文得到的S波速度结构(图4)相似。李永华等(2006a)对INDEPTH-III台站开展接收函数分析， 研究显示拉萨地体北部和羌塘地体下方中地壳有低速层存在， 也与本文结果有较好的一致性(如图4)。不同方法及不同数据集获得的低速层分布特征较为一致， 表明本文利用一维速度反演获得结果是可靠的。

图2 研究剖面的地形、台站位置名称及台站下方接收函数倾斜的叠加Fig. 2 Surface elevation, station position name and slant-stack of receiver function along the profile

利用三次样条插值方法对反演获得台站下方一维速度结构插值得到了测线下方二维速度结构(图 5b)， 显示在拉萨地体和羌塘地体普遍存在不连续的低速体。为了更好的突出低速层分布特征， 利用地震波层析成像方法思想(Zhao et al.， 1992， 1994；Zhao， 2009)， 即与参考模型进行对比并求取速度异常扰动。由于IASPEI91速度模型是全球均匀模型，不适合巨厚地壳的青藏高原研究(He et al.， 2010)。为更好突出显示青藏高原内部低速层分布特征， 我们对该全球标准速度模型做了改进， 即将利用前人通过H-k方法得到的地壳厚度(李永华等， 2006a； 刘国成等， 2014)， 换算为 Moho 深度， 修改 IASPEI91参考模型。再利用计算反演速度(图5b)与改进后参考模型之间的速度异常扰动的公式(4)， 获取剖面下的S波速度异常扰动， 突出显示低速层分布特征。

式中LVZ表示扰动大小， 用百分数表示，VS表示某一深度处的剪切波速度，Vref表示该深度处的经修改的参考模型剪切波速度。通过三次样条插值方法获得了剖面下方的速度异常结构， 突出显示了S波速度异常体分布特征(图5c)。

图3 23个台站接收函数波形拟合结果Fig. 3 The comparison result between original waveform and inverted waveform by nonlinear inversion for 23 stations

图4 沿INDEPTH-III剖面下方地壳S波速度结构Fig. 4 The 1-D S-wave velocity structure along the INDEPTH-III profile

图5 INDEPTH-III剖面下的S波速度剖面Fig. 5 S wave velocity profile under INDEPTH-III profile

图6 INDEPTH-III剖面所穿过的不同构造单元下的和整体的平均S波速度模型Fig. 6 Average shear wave velocity model beneath various tectonic subdivisions along the INDEPTH-III profile and the whole average shear-wave velocity model

INDEPTH-III测线下的伪层析成像S波速度异常结构(图5c)中蓝色部分为S波速度结构速度值比iaspei91模型速度低的区域， 与图 5b展示的INDEPTH-III台站二维速度剖面相比， 图5c更加地突出显示了低速层的分布特征。图5b和图5c， 综合显示了在侧线中广泛分布有低速区， 低速区分布大多集中于 20～40 km 的范围内， 在 32°～32.5°之间壳内低速层出现明显不连续， 刚好位于班公湖—怒江缝合带(吴珍汉等， 2016)之下。沿 INDEPTH-III开展的大地电磁测深结果(魏文博等， 2006)显示在羌塘地体存在相邻两个10～70 km的高导体， 这与本次研究结果(图5)显示的在32.5°～34°之间存在低速区分布特征相似， 即低速高导体。INDEPTH-III宽频带数据计算拉萨地体北部和羌塘地体波速比(Vp/Vs)整体偏高(李永华等， 2006b)同样得出了测线下方普遍存在低速层的结果。在缝合带南部的拉萨地体相邻测线深反射结果显示15～20 km处出现高亮点， 拉萨地体地壳中普遍存在低速层(Nelson et al.， 1996)。沿东经 88.5°剖面下的 S 波速度结构(吴蔚等， 2017)显示羌塘地体普遍存在低速层。通过与不同地球物理方法结果(Nelson et al.， 1996； 贺日政，2003； 魏文博等， 2006； 吴蔚等， 2017)对比， 显示拉萨地体和羌塘地体下存在低速层。

4 讨论

班公湖—怒江缝合带及其两侧的复杂速度结构特征(图 5)表明沿 INDEPTH-III剖面下的拉萨地体(北部)及羌塘地体(南部)物质结构差异明显， 如二者的平均速度结构(图6)、壳内低速层分布特征。

4.1 拉萨地体及羌塘地体平均速度模型

将台站下方速度模型分成拉萨地体、班公怒江缝合带和羌塘地体三个区域获得平均的速度模型(图 6a， b， c)， 与李永华速度模型(2006a)相比， 本文得出的低速层不仅比参考模型速度值低而且出现了速度翻转即低速区比上下区域速度值都要低。班公湖—怒江缝合带下的低速层厚度小于其两侧的拉萨地体和羌塘地体下的低速层厚度(图5)。拉萨地体北部和羌塘地体南部普遍存在 20～40 km 低速层， 细节上有所不同(图 5)， 羌塘地体低速区分布比拉萨地体内低速层较深， 可能与构造变化特征相符， 如中生代晚期就位的拉萨地体北部(潘桂棠等， 2006)、班公湖—怒江缝合带(邹长桥等， 2012； 吴珍汉等，2016； 吴蔚等， 2017)和经过后期新生代岩浆作用(邓万明等， 1996)强烈改造的羌塘地体(吴蔚等， 2017；严江勇等， 2019)。基于三个地体下方速度模型编制了青藏高原中部平均速度模型(图 6d)， 突出显示了INDEPTH-III剖面下的地壳整体上可分为与全球标准模型接近的上地壳(0～20 km)、相对低速的中地壳(20～40 km)及具有强反射特征(Ross et al.， 2004)的相对高速的下地壳(40～65 km)。

4.2 沿INDEPTH-III剖面下的低速层分布特征

INDEPTH-III剖面下S波速度结构反演结果表明低速层在拉萨地体北部和羌塘地体南部是普遍存在的， 如图5c所示在拉萨地体北部存在一个低速区，而羌塘地体南部存在两个低速区， 整体上与其他地球物理研究成果(Brown et al.， 1996； Kind et al.，1996； Nelson et al.， 1996； Wei et al.， 2001； Zhao et al.，2001； 郑洪伟等， 2007； 贺日政等， 2007)相一致。此外， 图5c也显示拉萨地体北测和羌塘地体南侧壳内低速层横向分布并不连通， 即在班公湖—怒江缝合带处出现不连续。INDEPTH-III项目实施的同剖面的大地电磁探测， 显示在地壳内部分布着不连续的高导体(魏文博等， 2006)。与邻近的E88.5°剖面结果(吴蔚等， 2017)对比， 本文结果也显示青藏高原地壳中普遍存在低速层且在羌塘盆地下低速层横向分布并不连续(吴蔚等， 2017)。这表明青藏高原中部地壳内部低速层广泛存在但不连续分布， 推测其动力学成因不仅与寄主地体固有属性(Chang et al.， 1986)有关， 而且也与青藏高原快速隆升(Yin and Harrison， 2000)相关。

地壳速度变化、物质组成与地壳温压状态有关。随着温度升高， 岩石物理状态会发生变化， 如各种干变质岩在1000 MPa(地下30 km)和室温条件下 S波速度平均值为 3.65 km/s， 当温度上升达到900℃以上时干地壳会塑性变化并发生熔融， 其 S波速度会下降至3.45 km/s(Yang et al.， 2012)。在这种条件下， 由公式(4)获得的相对于改进后 iaspei91模型的速度扰动LVZ=(3.45-3.75)/3.75= -8%； 若地壳中存在水时， 地壳的岩石发生熔融所需的温度会更低(即较低的S波速度和相对高的S波扰动)。据此推测， 图5c所展示的INDEPTH-III剖面下的中地壳扰动为-8%以下的蓝色区域发生了部分熔融。

沿INDEPTH-III剖面实施的深地震反射剖面调查(Haines et al.， 2003)显示拉萨和羌塘地体地下的20～30 km 处出现了跳跃的负震相， 亮点反射被解释为因部分熔融引起(Haines et al.， 2003)。Mechie et al.(2004)在深地震测深速度模型(Zhao et al.， 2001)基础上开展了地方震事件波形剖面速度分析， 认为该部分熔融层顶部温度突然增高(＞700℃)造成了石英 α-β 相变(Mechie et al.， 2004)。该低速层也与INDEPTH-III电磁测深剖面中的高导体(魏文博等，2006； 金胜等， 2010)和 Lg 波衰减特征(Fan and Lay，2002)相一致， 被解释为部分熔融层。

印度岩石圈地幔的北向俯冲(郑洪伟等， 2007)所形成的岩石圈增厚， 进而发生的断离或拆沉等构造活动产生多种生热效应(Haines et al.， 2003； 金胜等， 2019)， 使得地壳温度升高， 造成部分熔融(Haines et al.， 2003； Ross et al.， 2004)。现有的地震学研究显示北向俯冲的印度板块岩石圈地幔前缘已经俯冲至班公湖—怒江缝合带之下的 250 km深(Kosarev et al.， 1999； Zhao et al.， 2001； 郑洪伟等，2007； He et al.， 2010； Liang et al.， 2012； Zhang et al.，2012； Wang et al.， 2019)， 印度岩石圈地幔的深俯冲所产生的构造热促使班公湖—怒江缝合带南侧中下地壳温度升高导致了中地壳发生了部分熔融， 地球物理异常表现为如图 5c所示的低速(Fan and Lay，2002； Haines et al.， 2003； Mechie et al.， 2004)、高导(魏文博等， 2006； 金胜等， 2010)和亮点反射(Nelson et al.， 1996； Ross et al.， 2004； 卢占武等， 2014)等。

此外， INDEPTH-III剖面下的0～10 km的上地壳表现为高低速异常体相间， 这与附近的 E88.5°剖面结果(吴蔚等， 2017)相似。依据该区域的地表地质特征(王成善等， 2001； 吴珍汉等， 2016)， 近地表的该低速异常体与INDEPTH-III剖面所经过的班戈盆地和南羌塘盆地等沉积坳陷一致， 而高速异常体与班公湖—怒江缝合带内出露蛇绿岩带分布(吴珍汉等， 2016)有很好的对应关系。

5 结论

本文计算了INDEPTH-III计划中分布在班公怒江缝合带两侧的 23个宽频带地震仪采集的接收函数， 并利用时频域相位滤波(曲中党等， 2015)对数据进行滤波， 用接收函数复谱比非线性方法反演(刘启元等， 1996)获得台站下方的一维 S波速度结构，结合实际地质情况对低速层分布及其成因进行讨论，获得以下认识：

南北位于班公怒江缝合带两侧的拉萨块体和羌塘块体地壳内普遍存在低速层且在缝合带处出现间断， 有些地区还存在两个低速层。上部低速层，埋深较浅， 大多在 0～15 km处， 这可能与羌塘盆地发育有沉积坳陷相关。位于地壳中下部的低速层可能是由于在由青藏高原隆升所致的特定温压条件下岩石发生部分熔融形成。

由于印度岩石圈地幔的俯冲产生的构造热使拉萨地体中下地壳发生部分熔融， 而羌塘地体由于印度板块的俯冲前缘进入地幔使得岩石圈地幔活化，进而地幔热物质上涌使中下地壳的温度升高出现部分熔融。

致谢：本文研究所用的宽频带地震观测数据从www.iris.edu下载。在此特别感谢 INDEPTH-III组野外辛苦工作， 感谢中国地震局地质研究所刘启元老师提供的接收函数复谱比非线性反演程序包， 感谢审稿人的细致审阅和宝贵意见！

Acknowledgements:

This study was supported by National Natural Science Foundation of China (Nos. 41761134094，41574086 and 41274095)， National Key Research &Development Program of China (Nos. 2016YFC0600301 and 2018YFC0604102)， and China Geological Survey(No. DD20190015).