回转支承作为回转机械的重要组成部件，具有承受轴向载荷、径向载荷以及倾覆力矩的功能. 其结构尺寸较大，在实际使用过程中安装、拆卸不方便，并且对其进行维修所需的停机时间长、成本高. 因此，对回转支承进行运行状态在线监测与故障诊断尤为重要^[1-2]，不仅可以避免因设备故障而造成巨大损失，同时可以充分利用回转支承的使用价值.

针对回转支承等机械设备振动信号故障信息微弱、特征提取困难的问题，Jaouher等^[3-4]提出基于经验模态分解获取能量熵的特征提取方法，并结合人工神经网络实现轴承的故障识别. 魏永合等^[5]提取16个轴承振动信号统计特征并采用遗传算法降维，以实现对轴承退化状态特征的提取. 陆超等^[6-7]提取振动信号8个时域和8个时频域指标组成特征向量，结合概率主成分分析降维方法实现对回转支承振动信号多领域特征的提取. 一般时域、频域指标分析法未考虑振动信号的非线性问题，且在利用经验模态分解处理非线性非平稳信号时存在过包络、欠包络等问题^[8]，直接影响后续特征的提取效果.

研究表明，机械振动信号作为典型的非线性、非平稳信号，在一定的尺度范围内具有分形特征^[9]，这一特征能精确刻画信号的内在几何结构特点. Lin等^[10]利用多分形去趋势波动分析（multi-fractal detrended fluctuation analysis，MF-DFA）方法量化轴承数据的多分形性，并提取5个多分形特征用于表征轴承故障状态. Xiong等^[11]利用核主成分分析方法将5个MF-DFA特征与5个稳态分布特征进行特征融合，以提取滚动轴承故障特征. 在MF-DFA方法中，须经复杂的Legendre变换才能得到多分形特征值，计算量较大；在Wavelet leader方法中，采用Chhabra算法替代复杂的Legendre变换，计算简便. 因此，本研究利用Wavelet leader方法提取回转支承振动信号的高维多分形特征. 高维特征包含大量重复冗余信息，导致信息利用率较低且计算量较大. 等距映射（isometric mapping，ISOMAP）算法作为新型非线性数据压缩算法，在去除冗余信息的同时能较好地保持数据空间结构不变. 因此，为了提取回转支承振动信号的有效特征信息，本研究提出基于Wavelet leader和融合参数优化算法的ISOMAP算法的多分形自适应特征提取方法.

本研究通过Wavelet leader方法计算正常状态、螺栓破坏和外圈破坏3种状态下的多分形特征，挖掘数据几何结构信息，由此构造出高维多分形特征矩阵，并通过优化的ISOMAP方法对高维矩阵进行有效特征的自适应筛选. 为了定量分析特征提取效果，利用经遗传算法优化的最小二乘支持向量机（genetic algorithm-least squares support vector machine，GA-LSSVM）方法对筛选后的特征进行建模分析，并通过某型号回转支承的全寿命加速实验进行验证.

设所需处理的信号为 $X(t)$，假设母小波 ${\psi _0}(t)$满足： $\forall k = 0,1,\cdots,{N_\psi } - 1,$ $\int_{\bf{R}} {{t^k}{\psi _0}(t){\rm{d}}t \equiv 0} $， $\int_{\bf{R}} {{t^{{N_\psi }}}{\psi _0}(t){\rm{d}}t \ne } $0. 其中， ${N_\psi }$为 ${\psi _0}(t)$的消失矩（ ${N_\psi }$为正整数，且 ${N_\psi } \geqslant 1$），t为时间序列，k为整数^[12-13]. 对 ${\psi _0}(t)$进行伸缩和平移产生一系列函数簇 $\left\{ {{\psi _{j,k}}(t) = {2^{ - j/2}}{\psi _{0}} }\times\right.$ $\left.({2^{ - j}}t - k),\; j \in {\bf{Z}},k \in {\bf{Z}}\right\}$以构成Lebesgue空间 ${L^2}({\bf{R}})$上的正交基. $X(t)$的离散小波变换系数为

(1) $ {d_X}(j,k) = \int_{\bf{R}} {X(t){2^{ - j}}} {\psi _0}({2^{ - j}}t - k){\rm{d}}t. $

