基于关节几何误差的机器人绝对位置精度标定算法

Absolute Position Accuracy Calibration Algorithm for Robots Based on Joint Geometric Error

摘要

提出一种引入关节几何误差参数的工业机器人运动学模型和误差参数标定算法.首先，在DH（Denavit-Hartenberg）模型基础上，为每个关节引入6个几何误差参数，建立更为完善的标定误差模型，并实现该模型的正逆运动学求解算法.然后，建立包含关节误差、基坐标误差和工具坐标误差共45个参数的微分运动学雅可比矩阵，并采用小样本测试集实现基于迭代算法的误差参数求解.最后，在新松SR10C机器人上使用激光跟踪仪完成实验验证，将标定得到的误差参数补偿到模型中.结果表明，经过标定补偿后的机器人最大位置误差减小约80%，平均位置误差减小约80%，误差方差减小约97%，证明该方法可大幅提升工业机器人绝对位置精度和确定性.

Abstract

An industrial robot kinematic model with joint geometric error parameters and a calibration algorithm is proposed. Firstly， based on the DH model， six geometric error parameters are introduced for each joint to establish a more comprehensive error calibration model. The solutions of forward and inverse kinematic for the model are realized. Then， a differential kinematic Jacobian matrix containing 45 parameters， including joint errors， base coordinate errors， and tool coordinate errors is established. An iterative algorithm based on a small sample test set is used to solve the error parameters. Finally， experimental verification is carried out using a laser tracker on the SIASUN SR10C robot. The calibrated error parameters are then compensated into the model. Results show that， after calibration compensation， the maximum position error of the robot decreases by approximately 80%， the average position error decreases by approximately 80%， and the error variance decreases by approximately 97%， demonstrating that this method can significantly improve the absolute position accuracy and determinacy of industrial robots.

Graphical abstract

关键词

工业机器人 / 绝对位置精度 / 运动学模型 / 几何误差 / 标定

Key words

industrial robot / absolute position accuracy / kinematic model / geometric error / calibration

工业机器人的精度主要包括重复位置精度和绝对位置精度.随着高精度编码器和伺服技术的不断发展，国内外工业机器人的重复位置精度均已达到较高水平，能够充分满足示教再现任务的要求.但在离线编程、精密加工、精密检测等应用场景中，要求机器人实现精确点位控制和轨迹控制，则需要对机器人进行标定以提升绝对位置精度.工业机器人绝对精度标定算法研究的目的是通过分析、建模和校准机器人的误差源，来提高机器人的绝对位置精度.根据待辨识参数建模层级的不同，机器人参数校准研究分为关节层、运动学层和非几何参数层^［1］.关节层主要通过校准减速比和零位偏置来减小各关节转角实测值与理论值之间的偏差^［2-3］；运动学层主要校准机器人的结构参数与位移偏差^［4-5］；非几何参数层主要校准由于关节变形、齿轮间隙、运动控制等因素导致的定位与定向偏差^［6-7］.针对运动学层级的误差校准，需要选择合适的运动学模型，使其精确地表达各关节和连杆间的转换关系，进而构建相应的误差模型.目前构建机器人运动学模型主要采用DH模型^［8］，模型只使用4个参数来表示单关节的空间位置坐标.改进DH（MDH）^［9］模型通过在平行的奇异关节坐标系中加入绕y轴旋转的参数，解决了DH模型在某些姿态下不满足连续性的问题.基于旋量理论的指数积模型^［10］，每个关节的坐标系使用6个参数描述，避免了奇异性问题，但参数存在冗余.机器人运动学参数辨识是依据测量得到的末端位姿信息与自身关节角度信息，使用优化算法得到最优机器人运动学结构参数.当机器人的结构参数偏差较小时，可使用线性化方法构建运动学误差模型来描述机器人运动学参数偏差与位姿误差间的近似关系，并用各类迭代方法得到结构参数偏差的最小二乘解^［11］.另外，基于扩展卡尔曼滤波法、人工神经网络等^［12-13］经典最优化求解方法也应用于优化求解问题中.

本文基于标准DH模型，提出了一种引入关节几何误差的运动学模型，并实现了基于该运动学模型的正逆运动学求解算法.之后设计了同时包含关节误差、工具坐标误差和测量坐标误差的微分运动学雅可比矩阵，构建了绝对位置精度测量模型，并采用最优测试集和最小二乘法计算误差参数.最后通过工业机器人和激光跟踪仪验证了该算法的有效性.

