Python軸承故障診斷之連續(xù)小波變換CWT介紹

軸承故障數(shù)據(jù)的信號特性常隨時間變化，對于這種非平穩(wěn)信號分析，CWT 提供了一種同時在時間和頻率上定位信號特征的方法。

1 連續(xù)小波變換CWT原理介紹

1.1 CWT概述

連續(xù)小波變換（Continuous Wavelet Transform，CWT）是一種用于在時域和頻域上同時分析信號的方法，它通過使用不同尺度和位置的小波函數(shù)對信號進(jìn)行變換，以獲取信號的局部特性。

CWT的公式表示為：

其中，

信號x(t)經(jīng)過小波變換后，得到的結(jié)果是小波系數(shù)C，小波系數(shù)C是尺度a和位置b的函數(shù)。從物理意義上講，小波系數(shù)C中蘊含著信號在各個尺度a和位置b上的信息[1]。

不同尺度和位置下小波的形狀變換如圖所示：

1.2 CWT的原理和本質(zhì)

CWT的核心思想是在不同尺度（頻率）和位置上對信號進(jìn)行小波分解。為了達(dá)到這個目的，CWT使用一個小波函數(shù)（wavelet），通常稱為母小波或基本小波。這個小波函數(shù)是一個可調(diào)整尺度的波形。

• 尺度和平移：CWT使用可調(diào)整尺度的小波函數(shù)，這個尺度參數(shù)決定了小波的頻率，同時也使用平移參數(shù)，控制小波在時間軸上的位置。
• 小波函數(shù)：小波函數(shù)通常稱為母小波，是一種局部化的波形，通過縮放和平移，可以適應(yīng)信號的不同頻率和位置。
• 卷積過程：CWT通過在不同尺度和位置上對信號進(jìn)行小波函數(shù)的卷積，生成一系列的小波系數(shù)。

1.3 時頻圖譜

連續(xù)小波變換的結(jié)果是小波系數(shù)，提供了信號在時間和頻率上的局部信息。這些小波系數(shù)構(gòu)成了時頻平面上的圖像，被稱為時頻圖譜。時頻圖譜顯示了信號在時間和頻率上的局部特性。這對于定位故障信號中的異常事件以及了解信號的時頻結(jié)構(gòu)非常有用。

2 基于Python的CWT實現(xiàn)與參數(shù)對比

在 Python 中，使用 pywt 庫來實現(xiàn)連續(xù)小波變換（CWT）

2.1 代碼示例

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))


# 連續(xù)小波變換參數(shù)
# 采樣頻率
sampling_rate = 3000
# 尺度長度
totalscal = 128    
# 小波基函數(shù)
wavename = 'morl'
# 小波函數(shù)中心頻率
fc = pywt.central_frequency(wavename)
# 常數(shù)c
cparam = 2 * fc * totalscal  
# 尺度序列
scales = cparam / np.arange(totalscal, 0, -1)
 # 進(jìn)行CWT連續(xù)小波變換
coefficients, frequencies = pywt.cwt(signal, scales, wavename, 1.0/1000)
# 小波系數(shù)矩陣絕對值
amp = abs(coefficients)


# 根據(jù)采樣頻率 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('對應(yīng)時頻圖')


plt.show()

參數(shù)解釋

• totalscal：表示尺度長度，表征頻率的參數(shù)，選擇不同的尺度可以捕捉信號不同尺度上的特征；
• wavename：表示小波函數(shù)的類型，不同的小波適用于不同類型的信號；
• fc：表示小波函數(shù)的中心頻率，小波函數(shù)在頻率域中有一個中心頻率；這個中心頻率是與小波函數(shù)形狀和性質(zhì)有關(guān)的一個重要參數(shù)；
• scales：表示尺度序列，CWT本質(zhì)上是將你的信號與不同尺度的小波進(jìn)行相關(guān)，scales 參數(shù)確定尺度范圍；
• coefficients：表示信號變換后的小波系數(shù)；
• frequencies ：表示對應(yīng)的頻率信息

2.2 參數(shù)介紹和選擇策略

2.2.1 尺度長度：

