CN102736558A - 基于时间序列算法的数控机床热误差实时补偿建模方法 - Google Patents

Abstract

一种属于精密机床加工技术领域的基于时间序列算法的数控机床热误差实时补偿建模方法。包括以下步骤：(1)数据零均值化预处理，用逆序检验法和峰度偏度检验法判定数据的平稳性和正态性；(2)用自相关函数和偏相关函数及其截尾性结果判定准则进行热误差数学模型的模式识别；(3)采用最小二乘估计法或长自回归计算残差法实现热误差数学模型的参数估计；(4)热误差数学模型的定阶，采用AIC定阶准则、F检验定阶准则和白度检验定阶准则相结合的判定方法，完成热误差数学模型的定阶；(5)综合判定条件的整合处理，构建完整的预测数学模型公式。本发明的建模方法硬件需求较低，适用性广，建立的模型具有较高的预测精度和可靠性。

Description

基于时间序列算法的数控机床热误差实时补偿建模方法

技术领域

本发明涉及的是一种精密机床加工技术领域的热误差补偿建模方法，具体地说，是一种基于时间序列算法的数控机床热误差实时补偿建模方法。

背景技术

国产数控机床行业经过多年的发展，在生产规模上取得了很大的进步，成功跻身机床生产大国的行列中。目前我国生产的数控机床约占国内数控机床市场份额的35％左右，但主要集中在低端和中低端市场，而中高端的数控机床市场则绝大部分被进口机床占领，特别是在高端数控机床领域中，国内机床产品的市场占有率仅占4％左右。随着当前工业产品向着小型化、精细化方向发展，数控机床对加工精密度的要求也越来越高，从而使得现代制造业对于高精度数控机床的需求也越来越大。我国每年在进口高端数控机床这一项上都花费了大量的外汇，而且国外对高端机床进口采取了多种限制，这就使得高端数控机床在我国的大量应用受到了一定限制，对整个制造业今后的发展也带来了不利因素。因而，通过有效的手段提高国产数控机床的加工精密度、稳定性和可靠度对我国制造业乃至整个民族工业的发展具有重要的现实意义。

目前，国产数控机床迈入高端领域的“瓶颈”问题主要是加工精度和可靠性难以达到国外同类产品的水准，而这涉及到多方面的问题：一、高端数控机床本身各种装配部件的加工精度由于缺乏更精密的机床作为工作母机而难以保证；二、机床硬件结构在考虑误差补偿方面的设计以及制造途径上考虑误差因素的加工还没有达到国外机床生产厂商的同等水平；三、国产的数控系统还不具备软件补偿技术，相关的研究还处于实验室阶段，未能在实际生产中进行实施。通常数控机床的误差补偿方法包括：一、根据实际加工所测试的误差数据，对数控加工程序进行人工调整；二、利用数控系统可提供的参数设定方式的误差补偿功能，将可以预估的误差数据提前输入对应的误差补偿设置项(如背隙补偿、螺距补偿和刀杆补偿等)，在实际加工中，数控系统将这些预设的误差项纳入过程计算进行补偿。然而，对于由机床温度场的变化而产生的热误差波动，通常的误差补偿方法无法进行预估，也就不能在机床加工过程中实时地对热误差做出补偿。而在数控机床的整体误差中，热误差所占的比例很大，在相当程度上会直接影响到工件的加工精度。因此，研发适合多类型、多规格、多品种的数控机床热误差实时补偿方法对我国数控机床产业的发展是非常必要和有益的。

经对现有技术文献检索发现，中国专利申请号：200410093428.1，专利名称为：基于机床外部坐标系偏置的数控机床误差实时补偿器。根据该发明提供的实时补偿器的结构和原理来看，主要论述了一种基于机床外部坐标系偏置的数控机床误差实时补偿器的构架方案。该补偿器采用基于单片机结构的计算处理模块、CNC接口和运动控制模块、传感器及变送模块来实现热误差的计算和补偿，通过外部计算机进行建模分析并将误差模型放入计算处理模块。但是，由于数控机床实际加工中的热误差不但与机床结构相关联，而且受到各种加工条件和环境因素的影响，如切削速度、切深、冷却液类型、工件材料类型、加工周期、环境温度变化等，并且由于机床加工中温度变化的缓慢性，机床的热误差多表现为交叉作用和非线性。仅根据多元回归理论用最小二乘原理进行热误差建模分析，难以得到准确的实用化模型。中国专利申请号：200710045903.1，专利名称为：数控机床定位误差实时补偿装置。根据该发明提供的实时补偿装置的结构和原理来看，主要论述了一种集成了计算处理模块、温度传感器采集和变送模块、数控接口控制模块的机床定位误差补偿装置，采用和数控机床的PMC之间的数据交互，来实现定位误差的补偿工作，但是，该发明没有描述补偿装置中定位误差数学模型的建模方法，只是对模型的输入信号和输出内容进行了阐述。李永祥等在期刊《四川大学学报》第38卷第2期发表了《数控机床热误差的时序分析法建模及其应用》，该文章论述了利用实测的热误差序列进行时序分析识模、建模和预报的方法，并提出了部分的判别公式。但是，文中提出的时序建模方法缺少了时序建模的关键判定条件，并且采用单一判定准则难以保证模型的精度。第一，文中仅对热误差数据序列的平稳性采用统计方差的方法进行了检验，没有进行样本的零均值法处理，并且没有指出具体的平稳性检验的具体方法，同时也未对零均值法处理后的样本进行正态性的相关判定；第二，文中直接采用了AR模型作为热误差系统的拟合模型，没有根据样本序列先进行模型的模式识别，因而可能由于模型选错而使得最后的数学模型达不到补偿效果；第三，模型定阶式只采用了单一的FPE准则，而单一定阶准则判定时，因为所遇到的不确定因素较多，有时会出现模棱两可的局面，增加了阶数辨识的困难和不稳定性。

