0 引言

射线追踪是地震层析成像和一些偏移成像算法的关键。 射线追踪算法的精度和效率直接影响成像效果和计算效率。 有两类地质情况在实际的射线追踪工作中经常遇到,一是在起伏的地表面上做地震观测,二是在隧道或巷道超前地震勘探中,经常遇到地震波不能穿透的隧道或巷道、采空区、大型空洞等地质异常体。在隧道和煤矿巷道地震超前探测偏移成像处理中,成像区域近似全空间。在走时计算中,隧道或巷道本身作为一种非穿透体,以及采空区、空洞等的存在对地震波传播的影响不可忽略,因而会影响射线的弯曲和走时的准确计算。

起伏地形作为地震走时计算的边界处理,已经有不少学者研究过,但关于非穿透性地质体存在的情况下的走时计算方法鲜有报道。孙章庆等^[1-2]提出了两种复杂地表条件下的解决方案, 一种基于线性插值和窄带技术的地震波走时计算,另一种是不等距迎风差分法;这两个方法对起伏地表能够灵活处理,局部走时计算也具有较高的效率和稳定性。卢回忆等^[3]采用MSFM计算复杂近地表模型走时, 改善了标准FMM(fast marching method)存在对角方向误差大的缺陷。

从另一个角度来看,起伏地表模型边界实际上也是更大规模的非穿透性地质异常体的边界,如果能研究一种射线追踪方法能同时解决起伏边界与含非穿透性地质异常体的介质模型中的走时计算问题,将具有一定的实用价值。

传统的射线追踪方法有试射法和弯曲射线法^[4-6]。这两种方法在早期的地震勘探中被广泛应用, 但在局部非常复杂的介质模型中容易陷入局部最小值, 并且计算精度和效率都较低^[7]。Vidale^[8]最早提出了利用有限差分技术求解程函方程的射线追踪方法:通过矩形盒式扩展的方式递推计算整个区域的旅行时, 进而进行初至波追踪。 这种方法在介质局部速度变化强烈时算法不稳定, 且其扩展方式不符合波前传播规律,无法解决反射波及多值走时问题^[9]。 Podvin^[10]对Vidale方法进行了改进, 通过分析各计算网格点可能出现的不同的波传播形式, 确定旅行时面, 提高了旅行时计算精度, 但增加了计算量。

Sethian^[11]于1996年基于Level Set方法(Osher S.,et al,1988)提出了FMM方法。该方法遵循波前传播的熵守恒理论,采用有限差分迎风格式求解程函方程, 在实现过程中引入窄带技术和堆排序技术。Rawlinson和Sambridge^[12-13]提出了分区多步FMM界面处理方法, 实现了复杂层状介质反射波和多次波追踪计算。为了克服FMM方法在误差分布上对角线方向误差较大的缺点, Hassouna等^[14]提出了MSFM算法(multi-stencils fast marching methods)。改进的算法相较于标准FMM其计算精度显著提高。目前,FMM方法是公认效率最高、精度最高、且对任意复杂模型无条件稳定的射线追踪算法之一,已经被广泛用于地震层析成像和偏移成像的走时计算中^[3]。

本文主要研究基于FMM的含非穿透体和复杂地表模型走时计算。新算法在原始FMM算法的基础上, 通过对异常体区域和地表以上区域进行标记,当FMM窄带区计算到异常体区域时波前计算自动避开该区域,同时解决了非穿透体和起伏地表的走时计算问题。为此设计了一个二维规则模型进行旅行时计算,与理论时间进行比较, 分析了算法的精度问题。 最后用一个三维复杂模型展示了走时计算的实际效果。

1 常规FMM算法原理

在高频近似条件下,地震波传播时间满足程函方程:

|∇xT|=s(x),(1)

其中,∇_x表示梯度算子,T表示地震波的旅行时间,s(x)表示介质的慢度,是一个关于位置x的函数。 对式(1)近似程函方程中旅行时函数梯度项的迎风差分格式可写为^[15]:

