近年来，工业机器人以其低成本、高灵活性的特点被广泛应用于制孔、铣削、磨削等加工领域. 在飞机装配中，飞机结构件的连接需要大量孔. 人工钻孔耗时耗力且难以保证孔的质量，大型专用机床虽然具有较好的加工质量，但需要较大的安装空间和资金投入. 因此，越来越多的机器人制孔系统被开发用于航空制造领域，以适应低成本、高质量、高灵活性的自动化加工需求^[1-3].

与传统机床相比，工业机器人存在刚性差和易变形的问题^[4-5]. 在机器人制孔过程中，切削产生的动态切削力作为激振力易引发机器人颤振^[6-7]，须在机械臂末端加装压脚机构压紧工件表面以提升机器人系统动态刚度，利用压脚与工件接触面之间的摩擦力平衡动态切削力，达到抑制振动的目的. 此外，压脚还能补偿重力对末端执行器角度造成的影响以及消除叠层材料层与层之间的间隙^[8]. 但是，压脚产生的压紧力会使机器人关节发生变形，导致机器人末端位置产生不可忽略的偏差，从而降低最终制孔的位置精度.

在机器人运行过程中，由于存在负载变形、几何误差、关节磨损、热误差等原因，其绝对定位精度往往只有2~3 mm. 对此，实验室前期已对机器人自动化制孔中绝对定位误差进行分析研究，提出多种补偿方法用于提高机器人制孔定位精度，但均未考虑压脚压紧后产生的二次定位误差^[9-11]. 机器人在压脚压紧力作用下的变形会导致末端执行器压脚产生轴向变形和滑移变形. 其中轴向变形导致制孔过程中锪窝深度产生偏差，而滑移变形造成制孔定位误差. 针对轴向变形，费少华等^[12]提出通过基于压脚位移反馈的补偿技术进行窝深控制. 对于滑移变形，Olsson等^[13]利用6自由度力/力矩传感器，通过外力控制回路和相应的干扰补偿方案来控制压脚水平方向的力始终小于压脚与工件之间的最大静摩擦力，从而抑制压脚滑移. 然而该方法只能用于特定的加工系统，且须构建复杂的控制系统，在工业中的应用性相对较低. 考虑到机器人受力变形与其刚度密切相关，因此，解决这个问题的可行方法是根据机器人刚度模型对其变形进行离线预测和补偿. 为了实现该方法，首先须对机器人系统进行刚度建模以及相应的关节刚度辨识.

在机器人刚度建模方面，Salisbury^[14]最早提出机器人刚度模型. 在此基础上，Chen等^[15]通过引入补充刚度矩阵，提出增强型刚度模型，促进机器人刚度建模的发展. Dumas等^[16]在对6自由度关节式工业机器人进行关节刚度辨识时发现，增强型刚度模型中的补充刚度矩阵在大多数情况下都可以忽略. Bu等^[17]在研究机器人刚度特性时，为了简化计算而忽略末端的旋转变形和所受扭矩，将传统刚度模型简化为只有位移和力之间的关系，并通过与完整刚度模型下的实验变形进行对比，证明该简化刚度模型的有效性. 在机器人关节刚度辨识方面，大量学者基于各自的刚度模型和方法进行关节刚度辨识，但多数采用在机器人末端悬挂重物的加载方式，实验的难度和危险性较大^[18-20]. 曲巍崴等^[21]在对装有末端执行器的工业机器人系统进行关节刚度辨识时，将测力仪测得的压脚作用力作为机器人所受外载，简化关节刚度辨识的实验装置，但由于未考虑压脚与工件之间的摩擦力，导致实验结果存在误差.

本研究以机器人刚度分析为基础，提出适用于绝大多数机器人制孔系统的变形预测和离线补偿方法，用于进一步提高机器人制孔系统的定位精度. 以装有末端执行器的KUKA KR600型机器人为研究对象，建立机器人运动学模型，通过研究机器人末端平移变形与压脚压紧力的关系，建立压脚压紧力作用下的机器人刚度模型；通过基于L-M算法的机器人关节刚度辨识实验获取其关节刚度，并规划实验验证关节刚度的准确性；运用该机器人刚度模型对机器人滑移变形进行预测并补偿. 通过补偿前、后机器人制孔位置误差对比，证明该补偿方法能有效提高机器人末端定位精度.

KUKA KR600型机器人为6自由度关节式工业机器人，且所有关节均为旋转副，按照改进的Denavit-Hartenberg方法（MOD_DH）对机器人进行运动学建模，得到机器人连杆坐标系，如图1所示. 图中，坐标系{O}为基坐标系，其余坐标系{O_i}对应关节i，i=1，2，···，6，a_i为关节i的连杆长度，d_i为关节i的连杆偏置. MOD_DH参数如表1所示. 表中，α_i−1为关节i−1的连杆扭角，θ_i为关节i的转角，以逆时针旋转为正.

由文献[22]可知，相邻坐标系{O_i-1}和{O_i}之间的齐次变换矩阵为

(1) $ \begin{split} {}_i^{i - 1}{{T}} = &{\rm{Rot }}\;({X_{i - 1}},{\alpha _{i - 1}}){\rm{Trans }}\;({X_{i - 1}},{a_{i - 1}}){\rm{Rot }}\;({Z_i},{\theta _i})\times\\ &{\rm{Trans }}\;({Z_i},{d_i}). \end{split} $

式中： ${\rm{Rot }}\;({X_{i - 1}},{\alpha _{i - 1}})$ 表示绕 ${X_{i - 1}}$ 轴旋转 ${\alpha _{i - 1}}$， ${\rm{Trans }}\;({X_{i - 1}},{a_{i - 1}})$表示沿 ${X_{i - 1}}$轴移动 ${a_{i - 1}}$， ${\rm{Rot }}\;({Z_i},{\theta _i})$表示绕 ${Z_i}$轴旋转 ${\theta _i}$， ${\rm{Trans }}\;({Z_i},{d_i})$表示沿 ${Z_i}$轴移动 ${d_i}$.

将各坐标系之间的齐次变换矩阵顺序相乘，即可得到机器人末端坐标系相对于基坐标系的变换矩阵，即机器人运动学方程为

(2) ${}_6^0{{T}} = {}_1^0{{T}}{}_2^1{{T}}{}_3^2{{T}}{}_4^3{{T}}{}_5^4{{T}}{}_6^5{{T}}.$

至此，完成机器人MOD_DH运动学模型建立，为后续机器人刚度分析奠定基础.

为了验证该运动学模型的正确性，将表1中6个关节转角参数代入机器人运动学方程，得到机器人末端在基坐标系下的位姿矩阵：

(3) ${}_6^0{{T}} = \left[ {\begin{array}{*{20}{c}} 0&0&1&{1\;815} \\ 0&1&0&0 \\ { - 1}&0&0&{2\;290} \\ 0&0&0&1 \end{array}} \right].$

同时，使用Matlab建立相应机器人运动学模型，并采用表1所示的MOD_DH关节转角参数，得到如图2所示结果. 可以看出，仿真结果与式（3）中的机器人末端位姿完全相同.

机器人刚度主要反映机器人在外力作用下抵抗变形的能力，机器人刚度模型则表示外力与机器人变形之间对应的数值关系. 在机器人制孔中，末端执行器压脚产生的压紧力使机器人发生变形，导致压脚与工件表面发生相对滑移而产生滑动摩擦力，从而进一步影响机器人的变形情况. 因此，对于有压脚约束的机器人刚度模型应综合考虑压脚压紧力和摩擦力对机器人变形产生的影响.

Salisbury^[14]最早建立机器人关节刚度和笛卡尔刚度的映射关系：

(4) ${{K}} = {{{J}}^{ - {\rm{T}}}}{{{K}}_\theta }{{{J}}^{ - 1}}.$

式中： ${{K}}$为机器人笛卡尔刚度矩阵； ${{J}}$为机器人雅可比矩阵； ${{{K}}_\theta }$为机器人关节刚度矩阵， ${{{K}}_\theta }= {\rm{ diag}}\;\left[ {{k_1},{k_2}, \cdots ,{k_6}} \right]$， ${k_1},{k_2}, \cdots ,{k_6}$分别为机器人6个关节的刚度.

由式（4）可以求出机器人在笛卡尔空间下的刚度，进而得到传统机器人刚度模型：

(5) $ \Delta {{X}} = {{{K}}^{ - 1}}{{F}} = {{JK}}_\theta ^{ - 1}{{{J}}^{\rm{T}}}{{F}} = {{CF}}. $

式中： $\Delta {{X}}$为机器人末端广义变形量， $\Delta {{X}} = \left[ \Delta x, \Delta y,\Delta z,\Delta \alpha ,\Delta \beta ,\Delta \gamma \right]^{\rm{T}}$， $\Delta x$、 $\Delta y$、 $\Delta z$、 $\Delta \alpha $、 $\Delta \beta $、 $\Delta \gamma $分别为沿X、Y、Z轴的平移量和绕X、Y、Z轴的旋转量； ${{F}}$为机器人末端广义力矢量， ${{F}} = \left[ {F_x},{F_y},{F_z}, {M_x},{M_y}, {M_z} \right]^{\rm{T}}$， ${F_x}$、 ${F_y}$、 ${F_z}$、 ${M_x}$、 ${M_y}$、 ${M_z}$分别为沿X、Y、Z轴的力和绕X、Y、Z轴的力矩； ${{C}}$为机器人柔度矩阵.

由文献[17]可知，机器人制孔时所受的扭矩及其末端的旋转变形可以忽略不计，因此本研究主要研究机器人平移变形与外力（不含力矩）之间的关系. 假设机器人所受力矩和旋转变形均为零，则式（5）可以写为

(6) $\left[ {\begin{array}{*{20}{c}} {\Delta {{{X}}_{\rm{t}}}} \\ {\bf{0}} \end{array}} \right] = {{C}}\left[ {\begin{array}{*{20}{c}} {{{{F}}_{\rm{t}}}} \\ {\bf{0}} \end{array}} \right].$

式中： $\Delta {{{X}}_{\rm{t}}}$为机器人末端平移变形， $\Delta {{{X}}_{\rm{t}}} = {\left[ {\Delta x,\Delta y,\Delta z} \right]^{\rm{T}}}$； ${{{F}}_{\rm{t}}}$为机器人所受作用力， ${{{F}}_{\rm{t}}} = {\left[ {{F_x},{F_y},{F_z}} \right]^{\rm{T}}}$.

将柔度矩阵 ${{C}}$写成分块形式：

(7) ${{C}} = \left[ {\begin{array}{*{20}{c}} {{{{C}}_{{\rm{tt}}}}}&{{{{C}}_{{\rm{tr}}}}}\\ {{{C}}_{{\rm{tr}}}^{\rm{T}}}&{{{{C}}_{{\rm{rr}}}}} \end{array}} \right].$

式中： ${{{{C}}_{{\rm{tt}}}}}$为平移柔度子矩阵， ${{{C}}_{{\rm{rr}}}}$为旋转柔度子矩阵， ${{{C}}_{{\rm{tr}}}}$为耦合柔度子矩阵.

结合式（6）、（7），可以得到机器人平移变形和所受外力之间的关系，即简化后的传统机器人刚度模型：

(8) $ \Delta {{{X}}_{\rm{t}}} = {{{C}}_{{\rm{tt}}}}{{{F}}_{\rm{t}}}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {{{{c}}_1}}\\ {{{{c}}_2}}\\ {{{{c}}_3}} \end{array}} \right]{{{F}}_{\rm{t}}} = \left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{{c_{13}}}\\ {{c_{21}}}&{{c_{22}}}&{{c_{23}}}\\ {{c_{31}}}&{{c_{32}}}&{{c_{33}}} \end{array}} \right]{{{F}}_{\rm{t}}}. $

式中： ${{{c}}_1}$、 ${{{c}}_2}$、 ${{{c}}_3}$分别为平移柔度子矩阵 ${{{C}}_{{\rm{tt}}}}$的第1、2、3行向量.

在机器人制孔前，末端执行器压脚垂直压紧工件表面，产生压紧力 ${{{F}}_0} = {\left[ {0{\kern 1pt} ,0{\kern 1pt} ,{F_z}{\kern 1pt} } \right]^{\rm{T}}}$，如图3所示. 压紧力由压脚气缸的推力和压脚重力沿Z向分力组成，在实验中可以由工件后方的测力仪测得.

由式（8）可以求得机器人在初始压紧力作用下的理论末端变形 $\Delta {{{X}}_0} = {\left[ {\Delta {x_0},\Delta {y_0},\Delta {z_0}} \right]^{\rm{T}}}$. 机器人末端变形将驱使压脚沿工件表面，即X-Y平面发生滑移，从而在压脚和工件表面产生滑动摩擦力 ${{{F}}_{\rm{f}}}$，其模长可以由滑动摩擦公式求得：

(9) ${F_{\rm{f}}} = \mu {F_z} = \mu \left| {{{{F}}_0}} \right|.$

式中： $\mu $为2种材料的滑动摩擦系数.

滑动摩擦力的方向与滑移方向相反，如图4所示. 将其沿X、Y方向分解后，可以得到

(10) $ \left. \begin{array}{l} {F_x} = - {F_{\rm{f}}}\cos \; \alpha ,\\ {F_y} = - {F_{\rm{f}}}\sin \; \alpha . \end{array} \right\} $

式中： $\alpha = {\rm{arctan}}\;({{\Delta {y_0}} / {\Delta {x_0}}})$.

结合式（8）~（10）可以得到滑动摩擦力：

(11) $ {{{F}}_{\rm{f}}} = \left[ {\begin{array}{*{20}{c}} {{F_x}}\\ {{F_y}}\\ 0 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - \mu \left| {{{{F}}_0}} \right|\cos \; \left(\arctan \; \dfrac{{{{{c}}_2}{{{F}}_0}}}{{{{{c}}_1}{{{F}}_0}}}\right)}\\ { - \mu \left| {{{{F}}_0}} \right|\sin \; \left(\arctan \; \dfrac{{{{{c}}_2}{{{F}}_0}}}{{{{{c}}_1}{{{F}}_0}}}\right)}\\ 0 \end{array}} \right]. $

因此，压脚压紧状态下的机器人末端作用力为

(12) $\begin{split} {{{F}}_{{\rm{total}}}} =& {{{F}}_0} + {{{F}}_{\rm{f}}} = {\left[ {{F_x},{F_y},{F_z}} \right]^{\rm{T}}}=\\ & \left[ {\begin{array}{*{20}{c}} { - \mu \left| {{{{F}}_0}} \right|\cos \; \left(\arctan \; \dfrac{{{{{c}}_2}{{{F}}_0}}}{{{{{c}}_1}{{{F}}_0}}}\right)}\\ { - \mu \left| {{{{F}}_0}} \right|\sin \; \left(\arctan \; \dfrac{{{{{c}}_2}{{{F}}_0}}}{{{{{c}}_1}{{{F}}_0}}}\right)}\\ {\left| {{{{F}}_0}} \right|} \end{array}} \right]. \end{split}$

由 ${{{F}}_0} = {\left[ {0,0,{F_z}} \right]^{\rm{T}}}$可知， ${{{c}}_1}{{{F}}_0} = {c_{13}}{F_z}$， ${{{c}}_2}{{{F}}_0} = {c_{23}}{F_z}$，因此式（12）可以简化为

(13) $\begin{split}& {{{F}}_{{\rm{total}}}} = \left[ {\begin{array}{*{20}{c}} { - \mu \cos \; \left(\arctan \; \dfrac{{{c_{23}}}}{{{c_{13}}}}\right)}\\ { - \mu \sin \; \left(\arctan \; \dfrac{{{c_{23}}}}{{{c_{13}}}}\right)}\\ 1 \end{array}} \right] \;\; \left| {{{{F}}_0}} \right|{\rm{ = }}\\&\left[ {\begin{array}{*{20}{c}} 1&0&{ - \mu \cos \; \left(\arctan \; \dfrac{{{c_{23}}}}{{{c_{13}}}}\right)}\\ 0&1&{ - \mu \sin \; \left(\arctan \; \dfrac{{{c_{23}}}}{{{c_{13}}}}\right)}\\ 0&0&1 \end{array}} \right]{{{F}}_0}. \end{split}$

综合式（8）、（13），可以得到压脚约束下的机器人刚度模型为

