CN111624647B - 一种变偏移距vsp射线追踪积分叠前时间偏移方法和装置 - Google Patents

Abstract

本发明公开了一种变偏移距VSP射线追踪积分叠前时间偏移方法和装置，包括以下步骤：S1.输入变偏移距VSP共检波点数据、2维速度地质模型和偏移参数；S2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；S3.射线追踪计算炮点到模型网格节点的走时，保存各个炮点的走时表；S4.射线追踪计算检波点到模型网格节点的走时，保存各个检波点的走时表；S5.计算2维地质模型的时深表，并基于步骤S3和步骤S4中得到的走时表计算共散射走时，实现变偏移距VSP共检波点积分叠前时间偏移成像；S6.利用步骤S5的积分叠前时间偏移重排成共成像道集，通过共成像道集叠加，得到变偏移VSP积分叠前时间偏移叠加成像。本发明具有走时计算精度高，成像效果好的优势。

Description

一种变偏移距VSP射线追踪积分叠前时间偏移方法和装置

技术领域

本发明涉及地震数据处理，特别是涉及一种变偏移距VSP射线追踪积分叠前时间偏移方法和装置。

背景技术

近年来井中地震技术发展迅速，变偏移距VSP(即Walkaway VSP)已实现工业化应用。通常，变偏移距VSP是利用成像数据来研究井旁构造、储层描述的。VSP成像方法有VSP-CDP成像、单程波相移叠前深度偏移、单程波有限差分叠前深度偏移、双程波逆时偏移、高斯束偏移、克希霍夫积分偏移。

王楠等研究了《表驱Kirchhoff叠前时间偏移角度域成像方法》，文献利用两点射线追踪计算走时和反射张角，实现地面地震数据Kirchhoff叠前时间偏移角度域成像。王珺等研究了《克希霍夫法VSP多波联合成像》，文献实现了偏移距VSP反射P波、反射S波、透射P波、透射S波克希霍夫法联合成像。

现有的变偏移距VSPKirchhoff叠前时间偏移采用双曲线公式、弯曲射线追踪计算双程时，精度不够、不能适应复杂构造及横向变速；没有保存走时表，重复计算，效率低；没有上下行波分离，计算效率低。

发明内容

本发明的目的在于克服现有技术的不足，提供一种变偏移距VSP射线追踪积分叠前时间偏移方法和装置，具有走时计算精度高，成像效果好的优势。

本发明的目的是通过以下技术方案来实现的：一种变偏移距VSP射线追踪积分叠前时间偏移方法，包括以下步骤：

S1.输入变偏移距VSP共检波点数据、2维速度地质模型和偏移参数；

所述步骤S1包括：

