CN115409067A - 一种数控机床组件剩余寿命预测方法 - Google Patents

一种数控机床组件剩余寿命预测方法 Download PDF

Info

Publication number
CN115409067A
CN115409067A CN202211087850.0A CN202211087850A CN115409067A CN 115409067 A CN115409067 A CN 115409067A CN 202211087850 A CN202211087850 A CN 202211087850A CN 115409067 A CN115409067 A CN 115409067A
Authority
CN
China
Prior art keywords
machine tool
component
signal
formula
state monitoring
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202211087850.0A
Other languages
English (en)
Inventor
牟黎明
李丽娟
闫莉
于广伟
高晓兵
李雨菲
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Technological University
Original Assignee
Xian Technological University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Technological University filed Critical Xian Technological University
Priority to CN202211087850.0A priority Critical patent/CN115409067A/zh
Publication of CN115409067A publication Critical patent/CN115409067A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/20—Administration of product repair or maintenance

Landscapes

  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Engineering & Computer Science (AREA)
  • Economics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Strategic Management (AREA)
  • Tourism & Hospitality (AREA)
  • Physics & Mathematics (AREA)
  • General Business, Economics & Management (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Numerical Control (AREA)

Abstract

本发明涉及一种数控机床组件剩余寿命预测方法，包括：1、确定组件状态监测数据类型及监测部位，搭建组件状态监测平台，实时采集数控机床组件运行状态信息；2、对监测信号进行消除趋势项、降噪、特征提取以及PCA降维处理；3、根据数控机床组件状态监测信息，考虑机床运行工况，以主成分分析法降维后得到的新特征作为内部协变量，将运行载荷及速度作为外部协变量，建立一种充分考虑机床运行状态信息威布尔回归模型；4、基于威布尔回归模型，开展机床组件剩余寿命预测；本发明综合考虑机床加工过程中组件外部工况运行条件和内部振动特征，提高剩余寿命预测精度，准确预测剩余寿命，对于降低停机损失，提高生产效率，增加经济收益意义显著。

Description

一种数控机床组件剩余寿命预测方法
技术领域
本发明属于数控机床技术领域，涉及一种数控机床组件剩余寿命预测方法，具体涉及数控机床组件运行状态监测数据获取、数控机床组件状态监测信号处理、数控机床组件威布尔回归模型构建、基于威布尔回归模型的数控机床组件剩余寿命预测。
背景技术
在科技飞速发展的今天，剩余寿命预测技术俨然成为智能制造的重要组成部分，开展高效准确的剩余寿命预测，能够为制造装备预防性维护策略的制定提供指导，对于降低装备停机损失以及提升装备安全性和可靠性意义重大。
数控机床作为智能制造的“工作母机”，如若发生故障未能及时处理，那么势必会影响整个生产进度，轻则导致加工零部件报废、无法按期交付，项目合同终止，重则甚至会威胁到现场操作人员的人身安全。因此，开展机床组件剩余寿命预测研究，在机床组件发生故障之前，准确识别机床组件运行状态，对机床进行合理的维修保养，能够大大降低机床维护成本，提高机床利用率。
目前，关于数控机床的剩余寿命预测多是针对轴承、齿轮、刀具等组件展开，而对机床整机以及其它组件的剩余寿命预测研究相对较少；并且数控机床剩余寿命预测研究需要有大量、长周期的数据作为支撑，但随着机床结构功能复杂性的不断增加，加之变工况、多工况情形，致使数据收集有限，进而机床的剩余寿命预测也愈加困难；传统的数控机床剩余寿命预测模型多是基于机床组件单一工况运行条件展开，当工况发生变化时，无法及时更新预测模型，这就导致模型的适用性不强；或者仅基于机床组件的信号特征开展剩余寿命预测，未将机床实时加工的工况考虑进来，亦或是考虑了机床组件的运行工况，但未将二者进行很好的融合，这势必会导致机床组件剩余寿命预测结果存在较高的误差。
发明内容
针对传统的数控机床剩余寿命预测模型多是基于机床组件在某一单一工况下的运行展开，这就需要针对每种工况条件单独建模，繁琐又耗时；或者仅基于机床组件的信号特征开展剩余寿命预测，未将机床实时加工的工况考虑进来，亦或考虑了机床的运行工况，但未将二者进行融合，这势必会导致机床组件剩余寿命预测结果存在较高的误差。因此，本发明以数控机床组件为研究对象，对组件展开运行状态信息监测以及剩余寿命预测研究，以期探索更为合理准确的剩余寿命预测方法，为制造企业实际生产过程中机床组件的预测性维护策略制定以及机床组件健康管理提供一定的参考和指导。
为解决上述技术问题，本发明是采用如下技术方案实现的，结合附图说明如下：
一种数控机床组件剩余寿命预测方法，包括下述步骤：
步骤一、根据机床组件的基本结构以及工作原理，确定机床组件状态监测数据类型及监测部位，通过搭建组件状态监测平台，实时采集数控机床组件运行状态信息；
步骤二、对获得的组件状态监测信号进行消除趋势项、降噪、特征提取以及主成分分析(Principal Component Analysis，PCA)降维处理，获取趋势性更好的信号特征；
步骤三、根据数控机床组件状态监测信息，考虑机床的运行工况，以PCA降维后得到的信号特征作为内部协变量，将运行载荷及速度作为外部协变量，建立一种考虑机床运行状态信息的威布尔回归模型；
步骤四、基于威布尔回归模型的数控机床组件剩余寿命预测；
所述的步骤一中数控机床组件运行状态信息采集的具体方法如下：
(1)选择机床组件状态监测数据类型；
选取能够全面反映机床组件状态变化的信号进行监测；
(2)选择机床组件状态监测部位；
监测部位的选择不能影响机床正常的生产加工，并且易于传感器的安装；
(3)选择传感器；
对于高速旋转类组件，采用非接触式位移传感器监测回转精度误差，对于导轨、刀具这类组件使用加速度传感器采集振动信号；
(4)搭建机床组件状态监测平台；
搭建数控机床组件状态监测平台，主要包括被监测组件、传感器、数据采集仪及信号分析平台、连接线以及PC端；
(5)获取机床组件状态监测信息；
制定组件状态监测试验方案，基于搭建的机床组件状态监测平台，得到组件不同工况情形下的信号，并将信号测试分析系统采集到的信号，传输到信号分析平台，获取组件监测信息；
所述的步骤二中信号处理的具体方法如下：
(1)组件状态监测信号趋势项消除；
当出现原始特征信号在基线附近上下偏移并且偏移的程度随时间变化的现象时，需要根据采集到的信号特点，选择相应的趋势项消除方法，获取消除趋势项后的信号；
(2)组件状态监测信号降噪处理；
应用小波阈值降噪法对消除趋势项后的信号进行降噪处理；
对于传感器采集到的数控机床组件任意一段输出信号ω(t)，且该信号中含有噪声信号，则有：
ω(t)＝s(t)+σ(t)……………………………………(1)
式中，t为时间，s(t)为目标信号，或者是能够反映数控机床组件运行状态的有用信号；σ(t)则为噪声信号；
对于任意信号ω(t)，利用式(2)，得到它的连续小波变换Wf(ξ,δ)；
Figure BDA0003835903550000031
式中，ψξ,δ(t)为小波基函数；
Figure BDA0003835903550000032
为ψξ,δ(t)的共轭，利用式(3)得到ψξ,δ(t)的取值；
Figure BDA0003835903550000033
式中，ξ表示伸缩尺度参数，δ表示平移尺度参数，且有ξ,δ∈R,ξ≠0；
对式(2)、(3)进行小波变换离散化处理，取离散值ξ＝2m,δ＝2mn；m和n为自然数；利用式(4)得到任意信号的离散小波变换公式；
Figure BDA0003835903550000034
式中，W′f(m,n)即为ω(t)的离散小波变换；
根据机床组件信号规律，选择相应的小波基对组件原始信号进行小波分解降噪，进而得到降噪后的信号；
(3)组件状态监测信号特征提取；
对采集到的不同工况下的组件状态监测信号展开分析，根据组件状态监测信号特点，提取组件状态监测信号时域和频域特征；
1)组件状态监测信号时域特征提取；
数控机床组件状态监测信号的时域分析是基于时域幅值波形展开，机床加工运行过程中，不同工况情形下，组件信号各时域特征存在差异，各时域特征对不同工况的响应不同，利用式(5)-(8)提取时域特征：有效值Xrms、峰峰值Xv、偏度Xα和峭度XK；
Figure BDA0003835903550000035
Xv＝max(xi)-min(xi)……………………………(6)
Figure BDA0003835903550000041
Figure BDA0003835903550000042
式中，N为采样点数(1≤i≤N)，xi为一离散信号；
2)组件状态监测信号频域特征提取；
应用傅里叶变换将监测到的组件时域信号转换成频域信号，进而实现组件监测信号频域特征提取，利用式(9)-(11)提取频域特征：重心频率XFC、频率均方根XRMSF、频率标准差XRVF；
Figure BDA0003835903550000043
Figure BDA0003835903550000044
Figure BDA0003835903550000045
式中，A(k)为信号幅值，fk为频率，K表示谱线数目(1≤k≤K)；
3)组件状态监测信号特征降维；
考虑不同的特征对机床组件退化或故障的响应不同，应用主成分分析PCA对信号特征进行降维，去除相关度较低的特征，首先利用式(12)对g个特征样本Z＝[z1,z2,…,zg]进行信号特征去中心化处理，得到去中心化后的信号特征样本C；
Figure BDA0003835903550000046
式中，
Figure BDA0003835903550000047
为第g个样本均值，cg为第g个去中心化后的信号特征样本；
进而利用式(13)，得到特征样本集的协方差矩阵S；
Figure BDA0003835903550000051
式中，ci表示第i行元素，cj表示第j列元素，(ci)T为ci的转置；
假设存在等式(14)，即一个线性变换u使得uci的方差最大；
Figure BDA0003835903550000052
根据式(14)，为获得最大协方差，构建式(15)所示的拉格朗日条件极值；
Figure BDA0003835903550000053
则有等式(16)：
Figure BDA0003835903550000054
式中，λ为协方差矩阵S的特征值，u为特征向量；
在此基础上，求出全部特征值λ，其中λ1≥λ2≥…≥λg，及对应的特征向量u1,u2,…,ug，此时得到新的特征Y：
Y＝UT×C……………………………………(17)
式中，U＝(u1,u2,…,ug)。
定义第i个主分量yi的“方差贡献率”计算公式为：
Figure BDA0003835903550000055
则有前h个主分量y1,y2,…,yh的“累计方差贡献率”为：
Figure BDA0003835903550000056
当前h个主分量的累计方差贡献率大于95％时，取前h个主分量作为新特征，此时有：
Figure BDA0003835903550000057
其余(g-h)个新特征则被舍去；
至此在获得机床组件状态监测信息基础上，完成了信号处理，得到降维后的主成分PC；
本发明以PCA降维后的各主成分作为剩余寿命预测模型的内部协变量，将机床组件的负载和速度作为外部协变量，建立机床组件剩余寿命预测模型；
所述的步骤三中威布尔回归模型的构建具体方法如下：
机床组件的故障数据服从威布尔分布，威布尔分布模型的累积故障概率函数F(t)，如公式(21)所示：
Figure BDA0003835903550000061
式中，t是时间变量，t≥0；θ是尺度参数，θ＞0；γ是形状参数，γ＞0；
以机床组件运行过程中的内、外部协变量描述模型的尺度参数，将威布尔模型推广到威布尔回归模型WRM，此时，WRM模型的故障概率函数，如式(22)所示：
Figure BDA0003835903550000062
式中，R(t)表示可靠度函数；a0表示截距；I为内部协变量，q(1≤l≤q)表示内部协变量数目；E为外部协变量，p(1≤j≤p)表示外部协变量数目；a为内部协变量回归系数，al(1≤l≤q)表示第l个内部斜变量的回归系数，b为外部协变量回归系数，bj(1≤j≤p)表示第j个外部协变量的回归系数；
将式(22)推广，得到WRM模型故障概率函数F(t,I,E)为：
Figure BDA0003835903550000063
由此，得到WRM模型的故障率函数λ(t,I,E)，如公式(24)所示：
Figure BDA0003835903550000064
根据公式(23)和(24)，推出WRM模型的概率密度函数f(t,I,E)为：
Figure BDA0003835903550000065
至此完成了数控机床组件威布尔回归模型的建立；
所述的步骤四中基于威布尔回归模型的数控机床组件剩余寿命预测方法的具体步骤如下：
(1)明晰机床组件可靠度函数；
根据步骤三的威布尔回归模型的构建原理，得到机床组件在t时刻的可靠度函数为：
Figure BDA0003835903550000071
式中，a0为截距，γ为分布的形状参数，a1,a2,…,aq分别为内部协变量PC1,PC2,…,PCq的回归系数，b1和b2则分别为外部协变量负载L和速度S的回归系数；
(2)机床组件WRM模型参数估计；
应用极大似然估计法对WRM模型进行参数估计，根据公式(22)则有对数似然函数L(γ,θ′)如下：
Figure BDA0003835903550000072
式中，θ′＝exp(a0+a1×PC1+a2×PC2+…+aq×PCq+b1×L+b2×S)；
分别对式(27)中参数a0、a1,a2,…，aq、b1、b2求偏导，并将偏导数设为零，进而得到各参数的最优值，应用MATLAB编程软件中的fminsearch函数进行参数求解，进而得到不同工况下机床组件的WRM模型；
(3)机床组件剩余寿命预测；
数控机床组件的剩余寿命预测，是指机床关键组件从当前时刻t0运行到用户要求可靠性阈值或故障时刻t的时间差，即为该组件的剩余寿命，设剩余寿命函数为RUL(t,Z(t))，可靠度函数为R，t0时刻剩余寿命的几何意义为可靠度函数曲线在运行区间[t0,t]上的面积，则有：
Figure BDA0003835903550000073
式中，Z(t)代表机床组件运行工况信息；
因此，根据式(26)、(27)获取机床组件的可靠度，结合式(28)，实现机床组件在时间t0时的剩余寿命预测，采用留一法交叉验证展示WRM剩余寿命预测模型，将采集到的数据中的一组数据作为测试数据，剩余的全部数据用于WRM模型的参数估计，以此验证机床组件剩余寿命预测方法的准确性。
至此完成了基于威布尔回归模型的数控机床组件剩余寿命预测。
与现有技术相比本发明的有益效果是：
本发明的机床组件剩余寿命预测方法，以组件不同工况运行条件下状态监测信息为输入，综合考虑了机床加工过程中组件的工况运行条件和信号特征，应用PCA降维后的主成分作为内部协变量，将组件的负载和速度作为外部协变量，建立了威布尔回归模型，通过参数估计获得了模型最优参数值，进而实现机床组件剩余寿命预测，考虑机床组件的工况运行条件，提高剩余寿命预测精度，准确合理的剩余寿命预测结果，对于降低生产制造企业停机损失，提高企业生产效率，增加经济收益意义显著。
附图说明
图1是本发明所述的数控机床组件剩余寿命预测方法流程图；
图2是本发明的数控机床主轴监测点；
图3是本发明的数控机床组件状态监测平台示意图；
图4a是本发明的数控机床主轴加速度信号原始波形图；
图4b是本发明的数控机床主轴加速度信号消除趋势项后波形图；
图5a是本发明的小波降噪前的主轴振动信号波形图；
图5b是本发明的小波降噪后的主轴振动信号波形图；
图6a是本发明的数控机床主轴不同工况有效值变化趋势对比图；
图6b是本发明的数控机床主轴不同工况峰峰值变化趋势对比图；
图6c是本发明的数控机床主轴不同工况偏度变化趋势对比图；
图6d是本发明的数控机床主轴不同工况峭度变化趋势对比图；
图7a是本发明的数控机床主轴工况1各频域特征变化趋势对比图；
图7b是本发明的数控机床主轴工况2各频域特征变化趋势对比图；
图7c是本发明的数控机床主轴工况3各频域特征变化趋势对比图；
图7d是本发明的数控机床主轴工况4各频域特征变化趋势对比图；
图8是本发明的数控机床主轴各特征分量方差贡献率及累计方差贡献率；
图9是本发明的工况1运行条件下机床主轴剩余寿命预测结果；
图10a是本发明的工况2运行条件下机床主轴剩余寿命预测结果；
图10b是本发明的工况3运行条件下机床主轴剩余寿命预测结果；
图10c是本发明的工况4运行条件下机床主轴剩余寿命预测结果。
具体实施方式
下面结合附图对本发明作详细的描述：
参阅图1，本发明的数控机床组件剩余寿命预测方法包括下述步骤：根据机床组件的基本结构以及工作原理，确定组件状态监测数据类型及监测部位，搭建组件状态监测平台，采集数控机床组件运行状态信息；对采集到的信号进行消除趋势项、降噪、特征提取以及降维处理；威布尔回归模型构建；基于威布尔回归模型预测机床组件剩余寿命。
一、机床组件运行状态信息采集
1.确定机床组件状态监测数据类型及监测部位；
选取能够全面反映机床组件状态变化的信号进行监测；监测部位的选择不能影响机床正常的生产加工，并且易于传感器的安装；
2.传感器的选择；
对于高速旋转类组件，采用非接触式位移传感器监测回转精度误差，对于导轨、刀具这类组件使用加速度传感器采集振动信号；
3.搭建机床组件状态监测平台；
搭建数控机床组件状态监测平台，主要包括被监测组件、传感器、数据采集仪及信号分析平台、连接线以及PC端；
4.获取机床组件状态监测信息；
制定组件状态监测试验方案，基于搭建的机床组件状态监测平台，得到组件不同工况情形下的信号，并将信号测试分析系统采集到的信号，传输到信号分析平台，获取组件监测信息。
二、机床组件状态监测信号处理
1.组件状态监测信号趋势项消除；
当出现原始特征信号在基线附近上下偏移并且偏移的程度随时间变化的现象时，需要根据采集到的信号特点，选择相应的趋势项消除方法，获取消除趋势项后的信号；
2.组件状态监测信号降噪处理；
应用小波阈值降噪法对消除趋势项后的信号进行降噪处理；
对于传感器采集到的数控机床组件任意一段输出信号ω(t)，且该信号中含有噪声信号，则有ω(t)＝s(t)+σ(t)，其中t为时间，s(t)为目标信号，σ(t)则为噪声信号，对于任意信号ω(t)，它的连续小波变换表达式为
Figure BDA0003835903550000091
其中ψξ,δ(t)为小波基函数；
Figure BDA0003835903550000101
为ψξ,δ(t)的共轭，
Figure BDA0003835903550000102
其中ξ表示伸缩尺度参数，δ表示平移尺度参数，且有ξ,δ∈R,ξ≠0，由于采集到的机床组件信号为离散信号，因此，需要进行小波变换离散化处理，取离散值ξ＝2m,δ＝2mn(m和n为自然数)，利用式
Figure BDA0003835903550000103
得到任意信号的离散小波变换，其中W′f(m,n)为ω(t)的离散小波变换；然后根据机床组件信号规律，选择相应的小波基对组件原始信号进行小波分解降噪，进而得到降噪后的信号；
3.组件状态监测信号特征提取；
对采集到的不同工况下的组件状态监测信号展开分析，根据组件状态监测信号特点，提取组件状态监测信号时域特征和频域特征；
(1)组件状态监测信号时域特征提取；
数控机床组件状态监测信号的时域分析是基于时域幅值波形展开，机床加工运行过程中，不同工况情形下，组件信号各时域特征存在差异，各时域特征对不同工况的响应不同，利用式
Figure BDA0003835903550000104
提取有效值Xrms、利用式Xv＝max(xi)-min(xi)提取峰峰值Xv、利用式
Figure BDA0003835903550000105
提取偏度Xα、利用式
Figure BDA0003835903550000106
提取峭度XK，其中，N为采样点数(1≤i≤N)，xi为一离散信号；
(2)组件状态监测信号频域特征提取；
应用傅里叶变换将监测到的组件时域信号转换成频域信号，进而实现组件监测信号频域特征提取，利用式
Figure BDA0003835903550000107
提取频域特征重心频率XFC、利用式
Figure BDA0003835903550000108
提取频率均方根XRMSF、利用式
Figure BDA0003835903550000109
提取频率标准差XRVF；
(3)组件状态监测信号特征降维；
考虑不同的特征对机床组件退化或故障的响应不同，应用PCA对信号特征进行降维，去除相关度较低的特征，首先利用式
Figure BDA0003835903550000111
对g个特征样本Z＝[z1,z2,…,zg]进行信号特征去中心化处理，得到去中心化后的信号特征样本C；进而利用式
Figure BDA0003835903550000112
得到特征样本集的协方差矩阵S；假设存在等式
Figure BDA0003835903550000113
即一个线性变换u使得uci的方差最大；构建拉格朗日条件极值：
L(u,s,λ)＝uTSu+λ(1-|u|2)
s.t(|u|2＝1)
则有等式：
Figure BDA0003835903550000114
其中，λ为协方差矩阵S的特征值，u为特征向量；
在此基础上，求出全部特征值λ，其中λ1≥λ2≥…≥λg，及对应的特征向量u1,u2,…,ug，此时得到新的特征Y：Y＝UT×X，其中U＝(u1,u2,…,ug)；定义第i个主分量yi的“方差贡献率”计算公式
Figure BDA0003835903550000115
则有h个主分量y1,y2,…,yh的“累计方差贡献率”为
Figure BDA0003835903550000116
当前h个主分量的累计方差贡献率大于95％时，取前h个主分量作为新特征，此时有：
Figure BDA0003835903550000117
其余(g-h)个新特征则被舍去，进而完成信号的处理，得到降维后的主成分PC。
三、威布尔回归模型建立
机床组件的故障数据服从威布尔分布，威布尔分布模型的累积故障概率函数F(t)：
Figure BDA0003835903550000121
其中，t是时间变量，t≥0；θ是尺度参数，θ＞0；γ是形状参数，γ＞0；
以机床运行过程中的速度和负载作为外部协变量，以降维后的主成分PC作为内部协变量，将威布尔模型推广到威布尔回归模型WRM，此时，WRM模型的故障概率函数为：
Figure BDA0003835903550000122
式中，R(t)表示可靠度函数；a0表示截距；I为内部协变量，q(1≤l≤q)表示内部协变量数目；E为外部协变量，p(1≤j≤p)表示外部协变量数目；a为内部协变量回归系数，al(1≤l≤q)表示第l个内部斜变量的回归系数，b为外部协变量回归系数，bj(1≤j≤p)表示第j个外部协变量的回归系数；
由此，得到WRM模型的故障率函数λ(t,I,E)：
Figure BDA0003835903550000123
进而推出WRM模型的概率密度函数f(t,I,E)为：
Figure BDA0003835903550000124
四、基于威布尔回归模型的数控机床组件剩余寿命预测
1.机床组件可靠度函数；
根据构建的威布尔回归模型，得到机床组件在t时刻的可靠度函数R(t,PC,L,S)为：
Figure BDA0003835903550000125
式中，a0为截距，γ为分布的形状参数，a1,a2,…,aq分别为内部协变量PC1,PC2,…,PCq的回归系数，b1和b2则分别为外部协变量负载L和速度S的回归系数；
2.机床组件WRM模型参数估计；
应用极大似然估计法对WRM模型进行参数估计，根据WRM模型的故障率函数，则有对数似然函数：
Figure BDA0003835903550000131
式中，θ′＝exp(a0+a1×PC1+a2×PC2+…+aq×PCq+b1×L+b2×S)；
分别对上式中的参数a0、a1,a2,…，aq、b1、b2求偏导，并将偏导数设为零，进而得到各参数的最优值，应用MATLAB编程软件中的fminsearch函数进行参数求解，进而得到不同工况下机床组件的WRM模型；
3.机床组件剩余寿命预测；
数控机床组件的剩余寿命预测，是指机床组件从当前时刻t0运行到用户要求可靠性阈值或故障时刻t的时间差，即为该组件的剩余寿命。设剩余寿命函数为RUL(t,Z(t))，可靠度函数为R，t0时刻剩余寿命的几何意义为可靠度函数曲线在运行区间[t0,t]上的面积，则有
Figure BDA0003835903550000132
其中Z(t)代表机床组件运行工况信息；由此实现机床组件在时间t0时的剩余寿命预测。
采用留一法交叉验证展示WRM剩余寿命预测模型，将采集到的数据中的一组数据作为测试数据，剩余的全部数据用于WRM模型的参数估计，以此验证机床组件剩余寿命预测方法的准确性。
实施例
数控机床组件剩余寿命预测方法
本发明以数控机床关键组件主轴为例，进行实施例说明，以某型号数控机床主轴为研究对象，根据主轴的基本结构以及工作原理，选择振动信号作为主轴的状态监测指标，并确定主轴最佳监测点，参阅图2。
针对主轴这种高速旋转的组件，选择非接触式传感器，本发明选用5E101型电涡流位移传感器、1A313E型加速度传感器，监测数控机床主轴振动性能指标，选用DH5922D型动态信号测试分析系统采集主轴的状态信号，并将该系统与DHDAS信号分析平台相结合，通过连接传感器实现信号的实时采集与分析，搭建的组件状态监测平台示意图，参阅图3。
本发明仅以机床用户企业生产加工零件的一道工序为例，开展监测实验。加工过程中，采用B17R0.8刀具，对工件表面进行精加工，采用5E101型电涡流位移传感器、1A313E型加速度传感器测量机床主轴加工过程中的振动信号，采样频率为2kHz，为获取加工过程中机床主轴不同工况下的振动状态，以该工件正常加工工况为标准，合理的改变主轴转速，来获取不同转速下主轴的振动信号，具体的加工参数设置，如表1所示。
表1不同工况加工参数设置
Figure BDA0003835903550000141
根据采集到的信号特点，应用低通滤波法进行趋势项消除，消除趋势项前后的主轴加速度信号对比，参阅图4a、图4b。应用小波阈值降噪法对消除趋势项后的信号进行降噪处理，降噪前后的信号波形对比图，参阅图5a、图5b，通过对比可以发现，经过小波降噪处理后得到的信号在保留主轴原始振动信号基础上，去除了噪声信号的干扰，且降噪效果显著。
根据步骤二式(5)-(8)计算得到不同工况下主轴振动信号的各时域特征对比图，参阅图6a、图6b、图6c、图6d，由图6a、图6b、图6c、图6d可知，在机床的加工运行过程中，不同工况情形下，主轴的振动信号各时域特征存在差异，各时域特征对不同工况的响应不同。同理，根据步骤二中式(9)-(11)，能够得到不同工况下各频域特征变化趋势图，参阅图7a、图7b、图7c、图7d，由图7a、图7b、图7c、图7d可知，不同工况下，机床主轴振动信号频域特征存在差异，且随着时间的延长以及磨损的加剧，信号频域特征频率的幅值也会呈现一定的变化趋势。
根据步骤二中式(12)-(20)，对信号进行PCA特征降维，以工况1为例，计算得到该工况下各信号特征向量的特征值，如表2所示。
表2工况1振动信号特征向量的特征值
Figure BDA0003835903550000142
同时，计算各特征分量的方差贡献率与累计方差贡献率，得到前5个特征分量的累计方差贡献率为95.03％，参阅图8，因此，降维后的新特征向量能够包括原始特征向量的大部分信息。
根据步骤三的威布尔回归模型构建方法，结合式(21)-(25)，以主轴信号降维后的主成分作为内部协变量、以速度和负载作为外部协变量，构建主轴的威布尔回归模型。根据步骤四式(26)-(27)，应用MATLAB编程软件中的fminsearch函数进行参数求解，进而得到工况1运行条件下机床主轴WRM模型的可靠度函数：
Figure BDA0003835903550000151
θ′＝exp(0.054+0.060×PC1+0.036×PC2-10.220×PC3-1.357×PC4-19.420×PC5+2.417×L-0.095×S)
同理，可以获得不同工况下，机床组件主轴的WRM模型，相应的各参数值，如表3所示。
表3不同工况运行条件下机床主轴的WRM模型参数值
Figure BDA0003835903550000152
因此，根据步骤四式(28)，结合采集到的机床运行工况信息，获取任意时刻、任意工况条件下机床各组件的可靠度值，进而开展机床组件的剩余寿命预测研究。
本发明以工况1为例，将采集到的173组数据的前172组用于WRM模型的参数估计，剩余的1组数据作为测试数据，用于预测工况1运行条件下的剩余寿命，应用WRM模型预测机床主轴的剩余寿命，预测结果参阅图9。同理，可以得到不同工况运行条件下，机床主轴的剩余寿命预测结果，参阅图10a、图10b、图10c。通常20％范围内的剩余寿命预测误差，是可接受的，由图10a、图10b、图10c可以发现，应用WRM模型预测的机床组件剩余寿命，与组件实际剩余寿命非常接近，大部分落都在20％的误差范围内。
本发明根据数控机床组件的结构特征，选取了组件运行状态监测数据类型，并确定组件运行状态监测部位，搭建数控机床组件运行状态监测平台，获取了组件的状态监测数据；对获得的监测信号进行了趋势项消除、降噪以及特征提取等信号处理；考虑不同的特征对机床组件退化或故障的响应不同，应用PCA对信号特征进行降维，以降维后的主成分作为内部协变量，将组件的负载和速度作为外部协变量，建立了威布尔回归模型，通过参数估计获得了模型最优参数值，进而实现机床关键组件剩余寿命预测，准确合理的剩余寿命预测结果，对于降低生产制造企业停机损失，提高企业生产效率，增加经济收益意义显著。