发明内容

本发明的目的在于克服现有数控机床热误差补偿中数学模型建模方法的不足，提供一种基于时间序列算法的数控机床热误差实时补偿方法，使本发明的补偿方法可靠、全面。

由一串随机变量…，x₁，x₂，x₃，…构成的序列叫做随机序列，用x_t(t＝…，1，2，3，…)或{x_t}表示。如果下标t是整数变量，它代表着等间隔的时刻增长量，如第t时、第t天、第t次等，我们就称这种随机序列为时间序列，而整数变量t即认为是指某时刻。对于时间序列x₁，x₂，x₃，…，x_N，其数据的特点是承认观测值排列顺序的重要性，正是数据的顺序与大小反映了数据所包含的信息，反映了数据内部的相互联系。事实上，所获得的数据最为重要和有用的特性，就是观测值之间的依赖关系或相关性，正是这种相关性表征了产生这些数据的现象、过程、系统的“动态”或“记忆”。这种相关性一旦被定量地描述出来，就可以从系统的过去值预测其将来值。

用来分析各种相依有序的离散数据集合的方法称为时间序列算法。从表面上看，时序分析撇开了系统变量内在因果关系和结构关系的影响。事实上，时序中反映了曾经发生过的所有因果关系和结构关系的影响。时序分析是从总的方面进行考察，来说明各种元素的综合作用。因此，当我们所关心的影响因素错综复杂或有关的数据资料无法得到时，就直接采用时间t作为变量来综合地代替这些因素。时间作为一个明确的变量进入模型，其意义表面上是表示因变量随时间而自发的变化，而实际上是代表了决定因变量变化的诸因素的联合影响。时间序列算法实际上是一种处理动态数据的参数化时域分析方法，是指对观测数据拟合一个参数模型，再利用这个模型对观测数据及产生这一数据的系统进行分析，以便更本质地了解数据的内在结构和系统的动态特征，从而可以利用过去的观测数据对未来值进行预测和控制。

数控机床热误差形成的原因错综复杂，包括机床上各种热源影响、外部环境影响、人为影响等众多因素，而这些影响因素之间又相互干扰和耦合。通常对热误差进行建模时，除了要测试机床各运动轴的热误差数据外，还需要测试热源点的温度数据值，这就需要在机床的热源位置布置温度传感器，对于多轴联动的数控机床，由于热源点位置较多，因而布置的温度传感器也需要很多，这使得建模测试过程变得复杂。机床上温度传感器的布线一般需要较长的电缆，其容易受到机床内部的电、磁干扰而使得测试数据产生偏差；热电偶或热电阻式传感器本身的测试精度也容易受到如安装位置，环境温度以及人为因素的干扰而产生误差；由于温度测试需要经过变送器和A/D转换模块，因而整个温度测试环节中经过了多次信号转换，不可避免的会有系统误差存在；此外，由于热源之间存在一定的相关性，因此测试所得的温度数据值之间也就存在耦合，要区分这种耦合性，需要模型在处理温度数据时另行设置相应的分析方法，这就增加了建模的难度。而热源点温度变化与热误差之间的变化一致性，也是必须要予以考虑的问题，因为温度变化与热误差变化如果不一致，则会使得模型的准确性难以得到保证，而这种一致性问题也是通常热误差建模方法的难点。

本发明论述了基于时序法的热误差补偿建模策略，只需要采集机床各运动轴的热误差信息，然后通过对各控制热误差动态历史数据的分析和归纳来提取数据中所包含的信息，根据数据序列的变化规律来分析机床系统的特性，并通过有效的建模方法拟合热误差变化的规律性，从而能够有效推断和预测机床未来的热变形情况，这样既减少了建模分析和计算的复杂程度，也降低了采用温度采集方案的硬件成本，对于提高数控机床热误差补偿的建模精度和实用性推广的研究工作具有实际意义。

本发明是通过以下技术方案来实现的，本发明包括：

(1)数据预处理，采用零均值法导出建模用的时序数据，用逆序检验法判定序列的平稳性，用峰度偏度检验法判定热误差数据序列的正态性；

(2)热误差数学模型的模式识别，对于满足平稳性和正态性的样本用自相关函数和偏相关函数及其截尾性结果进行判定，完成热误差数学模型的选型；

(3)热误差数学模型的参数估计，采用最小二乘估计法或长自回归计算残差法，实现对模型展开式中系数的估计；

(4)热误差数学模型的定阶，采用AIC定阶准则、F检验定阶准则和白度检验定阶准则相结合的判定方法，完成模型阶次的判定；

(5)综合判定条件的整合处理，构建完整的预测数学模型公式。