S101.输入变偏移距VSP共检波点波场数据{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；

S102.输入2维网格速度地质模型{ModVel}，读取网格节点数MM×MN、网格大小dg、井口模型坐标MWMX信息；

S103.输入偏移孔径aper、成像起始坐标xmig1、成像结束坐标xmig2、成像网格dxmig参数。

S2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；

所述步骤S2包括：用步骤S1中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在2维地质模型中的坐标：

炮点在2维地质模型中的坐标为：

MSZ_i＝SZ_i

其中，MSX_i是第i个炮点在2维地质模型中的坐标，SX_i、SY_i是第i个炮点的大地坐标，SZ_i是第i个炮点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；

检波点在2维地质模型中的坐标为：

MRZ_j＝RZ_j

其中，MRX_j是第j个检波点在2维地质模型中坐标，RX_j、RY_j是第j个检波点的大地坐标，RZ_j是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。

S3.在步骤S2的模型坐标系中，采用Runge-Kutta射线追踪计算炮点到模型网格节点的走时，保存各个炮点的走时表；

所述步骤S3包括：

S301.输入射线追踪参数，时间步长tstep、网格节点坐标{RayX,RayZ}、射线角度间隔RayDa；

S302.采用Runge-Kutta射线追踪单炮点单方向走时：

A1、设第i个炮点MSX_i、MSZ_i，第is条射线的初始方向为TestAngle_is，源点向量为：

A2、源点r_ti到r_ti+1的Runge-Kutta射线追踪为：

其中，r_ti是源点的射线参数数组，r_ti(1)、r_ti(2)是源点的x、z坐标，r_ti(3)、r_ti(4)是源点的射线方向参数；

是源点坐标对应的速度的平方；

是速度的x方向的偏导数在源点的值，

是速度的z方向的偏导数在源点的值。

rk_A＝r_ti+0.5k₁

其中，rk_A是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_A是射线参数数组，rk_A(1)、rk_A(2)是射线的x、z坐标，rk_A(3)、rk_A(4)是射线方向参数；

是rk_A对应的速度的平方；

是速度的x方向的偏导数在rk_A处的值，

是速度的z方向的偏导数在rk_A点的值。

rk_B＝r_ti+0.5k₂

其中，rk_B是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_B是射线参数数组，rk_B(1)、rk_B(2)是射线的x、z坐标，rk_B(3)、rk_B(4)是射线方向参数；

是rk_B对应的速度的平方；

是速度的x方向的偏导数在rk_B处的值，

是速度的z方向的偏导数在rk_B点的值。

rk_C＝r_ti+k₃

其中，rk_C是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_C是射线参数数组，rk_C(1)、rk_C(2)是射线的x、z坐标，rk_C(3)、rk_C(4)是射线方向参数；

是rk_C对应的速度的平方；

是速度的x方向的偏导数在rk_C处的值，

是速度的z方向的偏导数在rk_C点的值。

r_ti+1＝r_ti+k₁/6+k₂/3+k₃/3+k₄/6

r_ti+1(3)＝A

r_ti+1(4)＝B

t_ti+1＝t_ti+tstep

其中，r_ti+1是新的源点的射线参数数组对应时间为t_ti+1，r_ti是源点的射线参数数组t_ti，r_ti+1(3)、r_ti+1(4)是新的源点的射线方向参数，

是在rk_C点的速度。

A3、循环执行步骤A2直至模型边界；

S303.循环执行步骤S302实现Runge-Kutta射线追踪单炮点多方向走时：

给定第i个炮点的坐标为MSX_i、MSZ_i；射线的初始方向为TestAngle∈[0:RayDa:2π]；循环执行步骤S302，得到第i个炮点的射线角度间隔为RayDa的射线追踪数组RayS_i＝[r(1),r(2),t]。

S304.网格化步骤S303得到的射线追踪数组RayS_i＝[r(1),r(2),t]，得到：

RayTS_i＝[RayX,RayZ,RayT]，保存为第i个炮点的走时表；

S305.并行循环步骤S302～S304实现全部炮点射线追踪。

S4.在步骤S2的模型坐标系中，采用Runge-Kutta射线追踪计算检波点到模型网格节点的走时，保存各个检波点的走时表；

所述步骤S4包括：

S401.输入射线追踪参数，时间步长tstep、网格节点坐标{RayX,RayZ}、射线角度间隔RayDa；

S402.采用Runge-Kutta射线追踪单检波点单方向走时：

B1、设第j个检波点MRX_j、MRZ_j，第jr条射线的初始方向为TestAngle_jr，源点向量为

B2、源点r_ti到r_ti+1的Runge-Kutta射线追踪为：

其中，r_ti是源点的射线参数数组，r_ti(1)、r_ti(2)是源点的x、z坐标，r_ti(3)、r_ti(4)是源点的射线方向参数；

是源点坐标对应的速度的平方；

是速度的x方向的偏导数在源点的值，

是速度的z方向的偏导数在源点的值。

rk_A＝r_ti+0.5k₁

其中，rk_A是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_A是射线参数数组，rk_A(1)、rk_A(2)是射线的x、z坐标，rk_A(3)、rk_A(4)是射线方向参数；

是rk_A对应的速度的平方；

是速度的x方向的偏导数在rk_A处的值，

是速度的z方向的偏导数在rk_A点的值。

rk_B＝r_ti+0.5k₂

其中，rk_B是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_B是射线参数数组，rk_B(1)、rk_B(2)是射线的x、z坐标，rk_B(3)、rk_B(4)是射线方向参数；

是rk_B对应的速度的平方；

是速度的x方向的偏导数在rk_B处的值，

是速度的z方向的偏导数在rk_B点的值。

rk_C＝r_ti+k₃

其中，rk_C是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_C是射线参数数组，rk_C(1)、rk_C(2)是射线的x、z坐标，rk_C(3)、rk_C(4)是射线方向参数；

是rk_C对应的速度的平方；

是速度的x方向的偏导数在rk_C处的值，

是速度的z方向的偏导数在rk_C点的值。

r_ti+1＝r_ti+k₁/6+k₂/3+k₃/3+k₄/6

r_ti+1(3)＝A

r_ti+1(4)＝B

t_ti+1＝t_ti+tstep

其中，r_ti+1是新的源点的射线参数数组对应时间为t_ti+1，r_ti是源点的射线参数数组t_ti，r_ti+1(3)、r_ti+1(4)是新的源点的射线方向参数，

是在rk_C点的速度。

B3、循环执行步骤B2直至模型边界。

S403.循环执行步骤S402实现Runge-Kutta射线追踪单检波点多方向走时

给定第j个检波点的坐标为MRX_j、MRZ_j；射线的初始方向为TestAngle∈[0:RayDa:2π]；循环执行步骤S402，得到第j个检波点的射线角度间隔为RayDa的射线追踪数组RayR_j＝[r(1),r(2),t]；

S404.网格化步骤S403得到的射线追踪数组RayR_j＝[r(1),r(2),t]为：

RayTR_j＝[RayX,RayZ,RayT]，保存为j个检波点的走时表；

S405.并行循环步骤S402～S403，实现全部检波点射线追踪。

S5.在步骤S2的模型坐标系中，计算2维地质模型的时深表，并基于步骤S3和步骤S4中得到的走时表计算共散射走时，实现变偏移距VSP共检波点积分叠前时间偏移成像；

所述步骤S5包括：

S501.计算步骤S1中2维地质模型的时深表

ModT_k,l＝ModT_k-1,l+2·dg/ModVel_k-1,l

ModH_k,l＝ModH_k-1,l+dg

其中，ModT是模型双程垂直时，ModH是模型深度，k∈[1,MM]是模型行数，l∈[1,MN]是模型列数；

S502.利用步骤S1的偏移参数定义成像坐标，选择步骤S1中输入数据中第j个检波点的波场数据VSPData_j

xmig_ktr＝xmig1+(ktr-1)·dxmig ktr＝1,2,…,(xmig2-xmig1)/dxmig+1

其中，xmig_ktr是第ktr道的成像坐标，xmig1是成像起始坐标，xmig2是成像结束坐标；

MigData_ktr,j＝{0}

其中，MigData_ktr,j第j个检波点第ktr道的成像数据，先置为零；

S503.利用步骤S1的偏移孔径、步骤S2的炮点模型坐标，结合成像坐标确定成像输入地震道范围

xmig_ktr-aper≤MSX_kaper,j≤xmig_ktr+aper

其中，xmig_ktr是第ktr道的成像坐标，aper是偏移孔径，MSX_kaper,j是第j个检波点第ktr个成像道孔径内第kaper个炮点的模型坐标；

S504.选择步骤S3第kaper个炮点xmig_ktr处的走时加上步骤S4第个j检波点xmig_ktr处的走时，即为xmig_ktr处kaper、j炮检对的共散射走时；

其中，TSR是xmig_ktr处kaper、j炮检对的共散射走时，

是第kaper个炮点xmig_ktr处的走时，

第j个检波点xmig_ktr处的走时；

S505.选择步骤S503中第j个检波点第kaper个炮点的波场数据VSPData_kaper,j，把TSR对应的振幅映射为步骤S501中时深表对应的双程垂直时，与xmig_ktr处的成像道MigData_ktr,j积分求和；

S506.循环执行步骤S503～S505，实现第j个检波点第ktr个成像道孔径内全部地震道积分叠前时间偏移成像，得到MigData_ktr,j；

S507.循环执行步骤S502～S506，实现第j个检波点全部成像道的积分叠前时间偏移成像得到MigData_j；

S508.并行循环步骤S502～S507，实现全部检波点全部成像道的积分叠前时间偏移成像,得到MigData。

S6.利用步骤S5的积分叠前时间偏移重排成共成像道集，通过共成像道集叠加，得到变偏移VSP积分叠前时间偏移叠加成像。

一种变偏移距VSP射线追踪积分叠前时间偏移装置，包括：

数据输入单元，用于输入变偏移距VSP共检波点数据、2维速度地质模型和偏移参数；

模型坐标建立单元，用于根据输入数据，建立炮点、检波点在2维地质模型中的坐标；

炮点走时表追踪计算单元，用于采用Runge-Kutta射线追踪计算炮点到模型网格节点的走时，保存各个炮点的走时表；

检波点走时表追踪计算单元，用于采用Runge-Kutta射线追踪计算检波点到模型网格节点的走时，保存各个检波点的走时表；

叠前时间偏移成像单元，用于计算2维地质模型的时深表，并基于得到的走时表计算共散射走时，实现变偏移距VSP共检波点积分叠前时间偏移成像；

变偏移VSP叠加成像单元，用于将积分叠前时间偏移重排成共成像道集，通过共成像道集叠加，得到变偏移VSP积分叠前时间偏移叠加成像。

本发明的有益效果是：本发明利用Runge-Kutta射线追踪求解程函方程计算走时；先计算炮点到地质模型各个网格节点的走时、存为走时表，再计算检波点到地质模型各个网格节点的走时、存为走时表；炮点的走时加上检波点走时，即为共散射点走时；将共散射点走时对应的地震数据映射为双程垂直时，偏移孔径内积分求和，实现共检波点波场积分叠前时间偏移成像；多个检波点的共散射成像叠加，实现变偏移距VSP积分叠前时间偏移成像叠加，具有走时计算精度高，成像效果好的优势。

附图说明

图1为本发明的方法流程图；

图2为实施例中输入变偏移距VSP共检波点波场示意图；

图3为实施例中输入2维速度地质模型示意图；

图4为实施例中炮点、检波点的模型坐标示意图；

图5为实施例中炮点到模型网格节点走时，等值线显示示意图；

图6为实施例中检波点到模型网格节点走时，等值线显示示意图；

图7为实施例中共散射的网格节点走时，等值线显示示意图；

图8为实施例中变偏移距VSP共检波积分叠前时间偏移成像剖面示意图；

图9为实施例中变偏移距VSP共检波积分叠前时间偏移成像共成像道集剖面示意图；

图10为变偏移距VSP共检波积分叠前时间偏移成像叠加剖面实施例中示意图；

图11为实施例中变偏移距VSP共检波积分叠前时间偏移成像叠加后去噪剖面示意图；

图12为本发明的装置原理示意图。

具体实施方式

下面结合附图进一步详细描述本发明的技术方案，但本发明的保护范围不局限于以下所述。

进行变偏移距VSP井中地震数据采集时，需要在地面布设一条由多个炮点所组成的炮线，在井中的各个检波点设置检波器，在炮点人工激发地震数据的条件下，实现变偏移距VSP井中数据采集，本申请结合VSP井中数据采集时得到的数据和相关参数，实现变偏移距VSP积分叠前时间偏移成像叠加，具体地：

如图1所示，一种变偏移距VSP射线追踪积分叠前时间偏移方法，包括以下步骤：

S1.输入变偏移距VSP共检波点数据、2维速度地质模型和偏移参数；

所述步骤S1包括：

S101.输入变偏移距VSP共检波点波场数据{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；

S102.输入2维网格速度地质模型{ModVel}，读取网格节点数MM×MN、网格大小dg、井口模型坐标MWMX信息；

S103.输入偏移孔径aper、成像起始坐标xmig1、成像结束坐标xmig2、成像网格dxmig参数。

在本申请的实施例中，输入的变偏移距VSP共检波点波场如图2所示，图2中，横坐标为道号，纵坐标为时间(单位：毫秒)。输入的2维速度地质模型如图3所示，图3中，横坐标为长度(单位：米)，纵坐标为深度(单位：米)。

S2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；

所述步骤S2包括：用步骤S1中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在2维地质模型中的坐标：

炮点在2维地质模型中的坐标为：

MSZ_i＝SZ_i

其中，MSX_i是第i个炮点在2维地质模型中的坐标，SX_i、SY_i是第i个炮点的大地坐标，SZ_i是第i个炮点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；

检波点在2维地质模型中的坐标为：

MRZ_j＝RZ_j

其中，MRX_j是第j个检波点在2维地质模型中坐标，RX_j、RY_j是第j个检波点的大地坐标，RZ_j是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。

在上述实施例中，建立的炮点、检波点模型坐标如图4所示，图4中横坐标为长度(单位：米)；纵坐标为深度(单位：米)；

S3.在步骤S2的模型坐标系中，采用Runge-Kutta射线追踪计算炮点到模型网格节点的走时，保存各个炮点的走时表；

所述步骤S3包括：

S301.输入射线追踪参数，时间步长tstep、网格节点坐标{RayX,RayZ}、射线角度间隔RayDa；

S302.采用Runge-Kutta射线追踪单炮点单方向走时：

A1、设第i个炮点MSX_i、MSZ_i，第is条射线的初始方向为TestAngle_is，源点向量为：

A2、源点r_ti到r_ti+1的Runge-Kutta射线追踪为：

其中，r_ti是源点的射线参数数组，r_ti(1)、r_ti(2)是源点的x、z坐标，r_ti(3)、r_ti(4)是源点的射线方向参数；

是源点坐标对应的速度的平方；

是速度的x方向的偏导数在源点的值，

是速度的z方向的偏导数在源点的值。

rk_A＝r_ti+0.5k₁

其中，rk_A是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_A是射线参数数组，rk_A(1)、rk_A(2)是射线的x、z坐标，rk_A(3)、rk_A(4)是射线方向参数；

是rk_A对应的速度的平方；

是速度的x方向的偏导数在rk_A处的值，

是速度的z方向的偏导数在rk_A点的值。

rk_B＝r_ti+0.5k₂

其中，rk_B是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_B是射线参数数组，rk_B(1)、rk_B(2)是射线的x、z坐标，rk_B(3)、rk_B(4)是射线方向参数；

是rk_B对应的速度的平方；

是速度的x方向的偏导数在rk_B处的值，

是速度的z方向的偏导数在rk_B点的值。

rk_C＝r_ti+k₃

其中，rk_C是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_C是射线参数数组，rk_C(1)、rk_C(2)是射线的x、z坐标，rk_C(3)、rk_C(4)是射线方向参数；

是rk_C对应的速度的平方；

是速度的x方向的偏导数在rk_C处的值，

是速度的z方向的偏导数在rk_C点的值。

r_ti+1＝r_ti+k₁/6+k₂/3+k₃/3+k₄/6

r_ti+1(3)＝A

r_ti+1(4)＝B

t_ti+1＝t_ti+tstep

其中，r_ti+1是新的源点的射线参数数组对应时间为t_ti+1，r_ti是源点的射线参数数组t_ti，r_ti+1(3)、r_ti+1(4)是新的源点的射线方向参数，

是在rk_C点的速度。

A3、循环执行步骤A2直至模型边界；

S303.循环执行步骤S302实现Runge-Kutta射线追踪单炮点多方向走时：

给定第i个炮点的坐标为MSX_i、MSZ_i；射线的初始方向为TestAngle∈[0:RayDa:2π]；循环执行步骤S302，得到第i个炮点的射线角度间隔为RayDa的射线追踪数组RayS_i＝[r(1),r(2),t]。

S304.网格化步骤S303得到的射线追踪数组RayS_i＝[r(1),r(2),t]，得到：

RayTS_i＝[RayX,RayZ,RayT]，保存为第i个炮点的走时表；

S305.并行循环步骤S302～S304实现全部炮点射线追踪。

在上述实施例中，炮点到模型网格节点走时、等值线显示如图5所示，图5中横坐标为长度(单位：米)；纵坐标为深度(单位：米)；

S4.在步骤S2的模型坐标系中，采用Runge-Kutta射线追踪计算检波点到模型网格节点的走时，保存各个检波点的走时表；

所述步骤S4包括：

S401.输入射线追踪参数，时间步长tstep、网格节点坐标{RayX,RayZ}、射线角度间隔RayDa；

S402.采用Runge-Kutta射线追踪单检波点单方向走时：

B1、设第j个检波点MRX_j、MRZ_j，第jr条射线的初始方向为TestAngle_jr，源点向量为

B2、源点r_ti到r_ti+1的Runge-Kutta射线追踪为：

其中，r_ti是源点的射线参数数组，r_ti(1)、r_ti(2)是源点的x、z坐标，r_ti(3)、r_ti(4)是源点的射线方向参数；

是源点坐标对应的速度的平方；

是速度的x方向的偏导数在源点的值，

是速度的z方向的偏导数在源点的值。

rk_A＝r_ti+0.5k₁

其中，rk_A是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_A是射线参数数组，rk_A(1)、rk_A(2)是射线的x、z坐标，rk_A(3)、rk_A(4)是射线方向参数；

是rk_A对应的速度的平方；

是速度的x方向的偏导数在rk_A处的值，

是速度的z方向的偏导数在rk_A点的值。

rk_B＝r_ti+0.5k₂

其中，rk_B是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_B是射线参数数组，rk_B(1)、rk_B(2)是射线的x、z坐标，rk_B(3)、rk_B(4)是射线方向参数；

是rk_B对应的速度的平方；

是速度的x方向的偏导数在rk_B处的值，

是速度的z方向的偏导数在rk_B点的值。

rk_C＝r_ti+k₃

其中，rk_C是射线参数数组，r_ti是源点的射线参数数组。

其中，rk_C是射线参数数组，rk_C(1)、rk_C(2)是射线的x、z坐标，rk_C(3)、rk_C(4)是射线方向参数；

是rk_C对应的速度的平方；

是速度的x方向的偏导数在rk_C处的值，

是速度的z方向的偏导数在rk_C点的值。

r_ti+1＝r_ti+k₁/6+k₂/3+k₃/3+k₄/6

r_ti+1(3)＝A

r_ti+1(4)＝B

t_ti+1＝t_ti+tstep

其中，r_ti+1是新的源点的射线参数数组对应时间为t_ti+1，r_ti是源点的射线参数数组t_ti，r_ti+1(3)、r_ti+1(4)是新的源点的射线方向参数，

是在rk_C点的速度。

B3、循环执行步骤B2直至模型边界。

S403.循环执行步骤S402实现Runge-Kutta射线追踪单检波点多方向走时

给定第j个检波点的坐标为MRX_j、MRZ_j；射线的初始方向为TestAngle∈[0:RayDa:2π]；循环执行步骤S402，得到第j个检波点的射线角度间隔为RayDa的射线追踪数组RayR_j＝[r(1),r(2),t]；

S404.网格化步骤S403得到的射线追踪数组RayR_j＝[r(1),r(2),t]为：

RayTR_j＝[RayX,RayZ,RayT]，保存为j个检波点的走时表；

S405.并行循环步骤S402～S403，实现全部检波点射线追踪。

在上述的实施例中，检波点到模型网格节点走时、等值线显示如图6所示，图6中横坐标为长度(单位：米)；纵坐标为深度(单位：米)。

S5.在步骤S2的模型坐标系中，计算2维地质模型的时深表，并基于步骤S3和步骤S4中得到的走时表计算共散射走时，实现变偏移距VSP共检波点积分叠前时间偏移成像；

所述步骤S5包括：

S501.计算步骤S1中2维地质模型的时深表

ModT_k,l＝ModT_k-1,l+2·dg/ModVel_k-1,l

ModH_k,l＝ModH_k-1,l+dg

其中，ModT是模型双程垂直时，ModH是模型深度，k∈[1,MM]是模型行数，l∈[1,MN]是模型列数；

S502.利用步骤S1的偏移参数定义成像坐标，选择步骤S1中输入数据中第j个检波点的波场数据VSPData_j

xmig_ktr＝xmig1+(ktr-1)·dxmig ktr＝1,2,…,(xmig2-xmig1)/dxmig+1

其中，xmig_ktr是第ktr道的成像坐标，xmig1是成像起始坐标，xmig2是成像结束坐标；

MigData_ktr,j＝{0}

其中，MigData_ktr,j第j个检波点第ktr道的成像数据，先置为零；

S503.利用步骤S1的偏移孔径、步骤S2的炮点模型坐标，结合成像坐标确定成像输入地震道范围

xmig_ktr-aper≤MSX_kaper,j≤xmig_ktr+aper

其中，xmig_ktr是第ktr道的成像坐标，aper是偏移孔径，MSX_kaper,j是第j个检波点第ktr个成像道孔径内第kaper个炮点的模型坐标；

S504.选择步骤S3第kaper个炮点xmig_ktr处的走时加上步骤S4第个j检波点xmig_ktr处的走时，即为xmig_ktr处kaper、j炮检对的共散射走时；

其中，TSR是xmig_ktr处kaper、j炮检对的共散射走时，

是第kaper个炮点xmig_ktr处的走时，

第j个检波点xmig_ktr处的走时；

在上述实施例中，共散射的网格节点走时、等值线显示如图7所示，图7中横坐标为长度(单位：米)；纵坐标为深度(单位：米)。

S505.选择步骤S503中第j个检波点第kaper个炮点的波场数据VSPData_kaper,j，把TSR对应的振幅映射为步骤S501中时深表对应的双程垂直时，与xmig_ktr处的成像道MigData_ktr,j积分求和；

S506.循环执行步骤S503～S505，实现第j个检波点第ktr个成像道孔径内全部地震道积分叠前时间偏移成像，得到MigData_ktr,j；

S507.循环执行步骤S502～S506，实现第j个检波点全部成像道的积分叠前时间偏移成像得到MigData_j；

S508.并行循环步骤S502～S507，实现全部检波点全部成像道的积分叠前时间偏移成像,得到MigData。

在上述实施例中，变偏移距VSP共检波积分叠前时间偏移成像剖面如图8所示，图8中横坐标为道号；纵坐标为时间(单位：毫秒)；

S6.利用步骤S5的积分叠前时间偏移重排成共成像道集，通过共成像道集叠加，得到变偏移VSP积分叠前时间偏移叠加成像。

在上述实施例中，变偏移距VSP共检波积分叠前时间偏移成像共成像道集剖面如图9所示。变偏移距VSP共检波积分叠前时间偏移成像叠加剖面如图10所示。变偏移距VSP共检波积分叠前时间偏移成像叠加后去噪剖面如图11所示；在图9～11中，横坐标为道号；纵坐标为时间(单位：毫秒)。

如图12所示，一种变偏移距VSP射线追踪积分叠前时间偏移装置，包括：

数据输入单元，用于输入变偏移距VSP共检波点数据、2维速度地质模型和偏移参数；

模型坐标建立单元，用于根据输入数据，建立炮点、检波点在2维地质模型中的坐标；

炮点走时表追踪计算单元，用于采用Runge-Kutta射线追踪计算炮点到模型网格节点的走时，保存各个炮点的走时表；

检波点走时表追踪计算单元，用于采用Runge-Kutta射线追踪计算检波点到模型网格节点的走时，保存各个检波点的走时表；

叠前时间偏移成像单元，用于计算2维地质模型的时深表，并基于得到的走时表计算共散射走时，实现变偏移距VSP共检波点积分叠前时间偏移成像；

变偏移VSP叠加成像单元，用于将积分叠前时间偏移重排成共成像道集，通过共成像道集叠加，得到变偏移VSP积分叠前时间偏移叠加成像。

综上，本发明利用Runge-Kutta射线追踪求解程函方程计算走时；先计算炮点到地质模型各个网格节点的走时、存为走时表，再计算检波点到地质模型各个网格节点的走时、存为走时表；炮点的走时加上检波点走时，即为共散射点走时；将共散射点走时对应的地震数据映射为双程垂直时，偏移孔径内积分求和，实现共检波点波场积分叠前时间偏移成像；多个检波点的共散射成像叠加，实现变偏移距VSP积分叠前时间偏移成像叠加，具有走时计算精度高，成像效果好的优势。

显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其他不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引申出的显而易见的变化或变动仍处于本发明创造的保护范围之中。

Claims (5)

1.一种变偏移距VSP射线追踪积分叠前时间偏移方法，其特征在于：包括以下步骤：

步骤S1.输入变偏移距VSP共检波点数据、2维速度地质模型和偏移参数；

步骤S2.根据输入数据，建立炮点、检波点在2维地质模型中的坐标；

步骤S3.在步骤S2的模型坐标系中，采用Runge-Kutta射线追踪计算炮点到模型网格节点的走时，保存各个炮点的走时表；

所述步骤S3包括：

S301.输入射线追踪参数，时间步长tstep、网格节点坐标{RayX,RayZ}、射线角度间隔RayDa；

S302.采用Runge-Kutta射线追踪单炮点单方向走时：

A1、设第i个炮点MSX_i、MSZ_i，第is条射线的初始方向为TestAngle_is，源点向量为：

是第i个炮点MSX_i、MSZ_i处的速度；

A2、源点r_ti到r_ti+1的Runge-Kutta射线追踪为：

其中，r_ti是源点的射线参数数组，r_ti(1)、r_ti(2)是源点的x、z坐标，r_ti(3)、r_ti(4)是源点的射线方向参数；

是源点坐标对应的速度的平方；

是速度的x方向的偏导数在源点的值，

是速度的z方向的偏导数在源点的值；

rk_A＝r_ti+0.5k₁

其中，rk_A是射线参数数组，r_ti是源点的射线参数数组；

其中，rk_A是射线参数数组，rk_A(1)、rk_A(2)是射线的x、z坐标，rk_A(3)、rk_A(4)是射线方向参数；

是rk_A对应的速度的平方；

是速度的x方向的偏导数在rk_A处的值，

是速度的z方向的偏导数在rk_A点的值；

rk_B＝r_ti+0.5k₂

其中，rk_B是射线参数数组，r_ti是源点的射线参数数组；

其中，rk_B是射线参数数组，rk_B(1)、rk_B(2)是射线的x、z坐标，rk_B(3)、rk_B(4)是射线方向参数；

是rk_B对应的速度的平方；

是速度的x方向的偏导数在rk_B处的值，

是速度的z方向的偏导数在rk_B点的值；

rk_C＝r_ti+k₃

其中，rk_C是射线参数数组，r_ti是源点的射线参数数组；

其中，rk_C是射线参数数组，rk_C(1)、rk_C(2)是射线的x、z坐标，rk_C(3)、rk_C(4)是射线方向参数；

是rk_C对应的速度的平方；

是速度的x方向的偏导数在rk_C处的值，

是速度的z方向的偏导数在rk_C点的值；

r_ti+1＝r_ti+k₁/6+k₂/3+k₃/3+k₄/6

r_ti+1(3)＝A

r_ti+1(4)＝B

t_ti+1＝t_ti+tstep

其中，r_ti+1是新的源点的射线参数数组，对应时间为t_ti+1，r_ti是源点的射线参数数组，对应时间为t_ti；r_ti+1(3)、r_ti+1(4)是新的源点的射线方向参数，

是在rk_C点的速度；

A3、循环执行步骤A2直至模型边界；

S303.循环执行步骤S302实现Runge-Kutta射线追踪单炮点多方向走时：

给定第i个炮点的坐标为MSX_i、MSZ_i；射线的初始方向为TestAngle∈[0:RayDa:2π]；循环执行步骤S302，得到第i个炮点的射线角度间隔为RayDa的射线追踪数组RayS_i＝[r(1),r(2),t]；r(1)是射线参数数组的元素1，r(2)是射线参数数组的元素2，t是走时；

S304.网格化步骤S303得到的射线追踪数组RayS_i＝[r(1),r(2),t]，得到：

RayTS_i＝[RayX,RayZ,RayT]，保存为第i个炮点的走时表；

RayT是网格坐标RayX、RayZ对应的网格走时；

S305.并行循环步骤S302～S304实现全部炮点射线追踪；

步骤S4.在步骤S2的模型坐标系中，采用Runge-Kutta射线追踪计算检波点到模型网格节点的走时，保存各个检波点的走时表；

所述步骤S4包括：

S401.输入射线追踪参数，时间步长tstep、网格节点坐标{RayX,RayZ}、射线角度间隔RayDa；

S402.采用Runge-Kutta射线追踪单检波点单方向走时：

B1、设第j个检波点MRX_j、MRZ_j，第jr条射线的初始方向为TestAngle_jr，源点向量为

是第j个检波点MRX_j、MRZ_j处的速度；

B2、源点r_ti到r_ti+1的Runge-Kutta射线追踪为：

其中，r_ti是源点的射线参数数组，r_ti(1)、r_ti(2)是源点的x、z坐标，r_ti(3)、r_ti(4)是源点的射线方向参数；

是源点坐标对应的速度的平方；

是速度的x方向的偏导数在源点的值，

是速度的z方向的偏导数在源点的值；

rk_A＝r_ti+0.5k₁

其中，rk_A是射线参数数组，r_ti是源点的射线参数数组；

其中，rk_A是射线参数数组，rk_A(1)、rk_A(2)是射线的x、z坐标，rk_A(3)、rk_A(4)是射线方向参数；

是rk_A对应的速度的平方；

是速度的x方向的偏导数在rk_A处的值，

是速度的z方向的偏导数在rk_A点的值；

rk_B＝r_ti+0.5k₂

其中，rk_B是射线参数数组，r_ti是源点的射线参数数组；

其中，rk_B是射线参数数组，rk_B(1)、rk_B(2)是射线的x、z坐标，rk_B(3)、rk_B(4)是射线方向参数；

是rk_B对应的速度的平方；

是速度的x方向的偏导数在rk_B处的值，

是速度的z方向的偏导数在rk_B点的值；

rk_C＝r_ti+k₃

其中，rk_C是射线参数数组，r_ti是源点的射线参数数组；

其中，rk_C是射线参数数组，rk_C(1)、rk_C(2)是射线的x、z坐标，rk_C(3)、rk_C(4)是射线方向参数；

是rk_C对应的速度的平方；

是速度的x方向的偏导数在rk_C处的值，

是速度的z方向的偏导数在rk_C点的值；

r_ti+1＝r_ti+k₁/6+k₂/3+k₃/3+k₄/6

r_ti+1(3)＝A

r_ti+1(4)＝B

t_ti+1＝t_ti+tstep

其中，r_ti+1是新的源点的射线参数数组，对应时间为t_ti+1，r_ti是源点的射线参数数组，对应时间为t_ti；r_ti+1(3)、r_ti+1(4)是新的源点的射线方向参数，

是在rk_C点的速度；

B3、循环执行步骤B2直至模型边界；

S403.循环执行步骤S402实现Runge-Kutta射线追踪单检波点多方向走时

给定第j个检波点的坐标为MRX_j、MRZ_j；射线的初始方向为TestAngle∈[0:RayDa:2π]；循环执行步骤S402，得到第j个检波点的射线角度间隔为RayDa的射线追踪数组RayR_j＝[r(1),r(2),t]；

S404.网格化步骤S403得到的射线追踪数组RayR_j＝[r(1),r(2),t]为：

RayTR_j＝[RayX,RayZ,RayT]，保存为j个检波点的走时表；

S405.并行循环步骤S402～S403，实现全部检波点射线追踪；

步骤S5.在步骤S2的模型坐标系中，计算2维地质模型的时深表，并基于步骤S3和步骤S4中得到的走时表计算共散射走时，实现变偏移距VSP共检波点积分叠前时间偏移成像；

步骤S6.利用步骤S5的积分叠前时间偏移重排成共成像道集，通过共成像道集叠加，得到变偏移VSP积分叠前时间偏移叠加成像。

2.根据权利要求1所述的一种变偏移距VSP射线追踪积分叠前时间偏移方法，其特征在于：所述步骤S1包括：

S101.输入变偏移距VSP共检波点波场数据{VSPData}，读取炮点坐标{SX,SY}、炮点深度SZ、检波点坐标{RX,RY}、检波点深度RZ、井口大地坐标{WMX,WMY}；

S102.输入2维网格速度地质模型{ModVel}，读取网格节点数MM×MN、网格大小dg、井口模型坐标MWMX信息；

S103.输入偏移孔径aper、成像起始坐标xmig1、成像结束坐标xmig2、成像网格dxmig参数。

3.根据权利要求2所述的一种变偏移距VSP射线追踪积分叠前时间偏移方法，其特征在于：所述步骤S2包括：用步骤S1中的炮点坐标、炮点深度、检波点坐标、检波点深度、井口大地坐标、井口模型坐标，建立炮点、检波点在2维地质模型中的坐标：

炮点在2维地质模型中的坐标为：

MSZ_i＝SZ_i

其中，MSX_i是第i个炮点在2维地质模型中的坐标，炮点在2维地质模型中的坐标也称为炮点模型坐标，SX_i、SY_i是第i个炮点的大地坐标，SZ_i是第i个炮点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数；

检波点在2维地质模型中的坐标为：

MRZ_j＝RZ_j

其中，MRX_j是第j个检波点在2维地质模型中坐标，RX_j、RY_j是第j个检波点的大地坐标，RZ_j是第j个检波点的深度，WMX、WMY是井口的大地坐标，MWMX是井口模型坐标，sign是符号函数。

4.根据权利要求3所述的一种变偏移距VSP射线追踪积分叠前时间偏移方法，其特征在于：所述步骤S5包括：

S501.计算步骤S1中2维地质模型的时深表

ModT_k,l＝ModT_k-1,l+2·dg/ModVel_k-1,l

ModH_k,l＝ModH_k-1,l+dg

其中，ModT是模型双程垂直时，ModH是模型深度，k∈[1,MM]是模型行数，l∈[1,MN]是模型列数；

S502.利用步骤S1的偏移参数定义成像坐标，选择步骤S1中输入数据中第j个检波点的波场数据VSPData_j

xmig_ktr＝xmig1+(ktr-1)·dxmig ktr＝1,2,…,(xmig2-xmig1)/dxmig+1

其中，xmig_ktr是第ktr道的成像坐标，xmig1是成像起始坐标，xmig2是成像结束坐标；

MigData_ktr,j＝{0}

其中，MigData_ktr,j第j个检波点第ktr道的成像数据，先置为零；

S503.利用步骤S1的偏移孔径、步骤S2的炮点模型坐标，结合成像坐标确定成像输入地震道范围

xmig_ktr-aper≤MSX_kaper,j≤xmig_ktr+aper

其中，xmig_ktr是第ktr道的成像坐标，aper是偏移孔径，MSX_kaper,j是第j个检波点第ktr个成像道孔径内第kaper个炮点的模型坐标；

S504.选择步骤S3第kaper个炮点xmig_ktr处的走时加上步骤S4第个j检波点xmig_ktr处的走时，即为xmig_ktr处kaper、j炮检对的共散射走时；

其中，TSR是xmig_ktr处kaper、j炮检对的共散射走时，

是第kaper个炮点xmig_ktr处的走时，

第j个检波点xmig_ktr处的走时；

S505.选择步骤S503中第j个检波点第kaper个炮点的波场数据VSPData_kaper,j，把TSR对应的振幅映射为步骤S501中时深表对应的双程垂直时，与xmig_ktr处的成像道MigData_ktr,j积分求和；

S506.循环执行步骤S503～S505，实现第j个检波点第ktr个成像道孔径内全部地震道积分叠前时间偏移成像，得到MigData_ktr,j；

S507.循环执行步骤S502～S506，实现第j个检波点全部成像道的积分叠前时间偏移成像得到MigData_j；

S508.并行循环步骤S502～S507，实现全部检波点全部成像道的积分叠前时间偏移成像,得到MigData。

5.一种变偏移距VSP射线追踪积分叠前时间偏移装置，采用权利要求1～4中任意一项所述的方法，其特征在于：包括：

数据输入单元，用于输入变偏移距VSP共检波点数据、2维速度地质模型和偏移参数；

模型坐标建立单元，用于根据输入数据，建立炮点、检波点在2维地质模型中的坐标；

炮点走时表追踪计算单元，用于采用Runge-Kutta射线追踪计算炮点到模型网格节点的走时，保存各个炮点的走时表；

检波点走时表追踪计算单元，用于采用Runge-Kutta射线追踪计算检波点到模型网格节点的走时，保存各个检波点的走时表；

叠前时间偏移成像单元，用于计算2维地质模型的时深表，并基于得到的走时表计算共散射走时，实现变偏移距VSP共检波点积分叠前时间偏移成像；

变偏移VSP叠加成像单元，用于将积分叠前时间偏移重排成共成像道集，通过共成像道集叠加，得到变偏移VSP积分叠前时间偏移叠加成像。

Images