CN104200063A - 机床空间加工误差的非确定性描述及预测方法 - Google Patents

Abstract

机床空间加工误差的非确定性描述及预测方法，是属于机床精度设计领域。首先，根据多体系统理论建立起机床的误差模型，在误差模型的基础上，对误差项进行合理的削减至三个方向的“当量误差”。当量误差中也同样存在着非确定性的波动，在本项发明中，加工平面时所拥有随机波动可以根据随机过程理论进行描述和预测。波动的范围也应该被限制在一定范围内；此外对加工误差波动有较大影响的关键误差项会被甄别出来，根据得到的结论，提出一些对于机床零件进行改善的地方。改善后的试验机，可以从测量结果中清晰的看到，精度的提升及波动范围的减小。这对于精密及超精密加工有着至关重要的指导意义。

Description

机床空间加工误差的非确定性描述及预测方法

技术领域

本发明是一种关于机床空间加工误差的非确定性描述及预测方法，属于机床精度设计领域。 

背景技术

随着科学技术和社会经济的快速发展，数控机床在现代加工制造和高性能装备制造的重要组成部分。如何更好的提高数控机床的加工精度成为国内外学术界的交流热点。影响机床加工精度的因素有很多，例如说几何误差、压力变形误差、热误差和动态误差等等。机床几何误差是影响加工精度的最重要部分，几乎占到所有误差的30％-40％，特别在精密及超精密的加工情况下。又因为几何误差很少受外界环境的影响，从而建立几何误差模型并进行分析对于提高加工精度有着很深的意义。 

机床的几何误差最主要来源于其导轨的制造精度还有安装精度及本身的直线度等误差。由于在安装过程中几何误差存在着一定的随机性，所以在不同的位置也会存在着一定的波动。几何误差都可以被分为两部分，一个确定量部分和符合某种概率特征分布的围绕着确定值波动的随机量部分。定义确定部分可以被补偿，随机部分也可以被控制在更小的范围内，这对加工精度的提高有着至关重要的意义。为了更好的提高数控机床的精度，误差模型的建立也是十分重要的，稳健精确的误差模型也是误差纠正和补偿的第一步。然而，随机波动的误差过大，在生产一些需要控制在特定范围内的零件，就会有面临超差报废的浪费现象。因此，如何更好的表达和分析几何误差的随机部分对于提高加工精度也是十分重要的。误差模型的建立好比绘制了一幅关于几何误差源的地图，这也是进行精度设计还有补偿的最初始最重要的部分。数十年前，科研人员主要解决的问题是机床的误差模型；在近十年，大多数的研究内容主要针对于几何误差模型建立的方法。以更稳健、简洁和精准的描述误差模型，是成为实现误差补偿最基础的要求。这些工作也大多在三轴机床上开展的。国内外专家学者一直在建立数控机床空间误差模型领域进行不懈的探索和研究，开展了多方面的工作。例如三角关系建模法、误差矩阵法、二次关系模型法、机构学建模法、刚体运动学法等。多体系统运动特征分析方法采用齐次列阵表示点的位置和矢量的姿态，在多体系统中建立广义坐标系，将机床抽象为多体系统，将在理想条件下和实际条件下的静态和动态过程中的体间的相对位置和姿态变化以及误差情况作了统一的、完整的描述，使多体系统误差的分析变得简单、迅速、明了和普遍适用，从而为实现计算机快速建模提供基础。然而，在加工过程中，由于机床加工环境具有着大量的不确定性因素，加之机床本身，零件材料等，因此一个零件拥有十分精确的尺寸是不现实的。所以，对于大多数生产过程而言，随机部分误差都被定义成可重复性的即在总体样本中取众数，而随机部分在围绕重复测量的均数以2-3倍的方差来控制随机部分的偏差。 

本项发明以经典的三轴高精度数控机床为例，对于机床的非确定性波动开展预测分析。首先，根据多体系统理论建立起机床的误差模型，在误差模型的基础上，对误差项进行合理的削减至三个方向的“当量误差”。当量误差中也同样存在着非确定性的波动，加工平面时的随机波动在本项发明中，可以根据随机过程理论进行描述和预测。波动的范围也应该被限制在一定范围内；此外对加工误差波动有较大影响的关键误差项会被甄别出来，根据得到的结论，提出一些对于机床零件进行改善的地方。改善后的试验机，可以从测量结果中清晰的看到，精度的提升及波动范围的减小。这对于精密及超精密加工有着至关重要的指导意义。 

发明内容

本发明的目的在于提供一种机床空间加工误差的非确定性描述及预测方法，现有研究技术方法中，机床的空间误差被分为两部分：确定性误差和围绕确定性误差的符合某种概率分 布的随机波动。虽然随机波动时不能被补偿的，但是尽可能减少波动范围这对于精密和超精密加工过程中，也是至关重要的。 

为实现上述目的，本发明采用的技术方案为机床空间加工误差的非确定性描述及预测方法，首先，根据多体系统理论建立起机床的误差模型，在误差模型的基础上，对误差项进行合理的削减至三个方向的“当量误差”。当量误差中也同样存在着非确定性的波动，在本发明中，加工平面时所拥有随机波动可以根据随机过程理论进行描述和预测。波动的范围也应该被限制在一定范围内；此外对加工误差波动有较大影响的关键误差项会被甄别出来，根据得到的结论，提出一些对于机床零件进行改善的地方。改善后的试验机，可以从测量结果中清晰的看到，精度的提升及波动范围的减小。这对于精密及超精密加工有着至关重要的指导意义。 

如图1所示，本方法的具体实施步骤如下， 

步骤一为三轴机床设置广义坐标系，并建立机床的空间误差模型。 

基于多体系统运动学理论，采用低序体阵列描述抽象机床系统的拓扑结构，在多体系统中建立广义坐标系，用矢量及其列向量表达位置关系，用齐次变换矩阵表示多体系统间的相互关系； 

步骤1.1建立三轴机床的拓扑结构 

分析机床的结构，定义三轴机床的各个组成部件，以及刀具和工件为“典型体”，用“B_j”表示，其中j＝0，1，2...n，j表示各典型体的序号，n-1表示机床所包含典型体的个数。 

典型体的编号规则如下： 

1)选定床身为典型体“B₀” 

2)将三轴机床分为刀具分支和工件分支，共两个分支。首先对刀具分支沿远离床身的方向，按照自然增长数列，对各典型体进行编号。再对工件分支沿远离床身的方向，按照自然增长数列，对各典型体进行编号，如图2，其中m表示刀具分支中典型体的个数，n+1表示机床总共包含的典型体的个数。 

3)任选系统中典型体B_j，体B_j的R阶低序体的序号定义为： 

L^r(j)＝i     (1) 

当B_j体为B_i体的r阶高序体(或B_j体为B_i体的相邻高序体)时，会满足： 

L^r(j)＝L(L^r-1(j))     (2) 

式中L——低序体算子； 

r、j——自然数 

且补充定义： 

L⁰(j)＝j，L^r(0)＝0     (3)，(4) 

步骤1.2建立三轴机床的特征矩阵。 

该方法所研究的三轴数控机床几何误差项(如图3所示)的几何意义及其表达式如表1所示 

表1：几何误差释义表 

在床身B₀和所有部件B_j上均建立起与其固定联接的右手直角笛卡尔三维坐标系O₀-X₀Y₀Z₀和O_j-X_jY_jZ_j，这些坐标系的集合称为广义坐标系，各体坐标系称为子坐标系，每个坐标系的三个正交基按右手定则分别取名为X，Y，Z轴；各个子坐标系的相对应的坐标轴分别对应平行；坐标轴的正方向与其所对应的运动轴的正方向相同。 

将各体之间的运动和静止情况，看作坐标系之间的运动和静止情况。根据两相邻典型体之间的静止和运动情况，在理想运动特征矩阵和误差特征矩阵表中选择相应的运动特征矩阵，如表2； 

表2：理想运动特征矩阵和运动误差特征矩阵表 

其中：T_ijS表示典型体B_j相对于典型体B_i运动的理想运动特征矩阵； 

ΔT_ijS表示典型体B_j相对于典型体B_i运动的运动误差特征矩阵； 

x_s表示沿X轴平移的距离； 

y_s表示沿Y轴平移的距离； 

z_s表示沿Z轴平移的距离； 

其余参数均已在表1(几何误差释义表)中列出。 

若相邻的典型体B_i与典型体B_j之间不存在相对运动，则理想运动特征矩阵T_ijS＝I_4×4，运动误差特征矩阵ΔT_ijS＝I_4×4，I_4×4表示4×4的单位矩阵。由于本发明是有关于机床空间加工误差的非确定性描述及预测方法，故使用过程中忽略除几何误差之外的所有误差因素，因此典型体间的体间静止特征矩阵均为T_ijP＝I_4×4。 

根据相邻典型体在静止状态下的实际位置关系，确定典型体间的体间静止误差特征矩阵ΔT_ijP

步骤1.3建立机床的空间误差模型 

刀具成型点实际运动位置与理想运动位置的偏差即为机床的空间误差。 

设刀具加工点在刀具坐标系中的坐标为： 

P_T＝[x_t，y_t，z_t，0]^T     (5) 

其中x_t表示刀具加工点在刀具坐标系中X轴方向的坐标值； 

y_t表示刀具加工点在刀具坐标系中Y轴方向的坐标值； 

z_t表示刀具加工点在刀具坐标系中Z轴方向的坐标值； 

下标t表示刀具 

机床在理想状态时成型点的运动位置： 

$\begin{array}{c}{P}_{\mathrm{wideal}}\\ ={\left[\underset{j=n,L^r\left(n\right)=0}{\overset{j=1}{\mathrm{\Π}}}{T}_{L^j\left(n\right)L^{j-1}\left(n\right)P}{T}_{L^j\left(n\right)L^{j-1}\left(n\right)S}\right]}^{-1}\left[\underset{u=r,L^r\left(m\right)=0}{\overset{u=1}{\mathrm{\Π}}}{T}_{L^u\left(m\right)L^{u-1}\left(m\right)P}{T}_{L^u\left(m\right)L^{u-1}\left(m\right)S}\right]P_T\end{array}---\left(6\right)$

式中T_ijP表示典型体B_j与典型体B_i之间的体间静止特征矩阵； 

T_ijS表示典型体B_j与典型体B_i之间的理想运动特征矩阵； 

P_T表示刀具加工点在刀具坐标系中的坐标； 

P_wideal表示理想条件下成型点在工件坐标系中的坐标， 

m+1表示刀具分支中典型体的个数； 

n+1表示三轴机床所包含的典型体的总个数。 

机床在实际状态时成型点的运动位置： 

$P_W={\left[\underset{u=n,L^r\left(n\right)=0}{\overset{u=1}{\mathrm{\Π}}}{T}_{L^u\left(n\right)L^{u-1}\left(n\right)}\right]}^{-1}\left[\underset{j=m,L^r\left(m\right)=0}{\overset{j=1}{\mathrm{\Π}}}{T}_{L^j\left(m\right)L^{j-1}\left(m\right)}\right]P_T---\left(7\right)$

其中T_ij＝T_ijP·ΔT_ijP·T_ijS·ΔT_ijS

T_ijP表示典型体B_j与典型体B_i之间的体间静止特征矩阵； 

ΔT_ijS表示典型体B_j与典型体B_i之间的体间静止误差特征矩阵； 

T_ijS表示典型体B_j与典型体B_i之间的理想运动特征矩阵； 

ΔT_ijS表示典型体B_j与典型体B_i之间的运动误差特征矩阵； 

P_T表示刀具加工点在刀具坐标系中的坐标。 

则机床的空间误差模型表示为： 

E＝P_wideal-P_W     (8) 

步骤1.4误差项合理削减及当量误差方程的建立 

本发明此步将以空间误差模型为基础，进一步将机床的所有误差项进行合理削减。机床的误差均值模型可以表示为： 

F＝F(E，G，P_W，U，U_W，U_t，G_V)     (9) 

其中： 

F＝[f₁，f₂，...，f_r]^T：r个独立方程组成的向量； 

E＝[E_x，E_y，E_z，0]^T：机床的空间误差向量； 

G＝[g₁，g₂，……，g_n]^T：n个机床各零部件几何误差组成的向量； 

G_v＝[Δγ_xy，Δβ_xz，Δα_yz，1]^T：三主轴间姿态形式误差； 

P_w＝[P_wx，P_wy，P_wz，1]^T：工件上成形点在工件坐标系中的坐标向量； 

U＝[x，y，z，B]^T：机床各运动轴的位置向量； 

U_w＝[x_w，y_w，z_w，1]^T：工件位置坐标向量； 

U_t＝[x_t，y_t，z_t，1]^T：刀具位置坐标向量； 

本发明中定义P_w，U，U_w，U_t是不存在误差。因此，可以进一步写为： 

F＝F(E，G，G_V)     (10) 

其中G的表达式可书写为： 

$G=\left[\begin{array}{cccccc}{\mathrm{\Δx}}_x& {\mathrm{\Δy}}_x& {\mathrm{\Δz}}_x& 0& 0& 0\\ {\mathrm{\Δx}}_y& {\mathrm{\Δy}}_y& {\mathrm{\Δz}}_y& 0& 0& 0\\ {\mathrm{\Δx}}_z& {\mathrm{\Δy}}_z& {\mathrm{\Δz}}_z& 0& 0& 0\\ 0& 0& 0& {\mathrm{\Δ\α}}_x& {\mathrm{\Δ\β}}_x& {\mathrm{\Δ\γ}}_x\\ 0& 0& 0& {\mathrm{\Δ\α}}_y& {\mathrm{\Δ\β}}_y& \mathrm{\Δ\γ}\\ 0& 0& 0& {\mathrm{\Δ\α}}_z& {\mathrm{\Δ\β}}_z& {\mathrm{\Δ\γ}}_z\end{array}\right]---\left(11\right)$

若存在空间误差项，可采用的方法利用激光干涉仪、球杆仪和五坐标测量仪工具来得出。其中对于机床测量方法而言，最常用的方法就是激光干涉仪。优点是可以通过一个轴的测量测出该方向上的6个误差项，总类可以分为直线度误差和线性误差，如果定义有一个与该轴运动形式趋势一致的激光干涉仪测量，此时产生的某些线性误差和直线度误差会有一定的相关性，因此，本发明中定义一个相关系数ρ表示其中的关系。 

例如：在激光干涉仪测量X向的六项基本误差时，并与此同时，在Y的方向，再加之一个激光干涉仪，运动趋势与X向运动一致，此时产生的Y项的6项基本误差就会产生一定的重叠项。X轴沿Y项的线性误差Δy_x与Y轴的定位误差Δy_y从空间来看二者是存在一定的关系，定义ρ＝Cov(Δy_x，Δy_y)便为二者的相关系数。一般情况下，设ρ＝Cov(ΔI_j，ΔJ_i)为误差与误差之间的相关系数，其中任意两项定位误差的相关系数为零。同理则可定义出其他误差项之间的相关性，矩阵： 

$U_{\mathrm{\ρ}}=\left[\begin{array}{cccccc}{\mathrm{\ρ}}_{11}& {\mathrm{\ρ}}_{12}& {\mathrm{\ρ}}_{13}& {\mathrm{\ρ}}_{14}& {\mathrm{\ρ}}_{15}& {\mathrm{\ρ}}_{16}\\ {\mathrm{\ρ}}_{21}& {\mathrm{\ρ}}_{22}& {\mathrm{\ρ}}_{23}& {\mathrm{\ρ}}_{24}& {\mathrm{\ρ}}_{25}& {\mathrm{\ρ}}_{26}\\ {\mathrm{\ρ}}_{31}& {\mathrm{\ρ}}_{32}& {\mathrm{\ρ}}_{33}& {\mathrm{\ρ}}_{34}& {\mathrm{\ρ}}_{35}& {\mathrm{\ρ}}_{36}\\ {\mathrm{\ρ}}_{41}& {\mathrm{\ρ}}_{42}& {\mathrm{\ρ}}_{43}& {\mathrm{\ρ}}_{44}& {\mathrm{\ρ}}_{45}& {\mathrm{\ρ}}_{46}\\ {\mathrm{\ρ}}_{51}& {\mathrm{\ρ}}_{52}& {\mathrm{\ρ}}_{53}& {\mathrm{\ρ}}_{54}& {\mathrm{\ρ}}_{55}& {\mathrm{\ρ}}_{56}\\ {\mathrm{\ρ}}_{61}& {\mathrm{\ρ}}_{62}& {\mathrm{\ρ}}_{63}& {\mathrm{\ρ}}_{64}& {\mathrm{\ρ}}_{65}& {\mathrm{\ρ}}_{66}\end{array}\right]---\left(12\right)$

当量误差，由于机床几何误差最终体现在定位精度上，本发明中定义一种新的误差含义：即将空间误差量，投影到在各个轴线上的误差分量。 

$\begin{array}{c}{U}_G=F\×U_{\mathrm{\ρ}}=F(E,G,G_V,U_w,U_t)\×\left[\begin{array}{cccccc}{\mathrm{\ρ}}_{11}& {\mathrm{\ρ}}_{12}& {\mathrm{\ρ}}_{13}& {\mathrm{\ρ}}_{14}& {\mathrm{\ρ}}_{15}& {\mathrm{\ρ}}_{16}\\ {\mathrm{\ρ}}_{21}& {\mathrm{\ρ}}_{22}& {\mathrm{\ρ}}_{23}& {\mathrm{\ρ}}_{24}& {\mathrm{\ρ}}_{25}& {\mathrm{\ρ}}_{26}\\ {\mathrm{\ρ}}_{31}& {\mathrm{\ρ}}_{32}& {\mathrm{\ρ}}_{33}& {\mathrm{\ρ}}_{34}& {\mathrm{\ρ}}_{35}& {\mathrm{\ρ}}_{36}\\ {\mathrm{\ρ}}_{41}& {\mathrm{\ρ}}_{42}& {\mathrm{\ρ}}_{43}& {\mathrm{\ρ}}_{44}& {\mathrm{\ρ}}_{45}& {\mathrm{\ρ}}_{46}\\ {\mathrm{\ρ}}_{51}& {\mathrm{\ρ}}_{52}& {\mathrm{\ρ}}_{53}& {\mathrm{\ρ}}_{54}& {\mathrm{\ρ}}_{55}& {\mathrm{\ρ}}_{56}\\ {\mathrm{\ρ}}_{61}& {\mathrm{\ρ}}_{62}& {\mathrm{\ρ}}_{63}& {\mathrm{\ρ}}_{64}& {\mathrm{\ρ}}_{65}& {\mathrm{\ρ}}_{66}\end{array}\right]\\ ={\left[\begin{array}{cccc}\mathrm{\Δ}{X}_x& \mathrm{\Δ}{Y}_y& \mathrm{\Δ}{Z}_z& 1\end{array}\right]}^T\end{array}---\left(13\right)$

其中： 

ΔX_x：X项上的当量误差； 

ΔY_y：Y项上的当量误差； 

ΔZ_z：Z项上的当量误差； 

最后得到当量误差方程。 

步骤二：数控机床各几何误差的测量及其测量数据的整理 

激光干涉仪被频繁的用于机床误差检测上，本发明通过定点多测的方法在X，Y，Z三个方向上进行测量。分别在各轴50-600mm的行程上，以每20mm为一节点，进行测量重复9次并计算均值。只保留误差值： 

t_r＝T_r-D   (14) 

D：目标点； 

T_r：激光干涉仪测量值； 

t_r：误差值； 

使用垂直度测量仪测量机床的三项垂直度误差。 

定义各项几何误差均符合t_r～N(μ，σ²)均符合高斯分布的独立同分布。 

$f_{\left(t_r\right)}=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\exp \{-\frac{{(t_r-\mathrm{\μ})}^2}{2{\mathrm{\σ}}^2}\}---\left(15\right)$

μ：为误差均值； 

σ²：为误差的方差； 

步骤三：计算当量误差并利用随机过程对加工轴及面的随机性波动进行描述及预测 

步骤3.1计算当量误差并进行线条拟合 

本发明中，认为ΔX_x，ΔY_y，ΔZ_z被设为独立同分布的。根据实验数据的均值，可以计算出三向的当量误差。利用B-spline曲线进行数据在位置点的拟合。拟合原理如下： 

$N_{i,p}\left(u\right)=\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}\left(u\right)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}\left(u\right)---\left(17\right)$

其中： 

u：表示当量误差； 

p：表示阶数(一般用三阶)； 

步骤3.2轴向随机性描述及预测原理 

针对其中一项误差的随机过程，可以将其称之为“高斯白噪声序列”，由白噪声过程定义可知，其中任意两点过程n₁，n₂两点的相关函数与其协方差函数相同均为σ²δ(n₁，n₂)，并且在移动过程中的任意时刻，均为不相关的，并且任意时刻均为N(0，σ²)，于是在这个过程中得到任一点的概率密度函数为： 

$f(\mathrm{\Δ}{X}_{x1},\mathrm{\Δ}{X}_{x2},...,\mathrm{\Δ}{X}_{\mathrm{xn}};n_1,n_2,....,n_n)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}f\left(\mathrm{\Δ}{X}_{\mathrm{xi}}\right)=\frac{1}{{\left(2\mathrm{\π}\right)}^{n/2}{\mathrm{\σ}}^n}\exp (-\frac{1}{2{\mathrm{\σ}}^2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\Δ}{X}_{\mathrm{xi}}^2)---\left(18\right)$

其中：ΔX_xi：为某一方向上的当量误差； 

n_n：为某一方向上的位置点； 

步骤3.3在平面上的随机性描述及预测原理 

任意两个当量误差((ΔX_x，ΔY_y)，(ΔY_y，ΔZ_z)和(ΔX_x，ΔZ_z))都是独立的随机变量，并且均符合N(0，σ²)分布。定义在X-Y平面上加工一平面，根据随机过程理论。可将平面上的任意点的误差点的误差预测为： 

{XY(n)＝ΔX_xcosωn+ΔY_ysinωn，n∈(-∞，+∞)}   (19) 

ΔX_x：X向当量误差； 

ΔY_y：Y向当量误差； 

ω：相对加工平面坐标系任一点与远点的矢量角； 

E_xy(n)属于联合高斯过程，从而也可以求出： 

E_xy(t)＝EΔX_x×cosωt+EΔY_y×sinωt＝0   (20) 

在机床运作过程中，任意两个过程点n₁，n₂时，可得到他们的相关函数与其协方差函数是相等的且均为： 

$\begin{array}{c}C(n_1,n_2)=R(n_1,n_2)=E\left[(\mathrm{\Δ}{X}_x\cos \mathrm{\ω}{n}_1+\mathrm{\Δ}{Y}_y\sin \mathrm{\ω}{n}_1)(\mathrm{\Δ}{X}_x\cos \mathrm{\ω}{n}_2+\mathrm{\Δ}{Y}_y\sin \mathrm{\ω}{n}_2)\right]\\ =\mathrm{E\Δ}{X}_x^2\×\cos \mathrm{\ω}{n}_1\cos \mathrm{\ω}{n}_2+\mathrm{E\Δ}{Y}_y\×\sin \mathrm{\ω}{n}_1\sin \mathrm{\ω}{n}_2+\\ \mathrm{E\Δ}{X}_x\×\mathrm{E\Δ}{Y}_y\×\cos \mathrm{\ω}{n}_1\sin \mathrm{\ω}{n}_2+\mathrm{E\Δ}{X}_x\×\mathrm{E\Δ}{Y}_y\×\sin \mathrm{\ω}{n}_1\cos \mathrm{\ω}{n}_2\\ ={\mathrm{\σ}}^2\cos \mathrm{\ω}(n_1-n_2)\end{array}---\left(21\right)$

由于、各个点为独立同分布的，因此对于可得到相关系数： 

$\mathrm{\ρ}=\frac{C(n_1,n_2)}{\mathrm{\σ}\left(n_1\right)\mathrm{\σ}\left(n_2\right)}=\cos \mathrm{\ω}(n_1-n_2)---\left(22\right)$

而ΔX_x，ΔY_y是服从N(0，σ²；0，σ²；cosω(n₁-n₂))，其二维密度函数为： 

$f_{\mathrm{XY}}(\mathrm{\Δ}{X}_x,\mathrm{\Δ}{Y}_y,n_1,n_2)=\frac{1}{2\mathrm{\π}\left|\sin \mathrm{\ω}(n_1-n_2)\right|}\exp [-\frac{\mathrm{\Δ}{X}_x^2-2\mathrm{\Δ}{X}_x\mathrm{\Δ}{Y}_y\cos \mathrm{\ω}(n_1-n_2)+\mathrm{\Δ}{Y}_y^2}{2{\mathrm{\σ}}^2{\sin }^2\mathrm{\ω}(n_1-n_2)}]---\left(23\right)$

依据本方法，同样可以到得到在Y-Z、X-Z面的联合概率密度函数。 

步骤四：关键误差识别与修改意见 

在本发明的前面步骤中，已经提及过当量误差及波动预测的求解方法。当量误差及其波动作为空间误差项的反应结果，如何将对空间误差项影响较大的误差甄别出来，并减少波动范围就成为此步骤的重点。控制波动范围，最直观的方法是控制影响该项的方差，根据步骤1.4提出的均值误差模型则有： 

$\begin{array}{c}{\mathrm{\σ}}_F^2={\left(\frac{\∂F}{\∂E}\right)}^2{\mathrm{\σ}}_E^2+{\left(\frac{\∂F}{\∂G}\right)}^2{\mathrm{\σ}}_G^2+{\left(\frac{\∂F}{\∂P_W}\right)}^2{\mathrm{\σ}}_{P_W}^2+{\left(\frac{\∂F}{\∂U}\right)}^2{\mathrm{\σ}}_U^2+{\left(\frac{\∂F}{\∂U_W}\right)}^2{\mathrm{\σ}}_{U_W}^2\\ +{\left(\frac{\∂F}{\∂U_t}\right)}^2{\mathrm{\σ}}_{U_t}^2+{\left(\frac{\∂F}{\∂G_V}\right)}^2{\mathrm{\σ}}_{G_V}^2\end{array}---\left(24\right)$

由于本发明只针对机床的几何误差项则有： 

$\begin{array}{c}{\mathrm{\σ}}_{G+G_V}^2={\left(\frac{\∂F}{\∂\mathrm{\Δ}{x}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{x}_x}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{y}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{y}_x}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{z}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{z}_x}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{x}_y}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{x}_y}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{y}_y}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{y}_y}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{z}_y}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{z}_y}^2+\\ {\left(\frac{\∂F}{\∂\mathrm{\Δ}{x}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{x}_z}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{y}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{y}_z}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{z}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{z}_z}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\α}}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\α}}_x}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\β}}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\β}}_x}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\γ}}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\γ}}_x}^2+\\ {\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\α}}_y}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\α}}_y}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\β}}_y}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\β}}_y}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\γ}}_y}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\γ}}_y}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\α}}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\α}}_z}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\β}}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\β}}_z}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\γ}}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\γ}}_z}^2+\\ {\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\α}}_{\mathrm{yz}}}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\α}}_{\mathrm{yz}}}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\β}}_{\mathrm{xz}}}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\β}}_{\mathrm{xz}}}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\α}}_{\mathrm{yz}}}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\α}}_{\mathrm{yz}}}^2\end{array}$

                                                                       (25) 

其中偏微分是用来具体识别出具体对于加工影响较大的误差项的，可将其就某一方向上展开归一化处理： 

$m_{\mathrm{ni}}=\frac{\left|M_{\mathrm{ni}}\right|}{\mathrm{\Σ}\left|M_{\mathrm{ni}}\right|},n=x,y,z---\left(26\right)$

m_ni的总量为1。在某一方向上的m_ni表示了该项误差对于结果影响的大小。并可以根据此原理来进行可削减波动范围的关键误差项识别工作。本发明中，为了验证预测及比较随机性效果，在各轴50-600mm的行程上，以每3mm为一节点随机记录下一组数据。证明描述预测方法的准确性与实用性。 

与现有技术相比，本发明具有如下有益效果。 

本发明是以机床的非确定性波动提供一种描述预测方法，以多体系统理论为基础建立空间误差模型，为了便于分析非确定性波动的特性及描述预测方法，提出了一种当量误差的概念；用当量误差取描述了机床加工时点、线、面产生非确定性波动的特点及预测效果；随后为了甄别出对于当量误差的数值及非确定性波动影响较大的原始误差项，提出了一种归一化处理甄别方法；最后通过修改试验机并比较，可以清楚的看到本发明提出的对于机床非确定性波动的描述预测方法对于精密及超精密加工有着实质性的生产指导意义。 

附图说明

图1为本方法实施流程图。 

图2为典型体的编号规则示意图。 

图3为一般机床误差项说明示意图。 

图4为三轴精密立式实验机床示意图。 

图5为三轴机床的拓扑结构图。 

图6为X向当量误差点及拟合图。 

图7为Y向当量误差点及拟合图。 

图8为Z向当量误差点及拟合图。 

图9为X向当量误差添加白噪声序列的随机性波动描述示意图。 

图10为Y向当量误差添加白噪声序列的随机性波动描述示意图。 

图11为Z向当量误差添加白噪声序列的随机性波动描述示意图。 

图12为在X-Y面上加工时对于非确定性波动描述示意图。 

图13为在X-Z面上加工时对于非确定性波动描述示意图。 

图14为在Y-Z面上加工时对于非确定性波动描述示意图。 

图15为在X向对加工误差影响较大的误差项分布图。 

图16为在Y向对加工误差影响较大的误差项分布图。 

图17为在Z向对加工误差影响较大的误差项分布图。 

图18为修改后试验机床X向当量误差点及拟合图。 

图19为修改后试验机床X向对行程内以3mm对的随机测量点的当量误差分布图。 

图20为修改后试验机床X向当量误差添加白噪声序列的随机性波动描述示意图。 

图21为未修改试验机床X向对行程内以3mm对的随机测量点的当量误差分布图。 

具体实施方式

本发明以三轴精密立式加工中心为例，对上述的机床空间加工误差的非确定性描述及预测方法进行验证。 

步骤一：为三轴机床设置广义坐标系，并建立机床的空间误差模型。 

基于多体系统运动学理论，采用低序体阵列描述抽象机床系统的拓扑结构，在多体系统中建立广义坐标系，用矢量及其列向量表达位置关系，用齐次变换矩阵表示多体系统间的相互关系； 

步骤1.1建立三轴机床的拓扑结构 

