CN114528879A - 轴承故障周期自动检测方法和系统 - Google Patents

轴承故障周期自动检测方法和系统 Download PDF

Info

Publication number
CN114528879A
CN114528879A CN202210143261.3A CN202210143261A CN114528879A CN 114528879 A CN114528879 A CN 114528879A CN 202210143261 A CN202210143261 A CN 202210143261A CN 114528879 A CN114528879 A CN 114528879A
Authority
CN
China
Prior art keywords
signal
peak
statistic
fault
signals
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
CN202210143261.3A
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.)
Shanghai Jiaotong University
Original Assignee
Shanghai 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 Shanghai Jiaotong University filed Critical Shanghai Jiaotong University
Priority to CN202210143261.3A priority Critical patent/CN114528879A/zh
Publication of CN114528879A publication Critical patent/CN114528879A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/12Classification; Matching
    • G06F2218/16Classification; Matching by matching signal segments
    • G06F2218/20Classification; Matching by matching signal segments by applying autoregressive analysis
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01MTESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M13/00Testing of machine parts
    • G01M13/04Bearings
    • G01M13/045Acoustic or vibration analysis
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08Feature extraction
    • G06F2218/10Feature extraction by analysing the shape of a waveform, e.g. extracting parameters relating to peaks

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Signal Processing (AREA)
  • General Engineering & Computer Science (AREA)
  • Acoustics & Sound (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明提供了一种轴承故障周期自动检测方法和系统,包括:通过传感器从故障轴承上收集振动信号;将振动信号进行划分重组,拼合成两段子信号;对每段子信号采用NRC算法处理,加权求和得到EnNRC统计量;构造缩减置信区间;从EnNRC统计量的最高峰开始,向之前的波峰构造波峰差统计量,并构建动态置信区间,检查区间内是否有其他波峰,若存在则从区间内的最高波峰开始下一次假设检验,迭代到区间内不存在其他波峰为止,此时的波峰为算法检测到的故障真实周期。本发明通过迭代假设检验比较存在干扰效应的周期性波峰及其他非周期波峰,逐步消除局部最优的误导效应,提高了周期检测的精度。

Description

轴承故障周期自动检测方法和系统
技术领域
本发明涉及轴承故障检测技术领域,具体地,涉及一种轴承故障周期自动检测方法和系统。
背景技术
轴承故障检测起源于工业发展初期的恒定工况,且已经积累了一定成果,该领域的研究已相对成熟。传统的检测技术大多基于振动信号的时域、频域等信息,通过提取并分析振动信号的相关特征,与健康状态下产生的振动信号的特征指标对比,分析设备当前的运行状态,判断是否存在潜在故障。
为了克服噪声的掩蔽效应所带来的困难,学者们在时域上设计了一些基准方法。Rabiner等人提出了一种自相关函数(ACF)来提取周期信号的周期波形特征,可以从时域变换信号中直观地识别周期。Huang等人提出了改进版的ACF——加权平均幅度差函数(AMDF)——用于声波基音提取。Variability Method和NRC算法通过多段子信号的相关性计算可以克服强背景噪声对周期信号周期检测的掩蔽效应。如果提供足够长的信号,这些方法下信号的周期性特征将不会被背景噪声所掩盖。基于这些时域上的变换信号,目前已经研究出一些自动估计周期的针对性方法。例如,Martin等人提出了一种使用AMDF周期图谱来检测周期的方法,Antoni等人提出了一种ACF和周期图谱方法结合的混合方法来自动识别周期,称为自动周期法(Autoperiod)。Wang等人还提出了一种正交子空间分解方法,利用经Variability Method变化后的特征信号来识别周期。在中强度背景噪声下,也有一些频域方法用以识别周期,比如周期图(periodogram)。基于Ramanujan和的变换最近也被设计出并用于研究频域信号。为了自动识别噪声信号的周期,一些利用最大似然估计(ML)的统计方法也逐渐兴起。Wise等人提出了一种带有惩罚系数的ML方法,即最大似然基音估计(MLPE),用于从语音信号中估计基音周期。Ramírez等人提出了一种带正则项的近似最大似然估计方法,用于循环平稳信号的周期估计。这些方法通过定位似然函数的最大峰值来估计周期。
专利文献CN108106838B(申请号:CN201611122045.1)公开了一种机械故障精密诊断的整周期时域回归及报警方法,利用广义共振/共振解调方法和转速跟踪检测方法检测轴承、齿轮等故障的振动、冲击信号,通过自动修正每转跟踪采样点数后,先对原始信号按照分析的故障类型进行整周期截取,消除现场数据长度有限、非整周期采样造成频谱泄露及其与傅里叶变换频率分辨率之矛盾,可大大提高频谱定性的精准度,再对信号进行时域回归处理,其他无关信号成分的能量均不会产生干扰,最后利用级差公式对故障进行定量计算,可提高故障分类定量诊断的精准度。
现有研究中也提出了一些统计假设检验方法用于估计信号周期。例如,Dandawate等人提出了一种基于循环协方差和频谱的卡方检验方法来检测信号循环平稳性的存在并估计出信号的周期。最近,Horstmann等人中提出了一种多重假设检验方法,利用Holm’s的次序拒绝检验法来估计近似循环平稳信号的周期。这些假设检验方法是基于一种渐近广义似然比检验。通过确定具有最小p值的二进制检验的索引位置来估计周期。时域特征提取算法方法主要是在时域内提供能够消除背景噪声掩蔽效应的特征信号,且周期通常是通过基于人工经验的主观判断来识别。还有一些方法通过找出提取出的周期特征中的第一个波峰所在位置来确定周期,但由于轴承故障信号受到强背景噪声的干扰,最大波峰可能出现在任何倍数周期处,故这种识别方式没有任何的稳定性和说服力。由于强背景噪声的掩蔽效应,原始信号的频谱容易被噪声所掩盖,因此在强背景噪声下定位最大幅值频谱的简易判断方法不再能有效识别周期。因此,在频域上设计的传统方法不适用于一些信号被强背景噪声污染的应用。对于通过定位似然函数的最大峰值来估计周期的方法,似然函数在真实周期及其倍数上具有相似的统计性质,呈现出的峰值也是相似的,导致对基于最大峰值的周期识别有一定的误导效应。
因此,简单地定位具有最大似然估计值的位置以识别信号周期可能会识别出位于倍数周期位置的错误峰值,最终造成周期误诊。通过确定具有最小p值的二进制检验的索引位置来估计周期的方法同样面临周期倍数相似统计性质的干扰,即周期信号在真实周期及其倍数处会产生相似的p值。因此,假设检验方法同样可能会错误地将倍数周期识别为信号的真实周期。
因此,时域变换后的信号(如自相关系数、似然值和p值)在真实周期及其倍数处通常具有相似的视觉特征和统计特征,呈现出具有误导效应的局部最优现象。这种局部最优现象的误导效应使得针对性的周期估计算法、最大似然周期估计算法和假设检验类方法容易出现很高的误报和漏报概率,且除了假设检验算法,其他算法的误报、漏报概率都是不可控的,使得周期估计非常困难。
发明内容
针对现有技术中的缺陷,本发明的目的是提供一种轴承故障周期自动检测方法和系统。
根据本发明提供的轴承故障周期自动检测方法,包括:
步骤1:通过传感器从故障轴承上收集振动信号;
步骤2:将振动信号进行划分,将奇、偶片段分别重组,拼合成两段子信号;
步骤3:对每段子信号采用NRC算法处理,加权求和得到EnNRC统计量;
步骤4:利用该统计量进一步构造波峰差统计量;
步骤5:构造基于动态显著性水平的缩减置信区间;
步骤6:从EnNRC统计量的最高峰开始,向之前的波峰构造波峰差统计量,并构建动态置信区间,检查区间内是否有其他波峰,若存在其他波峰,则从区间内的最高波峰开始下一次假设检验,一直迭代到区间内不存在其他波峰为止,此时的波峰为算法检测到的故障真实周期。
优选的,将收集的故障轴承运转产生的周期脉冲信号y(t)分解为真实故障信号x(t)和背景噪声∈(t),则有:y(t)=x(t)+∈(t),其中,∈(t)~N(0,σ2);
定义时间序列y=[y(1),y(2),...,y(l)]T为y(t)的离散观测值,l为样本总数;x=[x(1),x(2),...,x(l)]T,∈=[∈(1),∈(2),...,∈(l)]T
则信号y表示为:y=x+∈,∈~N(0,σ2Il),Il为大小为l×l的单位矩阵;
设故障信号x(t)的周期为T0,T0=p0/fs,其中,p0>1为一个周期内的样本数目,fs为采样频率;
通过NRC算法将信号y划分为
Figure BDA0003507442830000034
段,第k个信号片段yk=[y((k-1)n+1),…,y(kn)],其中k=1,2,…,m;n是每段信号中样本数目;
则EnNRC统计量表示为:
Figure BDA0003507442830000031
其中,i、j为序列号;若l%n≠0,则信号被划分为2m+1段;
奇数段的信号片段长度为:n1=l-mn,其中
Figure BDA0003507442830000035
偶数段的信号片段长度为:n2=n-n1
将奇数段信号重新拼合成新的子信号:
Figure BDA0003507442830000032
将偶数段信号重新拼合成新的子信号:
Figure BDA0003507442830000033
在yo和ye上分别构造NRC函数,得到:
Figure BDA0003507442830000041
其中,
Figure BDA00035074428300000415
EnNRC统计量为
Figure BDA0003507442830000042
Figure BDA0003507442830000043
的加权和,记作
Figure BDA0003507442830000044
表达式为:
Figure BDA0003507442830000045
其中,ωo=onn1/l,ωe=enn2/l。
优选的,EnNRC统计量的表达式为:
Figure BDA0003507442830000046
定义l%n=0时yj,1为空,yj,2=yj,并且两个空向量的内积为0,对于所有i和j,
Figure BDA0003507442830000047
故对于任意l,EnNRC统计量统一为:
Figure BDA0003507442830000048
Figure BDA0003507442830000049
的表达式改写成矩阵形式:
Figure BDA00035074428300000410
其中,
Figure BDA00035074428300000411
是大小为l×l的矩阵:
Figure BDA00035074428300000412
其中,Mn是由On和En构成的块对角矩阵,
Figure BDA00035074428300000416
且:
Figure BDA00035074428300000413
Figure BDA00035074428300000414
在真实信号x和背景噪声∈上用同样方式划分信号并构造EnNRC统计量:
Figure BDA0003507442830000051
从而
Figure BDA0003507442830000052
改写为:
Figure BDA0003507442830000053
Figure BDA0003507442830000054
Figure BDA0003507442830000055
的均值和方差为:
Figure BDA0003507442830000056
Figure BDA0003507442830000057
Figure BDA0003507442830000058
Figure BDA0003507442830000059
Figure BDA00035074428300000510
的期望不大于真实信号x的平均能量
Figure BDA00035074428300000511
根据柯西-施瓦茨不等式:
Figure BDA00035074428300000512
其中,
Figure BDA00035074428300000513
当且仅当n∈N时等式成立;当l趋向于无穷时,
Figure BDA00035074428300000514
依概率收敛于
Figure BDA00035074428300000515
优选的,最高波峰出现在故障真实周期
Figure BDA00035074428300000522
及其倍数处的概率随着信号长度l→∞收敛到1,表达式为:
Figure BDA00035074428300000516
其中,
Figure BDA00035074428300000517
对任意n≤K及
Figure BDA00035074428300000518
随着l→∞,概率表达式为:
Figure BDA00035074428300000519
噪声项
Figure BDA00035074428300000520
对任意待检测波峰p和备选波峰q服从近似正态分布,表达式为:
Figure BDA00035074428300000521
Figure BDA0003507442830000061
Figure BDA0003507442830000062
时,所提出的假设检验统计量
Figure BDA0003507442830000063
与未知信号x有关;而对于任意p和q,波峰差
Figure BDA0003507442830000064
都是有界的,表达式为:
Figure BDA0003507442830000065
Figure BDA0003507442830000066
为故障的真实周期波峰;
Figure BDA0003507442830000067
为目标周期波峰,当且仅当q∈N时等号成立;
根据这一不等式,构造任意q值下的零假设的有界接受域:
Figure BDA0003507442830000068
其中,δ是控制接受域大小的参数;
构造任意q值下的零假设的缩减接受域为:
Figure BDA0003507442830000069
缩减接受域由以下公式得到,对所有p和q,有:
Figure BDA00035074428300000610
等价表示为:
Tr(CpCp+CqCq-2CpCq)≥Tr(CpCp-CqCq)
当且仅当p%q=0时等号成立,验证该不等式等同于验证不等式:
Tr(CpCq)≤Tr(CqCq)
当且仅当p%q=0时等号成立;
Figure BDA00035074428300000611
的每个列向量都是周期为n的周期序列,而p%q=0时,有:
Figure BDA00035074428300000612
因此:
Figure BDA00035074428300000613
所以当p%q=0时,有:
Tr((Cp-Cq)Cq)=Tr(-(Cp-Cq)Dq)=0
不等式中等式成立的条件为p%q=0,而p%q≠0时,得到:
Tr(CpCq)<Tr(CkCq)=Tr(CqCq)
其中,k是p和q的公倍数。
优选的,在接受域中存在多个显著波峰时,选择接受p值最大的
Figure BDA0003507442830000071
表达式为:
Figure BDA0003507442830000072
Figure BDA0003507442830000073
一直迭代到:
Figure BDA0003507442830000074
Figure BDA0003507442830000075
时,接受域
Figure BDA0003507442830000076
里不存在新的波峰
Figure BDA0003507442830000077
通过NRC算法中的噪声估计法对σ2进行估计,表达式为:
Figure BDA0003507442830000078
基于动态显著性水平设计迭代假设检验IHTDSL用于故障周期自动检测,IHTDSL的接受域
Figure BDA0003507442830000079
从故障信号中迭代学习,置信水平δ在每次迭代中根据下式动态估计:
Figure BDA00035074428300000710
其中,γ∈(0,1),
Figure BDA00035074428300000711
为本次假设检验的待检测波峰;
终止条件简化为:
Figure BDA00035074428300000712
通过估计IHTDSL中的δ值来动态控制显著性水平,每次迭代中的自适应接受域都在第一类错误和第二类错误之间实现平衡,使两者都收敛于0,表现为故障周期的误判和漏判都在信号足够长时消失。
根据本发明提供的轴承故障周期自动检测系统,包括:
模块M1:通过传感器从故障轴承上收集振动信号;
模块M2:将振动信号进行划分,将奇、偶片段分别重组,拼合成两段子信号;
模块M3:对每段子信号采用NRC算法处理,加权求和得到EnNRC统计量;
模块M4:利用该统计量进一步构造波峰差统计量;
模块M5:构造基于动态显著性水平的缩减置信区间;
模块M6:从EnNRC统计量的最高峰开始,向之前的波峰构造波峰差统计量,并构建动态置信区间,检查区间内是否有其他波峰,若存在其他波峰,则从区间内的最高波峰开始下一次假设检验,一直迭代到区间内不存在其他波峰为止,此时的波峰为算法检测到的故障真实周期。
优选的,将收集的故障轴承运转产生的周期脉冲信号y(t)分解为真实故障信号x(t)和背景噪声∈(t),则有:y(t)=x(t)+∈(t),其中,∈(t)~N(0,σ2);
定义时间序列y=[y(1),y(2),...,y(l)]T为y(t)的离散观测值,l为样本总数;x=[x(1),x(2),...,x(l)]T,∈=[∈(1),∈(2),...,∈(l)]T
则信号y表示为:y=x+∈,∈~N(0,σ2Il),Il为大小为l×l的单位矩阵;
设故障信号x(t)的周期为T0,T0=p0/fs,其中,p0>1为一个周期内的样本数目,fs为采样频率;
通过NRC算法将信号y划分为
Figure BDA0003507442830000089
段,第k个信号片段yk=[y((k-1)n+1),…,y(kn)],其中k=1,2,…,m;n是每段信号中样本数目;
则EnNRC统计量表示为:
Figure BDA0003507442830000081
其中,i、j为序列号;若l%n≠0,则信号被划分为2m+1段;
奇数段的信号片段长度为:n1=l-mn,其中
Figure BDA00035074428300000810
偶数段的信号片段长度为:n2=n-n1
将奇数段信号重新拼合成新的子信号:
Figure BDA0003507442830000082
将偶数段信号重新拼合成新的子信号:
Figure BDA0003507442830000083
在yo和ye上分别构造NRC函数,得到:
Figure BDA0003507442830000084
其中,
Figure BDA00035074428300000811
EnNRC统计量为
Figure BDA0003507442830000085
Figure BDA0003507442830000086
的加权和,记作
Figure BDA0003507442830000087
表达式为:
Figure BDA0003507442830000088
其中,ωo=onn1/l,ωe=enn2/l。
优选的,EnNRC统计量的表达式为:
Figure BDA0003507442830000091
定义l%n=0时yj,1为空,yj,2=yj,并且两个空向量的内积为0,对于所有i和j,
Figure BDA0003507442830000092
故对于任意l,EnNRC统计量统一为:
Figure BDA0003507442830000093
Figure BDA0003507442830000094
的表达式改写成矩阵形式:
Figure BDA0003507442830000095
其中,
Figure BDA0003507442830000096
是大小为l×l的矩阵:
Figure BDA0003507442830000097
其中,Mn是由On和En构成的块对角矩阵,
Figure BDA0003507442830000098
且:
Figure BDA0003507442830000099
Figure BDA00035074428300000910
在真实信号x和背景噪声∈上用同样方式划分信号并构造EnNRC统计量:
Figure BDA00035074428300000911
从而
Figure BDA00035074428300000912
改写为:
Figure BDA00035074428300000913
Figure BDA00035074428300000914
Figure BDA00035074428300000915
的均值和方差为:
Figure BDA00035074428300000916
Figure BDA00035074428300000917
Figure BDA0003507442830000101
Figure BDA0003507442830000102
Figure BDA0003507442830000103
的期望不大于真实信号x的平均能量
Figure BDA0003507442830000104
根据柯西-施瓦茨不等式:
Figure BDA0003507442830000105
其中,
Figure BDA0003507442830000106
当且仅当n∈N时等式成立;当l趋向于无穷时,
Figure BDA0003507442830000107
依概率收敛于
Figure BDA0003507442830000108
优选的,最高波峰出现在故障真实周期
Figure BDA00035074428300001023
及其倍数处的概率随着信号长度l→∞收敛到1,表达式为:
Figure BDA0003507442830000109
其中,
Figure BDA00035074428300001010
对任意n≤K及
Figure BDA00035074428300001011
随着l→∞,概率表达式为:
Figure BDA00035074428300001012
噪声项
Figure BDA00035074428300001013
对任意待检测波峰p和备选波峰q服从近似正态分布,表达式为:
Figure BDA00035074428300001014
Figure BDA00035074428300001015
Figure BDA00035074428300001016
时,所提出的假设检验统计量
Figure BDA00035074428300001017
与未知信号x有关;而对于任意p和q,波峰差
Figure BDA00035074428300001018
都是有界的,表达式为:
Figure BDA00035074428300001019
Figure BDA00035074428300001020
为故障的真实周期波峰;
Figure BDA00035074428300001021
为目标周期波峰,当且仅当q∈N时等号成立;
根据这一不等式,构造任意q值下的零假设的有界接受域:
Figure BDA00035074428300001022
其中,δ是控制接受域大小的参数;
构造任意q值下的零假设的缩减接受域为:
Figure BDA0003507442830000111
缩减接受域由以下公式得到,对所有p和q,有:
Figure BDA0003507442830000112
等价表示为:
Tr(CpCp+CqCq-2CpCq)≥Tr(CpCp-CqCq)
当且仅当p%q=0时等号成立,验证该不等式等同于验证不等式:
Tr(CpCq)≤Tr(CqCq)
当且仅当p%q=0时等号成立;
Figure BDA0003507442830000113
的每个列向量都是周期为n的周期序列,而p%q=0时,有:
Figure BDA0003507442830000114
因此:
Figure BDA0003507442830000115
所以当p%q=0时,有:
Tr((Cp-Cq)Cq)=Tr(-(Cp-Cq)Dq)=0
不等式中等式成立的条件为p%q=0,而p%q≠0时,得到:
Tr(CpCq)<Tr(CkCq)=Tr(CqCq)
其中,k是p和q的公倍数。
优选的,在接受域中存在多个显著波峰时,选择接受p值最大的
Figure BDA0003507442830000116
表达式为:
Figure BDA0003507442830000117
Figure BDA0003507442830000118
一直迭代到:
Figure BDA0003507442830000119
Figure BDA00035074428300001110
时,接受域
Figure BDA00035074428300001111
里不存在新的波峰
Figure BDA00035074428300001112
通过NRC算法中的噪声估计法对σ2进行估计,表达式为:
Figure BDA00035074428300001113
基于动态显著性水平设计迭代假设检验IHTDSL用于故障周期自动检测,IHTDSL的接受域
Figure BDA0003507442830000121
从故障信号中迭代学习,置信水平δ在每次迭代中根据下式动态估计:
Figure BDA0003507442830000122
其中,γ∈(0,1),
Figure BDA0003507442830000123
为本次假设检验的待检测波峰;
终止条件简化为:
Figure BDA0003507442830000124
通过估计IHTDSL中的δ值来动态控制显著性水平,每次迭代中的自适应接受域都在第一类错误和第二类错误之间实现平衡,使两者都收敛于0,表现为故障周期的误判和漏判都在信号足够长时消失。
与现有技术相比,本发明具有如下的有益效果:
(1)本发明通过动态置信区间平衡假设检验过程中的第一类、第二类错误,也能在信号足够长时消除故障周期的误报和漏报;
(2)本发明通过迭代假设检验比较存在干扰效应的周期性波峰及其他非周期波峰,逐步消除局部最优的误导效应,提高了周期检测的精度;
(3)本发明可以将该算法应用于其他本质为周期检测的应用。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述,本发明的其它特征、目的和优点将会变得更明显:
图1为算法实施流程图;
图2为ACF典型波形图;
图3为NRC典型波形图;
图4为EnNRC信号划分方式示意图;
图5为迭代假设检验第一步的有界接受域和缩减接受域示意图(δ=3);
图6为迭代假设检验最后一步的有界接受域和缩减接受域示意图(δ=3);
图7a和图7b分别为仿真信号x(n)和y(n)的示意图;
图8为轴承振动信号故障周期检测精度示意图(CWRU实验)。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明,但不以任何形式限制本发明。应当指出的是,对本领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干变化和改进。这些都属于本发明的保护范围。
实施例:
轴承出现单点故障时最明显的特征就是匀速转动时,振动信号会呈现出周期性特征,即故障有着一定的频率。该频率可以利用轴承的参数以及电机的转速计算得到。然而,由于设备的整体性和转速的不稳定性,且轴承随着使用时间的增加会产生一定磨损,实际参数会发生变化,实际应用中理论上的故障频率无法通过这种方式计算得到。因此,设计算法及时通过振动信号捕捉轴承故障周期,对轴承周期性故障成分的提取及后期的诊断是十分必要的。
故障轴承运转产生的周期脉冲信号y(t)可分解为真实故障信号x(t)和背景噪声∈(t),即:
y(t)=x(t)+∈(t)
其中,∈(t)~N(0,σ2)。
时间序列y=[y(1),y(2),...,y(Ll)]T为y(t)的离散观测值,也即振动传感器在轴承上能收集到的信号形式,其中l为样本总数。类似地,定义x=[x(1),x(2),...,x(l)]T,以及∈=[∈(1),∈(2),...,∈(l)]T。信号y可以表示为:
y=x+∈
其中,∈~N(0,σ2Il),Il为大小为l×l的单位矩阵。假设故障信号x(t)的周期为T0,T0=p0/fs,其中p0>1为一个周期内的样本数目,fs为采样频率。本发明还假设故障信号x的平均能量
Figure BDA0003507442830000134
是有界的。
大多数自相关类函数可以写为以下形式:
Figure BDA0003507442830000131
其中,Cn是大小为l×l的矩阵。如果Cn不是对称矩阵,则可以用
Figure BDA0003507442830000132
来替代Cn,故不失普遍性,可以假设Cn为对称矩阵。例如,Rabiner等人提出的ACF算法可表示为:
Figure BDA0003507442830000133
其中l≥n。该等式可被进一步改写成二次型形式:
Figure BDA0003507442830000141
Figure BDA0003507442830000142
Rn为l×l的矩阵,该矩阵在第n个非对角线处为单位元素,其余均为0。ACF典型波形图如图2所示。
NRC算法将信号y划分为
Figure BDA00035074428300001410
段,第k个信号片段yk=[y((k-1)n+1),…,y(kn)],其中k=1,2,…,m,n是每段信号中样本数目。NRC也是基于自相关算法设计的统计量,表示为:
Figure BDA0003507442830000143
可以看出,NRC统计量中并未考虑信号分段后最后一段长度为n1=l-mn的残余片段,故NRC一定程度上并未完全利用原始故障振动信号中的全部信息。类似地,NRC统计量可以改写为二次型:
Figure BDA00035074428300001411
Figure BDA0003507442830000144
其中,In为n×n的单位矩阵,0n为n×n的零矩阵,
Figure BDA0003507442830000145
为n1×n的零矩阵。NRC典型波形图如图3所示。
由于信号y的周期性特征,
Figure BDA0003507442830000146
的波形图通常在周期位置会出现明显的波峰,即
Figure BDA0003507442830000147
在n∈N,N={p0,2p0,3p0,…}时出现局部最优值。n∈N时
Figure BDA0003507442830000148
的波峰称为周期峰,而
Figure BDA0003507442830000149
时出现的波峰称为非周期峰。然而,这些周期峰值统计性质的相似性使得周期估计算法很容易出现故障误报和漏报,因为算法可能会被这些局部最优周期峰值所误导。因此,识别正确峰值所处的真实周期是相当困难的。具体而言,传统优化算法擅长通过以下方式找到搜索区域的最高峰值[1,K]:
Figure BDA0003507442830000151
其中,K>2p0
Figure BDA0003507442830000152
是最高峰
Figure BDA0003507442830000153
的坐标。在n≤≤K时,如果
Figure BDA0003507442830000154
有两个相同的最高峰,则返回坐标更小的峰值所在位置,即返回第一个最大峰值所在的位置。
按照这一方法找到的图2和3中的最高峰,并不是信号的真实周期所在位置,可能是信号真实周期的倍数。
在MLPE中提出的似然函数也可以写成二次型形式。这个似然函数的波形与NRC函数非常相似,呈现周期性的波峰。尽管使用了惩罚项来增强坐标较小的周期峰,MLPE统计量呈现的这些周期峰仍然具有相当大的误导性。因此,由于周期峰的误导作用,正确识别似然函数真实周期处的波峰也是相当困难的。
这种误导效应使得许多实际轴承故障检测应用中严重依赖领域知识和主观判断来识别故障周期,反馈的信息不够及时和可靠。为满足实际生产中对轴承故障周期准确而快速的估计而非人工参与周期识别的迫切需求,设计一种稳健的周期估计方法解决上述局部最优导致的的误导效应,从而确定真实周期,是至关重要和紧迫的。
简单地识别最大峰值可能会导致倍数周期错误,因为真正的峰值可能位于
Figure BDA0003507442830000155
之前的位置(即
Figure BDA0003507442830000156
)。要确定最高峰
Figure BDA0003507442830000157
是否为真实周期处的波峰,需要将该波峰与其之前的波峰进行比较,否则容易造成很高的误检率。因此,为了降低误检率,需要将最大波峰
Figure BDA0003507442830000158
与任何位于其之前的波峰
Figure BDA0003507442830000159
进行统计上的比较。
从统计上比较两个波峰可以利用统计假设检验完成。本发明提出以下假设检验:
Figure BDA00035074428300001521
来比较待检测波峰
Figure BDA00035074428300001510
与其他备选波峰
Figure BDA00035074428300001511
如果所有的
Figure BDA00035074428300001512
都明显小于
Figure BDA00035074428300001513
就可以认定
Figure BDA00035074428300001514
是故障的真实周期。如果存在显著接近
Figure BDA00035074428300001515
的波峰
Figure BDA00035074428300001516
则可以拒绝
Figure BDA00035074428300001517
位于真实周期处的假设,并重复假设检验过程以验证
Figure BDA00035074428300001518
是否为目标周期波峰。图2中两个标记的波峰在统计上可认为是显著相等的,这表明第一步搜索到的最大峰
Figure BDA00035074428300001519
可能不是正确的周期波峰。
本发明将对故障特征统计量两个可能位于周期处的备选波峰进行比较,以消除倍数周期波峰的误导效应。构造新的统计量:
Figure BDA00035074428300001520
即两个波峰的差值。基于这一统计量,可以采用迭代假设检验方法来从统计学上找到真实周期。
本发明拟设计Enhanced Noise Resistant Correlation(EnNRC)算法对轴承故障特征进行提取,这一步是为了初步削弱背景噪声的掩蔽效应,增强故障特征。它不仅具有与NRC相似的统计性质,还具有方便构造迭代假设检验的特殊统计性质。与NRC类似,EnNRC也将传感器收集到的振动信号划分为若干信号片段。如果信号能正好被划分为m段,每段信号有n个点(即l%n=0),EnNRC函数的构造方式与NRC完全相同。在这种情况下,EnNRC统计量可以表示为:
Figure BDA0003507442830000161
如果信号不能被恰好划分(即l%n≠0),则信号被划分为2m+1段。图4展示了这种情况下信号的具体划分方式。由图可知,奇数段的信号片段长度为n1=l-mn,其中
Figure BDA00035074428300001610
偶数段信号片段长度为n2=n-n1
将奇数段信号重新拼合成新的子信号:
Figure BDA0003507442830000162
而将偶数段信号重新拼合成:
Figure BDA0003507442830000163
在yo和ye上分别构造NRC函数,得到:
Figure BDA0003507442830000164
其中,
Figure BDA00035074428300001611
(不小于l/n的最小整数),
Figure BDA00035074428300001612
EnNRC统计量为
Figure BDA0003507442830000165
Figure BDA0003507442830000166
的加权和,记作
Figure BDA0003507442830000167
即:
Figure BDA0003507442830000168
其中,ωo=onn1/l,ωe=enn2/l。在该权重下,EnNRC会具有一些特殊的统计性质,将在下一节中详细讨论。
EnNRC统计量的完整形式为:
Figure BDA0003507442830000169
为了便于表示,定义l%n=0时yj,1为空,yj,2=yj,且两个空向量的内积为0,即对于所有i和j,
Figure BDA0003507442830000171
故对于任意l,EnNRC统计量可以统一写成:
Figure BDA0003507442830000172
Figure BDA0003507442830000173
的表达式也可以改写成矩阵形式:
Figure BDA0003507442830000174
其中,
Figure BDA0003507442830000175
是大小为l×l的矩阵:
Figure BDA0003507442830000176
其中Mn是由On和En构成的块对角矩阵,
Figure BDA0003507442830000177
且:
Figure BDA0003507442830000178
Figure BDA0003507442830000179
在真实信号x和背景噪声∈上也可用同样方式划分信号并构造EnNRC统计量:
Figure BDA00035074428300001710
从而
Figure BDA00035074428300001711
可以进一步写成:
Figure BDA00035074428300001712
Figure BDA00035074428300001713
Figure BDA00035074428300001714
的均值和方差为:
Figure BDA00035074428300001715
Figure BDA00035074428300001716
Figure BDA00035074428300001717
Figure BDA00035074428300001718
尽管比起Qn,Cn有着更复杂的矩阵结构,但EnNRC统计量与NRC统计量仍有着相似的统计性质。表明
Figure BDA00035074428300001719
的期望不大于真实信号x的平均能量
Figure BDA00035074428300001720
即根据柯西-施瓦茨不等式:
Figure BDA0003507442830000181
其中,
Figure BDA0003507442830000182
当且仅当n∈N时等式成立。当l趋向于无穷时(即l→∞),
Figure BDA0003507442830000183
依概率收敛于
Figure BDA0003507442830000184
与NRC统计量类似,EnNRC统计量提取的故障特征也会呈现一些显著突出的波峰。最高波峰出现在故障真实周期及其倍数处的概率随着信号长度l→∞收敛到1,即:
Figure BDA0003507442830000185
其中,
Figure BDA0003507442830000186
对任意n≤K及
Figure BDA0003507442830000187
随着l→∞,则有:
Figure BDA0003507442830000188
表明最大波峰是周期波峰之一,然而最大波峰可能不位于p0处,因为所有周期波峰都有着相同的期望值,即:
Figure BDA0003507442830000189
且背景噪声导致的方差可能造成最大波峰出现在任意位置。这种在统计上具有误导性的相似性质表明,仅仅定位最大的波峰而不与其之前的波峰进行比较,识别出的周期可信度较低,很容易导致故障周期的误报和漏报。
考虑到真实故障信号x是未知的,从统计上比较任何基于相关性函数的统计量的波峰是非常困难的(如
Figure BDA00035074428300001810
Figure BDA00035074428300001811
),这使得假设检验的解析拒绝域或接受域几乎不可能实现。然而,从另一个角度来看,EnNRC统计量可以很好地解决这一问题,因为对于任意p,q∈N,有:
Figure BDA00035074428300001812
这说明对于p,q∈N,
Figure BDA00035074428300001813
Figure BDA00035074428300001814
这两个波峰的差值只由背景噪声决定,而故障信号的背景噪声可以通过一定方法估算出来。因此,当p,q∈N时,仅研究
Figure BDA00035074428300001815
而非研究
Figure BDA00035074428300001816
就足以从统计上比较周期波峰。本发明将构造出一个特别的接受域,用于在假设检验中比较待检测波峰
Figure BDA00035074428300001817
和备选波峰
Figure BDA00035074428300001818
Figure BDA00035074428300001819
的精确概率分布可以用多种方法来计算,然而这些方法在l较大时都会出现插值效果差等问题。此外,为了克服背景噪声的掩蔽效应,在实际应用中一般都会收集足够长的信号来检测故障,这些方法如果要用数值积分来精确计算分布,则需要大量的计算时间,无法用于故障的在线检测。
本发明提出一种
Figure BDA0003507442830000191
的近似概率分布计算方法来解决这些难题。噪声项
Figure BDA0003507442830000192
对任意p和q服从近似正态分布,即:
Figure BDA0003507442830000193
Figure BDA0003507442830000194
然而,当
Figure BDA0003507442830000195
时,所提出的假设检验统计量
Figure BDA0003507442830000196
依然不可避免地与未知信号x有关,所以无法构造出准确的拒绝域或接受域。而对于任意p和q,波峰差
Figure BDA0003507442830000197
Figure BDA0003507442830000198
都是有界的,即:
Figure BDA0003507442830000199
当且仅当q∈N时等号成立。
根据这一不等式,可以构造任意q值下的零假设的有界接受域:
Figure BDA00035074428300001910
其中δ是控制接受域大小的参数。如图5所示,图中的区域即为假设检验的有界接受域。
与精确接受域相比,有界接受域具有一些特殊而重要的统计性质。首先,利用有界接受域进行假设检验时,当q∈N时,由于接受域大小不变,所以第一类错误和第二类错误率保持不变。其中第一类错误为
Figure BDA00035074428300001911
时错误地拒绝了零假设,而第二类错误为
Figure BDA00035074428300001912
时错误地接受了原假设。其次,当
Figure BDA00035074428300001913
时,由于有界接受域比精确接受域更窄,利用有界接受域进行假设检验会提高第一类错误率,降低第二类错误率。值得注意的是,由于
Figure BDA00035074428300001914
时出现的波峰都是非周期处的干扰波峰,故提高第一类错误率能有效地减少非周期波峰的接受域,即提高
Figure BDA00035074428300001915
的拒绝率反而能降低故障检测的误报率。因此,这两个性质表明,利用有界接受域进行假设检验有助于克服背景噪声污染导致的非周期波峰干扰,从而比精确接受域具有更高的故障周期检测精度。
然而,计算有界接受域非常耗费时间和内存,因为信号长度l通常相当大。因此,本研究进一步提出缩减接受域,如图5中的区域,以加快接受域的计算速度。任意q值下的零假设的缩减接受域为:
Figure BDA0003507442830000201
缩减接受域主要由以下公式得到:对所有p和q,有:
Figure BDA0003507442830000202
或等价表示为:
Tr(CpCp+CqCq-2CpCq)≥Tr(CpCp-CqCq)
当且仅当p%q=0时等号成立。验证该不等式等同于验证不等式:
Tr(CpCq)≤Tr(CqCq)
当且仅当p%q=0时等号成立。可以看出,
Figure BDA0003507442830000203
的每个列向量都是周期为n的周期序列,而p%q=0时,有:
Figure BDA0003507442830000204
因此:
Figure BDA0003507442830000205
所以当p%q=0时,有:
Tr((Cp-Cq)Cq)=Tr(-(Cp-Cq)Dq)=0
即不等式中等式成立的条件为p%q=0。而p%q≠0时,可以得到:
Tr(CpCq)<Tr(CkCq)=Tr(CqCq)
其中k是p和q的公倍数。该不等式成立是因为Ck的非零元素比Cp对应位置的元素更大。不等式的取等条件说明基于缩减接受域的假设检验在
Figure BDA0003507442830000206
时依然保持相同的第一类错误率和第二类错误率。这一结论意味着至少在q=p0时假设检验还保持着相同的误报率,因为
Figure BDA0003507442830000207
也就是说,使用缩减接受域时接受故障真实周期波峰的概率保持不变。使用缩减接受域的另一个好处是减少了其他不满足
Figure BDA0003507442830000208
的周期波峰
Figure BDA0003507442830000209
的接受率,一定程度上削弱了倍数周期波峰的误导效应。此外,与有界接受域类似,在
Figure BDA00035074428300002010
时,由于基于缩减接受域的假设检验使用了较窄的接受域,可以在很大程度上缓解背景噪声的掩蔽效应。如图5所示,有界接受域和缩减接受域均可以覆盖真实周期处的波峰,但缩减接受域覆盖的非周期波峰比有界接受域少的多,从而提高了判断和计算效率。
然而,接受域中可能不止存在一个
Figure BDA00035074428300002011
如图4所示,在接受域中共有3个显著波峰。在这种情况下,选择接受p值最大的
Figure BDA00035074428300002012
Figure BDA0003507442830000211
与通过简单优化指标(如似然性或p值)确定周期的传统方法类似,通过上式得到的
Figure BDA0003507442830000212
可能不是真实周期,这是由周期波峰的误导效应导致的。因此,这种定位统计量最大值来识别周期的一次性算法不能有效地克服误导效应。为了进一步解决这一问题,新找到的波峰
Figure BDA0003507442830000213
也应该继续与其之前的波峰进行比较,故需要重复上节提出的假设检验过程,将该波峰视为新的最大波峰,即令
Figure BDA0003507442830000214
一直迭代到:
Figure BDA0003507442830000215
也就是
Figure BDA0003507442830000216
时接受域
Figure BDA0003507442830000217
里不存在新的波峰
Figure BDA0003507442830000218
表1给出了迭代假设检验过程的详细步骤。图6为满足上式中终止条件时最后一次假设检验的接受区域。从图中可以看出,在最后找到的波峰之前,没有其他波峰位于接受域,则该波峰位于故障真实周期,该迭代假设检验过程可以有效克服倍数周期波峰的误导效应。
表1迭代假设检验算法流程
Figure BDA0003507442830000219
一般来说,确定置信区间时,参数δ都是固定的,如选择δ=3,因为待测统计量一般都服从正态分布。使用固定置信水平(δ值固定)的迭代假设检验方法被称为IHTFSL(Iterative Hypothesis Testing on Fixed Significance Level)。由于实际应用中背景噪声大小σ2是未知的,故参考NRC算法中的噪声估计方法来估计σ2
Figure BDA00035074428300002110
利用估计的σ2,执行IHTFSL时可计算出接受域。
然而,使用固定的显著性水平不可避免地会导致固定的第一类错误率,即使信号足够长。而当信号长度较短时,较小的显著水平会导致较高的第二类错误率,反之亦然。此外,同时减少第一类和第二类错误几乎是不可能的,因为在l保持不变的情况下,一种错误的减少几乎肯定会增加另一种错误。这些矛盾表明,在第一类错误和第二类错误之间找到更好的平衡可以有效地降低故障误报和漏报率。
为了克服这一矛盾,本发明将改进常用的固定置信区间的假设检验,基于动态显著性水平来设计迭代假设检验(IHTDSL)用于故障周期自动检测。IHTDSL的接受域
Figure BDA0003507442830000221
从故障信号中迭代学习,其中δ在每次迭代中根据以下公式动态估计:
Figure BDA0003507442830000222
其中γ∈(0,1),
Figure BDA0003507442830000223
为本次假设检验的待检测波峰。大量仿真和实验研究表明,在大多数情况下,q<p0时的非周期波峰
Figure BDA0003507442830000224
的期望值通常小于
Figure BDA0003507442830000225
的一半。因此,在迭代假设检验过程中选择γ≤0.5是完全合理的。仿真和实验研究结果还表明,γ=0.5时IHTDSL通常具有更高的故障周期检测精度。
终止条件可以简化为:
Figure BDA0003507442830000226
可以看出,此时的终止条件不再取决于未知的σ2,这是IHTDSL优于IHTFSL的第一个优点。第二个优点是随着信号长度趋于无穷,周期检测的第一类错误率和第二类错误率都会收敛于0。在传统的假设检验方法中增加数据样本的数量通常有助于减少第二类错误,而手动指定的第一类错误率保持不变。然而,通过估计IHTDSL中的δ值来动态控制显著性水平,每次迭代中的自适应接受域都能在第一类错误和第二类错误之间实现很好的平衡,使两者都收敛于0,表现为故障周期的误判和漏判都会在信号足够长时消失。
算法构建流程如图1所示,步骤如下:
1、利用传感器从故障轴承上收集振动信号y(t);
2、以图4所示的方式将信号进行再划分,将奇、偶片段分别重组,拼合成两段子信号;
3、对每段子信号采用NRC算法处理,加权求和得到EnNRC统计量;
4、利用该统计量进一步构造波峰差统计量;
5、构造基于动态显著性水平的缩减置信区间,可令γ=0.5;
6、从EnNRC统计量的最高峰开始,向之前的波峰构造波峰差统计量,并构建动态置信区间,检查区间内是否有其他波峰;
若存在其他波峰,则从区间内的最高波峰开始下一次假设检验,一直迭代到区间内不存在其他波峰为止,此时的波峰便可认为是算法检测到的故障真实周期。
优选地,本发明还包括:仿真验证;
本部分将所提出的IHTDSL和IHTFSL的性能进行比较,使用周期瞬态信号来模拟周期信号x:
Figure BDA0003507442830000231
其中ζ=0.01为阻尼比,f=10Hz为固有频率,p0=400为一个周期内样本点数目,fs=200Hz为采样频率。此外,还添加了强度为σ2=4的背景噪声,合成被背景噪声污染的信号y(n),波形如图7a和图7b所示。该信号的信噪比约为-19dB。
合成的仿真信号长度在10000~150000个点不等,每种长度的信号进行200次测试,且每种方法的周期识别区间为[1,2500],即K=2500。根据200次测试中周期检测正确次数
Figure BDA0003507442830000233
计算每种信号长度下的周期检测精度,结果如表2所示。
表2不同γ和l组合下IHTDSL、IHTFSL、NRCPE和MLPE周期检测精度
Figure BDA0003507442830000232
Figure BDA0003507442830000241
从表2结果可以看出,当γ取值较大时,IHTDSL具有较高的周期检测精度。随着信号长度的增加,所有γ取值下IHTDSL的精度均收敛到1,这进一步证明了:信号长度趋于无穷时误检率将收敛到0。周期检测精度收敛到1的速度随着γ值的增大而增大,收敛速度在γ≥0.5时不会继续增大,这也解释了为什么γ=0.5时周期检测精度收敛速度反而不如γ=0.4时快。一般来说,由于γ=0.5在大多数应用场景下的周期检测精度已足够适用,在使用IHTDSL方法时建议选择γ=0.5。
同样地,如果δ值保持不变,IHTFSL在处理更长的信号时具有更高的周期检测精度,这是因为增加信号长度可以有效减小第二类错误率。然而,IHTFSL使用固定的显著性水平,因此即使信号足够长,第一类错误率始终是固定不变的,如表2中所示,在δ=2时,在信号长度足够长的情况下,IHTFSL的周期检测精度依然在一定水平上下波动。这种波动会随着δ值的增加而减小,因为δ值的增大降低了第一类错误率。
当信号长度固定时,增大δ并不一定能提高周期检测精度。例如,当信号长度足够长(l≥30000)时,周期检测精度随着δ的增大而增大,因为在这种情况下第二类错误率远远小于第一类错误率。然而,在信号长度较短(l≤20000)时,第二类错误率相对较大,此时的周期检测精度不仅仅受第一类错误率支配,而是两者共同影响,随着δ的增大,周期检测精度先增大后减小,并没有呈现出明显的上升趋势。周期检测精度在不同δ取值上的非递增趋势表明,根据不同情况动态选择合适的δ,即自适应平衡第一类和第二类错误,可以降低周期检测的误报率。
实验验证——轴承故障周期检测
轴承是旋转机械的关键部件。当轴承发生故障时,从故障轴承采集到的信号会呈现周期性特征,因此可以对信号进行周期分析来检测故障。本实验使用CWRU轴承实验数据来比较IHTDSL、IHTFSL轴承故障周期检测效果,采样频率为48kHz,用于驱动端轴承实验。实验中使用的测试故障轴承为62052RS JEM SKF深沟球轴承。轴承外圈故障直径为0.007英寸,该单点故障通过电火花加工产生。轴承具体工作参数如表3所示。轴承故障周期经计算约为9.3ms,即一个周期内的样本点约为446个。
表3测试轴承62052RS JEM SKF工作参数
Figure BDA0003507442830000242
Figure BDA0003507442830000251
表3(续)测试轴承62052RS JEM SKF工作参数
Figure BDA0003507442830000252
本案例研究从采集到的振动信号中随机抽取长度为l的子信号,其中l=2500,3000,3500,…,6000。对每种长度随机抽取200条子信号来检测不同长度上IHTDSL和IHTFSL的周期检测精度。周期检测范围设置为[1,1200],即K=1200。每种方法的周期检测正确率(即满足
Figure BDA0003507442830000253
)如图8所示。
从图8可以看出,当待测信号长度较短时,所有方法的周期检测精度都比较低,主要由于短信号的背景噪声掩蔽效应较强,如果随机产生的非周期波峰具有与真实周期波峰相似的高度(甚至比真实周期波峰更高),那么识别真实周期波峰几乎是不可能的,从而导致周期误报。然而,在这种情况下,本章提出的IHTDSL方法仍然优于传统方法。
随着信号长度逐渐增大,γ=0.4和γ=0.5时的IHTDSL方法的周期检测精度稳步上升,且明显优于其他方法。相比之下,IHTFSL的性能很大程度取决于δ的选择。采用δ=3的IHTFSL精度最低,而当δ设置为30时,IHTFSL的精度得到了显著提高,说明简单地固定δ值并不能很好地适应具有不同噪声条件的信号。
IHTFSL的性能相对较差可能是因为所检测的信号可能不是严格周期性的,两个周期波峰之间的差异会受到周期波动的影响。因此,理论接受域太窄,不能覆盖真实周期波峰,导致了相当高的第一类错误率。这也可以解释为什么增大δ到30会大大提高周期检测的精度。这一结果表明,在假设检验中使用动态置信区间是十分必要的,因为在不同的应用背景下,手动确定合适的显著性水平相当困难。
本发明设计出能够对比并逐步排除局部最优值的统计方法,是降低轴承故障周期估计误诊率的重要思路和迫切需求。为了填补现有成果中对这一思路的研究空白,本技术将首先借用NRC的类似思路,设计增强抗噪声相关(EnNRC)函数来抵抗强背景噪声的掩蔽效应。在此基础上设计出一种迭代假设检验方法,逐步比较提出的EnNRC统计量的两个峰值。在迭代假设检验的每一步中动态计算接受区域,从而大大降低误报和漏报率。该方法通过迭代假设检验,动态地将备选最优值与其他局部最优值进行比较,最终排除局部最优的误导,消除故障周期误报和漏报,正确定位周期。
本发明提出了一种带有动态显著性水平的迭代假设检验方法,用于轴承单点故障信号的周期自动检测。首先提出了基于自相关函数设计的EnNRC函数用于提取信号特征。由于EnNRC特征信号波峰差的特殊统计性质,该统计量被作为假设检验的待测统计量,利用其统计性质特别设计原假设的接受域。为了充分克服周期波峰的误导效应,利用动态学习的显著性水平反复进行迭代假设检验,直到不再有明显的备选波峰。通过迭代假设检验,IHTDSL实现了第I类错误和第II类错误的平衡,并在信号足够长时消除了故障周期误报和漏报。
作为NRC的改进算法,基于EnNRC统计量设计的波峰差统计量和IHTDSL方法也同样适用于NRC算法,且由于自相关类算法统计量结构的相似性,通过简单适配则可将IHTDSL算法推广应用到这一类算法的改进。因此,本发明提出的IHTDSL具有一定的普适性,为自相关类特征提取算法的周期自动检测算法设计提供了理论支持和研究方向。
本领域技术人员知道,除了以纯计算机可读程序代码方式实现本发明提供的系统、装置及其各个模块以外,完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统、装置及其各个模块以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同程序。所以,本发明提供的系统、装置及其各个模块可以被认为是一种硬件部件,而对其内包括的用于实现各种程序的模块也可以视为硬件部件内的结构;也可以将用于实现各种功能的模块视为既可以是实现方法的软件程序又可以是硬件部件内的结构。
以上对本发明的具体实施例进行了描述。需要理解的是,本发明并不局限于上述特定实施方式,本领域技术人员可以在权利要求的范围内做出各种变化或修改,这并不影响本发明的实质内容。在不冲突的情况下,本申请的实施例和实施例中的特征可以任意相互组合。

Claims (10)

1.一种轴承故障周期自动检测方法,其特征在于,包括:
步骤1:通过传感器从故障轴承上收集振动信号;
步骤2:将振动信号进行划分,将奇、偶片段分别重组,拼合成两段子信号;
步骤3:对每段子信号采用NRC算法处理,加权求和得到EnNRC统计量;
步骤4:利用该统计量进一步构造波峰差统计量;
步骤5:构造基于动态显著性水平的缩减置信区间;
步骤6:从EnNRC统计量的最高峰开始,向之前的波峰构造波峰差统计量,并构建动态置信区间,检查区间内是否有其他波峰,若存在其他波峰,则从区间内的最高波峰开始下一次假设检验,一直迭代到区间内不存在其他波峰为止,此时的波峰为算法检测到的故障真实周期。
2.根据权利要求1所述的轴承故障周期自动检测方法,其特征在于,将收集的故障轴承运转产生的周期脉冲信号y(t)分解为真实故障信号x(t)和背景噪声∈(t),则有:y(t)=x(t)+∈(t),其中,∈(t)~N(0,σ2);
定义时间序列y=[y(1),y(2),...,y(l)]T为y(t)的离散观测值,l为样本总数;x=[x(1),x(2),...,x(l)]T,∈=[∈(1),∈(2),...,∈(l)]T
则信号y表示为:y=x+∈,∈~N(0,σ2Il),Il为大小为l×l的单位矩阵;
设故障信号x(t)的周期为T0,T0=p0/fs,其中,p0>1为一个周期内的样本数目,fs为采样频率;
通过NRC算法将信号y划分为
Figure FDA0003507442820000015
段,第k个信号片段yk=[y((k-1)n+1),…,y(kn)],其中k=1,2,…,m;n是每段信号中样本数目;
则EnNRC统计量表示为:
Figure FDA0003507442820000011
其中,i、j为序列号;若l%n≠0,则信号被划分为2m+1段;
奇数段的信号片段长度为:n1=l—mn,其中
Figure FDA0003507442820000014
偶数段的信号片段长度为:n2=n—n1
将奇数段信号重新拼合成新的子信号:
Figure FDA0003507442820000012
将偶数段信号重新拼合成新的子信号:
Figure FDA0003507442820000013
在yo和ye上分别构造NRC函数,得到:
Figure FDA0003507442820000021
其中,
Figure FDA0003507442820000022
EnNRC统计量为
Figure FDA0003507442820000023
Figure FDA0003507442820000024
的加权和,记作
Figure FDA0003507442820000025
表达式为:
Figure FDA0003507442820000026
其中,ωo=onn1/l,ωe=enn2/l。
3.根据权利要求2所述的轴承故障周期自动检测方法,其特征在于,EnNRC统计量的表达式为:
Figure FDA0003507442820000027
定义1%n=0时yj,1为空,yj,2=yj,并且两个空向量的内积为0,对于所有i和j,
Figure FDA0003507442820000028
故对于任意l,EnNRC统计量统一为:
Figure FDA0003507442820000029
Figure FDA00035074428200000215
的表达式改写成矩阵形式:
Figure FDA00035074428200000214
其中,
Figure FDA00035074428200000210
是大小为l×l的矩阵:
Figure FDA00035074428200000211
其中,Mn是由On和En构成的块对角矩阵,
Figure FDA00035074428200000216
且:
Figure FDA00035074428200000212
Figure FDA00035074428200000213
在真实信号x和背景噪声∈上用同样方式划分信号并构造EnNRC统计量:
Figure FDA0003507442820000031
从而
Figure FDA00035074428200000321
改写为:
Figure FDA0003507442820000032
Figure FDA0003507442820000033
Figure FDA0003507442820000034
的均值和方差为:
Figure FDA0003507442820000035
Figure FDA0003507442820000036
Figure FDA0003507442820000037
Figure FDA0003507442820000038
Figure FDA0003507442820000039
的期望不大于真实信号x的平均能量
Figure FDA00035074428200000310
限据柯西-施瓦茨不等式:
Figure FDA00035074428200000311
其中,
Figure FDA00035074428200000312
当且仅当n∈N时等式成立;当l趋向于无穷时,
Figure FDA00035074428200000313
依概率收敛于
Figure FDA00035074428200000314
4.根据权利要求3所述的轴承故障周期自动检测方法,其特征在于,最高波峰出现在故障真实周期
Figure FDA00035074428200000322
及其倍数处的概率随着信号长度l→∞收敛到1,表达式为:
Figure FDA00035074428200000315
其中,
Figure FDA00035074428200000316
对任意n≤K及
Figure FDA00035074428200000320
随着l→∞,概率表达式为:
Figure FDA00035074428200000317
噪声项
Figure FDA00035074428200000318
对任意待检测波峰p和备选波峰q服从近似正态分布,表达式为:
Figure FDA00035074428200000319
Figure FDA0003507442820000041
Figure FDA0003507442820000042
时,所提出的假设检验统计量
Figure FDA0003507442820000043
与未知信号x有关;而对于任意p和q,波峰差
Figure FDA0003507442820000044
都是有界的,表达式为:
Figure FDA0003507442820000045
Figure FDA0003507442820000046
为故障的真实周期波峰;
Figure FDA0003507442820000047
为目标周期波峰,当且仅当q∈N时等号成立;
根据这一不等式,构造任意q值下的零假设的有界接受域:
Figure FDA0003507442820000048
其中,δ是控制接受域大小的参数;
构造任意q值下的零假设的缩减接受域为:
Figure FDA0003507442820000049
缩减接受域由以下公式得到,对所有p和q,有:
Figure FDA00035074428200000410
等价表示为:
Figure FDA00035074428200000411
当且仅当p%q=0时等号成立,验证该不等式等同于验证不等式:
Tr(CpCq)≤Tr(CqCq)
当且仅当p%q=0时等号成立;
Figure FDA00035074428200000412
的每个列向量都是周期为n的周期序列,而p%q=0时,有:
Figure FDA00035074428200000413
因此:
Figure FDA00035074428200000414
所以当p%q=0时,有:
Tr((Cp—Cq)Cq)=Tr(-(Cp—Cq)Dq)=0
不等式中等式成立的条件为p%q=0,而p%q≠0时,得到:
Tr(CpCq)<Tr(CkCq)=Tr(CqCq)
其中,k是p和q的公倍数。
5.根据权利要求4所述的轴承故障周期自动检测方法,其特征在于,在接受域中存在多个显著波峰时,选择接受p值最大的
Figure FDA00035074428200000512
表达式为:
Figure FDA0003507442820000051
Figure FDA0003507442820000052
一直迭代到:
Figure FDA0003507442820000053
Figure FDA0003507442820000054
时,接受域
Figure FDA0003507442820000055
里不存在新的波峰
Figure FDA0003507442820000056
通过NRC算法中的噪声估计法对σ2进行估计,表达式为:
Figure FDA0003507442820000057
基于动态显著性水平设计迭代假设检验IHTDSL用于故障周期自动检测,IHTDSL的接受域
Figure FDA0003507442820000058
从故障信号中迭代学习,置信水平δ在每次迭代中根据下式动态估计:
Figure FDA0003507442820000059
其中,γ∈(0,1),
Figure FDA00035074428200000510
为本次假设检验的待检测波峰;
终止条件简化为:
Figure FDA00035074428200000511
通过估计IHTDSL中的δ值来动态控制显著性水平,每次迭代中的自适应接受域都在第一类错误和第二类错误之间实现平衡,使两者都收敛于0,表现为故障周期的误判和漏判都在信号足够长时消失。
6.一种轴承故障周期自动检测系统,其特征在于,包括:
模块M1:通过传感器从故障轴承上收集振动信号;
模块M2:将振动信号进行划分,将奇、偶片段分别重组,拼合成两段子信号;
模块M3:对每段子信号采用NRC算法处理,加权求和得到EnNRC统计量;
模块M4:利用该统计量进一步构造波峰差统计量;
模块M5:构造基于动态显著性水平的缩减置信区间;
模块M6:从EnNRC统计量的最高峰开始,向之前的波峰构造波峰差统计量,并构建动态置信区间,检查区间内是否有其他波峰,若存在其他波峰,则从区间内的最高波峰开始下一次假设检验,一直迭代到区间内不存在其他波峰为止,此时的波峰为算法检测到的故障真实周期。
7.根据权利要求6所述的轴承故障周期自动检测系统,其特征在于,将收集的故障轴承运转产生的周期脉冲信号y(t)分解为真实故障信号x(t)和背景噪声∈(t),则有:y(t)=x(t)+∈(t),其中,∈(t)~N(0,σ2);
定义时间序列y=[y(1),y(2),...,y(l)]T为y(t)的离散观测值,l为样本总数;x=[x(1),x(2),...,x(l)]T,∈=[∈(1),∈(2),...,∈(l)]T
则信号y表示为:y=x+∈,∈~N(0,σ2Il),Il为大小为l×l的单位矩阵;
设故障信号x(t)的周期为T0,T0=p0/fs,其中,p0>1为一个周期内的样本数目,fs为采样频率;
通过NRC算法将信号y划分为
Figure FDA00035074428200000611
段,第k个信号片段yk=[y((k-1)n+1),…,y(kn)],其中k=1,2,…,m;n是每段信号中样本数目;
则EnNRC统计量表示为:
Figure FDA0003507442820000061
其中,i、j为序列号;若1%n≠0,则信号被划分为2m+1段;
奇数段的信号片段长度为:n1=l-mn,其中
Figure FDA00035074428200000610
偶数段的信号片段长度为:n2=n—n1
将奇数段信号重新拼合成新的子信号:
Figure FDA0003507442820000062
将偶数段信号重新拼合成新的子信号:
Figure FDA0003507442820000063
在yo和ye上分别构造NRC函数,得到:
Figure FDA0003507442820000064
其中,
Figure FDA0003507442820000065
EnNRC统计量为
Figure FDA0003507442820000066
Figure FDA0003507442820000067
的加权和,记作
Figure FDA0003507442820000068
表达式为:
Figure FDA0003507442820000069
其中,ωo=onn1/l,ωe=enn2/l。
8.根据权利要求7所述的轴承故障周期自动检测系统,其特征在于,EnNRC统计量的表达式为:
Figure FDA0003507442820000071
定义l%n=0时yj,1为空,yj,2=yj,并且两个空向量的内积为0,对于所有i和j,
Figure FDA0003507442820000072
故对于任意l,EnNRC统计量统一为:
Figure FDA0003507442820000073
Figure FDA0003507442820000074
的表达式改写成矩阵形式:
Figure FDA0003507442820000075
其中,
Figure FDA0003507442820000076
是大小为l×l的矩阵:
Figure FDA0003507442820000077
其中,Mn是由On和En构成的块对角矩阵,
Figure FDA00035074428200000715
且:
Figure FDA0003507442820000078
Figure FDA0003507442820000079
在真实信号x和背景噪声∈上用同样方式划分信号并构造EnNRC统计量:
Figure FDA00035074428200000710
从而
Figure FDA00035074428200000711
改写为:
Figure FDA00035074428200000712
Figure FDA00035074428200000713
Figure FDA00035074428200000714
的均值和方差为:
Figure FDA0003507442820000081
Figure FDA0003507442820000082
Figure FDA0003507442820000083
Figure FDA0003507442820000084
Figure FDA0003507442820000085
的期望不大于真实信号x的平均能量
Figure FDA0003507442820000086
根据柯西-施瓦茨不等式:
Figure FDA0003507442820000087
其中,
Figure FDA0003507442820000088
当且仅当n∈N时等式成立;当l趋向于无穷时,
Figure FDA0003507442820000089
依概率收敛于
Figure FDA00035074428200000810
9.根据权利要求8所述的轴承故障周期自动检测系统,其特征在于,最高波峰出现在故障真实周期
Figure FDA00035074428200000824
及其倍数处的概率随着信号长度l→∞收敛到1,表达式为:
Figure FDA00035074428200000811
其中,
Figure FDA00035074428200000812
对任意n≤K及
Figure FDA00035074428200000813
随着l→∞,概率表达式为:
Figure FDA00035074428200000814
噪声项
Figure FDA00035074428200000815
对任意待检测波峰p和备选波峰q服从近似正态分布,表达式为:
Figure FDA00035074428200000816
Figure FDA00035074428200000817
Figure FDA00035074428200000818
时,所提出的假设检验统计量
Figure FDA00035074428200000819
与未知信号x有关;而对于任意p和q,波峰差
Figure FDA00035074428200000820
都是有界的,表达式为:
Figure FDA00035074428200000821
Figure FDA00035074428200000822
为故障的真实周期波峰;
Figure FDA00035074428200000823
为目标周期波峰,当且仅当q∈N时等号成立;
根据这一不等式,构造任意q值下的零假设的有界接受域:
Figure FDA0003507442820000091
其中,δ是控制接受域大小的参数;
构造任意q值下的零假设的缩减接受域为:
Figure FDA0003507442820000092
缩减接受域由以下公式得到,对所有p和q,有:
Figure FDA0003507442820000093
等价表示为:
Figure FDA0003507442820000094
当且仅当p%q=0时等号成立,验证该不等式等同于验证不等式:
Figure FDA0003507442820000095
当且仅当p%q=0时等号成立;
Figure FDA0003507442820000096
的每个列向量都是周期为n的周期序列,而p%q=0时,有:
Figure FDA0003507442820000097
因此:
Figure FDA0003507442820000098
所以当p%q=0时,有:
Figure FDA0003507442820000099
不等式中等式成立的条件为p%q=0,而p%q≠0时,得到:
Tr(CpCq)<Tr(CkCq)=Tr(CqCq)
其中,k是p和q的公倍数。
10.根据权利要求9所述的轴承故障周期自动检测系统,其特征在于,在接受域中存在多个显著波峰时,选择接受p值最大的
Figure FDA00035074428200000910
表达式为:
Figure FDA00035074428200000911
Figure FDA00035074428200000912
一直迭代到:
Figure FDA0003507442820000101
Figure FDA0003507442820000109
时,接受域
Figure FDA0003507442820000102
里不存在新的波峰
Figure FDA0003507442820000103
通过NRC算法中的噪声估计法对σ2进行估计,表达式为:
Figure FDA0003507442820000104
基于动态显著性水平设计迭代假设检验IHTDSL用于故障周期自动检测,IHTDSL的接受域
Figure FDA0003507442820000105
从故障信号中迭代学习,置信水平δ在每次迭代中根据下式动态估计:
Figure FDA0003507442820000106
其中,γ∈(0,1),
Figure FDA0003507442820000107
为本次假设检验的待检测波峰;
终止条件简化为:
Figure FDA0003507442820000108
通过估计IHTDSL中的δ值来动态控制显著性水平,每次迭代中的自适应接受域都在第一类错误和第二类错误之间实现平衡,使两者都收敛于0,表现为故障周期的误判和漏判都在信号足够长时消失。
CN202210143261.3A 2022-02-16 2022-02-16 轴承故障周期自动检测方法和系统 Pending CN114528879A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210143261.3A CN114528879A (zh) 2022-02-16 2022-02-16 轴承故障周期自动检测方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210143261.3A CN114528879A (zh) 2022-02-16 2022-02-16 轴承故障周期自动检测方法和系统

Publications (1)

Publication Number Publication Date
CN114528879A true CN114528879A (zh) 2022-05-24

Family

ID=81622767

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210143261.3A Pending CN114528879A (zh) 2022-02-16 2022-02-16 轴承故障周期自动检测方法和系统

Country Status (1)

Country Link
CN (1) CN114528879A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116484308A (zh) * 2023-06-25 2023-07-25 火眼科技(天津)有限公司 一种基于边缘自适用计算的数据采集方法
CN116827055A (zh) * 2022-09-09 2023-09-29 东莞市智美生活电子科技有限公司 一种电机结构

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116827055A (zh) * 2022-09-09 2023-09-29 东莞市智美生活电子科技有限公司 一种电机结构
CN116827055B (zh) * 2022-09-09 2024-01-30 东莞市智美生活电子科技有限公司 一种电机结构
CN116484308A (zh) * 2023-06-25 2023-07-25 火眼科技(天津)有限公司 一种基于边缘自适用计算的数据采集方法
CN116484308B (zh) * 2023-06-25 2023-09-29 火眼科技(天津)有限公司 一种基于边缘自适用计算的数据采集方法

Similar Documents

Publication Publication Date Title
CN110132598B (zh) 旋转设备滚动轴承故障噪声诊断算法
Chen et al. Multi-fault condition monitoring of slurry pump with principle component analysis and sequential hypothesis test
CN109187025B (zh) 一种集成kelm的滚动轴承剩余使用寿命预测方法
CN114528879A (zh) 轴承故障周期自动检测方法和系统
Chen et al. Frequency-temporal-logic-based bearing fault diagnosis and fault interpretation using Bayesian optimization with Bayesian neural networks
Yang et al. Bearing fault automatic classification based on deep learning
CN110610035A (zh) 一种基于gru神经网络的滚动轴承剩余寿命预测方法
Jia et al. GTFE-Net: A gramian time frequency enhancement CNN for bearing fault diagnosis
US20230316051A1 (en) Pre-alarming method for rotary stall of compressors based on temporal dilated convolutional neural network
CN116226646B (zh) 轴承健康状态及剩余寿命的预测方法、系统、设备及介质
CN111881594B (zh) 一种核动力设备的非平稳信号状态监测方法及系统
CN113342597B (zh) 一种基于高斯混合隐马尔可夫模型的系统故障预测方法
Jiang et al. Dual-attention-based multiscale convolutional neural network with stage division for remaining useful life prediction of rolling bearings
CN107480386A (zh) 一种基于响应混叠性度量与遗传算法的测试激励优选方法
CN114417699A (zh) 泵阀故障检测方法
Sadoughi et al. A deep learning approach for failure prognostics of rolling element bearings
CN113256443B (zh) 核电水泵导轴承故障检测方法、系统、设备及可读存储介质
CN114064459A (zh) 基于生成对抗网络和集成学习的软件缺陷预测方法
Lv et al. A new feature extraction technique for early degeneration detection of rolling bearings
Yu et al. Rolling bearing degradation state identification based on LCD relative spectral entropy
CN110222390B (zh) 基于小波神经网络的齿轮裂纹识别方法
CN116758922A (zh) 一种用于变压器的声纹监测与诊断方法
Qi et al. A new deep fusion network for automatic mechanical fault feature learning
CN113051809A (zh) 一种基于改进受限玻尔兹曼机的虚拟健康因子构建方法
CN116625685B (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