定位精度是工业机器人性能的重要衡量指标之一^[1]，是决定机器人执行任务精确度、效率和可重复性的关键因素. 造成工业机器人的定位误差的因素可分为几何误差与非几何误差^[2-3]，其中，几何误差占总定位误差的80%^[4]. 因此，设计有效的校正方法降低几何误差进而提升工业机器人定位精度具有重要意义.

在降低误差的方法中，误差补偿法具有成本较低且效率较高的特点，该方法包括4个步骤：运动学建模、参数测量、参数辨识、误差补偿.

运动学建模是轨迹实现和精度预测的基础，常见模型包括D-H^[5]、CGA^[6]、CPC^[7] 和POE^[8]等. Long等^[9]建立的修正D-H模型解决了相邻关节轴平行时的奇异性问题. He等^[10]提出局部POE误差模型. 相较于D-H及M-DH模型需要24个运动学参数，POE模型仅需要关节旋量与转角，共12个参数^[11]，且指数积公式满足连续性与完整性^[12]，适用于机器人运动学模型建立.

误差参数测量为后续的误差补偿提供必要的数据. 双目视觉测量系统成本较低、可跟踪6自由度，但精度最高仅为0.1 mm，不适用于高精度场景^[13]. 激光跟踪仪测量空间大，精度达10~100 μm，但成本较高^[14]. 球杆仪实验成本相对较低且在精度方面可以达到1 μm^[15]. Ferrarini等^[16]用跟踪仪实验实现路径误差补偿. 黄鹏等^[17] 提出加权点匹配法实现激光跟踪仪坐标转换. Jiang等^[18]设计球杆仪“S”型轨迹简化误差辨识. 赵业和等^[19]结合多种群竞争松鼠搜索算法规划最短时间轨迹. 上述研究均未考虑测量空间中定位误差的相似特性，该特性会导致运动轨迹偏差，影响测量结果准确性.

误差参数辨识是提高机器人精度和稳定性的核心. 常用的方法为最小二乘法^[20]. Zhu等^[21]采用Levenberg-Marquardt算法辨识双臂串联机器人运动学误差参数. Cao等^[22]提出基于蝴蝶花授粉算法的人工神经网络标定方法以强化全局最优搜索能力. Zhao等^[23]提出自适应权重与阈值平衡的串联机械臂参数辨识方法. Zhang等^[24] 通过自适应布谷鸟搜索算法实现动态约束下机械臂关节空间运动用时最小化. 尽管上述研究均实现了误差参数辨识，但其所采用的算法普遍存在对初始值敏感、搜索效率低、收敛速度慢且易陷入局部最优等问题.

误差补偿是提升机器人精度和性能的关键环节. Li等^[25]基于神经网络建立误差预测模型实现补偿. Li等^[26]通过CDTW算法提出轨迹实时补偿方法. 张泽坤等^[27]采用高精度鲁棒控制提高了在线补偿精度与抗干扰能力. Zhang等^[28]利用绝对式直线光栅尺结合KUKA控制环境实现定位误差补偿. 以上补偿方法均取得了显著成效，但均未考虑测量设备重复定位误差导致的测量不确定度问题，该因素会导致补偿过程中机器人运动状态偏差，进而影响补偿效率.

针对上述问题，本研究以六自由度工业机器人为对象开展测量轨迹规划及考虑测量不确定度的误差参数辨识研究，具体贡献如下：1) 提出空间圆弧测量轨迹与测量模型；2) 建立改进的飞蛾火焰优化逆运动学求解方法；3) 提出考虑不确定度的补偿参数量化嵌套算法.

以6自由度工业机器人为研究对象，其零位构型及关节坐标系如图1所示.

基于POE模型对关节坐标系定义，z轴为关节轴线，x轴为相邻轴线公垂线，y轴依右手定则确定，原点为旋转关节中心. 图1中{A}为基础坐标系，{B}为工具坐标系，J_i (i=1,2,···,6)为旋转关节，箭头方向为旋转方向.

第i个旋转关节的关节旋量表达式如下：

(1)$ {\boldsymbol{\zeta}}_i=\left[-{\bf\textit{ω}}_i \times {\boldsymbol{q}}_i ,\quad {\bf\textit{ω}}_i\right]^{\mathrm{T}} . $