在連續(xù)小波變換（CWT）中，尺度參數(shù)是一個關(guān)鍵的選擇，因為它決定了小波函數(shù)的寬度，從而影響了頻率分辨率。尺度與頻率成反比，尺度反映了分析的頻率范圍，尺度越小，小波函數(shù)衰減越快，頻率越高；尺度越大，小波函數(shù)衰減越慢，頻率越低[1]。

選擇小波尺度的一般原則是：

• 高頻特征：如果關(guān)注信號的高頻特征，應(yīng)該選擇較小的尺度；
• 低頻特征：如果關(guān)注信號的低頻特征，應(yīng)該選擇較大的尺度；
• 覆蓋感興趣的頻率范圍：尺度參數(shù)的選擇應(yīng)該使小波函數(shù)能夠覆蓋感興趣的頻率范圍；如果期望信號有很高的頻率變化，可能需要選擇較小的尺度；
• 頻率分辨率 ： 較小的尺度提供更好的頻率分辨率 ，因為小波函數(shù)較窄，可以更精細(xì)地定位頻率；但是，這也意味著在時間上的分辨率較差，因此， 需要權(quán)衡時間分辨率和頻率分辨率 ；
• 信號持續(xù)時間 ：尺度參數(shù)的選擇還應(yīng)考慮信號的持續(xù)時間，如果信號是短暫的，可能需要較小的尺度；
• 尺度間隔： 在尺度參數(shù)上選擇合適的間隔，以確保在整個頻率范圍內(nèi)進(jìn)行了適當(dāng)?shù)牟蓸?；這取決于具體的應(yīng)用和信號特性。

2.2.2 小波函數(shù)（wavelet）：

小波函數(shù)（wavelet）的選擇也連續(xù)小波變換中的一個重要參數(shù)，它決定了小波基函數(shù)的形狀，不同的小波函數(shù)適用于不同類型的信號和應(yīng)用。

打印 Python pywt 包中的所有小波函數(shù)類型：

import pywt


# 獲取小波函數(shù)列表
wavelets = pywt.wavelist()


# 打印小波函數(shù)列表
# 按照12行，10列的形式打印數(shù)據(jù)
num_rows = 12
num_columns = 10


for i in range(0, len(wavelets), num_columns):
    row = wavelets[i:i + num_columns]
    print(row)

介紹常用的小波基函數(shù)：

• 'morl' ：Morlet小波是一種復(fù)雜的小波函數(shù)，它在頻率域和時域都有較好的局部化性質(zhì)；Morlet小波通常用于處理時頻局部化要求較高的信號，比如處理振動信號或某些生物醫(yī)學(xué)信號；
• 'cmor': Complex Morlet wavelet，復(fù)數(shù) Morlet 小波，是 Morlet 小波的一個變種；
• 'cgau' ：Complex Gaussian wavelet，復(fù)數(shù)高斯小波，用于近似高斯信號；
• 'db1'：Daubechies小波是離散小波變換（Discrete Wavelet Transform, DWT）中常用的小波函數(shù)。Daubechies小波是緊支撐的小波，適用于處理有限長度的信號。
• 'haar'：Haar小波是最簡單的小波函數(shù)之一，適用于對信號進(jìn)行基本的低通和高通分解；
• 'mexh'：Mexican Hat小波，也稱為Ricker小波，適用于處理具有尖峰或波包特性的信號；
• 'bior1.1'：Bior小波是一類雙正交小波，其特點是具有對稱和非對稱兩組濾波器，尤其適用于一些信號的多分辨率分析，如圖像處理；
• 'sym2'：Symlet小波（Symmetric Wavelets），是一類對稱的小波函數(shù)，它在某些方面類似于Daubechies小波，但是Symlet小波在設(shè)計上更加靈活。Symlet小波也是一種緊支撐小波，適用于有限長度的信號處理。

小波函數(shù)的選擇通常取決于處理的信號類型以及分析的目標(biāo)。在實際應(yīng)用中，可以嘗試不同的小波函數(shù)，觀察它們在信號上的效果，然后根據(jù)實驗結(jié)果選擇最適合的小波函數(shù)。在選擇小波函數(shù)時，也要考慮小波函數(shù)的性質(zhì)，如平滑性、局部化等。

2.3 凱斯西儲大學(xué)軸承數(shù)據(jù)的加載

選擇正常信號和 0.021英寸內(nèi)圈、滾珠、外圈故障信號數(shù)據(jù)來做對比

第一步，導(dǎo)入包，讀取數(shù)據(jù)

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英寸 內(nèi)圈
data3 = loadmat('21_2.mat') # 0.021英寸 滾珠
data4 = loadmat('21_3.mat') # 0.021英寸 外圈
# 注意，讀取出來的data是字典格式，可以通過函數(shù)type(data)查看。

第二步，數(shù)據(jù)集中統(tǒng)一讀取 驅(qū)動端加速度數(shù)據(jù)，取一個長度為1024的信號進(jìn)行后續(xù)觀察和實驗

# DE - drive end accelerometer data 驅(qū)動端加速度數(shù)據(jù)
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)
# 劃窗取值（大多數(shù)窗口大小為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]

第三步，進(jìn)行數(shù)據(jù)可視化

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('內(nèi)圈')
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與參數(shù)選擇對比

本實驗以某軸承故障診斷論文推薦'cgau8'小波，以及'morl'小波、'cmor1-1'小波、'cmor1.5-2'小波為實驗做對比，尺度先設(shè)定為128，來對比不同小波函數(shù)的影響。

2.4.1 基于尺度為128，選擇內(nèi)圈數(shù)據(jù)比較 CWT 的不同小波函數(shù)

import numpy as np
import matplotlib.pyplot as plt
import pywt
import pandas as pd


# 連續(xù)小波變換參數(shù)
# 采樣頻率
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)


# 進(jìn)行連續(xù)小波變換
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)


# 小波系數(shù)矩陣絕對值
amp1 = abs(coefficients1)
amp2 = abs(coefficients2)
amp3 = abs(coefficients3)
amp4 = abs(coefficients4)


# 根據(jù)采樣頻率 sampling_period 生成時間軸 t
t = np.linspace(0, 1.0/sampling_rate, sampling_rate, endpoint=False)

數(shù)據(jù)可視化

plt.figure(figsize=(20,10))
plt.subplot(2,2,1)
plt.contourf(t, frequencies1, amp1, cmap='jet')
plt.title('內(nèi)圈-cgau8')
plt.subplot(2,2,2)
plt.contourf(t, frequencies2, amp2, cmap='jet')
plt.title('內(nèi)圈-morl')
plt.subplot(2,2,3)
plt.contourf(t, frequencies3, amp3, cmap='jet')
plt.title('內(nèi)圈-cmor1-1')
plt.subplot(2,2,4)
plt.contourf(t, frequencies4, amp4, cmap='jet')
plt.title('內(nèi)圈-cmor1.5-2')
plt.show()

對比不同小波函數(shù)，對于內(nèi)圈故障信號來說，'cgau8'小波有著較高的頻率分辨率，需要對比其他類型故障數(shù)據(jù)，進(jìn)一步觀察。

2.4.2 根據(jù)正常數(shù)據(jù)和三種故障數(shù)據(jù)，對比不同小波函數(shù)的辨識度

比較來看，'cgau8'小波和'cmor1.5-2'小波對于不同故障都有著較高的辨識度，綜合考慮頻率分辨率和時間分辨率，選擇'cmor1.5-2'小波來進(jìn)一步分析。

2.4.3 基于'cmor1.5-2'小波，選擇滾珠故障數(shù)據(jù)比較 CWT 的不同尺度的變化：32、64、128、256；

同時尺度序列scales的常數(shù)項cparams均采用：cparam = 2 * fc * totalscal；'cmor1.5-2'小波中1 代表中心頻率參數(shù)，1.5代表帶寬參數(shù)；

注意，在連續(xù)小波變換中，影響小波系數(shù)的是尺度序列scales，但是仍能從上圖中看出尺度長度越大，對低頻特征關(guān)注越多（注意每幅小圖 的底部低頻的區(qū)別）；

為進(jìn)一步探索尺度序列scales的影響，選擇尺度長度為128，設(shè)置不同常數(shù)項cparams來觀察對時頻變換的影響：

import numpy as np
import matplotlib.pyplot as plt
import pywt
import pandas as pd


# 連續(xù)小波變換參數(shù)
# 采樣頻率
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)


