CN103358180B - 一种快速辨识机床主轴热态特性的方法 - Google Patents

Abstract

本发明公开了一种快速辨识机床主轴热态特性的方法，包括：在机床主轴上布置n个温度测量点，在采样间隔Δt下得到各个时间下n个温度测量点的温度值，采样次数为m，得到温度采样矩阵；设定延时时间，求出热特征值b_r(r=1,2,…,p)；选择温度测试点中温度最高的点为关键点，该关键点的温升曲线；计算该关键点的温度实测值与估计值的均方根误差σ；得到不同采样次数m下的均方根误差σ，直至均方根误差σ存在极小均方根值，即在极小均方根值对应的关键点的温升曲线即为能准确反映机床主轴热态特性的曲线。本发明方法只需30min左右即可实现快速辨识热态特性的目的，获取关键点的温升曲线，缩短辨识时间，结果准确，操作简单。

Description

一种快速辨识机床主轴热态特性的方法

技术领域

本发明涉及机床热态特性辨识方法领域，尤其涉及一种快速辨识机床主轴热态特性的方法。

背景技术

随着机床不断朝高速、高精度方向发展，热误差对加工精度的影响越来越大，甚至热误差所引起的加工误差占总误差的75%左右。如何有效减小热误差则成为研究重点。目前，主要有三方面的做法：一是减小温度变化，可以通过加强冷却、控制环境温度、加工前的预热处理等方式进行；二是减小机床结构对温度变化的敏感度，实施方法有采用对称结构设计、使用低热膨胀系数材料、隔离热源；三是补偿热误差。而无论采用上述哪种方式减小热误差，均首先需要了解机床热态特性，尤其是主轴部件的热态特性。因此，辨识热态特性便是设计和使用机床的关键因素。

目前，辨识机床热态特性方法有两种方式。一种是通过有限元分析进行数值仿真，模拟出机床各节点的温升曲线，计算热平衡时间及稳态温度；另一种是通过对机床进行实验，来测量出机床热特性。由于机床结构的复杂性、热源位置与强度的未知性、结合面热阻的不确定性等因素的影响，在进行有限元仿真时，热边界条件往往难以准确估计，使得最终的热态特性不符合实际情况。因此，直接对机床测量来辨识热态特性就成为可靠有效的手段。

但是现在的测量是通过机床空转，等待机床达到热平衡状态，然后对测量数据处理。这种方式一般需要3到6个小时才能完成测量，因此，费时费力。

申请公开号为CN 101972948A的中国发明专利申请公开了一种模拟工况载荷条件下机床热误差实验装置，包括主轴热误差测试试验系统，该实验装置还包括模拟工况载荷主轴加载装置，所述主轴加载装置包括力矩载荷加载装置和径向载荷加载装置，所述力矩载荷加载装置包括磁粉制动器，所述磁粉制动器安装在加载装置支架上，所述磁粉制动器通过导向键和检验棒连接，所述检验棒与所述主轴连接，所述磁粉制动器通过所述导向键和所述检验棒将模拟力矩载荷施加到所述主轴上；所述径向力载荷加载装置包括拉压力计，所述拉压力计安装在径向力加载装置支承体上，所述径向力加载装置与所述加载装置支架固定连接，所述拉压力计设有受拉端和受压端，所述拉压力计的受拉端设有径向力调整螺栓，所述拉压力计受压端设有径向力导杆，所述径向力导杆的远离拉压力计端安装有滚动轴承，所述滚动轴承与所述检验棒滚动连接，所述径向力调整螺栓通过挤压所述拉压力计产生模拟径向力载荷并通过所述径向力导杆、所述滚动轴承和所述检验棒滚动接触将模拟径向力载荷施加到所述主轴上。虽然该模拟工况载荷条件下机床热误差实验装置，使得主轴的热误差非常贴近于实际工况，但是其对机床主轴热特性的测试方法仍没有改进，仍存在费时费力的问题。

发明内容

本发明提供了一种快速辨识机床主轴热态特性的方法，通过对主轴上各温度测量点上温度数据的处理，实现快速辨识热态特性的目的，获取关键点的温升曲线，缩短辨识时间，结果准确，操作简单。

一种快速辨识机床主轴热态特性的方法，包括以下步骤：

（一）在机床主轴上布置n个温度测量点，在采样间隔Δt下得到各个时间下n个温度测量点的温度值，采样次数为m，得到温度采样矩阵为[{T(t₁)},{T(t₂)},……,{T(t_m)}]，具体为：

其中，m为采样次数，Δt为间隔时间，1,2,…,n代表不同的各温度测量点；

（二）设定延时时间τ，τ=gΔt，g为正整数，s=m-2τ；

构造两个温度采样序列：

[Z]_n×s=[{{T(t₁)}-{T(t₁+τ)}},{{T(t₂)}-{T(t₂+τ)}},…,{{T(t_s)}-{T(t_s+τ)}}]_n×s    (3)；

${\left[\hat{Z}\right]}_{n\×s}={\left[\right\{\left\{T(t_1+\mathrm{\τ})\right\}-\left\{T(t_1+2\mathrm{\τ})\right\}\},\{\left\{T(t_2+\mathrm{\τ})\right\}-\left\{T(t_2+2\mathrm{\τ})\right\}\},...,\{\left\{T(t_s+\mathrm{\τ})\right\}-\left\{T(t_s+2\mathrm{\τ})\right\}\left\}\right]}_{n\×s}---\left(4\right);$

矩阵[H]_n×n满足如下式(8)，

${\left[H\right]}_{n\×n}{\left[Z\right]}_{n\×s}={\left[\hat{Z}\right]}_{n\×s}---\left(8\right);$

则，将式(3)和(4)代入式(8)，求解得到[H]_n×n，

根据式(10)可知，

矩阵[H]_n×n的特征值矩阵为特征向量矩阵为对矩阵[H]_n×n进行特征值分解，便可求出热特征值b_r(r=1,2,…,p)，舍去热特征值b_r为虚数及负数的b_r，按保留的热特征值b_r个数重新确定p，p=保留的热特征值b_r个数，并将保留的热特征值b_r按热特征值b_r的大小由小到大重新排列；

（三）选择温度测试点中温度最高的点为关键点，采样次数为m，该关键点的温度列向量为 ${\left[Q\right]}_{m\×1}=\left[\begin{array}{c}T\left(t_1\right)\\ T\left(t_2\right)\\ .\\ .\\ .\\ T\left(t_m\right)\end{array}\right],$ 其中，T(t₁)为第1次采样得到的该关键点的温度，t₁为第1次采样对应的采样时间，依次类推，T(t_m)为第m次采样得到的该关键点的温度，t_m为第m次采样对应的采样时间；

根据式(15)，

[P]_m×(p+1)X_(p+1)×1=[Q]_m×1    (15)；

其中，

