CN105510032A - 基于谐躁比指导的解卷积方法 - Google Patents

Abstract

基于谐躁比指导的解卷积方法，首先对采集的振动信号进行截断和去均值处理；然后对截断后的信号直接进行解卷积处理，兼顾故障的冲击性和周期性特征，将峭度和谐躁比作为目标函数进行联立求迭代表达式，进而极大地减小了随机冲击的干扰；最后对解卷积后的信号进行包络分析，从包络谱中能提取故障特征频率，本发明在特征频率的提取过程中不需要人为参与，有利于实现故障特征提取和诊断监测的自动化，节约时间，效率更高。

Description

基于谐躁比指导的解卷积方法

技术领域

本发明涉及机械设备故障诊断技术领域，特别涉及基于谐躁比指导的解卷积(Harmonics-to-NoiseRatioGuidedDeconvolution，HNRGD)方法。

背景技术

振动分析为现阶段机械设备故障诊断最为有效途径之一，机械设备的状态劣化往往表现为振动信息的变化或异常。峭度是描述信号冲击性的有效指标，当振动信号表现为明显的冲击性时其峭度值较大，被广泛应用于机械系统监测诊断中。2007年最小熵解卷积(Minimumentropydeconvolution，MED)方法被Sawalhi和Randall等学者首次应用于滚动轴承的故障诊断领域，并取得了一定的效果。MED是一种不需要任何先验假设的信号时域盲解卷积技术，通过迭代选择一个(finiteimpulseresponse)FIR滤波器来最小化滤波信号的输出熵(即峭度最大化)，旨在最小化噪声的同时提取出故障冲击，因此其在高噪声下也能得到理想的诊断结果。虽然MED对冲击性的增强与提取具有明显的效果，但是其目标函数仅仅是寻求滤波信号峭度值最大化，所以易受随机孤立冲击的干扰。即当故障信号中既存在周期性冲击序列和随机孤立冲击时，通过MED增强技术有可能只能增强孤立冲击，而对真正反映故障特征的周期性冲击序列没有效果。鉴于上述问题，2012年McDonald等学者提出了最大相关峭度解卷积(Maximumcorrelatedkurtosisdeconvolution,MCKD)算法，提出了相关峭度的概念，兼顾了故障的冲击性和周期性特征，减小了随机冲击的干扰，并成功地将其应用到齿根裂纹故障诊断上。然而MCKD方法需要提前准确估计故障特征频率，而在工程实际中，测速计往往存在误差，而且复杂的工况使得设备转速不可能保持恒定。并且当转速精准确定的前提下，MCKD也很难提取带有随机滑动的轴承故障冲击信号。这些不足给MCKD的应用推广带来了诸多不便。

发明内容

为了克服上述现有技术的缺点，本发明的目的在于提供基于谐躁比指导的解卷积方法，实现机械设备不用预估故障特征频率等预处理，就能实现精确的故障诊断。

为实现上述目的，本发明方案采取的技术方案为：

基于谐躁比指导的解卷积方法，包括以下步骤：

步骤一：将振动加速度传感器吸附于对被测试滚动轴承的轴承座上，并对振动信号进行高频采样、截断和去均值处理，将振动信号记为x；接下来对振动信号x进行HNRGD处理；

步骤二：将峭度作为目标函数对滤波系数进行求偏导：

$O_k\[f\left(l\right)\]=\frac{\underset{n=1}{\overset{N}{\Σ}}{y}^4\left(n\right)}{{\[\underset{n=1}{\overset{N}{\Σ}}{y}^2\left(n\right)\]}^2}---\left(1\right)$

其中N为输出信号y的长度，使用其中f(l)为滤波器系数，l＝1,2,……,L,L为滤波器长度，求导后的结果为：

$\frac{\underset{n=1}{\overset{N}{\Σ}}{y}^2\left(n\right)\underset{n=1}{\overset{N}{\Σ}}{y}^3\left(n\right)x(n-l)}{\underset{n=1}{\overset{N}{\Σ}}{y}^4\left(n\right)}=\underset{n=1}{\overset{N}{\Σ}}x(n-l)y\left(n\right)---\left(2\right)$

再以谐躁比(Harmonics-to-NoiseRatio，HNR)为目标函数对滤波系数进行求偏导：

$O_k\[f\left(l\right)\]=\frac{\∫y\left(t\right)y(t+T)dt}{\∫y^2\left(t\right)dt-\∫y\left(t\right)y(t+T)dt}---\left(3\right)$

其中t为时间，T为周期，对目标函数进行离散处理，并使用 $\frac{\∂y\left(n\right)}{\∂f\left(l\right)}=x(n-l),$ 求导后的结果为:

$\begin{array}{c}\[\underset{n=1}{\overset{N}{\Σ}}x(n-l)y(n+T)+\underset{n=1}{\overset{N}{\Σ}}x(n+T-l)y\left(n\right)\]\\ =2\underset{n=1}{\overset{N}{\Σ}}x(n-l)y\left(n\right)\·\underset{n=1}{\overset{N}{\Σ}}y(n+T)y\left(n\right)/\underset{n=1}{\overset{N}{\Σ}}{y}^2\left(n\right)\end{array}---\left(4\right)$

联立(2)和(4)，使用得到解卷积迭代表达式：

将上式写成矩阵的形式：

Af＝b(6)

其中：

b--逆滤波器的输入振动信号x、输出y的互相关，b为L维列向量；

A--输入振动信号x的自相关，A为L×L维矩阵；

f--逆滤波器的滤波器系数，f为L维列向量；

首先计算自相关矩阵A；再假设逆滤波器的初始值f⁽⁰⁾，设置滤波器长度L＝100，并给定初始滤波器系数为[00……1-1……00]，用y⁽⁰⁾和x⁽⁰⁾计算列向量b⁽¹⁾；然后求解新的滤波器系数f⁽¹⁾＝A^-1b⁽¹⁾；通过对滤波器系数的计算进行收敛判断；如果通过比较后一次迭代与前一次迭代的谐躁比与峭度值，如果这两个参数都在增加那么迭代结束，如果并不是都在增加就更新滤波器系数进行下次迭代，直到满足收敛条件或者达到设定的最大迭代次数为30次，得到经过HNRGD处理后的信号y_k；

步骤三：对HNRGD处理后的信号y_k进行包络分析并得到包络谱，对包络谱进行分析，进而提取故障特征频率。

本发明相比于现有技术，具有以下有益效果：

a)本发明在传统MED方法的基础上，融入了谐躁比的概念，能兼顾故障的冲击性和周期性特征，极大地减小了随机冲击的干扰。

b)本发明不需要任何先验知识，也无需对系统故障特征频率进行预估，方法具有鲁棒性。

c)本发明在特征频率的提取过程中不需要人为参与，有利于实现故障特征提取和诊断监测的自动化，节约时间，效率更高。

附图说明

图1为实施例美国西储大学轴承实验台示意图。

图2为本发明方法流程图。

图3为实例中对原信号进行截断处理后的振动信号x。

图4为实例中对截断后的信号进行HNRGD方法处理后的信号y_k。

图5为实施例中经HNRGD方法处理后信号的包络谱。

图6为实施例中利用MCKD方法对截断信号进行改变采样率重采样之后的信号。

图7为实施例中从采样信号经MCKD滤波后的信号。

图8为实施例中经MCKD滤波后的信号的包络谱。

具体实施方式

下面结合附图和实施例对本发明做详细描述。

以下以美国西储大学轴承实验台为例进行说明，该试验台由电机、联轴器、转矩传感器、功率计以及轴承等组成，如图1所示，该实验台检测对象为驱动电机的驱动端轴承和风扇端轴承，轴承型号分别为：6205-2RSJEMSKF和6203-2RSJEMSKF，采样频率为12kHz。整个装置由电机驱动，将扭矩从电机沿转矩传感器传递到输出端的功率计。

本发明采用的数据是由美国西储大学提供的上述试验台测得的数据中编号为OR0146_3.mat的一组数据。根据W.A.SmithandR.B.Randall近期发表文章《RollingelementbearingdiagnosticsusingtheCaseWesternReserveUniversitydata:Abenchmarkstudy》描述，该组数据通过传统的方法无法进行故障信息的提取，因此该数据具有说服力。

利用基于谐躁比指导的解卷积方法对该试验台中驱动电机驱动端轴承进行诊断，应用本发明对实验数据进行分析并和传统MCKD方法进行对比。

