密集台阵背景噪声双聚束成像化龙断裂精细结构

吴晓阳， 谭俊卿， 郭震， 任鹏飞， 王力伟，叶秀薇， 陈永顺，4，5

1 哈尔滨工业大学， 哈尔滨 150000 2 南方科技大学海洋科学与工程系， 深圳 518055 3 广东省地震局， 广州 510070 4 南方海洋科学与工程广东省实验室(广州)深圳分部， 广州 511458 5 上海佘山地球科学国家野外科学观测研究站， 上海 201602

0 引言

Campillo和Paul(2003)建立了空间两点间的尾波互相关函数与格林函数的联系，拉开了地震学上利用地震背景噪声进行地球内部速度结构成像研究的序幕(Shapiro and Campillo，2004；Shapiro et al.，2005；Yao et al.，2006；Lin et al.，2007；Yang et al.，2007).由于环境噪声场中的能量以面波为主，面波的频散特性使不同周期的信号对地下不同深度的介质敏感，且越长周期的信号敏感深度越深，因此背景噪声面波成像可用于多尺度研究(上地幔，地壳及浅地表等)，有效地弥补了传统地震面波成像高频能量衰减过快的不足(Yang et al.，2008；Yang，2014).近年来，随着超级计算机、背景噪声成像算法的发展以及低成本的便携式节点地震仪的广泛应用，推进了背景噪声与密集台阵联合成像技术的发展(如，Lin et al.，2013；节点地震仪数量达到5200+).在地下浅层结构探测中，密集台阵背景噪声成像方法实施简单，对场地要求低且无破坏性，可在同一地区多次观测，摆脱了人工震源不易开展和天然地震时间空间不可控的限制，是一种安全、环保、经济的高分辨率浅层探测手段，在区域地质构造研究中有着其他方法无可比拟的优势.此外，相比于传统噪声成像方法，密集台阵背景噪声成像技术极大地提高了地震波场的空间和时间分辨率，而与高效成像方法的结合则实现了高分辨率的速度结构和波场空间特性的快速提取，如程函(Eikonal)或赫姆霍兹成像(Helmholtz tomography)，波场梯度法(wavefield gradiometry)，聚焦法(focal spots)等等(Liang and Langston，2009；Lin et al.，2009；Lin and Ritzwoller，2011；de Ridder and Biondi，2015；Hillers et al.，2016).基于射线理论的传统背景噪声成像方法一般间接利用相位速度测量提取台站间频散曲线(如，时频分析)，由于固有相位存在一定的不确定性，计算时可能会出现偏差；经典的Eikonal或Helmholtz方程成像方法，通过追踪波前的相位扰动和计算走时场梯度直接获得各点相速度，无需考虑反演问题，避免了传统成像由于各种假设和中间过程所带来的误差，同时也可以提取方位角各向异性信息(Lin et al.，2009；Mordret et al.，2013).

Beamforming也称单聚束阵列研究方法，通过信号叠加提取地震波能量在波场中传播的相干部分及其传播特性，进而估计入射波的慢度和反方位角范围等信息.Krüger等(1993，1996)将Beamforming方法与核动力源结合分析地震波传播的不对称和多路径效应，发现了核-幔边界的不均匀性.近年来，基于Beamforming的思想，利用背景噪声互相关成像源和接收点互易性的特点发展的双聚束Double beamforming(DBF)成像方法，能够进行微弱体波、面波信号的重构，通过叠加两个地震阵列之间互相关信号的能量相干部分达到增强信号相关性的目的，提取原本不清晰的微弱信号，因此DBF方法充当了两个地震阵列之间的空间滤波器.Nakata等(2016)在火山口附近利用DBF对三个子阵列的噪声互相关信号叠加实现了体波信号的重构，并识别出有助于提取两个阵列之间波场特征的反射面波信号，为火山区的结构监测研究提供重要依据.Castellanos等(2020)利用位于美国长滩的5200多个短周期密集台阵记录，通过DBF方法叠加并识别出互相关信号中微弱的体波信号，进行了互相关体波成像研究.

DBF方法同样适用于背景噪声面波成像研究.Boué等(2014)利用DBF完成美国中部地区的地壳尺度背景噪声成像，并讨论了利用DBF提供的方位角信息进行射线路径重构的问题；Wang等(2019a，b)首次将DBF方法用于局部尺度和区域尺度的线性密集阵列背景噪声面波成像研究；Gkogkas等(2021)利用DBF对人口密集的沉积盆地进行浅层高精度成像，为城市地震风险评估提供了重要依据.利用DBF进行成像研究，通过对密集阵列中的两个子阵列之间进行局部慢度(方位角)搜索可直接提取局部的波场特征，能反映方位分布良好的地震阵列中方位角各向异性信息，无需考虑任何层析反演参数(如正则化参数，对各向异性模型有重要影响).此外，DBF通过对相干波形的叠加，有效地应对了环境噪声场中的非扩散噪声带来的挑战(如，交通及其他人为活动).DBF方法同样适用于多尺度的背景噪声密集阵列成像研究，但与Eikonal成像类似，对台站的空间分布要求相对较高.

1 研究实例

1.1 研究区构造背景

珠江三角洲位于欧亚板块东南缘，属华南褶皱系的一部分，经历了海西—印支、燕山和喜山等复杂的地质构造运动，主要表现为继承性断裂活动和断块差异升降.珠江三角洲主要有三组断裂系统：北西向、东西向和北东向，多属于正断层.其中，北东向和东西向断裂形成于燕山运动时期，北西向断裂形成较晚且规模较小，可能成型于喜山运动末期.三组断裂系统互相交错，将三角洲基底切割成大小不等、运动速率不一的断陷和断隆(黄玉昆等，1983).与其他区域相比，珠江三角洲的地震活跃性相对较低，目前主要以微震为主且与活动断裂关系密切.值得注意的是震级大于4.8级的地震明显受北西向活动断裂控制，形成北西向成带的震中分布，这与三角洲形态是一致的.因此地震活动与北西向活动断裂及三角洲的形成、演化密切相关(钟建强等，1996).

番禺断隆作为珠江三角洲内以丘陵为地貌特征的隆起之一，在断裂围限作用下表现出一定的地震活动性，已有地震记录显示该区曾发生过中强地震(如1824年，市桥西北发生5级地震；姚衍桃等，2008).番禺断隆东侧的化龙断裂是珠江口水域西北向断裂之一(图1)，北起黄埔西北的吉山，往南南东经沙井村，切洪圣沙东缘地带，穿珠江至南岸化龙一带后继续向南南东延伸至南沙入伶仃洋，全长约43 km.已有关于化龙断裂浅地表结构研究多来自于地质剖面和钻孔等资料，钻孔资料获取地下50 m以内的岩土层信息，断层的破裂深度暂时没有详细资料可以考察，化龙断裂在第四纪是否活动尚不明确.化龙断裂在近地表的结构特征(断层形态和展布等)对于整个区域的地震危险性评估、地壳水文研究等均具有重要意义，断层周围高度破碎的岩石损伤区能圈闭地震能量并显著放大地面运动.因此本研究开展化龙断裂地下精细结构的地震学成像研究，对理解整个珠江三角洲的构造特点和发震机制有着重要意义.

图1 珠江三角洲构造地形图及横跨化龙断裂带线性密集台阵分布图 (a) 中红色虚线框标注了研究区所在位置，黄色圆圈标注了2008年以来珠三角及周边区域的地震分布，化龙断裂—HLF； (b) 化龙断裂密集台阵分布图，黄色三角为台站位置，红色实线为化龙断裂.Fig.1 Tectonic topographic map of the Pearl River Delta and the distribution map of linear dense seismic arrays across the Hualong fault zone The red dotted box in (a) shows the location of the study area, the yellow circle marks the distribution of earthquakes occurred in the Pearl River Delta and surrounding areas since 2008. Hualong fault—HLF；(b) Distribution map of dense seismic arrays in HLF. Seismic stations and HLF are denoted by yellow triangles and red line，respectively.

1.2 数据

