CN109581492A - 基于地震波模拟的岩石物理参数计算方法及系统 - Google Patents

Abstract

公开了一种基于地震波模拟的岩石物理参数计算方法及系统。该方法可以包括：根据岩心数据，建立初始大尺度储层模型，进而基于裂缝数据与孔隙度数据，获取大尺度储层模型；基于大尺度储层模型，建立小尺度网格，获取等效弹性参数模型；基于等效弹性参数模型，获取储层的各向异性参数。本发明通过数值计算方法获得复杂储层的岩石物理弹性参数，避免了获取不同地质特征岩芯的困难，得到小尺度非均匀导致的宏观物理参数变化。

Description

基于地震波模拟的岩石物理参数计算方法及系统

技术领域

本发明涉及油气勘探中岩石物理建模技术领域，更具体地，涉及一种基于地震波模拟的岩石物理参数计算方法及系统。

背景技术

岩石物理研究是沟通石油地质参数和地震响应的关键手段。由反射地震资料获取地下介质的速度、密度、泊松比、各向异性参数，进而由这些参数获得地下介质的孔隙度、裂缝发育和走向、孔隙和裂缝含流体(或气体)情况，就有望实现由反射地震资料直接识别储层，实现岩性勘探的目标。传统的岩石物理研究主要是实验测量和等效介质理论这两个途径。实验方法通过测量特定温压条件下岩芯样品的岩石物理参数特征及其变化规律，应用统计方法建立岩石物理弹性参数与储层参数关系的经验模型，给出由储层参数估算弹性参数的近似算法和公式以及相反。等效介质理论利用静力学方法，获取不同矿物组分的介质，以及介质包含裂缝、孔隙和孔隙含流体(或气体)情况下的等效物理模量，以此获得裂缝、孔隙和孔隙流体对岩石物理弹性参数影响的规律。两类方法已取得了很好的结果，但在实际应用中也遇到了一些困难。实验方法遇到的主要困难是岩芯样品和尺度的限制。实际采集的系列岩芯样品一般较难满足对特定地质特征，如孔隙度，在特定范围内变化的要求，因此由测量结果回归经验模型就存在数据不充分的问题。地震勘探中的地震波波长是百米量级，为准确获得地震勘探频率范围的岩石物理弹性参数，理想的情况是岩芯样品也是百米量级；而这无论从实验装置和岩芯采样都是不能实现的；因此强烈依赖于尺度的岩石物理弹性参数是难以通过实验直接获得的，如衰减特征。基于等效介质理论的理论方法的主要问题是理论模型的严格假设；现行的等效介质理论都对介质、孔隙和裂缝分布做了均匀或均匀分布的假设，因此将这些理论应用于具有明显小尺度非均匀特征的储层(例如页岩储层)时将产生偏差。

作为对实验和理论方法的补充，计算岩石物理方法近年来得到重视和发展。目前比较活跃的是岩芯成像与数值模拟相结合的方法(Dvorkin et al.，2011；Diaz et al.，2013)。这类方法首先利用束扫描电镜(FIB-SIM)和CT成像技术，获取岩芯的小尺度(微米级)精细构造，然后利用数值模拟获得岩石物理参数(渗透率、电阻率、弹性模量)。尽管相对实验测量这类方法有更多的灵活性，但由于电镜尺寸和计算能力的限制，这类方法所涉及的岩芯尺度是毫米量级，较实验测量的厘米量级又低一个量级。由于尺度太小，这类方法在获取弹性参数时利用的是静力学模拟方法。理论上，针对更大尺度岩芯应用数值模拟方法来获取岩石物理参数是可行的，但计算能力和获取大尺度精细地质模型的困难限制了这方面的应用。已有的工作(Deeks et al.，2012；Dvorkin et al.，2012)主要是针对完全的随机模型展开，由于计算机能力的限制，也还没有达到地震波波长的尺度。因此，有必要开发一种基于地震波模拟的岩石物理参数计算方法及系统。

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

发明内容