所述的数据预处理，是指：为了保证进行模式识别的动态数据为平稳数据，对采样获得的热误差数据先后进行零均值处理、平稳性和正态性判定，以得到满足模型识别和定阶要求的时序。

所述的零均值处理法，其方法为：估计测试所得样本的均值，再用每个样本数据减去样本均值，以差值序列作为建模用的时序。

所述的逆序检验法，其方法为：零均值处理后序列x_t＝{X(n)；n＝0，1，…，N-1}，当N充分大时，取N＝kM，M为一较大的正整数，按长度M将上述序列分为k个等长的子序列：

$\begin{array}{cccc}{x}_{11},& x_{12},& \·\·\·,& x_{1M};\\ x_{21},& x_{22,}& \·\·\·,& x_{2M};\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ x_{k1},& x_{k2},& \·\·\·,& x_{\mathrm{kM}};\end{array}$ (i＝1，2，…，k；j＝1，2，…，M)

对k个子序列计算均值

于是有序列：

定义随机变量：则统计量

给出了上式中按大小排列的逆序总次数。

如果总序列X(n)是平稳的，则

应是独立分布的，

这时应有：

，则统计量

渐进服从N(0，1)分布。

由此可对逆序数A进行统计检验，从而达到检验序列X(n)是否平稳的目的。对于给定的显著性水平α，查表得u_α/2＝1.96，如果|u|＜u_α/2，就可以认为总序列X(n)是平稳的。反之就是不平稳的。

所述的峰度偏度检验法，其方法为：对零均值处理后的序列检验其三阶矩和四阶矩的性质。因为对于正态随机变量X□N(μ，σ²)有：

偏度： $E{\left(\frac{X-\mathrm{\μ}}{\mathrm{\σ}}\right)}^3=0,$ 峰度： $E{\left(\frac{X-\mathrm{\μ}}{\mathrm{\σ}}\right)}^4=3$

偏度反映了概率密度函数的对称性，峰度反映了概率密度曲线的状态。

对于长度为N的样本序列x_t＝{X(n)；n＝0，1，…，N-1}，令：

$\stackrel{\‾}{x}=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(n\right),$ $S^2=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{(x\left(n\right)-\stackrel{\‾}{x})}^2$

$m_3=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{\left(\frac{x\left(n\right)-\stackrel{\‾}{x}}{S}\right)}^3,$ $m_4=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{\left(\frac{x\left(n\right)-\stackrel{\‾}{x}}{S}\right)}^4-3$

当N充分大时，如果m₃和m₄近似于零，即可以认为序列{x_t}具有正态性。

所述的热误差数学模型的模式识别，是指：对热误差数据序列选择匹配的拟合模型，再利用该模型对数据特征和产生这一数据的系统进行分析，从而可以利用过去的观测数据对未来值进行预测和控制。时间序列法的常用拟合模型包括自回归(AR)模型、滑动平均(MA)模型和自回归滑动平均(ARMA)模型，具体实现如下：

首先，求取数据序列的自相关函数ρ_k，并通过截尾性判断初步识别模型特征。用γ_k除以方差γ₀便得滞后k时的自相关函数，即ρ_k＝γ_k/γ₀，ρ₀≡1。

γ_k或ρ_k表示x_t对x_t-k的相关性。γ_k和ρ_k又称为理论自协方差函数和理论自相关函数。当我们根据一个样本的数据x₁，x₂，…，x_N去求γ_k和ρ_k的估计值时，则有如下定义：

${\hat{r}}_k={\hat{r}}_{-k}=\frac{1}{N}\underset{t=1}{\overset{N-k}{\mathrm{\Σ}}}{x}_t{x}_{t-k},$ k＝0，1，2，…，N-1

称为样本自协方差函数。其中x_t是减掉均值以后的观测数据。而

k＝0，1，2，…，N-1称为样本自相关函数。而AR(p)模型的自相关函数ρ_k是由指数衰减和衰减正弦波组成，即ρ_k的尾部不可能在延迟某步之后等于零，而是按负指数律衰减，具有“拖尾”的特性；MA模型的自相关函数ρ_k，从k＞q以后全部为零。具有“截尾”的特性，即若平稳时间序列的自相关函数为截尾，则其必是MA(q)序列；ARMA模型中的q是滑动平均的记忆部分，而当k＞q时，滑动平均部分就不起作用了，因此ARMA(p，q)模型的自相关函数ρ_k也具有拖尾的性质。

然后，求取数据序列的偏相关函数

进一步识别序列的模型特征。由于截尾性是MA序列特有的标志，据此可以推测出MA(q)模型的阶数。而AR(p)过程则没有这样的特性，采用偏相关函数可使AR(p)过程在滞后p步后截止。

设{x_t}为平稳、零均值序列，考虑x_t≡1，x_t≡2，…，x_t≡k对x_t的线性最小方差估计，即选择系数

使得：

达到极小。为此目的，对δ的右端分别求偏导数

j＝1，2，…，k，并令其为零，便得到

应该满足的Yule-Walker线性方程组，其系数矩阵为Toeplitz矩阵。

对k＝1，2，…，利用Cramer法则，可得：

(k＝1，2，…)即为x_t的偏相关函数。对于p阶自回归过程AR(p)，出现了序列的截尾特性：

即表明AR(p)模型的

是p步截尾的。

最后，根据前面数据序列自相关函数ρ_k和偏相关函数

的计算，可识别出系统可用哪种模型来拟合：(1)自相关函数ρ_k截尾而偏相关函数拖尾，系统可采用MA(q)模型拟合；(2)自相关函数ρ_k拖尾而偏相关函数是p步截尾的，系统可采用AR(p)模型拟合；(3)自相关函数ρ_k和偏相关函数

都为拖尾时，系统可采用ARMA(p，q)模型拟合。

所述的热误差数学模型的参数估计，是指：对于模型展开式

根据时序{x_t}按某一方法估计出

和

这n+1个参数。而由于对模型而言，{a_t}是模型的残差序列，其方差为：

可见，估计出

后，即可按上式估计出

所以，通过算法估计出

即可实现对模型全部参数的估计。

所述的最小二乘估计法，其方法为：将序列{x_t}(t＝1，2，…，N)代入AR(n)模型展开式可得到以下线性方程组：

用矩阵形式表示为：

$y=\left[\begin{array}{c}{x}_{n+1}\\ x_{n+2}\\ \·\\ \·\\ \·\\ y_N\end{array}\right],$

$a=\left[\begin{array}{c}{a}_{n+1}\\ a_{n+2}\\ \·\\ \·\\ \·\\ a_N\end{array}\right],$ $x=\left[\begin{array}{cccc}{x}_n& x_{n-1}& \·\·\·& x_1\\ x_{n+1}& x_n& \·\·\·& x_2\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ x_{N-1}& x_{N-2}& \·\·\·& x_{N-n}\end{array}\right]$

根据多元回归理论，参数矩阵的最小二乘估计为：

所述的长自回归计算残差法，用于对ARMA(n，m)模型进行参数估计，具体实现如下：第一，按AR模型的参数估计方法对序列{x_t}拟合高阶(长)AR{p}模型{p≥n+m}；第二，根据AR{p}模型计算残差：

(t＝p+1，p+2，…，N)

由此得到N-p个残差值a_p+1，a_p+2，…，a_N即{a_t}(t＝p+1，p+2，…，N)；第三，将{a_t}作为ARMA模型中滑动平均部分的观测值，连同{x_t}一起代入ARMA(n，m)模型，由于此时{x_t}、{a_t}均已知，则可按合适的AR模型的参数估计方法估计出ARMA模型参数：

写成矩阵方程为：Y＝Xβ+a

其中，

Y＝[x_p+n，x_p+n+1，…x_N]^T

$X=\left[\begin{array}{ccccccccc}{x}_{p+n-1}& x_{p+n-2}& \·\·\·& x_p& \begin{array}{c}\·\\ \·\\ \·\end{array}& a_{p+n-1}& a_{p+n-2}& \·\·\·& a_{p+1}\\ x_{p+n}& x_{p+n-1}& \·\·\·& x_{p+1}& \begin{array}{c}\·\\ \·\\ \·\end{array}& a_{p+n}& a_{p+n-1}& \·\·\·& a_{p+2}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \begin{array}{c}\·\\ \·\\ \·\end{array}& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ x_{N-1}& x_{N-2}& \·\·\·& x_{N-n}& \begin{array}{c}\·\\ \·\\ \·\end{array}& a_{N-1}& a_{N-2}& \·\·\·& a_{N-n+1}\end{array}\right]$

因此β的最小二乘估计为：β＝(X^TX)^-1X^TY。

所述的热误差数学模型的定阶，是指：通过结合3种定阶方法，采用AIC和F准则来为热误差模型初步定阶，再按照白度检验的定阶准则，通过χ²检验法检验残差序列是否为白噪声，进行模型的阶数分析和确认，具体实现如下：

首先，初步拟定模型迭代阶数为1，并根据AIC定阶准则求取AIC最小值对应的阶数，作为模型的初始定阶数；

然后，根据F检验定阶准则，判断模型初定的阶数是否合适，如F检验不显著，则表明初定的模型阶数是合适的；反之，如F检验显著，则表明模型阶数不合适，需要增加阶数，用AIC定阶准则继续往高阶识别；

最后，对于F检验不显著的模型，根据白度检验定阶准则，通过χ²检验法判定模型的残差序列{a_t}是否为白噪声，如残差序列不是白噪声，则需要增加模型阶数，用AIC定阶准则继续往高阶识别。

所述的AIC定阶准则，是指：在模型参数的极大似然数基础上，根据最小信息准则，对模型的阶数和相应的参数，同时给出一种最佳估计。定义AIC准则函数如下：

式中：N为样本容量，

为模型残差的方差，p为模型的阶数，对不同的p值计算出相应的AIC(P)值，AIC值最小的p就是适用模型的阶数。

所述的F检验定阶准则，其方法为：对序列{x_t，t＝1，2，…，N)所拟合的AR(p)或ARMA(p，q)模型的残差平方和：是阶数n＝p或n＝p+q的递减函数，S(n)随着阶数的增加而减小，最佳阶次是使S(n)减小不显著时的那个n值。

为了判定当模型的阶数改变时，相应的残差平方和的变化是否显著，引入检验的统计量：

式中：S(n₁)——参数个数为n₁模型的残差平方和；S(n₂)——参数个数为n₂模型的残差平方和；N——观测数据数目。

