三角网格曲面由于拓扑结构简单和计算效率高等优点^[1]，在汽车、船舶和模具等领域得到广泛应用。尤其随着数控加工领域针对网格曲面刀位轨迹生成方法研究的深入^[2-7]，针对基于三角网格曲面表示的工件磨削、平整需求也日益增多。然而，作为分片线性曲面，三角网格曲面缺少连续的法矢及曲率等局部微分几何信息；加之网格曲面作为真实模型的逼近曲面，存在精度损失的问题。因此，为了保证磨削过程的平稳性和磨削质量的均匀性，在利用机器人对三角网格曲面表示的工件进行磨削处理时，如何综合考虑网格曲面的局部微分特性，设计出相应的磨削工具姿态光顺优化算法，从而产生光顺性良好和过渡均匀的磨削工具姿态，成为目前亟需解决的问题。

在数控加工中，为解决三角网格曲面缺少连续的法矢及曲率的问题，不少研究学者给出了相对成熟的解决方案。在计算三角网格节点处的单位法矢时，常见的做法是对该节点周围的三角面片的法矢做面积加权平均^[1,3]，或考虑周围的三角面片形状等因素采用角度正弦加权等加权方式予以改进^[4,8]。另外，为了对曲面凹凸性判别及加工工具姿态进行合理规划，需要对节点处的曲率信息进行合理估计。相较于采用Taubin方法^[9]、基于多截形法的曲率估计^[3]即Voronoi方法^[10]等策略，本文中采用对局部顶点进行精确的参数曲面拟合的方式来进行曲率估计^[4]，该方式在操作和可控性上显得更为灵活，能够更好地为基于网格曲面的机器人磨削工具姿态光顺优化算法提供必要的曲面信息的支持。

另一方面，数控加工中在对加工工具姿态进行优化时，通常采用基于启发式算法^[11]、最短路径算法^[12]、样条拟合^[13]或径向基函数插值^[14]的方式来实现加工工具姿态间的平稳过渡。同样的，在使用机器人对网格曲面的工件进行磨削、平整的过程中，也需要对磨削工具的姿态变化的光顺性进行优化，从而避免磨削工具姿态变化剧烈导致的磨削表面的不均匀及机器人关节电机负荷过大。考虑到四元数插值在姿态过渡上的优势^[15]：能产生角速度均匀的姿态过渡、计算效率高及算法逻辑清晰，本文以四元数插值作为机器人磨削工具姿态光顺优化策略，实现对三角网格曲面磨削工具的姿态规划。

具体而言，本文将首先对网格曲面的法矢、曲率等局部微分几何信息进行计算，并基于所获得的的法矢和曲率信息，指定关键路径点的磨削工具姿态，最后通过四元数插值的方式对磨削工具姿态进行光顺，从而实现对磨削过程中机器人末端运动的优化。

1.   三角网格曲面的微分几何特性

1.1   法矢和曲率估计

曲面的法矢是加工过程中工具姿态调整的重要参考，然而三角网格曲面作为分片线性曲面，在节点处缺少连续的法矢和曲率信息。因此，对网格曲面的机器人磨削工具姿态光顺性优化而言，首先需要对网格曲面的节点处的法矢进行估算。

如图1所示，在三角网格下标为$ i $的节点$ {\upsilon _i} $处，所对应的法矢${{\boldsymbol{e}}_i}$可以由该节点星形邻域附近的面片的法矢做面积加权平均得到，记该节点附近的三角面片的面积和法矢分别为${{\boldsymbol{s}}}_{i,j}、{{\boldsymbol{e}}}_{i,j}，j=1,2,\cdots ,m$，具体公式如下。

图  1  网格曲面节点的星形邻域

下载: 全尺寸图片 幻灯片

除了法矢之外，三角网格曲面的曲率信息仍待估算。曲率反映了曲面或曲线在该点的弯曲程度，是姿态规划算法的重要依据。

对三角网格曲面的曲率估计，对计算效率要求较高，同时要求其估算的适应能力好，即能够依据精度的要求进行灵活调整。本文采用三次$ {\text{Be'zier}} $曲面进行对三角网格曲面的节点进行局部拟合，根据拟合曲面计算相应点的曲率。三次$ {\text{Be'zier}} $曲面的表达式

其中：$ {b_{i,j}} $为$ {\rm{Be'zier}} $曲面的控制顶点，$ u、v $为曲面的参数，$ {b}_{i,3}(u)、{B}_{j,3}(v) $为Bernstein基函数。

为完成三次$ {\rm{Be'zier}} $曲面的拟合过程，需要对$ {\upsilon _i} $邻域的节点$ {Q}_{j}({x}_{j},{y}_{j},{{\textit{z}}}_{j}) $进行参数化，然后利用最小二乘法求取拟合曲面。本文对节点进行参数化的过程如下。

（1） 计算点集$ {Q_j} $的最小二乘平面

记最小二乘平面的方程为

利用式（4）在极值条件下，计算得到相应的系数和最小二乘平面方程。

（2）将节点$ {Q_j} $投影至最小二乘平面上。

（3） 将所得投影点缩放至$ \left[ {0,1} \right] \times \left[ {0,1} \right] $的矩形域中，获得$ {Q_j} $对应的参数值$ \left[ {{u_j},{v_j}} \right] $。

然后，将节点$ {Q_j} $及对应的参数值$ \left[ {{u_j},{v_j}} \right] $代入式（2）中，获得线性方程组${{AX = B}}$，求解该方程，完成拟合过程。

接着，根据局部拟合得到的$ {\rm{Be'zier}} $曲面方程，若经过节点$ {\upsilon _i} $的磨削路径曲线为$ {{r}}(u(t),v(t)) $，依据式（5）计算节点$ {\upsilon _i} $处的曲率：

其中：K为曲面在该点处的法曲率，$ E、F、G $为曲面第一类基本量，$ L、M、N $为曲面第二类基本量。

1.2   关键路径点的姿态指定

利用等参数法，按步长$ \Delta t $，将磨削路径曲线$ {{r}}(u(t),v(t)) $分割成离散的路径点$ {{{r}}_k}(k = 1,2, \cdots ,m) $。随后利用1.1节中介绍的法矢和曲率估计方法计算每个离散路径点处的法矢${{{{\boldsymbol{e}}}}_k}(k = 1,2, \cdots, m)$和曲率$ {{{K}}_k}(k = 1,2, \cdots ,m) $。基于离散点处的曲率信息，按照以下流程指定关键路径点，并对关键路径点进行姿态指定：

（1）计算曲率序列$ {{{K}}_k}(k = 1,2, \cdots ,m) $的平均值。

（2）计算${K}_{c}=\overline{{K}}+\delta\left({K}_{\max }-\overline{{K}}\right) $作为曲率这一特征的阈值，其中${{{K}}_{{\rm{max}}}}$为曲率序列中的最大值，$ \delta \in \left[0，1\right] $，由人工指定，$ \delta $值越接近于1，表明设定的曲率阈值越大。

（3）遍历曲率序列，将曲率值大于曲率阈值$ {{{K}}_c} $的路径点区域视作关键区域$\left\{ {{r}_{{{k}_{1}}}},{{r}_{{{k}_{2}}}}, \cdots ,{{r}_{{{k}_{s}}}}, \cdots ,{{r}_{{{k}_{{{i}_{1}}}}}},{{r}_{{{k}_{{{i}_{2}}}}}}, \cdots ,{{r}_{{{k}_{{{i}_{t}}}}}} \right\}$，并将每个关键区域内曲率值局部最大的路径点视作关键路径点。

（4）将关键路径点的姿态指定为该点处的法矢。

至此，完成了关键路径点的机器人磨削工具姿态的指定，为后续机器人磨削工具姿态的光顺性优化做好了准备。

2.   四元数插值和姿态优化

本文采用四元数插值算法对机器人磨削工具姿态的光顺性进行优化，从而获得光顺性良好和过渡均匀的磨削工具姿态。下面介绍四元数插值的具体流程。

四元数是一种具有4个分量的、形如$ {\boldsymbol{q}}{\text{ = }}{q_1}i + {q_2}j + {q_3}k + {q_4} = \left[ {{q_1},{q_2},{q_3},{q_4}} \right] $的超复数^[16]，其中$ i,j,k $为虚轴单位，满足$ {i^2} = {j^2} = {k^2} = - 1 $，$ ij = - ji = k $，$ jk = - kj = i $，$ ik = - ki = - j $。

对于三维空间中的某一单位矢量$ {{\boldsymbol{e}}_n} = [ {e_{n,i}}, {e_{n,j}},{e_{n,k}}] $，在四元数空间内可以被表示为一个纯四元数$ {E}_{n}=\left[\left({e}_{n,i},{e}_{n,j},{e}_{n,k}\right)，0\right] $。若指定某一固定单位矢量$ {{\boldsymbol{e}}_{ref}} $作为参考矢量，则由参考矢量$ {{\boldsymbol{e}}_{ref}} $到单位矢量$ {{\boldsymbol{e}}_n} $的旋转变换可以被表示为

其中：$ {{\boldsymbol{E}}_{ref}} $为参考矢量${{\boldsymbol{e}}_{ref}}$所对应的纯四元数；${\boldsymbol{q}}_n^ *$是$ {{\boldsymbol{q}}_n} $的共轭四元数，$ {{\boldsymbol{q}}_n} $为由$ {{\boldsymbol{e}}_{ref}} $到$ {{\boldsymbol{e}}_n} $所经旋转变换对应的四元数，且该四元数$ {{\boldsymbol{q}}_n} $的计算公式为

下面利用四元数插值，对章节1中所指定的关键路径点序列$ {r}_{n}(n=1,2,\cdots ,N)$所对应的姿态序列（即该点的法矢）$ {e}_{n}(n=1,2,\cdots ,N) $进行插值，生成中间矢量作为过渡磨削工具姿态。具体做法如下：

（1） 选取参考矢量$ {{\boldsymbol{e}}_{ref}} $，并按照式（7）计算$ {{\boldsymbol{e}}_{ref}} $旋转至姿态序列中每一个姿态$ {{\boldsymbol{e}}}_{n}(n=1,2,\cdots ,N) $所对应的四元数$ {{\boldsymbol{q}}}_{n}(n=1,2,\cdots ,N) $。

（2）在所获得的四元数序列中，顺序选取两相邻四元数$ {{\boldsymbol{q}}}_{n},{{\boldsymbol{q}}}_{n+1}(n=1,2,\cdots ,N-1) $，按照式（9）进行四元数插值，得

式中：$ {m_n} = 1,2, \cdots ,{M_n};t \in \left[ {0,1} \right] $。$ {M_n} $为两关键路径点${{{r}}_n}$和${{{r}}_{n + 1}}$中间待插值的路径点数目，参数$ \left[ {0,1} \right] $可将区间$ \left[ {0,1} \right] $经$ {M_n} $等分得到。

（3）在获得所有路径点的姿态所对应的四元数后，需要将其转化为该点的姿态。具体做法是，在每个路径点处，按照插值得到的四元数$ {{\boldsymbol{q}}_{{m_n}}} $，利用式（7）对参考矢量$ {{\boldsymbol{e}}_{ref}} $做旋转变换，得到所有中间路径点的姿态$ {{\boldsymbol{e}}_{{m_n}}} $。

按照上述流程优化后的磨削工具姿态，如图2所示，可以看到生成的中间磨削工具姿态具有很好的均匀性。

图  2  四元数插值产生的磨削工具姿态

下载: 全尺寸图片 幻灯片

3.   实例验证

图3所示为MATLAB实现磨削姿态规划算法的可视化仿真环境，以及所使用的三角网格曲面。

图  3  算法的MATLAB界面及三角网格曲面

下载: 全尺寸图片 幻灯片

下面通过与未进行优化的磨削姿态进行对比，以验证本文算法的优越性。注意：未进行优化的磨削姿态指的是，将路径点处的法矢直接作为磨削工具姿态。基于上述算法流程展开仿真，仿真结果如图4所示。

图  4  仿真实例

下载: 全尺寸图片 幻灯片

图5为实验结果，通过对优化前后的磨削工具姿态运动学表现进行对比，以验证运动光顺性是否得到提高。其中，磨削工具姿态运动学表现为：整个运动过程中磨削工具姿态的角速度。

图  5  加工实验

下载: 全尺寸图片 幻灯片

利用向前差分，得到角速度的离散表达形式

式中：$ \omega $为角速度；$ f $为加工过程中的进给率；$ {\varphi _i} $为参数为$ {t_i} $时磨削工具姿态相对于绝对坐标系z轴的转角；$ {L_i} $为时间$ {t_i} $至$ {t_{i + 1}} $内磨削工具末端点移动距离。

图6所示为优化之前与优化之后的角速度曲线。可以看出，角速度的波动与优化之前相比明显降低。

图  6  磨削工具姿态的角速度

下载: 全尺寸图片 幻灯片

4.   结语

本文针对网格曲面磨削过程中遇到的问题，提出了相应的磨削姿态规划，所提出的具体算法具有以下特点：以局部拟合的方式对三角网格曲面进行法矢和曲率估计，并基于所获得的法矢和曲率信息，进行了关键路径点及其姿态的指定。然后，通过对磨削过程中机器人磨削工具姿态进行四元数插值，获得了更好的磨削工具姿态变化的均匀性，从而得到了更好的磨削工具姿态过渡的平稳性。经实例验证，本文所提算法优化后的磨削过程中，磨削工具姿态的角速度差起伏程度下降，波动减少，平稳性提高；同时，曲面的磨削效果更加均匀。