胡从军, 于广, 王立平
1. 清华大学 机械工程系, 北京 100084;
2. 精密超精密制造装备及控制北京市重点实验室, 北京 100084

收稿日期：2019-03-31
基金项目：博士后面上项目（2018M631455）；北京市青年基金资助项目（3194050）
作者简介：胡从军(1993-), 男, 硕士研究生
通信作者：王立平, 教授, E-mail:Lpwang@tsinghua.edu.cn

摘要：旋转刀轴中心控制（rotational tool center point，RTCP）精度是评价五轴混联机床的重要考核指标。该文以一台基于3P（4R）S主轴头的五轴混联机床为研究对象，对该机床终端的位置和角度精度进行标定，以提高混联机床的RTCP精度；提出了一种参数辨识算法，与最小二乘法、岭估计法相比，该算法缩短了辨识运算的时间。标定实验结果表明：该算法辨识得到的参数能有效提高机床的精度，验证了该算法的有效性，为五轴混联机床的铣削加工奠定了基础。
关键词：运动学标定参数辨识辨识算法旋转刀轴中心控制
Calibration of a 5-axis hybrid machine based on a 3P(4R)S spindle head
HU Congjun, YU Guang, WANG Liping
1.Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China;
2.Beijing Key Laboratory of Precision/Ultra-Precision Manufacturing Equipment and Control, Beijing 100084, China


Abstract: The rotational tool center point(RTCP) accuracy is an important evaluation index for evaluating five-axis hybrid machine tools. This study analyzed a 5-axis hybrid machine based on 3P(4R)S spindle head and calibrated the positioning and angular precision of the machine tool terminal to improve the RTCP precision of the hybrid machine tool. This paper presents an identification algorithm, which shortens the identification time compared with the least squares method and the ridge estimation method. Calibration tests show that the parameters identified by the algorithm can effectively improve the machine tool accuracy, which verifies the algorithm effectiveness and lays a foundation for milling using a 5-axis hybrid machine tool.
Key words: kinematic calibrationparameter identificationidentification algorithmrotational tool center point (RTCP)
近年来，中国航空航天事业得到飞速发展，其中高性能飞机发展尤为迅速，而高性能飞机的结构件通常具有尺寸大、曲面复杂、材料去除率高、加工易变形等特点，这要求加工机床切削量大、摆角灵活、速度高、精度高。为满足该加工要求，一类结合串并联机构的新型五轴混联机床被研发出来^[1-5]，其主轴头是一类并联机构，提供两转一移的3个自由度，再加上XY串联工作台提供的2个平动自由度，可共同实现5自由度任意曲面的大型飞机结构件加工^[6-8]。
在实际应用中，混联机床的零部件在加工过程中存在公差，使得机床装配完成后的实际几何参数与理论设计存在一定偏差，导致机床运动精度比较低。为了解决该问题，需要对五轴混联机床的精度进行研究，以提高机床的终端旋转刀轴中心控制(rotational tool center point，RTCP)精度。相关文献表明^[9-13]，影响混联机床静态精度的因素主要包括：零部件制造误差、零部件之间装配误差、关节间隙、自身重力以及外在负载引起的弹性变形、由热力等因素引起的结构件变形等。其中，有部分误差不随时间变化或者变化缓慢，称为准静态误差^[14-15]，大约占全部误差的70%，通过补偿消除准静态误差将显著提高混联机床终端的加工精度。五轴混联机床的精度标定方法主要为运动学标定^[16-25]，这种方法从理论模型出发，测量机床终端位置和姿态信息，再结合理论模型中机床各个部件的几何参数，利用辨识算法，计算机床各部件的误差，最后通过误差补偿方案进行补偿，从而提高终端RTCP精度。运动学标定方法主要分为误差建模、位姿测量、参数辨识、误差补偿4个部分。
误差建模和参数辨识是运动学标定的核心部分，参数辨识的结果将直接影响机床的精度，而参数辨识的关键是辨识算法。文[26-28]用最小二乘法对不同并联机构进行了运动学标定的研究，辨识出了几何参数误差，并通过补偿有效地提高了机床终端的精度。在理论上并联机构的终端可用6自由度描述当前位置和姿态，但是在实际测量时，存在某个自由度的数值如并联机床主轴刀具旋转的角度无法测量的问题，这会造成辨识矩阵奇异，从而无法进行参数辨识。Besnard^[29]在进行参数辨识时，用QR分解和奇异值分解的方法来解决辨识矩阵奇异的问题，但仍存在部分参数辨识不准确的现象。刘宇哲^[12]详细介绍了岭估计法和截断奇异值法，结合实验结果提出了综合岭估计法和截断奇异值法的奇异值局部修正算法，有效地提高了机床标定平面内的RTCP精度。但该修正算法中没有快速求得岭估计参数的方法，使得参数辨识过程复杂且耗时长。因此需要一种兼顾参数辨识精度和计算效率的辨识方法。
本文以一台基于3P(4R)S主轴头的五轴混联机床为研究对象，首先对混联机床进行误差建模，并设计合理的辨识算法；然后由测量的实验数据辨识出机构几何参数，并将其导入机床的运动学误差模型，最后完成运动学标定工作。
1 系统建模1.1 机构介绍本文研究的五轴混联机床由3P(4R)S主轴头和XY串联机构组合而成，混联机床整体模型如图 1所示。其中，3P(4R)S机构提供2个转动自由度(φ、θ)和Z方向上的平动自由度。由于设计需要，X和Y的平动自由度是分离的，3P(4R)S主轴头安装在立柱的滑块上，由滚珠丝杠驱动实现在竖直方向上的运动，提供Y方向上的平动自由度；工作台沿水平方向运动，提供X方向上的平动自由度，因此该混联机床具有5自由度。

