一种静态接触角和动态接触角的自动检测方法

摘要

本发明公开了一种自动的高精确度的液滴动静态接触角计算方法。本发明的技术方案是：采用水平集方法获得拍摄的清晰的水珠图像液滴的边缘；判断图像是静态接触角、前进角或后退角图像。若针对静态接触角图像且液滴体积大，选择椭圆拟合，而液滴体积和接触角均较小时，选择圆拟合算法。若计算前进角，则选择圆拟合或椭圆拟合算法；若计算后退角，则选择多项式拟合算法。获得拟合得到的曲线在三重线对应液滴边缘点处的切线进而获得左右2个接触角，二者平均即为最终的接触角。本方法可快速自动计算动静态接触角，准确率高。

说明书

技术领域

本发明涉及一种憎水性检测的静态接触角和动态接触角自动算法，属于材料性能测试领域。

背景技术

瓷和玻璃绝缘子由亲水性材料构成绝缘，在污湿情况下具有较小的表面电阻，导致污闪电压低。硅橡胶具有憎水性、憎水迁移特性，在污湿环境下仍然具有很高的表面电阻，因而具有较高的污闪电压。复合绝缘子和室温硫化(RoomTemperature Vulcanization，RTV)硅橡胶涂层由硅橡胶材料构成，它们具有较强的抗污闪能力，在电力系统外绝缘中有着广泛的应用。憎水性是硅橡胶材料的关键特性，在大气环境中硅橡胶表面容易遭受浸泡、电晕、染污，因此这些情况下其憎水性的丧失和恢复特性的研究具有重要的意义。

硅橡胶的憎水性表征方法主要有喷水定级法和接触角法，前者通过对硅橡胶表面喷水，然后根据水珠的大小、形状和分布来判断表面憎水性状况，它对测试条件要求不高，但判断的主观性稍大，主要针对现场的复合绝缘子或瓷和玻璃绝缘子上的RTV涂层。接触角法通过对硅橡胶表面滴去离子水，然后根据气、液、固三重线处液滴边缘切线与固液交界面的夹角(接触角)来判断材料的憎水性，该法对条件要求稍高，但对材料憎水性的表征更加细致，常用于规则平整的硅橡胶材料。目前针对硅橡胶材料憎水性的研究中常用静态接触角方法，该方法所需图像较少，处理较为方便，在材料表面比较光滑平整时取得了很好的效果。浸泡、电晕、染污后材料表面变得粗糙、有溶质沉淀，此时接触角在某个变化范围内变化时三重线固定不动，也就是说接触角测量存在滞后现象，开始时静态接触角往往偏大而接近前进角，其测试结果与表面憎水性状况未必非常吻合。考虑到一定体积的液滴在粗糙材料的表面，如果用微量进样器增加少量液滴则随着液滴的增加液滴底面积保持不变而接触角逐渐增大，直到液滴体积增加超过临界值后液滴的底面积变大，接触角有所减少，这一过程中得到的最大接触角即为前进角θ_A；如果用微量进样器抽取少量液滴则随着液滴的减少液滴底面积也保持不变而接触角逐渐减少，直到液滴体积减少超过临界值则液滴的底面积减少，接触角增大，这一过程中得到的最小接触角即为后退角θ_R。基于前进角/后退角的动态接触角法一定程度上能抑制粗糙、溶质沉淀等因素对接触角测量造成的干扰。

动态接触角测量时针对一个测试点需要拍摄多张图像才能获得准确的前进角和后退角，如果都通过肉眼选取边缘点然后拟合则工作量过大；而且受到针管的影响水珠的边缘不是完全的圆形或椭圆形，只有选择合适的边缘进行拟合才能获得准确的接触角，这对测试者提出了一定的要求。导致工作量和计算时间倍增。

发明内容

为解决现有接触角计算方法中存在的问题，本发明提供了一种材料憎水性检测时用的静态接触角和动态接触角自动算法，其特征在于，包括：

用数码相机或CCD相机结合图像采集卡拍摄得到水珠图像；

基于水平集的活动轮廓模型获得所述图像中的水珠边缘；

根据圆、椭圆拟合、一阶多项式以及二阶多项式拟合针对刨除针管附近液滴边缘所得误差的大小判断静态接触角、前进角图像或后退角图像；

采用圆拟合法或椭圆拟合算法计算静态接触角或前进角；

采用一阶和二阶多项式拟合算法计算后退角。

所述判断图像为静态接触角、前进角图像或后退角图像：即针对刨除针管附近后的液滴边缘，如果圆拟合算法与椭圆拟合算法的拟合误差小于一阶和二阶多项式拟合误差则为静态接触角或前进角图像；如果圆拟合算法与椭圆拟合算法的拟合误差大于一阶和二阶多项式拟合误差则为后退角图像。

计算静态接触角或前进角，根据拟合误差小于临界值选择待拟合的边缘，用椭圆算法先拟合，当所得接触角大于离心率对应接触角时，对应接触角即为最终的接触角，否则选用圆拟合算法计算接触角。

计算后退角，根据拟合误差小于临界值选择待拟合的边缘，分别选择一阶和二阶多项式拟合，若后者误差大于前者的1.5倍，则选择前者所得接触角为最终的接触角，否则选择后者所得接触角为最终的接触角。

所述各种拟合算法获得拟合误差小于临界值且最长的液滴边缘作为待拟合部分，该边缘在气、固、液交界线以上，包含气、固、液交界线附近的液滴边缘。

本发明的有益效果有：

(1)采用活动轮廓模型获得液滴边缘，抗干扰能力强

因为憎水性测量得到的水珠图像存在的干扰因素很多，常规的边缘算子得到的水珠边缘很容易受到干扰，且得到的边缘点不按顺序排列，不方便于后续处理。水平集方法基于求解偏微分方程，更好地利用了图像信息，能获得全局最优的液滴边缘，对初始轮廓不敏感，因此在各种干扰下均能较为准确获得水珠的边缘，为后续接触角的准确计算奠定了基础。

(2)能完全自动计算动静态接触角

根据圆拟合、椭圆拟合误差与多项式拟合误差大小的比较自动判断是静态接触角、前进角计算还是后退角计算，本发明无需提前知道待处理的图像是静态接触角图像还是动态接触角图像。针对静态接触角时的液滴图像，应选择这个液滴的整个边缘进行拟合，针对动态接触角图像应选择靠近三重线附近的液滴边缘进行拟合，这样才能更好利用液滴边缘信息准确获得接触角。本发明根据拟合误差的大小自适应选择待拟合的边缘，即选择拟合误差小于临界值时最长的待拟合边缘，这样准确性较高。本发明既无需知道待处理的图像属于静态接触角图像还是动态接触角图像，同时也无需提前获得液滴边缘信息。

(3)计算准确性高

针对静态接触角和前进角的液滴图像时由于本发明根据离心率与接触角的关系自动选择准确性更高的圆拟合算法或椭圆拟合算法，计算后退角时本发明根据拟合误差自适应选择线性拟合或二阶多项式拟合，这样能有效提高动静态接触角计算的准确性。

附图说明

下面结合附图对本发明作详细说明：

图1为根据本发明实施例的流程图；

图2为根据本发明实施例的圆拟合算法的流程图；

图3为根据本发明实施例的椭圆拟合算法的流程图；

图4为根据本发明实施例的圆与椭圆拟合算法选择图。

具体实施方式

参见附图1，应用本发明的步骤为：

(1)采用数码相机或CCD相机结合图像采集卡对材料上的水珠进行拍照，相机平面垂直于材料平面；

(2)将所得水珠图像采用基于水平集法获得水珠边缘；

(3)根据圆、椭圆拟合与多项式拟合针对刨除针管附近液滴边缘进行拟合，根据所得误差的大小判断为静态接触角、前进角图像还是后退角图像，如果前者小于后者则为静态接触角或前进角图像；否则为后退角图像；

(4)当计算静态接触角或前进角时选择圆拟合法或椭圆拟合法，根据拟合误差小于临界值选择待拟合的边缘，用椭圆算法先拟合，当所得接触角大于离心率对应接触角时对应接触角即为最终的接触角，否则选用圆拟合算法计算接触角。计算后退角时根据拟合误差小于临界值选择待拟合的边缘，分别选择一阶和二阶多项式拟合，考虑若后者误差大于前者的1.5倍，则选择前者所得接触角为最终的接触角，否则选择后者所得接触角为最终的接触角。

水平集函数曲面的演化遵循如下的Haminlton-Jacobi方程。

式中F表示曲线上各点的演化速度，方向沿着曲线的法线方向，通常与图像梯度和曲率有关；为水平集函数在二维平面上的梯度；不同的F表达式对应了不同的图像分割方法。

Chan-Vese模型的能量函数表示如下：

$>E(c_1,c_2,C)=\mathrm{\mu L}\left(C\right)+\mathrm{vA}\left(C\right)+{\lambda }_1\underset{I\left(C\right)}{\int }{|u_0-c_1|}^2\mathrm{dxdy}+{\lambda }_2\underset{O\left(C\right)}{\int }{|u_0-c_2|}^2\mathrm{dxdy}---\left(2\right)$>

式中A(C)为曲线C的内部区域面积；μ，v≥0，λ₁，λ₂＞0为权重系数；I(C)和O(C)分别为C内部和外部区域；c₁、c₂分别为I(C)和O(C)上图像灰度的均值。

和定义如下：

式中为水平集函数。

能量函数表示为

最小化能量函数将得到欧拉-拉格朗日方程，转换后得到水平集函数的梯度如下：

式(6)经过计算后即为所得图像的边缘。

在使用水平集获得图像中液滴边缘曲线后，如果是静态接触角图像，此时微量进样器的针管并没有插入到液滴中，所得图像中液滴的边缘为圆或椭圆，考虑到当液滴体积不是很大时液滴图像边缘满足圆形，使用基于圆的拟合算法，随着液滴体积的增加其边缘更接近于椭圆而选择椭圆拟合算法。如果是动态接触角图像，由于微量进样器的针管插入到水珠中，液滴的边缘不能完全满足圆或椭圆方程，当测量前进角时三重线(即气、液、固三态交界线)附近液滴边缘接近于圆或椭圆，根据离心率与接触角的关系选择圆拟合或椭圆拟合算法，然后根据误差小于临界值自适应选择最长的待拟合液滴边缘，根据拟合所得曲线获得最终接触角。如果测量后退角，三重线附近液滴边缘一般接近于多项式曲线，根据拟合误差小于临界值自适应选择最长的待拟合边缘，拟合获得曲线方程后最终获得接触角。

-基于圆的拟合方法

设获得的水珠边缘点横、纵坐标组成的数组为X(n)，n＝1，...，2i-1，2i，...2N，其中X(2i-1)、X(2i)分别为第i点的横、纵坐标；水珠边缘所在圆的圆心坐标为[X₀，Y₀]，半径为R，则定义第n点的误差如下：

$>e_n=\sqrt{{(X(2n-1)-X_0)}^2+{(X\left(2n\right)-Y_0)}^2}-R---\left(7\right)$>

则所有点总误差定义为：

$>E=\frac{1}{2}\underset{n=1}{\overset{N}{\Sigma }}{e}_n^2---\left(8\right)$>

式中X₀、Y₀、R为待求变量。Levenberg-Marquardt算法用一阶偏导的计算量获得了接近二阶偏导的计算速度，非常适合于非线性最小二乘问题，因此选用该方法。

设e＝[e₁，e₂，...，e_N]^T为误差列向量；W为非线性多元函数变量组成的列向量，W＝[X₀，Y₀，R]^T；J为雅克比矩阵，w_j为W中第j个元素；I为3维单位阵。可变参数迭代公式如下：

W(k+1)＝W(k)-(J(k)^TJ(k)+λI)^-1J(k)^Te(k)        (9)

式中k为迭代次数；λ根据前后两次计算误差的比较结果进行调整，如果误差增加则λ＝λ×10，如果误差减少则λ＝λ×0.1，其初值选择0.1效果不错。

对于静态接触角的拟合来说，雅克比矩阵计算如下：

$>J_{i,1}=\frac{{\partial e}_i}{\partial X_O}=\frac{-(X(2n-1)-X_O)}{\sqrt{{(X(2n-1)-X_O)}^2+{(X\left(2n\right)-Y_O)}^2}}---\left(10\right)$>

$>J_{i,2}=\frac{{\partial e}_i}{\partial Y_O}=\frac{-(X\left(2n\right)-Y_O)}{\sqrt{{(X(2n-1)-X_O)}^2+{(X\left(2n\right)-Y_O)}^2}}---\left(11\right)$>

$>J_{i,3}=\frac{{\partial e}_i}{\partial R}=-1---\left(12\right)$>

式中n为水珠边缘点的序号。

因为初始解与最优解的差距将严重影响收敛速度，甚至计算精度。气、液、固三态的交界线称为三重线，其在图像对应为水珠边缘中纵坐标最小的点，左、右侧各一个点，设序号分别为j、k。因硅橡胶表面水珠的接触角通常在90°周围分布，基于半圆的假设，算法使用的初始圆心和半径使用如下策略获得：

X₀＝(X(2j-1)+X(2k-1))/2；Y₀＝(X(2j)+X(2k))/2                    (13)

$>R=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}\sqrt{{(X(2i-1)-X_0)}^2+{(X\left(2i\right)-Y_0)}^2}---\left(14\right)$>

根据该方法能快速获得圆心和半径的初始值，且该值偏离准确值通常不大，实测结果表明该方式能保障测量的准确性和实时性。

设水珠左右两侧边缘与水平面交点坐标分别为(x₁，y₁)、(x₂，y₂)，拟合得到的圆心坐标为[X₀，Y₀]，则左右两侧斜率计算如下：

k₁＝-(x₁-X₀)/(y₁-Y₀)；k₂＝-(x₂-X₀)/(y₂-Y₀)                   (15)

式中k₁和k₂分别为三重线处圆弧上左、右切线的斜率。左、右两侧接触角计算公式如下。

θ₁＝atan(k₁)×180/π，k₁≥0；θ₁＝180+atan(k₁)×180/π，k₁＜0；

θ₂＝180-atan(k₂)×180/π，k₂≥0；θ₂＝-atan(k₂)×180/π，k₂＜0        (16)

式中θ₁和θ₂分别为左、右两侧的接触角，atan函数的取值范围为-π/2～π/2。

水珠的接触角θ计算如下：

θ＝(θ₁+θ₂)/2                      (17)

式中θ的单位为°。

Levenberg-Marquardt算法通过迭代来提高精确度，判断收敛的情况对算法的精确度和实时性有明显影响。算法根据迭代前后所得接触角的变化情况来决定是否停止计算。Levenberg-Marquardt算法迭代过程中连续三次计算结果往往能很好表征收敛情况。因此，收敛准则如下：设第N、N+1、N+2次迭代后所得接触角分别为A_N、A_N+1、A_N+2，如果满足

|A_N-A_N+1|≤C₁和|A_N+1-A_N+2|≤C₂              (18)

则迭代终止，其中C₁、C₂分别为设定的临界值，本发明选择为0.5°和0.01°，具体应用时可以根据实际情况适当调整，正常情况下使用该准则计算所得接触角的标准差仅为0.5°左右。算法流程图如附图2所示。

-基于椭圆的拟合方法

当液滴体积较大时圆拟合算法误差增加，应选择椭圆方程对液滴边缘进行拟合。常规最小二乘拟合计算量大，接触角计算的实时性将受到影响，选择一种直接计算椭圆参数的算法。椭圆的一般方程可表示为：

F(m，n)＝n·m＝ax²+bxy+cy²+dx+ey+f＝0                 (19)

上式中x、y分别表示横、纵坐标，m＝[a，b，c，d，e，f]^T，n＝[x²，xy，y²，x，y，1]。设n_i＝[x_i²，x_iy_i，y_i²，x_i，y_i，1]，B＝[n₁^T，n₂^T，...，n_N^T]^T，F(m，n_i)称之为平面上点(x_i，y_i)到曲线F(m，n)＝0的代数距离。当椭圆拟合所有离散数据点代数距离平方和最小时即可求解出相应的二次曲线，这属于非线性最小二乘问题，常规的最小二乘算法如Levenberg-Marquardt算法计算量和编程工作量都较大。对于N点的观测数据，拟合准则为：

$>E_{\min }=\min \left(\underset{i=1}{\overset{N}{\Sigma }}F{(m,n_i)}^2\right)---\left(20\right)$>

针对式(20)必须限定b²-4ac＜0才能保证拟合的结果为椭圆，否则拟合结果有可能是抛物线或者双曲线而非椭圆。因b²-4ac＜0不是一个等式限制条件，在实际求解时由于Kuhn-Tucker条件并不能保证有解，故引入限制条件b²-4ac＝-1，以矩阵形式表达为：

m^TCm＝1                (21)

式中$>C=\left(\begin{array}{cccccc}0& 0& 2& 0& 0& 0\\ 0& -1& 0& 0& 0& 0\\ 2& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right).$>

式(20)等价为

E_min＝min(|Bm|²)                  (22)

因为B是已知量，式(22)中待求量为m，问题的关键就是对m进行寻优。

引入Lagrange算子并求导得：

2B^TBm-2λCm＝0                 (23)

令S＝B^TB，式(23)改写为：

Sm＝λCm                       (24)

对于式(24)，可按广义特征值和广义特征向量的方法求解出6组对应的特征值λ_i和u_i。对于条件限制矩阵C，其特征值为[-2，-1，2，0，0，0]，仅有一个特征值为正。仅有惟一的广义特征值λ∈R⁺和广义特征向量u作为椭圆拟合解。对于任意k∈R⁺而言，应使(λ，ku)满足式(21)，即k²u^TCu＝1，则：

k＝[1/(u^TCu)]^1/2＝[1/(u^TSu)]^1/2                   (25)

m的唯一解：

$>\hat{m}=\mathrm{ku}---\left(26\right)$>

根据所得解即可获得椭圆的中心、长轴、短轴、旋转角度等信息，然后根据切线情况可获得接触角。

算法流程图如附图3所示。设根据求解的m获得椭圆的长、短半轴分别为L_L和L_S，中心为(X₀，Y₀)，倾斜角度为θ₀弧度；则椭圆上点(x₁，y₁)表达式如下：

x＝L_Lcosθ；y＝L_Ssinθ；              (27)

x₁＝xcosθ₀-ysinθ₀+X₀；

y₁＝xsinθ₀+ycosθ₀+Y₀；

液滴边缘上左右三重线对应点分别为(X_L，Y_L)和(X_R，Y_R)，考虑到样本可能不处于完全水平的状态，接触角计算思路如下，先获得左右三重线对应点对应的式(27)中的θ，记为θ_L1和θ_R1，angle(A)获得复数A的相位，函数的取值范围为-π～π。计算如下：

X_L1＝(X_L-X₀)cosθ₀+(Y_L-Y₀)sinθ₀；

Y_L1＝-(X_L-X₀)sinθ₀+(Y_L-Y₀)cosθ₀；

X_R1＝(X_R-X₀)cosθ₀+(Y_R-Y₀)sinθ₀；

Y_R1＝-(X_R-X₀)sinθ₀+(Y_R-Y₀)cosθ₀；

θ_L1＝angle(Y_L1/L_Sj+X_L1/L_L)；

θ_R1＝angle(Y_R1/L_Sj+X_R1/L_L).

根据θ_L1和θ_R1以及椭圆的倾斜角度θ₀可得不考虑固体平面倾斜时液滴左右两侧的接触角分别为θ_L2和θ_R2，单位为°，如式(28)所示。

$>\left(\begin{array}{c}{\theta }_{L2}=\frac{180}{\pi }[a\tan \left(\frac{-L_S}{L_L\tan {\theta }_{L1}}\right)+{\theta }_0],Y_{L1}\ge 0\\ {\theta }_{L2}=\frac{180}{\pi }[\pi +a\tan \left(\frac{-L_S}{L_L\tan {\theta }_{L1}}\right)+{\theta }_0],Y_{L1}<0\\ {\theta }_{R2}=\frac{180}{\pi }[a\tan \left(\frac{L_S}{L_L\tan {\theta }_{R1}}\right)+{\theta }_0],Y_{R1}\ge 0\\ {\theta }_{R2}=\frac{180}{\pi }[\pi +a\tan \left(\frac{L_S}{L_L\tan {\theta }_{R1}}\right)+{\theta }_0],Y_{R1}<0\end{array}\right)---\left(28\right)$>

设根据左右三重线对应点获得固体平面斜率为k₁，则左右接触角和最终接触角分别为θ_L、θ_R和θ，如式(29)、(30)所示，单位为°。

$>\left(\begin{array}{c}{\theta }_L={\theta }_{L2}-180a\tan \left(k_1\right)/\pi \\ {\theta }_R={\theta }_{R2}+180a\tan \left(k_1\right)/\pi \end{array}\right)---\left(29\right)$>

θ＝(θ_L+θ_R)/2                                              (30)

-多项式拟合法

多项式拟合法通过在三重线附近液滴边缘上选择几个点，对其进行多项式拟合，然后获得三重线处的切线进而获得接触角。设根据拟合获得左、右两侧多项式为而左右三重线处坐标为(X_L，Y_L)，(X_R，Y_R)，则根据多项式其对应在液滴上的切线的斜率分别为设根据左右三重线处获得固体平面斜率为k₁，则左右2侧接触角为

$>\left(\begin{array}{c}{\theta }_L=\frac{180}{\pi }[a\tan k_L-a\tan k_1],Y_{L1}\ge 0\\ {\theta }_L=\frac{180}{\pi }[\pi +a\tan k_L-a\tan k_1],Y_{L1}<0\\ {\theta }_R=\frac{180}{\pi }[a\tan k_R+a\tan k_1],Y_{R1}\ge 0\\ {\theta }_R=\frac{180}{\pi }[\pi +a\tan k_R+a\tan k_1],Y_{R1}<0\end{array}\right)$>

液滴的最终接触角

θ＝(θ_L+θ_R)/2                    (31)

-边缘与算法的自适应选择

当图像为静态接触角图像或前进角图像时除针管附近边缘外，整个液滴边缘可用圆或椭圆拟合算法较好表达；当图像为后退角图像则除针管附近边缘外，整个液滴边缘可用一阶或二阶多项式拟合算法较好表达。因此，根据圆、椭圆拟合与多项式拟合针对刨除针管附近液滴边缘误差的大小判断为静态接触角、前进角图像还是后退角图像，如果前者小于后者则图像为静态接触角或前进角图像；否则为后退角图像。当水珠体积较小时，当针管不在水珠内部则其边缘呈圆形，可将获得的边缘用圆方程拟合即可。当针管插入到水珠中，它对针管附近的水珠边缘有影响，水珠边缘不再满足圆形，但测量前进角时三重线附近的水珠边缘影响很小，如果拟合边缘选择合适则拟合误差较小，如果边缘选择过大则拟合误差变大，如果边缘选择过小则拟合误差虽小，但拟合得到曲线未必接近水珠的真实边缘，可以在误差小于临界值时选择尽量多的待拟合点。拟合曲线可以是圆形或多项式，测量前进角时接触角较大、待拟合的弧线较长时该拟合曲线比较接近于圆或椭圆，因此选择圆或椭圆拟合有较好效果。测量静态接触角和前进角时具体准则可根据附图4判断，当离心率确定时接触角大于对应的临界角时选择椭圆拟合算法否则选择圆拟合算法。测量后退角时接触角较小，弧线较短时，待拟合部分比较接近于低阶多项式，故此时选择多项式拟合，多项式一般根据实际情况选择二阶或一阶，考虑若后者误差大于前者的1.5倍，则选择前者，否则选择后者。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。