本发明提出了一种基于地震波模拟的岩石物理参数计算方法及系统，其能够通过数值计算方法获得复杂储层的岩石物理弹性参数，避免了获取不同地质特征岩芯的困难，得到小尺度非均匀导致的宏观物理参数变化。

根据本发明的一方面，提出了一种基于地震波模拟的岩石物理参数计算方法。所述方法可以包括：根据岩心数据，建立初始大尺度储层模型，进而基于裂缝数据与孔隙度数据，获取大尺度储层模型；基于所述大尺度储层模型，建立小尺度网格，获取等效弹性参数模型；基于所述等效弹性参数模型，获取储层的各向异性参数。

优选地，建立所述初始大尺度储层模型包括：设定单元类型，进而确定每个位置的单元类型的概率；对所述储层的底层进行初始化处理，确定所述底层的单元类型；基于所述每个位置的单元类型的概率与所述底层的单元类型，建立所述初始大尺度储层模型。

优选地，每个位置的单元类型的概率，即Z_i,j＝S_k的概率为：

其中，p_k表示每个位置的单元类型的概率，i表示二维空间的纵向变化，j表示二维空间的深度变化，Z_i,j表示(i，j)处的单元类型，S_k(k＝1,n)代表定义的单元类型，n为定义的单元类型总数，f为定义的单元类型编号，p_mk是纵向转移概率矩阵中(m，k)处的元素，r_lk是横向转移概率矩阵中(l，k)处的元素，N是横向的单元总数，代表模型的横向大小；l源于Z_i-1,j＝S_l，m源于Z_i,j-1＝S_m，q源于Z_N,j＝S_q，p_mf为纵向转移概率矩阵中(m，f)处的元素，r_kq、r_lf、r_fq分别为横向转移概率矩阵中(k，q)、(l，f)、(f，q)处的元素。

优选地，所述底层的单元类型的概率，即Z_i,1＝S_k的概率为：

优选地，获取所述弹性参数等效模型包括：基于所述大尺度储层模型，获取所述储层的矿物类型；根据自洽近似理论，计算每一种矿物类型的弹性模量；根据Gassmann理论，实现孔隙的流体替换，计算流体的弹性模量；根据所述每一种矿物类型的弹性模量与所述流体的弹性模量，获取所述储层的裂缝数据，进而获取所述弹性参数等效模型。

优选地，所述自洽近似理论为：

其中，表示矿物类型的等效体积模量，表示矿物类型的等效剪切模量，N代表共有N种矿物类型，x_i是第i种矿物的体积比，K_i表示第i种矿物的等效体积模量，μ_i表示第i种矿物的等效剪切模量，P^*i和Q^*i是第i种矿物在模量为和的背景介质中的几何系数。

优选地，所述各向异性参数为：

其中，δ、ε表示各向异性参数，(v_p)_v表示储层垂直方向的纵波速度，v_p(θ)表示储层沿角度θ方向的纵波速度。

根据本发明的另一方面，提出了一种基于地震波模拟的岩石物理参数计算系统，可以包括：存储器，存储有计算机可执行指令；处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：根据岩心数据，建立初始大尺度储层模型，进而基于裂缝数据与孔隙度数据，获取大尺度储层模型；基于所述大尺度储层模型，建立小尺度网格，获取等效弹性参数模型；基于所述等效弹性参数模型，获取储层的各向异性参数。

优选地，建立所述初始大尺度储层模型包括：设定单元类型，进而确定每个位置的单元类型的概率；对所述储层的底层进行初始化处理，确定所述底层的单元类型；基于所述每个位置的单元类型的概率与所述底层的单元类型，建立所述初始大尺度储层模型。

优选地，获取所述弹性参数等效模型包括：基于所述大尺度储层模型，获取所述储层的矿物类型；根据自洽近似理论，计算每一种矿物类型的弹性模量；根据Gassmann理论，实现孔隙的流体替换，计算流体的弹性模量；根据所述每一种矿物类型的弹性模量与所述流体的弹性模量，获取所述储层的裂缝数据，进而获取所述弹性参数等效模型。

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

附图说明

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

图1示出了根据本发明的基于地震波模拟的岩石物理参数计算方法的步骤的流程图。

