0 引言

碳酸盐岩油气藏是我国重要的战略接替资源。不同于国外中、新生代碳酸盐岩形成的大油气田，我国油气储层多处于古老的海相地质条件中，这意味着它们的沉积过程更有可能受到差异化生物、化学和机械作用的控制。碳酸盐岩在经历了多期次的地质演化和改造后，储层孔隙结构的复杂性提高，地球物理观测资料的多解性增强。这就要求数据处理解释工作中应尽力提升宏观观测资料的利用率，挖掘其中隐含的微观信息量，为储层地质解释提供有价值的参考数据和模型。

针对碳酸盐岩储层的岩石物理研究起步较晚。实验分析方面，Birch^[1]在20世纪60年代就开始针对碳酸盐岩进行变围压岩石物理实验; Nur^[2]与Mavko等^[3]根据不同地区样品的实验分析数据，发现碳酸盐岩的弹性性质变化主要与其内部所含的微观力学缺陷有关; Scotellar等^[4]认为压实、溶蚀、白云岩化和钙质胶结等过程通过对微观结构的影响极大程度上改变了碳酸盐岩的弹性性质; Wang^[5]进一步指出影响碳酸盐岩波速的主要因素为孔隙多少和孔隙类型。数值建模方面，碳酸盐岩的岩石物理模型通常将孔隙抽象成有不同包含物的椭球体来求得样品的等效弹性模量^[6-7]。

本研究旨在提供一套碳酸盐岩等效介质模型的改进方案，并形成与之配套的岩石微观参数提取方法和地层类型识别思路，为地质解释和储层预测提供参考数据和模型。

1 理论与方法 1.1 碳酸盐岩的岩石物理模型

有别于常规砂岩、泥岩储层比较单一的微观构造，碳酸盐岩储层的形成通常伴随着多期次的构造活动和多元的溶蚀改造，岩石内部往往同时存在3类以上的孔隙，它们的类型和含量是地质解释人员关心的问题^[8]。怎样从现有的地球物理资料中量化这些差异？这就需要借助岩石物理工具先从理论上认清微观特征与宏观观测数据的对应关系。碳酸盐岩中不同成因的孔隙往往具有不同的几何形态，因此岩石物理学中通常将复杂微观孔隙空间简化为3个部分：①以溶孔、窗格孔和铸模孔为代表的球状孔隙(硬孔隙)，孔隙度记为φ_SP；②以构造裂缝和溶蚀裂缝为代表的扁平孔隙(软孔隙), 孔隙度记为φ_M；③以晶间孔、粒间孔为代表的次球状孔隙, 孔隙度记为φ_IP。岩石的总孔隙度φ为三者之和：

等效介质建模是将物理现象转化为数学表达的过程。对于碳酸盐岩，商业软件通常通过微分等效介质模型(DEM)的迭代运算实现这一过程^[9]。该类方法预测的弹性参数与实验数据的可比性较强，但在实际生产应用中也有突出的局限性：①模型计算的是一种数学过程而不是实际物理过程；②涉及到大量不解耦参数和偏微分方程计算，需要通过迭代方法进行数值求解，难以兼顾计算效率；③数值计算具有路径依赖性，也就是说交换矿物背景和孔隙组分次序会改变预测结果，这与实际观测结果相悖；④建模过程涉及太多的岩石微观参数，如孔隙纵横比和体积分数，这些参数在井中难以获得。

为了克服这些不足，本研究在计算干岩石骨架(只包含孔隙和矿物)模量时借鉴了经典的临界孔隙度模型(critical porosity model, CPM)^[10]，即

式中：K_m和K_dry分别为矿物混合物和干岩石的基准体积模量；φ_c为临界孔隙度；μ_m和μ_dry分别为矿物混合物和干岩石的剪切模量。当造岩矿物单一时(深层岩石常认为是白云石与方解石的混合物)，K_m可以由单矿物组分含量表示。

临界孔隙度反映的是矿物骨架对外力的承受能力，这种能力的强弱与岩石的矿物成分和骨架微观结构相关(具体函数表达式参见Keys等^[11])。图 1是不同储集岩性的理论临界孔隙度随孔隙纵横比的变化规律，可以看出该物理量除了与组成岩石的矿物基质有关，更依赖于岩石中的微观孔隙成分。“临界”的含义为门槛值，它划分了饱和岩石两种截然不同的支撑机制：①在承载域(φ < φ_c)，介质整体呈固相，由矿物骨架承受压力起支撑作用；②在悬浮域(φ> φ_c)，介质整体呈流相，骨架颗粒在流体中呈悬浮态，流体承担起支撑作用。

大量实验数据表明：φ_c在砂岩这种孔隙结构简单的储层内变化很小(通常在0.16~0.22)，在碳酸盐岩地区则变化范围很大(0.02~0.50)。考虑到沉积环境的稳定性(岩性差异不大)，这种数值的变化可以归因为孔隙成分差异。在宏观尺度上看，φ_c反映了矿物混合物的坚硬程度(体积模量)和所包含的孔隙几何形态。坚硬矿物与硬孔隙组合而成的骨架临界孔隙度较大，反之较小。也就是说在承载域，干岩石的基准体积模量与临界孔隙度成正比，与孔隙度成反比：

真实地层条件下的岩石是处于流体饱和状态的，且通常满足承载域的假设。此时，碳酸盐岩的弹性参数可以使用Murphy^[12]改写的Gassmann方程形式进行描述：

式中：K_sat与μ_sat分别为饱和岩石的体积模量和剪切模量；K_p为孔隙流体引起的附加体积模量。

式中，K_fl为流体体积模量。当岩石由n种矿物或流体按照某种体积比例构成时，混合材料的等效弹性模量可以通过体积平均理论给出：

式中：下标V和R分别代表了矿物的Voigt和Reuss界限；K_l和v_l分别为第l种流体的体积模量和体积分数。在确定饱和岩石弹性模量(K_sat和μ_sat)、测定密度(通过密度测井获得，记为ρ)之后，可以进一步计算自然状态下饱和岩石的纵、横波速度(v_P、v_S)：

综上，新模型通过扩展经典的临界孔隙度模型建立了适用于碳酸盐岩的等效介质理论。优势在于：①规避迭代过程，大幅降低计算成本；②不存在孔隙添加的路径依赖性，提高鲁棒性；③临界孔隙度是具有物理意义的岩石参数，可以多级质控。

1.2 基于统计学的岩石特征参数提取方法

研究碳酸盐岩储层微观参数，传统方法(如声波地层因素公式或三孔隙度估算方法)的致命缺陷在于未考虑孔隙形态的因素。新策略为：结合岩石物理模型和概率统计学理论，尝试同时获取孔隙度(原生孔隙和次生溶蚀程度)和临界孔隙度(矿物组分和孔隙发育形态)两个物理量作为度量微观特征的参考指标。

概率统计学方法在模型参数估计的过程中利用了数据本身的分布特征，因此对随机噪声和野值有一定程度的免疫效果。虽然在物理意义上，φ和φ_c是彼此独立的物理量，但是从统计分析的结果来看，二者在分布上通常呈现正相关规律。也就是说，对于自然沉积和构造演化形成的碳酸盐岩来说，储层孔隙空间的大小取决于它的地质成因。例如，大规模的孔洞在概率上更可能是溶蚀孔洞或窗格孔，而不是粒间孔或裂缝等其他成因，此时岩石同时表现出较大的φ和φ_c观测值。

统计学反演的理论基础是贝叶斯公式：

式中：M为待反演参数向量[φ, φ_c]；d为宏观观测数据，建议使用测井中常见的弹性参数组合(如纵、横波速度和密度)；P(M)与P(d)为模型与数据的概率分布；P(M|d)与P(d|M)为条件概率分布。本次研究中后验概率的表达式为

式中：P(φ, φ_c)为先验概率；P(d|φ, φ_c)为似然函数。

从1.1节中我们得知，d是孔隙度、临界孔隙度、矿物基质和孔隙流体的函数：

式中，S_w为流体饱和度。不妨假设同一地区实验室岩心测试数据与测井数据具有统计可比性，那么先验概率P(φ, φ_c)可以从同一地层的拟合关系中推出。对于连续沉积的地层，M的联合概率满足高斯分布特征：

式中：(φ^(k), φ_c^(k))为当前地层k中待求参数的联合概率分布；N(·)为二元高斯分布，可由两对变量(μ_φ^(k)、μ^(k)_φc和σ_φ^(k)、σ^(k)_φc)的均值μ和标准差σ表示。

