反演裂缝密度的方法及系统与流程

本公开涉及油气地球物理勘探领域，更具体地，涉及一种反演裂缝密度的方法及系统。

背景技术：

在油气地球物理勘探领域，对裂缝性储层研究已成为日益重要的内容，而裂缝密度是裂缝性储层定量表征中的重要参数。

在现有的预测裂缝性质的方法中，最常见的是椭圆拟合方法，对于HTI(atransversely isotropic medium with a horizontal axis of symmetry(具有对称水平轴的横向各向同性介质))介质，纵波的地震属性(叠加速度和反射振幅)随方位变化呈椭圆形特征，长轴方向反映裂缝发育主要方位，长、短轴之比反映了裂缝发育程度。但此方法只能估计裂缝的发育程度，并不能定量的表征裂缝密度。Isabel Varela对含有AVAZ(Amplitude variation with incident angle and azimuth(振幅随入射角和方位角的变化))特征的数据运用SVD(sigular value decomposition奇异值分解(Singular Value Decomposition))反演裂缝密度。Morten Jakobsen等通过AVOZ分析估计裂缝储层渗透率。朱培民等提出用2条正交测线上的P波含有AVO特征的观测数据反演Thomsen参数，进而估计裂缝密度的方法。

发明人发现，上述几种有关的反演裂缝密度的方法均是建立在可以从反射信息中取得准确的反射系数的基础上，然而实际上地震勘探得到的是振幅信息，提取地震子波得到反射系数的过程不可避免也会对后续反演的精度产生影响，因此，有必要开发一种精确反演裂缝密度的方法。

公开于本公开背景技术部分的信息仅仅旨在加深对本公开的一般背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人 员所公知的现有技术。

技术实现要素：

本公开提出了一种反演裂缝密度的方法及系统，该方法可以基于相交测线上的反射振幅差值构造目标函数，基于目标函数进行反演得到与裂缝密度有关的Thomsen参数，以及基于Thomsen参数与裂缝密度的关系获得裂缝密度，从而实现了对裂缝密度的精确反演。

根据本公开的一方面，提出了一种反演裂缝密度的方法，该方法可以包括以下步骤：基于地震资料获取相交测线的纵波反射振幅差值，并基于拟合模型获得纵波反射系数差值；基于纵波反射系数差值和纵波反射振幅差值，构造目标函数fitness；以及基于目标函数fitness进行反演计算以获得裂缝密度，其中，目标函数fitness可以表示为：

其中，表示方位角；η_j表示相交测线的方位角差值；θ_k表示入射角；表示方位角分别为和入射角为θ_k时的纵波反射振幅差值；表示的最大值；表示反射系数差值或振幅差值最大时对应的入射角；Δγ、Δδ^(V)与Δε^(V)表示Thomsen参数，其中前置符号“Δ”代表上下两层Thomsen参数的差值；表示拟合模型中两条相交测线的纵波反射系数差值；表示的最大值；ω^ij表示权系数；I表示地震资料中方位角的测线总数量；J表示地震资料中方位角为时与其相交的测线的总数量；以及K 表示地震资料中入射角θ的总数量。

根据本公开的另一方面，提出了一种反演裂缝密度的系统，该系统可以包括以下单元：用于基于地震资料获取相交测线的纵波反射振幅差值，并基于拟合模型获得纵波反射系数差值的单元；用于基于纵波反射系数差值和纵波反射振幅差值，构造目标函数fitness的单元；以及用于基于目标函数fitness进行反演计算以获得裂缝密度的单元，其中，目标函数fitness可以表示为：

其中，表示方位角；η_j表示相交测线的方位角差值；θ_k表示入射角；表示方位角分别为和入射角为θ_k时的纵波反射振幅差值；表示的最大值；θ_km表示反射系数差值或振幅差值最大时对应的入射角；Δγ、Δδ^(V)与Δε^(V)表示Thomsen参数；表示拟合模型中两条相交测线的纵波反射系数差值；表示的最大值；ω^ij表示权系数；I表示地震资料中方位角的测线总数量；J表示地震资料中方位角为时与其相交的测线的总数量；以及K表示地震资料中入射角θ的总数量。