图 1 五轴混联机床模型

图选项

机床主轴头部分为3P(4R)S并联机构，机构模型如图 2所示，由动平台、3条相同的支链、固定导轨的移动箱体3部分组成。3条支链在空间上夹角120°分布，每条支链包括滑台、2个相同连杆和三角块，2个连杆通过4个转动副R组成平行四边形机构分别与滑台和三角块连接，滑台和导轨通过移动副P连接，三角块和动平台通过球副S连接。在工作时，每条支链通过电机驱动滚珠丝杠，进而带动滑台在导轨上移动，3个滑台的运动分别通过3个平行四边形机构传递到动平台上，使动平台实现两转一移的3自由度运动。

图 2 3P(4R)S主轴头模型

图选项

1.2 误差建模在对3P(4R)S机构进行运动学分析时，为了简化分析过程，每条支链的平行四边形机构可以等效为一根连杆，此时机构简化为3PRS机构，3PRS机构简图如图 3所示，其中，B₁、B₂、B₃分别为3条支链的输入零点，A₁、A₂、A₃分别为3条支链连接动平台球副的球铰中心。在移动箱体底面建立全局坐标系O-XYZ：OXY平面为3条支链输入零点B₁、B₂、B₃构成的平面，坐标原点O位于底面中心，X轴方向为OB₁方向，Z轴垂直于底面指向动平台方向，Y轴方向由右手定则确定。在动平台上建立动坐标系O′-X′Y′Z′：O′X′Y′平面与A₁A₂A₃平面重合，坐标原点O′位于动平台上3个球铰中心A₁、A₂、A₃构成的三角形的中心，Z′轴方向垂直于A₁A₂A₃平面指向刀尖位置，X′轴方向为O′A₁方向，Y′轴方向由右手定则确定。根据图 3建立3P(4R)S机构的闭环矢量环方程

图 3 3PRS机构简图

图选项

$\mathit{\boldsymbol{H}} + {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{a}}_i} - {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{L}}_d} = {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{b}}_i} + \mathit{\boldsymbol{R}}_{i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{q}}_i} + {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{R}}_{Ci}}{\mathit{\boldsymbol{l}}_i}.$

(1)

其中：R_Bi和R_Ci分别代表P副和R副的旋转矩阵；L_d、a_i、b_i、l_i分别表示动平台中心到刀尖矢量、动平台半径矢量、静平台半径矢量、连杆长度矢量；H表示动平台刀尖的位置矢量；q_i表示每个P副的驱动矢量；R_i表示各支链的旋转矩阵。R_i可表示为

$\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{R}}_i} = \left[ {\begin{array}{*{20}{c}}{\cos {\alpha _i}}&{ - \sin {\alpha _i}}&0\\{\sin {\alpha _i}}&{\cos {\alpha _i}}&0\\0&0&1\end{array}} \right],\;\;\;{\alpha _i} = \frac{{2\left( {i - 1} \right){\rm{ \mathsf{ π} }}}}{3},}\\{i = 1,2,3.}\end{array}$

R_T为动平台相对于全局坐标系的旋转矩阵，可用T-T(tilt torsion)角^[30]表示为

$\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{R}}_{\rm{T}}} = \left[ {\begin{array}{*{20}{c}}{\cos \varphi }&{ - \sin \varphi }&0\\{\sin \varphi }&{\cos \varphi }&0\\0&0&1\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{\cos \theta }&0&{\sin \theta }\\0&1&0\\{ - \sin \theta }&0&{\cos \theta }\end{array}} \right] \cdot }\\{\left[ {\begin{array}{*{20}{c}}{\cos \varphi }&{\sin \varphi }&0\\{ - \sin \varphi }&{\cos \varphi }&0\\0&0&1\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{\cos \psi }&{ - \sin \psi }&0\\{\sin \psi }&{\cos \psi }&0\\0&0&1\end{array}} \right].}\end{array}$

其中ψ表示动平台绕Z轴方向的角度。
机构在运动时，每条支链上连杆l_i和对应导轨的夹角是变化的，即R_Ci是未知的，此时不能直接用式(1)进行数值计算，可先将式(1)改写为

$\mathit{\boldsymbol{H}} + {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{a}}_i} - {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{L}}_d} - {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{b}}_i} - {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{q}}_i} + {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{R}}_{Ci}}{\mathit{\boldsymbol{l}}_i}.$

(2)

然后，在式(2)两边同时点乘连杆l_i的单位矢量得到3PRS机构的运动方程

${\left\| {\mathit{\boldsymbol{H}} + {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{a}}_i} - {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{L}}_d} - {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{b}}_i} - {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{q}}_i}} \right\|_2} = {l_i}.$

(3)

另外，每条支链上的连杆l_i和对应导轨上的滑台为转动副连接，则连杆l_i的转动平面应与转动副的转动轴垂直，由此可得到约束方程

$\left( {\mathit{\boldsymbol{H}} + {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{a}}_i} - {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{L}}_d} - {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{b}}_i} - {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{q}}_i}} \right){\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{R}}_{Ci}}\mathit{\boldsymbol{y}} = 0.$

(4)

每条支链均可改写得到运动方程和约束方程，整个机构一共有6个方程，可求出3PRS机构终端动平台的6个自由度数据。误差模型可由理想运动学模型推导而来，通过对式微分，忽略二阶及其以上的高阶小量，并通过对比误差项系数去掉冗余误差，最终得到每条支链包含δb_ix、δb_iy、δa_ix、δa_iy、δa_iz、δθ_Bix、δθ_Biy、δθ_Cix、δθ_Ciz、δq_i、δl_i共11项误差，于是3P(4R)S机构3条支链有33项误差。由于现有研究比较详细地阐释了3PRS机构误差模型的推导过程，误差模型基本一致，这里不再详细推导，只简要列出需要用到的模型方程^{[12-13, 31]}

$\begin{array}{*{20}{c}}{\mathit{\boldsymbol{H}} + {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{R}}_i}\left( {{\mathit{\boldsymbol{a}}_i} + \delta {\mathit{\boldsymbol{a}}_i}} \right) - {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{L}}_d} = }\\{{\mathit{\boldsymbol{R}}_i}\left( {{\mathit{\boldsymbol{b}}_i} + \delta {\mathit{\boldsymbol{b}}_i}} \right) + {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{R}}_{\theta Bi}}\left( {{\mathit{\boldsymbol{q}}_i} + \delta {\mathit{\boldsymbol{q}}_i}} \right) + }\\{{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{R}}_{\theta Bi}}{\mathit{\boldsymbol{R}}_{Ci}}{\mathit{\boldsymbol{R}}_{\theta Ci}}\left( {{\mathit{\boldsymbol{l}}_i} + \delta {\mathit{\boldsymbol{l}}_i}} \right).}\end{array}$

(5)

其中：δa_i表示动平台的几何误差矢量，δb_i表示定平台的几何误差矢量，δq_i、δl_i分别表示各个滑块位置误差矢量和连杆的长度误差矢量，R_θBi表示导轨角度误差矩阵，R_θCi表示滑块角度误差矩阵。
式(5)可改写为含几何误差项的运动方程和约束方程如下：

$\begin{array}{*{20}{c}}{\left\| {\mathit{\boldsymbol{H}} + {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{R}}_i}\left( {{\mathit{\boldsymbol{a}}_i} + \delta {\mathit{\boldsymbol{a}}_i}} \right) - {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{L}}_d} - } \right.}\\{{{\left. {{\mathit{\boldsymbol{R}}_i}\left( {{\mathit{\boldsymbol{b}}_i} + \delta {\mathit{\boldsymbol{b}}_i}} \right) + {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{R}}_{\theta Bi}}\left( {{\mathit{\boldsymbol{q}}_i} + \delta {\mathit{\boldsymbol{q}}_i}} \right)} \right\|}_2} = {l_i} + \delta {l_i},}\end{array}$

(6)

$\begin{array}{*{20}{c}}{\left( {\mathit{\boldsymbol{H}} + {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{R}}_i}\left( {{\mathit{\boldsymbol{a}}_i} + \delta {\mathit{\boldsymbol{a}}_i}} \right) - {\mathit{\boldsymbol{R}}_{\rm{T}}}{\mathit{\boldsymbol{L}}_d} - {\mathit{\boldsymbol{R}}_i}\left( {{\mathit{\boldsymbol{b}}_i} + \delta {\mathit{\boldsymbol{b}}_i}} \right) - } \right.}\\{\left. {{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{R}}_{\theta Bi}}\left( {{\mathit{\boldsymbol{q}}_i} + \delta {\mathit{\boldsymbol{q}}_i}} \right)} \right) \cdot {\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{R}}_{Bi}}{\mathit{\boldsymbol{R}}_{\theta Bi}}{\mathit{\boldsymbol{R}}_{Ci}}{\mathit{\boldsymbol{R}}_{\theta Ci}}{\mathit{\boldsymbol{y}}_{Di}} = 0.}\end{array}$

(7)

2 辨识算法由于混联机床部件的实际尺寸在机床装配完成后是未知的，因此在机床投入使用之前，需要用特定的方法得到机床结构参数。一个比较简单高效的方法是通过测量机床终端输出和驱动输入，利用基于模型设计的辨识算法求解出机床部件的结构参数。在测量数据保证一定精度的前提下，求解出的结构参数的精度完全取决于辨识算法，因此辨识算法的准确性和可靠性是参数辨识的关键。
对于混联机床的主轴头3PRS机构，驱动关节变量q、终端输出p、几何参数r之间的关系表示如下：

$\mathit{\boldsymbol{F}}\left( {\mathit{\boldsymbol{p}},\mathit{\boldsymbol{q}},\mathit{\boldsymbol{r}}} \right) = 0.$

(8)

该式为函数表达，其具体表达式为式(5)，对式(8)进行全微分可得

$\frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial \mathit{\boldsymbol{q}}}}{\rm{d}}\mathit{\boldsymbol{q}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial \mathit{\boldsymbol{p}}}}{\rm{d}}\mathit{\boldsymbol{p}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial \mathit{\boldsymbol{r}}}}{\rm{d}}\mathit{\boldsymbol{r}} = 0.$

(9)

在实际测量中，驱动关节变量q可直接在机床数控系统里面读取，因而dq=0，则可得几何误差参数dr与终端输出误差dp的之间的关系：

