CN104867136A - 一种基于距离方差最小的工件点云匹配算法 - Google Patents

Abstract

本发明公开了一种基于距离方差最小的工件点云匹配算法，包括以下步骤：(1)分别构建测量模型点云和带法矢的设计模型点云：(2)寻找最近点点集及最近点法矢集合；(3)根据余量要求，建构基于距离方差的目标函数：(4)通过目标函数求解测量模型与设计模型间的刚体变换参数，并利用刚体变化参数更新测量点云，实现点云匹配。本发明中的基于距离方差最小的点云匹配算法，对于工件存在余量的情况，不论凹面和凸面的余量是相等，还是呈差值或者比值关系，均能较好的得到点云匹配的结果，迭代效率高、匹配准确性好。

Description

一种基于距离方差最小的工件点云匹配算法

技术领域

本发明属于机器人砂带磨抛技术领域，更具体地，涉及一种基于距离方差最小的点云匹配算法，该算法可用于获取工件坐标系和校正工件装夹误差，尤其可针对叶片等其凹、凸面磨削余量不均匀的情况。

背景技术

叶片类复杂曲面零件在航空发动机、汽轮机等能源动力设备中被广泛应用。其型面加工质量将对动力设备的性能有着决定性的影响。目前叶片的磨抛工序基本依赖于手工方法，手工磨抛存在效率低、环境恶劣、一致性低等缺点，从未来发展趋势看，手工方法将会逐步被高柔性的机器人砂带磨抛取代。

在机器人砂带磨抛系统中，叶片类复杂曲面工件存在装夹误差、且没有特定的加工基准面。因此需校正装夹误差，确定基准面和加工余量。这些问题都依赖于叶片毛坯测量点云与设计模型的匹配，测量点云匹配效果对砂带磨抛质量有着关键影响。点云匹配一般通过求解设计模型和测量模型之间的最优刚体变化参数来改变二者的位姿使其整体距离最近。文献“Method for registration of 3-D shapes”(Robotics-DL tentative.International Society for Optics and Photonics,1992:586-606.)提出基于点-点距离的ICP算法。文献“Fitting B-spline curves to pointclouds by curvature-based squared distance minimization”.(ACMTransactions on Graphics,2006,25(2):214-238.)提出了基于点-切线距离的TDM算法。文献“复杂曲面零件数据拼合与精密加工技术研究”提出了基于自适应距离的ADF算法。

以上三篇文献中的匹配算法的目标函数都是基于距离平方和最小来求解刚体变化参数。设P＝{p₁,p₂,...p_i,...p_n}、Q＝{q₁,q₂,...q_i,...q_l}分别表示测量点云和设计模型的离散点云。基于距离平方和最小目标函数定义如下：

$\begin{array}{cc}\min & F=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{d}_i^2(t,R,p_i,q_i)\\ s.t.& q_i=q\left|\min \right||p_i-q||\end{array},$

其中t、R分别表示刚体变化参数中的平移变换矢量和旋转变换矩阵。q_i为p_i在设计模型Q中的欧式距离(对应欧几里德空间)最近点。根据距离定义的不同，衍生了不同的匹配算法：ICP算法采用点-点距离描述点-曲面最近点的几何度量，收敛速度取决于初始值。而TDM算法采用点-切面距离，收敛速度快但易陷入错误的局部最优值，ADF算法定义自适应距离，通过一个可变的修正系数的自适应的调整法向距离和切向距离的比重，以兼顾收敛的速度与稳定性。

以上ICP等基于距离最小化等算法只适合工件没有余量的情况，而当叶片的凹、凸面存在模型余量且相同时，若凹、凸面测量点云数量存在明显差异时或者叶片扭曲变形很大，则ICP、TDM等匹配算法将会使测量点云位姿向点云数量多的一面倾斜或凸面倾斜，因此无法保证测量点云均匀包裹设计模型，最终导致磨削时一面过切，一面欠切。

当叶片凹、凸面磨削余量不同时，匹配结果要求凹、凸面测量点云到设计模型之间的距离分别与凹、凸面设定的磨削余量相同。而ICP等匹配算法更无法达到预期的效果，其只适用于工件没有余量时点云匹配的情况，如工件的终检时的点云匹配。因此,针对机器人砂带磨抛中工件有余量的情况，现有的点云匹配算法并不能满足实际应用需要。

