弹性波在层状地质结构中的传播特性仿真研究

Simulation Study on Propagation Characteristics of Elastic Waves in Layered Geological Structures

摘要: 利用弹性波进行地质结构探测是地球物理勘探中最主流的方法之一，其主要利用弹性波在不同地质结构中的传播特性进行探测，由于地层介质按物性具有层状结构，因此研究弹性波在不同构成的层状地下半空间的传播特性是弹性波地质勘探的基础。针对弹性波在层状地质结构中的传播问题，本文首先推导了层状地层结构中的基本波动方程，然后基于有限元仿真软件COMSOL建立几种常见的地层结构模型并进行仿真分析。结果表明，只有两层高速覆盖层地质结构和速度递减型三层地质结构，其中弹性波能量能够传播到最下层，其余类型地质结构，弹性波在分界面会发生反射，弹性波基本都集中在表面层介质。本研究能够用于指导弹性波地质探测工程实践，具备一定的工程参考价值。

Abstract: Geological structure detection using elastic waves is one of the most mainstream methods in geophysical exploration. It mainly uses the propagation characteristics of elastic waves in different geological structures for detection. Because the formation medium has layered structure according to physical properties, studying the propagation characteristics of elastic wave in different layered underground half space is the basis of elastic wave geological exploration. For the propagation of elastic waves in layered geological structures, this paper firstly deduces the basic wave equation in the layered geological structure, and then establishes several common stratigraphic structure models based on the finite element simulation software COMSOL and conducts simulation analysis. The results show that there are only two-layer high-speed overburden geological structures and velocity-decreasing three-layer geological structures, in which the elastic wave energy can propagate to the lowest layer, and for other types of geological structures, elastic waves will be reflected at the interface, and elastic waves are basically concentrated in surface layer medium. This research can be used to guide the engineering practice of elastic wave geological exploration, and has certain engineering reference value.

1. 引言

近年来重大公路建设和房屋建设项目越来越多，路基以及地基勘探也成为大家越来越关注的问题。本质上来说，利用地下空间不同地质构造所带来的各种物理力学性质的差异对其进行探测属于地球物理探测的范畴，早期使用的方法主要有电法勘探、磁法勘探、重力勘探、以及地震勘探。其中地震勘探具有较强的分辨能力以及勘探精度，从而得到广泛的应用，是目前最高效和使用最多的探测方法。

地震勘探主要是利用人工产生的地震波(本质上是弹性波)在不同地质构造的地下半空间中传播时，由于地下介质弹性和密度之间的差异，弹性波会发生反射、衍射、折射等，通过地面布置的多个检波器拾取各类反射波或直达波信号，并进行一定的数据处理得出地下空间构造 [1]。地震勘探法目前较多的应用在工程地质一类的勘察、区域地质一类的调查、煤矿资源地质一类的勘测等领域 [2]。利用弹性波对层状地下空间进行探测大体上可分为两条路线，一种是利用瑞雷波在层状介质中的频散特性，对瑞雷波频散曲线进行反演来获取地层中的各层厚度以及横波速度 [3] [4] [5] [6]；一种是利用纵波或横波的反射波的时间信息或者对反射波进行反演来获取地层中的各层厚度和层速度 [7] [8]。近一二十年，许多研究人员对瑞雷波频散曲线的反演算法进行了深入研究和讨论，这其中包括早期出现的局部线性化方法，如杨成林和Forbriger使用最小二乘算法对瑞雷波频散曲线进行反演 [9] [10]。但是，局部线性化方法的缺点也很明显，其严重依赖初始模型，同时，雅可比矩阵的求解精度也会对算法涉及的偏导数的计算、反演结果的好坏产生很大的影响。所以，后续有研究人员提出了一种不同于局部线性化的优化方法，称为“非线性全局优化算法”，如遗传算法 [11]、粒子群算法 [12]、人工神经网络算法 [13]、以及其他非线性算法 [14] 等。使用反射波进行地震勘探研究较少，张银松等通过实验确定了地震反射波法的工作参数 [15] [16]，并研究得出了影响反射波分辨率的因素。总之，无论是瑞雷波法抑或是反射波法，其主要依据弹性波在层状地质结构中的传播特性，通过地面采集的多道弹性波信号来获取地质结构信息。因此，研究弹性波在不同构成的层状地下半空间中的传播特性是弹性波地质勘探的基础。

