CN114839675A - 一种建立三维速度模型的方法 - Google Patents

Abstract

本发明属于地震资料数据处理技术领域，提供了一种建立三维速度模型的方法，包括：对采集的原始三维地震资料进行预处理；在共炮点道集和共检波点道集拾取局部相关反射波走时和斜率数据，形成观测数据空间；建立离散速度模型，形成模型空间，并给定初始离散速度模型；在当前速度模型中，用拾取的斜率数据分别计算从炮点和检波点向地下传播的初始射线方向，通过求解三维程函方程获得分别从炮点与检波点的射线路径、走时和对应的射线参数；建立反演目标函数及反演方程，求解获得对速度模型的修正量，更新速度模型；判断当前速度模型是否符合要求，若是，则输出当前速度模型，本发明效率高、对复杂地质构造和低信噪比资料有很强适应性。

Description

一种建立三维速度模型的方法

技术领域

本发明属于地震资料数据处理技术领域，特别涉及一种建立三维速度模型的方法。

背景技术

地震勘探资料叠前深度偏移成像效果极大地依赖于所用的地下速度模型的准确度。地震走时层析反演是一种利用地震波走时资料建立地下低频速度模型的有效技术。常规反射波走时层析反演需要连续追踪各个同相轴，并拾取各个同相轴的走时数据，特别是需要事先建立这些同相轴与地下地层界面之间的对应关系，这对于复杂地质构造地区或低信噪比地震资料，往往难以实现。此外，反射波走时层析反演中，需要通过射线追踪计算当前速度模型的地震波走时以及反演方程中的核函数矩阵。确定一条从炮点到检波点的射线路径，除了要知道速度模型、炮点和检波点位置外，还需要知道射线的出射角和入射角。常规反射波走时层析反演通过若干次试验确定出射角和入射角，或者采用先求波前面再由波前面确定射线的方式，使射线追踪的计算量很大，导致走时层析反演的效率很低。

发明的内容

本发明的目的在于至少解决现有技术中存在的技术问题之一，提供一种使用地震波走时和斜率数据建立三维速度模型的方法，该方法易于计算机自动化实现，为复杂构造地区的地震成像提供了一种有效的三维速度建模方法。

本发明一种实施例提供的一种建立三维速度模型的方法，包括如下步骤：

步骤一，对原始三维地震资料进行预处理；

步骤二，拾取炮点和检波点的地震走时和斜率数据，形成观测数据空间，具体包括：利用基于曲波变换的局部相关同相轴的斜率和走时拾取算法，分别在共炮点道集和共检波点道集，追踪局部相关同相轴，并拾取相应的走时和斜率数据，对三维地震资料分别沿主测线(x方向)和交叉测线(y方向)进行操作后，获得如下的观测数据空间：

其中，T_{SR_obs}为观测的从炮点S到检波点R的地震波反射走时；P_{SX_obs}和P_{SY_obs}分别为观测的在炮点处沿x和y方向的地震斜率；P_{RX_obs},、P_{RY_obs}分别为观测的在检波点处沿x和y方向的地震斜率，N为拾取数据的炮点检波点对个数，即炮检对个数。

步骤三，建立离散速度模型，形成模型空间，并给定初始离散速度模型，具体包括：根据观测***资料的高程数据，确定起伏的观测面，再根据研究区的范围，确定三维速度模型在x、y和z方向的尺寸；根据地下构造的复杂程度和拟建立的速度模型分辨率要求，确定离散网格大小；对所述离散速度模型进行离散化，设置各个网格节点的速度值，所述离散速度模型中任意一点的速度值用所述网格节点速度值的三次样条函数表示，各个网格节点的速度值为模型参量，即未知量，所述模型空间为：

其中，V_j为第j个网格节点的速度值，M为网格节点数，给各网格点的速度赋值，形成初值速度模型；

步骤四，根据拾取的炮点斜率数据确定射线在炮点和在检波点的出射方向，即获得射线参数的初值，然后通过求解三维程函方程，获得当前速度模型的射线路径、走时和对应的射线参数数据，具体包括：

使用所述拾取的炮点斜率数据，用下式计算射线在所述炮点的出射方向角，以及所述炮点处的射线参数垂直分量：

其中，

为炮点射线方向与过x轴平面的夹角；θ_S为射线方向在过x轴平面的投影与z轴的夹角；V_0S为炮点处的速度；P_{SX_obs}和P_{SY_obs}分别为拾取的炮点处沿x和y方向的地震斜率；P_SZ0为炮点处垂直慢度分量；sin^-1表示反正弦函数；cos表示余弦函数，炮点处的初始射线参数P_S0为

P_S0＝[P_{SX_obs},P_{SY_obs},P_SZ0]

使用拾取的检波点斜率数据，用下式计算射线在检波点的出射方向角，以及检波点处的射线参数垂直分量：

其中，

为检波点射线方向与过x轴平面的夹角；θ_R为射线方向在过x轴平面的投影与z轴的夹角；V_0R为检波点处的速度；P_{RX_obs}和P_{RY_obs}分别为拾取的检波点处沿x和y方向的地震斜率；P_RZ0为检波点处垂直慢度分量；sin^-1表示反正弦函数；cos表示余弦函数。这时，检波点处的初始射线参数P_R0为

P_R0＝[P_{RX_obs},P_{RY_obs},P_RZ0]

分别从炮点S和检波点R处开始，使用在所述步骤S3和S4求得的炮点和检波点初始射线参数，采用Runge-Kutta法求解下列三维程函方程，

其中，x＝(X,Y,Z)为空间位置向量；p(x)＝[P_X(x),P_Y(x),P_Z(x)]为由三个坐标轴方向的射线参数或慢度分量组成的向量；V(x)为x处的速度；t为从炮点或检波点的地震波走时；▽表示梯度算子。

求解上述微分方程，不断求得从炮点S向地下传播的射线路径和单程走时t_S，以及从检波点R向地下传播的射线路径和单程走时t_R，当t_S与t_R之和等于该炮检对的观测双程走时T_{SR_obs}时，停止计算，这时炮点S的射线端点位置记为x_S，检波点R的射线端点位置记为x_R，存储炮点和检波点的射线在相应的各个网格单元的射线长度、射线参数；

步骤五，计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差，判断当前速度模型是否满足要求，确定反演过程是否终止。

进一步，所述步骤五计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差，判断当前速度模型是否满足要求，确定反演过程是否终止，具体包括以下步骤：

S1，如果x_S与x_R重合，且分别沿炮点和检波点的射线的单程走时之和与观测的双程反射走时一致时，说明所用的所述离散速度模型正确，用炮检对的炮点和检波点射线末端位置的距离表示从炮点和检波点确定的反射点误差，用下列式子分别计算反射点误差及反射走时误差：

其中，(X_Si，Y_Si，Z_Si)和(X_Ri，Y_Ri，Z_Ri)分别为从第i个炮检对的炮点和检波点的射线末端的坐标；t_Si和t_Ri分别为从第i个炮检对的炮点和检波点到各自射线末端的走时；T_{SRi obs}为观测的第i个炮检对的反射波走时；

S2，给定误差限，当计算的err_D和err_T小于给定的误差限时，输出当前速度模型，停止运行；否则，执行下列步骤；

S3，建立特定的目标函数和反演方程，求解得到对当前速度模型的修正量，更新速度模型；

S4，建立下列目标函数：

目

D_SRi＝[(X_Si-X_Ri)²+(Y_Si-Y_Ri)²+(Z_Si-Z_Ri)²]^1/2

T_SRi＝t_Si+t_Ri-T_{SRi_obs}

其中，N为炮检对的个数；t_Si和t_Ri分别为从第i个炮检对的炮点和检波点到各自射线末端的走时；T_{SRi_obs}为观测的第i个炮检对的反射波双程走时；T_SRi为第i个炮检对的两个单程走时之和与观测的双程走时之差；x_Si＝(X_Si，Y_Si，Z_Si)和x_Ri＝(X_Ri，Y_Ri，Z_Ri)分别为第i个炮检对的炮点和检波点的射线末端的坐标向量；D_SRi为射线末端点x_Si与x_Ri之间的距离；v为由网格节点速度组成的模型参数向量；λ为权系数；

S5，采用高斯-牛顿法求解上述非线性目标函数的极小值解，将所述目标函数进行泰勒级数展开取至一阶项，得到线性反演方程如下：

G·Δv＝Δd

其中，

其中，i＝1，2，...，N，j＝1，2，...，M；

