CN105865784A - Rolling bearing detection method based on LMD (Local Mean Decomposition) and gray correlation - Google Patents

Abstract

The invention discloses a rolling bearing detection method based on LMD (Local Mean Decomposition) and gray correlation, belongs to a rolling bearing fault detection method in the field of mechanical engineering, and relates to a rolling bearing detection method based on a fuzzy entropy algorithm of LMD (Local Mean Decomposition) and gray correlation. The method comprises the steps: employing an acceleration sensor to collect vibration acceleration signals of a rolling bearing during operation, wherein the vibration acceleration signals comprise a no-fault normal bearing vibration acceleration signal and inner ring, rolling body or outer ring fault bearing vibration acceleration signals. carrying out the LMD decomposition of the collected acceleration signals, and obtaining a plurality of PF (product function) components and residual errors; calculating the gray correlation degree of a test sample and a standard matrix through employing a gray correlation algorithm, and carrying out the fault mode recognition. The method can effectively carry out the feature extraction of the vibration signals, solves problems that the EMD decomposition is severe in excessively modal mixing and end effect and a PF component is large in data size and cannot serve as a characteristic vector, and achieves the effective recognition of the operation state of the rolling bearing.

Description

Decompose based on local mean value and the rolling bearing detection method of grey correlation

Technical field

The invention belongs to the fault detection method of the rolling bearing of mechanical engineering field, be specifically related to one The fuzzy entropy algorithm decomposed based on local mean value and the rolling bearing fault testing method of grey correlation.

Background technology

In plant equipment, rolling bearing is very important parts, daily life, commercial production with And the application of the every field such as national defense construction is extremely wide.The quality of the running status of rolling bearing directly shadow Ring stability, reliability and life-span that whole equipment runs.Therefore, the state prison of rolling bearing Survey and fault diagnosis technology, for safety in production, reduce economic loss, it is ensured that machine security is stable Operation has a very important role.

Owing to rolling bearing is affected by working environment, the fault-signal of rolling bearing be mostly non-stationary, Nonlinear properties, are difficult to extract the fault signature of rolling bearing.Currently, for non-stationary, non-thread Property the processing method of signal of rolling bearing have a lot, mainly have Wigner distribution, empirical mode decomposition (Empirical Mode Decomposition, EMD), but have obvious limitation, example As, Wigner can produce cross term when multicomponent data processing is analyzed by distribution；Wavelet transform process is believed Number lack adaptivity；Though it is adaptive signal processing method that EMD decomposes, but EMD decomposition can cause The phenomenon such as modal overlap, end effect.

It is a kind of new self adaptation that local mean value decomposes (Local Mean Decomposition, LMD) Nonstationary random response method.LMD decomposes compared with EMD decomposition, has the integrity of higher signal Holding capacity, decreases iterations, can preferably avoid the overshoot impact on signal decomposition simultaneously. Cheng Junsheng, high and Yang Yu etc. are at " vibrate and impact " volume 28 paper that the 5th phase delivered in 2009 LMD and EMD method is carried out by " local mean value decomposes the comparative study with empirical mode decomposition " Contrast, result shows to be better than EMD method in the aspect LMD method such as end effect, iterations. Owing to the acceleration signal of rolling bearing is decomposed into some multiplicative function component (Product through LMD Function, PF) containing substantial amounts of data volume, it is impossible to directly as the master sample of classification.Entropy It is a kind of index characterizing signal complexity, it is possible to effectively reduce the dimension of characteristic vector, fully Characterize the characteristic information of signal.Sample Entropy, approximate entropy are improved by fuzzy entropy (FuzzyEn), In terms of data independence and relative uniformity more prominent.Therefore use fuzzy entropy algorithm to calculate PF to divide The fuzzy entropy of amount, as simpler fault feature vector.

Summary of the invention

The present invention seeks to overcome the defect of prior art, invention is a kind of decomposes based on local mean value and ash The rolling bearing detection method of color association, by the analysis to bearing acceleration signal, and then identifies shaft The running status held.Use based on LMD fuzzy entropy algorithm and the rolling bearing detection side of grey correlation Method, overcomes EMD decomposition and had seriously modal overlap, the phenomenon of end effect and PF component The problem that data volume cannot function as greatly characteristic vector, it is achieved effective identification of rolling bearing running status.

