动态模式分解及其在轴承早期故障诊断中的应用

文 明， 党 章,3， 余 震， 吕 勇 魏国前

(1. 武汉科技大学 冶金装备与控制技术教育部重点实验室，武汉 430081；2. 武汉科技大学 机械传动与制造工程湖北省重点实验室，武汉 430081；3. 武汉科技大学 国家实验机械教学示范中心，武汉 430081)

作为机器中常见的零部件，轴承被广泛应用于各种机械设备中，大部分造成设备停机或者重大事故的滚动轴承在其服役早期通常含有不易发现的微弱故障[1]。对滚动轴承早期故障进行有效诊断不仅能够排除设备故障隐患，而且对于提高设备运行的经济性和可靠性意义重大[2]。传感器采集得到的原始系统动力学响应信号通常包含表征机器健康状况的特征信息，这些特征信息是滚动轴承健康状态判断和故障诊断的依据。早期故障轴承由于故障特征不明显，设备运行过程中激起的能量低，且信号在传递过程中与设备本征信号及外部环境噪声动态耦合，使得早期故障信号具有非线性、非平稳、故障特征微弱等特点,用传统的故障诊断方法难以提取[3]。

近年来，针对复杂非线性、非平稳信号，已经提出了各种有效的时频分析方法和模式分解方法，并应用于滚动轴承早期故障信号的故障诊断中[4]。短时傅里叶变换和小波分析方法基底具有构造简单、对称等特点，但其结构固定，缺乏变化，不适于分析复杂信号[5]。近十几年出现的经验模式分解(empirical mode decomposition, EMD)[6]、局部均值分解(local mean decomposition, LMD)[7]、变模式分解(variational mode decomposition, VMD)[8]及其变种模式分解方法[9-11]为信号分析带来了方便，但EMD和LMD缺乏严格的理论公式推导，且分解得到的本征模式函数中存在低秩成分和稀疏成分混叠现象。当轴承存在早期故障时，由于其故障成分能量不高，噪声相对于特征成分较强，其表征故障的低秩成分往往淹没在强背景噪声中，使得这些模式分解方法的有效性和效果受限，有时甚至不能提取出故障特征成分。VMD相对于EMD，LMD，其理论性更强，对噪声鲁棒性更好，但其分解效果依赖于参数的有效选取。这些信号分析方法在某些场景下具有良好的应用效果，但应用于强背景噪声的早期故障信号时，使用效果将不可避免受到影响。

相比上述模式分解方法，动态模式分解(dynamic mode decomposition, DMD)[12]作为一种基于奇异值分解(singular value decomposition, SVD)的无方程数据驱动频率分析方法，能够精确提取复杂动态系统的时空特征。其最早起源于流体动力学领域，近年来被广泛应用于流体场变量的特征提取、趋势预测以及复杂系统控制等方面[13]。DMD算法结合了SVD和多模式表征的优点，通过一个动态数据分析模型，能够将蕴含机械动力学特性的信号(数据)分解成一系列具有内在时空特征的单频非正交模式[14]，每个模式对应于原始信号的低秩成分(有用成分)或稀疏成分(噪声成分)，即使是具有强背景噪声的早期故障信号，也同样能将其分解为系列单频模式，避免了EMD和LMD方法中的模式混叠问题。DMD方法将提取出的单频模式与时间矩阵进行内积，得到的重构信号能够高概率地恢复出原始信号的动力学特性。基于数据驱动的DMD方法为机械故障诊断中的复杂信号分析提供了一种新的研究思路和手段。

1 动态模式分解与多尺度排列熵

1.1 动态模式分解

假设信号采集系统以等时间间隔(Δt=ti+1-ti)对原始连续系统进行采样,得到一维时间序列x=(x1,x2,…,xN)∈R1×N,以步长为1的定长窗口滑动选取x中元素并按行排列为Hankel矩阵，表示为

式中：m+n-1=N；X∈Rm×n。将X矩阵分解成两个连续汉克子矩阵Xt，Xt+1，将其表示为时间上前后的两个快照，表示为

Xt=[X1X2...Xn-1]
Xt+1=[X2X3...Xn]

(2)

估计一个最优线性算子A实现Xt到Xt+1的映射，使得

Xt+1=AXt

(3)

(4)

式中：U∈Rm×r,VT∈Rr×m分别为左右特征向量矩阵，且U×UT=I,V×VT=I;Δ为主对角线上含有不为零元素的矩阵，即特征值矩阵。通过最小化Xt+1与AXt之间的Frobenium范数误差求解最佳近似解,表示为

(5)

则最佳近似解可由式(5)解得

(6)

(7)

(8)

式中，φ为标准DMD模式。令ωi=ln(λi)/Δt，则原信号X的近似重构解为

(9)

(1)将原始信号X按照式(1),式(2)选取维度m进行相空间重构，得到式(2)中两个连续快照。

(4)按照式(3)～式(9)计算各模式和恢复信号。

1.2 截断秩选取与多尺度排列熵

1.2.1 截断秩的选取

DMD能够精确地将信号分解成多个代表系统低秩成分或稀疏成分的模式，但选择哪些模式用于恢复重构，即截断秩r的选取是个重点问题。截断秩r的大小决定了用于恢复的模式数量，如果r取值太大，可能将一些不必要的模式用于重构导致恢复信号包含大量噪声，不利于特征提取；如果r太小，那么恢复的信号将丢失代表早期故障特征的有用成分。许多学者曾提出并使用若干基于奇异值的硬阈值选取方法，如奇异值差分谱[16]、奇异值能量谱[17]及奇异值标准谱[18]等，但相比而言，根据信号自身特点的软阈值选取方法更加合理， Gavish等[19]提出了自适应截断秩r的选取方法。对于包含噪声和系统动力学信息的观测矩阵

Y=X+ηN

(10)

式中：N为独立分布的高斯白噪声;X为低秩矩阵，η为噪声水平系数。对于实际工程信号而言,信号噪声大小往往是未知的，Gavish等指出，其幅值可以根据Y的奇异值中值进行估计。Gavish等应用奇异值中值λmedian估计噪声矩阵ηN,以噪声矩阵Y的相关参数λ为阈值，小于阈值的特征值所对应的模式去除，大于阈值的模式保留。阈值可以定义为

λ=ω(β)λmedian

(11)

ω(β)≈0.56β3-0.95β2+1.82β+1.43

(12)

式中：β为Y的行与列之比;ω(β)为阈值系数计算公式;λmedian为奇异值中值;λ为最佳截断秩阈值。该文献提供了一种根据信号自身特点选取截断秩的方法，通过该方法确定截断秩并恢复信号后，信号频域可读性很高，但低能量特征频率(如内圈故障频率边频)容易被当作背景噪声去除。传统以奇异值能量作为硬阈值选取截断秩的方法往往不能准确的找到合适的阈值，有时为了避免部分特征频率被“丢弃”选择了过大的比值，使得故障识别过程中存在过多噪声干扰。为了克服这个不足，根据式(11)确定截断秩后，重新寻找一个自适应的奇异值能量阈值，把式(11)确定的奇异值能量与全部奇异值能量比值的均方根作为新的能量阈值，取得自适应能量截断秩，该方法能够更为准确的找到合理的比值，弥补了传统奇异值能量阈值选取不够精确的问题。假设通过式(11)得到的截断秩为r，则修正后的奇异值能量比值表示为

(13)

式中，p=min(2m,n-1)。

1.2.2 多尺度排列熵

排列熵[20]是一种衡量时间序列复杂度的算法,代表了系统的复杂程度，但其无法在多个时间尺度上进行数据复杂度的测量。多尺度排列熵(multi-scale permutation entropy, MPE)[21]是在排列熵基础上进行改进的一种反映信号不确定性和不规则性的非线性分析方法，具有更好的稳定性和更强的抗噪声能力[22]，当存在动态观测噪声时，这种方法很有用。MPE是一种在多个时间尺度上测量系统复杂度的方法，它具有算法简单，运算速度快等特点，其基本思想是将时间序列首先进行多尺度粗粒化，然后再进行排列熵值计算。熵值越小，表明系统时间序列趋于有序，呈现出周期性，如轴承、齿轮系统正常运转产生的周期振动信号，反之趋于无序，如机械系统中广泛存在的高斯白噪声。

对于一个给定的时间序列x={xk,k=1,2,…,N},首先对其进行粗粒化处理。

(14)

式中：s为尺度因子，s=1,2,…；当s=1时，粗粒化序列即为原始序列,MPE与排列熵值等价;[N/s]为对N/s往下取整；当s≥2时，原始序列被粗粒化长度为[N/s]的时间序列。

对向量中的每两个数据之间进行升序或降序考虑,每个重构向量有d!个可能的顺序模式。多尺度排列熵值H最终被定义为

(15)

式中，pi为某一个特定的顺序模式在所有重构向量中的频率，图1描述了针对一个时间序列pi的计算过程。正则参数ln(d!)保证了H(N,d,τ)的最大取值为1。

图1 MPE计算流程Fig.1 Calculation process for MPE

本文在改进的全局DMD基础上采用MPE算法计算DMD所得系列单频模式的复杂度，通过设定阈值，将大于该阈值的模式作为噪声过滤，小于阈值的模式认为包含有原始系统的动力学特征保留，并将其用于信号重构，以进行信号特征提取。本文技术路线如图2所示。

图2 基于MPE的DMD技术流程图Fig.2 Flow chat of DMD based on MPE

2 算法在仿真信号中的应用

将本文提出的基于自适应截断秩和MPE的DMD应用于Randall等[23]提出的轴承仿真信号模型，该模型由方程式(16)描述，根据轴承工作特点，当轴承内圈相对旋转时，不同位置的转子将承受不同的交变载荷，当轴承出现点蚀，胶合等缺陷时，轴承将产生瞬时冲击信号，且振动信号以自身高频振荡衰减。

假设以等时间间距Δt进行信号采样，S(t)为轴承的自然衰减频率函数，Ak为第k个冲击响应的幅值，N(t)为均值为零的背景噪声，仿真信号模型表示为

(16)

式中：ak为第k个冲击能量幅值;γ和φA为初始相位；fm为调制频率；B为轴承系统相关的衰减系数；τk为波动引起的滞后时间；cA为一个常数，轴承内圈故障为fi，转频为fr，式(16)中的参数选取如表1所示。在轴承早期故障信号中，故障特征信息能量较低，容易被背景噪声淹没，因此仿真信号增加了-1 dB的高斯白噪声，采样频率设置为fs=3 000 Hz。

表1 轴承仿真信号参数选择Tab.1 Parameter selection for bearing simulation signal

试验中设置m=2 000,文献[24]建议采样点数N>5d!=3 600,因此取N=4 000。轴承内圈故障仿真信号时域频域如图3所示。

图3 仿真信号时域与频域图Fig.3 Simulation signal in time domain and frequency domain

从图3可以看出，信号频域存在尖峰干扰，故障频率需要经过初步识别才能找到，且其转频和特征频率的边频寻找困难，能量较小的边频淹没在背景噪声中。根据1.1节式(1)、式(2)构造汉克时间序列矩阵，按照步骤1～步骤3将Xt,Xt+1矩阵进行投影降噪，截断秩根据提出的自适应能量软阈值确定，使得信号在初步降噪时保留住特征频率而减少噪声的干扰，然后根据截断秩对系统特征算子A进行特征估计。为了进一步滤除其中的一些噪声，获得更加可辨识的故障信号，采用1.2.2节中的MPE对模式进行阈值筛选，熵大于阈值的模式认为包含较多噪声滤除，熵小于阈值的模式进行保留，通过最佳截断秩筛选，共有r=309个模式，然后对筛选出的模式进行MPE值计算。

MPE是系统复杂性的度量，其值大小与尺度因子s和嵌入维度d有着明显的关系。若d太小，重构向量中包含状态太少，算法失去意义和有效性。若d太大相空间重构会均匀化时间序列，不仅导致计算费时，也不能反映出时间序列的微小变化。根据文献[25]的建议，MPE中嵌入维数d取值为3～7，时延τ对时间序列的影响较小[26]，一般取τ=1。并据此经验选取d=6,s=5,将选取的模式进行熵计算，得到各模式的熵分布如图4所示。

图4 各模式MPE值Fig.4 MPE value of each mode

取熵阈值为0.85，将熵值小于0.85的模式用于信号重构，共得259个模式，并对恢复的信号进行频域观察。恢复信号频域如图5所示。

图5 恢复信号频域Fig.5 Frequency spectrum of recovered signal

通过MPE筛选后，恢复信号频域噪声相比于原信号少很多，虽然频域依旧存在若干不明干扰成分，但能量较强的特征频率和能量较弱的边频都能够找到，说明了基于MPE的DMD算法在早期故障信号特征提取中的可行性。为了表明该方法的故障特征提取效果，将该方法与EMD和SVD去噪效果进行对比。对信号采用EMD后，将含有特征频率的本征模式函数进行恢复得到经验模式恢复信号，SVD方法的截断秩通过式(11)确定，通过各方法恢复的信号频域如图6所示。

图6 不同方法处理后频谱图Fig.6 Frequency spectrum of signal processed by different methods

从各恢复信号频域看，EMD恢复频域依旧存在很多噪声，且特征频率边频辨识困难，没有很好的从背景噪声中分离。SVD虽然去噪效果很好，但是也去除了淹没在噪声中的边频，对信号进行过度去噪，导致丢失了有用成分。相比于这两种传统的去噪方法，MPE算法对模式进行阈值降噪，其恢复效果是这两种传统方法的折衷，通过阈值筛选过滤掉噪声模式，进一步去除了背景噪声的干扰，从恢复的频域看，特征频率及其边频能够全部找到，且噪声干扰更小，识别难度更低，准确度更高。通过正交投影降噪和MPE去噪，信号得到了良好的去噪效果，恢复信号特征频率辨识更为简单方便，综合识别的准确性和难易程度看，本文方法在强背景噪声下的轴承故障成分提取中具有一定的优越性。

3 算法在试验信号中的应用

实际工程信号比仿真信号更为复杂，信号成分更加多样，将本文方法应用于实际信号中，以验证其工程应用效果。该节采用了Cincinnati大学公布的轴承试验数据[27]。

Cincinnati大学试验台结构如图7所示，4个Rexnord ZA-2115型双列轴承安装在一个由交流电机带动的2 000 r/min的轴上，轴承每列有16个滚动体，轴承节距和滚动体直径分别为28.15 mm，3.31 mm，轴承接触角为15.17°，信号由装在轴承箱的高灵敏度石英ICP加速度计(PCB 353B33)收集，采用NI DAQ Card 6062E收集信号。

图7 Cincinnati 大学测试试验台Fig.7 Cincinnati University test bench

采用第二个数据集文件(2004.02.12.19.32.39)第一个信道进行验证，该信道存在轴承外圈故障，由于噪声的存在，故障频率淹没在背景噪声中。原信号的采样频率fs=20 kHz,数据点数为N=20 480个，为了降低计算时间成本，在原信号基础上按照等间距对信号进行重新采样，采样率降为一半，其频域几乎不变，因此设置新采样频率fs=10 kHz,N=10 240。由文献[28]中的计算公式求解转频fr=33.3 Hz,外圈故障为fo=236.17 Hz。m取5 000，图8为其信号时域和频域图。

图8 Cincinnati试验信号时频域Fig.8 Cincinnati experimental signal in time domain and frequency domain

从图8所示信号看，信号在低频处出现若干尖峰，为转频所在频段，信号特征频率被背景噪声淹没无法识别。对构造的汉克矩阵进行初步投影降噪后，通过自适应截断秩方法共获得752个模式。利用MPE对这些模式进行熵计算，各模式熵值如图9所示。

图9 各模式MPE值Fig.9 MPE value of each mode

设置阈值为0.85，将高于阈值的模式进行滤除，剩余得到289个模式。将这些模式进行重构恢复，得到恢复信号，其频域如图10所示。

图10 MPE恢复频域Fig.10 Recovered frequency spectrum with MPE

虽然恢复频域中仍然存在噪声，但是通过r阶投影去噪和MPE过滤后，相比于原信号频域，噪声已经大为减少，频域中仅存在少量无法辨认的高能频率，且故障频率的二倍谐频也被找到，故障频率及其谐频已经易于辨识。因此结合自适应截断秩和MPE的动态模式分解方法在轴承故障提取中具有一定的可行性及有效性。同样将EMD，SVD应用于该信号进行提取效果对比，各频域如图11所示。

图11 各方法处理后频域Fig.11 Frequency spectrum of signal processed by different methods

从图11可知，相比于原始频域，EMD虽然能够提取特征频率，但降噪效果不明显，背景噪声围绕在特征频率附近,没有得到有效去除。在早期故障信号的强噪声背景下，通过SVD方法恢复的频域图中，特征频率及其谐波被当作噪声去除，导致故障特征无法提取。MPE能够根据熵值筛选出不同的单频模式分量，熵值大于阈值的模式被视为噪声去除，低于阈值的模式被视为有用成分保留，在去除噪声同时，能够较为方便的找到特征频率。虽然MPE不能去除全部噪声，从特征识别结果中能够看出该方法具有优于传统去噪方法的优势，具备一定的应用前景。总体来说，本文提出的基于自适应截断秩和MPE的DMD方法能够更为有效的提取出含有早期强背景噪声故障信号中的故障特征频率，具有较强的应用前景。

4 结 论

在轴承早期故障信号中，噪声往往充斥在整个频域，由于早期故障频率能量较弱，噪声相对较强，轴承故障成分常常淹没于本底噪声之中而难以直接识别。为有效剔除早期轴承故障信号中的强背景噪声，准确提取出故障频率，本文提出选取改进的全局算子DMD方法将信号分解成系列单频模式，克服了传统模式分解的一些局限，为后续噪声滤除带来方便。进一步提出通过自适应最佳截断秩方法和MPE阈值法对原始系统模式成分进行过滤筛选，对含噪信号中的有序模式进行了高概率恢复。

本文通过对汉克时间序列矩阵进行r阶正交投影，构造了去噪汉克矩阵，在初步去噪的同时解决自适应截断秩选取的问题，使得信号在保留特征频率的同时包含更少的噪声。对系列表征动力学特征信息的模式分量使用MPE进行阈值筛选，进一步提取出含噪更少的模式。通过对筛选出的模式进行恢复，能够更加方便简单地提取出原信号的特征频率。仿真信号和试验信号的试验结果表明，本文提的基于自适应截断秩和MPE的DMD方法在强背景噪声下的故障信号特征提取中具有良好的应用潜力。

虽然通过本文方法能够较好地对早期轴承故障信号进行降噪并提取特征频率，但重构信号依旧存在若干不明成分，同时，在MPE阈值的选取上，也需要针对主观恢复效果进行人工选取。因此，基于典型工况下不同轴承信号模式分量的熵值信息值得进一步研究，以提出自适应熵阈值选取方法。

Vol．41 No．12 2022