B_ij＝λ(L_Sij+L_Rij)，i＝1，2，...，N，j＝1，2，...，M；

Δd＝[D_SR1 D_SR2 … D_SRN T_SR1 T_SR2 … T_SRN]^T；

Δv＝[V₁ V₂ … V_M]^T；

D_sRi＝[(X_Si-X_Ri)²+(Y_Si-Y_Ri)²+(Z_Si-Z_Ri)²]^1/2，i＝1，2，...，N；

T_SRi＝t_Si+t_Ri-T_{SRi_obs}，i＝1，2，...，N；

t_Si和t_Ri分别为从第i个炮检对的炮点和检波点到各自射线末端的走时；T_{SRi_obs}为观测的第i个炮检对的反射波走时；(X_Si，Y_Si，Z_Si)和(X_Ri，Y_Ri，Z_Ri)分别为第i个炮检对的炮点和检波点的射线末端的坐标；V_i是第i个离散网格的速度值；V_Si和V_Ri分别为第i个炮检对的炮点和检波点射线末端所在单元的速度；P_SXi，P_SYi和P_SZi分别为第i个炮检对的炮点射线末端沿x，y和z方向的射线参数；P_RXi，P_RYi和P_RZi分别为第i个炮检对的检波点射线末端沿x，y和z方向的射线参数；L_Sij和L_Rij分别为第i个炮检对的炮点射线和检波点射线在第j个单元的射线段长度；L_SXij，L_SYij和LS_Zij分别为第i个炮检对的炮点射线在第j个单元的射线段在x，y和z方向的投影长度；L_RXij，L_RYij和L_RZij分别为第i个炮检对的检波点点射线在第j个单元的射线段在x，y和z方向的投影长度；N为炮-检对总数；M为离散模型的网格节点总数；

S6，使用LSQR算法求解上述线性反演方程，得到离散速度模型更新量Δv，用该模型更新量修改当前离散速度模型，就得到修正后的模型，计算公式如下：

v^(k+1)＝v^(k)+Δv

其中，v^(k)是当前离散速度模型，即上次(第k次)迭代的离散速度模型，Δv是本次迭代反演得到的离散速度模型修正量；v^(k+1)是本次迭代得到的修正后的离散速度模型，再把修正后的离散速度模型作为新的初始模型，转到步骤四。

本发明的有益效果是：采用本发明所述的使用地震反射波走时和斜率数据建立三维速度模型的方法，用局部相关反射波同相轴的斜率数据计算从炮点和检波点向地下传播方向和初始射线参数，提高了计算射线路径和走时的效率；使用的局部相关反射波同相轴易于自动化追踪，其走时和斜率数据易于自动化拾取，对复杂地下构造和低信噪比地震资料有很强的适应性。本发明实施例易于计算机自动化实现，为复杂构造地区的地震成像提供一种有效的三维速度建模方法。

附图说明

下面结合附图和实施例对本发明进一步地说明；

图1为本发明实施例一种建立三维速度模型的方法流程图；

图2为图1所示实施例利用走时和斜率数据建立的三维速度模型切片图；

图3为常规利用走时数据的层析反演方法建立的三维速度模型切片图；

图4为使用本发明实施例建立的三维速度模型进行偏移后的共成像点道集图；

图5为使用常规走时层析反演建立的三维速度模型进行偏移后的共成像点道集图。

具体实施方式

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

在地震数据中，共炮道集和共检波点道集容易追踪局部相关同相轴，局部相关同相轴的走时随偏移距的变化率被称为地震斜率。根据射线理论，地震斜率是射线参数，也即慢度矢量的水平分量。如果已知炮点和检波点处的速度值，就可用斜率参数计算出反射射线在炮点的出射角和在检波点的入射角。而且，局部相干同相轴易于追踪，其走时和斜率也易于拾取，对低信噪比地震资料有很强的适应性。若采用光滑速度模型，通过地震斜率对射线进行约束，只需要局部相关的同相轴走时和斜率数据，而不需要事先建立各个同相轴与地层界面之间的对应关系，从而提高反射波走时层析反演的自动化程度和实用性。因此，如果把地震斜率数据引入到反射波走时层析反演中，用地震斜率确定射线出射角和入射角，就可以提高射线追踪的效率，进而形成一种实用化的建立地下低频速度模型的方法，为叠前深度偏移速度分析和波形反演提供低频或初始速度模型。

本发明实施例把地震斜率数据引入到三维地震反射走时层析反演中，用地震斜率确定射线出射角和入射角，根据当速度模型准确时，分别从炮点和检波点向地下传播的射线将相遇在反射点，且两者的走时之和必等于实测的双程反射走时的情况，提出了一种使用地震反射波走时和斜率数据建立三维速度模型的方法，为复杂构造地区地震勘探或低信噪比地震资料的偏移成像，提供一种建立地下三维速度模型的技术。

参照如图1所示，为本发明实施例一种建立三维速度模型的方法，包括如下步骤：

步骤一，对原始三维地震资料进行预处理；

步骤二，拾取炮点和检波点的地震走时和斜率数据，形成观测数据空间，具体包括：利用基于曲波变换的局部相关同相轴的斜率和走时拾取算法，分别在共炮点道集和共检波点道集，追踪局部相关同相轴，并拾取相应的走时和斜率数据，对三维地震资料分别沿主测线(x方向)和交叉测线(y方向)进行操作后，获得如下的观测数据空间：

其中，T_{SR_obs}为观测的从炮点S到检波点R的地震波反射走时；P_{SX_obs}和P_{SY_obs}分别为观测的在炮点处沿x和y方向的地震斜率；P_{RX_obs},、P_{RY_obs}分别为观测的在检波点处沿x和y方向的地震斜率，N为拾取数据的炮点检波点对个数，即炮检对个数。

步骤三，建立离散速度模型，形成模型空间，并给定初始离散速度模型，具体包括：根据观测***资料的高程数据，确定起伏的观测面，再根据研究区的范围，确定三维速度模型在x、y和z方向的尺寸；根据地下构造的复杂程度和拟建立的速度模型分辨率要求，确定离散网格大小；对所述离散速度模型进行离散化，设置各个网格节点的速度值，所述离散速度模型中任意一点的速度值用所述网格节点速度值的三次样条函数表示，各个网格节点的速度值为模型参量，即未知量，所述模型空间为：

其中，V_j为第j个网格节点的速度值，M为网格节点数，给各网格点的速度赋值，形成初值速度模型；

步骤四，根据拾取的炮点斜率数据确定射线在炮点和在检波点的出射方向，即获得射线参数的初值，然后通过求解三维程函方程，获得当前速度模型的射线路径、走时和对应的射线参数数据，具体包括：

使用所述拾取的炮点斜率数据，用下式计算射线在所述炮点的出射方向角，以及所述炮点处的射线参数垂直分量：

其中，

为炮点射线方向与过x轴平面的夹角；θ_S为射线方向在过x轴平面的投影与z轴的夹角；V_0S为炮点处的速度；P_{SX_obs}和P_{SY_obs}分别为拾取的炮点处沿x和y方向的地震斜率；P_SZ0为炮点处垂直慢度分量；sin^-1表示反正弦函数；cos表示余弦函数，炮点处的初始射线参数P_S0为

P_S0＝[P_{SX_obs},P_{SY_obs},P_SZ0]

使用拾取的检波点斜率数据，用下式计算射线在检波点的出射方向角，以及检波点处的射线参数垂直分量：

其中，

为检波点射线方向与过x轴平面的夹角；θ_R为射线方向在过x轴平面的投影与z轴的夹角；V_0R为检波点处的速度；P_{RX_obs}和P_{RY_obs}分别为拾取的检波点处沿x和y方向的地震斜率；P_RZ0为检波点处垂直慢度分量；sin^-1表示反正弦函数；cos表示余弦函数。这时，检波点处的初始射线参数P_R0为

P_R0＝[P_{RX_obs},P_{RY_obs},P_RZ0]

分别从炮点S和检波点R处开始，使用在所述步骤S3和S4求得的炮点和检波点初始射线参数，采用Runge-Kutta法求解下列三维程函方程，

其中，x＝(X,Y,Z)为空间位置向量；p(x)＝[P_X(x),P_Y(x),P_Z(x)]为由三个坐标轴方向的射线参数或慢度分量组成的向量；V(x)为x处的速度；t为从炮点或检波点的地震波走时；

表示梯度算子。

求解上述微分方程，不断求得从炮点S向地下传播的射线路径和单程走时t_S，以及从检波点R向地下传播的射线路径和单程走时t_R，当t_S与t_R之和等于该炮检对的观测双程走时T_{SR_obs}时，停止计算，这时炮点S的射线端点位置记为x_S，检波点R的射线端点位置记为x_R，存储炮点和检波点的射线在相应的各个网格单元的射线长度、射线参数；

步骤五，计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差，判断当前速度模型是否满足要求，确定反演过程是否终止。

进一步，所述步骤五计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差，判断当前速度模型是否满足要求，确定反演过程是否终止，具体包括以下步骤：

S1，如果x_S与x_R重合，且分别沿炮点和检波点的射线的单程走时之和与观测的双程反射走时一致时，说明所用的所述离散速度模型正确，用炮检对的炮点和检波点射线末端位置的距离表示从炮点和检波点确定的反射点误差，用下列式子分别计算反射点误差及反射走时误差：

其中，(X_Si，Y_Si，Z_Si)和(X_Ri，Y_Ri，Z_Ri)分别为从第i个炮检对的炮点和检波点的射线末端的坐标；t_Si和t_Ri分别为从第i个炮检对的炮点和检波点到各自射线末端的走时；T_{SRi_obs}为观测的第i个炮检对的反射波走时；

S2，给定误差限，当计算的err_D和err_T小于给定的误差限时，输出当前速度模型，停止运行；否则，执行下列步骤；

S3，建立特定的目标函数和反演方程，求解得到对当前速度模型的修正量，更新速度模型；

S4，建立下列目标函数：

且

D_SRi＝[(X_Si-X_Ri)²+(Y_Si-Y_Ri)²+(Z_Si-Z_Ri)²]^1/2

T_SRi＝t_Si+t_Ri-T_{SRi_obs}

其中，N为炮检对的个数；t_Si和t_Ri分别为从第i个炮检对的炮点和检波点到各自射线末端的走时；T_{SRi_obs}为观测的第i个炮检对的反射波双程走时；T_SRi为第i个炮检对的两个单程走时之和与观测的双程走时之差；x_Si＝(X_Si，Y_Si，Z_Si)和x_Ri＝(X_Ri，Y_Ri，Z_Ri)分别为第i个炮检对的炮点和检波点的射线末端的坐标向量；D_SRi为射线末端点x_Si与x_Ri之间的距离；v为由网格节点速度组成的模型参数向量；λ为权系数；

S5，采用高斯-牛顿法求解上述非线性目标函数的极小值解，将所述目标函数进行泰勒级数展开取至一阶项，得到线性反演方程如下：

G·Δv＝Δd

其中，

其中，i＝1，2，...，N，j＝1，2，...，M；

B_ij＝λ(L_Sij+L_Rij)，i＝1，2，...，N，j＝1，2，...，M；

Δd＝[D_SR1 D_SR2 … D_SRN T_SR1 T_SR2 … T_SRN]^T；

Δv＝[V₁ V₂ … V_M]^T；

D_SRi＝[(X_Si-X_Ri)²+(Y_Si-Y_Ri)²+(Z_Si-Z_Ri)²]^1/2，i＝1，2，...，N；

T_SRi＝t_Si+t_Ri-T_{SRi_obs}，i＝1，2，...，N；

t_Si和t_Ri分别为从第i个炮检对的炮点和检波点到各自射线末端的走时；T_{SRi_obs}为观测的第i个炮检对的反射波走时；(X_Si，Y_Si，Z_Si)和(X_Ri，Y_Ri，Z_Ri)分别为第i个炮检对的炮点和检波点的射线末端的坐标；V_i是第i个离散网格的速度值；V_Si和V_Ri分别为第i个炮检对的炮点和检波点射线末端所在单元的速度；P_SXi，P_SYi和P_SZi分别为第i个炮检对的炮点射线末端沿x，y和z方向的射线参数；P_RXi，P_RYi和P_RZi分别为第i个炮检对的检波点射线末端沿x，y和z方向的射线参数；L_Sij和L_Rij分别为第i个炮检对的炮点射线和检波点射线在第j个单元的射线段长度；L_SXij，L_sYij和L_SZij分别为第i个炮检对的炮点射线在第j个单元的射线段在x，y和z方向的投影长度；L_RXij，L_RYij和L_RZij分别为第i个炮检对的检波点点射线在第j个单元的射线段在x，y和z方向的投影长度；N为炮-检对总数；M为离散模型的网格节点总数；

S6，使用LSQR算法求解上述线性反演方程，得到离散速度模型更新量Δv，用该模型更新量修改当前离散速度模型，就得到修正后的模型，计算公式如下：

v^(k+1)＝v^(k)+Δv

其中，v^(k)是当前离散速度模型，即上次(第k次)迭代的离散速度模型，Δv是本次迭代反演得到的离散速度模型修正量；v^(k+1)是本次迭代得到的修正后的离散速度模型，再把修正后的离散速度模型作为新的初始模型，转到步骤四。

作为本发明实施例的一个具体示例，用我国西南某地区的三维油气地震勘探资料对本发明实施例方法进行了测试，并与仅用反射走时数据的反射层析成像的结果进行比较。

步骤100，本工区长16km，宽8km。均匀布设25条炮线和59条检波点线，炮线间距300m，检波线间距150m，炮间距100m，检波点间距50m。每炮激发时，8条检波线且每条检波线240个检波点接收。对原始地震资料进行了压制噪声、振幅补偿、反褶积、静校正等预处理。

步骤110，利用基于曲波变换的局部相关同相轴斜率和走时拾取方法，分别沿检波线和炮线方向，在共炮点域和共检波点域拾取走时和斜率数据，获得236152个炮-检对的反射走时和斜率数据，形成数据空间。

步骤120，建立大小为16km×8km×4km的速度模型，并用250m×250m×250m的网格离散化该三维速度模型，形成模型空间。模型表面用对应的高程数据构建起伏地形。模型初始速度设置成随深度线性增加，表面最高处速度为1.0km/s，模型底界处速度4.5km/s，形成初始速度模型。

步骤130，利用初始速度和拾取的炮点、检波点的斜率数据，获得分别从炮点和检波点向地下传播的地震射线方向，给出射线出射的初始条件。

步骤140，根据初始速度模型和射线初始条件，采用Runge-Kutta法求解三维程函方程，得到各个炮-检对的射线末端位置、走时以及对应的各个网格单元的射线路径和射线参数。

步骤150，用各个炮-检对的炮点和检波点射线末端位置的距离计算反射点误差，用各个炮-检对的炮点和检波点射线末端走时及实测的反射波走时计算反射走时误差。

步骤160，把计算的反射点误差和反射走时误差分别与事先设置的误差限对比，判断当前速度模型是否满足要求。若计算的反射点误差和反射走时误差都小于设置的误差限，则输出当前速度模型，终止计算；否则，进行下列步骤。

步骤170，利用观测的反射走时和斜率数据、在S5步计算的各个炮-检对的射线末端位置、走时以及对应的各个网格单元的射线路径和射线参数，建立反演方程。求解反演方程，得到对当前速度模型的修正量。

步骤180，用求得的模型修正量对当前模型进行修正，得到更新的模型。把更新模型作为新的初始模型，转到步骤S4。

图2是本发明在迭代30次后输出的地下速度模型沿中间测线的垂直切片,图3是仅用反射走时数据建立的地下速度模型沿中间测线的垂直切片。对比这两个模型发现，本发明使用走时和斜率数据建立的速度模型与常规的仅用走时数据建立的速度模型，在横向上和在垂直方向上都有较大的差别。

为了进一步说明本发明的有益效果，分别利用本发明使用走时与斜率数据建立的速度模型和常规的仅用走时数据建立的速度模型，对地震数据进行叠前深度偏移成像。图4所示是用本发明建立的速速模型进行偏移，得到沿中间测线的共成像点道集，图5所示是用常规走时层析反演建立的速度模型进行偏移，得到的沿中间测线的共成像点道集。图5中的同相轴下拉与上翘现象较为明显，同相轴基本都呈现非水平形态，而图4中的同相轴则相对较为水平和光滑。在共成像点道集中，常常通过同相轴的水平程度来判断偏移时所用速度模型的准确性。同相轴越趋于水平，说明偏移所使用的速度模型越准确。比较图4与图5可见，本发明建立的速度模型更准确，表明了本发明的良好应用效果。

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

Claims (2)

1.一种利用地震数据建立三维速度模型的方法，其特征在于，包括以下步骤：

步骤一，对原始三维地震资料进行预处理；

