0 引言

在浅层介质结构探测中，瑞雷波的能量约占波场总能量的三分之二^[1-2]，因而地震记录中最突出的波形是直达瑞利波及其频散波组。瑞利波是由P波和SV波干涉形成，其频散波组中包含丰富的S波速度信息，利用瑞利波反演可以得到地下结构的S波速度。全波形反演技术因具有精度高、可进行多参数建模的特点，越来越受到勘探地球物理学家的重视。然而，全波形反演存在以下不足：其大多数应用集中于体波，面波常被视为噪声而去除；常规全波形反演通常缺乏低频信息，多数研究以反演P波速度为主。而瑞利波的低频特性恰好可以不受低频信息缺失的限制，且瑞利波相速度与S波速度密切相关，因此，利用瑞利波波形反演浅层S波速度结构较常规全波形反演具有一定优势。

早期瑞利波波形反演主要用于研究全球或区域大尺度的地球结构^[3-6]，由于研究中使用的面波波长大于浅层结构的尺寸，使该方法在浅层分辨率方面存在不足。随着研究的进展，出现了适用于浅层构造成像的面波全波形反演方法。Zeng等^[7]证明了瑞利波波形反演用于探测浅层速度结构的可行性；Tran等^[8]基于高斯—牛顿法在时域进行弹性全波形反演，成功重构了20m范围内的浅层S波速度结构；Groos等^[9-10]研究了介质黏弹性品质因子对反演结果的影响，随后实现了浅层瑞利波全波形反演，并通过模型试验分析其稳定性；Borisov等^[11]利用基于包络的目标函数反演浅层大尺度、不均匀结构，证明了基于包络的目标函数能有效减少周期跳跃问题。虽然瑞利波波形反演研究已经取得了很大进展，但仍存在震源子波估计问题，需要进一步研究。

由于野外资料的震源子波未知，在反演中使用不准确的地震子波将不可避免地影响波形反演的结果^[12]，因此，精确的子波估计是波形反演成功的关键。处理地震子波的方法主要有两类：①在反演的迭代过程中进行震源子波校正滤波^{[9, 13-16]}，即在正演模拟时首先采用预先选取的地震子波(如Ricker子波)获得合成数据(通常为共炮点道集)，然后利用最优校正滤波器对合成数据进行震源子波校正滤波。基于线性阻尼最小二乘理论，通过使观测数据和波形校正滤波后的合成数据的傅里叶系数的残差最小，在频域构造校正滤波器。虽然这种子波反演(或校正)方法可以获得较好的反演结果，但反演过程需要多次迭代，增加了运行时间。Park等^[17]采用阻尼函数对共中心点(CMP)道集进行子波校正，提高了计算效率。②将零炮检距地震道作为经验子波或采用虚拟震源法，通过特征道褶积消除震源子波的影响^[18-21]。

Schuster^[22]提出一种反射波加窗提取地震子波的方法，可快速获得近似震源子波，但窗口大小是人为确定的，窗口选取不合适会影响子波提取精度。为了更准确估计震源子波，本文采用自适应子波提取与子波校正联合的方法。首先，对观测数据进行自适应窗口时移叠加提取震源子波，子波提取时可自动确定窗口大小；其次，将提取的子波作为初始模型进行子波反演校正；再次，建立一组断层模型并采用交错网格有限差分法获得其对应的合成记录，进行波形反演试验，采用不同的子波处理方法得到各自的反演结果，对这些反演结果进行对比分析以验证各方法的应用效果；最后，通过一组实际数据进一步验证联合窗口自适应子波提取和子波校正波形反演方法的实际应用效果。

1 方法原理 1.1 瑞利波波形反演基本理论

全波形反演的基本原理是寻找一个能够对观测数据进行合理解释的“最优模型”，该模型能够拟合特定地震波场(直达、反射、折射、面波等)的相位和振幅。

设炮点x_s处激发、检波点x_r处接收到的t时刻地震波的第i个位移分量为u_i(x_s, x_r, t)。如果已知地下模型参数分布，则可用有限差分方法求解波动方程，得到的模型正演数据记为u^mod。将理论模型数据与实际观测数据u^obs进行比较，得到数据残差为δu=u^mod-u^obs。如果δu很小，则表明该模型能够很好地解释观测数据。波场的能量残差E可用残差向量的L₂范数表示

$ \boldsymbol{E}=|L|_{2}=\frac{1}{2} \delta \boldsymbol{u}^{\mathrm{T}} \delta \boldsymbol{u} $

(1)

瑞利波波形反演的基本步骤为：①给定一个初始模型，然后通过波场数值模拟计算初始模型或迭代修正模型的瑞利波场；②计算观测波场与模拟波场之间的波场残差即目标函数；③进行反演迭代，直到获得满足迭代精度要求的地下速度结构成像。

在进行初始模型或迭代修正模型的波场数值模拟计算时，震源子波的估计或校正是波形反演中的一个关键步骤，需给定与观测波场一致的震源子波，否则将会对波形反演结果带来负面影响。

1.2 窗口自适应子波提取

通过自适应设定一个窗口，在窗口内对瑞利波记录进行震源子波提取(图 1)。具体步骤如下。

(1) 统计每一道峰值大小以及对应的时刻，根据最小时刻和最大时刻之间的时间宽度各向两端扩大一定比例(如5%)，确定自适应窗口的时间长度T；再根据最大峰值所在的道号以及第一个峰值小于某个阈值(如最大峰值的20%)所在的道号确定自适应窗口的位置长度R(R为道数)，可以自适应地获得一个大小为T×R的局部窗口。

(2) 在局部窗口内通过第一道与其他道的互相关计算时间偏移，通过时间偏移对局部窗口中的所有道进行平移。

(3) 对窗口内所有道进行叠加以提高子波的信噪比。

(4) 根据需要调整子波长度。

1.3 子波校正滤波

波形反演中观测数据与合成数据的残差主要由真实模型与当前计算模型的差异造成，也有一些是由于观测数据对应的真实子波与合成记录使用的假定子波之间存在差异造成的。为了消除子波带来的影响，Groos等^[9]提出一种校正滤波方法处理波形反演中的子波问题，采用如下基本思路。

首先对当前模型采用同一个假定子波得到各炮合成记录，然后将合成记录与观测数据进行反褶积，得到子波校正滤波器，利用这些滤波器对前面正演步骤中使用的假定子波进行滤波，得到每一炮对应的最优子波。在多尺度波形反演过程中，每个频段只进行一次震源子波校正滤波。子波校正滤波器是在频率域构建的，利用线性阻尼最小二乘法使校正滤波后的观测数据与合成数据的傅里叶系数的残差最小，进而构造一个最优校正滤波器。

采用目标函数

$ \begin{aligned} F\left(\tilde{c}_{l} ; \varepsilon\right)=& \sum\limits_{l=0}^{N-1} \sum\limits_{k=1}^{M} f_{k}^{2}\left|\widetilde{d}_{l k}-\tilde{c}_{l} \widetilde{g}_{l k}\right|^{2}+\\ & M \overline{{E}} \varepsilon^{2} \sum\limits_{l=0}^{N-1}\left|\tilde{c}_{l}\right|^{2}=\chi^{2}+\psi^{2} \end{aligned} $

(2)

式中：$\tilde{d}_{l k} $、$ \widetilde{g}_{l k}$和$ \tilde{c}_{l} $分别为炮检距r_k处观测数据、合成数据以及待求滤波器的傅里叶系数；比例系数f_k是与炮检距相关的权重因子，用于补偿信号在不同炮检距处的幅值，消除炮检距的影响；ε是阻尼因子，对合成数据平均频谱能量进行加权；l为傅里叶系数个数，l=0, 1, …, N-1；k为检波器个数，k=1, 2, …, M。上式在任意频点处都成立。

$ \tilde{c}_{l} $是一个复数，包含实部和虚部

$ {\tilde c_l} = \tilde c_l^\prime + {\rm{i}}\tilde c_l^{\prime \prime } $

(3)

式(2)中

$ \chi^{2}=\sum\limits_{l=0}^{N-1} \sum\limits_{k=1}^{M} f_{k}^{2}\left|\widetilde{d}_{l k}-\tilde{c}_{l} \widetilde{g}_{l k}\right|^{2} $

(4)

表示用f_k加权后的数据残差。

$\psi^{2}=M \bar{E} \varepsilon^{2} \sum\limits_{l=0}^{N-1}\left|\tilde{c}_{l}\right|^{2} $

(5)

表示对数据进行归一化处理。其中

$ \bar{E}=\frac{1}{M N} \sum\limits_{l=0}^{N-1} \sum\limits_{k=1}^{M} f_{k}^{2}\left|\widetilde{g}_{l k}\right|^{2} $

(6)

通过对$\tilde{c}_{l} $的实部和虚部的目标函数最小化，即

$ \left\{\begin{array}{l} \frac{\partial F}{\partial \tilde{c}_{l}^{\prime}}=0 \\ \frac{\partial F}{\partial \tilde{c}_{l}^{\prime \prime}}=0 \end{array}\right. $

(7)

可得滤波器的傅里叶系数

$\tilde{c}_{l}=\frac{\sum\limits_{k=1}^{M} M f_{k}^{2} \widetilde{g}_{l k}^{*} \widetilde{d}_{l k}}{M \bar{E} \varepsilon^{2}+\sum\limits_{k=1}^{M} M f_{k}^{2}\left|\widetilde{g}_{l k}\right|^{2}} $

(8)

式中“*”号表示复共轭。

2 模型试算

为了验证子波提取与校正方法的有效性，建立如图 2所示断层模型。模型尺寸为40m×12m，断层上、下盘埋深分别为8、4m。上层介质P、S波速度分别为800、200m/s，下层介质P、S波速度分别为1200、400m/s。上、下层介质的密度均为2000kg/m³。数值模拟采用交错网格有限差分法。网格单元尺寸为0.2m×0.2m，时间步长为0.05ms。炮点位置分别为0、10、20、30、40m，炮间距为10m。每炮记录36道接收，道间距为1m，记录长度为0.6s，采样率为0.5ms。每炮的炮点、接收点位置固定。

采用图 3所示的S波速度递增模型作为初始模型，分别以真实震源子波、错误子波、错误子波联和子波校正、子波提取得到的子波、提取子波联和子波校正等5种子波测试瑞利波波形反演效果。为验证本文方法不受真实子波的影响，对于前三种情况，模型对应的合成记录均采用30Hz正弦衰减子波。对于模型试算中的其他情况，模型对应的合成记录则采用30Hz的Ricker子波。

第一种情况，当实际资料的震源子波已知时，正、反演过程均采用同一个真实子波(30Hz正弦衰减子波)，反演结果如图 4所示。图 4a为经过79次迭代后的断层模型速度结构反演结果。可以看出，断层界面在地下4、8m处有清晰的速度分层，基本可以准确重建断层模型。图 4b为真实模型、反演使用的初始模型和反演结果对应的第3炮记录的波形对比图(每隔2道显示1道)。为便于分析，图 4c中只显示了第10道。从图 4c可以看出，真实模型对应的波形与初始模型对应的波形差异较大，但与反演结果基本吻合。由此可见，当震源已知时，瑞利波波形反演结果是可靠的。

第二种情况，采用错误的地震子波进行波形反演，即真实子波为30Hz的正弦衰减子波，反演采用30Hz的Ricker子波。图 5为经过92次迭代得到的反演结果。可以看出，断层位置明显错误，分层边界不清晰，上、下层速度误差较大。由此可见，使用错误震源子波得到的反演结果明显不可靠。因此，为提高反演精度，有必要在反演之前先进行震源子波提取。

第三种情况，反演采用错误的震源子波，但同时进行了子波校正，亦是真实子波为30Hz的正弦衰减子波，反演采用30Hz的Ricker子波。图 6为经90次迭代得到的反演结果。从反演结果可以识别出断层构造的基本轮廓。与未进行子波校正滤波的反演结果(图 5)相比，反演精度得到了一定程度的提高，但仍存在断层位置错误、分层边界不明显、速度反演结果误差大等不足。由此可见，在真实震源子波未知的情况下，基于震源子波校正的波形反演方法可以在一定程度上提高反演精度。然而，对于更为复杂的横向非均匀模型，精度改进有限，无法满足高精度反演的要求。

第四种情况，真实子波为30Hz的Ricker子波，反演采用从不同炮记录中提取的子波，提取结果如图 7所示，其中红色曲线为真实的Ricker子波，其他颜色的曲线分别对应不同炮记录的子波提取结果。可以看出，虽然每一炮提取的子波有所不同，但都与真实子波比较接近，可以直接用于波形反演。

图 8为仅使用提取的子波进行反演对应的断层模型反演结果。图 8a为经97次迭代后的S波速度反演结果，反演使用的初始模型与图 3相同。可以看出，断层成像准确，速度分层较为清晰，上、下各层速度接近真实值。但在水平方向25~40m范围内，下层速度约为300m/s，与真实值400m/s有较大差异。这可能是提取子波与真实子波之间的误差引起的。图 8b为真实模型、初始模型和反演结果对应的第3炮记录的波形对比(每隔2道显示1道)，图 8c为图 8b中的第28道。从图中可以看出，真实模型对应的波形与初始模型相差较大，但与反演结果基本吻合。说明仅利用提取的子波进行波形反演，反演结果基本准确。

第五种情况，为了进一步提高反演精度，在子波提取的基础上加入子波校正，进一步对地震子波进行约束。联合使用提取的子波和子波校正进行反演，图 9为经过89次迭代后的反演结果。从反演结果可以看出，断层模型的成像效果较好，速度分层清晰(图 9a)。值得注意的是，8m以下的断层上盘速度比仅用提取的子波反演的速度(图 8a)更准确(接近真实值400m/s)。虽然浅部速度误差较大，且存在局部速度异常，但经过子波校正后，反演精度明显提高。图 9b为真实模型、初始模型和反演结果对应的第3炮记录的波形对比(每隔2道显示1道)，图 9c为图 9b中的第34道。可以看出，真实模型对应的波形与初始模型相差较大，但与反演结果基本吻合。因此，联合使用提取的子波和子波校正可以进一步提高反演精度，特别是提高深层的速度反演精度。另外，反演迭代次数比仅使用提取的子波进行反演相比减少了8次，说明联合使用子波提取与校正的方法可以节省一定的计算时间。

本文采用多尺度反演策略，在给定的频带内对相应的观测数据和提取的子波从低频到高频逐带进行反演。在反演过程中，每次使用一个新的带宽更大的数据子集时，子波校正滤波器都会估计出一个最优的震源子波用于波形反演，直到数据带宽再次增加。图 10为第3炮提取的子波在相应频带宽度内经过子波校正后得到的地震子波。可以看出，校正后的子波随着频带宽度的变化而变化。当频带宽度为0~70Hz时，采用提取的子波联合子波校正得到的子波(蓝色虚线)与真实子波(红色实线)及提取的子波(蓝色实线)基本重合。图 11是每一炮提取的子波经子波校正滤波后得到的最终子波。可以看出，最终的子波与真实子波相似。

由图 4~图 9不同子波处理的断层模型反演测试结果表明，基于子波提取与校正的波形反演可有效提高反演精度，特别是可以提高深层速度结构的反演精度。

3 实际资料测试

选取了一个人工隧道进行实际资料瑞利波探测试验。Wang等^[23]和Smith等^[24]开展的类似研究证明，二维或三维弹性全波形反演方法在浅层隧道探测和成像方面具有广阔的应用前景。本次试验利用基于二维弹性全波形反演方法验证震源子波对隧道探测效果的影响。现场测试点如图 12所示，隧道宽约6m，顶板埋深约3m，底板埋深约6m。测线垂直于隧道走向，位于隧道上方的平台上。从西向东0~47m之间放置48个间隔1m的低频(4.5Hz)检波器。炮点位于12、18、24、30、36m处，炮间距为6m。采用15磅锤击震源，记录长度为512ms，采样率为0.5ms。隧道位于测线水平位置27~33m处。

反演计算区域为59m×25m，网格单元尺寸为0.25m×0.25m，采用0.05ms时间步长。波场计算采用了30Hz的正弦衰减子波，合成记录长度为512ms，采样率为0.5ms。每炮记录共48道，道间距为1m，炮间距为6m。采用自由表面边界和完全匹配层吸收边界(厚度为5m)进行计算。除去吸收边界厚度，与实际数据匹配的模型尺寸为49m×20m。

图 13为现场实际数据的5个单炮记录，从图 13a和图 13b可以看出，由于隧道空腔的存在，第27~第33道之间(红色圆圈为空洞位置)同相轴的能量及相位发生了明显的变化。由于野外采集的面波数据具有较高的信噪比，对数据仅进行了5~70Hz的带通滤波预处理。图 14为初始模型(图 14a)仅使用提取的子波(图 14b)和联合使用子波提取与校正经过20次迭代得到的反演结果(图 14c)。从图 14b、图 14c可以看出，在27~33m处(红色圆圈为空洞实际位置)存在明显的低速异常，其位置与隧道实际位置几乎一致。隧道顶板埋深可以很好地成像，但底部界面反演深度略深于实际埋深，总体成像效果令人满意。图 15为第4炮的观测波形与预测波形对比。图 16为相应的反演拟合误差收敛曲线对比。对比可以看出，仅使用子波提取方法得到的隧道宽度反演误差比联合使用子波提取与校正方法的误差大得多。实测数据反演结果表明，联合子波提取与校正可以进一步提高瑞利波波形反演精度。

4 结论

震源子波的选取是影响地震波场正演精度的重要因素之一，进而影响瑞利波全波形反演结果的精度。理论模型的测试结果表明如下认识。

(1) 采用正确子波对应的反演结果最好，采用错误子波对应的反演结果最差。如果震源子波未知，直接进行瑞利波波形反演是不可靠的。

(2) 窗口自适应子波提取方法可以从观测记录中提取子波，并将提取的子波直接用于波形反演，能够获得较好的反演效果。

(3) 在子波提取方法基础上，对提取的子波进行子波校正滤波，进一步约束震源子波。试验结果表明，联合子波提取与子波校正能够有效地提高深部速度的反演精度，使反演速度更接近实际值。