结合连续小波变换和多约束非负矩阵分解的故障特征提取方法

唐曦凌，梁 霖，高慧中，罗爱玲

(1.广东电网公司电力科学研究院，广州 510080;2.西安交通大学 机械工程学院，西安 710049)

在机械故障诊断中，能否准确地提取故障特征是诊断的关键所在。由于机械振动信号常表现为复杂的非平稳信号，且早期故障特征往往会被大量噪声淹没，这为故障特征提取带来了很大的困难。

小波变换是目前分析非平稳信号的一种有效工具，而尺度连续变化的连续小波变换能够细致刻划信号的局部形态，采用与特征波形类似的小波基对信号进行连续小波变换，可通过时频分布图观测到故障特征的周期性变化。但是由于机械设备结构复杂，采集的振动信号复杂，且含有大量噪声，仅通过时频分布很难提取早期故障特征，且需要先验知识，对人员经验依赖性强，容易误判。

研究表明，故障特征信息主要存在于变换系数的部分尺度中，这使得变换的系数矩阵具有一定的稀疏性，同时，噪声也包含在系数矩阵中。因此，有学者提出利用SVD对系数矩阵进行特征约简［1－2］，并提取故障特征。然而SVD只能提取系数矩阵中的线性特征，当系数矩阵具有非线性分布时，则不能很好获得其低维本质结构。

而非负矩阵分解(Nonnegative Matrix Factorization，NMF)为提取数据非线性特征提供了一个新思路。它是一种以非负性为约束的矩阵分解方法，能够实现非线性降维和基于局部的特征提取，且分解结果更稀疏。该方法已逐渐被应用于多个领域，如Lee和Seung将NMF用于机器人对外界环境的感知，用来提取有意义的局部特征［3］;Kawamoto运用NMF分解音乐信号的幅度谱，最终得到单个音调的特征［4］;蔡蕾等［5］结合SNMF和支持向量机，对机械运行状态的时频图进行识别，进而对机器运行状态进行判断。文献［6］针对机械故障信号，利用正交性约束NMF在频域实现故障信号特征提取。

现有研究成果表明，由于非负性是较弱的约束性，这使得NMF分解结果常常并不唯一，因此需要根据实际应用中加入其他约束，从而获得满意的唯一解。

本文针对早期微弱故障的特征提取问题，根据故障特征在小波分解系数矩阵中的稀疏性以及信噪比较低的特点，将基于稀疏性和光滑性约束的多约束NMF(multiple constraints NMF，mcNMF)方法应用于连续小波变换(Continuous Wavelet Transform，CWT)的系数矩阵的分析中，提出了结合连续小波变换和多约束非负矩阵分解的特征提取方法CWT－mcNMF。通过仿真数据和实际案例表明该方法的有效性。

1 连续小波变换

设x(t)是具有有限能量的信号，则连续小波变换定义为x(t)与小波基函数的内积［7］，即:

式中:Wx(τ，a)为小波变换系数矩阵;τ为时间定位参数;a 为尺度参数，a ＞0;ψ(t)为小波函数;ψτ，a(t)为小波基函数，是由小波母函数ψ(t)进行伸缩和平移而形成的一组函数系列，即:

连续小波变换通过时间和尺度的连续变化将信号投影到时间—尺度平面上，系数矩阵Wx(τ，a)度量了信号与小波间的相似程度，反映了信号的特征信息。

2 非负矩阵分解方法原理

对一个T维的随机向量v进行了t次观测，得到观测矩阵 V=［v1，v2，…，vI］T，非负矩阵分解的问题可描述为:已知非负矩阵T，找出非负的W和非负的H，使式(3)成立［8］。

式中:V称为观测矩阵，是I×T的非负矩阵;W称为基矩阵，是I×J的非负矩阵;H称为系数矩阵，是J×T的非负矩阵。通常情况下I≫J，这表示可对观测矩阵进行维数约简，H则表示特征约简后的矩阵。

为了解决该问题，目前的算法中，标准的ALS算法较其他算法，效率更高，更适合大数据的处理。但该方法对噪声敏感，常获得次优解。而利用被处理信号的先验知识，对W和H施加更多的限制条件，如稀疏性、光滑性、正交性等，可获得较满意的解。加入相应的约束性后，非负矩阵分解最优化问题变为［9］:

式中:·F—矩阵的 Frobenius范数;αW为正则化系数，αW≥0;JW为对W施加的约束条件;αH为正则化系数，αH≥0;JH为对H施加的约束条件。