发明内容

针对现有技术的以上缺陷或改进需求，本发明的目的在于一种基于距离方差最小的工件点云匹配算法，其中通过对其关键目标函数的构建、目标函数的求解过程、以及求解过程中参量变换规则等进行改进，与现有技术相比能够有效解决有余量的测量模型与无余量的设计模型之间的匹配问题，并且该方法也适用于测量模型凹、凸面余量不相同情况下的模型匹配；此外，本发明中根据余量分配要求构建的误差判断函数，给出了匹配效果的评价指标，技术人员可根据该匹配效果评价指标，对该匹配方法的匹配效果得到精准认识。

为实现上述目的，按照本发明，提供了一种基于距离方差最小的工件点云匹配算法，其特征在于，包括以下步骤：

(1)分别构建测量模型点云和设计模型点云：

(1-1)扫描被加工的工件，获得测量模型点云P＝{p₁,p₂,...p_i,...p_n}，该测量模型点云含有n个点；

(1-2)根据工件设计模型，获得设计模型点云Q＝{q₁,q₂,...q_i,...q_l}及其法矢集合N＝{n₁,n₂,...n_i,...n_l}；该设计模型点云含有l个点，l>n；该法矢集合N中的第i个元素n_i为所述设计模型点云中的第i个元素q_i对应的单位法向矢量，i＝1,2,…,l；

(1-3)根据所述工件设计模型，依次对所述设计模型点云中的点进行划分；其中当q_i属于所述工件设计模型的凹面时，q_i为凹面点；当q_i属于所述工件设计模型的凸面时，q_i为凸面点；i＝1,2,…,l；

(2)依次根据所述测量模型点云P中的点p_i，i＝1,2,…,n，在所述设计模型点云Q中提取与p_i距离最近的设计模型点云中的点q′_i，q′_i∈Q，以及与q′_i对应的单位法向矢量n′_i，n′_i∈N，形成最近点点集Q′＝{q′₁,q′₂,...q′_i,...q′_n}和最近点法矢集合N′＝{n′₁,n′₂,...n′_i,...n′_n}；该最近点点集中的点与该最近点法矢集合中的点、以及所述测量模型点云中的点一一对应；

(3)根据余量要求，建构基于距离方差的目标函数：

(3-1)当凹面余量c_a与凸面余量c_t均为常数值时，目标函数为

$F_1=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(d_i-c_i)}^2,$

其中d_i＝||R·p_i+t-q′_i||，R、t分别为3×3的矩阵和3×1的矩阵；当q′_i属于设计模型凹面时，c_i＝c_a；当q′_i属于设计模型凸面时，c_i＝c_t；

然后执行步骤(4)；

(3-2)当凹面余量c_a与凸面余量c_t的比值为常数值c时，目标函数为

$F_2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{k}_i^2{d}_i^2-\mathrm{\λ}\frac{{\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{k}_i{d}_i\right)}^2}{n},$

其中λ∈(0,1]，d_i＝||R·p_i+t-q′_i||，R、t分别为3×3的矩阵和3×1的矩阵；当q′_i属于设计模型凹面时，k_i＝1；当q′_i属于设计模型凸面时，k_i＝c；

然后执行步骤(4)；

(4)根据所述测量模型点云P、最近点点集Q′和最近点法矢集合N′，求解步骤(3)中的R和t，使所述目标函数的数值最小，实现点云匹配。

作为本发明的进一步优选，所述步骤(4)包括以下步骤：

(4-1)当凹面余量c_a与凸面余量c_t均为常数值时，记

ξ＝(A^TA)^-1B，

其中， $A=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left[n_i^{\′T}\right|{(p_i\×n_i^{\′})}^T],$

$B=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left[\right[{(p_i-q_i^{\′\′})}^T\·n_i^{\′}\left]{\left[n_i^{\′T}\right|{(p_i\×n_i^{\′})}^T]}^T\right],$ 其中q″_i＝q′_i+c_in′_i，当q′_i属于设计模型凹面时，c_i＝c_a；当q′_i属于设计模型凸面时，c_i＝c_t；

接着，执行步骤(4-3)；

