三维地震射线追踪的双线性走时扰动插值方法与流程

本发明涉及一种三维地震射线追踪的双线性走时扰动插值方法，属于地震勘探技术领域。

背景技术：

地震射线追踪是研究地震波的一种快速有效技术，被广泛应用在地震层析成像、偏移成像、地震定位以及地震数据采集设计等领域。线性走时插值法包括二维走时插值法和三维走时插值法，是地震走时计算和和射线追踪的主要方法之一。该方法假设走时沿离散单元边界线性变化，单元边界上任意点的走时可通过相邻离散网格节点上的走时线性插值表示。

实际上，单元边界上的走时呈非线性变化，因而当离散单元较大时，线性走时插值法的线性走时假设会导致较大的走时和射线路径计算误差。为解决这一问题，提出了计算三维模型地震走时的基于规则六面体单元的线性走时扰动插值方法(zhangetal.，2015)。该方法将单元边界的实际走时分解为等效匀速介质中的参考走时和走时扰动，假设走时扰动呈线性变化，参考走时沿单元边界非线性变化，保留了实际走时沿单元边界非线性变化的特征，具有较高的走时计算精度。但该方法未涉及射线路径的追踪，并且规则六面体单元不利于模拟地形和地层界面的起伏以及各地层速度的变化。

由于在二维不规则四边形单元中，单元边界为直线，单元边界上任意一点的走时扰动可通过相邻两个网格节点处走时扰动的线性插值求得；类似地，在三维规则六面体单元中，单元边界为与某一坐标轴垂直的水平平面或垂直平面。但是，在不规则六面体单元中，单元边界上的四个节点可能不在同一平面上，即单元边界面可能为曲面，建立在单元边界为平面基础上的走时扰动插值方程不适用于这种情况。因此，需要推导适用于三维不规则单元的走时扰动插值方法。

技术实现要素：

针对现有技术存在的上述缺陷，本发明提出了一种三维地震射线追踪的双线性走时扰动插值方法，用于计算复杂三维模型的地震波前走时并追踪射线路径。

本发明所述的三维地震射线追踪的双线性走时扰动插值方法，包括如下步骤：

步骤s1：建立三维速度模型：根据各个地层分界面和速度变化确定各个地层中离散网格的大小，将三维模型离散为一系列不规则六面体单元，设置炮点与检波点位置；

步骤s2：从震源点开始，采用群波前扩展算法扩展波前：将震源点处波前走时值记为0，并且记k＝2，表示此处已获得波前走时；其余网格节点的波前走时值均记为无穷大，并令k＝0，表示此处网格节点未计算波前走时；

步骤s3：计算与炮点相邻的各个网格节点上的波传播时间：上述节点组成波前窄带l，令k＝1，表示该网格节点已计算出波传播时间，但还未确定波前时间；

步骤s4：在波前窄带l中按照一定规则找出次级震源点，各次级震源间的波传播能量不会相互干扰，将满足条件的网格节点从波前窄带l中移至次级震源集合b中，并令k＝2；

步骤s5：使用新提出的双线性走时扰动插值方法更新与次级震源集合b中的网格节点相邻的节点处的波传播时间，相邻节点需满足k≠2；若k＝0，则将该节点移至波前窄带l中并令k＝1；

步骤s6：若波前窄带集合l非空，则跳转至步骤s4；若集合l为空，则波前扩展结束；

步骤s7：在计算出模型中全部网格节点的波前走时后，根据费马原理和互换原理，从接收点出发，应用新提出的双线性走时扰动插值方法反向计算射线与检波点所在的离散单元各个边界的交点，并选取具有最小走时的交点作为射线与该离散单元边界的交点；

步骤s8：将该交点作为新的接收点：重复步骤s7计算射线与下一个离散单元边界的交点，若该离散单元为炮点所在单元则跳至步骤s9；

步骤s9：顺次连接炮点、射线与各单元边界的交点和接收点，便得到具有最小走时的直达波射线路径。

优选地，所述步骤s1中，模型离散过程具体按照以下步骤进行：

步骤s11：在x方向和y方向分别用一系列等间隔的垂向平面剖分模型；

步骤s12：由浅至深分别将每个地层划分成多个薄层：在同一地层中，薄层数目相等，并且在同一x、y坐标处，各薄层在z方向上的厚度相等；不同地层的薄层数目可以不同，这样就把三维速度模型离散成一系列不规则六面体单元。

优选地，所述步骤s4中，次级震源点集合b的选择需遵循如下规则：

b＝{(i,j,k)∈l,ti,j,k≤tl,min+δtl}(1)

其中，i，j，k分别表示离散网格节点在x，y，z三个方向的编号；l表示当前波前窄带；ti，j，k表示节点(i，j，k)处的波传播时间；tl，min表示波前窄带l中的最小波传播时间；δtl为波前窄带中波传播能量不会相互影响的两个相邻网格节点的波传播时间之差。

