轴承故障诊断——领域知识

林文韬0327

已于 2024-08-26 01:41:33 修改

阅读量1.8k

点赞数 10

文章标签：笔记

于 2024-08-02 20:51:45 首次发布

本文链接：https://blog.csdn.net/qq_58516607/article/details/140877960

版权

先把故障机理以及相应的基础的现代信号处理方法搞一搞，然后再开始上机器学习/深度学习算法。

然后，开始时域，频域分析，关于时域特征；了解下频谱，包络谱，平方包络谱，倒谱，频谱细化分析等等；还有各种熵值特征，分形特征。

发现时域，频域分析在面对（变速，变负载等等）复杂工况下不够用，开始进行时频谱分析，这方面的内容可太丰富了，找感兴趣的方向即可。

先说缺乏数学理论基础的，经验模态分解及其众多变体（EMD,EEMD,CEEMD，还有各种抑制模态混叠和端点效应的方法），局部均值分解LMD，局部特征尺度分解等等。然后开始分数阶傅里叶变换，短时傅里叶变换STFT，小波分析及众多变体（双树复小波，可调Q因子小波，多小波等等），WVD分布，变分模态分解，经验小波变换，辛几何模态分解，同步挤压变换，高阶同步挤压变换等等。

最后，就可以上机器学习/深度学习模型了，支持向量机SVM,随机森林，其他集成模型，线性回归，隐马尔可夫模型，字典学习，主成分分析，贝叶斯方法，堆栈自编码器及众多变体，深度信念网络DBN，卷积神经网络CNN，长短时记忆网络LSTM，生成对抗网络GAN等等，找感兴趣的方向即可。

作者：哥廷根数学学派
链接：https://zhuanlan.zhihu.com/p/550433862

绪论

机械故障诊断技术是一种利用各种检测方法，测取并分析、处理机械设备在运行中的状态信息，确定其整体或局部是否正常，早期发现故障及其原因，并能预报故障发展趋势的技术。

①机械运行过程是动态随机过程。即使同型号机械设备由于装配、安装及工作条件上的差异，也往往导致机器的工况状态及故障模式改变。

②从系统特性来看，机械设备都是由成百上千个零部件装配而成，零部件间相互偶合，决定了机械设备故障的多层次性。故障与现象之间没有简单的对应关系。

③如果仅从某一个参数或某一个侧面去分析而做出判断，一般很难做出正确的决策。应该从随机过程的理论出发，运用各种现代多学科融合的分析工具，综合判断机械的故障现象、属性、形成及其发展趋势。

研究内容

机械状态信号的测量、机器状态或失效信息的提取、状态识别、诊断决策等几个步骤.

诊断方法

1 基于信号处理的方法：对测得的振动量在时域、频域和时频域进行特征分析，用于确定机器各种故障的类型和性质。

2 基于知识的方法：通过处理测量到的输入输出信号来实现故障诊断。必须拥有大量的关于系统故障的经验知识，具有实测到的大量的各类故障样本数据。(机器学习)

3 基于解析模型的方法: 建立被诊断对象的较为精确的数学模型.

诊断手段

1 振动诊断，对机器主要部位的振动值（位移、速度、加速度、转速、相位值）进行测定，在时域、频域、时-频域进行特征分析。

2 噪声诊断

3 温度、压力等常规参数（价格便宜、形式多样）

4 无损诊断（超声波探伤、X射线探伤、渗透探伤、磁粉探伤），一般用于材料表面或内部。

5 油液分析（光谱分析、铁谱分析）

多种故障诊断方法融合的复合诊断，多元传感器信息的融合。

信号处理

时域平均方法对于滚动轴承的诊断基本上是没有用的。因为滚动轴承特征频率一般不在轴的旋转频率上，时域平均后振动特征被消除掉了。另外，由于滚动体与内、外圈之间存在相对滑动．故障点所产生的冲击振动重复性不好。因此，这种信号不适合时域平均方法。

倒谱是英文 cepstrum 的直译．也称二次谱和对数功率谱. 倒谱分析是检测复杂谱图中周期分拔的有效工具，可将振动信号功率谱图上的众多边带谱线简化为单根谱线，具有信息凝聚作用。

轴承相关概念

轴系是由轴、轴承、安装于轴上的传动体、密封件及定位组件组成.

滚动轴承

相关参数

轴的分类

分类

代号

载荷分布

轴承边缘的载荷和应力

失效形式

滑动轴承

数据集说明

包括十种轴承故障的振动加速度数据：正常（N）、内圈故障（IF）、外圈故障（of）、滚珠故障（BF）、保持架故障（CF）……

文件格式为.mat，每种类型的轴承故障数据都包含两个文件夹：

1.rpm1000_var_load：包含1000 rpm固定速度下七种不同负载力的数据。负载力为0、200、400、600、800、1000和1200 N。

2.var_speed_f600：包含固定负载力为600 N时七种不同速度的数据。速度为400、600，800、1000、1200、1400和1600 rpm。

特定操作条件下的每个文件夹包含10个.mat文件，代表10个数据集合。

该实验平台的数据采样频率为10kHz，采样时间为20至25秒。四个振动传感器用于数据采集，数据尺寸按1-4排的顺序排列，代表电机、轴承座、轴承平行上端和轴承平行侧端。

因此，每个.mat文件中的数据维度为（4，200000+），其中第一个维度对应于上述四个不同的振动测量位置，每个位置都有超过200000个数据点。

幅域分析

幅值概率密度密度函数

表示信号的幅值落在某一指定区间内的概率。

def interval_num_count(data, low, high):
    '''
    fun: 统计一维数据data落入某一个区间[low, high]内的数量
    param low: 区间下限
    param high: 区间上限
    return count_num: 落入某一个区间[low, high]内的数量
    '''
    count_num = 0
    for i in range(len(data)):
        if data[i]>low and data[i]<high:
            count_num += 1
    return count_num

def plt_amp_prob_density_fun(data, n):
    '''
    fun: 绘制幅值概率密度函数
    param data: 输入数据，1维array
    param n: 分割成段数的数量
    return: 绘制幅值概率密度函数
    '''
    xt = data - np.mean(data)
    max_value = np.abs( xt[np.argmax( np.abs(xt) )] ) #
    count_num = []
    for i in range(n):
        interval = max_value*2/n              # 区间长度为interval_len
        low = -max_value + i*interval         # 区间下限
        high = -max_value + (i+1)*interval    # 区间上限
        count = interval_num_count(data=xt, low=low, high=high)  # 统计落入该区间的幅值个数
        count_num.append(count)
    plt.bar(x=range( len(count_num) ), height=count_num)  # 绘制柱状图
    plt.show()

随机信号的概率论密度函数符合正态分布规律，确定校信号则呈现盆形曲线。一般故障信号多是随机信号和确定性信号的混合，当信号概率密度函数的正态分布曲线上端出现盆形漏斗时，往往预示着故障征兆。

统计参数

均值

${{\mu }_{\text{x}}}=\underset{T\to \infty }{\mathop{\lim }}\,\frac{1}{T}\int_{0}^{T}{x(t)dt}$

均方值

$\psi _{x}^{2}=\underset{T\to \infty }{\mathop{\lim }}\,\frac{1}{T}\int_{0}^{T}{{{x}^{2}}(t)dt}$

方差

$\sigma _{x}^{2}=\underset{T\to \infty }{\mathop{\lim }}\,\frac{1}{T}\int_{0}^{T}{{{\left( x(t)-{{\mu }_{x}} \right)}^{2}}dt}$

若信号的均值为零，则均值值等于方差。方差的正算术平方根成为标准差。

$\sigma _{x}^{2} = \psi _{x}^{2} - \mu _{x}^2$

实际工程信号长度有限，故采用：

有量纲幅域参数

均值

$\bar{x}=\frac{1}{N}\sum\limits_{i=1}^{N}{{{x}_{i}}}$

反应信号中静态部分，对诊断不起作用，但对计算其他参数有很大影响。零均值处理：一般在计算时应该先从数据中除去均值，剩下对诊断有用的动态部分。所有动态信号必须进行。

频谱分析等进一步处理时，也必须进行零均值处理。

均方根值（有效值）

${{x}_{rms}}=\sqrt{\frac{1}{N}\sum\limits_{i=1}^{N}{x_{i}^{2}}}$

表示信号的能量，多用于品佳振动等级或烈度。

最大值

${x_{\max }} = \max \left\{ {\left| {{x_i}} \right|} \right\}$

峰峰值

${x_{p - p}} = \max \left\{ {\left| {{x_i}} \right|} \right\} - \min \left\{ {\left| {{x_i}} \right|} \right\}$

歪度

$\alpha = {1 \over N}\sum\limits_{i = 1}^N {x_i^3}$

表示信号的幅值概率密度函数对纵坐标的不对称性。值越大，越不对称。

峭度

$\beta = {1 \over N}\sum\limits_{i = 1}^N {x_i^4}$

当K=3定义为分布曲线具有正常峰度（即零峭度）；当K＞3时，分布曲线具有正峭度。由此可知，当标准差 σt小于正常状态下的标准差，即观测值的分散程度较小时，K增大，此时正态分布曲线峰顶的高度高于正常正态分布曲线，故称为正峭度。当K＜3时，分布曲线具有负峭度，当标准差σt大于正常状态下的标准差，即观测值的分散程度较大时，K减小，此时正态分布曲线峰顶的高度低于正常正态分布曲线，故称为负峭度。

在轴承无故障运转时，由于各种不确定因素的影响，振动信号的幅值分布接近正态分布，峭度指标值K≈3;随着故障的出现和发展，振动信号中大幅值的概率密度增加，信号幅值的分布偏离正态分布，正态曲线出现偏斜或分散，峭度值也随之增大。峭度指标的绝对值越大，说明轴承偏离其正常状态，故障越严重，如当其K＞8时，则很可能出现了较大的故障。

无量纲幅域参数

有量纲参数的大小均与信号振动的绝对幅值有关，也就是和振动产生的工作条件有关，不同工作条件下的有量纲时域统计特征不可比、为此构造了无量纲时域统计特征。

波形指标

${s_f} = {{{x_{rms}}} \over {\left| {\bar x} \right|}}$

峰值指标

${C_f} = {{{x_{\max }}} \over {{x_{rms}}}}$

脉冲指标

${I_f} = {{{x_{\max }}} \over {\left| {\bar x} \right|}}$

裕度指标，是无量纲的歪度指标，表示信号的幅值概率密度函数对纵坐标的不对称性。

$C{L_f} = {{{x_{\max }}} \over {{x_r}}}$

峭度指标

${K_r} = {\beta \over {x_{rms}^4}}$

其中：

平均幅值

$\left| x \right| = {1 \over N}\sum\limits_{i = 1}^N {\left| {{x_i}} \right|}$

方根幅值

${x_{r}} = {\left[ {{1 \over N}\sum\limits_{i = 1}^N {\sqrt {\left| {{x_i}} \right|} } } \right]^2}$

同时用峭度指标与有效值(均方根值)进行故障监测令以兼顾敏感性与稳定性。
连续监测可发现峭度指标的变化趋势，当指标值上升到顶点开始下降时，就要密切注意是否有故障发生。

这些参数仅取决于信号的概率密度函数，而和频率和幅值无关。

信号的相关分析

自相关函数主要用于检测混淆在随机信号中的周期信号成分。

周期信号的自相关函数会按原频率重复出现，随机信号在时间位移 $\tau$ 稍大时，由于自身的相乘相消作用，自相关函数很快趋于0（假设均值为0）.

实际工程信号多数是随机噪声和确定性周期信号的混合体。一般情况下，周期信号和故障特征有关，随机噪声对诊断无用。

可以利用随机信号的自相关函数迅速衰减而周期函数不衰减的特性，在自相关图的右侧部分测取信号的周期。

自相关函数是从干扰噪声中找出周期信号或瞬时信号的重要手段，延长变量 $\tau$ 的取值，就可将信号中的周期分量暴露出来。

直接计算法计算量大，近代信号分析不再采用。利用自相关函数与功率谱密度函数的关系，采用快速傅里叶算法实现。

互相关函数的应用：

相关分析在机械设备故障诊断和振动控制中最直接的应用是传递问题，其中包括传递路径的识别和故障源的识别这两类问题，间接的应用是测速和定位。

频域分析

傅里叶级数

周期信号 $x(t)=x(t) \pm nT$ 可以展开成傅里叶级数——用无限个三角函数去拟合。

$f_{0}$ 为基频， ${f_0} = {1 \mathord{\left/ {\vphantom {1 {{T_0}}}} \right. \kern-\nulldelimiterspace} {{T_0}}}$ ， $T_{0}$ 为周期。

欧拉公式

$\cos \theta = {1 \over 2}\left( {{e^{j\theta }} + {e^{ - j\theta }}} \right)$

