CN103543637B - 一种重型机床环境温度解析建模方法 - Google Patents

Abstract

本发明提供一种重型机床车间环境温度解析建模方法，其结合时间序列分析，傅立叶级数分解方法，将实测的环境温度用傅立叶三角级数形式表达，获取环境温度的时间序列、波动频率、波动幅值、相位等信息，实时测量更新的温度数据和当前时间信号作为输入，实现温度的实时预测，代替实测温度用于热误差响应预测建模。本发明考虑了环境温度波动的周期性和非周期性，以及波动随季节和年度变化的特征，同时考虑了当前机床当前热状态既受当前环境温度影响又受历史温度状态影响的客观事实，有利于更准确地定量描述机床热变形误差的时滞响应特性，进而提高机床热误差预测的精度和鲁棒性。

Description

一种重型机床环境温度解析建模方法

技术领域

本发明涉及数控机床热误差建模补偿领域，更具体地，涉及放置于非恒温加工车间的重型数控机床环境温度解析建模方法。

背景技术

重型数控机床所处的环境温度对机床的热变形影响很大，在没有温控条件的车间，环境温度波动具有缓变、近似周期的特征，每日最低和最高温度出现的时间随季节波动而变化，同时重型数控机床对环境温度波动的热变形响应具有时间滞后，且滞后时间也随季节温度改变。对环境温度长期波动特征的数学描述，即对重型机床环境温度进行解析建模，成为准确建立温度与机床热误差之间鲁棒性模型的重要前提。

热误差建模通常是指建立机床的热变形与对应温度之间的映射关系，常见的热误差预测模型有：最小二乘模型、贝叶斯网络模型、支持向量机模型、灰色***模型，人工智能模型、模糊***模型等。对应的温度信息包括床身温度和环境温度等，通常确定建模所需温度测点的方法是：先根据经验在机床和环境布置大量温度测点，然后通过选择、分组、优化等，获得最优或最近线性布点。当前的研究中，环境温度对机床热误差影响的研究局限于实验统计和定性分析。在减少环境温度对机床的结构热变形的影响方面，主要采用的方法是对机床热结构优化、环境温度控制、改善散热条件、优化热源布置等。

随着最优布点识别方法和热误差非线性建模理论的发展，对机床内部热源的影响机理和热误差减少方法的研究日趋成熟，通过对机床本体温度测点的布置、分组、优化等，可以获得数量尽量少的最优或最近线性测点。

在内部热源引起的热变形机理分析逐步完善时，环境温度作用于机床引起的热误差的成分逐渐上升，而且没有得到有效解决。因此从机理和解析建模角度，准确建立机床热变形与环境温度之间的映射关系，是提高热误差补偿能力并进而提高机床加工精度的前提。

环境温度与机床本体的最优布点影响特征不同，机床对环境温度波动的响应一直存在滞后，不存在最近线性环境温度布点，而且滞后的时间随气候和季节变化而改变。但是当前的线性和非线性的热误差预测模型不能体现随温度滞后特征，导致当前的热误差预测模型在不同环境的通用性有限，甚至补偿精度失效。另一方面，当前对环境温度的研究和分析局限于实验测量、现象说明和定性讨论。目前的文献检索未发现将环境温度的非线性和时滞特征定量用于热误差建模或补偿的论述，也没有发现关于重型机床环境温度解析建模方面的研究。

发明内容

由上可知，当前需要建立能够反映重型机床环境温度非线性滞后特性的热误差模型，从而提高模型预测的精度和鲁棒性。作为模型输入的环境温度变量是实时测量的温度，不包含历史的温度信息和热信息，而机床的热误差既受到当前环境温度的影响，还同过去的机床状态以及过去的环境状态有关。环境温度是可以实时测量的，但是实测的环境温度直接用于时滞热误差预测会使得预测精度降低。

因此，当前需要解决的技术问题是：提供一种能够实时反映环境温度的解析模型，该模型具有明确的数学公式，能够准确描述车间环境温度的波动和随季节温度变化的特征，而且能够替代实测环境温度，用于时滞热误差响应的预测建模。

为了解决上述技术问题，本发明采用的技术方案是：

结合时间序列分析，傅立叶级数分解方法，将实测的环境温度用傅立叶三角级数形式表达，获取环境温度的时间序列、波动频率、波动幅值、相位等信息。实时测量更新的温度数据和当前时间信号作为输入，实现温度的实时预测，代替实测温度用于热误差响应预测建模。

