CN110210176B - 一种花岗岩微观裂缝结构识别和微观渗流分析方法 - Google Patents

Abstract

本发明公开了一种花岗岩微观裂缝结构识别和微观渗流分析方法，读取花岗岩单偏光显微图片M；对花岗岩单偏光显微图片M依次进行超像素化分割和滤波处理获得花岗岩显微数字图像M2；对花岗岩显微数字图像M2进行二值化处理，获得花岗岩显微二值化图像M3；分析获得花岗岩微观裂缝系统的流体流动速度分布图、流体压力分布图、等效渗透率；分析获得花岗岩微观裂缝系统‑花岗岩基质系统中的流体流动速度分布图和流体压力分布图。本发明实用简单，成本较低，能够重复利用花岗岩单偏光显微图片M，对花岗岩微观裂缝系统渗流规律进行可视化研究，并能够容易地计算花岗岩微观裂缝系统的等效渗透率。

Description

一种花岗岩微观裂缝结构识别和微观渗流分析方法

技术领域

本发明属于致密花岗岩中裂缝表征及渗流分析技术领域，具体涉及一种花岗岩微观裂缝结构识别和微观渗流分析方法。

背景技术

缝洞不发育的花岗岩是致密储层，孔隙度和渗透率极低。但花岗岩中油气藏的发现，表明受长期风化淋滤作用，局部花岗岩体内部微裂缝、溶洞和溶孔都很发育，使其具有较好的储集性能。并且其储层物性的好坏与否取决于岩性、风化程度、变质程度及花岗岩成份和热液溶蚀洞穴的数量等。

对于致密花岗岩而言，其基质渗透率极低，含非常少的基质溶蚀孔隙，渗流的通道主要是微细观的裂缝系统。致密花岗岩微观裂缝的分布特征制约着宏观裂缝的扩展方向，影响裂缝的交叉、转向等，因此对花岗岩内微观裂缝结构的精细刻画有利于揭示宏观裂缝的扩展规律。

但限于工业CT分辨率低且成本高，扫描电镜SEM等常规手段无法观察到微观裂缝的全貌，压汞法因所需泵注压力太高，使其无法获得致密花岗岩中孔隙和裂缝空间分布、孔径大小等参数。所以传统的这些方法无法深入地对微观裂缝的空间形态进行精细表征，从而很难捕捉到流体的微观渗流路径和渗流规律。而利用显微镜观测法，可以直观的揭示花岗岩中微观裂缝、孔隙及颗粒的空间分布特征。

认识和掌握裂缝系统中的渗流规律对花岗岩裂缝发育储层的流体开发具有重要意义。花岗岩中微观裂缝系统的渗流机理和渗流规律决定了宏观上的裂缝系统渗流特性。开展花岗岩微观裂缝系统渗流规律研究对认识花岗岩裂缝系统渗流特性具有关键作用。但目前关于花岗岩中微观裂缝系统渗流规律的研究较少，需要进一步的研究。

发明内容

为了克服现有技术方案的不足，本发明提供了一种花岗岩微观裂缝结构识别和微观渗流分析方法，有效实现背景技术中提到的花岗岩微观裂缝表征和含微裂缝系统的花岗岩渗流规律研究。

本发明采用的技术方案为：

一种花岗岩微观裂缝结构识别和微观渗流分析方法，包括以下步骤：

步骤1、读取花岗岩单偏光显微图片M；

步骤2、对花岗岩单偏光显微图片M依次进行超像素化分割和滤波处理获得花岗岩显微数字图像M2；

步骤3、对步骤2获得的花岗岩显微数字图像M2进行二值化处理，得到含花岗岩微观裂缝系统和花岗岩基质系统的花岗岩显微二值化图像M3；

步骤4、将花岗岩显微二值化图像M3导入到COMSOL Multiphysics软件中；

