在航空航天、船舶制造和能源化工等领域，常需采用电弧焊、激光焊和搅拌摩擦焊等方法实现复杂空间曲线轨迹构件的精密连接^[1-4]。对执行机构各关节运动参量的合理规划是实现高质量焊接的前提^[5]。

在复杂空间曲线轨迹焊接运动规划过程中，往往需满足以下焊接工艺要求：1)焊接速度可在焊前预设，且在焊接过程中保持恒定；2)在焊接过程中焊接姿态保持恒定，焊炬轴线与焊接区(熔池或搅拌区)法向可保持恒定的预设倾角，且焊接区与世界坐标系保持相对稳定的姿态，例如在电弧焊和激光焊中常需保持平焊位置，且焊炬轴线与熔池表面垂直或保持一定的倾角，又如在搅拌摩擦焊中常需使轴肩与工件表面成恒定角度以施加一定的顶锻压力；3)焊炬末端和焊接区距离可在焊前预设，且在焊接过程中保持恒定，在电弧焊中表现为弧长恒定，在激光焊中表现为离焦量恒定，在搅拌摩擦焊中表现为搅拌头插入深度恒定。

传统的焊接轨迹规划方法主要分为示教再现和离线编程两类，其基本原理是在工件表面选取一定数量的离散点，以运动时间^[6-7]或消耗能量^[8-9]等为优化准则，在关节坐标系或Descartes坐标系中进行运动规划，计算每个时刻各关节的运动参量，实现连续轨迹焊接。这种运动规划方法主要考虑机器人的运动学和动力学特征，难以保证焊接速度、焊接姿态、弧长、离焦量、搅拌头插入深度等焊接工艺需求。文[10-11]提出了一种用于立面内任意2维曲线焊接的3轴联动装置控制方法，该方法对轨迹上选取的离散点进行圆弧插补，在焊接过程中保证焊接速度、焊炬末端与焊接区距离保持恒定，且熔池始终位于水平位置，焊炬轴线始终垂直于熔池表面。但该技术只适用于平面曲线焊接运动规划，不能直接应用于复杂空间轨迹焊接，且不能适用于要求焊炬具有一定前倾角或后倾角的场合。

本文基于运动学建模的方法，针对焊接工艺对焊接速度、焊接姿态、焊炬末端与焊接区距离等的要求，提出一种适用于任意复杂空间轨迹焊接的运动规划方法，根据工件的计算机辅助设计(computer aided design, CAD)模型或实际测量数据自动计算每个时刻各关节位置、速度等运动参量，适用于复杂空间轨迹构件的自动焊接场合。

1 空间曲线轨迹样条插补方法

在进行运动规划前，必须首先获得待焊轨迹上若干个离散点的信息。记与待焊工件固结的工件坐标系为{P}，沿待焊轨迹依次选取N个离散点，测量其在工件坐标系{P}中的3维坐标X₁, X₂, …, X_N以及各点处的法向量n₁, n₂, …, n_N，这些参数可通过三坐标测量仪测量或CAD模型数据计算等方式获得。为实现连续焊接，必须首先对待焊轨迹进行样条插补。本文选取B样条曲线进行插补操作。选取B样条曲线进行插补的原因是易于通过控制样条曲线的阶次使其高阶连续，以保证焊接过程各关节运动速度、加速度等参量的连续性^[12]。

1.1 待焊轨迹曲线插补

设待焊轨迹在工件坐标系{P}中的参数方程为：

$ \boldsymbol{X}\left( u \right) = \sum\limits_{i = 1}^N {{\boldsymbol{C}_i} \cdot {B_{i,q}}\left( u \right),u \in \left[ {0,1} \right].} $

(1)

式中：u为B样条曲线参数方程的自变量；C_i为B样条曲线控制点的3维坐标；B_{i, q}(u)为q-1阶B样条曲线的基函数，其节点值为

$ \underbrace {0,0, \cdots ,0}_{q + 1个0},\underbrace {{u_1},{u_2}, \cdots ,{u_{N - 1 - q}}}_{N - 1 - q个},\underbrace {1,1, \cdots ,1}_{q + 1个1}, $

(2)

且满足：

$ {B_{i,1}}\left( u \right) = \left\{ \begin{array}{l} 1,{u_i} \le u \le {u_{i + 1}},\\ 0,其他 \end{array} \right. $

(3)

$ \begin{array}{l} {B_{i,q}}\left( u \right) = \frac{{u - {u_i}}}{{{u_{i + q - 1}} - {u_i}}}{B_{i,q - 1}}\left( u \right) + \\ \frac{{{u_{i + q}} - u}}{{{u_{i + q}} - {u_{i + 1}}}}{B_{i + 1,q - 1}}\left( u \right),q \ge 2. \end{array} $

(4)

由于点X₁, X₂, …, X_N位于待焊轨迹上，因此存在β_k∈[0, 1]，使得

$ \boldsymbol{X}\left( {{\beta _k}} \right) = {\boldsymbol{X}_k}. $

(5)

为求解式(1)所示B样条曲线的未知参数以实现曲线插补，必须首先确定式(2)中的B样条曲线节点值u_j和式(5)中的自变量取值β_k。本文采用累积弦长法确定β_k取值，即

$ {\beta _k} = \left\{ \begin{array}{l} 0,k = 1\\ \frac{{\sum\limits_{i = 1}^{k - 1} {\left\| {{\boldsymbol{X}_{i + 1}} - {\boldsymbol{X}_i}} \right\|} }}{{\sum\limits_{i = 1}^{N - 1} {\left\| {{\boldsymbol{X}_{i + 1}} - {\boldsymbol{X}_i}} \right\|} }} \end{array} \right.,2 \le k \le N. $

(6)

获得β_k后，B样条曲线的节点值由式(7)确定：

$ {u_j} = \frac{1}{q}\sum\limits_{i = j + 1}^{j + q} {{\beta _i},1} \le j \le N - 1 - q. $

(7)

将式(5)代入式(1)可得

$ \sum\limits_{i = 1}^N {{\boldsymbol{C}_i} \cdot {B_{i,q}}\left( {{\beta _k}} \right) = {\boldsymbol{X}_k}.} $

(8)

式(8)为含3N个未知数、3N个方程的线性方程组，求解之即可获得B样条曲线的控制点坐标C_i，完成待焊轨迹的曲线插补。

1.2 待焊轨迹法向量插补

根据1.1节中获得的插补轨迹方程，可计算轨迹上每一点处的切向量坐标

$ \boldsymbol{s}\left( u \right) = \frac{{{\rm{d}}\boldsymbol{X}\left( u \right)}}{{{\rm{d}}u}}. $

(9)

令

$ {\boldsymbol{s}_k} = \boldsymbol{s}\left( {{\beta _k}} \right) $

(10)

为点X_k处的切向量。然而，根据式(9)和(10)计算的切向量s_k不一定与法向量n_k正交。寻找合适的B样条曲线X(u)使其经过点X_k且满足切向量与n_k正交较为困难，因此本文首先采用1.1节中的方法对待焊轨迹进行曲线插补，然后利用切向量s_k对法向量n_k进行修正，获得与切向量s_k正交的修正的法向量m_k。修正后的法向量m_k与真实法向量n_k间应保持较小的角度偏差。假设修正的法向量m_k与切向量s_k和真实法向量n_k共面，且m_k与切向量s_k正交，则

$ {\boldsymbol{m}_k} = \frac{{{\boldsymbol{n}_k} - \frac{{\boldsymbol{s}_k^{\rm{T}}{\boldsymbol{n}_k}}}{{\boldsymbol{s}_k^{\rm{T}}{\boldsymbol{s}_k}}}{\boldsymbol{s}_k}}}{{\left\| {{\boldsymbol{n}_k} - \frac{{\boldsymbol{s}_k^{\rm{T}}{\boldsymbol{n}_k}}}{{\boldsymbol{s}_k^{\rm{T}}{\boldsymbol{s}_k}}}{\boldsymbol{s}_k}} \right\|}} $

(11)

令r_k为与修正法向量m_k和切向量s_k均垂直的单位向量，即

$ {\boldsymbol{r}_k} = {\boldsymbol{m}_k}{\rm{ \times }}\frac{{{\boldsymbol{s}_k}}}{{\left\| {{\boldsymbol{s}_k}} \right\|}}. $

(12)

使用与1.1节类似的方法对向量r_k进行B样条插补，获得插补方程r(u)，并令法向量插补函数为

$ \boldsymbol{m}\left( u \right) = \boldsymbol{s}\left( u \right) \times \boldsymbol{r}\left( u \right), $

(13)

则法向量插补函数m(u)在点X_k处的取值等于修正法向量m_k，且m(u)与切向量s(u)处处正交。

2 空间曲线轨迹焊接过程运动学建模

图 1为本文使用的空间曲线轨迹焊接运动机构示意图。焊炬安装在倾角调节机构上，可在焊前预先调整焊炬倾角；倾角调节机构安装于3自由度平移台上，由3自由度平移台带动焊炬进行3维平移运动；待焊工件安装在双自由度转台上，可进行2维旋转运动。

设在焊接过程中要求焊接速度恒为C、焊炬倾角恒为α、焊炬末端与焊接区的有向距离恒为h，且焊接区始终处于水平位置。在焊接前，预先调整倾角调节机构，使焊炬轴线与竖直方向成α角，并在焊接过程中驱动双自由度转台自动调整工件姿态，保证焊接区的法向始终位于竖直位置。

如图 2所示，建立世界坐标系{W}，其横轴、纵轴和竖轴方向分别与平移台的3轴运动方向平行，原点与焊炬末端重合；建立与工件固结的工件坐标系{P}，设初始时刻工件坐标系{P}和世界坐标系{W}的横轴、纵轴和竖轴方向相互平行。

设在t时刻，焊接区在工件坐标系{P}中的坐标为X(u)，焊炬轴线的方向向量为l，双自由度转台的转角为θ和γ，3自由度平移台的位移为G(定义为焊炬末端相对于初始时刻工件坐标系{P}原点的位移量)，t时刻工件坐标系{P}与初始时刻的工件坐标系的旋转转换矩阵和平移转换矩阵分别为R(θ, γ)和T(θ, γ)。

根据焊接区的法向始终位于竖直位置可得

$ \boldsymbol{R}\left( {\theta ,\gamma } \right) \cdot \frac{{\boldsymbol{m}\left( u \right)}}{{\left\| {\boldsymbol{m}\left( u \right)} \right\|}} = {\boldsymbol{e}_3}. $

(14)

式中e₃为世界坐标系{W}竖轴的单位向量。根据式(14)可获得转台转角θ和γ与u的关系。

根据坐标系转换关系，在t时刻，焊接区在世界坐标系{W}中的坐标为

$ {\boldsymbol{X}_w} = \boldsymbol{R}\left( {\theta ,\gamma } \right) \cdot \boldsymbol{X}\left( u \right) + \boldsymbol{T}\left( {\theta ,\gamma } \right). $

(15)

由于在焊接过程中焊炬仅进行平移运动，因此在世界坐标系{W}中焊炬轴线的方向向量为常向量，

$ {\boldsymbol{l}_w} = \left[ \begin{array}{l} 0\\ \sin \alpha \\ \cos \alpha \end{array} \right]. $

(16)

根据焊炬末端与焊接区保持恒定的有向距离h，可计算焊炬末端相对于初始时刻工件坐标系{P}原点的位移量，即3自由度平移台的位移量G，

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\boldsymbol{G} = {\boldsymbol{X}_w} + \boldsymbol{h} \cdot {\boldsymbol{l}_w} = \\ \boldsymbol{R}\left( {\theta ,\gamma } \right) \cdot \boldsymbol{X}\left( u \right) + \boldsymbol{T}\left( {\theta ,\gamma } \right) + \boldsymbol{h} \cdot {\boldsymbol{l}_w}. \end{array} $

(17)

根据式(16)和(17)可计算t时刻焊炬末端在世界坐标系{W}中的瞬时速度，

$ \begin{array}{l} \boldsymbol{v} = \frac{{{\rm{d}}\boldsymbol{G}}}{{{\rm{d}}t}} = \frac{{{\rm{d}}\boldsymbol{R}\left( {\theta ,\gamma } \right)}}{{{\rm{d}}t}} \cdot \boldsymbol{X}\left( u \right) + \\ \boldsymbol{R}\left( {\theta ,\gamma } \right) \cdot \boldsymbol{s}\left( u \right) \cdot \frac{{{\rm{d}}u}}{{{\rm{d}}t}} + \frac{{{\rm{d\boldsymbol{T}}}\left( {\theta ,\gamma } \right)}}{{{\rm{d}}t}}. \end{array} $

(18)

根据式(15)可计算t时刻焊接区在世界坐标系{W}中的瞬时速度

$ \boldsymbol{v'} = \frac{{{\rm{d}}\boldsymbol{R}\left( {\theta ,\gamma } \right)}}{{{\rm{d}}t}} \cdot \boldsymbol{X}\left( u \right) + \frac{{{\rm{d\boldsymbol{T}}}\left( {\theta ,\gamma } \right)}}{{{\rm{d}}t}}. $

(19)

为保持焊接速度恒定，必须使焊炬末端瞬时速度和焊接区瞬时速度在工件表面切向方向的投影之差恒定，且等于焊接速度C，即

$ \frac{{\boldsymbol{s}_w^{\rm{T}}}}{{\left\| {{\boldsymbol{s}_w}} \right\|}} \cdot v - \frac{{\boldsymbol{s}_w^{\rm{T}}}}{{\left\| {{\boldsymbol{s}_w}} \right\|}} \cdot \boldsymbol{v'} = C. $

(20)

式中，

$ {\boldsymbol{s}_w} = \boldsymbol{R}\left( {\theta ,\gamma } \right) \cdot \boldsymbol{s}\left( u \right) $

(21)

为t时刻工件表面切向量在世界坐标系{W}中的坐标。

结合式(18)-(21)可获得最终的焊速恒定方程，

$ \left\| {\boldsymbol{s}\left( u \right)} \right\| \cdot \frac{{{\rm{d}}u}}{{{\rm{d}}t}} = C, $

(22)

即

$ \frac{{{\rm{d}}u}}{{{\rm{d}}t}} = \frac{C}{{\left\| {\boldsymbol{s}\left( u \right)} \right\|}}. $

(23)

式(23)微分方程的解为

$ t = \frac{1}{C}\int_0^u {\left\| {\boldsymbol{s}\left( \xi \right)} \right\|} {\rm{d}}\xi . $

(24)

式(24)等号右边的积分式为弧长积分，可通过数值积分的方法求解。焊速恒定方程(24)的物理意义为焊缝长度等于焊接速度和焊接时间的乘积。根据式(24)可计算u随时间t的变化关系，即可计算每个时刻焊接区的位置X(u)，并结合式(14)和(17)计算转台转角θ和γ、平移台位移G随时间t的变化关系。

在焊接过程中，除了设定不同时刻运动机构每个关节的位移量外，还需要设定各关节的瞬时速度信息。式(14)对时间t求导可获得转台各关节的瞬时角速度ω_θ和ω_γ满足的关系式，

$ \begin{array}{l} \left( {\frac{{\partial \boldsymbol{R}}}{{\partial \theta }}{\omega _\theta } + \frac{{\partial \boldsymbol{R}}}{{\partial \gamma }}{\omega _\gamma }} \right) \cdot \frac{{\boldsymbol{m}\left( u \right)}}{{\left\| {\boldsymbol{m}\left( u \right)} \right\|}} = \\ - \boldsymbol{R}\left( {\theta ,\gamma } \right)\frac{{\rm{d}}}{{{\rm{d}}u}}\left[ {\frac{{\boldsymbol{m}\left( u \right)}}{{\left\| {\boldsymbol{m}\left( u \right)} \right\|}}} \right]\frac{{{\rm{d}}u}}{{{\rm{d}}t}}. \end{array} $

(25)

即

$ \begin{array}{l} \left[ {\frac{{\partial \boldsymbol{R}}}{{\partial \theta }}\frac{{\boldsymbol{m}\left( u \right)}}{{\left\| {\boldsymbol{m}\left( u \right)} \right\|}},\frac{{\partial \boldsymbol{R}}}{{\partial \gamma }}\frac{{\boldsymbol{m}\left( u \right)}}{{\left\| {\boldsymbol{m}\left( u \right)} \right\|}}} \right]\left[ \begin{array}{l} {\omega _\theta }\\ {\omega _\gamma } \end{array} \right] = \\ - \frac{C}{{\left\| {s\left( u \right)} \right\|}}\boldsymbol{R}\left( {\theta ,\gamma } \right)\frac{{\rm{d}}}{{{\rm{d}}u}}\left[ {\frac{{\boldsymbol{m}\left( u \right)}}{{\left\| {\boldsymbol{m}\left( u \right)} \right\|}}} \right]. \end{array} $

(26)

式中，

$ \begin{array}{l} \frac{{\rm{d}}}{{{\rm{d}}u}}\left[ {\frac{{\boldsymbol{m}\left( u \right)}}{{\left\| {\boldsymbol{m}\left( u \right)} \right\|}}} \right] = \frac{1}{{\left\| {\boldsymbol{m}\left( u \right)} \right\|}}\frac{{{\rm{d}}\boldsymbol{m}\left( u \right)}}{{{\rm{d}}u}} - \\ \;\;\;\;\;\;\;\;\;\;\;\frac{{\boldsymbol{m}\left( u \right){\boldsymbol{m}^{\rm{T}}}\left( u \right)}}{{{{\left\| {\boldsymbol{m}\left( u \right)} \right\|}^3}}}\frac{{{\rm{d}}\boldsymbol{m}\left( u \right)}}{{{\rm{d}}u}}. \end{array} $

(27)

通过线性二乘法求解式(26)的线性方程组即可计算出转台各关节的瞬时角速度ω_θ和ω_γ随时间t的变化关系。

平移台的瞬时速度等于焊炬末端在世界坐标系{W}中的瞬时速度v。根据式(18)和(23)可得

$ \begin{array}{l} \boldsymbol{v} = C \cdot \boldsymbol{R}\left( {\theta ,\gamma } \right) \cdot \frac{{\boldsymbol{s}\left( u \right)}}{{\left\| {\boldsymbol{s}\left( u \right)} \right\|}} + \left[ {\frac{{\partial \boldsymbol{\boldsymbol{R}}}}{{\partial \theta }} \cdot \boldsymbol{\boldsymbol{X}}\left( u \right) + \frac{{\partial \boldsymbol{\boldsymbol{T}}}}{{\partial \theta }}} \right]{\boldsymbol{\boldsymbol{\omega}} _\theta } + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {\frac{{\partial \boldsymbol{R}}}{{\partial r}} \cdot \boldsymbol{X}\left( u \right) + \frac{{\partial \boldsymbol{T}}}{{\partial \theta }}} \right]{\boldsymbol{\omega} _r}. \end{array} $

(28)

根据式(28)可计算平移台各关节的瞬时速度随时间t的变化关系。

在图 1所示的运动机构中：

$ \boldsymbol{R}\left( {\theta ,\gamma } \right) = \left[ {\begin{array}{*{20}{c}} {\cos \gamma }&0&{\sin \gamma }\\ {\sin \theta sin\gamma }&{\cos \theta }&{ - \sin \theta sin\gamma }\\ { - \cos \theta \sin \gamma }&{\sin \theta }&{\cos \theta \cos \gamma } \end{array}} \right], $

(29)

$ \boldsymbol{T}\left( {\theta ,\gamma } \right) = \left[ {\begin{array}{*{20}{c}} 0\\ { - L\sin \theta }\\ {L\sin \theta } \end{array}} \right]. $

(30)

式中L为双自由度转台两转轴间的距离。

将式(29)和(30)代入式(14)可求得：

$ \theta = {\sin ^{ - 1}}\left[ {\boldsymbol{e}_{\rm{2}}^{\rm{T}}\frac{{\boldsymbol{m}\left( u \right)}}{{\left\| {\boldsymbol{m}\left( u \right)} \right\|}}} \right], $

(31)

$ \gamma = {\tan ^{ - 1}}\frac{{ - \boldsymbol{e}_{\rm{1}}^{\rm{T}}\boldsymbol{m}\left( u \right)}}{{\boldsymbol{e}_3^{\rm{T}}\boldsymbol{m}\left( u \right)}}. $

(32)

式中：e₁, e₂, e₃分别为世界坐标系{W}横轴、纵轴和竖轴的单位向量。

至此，运动执行机构各平移和转动关节的位移量和瞬时速度随时间t的变化关系均可求得。

3 仿真计算实例分析

本文针对图 1和2所示的相贯线轨迹采用本文提出的运动规划规划方法进行仿真计算分析。该待焊轨迹由外径分别为150 mm和100 mm的圆管相贯而成，圆管厚度均为5 mm，对外相贯线轨迹进行焊接。双自由度转台两旋转轴之间的距离L=580 mm。要求在焊接过程中焊炬维持α=5°的摆角，焊接速度恒为C=6 mm/s，焊炬末端与焊接区保持h=5 mm的距离。待焊轨迹的法向量n_k定为垂直于外相贯线且与内相贯线相交的直线的方向向量。

在待焊外相贯线轨迹上选取N=100个离散点，其在工件坐标系{P}中的坐标为：

$ \begin{array}{l} {\boldsymbol{X}_k} = {\left[ {{x_k},{y_k},{x_k}} \right]^{\rm{T}}},\\ {x_k} = \left\{ \begin{array}{l} \sqrt {{{100}^2} - y_k^2} ,1 \le k \le 51\\ - \sqrt {{{100}^2} - y_k^2} ,52 \le k \le N = 100' \end{array} \right.\\ {y_k} = \left\{ \begin{array}{l} - 100 + 4\left( {k - 1} \right),1 \le k \le 51\\ 96 + 4\left( {52 - k} \right),52 \le k \le N = 100' \end{array} \right.\\ {z_k} = \sqrt {{{150}^2} - x_k^2} . \end{array} $

(33)

设过点X_k且与内相贯线相交的直线交内相贯线于点Q_k=[p_k, q_k, r_k]^T，则

$ {\boldsymbol{n}_k} = {\boldsymbol{X}_k} - {\boldsymbol{Q}_k} $

(34)

为待焊轨迹的法向量，其满足

$ \left\{ \begin{array}{l} p_k^2 + q_k^2 = {95^2}\\ p_k^2 + r_k^2 = {145^2}\\ - {p_k}/{x_k} + {q_k}/{y_k} + {r_k}/{z_k} = 1 \end{array} \right.. $

(35)

根据式(34)和(35)可计算真实法向量n_k。

获得待焊轨迹上若干个离散点的3维坐标X_k和法向量n_k后，利用1.1节中的方法对待焊轨迹进行四阶B样条曲线插补，并据此计算各离散点处的修正法向量m_k，获得法向量插补函数m(u)，最终利用本文建立的运动学模型求解每个平移和转动关节各个时刻的位移量和瞬时速度。计算结果如图 3-6所示。计算结果表明，修正法向量m_k与真实法向量n_k间的摆角偏差不大于0.08°，因此焊炬倾角偏差不超过0.08°，可满足实际焊接需求。

4 结论

本文提出了一种适用于任意复杂空间轨迹焊接过程的运动规划方法，得到了以下结论：

1)本文提出的焊缝轨迹点和法向量插补方法可对焊道轨迹进行高精度插补，仿真计算结果表明插补后的法向量角度偏差不大于0.08°。

2)本文建立的运动学模型可保证在整个焊接过程中焊接速度、焊炬倾角、焊炬末端与焊接区距离保持恒定，对保证焊缝成形质量和冶金质量具有重要意义。

3)本文提出的运动规划方法可准确计算各时刻执行机构各关节位移量和瞬时速度信息，适用于任意复杂空间轨迹电弧焊、激光焊和搅拌摩擦焊等多种焊接场合。