CN113029049A - 一种基于加权正负余量方差最小化算法的复杂工件光学测量方法 - Google Patents

一种基于加权正负余量方差最小化算法的复杂工件光学测量方法 Download PDF

Info

Publication number
CN113029049A
CN113029049A CN202110573411.XA CN202110573411A CN113029049A CN 113029049 A CN113029049 A CN 113029049A CN 202110573411 A CN202110573411 A CN 202110573411A CN 113029049 A CN113029049 A CN 113029049A
Authority
CN
China
Prior art keywords
negative
formula
point
positive
point cloud
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110573411.XA
Other languages
English (en)
Other versions
CN113029049B (zh
Inventor
朱大虎
吕睿
王周君
华林
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University of Technology WUT
Original Assignee
Wuhan University of Technology WUT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University of Technology WUT filed Critical Wuhan University of Technology WUT
Priority to CN202110573411.XA priority Critical patent/CN113029049B/zh
Publication of CN113029049A publication Critical patent/CN113029049A/zh
Application granted granted Critical
Publication of CN113029049B publication Critical patent/CN113029049B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00—Measuring arrangements characterised by the use of optical techniques
    • G01B11/30—Measuring arrangements characterised by the use of optical techniques for measuring roughness or irregularity of surfaces
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Abstract

本发明公开了一种基于加权正负余量方差最小化算法的复杂工件光学测量方法，首先由光学扫描仪获得工件测量点云，并由三维软件获得带有法向量的CAD点云；之后由Kdtree搜索算法获得最近点对集合，并将测量点云划分为正负余量测点；定义正负余量测点的权重系数函数后，定义加权正负偏差距离，由加权正负偏差距离生成新的目标函数；对获得的对应点对集合应用目标函数求解转换矩阵，并将转换矩阵应用于测量点云转换；将转换后的测量点云与标准CAD点云比对并生成色谱图。本发明有效解决传统算法在测量点云存在负余量和异常余量(离群点、离群噪声等)情况下的匹配失真问题，适用于余量分布复杂且未知的工件测量，具有精度高、效率高及测量稳定等优点。

Description

一种基于加权正负余量方差最小化算法的复杂工件光学测量 方法
技术领域
本发明属于复杂工件加工制造领域，涉及一种点云匹配算法及光学测量技术，具体涉及一种基于加权正负余量方差最小化算法（WPMAVM算法）的复杂工件光学测量方法。
背景技术
在复杂工件加工制造领域，为了保证工件的加工质量，需要对加工完成的工件进行光学测量并与标准CAD模型比对，以判断工件是否满足加工标准。以热端冲压工艺为例，在高温环境下，单次热冲压加工出的批量钢板并不能全部满足加工精度要求，当冲压强度过大时会导致钢板凹陷，而冲压强度过小时会导致钢板凸出，应对冲压完成后的钢板测量，对不满足加工精度的钢板冲压补偿，使其满足加工精度要求，因此热端冲压件存在分布未知的负余量，当冲压效果过差时会出现异常余量(错误的冲压方式使得工件局部区域严重偏离CAD模型标准件)。以复杂叶片为例，叶片通常具有前后缘薄、型面扭曲、易弯曲变形、材料难加工等特点，在前后缘等高曲率区域磨抛时材料去除量大易过磨，而在低曲率区域磨抛时材料去除量小易欠磨，单次叶片加工通常不能保证加工精度，尤其是叶片加工表面余量分布不均时，需要通过叶片表面测量工序获取叶片在当前加工工序完成时叶片的加工表面余量，以此判断叶片加工是否合格，而加工叶片在过磨区域会呈负余量，欠磨区域呈异常余量。上述复杂工件均属于余量分布复杂且未知的工件，传统匹配算法易受负余量和异常余量影响而匹配失真，难以获得工件的真实加工情况。文献“A method for registrationof 3-D shapes”(IEEE Transactions on Pattern Analysis and MachineIntelligence, 14 (1992) 239 - 256.)提出了一种ICP算法用于点云匹配，该算法以点到点的距离平方和作为目标函数，用于上述复杂工件测量时，未考虑余量分布和异常余量的存在，会导致匹配失真。申请号为CN201510226138.8的授权发明专利提出了一种基于距离方差最小的工件点云匹配算法，对余量已知的刚性工件，或余量已知但分布非均匀工件(如叶片凹面余量已知为d1mm，凸面已知为d2mm)可获得良好的匹配效果，但对于上述余量未知，存在负余量和异常余量的复杂工件，该方法未区分测点的正负性和异常性进而导致匹配失真，无法对加工后的复杂工件精确测量。本发明针对上述问题，提出了一种考虑测点正负性和异常性的WPMAVM算法，并应用于复杂工件光学测量，可有效解决传统算法因工件存在负余量测点和异常余量测点引起的匹配失真问题。
发明内容
本发明的目的在于提供一种基于加权正负余量方差最小化(Weighted Plus-and-Minus Allowance Variance Minimization，以下简称“WPMAVM”)算法的复杂工件光学测量方法，以解决现有算法在工件存在负余量、异常余量且余量分布复杂未知情况下的匹配失真问题，构造了一种新的加权正负偏差距离，并基于上述距离定义了新的目标函数，同时建立了算法的流程，并最终应用于复杂工件光学测量。
为了解决上述技术问题，本发明采用的技术方案为：
一种基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，包括如下步骤：
步骤1、由光学扫描仪扫描工件，获取工件的测量点云P，通过三维软件离散CAD模型获得点云并计算法向量，得到带有法向量的CAD模型点云Q n ，设置初始转换参数；
步骤2、由Kdtree搜索算法获得测量点云P在Q n 中的最近点及最近点处的法向量；根据测量点云和CAD模型点云的位置关系将测量点云P划分为正余量测点集合P 1 （p 1 ，p 2 ，p 3 .....p i ）和负余量测点集合P 2 （p 1 ，p 2 ，p 3 .....p l ），i和l均为正整数的顺序角标；
步骤3、定义正余量测点p i 和负余量测点p l 的权重系数函数，并计算单步转换后的加权正负均值；
步骤4、由步骤3定义的权重系数函数和加权正负均值定义加权正负偏差距离，由定义的加权正负偏差距离生成WPMAVM算法的目标函数；
步骤5、应用目标函数求解转换矩阵，并将转换矩阵应用于测量点云，重复步骤2-步骤5直到步骤1中转换参数满足收敛条件；
步骤6、定义当测量点云存在负余量点云和异常余量点云时的误差评价函数；
步骤7、将转换后的测量点云与标准CAD点云比对并生成色谱图，完成工件测量。
作为优选，步骤1中，初始转换参数包括：最大迭代次数N max ，初始迭代次数n=0，转换矩阵H，总体转换矩阵H F ，H F 的初始值设置为四阶单位转换矩阵。
作为优选，步骤2中，测量点云位于CAD模型上方的划分为正余量测点集合P 1 ，测量点云位于CAD模型下方的划分为负余量测点集合P 2 ，具体步骤如下：
正余量测点集合P 1 中的测点p i 和CAD点云中的最近点q j 及q j 处的法向量n j 满足公式(1)
Figure 100002_DEST_PATH_IMAGE001
公式(1)
负余量测点集合P 2 中的测点p l 和CAD点云中的最近点q j 及q j 处的法向量n j 满足公式(2)
Figure 100002_DEST_PATH_IMAGE002
公式(2)
作为优选，步骤3中，定义正余量测点和负余量测点的权重系数函数具体步骤如下：
正余量点云权重系数w i 的计算公式如下：
Figure 100002_DEST_PATH_IMAGE003
公式(3)
其中
Figure 100002_DEST_PATH_IMAGE004
为当前位置的正余量测点p i 到最近点q j 所在切平面之间的距 离；
Figure 100002_DEST_PATH_IMAGE005
为当前位置的正余量测点p i 到最近点q j 所在切平面加权距离正均 值；m为满足
Figure 100002_DEST_PATH_IMAGE006
的测点总数，
Figure 100002_DEST_PATH_IMAGE007
；
Figure 100002_DEST_PATH_IMAGE008
为自适应调整因 子，根据具体情况调整；
负余量点云权重w l 的计算公式如下：
Figure 100002_DEST_PATH_IMAGE009
公式(4)
其中
Figure 100002_DEST_PATH_IMAGE010
为当前位置的负余量测点p l 到最近点q j 所在切平面 之间的距离;
Figure 100002_DEST_PATH_IMAGE012
为当前位置的负余量测点p l 到最近点q j 所在切平面加权距 离负均值；n为满足
Figure 100002_DEST_PATH_IMAGE013
的测点总数，
Figure 100002_DEST_PATH_IMAGE014
；
Figure 100002_DEST_PATH_IMAGE015
为自适 应调整因子，根据具体情况调整。
作为优选，步骤3中，计算单步转换后的加权正负均值具体方法如下：
单步转换后的正余量测点p i 到相应最近点q j 所在切平面加权距离正均值如下：
Figure 100002_DEST_PATH_IMAGE016
公式(5)
其中d iTDM 表示单步转换后的正余量测点p i 到相应最近点q j 所在切平面距离如公式(6)。
Figure 100002_DEST_PATH_IMAGE017
公式(6)
其中p i+ 表示单步转换后的正余量测点，
Figure 100002_DEST_PATH_IMAGE018
表示单步转换的微 分旋转，δx、δy、δz分别表示测点（x，y，z）沿坐标系X、Y、Z轴的微分旋转运动量，
Figure 100002_DEST_PATH_IMAGE019
表示单步转换的微分平移，Δx、Δy、Δz分别表示测点（x，y，z） 沿坐标系X、Y、Z轴的微分平移运动量。
同理单步转换后的负余量测点p l 到相应最近点q j 所在切平面加权距离负均值如下：
Figure 100002_DEST_PATH_IMAGE020
公式(7)
其中d lTDM 表示单步转换后的负余量测点p l 到相应最近点q j 所在切平面距离，计算公式如下
Figure 100002_DEST_PATH_IMAGE021
公式(8)
其中p l+ 表示单步转换后的负余量测点。
作为优选，步骤4中，定义加权正负偏差距离并由定义的加权正负偏差距离生成WPMAVM算法目标函数的具体步骤如下：
当正余量测点p i 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure 100002_DEST_PATH_IMAGE022
时，定义加权正偏差距离
Figure 100002_DEST_PATH_IMAGE023
如下式：
Figure 100002_DEST_PATH_IMAGE024
公式(9)
当负余量测点p l 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure 100002_DEST_PATH_IMAGE025
时，定义加权负偏差距离
Figure 100002_DEST_PATH_IMAGE026
如下式：
Figure 100002_DEST_PATH_IMAGE027
公式(10)
由上述定义的加权正偏差距离和加权负偏差距离可建立目标函数如下：
Figure 100002_DEST_PATH_IMAGE028
公式(11)
作为优选，步骤5中，应用目标函数求解转换矩阵的具体步骤为：
由单步转换的微分旋转
Figure 100002_DEST_PATH_IMAGE029
和微分平移
Figure 100002_DEST_PATH_IMAGE030
可构成一组运动旋量
Figure 100002_DEST_PATH_IMAGE031
，δx、δy、δz分别表示测点（x，y，z）沿坐标系X、Y、Z轴的微分旋转运动量，
Figure 573956DEST_PATH_IMAGE030
表示单步转换的微分平移，Δx、Δy、Δz分别表示测点（x，y，z）沿坐 标系X、Y、Z轴的微分平移运动量；
则单步转换的微分运动可表示为:
Figure 100002_DEST_PATH_IMAGE032
公式(12)
首先讨论当正余量测点p i 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure 100002_DEST_PATH_IMAGE033
时，对加权正偏差距离平方和化简可得：
Figure 100002_DEST_PATH_IMAGE034
公式(13)
对上述公式(6)化简可得：
Figure 100002_DEST_PATH_IMAGE035
公式(14)
其中A i 为1x6矩阵，ξ为运动旋量6x1矩阵。
将
Figure 100002_DEST_PATH_IMAGE036
代入
Figure 100002_DEST_PATH_IMAGE037
得：
Figure 100002_DEST_PATH_IMAGE038
公式(15)
其中D 1为标量，E 1为6x6矩阵，F 1为6x1矩阵，
Figure 100002_DEST_PATH_IMAGE039
为加权正平均线矢量，计算公式如 下：
Figure 100002_DEST_PATH_IMAGE040
公式(16)
当负余量测点p l 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure 100002_DEST_PATH_IMAGE041
时，可得化简后的负余量测点p l 到最近点q j 所在切平面距离d lTDM 如 式 ：
Figure 100002_DEST_PATH_IMAGE042
公式(17)
同理将
Figure 100002_DEST_PATH_IMAGE043
代入
Figure 100002_DEST_PATH_IMAGE044
得：
Figure 100002_DEST_PATH_IMAGE045
公式(18)
故WPMAVM算法目标函数可化简为：
Figure 100002_DEST_PATH_IMAGE046
公式(19)
由公式（19）对转换矢量ξ求导可得其最小化条件为：
Figure 100002_DEST_PATH_IMAGE047
公式(20)
求解线性方程组可得：
Figure 100002_DEST_PATH_IMAGE048
公式(21)
经求解，转换矢量ξ的收敛结果为：
Figure 100002_DEST_PATH_IMAGE049
公式(22)
其中
Figure 100002_DEST_PATH_IMAGE050
为加权负平均线矢量如下：
Figure 100002_DEST_PATH_IMAGE051
公式(23)
由转换矢量ξ可得旋转矢量δ和平移矢量t，故旋转矩阵R和平移矩阵T可由下式求得：
Figure 100002_DEST_PATH_IMAGE052
公式(24)
公式（24）中，e为自然底数，
Figure 100002_DEST_PATH_IMAGE053
表示对旋转矢量δ求反对称矩阵；由上述旋转矩阵R 和平移矩阵T可构造单步转换矩阵H从而完成目标函数求解。
作为优选，步骤5中，将转换矩阵应用于测量点云，重复步骤2-步骤5直到步骤1中转换参数满足收敛条件的具体步骤为：
每次采用单步转换矩阵H更新测量点云P和总体转换矩阵H F 后重复步骤2-步骤5计算新的单步转换矩阵H，再次更新测量点云P和总体转换矩阵H F ，直至更新迭代次数n大于最大迭代次数，更新迭代公式如下：
Figure 100002_DEST_PATH_IMAGE054
公式(25)
当迭代次数n满足n>N max 时，迭代终止。
作为优选，步骤6中，定义当测量点云存在负余量点云和异常余量点云时的误差评价函数具体步骤如下：
当存在异常点云和负余量点云时，匹配的理想结果是除异常点云外，剩余正余量点云和负余量点云均匀分布于CAD模型表面互不失真，因此应削弱异常点云对误差评价的影响，并对正负余量加以区分，定义如下误差评价函数：
Figure DEST_PATH_IMAGE055
公式(26)
作为优选，步骤7中，将转换后的测量点云与标准CAD点云比对并生成色谱图，完成工件测量的具体步骤为：
将转换后的测量点云与CAD模型比对后生成色谱图，通过色谱图可得到工件各区域的余量分布情况，还可以截取部分点云显示型面色谱图，完成工件测量。
与现有的用于复杂工件光学测量的匹配算法作比较，本发明的优势在于：
本发明提出的方法充分考虑了测点的正负性问题，并在目标函数中做了具体区分，充分考虑了传统算法目标函数在存在异常测点时倾向于异常测点平方和最小化而引起的匹配失真问题，在目标函数中添加了权重系数，弱化异常测点对匹配的影响，相比较于传统匹配算法，本发明在复杂工件光学测量过程中具有更高的测量精度，可有效解决传统算法在测量点云存在负余量和异常余量(离群点、离群噪声等)情况下的匹配失真问题，适用于余量分布复杂且未知的工件测量，具有精度高、效率高及测量稳定等优点。
本发明适用于热端冲压、机械加工后的复杂工件测量，对于余量分布规则的工件同样适用，且在余量分布复杂且未知的工件测量中相较于传统算法优势明显。
附图说明
图1测点单步转换示意图。
图2 ICP算法对加工后的叶片匹配倾斜示意图。
图3距离方差最小算法对加工后的叶片匹配倾斜示意图。
图4 WPMAVM算法对加工后的叶片匹配示意图。
图5 WPMAVM算法对加工后的叶片匹配测量图。
图6 三种算法对加工后的叶片匹配测量误差图。
图7 基于WPMAVM算法的复杂工件光学测量流程图。
图8 本发明实施例中工件测量色谱图。
附图标记：1-切平面，2-CAD模型曲面，3-测点，4-均值曲面，5-加权正均值曲面，6-两倍加权正均值曲面，7-加权负均值曲面，8-两倍加权负均值曲面。
具体实施方式
下面结合附图并以加工叶片测量为例对本发明的实施方式作进一步详细描述。以下实施例用于说明本发明，但不能用来限制本发明的范围。
以下结合附图对本发明的实施过程做相关解释，涉及部分内容而非全部内容。
S1、由光学扫描仪扫描工件，获取工件上测点3的坐标，由所有测点3的坐标组成测量点云P，每个测点3的坐标记为（x，y，z），通过三维软件离散CAD模型获得点云并计算法向量，得到带有法向量的CAD模型点云Q n ，设置初始转换参数。
S1.1、步骤S1中，初始转换参数包括：最大迭代次数N max =15，初始迭代次数n=0，转换矩阵H，总体转换矩阵H F ，H F 的初始值设置为四阶单位转换矩阵。
S2、由Kdtree搜索算法获得P在Q中的最近点及最近点处的法向量；根据测量点云和CAD模型点云的位置关系将测量点云P划分为正余量测点集合P 1 （p 1 ，p 2 ，p 3 .....p i ）和负余量测点集合P 2 （p 1 ，p 2 ，p 3 .....p l ），i和l均为正整数的顺序角标。
S2.1、步骤S2中，测量点云位于CAD模型上方的划分为正余量测点集合P 1 ，测量点云位于CAD模型下方的负余量测点集合P 2 ，具体步骤如下：
正余量测点集合P 1 中的测点p i 和CAD点云中的最近点q j 及q j 处的法向量n j 满足公式(1)
Figure 737434DEST_PATH_IMAGE001
公式 (1)
于负余量测点集合P 2 中的测点p l 和CAD点云中的最近点q j 及q j 处的法向量n j 满足公式(2)
Figure 517171DEST_PATH_IMAGE002
公式(2)
图1为测点单步转换示意图，测点分为正余量测点p i (位于CAD模型曲面2以上)和负余量测点p l (位于CAD模型曲面2以下)，相应的最近点为q j ，q j 在CAD模型曲面2上的切面为切平面1，p t 为p i+ 和p l+ 在曲面上的投影点，p i+ 和p l+ 为F(R,T)转换后的更新点。图2为 ICP算法对加工后的叶片匹配倾斜示意图，图3为距离方差最小算法对加工后的叶片匹配倾斜示意图，所有测点3到CAD模型曲面2之间的距离求平均，得到均值距离，与CAD模型曲面2相距均值距离的曲面为均值曲面4。传统算法未区分正负测点，负余量测点会引起传统算法的匹配倾斜；未对异常测点加以限制，目标函数倾向于异常测点平方和最小化。综上所述，传统方法在存在负余量和异常余量的情况下会匹配失真，从而影响测量结果，最终影响加工工件的轮廓质量评估。
S3、定义正余量测点和负余量测点的权重系数函数，并计算单步转换后的加权正负均值。
S3.1、步骤S3中，定义正余量测点和负余量测点的权重系数函数具体步骤如下：
w i 为正余量点云权重系数函数如公式(3)，
Figure DEST_PATH_IMAGE057
为自适应调整因子，本发明在 加工叶片匹配实验中取2，可视具体情况调节。
Figure 286282DEST_PATH_IMAGE003
公式(3)
其中
Figure 162971DEST_PATH_IMAGE004
为当前位置的正余量测点p i 到最近点q j 所在切平面1的距离；
Figure 985433DEST_PATH_IMAGE005
为当前位置的正余量测点p i 到最近点q j 所在切平面1的加权距离正均 值，如图4所示，与CAD模型曲面2相距加权距离正均值的曲面为加权正均值曲面5，与CAD模 型曲面2相距加权距离正均值两倍的曲面为两倍加权正均值曲面6；m为满足
Figure 65516DEST_PATH_IMAGE006
的测点总数，
Figure 326733DEST_PATH_IMAGE007
。
w l 为负余量点云权重函数如公式(4)，
Figure DEST_PATH_IMAGE058
为自适应调整因子，本发明在 加工叶片匹配实验取2，可视具体情况调节。
Figure 392951DEST_PATH_IMAGE009
公式(4)
其中
Figure 199364DEST_PATH_IMAGE010
为当前位置的负余量测点p l 到最近点q j 所在切平面 1的距离；
Figure DEST_PATH_IMAGE059
为当前位置的负余量测点p l 到最近点q j 所在切平面1的加 权距离负均值，如图4所示，与CAD模型曲面2相距加权距离负均值的曲面为加权负均值曲面 7，与CAD模型曲面2相距加权距离负均值两倍的曲面为两倍加权负均值曲面8；n为满足
Figure 327595DEST_PATH_IMAGE013
的测点总数，
Figure 126923DEST_PATH_IMAGE014
。
如图4所示，为WPMAVM算法对加工后的叶片匹配示意图，由于对负余量测点和异常余量测点分别加以区分和限制，在匹配过程中不会因负余量测点和异常余量测点的存在而匹配失真。
S3.2、步骤3中，计算单步转换后的加权正负均值具体方法如下：
单步转换后的正余量测点到切平面加权距离正均值如下：
Figure 650309DEST_PATH_IMAGE016
公式(5)
其中d iTDM 表示单步转换后的正余量测点p i 到最近点q j 所在切平面1距离，计算如公式(6)。
Figure 893202DEST_PATH_IMAGE017
公式(6)
其中p i+ 表示单步转换后的正余量测点，
Figure 197145DEST_PATH_IMAGE018
表示单步转换的微 分旋转，δx、δy、δz分别表示测点（x，y，z）沿坐标系X、Y、Z轴的微分旋转运动量，
Figure 783853DEST_PATH_IMAGE019
表示单步转换的微分平移，Δx、Δy、Δz分别表示测点（x，y，z） 沿坐标系X、Y、Z轴的微分平移运动量。
同理单步转换后的负余量测点p l 到相应最近点q j 所在切平面1加权距离负均值如下：
Figure 489640DEST_PATH_IMAGE020
公式 (7)
其中d lTDM 表示单步转换后的负余量测点到切平面距离如公式(8)。
Figure 637856DEST_PATH_IMAGE021
公式(8)
其中p l+ 表示单步转换后的负余量测点。
S4、由步骤3定义的权重系数函数和加权正负均值定义加权正负偏差距离，由定义的加权正负偏差距离可生成WPMAVM算法的目标函数。
S4.1、步骤S4中，定义加权正负偏差距离并由定义的加权正负偏差距离生成WPMAVM算法目标函数的具体步骤如下：
当正余量测点p i 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure 429095DEST_PATH_IMAGE022
时，定义加权正偏差距离
Figure 819494DEST_PATH_IMAGE023
如下式：
Figure 379788DEST_PATH_IMAGE024
公式(9)
当负余量测点p l 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure 151435DEST_PATH_IMAGE025
时，定义加权负偏差距离
Figure 915123DEST_PATH_IMAGE026
如下式：
Figure 859945DEST_PATH_IMAGE027
公式(10)
由上述定义的加权正偏差距离和加权负偏差距离可建立目标函数如下：
Figure 524013DEST_PATH_IMAGE028
公式(11)
S5、应用目标函数求解转换矩阵，并将转换矩阵应用于测量点云和转换参数的更新，重复步骤S2-步骤S5直到步骤S1中转换参数满足收敛条件。
S5.1、步骤5中，应用目标函数求解转换矩阵的具体步骤为：
由单步转换的微分旋转
Figure 528878DEST_PATH_IMAGE029
和微分平移
Figure 779862DEST_PATH_IMAGE030
可构成一组运动旋量
Figure 262796DEST_PATH_IMAGE031
单步转换的微分运动可表示为:
Figure 735366DEST_PATH_IMAGE032
公式(12)
首先讨论当正余量测点p i 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure 160400DEST_PATH_IMAGE033
时，对加权正偏差距离平方和化简可得：
Figure 147947DEST_PATH_IMAGE034
公式(13)
对上述公式(6)化简可得：
Figure 185305DEST_PATH_IMAGE035
公式(14)
其中A i 为1x6矩阵，ξ为运动旋量6x1矩阵。
将
Figure 246802DEST_PATH_IMAGE036
代入
Figure 593469DEST_PATH_IMAGE037
得：
Figure 317581DEST_PATH_IMAGE038
公式(15)
其中D 1为标量，E 1为6x6矩阵，F 1为6x1矩阵，
Figure 407896DEST_PATH_IMAGE039
为加权正平均线矢量，计算公式如 下：
Figure 199266DEST_PATH_IMAGE040
公式(16)
当负余量测点p l 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure DEST_PATH_IMAGE060
时，可得化简后的负余量测点p l 到最近点q j 所在切平面1的距离d lTDM 如式 ：
Figure 966102DEST_PATH_IMAGE042
公式(17)
同理将
Figure 865925DEST_PATH_IMAGE043
代入
Figure 494353DEST_PATH_IMAGE044
得：
Figure 609070DEST_PATH_IMAGE045
公式(18)
故WPMAVM算法目标函数可化简为：
Figure 31962DEST_PATH_IMAGE046
公式(19)
由上式对转换矢量ξ求导可得其最小化条件为：
Figure 730665DEST_PATH_IMAGE047
公式(20)
求解线性方程组可得：
Figure 100466DEST_PATH_IMAGE048
公式(21)
经求解，转换矢量ξ的收敛结果为：
Figure 53379DEST_PATH_IMAGE049
公式(22)
其中
Figure 460221DEST_PATH_IMAGE050
为加权负平均线矢量如下：
Figure 334636DEST_PATH_IMAGE051
公式 (23)
由转换矢量ξ得旋转矢量δ和平移矢量t，故旋转矩阵R和平移矩阵T由下式求得：
Figure 819713DEST_PATH_IMAGE052
公式(24)
公式（24）中，e为自然底数，
Figure 830394DEST_PATH_IMAGE053
表示对旋转矢量δ求反对称矩阵。
由上述旋转矩阵R和平移矩阵T可构造单步转换矩阵H从而完成目标函数求解。
S5.2、步骤S5中，将转换矩阵应用于加工叶片测量点云，重复步骤S2-步骤S5直到步骤S1中转换参数满足收敛条件的具体步骤为：
每次采用单步转换矩阵H更新测量点云P和总体转换矩阵H F （H F 初始设置为四阶单位转换矩阵，更新一次后就是四阶转换矩阵）后重复步骤S2-步骤S5计算新的单步转换矩阵H，再次更新测量点云P和总体转换矩阵H F ，直至更新迭代次数n大于最大迭代次数，更新迭代公式如下：
Figure 860667DEST_PATH_IMAGE054
公式(25)
当迭代次数n满足n>N max 时，迭代终止，输出更新后的测量点云P并于工件点云比对，N max 根据需要设定，一般可以为10-50。图5为WPMAVM算法对加工后的叶片匹配测量图，黑色点云表示测点云，白色点云表示CAD点云，图中白色区域表示测点云位于CAD模型以下过磨为负余量测点，黑色区域表示测点云位于CAD模型以上为正余量测点。
S6、定义当测量点云存在负余量点云和异常余量点云时的误差评价函数。
S6.1、步骤6中，定义当测量点云存在负余量点云和异常余量点云时的误差评价函数具体步骤如下：
当存在异常点云和负余量点云时，匹配的理想结果是除异常点云外，剩余正余量点云和负余量点云均匀分布于CAD模型表面互不失真，因此应削弱异常点云对误差评价的影响，并对正负余量加以区分，定义如下误差评价函数：
Figure 769848DEST_PATH_IMAGE055
公式(26)
S7、步骤7中，如图8所示，将转换后的测量点云与标准CAD点云比对并生成色谱图，完成工件测量。
S7.1、将转换后的测量点云与标准CAD点云比对并生成色谱图，完成工件测量的具体步骤为：
将转换后的测量点云与CAD模型比对后生成色谱图，通过色谱图可得到工件各区域的余量分布情况，还可以截取部分点云显示型面色谱图，完成工件测量。
加工叶片根部存在未加工的区域，其余量分布是均匀的，可用于验证各算法匹配精度。对未加工区域截取部分测量点云分析误差色谱图，ICP算法出现了严重倾斜，匹配结束后叶片未加工区域的余量分布杂乱不均匀；距离方差最小算法出现了明显倾斜，匹配结束后叶片未加工区域的余量分布不均匀；而所提WPMAVM算法未出现明显倾斜，匹配结束后叶片未加工区域的余量分布均匀。应用公式(26)可计wRMSE误差，如图6为三种算法对加工后的叶片匹配测量误差图；WPMAVM算法的误差最小，距离方差最小算法的误差次之，ICP算法的误差最大，通过色谱图和误差计算的方式验证了所提算法相较于传统算法的优越性。
以上实施方式仅用于说明本发明，而非对本发明的限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行各种组合、修改或者等同替换，都不脱离本发明技术方案的精神和范围，均应涵盖在本发明的权利要求范围当中。

Claims (10)

1.一种基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，包括如下步骤：
步骤1、由光学扫描仪扫描工件，获取工件的测量点云P，通过三维软件离散CAD模型获得点云并计算法向量，得到带有法向量的CAD模型点云Q n ，设置初始转换参数；
步骤2、由Kdtree搜索算法获得测量点云P在Q n 中的最近点及最近点处的法向量；根据测量点云和CAD模型点云的位置关系将测量点云P划分为正余量测点集合P 1 （p 1 ，p 2 ，p 3 ..... p i ）和负余量测点集合P 2 （p 1 ，p 2 ，p 3 ..... p l ），i和l均为正整数的顺序角标；
步骤3、定义正余量测点p i 和负余量测点p l 的权重系数函数，并计算单步转换后的加权正负均值；
步骤4、由步骤3定义的权重系数函数和加权正负均值定义加权正负偏差距离，由定义的加权正负偏差距离生成WPMAVM算法的目标函数；
步骤5、应用目标函数求解转换矩阵，并将转换矩阵应用于测量点云和转换参数的更新，重复步骤2-步骤5直到步骤1中转换参数满足收敛条件；
步骤6、定义当测量点云存在负余量点云和异常余量点云时的误差评价函数；
步骤7、将转换后的测量点云与标准CAD点云比对并生成色谱图，完成工件测量。
2.根据权利要求1所述的基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，步骤1中，初始转换参数包括：最大迭代次数N max ，初始迭代次数n=0，转换矩阵H，总体转换矩阵H F ，H F 的初始值设置为四阶单位转换矩阵。
3. 根据权利要求2所述的基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，步骤2中，测量点云位于CAD模型上方的划分为正余量测点集合P 1 ，测量点云位于CAD模型下方的划分为负余量测点集合P 2 ，具体步骤如下：
正余量测点集合P 1 中的测点p i 和CAD点云中的最近点q j 及q j 处的法向量n j 满足公式(1)
Figure DEST_PATH_IMAGE001
公式(1)
负余量测点集合P 2 中的测点p l 和CAD点云中的最近点q j 及q j 处的法向量n j 满足公式(2)
Figure DEST_PATH_IMAGE002
公式(2)
上式中角标T表示转置矩阵。
4.根据权利要求3所述的基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，步骤3中，定义正余量测点和负余量测点的权重系数函数具体步骤如下：
正余量点云权重系数w i 的计算公式如下：
Figure DEST_PATH_IMAGE003
公式(3)
其中
Figure DEST_PATH_IMAGE004
为当前位置的正余量测点p i 到最近点q j 所在切平面之间的距离；
Figure DEST_PATH_IMAGE005
为当前位置的正余量测点p i 到最近点q j 所在切平面加权距离正均值；m 为满足
Figure DEST_PATH_IMAGE006
的测点总数，
Figure DEST_PATH_IMAGE007
；
Figure DEST_PATH_IMAGE008
为自适应调整因子，根据 具体情况调整；
负余量点云权重w l 的计算公式如下：
Figure DEST_PATH_IMAGE009
公式(4)
其中
Figure DEST_PATH_IMAGE010
为当前位置的负余量测点p l 到最近点q j 所在切平面之间 的距离；
Figure DEST_PATH_IMAGE011
为当前位置的负余量测点p l 到最近点q j 所在切平面加权距离 负均值；n为满足
Figure DEST_PATH_IMAGE012
的测点总数，
Figure DEST_PATH_IMAGE013
；
Figure DEST_PATH_IMAGE014
为自适应 调整因子，根据具体情况调整。
5.根据权利要求4所述的基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，步骤3中，计算单步转换后的加权正负均值具体方法如下：
单步转换后的正余量测点到最近点q j 所在切平面加权距离正均值如下：
Figure DEST_PATH_IMAGE015
公式(5)
其中d iTDM 表示单步转换后的正余量测点p i 到最近点q j 所在切平面距离，计算公式如下：
Figure DEST_PATH_IMAGE016
公式(6)
其中p i+ 表示单步转换后的正余量测点，
Figure DEST_PATH_IMAGE017
表示单步转换的微分旋 转，δx、δy、δz分别表示测点（x，y，z）沿坐标系X、Y、Z轴的微分旋转运动量，
Figure DEST_PATH_IMAGE018
表示单步转换的微分平移，Δx、Δy、Δz分别表示测点（x，y，z） 沿坐标系X、Y、Z轴的微分平移运动量；
同理单步转换后的负余量测点p l 到相应最近点q j 所在切平面加权距离负均值如下：
Figure DEST_PATH_IMAGE019
公式(7)
其中d lTDM 表示单步转换后的负余量测点p l 到最近点q j 所在切平面距离，计算公式如下：
Figure DEST_PATH_IMAGE020
公式(8)
其中p l+ 表示单步转换后的负余量测点。
6.根据权利要求5所述的基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，步骤4中，定义加权正负偏差距离并由定义的加权正负偏差距离生成WPMAVM算法目标函数的具体步骤如下：
当正余量测点p i 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure DEST_PATH_IMAGE021
时，定义加权正偏差距离
Figure DEST_PATH_IMAGE022
如下式：
Figure DEST_PATH_IMAGE023
公式(9)
当负余量测点p l 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure DEST_PATH_IMAGE024
时，定义加权负偏差距离
Figure DEST_PATH_IMAGE025
如下式：
Figure DEST_PATH_IMAGE026
公式(10)
由上述定义的加权正偏差距离和加权负偏差距离建立目标函数如下：
Figure DEST_PATH_IMAGE027
公式(11)。
7.根据权利要求6所述的基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，步骤5中，应用目标函数求解转换矩阵的具体步骤为：
由单步转换的微分旋转
Figure DEST_PATH_IMAGE028
和微分平移
Figure DEST_PATH_IMAGE029
构成 一组运动旋量
Figure DEST_PATH_IMAGE030
，则单步转换的微分运动表示 为:
Figure DEST_PATH_IMAGE031
公式(12)
当正余量测点p i 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure DEST_PATH_IMAGE032
时，对加权正偏差距离平方和化简得：
Figure DEST_PATH_IMAGE033
公式(13)
对上述公式(6)化简得：
Figure DEST_PATH_IMAGE034
公式(14)
其中A i 为1x6矩阵，ξ为运动旋量6x1矩阵；
将
Figure DEST_PATH_IMAGE035
代入
Figure DEST_PATH_IMAGE036
得：
Figure DEST_PATH_IMAGE037
公式(15)
其中D 1为标量，E 1为6x6矩阵，F 1为6x1矩阵，
Figure DEST_PATH_IMAGE038
为加权正平均线矢量，计算公式如下：
Figure DEST_PATH_IMAGE039
公式(16)
当负余量测点p l 与CAD模型点云中的最近点q j 及q j 处法向量n j 满足
Figure DEST_PATH_IMAGE040
时，得化简后的负余量测点p l 到最近点q j 所在切平面距离d lTDM 如 式 ：
Figure DEST_PATH_IMAGE041
公式(17)
同理将
Figure DEST_PATH_IMAGE042
代入
Figure DEST_PATH_IMAGE043
得：
Figure DEST_PATH_IMAGE044
公式(18)
故WPMAVM算法目标函数化简为：
Figure DEST_PATH_IMAGE045
公式(19)
由公式（19）对转换矢量ξ求导得其最小化条件为：
Figure DEST_PATH_IMAGE046
公式(20)
求解线性方程组得：
Figure DEST_PATH_IMAGE047
公式(21)
经求解，转换矢量ξ的收敛结果为：
Figure DEST_PATH_IMAGE048
公式(22)
其中
Figure DEST_PATH_IMAGE049
为加权负平均线矢量如下：
Figure DEST_PATH_IMAGE050
公式(23)
由转换矢量ξ得旋转矢量δ和平移矢量t，故旋转矩阵R和平移矩阵T由下式求得：
Figure DEST_PATH_IMAGE051
公式(24)
公式（24）中，e为自然底数，
Figure DEST_PATH_IMAGE052
表示对旋转矢量δ求反对称矩阵；
由上述旋转矩阵R和平移矩阵T构造单步转换矩阵H，从而完成目标函数求解。
8.根据权利要求7所述的基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，步骤5中，将转换矩阵应用于测量点云，重复步骤2-步骤5直到步骤1中转换参数满足收敛条件的具体步骤为：
每次采用单步转换矩阵H更新测量点云P和总体转换矩阵H F 后重复步骤2-步骤5计算新的单步转换矩阵H，再次更新测量点云P和总体转换矩阵H F ，直至更新迭代次数n大于最大迭代次数，更新迭代公式如下：
Figure DEST_PATH_IMAGE053
公式(25)
当迭代次数n满足n>N max 时，迭代终止。
9.根据权利要求8所述的基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，步骤6中，定义当测量点云存在负余量点云和异常余量点云时的误差评价函数具体步骤如下：
当存在异常点云和负余量点云时，匹配的理想结果是除异常点云外，剩余正余量点云和负余量点云均匀分布于CAD模型表面互不失真，因此应削弱异常点云对误差评价的影响，并对正负余量加以区分，定义如下误差评价函数：
Figure DEST_PATH_IMAGE054
公式(26)。
10.根据权利要求9所述的基于加权正负余量方差最小化算法的复杂工件光学测量方法，其特征在于，步骤7中，将转换后的测量点云与标准CAD点云比对并生成色谱图，完成工件测量的具体步骤为：
将转换后的测量点云与CAD模型比对后生成色谱图，通过色谱图得到工件各区域的余量分布情况，完成工件测量。
CN202110573411.XA 2021-05-25 2021-05-25 基于加权正负余量方差最小化算法的工件光学测量方法 Active CN113029049B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110573411.XA CN113029049B (zh) 2021-05-25 2021-05-25 基于加权正负余量方差最小化算法的工件光学测量方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110573411.XA CN113029049B (zh) 2021-05-25 2021-05-25 基于加权正负余量方差最小化算法的工件光学测量方法

Publications (2)

Publication Number Publication Date
CN113029049A true CN113029049A (zh) 2021-06-25
CN113029049B CN113029049B (zh) 2021-09-03

Family

ID=76456211

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110573411.XA Active CN113029049B (zh) 2021-05-25 2021-05-25 基于加权正负余量方差最小化算法的工件光学测量方法

Country Status (1)

Country Link
CN (1) CN113029049B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113674227A (zh) * 2021-08-02 2021-11-19 上海工程技术大学 一种用于离子推力器栅极组件的层间距检测方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH07182521A (ja) * 1993-10-26 1995-07-21 Gerber Syst Corp プリント基板の自動検査方法及び検査システム
US5815400A (en) * 1995-07-10 1998-09-29 Mitsubishi Denki Kabushiki Kaisha Machining method using numerical control apparatus
CN104867136A (zh) * 2015-05-06 2015-08-26 华中科技大学 一种基于距离方差最小的工件点云匹配算法
CN110990975A (zh) * 2019-12-11 2020-04-10 南京航空航天大学 基于实测数据的舱门边框轮廓铣切余量测算方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH07182521A (ja) * 1993-10-26 1995-07-21 Gerber Syst Corp プリント基板の自動検査方法及び検査システム
US5815400A (en) * 1995-07-10 1998-09-29 Mitsubishi Denki Kabushiki Kaisha Machining method using numerical control apparatus
CN104867136A (zh) * 2015-05-06 2015-08-26 华中科技大学 一种基于距离方差最小的工件点云匹配算法
CN110990975A (zh) * 2019-12-11 2020-04-10 南京航空航天大学 基于实测数据的舱门边框轮廓铣切余量测算方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
WEN-LONG LI 等: "3-D Shape Matching of a Blade Surface in robotic grinding applications", 《IEEE/ASME TRANSACTIONS ON MECHATRONICS》 *
李振华 等: "基于点云与CAD模型配准对齐的毛坯加工定位技术", 《制造业自动化》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113674227A (zh) * 2021-08-02 2021-11-19 上海工程技术大学 一种用于离子推力器栅极组件的层间距检测方法
CN113674227B (zh) * 2021-08-02 2023-08-08 上海工程技术大学 一种用于离子推力器栅极组件的层间距检测方法

Also Published As

Publication number Publication date
CN113029049B (zh) 2021-09-03

Similar Documents

Publication Publication Date Title
CN104867136B (zh) 一种基于距离方差最小的工件点云匹配算法
CN110069041B (zh) 一种基于在机测量的工件加工方法及系统
CN110688709B (zh) 一种基于工件点云模型的蒙皮工艺模型修正方法
CN109101741B (zh) 一种基于三角网格简化的复杂曲面检测自适应采样方法
CN112396690B (zh) 基于改进型向心参数化法的曲面高精重构方法
CN112033338B (zh) 一种叶片类曲面接触式扫描测量测头半径面补偿方法
CN111820545A (zh) 一种结合离线与在线扫描自动生成鞋底喷胶轨迹的方法
CN113029049B (zh) 基于加权正负余量方差最小化算法的工件光学测量方法
CN114972387B (zh) 基于三维实测的复材成型过程模具变形修配方法及系统
CN112817271A (zh) 基于在机测量的铸造机匣毛坯加工余量优化方法
CN113094964B (zh) 一种生成叶片加工坐标的方法和装置
CN110738728B (zh) 基于线性组合过渡算法的叶片修复模型重建方法
CN109033669B (zh) 基于万能运动参数驱动的螺旋锥齿轮仿真加工建模方法
CN116720268A (zh) 一种周期性描述的叶片叶型全局光顺重建方法
CN113970311A (zh) 一种航空发动机叶片矢量逼近迭代测量方法
KR100344917B1 (ko) 선박용 1차곡면 외판의 가공정보 산출방법
CN115056213A (zh) 一种面向大型复杂构件的机器人轨迹自适应修正方法
CN116718137A (zh) 面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法
CN113954102A (zh) 一种基于离线编程的百叶轮抛磨叶片路径规划方法
CN113393507A (zh) 无人机点云与地面三维激光扫描仪点云配准方法
CN116610905B (zh) 一种基于各向异性尺度修正的反距离权重数据插值方法
Makem et al. A virtual inspection technique for assessing the dimensional accuracy of forged compressor blades using FE modeling and CMM inspection
CN110385601B (zh) 一种汽车覆盖件模具研配型面的球形铣刀刀轨的生成方法
CN115994931A (zh) 基于去伪加权方差最小化算法的工件点云精配准方法
CN115235439B (zh) 基于网格划分和数据校正的地块测亩方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant