微米毛细管内气液两相动态驱替图像的接触角获取方法

摘要

本发明提供了微米毛细管内气液两相动态驱替图像的接触角获取方法，通过显微镜获取当前微米毛细管内气液两相动态驱替过程的图像再转化成二值图像，利用连通区域方法通过像素点的颜色提取出构成管道形边界图像的边界坐标点；初步确定相界面和两个壁面的坐标点分布趋势，进而初步划分出分别构成相界面和两个壁面的坐标点；对构成相界面的坐标点进行三次多项式拟合后得到三次函数，对两个壁面的坐标点进行一次函数拟和合得到两条对应的表达式，求解后即可得到相界面两个端点的精确坐标，即得到相界面两端点处的接触角。本发明的优点在于鲁棒性强,测量结果稳定，具有较强容错性，对于质量不好的图像也能提供准确的测量结果。

说明书

技术领域

本发明涉及地质领域，特别涉及一种微米毛细管内气液两相动态驱替图像的接触角获取方法。

背景技术

单个微米尺度的毛细管中的两相驱替机理是自然环境多孔介质内的两相驱替的基础，多孔介质中的两相驱替发生在众多的工业或自然过程中，比如致密油藏的水驱气、溶解气驱、水驱油等过程，上述的工业过程的特点是粘性力微弱(毛管数的范围10^-10<Ca<10^-5)，重力微弱(邦德数B_o小于10^-4)。影响驱替过程的关键因素为毛管力，而计算毛管力的关键参数是两相界面与管道壁面的接触角。毛管力的计算按照Young-Laplace公式计算，如下式所示，

γ是界面张力，接触角θ的实际范围为0～180°，当θ<90°时P_c为正值，当θ>90°时，P_c为负值，因此毛管力有可能是阻碍驱替的阻力也可能是推进驱替的动力。不仅如此，接触角评估的错误将会使毛管力计算存在严重的误差，从而导致驱替过程的评估失误。因此，获取上述工业条件下的动态接触角信息对于评估整个多孔介质内两相驱替过程至关重要。

现已较完善的测量相界面运动状态下(动态)接触角的技术，比如平板浸没法、毛细管液升法。平板浸没法是将一块待测固体平板，以一定速度、以平板垂直于液面的姿态进入或抽离液面，这个过程采用高速摄像机拍摄液面与平板接触点处相界面的接触角，在通过图像测量接触角的大小；或者采用测力计测量平板在这个运动过程中的受力，测算出界面张力反推出接触角。

毛细管液升法，是将一根已知半径的毛细管入口垂直伸入湿润相液体的液面下，通过高速摄像机拍摄湿润相液面上升的高度及对应的时间，通过Washburn公式推算出接触角。

但是采用上述两种方案都不能够测量毛管数10^-10<Ca<10^-5及邦德数Bo<10^-4工况下的动态接触角。

发明内容

本发明的目的是要提供一种微米毛细管内气液两相动态驱替图像的接触角获取方法。

特别地，本发明提供微米毛细管内气液两相动态驱替图像的接触角获取方法，包括如下步骤：

步骤100，利用显微镜获取当前微米毛细管内气液两相动态驱替过程的图像，并传送至控制系统，控制系统对该图像进行灰度转化和滤波处理，然后利用阈值分割算法将灰度图像转化成二值图像，利用连通区域方法通过像素点的颜色提取出构成管道形边界图像的边界点；

步骤200，识别边界点中代表管道形边界图像中表示左右两个相对界面的相界面和末端面及上下两个相对壁面的坐标点，并根据坐标值的大小对各坐标点进行排序和编号，依据构成未端面的坐标点特性，将所有未端面的坐标点删除，同时确定两个壁面与未端面相交的两个末端接触点坐标；

步骤300，对剩余坐标序列利用能量函数根据坐标值的变化对各坐标点进行排序，初步确定相界面和两个壁面的坐标点分布趋势，然后设定能量阈值选取变化低于该能量阈值的坐标点并进行划分和编号，通过取编号的平均值即可确定相界面两个端点的粗略坐标，进而初步划分出分别构成相界面和两个壁面的坐标点；

步骤400，对构成相界面的坐标点进行三次多项式拟合后得到三次函数，对两个壁面的坐标点进行一次函数拟合得到两条对应的表达式，求解后即可得到相界面两个端点的精确坐标，然后将其分别代入三次函数的一次导数中，即得到相界面两端点处的接触角。

在本发明的一个实施方式中，所述步骤100中提取出构成管道形边界图像的边界点的具体处理步骤如下：

步骤110，将二值图像中黑色区域的像素点记为数值0，白色区域的像素点记为数值1；

步骤120，将所有像素值I(x_i,y_i)＝1的坐标点同其相邻的像素点进行“与”运算；若运算结果为1则将像素点标记为内部点，若运算结果为0则将该像素点标记为边界点；

步骤130，运算结束后将所有边界点提出来，即得到所有构成管道形边界图像的边界点。

在本发明的一个实施方式中，构成末端面的坐标点特性如下：

由于未端面为一条竖线，因此其横坐标值变化一致或在一定范围内，而纵坐标变化为渐增，根据所有坐标点中横、纵坐标值的大小对各坐标点进行从小到大的排序，选择横坐标值变化在预定范围内的坐标点全部删除，即删除了未端面。

在本发明的一个实施方式中，所述初步确定相界面和两个壁面的坐标点分布趋势的过程如下：

将坐标点序列上的任意一个坐标点(x_i,y_i)，与其在序列上相距m位的坐标点(x_i+m,y_i+m)按照下式进行比较：

E_i＝||x_i-x_i+m|-|y_i-y_i+m||

可得到相界面上的坐标点在x方向上的变化较弱，在y方向上变化较强，因此E_i≈m；而两个壁面上的坐标在x方向上的变化较强，在y方向上变化较弱，会有E_i≈m；而从两个壁面过渡到相界面的这部分边界点其在x方向和y方向上的变化基本相同E_i≈0。

在本发明的一个实施方式中，初步划分出分别构成相界面和两个壁面的坐标点的过程如下：

步骤310，从任意一个坐标点开始，将另一个坐标点以或者的条件相对坐标点进行排列，待所有的坐标点排列成一个序列后，设定一个能量阈值E_T，其取值范围为1≤E_T≤m，在所有坐标点中搜索满足E_i≤E_T的所有坐标点，得到由两段连续编号坐标点组成的区域序列；

步骤320，将区域序列中编号最小者作为上壁面的一个端点，其与对应的上壁面的未端接触点坐标之间的坐标点为形成上壁面的坐标点编号范围；

将区域序列编号中最大者作为下壁面的一个端点，其与对应的下壁面的未端接触点坐标之间的坐标点为形成下壁面的坐标点编号范围；

对区域序列中所有编号求和后取平均值，则第一个小于平均值的编号为相界面的一个端点，而第一个大于平均值的编号则为相界面的另一个端点，两个端点之间的坐标点则为形成相界面的坐标点编号范围。

在本发明的一个实施方式中，还包括对坐标点排序时的顺序修正程序：

设置边界点序列的能量函数，由两个壁面的末端接触点坐标之间的区域且根据两个未端接触点坐标的纵坐标的差值确定极大值，即E_max＝|y_D-y_C|，再通过寻找单个最大值算法确定E_i＝E_max的坐标编号，设其编号为i＝imax，将坐标点序列中的第一坐标点分别与能量函数进行循环位移，将编号imax变成编号0，仍然是n+1个坐标点，从0开始编号，新的编号j和原编号i的关系为：

通过对排序的坐标点进行循环位移方式调整排序后，使得随着编号的递增，代表管道形边界图像的四个边界点依次出现。

在本发明的一个实施方式中，所述步骤400中，得到相界面两个端点的精确坐标的步骤如下：

步骤410，设构成相界面的第一个坐标点为P，最后一个坐标点为Q，用P至Q范围内所有坐标点进行三次多项式拟合，以横坐标x为因变量，纵坐标y为自变量，拟合得到弯月面形状的三次函数x＝f_PQ(y)；

步骤420，设构成上壁面的第一个坐标点为D，最后一个坐标点为O；构成下壁面的第一个坐标点为R，最后一个坐标点为C；采用一次函数拟合，以横坐标y为因变量，纵坐标x为自变量，得到两条直线的表达式y＝k_DOx+c_DO及y＝k_RCx+c_RC，其中k_DO和k_RC为直线的斜率；

