CN109635306A - Rotary machinery fault diagnosis method based on wavelet decomposition and spectrum kurtosis - Google Patents
Rotary machinery fault diagnosis method based on wavelet decomposition and spectrum kurtosis Download PDFInfo
- Publication number
- CN109635306A CN109635306A CN201811110099.5A CN201811110099A CN109635306A CN 109635306 A CN109635306 A CN 109635306A CN 201811110099 A CN201811110099 A CN 201811110099A CN 109635306 A CN109635306 A CN 109635306A
- Authority
- CN
- China
- Prior art keywords
- signal
- wavelet decomposition
- data
- layer
- envelope
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Withdrawn
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/17—Mechanical parametric or variational design
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
- G01M13/00—Testing of machine parts
- G01M13/04—Bearings
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/12—Classification; Matching
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Geometry (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Computational Mathematics (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
Abstract
The invention discloses the rotary machinery fault diagnosis methods based on wavelet decomposition and spectrum kurtosis, this method decomposes original signal first with wavelet-decomposing method, then noise component(s) and trend term in decomposition result are excluded using the rearrangement of data and substitution operation, then filtered signal for the first time is analyzed using spectrum kurtosis method again, obtain the centre frequency and bandwidth of optimal filter, then second of filtering is carried out again to filtered signal for the first time using the filter, then Envelope Analysis is carried out to second of filtered signal using cubic spline iteration smoothed envelope analysis method, the fault type of rotating machinery is finally determined according to envelope spectrum.The present invention is suitable for the complicated rotating machinery fault signal of processing, can accurately determine the fault type of rotating machinery, has good noise immunity and robustness, is convenient for engineer application.
Description
The present invention is application number 2016104924296, and on June 29 2016 applying date, " one kind is based on small echo to denomination of invention
The divisional application of the envelope Analysis Method of decomposition and spectrum kurtosis ".
Technical field
The present invention relates to condition monitoring for rotating machinery and fault diagnosis field, and in particular to based on wavelet decomposition and spectrum kurtosis
Rotary machinery fault diagnosis method.
Background technique
Envelope Analysis technology is widely used in the fault diagnosis of gear and rolling bearing.Existing Envelope Analysis technology has
Three defects below: 1. existing Envelope Analysis technology either directly analyzes original signal, or only to original
Signal analyzed again after simply filtering, therefore existing method is easy to be done by noise, trend and other ingredients
It disturbs, it is lower so as to cause the analysis precision of the prior art;2. existing Envelope Analysis technology is based on being converted by Hilbert,
And it must be the narrow band signal of simple component that Hilbert transformation, which requires analyzed signal, otherwise the frequency modulating section of signal will
The amplitude envelope analysis of signal is polluted as a result, still signal to be analyzed at present does not meet the item of simple component and narrowband strictly
Part will lead to the prior art because precision is not high in this way and be easy to appear erroneous judgement problem;3. the envelope spectrum obtained by conventional method
There is end effects.
Summary of the invention
The problem to be solved in the present invention is against the above deficiency, to propose the rotating machinery event based on wavelet decomposition and spectrum kurtosis
Hinder diagnostic method, after envelope Analysis Method of the invention, has analysis result accuracy and accuracy high, and can be accurately
The advantages of detecting rotating machinery fault type.
In order to solve the above technical problems, the technical solution adopted by the present invention is as follows: based on wavelet decomposition and spectrum kurtosis rotation
Turn mechanical failure diagnostic method, which comprises the following steps:
), step 1: measuring the vibration signal x(k of rotating machinery with sample frequency fs using acceleration transducer k=1,2 ...,
N, N are the length of sampled signal;
Step 2: signal x(k) being resolved by the sum of n component and a trend term using wavelet decomposition algorithm, i.e.,, wherein ci(k) i-th of the component obtained by wavelet decomposition algorithm, r are representedn
(k) trend term obtained by wavelet decomposition algorithm is represented;N=15 in this example;
Step 3: to ci(k) reordering operations and substitution operation are executed, it is rearranged to operate obtained data ci shuffle(k) it indicates,
Data c is obtained after substitution operationi FTran(k) it indicates;
Step 4: to ci(k), ci shuffle(k) and ci FTran(k) multi-fractal is executed respectively remove trend fluction analysis
(Multifractal Detrended Fluctuation Analysis, MFDFA), obtains generalized Hurst index curve, ci
(k) generalized Hurst index curve Hi(q) it indicates;ci shuffle(k) generalized Hurst index curve Hi shuffle(q) table
Show;ci FTran(k) generalized Hurst index curve Hi FTran(q) it indicates;
Step 5: if Hi(q) and Hi shuffle(q) or Hi(q) and Hi FTran(q) relative error between is less than 5% or Hi
(q), Hi shuffle(q) and Hi FTran(q) three does not change with q, then abandons corresponding ci(k) component;
Step 6: to remaining ci(k) component is summed, and by this and to be denoted as signal rearranged and substitute filtered result xf1(k);
Step 7: to xf1(k) spectrum kurtosis analysis is executed, centre frequency f corresponding to signal kurtosis maximum is found out0And bandwidth B;
Step 8: according to centre frequency f0With bandwidth B to xf1(k) bandpass filtering is carried out, x is obtainedf2(k);
Step 9: to signal xf2(k) analysis of cubic spline iteration smoothed envelope is executed, signal envelope eov(k is obtained);
Step 10: to obtained signal envelope eov(k it) executes discrete Fourier transform and obtains envelope spectrum, according to envelope spectrum signature
Frequency judges the fault type of machine;
MFDFA method in the step 4 the following steps are included:
1) x is constructedy(k) profileY(i), k=1,2 ..., N:
,
,
xy(k) c in the step 4 is representedi(k) or ci shuffle(k) or ci FTran(k);
2) by signal profileY(i) be divided into it is nonoverlappingN s Segment length issData, due to data lengthNIt generally can not divide exactlys,
So the remaining one piece of data of meeting cannot utilize;
In order to make full use of the length of data, then from the opposite direction of data 2 are obtained with identical length segmentation, such oneN s Section
Data;
3) it is fitted the polynomial trend of every segment data using least square method, then calculates the variance of every segment data:
;
y v (i) it is the of fittingvThe trend of segment data, if the polynomial trend of fitting ismRank then remembers that this goes the trend process to be
(MF-) DFAm;
4) the is calculatedqThe average value of rank wave function:
;
If 5) xy(k) there are self-similarity characteristics, thenqThe average value of rank wave functionF q (s) and lengthsBetween there are power laws
Relationship:
F q (s)~s H(q);
WhenqWhen=0, the formula in step 4) dissipates, at this momentH(0) it is determined by logarithmic mean process defined in following formula:
;
6) to the formula both sides in step 5) take logarithm can obtain ln [F q (s)]=H(q)ln(s)+c,cFor constant, it is possible thereby to obtain
Obtain the slope of straight lineH(q)。
Further, in the step 2 wavelet decomposition algorithm the following steps are included:
1) the 1st layer of wavelet decomposition: to x (t), t=1,2 ..., N, wavelet decomposition is executed
In above formula, g [k] and h [k] are respectively low-pass filter and high-pass filter, x1,L[n] and x1,H[n] is respectively signal x
(t) after low-pass filtered device and high pass filter filters as a result, x1,H[n] is known as the 1st layer of coefficient of wavelet decomposition, x1,L[n] is
The signal to be decomposed that 1st layer of wavelet decomposition obtains;In the 1st layer of wavelet decomposition, x (t) represents x (k) in the step 2;
2) the 2nd layer of wavelet decomposition: x1,L[n] is used as new data, executes above-mentioned steps 1 again)
x2,H[n] is known as the 2nd layer of coefficient of wavelet decomposition, x2,L[n] is the signal to be decomposed that the 2nd layer of wavelet decomposition obtains;
3) pth time wavelet decomposition: in pth time wavelet decomposition, xp-1,L[n] is used as new data, executes above-mentioned steps again
1) original signal finally, is decomposed into following form:, ci(t) it represents
I-th layer of coefficient of wavelet decomposition, rp(t) pth layer signal to be decomposed is represented;In pth time wavelet decomposition,;
xp-1,L[n] represents the signal to be decomposed that (p-1) layer wavelet decomposition obtains.
Further, in the step 3 data rearrangement operation the following steps are included:
Upset component c at randomi(k) put in order.
Further, in the step 3 data substitution operation the following steps are included:
1) to component ci(k) discrete Fourier transform is executed, component c is obtainedi(k) phase;
2) it is located at the pseudo- independent same distribution number in the section (- π, π) with one group to replace component ci(k) original phase;
3) inverse discrete Fourier transform is executed to the frequency domain data after phase substitutes and obtains data ci IFFT(k), it seeks counting
According to ci IFFT(k) real part.
Further, kurtosis method is composed in the step 7:
1) constructing a cutoff frequency is fcThe low-pass filter h (n) of=0.125+ ε;ε=0.175;
It 2) is the quasi- low-pass filter h of [0,0.25] based on h (n) construction passband0(n) and passband is [0.25,0.5]
Quasi- high-pass filter h1(n),
;
3) signal ci k(n) through h0(n)、 h1(n) it filters and resolves into low frequency part c after down-sampled2i k+1(n) and high frequency section
c2i+1 k+1(n), the down-sampled factor is 2, then the shaping filter tree after successive ignition filters, kth layer have 2kA frequency band, wherein
ci k(n) output signal of i-th of filter in expression filter tree on kth layer, i=0 ..., 2k- 1,0≤k≤K-1, K generation
The number of plies of table filter tree;c0(n) x in the step 7 is representedf1(k);
4) the centre frequency f of i-th of filter in decomposition tree on kth layerkiAnd bandwidth BkRespectively
;
5) each filter results c is calculatedi k(n) i=0,…, 2k- 1 kurtosis;
6) all spectrum kurtosis are summarized, obtains the total spectrum kurtosis of signal.
Further, the cubic spline iteration smoothed envelope analysis method in the step 9 the following steps are included:
1) signal is calculatedz(k) Jue Dui Zhi ∣z(k) the local extremum of ∣;In the 1st iteration,z(k) represent in the step 9
xf2(k);
2) envelope eov is obtained using cubic spline interpolation Local Extremum1(k);
3) rightz(k) be normalized to obtain;
4) the 2nd iteration:z 1(k) it is re-used as new data, repeat above-mentioned steps 1) ~ 3), it obtains;
5) i-th iteration:z i-1(k) it is re-used as new data, repeat above-mentioned steps 1) ~ 3), it obtains;
If 6)nWhat secondary iteration obtainedz n (k) amplitude be less than or equal to 1, then iterative process stop, finally obtaining signal
Envelope eov(k) be。
The invention adopts the above technical scheme, compared with prior art, the invention has the following advantages that
1) original signal is decomposed using wavelet decomposition, then excludes therein make an uproar using the rearrangement of data and substitution operation
Sound and trend component, the only useful component in stick signal component, so as to avoid noise and trend component to Envelope Analysis
As a result influence, analysis result accuracy and accuracy are high.
2) signal envelope and frequency modulating section are kept completely separate using cubic spline iteration smoothed envelope analysis method, energy
Influence of the frequency modulating section to signal envelope analysis result is enough avoided, to improve the precision of Envelope Analysis.
3) fault type of rotating machinery can be accurately detected.
4) there are end effects for the envelope spectrum obtained by conventional method, and can be avoided by the envelope spectrum that the present invention obtains
End effect.
The present invention will be further described with reference to the accompanying drawings and examples.
Detailed description of the invention
Attached drawing 1 is the flow chart of the method for the present invention in the embodiment of the present invention;
Attached drawing 2 is the signal for carrying out preliminary exposition in the embodiment of the present invention to signal using low-pass filter and high-pass filter
Figure;
Attached drawing 3 is the schematic diagram for quickly calculating spectrum kurtosis in the embodiment of the present invention using tree-shaped filter construction;
Attached drawing 4 is the bearing vibration signal in the embodiment of the present invention with inner ring failure;
Attached drawing 5 is the analysis in the embodiment of the present invention using traditional envelope Analysis Method to inner ring faulty bearing vibration signal
As a result;
Attached drawing 6 is the present invention in the embodiment of the present invention to the analysis result of inner ring faulty bearing vibration signal;
Attached drawing 7 is the bearing vibration signal in the embodiment of the present invention with outer ring failure;
Attached drawing 8 is the analysis in the embodiment of the present invention using traditional envelope Analysis Method to outer ring faulty bearing vibration signal
As a result;
Attached drawing 9 is the present invention in the embodiment of the present invention to the analysis result of outer ring faulty bearing vibration signal.
Specific embodiment
Embodiment, as shown in Figure 1, Figure 2, Figure 3 shows, the rotary machinery fault diagnosis method based on wavelet decomposition and spectrum kurtosis,
The following steps are included:
Step 1: measuring the vibration signal x(k of rotating machinery with sample frequency fs using acceleration transducer), (k=1,2,
..., N), N is the length of sampled signal;
Step 2: signal x(k) being resolved by the sum of n component and a trend term using wavelet decomposition algorithm, i.e.,, wherein ci(k) i-th point obtained by wavelet decomposition algorithm is represented
Amount, rn(k) trend term obtained by wavelet decomposition algorithm is represented;N=15 in this example;
Step 3: to ci(k) reordering operations and substitution operation are executed, it is rearranged to operate obtained data ci shuffle(k) it indicates,
Data c is obtained after substitution operationi FTran(k) it indicates;
Step 4: to ci(k), ci shuffle(k) and ci FTran(k) multi-fractal is executed respectively remove trend fluction analysis
(Multifractal Detrended Fluctuation Analysis, MFDFA), obtains generalized Hurst index curve, ci
(k) generalized Hurst index curve Hi(q) it indicates;ci shuffle(k) generalized Hurst index curve Hi shuffle(q) table
Show;ci FTran(k) generalized Hurst index curve Hi FTran(q) it indicates;
Step 5: if Hi(q) and Hi shuffle(q) or Hi(q) and Hi FTran(q) relative error between is less than 5% or Hi
(q), Hi shuffle(q) and Hi FTran(q) three does not change with q, then abandons corresponding ci(k) component;
Step 6: to remaining ci(k) component is summed, and by this and to be denoted as signal rearranged and substitute filtered result xf1(k);
Step 7: to xf1(k) spectrum kurtosis analysis is executed, centre frequency f corresponding to signal kurtosis maximum is found out0And bandwidth B;
Step 8: according to centre frequency f0With bandwidth B to xf1(k) bandpass filtering is carried out, x is obtainedf2(k);
Step 9: to signal xf2(k) analysis of cubic spline iteration smoothed envelope is executed, signal envelope eov(k is obtained);
Step 10: to obtained signal envelope eov(k it) executes discrete Fourier transform and obtains envelope spectrum, according to envelope spectrum signature
Frequency judges the fault type of machine.
Wavelet decomposition algorithm in step 2 the following steps are included:
1) the 1st layer of wavelet decomposition: to x (t), (t=1,2 ..., N) executes wavelet decomposition
In above formula, g [k] and h [k] are respectively low-pass filter and high-pass filter, x1,L[n] and x1,H[n] is respectively signal x
(t) after low-pass filtered device and high pass filter filters as a result, x1,H[n] is known as the 1st layer of coefficient of wavelet decomposition, x1,L[n] is
The signal to be decomposed that 1st layer of wavelet decomposition obtains;In the 1st layer of wavelet decomposition, x (t) represents x (k) in the step 2;
2) the 2nd layer of wavelet decomposition: x1,L[n] is used as new data, executes above-mentioned steps 1 again)
x2,H[n] is known as the 2nd layer of coefficient of wavelet decomposition, x2,L[n] is the signal to be decomposed that the 2nd layer of wavelet decomposition obtains;
3) pth time wavelet decomposition: in pth time wavelet decomposition, xp-1,L[n] is used as new data, executes above-mentioned steps again
1) original signal finally, is decomposed into following form:, ci(t) it represents
I-th layer of coefficient of wavelet decomposition, rp(t) pth layer signal to be decomposed is represented;;In pth time wavelet decomposition,;
xp-1,L[n] represents the signal to be decomposed that (p-1) layer wavelet decomposition obtains.
In step 3 data rearrangement operation the following steps are included:
Upset component c at randomi(k) put in order.
In step 3 data substitution operation the following steps are included:
1) to component ci(k) discrete Fourier transform is executed, component c is obtainedi(k) phase;
2) it is located at the pseudo- independent same distribution number in the section (- π, π) with one group to replace component ci(k) original phase;
3) inverse discrete Fourier transform is executed to the frequency domain data after phase substitutes and obtains data ci IFFT(k), it seeks counting
According to ci IFFT(k) real part.
MFDFA method in the step 4 the following steps are included:
1) x is constructedy(k) profileY(i), k=1,2 ..., N:
,
,
xy(k) c in the step 4 is representedi(k) or ci shuffle(k) or ci FTran(k);
2) by signal profileY(i) be divided into it is nonoverlappingN s Segment length issData, due to data lengthNIt generally can not divide exactlys,
So the remaining one piece of data of meeting cannot utilize;
In order to make full use of the length of data, then from the opposite direction of data 2 are obtained with identical length segmentation, such oneN s Section
Data;
3) it is fitted the polynomial trend of every segment data using least square method, then calculates the variance of every segment data:
;
y v (i) it is the of fittingvThe trend of segment data, if the polynomial trend of fitting ismRank then remembers that this goes the trend process to be
(MF-) DFAm;
4) the is calculatedqThe average value of rank wave function:
;
If 5) xy(k) there are self-similarity characteristics, thenqThe average value of rank wave functionF q (s) and lengthsBetween there are power laws
Relationship:
F q (s)~s H(q);
WhenqWhen=0, the formula in step 4) dissipates, at this momentH(0) it is determined by logarithmic mean process defined in following formula:
;
6) to the formula both sides in step 5) take logarithm can obtain ln [F q (s)]=H(q)ln(s)+c,cFor constant, it is possible thereby to obtain
Obtain the slope of straight lineH(q)。
Kurtosis method is composed in step 7:
1) constructing a cutoff frequency is fcThe low-pass filter h (n) of=0.125+ ε;ε > 0, ε=0.175;
It 2) is the quasi- low-pass filter h of [0,0.25] based on h (n) construction passband0(n) and passband is [0.25,0.5]
Quasi- high-pass filter h1(n),
;
3) signal ci k(n) through h0(n)、 h1(n) it filters and resolves into low frequency part c after down-sampled2i k+1(n) and high frequency section
c2i+1 k+1(n), the down-sampled factor is 2, then the shaping filter tree after successive ignition filters, kth layer have 2kA frequency band, wherein
ci k(n) output signal of i-th of filter in expression filter tree on kth layer, i=0 ..., 2k- 1,0≤k≤K-1, this example
Middle K=8;c0(n) x in the step 7 is representedf1(k);
4) the centre frequency f of i-th of filter in decomposition tree on kth layerkiAnd bandwidth BkRespectively
;
5) each filter results c is calculatedi k(n)( i=0,…, 2k- 1) kurtosis;
6) all spectrum kurtosis are summarized, obtains the total spectrum kurtosis of signal.
Cubic spline iteration smoothed envelope analysis method in step 9 the following steps are included:
1) signal is calculatedz(k) Jue Dui Zhi ∣z(k) the local extremum of ∣;In the 1st iteration,z(k) represent in the step 9
xf2(k);
2) envelope eov is obtained using cubic spline interpolation Local Extremum1(k);
3) rightz(k) be normalized to obtain;
4) the 2nd iteration:z 1(k) it is re-used as new data, repeat above-mentioned steps 1) ~ 3), it obtains;
5) i-th iteration:z i-1(k) it is re-used as new data, repeat above-mentioned steps 1) ~ 3), it obtains;
If 6)nWhat secondary iteration obtainedz n (k) amplitude be less than or equal to 1, then iterative process stop, finally obtaining signal
Envelope is。
Test 1, tests the performance of algorithm of the present invention using the bearing vibration data with inner ring failure
Card.
Testing bearing used is 6205-2RS JEM SKF, using electric discharge machining method on bearing inner race working depth
The groove for being 0.3556mm for 0.2794mm, width simulates bearing inner race failure, this experiment load is about 0.7457kW, driving
It is about 29.5Hz that motor, which turns frequency, and bearing inner race fault characteristic frequency is about 160Hz, sample frequency 4.8KHz, when signal sampling
A length of 1s.
Collected inner ring fault-signal is as shown in Figure 4.
Signal shown in Fig. 4 is analyzed using traditional envelope Analysis Method first, obtained analysis result such as Fig. 5
It is shown.From fig. 5, it can be seen that the fault signature of bearing is blanked completely, therefore traditional envelope Analysis Method cannot be effectively
Extract the fault signature of bearing;In addition, this explanation is by passing from fig. 5, it can be seen that the left end point of envelope spectrum is there is abnormal high level
There is end effects for the envelope spectrum that system method obtains.
Signal shown in Fig. 4 is analyzed using method proposed by the invention, obtained analysis result such as Fig. 6 institute
Show.From fig. 6, it can be seen that spectral line corresponding to 160Hz and 320Hz is apparently higher than other spectral lines, the two frequencies are respectively corresponded
1 frequency multiplication and 2 frequencys multiplication of bearing inner race fault characteristic frequency may determine that bearing has inner ring failure accordingly;It can from Fig. 6
Out, there is no end effect by the envelope spectrum that the present invention obtains.
Show that the present invention is capable of reliable recognition in the case where load and failure dimensional depth are constant through many experiments
Minimum inner ring failure dimension width is about 0.23mm, and conventional method is capable of the minimum inner ring failure dimension width of reliable recognition about
For 0.53mm, precision improves 56.6%.
Test 2, tests the performance of algorithm of the present invention using the bearing vibration data with outer ring failure
Card.
Testing bearing used is 6205-2RS JEM SKF, using electric discharge machining method on bearing outer ring working depth
The groove for being 0.5334mm for 0.2794mm, width simulates bearing outer ring failure, this experiment load is about 2.237 kW, driving
It is about 28.7Hz that motor, which turns frequency, and bearing outer ring fault characteristic frequency is about 103Hz, sample frequency 4.8KHz, when signal sampling
A length of 1s.
Collected outer ring fault-signal is as shown in Figure 7.
Signal shown in Fig. 7 is analyzed using traditional envelope Analysis Method first, obtained analysis result such as Fig. 8
It is shown.From figure 8, it is seen that the fault signature of bearing is blanked completely, therefore traditional envelope Analysis Method cannot be effectively
Extract the fault signature of bearing;In addition, this explanation is by passing from figure 8, it is seen that the left end point of envelope spectrum is there is abnormal high level
There is end effects for the envelope spectrum that system method obtains.
Signal shown in Fig. 7 is analyzed using method proposed by the invention, obtained analysis result such as Fig. 9 institute
Show.From fig. 9, it can be seen that spectral line corresponding to 103Hz and 206Hz is apparently higher than other spectral lines, the two frequencies are respectively corresponded
1 frequency multiplication and 2 frequencys multiplication of bearing outer ring fault characteristic frequency may determine that bearing has outer ring failure accordingly;It can from Fig. 9
Out, there is no end effect by the envelope spectrum that the present invention obtains.
Show that the present invention is capable of reliable recognition in the case where load and failure dimensional depth are constant through many experiments
Minimum outer ring failure dimension width is about 0.33mm, and conventional method is capable of the minimum outer ring failure dimension width of reliable recognition about
For 0.68mm, precision improves 51.5%.
According to test result, think after analysis:
1) traditional envelope Analysis Method directly carries out Envelope Analysis to original signal, or to after merely through simple process
Original signal carries out Envelope Analysis, and different from traditional envelope Analysis Method, the invention firstly uses wavelet decompositions to original letter
It number is decomposed, then excludes noise and trend component therein using the rearrangement of data and substitution operation, only stick signal
Useful component in component, the influence so as to avoid noise and trend component to Envelope Analysis result, improve accuracy and
Accuracy.
2) based on traditional envelope Analysis Method is converted by Hilbert, and Hilbert transformation requires analyzed letter
Number must be the narrow band signal of simple component, otherwise the frequency modulating section of signal will pollute signal Envelope Analysis as a result, but
It is the condition that signal to be analyzed does not meet simple component and narrowband strictly at present, will lead to the prior art in this way because of precision not
High and be easy to appear erroneous judgement problem, different from traditional envelope Analysis Method, the present invention utilizes cubic spline iteration smoothed envelope to divide
Signal envelope and frequency modulating section are kept completely separate by analysis method, be can be avoided frequency modulating section and are analyzed result to signal envelope
Influence, to improve the precision of Envelope Analysis.
3) fault type of rotating machinery can be accurately detected.
4) there are end effects for the envelope spectrum obtained by conventional method, and can be avoided by the envelope spectrum that the present invention obtains
End effect.
5) each step effect:
The 1) step: acquisition vibration signal;
Original signal: being resolved into the form of different component sums by the 2) step, and some of them component corresponds to noise and trend term, some
Component corresponds to useful signal;
The 3) ~ 5) step: reordering operations are executed to the signal that above-mentioned decomposition obtains and substitution operates, reject noise component(s) therein and
Trend term only retains useful signal;
The 6) step: remaining useful signal is summed, by this and as signal it is rearranged and substitute filtered result xf1(k);
7) step: to filtered signal xf1(k) spectrum kurtosis analysis is executed, corresponding center frequency at signal maximum kurtosis is found out
Rate f0And bandwidth B;
8) step: according to centre frequency f0With bandwidth B to xf1(k) bandpass filtering is carried out, signal x is obtainedf2(k);
The 9) step: signal x is calculatedf2(k) envelope eov (k);
The 10) step: discrete Fourier transform is executed to eov (k) and obtains envelope spectrum, the failure classes of bearing are judged according to envelope spectrum
Type.
One skilled in the art would recognize that the above specific embodiments are only exemplary, it is to make ability
Field technique personnel can better understand the content of present invention, should not be construed as limiting the scope of protection of the present invention, as long as
Technical solution improvements introduced according to the present invention each falls within protection scope of the present invention.
Claims (6)
1. the rotary machinery fault diagnosis method based on wavelet decomposition and spectrum kurtosis, which comprises the following steps:
), step 1: measuring the vibration signal x(k of rotating machinery with sample frequency fs using acceleration transducer k=1,2 ...,
N, N are the length of sampled signal;
Step 2: signal x(k) being resolved by the sum of n component and a trend term using wavelet decomposition algorithm, i.e.,, wherein ci(k) i-th of the component obtained by wavelet decomposition algorithm, r are representedn(k)
Represent the trend term obtained by wavelet decomposition algorithm;N=15 in this example;
Step 3: to ci(k) reordering operations and substitution operation are executed, it is rearranged to operate obtained data ci shuffle(k) it indicates,
Data c is obtained after substitution operationi FTran(k) it indicates;
Step 4: to ci(k), ci shuffle(k) and ci FTran(k) multi-fractal is executed respectively remove trend fluction analysis
(Multifractal Detrended Fluctuation Analysis, MFDFA), obtains generalized Hurst index curve, ci
(k) generalized Hurst index curve Hi(q) it indicates;ci shuffle(k) generalized Hurst index curve Hi shuffle(q) table
Show;ci FTran(k) generalized Hurst index curve Hi FTran(q) it indicates;
Step 5: if Hi(q) and Hi shuffle(q) or Hi(q) and Hi FTran(q) relative error between is less than 5% or Hi
(q), Hi shuffle(q) and Hi FTran(q) three does not change with q, then abandons corresponding ci(k) component;
Step 6: to remaining ci(k) component is summed, and by this and to be denoted as signal rearranged and substitute filtered result xf1(k);
Step 7: to xf1(k) spectrum kurtosis analysis is executed, centre frequency f corresponding to signal kurtosis maximum is found out0And bandwidth B;
Step 8: according to centre frequency f0With bandwidth B to xf1(k) bandpass filtering is carried out, x is obtainedf2(k);
Step 9: to signal xf2(k) analysis of cubic spline iteration smoothed envelope is executed, signal envelope eov(k is obtained);
Step 10: to obtained signal envelope eov(k it) executes discrete Fourier transform and obtains envelope spectrum, according to envelope spectrum signature
Frequency judges the fault type of machine;
MFDFA method in the step 4 the following steps are included:
1) x is constructedy(k) profileY(i), k=1,2 ..., N:
,
,
xy(k) c in the step 4 is representedi(k) or ci shuffle(k) or ci FTran(k);
2) by signal profileY(i) be divided into it is nonoverlappingN s Segment length issData, due to data lengthNIt generally can not divide exactlys,
So the remaining one piece of data of meeting cannot utilize;
In order to make full use of the length of data, then from the opposite direction of data 2 are obtained with identical length segmentation, such oneN s Section
Data;
3) it is fitted the polynomial trend of every segment data using least square method, then calculates the variance of every segment data:
;
y v (i) it is the of fittingvThe trend of segment data, if the polynomial trend of fitting ismRank then remembers that this goes the trend process to be
(MF-) DFAm;
4) the is calculatedqThe average value of rank wave function:
;
If 5) xy(k) there are self-similarity characteristics, thenqThe average value of rank wave functionF q (s) and lengthsBetween there are power laws
Relationship:
F q (s)~s H(q);
WhenqWhen=0, the formula in step 4) dissipates, at this momentH(0) it is determined by logarithmic mean process defined in following formula:
;
6) to the formula both sides in step 5) take logarithm can obtain ln [F q (s)]=H(q)ln(s)+c,cFor constant, it is possible thereby to obtain
The slope of straight lineH(q)。
2. the rotary machinery fault diagnosis method according to claim 1 based on wavelet decomposition and spectrum kurtosis, feature exist
In, wavelet decomposition algorithm in the step 2 the following steps are included:
1) the 1st layer of wavelet decomposition: to x (t), t=1,2 ..., N, wavelet decomposition is executed
In above formula, g [k] and h [k] are respectively low-pass filter and high-pass filter, x1,L[n] and x1,H[n] is respectively signal x
(t) after low-pass filtered device and high pass filter filters as a result, x1,H[n] is known as the 1st layer of coefficient of wavelet decomposition, x1,L[n] is
The signal to be decomposed that 1st layer of wavelet decomposition obtains;In the 1st layer of wavelet decomposition, x (t) represents x (k) in the step 2;
2) the 2nd layer of wavelet decomposition: x1,L[n] is used as new data, executes above-mentioned steps 1 again)
x2,H[n] is known as the 2nd layer of coefficient of wavelet decomposition, x2,L[n] is the signal to be decomposed that the 2nd layer of wavelet decomposition obtains;
3) pth time wavelet decomposition: in pth time wavelet decomposition, xp-1,L[n] is used as new data, executes above-mentioned steps again
1) original signal finally, is decomposed into following form:, ci(t) it represents
I-th layer of coefficient of wavelet decomposition, rp(t) pth layer signal to be decomposed is represented;In pth time wavelet decomposition,;
xp-1,L[n] represents the signal to be decomposed that (p-1) layer wavelet decomposition obtains.
3. the rotary machinery fault diagnosis method according to claim 1 based on wavelet decomposition and spectrum kurtosis, feature exist
In, in the step 3 data rearrangement operation the following steps are included:
Upset component c at randomi(k) put in order.
4. the rotary machinery fault diagnosis method according to claim 1 based on wavelet decomposition and spectrum kurtosis, feature exist
In: in the step 3 data substitution operation the following steps are included:
1) to component ci(k) discrete Fourier transform is executed, component c is obtainedi(k) phase;
2) it is located at the pseudo- independent same distribution number in the section (- π, π) with one group to replace component ci(k) original phase;
3) inverse discrete Fourier transform is executed to the frequency domain data after phase substitutes and obtains data ci IFFT(k), data are sought
ci IFFT(k) real part.
5. the rotary machinery fault diagnosis method according to claim 1 based on wavelet decomposition and spectrum kurtosis, feature exist
In: kurtosis method is composed in the step 7:
1) constructing a cutoff frequency is fcThe low-pass filter h (n) of=0.125+ ε;ε=0.175;
It 2) is the quasi- low-pass filter h of [0,0.25] based on h (n) construction passband0(n) and passband is [0.25,0.5]
Quasi- high-pass filter h1(n),
;
3) signal ci k(n) through h0(n)、 h1(n) it filters and resolves into low frequency part c after down-sampled2i k+1(n) and high frequency section c2i +1 k+1(n), the down-sampled factor is 2, then the shaping filter tree after successive ignition filters, kth layer have 2kA frequency band, wherein ci k
(n) output signal of i-th of filter in expression filter tree on kth layer, i=0 ..., 2k- 1,0≤k≤K-1, K are represented
The number of plies of filter tree;c0(n) x in the step 7 is representedf1(k);
4) the centre frequency f of i-th of filter in decomposition tree on kth layerkiAnd bandwidth BkRespectively
;
5) each filter results c is calculatedi k(n) i=0,…, 2k- 1 kurtosis;
6) all spectrum kurtosis are summarized, obtains the total spectrum kurtosis of signal.
6. the rotary machinery fault diagnosis method according to claim 1 based on wavelet decomposition and spectrum kurtosis, feature exist
In, cubic spline iteration smoothed envelope analysis method in the step 9 the following steps are included:
1) signal is calculatedz(k) Jue Dui Zhi ∣z(k) the local extremum of ∣;In the 1st iteration,z(k) represent in the step 9
xf2(k);
2) envelope eov is obtained using cubic spline interpolation Local Extremum1(k);
3) rightz(k) be normalized to obtain;
4) the 2nd iteration:z 1(k) it is re-used as new data, repeat above-mentioned steps 1) ~ 3), it obtains;
5) i-th iteration:z i-1(k) it is re-used as new data, repeat above-mentioned steps 1) ~ 3), it obtains;
If 6)nWhat secondary iteration obtainedz n (k) amplitude be less than or equal to 1, then iterative process stop, finally obtaining signal packet
Network eov(k) be。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811110099.5A CN109635306A (en) | 2016-06-29 | 2016-06-29 | Rotary machinery fault diagnosis method based on wavelet decomposition and spectrum kurtosis |
Applications Claiming Priority (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610492429.6A CN106096200B (en) | 2016-06-29 | 2016-06-29 | A kind of envelope Analysis Method based on wavelet decomposition and spectrum kurtosis |
CN201811110099.5A CN109635306A (en) | 2016-06-29 | 2016-06-29 | Rotary machinery fault diagnosis method based on wavelet decomposition and spectrum kurtosis |
Related Parent Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610492429.6A Division CN106096200B (en) | 2016-06-29 | 2016-06-29 | A kind of envelope Analysis Method based on wavelet decomposition and spectrum kurtosis |
Publications (1)
Publication Number | Publication Date |
---|---|
CN109635306A true CN109635306A (en) | 2019-04-16 |
Family
ID=57215285
Family Applications (2)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610492429.6A Active CN106096200B (en) | 2016-06-29 | 2016-06-29 | A kind of envelope Analysis Method based on wavelet decomposition and spectrum kurtosis |
CN201811110099.5A Withdrawn CN109635306A (en) | 2016-06-29 | 2016-06-29 | Rotary machinery fault diagnosis method based on wavelet decomposition and spectrum kurtosis |
Family Applications Before (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610492429.6A Active CN106096200B (en) | 2016-06-29 | 2016-06-29 | A kind of envelope Analysis Method based on wavelet decomposition and spectrum kurtosis |
Country Status (1)
Country | Link |
---|---|
CN (2) | CN106096200B (en) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110929789A (en) * | 2019-11-22 | 2020-03-27 | 北京理工大学 | Liver tumor automatic classification method and device based on multi-stage CT image analysis |
CN111769810A (en) * | 2020-06-29 | 2020-10-13 | 浙江大学 | Fluid mechanical modulation frequency extraction method based on energy kurtosis spectrum |
CN112101174A (en) * | 2020-09-09 | 2020-12-18 | 洛阳师范学院 | LOF-Kurtogram-based mechanical fault diagnosis method |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107220402B (en) * | 2017-04-14 | 2020-11-13 | 桂林理工大学 | Aluminum liquid interface simulation method |
CN112595514A (en) * | 2020-11-26 | 2021-04-02 | 上海航天控制技术研究所 | High-speed bearing assembly vibration signal noise reduction processing method |
CN113092114B (en) * | 2021-04-08 | 2024-05-10 | 陕西科技大学 | Bearing fault diagnosis method, device and storage medium |
Family Cites Families (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20130096848A1 (en) * | 2011-10-13 | 2013-04-18 | Charles Terrance Hatch | Methods and systems for automatic rolling-element bearing fault detection |
CN102620928A (en) * | 2012-03-02 | 2012-08-01 | 燕山大学 | Wind-power gear box fault diagnosis method based on wavelet medium-soft threshold and electronic-magnetic diaphragm (EMD) |
CN103837345B (en) * | 2014-03-25 | 2016-06-22 | 上海电机学院 | Fault Diagnosis of Gear Case method and device |
CN104019831B (en) * | 2014-06-20 | 2017-01-04 | 哈尔滨工业大学 | Gyroscope method for diagnosing faults based on EMD and entropy weight |
CN104729853B (en) * | 2015-04-10 | 2017-06-06 | 华东交通大学 | A kind of rolling bearing performance degradation assessment device and method |
-
2016
- 2016-06-29 CN CN201610492429.6A patent/CN106096200B/en active Active
- 2016-06-29 CN CN201811110099.5A patent/CN109635306A/en not_active Withdrawn
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110929789A (en) * | 2019-11-22 | 2020-03-27 | 北京理工大学 | Liver tumor automatic classification method and device based on multi-stage CT image analysis |
CN111769810A (en) * | 2020-06-29 | 2020-10-13 | 浙江大学 | Fluid mechanical modulation frequency extraction method based on energy kurtosis spectrum |
CN112101174A (en) * | 2020-09-09 | 2020-12-18 | 洛阳师范学院 | LOF-Kurtogram-based mechanical fault diagnosis method |
Also Published As
Publication number | Publication date |
---|---|
CN106096200B (en) | 2019-02-22 |
CN106096200A (en) | 2016-11-09 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106096313B (en) | A kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis | |
CN106096200B (en) | A kind of envelope Analysis Method based on wavelet decomposition and spectrum kurtosis | |
CN106198015B (en) | A kind of VMD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
CN106153339B (en) | A kind of envelope Analysis Method based on the filtering of variation Mode Decomposition | |
CN106198013B (en) | A kind of envelope Analysis Method based on empirical mode decomposition filtering | |
CN106096198B (en) | A kind of envelope Analysis Method based on variation Mode Decomposition and spectrum kurtosis | |
CN106096199B (en) | A kind of WT of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
CN105954031B (en) | A kind of envelope Analysis Method based on unusual spectral factorization filtering | |
CN106198009B (en) | A kind of EMD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
CN106053069B (en) | A kind of SSD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
CN106153333B (en) | A kind of envelope Analysis Method based on wavelet decomposition filtering | |
CN106053060B (en) | A kind of envelope Analysis Method that filtering is decomposed based on nonlinear model | |
CN106198010B (en) | A kind of envelope Analysis Method that filtering is decomposed based on local mean value | |
CN106198012B (en) | A kind of envelope Analysis Method for decomposing and composing kurtosis based on local mean value | |
CN106096201B (en) | A kind of EEMD and smoothed cubic spline envelope Analysis Method of rotating machinery | |
CN106198014B (en) | A kind of envelope Analysis Method based on empirical mode decomposition and spectrum kurtosis | |
CN106053059B (en) | It is a kind of based on it is interior grasp time scale decompose filtering envelope Analysis Method | |
CN106053061B (en) | A kind of envelope Analysis Method for decomposing and composing kurtosis based on nonlinear model | |
CN105954030B (en) | It is a kind of based on it is interior grasp time scale decompose and spectrum kurtosis envelope Analysis Method | |
CN106198018B (en) | A kind of EEMD of rotating machinery and smooth iteration envelope Analysis Method | |
CN106198017B (en) | A kind of LMD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
CN105973603B (en) | The EEMD and rational spline smoothed envelope analysis method of a kind of rotating machinery | |
CN106124200B (en) | A kind of ELMD of rotating machinery and smooth iteration envelope Analysis Method | |
CN106198016B (en) | A kind of NMD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
CN106153338B (en) | A kind of ELMD of rotating machinery and rational spline smoothed envelope analysis method |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
WW01 | Invention patent application withdrawn after publication |
Application publication date: 20190416 |
|
WW01 | Invention patent application withdrawn after publication |