CN105242307B - 复杂碳酸盐储层地震孔隙度获取方法及装置 - Google Patents

Abstract

本发明提供了复杂碳酸盐储层地震孔隙度获取方法及装置，涉及石油勘探领域。本发明所提供的方法通过先获取目标地震固结系数、目标地震临界孔隙度、地震横波速度、地震饱含流体岩石密度和岩石基质剪切模量来计算目标地域的地震孔隙度，不仅考虑了临界孔隙度现象，同时考虑了孔隙结构对地震波速度的影响，因此求取孔隙度物理意义明确，对于具有复杂孔隙结构的碳酸盐岩储层具有很好的适用性，与其它地震孔隙度预测方法相比更符合实际地质情况。并且利用经过钻井开采所得出的实际产液量进行标定，能够直观看出使用该方法所计算的地震孔隙度结果能够较为准确的反映碳酸盐岩储层的储集空间大小。

Description

复杂碳酸盐储层地震孔隙度获取方法及装置

技术领域

本发明涉及石油勘探领域，具体而言，涉及复杂碳酸盐储层地震孔隙度获取方法及装置。

背景技术

随着工业技术的发展，人们对能源的需求量越来越大，常见的能源有电能、风能、化石能源等。日常生活中，如汽车中所使用的石油，即为化石能源。不同种类能源的来源是不相同的，如电能通常是靠其他形式的能量转化得来，如风力发电、水力发电均是将动能转化为电能。与电能不同的是，化石能源是通过动植物深埋一定时间之后所转化的，并且需要通过勘探和挖掘来获得。

随着化石能源的开发，当前社会越来越多的使用了化石能源。为了供给充分的化石能源，高效开采化石能源的工作越发受到重视。化石能源的获取可以分为两个步骤，第一个步骤是勘探，来判断当地是否存在化石能源，或者说判断当地化石能源的富集程度是否足够；如果第一个步骤判断为是的话，就会执行第二个步骤，也就是通过开采的方式来得到化石能源。由此可见，如何准确的探明当地化石能源的储量(即富集程度)，是进行开采的前提条件。

描述化石能源储量的一个主要参数是孔隙度，该参数描述了储集空间的大小，因此，如何准确的确定指定地域的孔隙度，便成为了确定该地域富集程度的首要工作。

针对我国的具体情况来看，深层复杂碳酸盐具有巨大油气勘探潜力，该储层是我国重要的油气资源战略接替区。复杂碳酸盐岩储层年代古老、埋藏深度大，由于受到长期的深埋压实与成岩胶结作用，原生孔隙已消失殆尽。后期成岩与构造作用改造强烈，储集空间以次生溶蚀孔洞和裂缝为主，具有十分复杂的孔隙结构和流体分布特点。各缝洞系统连通关系复杂，无统一的油水界面，油水分布不受现今构造高低的限制，同一套地层中可能高部位含水，低部位出油气，普遍存在油、水同出或气、水同出的现象。油气分布受储集体发育程度和连通性、断裂系统分布、油气充注方向与程度、成藏时间及成藏期构造格局等共同控制。依靠叠后属性和叠后反演等常规手段很难对储层中流体的充填性质及富集程度做出准确判断，致使在实际勘探开发生产中储量动用率、采出程度、采油速度等关键指标都非常低，远远达不到合理的工业指标。在实际勘探中高产，而不稳产现象十分普遍，甚至是许多井虽然钻遇到了很好的储层，但却被泥质或是水充填，导致勘探失败。随着深层复杂碳酸盐岩储层勘探的不断深入，逐渐形成了将高效井作为主要勘探和开发目标的理念。高效井的划分将石油价格和开发成本作为主要考核指标，已开发的生产井，按照一定的计算标准，如果累积油气产量以国际石油价格计算超过规定数值时，即称为高效井。高效井一般具有油气潜力大、经济效益高、容易进行开发作业等特点，但是如果要准确识别高效井，必须建立在定量评价储层内油气富集程度的基础上。而油气富集程度的确定又是建立在准确获取孔隙度的基础上。

相关技术中，孔隙度的获取方法归纳起来主要有实验室岩心测定、测井资料解释和地震资料预测等。由于碳酸盐岩储层一般具有复杂的孔隙结构，尤其是针对次生碳酸盐岩储层来说，有效的储集空间为次生的溶蚀孔、溶蚀洞、溶蚀缝和构造缝等，储层表现为极强的非均质性特征，因此对于这种复杂的碳酸盐岩储层而言，若要在横向上准确预测孔隙度的发育情况，就必须借助地震资料，并且要考虑复杂孔隙结构的影响。

目前存在多种地震孔隙度预测的方法，主要包括：以Wyllie时间平均方程(1956)和改进型的Raymer公式(1980)为基础，由地震纵波速度求取孔隙度的方法；通过建立地震反演结果与孔隙度的统计关系，利用纵波阻抗、密度、泊松比等参数和孔隙度的交汇图(Doyen，1988，Anderson，1996，Pramalik，2004)来求取孔隙度。以上两种方法具有计算简单、方便实用的优点，在地质条件简单、横向变化不大的地区应用效果较为理想，但很难应用于储层孔隙结构复杂、横向变化快的碳酸盐岩地区。

综上，相关技术中，有多种确定孔隙度的方法，但是，这些方法均是建立正在各自不同的假设基础上，在实际应用中均存在各自的局限性，计算得到的结果并不能够准确的反映实际的孔隙度。

发明内容

本发明的目的在于提供复杂碳酸盐储层地震孔隙度获取方法及装置，以提高获取碳酸盐岩储层孔隙度的准确度。

第一方面，本发明实施例提供了复杂碳酸盐储层地震孔隙度获取方法，包括：

通过叠前反演法获取目标地域的地震纵波速度、地震横波速度和地震饱含流体岩石密度；

通过统计测井资料，获取目标地域的岩石基质体积模量和岩石基质剪切模量；

根据测井资料、地层压力、温度、地层水的矿化度因素，获取目标地域的流体体积模量；

根据目标地域的地震纵波速度、地震横波速度、地震饱含流体岩石密度，以及测井资料统计的岩石基质体积模量、岩石基质剪切模量、流体体积模量计算目标地域的目标地震临界孔隙度和目标地震固结系数；

根据目标地震固结系数、目标地震临界孔隙度、地震横波速度、地震饱含流体岩石密度和测井资料统计的岩石基质剪切模量计算目标地域的地震孔隙度。

结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，获取目标地域的目标地震临界孔隙度包括：

分别计算目标地域中每个测井的参考测井临界孔隙度；

将每个测井的参考测井临界孔隙度与每个测井的叠前弹性参数进行交汇，以确定多个叠前弹性参数中，目标叠前弹性参数与参考测井临界孔隙度的函数关系，所述目标叠前弹性参数与参考测井临界孔隙度的拟合程度要符合预设的要求；

根据目标叠前弹性参数与参考测井临界孔隙度的函数关系，利用叠前地震反演的弹性参数计算目标地震临界孔隙度。

结合第一方面，本发明实施例提供了第一方面的第二种可能的实施方式，其中，获取目标地域的目标地震固结系数包括：

使用如下公式计算目标地域的目标地震固结系数， $\sqrt{\frac{k_m\left(\frac{{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{sei}\right)}^2+{\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei}{\left(1+{\mathrm{\α}}^{sei}\right)}^2}{{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{sei}\right)}^2+{\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei}{\left({\mathrm{\α}}^{sei}\right)}^2}\right)}{\left(1-{\mathrm{\φ}}^{sei}+{\mathrm{\φ}}^{sei}\frac{K_m}{K_f}\right){\mathrm{\ρ}}^{sei}-\frac{{\left(V_s^{sei}\right)}^2\left({\mathrm{\ρ}}^{sei}\right)\left(2{\left({\mathrm{\α}}^{sei}\right)}^2+2{\mathrm{\α}}^{sei}+1\right)}{{\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei}{\left({\mathrm{\α}}^{sei}\right)}^2+{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{sei}\right)}^2}}+\frac{k_m{\left(V_s^{sei}\right)}^2\left(2{\left({\mathrm{\α}}^{sei}\right)}^2+2{\mathrm{\α}}^{sei}+1\right)}{{\left({\mathrm{\α}}^{sei}\right)}^2{\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei}+{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{sei}\right)}^2}+\frac{4}{3}{\left(V_s^{sei}\right)}^2}-V_p^{sei}=0$ 其中， ${\mathrm{\φ}}^{sei}=\frac{{\mathrm{\φ}}_c^{sei}(1+{\mathrm{\α}}^{sei})({\mathrm{\μ}}_m-{\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei})}{{\mathrm{\α}}^{sei}(1+2{\mathrm{\α}}^{sei}){\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei}+{\mathrm{\μ}}_m(1+{\mathrm{\α}}^{sei})},{\mathrm{\φ}}_c^{sei}$ 为计算的目标地震临界孔隙度，ρ^sei分别为地震纵波速度、地震横波速度和地震饱含流体岩石密度，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^sei为需要求解的目标地震固结系数。并且，使用如下方法计算目标地震固结系数，该公式当中ρ^sei、K_m、K_d、K_f、μ_m为已知量，以目标地震固结系数α^sei为未知量，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出目标地震固结系数。

结合第一方面，本发明实施例提供了第一方面的第三种可能的实施方式，其中，根据目标地震固结系数、目标地震临界孔隙度、地震横波速度、地震饱含流体岩石密度和岩石基质剪切模量计算目标地域的地震孔隙度包括：

使用如下公式计算目标地域的地震孔隙度，其中，φ^sei为地震孔隙度，α^sei为目标地震固结系数，为目标地震临界孔隙度，为地震横波速度，ρ^sei为地震饱含流体岩石密度，μ_m为测井资料统计的岩石基质剪切模量。

结合第一方面，本发明实施例提供了第一方面的第四种可能的实施方式，其中，步骤分别计算目标地域中每个测井的参考测井临界孔隙度包括：

使用如下公式计算指定测井的参考测井临界孔隙度， ${\mathrm{\φ}}_c^{\log }=\frac{{\mathrm{\φ}}^{log}\[{\mathrm{\α}}^{log}(1+2{\mathrm{\α}}^{log}){\left(V_s^{log}\right)}^2{\mathrm{\ρ}}^{log}+{\mathrm{\μ}}_m(1+{\mathrm{\α}}^{log})\]}{(1+{\mathrm{\α}}^{log})({\mathrm{\μ}}_m-{\left(V_s^{1og}\right)}^2{\mathrm{\ρ}}^{log})},$ 其中，为指定测井的参考测井临界孔隙度，α^log为指定测井的参考测井固结系数，φ^log为指定测井的孔隙度测井曲线数据，为指定测井的横波速度测井曲线数据，ρ^log为指定测井的饱含流体岩石密度测井曲线数据，μ_m为测井资料统计的岩石基质剪切模量。

结合第一方面，本发明实施例提供了第一方面的第五种可能的实施方式，其中，步骤指定测井的参考测井固结系数：

使用如下公式计算指定测井的参考测井固结系数， $\sqrt{\frac{k_m\left(\frac{{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{\log }\right)}^2+{\left(V_s^{\log }\right)}^2{\mathrm{\ρ}}^{\log }{\left(1+{\mathrm{\α}}^{\log }\right)}^2}{{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{\log }\right)}^2+{\left(V_s^{\log }\right)}^2{\mathrm{\ρ}}^{\log }{\left({\mathrm{\α}}^{\log }\right)}^2}\right)}{\left(1-{\mathrm{\φ}}^{\log }+{\mathrm{\φ}}^{\log }\frac{K_m}{K_f}\right){\mathrm{\ρ}}^{\log }-\frac{{\left(V_s^{\log }\right)}^2\left({\mathrm{\ρ}}^{\log }\right)\left(2{\left({\mathrm{\α}}^{\log }\right)}^2+2{\mathrm{\α}}^{\log }+1\right)}{{\left(V_s^{\log }\right)}^2{\mathrm{\ρ}}^{\log }{\left({\mathrm{\α}}^{\log }\right)}^2+{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{\log }\right)}^2}}+\frac{k_m{\left(V_s^{\log }\right)}^2\left(2{\left({\mathrm{\α}}^{\log }\right)}^2+2{\mathrm{\α}}^{\log }+1\right)}{{\left({\mathrm{\α}}^{\log }\right)}^2{\left(V_s^{\log }\right)}^2{\mathrm{\ρ}}^{\log }+{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{\log }\right)}^2}+\frac{4}{3}{\left(V_s^{\log }\right)}^2}-V_p^{\log }=0$ ，其中，ρ^log分别为指定测井的纵波速度、横波速度和饱含流体岩石密度的测井曲线数据，φ^log为指定测井的孔隙度测井曲线数据，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^log为需要求解的指定测井的参考测井固结系数。并且，使用如下方法计算参考测井固结系数，该公式当中ρ^log、φ^log、K_m、K_d、K_f、μ_m为已知量，以指定测井的参考测井固结系数α^log为未知量，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出指定测井的参考测井固结系数。

结合第一方面，本发明实施例提供了第一方面的第六种可能的实施方式，其中，多个叠前弹性参数包括纵波阻抗、横波阻抗、纵横波速度比和密度，目标叠前弹性参数优选为与参考测井临界孔隙度拟合程度符合预设要求的弹性参数；

第二方面，本发明实施例还提供了复杂碳酸盐储层地震孔隙度获取装置，包括：

第一获取模块，用于通过叠前反演方法获取目标地域的地震纵波速度、地震横波速度和地震饱含流体岩石密度；

第二获取模块，用于通过统计测井资料，获取目标地域的岩石基质体积模量和岩石基质剪切模量；

第三获取模块，用于根据测井资料、地层压力、温度、地层水的矿化度因素，获取目标地域的流体体积模量；

第一计算模块，用于根据目标地域的地震纵波速度、地震横波速度、地震饱含流体岩石密度，以及测井资料统计的岩石基质体积模量、岩石基质剪切模量、流体体积模量计算目标地域的目标地震临界孔隙度和目标地震固结系数；

第二计算模块，用于根据目标地震临界孔隙度、目标地震固结系数、地震横波速度、地震饱含流体岩石密度和测井资料统计的岩石基质剪切模量计算目标地域的地震孔隙度。

结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方式，其中，第一计算模块包括：

第一计算单元，用于使用如下公式计算每个测井的参考测井固结系数，

$\sqrt{\frac{k_m\left(\frac{{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{\log }\right)}^2+{\left(V_s^{\log }\right)}^2{\mathrm{\ρ}}^{\log }{\left(1+{\mathrm{\α}}^{\log }\right)}^2}{{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{\log }\right)}^2+{\left(V_s^{\log }\right)}^2{\mathrm{\ρ}}^{\log }{\left({\mathrm{\α}}^{\log }\right)}^2}\right)}{\left(1-{\mathrm{\φ}}^{\log }+{\mathrm{\φ}}^{\log }\frac{K_m}{K_f}\right){\mathrm{\ρ}}^{\log }-\frac{{\left(V_s^{\log }\right)}^2\left({\mathrm{\ρ}}^{\log }\right)\left(2{\left({\mathrm{\α}}^{\log }\right)}^2+2{\mathrm{\α}}^{\log }+1\right)}{{\left(V_s^{\log }\right)}^2{\mathrm{\ρ}}^{\log }{\left({\mathrm{\α}}^{\log }\right)}^2+{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{\log }\right)}^2}}+\frac{k_m{\left(V_s^{\log }\right)}^2\left(2{\left({\mathrm{\α}}^{\log }\right)}^2+2{\mathrm{\α}}^{\log }+1\right)}{{\left({\mathrm{\α}}^{\log }\right)}^2{\left(V_s^{\log }\right)}^2{\mathrm{\ρ}}^{\log }+{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{\log }\right)}^2}+\frac{4}{3}{\left(V_s^{\log }\right)}^2}-V_p^{\log }=0$ ，其中，ρ^log分别为指定测井的纵波速度、横波速度和饱含流体岩石密度的测井曲线数据，φ^log为指定测井的孔隙度测井曲线数据，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^log为需要求解的指定测井的参考测井固结系数。并且，使用如下方法计算指定测井的参考测井固结系数，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出指定测井的参考测井固结系数；

第二计算单元，用于使用如下公式计算每个测井的参考测井临界孔隙度， ${\mathrm{\φ}}_c^{\log }=\frac{{\mathrm{\φ}}^{log}\[{\mathrm{\α}}^{log}(1+2{\mathrm{\α}}^{log}){\left(V_s^{log}\right)}^2{\mathrm{\ρ}}^{log}+{\mathrm{\μ}}_m(1+{\mathrm{\α}}^{log})\]}{(1+{\mathrm{\α}}^{log})({\mathrm{\μ}}_m-{\left(V_s^{1og}\right)}^2{\mathrm{\ρ}}^{log})},$ 其中，为需要求解的指定测井的参考测井临界孔隙度，α^log为指定测井的参考测井固结系数，φ^log为指定测井的孔隙度测井曲线数据，为指定测井的横波速度测井曲线数据，ρ^log为指定测井的饱含流体岩石密度测井曲线数据，μ_m为测井资料统计的岩石基质剪切模量。

交汇模块，用于将每个测井的参考测井临界孔隙度与每个测井的叠前弹性参数进行交汇，以确定多个叠前弹性参数中，目标叠前弹性参数与参考测井临界孔隙度的函数关系，所述目标叠前弹性参数与参考测井临界孔隙度的拟合程度要符合预设的要求；

第三计算单元，用于根据目标叠前弹性参数与参考测井临界孔隙度的函数关系，利用叠前地震反演的目标叠前弹性参数计算目标地震临界孔隙度。

结合第二方面，本发明实施例提供了第二方面的第二种可能的实施方式，其中，第一计算模块还包括：

第四计算单元，用于使用如下公式计算目标地域的目标地震固结系数， $\sqrt{\frac{k_m\left(\frac{{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{sei}\right)}^2+{\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei}{\left(1+{\mathrm{\α}}^{sei}\right)}^2}{{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{sei}\right)}^2+{\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei}{\left({\mathrm{\α}}^{sei}\right)}^2}\right)}{\left(1-{\mathrm{\φ}}^{sei}+{\mathrm{\φ}}^{sei}\frac{K_m}{K_f}\right){\mathrm{\ρ}}^{sei}-\frac{{\left(V_s^{sei}\right)}^2\left({\mathrm{\ρ}}^{sei}\right)\left(2{\left({\mathrm{\α}}^{sei}\right)}^2+2{\mathrm{\α}}^{sei}+1\right)}{{\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei}{\left({\mathrm{\α}}^{sei}\right)}^2+{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{sei}\right)}^2}}+\frac{k_m{\left(V_s^{sei}\right)}^2\left(2{\left({\mathrm{\α}}^{sei}\right)}^2+2{\mathrm{\α}}^{sei}+1\right)}{{\left({\mathrm{\α}}^{sei}\right)}^2{\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei}+{\mathrm{\μ}}_m{\left(1+{\mathrm{\α}}^{sei}\right)}^2}+\frac{4}{3}{\left(V_s^{sei}\right)}^2}-V_p^{sei}=0$ 其中， ${\mathrm{\φ}}^{sei}=\frac{{\mathrm{\φ}}_c^{sei}(1+{\mathrm{\α}}^{sei})({\mathrm{\μ}}_m-{\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei})}{{\mathrm{\α}}^{sei}(1+2{\mathrm{\α}}^{sei}){\left(V_s^{sei}\right)}^2{\mathrm{\ρ}}^{sei}+{\mathrm{\μ}}_m(1+{\mathrm{\α}}^{sei})},{\mathrm{\φ}}_c^{sei}$ 为计算的目标地震临界孔隙度，分别为地震纵波速度、地震横波速度和地震饱含流体岩石密度，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^sei为需要求解的目标地震固结系数；

并且，使用如下方法求得目标地震固结系数，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出目标地震固结系数。

本发明实施例提供的复杂碳酸盐储层地震孔隙度获取方法，使用目标地震固结系数、目标地震临界孔隙度、地震横波速度、地震饱含流体岩石密度和岩石基质剪切模量来计算目标地域的地震孔隙度，不仅考虑了临界孔隙度现象，认为碳酸盐岩的临界孔隙度并不是固定不变的，而是受到地区岩性、沉积历史、成岩形式和岩溶改造诸多因素的影响而变化；同时通过自适应调整固结系数，考虑了孔隙结构对地震波速度的影响，因此求取孔隙度物理意义明确，对于具有复杂孔隙结构的碳酸盐岩储层具有很好的适用性，与其它地震孔隙度预测方法相比更符合实际地质情况。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法的细化流程图；

图2示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法中，塔里木盆地塔中地区某区块地震纵波速度反演连井剖面(图中测井曲线为纵波速度曲线)；

图3示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法中，塔里木盆地塔中地区某区块地震横波速度反演连井剖面(图中测井曲线为横波速度曲线)；

图4示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法中，塔里木盆地塔中地区某区块地震饱含流体岩石密度反演连井剖面(图中测井曲线为饱含流体岩石密度曲线)；

图5示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法中，基于测井资料的体积模量与孔隙度曲线交汇图；

图6示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法中，基于测井资料的剪切模量与孔隙度曲线交汇图；

图7示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法中，基于测井资料计算的参考测井临界孔隙度与各叠前弹性参数交汇图；

图8示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法中，塔里木盆地塔中地区某区块地震孔隙度计算结果连井剖面(图中测井曲线为孔隙度曲线)；

图9示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法中，塔里木盆地塔中地区某区块目的层地震孔隙度计算结果平面分布图；

图10示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法中，塔里木盆地塔中地区某区块各井地震预测孔隙度与累计产液量统计图；

图11示出了本发明实施例所提供的复杂碳酸盐储层地震孔隙度获取方法的基本流程图。

具体实施方式

下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

油气富集程度能够说明某个地域内的油气储量，如果该地区的油气储量超过预定的阈值后，才有必要对该地区的进行开采。可以认为，在开采之前，确定目标地域的油气储量是必须进行的一个步骤，而油气储量又是通过油气富集程度来表征的。因此，如何准确的确定油气富集程度便成为了一个重要的工作。

而精确的确定指定地域的油气富集程度又是依赖于准确的确定该地域的孔隙度。因此，如何准确的确定孔隙度便成为更为重要的一个环节。

相关技术中，已存在多种地震孔隙度预测的方法，主要包括：以Wyllie时间平均方程(1956)和改进型的Raymer公式(1980)为基础，由地震纵波速度求取孔隙度的方法；通过建立地震反演结果与孔隙度的统计关系，利用纵波阻抗、密度、泊松比等参数和孔隙度的交汇图(Doyen，1988，Anderson，1996，Pramalik，2004)来求取孔隙度。以上两种方法具有计算简单、方便实用的优点，在地质条件简单、横向变化不大的地区应用效果较为理想,但很难应用于储层孔隙结构复杂、横向变化快的碳酸盐岩地区。

Gassmann(1951)提出了饱和流体岩石弹性模量计算公式，从严格意义上建立了纵波速度与岩石固体颗粒的弹性模量，流体体积模量，干燥岩石体积模量及孔隙度、岩石密度的关系。由于该方程适用于充分低频的情况(孔隙流体处于静态)，所以它能很好地适用于地震资料(<100Hz)，具有非常重要的实用价值。

$V_p^2=\frac{1}{\mathrm{\&rho;}}\&lsqb;\frac{\left(K_m+K_d\right)}{K_m\left(1-\mathrm{\&phi;}-\frac{K_d}{K_m}+\mathrm{\&phi;}\frac{K_m}{K_f}\right)}+K_d+\frac{4}{3}{\mathrm{\&mu;}}_d\&rsqb;$ (式1)

${\mathrm{\&mu;}}_d=V_s^2\mathrm{\&rho;}$ (式2)