${\rm{d}}\mathit{\boldsymbol{p}} = \mathit{\boldsymbol{J}}{\rm{d}}\mathit{\boldsymbol{r}}.$

(10)

其中J为误差映射矩阵。J可表示为

$\mathit{\boldsymbol{J}} = - {\left( {\frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial \mathit{\boldsymbol{p}}}}} \right)^{ - 1}}\frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial \mathit{\boldsymbol{r}}}}.$

对于一个确定的终端位置(x，y，z，φ，θ，ψ)，根据运动方程和约束方程可以得到6个方程式，因此，式中的误差映射矩阵J为6×33的矩阵。为了能顺利求出几何误差参数dr，至少需要6组不同的终端数据。另外，在进行参数辨识计算求解几何误差参数dr时，还需要对误差映射矩阵进行堆叠处理，得到机床并联机构的误差辨识矩阵J⁺、J⁺与J的关系为

${\mathit{\boldsymbol{J}}^ + } = {\left[ {\begin{array}{*{20}{c}}{{\mathit{\boldsymbol{J}}_1}}\\ \vdots \\{{\mathit{\boldsymbol{J}}_n}}\end{array}} \right]_{mn \times 33}}.$

(11)

其中：n表示有n组终端数据，m表示测量数据的自由度数。为了能进行辨识计算，要求mn≥33；而为了减小测量误差的干扰，应使mn?33。理论上机床终端每一个位姿有6个自由度(x，y，z，φ，θ，ψ)，即m=6。但是在实验中ψ很难测量，而且ψ是机床动平台绕主轴方向旋转的角度，该方向同时也是机床切削过程中刀具旋转的方向，其精度不需要保证。因此，没必要花大代价来测量ψ，故机床终端只有5个自由度(x，y，z，φ，θ，ψ)的数据，此时误差映射矩阵J变为5×33的矩阵。
在测量若干组数据后，需要对误差映射矩阵进行堆叠处理得到误差辨识矩阵J⁺，此时几何误差参数dr与终端输出误差dp的之间的关系为

${\rm{d}}\mathit{\boldsymbol{p}} = {\mathit{\boldsymbol{J}}^ + }{\rm{d}}\mathit{\boldsymbol{r}}.$

(12)

2.1 最小二乘法在式(12)中，由于J⁺是mn×33(其中mn?33)的矩阵，因此无法直接求逆。为了能顺利求解出机床几何误差参数dr，需要用最小二乘法求解，式(12)可以改写为

${\mathit{\boldsymbol{J}}^{ + {\rm{T}}}}{\rm{d}}\mathit{\boldsymbol{p}} = {\mathit{\boldsymbol{J}}^{ + {\rm{T}}}}{\mathit{\boldsymbol{J}}^ + }{\rm{d}}\mathit{\boldsymbol{r}}.$

(13)

此时，如果J^+TJ⁺可逆，则可直接求出机床几何误差参数，即

${\rm{d}}\mathit{\boldsymbol{r}} = {\left( {{\mathit{\boldsymbol{J}}^{ + {\rm{T}}}}{\mathit{\boldsymbol{J}}^ + }} \right)^{ - 1}}{\mathit{\boldsymbol{J}}^{ + {\rm{T}}}}{\rm{d}}\mathit{\boldsymbol{p}}.$

(14)

式(14)中，要求矩阵J^+TJ⁺必须可逆，而J^+T·J⁺可逆的条件是测量时能完整地获取终端输出p的6个自由度数据，以保证误差映射矩阵J为6×33的矩阵。
在可完整获取终端输出信息的条件下，最小二乘法是一种简单且高效的辨识算法。
2.2 岭估计法由于在实际测量时只能获取机床终端的5个自由度(x，y，z，φ，θ，ψ)的数据，导致机床终端自由度信息获取不全，误差映射矩阵J为5×33的矩阵，辨识矩阵J^+TJ⁺不可逆，最小二乘法不再适用，此时可采用岭估计法。
岭估计法引入了岭估计参数λ，岭估计参数是该算法的核心。用(J^+TJ⁺+λI)代替式(14)中的J^+TJ⁺，其中I为与J^+TJ⁺同阶的单位矩阵，岭估计算法的表达式为

${\rm{d}}\mathit{\boldsymbol{r}} = {\left( {{\mathit{\boldsymbol{J}}^{ + {\rm{T}}}}{\mathit{\boldsymbol{J}}^ + } + \lambda \mathit{\boldsymbol{I}}} \right)^{ - 1}}{\mathit{\boldsymbol{J}}^{ + {\rm{T}}}}{\rm{d}}\mathit{\boldsymbol{p}}.$

(15)

