工業數據:基於統計和信號的特徵提取方法

在機械加工過程中，刀具銑削過程傳感器監控數據均為時序數據，如下圖所示，表現為一系列沒有固定規律的隨機離散值，不能有效識別數據特徵與磨損度的關係，因此需要進行特徵提取，把原有的低維時序數據映射到高維表示，挖掘和構造輸入數據的新特徵，進而找出最能幫助模型擬合的輸入(X)和輸出(Y)特徵。

在很多預測性維護的場景（如旋轉件、振動等）中，我們完全可以採用基於統計和信號的特徵提取方法就能得到一些較好的特徵。本文的特徵提取方法主要包含三個角度：

基於統計的時域特徵，就是使用統計學方法如均方根、方差、最大值、最小值、偏斜度、峰度、峰峰值等提取出的新數據。

1.1 均方根(Root Mean Square, RMS)

均方根(Root Mean Square, RMS)，是信號有效值的反映：

方差(Variance)，衡量隨機變量或一組數據離散程度，是源數據和期望值相差的度量：

最大值(Max)和最小值(Min)，表示在一定時間範圍內數據的最強和最弱的程度：

1.4 偏度(Skewness)

偏度(Skewness)又稱偏態係數，度量數據的概率密度曲線分布偏斜的方向和程度，即與平均值相比數據的非對稱的程度特徵。

其中，

1.5 峰峰值(Peak to Peak, PP)

峰峰值(Peak to Peak, PP)，即一個周期內數據最高值和最低值之間的差值，描述了數據的變化範圍的大小：

1.6 峰度(Kurtosis)

峰度(Kurtosis)又稱峰態係數，表示在平均值處數據的概率密度分布曲線的最大值大小。峰度直觀地反映了曲線的峰部尖度，峰度越高，說明曲線底部厚度越大，會導致方差越大，原因是數據存在一定數量的顯著比平均值大或小的差值。

其中，

頻域分析是信號領域中很重要的分析方法，具體包括了頻譜分析、能量譜分析、功率譜分析和倒頻譜分析等，其中以頻譜分析最為常用也最為重要。頻譜分析主要將信號數據通過傅立葉變換得到幅頻譜和相頻譜進行分析，有關傅立葉變換的內容網上有很多，就不贅述了。

通過計算機和傳感器採集得到的信號是離散信號，可以採用離散傅立葉變換：

信號的幅值

分別以幅值

2.1 譜偏態係數

譜偏態係數(Spectral Skewness)：統計頻域上幅值的分布偏斜的方向和程度：

2.2 譜峰態係數

譜峰態係數(Spectral Kurtosis)：表示頻域上幅值波形的尖銳程度：

2.3 譜功率

譜功率(Spectral Power)：表示單位頻帶內的信號功率：

小波變換（wavelet transform，WT）是一種新的變換分析方法，它繼承和發展了短時傅立葉變換局部化的思想，同時又克服了窗口大小不隨頻率變化等缺點，能夠提供一個隨頻率改變的「時間-頻率」窗口，是進行信號時頻分析和處理的理想工具。

小波變換過程簡述如下：

原始輸入信號X通過一組小波函數基經低通濾波器低頻分量（又稱為近似子帶）和高頻分量（又稱為細節子帶），將低頻分量作為輸入信號，又進行小波分解，得到下一層的低頻分量和高頻分量，以此類推，可以得到

小波包分解(Wavelet Packet Decomposition)又稱為最優子帶樹結構(Optimal Subband Tree Structuring)，如圖(b)所示，是對小波變換的進一步優化，具體為：在小波變換的基礎上，也對細節子帶進一步分解，最後通過最小化代價函數，計算出最優的信號分解路徑，並以此路徑對原始輸入信號進行分解。

選擇合適的小波函數基是進行信號分解的關鍵步驟。常用的小波函數基有哈爾小波函數基、多貝西小波和雙正交小波等。

多貝西小波作為稀疏基引入信號的光滑誤差不容易被察覺，因此具有較好的正則性，信號重構過程比較光滑，而且可以通過階次N的來控制頻域的局部化能力，使用較為靈活，因此本文的小波包分解過程採用多貝西小波函數基。

介紹完枯燥的原理後，我們通過實際操作來看看這麼常用的特徵提取方法最終的效果如何。

4.1 時域特徵提取效果

圖中展示了PHM2010刀具磨損預測數據集C1數據的X軸的振動信號在銑刀磨損周期內的6個統計學特徵與磨損度的關係圖。

4.2 頻域特徵提取效果

頻域特徵首先對預處理後的原始數據進行FFT處理，得到幅值譜圖，然後分別計算譜偏態係數、譜峰態係數和譜功率，獲得銑刀磨損度預測的頻域特徵。圖中(a)展示的是C1數據的X軸振動信號預處理後的數據經過快速傅立葉變換得到的幅值譜圖，圖中(b)(c)(d)分別是譜偏態係數、譜峰態係數、譜功率與磨損度的關係圖。

4.3 時頻聯合域特徵提取效果

在小波包分解過程中，小波函數基採用多貝西小波函數，採用分解層數

至此，原始輸入數據的基於統計的時域特徵、基於頻譜的頻域特徵和基於小波能量的時頻域特徵被提取出來，原來7個通道的數據擴展到70個通道，在更高維的數值空間表徵了更多的過程數據信息，為銑刀磨損度預測提供了豐富的樣本數據。

import pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport scipy.stats as stsimport timefrom pywt._wavelet_packets import WaveletPacket
start_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 = 70
print()print('-' * 30)print('         正在提取數據     ')print('-' * 30)print()
data = 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

data = 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

data = 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


點亮 ，是對我們最大的讚賞