多尺度阵列嵌套组合反演宾川气枪源区横波速度结构

孙楠, 潘磊, 王伟涛, 叶泵, 王彬, 陈晓非,5*

1 中国科学技术大学地球与空间科学学院， 合肥 230026 2 云南省地震局， 昆明 650224 3 南方科技大学地球与空间科学系 深圳 518055 4 中国地震局地球物理研究所， 北京 100081 5 深圳市深远海油气勘探技术重点实验室， 深圳 518055

0 引言

背景噪声成像是利用背景噪声分析提取面波频散信息，进而反演地下介质速度结构的方法，是近年来地震学研究的热点，无论在工程物理勘探还是地壳上地幔大尺度速度结构研究方面都得到了快速的发展，其中面波频散曲线的提取是这些研究的关键问题.基于不同的研究尺度，频散提取方法可以分为两种：(1)浅地表面波勘探尺度的空间自相关法(Aki, 1957)、F-K法(Capon, 1969；周云腾和张致付，2020)、相移法(Park et al., 1998；伍敦仕等，2017)、倾斜叠加法(Xia et al., 2007)以及拉东变换法(McMechan and Yedlin, 1981; Luo et al., 2008)等;(2)大地球物理尺度的双台法(Knopoff, 1964; Yao et al., 2006)和基于程函方程的相速度测量方法(Lin et al., 2009)等.随着面波成像技术的发展，高阶面波在地下结构分析中的影响越来越被重视，利用多阶面波联合反演可以降低反演结果的非唯一性，使反演结果更接近真实结构(罗银河, 2008；Xia et al., 2003；何耀峰, 2005)，而传统方法对高阶面波的反应不敏感，很难从实际背景噪声中提取到清晰的高阶面波信息.最近陈晓非课题组提出了基于台阵数据的频率-贝塞尔变换法(Frequency-Bessel, F-J方法)，能够有效地从背景噪声中提取分辨率高、频率范围广的基阶和高阶面波频散曲线(Wang et al., 2019；Pan et al., 2019；Wu et al., 2020)，无论在浅地表勘探(李雪燕，2019；吴华礼，2019)，还是地壳、上地幔等大尺度的速度结构成像研究中都取得了很好的应用成果(王建楠，2019；Zhan et al., 2020; Wu et al., 2020)，对被动源和主动源数据研究都有不错的效果(Li and Chen, 2020)，是一种应用前景非常广阔的新方法.

但是，对于横向结构变化非常大的区域，目前背景噪声成像方法存在一些矛盾.即小尺度阵列对于局部浅层结构的反演成像有效，能够得到分辨率较高的浅层速度结构，但对深部速度结构不敏感.而大尺度阵列对深部速度结构研究有效，对浅层结构却约束不够.因此，为得到分辨率较高地从浅层到深层的速度结构，本文中我们采用了一种多尺度阵列嵌套组合的方式提取Rayleigh面波频散曲线.即用小尺度阵列资料提取相对高频部分的频散曲线，而用大尺度阵列资料提取较低频率的频散曲线，融合不同尺度阵列资料的分析结果，获得宽频带频散曲线，通过对宽频带频散曲线的反演获得研究目标区域从浅至深(0～70 km)的横波速度结构.

云南位于青藏高原东南缘，受板块推挤作用，区域内横向结构复杂，中强地震持续活跃，地下构造特征备受关注.为探测和监测区域地下结构变化，2012年在云南宾川建成了水库激发气枪震源发射台，利用气枪震源可重复性、频率低、传播距离远的优势进行了一系列研究(陈颙等，2007；王彬等，2016；王宝善等，2011；陈蒙，2014；孙楠和孙耀充，2020；王伟涛等，2017；向涯等，2021).为得到气枪源发射台周边区域的浅层速度结构，2017年3—5月围绕气枪发射台开展了密集台阵实验，记录到大量的连续数据，本研究利用这一密集台阵记录到的连续背景噪声数据，采用F-J方法提取Rayleigh波频散曲线.由于该地区近地表横向结构变化大、地质地貌格局复杂，直接获得较宽频带频散曲线的难度较大，因此，我们采用从密集台阵到云南地区多尺度台阵嵌套组合的方式，来拓宽频散曲线的频段，进而反演得到研究区域不同深度的横波速度结构，为该区地下结构的探测和监测提供基础.

1 台站数据

2017年密集台阵实验围绕宾川气枪震源发射台布设，布设时间长度从两周到两个月不等，记录了大量的连续波形数据，台阵覆盖面积40 km×40 km，台间距为2～3 km，共使用了381个短周期地震计，频带范围为5 s～150 Hz，覆盖区内涵盖盆地和丘陵，地表起伏最大1.4 km，下文中该区用BA来表示.本文要反演的浅层区域为BA区的一部分，位于密集台阵内部气枪发射台的西南侧，覆盖区域长25 km宽10 km，包括57个台站，下文中该区域用PA来表示.

