一种复杂介质条件下射线追踪方法与流程

本发明属于陆上或海洋地震勘探领域，具体涉及一种复杂介质条件下射线追踪方法。

背景技术：

随着地震勘探程度的不断加深，地震勘探的难度也在不断加大，为了获取复杂地质条件下精细的地层构造和物性参数，往往需要开展地震反演的研究工作。而地震层析反演可以利用检波器接收的初至波信息求取介质内部的构造形态和物性参数，能够确定地球内部的局部不均匀性和精细结构，为全波形反演和地质解译奠定基础。传统的地震层析反演需要用到震源点和检波点之间的射线传播路径，而求取射线路径的方法包含基于初值问题的试射法和基于边值问题的弯曲法。在速度非均匀性较弱的情况下，这两种方法都能取得较好的效果。但是，在具有较强速度变化的介质条件下，试射法求取的射线路径稳定性较差且存在阴影区；而弯曲法射线追踪效率较低，且不能保证求取的走时为全局最小解。因此，传统的试射法和弯曲法无法实现复杂介质条件下射线路径的精确求取，严重制约了地震层析反演在实际应用中的效果。

因此,仍需要研究设计一种复杂介质条件下射线路径的精确求取方法,已解决上述现有技术不足。

技术实现要素：

本发明的目的针对现有技术而不足，提供一种复杂介质条件下射线追踪方法，用以实现复杂介质条件下射线路径的高效精确求解。

本发明的技术方案是：

一种复杂介质条件下射线追踪方法，包括如下步骤：

步骤一：对介质的慢度模型进行差分网格离散剖分，设定每个单元网格内部为常慢度，单元网格交点位置为待求取的初至波旅行时；

步骤二：初始化震源点位置；利用自适应程函方程求解方法求取单元网格交点位置的初至波旅行时t(i,j)，其中，i和j分别代表了x和z方向离散化后的网格节点坐标；

步骤三：基于初至波旅行时，利用中心差分公式，求取每个单元网格内部的入射角θ(i,j)；

步骤四：初始化检波点位置信息和入射角方向；

步骤五：从检波点出发，根据射线入射点位置信息和单元网格内的入射角方向，计算当前单元网格内的射线出射点坐标和进入的下个单元网格信息；

步骤六：将计算的当前网格射线出射点坐标作为进入的下一个单元网格的入射点坐标，并利用入射方向(flag)计算下一个单元网格的出射点坐标和之后进入的单元网格信息；

步骤七：重复步骤五和六，直至射线出射点坐标到达震源位置。

所述步骤三中每个单元网格内部的入射角θ(i,j)求取公式为tz(i,j)＝((t(i+1,j+1)+t(i+1,j))-(t(i,j+1)+t(i,j)))/(2dz)

tx(i,j)＝((t(i+1,j+1)+t(i,j+1))-(t(i+1,j)+t(i,j)))/(2dx)

其中，dz和dx分别代表了z方向和x方向的离散采样间隔，tz(i,j)和tx(i,j)分别代表了初至波旅行时间场沿z方向和x方向的方向导数，atan函数为反正切函数。

所述步骤四初始化检波点位置信息和入射角方向还包括，初始化公式，其初始化公式为：

其中，rz和rx分别代表了检波点在z方向和x方向的坐标信息；bdy_up、bdy_bm、bdy_lt和bdy_rt分别代表了慢度模型上、下、左和右的边界坐标信息；flag代表了入射角方向的选择。

所述步骤五从检波点出发，根据射线入射点位置信息和单元网格内的入射角方向，计算当前单元网格内的射线出射点坐标和进入的下个单元网格信息还包括：当入射角方向flag＝1时，射线从网格左上侧端点进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<θ<π/4时，出射点位置坐标(x1+h，z1+h*tan(θ))，进入的下一个单元网格信息是i+1，j，flag＝8；

当θ＝π/4时，出射点位置坐标(x1+h，z1+h)；进入的下一个单元网格信息是i+1，j+1，flag＝1；

