基于改进傅里叶级数轨迹的运动学参数辨识

Kinematic Parameter Identification Based on Improved Fourier Series Trajectory

摘要

针对高端制造领域焊接机器人焊接效率低、精度差、运行不平稳等多种问题，提出一种基于改进傅里叶级数轨迹的焊接机器人运动学参数辨识方法. 首先，依据MDH（modified denavit-hartenberg）准则建立焊接机器人的运动学模型，并基于微分运动学建立焊接机器人的误差模型；随后，通过模式搜索算法找到最优的改进傅里叶级数轨迹，从中选出50个条件数最小的位姿点作为实验位姿集；最后，提出一种递归最小二乘方-差分进化混合（recursive least squares and differential evolution hybrid，RLS-DEH）算法，能够有效提高参数辨识的精度. 以库卡（KUKA）的焊接机器人作为实验对象，验证上述方法的可行性.实验结果表明，将使用改进傅里叶级数轨迹方法得到的最优位姿集代入误差模型后，经过RLS-DEH算法辨识，焊接机器人的绝对位置平均误差从1 mm降低到0.05 mm以内，相比标定前精度提高了95%，证明了所提方法的高效性与实用性.

Abstract

Aiming at the problems of low welding efficiency， poor accuracy and unstable operation of welding robots in high-end manufacturing field， a method of kinematics parameter identification of welding robots based on improved Fourier series trajectory was proposed. Firstly， a kinematic model of the welding robot was established based on the MDH （modified denavit-hartenberg） criterion， and an error model of the welding robot was established based on differential kinematics. Secondly， the optimal improved Fourier series trajectory was found through pattern search algorithm， and 50 pose points with the smallest number of conditions were selected as the experimental pose set. Finally， a RLS-DEH （recursive least squares and differential evolution hybrid） algorithm was proposed， which can effectively improve the accuracy of parameter identification. The above method was validated using KUKA’s welding robot as the experimental object. The experimental results showed that after， substituting the optimal pose set obtained using the improved Fourier series trajectory method into the error model and identifying it using the RLS-DEH algorithm， the average absolute position error of the welding robot was reduced from 1 mm to within 0.05 mm， which improved the accuracy by 95% compared with the previous calibration， proving the efficiency and practicality of the proposed method.

Graphical abstract

关键词

运动学标定 / 改进傅里叶级数轨迹 / RLS-DEH / 焊接机器人

随着我国工业水平的快速提高，焊接机器人逐步进入工业市场并扮演着越来越重要的角色^［1］. 焊接机器人是主要从事焊接工作的工业机器人，一般由工业机器人在末端轴加装焊枪组成^［2］，可以替代人工在恶劣工作环境下进行高质量的焊接工作. 与国外相比，我国焊接机器人研究起步相对较晚，焊接精度相对较低，因此提高焊接机器人的绝对定位精度显得尤为重要^［3-5］.

研究表明，影响焊接机器人绝对定位精度的主要因素有运动学参数误差、关节参数误差、刚度参数误差以及装配误差等，其中影响最大的是运动学参数误差^［6］. 机器人标定是降低运动学参数误差最常用的方法之一，一般包括运动学模型建立、位姿测量、运动学参数辨识、误差补偿^［7］. 运动学建模的基础是建立在DH（denavit-hartenberg，DH）准则^［8］之上，该准则只用4个参数就可以确定连杆坐标系，但如果相邻两轴平行或者接近平行时运动学参数会发生突变，为了解决这一问题，Cao等^［9］增加了一个绕y轴旋转的参数，改进了DH准则，提出了MDH准则. 为了确保测量效率，保证后期辨识的准确率，位姿测量主要通过借助高精密仪器完成. 目前市场的测量设备主要有激光跟踪仪、自动经纬仪、三坐标测量仪、球杆仪、拉线传感器等^［10-12］. 激光跟踪仪是测量位姿的主流设备，具有精度高、可测范围大、动态跟踪、占地小的优点，但其操作较为复杂，价格昂贵. 刘保军^［13］提出了一种基于激光跟踪仪的机器人运动学参数补偿的方法，为后续更深入研究机器人奠定了基础. 王殿君等^［14］使用激光跟踪仪对机器人进行运动学标定，最终将机器人的绝对定位精度提升了69.6%. 运动学参数辨识是机器人标定中重要的环节，其中最经典的方法是最小二乘法辨识，但此方法容易造成矩阵的奇异导致无法求解. 高贯斌等^［15］克服了最小二乘法辨识噪声大、发散慢的缺点，提出一种卡尔曼滤波算法，实现了运动学参数的快速收敛，将绝对定位精度由1.039 mm提高到了0.226 mm. Huang等^［16］提出了一种混合非线性参数辨识的LM（levenberg-marquarat）和递推最小二乘法的混合辨识算法并进行了仿真试验，证明了算法的稳定性. Fang等^［17］提出了基于粒子分类和按维动态改变权重的改进粒子群算法，有效避免陷入局部最优，提高了算法的精度和收敛速度，具有更强的优化性能. Schnierle等^［18］融合了遗传算法和禁忌算法使得该算法既有全局优化能力又跳出了局部最优解的困局，实现了模型误差的精确辨识. 误差补偿^［19-20］即将辨识出的运动学参数误差补偿到运动学参数中. 梁兆东等^［21］设计了一种误差补偿的运动学算法，将绝对位置误差的均值从0.575 4 mm降低到0.277 9 mm. 高贯斌等^［15］提出了一种关节空间插值误差补偿方法，建立网格形式的误差补偿数据库，并通过实验验证了方法的可行性.