图2示出了根据本发明的一个实施例的储层模型的示意图。

图3a和图3b分别示出了根据本发明的一个实施例的3％TOC含量储层的垂直和水平速度分量炮记录的示意图。

图4示出了根据本发明的一个实施例的页岩储层生成的数值模型的纹层状部分的局部示意图。

图5示出了根据本发明的一个实施例的页岩储层生成的数值模型的层状部分的局部示意图。

图6示出了根据本发明的一个实施例的提取弹性参数的弹性介质纵波速度模型的示意图。

图7a和图7b分别示出了根据本发明的一个实施例的水平平面纵波入射0.1s和0.25s时刻的波场快照的示意图。

图8a和图8b分别示出了根据本发明的一个实施例的斜入射平面纵波0.1s和0.29s时刻的波场快照的示意图。

图9a和图9b分别示出了根据本发明的一个实施例的水平平面横波入射0.2s和0.38s时刻的波场快照的示意图。

图10示出了根据本发明的一个实施例的储层垂直方向纵波速度和横波速度随TOC含量变化曲线图。

图11示出了根据本发明的一个实施例的储层各向异性参数ε参数和δ参数随TOC含量变化曲线图。

具体实施方式

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

图1示出了根据本发明的基于地震波模拟的岩石物理参数计算方法的步骤的流程图。

在该实施例中，根据本发明的基于地震波模拟的岩石物理参数计算方法可以包括：

步骤101，根据岩心数据，建立初始大尺度储层模型，进而基于裂缝数据与孔隙度数据，获取大尺度储层模型；在一个示例中，建立初始大尺度储层模型包括：设定单元类型，进而确定每个位置的单元类型的概率；对储层的底层进行初始化处理，确定底层的单元类型；基于每个位置的单元类型的概率与底层的单元类型，建立初始大尺度储层模型。

在一个示例中，每个位置的单元类型的概率，即Z_i,j＝S_k的概率为：

其中，p_k表示每个位置的单元类型的概率，i表示二维空间的纵向变化，j表示二维空间的深度变化，Z_i,j表示(i，j)处的单元类型，S_k(k＝1,n)代表定义的单元类型，n为定义的单元类型总数，f为定义的单元类型编号，p_mk是纵向转移概率矩阵中(m，k)处的元素，r_lk是横向转移概率矩阵中(l，k)处的元素，N是横向的单元总数，代表模型的横向大小；l源于Z_i-1,j＝S_l，m源于Z_i,j-1＝S_m，q源于Z_N,j＝S_q，p_mf为纵向转移概率矩阵中(m，f)处的元素，r_kq、r_lf、r_fq分别为横向转移概率矩阵中(k，q)、(l，f)、(f，q)处的元素。

在一个示例中，底层的单元类型的概率，即Z_i,1＝S_k的概率为：

其中，r_lq为横向转移概率矩阵中(l，q)处的元素。

具体地，储层的精细数值建模是计算岩石物理方法应用的基础。为准确获得地震勘探频率范围的岩石物理弹性参数，特别是考虑储层非均匀导致的尺度效应，对百米量级的岩芯样品进行计算岩石物理的数值模拟是更合理的选择；而如何获得这样较大尺度的地质模型是这一思路能否实现的关键。页岩储层有别于其他种类储层的一个突出特点是其成层特征导致的小尺度非均匀性，这一非均匀性导致较大尺度的模型不能通过数值岩芯(毫米尺度)的简单放大“复制”来得到。为此，利用地层沉积的统计学模型，来获取地震波波长尺度的页岩储层精细(毫米级)数值模型；而岩芯样品，则是用来获得这一随机过程的关键控制因素的依据。

地质模型的非均匀是典型的分形(Fractal)结构，可以在不同尺度下无限细分。在建模中将垂向最小尺度定义在与页岩储层的纹层结构同量级(毫米)，而横向非均匀的尺度定义在“10米”量级。采用二维马尔科夫链模型生成地震波波长尺度的二维精细地质模型。为此，将以一定长度的“层”和“纹层”为单元，而每个“层”和“纹层”又可以由不同的岩相构成，不同类型岩相的“层”或“纹层”单元构成了组成地质模型的基本单元。由纵向和横向转移概率矩阵联合确定不同新事件(即基本单元)的发生概率，即可由这些基本单元随机构成二维地质模型。