数理统计的结论表明，当N足够大时，统计量t渐进地服从F(n₂-n₁，N-n₂)分布，自由度为(n₂-n₁，N-n₂)。根据以上原理，ARMA模型定阶的F检验准则为：

记S₀为(高阶)模型ARMA(2p，2p-1)的残差平方和，S₁为(低阶)模型ARMA(2p-2，2p-3)的残差平方和，而两模型参数个数之差为4。当N充分大时，对于比较ARMA(2p-2，2p-3)与ARMA(2p，2p-1)而言，近似地有：

给定显著性水平α(α＝0.05)，查F分布可得临界值F_α。若F＜F_α，则F检验不显著，所选模型的阶数合适；若F≥F_α，则F检验显著，说明需增加阶数，继续往高阶识别。

所述的白度检验定阶准则，是指：通过χ²检验法判定模型的残差序列{a_t}是否为白噪声序列，如所选模型的残差近似白噪声，说明模型接近于包含给定的时间序列所提供的全部信息。反之，残差不是白噪声，就意味着残差中含有有用信息，应增大模型阶次后重新估计所有参数，以便将残差中的有用信息挖掘出来，直至{a_t}成为白噪声序列为止。其方法为：设a_t＝{a(n)；n＝0，1，2，…，N-1}是白噪声序列，令

$\widehat{r_k}(N,a)=\frac{1}{N}\underset{t=1}{\overset{N-k}{\mathrm{\Σ}}}{a}_t{a}_{t+k}$ (k＝1，2，…，K)

$\widehat{{\mathrm{\ρ}}_k}(N,a)=\widehat{r_k}(N,a)/\widehat{r_0}(N,a)$ (k＝1，2，…，K)

为

的估计。当N□K时，可近似地看成是K个独立的N(0，1)分布所组成的随机向量。于是a(n)的独立性检验就转化为(k＝1，2，…，K)是否是N(0，1)分布总体的K个样本序列，这可用χ²检验法来检验。由于K个独立N(0，1)随机变量的平方和服从自由度为K的χ²分布，为此可作假设：H0：

k＝1，2，…K是K个独立的N(0，1)分布随机变量。也就是说：

$Q_k=N\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\widehat{{\mathrm{\ρ}}_k^2}(N,a)$

是自由度为K的χ²分布。因此，如果

则假设H₀成立；反之假设H₀不成立，因而否定k＝1，2，…K是K个独立N(0，1)分布的随机变量的假设，从而也就否定了序列{a_t}是白噪声序列。实际应用中，当N在数百以上时，K可取20～30，显著性水平α常用0.05。

所述的综合判定条件的整合处理，是指：根据模型定阶结果，将估计参数代入预选的模型系统，进行数据整合后得到基于时序算法分析和计算的数控机床热误差数学模型。

本发明将时间序列算法和计算机自动建模技术应用于数控机床的热误差补偿建模方法中，解决了单一最小二乘原理拟合建模由于热误差的非线性而产生的偏差和缺陷。同时为数控机床热误差补偿方面的研究人员进行热误差补偿模型优化和补偿实效性的深入研究提供了有效的分析手段和应用模型。该方法采用零均值法、逆序检验法和峰度偏度检验法进行热误差采集数据的预处理，以自相关函数和偏相关函数及其截尾性判定识别拟合模型的类型，采用最小二乘估计法或长自回归计算残差法实现模型的参数估计，并采用AIC定阶准则、F检验定阶准则和白度检验定阶准则相结合的判定方法，对模型进行精确定阶。该建模方法能够高效、准确地完成对数控机床热误差的建模计算，并且获得的数学模型公式能够精确地推断和预测机床热误差的变化情况，使机床的热误差补偿工作更为可靠和有效。

附图说明

图1为本发明建模方法的工作流程示意图

具体实施方式

下面结合附图对本发明的实施例进行说明。本实施例在以发明技术方案为前提下进行实施，给出了详细的实施方式和具体过程，但本发明的保护范围不限于下述的实施例。

本发明所述的热误差实时建模方法是一种基于时间序列算法的推理方法，具体实施过程如下：

对一台车削加工中心的主轴进行热误差补偿模型的建模分析，主轴X和Z向的热误差数据通过固定在刀架上的两个涡流位移传感器来测量，并通过数据自动采集平台按照10Hz的采样速率自动记录两个方向上的热漂移误差值。在空切削条件下实时进行测试：机床主轴旋转、拖板移动和冷却液流动而无切削加工，机床先运行4小时，接着模拟中午休息停机0.5小时，然后再运行4小时后停机0.5小时。在各种速度条件下重复测试过程，X和Z轴方向各取得10组数据。然后按照下面的步骤对取得的数据进行建模计算：

1.数据预处理

采用零均值法处理样本数据，然后用逆序检验法判定序列的平稳性，用峰度偏度检验法判定热误差数据序列的正态性，对不符合平稳性要求或正态性要求的数据序列分别进行平稳化处理或正态化处理，舍弃经处理后仍然达不到平稳性或正态性要求的数据序列。

(1)对X和Z轴的各10组热误差数据作零均值处理，使经处理的样本序列满足随机过程零均值的要求，即保证了自协方差函数与自相关函数是一致的。

(2)用逆序检验法判定序列的平稳性，每组数据的样本大小为32000左右，经检验X向有3组数据满足平稳性要求，Z向有5组数据满足平稳性要求。对于给定的显著性水平α(u_α/2＝1.96)，最终计算得到X方向3组数据的统计值为u＝0.25549，0.21417，0.29834，3组数据的u值都小于u_α/2，Z方向5组数据的统计值为u＝0.37132，0.3578，0.33312，0.4078，0.35269，5组数据的u值都小于u_α/2。X向其余7组数据和Z向其余5组数据经平稳化处理后分别有4组和2组可达到平稳性要求，不过计算出的统计值都较大，因而不适宜作为建模的数据样本。

(3)用峰度偏度检验法判定序列的正态性，经检验X向有2组数据满足正态性要求，Z向有2组数据满足正态性要求。而X向其余1组数据和Z向其余3组数据经正态化处理后各有1组数据可以达到正态化要求，但计算出的m₃和m₄值偏离0值较大，因此不适宜作为建模用的样本数据。经峰度偏度计算，X方向2组数据的峰度值m₃和偏度值m₄分别为：m₃＝0.0237138，0.100529；m₄＝0.0175181，0.062136。Z方向2组数据的峰度值m₃和偏度值m₄分别为：m₃＝0.0332945，0.041001；m₄＝0.024871，0.035924。m₃和m₄越接近于零，则表明序列的正态性越好，因此，最终选取峰度偏度均最小的X方向和Z方向数据各1组作为样本序列，再进行后续的模式识别。

2.热误差数学模型的模式识别

对经过预处理筛选后得到的样本序列用自相关函数和偏相关函数及其截尾性结果进行模式识别，完成热误差数学模型的选型。

(1)求取数据序列的自相关函数ρ_k，序列的自相关系数的极值接近于1，但是在整个极值右侧的动态数据分布都是在0值上下波动，没有明显下降的趋势。因而从动态数据的分布趋势来看，也不难发现其自相关系数具有拖尾的特性。

取M＝样本总量/10，即3200。且

从

开始到计算其中满足的个数判别如下：

①X方向序列：813个，813/M＝25.4％＜68.3％，所以不截尾，即具有拖尾性；

②Z方向序列：1077个，1077/M＝33.7％＜68.3％，所以不截尾，即具有拖尾性；以上实验数据说明序列具有拖尾性。因而不能用滑动平均模型MA(q)进行拟和。

(2)求取数据序列的偏相关函数

进一步识别序列的模型特征。同样取M＝样本总量/10，即3200。逐个检验

计算其中满足

的个数判别如下：

①X方向序列：1345个，1345/M＝42％＜68.3％，所以不截尾，即具有拖尾性；

②Z方向序列：982个，1077/M＝30.7％＜68.3％，所以不截尾，即具有拖尾性；

以上说明信号数据序列具有拖尾性。因而不能用自回归模型AR(p)模型进行拟和。根据上面计算的结果，可以看出对于机床主轴的X和Z方向热漂移误差都应该以自回归滑动平均模型ARMA(p，q)作为拟和模型。

3.热误差数学模型的定阶和参数估计

采用AIC定阶准则、F检验定阶准则和白度检验定阶准则相结合的判定方法，完成模型阶次的判定；并对模型采用长自回归计算残差法，实现对模型展开式中系数的估计。

(1)设迭代阶数1为模型的初始定阶数，考虑到数学模型的实用性，如模型阶数过高，则不利于实际使用中的在线计算，因此设定模型定阶上限为5。根据AIC定阶准则求取AIC最小值对应的阶数，AIC检验准则的特性是：p值在设定范围内逐渐增大时，AIC的计算值呈现由大变小，然后再变大的规律，因而在运用AIC准则时，它们的最小值对应的阶数，即为模型的最佳阶数。经第一次循环计算AIC值，得出的最小AIC值为-180.5，最终求得的自回归部分的阶数p为2，最终求得的滑动平均部分的阶数q为1，模型为ARMA(2，1)；

(2)根据F检验定阶准则，判断模型初定的阶数是否合适。根据F准则可知，对于给定显著性水平α＝0.05，查F分布可得临界值：F_0.05(4，4923)＝2.37，在初定的ARMA(2，1)模型的基础上，F分布值为：F＝3.2＞F_0.05＝2.37，据F准则判定原则，该模型初定的阶数是不适合的，因此需增加模型的阶数，继续往高阶识别。增加阶数为p+1，则增阶后的模型为ARMA(3，2)，然后继续进行F准则检验，计算得：F分布值F＝1.9815＜F_0.05＝2.37，说明模型：ARMA(3，2)是合适的(F检验不显著)；

(3)对模型展开式进行参数估计：根据F检验定阶准则增阶判定后得出的模型ARMA(3，2)，采用长自回归计算残差法，求取一组自回归系数

残差序列{a_t}以及残差的方差

(4)根据白度检验定阶准则，通过χ²检验法判定模型的残差序列{a_t}是否为白噪声序列，如模型ARMA(3，2)的残差近似白噪声，说明模型接近于包含给定的时间序列所提供的全部信息。反之，残差不是白噪声，就意味着残差中含有有用信息，应增大模型阶次后重新估计所有参数，以便将残差中的有用信息挖掘出来，直至{a_t}成为白噪声序列为止。对ARMA(3，2)模型进行白度检验所示。计算得：导入自相关系数的容量为32000，K值取

约等于57，index值为N-5＝31995，Q_K＝64.109。对显著性水平α＝0.05，

则有：因此，对前面已经F检验的结果，残差序列不是白噪声序列，这意味着残差中仍含有有用信息，应该继续增大模型阶次的估算范围重新进行定阶；

(5)第二轮分析计算得最小AIC值为-187.6，，最终求得的自回归部分的阶数p为4，滑动平均部分的阶数q为3。模型修正为ARMA(4，3)；

(6)第二轮进行F准则检验后F分布值的结果为：F＝1.4192＜F_0.05＝2.37。经F准则判定后说明由前面由AIC准则确立的模型：ARMA(4，3)是合适的；

(7)再次对模型展开式进行参数估计：根据第二轮F检验定阶准则判定后得出的模型ARMA(4，3)，采用长自回归计算残差法，求取相关的回归系数

残差序列{a_t}以及残差的方差

(8)同样对由AIC和F准则得出的模型ARMA(4，3)再进行白度检验。自相关系数的容量为32000，K值取

约等于57，index值为N-7＝31993，Q_K＝57.337。对显著性水平α＝0.05，

则有：

表明的残差序列是白噪声序列。因此，由上述检验初步确定模型的阶数，即模型可定为：ARMA(4，3)。

4.热误差数学模型的整合

根据模型定阶结果，将参数估计结果代入预选的模型，进行数据整合后得到基于时序算法分析和计算的数控机床热误差数学模型如下：

(1)X方向的数学模型

δ_x＝3.922×10^-6x⁴+1.051×10^-5x³+11.836×10^-4x²-1.243×10^-6x+0.6185+(-3.241×10^-2-2.439×10^-3×ΔT)x

式中，δ_x为X方向的热误差，ΔT为温度修正系数。

(2)Z方向的数学模型

δ_z＝1.356×10^-7z⁴+2.497×10^-6z³-2.331×10^-4z²-3.834×10^-3z-4.4072+(-1.5842×10^-2-7.1508×10^-3×ΔT)z

式中，δ_z为Z方向的热误差，ΔT为温度修正系数。

Claims (9)

1.一种基于时间序列算法的数控机床热误差实时补偿建模方法，其特征在于，包括以下步骤：(1)数据预处理，采用零均值法导出建模用的时序数据，用逆序检验法判定序列的平稳性，用峰度偏度检验法判定热误差数据序列的正态性；(2)热误差数学模型的模式识别，对于满足平稳性和正态性的样本用自相关函数和偏相关函数及其截尾性进行判定，完成热误差数学模型的选型；(3)热误差数学模型的参数估计，采用最小二乘估计法或长自回归计算残差法，实现对模型展开式中系数的估计；(4)热误差数学模型的定阶，采用AIC定阶准则、F检验定阶准则和白度检验定阶准则相结合的判定方法，完成模型阶次的判定；(5)综合判定条件的整合处理，构建完整的预测数学模型公式；

所述的数据预处理，是指：对采样获得的热误差数据先后进行零均值处理、平稳性和正态性判定，保证进行模式识别的动态数据为平稳数据；

所述的热误差数学模型的模式识别，是指：对热误差数据序列选择匹配的拟合模型，再利用该模型对数据特征和产生这一数据的系统进行分析；

所述的热误差数学模型的参数估计，是指：对于模型展开式

根据时序{x_t}估计出和这n+1个参数；

所述的热误差数学模型的定阶，是指：采用AIC和F准则对模型初步定阶，再用白度检验准则，通过χ²检验法检验残差序列是否为白噪声，完成对模型的定阶；

所述的综合判定条件的整合处理，是指：根据模型定阶结果，将估计参数代入预选的模型系统，进行数据整合后得到基于时序算法分析和计算的数控机床热误差数学模型。

2.根据权利要求1所述的逆序检验法，其特征是，对零均值处理后的误差序列，取N＝kM，按长度M将上述序列分为k个等长的子序列：对k个子序列计算均值

得到序列：定义随机变量：

统计量给出了上式中按大小排列的逆序总次数，如序列X(n)是平稳的，则

(i＝1，2，…，k)应是独立同分布的，这时应有：

则统计量

渐进服从N(0，1)分布，由此可对逆序数A进行统计检验，从而达到检验序列X(n)是否平稳的目的；对于给定的显著性水平α，查表得u_α/2＝1.96，如果|u|＜u_α/2，就可以认为总序列X(n)是平稳的；反之就是不平稳的。

3.根据权利要求1所述的峰度偏度检验法，其特征是，对零均值处理后的误差序列检验其三阶矩和四阶矩的性质。对于正态随机变量X□N(μ，σ²)有：

偏度：

峰度：

偏度反映了概率密度函数的对称性，峰度反映了概率密度曲线的状态；对于长度为N的样本序列x_t＝{X(n)；n＝0，1，…，N-1}，

$m_3=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{\left(\frac{x\left(n\right)-\stackrel{\‾}{x}}{S}\right)}^3,$ $m_4=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{\left(\frac{x\left(n\right)-\stackrel{\‾}{x}}{S}\right)}^4-3$

如果m₃和m₄近似于零，即可以认为序列{x_t}具有正态性。

4.根据权利要求1所述的自相关函数及其截尾性判定，其特征是，求取序列的自相关函数ρ_k，并通过截尾性判断识别模型特征；用γ_k除以方差γ₀便得滞后k时的自相关函数，即ρ_k＝γ_k/γ₀，ρ₀≡1。γ_k或ρ_k表示x_t对x_t-k的相关性；根据x₁，x₂，…，x_N去求γ_k和ρ_k的估计值，k＝0，1，2，…，N-1为样本自协方差函数，其中x_t是减掉均值以后的观测数据，而

k＝0，1，2，…，N-1为样本自相关函数；AR(p)模型的自相关函数ρ_k是由指数衰减和衰减正弦波组成，ρ_k的尾部按负指数律衰减，具有“拖尾”的特性；MA模型的自相关函数ρ_k，从k＞q以后全部为零。具有“截尾”的特性；ARMA模型中的q是滑动平均的记忆部分，而当k＞q时，ARMA(p，q)模型的自相关函数ρ_k也具有拖尾的性质。

5.根据权利要求1所述的偏相关函数及其截尾性判定，其特征是，求取数据序列的偏相关函数

进一步识别序列的模型特征；设{x_t}为平稳、零均值序列，考虑x_t≡1，x_t≡2，…，x_t≡k对x_t的线性最小方差估计，即选择系数

使得

达到极小；对δ的右端分别求偏导数

j＝1，2，…，k，并令其为零，便得到

应该满足的Yule-Walker线性方程组，其系数矩阵为Toeplitz矩阵。对k＝1，2，…，可得，

(k＝1，2，…)即为x_t的偏相关函数。对于p阶自归过程AR(p)，出现了序列的截尾特性，即表明AR(p)模型的

是p步截尾的。

6.根据权利要求1所述的长自回归计算残差法，其特征是，对ARMA(n，m)模型进行参数估计，按AR模型的参数估计方法对序列{x_t}拟合高阶(长)AR{p}模型{p≥n+m}；根据AR{p}模型计算残差：

(t＝p+1，p+2，…，N)得到N-p个残差值a_p+1，a_p+2，…，a_N，即{a_t}(t＝p+1，p+2，…，N)；将{a_t}作为ARMA模型中滑动平均部分的观测值，连同{x_t}一起代入ARMA(n，m)模型，由于此时{x_t}、{a_t}均已知，则可按合适的AR模型的参数估计方法估计出ARMA模型参数：Y＝Xβ+a，因此β的最小二乘估计为：β＝(X^TX)^-1X^TY。

7.根据权利要求1所述的AIC定阶准则，其特征是，在模型参数的极大似然数基础上，根据最小信息准则，对模型的阶数和相应的参数同时给出一种最佳估计；AIC准则函数

N为样本容量，

为模型残差的方差，p为模型的阶数，对不同的p值计算出相应的AIC(P)值，AIC值最小的p就是适用模型的阶数。

8.根据权利要求1所述的F检验定阶准则，其特征是，对误差序列所拟合的模型的残差平方和

是阶数n＝p或n＝p+q的递减函数，S(n)随着阶数的增加而减小，最佳阶次是使S(n)减小不显著时的那个n值；引入检验的统计量

式中S(n₁)是参数个数为n₁模型的残差平方和，S(n₂)是参数个数为n₂模型的残差平方和，N是观测数据数目；当N足够大时，统计量t渐进地服从F(n₂-n₁，N-n₂)分布；记S₀为(高阶)模型的残差平方和，S₁为(低阶)模型的残差平方和，而两模型参数个数之差为4。近似地有：

□F(4，N-(4p-1))给定显著性水平α(α＝0.05)，查F分布可得临界值F_α。若F＜F_α，则F检验不显著，所选模型阶数合适；若F≥F_α，则F检验显著，说明需增加阶数，继续往高阶识别。

9.根据权利要求1所述的白度检验定阶准则，其特征是，通过χ²检验法判定模型的残差序列{a_t}是否为白噪声序列，所选模型的残差近似白噪声，说明模型接近于包含给定的时间序列所提供的全部信息；反之，意味着残差中含有有用信息，应增大模型阶次后重新估计所有参数，直至{a_t}成为白噪声序列为止；其方法为设a_t＝{a(n)；n＝0，1，2，…，N-1}是白噪声序列，令 $\widehat{r_k}(N,a)=1/N\underset{t=1}{\overset{N-k}{\mathrm{\Σ}}}{a}_t{a}_{t+k}(k=\mathrm{1,2},\·\·\·,K)$ $\widehat{{\mathrm{\ρ}}_k}(N,a)=\widehat{r_k}(N,a)/\widehat{r_0}(N,a)(k=\mathrm{1,2},\·\·\·,K)$

为

的估计；当N□K时，

可近似地看成是K个独立的N(0，1)分布所组成的随机向量，a(n)的独立性检验就转化为

(k＝1，2，…，K)是否是N(0，1)分布总体的K个样本序列，这可用χ²检验法来检验；由于K个独立N(0，1)随机变量的平方和服从自由度为K的χ²分布，为此可作假设：

k＝1，2，…K是K个独立的N(0，1)分布随机变量；即

是自由度为K的χ²分布；如则假设H₀成立；反之假设H₀不成立，从而也就否定了序列{a_t}是白噪声序列。实际应用中，显著性水平α常用0.05。

Images