青藏高原东北缘瑞利面波成像：射线理论与程函方程结果的比较

郝识杰，黄周传，王良书，徐鸣洁，米 宁，于大勇

南京大学 地球科学与工程学院，南京 210023

1 引言

青藏高原是印度板块和欧亚板块碰撞形成的年轻的高原（图1；Yin and Harrison, 2000）。对于青藏高原的隆升机制存在热烈的讨论，前人提出许多端元模型，如地壳流模型（Clark and Royden,2000; Royden et al., 2008），刚性块体挤出模型（Tapponnier et al., 2001），岩石圈连续变形模型（Molnar and Tapponnier, 1975）等。获取精细的岩石圈地震波速度结构，对于讨论青藏高原的构造演化过程，具有重要的意义，广泛布设的密集地震观测台阵为获取更高精度的地球内部结构提供了可能。

“中国地震科学台阵南北地震带北段”观测项目（ChinArray 二期），在青藏高原东北缘及其相邻区域布设了675台流动地震仪。面波速度成像是利用密集台阵数据，对岩石圈结构进行约束的有效方法。传统的面波成像研究多基于射线假设，即假设地震台接收到的面波信号是沿震源与台站之间的大圆弧路径传播（Barmin et al., 2001）。在实际情况下，由于岩石圈波速结构存在横向变化，面波的传播方向会偏离大圆弧路径，导致基于射线假设的方法出现误差。Foster等（2014）利用类似聚束分析的方法分析了USArray台阵记录的地震事件的入射角度，发现入射角度的偏差可以达到15°，由此造成的相速度误差可以达到4%。Lin等（2009）提出了基于程函方程的成像方法，该方法首先计算面波的相位走时场，然后基于程函方程计算相速度的大小和方向。在该方法中，面波速度的方向从地震记录计算得来，理论上可以减少因面波传播偏离大圆弧路径导致的误差。在USArray数据集，基于程函方程的瑞利波相速度成像结果和Foster等（2014）矫正入射角度之后的成像结果有很好的一致性（Jin and Gaherty, 2015）。这一结果表明基于程函方程的成像方法可以有效地修正面波传播方向偏离大圆弧路径导致的误差。

在青藏高原东北缘，前人已经开展了许多基于面波的工作。如：基于射线的面波成像（潘佳铁等, 2017），基于程函方程的面波成像（钟世军等, 2017），面波反演三维S波速度结构（Li et al.,2017），面波与接收函数联合反演三维S波速度结构（Wang et al., 2017）。对于这些基于面波的研究工作而言，面波速度成像都是重要的基础。在这些研究中，面波速度成像的结果有相似的总体特征（青藏高原东北缘存在低速异常，鄂尔多斯块体西部上地壳存在低速异常，中下地壳和岩石圈表现出高速异常），在细节上有所区别（银川河套地堑在地壳是否存在低速异常）。由于这些研究使用了同一个数据集，所以成像结果的区别应该与使用的成像方法有关。

本研究分别使用基于射线的成像方法（Debayle and Sambridge, 2004）和基于程函方程的成像方法（Jin and Gaherty, 2015）计算青藏高原东北缘的相速度分布。通过对两种方法原理的阐述和两种方法计算结果的对比，讨论前人研究中存在争议的速度结构细节。

图1 青藏高原及周边地区地形图（红色矩形标识研究区域）Fig.1 Topography map of the Tibetan Plateau and adjacent region (red rectangle denotes the study region)

2 数据和方法

2.1 数据

本研究中使用的地震数据来自ChinArray二期项目，包含675台宽频地震仪（图2b），平均台间距约40 km，记录时间从2013年9月到2015年5月，观测的采样频率为100 Hz。两种成像方法使用了相同的天然地震瑞利波数据。为了挑选较为清晰的基阶面波波形，笔者选择了震中距在10°～160°范围内、矩震级大于5.5、震源深度小于50 km的地震事件，共351个（图2a）。

2.2 计算台站之间的相位走时差

