读研之掉进故障检测（八）—辛辛那提大学IMS轴承数据集数据分析（带代码）

---

---

为什么从凯斯西储大学的轴承数据集转到了辛辛那提大学IMS轴承数据集呢？
因为IMS的数据集是run-to-failure的，能够反映一个轴承整个生命周期。
关于凯斯西储大学的数据处理可以看我的这篇文章读研之掉进故障检测（一）—CWRU轴承数据集的三种已建立的诊断技术

前言

本文将对IMS轴承数据集进行简要的说明，以及对其进行时域分析、频谱分析等

---

一、IMS轴承数据集介绍

关于数据集，我是用CSDN里面老哥提供的下载
数据集中一个有三个文件夹，分别是三份数据。具体介绍如下。

简要的说明是，文件统一是Ascll格式的数据。文件每10分钟记录一次数据，每个文件由20,480个点组成，采样率设置为20 kHz（但是在一篇具体介绍论文中指出了其实实际频率是20.48KHZ）。因此一个文件也就是采集了1s的数据。

• 数据集1有4个轴承，每个轴承分别有两个通道，在试验-失效实验结束时，轴承3出现内圈缺陷，轴承4出现滚子元件缺陷。
• 数据集2有4个轴承，每个轴承一个通道。在试验到失效实验结束时，轴承1发生了外圈故障。
• 数据集3有4个轴承，每个轴承一个通道。在试验到失效实验结束时，轴承3发生了外圈故障。
  同时这里要给出轴承的特征频率。
  

二、一些分析

1.时域分析

之前有看到一篇博文，博主在其中介绍了一种添加后缀再读取的方法，其实不需要这么复杂。ASCII格式文件是以ASCII编码表示的文本文件。默认情况下，pd.read_csv()函数将假定输入文件是以逗号作为字段分隔符的CSV文件。然而ASCII文件不是使用逗号分隔字段，在使用txt文件中可以看到数据使用空白间隔进行的。

那么间隔设置为sep=‘\s+’，下面直接上代码吧。由于数据集很大，所以我只取了一个轴承的两个通道。

import os
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
config = {
    "font.family": 'serif', # 衬线字体
    "font.size": 10, # 相当于小四大小
    "font.serif": ['SimSun'], # 宋体
    "mathtext.fontset": 'stix', # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大
    'axes.unicode_minus': False # 处理负号，即-号
}
rcParams.update(config)

folder_path = "E:/下载物品处理/be2023_12/IMS/1st_test"  # 替换为你的文件夹路径

# 获取文件夹中的所有文件
files = os.listdir(folder_path)

# 创建空的DataFrame用于存储拼接后的数据
combined_data = pd.DataFrame()

# 遍历每个文件
for file in files:
    file_path = os.path.join(folder_path, file)

    # 解析文件名
    date_parts = file.split('.')
    year = int(date_parts[0])
    month = int(date_parts[1])
    day = int(date_parts[2])
    hour = int(date_parts[3])
    minute = int(date_parts[4])
    second = int(date_parts[5])

    # 读取文件数据
    data = pd.read_csv(file_path, sep='\s+', header=None, usecols=[0,1], names=['Amplitude1', 'Amplitude2'])

    # 构建日期时间列
    datetime_column = pd.to_datetime(f"{year}-{month:02d}-{day:02d} {hour:02d}:{minute:02d}:{second:02d}")

    # 添加日期时间列到数据中
    data['DateTime'] = datetime_column

    # 将当前文件的数据合并到总体数据中
    combined_data = pd.concat([combined_data, data], ignore_index=True)

# 绘制时域图
plt.plot(combined_data['DateTime'], combined_data['Amplitude1'])
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('轴承1通道一数据')
plt.xticks(rotation=45)
plt.show()


# 绘制时域图
plt.plot(combined_data['DateTime'], combined_data['Amplitude2'])
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('轴承1通道二数据')
plt.xticks(rotation=45)
plt.show()

这里对结果做相关解释中间那些一根线的是因为在那段时间内没有数据，原因是折线图所以会有线跟过去。然后我选择的是数据集1轴承1，所以从时域上看不出问题。因为他本来就没问题。一般来说从时域来讲，如果有问题出现，它的幅度会变大。

2.频谱分析（做傅里叶变换）

在论文Analysis of the Rolling Element Bearing data set of the Center for Intelligent Maintenance Systems of the University of Cincinnati中，有这样的频谱分析但是我不是很理解这个图是怎么画出来的。

在思考几天之后，我发现可以直接使用plt.specgram()函数来画。
同时仔细观察图中的数据可以发现这个数据集时间长度在6-7天之间，那么能够满足的数据很显然是辛辛那提数据集2的，因此可以得到这样一张图。

