基于多通道互相关的青藏高原东北缘背景噪声程函层析成像

马小军 吴庆举 潘佳铁 钟世军 徐 荟

1)中国地震局地球物理研究所，北京 100081

2)宁夏回族自治区地震局，银川 750001

3)北京市地震局，北京 100080

0 引言

解释青藏高原地壳缩短增厚成因的动力学模型包括：1)刚性块体模型，该模型认为高原下部的物质沿着大型断裂向E挤出(Molnaretal.，1975；Tapponnieretal.，1982)；2)下地壳管道流模型，即青藏高原周缘中下地壳广泛分布塑性通道流(Roydenetal.，1997；Clarketal.，2000)。因此，很多地震学者的研究都是围绕上述模型试图解释高原地壳增厚、壳幔变形的动力学机制(Shapiroetal.，2004；李永华等，2006；Yuanetal.，2013；Shenetal.，2015；Wangetal.，2016)。地震波走时成像是研究地球内部壳幔结构和动力学特征的重要手段，走时成像的分辨率与台站空间分布和射线交叉程度正相关。以往由于台站空间尺度大、分布不均匀且稀少，所得研究结果的精度和分辨率受到了限制。现今，开展密集台阵观测已成为一种研究壳幔精细结构的重要方法。“中国科学台阵探测”2期项目于2013—2016年在南北地震带北段布设了676个密集流动台站，台站之间平均距离为35km，所得观测结果为得到区域更精细的三维波速结构提供了良好的数据支撑。以这些数据为基础，目前已在体波成像(张风雪等，2018)、接收函数成像(王兴臣等，2017；Wangetal.，2017)和面波成像(Lietal.，2017；潘佳铁等，2017；钟世军等，2017；郑晨等，2018；Wangetal.，2020)等方面取得了很多研究成果。

背景噪声成像方法自创建以来，在壳幔结构研究方面得到了广泛应用，如全球尺度(Sageretal.，2020)、区域尺度(Sabraetal.，2005；Shapiroetal.，2005)和小尺度(Brenguieretal.，2007；Mordretetal.，2015)剪切波速度结构研究，通过噪声互相关提取体波并将其用于地幔转换带结构的研究(Polietal.，2012a，b)。传统的噪声面波成像以射线走时为基础，假设面波沿大圆路径传播。实际在横向结构十分不均匀的地区，面波的传播路径存在一定程度的弯曲，且传统的面波射线走时成像往往采用正则化反演方法，正则化过程中的模型平滑参数的选择容易引入虚假异常。

近年来，Lin等(2009)发展了基于程函方程的噪声面波成像方法。程函层析成像方法充分考虑射线的弯曲，通过追踪波前局部相位变化得到波场走时梯度，进而由波场梯度得到慢度矢量。虽然面波程函层析成像具有反演速度快、分辨率高等优点，但要求台站布局足够密集，其分辨率主要由台站孔径和格点尺寸决定。Lin等(2009)提出的噪声程函成像方法首先需提取频散曲线，目前广泛采用的频时分析(Frequency-time analysis，FTAN)相速度频散提取方法(Levshinetal.，1972，1992)是根据高斯窄带通滤波后解析信号的瞬时相位的定义得到的，其中需要相位模糊度参数，而相位模糊度参数由理论或经验相速度频散曲线来确定，无疑将带来系统误差(Bensenetal.，2007；Linetal.，2008；Bouéetal.，2013)。Jin等(2015)发展了基于多通道互相关相位延迟的天然地震面波自动成像方法，通过自动计算邻近所有台站对的相对相位走时延迟，采用程函成像方法反演得到面波频散的二维分布图像。多通道互相关相能够抑制通过相位延迟方法计算得到的面波频散的非相干噪声，并降低因多路径散射波导致的测量偏差。

图1 研究区域的构造背景与台阵分布Fig.1 Tectonic setting and stations in the study area.蓝色三角为科学台阵2期流动台站，灰色实线为主要的构造边界或断裂。HYF 海原断裂；SQS 南祁连缝合带；KF 昆仑断裂；HT 河套-吉兰泰盆地；ALS 阿拉善盆地；OD 鄂尔多斯地块；YCD 银川盆地；QL 祁连造山带； KQ 西秦岭造山带；SG 松潘-甘孜地块

本研究收集了“中国科学台阵探测” 2 期项目的676个流动台阵记录的连续波形数据(图1)，通过噪声互相关计算叠加2015年1—12月的互相关波形，得到ZZ、ZR、RZ、ZT、TZ、RT、TR、TT等9个分量的经验格林函数。本文将Jin等(2015)的方法首次发展应用到背景噪声层析成像中，首先采用多通道互相关计算相位的相对走时延迟；然后基于程函层析成像得到整个研究区域Rayleigh面波8～40s周期的相速度高分辨率分布结果，大部分研究区域的分辨率在40～100km之间；最后对地震面波程函成像结果和传统噪声程函成像结果进行了对比分析，并结合其他研究成果探讨了研究区域的构造变形和动力学过程。