步骤5、基于COMSOL Multiphysics软件中的层流模块对花岗岩显微二值化图像M3的花岗岩微观裂缝系统进行渗流数值试验，得到花岗岩微观裂缝系统中的流场分布图和等效渗透率K₀，花岗岩微观裂缝系统中的流场分布图包括花岗岩微观裂缝系统的流体流动速度分布图和流体压力分布图；

步骤6、基于COMSOL Multiphysics软件中的Brinkman方程模块对花岗岩显微二值化图像M3中的花岗岩微观裂缝系统-花岗岩基质系统在设定压差驱动下的渗流特征进行数值模拟试验，得到流体渗流规律，流体渗流规律包括花岗岩微观裂缝系统-花岗岩基质系统中的流体流动速度分布图和流体压力分布图。

如上所述的步骤5中的花岗岩微观裂缝系统的流体流动速度分布图和流体压力分布图通过以下步骤获得：

对导入到COMSOL Multiphysics软件中的花岗岩显微二值化图像M3做布尔差集，得到花岗岩微观裂缝系统的几何模型，采用三角形单元对该几何模型进行网格划分，并在裂缝边界处进行网格加密处理。然后对层流模块的物理场条件进行设定，即设置裂缝系统入口端和出口端的压差，设置裂缝壁面设为无滑移边界，通过求解层流场中流体流动的控制方程，得到花岗岩微观裂缝系统的流体流动速度和流体压力，从而生成花岗岩微观裂缝系统中的流体流动速度分布图和流体压力分布图；

层流场中流体流动的控制方程为：

式中，u-流体流动速度；I-单位矩阵；p-流体压力；ρ-流体密度；μ-流体粘度；T-矩阵的转置；

和·分别表示张量的梯度算子和散度符号。

如上所述的步骤5中的等效渗透率K₀通过以下步骤获得：

其中，A-流体过流面积；ΔL-花岗岩微观裂缝系统特征长度；Δp-裂缝系统入口端和出口端的压差；Q-花岗岩微观裂缝系统的流体流量。

如上所述的步骤6中的花岗岩微观裂缝系统-花岗岩基质系统中的流体流动速度分布图和流体压力分布图通过以下步骤获得：

利用COMSOL Multiphysics软件中的内置图像函数im1(x,y)对花岗岩显微二值化图像M3进行处理获得图像M4，并利用公式1-0.99×(1-im1(x,y))和K₀/(100×(1-im1(x,y))+0.1)分别获得图像M4的各点的孔隙度φ和渗透率k，其中，x，y分别为图像M4中各点的坐标，将图像M4导入到COMSOL Multiphysics软件中作为花岗岩微观裂缝系统-花岗岩基质系统的花岗岩微观二维几何模型，设定流体密度ρ和流体粘度μ分别为1000kg/m³和1×10^{- 3}Pa·s，利用Brinkman方程模块，对物理场进行条件设定，给定花岗岩微观二维几何模型左右两端压差，裂缝壁面边界设定为无滑移边界，通过求解Brinkman方程模块中流体运动的控制方程，得到基于Brinkman方程模块的花岗岩微观裂缝系统-花岗岩基质系统中的流体流动速度分布图和流体压力分布图，

Brinkman方程模块中流体运动的控制方程如下

式中，p-流体压力；u-流体流动速度；T-矩阵的转置；

和·分别表示张量的梯度算子和散度符号；ρ-流体密度；μ-流体粘度。

与现有技术相比，本发明具有以下优点及有益效果：

1、本发明能够方便的重构花岗岩显微图像，快速识别花岗岩显微图像中的花岗岩微观裂缝系统-花岗岩基质系统，得到可用于微观渗流分析的花岗岩显微二值化图像，提高图像处理效率；

2、本发明能够快速方便的计算得到花岗岩微观裂缝系统的等效渗透率；

3、本发明简便实用，成本较低，能够重复利用显微镜下得到的花岗岩单偏光显微图片开展数值渗流试验，对花岗岩微观裂缝系统渗流规律进行可视化研究。

附图说明

图1为本发明的流程图；

图2为本发明中使用的花岗岩单偏光显微图片M；