Claims (8)

1.一种数控机床组件剩余寿命预测方法，包括下述步骤：
步骤一、确定机床组件状态监测数据类型及监测部位，通过搭建组件状态监测平台，实时采集数控机床组件运行状态信息；
步骤二、对获得的组件状态监测信号进行消除趋势项、降噪、特征提取以及主成分分析降维处理，获取信号特征；
步骤三、根据数控机床组件状态监测信息，考虑机床的运行工况，以主成分分析降维后得到的信号特征作为内部协变量，将运行载荷及速度作为外部协变量，建立一种考虑机床运行状态信息的机床组件威布尔回归模型；
步骤四、基于威布尔回归模型的数控机床组件剩余寿命预测。
2.根据权利要求1所述的数控机床组件剩余寿命预测方法，其特征在于：
所述的步骤一中实时采集数控机床组件运行状态信息，具体方法如下：
(1)选择机床组件状态监测数据类型；
选取能够全面反映机床组件状态变化的信号进行监测；
(2)选择机床组件状态监测部位；
监测部位的选择不能影响机床正常的生产加工，并且易于传感器的安装；
(3)选择传感器；
对于高速旋转类组件，采用非接触式位移传感器监测回转精度误差，对于导轨、刀具这类组件使用加速度传感器采集振动信号；
(4)搭建机床组件状态监测平台；
搭建数控机床组件状态监测平台，主要包括被监测组件、传感器、数据采集仪及信号分析平台、连接线以及PC端；
(5)获取机床组件状态监测信息；
制定组件状态监测试验方案，基于搭建的机床组件状态监测平台，得到组件不同工况情形下的信号，并将信号测试分析系统采集到的信号，传输到信号分析平台，获取组件监测信息。
3.根据权利要求1所述的数控机床组件剩余寿命预测方法，其特征在于：
所述组件状态监测信号进行消除趋势项处理，具体是指：
当出现原始特征信号在基线附近上下偏移并且偏移的程度随时间变化的现象时，需要根据采集到的信号特点，选择相应的趋势项消除方法，获取消除趋势项后的信号；
所述组件状态监测信号进行降噪处理，具体是指：
应用小波阈值降噪法对消除趋势项后的信号进行降噪处理；
对于传感器采集到的数控机床组件任意一段输出信号ω(t)，且该信号中含有噪声信号，则有：
ω(t)＝s(t)+σ(t)……………………………………(1)
式中，t为时间，s(t)为目标信号，或者是能够反映数控机床组件运行状态的有用信号；σ(t)则为噪声信号；
对于任意信号ω(t)，利用式(2)，得到它的连续小波变换Wf(ξ,δ)；
Figure FDA0003835903540000021
式中，ψξ,δ(t)为小波基函数；
Figure FDA0003835903540000022
为ψξ,δ(t)的共轭，利用式(3)得到ψξ,δ(t)的取值；
Figure FDA0003835903540000023
式中，ξ表示伸缩尺度参数，δ表示平移尺度参数，且有ξ,δ∈R,ξ≠0；
对式(2)、(3)进行小波变换离散化处理，取离散值ξ＝2m,δ＝2mn；m和n为自然数；利用式(4)得到任意信号的离散小波变换公式；
Figure FDA0003835903540000024
式中，Wf′(m,n)即为ω(t)的离散小波变换；
根据机床组件信号规律，选择相应的小波基对组件原始信号进行小波分解降噪，进而得到降噪后的信号；
所述组件状态监测信号进行特征提取，具体是指：
对采集到的不同工况下的组件状态监测信号展开分析，根据组件状态监测信号特点，提取组件状态监测信号时域和频域特征。
4.根据权利要求3所述的数控机床组件剩余寿命预测方法，其特征在于：
所述提取组件状态监测信号时域特征，具体内容包括：
数控机床组件状态监测信号的时域分析是基于时域幅值波形展开，机床加工运行过程中，不同工况情形下，组件信号各时域特征存在差异，各时域特征对不同工况的响应不同，利用式(5)-(8)提取时域特征：有效值Xrms、峰峰值Xv、偏度Xα和峭度XK；
Figure FDA0003835903540000025
Xv＝max(xi)-min(xi)……………………………(6)
Figure FDA0003835903540000031
Figure FDA0003835903540000032
式中，N为采样点数(1≤i≤N)，xi为一离散信号。
5.根据权利要求3所述的数控机床组件剩余寿命预测方法，其特征在于：
所述提取组件状态监测信号频域特征，具体内容包括：
应用傅里叶变换将监测到的组件时域信号转换成频域信号，进而实现组件监测信号频域特征提取，利用式(9)-(11)提取频域特征：重心频率XFC、频率均方根XRMSF、频率标准差XRVF；
Figure FDA0003835903540000033
Figure FDA0003835903540000034
Figure FDA0003835903540000035
式中，A(k)为信号幅值，fk为频率，K表示谱线数目(1≤k≤K)。
6.根据权利要求1所述的数控机床组件剩余寿命预测方法，其特征在于：
所述组件状态监测信号进行主成分分析降维处理，具体包括：
考虑不同的特征对机床组件退化或故障的响应不同，应用主成分分析PCA对信号特征进行降维，去除相关度较低的特征，首先利用式(12)对g个特征样本Z＝[z1,z2,…,zg]进行信号特征去中心化处理，得到去中心化后的信号特征样本C；
Figure FDA0003835903540000036
式中，
Figure FDA0003835903540000037
为第g个样本均值，cg为第g个去中心化后的信号特征样本；
进而利用式(13)，得到特征样本集的协方差矩阵S；
Figure FDA0003835903540000041
式中，ci表示第i行元素，cj表示第j列元素，(ci)T为ci的转置；
假设存在等式(14)，即一个线性变换u使得uci的方差最大；
Figure FDA0003835903540000042
根据式(14)，为获得最大协方差，构建式(15)所示的拉格朗日条件极值；
Figure FDA0003835903540000043
则有等式(16)：
Figure FDA0003835903540000044
式中，λ为协方差矩阵S的特征值，u为特征向量；
在此基础上，求出全部特征值λ，其中λ1≥λ2≥…≥λg，及对应的特征向量u1,u2,…,ug，此时得到新的特征Y：
Y＝UT×C……………………………………(17)
式中，U＝(u1,u2,…,ug)；
定义第i个主分量yi的“方差贡献率”计算公式为：
Figure FDA0003835903540000045
则有前h个主分量y1,y2,…,yh的“累计方差贡献率”为：
Figure FDA0003835903540000046
当前h个主分量的累计方差贡献率大于95％时，取前h个主分量作为新特征，此时有：
Figure FDA0003835903540000051
其余(g-h)个新特征则被舍去；
至此在获得机床组件状态监测信息基础上，完成了信号处理，得到降维后的主成分PC。
7.根据权利要求1所述的数控机床组件剩余寿命预测方法，其特征在于：
步骤三中建立机床组件威布尔回归模型的具体方法如下：
机床组件的故障数据服从威布尔分布，威布尔分布模型的累积故障概率函数F(t)，如公式(21)所示：
Figure FDA0003835903540000052
式中，t是时间变量，t≥0；θ是尺度参数，θ＞0；γ是形状参数，γ＞0；
以机床运行过程中的速度和负载作为外部协变量，以降维后的主成分PC作为内部协变量，将威布尔模型推广到威布尔回归模型WRM，此时，WRM模型的故障概率函数，如式(22)所示：
Figure FDA0003835903540000053
式中，R(t)表示可靠度函数；a0表示截距；I为内部协变量，q(1≤l≤q)表示内部协变量数目；E为外部协变量，p(1≤j≤p)表示外部协变量数目；a为内部协变量回归系数，al(1≤l≤q)表示第l个内部斜变量的回归系数，b为外部协变量回归系数，bj(1≤j≤p)表示第j个外部协变量的回归系数；
将式(22)推广，得到WRM模型故障概率函数F(t,I,E)为：
Figure FDA0003835903540000054
由此，得到WRM模型的故障率函数λ(t,I,E)，如公式(24)所示：
Figure FDA0003835903540000055
根据公式(23)和(24)，推出WRM模型的概率密度函数f(t,I,E)为：
Figure FDA0003835903540000061
至此完成了数控机床组件威布尔回归模型的建立。
8.根据权利要求7所述的数控机床组件剩余寿命预测方法，其特征在于：
所述的步骤四中基于威布尔回归模型的数控机床组件剩余寿命预测的具体方法如下：
(1)明晰机床组件可靠度函数；
根据步骤三的威布尔回归模型的构建原理，得到机床组件在t时刻的可靠度函数R(t,PC,L,S)为：
Figure FDA0003835903540000062
式中，a0为截距，γ为分布的形状参数，a1,a2,…,aq分别为内部协变量PC1,PC2,…,PCq的回归系数，b1和b2则分别为外部协变量负载L和速度S的回归系数；
(2)机床组件WRM模型参数估计；
应用极大似然估计法对WRM模型进行参数估计，根据公式(22)则有对数似然函数L(γ,θ′)如下：
Figure FDA0003835903540000063
式中，θ′＝exp(a0+a1×PC1+a2×PC2+…+aq×PCq+b1×L+b2×S)；
分别对式(27)中参数a0、a1,a2,…，aq、b1、b2求偏导，并将偏导数设为零，进而得到各参数的最优值，应用MATLAB编程软件中的fminsearch函数进行参数求解，进而得到不同工况下机床组件的WRM模型；
(3)机床组件剩余寿命预测；
数控机床组件的剩余寿命预测，是指机床关键组件从当前时刻t0运行到用户要求可靠性阈值或故障时刻t的时间差，即为该组件的剩余寿命，设剩余寿命函数为RUL(t,Z(t))，可靠度函数为R，t0时刻剩余寿命的几何意义为可靠度函数曲线在运行区间[t0,t]上的面积，则有：
Figure FDA0003835903540000071
式中，Z(t)代表机床组件运行工况信息；
因此，根据式(26)、(27)获取机床组件的可靠度，结合式(28)，实现机床组件在时间t0时的剩余寿命预测。采用留一法交叉验证展示WRM剩余寿命预测模型，将采集到的数据中的一组数据作为测试数据，剩余的全部数据用于WRM模型的参数估计，以此验证机床组件剩余寿命预测方法的准确性。
CN202211087850.0A 2022-09-07 2022-09-07 一种数控机床组件剩余寿命预测方法 Pending CN115409067A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211087850.0A CN115409067A (zh) 2022-09-07 2022-09-07 一种数控机床组件剩余寿命预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211087850.0A CN115409067A (zh) 2022-09-07 2022-09-07 一种数控机床组件剩余寿命预测方法

