0.   引言

瑞雷面波法作为一种工程勘察手段，由于具有抗干扰能力强，垂向分辨率高和野外数据采集相对简便等优点，在探查覆盖层厚度、划分松散地层沉积层序、划分基岩的风化带、探测地下隐埋物体、古墓遗址、洞穴和采空区等方面得到广泛应用.铁路路基坍塌时有发生，面波在铁路路基浅层结构探测中得到了初步应用，据2005年到2020年的全国地质通报，地质灾害中滑坡占比高达71.1%（刘谢攀等，2024），面波在探测滑坡体的滑动带和滑动面起伏形态方面成为未来的研究重要领域之一，因此，开展面波高分辨率勘探研究具有重要的现实意义.20世纪60年代，美国人提出面波的半波长解释方法，并将稳态面波法首先用于地基勘察.较早将稳态面波法形成探测系统用于工程实践的是日本VIC公司，他们经过多年努力，于80年代初推出GR810佐藤式全自动地下勘查机，由于设备昂贵，我国讫今仅购置二台.80年代后期，稳态面波法试验研究在我国悄然兴起，地矿部、水利水电部、冶金部、中科院、浙江大学等均先后开展了应用开发研究.进入90年代，稳态面波法，特别是瞬态面波法，在硬件研制和软件开发两个方面，都相继取得了引人注目的进展.黄嘉正等（1995）采用河北物探队仪器研究所生产的RL-1型稳态激振系统，配备中科院武汉岩土力学所研制的RSM-16H动测仪，应用在深圳机场停机坪软土地基加固效果检测中.雷旭友等（2008）在2008年采用英国GDS公司研制并生产的一种新型稳态面波勘探设备连续表面波测试系统（CSWS），用于软土勘察中.由于稳态激振系统笨重且勘探深度较浅，于是从90年代中期开始，多道瞬态面波分析法（Multi-channel Analaysis of SurfaceWave，简称MASW）在国外开始兴起（Al⁃ Hunaidi，1994；Stockwell et al.，1996），同时国产的多道瞬态面波分析法的仪器和软件也得到了应用（刘云祯和王振东，1996），后来MASW方法逐步取代了稳态面波方法.MASW方法的实施过程主要分为3个步骤：（1）基于线性接收排列采集高信噪比的面波数据，对面波数据进行数学变换提取频散能量；（2）结合频散能量极值与趋势，拾取频散极值得到频散曲线；（3）反演频散曲线，获得一维横波速度剖面，对多个一维横波速度剖面进行插值，形成拟二维横波速度剖面.这种MASW方法在面波多模式分离、地层结构探测、古城墙勘察、堆石坝密实度检测、地基加固效果检测、隧底隐伏岩溶探测中得到应用（Karray and Lefebvre，2009；田宝卿和丁志峰，2021；张瑜鹏等，2021；克里木等，2020；于东凯和魏亮，2020；王先龙，2020）.上述第3步骤对每条频散曲线采用半波长解释方法得到一维面波相速度剖面，然后对多个一维剖面进行插值，形成拟二维面波相速度剖面.这种方法在斜坡内部结构探测、岩溶通道探测、浅部煤层采空区探测中得到应用（刘现锋等，2020；李巧灵等，2022；杨智等，2022）. MASW方法实施过程的第二步提取频散能量是所用接收排列下方地层结构的综合响应，因此接收排列越长，对地层结构的横向平均效应越大，即横向分辨率越低.

在瑞雷面波法勘探中，理想的情况是尽可能利用少的地震道（MASW方法对地震道数和排列长度是有要求的）数据提取到较为可靠的频散曲线，这样既能够满足垂向探测的需求，在横向上也有较高的分辨率，达到能准确探测地下地质结构的目的.相邻道互相关函数法（也称相邻道相位差法）最早在被动源面波法中得到应用，21世纪开始应用在瞬态面波勘探中，该方法的关键是怎么把体波对提取面波频散曲线的干扰降到最低.徐华全（2009）直接从原始采样的地震波记录中采用相邻道互相关函数法提取频散曲线，并展示了该方法的应用效果.然而采用该方法直接提取的频散曲线容易受到体波信号的干扰.马见青等（2010）做了基于S变换的相邻道提取频散曲线方法，对实测资料做了试验处理.尹晓菲等（2018）提取基阶模式面波，用互相关方法提取相邻两道之间的面波相速度频散曲线，然后反演纯路径相速度频散曲线得到每个目标网格的一维横波速度剖面，进而得到拟二维速度剖面，在水平层状介质模拟及实测资料应用中取得较好的效果.陈杰等（2018）基于多重滤波方法直接对地震波记录相邻道提取频散曲线，并将其应用到理论模拟和实测资料处理中. 谌强等（2021）通过对原始地震数据在F-K域进行面波分离提取，然后对分离后的两道面波数据分别利用相位差法和多重滤波法计算出相速度后进行叠加，通过对两层水平及倾斜介质速度递增理论模型的波动方程模拟地震记录试算，说明了该方法的有效性，同时，其对混凝土球物理模型采集的数据进行了处理，在二维相速度剖面中可见高速异常体的存在，也说明了该方法的有效性，相比传统的MASW方法又前进了一步.

