0 引言

含水饱和度是储层评价极为重要的参数之一。在非均质储层中，受储层微观孔隙结构、平均颗粒尺寸、胶结物类型及其含量、裂缝及泥质含量等因素的影响，含水饱和度的准确评价较常规储层难度更大。现有的储层评价中确定饱和度的方法主要包括利用密闭取心直接测定^[1]、实验室岩心测定方法^[2]及地球物理测井方法^[3]。相比而言，利用密闭取心直接测定的岩心含水饱和度是现场第一手资料，但成本高，推广难度大。实验室岩心测定法是获取含水饱和度较可靠的方法，但岩心分析数据的离散特征使储层定量表征存在明显的不连续性，因而通常只作为理论分析的验证依据。目前使用最为广泛的评价储层饱和度的方法是基于测井资料的含水饱和度评价，主要基于电法测井数据^[3]、核磁共振(NMR)测井数据^[4]及声波测井数据^[5]等。基于室内岩心电性实验和理论分析构建饱和度评价模型，进而结合测井曲线和多个影响因素进行储层含水饱和度评价，是目前储层含水饱和度评价的重要手段^[6]。

目前，基于实验室数据和测井数据的NMR技术正迅速成为评价储层质量的岩石物理手段。对常规储层，常用Archie经验公式评价其含水饱和度^[7]。Archie公式中的参数不但与地层岩性、孔隙结构有关，还与地层水矿化度、围压和温度相关，因此，对于低孔、低渗储层的含水饱和度评价，传统Archie公式无能为力^[8-9]。Waxman等^[10]研究发现黏土矿物具有附加导电性，提出基于阳离子交换容量计算的W-S模型；Clavier等^[11]将黏土矿物表面导电层分为受黏土矿物成分影响的、靠近黏土矿物的表层电子层及不受黏土影响的、远离黏土矿物的自由电子层，提出了改进的W-S模型。随着对油藏成因认识的不断深入，以复杂孔隙结构导电机理为基础开展饱和度模型研究逐渐成为油藏储层评价的热点。Shang等^[12]提出了等效岩石组分饱和度(EREM)模型，将组成孔隙空间的孔隙和喉道分成并联和串联孔隙，定义了孔隙结构系数和饱和度刻度因子，并给出了明确的参数计算方法，有效提高了孔隙结构复杂储层饱和度计算精度。邵才瑞等^[13]通过统计不同毛管压力曲线形态的物性特征，按孔隙度范围分类建立“J函数”与饱和度之间的统计关系。Hu等^[14]将致密砂岩储层的导电性等效为孔隙中截面积不变的直孔隙与截面积变化的梯形孔隙的并联模型，建立了考虑迂曲度、直孔隙比例和梯形因子等参数的梯形孔隙饱和度模型。

对于具有一定泥质含量及复杂储层孔隙结构的储层，孔隙结构及孔隙流体分布形态等非饱和度因素都会影响岩石的电阻率^[15-16]，而且这种影响在储层的不同部位是逐点变化的。现在通用的Archie饱和度模型^[7]、W-S模型^[10]都采用不同参数设法消除流体饱和度之外的因素对电阻率的影响。但是，实际的测井数据处理过程中，由于无法针对储层中不同部位所受影响的差异性实时地改变饱和度参数，导致复杂储层饱和度的计算精度明显受到限制^[17-19]。在实际的测井数据处理过程中，考虑到储层非均质性的影响^[20]，一般采用分层的方式确定饱和度模型的参数^[21]，这在一定的程度上可降低储层非均质性对饱和度计算精度的影响。但是，分层确定饱和度参数的方法工作量很大，也难以完全消除严重非均质性储层对于饱和度计算精度的影响^[22-23]。

为此，本文基于核磁岩心实验数据、含油饱和后的水驱替岩石导电实验数据、实际核磁测井数据及球管模型优化反演数据，建立储层孔隙结构与岩石导电模型参数之间的适配关系。测井数据处理过程中，依据这种关系实现逐点自动拾取饱和度模型参数，达到提高复杂储层饱和度计算精度的效果。

1 研究背景

冀东油田研究工区碎屑岩储层发育了复杂的孔隙储集及渗流系统，中深层油藏中低饱和度油气储层发育，储层埋藏较深、成岩后生作用强，孔隙结构复杂，几乎无法采用解析方法描述(图 1)。储层不动水饱和度高导致油、水层测井响应差异小，因而难以准确计算储层的含水饱和度。

本文基于岩心核磁实验数据、含油饱和后的水驱替岩石导电实验数据、核磁测井数据及常规测井数据，采用优化反演算法及相似对比方法，基于岩石孔隙结构分类计算得到相应的饱和度模型参数。在测井数据处理中，基于C_d(即孔隙中管形半径与球形半径之比)路径的相似性自动拾取与储层孔隙结构对应的饱和度模型参数，据此消除孔隙结构和孔隙流体分布对岩石导电性的影响，实现提高复杂储层饱和度计算精度的目的。

2 储层孔隙的球管模型刻画方法

真实储层岩石的孔隙结构非常复杂，几乎无法用解析方法描述。但是，在合理的近似条件下，可以采用某些特定模型实现岩石孔隙结构的近似描述。本文采用球管模型对储层岩石进行分类，建立球管模型的核心思想是：将岩石中的孔隙近似为一系列不同类型的球形孔(孔隙部分)和管形孔(孔喉部分)的组合。球形孔与管形孔的不同匹配方式代表不同类型的孔隙结构(图 2)^[24]。

图 2展示了不同半径的球形孔与管形孔的配置关系。根据球管模型储层评价方法可得到岩石孔隙中分组球管模型的管形孔半径^[24]。利用岩心核磁数据开展基于球管模型的优化反演，可得到T₂总谱^[25]、球形孔谱、管形孔谱和C_d路径(图 3)。其中，C_d路径反映每一个反演布点弛豫组分中管半径与球半径的配置关系

$ C_{\mathrm{d} i}=\frac{R_{\mathrm{c} i}}{R_{\mathrm{s} i}} $

(1)

式中：i=1，2，…，表示孔隙分组号；R_ci为第i组孔隙中管形孔隙半径；R_si为第i组孔隙中球形孔隙半径；C_di表示第i组孔隙中管形孔半径与球形孔半径之比，取值范围为[0, 1]。

因此，确定孔隙结构的关键问题转化为确定C_d的值。理论上，每个组分内的球管模型都有无数种，为了方便计算，本文提出利用扫描方法计算孔隙结构参数(图 4)。图中横轴是横向弛豫时间的分组编号i，这个分组编号对应的弛豫时间是T_2i，也就是反演核磁回波的T₂布点值。图中每一条折线都代表岩石中某一种特定的球形孔与管形孔的配置关系，即一种孔隙结构特征。沿每一条折线从左到右，都可以计算横向弛豫时间分组对应的T₂布点值T_2i、R_ci、R_si，利用这一组的T_2i布点值反演核磁回波。当某一条折线对应的布点值(T_2i)使核磁回波的拟合误差达到最小时，这个T₂布点值对应的C_d值C_di、球形孔半径R_si、管形孔半径R_ci即岩石孔隙结构的最优化近似描述。

本文对采自南堡11-L8X204井和南堡23-X2282井的13块岩心核磁数据进行优化反演，得到这13块岩样的C_d路径(图 5)。根据各条C_d路径之间的相似度，孔隙结构可以分为4组，分组算法如下。

首先，选定一个基本的C_d路径，根据

$ r=\frac{\sum\limits_{i=1}^N\left(C_{\mathrm{d} 1 i}-\overline{C_{\mathrm{d} 1}}\right) \sum\limits_{i=1}^N\left(C_{\mathrm{d} 2 i}-\overline{C_{\mathrm{d} 2}}\right)}{\sqrt{\sum\limits_{i=1}^N\left(C_{\mathrm{d} 1 i}-\overline{C_{\mathrm{d} 1}}\right)^2 \sum\limits_{i=1}^N\left(C_{\mathrm{d} 2 i}-\overline{C_{\mathrm{d} 2}}\right)^2}} $

(2)

计算该C_d路径与其他C_d路径间的相关系数。式中：N表示总路径数；C_d1i、C_d2i分别为第1条C_d路径和第2条C_d路径在第i个布点的C_d值；C_d1、C_d2——分别是第1条和第2条C_d路径的平均值。

然后，根据r最大化原则，把C_d路径相关系数最大的岩心分为一组(图 6)。

3 改进的EREM模型及含水饱和度计算方法

研究表明，对于复杂孔隙结构，岩石电阻率受岩石孔隙度的影响，且孔隙度与地层因子之间也不是简单的线性关系。Shang等^{[12, 24]}研究了复杂岩石孔隙结构对地层因子的影响，验证了孔隙度与地层因子之间的非线性关系EREM模型。本文基于文献[26-27]中的改进EREM模型，使孔隙结构参数参与岩石含水饱和度的计算过程，利用双孔隙(球形孔和管形孔)饱和度模型提高含水饱和度S_w的计算精度。相关计算公式如下

$ 1 \mathrm{~g} R_{\mathrm{F}}=c 1 \mathrm{~g} \phi+d(1 \mathrm{~g} \phi)^2 $

(3)

$ \lg R_{\mathrm{I}}=e 1 \mathrm{~g} S_{\mathrm{w}}+f\left(1 \mathrm{~g} S_{\mathrm{w}}\right)^2 $

(4)

$ F=\frac{R_0}{R_{\mathrm{w}}}=\frac{(1-\phi)^2}{C_{\mathrm{F}} \phi}+\frac{1}{\phi} $

(5)

$ I=\frac{R_{\mathrm{t}}}{R_0}=\frac{\left(1-S_{\mathrm{w}}\right)^2}{C_{\mathrm{I}} S_{\mathrm{w}}}+\frac{1}{S_{\mathrm{w}}} $

(6)

$ C_{\mathrm{F}}=\frac{C_0 R_{\mathrm{F}}}{1+\left(1-R_{\mathrm{F}}\right) C_0} $

(7)

$ C_{\mathrm{I}}=\frac{C_{\mathrm{F}} R_{\mathrm{I}}}{1+\left(1-R_{\mathrm{I}}\right) C_{\mathrm{F}}} $

(8)

$ C_0=\frac{V_{\mathrm{c}}}{V_{\mathrm{s}}} $

(9)

式中：R_F为孔隙度刻度因子，反映孔隙结构和分布对岩石导电性的影响；R_I表示饱和度刻度因子，反映含烃数量和分布形态对岩石导电性的影响；c、d是与孔隙结构相关的系数；e、f是与含烃饱和度及分布有关的系数；F、I分别为地层因素和电阻增大系数；R₀、R_t分别为饱含水和含油岩石样品的电阻率；R_w表示实验配置盐水的电阻率；ϕ表示地层的孔隙度；C_F表示初始孔隙结构优度；C_I是饱和度结构优度，反映含烃分布形态和孔隙结构形态对岩石导电性的影响；C₀为初始孔隙结构优度，反映孔隙形态结构对岩石导电性的影响，其值可根据球管模型的分布曲线得到，本文取值为管形孔、球形孔累计分布曲线的中位值的比值；V_c表示管形孔隙体积；V_s表示球形孔隙体积。

Shang等^[24]提出的EREM模型主要基于式(3)~式(9)讨论储层孔隙结构对岩石导电性的影响。本文主要从三个方面改进EREM模型，并将改进成果应用于复杂储层含水饱和度的计算。首先，提出实用的模型参数(c、d、e、f)确定方法；其次，将改进的EREM模型用于含水饱和度的计算；最后，根据孔隙结构参数的计算结果逐点自动拾取模型参数(c, d, e, f)。从式(6)可以看出，EREM模型的核心思想是放弃了胶结指数和饱和度指数为定值的概念^[7]，转而使用函数关系反映孔隙结构和饱和度分布对岩石导电性的影响，使参数具有了明确的物理含义，其核心参数就是c、d、e、f这4个系数。

