CN107831013A - A kind of Method for Bearing Fault Diagnosis for strengthening cyclic bispectrum using probability principal component analysis - Google Patents
A kind of Method for Bearing Fault Diagnosis for strengthening cyclic bispectrum using probability principal component analysis Download PDFInfo
- Publication number
- CN107831013A CN107831013A CN201710942216.3A CN201710942216A CN107831013A CN 107831013 A CN107831013 A CN 107831013A CN 201710942216 A CN201710942216 A CN 201710942216A CN 107831013 A CN107831013 A CN 107831013A
- Authority
- CN
- China
- Prior art keywords
- frequency
- bispectrum
- principal component
- signal
- cyclic
- 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
-
- 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
Abstract
The invention belongs to mechanical fault diagnosis field, and the invention discloses a kind of Method for Bearing Fault Diagnosis of probability principal component enhancing cyclic bispectrum, primary signal and the potential probability principal component model contacted of principal component vector can be reflected by initially setting up；Secondly the useful component of de-noising, intactly stick signal is carried out to primary signal using probability principal component model, and signal to noise ratio is significantly increased, solves the problems, such as fault-signal easily by noise jamming；Finally make cyclic bispectrum analysis to noise cancellation signal, and bearing fault is diagnosed to be from the contour map of simple subprogram frequency bispectrum.On the one hand the inventive method can effectively improve signal to noise ratio using probability principal component analysis de-noising, suppress noise, highlight the failure-frequency modulation phenomenon near bearing resonance frequency, so as to increase the hexagonal structure definition of cyclic bispectrum formation；On the other hand, bearing fault type can be visually detected from the contour map of simple subprogram frequency bispectrum.
Description
Technical Field
The invention belongs to the technical field of fault diagnosis of mechanical equipment, and particularly relates to a bearing fault diagnosis method for enhancing cyclic bispectrum by utilizing probability principal component analysis.
Background
The rotating machinery generally comprises a supporting part (bearing), a rotating part (shafting) and a transmission part (gear), is an indispensable important power device of modern equipment, has a complex and various structural form, and is widely applied to the fields of modern machining, nuclear power, wind energy, hydroelectric energy, large-scale unit equipment and aviation equipment. Statistically, the bearing is a vulnerable part in the fault of the rotating machine, and the fault accounts for 19% of the fault of the rotating machine. Therefore, it is valuable and realistic to monitor the operating condition of the bearing in real time to detect abnormal faults at an early stage.
However, early failure characteristics of rotating mechanical bearings are weak and are easily masked by various strong background noises. In addition, the fault signal of the bearing is usually a cyclostationary signal with an impact characteristic, having a non-gaussian and wide frequency band. To accurately extract fault features from strong background noise, it is almost impossible to use conventional spectral analysis, and modern signal processing methods that can handle nonlinear, non-gaussian, and non-white additive noise must be employed. One of the problems that the modern signal processing method urgently needs to solve in the field of mechanical fault diagnosis is a reliable bearing fault feature extraction method.
Probability Principal Component Analysis (PPCA) is essentially a factor Analysis hidden variable model, which expresses all variables of a signal by probability distribution, estimates relevant parameters of the model by constructing a probability Principal Component model and using an Expectation Maximization (EM) algorithm, and finally regenerates observation data by the probability Principal Component model, thereby realizing noise elimination and improving the signal-to-noise ratio. Cyclic bispectrum analysis requires a non-linear transformation of the cyclostationary signal and produces a sine wave of finite intensity. However, the signal itself does not contain any additive sinusoidal components of finite strength, and in order to generate a sinusoidal component, a sinusoidal decimation operation needs to be performed on the basis of nonlinear transformation, and the frequency of the generated sinusoidal wave is the cycle frequency. Because the cyclic bispectrum analysis is easily influenced by non-Gaussian noise, the success and the accuracy of fault diagnosis are directly influenced. Therefore, in order to overcome the defect that the cyclic bispectrum is susceptible to non-Gaussian noise, probability principal component analysis is adopted to enhance signals, a bearing fault diagnosis method for enhancing the cyclic bispectrum by using the probability principal component analysis is provided, and no report is found in the research on the aspect.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention provides a bearing fault diagnosis method for enhancing cycle bispectrum by utilizing probability principal component analysis, and the diagnosis method can accurately detect the bearing fault condition.
In order to achieve the purpose, the invention adopts the technical scheme that: a bearing fault diagnosis method for enhancing cyclic bispectrum by utilizing probabilistic principal component analysis is characterized by comprising the following steps of: (1) establishing a probability pivot model capable of reflecting potential relation between an original signal and a pivot vector; (2) denoising an original signal by using a probability principal component model to obtain a high signal-to-noise ratio signal, completely retaining useful components of the signal, greatly improving the signal-to-noise ratio and solving the problem that a fault signal is easily interfered by noise; (3) performing cyclic bispectrum analysis on the noise-removed signal, calculating third-order cumulant, sine extraction operation and two-dimensional FFT (fast Fourier transform) conversion to obtain a contour map of the single cyclic frequency bispectrum, and finally diagnosing a bearing fault by comparing regular hexagon vertex coordinates in the contour map of the single cyclic frequency bispectrum with fault characteristic frequency;
wherein, step (1) includes the following content:
the probability principal component model is based on the factor analysis method hidden variable model and assumes the variance of Gaussian noise
Isotropy, the generative model produced by the superposition of noise on the factor vectors (principal component vectors), whose mathematical form is:
X＝P·u+E (1)
in formula (1), X = { s = 1 ,s 2 ,…,s m }∈R n×m Is an n × m order matrix generated from an original signal s, wherein: n is the embedding dimension (number of original variables), m is the number of samples; p = { P 1 ,p 2 ,…,p k }∈R n×k Is the load matrix (pivot matrix) of order n × k to be solved, k is the pivot number, and k is<n；u＝{u 1 ,u 2 ,…,u m }∈R k×m As a principal component variable, u i (i =1,2, … m) as a component of the pivot variable; e is obedient N (0, sigma) 2 I) Distributed Gaussian noise, σ 2 Is the variance of the noise; thus, X is subject to N (0,PP) T +σ 2 I) (ii) a gaussian distribution of;
after the probability model is determined, the distribution of the principal component variables, the distribution of the original variables under the given principal component variables and the distribution of the principal component variables under the given original variables are respectively calculated, and then the probability principal component model can be obtained; the specific calculation process is as follows:
assuming that the pivot variable u follows a multivariate standard normal distribution, the probability distribution is:
the conditional distribution of the original variables under a given pivot variable is:
the distribution of the original variables is:
in formula (4), C = PP T +σ 2 I is the variance of the original variables of the order of n multiplied by n, and the conditional distribution of the principal component variables under the given original variables is obtained according to the Bayes theorem:
in formula (5), M = P T P+σ 2 I is a reduced dimension k × k order square matrix; as can be seen from the expressions (2) to (5), the parameters P and σ are determined 2 Then, the probability model can be obtained; estimating parameters by a hidden model classical parameter estimation EM algorithm:
in the formulae (6) and (7),the vector is the covariance matrix of the original variables, tr (-) is the trace of the matrix, i.e. the sum of all the elements on the main diagonal is taken, the equations (6) and (7) are iterated until convergence is reached, and the variance sigma of the load vector P and the Gaussian noise can be obtained 2 ；
The step (2) comprises the following steps:
according to a load vector P obtained by the established probability principal component model, decomposing an original signal X into the sum of outer products of P to obtain the principal component model of the original signal, wherein the mathematical expression of the principal component model is as follows:
in equation (8), each principal element is actually a projection of the original signal X in the direction of the load vector corresponding to this principal element component; the corresponding eigenvalues of the pivot vector reflect the original signal X at p i The degree of coverage of the direction; sorting each pivot vector according to the magnitude of its characteristic value, then its correspondent load vector can express the magnitude of original signal X change direction, and selecting front portion containing main information quantityThe (1,2, …, k) load vectors reconstruct data, and finally an original signal matrix with reduced dimension is obtained, so that the purpose of signal denoising is realized;
the probability pivot model carries out denoising on the original signal depending on the embedding dimension (original variable number) n and the pivot number k; k may be calculated from the Cumulative Variance Contribution Rate (CVCR); the single principal component contribution rate calculation formula is as follows:
in the formula (9), λ is an eigenvalue of a covariance matrix S of the original variables; therefore, the total principal component contribution rate calculation formula is:
in general, k isThe threshold value of (2 d) reaches the minimum value of 85% (document: I.T. Jolliffe, principal Component Analysis (2 d edition), springer-Verlag, new York, 2008);
when k is determined, determining the minimum embedding dimension (the number of original variables) by k < n;
the step (3) comprises the following steps:
by calculating the third-order cumulant, sine extraction operation and two-dimensional FFT transformation, a contour map of a single cycle frequency bispectrum is obtained, which specifically comprises the following steps:
calculating a third order cumulant formula in the frequency domain:
in the formula (11), u (f) represents Fourier transform of u (n), E (. Cndot.) represents expectation, and u * Denotes the conjugation, alpha being the cycle introducedFrequency; the calculation of the cyclic bispectrum from equation (11) can be divided into two steps: estimating a bispectrum; sine extraction operation;
the classical dual-spectrum estimation comprises a nonparametric method and a parametric method, and the nonparametric method comprises a direct method and an indirect method; the invention firstly adopts an indirect method to carry out bispectrum estimation, and the specific algorithm is as follows:
let u (1), u (2),.. U (m) be the sampled signal of a set of principal component variable components, and set f s Is the sampling frequency, N 0 Is the total number of frequency samples, Δ f = f s /N 0 Is the frequency sampling interval in the horizontal and vertical directions of the bispectral region;
step 1: divide the signal u (1), u (2),.. U (m) into k segments of length m, i.e. N = k × m, and subtract the mean of the segments;
step 2: let { u (i) (j) J =1,2,.. M } is the ith segment of data, and the third order cumulant for each segment is estimated:
in the formula (12), m 1 ＝max(0,-k,-l),m 2 ＝min(m-1,m-1-k,m-1-l)；
And 3, step 3: calculating the mean value of the third-order cumulants of all the sections as the estimated value of the third-order cumulants of the whole observation data group
And 4, step 4: performing bispectrum estimation:
in the formula (14), L is less than M-1, and omega (k, L) is a two-dimensional window function;
secondly, sine extraction operation is carried out;
as can be seen from the formula defined in formula (11), the cyclic bispectrum defined in the frequency domain is obtained by shifting the frequency of the bispectrum to the left by the cyclic frequency α; according to the Fourier frequency shift theorem, in the time domain, the third-order cyclic accumulation of the signal is multiplied by a factor e jαt Then, this is the sine extraction operation, and the statistic of the obtained cyclic bispectrum is defined as:
in the formula (15), T is a period, and the third-order cumulant is calculated in the time domainThe formula:
in formula (16), first, second and third momentsRespectively as follows:
the key of sine extraction operation is to estimate the cycle frequency alpha in advance, which directly determines the success of the bispectrum analysis and the accuracy of fault frequency diagnosis; when the cyclic bispectrum is actually applied, the cyclic frequency is determined by adopting a linear search method:
step 1: selecting a specific cyclic frequency set { alpha } in T = [0,2 pi), searching by a fixed step size delta alpha, and calculating each cyclic frequency alpha according to equation (16) j Third order cyclic accumulation of e { alpha }
Step 2: summing all the third-order cumulants calculated in the cyclic frequency set { alpha }, to obtain
And 3, step 3: to pairPerforming two-dimensional Fourier transform to obtain a cyclic frequency interval (alpha) 1 ,α 2 ) Cyclic bispectrum of (c):
the step (4) comprises the following steps: from the contour plot of the bispectrum of the single cycle frequency
And comparing the regular hexagonal vertex coordinates in the contour map of the rate bispectrum with the fault characteristic frequency, and finally diagnosing the bearing fault.
By adopting the scheme, the invention has the advantages that: on one hand, the method can effectively improve the signal-to-noise ratio and inhibit noise by utilizing the probability principal component analysis to eliminate the noise, so that the fault frequency modulation phenomenon near the bearing resonance frequency is highlighted, and the definition of a hexagonal structure formed by circulating double spectrums is increased; on the other hand, the bearing fault type is reliably detected from the contour map of the single cycle frequency bispectrum.
The invention is further described below with reference to the accompanying drawings.
Drawings
FIG. 1 is a bearing fault diagnosis flow chart of probability principal component analysis enhanced cycle bispectrum;
FIG. 2 is a waveform and spectrum diagram of an original simulation signal, (a) a waveform diagram, and (b) a spectrum diagram;
FIG. 3 is a diagram of a signal waveform and spectrum after noise cancellation in probabilistic principal component analysis (a), a waveform diagram, and a spectrum diagram (b);
fig. 4 is a dual-spectrum contour line and a three-dimensional graph of a noise-removed signal at α =1000Hz cycle, (a) a contour graph, and (b) a three-dimensional graph;
fig. 5 is a dual-spectrum contour line and a three-dimensional graph of the original signal at α =1000Hz cycle, (a) a contour graph, and (b) a three-dimensional graph;
FIG. 6 is a waveform diagram of an original signal of a bearing inner ring fault, (a) a time domain diagram, and (b) a frequency spectrum diagram;
FIG. 7 is a diagram of a signal waveform and spectrum after noise cancellation in probabilistic principal component analysis, (a) a time domain diagram, and (b) a spectrum diagram;
FIG. 8 is a circular bispectral contour line and three-dimensional map, (a) a contour map, (b) a partial magnified view of the contour line, and (c) a three-dimensional map;
FIG. 9 shows a cyclic bispectral contour line and a three-dimensional plot of an unpoised signal, (a) a contour plot, and (b) a three-dimensional plot.
Detailed Description
The protection scope of the present invention is not limited to the following embodiments, and those skilled in the art can implement the present invention in other embodiments according to the disclosure of the present invention, or make simple changes or modifications according to the design structure and idea of the present invention, and fall into the protection scope of the present invention.
An embodiment of the present invention is shown in FIG. 1=9.
Specific example 1: in order to verify the correctness of the bearing fault diagnosis method for enhancing the cyclic bispectrum by probability principal component analysis, a digital simulation model is established according to the bearing fault characteristics:
in equation (20), n (t) is a zero-mean stationary random signal, h (t) is a series of repetitive pulse bursts resulting from bearing failure, A i For randomly modulating amplitude, T is the impact interval of two adjacent pulses, tau i Showing the minor fluctuation of the ith impact with respect to the fault averaging period T, i being the number of impacts.
The response of the bearing to each pulse is defined as:
h(t)＝e -Bt sin 2πf n t (21)
in the formula (21), B is a coefficient determined by the resonance frequency, mass coefficient and stress relaxation time of the bearing, and f n Is the center frequency of the bearing.
The severity of the bearing failure is affected by random variables such as radial load distribution on each rolling element, dynamic stiffness of assembly, roundness errors of the rolling elements and the ring, and rolling element group difference, so that the amplitude A of the impact pulse i Should be a random modulation amplitude, expressed as:
A i ＝cos(2πf A t+φ A )+C A +randn(t) (22)
in the formula (22) < phi > A And C A Is an arbitrary constant, f A To modulate the frequency (outer ring failure f) A =0; at inner ring failure f A Equal to the shaft frequency; at fault of rolling body f A Equal to the cage transition frequency), gaussian random noise randn (t) models the randomness of the amplitude.
The zero-mean stationary random signal n (t) is modeled by a sinusoidal random phase signal, defined as:
n(t,φ i )＝Acos(ω 0 t+φ i ) φ i ∈(-π,π) (23)
in the formula (23), A, ω 0 Is a constant value of phi i Are random variables distributed uniformly over (-pi, pi).
Setting bearingFault impulse period T =0.02, i.e. fault frequency of 50Hz, modulation frequency f A =1480/60Hz, sampling frequency f s =4000Hz, small fluctuation τ from mean period of failure i =0.0155T, center frequency f of bearing n ＝1000Hz。
As can be seen from equation (23), the simulation signal contains both white gaussian noise, which appears in the amplitude-modulated portion, and stationary non-gaussian noise as additive noise of the entire signal. The embodiment calculates the cyclic bispectrum by analyzing signals before and after denoising through the probability principal component, and verifies the effectiveness of PPCA denoising and the accuracy of cyclic bispectrum diagnosis.
Selecting a sample number m =4076 from the total sampling points as preprocessed original data, and dividing the other samples into 20 segments, namely, the number of original variables n =20; when the number of pivot variables k =2 is set, the original data is represented as X 20×4076 In the form of a matrix. The waveform and spectrum of the original signal before noise cancellation by probabilistic principal component analysis is shown in fig. 2. The waveform and spectrum of the noise-removed signal by probabilistic principal component analysis are shown in fig. 3. Comparing the waveform diagrams of fig. 2 (a) and fig. 3 (a) shows that: the noise elimination of the probability principal component analysis can effectively inhibit the noise. The spectrograms of the signals before and after the noise cancellation are shown in fig. 2 (b) and fig. 3 (b), respectively. Comparing the spectrograms of the signals before and after noise cancellation, it can also be seen that: the P probability principal component analysis can effectively eliminate noise, so that the frequency spectrum characteristics of signals are clearer and more accurate.
Determining the cycle frequency by adopting a linear search method so as to determine the optimal cycle frequency alpha =1000Hz, and then performing sine extraction operation to obtain a contour map and a three-dimensional map calculation result of the bispectrum analysis, as shown in 4. As can be seen from FIG. 4, the vertex coordinates of the six deformations are (0, -50), (50, -50), (50,0), (0,50), (-50,50), (-50,0), respectively. This indicates that: for the cyclostationary signal, the signal-to-noise ratio of the central frequency position is highest, the cyclostationary signal is not easily interfered by noise, the relationship between the central frequency and other frequency spectrum components is relatively fixed, the fault frequency is exactly equal to 50Hz, and the stability is strongest. In order to verify the effect of enhancing the cyclic bispectrum by probability principal component analysis, the cyclic bispectrum is directly calculated for the original simulation signal, and the cyclic frequency is still alpha =1000Hz. The contour diagram and three-dimensional diagram are shown in fig. 5. As can be seen from fig. 5 (a), a regular cyclic bispectral hexagonal structure cannot be constructed from the failure frequency values due to the stationary non-gaussian noise interference.
Specific example 2: bearing inner race fault diagnosis
And taking a fault signal of a bearing inner ring of a certain real transmission system, wherein the outer ring of the bearing for experiment has a fault, the model is ER-12K, the number of the rolling bodies is 8, the diameter of the rolling bodies is 7.9375 mm, the diameter of a pitch circle of the bearing is 33.4772mm, and the theoretical calculation result of the fault characteristic frequency of the inner ring of the rolling bearing is 197.1Hz. The experimental bearing runs in an idle load mode, the rotating speed is 2990r/min, the sampling frequency is 25.6kHz, the number of sampling points is 327680, the number of the sampling points m =327660 is selected as preprocessed original data, and the preprocessed original data is divided into 20 sections, namely the number of original variables n =20; assuming that the number of pivot variables k =2, the original data is represented as X 20×327660 In the form of a matrix.
The original signal waveform and spectrogram are shown in fig. 6, the original data signal is subjected to probability principal component analysis denoising preprocessing, and the denoising signal waveform and spectrogram are shown in fig. 7. As can be seen from fig. 6, the useful information is completely buried in the noise. As is clear from fig. 7, the impact component is present in a prominent manner, and the noise is suppressed to a certain extent, but the failure characteristics cannot be diagnosed.
In order to further extract the fault characteristic information of the bearing inner ring, a cyclic bispectrum is calculated for the signal subjected to probability principal component analysis and noise elimination. Determining the cycle frequency by using a linear search method in order to determine the optimal cycle frequency α =2551Hz, a sine extraction operation is then performed to obtain a contour map and a three-dimensional map calculation result of the bispectrum analysis, as shown in fig. 8. It is clear from the enlarged view of fig. 8 (b) that a regular hexagonal structure is generated around the central coordinate (0,0), and six vertexes are respectively composed of pairs of fault characteristic frequencies, the frequency coordinates are respectively (0, -197), (197, -197), (197,0), (0,197), (-197,197), (-197,0), which indicates that the signal modulation frequency at the bearing central frequency α =2551Hz is 197Hz. The modulation frequency is close to the theoretical calculation result 197.1Hz of the fault characteristic frequency of the inner ring of the rolling bearing, so that the fault of the inner ring of the bearing is diagnosed.
In order to verify the effect of enhancing the cyclic bispectrum by probability principal component analysis, the cyclic bispectrum is directly calculated for the rolling bearing signals of the original data, and the cyclic frequency is still taken as the central frequency of alpha =2551Hz. The contour diagram and three-dimensional diagram are shown in fig. 9. As can be seen from fig. 9 (a), a regular cyclic bispectral hexagonal structure cannot be constructed from the failure frequency values due to the stationary non-gaussian noise interference.
Claims (1)
1. A bearing fault diagnosis method for enhancing cyclic bispectrum by using probabilistic principal component analysis is characterized by comprising the following steps of: (1) establishing a probability pivot model capable of reflecting potential relation between an original signal and a pivot vector; (2) denoising an original signal by using a probability principal component model to obtain a signal with a high signal-to-noise ratio, completely retaining useful components of the signal and greatly improving the signal-to-noise ratio; (3) performing cyclic bispectrum analysis on the noise-eliminating signal, calculating third-order cumulant, sine extraction operation and two-dimensional FFT (fast Fourier transform) conversion to obtain a contour map of the single cyclic frequency bispectrum, and finally diagnosing a bearing fault by comparing regular hexagonal vertex coordinates in the contour map of the single cyclic frequency bispectrum with fault characteristic frequency;
wherein, step (1) includes the following content:
the probability principal component model is a generation model generated by superposing noise by a factor vector (principal component vector) on the assumption of isotropic variance of Gaussian noise on the basis of a factor analysis hidden variable model, and the mathematical form of the probability principal component model is as follows:
X＝P·u+E (1)
in formula (1), X = { s = 1 ,s 2 ,…,s m }∈R n×m Is an n × m order matrix generated from an original signal s, wherein: n is the embedding dimension (number of original variables), m is the number of samples; p = { P 1 ,p 2 ,…,p k }∈R n×k Is the load matrix (pivot matrix) of order n × k to be solved, k is the pivot number, and k is<n；u＝{u 1 ,u 2 ,…,u m }∈R k×m Is mainly composed ofElement variable, u i (i =1,2, … m) as a component of the pivot variable; e is obedient N (0, sigma) 2 I) Distributed Gaussian noise, σ 2 Is the variance of the noise; thus, X is subject to N (0,PP) T +σ 2 I) A gaussian distribution of (d);
after the probability model is determined, the distribution of the principal component variables, the distribution of the original variables under the given principal component variables and the distribution of the principal component variables under the given original variables are respectively calculated, and then the probability principal component model can be obtained; the specific calculation process is as follows:
assuming that the pivot variable u follows a multivariate standard normal distribution, the probability distribution is:
the conditional distribution of the original variables under a given pivot variable is:
the distribution of the original variables is:
in formula (4), C = PP T +σ 2 I is the variance of the original variables of n multiplied by n orders, and the conditional distribution of the principal component variables under the given original variables is obtained according to the Bayes theorem:
in formula (5), M = P T P+σ 2 I is a k multiplied by k order square matrix with reduced dimensions; as can be seen from the expressions (2) to (5), the parameters P and σ are determined 2 Then, a probability model can be obtained; estimating parameters by a hidden model classical parameter estimation EM algorithm:
in the formulae (6) and (7),the vector is the covariance matrix of the original variables, tr (-) is the trace of the matrix, i.e. the sum of all the elements on the main diagonal is taken, the equations (6) and (7) are iterated until convergence is reached, and the variance sigma of the load vector P and the Gaussian noise can be obtained 2 ；
The step (2) comprises the following steps:
according to a load vector P obtained by the established probability principal component model, decomposing an original signal X into the sum of outer products of P to obtain the principal component model of the original signal, wherein the mathematical expression of the principal component model is as follows:
in equation (8), each principal element is actually a projection of the original signal X in the direction of the load vector corresponding to this principal element component; the corresponding eigenvalues of the pivot vector reflect the original signal X at p i The degree of coverage of the direction; sorting each principal component vector according to the magnitude of the characteristic value, then representing the magnitude of the original signal X change direction by the corresponding load vector, selecting the front (1,2, …, k) load vectors containing main information quantity to reconstruct data, finally obtaining a dimension-reduced original signal matrix, and realizing the purpose of signal denoising;
the probability pivot model denoises the original signal depending on the embedding dimension (original variable number) n and the pivot number k; k may be calculated from the Cumulative Variance Contribution Rate (CVCR); the single principal component contribution rate calculation formula is as follows:
in the formula (9), λ is an eigenvalue of a covariance matrix S of the original variables; therefore, the total principal component contribution rate calculation formula is:
in general, k isA minimum value of 85%;
when k is determined, determining the minimum embedding dimension (the number of original variables) by k < n;
the step (3) comprises the following steps:
by calculating the third-order cumulant, sine extraction operation and two-dimensional FFT transformation, a contour map of a single cycle frequency bispectrum is obtained, which specifically comprises the following steps:
calculating a third order cumulant formula in the frequency domain:
in the formula (11), u (f) represents Fourier transform of u (n), E (. Cndot.) represents expectation, and u * Representing the conjugate, alpha is the introduced cycle frequency; the calculation of the cyclic bispectrum from equation (11) can be divided into two steps: estimating a bispectrum; sine extraction operation;
the bispectrum estimation adopts an unparameterization method or a parameterization method, the unparameterization method comprises a direct method and an indirect method, and the indirect method adopts the specific algorithm for bispectrum estimation as follows:
let u (1), u (2),.. U (m) be the sampled signal of a set of principal component variable components, and set f s Is the sampling frequency, N 0 Is the total number of frequency samples, Δ f = f s /N 0 Is the frequency sampling interval in the horizontal and vertical directions of the bispectral region;
step 1: divide the signal u (1), u (2),.. U (m) into k segments of length m, i.e. N = k × m, and subtract the mean of the segments;
step 2: let { u (i) (j) J =1,2,.. M } is the ith segment of data, and the third order cumulant for each segment is estimated:
in the formula (12), m 1 ＝max(0,-k,-l),m 2 ＝min(m-1,m-1-k,m-1-l)；
And 3, step 3: calculating the mean value of the third-order cumulants of all the sections as the third-order cumulants estimated value of the whole observation data group
And 4, step 4: performing bispectrum estimation:
in the formula (14), L is less than M-1, and omega (k, L) is a two-dimensional window function;
secondly, sine extraction operation is carried out;
as can be seen from the formula defined in formula (11), the cyclic bispectrum defined in the frequency domain is obtained by shifting the frequency of the bispectrum to the left by the cyclic frequency α; according to the Fourier frequency shift theorem, in the time domain, the third-order cyclic accumulation of the signal is multiplied by a factor e jαt Then, this is the sine extraction operation, and the statistic of the obtained cyclic bispectrum is defined as:
in the formula (15), T is a period, and the third-order cumulant is calculated in the time domainThe formula is as follows:
in formula (16), first, second and third momentsRespectively as follows:
the key of sine extraction operation is to estimate the cycle frequency alpha in advance, which directly determines the success of the bispectrum analysis and the accuracy of fault frequency diagnosis; when the cyclic bispectrum is actually applied, the cyclic frequency is determined by adopting a linear search method:
step 1: selecting a specific cyclic frequency set { alpha } in T = [0,2 pi), searching by a fixed step size delta alpha, and calculating each cyclic frequency alpha according to equation (16) j Third order cyclic accumulation of e { alpha }
Step 2: summing all the third-order cumulants calculated in the cyclic frequency set { alpha }, to obtain
And 3, step 3: to pairPerforming two-dimensional Fourier transform to obtain a cyclic frequency interval (alpha) 1 ,α 2 ) Cyclic bispectrum of (c):
the step (4) comprises the following steps: and comparing the regular hexagonal vertex coordinates in the contour map of the single cycle frequency bispectrum with the fault characteristic frequency according to the contour map of the single cycle frequency bispectrum, and finally diagnosing the bearing fault.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710942216.3A CN107831013B (en) | 2017-10-11 | 2017-10-11 | Utilize the Method for Bearing Fault Diagnosis of probability principal component analysis enhancing cyclic bispectrum |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710942216.3A CN107831013B (en) | 2017-10-11 | 2017-10-11 | Utilize the Method for Bearing Fault Diagnosis of probability principal component analysis enhancing cyclic bispectrum |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107831013A true CN107831013A (en) | 2018-03-23 |
| CN107831013B CN107831013B (en) | 2019-11-22 |
Family
ID=61647836
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710942216.3A Active CN107831013B (en) | 2017-10-11 | 2017-10-11 | Utilize the Method for Bearing Fault Diagnosis of probability principal component analysis enhancing cyclic bispectrum |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107831013B (en) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108548957A (en) * | 2018-05-23 | 2018-09-18 | 西北工业大学 | The double-spectrum analysis method being combined based on circular modulating frequency spectrum and segmentation cross-correlation |
| CN108647707A (en) * | 2018-04-25 | 2018-10-12 | 北京旋极信息技术股份有限公司 | Probabilistic neural network creation method, method for diagnosing faults and device, storage medium |
| CN109029999A (en) * | 2018-09-19 | 2018-12-18 | 河北工业大学 | Fault Diagnosis of Roller Bearings based on enhanced modulation double-spectrum analysis |
| CN109374293A (en) * | 2018-10-29 | 2019-02-22 | 桂林电子科技大学 | A kind of gear failure diagnosing method |
| CN111207926A (en) * | 2019-12-27 | 2020-05-29 | 三明学院 | Fault diagnosis method based on rolling bearing, electronic device and storage medium |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1926413A (en) * | 2004-10-18 | 2007-03-07 | 日本精工株式会社 | Abnormality diagnosis system for machinery |
| CN102519726A (en) * | 2011-12-28 | 2012-06-27 | 昆明理工大学 | Acoustic-based diagnosis (ABD) method for compound fault of rolling bearing |
| US20140365411A1 (en) * | 2013-06-05 | 2014-12-11 | The Trustees Of Columbia University In The City Of New York | Monitoring Health of Dynamic System Using Speaker Recognition Techniques |
| CN105699080A (en) * | 2015-12-18 | 2016-06-22 | 华北电力大学（保定） | Wind turbine generator set bearing fault feature extraction method based on vibration data |
| CN106338385A (en) * | 2016-08-25 | 2017-01-18 | 东南大学 | Rotation machinery fault diagnosis method based on singular spectrum decomposition |
| CN107144428A (en) * | 2017-03-17 | 2017-09-08 | 北京交通大学 | A kind of rail traffic vehicles bearing residual life Forecasting Methodology based on fault diagnosis |
-
2017
- 2017-10-11 CN CN201710942216.3A patent/CN107831013B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1926413A (en) * | 2004-10-18 | 2007-03-07 | 日本精工株式会社 | Abnormality diagnosis system for machinery |
| CN102519726A (en) * | 2011-12-28 | 2012-06-27 | 昆明理工大学 | Acoustic-based diagnosis (ABD) method for compound fault of rolling bearing |
| US20140365411A1 (en) * | 2013-06-05 | 2014-12-11 | The Trustees Of Columbia University In The City Of New York | Monitoring Health of Dynamic System Using Speaker Recognition Techniques |
| CN105699080A (en) * | 2015-12-18 | 2016-06-22 | 华北电力大学（保定） | Wind turbine generator set bearing fault feature extraction method based on vibration data |
| CN106338385A (en) * | 2016-08-25 | 2017-01-18 | 东南大学 | Rotation machinery fault diagnosis method based on singular spectrum decomposition |
| CN107144428A (en) * | 2017-03-17 | 2017-09-08 | 北京交通大学 | A kind of rail traffic vehicles bearing residual life Forecasting Methodology based on fault diagnosis |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108647707A (en) * | 2018-04-25 | 2018-10-12 | 北京旋极信息技术股份有限公司 | Probabilistic neural network creation method, method for diagnosing faults and device, storage medium |
| CN108647707B (en) * | 2018-04-25 | 2022-09-09 | 北京旋极信息技术股份有限公司 | Probabilistic neural network creation method, failure diagnosis method and apparatus, and storage medium |
| CN108548957A (en) * | 2018-05-23 | 2018-09-18 | 西北工业大学 | The double-spectrum analysis method being combined based on circular modulating frequency spectrum and segmentation cross-correlation |
| CN108548957B (en) * | 2018-05-23 | 2020-08-07 | 西北工业大学 | Dual-spectrum analysis method based on combination of cyclic modulation spectrum and piecewise cross correlation |
| CN109029999A (en) * | 2018-09-19 | 2018-12-18 | 河北工业大学 | Fault Diagnosis of Roller Bearings based on enhanced modulation double-spectrum analysis |
| WO2020057066A1 (en) * | 2018-09-19 | 2020-03-26 | 河北工业大学 | Fault diagnosis method employing enhanced modulation signal bispectrum analysis for rolling bearing |
| CN109374293A (en) * | 2018-10-29 | 2019-02-22 | 桂林电子科技大学 | A kind of gear failure diagnosing method |
| CN109374293B (en) * | 2018-10-29 | 2020-07-24 | 珠海市华星装备信息科技有限公司 | Gear fault diagnosis method |
| CN111207926A (en) * | 2019-12-27 | 2020-05-29 | 三明学院 | Fault diagnosis method based on rolling bearing, electronic device and storage medium |
| CN111207926B (en) * | 2019-12-27 | 2022-02-01 | 三明学院 | Fault diagnosis method based on rolling bearing, electronic device and storage medium |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107831013B (en) | 2019-11-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107831013A (en) | A kind of Method for Bearing Fault Diagnosis for strengthening cyclic bispectrum using probability principal component analysis | |
| Wang et al. | Rolling element bearing fault diagnosis via fault characteristic order (FCO) analysis | |
| Hu et al. | An adaptive spectral kurtosis method and its application to fault detection of rolling element bearings | |
| CN110186682B (en) | Rolling bearing fault diagnosis method based on fractional order variation modal decomposition | |
| Du et al. | Feature identification with compressive measurements for machine fault diagnosis | |
| Liu et al. | Application of correlation matching for automatic bearing fault diagnosis | |
| Chen et al. | Wavelet denoising for the vibration signals of wind turbines based on variational mode decomposition and multiscale permutation entropy | |
| CN110926594B (en) | Method for extracting time-varying frequency characteristics of rotary machine signal | |
| CN109029999B (en) | Rolling bearing fault diagnosis method based on enhanced modulation bispectrum analysis | |
| Ding et al. | Sparsity-based algorithm for condition assessment of rotating machinery using internal encoder data | |
| Chen et al. | A fault pulse extraction and feature enhancement method for bearing fault diagnosis | |
| Ma et al. | An improved time-frequency analysis method for instantaneous frequency estimation of rolling bearing | |
| Xiang et al. | Comparison of Methods for Different Time-frequency Analysis of Vibration Signal. | |
| CN112668518A (en) | VMSST time-frequency analysis method for vibration fault signal | |
| Xiao et al. | Research on fault feature extraction method of rolling bearing based on NMD and wavelet threshold denoising | |
| Zheng et al. | Faults diagnosis of rolling bearings based on shift invariant K-singular value decomposition with sensitive atom nonlocal means enhancement | |
| Zhang et al. | Research on the fault diagnosis method for rolling bearings based on improved VMD and automatic IMF acquisition | |
| Liu et al. | Rotating machinery fault diagnosis under time-varying speeds: A review | |
| Yan et al. | A bearing fault feature extraction method based on optimized singular spectrum decomposition and linear predictor | |
| Zhang et al. | Generalized transmissibility damage indicator with application to wind turbine component condition monitoring | |
| Gong et al. | Application of optimized multiscale mathematical morphology for bearing fault diagnosis | |
| Huang et al. | A novel wheelset bearing fault diagnosis method integrated CEEMDAN, periodic segment matrix, and SVD | |
| Xu et al. | Rolling bearing fault feature extraction via improved SSD and a singular-value energy autocorrelation coefficient spectrum | |
| Wang et al. | Multi-domain extreme learning machine for bearing failure detection based on variational modal decomposition and approximate cyclic correntropy | |
| Liu et al. | A morphology filter-assisted extreme-point symmetric mode decomposition (MF-ESMD) denoising method for bridge dynamic deflection based on ground-based microwave interferometry |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| EE01 | Entry into force of recordation of patent licensing contract | ||
| EE01 | Entry into force of recordation of patent licensing contract |
Application publication date: 20180323 Assignee: ZHEJIANG GREATWALL MIXERS CO.,LTD. Assignor: Wenzhou University Contract record no.: X2023330000104 Denomination of invention: Bearing Fault Diagnosis Method Using Probabilistic Principal Component Analysis to Enhance Cyclic Bispectrum Granted publication date: 20191122 License type: Common License Record date: 20230311 |