CN111143990A - 一种敏感测点选择及融合的机床铣刀剩余寿命预测方法 - Google Patents

一种敏感测点选择及融合的机床铣刀剩余寿命预测方法 Download PDF

Info

Publication number
CN111143990A
CN111143990A CN201911353616.6A CN201911353616A CN111143990A CN 111143990 A CN111143990 A CN 111143990A CN 201911353616 A CN201911353616 A CN 201911353616A CN 111143990 A CN111143990 A CN 111143990A
Authority
CN
China
Prior art keywords
measuring point
state
measuring
milling cutter
point group
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.)
Granted
Application number
CN201911353616.6A
Other languages
English (en)
Other versions
CN111143990B (zh
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 Jiaotong University
Original Assignee
Xian Jiaotong 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 Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201911353616.6A priority Critical patent/CN111143990B/zh
Publication of CN111143990A publication Critical patent/CN111143990A/zh
Application granted granted Critical
Publication of CN111143990B publication Critical patent/CN111143990B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

一种敏感测点选择及融合的机床铣刀剩余寿命预测方法，首先建立机床铣刀退化的状态空间模型，其中状态传递函数表示隐藏系统状态值随时间的变化关系，观测函数表示不同传感器测点信号与系统状态之间的函数关系；其次利用单位极大似然估计法估计状态传递函数参数，最小二乘拟合法估计观测函数参数；然后利用斯皮尔曼相关系数绝对值选择最优测点组，粒子滤波算法对多测点数据进行融合；最后根据退化过程模拟算法实时预测铣刀剩余使用寿命的分布情况；本发明更加直观地解释了铣刀系统状态退化过程和多测点信号之间的关系，并对多测点数据进行融合，从而提高机床铣刀状态估计和剩余寿命预测的精度。

Description

一种敏感测点选择及融合的机床铣刀剩余寿命预测方法
技术领域
本发明属于机床刀具剩余寿命预测技术领域，具体涉及一种敏感测点选择及融合的机床铣刀剩余寿命预测方法。
背景技术
现代工业发展迅速，精密机床也朝着自动化、高精度、高效化的方向发展。然而在加工过程中，由于工件挤压、摩擦、冲击等因素影响，机床刀具会产生性能退化，出现崩刃、磨损等现象，轻则增大工件加工误差，重则导致设备损坏和人员伤亡。因此，对机床刀具进行剩余寿命预测以保障其运行可靠性变得越发重要。
随着工业物联网的飞速发展，越来越多的传感器被用于工业系统的健康状态监测中，构成了监测传感器网络，所采集的海量数据为剩余寿命预测带来了巨大的信息资源和技术挑战。如何选择传感器敏感测点以及如何融合多传感器测点数据以提高预测性能和精度成为急需解决的关键问题。现有的模型大多侧重于使用单一测点信息或是多个测点信号简单组合成的复合健康指标进行剩余寿命预测，缺乏通用性和扩展性，对于机床铣刀，不能够准确反映其退化过程。因此，采用敏感测点选择及融合的方式来提高机床铣刀剩余寿命预测的精度具有十分重要的意义。
发明内容
本发明的目的在于克服上述现有模型的缺点，提供一种敏感测点选择及融合的机床铣刀剩余寿命预测方法，通过一个状态传递函数表示隐藏系统状态值随时间的变化，一个观测函数表示不同测点信号之间的相关性以及它们的不同变化趋势，采用粒子滤波算法对预测性能优良的多个测点数据进行融合；从状态空间建模的角度反映多测点信号性能，并通过测点优选算法选择最优测点组，从而更加直观地解释了铣刀系统状态退化过程和多测点信号之间的关系，提高机床铣刀状态估计和剩余寿命预测的精度。
为达到上述目的，本发明采取的技术方案如下：
一种敏感测点选择及融合的机床铣刀剩余寿命预测方法，包括以下步骤：
1)建立铣刀系统状态退化模型：
建立如下多元状态空间模型：
xk＝xk-1+λωΔtk-1+εk-1 (1)
Yk＝AΦk+B+Vk (2)
式(1)为状态传递函数，用来描述系统状态的退化过程，其中，xk是tk时刻的系统状态值，表征铣刀的真实健康状态退化程度，即后刀面磨损量；λ是描述单位退化速率变化的参数；ω是与时间tk-1、状态xk-1和参数矩阵θ有关的退化速率函数；时间间隔Δtk-1＝tk-1-tk-2，随机噪声εk-1服从正态分布
Figure BDA0002335309520000021
式(2)为观测函数，用来描述后刀面磨损量与多测点信号之间的关系，假设共有P个传感器测点，则观测矩阵
Figure BDA0002335309520000022
Figure BDA0002335309520000023
表示第m个测点在tk时刻的观测值；比例系数矩阵A＝diag(α1，α2，…，αm，…，αP)，αm表示第m个测点的比例系数；位置系数矩阵B＝(β1，β2，…，βm，…，βP)T，βm表示第m个测点的位置系数；测点退化趋势矩阵Φk＝(φ(xk，γ1)，φ(xk，γ2)，…，φ(xk，γm)，…，φ(xk，γP))T，φ(xk，γm)表示第m个测点在状态xk下的退化趋势，γm是第m个测点监测信号的参数矩阵；噪声矩阵
Figure BDA0002335309520000024
Figure BDA0002335309520000025
服从P维多元正态分布N(0，∑)，
Figure BDA0002335309520000026
是第m个测点监测信号的随机噪声；
2)估计多元状态空间模型参数：
上述多元状态空间模型中，λ、θ、
Figure BDA0002335309520000027
A、B、γm、∑为待估计参数；
假设已获得N个训练样本及其退化信号，对于第n个样本，在离散时间序列
Figure BDA0002335309520000031
下，对应的后刀面磨损量观测序列记为
Figure BDA0002335309520000032
其中Kn表示第n个样本的采样总数，第m个测点的测量序列记为
Figure BDA0002335309520000033
使用训练样本数据进行多元状态空间模型的参数估计，过程如下：
2.1)利用单位极大似然估计法估计状态传递函数的参数(λ1，λ2，…，λN)、θ、
Figure BDA0002335309520000034
令
Figure BDA0002335309520000035
显然ΔXn服从多元正态分布
Figure BDA0002335309520000036
其中
Figure BDA0002335309520000037
ωn，k(θ)＝ω(xn，k，tn，k，θ)Δtn，k，
Figure BDA00023353095200000318
未知参数
Figure BDA0002335309520000038
的对数似然函数如下：
Figure BDA0002335309520000039
求式(3)关于λn和
Figure BDA00023353095200000310
的一阶偏导，令其等于零，得λn和
Figure BDA00023353095200000311
关于θ的函数，
Figure BDA00023353095200000312
Figure BDA00023353095200000313
将式(4)和式(5)代入式(3)，由式(6)求得θ的估计值
Figure BDA00023353095200000314
Figure BDA00023353095200000315
再将
Figure BDA00023353095200000316
代入式(4)和式(5)，得(λ1，λ2，…，λN)和
Figure BDA00023353095200000317
2.2)分步估计观测函数的参数A、B、γm、∑；
根据式(2)知，参数αm、βm和γm与每个测点自身有关，协方差矩阵∑与不同测点之间的相关性有关，为反映不同测点之间的变化和相关性，分两步进行参数估计；
2.2.1)单独分析每个测点信号，估计参数αm、βm和γm，
利用最小二乘拟合算法，根据式(7)分别估计参数αm、βm和γm：
Figure BDA0002335309520000041
2.2.2)联合分析所有测点信号，估计协方差矩阵∑，
将第n个样本在tk时刻的多测点测量向量记为
Figure BDA0002335309520000042
服从多元正态分布N(AΦn，k+B，∑)，其中Φn，k＝(φ(xn，k，γ1)，φ(xn，k，γ2)，…，φ(xn，k，γP))T，∑的对数似然函数如下：
Figure BDA0002335309520000043
求式(8)关于∑的一阶偏导，令其等于零，得∑的极大似然估计值，如式(9)所示：
Figure BDA0002335309520000044
3)最优测点组选择：
首先根据测点信号和后刀面磨损量之间的相关系数对P个测点进行优先级排序，然后根据优先级顺序生成P个测点组，最后从中选出最能准确估计后刀面磨损量的测点组作为最优测点组，具体方法如下：
3.1)定量估计后刀面磨损量与测点信号的相关性：
使用斯皮尔曼相关系数绝对值(ASC)进行非线性相关性度量，如式(10)所示，
Figure BDA0002335309520000045
其中
Figure BDA0002335309520000051
表示第m个测点的测量序列和刀具磨损厚度测量序列之间的斯皮尔曼相关系数绝对值，且
Figure BDA0002335309520000052
Figure BDA0002335309520000053
和
Figure BDA0002335309520000054
分别表示xn，k和
Figure BDA0002335309520000055
在各自序列中的排序；
3.2)生成测点组：
按N个训练样本的平均ASC值对P个测点进行降序排列，取前m个测点作为第m个测点组，共获得P个测点组；
3.3)选择最优测点组：
为不失一般性，假设测点组M是P个测点的任意组合，使用基于模糊重采样的粒子滤波算法，对每一个测点组进行训练样本退化速率和磨损厚度的估计；在初始时刻，生成一系列退化速率粒子和系统状态粒子，记为
Figure BDA0002335309520000056
其中Ns是粒子数，使用线性插值方法采样，从序列
Figure BDA0002335309520000057
和
Figure BDA0002335309520000058
中分别获得退化速率粒子
Figure BDA0002335309520000059
和系统状态粒子
Figure BDA00023353095200000510
每次系统状态粒子按照式(11)进行单步传送，
Figure BDA00023353095200000511
粒子权重初始化为
Figure BDA00023353095200000512
并按照式(12)进行更新，
Figure BDA00023353095200000513
其中，|M|表示测点组M的测点数量，YM，n，k是测点组M中，第n个样本在tk时刻的观测值，AM、BM、
Figure BDA00023353095200000514
和∑M分别是测点组M所对应的A、B、
Figure BDA00023353095200000515
和∑的子矩阵；
使用模糊重采样算法对粒子进行重采样，将粒子权重重置为
Figure BDA00023353095200000516
重采样的粒子记为
Figure BDA00023353095200000517
取退化速率粒子和系统状态粒子的中位数
Figure BDA00023353095200000518
和
Figure BDA00023353095200000519
作为各自的估计值，通过计算估计值和实际值之间的平均相对误差绝对值MARE来评估不同测点组的性能；
Figure BDA0002335309520000061
MARE值越低，状态估计值越准确，选择MARE值最低的测点组作为最优测点组；
4)剩余使用寿命预测：
首先，根据步骤3.3)中得到的最优测点组的监测信号，使用步骤3.2)中的式(11)和式(12)估计测试样本的系统状态；然后，将估计的模型参数和系统状态代入式(1)即得系统的状态退化过程；最后，根据退化过程模拟算法实时预测铣刀剩余使用寿命的分布情况。
本发明的有益效果为：
本发明针对机床铣刀，提出一种敏感测点选择及融合的机床铣刀剩余寿命预测方法，使用状态传递函数表示铣刀隐藏系统状态的退化过程，观测函数表示不同传感器测定信号之间的相关性和变化趋势，同时使用粒子滤波算法对预测性能优良的最优测点组数据进行融合。与使用单一测点信息或复合健康指标的传统方法相比，本发明更加直观地解释了铣刀系统状态退化过程和多测点信号之间的关系，并对多测点数据进行融合，从而提高机床铣刀状态估计和剩余寿命预测的精度。
附图说明
图1为本发明方法的流程图。
图2为实施例9个信号通道的RMS值。
图3为实施例测点选择结果过程，图(a)表示9个信号通道的ASC值，图(b)表示9个测点组的MARE值。
图4为实施例在三种预测方法下的预测结果。
具体实施方式
下面结合附图和实施例对本发明做进一步阐述。
参照图1，一种敏感测点选择及融合的机床铣刀剩余寿命预测方法，包括以下步骤：
1)建立铣刀系统状态退化模型：
建立如下多元状态空间模型：
xk＝xk-1+λωΔtk-1+εk-1 (1)
Yk＝AΦk+B+Vk (2)
式(1)为状态传递函数，用来描述系统状态的退化过程，其中，xk是tk时刻的系统状态值，表征铣刀的真实健康状态退化程度，即后刀面磨损量；λ是描述单位退化速率变化的参数；ω是与时间tk-1、状态xk-1和参数矩阵θ有关的退化速率函数；时间间隔Δtk-1＝tk-1-tk-2，随机噪声εk-1服从正态分布
Figure BDA0002335309520000071
式(2)为观测函数，用来描述后刀面磨损量与多测点信号之间的关系，假设共有P个传感器测点，则观测矩阵
Figure BDA0002335309520000072
Figure BDA0002335309520000073
表示第m个测点在tk时刻的观测值；比例系数矩阵A＝diag(α1，α2，…，αm，…，αP)，αm表示第m个测点的比例系数；位置系数矩阵B＝(β1，β2，…，βm，…，βP)T，βm表示第m个测点的位置系数；测点退化趋势矩阵Φk＝(φ(xk，γ1)，φ(xk，γ2)，…，φ(xk，γm)，…，φ(xk，γP))T，φ(xk，γm)表示第m个测点在状态xk下的退化趋势，γm是第m个测点监测信号的参数矩阵；噪声矩阵
Figure BDA0002335309520000074
Figure BDA0002335309520000075
服从P维多元正态分布N(O，∑)，
Figure BDA0002335309520000076
是第m个测点监测信号的随机噪声；
2)估计多元状态空间模型参数：
上述多元状态空间模型中，λ、θ、
Figure BDA0002335309520000077
A、B、γm、∑为待估计参数；
假设已获得N个训练样本及其退化信号，对于第n个样本，在离散时间序列
Figure BDA0002335309520000078
下，对应的后刀面磨损量观测序列记为
Figure BDA0002335309520000079
其中Kn表示第n个样本的采样总数，第m个测点的测量序列记为
Figure BDA0002335309520000081
使用训练样本数据进行多元状态空间模型的参数估计，过程如下：
2.1)利用单位极大似然估计法估计状态传递函数的参数(λ1，λ2，…，λN)、θ、
Figure BDA0002335309520000082
令
Figure BDA0002335309520000083
显然ΔXn服从多元正态分布
Figure BDA0002335309520000084
其中
Figure BDA0002335309520000085
ωn，k(θ)＝ω(xn，k，tn，k，θ)Δtn，k，
Figure BDA0002335309520000086
未知参数
Figure BDA0002335309520000087
的对数似然函数如下：
Figure BDA0002335309520000088
求式(3)关于λn和
Figure BDA0002335309520000089
的一阶偏导，令其等于零，得λn和
Figure BDA00023353095200000810
关于θ的函数，
Figure BDA00023353095200000811
Figure BDA00023353095200000812
将式(4)和式(5)代入式(3)，由式(6)求得θ的估计值
Figure BDA00023353095200000813
Figure BDA00023353095200000814
再将
Figure BDA00023353095200000815
代入式(4)和式(5)，得(λ1，λ2，…，λN)和
Figure BDA00023353095200000816
2.2)分步估计观测函数的参数A、B、γm、∑；
根据式(2)知，参数αm、βm和γm与每个测点自身有关，协方差矩阵∑与不同测点之间的相关性有关，为反映不同测点之间的变化和相关性，分两步进行参数估计；
2.2.1)单独分析每个测点信号，估计参数αm、βm和γm，
利用最小二乘拟合算法，根据式(7)分别估计参数αm、βm和γm：
Figure BDA0002335309520000091
2.2.2)联合分析所有测点信号，估计协方差矩阵∑，
将第n个样本在tk时刻的多测点测量向量记为
Figure BDA0002335309520000092
服从多元正态分布N(AΦn，k+B，∑)，其中Φn，k＝(φ(xn，k，γ1)，φ(xn，k，γ2)，…，φ(xn，k，γP))T，∑的对数似然函数如下：
Figure BDA0002335309520000093
求式(8)关于∑的一阶偏导，令其等于零，得∑的极大似然估计值，如式(9)所示：
Figure BDA0002335309520000094
3)最优测点组选择：
首先根据测点信号和后刀面磨损量之间的相关系数对P个测点进行优先级排序，然后根据优先级顺序生成P个测点组，最后从中选出最能准确估计后刀面磨损量的测点组作为最优测点组，具体方法如下：
3.1)定量估计后刀面磨损量与测点信号的相关性：
使用斯皮尔曼相关系数绝对值(ASC)进行非线性相关性度量，如式(10)所示，
Figure BDA0002335309520000095
其中
Figure BDA0002335309520000096
表示第m个测点的测量序列和刀具磨损厚度测量序列之间的斯皮尔曼相关系数绝对值，且
Figure BDA0002335309520000097
Figure BDA0002335309520000098
和
Figure BDA0002335309520000099
分别表示xn，k和
Figure BDA00023353095200000910
在各自序列中的排序；
3.2)生成测点组：
按N个训练样本的平均ASC值对P个测点进行降序排列，取前m个测点作为第m个测点组，共获得P个测点组；
3.3)选择最优测点组：
为不失一般性，假设测点组M是P个测点的任意组合，使用基于模糊重采样的粒子滤波算法，对每一个测点组进行训练样本退化速率和磨损厚度的估计；在初始时刻，生成一系列退化速率粒子和系统状态粒子，记为
Figure BDA0002335309520000101
其中NS是粒子数，使用线性插值方法采样，从序列
Figure BDA0002335309520000102
和
Figure BDA0002335309520000103
中分别获得退化速率粒子
Figure BDA0002335309520000104
和系统状态粒子
Figure BDA0002335309520000105
每次系统状态粒子按照式(11)进行单步传送，
Figure BDA0002335309520000106
粒子权重初始化为
Figure BDA0002335309520000107
并按照式(12)进行更新，
Figure BDA0002335309520000108
其中，|M|表示测点组M的测点数量，YM，n，k是测点组M中，第n个样本在tk时刻的观测值，AM、BM、
Figure BDA0002335309520000109
和∑M分别是测点组M所对应的A、B、
Figure BDA00023353095200001010
和∑的子矩阵；
使用模糊重采样算法对粒子进行重采样，将粒子权重重置为
Figure BDA00023353095200001011
重采样的粒子记为
Figure BDA00023353095200001012
取退化速率粒子和系统状态粒子的中位数
Figure BDA00023353095200001013
和
Figure BDA00023353095200001014
作为各自的估计值，通过计算估计值和实际值之间的平均相对误差绝对值MARE来评估不同测点组的性能；
Figure BDA00023353095200001015
MARE值越低，状态估计值越准确，选择MARE值最低的测点组作为最优测点组；
4)剩余使用寿命预测：
首先，根据步骤3.3)中得到的最优测点组的监测信号，使用步骤3.2)中的式(11)和式(12)估计测试样本的系统状态；然后，将估计的模型参数和系统状态代入式(1)即得系统的状态退化过程；最后，根据退化过程模拟算法实时预测铣刀剩余使用寿命的分布情况。
使用在某数控铣床铣削过程中采集到的刀具磨损监测数据对本发明方法进行验证，铣削具体参数如下：转速4500rpm，切削速度400mm/min，切削深度1.5mm，工件材料45钢，工件尺寸100×60×60mm，共记录7把刀具C1～C7在其整个寿命周期内的后刀面磨损量数据。采用的传感器共有5种不同类型，即三向加速度传感器、三向力传感器、电流传感器、声传感器和声发射传感器，其中三向加速度传感器安装在主轴上，三向力传感器安装在工件下方，两者都有X，Y和Z三个不同方向的通道，电流传感器用于监视主轴电机的电流，声传感器和声发射传感器安装在工件旁边，分别监测声音信号和声发射信号。5个传感器共产生9个信号通道，记为P1～P9，计算每个信号的均方根值(RMS)，用于预测刀具的剩余使用寿命，9个信号通道的RMS值如图2所示。
为了描述刀具的时变退化速率，使用θexp(θτ)来表示式(6)中的退化速率ω(x(τ)，τ，θ)，并用多项式函数
Figure BDA0002335309520000111
表示状态值和观测值之间的关系。图3中(a)和(b)分别给出了9个信号通道的ASC值和9个测点组的MARE值。从图3的(a)可以看出，力信号P4～P6、电流信号P7和声信号P8的ASC值比其他信号通道更高。在图2的RMS曲线中，这5个信号通道的增长趋势也比其他通道更明显。两者结果一致，因此证明ASC值可以有效评估测点信号的灵敏度。在图3的(b)中，第5个由{5、7、4、8、6}组成的测点组的MARE值最低，即该组合能够获得最准确的状态估计结果。因此，选择该测点组作为剩余寿命预测的最优测点组。
使用留一法交叉验证对实验数据进行分析，将六把刀具用作训练样本，一把刀具作为测试样本，选取上述最优测点组，使用本发明方法预测测试刀具的剩余使用寿命。为了进行比较，同时使用其他两种方法进行预测，状态空间模型相同，但选用不同的测点：一种只选择P5通道，记为“方法1”；另一种选择P1～P9全部通道，记为“方法2”。三种方法预测结果如图4所示，从图中可知，本发明方法明显优于另外两种方法。
本发明方法适用于各类机床铣刀的剩余寿命预测，在实际应用中，实施者可以利用本方法选择最优测点组，并进行多测点数据融合，更加直观地解释了系统状态退化过程和多测点信号之间的关系，有助于提高铣刀状态估计和剩余寿命预测的精度。应当指出，在不脱离本发明构想的前提下，对本方法所做的调整和变形，也应视为本发明的保护范围。

Claims (1)

1.一种敏感测点选择及融合的机床铣刀剩余寿命预测方法，其特征在于，包括以下步骤：
1)建立铣刀系统状态退化模型：
建立如下多元状态空间模型：
xk＝xk-1+λωΔtk-1+εk-1 (1)
Yk＝AΦk+B+Vk (2)
式(1)为状态传递函数，用来描述系统状态的退化过程，其中，xk是tk时刻的系统状态值，表征铣刀的真实健康状态退化程度，即后刀面磨损量；λ是描述单位退化速率变化的参数；ω是与时间tk-1、状态xk-1和参数矩阵θ有关的退化速率函数；时间间隔Δtk-1＝tk-1-tk-2，随机噪声εk-1服从正态分布
Figure FDA0002335309510000011
式(2)为观测函数，用来描述后刀面磨损量与多测点信号之间的关系，假设共有P个传感器测点，则观测矩阵
Figure FDA0002335309510000012
Figure FDA0002335309510000013
表示第m个测点在tk时刻的观测值；比例系数矩阵A＝diag(α1，α2，…，αm，…，αP)，αm表示第m个测点的比例系数；位置系数矩阵B＝(β1，β2，…，βm，…，βP)T，βm表示第m个测点的位置系数；测点退化趋势矩阵Φk＝(φ(xk，γ1)，φ(xk，γ2)，…，φ(xk，γm)，…，φ(xk，γP))T，φ(xk，γm)表示第m个测点在状态xk下的退化趋势，γm是第m个测点监测信号的参数矩阵；噪声矩阵Vk＝
Figure FDA0002335309510000014
服从P维多元正态分布N(0，∑)，
Figure FDA0002335309510000015
是第m个测点监测信号的随机噪声；
2)估计多元状态空间模型参数：
上述多元状态空间模型中，λ、θ、
Figure FDA0002335309510000016
A、B、γm、∑为待估计参数；
假设已获得N个训练样本及其退化信号，对于第n个样本，在离散时间序列
Figure FDA0002335309510000021
下，对应的后刀面磨损量观测序列记为
Figure FDA0002335309510000022
其中Kn表示第n个样本的采样总数，第m个测点的测量序列记为
Figure FDA0002335309510000023
使用训练样本数据进行多元状态空间模型的参数估计，过程如下：
2.1)利用单位极大似然估计法估计状态传递函数的参数(λ1，λ2，…，λN)、θ、
Figure FDA0002335309510000024
令
Figure FDA0002335309510000025
显然ΔXn服从多元正态分布
Figure FDA0002335309510000026
其中
Figure FDA0002335309510000027
ωn，k(θ)＝ω(xn，k，tn，k，θ)Δtn，k，
Figure FDA0002335309510000028
未知参数
Figure FDA0002335309510000029
的对数似然函数如下：
Figure FDA00023353095100000210
求式(3)关于λn和
Figure FDA00023353095100000211
的一阶偏导，令其等于零，得λn和
Figure FDA00023353095100000212
关于θ的函数，
Figure FDA00023353095100000213
Figure FDA00023353095100000214
将式(4)和式(5)代入式(3)，由式(6)求得θ的估计值
Figure FDA00023353095100000215
Figure FDA00023353095100000216
再将
Figure FDA00023353095100000217
代入式(4)和式(5)，得(λ1，λ2，…，λN)和
Figure FDA00023353095100000218
2.2)分步估计观测函数的参数A、B、γm、∑；
根据式(2)知，参数αm、βm和γm与每个测点自身有关，协方差矩阵∑与不同测点之间的相关性有关，为反映不同测点之间的变化和相关性，分两步进行参数估计；
2.2.1)单独分析每个测点信号，估计参数αm、βm和γm，
利用最小二乘拟合算法，根据式(7)分别估计参数αm、βm和γm：
Figure FDA0002335309510000031
2.2.2)联合分析所有测点信号，估计协方差矩阵∑，
将第n个样本在tk时刻的多测点测量向量记为
Figure FDA0002335309510000032
服从多元正态分布N(AΦn，k+B，∑)，其中Φn，k＝(φ(xn，k，γ1)，φ(xn，k，γ2)，…，φ(xn，k，γP))T，∑的对数似然函数如下：
Figure FDA0002335309510000033
求式(8)关于∑的一阶偏导，令其等于零，得∑的极大似然估计值，如式(9)所示：
Figure FDA0002335309510000034
3)最优测点组选择：
首先根据测点信号和后刀面磨损量之间的相关系数对P个测点进行优先级排序，然后根据优先级顺序生成P个测点组，最后从中选出最能准确估计后刀面磨损量的测点组作为最优测点组，具体方法如下：
3.1)定量估计后刀面磨损量与测点信号的相关性：
使用斯皮尔曼相关系数绝对值(ASC)进行非线性相关性度量，如式(10)所示，
Figure FDA0002335309510000035
其中
Figure FDA0002335309510000041
表示第m个测点的测量序列和刀具磨损厚度测量序列之间的斯皮尔曼相关系数绝对值，且
Figure FDA0002335309510000042
Figure FDA0002335309510000043
和
Figure FDA0002335309510000044
分别表示xn，k和
Figure FDA0002335309510000045
在各自序列中的排序；
3.2)生成测点组：
按N个训练样本的平均ASC值对P个测点进行降序排列，取前m个测点作为第m个测点组，共获得P个测点组；
3.3)选择最优测点组：
为不失一般性，假设测点组M是P个测点的任意组合，使用基于模糊重采样的粒子滤波算法，对每一个测点组进行训练样本退化速率和磨损厚度的估计；在初始时刻，生成一系列退化速率粒子和系统状态粒子，记为
Figure FDA0002335309510000046
其中Ns是粒子数，使用线性插值方法采样，从序列
Figure FDA0002335309510000047
和
Figure FDA0002335309510000048
中分别获得退化速率粒子
Figure FDA0002335309510000049
和系统状态粒子
Figure FDA00023353095100000410
每次系统状态粒子按照式(11)进行单步传送，
Figure FDA00023353095100000411
粒子权重初始化为
Figure FDA00023353095100000412
并按照式(12)进行更新，
Figure FDA00023353095100000413
其中，|M|表示测点组M的测点数量，YM，n，k是测点组M中，第n个样本在tk时刻的观测值，AM、BM、
Figure FDA00023353095100000414
和∑M分别是测点组M所对应的A、B、
Figure FDA00023353095100000415
和∑的子矩阵；
使用模糊重采样算法对粒子进行重采样，将粒子权重重置为
Figure FDA00023353095100000416
重采样的粒子记为
Figure FDA00023353095100000417
取退化速率粒子和系统状态粒子的中位数
Figure FDA00023353095100000418
和
Figure FDA00023353095100000419
作为各自的估计值，通过计算估计值和实际值之间的平均相对误差绝对值MARE来评估不同测点组的性能；
Figure FDA0002335309510000051
MARE值越低，状态估计值越准确，选择MARE值最低的测点组作为最优测点组；
4)剩余使用寿命预测：
首先，根据步骤3.3)中得到的最优测点组的监测信号，使用步骤3.2)中的式(11)和式(12)估计测试样本的系统状态；然后，将估计的模型参数和系统状态代入式(1)即得系统的状态退化过程；最后，根据退化过程模拟算法实时预测铣刀剩余使用寿命的分布情况。
CN201911353616.6A 2019-12-25 2019-12-25 一种敏感测点选择及融合的机床铣刀剩余寿命预测方法 Active CN111143990B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911353616.6A CN111143990B (zh) 2019-12-25 2019-12-25 一种敏感测点选择及融合的机床铣刀剩余寿命预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911353616.6A CN111143990B (zh) 2019-12-25 2019-12-25 一种敏感测点选择及融合的机床铣刀剩余寿命预测方法