步骤二，拾取炮点和检波点的地震走时和斜率数据，形成观测数据空间，具体包括：利用基于曲波变换的局部相关同相轴的斜率和走时拾取算法，分别在共炮点道集和共检波点道集，追踪局部相关同相轴，并拾取相应的走时和斜率数据，对三维地震资料分别沿主测线(x方向)和交叉测线(y方向)进行操作后，获得如下的观测数据空间：

其中，T_{SR_obs}为观测的从炮点S到检波点R的地震波反射走时；P_{SX_obs}和P_{SY_obs}分别为观测的在炮点处沿x和y方向的地震斜率；P_{RX_obs},、P_{RY_obs}分别为观测的在检波点处沿x和y方向的地震斜率，N为拾取数据的炮点检波点对个数，即炮检对个数。

步骤三，建立离散速度模型，形成模型空间，并给定初始离散速度模型，具体包括：根据观测***资料的高程数据，确定起伏的观测面，再根据研究区的范围，确定三维速度模型在x、y和z方向的尺寸；根据地下构造的复杂程度和拟建立的速度模型分辨率要求，确定离散网格大小；对所述离散速度模型进行离散化，设置各个网格节点的速度值，所述离散速度模型中任意一点的速度值用所述网格节点速度值的三次样条函数表示，各个网格节点的速度值为模型参量，即未知量，所述模型空间为：

其中，V_j为第j个网格节点的速度值，M为网格节点数，给各网格点的速度赋值，形成初值速度模型；

步骤四，根据拾取的炮点斜率数据确定射线在炮点和在检波点的出射方向，即获得射线参数的初值，然后通过求解三维程函方程，获得当前速度模型的射线路径、走时和对应的射线参数数据，具体包括：

使用所述拾取的炮点斜率数据，用下式计算射线在所述炮点的出射方向角，以及所述炮点处的射线参数垂直分量：

其中，

为炮点射线方向与过x轴平面的夹角；θ_S为射线方向在过x轴平面的投影与z轴的夹角；V_0S为炮点处的速度；P_{SX_obs}和P_{SY_obs}分别为拾取的炮点处沿x和y方向的地震斜率；P_SZ0为炮点处垂直慢度分量；sin^-1表示反正弦函数；cos表示余弦函数，炮点处的初始射线参数P_S0为

P_S0＝[P_{SX_obs},P_{SY_obs},P_SZ0]

使用拾取的检波点斜率数据，用下式计算射线在检波点的出射方向角，以及检波点处的射线参数垂直分量：

其中，

为检波点射线方向与过x轴平面的夹角；θ_R为射线方向在过x轴平面的投影与z轴的夹角；V_0R为检波点处的速度；P_{RX_obs}和P_{RY_obs}分别为拾取的检波点处沿x和y方向的地震斜率；P_RZ0为检波点处垂直慢度分量；sin^-1表示反正弦函数；cos表示余弦函数。这时，检波点处的初始射线参数P_R0为

P_R0＝[P_{RX_obs},P_{RY_obs},P_RZ0]

分别从炮点S和检波点R处开始，使用在所述步骤S3和S4求得的炮点和检波点初始射线参数，采用Runge-Kutta法求解下列三维程函方程，

其中，x＝(X,Y,Z)为空间位置向量；p(x)＝[P_X(x),P_Y(x),P_Z(x)]为由三个坐标轴方向的射线参数或慢度分量组成的向量；V(x)为x处的速度；t为从炮点或检波点的地震波走时；▽表示梯度算子。

求解上述微分方程，不断求得从炮点S向地下传播的射线路径和单程走时t_S，以及从检波点R向地下传播的射线路径和单程走时t_R，当t_S与t_R之和等于该炮检对的观测双程走时T_{SR_obs}时，停止计算，这时炮点S的射线端点位置记为x_S，检波点R的射线端点位置记为x_R，存储炮点和检波点的射线在相应的各个网格单元的射线长度、射线参数；

步骤五，计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差，判断当前速度模型是否满足要求，确定反演过程是否终止。

2.根据权利要求1所述的一种使用地震波走时和斜率数据建立三维速度模型的方法，其特征在于，所述步骤五计算分别从所述炮点和检波点确定的反射点误差及观测反射走时与计算反射走时的误差，判断当前速度模型是否满足要求，确定反演过程是否终止，具体包括以下步骤：

