三维层状介质中基于走时梯度的多次波射线追踪

三维层状介质中基于走时梯度的多次波射线追踪张婷婷;张东;邱达【摘要】本文结合分区多步思想改进初至波走时梯度射线追踪算法(MTG),使之能适用于三维层状介质中多次波射线追踪计算.在计算过程中按照射线类型和速度分界面将三维模型分为不同的计算区域,根据多次波的传播顺序从震源开始在对应的区域内采用网格剖分线性插值法计算波前走时场.然后从检波点开始根据多次波的类型在对应的计算区域内沿走时梯度方向分步追踪多次波射线路径.在计算多次波波前走时场的过程中,为提高计算效率加入了逐次网格剖分技术,使分界面上波前新的震源位置更加精确;在多次波射线追踪的过程中根据走时梯度方向追踪射线路径,在速度分界面区域内采用局部追踪策略求取下一个追踪点以提高多次波的射线精度.本文在多层层状介质模型及扩展到三维的Marmousi模型中进行试验,计算结果表明,改进的MTG算法结合分区多步技术计算多次波时具有常规MTG法计算初至波时的优点,计算精度高、算法简单、计算效率高,验证了MTG算法在复杂的三维模型中计算多次波的有效性.【期刊名称】《石油地球物理勘探》【年(卷),期】2014(049)006【总页数】9页(P1097-1105)【关键词】多次波追踪;分区多步;走时梯度算法;线性插值【作者】张婷婷;张东;邱达【作者单位】武汉大学物理科学与技术学院,湖北武汉430072;湖北民族学院信息工程学院,湖北恩施445000;武汉大学物理科学与技术学院,湖北武汉430072;湖北民族学院信息工程学院,湖北恩施445000【正文语种】中文【中图分类】P6311 引言如今的地震勘探技术正朝着精细化方向发展。

以往主要应用反射波地震勘探，多次波被当作干扰波看待，实际上多次波也传递了许多地下信息，应被加以有效应用：在多次波走时联合层析成像［1，2］中，利用反射层析弥补了大炮检距的折射层析盲区，同时又能更真实地反映地下不同位置以及不同深度、尺度更精细的结构分布；在多次波的压制方法研究中，需要实用的多次波正演模拟方法；在多波偏移成像等领域中也有应用。

因此，多次波的追踪计算变得越来越重要［3，4］。

最初的射线追踪算法用于计算地震初至波走时，包括传统的试射法和弯曲法、有限差分法［5～10］、最短路径法［11～14］、线性插值法［15～19］以及快速推进法［20～23］。

关于后续波的射线追踪有：Moser［24］利用最短路径方法实现了二维模型中反射波追踪计算；Podvin等［25］、Hole等［26］以及 Riahi 等［27］利用有限差分法计算速度分界面上的反射波；赵爱华等［28］实现了PS 转换波的计算；王童奎等［29］利用二次界面源算法实现了PS转换波射线路径追踪；赵改善等［30］利用互换原理采用走时线性插值计算非零炮检距情况下多个反射界面、多个接收点的反射波。

为了有效追踪多次波，Rawlinson等［31］提出了基于快速推进法（Fast Marching Method）的多次波计算技术，在二维模型中实现了多次波射线路径的追踪计算。

之后De Kool等［32］又将该算法推广至三维模型，实现了三维速度介质模型下的多次波波前模拟和射线追踪。

唐小平等［33，34］将分区多步思想与改进的最短路径法结合，实现了二维、三维复杂层状介质中多次反射、透射及转换波的追踪计算。

Huang等［35］利用改进后的双线性插值波前快速推进法计算三维模型中不规则网格下的一次反射和透射波。

在上述多次波射线追踪算法中大多采用的是有限差分法或最短路径法。

