唐晓红, 胡俊锋,2, 熊国良, 张 龙

(1. 华东交通大学机电与车辆工程学院 南昌，330013) (2. 中国铁路南昌局集团有限公司科学技术研究所 南昌，330002)

引 言

在旋转机械中，滚动轴承是运用最为广泛的关键零部件之一。随着机械内部结构复杂程度的不断提高，各零部件之间的联系显得极为紧密，而轴承的运转状态决定了其常处于高速高温的恶劣环境下，极易发生损伤类故障[1-2]。为了提高机械设备的安全性和可靠性，减少不必要的经济损失及人员伤亡，开展机械设备状态监测与故障诊断研究成为现代工业化建设亟待解决的热点问题[3]。近年来，对于机械设备的故障诊断研究已经取得了一定进展，一系列优异的消噪算法如谱峭度[4]、小波变换[5]和最小熵解卷积[6]等均得到了国内外学者的广泛关注并成功应用于轴承的故障诊断研究，然而这些算法均存在一定的局限性，如谱峭度只反映冲击性不反映周期性，小波变换基函数选取困难，这些不足均会导致诊断效果不佳。

非局部均值算法是近年来图像处理领域的新兴算法之一，是Buades等[7]为了解决传统消噪方法在图像信号处理出现的无法很好的保留图像细节及边缘的问题而提出的一种新的消噪方法。该算法在图像处理领域的应用极为广泛，在ECG信号[8]及叠前地震信号[9]的处理中也体现出优异的消噪效果。文献[10]将NLM算法应用于轴承振动数据中的测量噪声消除，然后将处理后的信号通过经验模态分解(empirical mode decomposition，简称EMD)来消除其他噪声干扰，最后对本征模函数(intrinsic mode function，简称IMF)包络谱进行分析,成功提取出故障特征频率。Lü等[11]对传统的NLM算法进行了适当改进，提出了快速NLM算法。文献[12-13]为了解决NLM算法处理低信噪比时出现的故障冲击均值化问题，利用加权运算得出的权重曲线以突显故障冲击特征。总体来说，NLM算法在轴承故障诊断研究中仍然处于初级阶段，存在许多问题有待解决，如NLM算法的滤波效果很大程度上取决于λ，M和P 3个参数的选择，而上述研究均未对该问题进行深入探讨。

笔者针对NLM参数难以选取的问题，提出采用粒子群算法(particle swarm optimization, 简称PSO)对NLM的3个决定性参数进行优化选择，实现NLM参数的自适应选取，并通过轴承故障仿真数据和实验数据对该方法进行了验证。

1 方法介绍

1.1 非局部均值算法

NLM算法的原理是利用图像中存在的众多相似结构的特性，以像素点所在领域为单元，在图像中搜索与该单元相似的成分，并对这些相似成分加权平均操作，以此消除图像中的异常噪声。Buades等[7]提出NLM算法的同时证明了其一致性，即NLM算法并非只对特定的像素点有效，而是适用于图像中的任一像素点，即在自然图像中，任一结构均存在与之类似的结构样本。这些相似特性不仅表现在二维图像中，在一维信号中也具有广泛体现。

假设u(t)为真实轴承故障信号，n(t)为背景噪声干扰，则实际采样得到的带噪声信号v(t)表示为

v(t)=u(t)+n(t)

(1)

对于v(t)而言，基于NLM的数据处理的目的是从包含噪信号v(t)中消除噪声干扰n(t)，从而恢复真实状态下的故障信号u(t)的过程。u(t)的NLM算法估计值K(t)为

(2)

其中：N为以点t为中心的搜索范围；N(t)为搜索区域内所有点的集合。

Z(t)为归一化因子，其计算公式为

Z(t)=∑s∈N(t)ω(t,s)

(3)

其中：ω(t,s)表示权重。

ω(t,s)计算公式为

(4)

其中：d2为t和s两点邻域块之间欧氏距离的平方和，欧式距离越小，两者的相似程度越大；参数h=λ(2LΔ)1/2，控制权重ω(t,s)的衰减速度，即直接决定信号的最终滤波程度；Δ表示以点t为中心的邻域块；LΔ为以点s为中心的邻域块；λ为带宽参数。

ω(t,s)的大小取决于以t和s为中心的两个邻域块之间的相似程度，两者相似度越高，则ω(t,s)的取值越大，反之亦然。ω(t,s)取值满足0≤ω(t,s)/Z(t)≤1和∑jω(t,s)=1两个基本条件。

参数λ，M和P为NLM算法的决定性参数，控制着NLM的最终处理效果。三者关系如图1所示，参数P控制着邻域块的大小，即LΔ=2P+1；参数M为搜索区域N(t)长度的一半，理论上来说，M越大，则N(t)中的相似邻域块越多，NLM加权平均后的处理效果越好，但是相应地需要更多的计算时间，严重影响计算效率。

图1 NLM参数关系Fig.1 Illustration of NLM parameter relationship

1.2 粒子群算法

一群鸟在区域内寻觅食物，若该区域仅有一块食物，那么找到食物的最简单有效的方案便是搜索距离食物最近的鸟的周围区域。受此启发，Kennedy等[14]提出粒子群算法，该算法中的粒子飞行的行为规则类似于鸟类运动，从而使整个粒子群的运动表现出与鸟类觅食类似的特性，并通过群体中个体之间的协作和信息共享来寻求最优解。该算法鲁棒性强、收敛速度快且并不依赖特定的求解模型，非常适合求解模型的优化问题[15]。

PSO在求解最优问题时，其最优解对应搜索区域中的一个粒子的位置。该区域中存在许多粒子，每个粒子都具有自身独特的位置和速度以决定最终的飞行距离和方向。在搜索最优粒子前，首先需要初始化一些随机粒子，每个随机粒子均表示区域内的候选解，这些候选解的优劣程度将由适应度值进行评估。粒子通过不断迭代来寻找最优解，在每一次的迭代中，粒子通过跟踪个体极值点Pbest(粒子本身找到的最优解)和全局极值点Gbest(整个种群目前找到的最优解)来更新自己。

设d维空间中粒子i的当前位置和速度分别为Xi=[xi,1, xi,2, …, xi,j]，Vi=[vi,1, vi,2, …, vi,j]；t时刻粒子i所经历过的个体最优位置Pi,j=[pi,1, pi,2, …, pi,j]，在此记为Pbest；则该时刻群体中所有粒子经历过的全局最优位置Gj=[g1, g2, …, gj]，记为Gbest。对于每一代粒子，其第j维粒子在寻觅到Pbest和Gbest两个极值后，将通过这两个最佳位置按照下式更新各粒子的速度和位置

vi,jt+1=ωvi,jt+c1r1pi,j-xi,jt+

c2r2pg,j-xi,jt

(5)

xi,jt+1=xi,jt+vi,jt+1

(j=1,2,…,d)

(6)

其中：ω为惯性权重；c1和c2为正的学习因子，亦被称为加速常数，一般来说c1=c2，取值在0～4之间；r1和r2为0～1之间均匀分布的随机数。

在PSO迭代寻优过程中，需要对位置变量Xi,j和速度变量Vi,j进行约束以确保算法的收敛性，约束条件如下：当Xi,jXmax时，Xi,j=Xmin+rand(Xmax-Xmin)；当Vi,jVmax时,Vi,j=Vmax。其中：Xmin和Xmax分别为粒子处于每一维时候的最小和最大位置；rand(·)为0～1之间随机变化的数值；Vmax为粒子允许的最大运行速度。

2 PSO优化NLM的滚动轴承故障诊断模型

在NLM算法中，λ，M和P控制着该算法的最终处理效果，是NLM算法的决定性参数。利用PSO算法的寻优特性对这3个参数进行最优求解。在PSO算法中，目标函数即适应度的选择决定了最终求解的准确性，对于滚动轴承故障数据而言，故障诊断的实质在于消除背景噪声干扰、从而突出能够反映轴承故障的相关信息，诊断精度一定程度上取决于故障信息的提取程度。峭度指标是常用的判别故障严重与否的评估指标之一，对异常冲击特别敏感，能够随着异常冲击的凸显而增大。笔者采用峭度的倒数作为PSO的目标函数，即滤波效果越好，冲击凸显越明显，则峭度值越大，其倒数(目标函数适应度)越小。PSO优化NLM的滚动轴承故障诊断模型如图2所示。首先，设置好PSO参数；其次，以原始振动数据为对象，通过PSO的参数迭代逐步逼近期望输出，直至实际输出满足阈值条件或达到最大迭代次数，获得最优参数组合；最后，以最优参数构成NLM滤波器对信号滤波处理，通过包络谱分析得出诊断结果。诊断步骤如下:

1) 输入原始信号(signal)；

2) 给定初始条件，设置学习因子c1=c2=1.496 18，权重ω=0.729 8，适应度值阈值为0.01，最大迭代次数K=50，种群个体数目N=50；

3) 构造N个长度为3(参数λ，M和P)的向量X和V，分别代表种群的位置和速度，并初始化；

4) 计算粒子的适应度值，并初始化Pbest和Gbest；

5) 更新粒子位置和速度；

6) 比较每个粒子的适应度与经历最优位置时的适应度Pbest，若当前更好，则更新Pbest；

7) 比较每个粒子的适应度和群体经历最优位置的适应度Gbest，若当前更好，则更新Gbest；

