城镇地区短期同步-异步背景噪声层析成像—

Engineering ›› 2025, Vol. 48 ›› Issue (5) : 309 -325. DOI: 10.1016/j.eng.2025.02.001

研究论文

城镇地区短期同步-异步背景噪声层析成像——岩溶调查中的应用

柳亚 ^a ,
夏江海 ^a ,
管博 ^a ,
席超强 ^b ,
宁玲 ^c ,
张浩 ^d

作者信息 +

Short-Term Synchronous and Asynchronous Ambient Noise Tomography in Urban Areas: Application to Karst Investigation

Ya Liu ^a ,
Jianghai Xia ^a^,^* ,
Bo Guan ^a ,
Chaoqiang Xi ^b ,
Ling Ning ^c ,
Hao Zhang ^d ,

Author information +

文章历史 +

Received	Accepted	Published
2023-08-16	2025-02-08
Issue Date
2025-04-22

PDF (10456K)

摘要

密集台阵背景噪声层析成像是实现高分辨率地下成像的强大工具，对地质灾害防控具有重要影响。然而，传统密集台阵研究需要大量台站同步观测以实现广泛覆盖。为了使用有限台站进行全面的岩溶特征调查，我们设计了一种新型同步-异步观测系统，以促进密集台阵观测。我们结合同步骨干台站，开展了两轮各持续约24 h的异步观测。我们利用197个节点式检波器实现了研究区域的广泛覆盖，平均台站间距为7 m。聚束分析结果显示了昼夜间噪声源分布的明显变化。我们估计了稳相区的震源强度，并使用加权方案对互相关函数（ $C 1$ 函数）进行叠加，以抑制非均匀噪声源分布的影响。叠加权重由与瑞利波相关的多分量 $C 1$ 函数之间的相似性系数导出。我们采用 $C 1$ 函数的互相关（ $C 2$ 方法）来获取异步台站之间的经验格林函数。为了消除由 $C 1$ 函数高阶模式面波引起的 $C 2$ 函数伪影，我们基于与多阶瑞利波相关的不同粒子运动对 $C 1$ 函数进行了滤波。从 $C 1$ 和 $C 2$ 函数获得的瑞利波频散测量被用于面波层析成像。反演得到的三维剪切波（S波）速度模型揭示了在40~60 m深度范围内存在两个显著的低速区，这与钻井数据中发现的岩溶洞穴吻合良好。短期同步-异步背景噪声层析成像方法为城镇地质灾害调查提供了一种经济高效的新途径。

Abstract

Dense-array ambient noise tomography is a powerful tool for achieving high-resolution subsurface imaging, significantly impacting geohazard prevention and control. Conventional dense-array studies, however, require simultaneous observations of numerous stations for extensive coverage. To conduct a comprehensive karst feature investigation with limited stations, we designed a new synchronous–asynchronous observation system that facilitates dense array observations. We conducted two rounds of asynchronous observations, each lasting approximately 24 h, in combination with synchronous backbone stations. We achieved wide-ranging coverage of the study area utilizing 197 nodal receivers, with an average station spacing of 7 m. The beamforming results revealed distinct variations in the noise source distributions between day and night. We estimated the source strength in the stationary phase zone and used a weighting scheme for stacking the cross-correlation functions ( $C^{1}$ functions) to suppress the influence of nonuniform noise source distributions. The weights were derived from the similarity coefficients between multicomponent $C^{1}$ functions related to Rayleigh waves. We employed the cross-correlation of $C^{1}$ functions ( $C^{2}$ methods) to obtain the empirical Green’s functions between asynchronous stations. To eliminate artifacts in $C^{2}$ functions from higher-mode surface waves in $C^{1}$ functions, we filtered the $C^{1}$ functions on the basis of different particle motions linked to multimode Rayleigh waves. The dispersion measurements of Rayleigh waves obtained from both the $C^{1}$ and $C^{2}$ functions were utilized in surface wave tomography. The inverted three-dimensional (3D) shear-wave (S-wave) velocity model reveals two significant low-velocity zones at depths ranging from 40 to 60 m, which align well with the karst caves found in the drilling data. The method of short-term synchronous–asynchronous ambient noise tomography shows promise as a cost-effective and efficient approach for urban geohazard investigations.

关键词

地震干涉法 / 面波层析成像 / 异步背景噪声 / 地质灾害 / 地震背景噪声

Key words

Seismic interferometry / Surface wave tomography / Asynchronous ambient noise / Geohazards / Seismic ambient noise

引用本文

引用格式 ▾

柳亚,夏江海,管博,席超强,宁玲,张浩. 城镇地区短期同步-异步背景噪声层析成像——岩溶调查中的应用[J]. 工程（英文）, 2025, 48(5): 309-325 DOI:10.1016/j.eng.2025.02.001

登录浏览全文

4963

注册一个新账户忘记密码

1 引言

21世纪初，地震干涉理论揭示，对背景噪声进行互相关计算（即C¹方法）可得到地下结构的经验格林函数（EGF）[1‒3]。EGF是格林函数的近似值，格林函数描述了脉冲源在某一点产生的信号在另一点的接收过程，反映了地球地下介质的脉冲响应。面波能量在背景噪声中占主导，研究人员可从背景噪声数据中提取面波[4‒6]。背景噪声面波研究涵盖多尺度和多应用场景，包括精细地壳与岩石圈结构刻画[7‒13]、全球尺度层析成像[14]、浅地表结构成像[15‒18]、火山结构刻画[19]、地热储层成像[20‒22]以及断裂构造勘察[23‒25]。在浅地表应用中，研究表明，通过交通诱发地面震动的较短时长（数小时至数分钟）背景噪声记录，可获得高质量面波成像结果[18,26‒31]，大幅缩短了传统背景噪声观测时长。背景噪声成像具有环境友好、成本效益高、抗噪声干扰能力强等优势，在人类活动复杂的人口密集城区具有广阔应用前景。

岩溶灾害是常见的地质灾害，常导致地面沉降、塌陷、地下水污染[32]，严重影响人类居住安全与稳定，造成重大经济损失和不良社会影响[33]。明确岩溶特征的发展规律与分布特征，对维护社会安全稳定、促进经济可持续发展至关重要。目前已有多种地球物理方法用于异常区域刻画，包括声波走时层析成像[34‒36]、折射波法[37‒39]、电阻率法[40‒41]、电磁法[42‒43]、地震波形反演[44]以及多道面波分析[45‒46]。但复杂的城市环境（如土地可达性受限、噪声干扰等）给地球物理成像方法带来了巨大挑战。本文采用背景噪声成像技术对地下岩溶特征进行刻画。

密集台阵观测与背景噪声面波层析成像相结合，可提供高分辨率的地下介质3D横波（S波）速度结构[47‒49]。但密集台阵需大量台站，成本较高且有时难以实现。近年来，地震干涉理论提出了通过对C¹函数进行二次互相关或卷积来得到EGFs [50‒52]。对C¹函数尾波进行互相关称为C³方法[50,53‒54]，而C²方法则对C¹函数直达面波进行互相关或卷积[55‒57]。C¹函数的二次互相关或卷积使异步观测利用成为可能[58]。

本文在某废弃停车场设计了新型同步-异步背景噪声观测系统（图1），利用有限台站（197台节点接收器）实现了短期（两轮观测，每轮约24 h）密集台阵观测（平均台站间距为7 m），以实现地下岩溶特征的高分辨率探测与成像。首先，分析了城市背景噪声复杂的时空特征，提出了新的C¹函数加权叠加方法，以抑制噪声源分布不均对C¹函数精度的影响。其次，通过利用基阶与高阶瑞利波质点运动的差异，可抑制C¹函数中的高阶面波，从而提高异步C²的精度。结合C¹和C²函数获取的瑞利波频散数据，研究采用背景噪声面波层析成像方法，得到了高分辨率的地下3D S波速度模型。40 m以下识别出的低速区与钻探数据指示的岩溶洞穴吻合良好。最后，本文分析了模型分辨率及异步观测带来的改进效果。

2 数据

2.1 同步与异步地震台阵设计

研究区域为杭州某废弃停车场，该场地上方规划建设一栋办公楼[图1（a）]。钻孔数据显示，地下存在大量岩溶洞穴[图1（b）]，对拟建建筑安全构成威胁。采用背景噪声面波层析成像技术可获取地下异常的空间分布，以应对岩溶威胁。值得注意的是，停车场区域内共钻探了36个岩性钻孔[图1（b）]，其中，三个钻孔还进行了横波速度测井。这些信息为反演及模型解释提供了辅助验证依据。

为充分利用有限的地震仪（本研究共197台节点检波器）实现研究区域的大范围与高密度覆盖，本文设计了包含同步与异步地震台阵的新型地震观测系统[图1（a）]。在停车场范围内开展两轮观测，每轮观测时长约24 h。每轮观测包含6条测线，测线间距约10 m，每条测线上台站间距约5 m。第一轮观测结束后，立即移动台站开展第二轮观测。同时，在停车场外额外部署29台台站，在整个观测期间持续记录背景噪声。这些额外台站在后续分析中起桥梁作用，用于重建两轮异步观测台站组之间的面波。

2.2 背景噪声特征

本研究采用197台三分量节点检波器（固有频率为5 Hz），采样率为500 Hz。在两轮观测中，各选取一个节点[图1（a）中用黑圈标出]进行时频分析。如图2所示，由于观测在连续两个工作日进行，其时频特征高度相似。背景噪声频带范围较宽，白天能量较高，夜间能量较低，尤其是在20 Hz以上频率时，夜间能量降低更为显著。夜间时频谱中观察到多条垂直白色谱线，推测其与夜间零星交通事件有关。此外，午间时段能量明显降低，推测与人类活动（尤其是午餐和午休期间交通活动减少）有关。时频分析结果表明，记录的背景噪声主要来源于人类活动，且具有明显的昼夜变化特征，这也反映了人工震源复杂的时变特性。

除时间和频率特征外，本文还分析了震源的空间变化特征。背景噪声源分布影响地震波提取[3,59‒62]。本研究利用停车场内台阵第一轮数据，采用波束形成分析研究震源分布。利用台阵数据，波束形成可扫描穿过台阵的地震波方位角和速度。波束形成会叠加台阵信号，叠加能量越强，表明地震波的对应方位角越真实。不同频率的背景噪声源可能来自不同方向。本文分析了3 Hz和10 Hz频率的噪声源分布，并基于连续短时窗得到了震源分布的时间变化（图3）。

3 Hz噪声源在整个观测期间的分布呈现出近似随机均匀的模式[图3（a）]，这可能是由于低频地震波衰减小，推测其可能来自远处震源。相比之下，10 Hz噪声源能量具有显著的方向性和昼夜变化[图3（b）]。白天，能量主要集中在35°～120°和200°～300°方向，对应场地附近的交叉路口位置。根据现场观测，这些方向除交通外无其他明显人工活动。场地周围为办公楼，观测到的能量集中主要是因为这些交叉路口的交通诱发震动。夜间能量分布更均匀，无明显方向性，表明10 Hz噪声源主要来自附近道路交通。值得注意的是，为更好地观察源分布特征，对每个时段的波束形成分析结果进行了归一化，因此不同时窗间的震源能量不具直接可比性。

基于记录的背景噪声时空频分析，研究发现城市环境中的背景噪声具有高度复杂性。这种复杂性可能对地震干涉法的基本假设（尤其是损害噪声源随机均匀分布的假设）构成挑战，这关乎能否可靠提取地震波信息。因此，在数据处理中必须格外谨慎。下一节将介绍从同步和异步记录中提取面波及进行面波频散测量的方法。

3 方法

3.1 C¹函数计算

稳相理论表明，稳相区（SPZs）内的噪声源对从背景噪声中提取地震波起关键作用。在SPZs内，噪声源产生相长干涉，而区外噪声源产生相消干涉。增强SPZs内的噪声源可显著提高提取面波的精度[63‒64]。本文数据的上述分析表明噪声源分布不均匀，所以继续采用传统处理方法并依赖长时间叠加获取可靠信号可能并不合适。因此，本文重点保留仅包含SPZ源的时段，以提高提取面波的精度。

面波的SPZ包括两台站连线方向及其周围区域。与以往采用台阵方法分析沿线方向源强的思路不同，本文利用交叉分量C¹函数的相似性分析沿线源强。

在理想远场平面波和噪声源均匀分布情况下，瑞利波的交叉分量C¹函数满足以下关系[65]：

ϕ Z Z = ϕ R R

(1)

ϕ Z R - ϕ R Z = π

(2)

式中，

ϕ

(:)表示频域中的相位项。Z和R分别代表垂直分量与径向分量。ZR表示一个接收器的Z分量与另一个接收器的R分量之间的C¹函数。该约定同样适用于ZZ、RR和RZ，其中，ZZ和RR表示两个接收器相同分量之间的C¹函数，RZ则表示一个接收器的R分量与另一个接收器的Z分量之间的C¹函数。需特别强调的一点是：式（1）、式（2）以及后文提及的所有C¹函数，均为C¹函数的对称分量（因果部分与时间翻转后的非因果部分的平均值），此举旨在降低潜在噪声源分布不均匀带来的影响[66]。瑞利波呈椭圆极化特征，其垂直分量与径向分量之间存在相位差。根据质点运动类型（逆向运动或顺向运动）的不同，瑞利波垂直分量与径向分量的相位差可为