定义二元区间 ${\rm{\lambda }} = {{\rm{\lambda }}_{j,k}} = [k{2^j},(k + 1){2^j}]$，用 $3{\rm{\lambda }}$表示 ${\rm{\lambda }}$和其2个相邻二元区间的并集，即 $3{{\rm{\lambda }}_{j,k}} = $ $ {{\rm{\lambda }}_{j,k - 1}} \cup {{\rm{\lambda }}_{j,k}} \cup {{\rm{\lambda }}_{j,k + 1}}$. 定义 ${L_{\rm{\lambda }}} \equiv {L_X}(j,k) = \mathop {\sup }\limits_{{\rm{\lambda '}} \subset 3{\rm{\lambda }}} \left| {{d_{X,{\rm{\lambda '}}}}} \right|$ 为 ${d_X}(j,k)$ 的Wavelet leader. 其中， ${d_{X,{\rm{\lambda '}}}} $为X在区间 ${\lambda '} $中的小波系数.

定义Wavelet leader的结构函数为 ${S_{\rm{L}}}(q,j)$， ${S_{\rm{L}}}(q,j) = \dfrac{1}{{{n_j}}}\displaystyle\sum\limits_{k = 1}^{{n_j}} {\mathop {\left| {{L_X}(j,k)} \right|}\nolimits^q } $， ${n_j}$为第j个尺度的Wavelet leader的个数，q为计算多分辨量矩的阶数. 定义对应的尺度指数 ${\zeta _{\rm{L}}}\left( q \right) = \mathop {\lim \;\inf }\limits_{j \to 0}\; \left( {{j}^{-1}{{\log }_2}\;{S_{\rm{L}}}\left( {q,j} \right)} \right)$. 结构函数 ${S_{\rm{L}}}(q,j)$被视为集合平均数的样本均值估计量，因此，尺度指数可以扩展为

(2) $ {\zeta _{\rm{L}}}(q) = \sum\limits_{p = 1}^\infty {{c_p}{{q^p}/{p!}}} . $

式中： ${c_p}$为对数累积量.

为了简化计算，采用Chhabra算法获得直接计算多重分形特征的经验公式^[14]:

(3) $ D(q) = \sum\limits_{j = {j_1}}^{{j_2}} {{w_j}{U^{\rm{L}}}(q,j)} , $

(4) $ h(q) = \sum\limits_{j = {j_1}}^{{j_2}} {{w_j}{V^{\rm{L}}}(q,j)} . $

式中：D（q）为多分形谱；h（q）为奇异指数； ${j_1}$、 ${j_2}$分别为参与估计的最小尺度与最大尺度； ${w_j}$为尺度 $j$的权重，且 $\displaystyle\sum\nolimits_{j = {j_1}}^{{j_2}} {j{w_j} \equiv 1} $， $\displaystyle\sum\nolimits_{j = {j_1}}^{{j_2}} {{w_j} \equiv 0} $， ${w_j} =$ $ {b_j}({V_0}j - {V_1})/ ({V_0}{V_2} - {V_1}^2)$，其中， ${V_i} = \displaystyle\sum\nolimits_{j = {j_1}}^{{j_2}} {{j^i}{b_j}}(i = 0,$1,2)， ${b_j}$为尺度 $j$的加权值，通常选取 ${b_j} = {n_j}$； ${U^{\rm{L}}}(q,j)$、 ${V^{\rm{L}}}(q,j)$为统计量，表达式如下：

(5) $ {U^{\rm{L}}}(q,j) = \sum\limits_{k = 1}^{{n_j}} {R_X^q(j,k){{\log }_2}\;R_X^q(j,k) + {{\log }_2}\;{n_j}} , $

(6) $ {V^{\rm{L}}}(q,j) = \sum\limits_{k = 1}^{{n_j}} {R_X^q(j,k){{\log }_2}}\; {L_X}(j,k), $

(7) $ R_X^q\left( {j,k} \right) = \left({L_X}{\left( {j,k} \right)}\right)^q{\Biggr /}\sum\limits_{k = 1}^{{n_j}} \left({L_X}{\left( {j,k} \right)}\right)^q . $

ISOMAP算法^[15]在多维尺度分析（multidimensional scaling，MDS）算法基础上加以改进，并基于两点间降维前后距离保持不变的准则，以测地线距离替代欧氏距离，较好地保留高维数据的空间特征，避免有效信息的流失. 具体实现过程如下.

1）设高维集 ${{X}} = \left\{ {{{{x}}_{{i}}}\left| {{{{x}}_{{i}}}} \right. \in {{\bf{R}}^{c}}, i = 1,2,\cdots,j,\cdots,n} \right\}$，c为空间维数，计算各点到x_i的欧式距离d（i，j）；

2）选取k个距离x_i最近的点作为近邻点，记为 ${N_i} = \left\{{ {{x}}_i^1,{{x}}_i^2,\cdots,{{x}}_i^l,\cdots,{{x}}_i^k}\right\} $；

3）以采样点x_i为节点，欧氏距离d（i，l）为边，构造采样点间邻域关系图G;

4）利用Dijkstra算法^[16]计算邻域图G中任意两点间的最短距离d_G（i，j），以逼近高维空间上测地线距离D_G，测地线距离矩阵D={d_G（i，j）}；

5）应用MDS算法^[17]对高维数据集X进行降维处理，过程如下.

设高维矩阵D对应的低维空间距离矩阵为 ${{Z}} \in {{\bf{R}}^{d'}}$，其中d'为低维空间维数. 根据任两点间欧式距离保持不变的原则，即 $\left\| {{{{z}}_i} - {{{z}}_j}} \right\| = {d_G}\left( {i,j} \right)$，构造矩阵B=Z^TZ， ${b_{ij}} = {{z}}_i^{\rm{T}}{{{z}}_j}$，其中，B为低维内积矩阵，则

(8) $ d_G^{\rm{2}}\left( {i,j} \right) = {\left\| {{{{z}}_i}} \right\|^2} + {\left\| {{{{z}}_j}} \right\|^2} - 2{{z}}_i^{\rm{T}}{{{z}}_j} = {b_{ii}} + {b_{jj}} - 2{b_{ij}}. $

对Z进行中心化可以得到

(9) $ \left. \begin{array}{l} \displaystyle\sum\limits_{i = 1}^m {d_G^{\rm{2}}\left( {i,j} \right) = {\rm{tr}}\;({{B}}) + m{b_{jj}}} ,\\ \displaystyle\sum\limits_{j = 1}^m {d_G^{\rm{2}}\left( {i,j} \right) = {\rm{tr}}\;({{B}}) + m{b_{ii}}} ,\\ \displaystyle\sum\limits_{i = 1}^m {\displaystyle\sum\limits_{j = 1}^m {d_G^{\rm{2}}\left( {i,j} \right)} } = 2m {\rm{tr}}\;({B}). \end{array} \right\} $

式中：m为距离矩阵D的维数；tr (∙)表示矩阵的迹，且 ${\rm{tr}} \;({{B}}) = {\sum\nolimits_{i = 1}^m {\left\| {{{{z}}_i}} \right\|} ^2}$. 令

(10) $ \left. \begin{array}{l} d_G^2(i, \cdot ) = \dfrac{1}{m}\displaystyle\sum\limits_{j = 1}^m {{d^2_G}(i,j)} , \\ d_G^2( \cdot ,j) = \dfrac{1}{m}\displaystyle\sum\limits_{i = 1}^m {{d^2_G}(i,j)} , \\ d_G^2( \cdot , \cdot ) = \dfrac{1}{{{m^2}}}\displaystyle\sum\limits_{i = 1}^m {\displaystyle\sum\limits_{j = 1}^{{m}} {{d^2_G}(i,j)} } . \\ \end{array} \right\} $

联立式（9）、（10），可以得到：

(11) $ {b_{ij}} = - \frac{1}{2}(d_G^2(i,j) - d_G^2(i, \cdot ) - d_G^2( \cdot ,j) + d_G^2( \cdot , \cdot )). $

根据式（11），可通过D求出B. 对B进行特征值分解： ${{B}} = {{V}}{{\varLambda}}{{{V}}^{\rm{T}}}$，其中 ${{\varLambda }} = {\rm{diag}}\;[{\lambda _1},{\lambda _2},\cdots,{\lambda _d}]$， ${\lambda _i}$、V分别为B的特征值和特征向量矩阵，且 ${\lambda _1} \geqslant {\lambda _2} \geqslant \cdots \geqslant $ ${\lambda _i} \geqslant \cdots \geqslant {\lambda _d} $. 由此得到高维数据集X的低维嵌入空间Y， $Y = \{ {{{y}}_i}\left| {{{{y}}_i} = \sqrt {{\lambda _i}} } \right.{{{v}}_i},{{{v}}_i}{\text{为组成矩阵}} $ ${{V}}\;{\text{的第}}\;i\;{\text{个特征}} $ ${\text{向量}}\} $.

由2.1节可知，低维空间集合Y能否有效表征高维空间X的结构特征，取决于近邻参数k的取值，且由残差e作为评价指标^[18]，e=1−Corr² （D，Z）. Corr (·)为标准线性相关系数. e越趋于0，表明数据嵌入效果越好. 为了获取最优参数k_opt，提出混合灰狼算法（hybrid grey wolf optimization algorithm，HGWO）进行ISOMAP参数寻优.

灰狼算法（GWO）^[19]是模拟狼群捕食的仿生智能算法. 随机产生一个狼群位置P={P₁，P₂， $\cdots $，P_N}作为候选解，灰狼个体（单个解）位置P_i=[p_i1，p_i2， $\cdots $，p_ij， $\cdots $，p_id]，其中i=1，2， $\cdots $，N，N为狼群数，d为灰狼个体位置维数，通过迭代更新方式确定猎物位置，即最优解. 实现过程如下：

(12) $ {D_{\rm{d}}} = \left| {{{C}}{{{P}}_{\rm{p}}}(t) - {{P}}(t)} \right|, $

(13) $ {{P}}(t + 1) = {{{P}}_{\rm{p}}}(t) - {{A}}{D_{\rm{d}}}. $

式中：D_d为灰狼与猎物之间的距离；P_p为猎物位置；P为灰狼位置；t为当前迭代次数；A、C为参数向量.

狼群狩猎（优化）行为主要由α狼（最优解）、β狼（次优解）和δ狼（第3优解）主导. 在狩猎时个体位置因猎物位置改变而发生变化，因此狼群根据（α，β，δ）更新后的位置重新确定猎物（最优解）位置，更新公式如下：

(14) $ \left. \begin{array}{l} {D_{\rm{\alpha }}} = \left| {{{{C}}_1}{{{P}}_{\rm{\alpha }}}\left( {{t}} \right) - {{P}}\left( {{t}} \right)} \right|,\\ {D_{\rm{\beta }}} = \left| {{{{C}}_2}{{{P}}_{\rm{\beta }}}\left( {{t}} \right) - {{P}}\left( {{t}} \right)} \right|,\\ {D_{\rm{\delta }}} = \left| {{{{C}}_3}{{{P}}_{\rm{\delta }}}\left( {{t}} \right) - {{P}}\left( {{t}} \right)} \right|, \end{array} \right\} $

(15) $ \left. \begin{array}{l} {{{P}}_1} = {{{P}}_{\rm{\alpha }}}\left( {{t}} \right) - {{{A}}_1}{D_{\rm{\alpha }}},\\ {{{P}}_2} = {{{P}}_{\rm{\beta }}}\left( {{t}} \right) - {{{A}}_2}{D_{\rm{\beta }}},\\ {{{P}}_3} = {{{P}}_{\rm{\delta }}}\left( {{t}} \right) - {{{A}}_3}{D_{\rm{\delta }}}, \end{array} \right\} $

(16) $ {{{P}}_{\rm{p}}}({{t}} + 1) = \left( {{{{P}}_1} + {{{P}}_2} + {{{P}}_3}} \right)/3. $

式中：D_α、D_β、D_δ分别为α、β、δ狼到狼群的距离，P_α、P_β、P_δ分别为α、β、δ狼的位置，A₁、A₂、A₃、C₁、C₂、C₃分别为参数向量.

GWO在寻找最优位置时依赖参数A，当1<|A|<2时，狼群会偏离猎物，导致GWO陷入局部最优解. 差分进化（differential evolutionary，DE）算法^[20]具备全局搜索能力和鲁棒性强的特点，因此可将DE算法与GWO算法相结合，构造混合灰狼优化算法HGWO，以获取ISOMAP算法的最优参数. 在DE算法中，通过变异、交叉和选择实现对（α，β，δ）的位置更新，避免GWO陷入局部最优，实现过程如下：

(17) $ {{{H}}_i}\left( {t + 1} \right) = {{{P}}_{r1}}\left( t \right) + F \left( {{{{P}}_{r2}}\left( t \right) - {{{P}}_{r3}}\left( t \right)} \right), $

(18) $ {u_{ij}}(t + 1) = \left\{ \begin{array}{l} {h_{ij}}\left( {t + 1} \right),\;\;\;{\rm{rand}} < C\;{\text{或}}\;j = {\rm{sn}};\\ {p_{ij}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{其他}, \end{array} \right. $

(19) $ {{{P}}_{\rm{p}}}\left( {t + 1} \right) = \left\{ \begin{array}{l} {{{U}}_i}\left( {t + 1} \right),\;\;f\left( {{{{U}}_i}\left( {t + 1} \right)} \right) \leqslant f\left( {{{{P}}_i}\left( t \right)} \right);\\ {{{P}}_i}\left( t \right),\;\;\;\;\;\;\;\;\text{其他}. \end{array} \right. $

式中：H_i=[h_i1，h_i2， $\cdots $，h_ij， $\cdots $，h_id]为变异个体的位置向量，U_i=[u_i1，u_i2， $\cdots $，u_ij， $\cdots $，u_id]为子代灰狼的位置向量，F∈[0，2]为缩放因子，r₁、r₂、r₃为互不相等的随机整数，P_r1、P_r2、P_r3为从狼群中任选的3个灰狼位置，C为交叉概率，rand∈[0，1]为随机数，sn为随机维数，f (·)为贪婪选择因子.

选择的适应度函数表达式为

(20) $ {f_{\rm{fit}}} = 1 - {\rm{Cor}}{{\rm{r}}^2}\;\left( {{{D}},{{Z}}} \right). $

ISOMAP的参数优化实现过程如下.

1）初始化灰狼参数（N，d，t）.

2）For i=1，2，···，t

根据式（14）计算狼群与α狼、β狼、δ狼的距离，并根据式（15）、（16）更新α狼、β狼、δ狼与猎物的位置.

根据式（20）计算父代灰狼的个体适应度，将前3个最优解标记为（P_α，P_β，P_δ），并将适应度最优值记为父代狼群最优解.

根据式（17）、（18）对父代狼群变异、交叉操作，获得子代狼群U_i.

根据式（20）计算每个子代灰狼的个体适应度，并将最优值记为子代狼群最优解.

根据式（19）筛选出子代和父代狼群中的优秀个体位置，由此确定猎物位置（最优解）.

END

3）求出狼群最优解，确定最优参数k_opt.

如图1所示为回转支承自适应特征提取方法的流程图，实施步骤如下. 1）从全寿命加速试验的回转支承振动信号中采集数据信息，确定能有效表征回转支承运行状态的数据，对数据进行集合经验模态分解-相关系数（ensemble empirical mode decomposition- correlation coefficient，EEMD-CC）^[21-22]降噪处理. 2）对实验过程中采集到的正常状态、螺栓破坏状态以及外圈破坏状态的数据进行Wavelet leader多分形分析，分别计算出3种状态下的多分形谱值D（q）、奇异指数h（q）以及尺度指数 ${\zeta _{\rm{L}}}(q)$，将其重组构造特征向量组. 3）利用HGWO算法优化ISOMAP算法的参数，利用满足要求的ISOMAP模型对重构特征向量组进行特征筛选. 4）利用步骤3）中筛选出的特征向量作为训练样本对经遗传算法参数优化后的最小二乘支持向量机进行训练建模. 5）将测试样本输入到训练好的GA-LSSVM模型中，进行故障状态识别.

为了验证所提方法的可行性，通过课题组自主研发的回转支承综合性能试验台模拟实际工况，实现振动信号采集. 为了加速试验，对待测设备施加100%极限载荷：轴向力为171.758 kN，倾覆力矩为482.618 kN∙m，径向力影响较小，忽略不计.

风电回转支承试验台主要由机械结构、液压驱动以及测控系统三大部分组成. 实验台的实物图如图2所示. 机械结构用于待测设备的安装、支撑以及固定；液压部分主要包括液压马达和液压缸，通过液压马达和液压缸联合作用模拟回转支承的实际工况；测控系统通过控制液压站完成试验台的动力输送，并对试验台的工作状态进行实时监测. 测控系统的流程图如图3所示. 测控系统采用美国NI公司LabVIEW软件平台及相应的cDAQ数采设备通过OPC协议实现与西门子公司S7-200控制系统的通信，以实现测控一体化功能.

针对某型号回转支承进行全寿命加速试验，回转支承结构参数如表1所示. 回转支承在工作时，承受轴向载荷、径向载荷以及倾覆力矩，回转支承不同部位的受力情况有所不同，因此，在回转支承上沿周向间隔90°均匀布置传感器. 传感器编号分别为1、2、3、4，设置采样频率为2 048 Hz. 完成为期7 d的回转支承全寿命加速试验，借助测试系统完成信息的采集. 为了反映回转支承不同时刻的运行状况，每隔2 h取一次样本，取样时间为20 s. 对4组加速度传感器采集的数据进行观察分析，发现4号传感器采集的信号变化最为明显，因此对4号传感器采集的信息进行分析. 回转支承全寿命振动信号（部分数据）如图4所示. 图中，A_m为加速度幅值，T为测试时间. 在实验结束后拆机检查，发现回转支承外圈滚道及螺栓出现严重破坏.

1）对采集的信号进行EEMD-CC降噪处理. 振动信号降噪前后的对比结果如图5所示. 可以看出，降噪能初步去除干扰信号，保留有用信息. 2）对降噪后的振动信号进行多分形分析. 通过Wavelet leader方法分别计算出正常状态、螺栓破坏以及外圈破坏3种状态下振动信号的尺度指数、奇异指数以及多分形谱值，计算结果如图6所示（以10个样本为例）. 如图6（a）所示，在同一状态下尺度指标与阶数存在线性递增关系，且不同阶数对应的尺度指标均聚集在同一数值附近；在不同状态下，同一阶数对应的尺度指标的聚集点有明显的差别（（0，0）点除外）. 如图6（b）所示，在同一状态下奇异指数随阶数线性递减，不同阶数对应的奇异指数聚集在同一数值附近，不同状态下的奇异指数的分布范围具有明显的差异性. 如图6（c）所示，多分形谱与奇异指数呈现凸形曲线关系，同一状态下的D (q)-h (q)图集聚在同一范围内，不同状态下的D (q)-h (q)图分布较分散.

如图7所示为正常状态下某一样本的D (q)-h (q)图. 可以看出，回转支承振动信号的多分形谱与奇异指数具有凸曲线关系，结合图6（c）可知，图7中11个点均分布在一定的范围内，并且不同状态下的分布范围具有明显的区别. 为了定量说明3种状态下各点的分布情况，对3种状态的样本计算结果进行统计. 统计结果如图8所示（以初始点、最高点及终止点3个特殊点为例）. 图中，N_u为频数，R_a为分布区间. 可以看出，在正常状态下初始点、最高点、终止点的主要分布范围分别为0.05~0.15、0.22~0.29、0.32~0.40. 在螺栓破坏状态下初始点、最高点、终止点的主要分布范围分别为0.04~0.07、0.07~0.10、0.12~0.20. 在外圈破坏状态下初始点、最高点、终止点的主要分布范围分别为0.35~0.54、0.57~0.65、0.72~0.88，如表2所示.

奇异指数反映不同测度回转支承振动信号的随机性，奇异指数越大，表明振动信号越不规则，随机性越强. 多分形谱表征奇异指数所在测度的差别（见图7）. 计算回转支承振动信号在3种状态下的多分形谱、奇异指数和尺度指数，构造1×33维特征矩阵T_j=[h_i(q)，D_i(q)， ${\zeta _{\rm{L}}}_i(q)$]，i为阶数点，i=1，2， $\cdots $，11，j为样本数.

高维特征能全面包含信号信息，但信息重叠问题必然存在. 为了去除冗余信息，降低特征维数，以解决数据处理工作量大、耗时长的问题，利用所提出的降维方法对高维向量特征进行筛选. 经HGWO优化后得到的ISOMAP近邻参数k=96，ISOMAP自适应筛选特征结果如图9所示. 由图可知，当特征维数M=4时，剩余方差e达到最小值，且随着维数的增加，剩余方差趋于稳定，由此可知，当特征维数降到4时，所得数据能基本包含原特征矩阵的全部信息. 对降维后的数据进行可视化处理，结果如图10所示. 图中，f₁为降维后第1维数据，f₂为降维后第2维数据，f₃为降维后第3维数据，通过标记大小表示降维后第4维数据f₄. 可以看出，3类状态的区分度较大，其中外圈故障状态与其他状态间的区分最为明显. 正常状态与螺栓破坏状态下的数据虽有部分重合，但整体区分度仍较大. 为了进一步分析降维后得到的4维数据的聚类效果，选取4维数据中任意3维数据进行分析，绘制聚类图，如图11所示. 结果表明，f₁-f₂-f₃和f₁-f₂-f₄聚类效果显著，类间的区分度大；f₂-f₃-f₄类间区分度较小，类间存在严重的重叠现象.

为了定量分析HGWO-ISOMAP的特征筛选效果，将降维后的特征矩阵分成训练样本和测试样本两部分，训练样本的输入为某状态降维后的特征矩阵，输出为该状态对应的状态函数值f，定义f={1，2，3}（f=1、2、3分别表示正常状态、螺栓破坏状态、外圈破坏状态），测试样本的输入为某一降维后的特征矩阵. 由此，总共形成了300组训练样本和150组测试样本. 通过训练样本对LSSVM进行训练建模，将测试样本输入训练好的LSSVM模型中，以实现定量评价特征筛选效果的目的. 为了提高模型的识别精度，采用遗传算法对LSSVM模型进行参数优化，初始种群大小为80，最大迭代次数为100，最终获得的最优参数 $\sigma {\rm{ = 0}}{\rm{.05,}}\;\gamma {\rm{ = 9}}{\rm{.07}}$；利用优化后的参数对LSSVM模型进行训练建模，最后利用测试样本检验模型的识别精度，识别结果如图12(a)所示. 如图12(b)~(e)所示分别为降维前分类结果和降维后任选取3维数据进行分类的结果. 图中，横坐标为测试样本点N_m，纵坐标为对应的状态值f. 由分类结果可知，采用HGWO-ISOMAP方法提取特征的识别效果最好，仅1个点发生误判. 5种分类结果如表3所示. 表中，R_c为识别精度，t为时间. 由识别精度可知，HGWO-ISOMAP筛选后的Wavelet leader特征具有较高的识别正确率，识别精度为99.33%；在此基础上任选择的三维特征也具有较高的识别率，除了f₂-f₃-f₄识别率为88%外，其他均大于95%，比未经特征筛选的33维Wavelet leader的特征识别效果更好. 由此可见，经过HGWO-ISOMAP筛选后的特征在保留原有特征有效信息的基础上，极大地剔除了冗余信息，减少了计算量，缩短了计算时间，提高了计算效率. 如表4所示为采用时域-频域特征^[23]、时频域特征^[24]以及BP神经网络作为分类器时的分类结果. 结合表3、4的分类结果可知，基于Wavelet leader提取的多分形特征在回转支承故障状态的识别上具有较好的效果. 同时，HGWO-ISOMAP自适应筛选特征方法能进一步改进故障状态识别性能.

（1）基于Wavelet leader方法挖掘的回转支承振动信号具有显著的多分形特征，在3种状态下，多分形谱-奇异指数关系图具有较大的区分度，可以作为回转支承故障状态识别的可靠特征指标.

（2）HGWO-ISOMAP方法能自适应性筛选出多分形特征矩阵中的有效信息，在保留原特征矩阵中全部有效信息的同时，剔除干扰信息，降低特征维数，显著提高计算效率和识别准确率.

（3）对本研究所提出的特征提取方法进行GA-LSSVM模型定量分析，进一步证明了所提方法的优越性. 所提方法为回转支承振动信号特征提取的新型有效分析方法. 这对回转支承运行状态在线监测和故障诊断具有一定的工程意义.