简介面向电机数据分析与风电系统研究的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按文本存储压缩比通常在520倍之间如果某个文件压缩比不到2倍很可能内部已经是二进制或内嵌压缩格式后面读取方式要区分对待。文件清单还能暴露非数据文件比如.ini、.log、.xlsx这些偶尔携带台架号、电机型号、环境温度等元信息先记下来后面做数据校验时能对得上。命名规律在这一步基本能判断出来。ch1_20231101_100000.csv这种带通道名和起始时间戳的说明记录程序按通道分文件Rec_001.csv这种顺序编号的说明按下电或实验阶段切段。前者时间拼接靠文件名解析后者往往需要从数据内部的时间列再对齐。2.2 字段体检与三分类映射波形量、统计量、工况量电机数据字段通常能按物理含义分成三类字段类别典型字段示例采样率特征分析用途波形量I_PhaseA、U_BC、Vibr_X1kHz50kHz频谱分析、瞬态冲击统计量Temp_Winding、Temp_Bearing1Hz10Hz热趋势、温升模型工况量Speed_Set、Torque_Fbk、Mode事件触发/低速工况分段、基线标定分类方法不用人工逐列看用我常用的体检函数跑一遍import pandas as pd def inspect_motor_csv(path: str, nrows: int 200) - None: df pd.read_csv(path, nrowsnrows, encoding_errorsreplace, skipinitialspaceTrue) for col in df.columns: uniq df[col].nunique(dropnaTrue) sample df[col].dropna().iloc[0] if uniq else None print(f{col:24s} dtype{str(df[col].dtype):12s} funique{uniq:8d} first{str(sample)[:20]}) inspect_motor_csv(ch1_20231101_100000.csv)参数说明nrows200只抽样前200行探查场景下足够判断字段类型skipinitialspaceTrue专门处理数采软件导出CSV时逗号后带空格的坏习惯不加这个参数数值列会被解析成字符串并产生隐性缺失encoding_errorsreplace防止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_rmslambda 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(methodffill)等于把停机区硬写成上一段运行值不仅制造假数据还会让后续均值、方差特征全部失真。我处理缺失值的固定套路是先按缺失段长度分级。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_areainside, limitshort) 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], default0) }) return out代码逻辑说明用isna().cumsum()把连续缺失区间分组这是处理缺失段的通用trick每变化一次就产生新的段号interpolate(limit_areainside, limitshort)只允许插值填补缺失段内部且最多补5个点超过就放弃缺失段长度超过50个采样点直接标记state2对应停机或断连状态返回的state列作为后续模型的辅助特征比单纯把缺失值去掉多保留了一层因果信息这样处理之后停机段不会参与稳态工况统计但能用于标记数据覆盖范围。3.2 滑动窗口构造RMS、峰值与温度变化率特征清洗完原始量后下一步是构造时序特征。电机数据分析最常用的三类滑窗特征和参数配置如下特征名滑动窗口步长计算公式作用vibr_rms_1s2000点100点sqrt(mean(x^2))振动能量水平current_peak_1s2000点100点max(abs(x))瞬时冲击电流temp_drate_60s60点1点diff(temp)/60绕组温升速率滑窗特征不建议全部用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, axis1)) vibr_rms sliding_rms(vibration_signal, 2000, 100)参数说明window2000对应2kHz采样率下的1秒物理窗口能量意义上的秒级RMSstep100是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, d1.0 / fs) mag np.abs(spec) return freq, mag参数说明fs是振动通道采样率比如2000Hznperseg4096对应约2秒窗长频率分辨率就是fs/nperseg≈0.49Hz能分辨转频附近的微小偏移除以nperseg/2是把幅度谱归一化到实际物理幅值量纲便于不同窗长之间互相比较频谱图出来后我会先用freq找特征频率。四极电机在转速1470rpm时转频约24.5Hz转频的整数倍位置出现能量峰是正常的如果转频附近出现边带或者在某些谐波倍数上能量异常升高说明可能存在轴承局部缺陷或转子不平衡。4.2 能量比指标比绝对幅值更稳的异常判定电机在不同的负载、温度下同频率的振动幅值差异很大直接用绝对幅值设阈值误报率偏高。我一般会看频带能量比也就是特征频带能量占整个分析频段能量的比例这个比例在不同工况间的波动远小于绝对幅值。指标计算方式正常范围参考异常倾向转频能量比sum(mag[23.5:25.5]) / sum(mag[5:1000])8%20%转子不平衡时升高轴承外圈特征频率能量sum(mag[fs*0.08:fs*0.12]) / total2%6%局部缺陷时出现峰值2倍转频比值mag[2*f0] / (mag[f0]1e-6)0.5轴系不对中时1表格里第三行加的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), axis1) rms np.sqrt(np.mean(windows ** 2, axis1)) return peak / (rms 1e-8) cf crest_factor_series(vibration_signal, win1000, step200) alert_idx np.where(cf 8)[0] # 波峰因子超过8视为冲击win1000在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, deg1) # [b, a] residual y - P.polyval(x, coeff)参数说明polyfit(deg1)返回常数项b和一次项aresidual是实际温度与模型预测的差值正常运行残差通常在±2°C内如果随运行时间单调走高说明同电流下温度在上升散热能力衰减5.2 残差滑动均值把时变漂移转成可持续观测的指标残差序列受环境温度、冷却液温度影响较大单点判断容易误报。我一般把残差做30分钟滑动均值再统计连续上升窗口数量。关键看持续而不是瞬时。resid_ma pd.Series(residual).rolling(1800, min_periods600).mean() trend_up resid_ma.diff().rolling(300).mean() 0.05代码假设对齐后数据是1Hzrolling(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这张配置表可以直接对应前面各章用到的参数清洗阶段的缺失阈值、对齐阶段的采样率、特征阶段的稳态段时长、频谱阶段的波峰因子阈值、回归阶段的温升残差报警线。参数化之后离线分析脚本和在线监测脚本共用同一份定义避免两侧阈值不一致导致同样的数据在不同环境里结论打架。本文还有配套的精品资源点击获取