基于波场上传试射的三维射线追踪方法

辛全圣，贾 雨，舒 引

(1.成都理工大学核技术与自动化工程学院，成都 610059;2.电子科技大学通信与信息工程学院，成都 611731)

引 言

射线追踪方法在地震勘探中有广泛应用，它在地震定位、地球物理层析成像、地震勘探资料的偏移处理、地震模拟及其它地震数字处理与反演中都起着极重要的作用。传统的射线追踪方法有试射法［1］和弯曲法［2］，后来发展起来的有，有限差分法［3-4］、旅行时线性插值法［5-7］、最短路径法［8-9］、波前重建方法［10］等。与其他方法相比，试射法精确性高，在全面搜索和适应复杂地质模型等方面有很大优势，是两点追踪问题的一种重要方法［11］。

试射法，又称打靶法，是基于斯奈尔(Snell)定理的计算方法。试射法是最早提出和使用的射线追踪方法，在数学上属于初值问题。试射法根据由震源出发的一束射线到达接收点的情况对射线出射角及其密度进行调整，最后由最靠近接收点的两条射线走时内插求出接收点走时［12］，从而实现两点间的射线追踪计算。试射法虽然是一种经典的射线追踪方法，但不能误认为它是低级的、行将淘汰的算法，相反，一直在不断的发展。试射法在部分问题上仍需进一步改进:其一是在计算精度方面存在难以解决每一个炮检对存在多条射线的问题［其二是在计算效率方面存在每一次观测系统变化都需要进行试射和迭代的完整流程而导致计算效率较慢的问题。为此，本文对CRP统计数据的准确度和三维射线追踪的计算效率进行改进。

1 方法原理

通过试射获得可能存在射线的区域，通过迭代计算准确的射线轨迹。首先以目标反射层为基准，计算反射层的每一个控制点到地表的全方位试射射线，并作为迭代第一步的初始数据，只要模型不变化，则试射结果一直有效;然后给定观测系统，对控制点到每一个炮点的射线进行迭代，并记录有效射线作为迭代第二步的初始数据，只要观测系统的炮位置不变化，则炮点到目标层控制点计算结果一直有效;最后遍历每个目标层三角形，迭代出该反射点在该三角形内部的所有射线。

1.1 试射

传统的试射方法主要采用从炮点进行试射，水平角度的范围为0到360度，垂直角度范围为0到180度，往反射层试射。到达反射层之前通过和子面求交点，然后发生折射，如果折射不成功，则这个角度试射失败，用下一角度继续进行试射。当射线到达反射层时，根据反射定理，利用当前方向和反射点的法向求得反射方向，从反射层往地表计算射线，简单的二维示意图如图1所示。本文采用的试射方法对传统的试射方法做了相应改进，本文采用的试射方法从反射层的反射点进行试射，水平角度的范围为0到180度，垂直角度范围为0到90度。简单的二维示意图如图2所示。

图1 从炮点进行试射二维示意图

图2 从反射点进行试射二维示意图

1.2 迭代

迭代部分是把试射部分的结果作为原始数据，通过三角形面积迭代法最终找出炮点到检波点的有效射线。先根据试射文件求出每一个炮点到反射层每一个顶点经反射后再到接收层的射线，记录文件，如图3所示。然后根据上一步求出的射线找出炮点经反射层反射到其关联检波器的射线，即为射线追踪最终射线结果，如图4所示。

1.3 三角形面积迭代法

图3 迭代第一步示意图

图4 迭代第二步示意图

如图5所示，已知目标点P和所在的出射三角形T1T2T3，面积为S，P将T1T2T3分成三个三角形，面积分别为S1，S2，S3。三个顶点对应的射线的方向向量分别为点在三角形中的面积坐标加权求和得出新的方向分量以该方向分量试射并得到出射点T，如果出射点T和接收器的距离小于给定的精度，则试射迭代结束。否则，用出射点T把原来的三角形分割成三个新的三角形 TT1T2、TT2T3和TT1T3。找到P点所在三角形，并作为下次试射迭代的初始三角形。按上述方法继续试射直到出射点和接收器的距离小于给定的精度为止。

图5 三角形面积迭代法示意图

2 基于波场上传试射算法

2.1 试射

试射是射线追踪核心算法的第一步，试射的结果以射线的格式存放文件，一条射线记录了该射线的方向、反射点坐标、射线路径中和每一个层面的交点的坐标和法向、以及接收层的入射点和出射点。试射的第一步是确定试射角度。理论上试射射线在接收层越密集，即角度间隔越小，射线追踪的结果越准确，但是角度间隔越小，试射部分所花的时间就越多，可根据工区的实际情况选择合适的角度间隔。

2.2 迭代第一步

(1)首先初始化观测系统，读取所有炮点和检波点，获得炮点和检波点的x值和y值，通过向接收层投影，给炮点和检波点赋z值。

(2)逐一处理每一个反射点，首先读取反射点的射线文件，根据每条射线的出射点x，y坐标，将该射线投入网格中，利用一个二维的vector变量存放对应射线的索引。

(3)遍历处理每一个炮点，首先根据炮点的x，y坐标找到炮点属于的网格，然后根据二维的vector变量数据，取出该网格以及周围一圈网格中的射线索引，通过这些射线索引获取每一条射线，从这些射线中选出三条出射点包围该炮点的射线，如果存在多条，则选出构成的三角形面积最小的三条射线。如图3所示，三条紫色射线r1，r2，r3为选出的三条满足要求的射线。