Publications (1)

Publication Number Publication Date
CN115409067A true CN115409067A (zh) 2022-11-29

Family

ID=84163990

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211087850.0A Pending CN115409067A (zh) 2022-09-07 2022-09-07 一种数控机床组件剩余寿命预测方法

Country Status (1)

Country Link
CN (1) CN115409067A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115713044A (zh) * 2023-01-09 2023-02-24 佰聆数据股份有限公司 一种多工况切换下的机电设备剩余寿命分析方法和装置
CN116187047A (zh) * 2023-02-16 2023-05-30 北京思维实创科技有限公司 设备寿命预测方法及相关装置
CN116880359A (zh) * 2023-09-07 2023-10-13 天津艺仕机床有限公司 一种可信数控系统的测试方法及系统

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115713044A (zh) * 2023-01-09 2023-02-24 佰聆数据股份有限公司 一种多工况切换下的机电设备剩余寿命分析方法和装置
CN115713044B (zh) * 2023-01-09 2023-04-25 佰聆数据股份有限公司 一种多工况切换下的机电设备剩余寿命分析方法和装置
CN116187047A (zh) * 2023-02-16 2023-05-30 北京思维实创科技有限公司 设备寿命预测方法及相关装置
CN116880359A (zh) * 2023-09-07 2023-10-13 天津艺仕机床有限公司 一种可信数控系统的测试方法及系统
CN116880359B (zh) * 2023-09-07 2023-11-10 天津艺仕机床有限公司 一种可信数控系统的测试方法及系统

