CN105760839A - 基于多特征流形学习与支持向量机的轴承故障诊断方法 - Google Patents

Abstract

本发明公开了一种基于多特征流形学习与支持向量机的轴承故障诊断方法，1)通过加速度传感器采集每种工况下滚动轴承在不同转速的振动加速度信号作为训练样本；2)对训练样本分别提取时域、频域和时频域特征参数；3)进行流行学习得到低维流形结构；4)通过传感器采集待测滚动轴承在转动时的振动加速度信号作为测试样本；5)提取测试样本的时域、频域和时频域特征参数；6)对测试样本进行流形学习得到低维流行结构；7)采用支持向量机分类方法对测试样本与训练样本进行匹配，将与测试样本最为匹配的训练样本所属的工况类别判定为测试样本的工况类别。本方法能够提高滚动轴承故障诊断的准确性和有效性。

Description

基于多特征流形学习与支持向量机的轴承故障诊断方法

技术领域

本发明涉及轴承机械故障诊断，尤其涉及一种基于多特征流形学习与支持向量机的滚动轴承故障诊断方法，属于机械故障诊断和计算机人工智能技术领域。

背景技术

滚动轴承是旋转机械中的关键零部件，具有摩擦小、精度高、成本低廉、互换性好等优点，在冶金、石油、化工、航天航空、煤炭电力等各个部门中得到广泛应用。但是，滚动轴承也是旋转机械中最容易损坏的零部件之一，滚动轴承抗冲击能力较弱，在冲击作用下很容易诱发故障产生。滚动轴承一旦发生故障极易导致整个机械系统瘫痪，因此对滚动轴承进行早期状态监测、分析与诊断在工程应用中具有重要意义。

滚动轴承故障诊断的关键是提取有效反应轴承运行状态与故障形式的特征参数。振动信号携带大量代表机械设备健康状况的信息，由此振动分析方法成为旋转机械故障诊断中最常用的方法之一。但滚动轴承早期故障信号微弱，加上工作环境中的噪声和其它干扰因素影响，信噪比很低，而且大多为非平稳非线性信号，特征提取困难很大，采用传统的振动分析方法难以做到可靠、高效的滚动轴承状态识别与故障诊断。

目前常用的机械故障特征提取方法主要有：傅里叶变换(FastFouriertransform，简称FFT)、稀疏表示提取瞬态特征、小波变换和经验模态分解(Empiricalmodedecomposition，简称EMD)、局部中值分解(Localmeandecomposition，简称LMD)、人工智能等。FFT方法能够直接从信号中找到感兴趣的频率成分，但是却无法同时兼顾信号在时域和频域中的全貌和局部化信息。稀疏表示提取瞬态特征，但是由于它的计算量大等特点，在一定程度上限制了它的发展。小波变换的结果较大程度上取决于小波基的选择，而通常小波基的选取是较难的。EMD方法和LMD方法均属于递归式模态分解，均存在模态混叠、端点效应、受采样频率的影响、同时还无法将两个频率相近的分量正确分离等缺点。基于人工智能的故障诊断方法中，目前应用较多的就是人工神经网络，但其缺点是推理过程解释性差，而且当待诊断样本不完备(数据有缺失)时，神经网络不能进行有效的推理工作，无法利用故障的早期特征对轴承进行相应诊断。

由于现有的机械故障诊断手段都在一定程度上存在缺陷，故而采用现有的机械故障诊断方法都难以对轴承故障进行较为准确的识别和诊断。

发明内容

针对现有技术存在的上述不足，本发明的目的在于提出一种基于多特征流形学习与支持向量机的轴承故障诊断方法，本方法能够提高滚动轴承故障诊断的准确性和有效性。

本发明的技术方案是这样实现的：

基于多特征流形学习与支持向量机的轴承故障诊断方法，步骤如下：

1)在四种不同工况下的滚动轴承转动工作时，通过加速度传感器分别采集每种工况下滚动轴承在不同转速工作的振动加速度信号，进行去噪预处理，并添加工况标签，将经过预处理并添加工况标签后的各种工况下的各个振动加速度信号作为训练样本；四种工况分别为正常运转、轴承内圈故障运转、轴承滚动体故障运转、轴承外圈故障运转；

2)对训练样本的四类工况数据分别提取它们的时域特征参数、频域特征参数和时频域特征参数；时域特征参数包括有量纲参数和无量纲参数，其中有量纲参数为均值、均方根值、方根幅值、绝对平均值、方差、最(小)大值、峰—峰值、峭度；无量纲参数为波形指标、峰值指标、脉冲指标、峭度指标、波峰因子、边缘因子、偏斜度、时域信息熵；频域特征参数为总功率谱和、莱斯频率、频率重心、频率方差、谐波指标、均方频率、频域信息熵；时频域特征参数为小波包能量；

3)对四类工况数据多特征提取之后，再进行流行学习，得到从高维数据特征集中提取隐藏其中的低维流形结构；

4)通过加速度传感器采集待测滚动轴承在转动工作时的振动加速度信号，并进行去噪预处理，作为测试样本；

5)提取测试样本的时域特征参数、频域特征参数和时频域特征参数；时域特征参数包括有量纲参数和无量纲参数，其中有量纲参数为均值、均方根值、方根幅值、绝对平均值、方差、最(小)大值、峰—峰值、峭度；无量纲参数为波形指标、峰值指标、脉冲指标、峭度指标、波峰因子、边缘因子、偏斜度、时域信息熵；频域特征参数为总功率谱和、莱斯频率、频率重心、频率方差、谐波指标、均方频率、频域信息熵；时频域特征参数为小波包能量；

6)按步骤3)对测试样本进行流形学习得到低维流行结构；

7)将测试样本的低维流行结构作为测试样本的匹配特征，将每种工况下各个训练样本对应的低维流行结构作为匹配基准，采用支持向量机分类方法对测试样本与训练样本进行匹配，将与测试样本最为匹配的训练样本所属的工况类别判定为测试样本的工况类别，从而得到待测滚动轴承的故障诊断结果。

其中，第2)步和第5)中各特征参数的提取方法为，

提取时域特征参数中的有量纲参数：将四类机械振动数据分割为g个周期数据y_i，i＝1，2，3，…，g，在时域分别按周期提取8个有量纲参数来反映信号的变化，其中：

均值 $\stackrel{\‾}{y}=\frac{1}{T}{\∫}_0^T{y}_{i}dt;$ 均方根值 $y_{rms}=\sqrt{\frac{1}{T}{\∫}_0^T{y}_i^{2}dt};$

方根幅值绝对平均值

方差最(小)大值y_max＝max(y_i)，y_min＝min(y_i)；

峰—峰值为y_ff＝y_max-y_min；峭度

在时域分别按周期提取8个无量纲参数来反映信号的变化，其中：

波形指标 $S_f=\frac{y_{rms}}{\stackrel{\‾}{y}};$ 峰值指标 $C_f=\frac{y_{\max }}{y_{rms}};$

脉冲指标 $I_f=\frac{y_{\max }}{\left|\stackrel{\‾}{y}\right|};$ 峭度指标 $K_v=\frac{\mathrm{\β}}{{y_{rms}}^4};$

波峰因子 $B_f=\frac{\max \left|y_i\right|}{y_{rms}};$ 边缘因子 $M_f=\frac{\max \left|y_i\right|}{{\[\frac{1}{T}{\∫}_0^T\sqrt{\left|y_i\right|}dt\]}^2};$

偏斜度 $P_f=\frac{1}{T}{\∫}_0^T{(y_i-\stackrel{\‾}{y})}^{3}dt/{y_{rms}}^3;$ 时域信息熵 $H_t=-{\mathrm{\Σ}}_{i=1}^M{p}_i{\mathrm{logp}}_i;$

其中为第i个奇异值在全部奇异值之和中所占的比例；δ_i为原信号高维重构相空间的奇异值；

提取频域特征参数：将四类机械振动数据分割为g个周期数据y_i，i＝1，2，3，…，g，在频域提取7个参数来反映信号的变化；

则其FFT变换的公式为：

$s\left(f\right)={\∫}_0^T{y}_i{e}^{-j2\mathrm{\π}ft}dt$

总功率谱和G_t＝∫s(f)df；莱斯频率

频率重心 $f_c={\∫}_0^{\mathrm{\∞}}fs{\left(f\right)}^{2}df/G_t;$ 频率方差 $v_f={\∫}_0^{\mathrm{\∞}}{(f-f_c)}^{2}s\left(f\right)df/G_t;$

谐波指标H＝f_x/f∫_x：；均方频率

频域信息熵 $H_f=-{\mathrm{\Σ}}_{i=1}^g{q}_i{\mathrm{logq}}_i;$

f∫_x：表示x(t)对时间积分再求莱斯频率；频域信息熵H_f定义式中q_i为第i个功率谱值与总功率谱的比值；

时频域特征参数提取过程：利用db4小波包函数对4类机械信号分别进行3层正交小波包分解用以进行时频域特征参数提取，得到这4类信号对应的由全频带均匀划分8个子频带的滤波信号，通过重构每个子频带的节点小波包系数，确保重构信号和原4类机械信号的长度一样，然后计算每个子频带滤波信号的各个采样点的幅值平方和作为其能量，将各个子频带滤波信号的能量与全频带重构信号总能量的比值作为其时频域统计特征，即(E₀/E，E₁/E，...，E₇/E)，其中E₀，E₁，...，E₇分别为8个子频带滤波信号的能量，E为全频带重构信号的总能量；由此针对四类工况信号分别构造了它们的时频域信号，即8个小波包能量信号；

针对获取的四类工况信号，分别构造包括8个有量纲参数和8个无量纲参数的16个时域特征参数、7个频域特征参数、8个时频域特征参数，总共31个特征参数构成滚动轴承信号的特征矩阵。

其中第3)步低维流形学习步骤为：

步骤一：根据特征矩阵的欧氏距离来计算k个近邻点，即计算x_p与其余数据点间的欧式距离，把距离最小的k个数据点作为近邻点；

步骤二：通过取以下函数的最小值计算出样本点的局部重建权值矩阵W；

$\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}\left|\right|x_p-\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}{x}_{pj}|{|}_2^2$

其中：x_pj(j＝1，2，...，k)为x_p的第j个近邻点，w_pj表示由k个近邻点的线性组合重构x_p时第j个近邻点的权值系数，且

(t为热核参数)，x_pj是x_p的近邻点，否则w_pj＝0；

上式可改写为：

$\begin{array}{c}\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}\left|\right|x_p\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}-\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}{x}_{pj}|{|}_2^2=\underset{p=1}{\overset{N}{\Σ}}|\left|\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}(x_p-x_{pj})\right|{|}_2^2\\ =\underset{p=1}{\overset{N}{\Σ}}\left|\right|\[\begin{array}{cccc}{x}_p-x_{p1}& x_p-x_{p2}& \mathrm{...}& x_p-x_{pk}\end{array}\]{\left[\begin{array}{cccc}{w}_{p1}& w_{p2}& \mathrm{...}& w_{pk}\end{array}\right]}^T|{|}_2^2\end{array}$

令w_p＝[w_p1w_p2…w_pk]^T，表示第p个样本点的局部重建权值向量，根据上式可改写为

$\begin{array}{c}\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}{(\[\begin{array}{cccc}{x}_p-x_{p1}& x_p-x_{p2}& \mathrm{...}& x_p-x_{pk}\end{array}\]w_P)}^T(\[\begin{array}{cc}{x}_p-x_{p1}& x_p\end{array}\\ \begin{array}{ccc}-x_{p2}& \mathrm{...}& x_p-x_{pk}\end{array}\]w_P)\\ =\underset{p=1}{\overset{N}{\Σ}}{w_p}^T{\left[\begin{array}{cccc}{x}_p-x_{p1}& x_p-x_{p2}& \mathrm{...}& x_p-x_{pk}\end{array}\right]}^T\[\begin{array}{cc}{x}_p-x_{p1}& x_p\end{array}\\ \begin{array}{ccc}-x_{p2}& \mathrm{...}& x_p-x_{pk}\end{array}\]w_P\end{array}$

令Z_p＝[x_p-x_p1x_p-x_p2...x_p-x_pk]^T[x_p-x_p1x_p-x_p2...x_p-x_pk]

表示第p个样本点的局部协方差，则：

$\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}{w_P}^T{Z}_p{w}_P;$

求解上式是一个带约束条件的最小值问题，可以利用拉格郎日乘数法求解，设

$L\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}\[{w_P}^T{Z}_p{w}_P+\mathrm{\λ}(\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}-1)\]$

对上式的两边分别对w_p求偏导数，可得：

$\frac{\∂L\left(W\right)}{\∂w_P}=2Z_p{w}_P+\mathrm{\λ}e;$

其中e是一个值全为1的k维向量，令可得代入可得λ，最后求得W_p；

步骤三：计算机械数据的低维嵌入流行结果Y，Y应满足如下条件：

$\min \mathrm{\φ}\left(Y\right)=\underset{p=1}{\overset{N}{\Σ}}\left|\right|Y_p-\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}{Y}_{pj}|{|}_2^2$

这里w_pj已知，输出Y应该满足条件：和

将w_pj(p＝1，2，...，N；j＝1，2，...，k)存放在稀疏矩阵W中，W的大小为N×N；如果x_j是x_p的近邻点，且是第t个近邻点(t＝1，2，...，k)，则W_jp＝w_pt；如果x_j不是x_p的近邻点，则W_jp＝0；用W_i表示W矩阵的第P列，I_i表示N×N的第P列，将上式改写为：

$\begin{array}{c}\min \mathrm{\φ}\left(Y\right)=\underset{P=1}{\overset{N}{\Σ}}\left|\right|{\mathrm{YI}}_P-{\mathrm{YW}}_P|{|}_2^2=|\left|Y(I-W)\right|{|}_2^2=tr\left(Y\right(I-W\left){\left(I-W\right)}^T{Y}^T\right)\\ =tr\left({\mathrm{YMY}}^T\right)\end{array}$

其中M＝(I-W)(I-W)^T；

使代价函数最小，存在两个约束条件，采用拉格郎日乘数法，设

L(Y)＝YMY^T+β(YY^T-NI)，

对上式两边求偏导可得：

$\frac{\∂L}{\∂Y}=2{\mathrm{MY}}^T+2{\mathrm{\βY}}^T,$

令可得

MY^T＝-βY^T

令c＝-β，得

MY^T＝cY^T；

为使代价函数最小，取矩阵M的最小d个特征值对应的特征向量作为所求的低维嵌入结果。

其中，步骤7)采用支持向量机分类方法对测试样本与训练样本进行匹配的过程为：

在四种工况训练样本中，将其中第k类训练样本看作正类，k∈{1,2,3,4}，将其它三类训练样本看作负类，通过两类支持向量机分类方法得到第k类的分类决策函数f_k(x)：

$f_k\left(x\right)=\underset{n=1}{\overset{N}{\Σ}}{\mathrm{\α}}_n^k{y}_{n}K(x,x_n)+b_k$

其中，为第k类分类决策函数f_k(x)中第n个训练样本重构信号x_n对应的拉格朗日系数；b_k为第k类分类决策函数f_k(x)的最优超平面位置系数；y_n表示第n个训练样本对应的分类标记，当第n个训练样本属于正类时，y_n＝1；当第n个训练样本属于负类时，y_n＝-1；n∈{1,2,…,N}，N为四种工况的训练样本总数；K(x，x_n)表示分类决策函数f_k(x)的输入量x相对于第n个训练样本重构信号x_n的高斯径向基核函数；

由此得到四种工况中每一种工况所对应的分类决策函数；

最后将测试样本重构信号分别作为四种工况对应的分类决策函数的输入量，计算出测试样本重构信号作为输入量的四个分类决策函数值，以其中最大的分类决策函数值所对应的工况类别判定为测试样本的工况类别，得到待测滚动轴承的故障诊断结果。

相比于现有技术，本发明具有如下有益效果：

1、本发明利用多特征的思想，提取机械数据的时域、频域、时率域特征，充分弥补单特征提取故障特征的不足，再利用非线性降维方法从高维数据特征集中提取隐藏其中的低维流形成分，有效消除冗余，提炼出原始信号内在本质特征，更为方便有效地描述故障特征。

2、本发明采用了支持向量机分类方法对测试样本进行分类识别，支持向量机分类方法中的学习过程可以被看成是一个优化寻找最优解的过程，因此可以采用之前设计好的有效方法去寻找和发现目标函数的全局最小值，提高故障识别的有效性。

3、与现有技术比较，本发明的滚动轴承故障诊断方法能够提高滚动轴承故障诊断的准确性和有效性，为解决滚动轴承故障诊断问题提供一种新的新的思路新的方法，可广泛应用于化工、冶金、电力、航空等领域的复杂系统中。

附图说明

图1-本发明基于多特征流形学习与支持向量机的轴承故障诊断方法流程图。

图2-滚动轴承内圈故障运转的原始振动加速度信号时域分布示例图(时域单位为ms)。

图3-滚动轴承外圈故障运转的原始振动加速度信号时域分布示例图(时域单位为ms)。

图4-LLE维数约减过程示意图。

图5-SVM规范化最优分类超平面的位置关系示意图。

具体实施方式

下面结合附图对本发明作进一步详细说明。

本发明基于多特征流形学习与支持向量机的滚动轴承故障诊断方法的思路是，首先获取机械滚动轴承数据的时域特征、频域特征、时频域特征，得到特征矩阵；之后采用流形学习算法从高维数据特征集中提取隐藏其中的低维流形成分，有效消除冗余，提炼出原始信号内在本质特征；最后采用支持向量机分类方法对测试样本进行分类识别，从而确定滚动轴承故障工况类别，实现对滚动轴承故障类别的诊断，以提高滚动轴承故障诊断的准确性和有效性。

针对时域特征随着故障程度的不断加深，其值会呈现出一种非正常的大幅度波动的特点，频域特征能够直接从信号中找到感兴趣频率成分，但是它们都是针对平稳理论提出来的，时频特征可以刻画非平稳信号统计特征随时间变化的情况，本发明利用这三种特征的特点，提出多特征流形学习，利用流形学习从高维数据特征集中提取隐藏其中的低维流形成分，有效消除冗余，提炼出原始信号内在本质特征，更为方便有效地描述故障特征。

支持向量机(SupportVectorMachine，简称SVM)分类器的核心思想是通过某种非线性映射(核函数)将输入向量映射到一个高维特征空间，并构造最优分类超平面，从而实现分类识别。在解决小样本、非线性及高维模式识别中具有独特优势，又能很好地限制过学习，特别适合于小样本集的数据处理，因此能够在故障诊断和故障预测方面得到应用。

