基于改进的双线性旅行时插值的三维射线追踪

梅胜全,邓 飞,钟本善,周熙襄

(成都理工大学 信息工程学院,四川成都 610059)

0 前言

三维地震射线走时计算是地震层析成像的基础,随着三维高精度、高密度地震勘探的广泛应用,其作用也越来越突出。而走时计算的精度与效率问题[1、2],也一直是三维地震层析成像难以有效解决的难点。地震射线走时计算的精度与效率是一个矛盾问题,一方面较高的射线精度能更准确地反映地震波场主能量传播轨迹,但计算复杂度高,计算量大,计算效率会降低;另一方面高效、快速的走时计算,却往往是以牺牲射线走时精度为代价的,这成为一个难以两全的技术难点[1]。

三维地震走时计算的精度是关键,射线精度越高,其反映地震波场主能量传播的射线路径越合理,地震波走时拟合累积误差越小,层析反演收敛就越快,反演的速度结构也越接近真实情况,这在我国西部复杂近地表层析成像中尤其突出。若射线精度达不到一定的要求,其层析反演解的不适定性就严重,地震资料处理的成像效果就差。同时,射线走时计算效率也是三维高精度三维地震勘探中的一个瓶颈问题,由于炮点、检波点加密,道间距变小,导致射线走时计算数据规模呈指数级增长,计算效率成为三维地震层析成像实际应用的一个巨大障碍。那么,如何有效地解决三维地震走时计算的精度与效率问题呢?

由于旅行时射线追踪算法已经很成熟了[5],大多以初值问题的试射法与边值问题的弯曲法为基础,目前应用广泛的有最短路径算法[6～8]和旅行时插值算法(LTI)[9]。众多的文献与应用表明,双线性旅行时插值算法能很好地模拟初至波传播过程,计算的初至时间具有较高的精度[10、11],但其计算量大,计算效率相对不高。

作者在三维高精度地震层析成像应用中,结合最短路径算法和旅行时插值算法,通过对双线性插值算法的改进,提出了一种三维射线追踪算法。采用最短路径算法寻找地震初至波的波前面,利用改进的双线性旅行时插值算法,计算网格单元的插值点,进行正演计算射线走时。这样既保证了射线的精度,同时又简化计算,提高了计算效率。

1 三维初至波射线追踪原理

在进行三维射线追踪时,首先应对模型进行离散化处理,将地层模型转化为以规则网格(一般为长方体)为基本单元的三维网格模型,假设单元内为常速,射线追踪算法分为二个步骤:

(1)向前处理计算单元结点的时间场分布,从激发点出发,基于波前面扩展,利用初至波旅行时插值算法计算整个模型节点上的最小走时。

图1 向前处理过程Fig.1 Forward processing procedure

(2)向后处理则利用互易原则,反向追踪射线路径,从接收点开始利用旅行时插值公式反向计算各单元的插值,直至激发点,得到一条完整的初至射线路径。

1.1 向前处理:计算模型节点时间场

向前处理过程用于计算激发点所覆盖的网格节点的最小旅行时间,这是最核心也最复杂的一步,基本的算法步骤如下:

(1)在激发点单元,直接计算网格单元节点的旅行时间,如图1(a)所示。

(2)以激发点单元为中心的,利用旅行时插值算法,逐一计算出六个相邻单元的网格节点的旅行时间,如图1(b)所示。

(3)使用最短路径射线追踪算法,求取新的次震源网格单元(假定次震源点为网格单元中心点),寻找最小走时的新的次震源点,继续向外插值计算或更新相邻的非次震源单元网格的节点旅行时,如图1(c)所示。

(4)重复执行步骤(3),直到整个模型的网格节点都作为次震源被计算完为止。

1.2 向后处理:反向追踪射线路径

向前处理计算了全部网格节点最小旅行时间场的分布。采用双线性插值算法计算射线的入射点,可以从长方体单元中任意穿过,而不是在固定的网格节点,或次震源路径上,因此不能简单通过保留次震源点,直接连接作为射线路径,还需要利用旅行时插值算法,根据互换原则,反向确定满足费马原理的最小旅行时节点位置,逐一插值计算出射线在网格单元的交点,直至回到激发点,连接而成为一条完整的射线路径。

具体步骤如下:

(1)利用旅行时插值算法,对检波点所在单元的六个面分别计算出初至射线穿出网格单元时与网格单元的交点,找出最小旅行时的插值点即为射线入射点,如图2(a),这个最小旅行时插值就是检波点的初至旅行时。

(2)根据当前交点的位置,利用旅行时插值算法,继续找出下一个网格单元的交点如图2(b)所示,需要对交点的位置进行判定。如果位于单元的一个面上,则由下一个网格单元决定;如果位于网格单元的棱线或者顶点上,则需要用周围多个网格单元共同确定下一个插值最小的交点。

(3)如果射线已经回到激发点所在网格单元,则已完成一条射线路径的追踪,否则继续执行步骤(2)。

图2 向后处理过程Fig.2 Backward processing procedure