瑞利波相速度成像的第一步是利用同一个地震事件的记录，计算两个台站之间的瑞利波相位走时差。为了客观地比较两种成像方法，使用同一套流程计算相位走时差（Jin and Gaherty, 2015），然后分别用于两种成像方法。首先进行预处理，通过群延迟时间计算窗口函数Ws，分离出基阶面波的信号。然后通过如下公式计算互相关函数C(t):

其中S1和S2分别代表两个台站记录的地震波形。对计算得到的互相关函数进行窄带通滤波，通过拟合如下公式得到相位走时差等参数（Gee and Jordan, 1992）：

等式左侧表示对互相关函数用滤波器Fi进行窄带通滤波，并用窗口函数Wc突出主要能量。等式右侧Ga表示高斯函数，用来拟合的五个参数分别为A、σ、ω、tg和tp，分别表示幅值参数，高斯窗口半宽度，窄带通滤波的中心频率，群延迟时间和相位延迟时间。利用拟合得到的相位走时差，可以进行下一步的成像工作。因为拟合的曲线是余弦曲线，拟合时可能产生误差，导致结果跳过一个或多个波长的范围。笔者选择间距200 km以内的台站对计算相位走时，这样面波波列较短，拟合更准确。

2.3 基于射线的面波相速度成像方法

如果震源和两个台站近似处于同一个大圆弧路径上，那么两个台站的相位走时差可以用于计算两个台站之间的平均相速度。在上一步计算的相位走时差中，挑选符合条件的数据，计算了台站对之间的平均相速度，射线密度如图3。在得到射线上的平均相速度之后，按照Debayle和Sambridge（2004）提出的面波成像流程进行相速度成像，得到二维的相速度结构。基于射线理论，平均速度可以表述为如下公式：

图2 （a）使用的地震事件分布图（b）本研究使用的青藏高原东北缘台站分布图（据Taylor and Yin, 2009修改）Fig.2 (a)Seismic events used in this research (b)Stations used in this research

图3 基于射线的成像方法中不同周期射线密度分布图Fig.3 The density of ray paths for tomography based on ray theory

其中T表示周期，c(T,rˆ)是面波的传播速度，rˆ表示地球表面的位置矢量，p表示面波的传播路径，cp(T)表示路径上的平均相速度，Lp表示路径的长度。面波传播路径近似为震源和台站之间的大圆弧路径，这一假设成立的条件是面波速度横向不均一性足够小，对射线路径造成的扰动小于想要达到的分辨率。通过阻尼线性反演，使模型正演得到的两点之间平均相速度逼近实际观测的台站之间的平均相速度。

2.4 基于程函方程的面波相速度成像方法

台阵记录到的，地震事件激发的单一频率的面波波场，可以由程函方程表示：

其中δτ为两个台站之间的相位走时差，是慢度矢量，表示两个台站之间的路径。同样地，通过阻尼线性反演，使模型正演的相位走时差逼近观测的相位走时差。由于基于程函方程的成像方法计算的是慢度（速度）矢量，所以任意两个台站对之间的相位走时差数据都能参与反演，从而更充分地利用数据。

3 结果

瑞利波相速度主要受岩石圈的P波速度，S波速度和密度影响，对S波速度变化敏感性最强，所以可以用来约束岩石圈S波速度结构。总体上，周期较长的相速度对更深部速度结构变化更敏感，速度结构变化对特定周期的相速度的影响可以通过敏感核量化地表示。基于LITHO1.0模型（Pasyanos et al., 2014）和ak135模型（Kennett et al., 1995）设置模型，地壳分为三层，莫霍面深度约为50 km，计算了20 s、30 s、60 s的瑞利波相速度敏感核，作为分析讨论成像结果的依据（图4）。从敏感核可以看出，20 s和30 s的瑞利波相速度主要反映中下地壳的结构，60 s的结果主要反映上地幔顶部的结构。

