川西坳陷致密砂岩气储层沉积、成岩作用复杂，导致孔渗相关性较差，且孔隙结构复杂，其中蓬莱镇组、沙溪庙组致密砂岩储层孔隙空间大部分为纳米—微米级孔隙，且具有较强的非均质性^[1]。针对致密砂岩储层特征，众多学者^[2-4]开展了大量微观孔隙结构研究工作，并对其孔隙类型、孔隙形态、孔隙大小和范围、孔隙空间分布等进行了研究，但对孔渗宏观参数与微观孔隙结构评价参数的相关研究稍有不足^[5]，且渗透率的计算也多在孔隙结构分类的基础上进行^[6]。在微观孔隙结构定量评价方面，众多学者^[7-9]通过多种实验和方法证实：沉积岩孔隙结构具有分形特征，而分形理论是描述物体复杂程度、非均质性的有效手段^[10-11]。

目前致密气储层孔隙结构实验方法主要有高压压汞、扫描电镜、小角度中子散射、核磁共振、氮气吸附、微—纳米CT等^[12]，其中高压压汞、氮气吸附、核磁共振是较常用的3种实验方法。氮气吸附所表征的孔隙通常小于100 nm，该范围难以描述致密气储层的微米级孔隙^[13-14]，分形处理得到的分形维数也仅代表纳米孔的孔隙结构特征^[15-16]。核磁共振与高压压汞均能表征纳米—微米级孔隙，但是在基于核磁共振资料的分形研究中，包含了横向弛豫时间(T₂)与毛管进汞压力(p_c)或孔隙半径(r)成线性关系这一假设^[17-19]。实际上T₂与p_c关系复杂，线性关系、幂函数关系均存在^[20-22]，直接利用线性关系进行致密砂岩储层核磁分形公式推导较为不妥。本次研究通过岩心高压压汞实验，并结合物性资料，利用分形理论对致密砂岩孔隙结构非均质性进行研究。

本次实验样品选自川西坳陷蓬莱镇组、沙溪庙组致密砂岩气藏12块代表性岩心。样品孔隙度为6.32%~19.90%，平均值为12.4%，渗透率为0.006~0.726 mD，平均值为0.124 mD。蓬莱镇组岩心样品深度为1 200~1 840 m，沙溪庙组岩心样品深度为2 130~2 550 m。本次实验利用9505型压汞仪对12块样品进行压汞测试，仪器体积分辨率高于0.01 mL，最高进汞压力可达200 MPa，实验温度为22 ℃，湿度为60%。实验方法和实验数据处理均按照《岩石毛管压力曲线的测定》行业标准^[23]来执行。

基于压汞资料进行分形处理的方法较多^[24-26]。贺承祖等^[25]最早推导了储层岩石孔隙(孔喉)分布的分形公式，进而推导了毛管压力和J函数等分形几何公式。在孔喉分布分形公式基础上，葛小波等^[27]直接利用拉普拉斯方程进行毛管压力分形公式推导。吴浩等^[12]将毛管束模型与拉普拉斯方程结合进行分形公式推导，推导结果差异在于分形维数与压力-累积进汞量曲线斜率的关系略有不同。司马立强等^[26]对J函数分形公式、毛管束分形公式、热力学分形公式的适用性进行了比较。结果表明，毛管束分形公式更能表征致密砂岩储层的孔隙结构特征。因此，本次实验选用毛管束分形公式对压汞资料进行处理，推导过程如下^{[12, 25]}。

储层孔喉分布具有分形结构，孔喉数量N(> r)与r必有如下幂函数关系

式中：N为孔喉数量；r为孔喉半径，m；a为常数；D

为孔喉分形维数，其理论值为2~3。同时，利用压汞实验和毛管束模型，毛管束数量N(> r)的另一种表达形式为

式中：V_Hg为孔喉半径达到r时的累计进汞体积，m³；l为毛细管长度，m。将拉普拉斯方程式与式(1)、式(2)相结合，同时将累计进汞体积转换为累计进汞饱和度，得如下关系式

式中：S_Hg为累计进汞饱和度，%；p_c为毛管压力，MPa；b为常数。

将式(3)两边取对数，可得

由上述可知，当储层孔喉具有分形结构时，其分形维数可通过lg S_Hg与lg p_c作图，由线性拟合确定直线斜率K，进而得到D，即

在实际压汞资料处理的过程中，川西致密砂岩储层分形曲线出现明显拐点，呈两段式线性关系，每一线性段的相关系数平方均大于0.9，分段点进汞压力值分布于2.0~27.6 MPa，表现出储层孔隙结构具有2种分形特征：第一分形特征段的维数称为D_b；第二分形特征段的维数称为D_s，其加权平均分形维数为D_ave。

