CN105425286A - 地震走时获取方法及基于其的井间地震走时层析成像方法 - Google Patents

Abstract

本申请实施例提供了一种地震走时获取方法及基于其的井间地震走时层析成像方法，该成像方法包括：首先利用初至提取算法的方式获得地震剖面的初至，然后建立一个均匀初始慢度模型，接着使用高阶多方向程函方程有限差分法计算慢度模型的地震走时，再用最速下降法追踪射线路径，进而构建出雅可比矩阵，最后进行走时层析成像来反演慢度场，得到一个新的慢度模型。通过迭代线性反演算法不断更新慢度模型并计算误差，最终得到一个误差足够小的慢度模型来对井间慢度场进行描述，从而获得比基于常规FMM算法的井间地震走时层析成像更好的重建结果和更高的反演精度。

Description

地震走时获取方法及基于其的井间地震走时层析成像方法

技术领域

本申请涉地震层析成像技术领域，尤其是涉及一种地震走时获取方法及基于其的井间地震走时层析成像方法。

背景技术

井间地震走时层析成像是石油开发、确定剩余油分布、探测未钻遇砂岩透镜体位置等方法的一种新的有效手段。其是由井中激发的地震波在井中接收，当地震波穿过井间地层中的异常结构体时其传播时间、振幅及波形均会发生较之正常的偏离变化，因此，可以利用这种偏离变化通过一些反演公式获得地层内部异常结构体的物性参数分布，从而揭示其地质构造、岩性分布和矿藏形态。

目前井间地震走时层析成像作为一种迭代反演公式，其层析反演过程需要反复实施走时正演计算，而正演算法的精度和效率则决定了整个反演流程的精度和效率。由此可见，地震走时计算是井间地震走时层析成像关键的步骤之一。然而，由于现有的地震走时计算方法精度不高，这导致井间地震走时层析成像结果也受到了影响。

发明内容

本申请实施例的目的在于提供一种地震走时获取方法及基于其的井间地震走时层析成像方法，以提高地震走时计算的精度，从而提高井间地震走时层析成像的质量。

为达到上述目的，一方面，本申请实施例提供了一种地震走时获取方法，包括以下步骤：

拾取地震记录的初至，得到其地震走时观测值；

依据所述地震走时观测值建立均匀网格化慢度模型；

确定所述均匀网格化慢度模型中源点邻点的地震走时；

基于高阶多方向程函方程有限差分法确定所述均匀网格化慢度模型中除所述源点邻点外其他网格节点的地震走时。

本申请实施例的地震走时获取方法，所述依据所述地震走时观测值建立均匀网格化慢度模型具体包括：

根据公式获取预设初始慢度模型的速度值v₀，其中，M为地震射线条数，l_i为第i条地震射线的长度，t_i为拾取的第i条地震射线的地震走时观测值；

取所述预设初始慢度模型的速度值v₀的倒数，获得所述预设初始慢度模型的慢度值s₀；

使用正方形网格将所述预设初始慢度模型均匀划分为多个网格节点，并将每个网格节点的慢度值赋值为s₀，得到均匀网格化慢度模型。

本申请实施例的地震走时获取方法，所述确定所述均匀网格化慢度模型中源点邻点的地震走时，具体包括：

在所述均匀网格化慢度模型中：

将震源所在位置作为源点，将源点标记为近点，并将其地震走时赋值为0；

将除所述近点之外的所有网格节点标记为远点并将其赋值为预设的地震走时初值；

根据一阶程函方程有限差分公式 $max{\left(\frac{1}{h}\right({T^{\′}}_{i,j}-{T^{\′}}_1),0)}^2+max(\frac{1}{h}({T^{\′}}_{i,j}-$ ${T^{\′}}_2),0{)}^2=S_{i,j}^{\′2}$ 获取所述近点的邻点的地震走时；

式中，T′_i,j为邻点处的地震走时，下角标i和j为该邻点在慢度模型内的索引位置，S′_i,j为邻点处的慢度值，T′₁＝min(T′_i-1,j,T′_i+1,j)为索引位置(i-1,j)和(i+2,j)处地震走时的最小值，T′₂＝min(T′_i,j-1,T′_i,j+1)为索引位置(i,j-1)和(i,j+2)处地震走时的最小值，h为网格间距。

本申请实施例的地震走时获取方法，所述基于高阶多方向程函方程有限差分法确定所述均匀网格化慢度模型中除所述源点邻点外其他网格节点的地震走时，具体包括：

将当前近点的所有邻点放入预设的波前窄带中并标记为窄带点；

在所述波前窄带中选取其中地震走时最小的窄带点并将其移出所述波前窄带作为新的近点；

根据公式 $max{\left(\frac{1}{\mathrm{\Δ}d}\right(T_{i,j}-T_1),0)}^2+max{\left(\frac{1}{\mathrm{\Δ}d}\right(T_{i,j}-T_2),0)}^2=S_{i,j}^2$ 获取或更新与所述新的近点相邻的窄带点或远点的地震走时，式中，T_i,j为与所述新的近点相邻的窄带点或远点的地震走时，S_i,j为其对应的慢度值，T₁和T₂为近点的地震走时，Δd为节点间距；

重复以上步骤，直至所述波前窄带中为空时停止，从而计算出所述均匀网格化慢度模型中除所述源点邻点外其他网格节点的地震走时。

另一方面，本申请实施例还提供了一种基于上述地震走时获取方法的井间地震走时层析成像方法，包括以下步骤：

基于高阶多方向程函方程有限差分法确定均匀网格化慢度模型中所有网格节点的地震走时；

利用最速下降法由各检波点向源点反向追踪对应的地震射线路径；

计算所述每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段的长度，并构建所述均匀网格化慢度模型的雅可比矩阵；

基于所述每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段的长度，计算每条地震射线路径的传播时间，得到所述均匀网格化慢度模型的走时向量；

根据所述雅可比矩阵、所述走时向量以及预设的迭代线性反演公式及算法对所述均匀网格化慢度模型进行迭代反演，以更新所述均匀网格化慢度模型；

重复以上步骤，直至当前更新后的均匀网格化慢度模型与所述地震走时观测值之间的关系满足设定条件，此时所述当前更新后的均匀网格化慢度模型为所求慢度模型。

本申请实施例的井间地震走时层析成像方法，均匀网格化慢度模型与地震走时观测值之间的误差小于设定阈值。

本申请实施例的井间地震走时层析成像方法，所述利用最速下降法由各检波点向源点反向追踪对应的地震射线路径中，每条地震射线路径在所述均匀网格化慢度模型中各个网格节点的最速下降方向角通过以下公式计算得到：

$\mathrm{\θ}={\tan }^{-1}\[(\∂T/\∂y)/(\∂T/\∂x)\]+\mathrm{\π};$

式中， $\frac{\∂T}{\∂x}=\frac{\[(T_{i+1,j}+T_{i+1,j+1})-(T_{i,j}+T_{i,j+1})\]}{2h}$ 和 $\frac{\∂T}{\∂y}=\frac{\[(T_{i,j+1}+T_{i+1,j+1})-(T_{i,j}+T_{i+1,j})\]}{2h}$ 分别为走时梯度向量的水平分量和竖直分量的中心差分格式，T_i,j为所述均匀网格化慢度模型中中索引位置(i,j)处的地震走时。

本申请实施例的井间地震走时层析成像方法，基于所述每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段的长度，计算每条地震射线路径的传播时间，具体包括：

根据公式d(s)＝∑_ms_ml_m计算每条地震射线路径的传播时间；

式中，m是每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段标记，l_m和s_m分别为第m条射线段的长度及该射线段所在网格节点内波传播的平均慢度。

本申请实施例的井间地震走时层析成像方法，所述预设的迭代线性反演公式包括：

$\left[\begin{array}{c}{W}_d{L}^k\\ {\mathrm{\λW}}_m\end{array}\right](s^{k+1}-s^k)=\left[\begin{array}{c}{W}_d\[t-d\left(s^k\right)\]\\ 0\end{array}\right],k=0,1,2,3,...$

式中，W_d和W_m是两个对称矩阵，λ为加权因子，L^k为第k次迭代的雅可比矩阵，d(s^k)为第k次迭代的均匀网格化慢度模型的走时向量，s^k为第k次迭代慢度向量。

本申请实施例的井间地震走时层析成像方法，所述预设的迭代线性反演算法包括最小二乘反演算法。

附图说明

此处所说明的附图用来提供对本申请实施例的进一步理解，构成本申请实施例的一部分，并不构成对本申请实施例的限定。在附图中：

图1为本申请一实施例的地震走时获取方法的流程图；

图2为本申请实施例的地震走时获取方法中，共炮集井间地震记录初至拾取示意图；

图3为本申请实施例的地震走时获取方法中高阶多方向程函方程有限差分算法中所定义的三种类型的网格节点示意图；

图4a-4d为本申请实施例的地震走时获取方法中使用一阶差分计算源点的相邻的走时示意图；其中，图4a为左侧邻点、图4b为右侧邻点、图4c为上侧邻点、图4d为下侧邻点；

图5a-5d为本申请实施例的地震走时获取方法中高阶多方向程函方程有限差分算法中所使用的四个方向；其中，图5a为方向1、图5b为方向2、图5c为方向3、图5d为方向4；

图6a-6c为本申请实施例的地震走时获取方法中均匀网格化慢度模型使用高阶多方向程函方程有限差分算法计算的走时与理论走时之间的绝对误差分布；其中，图6a为一阶差分绝对误差、图6b为二阶差分绝对误差、图6c为三阶差分绝对误差；

图7为本申请一实施例的井间地震走时层析成像方法的流程图；

图8为本申请实施例的井间地震走时层析成像方法中最速下降法中射线追踪示意图；

图9a-9d为本申请实施例的井间地震走时层析成像方法的走时层析成像反演结果与基于现有技术的走时层析成像反演结果的对比示意图；其中，图9a为真实模型、图9b为基于现有快速推进法(FMM)的走时层析成像演结果、图9c为二阶多方向有限差分走时层析成像反演结果、图9d为三阶多方向有限差分走时层析成像反演结果。

具体实施方式

为使本申请实施例的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本申请实施例做进一步详细说明。在此，本申请实施例的示意性实施例及其说明用于解释本申请实施例，但并不作为对本申请实施例的限定。

下面结合附图，对本申请实施例的具体实施方式作进一步的详细说明。

参考图1所示，本申请实施例的地震走时获取方法包括以下步骤：

S11，拾取地震记录的初至，得到其地震走时观测值。井间地震走时层析成像中所使用的观测数据是地震走时观测值，也就是检波点处接收到的直达波初至。而这些数据则需要对直达波的波形进行初至提取才能获得。因此，初至提取的精度将直接影响到走时层析成像的反演效果。对信噪比较高的数据，可使用初至提取算法提取初至，如传统的信噪比法和能量比法等，或最新的互相关法和数字图像分割法等；对于信噪比较低的数据，也可采用人工手动拾取的方式提取初至。如图2所示，为共炮集井间地震记录初至拾取示意图，图2中左侧所示的黑线即为所拾取的初至位置。

S12，依据所述地震走时观测值建立均匀网格化慢度模型。其中，依据地震走时观测值建立均匀网格化慢度模型具体包括：

根据公式获取预设初始慢度模型的速度值v₀，其中，M为地震射线条数，l_i为第i条地震射线的长度，t_i为拾取的第i条地震射线的地震走时观测值；

取所述预设初始慢度模型的速度值v₀的倒数，获得所述预设初始慢度模型的慢度值s₀；

使用正方形网格将所述预设初始慢度模型均匀划分为多个网格节点(即使用正方形网格将预设初始慢度模型划分成m行n列的离散模型)，并将每个网格节点的慢度值赋值为s₀，得到均匀网格化慢度模型。

S13，确定所述均匀网格化慢度模型中源点邻点的地震走时。其中，确定所述均匀网格化慢度模型中源点邻点的地震走时，具体包括：

在所述均匀网格化慢度模型中：

将震源所在位置作为源点，将源点标记为近点，并将其地震走时赋值为0；对于源点的邻点(即相邻网格节点)而言，其相邻的网格节点中只有源点的地震走时是已知的，因此源点的邻点的波前传播时间只能用一阶差分的形式通过源点的地震走时来计算(其中，源点的邻点如图4a-4d所示)。

将除所述近点之外的所有网格节点标记为远点并将其赋值为预设的地震走时初值；其中，该预设的地震走时初值一般赋予一个足够大的波前传播时间初值(如10⁵)。

根据一阶程函方程有限差分公式 $max{\left(\frac{1}{h}\right({T^{\′}}_{i,j}-{T^{\′}}_1),0)}^2+max(\frac{1}{h}({T^{\′}}_{i,j}-$ ${T^{\′}}_2),0{)}^2=S_{i,j}^{\′2}$ 获取所述近点的邻点的地震走时；

式中，T′_i,j为邻点处的地震走时，下角标i和j为该邻点在慢度模型内的索引位置，S′_i,j为邻点处的慢度值，T′₁＝min(T′_i-1,j,T′_i+1,j)为索引位置(i-1,j)和(i+2,j)处地震走时的最小值，T′₂＝min(T′_i,j-1,T′_i,j+1)为索引位置(i,j-1)和(i,j+2)处地震走时的最小值，h为网格间距。S14，基于高阶多方向程函方程有限差分法确定所述均匀网格化慢度模型中除所述源点邻点外其他网格节点的地震走时。其中，基于高阶多方向程函方程有限差分法确定所述均匀网格化慢度模型中除所述源点邻点外其他网格节点的地震走时，具体包括：

将当前近点的所有邻点放入预设的波前窄带中并标记为窄带点；

在所述波前窄带中选取其中地震走时最小的窄带点并将其移出所述波前窄带作为新的近点；

根据公式 $max{\left(\frac{1}{\mathrm{\Δ}d}\right(T_{i,j}-T_1),0)}^2+max{\left(\frac{1}{\mathrm{\Δ}d}\right(T_{i,j}-T_2),0)}^2=S_{i,j}^2$ 获取或更新与所述新的近点相邻的窄带点或远点的地震走时；当然，这里一般要求更新后的值比之前的值小，式中，T_i,j为与所述新的近点相邻的窄带点或远点的地震走时，S_i,j为其对应的慢度值，T₁和T₂为近点的地震走时，Δd为节点间距；

重复以上步骤，直至所述波前窄带中为空时停止，从而计算出所述均匀网格化慢度模型中除所述源点邻点外其他网格节点的地震走时。

在高阶多方向程函方程有限差分算法中，对于不同的方向以及所采用的不同阶数的有限差分，上式中的T₁、T₂以及Δd的值是不同的。如图5a-5d所示，为高阶多方向程函方程有限差分算法中所使用的四个方向。其中方向1为水平方向，方向2为垂直方向，方向3为主对角线方向，方向4为副对角线方向。计算中需要在这四个方向上分别对走时程函方程中的偏导进行有限差分离散，并且有限差分可以选择不同的阶数：阶数越高，差分中所使用的网格节点数越多，精度也就越高。下面所示的表1至表4给出了不同方向、不同阶数的有限差分中T₁、T₂以及Δd值的选择方法，具体的：

表1

表2

表3

表4

将方向1和方向2中分别对应的T₁、T₂以及Δd值代入上面的方程中并求解，可以得到一个解；将方向3和方向4中分别对应的T₁、T₂以及Δd值代入上面的方程中并求解，又可以得到一个解。由于T(i,j)的值是未知的，需要检查计算出的两个解T(i,j)是否大于邻近的T₁和T₂。通常只有大于T₁和T₂的解才是被接受的。

循环执行步骤S14，直到波前窄带中为空时停止。这样就计算出了均匀网格化慢度模型中所有网格节点的波前传播时间(即计算出了均匀网格化慢度模型中所有网格节点的地震走时)。

图6a-6c为均匀网格化慢度模型使用高阶多方向程函方程有限差分算法计算的走时与理论走时之间的绝对误差分布，其中图6a使用的是一阶差分，图6b使用的是二阶差分，图6c使用的是三阶差分。由此可以看出，差分阶数越高，绝对误差相对来说越小。

本申请实例中，上述S13和S14中出现了近点、远点和窄带点的概念，下面就这些概念进行解释：

在高阶多方向程函方程有限差分走时算法中，认为波是以波前扩散的方式向外传播的，而不是以传统假设的射线传播方式。因此，均匀网格化慢度模型中网格节点上的地震走时是波前传播时间。

高阶多方向程函方程有限差分算法中使用三种类型的点，如图3所示，其中第一种为近点(标记为M＝2)，表示已经计算过波前传播时间，并且波前传播时间不再被更新的节点；第二种为窄带点(标记为M＝1)，表示计算过波前传播时间，但是波前传播时间可能会被更新的节点；第三种为远点(标记为M＝0)，表示尚未计算波前传播时间的节点。

本申请实施例的地震走时获取方法依据实际观测的地震走时观测值建立均匀网格化慢度模型，并采用高阶多方向程函方程有限差分算法求解该均匀网格化慢度模型中网格节点的地震走时，由于差分阶数越高，绝对误差相对来说越小，因此，本申请实施例的地震走时获取方法可根据误差要求来设定差分阶数，因而保证了所获取的地震走时精度，从而为后续走时层析成像反演结果的奠定较好的基础。

结合图7所示，本申请实施例的基上述地震走时获取方法的井间地震走时层析成像方法包括以下步骤：

S71，基于高阶多方向程函方程有限差分法确定均匀网格化慢度模型中所有网格节点的地震走时。具体请参见上述地震走时获取方法的实施例，在此不再赘述。

S72，利用最速下降法由各检波点向源点反向追踪对应的地震射线路径。对于每条地震射线，以检波点所在位置为追踪地震射线的起点，首先计算该点所在网格节点内的最速下降方向角：

$\mathrm{\θ}={\tan }^{-1}\[(\∂T/\∂y)/(\∂T/\∂x)\]+\mathrm{\π};$

式中，和分别为地震走时梯度向量的水平分量和竖直分量的中心差分格式，T_i,j为所述均匀网格化慢度模型中中索引位置(i,j)处的地震走时。最速下降方向角θ沿正水平轴方向顺时针测量，对于一个确定的源，所有穿过网格节点(i,j)的地震射线具有此相同的方向角。

为了详细说明地震射线穿过一个网格节点时的过程，我们选取了如图8所示的一个典型情况进行讨论。在图8中，地震射线从网格节点右侧边界的A点进入，A点坐标为(x_a,y_a)。按照预先分配的该网格节点内最速下降角的值，射线可能从剩下的三条边中的任何一条穿出。借助两个辅助锐角α₁＝tan^-1[(y_j+1-y_a)/h]和α₂＝tan^-1[(y_a-y_j)/h]，通过分析不同的θ的取值范围，我们可以由图8判断出穿出点的坐标以及下一个网格的索引(如表5所示)。对于地震射线路径由其它位置穿入和穿出网格节点的情况，相应的穿出点坐标和下一个网格节点索引也可用相似的方法推导出。

表5

θ的范围	穿出点坐标	下个网格索引
			0<θ≤π/2	(xi+1,yj+1)	(i+1,j+1)
π/2<θ<π-α1	(xa+(ya-yj+1)cotθ,yj+1)	(i,j+1)
			θ＝π-α1	(xi,yj+1)	(i-1,j+1)
π-α1<θ<π+α2	(xi,ya+htanθ)	(i-1,j)
			θ＝π+α2	(xi,yj)	(i-1,j-1)
π+α2<θ<3π/2	(xa-(yj-ya)cotθ,yj)	(i,j-1)
			3π/2≤θ<2π	(xi+1,yj)	(i+1,j-1)

S73，计算所述每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段的长度，并构建所述均匀网格化慢度模型的雅可比矩阵。

追踪出每条地震射线在每个网格节点内的路径(即在该网格节点内的穿入点和穿出点坐标)之后，由于网格尺寸及其坐标点已知，该地震射线路径在该网格节点内的射线段的长度可以由三角函数关系很容易求得。

对于一个具有m行n列个网格节点的离散模型，其所有离散网格节点的个数为N＝(m+1)×(n+1)，如果共有M条地震射线，则其雅可比矩阵是一个M行N列的大型稀疏矩阵，其元素可表示如下：

$L_{ij}=\frac{\∂F_i}{\∂s_j}=l_{ij},i=1,2,3,...,M,j=1,2,3,...,N.$

式中，L_ij为雅可比矩阵L的第i行、第j列个元素，F_i为第i条地震射线的地震走时，s_j为第j个网格节点的慢度值，l_ij为第i条地震射线在第j个网格节点内的射线段的长度，M为总的地震射线数，N为总的网格节点数。

S74，基于所述每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段的长度，计算每条地震射线路径的传播时间d(s)，从而得到所述均匀网格化慢度模型的走时向量。具体包括：

根据公式d(s)＝∑_ms_ml_m计算每条地震射线路径的传播时间；

式中，m是每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段标记，l_m和s_m分别为第m条射线段的长度及该射线段所在网格节点内波传播的平均慢度。

S75，根据所述雅可比矩阵、所述走时向量以及预设的迭代线性反演公式及算法对所述均匀网格化慢度模型进行迭代反演，以更新所述均匀网格化慢度模型。其中，预设的迭代线性反演公式例如可以为：

$\left[\begin{array}{c}{W}_d{L}^k\\ {\mathrm{\λW}}_m\end{array}\right](s^{k+1}-s^k)=\left[\begin{array}{c}{W}_d\[t-d\left(s^k\right)\]\\ 0\end{array}\right],k=0,1,2,3,...$

式中，W_d和W_m是两个对称矩阵，分别称为数据加权矩阵和模型加权矩阵，它们可以包含一些先验信息，这些先验信息对降低多解性非常有用，W_d通常选择单位矩阵，W_m通常选择拉普拉斯算法子(二阶差分算子)；λ为加权因子，它的作用是用来权衡数据空间和模型空间之间的比例，λ的数值可通过L曲线法来选择；L^k为第k次迭代的雅可比矩阵；d(s^k)为第k次迭代的均匀网格化慢度模型的走时向量；s^k为第k次迭代慢度向量。

而预设的例如可以为最小二乘反演算法(比如LSQR算法、BICGSTAB算法等)。

S76，重复以上步骤，直至当前更新后的均匀网格化慢度模型与所述地震走时观测值之间的关系满足设定条件，此时所述当前更新后的均匀网格化慢度模型为所求慢度模型。

其中，当前更新后的均匀网格化慢度模型与所述地震走时观测值之间的关系例如可以为：

均匀网格化慢度模型与地震走时观测值之间的误差小于设定阈值。

为了更好的反映基于高阶多方向程函方程有限差分法的井间地震走时层析成像方法的反演效果，我们对一个合成数据模型进行了地震走时层析成像反演，图9a-9d为反演结果对比。其中图9a为真实模型，图9b为基于常规的FMM算法的井间地震走时层析成像的反演结果，图9c为本申请一个实施例使用二阶多方向程函方程有限差分算法的井间地震走时层析成像的反演结果，图9d本申请另一个实施例使用三阶多方向程函方程有限差分算法的井间地震走时层析成像的反演结果。可见相比较于常规的FMM算法，基于高阶多方向程函方程有限差分算法的井间地震走时层析成像对模型的反演效果更好，不但能够较好地重建尺寸较大的异常体，对于速度模型细节的重建效果也比较好。而在反演精度上，基于高阶多方向程函方程有限差分算法的井间地震走时层析成像的误差(二阶差分为0.0586ms，三阶差分为0.0554ms)也要明显小于基于FMM算法的走时层析成像误差(0.1447ms)。这说明基于高阶多方向程函方程有限差分算法的井间地震走时层析成像能获得比基于常规FMM算法的井间走时层析成像更好的重建结果和更高的反演精度。

本申请实施例的井间地震走时层析成像方法，首先利用初至提取算法的方式获得地震剖面的初至，然后建立一个均匀初始慢度模型，接着使用高阶多方向程函方程有限差分法计算慢度模型的地震走时，再用最速下降法追踪射线路径，进而构建出雅可比矩阵，最后进行走时层析成像来反演慢度场，得到一个新的慢度模型。通过迭代线性反演算法不断更新慢度模型并计算误差，最终得到一个误差足够小的慢度模型来对井间慢度场进行描述，从而获得比基于常规FMM算法的井间地震走时层析成像更好的重建结果和更高的反演精度。

本领域技术人员还可以了解到，本申请实施例列出的各种说明性步骤可以通过软件实现，也可以用硬件来实现，或者硬件、软件两者的结合来实现，至于是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本申请实施例保护的范围。

以上所述的具体实施例，对本申请的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本申请实施例的具体实施例而已，并不用于限定本申请的保护范围，凡在本申请的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims (10)

1.一种地震走时获取方法，其特征在于，包括以下步骤：

拾取地震记录的初至，得到其地震走时观测值；

依据所述地震走时观测值建立均匀网格化慢度模型；

确定所述均匀网格化慢度模型中源点邻点的地震走时；

基于高阶多方向程函方程有限差分法确定所述均匀网格化慢度模型中除所述源点邻点外其他网格节点的地震走时。

2.根据权利要求1所述的地震走时获取方法，其特征在于，所述依据所述地震走时观测值建立均匀网格化慢度模型具体包括：

根据公式获取预设初始慢度模型的速度值v₀，其中，M为地震射线条数，l_i为第i条地震射线的长度，t_i为拾取的第i条地震射线的地震走时观测值；

取所述预设初始慢度模型的速度值v₀的倒数，获得所述预设初始慢度模型的慢度值s₀；

使用正方形网格将所述预设初始慢度模型均匀划分为多个网格节点，并将每个网格节点的慢度值赋值为s₀，得到均匀网格化慢度模型。

3.根据权利要求1所述的地震走时获取方法，其特征在于，所述确定所述均匀网格化慢度模型中源点邻点的地震走时，具体包括：

在所述均匀网格化慢度模型中：

将震源所在位置作为源点，将源点标记为近点，并将其地震走时赋值为0；

将除所述近点之外的所有网格节点标记为远点并将其赋值为预设的地震走时初值；

根据一阶程函方程有限差分公式 $max{\left(\frac{1}{h}\right({T^{\′}}_{i,j}-{T^{\′}}_1),0)}^2+max\left(\frac{1}{h}\right({T^{\′}}_{i,j}-$ ${T^{\′}}_2),0{)}^2=S_{i,j}^{\′2}$ 获取所述近点的邻点的地震走时；

式中，T′_i,j为邻点处的地震走时，下角标i和j为该邻点在慢度模型内的索引位置，S′_i,j为邻点处的慢度值，T′₁＝min(T′_i-1,j,T′_i+1,j)为索引位置(i-1,j)和(i+2,j)处地震走时的最小值，T′₂＝min(T′_i,j-1,T′_i,j+1)为索引位置(i,j-1)和(i,j+2)处地震走时的最小值，h为网格间距。

4.根据权利要求3所述的地震走时获取方法，其特征在于，所述基于高阶多方向程函方程有限差分法确定所述均匀网格化慢度模型中除所述源点邻点外其他网格节点的地震走时，具体包括：

将当前近点的所有邻点放入预设的波前窄带中并标记为窄带点；

在所述波前窄带中选取其中地震走时最小的窄带点并将其移出所述波前窄带作为新的近点；

