地球物理学报  2021, Vol. 64 Issue (5): 1699-1709   PDF    
横向非均匀介质多道瑞雷波频散曲线正演
胡书凡1, 赵永辉1, 吴健生2,1, 葛双成3     
1. 同济大学海洋与地球科学学院, 上海 200092;
2. 九三学社上海市委员会, 上海 200041;
3. 浙江水利水电学院, 杭州 310018
摘要:作为近地表横波速度结构成像的主要手段之一,面波多道分析法的正问题研究对现场观测系统设计及后续反演计算具有重要意义.目前面波频散曲线的正演主要分为两类:一是对水平层状介质中面波的本征值问题进行求解,该类方法计算效率高但较难考虑地下介质在横向上的不均匀性;二是基于波动方程的全波场模拟,该类方法在理论上可考虑任意复杂的地质模型但计算成本相对较高.本文基于振幅归一化加权的聚束分析,提出了一种适用于横向非均匀介质模型的多道瑞雷波频散曲线正演方法.首先,基于聚束分析的计算公式推导得到了经振幅归一化加权后输出功率谱中相速度与局部相速度之间的关系,然后通过黄金分割极值搜索算法计算得到了多道瑞雷波数据的理论频散曲线.数值分析结果表明,该算法能够快速地实现横向非均匀介质中多道瑞雷波频散曲线的正演计算,所求取的频散曲线与采用二维弹性波时间域有限差分模拟分析得到的结果误差较小,这在一定程度上说明了该计算方法的可靠性,从而可为面波多道分析法中的观测系统快速优化设计以及横向非均匀介质中频散曲线的反演解释提供理论支撑.
关键词: 聚束分析      横向非均匀介质      面波多道分析法      频散曲线      正演模拟     
Forward modeling of multichannel Rayleigh wave dispersion curve in laterally inhomogeneous media
HU ShuFan1, ZHAO YongHui1, WU JianSheng2,1, GE ShuangCheng3     
1. School of Ocean and Earth science, Tongji University, Shanghai 200092, China;
2. The Shanghai Municipal Committee of Jiusan Society, Shanghai 200041, China;
3. Zhejiang University of Water Resources and Electric Power, Hangzhou 310018, China
Abstract: Multichannel analysis of surface wave (MASW) method has been one of the effective tools for retrieving near-surface shear-wave velocity. As with other geophysical methods, its forward problem is of great importance to the design of the field survey and the solution of the inversion. Currently, the forward algorithms of surface wave dispersion curve are generally classified into two categories: solution of surface wave eigenproblems in a horizontally layered model, which is time-efficient but difficult to consider lateral variations; and wave equation-based full wavefield simulation, which can theoretically deal with any complicated geological model but requires high computational costs. In this paper, we present a forward algorithm of multichannel Rayleigh wave dispersion curves in laterally inhomogeneous media based on amplitude-normalized beamforming. By analyzing the formula of beamforming, the relationship between the local phase velocity and the phase velocity in the power spectrum estimated by the amplitude-normalized beamforming is derived, and then the theoretical dispersion curve of multichannel Rayleigh wave data is calculated by utilizing the Golden Section search method. Numerical results show that the proposed algorithm can quickly and efficiently calculate the theoretical dispersion curve of multichannel Rayleigh waves in laterally inhomogeneous media. The difference between dispersion curves calculated by the proposed algorithm and produced by the 2D elastic wave finite-difference time-domain algorithm is small, which indicates the reliability of the proposed algorithm. The proposed forward algorithm can provide a theoretical foundation for the inversion of dispersion curves in laterally inhomogeneous media and an efficient tool for the optimization design of the MASW acquisition system.
Keywords: Beamforming    Laterally inhomogeneous media    MASW    Dispersion curve    Forward modeling    
0 引言

瑞雷波是由纵波和横波垂直分量在自由表面发生相长干涉而形成的一种沿自由表面传播的地震波.在垂向非均匀介质中,瑞雷波具有几何频散特性,即不同频率分量以不同的相速度传播,进而可用来探测地下介质的速度结构信息,其应用领域包括浅地表地球物理(曾求等,2020宋政宏等,2020林融冰等,2020)、勘探地球物理(孟小红和郭良辉,2007吴曲波等,2019)及区域与全球地震学(顾勤平等,2020付媛媛和肖卓,2020)等.在浅地表面波勘探中,早期采用的面波谱分析法(Spectral Analysis of Surface Waves, SASW)通过测量两道接收信号的相位差来获取频散曲线(Nazarian et al., 1986),但其计算结果存在易受随机噪声、体波以及多模式叠加等影响的问题(Neducza,2007Hashemi et al., 2019).而之后发展的面波多道分析法(Multichannel Analysis of Surface Waves, MASW)将接收到的时空域多道面波信号变换到另一个域中以此来构建频散图像(Park et al., 1999Xia et al., 1999),根据面波的强能量特征可从此图像中提取得到高质量的频散数据(Socco et al., 2010),常用的变换方法有:频率-波数变换(Gabriels et al., 1987Rix et al., 2002)、倾斜叠加变换(McMechan and Yedlin, 1981Xia et al., 2007)以及相移法(Park et al., 1998)等.当前,面波多道分析方法以其非侵入性、无损、高效及经济的特点而越来越受到学术界的重视(夏江海等,2015).

瑞雷波频散曲线的正演是其反演的基础,并可为实际工作调查的设计提供理论指导依据.当前,面波多道分析法中采用的理论频散曲线计算方法大多基于水平层状模型假设,主要有Thomson-Haskell法(Thomson,1950Haskell,1953)、Schwab-Knopoff法(Schwab and Knopoff, 19701972)、δ矩阵法(Watson,1970Buchen and Ben-Hador;1996)、Abo-Zena法(Abo-Zena,1979)、反射透射系数法(Kennett and Kerry, 1979)以及广义反射透射系数法(Chen,1993何耀锋等,2006)等.上述方法通过对频散方程进行搜根可快速实现瑞雷波理论相速度的计算,但实际浅地表介质中往往呈现较明显的横向不均匀性,接收排列下方是否能用水平层状模型来表达还有待商榷.而基于波动方程的全波场模拟理论上可考虑瑞雷波在任意横向变化介质中的传播,模拟过程通常采用有限元(刘雪明等;2009)、有限差分(周竹生等,2007袁士川等;2018)等数值计算方法来实现,其中又以交错网格有限差分算法(Virieux,1986)应用最为广泛.在交错网格有限差分算法中,自由表面边界的处理将对瑞雷波模拟精度产生较大的影响,早期常采用的真空近似法(Boore,1972)存在数值稳定性差和计算精度低的问题(Graves,1996).之后,通过许多学者的研究,分别提出了应力镜像法(Levander,1988)、横向各向同性介质替换法(Mittet,2002)、声学-弹性边界近似法(Xu et al., 2007)和改进的真空处理公式(Zeng et al., 2012)等自由表面边界的处理方法,有效地提高了瑞雷波的模拟精度.在模拟得到全波场记录后,可与实际处理一样采用面波多道分析法来提取瑞雷波频散曲线.然而,在实际正演模拟中往往需要对多组模型及其参数变化的影响进行分析,由于计算精度的要求以及稳定性条件的限制导致该类方法需要较长的计算时间.

鉴于现有频散曲线计算方法受限于水平层状模型假设或全波场模拟计算量较大等问题,本文从面波多道分析法中聚束分析的计算公式入手,通过推导振幅归一化后聚束分析输出功率谱中相速度与局部相速度之间的关系,采用类似于从频散图像中估计多道瑞雷波相速度的方式,利用黄金分割极值搜索算法求得理论相速度值,以实现横向非均匀介质中多道瑞雷波频散曲线的快速正演.并针对两类典型的横向不均匀介质模型,分析了局部频散曲线、平均频散曲线以及本文方法计算得到的理论频散曲线与二维弹性波场模拟分析得到的频散曲线之间的差异.

1 振幅归一化聚束分析

对于主动源瑞雷波勘探中采用的线性排列,假设共有L道接收且每道的偏移距由xl表示、每道接收信号的傅里叶谱由表示,传统聚束分析的输出可表示为(Johnson and Dugeon, 1993)

(1)

其中,j2=-1;ω为角频率,单位为rad/s;k为波数,单位为rad/m;wl为每道的权重.

上式也可写成矩阵形式,表示为

(2)

其中,H代表复共轭转置;e(k)=[e-jkx1, e-jkx2, …, e-jkxl]T为平面波导向矢量;S(ω)=[1(ω), 2(ω), …, L(ω)]T为由每道接收信号傅里叶谱组成的列向量;W为一对角矩阵,包含了每道的权重.

通过式(1)和式(2)可以看出,若取每道的权重wl为1,则利用传统聚束分析得到的频散图像与二维离散傅里叶变换(Gabriels et al., 1987)以及标准线性拉东逆变换(McMechan and Yedlin, 1981Sacchi and Ulrych, 1995潘冬明等,2010)的结果是等价的;若取每道的权重wl=1/||,即将每道信号的振幅作归一化从而补偿几何扩散和衰减的影响,则此情况下聚束分析等价于相移法(Park et al., 1998).

在频散分析过程中,通常计算的是聚束分析的输出功率谱(Rix et al., 2002Zywicki and Rix, 2005),表示为

(3)

其中,R(ω)为空间谱相关矩阵,表示为

(4)