依据排驱压力、中值半径、孔隙度、渗透率和进汞曲线形态差异等，将样品孔隙结构类型分为4类(图 1)。Ⅰ类储层孔隙结构表现为排驱压力低、渗透率高、中值半径较大；Ⅱ类储层孔隙结构表现为排驱压力较低、渗透率较高、中值半径较大；Ⅲ类储层孔隙结构表现为排驱压力较高、渗透率较低、中值半径较小；Ⅳ类储层孔隙结构表现为排驱压力高、渗透率低、中值半径小。由表 1可见，除渗透率外，孔隙度、排驱压力、中值压力和分形维数等单一参数无法区分孔隙结构类型，表明了孔隙结构的复杂性。

图 1

Fig. 1

下载eps/tif图 图 1 孔隙结构分类 Fig. 1 Pore structure classification

表 1(Table 1)

表 1 致密砂岩样品压汞和分形参数 Table 1 Mercury injection and fractal parameters in tight sandstone 样品 孔隙度/% 渗透率/mD 排驱压力/MPa 中值压力/MPa 中值半径/μm Db Ds Dave pf/MPa 分类 X41-5 16.63 0.726 0.265 6.31 0.274 2.59 2.14 2.34 2.00 Ⅰ类 F2-17 19.9 0.305 0.678 3.96 0.372 2.46 2.11 2.32 2.75 Ⅱ类 P4-2 13.04 0.151 0.674 8.43 0.186 2.60 2.15 2.36 2.75 Ⅱ类 X10-49 11.55 0.050 0.680 13.68 0.130 2.72 2.17 2.42 4.14 Ⅱ类 X41-2 13.87 0.096 0.674 25.49 0.106 2.52 2.19 2.33 2.75 Ⅱ类 X11-43 10.25 0.047 0.673 36.75 0.133 2.67 2.18 2.36 2.74 Ⅱ类 P3-3 13.9 0.027 0.677 7.94 0.170 2.75 2.13 2.45 6.00 Ⅲ类 P4-5 7.96 0.033 1.361 38.81 0.043 2.81 2.23 2.46 5.60 Ⅲ类 X10-45 14.67 0.028 0.674 33.53 0.062 2.85 2.23 2.41 2.75 Ⅲ类 J6-2 10.37 0.012 2.054 22.23 0.090 2.56 2.19 2.37 8.30 Ⅳ类 P3-1 6.32 0.008 2.738 47.67 0.032 2.95 2.28 2.52 8.30 Ⅳ类 X14-11 10.39 0.006 2.056 85.96 0.015 2.35 2.51 2.46 27.60 Ⅳ类

下载CSV 表 1 致密砂岩样品压汞和分形参数 Table 1 Mercury injection and fractal parameters in tight sandstone

基于压汞分析的4种孔隙结构类型划分结果，选取不同孔隙结构类型的4个典型样品，绘制其分形特征曲线(图 2)。结果显示，孔隙结构越好，分形曲线越平缓，在lg p_c - lg S_Hg直角坐标系的位置更趋于左上方，表明不同孔隙结构储层的分形曲线差异明显。虽然单独的分形维数无法区分致密砂岩储层孔隙结构类型，但分形曲线形态差异可作为储层孔隙结构类型划分的直接依据^[28]。

图 2

Fig. 2

下载eps/tif图 图 2 典型样品的分形特征曲线 Fig. 2 Fractal curves of typical samples

样品F2-17和J6-2分别对应于孔隙结构特征中的Ⅱ类和Ⅳ类，2个样品的分形特征曲线均有明显的两段式分段特征(图 3)，表明孔隙空间具有2种分形性质。图 3 (a)代表大孔的分形特征，曲线分段点在0.51，图 3(b)代表小孔的分形特征，曲线分段点在0.93。文献[15]、[29]中氮气吸附分形分界点基本固定在相对压力0.5处，表明氮气吸附的大、小孔界限固定，属于绝对的大孔和小孔，而研究区样品经压汞分形处理得到的分段点位置变化范围广，大、小孔界限不固定，属于相对的大孔和小孔。

图 3

Fig. 3

下载eps/tif图 图 3 部分样品分形特征曲线 Fig. 3 Fractal curves of some samples

由式(4)和式(5)可求出相对大孔、小孔的分形维数。同一样品大孔分形维数(D_b)明显大于小孔分形维数(D_s)，且不同样品当D_b大时，其D_s未必大(参见表 1)，表明大孔、小孔具有相对独立的孔隙结构。