Similar Documents

Publication Publication Date Title
Luo et al. Early fault detection of machine tools based on deep learning and dynamic identification
CN115409067A (zh) 一种数控机床组件剩余寿命预测方法
CN105718876B (zh) 一种滚珠丝杠健康状态的评估方法
Zhou et al. Tool wear condition monitoring in milling process based on current sensors
Wu Dynamic data system: a new modeling approach
US9508042B2 (en) Method for predicting machining quality of machine tool
CN106842922B (zh) 一种数控加工误差优化方法
Wang et al. A kMap optimized VMD-SVM model for milling chatter detection with an industrial robot
CN113188794B (zh) 一种基于改进pso-bp神经网络齿轮箱故障诊断方法及装置
Chen et al. A data-driven model for thermal error prediction considering thermoelasticity with gated recurrent unit attention
CN111460701B (zh) 一种故障诊断模型训练方法和装置
CN108629864A (zh) 一种基于振动的电主轴径向精度表征方法及其系统
Oberle et al. A use case to implement machine learning for life time prediction of manufacturing tools
CN112475410A (zh) 一种铣削温度与多元影响因子的关联分析系统及方法
CN114559298B (zh) 基于物理信息融合的刀具磨损监测方法
CN115741235A (zh) 基于五轴加工中心刀具的磨损预测与健康管理方法
CN116184928A (zh) 一种考虑刀具磨损的机床空间几何误差模型建模方法
Xie et al. Tool wear status recognition and prediction model of milling cutter based on deep learning
CN111289231B (zh) 基于不完全B-spline数据拟合的转子系统健康监测方法和系统
Hey et al. Sensor selection method to accurately model the thermal error in a spindle motor
Li et al. Tool wear classification in milling for varied cutting conditions: with emphasis on data pre-processing
CN116821828A (zh) 一种基于工业数据的多维时序预测方法
Pandhare et al. Ball screw health monitoring with inertial sensors
CN108415372B (zh) 精密机床热误差补偿方法
CN111241629A (zh) 基于数据驱动的飞机液压泵性能变化趋势智能预测方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination