Variable-speed bearing fault diagnosis method based on nonlinear mode decomposition

Abstract

The invention relates to a variable-speed bearing fault diagnosis method based on nonlinear modal decomposition, and the method comprises the steps: S1, obtaining a vibration signal of a rolling bearing at a variable speed, and calculating a fault feature order according to a bearing parameter; s2, performing low-pass filtering on the vibration signal to obtain a filtered signal; s3, extracting a non-linear mode component from the low-frequency part of the filtered vibration signal by using non-linear mode decomposition; s4, calculating the instantaneous frequency of the modal component by using a Hilbert method; s5, based on the instantaneous frequency, converting the vibration signal obtained in S1 into an angular domain pseudo-stationary signal by using an angular domain resampling method; and S6, performing envelope order analysis on the angular domain pseudo-stationary signal, solving a mode of an analysis signal, solving a Fourier spectrum for the mode of the analysis signal, comparing the Fourier spectrum with a fault feature order, and identifying a fault. According to the invention, the time-frequency analysis precision in a strong noise environment is improved, the accuracy of the rotating speed signal is ensured, and the accuracy of fault diagnosis is improved.

Description

Technical Field

[0001] The invention relates to the technical field of bearing fault diagnosis, and in particular to a variable speed bearing fault diagnosis method based on nonlinear modal decomposition. Background Art

[0002] Rolling bearings are commonly used components in rotating machinery and their health status is related to the reliability and stability of the equipment. However, when rolling bearings are running at variable speeds, they are more likely to fail due to the influence of working conditions, which in turn leads to serious safety problems and high maintenance costs. Therefore, it is extremely important to accurately diagnose the status of variable speed rolling bearings. Under variable speed operation, the fault characteristic frequency will be modulated by the speed, resulting in the inability to effectively judge the fault condition of the bearing by directly using the vibration signal envelope analysis.

[0003] The order tracking algorithm is one of the main methods for studying variable speed bearing fault diagnosis. Its core lies in tracking and obtaining the equal-angle sampling signal of the bearing. At present, the order tracking algorithm mainly includes hardware order tracking, computational order tracking and tachometer-free order tracking. The hardware order tracking method directly uses the photoelectric pulse angle encoder disk for equal-angle sampling to obtain the angular domain signal. This method has good real-time performance; however, its cost is high, the hardware equipment is complex, and it is not conducive to installation. Computational order tracking obtains the speed signal through the key phase device, and then calculates the equal-angle resampling time point according to the speed signal, and then obtains the angular domain signal through the interpolation algorithm; but the installation of the key phase device will increase the complexity of the equipment, especially in the transmission system. Tachometer-free order tracking mainly uses the time-frequency domain analysis method to directly extract the speed signal of the bearing from the vibration signal, and then obtains the angular domain signal by equal-angle sampling; this method does not require the installation of other measurement equipment, saving cost and equipment space. In a strong noise environment, most time-frequency analysis methods are difficult to guarantee the accuracy of the speed signal. Summary of the invention

[0004] In view of the shortcomings of the prior art, the present invention provides a variable speed bearing fault diagnosis method based on nonlinear mode decomposition (NMD), the purpose of which is to improve the accuracy of time-frequency analysis in a strong noise environment, ensure the accuracy of the speed signal, and improve the accuracy of fault diagnosis.

[0005] The technical solution adopted by the present invention is as follows:

[0006] A variable speed bearing fault diagnosis method based on nonlinear modal decomposition, comprising:

[0007] S1. Obtain the vibration signal of the rolling bearing under variable speed, and calculate the fault characteristic order according to the bearing parameters;

[0008] S2, performing low-pass filtering on the vibration signal to obtain a filtered signal;

[0009] S3, extracting nonlinear modal components from the low-frequency part of the filtered vibration signal using nonlinear modal decomposition;

[0010] S4. Calculate the instantaneous frequency of the modal component using the Hilbert method;

[0011] S5. Based on the instantaneous frequency, the vibration signal obtained in S1 is converted into a pseudo-stationary signal in the angular domain by using an angular domain resampling method;