max(Da-xT,-Db+xT,0)2+max(Dc-yT,-Dd+yT,0)2+max(De-zT,-Df+zT,0)2i,j,k12=si,j,k,(2)

其中,(i,j,k)表示直角坐标系(x,y,z)中网格点索引, a、b、c、d、e、f定义了迎风有限差分的阶数。 一般常用一阶或二阶梯度算子,例如:

x方向一阶梯度算子为:

D1-xTi,j,k=Ti,j,k-Ti-1,j,kΔx,D1+xTi,j,k=Ti+1,j,k-Ti,j,kΔx。(3)

x方向二阶梯度算子为:

D2-xTi,j,k=3Ti,j,k-4Ti-1,j,k+Ti-2,j,k2Δx,D2+xTi,j,k=4Ti+1,j,k-3Ti,j,k+Ti+2,j,k2Δx。(4)

y和z方向一阶和二阶算子具有与x方向类似结构, 可类比得出。

式(2)中各子项在计算时,具体采用的阶数根据计算的精度要求和网格节点计算所允许的阶数确定。 二阶精度比一阶要高。

FMM实现过程中采用窄带技术^[12-13] ,如图1所示,该技术将介质空间中所有网格点分为3种类型(用M标记网格点的类型) ^[7]:

1)近点(M=2),表示已经计算过的波前传播时间,在之后的计算过程中走时不再被更新的节点;

2)窄带点(M=1),表示计算过波前传播时间,但是波前传播时间可能会被更新的节点;

3)远点(M=0),表示尚未计算波前传播时间的节点。

由于震源附近波前面曲率非常大,走时计算在震源附近网格点会产生较大误差。 一种行之有效的方法是采用震源附近网格点加密策略^[12]。

图2给出了震源附近加密网格的波前推进以及映射到原始网格的波前。 首先以震源所在网格为中心确定合适的网格范围,在该网格范围内再将网格等间距细分。 图2a中选择震源周围3个网格点作为细分区域,每个网格再细分3个网格点。 然后计算细分网格走时,当计算到细分网格边界时, 把细分网格映射到原始粗网格中(图2b),再执行粗网格波前推进。

FMM实现流程如下:

步骤1:确定震源附近细分区域网格数,并细分网格。先在细分区域内计算走时。

步骤2:将震源点标记为近点(M=2),其他点标记为远点(M=0)。

步骤3:从震源点出发,利用式(2)计算震源周围点的波前传播时间,并将其表示为窄带点(M=1)。

步骤4:遍历窄带区,在其中找具有最小波前传播时间的点,将之标记为近点(M=2);

步骤5:将新的近点周围的远点标记为窄带点(M=1),并通过式(2)计算该点周围窄带点走时,如果得到的走时比原值小则更新该点,反之则保留。

步骤6:当计算到细分区域边界时,把细分区域走时映射到对应粗网格点上,得到新的窄带区。开始粗网格点波前推进。

步骤7:粗网格点波前推进过程和细分网格相同,即重复步骤4和步骤5,直到窄带区为空,最终计算出所有网格点走时。

2 含非穿透体及起伏地表FMM算法

在实际中,常常会碰到一些较为复杂的地质体, 如地下岩溶、空洞、煤矿中的巷道、采空区等。在这些异常体中地震波传播的速度为0,即非穿透体,程函方程在这些点处无意义。因而,在FMM波前推进过程中需要避开这些点的计算。地表以上区域网格点也不需要计算。实际上,复杂的地表区域也可以看作非穿透体的边界,因而,可以采用同一种方案同时处理这两种情况。

真实地表和非穿透体边界并不一定落在规则网格节点处,此时可以采用离边界最近的网格点近似代替边界。图3给出了含非穿透体和起伏地表的波前推进过程示意图。当波前推进的窄带区到达非穿透体时(图3a),窄带需要分裂,并沿各自的方向行进(图3b),当绕过异常体后,窄带区又会合并成一个整体并继续推进(图3c)。窄带区一端沿着地表推进。

对非穿透体和地表区域我们采用标记的方法,将其标记出来不参与旅行时计算。 因此,模型空间中的网格点可分为4种类型(用M标记网格点的类型):

1)异常点(M=3),表示不参与波前推进计算的点;

2)近点(M=2),表示已经计算过波前传播时间,在之后的计算过程中旅行时不再被更新的节点;

3)窄带点(M=1),表示计算过波前传播时间,但是波前传播时间可能会被更新的节点;

4)远点(M=0),表示尚未计算波前传播时间的节点。

含非穿透体及起伏地表FMM实现流程如下:

步骤1:确定震源附近细分区域网格数,并细分网格。 先在细分区域内计算旅行时。

步骤2:将异常体和地表以上区域所在网格点标记为3(M=3),将震源点标记为近点(M=2),其他点标记为远点(M=0)。

步骤3:从震源点出发,判断下一网格点标志,当M等于0时,利用式(2)计算波前传播时间, 并将其表示为窄带点(M=1);当M等于3时,跳过该点,不计算。

步骤4:遍历窄带区,在其中找具有最小波前传播时间的点,将之标记为近点(M=2);

步骤5:将新的近点周围的远点标记为窄带点(M=1),并通过式(2)计算该点周围窄带点旅行时,若得到的旅行时比原值小则更新该点,反之则保留。

步骤6:当计算到细分区域边界时,把细分区域旅行时映射到对应粗网格点上,得到新的窄带区,开始粗网格点波前推进。

步骤7:粗网格点波前推进过程和细分网格相同,即重复步骤4和步骤5,直到窄带区为空,最终计算出所有网格点走时。

3 数值算例

为了验证本文算法计算含非穿透体及地表模型的有效性,下面以一个二维和一个三维模型说明走时计算的效果。

3.1 模型1

模型1是一个含规则非穿透体二维均匀模型(图4)。 模型大小为50 m×50 m,背景速度1 000 m/s, 网格间距0.25 m。 非穿透体为正方形, 位于模型x和z方向均为20~30 m范围内。 震源位于点(25 m,0 m)。 采用一阶和二阶差分FMM进行了计算,图4是计算结果及其误差分布,图中等值线表示地震波的传播时间分布。 为了说明计算结果的准确性,图中标出了地震波传播的理论时间(图中虚线)。 由于模型速度为常量,模型各网格点理论时间可以通过传播距离和模型速度很容易直接计算得到。

由图4可以看出,计算时间场与理论时间场等值线吻合较好,说明算法对于非穿透体模型旅行时计算是正确的。 其中图4a为一阶FMM计算旅行时,其与理论值存在一定误差,但二阶FMM计算值与理论值等值线几乎重合(图4b),显然二阶FMM计算精度更高。 从误差结果来看(图4c和4d),二阶FMM误差分布总体比一阶FMM误差分布小很多, 总体均方根误差一阶是二阶的大约7倍(表1)。 无论是一阶还是二阶FMM其误差分布在震源点水平方向和垂直方向上误差都较小,45°方向误差较大。 当波场绕射到非穿透体后面, 绕射点形成新的震源, 也产生45°方向较大误差。 FMM算法误差分布存在明显的方向特性,即对角线方向较大。 本文算法的计算精度和误差分布规律与常规FMM保持一致。 计算效率上,二阶FMM计算时间略微比一阶长(表1)。 图5给出了震源到模型底端检波器各个位置的射线路径,由于模型背景速度为常量,射线为直线,但碰到非穿透体时,射线发生改变,射线计算结果符合实际情况。

3.2 模型2

模型2是一个三维复杂模型,包含一个隧道和两个空洞,并具有起伏地表的超前探测模型(图6)。

模型大小为40 m×20 m×20 m,隧道沿x方向掘进, 长15 m,宽和高均为5 m,其x方向中心线位于y=10 m,z=12.5 m处。在隧道前方存在两个不规则空洞(图6a)。模型上部为起伏地表面。震源位于隧道正前方(16 m,10 m,12.5 m)处,计算采用的网格间距为0.1 m。在起伏地表的处理上,把地表上方的区域也当作不可穿透介质。计算结果如图6b和6c所示。通过计算结果发现, 当地震波前碰到隧道、空洞,波场会自动绕射过去,符合地震波在地下传播的规律。同时,算法也能够解决起伏地表的走时计算问题,表现出了较好的适应性。此外,模型的计算结果也说明算法保持了常规FMM的稳定性和高效性。

4 结论

在原始FMM算法基础上,发展了含非穿透体和起伏地表的复杂模型地震波走时计算方法。理论研究和数值试验表明:①新算法采用标记非穿透体区域的方法, 当窄带区碰到非穿透异常体时,会自动避开或绕过异常体。新算法只需在原始FMM算法基础上改变较少量的代码即可实现,较为简单;②起伏地表上方的区域可以当作非穿透体对待,因而新算法同时适用于起伏地表模型,表现出较强的适应性;③新算法保持了原始FMM算法精度和效率,误差分布也与原算法类似。

(本文编辑:叶佩)

The authors have declared that no competing interests exist.

---

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子

[1]	孙章庆,孙建国,韩复兴. 复杂地表条件下基于线性插值和窄带技术的地震波走时计算 [J].地球物理学报,2009,52(11):2846-2853. [本文引用: 1]
[2]	孙章庆,孙建国,韩复兴. 三维起伏地表条件下地震波走时计算的不等距迎风差分法 [J].地球物理学报,2012,55(7):2441-2449. [本文引用: 1]
[3]	卢回忆,刘伊克,常旭. 基于MSFM的复杂近地表模型走时计算 [J].地球物理学报,2013,56(9):3100-3108. [本文引用: 2]
[4]	Julian B R,Gubbins D. Three-dimensional seismic ray tracing [J].Journal of Geophysics Research,1977,43(1):95-114. [本文引用: 1]
[5]	Um J,Thurber C. A fast algorithm for two-point seismic ray tracing [J].Bulletin of the Seismological Society of America,1987,77(3):972-986.
[6]	Prothero W A,Taylor W J,Eickemeyer J A. A fast,two-point,three-dimensional ray tracing algorithm using a simple step search method [J].Bulletin of the Seismological Society of America,1988,78(3):1190-1198. [本文引用: 1]
[7]	王飞,曲昕馨,刘四新,等. 一种新的基于多模板快速推进算法和最速下降法的射线追踪方法 [J].石油地球物理勘探,2014,49(6):1106-1114. [本文引用: 2]
[8]	Vidale J. Finite-difference calculation of travel-times [J].Bulletin of the Seismological Society of America,1988,78(6):2062-2076. [本文引用: 1]
[9]	李永博. VTI介质及复杂模型FMM射线追踪方法研究[D] .西安:长安大学,2012. [本文引用: 1]
[10]	Podvin P,Lecomte I. Finite difference computation of traveltimes in very contrasted velocity models:A massively parallel approach and its associated tools [J].Geophysical Journal International,1991,105(1):271-284. [本文引用: 1]
[11]	Sethian J A. A fast marching level set method for monotonically advancing fronts [J].Proc. Nat. Acad. Sci.,1996,93:1591-1595. [本文引用: 1]
[12]	Rawlinson N,Sambridge M. Multiple reflection and transmission phases in complex layered media using multistage fast marching method [J].Geophysics,2004,69:1338-1350. [本文引用: 3]
[13]	Rawlinson N,Sambridge M. Wave front evolution in strongly heterogeneous layered media using the fast marching method [J].Geohys. J. Int.,2004,156:631-647. [本文引用: 4]
[14]	Hassouna M S,Farag A A. Multistencils fast marching methods:A highly accurate solution to the eikonal equation on cartesian domains [J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2007,29(9):1563-1574. [本文引用: 1]
[15]	Sethian J A.Level set methods and fast marching methods[M].Cambridge:Cambridge University Press,1999. [本文引用: 1]