根据公式 $max{\left(\frac{1}{\mathrm{\Δ}d}\right(T_{i,j}-T_1),0)}^2+max{\left(\frac{1}{\mathrm{\Δ}d}\right(T_{i,j}-T_2),0)}^2=S_{i,j}^2$ 获取或更新与所述新的近点相邻的窄带点或远点的地震走时，式中，T_i,j为与所述新的近点相邻的窄带点或远点的地震走时，S_i,j为其对应的慢度值，T₁和T₂为近点的地震走时，Δd为节点间距；

重复以上步骤，直至所述波前窄带中为空时停止，从而计算出所述均匀网格化慢度模型中除所述源点邻点外其他网格节点的地震走时。

5.一种基于权利要求1所述地震走时获取方法的井间地震走时层析成像方法，其特征在于，包括以下步骤：

基于高阶多方向程函方程有限差分法确定均匀网格化慢度模型中所有网格节点的地震走时；

利用最速下降法由各检波点向源点反向追踪对应的地震射线路径；

计算所述每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段的长度，并构建所述均匀网格化慢度模型的雅可比矩阵；

基于所述每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段的长度，计算每条地震射线路径的传播时间，得到所述均匀网格化慢度模型的走时向量；

根据所述雅可比矩阵、所述走时向量以及预设的迭代线性反演公式及算法对所述均匀网格化慢度模型进行迭代反演，以更新所述均匀网格化慢度模型；

重复以上步骤，直至当前更新后的均匀网格化慢度模型与所述地震走时观测值之间的关系满足设定条件，此时所述当前更新后的均匀网格化慢度模型为所求慢度模型。

6.根据权利要求5所述的井间地震走时层析成像方法，其特征在于，所述当前更新后的均匀网格化慢度模型与所述地震走时观测值之间的关系包括：

均匀网格化慢度模型与地震走时观测值之间的误差小于设定阈值。

7.根据权利要求5所述的井间地震走时层析成像方法，其特征在于，所述利用最速下降法由各检波点向源点反向追踪对应的地震射线路径中，每条地震射线路径在所述均匀网格化慢度模型中各个网格节点的最速下降方向角通过以下公式计算得到：

$\mathrm{\θ}={\tan }^{-1}[(\∂T/\∂y)/(\∂T/\∂x)]+\mathrm{\π};$

式中， $\frac{\∂T}{\∂x}=\frac{\[(T_{i+1,j}+T_{i+1,j+1})-(T_{i,j}+T_{i,j+1})\]}{2h}$ 和 $\frac{\∂T}{\∂y}=\frac{\[(T_{i,j+1}+T_{i+1,j+1})-(T_{i,j}+T_{i+1,j})\]}{2h}$ 分别为走时梯度向量的水平分量和竖直分量的中心差分格式，T_i,j为所述均匀网格化慢度模型中中索引位置(i,j)处的地震走时。

8.根据权利要求5所述的井间地震走时层析成像方法，其特征在于，基于所述每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段的长度，计算每条地震射线路径的传播时间，具体包括：

根据公式d(s)＝∑_ms_ml_m计算每条地震射线路径的传播时间；

式中，m是每条地震射线路径在所述均匀网格化慢度模型中各个网格节点内的射线段标记，l_m和s_m分别为第m条射线段的长度及该射线段所在网格节点内波传播的平均慢度。

9.根据权利要求5所述的井间地震走时层析成像方法，其特征在于，所述预设的迭代线性反演公式包括：

$\left[\begin{array}{c}{W}_d{L}^k\\ {\mathrm{\λW}}_m\end{array}\right](s^{k+1}-s^k)=\left[\begin{array}{c}{w}_d\[t-d\left(s^k\right)\]\\ 0\end{array}\right],k=0,1,2,3,...$

式中，W_d和W_m是两个对称矩阵，λ为加权因子，L^k为第k次迭代的雅可比矩阵，d(s^k)为第k次迭代的均匀网格化慢度模型的走时向量，s^k为第k次迭代慢度向量。

10.根据权利要求5所述的井间地震走时层析成像方法，其特征在于，所述预设的迭代线性反演算法包括最小二乘反演算法。