Publications (2)

Publication Number Publication Date
CN111143990A true CN111143990A (zh) 2020-05-12
CN111143990B CN111143990B (zh) 2021-11-16

Family

ID=70519801

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911353616.6A Active CN111143990B (zh) 2019-12-25 2019-12-25 一种敏感测点选择及融合的机床铣刀剩余寿命预测方法

Country Status (1)

Country Link
CN (1) CN111143990B (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112855467A (zh) * 2021-03-22 2021-05-28 西安交通大学 一种风力发电机基准工况转换方法
CN112881518A (zh) * 2021-01-08 2021-06-01 东冶及策河北能源技术有限公司 一种动态滤波补偿器剩余寿命预测方法
CN113467375A (zh) * 2021-01-06 2021-10-01 南京航空航天大学 基于虚拟测量的智能数控加工刀具剩余可用寿命预测方法
CN113688534A (zh) * 2021-09-02 2021-11-23 江苏师范大学 一种基于多特征融合模型寻找最优铣削参数的研究方法
CN113847949A (zh) * 2021-09-23 2021-12-28 徐州万达回转支承有限公司 一种基于传感器信息融合的多工况砂带磨损状态在线检测方法
CN117829002A (zh) * 2024-03-05 2024-04-05 深圳市明谋科技有限公司 电力线缆的老化诊断监测方法及系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104792529A (zh) * 2015-04-12 2015-07-22 北京化工大学 基于状态空间模型的滚动轴承寿命预测方法
US20150349385A1 (en) * 2014-04-01 2015-12-03 Medtronic, Inc. Method and System for Predicting Useful Life of a Rechargeable Battery
CN109212966A (zh) * 2018-08-14 2019-01-15 西安交通大学 一种多工况动态基准化的机械设备剩余寿命预测方法
CN109376401A (zh) * 2018-09-29 2019-02-22 西安交通大学 一种自适应多源信息优选与融合的机械剩余寿命预测方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150349385A1 (en) * 2014-04-01 2015-12-03 Medtronic, Inc. Method and System for Predicting Useful Life of a Rechargeable Battery
CN104792529A (zh) * 2015-04-12 2015-07-22 北京化工大学 基于状态空间模型的滚动轴承寿命预测方法
CN109212966A (zh) * 2018-08-14 2019-01-15 西安交通大学 一种多工况动态基准化的机械设备剩余寿命预测方法
CN109376401A (zh) * 2018-09-29 2019-02-22 西安交通大学 一种自适应多源信息优选与融合的机械剩余寿命预测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
LI NAIPENG等: "Remaining useful life prediction of machinery under time-varying operating conditions based on a two-factor state-space model", 《RELIABILITY ENGINEERING AND SYSTEM SAFETY》 *
NAIPENG LI等: "Machine remaining useful life prediction considering unit-to-unit variability", 《2017 IEEE INTERNATIONAL CONFERENCE ON PROGNOSTICS AND HEALTH MANAGEMENT》 *
YAWEI HU等: "Remaining useful life model and assessment of mechanical products: a brief review and a note on the state space model method", 《CHINESE JOURNAL OF MECHANICAL ENGINEERING》 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113467375A (zh) * 2021-01-06 2021-10-01 南京航空航天大学 基于虚拟测量的智能数控加工刀具剩余可用寿命预测方法
CN112881518A (zh) * 2021-01-08 2021-06-01 东冶及策河北能源技术有限公司 一种动态滤波补偿器剩余寿命预测方法
CN112881518B (zh) * 2021-01-08 2022-07-12 东冶及策河北能源技术有限公司 一种动态滤波补偿器剩余寿命预测方法
CN112855467A (zh) * 2021-03-22 2021-05-28 西安交通大学 一种风力发电机基准工况转换方法
CN112855467B (zh) * 2021-03-22 2022-01-11 西安交通大学 一种风力发电机基准工况转换方法
CN113688534A (zh) * 2021-09-02 2021-11-23 江苏师范大学 一种基于多特征融合模型寻找最优铣削参数的研究方法
CN113688534B (zh) * 2021-09-02 2024-04-05 苏州莱库航空装备科技有限公司 一种基于多特征融合模型寻找最优铣削参数的研究方法
CN113847949A (zh) * 2021-09-23 2021-12-28 徐州万达回转支承有限公司 一种基于传感器信息融合的多工况砂带磨损状态在线检测方法
CN117829002A (zh) * 2024-03-05 2024-04-05 深圳市明谋科技有限公司 电力线缆的老化诊断监测方法及系统
CN117829002B (zh) * 2024-03-05 2024-05-14 深圳市明谋科技有限公司 电力线缆的老化诊断监测方法及系统

Also Published As

Publication number Publication date
CN111143990B (zh) 2021-11-16

Similar Documents

Publication Publication Date Title
CN111143990B (zh) 一种敏感测点选择及融合的机床铣刀剩余寿命预测方法
CN110647943B (zh) 基于演化数据聚类分析的切削刀具磨损监测方法
TWI640390B (zh) 刀具磨耗監測與預測方法
TWI481978B (zh) 工具機之加工品質的預測方法
Salonitis et al. Reliability assessment of cutting tool life based on surrogate approximation methods
Zhang et al. Modelling and prediction of tool wear using LS-SVM in milling operation
Zhao et al. An integrated prognostics method for failure time prediction of gears subject to the surface wear failure mode
CN109376401B (zh) 一种自适应多源信息优选与融合的机械剩余寿命预测方法
Jemielniak et al. Diagnosis of tool wear based on cutting forces and acoustic emission measures as inputs to a neural network
CN109333160B (zh) 高温合金钻削过程钻头磨损形式及磨损状态的在线监测方法
Brecher et al. Use of NC kernel data for surface roughness monitoring in milling operations
TWI584134B (zh) 製程異因分析方法與製程異因分析系統
CN113458873B (zh) 一种刀具磨损量和剩余寿命预测的方法
CN110209150B (zh) 基于多工序故障影响的作业车间调度方案鲁棒性测度方法
KR20190043038A (ko) 공작기계의 진동 적응제어를 위한 진동신호 처리 방법
Suhail et al. Workpiece surface temperature for in-process surface roughness prediction using response surface methodology
Zhang et al. Modeling of tool wear for ball end milling cutter based on shape mapping
Tansel et al. Genetic tool monitor (GTM) for micro-end-milling operations
CN113849901B (zh) 针对接触换热系数辨识的改进自适应优化方法及系统
El Ouafi et al. An ANN based multi-sensor integration approach for in-process monitoring of product quality in turning operations
Wiciak-Pikuła et al. Surface roughness and forces prediction of milling Inconel 718 with neural network
JP6809999B2 (ja) 導電性を備えた樹脂成形品の劣化箇所推定方法
Kayabaşi et al. On-line surface roughness prediction by using a probabilistic approach for end-milling
CN111708323A (zh) 一种具有热变形误差补偿的五轴小龙门数控加工中心
CN117961645B (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
GR01 Patent grant
GR01 Patent grant