CN111089726B - Rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition - Google Patents

Rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition Download PDF

Info

Publication number
CN111089726B
CN111089726B CN202010046098.XA CN202010046098A CN111089726B CN 111089726 B CN111089726 B CN 111089726B CN 202010046098 A CN202010046098 A CN 202010046098A CN 111089726 B CN111089726 B CN 111089726B
Authority
CN
China
Prior art keywords
singular
spectrum
decomposition
rolling bearing
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.)
Active
Application number
CN202010046098.XA
Other languages
Chinese (zh)
Other versions
CN111089726A (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.)
Southeast University
Original Assignee
Southeast 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 Southeast University filed Critical Southeast University
Priority to CN202010046098.XA priority Critical patent/CN111089726B/en
Publication of CN111089726A publication Critical patent/CN111089726A/en
Application granted granted Critical
Publication of CN111089726B publication Critical patent/CN111089726B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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

Abstract

The invention discloses a rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition, which comprises the following steps: step 1, mounting an acceleration sensor near a rolling bearing seat and collecting vibration signals; step 2, decomposing the acquired vibration signals by using an optimal dimension singular spectrum decomposition algorithm to obtain a plurality of singular spectrum components with relatively definite physical significance; step 3, selecting singular spectrum components containing rich fault characteristic information as principal component components according to a kurtosis criterion; step 4, calculating a 1.5-dimensional frequency domain weighted energy spectrum of the principal component; and 5, observing whether an obvious peak value appears at the fault characteristic frequency from the 1.5-dimensional frequency domain weighted energy spectrum, thereby realizing accurate diagnosis of the rolling bearing fault. The method is simple and easy to implement, and can more accurately realize fault diagnosis of the rolling bearing compared with other prior art.

Description

Rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition
Technical Field
The invention relates to the technical field of fault diagnosis in mechanical equipment, in particular to a rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition.
Background
Rolling bearings are key components of most rotating machines, and are widely used in various machines due to large load bearing capacity and remarkable energy saving effect. However, the rolling bearing is often in severe working conditions such as heavy load, high speed and the like during use, and various faults inevitably occur, so that the normal operation of the whole mechanical equipment is influenced. Therefore, the fault of the rolling bearing can be found in time, and the method has important significance for ensuring the safe and stable operation of mechanical equipment. In addition, due to transient characteristic changes of a mechanical system during operation, non-stationary characteristics (such as time-varying frequency characteristics) appear in the vibration response of the system, so that the conventional fast fourier method cannot be used for extracting the transient frequency characteristics sometimes. Therefore, it is necessary to find an effective fault feature extraction and rolling bearing fault detection method.
Mechanical systems are often subjected to multiple complex excitations simultaneously, so that the vibration signal contains many complex sub-signal components. Some representative time-frequency analysis techniques, such as Wavelet Transform (WT) and Empirical Mode Decomposition (EMD), can separate potential components with different frequency bands from the raw vibration data. However, these methods have some inherent drawbacks for non-stationary signal processing. For example, the WT has a higher spectral resolution at low frequencies, but its decomposition results depend on the choice of wavelet basis functions and may result in energy leakage associated with faulty signals. The decomposition process of EMD is limited by the local characteristics of the signal, so that some defects, such as end effects, modal aliasing, over-envelope and under-envelope, can occur.
In recent years, foreign researchers have proposed an adaptive signal Decomposition method based on Singular Spectrum Analysis (SSA), namely Singular Spectrum Decomposition (SSD), which adaptively determines the embedding dimension of each SSA process by an empirical criterion and reconstructs a signal from high frequency to low frequency. However, the embedding dimension in the singular spectral decomposition is chosen according to empirical criteria, which may lead to situations of modal aliasing and over-decomposition. Furthermore, some in-band noise and interference is more or less present in the singular spectral components obtained by the singular spectral decomposition. Therefore, there is a need for a technique to eliminate in-band noise and interference components while enhancing the strength of the fault.
Disclosure of Invention
The invention designs a novel rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition for realizing self-adaptive decomposition and fault characteristic enhancement of vibration signals of rolling bearings and other components, and aims to realize self-adaptive selection of embedded dimension in each iteration in a singular spectrum decomposition algorithm, improve the decomposition efficiency of the singular spectrum decomposition algorithm, improve the processing precision of the singular spectrum decomposition algorithm on signals, inhibit noise and harmonic interference in singular spectrum components through the idea of 1.5-dimensional spectrum and frequency domain weighted energy operator cooperation and enhance fault characteristic frequency.
The invention provides a rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition, which comprises the following steps of:
step 1: an acceleration sensor is arranged near the rolling bearing to collect vibration signals;
step 2: decomposing the acquired vibration signals by using an optimal dimension singular spectrum decomposition algorithm to obtain a plurality of singular spectrum components with relatively definite physical significance;
step 3, selecting singular spectrum components containing rich fault characteristic information as principal component components according to a kurtosis criterion;
step 4, calculating a 1.5-dimensional frequency domain weighted energy spectrum of the principal component;
and 5, observing whether an obvious peak value appears at the fault characteristic frequency from the 1.5-dimensional frequency domain weighted energy spectrum, thereby realizing accurate diagnosis of the rolling bearing fault.
The specific steps of decomposing the vibration signal by using the optimal embedding dimension singular spectrum decomposition in the step 2 are as follows:
step 2.1, adaptively determining the input sequence v of the iterative decomposition(j)(n) an optimal embedding dimension M when constructing a trajectory matrix, wherein j is an iteration number;
step 2.2, constructing a new track matrix expression for the vibration signals: taking the sequence X (n) {1,2,3,4,5} as an example, if the embedding dimension M is set to 3, the corresponding matrix X is
Figure GDA0003308682190000021
In the formula, the left half of the vertical line in the matrix X corresponds to a new trajectory matrix in the singular spectral decomposition;
inputting the iteration into a sequence v according to the track matrix construction mode(j)(N) constructing an M multiplied by N track matrix X, and performing singular value decomposition on the track matrix X to obtain:
Figure GDA0003308682190000022
wherein U is E.RP×PAnd u isk∈RP×1Is the kth column vector of U, V ∈ RN×NAnd v isk∈RN×1Is the kth column vector of V, Λ ∈ RP×NRepresenting a matrix of singular values, deltakIts k-th singular value;
step 2.3, in the first iteration, if a large trend item is detected, only the first singular value of the singular value matrix and the left and right eigenvectors thereof are used for obtaining g(1)(n) that is
Figure GDA0003308682190000031
And g is(1)(n) from X1Is obtained from the diagonal average of (1); when the number of iterations j>1 hour, component sequence g(j)(n) must have a physically well-defined time scale, whose frequency components are mainly concentrated in the frequency band [ f ]max-Δf,fmax+Δf]To (c) to (d); Δ f denotes the bandwidth and is determined by the residual component vj(n) the power spectral density function is obtained by gaussian interpolation; creating a subset I based on all feature groups with prominent dominant frequencies in the frequency band range for the left feature vector and one feature group with the largest contribution to the dominant peak energy of the selected modal componentsj＝{i1,…,ipThen through the matrix
Figure GDA0003308682190000032
Reconstructing the corresponding component sequence g by diagonal averaging(j)(n)；
Step 2.4, the component sequence g obtained in step 2.3 is subjected to(j)(n) separating from the original signal to obtain a residual component vj+1(n); calculating the normalized mean square error NMSE of the residual component and the original sequence, if NMSE is less than a given threshold th, stopping iteration, otherwise repeating steps 2.1-2.3, the normalized mean square error is defined as:
Figure GDA0003308682190000033
under adaptive determination in step 2.1Input sequence v at one iteration(j)(n) the specific steps of the optimal embedding dimension are:
step 2.1.1: at the j-th iteration, according to the formula
Figure GDA0003308682190000034
A range of embedding dimensions M to be selected is determined,
Figure GDA0003308682190000035
for the current input sequence v(j)(n) the frequency value corresponding to the largest peak in the power spectral density function,
Figure GDA0003308682190000036
representing rounded numbers;
step 2.1.2: taking each value in the value interval of the embedding dimension M as the embedding dimension, and respectively carrying out comparison on the current input sequence v(j)(n) constructing a trajectory matrix and performing decomposition and reconstruction, i.e. performing steps 2.2-2.3, obtaining a series of singular spectral components
Figure GDA0003308682190000037
l is the length of the embedding dimension interval;
step 2.1.3: an Over Decomposition Index (ODI) of the obtained singular spectral components is calculated, the ODI being defined as:
Figure GDA0003308682190000038
in the formula (I), the compound is shown in the specification,
Figure GDA0003308682190000039
representation of SSCjIn
Figure GDA00033086821900000310
The amplitude of (a) of (b) is,
Figure GDA00033086821900000311
representing the input signal v of the current iteration(j)(n) in
Figure GDA00033086821900000312
The amplitude of (c).
Step 2.1.4: using frequency band division method in empirical wavelet decomposition to input signal v of this iteration(j)The spectrum of (n) is adaptively band divided. After the frequency band division is finished, inputting a signal v according to the iteration(j)(n) main frequency
Figure GDA0003308682190000041
The frequency band of the channel is the singular spectral component SSC obtained in the iterationjClassifying into corresponding frequency bands, and then calculating a modal aliasing Index (MMI) of the obtained singular spectral component, wherein the MMI is defined as:
Figure GDA0003308682190000042
where I denotes the singular spectral component SSCjThe total number of spectral lines in the associated frequency band, T representing the singular spectral component SSCjTotal number of lines over the whole spectrum, AflAnd AftAre respectively the frequency flAnd ftThe corresponding amplitude.
Step 2.1.5: will be provided with
Figure GDA0003308682190000043
Respective over-decomposition index ODI is multiplied by the modal aliasing index MMI to obtain the corresponding objective function value
Figure GDA0003308682190000044
Namely;
Figure GDA0003308682190000045
step 2.1.6: and considering the corresponding embedding dimension when the objective function value is maximum as the optimal embedding dimension.
In step 3, when the singular spectrum components including rich fault feature information are selected as principal component components according to a kurtosis criterion, the kurtosis value of each singular spectrum component is calculated, and the singular spectrum component with the highest kurtosis is selected as the principal component, where the kurtosis is defined as:
Figure GDA0003308682190000046
wherein mu is the mean value of the sequence x (n), and sigma is the standard deviation of the sequence x (n).
The step 4 of calculating the 1.5-dimensional frequency domain weighted energy spectrum of the principal component comprises the following specific steps:
step 4.1, carrying out frequency domain weighted energy operator operation on the principal component x (t) to obtain an instantaneous energy signal thereof
Figure GDA0003308682190000047
The frequency domain weighted energy operator is defined as:
Figure GDA0003308682190000048
where H [. cndot. ] represents the Hilbert transform, and Γ [ x (t) ], represents the frequency domain weighted energy operator of x (t).
Step 4.2, calculating instantaneous energy signal
Figure GDA0003308682190000049
The third order cumulant diagonal slice of (a) can be written as:
Figure GDA00033086821900000410
where E {. cndot } represents the mathematical expectation and τ represents the time shift.
Step 4.3, for R3τ(τ, τ) is fast fourier transformed to obtain a 1.5-dimensional frequency-domain weighted energy spectrum, which can be written as:
Figure GDA0003308682190000051
in the step 5, if the fault characteristic frequency of the inner ring, the outer ring, the rolling body and the retainer of the rolling bearing which protrude exists, the corresponding fault can be regarded as occurring, otherwise, the rolling bearing operates normally.
Has the advantages that: the invention has the technical effects that:
1) according to the rolling bearing fault diagnosis method based on the optimal dimension singular spectrum decomposition, provided by the invention, the performance of a singular spectrum decomposition algorithm is improved by optimizing parameter selection in the singular spectrum decomposition, and a rolling bearing measurement signal can be adaptively decomposed into a plurality of singular spectrum components with definite physical significance. The defect that the embedded dimension is selected according to the empirical criterion in the original singular spectrum decomposition algorithm is overcome, and the decomposition quality is effectively improved.
2) According to the 1.5-dimensional frequency domain weighted energy spectrum analysis method provided by the invention, the in-band noise and interference components of singular spectrum components are effectively inhibited through the idea of cooperation of the frequency domain weighted energy operator and the 1.5-dimensional spectrum, so that the fault characteristic frequency is enhanced, and the fault type is more accurately judged.
In addition, the method is simple and easy to implement, and is suitable for online monitoring and fault diagnosis of the rolling bearing.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a schematic structural view of a rolling bearing tester;
FIG. 3 is a timing signal and an envelope spectrum of a vibration acceleration signal for a bearing outer race fault;
FIG. 4 is a time domain diagram of 8 singular spectral components obtained by an optimal dimension singular spectral decomposition;
FIG. 5 is a graph of the variation of the objective function when decomposing the singular spectral components 3;
fig. 6 is a 1.5-dimensional frequency-domain weighted energy spectrum of the singular spectral component 3 of fig. five.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in detail with reference to the accompanying drawings and specific embodiments.
A rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition is shown in fig. 1, and the steps can be summarized as follows:
step 1: an acceleration sensor is arranged near the rolling bearing to collect vibration signals;
step 2: decomposing the acquired vibration signals by using an optimal dimension singular spectrum decomposition algorithm to obtain a plurality of singular spectrum components with relatively definite physical significance;
step 3, selecting singular spectrum components containing rich fault characteristic information as principal component components according to a kurtosis criterion;
step 4, calculating a 1.5-dimensional frequency domain weighted energy spectrum of the principal component;
and 5, observing whether an obvious peak value appears at the fault characteristic frequency from the 1.5-dimensional frequency domain weighted energy spectrum, thereby realizing accurate diagnosis of the rolling bearing fault.
The specific method is described as follows:
the specific steps of decomposing the vibration signal by using the optimal embedding dimension singular spectrum decomposition in the step 2 are as follows:
step 2.1, adaptively determining the input sequence v of the iterative decomposition(j)(n) an optimal embedding dimension M when constructing a trajectory matrix, wherein j is an iteration number;
step 2.2, constructing a new track matrix expression for the vibration signals: taking the sequence X (n) {1,2,3,4,5} as an example, if the embedding dimension M is set to 3, the corresponding matrix X is
Figure GDA0003308682190000061
In the formula, the left half part of a vertical line in the matrix X corresponds to a new track matrix in singular spectral decomposition, and the new way of constructing the track matrix in a ring shape is to fully enhance oscillation components in an original sequence and simultaneously enable the sequence to meet the rule that signal energy is reduced along with the increase of iteration times. In addition, two parameters of dimension and window length are required to be embedded into the traditional track matrix, one parameter is reduced by the new construction method, the two-dimensional optimization problem is reduced to the one-dimensional optimization problem, and the difficulty of parameter selection is greatly reduced;
inputting the iteration into a sequence v according to the track matrix construction mode(j)(N) constructing an M multiplied by N track matrix X, and performing singular value decomposition on the track matrix X to obtain:
Figure GDA0003308682190000062
wherein U is E.RP×PAnd u isk∈RP×1Is the kth column vector of U, V ∈ RN×NAnd v isk∈RN×1Is the kth column vector of V, Λ ∈ RP×NRepresenting a matrix of singular values, deltakIts k-th singular value.
Step 2.3, in the first iteration, if a large trend item is detected, only the first singular value of the singular value matrix and the left and right eigenvectors thereof are used for obtaining g(1)(n) that is
Figure GDA0003308682190000063
And g is(1)(n) from X1Is obtained from the diagonal average of (c). When the number of iterations j>1 hour, component sequence g(j)(n) must have a physically well-defined time scale, whose frequency components are mainly concentrated in the frequency band [ f ]max-Δf,fmax+Δf]In the meantime. Δ f denotes the bandwidth and is determined by the residual component vjThe power spectral density function of (n) is obtained by gaussian interpolation. Creating a subset I based on all feature groups with prominent dominant frequencies in the frequency band range for the left feature vector and one feature group with the largest contribution to the dominant peak energy of the selected modal componentsj＝{i1,…,ipThen through the matrix
Figure GDA0003308682190000077
Reconstructing the corresponding component sequence g by diagonal averaging(j)(n)；
Step 2.4, the component sequence g obtained in step 2.3 is subjected to(j)(n) separating from the original signal to obtain a residual component vj+1(n) of (a). Computing residueAnd (3) normalizing the mean square error NMSE of the residual component and the original sequence, if the NMSE is less than a given threshold th, stopping iteration, otherwise, repeating the steps 2.1-2.3, wherein the normalized mean square error is defined as:
Figure GDA0003308682190000071
step 2.1 adaptively determines the input sequence v for the next iteration(j)(n) the specific steps of the optimal embedding dimension are:
step 2.1.1: at the j-th iteration, according to the formula
Figure GDA0003308682190000072
A range of embedding dimensions M to be selected is determined,
Figure GDA0003308682190000073
for the current input sequence v(j)(n) the frequency value corresponding to the largest peak in the power spectral density function,
Figure GDA0003308682190000074
representing rounded numbers;
fs/fmaxthe number of points of a fault period of the main peak frequency in the iteration is represented, and when the embedding dimension is lower than the value, more serious modal aliasing can be caused, so that the obtained singular spectral components do not have clear physical significance; the embedding dimension cannot be larger than 1/3 of the sequence length, otherwise over-decomposition can be caused, and the side frequency and amplitude of the resonance band are weakened, so that the fault characteristics are not obvious.
Step 2.1.2: taking each value in the value interval of the embedding dimension M as the embedding dimension, and respectively carrying out comparison on the current input sequence v(j)(n) constructing a trajectory matrix and performing decomposition and reconstruction, i.e. performing steps 2.2-2.3, obtaining a series of singular spectral components
Figure GDA0003308682190000075
l is the length of the embedding dimension interval;
step 2.1.3: an Over Decomposition Index (ODI) of the obtained singular spectral components is calculated, the ODI being defined as:
Figure GDA0003308682190000076
in the formula (I), the compound is shown in the specification,
Figure GDA0003308682190000081
representation of SSCjIn
Figure GDA0003308682190000082
The amplitude of (a) of (b) is,
Figure GDA0003308682190000083
representing the input signal v of the current iteration(j)(n) in
Figure GDA0003308682190000084
The amplitude of (c).
The over-decomposition index is mainly used for avoiding overlarge embedding dimension setting, and requires that the main frequency amplitude obtained by decomposition is consistent with the amplitude in the original signal, otherwise, part of main frequency components appear in the residual signal, so that modal aliasing is generated by subsequent decomposition.
Step 2.1.4: using frequency band division method in empirical wavelet decomposition to input signal v of this iteration(j)The spectrum of (n) is adaptively band divided. After the frequency band division is finished, inputting a signal v according to the iteration(j)(n) main frequency
Figure GDA0003308682190000085
The frequency band of the channel is the singular spectral component SSC obtained in the iterationjClassifying into corresponding frequency bands, and then calculating a modal aliasing Index (MMI) of the obtained singular spectral component, wherein the MMI is defined as:
Figure GDA0003308682190000086
where I denotes the singular spectral component SSCjThe total number of spectral lines in the associated frequency band, T representing the singular spectral component SSCjTotal number of lines over the whole spectrum, AflAnd AftAre respectively the frequency flAnd ftThe corresponding amplitude.
The modal aliasing index avoids the modal aliasing problem caused by excessively small embedding dimension and excessively small setting through a method of dividing a frequency band in advance, and a frequency band division result does not represent a final decomposition result. The singular spectrum decomposition algorithm sets the frequency band range to be separated in the next iteration through a Gaussian interpolation model. The modal aliasing index is used to check whether the final decomposition result is consistent with the theoretical result.
In addition, it is disclosed that a harmonic component having a fixed frequency is decomposed into singular values corresponding to two singular values having close magnitudes in a singular value matrix, and thus values in a preset embedding dimension selection interval are even numbers in order to sufficiently decompose a frequency component. This can reduce the amount of computation by half without the resolution accuracy.
Step 2.1.5: will be provided with
Figure GDA0003308682190000087
Respective over-decomposition index ODI is multiplied by the modal aliasing index MMI to obtain the corresponding objective function value
Figure GDA0003308682190000088
Namely;
Figure GDA0003308682190000089
step 2.1.6: and considering the corresponding embedding dimension when the objective function value is maximum as the optimal embedding dimension.
When the value of the objective function is maximized, it is actually required to avoid both modal aliasing and over-decomposition problems.
In step 3, when the singular spectrum components including rich fault feature information are selected as principal component components according to a kurtosis criterion, the kurtosis value of each singular spectrum component is calculated, and the singular spectrum component with the highest kurtosis is selected as the principal component, where the kurtosis is defined as:
Figure GDA0003308682190000091
wherein mu is the mean value of the sequence x (n), and sigma is the standard deviation of the sequence x (n).
The step 4 of calculating the 1.5-dimensional frequency domain weighted energy spectrum of the principal component comprises the following specific steps:
step 4.1, carrying out frequency domain weighted energy operator operation on the principal component x (t) to obtain an instantaneous energy signal thereof
Figure GDA0003308682190000092
The frequency domain weighted energy operator is defined as:
Figure GDA0003308682190000093
where H [. cndot. ] represents the Hilbert transform, and Γ [ x (t) ], represents the frequency domain weighted energy operator of x (t).
Step 4.2, calculating instantaneous energy signal
Figure GDA0003308682190000094
The third order cumulant diagonal slice of (a) can be written as:
Figure GDA0003308682190000095
where E {. cndot } represents the mathematical expectation and τ represents the time shift.
Step 4.3, for R3τ(τ, τ) is fast fourier transformed to obtain a 1.5-dimensional frequency-domain weighted energy spectrum, which can be written as:
Figure GDA0003308682190000096
when transient impact occurs, the vibration amplitude changes rapidly, the frequency domain weighted energy operator can quickly and accurately track the change of the total energy of the signal, the transient characteristics of the signal are strengthened, and interference is effectively inhibited. Because the third-order cumulant of the zero-mean Gaussian noise is equal to zero, the noise component can be effectively inhibited after the transient energy signal doped with the noise component is transformed by the three-range cumulant, and therefore the 1.5 FWEO spectrum is combined with the advantages of a frequency domain weighted energy operator and a 1.5 FWEO spectrum.
The effectiveness of the invention is verified by experiments through an ABLT-1A type bearing life strengthening tester. Fig. 2 is a schematic structural diagram of a rolling bearing testing machine, which mainly comprises a testing head, a testing head seat, a transmission system, a loading system, a lubricating system, an electrical appliance control system, a computer monitoring system and the like. The test head is arranged in the test head seat, the traditional system transmits the motion of the motor, and the test shaft rotates at a certain rotating speed through the coupler; the loading system provides load required by the test, and the lubricating system enables the test shaft to be fully lubricated under normal conditions for the test; the electric control system provides power and electric protection and controls the actions of a motor, a hydraulic oil cylinder and the like; and recording test temperature and vibration information by the computer, and monitoring the running condition of the machine.
In the experiment, the tested object is a single-row deep groove ball bearing with model 6205, a groove with the depth of 0.2mm and the width of 0.15mm is processed on the outer ring of the bearing through an electric spark technology to simulate damage, a fault bearing is installed at a channel of the sensor 4, and sensors of channels 1,2 and 3 are installed on other three normal bearings. Under the conditions that the rotating speed is 17.5HZ and the sampling frequency is 10240HZ, 5s of data is collected, an eddy current sensor is used for picking up vibration signals, an electric signal is converted into a digital signal and transmitted to a PC through a data acquisition card, data acquisition and signal analysis are carried out by means of software platforms such as Labview, MATLAB and the like, and specific parameters of an experiment are shown in a table 1. The calculated characteristic frequency of the fault of the outer ring of the rolling bearing is fo＝62.5Hz。
TABLE 1 Motor-bearing System Experimental parameters and Main technical indices
Parameter name Value of
Driving rate 1050r/min
Sampling frequency 10240Hz
Sensor with a sensor element PCB electric eddy current sensor
Data acquisition system NI 9234
Acquisition procedure Labview、Matlab
PC platform Intel(R)core(TM)i7-9750H CPU@2.60GH
Bearing model 6205 Single-row deep groove ball bearing
Load(s) 0kN
Vibration signal data length 10240
Ambient temperature 5-40 degree
Size of the whole machine 1500*720*1300mm
Fig. 3 shows a time-series waveform and a hilbert envelope spectrum of successive 10240 points in data acquired by the sensor 1. The signal waveform contains a large amount of noise and is influenced by a plurality of interference sources among transmission components, so that obvious periodic impact characteristics are difficult to find; only the fault characteristic frequency of the outer ring of the rolling bearing can be observed in the envelope spectrum, the amplitude is small, and the existence of frequency multiplication can hardly be observed, which is probably caused by the fact that the signal transmission path is long, so that the signal attenuation is caused.
The vibration signal is decomposed using an optimal embedding dimension singular spectrum decomposition method to obtain 8 singular spectrum components, as shown in fig. 4.
Fig. 5 is a graph showing the change in the objective function value when the singular spectral component 3 is decomposed.
The singular spectrum component 3 with the largest kurtosis is selected as a principal component, the frequency domain weighted energy spectrum is shown in fig. 6 as 1.5, and the obvious outer ring fault characteristic frequency and 2-8 harmonics thereof can be seen from the frequency domain weighted energy spectrum, so that the outer ring fault of the bearing is judged, and therefore the practicability and the accuracy of the method in fault diagnosis of the rolling bearing are demonstrated.
The above description is only for the specific embodiment of the present invention, but the scope of the present invention is not limited thereto, and any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope of the present invention are included in the scope of the present invention.

Claims (4)

1. A rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition is characterized by comprising the following steps:
step 1: an acceleration sensor is arranged near the rolling bearing to collect vibration signals;
step 2: decomposing the acquired vibration signals by using an optimal dimension singular spectrum decomposition algorithm to obtain a plurality of singular spectrum components with definite physical significance;
step 3, selecting singular spectrum components containing rich fault characteristic information as principal component components according to a kurtosis criterion;
step 4, calculating a 1.5-dimensional frequency domain weighted energy spectrum of the principal component;
step 5, observing whether an obvious peak value appears at the fault characteristic frequency from the 1.5-dimensional frequency domain weighted energy spectrum, thereby realizing accurate diagnosis of the rolling bearing fault;
in the step 2, the specific steps of decomposing the vibration signal by using the optimal dimension singular spectrum decomposition algorithm are as follows:
step 2.1, adaptively determining the input sequence v of the iterative decomposition(j)(n) an optimal embedding dimension M when constructing a trajectory matrix, wherein j is an iteration number;
step 2.2, constructing a new track matrix expression for the vibration signals: taking the sequence X (n) {1,2,3,4,5} as an example, if the embedding dimension M is set to 3, the corresponding matrix X is
Figure FDA0003308682180000011
In the formula, the left half of the vertical line in the matrix X corresponds to a new trajectory matrix in the singular spectral decomposition;
inputting the iteration into a sequence v according to the track matrix construction mode(j)(N) constructing an M multiplied by N track matrix X, and performing singular value decomposition on the track matrix X to obtain:
Figure FDA0003308682180000012
wherein U is E.RP×PAnd u isk∈RP×1Is the kth column vector of U, V ∈ RN×NAnd v isk∈RN×1Is the kth column vector of V, Λ ∈ RP×NRepresenting a matrix of singular values, deltakIts k-th singular value;
step 2.3, in the first iteration, if a large trend item is detected, only the first singular value of the singular value matrix and the left and right eigenvectors thereof are used for obtaining g(1)(n) that is
Figure FDA0003308682180000013
And g is(1)(n) from X1Is obtained from the diagonal average of (1); when the number of iterations j>1 hour, component sequence g(j)(n) must have a physically well-defined time scale, whose frequency components are mainly concentrated in the frequency band [ f ]max-Δf,fmax+Δf]To (c) to (d); Δ f denotes the bandwidth and is determined by the residual component vj(n) the power spectral density function is obtained by gaussian interpolation; creating a subset I based on all feature groups with prominent dominant frequencies in the frequency band range for the left feature vector and one feature group with the largest contribution to the dominant peak energy of the selected modal componentsj＝{i1,…,ipThen through the matrix
Figure FDA00033086821800000211
Reconstructing the corresponding component sequence g by diagonal averaging(j)(n)；
Step 2.4, the component sequence g obtained in step 2.3 is subjected to(j)(n) separating from the original signal to obtain a residual component vj+1(n); calculating the normalized mean square error NMSE of the residual component and the original sequence, if NMSE is less than a given threshold th, stopping iteration, otherwise repeating steps 2.1-2.3, the normalized mean square error is defined as:
Figure FDA0003308682180000021
step 2.1 adaptively determines the input sequence v for the next iteration(j)(n) the specific steps of the optimal embedding dimension are:
step 2.1.1: at the j-th iteration, according to the formula
Figure FDA0003308682180000022
A range of embedding dimensions M to be selected is determined,
Figure FDA0003308682180000023
for the current input sequence v(j)(n) the frequency value corresponding to the largest peak in the power spectral density function,
Figure FDA0003308682180000024
representing rounded numbers;
step 2.1.2: taking each value in the value interval of the embedding dimension M as the embedding dimension, and respectively carrying out comparison on the current input sequence v(j)(n) constructing a trajectory matrix and performing decomposition and reconstruction, i.e. performing steps 2.2-2.3, obtaining a series of singular spectral components
Figure FDA0003308682180000025
l is the length of the embedding dimension interval;
step 2.1.3: calculating an index ODI of the obtained singular spectral components, the ODI being defined as:
Figure FDA0003308682180000026
in the formula (I), the compound is shown in the specification,
Figure FDA0003308682180000027
representation of SSCjIn
Figure FDA0003308682180000028
The amplitude of (a) of (b) is,
Figure FDA0003308682180000029
representing the input signal v of the current iteration(j)(n) in
Figure FDA00033086821800000210
The amplitude of (d);
step 2.1.4: using empirical wavelet decompositionFor the current iteration input signal v(j)(n) performing adaptive band division on the frequency spectrum; after the frequency band division is finished, inputting a signal v according to the iteration(j)(n) main frequency
Figure FDA0003308682180000031
The frequency band of the channel is the singular spectral component SSC obtained in the iterationjClassifying the obtained singular spectral components into corresponding frequency bands, and then calculating a modal aliasing index MMI of the obtained singular spectral components, wherein the MMI is defined as:
Figure FDA0003308682180000032
where I denotes the singular spectral component SSCjThe total number of spectral lines in the associated frequency band, T representing the singular spectral component SSCjThe total number of spectral lines over the entire spectrum,
Figure FDA0003308682180000033
and
Figure FDA0003308682180000034
are respectively the frequency flAnd ftA corresponding amplitude value;
step 2.1.5: will be provided with
Figure FDA0003308682180000035
Respective over-decomposition index ODI is multiplied by the modal aliasing index MMI to obtain the corresponding objective function value
Figure FDA0003308682180000036
Namely;
Figure FDA0003308682180000037
step 2.1.6: and selecting the corresponding embedding dimension as the optimal embedding dimension when the objective function value is maximum.
2. The rolling bearing fault diagnosis method based on the optimal dimension singular spectrum decomposition as claimed in claim 1, wherein: in step 3, when the singular spectrum components including rich fault feature information are selected as principal component components according to a kurtosis criterion, the kurtosis value of each singular spectrum component is calculated, and the singular spectrum component with the highest kurtosis is selected as the principal component, where the kurtosis is defined as:
Figure FDA0003308682180000038
wherein mu is the mean value of the sequence x (n), and sigma is the standard deviation of the sequence x (n).
3. The rolling bearing fault diagnosis method based on the optimal dimension singular spectrum decomposition as claimed in claim 1, wherein: the step 4 of calculating the 1.5-dimensional frequency domain weighted energy spectrum of the principal component comprises the following specific steps:
step 4.1, carrying out frequency domain weighted energy operator operation on the principal component x (t) to obtain an instantaneous energy signal thereof
Figure FDA0003308682180000039
The frequency domain weighted energy operator is defined as:
Figure FDA00033086821800000310
wherein H [. cndot. ] represents Hilbert transform, and Γ [ x (t) ], represents the frequency domain weighted energy operator of x (t);
step 4.2, calculating instantaneous energy signal
Figure FDA00033086821800000311
The third order cumulant diagonal slice of (a) is written as:
Figure FDA00033086821800000312
where E {. denotes the mathematical expectation and τ denotes the time shift;
step 4.3, for R3τ(τ, τ) performing fast fourier transform to obtain a 1.5-dimensional frequency-domain weighted energy spectrum, which is recorded as:
Figure FDA00033086821800000313
4. the rolling bearing fault diagnosis method based on the optimal dimension singular spectrum decomposition as claimed in claim 1, wherein: in the step 5, if the fault characteristic frequency of the inner ring, the outer ring, the rolling body and the retainer of the rolling bearing which protrude exists, the corresponding fault is considered to occur, otherwise, the rolling bearing operates normally.
CN202010046098.XA 2020-01-16 2020-01-16 Rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition Active CN111089726B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010046098.XA CN111089726B (en) 2020-01-16 2020-01-16 Rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010046098.XA CN111089726B (en) 2020-01-16 2020-01-16 Rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition

Publications (2)

Publication Number Publication Date
CN111089726A CN111089726A (en) 2020-05-01
CN111089726B true CN111089726B (en) 2021-12-03

Family

ID=70399401

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010046098.XA Active CN111089726B (en) 2020-01-16 2020-01-16 Rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition

Country Status (1)

Country Link
CN (1) CN111089726B (en)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111444893A (en) * 2020-05-06 2020-07-24 南昌航空大学 Fault diagnosis method for main shaft device of mine hoist
CN112345184B (en) * 2020-09-28 2022-01-11 同济大学 Structural earthquake damage identification method based on real-time monitoring data
CN112213640B (en) * 2020-11-17 2024-01-26 润电能源科学技术有限公司 Motor fault diagnosis method and related equipment thereof
CN112595515A (en) * 2020-12-04 2021-04-02 中国船舶工业综合技术经济研究院 Power shafting bearing fault detection method and system
CN112857804B (en) * 2021-02-09 2022-06-17 广东海洋大学 Rolling bearing fault diagnosis method, device, medium and computer equipment
CN113378063B (en) * 2021-07-09 2023-07-28 小红书科技有限公司 Method for determining content diversity based on sliding spectrum decomposition and content sorting method
CN113865867A (en) * 2021-08-20 2021-12-31 北京工业大学 Bearing fault diagnosis method based on amplitude characteristic singular value decomposition
CN116226754B (en) * 2023-05-04 2023-06-27 北京锦源汇智科技有限公司 Equipment health state assessment method and system based on equipment modeling
CN117290640B (en) * 2023-11-27 2024-01-26 天津领语未来智能科技有限公司 Singular spectrum harmonic decomposition method for nonlinear signal processing
CN117407824B (en) * 2023-12-14 2024-02-27 四川蜀能电科能源技术有限公司 Health detection method, equipment and medium of power time synchronization device

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105954031A (en) * 2016-06-29 2016-09-21 潍坊学院 Envelopment analysis method based on singular spectrum decomposition filtering
CN106053069A (en) * 2016-06-29 2016-10-26 潍坊学院 SSD, spectral kurtosis and smooth iteration envelope analysis method of antifriction bearing
CN106338385A (en) * 2016-08-25 2017-01-18 东南大学 Rotation machinery fault diagnosis method based on singular spectrum decomposition
CN108489719A (en) * 2018-04-09 2018-09-04 常州湖南大学机械装备研究院 A kind of combined failure of rotating machinery diagnostic method based on the unusual spectral factorizations of G-P
CN109932179A (en) * 2019-04-09 2019-06-25 东南大学 A kind of rolling bearing fault testing method based on the reconstruct of DS Adaptive spectra

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180284758A1 (en) * 2016-05-09 2018-10-04 StrongForce IoT Portfolio 2016, LLC Methods and systems for industrial internet of things data collection for equipment analysis in an upstream oil and gas environment

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105954031A (en) * 2016-06-29 2016-09-21 潍坊学院 Envelopment analysis method based on singular spectrum decomposition filtering
CN106053069A (en) * 2016-06-29 2016-10-26 潍坊学院 SSD, spectral kurtosis and smooth iteration envelope analysis method of antifriction bearing
CN106338385A (en) * 2016-08-25 2017-01-18 东南大学 Rotation machinery fault diagnosis method based on singular spectrum decomposition
CN108489719A (en) * 2018-04-09 2018-09-04 常州湖南大学机械装备研究院 A kind of combined failure of rotating machinery diagnostic method based on the unusual spectral factorizations of G-P
CN109932179A (en) * 2019-04-09 2019-06-25 东南大学 A kind of rolling bearing fault testing method based on the reconstruct of DS Adaptive spectra

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于奇异谱分解的振动故障诊断方法研究;李琛;《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》;20200115(第1期);第8-11、22-23、37-38页 *
基于能量算子的滚动轴承故障诊断研究;刘晶;《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》;20190115(第12期);第19-20、57-58、61、68-74页 *

Also Published As

Publication number Publication date
CN111089726A (en) 2020-05-01

Similar Documents

Publication Publication Date Title
CN111089726B (en) Rolling bearing fault diagnosis method based on optimal dimension singular spectrum decomposition
CN108760327B (en) Diagnosis method for rotor fault of aircraft engine
CN110046476B (en) Ternary binary fractal wavelet sparse diagnosis method for rolling bearing faults
CN105784366A (en) Wind turbine generator bearing fault diagnosis method under variable speed
Dou et al. A rule-based intelligent method for fault diagnosis of rotating machinery
Zhang et al. A new feature extraction approach using improved symbolic aggregate approximation for machinery intelligent diagnosis
Li et al. Multiscale symbolic diversity entropy: a novel measurement approach for time-series analysis and its application in fault diagnosis of planetary gearboxes
Cui et al. An investigation of rolling bearing early diagnosis based on high-frequency characteristics and self-adaptive wavelet de-noising
CN113375939B (en) Mechanical part fault diagnosis method based on SVD and VMD
Han et al. Roller bearing fault diagnosis based on LMD and multi-scale symbolic dynamic information entropy
Attoui et al. Vibration-based bearing fault diagnosis by an integrated DWT-FFT approach and an adaptive neuro-fuzzy inference system
Huang et al. An improved empirical wavelet transform method for rolling bearing fault diagnosis
CN112487882B (en) Method for generating non-sparse index-guided enhanced envelope spectrum based on spectrum coherence
CN113607415A (en) Bearing fault diagnosis method based on short-time stochastic resonance under variable rotating speed
CN109710955A (en) Based on LCD-recurrence quantification analysis rolling bearing fault diagnosis and health evaluating method
Saini et al. Predictive monitoring of incipient faults in rotating machinery: a systematic review from data acquisition to artificial intelligence
Lu et al. Bearing fault diagnosis based on clustering and sparse representation in frequency domain
CN111811819A (en) Bearing fault diagnosis method and device based on machine learning
Lu et al. Early fault warning and identification in condition monitoring of bearing via wavelet packet decomposition coupled with graph
Fa-jun et al. Compound fault diagnosis of gearbox based on wavelet packet transform and sparse representation classification
CN113899548A (en) Rolling bearing fault diagnosis method based on harmonic kurtosis spectrum
CN110147637B (en) Rub-impact fault diagnosis method based on wavelet and harmonic component greedy sparse identification
Yu et al. Rolling bearing fault feature extraction and diagnosis method based on MODWPT and DBN
CN113283028A (en) Fault diagnosis method for gear of gear box
CN114136604A (en) Rotary equipment fault diagnosis method and system based on improved sparse dictionary

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