$\sin \theta = {1 \over {2j}}\left( {{e^{j\theta }} - {e^{ - j\theta }}} \right) = - {j \over 2}\left( {{e^{j\theta }} - {e^{ - j\theta }}} \right)$

*****************************************************************************************************

$x\left( t \right) = {{{a_0}} \over 2} + \sum\limits_{n = 1}^\infty {\left( {{a_n}\cos 2\pi n{f_0}t + {b_n}\sin 2\pi n{f_0}t} \right)}$

${{{a_0}} \over 2} + \sum\limits_{n = 1}^\infty {\left( {{{{a_n}} \over 2}\left( {{e^{2\pi jn{f_0}t}} + {e^{ - 2\pi jn{f_0}t}}} \right) + {{ - j{b_n}} \over 2}\left( {{e^{2\pi jn{f_0}t}} - {e^{ - 2\pi jn{f_0}t}}} \right)} \right)}$

${{{a_0}} \over 2} + \sum\limits_{n = 1}^\infty {{1 \over 2}\left( {{a_n} - j{b_n}} \right){e^{2\pi jn{f_0}t}}} + \sum\limits_{n = 1}^\infty {{1 \over 2}\left( {{a_n} + j{b_n}} \right){e^{ - 2\pi jn{f_0}t}}}$

*****************************************************************************************************

$\sum\limits_{n = 1}^\infty {{1 \over 2}\left( {{a_n} + j{b_n}} \right){e^{ - 2\pi jn{f_0}t}}} \mathop = \limits^{n = - n} \sum\limits_{n = - \infty }^1 {{1 \over 2}\left( {{a_{ - n}} + j{b_{ - n}}} \right){e^{2\pi jn{f_0}t}}} = \sum\limits_{n = - \infty }^1 {{1 \over 2}\left( {{a_n} - j{b_n}} \right){e^{2\pi jn{f_0}t}}}$

*****************************************************************************************************

${{{a_0}} \over 2} + \sum\limits_{n = 1}^\infty {{1 \over 2}\left( {{a_n} - j{b_n}} \right){e^{2\pi jn{f_0}t}}} + \sum\limits_{n = - \infty }^1 {{1 \over 2}\left( {{a_n} - j{b_n}} \right){e^{2\pi jn{f_0}t}}}$

${{{a_0}} \over 2} + \sum\limits_{n = - \infty }^\infty {{1 \over 2}\left( {{a_n} - j{b_n}} \right){e^{2\pi jn{f_0}t}}}$

*****************************************************************************************************

$x\left( t \right) = \sum\limits_{ - \infty }^\infty {{C_n}{e^{2\pi jn{f_0}t}}}$

${C_n} = {{{a_n} - j{b_n}} \over 2}$

${1 \over 2} \times {2 \over {{T_0}}}\int_0^{{T_0}} {x\left( t \right)\cos 2\pi n{f_0}t} dt + {{ - j} \over 2} \times {2 \over {{T_0}}}\int_0^{{T_0}} {x\left( t \right)\sin 2\pi n{f_0}t} dt$

${1 \over {{T_0}}}\int_0^{{T_0}} {x\left( t \right)\left( {{{{e^{2\pi jn{f_0}t}} + {e^{ - 2\pi jn{f_0}t}}} \over 2} - j{{{e^{2\pi jn{f_0}t}} - {e^{ - 2\pi jn{f_0}t}}} \over {2j}}} \right)dt}$

****************************************************************************************************

$x\left( t \right) = \sum\limits_{ - \infty }^\infty {{C_n}{e^{2\pi jn{f_0}t}}}$

${C_n} = {1 \over {{T_0}}}\int_0^{{T_0}} {x\left( t \right){e^{ - 2\pi jn{f_0}t}}dt}$

是个复变量，本身可以表示信号的幅值和相位，称为有限区间上信号的离散频谱。

取模值，傅里叶级数的幅值谱，每条谱线只出现在基波频率的整数倍上，不存在非整数倍的频率分量。可见，周期信号的频谱是离散的，且随着谐波次数的增加，谐波幅值是逐渐下降的。

根据这些特征，在频域判别信号是周期还是非周期的。

傅里叶积分（傅里叶变换）

零均值化处理的原因

采样定理

离散傅里叶变换

第一步，时域离散

第二步，时域截断