8)重复步骤5～7，N次后获得最优Pbest和Gbest；

9) 判断Pbest(实际输出)是否小于阈值0.01(期望输出)，若不满足条件则重复步骤5～8，直至满足阈值条件或达到最大迭代次数进入下一步；

10) 输出最优λ，M和P参数；

11) 将最优参数代入NLM算法获得最优滤波器；

12) 对原始振动信号滤波处理得到NLM滤波信号；

13) 对滤波信号包络计算后进行包络谱分析，从中提取出故障特征频率；

14) 提取出的故障特征频率与理论故障特征频率对比，判别是否存在故障，得到诊断结果。

图2 PSO优化NLM的滚动故障诊断模型Fig.2 Rolling bearing fault diagnosis model with NLM optimized by PSO

3 仿真信号分析

为了验证所提方法在轴承故障诊断中的可行性和有效性，建立轴承响应函数并通过“rand”函数加入高斯噪声干扰以模拟现实状态下的轴承外圈故障数据，设置仿真信号的采样频率fs=12 kHz，外圈故障特征频率为80 Hz。外圈故障仿真信号时域波形图如图3(a)所示。由于加入高斯噪声程度较大，原始信号中的极大一部分故障冲击特征被背景噪声掩盖，无法直观地判别故障冲击的位置，亦难以判别故障类型。因此需要对图3(a)进一步处理消除背景噪声以凸显故障冲击成分。图3(b)PSO收敛曲线可以看到随着参数的更替及迭代次数的递进，适应度值也在不断减小，即PSO寻优求解过程中目标函数适应度(滤波信号峭度的倒数)不断减小，相应地滤波信号后其噪声逐渐消除，冲击凸显程度不断增强，峭度随之增大。当迭代次数到达50次时寻优过程结束，得到最优参数组合λ= 0.340 6，M=5，P=5。将最优参数代入NLM得到最优NLM滤波器，并对原始信号进行处理，得到的滤波信号如图3(c)所示。相较图3(a)，图3(c)经过NLM滤波后的信号噪声干扰被压制到极小范围，冲击特征提取效果明显。图3(d)为对图3(c)包络计算后的包络谱图，发现能够成功提取出故障特征频率和倍频成分，而其他频率干扰较小，可见笔者提出的方法具有良好的消噪效果。

图3 仿真信号PSO-NLM诊断过程Fig.3 Diagnostic procedure on simulation Signal using PSO-NLM

4 实验数据分析

4.1 外圈故障分析

实验室外圈故障数据来源于美国辛辛那提大学智能维护中心(IMS)[16]，轴承疲劳试验台的整体结构和局部照片分别如图4(a)和图4(b)所示。电机提供动力，通过带传动驱动主轴旋转，主轴上依次安装有4个型号为Rexnord ZA-2115的双列圆锥滚子轴承。实验过程中通过安装于轴承座上的加速度传感器对轴承状态进行检测，实验时主轴转速设定为2 kr/min，采样频率为20 kHz，采样间隔为10 min。整个实验过程共持续164 h，共采集数据文件984个，每个数据文件包含4列长为20 480点的数据，分别对应主轴上4个轴承的状态。实验结束后拆解装置发现轴承1出现了严重外圈故障。Rexnord ZA-2115轴承参数及计算后得到的理论故障特征频率如表1所示。

图4 轴承疲劳试验台Fig.4 Fatigue test rig of rolling bearing

表1RexnordZA-2115轴承参数及故障特征频率

Tab.1RexnordZA-2115bearingparameterandfaultcharacteristicfrequency

轴承参数轴承节径/mm滚动体直径/mm滚珠个数接触角/(°)71.58.41615.17故障频率/Hz外圈内圈滚动体保持架236.43296.90140.0414.78

本节分析的故障数据为984个文件中的第976个第1列数据，此时接近实验尾声，外圈故障已经较为严重，代表性较强。因数据点过大，在此仅截取其中的3 000点数据进行分析。图5(a)为其时域波形图，可发现由于背景噪声较大，时域波形稍显杂乱无章，因轴承故障引发的脉冲冲击特征体现并不清晰，极大一部分故障脉冲隐藏于背景噪声中。NLM算法自适应过程中的目标函数收敛曲线如图5(b)所示，发现在优化过程中目标函数适应度逐步缩小，而NLM的峭度随之逐渐增大，直至迭代至50次时优化结束，其最优参数最终确定为λ=0.144 1，M=7，P=6。根据该参数代入NLM算法得到的最优滤波器并对原始信号图5(a)处理后的时域波形图如图5(c)所示。可以看到，处于图5(a)的噪声部分得到了极大程度的抑制，而能够代表轴承故障信息的脉冲冲击取得了极大的增强效果，相对于图5(a)来说，故障点和故障位置已经能够基本确定。图5(d)滤波包络谱能够清晰地提取出故障特征频率及其倍频成分，而频率的噪声干扰成分被压制到极小，因此可以判别出该轴承存在外圈故障，与实际情况相符。

