一种微米毛细管中液液两相驱替图像的接触角获取方法

技术领域

本发明涉及一种微米毛细管中液液两相驱替图像的接触角获取方法，适用于微尺度两相流动研究领域、石油强化开采领域的实验。

背景技术

单个微米尺度的毛细管中的两种不相融液体的驱替机理是自然环境多孔介质内的两相驱替的基础，多孔介质中的两相驱替发生在众多的工业或自然过程中，比如致密油藏的水驱油、聚合物驱油等过程，上述的工业过程的特点是粘性力微弱(毛管数的范围10^-10<Ca<10^-5)，重力微弱(邦德数Bo小于10^-4)。这两个无量纲数的定义为Ca＝μV/γ，μ为液体粘度，V是相界面速度，γ是界面张力；Bo＝ρgh²/γ，ρ是液体密度，g为重力加速度，h是流动区域的特征尺寸。影响驱替过程的关键因素为毛管力，而计算毛管力的关键参数是两相界面与管道壁面的接触角。毛管力的计算按照Young-Laplace公式计算，如式(1)所示。接触角θ的实际范围为0～180°，当θ<90°时P_c为正值，当θ>90°时，P_c为负值，因此毛管力有可能是阻碍驱替的阻力也可能是推进驱替的动力。不仅如此，接触角评估的错误将会使毛管力计算存在严重的误差，从而导致驱替过程的评估失误。因此，获取上述工业条件下的动态接触角信息对于评估整个多孔介质内液液驱替过程至关重要。

在微尺度毛细管中进行两相不相融液体的驱替实验是研究动态接触角的常用方法。已有的实验往往对相界面的显微图像进行人工分析，手动测量接触角。而已有的商业设备采用的图像处理方法也仅限与分析在固体表面的液滴形状(座滴法测量接触角)，液滴形状的图像分析算法包括对称形状分析方法(ADSA-P)、高宽比法、椭圆拟合法等。但是以上算法是针对平面上液滴的形状检测提出的，适用于对比度高的图像，不适用于微管道内液液驱替时低对比度的相界面形状。要测量出微管道内液液驱替相界面的接触角，需要从图像中识别出相界面。

发明内容

本发明的目的是提供一种微米毛细管中液液两相驱替图像的接触角获取方法，该方法能够从图像中识别出壁面、相界面和接触点，同时测量出接触角，接触点位置。

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

一种微米毛细管中液液两相驱替图像的接触角获取方法，将获得的显微图像进行坐标划分，设图像的像素尺寸宽为L_x和长L_y；然后进行以下步骤：

步骤一，首先通过多极值搜索技术沿着y方向搜索两个壁面最左端的两个点A和B，得到坐标分别为(0,y_A)和(0,y_B)；

步骤二，沿着A-B连线中点的E(0,y_E)的水平方向y＝x_E搜索单个局部最小值，获得相界面上某点F，坐标点为(x_F,y_F)；

步骤三，沿着y方向采用排序算法搜索单极值，进而搜索部分壁面，为了避免相界面对壁面的搜索，设置一个范围以避开相界面，令0<x_D1＝x_C1<x_F，在0≤x≤x_D1和0≤y≤y_E范围内搜索得到壁面A-D1，其中D1的坐标为(x_D1,y_D1)；在0<x_f<x_C1和y_E≤y≤L_y范围内搜索得到壁面B-C1，其中C1的坐标为(x_C1,y_C1)；

步骤四，沿着x方向采用排序算法搜索单极值，进而搜索部分相界面，为了避免两个壁面对相界面的搜索，设置一个范围以避开两个壁面，令0<y_D2<y_F<y_C2<L_y，在0≤x≤L_x和y_D2≤y≤y_C2范围内搜索得到部分相界面C2-D2，两个点的坐标分别是(x_D2,y_D2)、(x_C2,y_C2)；

步骤五，采用一次函数分别拟合A-D1和B-C1线段内所有的坐标点得到一次函数f_AD和f_BC，采用三次函数拟合C2-D2曲线内的所有点，得到三次函数f_CD；求解f_AD与f_CD的交点得到接触点D的坐标(x_D,y_D)，求解f_AD与f_CD的交点，得到接触点C的坐标(x_C,y_C)；线段B-C和A-D为壁面，曲线C-D为相界面；