然后，似然函数被定义为预测误差相反数的指数形式。也就是说，模型预测误差较小的决策变量输入对应较高的似然值，反之对应低值：

式中：s为归一化系数；ΔE为预测误差。进一步定义ΔE为模型速度预测算子g(φ, φ_c)与观测值d之差的绝对值：

总结以上，基于概率统计学的参数估计方法分为以下3个实施步骤：

1) 先验概率计算：进行岩心样品归位和测井层位标定工作。在实验室中获得样品的孔隙度、声波速度和微观结构信息。如图 2a，通过拟合φ和φ_c的联合分布得到决策变量的先验概率(岩心尺度的临界孔隙度计算方法将在下一节中给出)。

2) 似然函数计算：以临界孔隙度模型为正演算子，输入原始(或解释)测井数据，预测不同决策变量组合输入下的弹性参数响应，比较模型输出与实测曲线。如图 2b，根据误差函数获得似然函数。

3) 后验概率计算：如图 2c，遵循最大后验概率准则(maximum a posteriori，MAP)，输出决策变量的估计值M。

1.3 实验数据的验证

首先，给出在岩心尺度的临界孔隙度模型参数化方法。联合式(4)和式(7)，约去K_sat，得到

实验室的测试环境相对可靠(数据精度高且具有可重复性)，不妨将孔隙度测量值当做已知量对待。则式(14)即可进一步整理为φ_c的弹性参数表示形式：

现通过井中碳酸盐岩岩心样品(图 3)验证该参数提取方法。样品来自不同的井位、深度和地层，具有不同矿物组分和微观结构特征。如下准备工作已经完成：①微观结构扫描；②地层压力下的孔隙度测试；③饱水状态下的三轴超声波速度测试。④根据理论公式(15)量化的临界孔隙度值。以上数据被认为是岩石的真实属性(即M和d)。

参数提取结果用交会分析的方式进行评估，旨在考察反演结果(估计值)与实验观察(真实值)的一致性。如图 4所示，决策变量的所有样点均匀分布在对角线周围。回归系数(R²)达到0.891 3和0.969 1，平均相对误差(δ)小于10%，验证了新方法的可靠性。

2 应用效果 2.1 工区地质概况

YS1井位于我国西部地区的某盆地构造带，属局限台地潮坪、蒸发台地体系。本次研究的目标层为5 000~6 500 m埋深的超深层三叠系碳酸盐岩气藏。储层纵向变化快，在浅水环境下准同生白云化作用明显，在深水环境下则发育含泥质的灰岩。成

岩期间受到海平面扰动影响，部分区段夹杂了很薄的藻、砂屑灰质隔层^[13]。我们根据岩性纵向上的这种变化特征，将研究区段划分成了上储层段(5 702~5 730 m)、隔层(5 730~5 750 m)和下储层段(5 750~5 800 m)三部分。

此外，从储集空间上来看，YS1井有明显纵向产能差异，孔隙类型和成分比较复杂。早期铸模孔、粒内溶孔在埋藏成岩阶段被溶蚀扩大，流体溶蚀了颗粒间的胶结物及颗粒，形成了粒间溶孔、粒内溶孔、晶间溶孔及微裂隙。同时，储层在剪切和拉张的构造作用与溶蚀作用的共同影响下，局部常见构造裂缝及沿构造溶蚀的扩大缝, 形成油气运移的通道。

2.2 储层特征参数提取

从YS1井和邻近的2口井位收集碳酸盐岩岩心74块，经过基本的加工和测试分析步骤之后(薄片鉴定、CT扫描、X光衍射实验、覆压孔渗实验和三轴速度测试等手段)，掌握到了目的层段弹性参数和物性参数的基本分布特征。如图 5a、b所示，速度的变化趋势与微观参数呈现出明显的正相关分布特征。图 5c、d为临界孔隙度模型的预测结果，模型与数据之间达到较好的匹配效果。将这种分布规律进行定量统计分析，可以得出所有目的层的联合概率分布，进而得到各自的先验概率。

图 6举例说明了从测井数据中提取参数的标准流程。从图 6可以看出，受到沉积条件影响，碳酸盐岩储层即使在较小的深度范围内(5 702 ~5 800 m)孔隙类型也可能变化很大，导致先验概率(第1列)差异明显。本文提出的统计学方法通过将这种微观差异与模型预测进行匹配(第2列)，达到量化碳酸盐岩孔隙结构(第3列)的目的。