1 引入几何误差的机器人运动学模型

机器人坐标系设定如图1所示，{}_{}{}^{\mathrm{0}}F为机器人基坐标系，{}_{}{}^{\mathrm{m}}F为外部测量坐标系，机器人每个轴的固连坐标系为{}_{}{}^{i}F，均为x轴向前、z轴向上，各轴正方向均遵照重合轴方向的右手定则.由于机器人关节及连杆的加工、组装问题，实际机器人关节参数与理论不同，会导致这些固连坐标系发生偏移和旋转.此时若采用标准DH模型计算机器人末端位姿，会导致机器人绝对定位精度降低.

为使机器人具有更高的绝对位置精度，可在运动学模型中对上述误差进行补偿.在理论DH模型基础上，实际坐标系和理论坐标系间的偏移用\mathit{O}={\left[\begin{array}{ccc}{o}_x& o_y& o_z\end{array}\right]}^{\mathrm{T}}表示，旋转关系用俯仰角、航向角和横滚角表示，与机器人轴平行方向旋转量记为\gamma，其他两个方向依次记为\alpha和\beta，如图2所示.

1.1 正向运动学求解

引入几何误差的各关节坐标系相对于基坐标系的姿态矩阵：

\left.\begin{array}{c}{}_{\mathrm{1}}{}^{\mathrm{0}}\widehat{\mathit{R}}={\mathit{R}}_x\left({\alpha }_{\mathrm{1}}\right){\mathit{R}}_y\left({\beta }_{\mathrm{1}}\right){\mathit{R}}_z\left({\theta }_{\mathrm{1}}-{\gamma }_{\mathrm{1}}\right),\\ {}_{\mathrm{2}}{}^{\mathrm{1}}\widehat{\mathit{R}}={\mathit{R}}_x\left({\alpha }_{\mathrm{2}}\right){\mathit{R}}_z\left({\beta }_{\mathrm{2}}\right){\mathit{R}}_y\left({\theta }_{\mathrm{2}}-{\gamma }_{\mathrm{2}}\right),\\ {}_{\mathrm{3}}{}^{\mathrm{2}}\widehat{\mathit{R}}={\mathit{R}}_x\left({\alpha }_{\mathrm{3}}\right){\mathit{R}}_z\left({\beta }_{\mathrm{3}}\right){\mathit{R}}_y\left({\theta }_{\mathrm{3}}-{\gamma }_{\mathrm{3}}\right),\\ \begin{array}{l}{}_{\mathrm{4}}{}^{\mathrm{3}}\widehat{\mathit{R}}={\mathit{R}}_y\left({\alpha }_{\mathrm{4}}\right){\mathit{R}}_z\left({\beta }_{\mathrm{4}}\right){\mathit{R}}_x\left({\theta }_{\mathrm{4}}-{\gamma }_{\mathrm{4}}\right),\\ {}_{\mathrm{5}}{}^{\mathrm{4}}\widehat{\mathit{R}}={\mathit{R}}_x\left({\alpha }_{\mathrm{5}}\right){\mathit{R}}_z\left({\beta }_{\mathrm{5}}\right){\mathit{R}}_y\left({\theta }_{\mathrm{5}}-{\gamma }_{\mathrm{5}}\right),\\ {}_{\mathrm{6}}{}^{\mathrm{5}}\widehat{\mathit{R}}={\mathit{R}}_y\left({\alpha }_{\mathrm{6}}\right){\mathit{R}}_z\left({\beta }_{\mathrm{6}}\right){\mathit{R}}_x\left({\theta }_{\mathrm{6}}-{\gamma }_{\mathrm{6}}\right).\end{array}\end{array}\right\}(1)

式中：{\mathit{R}}_x，{\mathit{R}}_y和{\mathit{R}}_z分别为绕坐标轴x，y和z的旋转算子；{\alpha }_i，{\beta }_i，{\gamma }_i为关节i旋转误差角度；{\theta }_i为关节i旋转轴相对零位的角度值；{}_{i+\mathrm{1}}^i\widehat{\mathit{R}}为引入关节误差后相邻两关节i和i+1之间的旋转矩阵.

从而可得机器人末端姿态：

{}_{\mathrm{6}}{}^{\mathrm{0}}\widehat{\mathit{R}}={}_{\mathrm{1}}{}^{\mathrm{0}}\widehat{\mathit{R}}{}_{\mathrm{2}}{}^{\mathrm{1}}\widehat{\mathit{R}}{}_{\mathrm{3}}{}^{\mathrm{2}}\widehat{\mathit{R}}{}_{\mathrm{4}}{}^{\mathrm{3}}\widehat{\mathit{R}}{}_{\mathrm{5}}{}^{\mathrm{4}}\widehat{\mathit{R}}{}_{\mathrm{6}}{}^{\mathrm{5}}\widehat{\mathit{R}}.(2)

机器人工具中心点（tool center point，TCP）的位置计算：

\mathit{P}={\mathit{O}}_{\mathrm{1}}+{\mathit{L}}_{\mathrm{1}}+{}_{\mathrm{1}}{}^{\mathrm{0}}\widehat{\mathit{R}}\left({\mathit{O}}_{\mathrm{2}}+{\mathit{L}}_{\mathrm{2}}\right)+{}_{\mathrm{2}}{}^{\mathrm{0}}\widehat{\mathit{R}}\left({\mathit{O}}_{\mathrm{3}}+{\mathit{L}}_{\mathrm{3}}\right)+

{}_{\mathrm{3}}{}^{\mathrm{0}}\widehat{\mathit{R}}\left({\mathit{O}}_{\mathrm{4}}+{\mathit{L}}_{\mathrm{4}}\right)+{}_{\mathrm{4}}{}^{\mathrm{0}}\widehat{\mathit{R}}\left({\mathit{O}}_{\mathrm{5}}+{\mathit{L}}_{\mathrm{5}}\right)+{}_{\mathrm{5}}{}^{\mathrm{0}}\widehat{\mathit{R}}\left({\mathit{O}}_{\mathrm{6}}+{\mathit{L}}_{\mathrm{6}}\right).(3)

1.2 逆向运动学求解

引入几何误差后，机器人模型不再满足Pieper原则^［14］，因此该模型不存在解析解，需要通过数值法求解模型逆解.基于微分运动学原理，由雅可比矩阵获得关节变化量与机器人TCP变化量之间的关系，从而可通过迭代方式得到满足误差限值的机器人运动学逆解.雅可比矩阵\mathit{J}(\mathit{q})描述了机器人操作空间速度\dot{\mathit{x}}与关节空间速度\dot{\mathit{q}}之间的变换：

\dot{\mathit{x}}=\mathit{J}(\mathit{q})\dot{\mathit{q}}.(4)

式中，\mathit{q}=[{\theta }_{\mathrm{1}},{\theta }_{\mathrm{2}},...,{\theta }_{\mathrm{6}}{]}^{\mathrm{T}}为关节角度向量.雅可比矩阵可分为线速度（ v ）部分和角速度（ ω ）部分，迭代求解公式为

