CN115683629B - 一种轴承故障检测方法、装置及设备 - Google Patents

一种轴承故障检测方法、装置及设备 Download PDF

Info

Publication number
CN115683629B
CN115683629B CN202211398933.1A CN202211398933A CN115683629B CN 115683629 B CN115683629 B CN 115683629B CN 202211398933 A CN202211398933 A CN 202211398933A CN 115683629 B CN115683629 B CN 115683629B
Authority
CN
China
Prior art keywords
matrix
time
function
objective function
frequency domain
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.)
Active
Application number
CN202211398933.1A
Other languages
English (en)
Other versions
CN115683629A (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.)
Suzhou University
Original Assignee
Suzhou 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 Suzhou University filed Critical Suzhou University
Priority to CN202211398933.1A priority Critical patent/CN115683629B/zh
Publication of CN115683629A publication Critical patent/CN115683629A/zh
Application granted granted Critical
Publication of CN115683629B publication Critical patent/CN115683629B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及一种轴承故障检测方法包括：获取测试轴承的振动加速度时域信号，经短时傅里叶变换生成时频域系数矩阵；利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型的凸性目标函数；利用交替方向乘子法将凸性目标函数拆分成关于广义极小极大凹罚函数的第一子目标函数、关于截断核范数的第二子目标函数与辅助迭代函数；利用前向后向分裂算法求解第一子目标函数得关于Xk+1的迭代公式组；利用奇异值阈值算法求解第二目标子函数得关于Zk+1的迭代公式；初始化参数，将时频域系数矩阵输入辅助迭代函数中，迭代预设次数获取时频域稀疏低秩矩阵；对时频域稀疏低秩矩阵进行短时傅里叶逆变换得到重构时域信号后，进行平方包络谱分析，得到故障特征频率。

Description

