嵌入式Linux下Modbus RTU串口开发实战：从配置到轮询
2026/9/11 3:57:12
简介：针对轴承全寿命数据分析与故障诊断需求，这份资源提供了一套完整的时域与频域特征提取方案，覆盖17个时域统计特征和13个频域统计特征，适用于机械设备健康监测、剩余寿命预测等研究方向，也适合初涉特征工程的科研人员对照学习。包内包含11个文件，以9个MATLAB脚本（.m）为主，辅以2个自动保存文件（.asv），压缩包整体仅12KB。脚本模块划分清晰，时域统计与频域统计计算均有独立实现，另封装了希尔伯特解调、ZOOMFFT细化频谱等分析功能，可支撑从原始信号到特征向量的完整流程。目前已有1738人学习下载，说明该工具脚本在轴承数据挖掘场景中具有不错的参考价值。借助这些脚本，读者可快速复现常见时频特征提取逻辑，理解统计指标与包络谱、ZOOMFFT等方法的实现细节，并迁移到自己的数据集上进行二次开发。
1. 轴承全寿命数据时域频域特征提取：先用 30 个统计量把退化过程量化
做轴承全寿命（run-to-failure）退化预测，第一步通常不是建模，而是把原始振动信号压缩成一组能长期跟踪的特征。17 个时域特征加 13 个频域特征，是故障诊断里沿用多年的标准特征集，覆盖幅值、能量、冲击、频谱形状四个维度。RMS 和能量负责整体退化趋势，峭度、峰值因子、裕度因子负责早期点蚀冲击，重心频率等频域矩反映频谱重心随磨损的迁移。没有哪个单特征能走完全寿命，但 30 个组合起来，退化阶段、报警阈值、剩余寿命回归就都有据可依。下面把公式、实现、参数边界和筛选用法一次讲清，适合正在处理全寿命振动数据的工程师对照落地。
2. 17 个时域特征：公式清单、批量计算代码与无量纲因子的数值坑
2.1 先对齐 17 个时域特征：有量纲统计量与无量纲因子
不同论文里的“17 个时域特征”成员并不完全一致，按最常见的划分方式，前 10 个是有量纲统计量，后 7 个是无量纲因子和分布高阶矩。这里先给全清单，后面代码和章节都按这个编号来。
| 编号 | 特征名 | 公式 | 对退化的响应 |
|---|---|---|---|
| T1 | 均值 mean | mean(x) | 静态偏置，退化中变化小 |
| T2 | 绝对平均值 | mean(|x|) | 幅值水平，趋势与 RMS 接近 |
| T3 | 方差 var | var(x) | 波动能量 |
| T4 | 标准差 std | sqrt(var(x)) | 波动幅度 |
| T5 | 均方根 RMS | sqrt(mean(x²)) | 整体振动能量，全寿命主趋势 |
| T6 | 方根幅值 | (mean(sqrt(|x|)))² | 幅值水平，对孤立冲击不敏感 |
| T7 | 峰值 peak | max(|x|) | 冲击幅值 |
| T8 | 峰峰值 | max(x) − min(x) | 最大摆幅 |
| T9 | 最大值 | max(x) | 单向摆幅上限 |
| T10 | 最小值 | min(x) | 单向摆幅下限 |
| T11 | 偏度 skewness | E[(x−μ)³]/σ³ | 分布对称性 |
| T12 | 峭度 kurtosis | E[(x−μ)⁴]/σ⁴ | 早期冲击最敏感 |
| T13 | 峰值因子 | peak / RMS | 冲击相对能量 |
| T14 | 波形因子 | RMS / 绝对平均值 | 波形形状变化 |
| T15 | 脉冲因子 | peak / 绝对平均值 | 冲击相对幅值 |
| T16 | 裕度因子 | peak / 方根幅值 | 冲击敏感，比峰值因子更稳当 |
| T17 | 能量 | sum(x²) | 总能量，与 RMS 强相关 |
表格里 x 表示窗口内的一维振动序列，n 为窗口长度，μ 和 σ 分别是均值和标准差。有量纲特征依赖传感器灵敏度和安装位置，跨测点对比前要先归一化；T13～T16 四个无量纲因子消掉了幅值量纲，在载荷波动工况下表现更一致，这是它们能在早期故障里派上用场的原因之一。
2.2 用 NumPy 和 SciPy 批量计算 17 个时域特征
import numpy as np from scipy import stats def extract_td_features(x): """输入一段一维振动信号 x，返回 17 个时域特征字典。 参数 x : np.ndarray，长度建议 >= 1024，过短时偏度和峭度估计会失真。 """ x = np.asarray(x, dtype=np.float64) rms = np.sqrt(np.mean(x ** 2)) # T5 均方根 abs_mean = np.mean(np.abs(x)) # T2 绝对平均值 peak = np.max(np.abs(x)) # T7 峰值 sqrt_amp = np.mean(np.sqrt(np.abs(x))) ** 2 # T6 方根幅值 return { 'td_mean': np.mean(x), 'td_abs_mean': abs_mean, 'td_var': np.var(x), 'td_std': np.std(x), 'td_rms': rms, 'td_sqrt_amp': sqrt_amp, 'td_peak': peak, 'td_pk_pk': np.max(x) - np.min(x), 'td_max': np.max(x), 'td_min': np.min(x), 'td_skew': stats.skew(x), # T11 偏度 'td_kurt': stats.kurtosis(x, fisher=False), # T12 峭度，正态分布=3 'td_crest': peak / (rms + 1e-12), # T13 峰值因子 'td_shape': rms / (abs_mean + 1e-12), # T14 波形因子 'td_impulse': peak / (abs_mean + 1e-12), # T15 脉冲因子 'td_clearance': peak / (sqrt_amp + 1e-12), # T16 裕度因子 'td_energy': np.sum(x ** 2), # T17 能量 }几个参数和定义的取舍要讲清楚。np.var和np.std默认算总体值（除以 n），如果对标某些论文里的样本标准差（除以 n−1），小窗口下会差几个百分点，全寿命特征序列里必须固定一种口径。stats.skew返回的是有偏估计，stats.kurtosis(fisher=False)返回 Pearson 峭度，正态分布等于 3；很多文献直接用超额峭度（正态等于 0），跨文章对比前必须先确认。T13～T16 四个因子都要做分母保护，RMS、绝对平均值、方根幅值在轴承刚启动或传感器静默时可能接近 0，直接相除会出现极大的孤立尖峰，1e-12的作用是把异常值限制在合理范围，而不是掩盖真实冲击。
2.3 短窗口下的失真和峭度回落的误读
窗口长度小于 100 点时，偏度和峭度的估计方差会急剧变大，甚至返回 nan。全寿命切片时窗口至少要覆盖几百个采样点，工程上建议 4096 点起步，这不仅是为了统计可靠性，也是为了下一章频域特征的分辨率。另一个常见误读是峭度曲线：早期点蚀出现时峭度先冲高，随后剥落面扩大、冲击变成连续宽带振动，峭度反而回落——这不是退化结束，而是故障形态从“稀疏冲击”变成了“密集宽带”。所以峭度要配 RMS 一起看，一个看整体趋势，一个看冲击形态。
3. 13 个频域特征：幅值谱约定、频率矩公式与 Python 落地
3.1 频域特征统计的对象是幅值谱
13 个频域特征本质上是对频谱再做一次统计。第一步要定死统计对象：用幅值谱还是功率谱。幅值谱和时域幅值量纲一致，边频带结构看得直观；功率谱对幅值做了平方，会放大主频带的权重。我一般固定用单边幅值谱 S(k) 定义全部 13 个特征，P1～P5 和 P10～P13 都是对谱线的统计，P6～P9 是频率的一阶和二阶矩。常见论文里也有用功率谱算 P1～P6 的版本，数值不同但方向一致，关键是同一个项目里只允许一种定义。口径统一后，退化导致的频谱形态变化才能被这几个矩稳定捕捉。
3.2 13 个频域特征的公式与物理含义
| 编号 | 特征 | 公式 | 含义 |
|---|---|---|---|
| P1 | 幅值谱均值 | mean(S) | 谱线整体水平 |
| P2 | 幅值谱方差 | var(S) | 谱线离散程度 |
| P3 | 幅值谱标准差 | std(S) | 谱线离散程度 |
| P4 | 幅值谱偏度 | skew(S) | 谱分布的非对称性 |
| P5 | 幅值谱峭度 | kurt(S) | 谱峰尖锐程度 |
| P6 | 重心频率 FC | Σf·S / ΣS | 频谱重心位置（Hz） |
| P7 | 均方频率 MSF | Σf²·S / ΣS | 频率二阶矩（Hz²） |
| P8 | 频率方差 VF | MSF − FC² | 谱线围绕重心的分散度（Hz²） |
| P9 | 均方根频率 RMSF | sqrt(MSF) | 特征频率（Hz） |
| P10 | 谱能量 | ΣS² | 频谱总能量 |
| P11 | 峰值频率 | argmax(S) 对应 f | 主谱线频率 |
| P12 | 峰值幅值 | max(S) | 主谱线高度 |
| P13 | 谱总幅值 | ΣS | 幅值谱面积 |
P6～P9 四个频率矩对退化很敏感：磨损导致宽带能量上升时，高频段权重抬升，P6 和 P8 的上升往往早于 RMS 的明显变化。P11 和 P12 直接给出主频成分的位置，适合和转频、故障特征频率的理论值叠图对比。
3.3 Python 实现：加窗、单边谱与频率矩一条龙
import numpy as np from scipy import stats def extract_fd_features(x, fs): """输入一段一维振动信号 x 和采样率 fs，返回 13 个频域特征字典。""" x = np.asarray(x, dtype=np.float64) x = x - np.mean(x) # 去直流，避免 0 Hz 分量主导谱统计量 n = len(x) win = np.hanning(n) # Hann 窗抑制频谱泄露 xw = x * win X = np.fft.rfft(xw, n=n) # 单边 FFT，点数 n//2+1 amp = np.abs(X) / np.sum(win) * 2.0 # 幅值谱，补偿窗能量后单边乘 2 amp[0] = amp[0] / 2.0 # 直流分量的 2 要退回去 freq = np.fft.rfftfreq(n, d=1.0 / fs) S = amp denom = np.sum(S) + 1e-12 fc = np.sum(freq * S) / denom # P6 重心频率 msf = np.sum(freq ** 2 * S) / denom # P7 均方频率 vf = max(msf - fc ** 2, 0.0) # P8 频率方差，钳位到 0 peak_idx = np.argmax(S) return { 'fd_amp_mean': np.mean(S), # P1 'fd_amp_var': np.var(S), # P2 'fd_amp_std': np.std(S), # P3 'fd_amp_skew': stats.skew(S), # P4 'fd_amp_kurt': stats.kurtosis(S, fisher=False), # P5 'fd_fc': fc, # P6 'fd_msf': msf, # P7 'fd_vf': vf, # P8 'fd_rmsf': np.sqrt(msf), # P9 'fd_spec_energy': np.sum(S ** 2), # P10 'fd_peak_freq': freq[peak_idx], # P11 'fd_peak_amp': S[peak_idx], # P12 'fd_spec_sum': np.sum(S), # P13 }参数说明：去均值这步不能省，轴承信号里若有微小直流量，0 Hz 谱线会直接抬升 P1 和 P13，掩盖真实频谱形状。np.hanning的代价是主瓣展宽约一倍，换来旁瓣明显压低，适合做统计量，但不适合精确测单根谱线幅值。幅值标定中np.sum(win)是窗函数能量补偿，乘以 2 是把负频率能量并入正频，直流分量不乘 2。P8 用max(..., 0.0)钳位，因为浮点误差可能让 MSF − FC² 出现微小的负值，没有物理意义。
3.4 频率分辨率由窗口长度决定，边频带是检验尺
很多人习惯拿整段信号做一次 FFT，但全寿命特征提取里谱是逐窗口算的，频率分辨率 df = fs / n，n 是窗口长度而不是整段信号长度。以 25.6 kHz 采样率为例：窗口 4096 点时 df 是 6.25 Hz；窗口降到 1024 点，df 变成 25 Hz。滚动轴承外圈故障的调制边带间隔等于转频（约 33 Hz），25 Hz 的分辨率会把相邻边带直接糊在一起，P11 会在多个谱线间来回跳。窗口低于 2048 点时要警惕这类频域失真；反过来窗口太长又会拖慢特征序列的时间分辨率，报警变迟钝。4096～8192 点是全寿命数据常用的折中区间。
4. 全寿命特征序列：滑动窗口切片、特征矩阵构建与退化曲线
4.1 窗口长度和步长：按转频与退化速度两个约束定
滑动窗口参数是特征提取里最容易被忽略的环节。约束一是窗口至少覆盖 5～10 个转频周期，统计量才对相位不敏感；约束二是步长决定特征序列的时间分辨率，一般取窗口的 25%～50%。对一组典型的全寿命数据，按下面的参数起步基本不会跑偏。
| 工况参数 | 取值 | 依据 |
|---|---|---|
| 采样率 fs | 25600 Hz | 加速度传感器常见采集配置 |
| 转频 fr | 约 33 Hz | 对应 2000 RPM 工况 |
| 单圈采样点数 | 约 776 | 25600 / 33 |
| 窗口长度 | 4096 | 约 5.3 圈，统计上足够可靠 |
| 步长 | 1024 | 相邻窗口 75% 重叠，趋势平滑 |
| 频率分辨率 | 6.25 Hz | 25600 / 4096，可分辨 33 Hz 边带 |
有些全寿命数据集按固定间隔落盘（比如每 10 分钟一个文件），这时不需要滑动窗口，直接对每个文件算 30 个特征就是一条特征序列；滑动窗口用在连续采集的单一大文件上。两种做法的特征矩阵结构相同，后续退化趋势处理也一致。
4.2 批量切片并组装 30 维特征矩阵
import pandas as pd def build_feature_frame(x, fs, win_len=4096, step=1024): """把整段全寿命信号切成滑动窗口，返回 (窗口数, 30) 的特征 DataFrame。""" rows = [] n = len(x) for start in range(0, n - win_len + 1, step): seg = x[start:start + win_len] row = {'t_start': start / fs} # 窗口起始时间（秒） row.update(extract_td_features(seg)) row.update(extract_fd_features(seg, fs)) rows.append(row) return pd.DataFrame(rows) features_df = build_feature_frame(signal, fs=25600) print(features_df.shape) # (约900000, 31)，31 = 时间列 + 30 个特征时间开销主要在stats.skew和stats.kurtosis的纯 Python 循环上：10 小时、9000 万点的信号会切出约 90 万个窗口，单线程跑完大概几十分钟。常见做法是用multiprocessing或numba把两个特征函数并行化，或者用np.lib.stride_tricks.sliding_window_view生成窗口矩阵再向量化。内存上 90 万行 × 31 列 float64 约 220 MB，可以接受；内存紧张时就把特征行按块写入 parquet，不必一次全部驻留。
4.3 退化曲线长什么样：三条典型模式
特征序列生成后，先画三张图再谈建模。RMS 曲线通常呈现“平稳段—缓升段—剧升段”的三段式，对应正常磨损、裂纹扩展和严重剥落；峭度曲线的形态是早期冲高、中期回落的倒 V 型，这是冲击从稀疏到连续的过渡；重心频率 FC 和 RMSF 正常段基本平直，退化后期向高频漂移。如果出现异常——比如 RMS 初期就持续上涨、峭度全程无波动、FC 开局就往高频跑——优先检查传感器量程饱和、安装松动和工频干扰，不要急着改特征定义。曲线形态正常，特征才算真正进入可用状态。
5. 特征筛选与早期报警验证：两个指标加一条切分纪律
5.1 用单调性和趋势相关过滤 30 个特征
退化特征要进趋势模型，基本要求是单调。两个指标能快速完成初筛：
import numpy as np def monotonicity(s): """单调性，0~1，越接近 1 越单调。""" d = np.diff(np.asarray(s, dtype=np.float64)) return abs((d > 0).sum() - (d < 0).sum()) / len(d) def trend_corr(s): """与时间轴的线性相关，给出方向和强度。""" t = np.arange(len(s)) return np.corrcoef(t, s)[0, 1]对 30 列特征逐列计算，单调性低于 0.3 的（偏度、波形因子这类来回摆的特征）直接淘汰；trend_corr给出正负号，负相关表示特征在退化中下降，同样是有效信息。按经验，RMS、能量、峰峰值、谱能量的单调性通常排在最前，峭度的单调性反而不高，因为它中期回落。这两个指标跨数据集时排序基本一致，可以直接作为筛选依据。
5.2 早期报警点对比与时序切分纪律
把全寿命按 5% 分段，对每一段求峭度和 RMS 的均值，取前 5～10 段做正常基线，首次越限（基线均值 + 3σ）的时刻就是该特征的报警点。峭度通常比 RMS 提早 10%～20% 寿命报警，这正是早期点蚀的窗口期，也是把 30 个特征全部算出来的意义所在。建模划分训练集和验证集时，全寿命数据严禁随机打乱，要用 TimeSeriesSplit 保持时间顺序；标准化参数只从训练集拟合，验证集按同一套均值和方差变换。这一步漏掉，前 30 个特征里量级大的几个会把模型带偏，比任何调参问题都致命。
本文还有配套的精品资源，点击获取