尧图网站建设 
简介一套面向无人机早期故障检测的MATLAB完整代码包主要服务计算机、电子信息工程、数学等专业的学生适用于课程设计、期末大作业与毕业设计场景。代码采用参数化编程思路清晰、注释明细并配套可直接运行的飞行案例数据使用者可快速修改参数进行验证与二次开发。资源共26个文件以.mat数据文件为主20个配合5个.m脚本文件与1个说明文档整体大小18.51MB.mat用于存放飞行工况与故障数据.m脚本覆盖数据预处理、BiLSTM分类实验、LSTM自编码器构建等关键检测流程md文档提供环境配置与使用说明。已有84人学习下载适合需要开展时序数据分类、异常检测及深度模型对比实验的读者拿到即可复现实验并在此基础上扩展自己的算法思路。1. 无人机早期故障检测matlab代码.zip 到底能帮你解决什么问题无人机早期故障检测matlab代码.zip 这类压缩包在网上一搜一大把但很多人解压之后发现跑不通、改不动最后只能照抄一份毕业设计。真正能落地的代码应该回答三个问题故障在出现前多久能被发现靠什么特征发现报警阈值是怎么定出来的。这里说的“早期”不是等电机已经抖到飞不起来而是滚动轴承刚出现微米级剥落坑、螺旋桨出现细微裂纹时振动和电流信号里刚冒出来的那一丁点周期性冲击成分。对做毕设、准备技术答辩、给无人机巡检产线做预研的工程师来说这套流程的核心不在某个 .m 文件而在“特征提取→阈值生成→趋势预警”这条链路上。理解这条链路才算真正拿到这个压缩包里的东西。2. 早期故障检测难点拆解微弱信号、特征频率与 MATLAB 能做什么2.1 为什么无人机的早期故障特别难发现无人机和工业旋转机械最大的区别在于运行状态一直在变。工业电机可以稳定在额定转速下做长期监测而无人机从起飞到悬停再到机动电机转速可能有 30% 以上的浮动螺旋桨噪声、风噪、舵机换向噪声都是宽带成分。早期故障产生的冲击幅值往往只有正常振动幅值的十分之一甚至更低直接对原始振动信号做 FFT频谱上很难看出异常。另一个容易被忽略的点是采样率。飞控里陀螺仪和加速度计的采样率通常在 100Hz 到 500Hz那是给姿态控制用的根本覆盖不了轴承故障的特征频率。早期故障检测需要独立于飞控的数据采集链路振动通道采样率至少 12.8kHz电流通道至少 5kHz。这个前提不满足后面用什么算法都白搭。2.2 先把特征频率定下来电机电气故障与滚动轴承早期故障检测的第一步不是跑代码而是算出“要找什么频率”。对于无刷电机驱动的无人机主要看三类特征频率。电机电气故障特征频率为 2sfs 是转差率f 是供电频率。但在无人机场景里供电来自无刷电调而不是工频电网f 对应的是换相频率或 PWM 基频需要根据电机的极对数和电调控制方式折算。这个细节很多代码包没有处理好直接套 2sf 公式会出错。滚动轴承的特征频率计算公式相对固定。以外圈故障为例BPFO (n / 2) * fr * (1 - d/D * cos α)内圈故障BPFI (n / 2) * fr * (1 d/D * cos α)其中 n 是滚动体个数fr 是转频d 是滚动体直径D 是轴承节圆直径α 是接触角。内圈故障的特点是特征频率会被转频调制包络谱上可以看到 BPFI 及其边带。螺旋桨不平衡的特征频率就是 1×RPM但 RPM 随油门变化所以不能用固定频点去判断要用转速同步或者比例带宽追踪。2.3 MATLAB 在早期故障检测里的定位信号处理与自动阈值很多做深度学习的人会把早期故障检测理解为“训练一个分类网络”但在工程落地时深度学习在无人机场景里往往不实用早期故障样本太难采集正常数据倒是有一大堆。MATLAB 在这个场景里的定位是快速搭建一条可解释的信号处理链路从频域特征提取到统计阈值生成脚本改起来比 Python 还要顺手而且能直接和 Simulink 模型联合验证。在 Ubuntu 上搭过 PX4 无人机仿真的人会深有体会飞控仿真解决的是控制逻辑验证而故障检测需要的是真实传感器数据。MATLAB 里可以先离线把算法跑通再决定是否部署到机载处理器这个验证方式对毕设和产线预研都够用。2.4 关键对象与工具速查表监测对象敏感特征频率早期阶段表现MATLAB 工具箱电机滚动轴承BPFO / BPFI随转速漂移原始频谱出现边带包络谱出现清晰峰Signal Processing Toolbox电机电气换相2sf 谐波与边带电流谐波增大振动出现换相频率成分自建滤波与频谱分析螺旋桨叶片1×RPM 及 2 倍频基频幅值缓慢上升裂纹时 2 倍频非线性增长频域峰值追踪电池与供电电压纹波增大电流信号出现低频调制时域统计特征2.5 别把飞控控制和故障检测混在一起飞控里的串级 PID、LQR 控制器解决的是姿态稳定问题输入是 IMU 和磁力计数据采样率低、量程小本质上和故障检测不是一套体系。早期故障检测需要高频振动和电流数据处理的是机械和电气特征别指望从飞控日志里直接提取轴承故障特征。这个误区在不少开源代码里都存在下载数据包之前先确认数据来源。3. 从数据到特征矩阵MATLAB 读取、包络谱与特征提取实现3.1 数据目录约定与统一加载函数大多数人下载的 zip 包里文件命名和格式都不统一有 .mat、有 .csv列的排列顺序也不一致。我习惯先写一个统一入口把杂乱的原始文件整理成 MATLAB 结构体数组后续所有特征提取代码都只依赖这个结构体改动面积最小。function data load_uav_data(datasetDir, fs, channelNames) % 读取批量数据文件统一字段和标签 % 输入 % datasetDir 字符串数据文件夹路径 % fs 采样率单位 Hz建议 12800 或 51200 % channelNames 元胞数组例如 {vib_x,vib_y,current} % 输出 % data 结构体数组每一条包含 label、fs、signal files dir(fullfile(datasetDir, *.mat)); if isempty(files) files dir(fullfile(datasetDir, *.csv)); end for idx 1:length(files) [~, name, ~] fileparts(files(idx).name); raw load(fullfile(files(idx).folder, files(idx).name)); % 统一字段fault_type 为标签signal 为 nChannels x nSamples if isfield(raw, fault_type) data(idx).label raw.fault_type; elseif isfield(raw, label) data(idx).label raw.label; else data(idx).label normal; end data(idx).name name; data(idx).fs fs; data(idx).signal raw.signal(1:numel(channelNames), :); end end参数说明里最需要注意的是 channelNames 的顺序。有的数据包里 vibration 放第一行有的放第二行如果后面做特征提取时通道对应错了所有结论都会反转。建议加载后先打印 size(data(1).signal) 确认矩阵维度是“通道数 × 采样点数”如果方向反了用转置修正。3.2 用 FFT 与包络谱把特征频率“逼出来”早期故障的冲击成分信号很弱直接对原始信号做 FFT 往往看不到特征频率处的峰值。原因是冲击会激起结构的高频固有共振而共振频带里的信噪比比基频处高得多所以工程上一般不直接看原始频谱而是先做带通滤波再做包络解调最后对包络信号做 FFT这个流程叫高频共振解调。function [envSpectrum, freqAxis] envelope_spectrum(signal, fs, lowFreq, highFreq) % 高频共振解调提取包络谱 % 参数说明 % lowFreq / highFreq 带通滤波下边界和上边界单位 Hz % 常见做法是选在 1000 Hz 到 5000 Hz避开螺旋桨基频区域 signal signal - mean(signal); [b, a] butter(4, [lowFreq highFreq]/(fs/2), bandpass); band filtfilt(b, a, signal); % 零相位滤波避免相位偏移 envelope abs(hilbert(band)); % 希尔伯特变换求解析信号再取模 envSpectrum abs(fft(envelope)); freqAxis (0:length(envSpectrum)-1) * fs / length(envSpectrum); envSpectrum(1) 0; % 去掉直流分量 end代码里的 butter(4, ...) 是四阶巴特沃斯滤波器filtfilt 做零相位滤波保证冲击位置不发生偏移。hilbert 求解析信号后用 abs 得到包络最后对包络做 FFT 就得到包络谱。注意 lowFreq 的选择很关键如果带宽选到了 800 Hz 以下螺旋桨的正常转动噪声会大量进入如果选到 10kHz 以上可能把电调 PWM 开关噪声带进来。具体取值要先用 pwelch 看一下原始频谱的共振峰分布再定滤波边界。3.3 特征向量构建时域统计加频域峰值特征提取是把一段信号浓缩成一个固定长度的向量方便后续做阈值判断。下面是特征提取函数的核心部分。function feat extract_features(signal, fs, targetFreqs, bandWidth) % 从单通道信号中提取特征向量 % targetFreqs 已知特征频率向量例如 [BPFO, BPFI, 2*转频] % bandWidth 特征频率两侧搜索局部峰值的半宽建议 25~50 Hz feat []; feat(1) sqrt(mean(signal.^2)); % RMS 反映振动能量 feat(2) max(abs(signal)) / feat(1); % 峰值因子冲击性指标 [spec, fx] pwelch(signal, [], [], [], fs); % 功率谱密度 for k 1:length(targetFreqs) idx find(fx targetFreqs(k)-bandWidth fx targetFreqs(k)bandWidth); if isempty(idx) feat(2k) 0; else feat(2k) max(spec(idx)); % 特征频点峰值 end end % 边带能量比反映调制强度内圈故障很敏感 if length(targetFreqs) 0 f1 targetFreqs(1); sideLow find(fx f1-3*bandWidth fx f1-bandWidth); sideHigh find(fx f1bandWidth fx f13*bandWidth); feat(end1) sum(spec(sideLow)) sum(spec(sideHigh)); end end这段代码里 RMS 和峰值因子反映信号的能量和冲击强度特征频点峰值对应故障的特征成分。边带能量比是内圈故障的敏感指标因为内圈故障产生的冲击会被转频调制在包络谱上表现为 BPFI 附近出现等间隔边带。如果监测对象是螺旋桨建议把边带能量比替换成 2 倍频与基频的能量比裂纹初期这个比值呈非线性上升。4. 阈值不能拍脑袋自适应阈值与残差分析4.1 固定阈值的三个坑很多人在验证故障检测代码时发现误报率特别高第一个该背锅的就是固定阈值。无人机飞行环境温度变化大轴承温度升高后振动整体抬升固定阈值会在下午两点连续误报挂载负载变了电机转速脉动也会变特征值整体漂移再加上正常磨损会让特征值缓慢上升一周之后固定阈值就可能被正常数据击穿这就是运维里常说的“固定阈值失效”。4.2 滑动窗口自适应阈值算法自适应阈值的核心是用滑动窗口估计特征序列的均值和标准差再按均值加 k 倍标准差生成上下阈值。k 取 3 时对应正态分布的 99.7% 置信区间这是统计过程控制里最常用的经验值。function [thU, thL, mu] adaptive_threshold(featureSeq, windowLen, kFactor) % 滑动窗口自适应阈值 % featureSeq特征随时间变化的列向量单位是特征值本身 % windowLen窗口长度若每天一个特征点建议 24 或 48 % kFactor标准差倍数推荐 3误报敏感的场景可调到 3.5 n length(featureSeq); thU zeros(n, 1); thL zeros(n, 1); mu zeros(n, 1); for i windowLen:n seg featureSeq(i-windowLen1:i); mu(i) mean(seg); sigma std(seg); thU(i) mu(i) kFactor * sigma; thL(i) mu(i) - kFactor * sigma; end % 前 windowLen-1 个点历史不足用第一个有效窗口值填充 thU(1:windowLen-1) thU(windowLen); thL(1:windowLen-1) thL(windowLen); end窗口长度的选择要和数据采样间隔匹配。如果数据是每 10 分钟计算一个特征点窗口长度取 144 就是一天的滑动基线能适应昼夜温差带来的特征漂移。kFactor 调小的代价是漏报率降低但误报率升高实际调参时建议先固定 k3 跑一遍正常数据数一下一天内误报次数如果高于 2 次再逐步提高到 3.5。4.3 残差分析消除健康基线漂移自适应阈值只是解决了阈值随状态漂移的问题但真正触发报警的应当是特征值相对健康基线的“变化量”而不是绝对值。残差分析的做法是把当前特征值减去滑动窗口均值然后看残差是否持续为正且和零有显著距离。% 计算残差与连续报警判定 [thU, ~, mu] adaptive_threshold(featureSeq, windowLen, kFactor); residual featureSeq - mu; % 减去健康基线 limit kFactor * std(featureSeq(1:windowLen)); % 残差控制限 flags residual limit; % 单点超限标记 alarmIdx find(movmean(flags, 5) 1); % 连续 5 个点超限才报警这里的 limit 用的是初始健康段的波动范围反映的是“相对健康状态的偏差”而不是原始特征值的绝对值。连续 5 个点超限才报警的设计是为了过滤单点脉冲噪声飞控数据里偶发的射频干扰或者数传丢包不会触发误报。如果报警点集中在某几个时间点优先检查数据采集链路是否有间歇性接触不良而不是急着改算法。5. 从单点报警到趋势预警EWMA 控制图在线实现5.1 单点报警的局限上一章的自适应阈值虽然能适应基线漂移但本质还是“这一秒的特征值是否超过这一秒的阈值”。实际早期故障的特征是缓慢爬升单点超阈往往发生在故障已经比较明显之后。想要提前预警得把判定逻辑从“单点”变成“趋势”。指数加权移动平均EWMA是工业过程控制里最常用的趋势监控工具它给最近的观测值更高的权重能够在特征值还没有超过硬阈值时就捕捉到趋势变化。EWMA 的递推公式是 z_t λx_t (1-λ)z_{t-1}λ 为平滑系数越大表示越相信当前观测值越小表示越相信历史信息。5.2 EWMA 参数表λ、控制限因子与推荐取值参数含义推荐范围调参说明λ平滑系数0.1 ~ 0.3越小越平滑但响应越慢越大越灵敏但容易跟着噪声走L控制限因子2.7 ~ 3.0用于计算 UCL/LCLL 越大误报越少initLen健康基线长度50 ~ 200必须覆盖至少一个完整飞行起落周期采样间隔特征点计算周期5 ~ 15 分钟过短趋势不明显过长会延误报警这里要特别强调的是 initLen 的选择。EWMA 的初始值 z(1) 来自健康阶段的均值如果 initLen 里混入了已经发生早期的故障数据基线就会整体抬高后续报警灵敏度明显下降。正确做法是先画特征序列的时间曲线找到特征值平稳的下半段作为基线而不是直接取前 100 个点。5.3 MATLAB 实现 EWMA 控制图function [z, UCL, LCL, alarmIdx] ewma_control(featureSeq, lambda, L, initLen) % EWMA 控制图在线预警 % 输入 % featureSeq 特征值列向量例如 RMS 或特征频点峰值 % lambda 平滑系数推荐 0.15 % L 控制限因子推荐 2.8 % initLen 健康基线窗口长度 % 输出 % z 平滑后的特征序列 % UCL / LCL 上控制限 / 下控制限 % alarmIdx 报警位置索引 n length(featureSeq); z zeros(n, 1); UCL zeros(n, 1); LCL zeros(n, 1); init featureSeq(1:initLen); mu0 mean(init); sigma0 std(init); z(1) mu0; for t 2:n z(t) lambda * featureSeq(t) (1 - lambda) * z(t-1); % EWMA 方差随 t 收敛用标准公式计算控制限 sigmaZ sigma0 * sqrt(lambda / (2 - lambda) * (1 - (1-lambda)^(2*t))); UCL(t) mu0 L * sigmaZ; LCL(t) mu0 - L * sigmaZ; end alarmIdx find(z UCL | z LCL); end代码的核心逻辑在循环体里z(t) 是当前时刻的平滑值UCL 和 LCL 是围绕健康均值 mu0 的动态控制限。随着 t 增大(1-λ)^(2t) 趋近于零控制限会逐渐稳定到常值。如果报警idx出现在 t 很小的时候大概率是 initLen 选的健康基线有问题或者特征序列一开始就处于劣化状态。5.4 趋势预警前的数据质量检查EWMA 对数据质量的要求比自适应阈值更高。如果特征序列里有明显的跳跃点比如数传中断导致补零EWMA 会把跳变当作趋势信号拉响警报。跑控制图之前先做一次中值滤波把孤立野值剔除中值窗口取 3 或 5 即可太大反而把真实趋势抹平了。无人机巡检场景下一个架次 30 分钟特征值一分钟算一个的话一个架次只有 30 个特征点此时就不适合用 EWMA直接用单点加迟滞判定更实用。6. 用仿真信号校准整套检测脚本手头没有真实故障数据时可以自己生成一段仿真信号来验证代码链路是否通顺这个步骤能帮你避免把数据集里的偶然误差当作算法效果。fs 25600; t 0:1/fs:10-1/fs; fr 50; % 转频 50 Hz sigNormal 0.5*sin(2*pi*fr*t) 0.2*sin(2*pi*2*fr*t) 0.05*randn(size(t)); % 注入早期轴承外圈故障冲击幅值只有正常谐波的 6% impulseTrain zeros(size(t)); impulseIdx 1:round(fs/fr):length(t); for k 1:length(impulseIdx) if impulseIdx(k) length(t) - 1 break; end len min(30, length(t) - impulseIdx(k) 1); idxLen 0:len-1; impulseTrain(impulseIdx(k):impulseIdx(k)len-1) ... 0.03 * exp(-1000*idxLen/fs) .* sin(2*pi*3000*idxLen/fs); end sigFault sigNormal impulseTrain;这段代码里冲击幅值 0.03正常谐波幅值 0.5占比只有 6%放在时域图里肉眼几乎看不出来但包络谱里 BPFO 位置应当出现清晰峰值。把 sigFault 的前半段当作健康数据、后半段当作故障数据切成 1 秒一段跑一遍第 3 章和第 4 章的流程然后画特征值时间序列观察故障注入点之后特征值是否呈持续上升。调参顺序建议是先固定带通边界为 1000~5000 Hz然后调 windowLen 让阈值在健康段不误报最后调 λ 和 L让报警时刻比故障注入点提前至少 20 秒。如果 EWMA 报警太晚把 λ 从 0.15 调到 0.2报警提前量会明显改善如果开始误报优先确认是不是 initLen 里混入了注入故障后的样本。验证完成后在真实无人机数据集上跑最需要注意的还是转速同步问题——特征频率随转速漂移时直接固定频点搜索会丢峰这时把 targetFreqs 改成按转频比例计算的方式最稳妥。本文还有配套的精品资源点击获取 无人机早期故障检测MATLAB代码：特征提取与趋势预警实战
简介一套面向无人机早期故障检测的MATLAB完整代码包主要服务计算机、电子信息工程、数学等专业的学生适用于课程设计、期末大作业与毕业设计场景。代码采用参数化编程思路清晰、注释明细并配套可直接运行的飞行案例数据使用者可快速修改参数进行验证与二次开发。资源共26个文件以.mat数据文件为主20个配合5个.m脚本文件与1个说明文档整体大小18.51MB.mat用于存放飞行工况与故障数据.m脚本覆盖数据预处理、BiLSTM分类实验、LSTM自编码器构建等关键检测流程md文档提供环境配置与使用说明。已有84人学习下载适合需要开展时序数据分类、异常检测及深度模型对比实验的读者拿到即可复现实验并在此基础上扩展自己的算法思路。1. 无人机早期故障检测matlab代码.zip 到底能帮你解决什么问题无人机早期故障检测matlab代码.zip 这类压缩包在网上一搜一大把但很多人解压之后发现跑不通、改不动最后只能照抄一份毕业设计。真正能落地的代码应该回答三个问题故障在出现前多久能被发现靠什么特征发现报警阈值是怎么定出来的。这里说的“早期”不是等电机已经抖到飞不起来而是滚动轴承刚出现微米级剥落坑、螺旋桨出现细微裂纹时振动和电流信号里刚冒出来的那一丁点周期性冲击成分。对做毕设、准备技术答辩、给无人机巡检产线做预研的工程师来说这套流程的核心不在某个 .m 文件而在“特征提取→阈值生成→趋势预警”这条链路上。理解这条链路才算真正拿到这个压缩包里的东西。2. 早期故障检测难点拆解微弱信号、特征频率与 MATLAB 能做什么2.1 为什么无人机的早期故障特别难发现无人机和工业旋转机械最大的区别在于运行状态一直在变。工业电机可以稳定在额定转速下做长期监测而无人机从起飞到悬停再到机动电机转速可能有 30% 以上的浮动螺旋桨噪声、风噪、舵机换向噪声都是宽带成分。早期故障产生的冲击幅值往往只有正常振动幅值的十分之一甚至更低直接对原始振动信号做 FFT频谱上很难看出异常。另一个容易被忽略的点是采样率。飞控里陀螺仪和加速度计的采样率通常在 100Hz 到 500Hz那是给姿态控制用的根本覆盖不了轴承故障的特征频率。早期故障检测需要独立于飞控的数据采集链路振动通道采样率至少 12.8kHz电流通道至少 5kHz。这个前提不满足后面用什么算法都白搭。2.2 先把特征频率定下来电机电气故障与滚动轴承早期故障检测的第一步不是跑代码而是算出“要找什么频率”。对于无刷电机驱动的无人机主要看三类特征频率。电机电气故障特征频率为 2sfs 是转差率f 是供电频率。但在无人机场景里供电来自无刷电调而不是工频电网f 对应的是换相频率或 PWM 基频需要根据电机的极对数和电调控制方式折算。这个细节很多代码包没有处理好直接套 2sf 公式会出错。滚动轴承的特征频率计算公式相对固定。以外圈故障为例BPFO (n / 2) * fr * (1 - d/D * cos α)内圈故障BPFI (n / 2) * fr * (1 d/D * cos α)其中 n 是滚动体个数fr 是转频d 是滚动体直径D 是轴承节圆直径α 是接触角。内圈故障的特点是特征频率会被转频调制包络谱上可以看到 BPFI 及其边带。螺旋桨不平衡的特征频率就是 1×RPM但 RPM 随油门变化所以不能用固定频点去判断要用转速同步或者比例带宽追踪。2.3 MATLAB 在早期故障检测里的定位信号处理与自动阈值很多做深度学习的人会把早期故障检测理解为“训练一个分类网络”但在工程落地时深度学习在无人机场景里往往不实用早期故障样本太难采集正常数据倒是有一大堆。MATLAB 在这个场景里的定位是快速搭建一条可解释的信号处理链路从频域特征提取到统计阈值生成脚本改起来比 Python 还要顺手而且能直接和 Simulink 模型联合验证。在 Ubuntu 上搭过 PX4 无人机仿真的人会深有体会飞控仿真解决的是控制逻辑验证而故障检测需要的是真实传感器数据。MATLAB 里可以先离线把算法跑通再决定是否部署到机载处理器这个验证方式对毕设和产线预研都够用。2.4 关键对象与工具速查表监测对象敏感特征频率早期阶段表现MATLAB 工具箱电机滚动轴承BPFO / BPFI随转速漂移原始频谱出现边带包络谱出现清晰峰Signal Processing Toolbox电机电气换相2sf 谐波与边带电流谐波增大振动出现换相频率成分自建滤波与频谱分析螺旋桨叶片1×RPM 及 2 倍频基频幅值缓慢上升裂纹时 2 倍频非线性增长频域峰值追踪电池与供电电压纹波增大电流信号出现低频调制时域统计特征2.5 别把飞控控制和故障检测混在一起飞控里的串级 PID、LQR 控制器解决的是姿态稳定问题输入是 IMU 和磁力计数据采样率低、量程小本质上和故障检测不是一套体系。早期故障检测需要高频振动和电流数据处理的是机械和电气特征别指望从飞控日志里直接提取轴承故障特征。这个误区在不少开源代码里都存在下载数据包之前先确认数据来源。3. 从数据到特征矩阵MATLAB 读取、包络谱与特征提取实现3.1 数据目录约定与统一加载函数大多数人下载的 zip 包里文件命名和格式都不统一有 .mat、有 .csv列的排列顺序也不一致。我习惯先写一个统一入口把杂乱的原始文件整理成 MATLAB 结构体数组后续所有特征提取代码都只依赖这个结构体改动面积最小。function data load_uav_data(datasetDir, fs, channelNames) % 读取批量数据文件统一字段和标签 % 输入 % datasetDir 字符串数据文件夹路径 % fs 采样率单位 Hz建议 12800 或 51200 % channelNames 元胞数组例如 {vib_x,vib_y,current} % 输出 % data 结构体数组每一条包含 label、fs、signal files dir(fullfile(datasetDir, *.mat)); if isempty(files) files dir(fullfile(datasetDir, *.csv)); end for idx 1:length(files) [~, name, ~] fileparts(files(idx).name); raw load(fullfile(files(idx).folder, files(idx).name)); % 统一字段fault_type 为标签signal 为 nChannels x nSamples if isfield(raw, fault_type) data(idx).label raw.fault_type; elseif isfield(raw, label) data(idx).label raw.label; else data(idx).label normal; end data(idx).name name; data(idx).fs fs; data(idx).signal raw.signal(1:numel(channelNames), :); end end参数说明里最需要注意的是 channelNames 的顺序。有的数据包里 vibration 放第一行有的放第二行如果后面做特征提取时通道对应错了所有结论都会反转。建议加载后先打印 size(data(1).signal) 确认矩阵维度是“通道数 × 采样点数”如果方向反了用转置修正。3.2 用 FFT 与包络谱把特征频率“逼出来”早期故障的冲击成分信号很弱直接对原始信号做 FFT 往往看不到特征频率处的峰值。原因是冲击会激起结构的高频固有共振而共振频带里的信噪比比基频处高得多所以工程上一般不直接看原始频谱而是先做带通滤波再做包络解调最后对包络信号做 FFT这个流程叫高频共振解调。function [envSpectrum, freqAxis] envelope_spectrum(signal, fs, lowFreq, highFreq) % 高频共振解调提取包络谱 % 参数说明 % lowFreq / highFreq 带通滤波下边界和上边界单位 Hz % 常见做法是选在 1000 Hz 到 5000 Hz避开螺旋桨基频区域 signal signal - mean(signal); [b, a] butter(4, [lowFreq highFreq]/(fs/2), bandpass); band filtfilt(b, a, signal); % 零相位滤波避免相位偏移 envelope abs(hilbert(band)); % 希尔伯特变换求解析信号再取模 envSpectrum abs(fft(envelope)); freqAxis (0:length(envSpectrum)-1) * fs / length(envSpectrum); envSpectrum(1) 0; % 去掉直流分量 end代码里的 butter(4, ...) 是四阶巴特沃斯滤波器filtfilt 做零相位滤波保证冲击位置不发生偏移。hilbert 求解析信号后用 abs 得到包络最后对包络做 FFT 就得到包络谱。注意 lowFreq 的选择很关键如果带宽选到了 800 Hz 以下螺旋桨的正常转动噪声会大量进入如果选到 10kHz 以上可能把电调 PWM 开关噪声带进来。具体取值要先用 pwelch 看一下原始频谱的共振峰分布再定滤波边界。3.3 特征向量构建时域统计加频域峰值特征提取是把一段信号浓缩成一个固定长度的向量方便后续做阈值判断。下面是特征提取函数的核心部分。function feat extract_features(signal, fs, targetFreqs, bandWidth) % 从单通道信号中提取特征向量 % targetFreqs 已知特征频率向量例如 [BPFO, BPFI, 2*转频] % bandWidth 特征频率两侧搜索局部峰值的半宽建议 25~50 Hz feat []; feat(1) sqrt(mean(signal.^2)); % RMS 反映振动能量 feat(2) max(abs(signal)) / feat(1); % 峰值因子冲击性指标 [spec, fx] pwelch(signal, [], [], [], fs); % 功率谱密度 for k 1:length(targetFreqs) idx find(fx targetFreqs(k)-bandWidth fx targetFreqs(k)bandWidth); if isempty(idx) feat(2k) 0; else feat(2k) max(spec(idx)); % 特征频点峰值 end end % 边带能量比反映调制强度内圈故障很敏感 if length(targetFreqs) 0 f1 targetFreqs(1); sideLow find(fx f1-3*bandWidth fx f1-bandWidth); sideHigh find(fx f1bandWidth fx f13*bandWidth); feat(end1) sum(spec(sideLow)) sum(spec(sideHigh)); end end这段代码里 RMS 和峰值因子反映信号的能量和冲击强度特征频点峰值对应故障的特征成分。边带能量比是内圈故障的敏感指标因为内圈故障产生的冲击会被转频调制在包络谱上表现为 BPFI 附近出现等间隔边带。如果监测对象是螺旋桨建议把边带能量比替换成 2 倍频与基频的能量比裂纹初期这个比值呈非线性上升。4. 阈值不能拍脑袋自适应阈值与残差分析4.1 固定阈值的三个坑很多人在验证故障检测代码时发现误报率特别高第一个该背锅的就是固定阈值。无人机飞行环境温度变化大轴承温度升高后振动整体抬升固定阈值会在下午两点连续误报挂载负载变了电机转速脉动也会变特征值整体漂移再加上正常磨损会让特征值缓慢上升一周之后固定阈值就可能被正常数据击穿这就是运维里常说的“固定阈值失效”。4.2 滑动窗口自适应阈值算法自适应阈值的核心是用滑动窗口估计特征序列的均值和标准差再按均值加 k 倍标准差生成上下阈值。k 取 3 时对应正态分布的 99.7% 置信区间这是统计过程控制里最常用的经验值。function [thU, thL, mu] adaptive_threshold(featureSeq, windowLen, kFactor) % 滑动窗口自适应阈值 % featureSeq特征随时间变化的列向量单位是特征值本身 % windowLen窗口长度若每天一个特征点建议 24 或 48 % kFactor标准差倍数推荐 3误报敏感的场景可调到 3.5 n length(featureSeq); thU zeros(n, 1); thL zeros(n, 1); mu zeros(n, 1); for i windowLen:n seg featureSeq(i-windowLen1:i); mu(i) mean(seg); sigma std(seg); thU(i) mu(i) kFactor * sigma; thL(i) mu(i) - kFactor * sigma; end % 前 windowLen-1 个点历史不足用第一个有效窗口值填充 thU(1:windowLen-1) thU(windowLen); thL(1:windowLen-1) thL(windowLen); end窗口长度的选择要和数据采样间隔匹配。如果数据是每 10 分钟计算一个特征点窗口长度取 144 就是一天的滑动基线能适应昼夜温差带来的特征漂移。kFactor 调小的代价是漏报率降低但误报率升高实际调参时建议先固定 k3 跑一遍正常数据数一下一天内误报次数如果高于 2 次再逐步提高到 3.5。4.3 残差分析消除健康基线漂移自适应阈值只是解决了阈值随状态漂移的问题但真正触发报警的应当是特征值相对健康基线的“变化量”而不是绝对值。残差分析的做法是把当前特征值减去滑动窗口均值然后看残差是否持续为正且和零有显著距离。% 计算残差与连续报警判定 [thU, ~, mu] adaptive_threshold(featureSeq, windowLen, kFactor); residual featureSeq - mu; % 减去健康基线 limit kFactor * std(featureSeq(1:windowLen)); % 残差控制限 flags residual limit; % 单点超限标记 alarmIdx find(movmean(flags, 5) 1); % 连续 5 个点超限才报警这里的 limit 用的是初始健康段的波动范围反映的是“相对健康状态的偏差”而不是原始特征值的绝对值。连续 5 个点超限才报警的设计是为了过滤单点脉冲噪声飞控数据里偶发的射频干扰或者数传丢包不会触发误报。如果报警点集中在某几个时间点优先检查数据采集链路是否有间歇性接触不良而不是急着改算法。5. 从单点报警到趋势预警EWMA 控制图在线实现5.1 单点报警的局限上一章的自适应阈值虽然能适应基线漂移但本质还是“这一秒的特征值是否超过这一秒的阈值”。实际早期故障的特征是缓慢爬升单点超阈往往发生在故障已经比较明显之后。想要提前预警得把判定逻辑从“单点”变成“趋势”。指数加权移动平均EWMA是工业过程控制里最常用的趋势监控工具它给最近的观测值更高的权重能够在特征值还没有超过硬阈值时就捕捉到趋势变化。EWMA 的递推公式是 z_t λx_t (1-λ)z_{t-1}λ 为平滑系数越大表示越相信当前观测值越小表示越相信历史信息。5.2 EWMA 参数表λ、控制限因子与推荐取值参数含义推荐范围调参说明λ平滑系数0.1 ~ 0.3越小越平滑但响应越慢越大越灵敏但容易跟着噪声走L控制限因子2.7 ~ 3.0用于计算 UCL/LCLL 越大误报越少initLen健康基线长度50 ~ 200必须覆盖至少一个完整飞行起落周期采样间隔特征点计算周期5 ~ 15 分钟过短趋势不明显过长会延误报警这里要特别强调的是 initLen 的选择。EWMA 的初始值 z(1) 来自健康阶段的均值如果 initLen 里混入了已经发生早期的故障数据基线就会整体抬高后续报警灵敏度明显下降。正确做法是先画特征序列的时间曲线找到特征值平稳的下半段作为基线而不是直接取前 100 个点。5.3 MATLAB 实现 EWMA 控制图function [z, UCL, LCL, alarmIdx] ewma_control(featureSeq, lambda, L, initLen) % EWMA 控制图在线预警 % 输入 % featureSeq 特征值列向量例如 RMS 或特征频点峰值 % lambda 平滑系数推荐 0.15 % L 控制限因子推荐 2.8 % initLen 健康基线窗口长度 % 输出 % z 平滑后的特征序列 % UCL / LCL 上控制限 / 下控制限 % alarmIdx 报警位置索引 n length(featureSeq); z zeros(n, 1); UCL zeros(n, 1); LCL zeros(n, 1); init featureSeq(1:initLen); mu0 mean(init); sigma0 std(init); z(1) mu0; for t 2:n z(t) lambda * featureSeq(t) (1 - lambda) * z(t-1); % EWMA 方差随 t 收敛用标准公式计算控制限 sigmaZ sigma0 * sqrt(lambda / (2 - lambda) * (1 - (1-lambda)^(2*t))); UCL(t) mu0 L * sigmaZ; LCL(t) mu0 - L * sigmaZ; end alarmIdx find(z UCL | z LCL); end代码的核心逻辑在循环体里z(t) 是当前时刻的平滑值UCL 和 LCL 是围绕健康均值 mu0 的动态控制限。随着 t 增大(1-λ)^(2t) 趋近于零控制限会逐渐稳定到常值。如果报警idx出现在 t 很小的时候大概率是 initLen 选的健康基线有问题或者特征序列一开始就处于劣化状态。5.4 趋势预警前的数据质量检查EWMA 对数据质量的要求比自适应阈值更高。如果特征序列里有明显的跳跃点比如数传中断导致补零EWMA 会把跳变当作趋势信号拉响警报。跑控制图之前先做一次中值滤波把孤立野值剔除中值窗口取 3 或 5 即可太大反而把真实趋势抹平了。无人机巡检场景下一个架次 30 分钟特征值一分钟算一个的话一个架次只有 30 个特征点此时就不适合用 EWMA直接用单点加迟滞判定更实用。6. 用仿真信号校准整套检测脚本手头没有真实故障数据时可以自己生成一段仿真信号来验证代码链路是否通顺这个步骤能帮你避免把数据集里的偶然误差当作算法效果。fs 25600; t 0:1/fs:10-1/fs; fr 50; % 转频 50 Hz sigNormal 0.5*sin(2*pi*fr*t) 0.2*sin(2*pi*2*fr*t) 0.05*randn(size(t)); % 注入早期轴承外圈故障冲击幅值只有正常谐波的 6% impulseTrain zeros(size(t)); impulseIdx 1:round(fs/fr):length(t); for k 1:length(impulseIdx) if impulseIdx(k) length(t) - 1 break; end len min(30, length(t) - impulseIdx(k) 1); idxLen 0:len-1; impulseTrain(impulseIdx(k):impulseIdx(k)len-1) ... 0.03 * exp(-1000*idxLen/fs) .* sin(2*pi*3000*idxLen/fs); end sigFault sigNormal impulseTrain;这段代码里冲击幅值 0.03正常谐波幅值 0.5占比只有 6%放在时域图里肉眼几乎看不出来但包络谱里 BPFO 位置应当出现清晰峰值。把 sigFault 的前半段当作健康数据、后半段当作故障数据切成 1 秒一段跑一遍第 3 章和第 4 章的流程然后画特征值时间序列观察故障注入点之后特征值是否呈持续上升。调参顺序建议是先固定带通边界为 1000~5000 Hz然后调 windowLen 让阈值在健康段不误报最后调 λ 和 L让报警时刻比故障注入点提前至少 20 秒。如果 EWMA 报警太晚把 λ 从 0.15 调到 0.2报警提前量会明显改善如果开始误报优先确认是不是 initLen 里混入了注入故障后的样本。验证完成后在真实无人机数据集上跑最需要注意的还是转速同步问题——特征频率随转速漂移时直接固定频点搜索会丢峰这时把 targetFreqs 改成按转频比例计算的方式最稳妥。本文还有配套的精品资源点击获取