The technical solution used in the present invention is a kind of decomposition based on local mean value and the axis of rolling of grey correlation Holding detection method, it is characterized in that, detection method uses based on LMD fuzzy entropy algorithm and grey correlation phase In conjunction with computational methods, specifically comprising the following steps that of method

Step one: utilize acceleration transducer to gather the bearing vibration acceleration letter under running status Number, including trouble-free normal bearing and have inner ring fault, rolling element fault or the bearing of outer ring fault Vibration acceleration signal；

Step 2: the acceleration signal gathered is carried out LMD decomposition, obtains some multiplicative function PF and divide Amount and residual error；

1) all Local Extremum n of x (t) are determined_i, calculate adjacent extreme point n_iAnd n_i+1Meansigma methods m_i With envelope estimated value a_i, it may be assumed that

$m_i=\frac{n_i+n_{i+1}}{2},(i=1,2,...)---\left(1\right)$

$a_i=\frac{\left|n_i-n_{i+1}\right|}{2},(i=1,2,...)---\left(2\right)$

With straight line by meansigma methods m of all adjacent two extreme points_iCouple together, recycle moving average Line is smoothed by method, obtains local mean value function m_jk(t), (j=1,2 ...；K=1,2 ...)；Equally, By all adjacent envelope estimated values a of the bundle of lines_iCoupling together, line is put down by recycling moving average method Sliding process, obtains envelope estimation function a_jk(t), (j=1,2 ...；K=1,2 ...；K=)；

2) local mean value function m_jkT () separates from primary signal x (t), obtain function h_jk(t) be

h_jk(t)=x (t)-m_jk(t), (j=1,2 ...；K=1,2 ...；K=) (3)

3) h is used_jkT () is divided by envelope estimation function a_jkT () is to h_jkT () is demodulated, obtain FM signal s_jk(t) be

$s_{jk}\left(t\right)=\frac{h_{jk}\left(t\right)}{a_{jk}\left(t\right)},(j=1,2,...;k=1,2,...)---\left(4\right)$

Preferably s_jkT () is a pure FM signal, its local envelope function meets a_j(k+1)=1；If It is unsatisfactory for, then s_jkT () repeats above-mentioned step as primary signal, until obtaining pure FM signal s_jnT (), i.e. meets-1≤s_jn(t)≤1, its local envelope function a_j(n+1)T ()=1, has

$\left\{\begin{array}{c}{h}_{j1}\left(t\right)=x\left(t\right)-m_{j1}\left(t\right)\\ h_{j2}\left(t\right)=s_{j1}\left(t\right)-m_{j2}\left(t\right)\\ .\\ .\\ .\\ h_{jn}\left(t\right)=s_{jn}\left(t\right)-m_{jn}\left(t\right)\end{array}\right.,(j=1,2,\mathrm{...})---\left(5\right)$

$\left\{\begin{array}{c}{s}_{j1}\left(t\right)=\frac{h_{j1}\left(t\right)}{a_{j1}\left(t\right)}\\ s_{j2}\left(t\right)=\frac{h_{j2}\left(t\right)}{a_{j2}\left(t\right)}\\ .\\ .\\ .\\ s_{jn}\left(t\right)=\frac{h_{jn}\left(t\right)}{a_{jn}\left(t\right)}\end{array}\right.,(j=1,2,\mathrm{...})---\left(6\right)$

General stopping criterion for iteration is

$\underset{n\→\mathrm{\∞}}{\lim }{a}_{jn}\left(t\right)=1---\left(7\right)$

In actual applications, for avoiding too much decomposition number of times, if an amount of change Δ, meet 1-Δ≤a_jnT ()≤1+ Δ, iteration terminates；

(4) the local envelope function of generation is multiplied obtains envelope signal, i.e. the instantaneous width of PF component Value, i.e.

a_j(t)=a_j1(t)a_j2(t)…a_jn(t), (j=1,2 ...) (8)

5) by envelope signal a_j(t) and pure FM signal s_jnT () is multiplied, obtain first of primary signal PF component, i.e.

PF_j(t)=a_j(t)s_jn(t), (j=1,2 ...) (9)

PF_jT () contains highest frequency component in primary signal, be the AM/FM amplitude modulation/frequency modulation letter of a simple component Number, envelope signal a_jT () is exactly its instantaneous amplitude, its instantaneous frequency f_jT () is by pure FM signal s_jnT () is asked Go out, i.e.

$f_j\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{d\[\arccos \left(s_{jn}\right(t\left)\right)\]}{dt},(j=1,2,...)---\left(10\right)$

6) by PF₁T () separates from primary signal x (t), obtain time signal u₁(t), u₁T (), as primary signal repeat the above steps, is circulated p time, until u_pT () is a monotonic function.