本研究的密集地震台阵数据来自南方科技大学布设在跨化龙断裂带的125个三分量短周期地震记录(仪器型号为EPS-2-M6Q，频带范围为0.2～150 Hz).该地震台阵呈一维线性阵列排布，整条测线总长度为12.4 km，台间距为100 m左右(图1)，整条测线所跨区域的地势西南高东北低，西南部位于低丘地台，东北部位于冲积平原.本文研究用到了连续地震背景噪声记录，记录时间为2019年6月20日至7月10日(共21天)，采样率为200 Hz.

2 DBF背景噪声成像

图2 001台站与其他台站Z-Z分量互相关波形 (周期范围0.8～2 Hz)Fig.2 Waveforms of vertical-vertical cross-correlations between station 001 and all other stations (for a period range of 0.8～2 Hz)

DBF噪声成像的首要工作是计算所有台站之间的互相关函数.互相关函数的计算按照Bensen等(2007)提出的流程，包括：(1)将原始噪声记录的垂向分量(Z分量)截成1小时一段并重采样到20 Hz；(2)去除仪器响应、去趋势和平均值；(3)对信号在目标频带内进行带通滤波处理，滤波范围为0.5～3 Hz；(4)时域归一化(running-absolute-mean)和频域谱白化处理；(5)计算每一天的互相关函数，并通过相位加权叠加(tf-PWS)得到最终的互相关.第3步的滤波处理流程对于城市内短周期背景噪声互相关信号的有效提取非常关键.如果对未经滤波处理的噪声信号直接进行互相关计算，目标频带的能量会被强烈的高频噪声(城市活动)及波场中的非扩散噪声压制而无法提取到高质量的互相关信号，已有研究(Zhou et al.，2020)也证明了该步骤的重要性.图2为位于测线东北端的001台站与其他台站在0.8～2 Hz频带范围内的Z-Z分量互相关函数.在互相关的正负半轴均能看到基阶瑞利面波信号，但也夹杂着一些散射和反射信号等.

获取互相关结果后进行局部慢度搜索，流程与Wang等(2019a)总结的基本一致.首先沿测线两端按照200 m等间距选取63个格点作为聚束中心点(测线总长12.4 km).根据阵列的台间距分布(约100 m)，避免每一对聚束中叠加的波形数量过少，保证计算结果的稳定性，同时尽可能提高成像的空间分辨率，设定聚束宽度D=250 m，D也决定了成像的水平分辨率.理论上会产生62×62共3844对源点-接收点聚束对，但是为了避免两个聚束间的混叠并尽可能满足互相关远场条件，设定聚束对间的最小距离为3 km.传统的远场准则一般规定台间距至少大于目标频带波长的1～1.5倍(Yao et al.，2006；Luo et al.，2015)，在此我们使用了并非严格的远场条件，主要依据图2中互相关的计算结果(3 km内互相关信号杂乱或观察不到目标信号).图3展示了一对源-接收聚束示意图.

图3 源-接收聚束对示意图 蓝色三角形表示台站，Xsc和Xrc表示源-接收聚束中心点位置， D为聚束宽度.Fig.3 Illustration diagram of a source-receiver beam pair Blue triangles are seismic stations. Xsc and Xrc are source-receiver beam center points，respectively. D is the beam width.

之后，根据公式(1)完成互相关信号校正叠加计算(Thorson and Claerbout，1985；Rost and Thomas，2002).为避免采样率对时间域计算精度的限制，时间校正计算在频率域进行(Wang et al.，2019b).

+τ(Xr,j,ur,θr),Xs,i,Xr,j),(1)

其中，p表示处理后的互相关波形，b表示叠加后的互相关波形，ns和nr分别为每一对源-接收聚束内的台站数量，为保证计算结果的稳定性，当聚束对内满足条件的互相关数量小于5时，该聚束对将被舍弃.τ表示互相关中的源(接收)台站相对于源(接收)聚束中心点的校正时间，具体通过公式(2)计算获得，其中u和θ分别表示慢度和聚束对之间的方位角，(xc，yc)表示聚束中心点位置.

τ(X,u,θ)=u(x-xc)sinθ+u(y-yc)cosθ.

(2)

本研究中DBF方法基于平面波入射假设(故不做方位角搜索)，通过设定慢度搜索范围并展开每对聚束对的二维网格搜索，当叠加波形的包络线具有最大能量时所对应的慢度值即为源-接收聚束的局部最优慢度值.图5展示了图4中波形经不同慢度校正叠加后的波形和包络线能量.为提高搜索的效率，每次搜索采用两步法(粗网格和细网格)进行.如图6所示，首先选择粗网格间距(du1=0.05 s·km-1；图6a)在0.20～0.80 s·km-1范围内进行搜索，确定最优慢度ux1/y1后缩小网格间距在ux1/y1±du1范围内使用细网格间距(du2=0.005 s·km-1；图6b)继续搜索，其搜索效率明显高于直接采用细网格搜索.

图4 一对源-接收聚束对中所有满足条件的互相关信号 (a) 原始互相关波形； (b) 截取并归一化后的波形.所有波形按照台间距从下向上排列.Fig.4 The cross-correlation signals satisfying the conditions in a source-receiver beam pair (a) Original cross-correlation waveforms； (b) Intercepted and normalized waveforms. All waveforms are arranged from bottom to top according to interstation spacing.

图5 图4中互相关经慢度搜索偏移叠加后的波形和包络线 在us=0.4 s·km-1及ur=0.4 s·km-1时，包络线具有最大能量.每幅子图中标注了最大振幅值.Fig.5 Waveforms and envelopes of cross-correlations after slowness search and migration superposition in Fig.4 When us= 0.4 s·km-1 and ur=0.4 s·km-1, the envelope has the maximum energy. The maximum amplitude values are labeled in each subfigure.

图6 两步法网格搜索最大包络线能量示例 横纵坐标分别表示源-接收聚束中心点的局部慢度搜索范围，色条对应包络线能量强弱，绿色叉号对应搜索到的最佳慢度值的位置.Fig.6 An example of a two-step grid search for the maximum envelope energy Horizontal and vertical coordinates represent the local slowness search range of source and receiver beams respectively, and the color bar corresponds to the strength of envelope energy. The position of the green cross marks the location with the most optimal slowness.

计算过程中，不断移动源聚束和接收聚束中心点位置，因此对同一聚束中心点来说，将局部慢度的多次测量结果取平均即获得该点最终慢度值，慢度不确定度σs由公式(3)确定(Lin et al.，2009)，n表示每一个聚束中心点的有效测量个数，s0表示该点平均慢度值.因此，仅通过每个聚束中心点的局部慢度搜索直接获得相速度图，无需反演过程.

(3)

3 结果

3.1 相速度

按照前面介绍的步骤，利用DBF方法计算了整条测线在0.6～1.2 s周期范围内的慢度值，并对每个周期的慢度剖面进行平滑处理，每个位置的平滑慢度和不确定度为相邻三个点的平均值.图7展示了0.8 s周期在平滑处理前后相速度慢度剖面，误差棒表示了慢度不确定度(标准差)，可以看到经平滑处理的整条测线的结果更加稳定.

图7 0.8 s周期慢度剖面 (a) 未经处理； (b) 平滑处理.Fig.7 Slowness profiles at 0.8 s period (a) Profile before smoothing； (b) Profile after smoothing.

在提取了各周期的慢度和慢度不确定度后，通过插值计算得到0.6～1.2 s周期频带的相速度剖面(图8).结果显示，该测线在0.6～1.2 s周期范围内的相速度随周期增大而稳定增大，相速度介于2.2～3.25 km·s-1.测线的东北侧在0.6～0.7 s周期范围内出现了相对低速异常，低速异常约为2.2～2.4 km·s-1，且低速区的位置位于图中标注的化龙断裂位置附近.而在测线的西南侧，该周期的相速度达到2.4～2.6 km·s-1；在0.75～1.05 s周期范围，整条测线相速度值介于2.65～3.1 km·s-1之间，变化相对稳定；但在1.05～1.2 s周期，整条测线的相速度逐渐增大到3.25 km·s-1.

3.2 剪切波速度

接下来，我们提取了每个聚束中心点处的频散曲线，并利用Herrmann(2013)发展的迭代阻尼加权最小二乘算法单独反演每个聚束中心点下方的一维剪切波速度结构.由于在均匀半无限空间的泊松介质中，可近似估计1/3倍波长深度的剪切波速度为相速度的1.1倍(Fang et al.，2015)，因此基于上述关系对每个聚束中心点建立一个初始速度模型，并通过20次迭代反演得到最终的Vs速度模型.迭代过程中，通过Vp/Vs=1.75的恒定关系更新Vp速度模型，并根据Brocher(2005)给出的经验公式更新密度值.图9展示了聚束中心点A和B(图8)在0～1.2 km深度范围的一维S波速度结果和拟合频散曲线.然后将每个聚束中心点的一维S波速度结果合并，获得整条测线的二维S波速度剖面，如图10b所示.

图8 二维相速度剖面Fig.8 2-D phase velocity profile

S波速度结果显示，上地壳浅部(0～1.2 km)的速度分布反映了研究区的断裂及地质构造的总体特征，整条测线具有相对较高的剪切波速度，与该区域沉积层较薄、基底埋深较浅的特点相对应.并将其分为三层：层1，深度范围0～600 m，S波速度小于2.6 km·s-1，主要反映了浅层沉积层和风化破碎层的结构特征.值得注意的是，在测线东北侧化龙断裂附近表现出较强的低速异常，S波速度小于2.2 km·s-1.另外，速度异常的倾角与前人给出的化龙断裂倾角相吻合，约为60°(麦炜伦，2018)，但是其延伸到地表的位置与化龙断裂的位置略有偏差(～200 m，图10b)；层2，深度范围600～900 m，S波速度介于2.6～3.2 km·s-1，整条测线的S波速度横向变化相对稳定，垂向表现出相对低速向高速的过渡；层3，深度范围900～1200 m，S波速度大于3.2 km·s-1，区域内高速扰动相对平稳.总体上，S波速度结果展示了研究区内化龙断裂在浅地表(<600 m)出现水平方向上的差异，而在深部趋于一致的特征.

此外，对S波速度结果进行了垂向分辨率测试，参考该反演结果给定一个类似的速度模型(图11a 所示，该模型在同一深度处的速度相同，垂向上速度稳定变化，但在113.47°E—113.49°E的浅地表100 m深度范围内给定一个-6%的低速扰动)进行正演计算，获得0.6～1.2 s周期的频散数据，然后对其进行反演获得新的速度模型(图11b).测试结果显示，新模型可以较好地恢复输入模型中浅层的低速扰动.因此，本研究中S波结果在浅层的垂向分辨率可达到100 m，反演结果是可靠的.

图9 (a，b) A点处一维S波速度结果及拟合频散曲线； (c，d) B点处一维S波速度结果及拟合频散曲线Fig.9 (a，b) 1-D S wave velocity and synthetic dispersion curves at point A； (c，d) are the same as (a，b) but at point B

图10 (a) 研究区内断层和基岩分布图(麦炜伦，2018)； (b) 2-D S波速度结构剖面 珠江口断陷—PRE，番禺地隆—PY.图中红色箭头和黑色箭头分别指示了本文研究获得的速度异常边界和根据地质填图获得 的化龙断裂位置，黑色点线绘制了速度异常在浅层的形态.Fig.10 (a) Distribution of faults and bedrock in the study area (Mai, 2018)； (b) 2-D shear velocity profile Pearl River Estuary graben—PRE, Panyu uplift—PY. The red arrow shows the boundary of velocity anomalies obtained from this study and the black arrow indicates the location of HLF obtained from geological mapping. The black dotted line plots the pattern of velocity anomalies near the surface.

图11 垂向分辨率测试结果 (a) 输入剪切波速度模型； (b) 恢复剪切波速度模型.Fig.11 Vertical resolution test (a) Input shear wave velocity model； (b) Recovered shear wave velocity model.

4 讨论及结论

本研究基于密集台阵和DBF背景噪声成像方法建立了化龙断裂附近高分辨率的地下精细结构.提取了0.6～1.2 s频带范围的瑞利波相速度结果，并反演获得上地壳浅部0～1200 m深度范围的S波速度结构.研究结果表明：S波速度剖面准确指示了研究区的地质构造特征，最显著的特征是化龙断裂附近在浅层出现速度异常，东北侧相对于西南侧表现为低速异常，与北东倾向的正断层及所形成的断陷盆地特征相一致；低速异常的西南边界位于化龙断裂西南约200 m处，可能指示了断裂两侧薄弱的风化层边界.