(4)根据找出的三条射线和炮点利用三角形面积迭代法找出炮点到反射点的射线r，如果迭代过程中不收敛，即某一个方向的试射射线的出射点不在迭代的三角形中，则返回失败。

2.3 迭代第二步

(1)对每一个炮点逐一处理，读取该炮点到反射点经反射后再到接收层的所有射线。

(2)遍历反射层每一个三角形，取出三角形三个反射点，判断当前炮点到三个反射点是否都有射线，如果到某个点没有射线，则跳过当前三角形，处理下一个三角形。

(3)如果三角形三个顶点都有射线，则取出三条射线。如图4所示，r1，r2，r3为炮点到反射层当前三角形三个顶点的三条射线，三条射线在接收层的出射点构成的三角形为T1T2T3。遍历当前炮点的关联检波点，如果检波点在三角形内部，则用三角形面积迭代法找出最终射线r。需要注意的是，这里三角形迭代法找出最终射线的比例比迭代第一步找出射线的比例低，因为第一步中是利用三角形面积迭代法找出反射点到炮点的射线，而这一步是找出从炮点到检波点经过反射层反射的射线，随着折射反射越多，收敛比例越低。

3 数值模拟及分析

本文通过对复杂地质构造理论模型进行射线路径的计算，验证该算法的有效性。图6为射线路径图，炮点为2000个，25条线，每条线80个炮点，每炮关联检波器为1152个，红色排列片代表接收的检波点，绿色的点代表未接收的检波点。图7为水平界面的理论CMP面元情况与实际射线追踪统计的CRP面元情况对比图。可以看出，当界面为水平界面时，CRP方位角和炮检距统计图与基于CMP理论计算的完全一致，且不同偏移距范围的覆盖次数呈规则分布。CRP位于炮检点的中点，与理论CMP重合，射线路径对称，CRP覆盖次数分布规律与理论计算的覆盖次数分布规律完全一致。可以验证CRP统计数据的准确度明显提高。以冯家湾三维工区为例，试射单线程单机计算估算15天左右，迭代第一步计算量为试射的1/5左右，迭代第二步相对于前两步计算量很小，估计少一个数量级。计算效率的提高取决于实际的工区和观测系统，数据量越大、观测系统变化的次数越多，计算效率提高越明显。传统的试射方法每一次观测系统变化都需要进行试射和迭代的完整流程。本文的试射方法只要模型不变化，则试射结果一直有效;只要观测系统的炮位置不变化，则炮点到目标层控制点计算结果一直有效，提高了三维射线追踪的计算效率。

图6 射线路径图

图7 水平界面的理论CMP面元情况与实际射线追踪统计的CRP面元情况对比图

4 结束语

传统的试射方法从炮点出发，通过调整入射角使得射线逐步逼近检波点位置，有些检波点无法接收到射线，存在盲区。针对传统的试射方法难以解决每一个炮检对存在多条射线以及每一次观测系统变化都需要进行试射和迭代的完整流程而导致计算效率较慢的问题，本文提出了基于波场上传试射的三维射线追踪方法，从原理上解决了每一个炮检对存在多条射线的问题。同时，只要模型不变化，则试射结果一直有效;只要观测系统的炮位置不变化，则炮点到目标层控制点计算结果一直有效。数值模拟结果表明，该算法能够提高CRP统计数据的准确度，根据实际工区和观测系统的不同，计算效率有不同程度的提高，适于复杂地质构造三维射线追踪的工程应用。

［1］马争鸣，李衍达.二步法射线追踪［J］.地球物理学报，1991，34(4):501-508.

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

［3］ Vidale J E.Finite-difference calculation of traveltimes［J］.Bulletin of the seismological society of A-merica，1988，78(6):2067-2076.

［4］ Vidale J E.Finite-difference calculation of traveltimes in three dimension［J］.Geophysics，1990，55:521-526.

［5］张东，张婷婷，乔友锋，等.三维旅行时场B样条插值射线追踪方法［J］.石油地球物理勘探，2013，48(4):559-566.

［6］ Zhang D，Zhang T T，Zhang X L.A new 3-D ray tracing method based on LTI using successive partitioning of cell interfacesand traveltime gradients［J］.Journalof Applied Geophysics，2013，92:20-29.

［7］叶佩，李庆春.旅行时线性插值射线追踪提高计算精度和效率的改进方法［J］.吉林大学学报，2013，43(1):291-298.

［8］ Moser T J.Shortest path calculation of seismic rays［J］.Geophysics，1991，56(1):59-67.

［9］桑运云，李振春，张凯.抛物旅行时插值最短路径射线追踪［J］.石油地球物理勘探，2013，48(3):403-409.

［10］王庆林，白超英.相空间波前构造法地震多值射线追踪的改进［J］.石油地球物理勘探，2013，48(1):15-21.

［11］徐 涛，徐果明，高尔根，等.三维复杂介质的块状建模和试射射线追踪［J］.地球物理学报，2004，47(6):1118-1126.

［12］张钋，刘 洪，李幼铭.射线追踪方法的发展现状［J］.地球物理学进展，2000，15(1):36-45.