滚动轴承故障检测---基于希尔伯特变换解调（python实现）

一、基本原理

        请转向这里查看

二、方法实现

2.1 加载故障数据

        下面先以轴承外圈故障数据进行分析。原始数据为mat文件，是matlab文件，定义一个函数进行数据读取

from scipy.io import loadmat
import matplotlib.pyplot as plt
import math
import numpy as np
import scipy


plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False


def data_acquision(file_path):
    """
    fun: 从cwru mat文件读取加速度数据
    param file_path: mat文件绝对路径
    return accl_data: 加速度数据，array类型
    """
    file = loadmat(file_path)  # 加载mat数据
    file_keys = file.keys()
    for key in file_keys:
        if 'DE' in key:
            files = file[key].ravel()
    return files


if __name__ == "__main__":
    # step1: 获取测试数据
    path = r'D:\360安全浏览器下载\轴承故障检测\数据源\西储大学\CWRU轴承数据\cwru_data\0HP\12k_Drive_End_OR007@6_0_130.mat'
    data = data_acquision(path) #获取原始数据
    plt.figure(figsize=(12,3))
    plt.title('原始数据')
    plt.plot(data)
    print(data)

2.2 包络谱分析

    # step2: 截取数据 符合实际检测要求
    xt = data[:2000]
    plt.figure(figsize=(12,3))
    plt.title('截取数据')
    plt.plot(xt)
    print(xt)
    
    # step3: 做希尔伯特变换
    ht = scipy.fftpack.hilbert(xt)
    plt.figure(figsize=(12,3))
    plt.title('希尔伯特变换')
    plt.plot(ht)
    print(ht)
    
    # step4: 获取包络信号  
    # z(t)=x(t)+jh(t)=a(t)ejφt
    # 对z(t)取模，得到其幅值a(t). a(t)即为包络信号，也叫解析信号
    at = np.sqrt(ht**2+xt**2)   # at = sqrt(xt^2 + ht^2)  实部
    plt.figure(figsize=(12,3))
    plt.title('包络信号')
    plt.plot(at)
    print(at)
    
    # step5: 获取包络谱 
    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.title('包络谱')
    plt.plot(freq, am)
    plt.xlim(0,500)

2.3 特征频率计算

        2.3.1 轴承故障的特征频率

                内圈故障特征频率：

                外圈故障特征频率：

                滚动体故障特征频率：

                保持架故障特涨频率：

                其中， : 滚动体个数；轴转速；: 滚珠(子)直径；: 轴承节径

        2.3.2 轴承故障特征频率计算函数

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
        滚珠数量
    fr: float(r/min)
        转速
    d: float(mm)
        滚珠直径
    D: float(mm)
        轴承直径
    alpha: float(°)
        接触角
    fr:：float(r/min)
        转速
    Returns
    -------
    BPFI: float(Hz)
        内圈故障频率
    BPFO: float(Hz)
       外圈故障频率
    BSF: float(Hz)
        球体故障频率
    FTF: float(Hz)
        保持架频率
    '''
    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

        2.3.3 轴承特征频率计算

        轴承型号为：6205-2RSL JME SKF 深沟球滚珠轴承
        

        下面计算CWRU在转速1797rpm时的故障特征频率

    # step 6 
    bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1797)
    print('内圈故障特征频率',bpfi)
    print('外圈故障特征频率',bpfo)
    print('滚动体故障特征频率',bsf)
    print('保持架故障特征频率',ftf)
    print('当前转速',fr)

2.4 理论与实际验证

    plt.vlines(x=bpfo, ymin=0, ymax=0.2, colors='r')  # 一倍频
    plt.vlines(x=bpfo*2, ymin=0, ymax=0.2, colors='r')  # 二倍频