步骤430，根据拟合得到的三次函数及两条直线表达式，求解以下两个方程组；

即得到相界面上端点A的坐标(x_A,y_A)以及相界相下端点B的坐标(x_B,y_B)。

在本发明的一个实施方式中，获得接触角的过程如下：

根据相界面的拟合函数x＝f_PQ(y)，将上端点A的坐标(x_A,y_A)以及下端点B的坐标(x_B,y_B)分别代入f_PQ的一次导数中即可计算出接触角θ_A和θ_B：

在本发明的一个实施方式中，还包括倾斜修正程序：

设两个壁面的倾斜角度为上壁面和下壁面倾斜角度的平均值，即(arctank_DO+arctank_RC)/2；则修正后的接触角θ_A和θ_B为：

本发明的优点在于鲁棒性强,测量结果稳定，具有较强容错性，对于质量不好的图像也能提供准确的测量结果；采用该计算方法编程简单，程序运行效率高，采用常规处理器处理一张2300×1730像素的图片耗时小于1s。

附图说明

图1是本发明一个实施方式的接触角获取方法流程示意图；

图2是判断任意一个坐标点是否是边界点的示意图；

图3是识别出的形边界图像示意图；

图4是相界面、上下壁面和末端面的示意图；

图5是识别并提出的形边界图像示意图；

图6是形边界序列上的坐标点排列及各个顶点的示意图；

图7是能量函数算法的示意图；

图8循环位移前的坐标点序列示意图；

图9相界面和两个壁面的拟合以及接触点的确定示意图。

具体实施方式

如图1所示，本发明一个实施例的微米毛细管内气液两相动态驱替图像的接触角获取方法，包括如下步骤：

步骤100，利用显微镜获取当前微米毛细管内气液两相动态驱替过程的图像，并传送至控制系统，控制系统对该图像进行灰度转化和滤波处理，然后利用阈值分割算法将灰度图像转化成二值图像，利用连通区域方法通过像素点的颜色提取出构成管道形边界图像的边界坐标点；

图像中的白对应灰度值255，黑对应0，其余不同程度的灰色的灰度值在0和255之间。滤波处理是为了去除图像中像素的噪点、小面积的独立斑块。

其中提取出构成管道形边界图像的边界点的具体处理步骤如下：

步骤110，将二值图像中黑色区域的像素点记为数值0，白色区域的像素点记为数值1；

这里采用阈值分割算法将灰度图像转化成二值(黑白)图像，此时液体区域和固体区域被划分为黑色区域(数值为0)，空气所占据的区域为白色区域(数值为1)，因此图像中存在两个区域黑色区域及白色区域。

步骤120，将所有像素值I(x_i,y_i)＝1的坐标点同其相邻的像素点进行“与”运算；若运算结果为1则将像素点标记为内部点，若运算结果为0则将该像素点标记为边界点；

如图2所示，这里是通过识别连通区域来提取边界，先将所有像素值I(x_i,y_i)＝1的坐标点同其相邻的8个(若在图像边缘则为5个，若在图像四个角落中则为3个)像素点的像素值为1进行“与”运算；若运算结果为1则将像素点标记为内部点，若像素点相邻的8个坐标点中存在一个像素值I＝0(黑色点)的像素点，则运算的结果为0，该坐标点标记为边界点。

步骤130，运算结束后将所有边界点提出来，即得到所有构成管道形边界图像的边界点。如图3、4所示。

步骤200，识别边界坐标点中代表管道形边界图像中表示左右两个相对界面的相界面和末端面及上下两个相对壁面的坐标点，并根据坐标值的大小对各坐标点进行排序和编号，依据构成未端面的坐标点特性，将所有未端面的坐标点删除，同时确定两个壁面与未端面相交的两个末端接触点坐标；

对提取后各坐标点进行排序的方式如下：从任意一个边界点开始，将另一个边界点排列在之前或之后，排列条件为或者以此类推，将k个边界点排列成一个序列即获得了构成形边界的所有点的坐标，如图5所示，以方便后续对相界面、接触点和壁面等的识别。

