横向各向同性介质转换波角道集抽取方法和装置与流程

本公开涉及地震勘探数据处理领域，更具体地，涉及横向各向同性介质转换波角道集抽取方法和装置。

背景技术：

在叠前转换波共转换点道集上，转换波振幅随偏移距变化而变化(AVO)。叠前角度域AVO反演可以获得储层的纵波速度、横波速度、密度等，是储层弹性参数反演和储层预测、流体预测的重要技术。偏移距到角度的转换是叠前角度域AVO反演的关键技术。每个地震记录道都具有一个固定的炮检距，道上不同时刻采样点具有不同反射角度。为了获得反射波振幅随入射角变化的信息，必须把时间-偏移距域的固定偏移距共反射点道集转换为时间-角度域的固定反射角度的道集。

目前转换波入射角计算的方法主要有：1)直射线法；2)弯曲射线法；3)纵波等效纵波偏移距法。直射线法是假设反射界面以上为单一均匀介质，从地面到反射层位的速度为常数，忽略速度随深度变化，射线是反射点与出射点之间的直线，该方法计算的角度偏小，缺失大角度，不利于AVO反演,计算速度快，但是精度差。弯曲射线法是假设反射界面以上是多层介质，模型的速度随深度关系是阶跃式的，在已知纵波层速度、横波层速度和厚度情况下，用射线追踪求取转换波入射角。该模型需要人为设定速度界面的层数。由于折射使出射点偏移量随反射深度逆转，造成角度形态畸变。弯曲射线法精度较高，但是计算速度慢。纵波等效纵波偏移距法(参见参考文件1：中国专利申请号201210272650.2发明名称为“一种转换波角道集抽取方法”)是通过转化将转换 波炮检距、旅行时和速度映射为相应的纵波炮检距、旅行时和纵波速度，用纵波的射线参数法求取转换波入射角。该算法具有精度高的优点，速度介于直射线和弯曲射线法之间，但是在转化过程中，稳定性差，容易引人误差，且计算精度受纵波等效偏移距和纵波各向异常参数转换精度和映射后纵波旅行时计算精度约束

技术实现要素：

本公开针对横向各向同性介质提出了一种高精度的横向各向同性介质转换波角道集抽取方法。

根据一方面，提出了一种横向各向同性介质转换波角道集抽取方法，包括：基于叠前转换波的共转换点道集，通过速度分析得到转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff和转换波各向异性参数χ_eff；在第k个采样点处，针对共转换点道集中的各个地震道进行以下处理：

-根据转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff计算纵波叠加速度V_p2；

-根据转换波双程旅行时t_c和垂直速度比γ₀计算纵波双程旅行时t_p0；

-根据纵波双程旅行时t_p0和纵波叠加速度V_p2计算纵波层速度V_p,k；

-根据转换波双程旅行时t_c、转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff和转换波各向异性参数χ_eff计算射线参数p；

-根据纵波层速度V_p,k、射线参数p计算转换波入射角θ_k，其中，k为正整数，表示采样点编号且对应于层编号，转换波入射角θ_k为转换波在第k层处的入射角；

取k为1到n，并针对每个k值执行上述处理，以得到在每个采样点处针对各个地震道的转换波入射角，将转换波入射角相同的地震道归为同一入射角道集，得到转换波角道集，其中n为总采样点数。

另一方面提出了一种横向各向同性介质转换波角道集抽取装置，包括：用于基于叠前转换波的共转换点道集，通过速度分析得到转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff和转换波各向异性参数χ_eff的部件；用于在第k个采样点处，针对共转换点道集中的各个地震道进行以下处理的部件：

-根据转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff计算纵波叠加速度V_p2；

-根据转换波双程旅行时t_c和垂直速度比γ₀计算纵波双程旅行时t_p0；

-根据纵波双程旅行时t_p0和纵波叠加速度V_p2计算纵波层速度V_p,k；

-根据转换波双程旅行时t_c、转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff和转换波各向异性参数χ_eff计算射线参数p；

-根据纵波层速度V_p,k、射线参数p计算转换波入射角θ_k，其中，k为正整数，表示采样点编号且对应于层编号，转换波入射角θ_k为转换波在第k层处的入射角；

用于取k为1到n，并针对每个k值执行上述处理，以得到在每个采样点处针对各个地震道的转换波入射角，将转换波入射角相同的地震道归为同一入射角道集，得到转换波角道集的部件，其中n为总采样点数。

本公开提出了稳定的高精度的转换波入射角道集抽取方法，该方法能正确描述振幅随入射角变化关系，提高了叠前AVO反演地层弹性参数精度和地震资料储层含气性预测精度。

附图说明

通过结合附图对本公开示例性实施方式进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显，其中，在本公开示例性实施方式中，相同的参考标号通常代表相同部件。

