CN109408909A - 一种三维粗糙表面微凸体拟合方法 - Google Patents

Abstract

本发明公开了一种三维粗糙表面微凸体拟合方法，包括以下步骤：a.测量待拟合粗糙表面的表面形貌离散点高度矩阵z(i,j)；b.计算得到表面平均高度面h₀，取其为拟合参考面；c.对表面形貌离散点进行分割，划分为若干封闭区域；d.根据拟合参考面以及表面形貌离散点高度之间的相互关系，对每一个封闭区域内的离散点进行判断，将符合微凸体特征的离散点序列判定为微凸体；e.以步骤d中判定得到的微凸体，根据最小均方误差，用椭球面对每个微凸体的离散点进行拟合，得到拟合后椭球方程，据此计算微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等表面微凸体参数。本发明的一种三维粗糙表面微凸体拟合方法，具有测量精确的优点。

Description

一种三维粗糙表面微凸体拟合方法

技术领域

本发明涉及一种三维表面拟合方法，尤其涉及一种三维粗糙表面微凸体拟合方法。

背景技术

机械加工零件表面在微观尺度上都是凹凸不平的粗糙表面，当两表面接触的时候，实际产生接触的是一系列离散分布的接触区域，接触面积仅是光滑名义接触面积下的一部分，远小于名义接触面积，这会导致真实接触压力远大于名义接触压力。因此，接触表面的真实接触面积、接触压力等直接影响着传动部件的承载能力以及摩擦、磨损、疲劳等性能。

基于赫兹接触理论与数理统计原理，考虑微凸体间接触的粗糙表面接触模型目前广泛运用于粗糙表面接触分析，粗糙表面的表面微凸体参数：微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等对接触分析计算结果影响极大，然而考虑到光学仪器测量的表面形貌受到采样间距的影响，其微凸体参数计算同样将受到采样间距的巨大影响，通过经典三点法、五点法等计算微凸体参数得到的微凸体平均曲率半径通常过小，不能合适的用于粗糙表面微凸体接触模型计算，这极大影响了接触分析结果的准确性。为了提高参数计算的准确性，减小采样间距的影响，有研究者采用抛物线对表面微凸体进行拟合简化，再计算表面微凸体参数，在一定程度上提高了计算准确性。

然而现有拟合计算方法局限于二维粗糙轮廓的拟合简化，而粗糙表面接触模型是基于三维表面的接触发展而来，二维拟合简化得到的微凸体与实际三维微凸体存在显著的差别，简单的二维轮廓无法完全准确的表征粗糙表面的三维形貌，因此，现有方法的准确性还有待提高。

发明内容

本发明所要解决的技术问题是提供一种拟合误差小并且精确度高的三维粗糙表面微凸体拟合方法。

本发明是通过以下的技术方案实现的：

一种三维粗糙表面微凸体拟合方法，包括以下步骤，

a.测量待拟合粗糙表面的表面形貌离散点高度矩阵z(i,j)；

b.以步骤a中测量得到的表面形貌离散点高度序列z(i,j)为基础，计算得到表面平均高度面h₀，取其为拟合参考面；

c.基于分水岭算法，在Matlab软件中对表面形貌离散点进行分割，划分为若干封闭区域；

d.根据拟合参考面以及表面形貌离散点高度之间的相互关系，对每一个封闭区域内的离散点进行判断，将符合微凸体特征的离散点序列判定为微凸体；

e.以步骤d中判定得到的微凸体，根据最小均方误差，用椭球面对每个微凸体的离散点进行拟合，得到拟合后椭球方程，据此计算微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等表面微凸体参数。

优选的，a步骤中采用白光干涉仪进行测量。

优选的，b步骤中计算表面平均高度面h₀包括以下步骤：

(b1)根据表面形貌高度离散点矩阵z(i,j)计算得到平均高度线h₀，计算公式为

式中，N为表面形貌离散点垂直方向的个数，M为表面形貌离散点水平方向的个数；

