Python轴承故障诊断 (二)连续小波变换CWT

import numpy as np

import matplotlib.pyplot as plt

import pywt


# 生成三个不同频率成分的信号 ?3000个点

fs = 1000 ?# 采样率

time = np.linspace(0, 1, fs, endpoint=False) ?# 时间

# 第一个频率成分

signal1 = np.sin(2 * np.pi * 30 * time)

# 第二个频率成分

signal2 = np.sin(2 * np.pi * 60 * time)

# 第三个频率成分

signal3 = np.sin(2 * np.pi * 120 * time)

# 合并三个信号

signal = np.concatenate((signal1, signal2, signal3))


# 连续小波变换参数

# 采样频率

sampling_rate = 3000

# 尺度长度

totalscal = 128 ? ?

# 小波基函数

wavename = 'morl'

# 小波函数中心频率

fc = pywt.central_frequency(wavename)

# 常数c

cparam = 2 * fc * totalscal ?

# 尺度序列

scales = cparam / np.arange(totalscal, 0, -1)

?# 进行CWT连续小波变换

coefficients, frequencies = pywt.cwt(signal, scales, wavename, 1.0/1000)

# 小波系数矩阵绝对值

amp = abs(coefficients)


# 根据采样频率 sampling_period 生成时间轴 t

t = np.linspace(0, 1.0/sampling_rate, sampling_rate, endpoint=False)

# 绘制时频图谱

plt.figure(figsize=(20,10))

plt.subplot(2,1,1)

plt.plot(signal)

plt.title('30Hz和60Hz和120Hz的分段波形')

plt.subplot(2,1,2)

plt.contourf(t, frequencies, amp, cmap='jet')

plt.title('对应时频图')

plt.show()

参数解释

totalscal：表示尺度长度，表征频率的参数，选择不同的尺度可以捕捉信号不同尺度上的特征；
wavename：表示小波函数的类型，不同的小波适用于不同类型的信号；
fc：表示小波函数的中心频率，小波函数在频率域中有一个中心频率；这个中心频率是与小波函数形状和性质有关的一个重要参数；
scales：表示尺度序列，CWT本质上是将你的信号与不同尺度的小波进行相关，scales 参数确定尺度范围；

coefficients：表示信号变换后的小波系数；
frequencies ：表示对应的频率信息

2.2 参数介绍和选择策略

2.2.1 尺度长度：

在连续小波变换（CWT）中，尺度参数是一个关键的选择，因为它决定了小波函数的宽度，从而影响了频率分辨率。尺度与频率成反比，尺度反映了分析的频率范围，尺度越小，小波函数衰减越快，频率越高；尺度越大，小波函数衰减越慢，频率越低[1]。

选择小波尺度的一般原则是：

高频特征：如果关注信号的高频特征，应该选择较小的尺度；
低频特征：如果关注信号的低频特征，应该选择较大的尺度；

覆盖感兴趣的频率范围：尺度参数的选择应该使小波函数能够覆盖感兴趣的频率范围；如果期望信号有很高的频率变化，可能需要选择较小的尺度；
频率分辨率： 较小的尺度提供更好的频率分辨率，因为小波函数较窄，可以更精细地定位频率；但是，这也意味着在时间上的分辨率较差，因此，需要权衡时间分辨率和频率分辨率；
信号持续时间：尺度参数的选择还应考虑信号的持续时间，如果信号是短暂的，可能需要较小的尺度；
尺度间隔： 在尺度参数上选择合适的间隔，以确保在整个频率范围内进行了适当的采样；这取决于具体的应用和信号特性。

2.2.2 小波函数（wavelet）：

小波函数（wavelet）的选择也连续小波变换中的一个重要参数，它决定了小波基函数的形状，不同的小波函数适用于不同类型的信号和应用。

打印 Python pywt 包中的所有小波函数类型：

import pywt

# 获取小波函数列表

wavelets = pywt.wavelist()

# 打印小波函数列表

# 按照12行，10列的形式打印数据

num_rows = 12

num_columns = 10

for i in range(0, len(wavelets), num_columns):

? ? row = wavelets[i:i + num_columns]

? ? print(row)

介绍常用的小波基函数：

'morl' ：Morlet小波是一种复杂的小波函数，它在频率域和时域都有较好的局部化性质；Morlet小波通常用于处理时频局部化要求较高的信号，比如处理振动信号或某些生物医学信号；
'cmor': Complex Morlet wavelet，复数 Morlet 小波，是 Morlet 小波的一个变种；
'cgau' ：Complex Gaussian wavelet，复数高斯小波，用于近似高斯信号；
'db1'：Daubechies小波是离散小波变换（Discrete Wavelet Transform, DWT）中常用的小波函数。Daubechies小波是紧支撑的小波，适用于处理有限长度的信号。
'haar'：Haar小波是最简单的小波函数之一，适用于对信号进行基本的低通和高通分解；
'mexh'：Mexican Hat小波，也称为Ricker小波，适用于处理具有尖峰或波包特性的信号；

'bior1.1'：Bior小波是一类双正交小波，其特点是具有对称和非对称两组滤波器，尤其适用于一些信号的多分辨率分析，如图像处理；
'sym2'：Symlet小波（Symmetric Wavelets），是一类对称的小波函数，它在某些方面类似于Daubechies小波，但是Symlet小波在设计上更加灵活。Symlet小波也是一种紧支撑小波，适用于有限长度的信号处理。

小波函数的选择通常取决于处理的信号类型以及分析的目标。在实际应用中，可以尝试不同的小波函数，观察它们在信号上的效果，然后根据实验结果选择最适合的小波函数。在选择小波函数时，也要考虑小波函数的性质，如平滑性、局部化等。

2.3 凯斯西储大学轴承数据的加载

选择正常信号和?0.021英寸内圈、滚珠、外圈故障信号数据来做对比

第一步，导入包，读取数据

import numpy as np

from scipy.io import loadmat

import numpy as np

from scipy.signal import stft

import matplotlib.pyplot as plt

import matplotlib

matplotlib.rc("font", family='Microsoft YaHei')


# 读取MAT文件 ?

data1 = loadmat('0_0.mat') ?# 正常信号

data2 = loadmat('21_1.mat') # 0.021英寸 内圈

data3 = loadmat('21_2.mat') # 0.021英寸 滚珠

data4 = loadmat('21_3.mat') # 0.021英寸 外圈

# 注意，读取出来的data是字典格式，可以通过函数type(data)查看。

第二步，数据集中统一读取?驱动端加速度数据，取一个长度为1024的信号进行后续观察和实验

# DE - drive end accelerometer data 驱动端加速度数据

data_list1 = data1['X097_DE_time'].reshape(-1)

data_list2 = data2['X209_DE_time'].reshape(-1) ?

data_list3 = data3['X222_DE_time'].reshape(-1)

data_list4 = data4['X234_DE_time'].reshape(-1)

# 划窗取值（大多数窗口大小为1024）

data_list1 = data_list1[0:1024]

data_list2 = data_list2[0:1024]

data_list3 = data_list3[0:1024]

data_list4 = data_list4[0:1024]

第三步，进行数据可视化

plt.figure(figsize=(20,10))

plt.subplot(2,2,1)

plt.plot(data_list1)

plt.title('正常')

plt.subplot(2,2,2)

plt.plot(data_list2)

plt.title('内圈')

plt.subplot(2,2,3)

plt.plot(data_list3)

plt.title('滚珠')

plt.subplot(2,2,4)

plt.plot(data_list4)

plt.title('外圈')

plt.show()

2.4 CWT与参数选择对比

本实验以某轴承故障诊断论文推荐'cgau8'小波，以及'morl'小波、'cmor1-1'小波、'cmor1.5-2'小波为实验做对比，尺度先设定为128，来对比不同小波函数的影响。

2.4.1 基于尺度为128，选择内圈数据比较 CWT 的不同小波函数

import numpy as np

import matplotlib.pyplot as plt

import pywt

import pandas as pd


# 连续小波变换参数

# 采样频率

sampling_rate = 1024

# 尺度长度

totalscal = 128 ?


wavename1 = 'cgau8'

fc1 = pywt.central_frequency(wavename1)

cparam1 = 2 * fc1 * totalscal ?

scales1 = cparam1 / np.arange(totalscal, 0, -1)


wavename2 = 'morl' #

fc2 = pywt.central_frequency(wavename2)

cparam2 = 2 * fc2 * totalscal ?

scales2 = cparam2 / np.arange(totalscal, 0, -1)


wavename3 = "cmor1-1"

fc3 = pywt.central_frequency(wavename3)

cparam3 = 2 * fc3 * totalscal ?

scales3 = cparam3 / np.arange(totalscal, 0, -1)


wavename4 = 'cmor1.5-2'

fc4 = pywt.central_frequency(wavename4)

cparam4 = 2 * fc4 * totalscal ?

scales4 = cparam4 / np.arange(totalscal, 0, -1)


# 进行连续小波变换

coefficients1, frequencies1 = pywt.cwt(data_list2, scales1, wavename1, sampling_period)

coefficients2, frequencies2 = pywt.cwt(data_list2, scales2, wavename2, sampling_period)

coefficients3, frequencies3 = pywt.cwt(data_list2, scales3, wavename3, sampling_period)

coefficients4, frequencies4 = pywt.cwt(data_list2, scales4, wavename4, sampling_period)


# 小波系数矩阵绝对值

amp1 = abs(coefficients1)

amp2 = abs(coefficients2)

amp3 = abs(coefficients3)

amp4 = abs(coefficients4)


# 根据采样频率 sampling_period 生成时间轴 t

t = np.linspace(0, 1.0/sampling_rate, sampling_rate, endpoint=False)

数据可视化

plt.figure(figsize=(20,10))

plt.subplot(2,2,1)

plt.contourf(t, frequencies1, amp1, cmap='jet')

plt.title('内圈-cgau8')

plt.subplot(2,2,2)

plt.contourf(t, frequencies2, amp2, cmap='jet')

plt.title('内圈-morl')

plt.subplot(2,2,3)

plt.contourf(t, frequencies3, amp3, cmap='jet')

plt.title('内圈-cmor1-1')

plt.subplot(2,2,4)

plt.contourf(t, frequencies4, amp4, cmap='jet')

plt.title('内圈-cmor1.5-2')

plt.show()

对比不同小波函数，对于内圈故障信号来说，'cgau8'小波有着较高的频率分辨率，需要对比其他类型故障数据，进一步观察。

2.4.2 根据正常数据和三种故障数据，对比不同小波函数的辨识度

比较来看，'cgau8'小波和'cmor1.5-2'小波对于不同故障都有着较高的辨识度，综合考虑频率分辨率和时间分辨率，选择'cmor1.5-2'小波来进一步分析。

2.4.3 基于'cmor1.5-2'小波，选择滚珠故障数据比较 CWT 的不同尺度的变化：32、64、128、256；

同时尺度序列scales的常数项cparams均采用：cparam = 2 * fc * totalscal；'cmor1.5-2'小波中1 代表中心频率参数，1.5代表带宽参数；

注意，在连续小波变换中，影响小波系数的是尺度序列scales，但是仍能从上图中看出尺度长度越大，对低频特征关注越多（注意每幅小图的底部低频的区别）；

为进一步探索尺度序列scales的影响，选择尺度长度为128，设置不同常数项cparams来观察对时频变换的影响：

import numpy as np

import matplotlib.pyplot as plt

import pywt

import pandas as pd


# 连续小波变换参数

# 采样频率

sampling_rate = 1024

# 尺度长度

totalscal = 128 ?


wavename1 = 'cmor1.5-2'

fc1 = pywt.central_frequency(wavename1)

cparam1 = 1 * fc1 * totalscal

scales1 = cparam1 / np.arange(totalscal, 0, -1)


wavename2 = 'cmor1.5-2' #

fc2 = pywt.central_frequency(wavename2)

cparam2 = 2 * fc2 * totalscal

scales2 = cparam2 / np.arange(totalscal, 0, -1)


wavename3 = "cmor1.5-2"

fc3 = pywt.central_frequency(wavename3)

cparam3 = 3 * fc3 * totalscal ?

scales3 = cparam3 / np.arange(totalscal, 0, -1)


wavename4 = 'cmor1.5-2'

fc4 = pywt.central_frequency(wavename4)

cparam4 = 4 * fc4 * totalscal ?

scales4 = cparam4 / np.arange(totalscal, 0, -1)


# 进行连续小波变换

coefficients1, frequencies1 = pywt.cwt(data_list3, scales1, wavename1, sampling_period)

coefficients2, frequencies2 = pywt.cwt(data_list3, scales2, wavename2, sampling_period)

coefficients3, frequencies3 = pywt.cwt(data_list3, scales3, wavename3, sampling_period)

coefficients4, frequencies4 = pywt.cwt(data_list3, scales4, wavename4, sampling_period)


# 小波系数矩阵绝对值

amp1 = abs(coefficients1)

amp2 = abs(coefficients2)

amp3 = abs(coefficients3)

amp4 = abs(coefficients4)


# 根据采样频率 sampling_period 生成时间轴 t

t = np.linspace(0, 1.0/sampling_rate, sampling_rate, endpoint=False)

进行可视化

plt.figure(figsize=(20,10), dpi=300)

plt.subplot(2,2,1)

plt.contourf(t, frequencies1, amp1, cmap='jet')

plt.title('滚珠-32')

plt.subplot(2,2,2)

plt.contourf(t, frequencies2, amp2, cmap='jet')

plt.title('滚珠-64')

plt.subplot(2,2,3)

plt.contourf(t, frequencies3, amp3, cmap='jet')

plt.title('滚珠-128')

plt.subplot(2,2,4)

plt.contourf(t, frequencies4, amp4, cmap='jet')

plt.title('滚珠-256')

plt.show()

对于不同尺度序列scales的比较，更能说明之前的结论：

高频特征：如果关注信号的高频特征，应该选择较小的尺度；
低频特征：如果关注信号的低频特征，应该选择较大的尺度；

对于滚珠故障类型数据，从时频图结果来看，应该选择2倍的cparams参数，有着较高的频率分辨率，和我们感兴趣的频率区域。

2.4.4 比较cmor小波函数不同参数 ----中心频率，带宽参数

在轴承故障诊断中，中心频率和带宽的具体设置取决于多个因素，包括轴承类型、工作条件和故障特征等。由于每个应用场景和故障类型都有所不同，没有一个通用的固定数值。以下是一些常见的参考范围和建议：

中心频率（Center Frequency）：? ?

对于滚动轴承，常见的故障频率范围在几百赫兹到几千赫兹之间。可以选择在这个范围内的适当值作为中心频率。? ?
根据具体的故障类型，例如滚动体故障、内圈故障或外圈故障，可能需要设置不同的中心频率。
?

带宽（Bandwidth）：? ?

带宽的选择取决于信号的特征和所需的时频分辨率。? ?
通常情况下，较小的带宽值（如1-5）适用于较短时域特征，能够提供更好的时频局部化能力。??
?如果需要更好的频率分辨率，可以选择较大的带宽值（如10-20），但可能会牺牲一部分时频局部化能力。

这些值仅供参考，实际应用时需要进行实验和调整，根据信号的频谱特征和故障频率范围来确定最佳的参数配置。建议通过观察生成的时频图，确保故障频率和特征得到适当的捕捉和展示。同时，根据实际的故障案例和经验，不断调整参数以提高故障诊断的准确性和可靠性。

经过大量的对比实验和观察，本文得出最后的参数结论设置：

# 尺度长度

totalscal = 128 ? ?

# 小波基函数

wavename = 'cmor100-1'

# 小波函数中心频率

fc = pywt.central_frequency(wavename)

# 常数c

cparam = 2 * fc * totalscal ?

# 小波尺度序列

scales = cparam / np.arange(totalscal, 0, -1)

来实现对故障数据的诊断分类。

3?基于时频图像的轴承故障诊断分类

下面以连续小波变换（CWT）作为轴承故障数据的处理方法进行讲解：

数据介绍，凯斯西储大学（CWRU）轴承数据10分类数据集

train_set、val_set、test_set 均为按照7：2：1划分训练集、验证集、测试集，最后保存数据

3.1?生成时频图像数据集

如图所示为生成的时频图像数据集

3.2 定义数据加载器和VGG网络模型

制作数据标签，保存数据

定义VGG网络模型

3.3 设置参数，训练模型

30个epoch，准确率将近90%，继续调参可以进一步提高分类准确率

参考文献

[1]《小波分析及其工程应用》.机械工业出版社

文章来源:https://blog.csdn.net/qq_40949048/article/details/134758611
本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若内容造成侵权/违法违规/事实不符，请联系我的编程经验分享网邮箱：veading@qq.com进行投诉反馈，一经查实，立即删除！