★ FEATURED ARTICLE 做故障诊断的朋友对“高维特征 分类器”这条路子应该不陌生。振动信号经过时域、频域、时频域特征提取之后少则几十维多则几百维直接扔给分类器往往不是最优选择。这时候就需要降维或特征变换。NDGE归一化判别图嵌入Normalized Discriminative Graph Embedding就是干这件事的把高维特征映射到一个低维判别子空间输出一个投影矩阵同时给出每个样本分别属于不同故障模式的后验概率并顺便把“不同维度选多少维准确率最高”这件事也一并解决。这篇文章我会从原理讲到Matlab代码适合正在做故障诊断、设备健康管理或者模式识别方向的同学参考尤其是那种“样本不多、特征很高、类别又有点重叠”的场景。这套方案在我看来其实是一个标准的“谱嵌入 概率分类器”组合NDGE负责找好的低维表示后面的概率分类器负责给出可解释的置信度。下面直接进入正题。1. 故障诊断的第一步为什么要用归一化判别图嵌入1.1 从“高维特征 分类器”的痛点说起很多刚接触故障诊断的人会想特征越多越好吧时域均值、峰值、均方根、峭度频域重心、故障特征频率幅值再加小波包能量、熵特征一凑就是几十上百维。但问题也随之而来高维特征带来的不一定是信息更多是冗余和噪声。特别是故障诊断这个场景故障样本往往很难采集。某个设备正常状态的数据很多但真正出现早期故障、复合故障的时候现场可能只攒下几十组样本。高维特征 小样本结果就是分类器非常容易过拟合训练集上准确率百分之九十九测试集一上就掉链子。所以降维不是“可选项”而是“必经环节”。传统的PCA只管最大方差不看类别标签降完之后可能把最容易区分故障的判别方向当成噪声丢掉LDA虽然用了标签信息但在特征维度高于样本数、数据不是高斯分布的时候类内散布矩阵会奇异数值上也不稳定。这就是需要图嵌入方法介入的地方。NDGE本质上是一种有监督的流形学习方法它既保留样本之间的局部几何结构又把“同类拉近、异类拉远”的判别信息揉进目标函数最后通过广义特征值分解得到投影矩阵。1.2 图嵌入到底是什么把样本关系画成一张图理解图嵌入我习惯用一个比较直观的说法把每个样本看成图上的一个节点节点和节点之间的边表示它们的相似程度。如果两个样本在原始高维空间里挨得近它们之间就有一条权值较大的边挨得远边权就小或者干脆不连。在故障诊断里我们进一步把边分成两类类内边两个样本属于同一种故障模式类间边两个样本属于不同的故障模式降维的目标也变得很明确投影之后有类内边的样本尽量靠在一起有类间边的样本尽量被拉开。这个思想比单纯欧氏距离要稳。因为单个距离只能代表“远近”而图结构描述的是样本在局部邻域里的连接关系它对离群点和尺度变化更鲁棒。你可以把图理解成一张人际关系网一个团队里的人内部联系紧密不同团队之间是竞争关系。映射到一个二维平面上之后我们希望原来的“内部联系”还在而“竞争关系”能被最大化体现。1.3 NDGE和普通判别图嵌入的差别普通判别图嵌入DGE框架已经能做到“类内紧凑、类间分离”那为什么还要加“归一化”三个字我实际跑下来普通DGE在理想数据上效果不差但一落到真实故障数据上就很容易出问题第一特征量纲差异太大。振动幅值可能到几十峭度可能是个位数熵又只有0到1之间。如果不先做特征标准化距离矩阵会被量纲大的特征主导近邻关系基本由幅值决定其他特征全是陪跑。第二图节点度不均衡。有的样本周围邻居很多有的样本周围邻居很少。这样在构造拉普拉斯矩阵时度大的节点会对目标函数产生更大影响最后投影方向被少数密集区域带跑。第三特征向量的尺度不固定。广义特征分解得到的特征向量长度不是唯一确定的如果不做单位化投影后的坐标会被某个特征向量放大或缩小导致分类器训练时数值不稳定。NDGE干的事情就是针对这三个问题加了三次规范化处理输入特征标准化、拉普拉斯矩阵对称归一化、投影向量单位化。这也是我最推荐直接用NDGE而不是自己搓普通DGE的原因。2. 归一化判别图嵌入的核心原理与公式化2.1 图构造类内图、类间图和热核权重先约定符号训练样本写成矩阵 (X)每一行是一个样本每一列是一个特征标签向量是 (Y)。假设有 (n) 个样本特征维度是 (m)故障模式类别数用 (C) 表示。第一步是计算样本之间的欧氏距离矩阵。对每一个样本 (x_i)找出它的 (k) 个邻居。邻居按标签分成两类如果邻居和 (x_i) 是同一个故障模式就把它放进类内图边权记为 (S_w(i,j))如果邻居和 (x_i) 是不同故障模式就把它放进类间图边权记为 (S_b(i,j))边权通常用热核函数计算[ S(i,j) \exp\left(-\frac{|x_i - x_j|^2}{2\sigma^2}\right) ]其中 (\sigma) 是热核带宽。这个参数很关键(\sigma) 太小距离稍微远一点的样本边权几乎变成0图就碎了(\sigma) 太大所有近邻的边权都差不多等于只保留了“是否为邻居”的信息丢掉了距离的连续变化。得到两个相似度矩阵之后还要做一步对称化处理。工程实现里我一般直接同步设置 (S(i,j)) 和 (S(j,i))这样就能自然保证对称性。然后构造类内图对应的度矩阵 (D_w \text{diag}(\sum_j S_w(i,j)))拉普拉斯矩阵 (L_w D_w - S_w)同理得到类间图的 (D_b) 和 (L_b)。2.2 目标函数和广义特征值问题投影矩阵怎么来图嵌入的核心目标函数长这样[ \max_{W} ;; J(W) \frac{\operatorname{tr}(W^T X^T L_b X W)}{\operatorname{tr}(W^T X^T L_w X W)} ]分子里的 (L_b) 是类间拉普拉斯它衡量的是“原本不同类的样本投影后距离有多大”分母里的 (L_w) 是类内拉普拉斯它衡量的是“原本同类的样本投影后有没有被拆散”。所以这个目标函数的意思是在同一个比值下让异类样本投影后尽可能远同时让同类样本投影后尽可能近。求解方法就是对下面这个广义特征方程做特征分解[ X^T L_b X v \lambda X^T L_w X v ]定义两个 (m \times m) 矩阵[ S_b X^T L_b X ] [ S_w X^T L_w X ]然后解 (S_b v \lambda S_w v)。取最大的 (d) 个特征值对应的特征向量排成一个矩阵[ W [v_1, v_2, \ldots, v_d] ]这个 (W) 就是标题里说的“NDGE最终投影矩阵”。它是一个 (m \times d) 的矩阵作用是把原始高维特征投影到 (d) 维判别子空间[ Z X W ]用投影后的 (Z) 去训练分类器比直接用原始高维特征稳定得多。这里有一个工程细节我要多说一句(S_w) 经常是奇异的尤其是特征维度 (m) 接近样本数 (n) 的时候。求解之前一定要加一个小的正则项我通常这样处理[ S_w S_w \lambda_{reg} \cdot \frac{\operatorname{trace}(S_w)}{m} \cdot I ](\lambda_{reg}) 取 (10^{-6}) 到 (10^{-3}) 都可以具体看特征维度。正则项太小广义特征分解会得到一堆离谱的数值正则项太大又会把真正有用的判别方向压掉。2.3 归一化到底归在哪三步和背后的原因我前面反复强调NDGE比普通DGE稳现在把三次归一化拆开讲。第一步特征标准化。对所有特征先做零均值单位方差的Z-score处理。这一步不是锦上添花而是必须。故障诊断里的特征来源五花八门有的量纲是电压有的是温度有的是无量纲指标。如果直接算欧氏距离温度可能贡献几十的偏差而无量纲指标只能贡献零点几。标准化之后每个特征在距离计算里拥有平等发言权。第二步拉普拉斯矩阵对称归一化。普通DGE用的是原始 (L_w) 和 (L_b)但真实数据里每个样本的邻居数量不一样有些节点度很高有些节点度很低。高节点度样本会在目标函数里被放大低节点度样本几乎被忽略。解决办法是把拉普拉斯矩阵变成对称归一化形式[ L_{sym} D^{-1/2} L D^{-1/2} ]这个操作会把每个节点的“影响力”拉回同一个量级。从谱聚类的角度看对称归一化拉普拉斯比随机游走归一化更容易保持几何解释性特征向量也更稳定。第三步投影向量单位化。特征向量本身长度是可以任意伸缩的如果不做单位化一个范数为100的特征向量会在投影坐标中完全压制范数为0.1的特征向量。把每一列都除以其范数才能保证不同判别方向在坐标里权重一致。我自己的经验是这三步缺一不可。第一次跑NDGE时我只做了特征标准化没做对称归一化结果类内图稀疏的样本几乎全部被分错补上对称归一化之后准确率直接提了好几个百分点。3. Matlab代码实现一张图走完整个流程3.1 总体脚本结构训练、扫描、概率输出先说明一下整个代码的组织方式。NDGE不是单靠一个函数就能落地的东西我会拆成两部分compute_NDGE.m输入训练特征和标签输出投影矩阵 (W)、特征值、训练集的标准化参数ndge_dim_scan.m在不同维度下计算分类准确率并输出最优维度下每个测试样本的故障模式概率这样设计的好处是维度扫描不需要重复写构图和特征分解逻辑只需要循环调用核心函数。测试集的标准化参数也必须来自训练集绝对不能直接对测试集单独做Z-score否则等于偷偷用了测试集信息。3.2 compute_NDGE核心函数代码下面是完整可用的核心函数。相比网上的零碎片段我加了很多边界处理注释也写得比较详细。function [W, eigenvals, mu, stdv] compute_NDGE(Xtr, Ytr, opt) % compute_NDGE 归一化判别图嵌入 % 输入: % Xtr : n x m 训练特征矩阵, 每行一个样本 % Ytr : n x 1 训练标签, 建议为1,2,...,C % opt : 结构体, 字段: % k_in 类内近邻数, 默认5 % k_out 类间近邻数, 默认10 % sigma 热核带宽, 默认2 % reg 正则系数, 默认1e-6 % dim 输出投影维度, 默认30 % 输出: % W : m x dim 投影矩阵 % eigenvals: dim x 1 特征值降序 % mu, stdv : 训练集各维均值与标准差, 测试集投影时必须使用 if ~isfield(opt, k_in), opt.k_in 5; end if ~isfield(opt, k_out), opt.k_out 10; end if ~isfield(opt, sigma), opt.sigma 2; end if ~isfield(opt, reg), opt.reg 1e-6; end if ~isfield(opt, dim), opt.dim 30; end % 第一步: 特征标准化 mu mean(Xtr, 1); stdv std(Xtr, 0, 1); stdv(stdv eps) 1; % 防止零方差特征 X (Xtr - mu) ./ stdv; n size(X, 1); m size(X, 2); % 第二步: 欧氏距离矩阵 X2 sum(X.^2, 2); D2 X2 - 2 * (X * X) X2; D2(D2 0) 0; D sqrt(D2); % 第三步: 构造类内图S_w和类间图S_b S_w zeros(n, n); S_b zeros(n, n); [~, idx] sort(D, 2, ascend); for i 1:n same_cnt 0; diff_cnt 0; for j 1:n-1 nei idx(i, j 1); % 跳过自己, 所以j1 d D(i, nei); w exp(-d^2 / (2 * opt.sigma^2)); if Ytr(nei) Ytr(i) same_cnt opt.k_in S_w(i, nei) w; S_w(nei, i) w; same_cnt same_cnt 1; elseif Ytr(nei) ~ Ytr(i) diff_cnt opt.k_out S_b(i, nei) w; S_b(nei, i) w; diff_cnt diff_cnt 1; end end end % 第四步: 拉普拉斯矩阵及对称归一化 Dw diag(sum(S_w, 2)); Db diag(sum(S_b, 2)); Lw Dw - S_w; Lb Db - S_b; inv_sqrt_Dw diag(1 ./ sqrt(max(diag(Dw), eps))); inv_sqrt_Db diag(1 ./ sqrt(max(diag(Db), eps))); Lw_norm inv_sqrt_Dw * Lw * inv_sqrt_Dw; Lb_norm inv_sqrt_Db * Lb * inv_sqrt_Db; % 第五步: 广义特征分解 Sb X * Lb_norm * X; Sw X * Lw_norm * X; reg_val opt.reg * trace(Sw) / m; Sw Sw reg_val * eye(m); [V, E] eig(Sb, Sw); [eigenvals, order] sort(diag(E), descend); W V(:, order(1:opt.dim)); % 第六步: 投影向量单位化 for i 1:size(W, 2) nrm norm(W(:, i)); if nrm eps W(:, i) W(:, i) / nrm; end end end这段代码我测试过多次细节上有几个容易踩坑的地方第一stdv(stdv eps) 1必须加。有的特征可能所有样本取值都一样比如某个布尔型的开关状态特征。如果不处理标准化时会除以0整个矩阵变NaN。第二距离矩阵用展开式计算不要用两层循环套平方。在 (n) 不大时无所谓但 (n) 超过几百以后两层循环会慢到让你怀疑人生。上面这个写法利用矩阵乘法一次算完整张距离矩阵时间复杂度来自矩阵乘法比朴素循环快很多。第三构图的时候我用的是一趟遍历同时填 (S_w(i,nei)) 和 (S_w(nei,i))保证对称性。如果你用最大对称化 (S \max(S, S))效果一样但要注意遍历顺序否则会丢边。3.3 维度扫描和准确率统计不同故障诊断任务里最优投影维度并不会刚好等于类别数减一所以我们需要跑一条“维度-准确率”曲线。代码很简单循环调用核心函数就行但有几个工程细节训练集和测试集的划分要固定不能每个维度换一套数据分类器建议固定否则维度的对比会被分类器差异污染。function [acc_dim, prob_test, W_best, dim_best] ndge_dim_scan(... Xtr, Ytr, Xte, Yte, opt, dims) % dims: 需要扫描的维度向量, 例如 5:2:30 % 输出: % acc_dim : 每个维度对应的测试集准确率 % prob_test: 最优维度下, 每个测试样本的概率矩阵(n_te x C) % W_best : 最优维度对应的投影矩阵 % dim_best : 最优维度 acc_dim zeros(length(dims), 1); for ii 1:length(dims) opt.dim dims(ii); [W, ~, mu, stdv] compute_NDGE(Xtr, Ytr, opt); Ztr (Xtr - mu) ./ stdv * W; Zte (Xte - mu) ./ stdv * W; mdl fitcecoc(Ztr, Ytr, ... Learners, svm, ... FitPosterior, true); pred predict(mdl, Zte); acc_dim(ii) mean(pred(:) Yte(:)); end [~, best_idx] max(acc_dim); dim_best dims(best_idx); % 用最优维度再计算一次, 得到投影矩阵和概率 opt.dim dim_best; [W_best, eigenvals, mu, stdv] compute_NDGE(Xtr, Ytr, opt); Ztr (Xtr - mu) ./ stdv * W_best; Zte (Xte - mu) ./ stdv * W_best; mdl fitcecoc(Ztr, Ytr, ... Learners, svm, ... FitPosterior, true); [~, ~, ~, prob_test] predict(mdl, Zte); end这里用fitcecoc主要是为了直接拿到多分类后验概率。FitPosterior参数设置为true后predict返回的第四个输出是 (n \times C) 的概率矩阵每一行的和最接近1每一列对应一个故障模式。如果不想依赖统计和机器学习工具箱也可以用mnrfit/mnrval做多项逻辑回归或者自己实现一个softmax层。我的建议是先跑通上面这版再根据你自己的工具箱情况替换。3.4 每个样本的故障模式概率输出有了ndge_dim_scan拿到每个测试样本的概率矩阵之后我们可以做一件比单纯准确率更有意义的事输出“每个样本对不同故障模式的概率”而不是只给一个硬标签。故障诊断场景里硬标签有时会误导人。比如某个样本被强制分为“滚动体故障”但它的概率分布可能是“滚动体故障50%保持架故障40%外圈故障10%”。这种信息如果丢掉对现场维修决策非常危险。下面这段代码展示了如何把人读的结果格式化输出[acc_dim, prob_test, W_best, dim_best] ndge_dim_scan(... Xtr, Ytr, Xte, Yte, opt, 5:2:30); fprintf(最优维度: %d\n, dim_best); fprintf(测试集准确率: %.2f%%\n, acc_dim(best_idx) * 100); fprintf(\n样本序号 | 预测标签 | 各故障模式概率\n); for i 1:size(prob_test, 1) [maxp, pred] max(prob_test(i, :)); fprintf(%6d | %6d | , i, pred); fprintf(%.3f, , prob_test(i, :)); fprintf(\n); end % 把第一个样本的概率分布画成柱状图 figure; bar(prob_test(1, :)); set(gca, XTickLabel, 1:size(prob_test, 2)); xlabel(故障模式); ylabel(后验概率);这段代码跑完你会得到三个关键结果W_bestNDGE最终投影矩阵acc_dim不同维度下的准确率曲线prob_test每个样本对不同故障模式的概率这三个输出正好对应标题里的需求。我建议把W_best保存成.mat文件因为现场部署时如果要在新采集的数据上投影必须用完全相同的投影矩阵而不是重新训练。4. 实测记录和调参坑点4.1 一组实验记录我在一个滚动轴承故障模拟台采集的数据上跑过这套流程故障模式包括正常、内圈故障、外圈故障、滚动体故障每类训练样本大约30组测试样本大约20组。特征总共取了62维包括时域统计量、频域统计量和少量小波包能量特征。序号投影维度测试集准确率训练耗时1278.9%1.2s2588.4%1.3s3891.7%1.4s41295.6%1.6s51696.2%1.8s62096.1%1.9s73094.3%2.3s85091.5%3.1s从这张表能看出很典型的规律维度并不是越高越好。一开始增加维度判别信息越来越完整准确率上升到了16维附近达到峰值继续加维度反而把噪声和高频细节也保留下来准确率慢慢下降。这就是过拟合在降维层面的体现。所以我在实际项目里不会只看那个“最高点”而是找一段“平台期”。如果12维到20维准确率都在96%左右我会选择12维或14维而不是20维。维度越低分类器越稳定对现场传感器噪声的容忍度越高。4.2 三个最影响结果的参数第一个是类内近邻数k_in和类间近邻数k_out。k_in太小类内图可能断成多个小岛屿投影后同类样本不能被很好地聚拢k_in太大会把不同类别的样本也拉进类内图判别信息被污染。我的经验是k_in5~10k_out10~15。如果你的某个故障模式样本特别少可以按每个类别分别选近邻而不是全局统一一个数。第二个是热核带宽sigma。这个参数直接决定边的权重大小。sigma太小绝大多数边权跌破机器精度等于没有图sigma太大所有边的权重都趋近于1热核函数退化成0/1权重。一个可靠的初值是把所有样本距离排序取距离中位数作为sigma然后在0.5倍到2倍之间扫描一次。第三个是正则系数reg。Sw矩阵求逆或广义特征分解时正则项太小会数值崩溃太大会把判别方向压平。我一般通过观察特征值数量级来设置如果特征值出现NaN或超过1e10的异常值就把reg调大一个数量级如果投影结果在训练集上准确率很低就把reg调小一点。4.3 常见报错和排查速查表我在跑这种代码时积累了一个速查表现在分享出来。现象可能原因解决办法特征值全是NaN距离矩阵出现NaN或者特征里有NaN先检查Xtr是否含缺失值用any(isnan(Xtr(:)))排查广义特征分解报奇异性Sw奇异正则项太小增大reg或者先对原始特征做一次PCA预降维特征向量符号翻转特征分解每次结果符号不同投影矩阵训练后锁定保存部署时不再重算训练中可固定符号为正测试集准确率远低于训练集测试集没有使用训练集的mu/stdv必须保存compute_NDGE输出的mu,stdvfitcecoc不支持后验概率分类器设置或版本问题用templateSVM显式构造或改用softmax逻辑回归内存不足距离矩阵 (n \times n) 过大用分块计算或者对超大训练集采用锚点图某一类准确率极低该类样本太少图结构不完整增加类内近邻数或对少数类过采样后再构图最后再分享一个我的使用习惯在做维度扫描的时候不要只在单一测试集上比较准确率。我会把训练/测试划分跑五到十次随机划分然后把每个维度的平均准确率和波动范围画成误差棒图。选维度时不追求单次最高点而是选平均准确率高且波动小的区间。这一步看起来很浪费计算时间却能让NDGE的最终投影矩阵在不同批次数据上稳定得多。还有一点compute_NDGE里的图构造对样本顺序不敏感但对标签质量很敏感。如果训练集里存在误标注样本它会同时出现在类内图和类间图里造成图结构混乱。所以训练NDGE之前最好先做一个简单的异常样本筛查把明显错标的数据清理掉。这个预处理步骤往往比后面调任何参数都值钱。
阅读完成 · 觉得有帮助？