$\left\{\begin{array}{c}{u}_1\left(t\right)=x\left(t\right)-{\mathrm{PF}}_1\left(t\right)\\ u_2\left(t\right)=u_1\left(t\right)-{\mathrm{PF}}_2\left(t\right)\\ .\\ .\\ .\\ u_p\left(t\right)=u_{p-1}\left(t\right)-{\mathrm{PF}}_p\left(t\right)\end{array}\right.---\left(11\right)$

Primary signal is decomposed for p PF component and a monotonic function u_p(t) sum, i.e.

$x\left(t\right)=\underset{j=1}{\overset{p}{\Σ}}{\mathrm{PF}}_i\left(t\right)+u_p\left(t\right)---\left(12\right)$

Step 3:, inner ring fault, rolling element fault normal to bearing and four kinds of states of outer ring fault are each Take 5 groups of data as master sample, calculate the fuzzy entropy of 5 groups of front 3 PF components of sample, take it Average is as canonical matrix；

Step 4: four kinds of operating modes respectively take 3 groups of data as test sample, extract its front 3 PF components Fuzzy entropy, use Grey Relation Algorithm to calculate the grey relational grade of test sample and canonical matrix, enter And carry out Fault Pattern Recognition；

The general process of grey correlation analysis:

1) suitable characteristic parameter composition state model vector is selected

X=[x (1), x (2) ..., x (k) ..., x (n)] (k=1,2 ..., n) (13)

2) state number to be classified is determined, structure standard state pattern vector or matrix

X_i=[x_i(1), x_i(2) ..., x_i(n)] (i=1,2 ..., m) (14)

3) the state model vector of bearing to be checked is extracted

Y=[y (1), y (2) ..., y (n)] (15)

4) Y Yu m mode standard vector X of state model vector of bearing to be checked is calculated respectively_iAt k point Coefficient of association

$r\[y\left(k\right),x_i\left(k\right)\]=\frac{{\mathrm{\Δ}}_{\min }+{\mathrm{\ξ\Δ}}_{\max }}{{\mathrm{\Δ}}_i\left(k\right)+{\mathrm{\ξ\Δ}}_{\max }}---\left(16\right)$

In formula: ξ is resolution ratio, it is a previously selected constant, general desirable ξ≤0.5；Δ_min It is state model vector Y and mode standard vector X_iIn the minimum absolute difference value of each element, i.e.

${\mathrm{\Δ}}_{min}=\underset{i}{min}\underset{k}{min}|y\left(k\right)-x_i\left(k\right)|,(k=1,2,...,n;i=1,2,\mathrm{...},m)---\left(17\right)$

Δ_maxFor state model vector Y and mode standard vector X_iIn the maximum absolute difference of each element, I.e.

${\mathrm{\Δ}}_{\max }=\underset{i}{\max }\underset{k}{\max }|y\left(k\right)-x_i\left(k\right)|,(k=1,2,...,n;i=1,2,\mathrm{...},m)---\left(18\right)$

Δ_iK () is state model vector Y and mode standard vector X_iThe absolute difference of middle element k, i.e.

Δ_i(k)=| y (k)-x_i(k) |, (k=1,2 ..., n；I=1,2 ..., m) (19)

5) pattern vector Y to be checked and mode standard vector X is calculated respectively_iThe degree of association

$R_i=\frac{1}{n}\underset{k=1}{\overset{n}{\Σ}}r\[y\left(k\right),x_i\left(k\right)\]---\left(20\right)$