其中,Rmn(ω)= m(ωn*(ω)为第m与第n道之间的互功率谱,*表示复共轭.

若取每道的权重wl=1/||,并经过欧拉公式化简,式(3)可表示为

(5)

其中,Δφmn(ω)为第m与第n道之间的相位差.由式(3)可知,利用式(5)进行频散分析得到的相速度与利用相移法估计得到的相速度仍然是一致的,只是该式计算的是功率谱,而相移法生成的频散图像所表征的是幅度谱.

2 多道瑞雷波频散曲线正演

对于线性排列中的两道瑞雷波记录,经几何扩散校正后第二道信号的傅里叶谱n(ω)可以由第一道信号的傅里叶谱m(ω)表示(Kulesh et al., 2008Askari and Ferguson, 2012),其公式为

(6)

其中,dmn=xn-xm为两道之间的距离;kmn(ω)为两道之间的空间波数,确定了瑞雷波在两接收点之间的传播;Amn(ω)为与频率相关的衰减函数.

在正问题中可将地下模型离散为一系列的网格单元,网格单元的厚度可根据需求设置相同或不同(如图 1a所示),若只考虑弹性介质,则每个网格单元的弹性参数由横波速度、纵波速度以及密度表征.在此模型离散方式下,(6)式可表示为如下形式

(7)

图 1 正演算法流程图 (a) 模型网格剖分示例,黑色▽表示接收点位置;(b) 局部频散曲线;(c) 接收对之间的相位差;(d) 振幅归一化聚束分析的输出功率谱,黑色点线为通过黄金分割极值搜索得到的多道瑞雷波频散曲线. Fig. 1 Workflow for the forward algorithm (a) Example of model discretization. The black ▽ indicates the position of receivers; (b) Local dispersion curves; (c) Phase difference between a receiver pair; (d) Output power of the amplitude-normalized beamforming. The black dotted line represents the multichannel Rayleigh wave dispersion curve obtained by the Golden Section search.

其中,kq(ω)为局部空间波数,由mn道之间第q列网格单元的弹性参数确定;Δxqmn道之间第q列网格单元的宽度.

根据式(7),第m与第n道之间的相位差可表示为

(8)

将(8)式代入(5)式,并根据关系v=ω/k,可得到振幅归一化聚束分析输出功率谱中相速度与局部相速度之间的关系,表示为

(9)

其中,vq(ω)为局部相速度.

在实际多道瑞雷波数据的频散分析过程中,频散曲线通常是根据频散图像中的谱峰值来确定的.同样地,多道瑞雷波数据的理论频散曲线也可通过求取式(9)的局部极大值来求得,其具体实现方法:

(1) 利用传递矩阵算法求得由每一列模型网格单元构成的层状模型的局部相速度vq(ω).

(2) 将局部相速度vq(ω)转换为局部空间波数kq(ω),并根据(8)式计算每接收对之间的相位差Δφmn(ω).

(3) 利用黄金分割极值搜索算法求得式(9)的局部极大值v(ω),计算完所有频率点后即可得到多道瑞雷波数据的理论频散曲线.

正演算法流程如图 1所示.对于复杂地质模型式(9)可能存在多个局部极大值,因此在极值搜索过程中应充分考虑频散曲线的连续性,即将上一个频点的相速度值作为初值进行搜索,从而得到一条与实际频散曲线相匹配的理论频散曲线.

3 理论模型试验

为验证本文所提出的基于聚束分析的多道瑞雷波频散曲线正演算法的可靠性,对典型的横向不均匀介质模型(断层模型、含低速异常模型)进行了试验.首先通过二维交错网格时间域有限差分算法(Virieux,1986),对各向同性弹性模型进行全波场模拟. 全波场数值模拟中采用时间二阶、空间十二阶的差分精度,自由边界条件采用应力镜像法(Levander,1988)进行处理;模型大小为60 m(宽)×15 m(深),横纵向网格剖分为0.1 m×0.1 m;激发震源采用中心频率为20 Hz的雷克子波,时间采样间隔取0.25 ms;接收排列共48道,最小偏移距10 m,道间距1 m.在得到了瑞雷波的波形记录后,采用式(5)对多道瑞雷波数据进行频散分析并根据谱峰值提取相应的频散曲线.之后,采用基于聚束分析的多道瑞雷波频散曲线正演算法对相同模型进行计算得到理论频散曲线,并分析该理论频散曲线与波场模拟分析得到的频散曲线之间的误差.

由于在一定偏移距范围内瑞雷波以柱面波方式传播且此范围内体波能量较强,从而对频散分析造成影响并导致低频分量的相速度值被低估或高估,又称之为近场效应(Zywicki,1999),因此本文采用归一化排列中心距(Normalized Array Center Distance,NACD)和归一化瑞雷波速度(Normalized Rayleigh Wave Velocity,NRWV)来进行误差分析(Yoon and Rix, 2009),其计算公式分别为

(10)

(11)

其中,x为排列中所有接收道偏移距的平均值;λR为瑞雷波的波长;vR(ω)为波场模拟分析得到的瑞雷波相速度;vR, plane(ω)为平面瑞雷波的理论相速度,或称为参考相速度.

3.1 断层模型

图 2a所示,在该模型中间区域存在地层错断,模型中泊松比为常数0.333,每层密度均为1900 kg·m-3,在横向大于40 m的区域内地层横波速度由上到下依次为100、120、170 m·s-1和270 m·s-1.图中红色圆圈为激发震源位置,黑色三角形为接收点位置,黑色虚线代表接收排列中点.图 2b为利用全波场模拟得到的共炮点道集,可以看出在偏移距为20 m左右处瑞雷波记录的视速度发生了明显变化,意味着地下存在较强的横向不均匀性.

图 2 断层模型及相应的单炮记录 (a) 断层模型,红色○和黑色▽分别表示激发震源和接收点的位置,黑色虚线代表接收排列中点;(b) 模拟的共炮点道集. Fig. 2 Fault model and the corresponding single-shot record (a) Fault model. The red ○ and black ▽ represent the positions of source and receivers, respectively; (b) The simulated common-shot gather.

图 3a为利用式(9)计算得到的振幅归一化聚束分析的理论输出功率谱, 图中黑色虚线代表波数k=2π/L,其中L为接收排列长度.在模型离散过程中,每个网格单元的宽度固定为0.1 m、厚度根据模型设定.需要说明的是,在理论频散曲线的正演过程中并不需要像生成频散图像那样计算每个相速度值对应的理论功率谱,而只需进行极大值搜索即可.由图可见,在横向不均匀模型中理论功率谱将存在多个局部极大值,而利用考虑频散曲线连续性的策略可较好地搜索得到一条光滑变化的理论频散曲线(图 3a中白色点线).图 3b为利用式(5)对模拟单炮记录进行波场变换得到的频散图像,图中黑色虚线代表波数k=2π/L.相较于图 3a,此频散图像中的能量分布在一定程度上受到了体波及散射面波的影响.在频率小于5 Hz的范围内,由于低频分量的缺失以及泄露误差的影响,已难以准确地获取瑞雷波相速度,同时在高频段相速度趋于常速度值,因此在接下来的对比中只对5~50 Hz频段范围内的分量进行分析,从频散图像中提取的频散曲线如图 3b中的黑色点线所示.采用基于聚束分析的多道瑞雷波频散曲线正演算法计算得到的理论相速度值如图 3b中的白色点划线所示,可以看出在整个频段范围内该理论频散曲线与提取的频散曲线吻合较好,只是在5~6 Hz频段内相速度值要低于提取值,这主要是由于低频瑞雷波不发育以及近场效应所导致的,理论频散曲线与提取频散曲线之间的均方根误差为4.02 m·s-1,大于6 Hz频段内的均方根误差为0.94 m·s-1.接收排列中点处的局部频散曲线如图 3b中的黑色点划线所示,在所有频段内局部频散曲线的相速度值都要低于提取得到的相速度值,两者之间的均方根误差为6.83 m·s-1,大于6 Hz频段内的均方根误差为5.54 m·s-1,说明当地下存在较强的横向不均匀性时利用面波多道分析法提取的频散曲线并不满足中点假设,将提取得到的频散曲线置于排列中点会给之后的反演解释带入较大的误差.接收排列下方所有局部频散曲线的平均值如图 3b中的蓝色点划线所示,在大于6 Hz频段范围内平均相速度值皆高于提取得到的相速度值,两者之间的均方根误差为6.96 m·s-1,大于6 Hz频段内的均方根误差为6.16 m·s-1,说明面波多道分析法得到的频散曲线并不是排列下方局部频散曲线的简单算术平均.

图 3 断层模型频散分析结果 (a) 基于振幅归一化聚束分析得到的理论功率谱,白色点线为极值搜索得到的理论多道瑞雷波频散曲线,黑色虚线代表波数k=2π/L;(b) 模拟的单炮记录进行波场变换得到的频散图像,黑色点线为提取得到的频散曲线,白色点划线为正演得到的理论频散曲线,黑色点划线为排列中点处的局部频散曲线,蓝色点划线为排列下方的平均频散曲线,黑色虚线代表波数k=2π/L;(c) NACD-NRWV曲线图,黑色点划线、○标注线和×标注线分别为以理论相速度、局部相速度、以及平均相速度作为参考速度计算得到. Fig. 3 Results of dispersion analysis for fault model (a) Theoretical out power of the amplitude-normalized beamforming. The white dotted line represents the multichannel Rayleigh wave dispersion curve obtained by the local maximum search. The black dashed line indicates wavenumber k=2π/L; (b) Dispersion image generated from the simulated common-shot gather. The black dotted line indicates the extracted dispersion curve, the white dot-solid line represents the theoretical dispersion curve, the black dot-solid line is the local dispersion curve at the midpoint of receiver array, and the blue dot-solid line depicts the average phase velocity beneath the entire receiver array. The black dashed line indicates wavenumber k=2π/L; (c) The plot of normalized Rayleigh wave velocity versus normalized array center distance. The black dot-solid line, the ○ marked line, and the × marked line are calculated by the theoretical phase velocity, local phase velocity, and average phase velocity.

图 3c为将理论频散曲线、局部频散曲线和平均频散曲线作为参考相速度计算得到的归一化排列中心距-归一化瑞雷波速度曲线图(NACD-NRWV曲线图).根据水平层状模型中近场效应的研究(Yoon and Rix, 2009Bodet et al., 2009Roy and Jakka, 2017),当地层中横波速度呈速度递增型时,归一化排列中心距大于1的范围内相对误差小于10%;而对于所有类型地层,归一化排列中心距大于2的范围内相对误差小于5%.利用理论频散曲线作为参考值得到的曲线如图 3c中的黑色点划线所示,在归一化排列中心距大于1的范围内最大相对误差为7.8%,在归一化排列中心距大于2的范围内最大相对误差为2.1%,这与水平层状模型中的研究结果一致,进一步说明了本文正演方法的可靠性.利用排列中心点处局部频散曲线作为参考值得到的曲线如图 3c中的圆号标注线所示,在归一化排列中心距大于1和大于2的范围内最大相对误差皆为16.4%,可知在此情况下排列中点假设并不成立.利用排列下方局部频散曲线的平均值作为参考相速度得到的曲线如图 3c中的叉号标注线所示,在归一化排列中心距大于1和大于2的范围内最大相对误差皆为11.0%,相比于理论频散曲线的结果误差值较大.

本文所有算法的计算部分代码都是在Intel Core i3-4150处理器、Windows操作系统下利用Intel Fortran编译器进行编译执行的.对于二维弹性波时间域有限差分算法,由于需要对人工截断边界进行处理并考虑瑞雷波的影响深度(实际模型大小为80 m×60 m),此算例中的计算耗时为390 s,若将网格剖分为0.05 m×0.05 m并同时满足相同采样间隔的要求,则计算耗时上升至3192 s.而对于图 2a所示的层状模型,基于聚束分析的多道瑞雷波频散曲线正演算法无需对模型进行深度方向上的网格细分,在此算例中计算耗时仅为0.5 s,若将网格单元的宽度设置为0.05 m,计算耗时也仅为0.7 s,能够满足多组模型以及模型中参数变化影响的快速计算分析.

3.2 含低速异常模型

图 4a所示,在该模型20~40 m范围内存在一椭圆形低速异常体,该椭圆的长轴和短轴分别为10 m和1 m,模型中泊松比为常数0.333、密度均匀为1900 kg·m-3,地层中横波速度由上到下依次为120 m·s-1、170 m·s-1和270 m·s-1,低速异常体的横波速度为100 m·s-1.图 4b为利用全波场模拟得到的共炮点道集,由于此模型在横向上的不均匀性相较于上个模型要更弱,因此在单炮记录中瑞雷波视速度的变化也更缓一些.

图 4 含低速异常模型及相应的单炮记录 (a) 含低速异常模型,红色○和黑色▽分别表示激发震源和接收点的位置,黑色虚线代表接收排列中点;(b) 模拟的共炮点道集. Fig. 4 Low-velocity anomaly model and the corresponding single-shot record (a) Low-velocity anomaly model. The red ○ and black ▽ represent the positions of source and receivers, respectively; (b) The simulated common-shot gather.

图 5a为利用式(9)计算得到的理论输出功率谱以及利用黄金分割极值搜索得到的理论多道瑞雷波频散曲线(白色点线), 图中黑色虚线代表波数k=2π/L.在模型离散过程中,每个网格单元的宽度固定为0.1 m、厚度根据模型设定.图 5b为利用式(5)对模拟单炮记录进行波场变换得到的频散图像,图中黑色虚线代表波数k=2π/L.同样地,这里只对5~50 Hz频段范围内的分量进行分析,从频散图像中提取的频散曲线如图 5b中的黑色点线所示.基于振幅归一化聚束分析的多道瑞雷波频散曲线正演算法计算得到的理论相速度值如图 5b中的白色点划线所示,在整个频段范围内理论频散曲线与提取的频散曲线一致性较高,并且在此模型中频散图像受近场效应的影响要更弱,差异主要集中在8.5~10 Hz频段范围内,这是由于低速异常体引起的散射瑞雷波导致的,理论频散曲线与提取频散曲线之间的均方根误差为1.51 m·s-1.接收排列中点处的局部频散曲线如图 5b中的黑色点划线所示,由于模型浅层的横波速度无任何横向变化,局部频散曲线与提取频散曲线的差异主要集中在低于18 Hz的频段范围内,此范围内相速度值要低于从频散图像中提取得到的相速度值,两者之间的均方根误差为7.12 m·s-1.接收排列下方所有局部频散曲线的平均值如图 5b中的蓝色点划线所示,平均频散曲线与提取频散曲线的差异也主要是在频率低于18 Hz的范围内,两者之间的均方根误差为3.03 m·s-1.

图 5 含低速异常模型频散分析结果 (a) 基于振幅归一化聚束分析得到的理论功率谱,白色点线为极值搜索得到的理论多道瑞雷波频散曲线,黑色虚线代表波数k=2π/L;(b) 模拟的单炮记录进行波场变换得到的频散图像,黑色点线为提取得到的频散曲线,白色点划线为正演得到的理论频散曲线,黑色点划线为排列中点处的局部频散曲线,蓝色点划线为排列下方的平均频散曲线,黑色虚线代表波数k=2π/L;(c) NACD-NRWV曲线图,黑色点划线、○标注线和×标注线分别为以理论相速度、局部相速度、以及平均相速度作为参考速度计算得到;(d) 利用最小偏移距为1 m、5 m和10 m的单炮记录得到的叠加频散图像,图中线条表征与(b)一致. Fig. 5 Results of dispersion analysis for low-velocity anomaly model (a) Theoretical output power of the amplitude-normalized beamforming. The white dotted line represents the multichannel Rayleigh wave dispersion curve obtained by the local maximum search. The black dashed line indicates wavenumber k=2π/L; (b) Dispersion image generated from the simulated common-shot gather. The black dotted line indicates the extracted dispersion curve, the white dot-solid line represents the theoretical dispersion curve, the black dot-solid line is the local dispersion curve at the midpoint of receiver array, and the blue dot-solid line represents the average phase velocity beneath the entire receiver array. The black dashed line indicates wavenumber k=2π/L; (c) The plot of normalized Rayleigh wave velocity versus normalized array center distance. The black dot-solid line, the ○ marked line, and the × marked line are calculated by the theoretical phase velocity, local phase velocity, and average phase velocity. (d) Stacked dispersion image generated by common-shot gathers with the nearest offset of 1 m, 5 m and 10 m. The representation of plotted lines is the same as (b).

图 5c为将理论频散曲线、局部频散曲线和平均频散曲线作为参考相速度值计算得到的NACD-NRWV曲线图.利用理论频散曲线作为参考值得到的曲线如图 5c中的黑色点划线所示,在归一化排列中心距大于1和大于2的范围内最大相对误差都为4%,该误差对应的相速度点位于受散射瑞雷波影响的频段内.利用排列中心点处局部频散曲线作为参考值得到的曲线如图 5c中的圆号标注线所示,在归一化排列中心距大于1的范围内最大相对误差为15.9%,在大于2的范围内最大相对误差为14.8%,可知在此情况下排列中点假设仍然是不成立的.利用排列下方局部频散曲线的平均值作为参考相速度值得到的曲线如图 5c中的叉号标注线所示,在归一化排列中心距大于1的范围内最大相对误差为6.9%,在大于2的范围内最大相对误差为4.6%.

地下存在横向不均匀性产生的散射瑞雷波将严重影响频散图像的质量(Yilmaz and Kocaoglu, 2012Mi et al., 2017),而通过叠加不同偏移距数据的频散图像可以较好地压制非一致性噪声(Neducza,2007Socco et al., 2009Pasquet and Bodet, 2017).对于相同模型,保持接收排列固定不变,将最小偏移距分别设置为1 m和5 m并进行全波场模拟,然后将1 m、5 m和10 m三种偏移距单炮记录得到的频散图像进行叠加,叠加结果如图 5d所示,图中黑色点线为从叠加频散图像上提取的频散曲线,称之为叠加频散曲线.由图可见,叠加频散图像受散射瑞雷波的影响较弱,基于振幅归一化聚束分析的多道瑞雷波频散曲线正演算法计算的理论相速度值与提取得到的相速度值之间的一致性有所提升,两者之间的均方根误差为0.95 m·s-1,理论频散曲线与叠加频散曲线之间的平均相对误差为0.4%,最大相对误差为2.5%.在频率低于18 Hz范围内,排列中点局部频散曲线的相速度值仍低于提取得到的相速度值,两者之间的均方根误差为6.60 m·s-1,局部频散曲线与叠加频散曲线之间的平均相对误差为2.0%,最大相对误差为12.2%.平均相速度值与提取得到的相速度值之间的均方根误差有所上升,为3.38 m·s-1,两者之间的平均相对误差为1.1%,最大相对误差为7.6%.可以看出,基于振幅归一化聚束分析的多道瑞雷波频散曲线正演算法计算得到的相速度值与波场模拟分析得到的相速度值之间的误差较小,具有较高的可靠性.

4 结论

本文基于振幅归一化聚束分析,提出了一种适用于横向不均匀介质模型情况下的多道瑞雷波频散曲线快速正演算法,并通过与二维弹性波时间域有限差分模拟结果对比验证了该算法的可靠性.理论模型分析揭示,对于横向不均匀地下介质,利用面波多道分析法提取得到的频散曲线是接收排列下方介质的综合反映,但不是接收排列下方局部频散曲线的简单算术平均,更不能认为是接收排列中点处的介质响应.本文算法研究可为面波多道分析法的观测系统优化提供技术支撑,并且将该方法引入频散曲线的反演解释有望克服当前面波多道分析法在横向不均匀介质探测中水平分辨率的限制,提升反演结果的可靠性.从理论上讲,该算法也适用于其他类型导波(如勒夫波、Scholte波等)在横向不均匀介质中的多道频散曲线计算,而如何考虑多模式面波的传播,以及相应的多道频散曲线正演算法,将是更具理论意义及实际应用价值的新探索.

References
Abo-Zena A. 1979. Dispersion function computations for unlimited frequency values. Geophysical Journal International, 58(1): 91-105. DOI:10.1111/J.1365-246X.1979.TB01011.X
Askari R, Ferguson R J. 2012. Dispersion and the dissipative characteristics of surface waves in the generalized S-transform domain. Geophysics, 77(1): V11-V20. DOI:10.1190/geo2010-0330.1
Bodet L, Abraham O, Clorennec D. 2009. Near-offset effects on Rayleigh-wave dispersion measurements: Physical modeling. Journal of Applied Geophysics, 68(1): 95-103. DOI:10.1016/J.JAPPGEO.2009.02.012
Boore D. 1972. Finite difference methods for seismic wave propagation in heterogeneous materials. //Methods in Computational Physics: Advances in Research and Applications. New York: Academic Press, 11: 1-37, doi: 10.1016/B978-0-12-460811-5.50006-4.
Buchen P W, Ben-Hador R. 1996. Free-mode surface-wave computations. Geophysical Journal International, 124(3): 869-887. DOI:10.1111/J.1365-246X.1996.TB05642.X
Chen X F. 1993. A systematic and efficient method of computing normal modes for multilayered half-space. Geophysical Journal International, 115(2): 391-409. DOI:10.1111/J.1365-246X.1993.TB01194.X
Fu Y Y, Xiao Z. 2020. Ambient noise tomography of Rayleigh and Love wave in Northeast Tibetan plateau and adjacent regions. Chinese Journal of Geophysics (in Chinese), 63(3): 860-870. DOI:10.6038/cjg2020N0239
Gabriels P, Snieder R, Nolet G. 1987. In situ measurements of shear-wave velocity in sediments with higher-mode Rayleigh waves. Geophysical Prospecting, 35(2): 187-196. DOI:10.1111/J.1365-2478.1987.TB00812.X
Graves R W. 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences. Bulletin of the Seismological Society of America, 86(4): 1091-1106.
Gu Q P, Ding Z F, Kang Q Q, et al. 2020. Group velocity tomography of Rayleigh wave in the middle-southern segment of the Tan-Lu fault zone and adjacent regions using ambient seismic noise. Chinese Journal of Geophysics (in Chinese), 63(4): 1505-1522. DOI:10.6038/cjg2020N0117
Hashemi J M, Rahnema H, Boaga J, et al. 2019. Application of surface waves for detecting lateral variations: buried inclined plane. Near Surface Geophysics, 17(5): 501-531. DOI:10.1002/nsg.12059
Haskell N A. 1953. The dispersion of surface waves on multilayered media. Bulletin of the Seismological Society of America, 43(1): 17-34.
He Y F, Chen W T, Chen X F. 2006. Normal mode computation by the generalized reflection-transmission coefficient method in planar layered half space. Chinese Journal of Geophysics (in Chinese), 49(4): 1074-1081. DOI:10.3321/j.issn:0001-5733.2006.04.020
Johnson D H, Dudgeon D E. 1993. Array Signal Processing: Concepts and Techniques. Upper Saddle River: Prentice Hall.
Kennett B L N, Kerry N J. 1979. Seismic waves in a stratified half space. Geophysical Journal International, 57(3): 557-583. DOI:10.1111/J.1365-246X.1979.TB06779.X
Kulesh M, Holschneider M, Ohrnberger M, et al. 2008. Modeling of wave dispersion using continuous wavelet transforms Ⅱ: Wavelet-based frequency-velocity analysis. Pure and Applied Geophysics, 165(2): 255-270. DOI:10.1007/S00024-008-0299-7
Levander A R. 1988. Fourth-order finite-difference P-SV seismograms. Geophysics, 53(11): 1425-1436. DOI:10.1190/1.1442422
Lin R B, Zeng X F, Song Z H, et al. 2020. Distributed acoustic sensing for imaging shallow structure Ⅱ: Ambient noise tomography. Chinese Journal of Geophysics (in Chinese), 63(4): 1622-1629. DOI:10.6038/cjg2020N0272
Liu X M, Fan Y H, Zhai J Y, et al. 2009. Numerical simulation of Rayleigh wave zigzag dispersion curves of four typical strata. Chinese Journal of Geophysics (in Chinese), 52(12): 3042-3050. DOI:10.3969/j.issn.0001-5733.2009.12.013
McMechan G A, Yedlin M J. 1981. Analysis of dispersive waves by wave field transformation. Geophysics, 46(6): 869-874. DOI:10.1190/1.1441225
Meng X H, Guo L H. 2007. Using velocity inversion of seismic Rayleigh wave to compute S-wave statics of P-SV wave. Oil Geophysical Prospecting (in Chinese), 42(4): 448-453. DOI:10.3321/j.issn:1000-7210.2007.04.017
Mi B B, Xia J H, Shen C, et al. 2017. Horizontal resolution of multichannel analysis of surface waves. Geophysics, 82(3): EN51-EN66. DOI:10.1190/GEO2016-0202.1
Mittet R. 2002. Free-surface boundary conditions for elastic staggered-grid modeling schemes. Geophysics, 67(5): 1616-1623. DOI:10.1190/1.1512752
Nazarian S, Stokoe K H, Sheu J C, et al. 1986. Near-surface profiling of geotechnical sites by surface-wave method. //56th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 126-129.
Neducza B. 2007. Stacking of surface waves. Geophysics, 72(2): V51-V58. DOI:10.1190/1.2431635
Pan D M, Hu M S, Cui R F, et al. 2010. Dispersion analysis of Rayleigh surface waves and application based on Radon transform. Chinese Journal of Geophysics (in Chinese), 53(11): 2760-2766. DOI:10.3969/j.issn.0001-5733.2010.11.025
Park C B, Miller R D, Xia J H. 1998. Imaging dispersion curves of surface waves on multi-channel record. //68th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1377-1380.
Park C B, Miller R D, Xia J H. 1999. Multichannel analysis of surface waves. Geophysics, 64(3): 800-808. DOI:10.1190/1.1444590
Pasquet S, Bodet L. 2017. SWIP: An integrated workflow for surface-wave dispersion inversion and profiling. Geophysics, 82(6): WB47-WB61. DOI:10.1190/GEO2016-0625.1
Rix G J, Hebeler G L, Orozco M C. 2002. Near-surface Vs profiling in the New Madrid seismic zone using surface-wave methods. Seismological Research Letters, 73(3): 380-392. DOI:10.1785/GSSRL.73.3.380
Roy N, Jakka R S. 2017. Near-field effects on site characterization using MASW technique. Soil Dynamics and Earthquake Engineering, 97(1): 289-303. DOI:10.1016/J.SOILDYN.2017.02.011
Sacchi M D, Ulrych T J. 1995. High-resolution velocity gathers and offset space reconstruction. Geophysics, 60(4): 1169-1177. DOI:10.1190/1.1443845
Schwab F, Knopoff L. 1970. Surface-wave dispersion computations. Bulletin of the Seismological Society of America, 60(2): 321-344.
Schwab F A, Knopoff L. 1972. Fast surface wave and free mode computations. //Methods in Computational Physics: Advances in Research and Applications. New York: Academic Press, 11: 87-180, doi: 10.1016/B978-0-12-460811-5.50008-8.
Socco L V, Boiero D, Foti S, et al. 2009. Laterally constrained inversion of ground roll from seismic reflection records. Geophysics, 74(6): G35-G45. DOI:10.1190/1.3223636
Socco L V, Foti S, Boiero D. 2010. Surface-wave analysis for building near-surface velocity models-Established approaches and new perspectives. Geophysics, 75(5): 75A83-75A102. DOI:10.1190/1.3479491
Song Z H, Zeng X F, Xu S H, et al. 2020. Distributed acoustic sensing for imaging shallow structure Ⅰ: active source survey. Chinese Journal of Geophysics (in Chinese), 63(2): 532-540. DOI:10.6038/cjg2020N0184
Thomson W T. 1950. Transmission of elastic waves through a stratified solid medium. Journal of Applied Physics, 21(2): 89-93. DOI:10.1063/1.1699629
Virieux J. 1986. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics, 51(4): 889-901. DOI:10.1190/1.1442147
Watson T H. 1970. A note on fast computation of Rayleigh wave dispersion in the multilayered half-space. Bulletin of the Seismological Society of America, 60(1): 161-166.
Wu Q B, Pan Z Q, Chen J Y, et al. 2019. Detecting the three-dimensional distribution of shallow basalt by Rayleigh wave method——taking the Sabkhah Ad Dumathah area in Saudi Arabia as an example. Progress in Geophysics (in Chinese), 34(5): 1938-1944. DOI:10.6038/pg2019CC0401
Xia J H, Miller R D, Park C B. 1999. Estimation of near-surface shear-wave velocity by inversion of Rayleigh waves. Geophysics, 64(3): 691-700. DOI:10.1190/1.1444578
Xia J H, Xu Y X, Miller R D. 2007. Generating an image of dispersive energy by frequency decomposition and slant stacking. Pure and Applied Geophysics, 164(5): 941-956. DOI:10.1007/S00024-007-0204-9
Xia J H, Gao L L, Pan Y D, et al. 2015. New findings in high-frequency surface wave method. Chinese Journal of Geophysics (in Chinese), 58(8): 2591-2605. DOI:10.6038/cjg20150801
Xu Y X, Xia J H, Miller R D. 2007. Numerical investigation of implementation of air-earth boundary by acoustic-elastic boundary approach. Geophysics, 72(5): SM147-SM153. DOI:10.1190/1.2753831
Yilmaz O, Kocaoglu A. 2012. Effect of lateral heterogeneity in the soil column on shear-wave velocity estimation by Rayleigh-wave inversion. The Leading Edge, 31(7): 758-765. DOI:10.1190/TLE31070758.1
Yoon S, Rix G J. 2009. Near-field effects on array-based surface wave methods with active sources. Journal of Geotechnical and Geoenvironmental Engineering, 135(3): 399-406. DOI:10.1061/(ASCE)1090-0241(2009)135:3(399)
Yuan S C, Song X H, Zhang X Q, et al. 2018. Finite-difference modeling and characteristics analysis of Rayleigh waves in viscoelastic media. Chinese Journal of Geophysics (in Chinese), 61(4): 1496-1507. DOI:10.6038/cjg2018L0532
Zeng C, Xia J, Miller R D, et al. 2012. An improved vacuum formulation for 2D finite-difference modeling of Rayleigh waves including surface topography and internal discontinuities. Geophysics, 77(1): T1-T9. DOI:10.1190/GEO2011-0067.1
Zeng Q, Chu R S, Zeng X F, et al. 2020. Extraction of high-frequency surface waves on the Xishancun landslide from ambient seismic noises. Chinese Journal of Geophysics (in Chinese), 63(5): 1830-1843. DOI:10.6038/cjg2020N0208
Zhou Z S, Liu X L, Xiong X Y. 2007. A study of Richardson number and instability in moist saturated flow. Chinese Journal of Geophysics (in Chinese), 50(2): 567-573. DOI:10.3321/j.issn:0001-5733.2007.02.030
Zywicki D J. 1999. Advanced signal processing methods applied to engineering analysis of seismic surface waves[Ph. D. thesis]. Atlanta: Georgia Institute of Technology.
Zywicki D J, Rix G J. 2005. Mitigation of near-field effects for seismic surface wave velocity estimation with cylindrical beamformers. Journal of Geotechnical and Geoenvironmental Engineering, 131(8): 970-977. DOI:10.1061/(ASCE)1090-0241(2005)131:8(970)
付媛媛, 肖卓. 2020. 青藏高原东北缘及邻区Rayleigh和Love波背景噪声层析成像. 地球物理学报, 63(3): 860-870. DOI:10.6038/cjg2020N0239
顾勤平, 丁志峰, 康清清, 等. 2020. 郯庐断裂带中南段及邻区基于背景噪声的瑞利波群速度层析成像. 地球物理学报, 63(4): 1505-1522. DOI:10.6038/cjg2020N0117
何耀锋, 陈蔚天, 陈晓非. 2006. 利用广义反射-透射系数方法求解含低速层水平层状介质模型中面波频散曲线问题. 地球物理学报, 49(4): 1074-1081. DOI:10.3321/j.issn:0001-5733.2006.04.020
林融冰, 曾祥方, 宋政宏, 等. 2020. 分布式光纤声波传感系统在近地表成像中的应用Ⅱ: 背景噪声成像. 地球物理学报, 63(4): 1622-1629. DOI:10.6038/cjg2020N0272
刘雪明, 凡友华, 翟佳羽, 等. 2009. 四种典型地层的瑞雷波"之"字型频散曲线数值模拟研究. 地球物理学报, 52(12): 3042-3050. DOI:10.3969/j.issn.0001-5733.2009.12.013
孟小红, 郭良辉. 2007. 利用地震瑞利波速度反演求取P-SV波横波静校正量. 石油地球物理勘探, 42(4): 448-453. DOI:10.3321/j.issn:1000-7210.2007.04.017
潘冬明, 胡明顺, 崔若飞, 等. 2010. 基于拉东变换的瑞雷面波频散分析与应用. 地球物理学报, 53(11): 2760-2766. DOI:10.3969/j.issn.0001-5733.2010.11.025
宋政宏, 曾祥方, 徐善辉, 等. 2020. 分布式光纤声波传感系统在近地表成像中的应用Ⅰ: 主动源高频面波. 地球物理学报, 63(2): 532-540. DOI:10.6038/cjg2020N0184
吴曲波, 潘自强, 陈金勇, 等. 2019. 利用瞬态瑞雷面波法探测浅层玄武岩三维分布——以沙特Sabkhah Ad Dumathah地区为例. 地球物理学进展, 34(5): 1938-1944. DOI:10.6038/pg2019CC0401
夏江海, 高玲利, 潘雨迪, 等. 2015. 高频面波方法的若干新进展. 地球物理学报, 58(8): 2591-2605. DOI:10.6038/cjg20150801
袁士川, 宋先海, 张学强, 等. 2018. 黏弹性介质瑞雷波有限差分模拟与特性分析. 地球物理学报, 61(4): 1496-1507. DOI:10.6038/cjg2018L0532
曾求, 储日升, 曾祥方, 等. 2020. 利用地震背景噪声提取西山村滑坡高频面波信号. 地球物理学报, 63(5): 1830-1843. DOI:10.6038/cjg2020N0208
周竹生, 刘喜亮, 熊孝雨. 2007. 弹性介质中瑞雷面波有限差分法正演模拟. 地球物理学报, 50(2): 567-573. DOI:10.3321/j.issn:0001-5733.2007.02.030