图3为本发明中利用Matlab进行超像素化分割后的超像素化后的图像M1；

图4为本发明中利用Matlab进行二值化处理后的花岗岩显微二值化图像M3；

图5为COMSOL Multiphysics软件中花岗岩微观裂缝系统网格剖分；

图6为COMSOL Multiphysics软件中利用层流模块计算得到的花岗岩微观裂缝系统中的流体流动速度分布图；

图7为COMOSL中利用层流模块计算得到的花岗岩微观裂缝系统中的流体压力分布图。

图8为COMSOL中利用内置图像函数im1(x,y)将花岗岩显微二值化图像M3转化后的图像M4；

图9为COMSOL中求解Brinkman方程得到的花岗岩微观裂缝系统-花岗岩基质系统中的流体流动速度分布图；

图10为COMSOL中求解Brinkman方程得到的花岗岩微观裂缝系统-花岗岩基质系统中的流体压力分布图。

具体实施方式

为了更清楚的说明本发明所采用的技术手段和方法，下面结合具体实施例对本发明进行清楚、完整的介绍。显而易见的，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明专利保护范围。

本发明的具体实施步骤如下：

一种花岗岩微观裂缝结构识别和微观渗流分析方法，包括以下步骤：

步骤1、读取花岗岩单偏光显微图片M，具体为：

在Maltab中调用函数imread()读取花岗岩单偏光显微图片M；

步骤2、对花岗岩单偏光显微图片M依次进行超像素化分割和滤波处理获得花岗岩显微数字图像M2；

超像素化分割包括以下步骤：首先在Maltab中调用图像处理的函数superpixels(M,spix_numel,‘NumIterations’,spix_numIter,‘Method’,spix_meth{1})，对花岗岩单偏光显微照片M(长度为L，宽度为W)进行超像素化，超像素数spix_numel设置为5000，迭代次数NumIterations选择spix_numIter＝10，计算超像素的算法Method选择spix_meth{1}＝slic0，slic0的算法可以实现第一次线性迭代之后的图像紧密性，然后调用函数boundarymask(N,conn_2d(2))寻找分割区域边界获得超像素化后的图像M1，其中N为函数superpixels的矩阵返回值，conn_2d(2)＝8，是用来通过接触的边缘和角点建立一个像素和邻近像素区的连通性。

超像素分割处理通过将花岗岩单偏光显微图片M中相邻的具有相同性质或类似性质的像素或区域归并到一起，可将花岗岩单偏光显微图像M中具有基质属性和裂缝属性的像素集合到一起。

滤波处理包括以下步骤：在Matlab中调用函数imgaussfilt(M1,σ)对超像素化后的图像M1进行高斯滤波处理，其中高斯平滑的标准差设置为σ＝8，来获得高斯滤波后的花岗岩显微数字图像M2。

步骤3、对步骤2获得的花岗岩显微数字图像M2进行二值化处理，得到含花岗岩微观裂缝系统和花岗岩基质系统的花岗岩显微二值化图像M3；

二值化处理包括以下步骤：

调用rgb2gray()函数，将步骤2获得的花岗岩显微数字图像M2(RGB图)进行灰度图转换获得花岗岩灰度图Igray，接着调用函数imbinarize(Igray,0.99)对花岗岩灰度图进行二值化处理，得到含花岗岩微观裂缝系统和花岗岩基质系统的花岗岩显微二值化图像M3，使得裂隙部分被处理成为白色，基质部分被处理成为黑色；

步骤4、将花岗岩显微二值化图像M3导入到COMSOL Multiphysics软件中；

步骤5、基于COMSOL Multiphysics软件中的层流模块对花岗岩显微二值化图像M3的花岗岩微观裂缝系统进行渗流数值试验，得到花岗岩微观裂缝系统中的流场分布图和等效渗透率K₀。花岗岩微观裂缝系统中的流场分布图包括花岗岩微观裂缝系统的流体流动速度分布图和流体压力分布图。

具体包括以下步骤：

