0 引言

面波在早期地震勘探中被视为一种强干扰波，一般具有低速、高能量的特点，陆地勘探中的常见面波类型有瑞利波和勒夫波。随着地震勘探技术的发展，研究人员发现面波在层状介质中具有很强的频散特性，因此提出了一种新的近地表勘探方法——面波勘探，其基本思想是利用面波的频散特征勘查地质结构^[1-2]。面波勘探中的一个重要研究领域是面波的正演模拟，尤其是横向各向同性(TI)介质中的面波模拟^[3-5]。

一直以来，有限差分法是研究体波正演模拟的重要手段，现也被应用到面波的模拟研究中^[6]。相比体波波场的正演模拟，面波波场的正演模拟对自由表面有更严格的要求。自由表面的实现直接影响模拟波场的精度和准确性^[7-8]。

目前，基于标准交错网格(Standard Staggered Grid, SSG)的有限差分法在面波正演模拟中应用最广泛。Mittet^[9]实现了基于弹性波交错网格有限差分法的瑞利波波场模拟，并分析了其准确性和精度；Bohlen^[10]通过对比实验，给出了SSG方法下水平界面瑞利波模拟压制数值频散的频散准则；秦臻等^[11]实现了高精度的有限差分瑞利波模拟并分析了正演结果的频散特征；Zeng等^[12]将多轴完美匹配层技术(M-PML)应用到瑞利波数值模拟，解决了传统PML技术在高泊松比介质中不稳定的问题；邵广周等^[13]基于波场模拟方法分析了瑞利波的频散曲线及多模式能量分布；杜兴忠等^[14]利用GPU并行编程架构实现了基于MC-PML的二维瑞利波正演模拟。

当介质性质变化剧烈或地表起伏较大时，SSG有限差分法的稳定性和精度均较差，且需要内插介质参数影响模拟效果。对此，Saenger等^[15]提出一种稳定性更好且更适用于复杂介质情况的旋转交错网格(Rotated Staggered Grid, RSG)有限差分法。Saenger等^[16]又将RSG有限差分法应用到各向异性介质，证明了其在各向异性介质中进行地震波波场模拟的优势。Bohlen^[10]实现了基于RSG的瑞利波波场模拟，并与传统交错网格有限差分法的瑞利波波场模拟结果对比，认为二者在水平界面的模拟精度相当，在倾斜界面的RSG有限差分法的稳定性和精度更高。但是，Bohlen仅定性对比了波形，没有从多方面定量对比两种方法的模拟效果。

近年来，许多学者将RSG差分方法应用到各向异性介质尤其是TI介质的正演模拟研究，取得了相当不错的效果^[17-18]。在基于SSG的面波模拟中，瑞利波和勒夫波模拟一般分别进行P-SV波和SH波方程模拟。随着面波勘探方法的发展，单一的瑞利波和勒夫波方法已经不能满足勘探需求，研究人员开始发展联合反演方法^[19]和全波形反演方法^[20-21]。此外各向异性介质中的面波波场研究也开始被研究人员所重视^[22]。

本文基于TI介质一阶速度—应力二维三分量波动方程将P-SV波方程和SH波方程联合，利用RSG有限差分法进行二维三分量波动方程模拟，实现瑞利波和勒夫波的同步模拟。并从波前快照、波形曲线和频散曲线三个方面定性和定量地对比RSG和SSG有限差分法在二维各向同性介质均匀半空间中面波模拟的效果和精度。通过分析三种典型层状介质中RSG有限差分法模拟面波地震记录对应的频散曲线与理论频散曲线的吻合程度，以验证文中方法。进一步将基于RSG的面波模拟方法推广到TI介质，实现VTI、HTI和TTI介质的二维三分量面波模拟，并分析TI介质中的面波波形特征、频散曲线特征和地震记录特征，同时进行两层速度递增的TTI介质模型的二维三分量面波模拟。

1 二维三分量速度—应力波动方程

TI介质的二维三分量速度—应力波动方程为^[17]

$ \left\{\begin{array}{l}\frac{\partial {v}_{x}}{\partial t}=\frac{1}{\rho }\left(\frac{\partial {\tau }_{xx}}{\partial x}+\frac{\partial {\tau }_{xz}}{\partial z}\right)\\ \frac{\partial {v}_{y}}{\partial t}=\frac{1}{\rho }\left(\frac{\partial {\tau }_{yx}}{\partial x}+\frac{\partial {\tau }_{yz}}{\partial z}\right)\\ \frac{\partial {v}_{z}}{\partial t}=\frac{1}{\rho }\left(\frac{\partial {\tau }_{zx}}{\partial x}+\frac{\partial {\tau }_{zz}}{\partial z}\right)\\ \frac{\partial {\tau }_{xx}}{\partial t}={C}_{11}\frac{\partial {v}_{x}}{\partial x}+{C}_{16}\frac{\partial {v}_{y}}{\partial x}+{C}_{15}\left(\frac{\partial {v}_{z}}{\partial x}+\frac{\partial {v}_{x}}{\partial z}\right)+{C}_{14}\frac{\partial {v}_{y}}{\partial z}+{C}_{13}\frac{\partial {v}_{z}}{\partial z}\\ \frac{\partial {\tau }_{zz}}{\partial t}={C}_{13}\frac{\partial {v}_{x}}{\partial x}+{C}_{36}\frac{\partial {v}_{y}}{\partial x}+{C}_{35}\left(\frac{\partial {v}_{z}}{\partial x}+\frac{\partial {v}_{x}}{\partial z}\right)+{C}_{34}\frac{\partial {v}_{y}}{\partial z}+{C}_{33}\frac{\partial {v}_{z}}{\partial z}\\ \frac{\partial {\tau }_{yz}}{\partial t}={C}_{14}\frac{\partial {v}_{x}}{\partial x}+{C}_{46}\frac{\partial {v}_{y}}{\partial x}+{C}_{45}\left(\frac{\partial {v}_{z}}{\partial x}+\frac{\partial {v}_{x}}{\partial z}\right)+{C}_{44}\frac{\partial {v}_{y}}{\partial z}+{C}_{34}\frac{\partial {v}_{z}}{\partial z}\\ \frac{\partial {\tau }_{xz}}{\partial t}={C}_{15}\frac{\partial {v}_{x}}{\partial x}+{C}_{56}\frac{\partial {v}_{y}}{\partial x}+{C}_{55}\left(\frac{\partial {v}_{z}}{\partial x}+\frac{\partial {v}_{x}}{\partial z}\right)+{C}_{45}\frac{\partial {v}_{y}}{\partial z}+{C}_{35}\frac{\partial {v}_{z}}{\partial z}\\ \frac{\partial {\tau }_{yx}}{\partial t}={C}_{16}\frac{\partial {v}_{x}}{\partial x}+{C}_{66}\frac{\partial {v}_{y}}{\partial x}+{C}_{56}\left(\frac{\partial {v}_{z}}{\partial x}+\frac{\partial {v}_{x}}{\partial z}\right)+{C}_{46}\frac{\partial {v}_{y}}{\partial z}+{C}_{36}\frac{\partial {v}_{z}}{\partial z}\end{array}\right. $

(1)

式中：v_k（k=x, y, z)为不同方向的振动速度分量；$ {\tau }_{ij} $(i、j=x, y, z)为各个方向的应力分量，当i、j相同时为正应力，其余情况为剪应力；$ \rho $为密度；$ {C}_{uw} $(u、w=1, 2，…, 6)为弹性参数，由弹性参数矩阵确定，可由Thomsen参数及Bond变换计算得到。当Thomsen参数中的各向异性参数均等于零时，计算所得的弹性参数矩阵便是各向同性介质对应的弹性参数矩阵。所以各向同性介质也可被视为TI介质的特殊情况，其弹性参数矩阵结构更加简单，故式(1)同样适用于各向同性介质。

当介质为各向同性时，式(1)中的弹性参数组成的弹性参数矩阵为

$ \left[\begin{array}{cccccc}{C}_{11}& {C}_{12}& {C}_{13}& {C}_{14}& {C}_{15}& {C}_{16}\\ {C}_{21}& {C}_{22}& {C}_{23}& {C}_{24}& {C}_{25}& {C}_{26}\\ {C}_{31}& {C}_{32}& {C}_{33}& {C}_{34}& {C}_{35}& {C}_{36}\\ {C}_{41}& {C}_{42}& {C}_{43}& {C}_{44}& {C}_{45}& {C}_{46}\\ {C}_{51}& {C}_{52}& {C}_{53}& {C}_{54}& {C}_{55}& {C}_{56}\\ {C}_{61}& {C}_{62}& {C}_{63}& {C}_{64}& {C}_{65}& {C}_{66}\end{array}\right]=\left[\begin{array}{cccccc}\rho {V}_{\mathrm{P}}^{2}& \rho \left({V}_{\mathrm{P}}^{2}-2{V}_{\mathrm{S}}^{2}\right)& \rho \left({V}_{\mathrm{P}}^{2}-2{V}_{\mathrm{S}}^{2}\right)& 0& 0& 0\\ \rho \left({V}_{\mathrm{P}}^{2}-2{V}_{\mathrm{S}}^{2}\right)& \rho {V}_{\mathrm{P}}^{2}& \rho \left({V}_{\mathrm{P}}^{2}-2{V}_{\mathrm{S}}^{2}\right)& 0& 0& 0\\ \rho \left({V}_{\mathrm{P}}^{2}-2{V}_{\mathrm{S}}^{2}\right)& \rho \left({V}_{\mathrm{P}}^{2}-2{V}_{\mathrm{S}}^{2}\right)& \rho {V}_{\mathrm{P}}^{2}& 0& 0& 0\\ 0& 0& 0& \rho {V}_{\mathrm{S}}^{2}& 0& 0\\ 0& 0& 0& 0& \rho {V}_{\mathrm{S}}^{2}& 0\\ 0& 0& 0& 0& 0& \rho {V}_{\mathrm{S}}^{2}\end{array}\right] $

(2)

式中V_P、V_S分别为纵波速度和横波速度。

为证明基于式(1)实现P-SV波^[11]和SH波^[21]方程的联合是可行的，将式(2)代入式(1)可得各向同性介质P-SV波、SH波的一阶速度—应力方程