By the descending sequence of the degree of association of test pattern vector with m mode standard vector, according to association The size of degree determines that bearing to be checked belongs to the probability of which mode standard more greatly, and infers the event of bearing Barrier degree；By above-mentioned calculating, the cohesion of multiple parameters is become parameter, i.e. a grey relational grade；Ash Population mean degree of closeness between color degree of association quantitatively characterizing pattern vector to be checked and mode standard vector, Thus may be used for state classification and the fault diagnosis of rolling bearing.

The invention have the benefit that, non-gaussian non-linear for rolling bearing fault vibration signal and The feature of non-stationary, uses the LMD decomposition method can adaptive process signal, it is possible to obtain original The time-frequency distributions of signal integrity；Use the fuzzy entropy of PF component as characteristic vector, overcome PF Component data amount is big, it is impossible to as the shortcoming of characteristic vector；This recognition methods can effectively be shaken The feature extraction of dynamic signal, exactly to normal, inner ring fault, rolling element fault and outer ring fault four Plant the identification of bearing state.

Accompanying drawing explanation

Fig. 1 is the flow chart of detection method；

Fig. 2 a) it is the time domain beamformer of bearing inner race malfunction, Fig. 2 b) it is bearing roller fault shape The time domain beamformer of state, Fig. 2 c) it is the time domain beamformer of bearing normal condition, Fig. 2 d) it is bearing outer ring Malfunction time domain beamformer.

Detailed description of the invention

With technical scheme, the enforcement of the present invention is elaborated below in conjunction with the accompanying drawings.

Fig. 1 is the flow chart of detection method, and the power of motor used by this test is 1.5KW, tests institute It is SKF6205 by the model of bearing.Bearing rotating speed is 1750r/min, and sample frequency is 12KHz, therefore Hindering a diameter of 0.1778mm, the inner ring of bearing, outer ring and rolling element fault are all artificially processed with electric spark Form.The running status of bearing is divided into normally, inner ring fault, rolling element fault and outer ring fault, and four Plant the time domain beamformer of bearing state as shown in Figure 2.Specifically comprising the following steps that of detection method

The first step utilizes acceleration transducer to gather the vibration acceleration signal of rolling bearing, including without reason Barrier normal bearing and have the bear vibration acceleration of inner ring fault, rolling element fault or outer ring fault to believe Number.

The primary signal gathered is used LMD to decompose by second step, is decomposed into some PF components and Individual residual error, each PF component is to be multiplied gained by an envelope function and a pure FM Function. PF component owing to obtaining after LMD decomposes contains substantial amounts of data, it is impossible to be directly used as event Barrier characteristic vector, the fuzzy entropy therefore using fuzzy entropy method to calculate PF component is used as characteristic vector, Extract the state characteristic vector of rolling bearing.

Normal, inner ring fault, rolling element fault and four kinds of states of outer ring faulty bearings are respectively taken by the 3rd step 5 groups of data are as master sample, and often group 1024 points of data, calculating is often organized the PF component of data, taken The meansigma methods composition canonical matrix of 5 groups.

The fuzzy entropy of front 3 the PF components of 5 groups of samples of four-step calculation, takes its average as standard Matrix, table one is the fuzzy entropy average after LMD decomposes.

Fuzzy entropy average after table one LMD decomposition

5th step uses Grey Relation Algorithm to calculate the grey relational grade of test sample and canonical matrix, enters And carry out Fault Pattern Recognition；Calculate pattern vector Y to be checked and mode standard vector X respectively_iAssociation According to the size of the degree of association, degree, by descending for degree of association sequence, determines which mark bearing to be checked belongs to The probability of quasi-mode is bigger, and infers the fault degree of bearing；The cohesion of multiple parameters is become one Parameter, i.e. grey relational grade；Grey relational grade quantitatively characterizing pattern vector to be checked and mode standard vector Between population mean degree of closeness, thus may be used for state classification and the fault diagnosis of rolling bearing.

Four kinds of states of bearing respectively take 3 groups of data as test sample, extract the mould of its front 3 PF components Stick with paste entropy, use Grey Relation Algorithm to calculate the grey relational grade of test sample and canonical matrix, Jin Erjin Row Fault Pattern Recognition.If the degree of association is the biggest, illustrates that the similarity degree of the two is the highest, thus differentiate The fault mode of sample to be tested, table two grey relational grade result of calculation.