鉴于基于多特征流形学习与支持向量机具备的上述优点，本发明正是将多特征流形学习和支持向量机具备的上述优点加以整合，利用多特征流形学习与支持向量机进行滚动轴承故障工况的分类，实现对滚动轴承故障的识别和诊断，其具体操作流程如图1所示，包括如下步骤：

1)在四种不同工况下的滚动轴承转动工作时，通过加速度传感器分别采集每种工况下滚动轴承在不同转速工作的振动加速度信号，进行去噪预处理，并添加工况标签，将经过预处理并添加工况标签后的各种工况下的各个振动加速度信号数据作为训练样本；所述四种工况分别为正常运转、轴承内圈故障运转、轴承滚动体故障运转、轴承外圈故障运转；

滚动轴承在四种不同工况下转动工作的振动加速度信号相互之间存在一定的差异，例如，图2和图3分别示出了滚动轴承在内圈故障运转、外圈故障运转工况下的原始振动加速度信号时域分布图(时域单位为ms)，其信号差异较为明显。因此可以基于滚动轴承在不同工况下的振动加速度信号特征数据，对其故障情况进行识别。

2)对训练样本的四类工况数据分别提取它们的时域特征：有量纲参数(均值、均方根值、方根幅值、绝对平均值、方差、最(小)大值、峰—峰值、峭度)；无量纲参数(波形指标、峰值指标、脉冲指标、峭度指标、波峰因子、边缘因子、偏斜度、时域信息熵)；频域特征(总功率谱和、莱斯频率、频率重心、频率方差、谐波指标、均方频率、频域信息熵)；时频域特征(小波包能量)。

提取时域参数中的有量纲参数：将四类机械振动数据分割为g个周期数据y_i，i＝1，2，3，…，g，在时域分别按周期提取8个有量纲参数来反映信号的变化，其中：

均值 $\stackrel{\‾}{y}=\frac{1}{T}{\∫}_0^T{y}_{i}dt;$ 均方根值 $y_{rms}=\sqrt{\frac{1}{T}{\∫}_0^T{y}_i^{2}dt};$

方根幅值 $y_{\mathrm{\τ}}={\[\frac{1}{T}{\∫}_0^T\sqrt{\left|y_i\right|}dt\]}^2;$ 绝对平均值 $\left|\stackrel{\‾}{y}\right|=\frac{1}{T}{\∫}_0^T\left|y_i\right|dt;$

方差最(小)大值y_max＝max(y_i)，y_min＝min(y_i)；

峰—峰值为y_ff＝y_max-y_min；峭度

在时域分别按周期提取8个无量纲参数来反映信号的变化，其中：

波形指标为 $S_f=\frac{y_{rms}}{\stackrel{\‾}{y}};$ 峰值指标 $C_f=\frac{y_{\max }}{y_{rms}};$

脉冲指标为 $I_f=\frac{y_{\max }}{\left|\stackrel{\‾}{y}\right|};$ 峭度指标 $K_v=\frac{\mathrm{\β}}{{y_{rms}}^4};$

波峰因子为 $B_f=\frac{\max \left|y_i\right|}{y_{rms}};$ 边缘因子 $M_f=\frac{\max \left|y_i\right|}{{\[\frac{1}{T}{\∫}_0^T\sqrt{\left|y_i\right|}dt\]}^2};$

偏斜度为 $P_f=\frac{1}{T}{\∫}_0^T{(y_i-\stackrel{\‾}{y})}^{3}dt/{y_{rms}}^3;$ 时域信息熵为 $H_t=-{\mathrm{\Σ}}_{i=1}^M{p}_i{\mathrm{logp}}_i;$

其中为第i个奇异值在全部奇异值之和中所占的比例(δ_i为原信号高维重构相空间的奇异值)；

提取频域参数：将四类机械振动数据分割为g个周期数据y_i，i＝1，2，3，…，g，在频域提取7个参数来反映信号的变化。

则其FFT变换的公式为：

$s\left(f\right)={\∫}_0^T{y}_i{e}^{-j2\mathrm{\π}ft}dt$

总功率谱和为G_t＝∫s(f)df；莱斯频率为

频率重心为 $f_c={\∫}_0^{\mathrm{\∞}}fs{\left(f\right)}^{2}df/G_t;$ 频率方差 $v_f={\∫}_0^{\mathrm{\∞}}{(f-f_c)}^{2}s\left(f\right)df/G_t;$

谐波指标H＝f_x/f∫_x：；均方频率为

频域信息熵为 $H_f=-{\mathrm{\Σ}}_{i=1}^g{q}_i{\mathrm{logq}}_i;$

f∫_x：表示x(t)对时间积分再求莱斯频率；频域信息熵H_f定义式中：q_i为第i个功率谱值与总功率谱的比值。

利用db4小波包函数对4类机械信号分别进行3层正交小波包分解用以进行时频域特征提取，得到这4类信号对应的由全频带均匀划分8个子频带的滤波信号，通过重构每个子频带的节点小波包系数，确保重构信号和原4类机械信号的长度一样，然后计算每个子频带滤波信号的各个采样点的幅值平方和作为其能量，将各个子频带滤波信号的能量与全频带重构信号总能量的比值作为其时频域统计特征，即(E₀/E，E₁/E，...，E₇/E)，其中E₀，E₁，...，E₇分别为8个子频带滤波信号的能量，E为全频带重构信号的总能量；由此针对这四类的机械信号分别构造了它们的时频域信号，即8个小波包能量信号。

针对获取的以上四类机械信号，分别构造16个时域特征(包括8个有量纲参数和8个无量纲参数)、7个频域统计特征、8个时频域参数，构成以上滚动轴承信号的特征矩阵。

3)对四类工况数据的多特征提取之后，采用LLE(LocallyLinearEmbedding)进行流形学习，得到从高维数据特征集中提取隐藏其中的低维流形成分，有效消除冗余，提炼出原始信号内在本质特征，更为方便有效地描述故障特征。LLE维数约减过程如图4所示。

步骤一：根据样本特征矩阵的欧氏距离来计算k个近邻点；即计算x_p与其余数据点间的欧式距离，把距离最小的k个数据点作为近邻点。

步骤二：通过取以下函数的最小值计算出样本点的局部重建权值矩阵W；

$\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}\left|\right|x_p-\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}{x}_{pj}|{|}_2^2$

其中：x_pj(j＝1，2，...，k)为x_p的第j个近邻点，w_pj表示由k个近邻点的线性组合重构x_p时第j个近邻点的权值系数，且

(t为热核参数)，x_pj是x_p的近邻点，否则w_pj＝0；

上式可改写为，

$\begin{array}{c}\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}\left|\right|x_p\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}-\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}{x}_{pj}|{|}_2^2=\underset{p=1}{\overset{N}{\Σ}}|\left|\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}(x_p-x_{pj})\right|{|}_2^2\\ =\underset{p=1}{\overset{N}{\Σ}}\left|\right|\[\begin{array}{cccc}{x}_p-x_{p1}& x_p-x_{p2}& \mathrm{...}& x_p-x_{pk}\end{array}\]{\left[\begin{array}{cccc}{w}_{p1}& w_{p2}& \mathrm{...}& w_{pk}\end{array}\right]}^T|{|}_2^2\end{array}$

令w_p＝[w_p1w_p2…w_pk]^T，表示第p个样本点的局部重建权值向量，根据上式可改写为

$\begin{array}{c}\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}{(\[\begin{array}{cccc}{x}_p-x_{p1}& x_p-x_{p2}& \mathrm{...}& x_p-x_{pk}\end{array}\]w_P)}^T(\[\begin{array}{cc}{x}_p-x_{p1}& x_p\end{array}\\ \begin{array}{ccc}-x_{p2}& \mathrm{...}& x_p-x_{pk}\end{array}\]w_P)\\ =\underset{p=1}{\overset{N}{\Σ}}{w_p}^T{\left[\begin{array}{cccc}{x}_p-x_{p1}& x_p-x_{p2}& \mathrm{...}& x_p-x_{pk}\end{array}\right]}^T\[\begin{array}{cc}{x}_p-x_{p1}& x_p\end{array}\\ \begin{array}{ccc}-x_{p2}& \mathrm{...}& x_p-x_{pk}\end{array}\]w_P\end{array}$

令Z_p＝[x_p-x_p1x_p-x_p2...x_p-x_pk]^T[x_p-x_p1x_p-x_p2...x_p-x_pk]

表示第p个样本点的局部协方差，则：

$\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}{w_P}^T{Z}_p{w}_P$

求解上式是一个带约束条件的最小值问题，可以利用拉格郎日乘数法求解，设

$L\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}\[{w_P}^T{Z}_p{w}_P+\mathrm{\λ}(\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}-1)\]$

对上式的两边分别对w_p求偏导数，可得：

$\frac{\∂L\left(W\right)}{\∂w_P}=2Z_p{w}_P+\mathrm{\λ}e$

其中e是一个值全为1的k维向量，令可得代入可得λ，最后求得w_P；

第三步：计算机械数据的低维嵌入流行结果Y，Y应满足如下条件：

$\min \mathrm{\φ}\left(Y\right)=\underset{p=1}{\overset{N}{\Σ}}\left|\right|Y_p-\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}{Y}_{pj}|{|}_2^2$

这里w_pj已知，输出Y应该满足条件：和

可以将w_pj(p＝1，2，...，N；j＝1，2，...，k)存放在稀疏矩阵W中，W的大小为N×N。如果x_j是x_p的近邻点，且是第t个近邻点(t＝1，2，...，k)，则W_jp＝w_pt；如果x_j不是x_p的近邻点，则W_jp＝0。用W_i表示W矩阵的第P列，I_i表示N×N的第P列，可以将上式改写为：

$\begin{array}{c}\min \mathrm{\φ}\left(Y\right)=\underset{P=1}{\overset{N}{\Σ}}\left|\right|{\mathrm{YI}}_P-{\mathrm{YW}}_P|{|}_2^2=|\left|Y(I-W)\right|{|}_2^2=tr\left(Y\right(I-W\left){\left(I-W\right)}^T{Y}^T\right)\\ =tr\left({\mathrm{YMY}}^T\right)\end{array}$

其中M＝(I-W)(I-W)^T；

使代价函数最小，存在两个约束条件，同样可以采用拉格郎日乘数法，设

L(Y)＝YMYT+β(YYT-NI)，

对上式两边求偏导可得：

$\frac{\∂L}{\∂Y}=2{\mathrm{MY}}^T+2{\mathrm{\βY}}^T;$

令可得

MY^T＝-βY^T；

令c＝-β，得

MY^T＝cY^T；

为使代价函数最小，取矩阵M的最小d个特征值对应的特征向量作为所求的低维嵌入结果。

4)通过加速度传感器采集待测滚动轴承在转动工作时的振动加速度信号数据，并进行去噪预处理，作为测试样本。

5)利用同样的方法提取测试样本的时域特征：有量纲参数(均值、均方根值、方根幅值、绝对平均值、方差、最(小)大值、峰—峰值、峭度)；无量纲参数(波形指标、峰值指标、脉冲指标、峭度指标、波峰因子、边缘因子、偏斜度、时域信息熵)，频域特征(总功率谱和、莱斯频率、频率重心、频率方差、谐波指标、均方频率、频域信息熵)；时频域特征(小波包能量)。

6)对测试样本进行流形学习。

7)将测试样本的低维流形结构作为测试样本的匹配特征，将每种工况下的各个训练样本对应的的低维流形作为匹配基准，采用支持向量机分类方法对测试样本与训练样本进行匹配，将与测试样本最为匹配的训练样本所属的工况类别判定为测试样本的工况类别，从而得到待测滚动轴承的故障诊断结果。

支持向量机(SupportVectorMachines，简称SVM)于1963年在AT&TBell实验室由Vapnik等提出，它是以统计学中的VC维理论和结构风险最小原理为理论基础，根据有限的样本信息在模型的复杂度(即特定训练样本的学习精度)和学习能力(即正确识别任意样本的能力)之间寻求最佳折衷，以获得最好的推广能力。SVM将向量映射到一个更高维的空间里，在高维空间中建立一个最大分隔超平面，并在能够分开数据超平面的两边建立两个互相平行的超平面，分隔超平面使两个平行超平面的距离最大化，其距离越大，分类结果的误差越小。

图5为二维两类情况下的规范化最优超平面的位置关系示意图，H为分隔超平面，H1、H2为两互相平行的超平面，H1、H2为分类间隔d＝2/||w||。为保证数据的线性化，需把数据映射到核函数空间；同时，为能有效地把两类分开，应保证两类被正确分开使分类间隔最大，即有目标函数为：

$\underset{w,b}{min}\frac{1}{2}\left|\right|w|{|}^2+C\left({\mathrm{\Σ}}_{i=1}^N{\mathrm{\ϵ}}_i\right)=\frac{1}{2}(w\·w)+C\left({\mathrm{\Σ}}_{i=1}^N{\mathrm{\ϵ}}_I\right)---\left(1\right)$

对超平面H₁、H₂有：

上式可等价为：

其中：表示w与的内内为积，表示把x_i映射到核函数空间，b和C为常数；ε_i＞0为松弛变量，表示训练样本错分程度，其值越大表示错分样本越多。对(1)和(2)运用Lagrange乘子法，得到：

其中：ε_i＞0，β_i＞0为拉格朗日系数，L(w，b，α)为拉格朗日函数。

式(3)中对w、ε_i和b的偏导数为零，得到：

将(4)代入(3)，最优超平面的解等价于以下对偶问题的解。

$\begin{array}{c}\max Q\left(\mathrm{\α}\right)=\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\α}}_i-\frac{1}{2}\underset{i,j=1}{\overset{N}{\Σ}}{\mathrm{\α}}_i{\mathrm{\α}}_j{y}_i{y}_{j}K(x_i\·x_j)\\ s.t.{\Σ}_{i=1}^N{y}_i{\mathrm{\α}}_i=0,C>{\mathrm{\α}}_i>0\end{array}---\left(5\right)$

其中：

运用Lagrange乘子法，得到解为：

由式(6)得到最优分类面的分类规则函数为：

本发明选用高斯径向基(RBF)核函数：

$K(x\·x_j)=\exp (-\frac{\left|\right|x-x_j|{|}^2}{{\mathrm{\σ}}^2}))$

其中：σ为RBF核函数的参数。

SVM算法是以统计学为基础寻求数据间的最优分类面，通过把非线性数据映射到核函数空间，使其线性化，进而简化了计算复杂度，具有较好的分类效果。

步骤7)中，可以运用的支持向量机分类方法的具体分类方式有很多，例如一对一分类(OVO-SVM)、一对多分类(也称为一对余分类，OVR-SVM)、有向无环图分类(DAG-SVMS)、决策树分类、纠错输出编码分类等。但考虑到本发明方法中仅涉及正常运转、轴承内圈故障运转、轴承滚动体故障运转、轴承外圈故障运转这四种工况的故障工况分类识别，识别类别并不多，综合考虑识别效率和准确性的因素，采用一对多的SVM分类方法较为适合，因为采用一对多分类方法需要建立和进行识别运算的SVM分类器决策函数只有四个(每种故障工况类别对应一个)，并且滚动轴承在该四种不同工况下特征提取后经过流行学习，其本质特征的区别已经足以识别，能够较好的保证识别准确性。

在本发明中，步骤7)优选采用一对多支持向量机分类方法对测试样本与训练样本进行匹配的具体方式为：

在四种工况的训练样本中，针对其中第k类训练样本看作正类，k∈{1,2,3,4}，将其它3类训练样本看作负类，通过两类支持向量机分类方法得到第k类的分类决策函数：

$f_k\left(x\right)=\underset{n=1}{\overset{N}{\Σ}}{\mathrm{\α}}_n^k{y}_{n}K(x,x_n)+b_k$

其中，为第k类分类决策函数f_k(x)中第n个的训练样本重构信号x_n对应的拉格朗日系数；b_k为第k类分类决策函数f_k(x)的最优超平面位置系数；y_n表示第n个训练样本对应的分类标记，当第n个训练样本属于正类时，y_n＝1。

当第n个训练样本属于负类时y_n＝-1；n∈{1,2,…,N}，N为四种工况的训练样本的总数；K(x，x_n)表示分类决策函数f_k(x)的输入量x相对于第n个的训练样本重构信号x_n的高斯径向基核函数；

由此得到四种工况中每一种工况所对应的分类决策函数。

将测试样本重构信号分别作为四种工况对应的分类决策函数的输入量，计算出测试样本重构信号作为输入量的四个分类决策函数值，以其中最大的分类决策函数值所对应的工况类别判定为测试样本的工况类别，得到待测滚动轴承的故障诊断结果。

通过实验数据验证，采用本发明的基于多特征流行学习与支持向量机的滚动轴承故障诊断方法按上述流程进行故障诊断，在200个训练样本(每种工况50个训练样本)的条件下，随机抽样进行100次滚动轴承故障诊断识别，其识别准确率达到95.3％，完全能够满足实际应用需求。

综上所述，本发明基于多特征流形学习与支持向量机的滚动轴承故障诊断方法，利用多特征的思想，提取机械数据的时域、频域、时率域特征，充分弥补单特征提取故障特征的不足，再利用非线性降维方法从高维数据特征集中提取隐藏其中的低维流形成分，有效消除冗余，提炼出原始信号内在本质特征，更为方便有效地描述故障特征。本发明采用了支持向量机分类方法对测试样本进行分类识别，支持向量机分类方法中的学习过程可以被看成是一个优化寻找最优解的过程，因此可以采用之前设计好的有效方法去寻找和发现目标函数的全局最小值，提高故障识别的有效性。与现有技术比较，本发明的滚动轴承故障诊断方法能够提高滚动轴承故障诊断的准确性和有效性，为解决滚动轴承故障诊断问题提供一种新的思路，可广泛应用于化工、冶金、电力、航空等领域的复杂系统中。

最后需要说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管申请人参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，均应涵盖在本发明的权利要求范围当中。

Claims (4)

1.基于多特征流形学习与支持向量机的轴承故障诊断方法，其特征在于：步骤如下：

1)在四种不同工况下的滚动轴承转动工作时，通过加速度传感器分别采集每种工况下滚动轴承在不同转速工作的振动加速度信号，进行去噪预处理，并添加工况标签，将经过预处理并添加工况标签后的各种工况下的各个振动加速度信号作为训练样本；四种工况分别为正常运转、轴承内圈故障运转、轴承滚动体故障运转、轴承外圈故障运转；

2)对训练样本的四类工况数据分别提取它们的时域特征参数、频域特征参数和时频域特征参数；时域特征参数包括有量纲参数和无量纲参数，其中有量纲参数为均值、均方根值、方根幅值、绝对平均值、方差、最(小)大值、峰—峰值、峭度；无量纲参数为波形指标、峰值指标、脉冲指标、峭度指标、波峰因子、边缘因子、偏斜度、时域信息熵；频域特征参数为总功率谱和、莱斯频率、频率重心、频率方差、谐波指标、均方频率、频域信息熵；时频域特征参数为小波包能量；

3)对四类工况数据多特征提取之后，再进行流行学习，得到从高维数据特征集中提取隐藏其中的低维流形结构；

4)通过加速度传感器采集待测滚动轴承在转动工作时的振动加速度信号，并进行去噪预处理，作为测试样本；

5)提取测试样本的时域特征参数、频域特征参数和时频域特征参数；时域特征参数包括有量纲参数和无量纲参数，其中有量纲参数为均值、均方根值、方根幅值、绝对平均值、方差、最(小)大值、峰—峰值、峭度；无量纲参数为波形指标、峰值指标、脉冲指标、峭度指标、波峰因子、边缘因子、偏斜度、时域信息熵；频域特征参数为总功率谱和、莱斯频率、频率重心、频率方差、谐波指标、均方频率、频域信息熵；时频域特征参数为小波包能量；

6)按步骤3)对测试样本进行流形学习得到低维流行结构；

7)将测试样本的低维流行结构作为测试样本的匹配特征，将每种工况下各个训练样本对应的低维流行结构作为匹配基准，采用支持向量机分类方法对测试样本与训练样本进行匹配，将与测试样本最为匹配的训练样本所属的工况类别判定为测试样本的工况类别，从而得到待测滚动轴承的故障诊断结果。

2.根据权利要求1所述的基于多特征流形学习与支持向量机的轴承故障诊断方法，其特征在于，第2)步和第5)中各特征参数的提取方法为，

提取时域特征参数中的有量纲参数：将四类机械振动数据分割为g个周期数据y_i，i＝1，2，3，...，g，在时域分别按周期提取8个有量纲参数来反映信号的变化，其中：

均值 $\stackrel{\‾}{y}=\frac{1}{T}{\∫}_0^T{y}_{i}dt;$ 均方根值 $y_{rms}=\sqrt{\frac{1}{T}{\∫}_0^T{y}_i^{2}dt};$

方根幅值 $y_{\mathrm{\τ}}={\[\frac{1}{T}{\∫}_0^T\sqrt{\left|y_i\right|}dt\]}^2;$ 绝对平均值 $\[\stackrel{\‾}{y}|=\frac{1}{T}{\∫}_0^{\mathrm{\τ}}\left|y_i\right|dt;$

方差 $D_y=\frac{1}{T}{\∫}_o^T{(y_i-\stackrel{\‾}{y})}^{2}dt;$ 最(小)大值y_max＝max(y_i)，y_min＝min(y_i)；

峰—峰值为y_ff＝y_max-y_min；峭度

在时域分别按周期提取8个无量纲参数来反映信号的变化，其中：

波形指标 $S_f=\frac{y_{rms}}{\stackrel{\‾}{y}};$ 峰值指标 $C_f=\frac{y_{max}}{y_{rms}};$

脉冲指标 $I_f=\frac{y_{max}}{\left|\stackrel{\‾}{y}\right|};$ 峭度指标 $K_v=\frac{\mathrm{\β}}{{y_{rms}}^4};$

波峰因子 $B_f=\frac{max\left|y_i\right|}{y_{rms}};$ 边缘因子 $M_f=\frac{max\left|y_i\right|}{{\[\frac{1}{T}{\∫}_0^T\sqrt{\left|y_i\right|}dt\]}^2};$

偏斜度 $P_f=\frac{1}{T}{\∫}_0^T{(y_i-\stackrel{\‾}{y})}^{3}dt/{y_{rms}}^3;$ 时域信息熵 $H_t=-{\mathrm{\Σ}}_{i=1}^N{p}_i{\mathrm{logp}}_i;$

其中为第i个奇异值在全部奇异值之和中所占的比例；δ_i为原信号高维重构相空间的奇异值；

提取频域特征参数：将四类机械振动数据分割为g个周期数据y_i，i＝1，2，3，...，g，在频域提取7个参数来反映信号的变化；

则其FFT变换的公式为：

$s\left(f\right)={\∫}_0^T{y}_i{e}^{-j2\mathrm{\π}ft}dt$

总功率谱和G_t＝∫s(f)df；莱斯频率

频率重心 $f_c={\∫}_0^{\mathrm{\∞}}fs{\left(f\right)}^{2}df/G_t;$ 频率方差 $v_f={\∫}_0^{\mathrm{\∞}}{(f-f_c)}^{2}s\left(f\right)df/G_t;$

谐波指标H＝f_x/f∫_x：；均方频率

频域信息熵 $H_f=-{\mathrm{\Σ}}_{i=1}^g{q}_i{\mathrm{logq}}_i;$

f∫_x：表示x(t)对时间积分再求莱斯频率；频域信息熵H_f定义式中q_i为第i个功率谱值与总功率谱的比值；

时频域特征参数提取过程：利用db4小波包函数对4类机械信号分别进行3层正交小波包分解用以进行时频域特征参数提取，得到这4类信号对应的由全频带均匀划分8个子频带的滤波信号，通过重构每个子频带的节点小波包系数，确保重构信号和原4类机械信号的长度一样，然后计算每个子频带滤波信号的各个采样点的幅值平方和作为其能量，将各个子频带滤波信号的能量与全频带重构信号总能量的比值作为其时频域统计特征，即(E₀/E，E₁/E，...，E₇/E)，其中E₀，E₁，...，E₇分别为8个子频带滤波信号的能量，E为全频带重构信号的总能量；由此针对四类工况信号分别构造了它们的时频域信号，即8个小波包能量信号；

针对获取的四类工况信号，分别构造包括8个有量纲参数和8个无量纲参数的16个时域特征参数、7个频域特征参数、8个时频域特征参数，总共31个特征参数构成滚动轴承信号的特征矩阵。

3.根据权利要求1所述的基于多特征流形学习与支持向量机的轴承故障诊断方法，其特征在于，第3)步低维流形学习步骤为：

步骤一：根据特征矩阵的欧氏距离来计算k个近邻点，即计算x_p与其余数据点间的欧式距离，把距离最小的k个数据点作为近邻点；

步骤二：通过取以下函数的最小值计算出样本点的局部重建权值矩阵W；

$\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}\left|\right|x_p-\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}{x}_{pj}|{|}_2^2$

其中：x_pj(j＝1，2，...，k)为x_p的第j个近邻点，w_pj表示由k个近邻点的线性组合重构x_p时第j个近邻点的权值系数，且

(t为热核参数)，x_pj是x_p的近邻点，否则w_pj＝0；

上式可改写为：

$\begin{array}{c}\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}\left|\right|x_p\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}-\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}{x}_{pj}|{|}_2^2=\underset{p=1}{\overset{N}{\Σ}}|\left|\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}(x_p-x_{pj})\right|{|}_2^2\\ =\underset{P=1}{\overset{N}{\Σ}}\left|\right|\left[\begin{array}{cccc}{x}_p-x_{p1}& x_p-x_{p2}& ...& x_p-x_{pk}\end{array}\right]{\left[\begin{array}{cccc}{w}_{p1}& w_{p2}& ...& w_{pk}\end{array}\right]}^T|{|}_2^2\end{array}$

令w_p＝[w_p1w_p2...w_pk]^T，表示第p个样本点的局部重建权值向量，根据上式可改写为

$\begin{array}{c}\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}{\left(\left[\begin{array}{cccc}{x}_p-x_{p1}& x_p-x_{p2}& ...& x_p-x_{pk}\end{array}\right]w_P\right)}^T(\[\begin{array}{cc}{x}_p-x_{p1}& x_p\end{array}\\ \begin{array}{ccc}-x_{p2}& ...& x_p-x_{pk}\]w_P)\end{array}\\ =\underset{p=1}{\overset{N}{\Σ}}{w_p}^T{\left[\begin{array}{cccc}{x}_p-x_{p1}& x_p-x_{p2}& ...& x_p-x_{pk}\end{array}\right]}^T\[\begin{array}{cc}{x}_p-x_{p1}& x_p\end{array}\\ \begin{array}{ccc}-x_{p2}& ...& x_p-x_{pk}\]w_P)\end{array}\end{array}$

令Z_p＝[x_p-x_p1x_p-x_p2...x_p-x_pk]^T[x_p-x_p1x_p-x_p2...x_p-x_pk]

表示第p个样本点的局部协方差，则：

$\min \mathrm{\ϵ}\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}{w_P}^T{Z}_p{w}_P;$

求解上式是一个带约束条件的最小值问题，可以利用拉格郎日乘数法求解，设

$L\left(W\right)=\underset{p=1}{\overset{N}{\Σ}}\[{w_P}^T{Z}_p{w}_P+\mathrm{\λ}(\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}-1)\]$

对上式的两边分别对w_p求偏导数，可得：

$\frac{\∂L\left(W\right)}{\∂w_P}=2Z_p{w}_P+\mathrm{\λ}e;$

其中e是一个值全为1的k维向量，令 $\frac{\∂L\left(W\right)}{\∂w_P}=0,$ 可得 $Z_p{w}_P=-\frac{\mathrm{\λ}}{2}e,w_P=-\frac{\mathrm{\λ}}{2}{Z_p}^{-1}e,$ 代入可得λ，最后求得w_p；

步骤三：计算机械数据的低维嵌入流行结果Y，Y应满足如下条件：

$\min \mathrm{\φ}\left(Y\right)=\underset{P=1}{\overset{N}{\Σ}}\left|\right|Y_p-\underset{j=1}{\overset{k}{\Σ}}{w}_{pj}{Y}_{pj}|{|}_2^2$

这里w_pj已知，输出Y应该满足条件：和

将w_pj(p＝1，2，...，N；j＝1，2，...，k)存放在稀疏矩阵W中，W的大小为N×N；如果x_j是x_p的近邻点，且是第t个近邻点(t＝1，2，...，k)，则W_jp＝w_pt；如果x_j不是x_p的近邻点，则W_jp＝0；

用W_i表示W矩阵的第P列，I_i表示N×N的第P列，将上式改写为：

$\begin{array}{c}\min \mathrm{\φ}\left(Y\right)=\underset{P=1}{\overset{N}{\Σ}}\left|\right|{\mathrm{YI}}_P-{\mathrm{YW}}_P|{|}_2^2=|\left|Y(I-W)\right|{|}_2^2=tr\left(Y\right(I-W\left){\left(I-W\right)}^T{Y}^T\right)\\ =tr\left({\mathrm{YMY}}^T\right)\end{array}$

其中M＝(I-W)(I-W)^T；

使代价函数最小，存在两个约束条件，采用拉格郎日乘数法，设

L(Y)＝YMY^T+β(YY^T-NI)，

对上式两边求偏导可得：

$\frac{\∂L}{\∂Y}=2{\mathrm{MY}}^T+2{\mathrm{\βY}}^T,$

令可得

MY^T＝-βY^T

令c＝-β，得

MY^T＝cY^T；

为使代价函数最小，取矩阵M的最小d个特征值对应的特征向量作为所求的低维嵌入结果。

4.根据权利要求1所述的基于多特征流形学习与支持向量机的轴承故障诊断方法，其特征在于：步骤7)采用支持向量机分类方法对测试样本与训练样本进行匹配的过程为：

在四种工况训练样本中，将其中第k类训练样本看作正类，k∈{1,2,3,4}，将其它三类训练样本看作负类，通过两类支持向量机分类方法得到第k类的分类决策函数f_k(x)： $f_k\left(x\right)=\underset{n=1}{\overset{N}{\Σ}}{\mathrm{\α}}_n^k{y}_{n}K(x,x_n)+b_k$

其中，为第k类分类决策函数f_k(x)中第n个训练样本重构信号x_n对应的拉格朗日系数；b_k为第k类分类决策函数f_k(x)的最优超平面位置系数；y_n表示第n个训练样本对应的分类标记，当第n个训练样本属于正类时，y_n＝1；当第n个训练样本属于负类时，y_n＝-1；n∈{1,2,…,N}，N为四种工况的训练样本总数；K(x，x_n)表示分类决策函数f_k(x)的输入量x相对于第n个训练样本重构信号x_n的高斯径向基核函数；

由此得到四种工况中每一种工况所对应的分类决策函数；

最后将测试样本重构信号分别作为四种工况对应的分类决策函数的输入量，计算出测试样本重构信号作为输入量的四个分类决策函数值，以其中最大的分类决策函数值所对应的工况类别判定为测试样本的工况类别，得到待测滚动轴承的故障诊断结果。