两种方法在20 s、30 s、60 s的成像结果如图5所示，二者速度异常特征存在很好的一致性，且与主要构造单元的分布相关性较好。在青藏高原东北缘，两种方法在20 s、30 s、60 s的结果都表现为低速异常，且低速异常分布范围较为吻合。在20 s和30 s的结果中，松潘甘孜块体和祁连块体所在位置存在低速异常，且两块低速异常体互相分开。这一细节特征在两种方法的结果中都有体现，根据敏感核的计算结果，主要反映中下地壳的S波速度结构。用背景噪声面波反演得到的三维S波速度模型中，中下地壳松潘甘孜块体和祁连块体的低速异常表现为互相分离的两块（Bao et al., 2013; Li et al., 2014）。在60 s的相速度图像中，青藏高原东北缘又整体呈现为低速异常体。

图4 (a)计算敏感核使用的速度模型 (b)瑞利波相速度敏感核Fig.4 (a)Velocity model used in calculation of the kernel (b) sensitivity kernel of rayleigh wave phase velocity

在鄂尔多斯块体，相速度图像表现出显著的高速异常，这一特征在两种方法中都有体现，反映了鄂尔多斯块体作为稳定克拉通的特征。在鄂尔多斯块体西北部的银川河套地堑，两种方法的结果在20 s、30 s都表现出低速异常，在60 s基于程函方程的结果有低速异常，基于射线的结果射线较稀，约束效果不理想。基于程函方程的成像结果中，低速异常体的分布与银川河套地堑的轮廓更为吻合。

4 讨论

4.1 两种成像方法的对比

分别使用基于射线的成像方法和基于程函方程的成像方法得到了青藏高原东北缘的瑞利波相速度成像结果，两种方法使用了同一数据集计算的相位走时差。基于射线的成像方法只使用了与地震事件近似处于同一大圆弧的台站对之间的相位走时差，而基于程函方程的成像方法不受这一限制，利用数据更加充分。且后者同时反演面波速度的大小和方向，可以避免因偏离大圆弧路径传播导致的误差。从结果上看，两种方法取得了相似的相速度图像，但有一些细节区别。为了对比两种方法计算的瑞利波相速度差异，笔者用射线成像方法的结果减去程函方程成像方法的结果（图6）。两种方法得到的速度结构主要在研究区边缘及速度异常体边缘存在差异。笔者推测这是由于基于射线的成像利用的数据有限，对异常体形态约束不如程函方程方法充分，且在研究区边缘射线密度较小，表现出接近于平均速度的相速度特征。因此笔者认为对于研究区边缘的细节结构，基于程函方程的成像结果更加可靠。

4.2 青藏高原东北缘的低速异常结构

面波速度主要受S波速度影响，S波低速异常可能由于介质温度高、力学性质弱或包含流体等因素造成。此次研究的成像结果中，20 s和30 s瑞利波相速度在松潘甘孜块体和祁连块体表现出显著低速异常，且两块低速异常体互不相连。结合相速度敏感核，这反映了该区域中下地壳存在S波速度低速异常。

地壳流模型是解释青藏高原变形机制的一个端元模型，即青藏高原中下地壳较为软弱，在地形高程梯度形成的压强差作用下，由高原中心向高原边缘发生塑性流动（Clark and Royden, 2000;Royden et al., 2008）。地壳流模型解释了青藏高原内部平坦的地形特征和青藏高原边缘陡峭和平缓两种类型的地形变化（Clark and Royden, 2000）。青藏高原通过地壳流机制进行隆升和扩张的前提之一是存在软弱、可以流动的中下地壳。发生部分熔融可能改变中下地壳的流变性质，导致其可以流动。部分熔融会产生地球物理手段可以观测的异常现象。在松潘甘孜块体，接收函数和面波联合反演显示中下地壳存在较高的P波S波波速比（Ye et al., 2017），电阻率反演显示中下地壳存在低电阻率层（Wang et al., 2014），Lg波Q值成像反映出低Q值（Zhao et al., 2013）。这些地球物理研究结果支持松潘甘孜块体的中下地壳存在部分熔融。本研究结果中的松潘甘孜块体相速度低速异常反映了中下地壳S波低速异常，这些S波低速异常可能与中下地壳发生部分熔融有关。由于存在部分熔融，在重力作用下，松潘甘孜块体有产生地壳流的条件。同时，各向异性研究也在地壳观测到与地形梯度平行的快波方向，符合地壳流模型的预测（Huang et al., 2017）。

