用积分法求取岩石孔隙结构特征值

李 闽，陶正武，刘全稳，吴泽民，李 涛，肖文联，亢 鞠

（1.西南石油大学 油气藏地质及开发工程国家重点实验室，四川 成都 610500；2.广东石油化工学院 石油工程学院，广东 茂名 525000；3.中石油长庆油田 超低渗第四项目部，甘肃 庆阳 745000；4.中石化河南油田分公司 勘探开发研究院，河南 南阳 473132；5.中石油塔里木油田分公司 天然气事业部，新疆 库尔勒 841000）

1 引言

研究岩石微观孔隙结构能深刻揭示储集层的内部结构，加深对油气储集层的认识，同时微观孔隙结构是决定描述储层多相流过程中诸多宏观性质（毛管压力、相对渗透率等）的根本。目前，研究岩石孔隙结构的主要途径是用测定毛管压力的方法来确定孔隙喉道特征（定量角度）以及电镜扫描、铸体薄片（定性角度）。测取毛管压力曲线的方法主要有半渗透隔板法、离心法和压汞法。压汞法由于测定快速、准确，并且所加压力较高，便于对更小的孔隙进行测量，同时对岩样的形状、大小要求不严，因此，成为目前测定储集岩毛管压力的主要手段。压汞法测取的毛管压力主要用于求取孔隙结构特征值，常用的求取孔隙结构特征值方法有图解法、矩法。

（1）图解法：Chilingar等[1]曾提出碳酸盐岩的孔隙喉道大小（φ 标度）遵从正态分布，可以用图解法求解孔隙结构特征值。应用孔隙大小累积分布曲线上某些特征点数据，可以求取均值、分选系数、歪度等。

（2）矩法[2]：实际岩石的孔隙喉道分布是错综复杂的，在自然界中几乎没有单一的标准正态分布的孔隙结构。基于自然界岩石孔隙结构的复杂性，可以采用地质统计理论中的混合经验分布来描述岩样中的孔隙结构。根据地质统计理论，将储集岩的孔隙分布看成是在成岩及后生作用过程中几种成因所造成的孔隙分布的组合。这样，可以把储集岩的孔隙分布在统计学上使用地质混合经验分布的数字特征[2]。对于一个混合分布，使用矩法来求取孔隙分布的特征值更为合理，即用原点矩求均值、二阶中心距求方差等。除上述矩法参数外，还可以由它们进一步组合出其他类似参数，但因为它们都来源于原始测量数据，当不同岩样之间的同类参数相互对比时，其结果和上述矩法参数的对比所得到的结论不会有较大的原则差别，因为都是用于评价储集层的孔隙结构，实际应用中意义并不大。

以上两种求取孔隙结构特征值方法都是选取孔隙半径进行直接计算：其中图解法是选取几个具有代表性的点进行算术平均，所以为了提高精度，可以加密选样点；矩法可以看成是面积积分的饱和度加权平均化处理。由于选取点不同，计算结果便不同；而且图解法使用条件严格同时误差较大，矩法的缺陷在于饱和度范围受限于试验。由于这两种方法都存在不足，本文采用积分法求取特征值，即对孔隙半径进行积分求取特征值。考虑到面积积分的平均化处理可以转换成积分形式，这样只需得到孔隙半径的数学表达式，即毛管压力数学模型，就能积分求解孔隙结构特征值。原海涵[3]对积分法和矩法之间的理论关系作了证明，因而用积分法解释孔隙结构值具有坚实的理论基础。运用积分法来求取孔隙结构特征值更接近实际值，因为它采用积分的形式将整个饱和度范围进行平均求值，同时不需要考虑孔隙大小分布情形。

关于毛管压力数学模型的探究主要有Corey[4]、van Genuchten等[5]、Thomeer[6]和Donaldson等[7]、Tiab等[8]的模型，以上几种数学模型都存在一定的不足，比如表达式不形象[4－5]，参数的确定上繁琐[6－8]。文献[6－8]在分析大量试验数据的基础上，得出毛管压力曲线形态成双曲线：文献[6]得出了毛管压力与汞体积的近似双曲线模型，但毛管压力曲线一般指的是毛管压力与饱和度之间的关系曲线；文献[7－8]推导了毛管压力三参数双曲线模型：

式中：Pc为毛管压力（MPa）；Sw为湿相饱和度（%）；A、B、C 均为常数。因此，典型的毛管压力曲线成双曲线形态，下面推导更加形象准确的毛管压力幂函数数学模型。

2 积分法求取岩石孔隙结构特征值

2.1 幂函数模型推导

由于毛管压力曲线的开始陡段并不代表非湿相真正进入岩芯，此时非湿相饱和度的增加大多是由于岩样表面凹凸不平或切开较大孔隙引起的，非湿相是在排驱压力下才开始真正地进入岩芯[2，9]。所以，开始陡段的数据点可忽略（根据经验，忽略的饱和度范围在6%以内都可接受），探究数学模型时可以不用考虑这部分。本次研究选用的具有代表性岩芯取自中原油田，其孔隙度为25.51%，渗透率为230.74×10-3μm2，其压汞试验原始数据如图1(a)所示；将汞饱和度和毛管压力同时取对数，重新作图（见图1(b)），曲线成双曲线形态（其实质是反比例函数）。

考虑到反比例函数的两个边界条件：非湿相汞饱和度为0，即湿相饱和度 Sw取值为1时，毛管压力为排驱压力Pd（边界条件1）；非湿相汞饱和度为1，即湿相饱和度 Sw=0时，毛管压力为无穷大（边界条件2）。结合反比例函数的基本表达式，得出毛管压力反比例函数模型为

考虑到毛管压力曲线形态的复杂性，对反比例函数模型引入常数D，从而建立了毛管压力幂函数模型：

式中：Pd为排驱压力（MPa）；D为常数。

但在试验中，Sw=0的极值情况是极少出现的，即存在最小湿相饱和度 Smin，即湿相流体开始流动时的最小饱和度，主要原因是岩石本身的孔隙结构（如存在不连通的死孔隙）以及试验仪器的最大压力限制[9]。因此，毛管压力是在湿相饱和度为 Smin时达到最大（可以看成是无穷大），边界条件2应改为Sw=Smin，引入转换饱和度 Si：

图1 毛管压力曲线Fig.1 Capillary pressure curves

采用最小二乘法求取参数，可得

式中： Si为转换饱和度；Smin为最小湿相饱和度（%）；n为测试点个数。

将图1(a)压汞试验数据代入式（6），求出参数D=1.067。排驱压力Pd值按照Dullien[10]的做法确定为0.073 4 MPa。因此，毛管压力幂函数模型为

将推导的幂函数模型计算曲线作图如图1(c)所示，发现试验数据点与模型计算曲线能很好地拟合，没拟合上的4个点饱和度在6%以内，数学模型中可忽略这4个点。为了验证本文推导的毛管压力幂函数模型式（5），选取了具有代表性的Sami[11]图3（试样No 95）中的数据，通过上面的方法推导模型并确定参数，将模型计算曲线与试验数据进行对比，见图2。由于Sami采用隔板法测取气-盐水毛管压力，试验压力较低，图2(a)只显示了幂函数模型的一部分，因此，将其绘制在直角坐标系图2(b)中能更好的显示。其中：Sami图3（试样No 95）毛管压力幂函数模型为

图2 Sami数据(试样No.95)Fig.2 Sami's data(Sample No.95)

由推导出的毛管压力幂函数模型式（5），结合毛管压力定义，得到毛管半径r 与转换饱和度Si的关系式为

式中：σ为界面张力（mN/m；θ为接触角（°）。

2.2 积分法求解

用积分法求取孔隙结构特征值的关键在于确定毛管半径的数学表达式。有了公式，就可以仅用试验数据点确定整个岩样的半径表达式，进而求取特征值[3]。

（1）均值[9]ar

孔喉大小总平均数的量度。矩法求取平均半径表达式中分子表示对孔隙半径与饱和度差值面积求和，可以直接转换成孔隙半径对整个饱和度范围（0～100%）连续积分；分母表示对饱和度差值进行求和，积分法可以考虑到整个饱和度范围，因此，平均毛管半径可用如下积分式表示：

式中：ri为孔隙半径（μm）；ΔSi为间隔饱和度（%）。

由图1(c)可知，在湿相饱和度（0～ Smin）范围内，毛管压力Pc可以看成是无限大，此时的孔隙半径为0；在湿相饱和度（Smin～100%）范围内，毛管压力Pc为表达式（5），导出的毛管半径表达式为式（9），因此，式（10）实质为

代入式（9），求出平均半径ra为