为反演更深层的信息，我们将研究区域不断扩大，收集了宾川气枪源台网记录的连续波形数据，台网覆盖滇西地区(24.25°N—27.00°N，98.70°E—102.05°E)，共包括60个台站，台间距最大能到200 km以上，数据记录时间范围为2015—2018年，下文中该区用BC来表示.之后，我们将区域继续扩大到云南地区(20°N—30°N，96°E—107°E)，收集了云南区域台网2016—2019年3月共72台的连续波形数据，该区内台间距最大为927 km，下文中该区用YN来表示.本文所有区域的台站数据均选取了垂直分量的连续记录，采用多尺度阵列嵌套组合的方式进行分析计算，具体台站分布见图1，4个阵列具体信息见表1.

图1 多尺度阵列台站分布(a) 云南地区-YN； (b) 滇西地区-BC； (c) 2017年密集台阵区-BA，其中红色矩形框内为本文研究的最小区域-PA，台站用三角表示，黄色表示云南区域地震台网，绿色表示宾川气枪源台网，蓝色为2017年密集台阵，红色为研究小区域台站，黑色五角星为气枪发射台位置.Fig.1 Maps showing station distribution of multi-scale arrays(a) Yunnan area-YN； (b) Western Yunnan area-BC； (c) Dense array area in 2017-BA, where the red rectangle is the smallest study area in this paper-PA. Triangles represent stations, yellow for Yunnan regional seismic network, green for Binchuan airgun source seismic network, blue for dense array station in 2017, red for dense array of the smallest study area in this paper, and black pentagram indicates location of airgun source launcher.

表1 多尺度阵列信息表Tab1e 1 Information of multi-scale arrays

2 数据处理

文中数据处理主要包括4步：连续数据预处理、背景噪声互相关、提取频散信息和反演横波速度，具体流程见图2，对不同阵列数据分别作数据预处理，去除仪器响应的影响，经过多次筛选比较，YN区和BC区连续波形降采样到10 Hz，密集台阵降采样到100 Hz.背景噪声互相关计算采用姚华建提供的程序包(Yao et al., 2006, 2011)进行处理.针对连续波形记录中地震事件和气枪信号的影响，采用时间域one-bit的归一化操作进行消除(Bensen et al., 2007).

图2 数据处理流程图Fig.2 Flow chart of data processing

频散曲线提取方法采用陈晓非课题组(Wang et al., 2019)提出的F-J方法.该方法首先对噪声互相关函数的频谱C(r,ω)作贝塞尔变换，得到互相关系数的F-J频散谱：

(1)

其中，J0是零阶第一类贝塞尔函数.然后利用贝塞尔函数的同阶正交性，得到：

I(ω,k)=A·Im[g(ω,k)]，

(2)

其中，g(ω,k)是核函数，可以写为：

(3)

D(ω,k)具有频散特征，频散图上的点接近频散曲线时，D(ω,k)趋近于零，I(ω,k)会达到极值点，这些极值点就是频散点，而I(ω,k)可以利用背景噪声的互相关函数近似计算出来，这样就可以通过拾取频散谱图上的极值点来获得频散曲线信息.

提取频散曲线后，我们采用基于BFGS校正的拟牛顿法 (何耀峰，2005；田玥和陈晓非，2006；Pan et al., 2019)进行反演.首先以参考模型为中心，参考前人的研究经验，在±0.4 km·s-1扰动范围内随机生成50个初始速度模型(吴高雄，2020)，利用拟牛顿迭代进行反演，得到50个速度模型结果，然后将50个反演结果利用残差作为权重参数进行加权平均，计算得到最终的横波速度模型.

3 结果及分析

3.1 F-J频散谱

经过上述数据处理，得到4个尺度区域的F-J频散谱图(图3)，显示台阵区域越大，频散谱分布频段越低.在低频段，相同频率下的YN频谱分布明显比BC更窄，极值点的拾取也会更准确，PA中低频谱段噪声过大也暂不考虑，在高频段出现多阶信号，我们速度模型正演与频散谱对比，大致确定PA中基阶和1阶频散曲线的位置.基于以上考虑，分别拾取F-J频散谱上的极大值点，连线得到4个区域的面波频散曲线.图4显示YN区频散曲线频率范围为0.008～0.15 Hz，BC区频散曲线频段为0.046～0.33 Hz，BA区频段为0.18～0.56 Hz，PA区基阶频段为0.55～1.13 Hz，1阶频段为0.86～1.1 Hz.

图3 多尺度阵列的背景噪声数据提取到的F-J频散谱Fig.3 F-J dispersion spectra of multi-scale arrays from ambient noise data

图4 多尺度阵列的F-J频散曲线Fig.4 F-J dispersion curves of multi-scale arrays

将四个区域的频散谱按照区域大小顺序进行组合(图5a)，发现随着区域的逐渐扩大(PA→BA→BC→YN)，多尺度阵列的频散谱图能够很好的连接在一起，频散谱能量逐渐向低频端拓宽.图5b为组合后得到的频散曲线，更直观地表现出多尺度阵列的频散曲线在相同频段能够很好地吻合，表明越往深部区域介质横向均匀尺度越大，可以通过扩大区域尺度来拓宽低频信息，从而反演深层速度结构.

图5 多尺度阵列嵌套组合的F-J频散谱(a)和频散曲线(b)Fig.5 F-J dispersion spectrum (a) and dispersion curves (b) of combined multi-scale arrays

3.2 横波速度反演

背景噪声互相关结果经过F-J变换后，提取到多尺度阵列的频散曲线，下面将利用基于BFGS校正的拟牛顿迭代法反演不同深度的横波速度结构.反演中对于参考模型地选取，我们综合了中国地震局地球物理研究所吴建平研究员提供的速度模型和陈思文(2015)利用走时成像研究得到的滇西地区速度结构模型(图6)，参考模型共分42层，前两层层厚1 km，其余层厚2 km，最深为80 km.

图6 参考速度模型Fig.6 Reference velocity models

首先利用PA区频散曲线信息反演该区的浅层横波速度.图7为以参考模型为中心，在±0.4 km·s-1扰动范围内随机生成的50个初始速度模型分布.当只用基阶频散曲线反演的时候，50个模型的反演结果(图8a)显示：基阶频散曲线(红色实线)与实际值(黑点)拟合很好，而一阶频散曲线拟合不好，红色实线分布离散，图8c为其反演得到的横波速度结构，灰色代表50个模型的反演结果分布，颜色越深表示该速度值的比重越大，红色实线为50个反演结果加权平均得到的最终速度模型，结果显示横波速度的反演深度到4 km，0～2 km的横波速度变化特别大，从2.33 km·s-1增加到3.15 km·s-1(图8c)，罗睿洁等(2015)利用卫星遥感和地表调查认为该区域地层主要为泥盆系灰岩，那么横波速度0.82 km·s-1的变化梯度可能是从沉积层到基底层的转变.反演加入一阶频散曲线后(图8b)，50个模型的基阶和一阶频散曲线拟均能得到很好地拟合，反演模型结果(图8d)显示：50个速度模型在4 km以上更加收敛，也就是高阶面波的加入极大地降低反演结果的非唯一性，提高反演的稳定性，多个不同初始模型的反演结果更集中，减少了反演过程中浅层结果对初始模型的依赖，而且相对基阶，在同一频率下高阶面波的敏感深度更深，使得PA区横波速度的反演深度在8 km左右依然很收敛，表明该频段多阶联合反演的约束深度能达到8 km处.接下来的反演中，该区域浅层(<8 km)速度结构就采用PA区多阶反演结果，固定浅层反演结果不变，通过不断扩大台站分布区域来反演得到8 km以下的横波速度结构.

图7 扰动区间内50个初始速度模型的分布Fig.7 Initial velocity models within the disturbance interval

图8 分别利用PA区基阶和基阶加一阶的频散信息反演得到的拟合频散曲线和横波速度模型(a)(c) 利用基阶频散曲线得到的反演结果; (b)(d) 基阶加一阶频散曲线的反演结果；横波速度模型图(c)(d)中，两侧虚线为以参考模型为中心的初始模型扰动范围区间，灰色实线表示不同初始模型反演得到的结果分布，颜色越深表示该速度值比重越大，红色实线为最后反演得到的平均速度模型，下同.Fig.8 Fitting dispersion curves and shear wave velocity models inverted with PA region model(a) and (c) Using fundamental dispersion curves； (b) and (d) Fundamental plus first-order dispersion curves. In (c) and (d), dotted lines represent disturbance interval of initial models centered at the reference model, grey solid lines represent inversion results using different initial models, the darker the color, the greater the proportion of velocity value; and red solid lines represent the final average velocity model of inversion, the same below.

利用PA区频散信息反演得到浅层横波速度后，逐渐扩大区域尺度继续反演(图9).首先扩大到BA区，获取PA-BA组合后的频散曲线，可以发现基阶频段拓宽，低频段从0.55 Hz延伸到0.18 Hz，50个模型的频散曲线拟合较好(图9a)，横波速度的反演深度加深，最深达到12 km，深度10 km处出现明显的横波速度增加现象，12 km处稍微减小(图9d).

