代碼位置：https://github.com/HappyBoy-cmd/fault_diagnosis_signal_processing
參考資料：
機械故障診斷及典型案例解析（第2版，時獻江）
會議論文：Bearing Intelligent Fault Diagnosis Under Complex Working Condition Based on SK-ES-CNN，2021 Global Reliability and Prognostics and Health Management (PHM-Nanjing)
包括內圈、外圈、滾動體和正常數據，分別為一個。
1730_12k_0.007-InnerRace：內圈故障
1730_12k_0.007-OuterRace3：外圈故障
1730_12k_0.014-Ball：滾動體故障
1730_48k_Normal：正常數據
對CWRU軸承數據集不了解的同學見這裡：
CWRU數據集介紹 第1期
CWRU數據集介紹 第2期
CWRU數據集介紹 第3期
CWRU數據集介紹 第3期
下面先以軸承內圈故障數據進行分析。原始數據為mat文件，是matlab文件，定義一個函數進行數據讀取
def data_acquision(FilePath):
""" fun: 從cwru mat文件讀取加速度數據 param file_path: mat文件絕對路徑 return accl_data: 加速度數據，array類型 """
data = scio.loadmat(file_path) # 加載mat數據
data_key_list = list(data.keys()) # mat文件為字典類型，獲取字典所有的鍵並轉換為list類型
accl_key = data_key_list[3] # 獲取'X108_DE_time'
accl_data = data[accl_key].flatten() # 獲取'X108_DE_time'所對應的值，即為振動加速度信號,並將二維數組展成一維數組
return accl_data
data_acquision(FilePath)輸入參數FilePath，輸出一個1維array數據。下面進行演示
file_path = r'E:/研究生/pytorch/CSDN代碼/fault_diagnosis_signal_processing/第4篇-包絡譜/1730_12k_0.007-InnerRace.mat'
xt = data_acquision(file_path)
plt.figure(figsize=(12,3))
plt.plot(xt)
print(xt)
輸出結果：
[ 0.22269856 0.09323776 -0.14651649 ... -0.36125573 0.31138814
0.17055689]
希爾伯特解調法，亦叫包絡譜分析。
設 x ( t ) x(t) x(t)為一個實時域信號，其Hilbert變換定義為：
h ( t ) = 1 π ∫ − ∞ + ∞ x ( τ ) t − τ d τ = x ( t ) ∗ 1 π t h(t)=\frac{1}{\pi} \int_{-\infty}^{+\infty} \frac{x(\tau)}{t-\tau} \mathrm{d} \tau=x(t) * \frac{1}{\pi t} h(t)=π1∫−∞+∞t−τx(τ)dτ=x(t)∗πt1
則原始信號 x ( t ) x(t) x(t)和它的Hilbert變換信號 h ( t ) h(t) h(t)可以構建一個新的解析信號 z ( t ) z(t) z(t):
z ( t ) = x ( t ) + j h ( t ) = a ( t ) e j φ t z(t)=x(t)+j h(t)=a(t) e^{j \varphi t} z(t)=x(t)+jh(t)=a(t)ejφt
# step1: 做希爾伯特變換
ht = fftpack.hilbert(xt)
print(ht)
輸出結果：
[-0.02520403 -0.28707983 -0.00610516 ... 0.1100125 0.22821944
-0.11203138]
此時輸出的 h ( t ) h(t) h(t)是解析信號 a ( t ) a(t) a(t)的虛部系數
對 z ( t ) z(t) z(t)取模，得到其幅值 a ( t ) : a(t): a(t):
a ( t ) = ∣ z ( t ) ∣ = x 2 ( t ) + h 2 ( t ) {a(t)=|z(t)|=\sqrt{x^{2}(t)+h^{2}(t)}} a(t)=∣z(t)∣=x2(t)+h2(t)
注： a ( t ) a(t) a(t)即為包絡信號，也叫解析信號
sampling_rate = 12000
am = np.fft.fft(at) # 對希爾伯特變換後的at做fft變換獲得幅值
am = np.abs(am) # 對幅值求絕對值（此時的絕對值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # 獲取fft頻率，此時包括正頻率和負頻率
freq = freq[0:int(len(freq)/2)] # 獲取正頻率
plt.plot(freq, am)
觀察發現：
(1)在頻率為0Hz的地方幅值比較大
(2)在低頻部分貌似看到一倍頻和二倍頻
在0Hz的幅值比較高，使得其它頻率幅值較低，不便觀察。這種現象叫直流分量，去直流分量方法，y = y-mean(y)
sampling_rate = 12000
at = at - np.mean(at) # 去直流分量
am = np.fft.fft(at) # 對希爾伯特變換後的at做fft變換獲得幅值
am = np.abs(am) # 對幅值求絕對值（此時的絕對值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # 獲取fft頻率，此時包括正頻率和負頻率
freq = freq[0:int(len(freq)/2)] # 獲取正頻率
plt.plot(freq, am)
輸出結果：
sampling_rate = 12000
at = at - np.mean(at) # 去直流分量
am = np.fft.fft(at) # 對希爾伯特變換後的at做fft變換獲得幅值
am = np.abs(am) # 對幅值求絕對值（此時的絕對值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # 獲取fft頻率，此時包括正頻率和負頻率
freq = freq[0:int(len(freq)/2)] # 獲取正頻率
plt.figure(figsize=(12,3))
plt.plot(freq, am)
plt.xlim(0,500)
輸出結果：
內圈故障特征頻率： F B P F I = n f r 2 ( 1 + d D cos α ) F_{\mathrm{BPFI}}=\frac{n f_{r}}{2}\left(1+\frac{d}{D} \cos \alpha\right) FBPFI=2nfr(1+Ddcosα)
外圈故障特征頻率： F B P F O = n f r 2 ( 1 − d D cos α ) F_{\mathrm{BPFO}}=\frac{n f_{r}}{2}\left(1-\frac{d}{D} \cos \alpha\right) FBPFO=2nfr(1−Ddcosα)
滾動體故障特征頻率： F B S F = D f r 2 d [ 1 − ( d D cos α ) 2 ] F_{\mathrm{BSF}}=\frac{D f_{r}}{2 d}\left[1-\left(\frac{d}{D} \cos \alpha\right)^{2}\right] FBSF=2dDfr[1−(Ddcosα)2]
n n n: 滾動體個數， f r f_{r} fr: 軸轉速 d d d: 滾珠(子)直徑 D D D: 軸承節徑
軸承型號為：6205-2RSL JME SKF 深溝球滾珠軸承
d d d=7.94mm, D D D=39.04mm, α \alpha α=0, n n n=9
為了方便，定義了一個軸承故障特征頻率計算函數，只需輸入參數 d d d, D D D, α \alpha α, n n n和 f r f_{r} fr即可
def bearing_fault_freq_cal(n, d, D, alpha, fr=None):
''' 基本描述： 計算滾動軸承的故障特征頻率 詳細描述： 輸入4個參數 n, fr, d, D, alpha return C_bpfi, C_bpfo, C_bsf, C_ftf, fr 內圈 外圈 滾針 保持架 轉速 Parameters ---------- n: integer The number of roller element fr: float(r/min) Rotational speed d: float(mm) roller element diameter D: float(mm) pitch diameter of bearing alpha: float(°) contact angle fr:：float(r/min) rotational speed Returns ------- BPFI: float(Hz) Inner race-way fault frequency BPFO: float(Hz) Outer race-way fault frequency BSF: float(Hz) Ball fault frequency FTF: float(Hz) Cage frequency '''
C_bpfi = n*(1/2)*(1+d/D*np.math.cos(alpha))
C_bpfo = n*(1/2)*(1-(d/D)*np.math.cos(alpha))
C_bsf = D*(1/(2*d))*(1-np.square(d/D*np.math.cos(alpha)))
C_ftf = (1/2)*(1-(d/D)*np.math.cos(alpha))
if fr!=None:
return C_bpfi*fr/60, C_bpfo*fr/60, C_bsf*fr/60, C_ftf*fr/60, fr/60
else:
return C_bpfi, C_bpfo, C_bsf, C_ftf, fr
下面計算CWRU在轉速1730rpm時的故障特征頻率
bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1730)
print('內圈故障特征頻率',bpfi)
print('外圈故障特征頻率',bpfo)
print('滾動體故障特征頻率',bsf)
print(ftf)
print(fr)
sampling_rate = 12000
at = at - np.mean(at) # 去直流分量
am = np.fft.fft(at) # 對希爾伯特變換後的at做fft變換獲得幅值
am = np.abs(am) # 對幅值求絕對值（此時的絕對值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # 獲取fft頻率，此時包括正頻率和負頻率
freq = freq[0:int(len(freq)/2)] # 獲取正頻率
plt.figure(figsize=(12,3))
plt.plot(freq, am)
plt.xlim(0,500)
plt.vlines(x=156.13, ymin=0, ymax=0.2, colors='r') # 一倍頻
plt.vlines(x=156.13*2, ymin=0, ymax=0.2, colors='r') # 二倍頻
輸出結果：
紅色為理論內圈故障特征頻率，藍線為實際故障特征頻率。雖然沒有完全重合，但這個是在允許范圍內的。因此在實際情況中包絡譜出現該情況，可以該軸承出現了內圈故障。
那如果不希爾伯特變換解調，直接fft，能夠看到故障特征頻率嗎？下面進行對比分析
sampling_rate = 12000
am = np.fft.fft(xt) # 對希爾伯特變換後的at做fft變換獲得幅值
am = np.abs(am) # 對幅值求絕對值（此時的絕對值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(xt), d=1 / sampling_rate) # 獲取fft頻率，此時包括正頻率和負頻率
freq = freq[0:int(len(freq)/2)] # 獲取正頻率
plt.plot(freq, am)
輸出結果：
可見原始內圈故障信號的fft高頻較多
plt.plot(freq, am)
plt.xlim(0, 500)
plt.vlines(x=156.13, ymin=0, ymax=0.05, colors='r') # 一倍頻
plt.vlines(x=156.13*2, ymin=0, ymax=0.05, colors='r') # 二倍頻
可見直接fft的話，故障特征頻率幅值較低，被其它頻率幅值掩蓋了。反過來，希爾伯特解調可以更加方便觀察故障特征頻率低頻。
為了更加方便使用包絡譜，這裡封裝了一個包絡譜函數plt_envelope_spectrum
def plt_envelope_spectrum(data, fs, xlim=None, vline= None):
''' fun: 繪制包絡譜圖 param data: 輸入數據，1維array param fs: 采樣頻率 param xlim: 圖片橫坐標xlim，default = None param vline: 圖片垂直線，default = None '''
#----去直流分量----#
data = data - np.mean(data)
#----做希爾伯特變換----#
xt = data
ht = fftpack.hilbert(xt)
at = np.sqrt(xt**2+ht**2) # 獲得解析信號at = sqrt(xt^2 + ht^2)
am = np.fft.fft(at) # 對解析信號at做fft變換獲得幅值
am = np.abs(am) # 對幅值求絕對值（此時的絕對值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)] # 取正頻率幅值
freq = np.fft.fftfreq(len(at), d=1 / fs) # 獲取fft頻率，此時包括正頻率和負頻率
freq = freq[0:int(len(freq)/2)] # 獲取正頻率
plt.plot(freq, am)
if vline: # 是否繪制垂直線
plt.vlines(x=vline, ymax=0.2, ymin=0, colors='r') # 高度y 0-0.2，顏色紅色
if xlim: # 圖片橫坐標是否設置xlim
plt.xlim(0, xlim)
plt.xlabel('freq(Hz)') # 橫坐標標簽
plt.ylabel('amp(m/s2)') # 縱坐標標簽
file_path = r'E:/研究生/pytorch/CSDN代碼/fault_diagnosis_signal_processing/第4篇-包絡譜/1730_12k_0.007-OuterRace3.mat'
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bpfo)
輸出結果：
可見實際外圈故障特征頻率比較明顯
file_path = r'E:/研究生/pytorch/CSDN代碼/fault_diagnosis_signal_processing/第4篇-包絡譜/1730_12k_0.014-Ball.mat'
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bsf)
可見實際滾動體故障特征頻率不明顯
(1)包絡譜能夠檢測出內圈、外圈故障，滾動體比較困難
(2)直接使用fft難以檢測故障特征頻率，故障特征頻率易被高頻掩蓋