定义二维空间位置(i，j)，其中j代表代表纵向变化(即沉积历史)，j增加表明深度变化，令Z_i,j表示空间位置(i，j)处的基本单元类型；S_k(k＝1,n)代表定义的基本单元，p_mk是纵向转移概率矩阵中(m，k)处的元素，r_lk是横向转移概率矩阵中(l，k)处的元素，则Z_i,j＝S_k的概率为公式(1)。

在没有井约束情况，为完成二维建模，需先应用一维马尔科夫链模型生成两端，即i＝1,N处的单元类型Z_1,j和Z_N,j。为在此后的弹性波模拟计算中在横向应用循环边界条件，定义Z_1,j＝S_N,j。Z_1,j＝S_k的概率为公式(3)

p_k＝p_mk (3)

其中，m源于Z_1,j-1＝S_m。

在应用公式(1)构建二维模型时，需对底层初始化，即确定最底层的单元类型。再次应用一维马尔科夫链模型，Z_i,1＝S_k的概率为公式(2)。

由公式(1)-(3)计算空间位置(i，j)处S_k(k＝1,n)的发生概率，将这些概率顺序排列，组成一个长度为1的数轴；利用计算机产生随机数，随机数落在数轴的那个区间中，该处就是生成该区间对应的那类基本单元。

实际工作中，对多数不能明确给出横向转移概率矩阵的情况，可用将纵向转移概率矩阵中的对角元素放大一定的倍数来近似横向转移概率矩阵，这个倍数大小恰好反映了横向均匀的程度。

上述算法将产生的是一个由不同类型岩相的“层”或“纹层”单元构成的二维模型；而对每类基本单元，根据其岩相的矿物组分，由随机过程生成分辨率为0.001×0.001m的非均匀构造，每个0.001×0.001m的小尺度介质对应一类具体的矿物，如方解石、粘土、石英或TOC。这样就可以得到一个0.001×0.001m分辨率的、由基本矿物构成的初始大尺度储层模型。

生成基本单元的随机过程是由其岩相的矿物组分控制，根据矿物组分的百分比组成一个长度为1的数轴；对基本单元中的随机一个离散点(按0.001×0.001m离散)，利用计算机产生随机数，随机数落在数轴的那个区间中，该点就对应那个矿物；基本单元中某一点确定一类矿物后，使用的数轴也要相应调整(该类矿物的百分比减少)。TOC的分布还要考虑粘土的分布和“层”或“纹层”构造的情况。

生成初始大尺度储层模型后，还需进一步在模型上随机添加裂缝和孔隙。孔隙的添加将依据岩芯样品统计得到的不同层位、不同岩相的孔隙率，对模型中的离散点赋上孔隙，对TOC的离散点，其孔隙对应有机孔。裂缝总量的多少由统计得到的裂缝密度决定。垂直缝的方向将由主应力方向控制。

完成上述过程，就可得到包含孔隙、裂缝、有机孔的分辨率为0.001×0.001m的大尺度储层模型，该模型既考虑了结构特点(层与纹层结构)，也考虑了矿物组份比例(粘土、方解石、石英和TOC)，还包括了裂缝和孔隙。

步骤102，基于大尺度储层模型，建立小尺度网格，获取等效弹性参数模型；在一个示例中，获取等效弹性参数模型包括：基于大尺度储层模型，获取储层的矿物类型；根据自洽近似理论，计算每一种矿物类型的弹性模量；根据Gassmann理论，实现孔隙的流体替换，计算流体的弹性模量；根据每一种矿物类型的弹性模量与流体的弹性模量，获取储层的裂缝数据，进而获取等效弹性参数模型。

在一个示例中，自洽近似理论为：