对导入到COMSOL Multiphysics中的花岗岩显微二值化图像M3做布尔差集，得到花岗岩微观裂缝系统的几何模型，采用三角形单元对该几何模型进行网格划分，并在裂缝边界处进行网格加密处理。然后对层流模块的物理场条件进行设定，即设置裂缝系统入口端和出口端的压差为0.715Pa，裂缝壁面设为无滑移边界，通过求解层流场中流体流动的控制方程，得到花岗岩微观裂缝系统的流体流动速度(u)和流体压力(p)，从而生成花岗岩微观裂缝系统中的流体流动速度分布图和流体压力分布图。

层流场中流体流动的控制方程为：

式中，u-流体流动速度，单位为m/s；I-单位矩阵；p-流体压力，单位为Pa；ρ-流体密度，单位为kg/m³；μ-流体粘度，单位为Pa·s；T-矩阵的转置；

表示张量的梯度算子；·表示散度符号。

对花岗岩微观裂缝系统出口边界的流体流动速度积分，得到流过单位厚度(1m)的花岗岩微观裂缝系统的流体流量Q，本实施例中，Q＝1.12×10^-9m³/s。

根据以下公式获得花岗岩微观裂缝系统的等效渗透率K₀：

其中，A-流体过流面积(单位为m²)，ΔL-花岗岩微观裂缝系统特征长度(单位为m)，Δp-裂缝系统入口端和出口端的压差(单位为Pa)。

本实施例中，层流模块模拟所用流体为水，流体密度ρ被设定为1000kg/m³，流体粘度μ为1×10^-3Pa·s，流体过流面积A按单位厚度(1m)×花岗岩显微二值化图像M3的宽度(W)＝5.9017×10^-4m²，花岗岩微观裂缝系统特征长度ΔL取为花岗岩显微二值化图像M3的长度＝7.846×10^-4m，裂缝系统入口端和出口端的压差Δp设定为0.715Pa，因此通过计算可以得到花岗岩微观裂缝系统的等效渗透率K₀(单位为m²)为：

步骤6、基于COMSOL Multiphysics软件中的Brinkman方程模块对花岗岩显微二值化图像M3中的花岗岩微观裂缝系统-花岗岩基质系统在一定压差驱动下的渗流特征进行数值模拟试验，得到流体渗流规律，流体渗流规律包括花岗岩微观裂缝系统-花岗岩基质系统中的流体流动速度分布图和流体压力分布图。

具体包括如下步骤：

基于Matlab处理得到的花岗岩显微二值化图像M3，利用COMSOL Multiphysics中的内置图像函数im1(x,y)对花岗岩显微二值化图像M3进行处理获得图像M4，并利用公式1-0.99×(1-im1(x,y))和K₀/(100×(1-im1(x,y))+0.1)分别获得图像M4的各点的孔隙度φ和渗透率k，其中，x，y分别为图像M4中各点的坐标。将图像M4导入到COMSOL Multiphysics中作为花岗岩微观裂缝系统-花岗岩基质系统的花岗岩微观二维几何模型，设定流体密度ρ和流体粘度μ分别为1000kg/m³和1×10^-3Pa·s，利用Brinkman方程模块，对物理场进行条件设定，给定花岗岩微观二维几何模型左右两端压差为0.715Pa，裂缝壁面边界仍然设定为无滑移边界，通过求解Brinkman方程模块中流体运动的控制方程，得到基于Brinkman方程模块的花岗岩微观裂缝系统-花岗岩基质系统中的流体流动速度(u)分布图和流体压力(p)分布图。

Brinkman方程模块中流体运动的控制方程如下

式中，p-流体压力，单位为Pa；u-流体流动速度，单位为m/s；T-矩阵的转置；

和·分别表示张量的梯度算子和散度符号；ρ-流体密度，单位为kg/m³；μ-流体粘度，单位为Pa·s。

本文中所描述的具体实施例仅是对本发明作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或替代，但不会偏离本发明的精髓或者超越所附权利要求书定义的范围。

Claims (1)

1.一种花岗岩微观裂缝结构识别和微观渗流分析方法，其特征在于，包括以下步骤：

步骤1、读取花岗岩单偏光显微图片M；

步骤2、对花岗岩单偏光显微图片M依次进行超像素化分割和滤波处理获得花岗岩显微数字图像M2；

步骤3、对步骤2获得的花岗岩显微数字图像M2进行二值化处理，得到含花岗岩微观裂缝系统和花岗岩基质系统的花岗岩显微二值化图像M3；

步骤4、将花岗岩显微二值化图像M3导入到COMSOL Multiphysics软件中；

步骤5、基于COMSOL Multiphysics软件中的层流模块对花岗岩显微二值化图像M3的花岗岩微观裂缝系统进行渗流数值试验，得到花岗岩微观裂缝系统中的流场分布图和等效渗透率K₀，花岗岩微观裂缝系统中的流场分布图包括花岗岩微观裂缝系统的流体流动速度分布图和流体压力分布图；

步骤6、基于COMSOL Multiphysics软件中的Brinkman方程模块对花岗岩显微二值化图像M3中的花岗岩微观裂缝系统-花岗岩基质系统在设定压差驱动下的渗流特征进行数值模拟试验，得到流体渗流规律，流体渗流规律包括花岗岩微观裂缝系统-花岗岩基质系统中的流体流动速度分布图和流体压力分布图，

所述的步骤5中的花岗岩微观裂缝系统的流体流动速度分布图和流体压力分布图通过以下步骤获得：

对导入到COMSOL Multiphysics软件中的花岗岩显微二值化图像M3做布尔差集，得到花岗岩微观裂缝系统的几何模型，采用三角形单元对该几何模型进行网格划分，并在裂缝边界处进行网格加密处理，然后对层流模块的物理场条件进行设定，即设置裂缝系统入口端和出口端的压差，设置裂缝壁面为无滑移边界，通过求解层流场中流体流动的控制方程，得到花岗岩微观裂缝系统的流体流动速度和流体压力，从而生成花岗岩微观裂缝系统中的流体流动速度分布图和流体压力分布图；

层流场中流体流动的控制方程为：

式中，u-流体流动速度；I-单位矩阵；p-流体压力；ρ-流体密度；μ-流体粘度；T-矩阵的转置；

和·分别表示张量的梯度算子和散度符号，

所述的步骤5中的等效渗透率K₀通过以下步骤获得：

其中，A-流体过流面积；ΔL-花岗岩微观裂缝系统特征长度；Δp-裂缝系统入口端和出口端的压差；Q-花岗岩微观裂缝系统的流体流量，

所述的步骤6中的花岗岩微观裂缝系统-花岗岩基质系统中的流体流动速度分布图和流体压力分布图通过以下步骤获得：

利用COMSOL Multiphysics软件中的内置图像函数im1(x,y)对花岗岩显微二值化图像M3进行处理获得图像M4，并利用公式1-0.99×(1-im1(x,y))和K₀/(100×(1-im1(x,y))+0.1)分别获得图像M4的各点的孔隙度φ和渗透率k，其中，x，y分别为图像M4中各点的坐标，将图像M4导入到COMSOL Multiphysics软件中作为花岗岩微观裂缝系统-花岗岩基质系统的花岗岩微观二维几何模型，设定流体密度ρ和流体粘度μ分别为1000kg/m³和1×10^-3Pa·s，利用Brinkman方程模块，对物理场进行条件设定，给定花岗岩微观二维几何模型左右两端压差，裂缝壁面边界设定为无滑移边界，通过求解Brinkman方程模块中流体运动的控制方程，得到基于Brinkman方程模块的花岗岩微观裂缝系统-花岗岩基质系统中的流体流动速度分布图和流体压力分布图，

Brinkman方程模块中流体运动的控制方程如下

式中，p-流体压力；u-流体流动速度；T-矩阵的转置；

和·分别表示张量的梯度算子和散度符号；ρ-流体密度；μ-流体粘度。

Images