\left.\begin{array}{c}{}_{}{}^{i+\mathrm{1}}{\mathit{\omega }}_{i+\mathrm{1}}={}_i{}^{i+\mathrm{1}}\mathit{R}{}_{}{}^i\mathit{\omega }{\mathit{\omega }}_i+{\theta }_{i+\mathrm{1}}^{\mathrm{\text{'}}}{}_{}{}^{i+\mathrm{1}}\widehat{\mathit{Z}}{\widehat{\mathit{Z}}}_{i+\mathrm{1}},\\ {}_{}{}^{i+\mathrm{1}}{\mathit{v}}_{i+\mathrm{1}}={}_i{}^{i+\mathrm{1}}\mathit{R}({}_{}{}^i\mathit{v}{\mathit{v}}_i+{}_{}{}^i\mathit{\omega }{\mathit{\omega }}_i\times {}_{}{}^i\mathit{P}{\mathit{P}}_{i+\mathrm{1}}).\end{array}\right\}(5)

式中：{{}_{}{}^i}_i和{{}_{}{}^i}_i为连杆i的角速度和线速度相对于坐标系｛i｝的分量；{{}_{}{}^{i+\mathrm{1}}}_{i+\mathrm{1}}为关节i+1在坐标系｛i+1｝中的轴向量；{{}_{}{}^i}_{i+\mathrm{1}}为关节i轴向量与关节i+1轴向量间的距离.

由于机器人末端误差约90%来源于关节角度偏差 ^［15］，因此，对于引入误差的机器人模型，仅考虑在关节轴线上的旋转误差\mathrm{\Delta }\gamma，则关节变量可重新表示为

{\mathit{q}}_{\gamma }=\mathit{q}+\mathrm{\Delta }\gamma.(6)

将式（6）代入式（5），可得引入几何误差的雅可比矩阵\mathit{J}({\mathit{q}}_{\gamma }).由雅可比矩阵获得关节变化量与机器人TCP变化量之间的关系，通过迭代方式得到满足误差限值的机器人运动学逆解，流程如图3所示.

在迭代计算过程中，待求解逆运动学的机器人末端位置和姿态分别记为\mathrm{p}\mathrm{o}{\mathrm{s}}_{\mathrm{e}\mathrm{n}\mathrm{d}}和\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}{\mathrm{t}}_{\mathrm{e}\mathrm{n}\mathrm{d}}；当前迭代步对应的逆运动学解记为{\mathit{q}}_{\mathrm{c}\mathrm{u}\mathrm{r}}，根据正运动学计算得到的末端位置和姿态分别记为\mathrm{p}\mathrm{o}{\mathrm{s}}_{\mathrm{c}\mathrm{u}\mathrm{r}}和\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}{\mathrm{t}}_{\mathrm{c}\mathrm{u}\mathrm{r}}.上述定义确定后，则按照如下步骤计算给定末端位姿的逆运动学问题.

1） 依据DH模型解析求解方法，由\mathrm{p}\mathrm{o}{\mathrm{s}}_{\mathrm{e}\mathrm{n}\mathrm{d}}和\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}{\mathrm{t}}_{\mathrm{e}\mathrm{n}\mathrm{d}}得到标准模型反解{\mathit{q}}_{\mathrm{c}\mathrm{u}\mathrm{r}}^［14］，以此作为迭代求解的初始值；

2） 由{\mathit{q}}_{\mathrm{c}\mathrm{u}\mathrm{r}}按照式（2）和式（3）得到末端位置\mathrm{p}\mathrm{o}{\mathrm{s}}_{\mathrm{c}\mathrm{u}\mathrm{r}}和姿态\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}{\mathrm{t}}_{\mathrm{c}\mathrm{u}\mathrm{r}}，由此可得位置差值向量{\mathit{\delta }}_{\mathrm{p}\mathrm{o}\mathrm{s}}=\mathrm{p}\mathrm{o}{\mathrm{s}}_{\mathrm{e}\mathrm{n}\mathrm{d}}-\mathrm{p}\mathrm{o}{\mathrm{s}}_{\mathrm{c}\mathrm{u}\mathrm{r}}，姿态差值向量{\mathit{\delta }}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}}=\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}{\mathrm{t}}_{\mathrm{e}\mathrm{n}\mathrm{d}}-\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}{\mathrm{t}}_{\mathrm{c}\mathrm{u}\mathrm{r}}，并判断当前位置差值向量的二次范数是否小于迭代精度\epsilon.如果||{\mathit{\delta }}_{\mathrm{p}\mathrm{o}\mathrm{s}}||\le \epsilon，则迭代结束，当前关节值即为逆运动学解{\mathit{q}}_{\mathrm{e}\mathrm{n}\mathrm{d}}；如果||{\mathit{\delta }}_{\mathrm{p}\mathrm{o}\mathrm{s}}||>\epsilon，则进行步骤3）；

3） 由当前逆运动学解{\mathit{q}}_{\mathrm{c}\mathrm{u}\mathrm{r}}按照式（5）计算对应的雅可比矩阵\mathit{J}，并根据式（4）计算机器人末端微分变量{\mathit{\delta }}_{\mathrm{q}}={\mathit{J}}^{-\mathrm{1}}{\mathit{\delta }}_{\mathrm{p}\mathrm{o}\mathrm{s},\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}}，从而得到新的迭代关节值{\mathit{q}}_{\mathrm{c}\mathrm{u}\mathrm{r}}(k+\mathrm{1})={\mathit{q}}_{\mathrm{c}\mathrm{u}\mathrm{r}}(k)+{\mathit{\delta }}_{\mathrm{q}}.

重复步骤2），3），使位置向量差值的二次范数小于迭代精度\epsilon，得到机器人末端姿态\mathrm{p}\mathrm{o}{\mathrm{s}}_{\mathrm{e}\mathrm{n}\mathrm{d}}和\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{n}{\mathrm{t}}_{\mathrm{e}\mathrm{n}\mathrm{d}}的逆运动学解{\mathit{q}}_{\mathrm{e}\mathrm{n}\mathrm{d}}.如果迭代次数超出设定限值，则表明算法无法在预设迭代次数内收敛，即没有找到满足要求的逆解.

2 关节几何误差标定算法

2.1 测试集的选取

末端绝对位置误差服从正态分布规律，因此可以采用最小测量集N_{\mathrm{m}\mathrm{i}\mathrm{n}}进行误差标定：

N_{\mathrm{m}\mathrm{i}\mathrm{n}}={\left(\frac{z_{\alpha /\mathrm{2}}\sigma }{E}\right)}^{\mathrm{2}}.(7)

式中：z_{\alpha /\mathrm{2}}为误差分布置信区间；E为最大均差；{\sigma }^{}为方差.一般选取E不小于机器人的重复定位精度，当测试集数量不小于N_{\mathrm{m}\mathrm{i}\mathrm{n}}时，测量误差不大于E的置信度为(\mathrm{1}-\alpha )\times \mathrm{100}\mathrm{\%}.

测试集的选取采用拟随机法^［16］：基于蒙特卡罗法随机生成一组测量数据集，根据测量坐标系与基坐标系相对关系和可观测指数，计算最优的标定姿态值，并通过运动学逆解计算相应的关节角度值，最终，将测试集中的关节角度值作为控制命令发送给机器人控制器.选取O_{\mathrm{1}}可观测指数^［17］对测试集的预设位姿进行评估：

O_{\mathrm{1}}=\frac{({\sigma }_{\mathrm{1}}{\sigma }_{\mathrm{2}}...{\sigma }_m{)}^{\frac{\mathrm{1}}{m}}}{\sqrt[]{N}}.(8)

式中：{\sigma }_i为雅可比矩阵经奇异值分解（SVD）后对角线矩阵上第i个元素；m为待标定参数总数；N为测量集总数.

通过式（7）、式（8）的计算，既减少了测试样本数以提高标定效率，又保证了测量集覆盖整个机器人的工作空间和测试点观测性最优.

2.2 基坐标变换初始值计算

由图1可知，测量数据基于测量坐标系{}_{}{}^{\mathrm{m}}F，机器人模型基于基坐标系{}_{}{}^{\mathrm{0}}F，在标定计算之前需要获取一些已知位置或姿态的数据点，然后根据这些数据点计算得到基坐标系与测量坐标系之间的变换矩阵.该变换矩阵可以用来将测量坐标系下的位置和姿态值转换到基坐标系，实现测量位姿值与模型位姿值的匹配.选取k(k\ge \mathrm{3})个观测值构成观测数据集{}_{}{}^{\mathrm{m}}\mathit{p}=\{{{}_{}{}^{\mathrm{m}}\mathit{p}}_{{}_{\mathrm{1}}},{{}_{}{}^{\mathrm{m}}\mathit{p}}_{{}_{\mathrm{2}}},...,{{}_{}{}^{\mathrm{m}}\mathit{p}}_k\}，对应模型理论值构成理论数据集{}_{}{}^{\mathrm{0}}\mathit{p}=\{{{}_{}{}^{\mathrm{0}}\mathit{p}}_{{}_{\mathrm{1}}},{{}_{}{}^{\mathrm{0}}\mathit{p}}_{{}_{\mathrm{2}}},...,{{}_{}{}^{\mathrm{0}}\mathit{p}}_k\}，根据中心计算式（9）可得2个数据集的中心分别为{}_{}{}^{\mathrm{m}}\mathit{C}和{}_{}{}^{\mathrm{0}}\mathit{C}.

\mathit{C}=\frac{\mathrm{1}}{k}{\displaystyle \sum _{i=\mathrm{1}}^k}{\mathit{p}}_i.(9)

测试点到各自中心的向量为{{}_{}{}^{\mathrm{m}}}_{\mathrm{c}}=\left\{{}_{}{}^{\mathrm{m}}\mathit{p}-{}_{}{}^{\mathrm{m}}\mathit{C}\right\}， {{}_{}{}^{\mathrm{0}}}_{\mathrm{c}}=\left\{{}_{}{}^{\mathrm{0}}\mathit{p}-{}_{}{}^{\mathrm{0}}\mathit{C}\right\}，均为\mathrm{3}\times k向量，则观测坐标系到基坐标系的转换矩阵\mathit{R}可表示为