步骤六，根据步骤五得到的相界面的拟合函数x＝f_CD(y)，将C的坐标(x_C,y_C)以及D的坐标(x_D,y_D)分别代入f_CD的一次导数中，求得相界面在C点和D点的接触角θ_C和θ_D；f_CD的一次导数f′_CD为

f′_CD＝3a_CDy²+2b_CDy+c_CD(13)

将纵坐标y_C和y_D代入式(13)中，得到斜率，将斜率转换成角度得到接触角θ_C和θ_D；

本发明进一步的改进在于，步骤一的具体过程如下：

点A和点B是x＝0的直线上灰度极小值点，也是灰度梯度沿着y方向存在极大值的两个点，因此找到灰度梯度极大值的两个点的坐标并加以区分，确定点A和点B的坐标；

令任意坐标点(x,y)的灰度值为I(x,y)，则进行以下步骤：

(1)对x＝0的直线上的坐标点灰度值进行向前差分来表示变化梯度dI，dI存在两个极大值点；

(2)设x＝0的直线上所有点的y坐标按照灰度梯度值I(0,y)进行降序排列，得到一个坐标序列(y₁,y₂,...y_n)，坐标y₁是其中一个极大值点的位置，按顺序搜索坐标序列，搜索到第一个满足式(3)的坐标点y_k，其中1<k≤n；式(3)中α为小于1的系数，h为点A到点B的竖直距离h＝y_B-y_A；

|y_k-y₁|>αh(3)

(3)比较y₁和y_k的大小，令二者中较小值为y_min，较大值为y_max，则点A和点B的纵坐标确定为y_A＝y_min，y_B＝y_max。

本发明进一步的改进在于，步骤二的具体过程如下：

(1)求取A-B连线中点E的纵坐标y_E；

(2)将直线y＝y_E上所有坐标点的灰度值I(x,y_E)按式(5)求灰度的梯度，其中x∈[0,L_X]；得到直线y＝y_E上所有坐标点的灰度梯度dI(x,y_E)；

(3)由于dI(x,y_E)只有一个极大值，搜索dI(x,y_E)序列中的极大值点，确定相界面上点F的横坐标x_F，而点F的纵坐标x_F＝x_E。

本发明进一步的改进在于，步骤三的具体过程如下：

(1)为了沿着y方向搜索壁面时不受到相界面的影响，限定搜索的x范围为0≤x≤x_D1或0≤x≤x_C1，其中x_D1＝x_C1＝βx_F，β是一个小于1的系数；

(2)为了搜索壁面A-D1，横坐标x_i遍历0至x_D1，对于每个横坐标x_i，求取直线x＝x_i上灰度I(x_i,y)在y方向的梯度，如式(6)所示；搜索dI(x_i,y)上的极大值，获得对应壁面坐标点的纵坐标y_i，一系列坐标点(x_i,y_i)组成了壁面A-D1；

(3)为了搜索壁面B-C1，横坐标x_j遍历0至x_C1，对于每个横坐标x_j，求取直线x＝x_j上灰度I(x_j,y)在y方向的梯度，将式(6)中x_i替换成x_j计算dI(x_j,y)，搜索dI(x_j,y)上的极大值，获得对应壁面坐标点的纵坐标y_j，一系列坐标点(x_j,y_j)组成了壁面B-C1。

本发明进一步的改进在于，步骤四的具体过程如下：

(1)沿x方向搜索灰度的极小值确定相界面，需要避开两个壁面，设定y方向的搜索范围为y_D2<y<y_C2，其中y_D2和y_C2分别为

式(7)中，h为点A到点B的距离，即h＝y_B-y_A，ε是一个系数，其取值范围为0<ε<0.5；

(2)为了搜索相界面，y_i坐标遍历y_D2至y_C2，对于每个纵坐标y_i，求取直线y＝y_i上灰度I(x,y_i)在x方向的梯度dI(x,y_i)，如式(8)所示；搜索dI(x,y_i)上的极大值，获得对应壁面坐标点的纵坐标y_i，一系列坐标点(x_i,y_i)组成了壁面C2-D2；

本发明进一步的改进在于，步骤五的具体过程如下：

(1)采用一次函数拟合A-D1线段内所有点，以x坐标为自变量，y为因变量，得到的一次函数f_AD如式(9)所示，a_AD和b_AD为一次函数的系数，

y＝f_AD(x)＝a_ADx+b_AD(9)

采用一次函数拟合B-C1线段内所有点，以x坐标为自变量，y为因变量，得到的一次函数f_BC如式(10)所示，a_BC和b_BC为一次函数的系数

