1 引言

地震波走时计算在层析反演、速度分析、正演模拟和Kirchhoff积分偏移中都是关键和基础性的工作.在实现方法上分为“现算型”和“存储型”两种^[1].“现算型”算法公式简单，便于走时反复计算，常用的是Dix公式，不适应横向变速情况^[2].“存储型”算法存储量大，实际计算中皆通过抽稀目标点进行走时求取，然后插值来加密网格，求解时计算量大，插值时会损失精度^[1].在理论方面，走时算法层出不穷：主要有射线追踪及波前重建法、有限差分法等^[3~8].射线追踪和波前重建法适于速度横向变化不大的情况，存在阴影区和焦散区的问题^[3~6]，计算效率较低，难以并行加速^[9]；有限差分类型的方法^{[7, 8]}突破了计算效率的限制^[9]，利用有限差分求解程函方程进行走时计算，在速度复杂地区不稳定^[10]，有很多学者对此进行了研究和改进^[10~18].上述算法均在不同程度上对波场进行近似和离散化，在保结构^{[19, 20]}前提下存在一定的缺陷；均属于“存储型”算法，用于Kirchhoff叠前深度偏移时，需存储海量走时表，是应用于实际资料处理时的一个瓶颈.

理想的走时计算方法应具有如下特征：在横向变速情况下能减少对波场近似和离散化引起的误差，保证计算精度；具有“现算型”方法的特点，适于走时反复计算.地震波走时李代数积分算法^[21~27]在理论上满足这个要求.该算法基于保结构思想^{[19, 20]}，将李群算法和拟微分算子理论^[28~32]引入地震波传播规律研究，涉及到线性群的运算，通过李代数进行表达^{[28, 29]}.具体实现上包括单平方根算子的李代数积分^{[1, 21, 22]}和指数映射^[30~32]，利用稳相点原理^{[33, 34]}得到走时的解析表达式，属于“现算型”算法范畴.因此相对于以往走时计算方法，李代数积分法在理论计算和实现方式上均有很大的改进和提高，目前时间积分算法^{[21, 22]}已在Kirchhoff叠前时间偏移中取得了很好的应用效果^{[1, 23~27]}.

本文通过对李代数时间积分进行深化调整，得到适于Kirchhoff叠前深度偏移的地震波走时计算方法.李代数深度积分的难点在于：复杂介质中单平方根算子的透镜项是空变的，与后续项存在交换运算，收敛很慢^{[1, 21]}；象征表达式中包含波数一次项，高次项的系数在求积时亦需进行修正.本文将坐标变换与李代数积分法相结合，得到了地震波走时的解析表达式.多项式的表达形式适于不同源点和接收点走时的反复计算，包含非对称项，保证了计算精度.实际计算中仅需记录多项式各项的系数，不需要存储海量的走时表，有望改进Kirchhoff叠前深度偏移的精度和效率.

2 方法原理

本文在李代数时间积分算法的基础上，通过时间空间域到频率波数域和向量场到指数流形以及时间域到深度域的正反变换，提出了一种新的地震波走时计算方法.该方法通过坐标变换，使深度域单平方根算子透镜项转化为常数，避免与高阶项的交换运算.在射线坐标系上求解单平方根算子主象征，利用Magnus算法^[32]和BCH（Bacher-Campbell_Hausdorff）公式^[35~37]，通过李代数积分和指数映射的计算，由稳相点原理得到地震波走时的解析表达式.

2.1 坐标变换

空间频率域的声波方程可以表示为

(1)

其中，ω为频率，v为速度，f为地震波场.通过如下坐标变换：

(2)

将波场转换到射线坐标系（x′，y′，T）下.其中，x′，y′为射线坐标，后文统一用x，y表示；T为初始射线的单程走时，这里选取成像射线；T′_z=∂T/∂z为走时在z方向上的导数.

将方程（1）转换到象征域，可得：

(3)

其中，σ代表象征运算符号，其定义见附录；s=σ（）为单平方根算子的象征^{[21, 22]}；=σ（f）为象征域的波场；#代表Witt积运算，其定义见附录，该算子在指数映射的象征表达中起关键作用^[1].

令s_main=σ_main（）代表单平方根算子的主象征^{[21, 22]}，由（3）式可得：

(4)

其中，分别表示走时T在x，y方向的导数.与时间积分单平方根算子的主象征相比^[1]，单平方根算子的主象征s_main包含初始成像走时的横向导数.因此，欲求s_main，先据程函方程递推求解初始成像射线的走时T，将其转换到射线坐标系上，利用其横向导数，即可由（4）式进行计算.

通过数值算例来说明初始成像射线走时及其深时转换情况.图 1a背景代表速度深度模型，图 1a、1b和1c等值线分别代表横向变速模型下成像射线的走时、坐标转换后走时和走时横向导数的分布情况.我们基于此来求取单平方根算子的主象征，进而计算李代数积分和指数映射的表达式.