2 基于最短路径的三维次震源波前面扫描算法

在向前处理过程中,作者在本文采用最短射线路径算法不断计算次震源点,通过次震源点的最小旅行时排序,实现地震波前面的扩展。并以次震源为中心,推算临近网格节点的最小旅行时,直至计算出整个模型中所有网格节点的最小旅行时。作者采用堆排序算法进行波前面扫描,将网格单元分成二个集合,已经获得了最小旅行时的次震源集合与由其它网格插值计算旅行时,但还未成为次震源的集合。每一次从堆中取出最小旅行时的网格单元,作为次震源点,计算相邻的六个网格,并加入堆中。如此反复,直到所有的网格单元都成次震源点,计算完所有单元最小旅行时为止。

在向后处理插值反向追射线路径时,不需要简单地插值计算六个面的最小旅行时交点,而是先计算当前交点与其相邻次震源的最小旅行时,完成次震源传播路径的反推,从而指导射线的入射方向。再应用改进的双线性旅行时插值算法,快速求出下一个交点,直至整条射线路径追踪至激发点。

这里采用Moser计算公式计算插值交点与邻近网格次震源点的旅行时差[12],

其中 dij为插值交点与次震源点的欧氏距离;¯s为相邻网格计算单元的平均速度求得的平均慢度。

找出相邻网格单元的次震源点与当前插值交点的最小旅行时,则相应网格单元为射线下一个入射的方向。

通过这步次震源点扫描,可以最大程度地减少插值计算次数。

3 双线性旅行时插值原理及改进

旅行时插值算法的目的,是利用已知的最小旅行时间节点,计算未知节点的最小旅行时间。在整个向前及向后处理过程中,都需要多次使用旅行时插值计算次震源附近的网格节点的最小旅行时间,对旅行时插值算法的优化,可以在很大程度上提高整个算法的效率。作者在Asakawa提出的二维线性旅行时初至算法(LTI)基础上[13],进行算法推导,得出简化的插值公式与判定条件,下面是详细的算法说明。

在图3中,A、B、C、D为四个已知最小旅行时的节点,构成一个空间矩形,在X方向的宽度为w,在Y方向的高度为h,记A点坐标(0,0),最小旅行时记为t0,B点坐标(1,0),最小旅行时记为t1,C点坐标(1,1),最小旅行时记为t2,D点坐标(0,1),最小旅行时记为t3,F点是待求的未知最小旅行时节点,相对于A点的坐标记为(x0,y0,z0),网格单元的慢度为s。

图3 三维双线性旅行时插值示意图Fig.3 Three-dimensional bilinear travel-timeinterpolation scheme

在矩形ABCD上旅行时按双线性变化,满足

通过A、B、C、D四点,可以求得式(3)的系数:a=t0,b=t1-t0,c=t3-t0,d=t0+t2-t1-t3

在矩形ABCD内求一点E,坐标为(x,y),旅行时为Te,通过E点到达空间F点的旅行时为T。其中

要得到最小旅行时间,对式(4)求偏导得到方程组:

从式(5)中消去L后,可以得出:

对式(6)展开后,可得到:

从式(7)可见,最小旅行时点位于一个原点(-(cw-dx0)/2d,-(bh-dy0)/2d)的双曲线上。图4(见下页)是一个可能解的示意图,假设与矩形ABCD边界相交有二个点H、K,则最小旅行时点应该位于之内。这是一簇解,一般采用迭代法求解,如牛顿～拉斐森方法,可以在有限步迭代收敛于真实解附近。

图4 双线性旅行时插值可行解示意图Fig.4 The schematic diagram of bilinear travel-time interpolation feasible solution

但在三维射线追踪,尤其在三维高密度勘探中,炮点数据与检波点数据量大,这样插值计算将会成为射线追踪,耗费大量计算时间,严重影响插值计算效率。基于插值效率考虑,这里进行插值算法改进,通过必要的近似处理,可以有效地减化运算。

经变换可解得:

同理可得:

这就得到了一个快速求解双线性插值算法的公式,避免了迭代求解,可以直接快速得出插值点,提高了计算效率。当然与迭代法相比,在精度上有一定降低,但若控制在一个合理范围内,则射线精度是能得到保证的。

下面给出快速插值计算的步骤:

(1)由式(7)计算与矩形ABCD的交点,如果没有交点,则在边界上应用二维LTI得出二个解析解,比较求出最小的旅行时点作为插值点。

(2)应用式(10)、式(11)快速求出插值点,如果不满足判定条件,则直接比较计算矩形ABCD的交点(如H、K点)的最小旅行时点作为插值点。

改进的双线性插值算法的优点:式(10)、式(11)计算简单,通过波前面扫描算法与判定条件,能快速求出一个接近真实的解,最好情况是一个网格只需计算一次,最坏情况是计算边上插值5次。相比于牛顿～拉斐森方法迭代法而言,一个网格的每个面平均8次～12次的插值计算,六个面的计算最好也有50次,最差70次以上。由此可以看出,计算量大幅度降低,同时射线精度也能最大限度得到保证。

图5 一个复杂的三维山地资料Fig.5 A complex 3D mountain data example

图6 改进算法的射线追踪(局部)Fig.6 The improved algorithm of ray tracing(local)

4 资料试算及误差、效率分析

通过对一个复杂山地资料的射线追踪实际资料处理,采用纵向速度梯度变化建立模型,地表速度为500 m/s,每10 m速度增量为5 m/s,炮点7 000炮,检波点20 000个,炮水平覆盖范围11 000 m,模型最大深度1 500 m,道间距25 m,最大偏移距2 500 m,采用普通微机进行计算,平均单炮正演时间4 s左右。

为了验证改进的双线性插值算法的有效性和可靠性,通过牛顿～拉斐森方法迭代插值计算与改进的双线性插值算法快速计算的结果进行对比,这里以二种算法计算的走时对比基本上能反映其走时计算的精度(如图7及下页图8所示)。图9(见下页)是二种计算的运算效率(每炮正演计算用时)对比。

在图7中,通过对二种算法的单炮射线走时来反映二种计算方法的差异。图8中最大差值在20 ms,主要集中在近炮点附近,这说明在近炮点附近的计算误差较大。当偏移距越大,二者走时差值越小,平均在5 ms左右。同时可以看出,基于迭代计算的走时曲线比较光滑,而快速插值计算的走时曲线起伏较大,从炮点左右射线情况来看,炮点右边射线走时差值更小。经过分析发现,与炮点左右二边的地形起伏较大,网格间速度变化剧烈有较大的相关性。总的来说,二种算法计算的走时差值是在一个合理的范围内,通过快速插值计算的走时可以代替迭代计算走时,计算精度能得到保证。

图7 单炮的各条射线走时对比(ms)Fig.7 A shot's all rays travel-time contrast(ms)

图8 二种算法单炮射线走时差值对比(ms)Fig.8 A shot's all rays travel-time differential contrast with two algorithms(ms)

图9是选取800炮计算用时的一个对比情况。总的来说,每炮计算用时差异不大,但二种算法的单炮计算用时则非常明显,显然,改进的插值算法,计算效率得到了明显提升。

图9 800炮的射线追踪运算效率对比(s)Fig.9 The 800 shots'ray tracing operation efficiency contrast(in seconds)

5 结论

作者在本文提出的三维射线追踪算法,在一定程度上解决了射线追踪的精度与效率问题。在最短路径算法的基础上,提出了以次震源波前面扫描算法,代替了旅行时插值比较,从而减少了插值计算量。同时旅行时双线性插值改进算法原理简单,通过判定条件进行快速插值,既保留了旅行时插值算法较高精度的优点,同时简化计算,效率也得到了明显提升。通过数值计算表明,本方法可以适应各种复杂的网格模型,能够较准确地计算地震波的初至时间和射线路径,射线走时精度与效率能够满足层析反演需要。

在实际的复杂近地表地区资料处理中,采用此方法进行射线追踪正演模拟和反演成像,能获得较精确的近地表模型。从理论模型与实际资料试算数值计算对比分析,说明了方法的计算精度和速度是可靠和有效地,也说明了在现有计算机硬件条件下进行大规模的精确三维初至射线正演是可行的。

[1] JUL IAN B R,GUBBL INSD.Three-dimensional seismic ray tracing[J].J.Geophys,1977,43:95.

[2] N ISH I K.A three-d imensional robust seismic ray tracer for volcanic regions[J].Earth Planets Space,2001,53:101.

[3] 常旭,刘伊克.地震正反演与成像[M].北京:华文出版社,2001.

[4] 邓飞,周杲,王美平,等.基于剖面重构的三维地层建模[J].物探化探计算技术,2007,29(1):19.

[5] 吴军,贾雨.复杂介质的三维块状模型快速射线追踪[J].物探化探计算技术,2008,30(6):465.

[6] 张建中,陈世军,徐初伟.动态网络最短路径射线追踪[J].地球物理学报,2004,47(5):146.

[7] MOSER T J.Shortest path calculation of seismic rays[J].Geophysics,1991,56(1):59.

[8] 孔选林,李录明,罗省贤,等.各向异性介质中地震波射线正演[J].物探化探计算技术,2008,30(3):178.

[9] 张东,谢宝莲,杨艳,等.一种改进的线性走时插值射线追踪算法[J].地球物理学报,2009,52(1):200.

[10]高尔根,徐果明,蒋先艺,等.三维结构下逐段迭代射线追踪方法[J].石油地球物理勘探,2002,37(1):11.

[11]涂齐催,刘怀山.利用线性旅行时插值射线追踪计算地表模型初至波走时[J].2006,30(2):148.

[12]MOSER T J,VANECK T,NOLET G.Hypocenter determination in strongly heterogeneous earth models using the shortest path method[J].J.Geophys.Res,1992,97:6563.

[13]ASAKAWA E,KAWANAKA T.Seismic ray tracing using linear traveltime interpolation[J].Geophysical Prospecting,1993,41(1):99.