图1示出了根据本发明一个实施例的横向各向同性介质转换波角道集抽取 方法的流程图。

图2示出了典型的层状介质转换波射线路径的示意图。

图3示出了一个示例性的层状介质模型的纵波速度、横波速度和密度曲线示意图。

图4示出了正演模拟转换波道集经过动校正处理后的道集的示意图。

图5和图6分别示出了层状介质模型的第5层和第9层的本发明实施例计算的入射角与射线追踪的真正入射角的比较示意图。

图7示出了转换波入射角道集的一个示例的示意图。

具体实施方式

下面将参照附图更详细地描述本公开的优选实施方式。虽然附图中显示了本公开的优选实施方式，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

图1示出了根据本发明一个实施例的横向各向同性介质转换波角道集抽取方法的流程图，该方法包括：

步骤101，基于叠前转换波的共转换点道集，即CCP道集，通过速度分析得到转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff和转换波各向异性参数χ_eff；在一个示例中，本领域技术人员可根据任意已知的手段实现步骤101。

步骤102，在第k个采样点处，针对共转换点道集中的各个地震道进行以下处理：

步骤1021：根据转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff计算纵波叠加速度V_p2；

步骤1022：根据转换波双程旅行时t_c和垂直速度比γ₀计算纵波双程旅行时t_p0；

步骤1023：根据纵波双程旅行时t_p0和纵波叠加速度V_p2计算纵波层速度V_p,k；

步骤1024：根据转换波双程旅行时t_c、转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff和转换波各向异性参数χ_eff计算射线参数p；

步骤1025：根据纵波层速度V_p,k、射线参数p计算转换波入射角θ_k，其中，k为正整数，表示采样点编号且对应于层编号，转换波入射角θ_k为转换波在第k层处的入射角；

步骤103，取k为1到n，并针对每个k值执行步骤102，以得到在每个采样点处针对各个地震道的转换波入射角，将转换波入射角相同的地震道归为同一入射角道集，得到转换波角道集，其中n为总采样点数。

该实施例通过转化将转换波自激自收时间和叠加速度映射为相应的纵波自激自收时间和层速度，用转换波射线参数法求取转换波入射角，提高了叠前AVO反演地层弹性参数精度和地震资料储层含气性预测精度。

为便于理解该实施例，图2示出了典型的层状介质转换波射线路径的示意图，S是炮点，R是接收点，Xc是炮点到转换点距离，x是炮点到检波点距离，也即是偏移距x，纵波和横波层速度分别为V_p,1,…,V_p,n和V_s,1,…,V_s,n，各层纵波入射角和横波反射角分别为θ₁,…,θ_n和其中n是总采样点数即总层数。入射纵波(P)经过的每一层的射线参数为p(1)至p(n)，反射横波(S)经过的每一层的射线参数为p(n+1)至p(2n)。根据Snell定律，地震波传播路径的上每一点的射线参数为常数，即p(1)＝……＝p(n)＝p(n+1)＝……p(2n)＝p。

如本领域技术人员能够理解的，每个采样点对应着1-n层中的一层，第k个采样点对应于第k层，相同的入射角在不同采样点(即不同的层)对应的地震道不同(即偏移距不同)，所以这里在每个采样点处(即第k个采样点，k取1到n)，针对共转换点道集中的各个地震道(即各个偏移距)求对应于该采样点的转换波入射角θ_k，这样可以求出所有采样点处针对所有地震道的转换波入射角，然后再把相同转换波入射角的地震道放入同一道集，即形成入射角道集。

本实施例的理论依据如下。

地震波传播路径上每一层地震波入射角如射线参数的关系为：

因此

θ_k＝asin(V_p,kp) (1)

根据地震波旅行时理论，在共中心点道集内，射线参数p是旅行时t_c对偏移距x的导数。

$p=\frac{{\mathrm{dt}}_c}{dx}---\left(2\right)$

所以

${\mathrm{\θ}}_k=asin\left(V_{p,k}\frac{{\mathrm{dt}}_c}{dx}\right)---\left(3\right)$

根据各向异性介质理论，转换波共中心点道集不同偏移距x的旅行时t_c满足：

$t_c^2=t_{c0}^2+\frac{x^2}{V_{c2}^2}+\frac{A_4{x}^4}{1+A_5{x}^2}---\left(4\right)$

$A_4=-\frac{{({\mathrm{\γ}}_0{\mathrm{\γ}}_{eff}-1)}^2+8(1+{\mathrm{\γ}}_0){\mathrm{\γ}}_{eff}}{4t_{c0}^2{V}_{c2}^4{\mathrm{\γ}}_0{(1+{\mathrm{\γ}}_{eff})}^2}---\left(5\right)$