继续扩大区域到滇西地区，或取BC-BA-PA组合后的频散曲线进行反演，基阶低频段拓宽到0.046 Hz(图9b)，反演模型深度显著加深，最深能达到45 km左右(图9e)，前人的研究表明(陈思文，2015；白志明和王椿镛，2004)宾川地区莫霍面埋深约为45 km，也就是通过拓宽低频段信息，能反演得到该区从地表到地壳的横波速度结构.在深度10 km左右出现明显的高速-低速转变，这与张云鹏等(2020)利用天然地震研究该区上地壳速度结构的结果一致，该区域对应张云鹏等(2020)研究中x=-10，y=-5的区域，根据其计算结果该区域在10 km上下也存在明显的高低速变化现象.

最后将区域扩大到云南地区，获取YN-BC-BA-PA组合的频散曲线进行反演，显示低频最低能达到0.008 Hz(图9c)，反演深度能达70 km，70 km以上50个不同初始模型反演结果相对集中，10～45 km的反演结果比BC-BA-PA的反演结果更收敛(图9f)，表明低频信息的拓宽不仅能加深反演深度，对反演结果也具有一定的约束，降低了反演结果的非唯一性，在45 km左右存在横波速度迅速增加的现象，与该区域莫霍面的深度吻合.

图9 多尺度区域反演得到的频散曲线拟合图和横波速度模型图Fig.9 Fitting dispersion curves and shear wave velocity models inverted by multi-scale arrays

3.3 结果对比分析

将本文最终反演得到横波速度模型与参考模型进行对比，发现二者的深部速度结果基本吻合，而浅层结果存在明显差异，尤其是近地表处(图10a).将两种速度结构分别进行正演计算，与区域F-J频散谱分布进行比较，图10b中黄色虚线为参考模型的正演结果，红色虚线为反演模型的正演结果，显示反演模型的正演结果(红色虚线)与实际频散谱在所有频段都吻合.而参考模型的正演结果(黄色虚线)与实际频散曲线在低频段(<0.08 Hz，由YN区提取)基本重合，两者反映的深部的速度结构应该一致，图10a速度模型对比中证实了这一结果，二者在深部的速度结构基本一致.而0.08 Hz以上参考模型的正演结果(黄色虚线)与实际频散谱分布有明显区别，对应二者的速度结构在浅层存在明显差异.那么，对于该区域不同深度的横波速度结构，这种多尺度阵列嵌套组合的反演结果相对更准确.

图10 参考模型和反演模型的横波速度结构(a)与F-J频散谱正演模拟分布(b)Fig.10 Shear wave velocity (a) and F-J dispersion spectrum forward modeling (b) of reference model and inversion model

4 结论

本文基于多尺度阵列嵌套组合的方式，利用频率-贝塞尔变换方法，提取连续背景噪声数据中面波频散信息，将不同尺度的频散曲线融合进而反演，得到了宾川气枪发射台区域不同深度的横波速度结构.得到以下主要结论：

浅层横波速度反演结果显示，多阶面波反演结果更加收敛，反演深度更深.同样，利用多尺度阵列嵌套组合的方式(小区域PA—密集台阵BA—宾川BC—云南YN)，使低频信息不断拓宽(0.55 Hz—0.18 Hz—0.046 Hz—0.008 Hz)，对反演结果也能提供一定的约束，降低反演结果非唯一性.

从该区域反演得到0～70 km的横波速度变化特征可以看出：速度随深度曲折增长，浅层0～2 km处迅速增加，可能表明该区域沉积层的厚度在2 km以上的范围内.速度在深度10 km附近存在高低速的转变，说明该区域地壳内部存在速度异常体.在深度45 km处横波速度出现明显增大，与前人研究所得的该区域莫霍面的埋深相吻合，从结果上证实了本文利用多尺度阵列嵌套组合的方法是可行的.将反演得到的速度模型正演结果与实际频散谱进行对比，发现反演得到的区域横波速度的结果相对更准确，更能反映地下的真实介质情况，为该区地下结构的探测和监测提供了基础.

本文提出的多尺度阵列嵌套组合方式，为背景噪声成像提供了一种新的思路和方法.通过将多个阵列提取的频散曲线组合，得到频带更宽的频散曲线，既可以反演深度范围更广的速度结构，又可以对反演结果提供更好地约束，而且该方法操作简单方便，实用性极强，在区域结构研究中具有非常好的应用前景.

致谢感谢审稿专家为文章的修改提供的宝贵建议，感谢中国地震局地球物理研究所张云鹏博士、云南地震局彭关灵工程师和王光明工程师为本文提供资料.