基于截断高斯分布函数的岩心多峰分布孔隙结构表征方法

1.本发明涉及岩心孔隙结构建模技术领域，尤其涉及一种基于截断高斯分布函数的岩心多峰分布孔隙结构表征方法。


背景技术：

2.岩心孔隙结构复杂导致很难使用常规的几何方法对岩心孔隙结构进行描述。常规单分形模型得到的岩心孔隙分布为递减曲线，难以表示煤岩、页岩、碳酸盐岩等岩心的多峰分布特征。多重分形模型使用的维数过多，分形维数难以通过实验准确测试，并且函数通常不连续。所以，使用分形模型描述岩心孔隙结构特征具有一定的局限性。
3.根据岩心核磁实验可知，在半对数曲线上，页岩和碳酸盐岩毛管直径t2谱一般为双峰特征(branco et al 2017)，煤岩t2谱为双峰或者三峰特征(golsanami et al 2019)，砂岩t2谱一般为单峰或者双峰特征(xin et al 2022)，假设单个孔隙空间形态为孔径λ的圆柱体，t2弛豫时间可以转化为毛管直径：
4.λ＝10-3
t2ρ25.其中，t2为横向弛豫时间，ms；ρ2为横向表面弛豫率，μm/ms；λ为孔隙孔径， nm。
6.核磁测试曲线可以转化为孔隙分布曲线，曲线中x轴为毛管直径，y轴为孔隙占比。在整个孔隙尺寸范围内，毛管数量与毛管直径的分形关系式可由下式表示(yu and cheng 2002)：
[0007][0008]
其中，λ
min
—最小孔隙直径；λ
max
—最大孔隙直径。
[0009]
使用以上分形模型则可以计算得到的页岩孔隙分布和核磁测试结果转换的孔隙分布(ρ2＝0.02μm/ms)，可以看出，分形模型无法表征出孔隙的多峰分布特征，便需要建立一种用以描述岩心多峰分布特征的表征方法。


技术实现要素：

