高阶局部投影算法及其在滚动轴承故障诊断中的应用

吕 勇， 施 威， 易灿灿

(武汉科技大学 机械自动化学院 冶金装备及其控制教育部重点实验室，武汉 430081)

滚动轴承已经广泛地运用于各种旋转机械设备中，由于轴承等关键零部件发生故障导致整个设备停止运行的现象层出不穷，由此造成了巨大的经济损失。目前，由于一般很难知道机械系统具体的数学模型，因而对滚动轴承运行状态的监测主要是通过传感器测得反映其某一运行状态的振动信号进行的。这些振动信号本质上是由速度、加速度或位移等单一变量组成的一维时间序列。然而，测得的振动信号中含有很多与轴承运动状态无关的干扰信号，在信号的实际测量过程中，很容易受到环境噪声的干扰，最终获得的信号具有非平稳、非线性的特点，即呈现一定的混沌特性。由于在频域中混沌信号与噪声具有相似的宽频特征，因此传统的线性滤波方法具有很大的局限性。随着非线性科学的不断发展，基于相空间重构的降噪方法[1-5]逐渐成为研究非线性系统的重要工具之一。其中局部投影[6-9]是降噪效果比较好的一种方法，而且它不需要对系统模型有任何的先验知识，因而广泛地应用于各个领域[10-11]。

邻域质心的选取是局部投影降噪过程中十分关键的一步，它决定着降噪效果的好坏。Cawley等将邻域相点的平均值作为质心，这种方法由于局部线性化会产生较大的误差，对降噪效果影响较大。后来Sauer为了进一步地抑制噪声而利用二阶多项式对邻域质心进行了更精确地估计，使其与投影的超平面近似相切。最近，Moore等[12]提出对邻域质心进行更高阶地估计，提出了一种高阶局部投影降噪算法，并将其运用于医学信号处理中。虽然局部投影算法在机械故障诊断领域已有较为广泛的应用，但大多是针对齿轮信号而且都是从时域提取特征。本文提出将高阶局部投影单独地应用于滚动轴承的故障诊断中，首先通过仿真信号验证了该方法相比于传统的基于相空间重构的非线性时间序列降噪方法的优越性。最后利用高阶局部投影算法对从工业现场测得的风机轴承振动信号进行了降噪处理，降噪后从频域成功地提取出了风机轴承故障特征。

1 局部投影算法基本原理

一个从混沌系统测得的长度为N的时间序列s(1),s(2),….,s(N)∈R,选择合适的嵌入维数m和延迟时间τ就可以将其重构到m维的相空间中,则重构后相空间中的每个相点可以表示为

Xn=(sn,sn-τ,L,sn-(m-1)τ)

(1)

