简介一套面向工业设备监控与故障诊断的智能运维时序数据分析系统基于小波变换、傅里叶变换与核密度估计实现波形平稳性检测、周期性识别、异常程度转换与根因定位。适用于工业数据分析、设备健康管理和故障预测方向的开发者与研究者可帮助快速理解并搭建AI运维分析流程。资源共17个文件包含10个Python脚本、4个CSV数据集及说明文档整体约8.79MB脚本覆盖信号预处理、平稳性/周期性判断、异常检测与根因定位等核心环节CSV数据便于直接复现实验。目前已有94人学习下载。通过源码与配套说明既能学习小波多尺度分解与频域特征提取的具体实现也能掌握核密度估计在异常定位中的应用方式系统从数据读取、特征分析到异常评估给出了完整链路适合作为毕业设计、课题研究或实际工业监控项目的参考方案。1. 基于小波变换和傅里叶变换的智能运维时序数据分析系统解决的不只是“有没有故障”我在工业设备监控项目里见过太多这样的误判离心泵的振动频谱看起来完全正常转频分量也不突出但拆开轴承才发现已经严重磨损另一台电机电流波形到了某个时间段开始变“毛糙”FFT 峰值几乎没变化事后检查发现定子绝缘已经击穿。问题恰恰出在把傅里叶变换当成唯一的频谱工具平稳信号里那些瞬态成分被彻底摊平了。这套基于小波变换和傅里叶变换的智能运维时序数据分析系统核心做法是让两种变换分工傅里叶做全局周期识别和平稳性判断小波抓取局部突变与瞬时冲击再用核密度估计把重构残差映射成可告警的异常程度分数最后沿着频带能量的变化做根因定位。它面向的是设备状态监控工程师、故障诊断算法开发和 AI 运维平台建设者解决的是“看到异常但说不清异常从哪来”的落地问题。2. 傅里叶变换与小波变换的分工全局周期与局部突变的选型逻辑2.1 傅里叶变换在设备监控中能看见什么、看不见什么先明确一个基础事实傅里叶变换把一段时序信号整体变换到频域输出的是这段时长内各个频率分量的幅值和相位。对稳定运行的设备比如恒转速风机、齿轮箱、离心机频谱里基频和边频带是有规律的我通常会再叠加一层谐波乘积谱 HPS 来放大基频把“设备有没有周期”从人工看图变成自动判定。但代价也很直接一次 FFT 的频谱是整段窗口内的统计平均时间分辨率等于零。窗口里如果发生过一次短暂冲击它在频谱上的能量会被摊薄到整个观察时长里只在很宽的频段上留下一点不太起眼的凸起肉眼几乎无法识别。import numpy as np from scipy.fft import rfft, rfftfreq # 现场振动信号采样率 25600 Hz取 0.5 秒作为一次窗口 fs 25600 win_size int(fs * 0.5) # 12800 点 raw get_sensor_segment(device_idpump_01, start_ts..., lengthwin_size) # 去均值消除直流分量对低频段的压制 x raw - np.mean(raw) # 加汉宁窗限制截断处不连续产生的频谱泄漏 xw x * np.hanning(win_size) # 单边幅值谱rfft 只输出一半谱乘 2 恢复真实幅值 amp np.abs(rfft(xw)) * 2.0 / win_size freq rfftfreq(win_size, 1 / fs) # 切除 3 Hz 以下的极低频段只保留有物理意义的部分 keep freq 3 amp amp[keep] freq freq[keep]这段代码里有三个参数值得单独说明。第一个是加窗现场信号不可能正好截断在整数周期上端点不连续会让一个频率的能量泄漏到周围频段加汉宁窗是成本最低的缓解手段。第二个是频率分辨率它等于 fs 除以 win_size当前参数下分辨率是 2 Hz意味着两个间隔小于 2 Hz 的谱峰在频谱上只会显示成一个鼓包你要单独检测 0.5 Hz 的差异就必须把窗口拉长到 2 秒以上。第三个是单边谱的归一化rfft 输出的是复数频谱对实数信号只需要前半部分乘 2 再除以窗口长度才是物理上可读的幅值。这套参数是周期性识别和后续平稳性判断的共同起点换采样率时一定要同步改窗口长度不能只改其中一个。2.2 小波变换补上的局部时频能力与基函数选择小波变换不把信号拆成无限长的正弦波叠加而是用一组可伸缩、可平移的小波基去做内积输出不同尺度带上的系数。在设备监控场景下这相当于把信号按频带切开同时保留时间位置信息。细节系数对应高频瞬态近似系数对应低频趋势。我在这类系统里通常选 db4 或 db6 小波原因说得直白一点带点玄学成分但实测多年确实有效Daubechies 系小波的波形形态和机械冲击响应有相似性冲击能量能集中在少数几个系数上不像 Symlets 或 Coiflets 那样把一次冲击摊到好几层系数中给后续阈值判断带来额外麻烦。import pywt def decompose_signal(x, waveletdb4, level4): # wavedec 返回的元组第一位是近似系数后面依次是各层细节系数 coeffs pywt.wavedec(x, wavelet, modesymmetric, levellevel) names [fA{level}] [fD{i} for i in range(level, 0, -1)] return dict(zip(names, coeffs)) coeffs decompose_signal(x, waveletdb4, level4) # 采样率 fs25600Hz 时4 层分解各频带大致对应 # A4 约 0800HzD4 约 8001600Hz # D3 约 16003200HzD2 约 32006400HzD1 约 640012800Hz选层数时有几个边界要自己算清楚。每层系数长度减半D1 长度是原始信号的一半A4 长度约为窗口点数的十六分之一如果窗口本身只有 2048 点A4 只剩 128 点统计意义就很弱了。另一个是 Nyquist 频率小波每一层细节系数的频带是上一层的一半但原始采样率决定了最高只能识别到 fs/2也就是 12800 Hz。如果你关心的轴承故障特征频率落在 8001600 HzD4 就对应得很准如果关心的是 50100 Hz 的转频倍频就要把 level 加到 6 或 7或者直接回到傅里叶频谱去做。2.3 两种变换在一个系统里怎么串成管线我在这类监控平台上常把数据流拆成三个动作先用傅里叶做全局扫描用滑窗 FFT 算频谱一致性和周期主峰一旦平稳性判断触发「非平稳」标记就把后续特征计算切换到小波分解最后异常程度转换阶段再用小波重构残差作为输入。原因是两块计算成本差异很大FFT 一次窗口只做一组蝶形运算而小波分解要逐层卷积计算量大约是前者的 510 倍全部数据都跑小波在边缘网关上是浪费。只有被判定为存在瞬态或非平稳特征的数据段才值得用小波系数进一步分析。这样组织也方便加载历史故障样本做回放验证后面第六章会细说。下表是我在不同设备类型上的选型经验供你直接套用设备工况主要现象首选分析方式关注频带恒定转速轴承磨损周期冲击脉冲小波细节系数包络8003200Hz齿轮箱啮合异常边频带变化傅里叶高阶频谱啮合频率附近转子不平衡转频幅值稳定上升傅里叶基频趋势1倍转频变速工况下的喘振频率随时间漂移滑窗 FFT 频谱瀑布全频段观察漂移电机定子匝间短路电流波形毛糙小波细节系数能量突增高频噪声段3. 平稳性检测和周期性识别滑窗频谱相似度与谐波乘积谱3.1 平稳性不是统计书上的平稳而是频谱一致性的工程近似教科书里的平稳性定义是联合概率分布不随时间平移这个条件在工业现场几乎不可能满足也没有工程意义。我在智能运维平台里实际用的是另一套标准连续多个滑窗的幅度谱互相之间是否一致。设备平稳运转时每个窗口的频谱形状应该高度相似轴承出现局部剥落、叶片发生裂纹等瞬态源加入后某个窗口里会突然多出一段高频能量相邻窗口频谱相似度随即掉下去。这种近似虽然不像 ADF 检验那样有严格统计学证明但胜在能和频谱分析共用同一套滑窗机制计算开销极小。def spectral_similarity(amp_series): # 输入是多个窗口的幅值谱形状为 (num_windows, num_freq_bins) sims [] for i in range(1, len(amp_series)): a amp_series[i - 1] b amp_series[i] dot np.dot(a, b) norm_a np.linalg.norm(a) norm_b np.linalg.norm(b) # 加一个极小值避免除以 0 sims.append(dot / (norm_a * norm_b 1e-12)) return np.array(sims) # 假设已经用滑窗 FFT 拿到了 60 个窗口的幅值谱 # 窗口长度 0.5 秒重叠 50%覆盖 15 秒数据 window_size int(fs * 0.5) hop_size int(window_size * 0.5) amp_series [] # 每个元素是长度为窗口内的幅值谱数组 for start in range(0, len(x) - window_size, hop_size): seg x[start:start window_size] seg seg * np.hanning(window_size) amp_series.append(np.abs(rfft(seg))[: window_size // 2 1]) sims spectral_similarity(amp_series) # 平稳阈值SIM_TH 0.90 # 如果 sims 中连续 3 个值低于 0.90标记该时段为非平稳阈值 0.90 怎么定出来的多台设备离线数据回放后我发现正常定工况下频谱相似度通常在 0.93 以上偶发电源干扰会造成单点下降到 0.85但下一窗口立刻恢复真正进入故障发展期时相似度是持续走低并稳定在 0.85 以下。所以这里的判断逻辑不是看单点而是看连续低于阈值的窗口数量。如果你处理的信号是电流而非振动建议把阈值放宽到 0.88因为电流波形受负荷波动影响更大频谱变化比振动更频繁。3.2 周期性识别谐波乘积谱与倍频校验的组合手段周期性识别的目标是回答两个问题设备运行周期是多少这个周期是不是真实的机械周期。单靠 FFT 找峰值不够因为频谱里响应最大的是某次倍频而不是基频。我习惯在原始幅值谱基础上再加一层谐波乘积谱 HPS把谱分别下采样 1/2、1/3、1/4然后逐点相乘真实基频会在所有压缩谱中同时保留并叠加放大谐波分量则被相互错位抵消。def hps_spectrum(amp, num_products4): # amp: 原始幅值谱长度为 N max_idx len(amp) // num_products hps np.ones(max_idx) for k in range(1, num_products 1): # 每级下采样 k 倍截断到统一长度 compressed amp[::k][:max_idx] hps * compressed return hps # 峰值定位取 hps 中前 5 个局部极大值 peak_indices np.argsort(hps)[-5:][::-1] candidate_freqs freq[peak_indices] # 进一步用设备额定转频校验取最接近真实转速倍频的候选HPS 层数选 4 还是 5取决于频谱长度。num_products 越大基频权重越突出但有效频率范围会被压缩到原频谱的 1/5高频段信息丢失。我一般先限制检测范围在 3 Hz2000 Hz再确认 num_products 不大于 5。峰值定位后必须做一次倍频校验就是把候选频率和额定转频、叶片通过频率、轴承特征频率表对比能对上才写入周期特征库这种做法能把电源工频干扰误判成设备周期的问题挡在门外。3.3 平稳性结果怎么驱动小波分析接管平稳性检测结果不能只停留在曲线里要直接参与数据流选择。平台在实盘运行时按固定长度缓存原始信号每缓存一个窗口就做一次频谱相似度更新当相似度持续下降并触发非平稳标记时系统自动把这 1530 秒的原始缓存发送到小波分析模块而不是删除缓存继续等下一个窗口。这里有个常被忽略的工程细节小波分解需要足够长的序列来保证低层近似系数的统计意义如果触发发现得太晚缓存里只剩 2 秒信号分解到 A5 层几乎就剩不了几个点。我的做法是缓存区保持 60 秒、滑窗 0.5 秒标记触发时至少还有 30 秒数据可供小波分析。# 伪代码级的状态管理逻辑 STATE STABLE similarity_trend [] def on_window(amp_series): global STATE s spectral_similarity([amp_series[-2], amp_series[-1]])[0] similarity_trend.append(s) if len(similarity_trend) 10: similarity_trend.pop(0) # 连续 3 个窗口低于阈值才切换状态避免单点干扰 if all(v 0.90 for v in similarity_trend[-3:]): STATE UNSTABLE elif all(v 0.93 for v in similarity_trend[-3:]): STATE STABLE这里连续 3 个窗口的意义是消除偶然扰动。之前踩过一个坑只用单窗口相似度判定每次电压波动都会让平台误报“进入非平稳”小波分析模块被频繁唤醒网关 CPU 直接跑满。改成连续 3 次判定后误报率大幅下降这个参数也成了我以后做类似平台的默认值。4. 异常程度转换从重构残差到核密度估计的分数映射4.1 为什么要用重构残差定义异常程度直接拿振动幅值 RMS 做异常分级是最常见但最不可靠的做法。设备负载不同同一 RMS 值对应的健康状态完全不同同样幅值的信号在正常工况下可能只是负荷波动在重载工况下却可能是早期故障信号。我用的思路是先建立“正常模型”用小波重构逼近期望中的正常信号然后拿实际信号和期望信号的差当异常特征。这个差值叫重构残差幅值大、持续长才说明设备偏离正常状态比直接看绝对幅值靠谱得多。小波重构的做法是先对细节系数做阈值收缩把代表噪声的小系数置零或压缩再用处理后的系数反变换得到平滑信号。阈值的选择直接影响残差质量阈值太高重构信号过于干净残差里全是正常负荷波动后续异常程度会被夸大阈值太低阈值收缩等于没做残差和原始信号几乎相同。我一般用软阈值加自适应尺度参数尺度由分解层数和信号噪声水平共同决定。def wavelet_reconstruct_smooth(signal, waveletdb4, level4, scale1.2): coeffs pywt.wavedec(signal, wavelet, modesymmetric, levellevel) # 对每一层细节系数做软阈值收缩 thresholded [coeffs[0]] # 近似系数不动 for i, detail in enumerate(coeffs[1:]): sigma np.median(np.abs(detail)) / 0.6745 # 噪声方差估计 thresh scale * sigma * np.sqrt(2 * np.log(len(detail))) # 软阈值绝对值小于阈值的置零其余向零压缩 thresholded.append(pywt.threshold(detail, thresh, modesoft)) return pywt.waverec(thresholded, wavelet, modesymmetric) smooth wavelet_reconstruct_smooth(x, level4, scale1.2) residual x - smooth这段代码的核心在噪声方差估计 sigma 和阈值系数 scale。sigma 用中位数绝对偏差除以 0.6745是为了在受冲击信号影响时仍然稳健scale 是留给你的现场调节旋钮我通常在 1.0 到 1.5 之间调scale 越大重构信号越平滑检测灵敏度越高但误报风险也上升。注意 pywt.threshold 的 mode 用 soft 而不是 hard软阈值能避免重构信号在阈值点产生人工跳跃这个跳跃如果存在会直接形成伪冲击残差。4.2 核密度估计给残差建立分布模型拿到了残差序列下一步是回答“这个残差算不算异常”。最简单的方式是假设残差服从正态分布用均值加减 3 倍标准差做界限但残差往往带重尾和偏态特别是故障初期的瞬态冲击会让残差分布出现明显长尾。此时用核密度估计 KDE 更稳它不对分布形状做强假设通过把每个样本点变成一个小核函数并求和逼近真实概率密度。scipy.stats.gaussian_kde 提供现成实现核心要调的是带宽参数 bw_method。from scipy.stats import gaussian_kde # 用健康设备的历史正常残差训练 KDE # residual_train 是一维数组长度为 n kde_model gaussian_kde(residual_train, bw_methodsilverman) # 同时计算正常残差的 99% 分位数作为硬边界防止极端值干扰 abs_residual np.abs(residual_train) hard_limit np.percentile(abs_residual, 99) def residual_to_anomaly_score(residual_value, kde_model, hard_limit): # 求当前残差值在正常分布下的概率密度 density kde_model.evaluate([abs(residual_value)])[0] # 密度越低越位于分布尾部异常程度越高 # 用相对密度转成分数密度高于某个值记 0否则按低于比例分级 normal_density_99 kde_model.evaluate([hard_limit])[0] if abs(residual_value) hard_limit: return 0 if density normal_density_99 * 0.2: return 4 elif density normal_density_99 * 0.5: return 3 elif density normal_density_99: return 2 return 1这个分数映射里99% 分位数是工程上的安全垫。gaussian_kde 的带宽选 silverman 时是按样本量自动推导的适合大多数振动残差如果样本量小于 200建议手动传一个标量带宽值比如 0.05带宽太小会过拟合样本尖峰太大则会抹平双峰结构这是调整这个模型最需要耐心的地方。gaussian_kde.evaluate 输入要二维数组形式不写方括号直接传标量会报维度错误这是很多人第一次跑通代码时最容易卡住的地方。4.3 从密度到分数异常程度分级与告警联动异常程度转换最后要落到监控平台的告警联动。单独看一个窗口的分数意义有限因为瞬态冲击和电磁干扰也能造成单点高分但实际损害未必严重。我习惯把窗口分数送进一个 EWMA指数加权移动平均滤波器让分数平滑把瞬时抖动压下去只保留持续上升的趋势平滑系数 0.3 时一次孤立的冲击在经过 5 个窗口后影响几乎归零而持续上升的故障特征会稳定在一个高位。def ewma_smooth(scores, alpha0.3): smoothed [] # 平滑后的分数序列 prev scores[0] for s in scores: prev alpha * s (1 - alpha) * prev smoothed.append(prev) return np.array(smoothed) smoothed_scores ewma_smooth(raw_scores, alpha0.3) # 分级联动 # 01 级: 正常 2 级: 关注进入趋势跟踪 # 3 级: 预警通知现场检查 4 级: 停机建议平滑系数 alpha 是分级输出的关键参数。alpha 太大平滑几乎不起作用太小真实故障趋势要到很晚才浮现容易错过维护窗口。我一般把 alpha 控制在 0.250.35配合 10 秒一个输出周期使用如果你做的是慢变故障比如轴承润滑退化可以降到 0.1让趋势更稳定。分级阈值本身也不是固定不变的需要在平台上线后用历史故障样本回放校验这个落到第六章展开。5. 智能运维根因定位的排查与避坑频带映射、谐波误判与阈值漂移5.1 根因定位的事后链路异常时刻的频带能量反推异常程度分数触发预警后根因定位要回答“哪个部件出的问题”。我用的办法是频带能量反推原始信号经过小波分解后各层细节系数对应不同频带当异常分数升高时对比异常时段和正常基线的各频带能量占比找出能量变化最剧烈的频带再对照设备已知的特征频率表判断物理部件位置。轴承外圈故障的特征频率通常落在某个具体频邻齿轮啮合故障则集中在啮合频率及其边带电机转子断条表现为转频侧边带各有固定位置不需要靠猜。def locate_fault_subband(coeffs, baseline_energies): # coeffs: 小波分解系数baseline_energies 是健康时段各个子带的归一化能量 current_energies [] for c in coeffs: energy np.sum(np.square(c)) current_energies.append(energy) total np.sum(current_energies) current_norm current_energies / total # 计算每个子带相对基线的变化幅度 delta current_norm - np.array(baseline_energies) # 返回变化最大的子带索引和对应频带描述 max_idx np.argmax(np.abs(delta)) return max_idx, delta # 例出现第 2 层细节系数能量 30%而 1、3 层基本没变 # 配合轴承外圈故障特征频率落在 D2 频带可初步定位轴承外圈这段代码背后的逻辑是能量变化方向比单个频带绝对能量更重要。设备负载增加会让所有子带能量同步上升但只有一个子带显著变化时才更可能是局部缺陷。因此根因定位环节不要只看能量绝对值一定要带健康基线做比值。基线怎么建平台上线后连续采集 7 天正常工况信号按每个测点算平均归一化能量动态更新是另一个话题基线过老会导致误定位一般每 30 天要重算一次。5.2 五个避坑记录现象、原因、解决下面五条是本系统从原型到实盘最常见的坑每条都来自实际项目回访。按“现象 → 原因 → 解决”写方便你直接对照排查。坑一周期性识别把 50Hz 工频当成了设备周期现象是平台输出“设备主周期 50Hz”而现场设备实际转频是 25Hz诊断结论完全偏离机械真实状态。原因是传感器信号里混入了电磁工频干扰且幅值高于机械振动基频在 HPS 谱上 50Hz 的乘积能量盖过了真实基频。解决分两步入口处对 50Hz 做窄带陷波或用高通滤波把 30Hz 以下保留、50Hz 处挖掉同时用额定转频表校验候选频率凡是不能和转频或轴承特征频率对上的直接丢弃并标记“疑似电源干扰”。我后来在所有测点配置里都加了干扰频率列表这个坑基本不再出现。坑二小波分解层数过深导致重构信号过度平滑异常分数虚高现象是设备高负荷运行时平台频繁报出 3 级以上异常现场人员疲于确认但设备离线检查没有问题。原因是 level 取到 7 且阈值过大重构信号丢失了正常工况下应有的能量波动残差里混入的是正常结构噪声。解决方法是把 level 限制在 45 层同时把 tail 阈值从 1.5 回退到 1.01.2 之间直到正常工况下残差分布接近预期的正态形态再冻结参数。坑三滑窗长度太短周期性识别结果漂移不定现象是同一设备连续两次计算出的基频相差 0.8Hz无法做为告警基线。原因是窗口长度只有 2048 点频率分辨率等于 12.5Hz远高于转子转频间隔任何泄漏都会把峰值推偏。解决是把窗口加长到至少 1 秒频率分辨率降到 1Hz 以内当被分析的转频只有 8Hz 时窗口要拉到 5 秒才能稳定识别这也要在拿到测点真实转速后重新核算。坑四核密度估计带宽太小正常信号被识别为异常现象是训练好的 KDE 模型在正常工况下也频繁告警。原因是带宽取太小残差分布被拟合出大量尖峰每个训练样本都被突出为一个独立的密度峰新样本只要稍偏离训练点就被推到尾部。解决是改用 silverman 经验带宽或按 1.06σ 乘 n 的负五分之一次方手动计算每次更新训练样本集后也要重新检查带宽不要沿用旧参数。坑五根因定位时只看频带能量绝对值负载波动引发误定位现象是负载升高时定位模块输出“轴承磨损”实际只是进料量增大导致整体能量抬升。原因是把负载导致的整体能量变化误判成了局部缺陷信号。解决是强制使用归一化能量差同时引入转速和负荷信号作为辅助条件如果平台有 PLC 负荷输入识别出负荷波动和振动同步上升时应抑制根因定位输出并标记“负荷相关变化”。6. 进阶做法用历史故障样本回放验证整个链路与阈值调参顺序这一章分享一个我坚持了很久的习惯每次新接入一台设备我都会先找两三段明确故障历史的原始信号哪怕是拆机后确认的那种把信号切成一个标准化回放测试集喂给整套流程验证平稳性判断、周期识别、异常程度转换和根因定位能不能复现已知结论。这个回放测试集不做算法训练只做链路自检。具体做法是先把已标记的故障时间段单独分割保存输入系统后检查三步平稳性标记是否在故障段前后被置为非平稳异常程度分数是否在故障发生附近到达 3 级以上根因定位的子带编号是否落在已知部件特征频带内。三项里任何一项对不上就先重新检查前一个模块的参数而不是继续往下调。这个习惯能帮我快速锁定问题发生在哪一层。曾经有一次回放轴承内圈故障样本时异常分数和根因定位都正确但平稳性标记几乎没有触发追溯发现故障早期是非常微弱的高频调制频谱相似度还在阈值之上。这种情况的处理不是放宽平稳性阈值而是把周期性识别改到细节系数的包络谱上让早期冲击更容易暴露。回放测试集让这类调试变得可重复而不是每次等故障再临场分析。阈值调参的顺序我也固定在同一个框架里先调平稳性相似度阈值再调 HPS 乘积层数接着调 KDE 带宽和 EWMA 平滑系数最后才调根因定位能量差阈值。这个顺序严格依赖数据流的前后关系前面的阈值没定住后面的所有异常和定位结论都会跟着漂移。很多人喜欢一上来就调告警分数阈值结果后面每次更换设备都要整体重调等于进入一个无限循环。还有一个小技巧值得养成习惯所有阈值参数在部署时写进配置文件并带上来源备注标明“稳态数据 30 天统计”或“回放样本 X 确定”不要直接塞在代码里。我见过太多平台上线半年后没人能说清 0.9 和 1.2 这两个阈值当初是怎么定的改动参数全靠试探那是最痛苦的一种翻车方式。把参数来源和回放测试集一起管理相当于给阈值买了颗后悔药。这套结合傅里叶全局扫描与小波局部抓取的运维数据分析链路配上核密度估计的残差打分和频带能量反推已经在我经手的几个工业设备监控项目里稳定运行。不同设备的特征频率表、采样率、负荷特性都会有差异但链路骨架是可迁移的。希望你按这个思路搭起来之后能少踩几个我当年踩过的坑也希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/31 10:33:30 · 阅读 1.2k 小波变换与傅里叶变换融合的智能运维时序数据异常检测与故障定位
简介一套面向工业设备监控与故障诊断的智能运维时序数据分析系统基于小波变换、傅里叶变换与核密度估计实现波形平稳性检测、周期性识别、异常程度转换与根因定位。适用于工业数据分析、设备健康管理和故障预测方向的开发者与研究者可帮助快速理解并搭建AI运维分析流程。资源共17个文件包含10个Python脚本、4个CSV数据集及说明文档整体约8.79MB脚本覆盖信号预处理、平稳性/周期性判断、异常检测与根因定位等核心环节CSV数据便于直接复现实验。目前已有94人学习下载。通过源码与配套说明既能学习小波多尺度分解与频域特征提取的具体实现也能掌握核密度估计在异常定位中的应用方式系统从数据读取、特征分析到异常评估给出了完整链路适合作为毕业设计、课题研究或实际工业监控项目的参考方案。1. 基于小波变换和傅里叶变换的智能运维时序数据分析系统解决的不只是“有没有故障”我在工业设备监控项目里见过太多这样的误判离心泵的振动频谱看起来完全正常转频分量也不突出但拆开轴承才发现已经严重磨损另一台电机电流波形到了某个时间段开始变“毛糙”FFT 峰值几乎没变化事后检查发现定子绝缘已经击穿。问题恰恰出在把傅里叶变换当成唯一的频谱工具平稳信号里那些瞬态成分被彻底摊平了。这套基于小波变换和傅里叶变换的智能运维时序数据分析系统核心做法是让两种变换分工傅里叶做全局周期识别和平稳性判断小波抓取局部突变与瞬时冲击再用核密度估计把重构残差映射成可告警的异常程度分数最后沿着频带能量的变化做根因定位。它面向的是设备状态监控工程师、故障诊断算法开发和 AI 运维平台建设者解决的是“看到异常但说不清异常从哪来”的落地问题。2. 傅里叶变换与小波变换的分工全局周期与局部突变的选型逻辑2.1 傅里叶变换在设备监控中能看见什么、看不见什么先明确一个基础事实傅里叶变换把一段时序信号整体变换到频域输出的是这段时长内各个频率分量的幅值和相位。对稳定运行的设备比如恒转速风机、齿轮箱、离心机频谱里基频和边频带是有规律的我通常会再叠加一层谐波乘积谱 HPS 来放大基频把“设备有没有周期”从人工看图变成自动判定。但代价也很直接一次 FFT 的频谱是整段窗口内的统计平均时间分辨率等于零。窗口里如果发生过一次短暂冲击它在频谱上的能量会被摊薄到整个观察时长里只在很宽的频段上留下一点不太起眼的凸起肉眼几乎无法识别。import numpy as np from scipy.fft import rfft, rfftfreq # 现场振动信号采样率 25600 Hz取 0.5 秒作为一次窗口 fs 25600 win_size int(fs * 0.5) # 12800 点 raw get_sensor_segment(device_idpump_01, start_ts..., lengthwin_size) # 去均值消除直流分量对低频段的压制 x raw - np.mean(raw) # 加汉宁窗限制截断处不连续产生的频谱泄漏 xw x * np.hanning(win_size) # 单边幅值谱rfft 只输出一半谱乘 2 恢复真实幅值 amp np.abs(rfft(xw)) * 2.0 / win_size freq rfftfreq(win_size, 1 / fs) # 切除 3 Hz 以下的极低频段只保留有物理意义的部分 keep freq 3 amp amp[keep] freq freq[keep]这段代码里有三个参数值得单独说明。第一个是加窗现场信号不可能正好截断在整数周期上端点不连续会让一个频率的能量泄漏到周围频段加汉宁窗是成本最低的缓解手段。第二个是频率分辨率它等于 fs 除以 win_size当前参数下分辨率是 2 Hz意味着两个间隔小于 2 Hz 的谱峰在频谱上只会显示成一个鼓包你要单独检测 0.5 Hz 的差异就必须把窗口拉长到 2 秒以上。第三个是单边谱的归一化rfft 输出的是复数频谱对实数信号只需要前半部分乘 2 再除以窗口长度才是物理上可读的幅值。这套参数是周期性识别和后续平稳性判断的共同起点换采样率时一定要同步改窗口长度不能只改其中一个。2.2 小波变换补上的局部时频能力与基函数选择小波变换不把信号拆成无限长的正弦波叠加而是用一组可伸缩、可平移的小波基去做内积输出不同尺度带上的系数。在设备监控场景下这相当于把信号按频带切开同时保留时间位置信息。细节系数对应高频瞬态近似系数对应低频趋势。我在这类系统里通常选 db4 或 db6 小波原因说得直白一点带点玄学成分但实测多年确实有效Daubechies 系小波的波形形态和机械冲击响应有相似性冲击能量能集中在少数几个系数上不像 Symlets 或 Coiflets 那样把一次冲击摊到好几层系数中给后续阈值判断带来额外麻烦。import pywt def decompose_signal(x, waveletdb4, level4): # wavedec 返回的元组第一位是近似系数后面依次是各层细节系数 coeffs pywt.wavedec(x, wavelet, modesymmetric, levellevel) names [fA{level}] [fD{i} for i in range(level, 0, -1)] return dict(zip(names, coeffs)) coeffs decompose_signal(x, waveletdb4, level4) # 采样率 fs25600Hz 时4 层分解各频带大致对应 # A4 约 0800HzD4 约 8001600Hz # D3 约 16003200HzD2 约 32006400HzD1 约 640012800Hz选层数时有几个边界要自己算清楚。每层系数长度减半D1 长度是原始信号的一半A4 长度约为窗口点数的十六分之一如果窗口本身只有 2048 点A4 只剩 128 点统计意义就很弱了。另一个是 Nyquist 频率小波每一层细节系数的频带是上一层的一半但原始采样率决定了最高只能识别到 fs/2也就是 12800 Hz。如果你关心的轴承故障特征频率落在 8001600 HzD4 就对应得很准如果关心的是 50100 Hz 的转频倍频就要把 level 加到 6 或 7或者直接回到傅里叶频谱去做。2.3 两种变换在一个系统里怎么串成管线我在这类监控平台上常把数据流拆成三个动作先用傅里叶做全局扫描用滑窗 FFT 算频谱一致性和周期主峰一旦平稳性判断触发「非平稳」标记就把后续特征计算切换到小波分解最后异常程度转换阶段再用小波重构残差作为输入。原因是两块计算成本差异很大FFT 一次窗口只做一组蝶形运算而小波分解要逐层卷积计算量大约是前者的 510 倍全部数据都跑小波在边缘网关上是浪费。只有被判定为存在瞬态或非平稳特征的数据段才值得用小波系数进一步分析。这样组织也方便加载历史故障样本做回放验证后面第六章会细说。下表是我在不同设备类型上的选型经验供你直接套用设备工况主要现象首选分析方式关注频带恒定转速轴承磨损周期冲击脉冲小波细节系数包络8003200Hz齿轮箱啮合异常边频带变化傅里叶高阶频谱啮合频率附近转子不平衡转频幅值稳定上升傅里叶基频趋势1倍转频变速工况下的喘振频率随时间漂移滑窗 FFT 频谱瀑布全频段观察漂移电机定子匝间短路电流波形毛糙小波细节系数能量突增高频噪声段3. 平稳性检测和周期性识别滑窗频谱相似度与谐波乘积谱3.1 平稳性不是统计书上的平稳而是频谱一致性的工程近似教科书里的平稳性定义是联合概率分布不随时间平移这个条件在工业现场几乎不可能满足也没有工程意义。我在智能运维平台里实际用的是另一套标准连续多个滑窗的幅度谱互相之间是否一致。设备平稳运转时每个窗口的频谱形状应该高度相似轴承出现局部剥落、叶片发生裂纹等瞬态源加入后某个窗口里会突然多出一段高频能量相邻窗口频谱相似度随即掉下去。这种近似虽然不像 ADF 检验那样有严格统计学证明但胜在能和频谱分析共用同一套滑窗机制计算开销极小。def spectral_similarity(amp_series): # 输入是多个窗口的幅值谱形状为 (num_windows, num_freq_bins) sims [] for i in range(1, len(amp_series)): a amp_series[i - 1] b amp_series[i] dot np.dot(a, b) norm_a np.linalg.norm(a) norm_b np.linalg.norm(b) # 加一个极小值避免除以 0 sims.append(dot / (norm_a * norm_b 1e-12)) return np.array(sims) # 假设已经用滑窗 FFT 拿到了 60 个窗口的幅值谱 # 窗口长度 0.5 秒重叠 50%覆盖 15 秒数据 window_size int(fs * 0.5) hop_size int(window_size * 0.5) amp_series [] # 每个元素是长度为窗口内的幅值谱数组 for start in range(0, len(x) - window_size, hop_size): seg x[start:start window_size] seg seg * np.hanning(window_size) amp_series.append(np.abs(rfft(seg))[: window_size // 2 1]) sims spectral_similarity(amp_series) # 平稳阈值SIM_TH 0.90 # 如果 sims 中连续 3 个值低于 0.90标记该时段为非平稳阈值 0.90 怎么定出来的多台设备离线数据回放后我发现正常定工况下频谱相似度通常在 0.93 以上偶发电源干扰会造成单点下降到 0.85但下一窗口立刻恢复真正进入故障发展期时相似度是持续走低并稳定在 0.85 以下。所以这里的判断逻辑不是看单点而是看连续低于阈值的窗口数量。如果你处理的信号是电流而非振动建议把阈值放宽到 0.88因为电流波形受负荷波动影响更大频谱变化比振动更频繁。3.2 周期性识别谐波乘积谱与倍频校验的组合手段周期性识别的目标是回答两个问题设备运行周期是多少这个周期是不是真实的机械周期。单靠 FFT 找峰值不够因为频谱里响应最大的是某次倍频而不是基频。我习惯在原始幅值谱基础上再加一层谐波乘积谱 HPS把谱分别下采样 1/2、1/3、1/4然后逐点相乘真实基频会在所有压缩谱中同时保留并叠加放大谐波分量则被相互错位抵消。def hps_spectrum(amp, num_products4): # amp: 原始幅值谱长度为 N max_idx len(amp) // num_products hps np.ones(max_idx) for k in range(1, num_products 1): # 每级下采样 k 倍截断到统一长度 compressed amp[::k][:max_idx] hps * compressed return hps # 峰值定位取 hps 中前 5 个局部极大值 peak_indices np.argsort(hps)[-5:][::-1] candidate_freqs freq[peak_indices] # 进一步用设备额定转频校验取最接近真实转速倍频的候选HPS 层数选 4 还是 5取决于频谱长度。num_products 越大基频权重越突出但有效频率范围会被压缩到原频谱的 1/5高频段信息丢失。我一般先限制检测范围在 3 Hz2000 Hz再确认 num_products 不大于 5。峰值定位后必须做一次倍频校验就是把候选频率和额定转频、叶片通过频率、轴承特征频率表对比能对上才写入周期特征库这种做法能把电源工频干扰误判成设备周期的问题挡在门外。3.3 平稳性结果怎么驱动小波分析接管平稳性检测结果不能只停留在曲线里要直接参与数据流选择。平台在实盘运行时按固定长度缓存原始信号每缓存一个窗口就做一次频谱相似度更新当相似度持续下降并触发非平稳标记时系统自动把这 1530 秒的原始缓存发送到小波分析模块而不是删除缓存继续等下一个窗口。这里有个常被忽略的工程细节小波分解需要足够长的序列来保证低层近似系数的统计意义如果触发发现得太晚缓存里只剩 2 秒信号分解到 A5 层几乎就剩不了几个点。我的做法是缓存区保持 60 秒、滑窗 0.5 秒标记触发时至少还有 30 秒数据可供小波分析。# 伪代码级的状态管理逻辑 STATE STABLE similarity_trend [] def on_window(amp_series): global STATE s spectral_similarity([amp_series[-2], amp_series[-1]])[0] similarity_trend.append(s) if len(similarity_trend) 10: similarity_trend.pop(0) # 连续 3 个窗口低于阈值才切换状态避免单点干扰 if all(v 0.90 for v in similarity_trend[-3:]): STATE UNSTABLE elif all(v 0.93 for v in similarity_trend[-3:]): STATE STABLE这里连续 3 个窗口的意义是消除偶然扰动。之前踩过一个坑只用单窗口相似度判定每次电压波动都会让平台误报“进入非平稳”小波分析模块被频繁唤醒网关 CPU 直接跑满。改成连续 3 次判定后误报率大幅下降这个参数也成了我以后做类似平台的默认值。4. 异常程度转换从重构残差到核密度估计的分数映射4.1 为什么要用重构残差定义异常程度直接拿振动幅值 RMS 做异常分级是最常见但最不可靠的做法。设备负载不同同一 RMS 值对应的健康状态完全不同同样幅值的信号在正常工况下可能只是负荷波动在重载工况下却可能是早期故障信号。我用的思路是先建立“正常模型”用小波重构逼近期望中的正常信号然后拿实际信号和期望信号的差当异常特征。这个差值叫重构残差幅值大、持续长才说明设备偏离正常状态比直接看绝对幅值靠谱得多。小波重构的做法是先对细节系数做阈值收缩把代表噪声的小系数置零或压缩再用处理后的系数反变换得到平滑信号。阈值的选择直接影响残差质量阈值太高重构信号过于干净残差里全是正常负荷波动后续异常程度会被夸大阈值太低阈值收缩等于没做残差和原始信号几乎相同。我一般用软阈值加自适应尺度参数尺度由分解层数和信号噪声水平共同决定。def wavelet_reconstruct_smooth(signal, waveletdb4, level4, scale1.2): coeffs pywt.wavedec(signal, wavelet, modesymmetric, levellevel) # 对每一层细节系数做软阈值收缩 thresholded [coeffs[0]] # 近似系数不动 for i, detail in enumerate(coeffs[1:]): sigma np.median(np.abs(detail)) / 0.6745 # 噪声方差估计 thresh scale * sigma * np.sqrt(2 * np.log(len(detail))) # 软阈值绝对值小于阈值的置零其余向零压缩 thresholded.append(pywt.threshold(detail, thresh, modesoft)) return pywt.waverec(thresholded, wavelet, modesymmetric) smooth wavelet_reconstruct_smooth(x, level4, scale1.2) residual x - smooth这段代码的核心在噪声方差估计 sigma 和阈值系数 scale。sigma 用中位数绝对偏差除以 0.6745是为了在受冲击信号影响时仍然稳健scale 是留给你的现场调节旋钮我通常在 1.0 到 1.5 之间调scale 越大重构信号越平滑检测灵敏度越高但误报风险也上升。注意 pywt.threshold 的 mode 用 soft 而不是 hard软阈值能避免重构信号在阈值点产生人工跳跃这个跳跃如果存在会直接形成伪冲击残差。4.2 核密度估计给残差建立分布模型拿到了残差序列下一步是回答“这个残差算不算异常”。最简单的方式是假设残差服从正态分布用均值加减 3 倍标准差做界限但残差往往带重尾和偏态特别是故障初期的瞬态冲击会让残差分布出现明显长尾。此时用核密度估计 KDE 更稳它不对分布形状做强假设通过把每个样本点变成一个小核函数并求和逼近真实概率密度。scipy.stats.gaussian_kde 提供现成实现核心要调的是带宽参数 bw_method。from scipy.stats import gaussian_kde # 用健康设备的历史正常残差训练 KDE # residual_train 是一维数组长度为 n kde_model gaussian_kde(residual_train, bw_methodsilverman) # 同时计算正常残差的 99% 分位数作为硬边界防止极端值干扰 abs_residual np.abs(residual_train) hard_limit np.percentile(abs_residual, 99) def residual_to_anomaly_score(residual_value, kde_model, hard_limit): # 求当前残差值在正常分布下的概率密度 density kde_model.evaluate([abs(residual_value)])[0] # 密度越低越位于分布尾部异常程度越高 # 用相对密度转成分数密度高于某个值记 0否则按低于比例分级 normal_density_99 kde_model.evaluate([hard_limit])[0] if abs(residual_value) hard_limit: return 0 if density normal_density_99 * 0.2: return 4 elif density normal_density_99 * 0.5: return 3 elif density normal_density_99: return 2 return 1这个分数映射里99% 分位数是工程上的安全垫。gaussian_kde 的带宽选 silverman 时是按样本量自动推导的适合大多数振动残差如果样本量小于 200建议手动传一个标量带宽值比如 0.05带宽太小会过拟合样本尖峰太大则会抹平双峰结构这是调整这个模型最需要耐心的地方。gaussian_kde.evaluate 输入要二维数组形式不写方括号直接传标量会报维度错误这是很多人第一次跑通代码时最容易卡住的地方。4.3 从密度到分数异常程度分级与告警联动异常程度转换最后要落到监控平台的告警联动。单独看一个窗口的分数意义有限因为瞬态冲击和电磁干扰也能造成单点高分但实际损害未必严重。我习惯把窗口分数送进一个 EWMA指数加权移动平均滤波器让分数平滑把瞬时抖动压下去只保留持续上升的趋势平滑系数 0.3 时一次孤立的冲击在经过 5 个窗口后影响几乎归零而持续上升的故障特征会稳定在一个高位。def ewma_smooth(scores, alpha0.3): smoothed [] # 平滑后的分数序列 prev scores[0] for s in scores: prev alpha * s (1 - alpha) * prev smoothed.append(prev) return np.array(smoothed) smoothed_scores ewma_smooth(raw_scores, alpha0.3) # 分级联动 # 01 级: 正常 2 级: 关注进入趋势跟踪 # 3 级: 预警通知现场检查 4 级: 停机建议平滑系数 alpha 是分级输出的关键参数。alpha 太大平滑几乎不起作用太小真实故障趋势要到很晚才浮现容易错过维护窗口。我一般把 alpha 控制在 0.250.35配合 10 秒一个输出周期使用如果你做的是慢变故障比如轴承润滑退化可以降到 0.1让趋势更稳定。分级阈值本身也不是固定不变的需要在平台上线后用历史故障样本回放校验这个落到第六章展开。5. 智能运维根因定位的排查与避坑频带映射、谐波误判与阈值漂移5.1 根因定位的事后链路异常时刻的频带能量反推异常程度分数触发预警后根因定位要回答“哪个部件出的问题”。我用的办法是频带能量反推原始信号经过小波分解后各层细节系数对应不同频带当异常分数升高时对比异常时段和正常基线的各频带能量占比找出能量变化最剧烈的频带再对照设备已知的特征频率表判断物理部件位置。轴承外圈故障的特征频率通常落在某个具体频邻齿轮啮合故障则集中在啮合频率及其边带电机转子断条表现为转频侧边带各有固定位置不需要靠猜。def locate_fault_subband(coeffs, baseline_energies): # coeffs: 小波分解系数baseline_energies 是健康时段各个子带的归一化能量 current_energies [] for c in coeffs: energy np.sum(np.square(c)) current_energies.append(energy) total np.sum(current_energies) current_norm current_energies / total # 计算每个子带相对基线的变化幅度 delta current_norm - np.array(baseline_energies) # 返回变化最大的子带索引和对应频带描述 max_idx np.argmax(np.abs(delta)) return max_idx, delta # 例出现第 2 层细节系数能量 30%而 1、3 层基本没变 # 配合轴承外圈故障特征频率落在 D2 频带可初步定位轴承外圈这段代码背后的逻辑是能量变化方向比单个频带绝对能量更重要。设备负载增加会让所有子带能量同步上升但只有一个子带显著变化时才更可能是局部缺陷。因此根因定位环节不要只看能量绝对值一定要带健康基线做比值。基线怎么建平台上线后连续采集 7 天正常工况信号按每个测点算平均归一化能量动态更新是另一个话题基线过老会导致误定位一般每 30 天要重算一次。5.2 五个避坑记录现象、原因、解决下面五条是本系统从原型到实盘最常见的坑每条都来自实际项目回访。按“现象 → 原因 → 解决”写方便你直接对照排查。坑一周期性识别把 50Hz 工频当成了设备周期现象是平台输出“设备主周期 50Hz”而现场设备实际转频是 25Hz诊断结论完全偏离机械真实状态。原因是传感器信号里混入了电磁工频干扰且幅值高于机械振动基频在 HPS 谱上 50Hz 的乘积能量盖过了真实基频。解决分两步入口处对 50Hz 做窄带陷波或用高通滤波把 30Hz 以下保留、50Hz 处挖掉同时用额定转频表校验候选频率凡是不能和转频或轴承特征频率对上的直接丢弃并标记“疑似电源干扰”。我后来在所有测点配置里都加了干扰频率列表这个坑基本不再出现。坑二小波分解层数过深导致重构信号过度平滑异常分数虚高现象是设备高负荷运行时平台频繁报出 3 级以上异常现场人员疲于确认但设备离线检查没有问题。原因是 level 取到 7 且阈值过大重构信号丢失了正常工况下应有的能量波动残差里混入的是正常结构噪声。解决方法是把 level 限制在 45 层同时把 tail 阈值从 1.5 回退到 1.01.2 之间直到正常工况下残差分布接近预期的正态形态再冻结参数。坑三滑窗长度太短周期性识别结果漂移不定现象是同一设备连续两次计算出的基频相差 0.8Hz无法做为告警基线。原因是窗口长度只有 2048 点频率分辨率等于 12.5Hz远高于转子转频间隔任何泄漏都会把峰值推偏。解决是把窗口加长到至少 1 秒频率分辨率降到 1Hz 以内当被分析的转频只有 8Hz 时窗口要拉到 5 秒才能稳定识别这也要在拿到测点真实转速后重新核算。坑四核密度估计带宽太小正常信号被识别为异常现象是训练好的 KDE 模型在正常工况下也频繁告警。原因是带宽取太小残差分布被拟合出大量尖峰每个训练样本都被突出为一个独立的密度峰新样本只要稍偏离训练点就被推到尾部。解决是改用 silverman 经验带宽或按 1.06σ 乘 n 的负五分之一次方手动计算每次更新训练样本集后也要重新检查带宽不要沿用旧参数。坑五根因定位时只看频带能量绝对值负载波动引发误定位现象是负载升高时定位模块输出“轴承磨损”实际只是进料量增大导致整体能量抬升。原因是把负载导致的整体能量变化误判成了局部缺陷信号。解决是强制使用归一化能量差同时引入转速和负荷信号作为辅助条件如果平台有 PLC 负荷输入识别出负荷波动和振动同步上升时应抑制根因定位输出并标记“负荷相关变化”。6. 进阶做法用历史故障样本回放验证整个链路与阈值调参顺序这一章分享一个我坚持了很久的习惯每次新接入一台设备我都会先找两三段明确故障历史的原始信号哪怕是拆机后确认的那种把信号切成一个标准化回放测试集喂给整套流程验证平稳性判断、周期识别、异常程度转换和根因定位能不能复现已知结论。这个回放测试集不做算法训练只做链路自检。具体做法是先把已标记的故障时间段单独分割保存输入系统后检查三步平稳性标记是否在故障段前后被置为非平稳异常程度分数是否在故障发生附近到达 3 级以上根因定位的子带编号是否落在已知部件特征频带内。三项里任何一项对不上就先重新检查前一个模块的参数而不是继续往下调。这个习惯能帮我快速锁定问题发生在哪一层。曾经有一次回放轴承内圈故障样本时异常分数和根因定位都正确但平稳性标记几乎没有触发追溯发现故障早期是非常微弱的高频调制频谱相似度还在阈值之上。这种情况的处理不是放宽平稳性阈值而是把周期性识别改到细节系数的包络谱上让早期冲击更容易暴露。回放测试集让这类调试变得可重复而不是每次等故障再临场分析。阈值调参的顺序我也固定在同一个框架里先调平稳性相似度阈值再调 HPS 乘积层数接着调 KDE 带宽和 EWMA 平滑系数最后才调根因定位能量差阈值。这个顺序严格依赖数据流的前后关系前面的阈值没定住后面的所有异常和定位结论都会跟着漂移。很多人喜欢一上来就调告警分数阈值结果后面每次更换设备都要整体重调等于进入一个无限循环。还有一个小技巧值得养成习惯所有阈值参数在部署时写进配置文件并带上来源备注标明“稳态数据 30 天统计”或“回放样本 X 确定”不要直接塞在代码里。我见过太多平台上线半年后没人能说清 0.9 和 1.2 这两个阈值当初是怎么定的改动参数全靠试探那是最痛苦的一种翻车方式。把参数来源和回放测试集一起管理相当于给阈值买了颗后悔药。这套结合傅里叶全局扫描与小波局部抓取的运维数据分析链路配上核密度估计的残差打分和频带能量反推已经在我经手的几个工业设备监控项目里稳定运行。不同设备的特征频率表、采样率、负荷特性都会有差异但链路骨架是可迁移的。希望你按这个思路搭起来之后能少踩几个我当年踩过的坑也希望帮到你。本文还有配套的精品资源点击获取