设地面波场为f（x，y，0，ω），射线坐标系t₀处的波场为f（x，y，t₀，ω），若采用对波场递推的方式进行延拓可得：

(5)

根据李代数积分的定义^{[28, 29]}，的积分解满足，η为李代数积分运算符号，其定义见附录.可将（5）式转换为对算子进行积分的大步长延拓方式，其积分解为^[1]

(6)

其中，为大步长单程波算子.

在象征域进行计算，则（6）式转换为^[1]

(7)

令φ代表象征运算σ（）的复相位^{[21, 22]}，即exp[iφ（x，y，k_x，ky，ω）]，可知

(8)

其中，k_x，k_y分别为x，y方向的波数.从公式（5）~（8）可得^[1]：

(9)

若无横向变速，Witt积运算退化为简单的相乘，（9）式可得到简化，直接求得复相位φ，由稳相点原理即可得走时的解析表达式.横向变速情况下需增加高阶修正项.（9）式表明了从单程波算子到李代数积分及指数映射的对应关系，是整个算法的关键.

2.2 李代数积分

将（4）式展开，单平方根算子主象征s_main可以表示如下：

(10)

其中，s_i，j为波数各项的系数，下标代表k_x，k_y的幂次.

李代数积分的象征a亦是k_x，k_y的多项式，其系数可通过如下计算得到：根据单平方根算子的主象征s_main，利用Magnus算法进行级联式的递推，由低阶李代数积分的象征和其组成项的系数，通过交换运算和垂向积分，得到高阶修正及其组成项的系数.二维和三维情况下的具体表达式见附录，各项的详细推导参见文献[1]和^[22].

本算法中单平方根算子主象征s_main与时间积分^[1]中的表达式对应关系见表 1.可知s_main受初始射线走时横向导数的影响，且含波数的一次项.基于此，李代数积分和指数映射中各项的表达与时间积分中亦有所不同，通过数值算例来进行说明.图 2中红色和绿色曲线分别表示在图 1所示的速度模型中，本文算法各项系数与时间积分结果的对比，若为零即不显示.图 2a、2b和2d分别代表单平方根算子主象征一次项k_x的系数s_1，0、李代数积分象征一次项k_x系数a_1，0和单平方根算子主象征的二次项k_x²系数s_2，0，可见本算法中初始射线走时横向导数对李代数积分中各项系数的校正作用.

2.3 指数映射

若无横向变速，a#=a，根据（9）式可直接求解相位φ的表达式；在横向变速情况下，a#≠a，需要根据如下表达式进行高阶项修正^[1]：

(11)

在象征域利用Magnus算法和BCH分裂，通过级联式的递推，由低阶相位和其组成项的系数，通过交换运算得到高阶相位及其组成项的系数.二维和三维情况下相位的表达式见附录，详细推导参见文献[1]和^[22].

由于本算法中单平方根算子主象征包含波数的一次项（图 2a），使得相位系数也包含波数的一次项φ_1，0（该值在时间积分中为零）^{[1, 22]}，如图 2c所示，且相位其他项的系数亦得到了改进，以相位二次项k_x²的系数φ_2，0为例，其值如图 2e所示.

2.4 走时计算

根据稳相点原理，可得（8）式的高频渐进解，得到走时T_lie与相位φ的关系式，如下所示：

(12)

其中，Δx=（x_s-x），Δy=（y_s-y），（x_s，y_s）为地表源点的坐标，（12）式的具体推导见附录.对比T_lie²与（Δx+φ_1，0）²和（Δy+φ_0，1）²各项的对应关系，可得：

(13)

其中，c_n-i，i为走时多项式的系数，下标代表x，y方向距离的幂次；t₀为初始成像射线走时；深度z的信息通过走时系数反映出来.本算法中得到的走时二次项系数c_2，0与时间积分的对比如图 2f所示.

因此，对于一个给定的速度模型，通过坐标变换，从单平方根算子主象征s_main出发，李代数积分象征a涉及速度和走时横向导数的垂向积分和交换运算，相位多项式φ由李代数积分及其横向导数计算得到，最终求得走时的表达式.s_main、a和φ均是k_x，k_y的多项式，这里仅选取其中具有代表性的部分项系数来做图示，简要表明整个算法的计算流程.图 1是求取单平方根算子主象征之前的准备工作；图 2a，2b和2c代表各量一次项的系数，在李代数时间积分中均为零^[1]；图 2d，2e和2f对比了各量二次项的系数在本算法和时间积分中的不同，若无横向变速，即退化到时间积分的情况.

3 数值试验 3.1 点源响应

在线性横向变速介质中，地震波前面有解析解^[38]，可与之对比来验证本算法的有效性.设点源位于（x₀，z₀），速度函数为

(14)

其中，v₀为点源处的速度，v′为速度变化梯度，α为速度变化最快方向与垂向夹角，则t时刻的波前面方程为^[38]

(15)

可知线性横向变速介质中，波前面仍是圆，但圆心偏离点源的位置.将点源置于地表（750m，0m），v₀=2000m/s，dvdx=0.5s^-1，dvdz=0.5s^-1.本例中v′=，α=π/4.图 3a代表速度模型及理论走时曲线，图 3b是本算法结果与理论曲线的对比，可见两者相当吻合.

3.2 走时公式中各项对结果的影响

分析将走时公式（13）展开，可得

(16)

其中，X=Δx+φ_1，0，Y=Δy+φ_0，1.

可见，本算法是将传统Dix类型的公式推广到横向变速情况，用初始射线走时t₀加上其他各阶修正项得到最终地震波走时.针对上述模型分析各项对走时精度的影响：当t₀有10%的误差时，计算结果与理论波前的对比见图 4a；图 4b是仅包含二次校正项的计算结果；图 4c是加上了四次校正项即传统弯曲射线计算结果；图 4d是再加上三次校正项的计算结果，即非对称项的校正作用.

由图 4可知：（1）t₀的准确性对走时的影响最大，因此其准确求取对本算法的精度有决定性的作用.（2）若包含二次校正项，走时基本趋于真实值，证明了弱横向变速情况下Dix类型公式的有效性.（3）增加四次校正项，相对于图 4b，结果有所改进，但仍然偏离真实值，说明传统弯曲射线在横向变速情况下对于提高走时精度的贡献是有限的.（4）图 4d是增加了三次校正项的结果，与理论值最为吻合，表明横向变速情况下非对称项的校正作用不可忽略.

4 结语

（1）本文针对Kirchhoff叠前深度偏移中地震波走时求取计算量和存储量大的问题，基于保结构思想和拟微分算子理论，进行深入的理论研究，提出了一种新的地震波走时李代数积分算法.该算法利用坐标变换，从单平方根算子的主象征出发，在射线坐标系下进行李代数积分和指数映射的计算，运用稳相点原理，给出了地震波走时的解析表达式.

（2）李代数积分象征表达具有丰富的理论内涵和研究价值，有助于深入探索地震波场的传播规律.利用本文提出的算法在线性横向变速介质中进行了数值试验，结果表明本算法用于地震波走时求取是可行的，非对称项的加入保证了计算精度.通过对走时表达式中各项对最终结果的影响分析可知，横向变速情况下非对称项是不可忽略的.

（3）地震波走时李代数积分算法属于“现算型”方法范畴，数值运算过程是完全解耦的，便于并行加速，有利于走时的反复计算，提高了计算效率.实际计算中仅需记录走时各项的系数，由此可进行地震波走时的计算，不需要存储海量的走时数据体，节省了存储空间，为高效地进行Kirchhoff叠前深度偏移奠定了基础.

（4）同时指出，初始射线走时的准确求取和横向导数的合理计算是影响本算法数值计算精度的关键因素，需深入研究上述问题在复杂地区的具体实现，使李代数积分算法在Kirchhoff叠前深度偏移的走时求取中发挥更大的作用.

附录

算子的象征是拟微分算子作用于正弦波的响应，其表达式为

(A1)

若算子和的象征分别为σ（）和σ（），则算子乘法的象征为

(A2)

该运算称为Witt积，是运用莱布尼茨法则运算的结果.

根据保结构算法和李群理论，对于微分方程d′（ξ）=（ξ）d（ξ），其中（ξ）为拟微分算子，则其解d（ξ）满足：

(A3)

其中，d（ξ）是函数，和是算子，称为算子的李代数积分或李代数解.若用η代表李代数积分运算，即=η（），则满足如下Lie代数微分方程：

(A4)

其中，ad即代表李括号[, ]，是两个拟微分算子的交换运算结果，如下所示：

(A5)

将（A4）式展开，可表示为

(A6)

对（A6）式积分即可得到李代数积分解.

李代数积分象征表达式为二维情况下

(A7)

三维情况下

(A8)

通过指数映射求得相位的表达式为二维情况下

(A9)

三维情况下

(A10)

从二维出发，大步长波场延拓积分解的表达式为

(A11)

其中φ含有对深度z的积分信息.令p_x=k_x/ω，则

(A12)

(A12)式是关于慢度的积分，与频率无关，根据稳相点原理可知其高频渐进解为

(A13)

其中φ″_px为慢度的二阶导数.可得Kirchhoff偏移公式为

(A14)

其中，x_s是输入道地表的坐标.令Δx=x_s-x，b（p_x）=，则b（p_x0）代表地震波的走时，p_x0是∂b/∂p_x的零点.

由（A9）式可知

(A15)

由此可知

(A16)

根据相位φ可对（A16）式进行求解.同理，三维情况下即可得（12）式.