3 基于最优小波系数的多约束NMF特征提取

3.1 特征提取思路

由于通过连续小波变换，使信号的故障特征信息的变换系数与其他分量不同，并具有稀疏性，另外，变换系数中仍包含噪声分量，因此在利用非负矩阵分解对系数矩阵进行分析时，为了提取故障特征分量，可对结果，即H施加稀疏性和光滑性约束:施加稀疏性约束，可获得具有稀疏性的特征分量;施加光滑性约束，可提高特征分量的信噪比。最终实现系数矩阵的非线性降维;因缺乏W的先验知识，所以无法施加其他的约束性，本文采用以局部欧式距离为目标函数的HALS算法进行计算［10］。

小波变换中，为了最大程度的与故障特征匹配，需要根据不同的故障类型选择合理的小波基，变换后的系数矩阵才能准确的反映故障特征。在众多故障中，冲击类故障占有较高的比例，如缺陷轴承和齿轮产生的瞬时冲击等，因此本文针对冲击类故障，采用CWT－mcNMF方法提取微弱故障特征。

3.2 算法步骤

具体步骤如下:

步骤1:最优系数矩阵的获取。选择与冲击故障能够较好匹配的Morlet小波［11］，对于机械振动信号，常采用其实部作为小波基函数，即:

式中:β为波形参数，控制小波波形的衰减率。为了达到最优匹配，采用小波熵［12］确定波形参数β，并获得最优系数矩阵 Wβop(τ，a)。

步骤2:非负化。将系数矩阵每一行都加上一直流分量，使矩阵为非负。

步骤3:初始化。随机初始化W和H，输入约简维数J。

步骤4:分别迭代计算W和H，直到收敛。

(1)对H的迭代规则

对H加入稀疏性约束和平滑性约束的最优化问题如下:

式中:‖H‖2F=∑i‖hi‖22是对H施加光滑性约束，‖H‖l=∑jt|hjt|是对H施加稀疏性约束。

采用ALS获得H的迭代规则为:

式中:［x］+=max(ε，x)。

对带电粒子进行如下分析，带电粒子受到一个竖直向下的恒力电场力的作用，并且射入平行板电容器时带有一个水平的初速度v0。带电粒子在水平方向上做的是匀速直线运动，在竖直方向上做的是初速度为零的匀变速直线运动。这种运动类似于物体在重力作用下只受到一个水平初速度的平抛运动。所以可根据平抛运动的相关知识来解有关题目，可以求得以下有关物理量或相关结论：

(2)对W的迭代规则

采用HALS对W进行迭代，其获得的迭代规则为:

式中:V(j)=E+wjhTj，另外，每次迭代需对W的每一列进行归一化处理:wij←wij/‖wj‖2。

步骤5:获得特征约简结果H。

在上述算法步骤中，参数τ按采样周期离散，参数a的范围和步长根据特征频率和细化程度来选取;稀疏性正则化参数 αHsp∈［0.01，0.5］;光滑性正则化参数αHsm有 三 种 取 值 方 式［9］:常 数;αHcm(k)=20exp

4 仿真分析和实例分析

4.1 仿真分析

机械设备早期故障常常体现为微弱的冲击信号，常常被淹没在噪声或其他较大的振动中，为了验证本文特征提取方法的有效性，现模拟一带有噪声的周期性冲击衰减信号，信号波形如图1所示。

图1 含噪声冲击信号Fig.1 The impact noise simulation signal waveform

对该模拟信号采用CWT－mcNMF方法进行处理。相应参数为:尺度参数 a∈［1，10］，步长为 0.1，β =0.3(如图 2)，αHsp=0.09，αHsm，其中α0=0.001，J=2。考虑初始值对NMF算法的影响，重复运行算法多次，获得两个特征基向量，通过它们间的加权组合，能够构成原始的最优系数矩阵。其去均值后的波形如图3。

图2 冲击信号小波熵随波形参数的变化图Fig.2 Wavelet entropy change chart with wavelet shape factor

从图3可看出，第一个基向量含有包含冲击成分和幅值较大的噪声信息，而第二个基向量含有明显的冲击成分，说明被淹没的冲击特征被较好的提取出来。但是由于NMF属于矩阵分解问题，其结果存在幅值不确定性，导致两基向量的幅值与原始波形幅值相差较大。且由于NMF方法属于非线性特征提取方法，而在分解过程中并未施加正交性等约束，因此获得的基向量间具有冗余信息，在图中反映为在冲击脉冲较大的点，两个基向量的波形幅值都较大。但以上因素不影响故障特征的提取。

图3 CWT-mcNMF提取的微弱冲击信号Fig.3 The impact signal waveform extracted by CWT-mcNMF

4.2 实例分析

(1)轴承外环故障

现选取信噪比较低的轴承故障信号进行分析。实验对象为308型滚动轴承，模拟外环点蚀故障，采用加速度传感器采集振动信号，轴承转速为1 600 r/min，采样频率为40 k，采样点数为4 096，故障特征频率为82.05 Hz。信号时域波形图4，由图4可知，该故障信号信噪比较低，其轻微故障引起的周期性冲击信号被较大的噪声淹没，无法直接从时域观察到冲击周期分量。从低频区图5中也无法识别故障特征频率。

采用CWT－mcNMF提取波形中的周期性冲击成分。相应参数为:尺度参数 a∈［1，5］，步长为0.1;β =0.1(如图 6);αHsp=0.01;αHsm(k)=20exp(－ k/50);J=2。同样重复运行算法多次，获得两个特征基向量，去均值后波形如图7。

分析图7可知，第一个向量包含有冲击成分在内的轴承振动信息，第二个向量含有明显的冲击分量，且连续较大冲击间的时间间隔约为12.6 ms，对应特征频率为79.37 Hz，在此频率分辨率下，即为轴承故障特征频率。这表明CWT-mcNMF方法可有效提取微弱的冲击故障特征。同理，由于存在幅值不确定性，基向量与原始波形幅值相差较大;而第一个向量也存在微弱的冲击信息，说明NMF分解结果间具有冗余性。

图4 轴承外环故障原始波形Fig.4 The bearing outer ring fault signal waveform

图5 轴承外环故障原始波形频谱低频区Fig.5 The low-frequency region of the bearing outer ring fault signal frequency spectrum

图6 外环故障信号小波熵随波形参数的变化图Fig.6 Wavelet entropy change chart with wavelet shape factor

图7 CWT－mcNMF提取的轴承外环故障冲击信号Fig.7 The outer ring fault impact waveform extracted by CWT－mcNMF

(2)齿轮箱故障

为了进一步验证基于稀疏性和光滑性约束的非负矩阵分解提取冲击故障的能力，选择某二级传动的故障齿轮箱的输出轴信号进行分析。齿轮箱故障为中间轴的一级啮合齿轮断裂。电机转速约为1 500 r/min，采样频率为10 240 Hz，采样点数为4 096，从输入轴到输出轴间齿数分别为:26、64、40、85。经计算，输入轴转频z1=25 Hz，中间轴转频z2=10.15 Hz，输出轴 z3=5 Hz。故障特征频率为中间轴转频10.15 Hz。采集的输出轴加速度振动信号时域波形如图8。由于输出轴测点距断齿故障点较远，其冲击故障完全被淹没，而从其低频谱区图9中也无法直接找到故障特征频率。

采用CWT－mcNMF方法分析输出轴信号。相应参数为:尺度参数 a∈［1，10］，步长为 0.1，β =0.2(如图10)，αHsp=0.09;αHsm(k)=20exp(－ k/50);J=2，获得两个特征基向量去均值后的波形，如图11所示。

分析图11可知，在两个特征基向量中，第一个向量包含该输出轴的啮合等振动信息，而第二个向量具有明显的冲击分量，说明齿轮箱存在故障，但是仍无法直接判断其故障类型，对其进行FFT变换得到提取的故障波形频谱图，如图12所示，频谱图的低频区中非常清晰的显示10 Hz的故障特征频率，与齿轮箱中间轴一级啮合的齿轮断齿的故障吻合。

图8 齿轮箱断齿输出轴原始波形Fig.8 The output shaft fault signal waveform of gearbox tooth breaking

图9 齿轮箱断齿输出轴原始波形低频谱区Fig.9 The low-frequency region of theoutput shaftsignal frequency spectrum

从以上分析可知:① 由于非负矩阵分解实质上与独立分量分析相同，属于矩阵分解问题，因此得到的特征基向量具有幅值不确定性和顺序不确定性，但并不影响故障特征的提取;② 由于NMF属于非线性分解，且在分解过程中未施加正交性等约束，因此特征基向量间有冗余信息;③ 上述方法不仅适用于轴承和齿轮箱的故障信号，同样适用于其他机械振动信号。

图10 齿轮箱输出轴小波熵随波形参数的变化图Fig.10 Wavelet entropy change chart with wavelet shape factor

图11 CWT-mcNMF提取的齿轮箱输出轴故障冲击信号Fig.11 The tooth breaking fault impact waveform extracted by CWT-mcNMF

图12 CWT-mcNMF提取的故障波形频谱的低频区Fig.12 The low-frequency region of the gearbox tooth breaking fault signal frequency spectrum extracted by CWT-mcNMF

5 结论

本文针对小波变换在早期故障特征提取中存在的不足，提出结合连续小波变换和多约束非负矩阵分解的故障特征提取方法。该方法根据故障特征在系数矩阵中的表现特点，采用基于光滑性和稀疏性约束的非负矩阵分解方法对连续小波变换的最优系数矩阵进行非线性特征约简，获得信噪比较高的故障特征。仿真信号与实际的轴承和齿轮箱故障样本的应用表明，该方法能够获得基于局部的故障特征，实现非线性降维，能够在时域提取较明显的故障信息，很好的满足了故障诊断特征提取需求。

另外，非负矩阵分解方法是一种新兴的信号处理方法，在机械故障诊断的领域还比较少，因此，需要进一步研究适用于不同类型故障诊断的非负矩阵分解方法，和算法参数的控制，以及该方法与其他方法的结合等。另外，NMF的初始值对算法的收敛性也有较大影响，因此可结合初始化方法，以期得到最优结果。

［1］梁 霖，徐光华，侯成刚.基于奇异值分解的连续小波消噪方法［J］.西安交通大学学报，2004，38(9):904 －908.LIANG Lin，XU Guang-hua，HOU Cheng-gang.Continuous wavelet transform denoising method based on singualr value decomposition［J］.Journal of Xi’an Jiaotong University，2004，38(9):904 －908.

［2］Brenner M J.Nonstationary dynamics data analysis with wavelet-SVD filtering［J］.Mechanical Systems and Signal Processing，2003，17(4):765 －786.

［3］Lee D D，Seung H S.Learning in intelligent embedded systems［EB/OL］.http://hebb.mit.edu/people/seung/papers/lee_es99.pdf.1999.

［4］Kawamoto T，Hotta K，Mishima T，et al.Estimation of single tones from chord sounds using non-negative matrix factorization［J］.Neural Network World，2000，3:429－436.

［5］蔡 蕾，朱永生.基于稀疏性非负矩阵分解和支持向量机的时频图像识别［J］.自动化学报，2009，35(10):1272－1277.CAI Lei， ZHU Yong-sheng. Time-frequency spectra recognition based on sparse non-negative matrix factorization and support vector machine［J］.Acta Automatica Sinica，2009，35(10):1272－1277.

［6］栾美洁.基于非负矩阵分解的机械故障信号特征分析的研究［D］.南京，东南大学，2009.

［7］Mallat S.A theory for multiresolution signal decomposition:the wavelet representation［J］. IEEE Trans Pattern Analysisand Machine Intelligence，1989，11(7):674－693.

［8］Lee D，Seung H S.Learning the parts of objects by nonnegative matrix factorization［J］.Nature，1999，401:788－791.

［9］Cichocki A，Zdunek R，Phan A H，et al.Nonnogetive matrix and tensor factorizations［M］.United Kingdom:Wiley ＆Sons，Ltd，Publications，2009.

［10］Phan A H，Cichocki A.Multi-way nonnegative tensor factorization using fast hierarchical alternating leastsquares algorithm(HALS)［C］.Proc.of The 2008 International Symposium on NonlinearTheory and itsApplications，Budapest，Hungary，2008.

［11］栗茂林，梁 霖，王孙安，等.基于连续小波系数非线性流行学习的冲击特征提取方法［J］.振动与冲击，2012，31(1):106－111.LI Mao-lin，LIANG Lin，WANG Sun-an，et al.Mechanical impact feature extraction method based on nonlinearmanifold learning of continuous wavelet coefficients［J］.Journal of Vibration and Shock，2012，31(1):106 －111.

［12］林 京，屈梁生.基于连续小波变换的信号检测技术与故障诊断［J］.机械工程学报，2000，36(2):95－100.LING Jin，QU Liang-sheng.Feature detection and fault diagnosis based on continuous wavelet transform［J］.Journal of Mechanical Engineering，2000，36(2):95－100.