Table two grey relational grade result of calculation

As can be seen from Table II, decompose based on local mean value and the method for grey correlation identifies well Bearing fault state, can be widely applied to the fault detect of rolling bearing.

Claims (1)

1. decompose based on local mean value and a rolling bearing detection method for grey correlation, it is characterized in that, Detection method uses the computational methods combined based on LMD fuzzy entropy algorithm and grey correlation, method Specifically comprise the following steps that

Step one: utilize acceleration transducer to gather the bearing vibration acceleration letter under running status Number, including trouble-free normal bearing and have inner ring fault, rolling element fault or the bearing of outer ring fault Vibration acceleration signal；

Step 2: the acceleration signal gathered is carried out LMD decomposition, obtains some multiplicative function PF and divide Amount and residual error；

1) all Local Extremum n of x (t) are determined_i, calculate adjacent extreme point n_iAnd n_i+1Meansigma methods m_i With envelope estimated value a_i, that is:

$m_i=\frac{n_i+n_{i+1}}{2},(i=1,2,...)---\left(1\right)$

$a_i=\frac{|n_i-n_{i+1}|}{2},(i=1,2,...)---\left(2\right)$

With straight line by meansigma methods m of all adjacent two extreme points_iCouple together, recycle moving average Line is smoothed by method, obtains local mean value function m_jk(t), (j=1,2 ...；K=1,2 ...)；Equally, By all adjacent envelope estimated values a of the bundle of lines_iCoupling together, line is put down by recycling moving average method Sliding process, obtains envelope estimation function a_jk(t), (j=1,2 ...；K=1,2 ...)；

2) local mean value function m_jkT () separates from primary signal x (t), obtain function h_jk(t) be

h_jk(t)=x (t)-m_jk(t), (j=1,2 ...；K=1,2 ...) (3)

3) h is used_jkT () is divided by envelope estimation function a_jkT () is to h_jkT () is demodulated, obtain FM signal s_jk(t) be

$s_{jk}\left(t\right)=\frac{h_{jk}\left(t\right)}{a_{jk}\left(t\right)},(j=1,2,...;k=1,2,...)---\left(4\right)$

Preferably s_jkT () is a pure FM signal, its local envelope function meets a_j(k+1)=1；If It is unsatisfactory for, then s_jkT () repeats above-mentioned step as primary signal, until obtaining pure FM signal s_jnT (), i.e. meets-1≤s_jn(t)≤1, its local envelope function a_j(n+1)T ()=1, has