(14) $ \Delta {{{X}}_{\rm{t}}}{\rm{ = }}{{{C}}_{{\rm{tt}}}}{{{F}}_{{\rm{total}}}}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}&{ - \mu \cos \; \left(\arctan \; \dfrac{{{c_{23}}}}{{{c_{13}}}}\right){c_{11}} - \mu \sin \; \left(\arctan \; \dfrac{{{c_{23}}}}{{{c_{13}}}}\right){c_{12}} + {c_{13}}}\\ {{c_{21}}}&{{c_{22}}}&{ - \mu \cos \; \left(\arctan \; \dfrac{{{c_{23}}}}{{{c_{13}}}}\right){c_{21}} - \mu \sin \; \left(\arctan \; \dfrac{{{c_{23}}}}{{{c_{13}}}}\right){c_{22}} + {c_{23}}}\\ {{c_{31}}}&{{c_{32}}}&{ - \mu \cos \; \left(\arctan \; \dfrac{{{c_{23}}}}{{{c_{13}}}}\right){c_{31}} - \mu \sin \; \left(\arctan \; \dfrac{{{c_{23}}}}{{{c_{13}}}}\right){c_{32}} + {c_{33}}} \end{array}} \right]{{{F}}_0}{\rm{ = }}{{D}}{{{F}}_0}. $

式中： ${{D}}$为 $3 \times 3$的压脚约束下的平移柔度矩阵.

至此，压脚约束下的机器人刚度模型建立完成. 由式（13）可知，压脚作用力随机器人位姿变化而变化，因此利用传统机器人刚度模型预测不同位姿下的变形须测得所有位姿下的压脚作用力，无法在工程中实际应用. 利用式（14）所表示的机器人刚度模型只须测得压脚压紧力 ${{{F}}_0}$（取决于压脚气缸气压）便可预测任意位姿下的机器人末端变形，为后续机器人滑移变形补偿奠定理论基础.

由第2章的机器人刚度模型可知，压脚约束下的机器人平移柔度矩阵 ${{D}}$取决于机器人各个关节刚度和雅克比矩阵. 其中，雅克比矩阵由机器人的位姿决定，可以通过基于机器人运动学模型的微分变换法计算得到；关节刚度由机器人的驱动系统等决定，可以看作一组定值. 因此，可以基于压脚约束下的机器人刚度模型，通过实验辨识机器人的关节刚度.

由压脚约束下的机器人刚度模型可以得到

(15) $\Delta {{{X}}_{\rm{t}}} = {{D}}{{{F}}_0} = f({{q}},{k_i},{{{F}}_0}).$

式中： ${{q}} = {\left[ {{\theta _1},{\theta _2},{\theta _3},{\theta _4},{\theta _5},{\theta _6}} \right]^{\rm{T}}}$为机器人关节角； ${k_i}$为各关节刚度值， $i = 1,2,\cdots,6$.

由式（15）可知，机器人在压脚压紧力作用下产生的变形 $\Delta {{{X}}_{\rm{t}}}$与各关节刚度 ${k_i}$呈非线性关系，可由非线性最小二乘法求得各关节刚度. L-M算法因综合最速下降法和高斯牛顿法的优点，具有更快的求解速度和更好的全局收敛性，被广泛运用于求解非线性最小二乘问题^[23].

在不同机器人位姿 ${{q}}$下对机器人末端施加不同压紧力 ${{{F}}_0}$，测得机器人的实际平移变形 $\Delta {{{X}}_{\rm{a}}}$，同时由式（14）计算出平移变形的估计值 $\Delta {{{X}}_{\rm{e}}}$，将 $\Delta {{{X}}_{\rm{a}}}$和 $\Delta {{{X}}_{\rm{e}}}$差值的平方和最小作为目标函数，表达式为

(16) $\min \; S({{k}}) = \sum\limits_{i = 1}^n {{{{r}}_i}} {({{k}})^{\rm{T}}}{{{r}}_i}({{k}}) = {{R}}{({{k}})^{\rm{T}}}{{R}}({{k}}).$

式中： ${{k}}$为优化变量， ${{k}} = \left[ {{k_1},{k_2},{k_3},{k_4},{k_5},{k_6}} \right]$； ${{{r}}_i}({{k}})$为第 $i$组数据下实际变形和估算变形的差值，是 $3 \times 1$的向量， ${{{r}}_i}({{k}}){\rm{ = (}}\Delta {{{X}}_{\rm{a}}} - \Delta {{{X}}_{\rm{e}}}{)_i}$， $\Delta {{{X}}_{\rm{a}}} = {\left[ {\Delta {x_{\rm{a}}},\Delta {y_{\rm{a}}},\Delta {z_{\rm{a}}}} \right]^{\rm{T}}}$， $\Delta {{{X}}_{\rm{e}}} = {\left[ {\Delta {x_{\rm{e}}},\Delta {y_{\rm{e}}},\Delta {z_{\rm{e}}}} \right]^{\rm{T}}}$； ${{R}}({{k}})$为 $n$组数据的实际变形和估算变形的差值，是 $3n \times 1$的向量， ${{R}}({{k}}) = {\left[ {{{{r}}_1}({{k}}),{{{r}}_2}({{k}}),\cdots,{{{r}}_n}({{k}})} \right]^{\rm{T}}}$.

L-M算法求解关节刚度的算法流程如图5所示，具体步骤如下：1）选定初值 ${{{k}}^0}$，选取阻尼系数 $u$和放大系数 $v$，满足 $uv > 1$，选取精度 $\varepsilon > 0$，参数 $\beta \in (0,1.0)$，置 $j = 0$；2）计算函数 ${{R}}({{{k}}^j})$， $S({{{k}}^j})$；3）计算梯度 $\nabla {{R}}({{{k}}^j})$， $\nabla S({{{k}}^j})$，其中 $\nabla S({{{k}}^j}) = {(\nabla {{R}}({{{k}}^j}))^{\rm{T}}}{{R}}({{{k}}^j})$；4）令 ${{Q}} = {(\nabla {{R}}({{{k}}^j}))^{\rm{T}}}\nabla {{R}}({{{k}}^j})$，解方程 $({{Q}} + u{{I}})\Delta {{k}} = - \nabla S({{{k}}^j})$，得到迭代量 $\Delta {{k}}$，其中 ${{I}}$为单位矩阵；5）令 ${{{k}}^{j{\rm{ + }}1}} = {{{k}}^j} + {(\Delta {{k}})^{\rm{T}}}$，判断计算终止条件 $\left\| {\Delta {{k}}} \right\| < \varepsilon $是否满足，满足则 ${{{k}}^{j{\rm{ + }}1}}$即为最优解，不满足转步骤6）；6）若 $S({{{k}}^{j{\rm{ + }}1}}) < S({{{k}}^j}) + \beta {(\nabla {{S}}({{{k}}^j}))^{\rm{T}}}\Delta {{k}}$，说明 ${{k}}$的迭代方向正确，令 $u = u/v$，转步骤7），否则令 $u = uv$，转步骤4）；7）令 $j = j + 1$，转步骤2）.

如图6所示为关节刚度辨识实验系统，由KUKA KR600机器人、工件、带有气缸压脚系统的末端执行器、Kistler 9257B测力仪和Leica LTD640激光跟踪仪等组成，其中工件的材料为铝合金，压脚的材料为45钢.

在实验时，操作机器人将末端执行器定位到目标位置，并通过压脚上的4个激光位移传感器反馈调节机器人末端姿态，确保压脚垂直于工件表面. 压脚在气缸的作用下压紧工件表面，借助测力仪和激光跟踪仪测量机器人末端受力和变形. 逐渐改变压脚压力、变换机器人位姿进行多组实验. 本实验共选取5种不同的机器人位姿，在每种位姿下，通过调节气缸气压对机器人逐渐施加5次不同的压紧力. 各位姿下机器人关节角度和不同气缸气压P下的压脚压紧力分别如表2、3所示.

在由刚度辨识实验获得25组实验数据后，利用L-M算法求解KUKA KR600机器人的关节刚度. 1）根据式（16）编写目标函数的Matlab程序，程序中包含实验所用的机器人位姿、压脚压紧力、机器人末端变形的计算值和实测值等；2）确定优化变量 ${{{k}}^0}$的初值，合理的初值有助于提高算法的收敛速度与计算精度，根据Dumas等^[16]对KUKA KR240-2机器人关节刚度辨识的结果，将初值定为 ${{{k}}^0} = \left[ 1.0 \times {{10}^{10}},6.0 \times {{10}^9},6.0 \times {{10}^9},3.5 \times {{10}^9} \right. \left. 3.5 \times {{10}^9}, 3.5 \times {{10}^9} \right]$ N·mm/rad；3）在Matlab中调用L-M算法进行求解，得到优化变量的最优解：

