主–被动源面波联合反演的浅地表结构成像方法研究

Study on Imaging Method of Shallow Surface Structure Based on Joint Inversion of Active and Passive Source Surface Waves

摘要: 近年来，随着城市发展，人们开发城市地下空间的需求越来越大，城市地下灾害隐患也越来越多。如何有效、高精度和高分辨率地获取城市地下空间结构受到了更多的关注。本文利用主被动源面波联合成像方法对城市浅地表地下空间结构进行了探测。通过主被动源面波频散曲线的拼接和自动拾取，对某地下人行通道进行城市浅地表成像研究，清晰显示出埋深9~20 m的低速异常体，与实际资料中异常体位置高度吻合，具有较高的分辨率，为主被动源联合的面波探测方法在城市浅地表成像中的推广应用提供方法支持。

Abstract: In recent years, with the development of cities, people's demand for developing urban under-ground space is increasing, and the hidden dangers of urban underground disasters are also in-creasing. How to obtain urban underground spatial structure effectively, with high precision and high resolution has received more attention. In this paper, the combined imaging method of active and passive surface waves is used to detect the underground space structure of urban shallow surface. Through the splintering and automatic picking of surface wave dispersion curves of active and passive sources, the imaging research on urban shallow surface of an underground pedestrian passage is carried out, and the low-speed abnormal body with a buried depth of 9~20 m is clearly shown, which is highly consistent with the location of the abnormal body in the actual data, and has a high resolution. This method supports the promotion and application of the surface wave detection method combined with active and passive sources in urban shallow surface imaging.

1. 引言

随着近年来人们对城市地下空间的开发与利用，地质活动、开发工程、管线泄漏等造成了城市地下灾害隐患，需要对其进行探查和预防；同时在对城市地下空间开发前也需要对其结构进行探查。这些都使城市地下空间成像得到越来越多学者的关注和研究，目前城市地下空间探测的常用地球物理方法主要有浅层地震、面波勘探、高密度电法、探地雷达等，各个方法都在不同的应用条件下有着较好的探测效果，但面波勘探方法以其抗干扰强、探测深度大、应用场景多、布设简单等优势得到了更广泛的应用。

为充分利用面波高阶模式波在横波速度成像上的优势，陈晓非院士团队近五年来在高阶模式面波理论研究、频散能谱生成及应用上开展了一系列的研究，并应用与城市浅地表探测，取得了丰硕的成果 [1] [2] [3] 。夏江海团队近年来对背景噪声成像在城市地下空间探测的应用也展开了多方面研究，多次应用于实际，并取得颇多成果 [4] [5] [6] [7] 。杨抗等采用线性观测方式的主被动源面波联合方法对城市地下空洞进行探测，成像结果与钻孔揭露结果较为吻合 [8] 。杨志等对相移法进行改进，引入低频聚焦因子，其具有在短接收排列上提取低频频散能量的优势。该方法显著改善了频散能量在低频端的聚焦性，扩展了可拾取的频带范围，增大了探测深度，缩短了计算频散能量所需的排列长度，提高了面波多道分析方法的横向分辨率 [9] 。秦建军等在成都地区通过模型试验提出了研究地区城市地下灾害探测的最佳观测系统，并标定了对地下疏松体频散图谱的判识标准 [10] 。Shao J等利用光纤分布式代声学传感(DAS)记录城市交通噪声数据，通过反演面波频散曲线获取了测线下方的近地表横波速度剖面，证明了DAS在城市交通噪声近地表成像的有效性 [11] 。代新雲等利用高密度电法和连续面波法相结合对达州市城区某道路进行地下病害体探测，对其加以钻探验证，能基本确定地下病害体的位置及规模情况，为城市地下病害体的治理提供理论依据 [12] 。本文采用主被动源面波联合成像方法，实现了主被动源面波联合频率自动拾取，对某地下人行通道区域进行城市浅地表成像，清晰显示出埋深9~20 m的低速异常体，与实际资料中异常体位置高度吻合，误差小于1 m，对2 m尺寸的异常体分辨良好，具有着较高的分辨率。

2. 研究方法

主被动源面波联合浅层成像流程如图1所示，具体步骤如下：

1) 根据主被动源面波联合浅层成像原理和研究区域实际情况，进行观测系统设计和优化。前往研究区完成观测系统布设和实际数据采集。

2) 数据采集完成后，对原始数据进行筛选，降噪，归一化等预处理，分析数据质量。

3) 采用F-K法拾取主动源面波频散曲线，对主动源面波数据进行F-K变换，得到主动源F-K谱，并根据波数K与相速度V转换关系得到主动源面波F-V谱。对F-V谱分频率归一化，避免某一频率范围内能量过大导致其他频率范围F-V谱分辨率过低的情况。

4) 为实现自动拾取主被动源面波联合频率，对各频率下不同面波相速度对应能量幅值求取均值，观察分析不同能量收敛程度下均值的大小。设定阈值，若能量均值超过阈值则认为该频率下能量逸散，无法拾取频散曲线，将对应频率下F-V谱去除，实现主被动源面波联合频率求取和主动源F-V谱剪切的自动化。对剪切后主动源F-V谱各频率下能量最强点进行拾取，得到主动源频散曲线。

5) 采用空间自相关法拾取被动源面波频散曲线，以联合频率为上限，求取被动源面波自相关系数，拟合零阶贝塞尔函数，拾取被动源频散曲线。

6) 联合主被动源面波频散曲线，根据频率、面波相速度、波长关系由频散曲线求得不同速度面波波长，以半波长作为探测深度，可求得H-V曲线。

7) 为得到高质量初始模型，将目标深度内地层均匀划分为若干层，以各层范围内H-V曲线平均速度作为该层速度，得到后续反演所需初始模型。

8) 对频散曲线进行最小二乘反演，将反演结果与实际资料进行对比，计算误差，以验证该方法的精度和分辨率。

Figure 1. Combined imaging flow chart of active and passive source surface waves

图1. 主被动源面波联合成像流程图

3. 观测系统设计

本文采用多道瞬态面波分析法和空间自相关法的联合观测系统，即线性观测阵与同心圆观测阵嵌套观测系统，以达到主被动源面波数据同时采集的目的，并前往研究区进行实际调查，结合研究区实际情况和探测目的深度设定观测系统的炮检距，道间距，同心圆半径等参数。

研究区俯视图如图2所示，研究区紧邻主干道辅路，受场地限制，只可在人行道和绿化带范围进行检波器布设，经过实地测量，该研究区允许布设线性测线最大长度34 m，圆形观测阵列最大半径5 m。研究区地下人行道施工图纸如图3所示，红圈位置为观测点，测线走向与人行道走向垂直，图中可见观测点下方约9~11 m存在地下管线结构，约12~20 m为地下人行道，两者间距不足1 m。因此观测系统目的深度应大于20 m，故设定探测目表深度为25 m。

Figure 2. Top view of the study area

图2. 研究区俯视图

Figure 3. Construction drawing of underground pedestrian passage in research area

图3. 研究区地下人行通道施工图纸

根据经验公式，空间自相关法的同心圆观测阵探测深度约为同心圆半径的3~5倍 [13] ，场地限制最大半径为5 m，满足探测目标深度要求，故设定半径为5 m；主动源测线在场地限制内设置最长长度34 m，观测点18个，道间距2 m，以避免观测点数过低降低探测精度 [14] ，设定炮点5个，炮间距为4 m，最小炮检距4 m，最大炮检距20 m。最终设计观测系统如图4所示，观测系统实际布设和数据采集现场如图5所示。

Figure 4. Active and passive source surface wave joint observation system

图4. 主被动源面波联合观测系统

Figure 5. Observation system layout and data acquisition site

图5. 观测系统布设及数据采集现场

4. 数据处理

4.1. 主动源面波数据处理

对主动源面波原始数据进行筛选，降噪，归一化等预处理，分析数据质量，本文采用频率–波数法求取主动源面波F-V谱，该方法通过二维傅里叶变换将主动源面波数据由时间–空间域 $g\left(t,x\right)$ 转换至频率–波数域 $G\left(\omega ,k\right)$ ，转换关系如公式(1)所示 [15] ：

$G\left(\omega ,k\right)={\displaystyle {\int }_{-\infty }^{+\infty }{\displaystyle {\int }_{-\infty }^{+\infty }g\left(t,x\right){\text{e}}^{-j\left(\omega t+kx\right)}\text{d}t\text{d}x}}$ (1)

上式得到主动源面波的频率–波数谱，利用波数与相速度转换关系继续计算，如公式(2)所示 [15] ：

$k=\frac{1}{\lambda }=\frac{f}{V}$ (2)

其中k为波数，f为频率，v为相速度。从而得到主动源面波数据的F-V谱，并对其进行分频率归一化，避免某一频率范围内能量过大导致其他频率范围内能量分辨率过低的情况。

为实现主被动源面波频散曲线联合频率的自动求取，本文对各频率下不同面波相速度对应的能量大小求取均值，对比分析不同能量收敛程度下的均值大小，设定分辨能量是否足够收敛来拾取频散曲线的阈值。若均值超过阈值则认为该频率下能量逸散，无法拾取频散曲线，予以去除，以达到自动求取主被动源频散曲线拼接频率以及裁剪主动源F-V谱的目的。由于该方法建立在将各频率下能量归一化的基础上，理论上其具有良好的泛用性。