当π/4<θ<π/2时，出射点位置坐标(x1+h*tan(π/2-θ))；进入的下一个单元网格信息是i，j+1，flag＝2；

当π/2<＝θ<＝π时，出射点位置坐标(x1，z1+h)；进入的下一个单元网格信息是i-1，j+1，flag＝3；

当π<θ<3π/2时，出射点位置坐标无；进入的下一个单元网格信息无；

当3π/2<＝θ<＝2π时，出射点位置坐标(x1+h，z1)；进入的下一个单元网格信息为i+1，j-1，flag＝7。

其中z1和x1分别代表了入射点z方向和x方向的坐标，h代表了单元网格的剖分长度。

所述步骤五还包括：当入射角方向flag＝2时，射线从网格上侧进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<θ<atan(h/(i*h-x1))时，出射点位置坐标(i*h，z1+(i*h-x1)*tan(θ))；进入的下一个单元网格信息为i+1，j，flag＝8；

当θ＝atan(h/(i*h-x1))时，出射点位置坐标(i*h，z1+h)；进入的下一个单元网格信息为i+1，j+1，flag＝1；

当atan(h/(i*h-x1))<θ<π-atan(h/(x1-(i-1)*h))时；出射点位置坐标(x1+h*tan(π/2-θ)，z1+h)进入的下一个单元网格信息为i，j+1，flag＝2；

当θ＝π-atan(h/(x1-(i-1)*h))时；出射点位置坐标((i-1)*h，z1+h)进入的下一个单元网格信息为i-1，j+1，flag＝3；

当π-atan(h/(x1-(i-1)*h))<θ<π时；出射点位置坐标((i-1)*h，z1+(x1-(i-1)*h)tan(π-θ))进入的下一个单元网格信息为i-1，j，flag＝4；

当π<＝θ<3π/2时；出射点位置坐标((i-1)*h，z1)进入的下一个单元网格信息为i-1，j-1，flag＝5；

当3π/2<θ<＝2π时；出射点位置坐标(i*h，z1)进入的下一个单元网格信息为i+1，j-1，flag＝7。

所述步骤五还包括：当入射角方向flag＝3时，射线从网格右上侧端点进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<＝θ<＝π/2时；出射点位置坐标(x1，z1+h)，进入的下一个单元网格信息为i+1，j+1，flag＝1；

当π/2<θ<3π/4时；出射点位置坐标(x1-h*tan(θ-π/2)，z1+h)进入的下一个单元网格信息为i，j+1，flag＝2；

当θ＝3π/4时；出射点位置坐标(x1-h，z1+h)进入的下一个单元网格信息为i-1，j+1，flag＝3；

当3π/4<θ<π时；出射点位置坐标(x1-h，z1+h*tan(π-θ))进入的下一个单元网格信息为i-1，j，flag＝4；

当π<＝θ<＝3π/2时；出射点位置坐标(x1-h，z1)进入的下一个单元网格信息为i-1，j-1，flag＝5；

当3π/2<θ<2π时；出射点位置坐标无；进入的下一个单元网格信息无。

所述步骤五还包括：当入射角方向flag＝4时，射线从网格右侧进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<θ<＝π/2时；出射点位置坐标(x1，j*h)；进入的下一个单元网格信息为i+1，j+1，flag＝1；

当π/2<θ<π-atan((j*h-z1)/h)时；出射点位置坐标(x1-(j*h-z1)tan(θ-π/2)，j*h)；进入的下一个单元网格信息为i，j+1，flag＝2；

当θ＝π-atan((j*h-z1)/h)时；出射点位置坐标(x1-h，j*h)；进入的下一个单元网格信息为i-1，j+1，flag＝3；

当π-atan((j*h-z1)/h)<θ<π+atan((z1-(j-1)*h)/h)时；出射点位置坐标(x1+h，z1+h*tan(π-θ))；进入的下一个单元网格信息为i-1，j，flag＝4；

当θ＝π+atan((z1-(j-1)*h)/h)时；出射点位置坐标(x1-h，(j-1)*h)；进入的下一个单元网格信息为i-1，j-1，flag＝5；