一种轴承故障检测方法、装置及设备
技术领域
本发明涉及机械设备故障检测技术领域，尤其是指一种轴承故障检测方法、装置及设备。
背景技术
机械系统在国民经济发展中发挥着重要的作用，而轴承在机械系统的动力传输中又起着关键的作用，其健康状态会关系到整个设备的运行性能与安全。由于长期运行在复杂工况及环境下，轴承容易受到损坏，从而造成设备故障，由此可能引起重大的经济损失甚至人员伤亡。因此，对轴承进行准确、及时的故障诊断对于机械系统运行的可靠性具有重要意义。由于机械故障容易带来其振动信号上的变化，所以基于振动信号分析机械设备状态的方式已广泛应用于机械故障诊断领域。当轴承发生故障时，振动信号会激发脉冲响应，呈现周期性和非平稳性；但与此同时，复杂的工作环境也可能会引入大量的噪声，导致故障相关分量被噪声淹没，使故障特征不易被识别。
稀疏低秩矩阵估计是一种常用于图像识别与视频监控的去噪方法，其在这两方面的优异性能受到了广泛关注。从低秩矩阵的角度来看，信号矩阵的秩是矩阵行和列之间相关性的度量，对信号矩阵的秩优化可以有限地去除信号中的冗余信息，从而挖掘出信号中包含的更深层次的信息。由轴承故障引起的瞬态脉冲不仅在时域具有稀疏性，相应的，其时频域系数同样具有稀疏性，同时时频域系数的奇异值分布显示其具有低秩特性，而噪声信号的时频域系数不具备这两种特性。因此，稀疏低秩矩阵估计的方法同样适用于故障特征提取，其基本思想是：将一维含有噪声的信号通过短时傅里叶变换得到其时频域的系数矩阵，将从一维含噪信号中提取轴承故障脉冲信号的问题转化为一个从二维时频域系数矩阵中估计一个稀疏低秩矩阵的问题，再通过对所得到的稀疏低秩矩阵进行短时傅里叶逆变换，从而得到想要的故障瞬态脉冲。然而，由于核范数不能准确的逼近矩阵的秩，稀疏诱导范数不能准确地逼近L0范数，稀疏低秩矩阵估计的方法在故障诊断领域中的应用未取得广泛的发展。
传统稀疏低秩矩阵估计方法通过L1范数诱导时频系数矩阵的稀疏性，通过核范数来逼近矩阵的秩。由于L1范数本身的缺点，其保留信号幅值的能力差，使得最终得到的故障瞬态信号被明显低估。并且，核范数难以非常准确的逼近矩阵的秩，这导致其对于大的奇异值的保留结果会有较大的低估，小的奇异值没有及时的被舍去，由此也会影响最终得到的故障瞬态信号的幅值以及去噪效果。因此，现有的稀疏低秩矩阵估计方法稀疏诱导范数不能准确地逼近L0范数，容易引起故障瞬态脉冲的幅值低估；核范数不能准确地逼近矩阵的秩，奇异值保留效果差，影响故障瞬态脉冲的恢复效果；所得到的恢复信号幅值低估较为严重，平方包络谱幅值低，故障特征频率不明显。
发明内容
为此，本发明所要解决的技术问题在于克服现有技术中对轴承故障瞬态脉冲的重构时域信号幅值低估严重的问题。
为解决上述技术问题，本发明提供了一种轴承故障检测方法，包括：
S1：获取测试轴承的振动加速度时域信号，利用短时傅里叶变换生成时频域系数矩阵Y；
S2：利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型的目标函数，并根据其凸性条件转换为凸性目标函数；利用交替方向乘子法将所述凸性目标函数拆分成关于广义极小极大凹罚函数的第一子目标函数、关于截断核范数的第二子目标函数与辅助迭代函数；
S3：利用前向后向分裂算法，求解所述第一子目标函数，得到关于Xk+1的迭代公式组；
S4：根据所述关于Xk+1的迭代公式组，利用奇异值阈值算法求解所述第二目标子函数，得到关于Zk+1的迭代公式；
S5：利用关于Xk+1的迭代公式组与关于Zk+1的迭代公式更新所述辅助迭代函数，初始化参数，将所述时频域系数矩阵Y输入所述辅助迭代函数中，迭代预设次数Nit次获取对应的时频域稀疏低秩矩阵XNit；
S6：对所述时频域稀疏低秩矩阵XNit进行短时傅里叶逆变换得到轴承故障瞬态脉冲的重构时域信号；
S7：对所述轴承故障瞬态脉冲的重构时域信号进行平方包络谱分析，得到轴承故障特征频率。
在本发明的一个实施例中，所述利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型的目标函数，并根据其凸性条件转换为凸性目标函数包括：
所述稀疏低秩优化模型的目标函数：
Figure BDA0003934504310000031
所述凸性条件为ITI-λ2HTH≥0，
Figure BDA0003934504310000032
0≤γ≤1；
所述凸性目标函数：
Figure BDA0003934504310000033
其中，X初始化为与矩阵Y规格一致的零矩阵；L表示矩阵X的近似矩阵，X＝Rm×n；||X||F代表X的F范数，表示为
Figure BDA0003934504310000034
||X||1代表X的1范数，表示为
Figure BDA0003934504310000035
Tr(X)表示矩阵X的迹，第一辅助矩阵A与第二辅助矩阵B满足A∈Rm×n、BT∈Rm×n、AAT＝I、BBT＝I，I表示单位矩阵；
其中，将矩阵X的奇异值分解定义为X＝UΣVT：
Figure BDA0003934504310000041
Ak＝(Uk(:,[1:r]))T，Bk＝(Vk(:,[1:r]))T；
U＝(u1,…,um)∈Rm×m，Σ∈Rm×n，V＝(v1,…,vm)∈Rn×n；当第一辅助矩阵A与第二辅助矩阵B满足A＝(u1,…,ur)T，B＝(v1,…,vr)T时，存在||X||T＝||X||*-Tr(AXBT)，σi表示时频域系数矩阵X的第i个奇异值；r为截断阶数，表示将被保留的最大的r个非零奇异值；H表示将所述广义极小极大凹罚函数参数化的参数矩阵，满足HTH＝nI；λ1为第一正则化参数、λ2为第二正则化参数，γ为常量参数。
在本发明的一个实施例中，所述利用交替方向乘子法将所述凸性目标函数拆分成关于广义极小极大凹罚函数的第一子目标函数与关于截断核范数的第二子目标函数包括：
引入变量Z，Z满足Z＝X，带入所述凸性目标函数中得：
Figure BDA0003934504310000042
所述关于广义极小极大凹罚函数的第一子目标函数：
Figure BDA0003934504310000043
所述关于截断核范数的第二子目标函数：
Figure BDA0003934504310000044
辅助迭代函数：R＝R-(X-Z)；
其中，μ表示迭代步长，辅助矩阵R初始化为与所述时频域系数矩阵Y规格一致的零矩阵。
在本发明的一个实施例中，所述利用前向后向分裂算法，求解所述第一子目标函数，得到关于Xk+1的迭代公式组包括：
Figure BDA0003934504310000045
分别对X和L求偏导得到：
Figure BDA0003934504310000046
Figure BDA0003934504310000051
其中，符号函数
Figure BDA0003934504310000052
将对X和L求偏导得到的公式带入前向后向分裂算法中得关于Xk+1的迭代公式组：
Qk＝Xk-ξ((1+μ)Xk-Y-μ(Zk+Rk)-γ(Xk+Lk))；
Xk+1＝soft(Qk,λ2ξ)；
Tk＝Lk+ξγ(Xk+1-Lk)；
Lk+1＝soft(Tk,λ2ξ)；
其中，辅助矩阵R是初始化为与所述时频域系数矩阵Y规格一致的零矩阵，第一矩阵Qk为Xk梯度下降得到的矩阵，第二矩阵Tk为近似矩阵Lk梯度下降得到的矩阵,ξ表示步长，k表示迭代次数，soft代表软阈值算子，soft(x,C)＝x·max{0,1-C/|x|}。
在本发明的一个实施例中，所述根据所述关于Xk+1的迭代公式组，利用奇异值阈值算法求解所述第二目标子函数，得到关于Zk+1的迭代公式包括：
将
Figure BDA0003934504310000053
转化为：
Figure BDA0003934504310000054
利用奇异值阈值算法求解，得关于Zk+1的迭代公式：
Figure BDA0003934504310000055
Figure BDA0003934504310000056
其中，Ak为第一辅助矩阵，Bk为第二辅助矩阵，Pk为左奇异矩阵，
Figure BDA0003934504310000064
为右奇异矩阵，∑k为非负实对角矩阵。
在本发明的一个实施例中，所述利用奇异值阈值算法求解所述第二目标子函数，得到关于Zk+1的迭代公式还包括：
构建权值向量W与所述奇异值阈值算法中的阈值
Figure BDA0003934504310000061
相乘对阈值进行加权；
其中，权值向量W中的元素定义为：
Figure BDA0003934504310000062
σi＝∑(i,i)，σi表示奇异值，ε为常量常数；
将权值向量W引入所述Z的迭代公式中，更新的关于Zk+1的迭代公式表示为：
Figure BDA0003934504310000063
在本发明的一个实施例中，所述利用关于Xk+1的迭代公式组与关于Zk+1的迭代公式更新所述辅助迭代函数，初始化参数，将所述时频域系数矩阵Y输入所述辅助迭代函数中，迭代预设次数Nit次获取对应的时频域稀疏低秩矩阵XNit：
S51：利用关于Xk+1的迭代公式组与关于Zk+1的迭代公式更新所述辅助迭代函数：Rk+1＝Rk-(Xk+1-Zk+1)；
S52：设置迭代次数k＝0,1,2,3,…,Nit，输入时频域系数矩阵Y，初始化迭代步长μ、步长ξ、第一正则化参数λ1与第二正则化参数λ2、截断阶数r、常量参数γ；
S53：初始化迭代次数k＝0，初始化矩阵Xk、变量Zk、辅助矩阵Rk与近似矩阵Lk为与所述时频域系数矩阵Y规格一致的零矩阵：X0、Z0、R0与L0；
对X0进行梯度下降得到的第一矩阵Q0，利用软阈值函数，根据Q0迭代求出X1；对L0进行梯度下降得到第二矩阵T0，利用软阈值函数，根据T0迭代求解出L1；奇异值阈值算法求解得到Z1；根据R0、X1与Z1求解出R1；根据所述X1、L1、Z1与R1更新得到Q1；
S54：对于k＝1,2,3,…,Nit；
根据Xk、Zk、Rk与Lk更新得到Qk，利用软阈值函数，根据Qk求解Xk+1；
根据Xk+1更新Tk，利用软阈值函数，根据Tk更新Lk+1；
利用奇异值算法，根据Rk与Xk+1更新Zk+1；
根据Rk、Xk+1与Zk+1更新得到Rk+1；
利用Xk+1、Zk+1、Rk+1与Lk+1更新得到Qk+1；
S55：令k＝k+1，重复S54，直至k+1＝Nit，停止迭代，输出XNit。
本发明还提供了一种轴承故障检测装置，包括：
轴承故障信号预处理模块，用于获取测试轴承的振动加速度时域信号，利用短时傅里叶变换生成时频域系数矩阵；
稀疏低秩优化处理模块，用于利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型的目标函数，并利用交替方向乘子法求解获取目标函数迭代公式，将所述时频域系数矩阵输入所述目标函数迭代公式，迭代获取对应的时频域稀疏低秩矩阵；
时域信号重构模块，用于对所述时频域稀疏低秩矩阵进行短时傅里叶逆变换得到轴承故障瞬态脉冲的重构时域信号；
包络谱分析模块，用于对所述轴承故障瞬态脉冲的重构时域信号进行平方包络谱分析，得到轴承故障特征频率。
本发明还提供了一种轴承故障检测设备，包括：
采集装置，用于利用加速度传感器采集测试轴承的加速度时域信号；
轴承故障检测装置，与所述采集装置通讯连接，用于执行计算机程序实现如上述所述的轴承故障检测方法；
显示装置，与所述轴承故障检测装置通讯连接，用于显示所述重构时域信号的平方包络谱图。
在本发明的一个实施例中，所述采集装置包括：
驱动电机；
转动轴，一端通过联轴器连接所述驱动电机；
测试轴承安装座，套设于所述转动轴上，用于固定测试轴承；
加速度传感器，设置于所述转动轴的另一端，用于连接所述测试轴承并测量其振动加速度；
负载装置，用于添加可调节的负载至所述转动轴。
本发明的上述技术方案相比现有技术具有以下优点：
本发明所述的轴承故障检测方法将一维测试轴承振动加速度的时域信号进行短时傅里叶变换得到二维时频域系数矩阵；对二维时频域系数矩阵进行稀疏低秩优化处理生成稀疏低秩矩阵；对稀疏低秩矩阵进行短时傅里叶逆变换得到轴承故障瞬态脉冲的重构时域信号；对重构时域信号进行包络谱分析，得到故障特征频率。在进行稀疏低秩优化处理时，对于时频域系数矩阵的稀疏性，设计广义极小极大凹罚函数逼近L0范数，提升稀疏低秩优化模型的保幅值能力与去噪效果；对于时频域系数矩阵的低秩性，利用截断核范数逼近时频域系数矩阵的秩，并对传统奇异值阈值算法进行改进，设计了与奇异值大小呈负相关的权值向量，利用权值向量对奇异值阈值进行加权处理，提升了奇异值保留效果，从而提升了轴承故障检测的准确性，可信度高。
附图说明
为了使本发明的内容更容易被清楚的理解，下面根据本发明的具体实施例并结合附图，对本发明作进一步详细的说明，其中
图1是本发明实施例所提供的轴承故障检测方法流程图；
图2是本发明实施例所提供的轴承故障检测设备的采集装置示意图；
图3是本发明实施例所提供的故障轴承振动加速度的时域信号图；
图4是本发明实施例所提供的故障轴承振动加速度时域信号经过短时傅里叶变换得到的时频图；
图5是本发明实施例所提供的故障轴承振动加速度时域信号的时频图；
图6是本发明实施例所提供的轴承故障瞬态脉冲的重构时域信号的时频图；
图7是本发明实施例所提供的故障轴承振动加速度的时频域系数矩阵与轴承故障瞬态脉冲的重构时域信号的时频系数矩阵的奇异值分布对比图；
图8是本发明实施例所提供的轴承故障瞬态脉冲的重构时域信号图；
图9的(a)为本发明实施例所提供的轴承故障瞬态脉冲信号的平方包络谱图，图9的(b)为本发明实施例所提供的轴承故障瞬态脉冲信号的平方包络谱局部图。
具体实施方式
下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好地理解本发明并能予以实施，但所举实施例不作为对本发明的限定。
参照图1所示，本发明的轴承故障检测方法包括：
S1：获取测试轴承的振动加速度时域信号，利用短时傅里叶变换生成时频域系数矩阵Y：Y＝S(y)；
其中，Y表示时频域系数矩阵；y表示振动加速度时域信号；S表示短时傅里叶变换，其具体参数设计可根据使用需求进行调节；
S2：利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型的目标函数，并根据其凸性条件转换为凸性目标函数；利用交替方向乘子法将所述凸性目标函数拆分成关于广义极小极大凹罚函数的第一子目标函数、关于截断核范数的第二子目标函数与辅助迭代函数；
S21：利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型的目标函数，并根据其凸性条件转换为凸性目标函数：
广义极小极大凹罚函数：
Figure BDA0003934504310000101
截断核范数：
Figure BDA0003934504310000102
矩阵X的奇异值分解定义为X＝UΣVT，其中，U＝(u1,…,um)∈Rm×m，Σ∈Rm×n，V＝(v1,…,vm)∈Rn×n；当辅助矩阵A与B满足A＝(u1,…,ur)T，B＝(v1,…,vr)T时，存在||X||T＝||X||*-Tr(AXBT)；
构建的目标函数：
Figure BDA0003934504310000103
根据凸性条件ITI-λ2HTH≥0、
Figure BDA0003934504310000104
0≤γ≤1改写目标函数为凸性目标函数：
Figure BDA0003934504310000105
其中，X初始化为与矩阵Y规格一致的零矩阵；L表示矩阵X的近似矩阵，X＝Rm×n；||X||F代表X的F范数，表示为
Figure BDA0003934504310000106
||X||1代表X的1范数，表示为
Figure BDA0003934504310000107
Tr(X)表示矩阵X的迹，第一辅助矩阵A与第二辅助矩阵B满足A∈Rm×n、BT∈Rm×n、AAT＝I、BBT＝I，I表示单位矩阵；
将矩阵X的奇异值分解定义为X＝UΣVT：
Figure BDA0003934504310000108
Ak＝(Uk(:,[1:r]))T，Bk＝(Vk(:,[1:r]))T；
其中，U＝(u1,…,um)∈Rm×m，Σ∈Rm×n，V＝(v1,…,vm)∈Rn×n；当第一辅助矩阵A与第二辅助矩阵B满足A＝(u1,…,ur)T，B＝(v1,…,vr)T时，存在||X||T＝||X||*-Tr(AXBT)，σi表示时频域系数矩阵X的第i个奇异值；r为截断阶数，表示将被保留的最大的r个非零奇异值；H表示将所述广义极小极大凹罚函数参数化的参数矩阵，满足HTH＝nI；λ1为第一正则化参数、λ2为第二正则化参数，γ为常量参数。
S22：利用交替方向乘子法将所述凸性目标函数拆分成关于广义极小极大凹罚函数的第一子目标函数、关于截断核范数的第二子目标函数与辅助迭代函数：
引入变量Z，Z满足Z＝X，带入所述凸性目标函数中得：
Figure BDA0003934504310000111
所述关于广义极小极大凹罚函数的第一子目标函数：
Figure BDA0003934504310000112
所述关于截断核范数的第二子目标函数：
Figure BDA0003934504310000113
辅助迭代函数：R＝R-(X-Z)；
其中，μ表示迭代步长，辅助矩阵R初始化为与所述时频域系数矩阵Y规格一致的零矩阵。
S3：利用前向后向分裂算法，求解所述第一子目标函数，得到关于Xk+1的迭代公式组；
S31：
Figure BDA0003934504310000114
分别对X和L求偏导得到：
Figure BDA0003934504310000115
Figure BDA0003934504310000116
其中，符号函数
Figure BDA0003934504310000117
S32：将对X和L求偏导得到的公式带入前向后向分裂算法中得：
Qk＝Xk-ξ((1+μ)Xk-Y-μ(Zk+Rk)-γ(Xk+Lk))，
Xk+1＝soft(Qk,λ2ξ)，
Tk＝Lk+ξγ(Xk+1-Lk)，
Lk+1＝soft(Tk,λ2ξ)，
其中，第一矩阵Qk为Xk梯度下降得到的矩阵，第二矩阵Tk为近似矩阵Lk梯度下降得到的矩阵,ξ表示步长，k表示迭代次数，soft代表软阈值算子，soft(x,C)＝x·max{0,1-C/|x|}。
S4：根据所述关于Xk+1的迭代公式组，利用奇异值阈值算法求解所述第二目标子函数，得到关于Zk+1的迭代公式；
S41：根据所述关于Xk+1的迭代公式组，利用奇异值阈值算法求解所述第二目标子函数，得到关于Zk+1的迭代公式：
将
Figure BDA0003934504310000121
转化为：
Figure BDA0003934504310000122
利用奇异值阈值算法求解，得Z的迭代公式：
Figure BDA0003934504310000123
Figure BDA0003934504310000124
其中，Ak为第一辅助矩阵，Bk为第二辅助矩阵，Pk为左奇异矩阵，
Figure BDA0003934504310000125
为右奇异矩阵，∑k为非负实对角矩阵。
S42：对奇异值阈值算法进行加权求解所述第二目标子函数，更新关于Zk+1的迭代公式：
构建权值向量W与所述奇异值阈值算法中的阈值
Figure BDA0003934504310000126
相乘对阈值进行加权；
其中，权值向量W中的元素定义为：
Figure BDA0003934504310000127
σi＝∑(i,i)，σi表示奇异值，ε为常数；
将权值向量W引入所述Z的迭代公式中，更新的Z的迭代公式表示为：
Figure BDA0003934504310000131
S5：利用关于Xk+1的迭代公式组与关于Zk+1的迭代公式更新所述辅助迭代函数，初始化参数，将所述时频域系数矩阵Y输入所述辅助迭代函数中，迭代预设次数Nit次获取对应的时频域稀疏低秩矩阵XNit；
S51：利用关于Xk+1的迭代公式组与关于Zk+1的迭代公式更新所述辅助迭代函数：Rk+1＝Rk-(Xk+1-Zk+1)；
S52：设置迭代次数k＝0,1,2,3,…,Nit，输入时频域系数矩阵Y，初始化迭代步长μ、步长ξ、第一正则化参数λ1与第二正则化参数λ2、截断阶数r、常量参数γ；
S53：初始化迭代次数k＝0，初始化矩阵Xk、变量Zk、辅助矩阵Rk与近似矩阵Lk为与所述时频域系数矩阵Y规格一致的零矩阵：X0、Z0、R0与L0；
对X0进行梯度下降得到的第一矩阵Q0，利用软阈值函数，根据Q0迭代求出X1；对L0进行梯度下降得到第二矩阵T0，利用软阈值函数，根据T0迭代求解出L1；奇异值阈值算法求解得到Z1；根据R0、X1与Z1求解出R1；根据所述X1、L1、Z1与R1更新得到Q1；
S54：对于k＝1,2,3,…,Nit；
根据Xk、Zk、Rk与Lk更新得到Qk，利用软阈值函数，根据Qk求解Xk+1；
根据Xk+1更新Tk，利用软阈值函数，根据Tk更新Lk+1；
利用奇异值算法，根据Rk与Xk+1更新Zk+1；
根据Rk、Xk+1与Zk+1更新得到Rk+1；
利用Xk+1、Zk+1、Rk+1与Lk+1更新得到Qk+1；
S55：令k＝k+1，重复S54，直至k+1＝Nit，停止迭代，输出相对应的时频域稀疏低秩矩阵XNit；
S6：对所述时频域稀疏低秩矩阵XNit进行短时傅里叶逆变换得到轴承故障瞬态脉冲的重构时域信号：x*＝S-(XNit)；
其中，x*表示轴承故障瞬态脉冲的重构时域信号，S-表示短时傅里叶逆变换；XNit表示根据Rk+1的迭代公式求得的当k达到预设次数时的Xk；
S7：对所述轴承故障瞬态脉冲的重构时域信号进行平方包络谱分析，得到轴承故障特征频率。
基于上述实施例，本发明实施例还提供了一种轴承故障检测装置，包括：轴承故障信号预处理模块，用于获取测试轴承的振动加速度时域信号，利用短时傅里叶变换生成时频域系数矩阵；稀疏低秩优化处理模块，用于利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型的目标函数，并利用交替方向乘子法求解获取目标函数迭代公式，将所述时频域系数矩阵输入所述目标函数迭代公式，迭代获取对应的时频域稀疏低秩矩阵；时域信号重构模块，用于对所述时频域稀疏低秩矩阵进行短时傅里叶逆变换得到轴承故障瞬态脉冲的重构时域信号；包络谱分析模块，用于对所述轴承故障瞬态脉冲的重构时域信号进行平方包络谱分析，得到轴承故障特征频率。
基于上述实施例，本发明实施例还提供了一种轴承故障检测设备，包括：采集装置，用于利用加速度传感器采集测试轴承的加速度时域信号；轴承故障检测装置，与所述采集装置通讯连接，用于执行计算机程序实现如上述所述的轴承故障检测方法；显示装置，与所述轴承故障检测装置通讯连接，用于显示所述重构时域信号的平方包络谱图。
具体地，参照图2所示，采集装置包括：驱动电机；转动轴，一端通过联轴器连接所述驱动电机；测试轴承安装座，套设于所述转动轴上，用于固定测试轴承；加速度传感器，设置于所述转动轴的另一端，用于连接所述测试轴承并测量其振动加速度；负载装置，用于添加可调节的负载至所述转动轴。采集装置还可在转动轴上设置健康轴承，支撑转动轴，平衡施加负载后转动轴的偏心问题。
基于上述实施例，在本实施例中，对本发明实施例所提供的轴承故障检测方法进行验证；本发明实施例所采用的轴承故障的振动加速度时域信号是利用上述采集装置采集的，将加速度传感器放置在测试轴承座上采集振动加速度的时域信号；轴承故障模拟实验台的实验对象是型号为SKF6205-2RS的深沟球轴承，通过线切割加工的割痕来模拟轴承外圈故障，采样频率为10kHz，该轴承外圈故障特征频率的理论值为62.5Hz；本实施例检测的处理器为i5-11400，整个程序的耗时约为2.394s，检测结果准确，论证如下：
具体地，轴承故障检测的具体步骤包括：
将故障轴承振动加速度的时域信号y载入，设置采样频率为10kHz，采样点数为4000，将这4000个点所组成的信号段进行短时傅里叶变换得到其对应的时频域系数矩阵Y：Y＝S(y)；
其中，Y表示时频域系数矩阵；y表示振动加速度时域信号；S表示短时傅里叶变换，其设置的窗长为64，重合点数为窗长的50％，窗函数为余弦窗，快速傅里叶变换的点数为512；
如图3所示，为故障轴承振动加速度的时域信号图，如图4所示，为故障轴承振动加速度时域信号经过短时傅里叶变换得到的时频图；根据图3与图4可得，从故障轴承振动加速度的时域信号难以看出周期性的轴承故障瞬态脉冲，其时域信号没有明显的稀疏性，且其时频域也没有呈现稀疏性。
对所述时频域系数矩阵进行稀疏低秩优化处理，获取对应的时频域稀疏低秩矩阵；利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型；并利用交替方向乘子法将求解所述稀疏低秩优化模型的问题分解为利用前向后向分裂算法解决广义极小极大凹罚函数优化子问题与利用奇异值阈值算法解决截断核范数优化子问题，以获得稀疏低秩矩阵Xopt，具体步骤包括：
输入时频域系数矩阵Y，初始化变量X，V,Z,R为与时频域系数矩阵Y规格一致的零矩阵；
输入并初始化正则化参数λ1＝0.4,λ2＝0.1，步长μ＝1，ξ＝0.6，截断阶数r＝10，常量参数γ＝0.65,ε＝0.05,迭代次数Nit＝30；
对于k＝0,1,2,3,…,Nit，计算Ak和Bk：
Figure BDA0003934504310000161
Ak＝(Uk(:,[1:r]))T，Bk＝(Vk(:,[1:r]))T；
对于k＝0,1,2,3,…,Nit，计算Xk+1：
Qk＝Xk-ξ((1+μ)Xk-Y-μ(Zk+Rk)-γ(Xk+Lk))；
Xk+1＝soft(Qk,λ2ξ)；
对于k＝0,1,2,3,…,Nit，计算Vk+1：
Tk＝Lk+ξγ(Xk+1-Lk)；
Lk+1＝soft(Tk,λ2ξ)；
对于k＝0,1,2,3,…,Nit，计算Zk+1：
Figure BDA0003934504310000162
设置权值向量Wk，Wk中的元素为wi:i＝0,1,2,3,…,n，n＝min(size(Y))，
Figure BDA0003934504310000163
其中σi＝∑(i,i)；ε为常数，避免零分母情况，设为0.05；
将所述权值向量Wk与奇异值阈值算法中的阈值相乘：
Figure BDA0003934504310000164
对于k＝0,1,2,3,…,Nit，计算Rk+1：
Rk+1＝Rk-(Xk+1-Zk+1)
得到输出结果XNit，；
通过对时频域稀疏低秩矩阵XNit进行短时傅里叶逆变换，有x*＝S-(XNit)，从而得到轴承故障瞬态脉冲的重构时域信号x*。
参照图5所示，为故障轴承振动加速度时域信号的时频图；参照图6所示，为轴承故障瞬态脉冲的重构时域信号的时频图；通过对比图5与图6可知，通过稀疏低秩优化处理后成狗的轴承故障瞬态脉冲信号的时频图呈现了明显的稀疏性。参照图7所示，为故障轴承振动加速度的时频域系数矩阵与轴承故障瞬态脉冲的重构时域信号的时频系数矩阵的奇异值分布对比图，可以得知，相较于故障轴承振动加速度的时频域系数矩阵有近65个非零奇异值，轴承故障瞬态脉冲的重构时域信号的时频系数矩阵，仅有16个非零奇异值，呈现明显的低秩性，并且大的奇异值得到了较好的保留，提升了轴承故障检测的准确性。参照图8所示，为轴承故障瞬态脉冲的重构时域信号图，对比图3可知，本发明实施例所求得的轴承故障瞬态脉冲重构时域信号的特征，具有明显的稀疏性与周期性，处理效果好。
对所得轴承故障瞬态脉冲信号x*做平方包络谱分析，进一步可从轴承故障瞬态脉冲信号x*的平方包络谱图像中得到故障特征频率及其倍频等信息。
参照图9所示，图9的(a)为轴承故障瞬态脉冲信号的平方包络谱图，图9的(b)为轴承故障瞬态脉冲信号的平方包络谱局部图；根据图9可以明显观察到轴承故障特征频率及其倍频，可以得知本发明实施例所处理的轴承故障振动信号的故障特征频率为62.5Hz，与故障特征频率的理论值一致，证明了本发明所提供的方法的有效性与准确性。
本发明所提供的轴承故障检测方法，将含有噪声的一维轴承故障振动加速度的时域信号进行短时傅里叶变换得到二维时频域系数矩阵；对二维时频域系数矩阵进行改进的稀疏低秩优化处理生成稀疏低秩矩阵；对稀疏低秩矩阵进行短时傅里叶逆变换得到轴承故障瞬态脉冲的重构时域信号。在进行改进的稀疏低秩优化处理时，对于时频域系数矩阵的稀疏性，设计广义极小极大凹罚函数逼近L0范数，提升稀疏低秩优化模型的保幅值能力与去噪效果；对于时频域系数矩阵的低秩性，利用截断核范数逼近时频域系数矩阵的秩，并对传统奇异值阈值算法进行改进，设计了与奇异值大小呈负相关的权值向量，利用权值向量对奇异值阈值进行加权处理，提升了奇异值保留效果，从而提升了轴承故障检测的准确性。
显然，上述实施例仅仅是为清楚地说明所作的举例，并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明创造的保护范围之中。

