简介针对轴承全寿命数据分析与故障诊断需求这份资源提供了一套完整的时域与频域特征提取方案覆盖17个时域统计特征和13个频域统计特征适用于机械设备健康监测、剩余寿命预测等研究方向也适合初涉特征工程的科研人员对照学习。包内包含11个文件以9个MATLAB脚本.m为主辅以2个自动保存文件.asv压缩包整体仅12KB。脚本模块划分清晰时域统计与频域统计计算均有独立实现另封装了希尔伯特解调、ZOOMFFT细化频谱等分析功能可支撑从原始信号到特征向量的完整流程。目前已有1738人学习下载说明该工具脚本在轴承数据挖掘场景中具有不错的参考价值。借助这些脚本读者可快速复现常见时频特征提取逻辑理解统计指标与包络谱、ZOOMFFT等方法的实现细节并迁移到自己的数据集上进行二次开发。1. 轴承全寿命数据时域频域特征提取先用 30 个统计量把退化过程量化做轴承全寿命run-to-failure退化预测第一步通常不是建模而是把原始振动信号压缩成一组能长期跟踪的特征。17 个时域特征加 13 个频域特征是故障诊断里沿用多年的标准特征集覆盖幅值、能量、冲击、频谱形状四个维度。RMS 和能量负责整体退化趋势峭度、峰值因子、裕度因子负责早期点蚀冲击重心频率等频域矩反映频谱重心随磨损的迁移。没有哪个单特征能走完全寿命但 30 个组合起来退化阶段、报警阈值、剩余寿命回归就都有据可依。下面把公式、实现、参数边界和筛选用法一次讲清适合正在处理全寿命振动数据的工程师对照落地。2. 17 个时域特征公式清单、批量计算代码与无量纲因子的数值坑2.1 先对齐 17 个时域特征有量纲统计量与无量纲因子不同论文里的“17 个时域特征”成员并不完全一致按最常见的划分方式前 10 个是有量纲统计量后 7 个是无量纲因子和分布高阶矩。这里先给全清单后面代码和章节都按这个编号来。编号特征名公式对退化的响应T1均值 meanmean(x)静态偏置退化中变化小T2绝对平均值mean(|x|)幅值水平趋势与 RMS 接近T3方差 varvar(x)波动能量T4标准差 stdsqrt(var(x))波动幅度T5均方根 RMSsqrt(mean(x²))整体振动能量全寿命主趋势T6方根幅值(mean(sqrt(|x|)))²幅值水平对孤立冲击不敏感T7峰值 peakmax(|x|)冲击幅值T8峰峰值max(x) − min(x)最大摆幅T9最大值max(x)单向摆幅上限T10最小值min(x)单向摆幅下限T11偏度 skewnessE[(x−μ)³]/σ³分布对称性T12峭度 kurtosisE[(x−μ)⁴]/σ⁴早期冲击最敏感T13峰值因子peak / RMS冲击相对能量T14波形因子RMS / 绝对平均值波形形状变化T15脉冲因子peak / 绝对平均值冲击相对幅值T16裕度因子peak / 方根幅值冲击敏感比峰值因子更稳当T17能量sum(x²)总能量与 RMS 强相关表格里 x 表示窗口内的一维振动序列n 为窗口长度μ 和 σ 分别是均值和标准差。有量纲特征依赖传感器灵敏度和安装位置跨测点对比前要先归一化T13T16 四个无量纲因子消掉了幅值量纲在载荷波动工况下表现更一致这是它们能在早期故障里派上用场的原因之一。2.2 用 NumPy 和 SciPy 批量计算 17 个时域特征import numpy as np from scipy import stats def extract_td_features(x): 输入一段一维振动信号 x返回 17 个时域特征字典。 参数 x : np.ndarray长度建议 1024过短时偏度和峭度估计会失真。 x np.asarray(x, dtypenp.float64) rms np.sqrt(np.mean(x ** 2)) # T5 均方根 abs_mean np.mean(np.abs(x)) # T2 绝对平均值 peak np.max(np.abs(x)) # T7 峰值 sqrt_amp np.mean(np.sqrt(np.abs(x))) ** 2 # T6 方根幅值 return { td_mean: np.mean(x), td_abs_mean: abs_mean, td_var: np.var(x), td_std: np.std(x), td_rms: rms, td_sqrt_amp: sqrt_amp, td_peak: peak, td_pk_pk: np.max(x) - np.min(x), td_max: np.max(x), td_min: np.min(x), td_skew: stats.skew(x), # T11 偏度 td_kurt: stats.kurtosis(x, fisherFalse), # T12 峭度正态分布3 td_crest: peak / (rms 1e-12), # T13 峰值因子 td_shape: rms / (abs_mean 1e-12), # T14 波形因子 td_impulse: peak / (abs_mean 1e-12), # T15 脉冲因子 td_clearance: peak / (sqrt_amp 1e-12), # T16 裕度因子 td_energy: np.sum(x ** 2), # T17 能量 }几个参数和定义的取舍要讲清楚。np.var和np.std默认算总体值除以 n如果对标某些论文里的样本标准差除以 n−1小窗口下会差几个百分点全寿命特征序列里必须固定一种口径。stats.skew返回的是有偏估计stats.kurtosis(fisherFalse)返回 Pearson 峭度正态分布等于 3很多文献直接用超额峭度正态等于 0跨文章对比前必须先确认。T13T16 四个因子都要做分母保护RMS、绝对平均值、方根幅值在轴承刚启动或传感器静默时可能接近 0直接相除会出现极大的孤立尖峰1e-12的作用是把异常值限制在合理范围而不是掩盖真实冲击。2.3 短窗口下的失真和峭度回落的误读窗口长度小于 100 点时偏度和峭度的估计方差会急剧变大甚至返回 nan。全寿命切片时窗口至少要覆盖几百个采样点工程上建议 4096 点起步这不仅是为了统计可靠性也是为了下一章频域特征的分辨率。另一个常见误读是峭度曲线早期点蚀出现时峭度先冲高随后剥落面扩大、冲击变成连续宽带振动峭度反而回落——这不是退化结束而是故障形态从“稀疏冲击”变成了“密集宽带”。所以峭度要配 RMS 一起看一个看整体趋势一个看冲击形态。3. 13 个频域特征幅值谱约定、频率矩公式与 Python 落地3.1 频域特征统计的对象是幅值谱13 个频域特征本质上是对频谱再做一次统计。第一步要定死统计对象用幅值谱还是功率谱。幅值谱和时域幅值量纲一致边频带结构看得直观功率谱对幅值做了平方会放大主频带的权重。我一般固定用单边幅值谱 S(k) 定义全部 13 个特征P1P5 和 P10P13 都是对谱线的统计P6P9 是频率的一阶和二阶矩。常见论文里也有用功率谱算 P1P6 的版本数值不同但方向一致关键是同一个项目里只允许一种定义。口径统一后退化导致的频谱形态变化才能被这几个矩稳定捕捉。3.2 13 个频域特征的公式与物理含义编号特征公式含义P1幅值谱均值mean(S)谱线整体水平P2幅值谱方差var(S)谱线离散程度P3幅值谱标准差std(S)谱线离散程度P4幅值谱偏度skew(S)谱分布的非对称性P5幅值谱峭度kurt(S)谱峰尖锐程度P6重心频率 FCΣf·S / ΣS频谱重心位置HzP7均方频率 MSFΣf²·S / ΣS频率二阶矩Hz²P8频率方差 VFMSF − FC²谱线围绕重心的分散度Hz²P9均方根频率 RMSFsqrt(MSF)特征频率HzP10谱能量ΣS²频谱总能量P11峰值频率argmax(S) 对应 f主谱线频率P12峰值幅值max(S)主谱线高度P13谱总幅值ΣS幅值谱面积P6P9 四个频率矩对退化很敏感磨损导致宽带能量上升时高频段权重抬升P6 和 P8 的上升往往早于 RMS 的明显变化。P11 和 P12 直接给出主频成分的位置适合和转频、故障特征频率的理论值叠图对比。3.3 Python 实现加窗、单边谱与频率矩一条龙import numpy as np from scipy import stats def extract_fd_features(x, fs): 输入一段一维振动信号 x 和采样率 fs返回 13 个频域特征字典。 x np.asarray(x, dtypenp.float64) x x - np.mean(x) # 去直流避免 0 Hz 分量主导谱统计量 n len(x) win np.hanning(n) # Hann 窗抑制频谱泄露 xw x * win X np.fft.rfft(xw, nn) # 单边 FFT点数 n//21 amp np.abs(X) / np.sum(win) * 2.0 # 幅值谱补偿窗能量后单边乘 2 amp[0] amp[0] / 2.0 # 直流分量的 2 要退回去 freq np.fft.rfftfreq(n, d1.0 / fs) S amp denom np.sum(S) 1e-12 fc np.sum(freq * S) / denom # P6 重心频率 msf np.sum(freq ** 2 * S) / denom # P7 均方频率 vf max(msf - fc ** 2, 0.0) # P8 频率方差钳位到 0 peak_idx np.argmax(S) return { fd_amp_mean: np.mean(S), # P1 fd_amp_var: np.var(S), # P2 fd_amp_std: np.std(S), # P3 fd_amp_skew: stats.skew(S), # P4 fd_amp_kurt: stats.kurtosis(S, fisherFalse), # P5 fd_fc: fc, # P6 fd_msf: msf, # P7 fd_vf: vf, # P8 fd_rmsf: np.sqrt(msf), # P9 fd_spec_energy: np.sum(S ** 2), # P10 fd_peak_freq: freq[peak_idx], # P11 fd_peak_amp: S[peak_idx], # P12 fd_spec_sum: np.sum(S), # P13 }参数说明去均值这步不能省轴承信号里若有微小直流量0 Hz 谱线会直接抬升 P1 和 P13掩盖真实频谱形状。np.hanning的代价是主瓣展宽约一倍换来旁瓣明显压低适合做统计量但不适合精确测单根谱线幅值。幅值标定中np.sum(win)是窗函数能量补偿乘以 2 是把负频率能量并入正频直流分量不乘 2。P8 用max(..., 0.0)钳位因为浮点误差可能让 MSF − FC² 出现微小的负值没有物理意义。3.4 频率分辨率由窗口长度决定边频带是检验尺很多人习惯拿整段信号做一次 FFT但全寿命特征提取里谱是逐窗口算的频率分辨率 df fs / nn 是窗口长度而不是整段信号长度。以 25.6 kHz 采样率为例窗口 4096 点时 df 是 6.25 Hz窗口降到 1024 点df 变成 25 Hz。滚动轴承外圈故障的调制边带间隔等于转频约 33 Hz25 Hz 的分辨率会把相邻边带直接糊在一起P11 会在多个谱线间来回跳。窗口低于 2048 点时要警惕这类频域失真反过来窗口太长又会拖慢特征序列的时间分辨率报警变迟钝。40968192 点是全寿命数据常用的折中区间。4. 全寿命特征序列滑动窗口切片、特征矩阵构建与退化曲线4.1 窗口长度和步长按转频与退化速度两个约束定滑动窗口参数是特征提取里最容易被忽略的环节。约束一是窗口至少覆盖 510 个转频周期统计量才对相位不敏感约束二是步长决定特征序列的时间分辨率一般取窗口的 25%50%。对一组典型的全寿命数据按下面的参数起步基本不会跑偏。工况参数取值依据采样率 fs25600 Hz加速度传感器常见采集配置转频 fr约 33 Hz对应 2000 RPM 工况单圈采样点数约 77625600 / 33窗口长度4096约 5.3 圈统计上足够可靠步长1024相邻窗口 75% 重叠趋势平滑频率分辨率6.25 Hz25600 / 4096可分辨 33 Hz 边带有些全寿命数据集按固定间隔落盘比如每 10 分钟一个文件这时不需要滑动窗口直接对每个文件算 30 个特征就是一条特征序列滑动窗口用在连续采集的单一大文件上。两种做法的特征矩阵结构相同后续退化趋势处理也一致。4.2 批量切片并组装 30 维特征矩阵import pandas as pd def build_feature_frame(x, fs, win_len4096, step1024): 把整段全寿命信号切成滑动窗口返回 (窗口数, 30) 的特征 DataFrame。 rows [] n len(x) for start in range(0, n - win_len 1, step): seg x[start:start win_len] row {t_start: start / fs} # 窗口起始时间秒 row.update(extract_td_features(seg)) row.update(extract_fd_features(seg, fs)) rows.append(row) return pd.DataFrame(rows) features_df build_feature_frame(signal, fs25600) print(features_df.shape) # (约900000, 31)31 时间列 30 个特征时间开销主要在stats.skew和stats.kurtosis的纯 Python 循环上10 小时、9000 万点的信号会切出约 90 万个窗口单线程跑完大概几十分钟。常见做法是用multiprocessing或numba把两个特征函数并行化或者用np.lib.stride_tricks.sliding_window_view生成窗口矩阵再向量化。内存上 90 万行 × 31 列 float64 约 220 MB可以接受内存紧张时就把特征行按块写入 parquet不必一次全部驻留。4.3 退化曲线长什么样三条典型模式特征序列生成后先画三张图再谈建模。RMS 曲线通常呈现“平稳段—缓升段—剧升段”的三段式对应正常磨损、裂纹扩展和严重剥落峭度曲线的形态是早期冲高、中期回落的倒 V 型这是冲击从稀疏到连续的过渡重心频率 FC 和 RMSF 正常段基本平直退化后期向高频漂移。如果出现异常——比如 RMS 初期就持续上涨、峭度全程无波动、FC 开局就往高频跑——优先检查传感器量程饱和、安装松动和工频干扰不要急着改特征定义。曲线形态正常特征才算真正进入可用状态。5. 特征筛选与早期报警验证两个指标加一条切分纪律5.1 用单调性和趋势相关过滤 30 个特征退化特征要进趋势模型基本要求是单调。两个指标能快速完成初筛import numpy as np def monotonicity(s): 单调性0~1越接近 1 越单调。 d np.diff(np.asarray(s, dtypenp.float64)) return abs((d 0).sum() - (d 0).sum()) / len(d) def trend_corr(s): 与时间轴的线性相关给出方向和强度。 t np.arange(len(s)) return np.corrcoef(t, s)[0, 1]对 30 列特征逐列计算单调性低于 0.3 的偏度、波形因子这类来回摆的特征直接淘汰trend_corr给出正负号负相关表示特征在退化中下降同样是有效信息。按经验RMS、能量、峰峰值、谱能量的单调性通常排在最前峭度的单调性反而不高因为它中期回落。这两个指标跨数据集时排序基本一致可以直接作为筛选依据。5.2 早期报警点对比与时序切分纪律把全寿命按 5% 分段对每一段求峭度和 RMS 的均值取前 510 段做正常基线首次越限基线均值 3σ的时刻就是该特征的报警点。峭度通常比 RMS 提早 10%20% 寿命报警这正是早期点蚀的窗口期也是把 30 个特征全部算出来的意义所在。建模划分训练集和验证集时全寿命数据严禁随机打乱要用 TimeSeriesSplit 保持时间顺序标准化参数只从训练集拟合验证集按同一套均值和方差变换。这一步漏掉前 30 个特征里量级大的几个会把模型带偏比任何调参问题都致命。本文还有配套的精品资源点击获取