微米尺度平面孔隙网络内气液驱替图像的饱和度获取方法与流程

本发明属于油气藏强化开采领域，适用于二氧化碳地质封存领域的微尺度平面孔隙多孔介质中气液两相流动实验，具体涉及一种微米尺度平面孔隙网络内气液两相驱替图像的饱和度获取方法。

背景技术：

平面孔隙网络中的气液两相驱替研究是获取天然多孔介质中气液两相驱替机理的主要手段。天然多孔介质中的气液两相驱替过程发生在众多的工业过程中，例如致密油气藏的气驱油、水驱气、二氧化碳的地质封存等过程。毛管力与流体饱和度的关系曲线，即毛管力曲线，是油气储层中不混溶两相渗流的基本特征。流体相对渗透率与其饱和度的关系曲线，即相对渗透率曲线，是研究不混溶流体渗流的基础。不混溶的流体在多孔介质中驱替渗流时，每种流体会占据一定的孔隙体积，多孔介质孔隙中任意一点p处m相流体的饱和度定义为p点周围的特征体v^*内的m相流体所占据的孔隙体积的百分比，即

其中是特征体v^*内的孔隙体积，是特征体v^*内的m相流体的体积，特征体内各相流体饱和度之和为1。若取特征体为整个多孔介质，则多孔介质内m流体的饱和度sm为多孔介质孔隙vp内m相流体所占据的孔隙体积的百分比，即

其中vm是整个多孔介质孔隙内m相流体的体积，多孔介质内各相流体饱和度之和为1。毛管力曲线和相对渗透率曲线广泛应用于油气田开发参数计算、动态分析以及油藏数值模拟等，是油气田开发中的两个重要关系曲线。获得毛管力曲线和相对渗透率曲线的关键参数是确定各相流体的饱和度。

在微尺度平面孔隙网络中进行两相不混溶流体驱替实验是研究多孔介质渗流机理的常用方法。微米尺度平面孔隙网络内的渗流是小流量驱替工况，特别是进行准静态驱替实验时，此时难以通过实时测量各相流体的流量来确定孔隙内各相流体的体积，进一步地难以确定各相流体饱和度。油藏岩心样品中进行两相不混溶流体驱替实验时，往往先要测量多孔介质的孔隙度，再通过测量流经多孔介质的入口和出口流量之差来确定多孔介质内的流体的饱和度，该方法对多孔介质入口与出口处的流量传感器精度要求较高。以上方法不适用微米尺度的平面孔隙网络内的流体饱和度的测量。

技术实现要素：

本发明的目的在于提供微米尺度平面孔隙网络内气液驱替图像的饱和度获取方法，以克服现有技术的不足，本发明能够从平面孔隙网络气液两相驱替图像中识别出孔隙网络的固体区域、各流体流动区域、相界面，同时测量出各相流体流动区域的面积，各相流体的饱和度以及驱替流量。

为达到上述目的，本发明采用如下技术方案：

微米尺度平面孔隙网络内气液驱替图像的饱和度获取方法，包括以下步骤：

步骤一，获取平面孔隙网络模型图像，对该图像进行高斯滤波除去噪声后转化为灰度图像，对灰度图像进行坐标划分，并利用阈值分割方法将灰度图像转化像素尺寸为宽lx和高ly的黑白二值图像；

步骤二，裁剪黑白二值图像中入口与出口区域，保留中间孔隙网络区域，生成校验图像；

步骤三，统计校验图像中白色区域像素数，获得白色区域总面积，计算流动区域总面积atotal；

步骤四，获取平面孔隙网络气相驱替形态图像，裁剪平面孔隙网络气相驱替形态图像入口与出口区域，并进行灰度阈值分割生成黑白二值待测图像；

步骤五，统计黑白二值待测图像中白色区域像素数，获得白色区域总面积，计算气相和液相流动区域面积；

步骤六，对比校验图像与黑白二值待测图像，获得气液两相流动区域驱替形态图像，并求取气相饱和度sg和液相饱和度sl以及驱替流量qg和ql。

进一步的，步骤一中对灰度图像进行坐标划分：以平面孔隙网络模型灰度图像的左上角为坐标原点，建立直角坐标系，将图像分割为若干区块。

进一步的，求取每个区块的灰度阈值，将其转化为黑白二值图像，黑白二值图像中任意坐标(x,y)处的灰度值为i(x,y)，黑色区域的灰度i(x,y)＝255，白色区域的灰度i(x,y)＝0，坐标(x,y)处的像素值记为p(x,y)；

进一步的，裁剪黑白二值图像入口与出口区域的具体过程如下：

