考虑制造误差的工业机器人精度可靠性分析＊

李国发 韩良晟 何佳龙 王继利

（①吉林大学数控装备可靠性教育部重点实验室，吉林 长春 130015；②吉林大学机械与航空航天工程学院，吉林 长春 130015）

工业机器人的定位精度是衡量机器人性能的重要指标，包括（绝对）定位精度和重复定位精度。工业机器人的定位精度可靠性是指工业机器人末端执行器定位点落在精度域内的概率，是从概率的角度对定位精度进行量化[1]。工业机器人的运动学模型是用于描述工业机器人各关节运动与末端执行器之间的关系的模型，然而由于工业机器人结构复杂以及加工和装配误差等不确定因素的存在，使得机器人理想运动学模型和实际运动学模型之间存在误差，建立机器人运动学误差模型是机器人定位精度可靠性研究的基础。Mei B提出了一种针对五轴并联机械手弹性几何误差的建模方法[2]，揭示了结构误差、重力作用下的弹性变形和柔量参数误差对位姿偏差的综合影响；Li Y建立了考虑平行度误差的Delta机器人的运动学模型并进行了误差分析[3]；Fu Z基于Lie理论提出一种机器人统一误差模型[4]，并用于标定机器人的运动学参数和关节柔度；Yan Y基于矢量矩阵分析构建了六自由度并联机械手的误差模型[5]，并依据所建立的误差模型采用遗传算法对六自由度并联机器人的精度进行了优化。

建立机械产品的误差模型通常采用公差分析的思想。公差分析的数学模型主要包括实体漂移模型[6]、矩阵模型[7]、矢量环模型[8]、小位移旋量模型[9]、多面体模型[10]和公差变动矢量图模型[11]等。小位移旋量模型通过将公差变动划分为旋量与平移量两部分，利用机器人运动学方法将公差旋量转换成4×4的齐次变换矩阵，通过各个齐次变换矩阵建立三维空间的误差传递模型，并与机器人运动学模型结合建立机器人运动学误差模型。

目前，利用代理模型进行精度可靠性分析已成为热点。代理模型是指通过采样样本，利用插值或拟合等方法构建功能函数的代替模型，然后利用蒙特卡洛仿真对所构建的替代模型进行仿真分析估计失效概率，从而降低计算成本，提高计算效率，同时对于复杂的非线性多维度函数具有较高的计算精度。李国发[12]提出了一种面向多种代理模型的自适应加点策略用于结构可靠性分析。袁修开[13]提出了一种Kriging法与改进一次二阶矩融合的方法用于结构可靠性分析。常用的代理模型包括多项式响应曲面法[14]、Kriging[15]、梯度增强克里金法[16]、支持向量机[17]和人工神经网络[18]等。

本文基于小位移旋量理论（SDT）与MD-H模型相结合建立了考虑包括形状误差和尺寸误差的六自由度工业机器人关节误差的运动学模型，并基于蒙特卡洛法和Kriging代理模型对机器人进行了定位精度可靠性分析，通过灵敏度分析寻找到了影响工业机器人定位精度可靠性的主要误差来源。

1 机器人运动学模型

1.1 理想状态机器人正运动学模型

D-H模型[19]主要被用于在机器人连杆上建立坐标系，对于相邻连杆间的坐标变换的齐次变换矩阵为：

式中：将 sin简写为s，cos简写为c。其中α、d、a，θ表示机器人的DH参数。

则机器人的正运动学方程可以表示为：

1.2 考虑关节误差的机器人正运动学模型

SDT是利用机器人运动学中的齐次坐标变换方法描述三维实体各处误差的累积，同时与机器人的齐次变换矩阵相结合，能够用于建立考虑关节间隙的机器人运动学模型。小位移旋量D由两组矢量（即3个平动矢量和3个旋转矢量）精确描述。其对应的表达式和相应的齐次变换矩阵如下：

零部件之间的配合表面主要包含平面配合、圆环面配合和柱面配合，利用SDT建立的误差模型应符合设计阶段的公差要求，即相应的旋转和位移矢量所建立的最终的误差模型应满足设计公差要求。本文将机器人每个关节间的配合简化为圆环面和柱面配合，则对于圆环面配合的小位移旋量u、v、γ为零，对于柱面配合的小位移旋量 γ为零。

D-H模型描述了机器人相邻两个关节3个方向的平移及绕X轴和Z轴的旋转，不包括绕Y轴的转动，然而在实际运行中，机器人相邻关节间也存在绕Y轴的转动，此时传统的D-H模型就无法准确描述机器人的运动学模型。因此，在D-H模型的基础上，增加一个绕Y轴的旋转矩阵，构建了MD-H模型[20]，其对应的齐次变换矩阵如下式所示：