$ \left\{\begin{array}{l}\frac{\partial {v}_{x}}{\partial t}=\frac{1}{\rho }\left(\frac{\partial {\tau }_{xx}}{\partial x}+\frac{\partial {\tau }_{xz}}{\partial z}\right)\\ \frac{\partial {v}_{z}}{\partial t}=\frac{1}{\rho }\left(\frac{\partial {\tau }_{zx}}{\partial x}+\frac{\partial {\tau }_{zz}}{\partial z}\right)\\ \frac{\partial {\tau }_{xx}}{\partial t}=\rho \left[{V}_{\mathrm{p}}^{2}\frac{\partial {v}_{x}}{\partial x}+\left({V}_{\mathrm{p}}^{2}-2{V}_{\mathrm{S}}^{2}\right)\frac{\partial {v}_{z}}{\partial z}\right]\\ \frac{\partial {\tau }_{zz}}{\partial t}=\rho \left[\left({V}_{\mathrm{p}}^{2}-2{V}_{\mathrm{S}}^{2}\right)\frac{\partial {v}_{x}}{\partial x}+{V}_{\mathrm{p}}^{2}\frac{\partial {v}_{z}}{\partial z}\right]\\ \frac{\partial {\tau }_{xz}}{\partial t}=\rho {V}_{\mathrm{S}}^{2}\left(\frac{\partial {v}_{z}}{\partial x}+\frac{\partial {v}_{x}}{\partial z}\right)\end{array}\right. $

(3)

$ \left\{\begin{array}{l}\rho \frac{\partial {v}_{y}}{\partial t}=\frac{\partial {\tau }_{xy}}{\partial x}+\frac{\partial {\tau }_{zy}}{\partial z}\\ \frac{\partial {\tau }_{yx}}{\partial t}=\rho {V}_{\mathrm{S}}^{2}\frac{\partial {v}_{y}}{\partial x}\\ \frac{\partial {\tau }_{yz}}{\partial t}=\rho {V}_{\mathrm{S}}^{2}\frac{\partial {v}_{y}}{\partial z}\end{array}\right. $

(4)

式(1)包含P-SV波和SH波，式(3)、式(4)等价式(1)，所以通过数值求解式(1)能够实现包含瑞利波和勒夫波的波场模拟。

2 RSG有限差分法 2.1 RSG差分格式

相比于SSG速度分量和应力分量都在整网格节点(图 1)，RSG进行了网格旋转，将速度分量定义在整网格节点，应力分量定义在网格中心点（图 2），并将密度和弹性参数定义在计算所需的网格点，避免了计算过程中弹性参数和密度的插值，提高了计算精度^[15]。

RSG方法设计了两套坐标系计算微分，其一是对角线方向的$ \tilde{x}O\tilde{z} $坐标系，其二是坐标轴方向的$ xOz $坐标系，先计算$ \tilde{x}O\tilde{z} $坐标系的微分，再通过坐标转换求解$ xOz $坐标系对应的微分。两套坐标系微分算子之间的计算关系为^[15]

$ \left\{\begin{array}{l}\frac{\partial }{\partial x}=\frac{\sqrt{\mathrm{\Delta }{x}^{2}+\mathrm{\Delta }{z}^{2}}}{2\mathrm{\Delta }x}\left(\frac{\partial }{\tilde{z}}+\frac{\partial }{\tilde{x}}\right)\\ \frac{\partial }{\partial z}=\frac{\sqrt{\mathrm{\Delta }{x}^{2}+\mathrm{\Delta }{z}^{2}}}{2\mathrm{\Delta }z}\left(\frac{\partial }{\tilde{z}}-\frac{\partial }{\tilde{x}}\right)\end{array}\right. $

(5)

式中$ \Delta x、\Delta z $为坐标轴方向空间步长。

2.2 稳定性条件

Saenger等^[15]基于冯诺伊曼方法推导了时间二阶、空间2M阶RSG差分法的稳定性条件，并通过下式指出其在二维和三维情况下稳定性一致

$ \begin{array}{c}\frac{\mathrm{\Delta }t{V}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{\Delta }h}\sum\limits_{m=1}^{M}\left|{C}_{m}^{M}\right|\le 1\end{array} $

(6)

式中：$ \mathrm{\Delta }t $为时间步长；$ \mathrm{\Delta }h $为空间步长；V_max为模型速度最大值；$ {C}_{m}^{M} $为空间差分系数，M为差分系数的数目。

对应时间二阶、空间2M阶精度的SSG差分法的稳定性条件为

$ \begin{array}{c}\frac{\mathrm{\Delta }t{V}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{\Delta }h}\sum\limits_{m=1}^{M}\left|{C}_{m}^{M}\right|\le \end{array}\frac{1}{\sqrt{D}} $

(7)

式中D为空间维度。

对比式(6)和式(7)可以发现，理论上RSG方法比SSG方法的稳定性更强。

2.3 吸收边界

传统PML(Perfectly Matched Layer)方法在求解区域外镶嵌厚度有限的吸收衰减层时，将传播到镶嵌层中的波场按传播方向分为$ x $和$ z $两个部分，并引入衰减系数，衰减与边界垂直的部分，以$ {v}_{x} $分量为例^[22]

$ \left\{\begin{array}{l}{v}_{x}={v}_{x}^{\Vert }+{v}_{x}^{\perp }\\ \frac{\partial {v}_{x}^{\perp }}{\partial t}+{d}_{x}{v}_{x}^{\perp }=\frac{1}{\rho }\frac{\partial {\tau }_{xx}}{\partial x}\\ \frac{\partial {v}_{x}^{\Vert }}{\partial t}+{d}_{z}{v}_{x}^{\Vert }=\frac{1}{\rho }\frac{\partial {\tau }_{xz}}{\partial z}\end{array}\right. $

(8)

式中：$ {v}_{x} $分解为x方向分量$ {v}_{x}^{\perp } $和z方向分量$ {v}_{x}^{\Vert } $两个部分; $ {d}_{x}、{d}_{z} $为镶嵌层中$ x $方向和$ z $方向的衰减因子，其中$ {d}_{x} $在$ z $方向(即上、下边界)的镶嵌层中为零，$ {d}_{z} $在$ x $方向(即左、右边界)的镶嵌层中为零，两者在边角上均为非零^[23]。非零时$ {d}_{x}、{d}_{z} $计算公式为

$ \left\{\begin{array}{l}{d}_{x}\left({h}_{x}\right)=\mathrm{l}\mathrm{g}\frac{1}{R}\frac{2{V}_{\mathrm{P},\mathrm{m}\mathrm{a}\mathrm{x}}}{L}{\left(\frac{{h}_{x}}{L}\right)}^{4}\\ {d}_{z}\left({h}_{z}\right)=\mathrm{l}\mathrm{g}\frac{1}{R}\frac{2{V}_{\mathrm{P},\mathrm{m}\mathrm{a}\mathrm{x}}}{L}{\left(\frac{{h}_{z}}{L}\right)}^{4}\end{array}\right. $

(9)

式中：$ R $为理论反射系数，本文取R=0.000001;V_p，max为镶嵌层内的纵波最大速度；$ L $为吸收层厚度；h_x、h_z为距离求解区域分界的最小距离^[22]。

当地震波场变复杂时，PML边界的吸收效果会变差甚至不稳定。为解决这一问题，Meza等^[24]提出了一种多轴PML(M-PML)方法，分析了其在各向同性介质复杂波场和各向异性介质波场中的适用性和稳定性，并给出了数学意义上的稳定性分析结果，证明了M-PML边界的适用性；Zeng等^[12]验证了不同泊松比介质情况下M-PML吸收边界的稳定性和效果；刘东洋等^[25]将M-PML应用到TTI介质的正演模拟，取得了良好的效果。M-PML是PML方法的进一步发展，它在PML基础上对镶嵌层内分裂的波场进行$ x $和$ z $两个方向的衰减，$ x $和$ z $方向上的衰减因子定义为^[22]

$ \left\{\begin{array}{l}x:{d}_{x}={d}_{x}\left({h}_{x}\right),{d}_{z}={p}^{(z/x)}{d}_{x}\left({h}_{x}\right)\\ z:{d}_{x}={p}^{(x/z)}{d}_{z}\left({h}_{z}\right),{d}_{z}={d}_{z}\left({h}_{z}\right)\end{array}\right. $

(10)

式中：$ {p}^{\left(z/x\right)} $为$ x\mathrm{方}\mathrm{向} $(即左、右镶嵌层)的衰减因子比例系数；$ {p}^{\left(x/z\right)} $为$ z\mathrm{方}\mathrm{向} $(即上、下镶嵌层)的衰减因子比例系数。比例系数根据实际实验确定，文中取$ {p}^{\left(z/x\right)}={p}^{\left(x/z\right)}=0.5 $。面波模拟过程中，上边界不再采用吸收边界进行处理，而是采用自由表面边界。

2.4 自由表面边界条件

面波正演模拟研究中，自由表面边界条件会直接影响模拟波场的精度和准确性。Aki等^[8]给出了自由表面应力条件的理论表达式，二维自由表面边界条件的数学表达式为

$ \left\{\begin{array}{l}{\tau }_{zz}=0\\ {\tau }_{xz}=0\\ {\tau }_{yz}=0\end{array}\right. $

(11)

基于上式发展了诸多自由表面处理方法，其中有直接对自由表面应力进行处理的显式表达法，如应力镜像法(Stress Image Method, SIM)等。还有通过调整自由表面物性参数间接实现自由表面的隐式表达法，文中所使用的方法就属于这一类。Mittet^[9]在前人研究的基础上基于气—固界面物性变化给出了一种调整自由表面物性参数的方法，称为横向各向同性介质替换(MS)法，其二维表达式为

$ \left\{\begin{array}{l}{\tau }_{zz}=0\\ \rho =0.5{\rho }_{{}_{{}_{0}}}\\ \lambda =0\\ \mu =0.5{\mu }_{{}_{0}}\end{array}\right. $

(12)

式中：$ \lambda 、\mu $为自由表面的拉梅常数；ρ为自由表面的密度；$ {\rho }_{0}、{\mu }_{0} $为自由表面以下相邻介质的物性参数。

使用RSG差分法能更简单地实现MS法相应的自由表面条件。由于应力分量定义在同一网格点，自由表面位于速度分量采样点所在的整网格点，故只需要对自由表面的密度进行处理，而不需要调整其他物性参数，自由表面的密度满足^[10]

$ \rho =0.5{\rho }_{0} $

(13)

对于自由表面上方的波场分量，直接赋值为零。波场初始化后，从自由表面开始更新波场，隐式满足了自由表面的应力条件。后文通过数值测试证明了这一处理方式的正确性。

3 数值模拟分析

设计各向同性的均匀半空间模型和层状介质模型验证文中面波模拟方法，然后将模拟方法推广到TI介质，并分析TI介质的面波特征。正演模拟使用的震源函数为雷克子波，主频为40 Hz，延时为30 ms。

3.1 各向同性介质均匀半空间模型

建立各向同性介质均匀半空间模型(模型1)，弹性参数见表 1。针对该模型使用文中所提方法和SSG方法进行面波模拟。使用SSG方法模拟面波时，针对自由表面分别采用SIM、声学—弹性边界近似(Acoustic-elastic Boudary Approach, AEA)法，从波场快照、波形曲线和频散曲线三个角度对比分析三种方法面波模拟的精度。

理论上可以通过同时在式(1)中的$ {v}_{z} $和$ {v}_{y} $施加震源实现瑞利波和勒夫波的同时模拟。模拟区域大小为80 m×40 m，在自由表面的$ {v}_{z} $和$ {v}_{y} $施加震源子波，计算时间步长$ \mathrm{\Delta }t=0.1\mathrm{ }\mathrm{m}\mathrm{s} $，空间步长$ \mathrm{\Delta }h=0.1\mathrm{ }\mathrm{m}。$使用三种方法对模型1模拟，得到$ {v}_{x}\mathrm{、}{v}_{y} $和$ {v}_{z} $在t=0.15 s时的波场快照(图 3)。图中，$ {v}_{x} $和$ {v}_{z} $为瑞利波，$ {v}_{y} $为勒夫波，SSG-SIM和SSG-AEA的$ {v}_{x} $和$ {v}_{z} $基于P-SV波动方程正演模拟得到，$ {v}_{y} $基于SH波方程得到，RSG的三分量均是基于式(1)正演模拟得到。根据文中提出的RSG方法得到的波场快照能清楚观察到各类波形(图 3c)，包括瑞利波(R)、剪切首波(Head)、纵波(P)以及SV波和SH波，其形态符合勘探地球物理原理且与SSG-SIM(图 3a)和SSG-AEA(图 3b)得到的波场快照一致，证明本文提出的方法能够正确地实现含面波的二维三分量波场模拟。由于勒夫波不存在于均匀半空间模型中，后续将利用层状介质证明本文方法能同时实现勒夫波模拟。

面波模拟精度与网格剖分的精细程度有重要关系，引入参数ppw（points per minimum wavelength，面波最小波长$ {\lambda }_{\mathrm{m}\mathrm{i}\mathrm{n}} $与空间步长之比）衡量网格剖分的精度，其与空间步长之间的关系为

$ \begin{array}{c}\mathsf{ppw}=\frac{{\lambda }_{\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{\Delta }h}=\frac{{V}_{\mathrm{S},\mathrm{m}\mathrm{i}\mathrm{n}}^{}}{2{f}_{0}\mathrm{\Delta }h}\end{array} $

(14)

式中：$ {f}_{0} $为震源主频；$ {V}_{\mathrm{S},\mathrm{m}\mathrm{i}\mathrm{n}}^{} $为横波最小速度。为验证本文方法的模拟效果，针对模型1设计了ppw为10、25、50的三组网格剖分参数（表 2）。使用SSG-SIM、SSG-AEA和本文提出的RSG三种方法基于表 2的参数对模型1进行模拟。给定震源深度为1 m，接收点位于自由表面，最小炮检距为10 m，接收道数为70，道间距为1 m，得到三分量炮集地震记录。

以ppw=10，RSG为例，地震记录如图 4所示，可见低速强能量的瑞利波和体波，体波相比瑞利波能量很弱，在地震记录上波形不明显，由$ {v}_{y} $可见SH波(图 4a)，其中$ {v}_{x} $(图 4a)和$ {v}_{z} $(图 4b)为瑞利波炮集记录。抽取炮集中炮检距为30 m的单道记录进行对比，并与基于Cagniard-DeHoop^[26]技术得到的解析解进行对比，如图 5所示。为了定量描述数值解与解析解之间的差异，引入L₂范数误差对数值解的精度进行量化^[10]，计算方法为

$ E=\frac{\sum\limits_{l={l}_{\mathrm{s}}}^{{l}_{\mathrm{e}}}\left[f\right(l\mathrm{\Delta }t)-q{\left(l\mathrm{\Delta }t\right)]}^{2}}{\sum\limits_{l={l}_{\mathrm{s}}}^{l\mathrm{e}}{q}^{2}\left(l\mathrm{\Delta }t\right)} $

(15)

式中：f为数值解；q为解析解；$ {l}_{\mathrm{s}} $、$ {l}_{\mathrm{e}} $为计算的起始采样点和结束采样点。很明显，L₂范数误差对时间差异很敏感，而对振幅差异的敏感性较差。从图 5可看出依据RSG方法得到的波形与面波波形的拟合效果很好，随着网格划分精度的提高，拟合效果越来越好。计算图 5中波形曲线对应的数值解与解析解的L₂范数误差，如图 6所示，可发现RSG方法比SSG-AEA和SSG-SIM方法精度高很多，且在网格划分精度较低时也具有较高的精度。

在面波相关研究中，除了波场和波形特征，频散特征也是研究人员关注的重点，且频散特征是对地震记录的综合反应。通过高分辨率线性拉东变换方法(LRT)提取图 5中对应模拟地震记录的频散能量，再提取地震记录对应的频散曲线，并与快速矢量传递算法正演所得的理论频散曲线对比，计算二者的L₂范数误差(图 7)。由图可见，RSG方法模拟所得地震记录对应的频散能量分布(图 7a左、图 7b左、图 7c左)与理论频散曲线的拟合较好，可见三种方法在不同的网格剖分精度下频散曲线的误差都很小，且RSG方法的误差最小。

3.2 各向同性介质层状模型

各向同性介质均匀半空间模型中波场快照、波形曲线和频散曲线的定性—定量对比，已经证明RSG方法的模拟精度和稳定性明显优于SSG-AEA和SSG-SIM方法。为了进一步验证RSG方法，设计了速度递增、含低速软夹层和含高速硬夹层三种典型层状模型，模型参数分别如表 3~表 5所示。震源加载于自由表面的$ {v}_{z} $和$ {v}_{y} $，取时间步长$ \mathrm{\Delta }t=0.1\mathrm{ }\mathrm{m}\mathrm{s} $，空间步长$ \mathrm{\Delta }h=0.25\mathrm{ }\mathrm{m} $，接收线位于自由表面，最小炮检距为10 m，接收道数为100，道间距为1 m，记录时长为1 s。利用RSG方法模拟得到$ {v}_{x}\mathrm{、}{v}_{y} $和$ {v}_{z} $炮集记录，使用高分辨率LRT提取对应的频散能量图并与根据快速矢量传递算法和广义反射—透射系数算法计算得到的频散曲线进行对比，如图 8~图 10所示。可以发现RSG方法在三种典型层状模型中表现良好，能够同时实现瑞利波和勒夫波模拟，$ {v}_{x} $和$ {v}_{z} $炮集中瑞利波发育良好，十分清晰，频散明显，呈发散扫帚状，符合强能量、低速的特征，且伴有明显的体波，符合地球物理勘探规律，$ {v}_{y} $炮集中勒夫波呈扫帚状，同相轴发育清晰，频散明显，炮集记录对应的频散能量分布与理论频散曲线在基阶和高阶都吻合良好。

3.3 TI介质均匀半空间模型

各向同性介质模型的测试结果，证明了RSG方法在面波模拟中的适用性。本文将其推广到TI介质，设计了VTI、HTI和TTI介质均匀半空间模型，其模型参数如表 6所示，表中$ {V}_{\mathrm{P}0} $、$ {V}_{\mathrm{S}0} $、$ ε $、$ \delta $、$ \gamma $为Thomsen参数。$ {V}_{\mathrm{P}0} $、$ {V}_{\mathrm{S}0} $表示准纵波和准横波沿TI介质对称轴方向的传播速度；$ ε $、$ \delta $、$ \gamma $是描述TI介质各向异性强度的独立无量纲参数，$ ε $和$ \gamma $分别表征纵波和横波传播过程中各向异性强度，$ \delta $是控制垂直介质对称轴方向附近纵波各向异性的参数；$ \beta $、$ \psi $分别为极化角和方位角，是描述TI介质对称轴与观测坐标系角度的参数^[19]。当$ ε $、$ \delta $、$ \gamma $均为零时，TI介质退化为各向同性介质；当$ ε =\delta $时，TI介质为一种实际中常见且十分重要的椭圆各向异性介质，故取$ ε =\delta =0.25 $这一具有代表性的弱各向异性介质。

用RSG方法进行面波模拟，将震源分别加载在自由表面的$ {v}_{z} $和$ {v}_{y} $，取时间步长为$ \mathrm{\Delta }t=0.2\mathrm{ }\mathrm{m}\mathrm{s} $，空间步长为$ \mathrm{\Delta }h=0.25\mathrm{ }\mathrm{m} $，得到t=0.15 s的三分量波前快照（图 11~图 13）。由图 11可发现：VTI介质中波场特征和各向同性介质相似，x、z方向和y方向的波场彼此独立，$ {v}_{x} $和$ {v}_{z} $可见明显的瑞利波(R)、剪切首波(Head)、纵波(P)以及SV波，$ {v}_{y} $中见SH波；不同于各向同性介质波场的是纵波和横波波前面不再是半圆，而是椭圆形。由图 12和图 13可知，HTI介质和TTI介质中x、z方向和y方向的波场相互影响，震源加载于$ {v}_{z} $(图 12a、图 13a)和$ {v}_{y} $(图 12b、图 13b)得到的波场形态相似但能量分布不同，$ {v}_{x} $中瑞利波相比$ {v}_{z} $存在明显的横波分裂现象。

在自由表面上设置接收点，最小炮检距为10 m，道间距为1 m，接收道数为70，记录长度为1 s，得到模型7的三分量炮集地震记录，并提取各自对应的频散能量(图 14), 可见由于介质各向异性，导致地震记录复杂、横波分裂现象明显、横波能量更强，三分量的频散曲线相近，$ {v}_{y} $中存在明显的横波分裂现象，由其频散能量图可知慢横波能量的占比更大。

为进一步探究各向异性对频散特征的影响，设计了三个模型组，分别控制$ ε $、$ \delta $、$ \gamma $变化，模型组参数如表 7~表 9所示。将震源加载于自由表面的$ {v}_{z} $，取时间步长为$ \mathrm{\Delta }t=0.1\mathrm{ }\mathrm{m}\mathrm{s} $，空间步长$ \mathrm{\Delta }h=0.1\mathrm{ }\mathrm{m} $，接收线位于自由表面，最小炮检距为10 m，接收道数为70，道间距为1 m，记录时长为1s，对三个模型组进行面波模拟，得到地震记录后，提取对应的频散曲线并进行对比，如图 15~图 17所示。从图中可以看出TI介质中面波的频散特征主要受$ ε $、$ \delta $参数影响，受$ \gamma $参数的影响较小。

3.4 层状TI介质模型

为进一步研究TI介质中面波的特性，设计了一个两层的速度递增模型，其模型参数如表 10所示。将震源加载于自由表面的$ {v}_{z} $，取时间步长$ \mathrm{\Delta }t=0.1 $ ms，空间步长$ \mathrm{\Delta }h=0.1 $ m，接收线位于自由表面，最小炮检距为10 m，接收道数为100，道间距为1 m，记录时长为1 s，得到三分量地震记录并提取频散能量图（图 18）。图中地震记录中面波呈清晰的发散扫帚状，同相轴清晰，能量占比很高，频散特征明显（图 18上）。其中$ {v}_{y} $存在更复杂的面波与SH波（图 18上），对应的频散能量图中也存在疑似“模式接吻”的现象（图 18下中）。

由频散能量图（图 18下）可见层状介质中存在多阶面波，基阶面波的能量在其中占据主导地位，$ {v}_{y} $中高阶面波能量占比较$ {v}_{x} $和$ {v}_{z} $更高。

4 结论

本文基于TI介质一阶速度—应力二维三分量波动方程，通过旋转交错网格有限差分法结合M-PML边界法，并联合P-SV波方程和SH波方程，实现了包含瑞利波模拟和勒夫波模拟的二维三分量面波模拟。在证明文中模拟方案有效的基础上，将方案应用到TI介质二维三分量面波模拟，通过波场快照、地震记录和频散能量图分析TI介质中的面波特征和各向异性参数对面波频散特征的影响。结果表明：①TI介质面波波场相比各向同性介质变得复杂，普遍存在各向异性导致的极化现象；②VTI介质中面波波场特征与各向同性介质中面波波场特征相似，其中P-SV波和SH波的传播过程互不干扰；③HTI、TTI介质中P-SV波和SH波的传播过程相互影响，面波波场更复杂且与各向同性介质存在较大差异；④描述TI介质各向异性强度的三个Thomsen参数中，与纵波各向异性有关的$ ε $和$ \delta $对面波的频散特征影响较大；⑤在两层速度递增的TTI介质中，面波波场和频散特征变复杂，尤其在y方向上。

本文提出的二维三分量面波模拟方案为各向同性介质中瑞利波和勒夫波的模拟提供了精度更高、稳定条件更宽松、效率更高的模拟方法，实现了用同一个方程模拟瑞利波和勒夫波，为二维全波场地震波模拟甚至是全波场反演提供了重要参考，对TI介质中面波三分量波场的研究更接近真实面波勘探情况，对进一步认识各向异性介质中面波的传播特征具有重要意义。