0 引言

走时是地震波到达观测点的时间，是地震信号中最可靠的间接观测信息.在射线理论中，地震波走时被定义为射线走时，与地下慢度分布成线性关系.由于射线走时所反映的是地震波在大尺度非均匀介质中传播的运动学特征，所以通常被用于反演宏观速度或背景速度建模.通常而言，初至及折射波走时多用于近地表速度建模，而反射波的运动学信息(走时及斜率)则广泛应用于中深层速度建模.利用反射走时实现速度反演，主要需要回答两个问题：(1)如何测量反射波时差；(2)如何构造反射波路径并实现反射时差的反投影.在理论上，这两个问题既可以在数据域得到解决，也可以在成像域得到解决.因此，按照实现方式，反射波速度反演可分为成像域和数据域两类反演方法.

典型的成像域反演方法如基于射线理论的成像道集层析(Ray-based tomographic migration velocity analysis；Stork，1992；Nemeth，1995；Woodward et al., 1998)，用共成像点道集(common image gather，CIG)的剩余深度差(residual moveout，RMO)近似计算反射时差，避免了反射时差的显式测量，并用射线理论计算反射波路径.为突破射线理论的局限性，一些学者提出用单程波动方程(Xie and Yang, 2008；Bevc et al., 2008)代替射线理论，引入了波路径层析方法(Fliedner and Bevc, 2008)，可以适应强速度差异和多路径问题并考虑波传播的有限频带效应.Zhang等(2011, 2012)给出了双程波动方程构造的反射走时敏感度核函数，并用于成像域反射走时反演.然而，除射线理论的高频近似假设外，成像道集层析还引入了层析的线性近似、走时差计算和估计“真深度”时的常速假设等多种近似.此外，基于波动方程的偏移速度分析(Wave-equation migration velocity analysis，WEMVA)也是常用的成像域反演方法，通过极大化成像道集的能量(如叠加能量最大化目标函数)；(Soubaras and Gratacos, 2007)、或惩罚成像道集中像的差异性或聚焦性等(如差异相似优化，differential semblance optimization，DSO)；(Symes and Versteeg, 1993；Shen，2005；Symes，2008；Shen and Symes, 2008；Liu et al., 2014)，将“像差”向反射波路径反投影，实现背景速度反演.WEMVA用波动方程作为正演算子，克服了射线理论的高频近似假设；通过用波动方程产生成像道集，可以实现自动化反演(无需拾取地下反射点位置和倾角，以及扫描成像道集的RMO)；并且对成像道集的RMO没有任何假设，因此适用条件更广.然而，WEMVA方法的计算量巨大，且存在地下照明不均衡以及泛函梯度假象等问题(Fei and Williamson, 2010；Yang et al., 2013；Zhang and Shan, 2013；Shen and Symes, 2015)，影响了反演的收敛速度甚至导致反演不收敛.Luo等(2016)提出了全走时反演方法(Full traveltime inversion, FTI)，通过假定速度的变化仅引起地震波走时变化，构造了成像域的走时反演方法.然而，该方法仍然需要计算并存储成像道集，对计算和存储都有较高的需求.

对于数据域反演方法，通常需要测量反射波时差.在立体层析(Stereotomography，Billette and Lambaré，1998 ；Lambaré，2008 )中，通过关联数据(坐标、射线参数、双程走时)和模型(地下反射点、斜率、张角)参数，实现了反射时差的自动测量.同时，利用了炮、检处射线传播的局部方向约束了射线路径，在一定程度上避免了射线多路径现象.然而，由于正问题采用射线理论，仍然受制于射线理论固有的缺陷.此外，由于层析矩阵不仅包含模型对速度的敏感度矩阵，还引入了对坐标和角度等的敏感度矩阵，因此该方法还受多参数反演问题的制约，需要引入模型正则化等模型先验信息来帮助反演收敛.由于反射地震数据的频带宽度是有限的，所以若无关于模型的先验信息很难直接在数据域准确测量观测数据和模拟数据中来自地下同一反射界面的同相轴的时差.Ma和Hale(2013)用DTW(dynamic time warping)算法估计两个信号的时变时差，并构造了反射走时反演方法.然而DTW算法中关于目标函数(距离)的定义依赖地震波振幅，因而无法消除振幅对时差估计的影响.Feng和Wang(2015)基于波动方程反偏移和互相关时差测量方式，导出了反射波走时敏感度核函数表达，但没有给出反射波时差的自动测量方法.李辉等(2017)基于地震数据的特征高斯波包分解，实现了反射时差的测量并给出了高斯波包框架下的反射走时反演方法.然而，对于复杂地震数据的特征高斯波包分解仍然难以实现自动化.另一方面，基于聚焦准则的反演方法回避了时差的显式测量，试图用观测数据和模拟数据的互相关函数的某种加权范数直接构造目标泛函(van Leeuwen and Mulder, 2010)，利用泛函的聚焦特性反演背景速度.但Baek等(2013)指出，此类泛函存在梯度符号问题，仅对单个反射层成立.

