CN106353793A - 一种基于走时增量双线性插值射线追踪的井间地震层析反演方法 - Google Patents

Abstract

本发明公开了一种基于走时增量双线性插值射线追踪的井间地震层析反演方法，包括获取井间地震资料的直达波和反射波走时；建立离散初始模型；计算直达波和反射波波前走时；计算直达波和反射波射线路径；建立层析反演方程；求解该反演方程，获取迭代后的速度模型。该方法充分考虑地层先验信息，网格剖分采用不规则六面体单元网格，假设从炮点沿直线路径和等效慢度传到单元边界上各点的走时与实际走时的差值双线性变化，把走时表示成背景走时与走时增量之和，消除了射线追踪方法在近源场存在的射线路径和走时计算误差较大的问题，进而形成提高了走时层析成像的精度和反演稳定性，可以提高对变速起伏地层和速度异常体的适应性。

Description

一种基于走时增量双线性插值射线追踪的井间地震层析反演方法

技术领域

本发明属于油气勘探地震资料处理领域，是一种井间地震资料层析速度反演的有效方法。

现有技术

目前，井间地震走时层析成像应用广泛，但分辨率和精度不高，主要是射线追踪方法的精度和适应性严重地影响着旅行时层析成像的精度和可靠性。常规的射线追踪方法主要包括试射法、弯曲法、最短路径法和波前射线追踪等方法。试射法在介质速度结构较复杂时，难以确定震源点到所有接收点的射线路径。弯曲法对层状介质模型的适应性较强，但有时会收敛到一个局部最小旅行时的路径，而且计算成本太高。最短路径法计算出的走时比实际走时系统偏大，射线路径呈“之”字形，也比实际路径长。基于程函方程的波前射线追踪方法，效率较高，而且能够容易追踪到全局最小旅行时的路径；但常常基于矩形单元模型，对层状介质模型的适应性较差。

发明内容

本发明的目的是针对常规射线追踪方法存在的射线路径和走时计算存在误差、计算精度不高的问题，提出了一种基于走时增量双线性插值射线追踪的井间地震层析反演方法。

本发明的一种基于走时增量双线性插值射线追踪的井间地震层析反演方法包括:

(1)获取井间地震资料的直达波和反射波走时；

(2)建立离散初始模型；

(3)计算直达波和反射波波前走时；

(4)计算直达波和反射波射线路径；

(5)建立层析反演方程；

(6)求解该反演方程；

(7)将第(2)步的速度模型作为初始速度模型进行迭代反演，通过第(6)步解反演方程，获取迭代后的速度模型，重复第(2)-(6)步，直到得到最终的速度模型。

上述方案进一步细化方案：

(1)获取地震资料的直达波和反射波走时；

(2)建立不规则单元离散化模型，采用地层产状先验信息确定模型地层产状基本架构，然后沿垂直z方向离散化模型，水平方向布置等间距的节点，采用不规则网格剖分法把初始模型离散成若干个不规则六面体单元，用网格节点处的速度值表示整个模型速度的变化；

(3)计算直达波或反射波波前走时，采用基于走时增量双线性插值的三维波前走时计算方法同时计算直达波和反射波波前走时；

(4)计算直达波或反射波射线路径，采用基于波前走时增量双线性插值的射线追踪方法计算射线路径；

(5)建立层析反演方程，充分考虑地层先验信息约束，联合直达波走时与反射波走时，并引入平滑因子稳定反演条件，建立直达波与反射走时残差与网格节点慢度和深度修正量之间满足的线性方程；

(6)求解该反演方程，采用阻尼LSQR算法求解该反演方程，带阻尼LSQR算法求解的是A^TAx＝A^Tb，并能够在求解的过程中施加一定的阻尼，当数据误差较大时能进行有效的抑制；

(7)将第(2)步的速度模型作为初始速度模型进行迭代反演，通过第(6)步解反演方程，获取迭代后的速度模型，重复第(2)-(6)步6-10次，直到满足收敛条件和反演效果。

上述方案进一步包括：

所述步骤(4)计算直达波或反射波射线路径中，对于直达波，在计算出的不规则离散介质模型中所有网格节点上的从激发点传播的波前传播时间后，根据互换原理，从接收点开始，反向确定满足Fermat原理的射线在相应单元界面上的位置，直至追踪到源点所在单元为止，从而获得相应的射线路径；对于反射波，先利用上行波波前时间，从接收点开始，反向确定满足Fermat原理的射线，直至进行到速度界面上为止，界面上的该点，即为对应的反射点；然后再利用下行波时间，从反射点开始，反向追踪反射点到源点的射线，最后连接源点至反射点以及反射点到接收点的射线，便得到从源点到接收点的反射射线路径；两者的整个过程均需要不断地在每个不规则六面体网格单元内确定射线路径。

所述步骤(4)计算直达波或反射波射线路径中，

单个不规则六面体网格单元内确定射线路径的方法为：

设接收点R在该单元内部或某个界面上，且坐标为(x_R,y_R,z_R)，根据最小时间原理，确定该单元内的射线问题就转化为在该单元界面上求一点S′，使波从S′到R的传播时间达到最小；假设所求的点S′位于六面单元的左界面上，坐标记为，左界面上四个顶点A,B,C,D的坐标分别为(x₁,y₁,z₁),(x₂,y₂,z₂),(x₃,y₃,z₃),(x₄,y₄,z₄)，且它们对应的波前时间分别为t₁,t₂,t₃,t₄，假设在该左界面上，波前时间线性变化，则S′点的波前时间可由这四个顶点上波前时间的双线性插值函数表示为：

t_S′＝a·Δx+b·Δy+c·Δx·Δy+d (1)

波经过S′传播到E的时间为：

$t_E=t_{S^{\′}}+s\·\sqrt{{(x_0+\mathrm{\Δ}x-x_E)}^2+{(y_0+\mathrm{\Δ}y-y_E)}^2+{(z_0+a_1\mathrm{\Δ}x+b_1\mathrm{\Δ}y+c_1-z_E)}^2}---\left(2\right)$

其中： $a_1=-\frac{A_0}{C_0},b_1=-\frac{B_0}{C_0},c_1=-\frac{A_0(x_0-x_1)+B_0(y_0-y_1)}{C_0}+z_1-z_0;$