本文是基于相邻两道相位差法，通过在随机介质中建立低速异常体模型，进行波动方程数值模拟，采用线性拉冬（τ，p）变换提取面波，然后计算每道面波的相位谱，再计算每道初相位及相邻两道相位差，得到频散曲线，根据半波长理论得到拟二维面波相速度剖面，并将其应用到实际探测中.

1.   探测技术原理

当我们做瑞雷面波勘探时，野外采集的共炮点地震波记录里有纵波、横波和面波，当我们把面波作为有效波时，纵波和横波就是干扰波，它会严重影响面波频散曲线的计算.本文面波数据处理主要有5个步骤：①在地震波记录中分离提取面波；②利用提取的面波数据采用相邻两道相位差方法计算频散曲线；③根据半波长理论，将面波相速度随频率变化曲线转换为面波相速度随深度变化曲线；④将多个相邻两道的频散曲线集合起来，就得到面波相速度二维剖面；⑤剖面中相速度的大小反映地下介质的强弱.

从野外采集的地震波记录中分离提取面波，其核心算法就是线性拉冬变换，因为在时空（t，x）域面波和体波是叠加在一起的，当将地震波a（t，x）从时空域[t，x]变换到[τ，p]域后，面波和体波就基本上分布在不同的区域，然后将只有面波的区域再反变换到时空[t，x]域，这样就达到了分离提取面波的目的，下面给出具体的计算公式：定义p=dt/dx，p为视速度的倒数或称慢度，截距时间τ=t-px，t为地震波旅行时，x为偏移距.线性拉冬正变换公式（Thorson and Claerbout，1985）为：

假设a（t，x）有n道地震波记录，拉冬反变换公式为：

由于面波速度比体波小分布在一定的范围内，并且它主要沿着地表附近传播，传播方向与观测方向是一致的，所以其视速度与真速度相等，假设其在[τ，p]域对应的慢度区域为[p_k ，p_m]，截距时间区域为[τ₀，τ_j].在后面的数据处理中，就是根据地震波记录的实际情况设置适当的慢度区域和截距时间区域，用公式（1）和（2）进行正反两次变换提取面波的.

利用提取的面波记录ã（t，x），采用计算相邻两道相位差的方法得到n-1条频散曲线，以第k道与第k+1道两个相邻道为例，假设面波是从第k道向第k+1道传播，具体计算过程如下：

第k道的频谱的数值积分计算公式为：

m为采样点数，因为本文频率f是在对数域进行等间隔离散的，对低频进行了加密，所以在这里计算频谱的时候不能用快速傅里叶变换，直接用求和代替积分，其相位谱为：

注意，反正切函数的值域是[-π/2，π/2]，需要校正到[0，2π]，因为在复数域虚部是纵坐标，实部是横坐标，可以根据实部虚部的正负号将相位角校正到其对应的象限去.

同理，第k+1道的相位为φ_k+1（f），假设面波是从第1道向第n道传播的，则这两道的相位差为：

如果Δφ_k，_k+1（f）≤0，就加上2π，当道间距较大或频率较高时，相位差可能超过2π，就要补偿一个或多个周期.两道的时差为：

频散曲线为：

其中$ {\stackrel{-}{x}}_{k}=({x}_{k}+{x}_{k+1})/2 $为位置坐标，V_R为面波相速度.面波波长$ \lambda ={V}_{R}({\stackrel{-}{x}}_{k},f)/f $，根据半波长解释理论，探测深度h是波长的一半，即$ h=\lambda /2 $，这里就可以将频率f自变量换成深度h了，则面波相速度剖面在xoz面的表达式为：

2.   面波数值模拟

波动方程有限差分数值模拟是从理论上研究探讨解决实际问题的可行性，得到的结论可用于指导实际应用.这里着重介绍两种含低速体数学模型的建模，正演模拟，以及基于横向高分辨技术的数据处理结果分析.设计模型时考虑了城市工区场地的狭窄性，即不足以布置两个及以上排列.

本文是采用对二维弹性波动方程交错网格有限差分进行迭代计算的.数值模拟所使用的方程有两组，一是运动平衡方程，二是应力与应变的关系方程，将应变（位移）演算为质点振动速度，就可以建立应力与质点振动速度的关系式.以雷克（ricker）子波函数作为震源函数，将其加载到震源处的垂向质点振动速度上，先用现在时刻的质点振动速度和应力计算未来时刻的应力，再用现在时刻的应力和质点振动速度计算未来时刻的质点振动速度，这样就完成了一个时间间隔的迭代计算.模型的上边界（地表）采用自由边界条件，模型的左右及下边界采用完全匹配层（PML）条件.

这里建立数学模型的思路是先建立一个基于均匀介质背景的低速异常体模型，然后按照一定的随机性比例随机修改背景各个网格节点的物性参数，但保持低速异常体的物性参数不变，最后对模型参数进行3×3均滑处理.模型空间都设为23 m×25 m，网格间距dx=dz=0.5 m.在进行数值模拟计算时，时间步长0.000 1 s，记录长度0.8 s，震源子波的主频为25 Hz，以子波的周期40 ms作为延迟时，如图 1所示.

图  1  震源子波

Fig.  1.  The source wavelet

下载: 全尺寸图片 幻灯片

以下两个模型的初始数学模型背景介质的物性参数相同，异常体的物性参数也相同，只是异常体的形状大小不一样而已，背景介质的纵波速度Vp=2 000 m/s，横波速度Vs=800 m/s，密度ρ=1.9 g cm³，异常体的纵波速度Vp=1 200 m/s，横波速度Vs=300 m/s，密度ρ=1.2 g cm³.对背景介质各节点物性参数的随机修改是按照最大相对变化±20%生成随机介质模型的.

2.1   圆形低速异常模拟

设计圆形异常体的直径为2.5 m，其圆心在（11，5），初始模型如图 2a所示，受网格化影响圆形变为正方形了，在这个初始数学模型基础上建立的随机介质模型如图 2b所示，它是用3个物性参数中的横波速度绘制的.模拟计算了两炮地震波记录，震源点分别设在S₁（0，0）和S₂（23，0），道间距1 m，24道接收，地震波记录采样间隔0.5 ms，记录长度800 ms，为了显示主要波形，图像上只绘制了250 ms地震波记录，如图 3a、3c所示，由于排列短且面波能量强，波形记录中看不见体波直达波，低速异常体产生的绕射（散射）波通过波形振幅放大到一定程度依稀可见，主体波形是以背景介质产生的面波为主.

图  2  模型一的速度剖面

Fig.  2.  The velocity profile of model 1

下载: 全尺寸图片 幻灯片

图  3  模型一的地震波记录

a.左边炮原始记录；b.左边炮提取的面波记录；c.右边炮原始记录；d.右边炮提取的面波记录

Fig.  3.  Seismic wave recording of model 1

下载: 全尺寸图片 幻灯片

利用拉东变换公式（1）、公式（2）提取的面波如图 3b、3d所示，与原始波形对比剔除了绕射（散射）波，再利用公式（3）~（8）计算得到的瑞雷波相速度剖面如图 2c所示，与模型图 2b对比，反演的低速异常区域的主体与模型较吻合，低速区域的中心位置较准确，它的整体在水平方向有所拓展，没有其他低速区域干扰解译，虽然反演的相速度剖面与模型有许多的不同，但不影响对低速异常体的识别.图 4给出了两个正反向炮计算的23条叠加频散曲线中的3条，分别来自于8~9道（测点在8.5 m），11~12道（测点在11.5 m），14~15道（测点在14.5 m），对应图 4a、4b、4c，由于本文对频率是在对数域等间隔离散的，所以3个图的横坐标用频率的自然对数表示，对比3条曲线可以明显看出，中间速度低两边速度高，与低速异常体水平位置吻合.

图  4  模型一的三条频散曲线

Fig.  4.  Three dispersion curves of model 1

下载: 全尺寸图片 幻灯片

2.2   扁椭圆形低速异常模拟

考虑低速异常体的水平方向的贯通性，这里设计了一个长轴在水平方向的扁椭圆体型低速异常体，其宽度为12.5 m，高度为2.5 m，如图 5a所示，其物性参数与模型一的初始数学模型相同，受网格化影响其椭圆形态有所变化，由其为基础生成的随机介质模型如图 5b所示，模型中的低速区域为椭圆形.用该模型模拟计算了两炮地震波记录，震源点分别设在S₁（0，0）和S₂（23，0）处，道间距1 m，24道接收，地震波记录采样间隔0.5 ms，记录长度800 ms，为了显示主要波形，图像上只绘制了250 ms地震波记录，如图 6a，6c所示，由于排列短且面波能量强，波形记录中也看不见体波直达波，但低速异常体产生的绕射（散射）波清楚可见，主体波形还是以背景介质产生的面波为主.

图  5  模型二的速度剖面

Fig.  5.  The velocity profile of model 2

下载: 全尺寸图片 幻灯片

图  6  模型二的地震波记录

Fig.  6.  Seismic wave recording of model 2

下载: 全尺寸图片 幻灯片

利用公式（1）、（2）提取的面波如图 6b、6d所示，与原始波形对比剔除了绕射（散射）波，再利用公式（3）~（8）计算得到的瑞雷波相速度剖面如图 5c所示，与模型图 5b对比，反演的低速异常区域的主体与模型较吻合，低速区域的中心位置较准确.虽然模型图 5b与模型图 2b相比只是延伸了水平方向的低速区域，但反演的结果是低速区域在纵横向上都扩大了很多，中心低速异常区域成圆形状，与模型不符，这种现象给我们实际应用以很好的指导作用.

3.   实际应用

选择了两个具有代表性的工程实例，一个是在城市狭窄场地探测防空洞，洞顶埋深大约在8 m以上，电法、探地雷达等物探方法技术都不满足应用条件，多道瞬态面波分析法（MASW）也受到场地的限制，经济、高效、无损的方法只有面波横向高分辨技术.另一个是在宽阔场地探测地下渗水，可能渗漏区域较小，也是需要使用面波横向高分辨探测技术.

3.1   防空洞位置探测

该工区位于郑州市中心，场地狭窄只有23 m长，设计道间距1 m，24道接收，采样间隔0.2 ms，采样点数2 048，记录长度409.6 ms，检波器主频4.5 Hz，分别在第1道和第24道位置用落重锤激发了两炮地震波记录，如图 7a、7c所示，这两炮实测记录与前面的理论模拟记录相比较，频散严重波形复杂，说明地下介质的不均匀性更厉害，图 7a中隐约能看见绕射（散射）波的存在.提取的面波如图 7b、7d所示，与原始波形对比有较大的变化，去掉了绕射（散射）波，对面波处理得到的面波相速度剖面如图 7e所示，解释的防空洞顶埋深8.5 m，地面坐标X=18.5 m，防空洞的宽度和高度不好确定，钻探验证洞顶埋深8.1 m.

图  7  实例一地震波记录及其解译剖面

Fig.  7.  Example 1: seismic wave record and interpretation section

下载: 全尺寸图片 幻灯片

3.2   渗漏通道探测

该工区位于郑州市郊某处黄河大堤外堤脚，黄河洪水期在地面大约X=28 m为中心的区域有水渍，没有明显的水流，在此布置了3个排列，每个排列长度23 m，排列之间是连续衔接的，即第1个排列的第24道位置就是第2个排列的第1道，设计道间距1 m，24道接收，采样间隔0.2 ms，采样点数2 048，记录长度409.6 ms，检波器主频4.5 Hz，每个排列用落重锤激发两炮地震波记录，炮点分别在排列的两端外侧，偏移距都是8 m，地震波记录如图 8a~8c所示，提取的面波如图 8d~8f所示，与原始波形对比有较大的变化，去掉了干扰波，对面波处理得到的面波相速度剖面如图 9所示，图中可以看到明显的相对低速区域，建议在埋深h∈[0，16]，地面坐标X∈[22，33]的区域内进行封堵补强处理.

图  8  实例二地震波记录

Fig.  8.  Example 2: seismic wave record and interpretation section

下载: 全尺寸图片 幻灯片

图  9  瑞雷波相速度剖面

Fig.  9.  Rayleigh wave phase velocity profile

下载: 全尺寸图片 幻灯片

4.   结论

本文利用线性拉冬正反变换直接在[τ，p]域自动提取面波，采用相邻两道相位差法计算频散曲线，再用半波长解释方法将相速度与频率的关系转换为相速度与深度的关系，此时得到了相速度随测点及深度的变化关系，即拟二维面波相速度剖面，应用这项面波横向高分辨技术对两种模型数据及两个实测数据进行处理，得到如下的结论：

（1）理论模型试算表明，对一个排列两炮地震波记录，采用面波横向高分辨技术可以较准确地探测到低速异常体的中心位置，但成像结果中的低速异常区域与模型有一些的出入，对于圆形低速异常体横向被扩大的很明显，对于水平长条状低速异常体纵横向上都被扩大了，这对于实际数据的成像应用解释具有很好的指导意义.

（2）该技术用于两个典型工况的实际探测数据处理，经验证实取得了很好的效果.

（3）该技术适应于城市狭窄空间条件下的地下地质异常体探测，克服了传统的MASW方法的横向平均效应缺陷.