π / 2

或-

π / 2

[67]。逆向质点运动与顺向质点运动描述的是地震波中质点的运动方式：顺向运动指质点运动方向与面波传播方向一致，逆向运动则指质点运动方向与面波传播方向相反。当噪声源呈均匀分布或沿两个接收器的连线方向分布时，式（1）、式（2）适用于不同类型的质点运动。因此，本文提出将ZZ与RR的相似性、ZR与-RZ的相似性作为评估沿线震源强度的指标：

W = m a x P Z Z, R R, P Z R, - R Z, i f W < δ, W = 0

(3)

式中，P(..., ...)表示计算皮尔逊相关系数；W为两组相关系数（ZZ与RR、ZR与RZ）中的最大值；δ为阈值。若W低于阈值，则将其设为0。计算皮尔逊相关系数时需去除数据均值并按标准差（SD）归一化，使其成为评估信号相似性的合适指标。权重具有数据自适应性，仅需设置阈值。W越接近1，表明沿线震源越强；反之，W越接近0，表明震源更可能位于线外位置。研究可获取每个背景噪声时段的W，并将其作为后续C¹函数叠加的权重。这样可选择沿线震源强的数据段进行C¹函数叠加，以抑制噪声源分布不均匀对面波提取的影响，提高信号的可靠性和质量。

为验证加权叠加方法的有效性，本文给出了噪声源非均匀分布的模拟数据案例。采用双层模型作为地下介质（表1），模拟1 h的三分量（Z、N、E）背景噪声数据[68]，其中，N和E分别为南北分量和东西分量。1400个震源分布在内、外半径为500～1500 m的环形区域内（图4）。具体而言，研究特意在120°～130°和320°～330°方位角范围内增加了额外噪声源，以形成噪声源非均匀分布。中心区域布设25个接收器组成线性台阵，台站间距为5 m。所有点源的强度和方位随机生成，使得背景噪声中同时存在瑞利波和勒夫波。重要的是，仅模拟基阶面波，以确保波形和频散能量中观察到的任何额外信号均由噪声源非均匀分布引起，而非高阶面波。

研究对数据进行了多项标准预处理[66,69]，包括数据分段（3 s长，50%重叠）、去均值、去趋势和谱白化。随后以最左侧台站为参考台站，计算其与其他台站间的C¹函数。由于处理的是三分量数据，得到了ZN、ZE、NE等交叉分量C¹函数（ZN表示N分量与Z分量间的互相关，ZE表示E分量与Z分量间的互相关，NE表示N分量与E分量间的互相关）。通过分量旋转得到包含瑞利波（垂直分量和径向分量）的C¹函数，即ZZ、RR、ZR和RZ。计算数据段的权重W，并对ZZ分量C¹函数进行加权叠加。图4（a）显示了震源和接收器的分布，为数据模拟提供了空间背景。图4（b）对比了线性叠加和加权叠加得到的虚拟炮集（VSGs）。VSG根据台站与参考台站的距离对互相关对进行排序形成。与新叠加方法的结果相比，线性叠加的VSG存在虚假早到波。对VSGs进行多道频散成像[70]，结果表明原始叠加方法在频散谱中产生大量伪影[图4（c）]，而新叠加方法的频散趋势与理论频散曲线吻合良好，且无噪声干扰[图4（d）]。应用新叠加方法后，噪声源非均匀分布对面波提取的影响已被消除。

对于实际野外数据，采用新方法获取C¹函数。首先将数据下采样至100 Hz以提高计算效率，然后进行1～25 Hz带通滤波。数据预处理流程与模拟测试相同，但分段窗口设为30 s。采用阈值δ = 0.4计算叠加权重W，所有不同分量的C¹函数均按所得权重进行加权叠加。图5显示了以第一轮数据首道为参考道的VSG。传统线性叠加方法产生大量伪影，尤其是直达面波前的早到波；而应用加权叠加方法后，波形更清晰，噪声源分布不均匀引起的早到波明显被抑制。

对叠加后C¹函数不同分量间的相似性进行统计分析，相似性定义与式（3）相同，但不进行阈值处理。图6（a）、（b）表明，传统线性叠加方法未能抑制噪声源非均匀分布的影响，导致C¹函数相似性较低；而新叠加方法显著提高了叠加后C¹函数的相似性。本文还计算了C¹函数的信噪比（SNR）[图6（c）、（d）]，将SNR定义为信号窗口内最大振幅与信号窗口外数据均方根的比值。应用加权叠加方法后，SNR大于6的C¹函数数量显著增加，且这种改善在多个分量中均一致。

3.2 C²函数计算

上文已介绍了获取可靠C¹函数的方法，为计算C²函数奠定了坚实基础。在计算C²函数前，需考虑高阶面波的影响。C¹函数中的多模面波可能在C¹函数中产生基阶与高阶模式的交叉伪影。对图5（b）中的VSG计算其频散谱[图7（a）]，发现存在显著的高阶模式。假设基阶模式和高阶模式分别对应逆向和顺向椭圆质点运动，如以往研究所述[65,67,71]，通过对ZZ、RR、ZR和RZ分量进行适当相位校正后叠加，可抑制高阶模式。对ZR互相关施加正

π / 2

相位偏移，对RZ互相关施加负

π / 2

相位偏移。本文采用相位加权叠加方法[72]，充分利用其相位一致性，通过叠加ZZ、RR以及相位校正后的ZR和RZ互相关来抑制高阶面波[图7（b）]。为减轻高阶模式对C²函数的影响，利用滤波后的C¹函数计算C²函数。

图8为C ² 函数计算示意图。将停车场外的台站作为中间台站，分别计算其与停车场内第一轮、第二轮台站的C¹函数从而得到C¹，随后对这两组C¹函数进行互相关得到C²函数，通过选择性叠加得到最终C²函数的估计值，代表两个异步台站间的面波。选择依据为内部台站与外部台站连线形成的角度θ（图8）。如图8所示，C²函数的台站间距记为d₁，两组C¹函数的台站对间距分别记为d₂和d₃。对C¹函数进行互相关得到的信号实际传播距离为(d₃-d₂)，近似为d₁。θ越大，近似误差越大。本研究将θ设为30°，即仅使用与沿线方向夹角15°范围内的C¹台站对计算C²函数。根据本研究中的台站间距计算，对于处于最大角度边缘的台站对，其近似误差在5%以内。此外，对该角度范围内得到的C²信号进行叠加，可进一步减小近似误差。这种选择性叠加提高了C²函数的可靠性和质量[55]。