Claims (10)

1.一种轴承故障检测方法，其特征在于，包括：
S1：获取测试轴承的振动加速度时域信号，利用短时傅里叶变换生成时频域系数矩阵Y；
S2：利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型的目标函数，并根据其凸性条件转换为凸性目标函数；利用交替方向乘子法将所述凸性目标函数拆分成关于广义极小极大凹罚函数的第一子目标函数、关于截断核范数的第二子目标函数与辅助迭代函数；
S3：利用前向后向分裂算法，求解所述第一子目标函数，得到关于Xk+1的迭代公式组；
S4：根据所述关于Xk+1的迭代公式组，利用奇异值阈值算法求解所述第二子目标函数，得到关于Zk+1的迭代公式；
S5：利用关于Xk+1的迭代公式组与关于Zk+1的迭代公式更新所述辅助迭代函数，初始化参数，将所述时频域系数矩阵Y输入所述辅助迭代函数中，迭代预设次数Nit次获取对应的时频域稀疏低秩矩阵XNit；
S6：对所述时频域稀疏低秩矩阵XNit进行短时傅里叶逆变换得到轴承故障瞬态脉冲的重构时域信号；
S7：对所述轴承故障瞬态脉冲的重构时域信号进行平方包络谱分析，得到轴承故障特征频率。
2.根据权利要求1所述的轴承故障检测方法，其特征在于，所述利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型的目标函数，并根据其凸性条件转换为凸性目标函数包括：
所述稀疏低秩优化模型的目标函数：
Figure QLYQS_1
所述凸性条件为ITI-λ2HTH≥0，
Figure QLYQS_2
0≤γ≤1；
所述凸性目标函数：
Figure QLYQS_3
其中，X初始化为与矩阵Y规格一致的零矩阵；L表示矩阵X的近似矩阵，X＝Rm×n；||X||F代表X的F范数，表示为
Figure QLYQS_4
||X||1代表X的1范数，表示为
Figure QLYQS_5
Tr(X)表示矩阵X的迹，第一辅助矩阵A与第二辅助矩阵B满足A∈Rm×n、BT∈Rm×n、AAT＝I、BBT＝I，I表示单位矩阵；
其中，将矩阵X的奇异值分解定义为X＝UΣVT：
Figure QLYQS_6
Ak＝(Uk(:,[1:r]))T，Bk＝(Vk(:,[1:r]))T；
U＝(u1,…,um)∈Rm×m，Σ∈Rm×n，V＝(v1,…,vm)∈Rn×n；当第一辅助矩阵A与第二辅助矩阵B满足A＝(u1,…,ur)T，B＝(v1,…,vr)T时，存在||X||T＝||X||*-Tr(AXBT)，σi表示时频域系数矩阵X的第i个奇异值；r为截断阶数，表示将被保留的最大的r个非零奇异值；H表示将所述广义极小极大凹罚函数参数化的参数矩阵，满足HTH＝nI；λ1为第一正则化参数、λ2为第二正则化参数，γ为常量参数。
3.根据权利要求2所述的轴承故障检测方法，其特征在于，所述利用交替方向乘子法将所述凸性目标函数拆分成关于广义极小极大凹罚函数的第一子目标函数、关于截断核范数的第二子目标函数与辅助迭代函数包括：
引入变量Z，Z满足Z＝X，带入所述凸性目标函数中得：
Figure QLYQS_7
所述关于广义极小极大凹罚函数的第一子目标函数：
Figure QLYQS_8
所述关于截断核范数的第二子目标函数：
Figure QLYQS_9
辅助迭代函数：R＝R-(X-Z)；
其中，μ表示迭代步长，辅助矩阵R初始化为与所述时频域系数矩阵Y规格一致的零矩阵。
4.根据权利要求3所述的轴承故障检测方法，其特征在于，所述利用前向后向分裂算法，求解所述第一子目标函数，得到关于Xk+1的迭代公式组包括：
Figure QLYQS_10
分别对X和L求偏导得到：
Figure QLYQS_11
Figure QLYQS_12
其中，符号函数
Figure QLYQS_13
将对X和L求偏导得到的公式带入前向后向分裂算法中得关于Xk+1的迭代公式组：
Qk＝Xk-ξ((1+μ)Xk-Y-μ(Zk+Rk)-γ(Xk+Lk))；
Xk+1＝soft(Qk,λ2ξ)；
Tk＝Lk+ξγ(Xk+1-Lk)；
Lk+1＝soft(Tk,λ2ξ)；
其中，辅助矩阵R是初始化为与所述时频域系数矩阵Y规格一致的零矩阵，第一矩阵Qk为Xk梯度下降得到的矩阵，第二矩阵Tk为近似矩阵Lk梯度下降得到的矩阵,ξ表示步长，k表示迭代次数，soft代表软阈值算子，soft(x,C)＝x·max{0,1-C/|x|}。
5.根据权利要求4所述的轴承故障检测方法，其特征在于，所述根据所述关于Xk+1的迭代公式组，利用奇异值阈值算法求解所述第二目标子函数，得到关于Zk+1的迭代公式包括：
将
Figure QLYQS_14
转化为：
Figure QLYQS_15
利用奇异值阈值算法求解，得关于Zk+1的迭代公式：
Figure QLYQS_16
Figure QLYQS_17
其中，Ak为第一辅助矩阵，Bk为第二辅助矩阵，Pk为左奇异矩阵，
Figure QLYQS_18
为右奇异矩阵，∑k为非负实对角矩阵。
6.根据权利要求5所述的轴承故障检测方法，其特征在于，所述利用奇异值阈值算法求解所述第二目标子函数，得到关于Zk+1的迭代公式还包括：
构建权值向量W与所述奇异值阈值算法中的阈值
Figure QLYQS_19
相乘对阈值进行加权；
其中，权值向量W中的元素定义为：
Figure QLYQS_20
σi＝∑(i,i)，σi表示奇异值，ε为常量常数；
将权值向量W引入所述Z的迭代公式中，更新的关于Zk+1的迭代公式表示为：
Figure QLYQS_21
7.根据权利要求6所述的轴承故障检测方法，其特征在于，所述利用关于Xk+1的迭代公式组与关于Zk+1的迭代公式更新所述辅助迭代函数，初始化参数，将所述时频域系数矩阵Y输入所述辅助迭代函数中，迭代预设次数Nit次获取对应的时频域稀疏低秩矩阵XNit：
S51：利用关于Xk+1的迭代公式组与关于Zk+1的迭代公式更新所述辅助迭代函数：Rk+1＝Rk-(Xk+1-Zk+1)；
S52：设置迭代次数k＝0,1,2,3,…,Nit，初始化迭代步长μ、步长ξ、第一正则化参数λ1与第二正则化参数λ2、截断阶数r、常量参数γ；
S53：初始化迭代次数k＝0，输入时频域系数矩阵Y，初始化矩阵Xk、变量Zk、辅助矩阵Rk与近似矩阵Lk为与所述时频域系数矩阵Y规格一致的零矩阵：X0、Z0、R0与L0；
对X0进行梯度下降得到的第一矩阵Q0，利用软阈值函数，根据Q0迭代求出X1；对L0进行梯度下降得到第二矩阵T0，利用软阈值函数，根据T0迭代求解出L1；奇异值阈值算法求解得到Z1；根据R0、X1与Z1求解出R1；根据所述X1、L1、Z1与R1更新得到Q1；
S54：对于k＝1,2,3,…,Nit；
根据Xk、Zk、Rk与Lk更新得到Qk，利用软阈值函数，根据Qk求解Xk+1；
根据Xk+1更新Tk，利用软阈值函数，根据Tk更新Lk+1；
利用奇异值算法，根据Rk与Xk+1更新Zk+1；
根据Rk、Xk+1与Zk+1更新得到Rk+1；
利用Xk+1、Zk+1、Rk+1与Lk+1更新得到Qk+1；
S55：令k＝k+1，重复S54，直至k+1＝Nit，停止迭代，输出XNit。
8.一种轴承故障检测装置，其特征在于，包括：
轴承故障信号预处理模块，用于获取测试轴承的振动加速度时域信号，利用短时傅里叶变换生成时频域系数矩阵；
稀疏低秩优化处理模块，用于利用广义极小极大凹罚函数与截断核范数构建稀疏低秩优化模型的目标函数，并利用交替方向乘子法求解获取目标函数迭代公式，将所述时频域系数矩阵输入所述目标函数迭代公式，迭代获取对应的时频域稀疏低秩矩阵；
时域信号重构模块，用于对所述时频域稀疏低秩矩阵进行短时傅里叶逆变换得到轴承故障瞬态脉冲的重构时域信号；
包络谱分析模块，用于对所述轴承故障瞬态脉冲的重构时域信号进行平方包络谱分析，得到轴承故障特征频率。
9.一种轴承故障检测设备，其特征在于，包括：
采集装置，用于利用加速度传感器采集测试轴承的加速度时域信号；
轴承故障检测装置，与所述采集装置通讯连接，用于执行计算机程序实现如权利要求1至7任一项所述的轴承故障检测方法；
显示装置，与所述轴承故障检测装置通讯连接，用于显示所述重构时域信号的平方包络谱图。
10.根据权利要求9所述的轴承故障检测设备，其特征在于，所述采集装置包括：
驱动电机；
转动轴，一端通过联轴器连接所述驱动电机；
测试轴承安装座，套设于所述转动轴上，用于固定测试轴承；
加速度传感器，设置于所述转动轴的另一端，用于连接所述测试轴承并测量其振动加速度；
负载装置，用于添加可调节的负载至所述转动轴。
CN202211398933.1A 2022-11-09 2022-11-09 一种轴承故障检测方法、装置及设备 Active CN115683629B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211398933.1A CN115683629B (zh) 2022-11-09 2022-11-09 一种轴承故障检测方法、装置及设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211398933.1A CN115683629B (zh) 2022-11-09 2022-11-09 一种轴承故障检测方法、装置及设备

Publications (2)

Publication Number Publication Date
CN115683629A CN115683629A (zh) 2023-02-03
CN115683629B true CN115683629B (zh) 2023-06-27

Family

ID=85050580

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211398933.1A Active CN115683629B (zh) 2022-11-09 2022-11-09 一种轴承故障检测方法、装置及设备

Country Status (1)

Country Link
CN (1) CN115683629B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110991419A (zh) * 2019-12-23 2020-04-10 长安大学 基于稀疏低秩协同优化框架的齿轮箱局部故障诊断方法
CN112101082A (zh) * 2020-11-16 2020-12-18 华南理工大学 一种基于改进低秩稀疏分解的旋转机械故障诊断方法
CN112284727A (zh) * 2020-09-30 2021-01-29 华南理工大学 一种基于卷积极大极小凹罚算法的旋转机械故障诊断方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9317780B2 (en) * 2013-10-17 2016-04-19 Xerox Corporation Detecting multi-object anomalies utilizing a low rank sparsity model

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110991419A (zh) * 2019-12-23 2020-04-10 长安大学 基于稀疏低秩协同优化框架的齿轮箱局部故障诊断方法
CN112284727A (zh) * 2020-09-30 2021-01-29 华南理工大学 一种基于卷积极大极小凹罚算法的旋转机械故障诊断方法
CN112101082A (zh) * 2020-11-16 2020-12-18 华南理工大学 一种基于改进低秩稀疏分解的旋转机械故障诊断方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于低秩稀疏分解算法的航空锥齿轮故障诊断;陈礼顺;张晗;陈雪峰;程礼;曾林;;振动与冲击(12);108-117 *

