0 引言
对于数控机床热误差补偿，其补偿模型精度是影响最终补偿效果的主要因素[1-4]。许多研究发现，补偿模型的热误差拟合精度高，但并不代表模型能够长期准确地预测机床热误差[5-8]。
人们对数控机床热误差补偿控制技术的研究多集中在建模理论与应用方面。LO等[9]采用相关性分组、典型变量搜索和分组搜索等方法，优化了温度测点选择，保证了机床热误差多元线性回归模型的预测精度。ABDULSHAHED等[10]通过灰色模型和模糊C-均值聚类的温度测点优化选择方法，减少了温度测点数量，降低了温度传感器的不确定性及建模自变量之间的相关性，并应用自适应神经模糊推理系统建立了机床热误差模型。苗恩铭等[11-12]提出模糊聚类结合灰色关联度的温度测点选择方法，减少了温度测点数量，降低了建模自变量之间的相关性和多重共线性，建立了机床热误差支持向量机和多元线性回归模型，该模型具有良好的预测精度和稳健性；并进一步提出了机床热误差补偿建模最佳转速状态的选择方法[13]，比对分析了机床在主轴空转和实切状态下的热特性差异，且基于田口实验方法实现了实切状态下热误差补偿实验方案的优化设计[14]。张伟等[15]采用灰色关联分析和模糊聚类分析结合的方法对机床温度测点进行优化选择，保证了机床热误差模型的预测精度。李艳等[16]采用互信息和改进模糊聚类相结合的方法对机床温度测点进行优化选择，显著提高了机床热误差模型预测精度。上述研究通过对建模温度测点的优化选择，有效提升了模型准确反映机床热误差与温度之间变化规律的能力，从而保证了热误差模型的预测精度。通常，选出的建模温度测点称为温度敏感点。目前，模糊聚类结合灰色关联度算法是应用较为广泛的温度敏感点选择方法[5,11-12,15-16]。
我们通过大量实验研究发现，当机床运行参数发生变化时，热误差与温度之间的变化规律也会随之改变。当进行热误差预测时，机床很可能会在不同的参数设置条件下运行，热误差规律的差异性也会导致模型预测精度降低。为此，本文在目前传统热误差建模算法研究的基础上，提出了基于状态空间模型的数控机床热误差建模方法。该方法综合考虑了温度敏感点温升、主轴转速对机床热误差的影响，使模型可根据机床运行参数的变化而自动调整，且模型对机床运行参数变化具有良好的自适应性，从而为提升机床热误差模型的预测精度提供了有利条件。
1 机床运行参数变化对热误差模型预测精度的影响
为验证机床热误差状态空间模型的预测精度，本文以一台Leaderway V-450型数控机床为实验对象，在不同主轴转速条件下进行实验。实验结果表明：当采用传统算法(本文以多元线性回归算法为例)在一定转速下建立热误差模型,对不同转速下的实验数据进行预测时，相比于对相同转速下的实验数据进行预测，其模型预测精度降低，表明机床运行参数的变化确实会导致热误差规律发生变化。对机床热误差多元线性回归模型和状态空间模型的预测精度进行比对分析，实验结果表明：状态空间模型具有更高的预测精度和稳健性，从而验证了机床热误差状态空间模型的优越性。
1.1 模糊聚类结合灰色关联度算法
为降低各温度变量之间的相关性和多重共线性对模型预测精度的影响，在进行机床热误差建模时，首先应进行建模温度测点的优化选择(即温度敏感点选择)。本文采用模糊聚类结合灰色关联度的温度敏感点选择方法，有效降低了各温度变量之间的相关性和多重共线性，为机床热误差模型预测精度的提高提供了保障，具体如下。
(1)模糊聚类。①计算温度变量xi、xj之间的相关系数rij，确定模糊相似矩阵R=(rij)n×n，其中，i,j=1，2，…,n(n为温度测点总数)；②将模糊相似矩阵R换算成模糊等价矩阵，可采用平方法求出R的传递闭包
为一个模糊等价矩阵(即R→R2→(R2)2→…→R2k)，存在某个整数k使得R2k=R2k+1，则传递闭包
③根据精度要求，在[0,1]之间选取阈值λ，若传递闭包
中的元素
则该元素为1，否则，该元素为0，并得到布尔矩阵
若元素
则将温度变量xi、xj归为一类，从而达到分类目的，以降低不同类温度变量之间的多重共线性。相关系数rij按下式计算：
(1)
式中，t为不同时刻，是离散的测量时间点；L为实验数据长度。
(2)灰色关联度。采用邓氏关联度计算表达式获得机床热误差y与各温度变量xi(i=1,2,…,n)之间的灰色关联度γ(y,xi)，灰色关联度越大，温度测点对机床热误差的影响权重越大。选择出每类中与热误差关联度最大的温度测点作为温度敏感点，并参与建模，从而避免同类中温度变量间的强关联度对热误差模型预测精度的影响，同时，减少了使用的温度传感器数量。灰色关联度γ(y,xi)按下式计算：
(2)
r(y(t),xi(t))=
(3)
式中，ρ为分辨系数，ρ∈[0,1]，一般取ρ=0.5。
对温度敏感点进行选择后，将其作为自变量，通过多元线性回归算法即可建立机床热误差模型：
y=β0+β1x1+…+βmxm
(4)
式中，x1、x2、…、xm分别为不同的温度敏感点温升；β0、β1、…、βm为模型参数。
1.2 机床热误差实验
转速是机床主要的运行参数之一，为减少外界环境对实验结果的影响，在同一季度内，本文通过改变主轴转速来使机床处于不同的运行参数条件下。以Z向热误差为例，实验选择温度传感器DS18B20(测量精度为±0.5 ℃，最大分辨能力为0.062 5 ℃)采集温度数据，参考国际标准ISO 230-3:2007[17]给出的热误差测量方法，选择电涡流传感器(测量精度为±1 μm)采集机床主轴Z向热误差。温度传感器的贴放位置以影响机床主轴Z向热误差的主要热源附近为主，温度传感器T1～T9(对应的温度敏感点分别为T1～T9)和电涡流传感器S的安放位置见图1。温度传感器T10安放在机床外壳上以测量环境温度，未出现在图1中。各传感器的安放位置及作用见表1。实验装置如图2所示。
图1 传感器安放位置示意图
Fig.1 Sketch map of sensor placement
表1 传感器安放位置及作用
Tab.1 Installation locations and functions of sensors
传感器安放位置作用T1～T5主轴前端测量主轴发热T6,T9主轴套测量主轴发热T7打刀缸底座测量主轴发热T8主轴电机测量电机发热T10机床外壳测量环境温度SZ向正下端测量主轴热误差
图2 实验装置实物图
Fig.2 Physical map of experimental device
测量时，主轴以恒定转速(2 000 r/min、4 000 r/min、6 000 r/min)转动，每隔3 min主轴停转一次，以测量主轴Z向热误差，并通过温度测量系统采集该时刻的温度数据。测量实验持续时间均达4 h以上。
本研究共进行了6批次实验，分别记为K1～K6，各批次实验主轴转速的设置见表2，测得6批次实验数据的热误差曲线见图3。以K1批次实验数据为例，图4给出了各温度曲线。
表2 各批次实验主轴转速设置
Tab.2 Spindle speed setting of each batch test
批次主轴转速(r/min)K1,K22 000K3,K44 000K5,K66 000
图3 K1～K6批次实验数据的热误差曲线
Fig.3 Thermal error curves of K1～K6 batch experimental data
图4 K1批次实验数据的温度曲线
Fig.4 Temperature curves of K1batch experimental data
1.3 机床运行参数变化条件下热误差模型预测精度分析
采用模糊聚类结合灰色关联度算法对K1～K6批次实验数据的温度敏感点进行选择，温度敏感点的选择结果均为T1、T7。采用多元线性回归算法，分别获得K1～K6批次实验数据的机床Z向热误差多元线性回归模型，分别记为模型M1～M6。
M1：yZ=-1.67+9.82ΔT1-4.09ΔT7
M2：yZ=-3.03+8.65ΔT1-2.57ΔT7
⋮
M6：yZ=-4.01+7.15ΔT1-3.81ΔT7
式中，yZ为机床Z向热误差测量值；ΔT1、ΔT7分别为温度敏感点T1和T7的温升。
进而分别利用模型M1～M6对K1～K6批次实验数据进行拟合或预测，并计算出各模型对每批次实验数据的拟合或预测残余标准差S[18]。残余标准差S越小，表明拟合或预测误差越小，精度越高。残余标准差S的计算公式为
(5)
式中，
为机床Z向热误差估计值;h为自变量个数。
同一转速下的两个模型对同一转速下两批次数据的预测残余标准差S的平均值记为SCM，即
(6)
式中，I、J分别为转速相同的模型编号和实验批次。
SCM的计算结果见表3，并将表3中的计算结果绘制成条形图，如图5所示。
表3 不同转速下模型对各转速热误差预测残余标准差平均值
Tab.3 Average of predicted residual standard deviation of thermal error for each rotational speed by different speed models
预测批次主轴转速np(r/min)2 0004 0006 000建模批次主轴转速nm(r/min)2 0001.968.5621.624 0003.882.519.106 0006.826.011.66
图5 不同转速下模型对各转速热误差预测残余标准差平均值
Fig.5 Average of predicted residual standard deviation of thermal error for each rotational speed by different speed models
从图5中可以看出，机床热误差多元线性回归模型对相同转速(nm=np)下得到的热误差数据的预测残余标准差平均值，明显小于对不同转速(nm≠np)下得到的热误差数据的预测残余标准差平均值，这表明当数控机床转速发生变化时，会导致机床热误差规律发生改变，从而引起热误差模型预测精度降低。
2 机床热误差状态空间模型
状态空间模型的系统输出包含系统输入、系统自身特性、系统同外界相互联系三者信息，因此状态空间模型不仅可以描述系统输入与输出之间的关系，还能更好地揭示系统的内部结构和特性，对系统运行参数变化具有良好的适应能力[19]。为提升机床热误差模型在机床不同运行参数条件下的模型预测精度，本文建立以温度敏感点温升与主轴转速为输入的机床热误差状态空间模型，用于数控机床热误差预测。
2.1 机床热误差状态空间模型建立
将温度敏感点的温升ΔT1、ΔT7和该批次实验数据主轴转速nm作为输入，机床Z向热误差yZ(t)作为输出，建立机床热误差状态空间模型如下：
x(t+1)=Ax(t)+Bu(t)
(7)
yZ(t)=Cx(t)+Du(t)
(8)
u(t)=(ΔT1(t),ΔT7(t),nm(t))T
x(t)=(x1(t),x2(t),x3(t))T
C=(1,0,0) D=(d1,d2,d3)
式中，u(t)为系统的输入向量；x(t)为系统自身的状态向量;A为系统自身的结构矩阵；B为系统的参数矩阵；C、D均为系统的参数向量;a1、a2、a3为矩阵A中的元素数值；b11、b12、…、b33为矩阵B内的元素数值；d1、d2、d3为向量D中的元素数值。
2.2 机床热误差状态空间模型参数辨识
本文用最小二乘算法估计模型的参数，用状态观测器估计系统的状态[20]，基本过程如下。
(1)联合式(7)和式(8)可得到机床Z向热误差状态空间模型的参数辨识模型：
yZ(t)=x1(t)+Du(t)=HT(t)a+TT(t)b
(9)
a=(a1,a2,a3)T
b=(d1,d2,d3,b11,b21,b31,b12,b22,b32,b13,b23,b33)T
H(t)=[-x3(t-3) -x2(t-3) -x1(t-3)]T
T(t)=[ΔT1(t) ΔT7(t) nm(t) ΔT1(t-1)
ΔT1(t-2) ΔT1(t-3) ΔT7(t-1) ΔT7(t-2)
ΔT7(t-3) nm(t-1) nm(t-2) nm(t-3)]T
定义相关堆积矩阵Y=[yZ(1) yZ(2) …
yZ(L)]T，H=[H(1) H(2) … H(L)]T，T=[T(1) T(2) … T(L)]T，则式(9)可等价为
(10)
因此，可获得基于最小二乘算法的机床Z向热误差状态空间模型参数估计，即
(11)
由于式(11)对应的求逆矩阵在每次计算过程中均含不变的子矩阵和变化的子矩阵，为减少计算量，本文利用分块矩阵求逆引理，得到机床Z向热误差状态空间模型基于最小二乘算法的参数估计等价表达式，即
(12)
α=βTTY β=inv(TTT)
Q=HTMH M=I-TβTT
式中，I为L行L列的单位矩阵。
(2)本文通过设置初始状态向量x(1)=(1,1,1)T/p0(此处参数p0取106)，并利用上述热误差状态空间模型参数估计结果，依据式(7)(即状态观测器)获得热误差状态空间模型中未知的状态向量x(t)=(x1(t),x2(t),x3(t))T。
(3)在给定辨识精度q0条件下，利用迭代方案获得模型参数向量CO=(a,b)T，并建立机床Z向热误差状态空间模型。机床Z向热误差状态空间模型参数辨识流程见图6，其中q=‖CO,k-CO,k-1‖，CO,k-1、CO,k分别为第k-1次和第k次迭代热误差状态空间模型参数向量CO的估计结果。
图6 机床热误差状态空间模型参数辨识流程图
Fig.6 Machine thermal error state space model parameter identification flow chart
根据上述步骤，建立K1～K6批次实验数据的机床Z向热误差状态空间模型N1～N6，结果如下。
Ni：x(t+1)=Aix(t)+Biu(t)
yZ(t)=Cix(t)+Diu(t) i=1，2，…，6
⋮
⋮
C1=C2=…=C6=(1,0,0)
D1=(0.37，0.19，4.60×10-5)
D2=(0.38，0.26，2.26×10-4)
⋮
D6=(0.55，-0.21，1.07×10-3)
3 预测精度比对
为了对机床热误差多元线性回归模型、状态空间模型进行比对分析，分别利用多元线性回归模型M1～M6及状态空间模型N1～N6对K1～K6批次热误差数据进行拟合或预测。根据式(5)计算得到各模型对每批次实验数据的预测残余标准差S，并计算各模型对应所有预测残余标准差(可反映模型预测精度)S的均值参数SM和离散标准差参数SS，计算结果分别见图7和图8。SM、SS计算公式分别为
(13)
(14)
其中，对于M1～M6、N1～N6中的每个模型，Si为该模型对Ki(i=1,2,…,6)批次实验数据的预测残余标准差。模型预测精度均值参数SM小，表明模型的平均预测精度高；模型预测精度离散标准差参数SS小，表明模型对各批次热误差数据预测精度的分散程度小，各批次实验转速变化对预测精度的影响较小。
图7 机床热误差模型预测精度均值
Fig.7 Average of prediction accuracy of machine tool thermal error model
图8 机床热误差模型预测精度离散标准差
Fig.8 Discrete standard deviation of prediction accuracy of machine tool thermal error model
由图7可以看出，多元线性回归模型预测精度均值参数SM分布在4.75～12.33 μm，状态空间模型预测精度均值参数SM分布在2.63～3.30 μm。由图8可以看出，多元线性回归模型预测精度离散标准差参数SS分布在2.43～10.28 μm，状态空间模型预测精度离散标准差参数SS分布在0.64～1.62 μm。相对于多元线性回归模型，状态空间模型的预测精度均值和预测精度离散标准差均有所降低。
为综合对比分析两种模型的预测精度及稳健性，将机床热误差多元线性回归模型M1～M6和状态空间模型N1～N6的预测精度均值参数SM和预测精度离散标准差参数SS进行平均，得到两种模型的总预测精度均值参数
和预测精度离散标准差均值参数
见表4。
表4 两种模型的平均预测效果
Tab.4 Average of predicted residual standard deviation of thermal error for each rotational speed by different speed models μm
多元线性回归模型M1～M6状态空间模型N1～N6总预测精度均值SM6.902.89预测精度离散标准差均值SS5.040.97
由表4可知，在不同主轴转速下，相对于机床热误差多元线性回归模型，机床热误差状态空间模型的预测精度及稳健性得到了较好的提升。经计算，机床热误差状态空间模型的模型预测精度总体提升58.12%。由上述分析可知，相对于机床热误差多元线性回归模型，机床热误差状态空间模型实现了机床热误差、温度敏感点温升及主轴转速之间作用关系的综合描述，使模型可根据机床运行参数变化而自动调整，模型对机床运行参数变化具有良好的自适应性，较好地提升了模型的预测精度和稳健性。
4 实验验证
为验证数控机床热误差状态空间模型补偿效果，将模型N1存储到机床补偿卡系统中，分别在热误差补偿控制系统关闭、开启状态下，在Leaderway V-450型数控机床上，采用可转位面铣刀对铸铁进行平面切削加工实验。加工实验开始前，主轴以转速为6 000 r/min空转3 h；加工过程中，主轴转速1 000 r/min、进给速度300 mm/min、切削深度300 μm，切削过程持续15 min。加工试件见图9，其中，①为参考平面，②为未补偿加工平面，③为补偿加工平面。采用德国蔡司MC850三坐标测量机(测量范围为850 mm×700 mm×600 mm，分辨能力为0.1 μm)对切削加工试件进行均布9点检测，测试结果见表5。
图9 加工试件示意图
Fig.9 Processing specimen schematic
表5 平面切削加工实验补偿效果
Tab.9 The compensation effects of plane cutting experiments
μm
实际切削深度平均误差量关闭346.346.3开启306.76.7
由表5可知，相对于未补偿时的平均误差量，补偿时的平均误差量减少了85.5%。结果表明：基于状态空间模型的机床热误差补偿技术的实施，可获得很好的加工补偿效果，有效地提高了数控机床加工精度。
5 结论
(1)针对目前关于热误差建模的研究中忽略了机床运行参数变化对模型精度影响的问题，本文提出了基于状态空间模型的数控机床热误差建模方法。
(2)以不同主轴转速条件下的Leaderway V-450型数控加工中心为研究对象，分别建立机床热误差多元线性回归模型和状态空间模型，结果表明：与多元线性回归模型相比较，状态空间模型的预测精度及稳健性得到了有效提升。
(3)采用Leaderway V-450型数控加工中心进行了实际切削加工条件下的热误差补偿实验，结果表明：机床热误差状态空间模型的应用补偿效果显著。
[1] KIM B S, SONG Y C, PARK C H. Robust Thermal Error Modeling and Compensation for a Nano Level Thermal Drift in a High Precision Lathe[J]. International Journal of Precision Engineering and Manufacturing, 2011, 12(4): 657-661.
[2] WANG H T, WANG L P, LI T M. Thermal Sensor Selection for the Thermal Error Modeling of Machine Tool Based on the Fuzzy Clustering Method[J]. International Journal of Manufacturing Technology, 2013, 69: 121-126.
[3] HAN J, WANG L P, WANG H T, et al. A New Thermal Error Modeling Method for CNC Machine Tools[J]. International Journal of Manufacturing Technology, 2012, 62: 205-212.
[4] ZHANG T, YE W H, LIANG R J. Temperature Variable Optimization for Precision Machine Tool Thermal Error Compensation on Optimal Threshold[J]. Chinese Journal of Mechanical Engineering, 2013, 26(1): 158-165.
[5] MIAO E M, GONG Y Y, NIU P C， et al. Robustness of Thermal Error Compensation Modeling Models of CNC Machine Tools[J]. International Journal of Manufacturing Technology, 2013, 69: 2593-2603.
[6] ABDULSHAHED A M, LONGSTAFF A P, FLETCHER S， et al. Thermal Error Modeling of Machine Tools Based on ANFIS with Fuzzy C-means Clustering Using a Thermal Imaging Camera[J]. International Journal of Manufacturing Technology, 2015, 39(7): 1837-1852.
[7] YANG L, ZHAO W H, WU W W, et al. Thermal Error Modeling of the Spindle Based on Multiple Variables for the Precision Machine Tool[J]. International Journal of Manufacturing Technology, 2014, 72: 1415-1427.
[8] XIANG S T, LU H X, YANG J. Thermal Error Prediction Method for Spindles in Machine Tools Based on a Hybrid Model[J]. Journal of Engineering Manufacture, 2015, 229(1): 130-140.
[9] LO Chihhao, YUAN Jingxia, NI Jun. Optimal Temperature Variable Selection by Grouping Approach for Thermal Error Modeling and Compesation[J]. International Journal of Machine Tools & Manufacture, 1999, 39(9): 1383-1396.
[10] ABDULSHAHED A M, LONGSTAFF A P, FLETCHER S， et al. The Application of ANFIS Prediction Models for Thermal Error Compensation on CNC Machine Tools[J]. Applied Soft Computing, 2015, 27: 158-168.
[11] 苗恩铭，龚亚运，成天驹，等. 支持向量回归机在数控加工中心热误差建模中的应用[J]. 光学精密工程，2013, 21(4): 980-986.
MIAO Enming, GONG Yayun, CHENG Tianju, et al. Application of Support Vector Regression Machine to Thermal Error Modelling of Machine Tools[J]. Opt. Precision Eng., 2013, 21(4): 980-986.
[12] MIAO E M, GONG Y Y, DANG L C, et al. Temperature-sensitive Point Selection of Thermal Errors of Model of CNC Machine Center[J]. International Journal of Advanced Manufacturing Technology, 2014, 74:681-691.
[13] 苗恩铭，吕玄玄，苗继超，等. 数控机床热误差补偿最佳转速选择[J]. 光学精密工程，2015, 23(11): 3176-3182.
MIAO Enming, LYU Xuanxuan, MIAO Jichao, et al. Selection of Optimum Spindle Speed to Thermal Error Compensation of Machine Tools[J]. Opt. Precision Eng., 2015, 23(11): 3176-3182.
[14] MIAO E M, LIU H, FAN K C, et al. Analysis of CNC Machining Based on Characteristics of Thermal Errors and Optimal Design of Experimental Programs [J]. International Journal of Advanced Manufacturing Technology, 2017, 88(5): 1363-1371.
[15] 张伟，叶文华. 基于灰色关联和模糊聚类的机床温度测点优化[J]. 中国机械工程，2014, 25(4): 456-460.
ZHANG Wei, YE Wenhua. Optimization of Temperature Measuring Points for Machine Tools Based on Grey Correlation and Fuzzy Clustering Analysis[J]. China Mechanical Engineering, 2014, 25(4): 456-460.
[16] 李艳，李英浩，高峰，等. 基于互信息法和改进模糊聚类的温度测点优化[J]. 仪器仪表学报，2015, 36(11): 2466-2472.
LI Yan, LI Yinghao, GAO Feng, et al. Investigation on Optimization of Temperature Measurement Key Points Based on Mutual Information and Improved Fuzzy Clustering Analysis[J]. Chinese Journal of Scientific Instrument, 2015, 36(11): 2466-2472.
[17] International Standards Organization. ISO 230-3:2007. Test Code for Machine Tools—Part 3: Determination of Thermal Effects[S]. Geneva: ISO, 2007.
[18] 费业泰. 误差理论与数据处理[M]. 6版. 北京：机械工业出版社，2010.
FEI Yetai. Error Theory and Data Handle [M]. 6th ed. Beijing: China Machine Press, 2010.
[19] L K 蒂莫西, B E 博纳. 状态空间分析导论(下册)[M]. 胡钦训，刘颖，译. 北京：高等教育出版社，1986.
TIMOTHY L K, BONA B E. State Space Analysis: an Introduction(Volume Ⅱ)[M]. HU Qinxun, LIU Ying， Trans. Beijing: Higher Education Press, 1986.
[20] HU H, DING F. An iterative Least Squares Estimation Algorithm for Controlled Moving Average Systems Based on Matrix Decomposition [J]. Applied Mathematics Letters， 2012, 25: 2332-2338.