在传统的有限差分法中只能求网格点上的走时，而传统的最短路径法中需要解决入射角覆盖问题，否则射线路径始终是由炮点到接收点之间一系列网格节点直接连接而成的，射线出现“Z”字型的路径，射线精度和计算效率往往受到网格大小的限制，并且使得射线反射角（或折射角）不能随着入射角而连续变化。

为了快速求解三维模型，同时又不损失精度，张东等［36，37］提出了三维模型下利用线性插值或B样条插值沿最大走时梯度反方向追踪的射线追踪算法（the Maximum Traveltime Gradient Ray Tracing Method，简称 MTG）。

该算法的优点是由于B样条插值能够得到连续光滑的走时场，因此搜索步长不受网格大小的限制，在不增加网格密度的情况下，可以很大程度消除传统射线追踪方法由于单位网格划分给追踪路径带来的误差。

这种走时梯度射线追踪算法到目前为止只用于初至波射线追踪的计算，在本文中将其应用到多次波的射线追踪时对其做了改进。

首先采用逐次网格剖分法保证多次波波前走时场的计算效率；其次在多次波追踪过程中在速度变化剧烈的区域采用局部线性插值求下一个次级源点，在整个多次波射线追踪过程中沿走时梯度方向结合分区多步思想采用基于B样条／线性联合插值的方法计算三维层状模型下的多次透射、转换以及反射波射线路径。

2 模型分区和多步计算本文的分区多步思想源于Rawlinson等［31］提出的快速推进法多次波射线追踪算法。

第一步是模型分区，根据界面起伏的情况以及后续波的类型将模型分为不同的计算区域。

在计算过程中，速度模型至少分为三个区域：其一为有效区，即当前射线经过的区域；其二为边界区域，即速度分界面；其三为等待区，即当前还没有计算、波前还没有到达的区域。

第二步进行多步射线追踪计算，即首先根据后续波的类型计算波前走时场，然后追踪计算多次反射波、透射波及转换波。

图1为多次波分区多步计算示意图。

在一个两层速度模型中，起始下行P1波（在这里P、S分别表示纵波、横波；数字表示计算区域，上标、下标分别表示上行波、下行波）从炮点开始时在上层介质中传播直到遇到速度界面为止，此时有效区域为上层介质（浅灰色区域），等待区域为第二层介质（黑色区域）；如果在反射界面上发生反射（或反射转换），表示为P1P1（或P1S1），此时的计算区域不变，从速度分界上的最小走时点开始计算上行的反射波；如果在反射界面上发生了透射（或透射转换），则表示为P1P2（或P1S2），此时计算区域变为第二层介质区域，第一层变为等待区。

图1 多次波分区多步计算示意图图2为多次波前传播分区多步计算示意图。

分区多步计算多次波射线具体步骤如下。

（1）从炮点所在的区域开始往下进行波前扫描计算走时场，当当前区域波前扫描计算结束时，波前停留在速度分界面（或边界）上等待下一步的命令，图2a为下行波走时场。

（2）将分界面上的最小走时节点作为波前新的震源，根据进入新的计算区域所追踪射线的类型（如反射、透射、转换波）继续进行走时场的计算。

（新调区域必须与原计算区域相邻，在追踪反射波或反射转换波时，此时新调区域就是原计算区域。

）如果地震波在分界面处发生反射，则从新的震源开始计算当前区域内的上行波，如图2b所示；如果反射时追踪的是转换波，则调用转换波的速度计算上行波；如果是计算透射波，则从分界面上的新震源开始计算下层介质中的下行波走时场，如图2c所示；同样，如果是透射转换波，在计算下行波场时调用对应转换波的速度。

（3）如果是复杂的多次波，同样采用步骤（1）、（2）计算组合就可以计算多次反射波、透射波以及相应的转换波走时场并保存。

（4）根据步骤（1）、（2）、（3）计算出波前走时场后，按照多次波的传播规律从检波点开始在之前计算得到的多个独立走时场中沿走时梯度负方向进行射线追踪，就可以追踪多次波的射线路径。

