工业机器人以其低成本、高效率、高柔性的优势广泛应用于机械制造、航空航天、海洋探测等多个领域。随着现代飞机对性能和寿命，以及航空制造过程对飞机制造的高质量、高效率、低成本的要求越来越高，机器人自动钻铆技术以其独有的优势成为飞机数字化、柔性化装配的必然选择^[1]。然而，自动钻铆系统普遍采用的6自由度串联机器人存在机构弱刚性问题，很大程度影响着机器人钻铆加工精度和加工质量^[2]。为此，通过分析机器人的静刚度和优化机器人系统作业时的位姿来增强机器人刚度性能，对于改善系统加工质量、保证系统作业稳定性具有重要意义。

在机器人刚度建模方面，Chen等提出守恒合同变换法(conservation contract transformation，CCT)，建立了6自由度旋转关节机器人6×6的Descartes刚度矩阵^[3]。Abele等完成了机器人铣削系统工业机器人的刚度建模及参数辨识，分析了机器人末端滑移变形规律，并对其加工性能作出评价^[4]。曲巍葳等分析了KUKA KR360的奇异性和灵巧性，通过实验法获得机器人关节刚度，并利用遗传算法对机器人刚度性能进行了优化^[5-6]。布音研究了机器人制孔系统中机器人运动学性能、刚度建模以及刚度辨识，获取机器人关节刚度参数，分析了机器人作业典型位姿下刚度特性对制孔质量的影响规律^[7]。

机器人刚度是衡量机器人负载能力的重要指标，也是影响机器人钻铆加工质量的关键参数。本文以KUKA KR500L340型钻铆机器人为研究对象，建立了静刚度模型；然后利用机器人系统、API激光跟踪仪和Kistler测力仪搭建实验平台，完成了机器人关节刚度辨识实验，获取了机器人关节刚度值；在此基础上，结合刚度性能评价指标，采用Monte Carlo法仿真得到机器人工作空间内的刚度性能云图，为在机器人作业时选择较优刚度性能区域给出参考；进一步结合机器人钻铆系统作业时的具体约束条件，采用粒子群算法优化了机器人的位姿，通过对优化前后机器人末端承受相同载荷时的变形量进行对比，证明了通过优化机器人位姿可实现增强其刚度性能的目的，对改善机器人系统作业质量和稳定性具有重要意义。

1 机器人静刚度建模

KUKA KR500L340型钻铆机器人为6自由度旋转关节机器人。本文首先采用改进Denavit-Hartenberg(DH)参数法^[8-9]在机器人各关节处建立坐标系，进而得到其运动学模型并列出DH参数。该机器人具体结构如图 1所示，DH参数见表 1，运动学模型如图 2所示。

根据机器人运动学模型和DH参数推导出机器人运动学正解，采用微分变换法构造机器人运动学Jacobi矩阵(机器人关节速度与末端操作速度的映射关系式^[10])。由于重载机器人连杆是由铸铁和高强度钢等材料构成，关节处由电机、减速器等传动机构组成，在承受载荷时，机器人关节处的变形远大于连杆变形，因此机器人静刚度建模时假定连杆为刚性^[11]，且机器人末端变形全部由关节变形引起。根据Jacobi矩阵和Hooker定律有关公式，建立机器人静刚度模型为

$ \boldsymbol{F}=\boldsymbol{J}^{-\mathrm{T}} \boldsymbol{K}_{\theta} \boldsymbol{J}^{-1} \boldsymbol{X} . $

(1)

式中：F为机器人末端所受外力矢量，J为机器人当前位姿下的运动学Jacobi矩阵，K_θ为机器人关节刚度矩阵，X为机器人末端在外力矢量F下的变形。

机器人运动学Jacobi矩阵随位姿变化而变化，机器人末端变形不仅与外力有关，还与其位姿有关。

2 机器人关节刚度辨识

获取机器人关节刚度值的方法有2种：1) 通过关节处的结构、材料和尺寸等参数进行计算^[12]；2) 设计关节刚度辨识实验^[13]，通过处理实验数据得到。由于第1种方法过于理想化，关节处具体结构等参数很难完全掌握，因此本文选择通过关节刚度辨识实验的方法获取机器人关节刚度值。

2.1 关节刚度辨识原理

将静刚度模型(式(1))变形处理得到

$ \boldsymbol{X}=\boldsymbol{J} \boldsymbol{K}_{\theta}^{-1} \boldsymbol{J}^{\mathrm{T}} \boldsymbol{F} . $

(2)

式(2)详细展开可以得到

$ \left[\begin{array}{c} \boldsymbol{U}_{1} \\ \boldsymbol{U}_{2} \\ \boldsymbol{U}_{3} \\ \boldsymbol{U}_{\mathrm{R} 1} \\ \boldsymbol{U}_{\mathrm{R} 2} \\ \boldsymbol{U}_{\mathrm{R} 3} \end{array}\right]=\left[\begin{array}{ccc} \boldsymbol{J}_{11} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 1} \bf{F}_{i} & \cdots & \boldsymbol{J}_{16} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 6} \bf{F}_{i} \\ \vdots & & \vdots \\ \boldsymbol{J}_{61} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 1} \boldsymbol{F}_{i} & \cdots & \boldsymbol{J}_{66} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 6} \boldsymbol{F}_{i} \end{array}\right]\left[\begin{array}{c} \bf{K}_{\theta_{1}}^{-1} \\ \vdots \\ K_{\theta_{6}}^{-1} \end{array}\right] . $

(3)

其中：U₁、U₂、U₃、U_R1、U_R2和U_R3为变形向量X的6维分量。

进一步处理得到

$ \boldsymbol{X}=\boldsymbol{A} \boldsymbol{c} . $

(4)

式中：$ \boldsymbol{A}=\left[\begin{array}{ccc} \boldsymbol{J}_{11} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 1} \boldsymbol{F}_{i} & \cdots & \boldsymbol{J}_{16} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 6} \boldsymbol{F}_{i} \\ \bf{\vdots} & & \vdots \\ \boldsymbol{J}_{61} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 1} \boldsymbol{F}_{i} & \cdots & \boldsymbol{J}_{66} \sum\limits_{i=1}^{6} \boldsymbol{J}_{i 6} \boldsymbol{F}_{i} \end{array}\right]$，为Jacobi矩阵和力矢量所组成的系数矩阵；$ \boldsymbol{c}=\left[\begin{array}{c} \boldsymbol{K}_{\theta_{1}}^{-1} \\ \vdots \\ \boldsymbol{K}_{\theta_{6}}^{-1} \end{array}\right]$，为机器人关节刚度矩阵的变形形式。

通过关节刚度辨识实验测量机器人在不同位姿下的力与变形数据，即可求解得到c，进而得到关节刚度矩阵K_θ。

2.2 关节刚度辨识实验

机器人关节刚度辨识实验分3步：

步骤1   标定机器人基坐标系与激光跟踪仪坐标系的转换关系；

步骤2   标定测力仪坐标系与激光跟踪仪坐标系的转换关系；

步骤3   测量机器人末端的力与变形。

图 3所示为机器人关节刚度辨识实验现场。保持标定完成后的机器人系统、激光跟踪仪和测力仪等设备的位置不变，将激光跟踪仪靶镜安装在末端执行器靠近机器人末端法兰处。实验过程中，首先调节机器人各关节角来调整位姿使末端执行器靠近工件表面，记录此时的给关节角与末端位姿参数；然后控制压脚机构使其贴于工件表面并施加压紧力，记录测力仪与激光跟踪仪的数据；上述步骤重复6次，最终测得6组位姿，每组位姿可以得到16个测量点的力与变形数据，由于扭转变形量U_R1、U_R2、U_R3无法测得，因此取U₁、U₂、U₃，总计可以得到由288个方程构成的超定方程组。

2.3 辨识结果

首先求解超定方程组。采用最小二乘法可以解得c，进一步将其变换成矩阵形式求逆，得到机器人关节刚度矩阵K_θ，最终得到的KUKA KR500L340机器人的关节刚度矩阵为

$ \begin{aligned} &\boldsymbol{K}_{\theta}=\operatorname{diag}\left(7.65 \times 10^{6}, 5.76 \times 10^{6}, 2.43 \times 10^{7},\right. \\ &\left.8.49 \times 10^{6}, 3.43 \times 10^{6}, 3.88 \times 10^{6}\right) \mathrm{N} \cdot \mathrm{m} / \mathrm{rad} . \end{aligned} $

