Python軸承故障診斷 (二)連續小波變換CWT

前言

1?連續小波變換CWT原理介紹

1.1?CWT概述

1.2 CWT的原理和本質

2?基于Python的CWT實現與參數對比

2.1 代碼示例

2.2 參數介紹和選擇策略

2.2.1 尺度長度：

2.2.2 小波函數（wavelet）：

2.3 凱斯西儲大學軸承數據的加載

2.4 CWT與參數選擇對比

2.4.1 基于尺度為128，選擇內圈數據比較 CWT 的不同小波函數

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

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

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

3.1?生成時頻圖像數據集

3.2 定義數據加載器和VGG網絡模型

往期精彩內容：

Python房價分析（一）pyton爬蟲

時序預測：LSTM、ARIMA、Holt-Winters、SARIMA模型的分析與比較

Python-電能質量擾動信號數據介紹與分類

Python-凱斯西儲大學（CWRU）軸承數據解讀與分類處理

Python軸承故障診斷 (一)短時傅里葉變換STFT

前言

本文基于凱斯西儲大學（CWRU）軸承數據，進行連續小波變換的介紹與參數選擇，最后通過Python實現對故障數據的時頻圖像分類。凱斯西儲大學軸承數據的詳細介紹可以參考下文：

Python-凱斯西儲大學（CWRU）軸承數據解讀與分類處理

短時傅里葉變換STFT可以參考如下：? ?

Python軸承故障診斷 (一)短時傅里葉變換STFT

前期內容介紹了短時傅里葉變換（STFT），在其一定程度上改善了傅里葉變換（FT）無時間局部化能力的不足，但仍然存在以下不足：

在STFT中，選擇了固定的時間窗，時間和頻率的分辨率就保持不變；
由于STFT的基礎是傅里葉變換，所以它不適合于分析非平穩信號，可用于分析平穩信號和準平穩信號。

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

1?連續小波變換CWT原理介紹

1.1?CWT概述

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

CWT的公式表示為：

其中，

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

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

1.2 CWT的原理和本質

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

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

1.3 時頻圖譜

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

2?基于Python的CWT實現與參數對比

在 Python 中，使用 pywt 庫來實現連續小波變換（CWT）

2.1 代碼示例

import numpy as npimport matplotlib.pyplot as pltimport 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)# 常數ccparam = 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 生成時間軸 tt = 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 = 12num_columns = 10for 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 npfrom scipy.io import loadmatimport numpy as npfrom scipy.signal import stftimport matplotlib.pyplot as pltimport matplotlibmatplotlib.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 npimport matplotlib.pyplot as pltimport pywtimport 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 生成時間軸 tt = 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 npimport matplotlib.pyplot as pltimport pywtimport pandas as pd# 連續小波變換參數# 采樣頻率sampling_rate = 1024# 尺度長度totalscal = 128 ?wavename1 = 'cmor1.5-2'fc1 = pywt.central_frequency(wavename1)cparam1 = 1 * fc1 * totalscalscales1 = cparam1 / np.arange(totalscal, 0, -1)wavename2 = 'cmor1.5-2' #fc2 = pywt.central_frequency(wavename2)cparam2 = 2 * fc2 * totalscalscales2 = 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 生成時間軸 tt = 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)# 常數ccparam = 2 * fc * totalscal ?# 小波尺度序列scales = cparam / np.arange(totalscal, 0, -1)

來實現對故障數據的診斷分類。