其中，表示矿物类型的等效体积模量，表示矿物类型的等效剪切模量，N代表共有N种矿物类型，x_i是第i种矿物的体积比，K_i表示第i种矿物的等效体积模量，μ_i表示第i种矿物的等效剪切模量，P^*i和Q^*i是第i种矿物在模量为和的背景介质中的几何系数。

具体地，岩石物理的等效介质理论研究表明，等效弹性参数不仅决定于介质的矿物成分、孔隙和裂缝，还受矿物组分自身的几何形状影响。上一节给出了用矿物成分表达的地震波长尺度的数值模型，但没有给出各个矿物组分的几何形状和孔隙的形状；实际上，即使给出矿物组分的几何形状，若在弹性波模拟的数值离散中考虑这些形状，将需要远小于0.001×0.001m尺度的计算网格。就针对地震波长尺度的模型进行数值模拟而言，这样的网格尺度需要太多的计算资源和计算时间。为避免在生成数值模型时考虑各矿物组分的几何形状的困难，也为大幅减少数值模拟的计算规模，利用等效介质理论，发展了计算网格尺度的地球物理参数等效建模方法，即根据已有的计算资源，确定尽可能小的计算网格尺寸(远小于频散要求)；利用测井速度估计矿物的几何形状系数；针对每个计算网格包括的矿物成分、孔隙、裂缝、有机孔、孔隙流体(或气体)和得到的几何形状系数，生成该计算网格内的等效拉梅系数和密度。

矿物的几何形状系数估计的主要步骤是：针对数值模型的井旁部分，由等效介质方法得到测井尺度的纵波速度和横波速度，根据等效速度与实测速度的差异，调整几何形状系数，两者接近时的几何参数就可认为是实际的几何形状系数。

根据TOC中有机孔的孔隙率和含流体或气体情况，利用公式(4)的自洽近似(SCA)理论，求得TOC的等效弹性参数；根据计算网格内各类矿物成分，如TOC、方解石、石英、粘土的弹性参数，将孔隙作为真空，再次利用公式(4)的自洽近似理论求得等效的干岩的压缩和剪切模量；若孔隙中充满流体，可根据Brown-Korringa广义Gassmann理论(Brown andKorringa，1975)，实现孔隙的流体替换，得到网格内所含矿物的流体饱和的等效体积模量和等效剪切模量；根据得到的等效体积模量和等效剪切模量，根据各个计算网格中包含裂缝的数量，由Hudson理论(Hudson，1981)得到含有裂缝的等效弹性参数。

针对每个计算网格完成上述步骤，就可完成计算网格尺度的等效弹性参数模型，这一过程考虑了不同矿物组分、有机孔、孔隙和裂缝等介质特征。

等效化处理仅限于各个计算网格内，因此可保留页岩储层的层状非均匀特征；而常规的基于整个地质模型的等效方法均假定各类矿物成分、孔隙、裂缝是均匀分布的，不能有效地考虑储层的非均匀性。

步骤103，基于等效弹性参数模型，获取储层的各向异性参数。

在一个示例中，各向异性参数为：

其中，δ、ε表示各向异性参数，(v_p)_v表示储层垂直方向的纵波速度，v_p(θ)表示储层沿角度θ方向的纵波速度。

具体地，通过模拟平面纵波和横波在页岩储层中的传播，来提取储层的纵波速度、横波速度和各向异性参数。等效弹性参数是频率相关的，为获得服务于地震勘探的岩石物理参数，采用与实际地震勘探相同主频的震源。

通过生成一个100米厚、1000米宽的地质模型，设计较宽的模型是为了方便求取各向异性参数，通过对各计算网格进行等效计算，得到一个100米厚、1000米宽的非均匀弹性介质模型。通过拾取通过储层顶底两个时刻的波场快照上的波前，来计算储层的物理参数。因此将储层模型放置到一个均匀背景中，在模型的均匀介质部位设置平面波震源，应用弹性各向异性介质中弹性地震波模拟方法，模拟平面波传播；纪录两个时刻的波场快照，分别拾取两个波场快照中平面波波峰的位置，即可求得储层的弹性参数。若设置水平的平面纵波震源，在t₁和t₂两个时刻纪录波场快照，在两个波场快照上拾取得到平面波峰值的坐标为x₁和x₂，可得到储层的垂直方向纵波速度为公式(6)：