根据Takens[13]嵌入定理：当系统吸引子的分数维为d时，只需要满足最小嵌入维数m≥2d+1时，重构后的相空间与原来的动力学系统是微分同胚的，即具有相同的动力学特性。在m维相空间中，反映系统动力学行为的吸引子通常只局限在一个m0的低维子空间(d

步骤1选择合适的参数m和τ将含噪时间序列重构到m维相空间。

步骤2对于每个相点，确定其邻域Un。邻域的确定有固定邻域数目和固定邻域半径两种方式，本文采用的是前者。

步骤3计算每个邻域相点的质心。

步骤4计算每个相点邻域的协方差矩阵C并对其进行特征值分解，求出m-m0个较小特征值对应的特征向量aq，其中q=1,2,…,m-m0。

步骤5减去相点在噪声子空间上的投影，即

(2)

步骤6返回步骤2，直到所有数据处理完毕。

2 高阶局部投影算法

针对局部投影降噪步骤3，最初Cawley等将邻域相点的平均值

(3)

作为邻域质心，如图1(a)所示，这种方法由于局部线性化会产生较大的误差，各个邻域的中心与曲线都不是相切的，而向曲线内部发生了一定程度的偏移，对降噪效果影响较大。后来Sauer为了进一步地抑制噪声而利用二阶多项式

(4)

对邻域质心进行了估计，如图1(b)所示，本质上是通过将邻域中心向外移的方法使其与投影的超平面近似相切，本文将其称为标准局部投影算法。

(a)局部线性化(b)Sauer的改进方法

图1 质心估计方法示意图

Fig.1 Schematic diagram of centroid estimation methods

Moore等使用更高阶的多项式对邻域质心进行更准确地估计，并将其运用于医学信号处理中，取得了较好的效果。本文将其称为高阶局部投影算法，基本原理如下：

对于任意实数δ>0，定义连续移动平均算子Iδ

(5)

Sauer通过考虑将一阶和二阶的线性组合作为邻域的质心，提高了局部投影算法的降噪功效，依此类推，如果令

(6)

适用于所有的单项式

(7)

这个公式是用来证明对于k=1,2,…,20，由k个线性方程组成的系统

n=0,2,…,2(k-1)

(8)

的唯一解是

(9)

显然，由式(9)所确定的系数包括了Cawley等(i=1)和Sauer(i=2)提出的两种估算邻域质心的方法，而且还将邻域质心的估计推广到更高阶(i≥3)。

综上所述，高阶局部投影算法的实现步骤如下：

步骤1确定合适的嵌入维数m和延迟时间τ，将一维含噪信号重构到m维相空间。

步骤2给定邻域数目k，确定每个相点的邻域Un。

步骤3利用前面提到的高阶质心估计方法取代传统的质心估计方法，计算出每个邻域的质心。

步骤4计算每个相点的协方差矩阵C，并求得C的特征值和对应的特征向量。

步骤5确定m0个较大的特征值和剩余的m-m0个较小的特征值对应的特征向量。

步骤6利用式(2)进行修正，从而消除噪声。

步骤7返回步骤2，直至所有相点处理完后将其还原成一维时间序列，这样就得到了降噪后的信号。

3 仿真计算与工程应用

为了验证本文所提方法的有效性和实用性，本文首先采用调频信号和混沌信号进行数值仿真，并与奇异谱降噪、平滑局部子空间投影[14]和标准局部投影三种基于相空间重构的非线性时间序列降噪方法做了对比，最后对从工业现场测得的风机轴承数据进行了分析，并成功地提取出了风机轴承的故障特征。

由于这四种方法都是基于相空间重构，而嵌入维数m和延迟时间τ是其中两个非常关键的参数。目前针对嵌入维数的选取方法有：伪最小邻近法、关联维数法和奇异值分解法等。文献[15]中提出了一种嵌入维数自适应选取方法，并通过仿真和计算验证了该方法的有效性，本文将其作为嵌入维数的选择依据。而针对延迟时间的选取，通常使用平均互信息曲线第一次达到极小值时的τ值作为延迟时间。但是从数学理论的角度上讲，τ可以任意取值，本文经过多次仿真实验证明，τ值的选取对降噪效果影响很微弱。另外，τ取较小值时，邻域点能够较好地逼近吸引子的主流形，所以文中延迟时间τ均取值为1。

3.1 对调频信号的数值仿真

采用的仿真信号为

x=cos(70πt+sin(20πt))

(10)

其中采样点数N=8 000，采样频率fs=8 000 Hz，添加一定水平的高斯白噪声，使其添加噪声后(降噪前)的信噪比为3 dB。然后使用高阶局部投影算法进行降噪处理。降噪效果如图2所示，图2(a)是添加一定噪声后信噪比为3 dB的时域图，图2(b)是使用高阶局部投影算法降噪后的时域图。可以看到，降噪后的时域图上基本看不出噪声的影响，说明高阶局部投影算法具有较好的降噪效果。

(a)含噪信号的时域图(b)降噪后的时域图

图2 正弦信号降噪前后时域图

Fig.2 Time waveforms of sinusoidal signal before and after noise reduction

为了进一步说明高阶局部投影算法的降噪效果，将高阶局部投影算法与奇异谱降噪、平滑局部子空间投影和标准局部投影降噪算法进行对比，使用降噪前后信噪比(Signal-to-Noise Ratio,SNR)的增益作为衡量降噪效果的标准。其中SNR计算式为

(11)

同样采用之前提到的调频信号进行仿真，添加不同水平的高斯白噪声，使其添加噪声后(降噪前)的信噪比分别为0 dB、3 dB、6 dB、9 dB、12 dB、15 dB和18 dB。然后对这7组数据分别使用这四种方法进行降噪处理。采样点数N=8 000，采样频率fs=8 000 Hz。计算结果均是经过多次计算后取平均值以减小偶然误差，图3可以看出高阶局部投影算法相比其它几种非线性时间序列降噪方法具有更好的降噪效果。

图3 四种降噪方法的对比结果

3.2 对混沌信号的数值仿真

本文采用Lorenz混沌时间序列进行仿真分析。其方程为

(12)

式中：σ=10,r=28,b=8/3，此时Lorenz系统呈现出混沌特性。取Lorenz系统的x分量进行仿真分析，采样点数N=3 000，采样频率fs=300 Hz，设定积分步长为0.05，用四阶Runge-Kutta法进行计算。对Lorenz系统的x分量添加一定水平的高斯白噪声，使其添加噪声后(降噪前)的信噪比为3 dB。然后使用高阶局部投影算法进行降噪处理，降噪效果如图4～图6所示。

图4 不含噪Lorenz信号的相图

图5 含噪Lorenz信号的相图

图6 降噪后Lorenz信号的相图

图4是未添加噪声前的相图，图5是添加一定噪声后信噪比为3 dB的相图，图6是使用高阶局部投影降噪后的相图。可以看到，即使在噪声比较大的情况下，高阶局部投影仍然具有较好的降噪效果，并且能够很好地保持吸引子的几何形状。

为了进一步说明高阶局部投影算法的降噪效果，将高阶局部投影算法与奇异谱降噪、平滑局部子空间投影和标准局部投影算法对比。同样采用之前提到的混沌信号进行仿真，添加不同水平的高斯白噪声，使其添加噪声后(降噪前)的信噪比分别为0 dB、3 dB、6 dB、9 dB、12 dB、15 dB和18 dB。然后对这7组数据分别使用这四种方法进行降噪处理。采样点数N=8 000，采样频率fs=8 000 Hz，分别计算四种方法降噪后的信噪比。为了减小偶然误差，计算结果均是经过多次计算后取平均值得到的。图7可以看到高阶局部投影具有最大的信噪比增益，说明了高阶局部投影算法在保留吸引子确定性结构的基础上对非线性信号具有较强的处理能力。

3.3 工程应用

大脱硫风机作为某炼钢厂的关键设备，在环保除尘，治理环境污染方面起着重要的作用。在实际生产过程中一旦发生故障，将导致转炉内的灰尘无法及时清除，不仅会造成环境的污染，而且有可能造成巨大的经济损失和安全隐患。设备的结构简图如图8所示，现场检测人员使用美国CSI2130振动采集仪在轴承3的垂直方向进行信号采集，2006年9月发现该轴承出现异常的振动，将采集的信号分别使用高阶局部投影算法和标准局部投影算法进行分析。轴承型号为22344CA，采样频率为25 600 Hz，转速为740 r/min，这里取12 000个点进行分析，计算得到的故障特征频率如表1所示。

图7 四种降噪方法的对比结果

图8 大脱硫风机的结构简图

转频f/Hz内圈fi/Hz外圈fo/Hz滚动体fr/Hz保持架fc/Hz12.3396.2064.1328.864.93

(a)原始信号时域图(b)高阶局部投影降噪后的时域图

图9 轴承振动信号降噪前后时域图

图10 轴承振动信号降噪前后时域图

Fig.10 The spectrum of bearing vibration signal before and after noise reduction

从图9和图10，即时域图和频域图中可以看出，经过本文所提出的方法降噪处理后，背景噪声大大降低了。将频域图放大后，在没有使用高阶局部投影降噪处理前，由于受到大量噪声的干扰，很多故障特征频率被噪声所淹没，无法判断出轴承的故障类型。而经高阶局部投影降噪处理后，由图11(b)可以清晰地观察到轴的转频12.8 Hz，故障频率64 Hz及其二倍频128 Hz、三倍频192 Hz和四倍频256 Hz。图12是经过标准局部投影降噪后的频域放大图，可以看到故障频率64 Hz及其二倍频128 Hz的幅值明显低于图11(b)中相对应频率的幅值，因此高阶局部投影算法能够更加有效地提取故障特征。虽然与理论计算得到的外圈故障特征频率值有微小的偏差，但在误差允许的范围内，符合要求。因此可以判断是轴承外圈发生了故障，后来经检测证实为外圈故障。

(a)原始信号频谱放大图(b)高阶局部投影降噪后的频谱放大图

图11 轴承振动信号降噪前后频谱放大图

Fig.11 The spectrum zoom of bearing vibration signal before and after noise reduction

图12 标准局部投影降噪后的频谱放大图

通过对以上数值仿真信号和实测的风机轴承振动信号的分析表明：使用高阶多项式对邻域质心进行估计能够更好地逼近吸引子的几何结构，减少了吸引子的扭曲变形，进一步抑制了噪声，提升了局部投影算法的降噪效果。

4 结 论

针对局部投影降噪算法过程中邻域质心的选取对吸引子形状的影响，使用高阶多项式对邻域质心进行了更精确地估计，减少了吸引子的扭曲变形，进一步抑制了噪声，提升了局部投影算法的降噪效果。通过使用调频信号和Lorenz混沌信号进行数值仿真并和其它基于相空间重构的非线性降噪方法对比，凸显了高阶局部投影算法的优越性。最后将其运用于工业现场风机轴承的故障诊断中，有效降低了背景噪声，并从频域成功提取出了轴承的故障特征，证明了高阶局部投影算法的实用性和可行性。

[1] 阳建宏, 徐金梧, 吕勇, 等. 基于相重构和主流形识别的非线性时间序列降噪方法[J]. 北京科技大学学报, 2005,27(5):631-634.

YANG Jianhong, XU Jinwu, LÜ Yong, et al. Nonlinear time series noise reduction method based on phase reconstruction and principal manifold learning[J].Journal of Beijing University of Science and Technology, 2005,27(5):631-634.

[2] CAWLEY R, HSU G H. Local-geometric-projection method for noise reduction in chaotic maps and flows[J]. Physical Review A, 1992,46(6):3057-3082.

[3] SAUER T. A noise reduction method for signals from nonlinear systems[J].Physical D,1992,58(1):193-201.

[4] KANTZ H, SCHREIBER T. Nonlinear time series analysis[M]. Cambridge: Cambridge University Press, 2003.

[5] 吕志民, 张武军, 徐金梧. 基于奇异谱的降噪方法及其在故障诊断技术中的应用[J]. 机械工程学报, 2000,35(3):95-99.

LÜ Zhimin, ZHANG Wujun, XU Jinwu. A noise reduction method based singular spectrum and its application in machine fault diagnosis[J].Journal of Mechanical Engineering, 2000,35(3):95-99.

[6] 徐金梧, 吕勇, 王海峰. 局部投影算法及其在非线性时间序列分析中的应用[J]. 机械工程学报, 2003,39(9):146-150.

XU Jinwu, LÜ Yong, WANG Haifeng. Local projective method and its application on nonlinear time series[J].Journal of Mechanical Engineering,2003,39(9): 146-150.

[7] 阳建宏, 徐金梧, 杨德斌, 等. 邻域自适应选取的局部投影非线性降噪方法[J]. 振动与冲击, 2006,25(4):64-67.

YANG Jianhong, XU Jinwu, YANG Debin, et al. Nonlinear noise reduction method by local adaptive neighborhood selection[J].Journal of Vibration and Shock, 2006,25(4):64-67.

[8] 韩敏, 刘云侠. 基于非线性约束的局部投影降噪[J]. 电子与信息学报, 2009,31(2):400-404.

HAN Min, LIU Yunxia. Local projection noise reduction based on nonlinear constraints[J].Journal of Electronics and Information Technology,2009,31(2):400-404.

[9] MICHAEL T J, RICHARD J P. Generalized phase space projection for nonlinear noise reduction[J]. Physica D, 2005,201(3/4):306-317.

[11] KOTAS M. Projective filtering of time warped ECG beats[J]. Computers in Biology and Medicine, 2008,38(1):121-137.

[12] MOORE J M, SMALL M, KARRECH A. Improvements to local projective noise reduction through higher order and multiscale refinements[J]. Chaos, 2015,25(6): 653-671.

[13] TAKENS F. Detecting strange attractors in turbulence[J]. Lecture Notes in Mathematics,1981, 898(1):366-381.

[14] CHELIDZE D. Smooth local subspace projection for nonlinear noise reduction[J]. Chaos, 2014,24(1): 274-276.

[15] 余成义. 流形学习在混沌时间序列降噪中的应用[D]. 武汉: 武汉科技大学, 2012.