振动信号阶次分析

JoinApper

于 2025-02-07 15:10:13 发布

阅读量1.1k

点赞数 20

分类专栏：故障诊断信号与系统文章标签：机器学习信号处理 python

本文链接：https://blog.csdn.net/qq1234534215/article/details/145473754

版权

故障诊断同时被 2 个专栏收录

10 篇文章

订阅专栏

信号与系统

9 篇文章

订阅专栏

前言

之前介绍过轴承故障频率的计算方法，如果遇到转速变化缓慢的工况，我们可以直接使用FFT傅里叶变换方法计算其频谱，再将其中的故障频率能量提取出来，观察其一段时间内是否有持续上升趋势，从而判别是否出现故障。
但在实际工业环境下，大部分都是变转速的情况，由于FFT是对整个信号进行全局频率分析，因此会将不同时间段的频率成分混合在一起，导致频谱中出现频率混叠，无法准确反映真实的频率成分。
因此我们需要借助其他手段进行分析，常见的包括时频分析、阶次分析等，这篇文章就先来看看阶次分析是如何进行的。

基本原理

在介绍阶次分析的步骤之前，我们先要知道阶次这个概念。实际上阶次是用来描述振动频率与旋转机械转速之间关系的工具，就拿常见的轴承举例，轴承转速实际上就是一阶振动频率，而通过上篇文章的介绍，我们知道外圈损伤特征频率的计算公式如下：
$\begin{aligned} f_{ro}&=\frac{f_i-f_e}{2}(1-\frac{dcos(\alpha)}{D_m}) \end{aligned}$
其中 $f_i$ 是内圈旋转频率， $f_e$ 是外圈旋转频率，如果轴承是外圈固定、内圈旋转的方式，那么 $f_i$ 实际上就是轴承转速，而 $f_e$ 为0，那么外圈损伤特征频率的阶次为：
$\frac{f_{ro}}{f_i}=\frac{1}{2}(1-\frac{dcos(\alpha)}{D_m})$
那么为什么阶次可以用来做变转速情况下的分析呢？核心就在于虽然很多分析振动频率受到转速的影响，发生了变化，但是其相比于转速的倍数关系实际上是保持一致的，既阶次不变。因此阶次分析特别适用于变转速机械的故障诊断，它能够有效提取与转速相关的振动特征。

分析步骤

由上可知，只要我们将频率转换成阶次，就可以进行变转速下的分析工作了，那么如何才能转换呢？时域信号转换为频域后，每个频率代表的含义是一秒内完成周期性运动的次数，比如转频就是一秒内旋转的圈数，而阶次分析需要以转速作为基准，每个阶次代表的含义是旋转一圈完成周期性运动的次数，实际上就是要把信号从时域转换到角域。
角域采样理论告诉我们，对于变转速过程中的非稳态信号，若相对于旋转轴进行恒角度增量采样，则该时域非稳态信号在角度域是稳态信号，此时对角度域稳态信号进行FFT变换则可以得到清晰的图谱，即阶次谱。
阶次分析可以分为以下四个步骤：

采集转速信号：通过编码器或转速传感器获取转速脉冲信号（Tachometer Signal）。
计算角度信号：对转速信号进行积分，得到角度随时间的变化。
等角度重采样：根据角度信号将时域振动信号重采样为等角度间隔的信号。
角域分析：对重采样后的信号进行FFT分析，得到阶次谱。

当然目前也有很多研究集中在无转速脉冲信号下的转速提取工作上，后续再详细讲讲这个话题，此次就先假设我们有转速脉冲信号，下面从三个关键环节进行展开。

转速脉冲计算转速

我们知道转速脉冲来源于转速传感器，转速脉冲信号的频率与旋转机械的转速成正比，那么如何计算转速呢？实际上，我们只要知道旋转一周的脉冲数量N，以及两个脉冲之间的时间间隔t（脉冲信号频率f），就可以计算出转速了，具体的计算公式如下所示：
$RPM=\frac{60}{N*t}=\frac{60 * f}{N}$
下面我们用加拿大渥太华的数据集来演示一下，使用了H-A-1.mat数据，它是加速条件下健康轴承的振动数据。
在这里插入图片描述
原始数据总共采集了10s，为了方便演示，此次只选取了前100个数据点。为了计算出两个脉冲之间的时间间隔t，我们需要找到脉冲的上升沿或者下降沿位置。以上升沿为例，可以使用以下代码：

# 寻找上升沿
diff = np.diff((speed_data > np.mean(speed_data)).astype(int))
pulse = np.where(diff == 1)[0]

我们来看看其是否能够找出上升沿的位置：
在这里插入图片描述
基本上还是可以找到相应的位置，然后我们就可以使用之前所提公式在每个脉冲信号处计算瞬时转速了。

# 计算相邻脉冲之间的时间间隔，fs是采样时钟频率
t = np.diff(pulse) / fs
# 计算瞬时转速 (转/分钟)
speeds = 60.0 / (t * N)  # 转换为RPM
times = pulse[1:] / fs
# 转速滤波
window_size = 1000
smooth_speeds = signal.savgol_filter(speeds, window_size, 3)

