等时间采样无法对变转速设备进行故障分析,由于转速的变化,导致在等时间采样的情况下,无法采集到同相位的点,所以直接进行FFT转换后,信号往往会出现非常大的泄露,如图1所示。
图1 等时间采样频谱
阶次分析首先要确定参考轴的旋转频率。假设参考轴的转速为
实现阶次谱存在两个关键步骤:
(1).插值获取等角度间隔的时间点
(2).使用等角度间隔的时间点对振动信号重采样得到等角度间隔的振动信号
最后只需对重采样的振动信号FFT变换,即可得到阶次谱。需要注意角采样频率
采用渥太华升速轴承内圈损伤数据,数据为.mat数据,包含两类数据振动信号(Channel_1)和脉冲信号(Channel_2),采样频率200000,每圈脉冲数1024,内圈故障阶次5.43。
读取原始振动数据与转速脉冲(局部放大),如图2、图3。
图2 原始数据
图3 转速脉冲
2.根据转速脉冲计算转频时间曲线如图4,对给定时间内求取旋转角度的累积变化量如图5。
图4 转频时间曲线
图5 角度时间曲线
3.使用上述获取到的角度以及时间创建一个插值函数(python中使用interp1d实现),将角度进行等间隔划分并带入,即可获取到等角度间隔的时间点。使用原始振动数据和时间再次创建一个插值函数,将等角度间隔的时间点带入,即可获取到等角度间隔采样的振动信号,如图6。
图6 重采样振动信号
4.对重采样的振动信号进行FFT转换,得到阶次谱如图7(进行局部放大),但是故障阶次不够明显,进行包络解调得到包络阶次谱如图8(进行局部放大),此时故障阶次5.43非常明显。
图7 阶次谱
图8 阶次包络谱
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
from scipy.fftpack import fft
from scipy.integrate import cumulative_trapezoid, trapezoid
from scipy.fftpack import hilbert, fft, ifft
from matplotlib.pylab import mpl
from scipy import fftpack
import pandas as pd
import scipy
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
def trend_remove(x, y):
# 用2次多项式拟合x,y数组,返回多项式系数
param = np.polyfit(x, y, 2)
# 拟合完之后,用生成的多项式系数用来生成多项式函数
y_poly_func = np.poly1d(param)
# 生成多项式函数之后,就是获取x在这个多项式处的值
y_poly = y_poly_func(x)
# 原始信号减去拟合函数的值即为去除趋势项后的信号
de_trend_y = y - y_poly
return de_trend_y
file = r'C:\Users\Deng\Desktop\新建文件夹 (2)\I-A-1.mat'
dataset = scipy.io.loadmat(file)
vib = np.array([i[0] for i in dataset['Channel_1']])
pulse = np.array([i[0] for i in dataset['Channel_2']])
fs = 200000
t = (1 / fs) * np.arange(0, len(vib))
# 振动信号去除趋势项
vib = trend_remove(t, vib)
# 转速计算:SpeedFre
PulseTime = []
for i in range(len(pulse)-1):
if pulse[i + 1] < 2 < pulse[i]:
PulseTime.append(t[i])
PulseTime = np.array(PulseTime)
SpeedFre = []
speedTime = []
for i in range(0, len(PulseTime) + 1, 1024):
if i == 0:
SpeedFre.append(0)
else:
Fre = 1 / (PulseTime[i] - PulseTime[i - 1024])
SpeedFre.append(Fre)
speedTime.append(PulseTime[i])
SpeedFre = np.array(SpeedFre)
SpeedFre = SpeedFre[1:]
PulseTime = PulseTime[1:]
# 积分求取转过的总角度Angle
Angle = cumulative_trapezoid(SpeedFre*360, speedTime, initial=0)
plt.figure(figsize=(11, 4))
plt.plot(speedTime, SpeedFre)
plt.title('转速')
plt.xlabel('时间')
plt.ylabel('hz')
plt.figure(figsize=(11, 4))
plt.plot(speedTime, SpeedFre*360)
plt.title('角度')
plt.xlabel('时间')
plt.tight_layout()
plt.show()
# 获取等角度间隔的时间点:EqualAngleTime
angleStep = 0.1
EqualAngle = np.arange(0, trapezoid(SpeedFre*360, speedTime), angleStep)
FunctionA = interp1d(Angle, speedTime)
EqualAngleTime = FunctionA(EqualAngle)
# 用等角度间隔的时间点对振动信号重新采样,得到了等角度间隔的振动信号:EqualAngleVib
FunctionV = interp1d(t, vib)
EqualAngleVib = FunctionV(EqualAngleTime)
# 对重采样信后进行fft变换
VibFft = (np.abs(fft(EqualAngleVib)) * 2) / len(EqualAngleVib)
VibFftHalf = VibFft[1:int(np.round(len(EqualAngleVib) / 2))]
order = np.arange(0, len(EqualAngleVib), dtype='int64') * (1/angleStep*360) / len(EqualAngleVib)
orderHalf = order[1:int(np.round(len(EqualAngleVib) / 2))]
plt.figure(figsize=(11, 4))
plt.plot(EqualAngle, EqualAngleVib)
plt.title('重采样之后的振动信号')
plt.xlabel('角度')
plt.ylabel('振动幅值')
plt.figure(figsize=(11, 4))
plt.plot(orderHalf, VibFftHalf)
plt.title('阶次谱')
plt.xlabel('阶次')
plt.ylabel('振动幅值')
plt.xlim(0, 100)
plt.show()
angleFs = (1/angleStep)*360
orderMin = 1
orderMax = 500
# 带通滤波后逆变换
signal_fft = fftpack.fft(EqualAngleVib)
sig_fre = fftpack.fftfreq(n=len(signal_fft), d=1 / angleFs)
signal_fft[(abs(sig_fre) < orderMin) | (abs(sig_fre) > orderMax)] = 0
signal_filter = fftpack.ifft(signal_fft).real
# 包络变换
H = np.abs(hilbert(signal_filter) - np.mean(signal_filter))
HP = np.abs(fft(H - np.mean(H))) * 2 / len(EqualAngleVib)
HP_half = HP[1:int(np.round(len(EqualAngleVib) / 2))]
plt.figure(figsize=(11, 4))
plt.plot(orderHalf, HP_half)
plt.title('阶次包络谱')
plt.xlabel('阶次')
plt.ylabel('振动幅值')
plt.xlim(0, 100)
plt.show()
校核:李正平、陈凯歌、赵栓栓、曹希铭、赵学功、白亮、任超、陈宇航、Tina、王金、陈莹洁
该文资料搜集自网络,仅用作学术分享,不做商业用途,若侵权,后台联系小编进行删除