基于多域流形的行星齿轮箱局部故障识别

工程科学学报，第39卷，第5期：769-777.2017年5月 Chinese Journal of Engineering,Vol.39,No.5:769-777,May 2017 D0I:10.13374/j.issn2095-9389.2017.05.016;htp://journals.ustb.edu.cn 基于多域流形的行星齿轮箱局部故障识别 赵 川，冯志鹏区 北京科技大学机械工程学院，北京100083 ☒通信作者，E-mail:fengzp@ustb.cdu.cn 摘要行星齿轮箱振动信号包含多种频率成分和噪声干扰，频谱具有复杂的边带结构，容易对故障识别造成误导甚至引起 错判.在不同故障状态下，行星齿轮箱振动信号的多域特征量将偏离正常范围且偏离程度不同，根据这一特点，提取振动信 号的时域，频域特征参量用于故障识别.为了避免传统分析方法中负频率及虚假模态问题，增强对噪声干扰的鲁棒性，采用 局部均值分解法将信号自适应地分解为单分量之和，提取时频域单分量瞬时幅值能量.针对多域特征空间构造过程中出现 的高维及非线性问题，采用流形学习对数据进行降维处理.提出基于改进的虚假近邻点的本征维数估计及最优k邻域确定方 法，并通过等距映射对多域特征空间进行降维分析.对于行星齿轮箱实验信号，根据样本流形特征聚类结果，分别识别出了 太阳轮、行星轮和齿圈的局部故障，从而验证了上述方法的有效性 关键词行星齿轮箱；故障识别；局部均值分解；本征维数估计：多域流形 分类号TH17 Localized fault identification of planetary gearboxes based on multiple-domain manifold ZHAO Chuan,FENG Zhi-peng School of Mechanical Engineering,University of Science and Technology Beijing,Beijing 100083,China Corresponding author,E-mail:fengzp@ustb.edu.cn ABSTRACT The vibration signals of planetary gearboxes are composed of complex frequency components and interfering noises,and their spectra have intricate sidebands,which cause difficulty in and even misleading fault identification.In different fault cases,the vibration signatures in multiple domains typically differ from normal states with different discrepancies.Based on this hypothesis,time and frequency domain features are extracted for the purposes of fault identification.The vibration signal is adaptively decomposed into a set of mono-components,and the instantaneous energy of each mono-component is calculated in time-frequency domain by exploiting the merits of local mean decomposition,including its better robustness to noise and freedom from pseudo-mode and negative frequency problems.Manifold learning is utilized to tackle the high-dimensionality and non-linearity aspects of multiple-domain feature space construction.A new method is proposed for estimating the intrinsic dimension and selecting the k-nearest neighborhood based on the improved pseudo-nearest neighbor.In addition,isometric feature mapping (ISOMAP)is utilized to reduce the dimensions of the mul- tiple-domain feature space.The proposed method is validated by analyzing the planetary gearbox lab experimental dataset.Based on the clustering analysis results of the extracted manifold features,the localized faults on the sun,planet,and ring gears are successfully identified. KEY WORDS planetary gearboxes;fault identification;local mean decomposition;estimation of intrinsic dimension;multiple-do- main manifold 行星齿轮箱作为传动系统中的关键部件在风力发电机组中应用广泛.风力发电机组的工作环境造成行 收稿日期：2016-06-24 基金项目：国家自然科学基金资助项目(11272047,51475038)

工程科学学报,第 39 卷,第 5 期:769鄄鄄777,2017 年 5 月 Chinese Journal of Engineering, Vol. 39, No. 5: 769鄄鄄777, May 2017 DOI: 10. 13374 / j. issn2095鄄鄄9389. 2017. 05. 016; http: / / journals. ustb. edu. cn 基于多域流形的行星齿轮箱局部故障识别 赵 川, 冯志鹏苣 北京科技大学机械工程学院, 北京 100083 苣 通信作者, E鄄mail: fengzp@ ustb. edu. cn 摘 要 行星齿轮箱振动信号包含多种频率成分和噪声干扰,频谱具有复杂的边带结构,容易对故障识别造成误导甚至引起 错判. 在不同故障状态下,行星齿轮箱振动信号的多域特征量将偏离正常范围且偏离程度不同,根据这一特点,提取振动信 号的时域、频域特征参量用于故障识别. 为了避免传统分析方法中负频率及虚假模态问题,增强对噪声干扰的鲁棒性,采用 局部均值分解法将信号自适应地分解为单分量之和,提取时频域单分量瞬时幅值能量. 针对多域特征空间构造过程中出现 的高维及非线性问题,采用流形学习对数据进行降维处理. 提出基于改进的虚假近邻点的本征维数估计及最优 k 邻域确定方 法,并通过等距映射对多域特征空间进行降维分析. 对于行星齿轮箱实验信号,根据样本流形特征聚类结果,分别识别出了 太阳轮、行星轮和齿圈的局部故障,从而验证了上述方法的有效性. 关键词 行星齿轮箱; 故障识别; 局部均值分解; 本征维数估计; 多域流形 分类号 TH17 收稿日期: 2016鄄鄄06鄄鄄24 基金项目: 国家自然科学基金资助项目 (11272047,51475038) Localized fault identification of planetary gearboxes based on multiple鄄domain manifold ZHAO Chuan, FENG Zhi鄄peng 苣 School of Mechanical Engineering, University of Science and Technology Beijing, Beijing 100083, China 苣 Corresponding author, E鄄mail: fengzp@ ustb. edu. cn ABSTRACT The vibration signals of planetary gearboxes are composed of complex frequency components and interfering noises, and their spectra have intricate sidebands, which cause difficulty in and even misleading fault identification. In different fault cases, the vibration signatures in multiple domains typically differ from normal states with different discrepancies. Based on this hypothesis, time and frequency domain features are extracted for the purposes of fault identification. The vibration signal is adaptively decomposed into a set of mono鄄components, and the instantaneous energy of each mono鄄component is calculated in time鄄frequency domain by exploiting the merits of local mean decomposition, including its better robustness to noise and freedom from pseudo鄄mode and negative frequency problems. Manifold learning is utilized to tackle the high鄄dimensionality and non鄄linearity aspects of multiple鄄domain feature space construction. A new method is proposed for estimating the intrinsic dimension and selecting the k鄄nearest neighborhood based on the improved pseudo鄄nearest neighbor. In addition, isometric feature mapping (ISOMAP) is utilized to reduce the dimensions of the mul鄄 tiple鄄domain feature space. The proposed method is validated by analyzing the planetary gearbox lab experimental dataset. Based on the clustering analysis results of the extracted manifold features, the localized faults on the sun, planet, and ring gears are successfully identified. KEY WORDS planetary gearboxes; fault identification; local mean decomposition; estimation of intrinsic dimension; multiple鄄do鄄 main manifold 行星齿轮箱作为传动系统中的关键部件在风力发 电机组中应用广泛. 风力发电机组的工作环境造成行

·770· 工程科学学报，第39卷，第5期 星齿轮箱承受动态重载负荷且运行工况变化频繁，存 因及变化规律进行研究，揭示故障类别及零部件的损 在多种因素容易导致行星齿轮箱的太阳轮、行星轮及 伤过程.由于行星齿轮箱主要由太阳轮、行星轮、齿圈 齿圈形成点蚀、脱落及裂纹等局部故障.局部故障最 和行星架等构成，各组成部分间的运动及振动存在着 终导致齿轮箱失效，造成损失口.因而有必要对行星 相互影响，这必然使得局部故障对于整个箱体运行状 齿轮箱的局部故障识别进行研究. 态的影响包含繁杂的耦合性成分.时域信号为在不同 行星齿轮箱振动信号成分复杂，主要包括旋转频时刻测量到的状态参数，在进行信号分析时显得直观， 率、啮合频率、倍频成分、边带成分、特征频率成分等 易于理解和实现.当行星齿轮箱发生局部故障时，振 在行星齿轮箱运行过程中，某个轮齿上的局部故障将 动信号在时域的幅值与相位会偏离正常状态[).但 对啮合振动产生调幅调频作用：齿轮局部故障相对传 仅仅做时域分析，所获得的齿轮箱工作状态信息是有 感器的振动传递路径的时变现象将对局部损伤引起的 限的，常常只能判断运行状态的正常与异常，对于异常 冲击产生额外的调幅作用].对于行星齿轮箱的故障 状态的类别、发生位置及原因很难给出有效的判断 识别，研究人员进行了各种探索，提出了时间平均、约 频域分析需要对时域信号进行傅里叶变换完成从时域 束自适应、Winger Ville分布、循环平稳、谱峭度、能量特 到频域的转化.信号频域信息包含了频率、幅值及相 征及灰色关联度、迭代广义压缩变换等识别方法[3-]. 位等内容，由文献[2]知行星齿轮箱振动信号频谱结 冯志鹏等)构建了行星齿轮箱局部故障振动信号模 构复杂，包含了丰富的频率信息，且不同故障条件下的 型，推导了频谱解析表达式并给出了相应部件的故障频谱幅值和分布也会发生变化)，为进行齿轮箱运行 特征频率的计算公式，利用频谱分析方法进行了行星 状态的判断提供了更丰富的依据，但频域分析不能准 齿轮箱的局部故障识别. 确反映信号能量在空间各尺度的分布规律，且在故障 频谱分析是故障识别的常用方法.由于行星齿轮 早期阶段仅通过时域和频域分析很难识别故障模式. 箱的振动信号成分复杂，应用频谱分析法进行行星齿 当故障或异常运行状态发生时，被测振动信号常为非 轮箱的故障识别时，需要对振动信号进行复杂的频谱 线性信号，基于单分量的瞬时幅值能量分析可以弥补 结构分析，特征频率的计算.由于频率成分多样，且不 时域及频域分析的不足，为进一步研究齿轮箱运行状 同的结构参数将对频谱结构造成影响，因而给故障识 态提供了更完整的信息[4).对齿轮箱振动样本信号 别造成了一定的难度 的多域分析必然导致样本空间的维数增加.通常齿轮 故障识别本质上是对设备的运行模式进行识 箱运行状态信息分布于嵌入在多域高维空间中的非线 别].近年来，流形学习作为一种高维非线性数据的 性低维流形上，因而有效实现多域高维空间构建、降维 降维工具，在模式识别中得到了广泛的应用.汤宝平 及提高识别性能是进行行星齿轮箱局部故障识别的重 等)提出的增殖正交邻域保持嵌入流形学习算法实 要任务.基于多域流形的行星齿轮箱局部故障识别模 现了对于新增样本的动态增殖学习.陈法法等[]通 型主要分为四个步骤：信号采集、特征提取、降维和模 过时域统计指标和内禀模态分量能量构建了齿轮箱故 式识别.实现过程如图1所示. 障振动特征空间.栗茂林等]利用基于时域指标和 时 小波频带能量组成特征空间的局部切空间排列算法提 域 多 取了滚动轴承的早期故障特征.综上所述，由测试信 域高维 测试样本集 号构建的高维特征空间和与系统状态相关且由流形学 高维特 习发掘的特定流形相应.系统故障对低维流形的分布 析 趋势和点的聚合程度造成影响. 空间构建 针对上述情况，本文提出基于多域流形的故障识 空间降维 练样本集 诊断模型 结果判定 别方法，对观测的振动信号构造时域、频域参量，利用 析 局部均值分解提取时频域单分量瞬时幅值能量，构建 图1行星齿轮箱故障识别流程图 多域高维空间，通过等距映射法进行降维，识别不同状 Fig.1 Fault identification flow chart of planetary gearboxes 态样本低维流形的聚类结果，从而实现对行星齿轮箱 的故障识别. 2 多域高维特征空间构建及流形学习 1多域流形故障识别机理 2.1多域高维特征空间构建 行星齿轮箱局部故障识别旨在对其局部故障的起 近年来，基于统计量特征的故障识别得到广泛关