本公开的方法和装置具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施方案中将是显而易见的，或者将在并入本文中的附图和随后的具体实施方案中进行详细陈述，这些附图和具体实施方案共同用于解释本公开的特定原理。

附图说明

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

图1示出了根据本公开的一个实施方案的反演裂缝密度的方法的步骤的流程图。

图2a-2b为示出了根据本公开的一个示例的Thomsen参数与裂缝密度关系的示意图，其中，图2a为裂缝含水，图2b为裂缝含气。

图3a-3e为示出了根据本公开的一个示例的Z分量地震记录的示意图，其中，图3a为方位角0°；图3b为方位角30°；图3c为方位角45°；图3d为方位角60°；以及图3e为方位角90°。

图4a-4h为示出了根据本公开的一个示例的相交测线上HTI介质反射振幅差值的示意图，其中，图4a为方位角0°和90°；图4b为方位角30°和60°；图4c为方位角30°和90°；图4d为方位角0°和30°；图4e为方位角0°和45°；图4f为方位角45°和90°；图4g为方位角45°和60°；以及图4h为方位角30°和45°。

具体实施方式

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

实施方案1

图1示出了根据本公开的一个实施方案的反演裂缝密度的方法的流程图。根据本公开的实施方案的反演裂缝密度的方法可以包括以下步骤：步骤101，基于地震资料获取相交测线的纵波反射振幅差值和拟合模型中纵波反射系数差值；步骤102，基于纵波反射系数差值和纵波反射振幅差值，构造目标函数fitness； 以及步骤103，基于目标函数fitness进行反演计算以获得裂缝密度，其中，目标函数fitness可以表示为：

其中，表示方位角；η_j表示相交测线的方位角差值；θ_k表示入射角；表示方位角分别为和入射角为θ_k时的纵波反射振幅差值；表示的最大值；θ_km表示反射系数差值或振幅差值最大时对应的入射角；Δγ、Δδ^(V)与Δε^(V)表示Thomsen参数，其中前置符号“Δ”代表上下两层Thomsen参数的差值；表示拟合模型中两条相交测线的纵波反射系数差值；表示的最大值；ω^ij表示权系数；I表示地震资料中方位角的测线总数量；J表示地震资料中方位角为时与其相交的测线的总数量；以及K表示地震资料中入射角θ的总数量。

本实施方案可以通过利用纵波反射系数差值和纵波反射振幅差值来构造目标函数，并基于目标函数进行计算以获得裂缝密度，实现了对裂缝密度的精确反演。

获取纵波反射系数差值和纵波反射振幅差值

在一个示例中，可以基于拟合模型获得纵波反射系数差值。其中，可以使用Rüger反射系数近似公式来计算纵波反射系数差值。

根据Rüger反射系数近似公式，当方位角为和时，纵波反射系数差值可以用如下公式(1)表示：

其中，R可以表示纵波反射系数，可以表示方位角为时的测线的纵波反射系数，可以表示方位角为时的测线的纵波反射系数，ΔR可以表示方位角为和时的相交测线的纵波反射系数差值，α可以表示HTI介质各向同性面的纵波速度，β可以表示HTI介质各向同性面的SH波速度，可以表示HTI介质各向同性面的纵波速度的平均值，可以表示HTI介质各向同性面的SH波速度的平均值，θ可以表示入射角。其中，γ、δ^(V)与ε^(V)可以表示Thomsen参数，γ可以表示在水平方向传播的SH波与SV波速度的差别，，δ^(V)可以表示P波相速度在垂直入射时的二阶导数，ε^(V)可以表示纵波速度在垂直和水平方向的差别，Thomsen参数上角标V代表区分针对VTI介质的Thomsen参数，前置符号“Δ”代表上下两层Thomsen参数的差值。

本领域技术人员应该理解，获得纵波反射系数差值的方法并不局限于此，而是可以采用本领域技术人员所公知的各种反射系数近似方式来获得纵波反射系数的差值。

