CN105021399B - A kind of feature extracting method based on single channel signal blind separation rolling bearing - Google Patents
A kind of feature extracting method based on single channel signal blind separation rolling bearing Download PDFInfo
- Publication number
- CN105021399B CN105021399B CN201510363665.3A CN201510363665A CN105021399B CN 105021399 B CN105021399 B CN 105021399B CN 201510363665 A CN201510363665 A CN 201510363665A CN 105021399 B CN105021399 B CN 105021399B
- Authority
- CN
- China
- Prior art keywords
- mrow
- signal
- frequency
- msub
- msup
- 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
Links
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Complex Calculations (AREA)
Abstract
The invention discloses a kind of feature extracting method based on single channel signal blind separation rolling bearing, frequency slice wavelet transformation is carried out to original vibration signal, the energy spectrum of signal is obtained, is no longer rely on the selection of wavelet basis, the selection of section also overcomes wavelet transformation frequency band and chooses the problem of limited.Number of source is determined using the method for principal component analysis, solves the problems, such as to obtain cluster conclusion difficulty when sample is larger.The vector project matrix of dimensionality reduction is obtained using projection pattern according to principal component analysis, avoids the openness influence of signal, so as to which underdetermined problem is converted into well-posed problem, in addition, also can effectively realize that fault signature extracts to the source signal with certain correlation.
Description
Technical field
The present invention relates to technology for mechanical fault diagnosis field, it is related to a kind of based on single channel signal blind separation rolling bearing
Feature extracting method.
Background technology
Blind Signal Separation is the hot issue of field of signal processing research in the last few years.So-called Blind Signal Separation just refer to from
Some sensors observe the method that the primary signal that can not be directly observed is recovered in the mixed signal of multiple signals.Fanaticism
Number separation is general requires that number of probes be more than signal source number, and single channel blind separation then refers to that being used for observation signal mixing believes
Number sensor there was only one.This is a difficult point in Blind Signal Separation, but this is in mechanical fault diagnosis engineer applied
It is more closing to reality condition.
For owing to determine situation of the blind separation i.e. number of sensors less than number of source, main method has poly- based on intermediate value
Class algorithm blind separating method[1], deficient based on potential function determine blind separating method[2], the problem of these are present is to be based on source signal
Sparse Problems, and then separating effect is bad for openness poor signal；The also fanaticism number point based on wavelet decomposition
From[3], the problem of it is present is that the selection dependence to wavelet basis is very strong, and wavelet basis difference separating effect difference is very big；In addition
Have based on intrinsic mode function (EMD), the blind separating method of overall intrinsic mode function (EEMD)[4-5], ask existing for EMD methods
Topic is modal overlap phenomenon be present, and EEMD method problems are that amount of calculation is larger, and the calculating time is also longer, therefore should in Practical Project
It is not strong with middle instantaneity.
Bibliography
【1】FABIAN J T；CARLOS G P；ELMAR W L Median-based clustering for
Underdetermined blind signal processing [foreign language periodicals] 2006 (02)
【2】Zhang Yun；Li Benwei；Wang Yonghua；Deficient based on potential function determines blind source separating identifying and diagnosing method [J] aviations
Power journal, 2010,25 (01), 218-223.
【3】Wang Jiao；Liu Yulin；He Wei；Chao Zhichao；Wavelet decomposition single channel blind separation disturbance restraining method [J]
Chongqing Mail and Telephones Unvi's journal (natural science edition), 2014,26 (5), 648-653.
【4】Li Shun is dead drunk；Liu Xiaowei；Guo Haidong；A kind of blind separating method [P] Chinese patents of single channel vibration signal：CN
102288285 B, 2011-05-24.
【5】Meng Zong；Cai Long；Vibrated and rushed based on EEMD subband extraction associated mechanical vibration signal single channel blind separation [J]
Hit, 2014,33 (20), 40-46.
The content of the invention
For above-mentioned problems of the prior art and defect, it is an object of the present invention to provide one kind to be based on single-pass
The feature extracting method of road blind signal separation rolling bearing.
To achieve these goals, the present invention adopts the following technical scheme that：
A kind of feature extracting method based on single channel signal blind separation rolling bearing, specifically includes following steps：
Step 1：Selecting frequency section wavelet transform function, frequency slice small echo change is carried out to given original vibration signal
Change, obtain time-frequency figure, then inverse transformation is carried out to the signal after frequency slice wavelet transformation and obtains reconstruction signal, and draw reconstruct letter
Number energy spectrum；
Step 2：The time-frequency figure and energy spectrum obtained according to step 1, select the multiple areas for including energy peak in energy spectrum
Between original vibration signal is cut into slices, obtain the reconstruction signal of multiple sections；
Step 3：The reconstruction signal that step 2 is obtained carries out noise reduction and de-redundancy, it is determined that optimal number of source m and its correspondingly
M characteristic vector；
Step 4：The reconstruction signal that the Matrix Multiplication of m characteristic vector composition is obtained with step 2, obtain the dimensions of the m after dimensionality reduction
Matrix；
Step 5：The m dimension matrixes that step 4 obtains, m separation signal is obtained using independent component analysis, and seek its bag respectively
Network is composed；Envelope spectrum is normalized, superposition obtains equivalent envelope spectrum；Envelope spectrum is observed, extracts fault signature.
Specifically, the concrete methods of realizing of the step 1 is as follows：
Frequency slice wavelet transformation definition is：
In formula, σ is scale factor, and σ ≠ 0, σ are constant or ω and t function,Or σ=kt etc., wherein, k is tune
It is time domain variable to save time domain or frequency domain sensitivity, t, and ω is frequency domain variable；P () is frequency slice function, p*() is p ()
Conjugate function；F (τ) is original vibration signal, f (τ) ∈ L2(R)；
Inverse transformation is carried out to frequency slice wavelet transform signal, frequency slice wavelet transform signal W (t, ω, σ) is obtained and exists
Time-frequency region (t1,t2,ω1,ω2) component of signal, i.e. reconstruction signal：
Energy definition formula is：In formula, E (ω) EDF, E (ω)=| F
(ω)|2, wherein F (ω) is reconstruction signal fa(t) Fourier transformation.
Specifically, the concrete methods of realizing of the step 3 is as follows：
Step 3.1：Seek the covariance matrix A=(A of reconstruction signalij)pxp, wherein
Wherein, AijFor the element in covariance matrix A, xikFor xFSWT(t) the i-th row k column elements,For xFSWT(t) the i-th row
Average value, n are data length；xjkFor xFSWT(t) jth row k column elements,For xFSWT(t) jth row average value；
Step 3.2：Calculate above-mentioned covariance matrix A eigenvalue λ1≥λ2≥λ3≥...λn＞ 0 and its orthogonal unit
Change characteristic vector from left to right to arrangeAccording to formula
Variance contribution ratio is calculated, when its value is more than setting value, determines number of source m and its corresponding m characteristic vector.
Compared with prior art, the present invention has following technique effect：
1st, the present invention carries out frequency slice wavelet transformation to original vibration signal, obtains the energy spectrum of signal, is no longer rely on
The selection of wavelet basis, the selection for frequency band of cutting into slices also overcome the problem of wavelet transformation frequency band is limited.
2nd, number of source is determined using the method for principal component analysis, solves and cluster conclusion difficulty is obtained when sample is larger
The problem of.
3rd, the vector project matrix of dimensionality reduction is obtained using projection pattern according to principal component analysis, it is openness avoids signal
Influence, so as to which underdetermined problem is converted into well-posed problem.
4th, the present invention also can effectively realize that fault signature extracts to the source signal with certain correlation.
5th, the present invention is compared with conventional method, and amount of calculation is smaller, and calculating speed is very fast.
Brief description of the drawings
Fig. 1 is flow chart of the method for the present invention；
Fig. 2 is the source signal and mixed signal time domain and spectrogram of embodiment 1；
Fig. 3 is the time-frequency figure and reconstruction signal energy spectrogram of the mixed signal of embodiment 1；
Fig. 4 is the source signal time domain and envelope spectrum oscillogram after the separation of embodiment 1；
Fig. 5 is the equivalent envelope spectrum oscillogram of normalization of embodiment 1；
Fig. 6 is the source signal time-domain diagram and envelope spectrogram after contrast experiment's separation of embodiment 1；
Fig. 7 is that the source signal after contrast experiment's separation of embodiment 1 normalizes equivalent envelope spectrum oscillogram；
Fig. 8 is the observation signal time domain beamformer of embodiment 2；
Fig. 9 is the observation signal time-frequency figure of embodiment 2 and the energy spectrogram of reconstruction signal；
Figure 10 is the observation signal time-domain diagram and envelope spectrogram after the separation of embodiment 2；
Figure 11 is the equivalent envelope spectrum oscillogram of normalization of the observation signal after the separation of embodiment 2.
Explanation and illustration in further detail is done to the solution of the present invention with reference to the accompanying drawings and examples.
Embodiment
Above-mentioned technical proposal is deferred to, referring to Fig. 1, the feature of the invention based on single channel signal blind separation rolling bearing carries
Method is taken, specifically includes following steps：
Step 1：Selecting frequency slice function, frequency slice wavelet transformation is carried out to given original vibration signal, obtained
Time-frequency figure in its 0~fs/2 frequency band, then inverse transformation is carried out to the signal after frequency slice wavelet transformation and obtains reconstruction signal, its
In, fs is sample frequency, and draws the energy spectrum of reconstruction signal.Its concrete methods of realizing is as follows：
Frequency slice wavelet transformation definition is：
In formula, σ is scale factor, and σ ≠ 0, σ are constant or ω and t function,Or σ=kt etc., wherein, k is tune
Time domain or frequency domain sensitivity are saved, its value is the compromise selection in time domain and frequency domain resolution, and t is time domain variable, and ω is frequency domain
Variable；P () is frequency slice function, and p* () is p () conjugate function；F (τ) is original vibration signal, f (τ) ∈ L2
(R), above-mentioned calculation is had no backing the selection to wavelet basis.
Because signal time domain and frequency domain behavior are not separate, thus frequency slice wavelet transform result is redundancy.
Its inverse transformation can have multi-form in theory.Inverse transformation is carried out to frequency slice wavelet transform signal, it is small to obtain frequency slice
Wave conversion signal W (t, ω, σ) is in time-frequency region (t1,t2,ω1,ω2) component of signal, i.e. reconstruction signal：
As can be seen from the above equation, by changing time-frequency region (t1,t2,ω1,ω2), it be free in time frequency space
Component needed for extraction, wherein, t1,t2For time value, ω1,ω2For frequency values.
Energy definition formula is：In formula, E (ω) EDF, i.e. cell frequency
Signal energy in 1Hz, E (ω)=| F (ω) |2, wherein F (ω) is reconstruction signal fa(t) Fourier transformation.
If reconstruction signal fa(t) it is discrete signal, then the drafting of energy spectrum is used as slicing band using certain frequency bandwidth factor
Width, according to the amount of bandwidth under Whole frequency band, with the moving window of a fixed step size, cut into slices, the selection of bandwidth is characterized frequency
1 to 5 frequencys multiplication of rate, so ensure that the section containing main information there can be higher energy in Envelope Analysis, step-length selection is got over
Small, then precision is higher, can be very long but calculate the time, therefore the selection of step-length should be precision and calculate the compromise selection of time.
The relative scale energy of each section, using normalization relative energy monitoring, the relative energy of m-th of section is:
In formula：En(fa(t)) it is signal fa(t) gross energy,Cut into slices for m-th,For m-th
The energy of section, En(m) it is the relative energy values of m-th of section.All section relative energy values are drawn out, obtain energy
Amount spectrum.
Step 2：The time-frequency figure and energy spectrum obtained according to step 1, select energy spectrum in comprising energy peak be energy divide
The higher multiple sections of cloth are cut into slices to original vibration signal, obtain the reconstruction signal of multiple sections.Its concrete methods of realizing
It is as follows：
Signal transducer is gathered in engineering and is generally LPF, and the energy of slicing band is by each frequency in frequency separation
Rate multicomponent energy form, therefore window move when, due to frequency separation change, Energy maximum value can unexpected appearing and subsiding, because
This according to LPF principle selection energy spectrum in rise it is quick and precipitous, that is, cut into slices section choose should be 0 arrive comprising energy
No more than one section bandwidth of difference of the frequency at each peak value in spectrum, the frequency and neighbouring peak value, then to original letter
Number cut into slices, obtain the reconstruction signal of the higher section of n energy, and formed a n dimension matrixes xFSWT(t)=(c1,
c2..., cn)T, wherein n >=3, and n is positive integer.
Step 3：The reconstruction signal that step 2 is obtained, i.e. n dimension matrix xFSWT(t)=(c1,c2..., cn)T, using based on
Principal component analysis (PCA) carries out noise reduction and de-redundancy, makes the correlation between the dimension that remains as small as possible, while its variance
Value is as big as possible, so that it is determined that optimal number of source.Solved compared with clustering methodology sample it is larger when, obtain cluster conclusion
The problem of difficult.Compared with singular value decomposition method, what principal component was selected is the vector after projection, and requirement openness to original signal is more
It is low.Its concrete methods of realizing is as follows：
Step 3.1：Seek the covariance matrix A=(A of reconstruction signalij)pxp, wherein Wherein AijFor the element in covariance matrix A, xikFor xFSWT(t) the i-th row k column elements,For xFSWT(t) i-th
Row average value, n are data length；xjkFor xFSWT(t) jth row k column elements,For xFSWT(t) jth row average value.
Step 3.2：Calculate above-mentioned covariance matrix A eigenvalue λ1≥λ2≥λ3≥...λn＞ 0 and its orthogonal unit
Change characteristic vector from left to right to arrangeAccording to formula
Variance contribution ratio is calculated, when its value is more than setting value, then it is assumed that the main information of primary signal has been covered, so that it is determined that information source
Number m and its corresponding m characteristic vector.
Step 4：The number of source m obtained according to step 3 and its corresponding m characteristic vector, then calculate m characteristic vector
The Matrix Multiplication of composition is n dimension matrixes x with reconstruction signalFSWT(t)=(c1,c2..., cn)T, the dimension matrixes of the m after dimensionality reduction are obtained, i.e.,
N dimension matrixes are projected under m-dimensional space, obtain the vector project matrix of m dimensions, are asked surely so as to which underdetermined problem to be converted into fit
Topic.
Step 5：The m dimension matrixes obtained to step 4, using independent component analysis (ICA), m separation signal is obtained, and divide
Its envelope spectrum is not sought；Envelope spectrum is normalized, superposition obtains equivalent envelope spectrum, i.e.,
Envelope spectrum is observed, fault signature is extracted, realizes single channel Blind Signal Separation.Wherein FFj(t) it is the envelope of j-th of signal, peakj
For its peak value.
Embodiment 1：
The feature extracting method based on single channel signal blind separation rolling bearing of the present embodiment, specifically includes following step
Suddenly：
Step 1：Give two original vibration signal s1And s2, it is respectively：
s1=cos (2 π f1t+π/3)
s2=cos (2 π fbt)[1+βcos(2πfrt)]
Wherein, f1=25Hz, fr=25Hz, fb=115Hz, β=2, sampling number 1024, sample frequency fs are
1000Hz, if mixed signal model s=as1(t)+bs2(t)+n (t), wherein, a=1, b=1, n (t) believe for random white noise
Number.Emulate signal s1、s2With s time-domain diagram and spectrogram, as shown in Figure 2.
Selecting frequency slice function, then frequency slice wavelet transformation is carried out to above-mentioned mixed signal s, obtain its 0~fs/2
Time-frequency figure in frequency band, then inverse transformation is carried out to the signal after frequency slice wavelet transformation and obtains reconstruction signal, and draw reconstruct
The energy spectrum of signal, as shown in figure 3, its concrete methods of realizing is as follows：
Frequency slice wavelet transformation definition is：For discrete letter
Number, in formula：σ is scale factor, σ ≠ 0, p ()=exp (- C2/ 2), p* ()=(exp (- C2/2))*, p* (σ (τ-t))=
(exp(-(σ(τ-t))2/2))*；σ=sqrt (2)/2/0.025.
Frequency slice wavelet transform function W (t, ω, σ) is calculated in time-frequency region (t1,t2,ω1,ω2) component of signal,
That is reconstruction signal：
Wherein, (t1,t2,ω1,ω2) it is (0,1,0,500)；
Energy definition formula is：In formula, E (ω) EDF, i.e. cell frequency 1Hz
Interior signal energy.For discrete signal, energy definition isN is signal length, and F (ω) is fa(t)
Fourier transformation, Fi(ω) is i-th point of F (ω), with a width of 50Hz, is characterized 2 times of frequency, step-length 1Hz.
Step 2：The time-frequency figure and energy spectrum obtained according to step 1, the higher multiple sections of Energy distribution are selected to original
Signal is cut into slices, and obtains the reconstruction signal of multiple sections.Its concrete methods of realizing is as follows：
Referring to Fig. 3, can be seen that from whole energy spectrum, energy rise is quick at 30Hz, 90Hz, 120Hz, 140Hz,
Its elsewhere energy is almost 0, therefore selects to be cut comprising each Energy distribution upper section according to LPF principle
Piece, [0 40], [0 100], [0 130], [0 150] are selected respectively, obtain the reconstruction signal of 4 sections, and formed one
Individual 4 dimension matrix xFSWT(t)=(c1,c2,c3,c4)T, form reconstruction signal.Test as a comparison, select four Energy distributions relatively low
Section carry out experimental analysis, respectively select [40 80], [100 120], [160 200], [300 400], obtain 4 reconstruct
Signal.
Step 3：The reconstruction signal x that step 2 is obtainedFSWT(t)=(c1,c2,c3,c4)T, using based on principal component analysis
(PCA) source number estimation method, calculates its characteristic value and characteristic vector, determines number of source m.Characteristic value is as shown in table 1, letter
Source number is 2.The source number estimation method based on principal component analysis (PCA) is also used for the reconstruction signal in contrast experiment,
Its characteristic value and characteristic vector are calculated, characteristic value is as shown in table 2, number of source 4.Its concrete methods of realizing is as follows：
Step 3.1：Seek the covariance matrix A=(A of reconstruction signalij)4x4, wherein N is data length 1024.
Step 3.2：Calculate above-mentioned covariance matrix A eigenvalue λ1≥λ2≥...λ4＞ 0 and its orthogonal unitization spy
Sign vector from left to right arrangesPass through formulaCalculate
Variance contribution ratio, when its value is more than 0.90, then it is assumed that primary signal full detail has been covered, so that it is determined that number of source is 2.
The multidimensional signal x of table 1FSWT(t) characteristic value
The multidimensional signal x of table 2FSWT(t) characteristic value
Step 4：The number of source m=2 obtained according to step 3, then the Matrix Multiplication of preceding 2 characteristic vectors composition is calculated with weight
Structure signal xFSWT(t)=(c1,c2,c3,c4)T, 2 dimension matrixes after dimensionality reduction are obtained, so as to which underdetermined problem is converted into well-posed problem.
By the source signal battle array in the contrast experiment also vector project under the space that its characteristic vector is formed, the dimension matrix of composition 4.
Step 5：To the new signal battle array being made up of in step 4 vector project, i.e., 2 dimension matrixes, using independent component analysis
(ICA) 2 separation signals, are obtained, and seek its envelope spectrum respectively, as shown in Figure 4；Envelope spectrum is normalized, is superimposed
Equivalent envelope spectrum is obtained, as shown in figure 5, i.e.Envelope spectrum is observed, extracts fault signature, it is real
Existing single channel Blind Signal Separation.As can be seen from Figure 2 source signal characteristics frequency is 25Hz, it is seen from figure 5 that signal after separation
Occur high peaks at 25Hz, this coincide with source signal characteristics frequency, time-domain signal and its envelope spectrum and figure after being separated in Fig. 4
Source signal and its envelope spectrum are also extremely similar in 2, and it is 0.9940,0.8061 to compare its coefficient correlation, and similarity is strong.Therefore this is blind
Separation method can realize that characteristic frequency is extracted.The 4 dimensional signal battle arrays obtained in contrast experiment are used into independent component analysis again, obtained
To 4 separation signals, and its envelope is sought respectively, as shown in Figure 6, it can be seen that signal is disorderly and unsystematic after separation, believes with source in Fig. 2
Number and its envelope it is completely dissimilar, then envelope spectrum is normalized, superposition obtains equivalent envelope spectrum, as shown in fig. 7, nothing
Method observes source signal characteristics frequency, illustrates infeasible in this way.
Embodiment 2
Blind separation is carried out to a certain mechanical bearing fault-signal using the method for the present invention, and extracts fault characteristic frequency,
The bearing fault is that bearing roller damages.Vibrating sensor is installed in the bearing drive end, sample frequency fs is 12K
Hz, the equipment are loaded with 1HP loads, and its rotating speed is 1777r/min, i.e., its fundamental frequency is 29.6Hz, according to its component feature frequency
It is 118.1Hz that coefficient, which calculates rolling element fault characteristic frequency,.In bearing operation, bearing ball interacts with inner ring, outer ring,
Its source signal has certain correlation.
Step 1：Single channel observation signal is measured by sensor, its time domain beamformer is as shown in Figure 8.
Selecting frequency slice function, Zai Duigai roads signal carry out frequency slice wavelet transformation, obtained in its 0~fs/2 frequency band
Time-frequency figure, then inverse transformation is carried out to the signal after frequency slice wavelet transformation and obtains reconstruction signal, and draw reconstruction signal
Energy spectrum, as shown in figure 9, its concrete methods of realizing is as follows：
Frequency slice wavelet transformation definition is：In formula：σ is chi
Spend the factor, σ ≠ 0, p ()=exp (- C2/ 2), p* ()=(exp (- C2/2))*, p* (σ (τ-t))=(exp (- (σ (τ-t)
)2/2))*；σ=sqrt (2)/2/0.025.
Frequency slice wavelet transform function W (t, ω, σ) is calculated in time-frequency region (t1,t2,ω1,ω2) component of signal,
That is reconstruction signal：
Wherein, (t1,t2,ω1,ω2) it is (0,1,0,6000)；
For discrete signal, then energy definition isN is signal length, and F (ω) is fa(t) in Fu
Leaf transformation, Fi(ω) is i-th point of F (ω), and n is signal length 8192.With a width of 300Hz, 3 times or so of frequency is characterized,
Step-length is 1Hz.
Step 2：The time-frequency figure and energy spectrum obtained according to step 1, the higher multiple sections of Energy distribution are selected to original
Signal is cut into slices, and obtains the reconstruction signal of multiple sections.Its concrete methods of realizing is as follows：
Referring to Fig. 9, can be seen that from whole energy spectrum, in 300Hz, 1000Hz, 2600Hz, 3200Hz, 3500Hz and
4200Hz energy rises nearby are quick, and energy is extremely low by contrast for its elsewhere, therefore are included according to the selection of LPF principle
Each Energy distribution upper section is cut into slices, and selects [0 500], [0 1200], [0 2800], [0 3400], [0 respectively
3700] and [0 4400] are cut into slices, and obtain its reconstruction signal, are combined into new multidimensional signal xFSWT(t)=(x1,c1,
c2..., c6)T。
Step 3：The reconstruction signal x that step 2 is obtainedFSWT(t)=(x1,c1,c2..., c6)T, using based on principal component point
The source number estimation method of (PCA) is analysed, its characteristic value and characteristic vector is calculated, determines number of source m.Characteristic value is as shown in table 3,
Number of source is 4.Its concrete methods of realizing is as follows：
Step 3.1：Seek the covariance matrix A=(A of reconstruction signalij)6x6, wherein N is data length 8192.
Step 3.2：Calculate above-mentioned covariance matrix A eigenvalue λ1≥λ2≥...λ6＞ 0 and its orthogonal unitization spy
Sign vectorPass through formulaCalculate variance contribution
Rate, when its value is more than 0.90, then it is assumed that covered source signal full detail, calculated its variance contribution ratio, first four it is cumulative and
Reach 98.5%, so that it is determined that number of source is 4.
Step 4：The number of source m=4 obtained according to step 3, then the Matrix Multiplication of preceding 4 characteristic vectors composition is calculated with weight
Structure signal xFSWT(t)=(x1,c1,c2..., c6)T, 4 dimension matrixes after dimensionality reduction are obtained, it is suitable fixed so as to which underdetermined problem be converted into
Problem.
The multidimensional signal x of table 3FSWT(t) characteristic value
Step 5：The new signal battle array being made up of to step 4 vector, i.e., 4 dimension matrixes, using independent component analysis (ICA), is obtained
To 4 separation signals, and its envelope spectrum is sought respectively, as shown in Figure 10；Envelope spectrum is normalized, superposition obtains equivalent
Envelope spectrum, as shown in figure 11, i.e.,Envelope spectrum is observed, fault signature is extracted, realizes single-pass
Road Blind Signal Separation.
As can be seen from Figure 11 it can be seen that occurring larger peak value, wherein 30Hz and fundamental frequency at 30Hz, 59Hz, 118Hz
Relatively, 59Hz is approximately two frequencys multiplication of fundamental frequency to 29.6Hz, 118Hz and rolling element fault characteristic frequency be 118.1Hz extremely
It is close, thus judge it is that bearing roller breaks down, after being changed to ball, analyzed again, find 118Hz points peak
Value disappears, and illustrates the accuracy of analysis, also illustrate that this method to there is the source signal of certain correlation to realize failure spy
Sign extraction.
Claims (2)
1. a kind of feature extracting method based on single channel signal blind separation rolling bearing, it is characterised in that specifically include following
Step：
Step 1：Selecting frequency section wavelet transform function, frequency slice wavelet transformation is carried out to given original vibration signal,
Time-frequency figure is obtained, then inverse transformation is carried out to the signal after frequency slice wavelet transformation and obtains reconstruction signal, and draws reconstruction signal
Energy spectrum；
Step 2：The time-frequency figure and energy spectrum obtained according to step 1, select the multiple sections pair for including energy peak in energy spectrum
Original vibration signal is cut into slices, and obtains the reconstruction signal of multiple sections；
Step 3：The reconstruction signal that step 2 is obtained carries out noise reduction and de-redundancy, it is determined that optimal number of source m and its corresponding m
Individual characteristic vector；
Step 4：The reconstruction signal that the Matrix Multiplication of m characteristic vector composition is obtained with step 2, obtain the dimension matrixes of the m after dimensionality reduction；
Step 5：The m dimension matrixes that step 4 obtains, m separation signal is obtained using independent component analysis, and seek its envelope respectively
Spectrum；Envelope spectrum is normalized, superposition obtains equivalent envelope spectrum；Envelope spectrum is observed, extracts fault signature；
The concrete methods of realizing of the step 1 is as follows：
Frequency slice wavelet transformation definition is：
<mrow>
<mi>W</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>,</mo>
<mi>&omega;</mi>
<mo>,</mo>
<mi>&sigma;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msup>
<mi>&sigma;e</mi>
<mrow>
<mi>i</mi>
<mi>&omega;</mi>
<mi>t</mi>
</mrow>
</msup>
<msubsup>
<mo>&Integral;</mo>
<mrow>
<mo>-</mo>
<mi>&infin;</mi>
</mrow>
<mrow>
<mo>+</mo>
<mi>&infin;</mi>
</mrow>
</msubsup>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>e</mi>
<mrow>
<mi>i</mi>
<mi>&omega;</mi>
<mi>&tau;</mi>
</mrow>
</msup>
<msup>
<mi>p</mi>
<mo>*</mo>
</msup>
<mo>&lsqb;</mo>
<mi>&sigma;</mi>
<mrow>
<mo>(</mo>
<mi>&tau;</mi>
<mo>-</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mi>d</mi>
<mi>&tau;</mi>
<mo>,</mo>
</mrow>
In formula, σ is scale factor, and σ ≠ 0, σ are constant or ω and t function,Or σ=kt, wherein, k is regulation time domain
Or frequency domain sensitivity, t are time domain variable, ω is frequency domain variable；P () is frequency slice function, p*() is p () conjugation
Function；F (τ) is original vibration signal, f (τ) ∈ L2(R)；
Wherein,C is the variable of time domain vibration signal,σ=sqrt (2)/2/0.025；
Inverse transformation is carried out to frequency slice wavelet transform signal, obtains frequency slice wavelet transform signal W (t, ω, σ) in time-frequency
Region (t1,t2,ω1,ω2) component of signal, i.e. reconstruction signal：
<mrow>
<msub>
<mi>f</mi>
<mi>a</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<mi>&pi;</mi>
</mrow>
</mfrac>
<msubsup>
<mo>&Integral;</mo>
<mrow>
<mi>&omega;</mi>
<mn>1</mn>
</mrow>
<mrow>
<mi>&omega;</mi>
<mn>2</mn>
</mrow>
</msubsup>
<msubsup>
<mo>&Integral;</mo>
<mrow>
<mi>t</mi>
<mn>1</mn>
</mrow>
<mrow>
<mi>t</mi>
<mn>2</mn>
</mrow>
</msubsup>
<mi>W</mi>
<mrow>
<mo>(</mo>
<mi>&tau;</mi>
<mo>,</mo>
<mi>&omega;</mi>
<mo>,</mo>
<mi>&sigma;</mi>
<mo>)</mo>
</mrow>
<msup>
<mi>e</mi>
<mrow>
<mi>i</mi>
<mi>&omega;</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>-</mo>
<mi>&tau;</mi>
<mo>)</mo>
</mrow>
</mrow>
</msup>
<mi>d</mi>
<mi>&tau;</mi>
<mi>d</mi>
<mi>&omega;</mi>
</mrow>
Energy definition formula is：In formula, E (ω) EDF, E (ω)=| F (ω) |2,
Wherein F (ω) is reconstruction signal fa(t) Fourier transformation.
2. the feature extracting method as claimed in claim 1 based on single channel signal blind separation rolling bearing, it is characterised in that
The concrete methods of realizing of the step 3 is as follows：
Step 3.1：Seek the covariance matrix A=(A of reconstruction signalij)pxp, wherein
<mrow>
<msub>
<mi>A</mi>
<mrow>
<mi>i</mi>
<mi>j</mi>
</mrow>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mi>n</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mfrac>
<mo>&times;</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</munderover>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mrow>
<mi>i</mi>
<mi>k</mi>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mover>
<mi>x</mi>
<mo>&OverBar;</mo>
</mover>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mrow>
<mi>j</mi>
<mi>k</mi>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mover>
<mi>x</mi>
<mo>&OverBar;</mo>
</mover>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2...</mn>
<mi>p</mi>
</mrow>
1
Wherein, AijFor the element in covariance matrix A, xikFor xFSWT(t) the i-th row k column elements,For xFSWT(t) the i-th row is averaged
Value, n is data length；xjkFor xFSWT(t) jth row k column elements,For xFSWT(t) jth row average value, xFSWT(t) it is reconstruct letter
Number；
Step 3.2：Calculate above-mentioned covariance matrix A eigenvalue λ1≥λ2≥λ3≥...λn＞ 0 and its orthogonal unitization feature
Vector from left to right arrangesAccording to formulaMeter
Variance contribution ratio is calculated, when its value is more than setting value, determines number of source m and its corresponding m characteristic vector.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510363665.3A CN105021399B (en) | 2015-06-26 | 2015-06-26 | A kind of feature extracting method based on single channel signal blind separation rolling bearing |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510363665.3A CN105021399B (en) | 2015-06-26 | 2015-06-26 | A kind of feature extracting method based on single channel signal blind separation rolling bearing |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105021399A CN105021399A (en) | 2015-11-04 |
| CN105021399B true CN105021399B (en) | 2017-12-05 |
Family
ID=54411535
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510363665.3A Expired - Fee Related CN105021399B (en) | 2015-06-26 | 2015-06-26 | A kind of feature extracting method based on single channel signal blind separation rolling bearing |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105021399B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109580224A (en) * | 2018-12-28 | 2019-04-05 | 北京中科东韧科技有限责任公司 | Rolling bearing fault method of real-time |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105928666B (en) * | 2016-04-20 | 2018-12-21 | 中国石油大学(华东) | Leakage acoustic characteristic extracting method based on Hilbert-Huang transform and blind source separating |
| CN106769052B (en) * | 2017-03-21 | 2018-12-21 | 桂林电子科技大学 | A kind of mechanical system rolling bearing intelligent failure diagnosis method based on clustering |
| JP6848813B2 (en) * | 2017-10-25 | 2021-03-24 | 日本製鉄株式会社 | Information processing equipment, information processing methods and programs |
| CN108647667B (en) * | 2017-09-19 | 2019-09-27 | 长安大学 | A kind of implementation method of the frequency domain amplitude spectrum kurtosis figure based on signal Time-frequency Decomposition |
| CN109684898A (en) * | 2017-10-18 | 2019-04-26 | 中国航发商用航空发动机有限责任公司 | Aero-engine and its vibration signal blind separating method and device |
| CN108388232B (en) * | 2018-03-20 | 2020-07-24 | 江南大学 | Method for monitoring operation mode fault in crude oil desalting process |
| CN111503010B (en) * | 2019-09-18 | 2022-03-18 | 山东建筑大学 | Method for extracting broadband radiation noise signal characteristics of centrifugal pump |
| CN113074935B (en) * | 2021-04-01 | 2022-09-13 | 西华大学 | Acoustic separation and diagnosis method for impact fault characteristics of gearbox |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103439110A (en) * | 2013-07-31 | 2013-12-11 | 哈尔滨工程大学 | Early-stage weak fault diagnosis method of rolling bearing |
| CN103499445A (en) * | 2013-09-28 | 2014-01-08 | 长安大学 | Time-frequency slice analysis-based rolling bearing fault diagnosis method |
| CN104359674A (en) * | 2014-10-20 | 2015-02-18 | 广东电网有限责任公司电力科学研究院 | High-speed rolling bearing fault diagnosing method based on time domain and frequency domain state monitoring |
| CN104655423A (en) * | 2013-11-19 | 2015-05-27 | 北京交通大学 | Rolling bearing fault diagnosis method based on time-frequency domain multidimensional vibration feature fusion |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2007304031A (en) * | 2006-05-15 | 2007-11-22 | Nsk Ltd | Abnormality diagnostic apparatus |
-
2015
- 2015-06-26 CN CN201510363665.3A patent/CN105021399B/en not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103439110A (en) * | 2013-07-31 | 2013-12-11 | 哈尔滨工程大学 | Early-stage weak fault diagnosis method of rolling bearing |
| CN103499445A (en) * | 2013-09-28 | 2014-01-08 | 长安大学 | Time-frequency slice analysis-based rolling bearing fault diagnosis method |
| CN104655423A (en) * | 2013-11-19 | 2015-05-27 | 北京交通大学 | Rolling bearing fault diagnosis method based on time-frequency domain multidimensional vibration feature fusion |
| CN104359674A (en) * | 2014-10-20 | 2015-02-18 | 广东电网有限责任公司电力科学研究院 | High-speed rolling bearing fault diagnosing method based on time domain and frequency domain state monitoring |
Non-Patent Citations (3)
| Title |
|---|
| 一种基于时频峭度谱的滚动轴承损伤诊断方法;段晨东 等;《机械工程学报》;20150615;第51卷(第11期);第78-79页第0-1节 * |
| 基于小波分量奇异值分解的单通道盲分离算法;张纯 等;《电子测量与仪器学报》;20111130;第25卷(第11期);第991-994页第1-3节以及第997页第5节 * |
| 基于时频切片分析的故障诊断方法及应用;段晨东 等;《振动与冲击》;20110515;第30卷(第9期);第1-2页第1节 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109580224A (en) * | 2018-12-28 | 2019-04-05 | 北京中科东韧科技有限责任公司 | Rolling bearing fault method of real-time |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105021399A (en) | 2015-11-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105021399B (en) | A kind of feature extracting method based on single channel signal blind separation rolling bearing | |
| Yang et al. | Structural damage identification via a combination of blind feature extraction and sparse representation classification | |
| Amezquita-Sanchez et al. | Signal processing techniques for vibration-based health monitoring of smart structures | |
| CN105699080B (en) | A kind of Wind turbines bearing fault characteristics extracting method based on vibration data | |
| He et al. | Wayside acoustic diagnosis of defective train bearings based on signal resampling and information enhancement | |
| Wang et al. | An adaptive SK technique and its application for fault detection of rolling element bearings | |
| Guo et al. | Envelope extraction based dimension reduction for independent component analysis in fault diagnosis of rolling element bearing | |
| Kang et al. | Singular value decomposition based feature extraction approaches for classifying faults of induction motors | |
| CN105258947B (en) | A kind of Fault Diagnosis of Roller Bearings under operating mode disturbed conditions based on compressed sensing | |
| Amini et al. | Underdetermined blind modal identification of structures by earthquake and ambient vibration measurements via sparse component analysis | |
| Yang et al. | An improved EMD method for modal identification and a combined static-dynamic method for damage detection | |
| CN105862935B (en) | A kind of damnification recognition method for retaining wall structure system | |
| CN103499445A (en) | Time-frequency slice analysis-based rolling bearing fault diagnosis method | |
| Yao et al. | Blind modal identification using limited sensors through modified sparse component analysis by time‐frequency method | |
| Poozesh et al. | Modal parameter estimation from optically-measured data using a hybrid output-only system identification method | |
| CN107525674A (en) | Frequency method of estimation and detection means are turned based on crestal line probability distribution and localised waving | |
| KR20090078075A (en) | Fault diagnosis of inductirn motors by dft and wavelet | |
| Yan et al. | Frequency slice algorithm for modal signal separation and damping identification | |
| Huang et al. | A Fault Diagnosis Approach for Rolling Bearing Based on Wavelet Packet Decomposition and GMM-HMM. | |
| CN102128788A (en) | Improved natural excitation technology-based steel framework damage diagnosis method | |
| CN103473559A (en) | SAR image change detection method based on NSCT domain synthetic kernels | |
| Almasri et al. | Toward compressed sensing of structural monitoring data using discrete cosine transform | |
| Cheng et al. | The identification of a dam's modal parameters under random support excitation based on the Hankel matrix joint approximate diagonalization technique | |
| Ma et al. | Analysis and design of modified window shapes for S-transform to improve time–frequency localization | |
| Meng et al. | Adaptive sparse denoising and periodicity weighted spectrum separation for compound bearing fault diagnosis |
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: 20171205 Termination date: 20180626 |