下面以系数c、d为例，说明其值的确定方法。确定系数e、f的方法与其类似，不再赘述。

(1) 根据式(5)和式(7)计算岩心的C_F；

(2) 根据球管模型的计算结果(式(1))计算V_fi和V_si，进而得到C₀=V_ci/V_si，这里取V_ci/V_si比值的中位值；

(3) 根据式(3)计算R_F，代入下面的矩阵

$ \boldsymbol{Y}=\left(\begin{array}{c} 1 \mathrm{~g} R_{\mathrm{F}_1} \\ 1 \mathrm{~g} R_{\mathrm{F} 2} \\ \vdots \\ 1 \mathrm{~g} R_{\mathrm{F} n} \end{array}\right) $

(10)

式中n表示实验数据点数。

(4) 根据岩心的孔隙度ϕ计算矩阵

$ \boldsymbol{A}=\left(\begin{array}{cc} 1 \mathrm{~g} \phi_1 & 1 g^2 \phi_1 \\ 1 g \phi_2 & 1 g^2 \phi_2 \\ & \vdots \\ 1 g \phi_n & 1 g^2 \phi_n \end{array}\right) $

(11)

(5) 基于式(10)、式(11)得到参数c、d的矩阵

$ \left(\begin{array}{ll} c & d \end{array}\right)^{\mathrm{T}}=\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}+\lambda \boldsymbol{I}\right)^{-1} \boldsymbol{A}^{\mathrm{T}} \boldsymbol{Y} $

(12)

式中：I为2×2的单位对角矩阵；λ为阻尼系数，主要作用是增强计算结果的稳定性和合理性。如果计算得到的F-ϕ关系曲线随孔隙度ϕ的降低保持单调递增的趋势，则上式中λ取值为零；如果F-ϕ关系曲线出现中间凸起的现象(即非单调递增或递减)，则λ需从0开始逐渐增大进行试算，直到F-ϕ关系曲线呈单调递增趋势。

关于球管模型数值模拟及导电性关系的研究已经证明，岩石的导电性与孔隙结构有密切的关系，也就是说，只有在孔隙结构背景一致的前提下，讨论岩石的导电性与含油性之间的关系，对于提高含油饱和度的计算精度才有意义。因此，确定改进EREM模型参数(c、d、e、f)时，应按照孔隙结构背景基本一致的原则，把同一井区的储层数据分为不同的组，分别确定改进EREM模型的参数。

根据研究区已经取得的储层分类研究成果，将储层分为4个类型。对应地，根据C_d路径之间的相似性，把饱和度模型的参数也分为4个类型，即针对研究区的储层类型(图 6)，分别确定相应的饱和度模型参数。

对于球管模型，表征管形孔与球形孔匹配关系的关键参数是C_d，即C_d相同或者近似的岩石孔隙的孔隙结构也非常近似。因此，可根据C_d路径判断岩石孔隙结构的相似程度，并根据C_d路径特征对岩石孔隙类型进行分类。这个路径分类代表了各组内部管形孔的宽度与球形孔半径之间的比值关系，所以C_d路径及路径的分布方式可刻画岩石内部孔隙结构的特征。

测井中含水饱和度的计算是一个迭代计算过程，这个过程是收敛的。改进的EREM模型含水饱和度迭代计算过程^[27]如下。

(1) 采用典型Archie公式^[7]计算含水饱和度初值

$ S_{\mathrm{w}}^{(0)}=\sqrt{\frac{R_{\mathrm{w}}}{\phi^2 R_{\mathrm{t}}}} $

(13)

(2) 根据改进EREM模型(式(6)、式(8))，计算参数I和C_I；

(3) 求解改进EREM模型(式(6))中关于S_w的一元二次方程，得到S_w的计算公式

$ S_{\mathrm{w}}=\frac{\left(2+C_{\mathrm{I}} I\right)-\sqrt{\left(2+C_{\mathrm{I}} I\right)^2-4\left(C_{\mathrm{I}}+1\right)}}{2} $

(14)

(4) 计算更新值S_w^(k)与初值S_w⁽⁰⁾的误差ε，这里k表示迭代次数。如果ε大于给定的阈值，则赋值S_w⁽⁰⁾=S_w^(k)，并返回步骤(2)，重新计算S_w，直到ε满足要求。

总结以上储层饱和度计算过程，其具体步骤如下：首先，基于球管模型优化反演核磁岩石实验数据，计算孔隙结构参数，即C_d路径；其次，对多个岩石样本的C_d路径进行相关性分析，把具有相似C_d路径的样本分为一组，并对每组数据采用阻尼最小二乘法，确定改进EREM饱和度模型的参数；然后，通过把岩石C_d路径和对应的改进EREM饱和度模型参数应用于测井数值处理程序，对核磁测井数据进行优化处理，得到井筒剖面储层的C_d路径参数，逐条计算核磁测井反演得到的C_d路径与岩心分组C_d路径之间的相关系数，并拾取相关系数最大的岩石C_d路径对应的饱和度参数，作为该储层点的饱和度值；最后，基于已拾取的饱和度模型参数和常规测井资料(电阻率数据和孔隙度数据)，采用改进EREM模型计算该储层点的饱和度。

4 实际应用效果评价

将本文方法应用于冀东油田南堡油田井区实际测井资料的处理，结果见图 7。图 7直观地比较了测井C_d路径与基于核磁数据反演的C_d路径的相似性，验证了本文含水饱和度模型参数选取方法是正确的。本文方法将孔隙结构参数融入饱和度计算过程，相较于传统的储层饱和度计算过程采用固定参数的方式有明显的改进效果。

图 8是工区南堡23-2694井基于本文自动拾取饱和度模型参数对测井数据进行处理的结果，其中“解释结论”中关于油水的信息来自实际生产数据。由图可见，该井段一共解释了15个储层段，包括5个油层、2个油水同层、4个水层及4个干层。其中，含油层段中基于本文新模型得到的S_w为15%，而原模型的S_w为35%~40%。分析其原因，是新模型剔除了孔隙结构对于岩石电阻率测量值的影响，所以计算的含水饱和度更接近真实值。在含水层段，基于新模型和原模型计算的含水饱和度均为100%。分析C_d路径的形态，可以看到C_d路径在储层内部和非储层中都有比较剧烈的变化(图 8)，说明储层和非储层都有比较明显的非均质性。该井段电阻率曲线变化也比较剧烈，即电阻率的剧烈变化也表明该井段储层物性的非均质性特点。因此，核磁优化反演结果(图 8中C_d路径)与电测井响应均反映出该井段地层的非均质特征，而且前者与实际生产情况吻合较好，且明显提高了储层饱和度的计算精度。

5 结论

本文基于核磁测井回波数据优化反演结果，逐点计算储层孔隙结构参数，并根据孔隙结构的相似性，把储层孔隙结构分成4个不同的类型。不同类型之间的储层孔隙结构具有比较显著的差异；同类型的储层孔隙结构的形态具有一定的相似性。在此基础上，根据不同类型的孔隙结构，确定相应的饱和度模型参数。该项研究有利于针对强非均质性储层消除孔隙结构差异对于计算储层含水饱和度的影响，提高含水饱和度的计算精度。基于以上分析，可得出如下结论：

(1) 基于核磁岩心实验数据、含油饱和后的水驱替岩石导电实验数据、真实核磁测井数据及球管模型优化反演成果数据，建立储层孔隙结构与岩石导电模型参数之间的适配关系，即C_d路径，该参数反映了孔隙中球形孔与管形孔的适配关系；

(2) 依据C_d路径所代表的适配关系，在测井数据处理过程中，实现逐点自动拾取含水饱和度模型参数，作为该储层点的饱和度参数。

本文方法应用于冀东油田南堡油田井区实际测井资料，显著提高了储层饱和度的计算精度。