当π+atan((z1-(j-1)*h)/h)<θ<3π/2时；出射点位置坐标(x1-(z1-(j-1)*h)*tan(3π/2-θ)，(j-1)*h)进入的下一个单元网格信息为i，j-1，flag＝6；

当3π/2<＝θ<2π时；出射点位置坐标(x1，(j-1)*h)，进入的下一个单元网格信息为i+1，j-1，flag＝7。

所述步骤五还包括：当入射角方向flag＝5时，射线从网格右下侧端点进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<θ<π/2时；出射点位置坐标无，进入的下一个单元网格信息无；

当π/2<＝θ<＝π时；出射点位置坐标(x1-h，z1)，进入的下一个单元网格信息为i-1，j+1，flag＝3；

当π<θ<5π/4时；出射点位置坐标(x1-h，z1-h*tan(θ-π))，进入的下一个单元网格信息为i-1，j，flag＝4；

当θ＝5π/4时；出射点位置坐标(x1-h，z1-h)，进入的下一个单元网格信息为i-1，j-1，flag＝5；

当5π/4<θ<3π/2时；出射点位置坐标(x1-h*tan(3π/2-θ)，z1-h)，进入的下一个单元网格信息为i，j-1，flag＝6；

当3π/2<＝θ<＝2π时；出射点位置坐标(x1，z1-h)，进入的下一个单元网格信息为i+1，j-1，flag＝7。

所述步骤五还包括：当入射角方向flag＝6时，射线从网格下侧进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<＝θ<π/2时；出射点位置坐标(i*h，z1)，进入的下一个单元网格信息为i+1，j+1，flag＝1；

当π/2<θ<＝π时；出射点位置坐标((i-1)*h，z1)，进入的下一个单元网格信息为i-1，j+1，flag＝3；

当π<θ<3π/2-atan((x1-(i-1)*h)/h)时；出射点位置坐标((i-1)*h，z1-(x1-(i-1)*h)*tan(θ-π))，进入的下一个单元网格信息为i-1，j，flag＝4；

当θ＝3π/2-atan((x1-(i-1)*h)/h)时；出射点位置坐标((i-1)*h，z1-h)，进入的下一个单元网格信息为i-1，j-1，flag＝5；

当3π/2-atan((x1-(i-1)*h)/h)<θ<3π/2+atan((i*h-x1)/h)时；出射点位置坐标(x1-h*tan(3π/2-θ)，z1-h)，进入的下一个单元网格信息为i，j-1，flag＝6；

当θ＝3π/2+atan((i*h-x1)/h)时；出射点位置坐标(i*h，z1-h)，进入的下一个单元网格信息为i+1，j-1，flag＝7；

当3π/2+atan((i*h-x1)/h)<θ<2π时；出射点位置坐标(i*h，z1-(i*h-x1)*tan(2π-θ))进入的下一个单元网格信息为i+1，j，flag＝8。

所述步骤五包括：当入射角方向flag＝7时，射线从网格左下侧端点进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<＝θ<＝π/2时；出射点位置坐标(x1+h，z1)进入的下一个单元网格信息为i+1，j+1，flag＝1；

当π/2<θ<π时；出射点位置坐标无；进入的下一个单元网格信息为无；

当π<θ<＝3π/2时；出射点位置坐标(x1，z1-h)；进入的下一个单元网格信息为i-1，j-1，flag＝5；

当3π/2<θ<7π/4时；出射点位置坐标(x1+h*tan(θ-3π/2)，z1-h)；进入的下一个单元网格信息为i，j-1，flag＝6；

当θ＝7π/4时；出射点位置坐标(x1+h，z1-h)；进入的下一个单元网格信息为i+1，j-1，flag＝7；

当7π/4<θ<2π时；出射点位置坐标(x1+h，z1-h*tan(2π-θ))；进入的下一个单元网格信息为i+1，j，flag＝8；