如图2所示，基于谐躁比指导的解卷积方法，包括以下步骤：

步骤一：将振动加速度传感器吸附于对被测试滚动轴承的轴承座上，并对振动信号进行高频采样、截断和去均值处理，其中采样频率为12kHz，在使用数据时需要去掉起始噪声部分，截取整段信号中0.5-2.0s共1.5s的数据，如图所示3，将该1.5s振动信号记为x；接下来对振动信号x进行HNRGD处理；

步骤二：将峭度作为目标函数对滤波系数进行求偏导：

$O_k\[f\left(l\right)\]=\frac{\underset{n=1}{\overset{N}{\Σ}}{y}^4\left(n\right)}{{\[\underset{n=1}{\overset{N}{\Σ}}{y}^2\left(n\right)\]}^2}---\left(1\right)$

其中N为输出信号y的长度，使用其中f(l)为滤波器系数，l＝1,2,……,L,L为滤波器长度，求导后的结果为：

$\frac{\underset{n=1}{\overset{N}{\Σ}}{y}^2\left(n\right)\underset{n=1}{\overset{N}{\Σ}}{y}^3\left(n\right)x(n-l)}{\underset{n=1}{\overset{N}{\Σ}}{y}^4\left(n\right)}=\underset{n=1}{\overset{N}{\Σ}}x(n-l)y\left(n\right)---\left(2\right)$

再以谐躁比(Harmonics-to-NoiseRatio，HNR)为目标函数对滤波系数进行求偏导：

$O_k\[f\left(l\right)\]=\frac{\∫y\left(t\right)y(t+T)dt}{\∫y^2\left(t\right)dt-\∫y\left(t\right)y(t+T)dt}---\left(3\right)$

其中t为时间，T为周期，对目标函数进行离散处理，使用 $\frac{\∂y\left(n\right)}{\∂f\left(l\right)}=x(n-l),$ 求导后的结果为:

$\begin{array}{c}\[\underset{n=1}{\overset{N}{\Σ}}x(n-l)y(n+T)+\underset{n=1}{\overset{N}{\Σ}}x(n+T-l)y\left(n\right)\]\\ =2\underset{n=1}{\overset{N}{\Σ}}x(n-l)y\left(n\right)\·\underset{n=1}{\overset{N}{\Σ}}y(n+T)y\left(n\right)/\underset{n=1}{\overset{N}{\Σ}}{y}^2\left(n\right)\end{array}---\left(4\right)$

联立(2)和(4)使用得到解卷积迭代表达式：

将上式写成矩阵的形式：

Af＝b(6)

其中：

b--逆滤波器的输入振动信号x、输出y的互相关，b为L维列向量；

A--输入振动信号x的自相关，A为L×L维矩阵；

f--逆滤波器的滤波器系数，f为L维列向量；

首先计算自相关矩阵A；再假设逆滤波器的初始值f⁽⁰⁾，为与MCKD具有可比性，将滤波器长度设置与MCKD方法一样为100，并给定初始滤波器系数为[00……1-1……00]，用y⁽⁰⁾和x⁽⁰⁾计算列向量b⁽¹⁾；然后求解新的滤波器系数f⁽¹⁾＝A^-1b⁽¹⁾；通过对滤波器系数的计算进行收敛判断；如果通过比较后一次迭代与前一次迭代的谐躁比与峭度值，如果这两个参数都在增加那么迭代结束，如果并不是都在增加就更新滤波器系数进行下次迭代，直到满足收敛条件或者达到设定的最大迭代次数为30次；得到经过HNRGD处理后的信号y_k，如图4所示；

步骤三：对HNRGD处理后的信号y_k进行包络分析并得到包络谱，如图5所示，BPFO(ballpassfrequencyofouterrace)表示外圈故障特征频率，对包络谱进行分析，可以清楚的看到轴承外圈特征故障频率102.9435Hz以及其2倍频、4倍频、6倍频和7倍频，效果十分明显。