式中：ω_i为第i个旋转关节轴线方向的单位矢量，q_i为第i个旋转关节的轴线位置矢量. 机器人正运动学使用的参数如表1所示.

机器人正向运动学模型表达式如下：

(2)$ {\boldsymbol{g}}_{{AB}}(\theta)={\rm{e}}^{\hat{{\boldsymbol{\zeta}}}_1 \theta_1} {\rm{e}}^{\hat{{\boldsymbol{\zeta}}}_2 \theta_2} \cdots {\rm{e}}^{\hat{{\boldsymbol{\zeta}}}_i \theta_i} {\boldsymbol{g}}_{{AB}}(0). $

式中：g_AB(0)为机器人初始位姿下末端工具坐标系对于基础坐标系的刚体变换矩阵，$\theta_i $为第i个关节的关节变量，$ {{{\hat {\boldsymbol{\zeta}} }}_i} $为运动旋量矩阵. $ {{{\hat{ {\boldsymbol{\zeta} }}}}_i} $计算过程如下：

(3)$ {{\hat {\boldsymbol{\zeta}} }_i} = \left[ {\begin{array}{*{20}{c}} {{{{\hat {{\bf\textit{ω}}} }}_i}}&{ - {{{\bf\textit{ω}} }_i} \times {{q}_i}+h{{\bf\textit{ω}} }} \\ 0&{\bf{0}} \end{array}} \right]. $

式中：h为旋量的平移分量，当关节为旋转关节时，h=0，ω_i =[ω_iX, ω_iY, ω_iZ]^T，||ω_i||=1；$ {{\hat {\bf\textit{ω}} }_i} $为ω_i在反对称矩阵空间的映射. $ {{\hat {\bf\textit{ω}} }_i} $表达式如下：

(4)$ {{\hat {\bf\textit{ω}} }_i} = \left[ {\begin{array}{*{20}{c}} 0&{ - {{{{\omega}} }_{i{Z} }}}&{{{{{\omega}} }_{i{Y} }}} \\ {{{{{\omega}} }_{i{Z} }}}&0&{ - {{{{\omega}} }_{i{X} }}} \\ { - {{{{\omega}} }_{i{Y} }}}&{{{{{\omega}} }_{i{X} }}}&0 \end{array}} \right]. $

式(2)中的$ {{\rm{e}}^{{{{\hat {\boldsymbol{\zeta}} }}_i}{\theta _i}}} $表达式如下：

(5)$ {{\rm{e}}^{{{{\hat {\boldsymbol{\zeta}} }}_i}{\theta _i}}} = \left[ {\begin{array}{*{20}{c}} {{{\rm{e}}^{{{{\hat {\bf\textit{ω}} }}_i}{\theta _i}}}}&{({{\boldsymbol{I}}} - {{\rm{e}}^{{{{\hat {\bf\textit{ω}} }}_i}{\theta _i}}})({{{\bf\textit{ω}} }_i} \times {{{\boldsymbol{v}}}_i})+{{{\bf\textit{ω}} }_i}{{{\bf\textit{ω}} }_i}^{\mathrm{T}}{{{\boldsymbol{v}}}_i}{\theta _i}} \\ 0&{\bf{1}} \end{array}} \right]. $

式中：${\boldsymbol{v}}_i $为第i个关节旋量的平移向量.

由罗德里格斯公式可知，$ {{\rm{e}}^{{{{\hat {\bf\textit{ω}} }}_i}{\theta _i}}} $表达式如下：

(6)$ {{\rm{e}}^{{{{\hat {\bf\textit{ω}} }}_i}{\theta _i}}} = {{\boldsymbol{I}}}+(\sin \;{\theta _i}) {{{{\hat {\bf\textit{ω}} }}_i}} +(1 - \cos\; {\theta _i}){\left[ {{{{\hat {\bf\textit{ω}} }}_i}} \right]^2}. $

式中：I为单位矩阵. 由此可得机器人运动学正解方程如下：

(7)$ \begin{split} {{\boldsymbol{g}}_{AB}}(\theta ) =& {{\rm{e}}^{{{{\hat {\boldsymbol{\zeta}} }}_1}{\theta _1}}}{{\rm{e}}^{{{{\hat {\boldsymbol{\zeta}} }}_2}{\theta _2}}} \cdot \cdot \cdot {{\rm{e}}^{{{{\hat {\boldsymbol{\zeta}} }}_i}{\theta _i}}}{{\boldsymbol{g}}_{AB}}(0) =\\ & \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{R}}}(\theta )}&{{{\boldsymbol{p}}}(\theta )} \\ {\bf{0}}&1 \end{array}} \right]. \end{split} $

