基于快速步进法的三维反射波追踪

姜莉

摘 要： 引入了一种改进的快速步进算法，对三维复杂介质中反射波进行波前追踪，并利用QT和OpenGL相结合的方法对速度场以及射线路径进行可视化模拟。通过模型仿真，验证了该方法的效率、精度以及可靠性。

关键词： 快速步进法； 反射波； 射线追踪； 波前追踪

中图分类号： TN911.7?34； P314.3 文献标识码： A 文章编号： 1004?373X（2014）11?0132?03

Abstract： In the paper， an improved fast stepping algorithm is introduced to make wavefront tracing of the reflected wave in 3D heterogeneous media， and perform the visualization simulation of velocity field and ray path by using QT and OpenGL. The efficiency， accuracy and reliability of the newalgorithm were verified by model simulation.

Keywords： fast stepping method； reflected wave； ray tracing； wavefront tracing

费马原理指出，地震波沿着所需时间最短的路径传播，在高频近似条件下，地震波场的主能量沿射线轨迹传播。地震射线追踪的目的就是寻找震源到各个接收点的地震波传播的最短传播路径，从而研究地震波的传播规律，对地层探测及层析成像具有重要的意义。针对传统射线追踪方法仅能求出地震波的初至旅行时的缺点，本文在快速步进算法的基础上给出一种多级的方法，对反射波进行射线追踪。该方法将波前到达的每一层分割为一个独立计算的区域，把波前到达的每层界面重新作为波源重新计算，从而逐层推进式地追踪反射波或透射波的波前。

1 算法实现

快速步进算法的主要思想就是通过系统地用上风格式构建旅行时场来解程函方程。程函方程的上风差分格式实际上意味着信息传播是从走时的小值到大值。因此，快速步进算法重点在通过从[T（x，y，z）]的小值向外扩展来构建方程的解[1]。

1.1 旅行时场方程

在3D旅行时计算时，程函方程为：

1.2 模型建立

把目标区域离散为网格，将没被界面穿过的网格称作规则网格，每个网格顶点称作规则结点，如图1中大实心点所示。在规则网格部分，每一个网格都是具有八个网格顶点的区域，同时每个规则结点属于八个网格。如果规则网格被界面切割，称作不规则网格，网格边与界面相交，产生的交点称作不规则网格结点，如图1中小实心点所示，每个交点属于四个网格。如果要更新某一点的走时，认为每个网格相互独立。

1.3 算法实现

如图2所示，为了加快计算，把目标区域的所有网格结点划分为三个子集，首先把具有初始旅行时值的点记为确认值点（实心黑点），这些点的旅行时在计算中不再改变；然后把与确认值点相邻一个网格的点记为近点（实心绿点），近点需要计算一个旅行时的值，但它不一定是最小值，是可改变的；其他点记为远点（空心点），远点的旅行时值还没有计算，需要逐次推进计算。

计算最小旅行时值的过程如下：首先在近点中搜索旅行时最小的点，把这个点从近点子集中去除，加入确认值点子集，然后把此点的所有非确认值的临近点归入近点子集，如果临近点属于远点子集，则从中去除；接着重新计算该点临近点的旅行时值，如果比原值大，则保留原值，如果小，则取而代之。重复前面的步骤，直到确认值点集包含整个目标区域。本算法利用堆排序技术，可有效地在近点集合中搜索最小元素，使算法加快速度。

1.4 反射波前追踪

对反射波射线追踪时如果直接利用快速步进法计算时，会存在以下问题：使用快速步进方法的有限差分解必须服从绝对初至时间。对于反射波，波前至少穿过计算区域两次，因此只能通过多级方法来解决。多级方法的原理如图3所示。

第一步，波前传播到反射界面，然后将该区域进行复制。

第二步，将快速步进方法用于复制的区域。第二步开始时，窄带围绕着界面，这样做并不破坏因果性，因为初至反射波前只到达每个界面结点一次。

1.5 走时更新

1.5.1 规则网格