为进一步说明高阶面波对C²函数的影响，对比了C²函数的VSGs（图9）。将所有异步C²函数对按台站间距排序，并以3 m为间隔取平均值。所有波形均进行2~15 Hz带通滤波，以确保对比的公平性。观察发现，当未对C¹函数中的高阶模式进行滤波时，对应的C²函数会出现大量非物理波形[图9（a）]；而对C¹函数进行高阶面波滤波后，得到的C²函数仅显示清晰的基阶面波[图9（b）]，确保了后续频散测量的可靠性。

3.3 频散测量

经过上述步骤，获得了停车场内部台站间C¹函数的对称分量，包括C¹和C²函数。接下来，应用双台站频散分析法[7,73]获取任意互相关对的频散能量谱。该方法首先对两台站间的C¹函数进行时频分析，然后根据台站间距将时间轴转换为速度轴，从而得到台站间面波的相速度。在拾取双台站相速度频散曲线时，为解决频散谱中因相位模糊导致的多能量峰簇问题，本文利用多道频散分析法获取参考频散曲线，并在双台站频散谱中选取参考曲线附近的频散趋势。仅选择SNR大于6的C¹函数进行频散曲线拾取。

为研究C²方法获取的相速度测量值可能存在的误差，计算了第一轮观测中台站对间的C²函数并自动提取频散曲线，对比了同一路径上从C¹和C²函数中提取的相速度（V

C 1

和V

C 2

）。在所有有数据的频率点上，计算了相对速度偏差（

V C 2 - V C 1 V C 1

× 100%）。直方图[图10（a）]显示，80%的速度偏差在6%以内。同时观察到正速度偏差占比较高，这与前文提及的C²方法固有的台站路径近似误差有关。速度偏差的来源众多，包括信号质量、噪声源分布不完善等因素。误差结果表明速度偏差无显著系统偏差，这在一定程度上保证了结果的可靠性。尽管存在这些误差，但C²方法允许异步观测，可实现更密集的台站覆盖。图10（b）显示了通过C¹函数从两轮观测中获取的频散曲线以及通过C²函数从异步台站对中获取的频散曲线。C²函数的平均频散曲线与C¹函数的非常接近，表明C²结果可靠，进一步排除了与平均速度偏差超过4个SDs的数据，以确保数据在合理范围内。用于反演的最终频散曲线数据如图11所示。为满足平面波近似理论，研究仅保留波长不超过台站间距的数据[74]。C¹频散测量数量在低频段急剧减少[图11（a）]，主要原因是台站间距有限，难以获取长波长信号。异步台阵覆盖地理范围更广，台站对数量多且台站间距更大，有效增加了长波长频散测量数量。相速度在不同频率下存在显著扰动，且低频与高频间相速度变化较大[图11（b）]，表明研究区域地下介质复杂且不均匀。

3.4 面波层析成像

获取面波频散曲线后，通过面波层析成像可以得到地下三维横波（3D S波）速度结构。掌握面波走时数据后，可以选择任何可靠的面波层析成像算法，并比较不同成像结果，但此类比较超出了本文的研究范围。本研究采用面波频散直接反演法[75]进行面波层析成像。该方法采用快速行进算法进行射线追踪，考虑了面波沿非大圆路径传播的可能性，更适合复杂近地表介质成像。直接反演法通过正演计算初始模型的理论走时，然后迭代调整模型，以最小化观测走时与计算走时的差异。本研究利用面波射线路径和一维面波深度敏感核，建立3D S波速度扰动与面波走时残差之间的线性关系，从而可直接从面波频散测量数据中获取地下3D S波速度结构，省去了传统方法中获取二维相速度/群速度图的中间步骤。

3.5 工作流程总结

前文介绍了同步-异步背景噪声面波层析成像的关键步骤，图12展示了完整的工作流程图。传统背景噪声面波层析成像的主要步骤包括背景噪声数据预处理、C¹函数计算、频散测量和面波层析成像。本研究考虑了短期观测中噪声源分布不均匀对C¹函数的影响，增加了增强SPZs内噪声源的步骤；同时注意到C¹函数中存在的高阶面波会影响C¹频散分析和C²函数计算，因此加入了抑制高阶面波的步骤。为提高C²函数提取面波的精度，设计了SPZ区C²叠加步骤，排除SPZ区外的C²函数对面波提取的干扰。

需注意的是，本文设计的工作流程为一个框架，各步骤所用具体技术并未明确限定，这意味着每个步骤可采用不同方法。例如，本研究中使用的面波层析成像算法可替换为其他面波层析成像算法；同样，本研究通过极化滤波抑制高阶面波，也可探索其他方法，如基于高分辨率时频分析的抑制方法[76]。本研究的主要创新点在于设计了多个额外的数据处理步骤，这些步骤对于在短期观测中可靠获取C¹和C²函数至关重要。

4 结果

4.1 初始模型与深度敏感性

模型参数化通过在经纬度方向设置等距网格实现。但由于停车场并非严格沿东西方向分布，传统参数化方法会导致大量网格位于停车场外。这些网格无面波路径穿过，会降低反演系统的稳定性，并增加计算成本。为解决此问题，研究采用坐标旋转法，将测线近似旋转至东西方向，反演完成后恢复真实坐标。

场地36个钻孔中，3个包含S波速度测井结果[图13（a）]，以此为基础构建初始层状模型。模型由浅至深层厚逐渐增加，具体包括15个3 m厚的层、2个5 m厚的层和3个10 m厚的层。反演网格的水平间距设置为与台站间距相当，代表物理分辨率极限，具体经、纬度方向水平间距均为0.000 05°（约5.5 m）。

在层析成像中，密度和纵波（P波）速度均通过经验关系由横波速度确定[77]。初始层状模型的密度和纵波速度也采用相同经验关系推导，其对基阶瑞利波相速度的深度敏感核如图13（b）所示。结果表明，随着频率降低，面波对更深层介质的敏感性增强。本研究频散数据覆盖的频率范围为2～15 Hz，对60 m以上的介质敏感性良好，对更深层介质敏感性相对较弱。因此，将反演结果的解释限定在60 m以上深度。

4.2 棋盘格测试

为评估路径覆盖对地下结构的恢复能力，构建了两个三维棋盘格模型。第一个模型水平异常尺寸约为14 m × 17 m（图14），第二个模型水平异常尺寸约为23 m × 28 m（图15）。通过对实际数据反演所用初始模型施加10%扰动，得到棋盘格模型的速度。基于理论模型计算走时数据，其台站对和频率覆盖与实际数据相同。在模拟走时数据中加入1%高斯随机噪声以更接近真实情况。研究采用与实际反演相同的权重因子和迭代次数进行反演。恢复结果显示，在30 m深度处可清晰分辨出14 m × 17 m的异常，在大约55 m深度处最大可分辨出23 m × 28 m的异常，这与敏感核分析结果一致。

需注意的是，棋盘格测试结果受多种因素影响，如噪声水平、模型参数化、异常体分布等，因此不建议盲目依赖测试结果。应结合其他相关信息，对实际反演结果进行综合分析和解释。棋盘格测试结果可对分辨率进行初步评估，避免对反演结果过度解释。

4.3 S波速度模型

层析成像的目标函数包括数据残差和模型粗糙度正则化项，因此需确定合适的权重因子以平衡二者对目标函数的贡献。研究通过测试多个权重因子，对比反演结果的数据拟合度和模型平滑度，最终确定权重因子为1。经过10次迭代后，面波走时残差平均值从-0.0165 s降至-0.0006 s，且后续几乎不再下降。残差分布在0 s附近的集中度更高（图16），表明反演模型能有效捕捉并拟合频散测量数据，反演已成功收敛。

钻探结果表明，岩溶洞穴主要分布在40 m及以下区域。图17显示了反演3D模型在45 m、50 m、55 m深度的水平切片，场地西南和东北区域存在大范围的低速异常，与钻探结果确定的主要位置一致；而东南和西北区域为高速异常，表明基岩面较浅。结合敏感性分析和分辨率测试结果，模型难以精确识别水平和垂直尺寸均小于10 m的异常，尤其是在40 m以下深度。因此，期望速度模型完全揭示所有岩溶洞穴位置是不现实的。将钻探数据识别的洞穴位置叠加到速度切片上，除少数小洞穴未显示速度异常外，大多数洞穴位置均与低速异常相关。此前，施工团队曾尝试对废弃停车场的地下空洞进行灌浆处理，为建设新办公楼做准备，但处理效果未达到建筑安全标准。该区域的建设计划自2019年启动后一直停滞，截至2023年仍未开工。地下注入水泥浆对改善地下介质弹性强度的效果甚微，推测可能是由于速度模型中存在两个分布较广的低速区，这些低速区可能对应空洞和（或）高孔隙度砾石。这些观察表明，反演结果有效捕捉了地下岩溶特征的整体分布。

图18显示了两条速度剖面AA′ [图18（a）]和BB′ [图18（b）]，其在3D模型中的位置如图17所示，并叠加了与剖面相交区域的岩性钻探结果。钻探数据分析表明，停车场下伏地下介质可分为两层：沉积盖层和下伏基岩。在0～50 m深度范围内，S波速度从100 m·s^-1逐渐增加到400 m·s^-1；由于存在基岩，超过该深度范围，速度会显著增加。速度界面凹陷位置与钻探数据识别的岩溶洞穴位置高度吻合。部分岩溶洞穴未在速度剖面中显示的原因将在后续讨论中说明。

5 讨论

5.1 模型分辨率的限制因素

反演分辨率的一个重要方面是对不同尺寸速度异常的恢复能力。为研究台站密度对分辨率的影响，本文分析了从数据集中均匀排除一半台站后的反演结果，发现棋盘格模型的恢复效果显著下降（附录A中的图S1），稀疏台阵对地下低速异常的识别能力变差（附录A中的图S2）。另一方面，真实介质的复杂结构也会影响理想分辨率的实现。例如，图18（b）中#18号钻孔显示基岩深度约为40 m，而速度模型显示该位置的陡速度界面在50 m深度。垂直分辨率受台站密度影响较小，主要由面波频散反演的敏感性决定。面波波形反演有望提高其分辨率[77‒78]，此外，加入高阶面波数据也可提高反演的垂直分辨率[49,79‒82]。本文未讨论采用新的反演方法和新类型数据（如多阶瑞利波、勒夫波以及背景噪声水平分量与垂直分量谱比）来提高层析成像分辨率，未来研究将进一步探索这些方向。

考虑反演分辨率时还会产生一个问题：利用速度模型解释异常分布的准确性如何？场地有36个岩性钻孔，其中3个还包含S波速度测井结果。通过分析3个钻孔的数据，研究S波速度与岩性的关系。如图19所示，除砾石和石灰岩具有明显的高速特征、岩溶洞穴具有低速特征外，其他不同岩性间无显著速度差异。不同土层的速度存在波动，且这些速度差异并非恒定。考虑到上述反演分辨率的复杂性，本研究将结果解释限定在低速异常的空间分布和基岩深度变化上，以避免过度解释。

5.2 增加异步观测的优势

本研究提出的新型同步-异步观测系统，可利用有限台站数量实现研究区域的密集台阵观测。图20显示了本研究中2.0 Hz和5.4 Hz频率下面波的射线路径密度。为模拟传统同步观测，本研究特意将两轮观测的台站数量各减少一半。结果表明，加入异步观测后，面波射线路径密度显著增加，研究区域覆盖度提高。本研究仅通过两轮观测就成功实现了射线路径密度的大幅增加，理论上，增加观测轮次可进一步提高射线路径密度。新观测系统大幅减少了密集台阵观测所需的台站数量，从而显著降低了成本，这是传统同步观测几乎无法实现的。这种同步-异步观测系统有望成为未来密集台阵背景噪声成像研究的替代方案，尤其是在台站数量受限的情况下。

6 结论

本研究展示了短期同步-异步背景噪声层析成像在城市岩溶调查中的应用。通过同步-异步背景噪声观测系统，利用有限台站实现了研究区域的大范围空间覆盖。以C₁函数不同分量间的相似性作为叠加权重，成功抑制了噪声源分布不均匀的影响。采用多分量瑞利波叠加技术对C₁函数中的高阶模式进行滤波后，C₂方法可有效提取异步台站间的EGFs。通过背景噪声面波层析成像，获得了研究区域下方70 m深度内的3D S波速度结构。结果显示，东北部和西南部40 m以下区域存在显著的低速异常，结合钻孔资料推测这些异常与岩溶洞穴和（或）高孔隙度破碎砾石有关。结合C₁和C₂数据，比仅使用C₁函数，空间覆盖范围和分辨率都得到了显著提升。本文提出的同步-异步背景噪声层析成像方法即使在台站资源有限、复杂背景噪声条件下，也能实现快速、高分辨率的地震成像，有望为城市地质灾害防控提供有力支持。

参考文献

原文顺序 | 出版日期 | 本文引用

[1]	Weaver RL, Lobkis OI. Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies. Phys Rev Lett 2001;87(13):134301. . 10.1103/physrevlett.87.134301

[2]	Campillo M, Paul A. Long-range correlations in the diffuse seismic coda. Science 2003;299(5606):547‒9. . 10.1126/science.1078551

[3]	Snieder R. Extracting the Green’s function from the correlation of coda waves: a derivation based on stationary phase. Phys Rev E Stat Nonlin Soft Matter Phys 2004;69(4):046610. . 10.1103/physreve.69.046610