当入射角方向flag＝8时，射线从网格左侧进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<θ<atan((j*h-z1)/h)时；出射点位置坐标(x1+h，z1+h*tan(θ))，进入的下一个单元网格信息为i+1，j，flag＝8；

当θ＝atan((j*h-z1)/h)时；出射点位置坐标(x1+h，j*h)，进入的下一个单元网格信息为i+1，j+1，flag＝1；

当atan((j*h-z1)/h)<θ<π/2时；出射点位置坐标(x1+(j*h-z1)*tan(π/2-θ)，j*h)，进入的下一个单元网格信息为i，j+1，flag＝2；

当π/2<＝θ<π时；出射点位置坐标(x1，j*h)，进入的下一个单元网格信息为i-1，j+1，flag＝3；

当π<θ<＝3π/2时；出射点位置坐标(x1，(j-1)*h)，进入的下一个单元网格信息为i-1，j-1，flag＝5；

当3π/2<θ<2π-atan((z1-(j-1)*h)/h)时；出射点位置坐标(x1+(z1-(j-1)*h)*tan(θ-3π/2)，(j+1)*h)，进入的下一个单元网格信息为i，j-1，flag＝6；

当θ＝2π-atan((z1-(j-1)*h)/h)时；出射点位置坐标(x1+h，(j-1)*h)，进入的下一个单元网格信息为i+1，j-1，flag＝7；

当2π-atan((z1-(j-1)*h)/h)<θ<＝2π时；出射点位置坐标(x1+h，z1-h*tan(2π-θ))，进入的下一个单元网格信息为i+1，j，flag＝8。

本发明的有益效果是：

本发明设计的一种复杂介质条件下射线追踪方法，在有限差分网格剖分的基础上，利用自适应程函方程计算网格节点上的初至波旅行时，进而求取每个单元网格内部的入射角信息，再结合检波点和震源点的相对位置，实现复杂介质条件下射线路径的精确求解。

附图说明

图1为本发明设计的一种复杂介质条件下射线追踪方法中有限差分离散网格剖分示意图；

图2为本发明设计的一种复杂介质条件下射线追踪方法中射线从网格左上侧端点进入单元网格图示；

图3为本发明设计的一种复杂介质条件下射线追踪方法中射线从网格上侧进入单元网格图示；

图4为本发明设计的一种复杂介质条件下射线追踪方法中为射线从网格右上侧端点进入单元网格图示；

图5为本发明设计的一种复杂介质条件下射线追踪方法中射线从网格右侧进入单元网格图示；

图6为本发明设计的一种复杂介质条件下射线追踪方法中射线从网格右下侧端点进入单元网格图示；

图7为本发明设计的一种复杂介质条件下射线追踪方法中射线从网格下侧进入单元网格图示；

图8为本发明设计的一种复杂介质条件下射线追踪方法中射线从网格左下侧端点进入单元网格图示；

图9为本发明设计的一种复杂介质条件下射线追踪方法中射线从网格左侧进入单元网格图示；

图10为本发明设计的一种复杂介质条件下射线追踪方法中具有纵向常梯度变化的慢度模型；

图11为图10模型下求取的初至波旅行时等时曲线和不同检波点接收的射线路径图示。

具体实施方式

下面结合附图与实施例对本发明进行进一步的介绍：一种复杂介质条件下射线追踪方法，包括如下步骤：

步骤一：对介质的慢度模型进行差分网格离散剖分，设定每个单元网格内部为常慢度，单元网格交点位置为待求取的初至波旅行时。如图1中空心圆所示；

步骤二：初始化震源点位置；利用自适应程函方程求解方法求取单元网格交点位置的初至波旅行时t(i,j)，其中，i和j分别代表了x和z方向离散化后的网格节点坐标；

步骤三：基于初至波旅行时，利用中心差分公式，求取每个单元网格内部的入射角θ(i,j)；

步骤四：初始化检波点位置信息和入射角方向；

步骤五：从检波点出发，根据射线入射点位置信息和单元网格内的入射角方向，计算当前单元网格内的射线出射点坐标和进入的下个单元网格信息；