[0012] S6. Analyze the pseudo-stationary signal in the angular domain through envelope order to obtain the modulus of the analytical signal, then obtain the Fourier spectrum of the modulus of the analytical signal, compare the Fourier spectrum with the fault characteristic order, and identify the fault.

[0013] Further technical solutions are:

[0014] In step S3, the nonlinear modal decomposition includes:

[0015] S31, obtaining a time-frequency representation of the signal using short-time Fourier transform or wavelet transform;

[0016] S32, extracting the time-frequency ridge curve in the time-frequency representation;

[0017] S33, after the time-frequency ridge curve is extracted, harmonic reconstruction is performed: the main harmonic is reconstructed using the ridge line method, and the subharmonic is reconstructed using the ridge line method and the direct method with the larger inconsistency coefficient;

[0018] S34, subharmonic test: Use the harmonic reconstruction method to obtain the hth subharmonic r h (t), the truth or falsity of the null hypothesis is determined by independence analysis, h ≥ 2;

[0019] S35, noise check: when the remaining signal is judged to be noise, the iteration is stopped.

[0020] In step S32, under the short-time Fourier transform, the peak sequence that maximizes the path functional is used as the time-frequency ridge curve, and the peak sequence expression is as follows:

[0021]

[0022] Among them, G s (·) is the time-frequency representation of short-time Fourier transform, ω p (t) represents the time-frequency ridge curve, std[·] represents the standard deviation, mean[·] represents the mean, Δω p (t n )=ω p (t n)-ω p (t n-1 ), t n Represents the nth moment.

[0023] In step S33, the time-frequency ridge curve ω is extracted p (t), the amplitude A(t) and phase frequency Extraction of harmonics can be achieved through reconstruction;

[0024] For the case where short-time Fourier transform is used to obtain time-frequency representation, the ridge reconstruction method is:

[0025]

[0026] For the case where wavelet transform is used to obtain time-frequency representation, the ridge reconstruction method is:

[0027]

[0028] Among them, W s (·) is the time-frequency representation of wavelet transform; δv d and δlnv d is the correction value for the discretization effect, is the Fourier transform of the window function of the short-time Fourier transform, ω ψ is the wavelet peak frequency, ψ * (t) is the complex conjugate function of the lognormal wavelet basis function ψ(t) of the wavelet transform.

[0029] The independence analysis of the null hypothesis includes: given the maximum time shift M, determining the proxy parameters of the r1(t) harmonics:

[0030]

[0031] Where N is the sampling length; f s is the sampling frequency; d is the number of time shift points; ΔT d represents a time interval;

[0032] For r h (t) subharmonic, h≥2, first shift the time-frequency signal forward by ΔT d / 2, and then use the harmonic reconstruction method to obtain the optimal parameters For r h (t) proxy parameters;

[0033] Using statistical parameters Quantization h The consistency of r(t) and r1(t):

[0034]

[0035] Provides an overall measure of the mutual consistency of the harmonics:

[0036]

[0037] Among them, w A ,w v ,w φ is the weight coefficient;

[0038] like If the confidence level of is less than the preset confidence level, the harmonic is a pseudo harmonic, otherwise it is a true harmonic; is the overall measure of consistency when ΔT0=0, is the overall measure of consistency when the number of time-shift points is d.

[0039] Step S35, the noise test includes:

[0040] Calculate the main harmonic noise test statistical feature D0 of this iteration;

[0041] Generate N s Proxy data, calculate the time-frequency representation corresponding to each proxy data, extract components from them, and calculate the corresponding noise test statistical characteristics When D s If the proportion of proxy data with value greater than D0 is higher than the set significance level, the null hypothesis of test noise is rejected.

[0042] In step S4, the instantaneous frequency of the modal component is calculated using the Hilbert method, including:

[0043] The signal is transformed by Hilbert to obtain the corresponding analytical signal;

[0044] Calculate the amplitude and phase of the analytical signal;

[0045] The inverse of the phase with respect to time is calculated as the instantaneous frequency.