如图6所示为利用MCKD方法对振动信号x进行改变采样率重采样之后的效果，按照McDonald等提出的MCKD方法对重采样后的信号进行处理，对该方法的参数按照McDonald提出的推荐设置，其中参数设置为精确计算故障特征频率f＝102.9435Hz，移位距离T＝fs*1/f＝12000/102.9435，移位次数M＝3，滤波器长度设置为100，最大滤波次数设置为30，图7是经MCKD滤波后的信号，图8为其包络谱，从图8可以看出，在精确知道电机转速，轴承参数等一系列参数之后，该例中还必须知道轴承系统外圈故障以及精确计算出其故障特征频率的众多前提下，才能艰难地从包络谱中看到故障特征频率的6倍频(617.6610Hz)附近的616Hz的频率成分。然而，在工程实际中，系统经常处于变转速变载荷的复杂工况下，而且测速计往往存在误差，很难精确预估故障特征频率。因此相比于此，本发明提出的方法能避免以上提到的困难并能精确的提取故障特征频率，所以本发明优势明显。

Claims (1)

1.基于谐躁比指导的解卷积方法，其特征在于，包括以下步骤：

步骤一：将振动加速度传感器吸附于对被测试滚动轴承的轴承座上，并对振动信号进行高频采样、截断和去均值处理，将振动信号记为x；接下来对振动信号x进行HNRGD处理；

步骤二：将峭度作为目标函数对滤波系数进行求偏导：

$O_k\[f\left(l\right)\]=\frac{\underset{n=1}{\overset{N}{\Σ}}{y}^4\left(n\right)}{{\[\underset{n=1}{\overset{N}{\Σ}}{y}^2\left(n\right)\]}^2}---\left(1\right)$

其中N为输出信号y的长度，使用其中f(l)为滤波器系数，l＝1,2,……,L,L为滤波器长度，求导后的结果为：

$\frac{\underset{n=1}{\overset{N}{\Σ}}{y}^2\left(n\right)\underset{n=1}{\overset{N}{\Σ}}{y}^3\left(n\right)x(n-l)}{\underset{n=1}{\overset{N}{\Σ}}{y}^4\left(n\right)}=\underset{n=1}{\overset{N}{\Σ}}x(n-l)y\left(n\right)---\left(2\right)$

再以谐躁比(Harmonics-to-NoiseRatio，HNR)为目标函数对滤波系数进行求偏导：

$O_k\[f\left(l\right)\]=\frac{\∫y\left(t\right)y(t+T)dt}{\∫y^2\left(t\right)dt-\∫y\left(t\right)y(t+T)dt}---\left(3\right)$

其中t为时间，T为周期，对目标函数进行离散处理，使用 $\frac{\∂y\left(n\right)}{\∂f\left(l\right)}=x(n-l),$ 求导后的结果为:

$\begin{array}{c}\[\underset{n=1}{\overset{N}{\Σ}}x(n-l)y(n+T)+\underset{n=1}{\overset{N}{\Σ}}x(n+T-l)y\left(n\right)\]\\ =2\underset{n=1}{\overset{N}{\Σ}}x(n-l)y\left(n\right)\·\underset{n=1}{\overset{N}{\Σ}}y(n+T)y\left(n\right)/\underset{n=1}{\overset{N}{\Σ}}{y}^2\left(n\right)\end{array}---\left(4\right)$

联立(2)和(4)并使用得到解卷积迭代表达式：

将上式写成矩阵的形式：

Af＝b(6)

其中：

b--逆滤波器的输入振动信号x、输出y的互相关，b为L维列向量；

A--输入振动信号x的自相关，A为L×L维矩阵；

f--逆滤波器的滤波器系数，f为L维列向量；

首先计算自相关矩阵A；再假设逆滤波器的初始值f⁽⁰⁾，设置滤波器长度L＝100，并给定初始滤波器系数为[00……1-1……00]，用y⁽⁰⁾和x⁽⁰⁾计算列向量b⁽¹⁾；然后求解新的滤波器系数f⁽¹⁾＝A^-1b⁽¹⁾；通过对滤波器系数的计算进行收敛判断；如果通过比较后一次迭代与前一次迭代的谐躁比与峭度值，如果这两个参数都在增加那么迭代结束，如果并不是都在增加就更新滤波器系数进行下次迭代，直到满足收敛条件或者达到设定的最大迭代次数为30次，得到经过HNRGD处理后的信号y_k；

步骤三：对HNRGD处理后的信号y_k进行包络分析并得到包络谱，对包络谱进行分析，进而提取故障特征频率。