# 進(jìn)行連續(xù)小波變換
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)


# 小波系數(shù)矩陣絕對值
amp1 = abs(coefficients1)
amp2 = abs(coefficients2)
amp3 = abs(coefficients3)
amp4 = abs(coefficients4)


# 根據(jù)采樣頻率 sampling_period 生成時間軸 t
t = np.linspace(0, 1.0/sampling_rate, sampling_rate, endpoint=False)

進(jìn)行可視化

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的比較，更能說明之前的結(jié)論：

• 高頻特征：如果關(guān)注信號的高頻特征，應(yīng)該選擇較小的尺度；
• 低頻特征： 如果關(guān)注信號的低頻特征，應(yīng)該選擇較大的尺度；

對于滾珠故障類型數(shù)據(jù)，從時頻圖結(jié)果來看，應(yīng)該選擇2倍的cparams參數(shù)，有著較高的頻率分辨率，和我們感興趣的頻率區(qū)域。

2.4.4 比較cmor小波函數(shù) 不同參數(shù) ----中心頻率，帶寬參數(shù)

在軸承故障診斷中，中心頻率和帶寬的具體設(shè)置取決于多個因素，包括軸承類型、工作條件和故障特征等。由于每個應(yīng)用場景和故障類型都有所不同，沒有一個通用的固定數(shù)值。以下是一些常見的參考范圍和建議：

中心頻率（Center Frequency）：

• 對于滾動軸承，常見的故障頻率范圍在幾百赫茲到幾千赫茲之間?？梢赃x擇在這個范圍內(nèi)的適當(dāng)值作為中心頻率。
• 根據(jù)具體的故障類型，例如滾動體故障、內(nèi)圈故障或外圈故障，可能需要設(shè)置不同的中心頻率。

帶寬（Bandwidth）：

• 帶寬的選擇取決于信號的特征和所需的時頻分辨率。
• 通常情況下，較小的帶寬值（如1-5）適用于較短時域特征，能夠提供更好的時頻局部化能力。
• 如果需要更好的頻率分辨率，可以選擇較大的帶寬值（如10-20），但可能會犧牲一部分時頻局部化能力。

這些值僅供參考，實際應(yīng)用時需要進(jìn)行實驗和調(diào)整，根據(jù)信號的頻譜特征和故障頻率范圍來確定最佳的參數(shù)配置。建議通過觀察生成的時頻圖，確保故障頻率和特征得到適當(dāng)?shù)牟蹲胶驼故?。同時，根據(jù)實際的故障案例和經(jīng)驗，不斷調(diào)整參數(shù)以提高故障診斷的準(zhǔn)確性和可靠性。

經(jīng)過大量的對比實驗和觀察，本文得出最后的參數(shù)結(jié)論設(shè)置：

# 尺度長度
totalscal = 128    
# 小波基函數(shù)
wavename = 'cmor100-1'
# 小波函數(shù)中心頻率
fc = pywt.central_frequency(wavename)
# 常數(shù)c
cparam = 2 * fc * totalscal  
# 小波尺度序列
scales = cparam / np.arange(totalscal, 0, -1)

來實現(xiàn)對故障數(shù)據(jù)的診斷分類。

3 基于時頻圖像的軸承故障診斷分類

下面以連續(xù)小波變換（CWT）作為軸承故障數(shù)據(jù)的處理方法進(jìn)行講解：

數(shù)據(jù)介紹，凱斯西儲大學(xué)（CWRU）軸承數(shù)據(jù)10分類數(shù)據(jù)集

train_set、val_set、test_set 均為按照7：2：1劃分訓(xùn)練集、驗證集、測試集，最后保存數(shù)據(jù)

3.1 生成時頻圖像數(shù)據(jù)集

如圖所示為生成的時頻圖像數(shù)據(jù)集

3.2 定義數(shù)據(jù)加載器和VGG網(wǎng)絡(luò)模型

制作數(shù)據(jù)標(biāo)簽，保存數(shù)據(jù)

定義VGG網(wǎng)絡(luò)模型

3.3 設(shè)置參數(shù)，訓(xùn)練模型

30個epoch，準(zhǔn)確率將近90%，繼續(xù)調(diào)參可以進(jìn)一步提高分類準(zhǔn)確率。