上述回顾表明，对于反射地震信号，如何在数据域利用纯走时信息反演背景速度仍然是一个颇具挑战性的工作，而其中的核心问题则是如何不受振幅干扰并准确估计反射波走时及时差.因此，本文从地震数据的特征表达出发，通过构造立体数据空间，提出了一种自动且稳健的反射波时差测量方法，实现了基于特征波场分解的反射走时反演，用于背景速度的自动反演.

1 理论 1.1 基于特征波场分解的反射时差测量方法

王华忠等(2015)给出了如下叠前地震数据的特征波场分解方法(Characteristic wavefield decomposition，CWD)，

(1)

其中，d(t, s, r)为地震记录，s=(s_x, s_z)和r=(r_x, r_z)分别为炮点和检波器的坐标，s₀和r₀分别为参考炮点和检波器的坐标，p_s=(p_sx, p_sz)和p_r=(p_rx, p_rz)分别为局部平面波的入射和出射射线参数的水平分量，W_s和W_r分别为炮点和检波器端局部平面波分解的空间窗长度(使得地震波满足局部线性特征)，c(t, p_s, p_r; s₀, r₀)为特征波分解后得到的具有特定传播方向的时空局部化的子波.

在实际应用中，(1)式需要预先估计局部平面波的射线参数才能实现特征波场分解(王华忠等，2015).为了同时实现特征波分解及传播方向估计，本文首先构造了如下特征波场分解的伴随变换：

(2)

利用射线参数空间的稀疏性作为约束条件，可以在稀疏反演框架下将特征波场分解问题(1)转化为如下模型参数估计问题P₀:

(3)

其中，L₀范数表示向量c中非零元素的个数，d=d(t, s, r; s₀, r₀)，c=c(t, p_s, p_r; s₀, r₀)，算子T为(2)式中积分过程的算子形式，σ²为噪声能量阈值.

(3) 式可以用匹配追踪类的方法求解(Mallat and Zhang, 1993)，从而同时得到局部平面波子波及入射和出射射线参数.经过特征波场分解，可以直接将叠前地震数据投影到立体空间(炮/检点坐标，入射/出射射线参数及走时)中.对于带限地震信号的到达时，有多种定义方式，如起跳时刻，最大振幅时刻等.考虑到实际地震信号中总是存在一定的噪声，用起跳时刻定义的走时在实际数据中难以直接应用.因此，本文采用立体层析方法(Billette et al., 2003)中走时拾取的策略，用地震道的包络极值定义地震波走时：

(4)

其中，i, j分别是参考炮、检点的下标，k为特征波下标，E为包络算子，定义为

(5)

为实信号c(t)的希尔伯特变换， *代表褶积.

显然，通过特征波场分解(3)式及走时定义(4)式，可以自动生成如下立体数据空间：

(6)

值得注意的是，本文中的特征波场分解方法并不追求观测数据的精确重构，仅仅提取了观测数据中符合高维局部平面波特征的波现象，作为观测数据运动学信息的一种抽象与特征表达，进而建立立体数据空间，用于后续的速度反演.

为了测量反射波时差，以二维模型为例，给出如下策略.如图 1所示，可以在初始模型中从炮点和检波点处分别以入射和出射射线参数向地下传播两根射线.若这两根射线相交于Q点，此时有如下走时恒等式：

(7)

其中，T_{i, j, k}^obs为(4)式定义的观测数据走时，T^cal(s_i, q, p_s^k; m₀)和T^cal(r_j, q, p_r^k; m₀)分别为初始模型m₀中从炮点S到交点Q以及检波点R到交点Q的两根射线的单程走时，Δt_{i, j, k}为反射时差.P点为真实速度模型中反射点的位置.

此时，反射时差可以直接写为

(8)

显然，利用本文采用的策略，可以直接测量立体数据空间中每一个立体数据的反射时差，既避免了周期跳跃以及串层等可能性，又消除了振幅对于时差测量的影响.

1.2 基于特征波场分解的反射走时反演方法(Characteristic wavefield reflection traveltime inversion, CWRTI)

利用(8)式的走时残差，可以构造如下走时目标泛函：

(9)

其中，m为模型参数化方式.本文中采用慢度作为模型参数化方式，m(x)=v^-1(x).

目标泛函的梯度为

(10)

根据式(8)，走时残差对模型的Fréchet导数可以表示为

(11)

显然，是关于模型参数的非线性函数.根据走时正问题的具体表达形式，可以推导的线性近似.如射线理论中，为从x到y的射线路径；而在有限频层析理论中，为从x到y的有限频敏感度核函数(Xie and Yang, 2008)；对于波动方程走时反演，可以在一阶Born近似下用波场对模型的Fréchet导数表达(Luo and Schuster, 1991；Feng and Wang, 2015).考虑到计算效率问题，本文采用射线近似计算公式(11).

慢度模型可以用梯度导引的算法更新：

(12)

其中，α^k为步长，g_l(m^k)为梯度▽J(m^k)的低波数部分(背景部分).

通常情况下，无论采用何种形式的近似，泛函梯度中会存在各种假象，影响反演的收敛.在全波形反演或反射波形反演中，泛函梯度的低波数部分可以在成像前或成像后提取，如利用波场传播方向(Poynting矢量)构造梯度的低波数部分(Tang et al., 2013)；频率-波数域上下左右行波分解成像条件构造梯度(Qin et al., 2013；Wang et al., 2013；Wu and Alkhalifah, 2014)；时移成像条件+保留大角度成像结果(Alkhalifah，2014)；波数域中高-低波数的尺度分解(Almomin and Biondo, 2013)；梯度正交化方法(Albertin et al., 2013).

本文将泛函梯度视为图像，用全变差(Total variation，TV)正则化方法提取图像纹理g_h(m^k)(即梯度的高波数部分)：

(13)

其中，g(x)=▽J(m^k)为第k次迭代时目标泛函的梯度，x为地下空间变量，g_h(x)为梯度中变差较大的边界纹理部分(待求解)，算子其中L_m和L_n分别是水平和垂直方向的一阶导数.

求解(13)可得梯度的高波数部分g_h，将g_h从原始梯度g中减去，可得梯度的低波数部分：

(14)

求解(13)式的优点在于：梯度计算与梯度预处理相互独立.因此，泛函梯度可以用射线近似计算，也可以用波动方程线性化近似计算.完成梯度计算之后，再进行图像稀疏提升.与其他成像前梯度预处理方法相比，基于图像去噪理论的梯度预处理方法更加高效.

2 数值试验 2.1 二维合成数据测试

首先，用二维合成数据测试CWRTI方法.速度模型取自Sigsbee速度模型的一部分(X=[0~6700]m, Z=[2300~6100]m)，如图 2a所示.速度模型以沉积岩为主，并有数组断层发育，无盐丘侵入.X和Z方向的采样点数分别为801和601，网格间隔为7.5 m.观测系统设计为拖缆观测，最小偏移距300 m，最大偏移距3900 m，道间距为15 m.第一炮在X=0 m激发，炮间距30 m，共201炮.震源子波采用20 Hz主频的Ricker子波.记录时间5 s，采样间隔1 ms.

初始模型采用1.5 km·s^-1的常速模型.在反演中，CWRTI梯度用射线近似计算，用于更新大尺度的背景速度.其中，(13)式的正则化因子取λ=0.2.经过100次迭代后，相对走时残差降低至1%，反演停止.图 2b为100次迭代后输出的速度模型，与真实模型的光滑背景部分吻合程度较好.为进一步对比反演结果，将不同地表位置(X=1, 3, 5, 7 km)的纵向速度抽取，如图 3所示，其中，灰色虚线，灰色实线和黑色曲线分别为初始速度，CWRTI反演速度和真实速度.显然，CWRTI反演对速度与真实速度的变化趋势较为吻合，较好的恢复了中-大尺度的背景速度结构.为对比成像结果和道集，我们利用高斯束偏移(GBM)方法并输出角度域共成像点道集(ADCIG).图 4(a，b)分别是反演模型和真实模型中输出的角度道集.可以看出，经过走时反演之后，角度道集基本拉平.图 5(a，b)分别是反演模型和真实模型的GBM偏移结果.对于速度模型、成像结果和角度道集可以看出，由于深度速度横向变化较为剧烈(深度大于3500 m)且缺乏足够的照明角度，导致反演结果欠佳.同时，最深部的水平反射界面的成像深度及形态与真实模型存在一定偏差.这种小尺度的速度结构变化需要用后续的全波形反演以进一步提高速度模型的分辨率.

2.2 二维实际数据测试