式中：R(θ)为机器人末端姿态矩阵，p(θ)为机器人末端位置矩阵.

(8)$ {{\boldsymbol{R}}}(\theta ) = \left[ {\begin{array}{*{20}{c}} {{r_{11}}}&{{r_{12}}}&{{r_{13}}} \\ {{r_{21}}}&{{r_{22}}}&{{r_{23}}} \\ {{r_{31}}}&{{r_{22}}}&{{r_{33}}} \end{array}} \right], $

(9)$ {{\boldsymbol{p}}}(\theta ) = {\left[ {\begin{array}{*{20}{c}} {{p_{x} }},&{{p_y}},&{{p_{z} }} \end{array}} \right]^{\text{T}}}. $

将i=0代入式(7)中，即各关节变量理论值均为零时得到机器人处于零位时的刚体变换矩阵g_AB(0)表达式如下：

(10)$ {{\boldsymbol{g}}_{AB}}(0) = \left[ {\begin{array}{*{20}{c}} { - 1}&0&0&{374} \\ 0&1&0&0 \\ 0&0&{ - 1}&{630} \\ 0&0&0&1 \end{array}} \right]. $

对不同运动学模型在构建误差模型时使用的参数进行对比，如表2所示. 可以看出，POE模型参数少于D-H模型与M-DH模型参数. 旋量坐标偏差Δζ_i由关节轴线方向误差Δω_i与关节轴线位置误差Δq_i引起，而关节轴线位置参数q_i由连杆长度参数l_i确定. 关节轴线的平行度与垂直度误差区间为[0, 0.05] mm，关节角度误差区间为[0,0.3°]，l_i区间为[0, 0.5] mm^[32]，轴线偏转角度表达式如下：

(11)$ {S_i} = {w}/{{{l_i}}}. $

式中：w为平行度与垂直度误差.

在考虑不同误差参数的情况下，随机选取100组误差值代入机器人正运动学模型，得到末端定位误差平均值及影响比例，如表3所示. 其中，e_k为考虑参数时末端误差，e_b为不考虑参数时末端误差，ρ为影响比例. 可以看出，关节轴线方向参数误差对定位精度影响较小，因此本研究在建立运动学误差模型时，仅考虑连杆参数误差Δl_i和关节变量偏差Δθ_i. 考虑误差的POE模型参数如表4所示. 其中，$l_i^{\mathrm{e}} $为含有误差的机器人连杆长度参数，$\theta ^{\mathrm{e}}_i $为含有误差的机器人关节变量.

在考虑误差的情况下，旋量的指数表达式如下：

(12)$ {{\rm{e}}^{{{{\hat {\boldsymbol{\zeta}} }}_i}\Delta {\theta _i}}} = \left[ {\begin{array}{*{20}{c}} {{{\rm{e}}^{{{{\hat {\bf\textit{ω}} }}_i}\Delta {\theta _i}}}}&{({\boldsymbol{I}} - {{\rm{e}}^{{{{\hat {\bf\textit{ω}} }}_i}\Delta {\theta _i}}})({{{\bf\textit{ω}} }_i} \times {{{\boldsymbol{v}}}_i})+{{{\bf\textit{ω}} }_i}{{{\bf\textit{ω}} }_i}^{\mathrm{T}}{{{\boldsymbol{v}}}_i}\Delta {\theta _i}} \\ {\bf{ 0}}&{\bf{1}} \end{array}} \right]. $

由表4和式(11)可以得到基于POE模型的工业机器人误差模型：

(13)$ {{\boldsymbol{g}}_{AB}}({\theta ^{\mathrm{e}}}) = {{\rm{e}}^{{{{\hat {\boldsymbol{\zeta}} }}_1}\theta _1^{\mathrm{e}}}}{{\rm{e}}^{{{{\hat {\boldsymbol{\zeta}} }}_2}\theta _2^{\mathrm{e}}}} \cdot \cdot \cdot {{\rm{e}}^{{{{\hat {\boldsymbol{\zeta}} }}_i}\theta _i^{\mathrm{e}}}}{{\boldsymbol{g}}_{AB}}(0). $

式中：$\theta^{\rm{e}} $为含误差的关节变量.

运动学逆解是在已知末端工具坐标系{B}相对于基础坐标系{A}的理论姿态矩阵${\boldsymbol{g}}_0^6 $的前提下，求解各关节转角. 使用飞蛾火焰优化（moth-flame optimization, MFO）算法进行逆运动学求解. MFO具有收敛性强、兼容性高、鲁棒性好的特点. 在MFO算法中，飞蛾作为搜索代理，火焰代表历史最优位置，飞蛾围绕火焰搜索并更新更优解，确保不丢失最优解. 飞蛾搜索最优解的轨迹S以及第i个飞蛾到第j个火焰的距离D_i表达式分别如下：

(14)$ S\left( {{M_i},{F_j}} \right) = {D_i} {{{\mathrm{e}}}^{{bt}}} \cos\; (2{\text{π}} t)+{F_j}, $

(15)$ {D_i} = \left| {{F_j} - {M_i}} \right|. $

式中：M_i表示i个飞蛾，F_j表示第j个火焰，b为定义对数螺旋的常数，t为[−1.0, 1.0]的随机数. MFO算法的位置更新公式如下：

(16)$ {{{M}}^{(i,j)}} = D \times {{{\mathrm{e}}}^{{bt}}} \times \cos\; (2{\text{π}} T)+{{{S}}^{(i,j)}}. $

式中：M^(i,j)为第i只飞蛾在维度j的位置，S^(i,j)为适应度排序后的蛾群位置数组，参数D为第i只飞蛾在第j维度上对应火焰参数之间的绝对距离，T为位置更新随机参数. D和T的表达式分别如下：

(17)$ D = \left| {{{{S}}^{(i,j)}} - {{{M}}^{(i,j)}}} \right|, $

(18)$ T = (A - 1) \; {T_0}+1. $

式中：T₀为[0,1.0]的随机数，A为线性递减参数.

(19)$ A=-1+K\left({-1}/{K_{\max }}\right) . $

式中：K为当前迭代次数，K_max为最大迭代次数. 在对工业机器人进行逆解时，对已知的末端工具坐标系{B}，相对于基础坐标系{A}构建期望目标矩阵$ {}_B^A{\boldsymbol{g}} $，将6个关节转角带入正运动学模型中，得到含有误差参数的末端位姿矩阵$ {}_B^A{{\boldsymbol{g}}_\theta } $. 构建差值矩阵$ {}_B^A{{\boldsymbol{g}}_{\mathrm{R}}} $，表达式如下：

(20)$ {}_B^A{{\boldsymbol{g}}_{\mathrm{R}}} = {}_B^A{{\boldsymbol{g}}_\theta } - {}_B^A{\boldsymbol{g}}. $

构建运动学逆解目标函数R_t，将姿态矩阵与位置矩阵进行拆分，目标函数如下：

(21)$ {\boldsymbol{R}}_{\mathrm{t}}=\sum_{i=1}^3 \sum_{j=1}^3\left({ }_{{B}}^{{A}} \boldsymbol{\boldsymbol{g}}_{\mathrm{R}(i, j)}\right)^2+\sum_{i=1}^3\left({ }_{{B}}^{{A}} \boldsymbol{\boldsymbol{g}}_{\mathrm{R}(i, 4)}\right)^2 . $

式中：$ {}_B^A{{\boldsymbol{g}}_{R} }_{(i,j)} $为姿态矩阵，$ {}_B^A{{\boldsymbol{g}}_R}_{(i,4)} $为位置矩阵.

在使用MFO进行机器人逆解计算时，问题的维度为12，包含6个关节角度参数与6个连杆长度参数. 为了评估优化算法性能，对MFO算法、遗传算法（GA）和粒子群优化算法（PSO）进行比较，定义初始种群为500、最大迭代次数为200. 对比结果如图2所示. 其中，n为迭代次数，c为目标函数值. 可以看出，MFO算法在约50次迭代时达到收敛阈值，PSO需约80次，GA收敛性较弱. MFO迭代曲线平滑，较少陷入局部最优，在收敛速度和避免局部最优方面优于GA和PSO.

