CN105973603A - EEMD and rational spline smooth envelope analysis method for rotating machine - Google Patents
EEMD and rational spline smooth envelope analysis method for rotating machine Download PDFInfo
- Publication number
- CN105973603A CN105973603A CN201610492064.7A CN201610492064A CN105973603A CN 105973603 A CN105973603 A CN 105973603A CN 201610492064 A CN201610492064 A CN 201610492064A CN 105973603 A CN105973603 A CN 105973603A
- Authority
- CN
- China
- Prior art keywords
- signal
- envelope
- data
- eemd
- component
- 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.)
- Granted
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)
Abstract
The invention discloses an EEMD and rational spline smooth envelope analysis method for a rotating machine, and the method comprises the steps: firstly carrying out the decomposition of an original signal through employing an ensemble empirical mode decomposition method; secondly eliminating a noise component and a trend term in a decomposition result through employing the rearrangement and replacement operations of data; thirdly carrying out the analysis of the signal after first filtering through employing a spectral kurtosis method, and obtaining the central frequency and bandwidth of an optimal filter; fourthly carrying out the second filtering of the signal after first filtering through employing the filter; fifthly carrying out the envelope analysis of the signal after second filtering through employing a rational spline iteration smooth envelope analysis method; and finally determining the fault type of the rotating machine according to an envelope spectrum. The method is suitable for the processing of a complex fault signal of the rotating machine, can accurately judge the fault type of the rotating machine, is good in anti-noise performance and robustness, and is convenient for engineering application.
Description
Technical field
The present invention relates to condition monitoring for rotating machinery and fault diagnosis field, be specifically related to a kind of rotating machinery EEMD and
Rational spline smoothed envelope analyzes method.
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: the most existing Envelope Analysis technology or directly primary signal is analyzed, or only to original
Signal is analyzed after simply filtering again, and the most existing method is easily subject to the dry of noise, trend and other composition
Disturb, thus cause the analysis precision of prior art relatively low；The most existing Envelope Analysis technology is to be transformed to basis with Hilbert,
And Hilbert conversion requires that analyzed signal must be the narrow band signal of simple component, otherwise the frequency modulating section of signal will
The amplitude envelope analysis result of signal to be polluted, but signal the most to be analyzed the most strictly meets the bar of simple component and arrowband
Part, so may result in prior art and easily occurs erroneous judgement problem because precision is the highest；3. the envelope spectrum obtained by traditional method
There is end effect.
Summary of the invention
The problem to be solved in the present invention is for above not enough, and the EEMD and the rational spline that propose a kind of rotating machinery smooth
Envelope Analysis Method, after using the envelope Analysis Method of the present invention, has analysis result accuracy and degree of accuracy is high, and can be accurately
Detect rotating machinery fault type advantage.
For solving above technical problem, the technical scheme that the present invention takes is as follows: the EEMD of a kind of rotating machinery and reasonable
Spline smooth envelope Analysis Method, it is characterised in that comprise the following steps:
Step 1: utilize acceleration transducer to measure the vibration signal x(k of rotating machinery with sample frequency fs), (k=1,2,
..., N), N is the length of sampled signal；
Step 2: use set empirical mode decomposition (Ensemble Empirical Mode Decomposition, EEMD) to calculate
Method is by signal x(k) resolve into n component and a trend term sum, i.e., wherein, ciK () represents
The i-th component obtained by EEMD algorithm, rnK () represents the trend term obtained by EEMD algorithm；
Step 3: to ciK () performs reordering operations and substitutes operation, data c that rearranged operation obtainsi shuffleK () represents,
Data c are obtained after substituting operationi FTranK () represents；
Step 4: to ci(k), ci shuffle(k) and ci FTranK () performs multi-fractal respectively and removes trend fluction analysis
(Multifractal Detrended Fluctuation Analysis, MFDFA), obtains generalized Hurst index curve, ci
The generalized Hurst index curve H of (k)iQ () represents；ci shuffleThe generalized Hurst index curve H of (k)i shuffle(q) table
Show；ci FTranThe generalized Hurst index curve H of (k)i FTranQ () represents；
Step 5: if Hi(q) and Hi shuffle(q) or Hi(q) and Hi FTranQ the relative error between () is less than 5%, or Hi
(q), Hi shuffle(q) and Hi FTranQ () three do not change with q, then abandon the c of correspondencei(k) component；
Step 6: to remaining ci(k) component sue for peace, by this and be designated as signal rearranged and substitute filtered result xf1(k)；
Step 7: to xf1K () performs spectrum kurtosis analysis, obtain the mid frequency f corresponding to signal kurtosis maximum0And bandwidth B；
Step 8: according to mid frequency f0With bandwidth B to xf1K () carries out bandpass filtering, obtain xf2(k)；
Step 9: to signal xf2K () performs rational spline iteration smoothed envelope and analyzes, obtain signal envelope eov(k)；
Step 10: the signal envelope eov(k to obtaining) perform discrete Fourier transform obtain envelope spectrum, according to envelope spectrum feature
Frequency judges the fault type of machine.
A kind of prioritization scheme, gathers Empirical Mode Decomposition Algorithm and comprises the following steps in described step 2:
(1) to data x0K () is added white noise sequence and is produced a new data xj(k):
Std[x0(k)] represent data x0The standard deviation of (k), wnjK () represents wnjIn kth data, wnjRepresent jth with
The white noise sequence that machine produces, wnjAmplitude is 1,1≤j≤K；x0K () represents x (k) in step 2 described in claim 1；This example
In, K=100；
(2) to xjK () performs empirical mode decomposition, obtain n component and a trend term
cijK () represents xjK () performs the i-th component that empirical mode decomposition obtains, rnjK () represents xjK () performs experience
The trend term that Mode Decomposition obtains；
(3) meansigma methods of K decomposition result is calculated
ciK () represents x0K () carries out gathering the i-th component that empirical mode decomposition obtains, rnK () represents x0K () collects
Close the trend term that empirical mode decomposition obtains.
Further, in described step 3, data rearrangement operation comprises the following steps:
Upset component c at randomiPutting in order of (k).
Further, in described step 3, the operation of data replacement comprises the following steps:
1) to component ciK () performs discrete Fourier transform, it is thus achieved that component ciThe phase place of (k)；
2) component c is replaced with one group of pseudo-independent same distribution number being positioned in (-π, π) intervaliThe original phase of (k)；
3) frequency domain data after phase place substitutes is performed inverse discrete Fourier transform and obtain data ci IFFTK (), asks for number
According to ci IFFTThe real part of (k).
Further, in described step 4, MFDFA method comprises the following steps:
1) structure x (k) (k=1,2 ..., N) profile Y (i):
X (k) represents the c in step 4 described in claim 1i(k) or ci shuffle(k) or ci FTran(k)；
2) signal profile Y (i) is divided into nonoverlapping NSSegment length is the data of s, owing to data length N is generally not capable of dividing exactly s,
Can not utilize so one piece of data can be remained；
In order to make full use of the length of data, then from the opposite direction of data with identical length segmentation, obtain the most altogether 2NSSection
Data；
3) utilize the polynomial trend of the every segment data of least square fitting, then calculate the variance of every segment data:
yvI () is the trend of the v segment data of matching, if the polynomial trend of matching is m rank, then remember that this goes the trend process to be
(MF-) DFAm；In this example, m=1；
4) meansigma methods of q rank wave function is calculated:
；
5) if x (k) exists self-similarity characteristics, then meansigma methods F of q rank wave functionqExist between (s) and time scale s
Power law relation:
As q=0, the formula in step 4) dissipates, and at this moment H (0) is determined by logarithmic mean process defined in following formula:
6) are taken the logarithm in the formula both sides in step 5) and can obtain ln [Fq(s)]=H (q) ln (s)+c(c is constant), thus can obtain
Obtain slope H (q) of straight line.
Further, the spectrum kurtosis method in described step 7 comprises the following steps:
1) one cut-off frequency of structure is fcLow pass filter h (n) of=0.125+ ε；ε > 0, f in this examplec=0.3;
2) based on the quasi-low pass filter h that h (n) structure passband is [0,0.25]0N () and passband are [0.25,0.5]
Quasi-high pass filter h1(n),
;
3) signal ci kN () is through h0(n)、 h1N () resolves into low frequency part c after filtering and being down-sampled2i k+1(n) and HFS
c2i+1 k+1N (), the down-sampled factor is 2, then shaping filter tree after successive ignition filters, and kth layer has 2kIndividual frequency band, wherein
ci kThe output signal of the i-th wave filter on kth layer in (n) expression wave filter tree, i=0 ..., 2k-1,0≤k≤K-1, this example
Middle K=8；c0N () represents x in step 7 described in claim 1f1(k)；
4) the mid frequency f of the i-th wave filter on kth layer in decomposition treekiAnd bandwidth BkIt is respectively
;
5) each filter results c is calculatedi k(n)( i=0,…, 2k-1) kurtosis；
6) all of spectrum kurtosis is collected, obtain the spectrum kurtosis that signal is total.
Further, the analysis of the rational spline iteration smoothed envelope in described step 9 method comprises the following steps:
1) signal calculatedz(k) absolute valuez(k) local extremum；In the 1st iteration,z(k) represent claim 1 institute
State x in step 9f2(k)；
2) rational spline curve matching Local Extremum is used to obtain envelope eov1(k)；
3) rightz(k) be normalized and obtain；
4) the 2nd iteration:z 1(k) again as new data, repeat above-mentioned steps 1) ~ 3), obtain；
5) ith iteration:z i-1(k) again as new data, repeat above-mentioned steps 1) ~ 3), obtain；
6) ifnSecondary iteration obtainsz n (k) amplitude less than or equal to 1, then iterative process stops, and finally obtains signalz
(k) envelope be。
Further, in described step (2), Empirical Mode Decomposition Algorithm comprises the following steps:
1) first screening process: find out data x(k respectively) upper and lower Local Extremum, use cubic spline curve respectively
The upper and lower Local Extremum of matching, obtains signal x(k) local maximum envelope and local minimum envelope, then by this two
The value of the respective points of bar envelope is averaged, and obtains averaged curve m1；
Seek signal x(k again) and this averaged curve m1Difference, i.e. h10=x(k)-m1, so far first screening process terminates；
X(k) x in step described in claim 2 (2) is representedj(k)；
2) second screening process: h10Again new data, repeat the above steps 1 it are taken as), available h11= h10-m11, here
Parameter m11Represent h10Mean curve, repeat this process j time, < when 0.3, screening process stops, here until 0.2 < SD, now, h1j= h1(j-1)-m1j, at this moment it is believed that h1jMode function is grasped in being one
(Intrinsic Mode Function, IMF), defining the 1st IMF is c1= h1j；
3) from x(k) deduct c1, r can be obtained1=x(k)-c1, then by r1As new data, and repeat above-mentioned two step operations, this
Sample can obtain the 2nd IMF；
4) the available a series of IMF of step 3) operation are repeated, if rnHave changed into a monotonous curve, then screening process is stopped
Only, primary signal is decomposed into following form the most at last:。
The present invention uses above technical scheme, compared with prior art, the invention have the advantages that
1) utilize set empirical mode decomposition (EEMD) that primary signal is decomposed, then utilize the rearrangement of data and substitute behaviour
Make to get rid of noise therein and trend component, the only useful component in stick signal component, thus avoid noise and trend
The component impact on Envelope Analysis result, analysis result accuracy and degree of accuracy are high.
2) utilize rational spline iteration smoothed envelope to analyze method to be kept completely separate with frequency modulating section by signal envelope, energy
Enough avoid the frequency modulating section impact on signal envelope analysis result, thus improve the precision of Envelope Analysis.
3) fault type of rotating machinery can be detected exactly.
4) there is end effect in the envelope spectrum obtained by traditional method, and the envelope spectrum obtained by the present invention it can be avoided that
End effect.
The present invention will be further described with embodiment below in conjunction with the accompanying drawings.
Accompanying drawing explanation
Accompanying drawing 1 is the flow chart of envelope Analysis Method in the embodiment of the present invention；
Accompanying drawing 2 is to use low pass and high pass filter that signal is carried out the schematic diagram of preliminary exposition in the embodiment of the present invention；
Accompanying drawing 3 is the schematic diagram using tree-shaped filter construction quickly to calculate spectrum kurtosis in the embodiment of the present invention；
Accompanying drawing 4 is the bearing vibration signal in the embodiment of the present invention with inner ring fault；
Accompanying drawing 5 is to use the tradition envelope Analysis Method analysis to inner ring faulty bearing vibration signal in the embodiment of the present invention
Result；
Accompanying drawing 6 is to use envelope Analysis Method of the present invention to divide inner ring faulty bearing vibration signal in the embodiment of the present invention
Analysis result；
Accompanying drawing 7 is the bearing vibration signal in the embodiment of the present invention with outer ring fault；
Accompanying drawing 8 is to use the tradition envelope Analysis Method analysis to outer ring faulty bearing vibration signal in the embodiment of the present invention
Result；
Accompanying drawing 9 is to use envelope Analysis Method of the present invention to divide outer ring faulty bearing vibration signal in the embodiment of the present invention
Analysis result.
Detailed description of the invention
Embodiment, as shown in Figure 1, Figure 2, Figure 3 shows, the EEMD of a kind of rotating machinery and rational spline smoothed envelope analysis side
Method, comprises the following steps:
Step 1: utilize acceleration transducer to measure the vibration signal x(k of rotating machinery with sample frequency fs), (k=1,2,
..., N), N is the length of sampled signal；
Step 2: use set empirical mode decomposition (Ensemble Empirical Mode Decomposition, EEMD) to calculate
Method is by signal x(k) resolve into n component and a trend term sum, i.e., wherein, ci(k) generation
The i-th component that table is obtained by EEMD algorithm, rnK () represents the trend term obtained by EEMD algorithm；
Step 3: to ciK () performs reordering operations and substitutes operation, data c that rearranged operation obtainsi shuffleK () represents,
Data c are obtained after substituting operationi FTranK () represents；
Step 4: to ci(k), ci shuffle(k) and ci FTranK () performs multi-fractal respectively and removes trend fluction analysis
(Multifractal Detrended Fluctuation Analysis, MFDFA), obtains generalized Hurst index curve, ci
The generalized Hurst index curve H of (k)iQ () represents；ci shuffleThe generalized Hurst index curve H of (k)i shuffle(q) table
Show；ci FTranThe generalized Hurst index curve H of (k)i FTranQ () represents；
Step 5: if Hi(q) and Hi shuffle(q) or Hi(q) and Hi FTranQ the relative error between () is less than 5%, or Hi
(q), Hi shuffle(q) and Hi FTranQ () three do not change with q, then abandon the c of correspondencei(k) component；
Step 6: to remaining ci(k) component sue for peace, by this and be designated as signal rearranged and substitute filtered result xf1(k)；
Step 7: to xf1K () performs spectrum kurtosis analysis, obtain the mid frequency f corresponding to signal kurtosis maximum0And bandwidth B；
Step 8: according to mid frequency f0With bandwidth B to xf1K () carries out bandpass filtering, obtain xf2(k)；
Step 9: to signal xf2K () performs rational spline iteration smoothed envelope and analyzes, obtain signal envelope eov(k)；
Step 10: the signal envelope eov(k to obtaining) perform discrete Fourier transform obtain envelope spectrum, according to envelope spectrum feature
Frequency judges the fault type of machine.
Step 2 is gathered Empirical Mode Decomposition Algorithm comprise the following steps:
(1) to data x0K () is added white noise sequence and is produced a new data xj(k):
Std[x0(k)] represent data x0The standard deviation of (k), wnjK () represents wnjIn kth data, wnjRepresent jth with
The white noise sequence that machine produces, wnjAmplitude is 1,1≤j≤K；x0K () represents x (k) in step 2 described in claim 1；This example
In, K=100；
(2) to xjK () performs empirical mode decomposition, obtain n component and a trend term
cijK () represents xjK () performs the i-th component that empirical mode decomposition obtains, rnjK () represents xjK () performs experience
The trend term that Mode Decomposition obtains；
(3) meansigma methods of K decomposition result is calculated
ciK () represents x0K () carries out gathering the i-th component that empirical mode decomposition obtains, rnK () represents x0K () collects
Close the trend term that empirical mode decomposition obtains.
In step 3, data rearrangement operation comprises the following steps:
Upset component c at randomiPutting in order of (k).
In step 3, the operation of data replacement comprises the following steps:
1) to component ciK () performs discrete Fourier transform, it is thus achieved that component ciThe phase place of (k)；
2) component c is replaced with one group of pseudo-independent same distribution number being positioned in (-π, π) intervaliThe original phase of (k)；
3) frequency domain data after phase place substitutes is performed inverse discrete Fourier transform and obtain data ci IFFTK (), asks for number
According to ci IFFTThe real part of (k).
In step 4, MFDFA method comprises the following steps:
1) structure x (k) (k=1,2 ..., N) profile Y (i):
X (k) represents the c in step 4 described in claim 1i(k) or ci shuffle(k) or ci FTran(k)；
2) signal profile Y (i) is divided into nonoverlapping NSSegment length is the data of s, owing to data length N is generally not capable of dividing exactly s,
Can not utilize so one piece of data can be remained；
In order to make full use of the length of data, then from the opposite direction of data with identical length segmentation, obtain the most altogether 2NSSection
Data；
3) utilize the polynomial trend of the every segment data of least square fitting, then calculate the variance of every segment data:
yvI () is the trend of the v segment data of matching, if the polynomial trend of matching is m rank, then remember that this goes the trend process to be
(MF-) DFAm；In this example, m=1；
4) meansigma methods of q rank wave function is calculated:
；
5) if x (k) exists self-similarity characteristics, then meansigma methods F of q rank wave functionqExist between (s) and time scale s
Power law relation:
As q=0, the formula in step 4) dissipates, and at this moment H (0) is determined by logarithmic mean process defined in following formula:
6) are taken the logarithm in the formula both sides in step 5) and can obtain ln [Fq(s)]=H (q) ln (s)+c(c is constant), thus can obtain
Obtain slope H (q) of straight line.
Spectrum kurtosis method in step 7 comprises the following steps:
1) one cut-off frequency of structure is fcLow pass filter h (n) of=0.125+ ε；ε > 0, f in this examplec=0.3;
2) based on the quasi-low pass filter h that h (n) structure passband is [0,0.25]0N () and passband are [0.25,0.5]
Quasi-high pass filter h1(n),
;
3) signal ci kN () is through h0(n)、 h1N () resolves into low frequency part c after filtering and being down-sampled2i k+1(n) and HFS
c2i+1 k+1N (), the down-sampled factor is 2, then shaping filter tree after successive ignition filters, and kth layer has 2kIndividual frequency band, wherein
ci kThe output signal of the i-th wave filter on kth layer in (n) expression wave filter tree, i=0 ..., 2k-1,0≤k≤K-1, this example
Middle K=8；c0N () represents x in step 7 described in claim 1f1(k)；
4) the mid frequency f of the i-th wave filter on kth layer in decomposition treekiAnd bandwidth BkIt is respectively
;
5) each filter results c is calculatedi k(n)( i=0,…, 2k-1) kurtosis；
6) all of spectrum kurtosis is collected, obtain the spectrum kurtosis that signal is total.
Rational spline iteration smoothed envelope in step 9 is analyzed method and is comprised the following steps:
1) signal calculatedz(k) absolute valuez(k) local extremum；In the 1st iteration,z(k) represent claim 1 institute
State x in step 9f2(k)；
2) rational spline curve matching Local Extremum is used to obtain envelope eov1(k)；
3) rightz(k) be normalized and obtain；
4) the 2nd iteration:z 1(k) again as new data, repeat above-mentioned steps 1) ~ 3), obtain；
5) ith iteration:z i-1(k) again as new data, repeat above-mentioned steps 1) ~ 3), obtain；
6) ifnSecondary iteration obtainsz n (k) amplitude less than or equal to 1, then iterative process stops, and finally obtains signalz
(k) envelope be。
In step (2), Empirical Mode Decomposition Algorithm comprises the following steps:
1) first screening process: find out data x(k respectively) upper and lower Local Extremum, use cubic spline curve respectively
The upper and lower Local Extremum of matching, obtains signal x(k) local maximum envelope and local minimum envelope, then by this two
The value of the respective points of bar envelope is averaged, and obtains averaged curve m1；
Seek signal x(k again) and this averaged curve m1Difference, i.e. h10=x(k)-m1, so far first screening process terminates；
X(k) x in step described in claim 2 (2) is representedj(k)；
2) second screening process: h10Again new data, repeat the above steps 1 it are taken as), available h11= h10-m11, here
Parameter m11Represent h10Mean curve, repeat this process j time, < when 0.3, screening process stops, here until 0.2 < SD, now, h1j= h1(j-1)-m1j, at this moment it is believed that h1jMode function is grasped in being one
(Intrinsic Mode Function, IMF), defining the 1st IMF is c1= h1j；
3) from x(k) deduct c1, r can be obtained1=x(k)-c1, then by r1As new data, and repeat above-mentioned two step operations, this
Sample can obtain the 2nd IMF；
4) the available a series of IMF of step 3) operation are repeated, if rnHave changed into a monotonous curve, then screening process is stopped
Only, primary signal is decomposed into following form the most at last:。
Test 1, utilizes the bearing vibration data with inner ring fault to test the performance of algorithm of the present invention
Card.
Experiment bearing used is 6205-2RS JEM SKF, utilizes electric discharge machining method working depth on bearing inner race
For 0.2794mm, width be the groove of 0.3556mm to simulate bearing inner race fault, this experiment load is about 0.7457kW, drives
Motor turns frequency and is about 29.5Hz, and bearing inner race fault characteristic frequency is about 160Hz, and sample frequency is 4.8KHz, during signal sampling
A length of 1s.
The inner ring fault-signal collected is as shown in Figure 4.
Initially with traditional envelope Analysis Method, the signal shown in Fig. 4 is analyzed, the analysis result obtained such as Fig. 5
Shown in.From fig. 5, it can be seen that the fault signature of bearing is blanked completely, the most traditional envelope Analysis Method can not be effectively
Extract the fault signature of bearing；Additionally, from fig. 5, it can be seen that the left end point of envelope spectrum exists abnormal high level, this explanation is by tradition
There is end effect in the envelope spectrum that method obtains.
Use method proposed by the invention that signal shown in Fig. 4 is analyzed, the analysis result obtained such as Fig. 6 institute
Show.From fig. 6, it can be seen that the spectral line corresponding to 160Hz and 320Hz is apparently higher than other spectral line, the two frequency correspondence respectively
1 frequency multiplication of bearing inner race fault characteristic frequency and 2 frequencys multiplication, may determine that bearing has inner ring fault accordingly；Can from Fig. 6
Go out, the present invention envelope spectrum obtained does not has end effect.
Showing through many experiments, in the case of load and fault dimensional depth are constant, the present invention can reliable recognition
Minimum inner ring fault dimension width is about 0.19mm, and traditional method can reliable recognition minimum inner ring fault dimension width about
For 0.53mm, precision improves 64.2%.
Test 2, utilizes the bearing vibration data with outer ring fault to test the performance of algorithm of the present invention
Card.
Experiment bearing used is 6205-2RS JEM SKF, utilizes electric discharge machining method working depth on bearing outer ring
For 0.2794mm, width be the groove of 0.5334mm to simulate bearing outer ring fault, this experiment load is about 2.237 kW, drives
Motor turns frequency and is about 28.7Hz, and bearing outer ring fault characteristic frequency is about 103Hz, and sample frequency is 4.8KHz, during signal sampling
A length of 1s.
The outer ring fault-signal collected is as shown in Figure 7.
Initially with traditional envelope Analysis Method, the signal shown in Fig. 7 is analyzed, the analysis result obtained such as Fig. 8
Shown in.From figure 8, it is seen that the fault signature of bearing is blanked completely, the most traditional envelope Analysis Method can not be effectively
Extract the fault signature of bearing；Additionally, from figure 8, it is seen that the left end point of envelope spectrum exists abnormal high level, this explanation is by tradition
There is end effect in the envelope spectrum that method obtains.
Use method proposed by the invention that signal shown in Fig. 7 is analyzed, the analysis result obtained such as Fig. 9 institute
Show.From fig. 9, it can be seen that the spectral line corresponding to 103Hz and 206Hz is apparently higher than other spectral line, the two frequency correspondence respectively
1 frequency multiplication of bearing outer ring fault characteristic frequency and 2 frequencys multiplication, may determine that bearing has outer ring fault accordingly；Can from Fig. 9
Go out, the present invention envelope spectrum obtained does not has end effect..
Showing through many experiments, in the case of load and fault dimensional depth are constant, the present invention can reliable recognition
Minimum outer ring fault dimension width is about 0.29mm, and traditional method can reliable recognition minimum outer ring fault dimension width about
For 0.68mm, precision improves 57.4%.
According to result of the test, think after analysis:
1) traditional envelope Analysis Method directly carries out Envelope Analysis to primary signal, or to after merely through simple process
Primary signal carries out Envelope Analysis, and different from traditional envelope Analysis Method, primary signal is entered by the present invention first with EEMD
Row decomposes, and then utilizes the rearrangement of data and substitutes operation eliminating noise therein and trend component, only stick signal component
In useful component, thus avoid the impact on Envelope Analysis result of noise and trend component, improve accuracy and accurately
Degree.
2) traditional envelope Analysis Method is transformed to basis with Hilbert, and Hilbert conversion requires analyzed letter
Number must be the narrow band signal of simple component, otherwise the frequency modulating section of signal will pollute the Envelope Analysis result of signal, but
It is the signal the most to be analyzed condition that the most strictly meets simple component and arrowband, so may result in prior art because of precision not
High and erroneous judgement problem easily occur, different from tradition envelope Analysis Method, the present invention utilizes rational spline iteration smoothed envelope to divide
Signal envelope is kept completely separate by analysis method with frequency modulating section, it is possible to avoid frequency modulating section to signal envelope analysis result
Impact, thus improve the precision of Envelope Analysis.
3) fault type of rotating machinery can be detected exactly.
4) there is end effect in the envelope spectrum obtained by traditional method, and the envelope spectrum obtained by the present invention it can be avoided that
End effect.
5) each step effect:
1st) step: gather vibration signal；
2nd) step: primary signal is resolved into the form of different component sum, some of which component correspondence noise and trend term, some
Component correspondence useful signal；
3rd) ~ 5) step: the signal that obtains above-mentioned decomposition performs reordering operations and substitutes operation, reject noise component(s) therein and
Trend term, only retains useful signal；
6th) step: remaining useful signal is sued for peace, should and as signal rearranged and substitute filtered result xf1(k)；
7th) step: to filtered signal xf1K () performs spectrum kurtosis analysis, obtain center frequency corresponding at signal maximum kurtosis
Rate f0And bandwidth B；
8th) step: according to mid frequency f0With bandwidth B to xf1K () carries out bandpass filtering, obtain signal xf2(k)；
9th) step: signal calculated xf2Envelope eov (k) of (k)；
10th) step: eov (k) is performed discrete Fourier transform and obtains envelope spectrum, judge the failure classes of bearing according to envelope spectrum
Type.
One skilled in the art would recognize that above-mentioned detailed description of the invention is exemplary, be 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
According to technical solution of the present invention improvements introduced, each fall within protection scope of the present invention.
Claims (8)
1. EEMD and the rational spline smoothed envelope of a rotating machinery analyzes method, it is characterised in that comprise the following steps:
Step 1: utilize acceleration transducer to measure the vibration signal x(k of rotating machinery with sample frequency fs), (k=1,2,
..., N), N is the length of sampled signal；
Step 2: use set empirical mode decomposition (Ensemble Empirical Mode Decomposition, EEMD) to calculate
Method is by signal x(k) resolve into n component and a trend term sum, i.e., wherein, ciK () represents
The i-th component obtained by EEMD algorithm, rnK () represents the trend term obtained by EEMD algorithm；
Step 3: to ciK () performs reordering operations and substitutes operation, data c that rearranged operation obtainsi shuffleK () represents,
Data c are obtained after substituting operationi FTranK () represents；
Step 4: to ci(k), ci shuffle(k) and ci FTranK () performs multi-fractal respectively and removes trend fluction analysis
(Multifractal Detrended Fluctuation Analysis, MFDFA), obtains generalized Hurst index curve, ci
The generalized Hurst index curve H of (k)iQ () represents；ci shuffleThe generalized Hurst index curve H of (k)i shuffle(q) table
Show；ci FTranThe generalized Hurst index curve H of (k)i FTranQ () represents；
Step 5: if Hi(q) and Hi shuffle(q) or Hi(q) and Hi FTranQ the relative error between () is less than 5%, or Hi
(q), Hi shuffle(q) and Hi FTranQ () three do not change with q, then abandon the c of correspondencei(k) component；
Step 6: to remaining ci(k) component sue for peace, by this and be designated as signal rearranged and substitute filtered result xf1(k)；
Step 7: to xf1K () performs spectrum kurtosis analysis, obtain the mid frequency f corresponding to signal kurtosis maximum0And bandwidth B；
Step 8: according to mid frequency f0With bandwidth B to xf1K () carries out bandpass filtering, obtain xf2(k)；
Step 9: to signal xf2K () performs rational spline iteration smoothed envelope and analyzes, obtain signal envelope eov(k)；
Step 10: the signal envelope eov(k to obtaining) perform discrete Fourier transform obtain envelope spectrum, according to envelope spectrum feature
Frequency judges the fault type of machine.
The EEMD of a kind of rotating machinery the most according to claim 1 and rational spline smoothed envelope analyze method, its feature
It is, described step 2 is gathered Empirical Mode Decomposition Algorithm and comprises the following steps:
(1) to data x0K () is added white noise sequence and is produced a new data xj(k):
Std[x0(k)] represent data x0The standard deviation of (k), wnjK () represents wnjIn kth data, wnjRepresent jth random
The white noise sequence produced, wnjAmplitude is 1,1≤j≤K；x0K () represents x (k) in step 2 described in claim 1；
(2) to xjK () performs empirical mode decomposition, obtain n component and a trend term
cijK () represents xjK () performs the i-th component that empirical mode decomposition obtains, rnjK () represents xjK () performs experience
The trend term that Mode Decomposition obtains；
(3) meansigma methods of K decomposition result is calculated
ciK () represents x0K () carries out gathering the i-th component that empirical mode decomposition obtains, rnK () represents x0K () collects
Close the trend term that empirical mode decomposition obtains.
The EEMD of a kind of rotating machinery the most according to claim 1 and rational spline smoothed envelope analyze method, its feature
Being, in described step 3, data rearrangement operation comprises the following steps:
Upset component c at randomiPutting in order of (k).
The EEMD of a kind of rotating machinery the most according to claim 1 and rational spline smoothed envelope analyze method, its feature
It is: in described step 3, the operation of data replacement comprises the following steps:
1) to component ciK () performs discrete Fourier transform, it is thus achieved that component ciThe phase place of (k)；
2) component c is replaced with one group of pseudo-independent same distribution number being positioned in (-π, π) intervaliThe original phase of (k)；
3) frequency domain data after phase place substitutes is performed inverse discrete Fourier transform and obtain data ci IFFTK (), asks for data
ci IFFTThe real part of (k).
The EEMD of a kind of rotating machinery the most according to claim 1 and rational spline smoothed envelope analyze method, its feature
It is: in described step 4, MFDFA method comprises the following steps:
1) structure x (k) (k=1,2 ..., N) profile Y (i):
X (k) represents the c in step 4 described in claim 1i(k) or ci shuffle(k) or ci FTran(k)；
2) signal profile Y (i) is divided into nonoverlapping NSSegment length is the data of s, divides with identical length from the opposite direction of data
Section, obtains 2NSSegment data；
3) utilize the polynomial trend of the every segment data of least square fitting, then calculate the variance of every segment data:
yvI () is the trend of the v segment data of matching, if the polynomial trend of matching is m rank, then remember that this goes the trend process to be
(MF-) DFAm；
4) meansigma methods of q rank wave function is calculated:
；
5) if x (k) exists self-similarity characteristics, then meansigma methods F of q rank wave functionqPower is there is between (s) and time scale s
Rule relation:
As q=0, the formula in step 4) dissipates, and at this moment H (0) is determined by logarithmic mean process defined in following formula:
6) are taken the logarithm in the formula both sides in step 5) and can obtain ln [Fq(s)]=H (q) ln (s)+c(c is constant), thus can obtain
Obtain slope H (q) of straight line.
The EEMD of a kind of rotating machinery the most according to claim 1 and rational spline smoothed envelope analyze method, its feature
It is: the spectrum kurtosis method in described step 7 comprises the following steps:
1) one cut-off frequency of structure is fcLow pass filter h (n) of=0.125+ ε；
2) based on the quasi-low pass filter h that h (n) structure passband is [0,0.25]0N () and passband are [0.25,0.5]
Quasi-high pass filter h1(n),
;
3) signal ci kN () is through h0(n)、 h1N () resolves into low frequency part c after filtering and being down-sampled2i k+1(n) and HFS c2i +1 k+1N (), the down-sampled factor is 2, then shaping filter tree after successive ignition filters, and kth layer has 2kIndividual frequency band, wherein ci k
The output signal of the i-th wave filter on kth layer in (n) expression wave filter tree, i=0 ..., 2k-1,0≤k≤K-1；
4) the mid frequency f of the i-th wave filter on kth layer in decomposition treekiAnd bandwidth BkIt is respectively
;
5) each filter results c is calculatedi k(n)( i=0,…, 2k-1) kurtosis；
6) all of spectrum kurtosis is collected, obtain the spectrum kurtosis that signal is total.
The EEMD of a kind of rotating machinery the most according to claim 1 and rational spline smoothed envelope analyze method, its feature
Being, the rational spline iteration smoothed envelope in described step 9 is analyzed method and is comprised the following steps:
1) signal calculatedz(k) absolute valuez(k) local extremum；In the 1st iteration,z(k) represent claim 1 institute
State x in step 9f2(k)；
2) rational spline curve matching Local Extremum is used to obtain envelope eov1(k)；
3) rightz(k) be normalized and obtain；
4) the 2nd iteration:z 1(k) again as new data, repeat above-mentioned steps 1) ~ 3), obtain；
5) ith iteration:z i-1(k) again as new data, repeat above-mentioned steps 1) ~ 3), obtain；
6) ifnSecondary iteration obtainsz n (k) amplitude less than or equal to 1, then iterative process stops, and finally obtains signalz
(k) envelope be。
The EEMD of a kind of rotating machinery the most according to claim 2 and rational spline smoothed envelope analyze method, its feature
Being, in described step (2), Empirical Mode Decomposition Algorithm comprises the following steps:
1) first screening process: find out data x(k respectively) upper and lower Local Extremum, use cubic spline curve respectively
The upper and lower Local Extremum of matching, obtains signal x(k) local maximum envelope and local minimum envelope, then by this two
The value of the respective points of bar envelope is averaged, and obtains averaged curve m1；
Seek signal x(k again) and this averaged curve m1Difference, i.e. h10=x(k)-m1, so far first screening process terminates；
X(k) x in step described in claim 2 (2) is representedj(k)；
2) second screening process: h10Again new data, repeat the above steps 1 it are taken as), available h11= h10-m11, here
Parameter m11Represent h10Mean curve, repeat this process j time, < when 0.3, screening process stops, here until 0.2 < SD, now, h1j= h1(j-1)-m1j, at this moment it is believed that h1jMode function is grasped in being one
(Intrinsic Mode Function, IMF), defining the 1st IMF is c1= h1j；
3) from x(k) deduct c1, r can be obtained1=x(k)-c1, then by r1As new data, and repeat above-mentioned two step operations, so
The 2nd IMF can be obtained；
4) the available a series of IMF of step 3) operation are repeated, if rnHave changed into a monotonous curve, then screening process stops,
Primary signal is decomposed into following form the most at last:。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610492064.7A CN105973603B (en) | 2016-06-29 | 2016-06-29 | The EEMD and rational spline smoothed envelope analysis method of a kind of rotating machinery |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610492064.7A CN105973603B (en) | 2016-06-29 | 2016-06-29 | The EEMD and rational spline smoothed envelope analysis method of a kind of rotating machinery |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105973603A true CN105973603A (en) | 2016-09-28 |
| CN105973603B CN105973603B (en) | 2018-03-13 |
Family
ID=57019472
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610492064.7A Expired - Fee Related CN105973603B (en) | 2016-06-29 | 2016-06-29 | The EEMD and rational spline smoothed envelope analysis method of a kind of rotating machinery |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105973603B (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108288058A (en) * | 2018-04-12 | 2018-07-17 | 大连理工大学 | A kind of improved wavelet threshold knee joint swinging signal Denoising Algorithm |
| CN108399368A (en) * | 2018-01-31 | 2018-08-14 | 中南大学 | A kind of artificial source's electromagnetic method observation signal denoising method |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2013088431A (en) * | 2011-10-13 | 2013-05-13 | General Electric Co <Ge> | Method and system for automatically detecting rolling bearing fault |
| CN103900815A (en) * | 2014-04-02 | 2014-07-02 | 兰州交通大学 | Rolling bearing fault diagnosis method based on EEMD and distribution fitting testing |
| CN104897403A (en) * | 2015-06-24 | 2015-09-09 | 北京航空航天大学 | Self-adaption fault diagnosis method based on permutation entropy (PE) and manifold-based dynamic time warping (MDTW) |
-
2016
- 2016-06-29 CN CN201610492064.7A patent/CN105973603B/en not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2013088431A (en) * | 2011-10-13 | 2013-05-13 | General Electric Co <Ge> | Method and system for automatically detecting rolling bearing fault |
| CN103900815A (en) * | 2014-04-02 | 2014-07-02 | 兰州交通大学 | Rolling bearing fault diagnosis method based on EEMD and distribution fitting testing |
| CN104897403A (en) * | 2015-06-24 | 2015-09-09 | 北京航空航天大学 | Self-adaption fault diagnosis method based on permutation entropy (PE) and manifold-based dynamic time warping (MDTW) |
Non-Patent Citations (1)
| Title |
|---|
| 林近山等: "多重分形去趋势波动分析在滚动轴承损伤程度识别中的应用", 《中国机械工程》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108399368A (en) * | 2018-01-31 | 2018-08-14 | 中南大学 | A kind of artificial source's electromagnetic method observation signal denoising method |
| CN108399368B (en) * | 2018-01-31 | 2021-08-20 | 中南大学 | Artificial source electromagnetic observation signal denoising method |
| CN108288058A (en) * | 2018-04-12 | 2018-07-17 | 大连理工大学 | A kind of improved wavelet threshold knee joint swinging signal Denoising Algorithm |
| CN108288058B (en) * | 2018-04-12 | 2022-03-29 | 大连理工大学 | Improved wavelet threshold knee joint swing signal denoising algorithm |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105973603B (en) | 2018-03-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106198015A (en) | The VMD of a kind of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN106198013A (en) | A kind of envelope Analysis Method based on empirical mode decomposition filtering | |
| 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 | |
| CN106153339B (en) | A kind of envelope Analysis Method based on the filtering of variation Mode Decomposition | |
| CN106096198A (en) | A kind of envelope Analysis Method based on variation Mode Decomposition with spectrum kurtosis | |
| CN106168538A (en) | The ITD of a kind of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN106198009A (en) | The EMD of a kind 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 | |
| CN106096199A (en) | The WT of a kind of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN106053069A (en) | SSD, spectral kurtosis and smooth iteration envelope analysis method of antifriction bearing | |
| CN106053060B (en) | A kind of envelope Analysis Method that filtering is decomposed based on nonlinear model | |
| CN106198012B (en) | A kind of envelope Analysis Method for decomposing and composing kurtosis based on local mean value | |
| CN105973603A (en) | EEMD and rational spline smooth envelope analysis method for rotating machine | |
| CN106198010B (en) | A kind of envelope Analysis Method that filtering is decomposed based on local mean value | |
| CN106153333A (en) | A kind of envelope Analysis Method based on wavelet decomposition filtering | |
| CN106096201A (en) | The EEMD of a kind of rotating machinery and smoothed cubic spline envelope Analysis Method | |
| CN106198018A (en) | The EEMD of a kind of rotating machinery and smooth iteration envelope Analysis Method | |
| CN106198014A (en) | A kind of envelope Analysis Method based on empirical mode decomposition with spectrum kurtosis | |
| 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 | |
| CN106053059B (en) | It is a kind of based on it is interior grasp time scale decompose filtering envelope Analysis Method | |
| CN106198017A (en) | The LMD of a kind of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method | |
| CN106124200A (en) | The ELMD of a kind of rotating machinery and smooth iteration envelope Analysis Method | |
| CN106198016A (en) | The NMD of a kind of rolling bearing, spectrum kurtosis and smooth iteration 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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20180313 Termination date: 20210629 |