一种重型机床环境温度解析建模方法，其具体步骤是：

1）环境温度傅立叶级数分解

环境温度t_e(x)是连续物理量，具有基本的日周期波动特征（周期T₀=1440分钟），满足狄里赫利条件，根据傅立叶定律，环境温度可以展开为傅立叶级数形式，即：

$t_e\left(x\right)=A_0+\underset{n=1}{\overset{\∞}{\mathrm{\Σ}}}{A}_n\sin (n{\mathrm{\ω}}_{0}x+{\mathrm{\φ}}_n),n=\mathrm{1,2,3},L---\left(1\right)$

其中，x表示时间，单位为分钟，它由时间参考起点和当前时刻确定；基频ω₀=2π/T₀，A₀为温度的均值项，A_n是各个阶次的温度波的振幅，n=1,2,…，φ_n为相位角，其物理意义为各阶滞后时间；

实际上，由于受到户外不确定的气候变化的影响，重型机床所处的环境温度不具备绝对的日周期特征，而是包含年周期和非周期特性。因此公式（1）不能直接用于温度分解和预测，而要对A₀,A_n,φ_n等各项参数进行调整和替换，使得公式（1）具有明确的物理意义，并且能够用于实测环境温度的替代；

2）均值项A₀

均值项A₀的物理意义是当前温度点所在的周期的均值，由环境温度于不是绝对的日周期重复而是缓慢变化的，且当前温度所在周期的均值受历史温度状态的影响，因此可以通过滑动平均的方法，将当前测点之前一个周期的据历史温度数据的平均值来代替当前所在周期的均值。也就是说：

$A_0=A_0^{\left(1\right)}\left(x\right)=\frac{t_{\left[x\right]-N}+t_{\left[x\right]-(N-1)}+L+t_{\left[x\right]-1}}{N},\left[x\right]\≥N---\left(2\right)$

在此，[x]是对任意时刻x取整数，t_[x]是在时间点x的测量温度，N是一个周期的采样点数，环境温度是缓变信息，可以认为时间[x]≤x＜[x]+1期间的环境温度等同于t_[x]。在滑动平均的方法，要求至少有一个周期N个采样点的历史数据来预测当前的温度。新的测量值t_[x]-1引入，旧的t_[x]-N退出，于是预测的均值A₀可以在线更新。滑动平均方法保留了所有低于基频频率的波动信息；

3）幅值项A_n

A_n为分解后各阶次温度波动的幅值项，可以认为其与所考虑周期内波动最大值T_max与平均值A₀之差成正比，也就是：

$A_n={\mathrm{\β}}_n(T_{\max }\left(x\right)-A_0^{\left(1\right)}\left(x\right)),n=\mathrm{1,2,3},L---\left(3\right)$

其中：

$T_{\max }=\left(x\right)=\underset{\left[x\right]-N<u<\left[x\right]-1}{\max }\left(t_{\left[u\right]}\right)---\left(4\right)$

β_n含义为不同频率成分温度波对总的温度变化的贡献成分即权重，表示车间的固有特征。(T_max(x)-A₀ ⁽¹⁾(x))在不同的观测周期变化的，但对已知的环境温度历史观测值是确定的；

4）相位φ_n

实际温度的波动部分的最低值出现的时间随环境温度的变化而变化，夏季高温的时候，波动部分最低值出现的时间早，冬季低温时，波动部分最低值出现的时间晚。波动部分最低值出现的时间，是各个频率成分的相位共同作用的效果。因此，在傅立叶分解的形式中，相位由两部分组成：定值部分和波动部分，其中定值部分的叠加效果构成了夏季高温时的波动特征，而变化部分描述了最低温度出现的时间随环境温度整体波动的规律。综上分析，相位φ_n可以表示为：

${\mathrm{\&phi;}}_n={\mathrm{\&phi;}}_{0n}+{\mathrm{\&alpha;}}_n\&CenterDot;A_0^{\left(1\right)}\left(x\right),n=\mathrm{1,2,3},L---\left(5\right)$

在此α_n表示季节变换的滞后时间系数，代表了特定车间的固有热特性，可以认为是一个特定值；

5）总的温度预测模型

综上所述，任意x时刻的环境温度预测的解析模型可表示为：

$t_e\left(x\right)=A_0^{\left(1\right)}\left(x\right)+\underset{n=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_n(T_{\max }\left(x\right)-A_0^{\left(1\right)}\left(x\right))\sin (n{\mathrm{\&omega;}}_{0}x+({\mathrm{\&phi;}}_{0n}+{\mathrm{\&alpha;}}_n\&CenterDot;A_0^{\left(1\right)}\left(x\right))),n=\mathrm{1,2,3},L---\left(6\right)$

式中，均值A₀ ⁽¹⁾(x)和最大值T_max(x)可以通过至少一个周期的历史测量值计算得到，而且是可以实时更新的，基频ω₀是已知的，因此可以由实测温度数据用非线性最小二乘原理辨识上述解析模型中β_n,φ_0n和α_n等参数，进而实现环境温度的在线预测。该解析预测模型具有实际的物理意义，而且包含频率和相位成分，便于实时热误差响应分析。

本发明的有益效果是：本发明考虑了环境温度波动的周期性和非周期性，以及波动随季节和年度变化的特征，同时考虑了当前机床当前热状态既受当前环境温度影响又受历史温度状态影响的客观事实，提出新的研究思路，将结合时间序列分析，实现了环境温度的傅立叶级数分解，将具有明确时间、频率、相位信息的解析形式的预测温度用于热误差建模，替代实测的温度数据，该方法有利于更准确地定量描述机床热变形误差的时滞响应特性，进而提高机床热误差预测的精度和鲁棒性。

附图说明

图1示出一种XK26龙门移动镗铣床及其所在车间环境温度连续20天测量数据（测点间隔10分钟）及傅立叶变换得到的幅值谱；

图2示出环境温度时频特性及傅立叶级数分解的特征图；

图3示出XK26所处环境温度的在不同季节的实测及预测效果对比图。

具体实施方式

结合附图，以一台XK26龙门铣床的所处环境温度分解预测的实施过程，对本发明进行详细的说明。

附图1示出一种XK26龙门移动镗铣床及其所在车间环境温度连续20天测量数据（测点间隔10分钟）及傅立叶变换得到的幅值谱，其中傅立叶级数的均值项用时间序列的滑动平均数据替代，保留了温度信号的低频部分，基频和倍频部分信号通过不同频率、幅值、相位的谐波项叠加来替代。

1、环境温度温度测量

机床的环境温度测量***选用基于1-wire总线的DS18B20数字温度传感器，温度传感器布置于机床主轴前端附近。环境温度温度的记录与机床的开机和加工状态无关，因此可以连续监控，由于是长期测量，温度是缓变的且采集数据量大，因此测量时间间隔可设定稍长一些，在此设定10分钟测量间隔，即ΔT=10min，于是N=144，春夏秋冬四个季节分别进行各20天左右的温度测量和记录。从图1的连续20天测量数据及其傅立叶变换得到的幅值谱中可以看出，环境温度具有明显的日周期特性，其等效基频为1/day，该频率成分的幅值所占比重高，且存在幅值依次降低的倍频成分。而且频率低于基频的部分也比较明显，它体现出环境温度较长周期的波动成分。

2、环境温度模型辨识

经过计算，第四阶及以上的高频温度波动成分，机床的响应温度变化幅值低于0.065℃，可忽略不计，因此环境温度可用三阶傅立叶级数形式表达。

选用春夏秋冬各10天的数据进行拟合非线性最小二乘拟合。得到非线性参数拟合数据如表2：

表2温度预测模型参数辨识

将上述参数代入上述得到的总的温度预测模型 $t_e\left(x\right)=A_0^{\left(1\right)}\left(x\right)+\underset{n=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_n(T_{\max }\left(x\right)-A_0^{\left(1\right)}\left(x\right))\sin (n{\mathrm{\&omega;}}_{0}x+({\mathrm{\&phi;}}_{0n}+{\mathrm{\&alpha;}}_n\&CenterDot;A_0^{\left(1\right)}\left(x\right))),n=\mathrm{1,2,3},L,$

于是本实施例的环境温度预测模型为：

$\begin{array}{c}{t}_e\left(x\right)=A_0^{\left(1\right)}\left(x\right)+0.73(T_{\max }\left(x\right)-A_0^{\left(1\right)}\left(x\right))\sin (0.0436x-1.75+0.012A_0^{\left(1\right)}\left(x\right))\\ -0.23(T_{\max }\left(x\right)-A_0^{\left(1\right)}\left(x\right))\sin (0.0872x+0.45-0.008A_0^{\left(1\right)}\left(x\right))\\ -0.04(T_{\max }\left(x\right)-A_0^{\left(1\right)}\left(x\right))\sin (0.131x+0.77-0.007A_0^{\left(1\right)}\left(x\right))\end{array}---\left(6\right)$

3、温度预测效果验证

分别选用四个季节共30天的数据，对所辨识的模型式（6）进行预测效果对比。通过比较可以发现，在不同季节0-40℃的温度范围内，解析模型都能够很好的预测并复现当前温度，预测值与测量值之间的残差小于1℃，而且该预测模型是一种变系数的解析模型，能够方便的应用于机床的时滞和非线性热变形响应分析中。

以上所述，仅为本发明较佳的具体实施方式，但本发明保护的范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (1)

1.一种重型机床车间环境温度解析建模方法，其具体步骤是：

1)环境温度傅立叶级数分解

环境温度展开为傅立叶级数形式，即：

$t_e\left(x\right)=A_0+\underset{n=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}{A}_{n}sin({\mathrm{n\&omega;}}_{0}x+{\mathrm{\&phi;}}_n)---\left(1\right)$

其中，x表示时间，单位为分钟，它由时间参考起点和当前时刻确定；基频ω₀＝2π/T₀，A₀为温度的均值项，A_n是各个阶次的温度波的振幅，n＝1,2,…，φ_n为相位角，其物理意义为各阶滞后时间，下面对A₀,A_n,φ_n各项参数进行调整和替换；

2)均值项A₀

均值项A₀的物理意义是当前温度点所在的周期的均值，在此由历史温度的滑动平均值来替代,其由下式(2)得到：

$A_0=A_0^{\left(1\right)}\left(x\right)=\frac{t_{\&lsqb;x\&rsqb;-N}+t_{\&lsqb;x\&rsqb;-(N-1)}+...+t_{\&lsqb;x\&rsqb;-1}}{N},\&lsqb;x\&rsqb;\&GreaterEqual;N---\left(2\right)$

其中，[x]是对任意时刻x取整数，N是一个周期的采样点数，t_[x]是在时间点x的测量温度，环境温度是缓变信息，认为时间[x]≤x＜[x]+1期间的环境温度等同于t_[x]；在滑动平均的方法中，要求至少用一个周期N个采样点的历史数据来预测当前的温度；新的测量值t_[x]-1引入，旧的t_[x]-N退出，于是预测的均值A₀可以在线更新，滑动平均方法保留了所有低于基频频率的波动信息；

3)幅值项A_n

A_n为分解后各阶次温度波动的幅值项，其与所考虑周期内波动最大值T_max与平均值A₀之差成正比，也就是：

$A_n={\mathrm{\&beta;}}_n\left(T_{max}\right(x)-A_0^{\left(1\right)}(x\left)\right),n=1,2,3,...---\left(3\right)$

其中：

$T_{max}\left(x\right)=\underset{\&lsqb;x\&rsqb;-N<u<\&lsqb;x\&rsqb;-1}{max}\left(t_{\&lsqb;u\&rsqb;}\right)---\left(4\right)$

β_n含义为不同频率成分温度波对总的温度变化的贡献成分即权重，表示车间的固有特征；(T_max(x)-A₀ ⁽¹⁾(x))在不同的观测周期是变化的，但对已知的环境温度历史观测值是确定的；

4)相位φ_n

相位φ_n表示为：

${\mathrm{\&phi;}}_n={\mathrm{\&phi;}}_{0n}+{\mathrm{\&alpha;}}_n\&CenterDot;A_0^{\left(1\right)}\left(x\right),n=1,2,3,...---\left(5\right)$

在此α_n表示季节变换的滞后时间系数，代表了特定车间的固有热特性，其为特定值；

5)总的温度预测模型

综上所述，任意x时刻的环境温度预测的解析模型可表示为：

$t_e\left(x\right)=A_0^{\left(1\right)}\left(x\right)+\underset{n=1}{\overset{\mathrm{\&infin;}}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_n\left(T_{max}\right(x)-A_0^{\left(1\right)}(x\left)\right)sin({\mathrm{n\&omega;}}_{0}x+({\mathrm{\&phi;}}_{0n}+{\mathrm{\&alpha;}}_n\&CenterDot;A_0^{\left(1\right)}\left(x\right)\left)\right),n=1,2,3,...---\left(6\right)$

式(6)中，滑动平均值A₀ ⁽¹⁾(x)和最大值T_max(x)可以通过至少一个周期的历史测量值计算得到，而且是可以实时更新的，基频ω₀是已知的，因此可以由实测温度数据用非线性最小二乘原理辨识上述解析模型中β_n,φ_0n和α_n，进而实现环境温度的在线预测。