一种基于单通道信号盲分离滚动轴承的特征提取方法

摘要

本发明公开了一种基于单通道信号盲分离滚动轴承的特征提取方法，对原始振动信号进行频率切片小波变换，得到信号的能量谱，不再依靠小波基的选取，切片的选取也克服了小波变换频带选取受限的问题。采用主成分分析的方法确定信源数目，解决了当样本较大时获得聚类结论困难的问题。根据主成分分析采用投影方式得到降维的矢量投影矩阵，避免了信号稀疏性的影响，从而将欠定问题转化为适定问题，此外，对具有一定相关性的源信号也可有效实现故障特征提取。

权利要求

1.一种基于单通道信号盲分离滚动轴承的特征提取方法，其特征在于，具体包括以下步骤：步骤1：选取频率切片小波变换函数，对给定的原始振动信号进行频率切片小波变换，得到时频图，再对频率切片小波变换后的信号进行逆变换得到重构信号，并画出重构信号的能量谱；

步骤2：根据步骤1得到的时频图和能量谱，选择能量谱中包含能量峰值的多个区间对原始振动信号进行切片，得到多个切片的重构信号；

步骤3：将步骤2得到的重构信号进行降噪和去冗余，确定最佳信源数目m及其对应的m个特征向量；

步骤4：将m个特征向量组成的矩阵乘以步骤2得到的重构信号，得到降维后的m维矩阵；

步骤5：步骤4得到的m维矩阵，采用独立分量分析得到m个分离信号，并分别求其包络谱；对包络谱进行归一化处理，叠加得到等效包络谱；观测包络谱，提取故障特征；

所述步骤1的具体实现方法如下：

频率切片小波变换定义式为：

式中，σ为尺度因子，σ≠0，σ为常数或ω和t的函数， 或σ＝kt，其中，k为调节时域或频域灵敏度，t为时域变量，ω为频域变量；p(·)为频率切片函数，p*(·)是p(·)的共轭函数；f(τ)为原始振动信号，f(τ)∈L2(R)；

其中， C为时域振动信号的变量，

σ＝sqrt(2)/2/0.025；

对频率切片小波变换信号进行逆变换，得到频率切片小波变换信号W(t,ω,σ)在时频区域(t1,t2,ω1,ω2)的信号分量，即重构信号：能量定义公式为： 式中，E(ω)能量密度函数，E(ω)＝|F(ω)|2，其中F(ω)为重构信号fa(t)的傅里叶变换。

2.如权利要求1所述的基于单通道信号盲分离滚动轴承的特征提取方法，其特征在于，所述步骤3的具体实现方法如下：步骤3.1：求重构信号的协方差矩阵A＝(Aij)pxp，其中其中，Aij为协方差矩阵A中的元素，xik为xFSWT(t)第i行k列元素，为xFSWT(t)第i行平均值，n为数据长度；xjk为xFSWT(t)第j行k列元素， 为xFSWT(t)第j行平均值，xFSWT(t)为重构信号；

步骤3.2：计算上述协方差矩阵A的特征值λ1≥λ2≥λ3≥...λn＞0及其正交的单位化特征向量从左到右排列 按照公式 计算方差贡献率，当其值大于设定值时，确定信源数目m及其对应的m个特征向量。

说明书

一种基于单通道信号盲分离滚动轴承的特征提取方法

技术领域

[0001] 本发明涉及机械故障诊断技术领域，涉及一种基于单通道信号盲分离滚动轴承的特征提取方法。

背景技术

[0002] 盲信号分离是近些年来信号处理领域研究的热点问题。所谓盲信号分离就是指从若干传感器观测到多个信号的混合信号中恢复出无法直接观测到的原始信号的方法。盲信号分离一般要求传感器个数多于信号源个数，而单通道盲分离则是指用于观测信号混合信号的传感器只有一个。这是盲信号分离中的一个难点，但是这在机械故障诊断工程应用中却是更加贴近实际条件的。

