珠江三角洲地区经济高度发达，人口密集，开展地质调查工作非常必要，地质调查工作对于大湾区地下空间资源合理开发、防灾减灾和基础地质研究有着重要的现实意义。在城市发展过程中，已经进行了大量的浅部地质调查工作，近几年也大面积铺开了城市地质调查工作，但这些地质调查工作主要以100 m以浅为主，中深部的地质调查工作开展的非常少，因此，开展中深部地壳结构的调查研究工作显得非常有必要。

谭皓原^[1]利用国家地震科学数据共享中心的地震目录及临时台网资料对华南地壳及上地幔三维速度结构成像；尹周勋^[2]、赵明辉^[3]使用深地震方法研究华南地区地壳结构及壳内低速层；韩松^[4]、刘国兴^[5]、韩凯^[6]使用超长周期大地电磁测深，对华南地区深部电性结构进行研究；曾维军^[7]利用重磁资料对广州—巴拉望地学断面进行了综合研究。普遍认为：地壳分为上地壳、下地壳或者上地壳、中地壳和下地壳；珠江三角洲地区莫霍面深度了30 km左右，由陆地向海洋逐渐变浅；在15～20 km深的深度上存在异常低速层，由陆地向海洋逐渐尖灭。但以往的深地震、深电性剖面为了避开城区的人文干扰，测线大多数位于珠江三角洲边缘或外围，在核心地区，缺少专门研究中深部地壳结构的资料，在地壳结构复杂的情况下，是否还存在低速带层；以往的深地震剖面震源使用的是炸药震源，对于几十公里的探测深度，能量有限。

密集台阵勘探是基于地震事件的地震体波、面波等地震成像方法，是研究地壳及地球更深部结构的主要方法之一，由于该方法成本较低，并且对中深部地壳具有高分辨率，抗干扰能力，因而该方法被广泛应用。李敏娟^[8]利用密集台阵研究了青藏高原东北缘地壳精细结构及九寨沟地震震区结构特征。黄焱羚^[9]利用密集台阵研究龙门山断裂带南段地震空段的地震活动性，为地震空段及周围地震活动的时空分布特征和孕震风险性评估提供了丰富的信息。张云鹏^[10]利用密集台阵研究云南宾川上地壳速度结构，深入理解该区的主要控震构造及其深部结构特征提供了重要依据。

为了研究珠江三角洲核心地区中深部地壳结果速度分布特征和结构特征，笔者采用密集台阵勘探方法采集数据，使用接收函数方法处理反演得到沿线地壳速度结构，结合地质资料进行综合解释，为深入研究珠江三角洲地区地壳结构提供了丰富的信息。

本次研究测线起于广东省佛山市三水区芦苞镇，经狮山镇、中山市东凤镇，至中山市三乡镇，测线由西北向东南方向延伸，测线全长124 km，共64个测点，平均点距约2 km。数据采集时段为2019年10月7日至2019年11月16日，累计时长为38 d，保证了有足够时间记录到足够多天然地震事件。数据采集过程中使用仪器为QS-5B三分量宽频（0.2～200 Hz）数字地震仪，仪器之间使用内置高精度时钟同步，数据采样间隔为4 ms。

数据采集期间记录到的天然地震事件方位角、震中距分布都比较均匀，观测期间共记录到震级大于5级、震中距100度以内的天然地震事件168余个，见图3。图4为台阵记录到的一次7.1级地震事件的波形图，从左到右分别为垂直、南北、东西分量的地震波形。从三分量波形可以看出，包括初至P波、多次波的震相波形的相位和振幅具有良好的一致性，尤其是垂直分量波形，震相波形具有高度的一致性。对于南北和东西两个分量的波形，整体上也高度一致。对于P波初至的时间，在整个120多千米长度的测线上，不同台站的P波初至与理论的P波初至相比有一定偏离，反映出测线下方的地壳结构存在显著的横向不均匀性和差异。原始数据的信号内容丰富、震相清晰，转换波信号清晰，数据质量优良。

图3 远震天然地震事件震中分布图（2019-10-07—2019-11-16） Fig.3 Epicenter distribution map of teleseismic natural earthquake events

图4 台阵记录到的一次7.1级地震事件波形图Fig.4 Waveform diagram of a M7.1 earthquake event recorded by seismic array

本次研究使用较为成熟的接收函数方法对数据进行处理，通过从远震体波波形数据中提取的P接收函数和S波接收函数可以获取地壳内部速度结构、地壳厚度及物质成分组成、地幔过渡带的厚度变化及岩石圈地幔的间断面^[15]。其原理如下：接收函数是利用P波波形的台阵记录来反演沿线下方二维S波速度结构的波形反演方法，远震P波波形含有关于震源时间函数、震区介质结构、上地幔传播路径以及接收区介质结构的丰富信息。远震P波波形与这些影响机制的关系可表示成：D（t）=S （t）*MS（t）*MRay（t）*MR（t）*I（t）。其中，D（t）为所记录的远震P波波形数据；S（t）为震源时间函数；MS（t）为近源介质结构响应；MR（t）为台站下方接收介质的响应；MRay（t）为P波在地幔中传播的透射响应；I（t）为仪器响应。在以上因素中，只有台站下方介质的响应才是我们感兴趣的可用来反演台站下方地壳、上地幔速度结构的波形信息。因此要有一种方法将接收介质的响应从整个P波波形中分离出来，而接收函数方法就是一种行之有效的方法。 Langston^[16]基于等效震源假定，提出了从远震P波波形数据中分离出径向和切向接收函数方法。

