一种用于轴承故障诊断的特征选取方法
【技术领域】
[0001] 本发明涉及一种机械故障诊断方法，尤其是一种用于轴承故障诊断的特征选取方 法。
【背景技术】
[0002] 轴承是机械系统的重要部件之一，轴承故障诊断越来越受到重视，已成为保证生 产安全稳定运行和提高产品质量的重要手段和关键技术。当前轴承故障特征提取方法大多 依赖于某种（类）故障的先验知识，专门针对一种或几种故障，而由于故障类型、工况环境 以及故障成因的多样性，其故障特征也各有差异，故面对复杂多变的故障环境，需要多种特 征提取方法配合才有所成效，造成计算量与资源浪费，其根本原因是缺乏通用的，不随故障 转移的特征提取方法。
[0003] 同时轴承故障特征提取方法的偶然性导致了特征选择的不确定性，甚至无法实施 特征选择过程，目前特征选择方法大多基于后续诊断的正确率效果反馈，以剔除弱故障区 分能力的特征（或选择强故障区分能力的特征），并重新建模。故特征对诊断效果先期估计 与评判能力不足，且反复建模使得诊断时间代价较大，模型难以适应动态变化的参数和环 境变化，缺乏诊断实时性。
【发明内容】
[0004] 为了克服现有故障诊断特征提取和选择方法的不足，本发明提供了一种用于轴承 故障诊断的特征选取方法。所述方法利用奇异值分解（SVD)进行故障信号特征提取，并定 义特征区分度公式进行特征选择，具体采用以下技术方案：
[0005] -种用于轴承故障诊断的特征选取方法，包括如下步骤：
[0006] 步骤（1):采集轴承不同故障类型的多个振动信号；所述不同故障类型包括不同 故障直径和不同故障部位；根据所采集的振动信号构建各个振动信号自相关函数的矩阵， 并对所构建的矩阵进行奇异值分解，得到最大奇异值元素；
[0007] 步骤（2):逐步增加矩阵维数，直至矩阵的最大奇异值元素保持稳定；得到此时各 个振动信号的奇异值向量和矩阵维数值；
[0008] 步骤（3):根据步骤（2)得到的各个振动信号的奇异值向量和矩阵维数值得到不 同维数下各个振动信号的奇异值向量的平均曲率；根据各个振动信号的奇异值向量的平均 曲率定义反映不同类型的轴承振动信号差异的区分度公式：
[0010] 其中，dab为第a类故障类型和第b类故障类型的特征区分度（a e {1，2，···，γ}; b e {1，2,…，r})，r为故障类型总数；kJPkb为这两类故障特征向量曲线的平均曲率；q为 分类性能指标参数统计值；
[0011] 步骤（4):逐步减小矩阵维数，计算矩阵奇异值向量，并根据步骤（3)定义的区分 度公式计算不同类型的轴承振动信号之间的区分度；当所有轴承振动信号之间的区分度最 小值为1时，此时的矩阵维数即为特征选取维数，同时得到的该维数下的奇异值向量即为 选取的特征向量。
[0012] 所述步骤（1)采集的振动信号：
[0013] yid= {x 丄（i = 1，2,…，k)
[0014] 其中，yld为第id (id = 1,2, 3···，r)类故障的轴承振动信号序列；r为总故障类别 数；Xi= x(t J为h时刻振动加速度数据；t ;= t o+Q-l) Δ t(i = 1，2,…，k) ;t。为初始采 样时刻，△ t = 1/f为相邻采样时间间隔；f为采样频率；k为总采样数；i为采样计数。
[0015] 所述步骤（2)中，逐步增加矩阵维数，直至矩阵的最大奇异值元素保持稳定过程 如下；
[0016] 步骤（21):周期扩展振动信号yld，并计算扩展信号的自相关函数
η为扩展后的振动信号第η个序列点，m为 计算自相关函数时的平移点数；选取较小的维数构造初始自相关函数Hankel矩阵C。：
[0018] 其中，C。为初始自相关函数Hankel矩阵，M。〉。，％>0为二维矩阵的初始选取维数；
[0019] 步骤（22)，对C。进行奇异值分解：C。= U1AV1/;其中，矩阵 ^jQ - I U^,U2,......G 。为 Μ。XM。维的正父矩阵，w」，，.:，为矩阵兀素向 量；1/?二(??e兵她_为NqXNq维的正交矩阵，1??:·、, Vwq为矩阵元素向 量；D。为对角阵：
[0021] 其中，H " "，％为自相关矩阵C。的奇异值，且4 > A2 > ··· >心〇，故
[0022] 步骤（23)，令Mt= M〇+t，Ns= N〇+s，（t = 1，2, 3..·，s = 1，2, 3..·)，Qp= min{Mt，Ns}， p = min{t, s};其中，t与s为矩阵横向和矩阵纵向的扩展维数，每次交替增加 I进行扩 展；重复步骤（21)和步骤（22)，直到得到最大奇异值元素 A1达到稳定；此时奇异值向量 /二（H.……，A(?t)为提取的特征向量。
[0023] 所述步骤（23)中得到最大奇异值元素 λ ^勺条件为若维数Q p时最大奇异值元素 等于维数Qp+1时的最大奇异值元素，或两者之差小于〇. 01 λ ^直，此时Q p的最大奇异值元素 为稳定的最大奇异值元素；其中Qp是指扩展后的维数；Q P+1是指Q p再次扩展后的维数。
[0024] 本发明的有益效果是：本发明提出了一种用于轴承故障诊断的特征提取和选择方 法，所述方法利用奇异值分解（SVD)进行故障信号特征提取，能够有效克服依赖故障特征 先验知识的局限，并提升轴承故障诊断通用性；所述方法定义特征区分度公式，能够有效评 判方法的故障诊断能力；特征选取能够有效提升后续诊断正确率，平衡诊断环境适应性与 实时性要求。本发明为轴承故障诊断提供了低维的，有明确区分的特征，在变工况负载等轴 承故障诊断应用领域具有潜在经济和社会价值。
【附图说明】
[0025] 图1为本发明的用于轴承故障诊断的特征提取和选择算法流程图；
[0026] 图2为本发明的目标振动数据与辅助振动数采集示意图（a :故障部位；b :故障直 径）；
[0027] 图3为本发明的轴承振动信号示例图（健康轴承，12kHz，驱动端对象）；
[0028] 图4为本发明的最大特征值与自相关矩阵维数关系曲线（故障类型：轻微，严 重）；
[0029] 图5为本发明的SVD提取奇异值向量的对数曲线图（a :故障直径；b :故障部位）；
[0030] 图6为本发明的奇异值特征向量的平均曲率曲线（a :电机转速；b :故障直径；c : 故障部位；d :试验对象）；
[0031] 图7为本发明的区分度与分类性能指标参数比值曲线（故障直径）；
[0032] 图8为小波变换特征提取与选择方法的特征频率的信号能量曲线（a :145. 83Hz ; b :157. 94Hz ;c :315. 89Hz ;d :473. 83Hz)。
【具体实施方式】
[0033] 下面结合附图及具体实施例对本发明进行详细说明。
[0034] 参见图1，一种用于轴承故障诊断的特征提取和选择方法，包括如下步骤：
[0035] 步骤（1):对含不同故障类型的轴承进行振动信号（加速度）采集，用于特征提 取、选择以及故障区分性能测试；
[0036] 步骤（1)中，不同故障类型，包括不同故障直径与不同故障部位的信号采集，其中 轴承振动信号采集示意图如图2所示：
[0037] a)例如，故障直径包括：0·000_、0· 178mm、0. 356mm、0. 533mm直径，分别定义为健 康、轻微故障、中度故障、严重故障；
[0038] b)例如，故障部位包括：内圈、滚珠、外圈6:00、外圈3:00、外圈12:00(负载6:00 方向）；
[0039] c)实测中，采集频率f = 12kHz或48kHz条件的振动信号，在本发明中，采集频率 为轴承转速的200倍以上。
[0040] Yid= {x il, (i = 1, 2, ···, k)
[0041] 其中，yld为第id(id= 1,2,3···,γ)类故障的含观测噪音轴承振动数据序列；r为 总故障类别数；Xi= x(t J为1^时刻振动加速度数据；t ;= t o+Q-l) Δ t (i = 1, 2,…，k); t。为初始采样时刻，△ t = Ι/f为相邻采样时间间隔；k为总采样数；i为采样计数。
[0042] 轴承振动信号示例图如图3所示，可以看出，在变化工作环境下，振动信号呈现周 期非平稳特性，并且混有高斯噪声，使得轴承的故障特征易被正常振动引起干扰信号以及 白噪声掩盖，增大提取轴承故障特征提取的复杂性和困难度。
[0043] 步骤（2):利用Hankel矩阵构建振动信号自相关函数的自相关矩阵C，增大矩阵维 数，观察其奇异值分解（SVD)得到的最大特征值元素变化，直至保持恒定，并得到稳定时的 振动信号特征向量（奇异值向量）与矩阵维数值；
[0044] 步骤（2)中，由于轴承故障发生时，损伤点与其它表面元件碰撞产生信号冲击，可 能覆盖轴承系统和传感器的固有频率，故自相关矩阵奇异值分解的实施目的是，克服大多 轴承故障诊断方法依赖故障特征先验知识（故障频率、噪音干扰等）的局限，提供不依赖先 验知识，并具备物理意义的清晰特征，提升轴承故障诊断通用性。
[0045] 步骤⑵中，自相关矩阵C的奇异值分解步骤如下：
[0046] 步骤（21)，周期扩展振动信号yld，并计算扩展信号的自相关函数
η为扩展后的振动信号第η个序列点，m为 计算自相关函数时的平移点数；选取较小的维数构造初始自相关函数Hankel矩阵C。：
[0048] 其中，C。为初始自相关函数Hankel矩阵，M。〉。，％>0为二维矩阵的初始选取维数；
[0049] 步骤（22)，对C。进行奇异值分解：C。= UQDQVQT;其中，矩阵
为Μ。X M。维的正父矩阵，?%; ，.4?..为矩阵兀素向 量；
为Ν〇ΧΝ。维的正交矩阵，…，％。为矩阵元素向 量；D。为对角阵：
[0051 ] 其中，......，心。为自相关矩阵C。的奇异值，且4 > I2 > .1(?、，故
[0052] 步骤（23)，令Mt= M〇+t，Ns= N〇+s，（t = 1，2, 3...，s = 1，2, 3...)，Qp= min{Mt，Ns}， p = min {t, s};利用MjP N s维数重复步骤B-I和B-2,得到最大奇异值元素 λ i，若维数Qp时最大奇异值元素等于（