[0003] 针对欠定盲分离即传感器数目小于信源数目的情况，主要的方法有基于中值的聚类算法盲分离方法[1]、基于位势函数的欠定盲分离方法[2]，这些存在的问题是基于源信号的稀疏问题，而对于稀疏性较差的信号则分离效果不好；还有基于小波分解的盲信号分[3]离 ，它存在的问题是对小波基的选择依靠性很强，小波基不同分离效果差别很大；此外也有基于固有模式函数(EMD)、总体固有模式函数(EEMD)的盲分离方法[4-5]，EMD方法存在的问题是存在模态混叠现象，EEMD方法问题是计算量较大，计算时间也较长，因此在实际工程应用中即时性不强。

[0004] 参考文献

[0005] 【1】FABIAN J T；CARLOS G P；ELMAR W L Median-based clustering for underdetermined blind signal processing[外文期刊]2006(02).

[0006] 【2】张赟；李本威；王永华；基于位势函数的欠定盲源分离识别诊断方法[J].航空动力学报，2010，25(01)，218-223.

[0007] 【3】王娇；刘郁林；何为；晁志超；小波分解单通道盲分离干扰抑制方法[J].[0008] 重庆邮电大学学报(自然科学版)，2014，26(5)，648-653.

[0009] 【4】李舜酩；刘晓伟；郭海东；一种单通道振动信号的盲分离方法[P].中国专利：CN 102288285 B，2011-05-24.

[0010] 【5】孟宗；蔡龙；基于EEMD子带提取相关机械振动信号单通道盲分离[J].振动与冲击，2014，33(20)，40-46.

发明内容

[0011] 针对上述现有技术中存在的问题和缺陷，本发明的目的在于，提供一种基于单通道信号盲分离滚动轴承的特征提取方法。

[0012] 为了实现上述目的，本发明采用如下技术方案：

[0013] 一种基于单通道信号盲分离滚动轴承的特征提取方法，具体包括以下步骤：

[0014] 步骤1：选取频率切片小波变换函数，对给定的原始振动信号进行频率切片小波变换，得到时频图，再对频率切片小波变换后的信号进行逆变换得到重构信号，并画出重构信号的能量谱；

[0015] 步骤2：根据步骤1得到的时频图和能量谱，选择能量谱中包含能量峰值的多个区间对原始振动信号进行切片，得到多个切片的重构信号；

[0016] 步骤3：将步骤2得到的重构信号进行降噪和去冗余，确定最佳信源数目m及其对应的m个特征向量；

[0017] 步骤4：将m个特征向量组成的矩阵乘以步骤2得到的重构信号，得到降维后的m维矩阵；

[0018] 步骤5：步骤4得到的m维矩阵，采用独立分量分析得到m个分离信号，并分别求其包络谱；对包络谱进行归一化处理，叠加得到等效包络谱；观测包络谱，提取故障特征。

[0019] 具体地，所述步骤1的具体实现方法如下：

[0020] 频率切片小波变换定义式为：

[0021]

[0022] 式中，σ为尺度因子，σ≠0，σ为常数或ω和t的函数， 或σ＝kt等，其中，k为调节时域或频域灵敏度，t为时域变量，ω为频域变量；p(·)为频率切片函数，p*(·)是p(·)的共轭函数；f(τ)为原始振动信号，f(τ)∈L2(R)；

[0023] 对频率切片小波变换信号进行逆变换，得到频率切片小波变换信号W(t,ω,σ)在时频区域(t1,t2,ω1,ω2)的信号分量，即重构信号：

[0024]

[0025] 能量定义公式为： 式中，E(ω)能量密度函数，E(ω)＝|F(ω)|2，其中F(ω)为重构信号fa(t)的傅里叶变换。

[0026] 具体地，所述步骤3的具体实现方法如下：

[0027] 步骤3.1：求重构信号的协方差矩阵A＝(Aij)pxp，其中

[0028]

[0029] 其中，Aij为协方差矩阵A中的元素，xik为xFSWT(t)第i行k列元素， 为xFSWT(t)第i行平均值，n为数据长度；xjk为xFSWT(t)第j行k列元素， 为xFSWT(t)第j行平均值；

[0030] 步骤3.2：计算上述协方差矩阵A的特征值λ1≥λ2≥λ3≥...λn＞0及其正交的单位化特征向量从左到右排列 按照公式计算方差贡献率，当其值大于设定值时，确定信源数目m及其对应的m个特征向量。

