微动是指地球表面非地震引起的微弱振动(Okada和Suto, 2003)，包含体波和面波信号等成分，其中面波的能量占微动信号总能量的70%以上(Toksöz和Lacoss, 1968).研究表明，可以从微动信号中提取得到面波相速度或群速度频散信息，并使用它们对地下结构进行成像(Aki, 1957; Campillo和Paul, 2003).目前，从微动信号中提取面波频散信息的主要方法有：空间自相关法(Aki, 1957)、频率-波数法(Capon et al., 1969; Lacoss et al., 1969)、多道面波分析法(Park et al., 1998; 2007; 2008;Xia et al., 1997; 2000a)和τ-p变换法(McMechan和Yedlin, 1981)等.

利用微动信号提取面波频散信息反演地下S波速度结构的方法具有经济快捷、受场地限制小和环境无损等优点，被广泛应用于大尺度的地下结构探测和小尺度的工程勘探中.在大尺度的地下结构探测中，Shapiro等(2005)利用USArray台阵一个月的背景噪声信号进行了互相关计算并提取面波群速度信息，对美国加利福尼亚地区的地下结构进行了研究，成像结果显示速度差异明显，其中低速区域主要为沉积盆地，高速区域主要为山脉的火成岩地区；Yao等(2006)利用背景噪声信号计算经验格林函数，对青藏高原东南部的地下结构进行了成像，发现青藏高原东南部地壳和上地幔剪切波具有低速特点；Moschetti等(2007)使用USArray台阵28个月的背景噪声数据计算并提取出了周期8~40 s的瑞利波群速度频散曲线，并反映出美国西部地区地下速度异常，并且与已知地质信息相吻合；Zhao等(2016)使用密集的USArray台阵证明了利用背景噪声成像宽频带5~150 s瑞利波相速度的可行性.在小尺度的工程勘探应用方面，Liu等(2000)在加利福尼亚的两个地区分别利用微震台阵得到了百米深度范围的S波速度结构，并分别对比地区钻孔数据信息，结果显示微震台阵技术成像的速度结构与钻孔信息相吻合，表明该方法具有实际应用价值；何正勤等(2007)在北京地区多个场地进行微动台阵探测方法试验，探查了0~3 km深度的速度结构，验证了地脉动数据划分物性差异大的地层速度界面的有效性；徐佩芬等(2012)利用SPAC方法从微动数据中提取基阶频散曲线成像二维视速度剖面，并结合少量钻孔资料成功圈定出地下“孤石”，为工程提供了基岩面埋深和起伏信息.

目前大部分研究仅利用基阶频散信息进行结构反演，没有使用高阶频散信息.这是因为高阶振型面波能量小、振幅随周期变化大，传统方法不易提取出有效的高阶面波.同时，研究者发现各模态面波在不同频率范围内能量分布不同，当某一模态面波能量占主导时，其他模态面波则相对较弱(张碧星等, 1998; 王建楠, 2017)，并且当存在低速层和速度突变时，频散曲线可能会出现“之”字形形态(鲁来玉等, 2006; 刘雪枫等, 2009)和频散曲线“交叉”现象(De Nil, 2005; Levshin和Panza, 2006; Boaga et al., 2013; 张凯等, 2016)，使得高阶频散曲线提取更加困难，且无法准确判断频散曲线模态.

随着研究的深入，研究者们发现高阶模态频散曲线具有更大的穿透深度(Xia et al., 2000b; 2003;Luo et al., 2007)，且对地层参数变化更敏感，若能够提取出准确有效的高阶模态频散曲线，并联合基阶进行反演，可以极大地降低反演的多解性，从而获得更加准确的地下速度结构(Xia et al., 1999; 罗银河等, 2008; Pan et al., 2018; 杨振涛等, 2019).因此，提取准确有效的高阶模态频散曲线是微动信号成像的焦点和难点.近期，Wang et al.(2019)提出了提取面波高阶频散曲线的新方法，F-J方法，并从理论公式推导和合成背景噪声数据方面验证了该方法的有效性.

本文首先简单介绍了F-J方法的理论和处理流程，然后使用该方法从上海市苏州河地区的微动信号中提取了高分辨率多模态瑞利波频散曲线，最后利用粒子群算法进行反演，获得了二维S波速度结构剖面.

Wang等(2019)提出一种从微动信号中提取面波相速度频散曲线的新方法，频率-贝塞尔变换方法(F-J方法)，并给出了F-J方法的详细理论推导过程和数值验证，本文仅对该方法做简单的介绍.

空间上位于x₁、x₂的任意两个观测台站记录的微动信号垂直分量的互相关正比于格林函数虚部(如Sánchez-Sesma和Campillo, 2006)，将格林函数用水平层状半空间弹性介质各向同性源的格林函数计算公式表示(Chen, 1993, 1999; Hisada, 1994; Kennett, 1986)，得：

(1)