从公式(1)可以获得纵波反射系数差值，并可以使用表示两条相交测线的纵波反射系数差值，其中(i＝1,2,...,I；j＝1,2,...,J；k＝1,2,...,K)，I可以表示地震资料中方位角的测线总数量，J可以表示地震资料中方位角为时与其相交的测线的总数量，K可以表示地震资料中入射角θ的总数量。

在一个示例中，可以基于地震资料获取相交测线的纵波反射振幅差值。其中，可以使用观测地震数据中目的层位上的相交测线反射振幅来计算纵波反射 振幅差值，其可以使用如下公式(2)表示：

其中，A可以表示从反射信息中提取到的纵波反射振幅，可以表示方位角分别为和入射角为θ_k时的纵波反射振幅差值。本领域技术人员应当理解，本公开并不限于此，而可以使用本领域公知的任何地震数据处理与解释手段来计算纵波反射振幅差值。

构造目标函数

在一个示例中，可以基于纵波反射系数差值和纵波反射振幅差值，构造目标函数。所构建的目标函数可以用如下公式(3)表示：

其中，ω^ij可以表示权系数，可以表示在K个入射角中反射系数差值或振幅差值最大时对应的入射角，可以表示方位角分别为和入射角为θ_k时的纵波反射振幅差值的最大值，可以表示两条相交测线的纵波反射系数差值的最大值。

可以使用如下公式(4)表示：

可以使用如下公式(5)表示：

其中，从拟合模型中获得的纵波反射系数差值和从地震资料中获取的纵波反射振幅差值两者均可以是无量纲的物理量。由于反射系数与反射振幅呈线性 关系，因此可以利用反射振幅差值中的最大值点为参考，对应拟合模型中反射系数差值的最大值，可以计算出反射振幅差值与反射系数差值的比值关系，再根据此比值关系将其它反射系数差值转换为拟合模型的反射振幅差值。因而公式(3)中的后一项可以为拟合模型的反射振幅差值，与前一项(实际的纵波反射振幅差值)可以是同一个物理量。因此，可以使用与的差值最小的情况(也即目标函数fitness取最小值)下的Thomsen参数计算裂缝密度。

在一个示例中，权系数ω^ij可以用于控制一组不同方位角的纵波反射振幅差值的可信程度。其中，权系数ω^ij可以由用户根据纵波反射振幅差值的可信程度自行设定，例如ω^ij可以设定为0-1之间的值。某一入射角或方位角的观测数据的质量越好，也即可信程度越高，则权系数可以越大。具体地，可以通过地震数据的信号与噪音情况、某一地震面元上各个方位上的地震数据是否足够多来判断数据是否可信，权系数可以是主观参数。在实际情况下，如果该地区数据采集是规律的，即每个面元的炮点和检波点分布一致，则某一个i或者某一个j对应的权系数是可以一样的。

获得裂缝密度

在一个示例中，可以基于目标函数fitness进行反演计算以获得裂缝密度。通过非线性反演算法求得公式(3)所示的目标函数fitness的最小值，基于目标函数fitness的最小值可以获得与该最小值对应的一组Thomsen参数Δγ，Δδ^(V)和Δε^(V)的值。本领域技术人员应当理解，求得目标函数fitness的最小值的方法并不限于此，而可以使用本领域公知的任何方法求得目标函数fitness的最小值。

如果上层介质为均匀各向同性介质，下层为HTI介质，则Δγ，Δδ^(V)和Δε^(V) 均等于下层介质的Thomsen参数γ，δ^(V)和ε^(V)。因此，根据所获得的Thomsen参数Δγ，Δδ^(V)和Δε^(V)的值，可以基于如下公式(6)和(7)中的HTI介质Thomsen参数与裂缝密度e的关系，获得裂缝密度e：

$e=\frac{3(3-2g)\mathrm{\γ}}{8(1+2\mathrm{\γ})}---\left(6\right)$

$g=\frac{\mathrm{\μ}}{\mathrm{\λ}+2\mathrm{\μ}}=\frac{V_S^2}{V_P^2}---\left(7\right)$

其中，V_P和V_S可以为背景介质的纵波和横波速度，e可以表示裂缝密度，λ和μ是拉梅常数，g可以表示横波与纵波速度比值的平方。

应用示例

为便于理解本公开实施方案的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本公开，其任何具体细节并非意在以任何方式限制本公开。

图2a和图2b为示出了根据本公开的一个示例的Thomsen参数与裂缝密度关系的示意图。其中，图2a为裂缝含水，图2b为裂缝含气。在图2a和图2b中，纵波与横波背景速度分别为3800m/s和2000m/s，密度为2.4g/cm3，裂缝纵横比为10-3，纵轴可以表示Thomsen参数，图2a的横轴可以表示含水裂缝密度，图2b的横轴可以表示含气裂缝密度。

图3a-3e为示出了根据本公开的一个示例的Z分量地震记录的示意图，其中，图3a为方位角0°；图3b为方位角30°；图3c为方位角45°；图3d为方位角60°；以及图3e为方位角90°。其中，模型上层为均匀各向同性介质，下层为HTI介质。

图4a-4h为示出了根据本公开的一个示例的相交测线上HTI介质反射振幅差值的示意图，其中，图4a为方位角0°和90°；图4b为方位角30°和60°；图4c为方位角30°和90°；图4d为方位角0°和30°；图4e为方位角0°和 45°；图4f为方位角45°和90°；图4g为方位角45°和60°；以及图4h为方位角30°和45°。

如图4a-4h所示，反演得到的Thomsen参数结果：γ相对误差为1.3％；δ^(V)相对误差为28％；ε^(V)相对误差为26.6％。其中相对误差公式为：|(真实值-反演值)/真实值|×100％。显然，三个Thomsen参数反演值中，γ反演结果最接近真实值。根据公式(4)中裂缝密度e与γ的关系，可估计裂缝密度为相对误差为9％。数值反演结果表明：相交测线上反射振幅差值正演值与反演值拟合情况较好，基于纵波反射系数差值和纵波反射振幅差值来反演裂缝密度的方法是可行且稳定的。

本领域技术人员应理解，上面对本公开的实施方案的描述的目的仅为了示例性地说明本公开的实施方案的有益效果，并不意在将本公开的实施方案限制于所给出的任何示例。

实施方案2

在该实施方案中，提供了一种反演裂缝密度的系统，该系统可以包括以下单元：用于基于地震资料获取相交测线的纵波反射振幅差值，并基于拟合模型获得纵波反射系数差值的单元；用于基于纵波反射系数差值和纵波反射振幅差值，构造目标函数fitness的单元；以及用于基于目标函数fitness进行反演计算以获得裂缝密度的单元。其中，目标函数fitness可以表示为：

其中，表示方位角；η_j表示相交测线的方位角差值；θ_k表示入射角；表示方位角分别为和入射角为θ_k时的纵波反射振幅差值；表示的最大值；表示反射系数差值或振幅差值最大时对 应的入射角；Δγ、Δδ^(V)与Δε^(V)表示Thomsen参数；表示拟合模型中两条相交测线的纵波反射系数差值；表示的最大值；ω^ij表示权系数；I表示地震资料中方位角的测线总数量；J表示地震资料中方位角为时与其相交的测线的总数量；以及K表示地震资料中入射角θ的总数量。

本实施方案可以通过利用纵波反射系数差值和纵波反射振幅差值来构造目标函数，并基于目标函数进行计算以获得裂缝密度，实现了对裂缝密度的精确反演。

在一个示例中，可以使用Rüger反射系数近似公式来计算纵波反射系数差值。

在一个示例中，可以使用观测地震数据中目的层位上的相交测线的反射振幅来计算纵波反射振幅差值。

在一个示例中，目标函数fitness中的权系数可以用于控制一组不同方位角的纵波反射振幅差值的可信程度。

在一个示例中，基于所述目标函数fitness进行反演计算以获得裂缝密度可以包括：通过非线性反演算法获得目标函数fitness的最小值；基于目标函数fitness的最小值获得Thomsen参数；以及基于Thomsen参数进行计算以获得裂缝密度。

本领域技术人员应理解，上面对本公开的实施方案的描述的目的仅为了示例性地说明本公开的实施方案的有益效果，并不意在将本公开的实施方案限制于所给出的任何示例。

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

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

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

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

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

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

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

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

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