浅部地下洞穴的瑞利波散射数值分析

王启渊，柴华友，丁 玮

武汉工程大学资源与土木工程学院，湖北 武汉 430074

随着我国经济的飞速发展，基础建设也随之得到了巨大的发展，尤其是地铁、隧道等地下工程。近年来，因地下岩溶导致施工过程中地面塌陷，基础不均匀沉降，造成施工机械设备被埋、人员伤亡等事故时有发生。岩土工程勘查中对地下洞穴的探测通常是大范围的，常规的地质勘察手段难以同时兼顾效率和成本。波动测试是近年来新发展的一种物探方法［1-2］，Stokoe等［3］采用表面波谱 分 析 方 法（spectral analysis of surface wave，SASW）对高速公路路面及路基进行了探测，为瞬态瑞利波法在工程中的广泛应用奠定了基础。

对表面震源，瑞利波能量主导浅部波场，瑞利波遇遭浅部地下异质体会发生散射［4］。表面震源在介质中会产生体波和瑞利波，当介质中有洞穴存在时，体波和瑞利波都会发生散射［5-9］，散射波会产生频散和谱扰动现象，其扰动的程度及形态与散射波中各类型波的能量及传播特性有关。近年来，一些学者尝试利用瑞利波在弹性半空间介质中的传播特性，将瑞利波法应用于探测浅部地下洞穴。卞鹏等［10］将通过研究传感器位置和间距对地下障碍物探测的影响，分析含洞穴半空间介质频散曲线的特征，为探测地下障碍物时频散曲线参数的变化研究做工作。柴华友等［11］通过提取激发震源后的介质表面质点的速度响应，从能量变化的角度分析洞穴对相速度扰动机制。刘中宪等［12］使用间接边界积分方法对入射瑞利波在衬砌隧道中的散射现象进行了分析，发现了衬砌刚度对波散射的影响，同时发现随着隧道埋深的加大，波的散射现象逐渐减弱。邵广周等［13］基于洞穴衍射现象的研究提出了利用广义S转换的方法探测洞穴的埋深。蒋辉等［14］利用散射波现象检验地下洞穴的注浆效果，以波速的变化为判断依据，通过研究表明波速速度及变化幅度越大，注浆的效果越好。张献民等［15］通过有限元法进行数值模拟，总结了瑞利波信号受地下管道的形状和管径变化影响的相应规律。柴华友等［16］通过数值分析的方法分析散射波对表面波相速度的影响，由相速度曲线的变化确定异质体的参数。但相速度扰动影响因素很多，对各种干扰较敏感，相速度扰动与埋深规律难以判断。

基于以上分析，研究散射波能量及传播特性与洞穴埋深及几何参数间关系有助于分析波场扰动特征，推动瑞利波在浅部洞穴探测中的应用。研究洞穴的形态特征主要考虑波在传播过程中遇到洞穴时产生的入射、反射现象，得到散射波场的特征。通过散射波理论可以对介质中的洞穴进行分析，由于解析解以积分形式给出，需迭代求解，求解过程繁琐，不利于研究洞穴散射波基本特征，因此使用数值模拟方法建立二维模型，分析洞穴对瑞利波的散射特征，是十分必要的。将数值模拟二维模型结果进行分析，得到相关特征，用于三维现象的分析中。本文主要从二维模型的角度分析洞穴对瑞利波散射特征。针对这一目的，本文主要开展以下工作：（a）利用有限元模拟计算均匀半空间中存在不同埋深及几何参数洞穴情况下的散射波场；（b）将频率域转换成波长域，在偏移距-波长域振幅谱图中研究波场中的能量变化及传播特性与扰动特征；（c）建立波长与地下洞穴埋深及几何参数间的关系。

1 有限元模型的建立与验证

1.1 模型类型

采用平面及轴对称模型研究瑞利波的散射，复杂边界形状的洞穴可近似用一组小矩形边界洞穴组合建模，因此，采用矩形洞穴模型具有代表性。为了便于分析表面波场，取介质模型长25 m，宽20 m。在波的传播过程中，当传递至人工模型边界时，波会产生反射现象，有可能叠加到行进波中，对分析结果产生误差。为了减少模型边界产生的散射波，在有限元模型边界施加无反射边界，取模型土介质性质参数为：密度 ρ＝1 800 kg/m3，泊松比υ＝0.35，剪切波速Cs=250 m/s，具体模型如图1所示。

图1 模型示意图Fig.1 Schematic diagram of numerical model

1.2 震源

在地震波研究中，如果子波在时域上很窄，则很容易从波场中识别出散射波。因此，狭窄的子波在时域上具有很高的分辨率。在频域中，如果能量集中在一个狭窄的范围内，有利于分析散射波的频率行为。Ricker子波函数具有这些特性，取瞬态冲击荷载变化用Ricker表示，其时域及频率域表达式：

式中t0为峰值对应时间，f=1/t0为中心频率（图中箭头表示位置）。假设 t0=0.01 s，则 fm=100 Hz，Ricker源及其频谱如图2所示。可以看出，Ricker源是由一个正值域和两个负值域组成。

图 2 Ricker函数：（a）时域；（b）频域Fig.2 Ricker function：（a）Time domain；（b）Frequency domain

1.3 模型合理性验证

兰姆给出弹性均匀半空间的瑞利波位移解析解，通过数值模拟Lamb源下响应与解析解比较可验证有限元模型、网格划分及无反射边界有效性。Lamb点荷载在时域的表达式：

频域表达式为：

式中参数P及q用于调整Lamb源的振幅及频域能量分布。

Lamb源的瑞利波的水平（径向）和垂直位移近似表示为：

式中φ=arctan[(t-r/cR)/q]，K和H是与横波、纵波和瑞利波速相关的常数，cR是瑞利波速度，μ是剪切模量，r是距离源的水平（径向）距离（即偏移量）。

峰值t0的Lamb源可写为：

根据傅里叶变换的时间位移规律，对应于上述响应的频谱为：

式（5）、式（6）中的 φ变为：

在有限元模型中，Lamb源参数设置为：P=1，q=1×10-3，t0=16ms。解析解与数值解在偏移距为0.01 m、0.05 m、1 m的对比如图3所示。由图3可知，这两种解的误差很小，表明本文模型、网格划分及无反射边界是合理、有效的。

图3 不同偏移距质点速度响应的数值解与解析解比较Fig.3 Comparison of numerical and analytical solutions of particle velocity responses at different offsets

2 洞穴对瑞利波散射现象

激发表面震源在介质中会产生直达波和瑞利波（R波），介质中若没有洞穴，波不会发生散射现象。图4（a）为表面竖直向质点速度等值云图，可以看出表面波场由瑞利波主导。图4（b）为某一时刻截面质点速度矢量幅值云图，P、S波在介质体内传播，R波沿表面传播，主要能量分布在浅部区域。

图4 均匀半空间中质点速度响应：（a）竖直向质点速度时程图；（b）波场质点速度幅值等值云图Fig.4 Responses of particle velocity in homogenous half space：（a）Time history of vertical particle velocity；（b）Contour plot of particle velocityin wave filed

由于瑞利波能量分布较浅，浅部信号能量主要由瑞利波主导，当浅部出现洞穴时，瑞利波质点轨迹会发生改变，瑞利波发生散射，在时域及频域，散射波对波场扰动，其扰动的程度及形态与散射波中各类型波的能量及传播特性有关。瑞利波的能量通常在一个波长的长度内，且多数能量分布于半波长的长度内，图5给出高频（短波长）与低频（长波长）瑞利波能量在深度方向分布以及能量分布深度与洞穴埋深及几何参数对比。由图5可以分析洞穴对不同波长瑞利波的散射程度。瑞利波与洞穴相遇后，波长小于洞穴埋深的波被洞穴散射的能量很少，大部分能量继续以原有的速度在洞穴上方介质传播。而对于波长大于埋深的波，与洞穴相互作用，一部分能量发生反向散射，一部分能量在洞穴上方透射，在洞穴上方，透射波不仅传播能量会发生变化而且传播速度、衰减规律也会发生变化。均匀半空间中有洞穴的模拟示意图如图6所示，通过与前部分无洞穴情况相对比，讨论瑞利波在遭遇洞穴后产生的散射现象对波场的影响。由图 6（a），（b），（c）不同时刻质点速度幅值云图可以看出入射波在洞穴散射过程，对浅部洞穴，入射波主要是R波，在波场中R波会发生反射、波的模式转换、透射等现象。部分R波经过波的模式转换转化为横波和纵波，还有一部分R波会直接穿透洞穴，形成透射波。其表面竖直向质点速度幅值见图7，图7和图4相比，在洞穴近界面和远界面处有明显的变化，无洞穴的介质中，瑞利波沿传播方向传播至模型底部被无反射边界吸收，无明显能量衰减，而在有洞穴介质中传播，在洞穴边界处有明显的散射波出现，这说明使用散射瑞利波探测地下洞穴是可行的。

图5 瑞利波能量分布深度与波长示意图Fig.5 Illustration of energy distribution depth and wavelength of Rayleigh wave

图6 不同时刻竖直向质点速度幅值云图：（a）波传播至空洞近界面；（b）波传播至空洞上方；（c）波传播至空洞远界面Fig.6 Cloud contours of vertical particle velocity at different times：（a）Waves before the cavity；（b）Waves over the cavity；（c）Waves after the cavity

图7 表面竖直向质点速度云图Fig.7 Contour of surface vertical particle velocity

3 洞穴埋深及几何参数对波散射影响

3.1 散射波场分析

对表面各点竖直向质点速度响应作一维傅里叶变换得到频谱图，为了便于与深度量纲一致及消除波几何扩散导致波衰减，将每个质点响应的频率域谱转化为波长域谱并将幅值进行归一化得到偏移距-波长域振幅谱图，通过对振幅谱能量变化分析预估洞穴的形状，当均匀半空间中存在洞穴时，由前文可知瑞利波在遭遇洞穴后会发生散射现象，散射波在不同的区域表现出不同的特征。

图8（a）给出介质无洞穴情形下的表面波场谱图，不同位置谱能量相似，波场谱均匀，无扰动。当介质出现洞穴时，表面波场谱图如图8（b）所示，在此情况下，波在遇到洞穴前传播状态与均匀介质类似，当波传递至洞穴近界面处（b图中Ⅰ处）时，瑞利波产生反射、绕射、透射等现象，导致在洞穴的近界面有明显的能量区分界线。反射、透射使能量产生损失，同时洞穴近界面和顶面的相接拐点处（b图中Ⅱ处）会使能量产生二次散射，导致在拐点处能量紊乱形成不规则的波场。不仅是洞穴顶部的拐点，底部与远近界面的拐点也会产生这种现象，由于底部拐点，大部分能量向下散射，底部拐点附近散射波对表面波场的影响很小，导致很难判断洞穴的深度。当波行进至洞穴上方（b图中Ⅲ处），波在介质表面和洞穴顶层来回反射，反射波相互干涉，能量或加强或抵消。当上方透射波到达远交界的拐点（b图中Ⅳ处）同样会产生散射现象，波的散射和绕射使上方能量明显减少。在远界面附近（b图中Ⅴ处），透射波随着传播距离增加，传播特性趋于入射瑞利波，能量分布深度逐渐趋于未遭遇洞穴前。通过在振幅能量谱图中的特点，可以判断出洞穴的大致特征。

图8 偏移距-波长域表面竖直向质点振动速度振幅谱图：（a）无洞穴；（b）有洞穴Fig.8 Amplitude spectra of vertical particle velocity of surface in offset-wavelength domain ：（a）Without cavity；（b）With cavity

3.2 散射波对表面波场扰动分析

下面以洞穴埋深、位置及几何尺寸不同的几种算例来验证以上分析，主频 fm=100Hz时主要瑞利波能量分布在波长0.5 m～2.5 m范围，见图8（a），为了避免埋深较小或较大洞穴的微弱的散射波被数值振荡和频谱泄漏引起的噪声掩没，取洞穴的埋深1.5 m～3 m，洞穴具体参数如表1所示，根据不同算例绘制偏移距-波长域质点振动速度振幅谱图，如图9所示。

基于波可能产生的散射行为，将研究区域分为洞穴前方（震源至洞穴近界面）、洞穴上方、洞穴后方（洞穴远界面至远场）三个区域。通过前文所述的洞穴形态判断方法对由洞穴造成波场扰动的三个区域波场，建立波长与洞穴埋深的关系。

1）洞穴前方：波未传播至洞穴近界面时，其传播方式与均匀半空间中相同。当波与洞穴近界面遭遇，发生散射现象，洞穴近界面区域的反射波与

前行波相互干涉。算例 a、b、c、f研究洞穴高度、水平位置相同时埋深参数对前方散射波的影响，对应谱图中特征并结合前文对洞穴近界面洞穴的判断方法发现，前方出现条纹状衰减，产生形状规则的能量分界面。研究表明，不同埋深的洞穴前仍存在能量扰动，出现类似于能量加强或递减的变化，这表明，只要洞穴位于瑞利波的能量影响深度内，不同埋深都会产生散射现象。算例c和d为埋深、水平位置相同但洞穴高度发生变化时洞穴对瑞利波散射的影响，将对应谱图进行对比分析发现，洞穴高度增加，但是高度的下半部分散射向表面的能量较小。同时，在底部与近界面拐点处，波大部分能量向介质体内散射，一部分经过绕射至洞穴底部，向表面波场散射较少。因此，高度的影响对表面波场相对较小。算例c和e为同埋深同高度仅洞穴的水平位置变化时洞穴对瑞利波散射的影响，通过对应谱图中能量变化情况结合前文所述判断方法，谱图中能量的变化趋势是类似的，发现随着洞穴水平位置的变化，只有能量变化区域的水平位置产生变化。