选择一个窗函数（区间宽度），获得有限个采样点。这样做的结果是原有信号的频谱出现褶皱，称为频谱泄露。

第三步，频域离散

计算出的频谱是连续的，不能用计算机进行。必须进行采样离散化，间隔为截断函数区间宽度的倒数。

时域采样会引起频域的周期化，频域采样会引起时域的周期化。称为时域延拓。

只有前一半计算点（ N/2 ）有意义.

直接的离散傅里叶变换

频率分析精度

整周期采样

DFT 方法用一定长度的矩形窗函数截取任意一段非周期函数，将其左右时域延拓后，形成一个新的周期函数，然后对其进行傅里叶变换。

当这个窗长度与被测信号的周期成整数倍时。截取的信号的首尾点的幅值是相同的，如果左右延拓的话，就会衔接上，在整个时间域上形成一个和原来信号完全相同的信号，所以其幅值谱和理论值相同。

而当截取周期与信号的周期不成整数倍时，时域延拓的结果就会出现断点，引起的波形畸变，这样在频谱中必然会衍生出很多的频率成分去弥补这些断点。

窗函数

傅里叶频谱分析

确定性信号是指可以用明确的数学关系或者图表描述的信号。若信号被表示为一确定的时间函数，对于指定的某一时刻，可以确定一相应的函数值，这种信号被称为确定性信号。

确定性信号的傅里叶谱是个复数，包含实部、虚部，可绘制实频、虚频、幅频、相频、幅-相。

而工程和实际中接触的信号往往都是随机信号，随机信号不满足FFT变换的条件（随机信号不满足傅里叶变换绝对可积条件），不能直接进行FFT变换，顾一般不采用幅值谱和相位谱进行分析。

能量信号（能量有限，平均功率为0，代表波形：一个孤零零的方波）
功率信号（能量无限，平均功率有限，代表波形：一个无限延伸的正弦波）
非功率非能量信号（无穷能量+无穷功率，代表波形：一个无限延伸的单调波形）

若信号的总能量是有限的，可以用能量谱或幅值谱来考察。能量信号的能量是一个非无穷的正值，这时候就可以把能量作为考察能量信号的有效量纲，而且由于能量信号的能量有限，能量信号的平均功率肯定是无限接近于0的，这时候从平均功率的角度去考察能量信号就没有意义。

若信号的总能量是无限的，但单位时间内的能量是有限的，则用功率谱密度函数考察。功率信号的能量是无穷的，从能量角度去考察就没有意义了，但是功率信号的平均功率肯定是个非零值，这时候选择平均功率作为考察的量纲就是合理的。(信号做功消耗能量，平均功率是信号所做的功与无穷时间的比值）

随机信号的功率谱

当时间趋于无穷时，平均功率

${P_x} = \mathop {\lim }\limits_{T \to \infty } {1 \over T}\int_{ - {T \over 2}}^{{T \over 2}} {{x^2}(t)dt}$

${P_x} = \mathop {\lim }\limits_{T \to \infty } {1 \over T}\int_{ - {T \over 2}}^{{T \over 2}} {{{\left| {X(f)} \right|}^2}df} = \int_{ - \infty }^\infty {\mathop {\lim }\limits_{T \to \infty } {{{{\left| {X(f)} \right|}^2}} \over T}df}$

*************************************************************************************************************

平均功率谱密度函数，具有单位频率轴上的平均功率量纲，简称功率谱。

${S_x}(f) = \mathop {\lim }\limits_{T \to \infty } {{{{\left| {X(f)} \right|}^2}} \over T}$

*************************************************************************************************************

fs = 10000  #采样频率
length = len(dataarr)  #采样点数
t = np.linspace(0, length/fs, length)

fig, axs = plt.subplots(2, 1)
axs[0].plot(t, dataarr)
axs[0].set_xlabel('time(s)')
axs[0].set_ylabel('Amplitude')

p = fftpack.fft(dataarr.real)
power = np.abs(p) ** 2 / length  #计算功率谱
power_db = 10 * np.log10(power)  # 转化为分贝
f = np.fft.fftfreq(p.size, 1/fs)  #计算频率

axs[1].plot(f[0:length//2], power_db[0:length//2], 'r')
axs[1].set_xlabel('Frequency(Hz)')
axs[1].set_ylabel('Magnitude(dB)')

fig.tight_layout()  #调整子图间距

python实现