SAMPLE_RATE = 20.48e3

# Calculate the spectrogram
plt.specgram(combined_data['Amplitude1'], Fs=SAMPLE_RATE, cmap='viridis')
plt.xlabel('Time')
plt.ylabel('Frequency (Hz)')
plt.title('Spectrogram - Bearing 1 Channel 1 Data')
plt.xticks(rotation=45)
plt.colorbar()  # Add a colorbar for the intensity scale
plt.show()

# Calculate the spectrogram
plt.specgram(combined_data['Amplitude2'], Fs=SAMPLE_RATE, cmap='viridis')
plt.xlabel('Time')
plt.ylabel('Frequency (Hz)')
plt.title('Spectrogram - Bearing 1 Channel 2 Data')
plt.xticks(rotation=45)
plt.colorbar()  # Add a colorbar for the intensity scale
plt.show()

得到结果如下

同时在此提供一个分析单文件频谱的方法（FT）做以分享。

import os
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
from scipy import fft
import numpy as np

config = {
"font.family": 'serif', # 衬线字体
"font.size": 10, # 相当于小四大小
"font.serif": ['SimSun'], # 宋体
"mathtext.fontset": 'stix', # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大
'axes.unicode_minus': False # 处理负号，即-号
}
rcParams.update(config)

SAMPLE_RATE = 20.48e3

file_path = "E:/下载物品处理/be2023_12/IMS/1st_test/2003.11.25.23.39.56"
data = pd.read_csv(file_path, sep='\s+', header=None, usecols=[6,7], names=['Amplitude1', 'Amplitude2'])

N = data.shape[0]


yf = fft.fft(data['Amplitude1'].values)

xf = fft.fftfreq(N, 1 / SAMPLE_RATE)

plt.plot(xf, np.abs(yf))
plt.xlabel('频率')
plt.ylabel('振幅')
plt.title('2003.11.25.23.39.56文件中轴承7通道1的fft')
plt.show()

data1 = pd.read_csv(file_path, sep='\s+', header=None, usecols=[2,3], names=['Amplitude1', 'Amplitude2'])

N1 = data1.shape[0]


yf1 = fft.fft(data1['Amplitude1'].values)

xf1 = fft.fftfreq(N, 1 / SAMPLE_RATE)

plt.plot(xf1, np.abs(yf1))
plt.xlabel('频率')
plt.ylabel('振幅')
plt.title('2003.11.25.23.39.56文件中轴承2通道1的fft')
plt.show()

3、倒谱预白化（算包络线）

import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
from scipy import fft, signal
import numpy as np

config = {
    "font.family": 'serif',  # 衬线字体
    "font.size": 10,  # 相当于小四大小
    "font.serif": ['SimSun'],  # 宋体
    "mathtext.fontset": 'stix',  # matplotlib渲染数学字体时使用的字体，和Times New Roman差别不大
    'axes.unicode_minus': False  # 处理负号，即-号
}
rcParams.update(config)

SAMPLE_RATE = 20.48e3

file_path = "E:/下载物品处理/be2023_12/IMS/1st_test/2003.11.25.23.39.56"
data = pd.read_csv(file_path, sep='\s+', header=None, usecols=[6, 7], names=['Amplitude1', 'Amplitude2'])

N = data.shape[0]

# ----去直流分量----#
data['Amplitude1'] = data['Amplitude1'] - data['Amplitude1'].mean()

yf = fft.fft(data['Amplitude1'].values)

xf = fft.fftfreq(N, 1 / SAMPLE_RATE)

# ----Cepstrum prewhitening----#
cepstrum = np.fft.ifft(np.log(np.abs(yf)))
cepstrum[:5] = 0

# ----做希尔伯特变换----#
xt = np.real(np.fft.fft(cepstrum))
ht = signal.hilbert(xt)
at = np.sqrt(xt ** 2 + ht ** 2)
am = np.fft.fft(at)
am = np.abs(am)
am = am / len(am) * 2
am = am[0: int(len(am) / 2)]
freq = fft.fftfreq(len(at), d=1 / SAMPLE_RATE)
freq = freq[0:int(len(freq) / 2)]

plt.plot(freq, am)
plt.xlabel('频率')
plt.ylabel('振幅')
plt.title('2003.11.25.23.39.56文件中轴承4通道1的倒谱预白化后的包络谱')
plt.show()

---

小总结

本人当时看那篇文章是，摘要说道：为了实现这一目标，使用了几种诊断技术：时间分析、谱分析、盲反卷积、谱相干、包络谱。
今天先写到这儿，下次把盲反卷积、谱相干、包络谱的分析补上