表1 洞穴模型位置和预估位置参数Tab.1 Actual and estimated locations of cavities

图9 偏移距-波长域质点振动速度振幅谱图（矩形实线表示空洞实际大小，虚线表示预估形状，Ⅰ，Ⅱ，Ⅲ分别表示顶面，近界面，远界面）Fig.9 Amplitude spectra of particle velocity in offset-wavelength domain（solid lines：actual shapes of cavities；dotted lines：estimated shapes，Ⅰ，Ⅱ and Ⅲ represent the top，near and far boundaries respectively）

2）洞穴上方：瑞利波通过洞穴近界面区域后，到达洞穴上方，继续发生散射现象，散射波在自由表面和顶面之间来回反射。通过算例a、b、c、f研究同洞穴尺寸、水平位置但埋深不同时，洞穴上方波散射现象。当空洞在介质浅部时（算例a，埋深1.5 m），由于洞穴的埋深与介质表面距离较近波长较大的散射波无法通过，瑞利波具有波长与频率成反比的特性，故能通过的波具有较高的频率。由前文震源选取部分的描述可知，Ricker子波函数在高频部分的频率较集中，因此反射波在埋深与介质表面间发生干涉现象，在谱图a中表现为洞穴顶部后半部分的能量明显强于前半部分。当洞穴埋深加深后，与介质表面距离增大，干涉现象不再发生，故在算例b、c、f中这种特殊现象消失，在谱图中受洞穴影响发生的能量变化特征与前文所述判断方法一致。将算例c分别与算例d和e对比，发现洞穴的高度和质中水平位置的变化对洞穴导致的散射现象影响不大。

3）洞穴后方：上方透射波在洞穴远界面继续发生散射，远界面附近的表面波场同样产生扰动。当波离开洞穴远界面区域后，其传播行为逐渐趋于入射瑞利波传播特性。在埋深较低时（1.5 m处），洞穴远界面与顶部的拐点处没有出现前文中洞穴形态预估中拐点处能量衰减区域，而是出现条纹状的能量加深区域。这是由于顶部与介质表面间距离过短，散射波在洞穴上方发生干涉后在拐点处继续发生干涉，在拐点处积聚也会使拐点处的能量加强，能量在绕过拐点后部分绕射至洞穴底部，部分继续传播。随着埋深增加，没有干涉现象产生，因此在算例b、c、f中都没有出现算例a中的现象。将算例c和d相对比，当洞穴的埋深和水平位置不变高度变化时，与洞穴前的现象类似，洞穴的高度增加后，对表面波场造成的扰动不敏感，波在传播至远界面与底部交汇拐点时，大部分能量散射进洞穴内，一部分向下反射，一部分继续传播，传播模式逐渐变为均匀半空间传播方式。算例c和e为洞穴埋深和高度相同，水平位置变化的情况，对应谱图进行分析，能量变化区域随水平位移变化改变，但能量变化特征与前文所述一致。

3.3波场扰动与洞穴埋深及几何参数的关系

波场能量扰动与洞穴埋深、几何参数有关，根据能量变化程度可确定临界扰动，由临界扰动对应的波长可预估埋深。根据能量幅度变化形成的近界面、顶面和远界面，可以判断洞穴的形状，建立λ与h的关系。预估洞穴位置参数如表1所示。

从波散射引起的扰动和能量模式的变化，可定位洞穴的远近位置。但是从谱分析图中，对于洞穴的高度很难估计。部分情况洞穴的远近界面的估计位置比实际值大，这种情况可能是由于前向散射波在到达某个距离后干扰波的传播，使反向散射波相对于远近界面的前向散射波较弱。

通过洞穴模型位置和预估位置参数表及偏移距-波长域振幅谱图，发现在洞穴的实际埋深（h）与波长（λ）存在一定的对应关系，各算例中无量纲几何参数h/λ如表2所示。

表2 数值模拟中无量纲参数h/λ表Tab.2 Dimensionless parameters h/λ in numerical simulations

在埋深1.5 m～3 m的范围内，无量纲几何参数h/λ的变化范围为 1.367～1.657，算例 a、b、c、f的数据表明随深度变化，h/λ存在一定幅度的增长，但仍以1.5为中心浮动，算例c、d、e表明洞穴在大小和位置等影响因素下h/λ的值基本没有变化，可以判断洞穴的大小和在介质中的位置对h/λ没有影响。综上所述，在介质浅部存在空洞穴时，其扰动无量纲几何参数h/λ介于1.3～1.6之间，即洞穴的埋深与偏移距-波长域振幅谱图中洞穴对应波长的关系以1.5倍为中心上下浮动。

4 结 语

对均匀半空间介质中瑞利波的产生及传播特征进行了分析，运用有限元软件建立含洞穴均匀半空间模型，通过与解析解比较，验证了有限元模型可行性，通过对不同埋深、几何参数算例分析，得到如下结论：

1）当介质中存在洞穴时，瑞利波传播至洞穴位置会发生散射，散射波导致表面波场能量发生扰动。

2）通过数值模拟波场质点速度响应分析了瑞利波在洞穴前方、上方、后方的散射现象及散射波对波场造成的扰动，由此从偏移距-波长域振幅谱图的能量扰动变化判断洞穴的形态，建立了波长与洞穴埋深的关系。

3）在偏移距-波长域，由振幅谱扰动，根据埋深与扰动特征波长的倍数关系可以预测洞穴埋深。

［1］ 栾明龙，魏红，林万顺.瞬态瑞利波技术在工程勘察中的应用［J］. 物探与化探，2012，36（5）：878-883.

［2］ 孙笑，赵明阶，汪魁，等.波动理论在岩土工程测试中的应用研究进展［J］.重庆交通大学学报（自然科学版），2013，32（1）：58-62.

［3］ STOKOE K H，NAZARIAN S.Effectiveness of ground improvement from spectral analysis of surface waves［C］//Proceeding of the eighth European conference on soil mechanics and foundation engineering，Helsinki，Finland.Netherlauds：AA Balkema Publishers，1983：91-95.

［4］ 刘优平，龚敏，徐斌.半空间饱和土中输水管道对瑞利波的散射［J］.四川大学学报（工程科学版），2012，44（5）：71-77

［5］ 卞鹏，王媛，王笃强.有障碍物地基瑞利波频散曲线的特征研究［J］.山东大学学报（工学版），2012，42（10）：99-103.

［6］ CHAI H Y，GOH S H，PHOON K K，et al.Effects of source and cavity depths on wave-fields in layered media［J］.Journal of Applied Geophysics，2014，107（8）：163-170.

［7］ CHAI H Y，GOH S H，PHOON K K，et al.Effects of source and cavity depths on wave fields in homogeneous half spaces［J］.Journal of Applied Geophysics，2013，93（2）：52-59.

［8］ CHAI H Y，PHOON K K，GOH S H，et al.Some theoretical and numericalobservations on scattering of Rayleigh waves in media containing shallow rectangularcavities［J］.Journal of Applied Geophysics，2012，83（6）：107-119.

［9］ CHAI H Y，PHOON K K，WEI C F，et al.Analysis of effects of active sourceson observed phase velocity based on the thin-layer method［J］.Journal of Applied Geophysics，2011，73（1）49-58.

［10］ CHAI H Y，WEI C F.Effects of cavities on the spectra of Rayleigh waves numerical aspects［C/CD］//HUANG A B，MAYNE P W.Geotechnical and Geophysical Site Characterization.Proceeding of the third International Conference on Site Charactorization，Taipei，Taiwan，April，1-4，2008.London，UK：Taylor&Francis Group，2008：807-711.

［11］ 柴华友，卢应发，李祺，等.半无限体介质内异质体方位探测数值研究［J］. 岩土力学，2007，28（1）：188-192

［12］ 刘中宪，梁建文，张贺.弹性半空间中衬砌隧道对瑞利波的散射［J］.岩石力学与工程学报，2011，30（8）：1627-1637.

［13］ SHAO G Z，TSOFLIAS G P，LI C J.Detection of near-surface cavities by generalized S-transform of Rayleigh waves［J］.Journal of Applied Geophysics，2016，129：53-65.

［14］ 蒋辉，柳超，贾超，等.地震散射方法在注浆效果检测中的应用［J］.工程地球物理学报，2015，12（2）：260-265.

［15］ 张献民，支雁飞.地下管道瞬态瑞雷波信号响应分析［J］. 工程勘察，2015，43（6）：81-87.

［16］ 柴华友，崔玉军，卢应发，等.地下洞穴对表面波表观相速度影响的数值分析［J］.岩石力学与工程学报，2006，25（5）：956-962.