工程科学学报,第 39 卷,第 5 期 星齿轮箱承受动态重载负荷且运行工况变化频繁,存 在多种因素容易导致行星齿轮箱的太阳轮、行星轮及 齿圈形成点蚀、脱落及裂纹等局部故障. 局部故障最 终导致齿轮箱失效,造成损失[1] . 因而有必要对行星 齿轮箱的局部故障识别进行研究. 行星齿轮箱振动信号成分复杂,主要包括旋转频 率、啮合频率、倍频成分、边带成分、特征频率成分等. 在行星齿轮箱运行过程中,某个轮齿上的局部故障将 对啮合振动产生调幅调频作用;齿轮局部故障相对传 感器的振动传递路径的时变现象将对局部损伤引起的 冲击产生额外的调幅作用[2] . 对于行星齿轮箱的故障 识别,研究人员进行了各种探索,提出了时间平均、约 束自适应、WingerVille 分布、循环平稳、谱峭度、能量特 征及灰色关联度、迭代广义压缩变换等识别方法[3鄄鄄9] . 冯志鹏等[1] 构建了行星齿轮箱局部故障振动信号模 型,推导了频谱解析表达式并给出了相应部件的故障 特征频率的计算公式,利用频谱分析方法进行了行星 齿轮箱的局部故障识别. 频谱分析是故障识别的常用方法. 由于行星齿轮 箱的振动信号成分复杂,应用频谱分析法进行行星齿 轮箱的故障识别时,需要对振动信号进行复杂的频谱 结构分析,特征频率的计算. 由于频率成分多样,且不 同的结构参数将对频谱结构造成影响,因而给故障识 别造成了一定的难度. 故障识别本质上是对设备的运行模式进行识 别[10] . 近年来,流形学习作为一种高维非线性数据的 降维工具,在模式识别中得到了广泛的应用. 汤宝平 等[11]提出的增殖正交邻域保持嵌入流形学习算法实 现了对于新增样本的动态增殖学习. 陈法法等[12] 通 过时域统计指标和內禀模态分量能量构建了齿轮箱故 障振动特征空间. 栗茂林等[13] 利用基于时域指标和 小波频带能量组成特征空间的局部切空间排列算法提 取了滚动轴承的早期故障特征. 综上所述,由测试信 号构建的高维特征空间和与系统状态相关且由流形学 习发掘的特定流形相应. 系统故障对低维流形的分布 趋势和点的聚合程度造成影响. 针对上述情况,本文提出基于多域流形的故障识 别方法,对观测的振动信号构造时域、频域参量,利用 局部均值分解提取时频域单分量瞬时幅值能量,构建 多域高维空间,通过等距映射法进行降维,识别不同状 态样本低维流形的聚类结果,从而实现对行星齿轮箱 的故障识别. 1 多域流形故障识别机理 行星齿轮箱局部故障识别旨在对其局部故障的起 因及变化规律进行研究,揭示故障类别及零部件的损 伤过程. 由于行星齿轮箱主要由太阳轮、行星轮、齿圈 和行星架等构成,各组成部分间的运动及振动存在着 相互影响,这必然使得局部故障对于整个箱体运行状 态的影响包含繁杂的耦合性成分. 时域信号为在不同 时刻测量到的状态参数,在进行信号分析时显得直观, 易于理解和实现. 当行星齿轮箱发生局部故障时,振 动信号在时域的幅值与相位会偏离正常状态[14] . 但 仅仅做时域分析,所获得的齿轮箱工作状态信息是有 限的,常常只能判断运行状态的正常与异常,对于异常 状态的类别、发生位置及原因很难给出有效的判断. 频域分析需要对时域信号进行傅里叶变换完成从时域 到频域的转化. 信号频域信息包含了频率、幅值及相 位等内容,由文献[2] 知行星齿轮箱振动信号频谱结 构复杂,包含了丰富的频率信息,且不同故障条件下的 频谱幅值和分布也会发生变化[14] ,为进行齿轮箱运行 状态的判断提供了更丰富的依据,但频域分析不能准 确反映信号能量在空间各尺度的分布规律,且在故障 早期阶段仅通过时域和频域分析很难识别故障模式. 当故障或异常运行状态发生时,被测振动信号常为非 线性信号,基于单分量的瞬时幅值能量分析可以弥补 时域及频域分析的不足,为进一步研究齿轮箱运行状 态提供了更完整的信息[14] . 对齿轮箱振动样本信号 的多域分析必然导致样本空间的维数增加. 通常齿轮 箱运行状态信息分布于嵌入在多域高维空间中的非线 性低维流形上,因而有效实现多域高维空间构建、降维 及提高识别性能是进行行星齿轮箱局部故障识别的重 要任务. 基于多域流形的行星齿轮箱局部故障识别模 型主要分为四个步骤:信号采集、特征提取、降维和模 式识别. 实现过程如图 1 所示. 图 1 行星齿轮箱故障识别流程图 Fig. 1 Fault identification flow chart of planetary gearboxes 2 多域高维特征空间构建及流形学习 2郾 1 多域高维特征空间构建 近年来,基于统计量特征的故障识别得到广泛关 ·770·

赵川等：基于多域流形的行星齿轮箱局部故障识别 ·771· 注，统计量特征能够发掘出振动信号中与故障相关的 循环y次，直到u,(t)成为常量或单调函数. 重要信息，为进行多域流形故障识别提供了基础.研 对于每个分量P℉的瞬时幅值能量有 究表明多域统计特征能够完整描述机械设备的运行状 1a(5)12. (5) 态信息，因而如何从振动信号中提取有效的多域特征 成为基于多域流形故障识别的重要内容.在被定义的 其中，N表示离散单分量P℉的数据长度，j表示P℉分 统计参数中，时域特征包含了六个量纲参数：均值、均 量中第广个数据点，，表示第j个数据点对应的时刻 方根植、方根幅值、标准差、歪度、峰态和六个非量纲参 构建时频域特征向量V,其中z为分解后P℉分量 数：波形因数、蜂值因数、脉冲因数、余隙因数、歪度因 个数 数、峰态因数，频域特征包含了四个统计参数：平均频 V={E,E2,…,E. (6) 率、均方根频率、中心频率、根方差频率，具体计算公式 在构建多域高维空间时，对每个样本信号进行多 参见文献[10].在基于单分量的瞬时幅值能量分析 域分析，按本文所列述的特征参量顺序进行多域统计 中，小波变换和经验模态分解是常用的信号分解方法. 量计算，并针对每个样本构建g维特征向量ω，其中9 但小波变换的结果与所选择的小波基紧密相关，且缺 为多域统计量的个数即q=6+6+4+z,w的元素记为 乏自适应性，经验模态分解法对噪音敏感且存在模态 0,c=1,2,…,9.0.的取值为对应的第c个统计量特 混合问题).程军圣等[6]对比分析了局部均值分解 征值的大小，则多域特征向量ω如式(7)所示 (local mean decomposition,LMD)和经验模态分解 ω={01,02，…，0g} (7) (empirical mode decomposition,EMD),描述了局部均 记p个样本的高维特征矩阵空间为2，x显然 值分解的优越性，并将局部均值分解方法有效的应用 P。x,不能直观的进行故障识别且维数较高，对数据的 于齿轮故障识别.局部均值分解是由Smith提出的一 后续处理造成困难.为了进一步提取低维的有效特征 种新的非线性和非平稳信号分析方法.由于局部均值 同时抑制噪音干扰，需要构建用于处理非线性特征的 分解是依据信号本身的信息进行自适应分解的，不涉 维数缩减映射. 及小波基选择与模态混合问题，产生的乘积函数P℉ 2.2等距映射降维 (production function)分量具有真实的物理意义，由此 流形学习作为处理非线性高维数据的降维方法具 得到的时频分布能够清晰准确地反映出信号能量在空 有广泛的应用，其构建降维映射的基本思路主要有两 间各尺度上的分布规律.具体计算过程如下6 种：一是利用经典多尺度降维方法将数据点投影到低 (1)找出以：为时间变量的信号x(t)的所有局部 维空间同时保持测地距离不变；二是通过分析局部的 极值点n,计算相邻极值点n,和n,的均值m,及包络 流形结构来估计全局结构.本文采用等距映射(IS0- 估计a,其中下标i表示第i个极值点 MAP)降维方法对多域高维特征空间进行进一步地特 (2)计算连续局部均值函数m(t),连续包络估 征提取.ISOMAP的本质是多维尺度变换(multidimen- 计函数a(t).其中，I表示对应第l个P℉分量，b表示 sional scaling,MDS).与MDS不同的是，ISOMAP采用 迭代次数. 测地线距离代替欧氏距离作为数据差异度量.在MDS (3)计算连续调频信号估计S(t).将m(t)从 中，利用欧氏距离作为距离度量，但只能体现数据的线 原始信号x(t)分离得h,(t),计算S:().若S,(t)的 性结构.如果原始的高维数据是高度扭曲的流形，则 包络估计函数a2(t)未满足式(2)的要求则需要重复 这种度量不能体现数据的内在结构，因而采用测地线 上述的计算过程次，直到获取一个纯调频函数 距离作为距离度量，再进行MDS降维.对于高维特征 S.(t). 空间2=[w1,w2,…,w],w.∈R,其中R为q维空 Su(t)=hu(t)/an(t), (1) 间，e=1,2,…,p.ISOMAP的应用日标是从中提取能 a2(t)=1. (2) 够表征其本质特征的4=[61，δ2，…，6]，6。∈R,其 (4)计算第一个分量PF,:PF,的包络可由式(3) 中R为d维空间且q>d.具体的实现过程如下. 计算得出.则PF,的计算表达式为(4).PF的瞬时幅 (1)计算每个点。的k个近邻点.计算每个数据 值即包络信号a,(t). 点对w和w之间的欧氏距离，记为disd(e,g),其中 e,g∈[1,2，…，p]且e≠g,确定哪些点是近邻的，构建 a,(t)=a.(t)a2(t)…a.(t)=Πa(t),(3) 近邻图G. PF,(t)=a,(t)S.(t). (4) (2)计算近邻图距离矩阵.计算近邻图G中每两 (5)将PF从原始信号x()分离得到一个新的函 点间的最短路径disg(e,g),记所得的距离矩阵为 数u(t),把u,(t)作为原始数据重复步骤(1)~(5)， Dc=disg(e,g)

赵 川等: 基于多域流形的行星齿轮箱局部故障识别 注,统计量特征能够发掘出振动信号中与故障相关的 重要信息,为进行多域流形故障识别提供了基础. 研 究表明多域统计特征能够完整描述机械设备的运行状 态信息,因而如何从振动信号中提取有效的多域特征 成为基于多域流形故障识别的重要内容. 在被定义的 统计参数中,时域特征包含了六个量纲参数:均值、均 方根植、方根幅值、标准差、歪度、峰态和六个非量纲参 数:波形因数、峰值因数、脉冲因数、余隙因数、歪度因 数、峰态因数,频域特征包含了四个统计参数:平均频 率、均方根频率、中心频率、根方差频率,具体计算公式 参见文献[10]. 在基于单分量的瞬时幅值能量分析 中,小波变换和经验模态分解是常用的信号分解方法. 但小波变换的结果与所选择的小波基紧密相关,且缺 乏自适应性,经验模态分解法对噪音敏感且存在模态 混合问题[15] . 程军圣等[16] 对比分析了局部均值分解 (local mean decomposition, LMD) 和 经 验 模 态 分 解 (empirical mode decomposition, EMD), 描述了局部均 值分解的优越性,并将局部均值分解方法有效的应用 于齿轮故障识别. 局部均值分解是由 Smith 提出的一 种新的非线性和非平稳信号分析方法. 由于局部均值 分解是依据信号本身的信息进行自适应分解的,不涉 及小波基选择与模态混合问题,产生的乘积函数 PF (production function)分量具有真实的物理意义,由此 得到的时频分布能够清晰准确地反映出信号能量在空 间各尺度上的分布规律. 具体计算过程如下[16] . (1) 找出以 t 为时间变量的信号 x(t)的所有局部 极值点 ni,计算相邻极值点 ni和 ni + 1的均值 mi及包络 估计 ai, 其中下标 i 表示第 i 个极值点. (2) 计算连续局部均值函数 mlb ( t),连续包络估 计函数 alb(t). 其中,l 表示对应第 l 个 PF 分量,b 表示 迭代次数. (3) 计算连续调频信号估计 Slb ( t). 将 m11 ( t)从 原始信号 x(t)分离得 h11 (t),计算 S11 (t). 若 S11 (t)的 包络估计函数 a12 (t)未满足式(2)的要求则需要重复 上述的计算过程 w 次, 直到获取一个纯调频函数 S1w (t). S11 (t) = h11 (t) / a11 (t), (1) a12 (t) = 1. (2) (4)计算第一个分量 PF1 . PF1的包络可由式(3) 计算得出. 则 PF1的计算表达式为(4). PF1的瞬时幅 值即包络信号 a1 (t). a1 (t) = a11 (t)a12 (t)…a1w (t) = 仪 w b = 1 a1b(t), (3) PF1 (t) = a1 (t)S1w (t). (4) (5)将 PF1从原始信号 x( t)分离得到一个新的函 数 u1 (t),把 u1 (t)作为原始数据重复步骤(1) ~ (5), 循环 y 次,直到 uy(t)成为常量或单调函数. 对于每个分量 PF 的瞬时幅值能量有 El = 移 N j = 1 | al(t j) | 2 . (5) 其中,N 表示离散单分量 PF 的数据长度,j 表示 PF 分 量中第 j 个数据点,t j 表示第 j 个数据点对应的时刻. 构建时频域特征向量 V,其中 z 为分解后 PF 分量 个数. V = {E1 ,E2 ,…,Ez}. (6) 在构建多域高维空间时,对每个样本信号进行多 域分析,按本文所列述的特征参量顺序进行多域统计 量计算,并针对每个样本构建 q 维特征向量 棕,其中 q 为多域统计量的个数即 q = 6 + 6 + 4 + z,棕 的元素记为 oc, c = 1,2,…,q. oc的取值为对应的第 c 个统计量特 征值的大小,则多域特征向量 棕 如式(7)所示. 棕 = {o1 ,o2 ,…,oq}. (7) 记 p 个样本的高维特征矩阵空间为 赘p 伊 q . 显然 赘p 伊 q不能直观的进行故障识别且维数较高,对数据的 后续处理造成困难. 为了进一步提取低维的有效特征 同时抑制噪音干扰,需要构建用于处理非线性特征的 维数缩减映射. 2郾 2 等距映射降维 流形学习作为处理非线性高维数据的降维方法具 有广泛的应用,其构建降维映射的基本思路主要有两 种:一是利用经典多尺度降维方法将数据点投影到低 维空间同时保持测地距离不变;二是通过分析局部的 流形结构来估计全局结构. 本文采用等距映射( ISO鄄 MAP)降维方法对多域高维特征空间进行进一步地特 征提取. ISOMAP 的本质是多维尺度变换(multidimen鄄 sional scaling, MDS). 与 MDS 不同的是,ISOMAP 采用 测地线距离代替欧氏距离作为数据差异度量. 在 MDS 中,利用欧氏距离作为距离度量,但只能体现数据的线 性结构. 如果原始的高维数据是高度扭曲的流形,则 这种度量不能体现数据的内在结构,因而采用测地线 距离作为距离度量,再进行 MDS 降维. 对于高维特征 空间 赘 = [棕1 ,棕2 ,…,棕p], 棕e沂R q ,其中 R q为 q 维空 间,e = 1,2,…,p. ISOMAP 的应用目标是从中提取能 够表征其本质特征的 驻 = [ 啄1 ,啄2 ,…,啄p ],啄e沂R d ,其 中 R d为 d 维空间且 q > d. 具体的实现过程如下[17] . (1)计算每个点 棕e的 k 个近邻点. 计算每个数据 点对 棕e和 棕g之间的欧氏距离,记为 disd( e,g),其中 e,g沂[1,2,…,p]且 e屹g,确定哪些点是近邻的,构建 近邻图 G. (2)计算近邻图距离矩阵. 计算近邻图 G 中每两 点间的最短路径 disg ( e, g),记所得的距离矩阵为 DG = {disg(e,g)}. ·771·