(4-2)当凹面余量c_a与凸面余量c_t的比值为常数值c时，记 $\mathrm{\ξ}=-{[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^T{E}_i-\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^T\right)\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i\right)/n]}^{-1}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F}_i{E}_i^T-\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F}_i\right)\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^T\right)/n]$ 其中，F_i＝k_i(p_i-q′_i)^T·n′_i，当q′_i属于设计模型凹面时，k_i＝1；当q′_i属于设计模型凸面时，k_i＝c；

接着，执行步骤(4-3)；

(4-3)记矩阵v和矩阵ω，v矩阵和ω矩阵的行数相同，列数也相同，并且矩阵v和矩阵ω还满足 $\mathrm{\ξ}=\left[\begin{array}{c}v\\ \mathrm{\ω}\end{array}\right],$

则 $t,v,R=e^{\widehat{\mathrm{\ω}}},$

其中，为ω的反对称矩阵；

(4-4)更新测量模型点云P：

记p′_i＝R·p_i+t，

并使p_i＝p′_i；i＝1,2,…,n；

更新后的测量模型点云P即为匹配后的工件测量点云。

作为本发明的进一步优选，所述基于距离方差最小的工件点云匹配算法，其特征在于，还包括以下步骤：

(5)计算误差判断函数：

(5-1)依次根据所述更新后的测量模型点云P中的点p_i，i＝1,2,…,n，在所述设计模型点云Q中提取与p_i距离最近的设计模型点云中的点q′_i，q′_i∈Q，以及与q′_i对应的单位法向矢量n′_i，n′_i∈N，形成最近点点集Q′＝{q′₁,q′₂,...q′_i,...q′_n}和最近点法矢集合N′＝{n′₁,n′₂,...n′_i,...n′_n}；该最近点点集中的点与该最近点法矢集合中的点、以及所述测量模型点云中的点一一对应；

(5-2)当凹面余量c_a与凸面余量c_t均为常数值时，均方误差Std为：

$\mathrm{Std}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(h_i-c_i)}^2}{n}},$

其中h_i＝(q′_i-p_i)^Tn′_i，当q′_i属于设计模型凹面时，c_i＝c_a；当q′_i属于设计模型凸面时，c_i＝c_t；

其中当Std大于预先设定的均方误差值，则执行步骤(5-4)；否则，测量模型点云P即为匹配后的工件点云，该基于距离方差最小的工件点云匹配算法结束；

(5-3)当凹面余量c_a与凸面余量c_t的比值为常数值c时，{q′₁,q′₂,...q′_i,...q′_n}中有m个q′_i属于设计模型凹面,有k个q′_i属于设计模型凸面，

则，均方误差Std、凹凸面余量比例平均误差p分别为：

$\mathrm{Std}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{k}_i^2{h}_i^2-\mathrm{\λ}{\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{k}_i{h}_i\right)}^2/n}{n}},$

$p=|k\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{d}_{\mathrm{ia}}/m\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{d}_{\mathrm{it}}-c|$

其中，h_i＝(q′_i-p_i)^Tn′_i；当q′_i属于设计模型凹面时，k_i＝1；当q′_i属于设计模型凸面时，k_i＝c；

当q′_i属于设计模型凹面时，d_ia＝||R·p_i+t-q′_i||；当q′_i属于设计模型凸面时，d_it＝||R·p_i+t-q′_i||；

其中当Std大于预先设定的均方误差值，或者p大于预先设定的比例误差值，则执行步骤(5-4)；否则，测量模型点云P即为匹配后的工件点云，该基于距离方差最小的工件点云匹配算法结束；

(5-4)重复步骤(3)到步骤(5)。

作为本发明的进一步优选，所述步骤(1-2)中，测量模型点云中的点的总数n和设计模型点云中的点的总数l满足l≥5×n。

作为本发明的进一步优选，所述步骤(3-2)中，λ＝1。

通过本发明所构思的以上技术方案，与现有技术相比，具有以下有益效果：

1.本发明中的点云匹配算法是基于距离方差最小，尤其是对于工件存在余量的情况，不论凹面和凸面的余量是相等，还是差值或者比值关系，均能较好的得到点云匹配的结果。

本发明的基于距离方差最小的点云匹配算法，即使凹面与凸面测量点云数量存在较大差异，或者叶片扭曲变形较大，测量点云位姿也不会倾斜，确保测量点云均匀包裹设计模型。并且，本发明中基于方差最小的点云匹配算法，针对叶片凹面、凸面磨削余量不同的情况，也能使测量点云与设计模型之间的匹配结果与凹面和凸面设定的磨削余量保持一致，能够广泛应用于叶片加工过程中时时校准、检测余量，对整体加工过程进行全面监测。

本发明中的点云匹配算法在迭代计算过程中，能够很好的收敛，求出符合误差的点云匹配关系，与目前常用的ICP等算法相比，即使对工件凹、凸面的余量要求不同(包括差值情况和比值情况)，也能快速得到匹配结果；并且，由于本方法中的目标函数是基于方差的，评价最终匹配效果的均方差参量可以控制得更好。

2.本发明根据步骤(3)中基于距离方差的目标函数，在对目标函数最优化过程中，采用一系列具体的公式对刚体运动学中的旋转变换矩阵R和平移变换矢量t进行运算，运算的效率高、并且结果的匹配性好。

在运用刚体运动学原理进行点云匹配过程中，刚体变换参数的求解是十分关键的，直接影响着迭代计算的效率以及输出结论的准确性。本发明根据构建的基于距离方差的目标函数，通过控制迭代过程中的旋转变换矩阵R和平移变换矢量t，并利用R和t更新测量点云，在确保匹配结果准确性的基础上，能够保证迭代计算的效率，快速实现点云匹配。

3.本发明根据余量分配要求，提供了合理评价匹配效果的指标。

本发明的误差判断函数并非是简单的目标函数的重复，本发明构建的误差判断函数对应的是测量点云(更新后的测量点云)与设计模型点云的距离在设计模型法线方向上的投影，与目标函数完全对应距离相比，采用投影数值来判断误差使误差判断更为精细，评价匹配效果也更为合理。

4.本发明还定义了方差权重系数λ，能够根据需要灵活控制方差对目标函数的影响比重。当凹面余量与凸面余量的比值为常数值时，λ∈(0,1](λ不能为0；若λ＝0时，则目标函数对应距离的平方和)。本发明优选λ＝1，此时的目标函数即为距离的方差(该方差为标准形式的方差)。目标函数中的λ与对应误差判断函数中的λ数值保持一致，使最终的误差判断也能根据需要灵活控制方差所占比重。

综上，本发明的基于距离方差最小的点云匹配算法，对于机器人砂带磨抛中工件有余量的情况，能够很好的满足实际应用中的要求，使得有余量的测量模型与设计模型之间相互匹配。若针对叶片凹、凸面存在余量且余量的数值相同，则可保证测量点云均匀包裹设计模型，凹面、凸面均匹配；若叶片凹、凸面存在余量但余量的数值不同时，可使得凹面、凸面各自的测量点云到设计模型之间的距离与设定的余量相同。

附图说明

图1是本发明点云匹配算法的流程图；

图2是测量点云(右侧)和设计模型(左侧)的初始姿态；

图3a是凹、凸面余量差值为常数时的叶片截面示意图；

图3b是凹、凸面余量比例为常数时的叶片截面示意图；

图4是匹配后设计模型与测量点云(该测量点云已被更新)的空间位置；

图5a是在凹、凸面余量分布为常数情况下，STD、TDM、ICP、WSD等匹配算法的均方差随迭代次数的变化情况，其中的STD即对应本发明中基于距离方差的点云匹配算法；图中横坐标轴为迭代次数，纵坐标轴为均方差数值；

图5b是在凹、凸面余量比为常数的情况下，STD、TDM、ICP、WSD等匹配算法的均方差随迭代次数的变化情况，图中横坐标轴为迭代次数，纵坐标轴为凹面距离与凸面距离的比值。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例1

本实施例工件优选航空叶片，其中测量点云通过机器人磨抛系统中的扫描仪扫描待磨削的叶片毛坯获取，设计模型点云通过离散叶片的CAD设计模型获取。测量模型点云定义在机器人末端坐标系，设计模型点云定义在设计模型坐标系，点云匹配算法的最终目标是建立机器人末端坐标系与设计模型坐标系的对应关系。刚体运动学方法为常用的点云匹配方法，本发明优选采用旋转变换矩阵R和平移变换矢量t进行运算。

如图3(a)和3(b)所示，需通过匹配，将设计模型点云与测量模型点云之间的余量(即对应设计模型点云与测量模型点云相应点之间的距离)与实际要求磨削的余量保持相同。本实例中凹、凸面余量分配要求包含以下D1与D2两种情况：

D1：{p_a＝c_a＝0.2mm、p_t＝c_t＝0.5mm}

和

D2：{p_a/p_t＝c＝0.6}

两种情况，其中p_a为叶片凹面的实际要求磨削的余量，c_a为叶片凹面设计模型点云与测量模型点云之间的余量；p_t为叶片凸面的实际要求磨削的余量，c_t为叶片凸面设计模型点云与测量模型点云之间的余量。图1所示为对应D2的流程图，具体包括以下步骤：

S1：测量模型点云为P＝{p₁,p₂,...p_i,...p_n}，设计模型点云为Q＝{q₁,q₂,...q_i,...q_l}，测量模型点云点的总数约为设计模型点云点的总数的1/5，即n/l≈1/5；创建与设计模型点云各个点对应的单位法向矢量的集合{n₁,n₂,...n_i,...n_l}，该单位法向矢量的集合中的元素与Q中的点一一对应；然后，根据设计模型的凹面与凸面，对Q中的各个点进行区分。

S2：利用例如KD-tree方法(也可利用二叉树、八叉树等方法)，在Q中寻找与p_i(i＝1,2,…,n)距离最近的点q_j(q_j∈Q)，以及与点q_j对应的单位法向矢量n_j(n_j∈{n₁,n₂,...n_i,...n_l})、点q_j是属于凹面点还是凸面点(根据步骤S1中设计模型的凹面、凸面判断)；令q′_i＝q_j，n′_i＝n_j，则Q′＝{q′₁,q′₂,...q′_i,...q′_n}为最近点点集，该最近点点集中的点与测量模型点云P中的点一一对应；{n′₁,n′₂,...n′_i,...n′_n}为单位法向矢量集合，该法向矢量的集合与测量模型点云P中的点一一对应。

S3：根据余量分配要求D1和D2分别定义目标函数

(1)要求为D1时，记Q'中对应设计模型凹面的点的总数为m个，对应设计模型凸面的点的总数为k个，n＝k+m；目标函数为

$\begin{array}{cc}\min & F_1=\underset{i=1}{\overset{m+k}{\mathrm{\Σ}}}{(d_i-c_i)}^2\\ & =\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(d_{\mathrm{ia}}-c_a)}^2+\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{(d_{\mathrm{it}}-c_t)}^2\end{array}$   (公式I)

当q′_i属于设计模型凹面时，c_i＝c_a，当q′_i属于设计模型凸面时，c_i＝c_t。

其中d_i＝||R·p_i+t-q′_i||(i＝1,2,…,n)；d_ia(i＝1,2,…,m)和d_it(i＝1,2,…,k)分别对应凹面d_i(q′_i属于设计模型凹面)和凸面d_i(q′_i属于设计模型凸面)。

(2)要求为D2时，目标函数为

$F_2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{k}_i^2{d}_i^2-\mathrm{\λ}\frac{{\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{k}_i{d}_i\right)}^2}{n}$   (公式II)

其中，自定义方差权重系数λ＝1；d_i＝||R·p_i+t-q′_i||(i＝1,2,…,n)，q′_i为与p_i对应的最近点点集Q′中的点；当q′_i属于设计模型凹面时，k_i＝1；当q′_i属于设计模型凸面时，k_i＝c。

S4：根据目标函数获取微分变换矢量

微分平移矢量、微分旋转矢量分别为v、ω(v和ω均为3×1的列向量)，微分变换矢量为 $\mathrm{\ξ}=\left[\begin{array}{c}v\\ \mathrm{\ω}\end{array}\right]$ (ξ为6×1的列向量)，公式(I)与(II)中的距离采用点-切面距离表示。对应于D1、D2的ξ的具体求解情况如下：

(1)在要求D1的情况下，ξ可由下式获取：

ξ＝(A^TA)^-1B，

其中，(矩阵A为1×6的矩阵)，

$B=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\left[\right[{(p_i-q_i^{\′\′})}^T\·n_i^{\′}\left]{\left[n_i^{\′T}\right|{(p_i\×n_i^{\′})}^T]}^T\right],$ (矩阵B为6×1的矩阵)，其中q″_i＝q′_i+c_in′_i(i＝1,2,…,n)。

(2)在要求D2的情况下，ξ可由下式获取：

$\mathrm{\ξ}=-{[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^T{E}_i-\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^T\right)\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i\right)/n]}^{-1}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F}_i{E}_i^T-\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F}_i\right)\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^T\right)/n]$

其中，F_i＝k_i(p_i-q′_i)^T·n′_i(i＝1,2,…,n；F_i为常数；当q′_i属于设计模型凹面时，k_i＝1；当q′_i属于设计模型凸面时，k_i＝c)； $E_i=k_i\left[n_i^{\′T}\right|{(p_i\×n_i^{\′})}^T](i=\mathrm{1,2},..,n;)$ E_i为1×6的矩阵，当q′_i属于设计模型凹面时，k_i＝1；当q_i′属于设计模型凸面时，k_i＝c)。

S5：根据ξ以及 $\mathrm{\ξ}=\left[\begin{array}{c}v\\ \mathrm{\ω}\end{array}\right],$ 则刚体变化参数g(t,R)由下式确定，

$t=v,R=e^{\widehat{\mathrm{\ω}}},$

其中，为ω的反对称矩阵(为3×3的矩阵，R也为3×3的矩阵)，

S6：更新测量模型点云P：

令p′_i＝R·p_i+t，

再使p_i＝p′_i。

S7：判断是否结束计算：当均方误差Std(有时还包括凹凸面余量比例平均误差p)满足设定的阈值条件时(如Std<0.05且p<5％)，迭代结束；否则执行步骤S3。均方误差Std和凹凸面余量比例平均误差p的计算过程如下：

(1)首先更新最近点点集Q′和最近点法矢集合N′：

依次根据上述更新后的测量模型点云P中的点p_i，i＝1,2,…,n，在所述设计模型点云Q中提取与p_i距离最近的设计模型点云中的点q′_i，q′_i∈Q，以及与q′_i对应的单位法向矢量n′_i，n′_i∈N，形成最近点点集Q′＝{q′_i,q′₂,...q′_i,...q′_n}和最近点法矢集合N′＝{n′₁,n′₂,...n′_i,...n′_n}；该最近点点集中的点与该最近点法矢集合中的点、以及所述测量模型点云中的点一一对应；

(2)D1时，迭代以所述匹配的均方误差Std控制在一定范围内时终止，Std可表示为

$\mathrm{Std}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(h_i-c_i)}^2}{n}},$

其中h_i＝(q′_i-p_i)^Tn′_i，表示更新后的测量模型点云中第i点到设计模型最近点的距离在法线方向n′_i上的投影值。

(3)D2时，迭代以所述匹配误差的均方差Std和凹凸面余量比例平均误差p控制在一定范围内时终止，分别可表示为

$\mathrm{Std}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_i^2{h}_i^2-\mathrm{\&lambda;}{\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_i{h}_i\right)}^2/n}{n}},$

$p=|k\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ia}}/m\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{d}_{\mathrm{it}}-c|.$

其中h_i＝(q′_i-p_i)^Tn′_i；由于Q'中有m个q′_i属于设计模型凹面,有k个q′_i属于设计模型凸面(n＝k+m)，当q′_i属于设计模型凹面时，d_ia＝d_i＝||R·p_i+t-q′_i||；当q′_i属于设计模型凸面时，

经本实施例计算，匹配结果参数见表1和附图3。图3中设计模型与更新后的测量点云姿态已经重合。由表1可知凹面余量为0.2999mm，凸面余量为0.499mm，基本等于理想值0.3mm，0.5mm。

表1.凹、凸面匹配结果

	平均距离(mm)	均方差(mm)	正向比例(mm)
				凹面	0.2999	0.0058	99.98％
凸面	0.4999	0.0078	99.98％
				整体	0.4053	0.0052	99.98％

匹配过程中，各匹配方法的均方差随匹配过程的变化情况如图5(a)所示，凹面、凸面距离比值p随匹配次数的变化情况如图5(b)所示。其中，WSD表示基于距离权重的匹配方法，STD为本发明基于方差权重的匹配方法，从图5(a)中可看出本发明方法与TDM方法都能快速收敛，但发明的误差在上述匹配方法中误差最小。从图5(b)中可看出本发明方法随匹配次数的增加，凹、凸面余量比值渐渐趋向于理想值0.6。

另外，当凹面、凸面的余量要求均匀时，即p_a/p_t＝1，可按实施例1中D2的情况处理，只是此时{p_a/p_t＝c＝1}，同样也能求解得到与设计模型最为匹配的测量模型结果；当没有余量要求时，可按实施例1中D1的情况处理，只是此时{p_a＝c_a＝0.2mm、p_t＝c_t＝0.5mm}，同样也能求解得到与设计模型最为匹配的测量模型结果。

本发明中各种矩阵、向量的运算均遵循常规线性代数运算法则，例如A|B表示将A、B矩阵简单合并，若A和B均为1×3的行向量，则合并后的矩阵为1×6；任何三维矢量p＝[p_x,p_y,p_z]^T，其反对称矩阵按常规方式建构，即 $\hat{p}=\left[\begin{array}{ccc}0& -p_z& p_y\\ p_z& 0& -p_x\\ -p_y& p_x& 0\end{array}\right].$ 本发明中的三维向量，可视为3×1的列向量；求和算符如也按其常规数学含义理解。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于距离方差最小的工件点云匹配算法，其特征在于，包括以下步骤：

(1)分别构建测量模型点云和设计模型点云：

(1-1)扫描被加工的工件，获得测量模型点云P＝{p₁,p₂,...p_i,...p_n}，该测量模型点云含有n个点；

(1-2)根据工件设计模型，获得设计模型点云Q＝{q₁,q₂,...q_i,...q_l}及其法矢集合N＝{n₁,n₂,...n_i,...n_l}；该设计模型点云含有l个点，l>n；该法矢集合N中的第i个元素n_i为所述设计模型点云中的第i个元素q_i对应的单位法向矢量，i＝1,2,…,l；

(1-3)根据所述工件设计模型，依次对所述设计模型点云中的点进行划分；其中当q_i属于所述工件设计模型的凹面时，q_i为凹面点；当q_i属于所述工件设计模型的凸面时，q_i为凸面点；i＝1,2,…,l；

(2)依次根据所述测量模型点云P中的点p_i，i＝1,2,…,n，在所述设计模型点云Q中提取与p_i距离最近的设计模型点云中的点q′_i，q′_i∈Q，以及与q′_i对应的单位法向矢量n′_i，n′_i∈N，形成最近点点集Q′＝{q′₁,q′₂,...q′_i,...q′_n}和最近点法矢集合N′＝{n′₁,n′₂,...n′_i,...n′_n}；该最近点点集中的点与该最近点法矢集合中的点、以及所述测量模型点云中的点一一对应；

(3)根据余量要求，建构基于距离方差的目标函数：

(3-1)当凹面余量c_a与凸面余量c_t均为常数值时，目标函数为

$F_1=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(d_i-c_i)}^2,$

其中d_i＝||R·p_i+t-q′_i||，R、t分别为3×3的矩阵和3×1的矩阵；当q′_i属于设计模型凹面时，c_i＝c_a；当q′_i属于设计模型凸面时，c_i＝c_t；

然后执行步骤(4)；

(3-2)当凹面余量c_a与凸面余量c_t的比值为常数值c时，目标函数为

$F_2=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_i^2{d}_i^2-\mathrm{\&lambda;}\frac{{\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_i{d}_i\right)}^2}{n}$

其中λ∈(0,1]，d_i＝||R·p_i+t-q′_i||，R、t分别为3×3的矩阵和3×1的矩阵；当q′_i属于设计模型凹面时，k_i＝1；当q′_i属于设计模型凸面时，k_i＝c；

然后执行步骤(4)；

(4)根据所述测量模型点云P、最近点点集Q′和最近点法矢集合N′，求解步骤(3)中的R和t，使所述目标函数的数值最小，实现点云匹配。

2.如权利要求1所述基于距离方差最小的工件点云匹配算法，其特征在于，所述步骤(4)包括以下步骤：

(4-1)当凹面余量c_a与凸面余量c_t均为常数值时，记

ξ＝(A^TA)^-1B，

其中， $A=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left[{n_i^{\&prime;}}^T\right|{(p_i\&times;n_i^{\&prime;})}^T],$

$B=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left[\right[{(p_i-q_i^{\&prime;\&prime;})}^T\&CenterDot;n_i^{\&prime;}\left]{\left[{n_i^{\&prime;}}^T\right|{(p_i\&times;n_i^{\&prime;})}^T]}^T\right];$ q″_i＝q′_i+c_in′_i，当q′_i属于设计模型凹面时，c_i＝c_a；当q′_i属于设计模型凸面时，c_i＝c_t；