1 数据与方法

1.1 波形互相关计算相位延迟

Jin等(2015)提出的多通道互相关地震面波成像方法，是对同一震源、不同台站接收到的面波波形进行互相关，以自动测量相对相位和波包群的相对变化。对于噪声经验格林函数，我们以一个台站为“震源”，测量以这个台站为虚源的所有台站对的相位延迟。假设以台站1为震源，台站1与台站2、台站3之间的经验格林函数波形分别为S1、S2，则S1和S2之间的加窗互相关函数可表示为含5个参数的小波近似：

F(ω)*WsC(t)≈AGa[σ(t-tg)]cos[ω(t-tp)]

(1)

通过非线性最小二乘拟合可以得到群延迟和相延迟tg和tp。其中，F(ω)为对应频率ω的高斯窄带通滤波器，Ga为高斯函数，A为尺度因子，WsC(t)为加窗互相关函数，σ为半带宽。

1.2 程函层析成像

与传统的基于射线走时的成像方法不同，程函层析成像是从基于射线理论的程函方程发展而来的，标准形式的程函方程为(Shearer，2009)

(2)

(3)

式中，δτp为波场相对相位延迟时间，ri为2个台站之间的路径。

与Lin等(2009)通过求解式(2)中的波场梯度求解慢度空间不同，Jin等(2015)采用最小方差方法反演单个事件的慢度空间分布，反演的目标函数为

(4)

权重因子λ可控制反演的平滑度，δτ为相对相位延迟时间，可对观测波形与参考波形进行互相关计算得到。

1.3 数据处理和质量控制

进行单台数据预处理时，首先按照Bensen等(2007)提出的过程进行处理，即波形截取、去仪器响应、去均值和去趋势，并在频率域进行归一化；然后采用Poli等(2012b)提出的地震窗截断方法，使用长4h的时间窗在每天的连续波形上滑动，并进行50%的重叠，去除每个时间窗内大于3倍标准差的波形；最后计算互相关函数并经线性叠加得到最终的互相关函数。我们通过互相关函数的Hilbert变换计算得到经验格林函数，相较于计算互相关函数的负导数得到的经验格林函数，Hilbert变换具有不改变波形振幅的优点。为了减少噪声源非均匀分布对反演结果的影响，我们对因果信号和非因果信号进行反序平均，得到最终的经验格林函数(Linetal.，2008)。

按照Jin等(2015)的定义，波形互相关系数为

(5)

1.4 相速度计算

本文采用Jin等(2015)提出的程函层析成像方法(式(3))，对以单台为虚源的所有射线的相对相延迟进行反演，得到网格尺寸为0.2°×0.2°的波场相速度分布。

图2 为以台站15581和51555为虚震源的8s周期Rayleigh面波相速度结果，其中箭头为波场走时梯度方向，即波场的传播方向。8s周期面波的相速度对浅层地壳的S波速度较为敏感，受地表构造约束，鄂尔多斯地块西缘、银川盆地、祁连造山带西部和松潘-甘孜地块的东北部表现为低速异常，而西秦岭造山带东部为高速异常。

图2 以台站15581(a)和51555(b)为虚源(紫色五角星)的8s周期相速度及波场梯度方向Fig.2 Phase velocity and gradient direction at the 8s period of Rayleigh wave with the central station 15581(a) and 51555(b) as the virtual source.

图3 台站64016对应的所有射线的走时延迟与距离差的关系Fig.3 Travel time Delay with relative distance difference for the Station 64016.彩色叉为不同周期的走时延迟与均值之差<8s的结果，黑色圆圈为超过8s的结果

图4 平均频散曲线Fig.4 The mean dispersion curve.蓝色实线为利用频时分析方法(FTAN)计算的平均相速度频散，红色实线为本文结果；误差棒为测量的标准差

最终每个周期的相速度结果由所有单台结果加权叠加得到，主要分为3步：首先对所有单台相速度结果进行加权平均叠加得到平均相速度结果；之后对每个单台相速度进行选择，剔除单台走时超过8s的结果；最后对剔除后的结果进行加权平均叠加，得到最终的相速度分布。

1.5 反演分辨率估计与误差分析

程函层析成像反演充分考虑了射线的弯曲，反演的分辨率由台阵的分布、射线交叉程度和网格大小决定(Linetal.，2009；Qiuetal.，2019)。本文采用检测板测试方法给出了0.4°×0.4°、1°×1°尺度异常格点的检测板恢复结果。

图5 8s、25s和40s周期的检测板测试结果Fig.5 Checkerboard test image at 8s，25s and 40s periods.a 网格尺寸为0.4°×0.4°的结果；b 网格尺寸为1°×1°的结果

为检验结果的可靠性，我们还进行了误差分析检验，每个格点相速度的误差由所有单台源相速度的标准误差给出(Linetal.，2009)：

(6)

其中，n为单台源的数目。图6 为8～40s周期的相速度误差，除边缘地区外，研究区域的相速度误差都在40m/s以下，表明反演的精度较高，结果可靠。

图6 8～40s周期的相速度误差Fig.6 Phase velocity uncertainty at periods 8～40s.

2 结果

图7 不同周期的相速度深度敏感核Fig.7 Sensitivity kernel of the phase velocity at different periods.

短周期(8～12s)相速度主要反映了上地壳(深6～15km)的S波速度结构(图8)，上地壳的速度结构与区域地表地貌特征具有很好的相关性，鄂尔多斯盆地西缘、银川地堑、河套-吉兰泰盆地、阿拉善盆地南部表现为显著的低速异常，另外祁连造山带西部、松潘-甘孜地块东北部也显示出低速异常(约3.1km/s)。高速异常地区与造山带、山脉隆起的分布一致，如西秦岭造山带、祁连造山带东部呈现高速特征。以上结果与前人的研究一致(Shenetal.，2016；Lietal.，2017；潘佳铁等，2017；钟世军等，2017)。

中等周期(18～22s)的相速度主要对中地壳(深15～30km)的S波速度结构敏感(图8)。祁连造山带西部、松潘-甘孜地块东北部表现出明显的低速异常(约3.1km/s)，而河套-吉兰泰盆地呈现低速异常特征，随着周期增大低速异常逐渐减弱；阿拉善盆地南部呈现微弱的低速异常；蒙古褶皱系、祁连造山带东部、鄂尔多斯盆地西缘和西秦岭造山带为高速特征。

长周期(25～40s)相速度主要对下地壳—上地幔顶部的S波速度结构敏感(图9)，下地壳速度结构的横向分布与莫霍面的起伏相关(Wangetal.，2017；王兴臣等，2017)。青藏高原内部与其周缘地区的莫霍面深度具有显著的差异性，青藏地块内部的祁连造山带西部和松潘-甘孜地块东北部表现为显著的低速异常特征，而青藏地块周缘的阿拉善盆地、鄂尔多斯盆地西缘、西秦岭造山带等大范围区域显示出明显的整体一致的高速特征，但河套-吉兰泰盆地与周缘相比显示出显著的低速异常。先前的面波成像研究也显示同样的结果(Lietal.，2017；潘佳铁等，2017；钟世军等，2017；杨志高等，2019)。

图8 8s、12s、18s和22s周期的相速度Fig.8 Phase velocity at 8s，12s，18s and 22s periods.

图9 25s、30s、35s和40s周期的相速度Fig.9 Phase velocity at 25s，30s，35s and 40s periods.

3 讨论

3.1 与他人结果的对比

程函层析成像的优势在于能够反演局部小尺度的结构特征，具有较高的横向分辨率。Jin等(2015)创立的多通道互相关面波自动成像方法主要用于长周期(>20s)天然地震面波成像。钟世军等(2017)利用科学台阵2期天然地震面波数据，同样采用多通道互相关面波自动成像方法(Jinetal.，2015)得到了长周期地震面波的相速度频散。本文将该方法应用到短周期的噪声面波成像中。我们将本文的相速度纯路径频散结果与钟世军等(2017)的地震面波成像进行了对比。

图10 本文结果(a)与钟世军等(2017)结果(b)的差异Fig.10 The phase velocity difference at different periods between this study (a) and the result from ZHONG Shi-jun et al.，2017 (b).

图10 为本文与钟世军等(2017)结果的比较(14s、20s、40s周期)，左列(图10a)为本文的结果，右列(图10b)为钟世军等(2017)的结果。从图10 可以看出，对于14s、20s周期二者在大部分研究区域的一致性较好，相速度差异在误差范围之内，均可刻画局部小尺度速度结构的细节，本文结果的分辨率为20～40km，特别是短周期(8～25s)的分辨率可达20～30km。但祁连山西部和西秦岭的14s周期的结果具有较大的差异，钟世军等(2017)给出的祁连山西部的速度值则更低。由于天然地震的短周期面波射线密度较低，因此频散的可靠性较低。

在40s周期，两者在大尺度上特征一致，但在小区域上的差异较大，天然地震面波成像结果具有更高的分辨率，推测主要原因在于两者使用的数据存在差异，背景噪声源的分布特征决定了长周期的面波频散数目较少，而地震面波长周期频散数目更多，因此结果具有更高的分辨率。

3.2 与传统噪声程函层析成像方法的对比

Lin等(2009)提出的传统噪声程函层析成像方法首先需要提取面波频散曲线，然后再进行程函层析成像。目前，可使用自动FTAN方法提取Rayleigh波相速度频散曲线(Levshinetal.，1972，1992)，利用窄带通高斯滤波器对噪声互相关函数的包络信号进行滤波，测量中心频率对应的最大振幅和到时，进而提取面波群速度频散，通过相慢度、瞬时相位、群速度和相位模糊度参数的关系式导出相速度。其中，需要由参考模型频散曲线或已有频散得到相位模糊度参数(Bensenetal.，2007；Linetal.，2008)。

图11 本文方法反演的结果与Lin等(2009)的传统程函成像方法结果的对比(12s和30s周期)Fig.11 Comparison of phase velocity maps for periods of 12s and 30s.a 本文计算结果；b 传统程函成像结果；c 二者之差(图a减去图b的结果)

3.3 低速层成因的简单解释

3.3.1 河套-吉兰泰盆地的低速层

地壳浅部的低速层往往与其上的沉积盖层密切相关。研究区内的银川盆地、鄂尔多斯盆地西缘和河套-吉兰泰盆地的上地壳存在巨厚沉积层，河套-吉兰泰盆地的上地壳存在显著的低速异常并持续到18s周期。已有的地质研究表明，河套-吉兰泰盆地是新生代以来形成的大型断陷盆地，新生界最大沉积层可厚达10km以上(郑孟林等，2006)。而长周期(35～40s)相速度(相当于40～80km深度)结果表明河套-吉兰泰盆地的下地壳、上地幔也存在弱低速层。Chen等(2009)的S波接收函数研究和He等(2017)的面波成像研究表明河套-吉兰泰盆地的岩石圈厚度只有80km，远低于其南部的鄂尔多斯盆地的岩石圈厚度(160km)。同时，杨嵩等(2013)对华北地区的岩石圈厚度及上地幔温度进行了研究，发现河套-吉兰泰盆地下方的温度明显高于其周缘地区。本文推测，“大地幔楔”中的地幔软流圈热物质上涌可能是下地壳及岩石圈地幔出现弱低速异常的主要原因(郑晨等，2018；Leietal.，2019)。

3.3.2 松潘-甘孜地块和祁连山西部的低速层

有研究表明，地壳的径向各向异性是云母晶体的定向排列所致。正的径向各向异性与扩张型地壳(即水平向物质流)对应；反之，负的径向各向异性对应垂向的增厚。而松潘-甘孜地块东北部存在弱的负径向各向异性特征(汪洋等，2001；Baoetal.，2015；杨志高等，2019)，表明不太可能存在水平向的下地壳管道流。

4 结论

本文基于青藏高原东北缘及邻区密集台阵记录的连续波形得到了背景噪声经验格林函数，进而发展了基于多通道互相关的面波程函层析成像方法，首次将其应用于背景噪声层析成像中，反演得到了研究区域不同周期的相速度分布图像。

经对比，本文所得结果与已有的天然地震程函层析成像结果一致，表明本文结果是可靠的，且具有高分辨率和高精度的特点。与传统程函层析成像方法的研究结果对比，本文结果降低了长周期信号由于信噪比的降低和有效射线数目减少所带来的误差，压制了非相干多路径散射面波引起的测量误差，提高了反演的精度和稳定性。

反演得到的8～40s周期的相速度分布图像表明，河套-吉兰泰盆地的下地壳、上地幔显示出弱低速特征，地幔高温、软弱物质的侵入可能造成了该区域下地壳、上地幔出现弱低速。松潘-甘孜地块东北部中下地壳存在低速体异常，可能与下地壳的部分熔融有关。祁连造山带西段中下地壳广泛存在低速层，分析认为，祁连山西段中下地壳可能不具备部分熔融或地壳管道流，低速层可能与其地壳在外力作用下的缩短、增厚和解耦有关。

致谢中国地震局地球物理研究所“中国地震科学台阵探测中心”为本文提供了连续地震波形数据；Jin 和Gaherty博士为本研究提供了面波自动成像程序(ASWMS)(1)http：∥www.iris.edu/ds/products/aswms/。，冯力理博士为本研究提供了传统程函层析成像程序(2)https：∥github.com/NoisyLeon。；审稿专家为本文提供了很好的意见和建议。在此一并表示感谢！