基于背景噪声面波的吉林省东南部S波速度结构成像

2021-08-10董向欣安百州毕凤仪于平

世界地质 2021年2期

董向欣，安百州，毕凤仪，于平

吉林大学地球探测科学与技术学院，长春 130026

0 引言

吉林省地处东北亚地理中心位置，地貌形态差异明显，可分为东南部山地和西北部平原两大地貌。其地质构造较为复杂，分布有郯庐断裂带北段主干的依兰—伊通断裂和郯庐断裂带北段分支的敦化—密山断裂，松辽盆地以及长白山火山、龙岗火山等全新世活动火山。太平洋板块向西俯冲于欧亚板块下，俯冲带在上地幔的水平运移影响了长白山火山、镜泊湖火山的火山活动[1]。位于吉林省东南部、中朝边境的长白山火山属于巨型的板内层状复式火山[2]，由数以百计的火山锥和广袤的熔岩台地组成了总面积达12 000 km2的庞大火山群，已有多次大规模喷发的地质历史，今后可能有中等规模的喷发危险[3]。位于吉林省靖宇县与辉南县交界处、长白山脉西北侧的龙岗山脉，是长白山脉的支脉，由160余座火山组成了龙岗火山群，是中国最典型的单成因火山作用发育特征的地区之一[4]。龙岗火山群的火山活动始于新近纪，鼎盛于第四纪，休眠于距今约1 600 a±。龙岗火山群和长白山火山群分别展现了单成因和多成因的不同火山成因类型。20世纪90年代以来，前人在吉林省东南部地区利用不同的地球物理学方法技术开展了探测研究，包括人工地震[5--6]、天然地震接收函数[7--10]、大地电磁测深[11--15]、背景噪声成像等[16--21]，显示在长白山火山下方存在低速异常、高泊松比异常、低阻异常和重力负异常等特征，但对于长白山火山地壳岩浆房的空间分布位置和部分熔融程度等地质学问题的认识仍存在较大分歧。近年来，王武等[17]和范兴利等[18]在长白山地区的背景噪声成像研究，探讨了长白山火山岩浆房的位置和地壳熔融程度等争议问题，但所使用的成像方法均为需要生成群速度或相速度图的中间步骤的两步法，具有一定的局限性。

利用多个地震台站的背景噪声记录计算互相关函数并进行长时间叠加，获得面波的经验格林函数并提取群速度或相速度的频散曲线，最终进行成像获得地下速度结构，是近年来发展快速且应用广泛的地震成像方法。背景噪声的理论可以追溯到Aki[22]的文献。背景噪声成像方法同基于地震面波层析成像的传统方法相比，具有射线分布均匀、不需震源信息、不受短周期面波衰减影响和分辨率高等优点，更适用于地震活动较少的地区。Shapiro[23]最早使用背景噪声面波信号反演得到美国加州的高分辨率瑞利面波群速度分布图，随后该方法逐渐应用于全球多个地区不同尺度的地壳速度的研究。Fang et al.[24]提出的基于射线追踪的面波走时直接反演横波速度结构的方法，考虑了复杂介质下的射线路径弯曲的影响，减小了计算误差，有利于较小尺度的地质构造复杂地区的面波反演。为研究吉林省东南部地区S波速度结构，验证直接反演方法在多火山发育区的适用性，本研究使用国家测震台网的吉林、黑龙江和辽宁台网的连续波形数据，利用台站对互相关函数获得瑞利面波的经验格林函数，提取群速度和相速度频散曲线，用面波直接成像方法反演得到吉林省东南部深度为0～30 km的S波速度结构图像。

1 数据预处理及频散曲线提取

选取国家测震台网数据备份中心[25]共39个固定地震台站从2019年1月1日至2020年12月31日的连续地震背景噪声数据。其中吉林测震台网26个台站、黑龙江测震台网7个台站、辽宁测震台网6个台站。研究区内的主要地质构造和台站位置分布如图1所示。仪器类型为长周期仪器，波形采样率为100 Hz，实际数据处理使用了垂直分量的地震记录。

图1 研究区主要地质构造和台站分布图Fig.1 Distribution map of main geological structures and stations in study area

1.1 数据预处理

数据预处理采用了Bensen et al.[26]介绍的由背景噪声数据测量瑞利面波群速度和相速度频散曲线的方法和处理流程。首先对所有地震台站的单台数据按天截取为相同长度，重采样为1 Hz。为降低环境噪声和仪器噪声及人为干扰造成的波形记录倾斜和零飘现象的影响，对重采样后的数据做了去平均值和去线性趋势处理。然后去除仪器响应，目的是消除不同地震仪器固有工作频率的干扰。带通滤波周期范围为3～40 s。最后做了时域归一化和频谱白化处理，时域归一化的目的是除去波形记录中能量较强的地震和爆破等干扰，频谱白化的目的是抑制单频信号和微震的干扰，拓宽背景噪声的频带。

对任意一个台站每天的垂直分量波形数据进行互相关运算，选取时间长度为-200～200 s。将该台站对24个月的单日互相关运算结果线性叠加，生成最终的互相关函数。假定有台站A和台站B，波场由台站A激发，被台站B接收，台站A和台站B记录的背景噪声的互相关函数对时间求导可获得台站B的经验格林函数。

如果在空间上的噪声源是均匀分布的，那么背景噪声记录的互相关函数的正向和负向分支应该是严格对称的，这样求导得到的经验格林函数与精确的经验格林函数具有相同的相位信息。由于实际地球介质中噪声源的不均匀空间分布和随机性，不能保证由实际背景噪声记录的互相关函数的正向和负向分支的对称性。因此，用背景噪声记录的互相关函数的正向和负向分支反序平均的结果求导计算经验格林函数。这同时可以提高计算结果的信噪比。图2是所有台站从2019年1月叠加至2020年12月(共24个月)得到的经验格林函数。图中的结果显示了清晰的面波信号，具有较高的对称性和信噪比。

红色虚线之间为群速度为2～5 km的瑞利面波信号范围，紫色虚线外侧为噪声范围。图2 所有台站经验格林函数Fig.2 Empirical Green function of all stations

1.2 群速度和相速度频散曲线提取

本研究采用姚华建等[27]提出的一种基于图像分析的双台面波频散曲线快速提取方法和由此编写的Matlab交互式处理程序(EGFAnalysisTimeFreq_version_2015)，该程序设计了一个图形用户界面(GUI)，可快速准确地从经验格林函数或互相关函数中提取出双台瑞利面波群速度和相速度频散曲线。笔者以辽宁地震台网的西丰台(LN.XFN)和吉林地震台网的榆树台(JL.YST)为例，介绍了频散曲线提取的主要过程。

如图3所示，生成西丰台和榆树台的瑞利面波群速度频谱能量图。为满足远场假设，获得理论上真实可靠的频散曲线，应用了台站间沿大圆路径的距离至少为2倍瑞利面波波长的质量控制标准。图3中的红色直线表示台站间沿大圆路径的距离等于2倍波长时对应的周期，红色的实心圆表示满足台站间距不小于2倍波长的的质量控制标准的频散点。

同时，为提高测量精度，获得尽可能准确的频散曲线，应用了满足测量信噪比(SNR)的质量控制标准。信噪比的定义参考Bensen et al.[26]，为一定周期或频率的面波信号时窗内包络的最大振幅与一段时间后的噪声时窗内包络的振幅均方根之比。本研究的噪声时窗在信号时窗之后的100 s选取。图3中的蓝色圆环表示信噪比至少为5的瑞利面波群速度和相速度频散点。

a.瑞利面波群速度频散曲线测量图，黑色曲线为各周期的信号最大振幅的连线，红色圆点表示各周期的信号最大振幅所在的位置，蓝色圆环表示提取的满足信噪比要求的频散点，红色直线表示台站间沿大圆路径的距离等于2倍波长时对应的周期；b.瑞利面波相速度频散曲线测量图，黑色曲线、红色圆点、蓝色圆环和红色直线的意义与图3a中的相同。紫色曲线为图3a中提取的群速度频散曲线，绿色曲线为根据保留的相速度频散点计算得到的群速度频散曲线。图3 瑞利面波群速度和相速度频散曲线测量示例Fig.3 Example of Rayleigh wave group velocity and phase velocity dispersion curve measurement

根据以上两个质量控制标准，舍弃了小于2倍波长的台站间距对应的周期和信噪比小于5的周期对应的频散点，只保留了图3a中红色实心圆和蓝色圆环同时所在的群速度频散点，以保证群速度测量的可靠性。

然后，基于已保留的群速度频散点，提取相速度频散曲线。在相速度频谱能量图(图3b)中存在多个相速度频散的分支，只选择离群速度频散曲线最接近且能量最强的分支用来提取相速度频散曲线。同样地，只保留图3b中红色实心圆和蓝色圆环同时所在的相速度频散点。紫色曲线表示上一步提取的频散曲线，绿色曲线表示根据保留的相速度频散点计算得到的群速度频散曲线，两者较为接近，说明根据已保留的群速度频散点提取的相速度频散曲线可靠。

1.3 射线路径分布

本研究共选取研究区域内39个台站，按照n(n-1)/2 公式计算，理论上可以得到741条频散曲线。经过1.2节的数据筛选，实际反演中拾取674条瑞利面波相速度的混合路径频散曲线进行面波直接成像，射线路径分布如图4所示。

图4 射线分布图Fig.4 Ray distribution map

2 S波速度结构反演

2.1 反演方法

采用Fang et al.[24]提出的基于射线追踪的三维浅层地壳结构的面波频散直接反演方法，反演吉林省东南部地区三维地下S波速度结构。相比于传统的两步法，如Barmin et al.[28]提出的面波成像方法，该直接反演方法不需要生成群速度图或相速度图的步骤，使用基于频率相关的射线追踪和基于小波的稀疏约束反演。射线追踪使用Rawlinson et al.[29]提出的快速行进法，计算每个周期的面波旅行时和源与接收器之间的射线路径，避免了在大多数面波层析研究中使用的大圆传播的假设。基于小波的稀疏约束反演[30]考虑了具有平滑变化地震特性的准层状介质，通过网格点下的一维剖面来表示三维剪切波速度模型，这些网格点是使用基于小波的稀疏约束成像方法从所有频散数据中同时确定的。波速模型的小波系数用迭代加权最小二乘算法来估计，并且在迭代时，使用新获得的波速模型来更新面波射线路径和数据灵敏度矩阵。

2.2 检测板测试

为检查反演参数设置的合理性，评价反演模型的分辨能力，通常使用模型的检测板恢复性测试，对台站分布与射线路径覆盖程度的影响进行检测。设置的检测板棋盘模型深度范围为地表至62.5 km，沿深度方向划分为26层，每层厚度为2.5 km，每层的S波速度为1.4 km+深度×0.05 km，根据检测板分辨率测试的原理，使用正弦函数在每层的S波相同速度值上添加了一定的速度扰动量，速度扰动量的幅度范围为每层的S波相同速度值±0.4 km。如图5所示，棋盘模型的平面沿经度和纬度方向划分成15×14共210个网格点，网格大小为0.5°×0.5°。在正演时，先按照实际射线路径计算理论走时，由于受噪声源分布不均匀因素造成的相速度频散误差约为1%[31]，在理论数据中加入2%的随机误差生成合成数据，最后对合成数据反演，获得不同深度的检测板测试结果。如图5的红色矩形所示，反演结果的范围比模型的范围分别在经纬度和深度方向上少1°和1层深度。

图5给出了初始棋盘模型及不同深度的测试结果。检测结果显示深度为5～30 km时，反演的棋盘图案基本重现了给定的初始模型和速度扰动，特别对于研究区域的中心位置，成像结果的可信度较高。但模型边缘重建结果相对较差，这与缺乏射线路径覆盖有关。随着深度的增加，反演重建效果逐渐变差，在深度为35～40 km的图像中表现得较为明显，这是由于长周期的瑞利面波的相速度对相应深度的S波速度的敏感性降低导致。

2.3 反演和讨论

一般来说，瑞利面波的相速度对其三分之一个波长左右对应的地下深度处的S波速度结构最为敏感。笔者的反演初始速度模型根据拾取的瑞利面波相速度频散点一维平均S波速度模型设置，深度范围为地表至62.5 km，沿深度方向划分为26层，每层厚度为2.5 km。模型的平面略大于研究区平面，沿经度和纬度方向划分成15×14共210个网格点，网格大小为0.5°×0.5°。使用拾取的符合质量控制标准的相速度频散作为输入数据，经过反复试验，选取了合适的平滑参数和阻尼参数，迭代次数为10次，反演得到了吉林省东南部地区三维地下S波速度结构。

根据检测板可靠性分析结果，选取研究区深度为5 km、10 km、15 km、20 km、25 km、30 km的S波速度分布图像进行讨论(图6)。

由图6可以看出，研究区存在造山带、沉积盆地、断裂带及火山等多种地质构造，同一深度的S波速度值不均匀分布，且速度异常变化范围较大，表明吉林省东南部地区中上地壳S波速度结构存在较强的横向不均匀性。5 km的成像结果与地表地质构造相对应，速度异常呈现出受郯庐断裂带北延的依兰—伊通断裂和敦化—密山断裂等构造控制的区块化特征。依兰—伊通断裂西北方的松辽盆地在上地壳深度表现为整体的低速异常，可能是由于其具有厚度较大的松散沉积层；东南方的广大山脉地区表现为大面积的高速异常。

随着深度的增加(10～30 km)，松辽盆地在中下地壳的低速异常逐渐减弱并有向高速异常转变的趋势。可能是盆地中央具有较大的沉积厚度，边缘具有较浅的结晶基底的原因。

长白山火山的低速异常逐渐显现，在15～20 km表现明显(朝鲜一侧可能也具有低速异常)并有向下延伸的趋势，可能是火山下方中部地壳发育的岩浆房或者地壳深部岩浆向浅部运移的通道，这和范兴利等[18]的研究认识相一致。

龙岗火山在上地壳深度(5 km)不具有低速异常，在中下地壳深度(20～30 km)显现出微弱的低速异常，且位于低速异常范围的边缘。从15 km、25 km和30 km的图像上可以看出，龙岗火山和长白山火山的低速异常处于同一个较大的异常范围中，但与后者相比程度微弱，这可能是由于龙岗火山岩浆喷发后的残留物导致的。

在5～20 km深度，镜泊湖火山在中上地壳深度也不具有低速异常。但从25 km、30 km的图像上可以看出，镜泊湖火山及其东北方在下地壳深度显示出一定的低速异常，原因可能是火山下方存在程度较低的部分熔融。镜泊湖火山的速度结构，与龙岗火山、长白山火山的速度结构没有直接联系，而与敦化—密山断裂有关。从岩浆类型来看，镜泊湖火山的喷出物无地壳混染，与长白山火山的喷出物类型不同。

为更好地对比和说明长白山火山和龙岗火山的S波速度结构，笔者在长白山火山附近沿纬度42°N和经度128°E分别作了2条剖面，对应图6a中的AA’和BB’线段(图7)。

黑色矩形区域为棋盘模型的范围，红色矩形区域为反演结果的范围，黑色多边形区域为射线路径覆盖的范围。a.5 km深度；b.5 km深度；c.10 km深度；d.15 km深度；e.20 km深度；f.25 km深度；g.30 km深度；h.35 km深度；i.35 km深度。图5 检测板测试结果Fig.5 Test results of check board

黑色三角代表研究区内的新生代火山，红色线段代表图7的剖面位置。a.5 km深度；b.10 km深度；c.15 km深度；d.20 km深度；e.25 km深度；f.30 km深度。图6 S波速度层析成像结果Fig.6 S-wave velocity tomography results

图7 S波速度结构剖面图Fig.7 S-wave velocity profiles

由图7可以看出，长白山下方的低速异常在地下6～23 km深度范围，分为一个主要的低速异常体和一个较小的低速异常体，具有向下延伸的趋势。在12 km深度附近含有少量的相对高速异常。龙岗火山下方7～15 km深度范围也存在低速异常，与相邻的长白山火山下方的低速异常是连通的，但没有向下延伸的趋势，说明龙岗火山下方可能不存在地壳岩浆囊。

由于广泛且快速的熔融物分离发生于深部地壳的热流条件下，深部陆壳因此是岩浆演化和存储最理想的位置[32]。结合图6c、图6d和图7，笔者认为长白山火山的地壳岩浆房主要位于火山的正下方，深度为16～21 km。这支持了Kyong--Song et al.[33]提出的地壳的很大一部分已被岩浆作用改造，与部分熔融相关的地壳岩浆体可能直接位于火山下方的说法。另外，Zhang et al.[34]观察到岩浆从基性喷发到中基性喷发的成分变化。岩浆房中可能发生进一步的岩浆分离和化学分异，导致中酸性熔融物的产生。化学分异后，来自下地壳岩浆房的密度较低、长英质的岩浆可能会上升至浅部，形成横向和纵向延伸的低速体。因此，从下地壳岩浆房运移的差异岩浆可能是目前在长白山火山观察到丰富的近地表热液活动的原因。