python实战故障诊断之CWRU数据集（四）：线性回归模型的应用_对于故障预测,基于线性回归的

1. 概述

  从本章节开始，我们正式开启机器学习在轴承故障诊断算法中的应用系列分析章节，首先我们介绍的是线性回归模型。
  本章节通过分析轴承振动数据得到了系列信号特征，将各特征与轴承故障尺寸进行相关分析，获取了相关性最好的一个特征，最后采用线性回归模型对该特征与故障尺寸之间的关系进行了分析，相关分析结果及过程可有效指导实际故障诊断及故障预测。

2. 线性回归模型相关原理

  线性回归模型是机器学习中最基本的模型，其原理及计算方法建议参照以下文章：

• 【机器学习】一文读懂线性回归、岭回归和Lasso回归

• 一篇文章搞懂：岭回归和Lasso回归

• 机器学习算法系列（五）- Lasso回归算法（Lasso Regression Algorithm）

  我们需要重点关注一下，岭回归会得到很多参数接近0的值，而Lasso回归会尽可能使非零参数最少。

3. 代码实现过程

  第一步：在相同转速下，分析不同内圈故障尺寸下原始信号时域和频域特征，以及包络信号频域特征，特征的计算方法详见Pyhon在振动信号处理中的高级应用（一）：振动信号特征（时域、频域）提取

from scipy.io import loadmat
import numpy as np
from scipy import signal, fftpack, stats
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.preprocessing import MinMaxScaler
from sklearn import pipeline, linear_model
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error

plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.size'] = '12'
plt.rcParams['font.sans-serif'] = ['SimHei']

# 从.mat文件中解析数据
def get_data(source):
    data_DE = []
    data_FE = []
    data_BA = []
    data_Speed = []
    
    
    for i, key in enumerate(source.keys()):
        if i == 3:
            data_DE = source[key].flatten()
        elif i == 4:
            data_FE = source[key].flatten()
        elif i == 5:
            if len(source.keys()) == 6:
                data_Speed = source[key].flatten()
            else:
                data_BA = source[key].flatten()
        elif i == 6:
            data_Speed = source[key].flatten()
                
    return data_DE, data_FE, data_BA, data_Speed

def envelope(sig, fs, low_pass):
    # 信号解包络
    sig_envelope = sig ** 2
    # 求取包络信号的FFT变换结果
    sig_fft = fftpack.fft(sig_envelope)
    # 求取包络信号的FFT变换频率
    sig_fre = fftpack.fftfreq(len(sig_fft), d=1/fs)
    # 低通滤波器
    sig_fft[abs(sig_fre)>low_pass] = 0
    # 转换为滤波后时域信号
    sig_filter = fftpack.ifft(sig_fft).real
    # 求取包络曲线的FFT结果
    fre_sig, mag_sig = signal.welch(sig_filter, fs=fs, nfft=len(sig_filter), 
                                    nperseg=len(sig_filter), scaling = 'spectrum')
    mag_sig = mag_sig[fre_sig<=low_pass]
    fre_sig = fre_sig[fre_sig<=low_pass]
    mag_sig = np.sqrt(2*mag_sig)
    return fre_sig, mag_sig