$ \begin{split} {{{k}}_{{\rm{opt}}}} =& \left[ {1.324 \times {{10}^{10}},1.028 \times {{10}^{10}},8.109 \times {{10}^9},} \right.\;\;\;\;\;\;\;\;\; \\ & \left. {4.820 \times {{10}^9},4.603 \times {{10}^9},4.601 \times {{10}^9}} \right]\;{\rm{N \cdot mm/rad}}. \end{split} $

因此，机器人的各关节刚度k₁、k₂、k₃、k₄、k₅、k₆分别为1.324×10¹⁰、1.028×10¹⁰、8.109×10⁹、4.820×10⁹、4.603×10⁹、4.601×10⁹ N·mm/rad.

在获得KUKA KR600机器人关节刚度后，通过式（14）计算得到任意位姿下的机器人末端理论变形，从而对机器人制孔中的压脚滑移现象进行离线补偿. 通过关节刚度验证实验和压脚滑移变形补偿实验证明该补偿方案能有效提高机器人制孔定位精度.

压脚滑移变形补偿原理如图7所示. 图中，点A为理论位置，在未补偿情况下，机器人末端受压紧力作用滑移到点B处，即补偿前位置误差为 $\sqrt {{{(\Delta {x_{\rm{a}}})}^2} + {{(\Delta {y_{\rm{a}}})}^2}} $；点C为补偿后的位置，由理论位置减去预测的变形所得. 在此位置，机器人末端受压紧力作用滑移到点D，点A、D之间的距离即为补偿后的位置误差 $\sqrt {{{({x_1} - {x_0})}^2} + {{({y_1} - {y_0})}^2}} $.

压脚滑移变形补偿实验以铝合金壁板为工件，如图8所示. 在进行滑移补偿实验前，先对关节刚度的准确性进行验证. 选取压脚的气缸气压为0.5 MPa，选取6组机器人位姿，并利用激光跟踪仪测得6组位姿下的末端实际滑移，同时根据式（14）计算得到6组位姿下机器人理论末端滑移，如表4所示. 将不同位姿下机器人实际滑移和理论滑移制成图表，如图9所示. 图中，L为末端滑移，Q为不同组别. 可以看出，实际值与理论值较符合，最大偏差为0.04 mm，平均相对误差为11.70%，从而定量验证了该关节刚度的正确性，可以使用该组数值进行后续分析研究.

在验证完关节刚度的正确性后，参考上述6组位姿下的机器人理论末端滑移变形，分别对6个点的理论坐标进行补偿，而后操作机器人定位到补偿后的坐标下，在压紧压脚后，机器人末端受力滑移，测得此时的机器人末端位置，并与理论坐标相比，从而得到补偿后的位置误差，如表5所示. 表中，（x₀，y₀，z₀）为理论位置，（x₁，y₁，z₁）为补偿后的最终位置，[(x₁−x₀)²+(y₁−y₀)²]^1/2为补偿后位置误差.

将上述6组位姿下的机器人滑移变形补偿前、后的位置误差进行对比，如图10所示. 图中，δ为位置误差. 可以发现，采用本研究所述方法补偿机器人的滑移变形后，平均位置误差由原先的0.22 mm降低到0.05 mm.

（1）建立压脚约束下的机器人刚度模型，综合考虑压脚压紧力和摩擦力对机器人变形的影响，用于对机器人滑移变形进行预测和补偿.

（2）提出基于L-M算法的机器人关节刚度辨识方法，并采用测力仪测得的压脚作用力作为机器人所受外载，简化关节刚度辨识的实验装置.

（3）基于压脚约束下的机器人刚度模型，提出机器人滑移变形预测和补偿方法，使机器人制孔的平均位置误差由0.22 mm降为0.05 mm，提高机器人末端定位精度.

（4）在进行机器人关节刚度辨识实验时，由于测力仪面积较小且机器人相对工装无法移动，实验所选的机器人位姿较接近，由此计算的关节刚度可能存在一定误差，后续研究可以通过在机器人底部安装导轨或AGV车的方式，扩大实验中机器人的位姿变化.