\mathit{R}=({}_{}{}^{\mathrm{0}}\mathit{p}{)}^{\mathrm{T}}{}_{}{}^{\mathrm{m}}\mathit{p}.(10)

由于存在测量误差和截断误差，齐次矩阵的旋转矩阵部分不能满足单位正交矩阵的性质.为了保证坐标系旋转后的精度，本文采用四元数法表示坐标系旋转^［18］，并选取旋转误差最小的轴作为四元数旋转轴，从而得到满足不等式（11）的旋转轴i(i=\mathrm{1,2},\mathrm{3,4}).

\sqrt[]{||(\mathit{R}(i){}_{}{}^{\mathrm{m}}\mathit{p}-{}_{}{}^{\mathrm{0}}\mathit{p})||}\le \eta \frac{S+D}{\mathrm{2}}.(11)

式中：\eta为阈值比例系数；S为{}_{}{}^{\mathrm{m}}p各点连线所构成的图形面积；D为{}_{}{}^{\mathrm{m}}\mathit{p}到{}_{}{}^{\mathrm{0}}\mathit{p}的位移大小，如式（12）、式（13）所示.

S=\frac{\mathrm{1}}{k}\sqrt[]{{\displaystyle \sum _{i=\mathrm{1}}^{k-\mathrm{1}}}||{}_{}{}^{\mathrm{m}}\mathit{p}{\mathit{p}}_i-{}_{}{}^{\mathrm{m}}\mathit{p}{\mathit{p}}_{i+\mathrm{1}}||+||{}_{}{}^{\mathrm{m}}\mathit{p}{\mathit{p}}_k-{}_{}{}^{\mathrm{m}}\mathit{p}{\mathit{p}}_{\mathrm{1}}||}，(12)

D=\frac{\mathrm{1}}{k}\sqrt[]{{\displaystyle \sum _{i=\mathrm{1}}^k}||{}_{}{}^{\mathrm{m}}\mathit{p}{\mathit{p}}_i-{}_{}{}^{\mathrm{0}}\mathit{p}{\mathit{p}}_i||}.(13)

最终，获得观测坐标系到基坐标系的旋转及平移向量如式（14）所示，将其作为基坐标误差的初始值{\mathit{T}}_{\mathrm{0}}.

\left.\begin{array}{l}\mathit{r}=\mathrm{Q}\mathrm{U}\mathrm{A}\mathrm{T}(\mathit{R},i)，\\ \mathit{p}={}_{}{}^{\mathrm{0}}\mathit{C}-\mathit{r}{}_{}{}^{\mathrm{1}}\mathit{C}.\end{array}\right\}(14)

式中，QUAT为基于四元数的旋转轴向量计算公式.

2.3 误差参数计算

机器人微分运动学描述了关节结构参数微小变化与末端笛卡尔坐标变化的关系，其映射关系可用雅可比矩阵描述：

\mathrm{\Delta }\mathit{p}=\mathit{J}(\mathit{l})\mathrm{\Delta }\mathit{l}.(15)

式中：\mathit{l}为引入关节几何误差的机器人结构参数；包括从测量坐标系到机器人基坐标系的坐标变换参数{\mathit{l}}_{\mathrm{b}}；机器人各关节结构参数{\mathit{l}}_{\mathrm{r}}和工具坐标系变换参数{\mathit{l}}_{\mathrm{t}}. {\mathit{l}}_{\mathrm{b}}=[o_x,o_y,o_z{,\alpha ,\beta ,\gamma ]}^{\mathrm{T}}共6个参数，关节i的参数为{\mathit{l}}_{\mathrm{r}}(i)=[o_x,o_y,o_z{,\alpha ,\beta ,\gamma ]}^{\mathrm{T}}，机器人关节参数共36个，{\mathit{l}}_{\mathrm{t}}=[o_x,o_y,o_z{]}^{\mathrm{T}}共3个参数，因此，标定需要校准的参数共有45个； \mathrm{\Delta }\mathit{p}表示机器人末端由于参数变化产生的位置变化量.假设标定实验中有N个观测点，则雅可比矩阵\mathit{J}为\mathrm{3}N\times \mathrm{45}的矩阵：

\left.\begin{array}{c}\mathit{J}[\mathrm{3}n+\mathrm{1}][k]=\mathrm{\Delta }x(k)/\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p},\\ \mathit{J}[\mathrm{3}n+\mathrm{2}][k]=\mathrm{\Delta }y(k)/\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p},\\ \mathit{J}[\mathrm{3}n+\mathrm{3}][k]=\mathrm{\Delta }z(k)/\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}.\end{array}\right\}(16)

式中：n=\mathrm{0,2},...,N-\mathrm{1}；k=\mathrm{1,2},...,\mathrm{45}；\mathrm{\Delta }x(k)，\mathrm{\Delta }y(k)，\mathrm{\Delta }z(k)为第k个参数引入微小变化step时产生的末端位移偏差.

本文采用线性最小二乘法求解机器人结构误差参数.使用雅可比矩阵伪逆，根据末端坐标变化值计算求得参数变化值：

\mathrm{\Delta }\mathit{l}=({\mathit{J}}^{\mathrm{T}}{\mathit{J})}^{-\mathrm{1}}{\mathit{J}}^{\mathrm{T}}\mathrm{\Delta }\mathit{P}.(17)

式中：\mathrm{\Delta }\mathit{P}为N个测试点由于参数引入微小变化step所产生的偏差值\mathrm{\Delta }\mathit{p}构成的向量；用\mathrm{\Delta }\mathit{l}更新机器人结构参数\mathit{l}，并使用引入几何误差的机器人正运动学模型计算N个观测点的末端坐标{}_{\mathrm{b}}\mathit{P}，当||{}_{\mathrm{b}}{}^{\mathrm{0}}\mathit{P}-{}_{\mathrm{m}}{}^{\mathrm{0}}\mathit{P}||小于收敛阈值时，得到最终的关节参数误差；否则，继续对各参数增加step，直到收敛或超过迭代次数限值.在迭代运算中，{\mathit{l}}_{\mathrm{b}}的初始值由式（14）计算得到，{\mathit{l}}_{\mathrm{r}}的各元素初始值均为零，{\mathit{l}}_{\mathrm{t}}的初始值为理论值，标定流程如图4所示.

3 实验和分析

3.1 实验平台

本文实验平台如图5所示，机器人型号为新松SR10C，测量仪器为FARO激光测试仪，机器人控制器采用嵌入式工业主板，运动学计算及标定参数计算均在控制器中完成，具体参数如表1所示.实验中，首先获得测量数据，通过标定算法得到各误差参数，并将结果补偿到运动学模型中；比较补偿前后的位置理论值与测量值，验证算法性能.

3.2 绝对精度标定实验

由式（7）、式（8）获取70组测试姿态及其相应的关节值，并通过激光跟踪仪测量相应的观测数据，如图6～图7所示，数据涵盖了各关节的运动范围和机器人末端的工作空间.

按照图4所示的标定计算流程，通过观测值和对应的关节角度值，得到关节几何误差参数、基坐标系和工具坐标系如表2所示.

将标定得到的误差参数代入式（2）、式（3），可根据测试集中的关节角度值得到标定后的机器人末端理论位置.标定前后的理论位置与观测位置的误差数据如图8所示，标定前最大位置误差为5.23 mm，平均误差为2.03 mm；标定后最大位置误差为1.03 mm，减小约80.31%，平均误差为0.35 mm，减小约8.76%.由实验数据可知，标定后机器人末端绝对位置精度有大幅度提升.

标定前后位置误差概率分布如图9所示.标定前位置误差方差为1.14，标定后为0.03，减小了约97.37%，可见标定后位置误差变小的同时，分布更加集中，提高了机器人末端位置的确定性.

再次随机生成具有20个测试点的验证集数据，用于验证标定得到的误差参数是否适用于一般情况.测试结果如图10～图11所示，当采用标准DH模型时，理论位置与测量位置之间的最大误差为3.37 mm，均值为1.78 mm，分布方差为0.7.采用引入关节几何误差的运动学模型，各几何误差采用表2中的标定值，测得机器人末端位置误差最大为0.66 mm，均值为0.36 mm，分布方差为0.02.与标准DH模型相比，最大位置误差减小约80.42%，均值减小约79.78%，位置误差分布方差减小约97.14%.可见，基于本文方法构建的运动学模型及标定方法得到的参数适用于一般情况，可大幅提升机器人的绝对位置精度和位置确定性.

4 结 论

1） 本文提出了一种引入关节几何误差构建的运动学模型，能够更准确地描述工业机器人的运动学特性.该模型的正逆运动学解实时性强，满足工程应用需求.

