Analysis of Rayleigh Wave Dynamic Response and Propagation Characteristics in Layered Media
丁玮 DING Wei；陈灿 CHEN Can；聂田 NIE Tian
（武汉工程大学土木工程与建筑学院，武汉 430073）
（School of Civil Engineering and Architecture，Wuhan Institute of Technology，Wuhan 430073，China）
摘要：在均匀半空间介质中，当地表面存在一个竖向振源时，可能产生两种波：体波、面波。体波包括横波和纵波，面波分为瑞利波和拉夫波。当瑞利波垂直自由表面向下传播时，最多只能传播到一个波长的深度，其中大部分能量都集中在1/2波长深度范围之内，由此可利用不同波长的瑞利波来探测不同深度的地层。以瑞利波勘探为代表的浅层地质勘探方法，操作简单，经济高效，对于研究复杂地层条件具有指导性的意义。文中利用ANSYS软件建立了轴对称层状介质模型，并模拟瞬态点或面区域荷载作用下土层的动力响应，同时在计算模型的侧边和底边设置粘滞边界可有效减少反射波的能量。并用MATLAB软件进行数据后处理，采用二维傅里叶变换做频率波数域分析，再现了瑞利波在近地表的传播状态。
Abstract: In a uniform half-space medium, when there is a vertical vibration source on the local surface, two kinds of waves may be generated: body wave and surface wave. Body waves include transverse waves and longitudinal waves, and surface waves are divided into Rayleigh waves and Love waves. When the Rayleigh wave propagates downward perpendicular to the free surface, it can only propagate to a depth of at most one wavelength, and most of the energy is concentrated in the depth range of 1/2 wavelength, so that can use different wavelengths of Rayleigh waves to detect different layers. The shallow geological survey method, which is represented by Rayleigh wave exploration, is simple and economical, and it has guiding significance for the study of complex stratigraphic conditions.In this paper, ANSYS software establishes an axisymmetric layered media model, which simulates the dynamic response of soil under transient or surface area. At the same time, the viscous boundary is set at the sides and the bottom of the model to reduce the reflected wave energy. And it uses MATLAB software for data postprocessing, uses two-dimensional Fourier transform for frequency domain analysis to reproduce the propagation state of Rayleigh wave in the near surface.
关键词：瑞利波；数值分析；频率波数域
Key words: Rayleigh wave；numerical analysis；frequency wavenumber domain
中图分类号：TU435 文献标识码：A 文章编号：1006-4311（2019）27-0178-05
0 引言
十九世纪末，英国学者瑞利在研究半无限体弹性介质的波动问题中，发现了一种与压缩波及剪切波不同的波，后来这种波被称为瑞利波。其具有以下主要特性：①瑞利波的能量主要集中于介质的浅部，而且随着深度的增加，能量迅速衰减；②瑞利波的质点振动呈椭圆极化振动特性，不同于体波的线性极化振动；③在水平分层介质中，瑞利波有多个模态，同一频率，不同模态的相速度不同；同一模态，不同频率相速度也是不同的，即瑞利波相速度会随频率变化而变化，这就是瑞利波的频散特性；④瑞利波的传播特性受介质的分层厚度、数量以及各层剪切波速、泊松比等材料参数的影响。
Lamb研究了不同振源下半无限体弹性介质表面质点位移响应。Wood研究了泊松比为0.25弹性半无限体在简谐点荷载作用下，各类型波阵面位移及振动相对幅值。距振源一定距离后，瑞利波沿介质表面以r-1/2几何衰减，P波、S波以r-2沿表面几何衰减，瑞利波占整个波场能量的67%，而剪切波和压缩波分别占26%和7%，这表明地基土表面瑞利波能量占总能量的主要部分。1988年，吴世明等采用瞬态测试法测试了土层波速，系统地阐述了计算原理以及实际应用方法[1]。1998到2003年间张碧星等研究了层状空间中导波的传播和“之”字型频散曲线的形成机理[2]。2004年，柴华友等通过数值分析，对表面波谱分析计算及影响因素做了系统的研究[3]。2006年，杨生采用有限元方法模拟复杂地层条件下的瑞利波传播特性[4]。2009年，潘东明采用改进的传播矩阵法进行了典型地层的瑞利波频散曲线的计算[5]。2013年，万远收对瞬态瑞利波法研究频散曲线做了系统的探讨[6]。自上世纪80年代以来，我国专家学者对瑞利波的原理及应用做了深入的研究，而且取得了许多实用而有价值的成果。由于瑞利波的传播特性包含着丰富的地质信息，所以对瑞利波的研究也就更具有实际意义，运用实测的瑞利波频散曲线，通过定量分析解释，能够得到各地质层的厚度及弹性波的传播速度，而传播速度的大小直接反映了地层的软硬程度，从而确定地基的持力层。瑞利波勘探可方便地划分出该软弱层的埋深及范围；波速的大小与介质的物理力学参数，如密度、剪切模量、压缩模量、泊松比等密切相关[7]。本文利用ANSYS建立层状介质模拟，消除反射，用Ls-Dyna计算层状介质在表面圆盘荷载下的动力响应，探讨分析不同分层情况下的瑞利波传播特性。
1 理论
1.1 瑞利波的形成
均匀半空间面波的形成有各种解释，较新的观点是从波速度的变化观点论述。在均匀介质自由界面半空间传播，纵、横波因为■的缘故，它们总是分开且独立的传播。但在非均匀介质的自由半空间的条件下，纵、横波在将发生反、透、折射，并产生变异转换的现象。从而形成了非均匀的纵波，其波速 （Vs＜Cs）。此时这些波的波速也在发生变异，尤其是波在转换后产生变化和差异，（一般正常均匀介质中■），在非均匀介质的自由半空间的条件下可能出现■。此时，则就可能出现变异的纵、横波同时到达一个质点，质点的运动就会是一种合成，因而迭加为一种“新”的波，这就叫面波。其运动轨迹如图1所示。图2表示的是瑞利波振动矢量图（组合振动方向和幅度）。
在各向均匀半无限空间弹性介质表面上，当一个圆形基础上下运动时，R波占67%，亦就是说，R波的能量占全部激振能量的2/3，因此利用R波作为勘探方法，其信噪比会大大提高。
1.2 瑞利波的传播特性
面波的波动方程，有三种表达形式，其解都是用复变函数表示的[8]：
LR波传播方向为X，质点位移用u表示，地下方向为Z，质点位移用w表示，随Z而衰减，在Y方向上，质点的位移υ = 0。
质点振动仅在X—Z平面内进行，仅有 X 方向位移分量u，Z 方向位移分量w。令坐标X-Y 平面与自由平面重合，Z轴垂直自由界面下，为了简便，我们仅讨论平面波的情况，波沿X轴方向传播，与Y 轴无关，就成了X-Z的两维图像问题。在LR波传播方向为X的条件下，仅有X方向的位移Ux和Z方向的位移Uz。
为此，可以写出三种波动方程在均匀介质完全弹性半空间条件下面，波传播的表达式：
■（1）
此为一般形式，或
■
（2）
第二种形式，又因为■（3）
式中?自为泊松比，且令■（4）
则式（4）变成■（5）
第三种形式
式中■分别代表纵波波速，横波波速和面波速度。
为了研究问题，假设介质为特殊的泊松体，即λ=μ，则■可以求出■，代入方程可以求解：
■（6）
当?自有变化时，可以得到不同形式的方程
■（7）
式中，u为瑞利波水平位移，W为瑞利波垂直位移。
在均匀各向同性介质中传播的面波速度VR、横波速度VS与泊松比σ的关系为[9]：
■（8）
从式（8）可以看出，在均匀各向同性介质中面波速度与频率无关。当泊松比ν变化时，纵波速度、横波速度随泊松比ν的变化。
1.3 关于面波的频散
在均匀介质中传播的瑞利波速度CR与波的频率f（也可以说与波长λ无关），没有频散出现。
但在非均匀介质中传播，瑞利波得传播速度具有频散特性，在这种条件下传播，瑞利波速发生了频散。此时其波动方程可以记为[10]：
■（9）
由这一方程可以求解CR速度。但此时，CR不仅仅与介质的参数（ρ，CP，CS）有关，而且还与圆频率ω有关；也就是说，不同频率的瑞利波的传播速度是不相同的。这也就是介质表面的疏松覆盖层对瑞利波的影响所在。
在均匀介质中传播的瑞利波速度VR没有频散出现，在非均匀介质中传播具有频散特性，这种特点就是瑞利波工程地震勘察的物理基础。
所谓频散是指面波在介质中的传播速度CR是面波频率 f 的函数■，即面波传播速度CR，随频率 f 改变而变化[11]。
由于不同的频率有很多个相速度，于是各种分振动的不同的相位随波的传播而改变，结果由这些分振动干涉迭加形成的总振动（面波的总振动）的波速也将发生变化。
1.4 瑞利波模态
在分层剪切波速小于底部半无限体剪切波速的水平分层介质中，相同频率瑞利波可能具有不同的波长和传播速度，这表明相同频率的瑞利波有多个不同的传播模态，相速度越高，模态的阶次也越高。除了基阶模态（第一阶模态）外，其他模态存在截止频率，即只有高于截止频率，该模态才会出现。对半无限体上固体层，截止频率分布为
■（11）
这里H是层的厚度，cs1，cs2分别是层和半无限体的剪切波速，n是模态数。
瑞利波模态具有以下特点[12]：基阶模态沿深度衰减较快，影响深度主要在一个波长，波长越长，穿透深度也越深；高阶模态相速度比低阶模态高；在相同频率下，高阶模态的波长比低阶模态长；由于高阶模态穿透深度深，对深层土物理力学参数敏感，结合高阶模态有利于对深层介质参数分析；由于截止频率的存在，模态数量一般会随频率增加；在剖面刚度递增的场地，随着频率增加（波长减小），瑞利波传播局限于表层，模态相速度渐进于表层介质瑞利波速。
2 数值模拟
2.1 平面内瑞利波的数值模拟
本文采用ANSYS有限元分析软件来模拟地基土在圆盘荷载下的动力响应，通过地表面的动态响应来确定成层地基土的传播特性，并且与实测结果相比较，进而验证模型的合理性与准确性。
2.1.1 地层条件
在低应变条件下，土壤假定为弹性的，土层结构为标准的水平分层结构。本文主要采用三种不同地基土模型进行有限元动力数值模拟计算，相关模型参数见表1。
2.1.2 有限元模型及边界条件
R波以平面或圆柱波前展开。因此，可以通过在二维模型中使用平面或轴对称模型来研究R波的传播行为。在非均匀半空间中，正向R波在各个方向上被异质散射。可以根据二维模型的结果分析散射R波的波形和传播行为。例如，散射物体可以被认为是点源的结构。通常，由点或盘状源引起的R波在圆柱形前端扩散。因此，轴对称模型是本研究的首选。模型尺寸深为30m，轴对称半径为30m。在理论计算过程中，一般都是采用的半无限空间介质，但数值模拟无法实现此种情况，非均匀半空间介质在水平和深度方向实际上都是无限大的。然而，在建模过程中还是用有限大小的模型模拟非均质半无限体，因为计算机的存储有限，并且，计算过于耗费时间。 这意味着我们需要在模型中加入一些人工边界。为了模拟无限情况，波不能在这些人造边界处散射。因此，入射波可以传播而不被反射的无反射边界必须应用于模型。对于三维模型，可直接用ANSYS软件命令流中ednb命令来执行，但对于二维模型，使用此命令无效。经过研究，我们也找出了一种方法，具体操作如下：①先编写基于Ansys的指令；②Nsel 待施加的边界的节点component；③在待施加的边界施加作用力Edload；④关闭Slove指令；⑤在Solution中选择转换成*.k文件；⑥在ansys product launcher对*.k 文件进行修改，*，将Define_curve 删除（对施加无反射边界的部分）然后，将*set_node_list,按逆时针排列!可以用逗号，然后用*Boundary_non_reflecting_2D代替；⑦当边界部分节点的排序不是countclockwise，则这部分没有施加无反射边界。或者在边界取一部分节点；⑧生成结果文件，然后用Ansys的后处理；⑨开始新的计算时，须将先前相关的计算结果文件删除。
2.1.3 震源的选择
在实际工程测试中，常采用电机激振器等来激发震源，不同的振源在频率域能量分布不同，信噪比也不同。本文选用的是一种在地震正演模型计算中普遍采用的Ricker子波。Ricker波是高斯函数的归一化二阶导数，定义为：
■（12）
■（13）
式中，fM是峰值频率，tM是波形在时域中的时间偏移量，其中■。
图4中显示了Ricker来源及其频谱。 可以看出，Ricker源由一个正的主瓣和两个负的旁瓣组成。像包络余弦源一样，增加持续时间，峰值频率相应减小。由于能量集中在峰值频率附近，没有任何旁瓣。 在本研究中，选择Ricker源作为脉冲源。
2.2 数值计算结果及分析
均匀各向同性的弹性水平层状半无限体，对称点源位于第一层介质中激发P-SV场。采用传递矩阵方法来计算层状半无限体空间中的弹性波场及瑞利波模式分析。假定所有的接收点都位于水平界面上，而且接收到的信号都为垂直方向的位移分量。对于实验结果，利用多道时域波形采取二维傅里叶变换即可得到频率波数域的表示。在实际激发过程中，二维傅里叶变换的结果具有一定的误差。
选取三种层状介质模型，它们的参数取值见表1，模型1是一个规则的剖面，横波波速随着深度的增加而增加，模型2包含了一个低速层，模型3包含了一个高速层。图5分别记录了三种情况下某一时刻的质点速度幅值等值线图。
所有的接收点在自由界面上成线性排列以接收波形，对每个波形，采样周期为0.3ms。图（a）是三种分层情况下的频率波数域图，图（b）是三种情况下的时间偏移距图。
模型1是一个横波速度随着深度的增加而增加的三层模型。图5表明在源附近，所有三种类型的波（P，S和R波）被叠加在一起，特别是对于S波和R波。图6（b）表示的是60道接收器收到的时间偏移距图，两个相邻的接收点的间距是0.4米，图6（a）表示的是对图6（b）中所得到的信号进行二维Fourier变换得到的在f-k域的波场。可以看出最大值点所对应的频率和波数值。
模型2模型3是分别含有低速夹层和高速夹层的三层介质模型，与模型1不同，由于不同频率范围内各模式的激发强度不一致，会得到之字频散曲线。
3 结论
本文对三种不同分层情况下瑞利波传播模型进行了数值模拟计算，进行了详细分析，得出了如下结论：
①在上软下硬的土层结构中，瑞利波速随波长呈递增趋势；在上硬下软的土层结构中，瑞利波速随波长呈递减趋势。②对于横波波速随深度逐层增加的土层结构，瑞利波速度随波长呈递增趋势；对于横波波速随土层逐层减小的土层结构，瑞利波速随波长呈递减趋势；土层中含有硬[10]夹层的土层结构，瑞利波速先递增到极值，再开始递减。③对于存在软弱夹层的土层，会出现之字形频散曲线，多道瑞利波地震记录图中也表现出瑞利波的高阶模式。
参考文献：
[1]吴世明，曾国熙，陈云敏，叶贵如.利用表面波频谱分析测试土层波速[J].地震工程与工程振动，1988（04）：27-32.
[2]张碧星，肖柏勋，杨文杰，曹思远，牟永光.瑞利波勘探中“之”形频散曲线的形成机理及反演研究[J].地球物理学报，2000（04）：557-567.
[3]柴华友.表面波有效相速度分析方法研究[D].中国科学院研究生院（武汉岩土力学研究所），2004.
[4]杨生.复杂地层条件下瑞利波传播特性的数值分析研究[D].西安理工大学，2006.
[5]潘冬明.瑞雷面波频散分析与应用[D].中国矿业大学，2009.
[6]万远收.瞬态瑞利波法频散曲线的仿真研究[D].华中科技大学，2013.
[7]夏唐代，颜可珍，孙鸣宇.饱和土层中瑞利波的传播特性[J]. 水利学报，2004（11）：81-84.
[8]杨天春，朱自强，周勇.Matlab环境下瑞利波有限差分正演与曲线绘制[J].煤田地质与勘探，2010（01）：62-65，70.
[9]夏唐代，胡永生，杨顺群，陈昭友.道路结构瑞利波特性及动力响应分析[J]. 郑州工业大学学报，2000（01）：19-22.
[10]柴华友，白世伟，刘明贵，卢应发.瑞利波特性刚度矩阵分析方法[J].岩土力学，2006（02）：209-213，218.
[11]柴华友，张电吉，卢海林，杨典森，周春梅.层状饱和介质中瑞利波传播特性薄层分析方法[J].岩土工程学报，2015（06）：1132-1141.
[12]张碧星，鲁来玉.用频率-波数法分析瑞利波频散曲线[J]. 工程地球物理学报，2005（04）：245-255.
[13]陈龙珠，黄秋菊，夏唐代.饱和地基中瑞利波的弥散特性[J].岩土工程学报，1998（03）：6-9.
[14]夏唐代，陈龙珠，吴世明，朱少杰，丁狄刚.半空间饱和土中瑞利波特性[J].水利学报，1998（02）：48-54.
[15]柴华友，汪江波，周一勤，陈星烨.瑞利波分析方法及应用进展[J].岩石力学与工程学报，2002（01）：119-125.
[16]Hakan ■nar,Hamdi Alkan. Crustal S-wave structure beneath Eastern Black Sea Region revealed by Rayleigh-wave group velocities[J]. Journal of Asian Earth Sciences,2016,115.
[17]Laiyu Lu,Chenghao Wang,Bixing Zhang. Inversion of multimode Rayleigh waves in the presence of a low‐velocity layer: numerical and laboratory study[J]. Geophysical Journal International,2007,168(3).
[18]Sri Atmaja P. Rosidi. Comparison Between 2-D and 3-D Stiffness Matrix Model Simulation of Sasw Inversion for Pavement Structure[J]. Civil Engineering Dimension,2007,9(1).
[19]Takanori Harada,Hongze Wang,Shouji Saitou,Norihiko Yamashita,Genji Mori. Characteristics of P-SV Wave Propagation and Displacements by Harmonic Point Load on The Ground Surface of Elastic-Layered Half-Space[J]. Journal of applied mechanics,2005,8.
|