[0046] In step S5, based on the instantaneous frequency, the vibration signal obtained in step S1 is converted into an angular domain pseudo-stationary signal using an angular domain resampling method, including:

[0047] According to the equal-angle sampling interval Δθ, the original vibration signal x(t) is resampled in the angular domain based on cubic spline interpolation to obtain the total rotation angle at time t(i):

[0048]

[0049] Among them, f r(t) is the speed signal, i.e., the instantaneous frequency obtained by S4; i represents the number of the midpoint of the time series;

[0050] Use cubic spline interpolation to fit the angle θ and time t between t(i) and t(i+1):

[0051] θ i (t) = a3(tt(i)) 3 +a2(tt(i)) 2 +(tt(t))+a0

[0052] Among them, a0, a2, and a3 are fitting coefficients respectively;

[0053] When θ(i)≤kΔθ<θ(i+1), the time T(k) of the kth equal-angle sampling in the time domain vibration signal is calculated:

[0054]

[0055] Interpolation operation on x(t) can obtain the vibration amplitude A(k) at time point T(k), thereby obtaining an angular domain signal with an angle sequence of Δθ, 2Δθ, …, kΔθ and an amplitude sequence of A(1), A(2), …, A(k).

[0056] In step S2, the cut-off frequency during low-pass filtering is determined according to the speed operating range of the bearing.

[0057] In step S1, the fault feature order includes the inner ring fault feature order:

[0058]

[0059] Outer race fault characteristic order:

[0060]

[0061] Rolling element fault characteristic order:

[0062]

[0063] Among them, z is the number of rolling elements, d is the rolling element diameter, D is the pitch diameter, and α is the contact angle.

[0064] The beneficial effects of the present invention are as follows:

[0065] The present invention extracts low-frequency nonlinear modal components in the rotation process of the bearing by the NMD method, calculates the instantaneous frequency by the Hilbert method, obtains the speed information of the bearing, converts the variable speed non-stationary signal in the time domain into a pseudo-stationary signal in the angular domain by the angular domain resampling method, and realizes accurate judgment of the fault of the variable speed bearing by envelope order analysis. The method of the present invention has good noise resistance, does not need to install other measuring equipment, and saves cost and equipment space.

[0066] Other features and advantages of the present invention will be set forth in the following description or may be learned by practicing the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 It is a bearing speed curve corresponding to the simulation signal of the embodiment of the present invention.

[0068] Figure 2 This is a simulated impact waveform according to an embodiment of the present invention.

[0069] Figure 3 It is a simulation signal waveform diagram of an embodiment of the present invention.

[0070] Figure 4 It is a nonlinear modal analysis waveform diagram extracted by an embodiment of the present invention.

[0071] Figure 5 is the instantaneous frequency calculated by Hilbert in the embodiment of the present invention.

[0072] Figure 6 It is a comparison diagram of the speed curve extracted by the embodiment of the present invention and the actual speed curve.

[0073] Figure 7 is the angular domain waveform of the simulation signal according to the embodiment of the present invention.

[0074] Figure 8 It is the result of envelope order analysis of the simulation signal of the embodiment of the present invention.

[0075] Fig. 9 This is the time domain waveform of the inner ring fault of the Ottawa variable speed bearing, which is a verification example of the present invention.

[0076] Fig.10 It is the Ottawa variable speed bearing speed signal extracted by the verification example of the present invention.

[0077] Fig.11 It is the angular domain signal of the Ottawa variable speed bearing in the verification example of the present invention.

[0078] Fig.12 It is the envelope order analysis result of the Ottawa variable speed bearing of the verification example of the present invention.

[0079] Fig.13It is a schematic diagram of the process of the present invention. DETAILED DESCRIPTION

[0080] The specific implementation of the present invention is described below with reference to the accompanying drawings.

[0081] See also Fig.13 , a variable speed bearing fault diagnosis method based on nonlinear mode decomposition (NMD) of this embodiment includes:

[0082] S1. Obtain the vibration signal of the rolling bearing under variable speed, and calculate the fault characteristic order according to the bearing parameters;

