电机数据分析全流程：从zip解压到FFT异常检测

简介：面向电机数据分析与风电系统研究的MATLAB脚本资源，压缩包内包含一个nan_wt05.m脚本，可用于感应双馈发电机的数据预处理、特征提取与仿真建模。脚本重点引入Relief特征选择算法，通过计算分类权重，帮助研究者从电压、电流、转速等多项运行参数中筛选关键特征，支撑故障诊断或效率优化等工作，也适用于亚同步/超同步等不同运行状态下的数据对比。压缩包共1个文件，类型为m脚本，大小仅6KB，轻量易用，适合已有MATLAB基础、正在开展电机或新能源发电数据分析的工程师与学生参考，目前已有128人学习。通过学习此脚本，可了解完整的电机数据分析流程，包括缺失值处理、特征权重计算及双馈发电机运行状态仿真思路，为后续搭建自己的分析模型或风电系统优化提供直接参考。

1. 拿到nan_wt05.zip先别急着解压：电机数据分析的第一步不是pandas

电机数据分析这活，很多人拿到nan_wt05.zip第一反应是unzip -o解压，再import pandas开读，这个顺序其实反了。nan_wt05大概率对应测试工位或电机批次，zip里可能几十个CSV，也可能一两个超大日志，采样率从1Hz到50kHz混在一起。电机数据跟行情数据不一样，转速、电流、温度、振动强耦合，少个物理量或时间轴对不齐，后边特征工程全白做。我一般先做三件事：校验压缩包完整性、建字段级物理量清单、核对采样率与时间轴。三步走完才轮到pandas。常见事故是把单位字符串混进数值列，DataFrame类型静默变成object，到建模阶段才炸。先盘点再解析，是从nan_wt05.zip这类电机数据包开始的正确姿势。

2. 解压与数据探查：从nan_wt05.zip的文件清单到电机参数字段映射

2.1 先校验压缩包结构，再决定是unzip还是流式读取

电机数据包通常不是几十兆的小文件。48小时连续采集，振动2kHz、电流1kHz，一个CSV文本能逼近2GB。直接把zip全量解压，探查阶段就浪费好几分钟磁盘IO。我习惯先用unzip -l nan_wt05.zip看压缩包目录，或者用Python的zipfile模块列出文件清单和压缩比，从而判断里面是单个大文件分卷，还是按日期拆分的多个CSV。

import zipfile ZIP_PATH = "nan_wt05.zip" with zipfile.ZipFile(ZIP_PATH, "r") as zf: for info in zf.infolist(): if not info.is_dir(): ratio = info.file_size / info.compress_size if info.compress_size else 0 print(f"{info.filename}\t{info.file_size/1024:.0f} KB\t" f"压缩比 {ratio:.1f}x")

这段代码不解压任何文件内容，只读取zip中央目录，把文件名、原始大小、压缩比列出来。电机数据的CSV按文本存储，压缩比通常在5～20倍之间；如果某个文件压缩比不到2倍，很可能内部已经是二进制或内嵌压缩格式，后面读取方式要区分对待。文件清单还能暴露非数据文件，比如.ini、.log、.xlsx，这些偶尔携带台架号、电机型号、环境温度等元信息，先记下来，后面做数据校验时能对得上。

命名规律在这一步基本能判断出来。ch1_20231101_100000.csv这种带通道名和起始时间戳的，说明记录程序按通道分文件；Rec_001.csv这种顺序编号的，说明按下电或实验阶段切段。前者时间拼接靠文件名解析，后者往往需要从数据内部的时间列再对齐。

2.2 字段体检与三分类映射：波形量、统计量、工况量

电机数据字段通常能按物理含义分成三类：

分类方法不用人工逐列看，用我常用的体检函数跑一遍：

import pandas as pd def inspect_motor_csv(path: str, nrows: int = 200) -> None: df = pd.read_csv(path, nrows=nrows, encoding_errors="replace", skipinitialspace=True) for col in df.columns: uniq = df[col].nunique(dropna=True) sample = df[col].dropna().iloc[0] if uniq else None print(f"{col:24s} dtype={str(df[col].dtype):12s} " f"unique={uniq:8d} first={str(sample)[:20]}") inspect_motor_csv("ch1_20231101_100000.csv")

参数说明：

• nrows=200只抽样前200行，探查场景下足够判断字段类型
• skipinitialspace=True专门处理数采软件导出CSV时逗号后带空格的坏习惯，不加这个参数，数值列会被解析成字符串并产生隐性缺失
• encoding_errors="replace"防止BOM或编码异常让读取直接中断
• unique值很少的字段基本是Mode或状态码，可以用astype("category")压缩存储

体检结果里如果发现"Current(A)"这种列名，顺手用正则清洗成Current_A。真正危险的信号是电流列打印出来是object类型且sample值形如" 12.3"或"12.3 A"，这种需要在读取阶段就做定制转换，而不是事后遍历全表替换。事后替换在大数据量下慢且容易漏。

2.3 时间轴对齐：把高采样率通道重采样到统一基准

分文件存储的电机数据，最麻烦的是不同通道采样率不在一个量级。振动是2kHz、电流是1kHz、温度是1Hz。对齐的常规做法是以最低频信号为基准，把高频信号降采样聚合。

df_vib = df_vib.set_index("Time").resample("1S").agg( vibr_rms=lambda x: float((x**2).mean() ** 0.5) ) df_cur = df_cur.set_index("Time").resample("1S").mean() df_tmp = df_tmp.set_index("Time").resample("1S").mean() df_aligned = df_vib.join(df_cur, rsuffix="_cur").join(df_tmp, rsuffix="_tmp")

这里有个电机数据特有的坑：振动信号降采样不能用mean()。振动幅值正负对称，均值接近0，直接平均等于把振幅信息全部抹掉。正确聚合方式是先做RMS，也就是对窗口内所有样本求平方和均值再开方，代码里用lambda嵌套实现。而电流可以先平均，因为电流均值对应电磁转矩水平，是有物理意义的量。温度本身变化慢，1Hz重采样已经足够。

提示：如果后续要做FFT分析，时间对齐后的1s重采样数据已经丢失了振动信号的高频细节。FFT必须在原始2kHz或更高采样率的数据上做，重采样数据只看趋势和建立回归模型。这条边界要提前跟团队对齐，否则后面特征和标签对不上。

3. 清洗与特征构造：把nan_wt05.zip的原始信号变成可建模的电机状态量

3.1 缺失值处理：停机和通信中断不能用前向填充

电机数据里的NaN，含义比普通时序数据复杂得多。缺几个点是传感器毛刺；缺几秒是通信丢包；缺几分钟往往意味着电机进入停机状态或数采系统重新缓存。如果一股脑fillna(method="ffill")，等于把停机区硬写成上一段运行值，不仅制造假数据，还会让后续均值、方差特征全部失真。我处理缺失值的固定套路是先按缺失段长度分级。

import numpy as np def flag_missing_blocks(s: pd.Series, short: int = 5, long: int = 50): is_na = s.isna() grp = (is_na != is_na.shift()).cumsum() na_len = is_na.groupby(grp).transform("sum").where(is_na, 0) s_fill = s.interpolate(limit_area="inside", limit=short) s_fill = s_fill.fillna(0) mask_short = (na_len > 0) & (na_len <= short) mask_long = na_len > long out = pd.DataFrame({ "value": s_fill, "state": np.select([mask_short, mask_long], [1, 2], default=0) }) return out

代码逻辑说明：

• 用isna().cumsum()把连续缺失区间分组，这是处理缺失段的通用trick，每变化一次就产生新的段号
• interpolate(limit_area="inside", limit=short)只允许插值填补缺失段内部，且最多补5个点，超过就放弃
• 缺失段长度超过50个采样点，直接标记state=2，对应停机或断连状态
• 返回的state列作为后续模型的辅助特征，比单纯把缺失值去掉多保留了一层因果信息

这样处理之后，停机段不会参与稳态工况统计，但能用于标记数据覆盖范围。

3.2 滑动窗口构造RMS、峰值与温度变化率特征

清洗完原始量后，下一步是构造时序特征。电机数据分析最常用的三类滑窗特征和参数配置如下：

滑窗特征不建议全部用pandas的rolling().mean()做，数据量上来了会很慢。我常用numpy的索引矩阵构造窗口，这里给一个numpy实现：

import numpy as np def sliding_rms(x: np.ndarray, window: int, step: int) -> np.ndarray: n = len(x) out_len = 1 + (n - window) // step idx = np.arange(window)[None, :] + step * np.arange(out_len)[:, None] win = x[idx] # shape: (out_len, window) return np.sqrt(np.mean(win ** 2, axis=1)) vibr_rms = sliding_rms(vibration_signal, 2000, 100)

参数说明：

• window=2000对应2kHz采样率下的1秒物理窗口，能量意义上的秒级RMS
• step=100是0.05秒滑动一次，兼顾时间分辨率和计算量
• 通过np.arange构造索引矩阵，比逐窗口切片快一个数量级，处理几百万点数据也不会有明显卡顿

3.3 工况段划分：把非平稳信号切成稳态区间

电机不是始终在一个转速下运行，起停、加减速、加载卸载之间电压电流特性完全不同。全序列统一建阈值模型一定会误报。我一般的做法是根据转速设定值和转矩反馈做工况切分，只保留稳态段做特征统计。

speed_set = df_aligned["Speed_Set"] mode_change = (speed_set.diff().abs() > 50).astype(int) segment_id = mode_change.cumsum() df_aligned["seg"] = segment_id seg_stats = df_aligned.groupby("seg").agg( dur_s=("Time", "size"), speed_mean=("Speed_Set", "mean"), curr_rms=("Current_A", lambda x: float(np.sqrt(np.mean(x**2)))) ) steady_segs = seg_stats[seg_stats["dur_s"] >= 30]

逻辑说明：

• diff().abs()>50：转速设定值相邻采样变化超过50rpm视为工况切换
• mode_change.cumsum()把每次切换生成一个新段号
• 稳态段定义是至少连续30秒没有切换，组内统计才有统计学意义
• 稳态段提取出来后，FFT、回归、阈值统计都只在这些段上计算，能规避加减速过程瞬时电流冲击对整体特征的干扰

4. 频谱分析与异常检测：在nan_wt05.zip上做FFT实操与阈值设定

4.1 抽取稳态段做FFT：窗函数与频率分辨率的取舍

电机振动信号的FFT不是直接对整段数据做变换就完事。常见做法是从稳态段里抽取连续2048点或4096点，加窗后再变换。为什么要加窗？因为整段截取等价于矩形窗，频谱泄漏会把特征频率的能量弄散到旁边频点。工程上我首选汉宁窗（Hann），它能较好兼顾主瓣宽度和旁瓣衰减。

import numpy as np def fft_spectrum(x: np.ndarray, fs: float, nperseg: int = 4096): x = x[:nperseg] win = np.hanning(nperseg) xw = x * win spec = np.fft.rfft(xw) / (nperseg / 2) freq = np.fft.rfftfreq(nperseg, d=1.0 / fs) mag = np.abs(spec) return freq, mag

参数说明：

• fs是振动通道采样率，比如2000Hz
• nperseg=4096对应约2秒窗长，频率分辨率就是fs/nperseg≈0.49Hz，能分辨转频附近的微小偏移
• 除以nperseg/2是把幅度谱归一化到实际物理幅值量纲，便于不同窗长之间互相比较

频谱图出来后，我会先用freq找特征频率。四极电机在转速1470rpm时转频约24.5Hz，转频的整数倍位置出现能量峰是正常的；如果转频附近出现边带，或者在某些谐波倍数上能量异常升高，说明可能存在轴承局部缺陷或转子不平衡。

4.2 能量比指标：比绝对幅值更稳的异常判定

电机在不同的负载、温度下，同频率的振动幅值差异很大，直接用绝对幅值设阈值，误报率偏高。我一般会看频带能量比，也就是特征频带能量占整个分析频段能量的比例，这个比例在不同工况间的波动远小于绝对幅值。

表格里第三行加的1e-6是防止除零。阈值建议不是一个固定值，要基于正常历史数据取分位数。比如取正常运行数据该指标95分位数作为报警线，比拍脑袋的常数靠谱得多。计算频带能量比时，直接用numpy切片取幅度谱区间再做求和，注意区间边界的频率索引要换算成物理频率。

4.3 瞬态冲击的滑窗检测：捕捉轴承点蚀的早期信号

轴承点蚀初期的频谱特征不明显，但时域里会出现等间隔的瞬态冲击。检测方法是滑窗求波峰因子（crest factor），也就是窗口内峰值除以RMS值。

def crest_factor_series(x: np.ndarray, win: int = 1000, step: int = 200): n = len(x) n_out = 1 + (n - win) // step idx = np.arange(win)[None, :] + step * np.arange(n_out)[:, None] windows = x[idx] peak = np.max(np.abs(windows), axis=1) rms = np.sqrt(np.mean(windows ** 2, axis=1)) return peak / (rms + 1e-8) cf = crest_factor_series(vibration_signal, win=1000, step=200) alert_idx = np.where(cf > 8)[0] # 波峰因子超过8视为冲击

win=1000在2kHz采样率下是0.5秒窗口，cf > 8意味着窗口内存在远高于平均振动能量的冲击峰。这个方法在轴承早期点蚀检测中比整体频谱阈值更早发现异常，代价是无法定量缺陷尺寸，但足够用于触发离线诊断。

5. 从离线数据到监测参数：用nan_wt05.zip沉淀电机状态规则的几个技巧

5.1 建立温度-电流线性基线

电机绕组温度与电流平方近似线性关系。取前述稳态段数据，对绕组温度关于电流平方做一次线性回归：

from numpy.polynomial import polynomial as P x = current_rms ** 2 y = temp_winding coeff = P.polyfit(x, y, deg=1) # [b, a] residual = y - P.polyval(x, coeff)

参数说明：

• polyfit(deg=1)返回常数项b和一次项a
• residual是实际温度与模型预测的差值
• 正常运行残差通常在±2°C内，如果随运行时间单调走高，说明同电流下温度在上升，散热能力衰减

5.2 残差滑动均值：把时变漂移转成可持续观测的指标

残差序列受环境温度、冷却液温度影响较大，单点判断容易误报。我一般把残差做30分钟滑动均值，再统计连续上升窗口数量。关键看持续而不是瞬时。

resid_ma = pd.Series(residual).rolling(1800, min_periods=600).mean() trend_up = resid_ma.diff().rolling(300).mean() > 0.05

代码假设对齐后数据是1Hz，rolling(1800)即30分钟窗口。resid_ma.diff()判断残差均值是否在上升，再做一次300点平滑，连续多个窗口为正才报警。

5.3 把整个分析流程固化成可配置的监测模板

最后，把整个清洗、对齐、特征、频谱、回归流程的参数抽出来放到YAML配置里。后续拿到新的电机数据包，只改通道映射和阈值，脚本不动。这也是承担多台电机数据看护工作时减少重复劳动的常见做法。

motor: id: nan_wt05 fs_vibration: 2000 fs_current: 1000 steady_min_duration: 30 crest_factor_threshold: 8 temp_residual_alarm: 3.0 channels: vibration: Vibr_X current: I_PhaseA temp_winding: Temp_Winding speed_set: Speed_Set

这张配置表可以直接对应前面各章用到的参数：清洗阶段的缺失阈值、对齐阶段的采样率、特征阶段的稳态段时长、频谱阶段的波峰因子阈值、回归阶段的温升残差报警线。参数化之后，离线分析脚本和在线监测脚本共用同一份定义，避免两侧阈值不一致导致同样的数据在不同环境里结论打架。

本文还有配套的精品资源，点击获取