$A_5=\frac{A_4{V}_{c2}^2(1+{\mathrm{\γ}}_0){\mathrm{\γ}}_{eff}\[({\mathrm{\γ}}_0-1){\mathrm{\γ}}_{eff}^2+2{\mathrm{\γ}}_{eff}\]}{({\mathrm{\γ}}_0-1){\mathrm{\γ}}_{eff}^2(1-{\mathrm{\γ}}_0{\mathrm{\γ}}_{eff})-2(1+{\mathrm{\γ}}_0){\mathrm{\γ}}_{eff}{\mathrm{\χ}}_{eff}}---\left(6\right)$

其中t_c为转换波双程旅行时，t_c0为垂向转换波双程旅行时，V_c2为转换波叠加速度，γ₀为垂直速度比，γ_eff为等效速度比，χ_eff为转换波各向异性系数。上述参数满足以下关系：

${\mathrm{\γ}}_0=\frac{V_{p0}}{V_{s0}}=\frac{t_{s0}}{t_{p0}}=2(\frac{t_{c0}}{t_{p0}}-1)---\left(7\right)$

${\mathrm{\γ}}_{eff}=\frac{{\mathrm{\γ}}_2^2}{{\mathrm{\γ}}_0},{\mathrm{\γ}}_2=\frac{V_{p2}}{V_{s2}}---\left(8\right)$

$V_{c2}^2=\frac{1}{1+{\mathrm{\γ}}_0}{V}_{p2}^2+\frac{{\mathrm{\γ}}_0}{1+{\mathrm{\γ}}_0}{V}_{s2}^2---\left(9\right)$

整理得：

$t_{p0}=\frac{t_{c0}}{1+\frac{{\mathrm{\γ}}_0}{2}}---\left(10\right)$

$V_{p2}^2=V_{c2}^2\frac{{\mathrm{\γ}}_{eff}(1+{\mathrm{\γ}}_0)}{1+{\mathrm{\γ}}_{eff}}---\left(11\right)$

其中V_p0是纵波垂向速度，V_s0是横波垂向速度，V_p2为纵波叠加速度，V_s2是横波叠加速度，t_p0为纵波垂向双程旅行时，t_s0是横波垂向双程旅行时，γ₂是纵波和横波叠加速度比。通过DIX公式，基于本领域技术人员已知的手段可得纵波速度层速度V_p,k。

把式(5)代入式(3)得射线参数：

$p=\frac{1}{t_c}(\frac{x}{V_{c2}^2}+\frac{2A_4{x}^3(1+A_5{x}^2)+A_4{A}_5{x}^5}{{(1+A_5{x}^2)}^2})---\left(12\right)$

因此转换波入射角

${\mathrm{\θ}}_k=asin\left(V_{p,k}p\right)=asin\left(\frac{V_{p,k}}{t_c}\right(\frac{x}{V_{c2}^2}+\frac{2A_4{x}^3(1+A_5{x}^2)+A_4{A}_5{x}^5}{{(1+A_5{x}^2)}^2}\left)\right)---\left(13\right)$

其中k为正整数，代表所求得的转换波入射角所对应的层数，该层数与采样点编号相对应。在一个示例中，步骤1021中的根据转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff计算纵波叠加速度V_p2，可通过将这些参数带入上述公式(12)来实现。

在一个示例中，步骤1022中的根据垂向转换波双程旅行时t_c0和垂直速度比γ₀计算纵波双程旅行时t_p0可通过将这些参数带入上述公式(11)来实现。

在一个示例中，步骤1023中可基于DIX公式，利用本领域技术人员已知的方式根据纵波双程旅行时t_p0和纵波叠加速度V_p2计算纵波层速度V_p,k。

在一个示例中，步骤1024中的根据转换波双程旅行时t_c、转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff和转换波各向异性参数χ_eff计算射线参数p可通过将这些参数带入公式(13)来实现。

在一个示例中，步骤1025中的根据纵波层速度V_p,k和射线参数p计算转换波入射角θ_k可通过将这些参数带入公式(14)来实现。

在一个示例中，可以在每个采样点处遍历所有偏移距值来遍历所有地震道，并在遍历过程中针对每个地震道执行本实施例的步骤102。

根据本发明的另一实施例，提出了一种横向各向同性介质转换波角道集抽取装置，该装置可包括：用于基于叠前转换波的共转换点道集，通过速度分析得到转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff和转换波各向异性参数χ_eff的部件；用于在第k个采样点处，针对共转换点道集中的各个地震道进行以下处理的部件：

-根据转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff计算纵波叠加速度V_p2；

-根据转换波双程旅行时t_c和垂直速度比γ₀计算纵波双程旅行时t_p0；

-根据纵波双程旅行时t_p0和纵波叠加速度V_p2计算纵波层速度V_p,k；

-根据转换波双程旅行时t_c、转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff和转换波各向异性参数χ_eff计算射线参数p；

-根据纵波层速度V_p,k、射线参数p计算转换波入射角θ_k，其中，k为正整数，表示采样点编号且对应于层编号，转换波入射角θ_k为转换波在第k层处的入射角；以及

用于取k为1到n，并针对每个k值执行上述处理，以得到在每个采样点处针对各个地震道的转换波入射角，将转换波入射角相同的地震道归为同一入射角道集，得到转换波角道集的部件，其中n为总采样点数。

应用示例

为说明本发明实施例的有益效果，以下通过图2-5给出了一个具体应用示例。本领域技术人员应理解，该应用示例仅为了便于理解本发明实施例以及说明本发明实施例的效果，其涉及到的任何细节均不构成对本发明的限制。

图3示出了一个示例性的层状介质模型的纵波速度、横波速度和密度曲线示意图，虚点划线是纵波速度，实线是横波速度，虚点线是密度。速度单位是千米/秒，密度单位是克/立方厘米。通过射线追踪，得到不同层位不同偏移距转换波反射时间和入射角，再根据Zoeppritz计算不同偏移距反射振幅。通过转换波速度分析，得到转换波叠加速度V_c2、垂直速度比γ₀、等效速度比γ_eff和转换波各向异性参数χ_eff。图4是正演模拟转换波道集经过动校正处理后的道集的示意图。第一层的最大偏移距(2800米)对应的入射角达到60度，第九层的最大偏移距(2800米)对应的入射角为37度。在每个采样点对各个偏移距进行上述步骤102，可以得到每个采样点的各个偏移距的入射角。本发明实施例计算的入射角与射线追踪的真正入射角非常吻合，图5和图6分别是层状介质模型的第5层和第9层的本发明实施例计算的入射角与射线追踪的真正入射角的比较，可见本发明实施例计算的入射角完全正确。图7示出了根据第1层和第9层的入射角分布范围，抽取了入射角范围为5至30度，间隔5度的转换波入射角道集的一个示例的示意图。

本公开可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本公开的各个方面的计算机可读程序指令。

计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是――但不限于――电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携 式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身，诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波(例如，通过光纤电缆的光脉冲)、或者通过电线传输的电信号。

这里所描述的计算机可读程序指令可以从计算机可读存储介质下载到各个计算/处理设备，或者通过网络、例如因特网、局域网、广域网和/或无线网下载到外部计算机或外部存储设备。网络可以包括铜传输电缆、光纤传输、无线传输、路由器、防火墙、交换机、网关计算机和/或边缘服务器。每个计算/处理设备中的网络适配卡或者网络接口从网络接收计算机可读程序指令，并转发该计算机可读程序指令，以供存储在各个计算/处理设备中的计算机可读存储介质中。

用于执行本公开操作的计算机程序指令可以是汇编指令、指令集架构(ISA)指令、机器指令、机器相关指令、微代码、固件指令、状态设置数据、或者以一种或多种编程语言的任意组合编写的源代码或目标代码，所述编程语言包括面向对象的编程语言—诸如Smalltalk、C++等，以及常规的过程式编程语言—诸如“C”语言或类似的编程语言。计算机可读程序指令可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络—包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。在一些实施例中，通过利用计算机可读程序指令的状态信息来个性化定制电子电路，例如可编程逻辑电路、现场可编程门阵列(FPGA)或可编程逻辑阵列(PLA)，该电子电路可以执行计算机可读程序指令，从而实现本公开的各个方面。

这里参照根据本公开实施例的方法、装置(系统)和计算机程序产品的流程图和/或框图描述了本公开的各个方面。应当理解，流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合，都可以由计算机可读程序指令实现。

这些计算机可读程序指令可以提供给通用计算机、专用计算机或其它可编程数据处理装置的处理器，从而生产出一种机器，使得这些指令在通过计算机或其它可编程数据处理装置的处理器执行时，产生了实现流程图和/或框图中的一个或多个方框中规定的功能/动作的装置。也可以把这些计算机可读程序指令存储在计算机可读存储介质中，这些指令使得计算机、可编程数据处理装置和/或其他设备以特定方式工作，从而，存储有指令的计算机可读介质则包括一个制造品，其包括实现流程图和/或框图中的一个或多个方框中规定的功能/动作的各个方面的指令。

也可以把计算机可读程序指令加载到计算机、其它可编程数据处理装置、或其它设备上，使得在计算机、其它可编程数据处理装置或其它设备上执行一系列操作步骤，以产生计算机实现的过程，从而使得在计算机、其它可编程数据处理装置、或其它设备上执行的指令实现流程图和/或框图中的一个或多个方框中规定的功能/动作。

附图中的流程图和框图显示了根据本公开的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或指令的一部分，所述模块、程序段或指令的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

以上已经描述了本公开的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。