简介：一套面向无人机早期故障检测的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 先把特征频率定下来：电机电气故障与滚动轴承
早期故障检测的第一步不是跑代码，而是算出“要找什么频率”。对于无刷电机驱动的无人机，主要看三类特征频率。
电机电气故障特征频率为 2sf，s 是转差率，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(2+k) = 0; else feat(2+k) = max(spec(idx)); % 特征频点峰值 end end % 边带能量比：反映调制强度，内圈故障很敏感 if length(targetFreqs) > 0 f1 = targetFreqs(1); sideLow = find(fx > f1-3*bandWidth & fx < f1-bandWidth); sideHigh = find(fx > f1+bandWidth & fx < f1+3*bandWidth); feat(end+1) = 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-windowLen+1: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 调小的代价是漏报率降低但误报率升高，实际调参时建议先固定 k=3 跑一遍正常数据，数一下一天内误报次数，如果高于 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/LCL，L 越大误报越少 |
| 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 改成按转频比例计算的方式最稳妥。
本文还有配套的精品资源，点击获取