将该方法应用于上文求得的主动源面波F-V谱，自动判定17 Hz为主被动源面波频散曲线联合频率，将低于17 Hz的主动源面波F-V谱剪切去除(如图6所示)。

Figure 6. Active source F-V spectrum shear results

图6. 主动源F-V谱剪切结果

对剪切后主动源F-V谱各频率下能量最强点进行拾取，如图7所示，得到主动源频散曲线。

Figure 7. Active source dispersion curve

图7. 主动源频散曲线

4.2. 被动源面波数据处理

被动源面波数据采用空间自相关法进行处理，该方法由日本著名地震学家Aki在1957年提出 [16] ，他的理论奠定了该方法在未来的发展，其中最为核心的有两部分，分别为：

1) 随机振动在时空上具有稳定性；

2) 陆地的表面波动会进行频散性的传播。Aki通过计算推导得出了在一维和二维条件下的面波频散曲线的计算公式 [15] 。

一维条件下，假设角频率为 ${\omega }_0$ ，则在此方位的相速度 $c\left({\omega }_0\right)$ 与空间自相关系数满足公式(3)。

$\rho \left(\delta ,{\omega }_0\right)=\cos \left(\frac{{\omega }_0\delta }{c\left({\omega }_0\right)}\right)$ (3)

式中： $\delta$ 表示两检波器之间的距离。

公式(3)成立的条件是震源严格位于同一直线，但是在实际工作中无法得到震源的位置分布情况，因此Aki研究得出了在二维条件下计算频散曲线的方法，该方法具有更好的实用性，计算公式如公式(4)所示：

$\bar{\rho }\left(r,{\omega }_0\right)=J_0\left(\frac{{\omega }_0}{c}r\right)$ (4)

以上公式是在震源分布方式为随机分布的前提下求得的，其中：J₀表示第一类零阶贝塞尔函数， $\bar{\rho }$ 表示方位平均之后的空间自相关系数，r代表两个检波点之间的距离，c表示瑞雷面波的相速度， ${\omega }_0$ 代表角频率。但在实际采集工作当中，随机震源的分布方式往往为不均匀分布，为保证采集微动信号记录的质量，学者们进行了许多实验和研究，寻找减弱甚至是消除随机震源不均匀分布的影响，最终得到了将观测系统设计为圆形台阵的方法。

空间自相关方法的观测系统设计一般为圆形台阵，检波器数量至少为四个，其中圆心布设一个检波器，其余三个均匀分布在以某一特定半径的圆上，形成等边三角形的形状，半径的取值取决于目标探测深度，计算方法前文已有介绍。本研究需要快速傅里叶变换求功率谱和互功率谱，具体如下方公式所示，对每个台站微动信号记录进行此步处理，计算过程如公式(5)~(7)所示 [2] ：

$S\left(r,\theta ,\omega \right)={\displaystyle {\int }_{-\infty }^{+\infty }{R}_{1,2}\left(\tau \right){\text{e}}^{-i\omega t}\text{d}t}$ (5)

$S\left(0,\omega \right)={\displaystyle {\int }_{-\infty }^{+\infty }{R}_{1,1}\left(\tau \right){\text{e}}^{-i\omega t}\text{d}t}$ (6)

$S\left(r,\omega \right)={\displaystyle {\int }_{-\infty }^{+\infty }{R}_{2,2}\left(\tau \right){\text{e}}^{-i\omega t}\text{d}t}$ (7)

式中： $R_{1,2}={\displaystyle \sum {\mu }_0\left(t\right){\mu }_r\left(t+\tau \right)}$ ， $R_{1,1}={\displaystyle \sum {\mu }_0\left(t\right){\mu }_0\left(t+\tau \right)}$ ， $R_{2,2}={\displaystyle \sum {\mu }_r\left(t\right){\mu }_r\left(t+\tau \right)}$ ， $S\left(r,\theta ,\omega \right)$ 为圆心台站与半径为r的圆上某一台站的互功率谱，其中 ${\mu }_0$ 为圆心点处台站的微震记录， ${\mu }_r$ 为半径为r的圆上所布设台站的微震记录， $S\left(0,\omega \right)$ 为圆心点台站记录的自功率谱， $S\left(r,\omega \right)$ 为半径r的圆上一台站记录的自功率谱。空间自相关函可以用得到的功率谱表示如公式(8)所示 [17] ：

$\rho \left(r,\omega \right)=\frac{1}{\text{2π}}{\displaystyle {\int }_0^{\text{2π}}\frac{\mathrm{Re}\left[S\left(r,\theta ,\omega \right)\right]}{\sqrt{S\left(0,\omega \right)S\left(r,\omega \right)}}\text{d}\theta }$ (8)