(5)

由式(5)可知，该机器人关节刚度在10⁶ N·m/rad量级，关节刚度较大，能够承受较大关节负载。为验证所辨识出的关节刚度矩阵准确性，增加测得一组位姿下机器人末端的力与变形数据。通过已辨识出关节刚度矩阵和Jacobi矩阵以及测力仪的力矢量数据，可以计算得到机器人在该位姿下的理论预测变形量，并与激光跟踪仪测得的实际变形量进行对比，对比结果如图 4所示。预测值与实际值的平均误差仅为3.2%，辨识得到的关节刚度矩阵能够较好地预测机器人末端受力变形情况，证明了通过关节刚度辨识实验所得到的关节刚度模型的正确性。

3 机器人工作空间刚度性能分析

结合所辨识出的关节刚度矩阵和机器人处于不同位姿时的Jacobi矩阵，根据式(6)可以得到机器人在不同位姿下的末端操作刚度矩阵K_x，

$ \boldsymbol{K}_{x}=\boldsymbol{J}^{-\mathrm{T}} \boldsymbol{K}_{\theta} \boldsymbol{J}^{-1} . $

(6)

3.1 刚度性能评价指标

末端操作刚度矩阵K_x为六阶方阵，并不能直接判断出机器人在当前位姿下刚度性能优劣，需引入刚度性能评价指标^[14-15]。根据操作刚度矩阵的性质和元素量纲，结合式(1)可以得到

$ \left[\begin{array}{l} \boldsymbol{f} \\ \boldsymbol{m} \end{array}\right]=\left[\begin{array}{ll} \boldsymbol{k}_{f d} & \boldsymbol{k}_{f \hat{o}} \\ \boldsymbol{k}_{m d} & \boldsymbol{k}_{m \delta} \end{array}\right]\left[\begin{array}{l} \boldsymbol{d} \\ \boldsymbol{\delta} \end{array}\right] . $

(7)

式中：f为机器人末端力矢量，N；m为机器人末端力矩矢量，N·m；k_fd为力-线位移刚度矩阵，N/m；k_fδ为力-角位移刚度矩阵，N/rad；k_md为力矩-线位移刚度矩阵，N；k_mδ为力矩-角位移刚度矩阵，N·m/rad；d为机器人末端移动变形矢量，m；δ为机器人末端转动变形矢量，rad。

为重点研究力与平移变形关系，令δ=0，由式(7)可以得到

$ \boldsymbol{f}=\boldsymbol{k}_{f d} \cdot \boldsymbol{d} . . $

(8)

假定机器人末端承受单位外力作用，即f=1，则有

$\left|\boldsymbol{f}\right|^{2}=\boldsymbol{f}^{\mathrm{T}} \boldsymbol{f}=\boldsymbol{d}^{\mathrm{T}} \boldsymbol{k}_{f d}^{\mathrm{T}} \boldsymbol{k}_{f d} \boldsymbol{d}=1 \text {. } $

(9)

进一步分析可知，式(9)描述的为一个椭球体，可以称之为变形椭球^[16]，表示机器人末端在承受单位外力作用时的变形情况。椭球的长半轴对应着力-线位移刚度矩阵k_fd的最小特征值λ_min(k_fd)，表示机器人在当前位姿下末端承受单位外力时的最大变形量，即当前位姿的最差刚度性能，可以将其作为刚度性能评价指标，λ_min(k_fd)越小，表示机器人在当前位姿下末端承受单位外力时的最大变形量越小，机器人的刚度性能越好。

同理，将式(7)进行变形处理得到

$ \left[\begin{array}{l} \boldsymbol{d} \\ \boldsymbol{\delta} \end{array}\right]=\left[\begin{array}{cc} \boldsymbol{c}_{f d} & \boldsymbol{c}_{m d} \\ \boldsymbol{c}_{f \hat{o}} & \boldsymbol{c}_{m \delta} \end{array}\right]\left[\begin{array}{l} \boldsymbol{f} \\ \boldsymbol{m} \end{array}\right] . $

(10)