这里的编号是指从0开始按1、2、3……的顺序对排序后的坐标点依次进行编号。

构成末端面的坐标点特性如下：

如图6所示，由于未端面为一条竖线，理论上应该是一条垂直的直线，但实际获取的图像中其可能会产生一定的偏差，因此其横坐标值变化一致或在一定范围内，而纵坐标变化为渐增，根据前面获取的坐标点中横、纵坐标值的大小对各坐标点进行从小到大的排序，选择横坐标值变化在预定范围内的坐标点全部删除，即删除了未端面。

本发明中的壁面指显微图像中微毛细管道的内壁，相界面指气体和液体的交界面，接触点指相界面与上下两个壁面的交点，未端面则为与相界面相对的管道末端。具体地，设管道末端的横坐标为x＝W_s，对于任意边界点i，只要满足x_i＝W_s，则将点i从边界点的序列中剔除，则如图所示，线段c(末端面)会被剔除，而剩余的点中与上壁面(线段d)和下壁面(线段b)两端相邻的末端接触点坐标D和C则可以获得。

步骤300，对剩余坐标序列利用能量函数根据坐标值的变化对各坐标点进行排序，初步确定相界面和两个壁面的坐标点分布趋势，然后设定能量阈值选取变化低于该能量阈值的坐标点并进行划分和编号，通过取编号的平均值即可确定相界面两个端点的粗略坐标，进而初步划分出分别构成相界面和两个壁面的坐标点；

初步确定相界面和两个壁面的坐标点分布趋势的过程如下：

设边界点中除去线段c之后还剩余n+1个像素点，坐标组成的坐标点序列为(x₀,y₀)、(x₁,y₁)…(x_n,y_n)，在实际图像中，相界面并不是一条光滑的曲线，上壁面和下壁面同样也不是笔直的直线。通过能量函数来计算坐标点的变化趋势，可将相界面和上壁面、下壁面进行粗略区分。

将坐标点序列上的任意一个坐标点(x_i,y_i)，与其在序列上相距m位的坐标点(x_i+m,y_i+m)按照下式进行比较：

E_i＝||x_i-x_i+m|-|y_i-y_i+m||

可得到相界面上的坐标点在x方向上的变化较弱，在y方向上变化较强，因此E_i≈m；而两个壁面上的坐标在x方向上的变化较强，在y方向上变化较弱，会有E_i≈m；而从两个壁面过渡到相界面的这部分边界点其在x方向和y方向上的变化基本相同E_i≈0。

初步划分出分别构成相界面和两个壁面的坐标点的过程如下：

步骤310，从任意一个坐标点开始，将另一个坐标点以或者的条件相对坐标点进行排列，待所有的坐标点排列成一个序列后，设定一个能量阈值E_T，其取值范围为1≤E_T≤m，在所有坐标点中搜索满足E_i≤E_T的所有坐标点，得到由两段连续编号坐标点组成的区域序列；

比如编号10、11、12、20、21、22中，10、11和12是连续的，20、21和22是连续的，但是12到20是不连续的。

步骤320，将区域序列中编号最小者作为上壁面的一个端点，其与对应的上壁面的未端接触点坐标之间的坐标点为形成上壁面的坐标点编号范围；

将区域序列编号中最大者作为下壁面的一个端点，其与对应的下壁面的未端接触点坐标之间的坐标点为形成下壁面的坐标点编号范围；

对区域序列中所有编号求和后取平均值，则第一个小于平均值的编号为相界面的一个端点，而第一个大于平均值的编号则为相界面的另一个端点，两个端点之间的坐标点则为形成相界面的坐标点编号范围。

如图7所示，假设这两段连续的区域的首末两个坐标点的编号为总共u个坐标点(如编号分别是10、11、12、20、21、22，此时u＝6)；搜索序列中最小值(例子中的编号10)，令此点为O点；搜索序列中最大值(例子中的编号22)，令此点为R点；将序列所有编号进行求和运算并处以u(在例子中为(10+11+12+20+21+22)/6＝16)，令计算结果为令序列中第一个小于的编号为P(例子中的编号12)，令第一个大于的编号为Q(例子中的编号20)。其中编号P到Q范围内的坐标点是粗略的相界面，而编号D到O范围内的坐标点是粗略的下壁面b，编号R到C范围内的坐标点是粗略的上壁面d。

前述得到的坐标点序列不一定按照上壁面右端点(D)-相界面上端点(A)-相界面下端点(B)-下壁面右端点(C)的顺序排列，为了避免保证能按前述方法正确获得点O、P、Q和R，需要在采用循环位移的方式，调整坐标点点序列的排序，使得随着编号的递增，关键点D-A-B-C依次出现。具体的对坐标点排序时的顺序修正程序如下：

如图8所示，设置边界点序列的能量函数，由两个壁面的末端接触点坐标之间的区域且根据两个未端接触点坐标的纵坐标的差值确定极大值，即E_max＝|y_D-y_C|，再通过寻找单个最大值算法确定E_i＝E_max的坐标编号，设其编号为i＝imax，将坐标点序列中的第一坐标点分别与能量函数进行循环位移，将编号imax变成编号0，仍然是n+1个坐标点，从0开始编号，新的编号j和原编号i的关系为：

通过对排序的坐标点进行循环位移方式调整排序后，使得随着编号的递增，代表管道形边界图像的四个边界点依次出现。

步骤400，对构成相界面的坐标点进行三次多项式拟合后得到三次函数，对两个壁面的坐标点进行一次函数拟和合得到两条对应的表达式，求解后即可得到相界面两个端点的精确坐标，然后将其分别代入三次函数的一次导数中，即得到相界面两端点处的接触角。

如图9所示，得到相界面两个端点的精确坐标的步骤如下：

步骤410，设构成相界面的第一个坐标点为P，最后一个坐标点为Q，用P至Q范围内所有坐标点进行三次多项式拟合，以横坐标x为因变量，纵坐标y为自变量，拟合得到弯月面形状的三次函数x＝f_PQ(y)；

步骤420，设构成上壁面的第一个坐标点为D，最后一个坐标点为O；构成下壁面的第一个坐标点为R，最后一个坐标点为C；采用一次函数拟合，以横坐标y为因变量，纵坐标x为自变量，得到两条直线的表达式y＝k_DOx+c_DO及y＝k_RCx+c_RC，其中k_DO和k_RC为直线的斜率；

步骤430，根据拟合得到的三次函数及两条直线表达式，求解以下两个方程组；

即得到相界面上端点A的坐标(x_A,y_A)以及相界相下端点B的坐标(x_B,y_B)。

获得接触角的过程如下：

根据相界面的拟合函数x＝f_PQ(y)，将上端点A的坐标(x_A,y_A)以及下端点B的坐标(x_B,y_B)分别代入f_PQ的一次导数中即可计算出接触角θ_A和θ_B：

其中，将上端点A的坐标(x_A,y_A)带入f_PQ的一次导数中可以得到相界面的函数在上端点A的切线斜率，通过arctan函数将斜率转化为角度，由于测量的是位于液体那一侧的接触角，因此加上π/2就可以得到θ_A；将下端点B的坐标(x_B,y_B)分别代入f_PQ的一次导数中，可以求得相界面在下端点B的切线斜率，同样采用arctan将斜率转化成角度，由于下端点B处的斜率是负值，因此转化成角度后需要乘上一个负号，再加上π/2就可以得到θ_B。

还包括对两条壁面的倾斜修正程序：

设两个壁面的倾斜角度为上壁面和下壁面倾斜角度的平均值，即(arctank_DO+arctank_RC)/2；则修正后的接触角θ_A和θ_B为：

本发明的优点在于鲁棒性强,测量结果稳定，具有较强容错性，对于质量不好的图像也能提供准确的测量结果；采用该计算方法编程简单，程序运行效率高，采用常规处理器处理一张2300×1730像素的图片耗时小于1s。

至此，本领域技术人员应认识到，虽然本文已详尽示出和描述了本发明的多个示例性实施例，但是，在不脱离本发明精神和范围的情况下，仍可根据本发明公开的内容直接确定或推导出符合本发明原理的许多其他变型或修改。因此，本发明的范围应被理解和认定为覆盖了所有这些其他变型或修改。