Pyhon在振动信号处理中的高级应用（一）：振动信号特征（时域、频域）提取

• 一、概述
• 二、算法原理
• • 2.1 时域统计特征
  • 2.2 频域统计特征
• 三、python代码实现
• 四. 展示效果
• 五. tips

  振动信号故障诊断经历了多年的发展，已从老专家式的“线谱分析”逐步过渡到大数据故障诊断。巧妙地提取信号特征，利用各种智能算法对其进行学习，最终获取可靠的故障诊断模型，这一套流程成为了近些年故障诊断相关文章的标准格式。那么振动信号有哪些特征可供我们挑选呢？让我们一起来看看！

  这里我们借鉴西安交通大学机械学院雷教授论文《Fault diagnosis of rotating machinery based on multiple ANFIS combination with GAs》里的相关研究成果。雷教授在论文中共列举了2类共24种振动信号的时域和频域特征。

  第一类为时域统计特征，可直接通过振动信号的时间序列数据计算得出，如下列公式所示，其中： P 1 、 P 3 — P 5 P_{1}、P_{3}—P_{5} P1​、P3​—P5​反映了振动信号的时域幅值和能量特性； P 2 、 P 6 — P 11 P_{2}、P_{6}—P_{11} P2​、P6​—P11​反映了振动信号的时域分布特性。
P 1 = ∑ n = 1 N x ( n ) N P 2 = ∑ n = 1 N ( x ( n ) − P 1 ) 2 N − 1 P 3 = ( ∑ n = 1 N ∣ x ( n ) ∣ N ) 2 P_{1}=frac{sum_{n=1}^{N}x(n)}{N}quad P_{2}=sqrt{frac{sum_{n=1}^{N}left(x(n)-P_{1}right)^{2}}{N-1}}quad P_{3}=left(frac{sum_{n=1}^{N} sqrt{|x(n)|}}{N}right)^{2} P1​=N∑n=1N​x(n)​P2​=N−1∑n=1N​(x(n)−P1​)2​ ​P3​=(N∑n=1N​∣x(n)∣ ​​)2
P 4 = ∑ n = 1 N ( x ( n ) ) 2 N P 5 = max ⁡ ∣ x ( n ) ∣ P 6 = ∑ n = 1 N ( x ( n ) − P 1 ) 3 ( N − 1 ) P 2 3 P_{4}=sqrt{frac{sum_{n=1}^{N}(x(n))^{2}}{N}}quad P_{5}=max |x(n)|quad P_{6}=frac{sum_{n=1}^{N}left(x(n)-P_{1}right)^{3}}{(N-1) P_{2}^{3}} P4​=N∑n=1N​(x(n))2​ ​P5​=max∣x(n)∣P6​=(N−1)P23​∑n=1N​(x(n)−P1​)3​
P 7 = ∑ n = 1 N ( x ( n ) − P 1 ) 4 ( N − 1 ) P 2 4 P 8 = P 5 P 4 P 9 = P 5 P 3 P_{7}=frac{sum_{n=1}^{N}left(x(n)-P_{1}right)^{4}}{(N-1) P_{2}^{4}}quad P_{8}=frac{P_{5}}{P_{4}}quad P_{9}=frac{P_{5}}{P_{3}} P7​=(N−1)P24​∑n=1N​(x(n)−P1​)4​P8​=P4​P5​​P9​=P3​P5​​
P 10 = P 4 1 N ∑ n = 1 N ∣ x ( n ) ∣ P 11 = P 5 1 N ∑ n = 1 N ∣ x ( n ) ∣ P_{10}=frac{P_{4}}{frac{1}{N} sum_{n=1}^{N}|x(n)|}quad P_{11}=frac{P_{5}}{frac{1}{N} sum_{n=1}^{N}|x(n)|} P10​=N1​∑n=1N​∣x(n)∣P4​​P11​=N1​∑n=1N​∣x(n)∣P5​​

  第二类为频域统计特征，是对振动信号频域数据进行计算得到的，如下列公式所示，其中： P 12 P_{12} P12​反映了振动信号的频域幅值和能量特性； P 13 — P 15 、 P 17 、 P 21 — P 24 P_{13}—P_{15}、P_{17}、P_{21}—P_{24} P13​—P15​、P17​、P21​—P24​反映了振动信号的频域分布特性； P 16 、 P 18 — P 20 P_{16}、P_{18}—P_{20} P16​、P18​—P20​反映了振动信号的主要频率峰值。