珠江三角洲位于东南沿海中段的构造断陷盆地内，可细分为西、北江三角洲和东江三角洲，化龙断裂作为西、北江三角洲的东缘(东西两侧分别为珠江口断陷和番禺断隆；陈伟光，1998)，切割元古代变质岩、燕山期花岗岩及白垩系碎屑岩基底(图10a；麦炜伦，2018).此外，近代沉积是叠置在早期断块格局之上，随着晚更新世以来的地壳活动性增强珠江三角洲出现了大幅度沉降，由于受多组活动断裂交叉切割且不同断裂的活动性差异较大，所以断陷盆地的形成和发展与基底断裂有着很强的关联性，表现为断裂控制了断陷内部各个不同构造单元的沉积过程，沉积岩相变化和厚度差异均受到下覆断块差异升降运动的影响(黄玉昆等，1983；张馨予等，2019).因此位于沉积盆地边缘的化龙断裂两侧的速度异常是番禺断隆和珠江口断陷的沉积层厚度、压实率和岩性成分差异的体现，同时也反映了盆地东侧的珠江口断陷相对于番禺断隆的差异性沉降运动，但也需要更多的地面真实资料(如钻孔)来进行解释验证.

速度剖面中900 m及以下深度处的S波速度达到3.4 km·s-1，且在断裂两侧未表现出明显差异. 珠江三角洲的岩浆活动表现为早第三纪的酸性或碱性喷发及晚第三纪的基性喷发和岩脉侵入，岩浆主要沿北西向和东西向断裂活动.研究区内广泛存在着基性岩脉侵入花岗岩体，且发现玄武岩沿断面贯入断层，因此推测900 m处及以下深度的相对高速异常与第三纪岩浆活动形成的基岩层有关，900 m以上的速度结构可能标志着浅层未固结的第四纪沉积物和由岩浆活动组成的第三纪地层之间的过渡.

化龙断裂作为珠江三角洲北西向断裂之一，对三角洲地区的沉积、构造演化过程起到重要作用.成像结果显示测线东北侧相对薄弱的低速区(图10b)可能是具有应力调整的弱化层，当所积累的应变能量不断增大并突然释放时就会引起潜在的地震活动.

此前对珠江三角洲及邻区的地下结构了解主要依靠华南地区大尺度面波成像、珠三角海陆联测和浅层地震勘测结合钻孔验证等工作(曹敬贺等，2014；王敏玲等，2015；董好刚等，2016).断裂活动信息主要通过物探、钻探及地质剖面来获得，而高效快捷的高分辨率背景噪声成像研究开展相对较少.本文利用近年来发展的DBF方法展开化龙断裂背景噪声成像研究，直接提取了0.6～1.2 s频带范围的瑞利波相速度而无需考虑反演问题，避免了传统成像由于各种假设和中间过程中人为干扰所带来的误差，清晰刻画了化龙断裂的浅部速度结构特征，有效地提高了我们对珠江三角洲内断层系统的认识.此外，可通过断裂带附近局部台站加密的方式进一步提高成像分辨率，从而为揭示断裂带的精确位置、深部延伸形态以及探讨断裂带附近弱化层的复杂损伤模式提供帮助，达到密集台阵高精度探测地下结构的终极目标.

DBF方法同样适用于多尺度(浅层、中下地壳及上地幔)的背景噪声二维密集(大型)阵列速度模型的构建，能快捷地提取高分辨率的速度结构并对方位角信息具有很强敏感性，通过迭代搜索不同聚束间的慢度和入射方位角，获得方位分布良好的地震阵列中局部方位角各向异性信息，帮助建立更真实的射线路径.因此，本文研究不仅推进了对珠江三角洲地区断裂系统的认识，更为快速构建整个珠江三角洲的地下三维速度结构模型奠定了方法基础和技术保障，希望为进一步研究珠江三角洲的区域构造特征、发震机制及潜在地震风险性评估提供参考.

致谢感谢南方科技大学参与化龙野外台站布设的各位老师和同学，感谢审稿专家和编辑提出的宝贵意见和建议.