在機械加工過程中，刀具銑削過程傳感器監控資料均為時序資料，如下圖所示，表現為一系列沒有固定規律的随機離散值，不能有效識别資料特征與磨損度的關系，是以需要進行特征提取，把原有的低維時序資料映射到高維表示，挖掘和構造輸入資料的新特征，進而找出最能幫助模型拟合的輸入(X)和輸出(Y)特征。
在很多預測性維護的場景(如旋轉件、振動等)中，我們完全可以采用基于統計和信号的特征提取方法就能得到一些較好的特征。本文的特征提取方法主要包含三個角度：
- 基于統計的時域特征
- 基于頻譜分析的頻域特征
- 基于小波包能量的時頻聯合域特征
1 基于統計的時域特征
基于統計的時域特征，就是使用統計學方法如均方根、方差、最大值、最小值、偏斜度、峰度、峰峰值等提取出的新資料。
1.1 均方根(Root Mean Square, RMS)
均方根(Root Mean Square, RMS)，是信号有效值的反映：
1.2 方差(Variance)
方差(Variance)，衡量随機變量或一組資料離散程度，是源資料和期望值相差的度量：
1.3 最大值(Max)和最小值(Min)
最大值(Max)和最小值(Min)，表示在一定時間範圍内資料的最強和最弱的程度：
1.4 偏度(Skewness)
偏度(Skewness)又稱偏态系數，度量資料的機率密度曲線分布偏斜的方向和程度，即與平均值相比資料的非對稱的程度特征。
其中， 是3階中心距， 為标準差。
1.5 峰峰值(Peak to Peak, PP)
峰峰值(Peak to Peak, PP)，即一個周期内資料最高值和最低值之間的內插補點，描述了資料的變化範圍的大小：
1.6 峰度(Kurtosis)
峰度(Kurtosis)又稱峰态系數，表示在平均值處資料的機率密度分布曲線的最大值大小。峰度直覺地反映了曲線的峰部尖度，峰度越高，說明曲線底部厚度越大，會導緻方差越大，原因是資料存在一定數量的顯著比平均值大或小的內插補點。
其中， 是四階中心矩， 是标準差。
2 基于頻譜分析的頻域特征
頻域分析是信号領域中很重要的分析方法，具體包括了頻譜分析、能量譜分析、功率譜分析和倒頻譜分析等，其中以頻譜分析最為常用也最為重要。頻譜分析主要将信号資料通過傅裡葉變換得到幅頻譜和相頻譜進行分析，有關傅裡葉變換的内容網上有很多，就不贅述了。
通過計算機和傳感器采集得到的信号是離散信号，可以采用離散傅裡葉變換：
信号的幅值和相位可以表示為：
分别以幅值和相位為縱坐标，以頻率作為橫坐标，繪制出的二維坐标圖為幅值譜圖和相位譜圖，幅值譜圖表征信号的幅值随頻率的分布情況，相位譜圖表征了信号的相位随頻率的變化情況。在銑刀磨損度預測任務中，可以通過幅值譜提取以下統計特征：
2.1 譜偏态系數
譜偏态系數(Spectral Skewness)：統計頻域上幅值的分布偏斜的方向和程度：
2.2 譜峰态系數
譜峰态系數(Spectral Kurtosis)：表示頻域上幅值波形的尖銳程度：
2.3 譜功率
譜功率(Spectral Power)：表示機關頻帶内的信号功率：
3 基于小波包能量的時頻聯合域特征
小波變換(wavelet transform，WT)是一種新的變換分析方法，它繼承和發展了短時傅立葉變換局部化的思想，同時又克服了視窗大小不随頻率變化等缺點，能夠提供一個随頻率改變的“時間-頻率”視窗，是進行信号時頻分析和處理的理想工具。
小波變換過程簡述如下：
- 原始輸入信号X通過一組小波函數基經低通濾波器和高通濾波器分解後進行下采樣，得到低頻分量(又稱為近似子帶)和高頻分量(又稱為細節子帶)，
- 将低頻分量作為輸入信号，又進行小波分解，得到下一層的低頻分量和高頻分量，以此類推，可以得到層的小波分解，如圖 (a)所示。 随着小波分解層數的增加，其在頻域上的分辨率越高，這被稱為多分辨率分析(MultiResolution Analysis, MRA)。
小波包分解(Wavelet Packet Decomposition)又稱為最優子帶樹結構(Optimal Subband Tree Structuring)，如圖(b)所示，是對小波變換的進一步優化，具體為：在小波變換的基礎上，也對細節子帶進一步分解，最後通過最小化代價函數，計算出最優的信号分解路徑，并以此路徑對原始輸入信号進行分解。
選擇合适的小波函數基是進行信号分解的關鍵步驟。常用的小波函數基有哈爾小波函數基、多貝西小波和雙正交小波等。
多貝西小波作為稀疏基引入信号的光滑誤差不容易被察覺，是以具有較好的正則性，信号重構過程比較光滑，而且可以通過階次N的來控制頻域的局部化能力，使用較為靈活，是以本文的小波包分解過程采用多貝西小波函數基。
4 特征提取實際效果
介紹完枯燥的原理後，我們通過實際操作來看看這麼常用的特征提取方法最終的效果如何。
4.1 時域特征提取效果
圖中展示了PHM2010刀具磨損預測資料集C1資料的X軸的振動信号在銑刀磨損周期内的6個統計學特征與磨損度的關系圖。
4.2 頻域特征提取效果
頻域特征首先對預處理後的原始資料進行FFT處理，得到幅值譜圖，然後分别計算譜偏态系數、譜峰态系數和譜功率，獲得銑刀磨損度預測的頻域特征。圖中(a)展示的是C1資料的X軸振動信号預處理後的資料經過快速傅裡葉變換得到的幅值譜圖，圖中(b)(c)(d)分别是譜偏态系數、譜峰态系數、譜功率與磨損度的關系圖。
4.3 時頻聯合域特征提取效果
在小波包分解過程中，小波函數基采用多貝西小波函數，采用分解層數的小波包，則最後一層小波包數量為，下圖 (b)所示的是C1資料集X軸振動信号進行小波包分解提取的能量特征與磨損度關系的可視化結果，最後通過計算小波包分解結果的2-範數來提取能量特征，繪制出C1刀具磨損量與小波能量特征的關系如圖 (c)。
6 小結
至此，原始輸入資料的基于統計的時域特征、基于頻譜的頻域特征和基于小波能量的時頻域特征被提取出來，原來7個通道的資料擴充到70個通道，在更高維的數值空間表征了更多的過程資料資訊，為銑刀磨損度預測提供了豐富的樣本資料。
7 實驗代碼
import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport scipy.stats as stsimport timefrom pywt._wavelet_packets import WaveletPacketstart_time = time.time()# ****************************************# 特征提取函數# ****************************************def rms_fea(a): return np.sqrt(np.mean(np.square(a)))def var_fea(a): return np.var(a) # 方差def max_fea(a): return np.max(a)def skew_fea(a): # 偏度 return sts.skew(a)def kurt_fea(a): # 峰度 return sts.kurtosis(a)def pp_fea(a): # 峰峰值 return np.max(a) - np.min(a)def wave_fea(a): wp = WaveletPacket(a, 'db1', maxlevel=8) nodes = wp.get_level(8, 'freq') return np.linalg.norm(np.array([n.data for n in nodes]), 2)def spectral_skw(a): n = a.shape[0] mag = np.abs(np.fft.fft(a)) mag = mag[1:n//2]*2.00/n return sts.skew(mag)def spectral_kurt(a): n = a.shape[0] mag = np.abs(np.fft.fft(a)) mag = mag[1:n//2]*2.00/n return sts.kurtosis(mag)def spectral_pow(a): n = a.shape[0] mag = np.abs(np.fft.fft(a)) mag = mag[1:n//2]*2.00/n return np.mean(np.power(mag, 3))def extract_fea(data, num_stat=10): data_fea = [] dim_feature = 1 for i in range(dim_feature): data_slice = data data_fea.append(rms_fea(data_slice)) data_fea.append(var_fea(data_slice)) data_fea.append(max_fea(data_slice)) data_fea.append(pp_fea(data_slice)) data_fea.append(skew_fea(data_slice)) data_fea.append(kurt_fea(data_slice)) data_fea.append(wave_fea(data_slice)) data_fea.append(spectral_kurt(data_slice)) data_fea.append(spectral_skw(data_slice)) data_fea.append(spectral_pow(data_slice)) data_fea = np.array(data_fea) return data_fea.reshape((1, dim_feature*num_stat))# ****************************************# 資料處理 -->> 生成特征# ****************************************sample_num = 315sensor_num = 7feature_num = 70print()print('-' * 30)print(' 正在提取資料 ')print('-' * 30)print()# C1_normaldata = np.zeros((sample_num, feature_num), dtype=np.float32)for sample_index in range(0, sample_num): print("Processing C1 set:" + str(sample_index + 1)) for m in range(0, sensor_num): str1 = "%03d" % (sample_index + 1) str2 = m + 1 read_path = 'C:/Users/Kaifeng Guan/Desktop/My Tool Wear rediction/Code/new_data/dataset1/c1/c_1_' \ + str(str1) + '/f' + str(str2) + '.csv' pre_data = np.array(pd.read_csv(read_path)) pre_data = np.reshape(pre_data, (len(pre_data),)) data[sample_index, m*10:(m+1)*10] = extract_fea(pre_data)C1_normal = data# C4_normaldata = np.zeros((sample_num, feature_num), dtype=np.float32)for sample_index in range(0, sample_num): print("Processing C4 set:" + str(sample_index + 1)) for m in range(0, sensor_num): str1 = "%03d" % (sample_index + 1) str2 = m + 1 read_path = 'C:/Users/Kaifeng Guan/Desktop/My Tool Wear rediction/Code/new_data/dataset1/c4/c_4_' \ + str(str1) + '/f' + str(str2) + '.csv' pre_data = np.array(pd.read_csv(read_path)) pre_data = np.reshape(pre_data, (len(pre_data),)) data[sample_index, m*10:(m+1)*10] = extract_fea(pre_data)C4_normal = data# C6_normaldata = np.zeros((sample_num, feature_num), dtype=np.float32)for sample_index in range(0, sample_num): print("Processing C6 set:" + str(sample_index + 1)) for m in range(0, sensor_num): str1 = "%03d" % (sample_index + 1) str2 = m + 1 read_path = 'C:/Users/Kaifeng Guan/Desktop/My Tool Wear rediction/Code/new_data/dataset1/c6/c_6_' \ + str(str1) + '/f' + str(str2) + '.csv' pre_data = np.array(pd.read_csv(read_path)) pre_data = np.reshape(pre_data, (len(pre_data),)) data[sample_index, m*10:(m+1)*10] = extract_fea(pre_data)C6_normal = data
添加微信，和我一起讨論
掃一掃關注我們一起交流，共同進步
點亮
，是對我們最大的贊賞