步骤六：将计算的当前网格射线出射点坐标作为进入的下一个单元网格的入射点坐标，并利用入射方向(flag)计算下一个单元网格的出射点坐标和之后进入的单元网格信息；

步骤七：重复步骤五和六，直至射线出射点坐标到达震源位置。

所述步骤三中每个单元网格内部的入射角θ(i,j)求取公式为tz(i,j)＝((t(i+1,j+1)+t(i+1,j))-(t(i,j+1)+t(i,j)))/(2dz)

tx(i,j)＝((t(i+1,j+1)+t(i,j+1))-(t(i+1,j)+t(i,j)))/(2dx)

其中，dz和dx分别代表了z方向和x方向的离散采样间隔，tz(i,j)和tx(i,j)分别代表了初至波旅行时间场沿z方向和x方向的方向导数，atan函数为反正切函数。

所述步骤四初始化检波点位置信息和入射角方向还包括，初始化公式，其初始化公式为：

其中，rz和rx分别代表了检波点在z方向和x方向的坐标信息；bdy_up、bdy_bm、bdy_lt和bdy_rt分别代表了慢度模型上、下、左和右的边界坐标信息；flag代表了入射角方向的选择，如图1所示。

所述步骤五从检波点出发，根据射线入射点位置信息和单元网格内的入射角方向，计算当前单元网格内的射线出射点坐标和进入的下个单元网格信息还包括：当入射角方向flag＝1时，如图2所示，射线从网格左上侧端点进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<θ<π/4时，出射点位置坐标(x1+h，z1+h*tan(θ))，进入的下一个单元网格信息是i+1，j，flag＝8；

当θ＝π/4时，出射点位置坐标(x1+h，z1+h)；进入的下一个单元网格信息是i+1，j+1，flag＝1；

当π/4<θ<π/2时，出射点位置坐标(x1+h*tan(π/2-θ))；进入的下一个单元网格信息是i，j+1，flag＝2；

当π/2<＝θ<＝π时，出射点位置坐标(x1，z1+h)；进入的下一个单元网格信息是i-1，j+1，flag＝3；

当π<θ<3π/2时，出射点位置坐标无；进入的下一个单元网格信息无；

当3π/2<＝θ<＝2π时，出射点位置坐标(x1+h，z1)；进入的下一个单元网格信息为i+1，j-1，flag＝7。

其中z1和x1分别代表了入射点z方向和x方向的坐标，h代表了单元网格的剖分长度。

所述步骤五还包括：当入射角方向flag＝2时，如图3所示，射线从网格上侧进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<θ<atan(h/(i*h-x1))时，出射点位置坐标(i*h，z1+(i*h-x1)*tan(θ))；进入的下一个单元网格信息为i+1，j，flag＝8；

当θ＝atan(h/(i*h-x1))时，出射点位置坐标(i*h，z1+h)；进入的下一个单元网格信息为i+1，j+1，flag＝1；

当atan(h/(i*h-x1))<θ<π-atan(h/(x1-(i-1)*h))时；出射点位置坐标(x1+h*tan(π/2-θ)，z1+h)进入的下一个单元网格信息为i，j+1，flag＝2；

当θ＝π-atan(h/(x1-(i-1)*h))时；出射点位置坐标((i-1)*h，z1+h)进入的下一个单元网格信息为i-1，j+1，flag＝3；

当π-atan(h/(x1-(i-1)*h))<θ<π时；出射点位置坐标((i-1)*h，z1+(x1-(i-1)*h)tan(π-θ))进入的下一个单元网格信息为i-1，j，flag＝4；

当π<＝θ<3π/2时；出射点位置坐标((i-1)*h，z1)进入的下一个单元网格信息为i-1，j-1，flag＝5；

当3π/2<θ<＝2π时；出射点位置坐标(i*h，z1)进入的下一个单元网格信息为i+1，j-1，flag＝7。

所述步骤五还包括：当入射角方向flag＝3时，如图4所示，射线从网格右上侧端点进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<＝θ<＝π/2时；出射点位置坐标(x1，z1+h)，进入的下一个单元网格信息为i+1，j+1，flag＝1；