通过分析机器人定位误差在笛卡尔和关节坐标空间中的空间相似性，推导测量轨迹，以降低采样点数据自相关性和样本偏差，提高样本精度. 在n自由度的机器人活动范围内，机器人定位误差的相似程度可以用变差函数表征：

(22)$ \gamma(\theta, h)=\frac{1}{2} E[\Delta P(\theta)-\Delta P(\theta+h)]^2 . $

式中：h为关节空间中2组关节的欧氏距离，作为关节分割量，即2组关节角度的广义距离；ΔP(θ+h)为距离坐标θ分割量为h时的定位误差.

变差函数值反映机器人定位误差空间相似度，值越小表明相似度越高. 假设在误差增量期望为0且方差平稳的条件下，变差函数值仅与关节输入增量相关. 不同关节分割量对误差相似性作用效果不同. 2组不同关节角度θ⁽ⁱ⁾与θ^(j)关节分割量h表达式如下：

(23)$ h = \left[ {{\sum\limits_{k = 1}^n {\left[ {\theta _k^{(i)} - \theta _k^{(j)}} \right]} }^2}\right]^{1/2} . $

式中：$ \theta _k^{(i)} $与$ \theta _k^{(j)} $分别为2组不同的关节输入量，i和j分别表示第i、j组关节角度.

已知单侧测量轨迹关节输入，基于POE正运动学模型推导理论位置坐标. 给定随机误差参数$ \Delta l_{{\mathrm{s}}} ^i $与$ \Delta \theta _{{\mathrm{s}}} ^i $，如表5所示.

使用表5中的随机误差参数与式(23)计算每隔5°关节分割量的实验变差函数值$ \gamma $，得到X、Y、Z三向变差函数值，如图3所示. 可以看出，变差函数极值在关节分割量105°附近. 在半球空间内，通过水平旋转单侧测量轨迹得到与单侧轨迹对应的对比轨迹，轨迹间夹角由0°水平旋转至180°后减小回初始位置. 当夹角为180°时，轨迹沿半球直径对称，对应最大关节分割量为194.57°，当选取3、4、5条轨迹在测量空间内均匀分布时，关节分割量分别为134.54°、87.39°、60.43°. 随测量轨迹增多，相邻轨迹间关节分割量减小，关节角度相似度提高，定位误差差异减小. 当选择3条轨迹时，相邻轨迹关节分割量为134.54°，定位误差相似性低于4条、5条轨迹的. 继续增加轨迹数量，关节分割量减小，变差函数值降低，定位误差相似性增大.

将误差信息不充分定义为测量时无法有效捕捉机械臂末端误差场全貌，导致误差分析信息不完整、无法代表真实误差分布或变动，难以充分反映完整测量空间内机械臂末端误差空间特征的现象. 所选轨迹测量误差范围小于随机误差范围，且误差场分布与完整空间差异较大，表明误差信息不充分. 这种局限性使得基于有限轨迹的误差分析会偏离实际，影响精度和有效性. 为了确保信息充分准确，须多轨迹全面反映误差特征.

在完整测量空间随机生成200组关节角度值，使用表5给出的随机误差参数，比较随机采样与2、3、4条轨迹测量下机械臂末端X、Y、Z方向误差的均值μ和方差σ，结果如图4所示.

如图4(a)所示，采用随机采样方式测得的误差均值曲线，与使用3条轨迹进行测量所得到的误差均值曲线最为接近. 误差较低且稳定；采用2条和4条轨迹，误差较高且波动大，分别采用2条、3条、4条轨迹进行测量时，测得的误差均值变化率分别为5.77%、85.34%、15.19%. 图4(b)的方差分析呈现类似特征，采用随机采样方式测得的方差曲线，与使用3条轨迹进行测量所得到的方差曲线最接近，2、3、4条轨迹的方差变化率分别为61.76%、1.96%、4.91%. 以随机采样为参考，3条轨迹在控制误差方差方面表现最佳，能够有效降低不确定性；4条轨迹效果不稳定；2条轨迹表现最差，可能由信息不足导致. 综上所述，采用3条轨迹测量，误差信息更完备，离散程度低，测量更可靠. 3条测量轨迹的示意图如图5所示.

在3条轨迹上共取f个采样点，选择半球球心位置为球杆仪底座位置，在机器人按采样轨迹运动时，固定姿态不变，姿态矩阵R和位置误差矩阵p表达式分别如下：

(24)$ {{\boldsymbol{R}}} = \left[ {\begin{array}{*{20}{c}} { - 1}&0&0 \\ 0&1&0 \\ 0&0&{ - 1} \end{array}} \right], $

(25)$ {{\boldsymbol{p}}} = {\left[ {\begin{array}{*{20}{c}} {\Delta {p_{x} }},&{\Delta {p_{y} }},&{\Delta {p_{z} }} \end{array}} \right]^{\mathrm{T}}}. $

设半球球心位置为(X_o, Y_o, Z_o)，在实际情况中，球心位置存在偏差ΔX_o、ΔY_o、ΔZ_o，则球心实际坐标如下：

(26)$ \left.\begin{array}{l}X_{\mathrm{o}}^{\mathrm{a}}=X_{\mathrm{o}}+\Delta X_{\mathrm{o}}, \\Y_{\mathrm{o}}^{\mathrm{a}}=Y_{\mathrm{o}}+\Delta Y_{\mathrm{o}}, \\Z_{\mathrm{o}}^{\mathrm{a}}=Z_{\mathrm{o}}+\Delta Z_{\mathrm{o}},\end{array} \right\} $

在使用球杆仪测量时，球杆仪末端在半球球面上移动，设球杆仪理论长度为L_t，实际测量长度为L_a，理论长度如下：

(27)$ {L_{\mathrm{t}}} = \left[ {{{({X_{\mathrm{t}}} - X_{{\mathrm{o}}} ^{{\mathrm{a}}} )}^2}+{{({Y_{{\mathrm{t}}} } - Y_{{\mathrm{o}}} ^{{\mathrm{a}}} )}^2}+{{({Z_{{\mathrm{t}}} } - Z_{{\mathrm{o}}} ^{{\mathrm{a}}} )}^2}}\right]^{1/2} . $

通过测量可得到f个包含误差参数的理论长度， L_t包含表4中的12个参数误差以及3个球心偏差量(ΔX_o, ΔY_o, ΔZ_o)共15个误差参数.

采用圆弧插补算法确定空间轨迹初始点位，随机选取不共线的3点(P₁, P₂, P₃)，其初始值为P₁(x₁, y₁, z₁)，P₂(x₂, y₂, z₂)，P₃(x₃, y₃, z₃)，可以确定半径为r、圆心为o的圆弧，P_a、P_b的垂直平分面方程如下：

(28)$ \begin{split} ({x_{b} } - &{x_a})x+({y_b} - {y_{a} })y+({z_{b} } - {z_{a} })z+ \\ & \frac{1}{2}(x_{b} ^2+y_{b} ^2+z_{b} ^2 - x_{a} ^2 - y_{a} ^2 - z_{a} ^2) = 0. \end{split} $

由式(28)可以求得P₁与P₂的垂直平分面方程和P₂与P₃的垂直平分面方程. 平面P₁P₂P₃外接圆方程如下：

(29)$ \begin{split} (({y_1} -& {y_3})({z_2} - {z_3}) - ({y_2} - {y_3})({z_1} - {z_3}))x+ \\ & (({x_2} - {x_3})({z_1} - {z_3}) - ({x_1} - {x_3})({z_2} - {z_3}))y+ \\ & (({x_1} - {x_3})({y_2} - {y_3}) - ({x_2} - {x_3})({y_1} - {y_3})z+ \\ & ( - ({A_{31}}{x_3}+{B_{32}}{y_3}+{C_{33}}{z_3})) = 0. \end{split} $

由式(28)、(29)可求得外接圆圆心o(x₀, y₀, z₀)，外接圆半径可由半径计算公式求得.

圆弧插补角度计算公式如下：

(30)$ {\theta _i} = \left\{ {\begin{array}{*{20}{l}} \arctan \;2({y_3},{x_3})+2{\text{π}} ,&{y_3} < 0; \\ \arctan\; 2({y_3},{x_3}),& {y_3} \geqslant 0. \end{array}} \right. $

根据求得的圆心坐标、半径R以及插补角度可求得圆弧上点的坐标：

(31)$\left. \begin{split} {x_i} =& {x_0}+R \cos\;{\theta _i}, \\ {y_i} =& {y_0}+R \sin \;{\theta _i}, \\ {z_i} =& (1 - {t_i}) {z_1}+{t_i} {z_3}. \end{split}\right\} $

式中：t_i为第i个点的线性插值参数.

当待辨识的误差参数趋于准确值时，球杆仪理论长度L_t趋近于测量值L_a，构建误差参数辨识的目标函数H，使目标函数值等于理论长度与实测长度的差且趋近于无穷小：

(32)$ H=\dfrac{1}{f} {\displaystyle\sum_{n=1}^f\left|L_{\mathrm{t}}-L_{\mathrm{a}}\right|} . $

为了避免局部最优，加速误差参数辨识，位置更新公式如下：

(33)$ \begin{split}X_i^{(t+1)}=& X_i^{(t)}+s \left(\phi_i^{(t)}-X_i^{(t)}\right)+ \\& {\rm{e}}^{{bl}} \cos \;(2 {\text{π}} l) \times \left(X_1^{(t)}-X_i^{(t)}\right) .\end{split} $

式中：$X_i^{(t)} $为第i个飞蛾在t次迭代时的位置，$\phi_i^{(t)} $第i个飞蛾感知到的最优适应度值，$X_1^{(t)} $为头领飞蛾位置，s为步长，b为螺旋搜索常数，l为[−1.0, 1.0]的随机数. 算法流程如图6所示.

使用Compute Unified Device Architecture (CUDA)对辨识算法进行加速，具体过程如下.

1)分离目标函数中计算(X_t−$X_{\mathrm{o}}^{\mathrm{a}} $)²、(Y_t−$Y_{\mathrm{o}}^{\mathrm{a}} $)²、(Z_t−$Z_{\mathrm{o}}^{\mathrm{a}} $)²的部分，将其传入GPU中.

2) 随机生成误差赋给步骤1)中的公式.

3) 计算结果导出至CPU，计算目标函数值，当接近0.02时收敛速度显著下降，表明逼近最优解. 函数值小于0.02时迭代停止.

对于采样轨迹上的f个采样点进行多次测量求得各点的标准不确定度：

(34)$ u = \sqrt {\frac{1}{{n(n - 1)}} {{\sum\limits_{i = 1}^n {({a_i} - \bar a)} }}} . $

式中：$ \bar a $为采样点测量值的均值. 计算均值的不确定度区间为[−b，b]，b可以表达为

(35)$ b = \bar a \pm z u. $

工业机器人定位误差符合正态分布，在95%的置信水平下，z=1.96.

工业机器人重复定位误差及随机误差导致测量值具有不确定度，须对测量值进行动态修正，将测量值修正部分与辨识算法进行结合，获得修正后的误差辨识值，嵌套算法流程图如图7所示.

以采样点p为例，其不确定度范围为[−b, b]，p点测量值为a，评估a并对参数c动态调节. 当a位于[−b,−b/2]或[b/2,b]时，c取2；当a位于(−b/2, b/4]或(b/4, b/2)时，c取4；当|a|小于b/4时，不修正. 修正后的测量值a_n表达式如下：

(36)$ {a_{\mathrm{n}}} = a - {{\mathrm{sgn}}} \;(a - \bar a) \times \frac{b}{c} \times q. $

式中：a为初始测量值，$ \bar a $为测量平均值，b为不确定度区间，c为不确定度修正系数，q为取值范围为(0,1.0]的随机数.

嵌套补偿算法首先进行种群初始化，并评估待辨识参数的不确定度范围. 根据式(36)修正待辨识参数的不确定度. 选择头领飞蛾并计算虑及不确定度的适应度. 最后，判断迭代过程是否满足收敛条件，以及修正后的待辨识参数是否符合修正判断条件，即目标函数值小于0.02且修正后的待辨识参数a_n的绝对值小于b/4.

实验中的ABB-IRB-120型机器人重复定位精度为±0.01 mm，额定负载为3 kg. 雷尼绍QC20-W球杆仪的测量半径为100 mm，理论球心坐标为：X_o =350，Y_o =0，Z_o =394，测量轨迹以5°为间隔均匀取点. 每条轨迹间隔取5个点用于轨迹绘制，如图8所示.

实验中机器人末端速度为30 mm/s，采样率为2545.6 mm/min，每条轨迹测量10次. 按式(35)计算采样点的不确定度区间. 采样点对应的误差ε与不确定度区间如图9所示. 可以看出，在不确定性较弱的采样点处，不确定度区间长度约为0.05 mm，在不确定度较强位置处，不确定度区间长度约为0.50 mm.

