CN106096313B - A kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis - Google Patents

A kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis Download PDF

Info

Publication number
CN106096313B
CN106096313B CN201610492078.9A CN201610492078A CN106096313B CN 106096313 B CN106096313 B CN 106096313B CN 201610492078 A CN201610492078 A CN 201610492078A CN 106096313 B CN106096313 B CN 106096313B
Authority
CN
China
Prior art keywords
signal
data
envelope
spectrum
frequency
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.)
Expired - Fee Related
Application number
CN201610492078.9A
Other languages
Chinese (zh)
Other versions
CN106096313A (en
Inventor
林近山
窦春红
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Weifang University
Original Assignee
Weifang University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Weifang University filed Critical Weifang University
Priority to CN201610492078.9A priority Critical patent/CN106096313B/en
Publication of CN106096313A publication Critical patent/CN106096313A/en
Application granted granted Critical
Publication of CN106096313B publication Critical patent/CN106096313B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16ZINFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS, NOT OTHERWISE PROVIDED FOR
    • G16Z99/00Subject matter not provided for in other main groups of this subclass

Landscapes

  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis that the invention discloses a kind of, this method decomposes original signal first with singular spectrum decomposition 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

A kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis
Technical field
The present invention relates to condition monitoring for rotating machinery and fault diagnosis field, and in particular to one kind based on unusual spectral factorization and Compose the envelope Analysis Method of kurtosis.
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 a kind of envelope based on unusual spectral factorization and spectrum kurtosis Analysis method has analysis result accuracy and accuracy high, and can accurately examine after envelope Analysis Method of the invention The advantages of measuring rotating machinery fault type.
In order to solve the above technical problems, the technical solution adopted by the present invention is as follows: one kind based on unusual spectral factorization and compose it is high and steep The envelope Analysis Method of degree, 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 is the length of sampled signal;
Step 2: using unusual spectral factorization (Singular Spectrum Decomposition) algorithm by signal x(k) The sum of n component is resolved into, i.e.,, wherein ci(k) the obtained by singular spectrum decomposition algorithm is represented I component, unusual spectral factorization is it is well known that be shown in document
P. Bonizzi, J.M. KAREL, O. Meste, R.L. Peeters, Singular spectrum decomposition: A new method for time series decomposition, Advances in Adaptive Data Analysis, 2014,6 (4): 1-29;
Step 3: to ci(k) reordering operations and substitution operation are executed, it is rearranged to operate obtained data ci shuffle(k) table Show, obtains data c 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 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 executed, centre frequency f corresponding to signal kurtosis maximum is found out0And band Wide 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 Characteristic frequency judges the fault type of machine.
A kind of prioritization scheme, data rearrangement operation in the step 3 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 Data ci IFFT(k) real part.
Further, in the step 4 MFDFA method the following steps are included:
1) the profile Y (i) of x (k) (k=1,2 ..., N) is constructed:
,
X (k) represents the c in the step 4i(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, since data length N generally can not be whole Except s, so the remaining one piece of data of meeting cannot utilize;
In order to make full use of the length of data, then it is obtained with identical length segmentation, such one from the opposite direction of data 2NSSegment data;
3) it is fitted the polynomial trend of every segment data using least square method, then calculates the variance of every segment data:
yv(i) trend for the v segment data of fitting remembers that this goes trend process if the polynomial trend of fitting is m rank For (MF-) DFAm;In this example, m=1;
4) average value of q rank wave function is calculated:
5) if there are self-similarity characteristics, the average value F of q rank wave function by x (k)q(s) between time scale s There are power law relations:
As q=0, the formula in step 4) dissipates, and at this moment H (0) is by logarithmic mean process defined in following formula come really It is fixed:
6) take logarithm that can obtain ln [F on the formula both sides in step 5)q(s)]=H (q) ln (s)+c(c is constant), thus may be used To obtain the slope H (q) of straight line.
Further, the spectrum kurtosis method in the step 7 the following steps are included:
1) constructing a cutoff frequency is fcThe low-pass filter h (n) of=0.125+ ε;ε > 0, f in this examplec=0.3;
It 2) is the quasi- low-pass filter h of [0,0.25] based on h (n) construction passband0(n) be with passband [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 radio-frequency head Divide 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, Middle 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 K=8 in example;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 the step X in 9f2(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 Signalz(k) envelope 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 unusual spectral factorization, then excludes it using the rearrangement of data and substitution operation In noise and trend component, the only useful component in stick signal component, so as to avoid noise and trend component to packet Network analyzes the influence of result, and 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 to carry out preliminary exposition to signal using low-pass filter and high-pass filter in the embodiment of the present invention to show It is intended to;
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 using traditional envelope Analysis Method in the embodiment of the present invention to inner ring faulty bearing vibration signal Analyze 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 using traditional envelope Analysis Method in the embodiment of the present invention to outer ring faulty bearing vibration signal Analyze 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, a kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis, packet Include 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 is the length of sampled signal;
Step 2: using unusual spectral factorization (Singular Spectrum Decomposition) algorithm by signal x(k) Resolve into the sum of n component, that is, wherein, ci(k) i-th of the component obtained by singular spectrum decomposition algorithm, singular spectrum point are represented Solution is it is well known that be shown in document
P. Bonizzi, J.M. KAREL, O. Meste, R.L. Peeters, Singular spectrum decomposition: A new method for time series decomposition, Advances in Adaptive Data Analysis, 2014,6 (4): 1-29;
Step 3: to ci(k) reordering operations and substitution operation are executed, it is rearranged to operate obtained data ci shuffle(k) table Show, obtains data c 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 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 executed, centre frequency f corresponding to signal kurtosis maximum is found out0And band Wide 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 Characteristic frequency judges the fault type of machine.
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 Data ci IFFT(k) real part.
MFDFA method in step 4 the following steps are included:
1) the profile Y (i) of x (k) (k=1,2 ..., N) is constructed:
,
,
X (k) represents the c in the step 4i(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, since data length N generally can not be whole Except s, so the remaining one piece of data of meeting cannot utilize;
In order to make full use of the length of data, then it is obtained with identical length segmentation, such one from the opposite direction of data 2NSSegment data;
3) it is fitted the polynomial trend of every segment data using least square method, then calculates the variance of every segment data:
yv(i) trend for the v segment data of fitting remembers that this goes trend process if the polynomial trend of fitting is m rank For (MF-) DFAm;In this example, m=1;
4) average value of q rank wave function is calculated:
5) if there are self-similarity characteristics, the average value F of q rank wave function by x (k)q(s) between time scale s There are power law relations:
As q=0, the formula in step 4) dissipates, and at this moment H (0) is by logarithmic mean process defined in following formula come really It is fixed:
6) take logarithm that can obtain ln [F on the formula both sides in step 5)q(s)]=H (q) ln (s)+c(c is constant), thus may be used To obtain the slope H (q) of straight line.
Spectrum kurtosis method in step 7 the following steps are included:
1) constructing a cutoff frequency is fcThe low-pass filter h (n) of=0.125+ ε;ε > 0, f in this examplec=0.3;
It 2) is the quasi- low-pass filter h of [0,0.25] based on h (n) construction passband0(n) be with passband [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 radio-frequency head Divide 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, Middle 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 K=8 in example;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 the step X in 9f2(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 Signalz(k) envelope be
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.23 mm, and conventional method is capable of the minimum inner ring failure dimension width of reliable recognition About 0.53mm, precision improve 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 merely through simple process Original signal afterwards carries out Envelope Analysis, and different from traditional envelope Analysis Method, the invention firstly uses unusual spectral factorizations pair Original signal is decomposed, and then excludes noise and trend component therein, Jin Jinbao using the rearrangement of data and substitution operation The useful component in signal component is stayed, the influence so as to avoid noise and trend component to Envelope Analysis result improves standard Exactness 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 components correspond to useful signal;
The 3) ~ 5) step: reordering operations are executed to the signal that above-mentioned decomposition obtains and substitution operates, reject noise therein point Amount and trend term only retain 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 at signal maximum kurtosis is found out Frequency 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 of bearing is 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. a kind of envelope Analysis Method based on unusual spectral factorization 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: being decomposed signal x(k) using unusual spectral factorization (Singular Spectrum Decomposition) algorithm At the sum of n component, i.e.,, wherein ci(k) i-th point obtained by singular spectrum decomposition algorithm is represented Amount, unusual spectral factorization are known;
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.
2. a kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis according to claim 1, which is characterized in that In the step 3 data rearrangement operation the following steps are included:
Upset component c at randomi(k) put in order.
3. a kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis according to claim 1, it is characterised in that: 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.
4. a kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis according to claim 1, it is characterised in that: 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;m=1;
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 time scalesBetween exist Power law relation:
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)。
5. a kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis according to claim 1, it is characterised in that: Spectrum kurtosis method in the step 7 the following steps are included:
1) constructing a cutoff frequency is fcThe low-pass filter h (n) of=0.125+ ε;ε > 0, fc=0.3;
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=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.
6. a kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis according to claim 1, which is characterized in that 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 signalz (k) envelope be
CN201610492078.9A 2016-06-29 2016-06-29 A kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis Expired - Fee Related CN106096313B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610492078.9A CN106096313B (en) 2016-06-29 2016-06-29 A kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610492078.9A CN106096313B (en) 2016-06-29 2016-06-29 A kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis

Publications (2)

Publication Number Publication Date
CN106096313A CN106096313A (en) 2016-11-09
CN106096313B true CN106096313B (en) 2018-12-25

Family

ID=57215279

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610492078.9A Expired - Fee Related CN106096313B (en) 2016-06-29 2016-06-29 A kind of envelope Analysis Method based on unusual spectral factorization and spectrum kurtosis

Country Status (1)

Country Link
CN (1) CN106096313B (en)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107560718B (en) * 2017-07-25 2020-01-10 苏州微著设备诊断技术有限公司 Deconvolution method based on autocorrelation kurtosis maximization
CN108489719B (en) * 2018-04-09 2019-11-05 常州湖南大学机械装备研究院 A kind of combined failure of rotating machinery diagnostic method based on the unusual spectral factorization of G-P
CN108731945B (en) * 2018-08-02 2019-12-13 南昌航空大学 Method for extracting fault signal characteristic information of aircraft engine rotor system
CN109141885B (en) * 2018-09-04 2020-04-10 厦门理工学院 Bearing vibration signal envelope demodulation method, device and equipment based on MRSVD
CN109187023B (en) * 2018-09-04 2021-01-26 温州大学激光与光电智能制造研究院 Automobile generator bearing fault diagnosis method
CN109612726A (en) * 2018-11-20 2019-04-12 南京航空航天大学 A kind of multiple superstage analysis method extracted for vibration signal characteristics
CN112305380A (en) * 2020-09-01 2021-02-02 华南理工大学 Partial discharge white noise suppression method based on S transformation and spectral kurtosis
CN112595515A (en) * 2020-12-04 2021-04-02 中国船舶工业综合技术经济研究院 Power shafting bearing fault detection method and system

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104036493A (en) * 2014-05-21 2014-09-10 浙江大学 No-reference image quality evaluation method based on multifractal spectrum
CN104634571A (en) * 2015-02-06 2015-05-20 北京航空航天大学 Fault diagnosis method for rolling bearing based on LCD-MF (Local Characteristic Scale Decomposition )-(Multifractal)
CN105129109A (en) * 2015-09-30 2015-12-09 北京航空航天大学 Method for evaluating health of aircraft aileron actuator system based on multi-fractal theory and self-organizing map (SOM) network

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104036493A (en) * 2014-05-21 2014-09-10 浙江大学 No-reference image quality evaluation method based on multifractal spectrum
CN104634571A (en) * 2015-02-06 2015-05-20 北京航空航天大学 Fault diagnosis method for rolling bearing based on LCD-MF (Local Characteristic Scale Decomposition )-(Multifractal)
CN105129109A (en) * 2015-09-30 2015-12-09 北京航空航天大学 Method for evaluating health of aircraft aileron actuator system based on multi-fractal theory and self-organizing map (SOM) network

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Fault diagnosis of rolling bearings based on multifractal detrended fluctuation analysis and Mahalanobis distance criterion;Jinshan Lin等;《Mechanical systems and signal processing》;20130731;第38卷(第2期);515-533 *
基于时间序列标度分析的旋转机械故障诊断方法研究;林近山;《中国博士学位论文全文数据库工程科技Ⅱ辑》;20141215;C029-24 *
基于经验模式分解和谱峭度的滚动轴承故障诊断;林近山;《机械传动》;20120915;第36卷(第9期);76-79 *

Also Published As

Publication number Publication date
CN106096313A (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
CN106096198B (en) A kind of envelope Analysis Method based on variation Mode Decomposition and spectrum kurtosis
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
CN105954031B (en) A kind of envelope Analysis Method based on unusual spectral factorization filtering
CN106168538A (en) The ITD of a kind of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method
CN106096199B (en) A kind of WT 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
CN106198009B (en) A kind of EMD 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
CN106198010B (en) A kind of envelope Analysis Method that filtering is decomposed based on local mean value
CN106153333B (en) A kind of envelope Analysis Method based on wavelet decomposition filtering
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
CN106053061B (en) A kind of envelope Analysis Method for decomposing and composing kurtosis based on nonlinear model
CN106198014B (en) A kind of envelope Analysis Method based on empirical mode decomposition and spectrum kurtosis
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
CN106198017B (en) A kind of LMD of rolling bearing, spectrum kurtosis and smooth iteration envelope Analysis Method
CN106198018A (en) The EEMD of a kind 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
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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20181225

Termination date: 20210629

CF01 Termination of patent right due to non-payment of annual fee