CN105427383A - 一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法 - Google Patents

Abstract

本发明公开了图像处理技术领域中一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法，其主要步骤包括：基于CT技术扫描真实岩心并进行三维重建；对三维岩心进行孔喉分割，将所有孔喉进行编号排序，统计孔喉几何参数信息；利用数值实验得到四边形最大内角和形状因子的关系；根据凹凸性、形状因子确定四边形最大内角；从该角顶点出发做对角线，将最大内角分为两部分并且确定对角线长度；建立截面非线性特征方程，利用牛顿迭代法求解四边形边长；检验特征方程是否有正解且符合物理意义，否则重新构造，直至所有孔喉截面全部构造完毕。该方法充分考虑了多孔介质截面的凹凸性、水力半径的等价性，提高了孔隙网络模型表征真实多孔介质的精确性，从而能准确预测多孔介质渗流参数。

Description

一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法

技术领域

本发明涉及图像处理技术领域，尤其涉及一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法，适用于基于孔隙网络模型进行的多孔介质重建或微观渗流模拟。

背景技术

多孔介质孔隙空间的形状和连通性极不规则，拓扑结构复杂，其微观结构和物理特性决定了许多宏观渗流性质。为了深入了解渗流机制和规律，人们通常从孔隙水平甚至更微观的水平入手研究孔隙间的复杂渗流现象。孔隙网络模型作为现在常用的孔隙级建模方法，具有可重复、可控制、计算速度快的优点，定量研究渗流规律简便易行。

利用基于岩心CT切片可以对多孔介质进行三维重建，较完善的保留真实孔喉的截面形状、拓扑结构等信息。进而可以对三维岩心做适当简化，提取其拓扑信息，并将多孔介质抽象为具有理想几何形状的孔隙空间，建立孔隙网络模型。孔隙网络模型由喉道及其相连的孔隙体构成，每个孔喉都有固定的截面形状，最常用的形状截面是C-T-S设置，即圆形、三角形、正方形(PatzekTW,SilinDB.Shapefactorandhydraulicconductanceinnoncircularcapillaries:I.One-phasecreepingflow.Journalofcolloidandinterfacescience,2001,236(2):295-304.)。这种经典的截面形状设置保证了截面形状因子(周长的平方与面积的比)相等，但是只能利用凸多边形进行形状等价，不能精确表征多孔介质中大量存在的凹形截面。同时该经典的截面设置方案不能完全保证截面的周长和面积与真实孔喉截面相等，造成水力半径(面积与周长的比)不准确。C-T-S截面设置作为一种简单近似，可能丢失大量角隅信息，影响渗流机制的进行和渗流参数的计算。

发明内容

本发明的目的是提高利用孔隙网络模型表征真实多孔介质时的精确性，从而快速、准确地预测多孔介质渗流参数。为了达到以上目的本发明提供了一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法，该方法充分考虑了多孔介质截面的凹凸性、水力半径的等价性。

本发明技术方案具体步骤如下：

(1)利用CT成像技术对岩心进行扫描，然后基于CT切片图像利用移动立方体法对岩心进行三维重建。

(2)利用细化算法对三维数字岩心进行孔喉分割，将所有孔喉编号并排序，依次为1，2，3···N_max，同时初始化N使得N＝1；统计每个孔喉截面周长P、面积A、形状因子G及孔喉半径r_in。

(3)利用数值实验随机构建包括凹四边形及凸四边形在内的M个四边形，M取值一般为5000；画出M个四边形最大内角和形状因子的散点关系图，其轮廓曲线称之为截面最大内角极值曲线，拟合曲线得到函数表达式。

(4)采用无因次水力半径H判断第N个截面的凹凸性，根据最大内角极值曲线表达式确定一个四边形最大内角β₄。

(5)利用一条从角β₄顶点出发的对角线l₀将β₄分为角α₁和角α₂两部分，同时根据公式确定该对角线长度。

(6)建立以四条边长为未知数的截面非线性特征方程。利用牛顿迭代法求解四边形参数，其中每条边长的迭代初值选为截面周长的四分之一，迭代结束条件为迭代前后绝对值之差小于10^-5或者迭代次数大于1000。

(7)检验截面非线性特征方程的解，若方程无正解或内角和不等于360°，返回步骤(5)。若方程有符合物理意义的解且N<N_max，执行N＝N+1，并转至步骤(4)，否则构造结束。

其中所述步骤(3)中的数值实验方法为：

在直角坐标系中以原点为圆心构造一个单位圆，在单位圆范围内的四个象限分别随机选取一点，依次连接各个点，组成一个四边形。

所述步骤(4)中判断截面凹凸性的方法为：

H＝r_h/r_in,

其中r_h为水力半径，r_in是孔隙半径。如果H≥0.5，则该截面形状是凸边形，否则为凹边形。

