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

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

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

1.1 CWT概述

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

CWT的公式表示為：

其中，

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

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

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

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

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

1.3 時(shí)頻圖譜

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

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

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

2.1 代碼示例

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


# 生成三個(gè)不同頻率成分的信號(hào)  3000個(gè)點(diǎn)
fs = 1000  # 采樣率
time = np.linspace(0, 1, fs, endpoint=False)  # 時(shí)間
# 第一個(gè)頻率成分
signal1 = np.sin(2 * np.pi * 30 * time)
# 第二個(gè)頻率成分
signal2 = np.sin(2 * np.pi * 60 * time)
# 第三個(gè)頻率成分
signal3 = np.sin(2 * np.pi * 120 * time)
# 合并三個(gè)信號(hào)
signal = np.concatenate((signal1, signal2, signal3))


# 連續(xù)小波變換參數(shù)
# 采樣頻率
sampling_rate = 3000
# 尺度長(zhǎng)度
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ù)矩陣絕對(duì)值
amp = abs(coefficients)


# 根據(jù)采樣頻率 sampling_period 生成時(shí)間軸 t
t = np.linspace(0, 1.0/sampling_rate, sampling_rate, endpoint=False)
# 繪制時(shí)頻圖譜
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('對(duì)應(yīng)時(shí)頻圖')


plt.show()

參數(shù)解釋

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

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

2.2.1 尺度長(zhǎng)度：

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

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

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

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

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

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

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

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

第一步，導(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')  # 正常信號(hào)
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ū)動(dòng)端加速度數(shù)據(jù)，取一個(gè)長(zhǎng)度為1024的信號(hào)進(jìn)行后續(xù)觀察和實(shí)驗(yàn)

# DE - drive end accelerometer data 驅(qū)動(dòng)端加速度數(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ù)選擇對(duì)比

本實(shí)驗(yàn)以某軸承故障診斷論文推薦'cgau8'小波，以及'morl'小波、'cmor1-1'小波、'cmor1.5-2'小波為實(shí)驗(yàn)做對(duì)比，尺度先設(shè)定為128，來對(duì)比不同小波函數(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
# 尺度長(zhǎng)度
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ù)矩陣絕對(duì)值
amp1 = abs(coefficients1)
amp2 = abs(coefficients2)
amp3 = abs(coefficients3)
amp4 = abs(coefficients4)


# 根據(jù)采樣頻率 sampling_period 生成時(shí)間軸 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()

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

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

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

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

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

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

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

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


# 連續(xù)小波變換參數(shù)
# 采樣頻率
sampling_rate = 1024
# 尺度長(zhǎng)度
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ù)矩陣絕對(duì)值
amp1 = abs(coefficients1)
amp2 = abs(coefficients2)
amp3 = abs(coefficients3)
amp4 = abs(coefficients4)


# 根據(jù)采樣頻率 sampling_period 生成時(shí)間軸 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()

對(duì)于不同尺度序列scales的比較，更能說明之前的結(jié)論：

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

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

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

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

中心頻率（Center Frequency）：

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

帶寬（Bandwidth）：

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

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

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

# 尺度長(zhǎng)度
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)

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

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

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

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

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

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

如圖所示為生成的時(shí)頻圖像數(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個(gè)epoch，準(zhǔn)確率將近90%，繼續(xù)調(diào)參可以進(jìn)一步提高分類準(zhǔn)確率。