[0010]
基于上述现有技术无法表征出孔隙的多峰分布特征的缺点，本发明提供了基于截断高斯分布函数的岩心多峰分布孔隙结构表征方法，具体技术方案如下：
[0011]
基于截断高斯分布函数的岩心多峰分布孔隙结构表征方法，包括以下步骤：
[0012]
建立基于截断高斯分布叠加的多峰分布岩心孔隙结构模型；
[0013]
建立岩心孔隙多峰分布孔隙度模型；
[0014]
建立岩心孔隙多峰分布渗透率模型；
[0015]
使用单目标增强遗传算法计算模型参数。
[0016]
具体的，所述建立基于截断高斯分布叠加的多峰分布岩心孔隙结构模型包括以下子步骤：
[0017]
步骤11：对孔径尺寸对数归一化，计算式为：
[0018][0019]
其中，λs为归一化孔径，取值范围为[0,1]，无因次；
[0020]
步骤12：将多组截断高斯分布概率密度函数的线性叠加，得到毛细管数量函数fn(λ)，计算式为：
[0021][0022]
其中，μ为数学期望；σ为标准差；φ
trc
为截断高斯分布的概率密度函数； wk为.....。
[0023]
具体的，所述建立岩心孔隙多峰分布孔隙度模型，包括以下子步骤：
[0024]
步骤21：计算孔径为λ的毛管束的体积v
λ
，计算式为：
[0025]vλ
＝n
λaλ
l
a,λ
[0026]
其中，v
λ
为孔径为λ的毛管束总管内体积，m3；n
λ
为孔径为λ的毛细管数量，无因次；l
a,λ
为孔径为λ的毛细管实际长度，m；
[0027]
步骤22：计算孔径为λ的毛细管截面积，计算式为：
[0028][0029]
其中，λ为单根毛细管孔径，m；a
λ
为孔径为λ毛细管的截面积，m2；f
λ
为孔径为λ毛细管形状因子，无因次；
[0030]
步骤23：计算以孔径为λ的圆作为基准形状，质心重叠且各边与质心最近距离为λ的规则多边形面积与基准形状面积的比值，记为形状因子f
λ
，计算式为：
[0031][0032]
步骤24：计算毛细管实际长度l
a,λ
，计算式为：
[0033]
l
a,λ
＝lbδ
λ
[0034]
其中，lb为页岩基质表观长度，m；δ
λ
为孔径为λ的毛细管迂曲度，无因次；
[0035]
步骤25：采用迂曲度分形维数来计算迂曲度，计算式为：
[0036][0037]
其中，df为迂曲度分形维数，无因次；
[0038]
步骤26：计算符合n
peak
峰分布的孔径为λ的毛管束总体积v
λ
，计算式为：
[0039][0040]
其中，na为毛管总数；
[0041]
步骤27：计算符合n
peak
峰分布的孔径为λ的毛细管集合，毛管束总孔隙体积v
p
，计算式为：
[0042][0043]
步骤28：计算孔隙度，计算式为：
[0044][0045]
其中，r为岩心直径，m。
[0046]
具体的，所述建立岩心孔隙多峰分布渗透率模型包括以下子步骤：
[0047]
步骤31：计算单根毛管中的流量，计算式为：
[0048][0049]
步骤32：计算通过岩心的总流量，计算式为：
[0050][0051]
步骤33：计算岩心多峰分布渗透率k，计算式为：
[0052][0053]
具体的，所述使用单目标增强遗传算法计算模型参数，包含以下子步骤：
[0054]
步骤41：计算目标函数为孔隙分布的实测值与计算值的均方根误差，计算式为：
[0055][0056]
其中，n为实验测试数据的数量，无因次；f
nreal
(i)为第i个孔隙分布的实测值，m3；f
ncalc
(i)为第i个孔隙分布的计算值，m3；ε
fλ
为孔隙分布的均方根误差， m3；
[0057]
步骤42：使用单目标增强遗传算法对未知模型参数进行求解，所述模型参数包括毛管束总量，na；迂曲分形维数，df；左峰期望，μ1；右峰期望，μ2；左峰方差，σ1；右峰方差，σ2；左峰权重，ω1和右峰权重，ω2。
[0058]
具体的，步骤42所述的使用单目标增强遗传算法未知模型参数进行求解，还包括以下子步骤：
[0059]
步骤421：初始化单目标增强遗传算法；
[0060]
步骤422：生成参数方案；
[0061]
步骤423：将参数方案分为k个模型参数组；
[0062]
步骤424：对各个模型参数组求解f
nk
(λs)；
[0063]
步骤425：对各个模型参数组计算目标函数为孔隙分布的实测值与计算值的均方根误差；
[0064]
步骤426：将各个模型参数组的目标函数为孔隙分布的实测值与计算值的均方根
误差返回到单目标增强遗传算法；
[0065]
步骤427：判断是否达到截止条件，若是，则输出最优参数组；若否，则返回步骤422。
[0066]
本发明的有益效果：
[0067]
本发明使用高斯截断分布的统计学方法，基于岩心毛管束基本假设，通过两个高斯分布叠加的方法建立了岩心孔隙多峰分布函数，并推导了岩心孔隙度、渗透率模型。进一步使用遗传算法，结合岩心核磁实验测试数据，进行多参数拟合，建立分布函数中的相关参数求解方法。本发明建立的方法可以准确表征岩心孔隙体积的多峰分布特征，为描述多孔介质孔隙结构提供新思路。进一步为描述岩心气水两相流动规律，研究产水气井生产规律提供理论支撑。
附图说明
[0068]
图1所示是本发明采用单目标增强遗传算法进行毛管束模型参数求取的计算流程图；
[0069]
图2所示为本发明#a1岩心核磁共振测试测试结果图；
[0070]
图3所示为本发明#a2岩心核磁共振测试测试结果图；
[0071]
图4所示为本发明#a3岩心核磁共振测试测试结果图；
[0072]
图5所示为本发明的#a1岩样孔隙分布拟合结果图；
[0073]
图6所示为本发明的#a1岩样累积孔隙体积对比图；
[0074]
图7所示为本发明的#a2岩样孔隙分布拟合结果；
[0075]
图8所示为本发明的#a2岩样累积孔隙体积对比图；
[0076]
图9所示为本发明的#a3岩样孔隙分布拟合结果；
[0077]
图10所示为本发明的#a3岩样累积孔隙体积对比图；
[0078]
图11所示为本发明迂曲度对拟合结果的影响图。
具体实施方式
[0079]
为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图说明本发明的具体实施方式。
[0080]
基于截断高斯分布函数的岩心多峰分布孔隙结构表征方法，包括以下步骤：
[0081]
建立基于截断高斯分布叠加的多峰分布岩心孔隙结构模型；
[0082]
建立岩心孔隙多峰分布孔隙度模型；
[0083]
建立岩心孔隙多峰分布渗透率模型；
[0084]
使用单目标增强遗传算法计算模型参数。
[0085]
具体的，所述建立基于截断高斯分布叠加的多峰分布岩心孔隙结构模型包括以下子步骤：
[0086]
步骤11：对孔径尺寸对数归一化，计算式为：
[0087][0088]
其中，λs为归一化孔径，取值范围为[0,1]，无因次；
[0089]
步骤12：将多组截断高斯分布概率密度函数的线性叠加，得到毛细管数量函数fn(λ)，计算式为：
[0090][0091]
其中，μ为数学期望；σ为标准差；φ
trc
为截断高斯分布的概率密度函数； wk为.....。
[0092]
具体的，所述建立岩心孔隙多峰分布孔隙度模型，包括以下子步骤：
[0093]
步骤21：计算孔径为λ的毛管束的体积v
λ
，计算式为：
[0094]vλ
＝n
λaλ
l
a,λ
[0095]
其中，v
λ
为孔径为λ的毛管束总管内体积，m3；n
λ
为孔径为λ的毛细管数量，无因次；l
a,λ
为孔径为λ的毛细管实际长度，m；
[0096]
步骤22：计算孔径为λ的毛细管截面积，计算式为：
[0097][0098]
其中，λ为单根毛细管孔径，m；a
λ
为孔径为λ毛细管的截面积，m2；f
λ
为孔径为λ毛细管形状因子，无因次；
[0099]
步骤23：计算以孔径为λ的圆作为基准形状，质心重叠且各边与质心最近距离为λ的规则多边形面积与基准形状面积的比值，记为形状因子f
λ
，计算式为：
[0100][0101]
步骤24：计算毛细管实际长度l
a,λ
，计算式为：
[0102]
l
a,λ
＝lbδ
λ
[0103]
其中，lb为页岩基质表观长度，m；δ
λ
为孔径为λ的毛细管迂曲度，无因次；
[0104]
步骤25：采用迂曲度分形维数来计算迂曲度，计算式为：
[0105][0106]
其中，df为迂曲度分形维数，无因次；
[0107]
步骤26：计算符合n
peak
峰分布的孔径为λ的毛管束总体积v
λ
，计算式为：
[0108][0109]
其中，na为毛管总数；
[0110]
步骤27：计算符合n
peak
峰分布的孔径为λ的毛细管集合，毛管束总孔隙体积v
p
，计算式为：
[0111][0112]
步骤28：计算孔隙度，计算式为：
[0113][0114]
其中，r为岩心直径，m。
[0115]
具体的，所述建立岩心孔隙多峰分布渗透率模型包括以下子步骤：
[0116]
步骤31：计算单根毛管中的流量，计算式为：
[0117][0118]
步骤32：计算通过岩心的总流量，计算式为：
[0119][0120]
步骤33：计算岩心多峰分布渗透率k，计算式为：
[0121][0122]
具体的，所述使用单目标增强遗传算法计算模型参数，包含以下子步骤：
[0123]
步骤41：计算目标函数为孔隙分布的实测值与计算值的均方根误差，计算式为：
[0124][0125]
其中，n为实验测试数据的数量，无因次；f
nreal
(i)为第i个孔隙分布的实测值，m3；f
ncalc
(i)为第i个孔隙分布的计算值，m3；ε
fλ
为孔隙分布的均方根误差，m3；
[0126]
步骤42：使用单目标增强遗传算法对未知模型参数进行求解，所述模型参数包括毛管束总量，na；迂曲分形维数，df；左峰期望，μ1；右峰期望，μ2；左峰方差，σ1；右峰方差，σ2；左峰权重，ω1和右峰权重，ω2。
[0127]
如图1所示，具体的，步骤42所述的使用单目标增强遗传算法未知模型参数进行求解，还包括以下子步骤：
[0128]
步骤421：初始化单目标增强遗传算法；
[0129]
步骤422：生成参数方案；
[0130]
步骤423：将参数方案分为k个模型参数组；
[0131]
步骤424：对各个模型参数组求解f
nk
(λs)；
[0132]
步骤425：对各个模型参数组计算目标函数为孔隙分布的实测值与计算值的均方根误差；
[0133]
步骤426：将各个模型参数组的目标函数为孔隙分布的实测值与计算值的均方根误差返回到单目标增强遗传算法；
[0134]
步骤427：判断是否达到截止条件，若是，则输出最优参数组；若否，则返回步骤422。
[0135]
实例运用本发明技术方案与分析如下：
[0136]
#a1号岩样核磁测试结果
[0137]
如图2所示，岩心表观体积33.9cm3，饱和吸水质量1.3g，称重法孔隙度3.8％；核磁共振累积孔隙体积1.10cm3，核磁孔隙度为3.2％，主峰峰值对应的孔径约为 5.7nm，主峰峰体分布孔径范围为0.5～141.1nm；次峰峰值对应的孔径约为 670.5nm，次峰峰体分布孔径范围为214.5～865.8nm；孔径约162.2nm处的分布较小，将主、次峰的峰体间隔开。
[0138]
#a2号岩样核磁测试结果
[0139]
如图3所示，岩心表观体积26.0cm3，饱和吸水质量1.2g，称重法孔隙度4.6％；核磁共振累积孔隙体积1.09cm3，核磁孔隙度为4.2％；主峰峰值对应的孔径约为 4.3nm，主峰峰体分布孔径范围为0.4～30.4nm；次峰峰值对应的孔径约为135.2nm，次峰峰体分布孔径范围为53.1～214.5nm；孔径约40.2nm处的分布较小，将主、次峰的峰体间隔开。
[0140]
#a3号岩样核磁测试结果
[0141]
如图4所示，岩心表观体积37.3cm3，饱和吸水质量1.3g，称重法孔隙度4.3％；核磁共振累积孔隙体积1.23cm3，核磁孔隙度为3.3％；主峰峰值对应的孔径约为 3.2nm，主峰峰体分布孔径范围为0.4～46.2nm；次峰峰值对应的孔径约为162.2nm，次峰峰体分布孔径范围为61.1～865.8nm；孔径约53.1nm处的分布较小，将主、次峰的峰体间隔开。
[0142]
岩心孔隙多峰分布模型参数计算
[0143]
#a1号岩样
[0144]
通过对实验数据进行拟合，考虑、不考虑迂曲度的孔隙分布拟合均方根误差分别为1.91
×
10-3
和1.92
×
10-3
；总孔隙体积真实值为1.095cm3，考虑、不考虑迂曲度的总孔隙体积计算值分别为1.122cm3、1.118cm3，相对误差分别为 2.43％和2.02％；模型参数拟合结果见表1，孔隙分布计算结果见图5，累积孔隙体积对比见图6。
[0145][0146]
表1#a1模型结构参数拟合结果
[0147]
#a2号岩样
[0148]
通过对实验数据进行拟合，考虑、不考虑毛细管迂曲度的孔隙分布拟合均方根误差分别为1.31
×
10-3
和1.59
×
10-3
；总孔隙体积真实值1.087cm3，考虑、不考虑迂曲度的总孔隙体积计算值分别为1.135cm3、1.159cm3，相对误差分别为4.46％和6.70％；累积孔隙体积模型参数拟合结果见表2，孔隙分布计算结果见图7，累积孔隙体积对比见图8。
[0149][0150][0151]
表2#a2岩样参数拟合结果
[0152]
#a3号岩样
[0153]
通过对实验数据进行拟合，考虑、不考虑毛细管迂曲度的孔隙分布拟合均方根误差分别为1.16
×
10-3
和1.29
×
10-3
；总孔隙体积真实值1.231cm3，考虑、不考虑迂曲度的总孔隙体积计算值分别为1.282cm3、1.293cm3，相对误差分别为4.20％和5.03％；模型参数拟合结果见表3，孔隙分布计算结果见图9，累积孔隙体积对比见图10。
[0154][0155]
表3#a3岩样参数拟合结果
[0156]
1.1.1结果讨论
[0157]
(1)模型适应性
[0158]
孔径为λ的毛细管集合的v
λ
满足双峰分布的条件下，采用了考虑、不考虑迂曲度的两种策略进行了三组样本的孔隙分布拟合，拟合精度分析见表4。
[0159]
[0160][0161]
表4双峰毛管束结构模型拟合精度分析表
[0162]
从计算结果上看，孔隙分布v
λ
的均方根误差在1.16～1.92
×
10-3
之间，总孔隙体积误差在2.02％～6.70％之间，证明了本文建立的双峰毛管束结构模型能够良好表征页岩基质孔隙分布的双峰分布特征，且能准确定量计算不同毛管的孔隙体积。
[0163]
(2)迂曲度分析
[0164]
为评估迂曲度对孔隙分布拟合精度的影响，将考虑、不考虑迂曲度的拟合误差进行了对比，见图11。
[0165]
考虑迂曲度和不考虑迂曲度三组岩样孔隙分布的平均rmse分别为1.46
×ꢀ
10-3
、1.6
×
10-3
，考虑迂曲度使平均均方根误差降低了9.59％，平均孔隙体积的相对误差降低了32.58％，迂曲分形维数df有利于提升拟合效果，但对与n
λ
相关的na、ω1、μ1、σ1、ω2、μ2、σ2等7个参数的拟合结果有较大影响。考虑迂曲度的迂曲度分形维数分别为1.136、1.025和1.049，与分形毛管束模型通常采用的[1.5～1.7]取值范围(friesen and mikula 1987；angulo and alvarado etal 1992；)有一定的差异，导致这一现象的原因可能是：未考虑多峰分布时，迂曲度较大才能更好的体现岩心非均质性；考虑了多峰分布后，岩心非均质性主要由多峰分布体现，迂曲度相对较小。总而言之，迂曲度及迂曲分形维数df都有利于提高双峰毛管束结构模型的准确性。
[0166]
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。