(1)将黑白二值图像的左边界上所有坐标点的像素值p(0,yn)与其相邻的坐标点的像素值p(0,yn-1)和p(0,yn+1)分别进行“与”运算，搜索运算结果为0的坐标点ai(0,ya)和bi(0,yb)即为二值图像左边界上像素灰度变化点；

(2)将平面孔隙网络右边界上所有坐标点的像素值p(lx,ym)与其相邻的坐标点p(lx,ym-1)和p(lx,ym+1)的像素点进行“与”运算，搜索运算结果为0的坐标点co(lx,yc)、do(lx,yd)即为二值图像右边界上像素灰度变化点；

(3)分别沿直线y＝ya和直线y＝yb向x轴正方向按步骤(1)的方法搜索获取左边界上像素灰度变化点a(xa,ya)、b(xb,yb)；

(4)分别沿直线y＝yc和直线y＝yd向x轴负方向按步骤(2)的方法搜索获取右边界上像素灰度变化点c(xc,yc)、d(xd,yd)；

(5)以直线a-b和直线c-d为界限，裁剪黑白二值图像，去除平面孔隙网络入口和出口区域，只保留中间孔隙网络区域，获得像素尺寸为宽(xd-xa)和高ly或宽(xc-xb)和高ly的校验图像。

进一步的，步骤三中求取流动区域总面积的具体过程如下：

(1)从步骤二中生成的校验图像的左上角开始计数，搜索整个校验图像区域，白色区域的像素数为整个校验图像区域所有的像素值之和；

(2)按照校验图像像素尺寸与长度尺寸的比例，将像素数换算为白色区域总面积，即为平面孔隙网络中流动区域总面积atotal。

进一步的，步骤四中裁剪气相驱替形态图像：(1)以平面孔隙网络中气相驱替形态图像的左上角为坐标原点，建立直角坐标系；以步骤二中获得的左边界上像素灰度变化点a(xa,ya)、b(xb,yb)建立直线a-b，以右边界上像素灰度变化点c(xc,yc)、d(xd,yd)建立直线c-d；

(2)以直线a-b和直线c-d为界限，裁剪二值图像，去除气相驱替形态图像的入口和出口区域，只保留中间区域，获得像素尺寸为宽(xd-xa)和高ly或宽(xc-xb)和高ly的驱替形态图像。

进一步的，将裁剪完成的驱替形态图像分割为若干区块，取每个区块的灰度阈值，将其转化为黑白二值图像，生成驱替形态待测图像，黑白二值图像中任意坐标(x,y)处的灰度值为i(x,y)，黑色区域的灰度i(x,y)＝255，白色区域的灰度i(x,y)＝0，坐标(x,y)处的像素值记为p(x,y)；

进一步的，步骤五中获取气相和液相流动区域面积的具体过程如下：

(1)从步骤四生成的黑白二值待测图像的左上角开始计数，搜索整个黑白二值待测图像区域，白色区域的像素数为整个黑白二值待测图像区域所有的像素值之和；

(2)按照黑白二值待测图像像素尺寸与长度尺寸的比例，将白色区域像素数换算为白色区域总面积，即为气液两相驱替中气相流动面积ag，流动区域总面积atotal减去气相流动面积即得液相流动面积al：

al＝atotal-ag(5)。

进一步的，步骤六中获得气液两相流动区域驱替形态图像、气液两相饱和度以及驱替流量的具体过程如下：

(1)将黑白二值待测图像与校验图像进行异或操作，得到可同时显示气相和液相流动区域的驱替形态图像，对驱替形态图像进行图形开操作，消除噪音，使边界平滑，得到同时显示气相和液相驱替形态图像；

(2)将待测图像中气相流动面积ag除以校验图像中流动区域总面积atotal，即得到气相饱和度sg；

(3)将待测图像中液相流动面积al除以校验图像中流动区域总面积atotal，即得到液相饱和度sl；

(4)δt时间内的相饱和度的变化即为δt时间内气液相平均驱替流量qg和ql；

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

本发明微米尺度平面孔隙网络内气液驱替图像的饱和度获取方法，获取平面孔隙网络模型图像，对该图像进行高斯滤波除去噪声、坐标划分得到黑白二值图像，然后对黑白二值图像中入口与出口区域裁剪统计校验图像中白色区域像素数，获得白色区域总面积，同时平面孔隙网络气相驱替形态图像入口与出口区域，并进行灰度阈值分割生成黑白二值待测图像；统计黑白二值待测图像中白色区域像素数获得白色区域总面积，对比校验图像与黑白二值待测图像，获得气液两相流动区域驱替形态图像从而计算得到气相饱和度和液相饱和度以及驱替流量，本方法通过识别流动区域总面积以及气相、液相占据区域，计算气相和液相饱和度，并获得气相和液相驱替形态图像和驱替流量，测量结果稳定，鲁棒性强，本方法可高效、精确地处理大量图像数据，特别是对低驱替速度下的饱和度分析和驱替流量计算，更加便捷、有效。