[0083] Order analysis can eliminate the influence of speed on fault characteristics. The order is defined as:

[0084]

[0085] Where, f is the bearing rotation frequency, f r is the rotation frequency. By combining the order and the fault characteristic frequency of the bearing, the fault order characteristics of the bearing can be obtained as follows:

[0086] Inner race fault characteristic order:

[0087]

[0088] Outer race fault characteristic order:

[0089]

[0090] Rolling element fault characteristic order:

[0091]

[0092] Among them, z is the number of rolling elements, d is the rolling element diameter, D is the pitch diameter, and α is the contact angle.

[0093] Specifically, when the outer ring of the bearing fails under variable speed conditions, the effects of continuous unit pulse excitation are as follows:

[0094]

[0095] Among them, F(t) is the impact generated at time t, δ is the impact response, and the above formula means t = T k When a shock occurs; f r (t) is the bearing rotation frequency; k is the number of pulse excitations; T k is the time when the kth pulse excitation is generated; or is the fault order;

[0096] Unit impulse response of the system The time domain impact signal x(t) generated by the bearing fault is:

[0097]

[0098] Where, ξ is the damping ratio; f n is the resonance frequency; A k is the amplitude, when f r < <f n When A k The change is very small and can be regarded as a constant. The numerical simulation model is expressed as:

[0099]

[0100] The simulation signal y(t) includes a periodic transient impulse signal x(t), a harmonic component d(t) and a noise signal e(t). n =3000Hz; A k =1, ξ=0.23; or=5. d(t) is used to simulate the harmonic components in the signal, which is modulated by the speed. e(t) is a Gaussian white noise with a signal-to-noise ratio of -13dB generated based on x(t). The sampling frequency of the simulated signal y(t) is fs=10000Hz, and the number of sampling points N=10000.

[0101] In order to better simulate the fluctuation of bearing speed, f r =10+5sin(4πt), its changing trend is as follows Figure 1 As shown, the simulated impact waveform is as follows Figure 2 As shown, the simulation signal waveform is as follows Figure 3 shown.

[0102] S2, low-pass filtering the vibration signal to obtain a filtered signal;

[0103] Since the operating frequency of the transmission chain bearing is usually less than 50Hz, low-pass filtering is performed with a cutoff frequency of 50Hz. Determining the cutoff frequency according to the speed operating range of the bearing can effectively prevent the influence of noise.

[0104] S3, extracting nonlinear modal components from the low-frequency part of the filtered vibration signal using nonlinear modal decomposition, comprising the following steps:

[0105] S31. Obtaining the time-frequency representation of the signal using short-time Fourier transform (Windowed Fourier transform, WFT) or wavelet transform (Wavelet transform, WT);

[0106]

[0107] Among them, G s (ω,t),W s(ω,t) are the time-frequency representations of WFT and WT respectively; ω is the frequency; u and ξ are the representations of time t in the time domain and frequency domain respectively; s + (u), are the representations of the signal s(t) before transformation in the time domain and frequency domain respectively; g(t) is the window function of WFT, is the Fourier transform of g(t);

[0108] ω ψ is the wavelet peak frequency; ψ(t) is the lognormal wavelet basis function of WT, ψ * (t) represents the complex conjugate function of ψ(t), represents the complex conjugate function of the Fourier transform of ψ(t), that is, the "pointed hat" in the formula represents the Fourier transform, and the "*" represents the complex conjugate function, and both represent the complex conjugate function of the Fourier transform.

[0109] S32, extracting the time-frequency ridge curve in the time-frequency representation;

[0110] Under the short-time Fourier transform, the peak sequence that maximizes the path functional is taken as the time-frequency ridge curve, and the peak sequence expression is as follows:

[0111]

[0112] Among them, G s (·) is the time-frequency representation of short-time Fourier transform, ω p (t) represents the time-frequency ridge curve, std[·] represents the standard deviation, mean[·] represents the mean, Δω p (t n )=ω p (t n )-ω p (t n-1 ), t n Represents the nth moment.

[0113] Since the wavelet transform is a logarithmic frequency resolution, we can replace ω in the above formula p Change to lnω p , Δω p Change to Δlnω p , the corresponding time-frequency ridges can be extracted according to this standard.

[0114] S33. After the time-frequency ridge curve is extracted, harmonic reconstruction is performed: the main harmonics are reconstructed using the ridge line method, and the subharmonics are reconstructed using the ridge line method and the direct method with a larger inconsistency coefficient.

[0115] Specifically, extract the time-frequency ridge curve ω p (t), the amplitude A(t) and phase frequency The extraction of harmonics is realized; the reconstruction is divided into ridge method and direct method.

[0116] For the case where short-time Fourier transform is used to obtain time-frequency representation, the ridge reconstruction method is:

[0117]

[0118] For the case where wavelet transform is used to obtain time-frequency representation, the ridge reconstruction method is:

[0119]

[0120] Among them, W s (·) is the time-frequency representation of wavelet transform; δv d and δlnv d is the correction value of the discretization effect, which can be obtained by parabolic interpolation; ψ * is the complex conjugate function of ψ.

[0121] For the case where short-time Fourier transform is used to obtain time-frequency representation, the direct reconstruction method is:

[0122]

[0123] Among them, Cg is just a calculation parameter, and Re is the real part of the complex number;

[0124] For the case of using wavelet transform to obtain time-frequency representation, the direct reconstruction method is:

[0125]

[0126] For the above two reconstruction methods, the reconstruction method with a larger inconsistency coefficient is selected. The inconsistency coefficient of the harmonic parameters can be expressed as:

[0127]

[0128] In the formula, k A , k v To adjust the coefficient, k is summarized based on experience A =2, k v =4 or both are 1. In addition, since the ridge curve method has better noise robustness, and the direct method allows more accurate tracking of the time variability of the parameters at a low noise level, this embodiment uses the ridge curve method to reconstruct the main harmonic r1(t) and the subharmonic r2(t). h (t) Reconstructed using the method with a larger inconsistency coefficient between the ridge curve method and the direct method.

[0129] The two time-frequency representation methods, WFT and WT, have different characteristics in different environments. Therefore, in the process of obtaining time-frequency representation, you can select according to the following empirical formula: When R<1, use WFT; when R≥1, use WT.

[0130]

[0131] Among them, std is the calculated standard deviation.

[0132] S34, subharmonic test: Use the harmonic reconstruction method to obtain the hth subharmonic r h (t), the truth or falsity of the null hypothesis is determined by independence analysis, h ≥ 2;

[0133] After finding the time-frequency ridge curve of the main harmonic, its hth harmonic is generally located near the time-frequency ridge curve of the main harmonic. After obtaining the hth harmonic using the harmonic reconstruction method, considering the actual physical meaning, it is necessary to determine whether it is a real harmonic. The NMD algorithm uses the independence analysis of the null hypothesis to identify the subharmonic r h (t)(h≥2) is true or false.

[0134] The independence analysis of the null hypothesis includes:

[0135] Given the maximum time shift M, determine the proxy parameters for the r1(t) harmonics:

[0136]

[0137] Where N is the sampling length; f s is the sampling frequency; d is the number of time shift points; ΔT d is the time interval;

[0138] For r h (t) subharmonic, h≥2, first shift the time-frequency signal forward by ΔT d / 2, and then use the harmonic reconstruction method to obtain the optimal parameters For r h (t) proxy parameters;

[0139] Using statistical parameters Quantization h The consistency of r(t) and r1(t):

[0140]

[0141] Provides an overall measure of the mutual consistency of the harmonics:

[0142]

[0143] Among them, w A ,w v,w φ is the weight coefficient;

[0144] like If the confidence level of is less than the preset confidence level, the harmonic is a pseudo harmonic, otherwise it is a true harmonic;

[0145] is the overall measure of consistency when Δt0=0, is the overall measure of consistency when the number of time-shift points is d.

[0146] S35, noise check: when the remaining signal is judged as noise, the iteration is stopped;