S1，如果x_S与x_R重合，且分别沿炮点和检波点的射线的单程走时之和与观测的双程反射走时一致时，说明所用的所述离散速度模型正确，用炮检对的炮点和检波点射线末端位置的距离表示从炮点和检波点确定的反射点误差，用下列式子分别计算反射点误差及反射走时误差：

其中，(X_Si,Y_Si,Z_Si)和(X_Ri,Y_Ri,Z_Ri)分别为从第i个炮检对的炮点和检波点的射线末端的坐标；t_Si和t_Ri分别为从第i个炮检对的炮点和检波点到各自射线末端的走时；T_{SRi_obs}为观测的第i个炮检对的反射波走时；

S2，给定误差限，当计算的err_D和err_T小于给定的误差限时，输出当前速度模型，停止运行；否则，执行下列步骤；

S3，建立特定的目标函数和反演方程，求解得到对当前速度模型的修正量，更新速度模型；

S4，建立下列目标函数：

且

D_SRi＝[(X_Si-X_Ri)²+(Y_Si-Y_Ri)²+(Z_Si-Z_Ri)²]^1/2

T_SRi＝t_Si+t_Ri-T_{SRi_obs}

其中，N为炮检对的个数；t_Si和t_Ri分别为从第i个炮检对的炮点和检波点到各自射线末端的走时；T_{SRi_obs}为观测的第i个炮检对的反射波双程走时；T_SRi为第i个炮检对的两个单程走时之和与观测的双程走时之差；x_Si＝(X_Si,Y_Si,Z_Si)和x_Ri＝(X_Ri,Y_Ri,Z_Ri)分别为第i个炮检对的炮点和检波点的射线末端的坐标向量；D_SRi为射线末端点x_Si与x_Ri之间的距离；v为由网格节点速度组成的模型参数向量；λ为权系数；

S5，采用高斯-牛顿法求解上述非线性目标函数的极小值解，将所述目标函数进行泰勒级数展开取至一阶项，得到线性反演方程如下：

G·Δv＝Δd

其中，

其中，i＝1,2,…,N,j＝1,2,…,M；

B_ij＝λ(LS_ij+LR_ij)，i＝1,2,…,N,j＝1,2,…,M；

Δd＝[D_SR1 D_SR2 … D_SRN T_SR1 T_SR2 … T_SRN]^T；

Δv＝[V₁ V₂ … V_M]^T；

D_SRi＝[(X_Si-X_Ri)²+(Y_Si-Y_Ri)²+(Z_Si-Z_Ri)²]^1/2,i＝1,2,…,N；

T_SRt＝t_St+t_Rt-T_{SRt_obs},i＝1,2,…,N；

t_Si和t_Ri分别为从第i个炮检对的炮点和检波点到各自射线末端的走时；T_{SRi_obs}为观测的第i个炮检对的反射波走时；(X_Si,Y_Si,Z_Si)和(X_Ri,Y_Ri,Z_Ri)分别为第i个炮检对的炮点和检波点的射线末端的坐标；V_i是第i个离散网格的速度值；V_Si和V_Ri分别为第i个炮检对的炮点和检波点射线末端所在单元的速度；P_SXi，P_SYi和P_SZi分别为第i个炮检对的炮点射线末端沿x,y和z方向的射线参数；P_RXi，P_RYi和P_RZi分别为第i个炮检对的检波点射线末端沿x,y和z方向的射线参数；L_Sij和L_Rij分别为第i个炮检对的炮点射线和检波点射线在第j个单元的射线段长度；L_SXij，L_SYij和L_SZij分别为第i个炮检对的炮点射线在第j个单元的射线段在x,y和z方向的投影长度；L_RXij，L_RYij和L_RZij分别为第i个炮检对的检波点点射线在第j个单元的射线段在x,y和z方向的投影长度；N为炮-检对总数；M为离散模型的网格节点总数；

S6，使用LSQR算法求解上述线性反演方程，得到离散速度模型更新量Δv,用该模型更新量修改当前离散速度模型，就得到修正后的模型，计算公式如下：

v^(k+1)＝v^(k)+Δv

其中，v^(k)是当前离散速度模型，即上次(第k次)迭代的离散速度模型，Δv是本次迭代反演得到的离散速度模型修正量；v^(k+1)是本次迭代得到的修正后的离散速度模型，再把修正后的离散速度模型作为新的初始模型，转到步骤四。

Images