图5 青藏高原东北缘瑞利波相速度分布图Fig.5 Rayleigh wave phase velocity maps of the Northeastern Tibetan Plateau

图6 射线成像方法与程函方程成像方法结果的对比Fig.6 Comparison between phase velocities which are calculated based on ray theory and Eikonal equation

成像结果在祁连块体表现出和松潘甘孜块体相似的相速度低速异常，但是在祁连块体支持中下地壳存在部分熔融的证据较弱，该块体的纵横波速比低于松潘甘孜地体（Tian and Zhang, 2013），但电阻率高于松潘甘孜地体（Xiao et al., 2013）。祁连块体可能发育逆冲推覆断层，地壳增厚，成为含水流体运移的通道（Xiao et al., 2012），造成低速异常。

4.3 银川河套地堑的岩石圈结构

银川河套地堑位于鄂尔多斯块体的西北部，在印度板块和欧亚板块碰撞之后形成（Zhang et al., 1998）。基于前文的讨论，笔者认为基于程函方程的成像能更好地约束这一区域。用该方法得到的相速度分布图在20 s、30 s和60 s，在银川河套地堑表现出低速异常，且分布范围与地堑的轮廓较为一致，表明本区域中下地壳和上地幔顶部存在S波速度低速异常。体波成像也支持上地幔存在低速异常（Huang and Zhao, 2006）。接收函数研究显示银川河套地堑下方岩石圈厚度（约80 km）比相邻的鄂尔多斯块体（约200 km）薄（Chen et al., 2009）。上地幔顶部的低速异常可能是由于软流圈物质上涌造成的（Tian et al.,2011）。太平洋板块俯冲到华北克拉通下方，可能在上地幔产生了地幔楔（Zhao et al., 2007）。软流圈上涌可能与板块脱水和地幔楔中的拐角流有关（Huang and Zhao, 2006）。中下地壳存在低速异常表明软流圈上涌的影响到达了地壳，同时表现出低电阻率（Dong et al., 2014）和高的纵横波速比（Tian et al., 2011）。

5 总结

本研究使用ChinArray二期共675台宽频流动地震仪记录的远震瑞利波数据，分别使用基于射线的面波成像方法和基于程函方程的面波成像方法计算了瑞利波相速度分布图，揭示了地壳和上地幔顶部的结构横向变化特征。并通过两种方法的对比，讨论了前人研究中存在争议的细部结构特征。

在青藏高原东北缘，祁连块体和松潘甘孜块体在20 s、30 s、60 s存在低速异常。这一速度特征在两种方法取得的结果中都显著存在，反映了地壳较厚且存在S波速度低速异常。在20 s、30 s的相速度分布图中，祁连块体和松潘甘孜块体的低速异常相互分开，在60 s的结果中，两块低速异常体不再显著分开。这一细节特征在两种方法的计算结果中较为一致。支持两个块体存在不同的变形机制，松潘甘孜块体可能存在地壳流，祁连块体可能通过逆冲推覆断层发生增厚。在20 s、30 s、60 s，鄂尔多斯块体存在显著的高速异常，反映了稳定克拉通的特征。

在银川河套地堑，基于程函方程的成像结果在20 s和30 s显示出显著的低速异常，且分布范围与地堑的范围较为一致。基于射线的成像也显示出低速异常，但是低速体分布范围与构造体吻合较差。这可能由于基于程函方程的成像方法更充分地利用了台站之间地相位走时差数据，对较小的速度异常体约束更加准确。

致谢：感谢“中国地震科学探测台阵数据中心”为本文提供观测资料。本研究由自然科学基金项目（41674044；41674049）以及中国地震科学台阵探测项目（DQJB16A0306）资助，黄周传受南京大学“登峰人才支持计划”资助。文中多数图像是由GMT绘制，在此表示衷心感谢。