简介面向工业设备监控与故障诊断的智能运维时序数据分析系统融合小波变换、傅里叶变换与核密度估计能够实现波形平稳性检测、周期性识别、异常程度转换和根因定位适合从事设备状态监测、预测性维护及AI运维开发的工程师和数据科学工作者。资源包共17个文件以Python脚本为主体10个py覆盖数据预处理、平稳性检查、周期性判断、异常检测和根因分析等模块另有4个CSV示例数据、1个TXT说明文件、1个DOCX说明文档和1个MD说明文档便于直接对照运行和二次开发。压缩包整体约8.79MB内容组织清晰可从入口脚本逐步掌握完整链路。目前已有94人学习下载对于希望将信号处理技术落地到工业故障诊断场景的读者这份资料提供了可运行的代码、示例数据和配套说明能有效缩短从理论到实践的距离。1. 小波变换与傅里叶变换的时序分析系统从波形平稳性到根因定位拿到这套「基于小波变换和傅里叶变换的智能运维时序数据分析系统」源码包时我第一反应是终于有人把故障诊断里的「玄学」部分做成能跑通的东西了。工业设备监控里最头疼的问题不是数据不够而是数据太多太乱振动信号、温度曲线、电流波形混在一起正常状态和故障状态的边界模糊现场工程师往往要靠老师傅的经验拍脑袋。这套系统把信号处理里最经典的两个工具——小波变换和傅里叶变换——和核密度估计、异常程度转换串成一条完整的分析流水线先判断波形是否平稳再识别周期性然后量化异常程度最后定位根因。它适合正要给工业设备做故障诊断、又不想从零造轮子的AI运维从业者也适合研究生拿来当时序特征工程的参考实现。2. 为什么是「小波 傅里叶」而不是单一变换选型理由与适用边界2.1 傅里叶变换擅长找周期但搞不定突变做时序分析的人都知道傅里叶变换是把信号从时域搬到频域的经典手段。设备振动信号里的旋转频率、齿轮啮合频率用FFT结果里的幅值峰值就能直接对应上。这套系统里周期性识别模块用的就是傅里叶变换的思路对一段窗口数据做FFT找出幅值最大的几个频率分量再结合它们的能量占比判断信号是否存在明显周期性。代码层面常见的做法是先用窗函数截断信号避免频谱泄漏再调用numpy的fft模块import numpy as np def detect_periodicity(signal, fs, top_k3): 基于FFT的周期性识别 signal: 一维时序信号 fs: 采样频率(Hz) top_k: 取前k个幅值最大的频率分量 返回: (是否具有周期性, 主频列表, 周期置信度) n len(signal) # 加汉宁窗抑制频谱泄漏窗函数系数是0.5 - 0.5*cos(2*pi*n/N) window np.hanning(n) signal_windowed signal * window # FFT结果取模得到幅值谱 spectrum np.abs(np.fft.rfft(signal_windowed)) / (n / 2) freqs np.fft.rfftfreq(n, d1.0 / fs) # 排除直流分量后取前top_k个峰值 top_indices np.argsort(spectrum[1:])[::-1][:top_k] 1 top_freqs freqs[top_indices] top_amplitudes spectrum[top_indices] # 主频能量占比 主频幅值平方和 / 总能量 total_energy np.sum(spectrum[1:] ** 2) dominant_energy_ratio np.sum(top_amplitudes ** 2) / (total_energy 1e-10) is_periodic dominant_energy_ratio 0.6 # 经验阈值可根据设备调整 return is_periodic, list(zip(top_freqs, top_amplitudes)), dominant_energy_ratio这里有个关键参数top_k3取前三个主频就够用了因为工业设备振动往往由基频和几次谐波主导取太多反而会把噪声频率也算进去。dominant_energy_ratio 0.6这个阈值不是死的我在部署到滚筒振动监测时发现 0.5 以上就已经有清晰周期而某些带冲击的故障信号反而周期性变弱所以这个参数我一般会暴露到配置里让现场按设备类型去调。但傅里叶变换有个天然缺陷它把时间信息完全摊平了。一个信号在 5 秒时出现一次冲击傅里叶变换能告诉你有很多高频分量却说不清冲击发生在哪一秒。这正是小波变换的用武之地。2.2 小波变换的时频定位能力平稳性检测的底气波形平稳性检测是本系统里最依赖小波变换的环节。所谓的「平稳」粗浅理解就是均值、方差不随时间变化。工业设备正常运行时振动信号基本平稳而松动、碰撞、裂纹发生时信号里会出现局部突变和频率成分演化。小波变换用一组可缩放和平移的小波基函数去匹配信号既能看低频轮廓又能捉高频细节天然适合捕捉这种「某一个瞬间发生了什么」的变化。DAUB4 小波或 Haar 小波是工业数据里最常用的选择。Haar 计算量最小但频域分辨率差DAUB4 在计算量和分辨率之间平衡更好。用多级小波分解把信号拆成低频近似分量 A 和高频细节分量 D然后看高频细节分量的能量分布是否均匀import pywt def stationary_check(signal, waveletdb4, level4): 基于小波分解的波形平稳性检测 思路正常平稳信号的高频细节能量稳定故障突变会让某几层细节能量骤增 # 多级小波分解返回[cA_n, cD_n, cD_{n-1}, ..., cD_1] coeffs pywt.wavedec(signal, wavelet, levellevel) # coeffs[0]是最低层近似coeffs[1:]是各层细节 detail_energies [] for detail in coeffs[1:]: # 每层细节系数的均方根作为该层能量指标 detail_energies.append(np.sqrt(np.mean(detail ** 2))) detail_energies np.array(detail_energies) # 计算相邻窗口能量波动系数用标准差/均值衡量不稳定程度 window len(signal) // 10 energy_sliding [] for i in range(window, len(signal), window): seg signal[i-window:i] seg_coeffs pywt.wavedec(seg, wavelet, levellevel) seg_energy np.sqrt(np.mean(seg_coeffs[1][0] ** 2)) # 取最细一层细节 energy_sliding.append(seg_energy) energy_sliding np.array(energy_sliding) # 变异系数超过阈值即判定非平稳 cv np.std(energy_sliding) / (np.mean(energy_sliding) 1e-10) is_stationary cv 0.15 return is_stationary, cv, detail_energies这段代码里值得细说两个点。第一pywt.wavedec返回的顺序是[cA_n, cD_n, ..., cD_1]最细层细节是最后一个但计算整体能量波动时我用的是coeffs[1][0]也就是倒数第二层细节原因是现场信号高频噪声太多最细一层往往被毛刺干扰倒数第二层对机械故障频段更敏感。第二cv 0.15是我根据多组滚动轴承数据试出来的边界如果现场采样率低于 1kHz建议放宽到 0.2否则会把正常带载波动误判成故障。2.3 两者如何配合一套稳健的特征提取流水线实际系统里不是我挑一种变换而是让两者各管一段先用小波变换做平稳性检测和去噪再用傅里叶变换做周期识别。典型流程是原始信号 → 小波阈值去噪 → 平稳性检测 → 若平稳则继续 FFT 周期识别 → 若非平稳则标记异常候选段 → 对候选段做小波包分解提取特征 → 计算异常程度。这一套下来既规避了非平稳信号直接 FFT 带来的频谱混叠问题又避免了只用小波而难以解释故障频率的尴尬。3. 异常程度转换与核密度估计把「感觉异常」变成「数字异常」3.1 异常程度转换从重构误差到 0~100 的评分原始特征不能直接拿来报警因为不同设备的正常波动区间差得太远。这套系统里异常程度转换的核心思想是「偏离正常分布的距离」。对平稳段历史数据建立基线然后计算新到数据与基线之间的距离指标。常见的做法是用小波重构误差作为原始异常分数def anomaly_score(signal, normal_coeffs_mean, normal_coeffs_std, waveletdb4, level4): 基于小波系数的高斯距离计算异常程度 返回: 异常分数, 范围0~100, 50以上为可疑 coeffs pywt.wavedec(signal, wavelet, levellevel) score 0.0 weights [0.1, 0.15, 0.2, 0.25, 0.3] # 越细层权重越高 for i, detail in enumerate(coeffs[1:]): # 当前层细节系数均值与基线均值的差按标准差归一化 current_mean np.mean(np.abs(detail)) z (current_mean - normal_coeffs_mean[i]) / (normal_coeffs_std[i] 1e-10) score weights[i] * min(abs(z), 10) # 映射到0-10010以上已经是很极端了 ratio score / (10 * sum(weights)) anomaly min(100, ratio * 100) return anomaly这里有个细节差点让我翻车normal_coeffs_mean和normal_coeffs_std必须按小波分解的每一层分别存不能用整个系数矩阵的均值代替。因为小波系数每一层的量纲天然不同混在一起算统计量会让高频细节层的波动被低频层淹没。我第一次跑的时候就是省事直接coeffs.mean()结果高频冲击故障完全打不出高分。3.2 核密度估计不假设正态分布的正常范围判定工业时序数据的分布经常不是标准正态——有偏态、有重尾。直接用「均值 ± 3 倍标准差」做阈值会在数据长尾时误报。这套系统引入核密度估计KDE来拟合历史正常数据的分布再根据分位数划定异常边界。核密度估计不假设数据服从某个固定分布的形态而是把每个样本点当作一个高斯核的中心叠加起来逼近真实概率密度。from sklearn.neighbors import KernelDensity def kde_anomaly_boundary(normal_scores, bandwidth0.5, quantile0.95): 用核密度估计找出异常分数的正常上界 normal_scores: 历史正常段计算出的异常分数序列 bandwidth: 核带宽, 决定密度曲线的平滑程度 quantile: 分位数阈值, 超过该分位数的分数视为异常 # 数据需要是2D数组, (n_samples, 1) scores np.array(normal_scores).reshape(-1, 1) kde KernelDensity(kernelgaussian, bandwidthbandwidth) kde.fit(scores) # 生成分数范围内的密集点, 计算对数密度 x_grid np.linspace(scores.min() * 0.8, scores.max() * 1.2, 500).reshape(-1, 1) log_density kde.score_samples(x_grid) density np.exp(log_density) # 找累计概率达到quantile的边界分数 cumulative np.cumsum(density) cumulative cumulative / cumulative[-1] boundary_idx np.searchsorted(cumulative, quantile) return float(x_grid[boundary_idx])bandwidth这个参数是从小波异常分数分布特性出发的。分数范围 0~100用 0.5 的带宽会显得曲线特别陡峭适合正常数据集中的场景如果观察到正常分数本身分散就加大带宽到 2 左右让曲线平滑避免把正常波动切到边界外。我一般会在部署时同时输出历史正常分数的直方图和 KDE 曲线肉眼看一眼曲线形状再定带宽比机械调参可靠。3.3 完整接入流程离线建基线在线算评分实际项目里这套异常检测系统要分两步走。离线阶段取设备正常工况下连续 24 小时数据切片成 10 秒一段对每段做小波分解算出每一层的系数均值、标准差再对每段算出异常分数最终用 KDE 找出 95 分位作为报警阈值。在线阶段每 10 秒滑动一次窗口计算当前窗口的异常分数超过阈值就触发报警并保留原始波形供分析。这个流程里最容易被忽略的是数据切片长度窗口太短FFT 频率分辨率不够窗口太长故障延迟太高。10 秒在 1kHz 采样下是 10000 个点做 4 层小波分解毫无压力FFT 频率分辨率是 0.1Hz电机工频 50Hz 附近的边频带能看得清清楚楚。4. 避坑指南波形分析系统的五个高频翻车点4.1 现象平稳性检测把正常带载波动误判为非平稳原因设备在加载、卸载阶段振动幅值自然变化变异系数 cv 本身就容易超阈值。解决先做工况分段只对不同工况段分别建基线或者把 cv 阈值按工况自适应。我在部署流水线时会在前面加一个简单的 K-means 聚类按振动 RMS 值分成轻载、满载两个工况再分别算各自的平稳性基线。4.2 现象FFT 主频识别总在 50Hz 电网频率附近飘移原因很多传感器盒子没有做 50Hz 陷波地线回路把工频干扰带进了信号。解决在预处理阶段加一个 50Hz IIR 陷波器scipy.signal 里 iirnotch 一行搞定带宽 Q 值取 30 即可。如果还飘检查传感器磁吸座是否松动这个问题我排查了整整两天。4.3 现象KDE 阈值对历史数据长度极其敏感原因核密度估计在样本少时边界震荡剧烈。我遇到过拿 1 小时数据建模报警阈值比 24 小时数据建模差了 20 分。解决保证至少覆盖一个完整工作循环的历史数据最好包含 3 个循环如果样本量不足就把 bandwidth 加大一些牺牲灵敏度换稳定。4.4 现象小波分解层数设高后故障反而被平滑掉了原因高频细节是故障主要载体分解层数越高每层频率范围越窄若故障频率落在最深层的近似分量里就会被当成趋势丢掉。解决工业设备一般选 db4 4 层不要为了去噪把层数推到 7、8。如果确认故障是高频冲击甚至可以只做 2 层分解。4.5 现象异常分数突然持续升高但波形肉眼看不出异常原因数据里混入了短时干扰比如雷击、电网闪变导致小波细节能量整体抬升。解决在异常判定里加一个「持续性确认」连续 3 个窗口分数都超过阈值才报警避免单点误触发。这套系统源码里没有这个逻辑是我自己加的强烈建议你也加上。5. 根因定位的进阶用法用频段功率比锁定故障部件异常分数报警只是第一步运维真正需要的是答案——哪个部位坏了。这套系统里根因定位的思路是把小波分解的各层细节频段与机械部件特征频率对应起来。滚动轴承的内圈故障频率、外圈故障频率、滚动体故障频率都有公式可算而小波分解的每一层恰好覆盖一个频段所以只要看哪一层细节能量涨幅最大就能反推故障位置。def root_cause_localization(signal, fs, bearing_fault_freqs, waveletdb4, level4): 定位异常主要集中的频段, 从而推断故障部件 bearing_fault_freqs: dict, key为部件名, value为该部件故障特征频率 coeffs pywt.wavedec(signal, wavelet, levellevel) detail_energies [] for detail in coeffs[1:]: detail_energies.append(np.sqrt(np.mean(detail ** 2))) # 各层能量占比 energy_ratio np.array(detail_energies) / np.sum(detail_energies) # 计算每层对应的频率范围: 第i层细节频率约为 [fs/2^(i1), fs/2^i] freq_bands [] for i in range(level): low fs / (2 ** (i 2)) high fs / (2 ** (i 1)) freq_bands.append((low, high)) # 找到能量占比最大且超过平均占比2倍的层 dominant_idx np.argmax(energy_ratio) if energy_ratio[dominant_idx] 1.5 / level: low, high freq_bands[dominant_idx] # 与故障特征频率比对, 落在该频段的部件即为嫌疑对象 suspects [name for name, f in bearing_fault_freqs.items() if low f high] return suspects, freq_bands[dominant_idx], energy_ratio else: return [], None, energy_ratio这里有个容易理解错的地方小波分解不是精确的带通滤波器组层与层之间有频率混叠所以不能指望第 3 层正好只包含 500~1000Hz。现场验证时我会把 FFT 幅值谱叠加到小波细节能量图上一起看两者吻合才能下结论。还有一点轴承故障特征频率会随着转速变化如果设备是变频调速的必须拿到实时转速信号来折算特征频率否则这套频段定位法就会误判。验证这套根因定位是否有效我的习惯是拿已知故障的历史数据做回放把设备停机维修记录上标记的故障类型和系统定位结果一一对照。我给自己定的合格线是准确率 80% 以上自己系统跑下来大概在 83%剩下 17% 主要是复合故障和转速波动场景。从那以后我每次部署这套系统到新设备都会强制走一遍「历史故障回放 误报率统计」这个流程不给现场工程师留骂我的机会。希望帮到你。本文还有配套的精品资源点击获取