【“工业大数据预测”系列】——第3篇：基于统计和信号的特征提取

文章目录

1 基于统计的时域特征
- 1.1 均方根(Root Mean Square, RMS)
- 1.2 方差(Variance)
- 1.3 最大值(Max)和最小值(Min)
- 1.4 偏度(Skewness)
- 1.5 峰峰值(Peak to Peak, PP)
- 1.6 峰度(Kurtosis)
2 基于频谱分析的频域特征提取
- 2.1 谱偏态系数
- 2.2 谱峰态系数
- 2.3 谱功率
3 基于小波包能量的时频联合域特征
4 特征提取实际效果
- 4.1 时域特征提取效果
- 4.2 频域特征提取效果
- 4.3 时频联合域特征提取效果
5 小结
6 实验代码

机械加工过程中，刀具铣削过程传感器监控数据均为时序数据，如下图所示，为没有固定规律的随机离散值，不能有效识别数据特征与磨损度的关系，因此需要进行特征提取，把原有的低维时序数据映射到高维表示，挖掘和构造输入数据的新特征，进而找出最能帮助模型拟合的输入(X)和输出(Y)特征。

很多预测性维护的场景（如旋转件、振动等）我们完全可以采用基于统计和信号的特征提取方法就能得到一些较好的特征，本文的特征提取方法主要包含三个角度：

基于统计的时域特征
基于频谱分析的频域特征
基于小波包能量的时频联合域特征

1 基于统计的时域特征

基于统计的时域特征，就是使用统计学方法如均方根、方差、最大值、最小值、偏斜度、峰度、峰峰值等提取出的新数据。

1.1 均方根(Root Mean Square, RMS)

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

RMS=1n∑i=1nxi2R M S=\sqrt{\frac{1}{n} \sum_{i=1}^{n} x_{i}^{2}}RMS=n1i=1∑nxi2

1.2 方差(Variance)

方差(Variance)，衡量随机变量或一组数据离散程度，是源数据和期望值相差的度量：

Var⁡=1n∑i=1n(xi−xˉ)2\operatorname{Var}=\frac{1}{n} \sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2} Var=n1i=1∑n(xi−xˉ)2

1.3 最大值(Max)和最小值(Min)

最大值(Max)和最小值(Min)，表示在一定时间范围内数据的最强和最弱的程度：
Max⁡=max⁡(X)Min⁡=min⁡(X)\begin{array}{l} \operatorname{Max}=\max (X) \\ \operatorname{Min}=\min (X) \end{array} Max=max(X)Min=min(X)

1.4 偏度(Skewness)

偏度(Skewness)又称偏态系数，度量数据的概率密度曲线分布偏斜的方向和程度，即与平均值相比数据的非对称的程度特征。
Skew=E[(X−μδ)3]=EX3−3EXEX2+2E3X(EX2−E2X)3/2=μ3μ232=μ3δ3\text {Skew}=E\left[\left(\frac{X-\mu}{\delta}\right)^{3}\right]=\frac{E X^{3}-3 E X E X^{2}+2 E^{3} X}{\left(E X^{2}-E^{2} X\right)^{3 / 2}}=\frac{\mu_{3}}{\mu_{2}^{\frac{3}{2}}}=\frac{\mu_{3}}{\delta^{3}}Skew=E[(δX−μ)3]=(EX2−E2X)3/2EX3−3EXEX2+2E3X=μ223μ3=δ3μ3

其中， μ3\mu_3μ3是3阶中心距， δ\deltaδ为标准差。

1.5 峰峰值(Peak to Peak, PP)

峰峰值(Peak to Peak, PP)，即一个周期内数据最高值和最低值之间的差值，描述了数据的变化范围的大小：
PP=max⁡(X)−min⁡(X)P P=\max (X)-\min (X) PP=max(X)−min(X)

1.6 峰度(Kurtosis)