岭估计法能有效解决最小二乘算法中辨识矩阵J^+TJ⁺不可逆的问题，弥补了最小二乘法的缺陷，但是由于参数λ的引入，使求解结果严重依赖参数λ的取值。对于参数λ，目前并没有合适高效的取值算法，基本上只能通过试错法选取，而且针对不同实验数据，穷举法选取的参数λ也可能不同。岭估计法虽然解决了辨识矩阵奇异的问题，但是存在效率低、辨识计算不够准确的问题。
2.3 改进的最小二乘法辨识算法综合分析最小二乘法和岭估计法的优缺点，本文提出一种改进的辨识算法。对于缺少终端自由度数据造成辨识矩阵的奇异问题，可通过一定的算法把缺失的数据补上，此时终端自由度信息是完整的，式(14)中辨识矩阵J^+TJ⁺可逆，这就解决了因为辨识矩阵奇异而无法直接应用式(14)进行辨识计算的问题。
并联机床的参数辨识过程可转化为非线性方程组的求解过程，一般用Newton法对式(14)进行迭代求解。本算法在进行迭代计算时，将终端缺少的自由度ψ先设定为0，带入测量数据，此时终端的五自由度测量数据p(x，y，z，φ，θ，ψ)可以改写为p(x，y，z，φ，θ，ψ)，即终端位姿信息是完整的。于是误差映射矩阵J可写成6×33的矩阵，此时误差辨识矩阵J^+TJ⁺可逆，可直接用式(14)求解。为了使迭代计算能顺利进行，每一步迭代都需要对自由度ψ进行处理，为了简化计算过程，可以把每一步通过理论计算出来的ψ值当成测量值，即每迭代一次，都需要将测量数据中的ψ更新为理论计算出来的ψ值，此时终端输出残差可表示为dp=(dx，dy，dz，dφ，dθ，0)。ψ的处理依据是，当理论模型的几何参数和实际机床的几何参数相同时，理论计算的终端位姿应与实际测量数据相同，由于机床终端实测位姿信息缺少自由度ψ的数据，因此每次计算时把ψ计算值当成测量值，以保证辨识计算的顺利进行。
改进的最小二乘辨识算法有如下几个优点：
1) 新算法可直接用最小二乘法进行辨识计算，且避免了辨识矩阵的奇异问题，它与原始的最小二乘算法的区别在于对终端缺失自由度ψ的处理，使计算的残差由dp=(dx，dy，dz，dφ，dθ，dψ)变为dp=(dx，dy，dz，dφ，dθ，0)。
2) 没有引入多余的参数，不存在岭估计法面临的参数选择问题。
3) 与原始的最小二乘法相比，只需测量机床终端容易获取的5个自由度数据，即可进行参数辨识计算，提高了效率。
3 辨识算法仿真效果对比3.1 仿真参数设置及流程本节设计一个仿真算例，对前面3种辨识算法进行对比。在算例中，3PRS机构各个尺寸名义参数设定为：定平台OB_i长度b_i=220 mm，动平台O′A_i长度a_i=201 mm，连杆C_iA_i长度l_i=487 mm，动平台中心到刀尖的距离为150 mm。根据误差模型设置一组误差参数如表 1所示。本次仿真算例用的计算机为笔记本，其处理器型号为i7-7700HQ，主频为2.8 GHz，内存为24 G的DDR4内存。
表 1 误差参数

误差参数	$\frac{\delta b_{i x}}{\rm{mm}}$	$\frac{\delta b_{i y}}{\rm{mm}}$	$\frac{\delta a_{i x}}{\rm{mm}}$	$\frac{\delta a_{i y}}{\rm{mm}}$	$\frac{\delta a_{i z}}{\rm{mm}}$	$\frac{\delta \theta_{B i x}}{\left(^{\circ}\right)}$	$\frac{\delta \theta_{B i y}}{\left(^{\circ}\right)}$	$\frac{\delta \theta_{C i x}}{\left(^{\circ}\right)}$	$\frac{\delta \theta_{C i z}}{\left(^{\circ}\right)}$	$\frac{\delta q_{i}}{\rm{mm}}$	$\frac{\delta l_{i}}{\rm{mm}}$
	支链1	3.000	2.000	4.000	-1.000	2.000	0.010	0.019	0.013	0.009	3.000	2.000
支链2	2.000	1.000	2.000	-3.000	1.000	0.052	-0.020	0.020	0.010	-3.000	5.000
支链3	1.000	2.000	3.000	-5.000	-1.000	0.014	-0.013	-0.052	-0.015	1.000	-3.000

表选项






仿真流程如图 4所示，具体步骤如下：

图 4 仿真流程图

图选项

步骤1??初始化参数，设置一组误差参数δr₀和迭代初值δr⁰，迭代初值各元素通常为0。
步骤2??选择n个位姿点(z₀，φ₀，θ₀)，利用逆解方程计算各个关节变量，此时便完成了机床实际情况的模拟，得到机床输入输出信息q₀和p₀。考虑到实际测量中，测量误差无法避免，输出信息加入测量误差，则终端输出数据可表示为p₀+Δc₀，其中Δc₀表示随机测量误差，位置测量误差的取值范围为±0.01 mm，角度误差为±0.01°。
步骤3??将步骤2中得到的输入信息q₀带入运动学正解方程，求解出带误差参数δr^k的终端输出p^k，根据相关辨识算法计算堆叠求出误差辨识矩阵J⁺(δr^k)。
步骤4??利用p₀+Δc₀和p^k建立残差方程，利用辨识算法求解出误差参数δr^k+1。
步骤5??重复步骤3和4，使迭代结果满足迭代终止条件为止，得到辨识出来的误差参数δr₁。
步骤6??对比辨识得到的δr₁和设置的δr₀对终端输出的影响，评价补偿效果。
3.2 最小二乘法算例测量点所在Z面分别有Z₁=700 mm、Z₂=800 mm、Z₃=900 mm，每个Z面测点的姿态分布如图 5所示，一共选取27个位姿点(z，φ，θ)对应的数据，当求解精度为8×10^-9或者迭代次数为150次时计算终止。比较辨识算法补偿后终端在Z=700 mm处的输出误差，机构在参数补偿前的终端输出误差如图 6所示，最大输出误差为7.257 8 mm。