如果对各种响应做一个理想的假设，三种传播响应共同施加给震源时间函数，构成一个串联的线性系统，此外，要考虑仪器响应I（t），最终的系统输出产生位移即为u（t）。

这种串联的表示对应一个卷积模型，由于线性系统具有叠加效应，所以可以进行卷积运算。在时间域内，把这个系统表示成数学公式如下：

其中，U_V（t）、U_R（t）和U_T（t）分别表示地震记录的垂向分量、径向分量和切向分量。

MR_R（t）和MR_T（t）即为要求的径向和切向接收函数，使用时间域反褶积提取接收函数，时间域反褶积方法提取的接收函数具有较小的背景噪声且对地震波形的记录长度不限制。数据处理流程见图5。

图5 数据处理流程图Fig.5 Data processing flow chart

（1）远震事件波形截取：根据远震事件目录给出的地震发生的时间（UTC时间，世界协调时）、震中经度纬度和深度等信息，计算出远震P波初至的理论到时，从密集台阵测线数据的原始三分量波形中截取出远震事件波形，截取窗口为P波之前50 s，P波之后150 s，保证了远震波形中的地壳内部转换波和多次波信号均位于截取窗口内部。

（2）利用时间域反褶积方法^[17]，对密集台阵测线每一台地震仪记录的三分量远震事件波形进行数据处理，获得远震P波接收函数波形。

（3）对获得的接收函数波形进行质量控制，挑选出具有清晰莫霍面转换波震相的径向接收函数波形，并用于整个地壳结构的反演成像。

（4）基于远震事件目录给出的震中位置和深度，计算理论的射线参数，并根据射线参数，对密集台阵测线上同一台站记录到的接收函数波形进行叠加和动校正，获得每个台站的平均接收函数波形，并用于密集台阵测线下方地壳结构的反演与成像。

（5）速度结构反演。通过研究区域的速度结构初始模型构建、接收函数波形反演与地壳S波速度结构建立、接收函数理论地震图计算及与实际数据的对比分析采用阻尼最小二乘线性+迭代接收函数反演方法，经多次迭代更新，直至理论与实际接收函数波形差最小，最终获得能够和观测接收函数波形数据拟合最好的地壳S波速度结构参数模型。图6为观测原始数据波形（黑色）和反演后理论波形（红色）对比。理论波形能够较好拟合实际观测的数据波形，表明反演获得的测线下方的S波速度结构模型能够较好拟合观测数据。

图6 观测原始数据波形（黑色）和反演后理论波形（红色）对比Fig.6 Comparison between the observed original data waveform（black）and theoretical waveform after inversion（red）

对密集台阵剖面的反演成果进行了综合研究，讨论了测区地壳结构速度特征，结果表明：

（1）使用密集台阵勘探方法，利用5级及以上的远震天然地震作为震源，在干扰严重的城区也取得了良好的勘探效果，其具有抗干扰能力、信噪比高、分辨率高、施工便捷、成本低的优点，非常适合在干扰严重的城区勘探，可作为研究珠江三角洲地区地壳结构，特别是中深部地壳结构的勘探方法。

（2）根据密集台阵剖面速度结构特征划分了上地壳、中地壳、下地壳和上地幔四个大层位。上地壳厚度分布表明在珠江三角洲临海一侧，在地质构造的作用下，造成地壳下沉，接受沉积，导致上地壳厚度变厚，在中山三乡镇一带，上地壳厚度达14 km。

（3）根据密集台阵剖面速度结构特征获得了沿线中深部断裂构造特征及其空间展布形态，速度结构与断裂构造具有良好的一致性。

（4）本次密集台阵勘探推测三水盆地盆底深度在3 km，与该地区沉积层厚度基本一致。

（5）莫霍面波速特征明显，界面上下横波速度变化大，由3.6 km/s突变为4.1 km/s。

（6）本次密集台阵勘探取得的新成果：①传统认为重力高是莫霍面隆起所致，而本次密集台阵勘探成果表明，重力高是由莫霍面隆起、上地幔物质侵入下地壳和高密度岩浆侵入上地壳共同引起的；②沿测线地壳在15～20 km的深度上不存在低速层。

致谢：感谢徐行教授为本文提出宝贵的修改意见；感谢在密集台阵勘探数据采集、处理和总结工作中付出辛勤劳动的领导和同事；感谢李志伟老师提供数据处理服务工作；感谢韩思旭同志对本文的指导。