从时域到频域:用 DolphinDB 做工业振动信号的故障特征提取

📅 2026/7/10 11:54:02
从时域到频域:用 DolphinDB 做工业振动信号的故障特征提取
摘要有一句行业黑话振动数据里藏着设备的体检报告。一台电机、一台泵、一台风机的轴承、齿轮、转子一旦出现早期损伤振动信号里往往比温度、电流更早露出马脚。但问题在于振动信号不像温度那样看一眼数值就知道好坏。它的价值藏在波形里需要从时域走到频域才能看清楚——而这恰恰是传统时序数据库最不擅长的事数据存在库里分析得导到 MATLAB 或 Python 里跑 FFT算完再导回来链路一长实时性就没了。这篇就围绕 DolphinDB把工业振动信号分析的完整链路在库内跑一遍从时域统计特征RMS、峭度、峰值因子到频域 FFT 幅度谱再到工况对齐与轴承特征频率诊断。一、为什么振动分析必须从时域走到频域先把这件事讲清楚不然后面的代码就只是一堆函数调用。时域分析看的是波形本身——某个时刻振动有多大、整体波动剧烈不剧烈。它的指标RMS、峰值、峭度适合判断这台设备是不是出了问题但很难告诉你问题出在哪里。而且时域指标对早期故障不敏感轴承上一个针眼大小的点蚀在时域 RMS 上几乎看不出来但设备其实已经开始劣化了。频域分析看的是信号里各个频率成分的能量分布。机械设备的故障往往有明确的频率特征不平衡 → 1 倍转频1×能量升高不对中 → 2 倍转频2×明显轴承外圈点蚀 → 外圈故障特征频率BPFO处出现峰值齿轮磨损 → 啮合频率及其边带换句话说时域告诉你有事发生频域告诉你是什么事。这就是为什么振动诊断离不开 FFT。而 DolphinDB 的价值在于在 signal 插件中集成了功能强大且性能高的 FFTW3 开源傅里叶变换库如fft、ifft、小波变换等信号处理函数并且内置了fft、ifft、小波变换等信号处理函数且全部向量化优化过。这意味着整个时域→频域→诊断的链路可以在数据库内一行脚本完成不用把 GB 级的原始波形导出去再导回来。二、数据准备高频振动表与工况表振动分析需要两类数据高频振动波形和设备工况转速、负载。后者很重要因为同一个振动幅值在不同转速下含义完全不同——这点后面会展开。// 振动表高频原始波形按日期设备复合分区 db1 database(, VALUE, 2026.01.01..2026.12.31) db2 database(, HASH, [SYMBOL, 20]) db database(dfs://vib, COMPO, [db1, db2]) ​ vibSchema table(1:0, tsdeviceIdxAxisyAxiszAxis, [TIMESTAMP, SYMBOL, DOUBLE, DOUBLE, DOUBLE, DOUBLE]) db.createPartitionedTable(vibSchema, vibration, tsdeviceId) ​ // 工况表低频转速/负载1Hz 采样 condSchema table(1:0, tsdeviceIdrpmloadstatus, [TIMESTAMP, SYMBOL, DOUBLE, DOUBLE, INT]) db.createPartitionedTable(condSchema, condition, tsdeviceId)振动表存三轴加速度工业现场常用 ICP 振动传感器输出采样率假设 10kHz。工况表存转速和负载1Hz 就够。两张表采样频率差了一万倍——这正是后面要用asof join对齐的原因。三、时域特征一行脚本出体检报告时域特征是振动分析的第一道筛子。常用的几个指标每个对故障的敏感度不一样指标物理含义对什么故障敏感RMS均方根振动整体能量中晚期磨损、不平衡峰值瞬时最大冲击冲击类故障、剥落峰值因子Crest Factor峰值 / RMS早期轴承冲击敏感峭度Kurtosis分布的尖峭程度早期轴承故障最敏感歪度Skewness分布的偏斜不对称磨损经验上峭度和峰值因子对早期轴承故障最敏感——RMS 还没明显变化时峭度往往已经跳起来了。这也是为什么工业现场做预警峭度几乎是必算的。在 DolphinDB 里这些指标一行脚本 就能算出来// 提取每台设备在指定时段的时域特征 select deviceId, // RMS sqrt(mean(x^2))三轴合成 sqrt(avg(xAxis*xAxis yAxis*yAxis zAxis*zAxis)) as rms, // 峰值合成幅值 max(sqrt(xAxis*xAxis yAxis*yAxis zAxis*zAxis)) as peak, // 峰值因子 max(sqrt(xAxis*xAxis yAxis*yAxis zAxis*zAxis)) / sqrt(avg(xAxis*xAxis yAxis*yAxis zAxis*zAxis)) as crestFactor, // 峭度用合成幅值的峭度 kurtosis(sqrt(xAxis*xAxis yAxis*yAxis zAxis*zAxis)) as kurtosis, // 歪度 skew(sqrt(xAxis*xAxis yAxis*yAxis zAxis*zAxis)) as skewness from loadTable(dfs://vib, vibration) where date(ts) 2026.06.15 group by deviceId几个工程细节值得说一下kurtosis()和skew()是 DolphinDB 内置函数直接对一组值算高阶矩不用自己写公式。RMS 用sqrt(avg(x*x))而不是sqrt(avg(pow(x,2)))是因为x*x在向量化下比pow更快——这种细节在大数据量下会累积成可观的差距。这些聚合全部在存储节点上完成存算一体不需要把波形数据搬到计算节点。3000 万行数据算这五个指标单节点耗时大概在秒级。算完之后按峭度降序排排在前面的设备就是重点怀疑对象。但到这一步你还不知道它们为什么异常——这就需要频域分析了。四、频域分析FFT 提取幅度谱把一台怀疑异常设备的振动波形取出来做 FFT 看频谱// 取 DEV_003 在 10:00~10:01 这一分钟的 Z 轴波形 wave exec zAxis from loadTable(dfs://vib, vibration) where deviceId DEV_003 and ts between datetime(2026.06.15 10:00:00) and datetime(2026.06.15 10:01:00) order by ts ​ // 做 FFT取幅度谱 sampleRate 10000 // 10kHz n size(wave) spectrum abs(fft(wave)) // 幅度谱复数取模 halfN n \ 2 // 实信号频谱对称取前一半 halfSpec spectrum[0 : halfN] ​ // 频率轴 freqs (0 .. halfN-1) * (sampleRate \ n) ​ // 找主频幅度最大的频率分量 peakFreqIdx imax(halfSpec) peakFreq peakFreqIdx * (sampleRate \ n)几个关键点exec ... order by ts返回的是一个有序向量FFT 要求输入是等间隔采样的时域序列所以必须保证时间顺序。fft()返回复数向量abs()取模得到幅度谱。实信号的频谱关于 Nyquist 频率对称所以只看前一半。imax()返回最大值的索引乘以频率分辨率sampleRate / n就得到主频。这是定位哪个频率出问题的关键。拿到主频后对照设备的转频就能初步判断故障类型。比如这台 DEV_003 标称转速 3000 RPM转频 50Hz如果幅度谱在 50Hz 处有强峰可能是转子不平衡在 100Hz2×转频有强峰可能是不对中。但这里有个坑转速不是恒定的。设备实际运行转速会随负载波动固定按 50Hz 去找峰值在变工况场景下会错得离谱。这就引出了下一节。五、工况对齐用 asof join 解决变工况下的误报工业现场很少有设备恒速运行。一台电机的转速可能在 2400~3000 RPM 之间随负载波动。如果你用固定转频做频域诊断转速一变特征频率就漂移诊断结果全乱。解决办法是把振动分析和工况数据关联起来——按当前转速动态计算特征频率。这正是asof join的用武之地// 把每条振动数据和当时最近的工况转速对齐 select v.ts, v.deviceId, v.zAxis, c.rpm, c.load, c.status from aj( loadTable(dfs://vib, vibration) as v, loadTable(dfs://vib, condition) as c, deviceIdts ) where v.deviceId DEV_003 and v.ts between datetime(2026.06.15 10:00:00) and datetime(2026.06.15 10:01:00)ajasof join的逻辑是对振动表里的每条高频记录在工况表里找同一设备、时间不晚于该记录的最近一条转速数据。这样不用对低频工况数据做重采样、也不丢失高频波形的细节直接在原始层面完成对齐。对齐之后就可以按转速分段做频域分析。比如把转速 2800~2900 RPM 的振动片段挑出来单独做 FFT再在这个片段里按实际转频约 47.5Hz去找特征频率。这种按工况分桶的诊断方式比固定转频可靠得多。一个完整的按转速分桶诊断查询大概是这样// 把转速分桶每 50 RPM 一档每桶算时域特征 select floor(c.rpm \ 50) * 50 as rpmBucket, v.deviceId, sqrt(avg(v.zAxis * v.zAxis)) as rms, kurtosis(v.zAxis) as kurtosis, count(*) as samples from aj( loadTable(dfs://vib, vibration) as v, loadTable(dfs://vib, condition) as c, deviceIdts) where v.deviceId DEV_003 and date(v.ts) 2026.06.15 group by floor(c.rpm \ 50) * 50, v.deviceId order by rpmBucket结果会告诉你这台设备在不同转速下的 RMS 和峭度分别是多少。正常设备的峭度在各转速桶里应该比较平稳如果某个转速段峭度突然跳高那多半是共振或该转速下激发的故障——这是固定转频诊断根本发现不了的。六、轴承特征频率把频谱和机械几何对应起来走到这一步可以上更硬核的了轴承故障特征频率。滚动轴承的故障频率不是随机的它由轴承的几何参数和转速决定有标准公式外圈故障频率 BPFO N/2 × 1 \- d/D × cosα × f_r内圈故障频率 BPFI N/2 × 1 \ d/D × cosα × f_r滚动体故障频率 BSF D/2d × 1 \- \(d/D × cosα²) × f_r保持架故障频率 FTF 1/2 × 1 \- d/D × cosα × f_r其中 N 是滚动体数量d 是滚动体直径D 是节圆直径α 是接触角f_r 是转频rpm/60。知道轴承型号这些几何参数就查得到。结合第五节对齐到的实际转速可以算出当前工况下的特征频率再去幅度谱里核对这些频率处有没有异常峰值// 轴承参数示例SKF 6205 深沟球轴承N9, d7.94mm, D39.04mm, α0 def bearingFreqs(rpm) { N 9; d 7.94; D 39.04; alpha 0 fr rpm \ 60.0 cosA cos(alpha) return dict(BPFOBPFIBSFFTF, [ N\2 * (1 - d\D * cosA) * fr, N\2 * (1 d\D * cosA) * fr, D\(2*d) * (1 - pow(d\D * cosA, 2)) * fr, 1\2 * (1 - d\D * cosA) * fr ]) } ​ // 结合当前转速算特征频率 bearingFreqs(2975) // 输出类似BPFO107.4, BPFI148.7, BSF139.2, FTF11.9然后在幅度谱里看这些频率附近有没有能量峰值。一个实用做法是先对原始信号做包络分析demodulation把高频谐振带的能量搬到低频再做 FFT 得到包络谱——早期轴承故障在包络谱里比在原始频谱里清晰得多。包络分析的标准实现是解析信号法对信号做 FFT → 把负频率置零构造解析信号→ IFFT → 取模得到包络 → 再 FFT 得到包络谱。这套流程在 DolphinDB 里用fft/ifft/abs完全可以实现因为它本质就是 FFT 的组合操作。算出包络谱后在 BPFO107.4Hz处如果出现明显峰值及其谐波基本就可以判定外圈故障了。需要说明的是包络分析对参数带通中心频率、带宽比较敏感不同设备的谐振频带不同往往要先做一段频谱分析确定谐振带再算包络。这是工程经验问题不是算力问题——而在库内完成的好处恰恰是这些试错过程不需要把数据搬来搬去。七、从批量到实时时域特征上流频域走定时前面的分析都是基于历史数据的批量分析。但生产现场还需要实时性——传感器数据一进来就得判断这台设备是不是开始恶化。这里有个工程上务实的划分时域特征RMS、峭度走流计算计算量小、对延迟敏感适合挂到响应式状态引擎做亚毫秒级实时预警。频域分析FFT、包络谱走定时批量FFT 计算量大且通常需要一段完整的波形才有意义不适合逐点流式处理。实际做法是定时比如每 10 秒取一段波形做一次 FFT。时域特征的实时化很简单定义一个带状态的峭度预警函数挂到响应式状态引擎// 定义有状态的峭度预警函数 state def vibAlert(zAxis, windowSize, kurtThreshold) { kurt mstd(zAxis, windowSize) ... // 简化示意实际峭度需维护四阶矩 return iif(kurt kurtThreshold, 1, 0) } ​ // 建流表接收实时振动数据 share streamTable(1:0, tsdeviceIdzAxis, [TIMESTAMP, SYMBOL, DOUBLE]) as vibStream ​ // 输出告警表 alerts table(1:0, tsdeviceIdzAxisisAnomaly, [TIMESTAMP, SYMBOL, DOUBLE, INT]) ​ // 响应式状态引擎每条数据进来都更新窗口状态 createReactiveStateEngine( namevibMonitor, metrics[zAxis, vibAlert(zAxis, 1000, 4.0)], dummyTablevibStream, outputTablealerts, keyColumndeviceId ) ​ subscribeTable(tableNamevibStream, actionNamemonitor, handlertableInsert{getStreamEngine(vibMonitor)})注意state注解——它声明这是一个有状态函数引擎会自动维护窗口状态增量计算On→O1而不是每来一个点就把窗口里的数据重算一遍。这就是为什么峭度这种高阶统计量也能做到实时。频域的实时监控则用定时任务每隔 N 秒取最近一段波形跑 FFT在特征频率处看能量超阈值就告警。这种时域实时 频域准实时的分层架构在工程上比所有分析都走流更现实。八、局限与适用场景客观说库内做信号处理不是银弹有几个边界要清楚不是替代专业诊断软件。成熟的企业级振动诊断系统像 SKFs ptitude、Bently Nevada里有大量专家规则库和案例库DolphinDB 提供的是算力底座和函数库帮你把分析逻辑跑起来但诊断规则还得靠领域专家沉淀。FFT 假设信号平稳。对强非平稳信号启停机过程、冲击FFT 会有频谱泄漏需要配合 STFT 或小波变换。DolphinDB 提供了这些函数但用对用错还是要懂信号处理基础。包络分析等高级流程需要自己组合。DolphinDB 提供的是 fft/ifft/小波这些原子函数包络分析、阶次分析这类组合流程要自己用脚本拼。好处是灵活代价是要对信号处理流程比较熟。轻量场景没必要。如果只是几台风机看个振动幅值超限直接用阈值报警就行没必要在库里跑 FFT。库内信号处理的价值在设备多、要算的特征多、要做频域诊断的规模化场景。小结把这条链路串起来看高频振动波形 工况数据 │ ├── 时域特征RMS/峭度/峰值因子── 一行脚本秒级 │ └── 挂流计算引擎 → 实时预警 │ ├── 频域分析FFT 幅度谱── 定位故障频率 │ └── 配合转速动态算特征频率 │ └── 工况对齐asof join── 解决变工况误报 └── 按转速分桶诊断核心就几句话时域看有没有问题用峭度和峰值因子捕捉早期异常频域看是什么问题用 FFT 把特征频率找出来工况对齐是前提用 asof join 让诊断跟得上转速波动时域走流、频域走定时是工程上务实的实时化分层。振动分析这件事过去之所以难落地不是因为算法有多复杂FFT 中学物理就学过而是因为数据搬运成本太高——波形存在 A 系统转速在 B 系统FFT 要导到 C 系统跑结果再回到 D 系统告警。链路一长实时性、可靠性、开发效率全下来了。DolphinDB 把这条链路收拢到一个库里价值不在于它发明了什么新算法而在于它让你不用再为搬运数据折腾。对一个要在几百台设备上规模化部署振动诊断的团队来说这件事的意义远大于几个百分点的算法精度差异。参考链接DolphinDB 技术文档DolphinDB 内置函数fft / kurtosis / imax 等DolphinDB 流数据处理DolphinDB asof join 多频数据对齐