图5 实验室外圈信号PSO-NLM诊断过程Fig.5 Diagnostic procedure on laboratory outer ring signal using PSO-NLM

4.2 内圈故障分析

为了进一步说明本研究方法在轴承故障诊断的有效性，采用美国凯斯西储大学轴承数据中心[17]的轴承内圈故障数据进行分析。滚动轴承故障模拟实验台主要由左侧的三相感应电机、中间的力矩传感器、右端的测功机及内部控制模块4部分组成。实验台的整体结构如图6所示。

图 6 滚动轴承故障模拟实验台Fig.6 Simulation test-bed of rolling bearing failure

该实验台可测试两种型号的轴承，型号分别为SKF6205-2RS和SKF6203-2RS，分别安装在电机驱动端和风扇端。采用电火花技术为轴承添加不同尺寸的点蚀故障，故障直径依次为0.178，0.356和0.533 mm。轴承的工作状态由安装于电机驱动端、风扇端和底座的3个加速度传感器进行记录，由此采样可得到所需的故障振动数据。选取的分析数据为驱动端SKF6205-2RS轴承内圈故障数据，其故障直径为0.533 mm，实验时的相关参数如下：电机转速为1 797 r/min，采样频率为12 kHz。SKF6205-2RS轴承的相关参数及实验状态下的故障特征频率如表2所示。

表2SKF6205-2RS轴承参数及故障特征频率

Tab.2SKF6205-2RSbearingparameterandfaultcharacteristicfrequency

轴承参数轴承节径/mm滚动体直径/mm滚珠个数接触角/(°)39.047.9490故障频率/Hz外圈内圈滚动体保持架107.36162.1970.5811.93

图7(a)为安装在电机驱动端的加速度传感器采样得到的振动数据。为了凸显本研究方法的诊断效果，通过“rand”函数加入了一定程度的高斯白噪声干扰。由于噪声比例较大，原始信号中代表故障信息的循环冲击成分被背景噪声压制，很难从中提取出故障信息。通过PSO优化后获得最优参数λ，M，P的数值分别为0.509 9，62和6，将之代入NLM算法得到最优滤波器。图7(c)和图7(d)分别为NLM滤波信号及其包络谱。图7(c)滤波信号中能够明显的提取出故障冲击特征，而图7(d)中故障特征频率基频及倍频成分均清晰明了且谱线突出，倍频谱线两端存在明显的转频边带，为典型的内圈故障，与内圈故障特征频率理论值也完全对应。可见本研究方法在处理轴承内圈故障也同样适用。

图7 实验室内圈信号PSO-NLM诊断过程Fig.7 Diagnostic procedure on laboratory inner ring signal using PSO-NLM

5 方案对比

在文献[10-11]中，使用非局部均值算法对数据进行处理时，参数λ的取值通过SURE准则选择。由于λ的取值直接影响到消噪后信号的平滑度，因此在上述文献中根据SURE准则确定参数λ的取值为0.5σ，其中σ为信号的方差。由于参数M和P通过工作人员的经验进行选取，因此数据处理结果受工作人员的主观性思想影响较大，同样的数据经过不同的工作人员呈现出结果的多样性。为了便于直观地体现本研究方法的优异性，以图5(a)的外圈故障数据为例进行举例说明。通过SURE准则选取λ=0.5，σ=0.187 8，根据经验选择参数M和P的值为10和10。将参数代入NLM后得到滤波器，对图5(a)信号处理后得到的滤波信号如图8(a)所示，其在时域波形图中并未有明显改善，故障冲击仍处于噪声淹没状态，而图5(c)的噪声消除效果极为优异，故障冲击凸显效果优势较大。对比双方包络谱，图8(b)虽然能得到故障特征频率及其倍频成分，但是4倍频和5倍频成分却被噪声频率压制，而图5(d)则能清晰提取出4倍频和5倍频成分，且图5(d)频率的峰值能量也比图8(b)略高。可见，从噪声消除效果和包络谱分析效果看，本研究方法优于传统的NLM方法。

6 结束语

提出应用粒子群算法对参数进行寻优求解，实现了NLM算法的自适应化。针对轴承故障诊断领域建立了PSO优化NLM的滚动轴承故障诊断模型。运用仿真数据和实验轴承内、外圈故障数据对所提模型进行了验证。结果证明，本研究方法能够较大幅度地消除噪声，并提取出故障信息，对于算法的普及和应用有着十分积极的作用。