该机床的结构如图4所示。该机床包括X轴、刀具、工件、Y轴、Z轴、床身； 

该三轴数控机床的成型系统由X轴平动单元、Y轴平动单元、Z轴平动单元组成。在数控机床成型运动中，本发明考虑机床的几何误差。本机床共有21项几何误差，包括X，Y，Z轴的六项几何误差(Δx_x Δy_x Δz_x Δα_x Δβ_x Δγ_x Δx_y Δy_y Δz_y Δα_y Δβ_y Δγ_y Δx_z Δy_z Δz_z Δα_z Δβ_z Δγ_z)和三项垂直度误差(Δγ_XY Δβ_XZ Δα_YZ)。 

根据多体理论的基本原理将该机床抽象对多体系统，该机床主要由6个典型体组成，定义三轴机床的各个组成部件，以及刀具和工件为“典型体”，用“B_j”表示，其中j＝0，1，2，3，4，5，j表示各典型体的序号，n+1表示机床所包含典型体的个数。 

根据编号规则选定床身为典型体“B₀”，将三轴机床分为刀具分支和工件分支，共两个分支。首先对刀具分支沿远离床身的方向，按照自然增长数列，对各典型体进行编号。再对工件分支沿远离床身的方向，按照自然增长数列，对各典型体进行编号。编号结果如图5所示。 

步骤1.2建立三轴机床的特征矩阵。 

在床身B₀和所有部件B_j上均建立起与其固定联接的右手直角笛卡尔三维坐标系O₀-X₀Y₀Z₀和O_j-X_jY_jZ_j，这些坐标系的集合称为广义坐标系，各体坐标系称为子坐标系，每个坐标系的三个正交基按右手定则分别取名为X，Y，Z轴；各个子坐标系的相对应的坐标轴分别对应平行；坐标轴的正方向与其所对应的运动轴的正方向相同。 

将各体之间的运动和静止情况，看作坐标系之间的运动和静止情况。根据两相邻典型体之间的静止和运动情况，在理想运动特征矩阵和运动误差特征矩阵表(表2)中选择相应的运动特征矩阵。选择结果如表4。 

表4：该三轴机床的运动特征矩阵和运动误差特征矩阵表 

由于B₅相对于B₀无相对运动，则T_50S＝I_4×4ΔT_50S＝I_4×4； 

B₄相对于B₃无相对运动，则T_34S＝I_4×4ΔT_34S＝I_4×4。 

由于本发明是一种关于机床空间加工误差的非确定性描述及预测方法，在使用过程中忽略除几何误差之外的所有误差因素。根据相邻典型体在静止状态下的位置关系，确定典型体间静止特征矩阵和静止误差特征矩阵。结果如表5。 

表5：该三轴机床的静止特征矩阵和静止误差特征矩阵表 

步骤1.3建立机床的空间误差模型 

刀具成型点实际运动位置与理想运动位置的偏差即为机床的空间误差。 

设刀具加工点在刀具坐标系中的坐标为： 

P_T＝[x_t，y_t，z_t，0]^T    (27) 

其中x_t表示刀具加工点在刀具坐标系中X轴方向的坐标值； 

y_t表示刀具加工点在刀具坐标系中Y轴方向的坐标值； 

z_t表示刀具加工点在刀具坐标系中Z轴方向的坐标值； 

下标t表示刀具 

机床在理想状态时成型点的运动位置： 

P_wideal

＝[T_05P×T_05S]^-[T_01P×T_01S×T_12P×T_12S×T_23P×T_23S×T_34P×T_34S]P_T    （28) 

式中T_ijP表示典型体B_j与典型体B_i之间的体间静止特征矩阵； 

T_ijS表示典型体B_j与典型体B_i之间的理想运动特征矩阵； 

P_T表示刀具加工点在刀具坐标系中的坐标； 

P_widedl表示理想条件下成型点在工件坐标系中的坐标， 

机床在实际状态时成型点的运动位置： 

P_W＝[T₀₅]^-1[T₀₁×T₁₂×T₂₃×T₃₄]P_T    (29) 

其中T_ij＝T_ijP·ΔT_ijP·T_ijS·ΔT_ijS

T_ijP表示典型体B_j与典型体B_i之间的体间静止特征矩阵； 

ΔT_ijS表示典型体B_j与典型体B_i之间的体间静止误差特征矩阵； 

T_ijS表示典型体B_j与典型体B_i之间的理想运动特征矩阵； 

ΔT_ijS表示典型体B_j与典型体B_i之间的运动误差特征矩阵； 

P_T表示刀具加工点在刀具坐标系中的坐标。 

则机床的空间误差模型表示为： 

E＝P_wideal-P_W    (30) 

步骤1.4误差项合理削减及当量误差方程的建立 

本发明此步将以空间误差模型为基础，进一步将机床的所有误差项进行削减。机床的误差均值模型可以表示为： 

F＝F(E，G，P_W，U，U_W，U_t，G_V)    (31) 

其中： 

F＝[f₁，f₂，…，f_r]^T：r个独立方程组成的向量； 

E＝[E_x，E_y，E_z，0]^T：机床的空间误差向量； 

G＝[g₁，g₂，……，g_n]^T：n个机床各零部件几何误差组成的向量； 

G_v＝[Δγ_xy，Δβ_xz，Δα_yz，1]^T三主轴间姿态形式误差； 

P_w＝[P_wx，P_wy，P_wz，1]^T：工件上成形点在工件坐标系中的坐标向量； 

U＝[x，y，z，B]^T：机床各运动轴的位置向量； 

U_w＝[x_w，y_w，z_w，1]^T：工件位置坐标向量； 

U_t＝[x_t，y_t，z_t，1]^T：刀具位置坐标向量； 

由于在实际加工过程中，装夹误差及刀具装夹误差是必然要存在误差项的，故本发明中定义P_w，U是没有误差的。因此，可以进一步写为： 

F＝F(E，G，G_V，U_w，U_t)    (32) 

其中G的表达式可书写为： 

$G=\left[\begin{array}{cccccc}\mathrm{\Δ}{x}_x& \mathrm{\Δ}{y}_x& \mathrm{\Δ}{z}_x& 0& 0& 0\\ \mathrm{\Δ}{x}_y& \mathrm{\Δ}{y}_y& \mathrm{\Δ}{z}_y& 0& 0& 0\\ \mathrm{\Δ}{x}_z& \mathrm{\Δ}{y}_z& \mathrm{\Δ}{z}_z& 0& 0& 0\\ 0& 0& 0& \mathrm{\Δ}{\mathrm{\α}}_x& \mathrm{\Δ}{\mathrm{\β}}_x& \mathrm{\Δ}{\mathrm{\γ}}_x\\ 0& 0& 0& \mathrm{\Δ}{\mathrm{\α}}_y& \mathrm{\Δ}{\mathrm{\β}}_y& \mathrm{\Δ\γ}\\ 0& 0& 0& \mathrm{\Δ}{\mathrm{\α}}_z& \mathrm{\Δ}{\mathrm{\β}}_z& \mathrm{\Δ}{\mathrm{\γ}}_z\end{array}\right]---\left(33\right)$

空间误差项，利用激光干涉仪、球杆仪和五坐标测量仪来得出。其中对于机床测量方法而言，最常用的方法就是激光干涉仪。优点是可以通过一个轴的测量测出该方向上的6个误差项，总类可以分为直线度误差和线性误差，如果有一个与该轴运动形式趋势一致的激光干涉仪测量，此时产生的某些线性误差和直线度误差会有一定的相关性，因此，本发明中定义一个相关系数ρ表示其中的关系。 

在激光干涉仪测量X向的六项基本误差时，并与此同时，在Y的方向，再加之一个激光 干涉仪，运动趋势与X向运动一致，此时产生的Y项的6项基本误差就会产生一定的重叠项。X轴沿Y项的线性误差Δy_x与Y轴的定位误差Δy_y从空间来看二者是存在一定的关系，定义ρ＝Cov(Δy_x，Δy_y)便为二者的相关系数。一般情况下，设ρ＝Cov(ΔI_j，ΔJ_i)为误差与误差之间的相关系数，其中任意两项定位误差的相关系数为零。同理则可定义出其他误差项之间的相关性，矩阵： 

$U_{\mathrm{\ρ}}=\left[\begin{array}{cccccc}{\mathrm{\ρ}}_{11}& {\mathrm{\ρ}}_{12}& {\mathrm{\ρ}}_{13}& {\mathrm{\ρ}}_{14}& {\mathrm{\ρ}}_{15}& {\mathrm{\ρ}}_{16}\\ {\mathrm{\ρ}}_{21}& {\mathrm{\ρ}}_{22}& {\mathrm{\ρ}}_{23}& {\mathrm{\ρ}}_{24}& {\mathrm{\ρ}}_{25}& {\mathrm{\ρ}}_{26}\\ {\mathrm{\ρ}}_{31}& {\mathrm{\ρ}}_{32}& {\mathrm{\ρ}}_{33}& {\mathrm{\ρ}}_{34}& {\mathrm{\ρ}}_{35}& {\mathrm{\ρ}}_{36}\\ {\mathrm{\ρ}}_{41}& {\mathrm{\ρ}}_{42}& {\mathrm{\ρ}}_{43}& {\mathrm{\ρ}}_{44}& {\mathrm{\ρ}}_{45}& {\mathrm{\ρ}}_{46}\\ {\mathrm{\ρ}}_{51}& {\mathrm{\ρ}}_{52}& {\mathrm{\ρ}}_{53}& {\mathrm{\ρ}}_{54}& {\mathrm{\ρ}}_{55}& {\mathrm{\ρ}}_{56}\\ {\mathrm{\ρ}}_{61}& {\mathrm{\ρ}}_{62}& {\mathrm{\ρ}}_{63}& {\mathrm{\ρ}}_{64}& {\mathrm{\ρ}}_{65}& {\mathrm{\ρ}}_{66}\end{array}\right]---\left(34\right)$

当量误差，由于机床几何误差最终体现在定位精度上，本发明中定义一种新的误差含义：即将空间误差量，投影到在各个轴线上的误差分量。 

$\begin{array}{c}{U}_G=F\×U_{\mathrm{\ρ}}=F(E,G,G_V,U_w,U_t)\×\left[\begin{array}{cccccc}{\mathrm{\ρ}}_{11}& {\mathrm{\ρ}}_{12}& {\mathrm{\ρ}}_{13}& {\mathrm{\ρ}}_{14}& {\mathrm{\ρ}}_{15}& {\mathrm{\ρ}}_{16}\\ {\mathrm{\ρ}}_{21}& {\mathrm{\ρ}}_{22}& {\mathrm{\ρ}}_{23}& {\mathrm{\ρ}}_{24}& {\mathrm{\ρ}}_{25}& {\mathrm{\ρ}}_{26}\\ {\mathrm{\ρ}}_{31}& {\mathrm{\ρ}}_{32}& {\mathrm{\ρ}}_{33}& {\mathrm{\ρ}}_{34}& {\mathrm{\ρ}}_{35}& {\mathrm{\ρ}}_{36}\\ {\mathrm{\ρ}}_{41}& {\mathrm{\ρ}}_{42}& {\mathrm{\ρ}}_{43}& {\mathrm{\ρ}}_{44}& {\mathrm{\ρ}}_{45}& {\mathrm{\ρ}}_{46}\\ {\mathrm{\ρ}}_{51}& {\mathrm{\ρ}}_{52}& {\mathrm{\ρ}}_{53}& {\mathrm{\ρ}}_{54}& {\mathrm{\ρ}}_{55}& {\mathrm{\ρ}}_{56}\\ {\mathrm{\ρ}}_{61}& {\mathrm{\ρ}}_{62}& {\mathrm{\ρ}}_{63}& {\mathrm{\ρ}}_{64}& {\mathrm{\ρ}}_{65}& {\mathrm{\ρ}}_{66}\end{array}\right]\\ ={\left[\begin{array}{cccc}\mathrm{\Δ}{X}_x& \mathrm{\Δ}{Y}_y& \mathrm{\Δ}{Z}_z& 1\end{array}\right]}^T\end{array}---\left(35\right)$

其中： 

ΔX_x：X项上的当量误差； 

ΔY_y：Y项上的当量误差； 

ΔZ_z：Z项上的当量误差； 

最后得到当量误差方程： 

ΔX_x＝Δx_z-Δx_x-Δx_y-Δx_wd+zΔβ_x-zΔβ_wd+yΔγ_wd-zΔβ_y                (36) 

ΔY_y＝z[(Δα_x+Δα_y)-(Δy_x+Δy_y)]-x(Δγ_wd+Δγ_y+Δγ_xy)-Δy_wd+zΔα_wd (37) 

ΔZ_z＝x(Δβ_wd+Δβ_y)+Δz_z+Δz_t+Δy_z+Δy_t-zΔα_z+yΔα_wd-Δz_wd          (38) 

步骤二：数控机床各几何误差的测量及其测量数据的整理 

激光干涉仪被频繁的用于机床误差检测上，本发明中，通过定点多测的方法在X，Y，Z三个方向上进行测量。分别在各轴50-600mm的行程上，以每20mm为一节点，进行测量重复9次并计算均值。只保留误差值： 

t_r＝T_r-D       (39) 

D：目标点； 

T_r：激光干涉仪测量值； 

t_r：误差值； 

使用垂直度测量仪测量机床的三项垂直度误差。 

定义各项几何误差均符合t_r～N(μ，σ²)均符合高斯分布的独立同分布。 

$f_{\left(t_r\right)}=\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\exp \{-\frac{{(t_r-\mathrm{\μ})}^2}{{2\mathrm{\σ}}^2}\}---\left(40\right)$

μ：为误差均值； 

σ²：为误差的方差； 

本发明中，为了验证预测及比较随机性效果，在各轴50-600mm的行程上，以每3mm为一节点随机记录下一组数据。表6～9为50-600mm的行程上，以每20mm为一节点，测量9次并取均值。由于篇幅限制，只列举数据的一部分 

表6 X轴几何误差测量值均值(mm) 

表7 Y轴几何误差测量值(mm) 

表8 Z轴几何误差测量值(mm) 

表9 单元间误差测量值(mm) 

步骤三：计算当量误差并利用随机过程对加工轴及面的随机性波动进行描述及预测 

步骤3.1计算当量误差并进行线条拟合 

本发明中，定义ΔX_x，ΔY_y，ΔZ_z被设为独立同分布的。根据实验数据的均值，可以计算出三向的当量误差。利用B-spline曲线进行数据在位置点的拟合。拟合原理如下： 

$N_{i,p}\left(u\right)=\frac{u-u_i}{u_{i+p}-u_i}{N}_{i,p-1}\left(u\right)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}\left(u\right)---\left(42\right)$

其中： 

u：表示当量误差； 

p：表示阶数(一般用三阶)； 

为了更为直观的观察ΔX_x，ΔY_y，ΔZ_z当量误差及其拟合效果如图6～8所示 

步骤3.2轴向随机性描述及预测原理 

针对其中一项误差的随机过程，可以将其称之为“高斯白噪声序列”，由白噪声过程定义可知，其中任意两点过程n₁，n₂两点的相关函数与其协方差函数相同均为σ²δ(n₁，n₂)，并且在移动过程中的任意时刻，均为不相关的，并且任意时刻均为N(0，σ²)于是在这个过程中得到任一点的概率密度函数为： 

$f({\mathrm{\ΔX}}_{x1},\mathrm{\Δ}{X}_{x2},...,{\mathrm{\ΔX}}_{\mathrm{xn}};n_1,n_2,....,n_n)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}f\left({\mathrm{\ΔX}}_{\mathrm{xi}}\right)=\frac{1}{{\left(2\mathrm{\π}\right)}^{n/2}{\mathrm{\σ}}^n}\exp (-\frac{1}{{2\mathrm{\σ}}^2}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ΔX}}_{\mathrm{xi}}^2)---\left(43\right)$

其中：ΔX_xi：为某一方向上的当量误差； 

n_n：为某一方向上的位置点； 

本发明就ΔX_x，ΔY_y，ΔZ_z三向当量误差添加高斯白噪声序列，并以此来描述和预测机床的几何误差非确定性波动，其波动范围在±3σ之间。如图9～11所示 

步骤3.3在平面上的随机性描述及预测原理 

任意两个当量误差((ΔX_x，ΔY_y)，(ΔY_y，ΔZ_z)和(ΔX_x，ΔZ_z))都是独立的随机变量，并且均符合N(0，σ²)分布。定义在X-Y平面上加工一平面，根据随机过程理论。可将平面上的任意点的误差点的误差预测为： 

{XY(n)＝ΔX_xcosωn+ΔY_ysinωn，n∈(-∞，+∞)}     (44) 

ΔX_x：X向当量误差； 

ΔY_y：Y向当量误差； 

ω：相对加工平面坐标系任一点与远点的矢量角； 

E_xy(n)属于联合高斯过程，从而也可以求出： 

E_xy(t)＝EΔX_x×cosωt+EΔY_y×sinωt＝0       (45) 

在机床运作过程中，任意两个过程点n₁，n₂时，可得到他们的相关函数与其协方差函数是相等的且均为： 

$\begin{array}{c}C(n_1,n_2)=R(n_1,n_2)=E\left[({\mathrm{\ΔX}}_x\cos {\mathrm{\ωn}}_1+{\mathrm{\ΔY}}_y\sin {\mathrm{\ωn}}_1)({\mathrm{\ΔX}}_x\cos {\mathrm{\ωn}}_2+{\mathrm{\ΔY}}_y\sin {\mathrm{\ωn}}_2)\right]\\ ={\mathrm{E\ΔX}}_x^2\×\cos {\mathrm{\ωn}}_1\cos {\mathrm{\ωn}}_2+{\mathrm{E\ΔY}}_y\×\sin {\mathrm{\ωn}}_1\sin {\mathrm{\ωn}}_2+\\ {\mathrm{E\ΔX}}_x\×{\mathrm{E\ΔY}}_y\×\cos {\mathrm{\ωn}}_1\sin {\mathrm{\ωn}}_2+{\mathrm{E\ΔX}}_x\×{\mathrm{E\ΔY}}_y\×\sin {\mathrm{\ωn}}_1\cos {\mathrm{\ωn}}_2\\ ={\mathrm{\σ}}^2\cos \mathrm{\ω}(n_1-n_2)\end{array}---\left(46\right)$

由于、各个点为独立同分布的，因此对于可得到相关系数： 

$\mathrm{\ρ}=\frac{C(n_1,n_2)}{\mathrm{\σ}\left(n_1\right)\mathrm{\σ}\left(n_2\right)}=\cos \mathrm{\ω}(n_1-n_2)---\left(47\right)$

而ΔX_x，ΔY_y是服从N(0，σ²；0，σ²；cosω(n₁-n₂))，其二维密度函数为： 

$f_{\mathrm{XY}}({\mathrm{\ΔX}}_x,{\mathrm{\ΔY}}_y,n_1,n_2)=\frac{1}{2\mathrm{\π}\left|\sin \mathrm{\ω}(n_1-n_2)\right|}\exp [-\frac{{\mathrm{\ΔX}}_x^2-2{\mathrm{\ΔX}}_x{\mathrm{\ΔY}}_y\cos \mathrm{\ω}(n_1-n_2)+{\mathrm{\ΔY}}_y^2}{2{\mathrm{\σ}}^2{\sin }^2\mathrm{\ω}(n_1-n_2)}]---\left(48\right)$

依据本专利的方法，同样可以到得到在Y-Z、X-Z面的联合概率密度函数，其波动范围也应在±3σ之间。其在平面的随机波动情况如图12～14所示。 

步骤四：关键误差识别与修改意见 

在本项发明的前面步骤中，已经提及过当量误差及波动预测的求解方法。当量误差及其波动作为空间误差项的反应结果，如何将对空间误差项影响较大的误差甄别出来，并减少波动范围就成为此步骤的重点。控制波动范围，最直观的方法是控制影响该项的方差，根据步骤1.4提出的均值误差模型则有： 