$\left\{\begin{array}{c}{h}_{j1}\left(t\right)=x\left(t\right)-m_{j1}\left(t\right)\\ h_{j2}\left(t\right)=s_{j1}\left(t\right)-m_{j2}\left(t\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ h_{jn}\left(t\right)=s_{jn}\left(t\right)-m_{jn}\left(t\right)\end{array}\right.,\left(j=1,2,...\right)---\left(5\right)$

$\left\{\begin{array}{c}{s}_{j1}\left(t\right)=\frac{h_{j1}\left(t\right)}{a_{j1}\left(t\right)}\\ s_{j2}\left(t\right)=\frac{h_{j2}\left(t\right)}{a_{j2}\left(t\right)}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ s_{jn}\left(t\right)=\frac{h_{jn}\left(t\right)}{a_{jn}\left(t\right)}\end{array}\right.,\left(j=1,2,...\right)---\left(6\right)$

General stopping criterion for iteration is

$\underset{n\→\mathrm{\∞}}{\lim }{a}_{jn}\left(t\right)=1---\left(7\right)$

In actual applications, for avoiding too much decomposition number of times, if an amount of change Δ, meet 1-Δ≤a_jnT ()≤1+ Δ, iteration terminates；

(4) the local envelope function of generation is multiplied obtains envelope signal, i.e. the instantaneous width of PF component Value, i.e.

a_j(t)=a_j1(t)a_j2(t)…a_jn(t), (j=1,2 ...) (8)

5) by envelope signal a_j(t) and pure FM signal s_jnT () is multiplied, obtain first of primary signal PF component, i.e.

PF_j(t)=a_j(t)s_jn(t), (j=1,2 ...) (9)

PF_jT () contains highest frequency component in primary signal, be the AM/FM amplitude modulation/frequency modulation letter of a simple component Number, envelope signal a_jT () is exactly its instantaneous amplitude, its instantaneous frequency f_jT () is by pure FM signal s_jnT () is asked Go out, i.e.

$f_j\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{d\[\arccos \left(s_{jn}\right(t\left)\right)\]}{dt},(j=1,2,...)---\left(10\right)$

6) by PF₁T () separates from primary signal x (t), obtain time signal u₁(t), u₁T (), as primary signal repeat the above steps, is circulated p time, until u_pT () is a monotonic function；

$\left\{\begin{array}{c}{u}_1\left(t\right)=x\left(t\right)-{\mathrm{PF}}_1\left(t\right)\\ u_2\left(t\right)=u_1\left(t\right)-{\mathrm{PF}}_2\left(t\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ u_p\left(t\right)=u_{p-1}\left(t\right)-{\mathrm{PF}}_p\left(t\right)\end{array}\right.---\left(11\right)$

Primary signal is decomposed for p PF component and a monotonic function u_p(t) sum, i.e.

$x\left(t\right)=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{PF}}_i\left(t\right)+u_p\left(t\right)---\left(12\right)$

Step 3:, inner ring fault, rolling element fault normal to bearing and four kinds of states of outer ring fault are each Take 5 groups of data as master sample, calculate the fuzzy entropy of 5 groups of front 3 PF components of sample, take it Average is as canonical matrix；

Step 4: four kinds of operating modes respectively take 3 groups of data as test sample, extract its front 3 PF components Fuzzy entropy, use Grey Relation Algorithm to calculate the grey relational grade of test sample and canonical matrix, enter And carry out Fault Pattern Recognition；

The general process of grey correlation analysis:

1) suitable characteristic parameter composition state model vector is selected

X=[x (1), x (2) ..., x (k) ..., x (n)] (k=1,2 ..., n) (13)

2) state number to be classified is determined, structure standard state pattern vector or matrix

X_i=[x_i(1), x_i(2) ..., x_i(n)] (i=1,2 ..., m) (14)

3) the state model vector of bearing to be checked is extracted

Y=[y (1), y (2) ..., y (n)] (15)

4) Y Yu m mode standard vector X of state model vector of bearing to be checked is calculated respectively_iAt k point Coefficient of association

$r\[y\left(k\right),x_i\left(k\right)\]=\frac{{\mathrm{\Δ}}_{\min }+{\mathrm{\ξ\Δ}}_{\max }}{{\mathrm{\Δ}}_i\left(k\right)+{\mathrm{\ξ\Δ}}_{\max }}---\left(16\right)$

In formula: ξ is resolution ratio, it is a previously selected constant, general desirable ξ≤05 Δ_min It is state model vector Y and mode standard vector X_iIn the minimum absolute difference value of each element, i.e.

${\mathrm{\Δ}}_{\min }=\underset{i}{\min }\underset{k}{\min }|y\left(k\right)-x_i\left(k\right)|,(k=1,2,...,n;i=1,2,...,m)---\left(17\right)$

Δ_maxFor state model vector Y and mode standard vector X_iIn the maximum absolute difference of each element, I.e.

${\mathrm{\Δ}}_{max}=\underset{i}{max}\underset{k}{max}|y\left(k\right)-x_i\left(k\right)|,(k=1,2,...,n;i=1,2,...,m)---\left(18\right)$

Δ_iK () is state model vector Y and mode standard vector X_iThe absolute difference of middle element k, i.e.

Δ_i(k)=| y (k)-x_i(k) |, (k=1,2 ..., n；I=1,2 ..., m) (19)

5) pattern vector Y to be checked and mode standard vector X is calculated respectively_iThe degree of association

$R_i=\frac{1}{n}\underset{k=1}{\overset{n}{\Σ}}r\[y\left(k\right),x_i\left(k\right)\]---\left(20\right)$

By the descending sequence of the degree of association of test pattern vector with m mode standard vector, according to association The size of degree determines that bearing to be checked belongs to the probability of which mode standard more greatly, and infers the event of bearing Barrier degree.