def cal_feature(source_path, source_list, types, speeds, fs, resample=False):
    feature_list = []
    type_list = []
    speed_list = []
    for _type, _speed, _list in zip(types, speeds, source_list):

        # 获取待处理数据
        source = loadmat(source_path + str(_list) + '.mat')
        # 解析数据
        data_DE, data_FE, data_BA, data_Speed = get_data(source)
        # 重采样
        if resample:
            data_DE = data_DE[::4]
        # 定义切片长度
        split_N = fs
        # 获取信号特征
        for i in range(len(data_DE)//split_N):
            # 切分数据
            data = data_DE[i*split_N:(i+1)*split_N]
            # 原始信号特征提取
            # 获取数据时域特征
            feature_time1 = get_time_features(data)
            # 获取数据频域特征
            fre_sig, mag_sig = signal.welch(data, fs=fs, nfft=split_N, 
                                            nperseg=split_N, scaling='spectrum')
            mag_sig = np.sqrt(2*mag_sig)
            feature_fre1 = get_frequence_features(fre_sig, mag_sig)
            # 获取包络谱频域特征
            fre_sig, mag_sig = envelope(data, fs, low_pass=200)
            feature_fre2 = get_frequence_features(fre_sig, mag_sig)
            # 特征组合
            feature = np.concatenate((feature_time1, feature_fre1, feature_fre2), axis=0)
            feature_list.append(list(feature))
            type_list.append(_type)
            speed_list.append(_speed)
    # 构建特征Dataframe
    columns = []
    for i in range(len(feature_time1)):
        columns.append('feature_time1-'+str(i+1))
    for i in range(len(feature_fre1)):
        columns.append('feature_fre1-'+str(i+12))
    for i in range(len(feature_fre2)):
        columns.append('feature_fre2-'+str(i+12))
    df = pd.DataFrame(feature_list, columns=columns)
    df['types'] = type_list
    df['speed'] = speed_list
    return df

# 获取正常数据特征
source_path = './/data/Normal Baseline Data/'
source_list = [97, 98, 99, 100]
types = [0, 0, 0, 0]
speeds = [1797, 1772, 1750, 1730]
fs = 12000

df_normal = cal_feature(source_path, source_list, types, speeds, fs, resample=True)

# 获取电机端内圈故障信号特征
source_path = './/data/12k Drive End Bearing Fault Data/'
source_list = [105, 106, 107, 108, 169, 170, 171, 172, 
               209, 210, 211, 212, 3001, 3002, 3003, 3004]
types = [0.18, 0.18, 0.18, 0.18, 0.36, 0.36, 0.36, 0.36,
         0.53, 0.53, 0.53, 0.53, 0.71, 0.71, 0.71, 0.71]
speeds = [1797, 1772, 1750, 1730, 1797, 1772, 1750, 1730,
          1797, 1772, 1750, 1730, 1797, 1772, 1750, 1730]
fs = 12000
resample = True

df_break = cal_feature(source_path, source_list, types, speeds, fs, resample=False)

# 整合数据结果
df = pd.concat([df_normal, df_break], ignore_index=True)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

  第二步：将各特征与故障尺寸之间进行相关分析，并取相关系数最大的一组特征，进行回归分析。

# 获取排名前十的相关系数
cof_corr = df[df['speed']==1772].corr()['types'].sort_values()[-10:-2]
print(cof_corr)

# 提取不同转速下的数据进行线性回归分析
# 构建线性回归分析模型
pipe = pipeline.Pipeline(steps=[('poly', PolynomialFeatures(degree=1)),
                                ('reg', linear_model.LinearRegression(fit_intercept=False))])
for speed in speeds[:4]:
    # 提取该转速下的数据
    data = df[df['speed']==speed][[cof_corr.index[-1], 'types']]
    # 对线性回归模型进行训练
    pipe.fit(data.values[:, 1].reshape(-1, 1), data.values[:, 0].reshape(-1, 1))
    x = np.arange(0, 1, 0.01).reshape(-1, 1)
    y = pipe.predict(x)
    [coef1, coef2] = pipe.named_steps['reg'].coef_[0]
    print('the coefficient is %.4f,%.4f'%(coef1, coef2))
    # 画图展示
    plt.figure(figsize=(10, 3))
    plt.scatter(data['types'], data.iloc[:, 0], c='blue', marker='+')
    plt.plot(x, y, c='r', linewidth=2)
    plt.title(str(speed) + 'rpm-' + cof_corr.index[-1])
    plt.xlabel('Fault  width/mm')
    plt.ylabel(cof_corr.index[-1])
    plt.xlim(0, 1)
    plt.ylim(0, 0.15)
    plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

图3.1 1730rpm-feature_fre2-21

图3.2 1750rpm-feature_fre2-21

图3.3 1772rpm-feature_fre2-21

图3.4 1797rpm-feature_fre2-21

4. Tips

  从结果中可以发现，在不同转速下，线性回归模型均表现了比较好的拟合特性，而岭回归和Lasso回归如何应用起来则需要大家根据自己的需求进行分析应用了。相信大家通过本章节内容可以得到一定的启发，咱们点到为止，防止各位小文章撞车，哈哈。