2010, Vol. 53
2. 中国科学院研究生院, 北京 100049
2. Graduate University of Chinese Academy of Sciences, Beijing 100049, China
地震波传播过程有三种等价的物理描述体系：牛顿体系、拉格朗日体系和哈密顿体系.它们分别可用微分方程、变分原理和群论进行数学表达，并通过离散化得到有限差分法、有限元方法和群的保结构算法.研究大尺度简单地质问题时，前两类方法是有效的，但对于小尺度复杂的地质问题，则需要探索新的算法理论[1].目前研究表明，复杂地质条件下地震波传播规律涉及到拟微分算子理论[2, 3]和群的基本性质与规则[4]，其表达基于保结构思想[5]：即地震波传播过程中波前的数学表达是一个连续流形，其离散化后的过程描述亦应该在此连续流形上.保结构算法的目的是使离散化数值计算结果中尽可能多地保留原物理系统（如运动学系统或动力学系统）的性质，例如群的性质和频散关系.当前发展起来的保结构李群算法[6, 7]，是辛几何保结构思想的拓展[5, 8, 9]，应用于地震波研究可使地震波数值计算具有保真性、稳定性和快捷性.
长期以来，地球物理界努力引进拟微分算子[2, 3]和象征理论[3, 10]对地震波传播进行抽象研究，挖掘其深层次的规律，提高数值计算精度和速度[11~15].本文利用保结构李群算法[6, 7]研究单程波，即对算子进行大步长积分延拓[16, 17]，目的是在横向变速情况下，简化地震波走时求取，便于走时的反复计算，称为李代数积分走时算法.该方法将地震波场延拓方程转化为李代数微分方程，求其李代数积分解，通过指数映射得到波场延拓算子的复相位，在离散计算中减少递推过程的误差积累，有效保证了计算精度[18~20].因此，在大步长延拓研究中引入保结构算法，需要解决李代数积分[6, 21]和指数映射计算两大问题，前者通过将微分方程转换到象征域，使计算得到简化，后者通过BCH（Bacher-Campbell-Hausdorff）公式[22, 23]，利用Magnus算法[24, 25]进行求解，由稳相点原理[26, 27]得到地震波走时的表达式.
目前李代数时间积分走时算法已在Kirchhoff叠前时间偏移中取得了很好的应用效果，提高了计算精度和效率[28~30].与李代数时间积分相比，深度域中单平方根算子透镜项与速度有关，在横向变速情况下存在空间变化，且含有波数的一次项[16]，直接进行李代数深度积分需要高阶修正，收敛很慢.本文对该问题进行详细讨论，并利用坐标变换，对李代数深度积分算法进行改进，在射线坐标系下将透镜项的系数转换为常数，避免与后续项的交换运算，加快了收敛速度，在保证精度的同时减少了计算时间.
2 单程波李代数深度积分频率空间域单程波波场延拓方程可表示为
|
(1) |
其中，
根据保结构李群算法[6, 7]，方程（1）的解W（x，y，z，ω）满足
|
(2) |
其中，A称为单平方根算子S的李代数积分或李代数解，可表示为
|
(2) |
其中，η为李代数积分运算，[A，S]=AS-SA为算子的交换运算[4]，Bn为Bernoulli数[17].
令G=exp（iA）=exp（iη（S）），即W（x，y，z，ω）=GW（x，y，0，ω），其中G为大步长单程波算子，可通过上述李代数积分进行计算，对地震波的运动学和动力学多种参数进行研究[28~32]，本文主要利用李代数积分算法求取地震波的走时.
当介质存在横向变速时，上述算子均为拟微分算子，为了简化运算，将计算转换到象征域进行[10].象征是频谱概念的推广，即拟微分算子
|
(4) |
将（4）式转换到象征域，令σ（S）=s，σ（A）=a，a=a1，1 +a2，1 +a3，1 +a4，1 +a5，1 + …，其中，上标代表a的各阶校正项，由以下递推方式计算[16]：
|
(5) |
其中，a1，2=s，a2，2=i{a1，s}，a3，2=i{a2，1，a1，2}，a3，3=i{a1，1，a2，2}，a4，2=i{a3，1，a1，2}，a4，3=i{a1，1，a3，2}，a4，4=i{a2，1，a2，2}；s为单平方根算子S的象征，是关于ω，kx，ky的多项式，若只考虑其主象征[16]，可表示为
|
(6) |
其中，si，j为多项式的系数，下标i，j分别代表kx，ky的幂次；{}表示交换算子的象征，定义为{a，b}=σ（AB-BA），在横向变速情况下算子相乘的象征不等于象征简单的相乘，满足如下关系[10]：
|
(7) |
其中，a，b代表算子A，B的象征，#为Witt积[10].因此交换算子的象征可表示为{a，b}=σ（AB-BA）=a#b-b#a.
根据（5）式，利用单平方根算子的主象征，可求得象征域李代数深度积分的表达式：
|
(8) |
其中，ai，j为多项式的系数，下标i，j分别代表kx，ky的幂次.
3 存在的问题/精度分析将波场延拓方程的李代数深度积分解（2）式转换到象征域求解[16, 17]，得到
|
(9) |
其中，Ŵ=σ（W）表示波场的象征，深度z的信息包含在积分运算σ（G）中.令φ表示象征运算σ（G）的复相位，即σ（G）=exp[iφ（x，y，kx，ky，ω）]，可知
|
(10) |
根据象征理论，exp（iφ）=exp[σ（G）]=exp{σ[exp（iA）]}=exp（ia#）[10, 16]，若无横向变速，可直接求得复相位φ=a，横向变速情况下，需要根据Witt积加入修正项，表示为
|
(11) |
其中，φi，j为多项式的系数，下标i，j分别代表kx，ky的幂次.
对方程（10）利用稳相点原理[26, 27]，可得如下方程组：
|
(12) |
其中，px0，py0分别代表稳相点处x，y方向的慢度.通过对比走时Tlie2与（x+φ0，1）2和（y+φ1，0）2各阶项的对应关系，可知
|
(13) |
其中，cn-i，i代表走时多项式的系数，下标表示x，y方向距离的幂次，深度z的信息通过走时系数反映出来.在数值计算中若保留到走时的四次项，共需要12个系数，均与相位系数φ0，0，φ2，0有关[28].
综上所述，走时系数由相位系数得到，最关键的两个参数是φ0，0，φ2，0，相位通过φ=a#进行计算，因此最终走时计算受到单平方根算子李代数深度积分的影响.而在李代数深度积分中，由于横向变速情况下单平方根算子的透镜项存在空间变化，且含有波数的一次项，李代数积分a需要进行高阶修正.但实际数值计算中不可能无限次增加高阶项，同时由于高阶项涉及到速度的高阶导数，在复杂地质条件下会出现畸变，且计算时间随着校正阶数的增加而增加，如图 1所示.
|
图 1 李代数深度积分计算时间随校正阶数的变化 Fig. 1 The calculation time of Lie algebra depth integral with correction items |
为了提高运算效率和稳定性，在具体实施中，需保留有限的校正项参与计算，因此要求该积分能快速收敛.但由于透镜项的持续作用使得李代数深度积分的各幂次均需要高阶修正，收敛很慢，若仅保留低阶项，会由于校正不足影响后续相位和走时的计算精度.李代数积分各项的校正情况如表 1所示，这里仅显示前五阶校正项.
|
表 1 李代数深度积分方法中各项的校正情况 Table 1 The coefficients correction of items in Lie algebra depth integral |
由于线性横向变速介质中走时有解析解，与之对比来说明李代数深度积分存在的问题.设理论模型速度函数为
|
(14) |
其中，v0为源点处的速度，v′为速度变化梯度，α为速度变化最大方向与垂向夹角，点源位于（x0，z0），则t时刻波前面方程为[33]
|
(15) |
将点源放置在地表（750m，0m）处，v0=2000m/s，dvdx=0.5s-1，dvdz=0.5s-1，即
|
图 2 横向线性变速介质中李代数深度积分方法计算走时与理论值的对比和误差曲线 (a)计算走时(绿线)与理论值(红线)的对比(0.2 s); (b)两者的误差曲线(0.2 s). Fig. 2 Comparison between travel time by Lie algebra depth integral method and theoretical value and error curve in medium with linear lateral velocity variation (a) Comparison between computed travel time (0.2 s) (green) and theoretical value (red); (b) Error curve (0.2 s). |
针对李代数深度积分存在的问题，本文通过坐标变换，在射线坐标系下进行计算，将单平方根算子的透镜项转换为常数，消除其与后续项的作用，减少高阶修正，使改进后的算法收敛更快，进而利用指数映射和稳相点原理得到旁轴走时的表达式.
对频率空间域波动方程进行如下坐标变换：
|
(16) |
其中，T=t0为成像射线的走时.将其转换到射线坐标系下，可得
|
(17) |
其中，T′x，T′y，T′z分别代表走时T沿x，y，z方向的导数，x′，y′是射线坐标，用来惟一地确定一条射线，后文统一用x，y表示.将（17）式转换到象征域，表示为
|
(18) |
其中，s代表射线坐标系下单平方根算子的象征，若忽略高阶项的校正作用，可知主象征
|
(19) |
将（19）式展开可得
|
(20) |
从（20）式可知，坐标变换之后透镜项系数转化为常数，与后续项交换运算为零[16, 17]，从而在李代数积分计算中避免高阶修正.在弱横向变速介质中，Tx和Ty变化很小，忽略其对si，j的校正作用，则改进后李代数积分中kx，ky各幂次的校正情况如表 2所示.可以看出，系数a0，0，a2，0，a0，2不需要二阶以上的修正项，且低于四次幂的系数只需要三阶校正.实际计算时在si，j中加上Tx，Ty的修正对上述性质影响不大，从而保证了李代数积分能快速收敛.
|
|
表 2 改进的李代数深度积分方法各项系数的校正情况 Table 2 The coefficients correction of items in improved Lie algebra depth integral |
综上所述，改进后的算法通过坐标变换加快了李代数积分的收敛速度，通过指数映射和稳相点原理进行走时计算，不改变算法的保结构性质.李群算法指出，利用指数表示，可以将演化的函数限定在李群流形上，而李群可以使得映射前后的向量是保结构的.本文计算走时的基础是研究相位，它实际上是李代数指数映射（李群）象征，即纯虚数的指数，在相空间模为1，因而该映射保持了相空间（复空间）的向量模长，从而是保结构的（这里的结构，在复空间为酉结构，对应于实空间的辛结构）[1].
因此对于一个给定的速度模型，单程波李代数深度积分走时计算的详细步骤如下：首先利用式（16）的坐标变换将深度域单平方根算子的透镜项转化为常数；然后在射线坐标系下根据式（20）计算单平方根算子的主象征；再由递推的方式根据式（5）计算李代数积分的象征；通过指数映射由式（11）计算单程波延拓方程的复相位；最后根据稳相点原理利用式（13）得到地震波旁轴走时的解析表达式.
对于图 2的线性横向变速模型，利用改进后的算法计算走时，李代数积分保留三阶校正项，计算结果如图 3a中绿线所示，红色代表理论曲线.图 3b比较了0.2s处各阶校正项对走时结果的影响，绿线、红线和蓝线分别表示保留三阶校正、二阶校正和一阶校正计算结果与理论值的误差曲线，可见即使只保留一阶校正项，相对于改进前的算法，误差也有明显减小，表明改进后的方法使数值计算精度得到了很大的提高.
|
图 3 横向线性变速介质中改进后算法计算走时与理论值的对比和误差曲线 (a)计算走时(绿线)与理论值(红线)的对比；(b)两者的误差曲线(0.2 s). Fig. 3 Comparison between travel time by improved method and theoretical value and error curve in medium with linear lateral velocity variation (a) Comparison between computed travel time (green) and theoretical value (red); (b) Error curve (0.2 s). |
对于图 4背景所示的速度模型，含有不连续面和断层，地震波走时无解析解.将改进后的算法（LAI）与有限差分（FD）和射线追踪（RT）进行对比，如图 4中曲线所示，可见三者基本吻合，在边界和深部LAI与FD和RT的误差增大，主要由于在横向变速大的介质中，表 2所示的收敛性质会受到一定程度的影响，且旁轴走时算法在大角度时存在误差，有待进一步深入研究，将李代数积分法走时计算的实用化从Kirchhoff叠前时间偏移[34, 35]推广到叠前深度偏移.
|
图 4 LAI算法(黄色)与有限差分(红色)和射线追踪(绿色)方法计算的走时等值线对比 Fig. 4 Comparison of travel time contours computed by LAI (yellow), FD (red) and RT (green) |
若利用改进后的算法进行Kirchhoff叠前深度偏移走时计算，只需记录目标区域内走时多项式系数，提高了计算效率，且不需要存储海量的走时数据体[28~30]，节省了存储空间.对于图 4所示的速度模型，上述三种方法（LAI、FD和RT）的计算效率曲线如图 5所示，炮点均置于地表，样点数目为750× 750，RT算法中每间隔2°发射一条射线，每炮共追踪90条射线.可见LAI算法的计算时间与有限差分基本相当，远远小于射线追踪耗费的时间.对于存储量而言，理论情况下，若事先计算N炮走时用于偏移，不考虑孔径等因素，二维情况下，LAI算法计算到走时的四次项仅需要存储4个系数，对一个模型若计算nx·nz个样点的走时，需要开辟4nx·nz个数据存储空间，而有限差分和射线追踪方法均需存储N·nx·nz数据；三维情况下，LAI算法计算到走时的四次项仅需要存储12个系数，对一个模型若计算nx·ny·nz个样点的走时，需要开辟12nx· ny·nz个数据存储空间，而有限差分和射线追踪方法均需存储N·nx·ny·nz个数据.以上均为理论情况下的对比，在实际计算中，由于不同的精度要求以及偏移孔径的选择等因素，可能使得存储量的大小有别于上述结果，但LAI方法节省存储空间的优势依然存在.
|
图 5 LAI算法与有限差分和射线追踪方法计算效率对比 Fig. 5 Computing efficiency comparison of LAI, FD and RT with shot numbers |
保结构李群算法和象征运算具有丰富的理论内涵，可用于对地震波传播深层次规律的研究.本文主要推导了单平方根算子李代数深度积分的具体表达，详细分析了横向变速情况下，李代数深度积分中由于单平方根算子透镜项存在空间变化引起的问题.我们利用坐标变换将透镜项的系数转换为常数，在射线坐标系下进行计算，避免其与后续项的交换运算，使改进后的李代数积分算法能快速收敛.利用指数映射和稳相点原理，得到了旁轴走时的解析表达式，实际计算中只需要存储走时的系数，易于并行实现.通过对比线性横向变速模型中的计算误差，表明改进后的方法使计算精度得到了很大的提高.数值算例表明，与有限差分和射线追踪方法相比，在介质存在弱横向变速的情况下，改进后的算法在保证计算精度的同时，极大地提高了计算效率并节省了存储空间，对于提高Kirchhoff叠前深度偏移的精度和效率非常有利.但需要指出的是，单平方根算子的象征需要计算初始走时和速度的横向导数，在地质情况复杂的地区导致算法不稳定，且计算的旁轴走时在大角度时存在误差，因此需对李代数积分算法开展进一步的研究和探索.
致谢感谢匿名审稿专家细致且富有建设性的意见，本文编辑提出了详尽的修改建议，在此一并致以诚挚的谢意.
| [1] | 刘洪, 李幼铭. 油气勘探二次创业的油储地球物理方法研究回顾. 石油与天然气地质 , 2008, 29(5): 648–653. Liu H, Li Y M. Review on geophysical methodology of oil reservoirs in the second round of hydrocarbon exploration. Oil & Gas Geology (in Chinese) , 2008, 29(5): 648-653. |
| [2] | Cordes H O. The Technique of Pseudo Differential Operators. Cambridge: Cambridge University Press, 1995 . |
| [3] | 陈恕行. 拟微分算子理论(第二版). 北京: 高等教育出版社, 1995 . Chen S X. Pseudo Differential Operator (Edition 2) (in Chinese). Beijing: Higher Education Press, 1995 . |
| [4] | Naimark M A, Stern A I. Theory of Group Representations. New York: Springer-Verlag New York Inc., 1982 . |
| [5] | Feng K. On difference schemes and symplectic geometry. In:Feng K ed. Symplectic Differential Geometry and Differential Equations. Beijing:Science Press, 1985 http://www.oalib.com/references/17659643 |
| [6] | Iserles A, Munthe-Kaas H, Nfrstt S, et al. Lie-group methods. Acta Numerica , 2000, 9(1): 215-263. |
| [7] | 田畴. 李群及其在微分方程中的应用. 北京: 科学出版社, 2003 . Tian C. Applications of Lie Groups to Differential Equations (in Chinese). Beijing: Science Press, 2003 . |
| [8] | 陈景波, 秦孟兆. 射线追踪、辛几何算法与波场的数值模拟. 计算物理 , 2001, 18(6): 481–486. Chen J B, Qin M Z. Ray tracing, symplectic geometry algorithm and numerical modeling of wave field. Computational Physics (in Chinese) , 2001, 18(6): 481-486. |
| [9] | 秦孟兆, 陈景波. Maslov渐近理论与辛几何算法. 地球物理学报 , 2000, 43(4): 522–533. Qin M Z, Chen J B. Maslov asymptotic theory and symplectic algorithm. Chinese J. Geophys. (in Chinese) , 2000, 43(4): 522-533. DOI:10.1002/(ISSN)2326-0440 |
| [10] | 齐民友, 徐超江, 王维克. 现代偏微分方程引论. 武汉: 武汉大学出版社, 2005 . Qi M Y, Xu C J, Wang W K. Introduction of Modern Partial Differential Equation (in Chinese). Wuhan: Wuhan University Press, 2005 . |
| [11] | Zhang G Q. A new algorithm for finite-difference migration of steep dips. Geophysics , 1988, 53(2): 167-175. DOI:10.1190/1.1442451 |
| [12] | 张关泉. 一维波动方程的反演问题. 中国科学 , 1988, 31(7): 707–721. Zhang G Q. One-dimensional wave equation inverse problem. Science China (in Chinese) , 1988, 31(7): 707-721. |
| [13] | 张关泉. 波场分裂、平方根算子与偏移.反射地震学论文集. 上海: 同济大学出版社, 2000 . Zhang G Q. Wave field split, square root operator and migration. Memoir of Reflection Seismology (in Chinese). Shanghai: Tongji University Press, 2000 . |
| [14] | Chen J B, Munthe-Kaas H, Qin M Z. Square-conservative schemes for a class of evolution equations using Lie-group methods. SIAM J. Numer. Anal. , 2002, 39: 2164-2178. DOI:10.1137/S0036142901383971 |
| [15] | 刘超颖, 王成祥. 非稳态相移法叠前深度偏移. 地球物理学进展 , 2004, 19(3): 543–546. Liu C Y, Wang C X. Prestack depth migration by nonstationary phase shift method. Progress in Geophysics (in Chinese) , 2004, 19(3): 543-546. |
| [16] | 刘洪, 袁江华, 陈景波, 等. 大步长波场深度延拓的理论. 地球物理学报 , 2006, 49(6): 1779–1793. Liu H, Yuan J H, Chen J B, et al. Theory of large-step wave field depth extrapolation. Chinese J. Geophys. (in Chinese) , 2006, 49(6): 1779-1793. |
| [17] | 刘洪, 王秀闽, 曾锐, 等. 单程波算子积分解的象征表示. 地球物理学进展 , 2007, 22(2): 463–471. Liu H, Wang X M, Zeng R, et al. Symbol description to integral solution of one-way wave operator. Progress in Geophysics (in Chinese) , 2007, 22(2): 463-471. |
| [18] | 王秀闽. 积分偏移的计算几何与并行实现. 北京: 中国科学院地质与地球物理研究所, 2008 . Wang X M. Computational geometry and parallel implementation of integral migration (in Chinese). Beijing: Institute of Geology and Geophysics, Chinese Academy of Sciences, 2008 . |
| [19] | 刘国峰, 刘洪, 王秀闽, 等. Kirchhoff积分时间偏移两种走时计算及并行算法. 地球物理学进展 , 2009, 24(1): 131–136. Liu G F, Liu H, Wang X M, et al. Two kinds of traveling time computation and parallel computing methods of Kirchhoff migration. Progress in Geophysics (in Chinese) , 2009, 24(1): 131-136. |
| [20] | Liu G F, Liu H, Li B, et al. Getting pre-stack time migration travel times from the single square root operator. Applied Geophysics , 2009, 6(2): 129-137. DOI:10.1007/s11770-009-0016-z |
| [21] | Celledoni E, Iserles A. Methods for the approximation of the matrix exponential in a Lie-algebra setting. IMA J. Numer. Anal. , 2001, 21(2): 463-488. DOI:10.1093/imanum/21.2.463 |
| [22] | Moler C, Van Loan C F. Nineteen dubious ways to compute the exponential of a matrix. SIAM Review , 1978, 20(4): 801-836. DOI:10.1137/1020098 |
| [23] | Liu Y, Wu R S. A comparison between phase screen, finite difference and eigenfunction expansion calculation for scalar waves in inhomogeneous media. Bulletin of the Seismological Society of America , 1994, 84(4): 1154-1168. |
| [24] | McLachlan R I, Quispel G R W. Splitting methods. Acta Numerica , 2002, 11(1): 341-434. |
| [25] | Magnus W. On the exponential solution of differential equations for a linear operator. Communications on Pure and Applied Mathematics , 1954, 7(5): 649-673. |
| [26] | 李世雄. 波动方程的高频近似与辛几何. 北京: 科学出版社, 2001 . Li S X. High-Frequency Approximation and Symplectic Geometry (in Chinese). Beijing: Science Press, 2001 . |
| [27] | Chen J. Specular ray parameter extraction and stationary phase migration. Geophysics , 2004, 69(1): 249-256. DOI:10.1190/1.1649392 |
| [28] | 刘洪, 刘国峰, 李博, 等. 基于横向导数的走时计算方法及其在叠前时间偏移中的应用. 石油物探 , 2009, 48(1): 3–10. Liu H, Liu G F, Li B, et al. Travel time calculation by lateral derivative and its application in pre-stack time migration. Geophysical Prospecting for Petroleum (in Chinese) , 2009, 48(1): 3-10. |
| [29] | 李博, 刘国峰, 刘洪. 地震叠前时间偏移的一种图形处理器提速实现方法. 地球物理学报 , 2009, 52(1): 245–252. Li B, Liu G F, Liu H. A method of using GPU to accelerate seismic pre-stack time migration. Chinese J. Geophys. (in Chinese) , 2009, 52(1): 245-252. |
| [30] | 刘国峰, 刘钦, 李博, 等. 油气勘探地震资料处理GPU/CPU协同并行计算. 地球物理学进展 , 2009, 24(5): 1671–1678. Liu G F, Liu Q, Li B, et al. GPU/CPU co-processing parallel computation for seismic data processing in oil and gas exploration. Progress in Geophysics (in Chinese) , 2009, 24(5): 1671-1678. |
| [31] | 刘洪, 袁江华, 勾永锋, 等. 地震逆散射波场和算子的谱分解. 地球物理学报 , 2007, 50(1): 240–247. Liu H, Yuan J H, Gou Y F, et al. Spectral factorization of wave field and operator in seismic inverse scattering. Chinese J. Geophys. (in Chinese) , 2007, 50(1): 240-247. |
| [32] | 刘洪, 何利, 刘国锋等.地层滤波公式的李代数积分证明和推广.金翔龙, 秦蕴珊, 朱日祥等.中国地质与地球物理研究进展--庆贺刘光鼎院士80华诞.北京:科学出版社, 2008. 412~422 Liu H, He L, Liu G F, et al. Derivation and generalization of the stratigraphic filtering formula via Lie algebra integral (in Chinese). Jin X L, Qin Y S, Zhu R X, et al. The Research Development of Geology and Geophysics in China-The Collection in Memory of 80th Anniversary of Academician Liu Guangding. Beijing:Sciences Press, 2008. 412~422 |
| [33] | Nolet G.地震层析成像.王椿镛, 吴宁远, 李幼铭等译.北京:学术书刊出版社, 1989 Nolet G.Seismic Tomography (in Chinese).Translated by Wang C Y, Wu N Y, Li Y M, et al.Beijing:Academic Book Press, 1989 |
| [34] | 周纳, 刘继勇. 非对称走时叠前时间偏移方法及应用. 油气地球物理 , 2010, 8(2): 48–52. Zhou N, Liu J Y. The method and application of asymmetric travel-time pre-stack time migration. Petroleum Geophysics (in Chinese) , 2010, 8(2): 48-52. |
| [35] | 陈新荣. 非对称走时叠前时间偏移在青东4东地区的应用. 勘探地球物理进展 , 2010, 33(4): 279–283. Chen X R. The application of asymmetric travel-time pre-stack time migration in area of Qingdong 4 dong. Progress in Exploration Geophysics (in Chinese) , 2010, 33(4): 279-283. |