图 5 测量点姿态分布

图选项

图 6 机构在参数补偿前终端输出误差

图选项

用最小二乘法进行辨识计算时，终端位姿数据为6自由度数据(x，y，z，φ，θ，ψ)，辨识计算速度非常快，计算迭代了7次，耗时0.2 s。把辨识出来的误差参数导入运动学模型，计算出Z=700 mm处的输出误差如图 7所示，最大误差为0.014 9 mm(由于在模拟测量数据中加入了随机测量误差，因此辨识出来的误差参数和给定值之间存在一点偏差，该偏差大小由测量误差决定)，相比于补偿前的终端精度有很大的提高。

图 7 机构在最小二乘法辨识补偿后的终端输出误差

图选项

3.3 岭估计法算例岭估计辨识计算时，仿真流程里面模拟的终端输出数据p₀+Δc₀(x，y，z，φ，θ，ψ)要去掉ψ，将剩下5自由度(x，y，z，φ，θ，ψ)作为模拟的测量数据，通过差分的方法求误差辨识矩阵。岭估计法的辨识效果非常依赖岭估计参数λ的选取，所以在进行辨识计算之前，应该先选取合适的参数λ。目前并没有有效的方法来选取该参数，只能用试错法来进行选取。在仿真时，从数量级上来着手选取，使λ分别取1，10^-1，10^-2，…，10^-9，10^-10等11组数进行辨识计算。然后将辨识出来的误差参数导入运动学模型求出机构终端输出误差，并求出最大输出误差，得到终端最大输出误差和岭估计参数λ的关系如图 8所示。

图 8 终端最大输出误差和参数λ之间的关系

图选项

从图 8中不难发现，λ=10^-3时，机构终端输出误差最小，因此最合适的参数λ为10^-3，此时终端最大输出误差为0.010 4 mm。当参数λ取10^-4，…，10^-9，10^-10时，机构终端最大输出误差为0.015 2 mm，这比通过最小二乘法补偿后的结果大一点。另外，在仿真计算过程中，λ在分别取1，10^-1，10^-2，…，10^-9，10^-10这11个数据进行辨识计算时，迭代误差都大于0.02 mm且收敛非常慢，没有一组数据的求解精度达到了迭代停止条件8×10^-9。
为了进一步研究影响参数λ选取的因素，分别去掉测量误差和更改迭代次数为50次，再进行仿真计算，得到终端最大输出误差和参数λ的关系分别如图 9和10所示，最佳选择分别为λ=10^-4和λ=10^-3，最大输出误差分别为10^-5 mm和0.01 mm。

图 9 无测量误差时的终端最大输出误差和参数λ之间的关系

图选项

图 10 迭代50次时，终端最大输出误差和参数λ之间的关系

图选项

通过这几组仿真数据的对比可以发现，参数λ的最优值受测量误差、迭代次数等因素影响，这给参数λ的选取带来了一定的困难。综合这几组数据曲线，可选择λ=10^-3，将参数λ代入进行辨识计算，迭代了550次，迭代误差变化曲线如图 11所示，其收敛速度非常慢，计算时间为10 min左右。将辨识得到误差参数导入运动学模型，最后得到机构终端输出误差如图 12所示，其最大输出误差为0.014 4 mm。

图 11 迭代误差和迭代次数的关系

图选项

图 12 机构在岭估计法辨识补偿后的终端输出误差

图选项

3.4 改进的最小二乘法辨识算法算例改进的最小二乘法是在最小二乘算法的基础之上改进提出的，与最小二乘算法的区别在于本算法可以用于5自由度数据的辨识。在仿真计算时，仿真流程里面模拟的终端输出数据p₀+Δc₀(x，y，z，φ，θ，ψ)去掉ψ，剩下5自由度(x，y，z，φ，θ)作为模拟的测量数据，然后将每一次迭代求得的p^k中的ψ^k当作测量数据，则模拟的测量数据应为(x，y，z，φ，θ，ψ^k)(对应残差方程中Δψ=0)。最后再用最小二乘算法进行辨识计算，辨识计算迭代了23次，耗时0.8 s。将辨识结果导入运动学方程得到终端输出误差曲线如图 13所示，最大输出误差为0.015 2 mm。

图 13 新的辨识算法辨识补偿后的终端输出误差

图选项

对比3种辨识算法，在进行参数辨识时，3种算法都能明显提高机构终端输出精度，补偿后机构终端最大输出误差基本一致。总结3种算法的优缺点：最小二乘法收敛快、精度高，缺点是需要全6自由度的终端数据；岭估计法可以计算5自由度数据，由岭估计法辨识的误差参数，补偿精度相比最小二乘算法稍微高一些，但是岭估计算法需要通过试错的方法选取参数λ，而且存在收敛速度慢、参数λ受测量数据影响等问题，效率很低；改进的最小二乘法辨识算法综合了最小二乘法和岭估计法的优点，可以用5自由度数据计算，且收敛快，补偿后精度和最小二乘法算法相同。综合考虑，可以将本文提出的辨识算法作为五轴混联机床标定过程中的参数辨识算法。
4 机床标定实验对五轴混联机床进行标定实验，在工作空间中选择3个Z面进行位姿测量，每个Z面测量点的姿态分布如图 5所示，一共需要得到27个位姿点(z，φ，θ)对应的机床输入输出数据。实验分为角度测量和位置测量两部分进行，如图 14和15所示。

图 14 机床终端角度测量实验