优选地，所述步骤s5中，新提出的双线性走时扰动插值方法计算步骤如下：

步骤s51：计算由炮点出发到达待计算点的射线与离散单元边界的交点：假设射线由炮点s(xs，ys，zs)出发，穿过离散单元左边界到达待计算点f(xf，yf，zf)；离散单元左边界顶点坐标为ei(i＝1，2，3，4)，且顶点处波前走时已知为ti(i＝1，2，3，4)；射线与离散单元左边界的交点e(xe，ye，ze)可表示为：

其中，

ax＝x2-x1，bx＝x3-x1，cx＝x1-x2-x3+x4，dx＝x1(3)

ay＝y2-y1，by＝y3-y1，cy＝y1-y2-y3+y4，dy＝y1(4)

az＝z2-z1，bz＝z3-z1，cz＝z1-z2-z3+z4，dz＝z1(5)

公式(2)中，m，n为点ei(i＝1，2，3，4)坐标对e点坐标的贡献权重，可由下式表示：

式中，

at＝t2-t1-s1(h2-h1)，bt＝t3-t1-s1(h3-h1)，ct＝t1-t2-t3+t4-s1(h1-h2-h3+h4)(13)

式中，hi(i＝1，2，3，4)为点s至点ei(i＝1，2，3，4)的直线距离；s1为炮点至单元左边界的参考慢度，有s1＝(t1/h1+t1/h1+t1/h1+t1/h1)/4；

步骤s52：利用与e点相邻的网格节点处的走时扰动，通过插值计算e点处的走时扰动，计算公式如下：

δte＝atm+btn+ctmn+dt(14)

式中，m，n可通过公式(6)-(7)计算得到，at，bt，ct可由公式(13)计算得到，dt＝t1-s1h1；

步骤s53：计算走时待计算点f处的走时：假设离散单元的平均慢度为s2，射线路径在离散单元中为直线，则点f处的走时可表示为：

本发明的有益效果是：本发明所述的三维地震射线追踪的双线性走时扰动插值方法，克服了基于三维规则单元的走时扰动插值方法难以适应复杂的地层起伏和速度变化的问题，实现了三维复杂介质中地震走时和射线追踪的高精度计算；本发明易于实现，对复杂模型具有很强的适应性，且计算精度和计算效率高。。

附图说明

图1为本发明的流程框图；

图2为本发明的模型离散及波前扩展示意图；

图3为本发明的离散单元走时计算示意图；

图4(a)为新提出的双线性走时扰动插值方法计算的波前走时相对误差分布图(y＝1000m切片)；

图4(b)为三维线性走时插值法计算的波前走时相对误差分布图(y＝1000m切片)；

图5(a)为新提出的双线性走时扰动插值方法计算的射线路径与理论路径对比图；

图5(b)为三维线性走时插值法计算的射线路径与理论路径对比图；

图6为检波点接收到的射线走时的相对误差图。

具体实施方式

为了使本发明目的、技术方案更加清楚明白，下面结合实施例，对本发明作进一步详细说明。

实施例1：

如图1至图6所示，本发明所述的三维地震射线追踪的双线性走时扰动插值方法，包括如下步骤：

步骤s1：建立三维速度模型，根据各个地层分界面和速度变化确定各个地层中离散网格的大小，将三维模型离散为一系列不规则六面体单元，设置炮点与检波点位置；

步骤s2：如图2所示，从震源点开始，采用群波前扩展算法(gmm)扩展波前。将震源点处波前走时值记为0，并且记k＝2，表示该处已获得波前走时。其余网格节点的波前走时值均记为无穷大，并令k＝0，表示该处网格节点未计算波前走时；

步骤s3：计算与炮点相邻的各个网格节点上的波传播时间，这些节点组成波前窄带l，令k＝1，表示该网格节点已计算出波传播时间，但还未确定波前时间；

步骤s4：在波前窄带l中按照一定规则找出次级震源点，各次级震源间的波传播能量不会相互干扰，将满足条件的网格节点从波前窄带l中移至次级震源集合b中，并令k＝2；

步骤s5：使用新提出的双线性走时扰动插值方法更新与次级震源集合b中的网格节点相邻的节点处的波传播时间，相邻节点需满足k≠2。若k＝0，则将该节点移至波前窄带l中并令k＝1；

步骤s6：若波前窄带集合l非空，则跳转至步骤s4；若集合l为空，则波前扩展结束；

步骤s7：在计算出模型中全部网格节点的波前走时后，根据费马原理和互换原理，从接收点出发，应用新提出的双线性走时扰动插值方法反向计算射线与检波点所在的离散单元各个边界的交点，并选取具有最小走时的交点作为射线与该离散单元边界的交点；

步骤s8：将该交点作为新的接收点。重复步骤s7计算射线与下一个离散单元边界的交点，若该离散单元为炮点所在单元则跳至步骤s9；

步骤s9：顺次连接炮点、射线与各单元边界的交点和接收点，便得到具有最小走时的直达波射线路径。

进一步的，所述步骤s1中，如图2所示的模型离散过程具体按照以下步骤进行：

步骤s11：在x方向和y方向分别用一系列等间隔的垂向平面剖分模型；

步骤s12：由浅至深分别将每个地层划分成多个薄层。在同一地层中，薄层数目相等，并且在同一x、y坐标处，各薄层在z方向上的厚度相等。不同地层的薄层数目可以不同。这样,就把三维速度模型离散成一系列不规则六面体单元。

进一步的，所述步骤s4中，次级震源点集合b的选择需遵循如下规则：

b＝{(i,j,k)∈l,ti,j,k≤tl,min+δtl}(1)

其中，i,j,k分别表示离散网格节点在x,y,z三个方向的编号；l表示当前波前窄带；ti,j,k表示节点(i,j,k)处的波传播时间；tl，min表示波前窄带l中的最小波传播时间；δtl为波前窄带中波传播能量不会相互影响的两个相邻网格节点的波传播时间之差。

进一步的，所述步骤s5中，新提出的双线性走时扰动插值方法计算步骤如下：

步骤s51：计算由炮点出发到达待计算点的射线与离散单元边界的交点。为方便理解，假设射线由炮点s(xs,ys,zs)出发，穿过如图3所示的离散单元左边界到达待计算点f(xf,yf,zf)。离散单元左边界顶点坐标为ei(i＝1,2,3,4)，且顶点处波前走时已知为ti(i＝1,2,3,4)。射线与离散单元左边界的交点e(xe,ye,ze)可表示为：

其中，

ax＝x2-x1，bx＝x3-x1，cx＝x1-x2-x3+x4，dx＝x1(3)

ay＝y2-y1，by＝y3-y1，cy＝y1-y2-y3+y4，dy＝y1(4)

az＝z2-z1，bz＝z3-z1，cz＝z1-z2-z3+z4，dz＝z1(5)

公式(2)中，m,n为点ei(i＝1,2,3,4)坐标对e点坐标的贡献权重，可由下式表示：

式中，

at＝t2-t1-s1(h2-h1)，bt＝t3-t1-s1(h3-h1)，ct＝t1-t2-t3+t4-s1(h1-h2-h3+h4)(13)

式中，hi(i＝1,2,3,4)为点s至点ei(i＝1,2,3,4)的直线距离；s1为炮点至单元左边界的参考慢度，有s1＝(t1/h1+t1/h1+t1/h1+t1/h1)/4。

步骤s52：利用与e点相邻的网格节点处的走时扰动，通过插值计算e点处的走时扰动，计算公式如下：

δte＝atm+btn+ctmn+dt(14)

式中，m,n可通过公式(4)计算得到，at,bt,ct可由公式(13)计算得到，dt＝t1-s1h1。

步骤s53：计算走时待计算点f处的走时。假设离散单元的平均慢度为s2,射线路径在离散单元中为直线，则点f处的走时可表示为

实施例2：

为了进一步说明本方法的实现思路及实现过程并证明方法的准确性，使用三维模型进行测试，并和三维线性走时插值方法的计算结果进行比较。因匀速模型的理论射线为直线，因此选用匀速模型，便于说明射线追踪和走时计算误差。

s1：建立速度为2000m/s的匀速模型，模型尺寸为3500×3500×3500m³，离散网格大小为50×50×50m³，炮点位于(1000m,1000m,3200m)，检波点共9个，以400m间距均匀分布在检波点1(3000m,200m,0m)与检波点9(3000m,3400m,0m)连接的直线上；

s2：应用群波前扩展算法从炮点开始扩展波前，利用新提出的双线性走时扰动插值方法计算模型中各个节点处的波前走时，y＝1000m切片处的波前走时相对误差分布如图4(a)所示,最大相对误差为1.0×10^-11。图4(b)为三维线性走时插值方法的计算误差，最大相对误差达1.0×10^-1。可以看出使用新提出的双线性走时扰动插值方法计算走时的相对误差远小于线性走时插值方法的走时相对误差。

s3：从接收点出发使用新提出的双线性走时扰动插值方法反向逐步计算射线与相应单元边界的交点，直至追踪至炮点为止。顺次连接炮点、射线与各单元边界的交点、接收点，所得到的直达波射线路径与理论射线路径的比较如图5(a)所示。图5(b)为使用三维线性走时插值方法的射线路径追踪与理论射线路径的比较，图6所示是各个检波点接收到的射线走时与理论走时的相对误差。可以看出，使用新提出的双线性走时扰动插值方法计算的射线路径和射线走时其精度均远高于三维线性走时插值方法。

本发明可广泛运用于地震勘探场合。

以上所述仅为本发明的较佳实施例而己，并不以本发明为限制，凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等，均应包含在本发明的专利涵盖范围内。