当π/2<θ<3π/4时；出射点位置坐标(x1-h*tan(θ-π/2)，z1+h)进入的下一个单元网格信息为i，j+1，flag＝2；

当θ＝3π/4时；出射点位置坐标(x1-h，z1+h)进入的下一个单元网格信息为i-1，j+1，flag＝3；

当3π/4<θ<π时；出射点位置坐标(x1-h，z1+h*tan(π-θ))进入的下一个单元网格信息为i-1，j，flag＝4；

当π<＝θ<＝3π/2时；出射点位置坐标(x1-h，z1)进入的下一个单元网格信息为i-1，j-1，flag＝5；

当3π/2<θ<2π时；出射点位置坐标无；进入的下一个单元网格信息无。

所述步骤五还包括：当入射角方向flag＝4时，如图5所示，射线从网格右侧进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<θ<＝π/2时；出射点位置坐标(x1，j*h)；进入的下一个单元网格信息为i+1，j+1，flag＝1；

当π/2<θ<π-atan((j*h-z1)/h)时；出射点位置坐标(x1-(j*h-z1)tan(θ-π/2)，j*h)；进入的下一个单元网格信息为i，j+1，flag＝2；

当θ＝π-atan((j*h-z1)/h)时；出射点位置坐标(x1-h，j*h)；进入的下一个单元网格信息为i-1，j+1，flag＝3；

当π-atan((j*h-z1)/h)<θ<π+atan((z1-(j-1)*h)/h)时；出射点位置坐标(x1+h，z1+h*tan(π-θ))；进入的下一个单元网格信息为i-1，j，flag＝4；

当θ＝π+atan((z1-(j-1)*h)/h)时；出射点位置坐标(x1-h，(j-1)*h)；进入的下一个单元网格信息为i-1，j-1，flag＝5；

当π+atan((z1-(j-1)*h)/h)<θ<3π/2时；出射点位置坐标(x1-(z1-(j-1)*h)*tan(3π/2-θ)，(j-1)*h)进入的下一个单元网格信息为i，j-1，flag＝6；

当3π/2<＝θ<2π时；出射点位置坐标(x1，(j-1)*h)，进入的下一个单元网格信息为i+1，j-1，flag＝7。

所述步骤五还包括：当入射角方向flag＝5时，如图6所示，射线从网格右下侧端点进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<θ<π/2时；出射点位置坐标无，进入的下一个单元网格信息无；

当π/2<＝θ<＝π时；出射点位置坐标(x1-h，z1)，进入的下一个单元网格信息为i-1，j+1，flag＝3；

当π<θ<5π/4时；出射点位置坐标(x1-h，z1-h*tan(θ-π))，进入的下一个单元网格信息为i-1，j，flag＝4；

当θ＝5π/4时；出射点位置坐标(x1-h，z1-h)，进入的下一个单元网格信息为i-1，j-1，flag＝5；

当5π/4<θ<3π/2时；出射点位置坐标(x1-h*tan(3π/2-θ)，z1-h)，进入的下一个单元网格信息为i，j-1，flag＝6；

当3π/2<＝θ<＝2π时；出射点位置坐标(x1，z1-h)，进入的下一个单元网格信息为i+1，j-1，flag＝7。

所述步骤五还包括：当入射角方向flag＝6时，如图7所示，射线从网格下侧进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<＝θ<π/2时；出射点位置坐标(i*h，z1)，进入的下一个单元网格信息为i+1，j+1，flag＝1；

当π/2<θ<＝π时；出射点位置坐标((i-1)*h，z1)，进入的下一个单元网格信息为i-1，j+1，flag＝3；

当π<θ<3π/2-atan((x1-(i-1)*h)/h)时；出射点位置坐标((i-1)*h，z1-(x1-(i-1)*h)*tan(θ-π))，进入的下一个单元网格信息为i-1，j，flag＝4；