本文对每个圆环面配合和柱面配合进行标号。将机器人末端中心作为工具坐标点。机器人各关节间的误差传递路径包括串联路径和并联路径，由于机器人关节是由圆环面和柱面耦合配合，圆环面配合会产生绕X轴Y轴的转动误差以及绕Z轴的移动误差，而柱面配合会产生绕X轴转动的误差和沿着X轴和Y轴的移动误差，因此，本文对并联传递的部分旋转变量进行求交运算而对平移变量进行了求和运算。

2 机器人定位精度可靠性分析

2.1 六自由度工业机器人定位精度可靠性功能函数

针对建立包含关节误差的机器人正运动学模型，通过仿真对机器人的定位精度可靠性进行分析。先构建机器人定位精度可靠性分析的功能函数，如式（6）所示：

式中： [s]是 机器人末端点的允许误差，s(x)是机器人末端点的实际位置，s0是机器人末端点的实际位置，x=[Δx1,···,Δx6,Δy1,···,Δy6,Δz1,···,Δz6,Δα1,···,Δα6,Δβ1,···,Δβ6]是包含关节误差的向量，其中包括30个参数，各关节误差变量服从高斯分布，则机器人的定位的失效概率可以表示为：

2.2 试验设计

试验设计的核心内容是选择试验样本点的策略。样本点的选取应遵循样本点尽可能地反映整个样本空间的信息。常用的试验设计方法包括均匀设计、正交设计和拉丁超立方抽样设计等。拉丁超立方抽样方法因所抽取的样本点均匀的分布在整个样本空间内，被广泛使用，本文采用拉丁超立方抽样设计进行试验设计，其具体步骤为：首先确定抽样的变量个数n以及准备抽取的样本点个数m，然后将每个随机变量的区间m等分并在每个随机变量的每个区间内随机的取值，最终将取得的值随机组合形成最终的样本点用于后续Kriging代理模型的构建。

2.3 Kriging代理模型的构建

Kriging模型是依据协方差函数对随机过程或随机场进行空间建模和预测的回归算法。对于六自由度工业机器人的功能函数所涉及的计算量较大，而对于其定位精度的失效判定只需考虑其功能函数是否是小于0即可，Kriging代理模型可以根据给定的输入变量和响应拟合出代理模型，用于计算机器人定位精度可靠性的功能函数，进而进行定位精度可靠性分析。Kriging代理模型构建的步骤参如图1所示。

图1 Kriging代理模型构建流程

Kriging模型包括回归部分和随机部分，随机部分通常包括高斯过程和指数过程，本文中采用高斯过程。设输入变量为X=[x1,x2,···,xn]T，对应的真实响应值为G=[G(x1),G(x2),···,G(xk)]，Kriging代理模型的结构形式如下式所示：

式中：F(β,x)为 线性回归模型；fT(x)为多项式基函数；x为输入变量；G＾(x)为 对预测值的响应值；β是回归系数向量，其中任意两个样本点之间的协方差如下式所示：

对于给定的样本点X与对应的真实响应G，对于给定的相关性参数 θ和回归系数 β的估计值如下式所示：σ2

为保证在某点预测值无偏并且预测均方误差最小，在该点处的预测的均值和方差满足下式：

式中：r(x)=[R(θ,xi,x),R(θ,xi,x),···,R(θ,xi,x)]；u(x)=FTR-1r-f；Kriging模型中方差表明了模型在预测点处的误差。

3 机器人定位精度可靠性灵敏度分析

根据建立的机器人定位精度功能函数的Kriging代理模型，结合候选样本点，利用Monte Carlo抽样法获取随机参数均值与方差的灵敏度。

设随机参数向量为x=[x1,x2,···,xn]T，其联合概率密度为fx(x)，则失效概率为：

式中：F表示失效域，在本文中为功能函数小于零的部分，即F={x|G(x)≤0}。随机参数的可靠性灵敏度实质上是该系统失效概率对随机参数的均值以及方差求偏导运算，如式（16）所示：

式中：I[g(x)]如式（17）所示。

在实际工程中，常认为误差的各随机变量为独立的，所以其联合概率密度为各随机变量概率密度的乘积，如下式所示：

利用拉丁超立方抽样方法进行随机抽样，并将得到的样本点代入Kriging代理模型中，即得到相应的预测值，再对所得到的预测值进行分析和分类，根据大数定律可知，灵敏度计算可以表示为：

对于相互独立的正态随机变量， θxi包括均值μxi和标准差 σxi，可以得到下式：

代入求得系统对于参数均值和标准差失效概率灵敏度，如下式所示：

式中：NF表示利用联合概率密度函数fx(x)抽取的且满足功能函数小于0即发生失效的样本数量。对Sμxi和Sσxi取模得到：

随机参数灵敏度因子：

4 案例分析

4.1 机器人理想运动学模型

本文以RB08型六自由度工业机器人为例，建立了其理想状态下的运动学模型，机器人的具体参数见表1，建立的模型结果见图2。

表1 RB08机器人DH参数表

图2 理想机器人运动学模型

表1中：a1=170mm,a2=560mm,a3=155mm ，d4=630mm,d6=110mm。

4.2 包含关节误差的机器人运动学模型

将机器人关节间的配合表面统一简化为柱面配合和圆环面配合，并进行编号，详见表2。

表2 各个配合表面编号

建立机器人各关节的误差传递路径，如图3所示。根据1.2中所述，以 1→2误差传递为例。由1→2共包含2条误差传递路径，两条为并联路径。对于并联路径的旋量求交运算，平移量求和运算对于串联部分旋量和平移量都进行求和运算：

图3 机器人误差传递路径

式中：T表示误差；D表 示平移误差；R表示旋量误差。

由包含关节误差的机器人齐次变换矩阵式（4）和式（5）可以得到机器人的包含关节误差的正运动学模型，如下式所示：

以机器人各个关节间隙服从高斯分布进行了3 000次的仿真计算，得到的六自由度工业机器人相对理想位置的误差分布情况，其在空间的分布情况如图4。

图4 机器人末端点的位置变动

设置仿真时间为3 s，各关节匀速运动，各个关节变量见表3，一次随机仿真后机器人末端点在各方向上的理想运动与包含误差的运动对比结果见图5。

图5 时变状态下机器人末端点在各方向上的位置变动

表3 机器人关节初始角度和最终角度 （°）

4.3 构建功能函数代理模型

采用拉丁超立方抽样设计试验，获得了1 000组输入变量及其真实的响应，并用于对Kriging模型进行训练，将训练得到的代理模型用于预测六自由度工业机器人定位精度失效概率，对给定的2 000组输入变量的真实响应和预测响应结果见图6。在2 000次试验中基于功能函数及代理模型对于发生定位失效的真实判断结果和预测结果对比见图7。从图6可以看出，所建立的包含关节误差的机器人运动学代理模型很好的代替了机器人的运动学误差模型。从图7中可以看出，在2 000次仿真试验中，代理模型可以准确地识别出在试验过程中机器人是否发生定位精度失效。

图6 真实响应和预测响应对比

图7 2000次试验真实失效和预测结果对比

4.4 机器人定位精度可靠性灵敏度分析

对30个关节变量灵敏度分析的结果如图8和图9所示，其中图9中对各灵敏度结果进行了取正。结果表明六自由度工业机器人定位失效率对前18组变量的随机参数标准差比较敏感，经过计算得到前18组变量灵敏度因子占比为85.78%，而后12组变量灵敏度因子占比仅为14.22%。即在公差允许范围内，各关节的移动变量对机器人定位精度可靠性的影响比较明显。

图8 无量纲灵敏度条形图

图9 灵敏度因子条形图

然而针对前18组关节的移动变量，其灵敏度因子也存在较大的差异，这是由于机器人关节误差样本的随机产生以及机器人定位姿态影响的结果。如图10所示，利用建立的机器人运动学误差模型，随机指定一组关节运动，以第三和第四组误差变量为例，在2 s的运动内，第三个误差变量的灵敏度因子产生了显著的变化。即在同样一组误差变量的条件下，机器人的姿态同样会对机器人定位精度可靠性灵敏度分析结果产生影响。由于工业机器人通常从事重复性的工作，即机器人的动作通常是单一的，因此，可以使用本文方法对机器人在系列姿态下进行定位精度可靠性灵敏度的分析。

图10 不同姿态下第三、第四个误差灵敏因子结果

5 结语

针对工业机器人各关节加工和装配误差影响六自由度工业机器人末端精度可靠性的问题，本文考虑了机器人各关节处尺寸误差和形状误差。构建了机器人定位精度可靠性功能函数的Kriging代理模型。并开展了机器人定位精度可靠性分析。

（1）利用MD-H法和小位移旋量理论，建立了工业机器人包含关节尺寸误差和形状误差的运动学模型，提出了各关节误差的传递路径以及误差累计的方法，相对于传统方法更加接近实际情况。

（2）建立了基于蒙特卡洛和Kriging代理模型的工业机器人定位精度可靠性代理模型，通过蒙特卡洛仿真试验结果表明，建立的精度可靠性功能函数代理模型具有很高的预测精度。

（3）开展了基于代理模型的工业机器人定位精度灵敏度分析，确定了影响机器人定位精度可靠性的主要误差来源：各关节处的移动误差，验证了建立的工业机器人定位精度可靠性分析方法适用性强，具有较强的通用性。