本文主要对位姿测量中选取位姿点和运动学参数辨识中辨识算法进行研究，参照动力学参数辨识^［22-24］中最优激励设计方法，提出一种五次多项式替代傅里叶轨迹中常数项的改进傅里叶级数轨迹^［25-26］，由误差模型中雅可比矩阵的条件数作为优化函数，在确定了最优轨迹后从中选出50个条件数最小的位姿点作为实验的位姿集，即为最优位姿集. 本文提出的RLS-DEH算法提高了参数辨识的精度，最后通过实验验证了本文所提方法的可行性.

1 运动学建模

1.1 MDH准则建模

在工业机器人的运动学建模中，DH准则是最经典的建模方式. 然而DH准则遇到相邻两轴线平行或者近似平行时会出现奇异问题. 为了解决这一问题， Cao等^［9］对DH准则的奇异问题进行了简化处理，使两相邻平行轴平行即忽略平行轴之间的相对误差，提出了MDH准则，并成功标定了工业机器人的参数误差^［27］. 本文所研究的焊接机器人采用MDH准则建立运动学模型，其相邻连杆坐标系如图1所示.

单个连杆的齐次变换矩阵通式可以表示为：

{\mathit{A}}_i=\left[\begin{array}{cccc}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i& -\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i& \mathrm{0}& a_{i-\mathrm{1}}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{i-\mathrm{1}}& \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{i-\mathrm{1}}& -\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{i-\mathrm{1}}& -d_i\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{i-\mathrm{1}}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{i-\mathrm{1}}& \mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{i-\mathrm{1}}& \mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{i-\mathrm{1}}& d_i\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{i-\mathrm{1}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（1）

式中：θ_i 表示相邻坐标系之间x轴的夹角；α_i-1表示相邻坐标系之间z轴的夹角；d_i 表示相邻坐标系之间x轴的距离；a_i-1表示相邻坐标系之间z轴的距离.

本文所用KUKA焊接机器人的DH参数和三维模型分别如表1和图2所示.

基于DH参数，由正向运动学得到总的齐次变换矩阵为：

{\mathit{A}}_{\mathrm{6}}^{\mathrm{0}}=A_{\mathrm{1}}{A}_{\mathrm{2}}{A}_{\mathrm{3}}{A}_{\mathrm{4}}{A}_{\mathrm{5}}{A}_{\mathrm{6}}=\left[\begin{array}{cccc}{n}_x& o_x& a_x& p_x\\ n_y& o_y& a_y& p_y\\ n_z& o_z& a_z& p_z\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（2）

1.2 误差建模

由焊接机器人的正向运动学方程可知，末端位姿由4个连杆参数计算得到. 假设焊接机器人末端位姿的微小误差是由连杆参数的微小误差导致，得：

\mathrm{\Delta }{\mathit{A}}_i=\frac{\partial {\mathit{A}}_i}{\partial a_{i-\mathrm{1}}}\mathrm{\Delta }{a}_{i-\mathrm{1}}+\frac{\partial {\mathit{A}}_i}{\partial {\alpha }_{i-\mathrm{1}}}\mathrm{\Delta }{\alpha }_{i-\mathrm{1}}+\frac{\partial {\mathit{A}}_i}{\partial d_i}\mathrm{\Delta }{d}_i+\frac{\partial {\mathit{A}}_i}{\partial {\theta }_i}\mathrm{\Delta }{\theta }_i=A_i\delta {\mathit{A}}_i（3）

式中：\mathrm{\Delta }{\mathit{A}}_i为矩阵的偏移误差；Δa_i-1、Δα_i-1、Δd_i 、Δθ_i 分别为相应的连杆参数误差；δA_i 为微分矩阵.

结合式（1），可得：

\frac{\partial {\mathit{A}}_i}{\partial {\theta }_i}=\left[\begin{array}{cccc}-s{\theta }_i& -c{\theta }_i& \mathrm{0}& \mathrm{0}\\ c{\theta }_{i}c{\alpha }_{i-\mathrm{1}}& -s{\theta }_{i}c{\alpha }_{i-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ c{\theta }_{i}s{\alpha }_{i-\mathrm{1}}& s{\theta }_{i}s{\alpha }_{i-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]（4）

将式（4）转换成：

\frac{\partial {\mathit{A}}_i}{\partial {\theta }_i}={\mathit{A}}_i{V}_{\theta }（5）

求得：

V_{\theta }={\mathit{A}}_i^{-\mathrm{1}}\frac{\partial {\mathit{A}}_i}{\partial {\theta }_i}=\left[\begin{array}{cccc}\mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]（6）

其中 A_i^-1为：

{\mathit{A}}_i^{-\mathrm{1}}=\left[\begin{array}{cccc}c{\theta }_i& s{\theta }_{i}c{\alpha }_{i-\mathrm{1}}& s{\theta }_{i}s{\alpha }_{i-\mathrm{1}}& \mathrm{0}\\ -s{\theta }_i& c{\theta }_{i}c{\alpha }_{i-\mathrm{1}}& c{\theta }_{i}s{\alpha }_{i-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& -s{\alpha }_{i-\mathrm{1}}& c{\alpha }_{i-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]（7）

同理，可以得到V_a 、V_d 、V_α，分别将其代入式（3），可得微分矩阵δA_i .

\delta {\mathit{A}}_i=V_{\theta }\mathrm{\Delta }{\theta }_i+V_d\mathrm{\Delta }{d}_i+V_{\alpha }\mathrm{\Delta }{\alpha }_{i-\mathrm{1}}+V_a\mathrm{\Delta }{a}_{i-\mathrm{1}}（8）

将每个轴的偏移误差相加可以得到焊接机器人的末端位姿误差，展开并忽略二阶及高阶项可得：

\mathrm{\Delta }\mathit{A}\approx \mathit{A}{\displaystyle \sum _{i=\mathrm{1}}^I}(A_{i+\mathrm{1}}\cdots A_I{)}^{-\mathrm{1}}\delta {\mathit{A}}_i(A_{i+\mathrm{1}}\cdots A_I)=\mathit{A}\delta \mathit{A}（9）

用 B_i+1表示式（9）连乘矩阵：

A_{i+\mathrm{1}}\dots A_I={\mathit{B}}_{i+\mathrm{1}}=\left[\begin{array}{cccc}B{n}_{i+\mathrm{1}}& B{o}_{i+\mathrm{1}}& B{a}_{i+\mathrm{1}}& B{p}_{i+\mathrm{1}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（10）

最终可以得到焊接机器人各个方向的位置分量：

\begin{array}{l}\mathrm{d}{x}_I={\displaystyle \sum _{i=\mathrm{1}}^I}\left[{\vartheta }_i\cdot (B{p}_{i+\mathrm{1}}\times B{n}_{i+\mathrm{1}})+d_i\cdot B{n}_{i+\mathrm{1}}\right]\\ \mathrm{d}{y}_I={\displaystyle \sum _{i=\mathrm{1}}^I}\left[{\vartheta }_i\cdot (B{p}_{i+\mathrm{1}}\times B{o}_{i+\mathrm{1}})+d_i\cdot B{o}_{i+\mathrm{1}}\right]\\ \mathrm{d}{z}_I={\displaystyle \sum _{i=\mathrm{1}}^I}\left[{\vartheta }_i\cdot (B{p}_{i+\mathrm{1}}\times B{a}_{i+\mathrm{1}})+d_i\cdot B{a}_{i+\mathrm{1}}\right]\end{array}（11）

式中：{\vartheta }_i和d_i 分别为微分旋转分量与微分平移分量.

{\vartheta }_i=\left[\begin{array}{c}c{\theta }_i\mathrm{\Delta }{\alpha }_{i-\mathrm{1}}\\ -s{\theta }_i\mathrm{\Delta }{\alpha }_{i-\mathrm{1}}\\ \mathrm{\Delta }{\theta }_i\end{array}\right]（12）

d_i=\left[\begin{array}{c}c{\theta }_i\mathrm{\Delta }{\alpha }_{i-\mathrm{1}}-d_{i}s{\theta }_i\mathrm{\Delta }{\alpha }_{i-\mathrm{1}}\\ -s{\theta }_i\mathrm{\Delta }{\alpha }_{i-\mathrm{1}}-d_{i}c{\theta }_i\mathrm{\Delta }{\alpha }_{i-\mathrm{1}}\\ \mathrm{\Delta }{d}_i\end{array}\right]（13）

结合式（9）和式（11）可得焊接机器人末端位置误差与机器人运动学参数误差之间的线性关系为：

\mathrm{\Delta }P=\mathit{J}(q)\mathit{\epsilon }（14）

式中：q为关节变量； J（q）为参数误差的雅可比矩阵；ΔP为焊接机器人的末端误差； ε 为运动学参数误差矩阵， ε=［Δa₀ Δα₀ Δd₁ Δθ₁… Δa_I-1 Δα_I-1 Δd_I Δθ_I ］^-1.

2 基于改进傅里叶级数轨迹求位姿集

选取最优位姿集是标定工业机器人的重要环节.传统方法是在机器人活动的最大空间内随机选取50个位姿点作为实验位姿集. 但这种方法找到的位姿集并不是最优位姿集，可能会出现矩阵条件数过大、矩阵奇异的情况，进而影响运动学参数辨识的精度. 仿照动力学参数辨识中的激励轨迹设计，本文提出了改进傅里叶级数轨迹求最优位姿集的方法，即用五次多项式替代原傅里叶轨迹中的常数项. 具体表达形式如下：

q_i(t)={\displaystyle \sum _{k=\mathrm{1}}^N}\frac{a_k^i}{w_{\mathrm{f}}k}sin(w_{f}kt)-

{\displaystyle \sum _{k=\mathrm{1}}^N}\frac{b_k^i}{w_{\mathrm{f}}k}cos(w_{f}kt)+{\displaystyle \sum _{k=\mathrm{0}}^{\mathrm{5}}}{c}_k^i{t}^k

\stackrel{}{{\dot{q}}_i(t)}={\displaystyle \sum _{k=\mathrm{1}}^N}{a}_k^{i}cos(w_{f}kt)+

{\displaystyle \sum _{k=\mathrm{1}}^N}{b}_k^{i}sin(w_{f}kt)+{\displaystyle \sum _{k=\mathrm{1}}^{\mathrm{5}}}{c}_k^{i}k{t}^{k-\mathrm{1}}

{\ddot{q}}_i(t)=-{\displaystyle \sum _{k=\mathrm{1}}^N}{a}_k^i{w}_{f}ksin(w_{f}kt)+

{\displaystyle \sum _{k=\mathrm{1}}^N}{b}_k^i{w}_{f}kcos(w_{f}kt)+{\displaystyle \sum _{k=\mathrm{2}}^{\mathrm{5}}}{c}_k^{i}k(k-\mathrm{1})t^{k-\mathrm{2}}（15）

式中： t表示时间；q表示角度；w_f为傅里叶轨迹基频角频率；N为谐波阶数；a_k^i、b_k^i为待优化参数；c_k^i可由a_k^i、b_k^i表示.

将改进傅里叶级数轨迹作为寻找位姿点轨迹的表达式，轨迹参数a_k^i、b_k^i的取值应该受到各个关节的极限位置、速度、加速度的限制，由此可以得到该非线性优化问题的约束条件为：

\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n} {\mathit{J}}_{\mathrm{1}}=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\left(\mathit{J}\right)\\ \mathrm{s}.\mathrm{t}. \left\{\begin{array}{l}s\left(q\left(t\right)\right) \in S \\ q_{\mathrm{m}\mathrm{i}\mathrm{n}}\le q\left(t\right)\le q_{\mathrm{m}\mathrm{a}\mathrm{x}}\\ {\dot{q}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le \dot{q}\left(t\right)\le {\dot{q}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\\ {\ddot{q}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le \ddot{q}\left(t\right)\le {\ddot{q}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}\right.\end{array}（16）

式中：J₁为目标函数；cond（J）表示雅可比误差矩阵的条件数；q（t）表示机器人各个关节的角度； \dot{q}\left(t\right)表示角速度； \ddot{q}\left(t\right)表示角加速度；q_{\mathrm{m}\mathrm{i}\mathrm{n}}、 {\dot{q}}_{\mathrm{m}\mathrm{i}\mathrm{n}}、 {\ddot{q}}_{\mathrm{m}\mathrm{i}\mathrm{n}}分别表示角度、角速度、角加速度的下限；q_{\mathrm{m}\mathrm{a}\mathrm{x}}、{\dot{q}}_{\mathrm{m}\mathrm{a}\mathrm{x}}、{\ddot{q}}_{\mathrm{m}\mathrm{a}\mathrm{x}} 分别表示角度、角速度、角加速度的上限； s\left(q\left(t\right)\right)\in S表示关节角度在机器人的整个运动空间内.

从带有参数的改进傅里叶级数轨迹中选出50个点代入误差模型得到雅可比矩阵 J，则目标函数设置为寻找条件数最小的雅可比矩阵 J₁. 由于目标函数比较复杂，用传统的求解方法非常困难. 模式搜索算法是一种高效的智能寻优算法，这种算法具有运算速度快、计算参数少、容易上手等优点，所以本文使用模式搜索算法来求解上述约束优化问题. 根据经验值将模式搜索算法参数代入算法模型，可得到轨迹参数进而确定轨迹. 从确定好的轨迹中依据时间间隔选出200个位姿点，再从这200个位姿点依次选出50个位姿点进行迭代，将每50个位姿点代入误差模型算出条件数，最终找到条件数最小的50个位姿点，进而确定最优位姿集. 改进傅里叶级数轨迹求最优位姿集的流程图如图3所示.

在最终的改进傅里叶级数轨迹中，轨迹的运动周期T=25 s，基频w_{\mathrm{f}}=\mathrm{2}\mathrm{\pi }/T=\mathrm{0.08}\mathrm{\pi }，谐波项取N=5，经过模式搜索算法不断迭代寻优得到最终轨迹，如图4所示.

从图4中可以看出，角速度与角加速度的轨迹周期起止时刻都为0，这就说明了轨迹运行稳定，机械臂能够准确地沿着轨迹运动，有效避免了激光跟踪仪测量位姿点时产生的误差；同时轨迹中的位姿点在焊接机器人最大运动空间中都有分布，相比于传统方法在焊接机器人运动空间中随机寻找位姿点可能出现的局限性，该方法可以更好地降低误差模型中雅可比矩阵的条件数，提高辨识精度. 轨迹具有周期性，可以多次采集数据，提高实验数据的准确性.

3 基于RLS-DEH算法的运动学参数辨识

3.1 RLS-DEH优化算法

由第1节焊接机器人的误差模型可知，若已知焊接机器人末端定位误差和运动学参数雅可比矩阵，即可求出焊接机器人的运动学参数误差值. 依据第2节改进傅里叶级数轨迹中找到的最优位姿集，使用激光跟踪仪测得焊接机器人的实际末端位置，然后通过正向运动学模型可以计算得到焊接机器人的理论末端位置，即可获得焊接机器人末端定位误差. 运动学参数雅可比矩阵由带有误差的运动学参数依据误差模型计算得到，因此辨识运动学参数误差量这一问题可以转换成求解形如AX=B的非线性方程. 对于参数辨识中的非线性方程，一般采用最小二乘法进行求解，如递归最小二乘方（recursive least squares，RLS）算法. 另外对于参数辨识问题，也可以将其看成一个全局优化问题，然后利用智能优化算法进行求解. 差分进化（differential evolution，DE）算法是一种全局优化算法，与常见的遗传算法、粒子群算法等相似，都是受自然现象启发而发展出来的智能优化算法，可以应用于参数辨识问题.

RLS算法是最小二乘算法中的一类快速算法. RLS算法主要是在误差平方和最小原则的基础上提出一种解析拟合模型参数的迭代递推公式，可以实现在新的样本数据到来时，利用新的样本数据以及旧的最优模型参数便捷地计算新的最优模型参数，从而避免直接计算方法中的逆矩阵. 在第1节的误差模型中，将焊接机器人末端误差项展开并忽略二阶及高阶项，这样就把问题转化成了线性最小二乘问题. 应用到参数辨识中的公式可以写为：

\mathrm{\Delta }P=\mathit{J}\mathrm{\Delta }X（17）

式中：ΔP为焊接机器人末端位置误差，是RLS算法的输入；ΔX为运动学参数误差，是RLS算法的输出. 经过层层迭代，循环结束后即可得到次优的运动学参数误差值.

DE算法是在遗传算法的思想基础上提出的，属于一种基于群体差异的、多目标优化的算法，用于求解多维空间的最优解. DE算法相对于遗传算法（genetic algorithm，GA）而言，相同点是两者都包含种群初始化、变异、交叉、选择等步骤；不同点是DE算法变异向量是由父代差分向量生成，并与父代个体向量生成新的个体向量，且直接与其父代个体进行选择. 显然DE算法的搜索范围相对于GA算法更大且具有更广泛的适应性.

运用DE算法进行参数辨识的步骤如下.

步骤1：正确选择DE算法的参数并随机生成初始化种群. 在特定的搜索空间内随机生成N_P个个体X_i，g 以完成初始化种群.

X_{i,g}=\mathit{U}（18）

式中：i=1，2，3，…，N_P； U 是25×1的向量，每个分量值是所设置的上下界的随机值.

步骤2：差分变异操作.随机选取种群中一个个体进行差分向量缩放，与种群中相异个体相加得到变异向量.

V_{i,g}=X_{r\mathrm{1},g}+F(X_{r\mathrm{2},g}-X_{r\mathrm{3},g})（19）

式中： X_r1、X_r2、X_r3是3个随机个体，区间为［1， N_P］； F为变异因子，用于控制偏差的放大，主要影响算法的搜索步长，扩大搜索范围.

步骤3：对变异个体进行交叉，生成新的个体，交叉操作的方法为：

U_{i,j,g}=\left\{\begin{array}{l}{V}_{i,j,g}\mathrm{i}\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}(\mathrm{0,1})\le C_R\\ X_{i,j,g},\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.（20）

式中：C_R∈［0，1］为交叉因子.

步骤4：采用贪婪选择的策略，在目标向量与实验向量之间选择最优个体作为下一代种群个体.

（21）

在产生新一代种群之后，若满足终止条件则结束迭代过程，否则将继续执行变异、交叉、选择等操作. 迭代结束后获得的种群应该是评价函数值最小的个体，也就是最优的运动学参数误差值.

RLS算法和DE算法都可以辨识焊接机器人的运动学参数误差，但两者都存在一定的弊端. RLS算法在求解运动学参数误差时将非线性模型简化为线性模型，因此得到的只是一个近似解. DE算法是根据误差模型建立目标函数，拥有全局搜索能力，得到的结果更加精确，但DE算法的搜索范围庞大导致求解过程需要花很长时间. 为了兼顾效率和精度，本文提出一种更优的RLS-DEH算法，即将RLS算法得到的近似解作为DE算法种群初始化的中心值，合理地缩小了搜索范围，在保证DE算法辨识精度的同时提高辨识效率.

X_{i,g}=\mathrm{\Delta }X+U（22）

式（22）表示RLS-DEH算法的种群初始化.

算法改进示意图如图5所示，其中红色圆圈表示DE算法的搜索空间，红色的点表示DE算法种群初始化的中心值，由式（18）进行初始化. 蓝色圆圈表示RLS-DEH算法的搜索空间，绿色的点表示RLS-DEH算法种群初始化的中心值，由式（22）进行初始化. RLS-DEH算法相比于DE算法改进的本质在于，由RLS算法提供的近似解ΔX完成了种群初始化中心值从红色点到绿色点的转化，使初始化的中心值更加接近于最优值的附近，减小了搜索空间的范围，提高了辨识效率.

3.2 参数辨识实验

为了验证以上所提方法的可行性，本文采用KUKA焊接机器人、Leica激光跟踪仪和计算机进行标定实验. 焊接机器人标定实验平台如图6所示. 首先，利用改进的傅里叶级数轨迹方法得到焊接机器人的运动轨迹和最优位姿集，使用激光跟踪仪采集焊接机器人末端的实际位置，将采集到的数据代入误差模型并利用计算机辨识出运动学参数误差；然后完成对焊接机器人的误差补偿；最后采用激光跟踪仪验证补偿后的效果. 焊接机器人标定实验步骤如图7所示.

在参数辨识中采用RLS-DEH算法进行辨识. 同时为了验证所提方法的有效性，将RLS-DEH算法辨识后的定位误差和计算时效性分别与RLS算法、DE算法辨识后的定位误差和计算时效性进行比较.

图8为不同算法辨识结果分布图. 由图8可知，标定前焊接机器人X、Y、Z方向的平均位置误差均达到1 mm. 分别经过RLS算法、DE算法和RLS-DEH算法辨识后，平均位置误差分别降低到0.6 mm、 0.3 mm和0.05 mm以内. 辨识算法计算过程对比图如图9所示. 由图9可知，RLS算法、DE算法和RLS-DEH算法分别迭代300次、180次和50次时适应度值为0. 综上所述，相比传统算法，RLS-DEH算法辨识精度更高且时效性更好， 相比之前焊接机器人精度提高了95%.

3.3 补偿验证实验

由3.2节参数辨识实验可得，RLS-DEH算法的辨识效果最好，为了验证该算法的辨识结果在焊接机器人实际工作时的工作效果，分别设立了直线、斜线和空间补偿验证实验. 焊接机器人按照预定好的轨迹运动，经计算得到RLS-DEH算法误差补偿验证图如图10所示. 由图10可知，直线运动轨迹与斜线运动轨迹的平均位置误差在0.03 mm以内，空间运动轨迹的平均位置误差在0.05 mm以内，与标定前相比，焊接精度得到了大幅度提高，本文所提运动学标定方法达到实验预期效果.

4 结 论

本文主要对基于改进傅里叶级数轨迹的焊接机器人运动学参数辨识进行研究. 首先，建立焊接机器人的DH运动学模型，基于微分运动学方法构建焊接机器人的运动学误差模型；其次，参照动力学参数辨识中求取激励轨迹的方法得到焊接机器人的改进傅里叶级数轨迹，从该轨迹中筛选出用于参数辨识的最优位姿集；最后，提出了一种基于RLS-DEH算法对运动学参数进行辨识，通过实验证明了该混合优化算法拥有更高的辨识精度，焊接机器人的平均位置误差从1 mm降低到了0.05 mm以内. 综上，基于改进傅里叶级数轨迹的焊接机器人运动学模型具有更好的泛化能力.

参考文献

[1] 郭岩宝， 王斌， 王德国， 等 ．焊接机器人的研究进展与发展趋势［J］．现代制造工程，2021（5）：53-63．

[2] GUO Y B， WANG B， WANG D G，et al ．Research progress and development trend of welding robot［J］．Modern Manufacturing Engineering，2021（5）：53-63．（in Chinese）

[3] 黄通，刘志浩，高钦和，等 ．重载车辆轮胎模型参数辨识与灵敏度分析［J］．湖南大学学报（自然科学版），2022，49（8）：36-44．

[4] HUANG T， LIU Z H， GAO Q H，et al ．Parameter identification and sensitivity analysis of heavy-duty vehicle tire model［J］．Journal of Hunan University （Natural Sciences）， 2022， 49（8）：36-44．（in Chinese）

[5] 洪磊，杨小兰，王保升，等 ．基于球面均匀分布的焊接机器人TCP标定方法［J］．焊接学报，2020，41（8）：14-21．

[6] HONG L， YANG X L， WANG B S，et al ．TCP calibration method based on spherical uniform distribution for welding robot［J］．Transactions of the China Welding Institution，2020，41（8）：14-21．（in Chinese）

[7] 陈启军，周自强，刘成菊，等 ．基于视觉的六自由度机械臂运动学参数辨识［J］．同济大学学报（自然科学版），2020，48（7）：1023-1029．

[8] CHEN Q J， ZHOU Z Q， LIU C J，et al. Design of a vision calibration method of kinematic parameters for 6-DoF manipulator［J］.Journal of Tongji University （Natural Science）， 2020， 48（7）：1023-1029．（in Chinese）

[9] 孙炜，舒帅，李汶霖，等 ．基于阶梯形标定物的线激光传感器手眼标定算法［J］．湖南大学学报（自然科学版），2022，49（2）：83-92．

[10] SUN W， SHU S， LI W L，et al.Hand-eye calibration algorithm of line laser sensor based on step-shaped calibrators［J］.Journal of Hunan University （Natural Sciences），2022，49（2）：83-92．（in Chinese）

[11] 张颖，乔贵方，王保升，等 ．基于优化位姿集的工业机器人运动学参数辨识方法研究［J］．中国测试，2023，49（3）：91-95．

[12] ZHANG Y， QIAO G F， WANG B S，et al ．Research on kinematics parameter identification method of industrial robot based on optimized pose set［J］．China Measurement & Test，2023，49（3）：91-95．（in Chinese）

[13] 李楠，平雪良，王晨学，等 ．一种机器人的平面约束自标定技术研究［J］．机械科学与技术，2019，38（8）：1214-1222．

[14] LI N， PING X L， WANG C X，et al ．Exploring a planar constraints self-calibration technique of robots［J］．Mechanical Science and Technology for Aerospace Engineering，2019，38（8）：1214-1222．（in Chinese）

[15] DU G L， LIANG Y H， LI C Q，et al ．Online robot kinematic calibration using hybrid filter with multiple sensors［J］．IEEE Transactions on Instrumentation and Measurement，2020，69（9）：7092-7107．

[16] CAO H Q， NGUYEN H X， TRAN T N C，et al ．A robot calibration method using a neural network based on a butterfly and flower pollination algorithm［J］．IEEE Transactions on Industrial Electronics，2022，69（4）：3865-3875．

[17] ZHANG Y， GUO J W， LI X Y ．Study on redundancy in robot kinematic parameter identification［J］．IEEE Access，2022，10：60572-60584．

[18] CONTE J， SANTOLARIA J， MAJARENA A C，et al ．Modelling，kinematic parameter identification and sensitivity analysis of a Laser Tracker having the beam source in the rotating head［J］．Measurement，2016，89：261-272．

[19] CONTE J， MAJARENA A C， ACERO R，et al.Performance evaluation of laser tracker kinematic models and parameter identification［J］.The International Journal of Advanced Manufacturing Technology， 2015， 77（5/8）： 1353-1364.

[20] 刘保军 ．激光跟踪仪在机器人运动学参数补偿中的研究［J］．激光杂志，2016，37（7）：77-80．

[21] LIU B J.Study on laser tracker in robot kinematics parameter compensation［J］.Laser Journal，2016，37（7）：77-80．（in Chinese）

[22] 王殿君，王超星， 陈亚， 等 ．模块化六自由度机器人运动学标定与实验研究［J］．高技术通讯， 2016， 26（6）： 577-584．

[23] WANG D J， WANG C X， CHEN Y，et al ．Kinematic calibration and experimental study of a modular six-DOF robot［J］．Chinese High Technology Letters，2016，26（6）：577-584．（in Chinese）

[24] 高贯斌，张石文，那靖，等 ．基于标定和关节空间插值的工业机器人轨迹误差补偿［J］．机械工程学报，2021，57（21）：55-67．

[25] GAO G B， ZHANG S W， NA J，et al ．Compensation of trajectory error for industrial robots by interpolation and calibration method in the joint space［J］．Journal of Mechanical Engineering，2021，57（21）：55-67．（in Chinese）

[26] HUANG Y J， KE J H， ZHANG X M，et al ．Dynamic parameter identification of serial robots using a hybrid approach［J］．IEEE Transactions on Robotics，2023，39（2）：1607-1621．

[27] FANG L J， DANG P F. A step identification method of joint parame-ters of robots based on the measured pose of end-effector［J］. Proceedings of the Institution of Mechanical Engineers， Part C：Journal of Mechanical Engineering Science，2015，229（17）：3218-3233．

[28] SCHNIERLE M， RÖCK S. Latency and sampling compensation in mixed-reality-in-the-loop simulations of production systems［J］.Production Engineering，2023，17（3/4）：341-353．

[29] QI J D， CHEN B， ZHANG D H. Compensation for absolute positioning error of industrial robot considering the optimized measurement space［J］. International Journal of Advanced Robotic Systems， 2020， 17（2）： 172988142092164．

[30] 焦嘉琛， 田威， 张霖，等. 工业机器人作业误差分级补偿技术［J］.计算机集成制造系统，2022，28（6）：1627-1637．

[31] JIAO J C， TIAN W， ZHANG L，et al ．Hierarchical compensation technology for machining error of industrial robots［J］. Computer Integrated Manufacturing Systems， 2022，28（6）： 1627-1637．（in Chinese）

[32] 梁兆东，肖永强，游玮，等 ．工业机器人误差补偿运动学算法设计［J］．机械设计与制造，2021（5）：273-275．

[33] LIANG Z D， XIAO Y Q， YOU W，et al ．Design for errors compensation kinematic algorithm of industrial robot［J］．Machinery Design & Manufacture，2021（5）：273-275．（in Chinese）

[34] 杨浩，项丹，冯树辉，等 ．基于直接积分最小二乘法的同步发电机的参数辨识［J］．湖南大学学报（自然科学版），2013，40（7）：69-73．

[35] YANG H， XIANG D， FENG S H，et al ．Parameter identification of synchronous generator based on direct integration least square algorithm［J］．Journal of Hunan University （Natural Sciences），2013， 40（7）： 69-73．（in Chinese）

[36] ARGIN O F， BAYRAKTAROGLU Z Y ．Consistent dynamic model identification of the Stäubli RX-160 industrial robot using convex optimization method［J］．Journal of Mechanical Science and Technology，2021，35（5）：2185-2195．

[37] JIN H Z， LIU Z X， ZHANG H， et al.A dynamic parameter identification method for flexible joints based on adaptive control［J］. IEEE/ASME Transactions on Mechatronics， 2018， 23（6）： 2896-2908．

[38] 曾奎， 耿云海，陈炳龙 ．基于傅里叶级数的小推力航天器快速轨迹设计［J］．自动化学报，2016，42（11）：1641-1647．

[39] ZENG K， GENG Y H， CHEN B L ．Rapid design of low-thrust rendezvous trajectory with Fourier series［J］．Acta Automatica Sinica， 2016， 42（11）： 1641-1647．（in Chinese）

[40] 梅江平，方志炜，贺莹，等 ．基于傅里叶级数的最优能量运动轨迹优化方法［J］．机械设计，2018，35（10）：1-7．

[41] MEI J P， FANG Z W， HE Y，et al ．Method of optimal-energy motion trajectory planning based on Fourier series［J］．Journal of Machine Design，2018，35（10）：1-7．（in Chinese）

[42] 陆艺，葛文琦，郭斌 ．基于标准球距离约束的工业机器人参数标定［J］．计量学报，2020，41（9）：1048-1054．

[43] LU Y， GE W Q， GUO B. Calibration of industrial robot parameters based on standard ball distance constraints［J］. Acta Metrologica Sinica，2020，41（9）：1048-1054．（in Chinese）

基金资助

国家重点研发计划资助项目（2022YFB3404804）; 国家自然科学基金资助项目（52375041）