数十米至数百米深度的地壳浅层剪切波速度结构可为地震工程场地评估、强地面运动模拟、活动断裂、地热资源勘查等提供重要的参考.对于具有较厚沉积层地区而言，可靠的浅层波速结构可为地震定位、地壳上地幔深部结构成像等研究提供可靠的浅层结构模型(Xia et al., 1999; 何正勤等, 2007; 徐佩芬等, 2009, 2013; Xu et al., 2012; Li et al., 2014, 2016; Ni et al., 2014; 李志伟等, 2015).因此，地壳浅层剪切波速度结构对于大型工程建设和基础研究均具有重要意义.目前地壳浅层剪切波结构主要的研究包括主动源和被动源方法.主动源方法主要为采用人工震源(如爆炸震源，气枪，重锤，震源车)激发地震面波信号，提取面波频散信息并反演剪切波速度结构.主动源方法具有信噪比高、探测精度高等优点(Xia et al., 1999; Wang et al., 2016).然而，主动源探测方法的成本较高，同时基于炸药震源的方法对环境有一定破坏性，在城市等人口密集地区开展探测工作受到一定限制.被动源方法不依赖于人工震源，而是利用包括天然地震事件、地震背景噪声等在内的被动震源产生的信号探测浅层结构.例如，利用地震台站记录到的近震事件P波波形，基于初至P波垂向与径向振幅比，以及来自沉积基底的Ps、Sp转换震相的时差和振幅比，可以获得台站下方数十米至数百米深度的浅层剪切波速度结构(Li et al., 2014; Ni et al., 2014).

地球上时刻存在着由海洋、大气、地震等各种震源激发的地震波，携带了噪声源及地球内部结构的丰富信息.随着地震学理论和方法的发展，基于地震背景噪声数据的研究在地壳波速结构成像、地球内部波速结构随时间变化、地震水平位置测定等方面取得了丰富的成果(Shapiro et al., 2005; Stehly et al., 2006; Zhan et al., 2011).通过提取地震台阵记录的地震背景噪声数据中的面波信号，可以获得面波频散并反演地壳及上地幔的剪切波速度结构.而在地壳浅层结构研究中，也发展了包括空间自相关(Spatial Autocorrelation, SPAC)、拓展空间自相关(ESPAC)、频率波数法(F-K)、H/V谱比法等基于地震背景噪声数据的微动方法.在地壳浅部数十米至数百米深度的波速结构探测中，基于台阵的微动探测方法已在地热资源、重大工程场地等研究领域获得了广泛应用，取得了良好的应用效果(ToksÖz and Lacoss, 1968; Ling, 1994; Okada, 2003; Xu et al., 2012).基于地震背景噪声的地壳浅层结构成像方法不受地震事件时空分布的限制，施工成本较低，对环境无破坏，因此已经取得了较快发展，成为研究地壳浅层波速结构的重要地震学方法.

Aki(1957)最早提出了SPAC方法的概念，并给出了平稳随机波动理论假设下一维和二维的空间自相关理论和应用实例.Aki(1965), Asten(1976), Ferrazzini等(1991), Chouet等(1998), Satoh等(2001)又对微动探测的原理和技术方法进行了研究，推动了微动方法的发展.何正勤等(2007)利用微动方法研究了0~3 km深度的剪切波速度结构，为确定地热钻井位置提供了依据，并得到了后期钻孔资料的验证.徐佩芬等(2009, 2013)基于微动方法开展了煤矿陷落柱和地铁沿线地下“孤石”等浅层波速结构相关的探测研究.

除了空间自相关方法外，还有一些基于地震背景噪声的地壳浅层结构探测方法.Capon(1970)利用最大似然F-K法成功地从地震背景噪声中提取了面波信号.Nakamura(1989)发展出新的台阵数据处理方法，即H/V谱比法研究场地的卓越频率及放大因子等参数，用于估计沉积层厚度.Claerbout(1968)提出了地震背景噪声互相关理论的雏形，证明了在水平层状且底层为均匀无限半空间模型中，自由地表记录的底部透射地震波自相关波形与其自激自收波形是等价的.Xu等(2013)在三峡库区巴东县黄土坡滑坡体隧道内用47个Texans 4.5 Hz检波器组成的测线记录了28个小时微动数据，利用波形互相关方法得到了3~8 Hz面波信号，并采用多道面波分析方法(MASW)提取相速度频散曲线，最后反演获得滑坡体200 m深度内的剪切波速度结构.徐义贤和罗银河(2015)详细总结了噪声地震学在浅层勘探中的应用现状，并对噪声地震学的发展前景进行了展望.上述研究进一步表明微动方法在地壳浅层结构探测中的良好应用效果.

五大连池火山是我国最重要的活火山之一(其它重要的活火山包括腾冲火山、长白山火山及镜泊湖火山)，其最近一次喷发发生在约300年前，喷发的火山口位于尾山西南侧的老黑山和火烧山.尾山位于五大连池火山群的北部，靠近整个五大连池火山区的中间区域.准确可靠的地壳浅表层剪切波结构对于火山地震准确定位和火山下方地壳深部三维结构成像具有重要意义，能够提高深部结构成像的可靠性.我们在五大连池火山区尾山火山锥附近区域布设了由29个地震检波器组成的密集台阵，对记录的连续地震背景噪声波形数据，使用微动方法提取Rayleigh面波相速度频散曲线，并采用面波层析成像方法反演了研究区地壳浅层剪切波速度结构.在此基础上，我们对五大连池火山区尾山火山锥附近区域的波速异常结构及其相关地质意义进行了初步讨论.

基于空间平稳随机波动理论和基阶面波是微动信号中的主要能量两个假设下，Aki(1957)对空间自相关理论进行了研究.在规则台阵的观测系统下，设圆心处的垂向分量波场为u(x, y, t), 而圆周上任一台站的垂向分量波场为u(x+ξ, y+η, t)，则空间自相关函数Φ(ξ, η, t)可表示为

(1)

上划线表示时间平均.假设波场以单一速度c在圆内传播，则空间自相关函数的方位平均为

(2)

上式为极坐标表达形式，ψ为圆形观测台阵中半径为r处台站与中心点台站的方位角，ξ=rcos(ψ)，η=rsin(ψ).空间自相关函数方位平均与波场功率谱密度关系可表示为

(3)

(4)

Φ(ω)为波场的功率谱密度，ω为角频率，J₀为第一类零阶贝塞尔函数.对上述信号应用窄带滤波，功率谱密度可表示为

(5)

P(ω₀)为在频率为ω₀时的功率谱密度，δ(ω-ω₀)为狄拉克函数，此时方位平均的空间自相关函数可表示为

(6)

自相关系数定义为

(7)

由于功率谱密度与台站位置无关，因此可将方位平均后的自相关系数写为

(8)

Chávez-García等(2005)研究表明：长时间的地震时间序列可以替代台站的空间方位平均.即对于两个台站来说，通过延长观测记录时间，可以达到噪声源趋向于空间均匀分布的效果，因此长时间的观测与台站的空间方位平均可以达到相似的效果.对于台站间长时间的微动记录，自相关系数按照如下公式进行计算：

(9)

ρ_ij(f)是自相关系数，是频率f的函数；U_i(f)与U_j(f)表示台站i与台站j的频率域微动数据，*表示共轭，ω₀为角频率，r为台站i与台站j间的距离，c(ω₀)为Rayleigh面波相速度.

Asten(2006)提出当台站在圆周上充分多或为各向同性波场时，抑或台站间记录了足够长时间的微动信号时，自相关系数谱的虚部为零，其实部为第一类零阶贝塞尔函数.在实际计算中，虚部一般不为零，表明噪声源空间分布并非绝对均匀，或者地震噪声波场不是稳定不变的平面波，这会在自相关系数谱计算中引入一定误差.一般来说，将自相关系数谱的实部近似为第一类零阶贝塞尔函数，即

(10)

常规SPAC方法通常将观测系统布设成同心圆，以减小噪声源方位分布不均匀对面波频散曲线可靠性的影响.为了能够使空间自相关方法在较大区域开展研究，Ling和Okada(1993)发展出了基于不规则台阵的拓展空间自相关方法(Extend SPAC，即ESPAC).ESPAC方法布置观测系统比较灵活，可以布设成一字型、十字型或呈网状的任意形状，并通过两两台站间的数据处理获得台站间的面波频散信息.随着地震背景噪声方法理论的发展，以及对于台阵布设无严格限制的优点，基于不规则台阵的ESPAC方法在地壳浅层结构研究得到广泛应用.Ohori等(2002)利用ESPAC和FK方法处理了T型台阵数据，并对获得的Rayleigh相速度频散曲线进行对比分析，发现基于不规则台阵ESPAC方法得到的Rayleigh面波相速度频散曲线更为可靠.Chávez-García等(2005, 2006)和EkstrÖm等(2009)研究表明在噪声来源方向随机分布的情况下，由两个台站获取的空间自相关函数为贝塞尔函数，与同心圆台站得到的结果相同，即台站对间长时间的地震时间序列可以替代Aki(1957)提出的规则台阵空间方位平均；而在新西兰开展的多次不规则台阵观测也表明反演得到的浅层剪切波速度结构与已有研究成果一致.不仅如此，Chávez-García和Rodríguez(2007)基于不规则台阵背景噪声数据的研究表明：频率域互相关(SPAC)与时间域波形互相关方法得到的Rayleigh相速度频散曲线是一致的.Tsai和Moschetti(2010)从理论上也证明了基于频率域互相关的SPAC方法与时间域波形互相关方法的一致性.Nakata等(2011)和何正勤等(2013)利用ESPAC方法从线性台阵记录的地震背景噪声中恢复了面波信号，测量了面波相速度频散曲线，发现与主动源多道面波分析方法的结果相比具有较好一致性.上述研究表明：基于不规则台阵的微动方法可以有效探测地壳浅层剪切波速度结构.

五大连池火山群由14个火山锥组成，实心正三角为各火山口位置，尾山位于该火山群最北侧，整体呈北西-南东走向，其东侧及南侧为五大连池堰塞湖(图 1a).研究区域广布农田和树林，高程在320~520 m之间，西侧和北侧海拔较高，东南侧较低，整体地势较为平坦.尾山火山锥高程520 m左右.我们将29个主频为2 Hz的地震检波器组成密集观测台阵，布设在五大连池火山区尾山火山锥附近(图 1b).由于尾山火山口附近森林密布、地形复杂，地震仪器主要布设在尾山火山口的西侧.卢川等(2015)利用相同类型的地震检波器开展了沉积盆地的微动探测并取得了良好效果，表明此类地震检波器适用于微动探测研究.台站位置及时间服务则依据仪器内置GPS获得，采样率设为200 Hz.仪器采用内置锂电池供电，每个仪器都记录了至少15 h微动信号.由于白天人文干扰较大，因此挑选出2014年7月7号晚上6点至7月8号早上9点共15 h数据进行进一步处理.

图 1

Fig. 1

图 1 (a) 五大连池火山区，上覆白色烟圈的实心正三角为14座火山分布图，图右下方为研究区的地理位置.从左至右依次为南北格拉球山、卧虎山、笔架山、老黑山、药泉山、老黑山、火烧山、东西焦德布山、尾山、小孤山、东西龙门山、莫拉布山，实线方框为此次研究区域尾山周边；(b) 研究区域台站分布图，实心倒三角为29个台站分布，虚线OA, OB, OC, DD′是图 6，图 8和图 9所示的四组剖面；(c)  研究区域三维地形分布图.图 1a中实线方框圈定位置为研究区域 Fig. 1 (a) Locations of 14 volcanos (triangles with smoke rings) around the Wudalianchi City (In north eastern China), the study region is in the solid box area. (b) Locations of 29 stations (inverted triangles) around the study region of Weishan. The 4 dash lines in (b) delineate the positions of 4 profiles of shear velocity, which are plotted in Fig. 6, Fig. 8 and Fig. 9. (c) Three-dimensional topographic map of the study region

在计算台站对空间自相关系数前，我们分别对单个台站的原始数据进行了预处理，包括去均值、去倾斜、去尖脉冲等操作，并对数据进行频谱分析(图 2).数据振幅谱在1~5 Hz频段内能量较强，带通滤波窗口选为0.1~10 Hz.将每个台站垂向分量微动数据按60 s分段，共计900段.首先对每段数据进行傅里叶变换，按公式(9)用一个台站的傅里叶谱与另一个台站的谱的共轭计算自相关系数，并按公式(10)取实部后进行叠加.然后采用滑动绝对平均方法对叠加后的自相关系数波形进行平滑处理，得到每个台站对的自相关系数曲线，并剔除了信噪比较低、未能呈现第一类零阶贝塞尔曲线形态的自相关系数曲线.从不同台站距的自相关系数波形可以看出，自相关系数曲线与第一类零阶贝塞尔曲线类似，且随着台站间距的增大，第一个零点对应的频率呈下降趋势(图 2c).

图 2

Fig. 2

图 2 (a) 4个台站原始数据波形图；(b) 4个台站原始数据频谱；(c) 4个台站对的空间相关系数波形 Fig. 2 (a) Observed noise record in 4 stations; (b) Power spectrum of the data in panel (a); (c) Cross-correlation coefficients as a function of frequency for 4 different station pairs