下一步便是将所有位于同一圆形上的检波器采集数据求得的自相关系数进行方位平均，获得该圆形台阵微动信号记录的自相关系数。

根据Aki的理论假设，微动台阵记录的空间自相关系数可以用第一类零阶贝塞尔函数表示，如公式(9)所示：

$\rho \left(r,f\right)=J_0\left(\frac{\text{2π}f}{c\left(f\right)}r\right)$ (9)

式中， $c\left(f\right)$ 是频率f对应的瑞雷波相速度，r代表台站之间的距离，J₀为第一类零阶贝塞尔函数，f代表频率，ρ为方位平均后的自相关系数。

Figure 8. Autocorrelation coefficient

图8. 自相关系数

以上文得到的主被动源面波频散曲线合并频率17 Hz为上限，求取被动源面波数据自相关系数并做平滑处理，如图8所示。将自相关系数与第一类零阶贝塞尔函数进行拟合，根据拟合残差极小点拾取被动源频散曲线，如图9所示，得到被动源频散曲线。

Figure 9. Passive source surface wave dispersion curve

图9. 被动源面波频散曲线

4.3. 主被动源联合反演

将上文得到的主被动源频散曲线进行联合，得到0~50 Hz范围的主被动源面波联合频散曲线，如图10所示。

Figure 10. Combined dispersion curve of active and passive source surface waves

图10. 主被动源面波联合频散曲线

为得到后续反演所需的初始模型，还需进一步建立相速度与穿透深度的关系。波长、频率和相速度存在如公式(10)所示关系 [18] ：

${\lambda }_R=V_R/f$ (10)

式中 ${\lambda }_R$ 为面波波长，V_R为面波相速度，f为频率。

一般认为面波的穿透深度为其半波长 [19] [20] ，即公式(11)所示关系：

$H={\lambda }_R/2$ (11)

按公式(9)和公式(10)计算可由频散曲线得到H-V曲线，结算结果如图11所示。

Figure 11. H-V curve

图11. H-V曲线

Figure 12. Initial model

图12. 初始模型

将深度0~30 m的地层均匀分为30层，每层1 m，按层对H-V曲线进行平均，得到后续反演所需初始模型，如图12所示。

本文采用最小二乘法进行反演，将上文得到的频散曲线和初始模型输入，如图13所示，经过31次迭代，拟合误差达到最小值且趋于稳定，为2.94%，误差较小;频散曲线拟合结果如图14，可以看出拟合结果与实际频散曲线较吻合。

Figure 13. Iterative error variation

图13. 迭代误差变化

Figure 14. Dispersion curve fitting results

图14. 频散曲线拟合结果

最后如图15所示，将主被动源面波联合成像得到的横波速度结构与前文展示的地下人行道施工图纸进行比对，图中左侧为反演结果，右侧为施工图纸，褐色虚线为对比辅助线。可以看出横波速度结构所显示的低速异常体顶板位于地下9 m处，底板位于地下20 m处，在12 m深度处出现横波速度跳跃式下降，故可分为两部分，以12 m为分界点；图纸显示观测点下方约9~11 m存在地下管线结构，约12~20 m为地下人行道，两者间距不足1 m。由于管线结构内存在内容物，并非空洞，而地下人行通道可视为空洞，理论上管线结构对应横波速度应大于地下人行通道对应横波速度，该情况与成像结果高度吻合，且各对应位置误差小于1 m，且对埋深9~11 m的2 m尺寸地下管线异常体分辨良好，可证明主被动源面波联合成像方法具有着较高的垂向分辨率。

Figure 15. Comparison between the surface wave imaging results of active and passive sources and the actual data

图15. 主被动源面波浅层成像结果与实际资料对比图

5. 结论

本文将主被动源面波联合成像方法实际应用于某地下人行道区域，开展浅地表成像研究，得到以下结论：

1) 实现了主被动源面波频散曲线联合频率的自动拾取，为主被动源面波联合反演的自动化提供了方法支持。

2) 将成像结果与实际地下人行通道施工图纸进行比对，发现主被动源面波联合反演的浅地表结构成像方法所得横波速度结构与实际地下结构高度吻合，误差小于1 m，对2 m尺寸的地下异常体分辨良好，证明了该方法具有较高的垂向分辨率，在城市浅地表探测领域有着较高的应用价值，为主被动联合的面波探测方法在城市浅地表成像的推广应用提供了方法支持。

致谢

感谢北京市勘察设计研究院有限公司提供实际资料。

基金项目

国家重点研发计划(项目编号：2021YFC3090302)；

北京市自然科学基金资助项目(项目编号：8232012)；

中央高校基本科研业务费资助(项目编号：2021JCCXDC02)；

煤炭精细勘探与智能开发全国重点实验室开放基金资助(项目编号：SKLCRSM21LH01)。