本发明涉及数控加工技术领域，具体为一种基于测量点与自适应差分进化算法的工件姿态调整方法。
背景技术：
工件精确定位与余量分布均化是保证加工精度的关键，通过实时调整工件姿态，使其与理论姿态保持一致是实现精确定位与余量均布的直接有效的方法。余量分布对保证加工精度至关重要，而工件姿态对余量分布影响显著。在机械加工中，由于多次装夹等原因，工件姿态容易发生变化，因此有必要寻找一种方法来保证工件姿态的准确性。
关于匹配问题，besl最先提出了迭代最近点(icp)算法，该算法被广泛应用于点云配准，在特定的初始条件下能够较好地实现曲面匹配；在传统icp算法基础上，研究者提出了新的匹配算法，liu提出了一种新颖的自由曲面匹配算法，该方法通过基于共线性和接近性约束条件建立的传统icp准则来直接处理可能的点匹配，而不需要任何特征提取、图像预处理或异常数据的运动估计；du提出了一种基于期望最大化(em)估计的有界尺度的概率迭代最近点算法，用于噪声点集的各向同性尺度配准；针对自由曲面及其他类型曲面的匹配问题，zhang提出了一种点到最近点的迭代匹配算法，该算法采用基于距离分布的统计方法实现异常点的处理。然而，就icp算法而言，其鲁棒性和收敛速度有待提高，且采用该算法获得的解很可能是局部最优解。为了提高icp算法的收敛速度与运算精度，lena提出了收敛迭代最近点算法，该方法使得局部曲面匹配的精度大幅提升；针对3d点云，ge提出了基于非线性gauss-helmert模型和非线性最小二乘的点云数据配准方法，并完成相关三维模型的重构，该方法在收敛速度和配准精度上优于icp算法。在机械加工中，曲面匹配主要用于复杂零件的定位与余量优化，在曲面匹配的实际应用中，众多研究者也进行了相关的研究工作。yan为解决数控加工中复杂曲面类零件的定位问题，提出了两步匹配法，算法分为粗匹配与精确匹配，在精确匹配中引入了遗传算法(ga)；ko提出了一种针对自由曲面造型三维物体的匹配方法，该方法将匹配问题归结为采用区间投影法(ipp)求解非线性方程系统，并将曲面的高斯曲率和平均曲率作为匹配特征，该方法可实现全局及局部特征的匹配；mehrad通过粗定位与精确定位，实现了基于曲率和距离相似性的点到设计模型的匹配，该方法稳定性较好，但由于b样条曲面的重构，该方法在运行效率上有改进空间；shen针对大型的复杂自由曲面类毛坯件的最佳定位，提出了最佳匹配算法，该算法极大地改善了毛坯件余量分布，提高了加工效率；zhang提出了基于三坐标测量的加工余量优化模型，该方法通过坐标对齐以确保毛坯件具有足够的加工余量；sun针对复杂曲面加工提出了一种统一定位模型，并通过bfgs算法实现余量优化，该方法也涉及到了匹配过程，以满足用户定义的需优先保证加工余量的特定表面的约束。
上述匹配算法大多是针对点云与设计模型的匹配，且大多算法需要在已知曲面参数方程或其他曲面特性参数的情况下实现。对于不规则曲面，其曲面参数方程的确定较为困难。
技术实现要素：
针对上述问题，本发明的目的在于提供一种基于测量点与自适应差分进化算法(self-adaptivedifferentialevolutionalgorithm，后面简称sade)的工件姿态调整方法，通过对工件姿态的实时调整，辅助完成工件的精确定位，以确保工件每一工序余量分布均匀，提高加工精度；该方法只需少数测量点，因此对规则与非规则曲面都具有良好的适用性。技术方案如下：
一种基于测量点与自适应差分进化算法的工件姿态调整方法，包括：
进行点集匹配：设工件特征面实际测量点集p＝{pi|i＝1,2,…,n}，理论测量点集q＝{qi|i＝1,2,…,n}，则p到q的匹配表示为：
式中，r为空间旋转变换矩阵，α、β、γ分别表示实际测量点相对于xyz轴的旋转角度；t为平移变换矩阵，tx、ty、tz分别表示实际测量点相对于xyz轴的平移量；匹配流程如下：
step1：获取测量点集q和p；
step2：在理论曲面与对应的实际曲面上分别选取3个对应的参考点，建立理论与实际局部坐标系；
step3：初始匹配：通过局部坐标系的对齐实现初始匹配，并计算获取初始匹配矩阵r0和t0；
step4：点集变换：利用下式对实际点集p进行变换，获得新的实际点集p'；
p′＝pr0+t0(2)
step5：精确匹配：采用sade算法搜索最佳精确匹配矩阵r1和t1；
step6：按以下公式计算最终的变换矩阵r和t；
step7：计算机床调整参数并对机床进行调整；
step8：重新测量并获取新的点集p”，计算点集p”与理论点集q之间的平均距离dist_ad，并与工件加工面容差ε进行比较，如果dist_ad≤|ε|，结束流程，否则令p'＝p”，转至step5重新计算。
进一步的，所述局部坐标系建立的方法为：
在理论曲面上选取不共线的q0，q1，q2作为参考点，在实际曲面选取不共线的p0，p1，p2作为参考点；根据上述参考点确定各局部坐标系xyz轴单位方向矢量，据此建立局部坐标系，具体如下：
设ia，ja，ka为实际局部坐标系xyz轴单位方向矢量，p0为原点；it，jt，kt为理论局部坐标系xyz轴单位方向矢量，q0为原点；则有：
根据公式(4)、(5)求得坐标系各轴的单位方向矢量，建立实际局部坐标系p0-xayaza和理论局部坐标系q0-xtytzt。
更进一步的，所述初始匹配矩阵的求解方法包括：
初始匹配通过局部坐标系p0-xayaza与q0-xtytzt的对齐实现，两坐标系与初始匹配矩阵之间存在以下关系：
式中r0为初始匹配矩阵，其形式如下：
根据r0，实际曲面与理论曲面之间的方向误差，即实际测量点相对于理论坐标轴xt,yt和zt的旋转角度α0,β0和γ0，如下式所示：
式中，r0(i,j)为矩阵r0第i行第j列元素(i,j＝1,2,3)；
坐标系p0-xayaza相对于坐标系q0-xtytzt的平移量由下式计算得到：
t0＝μ-μ0′r0＝[tx0ty0tz0](8)
式中，
更进一步的，所述实际曲面与理论曲面之间的初始姿态误差通过下式所示姿态误差矩阵pose0描述：
通过pose0能够直观地表示实际曲面相对于理论曲面的旋转角度值与平移量，理想情况下，实际曲面与理论曲面完全重合，因此pose0为零矩阵。
更进一步的，所述精确匹配矩阵的求解方法为：
精确匹配旋转变换与平移变换矩阵分别为r1和t1，其中：
t1＝[tx1ty1tz1]
初始匹配后，实际点集由p转换为p′＝{p′i|i＝1,2,…,n}，将精确匹配归结为一个最小化问题，依据最小二乘原理确定其目标函数如下：
由式(10)知，目标函数包含α1,β1,γ1和tx1,ty1,tz16个变量，由式(8)可知，初始匹配后的实际点集p'相对于理论点集q的平移量可表示为：t1＝μ-μ1′r1＝[tx1ty1tz1]，其中p'的质心坐标
更进一步的，在所述精确匹配中，引入sade算法对目标函数进行求解，其中的自适应变异算子如下式所示：
式中，f为变异算子，f0为初始变异算子，gm为最大进化代数，g表示当前进化代数。
更进一步的，所述step7具体为：
通过最终的空间变换矩阵r和t，根据五轴数控机床结构类型求解出对应机床的调整量，包括旋转轴的旋转角度以及平移轴的平移量；根据调整量对机床进行相应的调整，以使工件姿态与匹配结果一致，进而保证余量分布均匀；
在czfxyb结构的五轴数控机床中，设绕c轴转角为δc，绕b轴转角δb，沿xyz轴的平移量分别为tx、ty、tz，则有：
t＝t0+t1＝[txtytz](14)
由式(13)、(14)求得机床调整量δb、δc及tx、ty、tz，实现机床的调整，进而完成工件姿态的调整。
本发明的有益效果是：本发明的方法分初始匹配与精确匹配，通过两次匹配并结合数控机床的结构，可以求得最佳机床调整参数，最后对机床进行调整以实现工件姿态的调整。可成功用于iso五轴数控机床典型加工测试件(isocd10791-7:2017,testconditionsformachiningcentres-part7:accuracyoffinishedtestpiece)姿态的调整，确保其定位的准确性与余量分布的均匀性；试验证明，本方法不需要已知工件曲面参数方程或其他参数，只需少数测量点即可完成匹配，实现姿态调整，操作简单可靠，具有更好的适用性；且相对于ga和bfgs算法，sade算法更加简便、可靠和稳定，其收敛速度较其他进化算法更快，且不受初值的影响；此外，对大空间、非线性和不可求导的连续问题，其求解效率较其他进化算法好。
附图说明
图1是匹配流程图。
图2是局部坐标系示意。
图3是sade算法流程图。
图4是测试件模型。
图5是初始姿态示意。
图6是sade优化曲线图。
图7是姿态调整前后点集间的距离。
具体实施方式
下面结合附图和具体实施例对本发明做进一步详细说明。本发明提出了一种新的匹配算法以实现对工件姿态的实时调整。所提算法包括初始匹配与精确匹配，在精确匹配环节引入自适应差分进化算法(sade)来搜索最佳精确匹配矩阵，通过分步匹配求得实际点集相对于理论点集的最佳空间旋转变换矩阵r和平移变换矩阵t。最后，根据最佳变换矩阵结合机床结构求得机床调整参数，并实现工件姿态的调整。
进行点集匹配：
设工件特征面实际测量点集p＝{pi|i＝1,2,…,n}，理论测量点集q＝{qi|i＝1,2,…,n}，则p到q的匹配可以表示为：
式中，r为空间旋转变换矩阵，α、β、γ分别表示实际测量点相对于xyz轴的旋转角度；t为平移变换矩阵，tx、ty、tz分别表示实际测量点相对于xyz轴的平移量。图1所示为匹配流程图，具体流程如下：
step1.获取测量点集q和p。
step2.在理论曲面与对应的实际曲面上分别选取3个对应的参考点，建立理论与实际局部坐标系。
如图2所示，理论曲面与实际曲面之间存在姿态误差。在理论曲面与实际曲面上分别选取q0，q1，q2和p0，p1，p2作为参考点。通常情况下，参考点由3个不共线的点构成，且尽可能地分布其他测量点的外侧。根据参考点可以确定各局部坐标系xyz轴单位方向矢量，据此建立局部坐标系，具体方法如下：
设ia，ja，ka为实际局部坐标系xyz轴单位方向矢量，p0为原点；it，jt，kt为理论局部坐标系xyz轴单位方向矢量，q0为原点。则有：
根据公式(2)、(3)可以求得坐标系各轴的单位方向矢量，建立如图2所示的实际局部坐标系p0-xayaza和理论局部坐标系q0-xtytzt。
step3.初始匹配：通过局部坐标系的对齐实现初始匹配，并计算获取初始匹配矩阵r0和t0。
初始匹配通过局部坐标系p0-xayaza与q0-xtytzt的对齐实现，经过初始匹配，可近似确定后续精确匹配中变量的边界值。理论上，初始匹配后实际局部坐标系与理论局部坐标系完全重合，因此两局部坐标系与初始匹配矩阵之间存在式(4)所示关系：
式中r0为初始匹配矩阵，其形式如下：
根据r0，实际曲面与理论曲面之间的方向误差(即实际测量点相对于理论坐标轴xt,yt和zt的旋转角度α0,β0和γ0可以求得，其结果如下式(5)所示：
式中，r0(i,j)为矩阵r0第i行第j列元素(i,j＝1,2,3)。
坐标系p0-xayaza相对于q0-xtytzt的平移量可由式(6)计算得到。
t0＝μ-μ′0r0＝[tx0ty0tz0](6)
式中，
求得α0,β0,γ0,tx0,ty0,tz0后，实际曲面与理论曲面之间的初始姿态误差可以通过式(7)所示姿态误差矩阵pose0描述：
通过pose0可以直观地表示实际曲面相对于理论曲面的旋转角度值与平移量，理想情况下，实际曲面与理论曲面完全重合，因此pose0为零矩阵。
step4.点集变换：利用公式(8)对实际点集p进行变换，获得新的实际点集p'；
p′＝pr0+t0(8)
step5.精确匹配：采用sade算法搜索最佳精确匹配矩阵r1和t1。
设精确匹配旋转变换与平移变换矩阵分别为r1和t1，其中：
t1＝[tx1ty1tz1](9)
初始匹配后，实际点集由p转换为p′＝{pi′i＝1,2,…,n}，精确匹配是为了寻找最佳匹配矩阵r1和t1，使匹配后的点集到理论点集q距离的均方误差最小，因此，精确匹配可以归结为一个最小化问题，依据最小二乘原理可确定其目标函数如式(10)所示。
由式(10)可知，目标函数包含α1,β1,γ1和tx1,ty1,tz16个变量，且属于非线性问题。为了提高运算效率，有必要对该函数进行降元处理，由式(6)可知，初始匹配后的实际点集p'相对于理论点集q的平移量可以表示为：t1＝μ-μ′1r1＝[tx1ty1tz1]，其中p'的质心坐标
在精确匹配中，本发明首次引入sade算法对目标函数进行求解，以获得最佳矩阵参数α1,β1,γ1和tx1,ty1,tz1。sade是传统差分进化算法(de)的改进，其操作主要包括交叉、变异、选择和边界条件处理。在sade中，可以灵活设计自适应变异算子，能够有效避免因变异率过大造成的搜索效率低下和全局最优解精度低，以及因变异率过小导致的种群多样性降低进而出现早熟等现象，自适应变异算子设计如式(12)所示。
式中，f为变异算子，f0为初始变异算子，gm为最大进化代数，g表示当前进化代数。在算法开始时自适应变异算子为2f0，具有较大值，在初期可以保持个体多样性，避免早熟；随着算法的进展，变异算子逐步降低，到后期其值接近f0，保留优良信息，避免最优解遭到破坏，增加搜索到全局最优解的概率。图3为sade流程图。
step6.按公式计算最终的变换矩阵r和t；
step7.计算机床调整参数并对机床进行调整。
根据初始匹配矩阵r0和t0及精确匹配矩阵r1和t1，最终的空间变换矩阵r和t可由式(13)计算获得。
通过r和t，根据五轴数控机床结构类型可以求解出对应机床的调整量，包括旋转轴的旋转角度以及平移轴的平移量。根据调整量对机床进行相应的调整，以使工件姿态与匹配结果一致，进而保证余量分布均匀。以czfxyb(工作台主轴双回转)结构的五轴数控机床为例，设绕c轴转角为δc，绕b轴转角δb，沿xyz轴的平移量分别为tx、ty、tz，则有[20]：
t＝t0+t1＝[txtytz](21)
由式(14)、(15)可以求得机床调整量δb、δc及tx、ty、tz，实现机床的调整，进而完成工件姿态的调整。
以iso五轴数控机床典型加工测试件为例，对本发明所开发的算法进行验证。图4所示为测试件模型，加工曲面容差为ε＝±0.035mm，加工曲面初始姿态如图5所示。具体实施步骤及相关测量数据和试验结果如下。
1)初始匹配
采用三坐标测量机对加工曲面进行测量，获取理论测量点集q与对应的实际测量点集p，选取参考点建立局部坐标系，测量数据如表1所示。点集初始平均距离dist_in＝2.02194mm，初始匹配结果如表2所示。
2)精确匹配
首先，根据初始匹配后姿态误差矩阵pose0中α0,β0,γ0的值确定精确匹配中α1,β1,γ1的边界值，理论上，其边界值不会超过α0,β0,γ0，但考虑到匹配中存在误差，可在α0,β0,γ0的基础上适当扩大其边界值。则α1,β1,γ1的边界值可分别定义为：[-α0-k1,α0+k1],[-β0-k2,β0+k2],[-γ0-k3,γ0+k3]，其中ki(i＝1,2,3)为[0°,3°]的任意常数。在实例中优化变量边界值如下：
[α1(l),α1(u)]＝[-5°,5°]
[β1(l),β1(u)]＝[-5.5°,5.5°]
[γ1(l),γ1(u)]＝[-3°,3°]
根据本发明确定的目标函数及优化变量，结合文献[13][17]确定sade算法相关的参数值如表3所示。图6所示为sade优化曲线，经过200次进化获得全局最优解α1＝-0.02994048°,β1＝0.00881839°，γ1＝0.00256112°。精确匹配矩阵r1、t1如下所示。
t1＝[0.3378890150.037944696-0.072290709]
3)机床调整
根据公式(3)可求得最终空间变换矩阵r和t，再结合公式(14)、(15)可以求得机床调整参数δc,δb,tx,ty,tz。最后对机床进行调整，以实现工件姿态的调整。dist_ad为姿态调整后理论点集与实际点集之间的平均距离，用以评估调整效果。调整参数及调整后的测量结果分别见表4、表5。
4)结果分析
图7所示为匹配前、初始匹配后及姿态调整后实际点集与理论点集之间的平均距离，从中可知，经过两次匹配，实际点集与理论点集间的误差大幅降低，最终收敛到加工曲面容差范围内，调整后点集间距离dist_ad＝0.028077mm<|ε|。
表1初始测量数据(dist_in表示初始点集间平均距离/mm)
表2初始匹配结果
表3sade算法参数值
表4调整参数及计算结果
表5初始匹配与调整后的测量点集