而 $\left\{\begin{array}{c}{A}_0=(y_2-y_1)(z_3-z_1)-(y_3-y_1)(z_2-z_1)\\ B_0=(x_3-x_1)(z_2-z_1)-(x_2-x_1)(z_3-z_1)\\ C_0=(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)\end{array}\right.;$

(x₀,y₀,z₀)为界面的中心点坐标，即

x₀＝(x₁+x₂)/2，y₀＝(y₁+y₂)/2，z₀＝(z₁+z₃)/2

即可根据算式(2)计算出波经过六面体各个界面到达接收点R的传播时间和相应的S′点位置；

在这些传播时间中选取最小值，该最小值所对应的插值边界点S′便是射线路径的通过点，该点与接收点的连线即为该单元内的射线路径。

在此基础上，从源点到接收点的直达波射线追踪步骤如下：

步骤一：在接收点所在单元内，利用该单元网格节点上计算出的波前传播时间，按照上述方法确定到达接收点的射线与该单元界面的交点，如果接收点位于单元边界或网格节点上，这时要同时考虑接收点所在的所有单元的界面；

步骤二：以该交点为新的接收点，重复步骤一,以确定第二个单元网格上的射线与单元边界的交点，依次类推,确定射线与其穿过的所有单元界面的交点。

步骤三：若追踪到的交点位于源点所在单元内，则终止射线追踪；否则，返回步骤二继续追踪。

步骤四：把源点和确定的所有交点及接收点顺次相连，就得到了源点到该接收点的直达波射线路径；

从源点到接收点的反射波射线追踪分为反射路径和入射路径两部分,具体步骤如下：

(1)利用来自某反射面的反射波前时间，计算反射路径部分：

(1.1)在接收点所在单元内，利用该单元网格节点上计算出的反射波前传播时间，按照上述直达波射线路径确定方法确定到达接收点的射线与该单元界面的交点，如果接收点位于单元边界或网格节点上，这时要同时考虑接收点所在的所有单元的界面；

(1.2)若确定的交点位于反射面上，此交点即为该反射面上的反射点，即完成了反射路径部分的射线追踪，则可进行入射路径的射线追踪,其步骤按照下文中步骤(2),否则若仍未达到反射面，则按照步骤(1.3)继续进行反射路径的追踪；

(1.3)以该交点为新的接收点，重复步骤(1.1)，依次确定射线与其穿过的各个单元界面的交点,直到确定完反射路径的射线与其穿过的所有单元界面的交点。

(2)利用来自源点的入射波前时间，计算入射路径：

(2.1)以反射点为新的接收点，利用入射波前时间，重复步骤(1.1)；

(2.2)若确定的交点位于源点所在单元内，终止射线追踪；否则，跳到(2.3)；

(2.3)以该交点为新的接收点，重复步骤(2.1)，依次确定射线与其穿过的各个单元界面的交点；

(3)把源点和确定的所有交点及接收点顺次相连，就得到了源点到该接收点的反射射线路径。

所述步骤(5)建立层析反演方程中，层析反演的目标函数的表达式为：

$\begin{array}{c}\mathrm{\φ}(s,d)={\mathrm{\ρ}}_1{\left(t_{tran}^{cal}\right(s,d)-t_{tran}^{obs})}^T{C}_{tran}^{-1}\left(t_{tran}^{cal}\right(s,d)-t_{tran}^{obs})\\ +{\mathrm{\ρ}}_2{\left(t_{refl}^{cal}\right(s,d)-t_{refl}^{obs})}^T{C}_{refl}^{-1}\left(t_{refl}^{cal}\right(s,d)-t_{refl}^{obs})\\ +{\mathrm{\λ}}_{s1}\Σs^T{L}^T{C}_{slow}^{-1}Ls\\ +{\mathrm{\λ}}_{d1}\Σd^T{L}^T{C}_{depth}^{-1}Ld\\ +{\mathrm{\λ}}_{s2}{(s-s_{prior})}^T{C}_{slow}^{-1}(s-s_{prior})\\ +{\mathrm{\λ}}_{d2}{(d-d_{prior})}^T{C}_{depth}^{-1}(d-d_{prior})\end{array}---\left(3\right)$

其中，s表示模型慢度向量，d表示界面深度向量，C_tran、C_refl、C_slow、C_depth分别表示透射波旅行时、反射波旅行时、模型速度、界面深度空间对应的先验协方差矩阵；λs1、λd1、λs2、λd2分别为控制各项作用大小的权系数，ρ1和ρ2是取值为0和1的常系数，用于控制使用不同类型的旅行时数据；

上式第一项表示透射波旅行时实测值与计算值之间的差异；

第二项表示反射波旅行时实测值与计算值之间的差异；

第三项和号内表示一层介质慢度的平滑程度，求和后表示所有层慢度值的光滑程度，L为光滑算子，可以采用二阶差分算子，即

$s^T{L}^{T}Ls=\underset{k=1}{\overset{K}{\Σ}}\[{\left(\frac{{\∂}^{2}s}{\∂x^2}\right)}^2+{\left(\frac{{\∂}^{2}s}{\∂z^2}\right)}^2+2{\left(\frac{{\∂}^{2}s}{\∂x\∂z}\right)}^2\]$

用二阶差分计算；

第四项和号内表示一个界面的深度的平滑程度，求和后表示所有界面深度的光滑程度，L为光滑算子，可以采用二阶差分算子，即

$s^T{L}^{T}Ls=\underset{k=1}{\overset{K}{\Σ}}{\left(\frac{{\∂}^{2}d}{\∂x^2}\right)}^2$

用二阶差分计算；

第五项表示反演出的模型慢度值与已知的模型慢度值之间的差异；第六项表示反演出的界面深度值与已知的界面深度值之间的差异；这两项是用已知模型数据对反演结果进行约束；系数λs2或λd2是相应的权重因子。

发明效果

本发明的方法能够较好地获取两井之间的速度模型，有着其他技术不具备的优势，其具体优势和特点表现在以下几个方面：

第一、技术效果的可靠性。该方法充分考虑地层先验信息，网格剖分采用不规则六面体单元网格，使初始模型作为反演的约束条件，降低了反演的多解性，射线追踪时采用基于走时增量的双线性插值波前走时计算，把走时表示成背景走时与走时增量之和，假设走时增量在单元边界上双线性变化，消除了射线追踪方法在近源场存在的射线路径和走时计算误差较大的问题，相比单独射线方法更符合实际地震波传播路径，同时联合直达波与反射波两种有效波场走时信息，所得到的结果可靠性更强，效果明显。

第二、射线追踪效率和反演算法效率高，群进式波前扩展走时计算方法和走时增量双线性插值的射线追踪方法，直达波与反射波走时同时计算，计算效率高；阻尼LSQR全局最有解反演算法反演效率高，节省大量的内存占用，稳定性强。该方法流程及参数设置简单，运算速度快，适合应用于射线条数非常大的井间地震资料层析反演处理。

附图说明

图1一种基于走时增量双线性插值射线追踪的井间地震层析反演方法的处理流程图。

图2为不规则单元内射线穿过某一界面ABCD到达某一点R的几何关系图。图中A、B、C、D为单元左界面的顶点。S为射线与ABCD界面的交点，R为接收点。

图3为常规方法计算的走时误差与本发明方法计算的走时误差对比。图中左侧视图为常规方法计算的走时误差，右侧为本发明方法计算的走时误差。

图4为常规方法和本发明的基于走时增量双线性插值波前走时计算方法计算的直达波时间与理论值的比较。

图5为起伏地层模型的层析反演效果。图中：左侧为理论模型切片；中间为用常规射线追踪方法的层析反演切片；右侧用本发明方法的层析反演切片。

具体实施方式

下面结合附图1对本发明的技术方案做进一步说明。

一种基于走时增量双线性插值射线追踪的井间地震层析反演方法，包括：

(1)获取地震资料的直达波和反射波走时。

(2)建立不规则单元离散化模型。采用地层产状先验信息确定模型地层产状基本架构，然后沿垂直z方向离散化模型，水平方向布置等间距的节点，采用不规则网格剖分法把初始模型离散成若干个不规则六面体单元，用网格节点处的速度值表示整个模型速度的变化。

(3)计算直达波或反射波波前走时。采用基于走时增量双线性插值的三维波前走时计算方法同时计算直达波和反射波波前走时。

(4)计算直达波或反射波射线路径。采用基于波前走时增量双线性插值的射线追踪方法计算射线路径。对于直达波，在计算出的不规则离散介质模型中所有网格节点上的从激发点传播的波前传播时间后，根据互换原理，从接收点开始，反向确定满足Fermat原理的射线在相应单元界面上的位置，直至追踪到源点所在单元为止，从而获得相应的射线路径，提高射线追踪精度；对于反射波，先利用上行波波前时间，从接收点开始，反向确定满足Fermat原理的射线，直至进行到速度界面上为止，界面上的该点，即为对应的反射点；然后再利用下行波时间，从反射点开始，反向追踪反射点到源点的射线，最后连接源点至反射点以及反射点到接收点的射线，便得到从源点到接收点的反射射线路径。两者的整个过程均需要不断地在每个不规则六面体网格单元内确定射线路径。

(5)建立层析反演方程。充分考虑地层先验信息约束，联合直达波走时与反射波走时，并引入平滑因子稳定反演条件，建立直达波与反射走时残差与网格节点慢度和深度修正量之间满足的线性方程。

(6)求解该反演方程。采用阻尼LSQR算法求解该反演方程，带阻尼LSQR算法求解的是A^TAx＝A^Tb，并能够在求解的过程中施加一定的阻尼，当数据误差较大时能进行有效的抑制，提高了反演的抗噪声能力。与其他求解方法相比，可以节省较多的内存占用，同时提高计算效率。

(7)将第(2)步的速度模型作为初始速度模型进行迭代反演。通过第(6)步解反演方程，获取迭代后的速度模型，重复第(2)-(6)步6-10次，直到满足收敛条件和反演效果。

上述实施例的优化方案是：

(1)获取地震资料的直达波和反射波走时。

(2)建立不规则单元离散化模型。采用地层产状先验信息确定模型地层产状基本架构，然后沿垂直z方向离散化模型，水平方向布置等间距的节点，采用不规则网格剖分法把初始模型离散成若干个不规则六面体单元，用网格节点处的速度值表示整个模型速度的变化。

(3)计算直达波或反射波波前走时。采用基于走时增量双线性插值的三维波前走时计算方法同时计算直达波和反射波波前走时。先从源点开始，按照波前扩展思路采用效率非常高的GMM群快速步进方法计算波前走时至介质分界面，再从介质分界面上的最小走时节点开始向下，计算界面以下介质节点上的波前时间。

(4)计算直达波或反射波射线路径。采用基于波前走时增量双线性插值的射线追踪方法。对于直达波，在计算出的不规则离散介质模型中所有网格节点上的从激发点传播的波前传播时间后，根据互换原理，从接收点开始，反向确定满足Fermat原理的射线在相应单元界面上的位置，直至追踪到源点所在单元为止，从而获得相应的透射波射线路径，提高射线追踪精度；对于反射波，先利用上行波波前时间，从接收点开始，反向确定满足Fermat原理的射线，直至进行到速度界面上为止，界面上的该点，即为对应的反射点；然后，再利用下行波时间，从反射点开始，反向追踪反射点到源点的射线，最后连接源点至反射点以及反射点到接收点的射线，便得到从源点到接收点的反射射线路径。两者的整个过程均需要不断地在每个不规则六面体网格单元内确定射线路径。

其中单个不规则六面体网格单元内确定射线路径的方法为：

设接收点R在该单元内部或某个界面上，且坐标为(x_R,y_R,z_R)，如图2所示，根据最小时间原理，确定该单元内的射线问题就转化为在该单元界面上求一点S′，使波从S′到R的传播时间达到最小。假设所求的点S′位于图2的六面单元的左界面上，坐标记为左界面上四个顶点A,B,C,D的坐标分别为(x₁,y₁,z₁),(x₂,y₂,z₂),(x₃,y₃,z₃),(x₄,y₄,z₄)，且它们对应的波前时间分别为t₁,t₂,t₃,t₄。假设在该左界面上，波前时间线性变化，则S′点的波前时间可由这四个顶点上波前时间的双线性插值函数表示为：

t_S′＝a·Δx+b·Δy+c·Δx·Δy+d (1)

波经过S′传播到E的时间为：

$t_E=t_{S^{\′}}+s\·\sqrt{{(x_0+\mathrm{\Δ}x-x_E)}^2+{(y_0+\mathrm{\Δ}y-y_E)}^2+{(z_0+a_1\mathrm{\Δ}x+b_1\mathrm{\Δ}y+c_1-z_E)}^2}---\left(2\right)$

其中： $a_1=-\frac{A_0}{C_0},b_1=-\frac{B_0}{C_0},c_1=-\frac{A_0(x_0-x_1)+B_0(y_0-y_1)}{C_0}+z_1-z_0;$

而 $\left\{\begin{array}{c}{A}_0=(y_2-y_1)(z_3-z_1)-(y_3-y_1)(z_2-z_1)\\ B_0=(x_3-x_1)(z_2-z_1)-(x_2-x_1)(z_3-z_1)\\ C_0=(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)\end{array}\right.;$

(x₀,y₀,z₀)为界面的中心点坐标，即

x₀＝(x₁+x₂)/2，y₀＝(y₁+y₂)/2，z₀＝(z₁+z₃)/2

即可根据算式(2)计算出波经过六面体各个界面到达接收点R的传播时间和相应的S′点位置。

在这些传播时间中选取最小值，该最小值所对应的插值边界点S′便是射线路径的通过点，该点与接收点的连线即为该单元内的射线路径。在此基础上，从源点到接收点的直达波射线追踪步骤如下：

步骤一：在接收点所在单元内，利用该单元网格节点上计算出的波前传播时间，按照上述方法确定到达接收点的射线与该单元界面的交点，如果接收点位于单元边界或网格节点上，这时要同时考虑接收点所在的所有单元的界面；

步骤二：以该交点为新的接收点，重复步骤一,以确定第二个单元网格上的射线与单元边界的交点，依次类推,确定射线与其穿过的所有单元界面的交点。

步骤三：若追踪到的交点位于源点所在单元内，则终止射线追踪；否则，返回步骤二继续追踪。

步骤四：把源点和确定的所有交点及接收点顺次相连，就得到了源点到该接收点的直达波射线路径；

从源点到接收点的反射波射线追踪分为反射路径和入射路径两部分,具体步骤如下：

(1)利用来自某反射面的反射波前时间，计算反射路径部分：

(1.1)在接收点所在单元内，利用该单元网格节点上计算出的反射波前传播时间，按照上述直达波射线路径确定方法确定到达接收点的射线与该单元界面的交点，如果接收点位于单元边界或网格节点上，这时要同时考虑接收点所在的所有单元的界面；

(1.2)若确定的交点位于反射面上，此交点即为该反射面上的反射点，即完成了反射路径部分的射线追踪，则可进行入射路径的射线追踪,其步骤按照下文中步骤(2),否则若仍未达到反射面，则按照步骤(1.3)继续进行反射路径的追踪；

(1.3)以该交点为新的接收点，重复步骤(1.1)，依次确定射线与其穿过的各个单元界面的交点,直到确定完反射路径的射线与其穿过的所有单元界面的交点。

(2)利用来自源点的入射波前时间，计算入射路径：

(2.1)以反射点为新的接收点，利用入射波前时间，重复步骤(1.1)；

(2.2)若确定的交点位于源点所在单元内，终止射线追踪；否则，跳到(2.3)；

(2.3)以该交点为新的接收点，重复步骤(2.1)，依次确定射线与其穿过的各个单元界面的交点；

(3)把源点和确定的所有交点及接收点顺次相连，就得到了源点到该接收点的反射射线路径。

(5)建立层析反演方程，充分考虑地层先验信息约束，联合直达波走时与反射波走时，并引入平滑因子稳定反演条件，建立直达波与反射走时残差与网格节点慢度和深度修正量之间满足的线性方程；

层析反演的目标函数的表达式为：

$\begin{array}{c}\mathrm{\φ}(s,d)={\mathrm{\ρ}}_1{\left(t_{tran}^{cal}\right(s,d)-t_{tran}^{obs})}^T{C}_{tran}^{-1}\left(t_{tran}^{cal}\right(s,d)-t_{tran}^{obs})\\ +{\mathrm{\ρ}}_2{\left(t_{refl}^{cal}\right(s,d)-t_{refl}^{obs})}^T{C}_{refl}^{-1}\left(t_{refl}^{cal}\right(s,d)-t_{refl}^{obs})\\ +{\mathrm{\λ}}_{s1}\Σs^T{L}^T{C}_{slow}^{-1}Ls\\ +{\mathrm{\λ}}_{d1}\Σd^T{L}^T{C}_{depth}^{-1}Ld\\ +{\mathrm{\λ}}_{s2}{(s-s_{prior})}^T{C}_{slow}^{-1}(s-s_{prior})\\ +{\mathrm{\λ}}_{d2}{(d-d_{prior})}^T{C}_{depth}^{-1}(d-d_{prior})\end{array}---\left(3\right)$

其中，s表示模型慢度向量，d表示界面深度向量，C_tran、C_refl、C_slow、C_depth分别表示透射波旅行时、反射波旅行时、模型速度、界面深度空间对应的先验协方差矩阵；λs1、λd1、λs2、λd2分别为控制各项作用大小的权系数。ρ1和ρ2是取值为0和1的常系数，用于控制使用不同类型的旅行时数据。

上式第一项表示透射波(直达波)旅行时实测值与计算值之间的差异；

第二项表示反射波旅行时实测值与计算值之间的差异；

第三项和号内表示一层介质慢度的平滑程度，求和后表示所有层慢度值的光滑程度。L为光滑算子，可以采用二阶差分算子，即

$s^T{L}^{T}Ls=\underset{k=1}{\overset{K}{\Σ}}\[{\left(\frac{{\∂}^{2}s}{\∂x^2}\right)}^2+{\left(\frac{{\∂}^{2}s}{\∂z^2}\right)}^2+2{\left(\frac{{\∂}^{2}s}{\∂x\∂z}\right)}^2\]$

用二阶差分计算。

第四项和号内表示一个界面的深度的平滑程度，求和后表示所有界面深度的光滑程度。L为光滑算子，可以采用二阶差分算子，即

$s^T{L}^{T}Ls=\underset{k=1}{\overset{K}{\Σ}}{\left(\frac{{\∂}^{2}d}{\∂x^2}\right)}^2$

用二阶差分计算。

第五项表示反演出的模型慢度值与已知的模型慢度值之间的差异；第六项表示反演出的界面深度值与已知的界面深度值之间的差异。这两项是用已知模型数据对反演结果进行约束。系数λs2或λd2是相应的权重因子。

(6)求解该反演方程，采用阻尼LSQR算法求解该反演方程，带阻尼LSQR算法求解的是A^TAx＝A^Tb，并能够在求解的过程中施加一定的阻尼，当数据误差较大时能进行有效的抑制，提高了反演的抗噪声能力。与其他求解方法相比，可以节省较多的内存占用，同时提高计算效率；

(7)将第(2)步的速度模型作为初始速度模型进行迭代反演，通过第(6)步解反演方程，获取迭代后的速度模型，重复第(2)-(6)步6-10次，直到满足收敛条件和反演效果，收敛条件为使迭代均方根残差小于一个采样点，井间地震一般为0.5ms采样；反演效果要求为井旁层析反演的速度曲线与声波测井曲线速度误差小于10％，即得到最终的速度模型。

本发明的实施首先对理论变速模型的走时进行了波前走时计算：图3所示是用常规插值方法和本发明插值方法计算的波前走时与理论走时的比较。模型长度100米，厚度10米，深度方向500米，速度随深度线性变化，在地表速度2000m/s，模型底部速度3000m/s；炮点位于左边界25米深度上。离散单元大小10m×10m×10m。其中，左图是常规插值方法计算的波前走时与理论走时误差；右图是本发明插值方法计算的波前走时与理论走时误差；从结果可以看出，常规插值方法计算的走时误差在近场位置达到本发明插值0.05ms以上，本发明方法的走时误差基本在0.01ms以内，本发明的走时计算方法比常规方法计算的波前走时的误差减小了2倍以上，特别是在近场区，原方法的误差较大，而新方法误差小得多。

图4所示的是不同方法计算的直达波走时比较。图形中最上面的跳跃较多的曲线是常规方法计算的直达波时距曲线，其下为本发明方法计算的时距曲线，最下方曲线为理论走时曲线，结果表明：常规方法计算的直达波时距曲线有较大的误差，且偏移距越小，误差越大；而本发明方法计算的直达波走时误差较小，与理论值吻合得较好。

试验模型二为包含倾斜和下凹界面的起伏地层模型：井间距100m，井深960m，模型在垂直于两井连线的方向上的厚度为10m，由五层匀速介质组成，从上至下，速度分别为2000m/s，2300m/s，2000m/s和2200m/s；4个反射界面，从上到下分别是斜界面、下凹界面、斜界面和水平界面，该模型截面如图5左图所示。图5中间为使用常规射线追踪方法的层析反演结果，图5右侧为使用本发明方法的层析反演结果。可以看出，本发明方法结果远远好于原方法结果。每个分层速度与理论模型基本一致，并且界面反映准确，起伏地层的层间无假速度异常体存在，反演精度高。

Claims (5)

1.一种基于走时增量双线性插值射线追踪的井间地震层析反演方法，其特征在于包括:

（1）获取井间地震资料的直达波和反射波走时；

（2）建立离散初始模型；

（3）计算直达波和反射波波前走时；

（4）计算直达波和反射波射线路径；

（5）建立层析反演方程；

（6）求解该反演方程；

（7）将第（2）步的速度模型作为初始速度模型进行迭代反演，通过第（6）步解反演方程，获取迭代后的速度模型，重复第（2）-（6）步，直到得到最终的速度模型。

2.根据权利要求1所述的基于走时增量双线性插值射线追踪的井间地震层析反演方法，其特征在于包括：

（1）获取地震资料的直达波和反射波走时；

（2）建立不规则单元离散化模型，采用地层产状先验信息确定模型地层产状基本架构，然后沿垂直z方向离散化模型，水平方向布置等间距的节点，采用不规则网格剖分法把初始模型离散成若干个不规则六面体单元，用网格节点处的速度值表示整个模型速度的变化；

（3）计算直达波或反射波波前走时，采用基于走时增量双线性插值的三维波前走时计算方法同时计算直达波和反射波波前走时；

（4）计算直达波或反射波射线路径，采用基于波前走时增量双线性插值的射线追踪方法计算射线路径；

（5）建立层析反演方程，充分考虑地层先验信息约束，联合直达波走时与反射波走时，并引入平滑因子稳定反演条件，建立直达波与反射走时残差与网格节点慢度和深度修正量之间满足的线性方程；

（6）求解该反演方程，采用阻尼LSQR算法求解该反演方程，带阻尼LSQR算法求解的是A^TAx=A^Tb，并能够在求解的过程中施加一定的阻尼，当数据误差较大时能进行有效的抑制；

（7）将第（2）步的速度模型作为初始速度模型进行迭代反演，通过第（6）步解反演方程，获取迭代后的速度模型，重复第（2）-（6）步6-10次，直到满足收敛条件和反演效果。

3.根据权利要求2所述的基于走时增量双线性插值射线追踪的井间地震层析反演方法，其特征在于所述步骤（4）计算直达波或反射波射线路径中，对于直达波，在计算出的不规则离散介质模型中所有网格节点上的从激发点传播的波前传播时间后，根据互换原理，从接收点开始，反向确定满足Fermat原理的射线在相应单元界面上的位置，直至追踪到源点所在单元为止，从而获得相应的射线路径；对于反射波，先利用上行波波前时间，从接收点开始，反向确定满足Fermat原理的射线，直至进行到速度界面上为止，界面上的该点，即为对应的反射点；然后再利用下行波时间，从反射点开始，反向追踪反射点到源点的射线，最后连接源点至反射点以及反射点到接收点的射线，便得到从源点到接收点的反射射线路径；两者的整个过程均需要不断地在每个不规则六面体网格单元内确定射线路径。

4.根据权利要求3所述的基于走时增量双线性插值射线追踪的井间地震层析反演方法，其特征在于所述步骤（4）计算直达波或反射波射线路径中，

其中单个不规则六面体网格单元内确定射线路径的方法为：

设接收点R在该单元内部或某个界面上，且坐标为(x_R,y_R,z_R)，根据最小时间原理，确定该单元内的射线问题就转化为在该单元界面上求一点S′，使波从S′到R的传播时间达到最小；假设所求的点S′位于六面单元的左界面上，坐标记为左界面上四个顶点A,B,C,D的坐标分别为(x₁,y₁,z₁),(x₂,y₂,z₂),(x₃,y₃,z₃),(x₄,y₄,z₄)，且它们对应的波前时间分别为t₁,t₂,t₃,t₄，假设在该左界面上，波前时间线性变化，则S′点的波前时间可由这四个顶点上波前时间的双线性插值函数表示为：

t_S′=a·Δx+b·Δy+c·Δx·Δy+d (1)

波经过S′传播到E的时间为：

$t_E=t_{S^{\′}}+s\·\sqrt{{(x_0+\mathrm{\Δ}x-x_E)}^2+{(y_0+\mathrm{\Δ}y-y_E)}^2+{(z_0+a_1\mathrm{\Δ}x+b_1\mathrm{\Δ}y+c_1-z_E)}^2}---\left(2\right)$

其中： $a_1=-\frac{A_0}{C_0},b_1=-\frac{B_0}{C_0},c_1=-\frac{A_0(x_0-x_1)+B_0(y_0-y_1)}{C_0}+z_1-z_0;$

而 $\{\begin{array}{c}{A}_0=(y_2-y_1)(z_3-z_1)-(y_3-y_1)(z_2-z_1)\\ B_0=(x_3-x_1)(z_2-z_1)-(x_2-x_1)(z_3-z_1)\\ C_0=(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)\end{array};$

(x₀,y₀,z₀)为界面的中心点坐标，即

x₀=（x₁+x₂)/2，y₀=(y₁+y₂)/2，z₀=(z₁+z₃)/2

即可根据算式（2）计算出波经过六面体各个界面到达接收点R的传播时间和相应的S′点位置；

在这些传播时间中选取最小值，该最小值所对应的插值边界点S′便是射线路径的通过点，该点与接收点的连线即为该单元内的射线路径；在此基础上，从源点到接收点的直达波射线追踪步骤如下：

步骤一：在接收点所在单元内，利用该单元网格节点上计算出的波前传播时间，按照上述方法确定到达接收点的射线与该单元界面的交点，如果接收点位于单元边界或网格节点上，这时要同时考虑接收点所在的所有单元的界面；

步骤二：以该交点为新的接收点，重复步骤一,以确定第二个单元网格上的射线与单元边界的交点，依次类推,确定射线与其穿过的所有单元界面的交点。

步骤三：若追踪到的交点位于源点所在单元内，则终止射线追踪；否则，返回步骤二继续追踪。

步骤四：把源点和确定的所有交点及接收点顺次相连，就得到了源点到该接收点的直达波射线路径；

从源点到接收点的反射波射线追踪分为反射路径和入射路径两部分,具体步骤如下：

（1）利用来自某反射面的反射波前时间，计算反射路径部分：

(1.1)在接收点所在单元内，利用该单元网格节点上计算出的反射波前传播时间，按照上述直达波射线路径确定方法确定到达接收点的射线与该单元界面的交点，如果接收点位于单元边界或网格节点上，这时要同时考虑接收点所在的所有单元的界面；

(1.2)若确定的交点位于反射面上，此交点即为该反射面上的反射点，即完成了反射路径部分的射线追踪，则可进行入射路径的射线追踪,其步骤按照下文中步骤(2),否则若仍未达到反射面，则按照步骤(1.3)继续进行反射路径的追踪；

(1.3)以该交点为新的接收点，重复步骤(1.1)，依次确定射线与其穿过的各个单元界面的交点,直到确定完反射路径的射线与其穿过的所有单元界面的交点。

（2）利用来自源点的入射波前时间，计算入射路径：

(2.1)以反射点为新的接收点，利用入射波前时间，重复步骤(1.1)；

(2.2)若确定的交点位于源点所在单元内，终止射线追踪；否则，跳到(2.3)；

(2.3)以该交点为新的接收点，重复步骤(2.1)，依次确定射线与其穿过的各个单元界面的交点；

（3）把源点和确定的所有交点及接收点顺次相连，就得到了源点到该接收点的反射射线路径。

5.根据权利要求3或4所述的基于走时增量双线性插值射线追踪的井间地震层析反演方法，其特征在于所述步骤（5）建立层析反演方程中，层析反演的目标函数的表达式为：

$\begin{array}{c}\mathrm{\φ}(s,d)={\mathrm{\ρ}}_1{(t_{tran}^{cal}(s,d)-t_{tran}^{obs})}^T{C}_{tran}^{-1}(t_{tran}^{cal}(s,d)-t_{tran}^{obs})\\ +{\mathrm{\ρ}}_2{(t_{refl}^{cal}(s,d)-t_{refl}^{obs})}^T{C}_{refl}^{-1}(t_{refl}^{cal}(s,d)-t_{refl}^{obs})\\ +{\mathrm{\λ}}_{s1}\Σs^T{L}^T{C}_{slow}^{-1}Ls\\ +{\mathrm{\λ}}_{d1}\Σd^T{L}^T{C}_{depth}^{-1}Ld\\ +{\mathrm{\λ}}_{s2}{\left(s-s_{prior}\right)}^T{C}_{slow}^{-1}\left(s-s_{prior}\right)\\ +{\mathrm{\λ}}_{d2}{\left(d-d_{prior}\right)}^T{C}_{depth}^{-1}\left(d-d_{prior}\right)\end{array}---\left(3\right)$

其中，s表示模型慢度向量，d表示界面深度向量，C_tran、C_refl、C_slow、C_depth分别表示透射波旅行时、反射波旅行时、模型速度、界面深度空间对应的先验协方差矩阵;λs1、λd1、λs2、λd2分别为控制各项作用大小的权系数，ρ1和ρ2是取值为0和1的常系数，用于控制使用不同类型的旅行时数据；

上式第一项表示透射波旅行时实测值与计算值之间的差异；

第二项表示反射波旅行时实测值与计算值之间的差异；

第三项和号内表示一层介质慢度的平滑程度，求和后表示所有层慢度值的光滑程度，L为光滑算子，可以采用二阶差分算子，即

$s^T{L}^{T}Ls=\underset{k=1}{\overset{K}{\Σ}}\[{\left(\frac{{\∂}^{2}s}{\∂x^2}\right)}^2+{\left(\frac{{\∂}^{2}s}{\∂z^2}\right)}^2+2{\left(\frac{{\∂}^{2}s}{\∂x\∂z}\right)}^2\]$

用二阶差分计算；

第四项和号内表示一个界面的深度的平滑程度，求和后表示所有界面深度的光滑程度，L为光滑算子，可以采用二阶差分算子，即

$s^T{L}^{T}Ls=\underset{k=1}{\overset{K}{\Σ}}{\left(\frac{{\∂}^{2}d}{\∂x^2}\right)}^2$

用二阶差分计算；

第五项表示反演出的模型慢度值与已知的模型慢度值之间的差异；第六项表示反演出的界面深度值与已知的界面深度值之间的差异；这两项是用已知模型数据对反演结果进行约束；系数λs2或λd2是相应的权重因子。