[0147] In the NMD decomposition process, the condition for stopping the iterative decomposition is that the remaining signal is judged to be noise. A ,α v ) is calculated as follows:

[0148]

[0149] In the formula, is the Fourier entropy spectrum of the amplitude; is the Fourier entropy spectrum of frequency.

[0150] During the inspection:

[0151] Calculate the main harmonic noise test statistical characteristics D0(α A ,α v );

[0152] Generate N s Proxy data, calculate the time-frequency representation corresponding to each proxy data, extract components from them, and calculate the corresponding noise test statistical characteristics When D s If the proportion of proxy data with value greater than D0 is higher than the set significance level, the null hypothesis of test noise is rejected.

[0153] This embodiment uses NMD to extract nonlinear modal components in the low-frequency part of the signal. The simulation signal is processed using NMD, and the obtained low-frequency modal components are composed of one main harmonic and two subharmonics. The extracted nonlinear modal components are the sum of three harmonics. The results are as follows: Figure 4 shown.

[0154] S4. Calculate the instantaneous frequency of the modal component using the Hilbert method, including:

[0155] The modal component signal x(t) is transformed by Hilbert:

[0156]

[0157] Get the corresponding analysis signal That is, the imaginary part of z(t) is the Hilbert transform of its real part. Where A is the amplitude and φ is the phase, which are calculated as:

[0158]

[0159] Calculate the instantaneous frequency as the inverse of the phase with respect to time:

[0160]

[0161] The result of calculating the instantaneous frequency according to the Hilbert method in this embodiment is as follows: Figure 5 As shown, compare it with the actual transfer frequency f r The comparison results are as follows Figure 6 As shown, it can be seen that NMD can effectively extract the speed information from the bearing signal.

[0162] S5. Based on the instantaneous frequency, the original vibration signal obtained in S1 (a non-stationary signal with variable speed in the time domain) is converted into a pseudo-stationary signal in the angular domain by using an angular domain resampling method, including:

[0163] According to the equal-angle sampling interval Δθ, the original vibration signal x(t) is resampled in the angular domain based on cubic spline interpolation to obtain the total rotation angle at time t(i):

[0164]

[0165] Among them, f r (t) is the speed signal, i represents the number of the time series point;

[0166] Use cubic spline interpolation to fit the angle θ and time t between t(i) and t(i+1):

[0167] θ i (t) = a3(tt(i)) 3 +a2(tt(i)) 2 +(tt(i))+a0

[0168] Among them, a0, a2, and a3 are fitting coefficients respectively;

[0169] When θ(i)≤kΔθ<θ(i+1), the time T(k) of the kth equal-angle sampling in the time domain vibration signal is calculated:

[0170]

[0171] Interpolation operation on the vibration signal x(t) can obtain the vibration amplitude A(k) at time point T(k), thereby obtaining an angular domain signal with an angle sequence of Δθ, 2Δθ, …, kΔθ and an amplitude sequence of A(1), A(2), …, A(k).

[0172] This embodiment uses the rotation frequency signal extracted by NMD to convert the simulation signal into a pseudo-stationary signal in the angular domain. The number of resampling points for one rotation of the bearing is 360. The angular domain signal is as follows: Figure 7 shown.

[0173] S6. Analyze the pseudo-stationary signal in the angular domain obtained in S5 through envelope order: perform Hilbert transform on the angular domain signal to obtain the modulus of the analytical signal, then obtain the Fourier spectrum of the modulus of the analytical signal, compare the Fourier spectrum with the fault characteristic order of step S1, and identify the fault.

[0174] This embodiment performs envelope order analysis on the bearing angular domain signal. The results are as follows: Figure 8 The envelope order spectrum can also clearly find the characteristic order of the transient impact signal, which is consistent with the fault characteristic order or=5 of the simulation signal, and effectively diagnoses the bearing fault.

[0175] Verification Example

[0176] In order to further verify the method of this embodiment, verification is performed using a measured signal, including the following steps:

