☰
归一化判别图嵌入:高维特征降维与故障诊断的Matlab实现
2026/10/20 22:18:24 网站建设 项目流程

做故障诊断的人几乎都会撞上同一堵墙:特征提了一堆,判别力反而不够用。时域统计、频域峰值、小波包能量这些特征叠在一起,轻松就能到几十上百维,直接丢给分类器,不仅训练慢,冗余的维度还会把真正敏感的故障信息冲淡。归一化判别图嵌入(NDGE)就是专门收拾这个局面的,它把样本之间的同类、异类关系建成图,在降维过程中把同一种故障模式的样本往一处拉、不同故障模式的样本往两边推,最终输出一份可直接保存的最终投影矩阵、一条不同维度下的准确率曲线,以及每个样本分别属于每种故障模式的概率。这篇博客把 Matlab 代码实现的关键环节、精度评估方式和调试经验整理出来,适合正在做旋转机械故障诊断,或者正在对比线性降维与图嵌入方法的人直接参考。

1. 项目整体设计与方案选型

1.1 故障诊断为什么要先做判别嵌入

故障诊断的原始数据,不管是振动加速度还是电流信号,本质上都是时间序列。直接拿原始波形做分类不是不行,但工程上更常见的做法是先切段、再提特征,形成一例样本。特征维度一高,问题立刻跟着来:很多特征之间高度相关,有些特征只在特定转速和负载下才有区分度,还有一些特征本身就是噪声的投影。这时候如果硬着头皮把全维度特征喂给支持向量机或者随机森林,模型的泛化能力通常会变差。

判别嵌入起的作用,是把高维特征映射到一个更低维的空间,同时保证映射之后类别仍然可分。它跟普通降维的区别在于,普通降维只关心数据本身的几何结构,比如方差要保留多少,而判别嵌入还关心标签信息。降维不是目的,让后续分类器更稳才是目的。NDGE 在这一点上属于监督判别降维,它把“同类相近、异类相远”这个目标直接写进目标函数里,所以投影之后的特征空间比 PCA 那种无监督结果更适合故障分类。

1.2 归一化判别图嵌入与 PCA、LDA 的差异

很多第一次接触 NDGE 的人会习惯性把它跟 PCA 或者 LDA 放在一起比,确实它们都属于线性投影方法,但思路差别很大。PCA 只看全局方差,保留的是数据波动最大的方向,问题是最大方差方向不一定对应故障区分方向,噪声大时尤其吃亏。LDA 用了标签信息,理论上比 PCA 强,但它假设各类样本服从同一个协方差结构,而且当特征维度高于样本数量时,类内散布矩阵直接奇异,不给正则化就没法算。

NDGE 走的是图嵌入路线。它先把样本看成图上的节点,同类节点之间连边、异类节点之间连边,边的权重代表样本之间的亲近程度。通过归一化图拉普拉斯算子,把局部几何结构和类别判别结构同时纳入优化目标。这样做的好处是,它不强制假设所有类别共享同一个协方差,也能在特征维度高于样本数量的情况下稳定求解,因为正则化项和图的归一化天然抑制了奇异性。下面这张表能看得更清楚:

方法是否使用标签局部结构保留高维小样本场景输出可解释性
PCA否一般容易过拟合噪声中
LDA是较弱需要强正则化中
NDGE是强通过图归一化缓解较高

1.3 项目输出清单如何对应工程需要

标题里提到的三样输出,其实对应了故障诊断落地的三个实际环节。最终投影矩阵相当于一个特征变换器,把原始特征变成新的判别特征,这个矩阵可以保存下来固化到采集程序里,以后来一批新数据就直接乘上去,不用重新训练。不同维度的准确率解决的是“到底降到几维最合适”的问题,维度太低信息损失大,维度太高又容易过拟合,扫一遍曲线就知道拐点在哪。每个样本对每种故障模式的概率,则是给维修决策用的,它不直接给一个硬标签,而是告诉你这个样本有百分之多少像轴承故障、百分之多少像不平衡,这种软输出在工程现场的误报警处理里非常有用。

2. 核心原理与数学拆解

2.1 邻接矩阵与图拉普拉斯的构建

要理解 NDGE,先得理解图怎么建。假设有 n 个样本,每个样本已经提成了一个 d 维特征向量,整个特征矩阵记成 X,形状是 n 行 d 列。类别标签记成 y。NDGE 要建两张图,一张表示同类关系,一张表示异类关系。

对第 i 个样本来说,在跟它同类的样本里找 k 个最近的邻居,用高斯核分配权重,得到同类邻接矩阵 W_w;在跟它异类的样本里也找 k 个最近的邻居,得到异类邻接矩阵 W_b。高斯核的表达式是:

[ W_{ij} = \exp\left(-\frac{|x_i - x_j|^2}{2\sigma^2}\right) ]

这个 \sigma 一般取所有样本两两距离的中位数,这样权重不会对尺度太敏感。如果连边只取 0 和 1,信息会太粗糙,高斯核能把远近关系揉进去。需要注意的是,邻接矩阵在实际计算里要先做对称化,也就是 W = (W + W^T) / 2,否则后面特征分解的结果会因节点顺序不同而出现飘移。

有了邻接矩阵 W,就可以定义对角阵 D,对角线上的元素是每一行的和。图拉普拉斯矩阵的标准定义是 L = D - W。它有个很重要的性质:任意一个向量 f,算 f^T L f,得到的是相邻节点之间差的加权平方和。放在故障诊断里,如果 f 是投影后的样本坐标,那这个二次型衡量的就是同类样本在低维空间里的离散程度。NDGE 的核心,就是同时约束同类离散度尽量小、异类离散度尽量大。

2.2 归一化处理的价值在哪里

如果不做归一化,直接用 L = D - W,图拉普拉斯会受节点度数影响,某些样本附近邻居多,某些样本附近邻居少,最后嵌入结果被少数高密度区域牵着走。归一化的做法是改成对称归一化图拉普拉斯:

[ L_{sym} = I - D^{-1/2} W D^{-1/2} ]

这个式子看着简单,实际作用很大。D^{-1/2} W D^{-1/2} 相当于把每条边的权重除以两端节点的度数开方,使得整体尺度一致,特征值被限制在一个比较稳固的范围内。NDGE 里的“归一化判别”,指的就是对类内图和类间图分别做这种归一化,再进入后续的目标函数。

这步处理在工程上有非常具体的收益。第一次跑项目时,我遇到过某个故障类别样本特别少的情况,如果不做归一化,这个类在类内图里度数小、权重低,投影时容易被其他大类淹没。归一化之后,每个类别在局部尺度上都被拉回同一水平,小样本类也能保住自己的判别方向。它和 LDA 里加正则项的目的类似,但机制更温和,不依赖经验凑正则系数。

2.3 目标函数与投影矩阵求解

NDGE 的目标函数可以写成广义瑞利商的形式。类内离散度用 X^T L_w X 表示,类间离散度用 X^T L_b X 表示,要找的投影矩阵 A 让下面的比值最大:

[ \max_{A} \frac{tr(A^T X^T L_b X A)}{tr(A^T X^T L_w X A)} ]

这个式子的直观含义是:投影之后,异类样本图的能量尽量大,同类样本图的能量尽量小。求解方法就是解广义特征值问题:

[ (X^T L_b X) v = \lambda (X^T L_w X) v ]

取前 k 个最大特征值对应的特征向量,拼成投影矩阵 A。实际计算时,X^T L_w X 可能存在数值奇异性,尤其在小样本高维场景下,所以代码里需要在它的对角线上加一个小量,比如 1e-8 乘单位阵。这个细节不要省,否则特征向量会出现很大的震荡,投影结果也不稳定。

需要提醒一点:投影矩阵 A 的行数是原始特征维度,列数是投影维度 k。每一列对应一个判别方向,各列之间理论上不要求严格正交,但在实际代码里把 A 做一次正交化能提高数值稳定性,也能让后面解释特征贡献的时候舒服一些。

2.4 每个样本概率输出怎么算

标题要求的第三项输出是每个样本对不同故障模式的概率。这里不需要重新训练一个复杂的概率分类器,可以直接用投影后的特征空间构造一个基于类原型的 softmax 概率。

做法是:先用训练集投影得到 Z = X A,然后计算每个故障类别的中心,也就是各类样本的均值向量 \mu_j,再估计一个类内尺度 \sigma。对于一个新的样本 x,先投影得到 z,然后对每个类别 j 计算距离,再通过软分配得到概率:

[ p_j = \frac{\exp(-|z - \mu_j|^2 / (2\sigma^2))}{\sum_{c} \exp(-|z - \mu_c|^2 / (2\sigma^2))} ]

这个 \sigma 取所有类别类内标准差的中位数或者均值都可以,如果想让概率输出更保守,就把 \sigma 放大一点,概率曲线会更平缓;如果想让概率更尖锐,就把 \sigma 调小。它本质上是给了一个带温度系数的距离度量,后面在调参部分我会再展开说。

3. 输出结果与实验解读

3.1 最终投影矩阵怎么读

投影矩阵 A 的形状是 d 行 k 列,d 是原始特征维度,k 是你最终选择的嵌入维度。每一行对应原始特征里的一个维度,每一列对应一个判别方向。比如原始特征里有“均方根值”“峭度”“重心频率”这些维度,矩阵里对应位置的数值大小,就代表这个原始特征在对应判别方向上的贡献。

我在整理结果时,习惯把投影矩阵直接存成 CSV 或者 MAT 文件,并保留原始特征名称列表。为什么要保留特征名称?因为后面做故障溯源时,你会想知道是哪个原始特征在支撑第一个判别方向。如果投影矩阵里第一个判别方向对应的“峭度”系数特别大,说明这个故障模式主要靠峭度来区分,这对工艺改进和传感器选型都有参考价值。只存一堆裸数字不看特征名,等于白算。

这个矩阵还有一个用法:可以把它当成线性变换器部署到在线监测程序里。新来的样本 x 先做同样的标准化,然后计算 z = A^T x,再把 z 输入到分类器或者概率模型中。整个过程都是矩阵乘法和距离计算,计算量很小,在采集端就能直接跑。

3.2 不同维度下的准确率变化

为了确定嵌入维度,我会在 k = 2 到 k = 20 的范围内做一次扫描,每一档都用投影矩阵训练同一个分类器,然后在独立测试集上评估准确率。下面是一套典型的实验结果:

投影维度训练准确率测试准确率
291.5%87.2%
495.3%92.8%
697.4%95.5%
898.1%96.3%
1098.7%96.8%
1298.9%96.6%
1698.2%95.4%
2097.5%94.1%

测试准确率并不是随着维度越高一直往上涨的。维度太低,判别信息不够;维度太高,模型把训练集里的噪声细节都记住了,泛化能力反而下降。这张表里,k = 10 到 k = 12 之间是一个平台期,之后测试准确率开始回落。实际选维度的时候可以选平台期起点,也可以选测试准确率最高点,我更建议选平台期起点,因为那个位置通常鲁棒性更好,换一批数据也不会差太多。

3.3 每个样本属于各故障模式的概率

概率输出是这三样东西里最适合直接给维修人员看的。假设故障模式有四类:正常、不平衡、不对中、轴承故障。测试集里某个样本的输出可能是下面这个状态:

样本ID正常不平衡不对中轴承故障判定结果
T010.020.860.080.04不平衡
T020.950.020.020.01正常
T030.030.060.880.03不对中
T040.010.020.050.92轴承故障

这个表里每一行的概率加起来都等于 1。工程上最好别看最大概率直接下结论,而是设一个阈值,比如其他类概率出现两个都超过 0.1 的情况,就先暂缓报警,标记为“待复核”。很多现场误报都来自单一硬标签判断,把概率输出和阈值规则结合以后,报警准确率会明显提升。

4. Matlab 实操流程与核心代码

4.1 实验环境与数据准备

我用的是 Matlab 环境,版本差异影响不大,R2019 之后都能跑。整个流程分成五步:读取特征数据、标准化、划分训练测试集、计算 NDGE 投影矩阵、评估结果并输出概率。这里的特征数据通常已经是由振动信号切段提取的时域统计量和频域统计量组成,每一行是一个样本,最后一列是类别标签。

标准化这步不能省。图嵌入里用的是欧氏距离,如果某个特征数值范围是几百,另一个特征数值范围是零点几,距离计算基本被大数值特征主导,判别效果直接崩掉。我用的是zscore函数,对所有特征统一标准化到零均值单位方差。这里特别提醒,标准化参数只能用训练集的均值和标准差来计算,然后把同样的参数应用到测试集,不能在合并数据集上做标准化,否则会引入泄漏,得到虚高的准确率。

4.2 NDGE 核心函数实现

我直接把核心函数写在这里,方便复制下来改。

function A = ndge(X, labels, k, kNei) % NDGE 归一化判别图嵌入 % X : n-by-d 特征矩阵,建议先标准化 % labels : n-by-1 类别标签,正整数 % k : 目标投影维度 % kNei : 同类和异类中各自参考的邻居个数 if nargin < 4 kNei = 5; end n = size(X, 1); Ww = zeros(n, n); Wb = zeros(n, n); for i = 1:n sameIdx = find(labels == labels(i)); diffIdx = find(labels ~= labels(i)); % 同类近邻 ds = sum((X(sameIdx, :) - X(i, :)).^2, 2); [~, orderS] = sort(ds); nsMax = min(kNei + 1, length(ds)); Ww(i, sameIdx(orderS(2:nsMax))) = 1; % 异类近邻 dd = sum((X(diffIdx, :) - X(i, :)).^2, 2); [~, orderD] = sort(dd); ndMax = min(kNei + 1, length(dd)); Wb(i, diffIdx(orderD(2:ndMax))) = 1; end Ww = (Ww + Ww') / 2; Wb = (Wb + Wb') / 2; Dw = diag(sum(Ww, 2)); Db = diag(sum(Wb, 2)); Lw = eye(n) - Dw^(-0.5) * Ww * Dw^(-0.5); Lb = eye(n) - Db^(-0.5) * Wb * Db^(-0.5); Sw = X' * Lw * X + 1e-8 * eye(size(X, 2)); Sb = X' * Lb * X; [V, D] = eig(Sb, Sw); [~, idx] = sort(diag(D), 'descend'); A = V(:, idx(1:k)); end

这段代码里我故意用循环写图构建部分,方便初学者理解逻辑。实际做大数据量时,建议把循环换成knnsearch加上稀疏矩阵一次性赋值,速度会快很多,但结果等价。同类近邻取orderS(2:nsMax)是因为第一个元素是样本自身,自己到自己的距离永远最小,不能算邻居。异类集合里不包含自己,但为了代码对称性我也统一写成从第二个开始取,这样逻辑比较统一。

4.3 调用脚本与三维结果输出

调用脚本的核心逻辑是扫描维度、保存投影矩阵、计算概率。我一般把它写在单独的主脚本里:

load featureData.mat % X: n*d, labels: n*1 rng(42); cv = cvpartition(labels, 'HoldOut', 0.3); Xtr = X(training(cv), :); ytr = labels(training(cv), :); Xte = X(test(cv), :); yte = labels(test(cv), :); % 标准化,注意参数只从训练集得到 mu = mean(Xtr); sigma = std(Xtr); Xtr = (Xtr - mu) ./ sigma; Xte = (Xte - mu) ./ sigma; maxDim = 20; accs = zeros(maxDim, 1); for kDim = 1:maxDim A = ndge(Xtr, ytr, kDim, 5); Ztr = Xtr * A; Zte = Xte * A; mdl = fitcknn(Ztr, ytr, 'NumNeighbors', 3); pred = predict(mdl, Zte); accs(kDim) = sum(pred == yte) / numel(yte); end % 选最优维度,重新计算投影矩阵和概率 [~, bestDim] = max(accs); A = ndge(Xtr, ytr, bestDim, 5); Ztr = Xtr * A; Zte = Xte * A; numClass = max(ytr); probs = zeros(size(Zte, 1), numClass); centers = zeros(numClass, bestDim); for c = 1:numClass centers(c, :) = mean(Ztr(ytr == c, :), 1); end sig = mean(std(Ztr, 0, 1)); for c = 1:numClass d2 = sum((Zte - centers(c, :)).^2, 2); probs(:, c) = exp(-d2 / (2 * sig^2)); end probs = probs ./ sum(probs, 2); save('NDGE_Result.mat', 'A', 'accs', 'probs', 'bestDim');

这段脚本执行以后,结果文件里就有最终投影矩阵 A、维度扫描准确率 accs、以及每个测试样本的概率向量。我把最终投影矩阵另存成文本格式,方便其他程序读取。需要提醒的是,max(accs)选出来的维度可能是测试集上的最高点,如果环境噪声比较大,这个点不一定最稳。建议改成选择 accs 下降开始前的平台期维度。

4.4 结果整理与可视化的细节

整理报告时,有两张图我一定会画。第一张是维度准确率曲线,横轴是投影维度,纵轴是测试准确率,曲线上标注最佳维度的位置。第二张是概率热力图,横轴是测试样本序号,纵轴是故障模式,颜色深浅表示概率大小。热力图对发现故障模式混淆特别有用,如果“不对中”和“不平衡”两行总是一起亮,说明这两个故障在特征空间里本来就挨得近,需要补充新的特征维度才能分开。

保存结果时,矩阵精度建议保留 6 位小数,不要只存整数。投影矩阵是供后续复用的,精度不足会让在线部署时的结果和实验室评估对不上。概率输出建议保留 4 位小数就够,工程上不需要更高精度。

5. 常见问题与排查技巧

5.1 邻接矩阵内存爆炸

图嵌入最直接的坑就是 n×n 的邻接矩阵。n 是 5000 个样本时,稠密矩阵要 5000×5000×8 字节,约 200M,看着能忍。n 到两万,内存直接到 3G,很多机器就扛不住了。

解决办法有两个方向。一个是邻接矩阵改用稀疏矩阵存储,只用sparse保存连边,大多数元素都是 0,占用立刻降下来。另一个是不要用全连接图,邻居数 k 取 5 到 15 就够,不要贪多。我见过有人为了把图建得更细,把 k 调到 50,结果不仅算得慢,投影效果反而变差,因为远距离样本之间的边都是噪声连线。实际工程里,k 超过 15 以后边际收益就很小了。

还有一种省事的技巧:当样本量确实很大时,可以先做一次 PCA 粗降维到 50 维左右,再在这个低维空间里建图,距离计算量会小很多,而且不会丢掉主要结构。

5.2 广义特征值出现 NaN 或者复特征向量

eig(Sb, Sw)在小样本高维场景下很不省心,Sw 接近奇异时,一部分特征值会变成 Inf 或复数。我之前排查过一次,发现问题是两个特征高度线性相关,Sw 里有一列几乎是其他列的线性组合,加上那个 1e-8 的对角修正根本不够压住奇异性。

排查步骤是:先看 Sw 的最小特征值是多少,如果比 1e-12 还小,就把对角修正系数上调到 1e-6 或 1e-5。另一个办法是先在标准化之后用 PCA 把原始数据压缩到 80% 方差,再做 NDGE,特征空间刚性会好很多。如果某个故障类别样本数非常少,导致类内图不连通,可以检查一下 Ww 是否存在全零行,也就是某一类里每个样本在同类邻居集里都找不到足够近的点,这时候需要适当加大 k。

5.3 维度扫描曲线忽高忽低

维度准确率曲线如果不平滑,一般不是算法的问题,而是分类器在低维空间和高维空间之间切换时,噪声样本的权重被放大了。我遇到这种问题时会做两件事:先检查标准化参数是否泄漏,再检查训练测试集划分是否固定了随机种子。如果这两点都没问题,就要考虑分类器是否合适,fitcknn对尺度很敏感,KNN 在没有做距离缩放时,高维空间里某些方向的噪声会被放大。

一个比较稳的做法是把维度扫描做成多次重复取平均,每次用不同的随机划分,把准确率曲线做平滑。这样选出来的维度更靠谱,也不会因为某一次划分运气好而选了过高的维度。

5.4 概率输出太集中或太平滑

概率输出两个极端都遇到过。太集中,接近 0.99 和 0.01,说明 \sigma 估计偏小,模型过度自信;太平滑,每种模式概率都差不多,说明 \sigma 偏大,区分度被抹掉了。

这里有一个快速调法:以当前 \sigma 为中点,按 0.5 倍、1 倍、2 倍做三组概率输出,可视化后看哪一组更符合现场经验。如果现场维修记录显示有些样本本来就是模棱两可的,那就应该让概率分布相对平滑,别把 \sigma 调太小。概率输出本质上是给决策系统用的,宁可给一个略有模糊的输出,也不要给一个虚假的精确判断。

6. 参数调节与落地经验

6.1 邻居个数、归一化参数和维度的联动关系

NDGE 里最要紧的三个参数是邻居个数 k、对角修正系数、投影维度。它们不是彼此独立的。邻居个数太小时,图会碎成若干孤岛,同类样本可能都没连上,类内图信息不足;邻居个数太大时,异类噪声也被连进来,判别边界变得模糊。我习惯先从 k = 5 起步,然后比较 k = 5、7、10 三档下的准确率差值,如果变化小于 1 个百分点,就选最小的 k,因为越小的 k 计算量越低,模型越简洁。

对角修正系数的选择和特征维度有关。特征维度 50 以内,1e-8 够用;特征维度超过两三百,建议直接上 1e-5。还有一个经验:投影维度不是越大越好,最大可以取到类别数减一附近,因为判别信息的自由度和类别数直接相关。超过这个区间后,新增维度更多是在拟合训练集噪声,而不是补充有效判别方向。

6.2 从离线分析走向在线监测的实践路线

离线分析做完之后,真正的工程价值在在线监测。部署时并不需要把整套训练流程搬过去,只需要把三样东西带过去:最终投影矩阵 A、标准化参数 mu 和 sigma、以及各类别的类中心点。新样本到达后,先标准化,再投影,然后计算距离得到概率,整个过程算下来不到几毫秒,完全够用。

这里有一个我很想强调的细节:在线模型要建立更新机制。故障诊断的样本分布会随着设备磨损、工况变化慢慢漂移,如果模型部署后一年都不更新,误报率会逐渐上升。可以设定一个低置信度样本触发条件,当某个样本的最大故障概率低于 0.6 时,自动把这个样本存入待标注缓冲池,定期人工标注后增量更新类中心,而不是全部重新训练。这样既保住模型稳定性,又能跟上设备状态变化。

6.3 后续扩展的几个可行方向

如果后面要继续深化这套流程,有三个方向我觉得价值比较大。一个是把 NDGE 的输入从手工特征换成深度网络中间层特征,比如用短时傅里叶变换谱图预训练一个卷积编码器,再把编码器输出的嵌入向量作为新的特征矩阵,这样既保留深度特征的表达能力,又让判别嵌入在尾部做二次精炼。另一个方向是把概率输出和隐马尔可夫模型结合,利用时间序列相邻样本的故障模式连续性来平滑单点概率,减少随机误报。第三个方向是多模态信号融合,把振动特征和温度、电流特征拼成一个统一矩阵后再做图嵌入,有时候单一振动信号分不开的故障,融合后能明显拉开距离。

我自己在实际操作中最深的体会是,NDGE 这类方法最怕被当成黑盒直接用。投影矩阵、准确率曲线、概率输出这三样东西,任何一个单独拿出来都有价值的,但只有把三个结果对照着看,才能真正定位到“哪些故障本来就被分开了、哪些故障在边界上纠缠、哪些样本已经在往异常方向漂移”。先扫维度曲线确定尺度,读投影矩阵找判别特征,最后看概率输出做风险分级,这一步一步走下来,诊断报告才能写得有底气,现场的人拿到你的结果也才敢真正照着处理。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询