接着，执行步骤(4-3)；

(4-2)当凹面余量c_a与凸面余量c_t的比值为常数值c时，记

$\mathrm{\&xi;}=-{[\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i^T{E}_i-\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i^T\right)\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i\right)/n]}^{-1}[\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{F}_i{E}_i^T-\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{F}_i\right)\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i^T\right)/n]$

其中，F_i＝k_i(p_i-q′_i)^T·n′i， $E_i=k_i\left[{n_i^{\&prime;}}^T\right|{(p_i\&times;n_i^{\&prime;})}^T],$ 当q′_i属于设计模型凹面时，k_i＝1；当q′_i属于设计模型凸面时，k_i＝c；

接着，执行步骤(4-3)；

(4-3)记矩阵v和矩阵ω，v矩阵和ω矩阵的行数相同，列数也相同，并且矩阵v和矩阵ω还满足 $\mathrm{\&xi;}=\left[\begin{array}{c}v\\ \mathrm{\&omega;}\end{array}\right],$

则t＝v， $R=e^{\widehat{\mathrm{\&omega;}}},$

其中，为ω的反对称矩阵；

(4-4)更新测量模型点云P：

记p′_i＝R·p_i+t，

并使p_i＝p′_i；i＝1,2,…,n；

更新后的测量模型点云P即为匹配后的工件测量点云。

3.如权利要求2所述基于距离方差最小的工件点云匹配算法，其特征在于，该基于距离方差最小的工件点云匹配算法，还包括以下步骤：

(5)计算误差判断函数：

(5-1)依次根据所述更新后的测量模型点云P中的点p_i，i＝1,2,…,n，在所述设计模型点云Q中提取与p_i距离最近的设计模型点云中的点q′_i，q′_i∈Q，以及与q′_i对应的单位法向矢量n′_i，n′_i∈N，形成最近点点集Q′＝{q′₁,q′₂,...q′_i,...q′_n}和最近点法矢集合N′＝{n′₁,n′₂,...n′_i,...n′_n}；该最近点点集中的点与该最近点法矢集合中的点、以及所述测量模型点云中的点一一对应；

(5-2)当凹面余量c_a与凸面余量c_t均为常数值时，均方误差Std为：

$\mathrm{Std}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(h_i-c_i)}^2}{n}},$

其中h_i＝(q′_i-p_i)^Tn′_i；当q′_i属于设计模型凹面时，c_i＝c_a；当q′_i属于设计模型凸面时，c_i＝c_t；

其中当Std大于预先设定的均方误差值，则执行步骤(5-4)；否则，测量模型点云P即为匹配后的工件点云，该基于距离方差最小的工件点云匹配算法结束；

(5-3)当凹面余量c_a与凸面余量c_t的比值为常数值c时，{q′₁,q′₂,...q′_i,...q′_n}中有m个q′_i属于设计模型凹面,有k个q′_i属于设计模型凸面，

则，均方误差Std、凹凸面余量比例平均误差p分别为：

$\mathrm{Std}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_i^2{h}_i^2-\mathrm{\&lambda;}{\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_i{h}_i\right)}^2/n}{n}},$

$p=|k\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ia}}/m\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{d}_{\mathrm{it}}-c|,$

其中，h_i＝(q′_i-p_i)^Tn′_i；当q′_i属于设计模型凹面时，k_i＝1；当q′_i属于设计模型凸面时，k_i＝c；

当q′_i属于设计模型凹面时，d_ia＝||R·p_i+t-q′_i||；当q′_i属于设计模型凸面时，d_it＝||R·p_i+t-q′_i||；

其中当Std大于预先设定的均方误差值，或者p大于预先设定的比例误差值，则执行步骤(5-4)；否则，测量模型点云P即为匹配后的工件点云，该基于距离方差最小的工件点云匹配算法结束；

(5-4)重复步骤(3)到步骤(5)。

4.如权利要求1-3任意一项所述基于距离方差最小的工件点云匹配算法，其特征在于，所述步骤(1-2)中，测量模型点云中的点的总数n和设计模型点云中的点的总数l满足l≥5×n。

5.如权利要求1-4任意一项所述基于距离方差最小的工件点云匹配算法，其特征在于，所述步骤(3-2)中，λ＝1。