y＝f_BC(x)＝a_BCx+b_BC(10)

采用三次函数拟合相界面C2-D2，以y坐标为自变量，x为因变量，得到的三次函数f_CD，如式(11)所示；

x＝f_BC(y)＝a_CDy³+b_CDy²+c_CDy+d_CD(11)

(2)接触点C和D分别是两个壁面与相界面的交点，因此通过求解f_CD与f_AD的交点，及f_CD与f_BC交点，获得接触点C和D的坐标(x_C,y_C)和(x_D,y_D)，上述过程是求解式(12)所示的方程组；

(3)相界面、壁面上所有坐标点的确定：要求取壁面A-D上的所有坐标点，将每个满足x_A≤x≤x_D的横坐标x带入式(9)中进行计算得到对应y坐标，从而获得壁面A-D上所有点的坐标；要求取壁面B-C上的所有坐标点，将每个满足x_B≤x≤x_C的横坐标x带入式(10)中进行计算得到对应y坐标，从而获得壁面B-C上所有点的坐标；对于相界面C-D上的所有坐标点，将每个满足y_C≤y≤y_D的纵坐标y带入式(11)中进行计算得到对应x坐标，从而获得相界面C-D上所有点的坐标。

本发明进一步的改进在于，还包括步骤七：当实际获取得到的壁面不是水平时，对式(14)获得的接触角进行修正，具体过程如下：壁面的倾斜角度为壁面A-D1和B-C1的倾斜角度平均值，即(arctan>AD+arctan>BC)/2，则接触角的最终值为：

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

本发明采用三次多项式拟合相界面并计算接触角，克服了传统算法ADSA-P只能拟合对称相界面的局限，该发明可以适应不对称相界面，更适合用于实际实验中获得的测量图像；该发明基于图像局部灰度极值的方式，实现了微管道中液液驱替的低对比图像的相界面，进而得到接触角。本发明基于图像局部灰度极值的方式，实现了微管道中液液驱替的低对比图像的相界面、壁面和接触点的识别，算法容错率高，可以适应存在噪声的图像，可以正确、高效识别出相界面、壁面和接触点，以上信息可用于正确测量接触角。该发明实现了在低毛管数、低邦德数状态下的液液驱替动态接触角测量。采用三次多项式对相界面拟合并测量接触角，该方法可适应不对称相界面的拟合，算法鲁棒性强，极大提高了数据处理的效率及测量精度。

附图说明

图1为液液驱替的显微图像。

图2为液液驱替的灰度图像中的灰度分布图。

图3为液液驱替图像的坐标设定。

图4为图像处理中关键点的位置。

图5为灰度及灰度梯度的分布在x＝0的直线上存在两个极小值点。

图6为灰度及灰度梯度的分布在x＝0的直线上存在两个极大值点。

图7为两个接触角为两个接触点上相界面与壁面形成的夹角。

图8为三次多项式导数与接触角之间的关系。

具体实施方式

下面结合附图对本发明进行详细说明。

本发明提供了在图像中寻找局部极值的方法，确定边界点，再从边界点识别壁面、相界面以及接触点。

液液驱替图像特征具体如下：

液液驱替过程的显微图像采用的是入射光(显微镜的底光源照射)观察方式，因此可以获得壁面、相界面与体相液体明显区分的图像。液液驱替的显微图像属于低对比度图像，如图1和图2所示。对比度低的图像的直方图分布特点就是灰度分布集中在一定区域内，高灰度和低灰度区域的像素数量较少。图2中两条直线为微孔道的壁面，曲线为相界面。利用这一特点，通过计算，从一系列显微图像中识别出相界面、测量接触角、相界面位移等其他参数。本方法可高效、精确地处理数量庞大的图像数据，极大提高了微毛细管中动态接触角研究的效率。

对比度低的图像的直方图分布特点就是灰度分布集中在一定区域内。液液驱替图像中，两条直线为微孔道的壁面，而曲线为相界面。

微孔道壁面和相界面的识别

此处边界的含义为显微图像中的微孔道壁面以及不相融液相之间的相界面。将获得的图像按照图3所示方法进行坐标划分。设图像的像素尺寸宽为L_x和长L_y。

图像上，壁面由灰度值较低的一系列像素点组成的两条近似直线。这两条直线基本是呈水平(x方向)延申的，因此需要沿着竖直方向(y方向)搜索灰度的局部最小值来确定壁面的。搜索一个数组序列中单个极值及其位置已经有成熟的算法，但是同时搜索多个极值及多个位置，还未有成熟算法。下面将对壁面和相界面的识别分别进行阐述。

部分壁面的识别

步骤一，首先通过多极值搜索技术沿着y方向搜索两个壁面最左端(x＝0)的两个点A和B，得到坐标分别为(0,y_A)和(0,y_B)，也就获得了两条壁面的大致位置；

步骤一的具体过程如下：

点A和点B是x＝0的直线上灰度极小值点，也就是灰度梯度沿着y方向存在极大值的两个点，因此找到灰度梯度极大值的两个点的坐标并加以区分，就可以确定点A和点B的坐标。如图5和图6所示；

令任意坐标点(x,y)的灰度值为I(x,y)，则进行以下步骤：

(1)对x＝0的直线上的坐标点灰度值进行向前差分来表示变化梯度dI，dI存在两个极大值点，如图6所示。

(2)设x＝0的直线上所有点的y坐标按照灰度梯度值I(0,y)进行降序排列，得到一个坐标序列(y₁,y₂,...y_n)，坐标y₁是其中一个极大值点的位置，按顺序搜索坐标序列，搜索到第一个满足式(3)的坐标点y_k，其中1<k≤n；式(3)中α为小于1的系数，h为点A到点B的竖直距离h＝y_B-y_A。

|y_k-y₁|>αh(3)

(3)比较y₁和y_k的大小，令二者中较小值为y_min，较大值为y_max，则点A和点B的纵坐标可确定为y_A＝y_min，y_B＝y_max。

步骤二，沿着A-B连线中点的E(0,y_E)的水平方向y＝x_E搜索单个局部最小值，获得相界面上某点F，坐标点为(x_F,y_F)；

步骤二的具体过程如下：

(1)求取A-B连线中点E的纵坐标y_E(横坐标已知x_E＝0)，

(2)将直线y＝y_E上所有坐标点的灰度值I(x,y_E)按式(5)求灰度的梯度，其中x∈[0,L_X]。得到直线y＝y_E上所有坐标点的灰度梯度dI(x,y_E)；

(3)由于dI(x,y_E)只有一个极大值，搜索dI(x,y_E)序列中的极大值点，就可以确定相界面上点F的横坐标x_F，而点F的纵坐标x_F＝x_E。

步骤三，沿着y方向采用已有的排序算法搜索单极值，进而搜索部分壁面，为了避免相界面对壁面的搜索，设置一个范围以避开相界面，令0<x_D1＝x_C1<x_F。在0≤x≤x_D1和0≤y≤y_E范围内搜索得到壁面A-D1，其中D1的坐标为(x_D1,y_D1)；在0<x_f<x_C1和y_E≤y≤L_y范围内搜索得到壁面B-C1，其中C1的坐标为(x_C1,y_C1)，参见图4；

步骤三的具体过程如下：

(1)为了沿着y方向搜索壁面时不受到相界面的影响，限定搜索的x范围为0≤x≤x_D1或0≤x≤x_C1，其中x_D1＝x_C1＝βx_F，β是一个小于1的系数；

(2)为了搜索壁面A-D1，横坐标x_i遍历0至x_D1，对于每个横坐标x_i，求取直线x＝x_i上灰度I(x_i,y)在y方向的梯度，如式(6)所示；搜索dI(x_i,y)上的极大值，可以获得对应壁面坐标点的纵坐标y_i，而一系列坐标点(x_i,y_i)就组成了壁面A-D1；

(3)为了搜索壁面B-C1，横坐标x_j遍历0至x_C1，对于每个横坐标x_j，求取直线x＝x_j上灰度I(x_j,y)在y方向的梯度，将式(6)中x_i替换成x_j计算dI(x_j,y)，搜索dI(x_j,y)上的极大值，可以获得对应壁面坐标点的纵坐标y_j，而一系列坐标点(x_j,y_j)就组成了壁面B-C1。

相界面和接触点的确定：

根据以上步骤得到的结果，可以进一步确定相界面和接触点。紧接着上述的步骤三：

步骤四，沿着x方向采用已有的排序算法搜索单极值，进而搜索部分相界面，为了避免两个壁面对相界面的搜索，设置一个范围以避开两个壁面，令0<y_D2<y_F<y_C2<L_y。在0≤x≤L_x和y_D2≤y≤y_C2范围内搜索得到部分相界面C2-D2，两个点的坐标分别是(x_D2,y_D2)、(x_C2,y_C2)。

步骤四的具体过程如下：

(1)沿x方向搜索灰度的极小值确定相界面，需要避开两个壁面，设定y方向的搜索范围为y_D2<y<y_C2，其中y_D2和y_C2分别为

式(7)中，h为点A到点B的距离，即h＝y_B-y_A，ε是一个系数，其取值范围为0<ε<0.5。

(2)为了搜索相界面，y_i坐标遍历y_D2至y_C2，对于每个纵坐标y_i，求取直线y＝y_i上灰度I(x,y_i)在x方向的梯度dI(x,y_i)，如式(8)所示；搜索dI(x,y_i)上的极大值，可以获得对应壁面坐标点的纵坐标y_i，而一系列坐标点(x_i,y_i)就组成了壁面C2-D2。

步骤五，采用一次函数分别拟合A-D1和B-C1线段内所有的坐标点得到一次函数f_AD和f_BC，采用三次函数拟合C2-D2曲线内的所有点，得到三次函数f_CD。求解f_AD与f_CD的交点可得接触点D的坐标(x_D,y_D)，求解f_AD与f_CD的交点可得接触点C的坐标(x_C,y_C)。线段B-C和A-D则为壁面，曲线C-D为相界面。

步骤五的具体过程如下：

(1)采用一次函数拟合A-D1线段内所有点，以x坐标为自变量，y为因变量，得到的一次函数f_AD如式(9)所示，a_AD和b_AD为一次函数的系数，

y＝f_AD(x)＝a_ADx+b_AD(9)

采用一次函数拟合B-C1线段内所有点，以x坐标为自变量，y为因变量，得到的一次函数f_BC如式(10)所示，a_BC和b_BC为一次函数的系数

y＝f_BC(x)＝a_BCx+b_BC(10)

采用三次函数拟合相界面C2-D2，以y坐标为自变量，x为因变量，得到的三次函数f_CD，如式(11)所示。

x＝f_BC(y)＝a_CDy³+b_CDy²+c_CDy+d_CD(11)

(2)接触点C和D分别是两个壁面与相界面的交点，因此通过求解f_CD与f_AD的交点，及f_CD与f_BC交点，就可以获得接触点C和D的坐标(x_C,y_C)和(x_D,y_D)，上述过程就是求解式(12)所示的方程组。

(3)相界面、壁面上所有坐标点的确定。要求取壁面A-D上的所有坐标点，将每个满足x_A≤x≤x_D的横坐标x带入式(9)中进行计算就可以得到对应y坐标，从而可以获得壁面A-D上所有点的坐标；要求取壁面B-C上的所有坐标点，将每个满足x_B≤x≤x_C的横坐标x带入式(10)中进行计算就可以得到对应y坐标，从而可以获得壁面B-C上所有点的坐标；对于相界面C-D上的所有坐标点，将每个满足y_C≤y≤y_D的纵坐标y带入式(11)中进行计算就可以得到对应x坐标，从而可以获得相界面C-D上所有点的坐标。

接触角的测量：

接触角的定义为相界面在接触点位置的与壁面形成的角度。下面求解相界面在接触点位置的接触角，根据上一步得到的相界面的拟合函数x＝f_CD(y)，将C的坐标(x_C,y_C)以及D的坐标(x_D,y_D)分别代入f_CD的一次导数中，可以求得相界面在C点和D点的接触角θ_C和θ_D。f_CD的一次导数f′_CD为，

f′_CD＝3a_CDy²+2b_CDy+c_CD(13)

将纵坐标y_C和y_D代入式(13)中，可以得到斜率，将斜率转换成角度就可以得到如图6所示的接触角θ_C和θ_D。

另外，考虑实际获取得到的壁面不一定是水平的，可能带有一定倾斜角度，因此要对式(14)获得的接触角进行修正。壁面的倾斜角度为壁面A-D1和B-C1的倾斜角度平均值，即(arctan>AD+arctan>BC)/2。则接触角的最终值为，

本发明提出的方法与已有的商业设备Krüss DSA100接触角仪器的测试结果对比，如下表1所示。

表1本发明的方法与已有的商业设备Krüss DSA100接触角仪器的测试结果对比

从表1可以看出，接触角处在40°至120°范围内的测量偏差小于0.5°，说明本发明的方法测量准确。