·772· 工程科学学报，第39卷，第5期 (3)求低维流形.将矩阵D。作为MDS的输入从 (4)确定本征维数及最优k值.遍历r值，步长为 而求得低维流形△，x, 1,且在r值一定时，遍历d重复计算步骤(1)、(2)和 2.3本征维数估计与k近邻点选择 (3)直到H(d)以单调递增的方式接近1且在d值大 在流形学习过程中，涉及两个重要参数即本征维 于某个特定值d'时H,(d)趋于平稳，此时d'即为本征 数d与k近邻中样本个数k的确定.d是降维的目标 维数，对应的k=d”+r. 维数，是拘束样本点保证数据内在结构所需坐标的最 3实验信号分析 小个数.k值的选取和确定对流形学习的效果产生影 响.k值过大容易导致局部几何结构的损失，而k值过3.1实验说明与信号采集 小则会引起计算错误.通常k应满足范围k≥d+1且 图2所示实验系统由电机及其控制器、一级行星 d≥1.本文基于改进的虚假近邻点法，同时对本征维 齿轮箱及振动测试系统组成.表1列出了齿轮箱的齿 数d和近邻点数k进行估计和选择.虚假近邻点法的 轮参数.多个传感器安装在齿轮箱的基座及箱体顶 基本思想是高维空间中并不相邻的两点投影到低维空 部，恰位于齿圈上方.本文分析通过箱体顶部传感器 间中可能会成为相邻的两点，而随着低维空间维数的 采集的振动信号，采集频率16384Hz,采样时间60s. 增加，虚假邻点也会被逐步剔除，从而使内嵌轨迹得到 连接齿轮箱太阳轮输入轴的旋转频率设为15.95Hz. 恢复[).方法的具体实现过程如下. 图3所示为用于仿真的带有局部磨损剥落损伤的齿轮 (1)设定本征维数为d(d≥1)，选取近邻点k= 实物照片.每个振动条件下采集100个样本，每个样 d+r,r为自然数且r≥2，通过ISOMAP对2进行降 本包含5000个离散数据点.随机选择50个样本作为 维获取低维矩阵Apx和Apxa) 训练数据，其他样本为测试数据，则四种工况下（正 (2)记⑧为4。xd中样本点6，的最近邻点，距 离为 R(d,e)=‖8-8rot‖. (8) (3)确定判断指标.定义虚假近邻点判断公式 a(d.) (9) 当a(d,e)大于阀值时，可判断虚假近邻点.但阀 值的确定涉及较多主观性.为避免主观性对判断结果 图2实验系统 的影响，定义 Fig.2 Experimental system H(d)=1R(d+1,e (10) p R(d,e) 表1行星齿轮箱参数 并按式(11)计算指标函数H(d) Table 1 Planetary gearbox configuration parameters H(d)=Hd+1) 齿轮 太阳轮行星轮（数量）齿圈 (11) H(d) 齿数 13 38(3) 92 2010/12101 图3齿轮局部损伤.(a)太阳轮：(b)行星轮：(c)齿圈 Fig.3 Localized gear damage:(a)sun gear;(b)planet gear;(c)ring gear