${\left[P\right]}_{m\×(p+1)}=\left[\begin{array}{cccc}1& e^{-b_1{t}_1}& ...& e^{-b_p{t}_1}\\ 1& e^{-b_1{t}_2}& ...& e^{-b_p{t}_2}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& e^{-b_1{t}_m}& ...& e^{-b_p{t}_m}\end{array}\right],$ $X_{(p+1)\×1}=\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ .\\ .\\ .\\ {\mathrm{\γ}}_p\end{array}\right],$ ${\left[Q\right]}_{m\×1}=\left[\begin{array}{c}T\left(t_1\right)\\ T\left(t_2\right)\\ .\\ .\\ .\\ T\left(t_m\right)\end{array}\right],$

计算得到 $\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ .\\ .\\ .\\ {\mathrm{\γ}}_p\end{array}\right];$

该关键点的温升曲线表达式如式(14)所示：

$T\left(t_k\right)={\mathrm{\γ}}_0+\underset{r=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\γ}}_r{e}^{-b_r{t}_k}=\left[\begin{array}{cccc}1& e^{-b_1{t}_k}& ...& e^{-b_p{t}_k}\end{array}\right]\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ .\\ .\\ .\\ {\mathrm{\γ}}_p\end{array}\right]---\left(14\right);$

其中，T(t_k)表示t_k时刻的温度，t_k为第k次采样的时间，γ_r（r=0,1,2,…,p）表示温升曲线表达式系数，b_r(r=1,2,…,p)表示热特征值；

（四）计算该关键点的温度实测值与估计值的均方根误差σ，均方根误差表达式(16)如为

$\mathrm{\σ}=\sqrt{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{[T_e\left(k\right)-T_o\left(k\right)]}^2/(m-1)}---\left(16\right);$

其中，m为采样次数，T_e(k)为在第k次采样时刻估计的温度值，T_e(k)通过该关键点的温升曲线得到，T_o(k)为第k次采样时刻测量的温度值；

（五）在某一采样次数m下得到均方根误差σ，采样次数m依次增大，返回步骤（一），得到不同采样次数m下的均方根误差σ，直至均方根误差σ出现先减小后增大的情形，即均方根误差σ存在极小均方根值，即在极小均方根值对应的采样次数m所对应的时间为辨识时间，该辨识时间下关键点的温升曲线即为能准确反映机床主轴热态特性的曲线。

通过上述方法，得到的辨识时间较短，一般在1小时以内，30min左右，从而通过在较短的辨识时间下，就可以预判机床主轴在之后时间所反映的热特性，从而快速辨识机床主轴热态特性。

步骤（一）中，温度采样矩阵为：

其中，m为采样次数，Δt为间隔时间，1,2,…,n代表不同的各温度测量点；

温度采样矩阵{T(t)}中的每一列即为同一采样时间不同温度测量点下的温度值，{T(t₁)}为温度采样矩阵{T(t)}第一列的温度数据，即第1次采样时间下各温度测量点的温度值；{T(t₂)}为温度采样矩阵{T(t)}第二列的温度数据，第2次采样时间下各温度测量点的温度值；{T(t_m)}为温度采样矩阵{T(t)}第m列的温度数据，第m次采样时间下各温度测量点的温度值；

如{T(t₁)}为：

如{T(t_m)}为：

温度采样矩阵{T(t)}也可写成{T(t)}=[{T(t₁)},{T(t₂)},……,{T(t_m)}]；

作为优选，所述各温度测量点沿机床主轴的轴向布置，从而使得机床主轴的热特性得到更完整的表达；

步骤（二）中，设定延时时间τ，τ=gΔt，g为正整数，g可以根据需要从正整数中取一个，m=s+2τ，s=m-2τ，s为带入计算的次数，；

现有技术中，机床在空转时的自由温升表达式{T(t)}为：

式中为p维列向量，b_r(r=1,2,…,p)为热特征值；

构造如下式(1)和式(2)两式；

构造两个温度采样序列

[Z]_n×s=[{{T(t₁)}-{T(t₁+τ)}},{{T(t₂)}-{T(t₂+τ)}},…,{{T(t_s)}-{T(t_s+τ)}}]_n×s    (3)；

${\left[\hat{Z}\right]}_{n\×s}={\left[\right\{\left\{T(t_1+\mathrm{\τ})\right\}-\left\{T(t_1+2\mathrm{\τ})\right\}\},\{\left\{T(t_2+\mathrm{\τ})\right\}-\left\{T(t_2+2\mathrm{\τ})\right\}\},...,\{\left\{T(t_s+\mathrm{\τ})\right\}-\left\{T(t_s+2\mathrm{\τ})\right\}\left\}\right]}_{n\×s}---\left(4\right);$

作为优选，采样次数m要大于温度测量点的个数n，能够提高计算的准确性；

将式(1)、式(2)代入式(3)、式(4)，得到式(5)、式(6)；

其中，式(5)和式(6)中[E]_p×s如式(7)所示，

${\left[E\right]}_{p\×s}=\left[\begin{array}{cccc}{e}^{-b_1{t}_1}& e^{-b_1{t}_2}& ...& e^{-b_1{t}_s}\\ e^{-b_2{t}_1}& e^{-b_2{t}_2}& ...& e^{-b_2{t}_s}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e^{-b_p{t}_1}& e^{-b_p{t}_2}& ...& e^{-b_p{t}_s}\end{array}\right]---\left(7\right);$

假设存在一个矩阵[H]_n×n满足如下式(8)，

${\left[H\right]}_{n\×n}{\left[Z\right]}_{n\×s}={\left[\hat{Z}\right]}_{n\×s}---\left(8\right);$

将式（5）、（6）代入式（8），得到下式(9)，

一般地，是非奇异矩阵，式（9）可变为式(10)，

根据线性代数方阵特征值及特征向量定义，设A是n阶矩阵，即方阵A，如果数λ和n维非零列向量x使关系式Ax=λx成立，那么，数λ成为方阵A的特征值，非零向量x成为A的对应于特征值λ的特征向量，将A所有的特征值和特征向量合在一起，即为

AX=XD    (11)；

其中，特征向量矩阵X=[x₁ x₂ … x_n]_n×n，特征值矩阵

式（11）也可写成

X^-1AX=D    (12)；

对比式（10）和式（12），可以看出矩阵[H]_n×n的特征值矩阵为特征向量矩阵为因此，只需对矩阵[H]_n×n进行特征值分解，便可求出热特征值br(r=1,2,…,p)，舍去热特征值b_r为虚数及负数的b_r，按保留的热特征值b_r个数重新确定p，p=保留的热特征值b_r个数，并将保留的热特征值b_r按热特征值b_r的大小由小到大重新排列。

将式（8）转置，得

${\left[Z\right]}_{s\×n}^T{\left[H\right]}_{n\×n}^T={\left[\hat{Z}\right]}_{s\×n}^T---\left(13\right);$

对矩阵[H]_n×n进行特征值分解，求出热特征值b_r(r=1,2,…,p)，在舍去热特征值b_r为虚数及负数的b_r，按保留的热特征值b_r个数重新确定p，p=保留的热特征值b_r个数，并将保留的热特征值b_r按热特征值b_r的大小由小到大重新排列。

步骤（三）中，计算关键点的温升表达式中的{γ_r}列向量，该关键点的温升曲线表达式如式(14)所示：

$T\left(t_k\right)={\mathrm{\γ}}_0+\underset{r=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\γ}}_r{e}^{-b_r{t}_k}=\left[\begin{array}{cccc}1& e^{-b_1{t}_k}& ...& e^{-b_p{t}_k}\end{array}\right]\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ .\\ .\\ .\\ {\mathrm{\γ}}_p\end{array}\right]---\left(14\right)$

其中，T(t_k)表示t_k时刻的温度，γ_r（r=0,1,2,…,p）表示温升曲线表达式系数，b_r(r=1,2,…,p)表示热特征值，{γ_r}中的p与b_r中的p具有相同含义。

进行m次采样后，得到如式(15)所示：

[P]_m×(p+1)X_(p+1)×1=[Q]_m×1    (15)

其中，

${\left[P\right]}_{m\×(p+1)}=\left[\begin{array}{cccc}1& e^{-b_1{t}_1}& ...& e^{-b_p{t}_1}\\ 1& e^{-b_1{t}_2}& ...& e^{-b_p{t}_2}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& e^{-b_1{t}_m}& ...& e^{-b_p{t}_m}\end{array}\right],$ $X_{(p+1)\×1}=\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ .\\ .\\ .\\ {\mathrm{\γ}}_p\end{array}\right],$ ${\left[Q\right]}_{m\×1}=\left[\begin{array}{c}T\left(t_1\right)\\ T\left(t_2\right)\\ .\\ .\\ .\\ T\left(t_m\right)\end{array}\right],$

求解式（15）的X_(p+1)×1，即可得到一点的γ_r（r=0,1,2,…,p），那么该关键点的温升表达式即可求得；

步骤（四）中，计算温度实测值与估计值的均方根误差σ，均方根误差表达式(16)为

$\mathrm{\σ}=\sqrt{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{[T_e\left(k\right)-T_o\left(k\right)]}^2/(m-1)}---\left(16\right)$

式(16)中m为采样次数，T_e(k)为在第k次采样时刻估计的温度值，T_o(k)为第k次采样时刻测量的温度值。

步骤（五）中，在某一采样次数m（即某一采样时间）下，代入式(16)可得到一均方根误差σ，随着采样次数m依次增加，即随着采样时间的增加，返回步骤（一），经过步骤（一）至（四）后，可再得到一均方根误差，即随着采样次数m的增加，得到不同的采样次数m下对应的不同的均方根误差，也即不同的采样时间下对应的不同的均方根误差，均方根误差随着采样次数m（也即采样时间）的增加，呈现一定的变化，当采样次数m较小时，得到的温升曲线与实际温度值拟合度差，均方根误差也大，均方根误差σ会随着采样次数m的增加出现先减小后增大的情形，即均方根误差σ存在极小均方根值，即在极小均方根值对应的采样次数m所对应的时间为辨识时间，该辨析时间下关键点的温升曲线即为能准确反映机床主轴热态特性的曲线。

热平衡时间定义为当温度达到稳态温度的95%时所处的时间。

与现有技术相比，本发明具有如下优点：

本发明快速辨识机床主轴热态特性的方法，通过对主轴上各温度测量点上温度数据的处理，实现快速辨识特态特性的目的，获取关键点的温升曲线，大幅减少实际辨识热特性操作的时间，缩短生产准备的时间，且此发明方法只需温度采样数据，无需热激励即可达到辨识目的，结果准确，操作简单

通过上述方法，得到的辨识时间较短，一般在1小时以内，只需30min左右即可到达辨识目的，在较短的辨识时间下，就可以预判机床主轴在之后时间所反映的热特性，从而快速辨识机床主轴热态特性。

附图说明

图1为本发明快速辨识机床主轴热态特性的方法流程图；

图2为本发明中机床主轴温度传感器布置的示意图；

图3为本发明辨识温升曲线与实测温升曲线对比图，其中，a为26.35min辨识时间的关键点的温升曲线，b为该关键点的实测温度曲线。

具体实施方式

下面结合附图和实施例对本发明方法做进一步详细说明，以下实施例不构成对本发明的限定。

如图1所示，一种快速辨识机床主轴热态特性的方法，包括以下步骤：

（一）在机床主轴上布置n个温度测量点，在采样间隔Δt下得到各个时间下n个温度测量点的温度值，采样次数为m，得到温度采样矩阵为[{T(t₁)},{T(t₂)},……,{T(t_m)}]，具体为：

其中，m为采样次数，Δt为间隔时间，1,2,…,n代表不同的各温度测量点；

（二）设定延时时间τ，τ=gΔt，g为正整数，s=m-2τ，s为带入计算的次数；

构造两个温度采样序列：

[Z]_n×s=[{{T(t₁)}-{T(t₁+τ)}},{{T(t₂)}-{T(t₂+τ)}},…,{{T(t_s)}-{T(t_s+τ)}}]_n×s    (3)；

${\left[\hat{Z}\right]}_{n\×s}={\left[\right\{\left\{T(t_1+\mathrm{\τ})\right\}-\left\{T(t_1+2\mathrm{\τ})\right\}\},\{\left\{T(t_2+\mathrm{\τ})\right\}-\left\{T(t_2+2\mathrm{\τ})\right\}\},...,\{\left\{T(t_s+\mathrm{\τ})\right\}-\left\{T(t_s+2\mathrm{\τ})\right\}\left\}\right]}_{n\×s}---\left(4\right);$

矩阵[H]_n×n满足如下式子(8)，

${\left[H\right]}_{n\×n}{\left[Z\right]}_{n\×s}={\left[\hat{Z}\right]}_{n\×s}---\left(8\right);$

则，将式(3)和(4)代入式(8)，求解得到[H]_n×n，

根据式(10)可知，

矩阵[H]_n×n的特征值矩阵为特征向量矩阵为对矩阵[H]_n×n进行特征值分解，便可求出热特征值b_r(r=1,2,…,p)，舍去热特征值b_r为虚数及负数的b_r，按保留的热特征值b_r个数重新确定p，p=保留的热特征值b_r个数，并将保留的热特征值b_r按热特征值b_r的大小由小到大重新排列；

（三）选择温度测试点中温度最高的点为关键点，采样次数为m，该关键点的温度列向量为 ${\left[Q\right]}_{m\×1}=\left[\begin{array}{c}T\left(t_1\right)\\ T\left(t_2\right)\\ .\\ .\\ .\\ T\left(t_m\right)\end{array}\right],$ 其中，T(t₁)为第1次采样得到的该关键点的温度，t₁为第1次采样对应的采样时间，依次类推，T(t_m)为第m次采样得到的该关键点的温度，t_m为第m次采样对应的采样时间；

根据式(15)，

[P]_m×(p+1)X_(p+1)×1=[Q]_m×1    (15)；

其中，

${\left[P\right]}_{m\×(p+1)}=\left[\begin{array}{cccc}1& e^{-b_1{t}_1}& ...& e^{-b_p{t}_1}\\ 1& e^{-b_1{t}_2}& ...& e^{-b_p{t}_2}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& e^{-b_1{t}_m}& ...& e^{-b_p{t}_m}\end{array}\right],$ $X_{(p+1)\×1}=\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ .\\ .\\ .\\ {\mathrm{\γ}}_p\end{array}\right],$ ${\left[Q\right]}_{m\×1}=\left[\begin{array}{c}T\left(t_1\right)\\ T\left(t_2\right)\\ .\\ .\\ .\\ T\left(t_m\right)\end{array}\right];$

计算得到 $\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ .\\ .\\ .\\ {\mathrm{\γ}}_p\end{array}\right];$

该关键点的温升曲线表达式如式(14)所示：

$T\left(t_k\right)={\mathrm{\γ}}_0+\underset{r=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\γ}}_r{e}^{-b_r{t}_k}=\left[\begin{array}{cccc}1& e^{-b_1{t}_k}& ...& e^{-b_p{t}_k}\end{array}\right]\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ .\\ .\\ .\\ {\mathrm{\γ}}_p\end{array}\right]---\left(14\right);$

其中，T(t_k)表示t_k时刻的温度，t_k为第k次采样的时间，γ_r（r=0,1,2,…,p）表示温升曲线表达式系数，b_r(r=1,2,…,p)表示热特征值；

（四）计算该关键点的温度实测值与估计值的均方根误差σ，均方根误差表达式(16)如为

$\mathrm{\σ}=\sqrt{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{[T_e\left(k\right)-T_o\left(k\right)]}^2/(m-1)}---\left(16\right);$

其中，m为采样次数，T_e(k)为在第k次采样时刻估计的温度值，T_e(k)通过该关键点的温升曲线得到，T_o(k)为第k次采样时刻测量的温度值；

（五）在某一采样次数m下得到均方根误差σ，采样次数m依次增大，返回步骤（一），得到不同采样次数m下的均方根误差σ，直至均方根误差σ出现先减小后增大的情形，即均方根误差σ存在极小均方根值，即在极小均方根值对应的采样次数m所对应的时间为辨识时间，该辨识时间下关键点的温升曲线即为能准确反映机床主轴热态特性的曲线。

实施例1

如图2所示，对一立式加工中心机床主轴做热特性辨识，在沿机床主轴的轴向布置温度传感器1、2、3、4、5、6、7和8，形成不同的温度测量点，即设置8个温度测量点，其中，9为心轴，温度传感器1位于轴承座10前端面，温度传感器2、3、4位于轴承座10外圆周面，温度传感器5位于法兰11的外圆周面，温度传感器6、7、8位于主轴箱12外表面。温度传感器3位于轴承座10外圆周面上（主轴前轴承处），温升最高，且对刀尖点热变形误差影响较大，温度传感器3所对应的温度测量点作为关键点，衡量机床主轴热态特性。

机床空载，主轴转速为5000r/min，环境温度为10.4℃，采样间隔Δt为1s，延时时间τ为8s，采用具体实施方式中的方法计算不同采样次数m（也即不同采样时间）下的均方根误差σ，其部分结果（关键部分）如表1所示。

表1

1、在采样时间5min下，采样次数m=301下计算均方根误差σ：

（一）如图2在机床主轴上布置8个温度测量点，在采样间隔1s下得到各个时间下8个温度测量点的温度值，进行5min采样，采样次数m=301，得到温度采样矩阵为[{T(t₁)},{T(t₂)},……,{T(t₃₀₁)}]，具体为：

$\left\{\begin{array}{cccccc}10.45& 10.45& 10.46& 10.46& ......& 11.10\\ 10.16& 10.17& 10.17& 10.17& ......& 10.81\\ 10.49& 10.50& 10.51& 10.51& ......& 11.23\\ 10.31& 10.32& 10.32& 10.33& ......& 10.71\\ 10.24& 10.24& 10.24& 10.25& ......& 10.42\\ 10.37& 10.37& 10.37& 10.38& ......& 10.43\\ 10.41& 10.42& 10.42& 10.43& ......& 10.65\\ 10.35& 10.36& 10.37& 10.37& ......& 10.74\end{array}\right\}$

（二）设定延时时间τ=8s，τ=gΔt，g=8，s为带入计算的次数，s=m-2τ，s=285，n为温度传感器的数目，n=8；

构造两个温度采样序列：

[Z]_n×s=[{{T(t₁)}-{T(t₁+τ)}},{{T(t₂)}-{T(t₂+τ)}},…,{{T(t_s)}-{T(t_s+τ)}}]_n×s    (3)；

${\left[\hat{Z}\right]}_{n\×s}={\left[\right\{\left\{T(t_1+\mathrm{\τ})\right\}-\left\{T(t_1+2\mathrm{\τ})\right\}\},\{\left\{T(t_2+\mathrm{\τ})\right\}-\left\{T(t_2+2\mathrm{\τ})\right\}\},...,\{\left\{T(t_s+\mathrm{\τ})\right\}-\left\{T(t_s+2\mathrm{\τ})\right\}\left\}\right]}_{n\×s}---\left(4\right);$

${\left[Z\right]}_{8\×285}={\left[\begin{array}{cccc}-0.01& -0.01& ...& -0.15\\ -0.01& -0.01& ...& -0.08\\ -0.02& -0.02& ...& -0.19\\ -0.02& -0.01& ...& -0.06\\ -0.01& -0.01& ...& -0.08\\ -0.01& -0.02& ...& -0.12\\ -0.02& -0.01& ...& -0.11\\ -0.02& -0.02& ...& -0.05\end{array}\right]}_{8\×285}---\left(3\right);$

${\left[\hat{Z}\right]}_{8\×285}={\left[\begin{array}{cccc}-0.01& -0.01& ...& -0.12\\ -0.01& -0.01& ...& -0.1\\ -0.02& -0.02& ...& -0.17\\ -0.02& -0.02& ...& -0.08\\ -0.02& -0.01& ...& -0.05\\ -0.01& -0.02& ...& -0.1\\ -0.02& -0.01& ...& -0.11\\ -0.02& -0.02& ...& -0.12\end{array}\right]}_{8\×285}---\left(4\right);$

矩阵[H]_8×8满足如下式子(8)，

${\left[H\right]}_{8\×8}{\left[Z\right]}_{8\×285}={\left[\hat{Z}\right]}_{8\×285}---\left(8\right);$

则，将式(3)和(4)代入式(8)，求解得到[H]_8×8，

${\left[H\right]}_{8\×8}=\left[\begin{array}{cccccccc}0.0143& 0.6921& 0.3484& 0.0448& -0.1842& 0.0347& -0.1066& -0.0689\\ 0.4049& 0.7790& 0.3217& -0.4174& -0.2996& 0.2874& 0.1633& -0.2750\\ 0.4378& 0.6830& -0.1219& -0.0810& -0.0595& 0.0458& 0.0965& 0.0228\\ 0.2335& 0.0082& 0.2957& -0.1615& 0.2516& 0.0718& -0.0322& 0.3162\\ -0.0834& -0.1002& 0.1194& 0.2502& -0.1206& 0.2525& 0.3460& 0.3316\\ -0.1509& 0.0824& -0.0510& -0.1001& 0.3229& -0.0630& 0.4926& 0.3012\\ -0.1618& 0.2259& -0.1089& -0.0433& 0.2049& 0.5079& 0.0594& 0.1792\\ 0.1166& 0.1226& 0.0715& 0.1959& 0.2833& 0.1615& 0.4203& -0.1765\end{array}\right]$

根据式(10)可知，矩阵[H]_8×8的特征值矩阵为特征向量矩阵为对[H]_8×8进行特征值分解，得到

$\left[\begin{array}{cccccccc}1.0277+0.0275i& 0& 0& 0& 0& 0& 0& 0\\ 0& 1.0277-0.0275i& 0& 0& 0& 0& 0& 0\\ 0& 0& 0.5726& 0& 0& 0& 0& 0\\ 0& 0& 0& -0.5142+0.147i& 0& 0& 0& 0\\ 0& 0& 0& 0& -5142-0.147i& 0& 0& 0\\ 0& 0& 0& 0& 0& -0.4527+0.0327i& 0& 0\\ 0& 0& 0& 0& 0& 0& -0.4527-0.0327i& 0\\ 0& 0& 0& 0& 0& 0& 0& -0.4851\end{array}\right];$

则b₁=-0.2076-0.2007i，i为虚数单位；

则b₂=-0.2076+0.2007i；

则b₃=4.1818；

则b₄=4.6940-21.4735i；

则b₅=4.6940+21.4735i；

则b₆=5.9244-23.0211i；

则b₇=5.9244+23.0211i；

则b₈=5.4255-23.5619i；

由于b_r是正实数，因此，需要舍去虚数及负数，得到最终特征值，按b_r个数确定p，p=1，并按升序重新排列，b₁=4.1818min^-1。

（三）选择温度测试点中温度最高的点为关键点，温度传感器3所对应的温度测量点作为关键点，采样次数为m，该关键点的温度列项量为 ${\left[Q\right]}_{m\×1}=\left[\begin{array}{c}T\left(t_1\right)\\ T\left(t_2\right)\\ .\\ .\\ .\\ T\left(t_m\right)\end{array}\right]=\left[\begin{array}{c}10.49\\ 10.50\\ .\\ .\\ .\\ 11.23\end{array}\right],$ 即将温度采样矩阵中第三行的数据转化为现在的温度列项量数据，其中，T(t₁)为第1次采样得到的该关键点的温度，t₁为第1次采样对应的采样时间，依次类推，T(t_m)为第m次采样得到的该关键点的温度，t_m为第m次采样对应的采样时间；

根据式(15)，[P]_m×(p+1)X_(p+1)×1=[Q]_m×1    (15)；

其中，m=301，p=1，

${\left[P\right]}_{301\×2}=\left[\begin{array}{cc}1& e^{-4.1818\×t_1}\\ 1& e^{-4.1818\×t_2}\\ .& .\\ .& .\\ .& .\\ 1& e^{-4.1818\×t_m}\end{array}\right],$ $X_{2\×1}=\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\end{array}\right],$ ${\left[Q\right]}_{301\×1}=\left[\begin{array}{c}10.49\\ 10.50\\ .\\ .\\ .\\ 11.23\end{array}\right],$

可以计算得到 $\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\end{array}\right]=\left[\begin{array}{c}10.8409\\ -0.4770\end{array}\right];$

该关键点的温升曲线表达式如式(14)所示：

$\left(t_k\right)={\mathrm{\γ}}_0+{\mathrm{\γ}}_1{e}^{-b_1{t}_k}=10.8409-0.4770\×e^{-4.1818\×t_k}---\left(14\right)$

其中，T(t_k)表示t_k时刻的温度，t_k为第k次采样的时间，单位min，γ_r（r=0,1,2,…,p）表示温升曲线表达式系数，b_r(r=1,2,…,p)表示热特征值，p=1；

（四）计算该关键点的温度实测值与估计值的均方根误差σ，均方根误差表达式(16)如：

$\mathrm{\σ}=\sqrt{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{[T_e\left(k\right)-T_o\left(k\right)]}^2/(m-1)}---\left(16\right);$

其中，m为采样次数，为301，T_e(k)为在第k次采样时刻估计的温度值，T_e(k)通过该关键点的温升曲线得到，T_o(k)为第k次采样时刻测量的温度值，计算得到σ=0.7424℃。

2、在采样时间26.35min下，采样次数m=1582下计算均方根误差σ：

（一）如图2在机床主轴上布置8个温度测量点，在采样间隔1s下得到各个时间下8个温度测量点的温度值，进行26.35min采样，采样次数m=1582，得到温度采样矩阵为[{T(t₁)},{T(t₂)},……,{T(t₁₅₈₂)}]，具体为：

$\left\{\begin{array}{cccccc}10.45& 10.45& 10.46& 10.46& ......& 12.80\\ 10.16& 10.17& 10.17& 10.17& ......& 12.81\\ 10.49& 10.50& 10.51& 10.51& ......& 13.25\\ 10.31& 10.32& 10.32& 10.33& ......& 12.61\\ 10.24& 10.24& 10.24& 10.25& ......& 12.25\\ 10.37& 10.37& 10.37& 10.38& ......& 12.04\\ 10.41& 10.42& 10.42& 10.43& ......& 12.14\\ 10.35& 10.36& 10.37& 10.37& ......& 12.65\end{array}\right\};$

（二）设定延时时间τ=8s，τ=gΔt，g=8，s为带入计算的次数，s=m-2τ，s=1566，n为温度传感器的数目，n=8；

构造两个温度采样序列：

[Z]_n×s=[{{T(t₁)}-{T(t₁+τ)}},{{T(t₂)}-{T(t₂+τ)}},…,{{T(t_s)}-{T(t_s+τ)}}]_n×s    (3)；

${\left[\hat{Z}\right]}_{n\×s}={\left[\right\{\left\{T(t_1+\mathrm{\τ})\right\}-\left\{T(t_1+2\mathrm{\τ})\right\}\},\{\left\{T(t_2+\mathrm{\τ})\right\}-\left\{T(t_2+2\mathrm{\τ})\right\}\},...,\{\left\{T(t_s+\mathrm{\τ})\right\}-\left\{T(t_s+2\mathrm{\τ})\right\}\left\}\right]}_{n\×s}---\left(4\right);$

${\left[Z\right]}_{8\×1566}={\left[\begin{array}{cccc}-0.01& -0.01& ...& -0.12\\ -0.01& -0.01& ...& -0.08\\ -0.02& -0.02& ...& -0.21\\ -0.02& -0.01& ...& -0.08\\ -0.01& -0.01& ...& -0.06\\ -0.01& -0.02& ...& -0.15\\ -0.02& -0.01& ...& -0.14\\ -0.02& -0.02& ...& -0.10\end{array}\right]}_{8\×1566}---\left(3\right);$

${\left[\hat{Z}\right]}_{8\×1566}={\left[\begin{array}{cccc}-0.01& -0.01& ...& -0.12\\ -0.01& -0.01& ...& -0.09\\ -0.02& -0.02& ...& -0.15\\ -0.02& -0.02& ...& -0.09\\ -0.02& -0.01& ...& -0.07\\ -0.01& -0.02& ...& -0.12\\ -0.02& -0.01& ...& -0.10\\ -0.02& -0.02& ...& -0.15\end{array}\right]}_{8\×1566}---\left(4\right);$

矩阵[H]_8×8满足如下式(8)，

${\left[H\right]}_{8\×8}{\left[Z\right]}_{8\×1566}={\left[\hat{Z}\right]}_{8\×1566}---\left(8\right);$

则，将式(3)和(4)代入式(8)，求解得到[H]_8×8，

${\left[H\right]}_{8\×8}=\left[\begin{array}{cccccccc}0.8515& 0.1394& 0.0534& -0.1132& 0.2779& -0.1586& -0.0317& -0.0591\\ 0.1110& 0.9226& 0.0136& 0.0297& -0.2002& 0.1425& 0.0446& -0.0338\\ 0.1205& 0.0977& 0.8820& -0.1700& 0.0122& 0.1291& -0.0322& -0.0088\\ 0.0343& 0.1800& -0.0735& 0.6899& 0.2499& -0.0546& 0.0091& -0.0260\\ 0.2439& -0.2201& -0.0210& 0.1812& 0.3912& 0.3868& 0.0580& 0.0538\\ -0.1344& 0.1295& 0.0194& -0.1077& 0.3327& 0.7534& -0.0052& -0.0265\\ 0.0336& 0.0398& -0.0440& -0.0532& -0.0046& 0.0540& 0.9650& 0.0182\\ 0.0285& 0.0844& -0.0387& -0.1313& 0.0994& -0.0294& 0.0231& 0.9720\end{array}\right]$

根据式(10)可知，矩阵[H]_8×8的特征值矩阵为特征向量矩阵为对[H]_8×8进行特征值分解，得到

$\left[\begin{array}{cccccccc}-0.1060& 0& 0& 0& 0& 0& 0& 0\\ 0& 0.6926& 0& 0& 0& 0& 0& 0\\ 0& 0& 0.9186& 0& 0& 0& 0& 0\\ 0& 0& 0& 0.9470& 0& 0& 0& 0\\ 0& 0& 0& 0& 0.9887& 0& 0& 0\\ 0& 0& 0& 0& 0& 0.9884& 0& 0\\ 0& 0& 0& 0& 0& 0& 0.9941+0.0114i& 0\\ 0& 0& 0& 0& 0& 0& 0& 0.9941-0.0114i\end{array}\right];$

则b₁=16.8322-23.5619i；

则b₂=2.7548；

则b₃=0.6376；

则b₄=0.4084；

则b₅=0.0852；

则b₆=0.0120；

则b₇=0.0439-0.0860i；

则b₈=0.0439+0.0860i；

舍去为虚数及负数的b_r，按b_r个数确定p，p=5，并按升序（即b_r热特征值的大小由小到大）重新排列，得到b₁=0.0120min^-1，b₂=0.0852min^-1，b₃=0.4084min^-1，b₄=0.6376min^-1，b₅=2.7548min^-1；

（三）选择温度测试点中温度最高的点为关键点，温度传感器3所对应的温度测量点作为关键点，采样次数为m，该关键点的温度列项量为 ${\left[Q\right]}_{m\×1}=\left[\begin{array}{c}T\left(t_1\right)\\ T\left(t_2\right)\\ .\\ .\\ .\\ T\left(t_m\right)\end{array}\right]=\left[\begin{array}{c}10.49\\ 10.50\\ .\\ .\\ .\\ 13.25\end{array}\right],$ 即将温度采样矩阵中第三行的数据转化为现在的温度列项量数据，其中，T(t₁)为第1次采样得到的该关键点的温度，t₁为第1次采样对应的采样时间，依次类推，T(t_m)为第m次采样得到的该关键点的温度，t_m为第m次采样对应的采样时间；

根据式(15)，

[P]_m×(p+1)X_(p+1)×1=[Q]_m×1    (15)；

其中，m=1582，p=5

${\left[P\right]}_{1582\×6}=\left[\begin{array}{cccccc}1& e^{-0.0120\×t_1}& e^{-0.0852\×t_1}& e^{-0.4084\×t_1}& e^{-0.6376\×t_1}& e^{-2.7548\×t_1}\\ 1& e^{-0.0120\×t_2}& e^{-0.0852\×t_2}& e^{-0.4084\×t_2}& e^{-0.6376\×t_2}& e^{-2.7548\×t_2}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 1& e^{-0.0120\×t_m}& e^{-0.0852\×t_m}& e^{-0.4084\×t_m}& e^{-0.6376\×t_m}& e^{-2.7548\×t_m}\end{array}\right],$

$X_{6\×1}=\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ {\mathrm{\γ}}_2\\ {\mathrm{\γ}}_3\\ {\mathrm{\γ}}_4\\ {\mathrm{\γ}}_5\end{array}\right],$ ${\left[Q\right]}_{1582\×1}=\left[\begin{array}{c}10.49\\ 10.50\\ .\\ .\\ .\\ 13.25\end{array}\right];$

可以计算得到 $\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ {\mathrm{\γ}}_2\\ {\mathrm{\γ}}_3\\ {\mathrm{\γ}}_4\\ {\mathrm{\γ}}_5\end{array}\right]=\left[\begin{array}{c}19.9723\\ -9.3344\\ 0.2691\\ -2.3818\\ 2.4242\\ -0.4552\end{array}\right];$

其温升表达式为：

$T\left(t_k\right)={\mathrm{\γ}}_0+\underset{r=1}{\overset{5}{\mathrm{\Σ}}}{\mathrm{\γ}}_r{e}^{-b_r{t}_k}=\left[\begin{array}{cccccc}1& e^{-b_1{t}_k}& e^{-b_2{t}_k}& e^{-b_3{t}_k}& e^{-b_4{t}_k}& e^{{-b}_5{t}_k}\end{array}\right]\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ {\mathrm{\γ}}_2\\ {\mathrm{\γ}}_3\\ {\mathrm{\γ}}_4\\ {\mathrm{\γ}}_5\end{array}\right]$

$=19.9723-9.3344\×e^{-0.0120\×t_k}+0.2691\×e^{-0.0852\×t_k}-2.3181\×e^{-0.4084\×t_k}$

$+2.4242\×e^{-0.6376\×t_k}-0.4552\×e^{-2.7548\×t_k};$

其中，t_k为第k次采样的时间，单位min；

（四）计算该关键点的温度实测值与估计值的均方根误差σ，均方根误差表达式(16)如：

$\mathrm{\σ}=\sqrt{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{[T_e\left(k\right)-T_o\left(k\right)]}^2/(m-1)}---\left(16\right);$

其中，m为采样次数，为1582，T_e(k)为在第k次采样时刻估计的温度值，T_e(k)通过该关键点的温升曲线得到，T_o(k)为第k次采样时刻测量的温度值，计算得到σ=0.0246℃。

3、在采样时间30min下，采样次数m=1801下计算均方根误差σ：

（一）如图2在机床主轴上布置8个温度测量点，在采样间隔1s下得到各个时间下8个温度测量点的温度值，进行30min采样，采样次数m=1801，得到温度采样矩阵为[{T(t₁)},{T(t₂)},……,{T(t₁₈₀₁)}]，具体为：

$\left\{\begin{array}{cccccc}10.45& 10.45& 10.46& 10.46& ......& 13.10\\ 10.16& 10.17& 10.17& 10.17& ......& 13.02\\ 10.49& 10.50& 10.51& 10.51& ......& 13.45\\ 10.31& 10.32& 10.32& 10.33& ......& 12.91\\ 10.24& 10.24& 10.24& 10.25& ......& 12.52\\ 10.37& 10.37& 10.37& 10.38& ......& 12.25\\ 10.41& 10.42& 10.42& 10.43& ......& 12.36\\ 10.35& 10.36& 10.37& 10.37& ......& 12.97\end{array}\right\};$

（二）设定延时时间τ=8s，τ=gΔt，g=8，s为带入计算的次数，s=m-2τ，s=1785，n为温度传感器的数目，n=8；

构造两个温度采样序列：

[Z]_n×s=[{{T(t₁)}-{T(t₁+τ)}},{{T(t₂)}-{T(t₂+τ)}},…,{{T(t_s)}-{T(t_s+τ)}}]_n×s    (3)；

${\left[\hat{Z}\right]}_{n\×s}={\left[\right\{\left\{T(t_1+\mathrm{\τ})\right\}-\left\{T(t_1+2\mathrm{\τ})\right\}\},\{\left\{T(t_2+\mathrm{\τ})\right\}-\left\{T(t_2+2\mathrm{\τ})\right\}\},...,\{\left\{T(t_s+\mathrm{\τ})\right\}-\left\{T(t_s+2\mathrm{\τ})\right\}\left\}\right]}_{n\×s}---\left(4\right)$

${\left[Z\right]}_{8\×1785}={\left[\begin{array}{cccc}-0.01& -0.01& ...& -0.15\\ -0.01& -0.01& ...& -0.09\\ -0.02& -0.02& ...& -0.20\\ -0.02& -0.01& ...& -0.1\\ -0.01& -0.01& ...& -0.07\\ -0.01& -0.02& ...& -0.18\\ -0.02& -0.01& ...& -0.12\\ -0.02& -0.02& ...& -0.10\end{array}\right]}_{8\×1785}---\left(3\right);$

${\left[\hat{Z}\right]}_{8\×1785}={\left[\begin{array}{cccc}-0.01& -0.01& ...& -0.16\\ -0.01& -0.01& ...& -0.10\\ -0.02& -0.02& ...& -0.18\\ -0.02& -0.02& ...& -0.08\\ -0.02& -0.01& ...& -0.10\\ -0.01& -0.02& ...& -0.15\\ -0.02& -0.01& ...& -0.09\\ -0.02& -0.02& ...& -0.16\end{array}\right]}_{8\×1785}---\left(4\right);$

矩阵[H]_8×8满足如下式(8)，

${\left[H\right]}_{8\×8}{\left[Z\right]}_{8\×1785}={\left[\hat{Z}\right]}_{8\×1785}---\left(8\right);$

则，将式(3)和(4)代入式(8)，求解得到[H]_8×8，

${\left[H\right]}_{8\×8}=\left[\begin{array}{cccccccc}0.9530& 0.0743& 0.0281& -0.0820& 0.0913& -0.0220& -0.0106& -0.0430\\ -0.0023& 0.9910& 0.0551& -0.0193& 0.0087& -0.0225& 0.0422& -0.0544\\ 0.0238& 0.1511& 0.9067& -0.1796& 0.1714& 0.0050& -0.0521& -0.0236\\ 0.0350& 0.1192& -0.0700& 0.8285& 0.1230& 0.0248& -0.0409& -0.0128\\ -0.0019& -0.0672& 0.0600& 0.0841& 0.8291& 0.0492& 0.0386& 0.0139\\ 0.0003& 0.0119& -0.0158& 0.0083& 0.0324& 0.9685& -0.0090& 0.0006\\ -0.0132& 0.0522& -0.0307& -0.0279& 0.0436& 0.0138& 0.9425& 0.0145\\ -0.0049& 0.0613& -0.0199& -0.0511& 0.0819& -0.0430& 0.0073& 0.9662\end{array}\right]$

根据式(10)可知，矩阵[H]_8×8的特征值矩阵为特征向量矩阵为对[H]_8×8进行特征值分解，得到

$\left[\begin{array}{cccccccc}0.6200& 0& 0& 0& 0& 0& 0& 0\\ 0& 0.8885& 0& 0& 0& 0& 0& 0\\ 0& 0& 0.9382& 0& 0& 0& 0& 0\\ 0& 0& 0& 0.9658& 0& 0& 0& 0\\ 0& 0& 0& 0& 0.9932+0.0093i& 0& 0& 0\\ 0& 0& 0& 0& 0& 0.9932-0.0093i& 0& 0\\ 0& 0& 0& 0& 0& 0& 0.9981& 0\\ 0& 0& 0& 0& 0& 0& 0& 0.9887\end{array}\right];$

则b₁=3.5853；

则b₂=0.8867；

则b₃=0.4784；

则b₄=0.261；

则b₅=0.0510-0.0703i；

则b₆=0.0510+0.0703i；

则b₇=0.0143；

则b₈=0.0852；

舍去为虚数及负数的b_r，按b_r个数确定p，p=6，并按升序（即b_r热特征值的大小由小到大）重新排列，得到b₁=0.0143min^-1，b₂=0.0852min^-1，b₃=0.261min^-1，b₄=0.4784min^-1，b₅=0.8867min^-1；b₆=3.5853min^-1；

（三）选择温度测试点中温度最高的点为关键点，温度传感器3所对应的温度测量点作为关键点，采样次数为m，该关键点的温度列项量为 ${\left[Q\right]}_{m\×1}=\left[\begin{array}{c}T\left(t_1\right)\\ T\left(t_2\right)\\ .\\ .\\ .\\ T\left(t_m\right)\end{array}\right]=\left[\begin{array}{c}10.49\\ 10.50\\ .\\ .\\ .\\ 13.45\end{array}\right],$ 即将温度采样矩阵中第三行的数据转化为现在的温度列项量数据，其中，T(t₁)为第1次采样得到的该关键点的温度，t₁为第1次采样对应的采样时间，依次类推，T(t_m)为第m次采样得到的该关键点的温度，t_m为第m次采样对应的采样时间；

根据式(15)，

[P]_m×(p+1)X_(p+1)×1=[Q]_m×1    (15)；

其中，m=1801，p=6，

${\left[P\right]}_{1801\×7}\left[\begin{array}{ccccccc}1& e^{-0.0143\×t_1}& e^{-0.0852\×t_1}& e^{-0.261\×t_1}& e^{-0.4784\×t_1}& e^{-0.8867\×t_1}& e^{-3.5853\×t_1}\\ 1& e^{-0.0143\×t_2}& e^{-0.0852\×t_2}& e^{-0.261\×t_2}& e^{-0.4784\×t_2}& e^{-0.8867\×t_2}& e^{-3.5853\×t_2}\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ 1& e^{-0.0143\×t_m}& e^{-0.0852\×t_m}& e^{-0.261\×t_m}& e^{-0.4784\×t_m}& e^{-0.8867\×t_m}& e^{-3.5853\×t_m}\end{array}\right],$

$X_{7\×1}=\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ {\mathrm{\γ}}_2\\ {\mathrm{\γ}}_3\\ {\mathrm{\γ}}_4\\ {\mathrm{\γ}}_5\end{array}\right],$ ${\left[Q\right]}_{1801\×1}=\left[\begin{array}{c}10.49\\ 10.50\\ .\\ .\\ .\\ 13.45\end{array}\right];$

可以计算得到 $\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ {\mathrm{\γ}}_2\\ {\mathrm{\γ}}_3\\ {\mathrm{\γ}}_4\\ {\mathrm{\γ}}_5\\ {\mathrm{\γ}}_6\end{array}\right]=\left[\begin{array}{c}19.3749\\ -9.2617\\ 1.6340\\ -3.9965\\ 4.3417\\ -1.7083\\ 0.2298\end{array}\right];$

得到的温升表达式为：

$T\left(t_k\right)={\mathrm{\γ}}_0+\underset{r=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\γ}}_r{e}^{-b_r{t}_k}=\left[\begin{array}{ccccccc}1& e^{-b_1{t}_k}& e^{-b_2{t}_k}& e^{-b_3{t}_k}& e^{-b_4{t}_k}& e^{-b_5{t}_k}& e^{-b_6{t}_k}\end{array}\right]\left[\begin{array}{c}{\mathrm{\γ}}_0\\ {\mathrm{\γ}}_1\\ {\mathrm{\γ}}_2\\ {\mathrm{\γ}}_3\\ {\mathrm{\γ}}_4\\ {\mathrm{\γ}}_5\\ {\mathrm{\γ}}_6\end{array}\right]$

$=19.3749-9.2617\×e^{-0.0143\×t_k}+1.634\×e^{-0.0852\×t_k}-3.9965\×e^{-0.261\×t_k}$

$+4.3417\×e^{-0.4784\×t_k}-1.7083\×e^{-0.8876\×t_k}+0.2298\×e^{-3.5853\×t_k};$

其中，t_k为第k次采样的时间，单位min。

（四）计算该关键点的温度实测值与估计值的均方根误差σ，均方根误差表达式(16)如：

$\mathrm{\σ}=\sqrt{\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{[T_e\left(k\right)-T_o\left(k\right)]}^2/(m-1)}---\left(16\right);$

其中，m为采样次数，为1801，T_e(k)为在第k次采样时刻估计的温度值，T_e(k)通过该关键点的温升曲线得到，T_o(k)为第k次采样时刻测量的温度值，计算得到σ=0.775℃。

如表1可以清晰地观察到在30min内，均方根误差存在极小值，在极小均方根值对应的采样时间为26.35min，作为辨识时间，该辨析时间下关键点的温升曲线即为能准确反映机床主轴热态特性的曲线，如图3所示a曲线所示。

如图3所示，将26.35min辨识时间的关键点的温升曲线a与该关键点的实测温度曲线b比较，两者呈现高度的拟合。辨识的稳态温度为19.955℃，辨识的热平衡时间为185.4min，而实测的稳态温度为19.8℃，实测的热平衡时间为184.2min。辨识的热态特性与实测相比误差在1%以内，因此，本发明结果准确可靠，此发明方法能够在快速辨识主轴热特性上取得很好效果。

Claims (3)

1.一种快速辨识机床主轴热态特性的方法，其特征在于，包括以下步骤：

(一)在机床主轴上布置n个温度测量点，在采样间隔Δt下得到各个时间下n个温度测量点的温度值，采样次数为m，得到温度采样矩阵为[{T(t1)},{T(t2)},……,{T(tm)}]，具体为：

其中，1,2,...,n代表不同的各温度测量点；

(二)设定延时时间τ，τ＝gΔt，g为正整数，s＝m-2τ；

构造两个温度采样序列：

[Z]_n×s＝[{{T(t₁)}-{T(t₁+τ)}},{{T(t₂)}-{T(t₂+τ)}},…,{{T(t_s)}-{T(t_s+τ)}}]_n×s       (3)；

矩阵[H]_n×n满足如下式(8)，

则，将式(3)和(4)代入式(8)，求解得到[H]_n×n，

根据式(10)可知，

矩阵[H]_n×n的特征值矩阵为特征向量矩阵为 对矩阵[H]_n×n进行特征值分解，求出热特征值b_r(r＝1,2,...,p)，舍去热特征值b_r为虚数及负数的b_r，按保留的热特征值b_r个数重新确定p，并将保留的热特征值按热特征值b_r的大小由小到大重新排列；

(三)选择温度测试点中温度最高的点为关键点，采样次数为m，该关键点的温度列向量为其中，T(t₁)为第1次采样得到的该关键点的温度，t₁为第1次采样对应的采样时间，依次类推，T(t_m)为第m次采样得到的该关键点的温度，t_m为第m次采样对应的采样时间；

根据式(15)，

[P]_m×(p+1)[X]_(p+1)×1＝[Q]_m×1       (15)；

其中，

计算得到

该关键点的温升曲线表达式如式(14)所示：

其中，T(t_k)表示t_k时刻的温度，t_k为第k次采样的时间，γ_r(r＝0,1,2,...,p)表示温升曲线表达式系数，b_r(r＝1,2,...,p)表示热特征值；

(四)计算该关键点的温度实测值与估计值的均方根误差σ，均方根误差表达式(16)为：

其中，m为采样次数，T_e(k)为在第k次采样时刻估计的温度值，T_e(k)通过该关键点的温升曲线得到，T_o(k)为第k次采样时刻测量的温度值；

(五)在某一采样次数m下得到均方根误差σ，采样次数m依次增大，返回步骤(一)，得到不同采样次数m下的均方根误差σ，直至均方根误差σ出现先减小后增大的情形，即均方根误差σ存在极小均方根值，即在极小均方根值对应的采样次数m所对应的时间为辨识时间，该辨识时间下关键点的温升曲线即为能准确反映机床主轴热态特性的曲线。

2.根据权利要求1所述的快速辨识机床主轴热态特性的方法，其特征在于，步骤(一)中，所述各温度测量点沿机床主轴的轴向布置。

3.根据权利要求1所述的快速辨识机床主轴热态特性的方法，其特征在于，步骤(一)中，采样次数m要大于温度测量点的个数n。