NDGE归一化判别图嵌入:滚动轴承故障诊断降维实战与Matlab实现 干故障诊断这行的朋友应该都有体会振动信号的特征好提但特征一多就麻烦。时域的均值、峰值、峭度、波形因子频域的幅值谱特征、边带能量再加上小波包分解出来的节点能量轻轻松松攒出三四十维。把这些全部丢给分类器很多时候准确率反而上不去——这就是典型的高维小样本问题。我这段时间在一个滚动轴承故障诊断的模拟项目里核心做的就是“归一化判别图嵌入”NDGE监督降维。这个方法的输出很明确一共三样东西NDGE的最终投影矩阵、不同投影维度下的分类准确率、每个测试样本对不同故障模式的概率。我直接用Matlab跑通了整套流程今天把原理、代码实现和踩过的坑一起整理出来给同样在做故障诊断、受困于特征维数爆炸的同行做个参考。这套代码量不大核心函数几十行就能写完但每一步背后的取舍值得好好说清楚。1. 故障诊断流程中降维环节为什么绕不开1.1 高维特征带来的真实麻烦先说一个我在模拟项目里反复遇到的现象。假设你从轴承振动信号里提取了40维特征每个工况只有几十个样本。样本数远小于特征维度时样本在高维空间里会变得非常稀疏样本之间的欧氏距离几乎失去区分意义。K近邻、支持向量机这类依赖距离的分类器在这种空间里很容易把训练集的个别噪声点当规律学进去测试集上一测就露馅。这不是分类器不行而是特征空间本身不适合判别任务。更麻烦的是这些特征之间有大量冗余。峰值因子和峰值高度相关峭度和标准差之间也互相纠缠直接用原始特征训练模型计算开销大、可解释性差还容易过拟合。所以在故障诊断的标准流程里降维从来不是锦上添花而是连接特征提取和分类器之间绕不开的一环。1.2 主流降维方法各自的“水土不服”我最早接触的是主成分分析PCA。PCA只找方差最大的方向完全不看标签用来做无监督压缩没问题但压缩后的特征不一定对分类有利。一个典型的例子是两类样本在某个方向上方差很大但类别完全重叠PCA一定会优先保留这个方向结果最利于重构的方向恰恰不含判别信息。线性判别分析LDA把标签信息用起来了但它的假设比较强要求每类样本近似服从单一高斯分布、各类协方差相差不大。实际故障数据往往不是这种规整分布同一类样本可能存在多个局部簇LDA的全局散布矩阵很容易把局部的判别结构抹平。局部保持投影LPP擅长保留局部流形结构但它完全不理会类标签同类近邻和异类近邻在目标函数里被一视同仁降维结果同样可能损失判别力。痛点其实就一句话能不能把“保留局部结构”和“利用类别信息”这两件事同时做掉。NDGE就是冲着这个目标去的。1.3 几种降维方法的直观对比方法是否监督是否保留局部结构是否显式最大化类间可分性对非高斯分布的容忍度PCA否否否一般LDA是否是较差LPP否是否较好NDGE是是是较好表格里最后一列“较好”是我在多次实验中得到的经验判断不是严格的理论结论。NDGE本质上是用类标号约束图结构所以它比LDA更宽容数据分布形态又比LPP多了一股“刻意拉开异类”的劲头。2. NDGE的核心思想两张图和一个比值目标函数2.1 类内图与类间图是怎么构造的图嵌入的思路是把每个样本看成图上的一个节点样本之间的相似关系用边来表示。NDGE要构造两张图。第一张叫类内图也叫固有图。对每个样本在同类样本里找出最近的k1个邻居连上边边的权重用热核函数计算w_ij exp(-||x_i - x_j||² / 2σ²)。这张图的作用是拉近同类样本——在降维后的空间里原本连在一条边上的样本应该尽量靠在一起。第二张叫类间图也叫惩罚图。对每个样本在异类样本里找出最近的k2个邻居同样按热核赋权。这张图起的作用正好相反是推远异类样本——原本有边相连的异类样本在低维空间里应该尽量分开。细看这两张图的设计就会发现它同时兼顾了两件事热核权重刻画的是样本间的局部几何关系这是从LPP那里继承来的而同类相连、异类分开的规则又带着明显的判别意图这是LDA那一路的思路。两张图共同决定投影方向这就是NDGE能兼顾局部保持和判别分离的根本原因。2.2 目标函数一个比值一次广义特征值求解有了两张图下一步是定义目标函数。投影后的样本写成y_i Wᵀx_i类内离散程度用固有图的拉普拉斯矩阵来度量类间分离程度用惩罚图的拉普拉斯矩阵来度量。NDGE最终要最大化下面这个比值J(W) tr(Wᵀ Xᵀ L_b X W) / tr(Wᵀ Xᵀ L_w X W)其中L_w是类内图的拉普拉斯矩阵L_b是类间图的拉普拉斯矩阵。分子衡量投影后异类样本被拉开的程度分母衡量同类样本被压缩的紧密程度。比值越大说明投影空间里“类间距离”相对“类内距离”越大判别效果越好。这个优化问题求解析解最后落成一个广义特征值问题Xᵀ L_b X w λ Xᵀ L_w X w把特征值从大到小排列取前p个特征向量拼成投影矩阵W就是NDGE的最终输出之一。每个特征值对应一个投影方向特征值越大说明这个方向对判别贡献越大。2.3 “归一化”这三个字为什么值钱NDGE里的Normalized我这套实现里体现在两个地方。第一个是数据标准化。计算距离之前先用zscore把每个特征变成零均值、单位方差避免量纲大的特征垄断距离度量。第二个是拉普拉斯矩阵的归一化处理。普通的拉普拉斯矩阵是L D - W其中D是度矩阵。归一化版本是L I - D^(-1/2) W D^(-1/2)。这两者的差别可以打个比方普通拉普拉斯像“按人头投票”一个簇里样本特别多它的意见在目标函数里就被放得很大归一化拉普拉斯像“按人均发言”每个簇的贡献相对均衡不容易让样本密度差异带偏整体投影方向。归一化还有一个实打实的好处它的特征值有界数值性质更好广义特征值问题求解起来更稳定对数据里的噪声没那么敏感。我在实验里对比过同一批数据用未归一化版本换一次训练集划分投影矩阵的抖动明显比归一化版本大。3. Matlab代码实现从特征矩阵到三个核心输出3.1 输入数据的组织方式NDGE的输入很干净就是一个特征矩阵X和一个标签向量labels。X按行放样本列是特征labels从1开始编号。至于特征怎么提取不在今天讨论范围内你可以用时域统计量、频域包络谱特征或者小波包能量特征拼起来。我代码里用了一个模拟数据来演示流程% 模拟数据4类故障每类60个样本40维特征 rng(2024); N 240; d 40; C 4; X randn(N, d); labels repelem([1;2;3;4], 60);实际工程里X换成你从振动信号里提取的真实特征矩阵就行后面的流程完全一样。注意标签一定要是整数类别编号后面构造类内类间图时全靠它来区分同类异类。3.2 NDGE核心函数完整实现下面是整个项目的核心函数我把关键步骤都写了注释function [W, eigVal] NDGE(X, labels, k1, k2, sigma) % 归一化判别图嵌入核心函数 % 输入 % X - N×d 特征矩阵每行一个样本 % labels - N×1 标签向量从1开始编号 % k1 - 类内图近邻数 % k2 - 类间图近邻数 % sigma - 热核宽度 % 输出 % W - d×d 投影矩阵各列按特征值从大到小排列 % eigVal - 对应的特征值 [N, d] size(X); % 第0步数据标准化量纲均衡非常重要 X zscore(X); % 第1步计算所有样本两两之间的距离 D pdist2(X, X); % 第2步构造类内邻接矩阵Ww和类间邻接矩阵Bb Ww zeros(N, N); Bb zeros(N, N); for i 1:N % 类内图同类样本中最近的k1个 sameIdx find(labels labels(i)); sameIdx(sameIdx i) []; [~, ord] sort(D(i, sameIdx)); nn1 sameIdx(ord(1:min(k1, numel(ord)))); for j nn1 Ww(i, j) exp(-D(i, j)^2 / (2 * sigma^2)); end % 类间图异类样本中最近的k2个 diffIdx find(labels ~ labels(i)); [~, ord2] sort(D(i, diffIdx)); nn2 diffIdx(ord2(1:min(k2, numel(ord2)))); for j nn2 Bb(i, j) exp(-D(i, j)^2 / (2 * sigma^2)); end end % 第3步度矩阵与归一化拉普拉斯矩阵 Dw diag(sum(Ww, 2)); Db diag(sum(Bb, 2)); sWw Dw^(-0.5) * Ww * Dw^(-0.5); sBb Db^(-0.5) * Bb * Db^(-0.5); Lw eye(N) - sWw; Lb eye(N) - sBb; % 第4步在数据空间里构造矩阵并求解广义特征值问题 A X. * Lb * X; % 类间分离项 B X. * Lw * X; % 类内紧致项 % B矩阵可能奇异加一小点正则项防止求解失败 B B 1e-6 * eye(d); [V, D] eig(A, B); eigVal diag(D); [eigVal, idx] sort(eigVal, descend); W V(:, idx); end这段代码里的关键点我拆开说。先看第2步的图构造。代码用两层循环找近邻逻辑很朴素但要注意排序时把自身样本先排除掉否则每个样本都会把自己当成最近的邻居把权重矩阵对角线污染掉。热核权重里D(i,j)用的是欧氏距离的平方sigma控制权重随距离衰减的快慢。再看第3步的归一化拉普拉斯。直接用D^(-0.5)WD^(-0.5)的写法在样本量不大时没有问题但要注意度矩阵Dw如果存在全零行也就是某个样本没有任何同类近邻取倒数的开方就会算出Inf。这种情况在类数少、样本密的数据里不常见但越小样本越容易碰到。最后看第4步的广义特征值问题。Matlab的eig(A,B)直接解Ax lambdaB*x返回的D是广义特征值对角阵。我之所以在B上加1e-6的单位阵正则项是因为XᵀLwX在特征维度高、样本量不足时几乎必然奇异不加正则项eig经常报错或者给出复数特征值。3.3 输出一最终投影矩阵W怎么用W是d×d的方阵每一列是一个投影方向列的顺序按广义特征值从大到小排好了。投影的时候用前p列就够了% 把全部样本投影到 p 维空间 p 12; Y X * W(:, 1:p); % Y 的每一行是原样本的p维投影特征这里有个细节值得注意。W的每一列并不是普通意义下的正交归一向而是满足广义正交条件w_iᵀ B w_j 0i≠j。所以你检查Wᵀ*W是不是单位阵时会发现它不是这是正常的不要当成代码bug去改。实际使用中投影后的Y可以直接替代原始特征X送进任何分类器。我做故障诊断时习惯用最近邻或者线性SVM效果都不错。W矩阵本身也可以保存下来作为特征提取模型的一部分部署到现场现场来了一组新数据先提特征再乘W最后进分类器决策。3.4 输出二不同维度的分类准确率怎么算NDGE的投影方向有判别能力但不是说p越大越好。不同维度下的准确率曲线就是用来确定最优p的依据。这段代码很简单maxDim min(d, 30); % 最多算到30维再高意义不大 accs zeros(1, maxDim); for p 1:maxDim projTrain Xtrain * W(:, 1:p); projTest Xtest * W(:, 1:p); % 用最近类中心分类避免引入分类器自身的超参数 pred classifyByNearestCentroid(projTrain, trainLabels, projTest); accs(p) mean(pred testLabels); end [bestAcc, bestP] max(accs); fprintf(最优维度: %d, 最优准确率: %.2f%%\n, bestP, bestAcc * 100);配合的分类函数function pred classifyByNearestCentroid(Xtrain, labels, Xtest) classes unique(labels); C numel(classes); mu zeros(C, size(Xtrain, 2)); for c 1:C mu(c, :) mean(Xtrain(labels classes(c), :), 1); end M size(Xtest, 1); pred zeros(M, 1); for i 1:M dists zeros(C, 1); for c 1:C dists(c) sum((Xtest(i, :) - mu(c, :)).^2); end [~, pred(i)] min(dists); pred(i) classes(pred(i)); end end为什么这里刻意选最近类中心做分类器而不是直接上SVM或者神经网络两个原因。第一最近类中心本身没有超参数准确率曲线完全反映NDGE投影的质量不会被分类器参数干扰。第二NDGE的目标函数本质上就是在“压缩同类中心距离、扩大异类中心距离”最近类中心规则和它的假设高度一致逻辑上自洽。3.5 输出三每个样本对不同故障模式的概率最后这个输出是我觉得最实用的一个。它不只给一个硬标签而是给出“这个样本有多大概率属于每一种故障模式”。我的做法是选定维度p后在投影空间里计算每个类别的质心然后用softmax把样本到各类质心的距离转成概率function proba softmaxProb(Ytrain, trainLabels, Ytest, tau) % tau 是温度参数控制概率分布的锐度 classes unique(trainLabels); C numel(classes); proba zeros(size(Ytest, 1), C); mu zeros(C, size(Ytest, 2)); for c 1:C mu(c, :) mean(Ytrain(trainLabels classes(c), :), 1); end for i 1:size(Ytest, 1) for c 1:C proba(i, c) exp(-sum((Ytest(i, :) - mu(c, :)).^2) / tau); end proba(i, :) proba(i, :) / sum(proba(i, :)); end endtau的初值可以取训练样本到自身类质心距离平方的平均值。如果tau太大所有概率会趋近均匀分布tau太小概率几乎是0和1失去了置信度参考价值。实际用的时候可以画出来看效果一般把tau设在平均距离平方的0.5到2倍之间就够用。这个概率矩阵是M×C的M是测试样本数C是故障类别数。每一行加起来等于1最大概率对应的类别就是该样本的判别结果而概率值本身可以作为置信度若最大概率只有0.4说明样本落在类边界附近需要人工复核。4. 一个完整的故障诊断实验实例4.1 实验配置说明为了验证流程我在模拟项目X里用了一组滚动轴承振动数据4种状态分别是正常、内圈故障、外圈故障、滚动体故障。每种状态80个样本每个样本提取35维特征包括时域统计量14个、频域特征10个、小波包节点能量11个。训练集占四分之三测试集占四分之一。NDGE的参数设置如下参数设定值样本总数320特征维度35故障类别数4训练 / 测试样本240 / 80类内近邻 k15类间近邻 k25热核宽度 sigma由训练集距离中位数自动估计4.2 不同维度下的准确率曲线跑完准确率扫描结果大致如下投影维度 p准确率 (%)278.75590.00893.751095.001296.251595.002093.7535原始特征91.25这张表很能说明问题。准确率并不是随着维度单调上升的最佳区间在10到12维附近。维度太低判别信息不够维度太高把噪声也放进来了过拟合开始占上风。原始35维特征的准确率只有91.25%反而比12维投影后的96.25%低这正是降维价值的直接证据。所以这条准确率曲线不是跑着玩的它直接决定了最终部署时p的取值。我在项目里会先把bestP找出来再用全部训练数据重新跑一遍NDGE把投影矩阵和概率模型固化下来供测试使用。4.3 概率矩阵的判读方法概率矩阵输出后随机抽了三个测试样本结果长这样样本编号P(正常)P(内圈故障)P(外圈故障)P(滚动体故障)判别结果测试#70.0020.9510.0310.016内圈故障测试#230.8760.0640.0450.015正常测试#410.2100.2380.3050.247边界样本前两行的置信度很高可以直接采信。第三行就比较有意思了最大概率只有0.305样本明显落在多个故障模式的重叠区域。实际运营中遇到这种情况我不会立刻下结论而是建议补充时频分析或者重新采样验证。硬标签模式会给一个确定的类但那种确定性往往是假的概率输出把不确定性摆在明面上反而更接近工程决策的真实需求。5. 调参、踩坑与常见问题排查5.1 热核宽度sigma图结构的灵魂参数sigma选不好NDGE的效果会急转直下。sigma太小热核权重迅速衰落到零图近乎不连通目标函数里全是数值噪声sigma太大权重趋近均匀局部结构全部丢失NDGE退化成一种没什么判别力的全局方法。我的做法是先用一个启发式初值计算同类近邻对距离的中位数让sigma等于这个中位数或它的两倍然后在这个值附近扫2到3个档观察准确率曲线的变化。不少文献把sigma写成一个固定的小数那是针对特定数据集调出来的直接抄到自己的数据上很容易翻车。5.2 近邻数k1、k2怎么定k1和k2的一般起点是3到7。k1太小类内图可能不连通目标函数求解不稳定k1太大会把较远的同类样本也拉进来局部结构的优势被稀释。k2同理但它的作用方向是推开异类k2太小时类间排斥力不足分隔不彻底。比较稳的经验法则是k1和k2取每类样本数的10%到20%。比如每类80个样本就从8到16这个区间去试。不过近邻数的敏感度远不如sigma只要不是太极端对准确率的影响通常在两个百分点以内不需要过度调优。5.3 忘记数据标准化一切权重都白搭这是我最常帮别人排查的问题。有些特征量纲差异极大峰值可以到几十峭度在3附近小波包能量可能上千。不标准化直接算距离量纲大的特征几乎主导了所有近邻关系的判定热核权重也跟着偏。NDGE里的拉普拉斯归一化能缓解一部分度数不平衡的问题但它管不了特征量纲这一步必须单独做。代码里用zscore有一点要注意训练集的均值和标准差要保存下来应用到测试集上而不是对测试集重新zscore。否则测试集的特征缩放和训练集不一致投影结果会系统性偏移。5.4 广义特征值求解失败的几类原因如果eig(A,B)报错或者特征值出现复数优先检查三件事。第一B矩阵是否奇异XᵀLwX在样本量不够时很容易奇异加正则项是常规操作。第二Lw和Lb的构造是否正确尤其是度矩阵有没有全零行。第三特征维度d太大时可以考虑先跑PCA把维度压到几十维再进NDGE我试过这种级联方案准确率损失很小但稳定性提升明显。5.5 常见问题速查表现象可能原因解决方法低维准确率上不去sigma取值不当图结构不合理用距离中位数启始扫描两三个sigma档高维准确率反而下降维度超过了有效信息容量开始过拟合用准确率曲线选最优p投影矩阵出现NaN数据里存在NaN或Inf清洗数据后再zscore检查特征矩阵eig报正定矩阵错误B矩阵奇异加1e-6正则项或先PCA再NDGE各类概率非常接近tau过大softmax被拉平用平均距离平方设定tau并适当缩小图不连通求解不稳定k1太小或单类样本太少增大k1或增加同类样本数量5.6 概率输出里tau的调法tau的选择没有绝对标准。我的经验是先计算出所有训练样本到自身类质心的距离平方均值把它作为tau_0然后跑几个tau值比如0.5tau_0、tau_0、2tau_0看概率分布的形态。做故障诊断时我们希望正常样本的概率输出接近“确定无疑是正常的”边界样本保留一定程度的模糊tau_0附近通常能满足这个需求。6. 我的一点实际体会我在几个故障诊断模拟项目里横向对比过PCA、LPP和NDGE最直观的感受是当各类样本分布非常紧、类间分离得很开时三种方法差距不大一旦样本分布有交叠、数据里噪声偏多NDGE的稳定性优势就体现出来了。它把局部流形的保持和类别的判别信息揉进了同一个目标函数这比LDA那种全局假设要皮实很多也比LPP那种完全不管标签的纯无监督方式更有指向性。概率输出是这个流程里我最舍不得丢掉的部分。实际操作里它把一个“判为什么故障”的结果变成了一组“有多大把握是哪种故障”的证据链。现场工程师看到高置信度可以放心检修看到低置信度会主动补做验证整个决策过程从“黑盒给答案”变成了“有依据的推断”这在实际项目里的价值比那一两个百分点的准确率提升要重要得多。如果后续要扩展我的建议有两个方向。一是把NDGE和深度特征提取结合起来先用自编码器或一维卷积网络从原始波形里学非线性特征再把这些特征送进NDGE做判别压缩效果往往能再上一个台阶。二是把概率输出和决策风险结合起来给不同误判设置不同代价用贝叶斯决策做在线报警逻辑这会让故障诊断系统更贴近实际运维场景。先把今天这套流程跑通拿到你自己的投影矩阵、准确率曲线和概率矩阵后面的事情就好办多了。