附图说明

图1是本发明一个实施方式的饱和度获取方法的流程示意图。

图2是平面孔隙网络模型示意图。

图3是黑白二值图像以及裁剪过程示意图。

图4是校验图像示意图。

图5是气相驱替形态图像以及裁剪过程示意图。

图6是气相驱替形态图像裁剪结果示意图。

图7是气相驱替形态待测图像示意图。

图8是同时显示气相和液相流动区域驱替形态图像示意图。

具体实施方式

下面结合附图对本发明做进一步详细描述：

气液驱替图像特征具体如下：气液驱替过程，显微镜光源采用的是反射光观察方式，利用显微镜反射光拍摄平面孔隙网络和气液两相驱替形态时，可以获得清晰的气固界面以及气液界面，但难以获得清晰的液固界面。气相与液相和固相具有很高的对比度，高灰度区域为平面孔隙网络固体和液相，低灰度区域为气相占据区域，如图2和图5所示。图2中圆孔连线区域(白色区域)为可流动区域，被驱替实验开始之前，被空气占据或抽为真空，其余区域(灰色区域)为平面孔隙网络固体。图5中灰色区域为平面孔隙网络固体和液相占据区域，白色区域被空气或其他气体占据。利用这一特点，本方法通过识别流动区域总面积以及气相、液相占据区域，计算气相和液相饱和度，并获得气相和液相驱替形态图像，可高效、精确地处理大量图像数据，特别是对低驱替速度下的饱和度分析，更加便捷、有效。

平面孔隙网络模型图像和气相驱替形态图像裁剪：

显微镜拍摄的平面孔隙网络模型和气相驱替形态图像包含了孔隙网络之外的入口和出口区域，这对饱和对分析带来的干扰，需要将其裁剪。如图3和图5所示，以直线a-b为界限，点(0,0)、(0,ly)、(xb,ly)和(xa,0)围城的区域为入口区域，以直线c-d为界限，点(0,0)、(0,ly)、(xb,ly)和(xa,0)围城的区域为出口区域，需要将入口区域和出口区域裁剪，再求取流动区域总面积以及气液两相占据面积。

如图1所示，微米尺度平面孔隙网络内气液驱替图像的饱和度获取方法，利用显微镜反射光拍摄平面孔隙网络和气相驱替形态的高对比度图像，包括以下步骤，以获取气液两相饱和度：

步骤一，读取平面孔隙网络模型图像，对该图像进行高斯滤波除去噪声，转化为灰度图像，对灰度图像进行坐标划分，并利用阈值分割方法将灰度图像转化像素尺寸为宽lx和高ly的黑白二值图像，平面孔隙网络模型如图2所示，黑白二值图像如图3所示；

步骤二，裁剪黑白二值图像中入口与出口区域，保留中间孔隙网络区域，生成校验图像，具体校验图像如图4所示；

步骤三，统计校验图像中白色区域像素数，获得白色区域总面积，计算流动区域总面积atotal；

步骤四，读取平面孔隙网络气相驱替形态图像，裁剪入口与出口区域，并进行灰度阈值分割生成黑白二值待测图像；

步骤五，统计黑白二值待测图像中白色区域像素数，获得白色区域总面积，计算气相和液相流动区域面积；

步骤六，对比校验图像与黑白二值待测图像，获得气液两相流动区域驱替形态图像，并求取气相饱和度sg和液相饱和度sl以及驱替流量qg和ql。

步骤一中对灰度图像进行坐标划分和灰度阈值分割的具体过程如下：(1)以平面孔隙网络模型灰度图像的左上角为坐标原点，建立直角坐标系，将图像分割为若干区块；

(2)求取每个区块的灰度阈值，将其转化为黑白二值图像，黑白二值图像中任意坐标(x,y)处的灰度值为i(x,y)，黑色区域的灰度i(x,y)＝255，白色区域的灰度i(x,y)＝0，坐标(x,y)处的像素值记为p(x,y)；

步骤二中裁剪黑白二值图像入口与出口区域的具体过程如下：