其中：c_fd为力-线位移柔度矩阵，N/m；c_md为力矩-线位移柔度矩阵，N；c_fδ为力-角位移柔度矩阵，N/rad；c_mδ为力矩-角位移柔度矩阵，N·m/rad。

可见，力-线位移柔度矩阵c_fd也对应着一个椭球体，可以称之为力椭球，表示使机器人末端产生单位平移变形所需的外力。椭球的短半轴对应着矩阵c_fd的最大特征值λ_max(c_fd)，表示引起机器人末端单位平移变形所需的最小外力，可以将其作为刚度性能评价指标。λ_max(c_fd)越大，说明使机器人末端产生单位平移变形所需的外力越大，机器人的刚度性能越好。

3.2 机器人工作空间刚度性能分析

本文利用Monte Carlo法和机器人DH参数仿真分析机器人工作空间^[17]，根据机器人工作空间内每一点位姿和操作刚度矩阵，得到机器人工作空间内的刚度性能云图，图 5和6所示为分别采用λ_min(k_fd)和λ_max(c_fd)为评价指标得到的该机器人的工作空间刚度性能云图和XOZ平面投影。

当以λ_min(k_fd)为评价指标时，深蓝色区域代表刚度性能较好区域，橙黄色区域代表刚度性能较差区域，而以λ_max(c_fd)为评价指标时，橙黄色区域代表刚度性能较好区域，深蓝色代表刚度性能较差区域。由图 5和6的分析结果可知，两种评价指标下刚度性能较好和较差区域完全一致。

4 基于刚度性能的位姿优化 4.1 优化问题数学模型

机器人钻铆系统作业要求刀具刀尖点与目标加工点重合，刀具进给方向与目标点法线方向一致，此时约束了机器人5个自由度，还存在绕目标点法线方向旋转自由度没有约束。此时，对于6自由度机器人，满足加工条件的位姿有无穷多个，为此以刚度性能最大化为目标对机器人作业位姿进行优化。

优化数学模型包括目标函数、优化变量和约束条件。选择力-线位移刚度矩阵的最小特征值λ_min(k_fd)作为刚度性能评价指标和优化目标函数；优化变量为机器人关节角θ，对应着机器人位姿；约束条件为机器人各个关节角限位和加工时位姿约束，即刀尖点与目标点重合、刀具进给方向与目标点法线方向一致。最终得到优化问题的数学模型为

$ \left\{ \begin{array}{l} \min \;\;\;\;\;\;\;\;\;\;\;\;\;\;f(\theta )\\ {\rm{s}}{\rm{. t}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;H(\theta ) = 0,\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\theta _{i\min }} \leqslant {\theta _i} \leqslant {\theta _{{i \mathop{\rm max}\nolimits} }}\quad (i = 1,2, \cdots ,6). \end{array} \right. $

(11)

4.2 优化算法

粒子群优化(particle swarm optimization，PSO)算法^[18]具有较强全局搜索能力，较快收敛速度，且参数设置相对较少，易于编程求解。

4.3 优化结果

选取一组机器人常用作业位姿，根据优化模型和算法，利用MATLAB求解得到优化后机器人位姿。表 2为位姿优化前后机器人各关节角和刚度性能指标对比，机器人位姿状态对比见图 7。

4.4 优化结果验证

假设在优化前后两组位姿下机器人末端执行器上刀尖点处，沿进给反方向施加400~1 000 N的外力，结合两组位姿下Jacobi矩阵和关节刚度矩阵，计算得到优化前后两组位姿下机器人末端理论变形量，结果对比如图 8所示。可见，优化后的变形量明显降低。

5 结论

结合机器人Jacobi矩阵和Hooker定律建立了某型钻铆机器人静刚度模型；设计并完成了机器人关节刚度辨识实验，获取了机器人关节刚度值，验证了所辨识出的关节刚度值的准确性；以力-线位移柔度矩阵的最小特征值和力-线位移柔度矩阵的最大特征值两种刚度性能评价指标，结合Monte Carlo法得到机器人工作空间内刚度性能云图。以机器人刚度性能为目标，采用粒子群算法优化了机器人位姿。优化前后机器人末端受力变形情况分析结果表明，机器人位姿优化可提高机器人作业刚度，对保证机器人系统作业质量和稳定性都具有重要意义。