Also Published As

Publication number Publication date
CN115683629A (zh) 2023-02-03

Similar Documents

Publication Publication Date Title
CN105827250B (zh) 一种基于自适应字典学习的电能质量数据压缩重构方法
CN110361778B (zh) 一种基于生成对抗网络的地震数据重建方法
Xiao et al. Multivariate global sensitivity analysis for dynamic models based on wavelet analysis
Gribonval et al. Blind calibration for compressed sensing by convex optimization
CN106644464A (zh) 一种基于载荷谱分析的轧机传动系统关键零部件的疲劳寿命预警方法
CN104748961A (zh) 基于svd分解降噪和相关性eemd熵特征的齿轮故障诊断方法
Hou et al. Fault diagnosis for rolling bearings under unknown time-varying speed conditions with sparse representation
CN111582137B (zh) 一种滚动轴承信号重构方法及系统
Dong et al. A repeated single-channel mechanical signal blind separation method based on morphological filtering and singular value decomposition
CN111324861A (zh) 一种基于矩阵分解的深度学习磁共振波谱重建方法
Medina et al. Deep learning-based gear pitting severity assessment using acoustic emission, vibration and currents signals
CN115683629B (zh) 一种轴承故障检测方法、装置及设备
CN113935246A (zh) 一种信号鲁棒稀疏时频分析方法、终端设备及存储介质
Ma et al. Sparse low-rank matrix estimation with nonconvex enhancement for fault diagnosis of rolling bearings
Wang et al. Wire rope damage detection signal processing using K-singular value decomposition and optimized double-tree complex wavelet transform
CN114863209B (zh) 类别比例引导的无监督领域适应建模方法、系统、设备及介质
Chen et al. Application of EMD-AR and MTS for hydraulic pump fault diagnosis
Cuomo et al. An inverse Bayesian scheme for the denoising of ECG signals
Prichard et al. Is the AE index the result of nonlinear dynamics?
Alcaraz et al. Signal noise filtering using wavelet coefficient temporal correlation techniques
Chen et al. Bearing fault detection based on SVD and EMD
CN112329626A (zh) 调制与深度学习融合的设备故障诊断方法、系统及介质
Guo et al. Application of EEMD singular value energy spectrum in gear fault identification
Paul et al. Application of discrete domain wavelet filter for signal denoising
CN113239608B (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