用一条二维陆上实际数据进一步测试本文方法.该数据共450炮，道长1125，采样间隔为4 ms.最小和最大偏移距分别为100 m和3200 m.炮间距50 m，道间距25 m.初始模型采用1.8 km·s^-1的常速模型.经过140次迭代后，相对走时残差降低至3%.反演的速度模型如图 6所示.分别用初始模型和反演模型进行高斯束偏移，输出的角度道集如图 7(a，b)所示.由于初始模型为常速模型，导致成像道集中存在较大的剩余曲率.相比之下，经过走时反演后成像道集基本拉平.图 8(a, b)中由于初始速度偏低，导致成像剖面画弧明显.而反演速度得到的偏移剖面的断层清晰，绕射波收敛干净，但局部区域(X=9~12 km)的成像质量较差.

3 讨论

本文重点讨论的是如何在数据域进行利用反射走时信息进行速度反演，因而没有给出特征波分解((3)式)的求解算法.详细的求解算法、特征波分解的影响因素分析及参数选择等问题将另行讨论.

本文所采用的包络极值方法定义的地震波走时，通常而言应该大于地震波的真实到达时.只有对于无限频带的数据(脉冲地震子波)二者才相同.对于有限频带的地震数据，若地震子波未知，很难给出走时的准确定义.然而本文的定位是反演中-大尺度的背景速度，虽然将包络极值时间用于走时反演存在一定的系统误差(走时误差小于一个子波周期)，但在实际应用中仍然是可以接受的.此外，业界广泛应用的立体层析方法中，在数据域的走时拾取方式也采用了相同的策略(拾取归一化包络平方τ-p谱的极值位置).

本文提出的反射走时计算策略不同于立体层析.采用观测数据的斜率作为地表约束条件，从地表向地下传播两根单程射线，并以射线相交为判定条件从而计算反射波时差.而立体层析则是从模型空间中的地下散射点向地表传播两根单程射线.因此，即使初始模型一致，本文方法与立体层析得出的反射时差也不相同.

此外，在计算效率方面，由于本文反射时差计算策略只需要两次初值射线追踪，因此具有极高的计算效率.相比之下，成像域射线层析反演方法需要首先计算并存储地下共成像点道集(如角度或偏移距道集)，通过(自动)拾取RMO进而计算反射波时差，因此需要较大的计算量和存储量.

射线理论或波动方程线性化(Born近似或者Rytov近似)皆可用于本文目标泛函的梯度计算，但二者隐含的假设及成立条件却不相同.基于射线理论的泛函梯度，不存在弱散射假设，因而更适于反演大尺度的背景速度扰动，尤其是当初始速度与真实速度相差较远时.另一方面，若基于波动方程线性化(如Born近似等)计算梯度，理论上对速度扰动的分辨率更好，但计算量巨大.因此，综合考虑到反演精度及计算效率，本文在数值实验中均采用射线近似计算泛函梯度和对角Hessian矩阵.

当强速度反差存在时(如盐丘嵌入沉积岩中)，在光滑速度模型中进行射线追踪可能会有较大误差.因此，今后的研究中会进一步分析强速度反差对本文方法的影响.此外，本文提出的时差计算策略难以直接推广至三维问题.对于三维问题，更好的解决方案是将“射线相交”条件改为“走时相等”(即在地下某个深度时，两根单程射线的走时之和等于观测走时)，并测量两根射线的空间距离.通过将时差极小化准则改为距离最小化准则，克服三维情况下射线不共面的困难.

4 结论

本文中提出了一种基于特征波场分解的立体数据空间估计方法，实现了叠前地震数据运动学信息的特征表达.通过局部化的波传播算子及特征波聚焦成像条件，本文提出了一种新的反射波时差测量方法，既避免了周期跳跃问题以及串层等可能性，又消除了振幅因素对时差测量的影响.基于特征波反射时差，本文给出了一种自动化反射走时反演方法，并提出一种基于L₁范数约束的梯度正则化方法，用于反演背景速度模型，可以为后续的高精度反演方法(如波形反演)提供较好的初始速度模型.

理论分析和数值实验均表明，本文提出的特征波走时反演方法具有如下优点：(1)无需长偏移距观测数据或低频信息；(2)对初始模型依赖性低(文中的所有数值实验均用常速作为初始模型)；(3)计算效率高(不需要计算及输出共成像点道集)且无需显式走时拾取；(4)整个反演流程可以实现全自动化.

致谢  感谢同济大学数学系李辉博士在射线理论方面的讨论并提供射线追踪程序.感谢两位匿名审稿人提出的中肯建议.感谢中石油勘探开发研究院及西北分院、中海油研究院和湛江分公司、中石化物探技术研究院和胜利油田分公司对波现象与反演成像研究组(WPI)研究工作的资助与支持.