其中，，F表示傅里叶变换，A为常数，g_z(z=0, k, ω)为与r无关的核函数，J₀(kr)为第一类零阶贝塞尔函数，k为波数，r为台站距离.根据贝塞尔函数的正交性以及各向同性条件下的频率矢量波数变换方法(Wang et al., 2019)，在方程(1)基础上应用贝塞尔变换，得：

(2)

公式(2)左边可以通过计算直接得到，根据核函数的性质(如Kennett, 1986)，当核函数趋于无穷大值时，该无穷大值就是频散曲线所对应的点，所以计算左部分提取最大值，便可以得到频散曲线.理论公式(2)左边中无限积分在实际计算中将转换为有限项的累加，将用线性方程近似：

(3)

其中， ，r_j为第j个互相谱的台站对的距离.

根据贝塞尔函数的递推性，

(4)

最终可以得到：

(5)

利用公式(5)对互相关谱在台站间距上进行线性插值，可以得到频率-相速度谱.

为了验证F-J方法可以从城市微动信号中提取出对城市浅地表勘探有效的高阶频散信息，本研究对上海市苏州河地区开展了实地微动信号采集并进行处理.

研究区域苏州河地区位于上海市，长江三角洲东南前缘，地形平坦，地下结构以层状沉积为主，除少数基岩露头外，覆盖了巨厚的第四系松散沉积物，厚度可达200~300 m.苏州河地区属滨海平原，在深度范围3~20 m普遍分布淤泥质黏性土，呈流塑性，是主要的软弱层.随深度的增加依次分布，软塑性灰色黏性土、硬土层、密实粉性土粉砂互层和致密粉细砂层，受古河道影响局部范围存在灰色黏性土厚度增大、土性差异增大和硬土层缺失的情况.深部致密粉细砂层局部夹中粗砂，分布较为稳定.潜水层水位埋深一般为0.3~1.5 m，渗透性差，当场地大面积填土后水位上升.承压水不连续分布，水位低于潜水位，埋深约3~11 m.本次测区西南方向约200 m钻孔地质信息见表 1，测区南侧紧邻苏州河河道，北侧为城市公园，四周分布公路、居民区和工地等，微动信号震源复杂，且无明显方向性.

表 1 (Table 1)

表 1 钻孔地质信息 Table 1 Geological information of borehole 序号 岩性信息 地层深度 /m 密度 /(g·cm-3) P波速度 /(m·s-1) S波速度 /(m·s-1) 1 淤泥 0.8 1.93 795 100 2 黏质粉土 4.0 1.77 1000 110 3 淤泥质 粉质黏土 5.5 1.77 1078 113 4 淤泥质 黏土 15.1 1.81 1177 133 5 黏土 25.0 1.85 1113 184 6 粉质黏土 夹粉砂 37.6 1.91 1302 230 7 粉质黏土 39.3 1.92 1655 272 8 粉砂夹 粉质黏土 43.5 1.96 1355 278 9 粉质黏土 49.0 1.86 1281 236 10 粉质黏土 64.0 1.99 1525 278 11 粉砂夹 粉质黏土 70.0 1.97 1586 380 12 粉质黏土 夹粉砂 73.7 1.97 1586 380 13 粉细砂夹 中粗砂 80.0 2.00 1620 445

表 1 钻孔地质信息 Table 1 Geological information of borehole

F-J方法对观测系统无特定排列要求，随机分布、线性排列等均适用.本次微动信号采集于2018年3月25日上午10—12点进行，沿绿化带边缘采用近线性排列观测系统，如图 1a所示，使用38台宽频带数字地震仪，测线总展布约100 m，道间距约2 m，采集频率200 Hz，静置采集1 h.

图 1

Fig. 1

图 1 测区台站分布和微动记录图 (a)测区台站分布图; (b)测区8个台站预处理后的微动波形记录. Fig. 1 Nearly-linear observation array (a) and preprocessed microtremor recordings (b)

对微动信号进行预处理，包括去除均值、去除线性趋势和带通滤波等，具体流程如图 2.图 1b给出了预处理后的8个台站的10 min微动波形记录.

图 2

Fig. 2

图 2 F-J方法流程图 Fig. 2 Flowchart of the F-J method

应用F-J方法从微动信号中提取频散曲线可分为三个步骤：(1)将记录按照一定窗长分成若干段，任意台站两两组合计算空间互相关系数得到互相关谱.窗长的选择需要保证微动信号在有效的频带范围内频谱基本一致，且能量稳定，满足时间、空间平稳随机的条件.本文将记录按照窗长30 s分成120段.(2)将所有的互相关谱进行叠加平均，随着叠加次数的增加可以快速提高后续积分成像的分辨率，但当叠加次数增加到一定程度后，分辨率的提高趋于平缓.(3)叠加后的互相关谱进行积分计算，得到频率-相速度谱，进而提取频散曲线.

F-J方法基于平层结构假设，但实际地层中无法严格满足该假设，因此计算结果包含地下结构的平均效应.台阵展布越小，越易满足完全水平层假设条件，同时展布的缩小会降低低频面波的成像效果和探测深度.本实验采用21个邻近台站为一组，分为18组(Arr.01—Arr.18)，组1(Arr.01)包括台站ST.01—ST.21(见图 1)，组2(Arr.02)包括台站ST.02—ST.22，依次类推.使用F-J方法对18组的城市微动数据进行处理，计算得到频率范围0~30 Hz、相速度范围50~450 m·s^-1的频率-相速度谱.

图 3给出了SPAC方法和F-J方法计算的频率-相速度谱，通过对比可以发现：在低频区域(< 13 Hz)两种方法基阶频散曲线成像分辨率相差不大，但F-J方法在小于5 Hz部分频散曲线的宽度展布略大于SPAC法.频率13~25 Hz区域F-J方法具有更高的分辨率，SPAC法仅能成像一条高阶频散曲线且宽度展布模糊.在高频区域(25 Hz左右)，SPAC方法有能量集中点，但最终未能形成连续频散曲线，F-J方法以高分辨率成像多条高模态频散曲线.图 3b中的点线是使用广义反射-透射法系数法(Chen, 1993; 何耀峰等, 2006)对该地区钻孔数据(表 1)建立的速度模型计算得到的基阶和1~7阶理论频散曲线.

图 3

Fig. 3

图 3 频率-相速度谱 (a) SPAC法; (b) F-J方法(点线为钻孔数据理论频散曲线). Fig. 3 Frequency-phase velocity spectrograms (a) SPAC method; (b) F-J method (Dotted lines are theoretical dispersion curves by borehole information).

图 4分别为组Arr.01—Arr.18的频率-相速度谱.低频部分(< 13 Hz)，基阶频散曲线组间形态相似，与钻孔数据理论频散曲线相比无太大偏差.高频部分(>13 Hz)，高阶模态频散曲线组间形态相差较大，在频率范围15~30 Hz，Arr.01—Arr.18组高阶频散曲线的模式融合逐渐分离.主要原因是，低频基阶模态与深部结构有关，深部地层一般分布较为稳定，高频高阶模态部分与浅层结构有关，城市浅层受到建筑垃圾、地下空间开发等人为影响结构错综复杂，因而频散曲线形态不同.基阶频散曲线的相似性也决定了仅使用基阶频散信息进行反演，会增强反演结果的多解性，加入多模态频散信息可以有效地约束反演结果.

图 4

Fig. 4

图 4 Arr.01-Arr.18频率-相速度谱 点线为反演S波速度结构的理论频散曲线. Fig. 4 Frequecy-phase velocity spectrograms of Arr.01-Arr.18 Dotted line are dispersion curves calculated from inverted S-wave velocity.

本节使用钻孔数据建立初始模型，联合基阶、高阶多模态的频散信息，利用粒子群算法对S波速度和层厚进行反演.在反演过程中，设定密度和P波速度已知.反演得到的S波速度结构的理论频散曲线(图 4中点线)和F-J方法计算的频散谱可以很好地吻合.将不同组得到的反演S波速度进行了Kriging插值，获得二维S波速度的剖面图(图 5).

图 5

Fig. 5

图 5 S波速度剖面图 (a)钻孔S波速度模型; (b) F-J法S波速度剖面. Fig. 5 S-wave velocity profile (a) S-wave velocity model from borehole information (scale is on right); (b) S-wave velocity profile by the F-J method.

结合钻孔数据岩性信息，可以将图 5中的速度剖面分为三层：(1)淤泥质黏土层：深度范围0~23 m，S波速度范围110~160 m·s^-1，为表层软弱层，岩性较均匀.(2)粉质黏土层：深度范围23~55 m，S波速度范围220~290 m·s^-1.该深度范围存在低速区(图 5中白色虚线框)，本研究综合地质信息解释为古河道切割导致的地层缺失.(3)粉砂黏土层：深度范围大于55 m，S波速度范围300~500 m·s^-1，随着深度的增加S波速度平稳增加，土层横向分布稳定且土质致密，可作为承压层.

本文使用F-J方法成功地从城市微动信号中提取出了高分辨率多模态瑞利波频散曲线，利用粒子群算法联合多模态频散曲线对浅地表的S波速度结构进行了成像，直观地显示了城市浅地表S波横向和纵向速度结构变化，揭示了局部地层中因古河道的切割造成的结构缺失.研究结果表明：F-J方法能够准确地提取多模态频散曲线，为反演计算增加约束信息，提高了反演过程的稳定性和反演结果的准确性.F-J方法对观测系统的布设要求低，而且需要台站的数量小，抗干扰能力强，适用于交通繁忙、建筑物密集、各种场源干扰严重的城区，在城市浅地表精细结构成像方面具有较好的应用前景.