机械故障诊断信号的幅域分析 - 幅值概率密度函数 |在CWRU轴承数据上实战
很多初学者刚接触故障诊断可能觉得很简单，套用深度学习模型进行训练，分类准确率达到99%即可。
在写论文时，这样的确没问题。但是在工作或者在解决实际问题时，几乎是用不上。
最近在读时献江老师等人的《机械故障诊断及典型案例解析》这本书时，发现其内容真的很丰富，故障诊断方法讲解的详细全面，且有实际丰富的案例，更加容易消化理解各方法特点及其应用场景。
信号波形是某种物理量随时间变化的关系。机械振动信号处理的基本方法有幅域分析、时域分析和频域分析。
在幅值上的各种统计处理通常称为幅域分析。仅对波形的幅值用简单的方法进行统计分析， 如计算波形的最大值、平均值和有效值，或研究时域波形幅值的概率分布形式等，
信号在时间域内的变换或统计分析称为时域分析。频域分析是确定信号的频率结构，即信号中包含哪些频率成分，分析的结果是以频率为自变量的各种物理量的谱线或曲线。
不同的分析方法是从不同的角度观察、分析信号，使信号处理的结果更加利于故障分析与诊断。
本节将介绍幅域分析方法中的幅值概率密度函数，包括其数学公式定义及CWRU轴承数据实战案例。
1. 幅值概率密度函数定义
2. 基于CWRU数据代码实战准备和数据介绍
2.1 导入包
2.2 定义CWRU数据读取函数
3. 轴承内圈故障幅值概率密度函数案例分析
3.1 时域图绘制
3.2 编程思路分析
4. 封装成一个plt_amp_prob_density_fun()函数
4.1 轴承滚动体轴承故障幅值概率密度函数案例分析
4.2 轴承外圈轴承故障幅值概率密度函数案例分析
4.3 正常状态轴承故障幅值概率密度函数案例分析
5. 总结
6. 课后思考
1、幅值概率密度函数定义
随机信号的幅值概率密度函数表示信号的幅值落在某一个指定区间内的概率，幅值概率密度函数提供了随机信号沿幅值域分布的信息，是随机信号的主要统计特性参数之一。在图1中，
当
图1 时域波形及幅值概率密度函数
典型信号的时域波形和幅值概率密度函数如图2所示。根据随机过程理论，随机信号的幅值概率密度函数符合正态分布规律，而确定性信号如简谐信号的幅值概率密度函数则呈盆形曲线，如图2(a)所示。一般故障信号多是随机信号和简谐信号的混合体，所以当信号幅值概率密度函数的正态分布曲线上端出现盆型漏斗时［图2(b)]，往往预示着系统存在故障征兆。
图2典型信号的时域波形和幅值概率密度函数
可以采用幅值概率密度函数直接进行设备状态的检测与诊断。图3是新旧机床变速箱的噪声分布规律，可见，新旧两个变速箱的分布规律有着明显的差异。这是因为机床齿轮箱中的零件由于磨损等原因使得配合间隙增大，在机床噪声中会出现大幅值周期性冲击成分，使得噪声信号的方差增加，分散度加大，甚至使曲线的顶部变平或出现局部的凹形。
图3 新旧机床变速箱噪声分布规律
2、基于CWRU数据代码实战准备和数据介绍
2.1 导入包
import scipy.io as scioimport numpy as npimport matplotlib.pyplot as pltfrom scipy import signal, fftpack, stats
from matplotlib import rcParamsconfig = {"font.family": 'serif', # 衬线字体"font.size": 10, # 相当于小四大小"font.serif": ['SimSun'], # 宋体"mathtext.fontset": 'stix', # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大'axes.unicode_minus': False # 处理负号，即-号}rcParams.update(config)
2.2 定义CWRU数据读取函数
使用CWRU轴承数据进行分析，选取了4个mat文件，包括内圈故障、外圈故障、滚动体故障和正常状态。
文件名解释（以“1730_12k_0.007-InnerRace”为例）：
1730：转速，单位rpm
12k：采样频率，Hz
0.007：故障大小，单位inch，1inch=25.4mm
InnerRace：代表为内圈故障
更多CWRU数据集介绍见
下面定义一个mat文件数据读取函数
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
3、内圈故障幅值概率密度函数案例分析
3.1 时域图绘制
file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第4篇-包络谱/1730_12k_0.007-InnerRace.mat'xt = data_acquision(file_path) #从mat文件获取振动加速度数据##--------绘制时域图-------##plt.figure(figsize=(12,3))plt.plot(xt)plt.show()print(xt)
[ 0.22269856 0.09323776 -0.14651649 ... -0.36125573 0.311388140.17055689]
面对这10s的时域信号，如何实现获取绝对幅值概率密度函数呢？
3.2 编程思路分析
step1：确定幅值绝对值最大值max_value
step2：在幅值区间[-max_value, max_value]内等分划分成n个小区间，区间长度为interval_len
step3：在第i个幅值小区间[-max_value+i * interval_len, max_value+(i+1) * interval_len]内，统计落入该区间的幅值个数count_num
step5：将n个区间统计的个数count_num形成一个列表count_num_list
step6：绘制柱状图，即得到幅值概率密度图
def interval_num_count(data, low, high):'''fun: 统计一维数据data落入某一个区间[low, high]内的数量param low: 区间下限param high: 区间上限return count_num: 落入某一个区间[low, high]内的数量'''count_num = 0for i in range(len(data)):if data[i]>low and data[i]<high:count_num += 1return count_num
n = 10 # 分段数量xt = xt - np.mean(xt) # 去直流分量（也叫零均值化处理）max_value = np.abs( xt[np.argmax( np.abs(xt) )] ) # 确定幅值绝对值最大值max_valuecount_num_list = []for i in range(n):interval_len = max_value*2/n # 在幅值区间[-max_value, max_value]内等分划分成n个小区间，区间长度为interval_lenlow = -max_value + i*interval_len # 第i个幅值小区间[-max_value+i * interval_len, max_value+(i+1) * interval_len]的下限high = -max_value + (i+1)*interval_len # 第i个幅值小区间[-max_value+i * interval_len, max_value+(i+1) * interval_len]的上限count_num = interval_num_count(data=xt, low=low, high=high) # 统计落入该区间的幅值个数count_numcount_num_list.append(count_num) # 将n个区间统计的个数count_num形成一个列表count_num_listplt.bar(x=range( len(count_num_list) ), height=count_num_list) # 绘制柱状图，即得到幅值概率密度图
是不是已成雏形，有正态分布形状了。只是原始区间分成了10份，看不太出来
下面分成100份
n = 100xt = xt - np.mean(xt) # 去直流分量（也叫零均值化处理）max_value = np.abs( xt[np.argmax( np.abs(xt) )] ) # 确定幅值绝对值最大值max_valuecount_num_list = []for i in range(n):interval_len = max_value*2/n # 在幅值区间[-max_value, max_value]内等分划分成n个小区间，区间长度为interval_lenlow = -max_value + i*interval_len # 第i个幅值小区间[-max_value+i * interval_len, max_value+(i+1) * interval_len]的下限high = -max_value + (i+1)*interval_len # 第i个幅值小区间[-max_value+i * interval_len, max_value+(i+1) * interval_len]的上限count_num = interval_num_count(data=xt, low=low, high=high) # 统计落入该区间的幅值个数count_numcount_num_list.append(count_num) # 将n个区间统计的个数count_num形成一个列表count_num_listplt.bar(x=range( len(count_num_list) ), height=count_num_list) # 绘制柱状图，即得到幅值概率密度图
这下有正态分布那味了。不得不说这个形状还挺优美的，看着挺舒服。
4、封装成一个plt_amp_prob_density_fun()函数
def plt_amp_prob_density_fun(data, n):'''fun: 绘制幅值概率密度函数param data: 输入数据，1维arrayparam n: 分割成段数的数量return: 绘制幅值概率密度函数'''max_value = np.abs( xt[np.argmax( np.abs(xt) )] ) #count_num = []for i in range(n):interval = max_value*2/n # 区间长度为interval_lenlow = -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()
4.1 滚动体故障轴承幅值概率密度函数案例分析
此时的正态分布形状比较瘦小，峭度K值大于3。K值大于3时，说明信号中冲击成分幅值增大。
4.2 外圈故障轴承幅值概率密度函数案例分析
file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第4篇-包络谱/1730_12k_0.007-OuterRace3.mat'xt = data_acquision(file_path)plt.figure(figsize=(12,3))plt.plot(xt)plt.show()n = 100 # 设定分成份数plt_amp_prob_density_fun(data=xt, n=n)
在顶部出现了第1节介绍时说的盆型漏斗。
4.3 正常状态轴承幅值概率密度函数案例分析
file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第4篇-包络谱/1730_48k_Normal.mat'xt = data_acquision(file_path)plt.figure(figsize=(12,3))plt.plot(xt)plt.show()n = 100 # 设定分成份数plt_amp_prob_density_fun(data=xt, n=n)
看着与标准正态分布挺像的。
5、总结
6、课后思考
该方法的不足：
采集的样本点数数量与幅值概率密度函数成正比。因此对采样样本数量有有一定要求。
除了这个不足，你还发现哪些不足呢？
该程序存在两个不足：
第一个是纵坐标没有归一化，因为概率之和应该为1。
第二个是横坐标应该[-max_value, max_value]，也就是概率密度函数图应该分布在纵坐标两边。
发动你的小脑袋来完善这个程序吧。
注明：
本文内容摘抄自《机械故障诊断及典型案例解析》，仅供学习交流，若有侵权，烦请联系处理！
代码链接：https://github.com/HappyBoy-cmd/fault_diagnosis_signal_processing
参考文献
[1] 《机械故障诊断及典型案例解析，第2版，时献江等，2020年第一次印刷》
