做设备状态监测的工程师应该都有过这种体验传感器一多、特征一算几十上百维的数据摆在面前分类器反而不知道该听谁的。我前阵子处理一组齿轮箱振动信号提取了包括时域统计量、频域峰值、小波包能量在内共40维特征6类工况直接扔给线性判别器交叉验证准确率卡在91%上下怎么调都上不去。后来换成归一化判别图嵌入NDGE先降维再在低维空间里做分类准确率才明显松动。这篇就把NDGE的原理、Matlab实现、三个核心输出的解读和调试经验一次性说清楚给正在做故障诊断分类的同学一个能直接落地的参照。NDGE这个东西核心价值一句话就能讲完在降维的同时把不同故障模式之间的边界拉得更开把同一故障模式内部的样本聚得更紧。它跟PCA、LDA这类经典方法最大的区别是它用“图”来刻画样本之间的关系既用了样本的类别标签又考虑了样本在原始特征空间里的局部密度结构。所以它特别适合故障诊断里那种“类内变体多、类间差异细”的数据比如轴承的早期点蚀和剥落特征空间里往往靠得很近普通降维方法很容易把它们揉在一起。我写这篇文章不会只给一个能跑的代码包那没有意义。重点是把三件事拆开讲透第一NDGE的目标函数是怎么设计出来的它的“归一化”到底干了什么第二在Matlab里怎么一步步把它落地从构图到广义特征值分解再到三个输出文件——最终投影矩阵、不同维度准确率、每个样本的故障概率第三实际调试时最容易踩的坑是什么以及怎么根据概率输出做更细的故障判断。1. 故障诊断里的老问题特征越多分类反而越不准1.1 高维特征里的“信息冗余”和“判别冲突”做故障诊断的人应该都经历过一个反直觉的现象特征提得越多效果反而越差。有人提取时域里的均值、峰值、峭度又加频域里的重心频率、边带幅值最后几十个特征并在一起。你以为是信息丰富了实际上里面塞满了三样东西和故障类型完全无关的工况变量、彼此高度相关的冗余特征、以及少量真正有判别力的线索。我举一个生活里的类比。你要在一群人里分辨谁是运动员如果只告诉你身高、体重、鞋码其实判断已经比较靠谱但如果再把头发长度、指甲形状、手机品牌也当成判断依据反而会被这些无关信息带偏。高维特征就是这种“把噪音当信息”的典型场景。分类器在高维空间里很容易找到一些虚假的规律——比如某几个特征组合恰好能隔开当前这批样本但换一批数据就完全失效这就是过拟合。所以在故障诊断的标准流程里特征降维几乎是一个绕不开的环节。传统的做法是PCA它找的是“方差最大的方向”目标是尽可能保留数据本身的分散程度。但PCA有一个致命的盲区它不知道你要分的是哪几类故障它只管方差最大。有时候方差最大的方向恰好是工况变化的方向而不是故障类型的方向降完维故障信息反而被丢掉了。1.2 从PCA到图嵌入降维为什么不能只看方差LDA比PCA稍微聪明一点它知道利用类别标签目标是最小化类内距离、最大化类间距离。这个思路本身很对但它在小样本高维场景下有一个硬伤类内散布矩阵通常是奇异的求逆都不稳定Matlab里一跑就给你一堆NaN。图嵌入的思路是绕开了“全局散布矩阵”这一套改从局部关系入手。它把每个样本看成图里的一个节点把样本之间的近邻关系看成边。如果两个样本属于同一类且离得很近就连一条“类内边”如果属于不同类但离得很近就连一条“类间边”。降维的目标变成让类内边的两个端点尽量靠近让类间边的两个端点尽量拉远。这么做的好处是它不要求整个类服从一个高斯分布也不要求数据整体是线性的它只关心“每一条边两端的局部关系”。这就很符合故障数据的实际形态——同一种故障在早期、中期、后期的特征表现差异很大类内分布可能是多峰的、条状的硬套一个全局散布矩阵很容易失真而局部图结构能把这个形态描述得更好。1.3 NDGE和LDA、LPP的区别在哪里熟悉降维方法的人可能听说过LPP局部保留投影它也是图嵌入。LPP和NDGE的差别在于LPP是无监督的它只保持近邻结构不关心近邻的标签NDGE是有监督的它同时构造“类内近邻”和“类间近邻”两张图把判别信息直接织进图结构里。下面这张表可以快速看出几种方法的分工差异。方法是否用标签约束目标适合的故障诊断场景典型问题PCA否最大化全局方差数据探索、工况变化分析不关心类别边界LDA是最大化类间/类内散布比线性可分、样本量充足小样本高维时矩阵奇异LPP否保持邻域局部结构无标签异常检测没有利用故障类别NDGE是类内图拉近、类间图推远多分类故障、类间重叠大需要调构图参数NDGE里的“判别”二字本质就是LDA思想的图化版本LDA用全局均值来定义类内和类间散布NDGE用局部近邻连接来定义。这种局部化的好处前面已经说了代价则是多了一堆构图参数需要调比如近邻个数、核宽度、正则系数这也是我后面专门用一整章来讲调参的原因。2. 归一化判别图嵌入的目标函数三股劲儿拧成一股绳2.1 类内图与类间图图嵌入的构图逻辑先引入一套符号后面写代码和推导都用它。假设有n个样本每个样本是d维特征向量组成矩阵X维度是d×n。标签向量y的取值是1到CC是故障模式总数。第一步是构图。对每个样本xi在同类样本里找k1个最近邻在异类样本里找k2个最近邻。然后定义两个亲和矩阵类内亲和矩阵S_w如果xj和xi是同类近邻那么S_w(i,j) exp(-||xi-xj||² / 2σ²)否则为0。 类间亲和矩阵S_b如果xj和xi是异类近邻那么S_b(i,j) exp(-||xi-xj||² / 2σ²)否则为0。这里用的是高斯核作用是距离越近的样本亲和作用越强距离拉远之后亲和作用指数衰减。σ是核宽度控制“多远算近”这个参数很敏感后面会讲怎么选。为什么要用两个图而不是一个图因为两个图承担的任务相反类内图负责“吸引力”把同类样本往一起拉类间图负责“排斥力”把异类样本往两边推。降维的投影方向最终就是这两股力的折中结果。2.2 归一化到底归一的是什么度矩阵与Laplacian这是NDGE最容易让人迷糊的地方。很多人在代码里看到“归一化”三个字以为只是对数据做了标准化其实在NDGE里有两层归一化第一层确实是数据标准化第二层是图Laplacian的对称归一化这才是方法名字里“归一化”的主要来源。在图上定义度矩阵DD(i,i)等于第i个节点连出边的权重之和。通俗理解度矩阵表示“这个样本在图上被多少邻居以及多强的连接包围”。图Laplacian矩阵L D - S是图论里研究图结构的基本算子。对Laplacian做对称归一化之后得到L_norm D^(-1/2) L D^(-1/2)为什么要多此一举因为原始图中每个节点的度可能差异很大有的样本处在密集区连接一大堆近邻有的样本处在稀疏区身边几乎没人。如果不做归一化投影方向会被密集区样本牵着走稀疏区的判别结构则几乎被忽略。这个情况在故障诊断里很常见正常状态样本往往聚集一大堆各类早期故障样本则稀稀拉拉分布在边缘。归一化相当于给每个节点的作用力做了“按连接数摊销”让密集区和稀疏区对投影方向有近似平等的话语权。第二层归一化才是数据层面的也就是对原始特征的每个维度做z-score标准化。这一步很多人会忘忘了之后构图用的欧氏距离就会被量纲大的特征主导比如幅值特征动辄上千频谱能量特征可能是几万而某些无量纲指标只有个位数结果距离计算基本等于只看那几个量纲大的特征图中的近邻关系完全失真。2.3 目标函数与广义特征值分解的数学推导两张图和对应的归一化Laplacian都构造好之后需要定义一个目标函数。我们用Lw_norm表示类内图的归一化LaplacianLb_norm表示类间图的归一化Laplacian。要找的是一个投影矩阵W它把原始特征x投影成y W^T x维度从d降到m。投影之后希望类内图连接点尽量靠近类间图连接点尽量远离。这个想法落到公式上就是最小化J_w(W) tr(W^T X Lw_norm X^T W)最大化J_b(W) tr(W^T X Lb_norm X^T W)两个目标合并经典的写法是最大化两者的比值max J_b(W) / (J_w(W) αI)这里的I是单位矩阵α是一个很小的正数作用是防止分母奇异同时起到正则化作用。如果J_w是零矩阵那这个优化问题没有意义α就是兜底的。这个带约束的优化问题可以转化为广义特征值分解问题(X Lb_norm X^T) w λ (X Lw_norm X^T αI) w求出广义特征值和特征向量之后取最大的m个特征值对应的特征向量按列拼起来就得到投影矩阵W。特征值的大小代表对应投影方向在判别意义上的“能力”——特征值越大说明在这个方向上类间推开的效果相对类内聚合的效果越强也就是判别信息越集中。我在给刚接触这个方法的同学解释的时候喜欢用“两股力拔河”来类比。类间图是往外拉的那股力类内图是往里拽的那股力。广义特征值分解就是找一个方向让往外拉的力最大、往里拽的力最小。每个投影方向的“特征值”就是这条方向上两股力量对比的比值把特征值从大到小排个序你就知道该往哪些方向上看。3. Matlab实现落地从训练矩阵到三个输出文件3.1 先定输入特征矩阵和标签的格式约定在动手写代码之前先把接口想清楚后面用起来会省很多事。我把主函数定成下面这样function [W, accs, probs] NDGE_demo(Xfeat, labels, opts) % 功能归一化判别图嵌入用于故障诊断降维与分类 % 输入 % Xfeat : nSample x dFeature 的特征矩阵每一行是一个样本 % labels : nSample x 1 的标签向量取值为 1,2,...,C % opts : 结构体包含 k1, k2, sigma, alpha, maxDim % 输出 % W : dFeature x maxDim 的投影矩阵按判别力从高到低排列 % accs : maxDim x 1 的准确率向量accs(m) 表示降到 m 维的交叉验证准确率 % probs : nSample x C 的概率矩阵每行是该样本对每一类故障的概率这个函数签名本身就对应了标题里的三个输出最终投影矩阵W、不同维度的准确率accs、每个样本对不同故障模式的概率probs。用户只需要喂进来特征矩阵和标签剩下的全部由函数内部处理。3.2 构图和亲和矩阵的代码骨架第一步永远是对特征做标准化再算两两距离。标注代码之前提醒一下计算n×n的欧氏距离矩阵在样本量很大的时候会占内存n5000时矩阵就有约200MB所以要留意样本规模。X zscore(Xfeat); % 第一层归一化均值0标准差1 n size(X, 1); % 样本数 d size(X, 2); % 原始特征维度 k1 opts.k1; % 同类近邻个数 k2 opts.k2; % 异类近邻个数 sigma opts.sigma; % 高斯核宽度 alpha opts.alpha; % 正则系数 Sw zeros(n, n); % 类内亲和矩阵 Sb zeros(n, n); % 类间亲和矩阵 % 逐样本构建局部近邻关系 for i 1:n other setdiff(1:n, i); dist2 sum((X(other, :) - X(i, :)).^2, 2); [~, ord] sort(dist2); knn_idx other(ord(1:max(k1, k2))); % 取前max(k1,k2)个候选邻居 % 同类近邻连边权重用高斯核 cnt 0; for jj knn_idx if labels(jj) labels(i) cnt k1 Sw(i, jj) exp(-dist2(other jj) / (2 * sigma^2)); cnt cnt 1; end end % 异类近邻连边 cnt 0; for jj knn_idx if labels(jj) ~ labels(i) cnt k2 Sb(i, jj) exp(-dist2(other jj) / (2 * sigma^2)); cnt cnt 1; end end end % 强制对称保证后续特征分解稳定 Sw max(Sw, Sw); Sb max(Sb, Sb);这里的max操作在Matlab里的写法是Sw max(Sw, Sw)等价于取上三角和下三角的最大值构造对称矩阵。为什么要对称因为后边的图Laplacian和广义特征值分解都要求矩阵是对称的如果构图过程因为数值精度出现不对称就会在后面冒出复数特征值这是非常常见的坑。3.3 广义特征分解求解投影矩阵W有了亲和矩阵下一步转成度矩阵和Laplacian再做对称归一化最终求解广义特征值问题。% 度矩阵与Laplacian Dw diag(sum(Sw, 2)); Db diag(sum(Sb, 2)); Lw Dw - Sw; Lb Db - Sb; % 对称归一化 Laplacian inv_sqrt_dw diag(1 ./ sqrt(sum(Sw, 2) eps)); inv_sqrt_db diag(1 ./ sqrt(sum(Sb, 2) eps)); Lw_norm inv_sqrt_dw * Lw * inv_sqrt_dw; Lb_norm inv_sqrt_db * Lb * inv_sqrt_db; % 把图散度转回特征空间 P_b X * Lb_norm * X; % d x d P_w X * Lw_norm * X; % d x d P_b (P_b P_b) / 2; P_w (P_w P_w) / 2; % 求解广义特征值问题 [V, D] eig(P_b, P_w alpha * eye(d)); eigvals real(diag(D)); [~, sortIdx] sort(eigvals, descend); W_all real(V(:, sortIdx));需要注意Matlab的eig函数用于广义特征值分解时A和B都应该是方阵且B尽量可逆。这里B已经加了αI也就是α乘以单位阵所以数值上通常是可逆的。如果你用的是较小的数据集不要用eigs来求前几个特征因为eigs默认用迭代法对非对称问题或矩阵条件数大的情况容易不收敛直接在完整特征分解里挑前m个反而更稳。3.4 不同维度准确率怎么算循环降维加分类器投影矩阵W_all算出来后从1维到maxDim每一档都要做一次交叉验证实验得到准确率曲线。这是NDGE实验里最直观的一张图它告诉你降维到多少维信息最划算。maxDim min(d, max(1, opts.maxDim)); maxDim min(maxDim, C - 1); % 通常不超过类别数减1 accs zeros(maxDim, 1); rng(1); % 固定随机种子保证结果可复现 cvp cvpartition(labels, KFold, 5); % 5折交叉验证 for m 1:maxDim Proj X * W_all(:, 1:m); % 降维后的特征 acc_fold zeros(cvp.NumTestSets, 1); for f 1:cvp.NumTestSets trIdx cvp.training(f); teIdx cvp.test(f); % 用最简单的K近邻分类器k5 model fitcknn(Proj(trIdx, :), labels(trIdx), NumNeighbors, 5); pred predict(model, Proj(teIdx, :)); acc_fold(f) sum(pred labels(teIdx)) / numel(teIdx); end accs(m) mean(acc_fold); end % 找最优维度 [bestAcc, bestDim] max(accs); fprintf(最优维度: %d, 交叉验证准确率: %.4f\n, bestDim, bestAcc);维度上限设成类别数减1是一般做法因为判别信息在大多数情况下集中在前C-1个方向上继续增加维度只会让分类器在有限样本下更容易过拟合。当然这也不是绝对的如果你有大量样本或者特征内部结构很复杂可以适当放宽上限但至少要把C-1当成分界线去观察准确率曲线的走向。3.5 每个样本对不同故障模式的概率输出方案概率输出是整个标题里最容易被忽略但实际价值最高的部分。降维完成后要给出每个样本属于每一类故障的概率可以用一个软的分类器而不是只给一个硬分类结果。Matlab里最简单的做法是用多分类逻辑回归也就是softmax回归。bestProj X * W_all(:, 1:bestDim); % 方案一fitcecoc 逻辑回归学习器可以输出后验概率 mdl fitcecoc(bestProj, labels, ... Learners, templateLinear(Learner, logistic), ... FitPosterior, true); [~, score] predict(mdl, bestProj); % score 是 nSample x C 的后验概率矩阵每一行之和接近1 probs score; % 对概率做一个简单校准防止模型过于自信 T 2; % 温度系数 probs exp(score / T) ./ sum(exp(score / T), 2);温度系数T是一个老概念softmax里面的“温度”其实就是把分数除以一个常数T再归一化。T越大输出的概率分布越平缓越能反映“不确定”的状态。我一般取T1.5到3之间具体看你对误诊的容忍度。如果T1模型往往会对错误类别也给出很高的自信度这在故障诊断里是很危险的——一个0.98的概率如果本身是错的就可能直接导致错误停机或者错过隐患。4. 读懂实验输出投影矩阵、准确率曲线和概率结果各自在说什么4.1 投影矩阵W的物理含义与可视化检查投影矩阵W是一个d×m的矩阵它每一列都是一个d维向量含义是“在这个判别方向下每个原始特征被分配了多少权重”。这句话展开来说如果W中第一列的第5个元素数值很大说明原始特征里的第5个维度对这个判别方向的贡献很大如果某一行的所有列数值都很小说明这个原始特征在NDGE眼里基本没什么判别力。这是一个非常实用的副产品——做特征筛选。我在实际项目里经常先跑一遍NDGE然后看W矩阵里每一行元素的绝对值之和按降序排序排在前面的特征就是有效的故障敏感特征排在后面的可以直接砍掉。这样做的好处是它比单独算每个特征和标签的相关系数更靠谱因为它考虑的是特征之间的组合关系而不是单打独斗。拿到W之后还要做一步可视化验证把原始特征投影到前两维或者前三维画散点图用不同颜色标不同故障。这一步虽然不直接出现在标题的三个输出里但我强烈建议每个做NDGE的人都要做因为准确率是个标量它只能告诉你“好不好”却没法告诉你“为什么好”。散点图能直接看到各类别的重叠程度、异常点位置、聚类形状这些信息对下一步调参和故障机理分析都有帮助。4.2 维度准确率曲线先升后降的另一面很多人看到“不同维度的准确率”这个输出只关心最高点在哪这是一个误区。准确率曲线的形状信息量远超那个最高点。理想情况下曲线会先快速上升因为前几个判别方向把主要类间结构带出来了然后进入平台期准确率基本不再变化接着可能缓慢下降因为后面特征向量对应的特征值很小里面装的更多是噪声分类器在这些维度上的拟合只会加重过拟合。这条曲线的第一个拐点通常对应数据的本征判别维度。比如在40维特征里如果曲线在第5维就进入平台期说明有效判别信息大概集中在5个方向上其他35维基本是冗余。这个结论比单纯选一个“最优维度”更有价值它直接指导你后续特征工程的方向如果维度已经降到很干净还分不开某些故障那问题不在降维而在原始特征本身没有捕获到那些故障的差异。还有一点要提醒如果准确率曲线从头到尾都没有明显上升甚至在低维时非常低那说明要么构图参数出了问题要么原始特征本身质量太差。这时候不要继续堆维度先回头检查第5章的调试项。4.3 样本概率是“软诊断”可信度与误诊排查每个样本对各类故障的概率是把“硬分类”升级成“软诊断”的关键。硬分类只会告诉你“这个是第3类故障”软概率会告诉你“这个有60%可能第3类30%可能第4类10%可能第2类”。在后者的信息量里你能做三件额外的事情。第一发现可疑样本。如果某个样本对某一类的概率接近1那它基本是稳定样本。如果某个样本的最高概率只有0.5左右说明它落在两个类的边界区域这种样本可能是早期故障的过渡形态也可能是特征缺失导致的信息不足需要重点排查。第二融合多源信息。振动、电流、温度多个传感器各来一组特征每组特征各自跑一遍NDGE得到各自的概率输出最后用贝叶斯融合或者简单的加权平均比直接拼接所有特征后一次性降维要稳得多。这里的前提就是每个通道必须先给出概率而不是硬标签。第三指导主动补充实验。概率均分在多个类之间的样本恰恰是你最应该关心的“灰色地带”。在实际诊断里我会专门把这些样本挑出来去看原始信号波形或者补做一个变工况实验看它到底随工况漂移还是稳定地偏向某一类。这种做法比单纯追求平均准确率更能解决工程里的实际问题。5. 复现和调参中的几个坎这些坑网上不会写5.1 特征尺度敏感为什么先标准化前面提过数据标准化但这里必须再强调一次因为它是所有问题里最容易犯、也最隐蔽的一个。我见过一些人直接在原始特征上跑NDGE然后发现准确率比PCA还差跑过来问怎么回事一检查特征里有一个“峰值因子”的范围是2到8另一个“均方根值”的范围是20到1500构图时欧氏距离几乎完全被均方根值掌控峰值因子那点判别信息全被淹没了。解决办法很简单训练前务必执行zscore对所有特征列都归一到均值0、标准差1。如果特征本身存在明显的离群点可以换成RobustScaler之类的稳健标准化按中位数和四分位距缩放。但最基础的还是zscore闭眼先用它。5.2 类间图近邻个数k2过大导致的问题k2是异类近邻个数它控制“排斥力”的来源数量。这个参数很有意思太小了排斥力样本太少投影方向可能只分得开最近的那几个类分不开全局太大了会把一些本来就远的异类样本也拉进来当“近邻”排斥力矩阵里塞满噪声最后投影方向被各种随机连接带偏。我个人的经验是k2取3到8之间且不要超过每类样本数的三分之一。如果一类故障只有10个样本k2取10以上就要小心——它可能把整个异类样本团都给拉进来了判别边界会被整体位置主导局部结构全丢。k1的取值范围类似但k1稍微大一点问题不大同类近邻多连几条边不会引入太多噪声只会在数值上稀释掉最核心的“紧邻关系”。一般的起手式是k15、k25然后根据类间重叠情况再调。5.3 小样本高维下的正则选择故障诊断里样本量通常很紧张特别是早期故障样本能凑到每类几十个就算不错。这种情况下d n也就是特征维度远大于样本数P_w矩阵一定奇异广义特征值分解会非常不稳定甚至出现复数特征值。解决思路有两条。第一条是正则化opts.alpha不要太小在特征数几十、样本数几百的时候alpha取1e-4到1e-2之间都算合理区间具体可以用交叉验证在小范围里搜索一下。第二条是先用PCA把维度压到20维以内再跑NDGE相当于一个两阶段降维。两阶段的做法在工程上更稳妥因为PCA只做无损压缩不会破坏前20个主成分里的判别结构。还有一种是血泪教训特征值出现了复数别慌先检查P_w加正则之后是否依然对称。我的代码里专门加了(P_wP_w)/2这步对称化处理就是在N次遇到复数特征值之后总结出来的。如果你拿到网上的代码里面没有这一步建议直接补上。5.4 概率输出的校准问题概率输出还有一个容易忽略的问题分类器给的概率往往没有经过校准也就是说“0.9的概率”并不意味着真有90%的正确率。特别是用fitcecoc的时候不同学习器之间的分数可能不处于同一量纲直接转成概率会失真。我建议在得到score之后做一次简单的温度缩放或者Platt校准。Platt校准用逻辑回归把概率和真实标签之间再拟合一次Matlab里可以用fitglm做一维校准。温度缩放更简单前面代码里那个T就是干这个的。在实际调的时候我会先看验证集上的可靠性图——把预测概率分桶每个桶里算实际正例占比如果点都落在大约45度对角线上说明校准良好如果两头都太自信就加大T。我跑过几套振动故障数据之后最深的体会是NDGE这类图嵌入方法难点从来不在数学推导而在构图参数是否和数据本身的流形结构匹配。参数对不上准确率就是上不去参数一但合适效果往往比瞎调PCA、LDA好出一大截。推荐拿到代码后先不急着上自己的数据把维度准确率曲线和第一阶梯形近邻的散点图打出来看一眼确认图里两类界线清晰、投影矩阵没有明显集中在某几个原特征上再正式套到故障诊断流程里。这个“先看结构、再跑分类”的习惯能帮你省下大量在参数泥潭里打转的时间。 
长河君 · 2026/10/21 15:57:18 · 阅读 1.2k NDGE归一化判别图嵌入在故障诊断中的原理与Matlab实现
做设备状态监测的工程师应该都有过这种体验传感器一多、特征一算几十上百维的数据摆在面前分类器反而不知道该听谁的。我前阵子处理一组齿轮箱振动信号提取了包括时域统计量、频域峰值、小波包能量在内共40维特征6类工况直接扔给线性判别器交叉验证准确率卡在91%上下怎么调都上不去。后来换成归一化判别图嵌入NDGE先降维再在低维空间里做分类准确率才明显松动。这篇就把NDGE的原理、Matlab实现、三个核心输出的解读和调试经验一次性说清楚给正在做故障诊断分类的同学一个能直接落地的参照。NDGE这个东西核心价值一句话就能讲完在降维的同时把不同故障模式之间的边界拉得更开把同一故障模式内部的样本聚得更紧。它跟PCA、LDA这类经典方法最大的区别是它用“图”来刻画样本之间的关系既用了样本的类别标签又考虑了样本在原始特征空间里的局部密度结构。所以它特别适合故障诊断里那种“类内变体多、类间差异细”的数据比如轴承的早期点蚀和剥落特征空间里往往靠得很近普通降维方法很容易把它们揉在一起。我写这篇文章不会只给一个能跑的代码包那没有意义。重点是把三件事拆开讲透第一NDGE的目标函数是怎么设计出来的它的“归一化”到底干了什么第二在Matlab里怎么一步步把它落地从构图到广义特征值分解再到三个输出文件——最终投影矩阵、不同维度准确率、每个样本的故障概率第三实际调试时最容易踩的坑是什么以及怎么根据概率输出做更细的故障判断。1. 故障诊断里的老问题特征越多分类反而越不准1.1 高维特征里的“信息冗余”和“判别冲突”做故障诊断的人应该都经历过一个反直觉的现象特征提得越多效果反而越差。有人提取时域里的均值、峰值、峭度又加频域里的重心频率、边带幅值最后几十个特征并在一起。你以为是信息丰富了实际上里面塞满了三样东西和故障类型完全无关的工况变量、彼此高度相关的冗余特征、以及少量真正有判别力的线索。我举一个生活里的类比。你要在一群人里分辨谁是运动员如果只告诉你身高、体重、鞋码其实判断已经比较靠谱但如果再把头发长度、指甲形状、手机品牌也当成判断依据反而会被这些无关信息带偏。高维特征就是这种“把噪音当信息”的典型场景。分类器在高维空间里很容易找到一些虚假的规律——比如某几个特征组合恰好能隔开当前这批样本但换一批数据就完全失效这就是过拟合。所以在故障诊断的标准流程里特征降维几乎是一个绕不开的环节。传统的做法是PCA它找的是“方差最大的方向”目标是尽可能保留数据本身的分散程度。但PCA有一个致命的盲区它不知道你要分的是哪几类故障它只管方差最大。有时候方差最大的方向恰好是工况变化的方向而不是故障类型的方向降完维故障信息反而被丢掉了。1.2 从PCA到图嵌入降维为什么不能只看方差LDA比PCA稍微聪明一点它知道利用类别标签目标是最小化类内距离、最大化类间距离。这个思路本身很对但它在小样本高维场景下有一个硬伤类内散布矩阵通常是奇异的求逆都不稳定Matlab里一跑就给你一堆NaN。图嵌入的思路是绕开了“全局散布矩阵”这一套改从局部关系入手。它把每个样本看成图里的一个节点把样本之间的近邻关系看成边。如果两个样本属于同一类且离得很近就连一条“类内边”如果属于不同类但离得很近就连一条“类间边”。降维的目标变成让类内边的两个端点尽量靠近让类间边的两个端点尽量拉远。这么做的好处是它不要求整个类服从一个高斯分布也不要求数据整体是线性的它只关心“每一条边两端的局部关系”。这就很符合故障数据的实际形态——同一种故障在早期、中期、后期的特征表现差异很大类内分布可能是多峰的、条状的硬套一个全局散布矩阵很容易失真而局部图结构能把这个形态描述得更好。1.3 NDGE和LDA、LPP的区别在哪里熟悉降维方法的人可能听说过LPP局部保留投影它也是图嵌入。LPP和NDGE的差别在于LPP是无监督的它只保持近邻结构不关心近邻的标签NDGE是有监督的它同时构造“类内近邻”和“类间近邻”两张图把判别信息直接织进图结构里。下面这张表可以快速看出几种方法的分工差异。方法是否用标签约束目标适合的故障诊断场景典型问题PCA否最大化全局方差数据探索、工况变化分析不关心类别边界LDA是最大化类间/类内散布比线性可分、样本量充足小样本高维时矩阵奇异LPP否保持邻域局部结构无标签异常检测没有利用故障类别NDGE是类内图拉近、类间图推远多分类故障、类间重叠大需要调构图参数NDGE里的“判别”二字本质就是LDA思想的图化版本LDA用全局均值来定义类内和类间散布NDGE用局部近邻连接来定义。这种局部化的好处前面已经说了代价则是多了一堆构图参数需要调比如近邻个数、核宽度、正则系数这也是我后面专门用一整章来讲调参的原因。2. 归一化判别图嵌入的目标函数三股劲儿拧成一股绳2.1 类内图与类间图图嵌入的构图逻辑先引入一套符号后面写代码和推导都用它。假设有n个样本每个样本是d维特征向量组成矩阵X维度是d×n。标签向量y的取值是1到CC是故障模式总数。第一步是构图。对每个样本xi在同类样本里找k1个最近邻在异类样本里找k2个最近邻。然后定义两个亲和矩阵类内亲和矩阵S_w如果xj和xi是同类近邻那么S_w(i,j) exp(-||xi-xj||² / 2σ²)否则为0。 类间亲和矩阵S_b如果xj和xi是异类近邻那么S_b(i,j) exp(-||xi-xj||² / 2σ²)否则为0。这里用的是高斯核作用是距离越近的样本亲和作用越强距离拉远之后亲和作用指数衰减。σ是核宽度控制“多远算近”这个参数很敏感后面会讲怎么选。为什么要用两个图而不是一个图因为两个图承担的任务相反类内图负责“吸引力”把同类样本往一起拉类间图负责“排斥力”把异类样本往两边推。降维的投影方向最终就是这两股力的折中结果。2.2 归一化到底归一的是什么度矩阵与Laplacian这是NDGE最容易让人迷糊的地方。很多人在代码里看到“归一化”三个字以为只是对数据做了标准化其实在NDGE里有两层归一化第一层确实是数据标准化第二层是图Laplacian的对称归一化这才是方法名字里“归一化”的主要来源。在图上定义度矩阵DD(i,i)等于第i个节点连出边的权重之和。通俗理解度矩阵表示“这个样本在图上被多少邻居以及多强的连接包围”。图Laplacian矩阵L D - S是图论里研究图结构的基本算子。对Laplacian做对称归一化之后得到L_norm D^(-1/2) L D^(-1/2)为什么要多此一举因为原始图中每个节点的度可能差异很大有的样本处在密集区连接一大堆近邻有的样本处在稀疏区身边几乎没人。如果不做归一化投影方向会被密集区样本牵着走稀疏区的判别结构则几乎被忽略。这个情况在故障诊断里很常见正常状态样本往往聚集一大堆各类早期故障样本则稀稀拉拉分布在边缘。归一化相当于给每个节点的作用力做了“按连接数摊销”让密集区和稀疏区对投影方向有近似平等的话语权。第二层归一化才是数据层面的也就是对原始特征的每个维度做z-score标准化。这一步很多人会忘忘了之后构图用的欧氏距离就会被量纲大的特征主导比如幅值特征动辄上千频谱能量特征可能是几万而某些无量纲指标只有个位数结果距离计算基本等于只看那几个量纲大的特征图中的近邻关系完全失真。2.3 目标函数与广义特征值分解的数学推导两张图和对应的归一化Laplacian都构造好之后需要定义一个目标函数。我们用Lw_norm表示类内图的归一化LaplacianLb_norm表示类间图的归一化Laplacian。要找的是一个投影矩阵W它把原始特征x投影成y W^T x维度从d降到m。投影之后希望类内图连接点尽量靠近类间图连接点尽量远离。这个想法落到公式上就是最小化J_w(W) tr(W^T X Lw_norm X^T W)最大化J_b(W) tr(W^T X Lb_norm X^T W)两个目标合并经典的写法是最大化两者的比值max J_b(W) / (J_w(W) αI)这里的I是单位矩阵α是一个很小的正数作用是防止分母奇异同时起到正则化作用。如果J_w是零矩阵那这个优化问题没有意义α就是兜底的。这个带约束的优化问题可以转化为广义特征值分解问题(X Lb_norm X^T) w λ (X Lw_norm X^T αI) w求出广义特征值和特征向量之后取最大的m个特征值对应的特征向量按列拼起来就得到投影矩阵W。特征值的大小代表对应投影方向在判别意义上的“能力”——特征值越大说明在这个方向上类间推开的效果相对类内聚合的效果越强也就是判别信息越集中。我在给刚接触这个方法的同学解释的时候喜欢用“两股力拔河”来类比。类间图是往外拉的那股力类内图是往里拽的那股力。广义特征值分解就是找一个方向让往外拉的力最大、往里拽的力最小。每个投影方向的“特征值”就是这条方向上两股力量对比的比值把特征值从大到小排个序你就知道该往哪些方向上看。3. Matlab实现落地从训练矩阵到三个输出文件3.1 先定输入特征矩阵和标签的格式约定在动手写代码之前先把接口想清楚后面用起来会省很多事。我把主函数定成下面这样function [W, accs, probs] NDGE_demo(Xfeat, labels, opts) % 功能归一化判别图嵌入用于故障诊断降维与分类 % 输入 % Xfeat : nSample x dFeature 的特征矩阵每一行是一个样本 % labels : nSample x 1 的标签向量取值为 1,2,...,C % opts : 结构体包含 k1, k2, sigma, alpha, maxDim % 输出 % W : dFeature x maxDim 的投影矩阵按判别力从高到低排列 % accs : maxDim x 1 的准确率向量accs(m) 表示降到 m 维的交叉验证准确率 % probs : nSample x C 的概率矩阵每行是该样本对每一类故障的概率这个函数签名本身就对应了标题里的三个输出最终投影矩阵W、不同维度的准确率accs、每个样本对不同故障模式的概率probs。用户只需要喂进来特征矩阵和标签剩下的全部由函数内部处理。3.2 构图和亲和矩阵的代码骨架第一步永远是对特征做标准化再算两两距离。标注代码之前提醒一下计算n×n的欧氏距离矩阵在样本量很大的时候会占内存n5000时矩阵就有约200MB所以要留意样本规模。X zscore(Xfeat); % 第一层归一化均值0标准差1 n size(X, 1); % 样本数 d size(X, 2); % 原始特征维度 k1 opts.k1; % 同类近邻个数 k2 opts.k2; % 异类近邻个数 sigma opts.sigma; % 高斯核宽度 alpha opts.alpha; % 正则系数 Sw zeros(n, n); % 类内亲和矩阵 Sb zeros(n, n); % 类间亲和矩阵 % 逐样本构建局部近邻关系 for i 1:n other setdiff(1:n, i); dist2 sum((X(other, :) - X(i, :)).^2, 2); [~, ord] sort(dist2); knn_idx other(ord(1:max(k1, k2))); % 取前max(k1,k2)个候选邻居 % 同类近邻连边权重用高斯核 cnt 0; for jj knn_idx if labels(jj) labels(i) cnt k1 Sw(i, jj) exp(-dist2(other jj) / (2 * sigma^2)); cnt cnt 1; end end % 异类近邻连边 cnt 0; for jj knn_idx if labels(jj) ~ labels(i) cnt k2 Sb(i, jj) exp(-dist2(other jj) / (2 * sigma^2)); cnt cnt 1; end end end % 强制对称保证后续特征分解稳定 Sw max(Sw, Sw); Sb max(Sb, Sb);这里的max操作在Matlab里的写法是Sw max(Sw, Sw)等价于取上三角和下三角的最大值构造对称矩阵。为什么要对称因为后边的图Laplacian和广义特征值分解都要求矩阵是对称的如果构图过程因为数值精度出现不对称就会在后面冒出复数特征值这是非常常见的坑。3.3 广义特征分解求解投影矩阵W有了亲和矩阵下一步转成度矩阵和Laplacian再做对称归一化最终求解广义特征值问题。% 度矩阵与Laplacian Dw diag(sum(Sw, 2)); Db diag(sum(Sb, 2)); Lw Dw - Sw; Lb Db - Sb; % 对称归一化 Laplacian inv_sqrt_dw diag(1 ./ sqrt(sum(Sw, 2) eps)); inv_sqrt_db diag(1 ./ sqrt(sum(Sb, 2) eps)); Lw_norm inv_sqrt_dw * Lw * inv_sqrt_dw; Lb_norm inv_sqrt_db * Lb * inv_sqrt_db; % 把图散度转回特征空间 P_b X * Lb_norm * X; % d x d P_w X * Lw_norm * X; % d x d P_b (P_b P_b) / 2; P_w (P_w P_w) / 2; % 求解广义特征值问题 [V, D] eig(P_b, P_w alpha * eye(d)); eigvals real(diag(D)); [~, sortIdx] sort(eigvals, descend); W_all real(V(:, sortIdx));需要注意Matlab的eig函数用于广义特征值分解时A和B都应该是方阵且B尽量可逆。这里B已经加了αI也就是α乘以单位阵所以数值上通常是可逆的。如果你用的是较小的数据集不要用eigs来求前几个特征因为eigs默认用迭代法对非对称问题或矩阵条件数大的情况容易不收敛直接在完整特征分解里挑前m个反而更稳。3.4 不同维度准确率怎么算循环降维加分类器投影矩阵W_all算出来后从1维到maxDim每一档都要做一次交叉验证实验得到准确率曲线。这是NDGE实验里最直观的一张图它告诉你降维到多少维信息最划算。maxDim min(d, max(1, opts.maxDim)); maxDim min(maxDim, C - 1); % 通常不超过类别数减1 accs zeros(maxDim, 1); rng(1); % 固定随机种子保证结果可复现 cvp cvpartition(labels, KFold, 5); % 5折交叉验证 for m 1:maxDim Proj X * W_all(:, 1:m); % 降维后的特征 acc_fold zeros(cvp.NumTestSets, 1); for f 1:cvp.NumTestSets trIdx cvp.training(f); teIdx cvp.test(f); % 用最简单的K近邻分类器k5 model fitcknn(Proj(trIdx, :), labels(trIdx), NumNeighbors, 5); pred predict(model, Proj(teIdx, :)); acc_fold(f) sum(pred labels(teIdx)) / numel(teIdx); end accs(m) mean(acc_fold); end % 找最优维度 [bestAcc, bestDim] max(accs); fprintf(最优维度: %d, 交叉验证准确率: %.4f\n, bestDim, bestAcc);维度上限设成类别数减1是一般做法因为判别信息在大多数情况下集中在前C-1个方向上继续增加维度只会让分类器在有限样本下更容易过拟合。当然这也不是绝对的如果你有大量样本或者特征内部结构很复杂可以适当放宽上限但至少要把C-1当成分界线去观察准确率曲线的走向。3.5 每个样本对不同故障模式的概率输出方案概率输出是整个标题里最容易被忽略但实际价值最高的部分。降维完成后要给出每个样本属于每一类故障的概率可以用一个软的分类器而不是只给一个硬分类结果。Matlab里最简单的做法是用多分类逻辑回归也就是softmax回归。bestProj X * W_all(:, 1:bestDim); % 方案一fitcecoc 逻辑回归学习器可以输出后验概率 mdl fitcecoc(bestProj, labels, ... Learners, templateLinear(Learner, logistic), ... FitPosterior, true); [~, score] predict(mdl, bestProj); % score 是 nSample x C 的后验概率矩阵每一行之和接近1 probs score; % 对概率做一个简单校准防止模型过于自信 T 2; % 温度系数 probs exp(score / T) ./ sum(exp(score / T), 2);温度系数T是一个老概念softmax里面的“温度”其实就是把分数除以一个常数T再归一化。T越大输出的概率分布越平缓越能反映“不确定”的状态。我一般取T1.5到3之间具体看你对误诊的容忍度。如果T1模型往往会对错误类别也给出很高的自信度这在故障诊断里是很危险的——一个0.98的概率如果本身是错的就可能直接导致错误停机或者错过隐患。4. 读懂实验输出投影矩阵、准确率曲线和概率结果各自在说什么4.1 投影矩阵W的物理含义与可视化检查投影矩阵W是一个d×m的矩阵它每一列都是一个d维向量含义是“在这个判别方向下每个原始特征被分配了多少权重”。这句话展开来说如果W中第一列的第5个元素数值很大说明原始特征里的第5个维度对这个判别方向的贡献很大如果某一行的所有列数值都很小说明这个原始特征在NDGE眼里基本没什么判别力。这是一个非常实用的副产品——做特征筛选。我在实际项目里经常先跑一遍NDGE然后看W矩阵里每一行元素的绝对值之和按降序排序排在前面的特征就是有效的故障敏感特征排在后面的可以直接砍掉。这样做的好处是它比单独算每个特征和标签的相关系数更靠谱因为它考虑的是特征之间的组合关系而不是单打独斗。拿到W之后还要做一步可视化验证把原始特征投影到前两维或者前三维画散点图用不同颜色标不同故障。这一步虽然不直接出现在标题的三个输出里但我强烈建议每个做NDGE的人都要做因为准确率是个标量它只能告诉你“好不好”却没法告诉你“为什么好”。散点图能直接看到各类别的重叠程度、异常点位置、聚类形状这些信息对下一步调参和故障机理分析都有帮助。4.2 维度准确率曲线先升后降的另一面很多人看到“不同维度的准确率”这个输出只关心最高点在哪这是一个误区。准确率曲线的形状信息量远超那个最高点。理想情况下曲线会先快速上升因为前几个判别方向把主要类间结构带出来了然后进入平台期准确率基本不再变化接着可能缓慢下降因为后面特征向量对应的特征值很小里面装的更多是噪声分类器在这些维度上的拟合只会加重过拟合。这条曲线的第一个拐点通常对应数据的本征判别维度。比如在40维特征里如果曲线在第5维就进入平台期说明有效判别信息大概集中在5个方向上其他35维基本是冗余。这个结论比单纯选一个“最优维度”更有价值它直接指导你后续特征工程的方向如果维度已经降到很干净还分不开某些故障那问题不在降维而在原始特征本身没有捕获到那些故障的差异。还有一点要提醒如果准确率曲线从头到尾都没有明显上升甚至在低维时非常低那说明要么构图参数出了问题要么原始特征本身质量太差。这时候不要继续堆维度先回头检查第5章的调试项。4.3 样本概率是“软诊断”可信度与误诊排查每个样本对各类故障的概率是把“硬分类”升级成“软诊断”的关键。硬分类只会告诉你“这个是第3类故障”软概率会告诉你“这个有60%可能第3类30%可能第4类10%可能第2类”。在后者的信息量里你能做三件额外的事情。第一发现可疑样本。如果某个样本对某一类的概率接近1那它基本是稳定样本。如果某个样本的最高概率只有0.5左右说明它落在两个类的边界区域这种样本可能是早期故障的过渡形态也可能是特征缺失导致的信息不足需要重点排查。第二融合多源信息。振动、电流、温度多个传感器各来一组特征每组特征各自跑一遍NDGE得到各自的概率输出最后用贝叶斯融合或者简单的加权平均比直接拼接所有特征后一次性降维要稳得多。这里的前提就是每个通道必须先给出概率而不是硬标签。第三指导主动补充实验。概率均分在多个类之间的样本恰恰是你最应该关心的“灰色地带”。在实际诊断里我会专门把这些样本挑出来去看原始信号波形或者补做一个变工况实验看它到底随工况漂移还是稳定地偏向某一类。这种做法比单纯追求平均准确率更能解决工程里的实际问题。5. 复现和调参中的几个坎这些坑网上不会写5.1 特征尺度敏感为什么先标准化前面提过数据标准化但这里必须再强调一次因为它是所有问题里最容易犯、也最隐蔽的一个。我见过一些人直接在原始特征上跑NDGE然后发现准确率比PCA还差跑过来问怎么回事一检查特征里有一个“峰值因子”的范围是2到8另一个“均方根值”的范围是20到1500构图时欧氏距离几乎完全被均方根值掌控峰值因子那点判别信息全被淹没了。解决办法很简单训练前务必执行zscore对所有特征列都归一到均值0、标准差1。如果特征本身存在明显的离群点可以换成RobustScaler之类的稳健标准化按中位数和四分位距缩放。但最基础的还是zscore闭眼先用它。5.2 类间图近邻个数k2过大导致的问题k2是异类近邻个数它控制“排斥力”的来源数量。这个参数很有意思太小了排斥力样本太少投影方向可能只分得开最近的那几个类分不开全局太大了会把一些本来就远的异类样本也拉进来当“近邻”排斥力矩阵里塞满噪声最后投影方向被各种随机连接带偏。我个人的经验是k2取3到8之间且不要超过每类样本数的三分之一。如果一类故障只有10个样本k2取10以上就要小心——它可能把整个异类样本团都给拉进来了判别边界会被整体位置主导局部结构全丢。k1的取值范围类似但k1稍微大一点问题不大同类近邻多连几条边不会引入太多噪声只会在数值上稀释掉最核心的“紧邻关系”。一般的起手式是k15、k25然后根据类间重叠情况再调。5.3 小样本高维下的正则选择故障诊断里样本量通常很紧张特别是早期故障样本能凑到每类几十个就算不错。这种情况下d n也就是特征维度远大于样本数P_w矩阵一定奇异广义特征值分解会非常不稳定甚至出现复数特征值。解决思路有两条。第一条是正则化opts.alpha不要太小在特征数几十、样本数几百的时候alpha取1e-4到1e-2之间都算合理区间具体可以用交叉验证在小范围里搜索一下。第二条是先用PCA把维度压到20维以内再跑NDGE相当于一个两阶段降维。两阶段的做法在工程上更稳妥因为PCA只做无损压缩不会破坏前20个主成分里的判别结构。还有一种是血泪教训特征值出现了复数别慌先检查P_w加正则之后是否依然对称。我的代码里专门加了(P_wP_w)/2这步对称化处理就是在N次遇到复数特征值之后总结出来的。如果你拿到网上的代码里面没有这一步建议直接补上。5.4 概率输出的校准问题概率输出还有一个容易忽略的问题分类器给的概率往往没有经过校准也就是说“0.9的概率”并不意味着真有90%的正确率。特别是用fitcecoc的时候不同学习器之间的分数可能不处于同一量纲直接转成概率会失真。我建议在得到score之后做一次简单的温度缩放或者Platt校准。Platt校准用逻辑回归把概率和真实标签之间再拟合一次Matlab里可以用fitglm做一维校准。温度缩放更简单前面代码里那个T就是干这个的。在实际调的时候我会先看验证集上的可靠性图——把预测概率分桶每个桶里算实际正例占比如果点都落在大约45度对角线上说明校准良好如果两头都太自信就加大T。我跑过几套振动故障数据之后最深的体会是NDGE这类图嵌入方法难点从来不在数学推导而在构图参数是否和数据本身的流形结构匹配。参数对不上准确率就是上不去参数一但合适效果往往比瞎调PCA、LDA好出一大截。推荐拿到代码后先不急着上自己的数据把维度准确率曲线和第一阶梯形近邻的散点图打出来看一眼确认图里两类界线清晰、投影矩阵没有明显集中在某几个原特征上再正式套到故障诊断流程里。这个“先看结构、再跑分类”的习惯能帮你省下大量在参数泥潭里打转的时间。