CN106198013B - A kind of envelope Analysis Method based on empirical mode decomposition filtering - Google Patents
A kind of envelope Analysis Method based on empirical mode decomposition filtering Download PDFInfo
- Publication number
- CN106198013B CN106198013B CN201610492435.1A CN201610492435A CN106198013B CN 106198013 B CN106198013 B CN 106198013B CN 201610492435 A CN201610492435 A CN 201610492435A CN 106198013 B CN106198013 B CN 106198013B
- Authority
- CN
- China
- Prior art keywords
- signal
- envelope
- data
- empirical mode
- mode decomposition
- 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.)
- Active
Links
Classifications
-
- 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
- G01M13/045—Acoustic or vibration analysis
Landscapes
- Physics & Mathematics (AREA)
- Acoustics & Sound (AREA)
- General Physics & Mathematics (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
Abstract
The invention discloses a kind of envelope Analysis Methods based on empirical mode decomposition filtering, this method decomposes original signal first with ensemble empirical mode decomposition method, then the noise component(s) and trend term in decomposition result are excluded using the rearrangement and replacement operation of data, 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 wave filter, then Envelope Analysis is carried out to second of filtered signal using rational 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 determine the fault type of rotating machinery exactly, have good noise immunity and robustness, convenient for engineer application.
Description
Technical field
The present invention relates to condition monitoring for rotating machinery and fault diagnosis field, and in particular to one kind is based on empirical mode decomposition
The envelope Analysis Method of filtering.
Background technology
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 directly analyzes original signal or only to original
Signal analyzed again after simply filtering, therefore existing method is easily done be subject to noise, trend and other ingredients
It disturbs, it is relatively low so as to cause the analysis precision of the prior art；2. existing Envelope Analysis technology is based on Hilbert is converted,
And the signal that Hilbert conversion requirements are analyzed must be the narrow band signal of simple component, otherwise the frequency modulating section of signal will
The amplitude envelope analysis result of signal is polluted, but signal to be analyzed at present does not meet the item of simple component and narrowband strictly
Part so may result in the prior art and erroneous judgement problem be susceptible to due to precision is not high；3. the envelope spectrum obtained by conventional method
There are end effects.
The content of the invention
The problem to be solved in the present invention is for above deficiency, proposes a kind of envelope point based on empirical mode decomposition filtering
Analysis method after envelope Analysis Method using the present invention, has that analysis result accuracy and accuracy are high, and can detect exactly
The advantages of going out rotating machinery fault type.
For solution more than technical problem, the technical solution that the present invention takes is as follows：One kind is filtered based on empirical mode decomposition
Envelope Analysis Method, which is characterized in that comprise the following steps：
Step 1：The vibration signal x of rotating machinery is measured with sample frequency fs using acceleration transducer（k）, （k=1,
2, …,N）, N is the length of sampled signal；
Step 2：Using empirical mode decomposition（Empirical Mode Decomposition, EMD）Algorithm is by signal x
（k）The sum of n component and trend term are resolved into, i.e.,, wherein, ci（k）It represents
I-th of the component obtained by EMD algorithms, rn（k）Represent the trend term obtained by EMD algorithms；
Step 3：To ci（k）It performs reordering operations and substitutes and operate, it is rearranged to operate obtained data ci shuffle（k）Table
Show, data c is obtained after substituting operationi FTran（k）It represents；
Step 4：To ci（k）、ci shuffle（k）And ci FTran（k）Multi-fractal is performed respectively removes trend fluction analysis
（Multifractal Detrended Fluctuation Analysis, MFDFA）, obtain generalized Hurst index curve, ci
（k）Generalized Hurst index curve Hi（q）It represents；
Step 5：If Hi（q）With Hi shuffle（q）Or Hi（q）With Hi FTran（q）Between relative error be less than 5% or
Person 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 performed, the centre frequency f corresponding to signal kurtosis maximum is obtained0And band
Wide B；
Step 8：According to centre frequency f0With bandwidth B to xf1（k）Bandpass filtering is carried out, obtains xf2（k）；
Step 9：To signal xf2（k）The analysis of rational spline iteration smoothed envelope is performed, obtains signal envelope eov（k）；
Step 10：To obtained signal envelope eov（k）It performs discrete Fourier transform and obtains envelope spectrum, according to envelope spectrum
Characteristic frequency judges the fault type of machine.
A kind of prioritization scheme, empirical mode decomposition in the step 2（EMD）Algorithm comprises the following steps：
1）First screening process：Data x is found out respectively（k）Upper and lower Local Extremum, using cubic spline curve
Upper and lower Local Extremum is fitted respectively, obtains signal x（k）Local maximum envelope and local minimum envelope, then will
The value of the respective points of this two envelopes is averaged, and obtains an averaged curve m1；
Original signal x is sought again（k）With this averaged curve m1Difference, i.e. h10=x(k)-m1, so far first screening process knot
Beam；
2）Second screening process：h 10Again original signal is taken as, repeats the above steps 1）, can obtainh 11= h 10-m 11, parameter herem 11It representsh 10Mean curve, this process j times is repeated, until 0.2<SD<Screening process stops when 0.3,
Here, at this point,h 1j = h 1(j-1)-m 1j , at this moment can considerh 1j It is to grasp mould in one
State function（Intrinsic Mode Function, IMF）, defining the 1st IMF isc 1=h 1j ；
3）It is subtracted from original signalc 1, can obtainr 1=x(k)-c 1, then willr 1As new data, and repeat above-mentioned two steps behaviour
Make, can so obtain the 2nd IMF；
4）Repeat step 3）Operation can obtain a series of IMF, ifr n A monotonous curve is had changed into, then screening process
Stop, most original signal is decomposed into following form at last：。
Further, data rearrangement operation comprises the following steps in the step 3：
Upset component c at randomi（k）Put in order.
Further, the operation of data replacement comprises the following steps in the step 3：
1）To component ci（k）Discrete Fourier transform is performed, obtains component ci（k）Phase；
2）It is located at the pseudo- independent same distribution number in (- π, π) section with one group to replace component ci（k）Original phase；
3）Inverse discrete Fourier transform is performed to the frequency domain data after phase substitutes and obtains data ci IFFT（k）, ask for
Data ci IFFT（k）Real part.
Further, MFDFA methods comprise the following steps in the step 4：
1）Constructionx(k)(k=1,2,…,N) profileY(i)：
；
x(k) represent c in step 4 described in claim 1i（k）Or ci shuffle（k）Or ci FTran（k）；
2）By signal profileY(i) be divided into it is nonoverlappingN s Segment length issData, from the negative directions of data with identical
Length segmentation obtains 2N s Segment data；
3）Using polynomial trend of the least square fitting per segment data, the variance per segment data is then calculated：
；
y v (i) for the of fittingvThe trend of segment data, if the polynomial trend of fitting ismRank then remembers that this goes trend process
For（MF-）DFAm；
4）Calculate theqThe average value of rank wave function：
；
5）Ifx(k) there are self-similarity characteristics, thenqThe average value of rank wave functionF q (s) and time scalesBetween
There are power law relations：
F q (s)~s H(q)；
WhenqWhen=0, step 4）In formula diverging, at this momentH(0) by logarithmic mean process defined in following formula come really
It is fixed：
；
6）To step 5）In formula both sides take the logarithm can obtain ln [F q (s)]=H(q)ln(s)+c (cFor constant), thus
The slope of straight line can be obtainedH(q)。
Further, the spectrum kurtosis method in the step 7 comprises the following steps：
1）Constructing a cutoff frequency isf c =0.125+εLow-pass filterh(n)；ε>0；
2）It is based onh(n) passband is constructed as [0;0.25] quasi- low-pass filterh 0(n) and passband be [0.25;0.5]
Quasi- high-pass filterh 1 (n),
；
3）SignalThroughh 0(n)、h 1 (n) filter and resolve into low frequency part after down-sampledAnd radio-frequency head
Point, the down-sampled factor is 2, then the shaping filter tree after successive ignition filters, kth layer have 2 k A frequency band,
InIt represents the in wave filter treekOn layeriThe output signal of a wave filter,i=0,…,2k- 1,0≤k≤K-1；
4）In decomposition treekOn layeriThe centre frequency of a wave filterf ki And bandwidthB k Respectively
f ki =(i+2-1)2-k-1
B k =2-k-1；
5）Calculate each filter results ( i=0,…,2k- 1) kurtosis；
6）All spectrum kurtosis are summarized, obtain the total spectrum kurtosis of signal.
Further, the rational spline iteration smoothed envelope analysis method in the step 9 comprises the following steps：
1）Calculate signalz(k) Jue Dui Zhi ∣z(k) the local extremum of ∣；In the 1st iteration,z(k) represent claim
X in 1 step 9f2（k）；
2）Envelope eov is obtained using rational spline curve fitting Local Extremum1（k）；
3）It is rightz(k) be normalized to obtain；
4）2nd iteration：z 1(k) new data is re-used as, repeat above-mentioned steps 1）~3）, obtain；
5）Ith iteration：z i-1(k) new data is re-used as, repeat above-mentioned steps 1）~ 3), obtain；
If 6)nWhat secondary iteration obtainedz n (k) amplitude be less than or equal to 1, then iterative process stop, finally obtaining
Signalz(k) envelope be。
The present invention is using above technical scheme, and compared with prior art, the present invention has the following advantages：
1) original signal is decomposed using EMD, then excludes therein make an uproar using the rearrangement and replacement operation of data
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 rational spline iteration smoothed envelope analysis method, energy
Influence of the frequency modulating section to signal envelope analysis result is enough avoided, so as to improve the precision of Envelope Analysis.
3) fault type of rotating machinery can be detected exactly.
4) envelope spectrum obtained can avoid end effect.
The present invention will be further described with reference to the accompanying drawings and examples.
Description of the drawings
Attached drawing 1 is the flow chart of envelope Analysis Method in the embodiment of the present invention；
Attached drawing 2 is the schematic diagram for carrying out preliminary exposition in the embodiment of the present invention to signal using low pass and high-pass filter；
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 for having in the embodiment of the present invention inner ring failure；
Attached drawing 5 is to inner ring faulty bearing vibration signal in the embodiment of the present invention using traditional envelope Analysis Method
Analysis result；
Attached drawing 6 is to inner ring faulty bearing vibration signal in the embodiment of the present invention using envelope Analysis Method of the present invention
Analysis result；
Attached drawing 7 is the bearing vibration signal for having in the embodiment of the present invention outer ring failure；
Attached drawing 8 is to outer ring faulty bearing vibration signal in the embodiment of the present invention using traditional envelope Analysis Method
Analysis result；
Attached drawing 9 is to outer ring faulty bearing vibration signal in the embodiment of the present invention using envelope Analysis Method of the present invention
Analysis result.
Specific embodiment
Embodiment, as shown in Figure 1, Figure 2, Figure 3 shows, a kind of envelope Analysis Method based on empirical mode decomposition filtering, including
Following steps：
Step 1：The vibration signal x of rotating machinery is measured with sample frequency fs using acceleration transducer（k）, （k=1,
2, …,N）, N is the length of sampled signal；
Step 2：Using empirical mode decomposition（Empirical Mode Decomposition, EMD）Algorithm is by signal x
（k）The sum of n component and trend term are resolved into, i.e.,, wherein, ci（k）Generation
I-th of component that table is obtained by EMD algorithms, rn（k）Represent the trend term obtained by EMD algorithms；
Step 3：To ci（k）It performs reordering operations and substitutes and operate, it is rearranged to operate obtained data ci shuffle（k）Table
Show, data c is obtained after substituting operationi FTran（k）It represents；
Step 4：To ci（k）、ci shuffle（k）And ci FTran（k）Multi-fractal is performed respectively removes trend fluction analysis
（Multifractal Detrended Fluctuation Analysis, MFDFA）, obtain generalized Hurst index curve, ci
（k）Generalized Hurst index curve Hi（q）It represents；
Step 5：If Hi（q）With Hi shuffle（q）Or Hi（q）With Hi FTran（q）Between relative error be less than 5% or
Person 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 performed, the centre frequency f corresponding to signal kurtosis maximum is obtained0And band
Wide B；
Step 8：According to centre frequency f0With bandwidth B to xf1（k）Bandpass filtering is carried out, obtains xf2（k）；
Step 9：To signal xf2（k）The analysis of rational spline iteration smoothed envelope is performed, obtains signal envelope eov（k）；
Step 10：To obtained signal envelope eov（k）It performs discrete Fourier transform and obtains envelope spectrum, according to envelope spectrum
Characteristic frequency judges the fault type of machine.
Empirical mode decomposition in step 2（EMD）Algorithm comprises the following steps：
1）First screening process：Data x is found out respectively（k）Upper and lower Local Extremum, using cubic spline curve
Upper and lower Local Extremum is fitted respectively, obtains signal x（k）Local maximum envelope and local minimum envelope, then will
The value of the respective points of this two envelopes is averaged, and obtains an averaged curve m1；
Original signal x is sought again（k）With this averaged curve m1Difference, i.e. h10=x(k)-m1, so far first screening process knot
Beam；
2）Second screening process：h 10Again original signal is taken as, repeats the above steps 1）, can obtainh 11= h 10-m 11, parameter herem 11It representsh 10Mean curve, this process j times is repeated, until 0.2<SD<Screening process stops when 0.3,
Here, at this point,h 1j = h 1(j-1)-m 1j , at this moment can considerh 1j It is to grasp mould in one
State function（Intrinsic Mode Function, IMF）, defining the 1st IMF isc 1=h 1j ；
3）It is subtracted from original signalc 1, can obtainr 1=x(k)-c 1, then willr 1As new data, and repeat above-mentioned two steps behaviour
Make, can so obtain the 2nd IMF；
4）Repeat step 3）Operation can obtain a series of IMF, ifr n A monotonous curve is had changed into, then screening process
Stop, most original signal is decomposed into following form at last：。
Data rearrangement operation comprises the following steps in step 3：
Upset component c at randomi（k）Put in order.
Data substitute operation and comprise the following steps in step 3：
1）To component ci（k）Discrete Fourier transform is performed, obtains component ci（k）Phase；
2）It is located at the pseudo- independent same distribution number in (- π, π) section with one group to replace component ci（k）Original phase；
3）Inverse discrete Fourier transform is performed to the frequency domain data after phase substitutes and obtains data ci IFFT（k）, ask for
Data ci IFFT（k）Real part.
MFDFA methods comprise the following steps in step 4：
1）Constructionx(k)(k=1,2,…,N) profileY(i)：
；
x(k) represent c in step 4 described in claim 1i（k）Or ci shuffle（k）Or ci FTran（k）；
2）By signal profileY(i) be divided into it is nonoverlappingN s Segment length issData, from the negative directions of data with identical
Length segmentation obtains 2N s Segment data；
3）Using polynomial trend of the least square fitting per segment data, the variance per segment data is then calculated：
；
y v (i) for the of fittingvThe trend of segment data, if the polynomial trend of fitting ismRank then remembers that this goes trend process
For（MF-）DFAm；
4）Calculate theqThe average value of rank wave function：
；
5）Ifx(k) there are self-similarity characteristics, thenqThe average value of rank wave functionF q (s) and time scalesBetween
There are power law relations：
F q (s)~s H(q)；
WhenqWhen=0, step 4）In formula diverging, at this momentH(0) by logarithmic mean process defined in following formula come really
It is fixed：
；
6）To step 5）In formula both sides take the logarithm can obtain ln [F q (s)]=H(q)ln(s)+c (cFor constant), thus
The slope of straight line can be obtainedH(q)。
Spectrum kurtosis method in step 7 comprises the following steps：
1）Constructing a cutoff frequency isf c =0.125+εLow-pass filterh(n)；ε>0；
2）It is based onh(n) passband is constructed as [0;0.25] quasi- low-pass filterh 0(n) and passband be [0.25;0.5]
Quasi- high-pass filterh 1 (n),
；
3）SignalThroughh 0(n)、h 1 (n) filter and resolve into low frequency part after down-sampledAnd high frequency
Part, the down-sampled factor is 2, then the shaping filter tree after successive ignition filters, kth layer have 2 k A frequency band,
WhereinIt represents the in wave filter treekOn layeriThe output signal of a wave filter,i=0,…,2k- 1,0≤k≤K-
1；
4）In decomposition treekOn layeriThe centre frequency of a wave filterf ki And bandwidthB k Respectively
f ki =(i+2-1)2-k-1
B k =2-k-1；
5）Calculate each filter results ( i=0,…,2k- 1) kurtosis；
6）All spectrum kurtosis are summarized, obtain the total spectrum kurtosis of signal.
Rational spline iteration smoothed envelope analysis method in step 9 comprises the following steps：
1）Calculate signalz(k) Jue Dui Zhi ∣z(k) the local extremum of ∣；In the 1st iteration,z(k) represent claim
X in 1 step 9f2（k）；
2）Envelope eov is obtained using rational spline curve fitting Local Extremum1（k）；
3）It is rightz(k) be normalized to obtain；
4）2nd iteration：z 1(k) new data is re-used as, repeat above-mentioned steps 1）~3）, obtain；
5）Ith iteration：z i-1(k) new data is re-used as, repeat above-mentioned steps 1）~ 3), obtain；
If 6)nWhat secondary iteration obtainedz n (k) amplitude be less than or equal to 1, then iterative process stop, finally obtaining
Signalz(k) envelope be。
Experiment 1, tests the performance of algorithm of the present invention using the bearing vibration data with inner ring failure
Card.
Bearing used in experiment be 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, during signal sampling
A length of 1s.
The inner ring fault-signal collected 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, from fig. 5, it can be seen that the left end point of envelope spectrum has abnormal high level, this explanation is by tradition
There are end effects for the envelope spectrum that method obtains.
Signal shown in Fig. 4 is analyzed using method proposed by the invention, obtained analysis result such as Fig. 6 institutes
Show.From fig. 6, it can be seen that the spectral line corresponding to 160Hz and 320Hz, apparently higher than other spectral lines, the two frequencies correspond to respectively
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
Go out, the envelope spectrum obtained by the present invention does not have end effect.
Show that the present invention is capable of reliable recognition in the case where loading and failure dimensional depth being constant through many experiments
Minimum inner ring failure dimension width is about 0.20mm, and conventional method is capable of the minimum inner ring failure dimension width of reliable recognition about
For 0.53mm, precision improves 62.3%.
Experiment 2, tests the performance of algorithm of the present invention using the bearing vibration data with outer ring failure
Card.
Bearing used in experiment be 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, during signal sampling
A length of 1s.
The outer ring fault-signal collected 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, from figure 8, it is seen that the left end point of envelope spectrum has abnormal high level, this explanation is by tradition
There are end effects for the envelope spectrum that method obtains.
Signal shown in Fig. 7 is analyzed using method proposed by the invention, obtained analysis result such as Fig. 9 institutes
Show.From fig. 9, it can be seen that the spectral line corresponding to 103Hz and 206Hz, apparently higher than other spectral lines, the two frequencies correspond to respectively
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
Go out, the envelope spectrum obtained by the present invention does not have end effect.
Show that the present invention is capable of reliable recognition in the case where loading and failure dimensional depth being constant through many experiments
Minimum outer ring failure dimension width is about 0.30mm, and conventional method is capable of the minimum outer ring failure dimension width of reliable recognition about
For 0.68mm, precision improves 55.9%.
According to result of the test, think after analysis：
1) traditional envelope Analysis Method directly carries out original signal Envelope Analysis or to merely through simple process
Original signal afterwards carries out Envelope Analysis, different from traditional envelope Analysis Method, and the present invention is first with EMD to original signal
It is decomposed, then excludes noise and trend component therein using the rearrangement and replacement operation of data, only stick signal point
Useful component in amount so as to avoid the influence of noise and trend component to Envelope Analysis result, improves accuracy and essence
Exactness.
2) traditional envelope Analysis Method is based on Hilbert is converted, and the letter that Hilbert conversion requirements are analyzed
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, so may result in the prior art because of precision not
High and be susceptible to erroneous judgement problem, different from traditional envelope Analysis Method, the present invention is divided using rational spline iteration smoothed envelope
Signal envelope and frequency modulating section are kept completely separate by analysis method, can avoid frequency modulating section to signal envelope analysis result
Influence, so as to improve the precision of Envelope Analysis.
3) fault type of rotating machinery can be detected exactly.
4) envelope spectrum obtained by conventional method is there are end effect, and can be avoided by the envelope spectrum that the present invention obtains
End effect.
5）Each step effect：
1) step：Gather vibration signal；
2) step：Original signal is resolved into the form of different component sums, some of which component corresponds to noise and trend term,
Some components correspond to useful signal；
3) ~ 5) step：Is performed by reordering operations and is substituted for the signal that above-mentioned decomposition obtains and is operated, rejects noise therein point
Amount and trend term only retain useful signal；
6) step：Remaining useful signal is summed, using this and it is rearranged as signal and substitute filtered result xf1
(k)；
7) step：To filtered signal xf1(k) spectrum kurtosis analysis is performed, corresponding center at signal maximum kurtosis is obtained
Frequency f0And bandwidth B；
8) step：According to centre frequency f0With bandwidth B to xf1(k) bandpass filtering is carried out, obtains signal xf2(k)；
9) step：Calculate signal xf2(k) envelope eov (k)；
10) step：Discrete Fourier transform is performed to eov (k) and obtains envelope spectrum, the failure of bearing is judged according to envelope spectrum
Type.
One skilled in the art would recognize that above-mentioned specific embodiment is only exemplary, it is to make ability
Field technique personnel can be better understood from present invention, should not be understood as limiting the scope of the invention, as long as
Technical solution improvements introduced according to the present invention each falls within protection scope of the present invention.
Claims (7)
1. a kind of envelope Analysis Method based on empirical mode decomposition filtering, which is characterized in that comprise the following steps：
Step 1：The vibration signal x of rotating machinery is measured with sample frequency fs using acceleration transducer（k）, k=1, 2, …,
N, N are the length of sampled signal；
Step 2：Using empirical mode decomposition（Empirical Mode Decomposition, EMD）Algorithm is by vibration signal x
（k）The sum of n component and trend term are resolved into, i.e.,, wherein, ci（k）Represent by
I-th of component that EMD algorithms obtain, rn（k）Represent the trend term obtained by EMD algorithms；
Step 3：To ci（k）It performs reordering operations and substitutes and operate, it is rearranged to operate obtained data ci shuffle（k）It represents,
Data c is obtained after substituting operationi FTran（k）It represents；
Step 4：To ci（k）、ci shuffle（k）And ci FTran（k）Multi-fractal is performed respectively removes trend fluction analysis
（Multifractal Detrended Fluctuation Analysis, MFDFA）, obtain generalized Hurst index curve, ci
（k）Generalized Hurst index curve Hi（q）It represents；
Step 5：If Hi（q）With Hi shuffle（q）Or Hi（q）With Hi FTran（q）Between relative error be 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 performed, the centre frequency f corresponding to signal kurtosis maximum is obtained0And bandwidth B；
Step 8：According to centre frequency f0With bandwidth B to xf1（k）Bandpass filtering is carried out, obtains xf2（k）；
Step 9：To signal xf2（k）The analysis of rational spline iteration smoothed envelope is performed, obtains signal envelope eov（k）；
Step 10：To obtained signal envelope eov（k）It performs discrete Fourier transform and obtains envelope spectrum, according to envelope spectrum signature
Frequency judges the fault type of machine.
A kind of 2. envelope Analysis Method based on empirical mode decomposition filtering according to claim 1, which is characterized in that institute
State empirical mode decomposition in step 2（EMD）Algorithm comprises the following steps：
1）First screening process：Vibration signal x is found out respectively（k）Upper and lower Local Extremum, using cubic spline curve
Upper and lower Local Extremum is fitted respectively, obtains vibration signal x（k）Local maximum envelope and local minimum envelope,
The value of the respective points of this two envelopes is averaged again, obtains an averaged curve m1；
Vibration signal x is sought again（k）With this averaged curve m1Difference, i.e. h10=x(k)-m1, so far first screening process terminate；
2）Second screening process：h 10Again vibration signal is taken as, repeats the above steps 1）, can obtainh 11= h 10-m 11, this
In parameterm 11It representsh 10Mean curve, this process j times is repeated, until 0.2<SD<Screening process stops when 0.3, here, at this point,h 1j = h 1(j-1)-m 1j , at this moment thinkh 1j It is to grasp mode function in one
（Intrinsic Mode Function, IMF）, defining the 1st IMF isc 1=h 1j ；
3）It is subtracted from vibration signalc 1, can obtainr 1=x(k)-c 1, then willr 1As new data, and above-mentioned two steps operation is repeated,
So obtain the 2nd IMF；
4）Repeat step 3）Operation can obtain a series of IMF, ifr n Having changed into a monotonous curve, then screening process stops,
Most divided oscillation signal solution is following form at last：。
A kind of 3. envelope Analysis Method based on empirical mode decomposition filtering according to claim 1, which is characterized in that institute
Data rearrangement operation in step 3 is stated to comprise the following steps：
Upset component c at randomi（k）Put in order.
4. a kind of envelope Analysis Method based on empirical mode decomposition filtering according to claim 1, it is characterised in that：Institute
It states data in step 3 and substitutes to operate and comprise the following steps：
1）To component ci（k）Discrete Fourier transform is performed, obtains component ci（k）Phase；
2）It is located at the pseudo- independent same distribution number in (- π, π) section with one group to replace component ci（k）Original phase；
3）Inverse discrete Fourier transform is performed to the frequency domain data after phase substitutes and obtains data ci IFFT（k）, ask for data
ci IFFT（k）Real part.
5. a kind of envelope Analysis Method based on empirical mode decomposition filtering according to claim 1, it is characterised in that：Institute
MFDFA methods in step 4 are stated to comprise the following steps：
1）Construct xy(k) profileY(i), k=1,2 ..., N：
,
,
xy(k) c in step 4 described in claim 1 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 negative direction of data 2 are obtained with identical length segmentation, such oneN s Section
Data；
3）Using polynomial trend of the least square fitting per segment data, the variance per segment data is then calculated：
4P`F)18O2626ZX3J9SV283R
Y$[`H}@AV[I7X%IJ6A_JQWC；
y v (i) for 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；In this example, m=1；
4）Calculate theqThe average value of rank wave function：
0R(8CJGG2TXCM$2NT2@((TH；
5）If xy(k) there are self-similarity characteristics, thenqThe average value of rank wave functionF q (s) and time scalesBetween exist
Power law relation：
F q (s)~s H(q)；
WhenqWhen=0, step 4）In formula diverging, at this momentH(0) determined by logarithmic mean process defined in following formula：
S42M8HF%L8]B11YUNN21RY7；
6）To step 5）In formula both sides take the logarithm can obtain ln [F q (s)]=H(q)ln(s)+c,cFor constant, straight line is derived from
SlopeH(q)。
6. a kind of envelope Analysis Method based on empirical mode decomposition filtering according to claim 1, it is characterised in that：Institute
The spectrum kurtosis method stated in step 7 comprises the following steps：
1）Constructing a cutoff frequency isf c =0.125+εLow-pass filterh(n)；ε>0；
2）It is based onh(n) passband is constructed as [0;0.25] quasi- low-pass filterh 0(n) and passband be [0.25;0.5] standard
High-pass filterh 1 (n),
；
3）SignalThroughh 0(n)、h 1 (n) filter and resolve into low frequency part after down-sampledAnd high frequency section, the down-sampled factor is 2, then the shaping filter tree after successive ignition filters, kth layer have 2 k A frequency band, whereinIt represents the in wave filter treekOn layeriThe output signal of a wave filter,i=0,…,2k- 1,0≤k≤K-1；
4）In decomposition treekOn layeriThe centre frequency of a wave filterf ki And bandwidthB k Respectively
f ki =(i+2-1)2-k-1
B k =2-k-1；
5）Calculate each filter results i=0,…,2k- 1 kurtosis；
6）All spectrum kurtosis are summarized, obtain the total spectrum kurtosis of signal.
A kind of 7. envelope Analysis Method based on empirical mode decomposition filtering according to claim 1, which is characterized in that institute
The rational spline iteration smoothed envelope analysis method stated in step 9 comprises the following steps：
1）Calculate signalz(k) Jue Dui Zhi ∣z(k) the local extremum of ∣；In the 1st iteration,z(k) represent claim 1 institute
State x in step 9f2（k）；
2）Envelope eov is obtained using rational spline curve fitting Local Extremum1（k）；
3）It is rightz(k) be normalized to obtain；
4）2nd iteration：z 1(k) new data is re-used as, repeat above-mentioned steps 1）~3）, obtain；
5）Ith iteration：z i-1(k) new data is re-used as, repeat above-mentioned steps 1）~ 3), obtain；
If 6)nWhat secondary iteration obtainedz n (k) amplitude be less than or equal to 1, then iterative process stop, finally obtaining signalz
(k) envelope be。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610492435.1A CN106198013B (en) | 2016-06-29 | 2016-06-29 | A kind of envelope Analysis Method based on empirical mode decomposition filtering |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610492435.1A CN106198013B (en) | 2016-06-29 | 2016-06-29 | A kind of envelope Analysis Method based on empirical mode decomposition filtering |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106198013A CN106198013A (en) | 2016-12-07 |
| CN106198013B true CN106198013B (en) | 2018-05-29 |
Family
ID=57461056
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610492435.1A Active CN106198013B (en) | 2016-06-29 | 2016-06-29 | A kind of envelope Analysis Method based on empirical mode decomposition filtering |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106198013B (en) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10361895B1 (en) * | 2018-01-22 | 2019-07-23 | Huawei Technologies Co., Ltd. | Systems and methods for hybrid multi-layer signal decomposition |
| CN109522522B (en) * | 2018-12-14 | 2020-11-06 | 上海递缇智能系统有限公司 | Method for generating and applying envelope curve of press-fitting and screwing curve |
| CN112099100B (en) * | 2020-08-25 | 2021-08-27 | 吉林大学 | Magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering |
| CN112633121A (en) * | 2020-12-18 | 2021-04-09 | 国网浙江省电力有限公司武义县供电公司 | Radiation source identification method based on Hilbert transform and multi-fractal dimension characteristics |
| CN113435281A (en) * | 2021-06-18 | 2021-09-24 | 秦皇岛北方管业有限公司 | Ripple compensator fault diagnosis method based on improved HHT conversion |
| CN113567561A (en) * | 2021-07-20 | 2021-10-29 | 宝宇(武汉)激光技术有限公司 | Metal nondestructive test equipment |
| CN116750290B (en) * | 2023-08-14 | 2023-10-31 | 青岛义龙包装机械有限公司 | Intelligent monitoring method for running state of bag type packaging machine |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2581725A2 (en) * | 2011-10-13 | 2013-04-17 | General Electric Company | Methods and systems for automatic rolling-element bearing fault detection |
| CN105004523A (en) * | 2015-08-04 | 2015-10-28 | 潍坊学院 | Rolling bearing state monitoring method based on weighted similarity measure |
| CN105067262A (en) * | 2015-08-04 | 2015-11-18 | 潍坊学院 | Rolling bearing state monitoring method |
| CN105092239A (en) * | 2014-05-09 | 2015-11-25 | 潍坊学院 | Method for detecting early stage fault of gear |
| CN105588717A (en) * | 2015-12-10 | 2016-05-18 | 潍坊学院 | Gearbox fault diagnosis method |
-
2016
- 2016-06-29 CN CN201610492435.1A patent/CN106198013B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2581725A2 (en) * | 2011-10-13 | 2013-04-17 | General Electric Company | Methods and systems for automatic rolling-element bearing fault detection |
| CN105092239A (en) * | 2014-05-09 | 2015-11-25 | 潍坊学院 | Method for detecting early stage fault of gear |
| CN105004523A (en) * | 2015-08-04 | 2015-10-28 | 潍坊学院 | Rolling bearing state monitoring method based on weighted similarity measure |
| CN105067262A (en) * | 2015-08-04 | 2015-11-18 | 潍坊学院 | Rolling bearing state monitoring method |
| CN105588717A (en) * | 2015-12-10 | 2016-05-18 | 潍坊学院 | Gearbox fault diagnosis method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106198013A (en) | 2016-12-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106198013B (en) | A kind of envelope Analysis Method based on empirical mode decomposition filtering | |
| CN106198015B (en) | A kind of VMD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN106096313B (en) | A kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis | |
| CN106153339B (en) | A kind of envelope Analysis Method based on the filtering of variation Mode Decomposition | |
| CN106096200B (en) | A kind of envelope Analysis Method based on wavelet decomposition and spectrum kurtosis | |
| CN106168538B (en) | A kind of ITD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN106096198B (en) | A kind of envelope Analysis Method based on variation Mode Decomposition and spectrum kurtosis | |
| CN106198009B (en) | A kind of EMD 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 | |
| CN106198012B (en) | A kind of envelope Analysis Method for decomposing and composing kurtosis based on local mean value | |
| CN106096199B (en) | A kind of WT of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN106053060B (en) | A kind of envelope Analysis Method that filtering is decomposed based on nonlinear model | |
| CN106153333B (en) | A kind of envelope Analysis Method based on wavelet decomposition filtering | |
| CN106198010B (en) | A kind of envelope Analysis Method that filtering is decomposed based on local mean value | |
| CN106198014B (en) | A kind of envelope Analysis Method based on empirical mode decomposition and spectrum kurtosis | |
| CN106198018B (en) | A kind of EEMD of rotating machinery and smooth iteration envelope Analysis Method | |
| CN106096201B (en) | A kind of EEMD and smoothed cubic spline envelope Analysis Method of rotating machinery | |
| CN105973603B (en) | The EEMD and rational spline smoothed envelope analysis method of a kind of rotating machinery | |
| CN106053061B (en) | A kind of envelope Analysis Method for decomposing and composing kurtosis based on nonlinear model | |
| CN106198017B (en) | A kind of LMD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN106053059B (en) | It is a kind of based on it is interior grasp time scale decompose filtering envelope Analysis Method | |
| CN105954030B (en) | It is a kind of based on it is interior grasp time scale decompose and spectrum kurtosis envelope Analysis Method | |
| CN106198016B (en) | A kind of NMD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN106124200B (en) | A kind of ELMD of rotating machinery 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 |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |