随着我国基础设施建设的不断发展以及城市化进程的快速推进，准确描绘出浅地表天然存在或人工制造的不稳定区域(通常表现为低速、低密度特征)，如地下空洞、道路塌陷区、煤层采空区等，对于保证陆地交通等基础设施建设与人们安全生产生活具有重要意义^[1-2]。同时，构建高精度浅地表模型，对解决复杂地质条件下煤田地震勘探资料处理中的静校正问题至关重要^[3-6]。然而近地表岩土层性质复杂且探测目标尺寸小，对地球物理方法成像的分辨率与精度提出极高的要求，加剧了探测难度。利用对横波速度相当敏感且能量主要集中在近地表的地震面波进行近地表速度建模自本世纪以来在勘探地球物理界获得了极大关注，面波多道分析法(Multichannel Analysis of Surface Waves，MASW)^[7-8]即是其中的代表性方法之一，在探测如空洞等浅地表低速异常体方面得到了广泛的应用^[9-11]，取得了较为满意的结果。近年来出现的逆散射面波分析方法^[12]与多窗口面波分析方法^[13]等应用于相似的实际案例，成像分辨率与精度得到进一步提升。然而这类基于频散曲线反演方法的假设条件较为苛刻，作为输入数据的频散曲线是当前空间窗口介质的综合响应，当低速体尺寸较小或空间窗口较长时，频散曲线对异常体的敏感度降低而难以对其重构。另外，反演寻优基于层状模型，这进一步增加了反演结果的不确定性。

通过利用地震波运动学与动力学在内的全波场信息来直接进行二维或三维反演建模的全波形反演(Full Waveform Inversion，FWI)方法自A. Tarantola^[14]提出以来，不断发展，但多聚焦于体波信号的反演^[15-16]。利用面波信号进行反演的面波FWI方法在浅地表地球物理界的研究与应用在近10年来发展迅速^[17-18]，已被证明相较于上述基于频散分析的面波方法，成像精度与分辨率有大幅提升^[19]。根据波形形成机制的不同，面波可分为瑞雷波(P波与SV波相互干涉)与勒夫波(SH波相互干涉)两类。而勒夫波与纵波无关的特点，使得勒夫波FWI的系统更为简单，不仅能够保证计算效率还能减弱不同参数之间的串扰^[20]，但由于其水平激发源的限制，只有极少的实际应用案例。Pan Yudi等^[21]采用该方法处理了采集于美国怀俄明州的一单炮SH波记录，证明了测区高速表层的存在，E. Dokter^[20]、Yan Yingwei^[22]等分别基于该技术来判断斯洛伐克的某村庄处是否存在古遗址地基，均给出了否定的结论；M.I. Alam ^[23]则完成了对俄克拉荷马大学校园内埋深约1.5 m、直径约0.8 m管道的探测任务；R. Mecking等^[24]通过增强勒夫波数据中的SH型折射信息成功探测到了位于土耳其佩加蒙埋深约30 m 的一个空洞；Guan Jianbo等^[25]将该方法应用于考古地球物理中，清晰揭示了一个高1.6 m、长4.7 m的古墓。需要注意的是，当空洞或地下隧道等周围存在诸如混凝土等特高速介质时，勒夫波 FWI极大可能会失效^[24]，这在实际案例中得到了验证^[26]。瑞雷波由于垂向激发方式简单且不受上述限制的影响，针对瑞雷波FWI的研究及其在低速异常体探测中的实际应用更为广泛。K.T. Tran等^[27]采用该技术成功地完成了位于弗罗里达大学与美国纽贝里的两个地下空洞探测任务，横波速度与纵波速度反演剖面清晰地展现了空洞的位置与尺寸；之后，K.T. Tran等^[28]又成功表征了美国拉斯维加斯山谷附近的一处低速沉积夹层；Wang Yao等^[29]对美国亚利桑那州索诺兰沙漠的一处人工隧道进行了探测，反演结果与实际资料匹配较好。但上述研究或将密度、或将纵波速度设置为常值而忽略它们的影响。F. Wittkamp等^[30]在对卡尔斯鲁厄机场附近采集的数据进行处理时，对横、纵波速度与密度进行同步更新，横波速度模型在地表处清晰显示了一个回填渠道；之后，Pan Yudi等^[31]开发了一种随机目标函数瑞雷波FWI策略并应用于同一野外数据，同样对三参数进行更新，结果显示回填渠道的形态与边界更为清晰。另外，三维瑞雷波 FWI的研究与应用也已出现^[32-33]，但其巨大的计算消耗与内存需求使其在实际工程中的推广应用十分受限。

面波FWI作为一种高精度的浅地表成像方法，还未在应用地球物理界获得推广应用，尤其是国内的应用研究较少。本文通过数值测试与实际数据应用来展示瑞雷波多参数FWI在浅地表空洞探测中的可行性和有效性，期望能够为相同尺度的工程实践应用提供参考。

1   方法原理

FWI的核心理念是期望通过最小化目标函数$\varPhi$找到一个可以对记录的地震信号(包括波形与相位)进行准确描述的介质模型，该模型被认为是对应于观测信号的真实模型。与此同时，在实际情况中通常采集多偏移距数据以实现对观测区介质的充分覆盖，由此$\varPhi$可定义为：

式中：${\boldsymbol{m}}$为模型参数；$s$与N分别为炮序号与总炮数；$ {{\boldsymbol{d}}^{{\text{est}}}} $与$ {{\boldsymbol{d}}^{{\text{rec}}}} $分别为单炮预测数据与记录数据。采用计算量较低、基于梯度下降的局部优化算法来应对FWI这个大型非线性寻优问题：

式中：i+1与i为迭代次数；${\alpha _i}$与$ \Delta {{\boldsymbol{m}}_i} $分别为第i次迭代时模型的扰动步长与模型更新方向。由预处理共轭梯度算法^[13-14]确定的模型更新方向：

式中：$ {{\boldsymbol{H}}^{ - 1}}({{\boldsymbol{m}}_i}) $、$\nabla \varPhi {\text{(}}{{\boldsymbol{m}}_i}{\text{)}}$和${\beta _i} $分别为第i次迭代时的拟Hessian矩阵的逆、目标函数梯度与模型更新方向。伴随状态法^[34]被用来计算目标函数的梯度以避免显式计算方式^[21]的巨大计算量，其仅需要以数据残差作为震源额外多做一次正演，并由该次正演获得的所谓伴随波场与正传波场进行互相关即完成计算任务。这里直接给出$\nabla \varPhi {\text{(}}{\boldsymbol{m}}{\text{)}}$关于横、纵波速度与密度的具体形式，推导过程请参阅C. Castellanos等^[35]的工作：

式中：${v_{\text{S}}}$、${v_{\text{P}}}$与$\rho $分别为模型横波、纵波速度与密度；$T$为记录数据采样点数；$ v_x^{{\text{for}}} $与$v_{\textit{z}}^{{\text{for}}}$为正传波场的各速度分量；$ v_x^{{\text{adj}}} $、$v_{\textit{z}}^{{\text{adj}}}$与$ \sigma _{xx}^{{\text{adj}}} $、$\sigma _{x{\textit{z}}}^{{\text{adj}}}$、$\sigma _{{\textit{z}}{\textit{z}}}^{{\text{adj}}}$分别为反传波场的各速度分量与各应力分量。

目标函数的Hessian矩阵具有均衡梯度振幅、加强模型深部照明、压制地表假象以及提升小尺度异常的表征能力^{[17-18, 22, 25]}，为了避免直接对其计算而导致的计算及存储溢出的问题，可通过二阶近似构造的方式，对正传波场进行互相关来获取一个功能十分相近的预处理算子，即拟Hessian矩阵。这里同样直接给出$ {\boldsymbol{H}}({\boldsymbol{m}}) $关于横、纵波速度与密度的具体形式：

其中，${\boldsymbol{{H}}}{\text{(}}{\boldsymbol{m}}{\text{)}}{|_{{\boldsymbol{m}} = {v_{\rm{S}}}}}$、${\boldsymbol{{ H}}}{\text{(}}{\boldsymbol{m}}{\text{)}}{|_{{\boldsymbol{m}} = {v_{\rm{P}}}}}$与${\boldsymbol{{ H}}}{\text{(}}{\boldsymbol{m}}{\text{)}}{|_{{\boldsymbol{m}} = \rho }}$分别为横波速度、纵波速度与密度的拟Hessian算子。具体推导过程请参阅C. Shin等^[36]的工作。之后，直接对$ {\boldsymbol{H}}({\boldsymbol{m}}) $取倒数，便可得到$ {{\boldsymbol{H}}^{ - 1}}({\boldsymbol{m}}) $。为了避免${\boldsymbol{{ H}}}{\text{(}}{\boldsymbol{m}}{\text{)}}$中某些近零值导致$ {{\boldsymbol{H}}^{ - 1}}({\boldsymbol{m}}) $量级的陡增现象，对其进行调整：

式中：$\gamma $为取值范围为${10^{ - 1}}$~${10^{ - 5}}$的超参数，在后文的合成测试及野外案例中其值均为${10^{ - 5}}$。

另外，FWI技术在实际数据中应用时，存在震源子波未知、数据包含低、高频噪声以及波场正演程序与实际数据采集存在维数差等问题，相应的解决技术策略将在第3节实际案例中进行详细介绍。

2   合成测试

本节通过一个合成测试验证瑞雷波多参数FWI在探测地下小尺度低速异常体中的有效性。所建立的合成模型如图1a、图1b与图1c所示，其尺寸为10 m(垂向)×28 m(横向)，以0.5 m的方形网格步长被离散为20×56的网格维度。模型分为两层，表层厚度2 m，底层为半无限空间。表层与底层的横波速度分别为150 m/s与230 m/s(图1a)；纵波速度为横波速度的2倍 (图1b)；表层与底层的密度分别为1.2 g/cm³、1.84 g/cm³ (图1c)。在底层埋深4 m处嵌入一个尺寸2.5 m(垂向)×4 m(横向)，横、纵波速度与密度分别为80 m/s，160 m/s与0.64 g/cm³的低速体，即为探测目标空洞(图1a、图1b与图1c中白色实线圈定的区域)。在整个模型空间中，密度、纵波速度均与横波速度呈现正向相关，这符合近地表介质的一般属性。接收排列布设于地表，以1 m道间距与2 m炮间距均匀布设25道垂直分量检波器与15个激发点。地震记录的正演模拟采用空间精度10阶、时间精度2阶的交错网格差分算法^[37]，震源子波为主频20 Hz且延迟时间0.02 s的Ricker子波，采样率0.25 ms，记录时长0.4 s。

图  1  合成测试中的真实模型与初始模型

Figure  1.  Actual and initial models in the synthetic test

下载: 全尺寸图片 幻灯片

图2a和图2b分别给出了炮点位于模型左部(炮序号1)与右部(炮序号15)的观测地震数据，从中可观察到面波强烈的频散特性，此特性在反演过程中极易引起周波跳跃问题，因此本文在FWI中加入了多尺度反演策略^[38]。图3a为所有地震记录各道数据进行频谱分析后得到的叠加振幅谱，在第10—第15道处，低频能量明显增强，指示异常体的存在。振幅谱中频带能量主要集中于5~60 Hz，在此反演策略分为中低频5~35 Hz、全频带5~65 Hz两个尺度。图3b为第1炮数据的面波频散谱，基于半波长法则建立的1D梯度型横波速度初始模型如图1d所示，初始横波速度由地表的150 m/s以4 m/s的步长线性增加至模型底部的226 m/s；与真实模型一样，初始纵波为初始横波速度的2倍 (图1e)；初始密度由地表的1.2 g/cm³以0.032 g/cm³的步长线性增加至模型底部的1.808 g/cm³ (图1f)。各初始模型未包含空洞的任何先验信息。在反演进程中，认为震源子波已知，当迭代的次数达到31次或当相邻两次迭代的目标函数的下降值小于初始值的0.001倍时，迭代停止。此次合成测试完整的FWI代码是在普通笔记本电脑(Intel(R) Core(TM) i7-9750H，6 cores 2.6 GHz, 18-GB RAM)运行的，整个反演进程共耗时约1.52 h。

图  2  合成测试中的记录数据

Figure  2.  Recorded data in the synthetic test

下载: 全尺寸图片 幻灯片

图  3  合成测试中记录数据的频谱与频散谱

Figure  3.  Spectrum and dispersion spectrum images of recorded data in the synthetic test

下载: 全尺寸图片 幻灯片

图4展示了归一化目标函数下降情况，其收敛速度很快，经10次迭代便已下降至0.2左右，最终收敛至0.1左右，而未收敛至0的原因将在后面进行分析。图5a分别对比了第1炮与第15炮的观测波形与最终反演结果获得的预测波形，二者之间的拟合程度极高，与基于初始模型的波形拟合情况相比(图5b)，主要震相的匹配度(包括振幅与相位)得到了很大程度的提高。

图  4  合成测试中的目标函数下降曲线

Figure  4.  Descending curve of the objective function in the synthetic test

下载: 全尺寸图片 幻灯片

图  5  合成测试中FWI反演结果和初始模型的预测数据与观测数据的波形对比

Figure  5.  Comparison of observed waveforms with waveforms predicted using the final FWI results and initial models in the synthetic test

下载: 全尺寸图片 幻灯片

图6分别给出了横波速度、纵波速度以及密度的反演结果，与初始模型相比，3个参数均得到了明显更新，均较好地恢复了表层的厚度及参数值，虽然加入了预处理技术，但由于震源、检波器与介质之间强烈的耦合效应导致地表附近产生一些明显假象。与此同时，横波速度最为准确地刻画出了空洞的位置与尺寸，恢复出真实速度，纵波速度模型次之，密度模型最差。这也说明了在瑞雷波FWI中，面波对于横波速度最为敏感。因此，在实际案例研究中应更多关心横波速度的反演结果。对于底层，三者也均得到了较好的恢复，但由于低速体的存在加剧了反演多解性程度，导致其下方各参数均出现了过度估计，影响了周围参数的准确更新，使横波速度与密度中均出现了一些小尺度的低速或高速体的假象，这也是波形未能完全拟合(图5a为波形隔道抽稀显示结果)、目标函数(图4)未能收敛至0的原因。另外，由于有限的炮检照明，导致在反演结果的左、右边界出现了“拖尾”现象^[19]。

图  6  合成测试中FWI方法的反演结果

Figure  6.  FWI results of in the synthetic test

下载: 全尺寸图片 幻灯片

图7详细比较了位于模型中部第13道处横波速度、纵波速度与密度的真实模型、初始模型和反演结果的1D剖面，便于更为直观地展示以上观点。初始模型除空洞所在区域外与真实模型较为接近，因此这些区域反演结果的参数更新并不明显，而空洞处的横波速度得到了很高程度的更新并与真实模型十分接近(图7a)，说明了所发展的方法具有在小范围模型偏差内对参数进行准确重构的能力。虽然密度与纵波速度受面波敏感度的影响对异常体的位置及参数的恢复程度并不十分令人满意，但并不能忽视它们在反演进程中的同步更新对横波速度模型中异常体的重建精度所作出的贡献^[39]。

图  7  合成测试中FWI反演结果的1D剖面

Figure  7.  1D profiles of the FWI results in the synthetic test

下载: 全尺寸图片 幻灯片

为了更加有力地展示FWI方法的优越性，将其与基于面波频散曲线的MASW方法的反演结果进行比较。在模型中保持空洞位于中心位置，将速度模型(图1a)向左右两侧分别扩展10 m，以适应MASW方法数据采集所需的接收排列滚动模式。为避免近场效应与确保探测深度^[7]，设置最小和最大偏移距分别为10 m和30 m，即接收排列为20 m；同时，考虑到空洞尺寸较小，设置激发点距尽可能小，每次向前滚动距离为1 m，充分观测空洞。基于图1a速度模型，模拟地震记录中没有明显的高阶频散能量(图3b)，因此选用对基阶模态具有良好效果的相移法^[40]提取频散能量。根据频散能量聚焦性与空间采样定理，图8a展示了全部28条频散曲线。从图中可以看出，在第8—第20个频散测量点范围，相速度呈现出低值特征，约为200 m/s，暗示存在低速异常。初始横波速度模型为图1d，采用阻尼最小二乘方法^[41]，反演结果如图8b所示，揭示了一个中心位于(13.5 m, 3.5 m)的低速区，与真实空洞水平中心位置一致，但比真实空洞浅约2 m。由于低速异常体导致面波波场复杂化以及反演的多解性，表层并没有被完全揭示。同时MASW方法具有横向平均效应，其反演结果模糊了空洞形态，并对低速的反馈极其有限，仅在正常值230 m/s基础上下降了约12%(图8c)。这与频散谱代表的计算频散所用排列长度(20 m)下方介质的平均效应以及空洞尺寸密切相关。此外，这也说明，空洞尺寸越小，使用MASW方法探测难度越大。上述对比结果展示了瑞雷波多参数FWI方法对于重建浅地表空洞等小尺度低速异常体的有效性。因此，将进一步应用于实际案例，以验证其在解决实际问题中的可靠性。

图  8  合成测试中MASW方法反演结果

Figure  8.  Inversion results of the MASW method in the synthetic test

下载: 全尺寸图片 幻灯片

3   实际案例

3.1   数据采集与预处理

案例研究区位于陕西省宝鸡市，探测目标为一处人工空洞(图9a)。铺设的接收排列如图9b所示，在空洞顶上方约4.5 m的高度布设一条线性测线阵列，包含92个4.5 Hz垂直分量检波器，道间距为1 m，总接收阵列长度91 m，同时以2 m的炮间距共设置46个炮点 (图9c)。空洞宽3.8 m、高度2.6 m，空洞中心位于第45道下方。针对目标位置，抽取了第16—第30炮的第33—第57道数据进行反演(图9c)，因此，在后文中分别将第16炮与第33道称为第1炮与第1道。采用10.8 kg大锤作为人工震源，每个炮点均重复激发并对获得的地震记录进行垂直叠加来提高地震记录的信噪比。图10展示了第3炮与第13炮归一化的原始地震记录，除了近源地震道受震源近场效应影响严重外，记录整体具有较高的信噪比。

图  9  野外实验工区与观测系统

Figure  9.  Work area and deployed observation system in the field experiment

下载: 全尺寸图片 幻灯片

图  10  野外实验中采集的原始记录数据

Figure  10.  Raw recorded data in the field experiment

下载: 全尺寸图片 幻灯片

在FWI反演前，对单炮记录进行预处理以保证反演结果的可靠性。首先切除早于初至信号到达的噪声信号，并设置有效信号时间窗口以衰减其后的信号能量。由频谱分析的振幅谱可以发现(图11a)，数据能量集中于10~60 Hz，因此对各道信号进行5~65 Hz的带通滤波以消除来自地层深部天然源产生的低频噪声以及来自周边人为活动产生的高频噪声。随后各道数据与褶积因子$ \sqrt{{t}^{-1}} $进行褶积运算实现3D波场到2D波场的转换，同时各远偏道(偏移距$r$大于11 m)数据乘以$ r\sqrt{2{t}^{-1}} $进行振幅补偿来校正几何扩散的影响^[42]。最后各道信号均延迟0.05 s以避免估计震源子波受到非因果效应影响。预处理后的地震数据(图12)与原始的相比，信噪比得到进一步提升，有效震相的信号能量明显增强。

图  11  野外实验中记录数据的频谱与频散谱

Figure  11.  Spectrum and dispersion spectrum images of recorded data in the field experiment

下载: 全尺寸图片 幻灯片

图  12  野外实验中预处理后的记录数据

Figure  12.  Recorded data after preprocessing in the field experiment

下载: 全尺寸图片 幻灯片