式中：V_p、V_s分别代表纵、横波速度，ρ为饱含流体的岩石密度，φ为孔隙度，K_m、K_d、K_f分别代表岩石基质等效体积模量、干岩石的体积模量和孔隙内流体的体积模量，μ_d代表干岩石的剪切模量。公式1和公式2为统一的岩石物理模型公式，即可用于测井资料的计算，也可用于地震资料的计算，纵波速度Vp、横波速度Vs、饱含流体的岩石密度ρ和孔隙度φ即可代表测井资料数据，也可代表地震预测参数，仅是输入参数不同而已，如φ可代表测井资料孔隙度φ^log，也可代表地震预测的孔隙度φ^sei，本专利其它公式下同。

通过上式可知，Gassmann方程中存在7个参数：V_p、V_s、ρ、K_m、K_f、K_d、φ，当已知或合理估计出其中前6个参数后，就可以计算岩石的孔隙度。实际上，利用Gassmann方程进行测井资料的孔隙度计算，只有理论意义，孔隙度数据可利用测井数据直接得到。而利用地震资料来估计孔隙度则具有重要的应用价值，因为地震资料覆盖面广，并有较高的横向分辨率。

基于Gassmann方程进行孔隙度计算是目前较为流行的一种方法，该方法的优点是获得了计算孔隙度的解析式，有坚实的理论基础，缺点是需要预先提供的参数很多，尤其是需要事先获取纵波速度、横波速度和密度等作为输入参数，并且密度反演具有较强的不稳定性，利用常规的反演方法很难准确获取密度参数；此外，常规利用Gassmann方程求取孔隙度的方法没有考虑孔隙结构，因此限制了该类方法在碳酸盐岩储层中的应用。

利用Gassmann方程研究饱和岩石的弹性性质时，需要预先提供干岩石的体积模量。目前干岩石体积模量获取的方法可以分为三类：一是岩石样品的实验室直接测量，利用测量得到的岩石骨架的纵波速度、横波速度和密度来计算获得，该方法获得的干岩石体积模量精度高，但成本昂贵；二是等效介质理论，就是在已知岩石各组成成分相对含量和弹性模量情况下，通过调整岩石物理模型中的自由参数(往往代表岩石中组成成分的几何特征，比如孔隙结构参数)来定量求取干岩石的体积模量，其中适用于碳酸盐岩的模型主要有模型(1974)、Xu-Payne模型(2009)和DEM-Gassmann模型(2012)等，但这些模型包括参数较多，很难应用于实际地震资料；三是经验模型，其中常见的经验模型是临界孔隙度模型(1992,1998)、Krief模型(1990)、Pride模型(2004)和临界孔隙度Pride模型(2012)等，这些经验模型在计算干岩石体积模量过程中简单实用，但都有各自不同的假设，在实际应用中均存在各自的局限性。

侯波等(2012)基于Pride模型，通过引入临界孔隙度提出了临界孔隙度Pride模型。该模型兼取临界孔隙度模型和Pride模型的优点，建立了如下临界孔隙度Pride模型：

$K_d=\frac{K_m(1-\mathrm{\&phi;}/{\mathrm{\&phi;}}_c)}{1+\mathrm{\&alpha;}\mathrm{\&phi;}/{\mathrm{\&phi;}}_c}$ (式3)

${\mathrm{\&mu;}}_d=\frac{{\mathrm{\&mu;}}_m(1-\mathrm{\&phi;}/{\mathrm{\&phi;}}_c)}{1+\mathrm{\&gamma;}\mathrm{\&alpha;}\mathrm{\&phi;}/{\mathrm{\&phi;}}_c},\mathrm{\&gamma;}=\frac{1+2\mathrm{\&alpha;}}{1+\mathrm{\&alpha;}}$ `(式4)

式中：K_m、K_d分别代表岩石基质等效体积模量和干岩石的体积模量，μ_m、μ_d分别代表岩石基质和干岩样的剪切模量，φ为孔隙度，φ_c为临界孔隙度，α为固结系数，γ为关于固结系数的关系经验关系式。公式3和公式4为统一的岩石物理模型公式，即可用于测井资料的计算，也可用于地震资料的计算，公式中的临界孔隙度φ_c和固结系数α即可代表测井资料计算参数，也可代表地震资料预测参数，仅是输入参数不同而已，如φ_c可代表参考测井临界孔隙度也可代表地震预测的临界孔隙度本专利其它公式下同。

该模型假设同一岩性岩石的临界孔隙度是固定不变的，并参考Mavko(1998)的观点，将碳酸盐岩临界孔隙度取值为60％。但是由于碳酸盐岩储层成岩形式多样，后期岩溶改造作用强烈，导致同一岩性的临界孔隙度并不是固定不变的，而将其假设为常数显然并不适用于碳酸盐岩复杂的地质条件，求得的孔隙度结果精度不高。

有鉴于此，本申请提供了复杂碳酸盐储层地震孔隙度获取方法，如图11所示，包括如下步骤：

S101，通过叠前反演方法获取目标地域的地震纵波速度、地震横波速度和地震饱含流体岩石密度；

S102，通过统计测井资料，获取目标地域的岩石基质体积模量和岩石基质剪切模量；

S103，根据测井资料、地层压力、温度、地层水的矿化度因素，获取目标地域的流体体积模量；

S104，根据目标地域的地震纵波速度、地震横波速度、地震饱含流体岩石密度，以及测井资料统计的岩石基质体积模量、岩石基质剪切模量、流体体积模量计算目标地域的目标地震临界孔隙度和目标地震固结系数；

S105，根据目标地震临界孔隙度、目标地震固结系数、地震横波速度、地震饱含流体岩石密度和测井资料统计的岩石基质剪切模量计算目标地域的地震孔隙度。

其中，计算目标地域地震孔隙度是依赖于获取准确的目标地震固结系数、目标地震临界孔隙度，下面，说明如何获得这两个参数。

由于目标地域中通常同时存在有多个测井，因此可以对目标地震临界孔隙度和目标地震固结系数的获取进行调整，即，步骤S104，计算目标地域的目标地震临界孔隙度和目标地震固结系数包括：

分别计算目标地域中每个测井的参考测井临界孔隙度；

将每个测井的参考测井临界孔隙度与每个测井的叠前弹性参数进行交汇，以确定多个叠前弹性参数中，目标叠前弹性参数与参考测井临界孔隙度的函数关系，所述目标叠前弹性参数与参考测井临界孔隙度的拟合程度要符合预设的要求；

根据目标叠前弹性参数与参考测井临界孔隙度的函数关系，利用叠前地震反演的目标叠前弹性参数计算目标地震临界孔隙度。

其中，步骤多个叠前弹性参数包括纵波阻抗、横波阻抗、纵横波速度比和密度，目标叠前弹性参数为优选的与参考测井临界孔隙度拟合程度符合预设要求的弹性参数，本发明实施例当中通过实际计算与对比分析，目标叠前弹性参数优选为地震纵波阻抗。

进一步，使用如下公式计算目标地域的目标地震固结系数， $\sqrt{\frac{k_m\left(\frac{{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{sei}\right)}^2+{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}{\left(1+{\mathrm{\&alpha;}}^{sei}\right)}^2}{{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{sei}\right)}^2+{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}{\left({\mathrm{\&alpha;}}^{sei}\right)}^2}\right)}{\left(1-{\mathrm{\&phi;}}^{sei}+{\mathrm{\&phi;}}^{sei}\frac{K_m}{K_f}\right){\mathrm{\&rho;}}^{sei}-\frac{{\left(V_s^{sei}\right)}^2\left({\mathrm{\&rho;}}^{sei}\right)\left(2{\left({\mathrm{\&alpha;}}^{sei}\right)}^2+2{\mathrm{\&alpha;}}^{sei}+1\right)}{{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}{\left({\mathrm{\&alpha;}}^{sei}\right)}^2+{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{sei}\right)}^2}}+\frac{k_m{\left(V_s^{sei}\right)}^2\left(2{\left({\mathrm{\&alpha;}}^{sei}\right)}^2+2{\mathrm{\&alpha;}}^{sei}+1\right)}{{\left({\mathrm{\&alpha;}}^{sei}\right)}^2{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}+{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{sei}\right)}^2}+\frac{4}{3}{\left(V_s^{sei}\right)}^2}-V_p^{sei}=0$ 其中， ${\mathrm{\&phi;}}^{sei}=\frac{{\mathrm{\&phi;}}_c^{sei}(1+{\mathrm{\&alpha;}}^{sei})({\mathrm{\&mu;}}_m-{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei})}{{\mathrm{\&alpha;}}^{sei}(1+2{\mathrm{\&alpha;}}^{sei}){\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}+{\mathrm{\&mu;}}_m(1+{\mathrm{\&alpha;}}^{sei})},{\mathrm{\&phi;}}_c^{sei}$ 为计算的目标地震临界孔隙度，ρ^sei分别为地震纵波速度、地震横波速度和地震饱含流体岩石密度，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^sei为需要求解的目标地震固结系数。该公式当中ρ^sei、K_m、K_d、K_f、μ_m为已知量，以目标地震固结系数α^sei为未知量，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出目标地震固结系数。

步骤S105，根据目标地震临界孔隙度、目标地震固结系数、地震横波速度、地震饱含流体岩石密度和测井资料统计的岩石基质剪切模量计算目标地域的地震孔隙度包括：

使用如下公式计算目标地域的地震孔隙度， ${\mathrm{\&phi;}}^{sei}=\frac{{\mathrm{\&phi;}}_c^{sei}(1+{\mathrm{\&alpha;}}^{sei})({\mathrm{\&mu;}}_m-{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei})}{{\mathrm{\&alpha;}}^{sei}(1+2{\mathrm{\&alpha;}}^{sei}){\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}+{\mathrm{\&mu;}}_m(1+{\mathrm{\&alpha;}}^{sei})},$ 其中，φ^sei为地震孔隙度，α^sei为目标地震固结系数，为目标地震临界孔隙度，V_s ^sei为地震横波速度，ρ^sei为地震饱含流体岩石密度，μ_m为测井资料统计的岩石基质剪切模量。

进一步，步骤分别计算目标地域中每个测井的参考测井临界孔隙度包括：

使用如下公式计算指定测井的参考测井临界孔隙度， ${\mathrm{\&phi;}}_c^{\log }=\frac{{\mathrm{\&phi;}}^{log}\&lsqb;{\mathrm{\&alpha;}}^{log}(1+2{\mathrm{\&alpha;}}^{log}){\left(V_s^{log}\right)}^2{\mathrm{\&rho;}}^{log}+{\mathrm{\&mu;}}_m(1+{\mathrm{\&alpha;}}^{log})\&rsqb;}{(1+{\mathrm{\&alpha;}}^{log})({\mathrm{\&mu;}}_m-{\left(V_s^{1og}\right)}^2{\mathrm{\&rho;}}^{log})},$ 其中，为指定测井的参考测井临界孔隙度，α^log为指定测井的参考测井固结系数，φ^log为指定测井的孔隙度测井曲线数据，为指定测井的横波速度测井曲线数据，ρ^log为指定测井的饱含流体岩石密度测井曲线数据，μ_m为测井资料统计的岩石基质剪切模量。

更进一步，对指定测井的参考测井固结系数的获取进行说明，还包括：

使用如下公式计算指定测井的参考测井固结系数， $\sqrt{\frac{k_m\left(\frac{{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{\log }\right)}^2+{\left(V_s^{\log }\right)}^2{\mathrm{\&rho;}}^{\log }{\left(1+{\mathrm{\&alpha;}}^{\log }\right)}^2}{{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{\log }\right)}^2+{\left(V_s^{\log }\right)}^2{\mathrm{\&rho;}}^{\log }{\left({\mathrm{\&alpha;}}^{\log }\right)}^2}\right)}{\left(1-{\mathrm{\&phi;}}^{\log }+{\mathrm{\&phi;}}^{\log }\frac{K_m}{K_f}\right){\mathrm{\&rho;}}^{\log }-\frac{{\left(V_s^{\log }\right)}^2\left({\mathrm{\&rho;}}^{\log }\right)\left(2{\left({\mathrm{\&alpha;}}^{\log }\right)}^2+2{\mathrm{\&alpha;}}^{\log }+1\right)}{{\left(V_s^{\log }\right)}^2{\mathrm{\&rho;}}^{\log }{\left({\mathrm{\&alpha;}}^{\log }\right)}^2+{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{\log }\right)}^2}}+\frac{k_m{\left(V_s^{\log }\right)}^2\left(2{\left({\mathrm{\&alpha;}}^{\log }\right)}^2+2{\mathrm{\&alpha;}}^{\log }+1\right)}{{\left({\mathrm{\&alpha;}}^{\log }\right)}^2{\left(V_s^{\log }\right)}^2{\mathrm{\&rho;}}^{\log }+{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{\log }\right)}^2}+\frac{4}{3}{\left(V_s^{\log }\right)}^2}-V_p^{\log }=0$ ，其中，ρ^log分别为指定测井的纵波速度、横波速度和饱含流体岩石密度的测井曲线数据，φ^log为指定测井的孔隙度测井曲线数据，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^log为需要求解的指定测井的参考测井固结系数。该公式当中ρ^log、φ^log、K_m、K_d、K_f、μ_m为已知量，以指定测井的参考测井固结系数α^log为未知量，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出指定测井的参考测井固结系数。

本发明提出了一种新的考虑复杂孔隙结构的碳酸盐岩储层地震孔隙度计算方法，推导了地震孔隙度计算公式，不仅考虑了临界孔隙度现象，同时考虑了孔隙结构对地震波速度的影响，因此求取孔隙度物理意义明确，对于具有复杂孔隙结构的碳酸盐岩储层具有很好的适用性，与其它地震孔隙度预测方法相比更符合实际地质情况。由于碳酸盐岩储层成岩形式多样，后期岩溶改造作用强烈，导致同一岩性的临界孔隙度并不是固定不变的，而目前的岩石物理模型和孔隙度预测方法当中为分析方便，均假设同一岩性岩石的临界孔隙度都是相同的，显然这一假设并不适用于碳酸盐岩复杂的地质条件。本发明认为碳酸盐岩的临界孔隙度并不是固定的常数，而是受到地区岩性、沉积历史、成岩形式和岩溶改造等诸多因素的影响而变化，更加符合客观地质实际，并通过碳酸盐岩地区实际测井资料分析，利用叠前弹性参数来拟合地震临界孔隙度参数。另外许多研究都已经证实在碳酸盐岩储层中孔隙结构对声波速度影响很大，孔隙结构是影响声波速度不容忽视的重要因素，尤其是针对以次生储集空间主导的复杂碳酸盐岩储层来说，计算孔隙度时必须考虑孔隙结构的影响。本发明认为固结系数与岩石孔隙形状有关，被认为是带有孔隙形态调节作用的参数，本发明运用混沌量子粒子群非线性优化算法自适应调整固结系数，因此适用于具有复杂孔隙结构的碳酸盐岩储层。

下面，对本申请的地震孔隙度计算方法中所使用到的公式进行说明：

首先将临界孔隙度Pride模型和Gassmann进行整理，将公式2代入公式4中进行整理，推导出公式5：

$\frac{\mathrm{\&phi;}}{{\mathrm{\&phi;}}_c}=\frac{(1+\mathrm{\&alpha;})({\mathrm{\&mu;}}_m-{V_s}^2\mathrm{\&rho;})}{\mathrm{\&alpha;}(1+2\mathrm{\&alpha;}){V_s}^2\mathrm{\&rho;}+{\mathrm{\&mu;}}_m(1+\mathrm{\&alpha;})}$ (式5)

式中：φ为孔隙度，α为目标固结系数，φ_c为目标临界孔隙度，μ_m为岩石基质剪切模量，V_s为横波速度，ρ为饱含流体的岩石密度。

对公式5进行整理，可推导出公式6，即本申请所提供的孔隙度计算公式：

$\mathrm{\&phi;}=\frac{{\mathrm{\&phi;}}_c(1+\mathrm{\&alpha;})({\mathrm{\&mu;}}_m-{V_s}^2\mathrm{\&rho;})}{\mathrm{\&alpha;}(1+2\mathrm{\&alpha;}){V_s}^2\mathrm{\&rho;}+{\mathrm{\&mu;}}_m(1+\mathrm{\&alpha;})}$ (式6)

式中：φ为孔隙度，α为目标固结系数，φ_c为目标临界孔隙度，μ_m为岩石基质剪切模量，V_s为横波速度，ρ为饱含流体的岩石密度。

其中，岩石基质剪切模量μ_m可以通过测井资料统计分析获得，横波速度V_s和饱含流体的岩石密度ρ可通过叠前反演获得。则孔隙度计算公式(式6)右侧中存在两个未知量，分别是目标固结系数α和目标临界孔隙度φ_c。

将公式5代入公式3中进行整理，推导出公式7：

$K_d=\frac{K_m{V_s}^2\mathrm{\&rho;}(2{\mathrm{\&alpha;}}^2+2\mathrm{\&alpha;}+1)}{{\mathrm{\&alpha;}}^2{V_s}^2\mathrm{\&rho;}+{\mathrm{\&mu;}}_m{(1+\mathrm{\&alpha;})}^2}$ (式7)

式中：K_m、K_d分别代表岩石基质等效体积模量和干岩石的体积模量，μ_m为岩石基质剪切模量，V_s为横波速度，ρ为饱含流体的岩石密度，α为目标固结系数。

将公式7代入Gassmann方程公式1中进行整理，推导出公式8，即本申请的固结系数的隐性计算公式：

$\sqrt{\frac{k_m{\left(\frac{{\mathrm{\&mu;}}_m{(1+\mathrm{\&alpha;})}^2+{V_s}^2\mathrm{\&rho;}{(1+\mathrm{\&alpha;})}^2}{{\mathrm{\&mu;}}_m{(1+\mathrm{\&alpha;})}^2+{V_s}^2{\mathrm{\&rho;\&alpha;}}^2}\right)}^2}{(1-\mathrm{\&phi;}+\mathrm{\&phi;}\frac{K_m}{K_f})\mathrm{\&rho;}-\frac{{V_s}^2{\mathrm{\&rho;}}^2(2{\mathrm{\&alpha;}}^2+2\mathrm{\&alpha;}+1)}{{V_s}^2{\mathrm{\&rho;\&alpha;}}^2+{\mathrm{\&mu;}}_m{(1+\mathrm{\&alpha;})}^2}}+\frac{k_m{V_s}^2(2{\mathrm{\&alpha;}}^2+2\mathrm{\&alpha;}+1)}{{\mathrm{\&alpha;}}^2{V_s}^2\mathrm{\&rho;}+{\mathrm{\&mu;}}_m{(1+\mathrm{\&alpha;})}^2}+\frac{4}{3}{V_s}^2}-V_p=0$ (式8)

式中：V_p、V_s分别代表纵、横波速度，ρ为饱含流体的岩石密度，K_m、K_d、K_f分别代表岩石基质的体积模量、干岩石的体积模量和孔隙内流体的体积模量，μ_m代表岩石基质的剪切模量，φ为孔隙度，α为目标固结系数。

为了研究方便和简化计算过程，目前的岩石物理模型当中均假设同一岩性岩石的临界孔隙度都是相同的，如碳酸盐岩临界孔隙度取值为60％，砂岩临界孔隙度取值为40％(Mavko,1998)。但是由于碳酸盐岩储层成岩形式多样，后期岩溶改造作用强烈，导致同一岩性的临界孔隙度并不是固定不变的，而将其假设为常数显然并不适用于碳酸盐岩复杂的地质条件。本发明认为碳酸盐岩的临界孔隙度并不是固定的常数，而是受到地区岩性、沉积历史、成岩形式和岩溶改造等诸多因素的影响而变化，更加符合客观地质实际，进而可以提高孔隙度的计算精度。

对公式5进行整理，可推导出公式9，即本发明的临界孔隙度计算公式：

${\mathrm{\&phi;}}_c=\frac{\mathrm{\&phi;}\&lsqb;\mathrm{\&alpha;}(1+2\mathrm{\&alpha;}){V_s}^2\mathrm{\&rho;}+{\mathrm{\&mu;}}_m(1+\mathrm{\&alpha;})\&rsqb;}{(1+\mathrm{\&alpha;})({\mathrm{\&mu;}}_m-{V_s}^2\mathrm{\&rho;})}$ (式9)

式中：φ为孔隙度，φ_c为临界孔隙度，α为固结系数，μ_m为岩石基质剪切模量，V_s为横波速度，ρ为饱含流体的岩石密度。

下面，以一个具体的实例还来说明本申请所提供的复杂碳酸盐储层地震孔隙度获取方法。该方法可以分为如下6个步骤，如图1所示。

步骤1：将研究区中，钻井的纵波速度测井曲线数据横波速度测井曲线数据密度测井曲线数据ρ^log、孔隙度测井曲线数据φ^log以及利用测井资料统计获得的岩石基质的体积模量K_m、干岩石的体积模量K_d、孔隙内流体的体积模量K_f和岩石基质的剪切模量μ_m输入到公式8中，则公式中固结系数α为唯一待求解未知量，可见公式8为非线性求解方程。认为固结系数与岩石孔隙形状有关，被认为是带有孔隙形态调节作用的参数，利用混沌量子粒子群非线性优化算法进行全局搜索，使公式8计算误差最小，当计算误差小于设定的阈值时，则输出指定测井的参考测井固结系数α^log；

步骤2：将研究区已钻井的横波速度测井曲线数据密度测井曲线数据ρ^log、孔隙度测井曲线数据φ^log、利用测井资料统计获得的岩石基质剪切模量μ_m以及步骤1计算的指定测井的参考测井固结系数α^log输入到临界孔隙度计算公式(式9)中，就可以计算出指定测井的参考测井临界孔隙度

步骤3：将步骤2计算的各指定测井的参考测井临界孔隙度与测井资料计算的多种叠前弹性参数(纵波阻抗I_p ^log、横波阻抗I_s ^log、纵横波速度比(V_p/V_s)^log、密度ρ^log)进行交汇，建立目标叠前弹性参数与参考测井临界孔隙度的函数关系，所述目标叠前弹性参数与参考测井临界孔隙度的拟合程度要符合预设的要求：。

φ_c ^log＝A·I^log+B (式10)

式中：A和B为利用实际测井资料拟合得到的系数，I^log为测井资料计算的目标叠前弹性参数。

在本发明所记载的实施例中，纵波阻抗测井曲线数据与参考测井临界孔隙度线性拟合关系最好，符合预设要求。

步骤4：将步骤3获得的目标叠前弹性参数与参考测井临界孔隙度的函数关系应用到地震数据当中，则得到公式11，进而利用叠前地震反演获得叠前弹性参数获得目标地震临界孔隙度；

φ_c ^sei＝A·I^sei+B (式11)

式中：为目标地震临界孔隙度，I^sei为地震资料计算的叠前弹性参数(本发明所记载的实施例中指地震纵波阻抗I_p ^sei)。

步骤5：将公式6代入到公式8中，以地震资料反演获得的纵波速度V_p ^sei、横波速度V_s ^sei、饱含流体岩石密度ρ^sei以及利用测井资料统计获得的岩石基质体积模量K_m、孔隙内流体体积模量K_f、岩石基质剪切模量μ_m、步骤4计算的目标地震临界孔隙度φ_c ^sei作为输入参数，则公式中目标地震固结系数α^sei为唯一待求解未知量，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出目标地震固结系数α^sei；

步骤6：将步骤4计算的目标地震临界孔隙度φ_c ^sei、步骤5计算的目标地震固结系数α^sei、地震资料反演获得的横波速度V_s ^sei和饱含流体岩石密度ρ^sei、利用测井资料统计获得的岩石基质剪切模量μ_m输入到本发明的孔隙度计算公式(式6)中，则计算出基于地震资料的孔隙度φ^sei。

图2-图4分别是选取塔里木盆地塔中地区某区块为研究区，基于叠前地震资料利用叠前AVO反演和叠前压噪密度反演得到的纵波速度、横波速度和饱含流体岩石密度反演结果连井剖面，并利用纵波速度、横波速度和密度测井曲线进行反演结果的标定，其地震反演结果与测井曲线相吻合，证明了叠前弹性参数地震反演结果的正确性，从而为地震孔隙度计算提供准确的输入参数。

图5和图6分别是研究区基于测井资料的体积模量与孔隙度曲线交汇图和剪切模量与孔隙度曲线交汇图。假设研究区碳酸盐岩地层岩性比较单一，岩性的变化对孔隙度预测的影响较小，对测井资料统计分析，将孔隙度与体积模量、孔隙度和剪切模量进行交汇，并进行关系拟合，当孔隙度为0时对应的模量参数即为碳酸盐岩基质的体积模量和剪切模量。

图7为步骤3基于测井资料计算的参考测井临界孔隙度与各叠前弹性参数交汇图，通过对比分析发现，其中纵波阻抗与参考测井临界孔隙度线性拟合关系最好，线性拟合度达到了0.9168，因此利用纵波阻抗可以准确预测临界孔隙度。

图8为研究区地震孔隙度计算结果的连井剖面，图中测井曲线为孔隙度测井曲线。发现高孔隙度的储层集中发育在鹰山组碳酸盐岩地层的顶界面附近，并具有很强的非均质性，这是由于研究区在地质历史时期经历了多期次的构造运动，碳酸盐岩地层多次抬升出露地表，其中鹰山组沉积末期受到加里东运动的影响，遭到风化、剥蚀、大气淡水淋滤等作用最为强烈，与上覆地层形成了角度不整合，发育了次生岩溶风化壳储层，预测结果与地质规律相符合，并且地震孔隙度计算结果与测井解释的孔隙度相吻合，证明了该方法的可行性。

图9为研究区目的层地震孔隙度计算结果平面分布图。并选取研究区内投产时间超过1年以上的井，统计其累计产液量(包括产油量和产水量)，以保证统计的累计产液量能够一定程度代表储集空间的大小，并将其与地震预测的孔隙度进行对比(如图10)，发现二者存在一定的相关性，即孔隙度高值的井一般产液量较高，孔隙度低值的井一般产液量较低，说明地震预测的孔隙度与钻井结果相一致，能够较好的反映储层的发育程度，即储集空间的大小，证明本发明提出的地震孔隙度计算方法在碳酸盐岩储层中效果明显，能够利用地震资料在横向上对孔隙度进行相对准确的预测。

本申请实施例还提供了复杂碳酸盐储层地震孔隙度获取装置，包括：

第一获取模块，用于通过叠前反演方法获取目标地域的地震纵波速度、地震横波速度和地震饱含流体岩石密度；

第二获取模块，用于通过统计测井资料，获取目标地域的岩石基质体积模量和岩石基质剪切模量；

第三获取模块，用于根据测井资料、地层压力、温度、地层水的矿化度因素，获取目标地域的流体体积模量；

第一计算模块，用于根据目标地域的地震纵波速度、地震横波速度、地震饱含流体岩石密度，以及测井资料统计的岩石基质体积模量、岩石基质剪切模量、流体体积模量计算目标地域的目标地震临界孔隙度和目标地震固结系数；

第二计算模块，用于根据目标地震临界孔隙度、目标地震固结系数、地震横波速度、地震饱含流体岩石密度和测井资料统计的岩石基质剪切模量计算目标地域的地震孔隙度。

优选的，第一计算模块包括：

第一计算单元，用于使用如下公式计算每个测井的参考测井固结系数， $\sqrt{\frac{k_m\left(\frac{{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{\log }\right)}^2+{\left(V_s^{\log }\right)}^2{\mathrm{\&rho;}}^{\log }{\left(1+{\mathrm{\&alpha;}}^{\log }\right)}^2}{{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{\log }\right)}^2+{\left(V_s^{\log }\right)}^2{\mathrm{\&rho;}}^{\log }{\left({\mathrm{\&alpha;}}^{\log }\right)}^2}\right)}{\left(1-{\mathrm{\&phi;}}^{\log }+{\mathrm{\&phi;}}^{\log }\frac{K_m}{K_f}\right){\mathrm{\&rho;}}^{\log }-\frac{{\left(V_s^{\log }\right)}^2\left({\mathrm{\&rho;}}^{\log }\right)\left(2{\left({\mathrm{\&alpha;}}^{\log }\right)}^2+2{\mathrm{\&alpha;}}^{\log }+1\right)}{{\left(V_s^{\log }\right)}^2{\mathrm{\&rho;}}^{\log }{\left({\mathrm{\&alpha;}}^{\log }\right)}^2+{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{\log }\right)}^2}}+\frac{k_m{\left(V_s^{\log }\right)}^2\left(2{\left({\mathrm{\&alpha;}}^{\log }\right)}^2+2{\mathrm{\&alpha;}}^{\log }+1\right)}{{\left({\mathrm{\&alpha;}}^{\log }\right)}^2{\left(V_s^{\log }\right)}^2{\mathrm{\&rho;}}^{\log }+{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{\log }\right)}^2}+\frac{4}{3}{\left(V_s^{\log }\right)}^2}-V_p^{\log }=0$ ，其中，ρ^log分别为指定测井的纵波速度、横波速度和饱含流体岩石密度的测井曲线数据，φ^log为指定测井的孔隙度测井曲线数据，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^log为需要求解的指定测井的参考测井固结系数。并且，使用如下方法计算指定测井的参考测井固结系数，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出指定测井的参考测井固结系数；

第二计算单元，用于使用如下公式计算每个测井的参考测井临界孔隙度， ${\mathrm{\&phi;}}_c^{\log }=\frac{{\mathrm{\&phi;}}^{log}\&lsqb;{\mathrm{\&alpha;}}^{log}(1+2{\mathrm{\&alpha;}}^{log}){\left(V_s^{log}\right)}^2{\mathrm{\&rho;}}^{log}+{\mathrm{\&mu;}}_m(1+{\mathrm{\&alpha;}}^{log})\&rsqb;}{(1+{\mathrm{\&alpha;}}^{log})({\mathrm{\&mu;}}_m-{\left(V_s^{1og}\right)}^2{\mathrm{\&rho;}}^{log})},$ 其中，为需要求解的指定测井的参考测井临界孔隙度，α^log为指定测井的参考测井固结系数，φ^log为指定测井的孔隙度测井曲线数据，V_s ^log为指定测井的横波速度测井曲线数据，ρ^log为指定测井的饱含流体岩石密度测井曲线数据，μ_m为测井资料统计的岩石基质剪切模量。

交汇模块，用于将每个测井的参考测井临界孔隙度与每个测井的叠前弹性参数进行交汇，以确定多个叠前弹性参数中，目标叠前弹性参数与参考测井临界孔隙度的函数关系，所述目标叠前弹性参数与参考测井临界孔隙度的拟合程度要符合预设的要求；

第三计算单元，用于根据目标叠前弹性参数与参考测井临界孔隙度的函数关系，利用叠前地震反演的目标叠前弹性参数计算目标地震临界孔隙度。

优选的，第一计算模块还包括：

第四计算单元，用于使用如下公式计算目标地域的目标地震固结系数， $\sqrt{\frac{k_m\left(\frac{{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{sei}\right)}^2+{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}{\left(1+{\mathrm{\&alpha;}}^{sei}\right)}^2}{{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{sei}\right)}^2+{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}{\left({\mathrm{\&alpha;}}^{sei}\right)}^2}\right)}{\left(1-{\mathrm{\&phi;}}^{sei}+{\mathrm{\&phi;}}^{sei}\frac{K_m}{K_f}\right){\mathrm{\&rho;}}^{sei}-\frac{{\left(V_s^{sei}\right)}^2\left({\mathrm{\&rho;}}^{sei}\right)\left(2{\left({\mathrm{\&alpha;}}^{sei}\right)}^2+2{\mathrm{\&alpha;}}^{sei}+1\right)}{{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}{\left({\mathrm{\&alpha;}}^{sei}\right)}^2+{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{sei}\right)}^2}}+\frac{k_m{\left(V_s^{sei}\right)}^2\left(2{\left({\mathrm{\&alpha;}}^{sei}\right)}^2+2{\mathrm{\&alpha;}}^{sei}+1\right)}{{\left({\mathrm{\&alpha;}}^{sei}\right)}^2{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}+{\mathrm{\&mu;}}_m{\left(1+{\mathrm{\&alpha;}}^{sei}\right)}^2}+\frac{4}{3}{\left(V_s^{sei}\right)}^2}-V_p^{sei}=0$ 其中， ${\mathrm{\&phi;}}^{sei}=\frac{{\mathrm{\&phi;}}_c^{sei}(1+{\mathrm{\&alpha;}}^{sei})({\mathrm{\&mu;}}_m-{\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei})}{{\mathrm{\&alpha;}}^{sei}(1+2{\mathrm{\&alpha;}}^{sei}){\left(V_s^{sei}\right)}^2{\mathrm{\&rho;}}^{sei}+{\mathrm{\&mu;}}_m(1+{\mathrm{\&alpha;}}^{sei})},{\mathrm{\&phi;}}_c^{sei}$ 为计算的目标地震临界孔隙度，ρ^sei分别为地震纵波速度、地震横波速度和地震饱含流体岩石密度，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^sei为需要求解的目标地震固结系数；

并且，使用如下方法求得目标地震固结系数，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出目标地震固结系数。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的装置、模块和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (10)

1.复杂碳酸盐储层地震孔隙度获取方法，其特征在于，包括：

通过叠前反演方法获取目标地域的地震纵波速度、地震横波速度和地震饱含流体岩石密度；

通过统计测井资料，获取目标地域的岩石基质体积模量和岩石基质剪切模量；

根据测井资料、地层压力、温度、地层水的矿化度因素，获取目标地域的流体体积模量；

根据目标地域的地震纵波速度、地震横波速度、地震饱含流体岩石密度，以及测井资料统计的岩石基质体积模量、岩石基质剪切模量、流体体积模量计算目标地域的目标地震临界孔隙度和目标地震固结系数；

根据目标地震临界孔隙度、目标地震固结系数、地震横波速度、地震饱含流体岩石密度和测井资料统计的岩石基质剪切模量计算目标地域的地震孔隙度。

2.根据权利要求1所述的复杂碳酸盐储层地震孔隙度获取方法，其特征在于，所述计算目标地域的目标地震临界孔隙度包括：

分别计算目标地域中每个测井的参考测井临界孔隙度；

将每个测井的参考测井临界孔隙度与每个测井的叠前弹性参数进行交汇，以确定多个叠前弹性参数中，目标叠前弹性参数与参考测井临界孔隙度的函数关系，所述目标叠前弹性参数与参考测井临界孔隙度的拟合程度要符合预设的要求；

根据目标叠前弹性参数与参考测井临界孔隙度的函数关系，利用叠前地震反演的目标叠前弹性参数计算目标地震临界孔隙度。

3.根据权利要求1所述的复杂碳酸盐储层地震孔隙度获取方法，其特征在于，所述计算目标地域的目标地震固结系数包括：

使用如下公式计算目标地域的目标地震固结系数，其中， 为计算的目标地震临界孔隙度，ρ^sei分别为地震纵波速度、地震横波速度和地震饱含流体岩石密度，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^sei为需要求解的目标地震固结系数；

并且，使用如下方法求得目标地震固结系数，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出目标地震固结系数。

4.根据权利要求1所述的复杂碳酸盐储层地震孔隙度获取方法，其特征在于，根据目标地震临界孔隙度、目标地震固结系数、地震横波速度、地震饱含流体岩石密度和测井资料统计的岩石基质剪切模量计算目标地域的地震孔隙度包括：

使用如下公式计算目标地域的地震孔隙度，其中，φ^sei为地震孔隙度，α^sei为目标地震固结系数，为目标地震临界孔隙度，为地震横波速度，ρ^sei为地震饱含流体岩石密度，μ_m为测井资料统计的岩石基质剪切模量。

5.根据权利要求2所述的复杂碳酸盐储层地震孔隙度获取方法，其特征在于，步骤所述分别计算目标地域中每个测井的参考测井临界孔隙度包括：

使用如下公式计算指定测井的参考测井临界孔隙度，其中，为指定测井的参考测井临界孔隙度，α^log为指定测井的参考测井固结系数，φ^log为指定测井的孔隙度测井曲线数据，为指定测井的横波速度测井曲线数据，ρ^log为指定测井的饱含流体岩石密度测井曲线数据，μ_m为测井资料统计的岩石基质剪切模量。

6.根据权利要求5所述的复杂碳酸盐储层地震孔隙度获取方法，其特征在于，

使用如下公式计算指定测井的参考测井固结系数，，其中，ρ^log分别为指定测井的纵波速度、横波速度和饱含流体岩石密度的测井曲线数据，φ^log为指定测井的孔隙度测井曲线数据，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^log为需要求解的指定测井的参考测井固结系数；

并且，使用如下方法计算指定测井的参考测井固结系数，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出指定测井的参考测井固结系数。

7.根据权利要求2所述的复杂碳酸盐储层地震孔隙度获取方法，其特征在于，所述多个叠前弹性参数包括纵波阻抗、横波阻抗、纵横波速度比和密度。

8.复杂碳酸盐储层地震孔隙度获取装置，其特征在于，包括：

第一获取模块，用于通过叠前反演方法获取目标地域的地震纵波速度、地震横波速度和地震饱含流体岩石密度；

第二获取模块，用于通过统计测井资料，获取目标地域的岩石基质体积模量和岩石基质剪切模量；

第三获取模块，用于根据测井资料、地层压力、温度、地层水的矿化度因素，获取目标地域的流体体积模量；

第一计算模块，用于根据目标地域的地震纵波速度、地震横波速度、地震饱含流体岩石密度，以及测井资料统计的岩石基质体积模量、岩石基质剪切模量、流体体积模量计算目标地域的目标地震临界孔隙度和目标地震固结系数；

第二计算模块，用于根据目标地震临界孔隙度、目标地震固结系数、地震横波速度、地震饱含流体岩石密度和测井资料统计的岩石基质剪切模量计算目标地域的地震孔隙度。

9.根据权利要求8所述的复杂碳酸盐储层地震孔隙度获取装置，其特征在于，所述第一计算模块包括：

第一计算单元，用于使用如下公式计算每个测井的参考测井固结系数，

，其中，ρ^log分别为指定测井的纵波速度、横波速度和饱含流体岩石密度的测井曲线数据，φ^log为指定测井的孔隙度测井曲线数据，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^log为需要求解的指定测井的参考测井固结系数，且，使用如下方法计算指定测井的参考测井固结系数，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出指定测井的参考测井固结系数；

第二计算单元，用于使用如下公式计算每个测井的参考测井临界孔隙度，其中，为需要求解的指定测井的参考测井临界孔隙度，α^log为指定测井的参考测井固结系数，φ^log为指定测井的孔隙度测井曲线数据，为指定测井的横波速度测井曲线数据，ρ^log为指定测井的饱含流体岩石密度测井曲线数据，μ_m为测井资料统计的岩石基质剪切模量；

交汇模块，用于将每个测井的参考测井临界孔隙度与每个测井的叠前弹性参数进行交汇，以确定多个叠前弹性参数中，目标叠前弹性参数与参考测井临界孔隙度的函数关系，所述目标叠前弹性参数与参考测井临界孔隙度的拟合程度要符合预设的要求；

第三计算单元，用于根据目标叠前弹性参数与参考测井临界孔隙度的函数关系，利用叠前地震反演的目标叠前弹性参数计算目标地震临界孔隙度。

10.根据权利要求9所述的复杂碳酸盐储层地震孔隙度获取装置，其特征在于，所述第一计算模块还包括：

第四计算单元，用于使用如下公式计算目标地域的目标地震固结系数，

其中， 为计算的目标地震临界孔隙度，ρ^sei分别为地震纵波速度、地震横波速度和地震饱含流体岩石密度，K_m、K_f、μ_m分别为测井资料统计获得的岩石基质体积模量、孔隙内流体体积模量和岩石基质的剪切模量，α^sei为需要求解的目标地震固结系数；

并且，使用如下方法求得目标地震固结系数，利用混沌量子粒子群非线性优化算法进行全局搜索，当计算误差小于设定的阈值，则输出目标地震固结系数。