工程科学学报,第 39 卷,第 5 期 (3)求低维流形. 将矩阵 DG作为 MDS 的输入从 而求得低维流形 驻p 伊 q . 2郾 3 本征维数估计与 k 近邻点选择 在流形学习过程中,涉及两个重要参数即本征维 数 d 与 k 近邻中样本个数 k 的确定. d 是降维的目标 维数,是拘束样本点保证数据内在结构所需坐标的最 小个数. k 值的选取和确定对流形学习的效果产生影 响. k 值过大容易导致局部几何结构的损失,而 k 值过 小则会引起计算错误. 通常 k 应满足范围 k逸d + 1 且 d逸1. 本文基于改进的虚假近邻点法,同时对本征维 数 d 和近邻点数 k 进行估计和选择. 虚假近邻点法的 基本思想是高维空间中并不相邻的两点投影到低维空 间中可能会成为相邻的两点,而随着低维空间维数的 增加,虚假邻点也会被逐步剔除,从而使内嵌轨迹得到 恢复[18] . 方法的具体实现过程如下. (1)设定本征维数为 d( d逸1), 选取近邻点k = d + r, r 为自然数且 r逸2, 通过 ISOMAP 对 赘 进行降 维获取低维矩阵 驻p 伊 d和 驻p 伊 (d + 1) . (2)记 啄 nearest i 为 驻p 伊 d 中样本点 啄i 的最近邻点,距 离为 R(d,e) = 椰啄e - 啄 nearest e 椰. (8) (3)确定判断指标. 定义虚假近邻点判断公式 琢(d,e) = R(d + 1,e) R(d,e) . (9) 当 琢(d,e)大于阀值时,可判断虚假近邻点. 但阀 值的确定涉及较多主观性. 为避免主观性对判断结果 的影响,定义 H(d) = 1 p 移 p e = 1 R(d + 1,e) R(d,e) . (10) 并按式(11)计算指标函数 H1 (d). H1 (d) = H(d + 1) H(d) . (11) (4)确定本征维数及最优 k 值. 遍历 r 值,步长为 1,且在 r 值一定时,遍历 d 重复计算步骤(1)、(2) 和 (3)直到 H1 (d)以单调递增的方式接近 1 且在 d 值大 于某个特定值 d忆时 H1 ( d)趋于平稳,此时 d忆即为本征 维数,对应的 k = d忆 + r. 3 实验信号分析 3郾 1 实验说明与信号采集 图 2 所示实验系统由电机及其控制器、一级行星 齿轮箱及振动测试系统组成. 表 1 列出了齿轮箱的齿 轮参数. 多个传感器安装在齿轮箱的基座及箱体顶 部,恰位于齿圈上方. 本文分析通过箱体顶部传感器 采集的振动信号,采集频率 16384 Hz,采样时间 60 s. 连接齿轮箱太阳轮输入轴的旋转频率设为 15郾 95 Hz. 图 3 所示为用于仿真的带有局部磨损剥落损伤的齿轮 实物照片. 每个振动条件下采集 100 个样本,每个样 本包含 5000 个离散数据点. 随机选择 50 个样本作为 训练数据,其他样本为测试数据,则四种工况下(正 图 2 实验系统 Fig. 2 Experimental system 表 1 行星齿轮箱参数 Table 1 Planetary gearbox configuration parameters 齿轮 太阳轮 行星轮(数量) 齿圈 齿数 13 38 (3) 92 图 3 齿轮局部损伤. (a) 太阳轮; (b) 行星轮; (c) 齿圈 Fig. 3 Localized gear damage: (a) sun gear; (b) planet gear; (c) ring gear ·772·

赵川等：基于多域流形的行星齿轮箱局部故障识别 ·773· 常、太阳轮损伤、行星轮损伤、齿圈损伤)训练样本为 局部故障、行星轮局部故障、齿圈局部故障等运行状态 200个，测试样本为200个.行星齿轮箱正常、太阳轮 对应的采集信号时域波形如图4所示. 3 (a) 日 10 15 0 15 0 时间/s 时间s 0.5 d -1 0.5 10 15 20 10 15 20 时间/s 时间s 图4行星齿轮箱振动信号波形.(a)正常；(b)太阳轮损伤：(c)行星轮损伤；(d)齿圈损伤 Fig.4 Waveforms of vibration signals:(a)normal;(b)sun gear fault;(c)planet gear fault;(d)ring gear fault 3.2信号分析 理.可根据文献[10]直接计算12个时域统计量与4 3.2.1提取多域特征 个频域统计量.图5和图6分别给出了四种工况下 提取多域特征需要对每个振动条件下所采集的 (正常、太阳轮损伤、行星轮损伤、齿圈损伤)样本对应 50个样本中的每个样本进行时域、频域特征量及单分 的六个时域量纲统计特征值及六个时域量纲一统计特 量瞬时幅值能量的计算.由于不同域的特征量的量纲 征值.图7给出了相同样本对应的四个频域统计特 不同，为统一数据，需要将样本数据进行量纲归一化处 征值 1.0 (a 15 0.8 iFiwl.oowwwwwmi 0.6 whw 盖人 0.4 0. 50 100 150 200 0 50 10 150 200 样本点 样本点 0.8 c 0.15 nhyh-o-kmthyi 0.10 (d) 0.05 0.4--W 0 www- 0.2 0.05 0 50 100 150 200 50 100 150 200 样本点 样本点 3 200 (e) (f) 50 10 150 200 50 100 150 200 样本点 样本点 图5六个量纲时域特征.(a)绝对均值：(b)均方根植：(c)方根幅值：(d)标准差：(e)歪度：()蜂态 Fig.5 Six dimensional time-domain features:(a)mean;(b)root mean square;(c)root;(d)standard deviation;(e)skewness;(f)kurtosis 以一个正常工况信号样本为例进行局部均值分 同一状态下的样本在相同的样本空间中具有相似的分 解，所得5个P℉分量如图8所示 布属性和几何结构1o].本文应用SOMAP对行星齿轮 采用相同方法对其他工况典型样本进行分解，可 箱不同状态的训练样本多域高维特征空间进行降维 最终确定P℉分量个数为5.分别计算四种工况样本集 图10给出了特征空间的本征维数估计及k近邻点选 每个样本5个P℉分量的瞬时幅值能量，结果如图9 择结果 所示 由图10可知，(a)、(d)分别为k值过小与过大的 3.2.2多域流形学习 临界情况，考虑(c)的H,(d)值大于1，因而确定本征 每种运行条件下选取50个样本点作为训练样本 维数为2，最优k值为31，降维结果记为4mx2:定义 构建多域高维特征矩阵Pmx21·根据流形学习理论， 矩阵4ox2的第一列为维度V,第二列为维度V2

赵 川等: 基于多域流形的行星齿轮箱局部故障识别 常、太阳轮损伤、行星轮损伤、齿圈损伤) 训练样本为 200 个,测试样本为 200 个. 行星齿轮箱正常、太阳轮 局部故障、行星轮局部故障、齿圈局部故障等运行状态 对应的采集信号时域波形如图 4 所示. 图 4 行星齿轮箱振动信号波形 郾 (a) 正常; (b) 太阳轮损伤; (c) 行星轮损伤; (d) 齿圈损伤 Fig. 4 Waveforms of vibration signals: (a) normal; (b) sun gear fault; (c) planet gear fault; (d) ring gear fault 3郾 2 信号分析 3郾 2郾 1 提取多域特征 提取多域特征需要对每个振动条件下所采集的 50 个样本中的每个样本进行时域、频域特征量及单分 量瞬时幅值能量的计算. 由于不同域的特征量的量纲 不同,为统一数据,需要将样本数据进行量纲归一化处 理. 可根据文献[10]直接计算 12 个时域统计量与 4 个频域统计量. 图 5 和图 6 分别给出了四种工况下 (正常、太阳轮损伤、行星轮损伤、齿圈损伤)样本对应 的六个时域量纲统计特征值及六个时域量纲一统计特 征值. 图 7 给出了相同样本对应的四个频域统计特 征值. 图 5 六个量纲时域特征 郾 (a) 绝对均值; (b) 均方根植; (c) 方根幅值; (d) 标准差; (e) 歪度; (f) 峰态 Fig. 5 Six dimensional time鄄domain features: (a) mean; (b) root mean square; (c) root; (d) standard deviation; (e) skewness; (f) kurtosis 以一个正常工况信号样本为例进行局部均值分 解,所得 5 个 PF 分量如图 8 所示. 采用相同方法对其他工况典型样本进行分解,可 最终确定 PF 分量个数为 5. 分别计算四种工况样本集 每个样本 5 个 PF 分量的瞬时幅值能量,结果如图 9 所示. 3郾 2郾 2 多域流形学习 每种运行条件下选取 50 个样本点作为训练样本 构建多域高维特征矩阵 赘200 伊 21 . 根据流形学习理论, 同一状态下的样本在相同的样本空间中具有相似的分 布属性和几何结构[10] . 本文应用 ISOMAP 对行星齿轮 箱不同状态的训练样本多域高维特征空间进行降维. 图 10 给出了特征空间的本征维数估计及 k 近邻点选 择结果. 由图 10 可知,(a)、(d)分别为 k 值过小与过大的 临界情况,考虑(c)的 H1 ( d)值大于 1,因而确定本征 维数为 2,最优 k 值为 31,降维结果记为 驻200 伊 2 . 定义 矩阵 驻200 伊 2 的第一列为维度 V1 , 第二列为维度 V2 . ·773·