[4]	Sabra KG, Gerstoft P, Roux P, Kuperman WA, Fehler MC. Extracting time domain Green’s function estimates from ambient seismic noise. Geophys Res Lett 2005;32(3):L03310. . 10.1029/2004gl021862

[5]	Shapiro NM, Campillo M, Stehly L, Ritzwoller MH. High-resolution surface wave tomography from ambient seismic noise. Science 2005;307(5715):1615‒8. . 10.1126/science.1108339

[6]	Wapenaar K, Draganov D, Snieder R, Campman X, Verdel A. Tutorial on seismic interferometry: part 1—basic principles and applications. Geophysics 2010;75(5):75A195‒209. . 10.1190/1.3457445

[7]	Yao H, van der Hilst RD, de Hoop MV. Surface-wave array tomography in southeast Xizang from ambient seismic noise and two-station analysis—I. Phase velocity maps. Geophys J Int 2006;166(2):732‒44. . 10.1111/j.1365-246x.2006.03028.x

[8]	Yang Y, Ritzwoller MH, Levshin AL, Shapiro NM. Ambient noise Rayleigh wave tomography across Europe. Geophys J Int 2007;168(1):259‒74. . 10.1111/j.1365-246x.2006.03203.x

[9]	Rawlinson N, Salmon M, Kennett BLN. Transportable seismic array tomography in southeast Australia: illuminating the transition from Proterozoic to Phanerozoic lithosphere. Lithos 2014;189:65‒76. . 10.1016/j.lithos.2013.06.001

[10]	Bao X, Song X, Li J. High-resolution lithospheric structure beneath China’s mainland from ambient noise and earthquake surface-wave tomography. Earth Planet Sci Lett 2015;417:132‒41. . 10.1016/j.epsl.2015.02.024

[11]	Luo S, Yao H, Li Q, Wang W, Wan K, Meng Y, et al. High-resolution 3D crustal S wave velocity structure of the Middle-Lower Yangtze River metallogenic belt and implications for its deep geodynamic setting. Sci China Earth Sci 2019;62(9):1361‒78. . 10.1007/s11430-018-9352-9

[12]	Nimiya H, Ikeda T, Tsuji T. Three-dimensional S wave velocity structure of central Japan estimated by surface-wave tomography using ambient noise. J Geophys Res Solid Earth 2020;125:e2019JB019043. . 10.1029/2019jb019043

[13]	Gaite B, Iglesias A, Villaseñor A, Herraiz M, Pacheco JF. Crustal structure of Mexico and surrounding regions from seismic ambient noise tomography. Geophys J Int 2012;188(3):1413‒24. . 10.1111/j.1365-246x.2011.05339.x

[14]	Haned A, Stutzmann E, Schimmel M, Kiselev S, Davaille A, Yelles-Chaouche A. Global tomography using seismic hum. Geophys J Int 2016;204(2):1222‒36. . 10.1093/gji/ggv516

[15]	Li C, Yao H, Fang H, Huang X, Wan K, Zhang H, et al. 3D near-surface shear wave velocity structure from ambient-noise tomography and borehole data in the Hefei urban area, China. Seismol Res Lett 2016;87(4):882‒92. . 10.1785/0220150257

[16]	Sun R, Li J, Yan Y, Liu H, Bai L, Chen Y. Three-dimensional urban subsurface space tomography with dense ambient noise seismic array. Surv Geophys 2024;45(3):819‒43. . 10.1007/s10712-023-09819-3

[17]	Nakata N. Near-surface S-wave velocities estimated from traffic-induced Love waves using seismic interferometry with double beamforming. Interpretation 2016;4(4):SQ23‒31. . 10.1190/int-2016-0013.1

[18]	Xia S, Zhang C, Cao J. Ambient noise tomography for coral islands. Engineering 2023;25 Jun:182‒93. . 10.1016/j.eng.2021.09.022

[19]	Carvalho J, Silveira G, Dumont S, Ramalho R. 3D-ambient noise surface wave tomography of Fogo volcano, Cape Verde. J Volcanol Geotherm Res 2022;432:107702. . 10.1016/j.jvolgeores.2022.107702

[20]	Lehujeur M, Vergne J, Schmittbuhl J, Zigone D, LeChenadec A. EstOF Team. Reservoir imaging using ambient noise correlation from a dense seismic network. J Geophys Res Solid Earth 2018;123(8):6671‒86. . 10.1029/2018jb015440

[21]	Cheng F, Xia J, Ajo-Franklin JB, Behm M, Zhou C, Dai T, et al. High-resolution ambient noise imaging of geothermal reservoir using 3C dense seismic nodal array and ultra-short observation. J Geophys Res Solid Earth 2021;126(8): e2021JB021827. . 10.1029/2021jb021827

[22]	Zhou C, Xia J, Pang J, Cheng F, Chen X, Xi C, et al. Near-surface geothermal reservoir imaging based on the customized dense seismic network. Surv Geophys 2021;42(3):673‒97. . 10.1007/s10712-021-09642-8

[23]	Gu N, Wang K, Gao J, Ding N, Yao H, Zhang H. Shallow crustal structure of the Tanlu fault zone near Chao Lake in eastern China by direct surface wave tomography from local dense array ambient noise analysis. Pure Appl Geophys 2019;176(3):1193‒206. . 10.1007/s00024-018-2041-4

[24]	Guan B, Mi B, Zhang H, Liu Y, Xi C, Zhou C. Selection of noise sources and short time passive surface wave imaging—a case study on fault investigation. J Appl Geophys 2021;194:104437. . 10.1016/j.jappgeo.2021.104437

[25]	Ning L, Dai T, Liu Y, Xi C, Zhang H, Zhou C. Application of multichannel analysis of passive surface waves method for fault investigation. J Appl Geophys 2021;192:104382. . 10.1016/j.jappgeo.2021.104382

[26]	Park CB, Miller RD. Roadside passive multichannel analysis of surface waves (MASW). J Eng Environ Geophys 2008;13(1):1‒11. . 10.2113/jeeg13.1.1

[27]	Nakata N, Snieder R, Tsuji T, Larner K, Matsuoka T. Shear wave imaging from traffic noise using seismic interferometry by cross-coherence. Geophysics 2011;76(6):SA97‒SA106. . 10.1190/geo2010-0188.1

[28]	Cheng F, Xia J, Xu Z, Hu Y, Mi B. Frequency-wavenumber (FK)-based data selection in high-frequency passive surface wave survey. Surv Geophys 2018;39(4):661‒82. . 10.1007/s10712-018-9473-3