通过分析自相关系数波形的零点和极点，并与第一类零阶贝塞尔函数零点和极点进行拟合，得到了Rayleigh面波相速度频散曲线(图 3).即已知第一类零阶贝塞尔函数y=J₀(x)的零点及极点对应的x值，由相关系数曲线的零点及极点可得其对应频率f，由公式即可计算得到每个台站对的Rayleigh相速度频散曲线c(f).最终经过质量控制和挑选，共得到60条较高质量的Rayleigh波相速度频散曲线(图 3).

图 3

Fig. 3

图 3 (a) 所有台站对实测频散曲线，其中灰色短虚线为由(b)中平均速度模型计算得到的理论频散曲线；(b) 研究区一维剪切波平均速度模型；(c) 周期为3 Hz, 2.5 Hz, 2 Hz, 1.25 Hz面波相速度敏感核曲线 Fig. 3 (a) The measured Rayleigh wave phase velocity dispersion curves base on the microtremor survey. The gray dashed line is theoretical dispersion curve based on the model in (b).(b) The 1-D model of average shear wave velocity in the study region.(c) Normalized depth sensitivity kernels for 1.25~3 Hz Rayleigh wave phase velocity

本文采用面波层析成像方法对2~5 Hz的相速度频散数据进行了层析成像反演(Yanovskaya and Ditmar, 1990).面波射线路径分布显示频散数据对研究区有较好覆盖(图 4).反演采用0.1 km×0.1 km网格，并经过多次成像测试确定合适的正则化参数和阻尼系数，从而保证成像结果相对平滑且能够较好拟合观测数据(图 5).为了更清晰地分析层析成像结果随空间和周期的变化，以尾山火山口为原点O，分别沿西向OA，西南向OB，南向OC以及尾山南侧东西走向DD′(见图 1b)获得4组相速度频散曲线变化剖面图(图 6).在研究区西侧、西南侧及南侧靠近尾山火山口的区域，在2.8~4.0 Hz频段内相速度相对较高.

图 4

Fig. 4

图 4 (a) 射线路径分布图；(b) 网格射线覆盖数统计 Fig. 4 (a)Ray path coverage; (b)Ray path number in each grid

图 5

Fig. 5

图 5 多个周期的相速度分布图.右上角为周期，红色实心正三角为尾山火山锥位置.图中给出研究区周期为4.2 Hz, 3.4 Hz, 2.8 Hz, 2.2 Hz的面波相速度水平层析成像结果 Fig. 5 Tomographic images for Rayleigh wave phase velocity at 4.2, 3.4, 2.8 and 2.2 Hz. The period for each subfigure is labelled on the upper-right corner. Red triangle shows the location of Weishan volcano cone

图 6

Fig. 6

图 6 以尾山火山口为原点O，沿西向OA，西南向OB，南向OC以及尾山南侧东西走向DD′(见图 1b)等4组测线的相速度频散曲线变化剖面图，四组剖面的色标相同 Fig. 6 (a—d)Cross sections of phase velocity dispersion curves along OA, OB, OC and DD′ profiles delineated in Fig. 1b. The solid line on the top is elevation along OA, OB, OC and DD′ profiles delineated in Fig. 1b

面波相速度频散数据不同周期的敏感核测试表明： 2.0~3.0 Hz频散数据对0~700 m结构均有一定分辨能力，对100~400 m结构较为敏感，其中2 Hz频散对100~500 m较浅深度具有较好分辨能力(图 3).为了获得研究区三维剪切波速度结构，将每个网格的相速度频散数据采用线性方法进行反演(Herrmann and Ammon, 2004)，得到1-D随深度变化的剪切波速度结构，并通过插值平滑便可获得研究区0~700 m的三维剪切波速度结构(图 7和8).100 m深度的剪切波速度图像显示：剪切波波速变化范围为600~950 m·s^-1，在靠近尾山火山锥和研究区西北部为高速异常，研究区东南部则为低速异常；至200 m以下深度，靠近尾山火山锥区域和研究区西北部则显示为低速异常，研究区东南部为高速异常.为了更好反映研究区波速结构变化特征，图 8和图 9分别给出了绝对和相对剪切波速度变化剖面.图 8给出了由尾山火山口向西1.7 km(OA)，向东南2.0 km(OB)，向南1 km(OC)及尾山南侧东西方向长度为2 km的测线DD′等四组剖面，清晰地展示了0~700 m深度绝对剪切波速度变化特征.OA，OB，OC三组由火山口为原点的测线海拔由480 m逐渐降为360 m，尾山近处100 m以上浅层剪切波速度较远处高，并呈现出地表低速覆盖层越靠近尾山越薄的特征.DD′测线由东向西海拔由400 m降至360 m，整体较为平坦，同样地，100 m深度以上的浅层，靠近尾山火山锥处的剪切波速度要较远处高，并呈现出低速覆盖层越靠近尾山越薄的特征.图 9同样反映了OA，OB，OC，DD′四组剖面在150 m深度以上的浅层地表尾山火山锥近处速度较远处偏高，而在150~700 m深度范围内，尾山火山锥近处速度较远处偏低.

