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

原创

Xxtaoaooo

发布于 2026-07-08 22:14:00

120

摘要

有一句行业黑话："振动数据里藏着设备的体检报告"。一台电机、一台泵、一台风机的轴承、齿轮、转子一旦出现早期损伤，振动信号里往往比温度、电流更早露出马脚。
但问题在于，振动信号不像温度那样"看一眼数值就知道好坏"。它的价值藏在波形里，需要从时域走到频域才能看清楚——而这恰恰是传统时序数据库最不擅长的事：数据存在库里，分析得导到 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, `ts`deviceId`xAxis`yAxis`zAxis,
                  [TIMESTAMP, SYMBOL, DOUBLE, DOUBLE, DOUBLE, DOUBLE])
db.createPartitionedTable(vibSchema, "vibration", `ts`deviceId)

// 工况表：低频转速/负载，1Hz 采样
condSchema = table(1:0, `ts`deviceId`rpm`load`status,
                   [TIMESTAMP, SYMBOL, DOUBLE, DOUBLE, INT])
db.createPartitionedTable(condSchema, "condition", `ts`deviceId)

振动表存三轴加速度（工业现场常用 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 处有强峰，可能是转子不平衡；在 100Hz（2×转频）有强峰，可能是不对中。

但这里有个坑：转速不是恒定的。设备实际运行转速会随负载波动，固定按 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,
    `deviceId`ts
)
where v.deviceId = "DEV_003"
  and v.ts between datetime(2026.06.15 10:00:00) and datetime(2026.06.15 10:01:00)

aj（asof 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,
    `deviceId`ts)
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 深沟球轴承，N=9, d=7.94mm, D=39.04mm, α=0）
def bearingFreqs(rpm) {
    N = 9; d = 7.94; D = 39.04; alpha = 0
    fr = rpm \ 60.0
    cosA = cos(alpha)
    return dict(`BPFO`BPFI`BSF`FTF, [
        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)
// 输出类似：BPFO=107.4, BPFI=148.7, BSF=139.2, FTF=11.9

然后在幅度谱里看这些频率附近有没有能量峰值。一个实用做法是：先对原始信号做包络分析（demodulation），把高频谐振带的能量"搬"到低频，再做 FFT 得到包络谱——早期轴承故障在包络谱里比在原始频谱里清晰得多。

包络分析的标准实现是解析信号法：对信号做 FFT → 把负频率置零（构造解析信号）→ IFFT → 取模得到包络 → 再 FFT 得到包络谱。这套流程在 DolphinDB 里用 fft / ifft / abs 完全可以实现，因为它本质就是 FFT 的组合操作。算出包络谱后，在 BPFO（107.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, `ts`deviceId`zAxis, [TIMESTAMP, SYMBOL, DOUBLE]) as vibStream

// 输出告警表
alerts = table(1:0, `ts`deviceId`zAxis`isAnomaly,
               [TIMESTAMP, SYMBOL, DOUBLE, INT])

// 响应式状态引擎：每条数据进来都更新窗口状态
createReactiveStateEngine(
    name="vibMonitor",
    metrics=<[zAxis, vibAlert(zAxis, 1000, 4.0)]>,
    dummyTable=vibStream,
    outputTable=alerts,
    keyColumn="deviceId"
)

subscribeTable(tableName="vibStream", actionName="monitor",
              handler=tableInsert{getStreamEngine("vibMonitor")})

注意 @state 注解——它声明这是一个有状态函数，引擎会自动维护窗口状态（增量计算，O(n)→O(1)），而不是每来一个点就把窗口里的数据重算一遍。这就是为什么峭度这种高阶统计量也能做到实时。

频域的实时监控则用定时任务：每隔 N 秒取最近一段波形，跑 FFT，在特征频率处看能量，超阈值就告警。这种"时域实时 + 频域准实时"的分层架构，在工程上比"所有分析都走流"更现实。

八、局限与适用场景

客观说，库内做信号处理不是银弹，有几个边界要清楚：

不是替代专业诊断软件。成熟的企业级振动诊断系统（像 SKFs @ptitude、Bently Nevada）里有大量专家规则库和案例库，DolphinDB 提供的是"算力底座和函数库"，帮你把分析逻辑跑起来，但诊断规则还得靠领域专家沉淀。
FFT 假设信号平稳。对强非平稳信号（启停机过程、冲击），FFT 会有频谱泄漏，需要配合 STFT 或小波变换。DolphinDB 提供了这些函数，但用对用错还是要懂信号处理基础。
包络分析等高级流程需要自己组合。DolphinDB 提供的是 fft/ifft/小波这些原子函数，包络分析、阶次分析这类组合流程要自己用脚本拼。好处是灵活，代价是要对信号处理流程比较熟。
轻量场景没必要。如果只是几台风机看个振动幅值超限，直接用阈值报警就行，没必要在库里跑 FFT。库内信号处理的价值在"设备多、要算的特征多、要做频域诊断"的规模化场景。

小结

把这条链路串起来看：

高频振动波形 + 工况数据
        │
        ├── 时域特征（RMS/峭度/峰值因子）── 一行脚本，秒级
        │         └── 挂流计算引擎 → 实时预警
        │
        ├── 频域分析（FFT 幅度谱）── 定位故障频率
        │         └── 配合转速动态算特征频率
        │
        └── 工况对齐（asof join）── 解决变工况误报
                  └── 按转速分桶诊断

核心就几句话：