图选项

图 15 机床终端位置测量实验

图选项

在完成数据测量之后，利用改进的最小二乘法辨识算法进行参数辨识计算，辨识出混联机床的几何误差参数，将这些参数导入机构模型，得到标定补偿后的并联机构模型。之后通过RTCP检测来验证补偿效果，在混联机床处于准静态或者低速运动的状态下，检测机床8个位姿的RTCP精度，这8个位姿的摆角均为θ=40°，补偿前后2组数据中，并联机构3个驱动轴数据相同，对比终端输出，得到补偿前后的精度分别如图 16和17所示。

图 16 标定前的RTCP精度

图选项

图 17 参数补偿后的RTCP精度

图选项

检测结果显示，混联机床在未标定时的精度误差都在毫米级，RTCP精度在X、Y、Z方向的最大误差分别为6.461、6.275和7.146 mm，精度非常差。在进行参数补偿之后，检测到RTCP精度在X、Y、Z方向的最大误差分别为0.041、0.039、0.040 mm。通过补偿前后的对比可以看出，通过参数补偿能显著提高机床终端精度，这也证明了前面辨识算法和测量方案的有效性和正确性。
5 结论本文对基于3P(4R)S主轴头的五轴混联机床的精度进行标定算法研究工作，通过对机床终端难以测量的ψ进行适当处理，避免了辨识矩阵奇异的问题，克服了岭估计法、截断奇异值法等算法计算结果依赖参数选择的不足，可直接利用最小二乘法估计方法进行参数辨识，提高了参数辨识的效率和准确率。

参考文献