[0177] Step S1, obtain the vibration signal of the rolling bearing under variable speed, and calculate the fault feature order according to the bearing parameters. The University of Ottawa data set is a variable speed bearing data set, the bearing model is ER16K, and its fault order is calculated as:

[0178] Inner race fault:

[0179]

[0180] Outer race fault:

[0181]

[0182] Rolling element failure:

[0183]

[0184] Vibration data is measured by an acceleration sensor with a sampling frequency of 200kHz and a single sample sampling duration of 10 seconds. The data set consists of two conditions: speed change mode and fault type. Bearing inner race fault data is selected for analysis. Since the high-frequency part of the data is mainly noise components, and too high a sampling frequency will result in a large amount of data and consume a lot of computing time. Therefore, the original data is resampled with a resampling frequency of 10000Hz, and 20,000 sampling points between 5 and 7 seconds are intercepted as the analysis data for this verification experiment. The time domain waveform is as follows: Fig. 9 shown.

[0185] Step S2: Since the operating rotation frequency of the transmission chain bearing of the University of Ottawa is less than 50 Hz, low-pass filtering is performed with a cut-off frequency of 50 Hz.

[0186] Step S3: Using NMD to extract nonlinear modal components in the low-frequency part of the signal.

[0187] Step S4, and calculate the instantaneous frequency according to the Hilbert method, the result is as follows Fig.10 shown.

[0188] Step S5: Use the angle domain resampling method to convert the time domain variable speed non-stationary signal into an angle domain pseudo-stationary signal. Use the rotation frequency signal extracted by NMD to convert the simulation signal into an angle domain pseudo-stationary signal. The number of resampling points for one rotation of the bearing is 180. The angle domain signal is as follows: Fig.11 shown.

[0189] Step S6: Mining fault features through envelope order analysis. Envelope order analysis is performed on the bearing angular domain signal. The results are as follows: Fig.12 The envelope order spectrum can also clearly reveal the characteristic order of the transient impact signal, which is consistent with the characteristic order of the inner ring fault of the simulation signal or i =5.43 is very consistent, which effectively diagnoses the bearing fault.

[0190] Those skilled in the art can understand that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention is described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A variable speed bearing fault diagnosis method based on nonlinear modal decomposition, characterized in that: include: S1. Obtain the vibration signal of the rolling bearing under variable speed, and calculate the fault characteristic order according to the bearing parameters; S2, performing low-pass filtering on the vibration signal to obtain a filtered signal; S3, extracting nonlinear modal components from the low-frequency part of the filtered vibration signal using nonlinear modal decomposition; S4. Calculate the instantaneous frequency of the modal component using the Hilbert method; S5. Based on the instantaneous frequency, the vibration signal obtained in S1 is converted into a pseudo-stationary signal in the angular domain by using an angular domain resampling method; S6. Analyze the pseudo-stationary signal in the angular domain through envelope order to obtain the modulus of the analytical signal, then obtain the Fourier spectrum of the modulus of the analytical signal, compare the Fourier spectrum with the fault characteristic order, and identify the fault.

2. The method according to claim 1, characterized in that In step S3, the nonlinear modal decomposition includes: S31, obtaining a time-frequency representation of the signal using short-time Fourier transform or wavelet transform; S32, extracting the time-frequency ridge curve in the time-frequency representation; S33, after the time-frequency ridge curve is extracted, harmonic reconstruction is performed: the main harmonic is reconstructed using the ridge line method, and the subharmonic is reconstructed using the ridge line method and the direct method with the larger inconsistency coefficient; S34, subharmonic test: Use the harmonic reconstruction method to obtain the hth subharmonic r h (t), the truth or falsity of the null hypothesis is determined by independence analysis, h ≥ 2; S35, noise check: when the remaining signal is judged to be noise, the iteration is stopped.

3. The method according to claim 2, characterized in that In step S32, under the short-time Fourier transform, the peak sequence that maximizes the path functional is used as the time-frequency ridge curve, and the peak sequence expression is as follows: Among them, G s (·) is the time-frequency representation of short-time Fourier transform, ω p (t) represents the time-frequency ridge curve, std[·] represents the standard deviation, mean[·] represents the mean, Δω p (t n )=ω p (t n )-ω p (t n-1 ), t n Represents the nth moment.

4. The method according to claim 3, characterized in that In step S33, the time-frequency ridge curve ω is extracted p (t), the amplitude A(t) and phase frequency Extraction of harmonics can be achieved through reconstruction; For the case where short-time Fourier transform is used to obtain time-frequency representation, the ridge reconstruction method is: For the case where wavelet transform is used to obtain time-frequency representation, the ridge reconstruction method is: Among them, W s (·) is the time-frequency representation of wavelet transform; δv d and δlnv d is the correction value for the discretization effect, is the Fourier transform of the window function of the short-time Fourier transform, ω ψ is the wavelet peak frequency, ψ * (t) is the complex conjugate function of the lognormal wavelet basis function ψ(t) of the wavelet transform.

5. The method according to claim 4, characterized in that The independence analysis of the null hypothesis includes: given the maximum time shift M, determining the proxy parameters of the r1(t) harmonics: Where N is the sampling length; f s is the sampling frequency; d is the number of time shift points; ΔT d represents a time interval; For r h (t) subharmonic, h≥2, first shift the time-frequency signal forward by ΔT d / 2, and then use the harmonic reconstruction method to obtain the optimal parameters For r h (t) proxy parameters; Using statistical parameters Quantization h The consistency of r(t) and r1(t): Provides an overall measure of the mutual consistency of the harmonics: Among them, w A ,w v ,w φ is the weight coefficient; like If the confidence level of is less than the preset confidence level, the harmonic is a pseudo harmonic, otherwise it is a true harmonic; is the overall measure of consistency when ΔT0=0, is the overall measure of consistency when the number of time-shift points is d.

6. The method according to claim 2, characterized in that Step S35, the noise test includes: Calculate the main harmonic noise test statistical feature D0 of this iteration; Generate N s Proxy data, calculate the time-frequency representation corresponding to each proxy data, extract components from them, and calculate the corresponding noise test statistical characteristics When D s If the proportion of proxy data with value greater than D0 is higher than the set significance level, the null hypothesis of test noise is rejected.

7. The method according to claim 1, characterized in that In step S4, the instantaneous frequency of the modal component is calculated using the Hilbert method, including: The signal is transformed by Hilbert to obtain the corresponding analytical signal; Calculate the amplitude and phase of the analytical signal; The inverse of the phase with respect to time is calculated as the instantaneous frequency.

8. The method according to claim 1, characterized in that: In step S5, based on the instantaneous frequency, the vibration signal obtained in step S1 is converted into an angular domain pseudo-stationary signal using an angular domain resampling method, including: According to the equal-angle sampling interval Δθ, the original vibration signal x(t) is resampled in the angular domain based on cubic spline interpolation to obtain the total rotation angle at time t(i): Among them, f r (t) is the speed signal, i.e., the instantaneous frequency obtained by S4; i represents the number of the midpoint of the time series; Use cubic spline interpolation to fit the angle θ and time t between t(i) and t(i+1): θ i (t)=a3(t-t(i)) 3 +a2(t-t(i)) 2 +(t-t(i))+a0 Among them, a0, a2, and a3 are fitting coefficients respectively; When θ(i)≤kΔθ<θ(i+1), the time T(k) of the kth equal-angle sampling in the time domain vibration signal is calculated: Interpolation operation on x(t) can obtain the vibration amplitude A(k) at time point T(k), thereby obtaining an angular domain signal with an angle sequence of Δθ, 2Δθ, …, kΔθ and an amplitude sequence of A(1), A(2), …, A(k).

9. The method according to claim 1, characterized in that: In step S2, the cut-off frequency during low-pass filtering is determined according to the speed operating range of the bearing.

10. The method according to claim 1, characterized in that In step S1, the fault feature order includes the inner ring fault feature order: Outer race fault characteristic order: Rolling element fault characteristic order: Among them, z is the number of rolling elements, d is the rolling element diameter, D is the pitch diameter, and α is the contact angle.

Images