当θ＝3π/2-atan((x1-(i-1)*h)/h)时；出射点位置坐标((i-1)*h，z1-h)，进入的下一个单元网格信息为i-1，j-1，flag＝5；

当3π/2-atan((x1-(i-1)*h)/h)<θ<3π/2+atan((i*h-x1)/h)时；出射点位置坐标(x1-h*tan(3π/2-θ)，z1-h)，进入的下一个单元网格信息为i，j-1，flag＝6；

当θ＝3π/2+atan((i*h-x1)/h)时；出射点位置坐标(i*h，z1-h)，进入的下一个单元网格信息为i+1，j-1，flag＝7；

当3π/2+atan((i*h-x1)/h)<θ<2π时；出射点位置坐标(i*h，z1-(i*h-x1)*tan(2π-θ))进入的下一个单元网格信息为i+1，j，flag＝8。

所述步骤五包括：当入射角方向flag＝7时，如图8所示，射线从网格左下侧端点进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<＝θ<＝π/2时；出射点位置坐标(x1+h，z1)进入的下一个单元网格信息为i+1，j+1，flag＝1；

当π/2<θ<π时；出射点位置坐标无；进入的下一个单元网格信息为无；

当π<θ<＝3π/2时；出射点位置坐标(x1，z1-h)；进入的下一个单元网格信息为i-1，j-1，flag＝5；

当3π/2<θ<7π/4时；出射点位置坐标(x1+h*tan(θ-3π/2)，z1-h)；进入的下一个单元网格信息为i，j-1，flag＝6；

当θ＝7π/4时；出射点位置坐标(x1+h，z1-h)；进入的下一个单元网格信息为i+1，j-1，flag＝7；

当7π/4<θ<2π时；出射点位置坐标(x1+h，z1-h*tan(2π-θ))；进入的下一个单元网格信息为i+1，j，flag＝8；

当入射角方向flag＝8时，如图9所示，射线从网格左侧进入单元网格内，计算单元网格内的射线出射点坐标和进入的下个单元网格信息有：

当0<θ<atan((j*h-z1)/h)时；出射点位置坐标(x1+h，z1+h*tan(θ))，进入的下一个单元网格信息为i+1，j，flag＝8；

当θ＝atan((j*h-z1)/h)时；出射点位置坐标(x1+h，j*h)，进入的下一个单元网格信息为i+1，j+1，flag＝1；

当atan((j*h-z1)/h)<θ<π/2时；出射点位置坐标(x1+(j*h-z1)*tan(π/2-θ)，j*h)，进入的下一个单元网格信息为i，j+1，flag＝2；

当π/2<＝θ<π时；出射点位置坐标(x1，j*h)，进入的下一个单元网格信息为i-1，j+1，flag＝3；

当π<θ<＝3π/2时；出射点位置坐标(x1，(j-1)*h)，进入的下一个单元网格信息为i-1，j-1，flag＝5；

当3π/2<θ<2π-atan((z1-(j-1)*h)/h)时；出射点位置坐标(x1+(z1-(j-1)*h)*tan(θ-3π/2)，(j+1)*h)，进入的下一个单元网格信息为i，j-1，flag＝6；

当θ＝2π-atan((z1-(j-1)*h)/h)时；出射点位置坐标(x1+h，(j-1)*h)，进入的下一个单元网格信息为i+1，j-1，flag＝7；

当2π-atan((z1-(j-1)*h)/h)<θ<＝2π时；出射点位置坐标(x1+h，z1-h*tan(2π-θ))，进入的下一个单元网格信息为i+1，j，flag＝8。

如图10所示在该模型的基础上，利用步骤一至七计算出基于初至波旅行时的射线路径，结果如图11所示，其中b线代表了初至波旅行时的等时曲线，a线代表了从震源出发，不同检波点接收的射线路径，可见该射线路径符合潜波传播的特征。

本发明提出的一种复杂介质条件下射线追踪方法并不限于以上所述的实施例，本领域的技术人员根据本发明的技术方案而得出的其他实施方式，满足本发明的原理方式，同样属于本发明的技术创新范畴。