图 7

Fig. 7

图 7 剪切波速度不同深度剖面图.右上角为深度，红色实心三角为尾山火山锥，色标为剪切波速度值变化 Fig. 7 Tomographic images for shear-wave velocities at 0.1~0.5 km depths. The depth for each subfigure is labelled on the upper-right corner. Red triangle shows the location of Weishan volcano cone

图 8

Fig. 8

图 8 (a)  OA测线剪切波速度剖面；(b) OB测线剪切波速度剖面；(c) OC测线剪切波速度剖面；(d) DD′测线剪切波速度剖面；图形上部为该测线地形变化，四组剖面色标相同，均为剪切波速度变化.四组剖面方位分布如图 1b所示 Fig. 8 (a—d) Cross sections of shear wave velocity along OA, OB, OC and DD′ profiles delineated in Fig. 1 (b). The solid line on the top is elevation along OA, OB, OC and DD′ profiles delineated in Fig. 1b

图 9

Fig. 9

图 9 (a) OA测线区域平均剪切波速度扰动剖面；(b) OB测线区域平均剪切波速度扰动剖面；(c) OC测线区域平均剪切波速度扰动剖面；(d) DD′测线区域平均剪切波速度扰动剖面；图形上部为该测线地形变化，四组剖面色标相同，均为区域平均剪切波速度±20%速度扰动.四组剖面方位分布如图 1b所示 Fig. 9 (a—d) Cross sections of the average shear wave velocity perturbations along OA, OB, OC and DD′ profiles delineated in Fig. 1b. The solid line on the top is elevation along OA, OB, OC and DD′ profiles delineated in Fig. 1b

鉴于面波频散数据覆盖范围的限制，本文只展示和讨论了研究区主要的异常结构特征，而更小尺度异常特征可能具有较大不确定性不进行讨论.野外地质调查发现在尾山火山锥附近区域地表覆盖有松散沉积土壤层，且靠近火山锥附近广泛出露有多孔隙玄武岩.因此，推测尾山火山锥附近的剪切波速度异常特征可能由以下因素造成：尾山为整体走向为北西-南东，尾山火山口近处及研究区西北部地势较高，有黑色多孔隙玄武岩出露，东南部地势较低，地表分布为松散沉积层，推测靠近尾山火山锥处及西北部深部为火山喷发的玄武岩，在波速结构成像结果中显示为高速异常；而在远离尾山火山锥处，其浅部的低速异常反映了地表松散沉积层的影响；并且越靠近尾山火山锥，松散沉积层的厚度越薄(图 7—9).至200 m以下深度，剪切波速度分布特征与火山喷发时火山口附近的岩浆分布密切相关：在靠近尾山火山锥区域，火山作用造成火山通道周围较为破碎多裂隙发育，加之火山岩为多孔隙玄武岩富含孔隙流体，使之显示为低速异常；而在远离尾山火山锥区域，其异常特征可能反映了具有较高剪切波速度的结晶变质岩的影响.

本文利用由29台检波器组成的无线地震检波器密集台阵，在五大连池火山区尾山火山锥附近开展了连续地震背景噪声观测，并基于微动和面波层析成像方法获得了研究区地表至700 m深度的地壳浅层三维剪切波速度结构.五大连池火山区尾山火山锥附近的三维剪切波速度成像结果显示：在0~150 m的浅部靠近尾山火山锥的区域显示为相对高速异常，远离火山锥的区域显示为相对低速异常.至150~700 m深部，波速异常特征则明显与0~150 m深度的波速特征相反，即靠近尾山火山锥的区域显示为相对低速异常，远离火山锥的区域显示为相对高速异常.在远离尾山火山锥区域，浅层的相对低速异常可能与松散沉积层有关，深部的高速异常则反映了结晶变质岩的影响；在靠近尾山火山锥区域，浅部的相对高速异常应该反映了出露地表的玄武岩，而深部的相对低速异常则可能反映了火山通道周围广泛发育的破碎裂隙结构及其火山喷发后孔隙流体填充的影响.

感谢审稿专家宝贵的修改意见和建议.感谢倪四道教授的指导和帮助.