[0031] 与现有技术相比，本发明具有以下技术效果：

[0032] 1、本发明对原始振动信号进行频率切片小波变换，得到信号的能量谱，不再依靠小波基的选取，切片频带的选取也克服了小波变换频带受限的问题。

[0033] 2、采用主成分分析的方法确定信源数目，解决了当样本较大时获得聚类结论困难的问题。

[0034] 3、根据主成分分析采用投影方式得到降维的矢量投影矩阵，避免了信号稀疏性的影响，从而将欠定问题转化为适定问题。

[0035] 4、本发明对具有一定相关性的源信号也可有效实现故障特征提取。

[0036] 5、本发明与传统方法相比，计算量较小，计算速度较快。

附图说明

[0037] 图1是本发明的方法流程图；

[0038] 图2是实施例1的源信号及混合信号时域及频谱图；

[0039] 图3是实施例1的混合信号的时频图及重构信号能量谱图；

[0040] 图4是实施例1的分离后的源信号时域及包络谱波形图；

[0041] 图5是实施例1的归一化等效包络谱波形图；

[0042] 图6是实施例1的对比实验分离后的源信号时域图及包络谱图；

[0043] 图7是实施例1的对比实验分离后的源信号归一化等效包络谱波形图；

[0044] 图8是实施例2的观测信号时域波形图；

[0045] 图9是实施例2的观测信号时频图及重构信号的能量谱图；

[0046] 图10是实施例2的分离后的观测信号时域图及包络谱图；

[0047] 图11是实施例2的分离后的观测信号的归一化等效包络谱波形图。

[0048] 下面结合附图和实施例对本发明的方案做进一步详细地解释和说明。

具体实施方式

[0049] 遵从上述技术方案，参见图1，本发明的基于单通道信号盲分离滚动轴承的特征提取方法，具体包括以下步骤：

[0050] 步骤1：选取频率切片函数，对给定的原始振动信号进行频率切片小波变换，得到其0～fs/2频带内的时频图，再对频率切片小波变换后的信号进行逆变换得到重构信号，其中，fs为采样频率，并画出重构信号的能量谱。其具体实现方法如下：

[0051] 频率切片小波变换定义式为：

[0052]

[0053] 式中，σ为尺度因子，σ≠0，σ为常数或ω和t的函数， 或σ＝kt等，其中，k为调节时域或频域灵敏度，其取值为在时域和频域分辨率的折中选择，t为时域变量，ω为频域变量；p(·)为频率切片函数，p*(·)是p(·)的共轭函数；f(τ)为原始振动信号，f(τ)∈L2(R)，上述计算方式没有依靠对小波基的选取。

[0054] 由于信号时域和频域行为并非相互独立，因而频率切片小波变换结果是冗余的。理论上其逆变换可以有不同形式。对频率切片小波变换信号进行逆变换，得到频率切片小波变换信号W(t,ω,σ)在时频区域(t1,t2,ω1,ω2)的信号分量，即重构信号：

[0055]

[0056] 由上式可以看出，通过改变时频区域(t1,t2,ω1,ω2)，可以自由的在时频空间上提取所需分量，其中，t1,t2为时间值，ω1,ω2为频率值。

[0057] 能量定义公式为： 式中，E(ω)能量密度函数，即单位频率1Hz内的信号能量，E(ω)＝|F(ω)|2，其中F(ω)为重构信号fa(t)的傅里叶变换。

[0058] 若重构信号fa(t)为离散信号，则能量谱的绘制以一定频带宽度因子作为切片带宽，在全频带下按照该带宽大小，以一定步长的移动窗口，进行切片，带宽的选择为特征频率的1到5倍频，这样保证含有主要信息的切片在包络分析中能有较高的能量，步长选择越小，则精度越高，但是计算时间会很长，因此步长的选取应是精度和计算时间的折中选择。各个切片的相对比例能量，采用归一化相对能量监测，第m个切片的相对能量为:

[0059]

[0060] 式中：En(fa(t))为信号fa(t)的总能量， 为第m个切片， 为第m个切片的能量，En(m)为第m个切片的相对能量值。将所有切片相对能量值绘制出来则得到能量谱。

[0061] 步骤2：根据步骤1得到的时频图和能量谱，选择能量谱中包含能量峰值即能量分布较高的多个区间对原始振动信号进行切片，得到多个切片的重构信号。其具体实现方法如下：

[0062] 工程中采集信号传感器一般为低通滤波，而切片带的能量是由频率区间内各个频率成分能量组成，因此窗口移动时，由于频率区间变化，能量最大值会突然出现和消失，因此按照低通滤波的原则选择能量谱中上升快速且陡峭，即切片区间选取应为0到包含能量谱中的每一个峰值处的频率，该频率与邻近峰值之差不大于一个切片带宽，然后对原始信号进行切片，得到n个能量较高的切片的重构信号，并将其组成一个n维矩阵xFSWT(t)＝(c1,c2,...，cn)T，其中n≥3，且n为正整数。

[0063] 步骤3：将步骤2得到的重构信号，即n维矩阵xFSWT(t)＝(c1,c2,...，cn)T，采用基于主成分分析(PCA)进行降噪和去冗余，使保留下来的维度间的相关性尽可能小，同时其方差值尽可能大，从而确定最佳信源数目。与聚类分析法相比解决了样本较大时，获得聚类结论困难的问题。与奇异值分解法相比，主分量选用的是投影后的矢量，对原信号稀疏性要求更低。其具体实现方法如下：

[0064] 步骤3.1：求重构信号的协方差矩阵A＝(Aij)pxp，其中其中Aij为协方差矩阵A中的元素，xik为xFSWT(t)第i行k列元素， 为xFSWT(t)第i行平均值，n为数据长度；xjk为xFSWT(t)第j行k列元素， 为xFSWT(t)第j行平均值。

[0065] 步骤3.2：计算上述协方差矩阵A的特征值λ1≥λ2≥λ3≥...λn＞0及其正交的单位化特征向量从左到右排列 按照公式计算方差贡献率，当其值大于设定值时，则认为已覆盖原始信号的主要信息，从而确定信源数目m及其对应m个特征向量。

[0066] 步骤4：根据步骤3得到的信源数目m及其对应的m个特征向量，再计算m个特征向量组成的矩阵乘以重构信号即n维矩阵xFSWT(t)＝(c1,c2,...，cn)T，得到降维后的m维矩阵，即将n维矩阵在m维空间下进行投影，得到m维的矢量投影矩阵，从而将欠定问题转化为适定问题。

[0067] 步骤5：对步骤4得到的m维矩阵，采用独立分量分析(ICA),得到m个分离信号，并分别求其包络谱；对包络谱进行归一化处理，叠加得到等效包络谱，即观测包络谱，提取故障特征，实现单通道盲信号分离。其中FFj(t)为第j个信号的包络，peakj为其峰值。

[0068] 实施例1：

[0069] 本实施例的基于单通道信号盲分离滚动轴承的特征提取方法，具体包括以下步骤：

[0070] 步骤1：给定两个原始振动信号s1和s2，分别为：

[0071] s1＝cos(2πf1t+π/3)

[0072] s2＝cos(2πfbt)[1+βcos(2πfrt)]

[0073] 其中，f1＝25Hz，fr＝25Hz，fb＝115Hz，β＝2，采样点数为1024，采样频率fs为1000Hz，设混合信号模型s＝as1(t)+bs2(t)+n(t)，其中，a＝1，b＝1，n(t)为随机白噪声信号。仿真信号s1、s2和s的时域图和频谱图，如图2所示。

[0074] 选取频率切片函数，再对上述混合信号s进行频率切片小波变换，得到其0～fs/2频带内的时频图，再对频率切片小波变换后的信号进行逆变换得到重构信号，并画出重构信号的能量谱，如图3所示，其具体实现方法如下：

[0075] 频率切片小波变换定义式为： 针对离散信号，式中：σ为尺度因子，σ≠0，p(·)＝exp(-C2/2)，p*(·)＝(exp(-C2/2))*，p*(σ(τ-t))＝(exp(-(σ(τ-t))2/2))*；σ＝sqrt(2)/2/0.025。

[0076] 计算频率切片小波变换函数W(t,ω,σ)在时频区域(t1,t2,ω1,ω2)的信号分量，即重构信号：

[0077]

[0078] 其中，(t1,t2,ω1,ω2)为(0,1,0,500)；

[0079] 能量定义公式为： 式中，E(ω)能量密度函数，即单位频率1Hz内的信号能量。针对离散信号，能量定义为 n为信号长度，F(ω)为fa(t)的
傅里叶变换，Fi(ω)为F(ω)的第i个点，带宽为50Hz，为特征频率的2倍，步长为1Hz。

[0080] 步骤2：根据步骤1得到的时频图和能量谱，选择能量分布较高的多个区间对原始信号进行切片，得到多个切片的重构信号。其具体实现方法如下：

[0081] 参见图3，从整个能量谱中可看出，在30Hz，90Hz，120Hz，140Hz处能量升高快速，在其他处能量均几乎为0，因此按照低通滤波原则选择包含每一个能量分布较高部分进行切片，分别选择[0 40]、[0 100]、[0 130]、[0 150]，得到4个切片的重构信号，并将其组成一个4维矩阵xFSWT(t)＝(c1,c2,c3,c4)T，形成重构信号。作为对比试验，选择四个能量分布较低的切片进行实验分析，分别选择[40 80]、[100 120]、[160 200]、[300 400]，得到4个重构信号。

[0082] 步骤3：将步骤2得到的重构信号xFSWT(t)＝(c1,c2,c3,c4)T，采用基于主成分分析(PCA)的源数目估计方法，计算其特征值及特征向量，确定信源数目m。特征值如表1所示，信源数目为2。针对对比实验中的重构信号也采用基于主成分分析(PCA)的源数目估计方法，计算其特征值及特征向量，特征值如表2所示，信源数目为4。其具体实现方法如下：

[0083] 步骤3.1：求重构信号的协方差矩阵A＝(Aij)4x4，其中n为数据长度1024。

[0084] 步骤3.2：计算上述协方差矩阵A的特征值λ1≥λ2≥...λ4＞0及其正交的单位化特征向量从左到右排列 通过公式 计算方差贡献率，当其值大于0.90时，则认为已覆盖原始信号全部信息，从而确定信源数目为2。

[0085] 表1 多维信号xFSWT(t)特征值

[0086]

[0087] 表2 多维信号xFSWT(t)特征值

[0088]

[0089] 步骤4：根据步骤3得到的信源数目m＝2，再计算前2个特征向量组成的矩阵乘以重T构信号xFSWT(t)＝(c1,c2,c3,c4) ，得到降维后的2维矩阵，从而将欠定问题转化为适定问题。
将对比实验中的源信号阵也在其特征向量构成的空间下矢量投影，组成4维矩阵。

[0090] 步骤5：对步骤4中由矢量投影组成的新信号阵，即2维矩阵，采用独立分量分析(ICA),得到2个分离信号，并分别求其包络谱，如图4所示；对包络谱进行归一化处理，叠加得到等效包络谱，如图5所示，即 观测包络谱，提取故障特征，实现单通道盲信号分离。从图2中可以看出源信号特征频率为25Hz，从图5中看到分离后信号在25Hz处出现较高峰值，这与源信号特征频率吻合，图4中分离后时域信号及其包络谱与图
2中源信号及其包络谱也极其相似，比较其相关系数为0.9940、0.8061，相似度强。因此该盲分离方法能够实现特征频率提取。再将对比实验中得到的4维信号阵采用独立分量分析，得到4个分离信号，并分别求其包络，如图6所示，可以看出分离后信号杂乱无章，与图2中源信号及其包络完全不相似，再对包络谱进行归一化处理，叠加得到等效包络谱，如图7所示，无法观察出源信号特征频率，说明了采用这种方法不可行。

[0091] 实施例2

[0092] 利用本发明的方法对某一机械轴承故障信号进行盲分离，并提取故障特征频率，该轴承故障为轴承滚动体出现损伤。在该轴承驱动端安装振动传感器，采样频率fs为12K Hz，该设备加载有1HP负载，其转速为1777r/min,即其基频为29.6Hz，根据其部件特征频率系数算出滚动体故障特征频率为118.1Hz。在轴承运转时，轴承滚珠与内圈、外圈相互作用，其源信号具有一定相关性。

[0093] 步骤1：通过传感器测得单路观测信号，其时域波形图如图8所示。

[0094] 选取频率切片函数，再对该路信号进行频率切片小波变换，得到其0～fs/2频带内的时频图，再对频率切片小波变换后的信号进行逆变换得到重构信号，并画出重构信号的能量谱，如图9所示，其具体实现方法如下：

[0095] 频率切片小波变换定义式为： 式中：σ为尺度因子，σ≠0，p(·)＝exp(-C2/2)，p*(·)＝(exp(-C2/2))*，p*(σ(τ-t))＝(exp(-(σ(τ-t)
2 *
) /2)) ；σ＝sqrt(2)/2/0.025。

[0096] 计算频率切片小波变换函数W(t,ω,σ)在时频区域(t1,t2,ω1,ω2)的信号分量，即重构信号：

[0097]

[0098] 其中，(t1,t2,ω1,ω2)为(0,1,0,6000)；

[0099] 针对离散信号则能量定义为 n为信号长度，F(ω)为fa(t)的傅里叶变换，Fi(ω)为F(ω)的第i个点，n为信号长度8192。带宽为300Hz，为特征频率3倍左右，步长为1Hz。

[0100] 步骤2：根据步骤1得到的时频图和能量谱，选择能量分布较高的多个区间对原始信号进行切片，得到多个切片的重构信号。其具体实现方法如下：

[0101] 参见图9，从整个能量谱中可看出，在300Hz，1000Hz，2600Hz，3200Hz、3500Hz及4200Hz附近能量升高快速，而其他处相比之下能量极低，因此按照低通滤波原则选择包含每一个能量分布较高部分进行切片，分别选择[0 500]、[0 1200]、[0 2800]、[0 3400]、[0 
3700]及[0 4400]进行切片，得到其重构信号，组合成新的多维信号xFSWT(t)＝(x1,c1,c2,...，c6)T。

[0102] 步骤3：将步骤2得到的重构信号xFSWT(t)＝(x1,c1,c2,...，c6)T，采用基于主成分分析(PCA)的源数目估计方法，计算其特征值及特征向量，确定信源数目m。特征值如表3所示，信源数目为4。其具体实现方法如下：

[0103] 步骤3.1：求重构信号的协方差矩阵A＝(Aij)6x6，其中n为数据长度8192。

[0104] 步骤3.2：计算上述协方差矩阵A的特征值λ1≥λ2≥...λ6＞0及其正交的单位化特征向量 通过公式 计算方差贡献率，当其值大于0.90时，则认为已覆盖源信号全部信息，计算其方差贡献率，前四项累加和达到98.5％，从而确定信源数目为4。

[0105] 步骤4：根据步骤3得到的信源数目m＝4，再计算前4个特征向量组成的矩阵乘以重构信号xFSWT(t)＝(x1,c1,c2,...，c6)T，得到降维后的4维矩阵，从而将欠定问题转化为适定问题。

[0106] 表3 多维信号xFSWT(t)特征值

[0107]

[0108] 步骤5：对步骤4有矢量组成的新信号阵，即4维矩阵，采用独立分量分析(ICA),得到4个分离信号，并分别求其包络谱，如图10所示；对包络谱进行归一化处理，叠加得到等效包络谱，如图11所示，即 观测包络谱，提取故障特征，实现单通道盲信号分离。

[0109] 从图11中可以看出可以看到30Hz、59Hz、118Hz处出现较大峰值，其中30Hz与基频29.6Hz比较接近，59Hz近似为基频的二倍频，118Hz与滚动体故障特征频率为118.1Hz极为接近，由此判断是轴承滚动体出现故障，对滚珠进行更换后，再次进行分析，发现118Hz点峰值消失，说明了分析的准确性，也说明了该方法对有一定相关性的源信号可以实现故障特征提取。