$\begin{array}{c}{\mathrm{\σ}}_F^2={\left(\frac{\∂F}{\∂E}\right)}^2{\mathrm{\σ}}_E^2+{\left(\frac{\∂F}{\∂G}\right)}^2{\mathrm{\σ}}_G^2+{\left(\frac{\∂E}{\∂P_W}\right)}^2{\mathrm{\σ}}_{P_W}^2+{\left(\frac{\∂F}{\∂U}\right)}^2{\mathrm{\σ}}_U^2+{\left(\frac{\∂F}{\∂U_W}\right)}^2{\mathrm{\σ}}_{U_W}^2\\ +{\left(\frac{\∂F}{\∂U_t}\right)}^2{\mathrm{\σ}}_{U_t}^2+{\left(\frac{\∂F}{\∂G_V}\right)}^2{\mathrm{\σ}}_{G_V}^2\end{array}---\left(24\right)$

由于本项发明只针对机床的几何误差项则有： 

$\begin{array}{c}{\mathrm{\σ}}_{G+G_V}^2={\left(\frac{\∂F}{\∂\mathrm{\Δ}{x}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{x}_x}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{y}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{y}_x}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{z}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{z}_x}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{x}_y}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{x}_y}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{y}_y}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{y}_y}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{z}_y}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{z}_y}^2+\\ {\left(\frac{\∂F}{\∂\mathrm{\Δ}{x}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{x}_z}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{y}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{y}_z}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{z}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{z}_z}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\α}}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\α}}_x}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\β}}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\β}}_x}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\γ}}_x}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\γ}}_x}^2+\\ {\left(\frac{\∂F}{{\∂\mathrm{\Δ\α}}_y}\right)}^2{\mathrm{\σ}}_{{\mathrm{\Δ}\∂}_y}^2+{\left(\frac{\∂F}{{\∂\mathrm{\Δ\β}}_y}\right)}^2{\mathrm{\σ}}_{{\mathrm{\Δ\β}}_y}^2+{\left(\frac{\∂F}{{\∂\mathrm{\Δ\γ}}_y}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\γ}}_y}^2+{\left(\frac{\∂F}{{\∂\mathrm{\Δ\α}}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\α}}_z}^2+{\left(\frac{\∂F}{{\∂\mathrm{\Δ\β}}_z}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\β}}_z}^2+{\left(\frac{\∂F}{{\∂\mathrm{\Δ\γ}}_z}\right)}^2{\mathrm{\σ}}_{{\mathrm{\Δ\γ}}_z}^2+\\ {\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\α}}_{\mathrm{yz}}}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\α}}_{\mathrm{yz}}}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\β}}_{\mathrm{xz}}}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\β}}_{\mathrm{xz}}}^2+{\left(\frac{\∂F}{\∂\mathrm{\Δ}{\mathrm{\α}}_{\mathrm{yz}}}\right)}^2{\mathrm{\σ}}_{\mathrm{\Δ}{\mathrm{\α}}_{\mathrm{yz}}}^2\end{array}---\left(25\right)$

其中偏微分是用来具体识别出具体对于加工影响较大的误差项的，可将其就某一方向上展开归一化处理： 

$m_{\mathrm{ni}}=\frac{\left|M_{\mathrm{ni}}\right|}{\mathrm{\Σ}\left|M_{\mathrm{ni}}\right|},n=x,y,z---\left(26\right)$

m_ni的总量为1。在某一方向上的m_ni表示了该项误差对于结果影响的大小。图15～17分别表示了，在各个方向上对误差结果即非确定性波动范围影响较大的误差项。分别是在X向上，Δx_z，Δx_x，Δx_y，Δβ_x，Δβ_y；Y方向上，Δy_x，Δy_y，Δα_x，Δα_y，Δγ_xy；Z方向上， Δz_z，Δy_z，Δα_z，Δβ_y对加工结果有着较大影响，本发明为了更为直观的看到的其影响，图15～17表示了其影响程度。并可以根据此原理来进行可削减波动范围的关键误差项识别工作。 

本发明为了更进一步的证明方法的准确性与实用性。根据关键误差源对机床的关键部位零件进行修改，其中Δx_x，Δy_y，Δz_z三项误差源于三个主轴向的移动螺母、螺栓的制造精度和积累误差；Δx_z，Δy_z两项误差源于机床导轨的垂直面的直线度误差；Δx_y，Δy_x两项误差源于机床导轨水平面的直线度误差；Δα_x，Δβ_y两项误差取决于导轨的平行度误差；Δα_y，Δβ_y两项误差取决于机床导轨的垂直面的直线度误差及导轨长度。根据以上建议对试验机进行改善，即使用更高精度的导轨替换。为了验证更换效果，对未改善的试验机通过定点多测的方法在X，Y，Z三个方向上进行测量，分别在各轴50-600mm的行程上以每3mm为一节点随机记录下一组数据并计算其当量误差(图21)。 

对于试验机进行再次测量，通过定点多测的方法在X，Y，Z三个方向上进行测量。分别在各轴50-600mm的行程上，以每20mm为一节点，进行测量重复9次并计算均值，并以X向的当量误差为例(图18所示)。可以清楚的看到相较之前图(图6所示)，当量误差值明显降低。为了验证预测及比较随机性效果，与前面在以50-600mm的行程上为例，以每3mm为一节点随机记录下一组数据并计算其当量误差(图19所示)；并以本发明步骤3.2的描述预测波动方法，也在每3mm下进行添加白噪声序列(图20)，都以X向为例，可以清楚看到实测的当量误差(图19)与描述预测方法的步骤3.2产生的波动效果(图20)是十分相似的，证明了本发明描述预测方法的实用性。 

另外通过图19与图21的比较，也可以明晰的看到改善后的试验机波动量范围在[-0.0014mm，0.0013mm]之间，而改善前的机床的波动范围在[-0.0025mm，0.0022mm]之间。波动范围明显降低，可证明本发明提出的非确定性波动的描述预测方法对非确定性误差波动范围的减少有着实际价值，这对精密及超精密加工有着很深刻的指导意义。 

Claims (2)

1.机床空间加工误差的非确定性描述及预测方法，其特征在于：首先，根据多体系统理论建立起机床的误差模型，在误差模型的基础上，对误差项进行合理的削减至三个方向的“当量误差”；当量误差中也同样存在着非确定性的波动，在本发明中，加工平面时所拥有随机波动可以根据随机过程理论进行描述和预测；波动的范围也应该被限制在一定范围内；此外对加工误差波动有较大影响的关键误差项会被甄别出来，根据得到的结论，提出一些对于机床零件进行改善的地方； 

本方法的具体实施步骤如下， 

步骤一为三轴机床设置广义坐标系，并建立机床的空间误差模型； 

基于多体系统运动学理论，采用低序体阵列描述抽象机床系统的拓扑结构，在多体系统中建立广义坐标系，用矢量及其列向量表达位置关系，用齐次变换矩阵表示多体系统间的相互关系； 

步骤1.1建立三轴机床的拓扑结构 

分析机床的结构，定义三轴机床的各个组成部件，以及刀具和工件为“典型体”，用“B_j”表示，其中j＝0,1,2...n,j表示各典型体的序号,n-1表示机床所包含典型体的个数； 

典型体的编号规则如下： 

1)选定床身为典型体“B₀” 

2)将三轴机床分为刀具分支和工件分支，共两个分支；首先对刀具分支沿远离床身的方向，按照自然增长数列，对各典型体进行编号；再对工件分支沿远离床身的方向，按照自然增长数列，对各典型体进行编号，其中m表示刀具分支中典型体的个数，n+1表示机床总共包含的典型体的个数； 

3)任选系统中典型体B_j，体B_j的R阶低序体的序号定义为： 

L^r(j)＝i   (1) 

当B_j体为B_i体的r阶高序体(或B_j体为B_i体的相邻高序体)时，会满足： 

L^r(j)＝L(L^r-1(j))   (2) 

式中L——低序体算子； 

r、j——自然数 

且补充定义： 

L⁰(j)＝j,L^r(0)＝0   (3),(4) 

步骤1.2建立三轴机床的特征矩阵； 

该方法所研究的三轴数控机床几何误差项的几何意义及其表达式如表1所示 

表1：几何误差释义表 

在床身B₀和所有部件B_j上均建立起与其固定联接的右手直角笛卡尔三维坐标系O₀-X₀Y₀Z₀和O_j-X_jY_jZ_j，这些坐标系的集合称为广义坐标系，各体坐标系称为子坐标系，每个坐标系的三个正交基按右手定则分别取名为X,Y,Z轴；各个子坐标系的相对应的坐标轴分别对应平行；坐标轴的正方向与其所对应的运动轴的正方向相同； 

将各体之间的运动和静止情况，看作坐标系之间的运动和静止情况；根据两相邻典型体之间的静止和运动情况，在理想运动特征矩阵和误差特征矩阵表中选择相应的运动特征矩阵，如表2； 

表2：理想运动特征矩阵和运动误差特征矩阵表 

其中：T_ijS表示典型体B_j相对于典型体Bi运动的理想运动特征矩阵； 

ΔT_ijS表示典型体B_j相对于典型体B_i运动的运动误差特征矩阵； 

x_s表示沿X轴平移的距离； 

y_s表示沿Y轴平移的距离； 

z_s表示沿Z轴平移的距离； 

其余参数均已在表1中列出； 

若相邻的典型体B_i与典型体B_j之间不存在相对运动，则理想运动特征矩阵T_ijS＝I_4×4，运动误差特征矩阵ΔT_ijS＝I_4×4，I_4×4表示4×4的单位矩阵；由于本发明是有关于机床空间加工误差的非确定性描述及预测方法，故使用过程中忽略除几何误差之外的所有误差因素，因此典型体间的体间静止特征矩阵均为T_ijP＝I_4×4； 

根据相邻典型体在静止状态下的实际位置关系，确定典型体间的体间静止误差特征矩阵ΔT_ijP

步骤1.3建立机床的空间误差模型 

刀具成型点实际运动位置与理想运动位置的偏差即为机床的空间误差； 

设刀具加工点在刀具坐标系中的坐标为： 

P_T＝[x_t,y_t,z_t,0]^T   (5) 

其中x_t表示刀具加工点在刀具坐标系中X轴方向的坐标值； 

y_t表示刀具加工点在刀具坐标系中Y轴方向的坐标值； 

z_t表示刀具加工点在刀具坐标系中Z轴方向的坐标值； 

下标t表示刀具 

机床在理想状态时成型点的运动位置： 

式中T_ijP表示典型体B_j与典型体B_i之间的体间静止特征矩阵； 

T_ijS表示典型体B_j与典型体B_i之间的理想运动特征矩阵； 

P_T表示刀具加工点在刀具坐标系中的坐标； 

P_wideal表示理想条件下成型点在工件坐标系中的坐标， 

m+1表示刀具分支中典型体的个数； 

n+1表示三轴机床所包含的典型体的总个数； 

机床在实际状态时成型点的运动位置： 

其中T_ij＝T_ijP·ΔT_ijP·T_ijS·ΔT_ijS

T_ijP表示典型体B_j与典型体B_i之间的体间静止特征矩阵； 

ΔT_ijS表示典型体B_j与典型体B_i之间的体间静止误差特征矩阵； 

T_ijS表示典型体B_j与典型体B_i之间的理想运动特征矩阵； 

ΔT_ijS表示典型体B_j与典型体B_i之间的运动误差特征矩阵； 

P_T表示刀具加工点在刀具坐标系中的坐标； 

则机床的空间误差模型表示为： 

E＝P_wideal-P_W   (8) 

步骤1.4误差项合理削减及当量误差方程的建立 

本发明此步将以空间误差模型为基础，进一步将机床的所有误差项进行合理削减；机床的误差均值模型可以表示为： 

F＝F(E,G,P_W,U,U_W,U_t，G_V)   (9)其中： 

F＝[f₁,f₂,...,f_r]^T：r个独立方程组成的向量； 

E＝[E_x,E_y,E_z,0]^T：机床的空间误差向量； 

G＝[g₁,g₂,……,g_n]^T:n个机床各零部件几何误差组成的向量； 

G_v＝[Δγ_xy,Δβ_xz,Δα_yz,1]^T：三主轴间姿态形式误差； 

P_w＝[P_wx,P_wy,P_wz,1]^T：工件上成形点在工件坐标系中的坐标向量； 

U＝[x,y,z,B]^T：机床各运动轴的位置向量； 

U_w＝[x_w,y_w,z_w,1]^T：工件位置坐标向量； 

U_t＝[x_t,y_t,z_t,1]^T：刀具位置坐标向量； 

本发明中定义P_w,U,U_w,U_t是不存在误差；因此，可以进一步写为： 

F＝F(E,G,G_V)   (10) 

其中G的表达式可书写为： 

若存在空间误差项，可采用的方法利用激光干涉仪、球杆仪和五坐标测量仪工具来得出；其中对于机床测量方法而言，最常用的方法就是激光干涉仪；优点是可以通过一个轴的测量测出该方向上的6个误差项，总类可以分为直线度误差和线性误差，如果定义有一个与该轴运动形式趋势一致的激光干涉仪测量，此时产生的某些线性误差和直线度误差会有一定的相关性，因此，本发明中定义一个 相关系数ρ表示其中的关系； 

在激光干涉仪测量X向的六项基本误差时，并与此同时，在Y的方向，再加之一个激光干涉仪，运动趋势与X向运动一致，此时产生的Y项的6项基本误差就会产生一定的重叠项；X轴沿Y项的线性误差Δy_x与Y轴的定位误差Δy_y从空间来看二者是存在一定的关系，定义ρ＝Cov(Δy_x,Δy_y)便为二者的相关系数；一般情况下，设ρ＝Cov(ΔI_j,ΔJ_i)为误差与误差之间的相关系数，其中任意两项定位误差的相关系数为零；同理则可定义出其他误差项之间的相关性，矩阵： 

当量误差，由于机床几何误差最终体现在定位精度上，本发明中定义一种新的误差含义：即将空间误差量，投影到在各个轴线上的误差分量； 

其中： 

ΔX_x：X项上的当量误差； 

ΔY_y：Y项上的当量误差； 

ΔZ_z：Z项上的当量误差； 

最后得到当量误差方程； 

步骤二：数控机床各几何误差的测量及其测量数据的整理 

激光干涉仪被频繁的用于机床误差检测上，本发明通过定点多测的方法在 X，Y，Z三个方向上进行测量；分别在各轴50-600mm的行程上，以每20mm为一节点，进行测量重复9次并计算均值；只保留误差值： 

t_r＝T_r-D   (14) 

D：目标点； 

T_r：激光干涉仪测量值； 

t_r：误差值； 

使用垂直度测量仪测量机床的三项垂直度误差； 

定义各项几何误差均符合t_r～N(μ,σ²)均符合高斯分布的独立同分布； 

μ:为误差均值； 

σ²：为误差的方差； 

步骤三：计算当量误差并利用随机过程对加工轴及面的随机性波动进行描述及预测 

步骤3.1计算当量误差并进行线条拟合 

本发明中，认为ΔX_x，ΔY_y，ΔZ_z被设为独立同分布的；根据实验数据的均值，可以计算出三向的当量误差；利用B-spline曲线进行数据在位置点的拟合；拟合原理如下： 

其中： 

u：表示当量误差； 

p:表示阶数； 

步骤3.2轴向随机性描述及预测原理 

针对其中一项误差的随机过程，可以将其称之为“高斯白噪声序列”，由白噪声过程定义可知，其中任意两点过程n₁,n₂两点的相关函数与其协方差函 数相同均为σ²δ(n₁,n₂)，并且在移动过程中的任意时刻，均为不相关的，并且任意时刻均为N(0,σ²)，于是在这个过程中得到任一点的概率密度函数为： 

其中：ΔX_xi：为某一方向上的当量误差； 

n_n:为某一方向上的位置点； 

步骤3.3在平面上的随机性描述及预测原理 

任意两个当量误差((ΔX_x,ΔY_y),(ΔY_y,ΔZ_z)和(ΔX_x,ΔZ_z))都是独立的随机变量，并且均符合N(0,σ²)分布；定义在X-Y平面上加工一平面，根据随机过程理论；可将平面上的任意点的误差点的误差预测为： 

{XY(n)＝ΔX_xcosωn+ΔY_ysinωn,n∈(-∞,+∞)}   (19) 

ΔX_x：X向当量误差； 

ΔY_y：Y向当量误差； 

ω：相对加工平面坐标系任一点与远点的矢量角； 

E_xy(n)属于联合高斯过程，从而也可以求出： 

E_xy(t)＝EΔX_x×cosωt+EΔY_y×sinωt＝0   (20) 

在机床运作过程中，任意两个过程点n₁,n₂时，可得到他们的相关函数与其协方差函数是相等的且均为： 

由于、各个点为独立同分布的，因此对于可得到相关系数： 

而ΔX_x，ΔY_y是服从N(0,σ²；0,σ²；cosω(n₁-n₂))，其二维密度函数为： 

依据本方法，同样可以到得到在Y-Z、X-Z面的联合概率密度函数； 

步骤四：关键误差识别与修改意见 

当量误差及其波动作为空间误差项的反应结果，如何将对空间误差项影响较大的误差甄别出来，并减少波动范围就成为此步骤的重点；控制波动范围，最直观的方法是控制影响该项的方差，根据步骤1.4提出的均值误差模型则有： 

由于本发明只针对机床的几何误差项则有： 

其中偏微分是用来具体识别出具体对于加工影响较大的误差项的，可将其就某一方向上展开归一化处理： 

m_ni的总量为1；在某一方向上的m_ni表示了该项误差对于结果影响的大小； 

并可以根据此原理来进行可削减波动范围的关键误差项识别工作； 

本发明中，为了验证预测及比较随机性效果，在各轴50-600mm的行程上，以每3mm为一节点随机记录下一组数据。 

2.根据权利要求1所述的机床空间加工误差的非确定性描述及预测方法，其特征在于：本发明以三轴精密立式加工中心为例，对上述的机床空间加工误差的非确定性描述及预测方法进行验证； 

步骤一：为三轴机床设置广义坐标系，并建立机床的空间误差模型； 

基于多体系统运动学理论，采用低序体阵列描述抽象机床系统的拓扑结构，在多体系统中建立广义坐标系，用矢量及其列向量表达位置关系，用齐次变换矩阵表示多体系统间的相互关系； 

步骤1.1建立三轴机床的拓扑结构 

该机床包括X轴、刀具、工件、Y轴、Z轴、床身； 

该三轴数控机床的成型系统由X轴平动单元、Y轴平动单元、Z轴平动单元组成；在数控机床成型运动中，本发明考虑机床的几何误差；本机床共有21项几何误差，包括X,Y,Z轴的六项几何误差(Δx_xΔy_xΔz_xΔα_xΔβ_xΔγ_xΔx_yΔy_yΔz_yΔα_yΔβ_yΔγ_yΔx_zΔy_zΔz_zΔα_zΔβ_zΔγ_z)和三项垂直度误差(Δγ_XYΔβ_XZΔα_YZ)； 

根据多体理论的基本原理将该机床抽象对多体系统，该机床主要由6个典型体组成，定义三轴机床的各个组成部件，以及刀具和工件为“典型体”，用“B_j”表示，其中j＝0,1,2,3,4,5,j表示各典型体的序号,n+1表示机床所包含典型体的个数； 

根据编号规则选定床身为典型体“B₀”，将三轴机床分为刀具分支和工件分支，共两个分支；首先对刀具分支沿远离床身的方向，按照自然增长数列，对各典型体进行编号；再对工件分支沿远离床身的方向，按照自然增长数列，对各典型体进行编号； 

步骤1.2建立三轴机床的特征矩阵； 

在床身B₀和所有部件B_j上均建立起与其固定联接的右手直角笛卡尔三维坐标系O₀-X₀Y₀Z₀和O_j-X_jY_jZ_j，这些坐标系的集合称为广义坐标系，各体坐标系称为子坐标系，每个坐标系的三个正交基按右手定则分别取名为X,Y,Z轴；各个子坐标系的相对应的坐标轴分别对应平行；坐标轴的正方向与其所对应的运动轴的正方向相同； 

将各体之间的运动和静止情况，看作坐标系之间的运动和静止情况；根据两相邻典型体之间的静止和运动情况，在理想运动特征矩阵和运动误差特征矩阵表 中选择相应的运动特征矩阵；选择结果如表4； 

表4：该三轴机床的运动特征矩阵和运动误差特征矩阵表 

由于B₅相对于B₀无相对运动，则T_50S＝I_4×4ΔT_50S＝I_4×4； 

B₄相对于B₃无相对运动，则T_34S＝I_4×4ΔT_34S＝I_4×4； 

由于本发明是一种关于机床空间加工误差的非确定性描述及预测方法，在使用过程中忽略除几何误差之外的所有误差因素；根据相邻典型体在静止状态下的位置关系，确定典型体间静止特征矩阵和静止误差特征矩阵；结果如表5； 

表5：该三轴机床的静止特征矩阵和静止误差特征矩阵表 

步骤1.3建立机床的空间误差模型 

刀具成型点实际运动位置与理想运动位置的偏差即为机床的空间误差； 

设刀具加工点在刀具坐标系中的坐标为： 

P_T＝[x_t,y_t,z_t,0]^T   (27) 

其中x_t表示刀具加工点在刀具坐标系中X轴方向的坐标值； 

y_t表示刀具加工点在刀具坐标系中Y轴方向的坐标值； 

z_t表示刀具加工点在刀具坐标系中Z轴方向的坐标值； 

下标t表示刀具 

机床在理想状态时成型点的运动位置： 

P_wideal

                                                       (28) 

＝[T_05P×T_05S]-[T_01P×T_01S×T_12P×T_12S×T_23P×T_23S×T_34P×T_34S]P_T

式中T_ijP表示典型体B_j与典型体B_i之间的体间静止特征矩阵； 

T_ijS表示典型体B_j与典型体B_i之间的理想运动特征矩阵； 

P_T表示刀具加工点在刀具坐标系中的坐标； 

P_wideal表示理想条件下成型点在工件坐标系中的坐标， 

机床在实际状态时成型点的运动位置： 

P_W＝[T₀₅]^-1[T₀₁×T₁₂×T₂₃×T₃₄]P_T   (29) 

其中T_ij＝T_ijP·ΔT_ijP·T_ijS·ΔT_ijS

T_ijP表示典型体B_j与典型体B_i之间的体间静止特征矩阵； 

ΔT_ijS表示典型体B_j与典型体B_i之间的体间静止误差特征矩阵； 

T_ijS表示典型体B_j与典型体B_i之间的理想运动特征矩阵； 

ΔT_ijS表示典型体B_j与典型体B_i之间的运动误差特征矩阵； 

P_T表示刀具加工点在刀具坐标系中的坐标； 

则机床的空间误差模型表示为： 

E＝P_wideal-P_W   (30) 

步骤1.4误差项合理削减及当量误差方程的建立 

本发明此步将以空间误差模型为基础，进一步将机床的所有误差项进行削减；机床的误差均值模型可以表示为： 

F＝F(E,G,P_W,U,U_W,U_t，G_V)   (31) 

其中： 

F＝[f₁,f₂,...,f_r]^T：r个独立方程组成的向量； 

E＝[E_x,E_y,E_z,0]^T：机床的空间误差向量； 

G＝[g₁,g₂,……,g_n]^T:n个机床各零部件几何误差组成的向量； 

G_v＝[Δγ_xy,Δβ_xz,Δα_yz,1]^T：三主轴间姿态形式误差； 

P_w＝[P_wx,P_wy,P_wz,1]^T：工件上成形点在工件坐标系中的坐标向量； 

U＝[x,y,z,B]^T：机床各运动轴的位置向量； 

U_w＝[x_w,y_w,z_w,1]^T：工件位置坐标向量； 

U_t＝[x_t,y_t,z_t,1]^T：刀具位置坐标向量； 

由于在实际加工过程中，装夹误差及刀具装夹误差是必然要存在误差项的，故本发明中定义P_w,U是没有误差的；因此，可以进一步写为： 

F＝F(E,G,G_V,U_w,U_t)   (32) 

其中G的表达式可书写为： 

空间误差项，利用激光干涉仪、球杆仪和五坐标测量仪来得出；其中对于机床测量方法而言，最常用的方法就是激光干涉仪；优点是可以通过一个轴的测量测出该方向上的6个误差项，总类可以分为直线度误差和线性误差，如果有一个与该轴运动形式趋势一致的激光干涉仪测量，此时产生的某些线性误差和直线度误差会有一定的相关性，因此，本发明中定义一个相关系数ρ表示其中的关系； 

在激光干涉仪测量X向的六项基本误差时，并与此同时，在Y的方向，再加之一个激光干涉仪，运动趋势与X向运动一致，此时产生的Y项的6项基本误差就会产生一定的重叠项；X轴沿Y项的线性误差Δy_x与Y轴的定位误差Δy_y从空间来看二者是存在一定的关系，定义ρ＝Cov(Δy_x,Δy_y)便为二者的相关系数；一般情况下，设ρ＝Cov(ΔI_j,ΔJ_i)为误差与误差之间的相关系数，其中任意两项定位误差的相关系数为零；同理则可定义出其他误差项之间的相关性，矩阵： 

当量误差，由于机床几何误差最终体现在定位精度上，本发明中定义一种新的误差含义：即将空间误差量，投影到在各个轴线上的误差分量； 

其中： 

ΔX_x：X项上的当量误差； 

ΔY_y：Y项上的当量误差； 

ΔZ_z：Z项上的当量误差； 

最后得到当量误差方程： 

ΔX_x＝Δx_z-Δx_x-Δx_y-Δx_wd+zΔβ_x-zΔβ_wd+yΔγ_wd-zΔβ_y   (36) 

ΔY_y＝z[(Δα_x+Δα_y)-(Δy_x+Δy_y)]-x(Δγ_wd+Δγ_y+Δγ_xy)-Δy_wd+zΔα_wd  (37) 

ΔZ_z＝x(Δβ_wd+Δβ_y)+Δz_z+Δz_t+Δy_z+Δy_t-zΔα_z+yΔα_wd-Δz_wd   (38) 

步骤二：数控机床各几何误差的测量及其测量数据的整理 

激光干涉仪被频繁的用于机床误差检测上，本发明中，通过定点多测的方法在X，Y，Z三个方向上进行测量；分别在各轴50-600mm的行程上，以每20mm为一节点，进行测量重复9次并计算均值；只保留误差值： 

t_r＝T_r-D   (39) 

D：目标点； 

T_r：激光干涉仪测量值； 

t_r：误差值； 

使用垂直度测量仪测量机床的三项垂直度误差； 

定义各项几何误差均符合t_r～N(μ,σ²)均符合高斯分布的独立同分布； 

μ:为误差均值； 

σ²：为误差的方差； 

本发明中，为了验证预测及比较随机性效果，在各轴50-600mm的行程上，以每3mm为一节点随机记录下一组数据；表6～9为50-600mm的行程上，以每20mm为一节点，测量9次并取均值；由于篇幅限制，只列举数据的一部分 

表6 X轴几何误差测量值均值(mm) 

表7 Y轴几何误差测量值(mm) 

表8 Z轴几何误差测量值(mm) 

表9 单元间误差测量值(mm) 

步骤三：计算当量误差并利用随机过程对加工轴及面的随机性波动进行描述及预测 

步骤3.1计算当量误差并进行线条拟合 

本发明中，定义ΔX_x，ΔY_y，ΔZ_z被设为独立同分布的；根据实验数据的均值，可以计算出三向的当量误差；利用B-spline曲线进行数据在位置点的拟合；拟合原理如下： 

其中： 

u：表示当量误差； 

p:表示阶数； 

步骤3.2轴向随机性描述及预测原理 

针对其中一项误差的随机过程，可以将其称之为“高斯白噪声序列”，由白噪声过程定义可知，其中任意两点过程n₁,n₂两点的相关函数与其协方差函数相同均为σ²δ(n₁,n₂)，并且在移动过程中的任意时刻，均为不相关的，并且任意时刻均为N(0,σ²)于是在这个过程中得到任一点的概率密度函数为： 

其中：ΔX_xi：为某一方向上的当量误差； 

n_n：为某一方向上的位置点； 

本发明就ΔX_x，ΔY_y，ΔZ_z三向当量误差添加高斯白噪声序列，并以此来描 述和预测机床的几何误差非确定性波动，其波动范围在±3σ之间； 

步骤3.3在平面上的随机性描述及预测原理 

任意两个当量误差((ΔX_x,ΔY_y),(ΔY_y,ΔZ_z)和(ΔX_x,ΔZ_z))都是独立的随机变量，并且均符合N(0,σ²)分布；定义在X-Y平面上加工一平面，根据随机过程理论；可将平面上的任意点的误差点的误差预测为： 

{XY(n)＝ΔX_xcosωn+ΔY_ysinωn,n∈(-∞,+∞)}   (44) 

ΔX_x：X向当量误差； 

ΔY_y：Y向当量误差； 

ω：相对加工平面坐标系任一点与远点的矢量角； 

E_xy(n)属于联合高斯过程，从而也可以求出： 

E_xy(t)＝EΔX_x×cosωt+EΔY_y×sinωt＝0   (45) 

在机床运作过程中，任意两个过程点n₁,n₂时，可得到他们的相关函数与其协方差函数是相等的且均为： 

由于、各个点为独立同分布的，因此对于可得到相关系数： 

而ΔX_x，ΔY_y是服从N(0,σ²；0,σ²；cosω(n₁-n₂))，其二维密度函数为： 

依据本专利的方法，同样可以到得到在Y-Z、X-Z面的联合概率密度函数，其波动范围也应在±3σ之间； 

步骤四：关键误差识别与修改意见 

在本项发明的前面步骤中，已经提及过当量误差及波动预测的求解方法；当量误差及其波动作为空间误差项的反应结果，如何将对空间误差项影响较大的误差甄别出来，并减少波动范围就成为此步骤的重点；控制波动范围，最直观的方法是控制影响该项的方差，根据步骤1.4提出的均值误差模型则有： 

由于本项发明只针对机床的几何误差项则有： 

其中偏微分是用来具体识别出具体对于加工影响较大的误差项的，可将其就某一方向上展开归一化处理： 

m_ni的总量为1；在某一方向上的m_ni表示了该项误差对于结果影响的大小；图15～17分别表示了，在各个方向上对误差结果即非确定性波动范围影响较大的误差项；分别是在X向上，Δx_z,Δx_x,Δx_y,Δβ_x,Δβ_y；Y方向上，Δy_x,Δy_y,Δα_x,Δα_y,Δγ_xy；Z方向上，Δz_z,Δy_z,Δα_z,Δβ_y对加工结果有着较大影响，本发明为了更为直观的看到的其影响，并可以根据此原理来进行可削减波动范围的关键误差项识别工作。