对于规则网格点的走时更新，用上风格式近似表示式（1）的梯度函数，可得到如下的表达式：

1.5.2 不规则网格

如果被更新结点处于不规则网格中，使用相邻的不规则点的走时来更新某一结点时，假设波前近似为平面。当使用不规则网格中的确认值点更新时，此网格单元中的确认值点数可能是一个，两个，三个或大于三个。如果由三个确认值点的信息来更新此点走时，这三点可以惟一确定一个波阵面，准确计算该点走时。如果仅有一个或两个点是确认值点，必须增加额外的约束信息，找到波阵面的解。当确认值点多于三个时，将产生冗余信息。下面分别考虑这四种情况。

（1） 由三个确认值点更新

[x，][y，][z]分别为[x，][y]和[z]三个方向的单位矢量。由以上的计算公式可以看出，每个网格的慢度矢量在计算旅行时的过程中，就可求出，不用特别计算，从而可节省计算量。根据射线理论，射线与慢度矢量的方向一致，射线在模型中某一点的传播方向，用该点周围网格点的慢度矢量插值求出，在本文中采用双线性插值。

2 实验及仿真

本文采用层状介质模型进行实验仿真，模型为一个两层的三维介质模型。如图4所示，上层用单一颜色表示其传播速度为常数，下层的速度随深度变化，颜色越深表示传播速度越快。目标区域长宽深分别为2 000 m，2 000 m和1 000 m，在[x，][y]和[z]方向分别划分为81，81和41个网格，网格间距在各个方向上都是25 m，炮点位于地面上（0，1 000，1 000）坐标处。图4是模型速度场的剖面图，拖动右边滑块可以控制各方向剖面的移动。图5是射线及界面模拟图，拖动滑块可以对图像进行旋转。从图中可以看出，模拟结果与理论结果一致，快速步进法确实是一种快速稳定，适用于实际的方法。

3 结 论

本文采用Fortran和QT类库以及OpenGL相结合，的方法很好地实现了快速步进法射线追踪，不仅利用前人已有的成果进一步降低开发难度，提高开发效率。同时易于向Unix移植，以适应一些用户的特殊需要。由模型实验看，所得结果比较客观地反映出了实际情况。所以，此法是一种很好的三维正演方法。

参考文献

[1] KIM Seongjai. 3?Deikonalsolvers： first?arrival traveltimes [J]. Geophysical， 2002， 67（4）： 1225?1231.

[2] 许琨，吴律，王妙月.改进Moser法射线追踪[J].地球物理学进展，1998（4）：60?64.

[3] 张云姝，顾汉明，师学明.基于MOSER曲射线追踪的SIRT声波层析成像[J].工程地球物理学报，2005（3）：167?169.

[4] 张霖斌，刘迎曦，赵振峰，等.有限差分法射线追踪[J].石油地球物理勘探，1993，42（6）：198?200.

[5] 张美根，贾豫葛，王妙月，等.界面二次源波前扩展法全局最小走时射线追踪技术[J].地球物理学报，2006，49（4）：1169?1175.

[6] 李艳.地震勘探中二维波动方程的小波数值正演模拟[D].哈尔滨：哈尔滨工业大学，2002.

[7] 周竹生，刘喜亮，熊孝雨.弹性介质中瑞雷面波有限差分法正演模拟[J].地球物理学报，2007，50（2）：567?573.

[8] 刘恩儒，岳建华，刘彦.具有离散裂缝空间分布的二维固体中地震波传播的有限差分模拟[J].地球物理学报，2006，49（1）：180?188.

[9] 朱亚平，杨慧珠.OpenGL技术在地震数据可视化中的应用[J].石油地球物理勘探，2000，35（4）：403?414.

[10] WOO Mason， NEIDER Jaekie， DAVIS Tom， et al. OpenGL编程权威指南[M].北京：中国电力出版社，2001.

[11] 王博桦，陈振一，王关桥，等.地震信号检测新机理研究[J].现代电子技术，2013，36（21）：111?112.