[1]	张曙. 航空结构件加工的新一代数控机床——解读Ecospeed领悟机床设计之道[J]. 金属加工(冷加工), 2012(3): 2-5. ZHANG S. A new generation of CNC machine tools for aeronautical structural parts processing:interpretation of ecospeed and understanding of machine tool design[J]. Machinist Metal Cutting, 2012(3): 2-5. (in Chinese)
[2]	王立平, 汪劲松. 新型并联机床的研究状况及应用前景[J]. 航空制造技术, 2004(6): 26-30. WANG L P, WANG J S. Research status and application prospect of new parallel machine tools[J]. Aeronautical Manufacturing Technology, 2004(6): 26-30. DOI:10.3969/j.issn.1671-833X.2004.06.001 (in Chinese)
[3]	晋娆, 翊舟. Ecospeed F系列产品进军中国航空业[J]. 航空制造技术, 2016(11): 108-109. JIN R, YI Z. Ecospeed F series products entering military china aviation industry[J]. Aeronautical Manufacturing Technology, 2016(11): 108-109. (in Chinese)
[4]	LIU X J, WANG J S, PRITSCHOW G. A new family of spatial 3-DoF fully-parallel manipulators with high rotational capability[J]. Mechanism and Machine Theory, 2005, 40(4): 475-494. DOI:10.1016/j.mechmachtheory.2004.10.001
[5]	HUANG T, LI M, ZHAO X M, et al. Conceptual design and dimensional synthesis for a 3-DOF module of the trivariant-A novel 5-DOF reconfigurable hybrid robot[J]. IEEE Transactions on Robotics, 2005, 21(3): 449-456. DOI:10.1109/TRO.2004.840908
[6]	SON S, KIM T, SARMA S E, et al. A hybrid 5-axis CNC milling machine[J]. Precision Engineering, 2009, 33(4): 430-446. DOI:10.1016/j.precisioneng.2008.12.001
[7]	KANAAN D, WENGER P, CHABLAT D. Kinematic analysis of a serial-parallel machine tool:The verne machine[J]. Mechanism and Machine Theory, 2009, 44(2): 487-498. DOI:10.1016/j.mechmachtheory.2008.03.002
[8]	张京雷. SMC35混联机床的运动学标定及其专用数控系统的研究[D].北京: 清华大学, 2017. ZHANG J L. Research on kinematic calibration of the SMC35 hybrid machine tool and its numerical control system[D]. Beijing: Tsinghua University, 2017. (in Chinese)
[9]	BI Q Z, HUANG N D, SUN C, et al. Identification and compensation of geometric errors of rotary axes on five-axis machine by on-machine measurement[J]. International Journal of Machine Tools and Manufacture, 2015, 89: 182-191. DOI:10.1016/j.ijmachtools.2014.11.008
[10]	刘焕牢.数控机床几何误差测量及误差补偿技术的研究[D].武汉: 华中科技大学, 2005. LIU H L. Research on the geometric error measurement and error compensation of the numerical control machine tools[D]. Wuhan: Huazhong University of Science and Technology, 2005. (in Chinese)
[11]	张云峰, 李春雷. 数控机床几何误差测量及误差补偿技术的研究[J]. 科技传播, 2016, 8(5): 179, 185. ZHANG Y F, LI C L. Research on the geometric error measurement and error compensation of the numerical control machine Tools[J]. Public Communication of Science & Technology, 2016, 8(5): 179, 185. (in Chinese)
[12]	刘宇哲.一类并联式摆角头的几何误差建模及标定研究[D].北京: 清华大学, 2017. LIU Y Z. Research on the geometrical error modeling and calibration of a class of swing tool heads with parallel kinematics[D]. Beijing: Tsinghua University, 2017. (in Chinese)
[13]	陈禹臻.一种空间三自由度关节型摆角头的误差补偿研究[D].北京: 清华大学, 2015. CHEN Y Z. Research on the error compensation of a spacial manipulator with 3-DOF[D]. Beijing: Tsinghua University, 2015. (in Chinese)
[14]	RAMESH R, MANNAN M A, POO A N. Error compensation in machine tools-A review Part I:Geometric, cutting-force induced and fixture-dependent errors[J]. International Journal of Machine Tools and Manufacture, 2000, 40(9): 1235-1256. DOI:10.1016/S0890-6955(00)00009-2
[15]	何小妹, 丁洪生, 付铁, 等. 并联机床运动学标定研究综述[J]. 机床与液压, 2004(10): 9-11, 31. HE X M, DING H S, FU T, et al. Review on kinematic calibration of parallel kinematic machine[J]. Machine Tool & Hydraulics, 2004(10): 9-11, 31. DOI:10.3969/j.issn.1001-3881.2004.10.004 (in Chinese)
[16]	SKOPEC T, ?IKA Z, VALá?EK M. Calibration using adaptive model complexity for parallel and fiber-driven mechanisms[J]. Robotica, 2016, 34(6): 1416-1435. DOI:10.1017/S0263574715000429
[17]	刘宇哲, 吴军, 王立平, 等. 5轴混联机床运动学标定的测量轨迹评价及误差补偿[J]. 清华大学学报(自然科学版), 2016, 56(10): 1047-1054. LIU Y Z, WU J, WANG L P, et al. Measurement trajectory evaluation and error compensation for kinematic calibration of a 5-axis hybrid machine tool[J]. Journal of Tsinghua University (Science and Technology), 2016, 56(10): 1047-1054. (in Chinese)
[18]	粟时平.多轴数控机床精度建模与误差补偿方法研究[D].长沙: 中国人民解放军国防科学技术大学, 2002. SU S P. Study on the methods of precision modeling and error compensation for multi-axis CNC machine tools[D]. Changsha: National University of Defense Technology, 2002. (in Chinese)
[19]	TIAN W J, GAO W G, CHANG W F, et al. Error modeling and sensitivity analysis of a five-axis machine tool[J]. Mathematical Problems in Engineering, 2014, 2014: 745250.
[20]	ZI B, DING H F, WU X, et al. Error modeling and sensitivity analysis of a hybrid-driven based cable parallel manipulator[J]. Precision Engineering, 2014, 38(1): 197-211. DOI:10.1016/j.precisioneng.2013.06.002
[21]	HUANG T, LI Y, TANG G B, et al. Error modeling, sensitivity analysis and assembly process of a class of 3-DOF parallel kinematic machines with parallelogram struts[J]. Science in China Series E:Technological Sciences, 2002, 45(5): 467-476. DOI:10.1360/02ye9054
[22]	LI J, YU L D, SUN J Q, et al. A kinematic model for parallel-joint coordinate measuring machine[J]. Journal of Mechanisms and Robotics, 2013, 5(4): 044501. DOI:10.1115/1.4025121
[23]	CHEBBI A H, AFFI Z, ROMDHANE L. Prediction of the pose errors produced by joints clearance for a 3-UPU parallel robot[J]. Mechanism and Machine Theory, 2009, 44(9): 1768-1783. DOI:10.1016/j.mechmachtheory.2009.03.006
[24]	PINTO J M, ARRIETA C, ANDIA M E, et al. Sensitivity analysis of geometric errors in additive manufacturing medical models[J]. Medical Engineering & Physics, 2015, 37(3): 328-334.
[25]	沈斌, 邓丽芬, 劳黎露, 等. 数控机床几何误差测量及误差补偿研究[J]. 机床与液压, 2016, 44(5): 80-83. SHEN B, DENG L F, LAO L L, et al. Research of geometric error measurement and error compensation for CNC machine tools[J]. Machine Tool & Hydraulics, 2016, 44(5): 80-83. DOI:10.3969/j.issn.1001-3881.2016.05.020 (in Chinese)
[26]	MAJARENA A C, SANTOLARIA J, SAMPER D, et al. Analysis and evaluation of objective functions in kinematic calibration of parallel mechanisms[J]. The International Journal of Advanced Manufacturing Technology, 2013, 66(5-8): 751-76. DOI:10.1007/s00170-012-4363-1
[27]	WU J F, ZHANG R, WANG R H, et al. A systematic optimization approach for the calibration of parallel kinematics machine tools by a laser tracker[J]. International Journal of Machine Tools and Manufacture, 2014, 86: 1-11. DOI:10.1016/j.ijmachtools.2014.06.003
[28]	IBARAKI S, YOKAWA T, KAKINO Y, et al. Kinematic calibration on a parallel kinematic machine tool of the stewart platform by circular tests[C]//Proceedings of the 2004 American Control Conference. Boston USA: IEEE, 2004: 1394-1399.
[29]	BESNARD S, KHAIL W. Identifiable parameters for parallel robots kinematic calibration[C]//Proceedings 2001 ICRA. IEEE International Conference on Robotics and Automation (Cat. No.01CH37164). Seoul, South Korea: IEEE, 2001: 2859-2866.
[30]	BONEV, I A. Geometric analysis of parallel mechanisms[D]. Canada: Laval University, 2002.
[31]	黄鹏.一类三自由度空间并联机器精度保证研究[D].北京: 清华大学, 2011. HUANG P. Research on the accuracy assurance of a class of 3-DOF spacial parallel manipulators[D]. Beijing: Tsinghua University, 2011. (in Chinese)