以最后一次测量结果作为补偿实验误差测量结果. 基于测量数据构建目标函数，分别用传统MFO、改进MFO和鲸鱼优化算法(WOA)进行10次迭代求解，比较总用时和迭代次数的平均值. 迭代条件为全局最优值小于0.02，种群数量分别为500、1500、3000，结果如图10所示. 其中，$T_{\mathrm{s}} $为总用时，$C_{\mathrm{s}} $为迭代次数. 可以看出，改进MFO算法在种群数量为500、1500、3000时，总计算时长较WOA算法平均降低62.03%，较传统MFO算法降低42.46%；迭代次数较WOA算法平均提升70.45%，较传统MFO算法提升42.33%.

对不考虑和考虑不确定度的误差参数进行辨识，通过优化算法求解，修正POE运动学模型参数，逆解轨迹点位得到修正关节角度，辨识及修正后的参数如表6所示. 可以看出，在不考虑不确定度精度的影响时，连杆长度波动为−0.263~0.269 mm，角度波动为−0.166°~0.181°；考虑不确定度后，连杆长度波动缩小至−0.064~0.158 mm，角度波动为−0.211°~0.044°. 机器人理论运行轨迹与补偿前、补偿后、虑及不确定度补偿的轨迹，如图11所示.

如表7所示为补偿效果对比. 其中，${{e}}_{\mathrm{p}} $、${{e}}_{\mathrm{j}} $分别为绝对平均误差、均方根误差，$\varphi $为补偿效率. 不考虑不确定度时绝对平均误差相较于补偿前减少49.83%，均方根误差相较于补偿前减少45.65%. 考虑不确定度时绝对平均误差相较于补偿前减少53.47%，相较于不考虑不确定度时减少7.31%. 考虑不确定度时均方根误差相较于补偿前减少48.89%，相较于不考虑不确定度时减少7.09%.

为了验证所提出的精度校准方法在工业应用中的性能与效果，在ABB IRB 120工业机器人上进行实验，选用Ф2 mm HRC45 QXD球头铣刀，末端执行器转速为3 000 r/min，以铝合金材质棒料为加工对象. 加工叶轮工序如下. 1) 叶轮粗加工. 用球形铣刀铣削叶片间流道，从轮盖侧至轮毂侧逐步进行. 2) 叶轮精加工. 精修叶片表面及流道，确保型线精确、表面光滑，以减少阻力，提升效率. 3) 对叶轮进行外观修整，去除毛刺、锐角，确保叶轮表面质量.

使用三坐标测量机对补偿前后的叶轮半径、叶轮高度、叶片高度、上下平面平行度与平面度进行检测，实验现场及工件对比如图12所示. 检测结果显示，补偿前叶轮半径、叶片高度、叶轮高度、平行度、平面度各自的误差分别为0.086、0.109、0.093、0.015、0.008 mm；补偿后分别为0.054、0.069、0.066、0.009、0.005 mm. 尺寸误差平均减小32.3%，形位误差平均减小38.9%.

针对提高工业机器人定位精度的问题，提出误差参数校正方法，并进行实验验证，本研究的主要结论如下.

(1)通过工业机器人关节与空间定位误差相似特性分析，设计了基于球杆仪的空间去冗余测量轨迹，实现了采样点的低相似性及去冗余测量.

(2) 建立了基于包围策略与猎物搜索机制的改进MFO优化逆解、辨识方法，改进的MFO算法相较于传统的飞蛾火焰优化算法与鲸鱼优化算法，收敛速率分别提升了42.46%、62.03%.

(3)构建了补偿参数嵌套算法并实现了辨识参数的不确定度修正. 实验结果表明，不考虑不确定度补偿，绝对平均误差相较于补偿前减少49.83%，均方根误差相较于补偿前减少45.65%；考虑不确定度补偿，绝对平均误差相较于补偿前减少53.47%，均方根误差相较于补偿前减少48.89%. 补偿后加工的叶轮尺寸误差相较于未补偿时加工的叶轮平均减小32.3%，形位误差相较于未补偿时加工的叶轮平均减小38.9%.

为了实现在线精度校准，后续将致力基于机器视觉的在线补偿策略研究，以实现工业机器人在复杂动态环境中定位精度的提升.