P 12 = ∑ k = 1 K s ( k ) K P 13 = ∑ k = 1 K ( s ( k ) − P 12 ) 2 K − 1 P 14 = ∑ k = 1 K ( s ( k ) − P 12 ) 3 K ( P 13 ) 3 P_{12}=frac{sum_{k=1}^{K} s(k)}{K} quad P_{13}=frac{sum_{k=1}^{K}left(s(k)-P_{12}right)^{2}}{K-1} quad P_{14}=frac{sum_{k=1}^{K}left(s(k)-P_{12}right)^{3}}{Kleft(sqrt{P_{13}}right)^{3}} P12​=K∑k=1K​s(k)​P13​=K−1∑k=1K​(s(k)−P12​)2​P14​=K(P13​ ​)3∑k=1K​(s(k)−P12​)3​
P 15 = ∑ k = 1 K ( s ( k ) − P 12 ) 4 K P 13 2 P 16 = ∑ k = 1 K f k s ( k ) ∑ k = 1 K s ( k ) P 17 = ∑ k = 1 K ( f k − P 16 ) 2 s ( k ) K P_{15}=frac{sum_{k=1}^{K}left(s(k)-P_{12}right)^{4}}{K P_{13}^{2}} quad P_{16}=frac{sum_{k=1}^{K} f_{k} s(k)}{sum_{k=1}^{K} s(k)} quad P_{17}=sqrt{frac{sum_{k=1}^{K}left(f_{k}-P_{16}right)^{2} s(k)}{K}} P15​=KP132​∑k=1K​(s(k)−P12​)4​P16​=∑k=1K​s(k)∑k=1K​fk​s(k)​P17​=K∑k=1K​(fk​−P16​)2s(k)​ ​
P 18 = ∑ k = 1 K f k 2 s ( k ) ∑ k = 1 K s ( k ) P 19 = ∑ k = 1 K f k 4 s ( k ) ∑ k = 1 K f k 2 s ( k ) P 20 = ∑ k = 1 2 f k 2 s ( k ) ∑ k = 1 K f k 2 s ( k ) P_{18}=sqrt{frac{sum_{k=1}^{K} f_{k}^{2} s(k)}{sum_{k=1}^{K} s(k)}} quad P_{19}=sqrt{frac{sum_{k=1}^{K} f_{k}^{4} s(k)}{sum_{k=1}^{K} f_{k}^{2} s(k)}}quad P_{20}=frac{sum_{k=1}^{2} f_{k}^{2} s(k)}{sum_{k=1}^{K} f_{k}^{2} s(k)} P18​=∑k=1K​s(k)∑k=1K​fk2​s(k)​ ​P19​=∑k=1K​fk2​s(k)∑k=1K​fk4​s(k)​ ​P20​=∑k=1K​fk2​s(k)∑k=12​fk2​s(k)​
P 21 = P 17 P 16 P 22 = ∑ k = 1 K ( f k − P 16 ) 3 s ( k ) K P 17 3 P 23 = ∑ k = 1 K ( f k − P 16 ) 4 s ( k ) K P 17 4 P_{21}=frac{P_{17}}{P_{16}} quad P_{22}=frac{sum_{k=1}^{K}left(f_{k}-P_{16}right)^{3} s(k)}{K P_{17}^{3}} quad P_{23}=frac{sum_{k=1}^{K}left(f_{k}-P_{16}right)^{4} s(k)}{K P_{17}^{4}} P21​=P16​P17​​P22​=KP173​∑k=1K​(fk​−P16​)3s(k)​P23​=KP174​∑k=1K​(fk​−P16​)4s(k)​
P 24 = ∑ k = 1 K ( f k − P 16 ) 1 / 2 s ( k ) K P 17 P_{24}=frac{sum_{k=1}^{K}left(f_{k}-P_{16}right)^{1 / 2} s(k)}{K sqrt{P_{17}}} P24​=KP17​ ​∑k=1K​(fk​−P16​)1/2s(k)​