由表 2可见，孔隙度与D_b和D_s相关性均较低，而与D_ave具有较好的负相关性，表明孔隙度更多受孔隙结构整体分形特征的影响，D_ave越大，孔隙度越低。

表 2(Table 2)

表 2 分形维数与孔隙结构参数关系 Table 2 Relationships between fractal dimensions and pore structure parameters 拟合参数 Db Ds Dave 孔隙度 y = 108.69 x-2.275，R2 = 0.22 y = 174.87 x-3.398，R2 = 0.22 y = 36 061 x-9.163，R2 = 0.58 渗透率 — y = 1E + 6 x-21.68，R2 = 0.48 y = 3E + 14 x-41.47，R2 = 0.57 中值半径 — y = 1E + 5 x-18.36，R2 = 0.84 y = 8E + 8 x-26.11，R2 = 0.55

下载CSV 表 2 分形维数与孔隙结构参数关系 Table 2 Relationships between fractal dimensions and pore structure parameters

渗透率与D_b无相关性，而与D_s和D_ave均具有明显的负相关性，且与D_ave相关性更好，表明渗透率同样更多受孔隙结构整体分形特征影响。比较大、小孔分形维数与渗透率的关系可知，小孔分形特征是制约渗透率大小的主要因素。中值半径与D_b无相关性，而与D_s和D_ave均具有明显的负相关性，且与D_s相关性更好，相关系数平方达到0.8以上，表明中值半径受小孔与孔隙结构整体分形特征的影响。

渗透率计算模型通常包括渗透率与孔隙度的幂函数拟合、多项式拟合以及指数拟合等。利用岩心孔渗测试资料，建立渗透率与孔隙度的多种拟合关系(图 4)。

图 4

Fig. 4

下载eps/tif图 图 4 渗透率与孔隙度的多种拟合关系 Fig. 4 Several fitted curves between permeability and porosity

渗透率与孔隙度的线性、多项式拟合关系较差，相关系数平方均小于0.5，而渗透率与孔隙度的幂函数、指数拟合相关系数平方虽然达到0.5，但是无法满足储层精细评价的要求。在考虑孔隙结构对渗透率影响的基础上，众多学者^[30-33]从理论上对多孔介质渗透率的分形公式进行推导，取得了较好的效果，但是这些公式实际上包含了难以确定的迂曲度等参数，实际应用时具有一定局限性。

本次研究尝试利用多元回归的方法建立渗透率计算模型。由前述可知，分形特征对渗透率大小具有明显影响，分形特征包含D_ave与p_f等参数；渗透率与R₅₀等参数通常具有较好的相关性。通过理论分析与多次试算，最终选取D_ave，p_f与R₅₀进行渗透率多元回归计算，回归公式为

式中：K为渗透率，mD；D_ave为加权平均分形维数；p_f为分形曲线分段点处压力，MPa；R₅₀为中值半径，μm。

由图 5可见，回归公式相关系数平方(R²)大于0.9，实测渗透率与计算渗透率处于y = x的直线上，误差较小，明显优于由孔隙度拟合的渗透率模型，且避免了迂曲度等参数的确定。由此表明，基于D_ave，p_f与R₅₀的回归公式可以更准确地计算致密砂岩储层的渗透率，为致密砂岩储层渗透率计算提供了新思路。考虑R₅₀与D_s的相关性较好，实际表达式为

图 5

Fig. 5

下载eps/tif图 图 5 实测渗透率与计算渗透率交会图 Fig. 5 Cross plot between measured and calculation permeability

该表达式与李留仁等^[34]提出的分形多孔介质渗透率理论关系模型K = f(D)具有较好的一致性，本文模型仅增加p_f这一参数，增加原因在于样品孔隙结构具有2种分形特征。

(1) 川西坳陷致密砂岩储层孔隙结构复杂，依据排驱压力、中值半径、孔隙度、渗透率和进汞曲线形态差异等孔隙结构特征，致密砂岩储层可划分为Ⅰ类、Ⅱ类、Ⅲ类、Ⅳ类共4种孔隙结构类型。

(2) 样品分形曲线具有分段特征，不同样品分段点压力不同，孔隙结构越好，分形曲线越平缓，在lg p_c - lg S_Hg直角坐标系的位置更趋于左上方，利用分形曲线形态差异可区分孔隙结构类型；孔隙具有相对独立的孔隙结构，孔隙度、渗透率主要受加权平均分形维数D_ave制约，R₅₀主要受小孔分形特征制约。

(3) 基于加权平均分形维数、分段压力、中值半径的多元回归公式可以更准确地计算渗透率，回归公式与理论关系模型表达式具有较好的一致性，为致密砂岩储层渗透率计算提供了新思路。