[29]	Liu Y, Xia J, Xi C, Dai T, Ning L. Improving the retrieval of high-frequency surface waves from ambient noise through multichannel-coherency-weighted stack. Geophys J Int 2021;227(2):776‒85. . 10.1093/gji/ggab253

[30]	Mi B, Xia J, Tian G, Shi Z, Xing H, Chang X, et al. Near-surface imaging from traffic-induced surface waves with dense linear arrays: an application in the urban area of Hangzhou, China. Geophysics 2022;87(2):B145‒58. . 10.1190/geo2021-0184.1

[31]	Cheng F, Xia J, Shen C, Hu Y, Xu Z, Mi B. Imposing active sources during high frequency passive surface-wave measurement. Engineering 2018;4(5):685‒93. . 10.1016/j.eng.2018.08.003

[32]	Waltham T, Bell FG, Culshaw MG, Knez M, Slabe T. Sinkholes and subsidence: karst and cavernous rocks in engineering and construction. Berlin: Springer; 2005. . 10.3986/ac.v35i1.551

[33]

Weary D. The cost of karst subsidence and sinkhole collapse in the United States compared with other natural hazards. In: Proceedings of the Fourteenth Multidisciplinary Conference; 2015 Oct 5‒6; Rochester, MN, USA. Tampa: University of South Florida Tampa Library; 2015. p. 433‒46. . 10.5038/9780991000951.1062

[34]	Zhang R, Ai T, Ren L, Li G. Failure characterization of three typical coal-bearing formation rocks using acoustic emission monitoring and X-ray computed tomography techniques. Rock Mech Rock Eng 2019;52(6):1945‒58. . 10.1007/s00603-018-1677-9

[35]	Dong L, Tong X, Ma J. Quantitative investigation of tomographic effects in abnormal regions of complex structures. Engineering 2021;7(7):1011‒22. . 10.1016/j.eng.2020.06.021

[36]	Dong L, Pei Z, Xie X, Zhang Y, Yan X. Early identification of abnormal regions in rock-mass using traveltime tomography. Engineering 2023;22(3):191‒200. . 10.1016/j.eng.2022.05.016

[37]	Hiltunen DR, Cramer BJ. Application of seismic refraction tomography in karst terrane. J Geotech Geoenviron Eng 2008;134(7):938‒48. . 10.1061/(asce)1090-0241(2008)134:7(938)

[38]	Aati AHA, Shabaan SH. Detection of karstic limestone bedrock by shallow seismic refraction in an area west of Assiut, middle Egypt. Lead Edge 2013;32(3):316‒22. . 10.1190/tle32030316.1

[39]	Peng D, Cheng F, Liu J, Zong Y, Yu M, Hu G, et al. Joint tomography of multi cross-hole and borehole-to-surface seismic data for karst detection. J Appl Geophys 2021;184:104252. . 10.1016/j.jappgeo.2020.104252

[40]	Zhang J, Liu S, Chen Q, Wang B, Ren C. Application of cross-borehole integrated geophysical methods for the detailed investigation of karst in urban metro construction. JEEG 2019;24(4):525‒36. . 10.2113/jeeg24.04.03.525

[41]	Shangxin F, Yufei Z, Yujie W, Shanyong W, Ruilang C. A comprehensive approach to karst identification and groutability evaluation—a case study of the Dehou reservoir, SW China. Eng Geol 2020;269:105529. . 10.1016/j.enggeo.2020.105529

[42]	Gan F, Han K, Lan F, Chen Y, Zhang W. Multi-geophysical approaches to detect karst channels underground—a case study in Mengzi of Yunnan Province, China. J Appl Geophys 2017;136:91‒8. . 10.1016/j.jappgeo.2016.10.036

[43]	Foudili D, Bouzid A, Berguig MC, Bougchiche SS, Abtout A, Guemache MA. Investigating karst collapse geohazards using magnetotellurics: a case study of M’rara basin, Algerian Sahara. J Appl Geophys 2019;160:144‒56. . 10.1016/j.jappgeo.2018.11.011

[44]	Kiernan M, Jackson D, Montgomery J, Anderson JB, McDonald BW, Davis KC. Characterization of a karst site using electrical resistivity tomography and seismic full waveform inversion. JEEG 2021;26(1):1‒11. . 10.32389/jeeg20-045

[45]	Xia J, Chen C, Li P, Lewis M. Delineation of a collapse feature in a noisy environment using a multichannel surface wave technique. Geotechnique 2004;54(1):17‒27. . 10.1680/geot.54.1.17.36326

[46]	Parker Jr EH, Hawman RB. Multi-channel analysis of surface waves (MASW) in karst terrain, southwest Georgia: implications for detecting anomalous features and fracture zones. J Environ Eng Geophys 2012;17(3): 129‒50. . 10.2113/jeeg17.3.129

[47]	Chang JP, de Ridder SAL, Biondi BL. High-frequency Rayleigh-wave tomography using traffic noise from Long Beach, California. Geophysics 2016;81(2):B43‒53. . 10.1190/geo2015-0415.1

[48]	Wang Y, Lin FC, Schmandt B, Farrell J. Ambient noise tomography across Mount St. Helens using a dense seismic array. J Geophys Res Solid Earth 2017;122(6):4492‒508. . 10.1002/2016jb013769

[49]	Fu L, Pan L, Li Z, Dong S, Ma Q, Chen X. Improved high-resolution 3D vs model of Long Beach, CA: inversion of multimodal dispersion curves from ambient noise of a dense array. Geophys Res Lett 2022;49:e2021GL097619. . 10.1029/2021gl097619

[50]	Stehly L, Campillo M, Froment B, Weaver RL. Reconstructing Green’s function by correlation of the coda of the correlation (C³) of ambient seismic noise. J Geophys Res 2008;113(B11):B11306. . 10.1029/2008jb005693

[51]	Curtis A, Halliday D. Source-receiver wave field interferometry. Phys Rev E Stat Nonlin Soft Matter Phys 2010;81(4):046601. . 10.1103/physreve.81.046601

[52]	Froment B, Campillo M, Roux P. Reconstructing the Green’s function through iteration of correlations. C R Geosci 2011;343(8‒9):623‒32.

[53]	Spica Z, Perton M, Calò M, Legrand D, Córdoba-Montiel F, Iglesias A. 3-D shear wave velocity model of Mexico and south US: bridging seismic networks with ambient noise cross-correlations (C¹) and correlation of coda of correlations (C³). Geophys J Int 2016;206(3):1795‒813. . 10.1093/gji/ggw240

[54]	Ansaripour M, Rezapour M, Saygin E. Shear wave velocity structure of Iranian plateau: using combination of ambient noise cross-correlations (C¹) and correlation of coda of correlations (C³). Geophys J Int 2019;218 (3):1919‒38. . 10.1093/gji/ggz257

[55]	Chen Y, Saygin E. Empirical Green’s function retrieval using ambient noise source-receiver interferometry. J Geophys Res Solid Earth 2020;125(2): e2019JB018261. . 10.1029/2019jb018261

[56]	Zhang S, Feng L, Ritzwoller MH. Three-station interferometry and tomography: coda versus direct waves. Geophys J Int 2020;221(1):521‒41. . 10.1093/gji/ggaa046

[57]	Rao H, Luo Y, Zhao K, Yang Y. Extracting surface wave dispersion curves from asynchronous seismic stations: method and application. Geophys J Int 2021;226(2):1148‒58. . 10.1093/gji/ggab153

[58]	Chen Y, Saygin E, Kennett B, Qashqai MT, Hauser J, Lumley D, et al. Next generation seismic model of the Australian crust from synchronous and asynchronous ambient noise imaging. Nat Commun 2023;14(1):1192. . 10.1038/s41467-023-36514-z

[59]	Boschi L, Weemstra C. Stationary-phase integrals in the cross correlation of ambient noise. Rev Geophys 2015;53(2):411‒51. . 10.1002/2014rg000455

[60]	Liu Y, Xia J, Cheng F, Xi C, Shen C, Zhou C. Pseudo-linear-array analysis of passive surface waves based on beamforming. Geophys J Int 2020;221(1):640‒50. . 10.1093/gji/ggaa024

[61]	Cheng F, Xia J, Luo Y, Xu Z, Wang L, Shen C, et al. Multichannel analysis of passive surface waves based on crosscorrelations. Geophysics 2016;81(5): EN57‒66. . 10.1190/geo2015-0505.1

[62]	Zhang H, Mi B, Xi C, Liu Y, Guan B, Ning L. Extraction of surface-wave phase velocities from ambient noise in the presence of local noise sources based on matched-field processing. J Appl Geophys 2022;204:104755. . 10.1016/j.jappgeo.2022.104755

[63]	Ning L, Xia J, Dai T, Liu Y, Zhang H, Xi C. High-frequency surface-wave imaging from traffic-induced noise by selecting in-line sources. Surv Geophys 2022;43(6):1873‒99. . 10.1007/s10712-022-09723-2

[64]	Liu Y, Xia J, Xi C, Zhang H, Guan B, Dai T, et al. Enhancing noise sources in stationary-phase zones for accurate phase-velocity estimation of high frequency surface waves. Geophysics 2023;88(1):L1‒9. . 10.1190/geo2021-0696.1

[65]	Nayak A, Thurber CH. Using multicomponent ambient seismic noise cross correlations to identify higher mode Rayleigh waves and improve dispersion measurements. Geophys J Int 2020;222(3):1590‒605. . 10.1093/gji/ggaa270

[66]	Lin FC, Moschetti MP, Ritzwoller MH. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps. Geophys J Int 2008;173(1):281‒98. . 10.1111/j.1365-246x.2008.03720.x

[67]	Gribler G, Mikesell TD. Methods to isolate retrograde and prograde Rayleigh wave signals. Geophys J Int 2019;219(2):975‒94. . 10.1093/gji/ggz341

[68]	Herrmann RB. Computer programs in seismology: an evolving tool for instruction and research. Seismol Res Lett 2013;84(6):1081‒8. . 10.1785/0220110096

[69]	Bensen GD, Ritzwoller MH, Barmin MP, Levshin AL, Lin F, Moschetti MP, et al. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophys J Int 2007;169(3):1239‒60. . 10.1111/j.1365-246x.2007.03374.x

[70]	Park CB, Miller RD, Xia J. Imaging dispersion curves of surface waves on multi channel record. In: SEG technical program expanded abstracts 1998. Houston: Society of Exploration Geophysicists; 1998. p. 1377‒80. . 10.1190/1.1820161

[71]	Gribler G, Liberty LM, Mikesell TD, Michaels P. Isolating retrograde and prograde Rayleigh-wave modes using a polarity mute. Geophysics 2016;81(5): V379‒85. . 10.1190/geo2015-0683.1

[72]	Schimmel M, Paulssen H. Noise reduction and detection of weak, coherent signals through phase-weighted stacks. Geophys J Int 1997;130(2):497‒505. . 10.1111/j.1365-246x.1997.tb05664.x

[73]	Yao H, Gouédard P, Collins JA, McGuire JJ, van der Hilst RD. Structure of young East Pacific rise lithosphere from ambient noise correlation analysis of fundamental- and higher-mode Scholte-Rayleigh waves. C R Geosci 2011;343(8‒9):571‒83.

[74]	Luo Y, Yang Y, Xu Y, Xu H, Zhao K, Wang K. On the limitations of interstation distances in ambient noise tomography. Geophys J Int 2015;201(2):652‒61. . 10.1093/gji/ggv043

[75]	Fang H, Yao H, Zhang H, Huang YC, van derHilst RD. Direct inversion of surface wave dispersion for three-dimensional shallow crustal structure based on ray tracing: methodology and application. Geophys J Int 2015;201(3):1251‒63. . 10.1093/gji/ggv080

[76]	Luo Y, Xia J, Miller RD, Xu Y, Liu J, Liu Q. Rayleigh-wave mode separation by high-resolution linear Radon transform. Geophys J Int 2009;179(1):254‒64. . 10.1111/j.1365-246x.2009.04277.x

[77]	Brocher TM. Empirical relations between elastic wavespeeds and density in the earth’s crust. Bull Seismol Soc Amer 2005;95(6):2081‒92. . 10.1785/0120050077

[78]	Li J, Feng Z, Schuster G. Wave-equation dispersion inversion. Geophys J Int 2017;208(3):1567‒78. . 10.1093/gji/ggw465

[79]	Pan Y, Gao L, Bohlen T. Random-objective waveform inversion of 3D‒9C shallow-seismic field data. J Geophys Res Solid Earth 2021;126: e2021JB022036. . 10.1029/2021jb022036

[80]	Xia J, Miller RD, Park CB, Tian G. Inversion of high frequency surface waves with fundamental and higher modes. J Appl Geophys 2003;52(1):45‒57. . 10.1016/s0926-9851(02)00239-2

[81]	Pan L, Chen X, Wang J, Yang Z, Zhang D. Sensitivity analysis of dispersion curves of Rayleigh waves with fundamental and higher modes. Geophys J Int 2019;216(2):1276‒303. . 10.1093/gji/ggy479

[82]	Liu Q, Chen X, Gao L, Yu Z, Chen J. Direct image dissimilarity inversion of ambient noise multimodal dispersion spectrograms. Bull Seismol Soc Amer 2023;113(5):1960‒81. . 10.1785/0120230012