(b2)令拟合参考面为高度为h＝h₀的平面。

优选的，c步骤中，包括以下步骤：

(c1)计算表面十点高度S_z，S_z定义为表面形貌离散点中，5个最高顶点的高度和5个最深谷点的深度的平均值，计算公式为

式中，z_maxk和z_mink(k＝1,2,...,5)为5个最高顶点的高度和5个最深谷点的深度，若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最大值，则该点为顶点，同理若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最小值，则该点为谷点；

(c2)根据拟合参考面对表面形貌离散点矩阵z(i,j)进行处理，仅对高于拟合参考面的表面进行分割处理，若z(i,j)>h，则令z(i,j)＝z(i,j)-h，若z(i,j)<＝h，则令z(i,j)＝h；

(c3)使用Matlab对表面形貌离散点矩阵z(i,j)进行垂直方向和水平方向的Sobel算子滤波，然后求取模值，得到表面形貌离散点的梯度幅值矩阵gradz(i,j)，计算公式为：

式中，I_x(i,j)为使用Sobel算子对表面形貌离散点矩阵进行水平方向边缘滤波的卷积结果，I_y(i,j)为使用Sobel算子对表面形貌离散点矩阵进行垂直方向边缘滤波的卷积结果；

(c4)根据表面十点高度S_z对梯度幅值矩阵gradz(i,j)进行处理，若gradz(i,j)<0.05S_z，则令gradz(i,j)＝0；

(c5)使用Matlab中的watershed分水岭算法函数对梯度幅值矩阵进行分水岭分割，将表面分割为若干封闭区域。

优选的，d步骤中，包括以下步骤：

(d1)使用Matlab中的regiongrops图像处理函数对分水岭算法分割后的结果进行统计处理，得到每个封闭区域内所有离散点的坐标和区域边界上离散点的坐标；

(d2)针对单个封闭区域，若该区域内所有离散点高度的平均值大于该区域边界上离散点高度的平均值，且该区域内所有离散点高度的平均值大于拟合参考面高度h，则将该封闭区域内的离散点判定为组成一个微凸体。

优选的，e步骤中，包括以下步骤：

(e1)针对判定后单个微凸体的封闭区域，根据Matlab中的regiongrops图像处理函数可确定与该封闭区域具有相同标准二阶中心距的椭圆，令a为与该区域具有相同标准二阶中心矩的椭圆的长轴长度，b为与该区域具有相同标准二阶中心矩的椭圆的短轴长度，O点为该区域的质心；

(e2)针对(e1)中确定的椭圆区域，以O点为坐标系原点，椭圆长轴方向为x轴，椭圆短轴方向为y轴，高度方向为z轴建立坐标系，设区域内共有n个离散点，离散点高度分别为z(x₁,y₁)，…，z(x_i,y_i)，…，z(x_n,y_n)；

(e3)用椭球面对该椭圆区域内离散点进行拟合，令椭球面方程为：

式中，c为椭球面系数；

(e4)运用最小均方根误差来判定最适合的拟合椭球面，拟合椭球面最小均方误差Er由以下公式确定：

令即可计算得到每个椭球面的系数c：

将得到的椭球面系数c代入椭球面方程，即可得到微凸体离散点的拟合椭球面方程，确定拟合形貌。

优选的，e步骤中，包括以下步骤：

(e5)每一个椭球顶端曲率半径即为微凸体顶峰曲率半径，椭球在x轴方向等效曲率半径为R_x，在y轴方向等效曲率半径为R_y，计算公式如下：

设共有m个微凸体，则微凸体平均曲率半径R为：

(e6)微凸体峰点高度标准偏差σ、微凸体密度η为：

η＝m/A (11)

式中，A为整个粗糙表面的名义面积。

有益效果是：

