横向非均匀介质多道瑞雷波频散曲线正演

胡书凡， 赵永辉* ， 吴健生2,， 葛双成

1 同济大学海洋与地球科学学院， 上海 200092 2 九三学社上海市委员会， 上海 200041 3 浙江水利水电学院， 杭州 310018

0 引言

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

瑞雷波频散曲线的正演是其反演的基础，并可为实际工作调查的设计提供理论指导依据.当前，面波多道分析法中采用的理论频散曲线计算方法大多基于水平层状模型假设，主要有Thomson-Haskell法(Thomson，1950；Haskell，1953)、Schwab-Knopoff法(Schwab and Knopoff，1970，1972)、δ矩阵法(Watson，1970；Buchen 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 振幅归一化聚束分析

(1)

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

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

Z(ω,k)=eH(k)WS(ω),

(2)

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

P(ω,k)=Z(ω,k)Z(ω,k)H

=eH(k)WS(ω)SH(ω)WHe(k)

=eH(k)WR(ω)WHe(k),

(3)

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

(4)

(5)

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

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

(6)

其中，dm n=xn-xm为两道之间的距离；km n(ω)为两道之间的空间波数，确定了瑞雷波在两接收点之间的传播；Am n(ω)为与频率相关的衰减函数.

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

(7)

其中，kq(ω)为局部空间波数，由m与n道之间第q列网格单元的弹性参数确定；Δxq为m与n道之间第q列网格单元的宽度.

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

(8)

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

(9)

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

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

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

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

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

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

图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.

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)

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，2009；Bodet et al.，2009；Roy 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，2012；Mi et al.，2017)，而通过叠加不同偏移距数据的频散图像可以较好地压制非一致性噪声(Neducza，2007；Socco et al.，2009；Pasquet 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波等)在横向不均匀介质中的多道频散曲线计算，而如何考虑多模式面波的传播，以及相应的多道频散曲线正演算法，将是更具理论意义及实际应用价值的新探索.