0 引言
利用微流控模型进行水驱油实验，可以直观呈现油水两相在驱替过程中的流体流动特征及剩余油分布规律，对认识复杂油藏的开发规律具有重要的参考意义。目前大部分微流控实验仅简单地观察了多相流体在不同环境条件下的渗流现象，而对其内在机理缺乏定量化的认识[1⇓-3]。不同于实验模拟、三维模型重建等分析方法，数字图像处理技术可以更真实地定量描述流体的流动状态，这促使微流控实验技术由现象观察向定量描述发展[4⇓-6]。王川等[7]基于剩余油分布的图像特征，实现了对微观剩余油的分类识别和定量统计。吴飞鹏等[8]、高文彬等[9]利用微观模型水驱油实验，分别在不同条件下对微观剩余油进行分类和统计。也有学者利用Micro CT扫描实验进行了微观剩余油分布规律的研究[10⇓-12]，可以看到上述研究均针对静态条件下微观剩余油的状态，而非流动条件下的分布规律。
除微观剩余油分布规律外，“油-水-岩石”的三相接触角也是油水渗流的重要研究对象。仅受毛细管力作用条件下的三相接触角能够反映岩石的润湿性，认清水驱油过程中三相接触角的动态变化，对认识油水两相的渗流规律具有重要意义[13]。三相接触角的测量最初主要利用接触角测量仪测量开放空间内单个液滴的接触角[14⇓-16]，这主要以单个液滴为研究对象，主要注重从图像角度表征流固接触界面，再计算接触角，并没有考虑液滴在孔喉环境下的分布状态。微观CT扫描图像技术也被应用于孔喉系统内接触角的测量，Andrew等[17]、Klise等[18]针对微观CT三维扫描图像，建立了像素级的原位接触角自动测量方法，研究了不同岩石结构下接触角的分布规律。Scanziani等[19]通过分离三维岩心重构数据中的岩石、水和油三相，对岩石边缘进行直线拟合，并对油水界面进行圆弧拟合并取切线，通过测量两条直线的夹角完成三相接触角的自动测量。该方法实现了微观图像上较为精准的测量，但其反映的是静态条件下三相接触角的特征。Mirzaei等[20]依据边缘曲线，通过Harris角点检测函数确定接触点的精确位置，实现了液滴接触角的测量。AlRatrout等[21]利用“流体-流体”恒定的界面曲率寻找接触线和法向量，计算接触角，实现了不互溶流体原位接触角的测量。Ibekwe等[22]通过识别接触点、接触角粗估计、角度精确化等方法实现了图像中微观三相接触角的精确测量。可以看到上述成果对接触角的抽象化表征方法考虑较少，且难以对图像中的大量接触角进行统计分析，但其对三相接触角的图像处理分析策略具有借鉴意义。考虑到微流控实验图像具有阈值分明、曝光不均、相态较少的视觉特点，以及油、水和岩石三相接触界面测量的相关要求，本文借鉴Scanziani等[19]与AlRatrout等[21]采用CT重构数据测量微观三相接触角的思路，建立基于微流控实验数字图像的三相接触角自动测量方法。
1 微流控驱替实验及图像预处理
在微流控模型上进行水驱油实验，使用CCD（电荷耦合元件）相机拍摄实验过程并保存视频，视频与图像分辨率为（1 236×1 624）像素。将获取的饱和油图像用于提取岩石骨架，水驱达到残余油状态时的剩余油图像分别用于剩余油分类识别和接触角自动测量（见图1 ）。实验微流控模型润湿性为水湿，孔隙度20.63%，驱替至残余油状态时的含水饱和度为67.78%。实验用原油黏度为4 mPa·s，模型总流动区域约为2 cm×1.5 cm。
2 三相接触角自动测量方法
2.1 邻域检测法识别接触角顶点（接触点）
识别接触角关键点无需检测岩石相的全部像素点，检测图像中岩石内部的像素点会导致接触角检测效率降低，使用Canny边缘检测方法提取岩石边缘像素点[26]，可大幅提升岩石相的检测效率。在剩余油三值图像矩阵中，Canny边缘检测方法可清楚地区分油、水和岩石三相的灰度值，而需要寻找的三相接触角顶点（接触点）是油、水和岩石3种不同灰度值像素点的交点。
2.2 提取与拟合三相接触界面
检测得到图像中的全部有效接触点后，为了避免复杂孔喉环境对接触角测量的影响，需提取油、水和岩石三相之间的界面。基于三值图像矩阵，以接触点为起点，提取相邻岩石边缘点集和油水界面曲线点集，点集将用于接触角的测量，能被有效测量的接触角称为有效接触角。在提取过程中，为减小误差，将舍去接触点附近油相或水相像素点数过少的接触点。
在检测油水界面曲线的过程中，实时计算曲线检测窗口附近像素点集的平均斜率差，并通过调整邻域平均斜率差的值来量化岩石边缘或油水界面曲线的拐动变化，以确定何时停止检测界面曲线。
微观图像中的三相接触角是微尺度下的变量，实验拍摄时的图像分辨率会对测量结果产生误差。本文用三次样条插值法对原图像进行放大处理，三次样条插值法保留了图像中的3种阈值，可实现对油水界面、岩石边缘等曲线图像意义上的平滑拟合，从图像层面提升了三值化图像的计算精度。
三相界面经过提取和拟合后，还需将角的两边表示为可用于计算的两个端点，接触点、两个端点的抽象位置关系如图5 所示。图像坐标系中接触点坐标标记为（xc，yc），岩石边缘拟合直线端点坐标标记为（x0，y0），圆弧切线端点坐标标记为（qx，qy），这3个点将用于测量接触角。对于岩石边缘可直接在拟合直线上任取一点作为端点;对于油水界面拟合圆弧，则可依据接触点和拟合圆弧的圆心（标记为（x1，y1））的连线，作过接触点的垂线得拟合圆弧的切线并在其上任取一点作为端点。
2.3 三相接触角的自动测量
3 三相接触角的测量结果
3.1 接触角测量结果误差
人工测量的接触角结果往往存在主观性，无法体现测量方法对多组测量结果的统一约束，虽难以达到标准统一的要求，但其仍可作为对本方法自动测量结果的辅助验证。微流控实验在驱替至40 min 10 s时达到残余油状态，选取该时刻的残余油图像，在三值图像上人工测量了47个接触角数据，随后将接触角由小到大排序，与本文方法计算的接触角（同样排序）进行对比（见图7 ）。本文方法去除一部分无效接触角（数据点过少、图像噪点、图像边界等）后，共测量得到2 373个有效接触角。
假设人工测量值为真值，经过计算，本文方法测量结果的绝对误差平均值为5.40°，相对误差平均值为1.50%，可见本文方法测量结果与人工测量结果基本一致，说明三相接触角自动测量方法所得结果可靠。自动测量方法的误差主要来源于计算误差，与三相阈值分割、岩石边缘和油水界面提取、数据拟合算法等有关。
3.2 三相接触角分布特征
“油-水-岩石”三相接触角的大小与模型润湿性、毛管压力、驱替参数等相关。而传统实验室测量的“空气-油-岩石”三相接触角为静态结果，与油水驱替条件下的微观三相动态接触角存在明显差异，微流控模型油水两相流动中，三相接触角体现了界面的复杂受力变化情况，随着驱替时间的延长，同一位置的三相接触角也在不断地变化。
为了更加直观地观测接触角测量结果的二维空间分布，采用不同颜色实心圆标记了4个驱替时刻原三值图像中的所有有效接触角（见图8 ）。分析发现驱替至16 min 5 s时，注入端附近接触角以锐角为主（见图8 a）;驱替至21 min 10 s时采出端见水，采出端附近存在大量钝角接触角（见图8 b右侧）;驱替至30 min 14 s时，接触角数量逐渐增多，钝角数量也大幅提升（见图8 c）;驱替至40 min 10 s时驱替过程达到水驱残余油状态，接触角在平面上分布趋于均匀，钝角数量略有提升（见图8 d）。
为了实现对接触角测量结果数据的概率密度分析，对这4组接触角数据分别作正态拟合概率密度曲线和高斯核密度估计概率曲线。其中，图9 a—图9 d为正态分布拟合后的结果，可分析接触角数据分布的均匀、集中程度，可以看到，随着驱替的进行，接触角均值逐渐增大，接触角峰值由锐角向钝角偏移，同时数据更加集中，表明模型整体逐步由水湿向油湿转变。图9 e—图9 h为接触角数据的高斯核密度估计概率曲线，主要体现接触角数据在不同角度范围内的分布特点，曲线上每个向上凸起的峰代表接触角数据在某个角度范围内的聚集。可以看到水驱前期（16 min 5 s），核密度估计曲线包含一个锐角峰和一个钝角峰，锐角峰的峰值更高，表明数据更加集中于锐角，接触角均值为锐角，推测此时微流控模型整体呈现水湿特征;随着驱替的进行，曲线高峰逐渐右移，驱替至21 min 10 s时接触角均值已变为钝角，核密度估计曲线仍然包含一个锐角峰和一个钝角峰，但钝角峰更高，数据更加集中于钝角;驱替至30 min 14 s时接触角数量开始增多，对应图8 c中接触角二维空间分布数量增多的情形，此时接触角均值与21 min 10 s时较为接近，核密度估计曲线在钝角区域包含2个峰;驱替至40 min 10 s时均值继续增大，数据比水驱前期更为集中，核密度估计曲线在钝角区域的2个峰与30 min 14 s时相比向右偏移，数据集中于两处更大的钝角。
3.3 不同驱替阶段接触角动态变化特征
微流控模型的整体润湿性变化是复杂的动态过程，且较难从微观角度研究单个“油-水-岩石”三相接触角处体现的润湿性，研究不同驱替阶段接触角数据的动态统计特征，更能有效表征微流控模型的整体润湿性。
进一步观察整个驱替过程中三相接触角和油水接触关系的动态变化，需要对视频取帧分析，计算各时刻的图像含水饱和度，测量各含水阶段的三相接触角。微流控驱替实验视频为25帧/s，总时长40 min 13 s，水驱约40 min后，模型中基本呈现为残余油状态，计算得该时刻含水饱和度为67.78%。视频分析过程仅测量水相一侧的接触角，测量数据包括水驱开始至达到残余油状态时的所有接触角数据，基本可以体现微观水驱油实验中油水接触关系的动态变化过程。
水驱油过程中微观三相接触角的变化，反映了流动条件下毛细管力与黏滞力间的相互作用。驱替结束后，黏滞力不再影响流体的界面分布，其接触角测量结果可用于模型的润湿性分析。图10 为不同含水阶段接触角均值变化曲线，主要体现三相接触角均值由锐角到钝角的动态变化过程。当含水饱和度小于12.02%时，所测三相接触角的均值几乎全部为锐角;当含水饱和度大于等于12.02%时，所测三相接触角的均值均为钝角，且基本分布在100°~120°;达到残余油状态时（含水饱和度约67.78%），接触角均值约为113.23°。
水驱油过程中，油水界面始终受黏滞力与毛细管力的相互作用，接触角变化反映的是毛细管力与黏滞力的相互作用结果。驱替初期，黏滞力占主导，接触角普遍为锐角;油水分布逐渐向残余油状态转化过程中，黏滞力对剩余油的作用逐步减弱，接触角也逐渐转化为钝角，这反映以黏滞力为主导逐步向毛细管力为主导过渡的过程;驱替过程结束，残余油仅受毛细管力作用，接触角以钝角为主，实验模型呈偏油湿特征，此时微流控模型的润湿性发生反转，分析认为这与原油中酸性物质使二氧化硅表面性质发生变化有关。
3.4 不同流态剩余油接触角变化特征
4 结论
对于水湿岩石，水驱油过程中，驱替前期接触角均值以锐角为主，随驱替时间的延长，油相受力逐渐由黏滞力为主导向毛细管力为主导过渡，接触角均值不断增大，至驱替中后期以钝角为主。
多孔流和柱状流剩余油在整个驱替过程中接触角变化幅度较小且均值皆为钝角，该类剩余油形成后流动性差，动用难度较大;簇状流在驱替过程中三相接触角均值整体呈上升趋势，簇状流剩余油动用难度逐渐升高;膜状流整体呈钝角，水驱后期接触角均值最大，动用难度最大。驱替结束后，不同流态剩余油油相仅受毛细管力作用的影响，接触角均为钝角，在原油作用下模型表现出偏油湿的特征。
符号注释：
n——需要计算的接触角的个数;（xc，yc）——接触点坐标;（x0，y0）——岩石边缘拟合线段端点坐标;（qx，qy）——圆弧切线端点坐标;（x1，y1）——圆心坐标。