（2）分选系数[9]Sp

孔喉大小的标准偏差，是分散程度的量度。根据均值积分形式算法，分选系数 Sp可由下式求出：

同理，式（13）主要计算过程为

计算结果为

（3）歪度[9]Skp

描述分布不对称的量度。根据均值积分形式算法，求取歪度的积分式为

同理，式（16）主要计算过程为

计算结果为

3 结果对比与分析

根据图1岩样压汞试验数据绘制出孔隙大小分布曲线（判定为类似正态分布）和孔隙大小累计分布曲线，选取特征点代入公式[1]，即可求取孔隙结构特征值（表1）。矩法计算稍显复杂，涉及分组问题，具体参照文献[2]方法；由前面推导的毛管半径数学表达式可以得出参数值，代入积分法求取特征值公式，求取结果见表1。

表1 3种方法求取结果对比表（φ 标度）Table 1 Comparison of the calculated results of the three methods(φ value)

根据同样的方法，Sami[11]图3（试样No.95）的特征值求解结果也列入表1。

由于图解法是选取其中的几个具有代表性数据点（某些点还需要外推求解，更加大了误差，如表1求取的均值为负数，明显错误；这也正凸显了图解的巨大误差性）进行算术平均，用有限的数据点代表整个范围的半径本身就不妥；因此，一般认为，图解法求取误差最大。考虑到矩法求解的饱和度范围受到限制，致使非湿相饱和度最大值无法达到100%，如矩法求取均值的表达式，见表2，由于分母饱和度取值比实际偏小，其求解结果比实际值偏大，因而求解的特征值只表示试验流体能注入的孔隙空间量度，而不是全孔隙量度，这说明矩法存在缺陷。积分法是由矩法转换而成，其采用积分的形式将整个饱和度范围进行平均求值，更加接近真实值。因而积分法准确性高，算法严谨。

表2 图解法[1，9]和矩法[2]求取公式Table 2 Formulations of graphical method[1，9]and rectangular method[2]

4 结论

（1）图解法要求岩石孔隙大小满足正态分布，使用条件严格，同时其求取的特征值是选取的几个具有代表性数据点，误差较大。

（2）矩法求解的饱和度范围受到限制，求解的特征值只表示试验流体能注入的孔隙空间量度，而不是全孔隙量度。

（3）推导出毛管压力幂函数数学模型，即Pc=，采用最小二乘法求取参数，并用公开发表的文章数据对本文推导的模型进行了验证。

（4）详细介绍了积分法求取岩石孔隙结构特征值的基本思想。对比3种求取孔隙结构特征值方法，得出积分法算法严谨；因为积分法是采用积分的形式将整个饱和度范围进行平均求值，同时不必考虑孔隙大小分布情形，更具有普遍意义。

[1]CHILINGAR G V，MANNON R W，RISKE H H.Oil and gas production from carbonate rocks[M].New York:American Elsevier Publishing Company，1972.

[2]罗蛰潭，王允诚.油气储集层的孔隙结构[M].北京:科学出版社，1986:114－118.

[3]原海涵.毛管理论在测井解释中的应用：毛管电动力学与多孔性岩石[M].北京:石油工业出版社，1995:11－14.

[4]COREY A T.The interrelation between gas and oil relative permeabilities[J].Producers Monthly，1954，19(1):38－41.

[5]VAN GENUCHTEN，TH M.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J].Soil Science Society of America Journal，1980，44(5):892－898.

[6]THOMEER J H M.Introduction of a pore geometrical factor defined by the capillary pressure curve[J].Journal of Petroleum Technology，1960，12(3):73－77.

[7]DONALDSON E C，NINA EWALL，BALJIT SINGH.Characteristics of capillary pressure curves[J].Journal of Petroleum Science and Engineering，1991(6):249－261.

[8]TIAB D，DONALDSON E C.Petrophysics[M].New York:Elsevier Inc.，2012:351－356.

[9]何更生.油层物理[M].北京:石油工业出版社，1994.

[10]DULLIEN F A L.Porous media fluid transport and pore structure[M].New York:Academic Press，Inc.，1979.

[11]SAMI M ABOUJAFAR，MOHAMED A AMARA.A comparison between capillary and electrical properties of rock samples obtained at ambient conditions and reservoir conditions[C]//The North Africa Technical Conference &Exhibition.Cairo:Society of Petroleum Engineers，2013.