图2 多次波前传播分区多步计算示意图（a）下行波P1走时场；（b）反射上行波P1P1或P1S1走时场；（c）透射下行波P1P2或P1S2走时场3 改进后的MTG方法多次波射线追踪计算和初至波追踪一样包括两个步骤：向前处理计算波前走时场以及向后处理计算接收点到炮点射线路径。

在MTG方法中走时场的计算采用基于走时线性插值的网格剖分法。

图3为逐次网格剖分示意图。

网格界面剖分假设网格单元边界面上走时线性分布，网格边界面上任意点的走时可由相邻四个网格顶点走时线性插值得到，通过线性插值计算选取网格界面上的剖分点中满足最小走时的点作为入射点，如图3左所示。

由于篇幅的原因具体计算公式见文献［37，38］。

为了进一步提高计算效率，在找到具有最小走时的入射点后，在这个点附近继续添加剖分点阵，然后在这些次级剖分点上继续做插值计算如图3右图所示，如果遇到更小走时的点，则将该点作为新的入射点。

按照文献［17］中介绍的扫描顺序，就能快速有效地计算出有效区域内各节点波前的走时。

在此需指出模型中网格节点的个数并没有增加，只是在单个网格单元内计算时选取了不同的剖分点计算。

在MTG方法应用到多次波波前走时场的计算中时，在分区多步计算中采用逐次网格剖分法计算边界区域上波前走时，由于逐次剖分点的引入能够在界面上找出更加精确的最小走时点作为透射波、反射波波前新的震源，所以此法可提高计算精度和效率。

图3 逐次网格剖分示意图多次波射线追踪向后处理依据地震波传播路径垂直于波前面的原理在各个独立的波前走时场内，通过3次B样条插值计算走时场梯度方向即波前面的法线方向，从检波点出发在相应的计算区域内沿走时场的梯度负方向追踪，即可得到从震源到检波点的多次波射线路径。

利用B样条插值计算追踪点走时需要用到周围4×4×4个节点的走时数据，计算公式如下其中：Tlmn是网格节点（l，m，n）的走时；i，j，k为追踪点所在的网格单元；s＝x－i；t＝y－j；u＝z－k；Bl（s）、Bm（t）、Bn（u）为B样条基函数。

与追踪初至波不同的是在利用走时梯度追踪多次波时需要分两种情况。

（1）当追踪到速度分界面上时如果追踪的是反射波，那么下一个计算区域还是当前计算区域，此时的追踪点位于有效区域的边界处。

B样条插值需要利用待插值点周围27个网格上64个网格节点的走时数据作为控制点进行插值，而此时一部分数据点位于等待区内，而等待区内的节点不能参与计算。

当参与插值的数据点不足时，如图4所示，这时的走时插值公式为图4 追踪点位于有效区域边界时的情况（2）当地震波到达速度突变的分界面发生折射时，根据斯内尔定律，Sinα／Sinβ＝v1／v2（v1＜v2），分界面两侧的走时场并不是连续的。

而根据三次样条函数的定义［39］，在整个区间［a，b］上，三次插值样条函数S（x）为二阶连续可导函数，S（k）（xi－0）＝S（k）（xi＋0），k＝0，1，2，其中xi（i＝0，1，…，n）为样条函数S（x）的节点。

而介质Ⅰ与介质Ⅱ的走时场在分界面上的两侧并不是连续的，不满足样条函数使用的条件，不能继续利用B样条插值公式计算走时。

在这种情况下采用局部射线追踪策略，即在射线追踪过程中沿梯度方向每向前追踪一个步长都判断追踪点是否到达速度突变区域。

当射线追踪到速度突变区域时，与线性插值计算走时场时类似，假设网格单元边界面上走时呈线性分布，追踪点所在网格中的界面上所有的点都有可能是次级震源，将网格界面进行剖分通过线性插值计算选取网格界面上的剖分点中满足最小走时的点作为次级源点。