机械故障诊断及工业工程故障诊断若干例子（第二篇）

MATLAB环境下基于平方包络谱的滚动轴承故障诊断
算法程序运行环境为MATLAB R2018a，执行基于平方包络谱的滚动轴承故障诊断，也可用于金融时间序列，地震信号，机械振动信号，语音信号，声信号等一维时间序列信号。

Python环境下基于麦克风信号与随机森林的机器轴承运行状态识别

以往的轴承故障诊断都是围绕内外圈、滚动体等进行的，而该算法程序使用麦克风阵列信号和随机森林模型进行机器轴承运行状态识别，轴承运行状态包括正常状态（good bearing），外圈单点损伤状态（bad bearing），低速运行状态（slow speed bearing），过度润滑状态（over lubrication），轴承润滑不足状态（under lubrication），轴承烧蚀故障状态（corona），轴承轻微变形（上拱，arching）等状态，看一下所使用的麦克风音频文件。

本项目要用到librosa模块，librosa是一个用于音乐和音频分析的python包，可以pip install librosa

首先导入相关模块

import librosa
import librosa.display
from scipy.io import wavfile as wav 
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

from sklearn.ensemble import RandomForestClassifier
#pip install xgboost
#from xgboost import XGBClassifier

导入音频样本并部分可视化

audio_file = "over lubrication.wav" #轴承过度润滑状态
data, sample_rate = librosa.load(audio_file)
plt.figure(figsize=(20,5)) # set size of visualization
librosa.display.waveshow(data, sr= sample_rate)
plt.xlabel('Time (samples)')
plt.ylabel('Amplitude')

audio_file = "under lubrication.wav" #轴承润滑不足
data, sample_rate = librosa.load(audio_file)
plt.figure(figsize=(20,5)) # set size of visualization
librosa.display.waveshow(data, sr= sample_rate)
plt.xlabel('Time (samples)')
plt.ylabel('Amplitude')

audio_file = "corona.wav" #轴承烧蚀故障
data, sample_rate = librosa.load(audio_file)
plt.figure(figsize=(20,5)) # set size of visualization
librosa.display.waveshow(data, sr= sample_rate)
plt.xlabel('Time (samples)')
plt.ylabel('Amplitude')

看一下某个样本音频的直方图

# 采样频率
print('The sampling frequency is: ' + str(sample_rate) + ' Hz')
# 音频样本的直方图
plt.figure(figsize=(20,5))
plt.hist(data, bins=100)
plt.xlabel('Amplitude')
plt.ylabel('Frequency')

大多数值接近0。

特征提取

对一个样本音频文件进行特征提取(梅尔频率倒谱系数MFCC)，MFCC一般是用于表示声音信号的特征集

mfccs = librosa.feature.mfcc(y = data, sr = sample_rate, n_mfcc=40)
mfccs.shape
#有2371 组不同的 40 个值的原因是因为音频样本的长度
#查看通过 mfcc 生成的特征
mfccs[0]

对所有音频文件进行特征提取

def feature_extractor(file):
    audio, s_r = librosa.load(file)
    mfcc_feat = librosa.feature.mfcc(y=audio, sr = s_r, n_mfcc=40)
    mfcc_scaled_feat = np.mean(mfcc_feat.T, axis=0)
    return mfcc_scaled_feat
 #输出标签
 labels = ['arching', 'bad bearing', 'good bearing', 'corona', 'over lubrication', 'under lubrication', 'tracking', 'slow speed bearing']
extracted_feat = []
for i in labels:
    file_path =  '' + i + ".wav"
    data = feature_extractor(file_path)
    extracted_feat.append(data)

看一下特征长度

len(extracted_feat)

显示提取特征的样本及其相应的标签

print(str(str(labels[0]) + " : " + str(extracted_feat[0])))

创建pandas数据框架，数据框架的的特征为列，标签为行

data = {'Features':extracted_feat,
        'Class':labels}
df = pd.DataFrame(data)
# 显示数据框架
df = df.convert_dtypes()
df

初始化 RandomForestClassifier分类器

rf = RandomForestClassifier()

创建输入和输出数组 ，输入数组是特征，输出数组是标签

X = df['Features'].tolist()
y = df['Class'].tolist()

训练模型

rf.fit(X, y)

准确率

rf.score(X,y)

准确率为100%

MATLAB环境下基于双树复小波变换的轴承故障诊断

算法程序运行环境为MATLAB R2021b，执行基于双树复小波变换的轴承故障诊断，也可用于金融时间序列，地震信号，机械振动信号，语音信号，声信号等一维时间序列信号。

双树复小波变换DTCWT 为两个独立的两通道滤波器组，在实际应用用，不能随意选择两棵树中使用的尺度小波滤波器。第一棵树 {h0 ,h1 } 的低通（尺度）和高通（小波）滤波器生成一个尺度函数和小波，另一棵树是由第一棵树的尺度函数的近似希尔伯特变换以及相应的高通滤波器生成的小波函数组成，记作{g0，g1 }。因此由两棵树形成的复值尺度函数和小波函数是近似解析的，DTCWT的冗余度明显小于未抽取的DWT 的冗余度。

双树复小波变换DTCWT基本理论

双树复小波变换DTCWT采用二叉树结构的两路滤波器组进行信号的分解和重构，第一棵树生成实部，第二棵生成虚部，合理设计实、虚部树低通滤波器，满足半采样延迟条件，具有近似平移不变性。两树的滤波器采样频率相同，但是它们之间的延迟恰好是一个采样间隔，这样虚部树中第１层的二抽取恰好采到实部树中二抽取所丢掉的采样值，在获得了复小波变换的平移不变性的同时避免了大量的计算并且具有容易实现的优势。下图为３层双树复小波的分解和重构过程。

虚线上方实部树变换的小波系数和尺度系数可由式（2）和（3）计算

相关的参考文献如下：

[1]Huang Tongyuan,Xu Jia,Yang Yuling,Han Baoru. Robust Zero-Watermarking Algorithm for Medical Images Using Double-Tree Complex Wavelet Transform and Hessenberg Decomposition[J]. Mathematics,2022,10(7).

[2]Zhou Yilu,Fu Xiaojin. Image Denoising Based on Dual-tree Complex Wavelet Transform and Convolutional Neural Network[J]. Journal of Physics: Conference Series,2021,1995(1).

[3]Lei Wang,Zhiwen Liu,Hongrui Cao,Xin Zhang. Subband averaging kurtogram with dual-tree complex wavelet packet transform for rotating machinery fault diagnosis[J]. Mechanical Systems and Signal Processing,2020,142(C).

下面开始进行轴承故障诊断

首先以70Hz轴承外圈故障为例，故障特征频率因子为BPFI=6.587，采样频率为10240HZ,首先看一下原始信号的包络谱

可见包络谱的故障特征频率虽然出现了，但是幅值相对较低。

下面进行DTCWT分解，分解层数为5

看一下实部树所对应的波形

实部树所对应的包络谱

重点看一下Level 2的包络谱

故障特征频率的相对幅值较高

然后看一下虚部树所对应的分解波形

虚部树所对应的包络谱

然后看一下50Hz的情形，包络谱如下

实部树分解所得到的包络谱如下

详细看一下Level 1的包络谱

可见轴承故障特征频率的相对幅值较高，虚部树与此类似

面包多代码

https://mbd.pub/o/GeBENHAGEN