其中，L是重叠储层的厚度，v_p0是均匀背景的纵波速度；为减少在储层界面处的反射，可令1/v_p0等于储层中各单元纵波速度倒数的平均值。

若设置与水平方向夹角为θ的平面纵波震源，可求得储层的沿角度θ方向的纵波速度为公式(7)：

其中，h₂-h₁是两个波前面沿垂直波前面方向的距离。

若设置水平的平面横波震源，在t₁和t₂两个时刻纪录波场快照，在两个波场快照上拾取得到平面波峰值的坐标为x₁和x₂，可得到储层的垂直方向横波速度(v_s)_v。

利用已求得的垂直、平行、斜入射的纵波速度，可利用公式(5)的VTI介质相速度公式，由最小二乘方法求得储层的各向异性参数。

本方法利用储层的统计数据，获得储层岩石物理弹性参数的变化规律，考虑实际储层的非均匀特征，可得到岩芯测试难以求取的与尺寸效应高度相关的弹性参数，也避免了求取弹性参数变化规律时获取不同地质特征岩芯的困难，为计算岩石物理面临的大尺度地质建模和计算能力限制问题提供了有效的解决方案。

应用示例

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

图2示出了根据本发明的一个实施例的储层模型的示意图。将100米后的页岩储层置于一个均匀的弹性介质中，爆炸源在储层上方800米；在与爆炸源相同的深度位置，设置检波器，接收储层的反射波，储层模型如图2所示。

图3a和图3b分别示出了根据本发明的一个实施例的3％TOC含量储层的垂直和水平速度分量炮记录的示意图。通过由计算岩石物理方法得到的弹性参数和各向异性参数(v_p,v_s,ε,δ,ρ)赋给图2模型中的储层，数值模拟得到纵波震源入射时的地震响应，如图3a和图3b所示。

根据某页岩储层岩芯分析结果，给出一个页岩储层数值建模例子。根据岩芯分析，确定3种不同岩相的“层”基本单元：含泥质灰岩层、泥质灰岩层、灰质泥岩层，其尺度定义为0.1×10m；确定2种不同岩相的“纹层”单元：灰质纹层、泥质纹层，其尺度定义为0.001×10m；共5个基本单元。储层的岩相构成、矿物占比如表1，“层”与“纹层”占比为4：1。

表1

图4示出了根据本发明的一个实施例的页岩储层生成的数值模型的纹层状部分的局部示意图。

图5示出了根据本发明的一个实施例的页岩储层生成的数值模型的层状部分的局部示意图。

图6示出了根据本发明的一个实施例的提取弹性参数的弹性介质纵波速度模型的示意图。

图7a和图7b分别示出了根据本发明的一个实施例的水平平面纵波入射0.1s和0.25s时刻的波场快照的示意图。

图8a和图8b分别示出了根据本发明的一个实施例的斜入射平面纵波0.1s和0.29s时刻的波场快照的示意图。

图9a和图9b分别示出了根据本发明的一个实施例的水平平面横波入射0.2s和0.38s时刻的波场快照的示意图。

然后，在页岩储层数值建模基础上，建立小尺度网格(计算网格尺度)，对网格内包括的非均匀矿物，考虑各种孔隙结构，利用等效介质理论生成等效的物理参数(纵波速度、横波速度和密度)。在小尺度网格上进行等效，充分考虑到了实际页岩的纵向强非均匀性。在地质与地球物理参数建模的基础上，根据物理参数和小尺度网格，进行不同角度平面地震波传播模拟，获取波长尺度的页岩模型的岩石物理参数。为提高速度拾取的精度，需建立远大于波长的模型，生成的页岩模型重复铺设，放置到一个均匀介质中。首先生成一个100米厚、1000米宽的地质模型(设计较宽的模型是为了方便求取各向异性参数)，通过对各计算网格进行等效计算，得到一个100米厚、1000米宽的非均匀弹性介质模型；较大的模型可减少拾取地震波波峰的误差，因此通过重叠该模型，可得到一个更大的，如400×1000m的模型。通过拾取两个时刻的波场快照上的波前，来计算储层的物理参数，因此将储层模型放置到一个均匀背景中，如图6所示。在均匀介质中，既容易激发平面波，也容易准确拾取平面波的波峰。图6中的储层模型，就是由图4和图5通过等效算法得到的，其计算网格的尺寸是0.02×0.02m。图7a和图7b、图8a和图8b、图9a和图9b分别展示了水平平面纵波入射、斜入射平面纵波以及水平平面横波入射不用时刻的波场快照。

以表1所示的胜利罗家的页岩储层为模板，生成了一组TOC含量从3％变化到21％的系列储层模型，系列储层模型的矿物占比、孔隙度、裂缝密度、纹层和层的比例见表2，地球物理参数等效建模所使用的各类矿物的岩石物理参数如表3所示。

表2

表3

图10示出了根据本发明的一个实施例的储层垂直方向纵波速度和横波速度随TOC含量变化曲线图，其中，实线代表纵波速度，虚线代表横波速度，从图中可看出，随TOC含量的增加，纵波速度迅速降低，横波速度也降低，但降低的幅度少于纵波速度。

图11示出了根据本发明的一个实施例的储层各向异性参数ε参数和δ参数随TOC含量变化曲线图，其中，实线代表ε参数，虚线代表δ参数，可见随TOC含量的增加，ε参数和δ参数也快速增加，这表明TOC含量增加导致储层的各向异性迅速增加。

综上所述，本发明利用储层的统计数据，获得储层岩石物理弹性参数的变化规律，考虑实际储层的非均匀特征，可得到岩芯测试难以求取的与尺寸效应高度相关的弹性参数，也避免了求取弹性参数变化规律时获取不同地质特征岩芯的困难，为计算岩石物理面临的大尺度地质建模和计算能力限制问题提供了有效的解决方案。

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

根据本发明的实施例，提供了一种基于地震波模拟的岩石物理参数计算系统，可以包括：存储器，存储有计算机可执行指令；处理器，处理器运行存储器中的计算机可执行指令，执行以下步骤：根据岩心数据，建立初始大尺度储层模型，进而基于裂缝数据与孔隙度数据，获取大尺度储层模型；基于大尺度储层模型，建立小尺度网格，获取等效弹性参数模型；基于等效弹性参数模型，获取储层的各向异性参数。

在一个示例中，建立初始大尺度储层模型包括：设定单元类型，进而确定每个位置的单元类型的概率；对储层的底层进行初始化处理，确定底层的单元类型；基于每个位置的单元类型的概率与底层的单元类型，建立初始大尺度储层模型。

在一个示例中，获取弹性参数等效模型包括：基于大尺度储层模型，获取储层的矿物类型；根据自洽近似理论，计算每一种矿物类型的弹性模量；根据Gassmann理论，实现孔隙的流体替换，计算流体的弹性模量；根据每一种矿物类型的弹性模量与流体的弹性模量，获取储层的裂缝数据，进而获取弹性参数等效模型。

本发明利用储层的统计数据，获得储层岩石物理弹性参数的变化规律，考虑实际储层的非均匀特征，可得到岩芯测试难以求取的与尺寸效应高度相关的弹性参数，也避免了求取弹性参数变化规律时获取不同地质特征岩芯的困难，为计算岩石物理面临的大尺度地质建模和计算能力限制问题提供了有效的解决方案。

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

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。

Claims (10)

1.一种基于地震波模拟的岩石物理参数计算方法，包括：

根据岩心数据，建立初始大尺度储层模型，进而基于裂缝数据与孔隙度数据，获取大尺度储层模型；

基于所述大尺度储层模型，建立小尺度网格，获取等效弹性参数模型；

基于所述等效弹性参数模型，获取储层的各向异性参数。

2.根据权利要求1所述的基于地震波模拟的岩石物理参数计算方法，其中，建立所述初始大尺度储层模型包括：

设定单元类型，进而确定每个位置的单元类型的概率；

对所述储层的底层进行初始化处理，确定所述底层的单元类型；

基于所述每个位置的单元类型的概率与所述底层的单元类型，建立所述初始大尺度储层模型。

3.根据权利要求2所述的基于地震波模拟的岩石物理参数计算方法，其中，每个位置的单元类型的概率，即Z_i,j＝S_k的概率为：

其中，p_k表示每个位置的单元类型的概率，i表示二维空间的纵向变化，j表示二维空间的深度变化，Z_i,j表示(i，j)处的单元类型，S_k(k＝1,n)代表定义的单元类型，n为定义的单元类型总数，f为定义的单元类型编号，p_mk是纵向转移概率矩阵中(m，k)处的元素，r_lk是横向转移概率矩阵中(l，k)处的元素，N是横向的单元总数，代表模型的横向大小；l源于Z_i-1,j＝S_l，m源于Z_i,j-1＝S_m，q源于Z_N,j＝S_q，p_mf为纵向转移概率矩阵中(m，f)处的元素，r_kq、r_lf、r_fq分别为横向转移概率矩阵中(k，q)、(l，f)、(f，q)处的元素。

4.根据权利要求2所述的基于地震波模拟的岩石物理参数计算方法，其中，所述底层的单元类型的概率，即Z_i,1＝S_k的概率为：

其中，r_lq为横向转移概率矩阵中(l，q)处的元素。

5.根据权利要求1所述的基于地震波模拟的岩石物理参数计算方法，其中，获取所述弹性参数等效模型包括：

基于所述大尺度储层模型，获取所述储层的矿物类型；

根据自洽近似理论，计算每一种矿物类型的弹性模量；

根据Gassmann理论，实现孔隙的流体替换，计算流体的弹性模量；

根据所述每一种矿物类型的弹性模量与所述流体的弹性模量，获取所述储层的裂缝数据，进而获取所述弹性参数等效模型。

6.根据权利要求5所述的基于地震波模拟的岩石物理参数计算方法，其中，所述自洽近似理论为：

其中，表示矿物类型的等效体积模量，表示矿物类型的等效剪切模量，N代表共有N种矿物类型，x_i是第i种矿物的体积比，K_i表示第i种矿物的等效体积模量，μ_i表示第i种矿物的等效剪切模量，P^*i和Q^*i是第i种矿物在模量为和的背景介质中的几何系数。

7.根据权利要求1所述的基于地震波模拟的岩石物理参数计算方法，其中，所述各向异性参数为：

其中，δ、ε表示各向异性参数，(v_p)_v表示储层垂直方向的纵波速度，v_p(θ)表示储层沿角度θ方向的纵波速度。

8.一种基于地震波模拟的岩石物理参数计算系统，其特征在于，该系统包括：

存储器，存储有计算机可执行指令；

处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：

根据岩心数据，建立初始大尺度储层模型，进而基于裂缝数据与孔隙度数据，获取大尺度储层模型；

基于所述大尺度储层模型，建立小尺度网格，获取等效弹性参数模型；

基于所述等效弹性参数模型，获取储层的各向异性参数。

9.根据权利要求8所述的基于地震波模拟的岩石物理参数计算系统，其中，建立所述初始大尺度储层模型包括：

设定单元类型，进而确定每个位置的单元类型的概率；

对所述储层的底层进行初始化处理，确定所述底层的单元类型；

基于所述每个位置的单元类型的概率与所述底层的单元类型，建立所述初始大尺度储层模型。

10.根据权利要求8所述的基于地震波模拟的岩石物理参数计算系统，其中，取所述弹性参数等效模型包括：

基于所述大尺度储层模型，获取所述储层的矿物类型；

根据自洽近似理论，计算每一种矿物类型的弹性模量；

根据Gassmann理论，实现孔隙的流体替换，计算流体的弹性模量；

根据所述每一种矿物类型的弹性模量与所述流体的弹性模量，获取所述储层的裂缝数据，进而获取所述弹性参数等效模型。