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

技术特征：

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是相应的权重因子。