所述步骤(4)中的四边形最大内角公式为：

其中β_4min表示四边形形状因子为G时最大内角的最小值，β_4max表示四边形形状因子为G时最大内角的最大值。

所述步骤(5)中的对角线长度确定公式为：

l₀＝5A/P+k(0.5P-5A/P),k∈(0,1)

其中k为符合均匀分布的随机数。

所述步骤(6)中的截面非线性特征方程为：

$\left\{\begin{array}{c}{l}_1+l_2+l_3+l_4-p=0\\ \left(0.5l_0{\mathrm{sin\&alpha;}}_2\right)l_3+\left(0.5l_0{\mathrm{sin\&alpha;}}_1\right)l_4-A=0\\ l_1^2-l_4^2+\left(2l_0{\mathrm{cos\&alpha;}}_1\right)l_4-{l_0}^2=0\\ l_2^2-l_3^2+\left(2l_0{\mathrm{cos\&alpha;}}_2\right)l_3-{l_0}^2=0\end{array}\right..$

其中l₁、l₂、l₃、l₄分别表示四边形四条边的边长。

本发明具有以下有益效果及优点：

(1)构造过程考虑了真实岩心中的凹形截面，使得多孔介质表征更准确。

(2)所构造截面不仅保证了与真实岩心截面形状因子等价，而且保证了水力半径等价。

(3)精细的孔喉截面构造方法使得孔隙网络模型模拟得到的渗流参数更加精确。

附图说明

图1为本发明的步骤流程图。

图2为CT切片、孔隙岩石分割结果及三维重建结果。

图3为四边形数值实验构造方法及形状因子与最大内角的散点关系图。

图4为凸形及凹形孔喉截面示意图。

图5为截面构造示意图。

图6为H值不同时的截面构造结果。

具体实施方式

结合附图及实施例对本发明作进一步说明：

如图1所示，一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法，其步骤如下：

(1)利用CT成像技术对岩心进行扫描，基于图像分割技术和移动立方体法对真实岩心进行三维重建。图2依次为CT切片、孔隙岩石分割结果及三维重建结果。

(2)利用细化算法对三维数字岩心进行孔喉分割，将所有孔喉编号并排序，依次为1，2，3···N_max，同时初始化N使得N＝1；统计每个孔喉截面周长P、面积A、形状因子G及孔喉半径r_in。

(3)如图3(a)所示，在直角坐标系中以原点为圆心构造一个单位圆，在单位圆范围内的四个象限分别随机选取一点A、B、C、D，依次连接各个点可组成一个四边形。如此往复，构建5000个四边形，画出最大内角和形状因子的散点关系图如图3(b)所示。在形状因子G一定的条件下，最大内角存在取值范围，其最大值和最小值对应的散点组成的轮廓曲线称之为截面最大内角极值曲线，曲线函数拟合结果为：

$\left\{\begin{array}{c}{\mathrm{\&beta;}}_{4\min }=-7645.4G^2+10.341G+119.9\\ {\mathrm{\&beta;}}_{4\max }=6\&times;{10}^7{G}^4-9\&times;{10}^6{G}^3+303284G^2-3535.7G+274.32\end{array}\right.$

其中G为形状因子，β_4min表示四边形形状因子为G时最大内角的最小值，β_4max表示四边形形状因子为G时最大内角的最大值。

(4)采用无因次水力半径H判断第N个截面的凹凸性：

H＝r_h/r_in,

其中r_h为水力半径，r_in是孔隙半径。如果H≥0.5，则该截面形状是凸边形，否则为凹边形。图4为2个孔喉截面示例，其孔隙半径为图中的最大内切圆半径。根据计算图4(a)及图4(b)的无因次水力半径分别为0.554及0.425，说明这2个截面分别为凸边形及凹边形。进一步根据最大内角极值曲线表达式确定一个四边形最大内角β₄：

该公式保证了凹边形截面的最大内角一定大于180°，凸边形截面的最大内角一定小于180°。

(5)如图5所示，利用一条从角β₄顶点出发的对角线l₀将β₄分为角α₁和角α₂两部分，同时在合理的范围内随机确定对角线长度：

l₀＝5A/P+k(0.5P-5A/P),k∈(0,1)

(6)现已知四边形周长、面积、最大内角及从该角出发的对角线长度及对角线与两边的夹角。利用正弦定理推论及余弦定理，建立以四条边长为未知数的截面非线性特征方程如下：

$\left\{\begin{array}{c}{l}_1+l_2+l_3+l_4-p=0\\ \left(0.5l_0{\mathrm{sin\&alpha;}}_2\right)l_3+\left(0.5l_0{\mathrm{sin\&alpha;}}_1\right)l_4-A=0\\ l_1^2-l_4^2+\left(2l_0{\mathrm{cos\&alpha;}}_1\right)l_4-{l_0}^2=0\\ l_2^2-l_3^2+\left(2l_0{\mathrm{cos\&alpha;}}_2\right)l_3-{l_0}^2=0\end{array}\right..$

然后利用牛顿迭代法求解该截面非线性特征方程组，迭代初值为(0.25P，0.25P，0.25P，0.25P)，迭代结束条件为迭代前后绝对值之差小于10^-5或者迭代次数大于1000。以G＝0.04、A/P＝10(无量纲)为例，在H值不同时，构造的图形如图6所示。图6(a)对应H＝0.5的情况,此时为正凸多边形；图6(b)对应H>0.5的情况，此时为凸多边形；图6(c)对应H<0.5的情况,此时为凹多边形。其关键参数取值在表1中列出。

表1图6中满足G＝0.04,A/p＝10的求解结果

(7)由于边长必须为正值，因此如果方程的根是负值时是不符合物理意义的；另外由于迭代结束条件的限制，正解在某些情况下可能由于精度不够导致不能拼接成四边形，此时可以分别求出四边形内角，看内角和是否为360°。求解方程后进一步检验截面非线性特征方程的解，若方程无正解或内角和不等于360°，返回步骤(5)。若方程有符合物理意义的解且N<N_max，执行N＝N+1，并转至步骤(4)，否则构造结束。

Claims (6)

1.一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法，其特征在于以下步骤：

(1)利用CT成像技术对岩心进行扫描，然后基于CT切片图像利用移动立方体法对岩心进行三维重建。

(2)利用细化算法对三维数字岩心进行孔喉分割，将所有孔喉编号并排序，依次为1，2，3···N_max，同时初始化N使得N＝1；统计每个孔喉截面周长P、面积A、形状因子G及孔喉半径r_in。

(3)利用数值实验随机构建包括凹四边形及凸四边形在内的M个四边形，M取值一般为5000；画出M个四边形最大内角和形状因子的散点关系图，其轮廓曲线称之为截面最大内角极值曲线，拟合曲线得到函数表达式。

(4)采用无因次水力半径H判断第N个截面的凹凸性，根据最大内角极值曲线表达式确定一个四边形最大内角β₄。

(5)利用一条从角β₄顶点出发的对角线l₀将β₄分为角α₁和角α₂两部分，同时根据公式确定该对角线长度。

(6)建立以四条边长为未知数的截面非线性特征方程。利用牛顿迭代法求解四边形参数，其中每条边长的迭代初值选为截面周长的四分之一，迭代结束条件为迭代前后绝对值之差小于10^-5或者迭代次数大于1000。

(7)检验截面非线性特征方程的解，若方程无正解或内角和不等于360°，返回步骤(5)。若方程有符合物理意义的解且N<N_max，执行N＝N+1，并转至步骤(4)，否则构造结束。

2.如权利要求1所述的一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法，其所述步骤(3)中的数值实验，其特征在于：

在直角坐标系中以原点为圆心构造一个单位圆，在单位圆范围内的四个象限分别随机选取一点，依次连接各个点，组成一个四边形。

3.如权利要求1所述的一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法，其所述步骤(4)中采用无因次水力半径H判断第N个截面的凹凸性，其特征在于：

H＝r_h/r_in,

其中r_h为水力半径，r_in是孔隙半径。如果H≥0.5，则该截面形状是凸边形，否则为凹边形。

4.如权利要求1所述的一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法，其所述步骤(4)中的四边形最大内角β₄，其特征在于由以下公式决定：

其中β_4min表示四边形形状因子为G时最大内角的最小值，β_4max表示四边形形状因子为G时最大内角的最大值。

5.如权利要求1所述的一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法，其所述步骤(5)中的对角线l₀，其特征在于该对角线长度确定公式为：

l₀＝5A/P+k(0.5P-5A/P),k∈(0,1)

其中k为符合均匀分布的随机数。

6.如权利要求1所述的一种考虑凹凸性的岩石孔隙网络模型的孔喉截面构造方法，其所述步骤(6)中的截面非线性特征方程，其特征在于满足以下形式：

$\left\{\begin{array}{c}{l}_1+l_2+l_3+l_4-p=0\\ \left(0.5l_0{\mathrm{sin\&alpha;}}_2\right)l_3+\left(0.5l_0{\mathrm{sin\&alpha;}}_1\right)l_4-A=0\\ l_1^2-l_4^2+\left(2l_0{\mathrm{cos\&alpha;}}_1\right)l_4-{l_0}^2=0\\ l_2^2-l_3^2+\left(2l_0{\mathrm{cos\&alpha;}}_2\right)l_3-{l_0}^2=0\end{array}\right..$

其中l₁、l₂、l₃、l₄分别表示四边形四条边的边长。