CN104677632B - 利用粒子滤波与谱峭度的滚动轴承故障诊断方法 - Google Patents
利用粒子滤波与谱峭度的滚动轴承故障诊断方法 Download PDFInfo
- Publication number
- CN104677632B CN104677632B CN201510031512.9A CN201510031512A CN104677632B CN 104677632 B CN104677632 B CN 104677632B CN 201510031512 A CN201510031512 A CN 201510031512A CN 104677632 B CN104677632 B CN 104677632B
- Authority
- CN
- China
- Prior art keywords
- signal
- kurtosis
- filtering
- filter
- fault
- 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
Links
Landscapes
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
Abstract
本发明公开了一种利用粒子滤波与谱峭度的滚动轴承故障诊断方法，涉及粒子滤波降噪处理以及谱峭度计算。该方法提出了在快速谱峭度的基础上，利用粒子滤波降噪处理提高信噪比，解决了快速谱峭度在低信噪比情况下可行性差的问题。首先建立振动信号的状态方程；然后提取背景噪声，将其与状态方程之和作为观测方程；最后联立状态方程与观测方程建立状态空间模型；采用粒子滤波算法对信号重新估计，得到新的序列即是降噪之后的信号；最后用快速谱峭度的方法，获取最佳分析频带，得出故障频率。本发明的滚动轴承故障诊断方法，降低故障信号中的噪声干扰，提高信噪比，实现了滚动轴承早期微弱故障的诊断。
Description
技术领域
本发明涉及一种利用粒子滤波与谱峭度的滚动轴承故障诊断方法，涉及粒子滤波降噪处理以及谱峭度计算。
背景技术
滚动轴承是机械设备尤其是旋转机械中重要的零部件之一，但是其寿命的随机性较大，容易损坏，目前还无法准确预测其寿命的长短。近年来，科技和工业生产都在快速的发展，机械设备逐渐向高速化、大型化和自动化发展，这在提高生产力的同时对设备的安全维护也提出了更高的要求，在某一环节出现一个微小的问题都将造成不可挽回的巨额损失。鉴于滚动轴承在机械设备中的地位，可以得出滚动轴承能否正常运行工作关系着机械设备乃至生产线的正常运转，因而掌握滚动轴承运行的工作状态以及故障的形成和发展，是目前机械故障诊断领域中所研究的重要课题之一。
根据故障诊断技术机理的不同，诊断滚动轴承故障的常用方法有振动分析法、噪声分析法、油样分析法、温度分析法、油膜电阻法、声发射诊断等。这些分析方法中，振动分析法是运用最广泛的一种方法。振动信号包含丰富的设备异常或故障的信息，通过对振动信号的实时在线监测和分析，可以判断机械系统及其部件是否正常运行，并以此确定机械故障产生的原因，判断故障产生的具体部位和损伤程度。较常用的振动诊断方法有倒频谱分析、特征参数分析法、冲击脉冲法、包络分析法、小波分析等。虽然诊断方法很多，但是对早期的微弱故障诊断仍然是一个难题。因为在故障早期，所需频段信号能量非常微弱、特征信号成分极易被噪声和其他频段信号俺没，未经任何处理的原始频段微弱信号中很难提取出有用的特征信号成分，如何提取被淹没在噪声无用频段内的微弱故障特征信号是目前国内外研究焦点。
发明内容
振动信号是典型的非平稳、非线性的时序信号，粒子滤波在解决非线性滤波问题时，随机量不会受到高斯分布的制约，比高斯模型具有更广泛的分布，对非线性特性有很强的建模能力。本发明提供了一种利用粒子滤波与谱峭度的滚动轴承故障诊断方法，解决了普通的快速谱峭度方法难以实现在低信噪比的情况下诊断出轴承故障的问题。
一种利用粒子滤波与谱峭度的滚动轴承故障诊断方法，具体过程如下：
粒子滤波是以状态空间方程为主对信号滤波，所以首先要建立轴承振动信号的状态空间方程。对于一个非线性、非高斯系统，状态空间模型的一般表达式为：
其中，vk和nk分别是过程噪声和观测噪声，并且都是相互独立的，协方差分别是Qk和Rk的零均值噪声序列，f和h为状态转移方程和观测方程。
众所周知，轴承振动信号是典型的非线性、非高斯时间序列。所以用AR时序模型作为状态方程，模型的残差方差即为过程噪声。提取原始信号的背景噪声作为观测噪声，将其与状态方程加载一起作为观测方程。综上所述，振动信号的状态空间表达式可以写为：
式中，p为AR模型阶数，at(t)是模型参数序列，为模型残差方差，即过程噪声，为信号背景噪声，具体获得方法为：用小波分解函数对采集的振动信号进行多尺度分解,得到小波变换系数。再设定合适的阈值,将大于阈值的小波系数置为零,小于阈值的小波系数保留下来。对经过处理的小波系数进行小波重构,就得到所需的噪声信号，即观测噪声。et和nt都是服从N～(0,1)的标准正太分布。模型参数估计方法为Yule-Walker法，模型定界准则是AIC准则。
粒子滤波的本质是利用蒙特卡罗积分法求解贝叶斯积分运算，它的关键是核心粒子的权值大小要能不断被修正和调整，最后获得精确的估计值，下文将标准粒子滤波的具体流程。
首先用一组随机粒子和这些粒子的权值大小来近似表达后验概率密度函数p(xk|z1:k)：
式中，δ(·)是狄拉克函数，z1:k是观测量。如果这些粒子是通过真实概率密度函数p(xk|z1:k)获得，那么式(3)可写为：
然而，在很多情况下获得真实概率密度函数粒子是非常困难的，所以就利用重要性概率函数q(xk|z1:k)获取粒子。假设p(xk|z1:k)与一个分析函数π(xk|z1:k)是线性关系，则是(3)中的权值大小可以通过下式计算：
假设重要性函数可以分解为：
将式(5)和时(6)代入到式(4)中更新粒子权值大小
当满足条件:的时候，式(7)可以写为：
又因为先验概率密度函数通常选择为：所以式(8)可以写为：
粒子滤波结束后结合快速谱峭度法提取特征故障频率。谱峭度的定义式为：
式中，为信号的n阶谱矩，其中，E<·>——均值运算，E|·|取模运算。X(f,t)是信号x(t)在频率f处的复包络。
快速谱峭度算法采用的是在全频带范围内1/3和二叉树结构设计的带通滤波器分析信号。算法如下：
(1)构造一个截止频率fc＝1/8+ε的低通滤波器h(n)，其中fs＝1,ε＞0。低通滤波器的截止频率增大ε是为了消除过度带的影响。并基于h(n)构造准解析低通滤波器h0(n)和准解析高通滤波器h1(n)，其分析频带分别为[0,1/4]、[1/4,1/2]，即：
(2)分别以h0(n)和h1(n)低通和高通滤波，对滤波结果做2倍降采样，如图1所示。如此迭代进行，可以获得如图2所示的滤波器树及相应的滤波结果。其中表示第k层第i个滤波器产生的滤波结果，其中i取值为0到2k-1；同时也是信号在中心频带fi＝(i+2-1)2k-1和带宽(Δf)k＝2k-1的复包络，2倍降采样的目的是保证滤波器每一层中的数据长度与原始数据相同。
(3)对于每一个滤波器结果，按照式(10)计算峭度：
(4)将所有的谱峭度汇总，从而得到“快速谱峭度图”。
根据峭度图选择峭度最大值频带区间解调分析，计算功率谱，对功率谱做FFT变换得到频谱图，根据频谱信息做出故障类型的诊断。
本发明提供了一种利用粒子滤波与谱峭度的滚动轴承故障诊断方法，降低故障信号中的噪声干扰，提高信噪比，实现了滚动轴承早期微弱故障的诊断。
附图说明
图1为低通/高通滤波器分解图。
图2为二叉树滤波器组及相应的滤波结果。
图3内环故障时域振动信号。
图4外环故障时域振动信号。
图5内环振动信号的AIC曲线图。
图6外环振动信号的AIC曲线图。
图7滤波之后的内环时域信号。
图8滤波之后的外环时域信号。
图9为内环信号降噪后频谱图。
图10为外环信号降噪后频谱图。
图11内环信号滤波之前的频谱图。
图12外环信号滤波之前的频谱图。
具体实施方式
以下结合技术方案和附图详细叙述本发明的实施例。
实施例的数据来自美国西储大学，数据为电机驱动端轴承故障数据,轴承型号为6205-2RS JEM SKF，深沟球轴承。轴承参数如下表所示：
表1 轴承参数表
| 内圈直径 | 外圈直径 | 厚度 | 滚动体直径Dd | 节经dm | 接触角α | 滚珠数Z |
| 25mm | 52mm | 15mm | 7.94mm | 139.04mm | 60° | 9 |
实验电机转速1797rpm/m，故障尺寸为0.1778mm，内外环的故障特征频率计算公式如下：内环
外环
观察图3和图4外所示的内环、外环时域振动信号，无法看出冲击信号，所以不能判断是否存在故障。要得到故障特征频率只有作进一步的分析处理。
图5和图6分别是内环和外环的AIC曲线图。AIC准则认为：AIC值越小，模型精度越高。但是这样模型的阶数就会变大，计算量就会上升。在上升到一定阶段时，模型的精确度度提高非常缓慢。综合考虑两者因素，我们选定当AIC值下降缓慢时候的阶数为模型阶数。内环的AIC曲线在阶数为20时，值已经趋于平缓，下降幅度很小；而外环在第8阶时下降幅度最小。所以模型阶数分别定位20阶、8阶。
其次采用小波降噪方法得到信号的背景噪声，即为观察噪声，观察噪声与AR模型之和即为观察方程；最后联立AR模型与观察方程便是粒子滤波所需的状态空间方程。利用粒子滤波方法对信号重新估计滤波，得到滤波后的时域信号如图7、图8所示，对比观察滤波前后的信号，可以发现滤波之后信号幅值下降的明显。结合谱峭度方法计算峭度值，根据峭度图选取峭度值最大子带为最优解调分析频带。提取特征频率如图9、图10所示。
在图9、图10频谱图中，很明显的观察到故障特征频率以及它们的倍频信号，可以说明本方法能成功地实现轴承故障诊断。为了进一步体现所提出方法的优越性，对比分析了粒子滤波之前的处理结果，如图11、图12所示。
图11是内环滤波之前的频域图，图中只能观察到轴频信号(因为实验转速为1797rpm/m，所以轴频是29.95Hz)及其倍频，没有故障频率。说明滤波之前谱峭度无法得到内环的故障信息。图12是外环振动信号滤波之前的频域图，频谱图中虽然可以看到故障频率107.7Hz，但是存在着边频带谱线，干扰信号太多，可视效果很差。再反观滤波之后频谱如图10所示，边频带频率被成功滤除，没有干扰信号，结果清晰一目了然。
上述对比分析再一次证明了所提出的方法的优越性。
Claims (1)
1.利用粒子滤波与谱峭度的滚动轴承故障诊断方法，其特征在于如下步骤：
(一)滤波降噪，过程如下：
(1)初始化；由先验概率p(x0)产生粒子群所有粒子权值为1/N；
(2)更新；在k时刻，更新粒子权值，
<mrow>
<msubsup>
<mi>w</mi>
<mi>k</mi>
<mi>i</mi>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>w</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mi>i</mi>
</msubsup>
<mi>p</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>z</mi>
<mi>k</mi>
</msub>
<mo>|</mo>
<msubsup>
<mi>x</mi>
<mi>k</mi>
<mi>i</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msubsup>
<mi>w</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mi>i</mi>
</msubsup>
<mi>p</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>z</mi>
<mi>k</mi>
</msub>
<mo>-</mo>
<mi>h</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>x</mi>
<mi>k</mi>
<mi>i</mi>
</msubsup>
<mo>)</mo>
</mrow>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mi>i</mi>
<mo>=</mo>
<mn>1,2</mn>
<mo>,</mo>
<mo>&CenterDot;</mo>
<mo>&CenterDot;</mo>
<mo>&CenterDot;</mo>
<mi>N</mi>
<mo>,</mo>
</mrow>
并且归一化权值则得k时刻未知参数x最小均方估计为：
(3)重采样；得到新的粒子集合
(4)预测；利用状态方程f预测未知参数
(5)时刻k＝k+1，转到步骤(2)；
(二)滤波降噪之后，采用快速谱峭度的算法计算得到峭度图，过程如下：
(6)构造一个截止频率fc＝1/8+ε的低通滤波器h(n)，其中fs＝1,ε＞0；低通滤波器的截止频率增大ε是为了消除过度带的影响；基于h(n)构造准解析低通滤波器h0(n)和准解析高通滤波器h1(n)，其分析频带分别为[0,0.25]、[0.25,0.5]，即：
<mfenced open='{' close=''>
<mtable>
<mtr>
<mtd>
<msub>
<mi>h</mi>
<mn>0</mn>
</msub>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>h</mi>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>e</mi>
<mrow>
<mi>j&pi;n</mi>
<mo>/</mo>
<mn>4</mn>
</mrow>
</msup>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>h</mi>
<mn>1</mn>
</msub>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>h</mi>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>e</mi>
<mrow>
<mi>j&pi;n</mi>
<mo>/</mo>
<mn>4</mn>
</mrow>
</msup>
</mtd>
</mtr>
</mtable>
</mfenced>
(7)分别以h0(n)和h1(n)进行低通和高通滤波，对滤波结果做2倍降采样；如此迭代进行，获得滤波器树及相应的滤波结果；其中表示第k层第i个滤波器产生的滤波结果，其中i取值为0到2k-1；同时也是信号在中心频带fi＝(i+2-1)2k-1和带宽(Δf)k＝2k-1的复包络，2倍降采样是保证滤波器每一层中的数据长度与原始数据相同；
(8)对于每一个滤波器结果，按照式f≠0计算峭度：
<mrow>
<mi>K</mi>
<mo>=</mo>
<mfrac>
<mrow>
<mi>E</mi>
<mo><</mo>
<msup>
<mrow>
<mo>|</mo>
<msubsup>
<mi>c</mi>
<mi>k</mi>
<mi>i</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>)</mo>
</mrow>
<mo>|</mo>
</mrow>
<mn>4</mn>
</msup>
<mo>></mo>
</mrow>
<mrow>
<mo>{</mo>
<mi>E</mi>
<mo><</mo>
<msup>
<mrow>
<mo>|</mo>
<msubsup>
<mi>c</mi>
<mi>k</mi>
<mi>i</mi>
</msubsup>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>)</mo>
</mrow>
<mo>|</mo>
</mrow>
<mn>2</mn>
</msup>
<msup>
<mrow>
<mo>></mo>
<mo>}</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>-</mo>
<mn>2</mn>
</mrow>
(9)将所有的谱峭度汇总，从而得到“快速谱峭度图”；
(10)选择角度最大的频带做解调分析，计算功率谱得到频谱图。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510031512.9A CN104677632B (zh) | 2015-01-21 | 2015-01-21 | 利用粒子滤波与谱峭度的滚动轴承故障诊断方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510031512.9A CN104677632B (zh) | 2015-01-21 | 2015-01-21 | 利用粒子滤波与谱峭度的滚动轴承故障诊断方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104677632A CN104677632A (zh) | 2015-06-03 |
| CN104677632B true CN104677632B (zh) | 2018-04-10 |
Family
ID=53312984
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510031512.9A Active CN104677632B (zh) | 2015-01-21 | 2015-01-21 | 利用粒子滤波与谱峭度的滚动轴承故障诊断方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104677632B (zh) |
Families Citing this family (16)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106017926A (zh) * | 2016-05-13 | 2016-10-12 | 山东理工大学 | 基于变模态分解的滚动轴承故障诊断方法 |
| CN106168538B (zh) * | 2016-06-29 | 2018-07-03 | 潍坊学院 | 一种滚动轴承的itd、谱峭度和平滑迭代包络分析方法 |
| CN106053059B (zh) * | 2016-06-29 | 2018-03-27 | 潍坊学院 | 一种基于内秉时间尺度分解滤波的包络分析方法 |
| CN106053061B (zh) * | 2016-06-29 | 2018-03-13 | 潍坊学院 | 一种基于非线性模式分解和谱峭度的包络分析方法 |
| CN105954030B (zh) * | 2016-06-29 | 2018-03-23 | 潍坊学院 | 一种基于内秉时间尺度分解和谱峭度的包络分析方法 |
| CN106441893B (zh) * | 2016-09-22 | 2018-08-10 | 北京邮电大学 | 火车滚动轴承故障与杂质振动区分方法 |
| CN108646296A (zh) * | 2018-05-16 | 2018-10-12 | 吉林大学 | 基于自适应谱峭度滤波器的沙漠地震信号噪声消减方法 |
| CN109029977B (zh) * | 2018-07-12 | 2019-12-31 | 福州大学 | 一种基于vmd-amckd的行星齿轮箱早期故障诊断方法 |
| CN109446625B (zh) * | 2018-10-22 | 2022-07-26 | 南昌航空大学 | 一种基于贝叶斯推理的直升机动部件动态阈值计算方法 |
| CN109670243B (zh) * | 2018-12-20 | 2020-11-24 | 华中科技大学 | 一种基于勒贝格空间模型的寿命预测方法 |
| CN111122161B (zh) * | 2019-12-24 | 2021-06-08 | 北京交通大学 | 一种基于fast kurtogram和深度残差学习的变工况轴承故障诊断方法 |
| CN110895687B (zh) * | 2019-12-24 | 2023-04-07 | 安徽工业大学 | 一种基于均变相位经验模态分解和快速谱峭度的轴承故障诊断方法 |
| CN111444893A (zh) * | 2020-05-06 | 2020-07-24 | 南昌航空大学 | 一种矿井提升机主轴装置的故障诊断方法 |
| CN113340995B (zh) * | 2021-05-11 | 2024-05-07 | 西安交通大学 | 一种激光冲击强化缺陷实时检测的声发射信号频段选择方法 |
| CN114295365B (zh) * | 2022-03-09 | 2022-07-08 | 深圳市信润富联数字科技有限公司 | 插针机精度异常监测方法、系统、设备及存储介质 |
| CN117056675A (zh) * | 2023-08-16 | 2023-11-14 | 西南石油大学 | 一种基于小波变换和粒子滤波结合的漏磁信号降噪方法 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| DE102012220222A1 (de) * | 2012-11-07 | 2014-05-08 | Siemens Aktiengesellschaft | Vorrichtung und Verfahren zur Zustandsüberwachung eines Wälzlagers |
| CN103424258A (zh) * | 2013-08-06 | 2013-12-04 | 昆明理工大学 | 一种滚动轴承的故障诊断方法 |
| CN104236905A (zh) * | 2014-08-26 | 2014-12-24 | 中国直升机设计研究所 | 一种轴承故障诊断方法 |
-
2015
- 2015-01-21 CN CN201510031512.9A patent/CN104677632B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN104677632A (zh) | 2015-06-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104677632B (zh) | 利用粒子滤波与谱峭度的滚动轴承故障诊断方法 | |
| CN108760327B (zh) | 一种航空发动机转子故障的诊断方法 | |
| CN109145727A (zh) | 一种基于vmd参数优化的轴承故障特征提取方法 | |
| Liao et al. | An automatic filtering method based on an improved genetic algorithm—with application to rolling bearing fault signal extraction | |
| CN112183259B (zh) | 一种基于ceemd与峭度加权平均阈值去噪的滚动轴承故障诊断方法 | |
| CN112945546B (zh) | 一种齿轮箱复杂故障精确诊断方法 | |
| CN110940522A (zh) | 强背景噪声下轴承故障周期性脉冲稀疏分离与诊断方法 | |
| Kim et al. | A new prognostics approach for bearing based on entropy decrease and comparison with existing methods | |
| Li et al. | The contact state monitoring for seal end faces based on acoustic emission detection | |
| CN110987433B (zh) | 一种基于高频信号特征幅值的轴承故障预警方法 | |
| Liao et al. | A vibration signal filtering method based on KL divergence genetic algorithm–with application to low speed bearing fault diagnosis | |
| Wang et al. | The diagnosis of rolling bearing based on the parameters of pulse atoms and degree of cyclostationarity | |
| Rehab et al. | Bearings fault detection using hidden Markov models and principal component analysis enhanced features | |
| CN115452377B (zh) | 一种基于ceemdan与自适应去噪的滚动轴承故障诊断方法 | |
| Liu et al. | Fault diagnosis method of rolling bearing based on the multiple features of LMD and random forest | |
| CN113639985B (zh) | 一种基于优化故障特征频谱的机械故障诊断与状态监测方法 | |
| CN114021275B (zh) | 一种基于深度卷积模糊系统的滚动轴承故障诊断方法 | |
| Wang et al. | A two-stage compression method for the fault detection of roller bearings | |
| CN113869358A (zh) | 基于循环相关熵和一维浅卷积神经网络的轴承故障诊断方法 | |
| Aijun et al. | A novel approach of impulsive signal extraction for early fault detection of rolling element bearing | |
| Meng et al. | Rolling bearing fault diagnosis method based on MCMF and SAIMFE | |
| Duan et al. | Bearing fault diagnosis based on state-space principal component tracking filter algorithm | |
| CN112836583B (zh) | 一种风力机故障诊断方法 | |
| Lin et al. | Fault diagnosis of rolling bearings using multifractal detrended fluctuation analysis and Mahalanobis distance criterion | |
| Zhao et al. | A new intelligent weak fault recognition framework for rotating machinery |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |