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

1.本发明属于数控机床技术领域，涉及一种数控机床组件剩余寿命预测方法，具体涉及数控机床组件运行状态监测数据获取、数控机床组件状态监测信号处理、数控机床组件威布尔回归模型构建、基于威布尔回归模型的数控机床组件剩余寿命预测。


背景技术：

2.在科技飞速发展的今天，剩余寿命预测技术俨然成为智能制造的重要组成部分，开展高效准确的剩余寿命预测，能够为制造装备预防性维护策略的制定提供指导，对于降低装备停机损失以及提升装备安全性和可靠性意义重大。
3.数控机床作为智能制造的“工作母机”，如若发生故障未能及时处理，那么势必会影响整个生产进度，轻则导致加工零部件报废、无法按期交付，项目合同终止，重则甚至会威胁到现场操作人员的人身安全。因此，开展机床组件剩余寿命预测研究，在机床组件发生故障之前，准确识别机床组件运行状态，对机床进行合理的维修保养，能够大大降低机床维护成本，提高机床利用率。
4.目前，关于数控机床的剩余寿命预测多是针对轴承、齿轮、刀具等组件展开，而对机床整机以及其它组件的剩余寿命预测研究相对较少；并且数控机床剩余寿命预测研究需要有大量、长周期的数据作为支撑，但随着机床结构功能复杂性的不断增加，加之变工况、多工况情形，致使数据收集有限，进而机床的剩余寿命预测也愈加困难；传统的数控机床剩余寿命预测模型多是基于机床组件单一工况运行条件展开，当工况发生变化时，无法及时更新预测模型，这就导致模型的适用性不强；或者仅基于机床组件的信号特征开展剩余寿命预测，未将机床实时加工的工况考虑进来，亦或是考虑了机床组件的运行工况，但未将二者进行很好的融合，这势必会导致机床组件剩余寿命预测结果存在较高的误差。


技术实现要素：

5.针对传统的数控机床剩余寿命预测模型多是基于机床组件在某一单一工况下的运行展开，这就需要针对每种工况条件单独建模，繁琐又耗时；或者仅基于机床组件的信号特征开展剩余寿命预测，未将机床实时加工的工况考虑进来，亦或考虑了机床的运行工况，但未将二者进行融合，这势必会导致机床组件剩余寿命预测结果存在较高的误差。因此，本发明以数控机床组件为研究对象，对组件展开运行状态信息监测以及剩余寿命预测研究，以期探索更为合理准确的剩余寿命预测方法，为制造企业实际生产过程中机床组件的预测性维护策略制定以及机床组件健康管理提供一定的参考和指导。
6.为解决上述技术问题，本发明是采用如下技术方案实现的，结合附图说明如下：
7.一种数控机床组件剩余寿命预测方法，包括下述步骤：
8.步骤一、根据机床组件的基本结构以及工作原理，确定机床组件状态监测数据类型及监测部位，通过搭建组件状态监测平台，实时采集数控机床组件运行状态信息；
9.步骤二、对获得的组件状态监测信号进行消除趋势项、降噪、特征提取以及主成分
分析(principal component analysis，pca)降维处理，获取趋势性更好的信号特征；
10.步骤三、根据数控机床组件状态监测信息，考虑机床的运行工况，以pca降维后得到的信号特征作为内部协变量，将运行载荷及速度作为外部协变量，建立一种考虑机床运行状态信息的威布尔回归模型；
11.步骤四、基于威布尔回归模型的数控机床组件剩余寿命预测；
12.所述的步骤一中数控机床组件运行状态信息采集的具体方法如下：
13.(1)选择机床组件状态监测数据类型；
14.选取能够全面反映机床组件状态变化的信号进行监测；
15.(2)选择机床组件状态监测部位；
16.监测部位的选择不能影响机床正常的生产加工，并且易于传感器的安装；
17.(3)选择传感器；
18.对于高速旋转类组件，采用非接触式位移传感器监测回转精度误差，对于导轨、刀具这类组件使用加速度传感器采集振动信号；
19.(4)搭建机床组件状态监测平台；
20.搭建数控机床组件状态监测平台，主要包括被监测组件、传感器、数据采集仪及信号分析平台、连接线以及pc端；
21.(5)获取机床组件状态监测信息；
22.制定组件状态监测试验方案，基于搭建的机床组件状态监测平台，得到组件不同工况情形下的信号，并将信号测试分析系统采集到的信号，传输到信号分析平台，获取组件监测信息；
23.所述的步骤二中信号处理的具体方法如下：
24.(1)组件状态监测信号趋势项消除；
25.当出现原始特征信号在基线附近上下偏移并且偏移的程度随时间变化的现象时，需要根据采集到的信号特点，选择相应的趋势项消除方法，获取消除趋势项后的信号；
26.(2)组件状态监测信号降噪处理；
27.应用小波阈值降噪法对消除趋势项后的信号进行降噪处理；
28.对于传感器采集到的数控机床组件任意一段输出信号ω(t)，且该信号中含有噪声信号，则有：
29.ω(t)＝s(t) σ(t)
……………………………………
(1)
30.式中，t为时间，s(t)为目标信号，或者是能够反映数控机床组件运行状态的有用信号；σ(t)则为噪声信号；
31.对于任意信号ω(t)，利用式(2)，得到它的连续小波变换wf(ξ,δ)；
[0032][0033]
式中，ψ
ξ,δ
(t)为小波基函数；为ψ
ξ,δ
(t)的共轭，利用式(3)得到ψ
ξ,δ
(t)的取值；
[0034]
[0035]
式中，ξ表示伸缩尺度参数，δ表示平移尺度参数，且有ξ,δ∈r,ξ≠0；
[0036]
对式(2)、(3)进行小波变换离散化处理，取离散值ξ＝2m,δ＝2mn；m和n为自然数；利用式(4)得到任意信号的离散小波变换公式；
[0037][0038]
式中，w
′f(m,n)即为ω(t)的离散小波变换；
[0039]
根据机床组件信号规律，选择相应的小波基对组件原始信号进行小波分解降噪，进而得到降噪后的信号；
[0040]
(3)组件状态监测信号特征提取；
[0041]
对采集到的不同工况下的组件状态监测信号展开分析，根据组件状态监测信号特点，提取组件状态监测信号时域和频域特征；
[0042]
1)组件状态监测信号时域特征提取；
[0043]
数控机床组件状态监测信号的时域分析是基于时域幅值波形展开，机床加工运行过程中，不同工况情形下，组件信号各时域特征存在差异，各时域特征对不同工况的响应不同，利用式(5)-(8)提取时域特征：有效值x
rms
、峰峰值xv、偏度x
α
和峭度xk；
[0044][0045]
xv＝max(xi)-min(xi)
……………………………
(6)
[0046][0047][0048]
式中，n为采样点数(1≤i≤n)，xi为一离散信号；
[0049]
2)组件状态监测信号频域特征提取；
[0050]
应用傅里叶变换将监测到的组件时域信号转换成频域信号，进而实现组件监测信号频域特征提取，利用式(9)-(11)提取频域特征：重心频率x
fc
、频率均方根x
rmsf
、频率标准差x
rvf
；
[0051][0052]
[0053][0054]
式中，a(k)为信号幅值，fk为频率，k表示谱线数目(1≤k≤k)；
[0055]
3)组件状态监测信号特征降维；
[0056]
考虑不同的特征对机床组件退化或故障的响应不同，应用主成分分析pca对信号特征进行降维，去除相关度较低的特征，首先利用式(12)对g个特征样本z＝[z1,z2,
…
,zg]进行信号特征去中心化处理，得到去中心化后的信号特征样本c；
[0057][0058]
式中，为第g个样本均值，cg为第g个去中心化后的信号特征样本；
[0059]
进而利用式(13)，得到特征样本集的协方差矩阵s；
[0060][0061]
式中，ci表示第i行元素，cj表示第j列元素，(ci)
t
为ci的转置；
[0062]
假设存在等式(14)，即一个线性变换u使得uci的方差最大；
[0063][0064]
根据式(14)，为获得最大协方差，构建式(15)所示的拉格朗日条件极值；
[0065][0066]
则有等式(16)：
[0067][0068]
式中，λ为协方差矩阵s的特征值，u为特征向量；
[0069]
在此基础上，求出全部特征值λ，其中λ1≥λ2≥
…
≥λg，及对应的特征向量u1,u2,
…
,ug，此时得到新的特征y：
[0070]
y＝u
t
×c……………………………………
(17)
[0071]
式中，u＝(u1,u2,
…
,ug)。
[0072]
定义第i个主分量yi的“方差贡献率”计算公式为：
[0073][0074]
则有前h个主分量y1,y2,
…
,yh的“累计方差贡献率”为：
[0075][0076]
当前h个主分量的累计方差贡献率大于95％时，取前h个主分量作为新特征，此时有：
[0077][0078]
其余(g-h)个新特征则被舍去；
[0079]
至此在获得机床组件状态监测信息基础上，完成了信号处理，得到降维后的主成分pc；
[0080]
本发明以pca降维后的各主成分作为剩余寿命预测模型的内部协变量，将机床组件的负载和速度作为外部协变量，建立机床组件剩余寿命预测模型；
[0081]
所述的步骤三中威布尔回归模型的构建具体方法如下：
[0082]
机床组件的故障数据服从威布尔分布，威布尔分布模型的累积故障概率函数f(t)，如公式(21)所示：
[0083][0084]
式中，t是时间变量，t≥0；θ是尺度参数，θ＞0；γ是形状参数，γ＞0；
[0085]
以机床组件运行过程中的内、外部协变量描述模型的尺度参数，将威布尔模型推广到威布尔回归模型wrm，此时，wrm模型的故障概率函数，如式(22)所示：
[0086][0087]
式中，r(t)表示可靠度函数；a0表示截距；i为内部协变量，q(1≤l≤q)表示内部协变量数目；e为外部协变量，p(1≤j≤p)表示外部协变量数目；a为内部协变量回归系数，a
l
(1≤l≤q)表示第l个内部斜变量的回归系数，b为外部协变量回归系数，bj(1≤j≤p)表示第j个外部协变量的回归系数；
[0088]
将式(22)推广，得到wrm模型故障概率函数f(t,i,e)为：
[0089][0090]
由此，得到wrm模型的故障率函数λ(t,i,e)，如公式(24)所示：
[0091][0092]
根据公式(23)和(24)，推出wrm模型的概率密度函数f(t,i,e)为：
[0093][0094]
至此完成了数控机床组件威布尔回归模型的建立；
[0095]
所述的步骤四中基于威布尔回归模型的数控机床组件剩余寿命预测方法的具体步骤如下：
[0096]
(1)明晰机床组件可靠度函数；
[0097]
根据步骤三的威布尔回归模型的构建原理，得到机床组件在t时刻的可靠度函数为：
[0098][0099]
式中，a0为截距，γ为分布的形状参数，a1,a2,
…
,aq分别为内部协变量pc1,pc2,
…
,pcq的回归系数，b1和b2则分别为外部协变量负载l和速度s的回归系数；
[0100]
(2)机床组件wrm模型参数估计；
[0101]
应用极大似然估计法对wrm模型进行参数估计，根据公式(22)则有对数似然函数l(γ,θ
′
)如下：
[0102][0103]
式中，θ
′
＝exp(a0 a1×
pc1 a2×
pc2
…
aq×
pcq b1×
l b2×
s)；
[0104]
分别对式(27)中参数a0、a1,a2,
…
，aq、b1、b2求偏导，并将偏导数设为零，进而得到各参数的最优值，应用matlab编程软件中的fminsearch函数进行参数求解，进而得到不同工况下机床组件的wrm模型；
[0105]
(3)机床组件剩余寿命预测；
[0106]
数控机床组件的剩余寿命预测，是指机床关键组件从当前时刻t0运行到用户要求可靠性阈值或故障时刻t的时间差，即为该组件的剩余寿命，设剩余寿命函数为rul(t,z(t))，可靠度函数为r，t0时刻剩余寿命的几何意义为可靠度函数曲线在运行区间[t0,t]上的面积，则有：
[0107][0108]
式中，z(t)代表机床组件运行工况信息；
[0109]
因此，根据式(26)、(27)获取机床组件的可靠度，结合式(28)，实现机床组件在时间t0时的剩余寿命预测，采用留一法交叉验证展示wrm剩余寿命预测模型，将采集到的数据中的一组数据作为测试数据，剩余的全部数据用于wrm模型的参数估计，以此验证机床组件剩余寿命预测方法的准确性。
[0110]
至此完成了基于威布尔回归模型的数控机床组件剩余寿命预测。
[0111]
与现有技术相比本发明的有益效果是：
[0112]
本发明的机床组件剩余寿命预测方法，以组件不同工况运行条件下状态监测信息
为输入，综合考虑了机床加工过程中组件的工况运行条件和信号特征，应用pca降维后的主成分作为内部协变量，将组件的负载和速度作为外部协变量，建立了威布尔回归模型，通过参数估计获得了模型最优参数值，进而实现机床组件剩余寿命预测，考虑机床组件的工况运行条件，提高剩余寿命预测精度，准确合理的剩余寿命预测结果，对于降低生产制造企业停机损失，提高企业生产效率，增加经济收益意义显著。
附图说明
[0113]
图1是本发明所述的数控机床组件剩余寿命预测方法流程图；
[0114]
图2是本发明的数控机床主轴监测点；
[0115]
图3是本发明的数控机床组件状态监测平台示意图；
[0116]
图4a是本发明的数控机床主轴加速度信号原始波形图；
[0117]
图4b是本发明的数控机床主轴加速度信号消除趋势项后波形图；
[0118]
图5a是本发明的小波降噪前的主轴振动信号波形图；
[0119]
图5b是本发明的小波降噪后的主轴振动信号波形图；
[0120]
图6a是本发明的数控机床主轴不同工况有效值变化趋势对比图；
[0121]
图6b是本发明的数控机床主轴不同工况峰峰值变化趋势对比图；
[0122]
图6c是本发明的数控机床主轴不同工况偏度变化趋势对比图；
[0123]
图6d是本发明的数控机床主轴不同工况峭度变化趋势对比图；
[0124]
图7a是本发明的数控机床主轴工况1各频域特征变化趋势对比图；
[0125]
图7b是本发明的数控机床主轴工况2各频域特征变化趋势对比图；
[0126]
图7c是本发明的数控机床主轴工况3各频域特征变化趋势对比图；
[0127]
图7d是本发明的数控机床主轴工况4各频域特征变化趋势对比图；
[0128]
图8是本发明的数控机床主轴各特征分量方差贡献率及累计方差贡献率；
[0129]
图9是本发明的工况1运行条件下机床主轴剩余寿命预测结果；
[0130]
图10a是本发明的工况2运行条件下机床主轴剩余寿命预测结果；
[0131]
图10b是本发明的工况3运行条件下机床主轴剩余寿命预测结果；
[0132]
图10c是本发明的工况4运行条件下机床主轴剩余寿命预测结果。
具体实施方式
[0133]
下面结合附图对本发明作详细的描述：
[0134]
参阅图1，本发明的数控机床组件剩余寿命预测方法包括下述步骤：根据机床组件的基本结构以及工作原理，确定组件状态监测数据类型及监测部位，搭建组件状态监测平台，采集数控机床组件运行状态信息；对采集到的信号进行消除趋势项、降噪、特征提取以及降维处理；威布尔回归模型构建；基于威布尔回归模型预测机床组件剩余寿命。
[0135]
一、机床组件运行状态信息采集
[0136]
1.确定机床组件状态监测数据类型及监测部位；
[0137]
选取能够全面反映机床组件状态变化的信号进行监测；监测部位的选择不能影响机床正常的生产加工，并且易于传感器的安装；
[0138]
2.传感器的选择；
[0139]
对于高速旋转类组件，采用非接触式位移传感器监测回转精度误差，对于导轨、刀具这类组件使用加速度传感器采集振动信号；
[0140]
3.搭建机床组件状态监测平台；
[0141]
搭建数控机床组件状态监测平台，主要包括被监测组件、传感器、数据采集仪及信号分析平台、连接线以及pc端；
[0142]
4.获取机床组件状态监测信息；
[0143]
制定组件状态监测试验方案，基于搭建的机床组件状态监测平台，得到组件不同工况情形下的信号，并将信号测试分析系统采集到的信号，传输到信号分析平台，获取组件监测信息。
[0144]
二、机床组件状态监测信号处理
[0145]
1.组件状态监测信号趋势项消除；
[0146]
当出现原始特征信号在基线附近上下偏移并且偏移的程度随时间变化的现象时，需要根据采集到的信号特点，选择相应的趋势项消除方法，获取消除趋势项后的信号；
[0147]
2.组件状态监测信号降噪处理；
[0148]
应用小波阈值降噪法对消除趋势项后的信号进行降噪处理；
[0149]
对于传感器采集到的数控机床组件任意一段输出信号ω(t)，且该信号中含有噪声信号，则有ω(t)＝s(t) σ(t)，其中t为时间，s(t)为目标信号，σ(t)则为噪声信号，对于任意信号ω(t)，它的连续小波变换表达式为其中ψ
ξ,δ
(t)为小波基函数；为ψ
ξ,δ
(t)的共轭，其中ξ表示伸缩尺度参数，δ表示平移尺度参数，且有ξ,δ∈r,ξ≠0，由于采集到的机床组件信号为离散信号，因此，需要进行小波变换离散化处理，取离散值ξ＝2m,δ＝2mn(m和n为自然数)，利用式得到任意信号的离散小波变换，其中w
′f(m,n)为ω(t)的离散小波变换；然后根据机床组件信号规律，选择相应的小波基对组件原始信号进行小波分解降噪，进而得到降噪后的信号；
[0150]
3.组件状态监测信号特征提取；
[0151]
对采集到的不同工况下的组件状态监测信号展开分析，根据组件状态监测信号特点，提取组件状态监测信号时域特征和频域特征；
[0152]
(1)组件状态监测信号时域特征提取；
[0153]
数控机床组件状态监测信号的时域分析是基于时域幅值波形展开，机床加工运行过程中，不同工况情形下，组件信号各时域特征存在差异，各时域特征对不同工况的响应不同，利用式提取有效值x
rms
、利用式xv＝max(xi)-min(xi)提取峰峰值xv、利用式提取偏度x
α
、利用式提取峭度xk，其中，n为采样点
数(1≤i≤n)，xi为一离散信号；
[0154]
(2)组件状态监测信号频域特征提取；
[0155]
应用傅里叶变换将监测到的组件时域信号转换成频域信号，进而实现组件监测信号频域特征提取，利用式提取频域特征重心频率x
fc
、利用式提取频率均方根x
rmsf
、利用式提取频率标准差x
rvf
；
[0156]
(3)组件状态监测信号特征降维；
[0157]
考虑不同的特征对机床组件退化或故障的响应不同，应用pca对信号特征进行降维，去除相关度较低的特征，首先利用式对g个特征样本z＝[z1,z2,
…
,zg]进行信号特征去中心化处理，得到去中心化后的信号特征样本c；进而利用式得到特征样本集的协方差矩阵s；假设存在等式即一个线性变换u使得uci的方差最大；构建拉格朗日条件极值：
[0158]
l(u,s,λ)＝u
t
su λ(1-|u|2)
[0159]
s.t(|u|2＝1)
[0160]
则有等式：
[0161][0162]
其中，λ为协方差矩阵s的特征值，u为特征向量；
[0163]
在此基础上，求出全部特征值λ，其中λ1≥λ2≥
…
≥λg，及对应的特征向量u1,u2,
…
,ug，此时得到新的特征y：y＝u
t
×
x，其中u＝(u1,u2,
…
,ug)；定义第i个主分量yi的“方差贡献率”计算公式则有h个主分量y1,y2,
…
,yh的“累计方差贡献率”为当前h个主分量的累计方差贡献率大于95％时，取前h个主分量作为新特征，此时有：
[0164][0165]
其余(g-h)个新特征则被舍去，进而完成信号的处理，得到降维后的主成分pc。
[0166]
三、威布尔回归模型建立
[0167]
机床组件的故障数据服从威布尔分布，威布尔分布模型的累积故障概率函数f(t)：
[0168][0169]
其中，t是时间变量，t≥0；θ是尺度参数，θ＞0；γ是形状参数，γ＞0；
[0170]
以机床运行过程中的速度和负载作为外部协变量，以降维后的主成分pc作为内部协变量，将威布尔模型推广到威布尔回归模型wrm，此时，wrm模型的故障概率函数为：
[0171][0172]
式中，r(t)表示可靠度函数；a0表示截距；i为内部协变量，q(1≤l≤q)表示内部协变量数目；e为外部协变量，p(1≤j≤p)表示外部协变量数目；a为内部协变量回归系数，a
l
(1≤l≤q)表示第l个内部斜变量的回归系数，b为外部协变量回归系数，bj(1≤j≤p)表示第j个外部协变量的回归系数；
[0173]
由此，得到wrm模型的故障率函数λ(t,i,e)：
[0174][0175]
进而推出wrm模型的概率密度函数f(t,i,e)为：
[0176][0177]
四、基于威布尔回归模型的数控机床组件剩余寿命预测
[0178]
1.机床组件可靠度函数；
[0179]
根据构建的威布尔回归模型，得到机床组件在t时刻的可靠度函数r(t,pc,l,s)为：
[0180][0181]
式中，a0为截距，γ为分布的形状参数，a1,a2,
…
,aq分别为内部协变量pc1,pc2,
…
,pcq的回归系数，b1和b2则分别为外部协变量负载l和速度s的回归系数；
[0182]
2.机床组件wrm模型参数估计；
[0183]
应用极大似然估计法对wrm模型进行参数估计，根据wrm模型的故障率函数，则有对数似然函数：
[0184][0185]
式中，θ
′
＝exp(a0 a1×
pc1 a2×
pc2
…
aq×
pcq b1×
l b2×
s)；
[0186]
分别对上式中的参数a0、a1,a2,
…
，aq、b1、b2求偏导，并将偏导数设为零，进而得到各参数的最优值，应用matlab编程软件中的fminsearch函数进行参数求解，进而得到不同工况下机床组件的wrm模型；
[0187]
3.机床组件剩余寿命预测；
[0188]
数控机床组件的剩余寿命预测，是指机床组件从当前时刻t0运行到用户要求可靠性阈值或故障时刻t的时间差，即为该组件的剩余寿命。设剩余寿命函数为rul(t,z(t))，可靠度函数为r，t0时刻剩余寿命的几何意义为可靠度函数曲线在运行区间[t0,t]上的面积，则有其中z(t)代表机床组件运行工况信息；由此实现机床组件在时间t0时的剩余寿命预测。
[0189]
采用留一法交叉验证展示wrm剩余寿命预测模型，将采集到的数据中的一组数据作为测试数据，剩余的全部数据用于wrm模型的参数估计，以此验证机床组件剩余寿命预测方法的准确性。
[0190]
实施例
[0191]
数控机床组件剩余寿命预测方法
[0192]
本发明以数控机床关键组件主轴为例，进行实施例说明，以某型号数控机床主轴为研究对象，根据主轴的基本结构以及工作原理，选择振动信号作为主轴的状态监测指标，并确定主轴最佳监测点，参阅图2。
[0193]
针对主轴这种高速旋转的组件，选择非接触式传感器，本发明选用5e101型电涡流位移传感器、1a313e型加速度传感器，监测数控机床主轴振动性能指标，选用dh5922d型动态信号测试分析系统采集主轴的状态信号，并将该系统与dhdas信号分析平台相结合，通过连接传感器实现信号的实时采集与分析，搭建的组件状态监测平台示意图，参阅图3。
[0194]
本发明仅以机床用户企业生产加工零件的一道工序为例，开展监测实验。加工过程中，采用b17r0.8刀具，对工件表面进行精加工，采用5e101型电涡流位移传感器、1a313e型加速度传感器测量机床主轴加工过程中的振动信号，采样频率为2khz，为获取加工过程中机床主轴不同工况下的振动状态，以该工件正常加工工况为标准，合理的改变主轴转速，来获取不同转速下主轴的振动信号，具体的加工参数设置，如表1所示。
[0195]
表1不同工况加工参数设置
[0196][0197]
根据采集到的信号特点，应用低通滤波法进行趋势项消除，消除趋势项前后的主
轴加速度信号对比，参阅图4a、图4b。应用小波阈值降噪法对消除趋势项后的信号进行降噪处理，降噪前后的信号波形对比图，参阅图5a、图5b，通过对比可以发现，经过小波降噪处理后得到的信号在保留主轴原始振动信号基础上，去除了噪声信号的干扰，且降噪效果显著。
[0198]
根据步骤二式(5)-(8)计算得到不同工况下主轴振动信号的各时域特征对比图，参阅图6a、图6b、图6c、图6d，由图6a、图6b、图6c、图6d可知，在机床的加工运行过程中，不同工况情形下，主轴的振动信号各时域特征存在差异，各时域特征对不同工况的响应不同。同理，根据步骤二中式(9)-(11)，能够得到不同工况下各频域特征变化趋势图，参阅图7a、图7b、图7c、图7d，由图7a、图7b、图7c、图7d可知，不同工况下，机床主轴振动信号频域特征存在差异，且随着时间的延长以及磨损的加剧，信号频域特征频率的幅值也会呈现一定的变化趋势。
[0199]
根据步骤二中式(12)-(20)，对信号进行pca特征降维，以工况1为例，计算得到该工况下各信号特征向量的特征值，如表2所示。
[0200]
表2工况1振动信号特征向量的特征值
[0201][0202]
同时，计算各特征分量的方差贡献率与累计方差贡献率，得到前5个特征分量的累计方差贡献率为95.03％，参阅图8，因此，降维后的新特征向量能够包括原始特征向量的大部分信息。
[0203]
根据步骤三的威布尔回归模型构建方法，结合式(21)-(25)，以主轴信号降维后的主成分作为内部协变量、以速度和负载作为外部协变量，构建主轴的威布尔回归模型。根据步骤四式(26)-(27)，应用matlab编程软件中的fminsearch函数进行参数求解，进而得到工况1运行条件下机床主轴wrm模型的可靠度函数：
[0204][0205]
θ
′
＝exp(0.054 0.060
×
pc1 0.036
×
pc
2-10.220
×
pc
3-1.357
×
pc
4-19.420
×
pc5 2.417
×
l-0.095
×
s)
[0206]
同理，可以获得不同工况下，机床组件主轴的wrm模型，相应的各参数值，如表3所示。
[0207]
表3不同工况运行条件下机床主轴的wrm模型参数值
[0208][0209]
因此，根据步骤四式(28)，结合采集到的机床运行工况信息，获取任意时刻、任意工况条件下机床各组件的可靠度值，进而开展机床组件的剩余寿命预测研究。
[0210]
本发明以工况1为例，将采集到的173组数据的前172组用于wrm模型的参数估计，
剩余的1组数据作为测试数据，用于预测工况1运行条件下的剩余寿命，应用wrm模型预测机床主轴的剩余寿命，预测结果参阅图9。同理，可以得到不同工况运行条件下，机床主轴的剩余寿命预测结果，参阅图10a、图10b、图10c。通常20％范围内的剩余寿命预测误差，是可接受的，由图10a、图10b、图10c可以发现，应用wrm模型预测的机床组件剩余寿命，与组件实际剩余寿命非常接近，大部分落都在20％的误差范围内。
[0211]
本发明根据数控机床组件的结构特征，选取了组件运行状态监测数据类型，并确定组件运行状态监测部位，搭建数控机床组件运行状态监测平台，获取了组件的状态监测数据；对获得的监测信号进行了趋势项消除、降噪以及特征提取等信号处理；考虑不同的特征对机床组件退化或故障的响应不同，应用pca对信号特征进行降维，以降维后的主成分作为内部协变量，将组件的负载和速度作为外部协变量，建立了威布尔回归模型，通过参数估计获得了模型最优参数值，进而实现机床关键组件剩余寿命预测，准确合理的剩余寿命预测结果，对于降低生产制造企业停机损失，提高企业生产效率，增加经济收益意义显著。