与现有技术相比，本发明的一种三维粗糙表面微凸体拟合方法通过实测工件的表面形貌，提取表面三维形貌离散点，基于分水岭算法对表面区域进行分割，通过离散点高度与初始参考面对分割区域进行判定，进一步确定拟合微凸体的区域，根据最小均方误差对每个微凸体进行椭球面拟合，方便计算微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等微凸体参数，拟合得到的椭球微凸体更加贴近实际接触的微凸体，使得拟合误差小，能更加准确的表征三维接触计算中的表面形貌，在一定程度上避免了采样间距对微凸体参数计算的巨大影响，使微凸体参数计算更加合理，以此计算得到微凸体参数能更加准确的表征接触微凸体的几何特征，提高了粗糙表面接触分析的准确性。

附图说明

以下结合附图对本发明的具体实施方式作进一步的详细说明，其中：

图1为工件表面形貌示意图；

图2为分割结果局部图；

图3为单个微凸体拟合前原始形貌示意图；

图4为单个微凸体拟合后形貌示意图；

图5为工件表面形貌拟合结果图。

具体实施方式

一种三维粗糙表面微凸体拟合方法，包括以下步骤：

a.测量待拟合粗糙表面的表面形貌离散点高度矩阵z(i,j)；

b.以步骤a中测量得到的表面形貌离散点高度序列z(i,j)为基础，计算得到表面平均高度面h₀，取其为拟合参考面；

c.基于分水岭算法，在Matlab软件中对表面形貌离散点进行分割，划分为若干封闭区域；

d.根据拟合参考面以及表面形貌离散点高度之间的相互关系，对每一个封闭区域内的离散点进行判断，将符合微凸体特征的离散点序列判定为微凸体；

e.以步骤d中判定得到的微凸体，根据最小均方误差，用椭球面对每个微凸体的离散点进行拟合，得到拟合后椭球方程，据此计算微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等表面微凸体参数。

较佳的实施方式，a步骤中采用白光干涉仪进行测量。

较佳的实施方式，b步骤中计算表面平均高度面h₀包括以下步骤：

(b1)根据表面形貌高度离散点矩阵z(i,j)计算得到平均高度线h₀，计算公式为

式中，N为表面形貌离散点垂直方向的个数，垂直方向也即是x方向，M为表面形貌离散点水平方向的个数，水平方向也即是y方向，表面形貌图如图1所示；

(b2)令拟合参考面为高度为h＝h₀的平面。

较佳的实施方式，c步骤中，包括以下步骤：

(c1)计算表面十点高度S_z，S_z定义为表面形貌离散点中，5个最高顶点的高度和5个最深谷点的深度的平均值，计算公式为

式中，z_maxk和z_mink(k＝1,2,...,5)为5个最高顶点的高度和5个最深谷点的深度，若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最大值，则该点为顶点，同理若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最小值，则该点为谷点；

(c2)根据拟合参考面对表面形貌离散点矩阵z(i,j)进行处理，仅对高于拟合参考面的表面进行分割处理，若z(i,j)>h，则令z(i,j)＝z(i,j)-h，若z(i,j)<＝h，则令z(i,j)＝h；

(c3)使用Matlab对表面形貌离散点矩阵z(i,j)进行垂直方向和水平方向的Sobel算子滤波，然后求取模值，得到表面形貌离散点的梯度幅值矩阵gradz(i,j)，计算公式为：

式中，I_x(i,j)为使用Sobel算子对表面形貌离散点矩阵进行水平方向边缘滤波的卷积结果，I_y(i,j)为使用Sobel算子对表面形貌离散点矩阵进行垂直方向边缘滤波的卷积结果；

(c4)根据表面十点高度S_z对梯度幅值矩阵gradz(i,j)进行处理，若gradz(i,j)<0.05S_z，则令gradz(i,j)＝0；

(c5)使用Matlab中的watershed分水岭算法函数对梯度幅值矩阵进行分水岭分割，将表面分割为若干封闭区域，表面分割结果局部如图2所示。

较佳的实施方式，d步骤中，包括以下步骤：

(d1)使用Matlab中的regiongrops图像处理函数对分水岭算法分割后的结果进行统计处理，得到每个封闭区域内所有离散点的坐标和区域边界上离散点的坐标；

(d2)针对单个封闭区域，若该区域内所有离散点高度的平均值大于该区域边界上离散点高度的平均值，且该区域内所有离散点高度的平均值大于拟合参考面高度h，则将该封闭区域内的离散点判定为组成一个微凸体。

较佳的实施方式，e步骤中，包括以下步骤：

(e1)针对判定后单个微凸体的封闭区域，根据Matlab中的regiongrops图像处理函数可确定与该封闭区域具有相同标准二阶中心距的椭圆，令a为与该区域具有相同标准二阶中心矩的椭圆的长轴长度，b为与该区域具有相同标准二阶中心矩的椭圆的短轴长度，O点为该区域的质心；

(e2)针对(e1)中确定的椭圆区域，以O点为坐标系原点，椭圆长轴方向为x轴，椭圆短轴方向为y轴，高度方向为z轴建立坐标系，设区域内共有n个离散点，离散点高度分别为z(x₁,y₁)，…，z(x_i,y_i)，…，z(x_n,y_n)；

(e3)用椭球面对该椭圆区域内离散点进行拟合，令椭球面方程为：

式中，c为椭球面系数；

(e4)运用最小均方根误差来判定最适合的拟合椭球面，拟合椭球面最小均方误差Er由以下公式确定：

令即可计算得到每个椭球面的系数c：

将得到的椭球面系数c代入椭球面方程，即可得到微凸体离散点的拟合椭球面方程，确定拟合形貌，表面拟合结果如图5所示。

较佳的实施方式，e步骤中，包括以下步骤：

(e5)每一个椭球顶端曲率半径即为微凸体顶峰曲率半径，椭球在x轴方向等效曲率半径为R_x，在y轴方向等效曲率半径为R_y，计算公式如下：

设共有m个微凸体，则微凸体平均曲率半径R为：

(e6)微凸体峰点高度标准偏差σ、微凸体密度η为：

η＝m/A (11)

式中，A为整个粗糙表面的名义面积。

以下，根据上述的方法提供具体算例：

用白光干涉仪Wyko NT9100对平板工件表面形貌进行测量，选择5×镜头倍数，采样间隔1μm，得到表面形貌离散点高度矩阵z(i,j)，共472×632个离散点；

一、确定拟合参考面

根据表面形貌高度离散点矩阵z(i,j)计算表面平均高度h₀，数据代入式(1)：

二、分割表面区域

若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最大值，则该点为顶点，同理若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最小值，则该点为谷点，以此得到所有顶点与谷点深度，计算表面十点高度S_z，S_z为5个最高顶点的高度和5个最深谷点的深度的平均值，数据代入式(2)：

再根据拟合参考面对表面形貌离散点矩阵z(i,j)进行处理，若z(i,j)>0.0075，则令z(i,j)＝z(i,j)-0.0075，若z(i,j)<＝0.0075，则令z(i,j)＝0.0075；

然后，使用Matlab对表面形貌离散点矩阵z(i,j)进行垂直方向和水平方向的Sobel算子滤波，然后求取模值，得到表面形貌离散点的梯度幅值矩阵gradz(i,j)：

式中，I_x(i,j)为使用Sobel算子对表面形貌离散点矩阵进行水平方向边缘滤波的卷积结果，I_y(i,j)为使用Sobel算子对表面形貌离散点矩阵进行垂直方向边缘滤波的卷积结果；

再根据表面十点高度S_z对梯度幅值矩阵进行处理，若gradz(i,j)<0.05S_z，即如果gradz(i,j)<0.275,则令gradz(i,j)＝0；

最后使用Matlab中的watershed分水岭算法函数对梯度幅值矩阵进行分水岭分割，将表面分割为若干封闭区域，分割结果局部如图2所示。

三、微凸体判定

使用Matlab中的regiongrops图像处理函数对分水岭算法分割后的结果进行统计处理，得到每个封闭区域内所有离散点的坐标和区域边界上离散点的坐标。

然后针对单个封闭区域，若该区域内所有离散点高度的平均值大于该区域边界上离散点高度的平均值，且该区域内所有离散点高度的平均值大于拟合参考面高度h，则将该封闭区域内的离散点判定为组成一个微凸体。

四、微凸体拟合

首先，针对判定后单个微凸体的封闭区域，根据Matlab中的regiongrops图像处理函数可确定与该封闭区域具有相同标准二阶中心距的椭圆，令a为与该区域具有相同标准二阶中心矩的椭圆的长轴长度，b为与该区域具有相同标准二阶中心矩的椭圆的短轴长度，O点为该区域的质心。任取一微凸体为例，如一区域内表面微凸体形貌如图3所示，该区域质心为点O(279,37)，与该区域具有相同标准二阶中心矩的椭圆的长轴长度a＝8.7μm，与该区域具有相同标准二阶中心矩的椭圆的短轴长度b＝6.4μm；

然后，针对上段中确定的椭圆区域，以点O(279,37)为坐标系原点，椭圆长轴方向为x轴，椭圆短轴方向为y轴，高度方向为z轴建立坐标系，则区域内共有46个离散点，离散点高度分别为z(x₁,y₁)，…，z(x_i,y_i)，…，z(x₄₆,y₄₆)；

再用椭球面对该椭圆区域内离散点进行拟合，数据代入式(4)，令椭球面方程为：

式中，c为椭球面系数；

然后，运用最小均方根误差来判定最适合的拟合椭球面，数据代入式(5)，拟合椭球面最小均方误差Er由以下公式确定：

令即可计算得到每个椭球面的系数c，数据代入式(6)，可得：

将得到的椭球面系数c代入椭球面方程式(15)，即可得到微凸体离散点的拟合椭球面方程，确定拟合形貌，单个微凸体拟合后形貌如图4所示，整体表面拟合结果如图5所示。

五、微凸体参数的计算

每一个椭球顶端曲率半径即为微凸体顶峰曲率半径，椭球在x轴方向等效曲率半径为Rx，在y轴方向等效曲率半径为R_y，数据分别代入式(7)、(8)，如上述微凸体拟合后顶峰曲率半径如下：

分割并判定后共有859个微凸体，数据代入式(9)，则微凸体平均曲率半径R为：

将数据分别代入式(10)、(11)，即可获得微凸体峰点高度标准偏差σ、微凸体密度η为：

η＝859/(472×632)＝0.00299/μm² (22)。

以上实施例仅用以说明本发明的技术方案而并非对其进行限制，凡未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明技术方案的范围内。

Claims (7)

1.一种三维粗糙表面微凸体拟合方法，其特征在于：包括以下步骤，

a.测量待拟合粗糙表面的表面形貌离散点高度矩阵z(i,j)；

b.以步骤a中测量得到的表面形貌离散点高度序列z(i,j)为基础，计算得到表面平均高度面h₀，取其为拟合参考面；

c.基于分水岭算法，在Matlab软件中对表面形貌离散点进行分割，划分为若干封闭区域；

d.根据拟合参考面以及表面形貌离散点高度之间的相互关系，对每一个封闭区域内的离散点进行判断，将符合微凸体特征的离散点序列判定为微凸体；

e.以步骤d中判定得到的微凸体，根据最小均方误差，用椭球面对每个微凸体的离散点进行拟合，得到拟合后椭球方程，据此计算微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等表面微凸体参数。

2.根据权利要求1所述的三维粗糙表面微凸体拟合方法，其特征在于：a步骤中采用白光干涉仪进行测量。

3.根据权利要求1所述的三维粗糙表面微凸体拟合方法，其特征在于，b步骤中计算表面平均高度面h₀包括以下步骤：

(b1)根据表面形貌高度离散点矩阵z(i,j)计算得到平均高度线h₀，计算公式为：

式中，N为表面形貌离散点垂直方向的个数，M为表面形貌离散点水平方向的个数；

(b2)令拟合参考面为高度为h＝h₀的平面。

4.根据权利要求1所述的三维粗糙表面微凸体拟合方法，其特征在于，c步骤中，包括以下步骤：

(c1)计算表面十点高度S_z，S_z定义为表面形貌离散点中，5个最高顶点的高度和5个最深谷点的深度的平均值，计算公式为

式中，z_maxk和z_mink(k＝1,2,...,5)为5个最高顶点的高度和5个最深谷点的深度，若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最大值，则该点为顶点，同理若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最小值，则该点为谷点；

(c2)根据拟合参考面对表面形貌离散点矩阵z(i,j)进行处理，仅对高于拟合参考面的表面进行分割处理，若z(i,j)>h，则令z(i,j)＝z(i,j)-h，若z(i,j)<＝h，则令z(i,j)＝h；

(c3)使用Matlab对表面形貌离散点矩阵z(i,j)进行垂直方向和水平方向的Sobel算子滤波，然后求取模值，得到表面形貌离散点的梯度幅值矩阵gradz(i,j)，计算公式为：

式中，I_x(i,j)为使用Sobel算子对表面形貌离散点矩阵进行水平方向边缘滤波的卷积结果，I_y(i,j)为使用Sobel算子对表面形貌离散点矩阵进行垂直方向边缘滤波的卷积结果；

(c4)根据表面十点高度S_z对梯度幅值矩阵gradz(i,j)进行处理，若gradz(i,j)<0.05S_z，则令gradz(i,j)＝0；

(c5)使用Matlab中的watershed分水岭算法函数对梯度幅值矩阵进行分水岭分割，将表面分割为若干封闭区域。

5.根据权利要求1所述的三维粗糙表面微凸体拟合方法，其特征在于，d步骤中，包括以下步骤：

(d1)使用Matlab中的regiongrops图像处理函数对分水岭算法分割后的结果进行统计处理，得到每个封闭区域内所有离散点的坐标和区域边界上离散点的坐标；

(d2)针对单个封闭区域，若该区域内所有离散点高度的平均值大于该区域边界上离散点高度的平均值，且该区域内所有离散点高度的平均值大于拟合参考面高度h，则将该封闭区域内的离散点判定为组成一个微凸体。

6.根据权利要求1所述的三维粗糙表面微凸体拟合方法，其特征在于，e步骤中，包括以下步骤：

(e1)针对判定后单个微凸体的封闭区域，根据Matlab中的regiongrops图像处理函数可确定与该封闭区域具有相同标准二阶中心距的椭圆，令a为与该区域具有相同标准二阶中心矩的椭圆的长轴长度，b为与该区域具有相同标准二阶中心矩的椭圆的短轴长度，O点为该区域的质心；

(e2)针对(e1)中确定的椭圆区域，以O点为坐标系原点，椭圆长轴方向为x轴，椭圆短轴方向为y轴，高度方向为z轴建立坐标系，设区域内共有n个离散点，离散点高度分别为z(x₁,y₁)，…，z(x_i,y_i)，…，z(x_n,y_n)；

(e3)用椭球面对该椭圆区域内离散点进行拟合，令椭球面方程为：

式中，c为椭球面系数；

(e4)运用最小均方根误差来判定最适合的拟合椭球面，拟合椭球面最小均方误差Er由以下公式确定：

令即可计算得到每个椭球面的系数c：

将得到的椭球面系数c代入椭球面方程，即可得到微凸体离散点的拟合椭球面方程，确定拟合形貌。

7.根据权利要求6所述的三维粗糙表面微凸体拟合方法，其特征在于：e步骤中，包括以下步骤：

(e5)每一个椭球顶端曲率半径即为微凸体顶峰曲率半径，椭球在x轴方向等效曲率半径为R_x，在y轴方向等效曲率半径为R_y，计算公式如下：

设共有m个微凸体，则微凸体平均曲率半径R为：

(e6)微凸体峰点高度标准偏差σ、微凸体密度η为：

η＝m/A (11)

式中，A为整个粗糙表面的名义面积。