import numpy as np
from scipy import stats
from scipy.signal import welch

# signal为一维信号
def get_features(signal, fs):
    # 一、获取时域相关特征
    # 1.均值(实际振动数据均值常常为0，为更广泛使用，采用绝对值进行计算)
    p_1 = np.mean(abs(signal))
    
    # 2.标准差
    p_2 = np.std(signal)
    
    # 3.方根幅值
    p_3 = ((np.mean(np.sqrt(np.abs(signal)))))**2
    
    # 4.均方根值
    p_4 = np.sqrt((np.mean(signal**2)))
    
    # 5.最大峰值
    p_5 = np.max(abs(signal))
    
    # 6.偏度
    p_6 = stats.skew(signal)
    
    # 7.峭度
    p_7 = stats.kurtosis(signal)
    
    # 8.峰值因子
    p_8 = p_5/p_4
    
    # 9.裕度因子
    p_9  = p_5/p_3  
    
    # 10.波形因子
    p_10 = p_4/p_1
    
    # 11.脉冲因子
    p_11 = p_5/p_1

    # 二、获取频域相关序特征
    # 求取频谱图（1Hz分辨率）
    fre_sig, mag_sig = welch(sig, fs, nperseg=fs, noverlap = fs // 2, scaling='spectrum')
    mag_sig = np.sqrt(mag_sig)*np.sqrt(2)
    k = len(mag_sig)
    
    p_12 = np.mean(mag_sig)

    p_13 = np.var(mag_sig)

    p_14 = (np.sum((mag_sig-p_12)**3))/(k*((np.sqrt(p_13))**3))

    p_15 = (np.sum((mag_sig-p_12)**4))/(k*((p_13)**2))

    p_16 = (np.sum(fre_sig*mag_sig))/(np.sum(mag_sig))

    p_17 = np.sqrt((np.mean(((fre_sig-p_16)**2)*(mag_sig))))

    p_18 = np.sqrt((np.sum((fre_sig**2)*mag_sig))/(np.sum(mag_sig)))

    p_19 = np.sqrt((np.sum((fre_sig**4)*mag_sig))/(np.sum((fre_sig**2)*mag_sig)))

    p_20 = (np.sum((fre_sig**2)*mag_sig))/(np.sqrt((np.sum(mag_sig))*(np.sum((fre_sig**4)*mag_sig))))

    p_21 = p_17/p_16

    p_22 = (np.sum(((fre_sig - p_16)**3)*mag_sig))/(k*(p_17**3))

    p_23 = (np.sum(((fre_sig - p_16)**4)*mag_sig))/(k*(p_17**4))
    
    p_24 = (np.sum((np.sqrt(abs(fre_sig - p_16)))*mag_sig))/(k*np.sqrt(p_17))
    
    p = np.array([p_1, p_2, p_3, p_4, p_5, p_6, p_7, p_8, p_9, p_10, p_11, p_12, 
                p_13, p_14, p_15, p_16, p_17, p_18, p_19, p_20, p_21, p_22, p_23, p_24]) 
    return p

# 构造模拟信号
fre = 10
fs = 1000
t = np.arange(0, 1, 1/fs)
signal = np.sin(2*np.pi*fre*t) 

p = get_features(signal, fs)
print(p)

  笔者以振动信号特征的提取作为本系列文章的开头，一是想蹭一波大数据故障诊断的热度；第二呢，嘿嘿，挖坑。别以为提取了特征就能做好故障诊断，不信你试试，没有一个能打得。没有好的、特征明显的原始振动信号，休想玩好大数据故障诊断，接下来咱们一起老老实实研究如何提取特征明显的原始振动信号吧！EMD、小波、盲源分析吧啦吧啦，走起。