2） 基于误差参数的微分运动学雅可比矩阵可将关节误差、基坐标系误差和工具坐标系误差纳入标定范围.采用基于最小二乘法的绝对位置精度标定算法，可以成功对工业机器人的45个参数进行标定.

3） 根据本文设计的运动学模型和标定方法，在新松SR10C机器人的标定结果中，最大位置误差减小了约80%，平均误差减小了约80%，误差分布方差减小了约97%.结果表明，该算法能够显著提升工业机器人的绝对位置精度和确定性.

参考文献

[1] Roth Z， Mooring B， Ravani B. An overview of robot calibration［J］. IEEE Journal on Robotics and Automation， 1987， 3（5）： 377-385.

[2] Dumas C， Caro S， Chérif M， et al. A methodology for joint stiffness identification of serial robots［C］//2010 IEEE/RSJ International Conference on Intelligent Robots and Systems. Taipei， 2010： 464-469.

[3] Yang K， Yang W Y， Cheng G D， et al. A new methodology for joint stiffness identification of heavy duty industrial robots with the counterbalancing system［J］. Robotics and Computer-Integrated Manufacturing， 2018， 53： 58-71.

[4] Ginani L S， Motta J M S T. Theoretical and practical aspects of robot calibration with experimental verification［J］. Journal of the Brazilian Society of Mechanical Sciences and Engineering， 2011， 33（1）： 15-21.

[5] Ye S H， Wang Y， Ren Y J， et al. Robot calibration using iteration and differential kinematics［J］. Journal of Physics： Conference Series， 2006， 48： 1-6.

[6] Jang J H， Kim S H， Kwak Y K. Calibration of geometric and non-geometric errors of an industrial robot［J］. Robotica， 2001， 19（3）： 311-321.

[7] Gong C H， Yuan J X， Ni J. Nongeometric error identification and compensation for robotic system by inverse calibration［J］. International Journal of Machine Tools and Manufacture， 2000， 40（14）： 2119-2137.

[8] Denavit J， Hartenberg R S. A kinematic notation for lower-pair mechanisms based on matrices［J］. Journal of Applied Mechanics， 1955， 22（2）： 215-221.

[9] Hayati S， Mirmirani M. Improving the absolute positioning accuracy of robot manipulators［J］. Journal of Robotic Systems， 1985， 2（4）： 397-413.

[10] Brockett R W. Robotic manipulators and the product of exponentials formula［C］// Mathematical Theory of Networks and Systems. Berlin： Springer Berlin Heidelberg， 1984： 120-129.

[11] 高雷阜. 最优化理论与方法［M］. 沈阳： 东北大学出版社， 2005： 102-130.

[12] Gao Lei-fu. Theory and methods of optimization［M］. Shenyang： Northeastern University Press， 2005： 102-130.

[13] Omodei A， Legnani G， Adamini R. Three methodologies for the calibration of industrial manipulators： experimental results on a SCARA robot［J］. Journal of Robotic Systems， 2000， 17（6）： 291-307.

[14] Nguyen H N， Zhou J， Kang H J. A calibration method for enhancing robot accuracy through integration of an extended Kalman filter algorithm and an artificial neural network［J］. Neurocomputing， 2015， 151： 996-1005.

[15] 约翰 J. 克雷格. 机器人学导论：第4版［M］. 北京：机械工业出版社，2018：71-89.

[16] John J. Craig. Introduction to robotics［M］.4th ed. Beijing： China Machine Press， 2018： 71-89.

[17] Nubiola A， Bonev I A. Absolute calibration of an ABB IRB 1600 robot using a laser tracker［J］. Robotics and Computer-Integrated Manufacturing， 2013， 29（1）： 236-245.

[18] 毛晨涛. 工业机器人的运动学校准方法研究［D］. 杭州： 浙江大学， 2020.

[19] Mao Chen-tao. Research on kinematics calibration method of industrial robot［D］. Hangzhou： Zhejiang University， 2020.

[20] Sun Y， Hollerbach J M. Observability index selection for robot calibration［C］//2008 IEEE International Conference on Robotics and Automation. Pasadena， 2008： 831-836.

[21] Wang W， Liu F， Yun C. Calibration method of robot base frame using unit quaternion form［J］. Precision Engineering， 2015， 41： 47-54.

基金资助

国家自然科学基金资助项目(U20A20197)