弹性波在层状地下空间中的传播规律由弹性波波动方程控制，并且有大量学者对其进行了研究。但层状地下空间中的弹性波波动方程形式过于复杂，难以从中看出规律，用于指导工程实践，常用的做法是通过有限差分数值模拟或有限元模拟，得出弹性波传播规律 [17] [18]。袁士川 [19] 等通过数值模拟研究了均匀半空间模型和两层速度递增模型中瑞雷波在黏弹性介质中的传播特性。杨天春 [20] 等采用4 × 10阶高阶交错网格有限差分方法对均匀半空间介质、三层速度递增介质、含软夹层介质进行全波场模拟，有效地压制了数值频散影响，使模拟结果更加真实反映瑞利波本身的传播特性。邹延延 [21] 等基于完全匹配层法吸收边界条件的高阶交错网格有限差分方法，对两层速度递增介质、含硬夹层介质和含软夹层介质的模型进行了数值模拟，解决了网格空间截断边界处的强反射问题以及减少了频散误差。刘超 [22] 采用有限元数值模拟方法分析了含软夹层介质时的弹性波传播规律。李伟华 [23] 通过显式有限元方法研究了含软夹层的层状沉积河谷场地的地震动力响应。

综上可知，已有大量学者采用有限差分法研究了均匀地层结构、二层地层结构以及三层地层结构中的弹性波传播规律。而有限元法通常只对其中一种地层结构进行了弹性波数值模拟研究，因此采用有限元法全面分析多种层状地质结构中弹性波的传播特性就显得尤为重要。本文首先推导层状地下半空间中的波动方程，然后基于有限元仿真软件建立常见的几种层状地质模型进行仿真研究，分析其中不同类型的弹性波的传播规律。本研究对于工程实践中采用瑞雷波或反射波进行地质勘探而言，具有一定的指导意义。

2. 层状地下半空间波动方程

2.1. 基本波动方程

由固体力学的运动微分方程、几何方程(应力–应变关系)以及本构方程(应力–应变关系)可以得到纳维–斯托克斯方程即式(1) [24]：

$\mu {\nabla }^{2}u+\left(\lambda +\mu \right)\nabla \left(\nabla \cdot u\right)+\rho f=\rho \ddot{u}$, (1)

其中，λ、μ为Lamé系数，其只与介质的杨氏模量和泊松比有关；ρ为密度；f为体力；u为位移向量。

根据亥姆霍兹定理，向量场 $u\left(x,t\right)$ 可分解为一标量势函数 $\varphi \left(x,t\right)$ 的梯度与一向量势函数 $\psi \left(x,t\right)$ 的旋度之和 [24]，得到式(2)：

$u\left(x,t\right)=\nabla \varphi \left(x,t\right)+\nabla \times \psi \left(x,t\right),\nabla \cdot \psi =0$, (2)

其中

$\{\begin{array}{l}\varphi \left(x,t\right)=\nabla \cdot W\\ \psi \left(x,t\right)=-\nabla \times W\end{array}$ (3)

依据式(2)，可以将位移场分为两部分，分别用 $u^{\left(1\right)}$ 及 $u^{\left(2\right)}$ 表示，并且有：

$\{\begin{array}{ll}{u}^{\left(1\right)}=\nabla \varphi ,\hfill & u^{\left(2\right)}=\nabla \times \psi \hfill \\ \nabla \times u^{\left(1\right)}=0,\hfill & \nabla \cdot u^{\left(2\right)}=0\hfill \end{array}$ (4)

式(4)第二式表明， $u^{\left(1\right)}$ 、 $u^{\left(2\right)}$ 分别表示位移场的无旋部分、等体积部分，因此也将两部分位移对应的弹性波称为无旋波和等体积波。无旋波也叫纵波(P波)，其波阵面上质点位移方向与波传播方向平行；等体积波也叫横波(S波)，其波阵面上质点位移方向与波传播方向垂直。并且在三维空间中，S波可分为水平偏振的横波(SH波)和垂直偏振的横波(SV波)。选择弹性介质为球对称区域，并采用球面坐标。此时，无旋波的位移为 $u=\nabla \varphi$，将其带入不计体力的纳维方程式(1)并注意到 $\nabla \times u=0$，得式(5)：

${\nabla }^2\varphi =\frac{1}{c_1^2}\ddot{\varphi }$. (5)

其中， $\ddot{\varphi }=\frac{{\partial }^2\varphi }{\partial t^2}$， $c_1=\sqrt{\frac{\lambda +2\mu }{\rho }}$。等体积波的位移可表示为 $u=\nabla \times \psi$，将其带入不计体力的纳维方程式(1)并注意到 $\nabla \cdot u=0$，得式(6)：

${\nabla }^2\psi =\frac{1}{c_2^2}\ddot{\psi }$, (6)

其中， $\ddot{\psi }=\frac{{\partial }^2\psi }{\partial t^2}$， $c_2=\sqrt{\frac{\mu }{\rho }}$。式(5)和式(6)就是均匀线弹性材料中的基本波动方程。

2.2. 层状地下半空间中波动方程

层状地下半空间介质模型如图1所示。第一层上界面为自由表面，第n层为半空间，对于每一层介质，均可以导出式(5)和式(6)形式的基本波动方程。

Figure 1. Layered media model

图1. 层状介质模型

柱坐标系下，每一层的位移场和应力场可以表示为 [25]：

$\{\begin{array}{l}{u}_r=\left[k\phi \left(z\right)+\partial \psi \left(z\right)/\partial z\right]H_0^{\left(1\right)}{}^{\text{'}}\left(kr\right){\text{e}}^{-iwt}\\ u_{\theta }=0\\ u_z=\left[k\psi \left(z\right)+\partial \phi \left(z\right)/\partial z\right]H_0^{\left(1\right)}\left(kr\right){\text{e}}^{-iwt}\end{array}$ (7)

$\{\begin{array}{l}{\tau }_{rz}=2\mu k\left(tk\psi \left(z\right)+\partial \phi \left(z\right)/\partial z\right)H_0^{\left(1\right)}{}^{\text{'}}\left(kr\right){\text{e}}^{-iwt}\\ {\tau }_{\theta z}=0\\ {\tau }_{zz}=2\mu k\left(tk\phi \left(z\right)+\partial \psi \left(z\right)/\partial z\right)H_0^{\left(1\right)}\left(kr\right){\text{e}}^{-iwt}\end{array}$ (8)

其中， $\phi \left(z\right)=A{\text{e}}^{i{\gamma }_{p}kz}+B{\text{e}}^{-i{\gamma }_{p}kz}$， $\psi \left(z\right)=C{\text{e}}^{i{\gamma }_{s}kz}+D{\text{e}}^{-i{\gamma }_{s}kz}$， ${\gamma }_p=\sqrt{\frac{c^2}{V_p^2}-1}$， ${\gamma }_s=\sqrt{\frac{c^2}{V_s^2}-1}$， $t=1-c^2/\left(2V_s^2\right)$。

V_p，V_s：纵波波速和横波波速；φ，χ，ψ：质点的位移势；k：波数；u_r，u_θ，u_z：位移在柱坐标系(r, θ, z)下的分量； ${\tau }_{rz},{\tau }_{\theta z},{\tau }_{zz}$ ：应力在柱坐标系(r, θ, z)下的分量； $H_0^{\left(1\right)}$， $H_0^{\left(1\right)}{}^{\text{'}}$ ：第一类零阶汉开尔函数及其导数；A、B、C、D：待定常数。

推导多层介质中弹性波波动方程的关键是首先得到某一层的位移应力响应，然后依据各层之间各变量的传递关系，推导出其余各层的位移应力响应。其中，应用最广泛的是凡友华 [25] 等人提出的快速矢量传递算法。其基本思想为，依据各层上下界面的位移应力的连续性条件，得到相邻两层介质上界面的位移应力矢量的传递矩阵，最后发现，传递矩阵中仅含有上层介质的参数，因此可以由第一层的位移应力矢量逐层推导，得到任意第n层的位移应力矢量。

3. 仿真模型的建立

3.1. 仿真模型的建立

本研究主要关心地表以下几十米范围内，不同地层结构的弹性波传播情况，在此范围内，地层大致可分为三种结构，即均匀地层结构、两层地层结构以及三层地层结构。对于两层地层结构，本研究建立低速覆盖层以及高速覆盖层共两种地层结构。对于三层地层结构，建立含软夹层、含硬夹层、速度递增型以及速度递减型地层结构共四种地层模型，分别进行仿真模拟，依据仿真结果分析其中弹性波传播情况。其中均匀介质模型几何如图2(a)所示；两层介质模型如图2(b)所示(这里只给出高速覆盖两层介质情况，其余情况模型几何一致，只是介质参数有所不同)；三层介质如图2(c)所示(这里只给出含软夹层三层介质模型几何，其余情况模型几何一致，只是介质参数不同)，模型几何尺寸为300 m × 100 m。

(a) (b) (c)

Figure 2. Simulation model. (a) Homogeneous media model; (b) two-layer media model; (c) three-layer media model

图2. 仿真模型。(a) 均匀介质模型；(b) 两层介质模型；(c) 三层介质模型

3.2. 仿真参数设置

仿真参数包括模型中介质的物性参数以及有限元模拟的参数，后者包括震源定义、时间步长、网格大小等。仿真用到的均匀介质模型参数如表1所示。两层介质模型参数和三层介质模型参数分别如表2、表3和表4所示。

Table 1. Homogeneous medium parameters

表1. 均匀介质参数

Table 2. Two-layer dielectric parameters

表2. 两层介质参数

Table 3. Parameters of three-layer dielectric with soft (hard) interlayer

表3. 含软(硬)夹层三层介质参数

Table 4. Speed increasing (minus) three-layer dielectric parameters

表4. 速度递增(减)三层介质参数

网格划分大小为1 m，采用瞬态研究，时间步长为0.0005 s，总计算时间为1 s，震源采用20 Hz的雷克子波。

4. 仿真结果与分析

4.1. 均匀介质仿真结果

均匀地层中的瑞雷波不发生频散现象。震源用雷克子波激发后，波场快照中可以清楚的看到横波、纵波以及瑞雷面波(DR波)，其中沿水平方向界面稳定传播的是瑞雷面波，体波(S波和P波)则在整个均匀介质内传播。

图3为均匀介质模型的水平分量波场快照图，分别表现的是传播时间为0.1 s、0.2 s、0.3 s时的质点振动情况。图中能清晰看到体波和瑞雷面波，纵波传播速度最快，瑞雷面波传播速度最慢，与均匀模型假设的速度一致。从图中可以看出，在水平分界面处面波的能量最强，体波能量较弱且迅速衰减。图3(b)中，水平分量中瑞雷面波的质点振动在到达某一个深度后会发生翻转，与表面瑞雷面波质点振动方向相反。且在传播时间为0.2 s时，由数值模拟得到的面波传播距离大约为50 m，理论计算得到的传播距离为48.4 m，二者基本一致。

(a) t = 0.1 s (b) t = 0.2 s (c) t = 0.3 s

Figure 3. Snapshot of the numerical simulation wave field of the homogeneous medium model (DP is direct longitudinal wave, S means transverse wave, DR is the direct rayleigh surface wave)

图3. 均匀介质模型数值模拟波场快照(DP为直达纵波，S表示横波，DR为直达瑞雷面波)

4.2. 两层介质仿真结果

图4为低速覆盖两层模型的水平分量波场快照图。上层介质的横波波速要比下层介质的横波波速低，上层介质瑞雷面波波速为242 m/s，其波长为12 m，上层介质的厚度为10 m，大约为一个瑞雷波波长。当上层介质的厚度在一个瑞雷波波长范围内时，可以明显观察到频散现象，从波场快照图可以看到能量大多位于上层介质中，透射到下层介质中的能量很少。由于分层介质中的瑞雷波不能像在均匀介质中一样以同一速度传播，因此在介质分界面处会发生波场转换，使得大部分能量以反射波的形式传播回自由界面不断产生新的瑞雷面波，少部分能量转换成体波透射到下层模型中传播。从图4中可以看到瑞雷面波在转换体波(RSP)和转换瑞雷面波(RSR)出现后有明显的衰减，直达瑞雷面波(DR)之后的转换瑞雷面波主要分布在上层介质中。

(a) t = 0.1 s (b) t = 0.2 s (c) t = 0.3 s

Figure 4. Snapshot of the numerical simulation wave field of the two-layer low-velocity coverage model (RSR is the Rayleigh surface wave generated by the body wave conversion)

图4. 两层低速覆盖模型数值模拟波场快照(RSR为体波转换生成的瑞雷面波)

图5是高速覆盖两层模型的水平分量波场快照图，上层介质瑞雷面波速度879 m/s，其波长约为40 m，上层介质厚度为20 m，小于一个瑞雷波波长，可以明显观测到速度频散。由前面的理论推导可得，当覆盖层波速大于下层介质波速时，瑞雷面波在介质分界面处也会发生波场转换。与低速覆盖两层模型不同的是，高速覆盖两层模型在分界面处发生波场转换后大部分能量都透射到下层介质中传播。图中可以看到瑞雷面波能量迅速减小，大部分能量都分布在下层介质中。此外，由于瑞雷面波与横波在各方面特征都很相近，因此直达瑞雷面波在介质分界面处主要转换成透射横波(CS)在下层介质中传播。

(a) t = 0.1 s (b) t = 0.2 s (c) t = 0.3 s

Figure 5. Snapshot of the numerical simulation wave field of the two-layer high-speed coverage model (CS is transmitted shear wave)

图5. 两层高速覆盖模型数值模拟波场快照(CS为透射横波)

4.3. 三层介质仿真结果

图6为三层软夹层模型的水平分量波场快照图，在第一层与第二次介质分界面上(深度10 m)，由于第一层介质的瑞雷面波波速比第二层介质波速高，类似于高速覆盖两层介质，大部分能量在界面处发生波场转换后都透射到软弱夹层中传播，自由表面瑞雷面波能量不断减弱。在第二层与第三层介质分界面上(深度20 m)，由于第二层介质的瑞雷面波波速比第三层介质波速低，类似于低速覆盖两层介质，大部分能量在界面处发生波场转换并反射回上一层介质继续在软弱夹层中传播，只有少部分能量透射到下层介质中以体波的形式传播。

(a) t = 0.1 s(b) t = 0.2 s(c) t = 0.3 s

Figure 6. Snapshot of the numerical simulation wave field of the soft sandwich model

图6. 软夹层型模型数值模拟波场快照

图7为硬夹层型模型的水平分量波场快照图。由于第一层介质的波速比第二层介质波速低，类似于低速覆盖两层介质，由震源激发的地震波在遇到第一层与第二层介质分界面(深度10 m)时会发生波场转换，能量大部分反射回上层介质又不断产生新的转换瑞雷面波。因此自由界面上直达瑞雷面波的能量逐渐减弱，直达瑞雷面波后是不断生成的转换瑞雷面波。第二层与第三层介质分界面(深度20 m)上，由于第二层介质的波速比第三层介质波速大，大部分能量在界面处发生波场转换后都透射到第三层介质中以体波的形式传播。因此能量主要分布在第一层和第三层介质中，在硬夹层中传播的能量很少。

(a) t = 0.1 s (b) t = 0.2 s(c) t = 0.3 s

Figure 7. Snapshot of the numerical simulation wave field of the hard sandwich model

图7. 硬夹层型模型数值模拟波场快照

图8为三层速度递增型模型的水平分量波场快照图。由于第一层介质的波速比第二层介质波速低，类似于低速覆盖两层介质，由震源激发的地震波能量在遇到第一层与第二层介质分界面(深度10 m)时会发生波场转换，能量大部分反射回上层介质又不断产生新的转换瑞雷面波(RSR)，因此自由界面上直达瑞雷面波的能量逐渐减弱，直达瑞雷面波后是不断生成的转换瑞雷面波。在第二层与第三层介质分界面上(深度20 m)，同样由于第二层介质的波速比第三层介质波速低，传播到第二三层介质分界面处的能量在界面处发生波场转换并反射回第二层介质中传播，在第三层介质中传播的能量很少。

(a) t = 0.1 s (b) t = 0.2 s(c) t = 0.3 s

Figure 8. Snapshot of the numerical simulation wave field of the three-layer velocity-increasing model

图8. 三层速度递增型模型数值模拟波场快照

图9为三层速度递减型模型的水平分量波场快照图。在第一层与第二层介质分界面上(深度10 m)，由于第一层介质的波速比第二层介质波速高，类似于高速覆盖两层介质，由震源激发的地震波大部分透射到第二层介质中传播，自由表面瑞雷面波能量逐渐减弱。在第二层与第三层介质分界面上(深度20 m)，同样由于第二层介质的波速比第三层介质波速高，透射到第二层介质中的能量在分界面处继续透射到第三层介质中传播。

(a) t = 0.1 s (b) t = 0.2 s(c) t = 0.3 s

Figure 9. Snapshot of the numerical simulation wave field of the three-layer velocity decreasing model

图9. 三层速度递减型模型数值模拟波场快照

5. 结论

目前利用人工在地面激发的地震波(弹性波)来进行浅层地质勘探已经成为了地球物理勘探领域的主要方法之一。无论是利用反射波还是瑞雷面波的频散特性来进行地质构造探测，其本质都是利用弹性波在不同的地层结构中传播特性的差异。为研究在地面激发弹性波时，各种不同类型的弹性波在层状介质中的传播特性，本文做了以下几方面的工作：

1) 推导建立了基本波动方程和层状地下半空间中的弹性波传播方程；

2) 基于有限元仿真软件COMSOL建立了几种常见的两层和三层地质模型；

3) 仿真分析得出了弹性波的传播规律，结果表明，对于两层介质，当上覆层为高速覆盖层时，大部分弹性波能够传播到下层介质，对于下层介质有较好的探测分辨率；当上覆层为低速层时，弹性波能量基本都集中在第一层，对于第一层的探测分辨率较好。对于三层介质，当中间层为软夹层介质时，弹性波能量基本都集中在中间的软夹层；当中间层为硬夹层介质使，能量主要集中在第一层和第三层介质；在速度递增介质中，能量主要集中在第一层介质中；在速度递减介质中，弹性波能量可以传播到各层介质中，其三层介质均能较好的分辨出来。

基金项目

国家重点研发计划课题“红层地区典型地质灾害多维多场探测监测装备及预警技术”(2018YFC1504903)；国家山区公路工程技术研究中心开放基金项目“基于弹性波反射法的浅层介质动力学模型及成像算法研究”(GSGZJ-2020-01)。

NOTES

^*通讯作者。