模型尺寸及网格步长与合成测试中的相同，并同样基于第1炮记录的面波频散谱(图11b)建立初始模型。初始横波速度由地表的150 m/s以4 m/s的步长线性增加至模型底部的226 m/s(图13a)。初始纵波速度(图13b)与初始密度(图13c)仍与横波模型保持正相关，初始纵波速度为初始横波速度的2倍，初始密度由地表的1.2 g/cm³以0.032 g/cm³的步长线性增加至模型底部的1.808 g/cm³ (图13d)。对数据采集中人工激发的子波进行合理估计是实际数据处理中的一项关键步骤。这里笔者采用校正滤波方法^[36]，即通过将观测的地震记录与基于特定参数模型使用已知子波所模拟的地震记录进行反卷积而构建一个校正滤波器，之后采用该滤波器与模拟地震记录进行卷积运算从而估计出震源子波。由于估计的震源子波是一个使模拟数据与观测数据残差最小化的最优结果，为了减弱其受到不准确的特定参数模型的影响，在反演的每次迭代进程中均对其进行动态估计。图14展示了基于初始模型估计的所有震源子波，它们具有高度一致性，这对于保证反演结果的可靠度具有重要意义。

图  13  野外实验中的初始模型与FWI最终反演结果

Figure  13.  Initial models and final FWI results in the field experiment

下载: 全尺寸图片 幻灯片

图  14  野外实验中估计的子波

Figure  14.  Estimated source wavelets in the field experiment

下载: 全尺寸图片 幻灯片

3.2   反演结果与分析

反演过程中的约束条件、超参数的设置以及反演尺度分割与合成测试中的均相同。两个尺度的反演进行了22次迭代，归一化目标函数值收敛至约0.55(图15)，考虑到数据中存在噪声、地层衰减以及参数串扰等影响，这个误差值是可以接受的。反演结果如图13d、图13e与图13f所示，相较于初始模型，反演模型展现出较强的横向非均匀性。反演模型在深度2 m处可大致分为两层，地表处存在未被预处理算子完全压制的局部假象，这种假象可以采取平滑滤波予以缓解。横波速度模型得到了较高程度的更新(图13d)，在模型中部包含了一个带有模糊边界的低速体(图13d黑色虚线标注的区域)，位于第11道—第15道下方，顶板与底板埋深分别约为4.5 m与7.5 m，最大宽度与高度分别约为4 m与3 m，与人工空洞的实际位置与尺寸十分相近。更新后的密度模型形态与横波模型的较为相似(图13f)，但中部的低速体位置略为偏右，这与面波对密度的敏感度较低且受横波的串扰影响有关。波场中面波能量极强导致体波相对被压制，致使纵波模型的更新有限(图13e)，未完全揭露出与横波模型、密度模型中相同的低速异常体。同样由于在模型边界处波场照明严重不足，导致反演结果存在较强的边界效应与“拖尾”现象^[19]。图16展示了最终反演结果对应的预测数据与观测数据之间的拟合度(波形隔道抽稀显示)。由于浅地表地层介质的复杂性、地层衰减、震源近场效应、背景噪声等因素导致某些地震道仍然存在一定程度的拟合残差，然而二者的波形拟合度整体性较好(图16a、图16b)，从而保证了最终反演结果的可靠性。

图  15  野外实验中的目标函数下降曲线

Figure  15.  The descending curve of the objective function in the field experiment

下载: 全尺寸图片 幻灯片

图  16  野外实验中最终预测数据与记录数据的波形对比

Figure  16.  Comparison of observed waveforms with waveforms predicted through the final FWI results in the field experiment

下载: 全尺寸图片 幻灯片

以上这些结果有力地说明，提出的方法具有对小尺度低速异常体(如空洞)的探测能力。另外，值得注意的是，整套反演流程在常规笔记本电脑上采用单核运算耗时不足1 h，展示了该方法计算的高效性。

4   结 论

a. 开发了一套完整的二维瑞雷波多参数全波形反演方法与流程，包含FWI目标函数的拟Hessian矩阵预处理算子、波场正演与实测数据维数转换、校正滤波方法动态估计震源子波等关键技术与处理流程。模型和实际数据检验了该方法与流程的可行性。

b. 经测试含有空洞的模型数据和实际数据，验证了该方法具有精细刻画浅地表小尺度诸如空洞的低速异常体的能力。实测数据中的横波速度模型结果清晰展示出一个4 m×3 m的人工空洞，与其实际位置与尺寸十分接近。

c. 二维瑞雷波多参数FWI方法获得的横波、纵波速度模型、密度模型具有总体一致性，横波速度模型最为准确地刻画出空洞的位置与尺寸。实测数据反演结果的评价应以横波速度模型为主，辅助参考纵波速度模型与密度模型。