在这里插入图片描述

转速计算角度

当我们计算好转速后，对之进行积分就可以得到角度数据了。但在这之前，我们还需要将脉冲转速进行插值，以便得到每个采样时刻的转速。

# 将 RPM 转换为 Hz
speeds_hz = smooth_speeds / 60.0

# 插值转速，获得完整时间序列上的连续转速
full_time = np.arange(len(speed_data)) / fs
interp_speeds = interp1d(times, speeds_hz, kind='linear', fill_value="extrapolate")
full_speeds = interp_speeds(full_time)

# 计算角度变化曲线 (累积角度，单位: rad)
full_angle = np.cumsum(full_speeds * 2 * np.pi / fs)

等角度重采样

通过上面的步骤，我们已经得到了采样时间与角度的关系，由于采样时间与振动信号幅值是一一对应的，因此角度与振动数据的关系也确定了。此时再通过插值的方式，我们就可以得到等角度重采样后的振动数据了。

# 重采样的时候，一圈设定采样points_per_rotation个数据
points_per_rotation = 360

# 创建均匀角度序列
n_rotations = int(full_angle[-1] / (2 * np.pi))
uniform_angle = np.linspace(0, n_rotations * 2 * np.pi, n_rotations * points_per_rotation)

# 重采样振动信号
vib_interp = interp1d(full_angle, vibration_data, kind='linear', bounds_error=False, fill_value=0)
resampled_vib = vib_interp(uniform_angle)

最后对重采样后的振动数据做FFT变换，就可以得到阶次谱了，下面我们还是用H-A-1.mat数据来进行演示。

案例演示

代码示例

我们定义三个函数，load_data负责加载mat数据：

def load_data(file_path):
    """加载振动和脉冲数据"""
    data = loadmat(file_path)
    vibration_data = data['Channel_1'].flatten()
    speed_data = data['Channel_2'].flatten()
    return vibration_data, speed_data

equal_angle_resample_vibration_signal用于对原始振动数据进行重采样。

def equal_angle_resample_vibration_signal(speed_data, vibration_data, N=1024, fs=200000, points_per_rotation=360):
    """
    等角度重采样数据
    :param speed_data: 转速脉冲数据
    :param vibration_data: 振动数据
    :param N: 旋转一圈的脉冲数
    :param fs: 采样率
    :param points_per_rotation: 重采样时设定每圈采样点数
    :return:
    """
    # 寻找上升沿
    diff = np.diff((speed_data > np.mean(speed_data)).astype(int))
    pulse = np.where(diff == 1)[0]

    # 计算相邻脉冲之间的时间间隔，fs是采样时钟频率
    t = np.diff(pulse) / fs
    # 计算瞬时转速 (转/分钟)
    speeds = 60.0 / (t * N)  # 转换为RPM
    times = pulse[1:] / fs
    # 转速滤波
    window_size = 1000
    smooth_speeds = signal.savgol_filter(speeds, window_size, 3)

    # 将 RPM 转换为 Hz
    speeds_hz = smooth_speeds / 60.0

    # 插值转速，获得完整时间序列上的连续转速
    full_time = np.arange(len(speed_data)) / fs
    interp_speeds = interp1d(times, speeds_hz, kind='linear', fill_value="extrapolate")
    full_speeds = interp_speeds(full_time)

    # 计算角度变化曲线 (累积角度，单位: rad)
    full_angle = np.cumsum(full_speeds * 2 * np.pi / fs)

    # 创建均匀角度序列
    n_rotations = int(full_angle[-1] / (2 * np.pi))
    uniform_angle = np.linspace(0, n_rotations * 2 * np.pi, n_rotations * points_per_rotation)

    # 重采样振动信号
    vib_interp = interp1d(full_angle, vibration_data, kind='linear', bounds_error=False, fill_value=0)
    resampled_vib = vib_interp(uniform_angle)
    return resampled_vib, uniform_angle

order_spectrum用于计算阶次谱。

def order_spectrum(resampled_data, points_per_rotation=360, max_order=50):
    """
    计算阶次谱，并限制阶次范围
    :param resampled_data: 重采样后的角域信号
    :param points_per_rotation: 重采样时设定每圈采样点数
    :param max_order: 最大阶次
    :return:
    """
    # 应用窗函数
    window = np.hanning(len(resampled_data))
    windowed_data = resampled_data * window

    # 计算FFT
    spectrum = np.fft.fft(windowed_data)
    spectrum = np.abs(spectrum)[:len(spectrum) // 2]  # 单边谱

    # 计算阶次轴，截取与 spectrum 一致的部分
    orders = np.fft.fftfreq(len(resampled_data), d=1 / points_per_rotation)[:len(spectrum)]

    # 限制阶次范围
    mask = (orders > 0) & (orders <= max_order)
    return orders[mask], spectrum[mask]