首页/文章/ 详情

机械故障诊断信号的幅域分析 - 时域统计特征 | 在CWRU和IMF轴承数据集上实战

故障诊断与python学习

1年前浏览1353

m脑袋chegnz

背景及摘要

很多初学者刚接触故障诊断可能觉得很简单，套用深度学习模型进行训练，分类准确率达到99%即可。

在写论文时，这样的确没问题。但是在工作或者在解决实际问题时，几乎是用不上。

最近在读时献江老师等人的《机械故障诊断及典型案例解析》这本书时，发现其内容真的很丰富，故障诊断方法讲解的详细全面，且有实际丰富的案例，更加容易消化理解各方法特点及其应用场景。

信号波形是某种物理量随时间变化的关系。机械振动信号处理的基本方法有幅域分析、时域分析和频域分析。

在幅值上的各种统计处理通常称为幅域分析。仅对波形的幅值用简单的方法进行统计分析，如计算波形的最大值、平均值和有效值，或研究时域波形幅值的概率分布形式等，

信号在时间域内的变换或统计分析称为时域分析。频域分析是确定信号的频率结构，即信号中包含哪些频率成分，分析的结果是以频率为自变量的各种物理量的谱线或曲线。

不同的分析方法是从不同的角度观察、分析信号，使信号处理的结果更加利于故障分析与诊断。

本节将介绍幅域分析方法中的时域统计特征，包括其数学公式定义及CWRU和IMF轴承数据实战案例。

代码位置：

‍https://github.com/HappyBoy-cmd/fault_diagnosis_signal_processing

参考资料：

书籍：机械故障诊断及典型案例解析（第2版，时献江）
学位论文：细纱机罗拉轴承故障诊断方法研究

1. 摘要

2. 有量纲时域统计特征计算公式及物理意义

3. 无量纲时域统计特征计算公式及物理意义

4. 模拟数据代码实战

4.1 导入包

4.2 生成模拟正弦数据

4.3 绘制幅值概率密度函数

4.4 计算时域统计特征

5. CWRU数据实战

5.1 内圈故障

5.2 正常状态

5.3 外圈故障

5.4 滚动体故障

6. IMF加速寿命实验轴承数据

6.1 IMF数据简介

6.2 单次数据尝试

6.3 整个寿命周期数据分析

7. 总结

8. 课后思考

1、摘要

信号的幅域分析也称统计特征分析（也叫幅域参数分析），主要利用振动信号的幅值统计特征来进行分析和诊断。应用比较广泛的有均方根值、峰值指标、波形指标和峭度等指标。信号的幅域分析也属于时域分析，和相关分析等时域分析方法不同，时域分析不考虑原始信号的时序，仅与信号的幅值大小及分布有关。时域统计特征包括有量纲时域统计特征和无量纲时域统计特征两大类。本文将介绍有量纲时域统计特征和无量纲时域统计特征，并在CWRU轴承和IMF(辛辛那提大学轴承数据)全寿命周期轴承实验数据上进行实战。

2、有量纲时域统计特征计算公式及物理意义

时域统计特征参数可以通过3种方式计算。

第1种：随机信号的时域统计特征与幅值概率密度函数有密切关系，对于各态历经的平稳信号，可以由幅值概率密度函数计算如下统计参数

均值：

均方根值：

方差：

绝对平均值：

方根幅值：

歪(偏)度：

峭度：

第2种：以上参数计算需要用到幅值概率密度函数，不易计算。实际上对于各态历经的平稳随机信号，可以直接利用单个样本进行计算，公式如下：（其实该种计算公式适合当为连续函数时）

均值：

均方根值：

方差：

绝对平均值：

方根幅值：

歪(偏)度：

峭度：

其中：

信号的均值反映信号中的静态部分，多数情况下表示振动的平衡位置；
均方根值反映信号的能量大小，相当于电学中的有效值，多用于评价振动等级或烈度；
方差仅反映了信号中的动态部分，反映振动信号以平衡位置为中心的幅值变化程度，若信号的均值为零，则均方值（均方根值的平方）等于方差。
歪度表示信号的幅值概率密度函数对纵坐标的不对称性，越大，越不对称（可参见后面的图1)。
峭度表示正态分布曲线的性状，当较小时表示分布曲线瘦而高，成为正峭度；当较大时，分布曲线具有负峭度，此时正态分布曲线峰顶的高度低于正常正态分布曲线（可参见后面的图2)。

图1 裕度指标的物理意义

图2 峭度指标的物理意义

第3种：以上参数是理论上的统计真值，实际工程信号采样长度有限，一般采用下述的估计值。以后不提示统计真值和估计值的区别，实际计算过程均为有限长度的估计值，有时为了说明问题方便，也常常使用统计真值的理论公式。

对于时间序列信号，有量纲时域统计指标估计值的计算公式如下：（其实该种计算公式适合当为离散函数时）

均值：

方差：

均方根值（有效值）：

均方值：

均对平均幅值：

方根幅值：

最大值：

最小值：

峰值：

峰峰值：

歪(偏)度：

峭度：

3、无量纲时域统计特征计算公式及物理意义

有量纲时域统计特征的大小与信号（振动）绝对幅值有关，也就是和振动产生的工作条件有关，不同工作条件下的有量纲时域统计特征不可比、为此构造了无量纲时域统计特征。对于时间序列信号，无量纲时域统计特征的计算公式如下：

波形指标：

峰值指标：

脉冲指标：

裕度指标：

峭度指标：

偏度指标：

其中，

裕度指标是无量纲的歪度指标，表示信号的幅值概率密度函数对纵坐标的不对称性。如果越大，越不对称，且不对称有正（右偏移）负（左偏移）之分，如图1所示。旋转机械等设备的振动信号由于存在某一方向的摩擦或碰撞，或者某一方向的支撑刚度较弱，会造成振动波形的不对称，使裕度指标增大。
峭度指标对大幅值最为敏感，当大幅值出现的概率增加时，值会迅速增加，这对探测信号中含有脉冲的故障特别有效。峭度指标的物理意义如图3所示。时定义分布曲线具有正常峰度（即零峭度）；当时，分布曲线具有正峭度，此时正态分布曲线峰顶的高度高于正常正态分布曲线，故称为正峭度。当时，分布曲线具有负峭度，此时正态分布曲线峰顶的高度低于正常正态分布曲线，故称为负峭度。

图3 滚动轴承振动幅值概率密度分布图

图3 为滚动轴承的振动幅值概率密度分布图。

实线为正常时，幅值概率密度函数近似为正态分布；
虚线为发生剥落时，此时幅值概率密度函数呈现头部窄、底部宽的形式，值较大表明信号中冲击成分幅值增大（底部宽），但是能量不大（值小），即系统处于剥落故障开始发生的时刻。

另外，必须着重指出，信号的均值反映信号中的静态部分，一般对诊断不起作用，但对计算上述参数有很大影响。所以，一般在计算时应先从数据中去除均值，保留对诊断有用的动态部分，这一过程称为零均值化处理（也叫去直流分量），其计算方法如下：

假设原始时间序列信号，其均值，则去均值后的新时间序列计算式为：
表1 为齿轮振动信号的无量纲幅域诊断参数。新齿轮经过运行产生了疲劳剥落故障，振动信号中有明显的冲击脉冲，各幅域参数中除了波形参数外，均有明显上升。

表1 齿轮振动信号的无量纲时域特征

峭度指标、裕度指标和脉冲指标对冲击脉冲型故障比较敏感。
当早期发生故障时，大幅值的脉冲还不是很多，均方根值变化不大，上述参数已有增加。
当故障逐步发展时，它们上升较快；但上升到一定程度后，由于分母上的有效值增大，这些指标反而会逐步下降。
这表明这些参数对早期故障有较高敏感性，但稳定性不很好。
均方根值则相反，虽然对早期故障不敏感．但稳定性好，随着故障发展单调上升。

图4为某滚动轴承振动信号的峭度指标和有效值随轴承疲劳试验时间的变化过程，可见，两个指标的变化符合上述规律。因此，要想取得较好的故障监测效果，一般可以采取以下措施：

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

图4 峭度指标和有效值随轴承疲劳试验时间的变化过程

4、模拟数据代码实战

4.1 导入包

import scipy.io as scioimport numpy as npimport matplotlib.pyplot as pltfrom scipy import signal, fftpack, stats

from matplotlib import rcParams
config = {    "font.family": 'serif', # 衬线字体    "font.size": 10, # 相当于小四大小    "font.serif": ['SimSun'], # 宋体    "mathtext.fontset": 'stix', # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大    'axes.unicode_minus': False # 处理负号，即-号}rcParams.update(config)

4.2 生成模拟正弦数据

fs = 1000f = 10n = np.arange(5120)t = n/fsx = 3 * np.sin(2 * np.pi * n * f/fs) + np.random.uniform(0,1, len(n)) #正弦信号3*sin(2*pi*10*t)+噪声信号

plt.figure(figsize=(12,2))plt.plot(t, x)

上图为加了噪声的正弦信号。

4.3 绘制幅值概率密度函数

首先定义概率密度函数函数，具体见机械故障诊断信号的幅域分析 - 幅值概率密度函数 | 基于python的代码实现，在CWRU轴承数据上实战

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    max_value = np.abs( xt[np.argmax( np.abs(xt) )] ) #    count_num = []    x = []    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)        x.append( (high+low)/2 )    count_num = count_num/np.sum(count_num)    plt.bar(np.arange(len(count_num)), height=count_num)  # 绘制柱状图    plt.show()

plt_amp_prob_density_fun(data=x, n=100)

可知概率密度函数顶部出现一个“碗状”。这是由于正弦函数的概率密度函数是向上抛物线状，而随机信号的概率密度函数是正态分布曲线。

4.3 计算时域特征

首先定义计算时域特征的函数

def get_time_domain_features(data):    x = data    #------有量纲指标------#    X_mean = np.mean(x)   # 1.平均值    X_std = np.std(x)     # 2.方差    x = x - X_mean        # 零均值化（去直流分量）    X_rms2 = np.sum(x**2)/len(x)  # 3.均方值    X_rms = np.sqrt(X_rms2)   # 4.均方根值(有效值)    X_max = np.max(x)     # 5.最大值    X_min = np.min(x)     # 6.最小值    X_p = max(abs(X_max), abs(X_min)) # 7.峰值    X_pp = X_max - X_min  # 8.峰峰值    X_avg = np.mean(np.abs(x))  # 9.平均绝对幅值    X_r = np.mean( np.sqrt(np.abs(x)) )**2  # 10.方根幅值    X_alpha = np.mean( np.power(x, 3) )     # 11.偏度(歪度)    X_beta = np.mean( np.power(x, 4) )   # 12.峭度
    #-------无量纲指标-------#    X_wf = X_rms/X_avg  # 13.波形指标    X_cf = X_pp/X_rms   # 14.峰值指标    X_if = X_pp/X_avg   # 15.脉冲指标    X_clf = X_pp / X_r  # 16 裕度指标    X_pf = X_alpha / X_rms ** 3  # 17.偏度指标    X_kf = X_beta/X_rms ** 4     # 18.峭度指标
    time_domain_features_list = [X_mean, X_std, X_rms2, X_rms, X_max, X_min, X_p, X_pp, X_avg, X_r, X_alpha, X_beta, X_wf, X_cf, X_if, X_clf, X_pf, X_kf]    time_domain_names_list = ['平均值', '方差', '均方值', '均方根值', '最大值', '最小值', '峰值', '峰峰值', '平均绝对幅值',                               '方根幅值', '偏度', '峭度', '波形指标', '峰值指标', '脉冲指标', '裕度指标', '偏度指标', '峭度指标']    for i in range(len(time_domain_names_list)):        print(time_domain_names_list[i],':',time_domain_features_list[i])    return time_domain_features_list, time_domain_names_list

然后调用函数

time_domain_features_list, time_domain_names_list = get_time_domain_features(data=x)

得到结果

>>>输出结果平均值 : 0.506101294634079方差 : 2.1422782642540135均方值 : 4.589356161495189均方根值 : 2.1422782642540135最大值 : 3.4797421092734804最小值 : -3.4874981121547806峰值 : 3.4874981121547806峰峰值 : 6.967240221428261平均绝对幅值 : 1.920591101461503方根幅值 : 1.752480189593869偏度 : -0.07538076391033736峭度 : 32.72574911884685波形指标 : 1.1154265281265827峰值指标 : 3.2522573456881902脉冲指标 : 3.627654119675153裕度指标 : 3.9756456379931424偏度指标 : -0.007667131112382436峭度指标 : 1.553767635488039

5、CWRU数据实战

5.1 内圈故障

定义CWRU数据读取函数

def data_acquision(FilePath):    """    fun: 从cwru mat文件读取加速度数据    param file_path: mat文件绝对路径    return accl_data: 加速度数据，array类型    """    data = scio.loadmat(file_path)  # 加载mat数据    data_key_list = list(data.keys())  # mat文件为字典类型，获取字典所有的键并转换为list类型    accl_key = data_key_list[3]  # 获取'X108_DE_time'    accl_data = data[accl_key].flatten()  # 获取'X108_DE_time'所对应的值，即为振动加速度信号,并将二维数组展成一维数组    return accl_data

使用CWRU轴承数据进行分析，选取了4个mat文件，包括内圈故障、外圈故障、滚动体故障和正常状态。

文件名解释（以“1730_12k_0.007-InnerRace”为例）：

1730：转速，单位rpm

12k：采样频率，Hz

0.007：故障大小，单位inch，1inch=25.4mm

InnerRace：代表为内圈故障

文件名解释（以“1730_12k_0.007-InnerRace”为例）：

1730：转速，单位rpm

12k：采样频率，Hz

0.007：故障大小，单位inch，1inch=25.4mm

InnerRace：代表为内圈故障

更多CWRU数据集介绍见

绘制内圈故障的时域图和幅值概率密度函数

file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第6篇-时频域指标/1730_12k_0.007-InnerRace.mat'data = data_acquision(FilePath=file_path)print(data)plt.figure(figsize=(12,2))plt.plot(data)plt.show()plt_amp_prob_density_fun(data=data, n=100)

>>>输出结果[ 0.22269856  0.09323776 -0.14651649 ... -0.36125573  0.31138814  0.17055689]

计算内圈故障的各时域统计特征大小

get_time_domain_features(data=data)

>>>输出结果平均值 : 0.004718397126522763方差 : 0.3135712339852505均方值 : 0.09832691878303272均方根值 : 0.3135712339852505最大值 : 1.6667390879034174最小值 : -1.5402176785636486峰值 : 1.6667390879034174峰峰值 : 3.2069567664670657平均绝对幅值 : 0.22413171955956412方根幅值 : 0.1788814956202537偏度 : -0.00040761964834175015峭度 : 0.05115487013606729波形指标 : 1.3990488923274307峰值指标 : 10.227203323816086脉冲指标 : 14.308357481792312裕度指标 : 17.92782845060225偏度指标 : -0.01322045690393005峭度指标 : 5.291053175312336

可发现其峭度指标大于3，幅值概率密度函数偏瘦尖。

5.2 正常状态

绘制正常状态的时域图和幅值概率密度函数图。

file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第6篇-时频域指标/1730_48k_Normal.mat'data = data_acquision(FilePath=file_path)print(data)plt.figure(figsize=(12,2))plt.plot(data)plt.show()plt_amp_prob_density_fun(data=data, n=100)

>>>输出结果[ 0.01460308  0.05444862  0.10764554 ... -0.02357354  0.00521538  0.04777292]

计算正常状态轴承的各时域统计特征大小

get_time_domain_features(data=data)

>>>输出结果平均值 : 0.01245851422293584方差 : 0.06469520385023537均方值 : 0.004185469401223511均方根值 : 0.06469520385023537最大值 : 0.2712584088539872最小值 : -0.3189145142229358峰值 : 0.3189145142229358峰峰值 : 0.590172923076923平均绝对幅值 : 0.05170725167996645方根幅值 : 0.04384435047720211偏度 : -3.4532433523248354e-05峭度 : 5.180413666679372e-05波形指标 : 1.2511824115243202峰值指标 : 9.12235974158347脉冲指标 : 11.41373606026678裕度指标 : 13.460637839390442偏度指标 : -0.1275295794513686峭度指标 : 2.9571686803135426

可发现其峭度指标在3左右，幅值概率密度函数与正态分布差不多。

5.3 外圈故障

绘制外圈故障的时域图和幅值概率密度函数图


file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第6篇-时频域指标/1730_12k_0.007-OuterRace3.mat'data = data_acquision(FilePath=file_path)print(data)plt.figure(figsize=(12,2))plt.plot(data)plt.show()plt_amp_prob_density_fun(data=data, n=100)

>>>输出结果[ 0.10436457  0.62537525  0.12304461 ...  0.88689581  0.84588094 -0.5636499 ]

计算外圈故障的各幅域特征参数大小

get_time_domain_features(data=data)

>>>输出结果平均值 : 0.006673843393593249方差 : 0.7975703939021197均方值 : 0.6361185332291823均方根值 : 0.7975703939021197最大值 : 3.612380847225169最小值 : -3.322787017046288峰值 : 3.612380847225169峰峰值 : 6.935167864271457平均绝对幅值 : 0.6080145950947324方根幅值 : 0.5062975339786826偏度 : 0.04229527425007543峭度 : 1.6551101479225734波形指标 : 1.3117619220601988峰值指标 : 8.695367728409641脉冲指标 : 11.406252284438855裕度指标 : 13.697810869771013偏度指标 : 0.08336519532029546峭度指标 : 4.090258951031924

可发现其峭度指标大于3，幅值概率密度函数顶部出现凹坑。

5.4 滚动体故障

绘制滚动体故障的时域图和幅值概率密度函数图

file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第6篇-时频域指标/1730_12k_0.014-Ball.mat'data = data_acquision(FilePath=file_path)print(data)plt.figure(figsize=(12,2))plt.plot(data)plt.show()plt_amp_prob_density_fun(data=data, n=100)

>>>输出结果[ 0.1054204  -0.10736962 -0.16340974 ... -0.5821675   0.07488259  0.56202555]

计算滚动体故障的各时域统计特征大小

get_time_domain_features(data=data)

>>>输出结果平均值 : 0.004535467952873973方差 : 0.13361739428910718均方值 : 0.017853608056610733均方根值 : 0.13361739428910718最大值 : 1.9846451308495212最小值 : -1.7174139110666464峰值 : 1.9846451308495212峰峰值 : 3.7020590419161676平均绝对幅值 : 0.09177927629408723方根幅值 : 0.07371838887961654偏度 : 0.0003906736827228592峭度 : 0.004736237523606714波形指标 : 1.4558558280734157峰值指标 : 27.70641548289771脉冲指标 : 40.33654645580015裕度指标 : 50.21893584735956偏度指标 : 0.16376653561183255峭度指标 : 14.858722825401527

可发现其峭度指标远远大于3，幅值概率密度函数也十分瘦尖。

6、IMF加速寿命实验轴承数据实战

6.1 IMF数据介绍

IMF是一个轴承全寿命周期数据，即在试验台上从正常跑到坏。IMF数据是每隔10分钟采集一次。

文件名是采集时间

如“2004.02.12.10.32”为2004年2月12日10:32采集。

本次数据选取的是一个最终跑出外圈故障的数据。

6.2 单个数据尝试

定义IMF数据读取函数

def IMF_data_acquision(FilePath, n):    """    fun: 从IMF文件读取加速度数据    param file_path: mat文件绝对路径    return accl_data: 加速度数据，array类型    """    data = np.loadtxt(fname=file_path)    accl_data = data[:, n]    return accl_data

file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第6篇-时频域指标/IMF_2nd_test_OF1/2004.02.12.10.32.39'data = IMF_data_acquision(FilePath=file_path, n=0)

绘制时域图和幅值概率密度函数图、计算时域统计特征大小

print(data)plt.figure(figsize=(12,2))plt.plot(data)plt.show()plt_amp_prob_density_fun(data=data, n=100)time_domain_features_list, time_domain_names_list = get_time_domain_features(data=data)

>>>输出结果[-0.049 -0.042  0.015 ... -0.012 -0.012  0.02 ]

>>>输出结果平均值 : -0.010195996093750001方差 : 0.07347493104305192均方值 : 0.005398565491781235均方根值 : 0.07347493104305192最大值 : 0.46419599609375最小值 : -0.37580400390625峰值 : 0.46419599609375峰峰值 : 0.8400000000000001平均绝对幅值 : 0.057688910923004155方根幅值 : 0.04853102183425118偏度 : 3.331677118397008e-05峭度 : 0.00010575850683660837波形指标 : 1.2736404599684148峰值指标 : 11.432470749891692脉冲指标 : 14.5608573044675裕度指标 : 17.308516661134114偏度指标 : 0.08399343541253229峭度指标 : 3.628762642644256

6.2 整个寿命周期数据分析

接下来遍历读取全寿命周期数据。并获取每次采集数据的时域统计特征大小。将所有文件的时域统计特征保存为二维list类型。

base_dir = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第6篇-时频域指标/IMF_2nd_test_OF1'file_name_list = os.listdir(base_dir)all_time_domain_features_list = []for file_name in file_name_list:    file_path = os.path.join(base_dir, file_name)    data = IMF_data_acquision(FilePath=file_path, n=0)    time_domain_features_list, time_domain_names_list = get_time_domain_features(data=data)    all_time_domain_features_list.append(time_domain_features_list)

>>>输出结果平均值 : -0.010195996093750001方差 : 0.07347493104305192均方值 : 0.005398565491781235均方根值 : 0.07347493104305192最大值 : 0.46419599609375最小值 : -0.37580400390625峰值 : 0.46419599609375峰峰值 : 0.8400000000000001...峰值指标 : 7.003117715378994脉冲指标 : 7.147508219574274裕度指标 : 7.207950142678743偏度指标 : 0.31703227680043844峭度指标 : 1.3902259645328505

将二维list转为二维array

all_time_domain_features_arr = np.array(all_time_domain_features_list)all_time_domain_features_arr

>>>结果输出array([[-1.01959961e-02,  7.34749310e-02,  5.39856549e-03, ...,         1.73085167e+01,  8.39934354e-02,  3.62876264e+00],       [-2.58500977e-03,  7.53381213e-02,  5.67583252e-03, ...,         1.52497399e+01,  5.21419073e-02,  3.64829108e+00],       [-2.48398437e-03,  7.61892928e-02,  5.80480834e-03, ...,         1.77537992e+01,  3.28079981e-02,  3.51347531e+00],       ...,       [-1.70307617e-03,  4.83832398e-01,  2.34093789e-01, ...,         2.50246000e+01, -3.77095385e-01,  7.89175523e+00],       [ 1.85683594e-03,  9.87327986e-04,  9.74816551e-07, ...,         1.46340168e+01,  5.79698251e-01,  6.63751264e+00],       [-1.16196289e-03,  9.99554810e-04,  9.99109819e-07, ...,         7.20795014e+00,  3.17032277e-01,  1.39022596e+00]])```

绘制整个周期的每个时域特征参数大小变化

for i in range(len(time_domain_features_list)):    plt.figure(figsize=(12,3))    plt.plot(all_time_domain_features_arr[:,i])    plt.title(time_domain_names_list[i])    plt.show()