简介面向信号处理、状态监测与故障诊断方向的学习者这份压缩包内仅含一个脚本文件体积约2KB代码轻量、易于阅读和修改。脚本围绕振动信号的特征提取展开既包括均值、标准差、峰值、峭度等时域指标也利用快速傅里叶变换计算频谱功率、频率峰值与带宽等频域参数并集成了前六阶CEEMDAN分解用于处理非平稳、非线性信号从而捕捉不同时间尺度上的故障特征。对于机械故障诊断、振动台试验数据分析等场景该脚本提供了一套完整可运行的参考流程从数据预处理到特征计算一步到位适合正在开展相关课题的学生或需要快速验证算法的工程师借鉴复用。目前已有416人学习浏览凭借简洁的代码结构和清晰的注释可直接套用或扩展为自定义的信号分析工具。1. 振动信号特征提取要解决的不是“算特征”而是“特征是否真的可用”在设备状态监测和故障诊断的现场振动信号特征提取常常被误读成“对波形做一次统计计算”。从 data_processing.rar 这类打包程序里跑出十几个特征值并不难难的是搞明白哪几个指标能区分正常与早期故障哪几个指标会被转速波动和负载变化直接淹没。时域统计量描述幅值分布与能量总量频域特征描述能量在频率轴上的位置两者互补但不冗余必须放到同一个程序流程里统一核算。接下来的内容按时域特征、频域特征、特征提取程序、特征校验四个环节推进。每个特征都会写出数学定义、代码实现、参数边界和常见误用第 4 章提供一个可以直接替换数据源的 Python 骨架第 5 章给出投入模型前的三个校验开关。目标不是让你多攒出几百个特征而是让写进监测系统里的每一个特征都知道自己负责什么失效时又会以什么形式暴露出来。2. 振动信号时域特征先做数据清洗再算统计指标2.1 data_processing.rar 这类程序为什么强制“先去均值”拿到 data_processing.rar 这类打包的振动信号特征提取程序时域特征模块的第一行代码十有八九是x x - np.mean(x)。它不是为了凑篇幅而是采集链路本身就往信号里塞了直流偏置压电加速度计内置电荷放大器的零点漂移、采集卡输入通道的偏置电压、A/D 转换器的失调误差都会在原始波形上叠加一个常数。一个 0.01 量级的直流偏置叠加在 0.05 的振动幅值上均方根值被整体抬高峭度被压扁峰值因子失真后端的阈值判断会全面跑偏。去均值之后还要看趋势项。趋势项是信号在较长窗口内缓慢滑动的低频分量多由温度漂移和采样电路响应造成它让信号均值不再稳定在零附近也让频谱 0Hz 附近的能量异常抬高。常见处理有两条路一阶差分去趋势实现最简单但会把低频段的高频噪声一起放大最小二乘拟合基线再相减适合趋势近似线性的长记录但遇到突发瞬态时容易过拟合。我的习惯是先用短时傅里叶变换看一眼低频能量有没有随时间移动只有 0.5Hz 以下的频带持续抬升才做去趋势如果只是固定直流偏置去均值已经足够。预处理做完时域特征才有意义。均方根值描述振动能量总量峰值描述瞬时冲击的最坏情况峰值因子是峰值与 RMS 的比值用来突出冲击成分峭度是四阶中心矩与方差平方的比值专门刻画脉冲型故障。波形因子、脉冲因子和裕度因子在单测点诊断里出场频率低一些但做机器学习多维特征时仍然值得保留它们对幅值尺度变化的敏感度不同能提供一定的互补信息。2.2 时域特征数学定义与适用条件特征数学定义诊断倾向使用注意均方根值 RMSsqrt(mean(x^2))振动能量总量对早期局部缺陷不敏感峰值 Peakmax(abs(x))瞬时冲击幅值易被毛刺污染峰值因子 Crest FactorPeak / RMS冲击性故障扩展后可能回落峭度 KurtosisE[(x-μ)^4] / σ^4脉冲型早期故障样本太短时离散度极大波形因子RMS / mean(abs(x))波形形状波形近似正弦时变化不大裕度因子Peak / (mean(abs(x)))^2早期磨损信号接近零时分母爆炸峭度值得单独说。正常振动接近高斯分布峭度在 3 附近早期剥落产生的冲击稀疏但幅值大峭度可以升到 4 到 10 以上。故障进入连续磨损阶段后冲击变密信号逐渐趋近窄带随机过程峭度反而回落。这就是典型的“早期敏感、晚期钝化”。峭度单独用会漏掉晚期故障所以它要和 RMS 配合RMS 持续上升而峭度同时回落往往意味着故障已经从点蚀扩散成面磨损维修策略要从监测转为计划停机。2.3 用 Python 计算振动信号时域特征的最小实现import numpy as np def time_domain_features(signal: np.ndarray, fs: float 10000.0) - dict: signal: 一段去均值后的原始振动波形 fs: 采样率用于记录不参与统计计算 x signal - np.mean(signal) # 去直流偏置 n len(x) rms np.sqrt(np.mean(x ** 2)) peak np.max(np.abs(x)) mean_abs np.mean(np.abs(x)) std np.std(x) kurtosis np.mean((x - np.mean(x)) ** 4) / (std ** 4 1e-8) features { rms: rms, peak: peak, crest_factor: peak / rms if rms 0 else 0.0, kurtosis: kurtosis, waveform_factor: rms / mean_abs if mean_abs 0 else 0.0, peak2peak: np.ptp(x), # 峰峰值用于描述轴承冲击范围 } return features逻辑并不复杂先减均值消除直流偏置再算 RMS、峰值、峭度和三个比值。峭度除以std**4是为了归一化幅值尺度不同转速、不同负载下的信号可以互相比较crest_factor和waveform_factor都做了分母保护防止设备停机时出现除零。实际工程版本里我会在返回字典里再加一个sample_length字段因为后续做滑窗批处理时窗长不同算出来的峭度离散度差异很大记录样本数能帮你回溯特征波动是不是窗口长度引起的。2.4 窗口长度与采样率对时域特征的影响时域特征看似简单窗口长度却是第一个坑。峭度是四阶统计量样本越少波动越大经验值是样本点数低于 2000 时峭度就开始变得不可靠低于 500 时基本没有诊断价值。按 10kHz 采样率估算2000 点对应 0.2 秒而 1500r/min 的设备每转 0.04 秒0.2 秒只有 5 个转频周期样本依然偏少。我通常会取 1 秒窗口也就是 25 个转动周期这样峭度的逐窗波动能压到可接受范围。RMS 的稳定性要求低得多50ms 就能得到工程上可用的估计但峰值因子对窗口长度更敏感窗口太短间隔较长的冲击会被漏掉。所以窗口长度的下限应该服从峭度和峰值因子而不是 RMS。写成经验公式就是窗口长度 T K / (RPM/60)K 取 10 到 30。K 取 10 适合快速巡检K 取 30 适合长期监测的基线统计。这个选择会直接传导到第 4 章的滑窗参数里先在这里定下来。提示设备振动监测里峭度阈值常取 3.5超过 3.5 的窗口先标成“疑似冲击”再结合 RMS 和频域边带比决定是否报警不要单点触发。3. 振动信号频域特征从 FFT 到功率谱、边频带与谱峭度3.1 时域特征给不了“能量分布”频域特征补上时域特征能回答“振动有多大、冲击有多强”回答不了“这些能量集中在哪些频率”。两个轴承故障的 RMS 完全一样一个是外圈正常磨损一个是外圈出现剥落坑能量总和相同但频率分布完全不同。齿轮断齿会在啮合频率两侧激起边频带轴承早期剥落会激起某个结构共振频段这些差异在时域波形里肉眼很难分辨频域谱上一眼就能看出来。频域提取的第一步是 FFT但直接用np.abs(np.fft.fft(x))得到的不是工程意义上的幅值。FFT 结果是复数谱实信号频谱关于奈奎斯特频率对称通常只保留正频率部分同时 FFT 的幅值受窗长和窗函数影响不做修正就无法跨窗比较。很多已经跑起来的振动信号特征提取程序问题恰恰出在这里换了窗口长度同一测点的同一特征频带幅值出现 30% 以上的偏差却找不到原因。3.2 幅值谱、功率谱、功率谱密度怎么选谱类型计算方式物理含义适用场景幅值谱abs(FFT(x))单位频率分量的幅值看特征频率幅值的绝对值功率谱abs(FFT(x))^2频率上的能量识别主频带、边频带功率谱密度功率谱 / 频率分辨率单位频带功率跨采样率、跨窗长比较幅值谱直观幅值单位就是振动单位功率谱对小幅值分量更敏感因为平方运算把小信号压到更接近噪声底动态范围更大找边频带一般用功率谱。功率谱密度多一个除以频率分辨率的操作当采样率和窗长不一致时不同记录之间的谱值可以互相比较。如果设备状态监测系统长期固定采样率和窗长功率谱就够用如果要横向对比不同设备、不同工况最好统一存功率谱密度。还有一个取舍特征提取程序是保存整条频谱还是只抽取频带特征。整条频谱信息全但占用空间大直接丢给机器学习模型容易过拟合到噪声线上。我常用的策略是保存三个量主轴频及其整数倍频的幅值、边频带宽带的能量、若干预设频段的能量占比。这样既压缩了存储又保留了诊断和训练模型所需的主要信息。3.3 边频带提取齿轮轴承故障的频域特征怎么定齿轮箱和滚动轴承的故障信号一般都带调制成分。齿轮啮合频率 fm z × fr 两侧会出现间隔为转频 fr 的边频带轴承外圈故障特征频率 BPFO 附近同样会出现转频间隔的边带族。频域特征提取的关键不是看主频幅值而是看边频带能量与中心频率幅值的比值。这个边带比能抵消转速小幅波动和载荷变化的影响比单独的幅值稳定得多。边频带提取的实现里窗函数的选择直接决定边带比值可不可信。矩形窗的旁瓣只有 -13dB主频能量会泄漏到两侧边带边带比被系统性抬高汉宁窗旁瓣约 -31.5dB泄漏压小了一个量级。所以提取边频带特征时必须用汉宁窗或旁瓣更低的布莱克曼-哈里斯窗。矩形窗只有在一个窗口内正好包含整数个信号周期时才安全而振动信号几乎不可能满足这个条件。3.4 用 Python 提取频域特征加窗、修正幅值、算边带比import numpy as np def frequency_features(signal: np.ndarray, fs: float, rpm: float 1500.0) - dict: signal: 去均值后的时域信号 fs: 采样率单位 Hz rpm: 转速用于把转频换算成谱线条数 n len(signal) win np.hanning(n) # 汉宁窗压低旁瓣 signal_win signal * win spec np.fft.rfft(signal_win) # 实信号只算正频段 freq np.fft.rfftfreq(n, d1.0 / fs) # 每条谱线对应的频率 # 单边幅值修正除以 n、乘 2再补偿窗函数能量增益 amp 2.0 * np.abs(spec) / n amp amp / np.mean(win) power amp ** 2 # 单边功率谱 start_idx max(5, int(10.0 * n / fs)) # 忽略 10Hz 以下噪声平台 peak_idx start_idx np.argmax(power[start_idx:]) peak_freq freq[peak_idx] # 边带带宽取半转频带宽转频换算成谱线数 rotation_freq rpm / 60.0 sideband_half int(0.5 * rotation_freq / (fs / n)) lo max(start_idx, peak_idx - sideband_half) hi min(n // 2, peak_idx sideband_half) sideband_energy (np.sum(power[lo:peak_idx - 1]) np.sum(power[peak_idx 2:hi 1])) features { peak_freq: peak_freq, peak_amp: amp[peak_idx], total_energy: np.sum(power), sideband_ratio: sideband_energy / (power[peak_idx] 1e-8), } return features代码逻辑分四块加窗、FFT、主峰搜索、边带能量计算。np.hanning(n)加窗是为了压低频谱旁瓣amp 2.0 * np.abs(spec) / n把 FFT 的归一化修正成单边幅值再除以np.mean(win)补偿加窗造成的能量损失两步都缺一不可。rfftfreq返回正频段谱线频率最后一个索引是n // 2对应奈奎斯特频率。sideband_energy分别从两个方向跳过紧挨主峰的谱线避免汉宁窗主瓣泄漏灌进边带统计。sideband_half是本函数最需要注意的参数。它把半转频带宽换算成谱线条数换算公式是带宽除以频率分辨率fs / n。采样率 10kHz、窗长 1 秒时分辨率是 1Hz25Hz 的转频对应 25 条谱线边带搜索范围足够窗长压到 0.2 秒时分辨率变成 5Hz转频只对应 5 条谱线边带比仍然可用但波动加大。实际程序里 rpm 不应该当成固定标量最好从转速通道同步读入否则转频波动会直接扭曲边带比。3.5 谱峭度定位共振频段的时域峭度升级版时域峭度只告诉你“信号有冲击”不告诉你冲击能量在哪个频段。谱峭度把峭度按频率展开能定位到局部缺陷激励起的共振频带。常见实现是短时傅里叶变换对信号做 STFT 得到时频矩阵再沿时间轴对每个频率点计算四阶中心矩与二阶矩平方的比值整张时频图转化成一维谱峭度谱。谱峭度值最大的频带通常就是缺陷激励起的高频共振带。下一步对这个频带做带通滤波再做希尔伯特变换求包络包络谱里的特征频率会清晰得多。这是轴承故障诊断的标准链路也是许多振动信号特征提取程序里“包络分析”按钮背后的实现逻辑。它比直接对原始信号做包络谱更稳健因为它先排除了其他频段的干扰。4. 振动信号特征提取程序落地滑窗、批处理与特征表设计4.1 定长窗口还是转速跟踪窗口特征提取程序里第一组要决定的是窗口长度和滑动步长。定长窗口实现最省事但对变转速设备不友好转速从 600 到 3000r/min 变化时固定 1 秒窗口对应的转数从 10 圈到 50 圈峭度和边频带的统计对象完全不同特征值不再可比。转速跟踪窗口则根据当前转频动态计算窗口长度保证每个窗口包含相同的转数适合风电、机床主轴这类转速多变的场景。没有键相通道时可以用瞬时频率估计代替。希尔伯特变换求瞬时相位再对相位求导得到瞬时频率反推窗口需要的采样点数。这个办法不需要额外硬件但低转速区间瞬时频率估计误差偏大窗口长度会跟着抖。我通常按设备的转速波动率做判断没有键相、转速波动小于 ±5%用定长窗口同时把每个窗口的 RMS 存下来建模时按 RMS 分桶处理。4.2 构建“时域 频域”联合特征表的 Python 骨架import numpy as np import pandas as pd def build_feature_table(waveforms: np.ndarray, fs: float, rpm: float, win_sec: float 1.0, step_sec: float 0.5) - pd.DataFrame: waveforms: (测点, 样本) 二维数组 fs: 采样率 rpm: 名义转速用于边频带计算 win_sec: 窗口长度 step_sec: 滑动步长 win_len int(fs * win_sec) step_len int(fs * step_sec) rows [] for ch_idx, channel in enumerate(waveforms): for start in range(0, len(channel) - win_len 1, step_len): seg channel[start:start win_len] tfeat time_domain_features(seg, fs) ffeat frequency_features(seg, fs, rpmrpm) row {channel: ch_idx, start: start} row.update(tfeat) row.update(ffeat) rows.append(row) return pd.DataFrame(rows)这个骨架把第 2 章的时域特征函数和第 3 章的频域特征函数组装成一张特征表。win_sec1.0和step_sec0.5表示相邻窗口重叠 50%特征行数接近不重叠时的两倍适合训练机器学习模型但重叠窗口的样本不独立模型评估时要做分组交叉验证否则测试集里混进强相关的训练样本准确率高得没有意义。win_sec step_sec时窗口不重叠样本独立性最好适合长期趋势监测。注意frequency_features里的rpm在骨架里传的是标量名义转速。实际工况下转频会随负载波动更合理的做法是从同一时段的转速通道逐一读取瞬时转速再按窗口中值转速计算边带带宽。如果设备没有转速通道退而求其次用名义转速但要记录rpm_valid标记转速偏离超过 10% 的窗口打上不可用标签。注意滑动重叠窗口产生的特征行不是相互独立的。做模型验证时先按“通道 起始时间”分组再执行分组交叉验证能避免把同一个大窗口的孪生样本同时分进训练集和测试集。4.3 多文件批量处理与并行度设置data_processing.rar 这类程序被打开得最多的场合是几十个测点、几百个历史数据文件的状态评估。单线程计算其实不会太慢1 秒窗口做一次时域特征加 FFT耗时在毫秒量级真正的瓶颈在文件读取、解压和特征落盘。所以并行优化的重点不是疯狂调用 FFT而是合理调度 I/O。以 Python 为例concurrent.futures.ProcessPoolExecutor按文件分块即可每个子进程只接收文件路径和参数不传波形数据本身。进程数不是越多越好机械硬盘下多进程并发读反而会因磁头争用变慢SSD 上进程数可以接近物理核数。如果数据已经加载进内存用joblib或multiprocessing按测点分块也能拿到接近线性的加速比。4.4 特征存储格式与缺失值策略特征表推荐用 Parquet 格式按channel date分区存储。Parquet 是列式存储按特征列做统计分析时只扫描需要的列块比 CSV 快 5 到 10 倍CSV 的唯一优势是肉眼可直接打开但特征表动辄几十万行这个优势基本用不上。落盘时给文件加多个时间戳字段而不是把时间编码在文件名里方便后续按时间段切片。缺失值主要来自三类情况设备停机时振动幅值接近零比值类特征分母为零转速信号丢失导致无法计算转数传感器断线导致整段数据无效。停机段落直接用 RMS 阈值过滤不进入特征表转速标记为 0 的窗口保留在表里但加rpm_valid0标签比值类特征计算时返回np.nan不要返回 0。用 0 填充缺失值会把“无数据”和“真实为 0”混淆模型学到的边界是错的。5. 验证振动信号特征可用的三个开关类间距、漂移与迁移5.1 用 Fisher 类间距过滤无效特征特征表生成后第一件事不是训练模型而是看正常样本和故障样本的特征分布有没有重叠。Fisher 判别比等于 (μ1 - μ2)^2 / (σ1^2 σ2^2)大于 2 说明两类分布分得开小于 1 说明重叠严重即使分类器强行拟合换一段数据也会失效。可视化时画分布直方图而不是只画散点图直方图的尾部交叉情况比二维投影更真实也更容易看出哪些特征存在多模态分布。5.2 用 Spearman 单调性看退化过程把特征按设备投运时间排序计算 Spearman 秩相关系数数值接近 ±1 说明特征随退化过程单调上升或下降这是趋势预警类应用的基础。转速突变、负载跃迁会造出虚假趋势所以校验前先按工况分桶每个桶内单独算单调性。样本量少时改用 Mann-Kendall 趋势检验它对异常值不敏感伪趋势不容易通过检验。5.3 用工况切换做迁移测试换一段不同负载、不同转速的故障数据重跑一次同一组特征看峭度、边带比是否还保持同样的变化方向。失效时先查 RMS 是否淹没了冲击成分再查边频带带宽是否被转频变化拉偏。如果换一个工况就失去区分度这个特征在模型里的权重无论多高都不该被信任。提示迁移测试最怕“标签污染”。训练调参时已经看过验证工况的数据分布再拿同一批工况验证就会虚高。保留一段从未参与建模的数据放到最后只做一次测试。5.4 一个最小校验脚本import numpy as np from scipy.stats import spearmanr def quick_feature_check(normal: np.ndarray, fault: np.ndarray, trend: np.ndarray) - dict: fdr ((np.mean(fault) - np.mean(normal)) ** 2 / (np.var(fault) np.var(normal) 1e-8)) rho, p spearmanr(trend, np.arange(len(trend))) return {fisher: fdr, spearman_rho: rho, p_value: p}normal和fault分别传正常、故障窗口的特征向量trend传按时间排序的退化过程特征值。fisher小于 1 的特征直接剔掉spearman_rho绝对值接近 1 且 p 小于 0.05说明特征有趋势保留价值。最后一步是把这个脚本放到不同工况的分桶数据上各跑一遍凡是只在单一工况下通过的特征都值得回看原始波形确认物理含义再决定要不要进模型。本文还有配套的精品资源点击获取