峰度(Kurtosis)又称峰态系数，表示在平均值处数据的概率密度分布曲线的最大值大小。峰度直观地反映了曲线的峰部尖度，峰度越高，说明曲线底部厚度越大，会导致方差越大，原因是数据存在一定数量的显著比平均值大或小的差值。
Kurt=μ4δ4=E[(X−μ)4](E[(X−μ)2])2\text {Kurt}=\frac{\mu_{4}}{\delta^{4}}=\frac{E\left[(X-\mu)^{4}\right]}{\left(E\left[(X-\mu)^{2}\right]\right)^{2}} Kurt=δ4μ4=(E[(X−μ)2])2E[(X−μ)4]
其中， μ4\mu_4μ4是四阶中心矩， δ4\delta^4δ4是标准差。

2 基于频谱分析的频域特征提取

频域分析是信号领域中很重要的分析方法，具体包括了频谱分析、能量谱分析、功率谱分析和倒频谱分析等，其中以频谱分析最为常用也最为重要。频谱分析主要将信号数据通过傅里叶变换得到幅频谱和相频谱进行分析，有关傅里叶变换的内容网上有很多，就不赘述了。

通过计算机和传感器采集得到的信号是离散信号，可以采用离散傅里叶变换：
F(ω)=∣F(ω)∣ejφ(ω)=Re⁡(F(ω))2+Im⁡(F(ω))2ejφ(ω)F(\omega)=|F(\omega)| e^{j \varphi(\omega)}=\sqrt{\operatorname{Re}(F(\omega))^{2}+\operatorname{Im}(F(\omega))^{2}} e^{j \varphi(\omega)} F(ω)=∣F(ω)∣ejφ(ω)=Re(F(ω))2+Im(F(ω))2ejφ(ω)
信号的幅值∣F(ω)∣|F(\omega)|∣F(ω)∣和相位φ(ω)\varphi(\omega)φ(ω)可以表示为：
∣F(ω)∣=Re⁡(F(ω))2+Im⁡(F(ω))2φ(ω)=arctan⁡(Im⁡(F(ω))Re⁡(F(ω)))\begin{array}{c} |F(\omega)|=\sqrt{\operatorname{Re}(F(\omega))^{2}+\operatorname{Im}(F(\omega))^{2}} \\ \varphi(\omega)=\arctan \left(\frac{\operatorname{Im}(F(\omega))}{\operatorname{Re}(F(\omega))}\right) \end{array} ∣F(ω)∣=Re(F(ω))2+Im(F(ω))2φ(ω)=arctan(Re(F(ω))Im(F(ω)))
分别以幅值∣F(ω)∣|F(\omega)|∣F(ω)∣和相位φ(ω)\varphi(\omega)φ(ω)为纵坐标，以频率作为横坐标，绘制出的二维坐标图为幅值谱图和相位谱图，幅值谱图表征信号的幅值随频率的分布情况，相位谱图表征了信号的相位随频率的变化情况。在铣刀磨损度预测任务中，可以通过幅值谱提取以下统计特征：

2.1 谱偏态系数

谱偏态系数(Spectral Skewness)：统计频域上幅值的分布偏斜的方向和程度：
Spec Skew=∑i=1k(Fi(ω)−Fˉ(ω)δ)3S(Fi(ω))\text {Spec Skew}=\sum_{i=1}^{k}\left(\frac{F_{i}(\omega)-\bar{F}(\omega)}{\delta}\right)^{3} S\left(F_{i}(\omega)\right) Spec Skew=i=1∑k(δFi(ω)−Fˉ(ω))3S(Fi(ω))

2.2 谱峰态系数

谱峰态系数(Spectral Kurtosis)：表示频域上幅值波形的尖锐程度：
Spec Kurt=∑i=1k(Fi(ω)−Fˉ(ω)δ)4S(Fi(ω))\text {Spec Kurt}=\sum_{i=1}^{k}\left(\frac{F_{i}(\omega)-\bar{F}(\omega)}{\delta}\right)^{4} S\left(F_{i}(\omega)\right) Spec Kurt=i=1∑k(δFi(ω)−Fˉ(ω))4S(Fi(ω))

2.3 谱功率

谱功率(Spectral Power)：表示单位频带内的信号功率：
Spec Pow =∑i=1k(Fi(ω))3S(Fi(ω))\text {Spec Pow }=\sum_{i=1}^{k}\left(F_{i}(\omega)\right)^{3} S\left(F_{i}(\omega)\right) Spec Pow =i=1∑k(Fi(ω))3S(Fi(ω))

3 基于小波包能量的时频联合域特征

小波变换（wavelet transform，WT）是一种新的变换分析方法，它继承和发展了短时傅立叶变换局部化的思想，同时又克服了窗口大小不随频率变化等缺点，能够提供一个随频率改变的“时间-频率”窗口，是进行信号时频分析和处理的理想工具。

小波变换过程简述如下：

原始输入信号X通过一组小波函数基经低通滤波器hkh_khk和高通滤波器gkg_kgk分解后进行下采样，得到低频分量（又称为近似子带）和高频分量（又称为细节子带），
将低频分量作为输入信号，又进行小波分解，得到下一层的低频分量和高频分量，以此类推，可以得到NNN层的小波分解，如图 (a)所示。
随着小波分解层数的增加，其在频域上的分辨率越高，这被称为多分辨率分析(MultiResolution Analysis, MRA)。

小波包分解(Wavelet Packet Decomposition)又称为最优子带树结构(Optimal Subband Tree Structuring)，如图(b)所示，是对小波变换的进一步优化，具体为：在小波变换的基础上，也对细节子带进一步分解，最后通过最小化代价函数，计算出最优的信号分解路径，并以此路径对原始输入信号进行分解。

选择合适的小波函数基是进行信号分解的关键步骤。常用的小波函数基有哈尔小波函数基、多贝西小波和双正交小波等。

多贝西小波作为稀疏基引入信号的光滑误差不容易被察觉，因此具有较好的正则性，信号重构过程比较光滑，而且可以通过阶次N的来控制频域的局部化能力，使用较为灵活，因此本文的小波包分解过程采用多贝西小波函数基。

4 特征提取实际效果

介绍完枯燥的原理后，我们通过实际操作来看看这么常用的特征提取方法最终的效果如何。

4.1 时域特征提取效果

图中展示了PHM2010刀具磨损预测数据集C1数据的X轴的振动信号在铣刀磨损周期内的6个统计学特征与磨损度的关系图。

4.2 频域特征提取效果

频域特征主要通过3.1.2节所述的基于频谱分析的方法来进行特征提取，首先对预处理后的原始数据进行FFT处理，得到幅值谱图，然后分别计算谱偏态系数、谱峰态系数和谱功率，获得铣刀磨损度预测的频域特征。
图中(a)展示的是C1数据的X轴振动信号预处理后的数据经过快速傅里叶变换得到的幅值谱图，图中(b)©(d)分别是谱偏态系数、谱峰态系数、谱功率与磨损度的关系图。

4.3 时频联合域特征提取效果

在小波包分解过程中，小波函数基采用多贝西小波函数，采用分解层数N=8N=8N=8的小波包，则最后一层小波包数量为2N=28=2562^N=2^8=2562N=28=256，下图 (b)所示的是C1数据集X轴振动信号进行小波包分解提取的能量特征与磨损度关系的可视化结果，最后通过计算小波包分解结果的2-范数来提取能量特征，绘制出C1刀具磨损量与小波能量特征的关系如图 ©。

5 小结

至此，原始输入数据的基于统计的时域特征、基于频谱的频域特征和基于小波能量的时频域特征被提取出来，原来7个通道的数据扩展到70个通道，在更高维的数值空间表征了更多的过程数据信息，为铣刀磨损度预测提供了丰富的样本数据。

6 实验代码

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sts
import time
from 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/nreturn 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/nreturn 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/nreturn np.mean(np.power(mag, 3))def extract_fea(data, num_stat=10):data_fea = []dim_feature = 1for i in range(dim_feature):data_slice = datadata_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 = 315
sensor_num = 7
feature_num = 70print()
print('-' * 30)
print('         正在提取数据     ')
print('-' * 30)
print()# C1_normal
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 + 1read_path = 'C:/Users/Kaifeng Guan/Desktop/My Tool Wear Prediction/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_normal
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 + 1read_path = 'C:/Users/Kaifeng Guan/Desktop/My Tool Wear Prediction/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_normal
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 + 1read_path = 'C:/Users/Kaifeng Guan/Desktop/My Tool Wear Prediction/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

如果您觉得文章对您有用，请给我点个赞吧！
您的肯定是对我最大的鼓励。