图 7是岩石物理研究中常用的孔隙度、速度交会分析图，通过这样交会分析可以有效地检测数据质量和分析参数聚类特征。数据颜色为临界孔隙度提取结果，参考线为按典型岩性和孔隙类型输入时的模型预测结果。图 7a、c为工业中常用的“三孔度测井”分析法，这种方法没有考虑岩石物理规律，不同类型的孔隙在交会图上没有表现出很好的聚类特征。图 7b、d为新方法得到的预测结果。在速度信息相同的情况下(仅对孔隙度做重新估算)，可以观察到数据规律性的明显增强(数据成簇分布且表现出差异化的孔隙依赖特征)。不同类型的孔隙在速度变化趋势上具有明显的区分度，刚性孔隙(临界孔隙度大的岩石)有较大的孔隙度和较缓的速度变化趋势，柔性孔隙(临界孔隙度小的岩石)表现为较小的孔隙度和较快的变化趋势，这与理论和实验岩石物理规律是一致的。即不同类型储层的弹性响应特征差异更易于识别，有效降低了资料的多解性。

2.3 岩石物理分析及沉积解释

如上文所述，碳酸盐岩气藏的储集空间类型多样，既有受组构控制的粒间溶孔、粒内溶孔、格架孔和晶间溶孔等，又有不受组构控制的溶洞、溶缝和构造缝^[14]。碳酸盐岩储层中临界孔隙度对微观结构的变化十分敏感，其纵向分布特征能为储层地质解释服务。以YS1井为例，假设不同孔隙成分的岩石对应参数的变化范围不同，通过实验室数据和岩心观察结果可以找到这个范围并赋值序号标签，得到量化识别标准(如表 1所示，典型的岩性-孔隙类型组合被分别定义标签管理)。从岩石物理的角度上：①脆性较强的灰岩中普遍发育微裂缝且难以保存球状孔隙，数值上为较低的φ_c响应。这样的微观构造不是油气的主要储集空间，但是可以增强储层整体的沟通能力，是值得关注的特征规律。②刚性较强的白云岩对应着出色的抗压能力，有能力保存岩石中的一些大孔隙。坚硬的骨架在岩石物理响应中表现为较高的φ_c响应，是油气的主要储集空间。

此外，利用地球物理参数纵向上的变化可以推测地区沉积历史和地质运动。图 8为研究区原始测井曲线与岩石物理分析结果，根据数据特征并结合地质背景资料，合理地得到以下推断：

1) 上储层段：潮间带的灰坪、云坪沉积环境，形成了大量层状分布的准同生泥、粉晶白云岩，其规模主要受局部水体咸化程度控制。表生期的淡水溶蚀作用使得白云岩储层中主要发育晶间(溶)孔及少量粒间溶孔。后期构造运动在脆性较强的灰岩中形成微裂缝。

2) 隔层：潮下带的藻灰坪沉积环境，砂屑灰岩夹杂白云质灰岩形成致密隔层，孔隙度很低，且以微裂缝类型为主。其中，在构造和有机酸作用下形成的早期裂缝是油气运移的有利通道，而高温热液造成的晚期裂缝则可以对溶蚀过程起到疏导作用。

3) 下储层段：潮间带的云坪、藻云坪沉积环境，形成大规模层状分布的白云岩，有利于孔隙发育和保存。岩石中顺层理方向发育了大量的鸟眼、窗格孔以及藻间溶孔，是极佳的油气储存空间。

3 结论

1) 改进的岩石物理建模方法可以有效地描述复杂微观结构和宏观响应之间的对应关系，降低地球物理资料的多解性，为弹性参数响应差异提供机理解释。

2) 碳酸盐岩临界孔隙度估算方法充分兼顾了岩石物理普遍规律与地区差异性特征，具有良好的计算精度。新方法可以提高地球物理资料对储层微观特征的描述能力，是分析储层沉积和改造历史的有利工具。

3) 以西部地区的YS1井为实例展示了新的孔隙微观特征分析流程，即岩石物理建模—特征参数提取—储层地质解释。应用效果充分证明了新方法的实施有助于碳酸盐岩储层微观特征分析，具有进一步完善和推广的价值。