(1)将黑白二值图像的左边界(直线x＝0)上所有坐标点的像素值p(0,yn)与其相邻的坐标点的像素值p(0,yn-1)和p(0,yn+1)分别进行“与”运算，搜索运算结果为0的坐标点ai(0,ya)和bi(0,yb)即为二值图像左边界上像素灰度变化点；

(2)将平面孔隙网络右边界(直线x＝lx)上所有坐标点的像素值p(lx,ym)与其相邻的坐标点p(lx,ym-1)和p(lx,ym+1)的像素点进行“与”运算，搜索运算结果为0的坐标点co(lx,yc)、do(lx,yd)即为二值图像右边界上像素灰度变化点；

(3)分别沿直线y＝ya和直线y＝yb向x轴正方向按步骤(1)的方法搜索左边界上像素灰度变化点a(xa,ya)、b(xb,yb)；

(4)分别沿直线y＝yc和直线y＝yd向x轴负方向按步骤(2)的方法搜索右边界上像素灰度变化点c(xc,yc)、d(xd,yd)；

(5)以直线a-b和直线c-d为界限，裁剪黑白二值图像，去除平面孔隙网络入口和出口区域，只保留中间孔隙网络区域，获得像素尺寸为宽(xd-xa)和高ly或宽(xc-xb)和高ly的校验图像，校验图像如图4所示。

步骤三中求取流动区域总面积的具体过程如下：

(1)从步骤二中生成的校验图像的左上角开始计数，搜索整个校验图像区域，白色区域的像素数为整个校验图像区域所有的像素值之和；

(2)按照校验图像像素尺寸与长度尺寸的比例，将像素数换算为白色区域总面积，即为平面孔隙网络中流动区域总面积atotal。

步骤四中裁剪气相驱替形态图像以及图像二值化处理的具体过程如下，气相驱替形态图像如图5所示：

(1)以平面孔隙网络中气相驱替形态图像的左上角为坐标原点，建立直角坐标系；以步骤二中获得的左边界上像素灰度变化点a(xa,ya)、b(xb,yb)建立直线a-b，以右边界上像素灰度变化点c(xc,yc)、d(xd,yd)建立直线c-d；

(2)以直线a-b和直线c-d为界限，裁剪二值图像，去除气相驱替形态图像的入口和出口区域，只保留中间区域，获得像素尺寸为宽(xd-xa)和高ly或宽(xc-xb)和高ly的驱替形态图像，驱替形态图像如图6所示；

(3)将裁剪完成的驱替形态图像分割为若干区块，取每个区块的灰度阈值，将其转化为黑白二值图像，生成驱替形态待测图像，即黑白二值待测图像，驱替形态待测图像如图7所示，黑白二值图像中任意坐标(x,y)处的灰度值为i(x,y)，黑色区域的灰度i(x,y)＝255，白色区域的灰度i(x,y)＝0，坐标(x,y)处的像素值记为p(x,y)；

步骤五中获取气相和液相流动区域面积的具体过程如下：

(1)从步骤四生成的驱替形态的黑白二值待测图像的左上角开始计数，搜索整个黑白二值待测图像区域，白色区域的像素数为整个黑白二值待测图像区域所有的像素值之和；

(2)按照黑白二值待测图像像素尺寸与长度尺寸的比例，将白色区域像素数换算为白色区域总面积，即为气液两相驱替中气相流动面积ag，步骤三获得的流动区域总面积atotal减去气相流动面积即得液相流动面积al：

al＝atotal-ag(5)

气液两相流动区域驱替形态图像如图8所示，步骤六中获得气液两相流动区域驱替形态图像、气液两相饱和度以及驱替流量的具体过程如下：

(1)将黑白二值待测图像与校验图像进行异或操作，得到可同时显示气相和液相流动区域的驱替形态图像，对驱替形态图像进行图形开操作，消除噪音，使边界平滑，得到图8的同时显示气相和液相驱替形态图像；

(2)步骤五得到待测图像中气相流动面积ag除以校验图像中流动区域总面积atotal，即得到气相饱和度sg；

(3)步骤五得到待测图像中液相流动面积al除以校验图像中流动区域总面积atotal，即得到液相饱和度sl；

(4)δt时间内的相饱和度的变化即为δt时间内气液相平均驱替流量qg和ql；

本发明测量结果稳定，鲁棒性强；利用该方法编程简单，程序运行效率高，采用常规处理器单线程处理像素2304×1728的1张校验图像和1张待测图像的图片耗时大约20秒；可并行处理，采用常规处理器10线程并行处理像素2304×1728的1张校验图像和2000张待测图像的图片耗时大约5分钟。