Variable speed fault feature enhancement method based on sparse feature adaptive convolution spectrum amplitude modulation

Abstract

The invention provides a variable speed fault feature enhancement method based on sparse feature adaptive convolution spectrum amplitude modulation. According to the method, parameters in an interval are selected to construct wavelets, cross-correlation analysis is carried out on the wavelets and signals, and a single wavelet with the maximum correlation coefficient is selected as a feature dictionary; tracking and denoising the signal convolution basis by using the feature dictionary, and analyzing to obtain a sparse signal; amplitude modulation processing is carried out on the sparse signal spectrum to obtain an amplitude correction component; and constructing an order envelope spectrum of the amplitude correction component, and realizing the diagnosis of the variable-speed operation fault of the bearing by comparing the feature order in the order envelope spectrum with the theoretical fault feature order of the bearing.

Description

Technical Field

[0001] The present invention relates to the field of rolling bearing fault diagnosis, and proposes an amplitude modulation algorithm based on sparse feature adaptive convolution spectrum, which can enhance the fault-related characteristic components of the vibration signal and is suitable for the diagnosis of bearing faults in variable speed operation. The present invention constructs a wavelet by selecting parameters in an interval and performs cross-correlation analysis with the signal, selects a single wavelet with the maximum correlation coefficient as a feature dictionary; uses the feature dictionary to perform convolution basis pursuit denoising on the signal and obtains a sparse signal by analysis; performs amplitude modulation processing on the sparse signal spectrum to obtain an amplitude correction component; constructs an order envelope spectrum of the amplitude correction component to achieve fault diagnosis of bearings in variable speed operation. Background Art

[0002] Rolling bearings are core components that support rotating mechanical equipment and are crucial to the reliable and efficient operation of the equipment. In industrial scenarios, bearings usually operate in harsh and complex working conditions and are prone to failure after long periods of operation. They need to be monitored to detect potential failures in time to avoid affecting production and causing losses. The state of the bearing is often identified by collecting vibration signals, but due to the interference of background noise, the internal fault information of the time domain signal may be submerged by redundant components, which poses a huge challenge to equipment fault diagnosis. Therefore, effectively mining fault information from complex bearing signals is the technical key to equipment variable speed fault diagnosis. Summary of the invention

[0003] When a local fault occurs in the bearing, the surface fault defects of the bearing will trigger a series of periodic transient shocks in the collected vibration signal. The fault pulse characteristics contained in the collected signal can be analyzed to determine the type of bearing fault. However, in actual industrial scenarios, due to interference from other vibrating components and signal transmission paths in the equipment, the internal vibration shock characteristics of the signal are submerged by noise, affecting the recognition effect of the fault pulse characteristics. In addition, rotating mechanical equipment is in a variable speed operation state during the start-up and shutdown process, and the fault shock mode in the signal may be modulated and aliased, which makes it more difficult to identify and extract the shock characteristics.

[0004] In order to solve the problem of difficulty in identifying fault features of variable speed bearing signals, the present invention proposes a variable speed fault feature enhancement method based on sparse feature adaptive convolution spectrum modulation to capture the time domain fault feature information of variable speed signals. The main difficulty of the present invention is to realize matching optimization of double damped wavelet feature dictionary for convolution basis pursuit denoising and combine spectrum modulation to realize fault feature enhancement.

[0005] The steps of the invention are as follows

[0006] S1: Bearing signal acquisition: Use an acceleration sensor to collect the bearing's variable speed vibration signal; use an incremental encoder to synchronously collect the bearing's key phase pulse signal for angular domain resampling of subsequent order analysis.

[0007] S2: Feature dictionary construction: The double-damped wavelet has the characteristics of bilateral attenuated oscillation, which can effectively characterize the time-domain impact waveform caused by defects. Set the key parameter selection interval that affects the wavelet waveform, use the cross-correlation analysis method combined with the particle swarm optimization algorithm to match the maximum correlation coefficient of the wavelet and the vibration signal, search for the optimal parameters of the double-damped wavelet and construct the feature dictionary.

[0008] S3: Convolutional basis pursuit denoising: A convolutional basis pursuit denoising model is constructed using the constructed feature dictionary and vibration signal, and the model is iteratively solved using the alternating direction multiplication method to obtain sparse coefficients.

[0009] S4: Spectral amplitude modulation: The obtained sparse coefficients are subjected to spectral amplitude modulation processing to obtain amplitude correction components of the sparse coefficients.

[0010] S5: Amplitude correction component characteristic analysis: The amplitude correction component is resampled in the angular domain according to the key phase pulse signal and envelope analysis is performed. The fault type is determined according to the main characteristic components displayed in the drawn order envelope spectrum. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 It is a schematic diagram of the analysis process of the invention;

[0012] Figure 2 It is the time domain diagram of the outer ring vibration signal of the faulty bearing;

[0013] Figure 3 It is the time domain diagram of the key phase pulse of the outer ring of the faulty bearing;

[0014] Figure 4 It is the order envelope spectrum analysis result of the original signal of the outer ring of the faulty bearing;

[0015] Figure 5 It is the result of adjusting the key parameters on the double-damped wavelet waveform;

[0016] Figure 6 The present invention obtains a sparse coefficient time domain graph;

[0017] Figure 7 This is the final order envelope spectrum analysis result obtained by the present invention. DETAILED DESCRIPTION

[0018] The present invention proposes a variable speed fault feature enhancement method based on sparse feature adaptive convolution spectrum amplitude modulation. Figure 1 Schematic diagram of the analysis process.

[0019] S1: Vibration signal acquisition: Use an acceleration sensor to collect the variable speed bearing vibration signal, expressed as Figure 2 The dimension is 1×n time series x, n is the number of sampling points, and the sampling frequency of the vibration signal is 200kHz. Use incremental encoder to synchronously collect Figure 3 The bearing key phase pulse signal is a time series S with a dimension of 1×n. Figure 4 No obvious characteristic components can be observed in the order envelope spectrum analysis results of the vibration signal.

[0020] S2: Feature dictionary construction: The double damped wavelet g(t) is expressed as:

[0021]

[0022] Where t is time, τ is the time shift, Sgn() is the sign function, ω is the angular frequency, ζ1 is the attenuation damping coefficient, ζ2 is the enhancement damping coefficient, g1(t) is the attenuation section of the wavelet, g2(t) is the enhancement section of the wavelet, and f d is the natural frequency. Among them, the attenuation damping coefficient ζ1, the enhancement damping coefficient ζ2 and the natural frequency f d are the key parameters that affect the wavelet waveform, such as Figure 5 Selecting different combinations of key parameters to construct double-damped wavelets has different forms. The pre-selected intervals of the attenuation damping coefficient ζ1 and the enhancement damping coefficient ζ2 of the double-damped wavelet are set to (0,1], and the natural frequency f d The pre-selected interval is [500,3000]. The time t interval is set to [0,0.025], and the time shift τ is set to 0.0125.

[0023] The cross-correlation analysis method is used to evaluate the similarity between the double damped wavelet and the internal impact component of the vibration signal. The cross-correlation function value is calculated by sliding the double damped wavelet g(t) as a window on the vibration signal at different positions of the vibration signal time series x:

[0024]

[0025] Where τ0 is the lag position of the window sliding of the double damped wavelet g(t) on the time series x. The lag position of the maximum cross-correlation function obtained by window sliding is used to calculate the normalized correlation coefficient:

[0026]

[0027]

[0028] Where τ1 is the lag position corresponding to the maximum cross-correlation function, R gx (τ1) is the maximum cross-correlation function value, is the calculated normalized correlation coefficient, Rxx (0) and R gg (0) are the values ​​of the vibration signal x and the double damped wavelet g(t) when the lag position is 0. The normalized correlation coefficient calculated is The larger the value is, the higher the matching degree between the double damped wavelet g(t) and the internal impact component of the vibration signal x is. Therefore, the optimization goal is to search for the key parameters used to construct the wavelet with the maximum normalized correlation coefficient within the key parameter selection interval.

[0029] Use particle swarm optimization algorithm to select the key parameters of double damped wavelet. Set the fitness value to the normalized correlation coefficient The fitness value of the double damped wavelet constructed by randomly initializing the key parameter particles in the parameter selection interval is calculated, and the size of the fitness value is compared to update the individual optimal solution and the group optimal solution of the particle swarm, thereby determining the position coordinates of the next movement of the particle in the interval. When the absolute value of the difference between the optimal group fitness values ​​of two consecutive iterations of the particle is less than 1×10 -5 When it is judged to be converged and reaches the maximum value, the attenuation damping coefficient ζ1, the enhanced damping coefficient ζ2 and the natural frequency f of the last iteration d is determined as the optimal key parameter and used to construct the double-damped wavelet feature dictionary D.

[0030] S3: Convolutional basis pursuit denoising: Use the constructed double-damped wavelet feature dictionary D to perform convolutional basis pursuit denoising on the vibration signal x and construct a sparse representation model:

[0031]

[0032] Among them, D is the feature dictionary constructed in the previous step, α is the sparse coefficient that needs to be solved with the same dimension as the vibration signal x, * represents the convolution operation, y is the variable used to assist in solving the sparse coefficient, and λ is the regularization parameter with a value of 0.1. The constructed optimization problem is solved efficiently by the alternating direction multiplier method:

[0033]

[0034]

[0035] u (i+1) =u (i) +α (i+1) -y (i+1) (Formula 8)

[0036] Among them, i is the number of iterations, α (i) ,y (i) and u (i) Respectively represent the sparse coefficients obtained at the i-th iteration, the auxiliary variables and the Lagrange multipliers at the i-th iteration, and the initialized sparse coefficients α(0) , auxiliary variable y (0) and the Lagrange multiplier u (0) are all zero vectors with the same dimension as the vibration signal x, ρ is the penalty parameter, set ρ = 100λ + 0.5. The algorithm iterates until ||α (i+1) -y (i+1) ||2 is less than 1×10 -5 When it is set to converge, the last iteration solution is obtained Figure 6 The sparse coefficient α is used for subsequent spectrum amplitude modulation processing.

[0037] S4: Spectral amplitude modulation: Spectral amplitude modulation processing of the acquired sparse coefficient α:

[0038] x SAM (MO) = IFT{A(α) MO e jφ(α)} (Formula 9)

[0039] Among them, MO is the order of magnitude of the correction component, j represents the imaginary part, IFT{} is the inverse Fourier transform, A(α) represents the amplitude of the coefficient α, and Φ(α) represents the phase of the coefficient α. The order of magnitude interval is set to [-0.5, 1.5], and the interval step is set to 0.1. According to the order of magnitude interval and step size, it can be confirmed that the number of order of magnitude MOs m is 21. MO corresponds to a total of m amplitude correction components x with a dimension of 1×n SAM (MO) for subsequent order envelope spectrum characteristic analysis.

[0040] S5: Amplitude correction component feature analysis: For non-stationary operation signals, it is necessary to perform envelope analysis after angular domain resampling to identify the main feature information contained in the signal. Assume that the axis moves at a uniform speed in a short period of time during the rotation process. The rotation angle is expressed as:

[0041] β(t)=b0+b1t+b2t 2 (Formula 10)

[0042] Among them, β(t) is the rotation angle in a period of time t, b0, b1, b2 are the unknown speed related coefficients. The shaft rotation angle increment between two consecutive key phase pulse generation times is 2π, and the unknown coefficients b0, b1, b2 are solved by fitting the generation times of three consecutive key phase pulses:

[0043]

[0044] Among them, t1, t2, and t3 are the occurrence times of three consecutive key phase pulses. According to the unknown coefficients b0, b1, and b2 obtained by solving, the time corresponding to equal angle sampling is obtained:

[0045]

[0046] Among them, t β is the time corresponding to the change of the rotation angle β. According to the solved time point, the vibration signal is interpolated by using the interpolation algorithm to construct the amplitude of the vibration signal in the angular domain corresponding to the sampling time, and the signal is resampled from the time domain to the angular domain. SAM (MO) respectively performs angular domain resampling to obtain m amplitude correction components with a dimension of 1×r, and the angular domain vector x SAMA (MO). Where r represents the number of data points after the time domain amplitude correction component is resampled in the angular domain.

[0047] Set the amplitude of each order in the final order envelope spectrum to m amplitude correction components angular domain vector x SAMA The maximum value of the corresponding order amplitude in the envelope spectrum of (MO). Figure 7 It is the final order envelope spectrum analysis result obtained by the present invention. The characteristic order O of the final order envelope spectrum is compared with the theoretical fault order of the bearing. When the absolute value of the error between the main characteristic order of the pulse enhancement component and the theoretical fault order of the bearing is within 0.25, it is judged as the corresponding bearing fault type.

[0048] Outer ring fault order:

[0049]

[0050] Inner race fault order:

[0051]

[0052] Rolling element failure order:

[0053]

[0054] Among them, N is the number of rolling elements, d1 is the rolling element diameter, d2 is the pitch circle diameter, and θ is the contact angle.

Claims

1. A variable speed fault feature enhancement method based on sparse feature adaptive convolution spectrum amplitude modulation, characterized in that: The steps are as follows: S1: Vibration signal acquisition: Use an acceleration sensor to collect the variable speed bearing vibration signal, which is expressed as a time series x with a dimension of 1×n, where n is the number of sampling points and the sampling frequency of the vibration signal is 200kHz; use an incremental encoder to synchronously collect the bearing key phase pulse signal, which is a time series S with a dimension of 1×n; S2: Feature dictionary construction: The double damped wavelet g(t) is expressed as: Where t is time, τ is the time shift, Sgn() is the sign function, ω is the angular frequency, ζ1 is the attenuation damping coefficient, ζ2 is the enhancement damping coefficient, g1(t) is the attenuation section of the wavelet, g2(t) is the enhancement section of the wavelet, and f d is the natural frequency; where the attenuation damping coefficient ζ1, the enhancement damping coefficient ζ2 and the natural frequency f d It is the key parameter that affects the wavelet waveform; the pre-selected intervals of the attenuation damping coefficient ζ1 and the enhanced damping coefficient ζ2 of the double-damped wavelet are set to (0,1], and the natural frequency f d The pre-selected interval of is [500,3000], the time t interval is set to [0,0.025], and the time shift τ is set to 0.0125; The cross-correlation analysis method is used to evaluate the similarity between the double-damped wavelet and the internal impact component of the vibration signal; the cross-correlation function value is calculated by sliding the double-damped wavelet g(t) as a window on the vibration signal at different positions of the vibration signal time series x: Where τ0 is the lag position of the window sliding of the double damped wavelet g(t) on the time series x; the lag position of the maximum cross-correlation function obtained by the window sliding is used to calculate the normalized correlation coefficient: Where τ1 is the lag position corresponding to the maximum cross-correlation function, R gx (τ1) is the maximum cross-correlation function value, is the calculated normalized correlation coefficient, R xx (0) and R gg (0) are the values ​​of the vibration signal x and the double damped wavelet g(t) when the lag position is 0; the calculated normalized correlation coefficient is The larger the value is, the higher the matching degree between the double-damped wavelet g(t) and the internal impact component of the vibration signal x is; therefore, the optimization goal is to search for the key parameters for constructing the wavelet with the maximum normalized correlation coefficient within the key parameter selection interval; Use particle swarm optimization algorithm to select key parameters of double damped wavelet; set fitness value to normalized correlation coefficient The fitness value of the double damped wavelet constructed by randomly initializing the key parameter particles in the parameter selection interval is calculated, and the size of the fitness value is compared to update the individual optimal solution and the group optimal solution of the particle swarm, and then determine the position coordinates of the next movement of the particle in the interval; when the absolute value of the difference between the optimal group fitness values ​​of two consecutive iterations of the particle is less than 1×10 -5 When it is judged to be converged and reaches the maximum value, the attenuation damping coefficient ζ1, the enhanced damping coefficient ζ2 and the natural frequency f of the last iteration d It is determined as the optimal key parameter and used to construct the double damped wavelet feature dictionary D; S3: Convolutional basis pursuit denoising: Use the constructed double-damped wavelet feature dictionary D to perform convolutional basis pursuit denoising on the vibration signal x and construct a sparse representation model: Where D is the feature dictionary constructed in the previous step, α is the sparse coefficient to be solved with the same dimension as the vibration signal x, * represents the convolution operation, y is the variable used to assist in solving the sparse coefficient, and λ is the regularization parameter with a value of 0.1; the constructed optimization problem is solved efficiently by the alternating direction multiplier method: u (i+1) = u (i) + α (i+1) - y (i+1) (Equation 8) Among them, i is the number of iterations, α (i) ,y (i) and u (i) Respectively represent the sparse coefficients obtained at the i-th iteration, the auxiliary variables and the Lagrange multipliers at the i-th iteration, and the initialized sparse coefficients α (0) , auxiliary variable y (0) and the Lagrange multiplier u (0) are all zero vectors with the same dimension as the vibration signal x, ρ is the penalty parameter, set ρ = 100λ + 0.5; the algorithm iterates until ||α (i+1) -y (i+1) ||2 is less than 1×10 -5 When , it is set to converge, and the sparse coefficient α obtained by the last iteration is used for subsequent spectrum amplitude modulation processing; S4: Spectral amplitude modulation: Spectral amplitude modulation processing of the acquired sparse coefficient α: x SAM (MO) = IFT{A(α) MO e jφ(α)}} (Formula 9) Among them, MO is the order of magnitude of the correction component, j represents the imaginary part, IFT{} is the inverse Fourier transform, A(α) represents the amplitude of the coefficient α, Φ(α) represents the phase of the coefficient α; the order of magnitude interval is set to [-0.5, 1.5], and the interval step is set to 0.

1. According to the order of magnitude interval and step size, it can be confirmed that the number of order of magnitude MOs m is 21; MOs correspond to a total of m amplitude correction components x with a dimension of 1×n SAM (MO), used for subsequent order envelope spectrum characteristic analysis; S5: Amplitude correction component feature analysis: For non-stationary operation signals, it is necessary to perform envelope analysis after angular domain resampling to identify the main feature information contained in the signal; it is assumed that the axis moves at a uniform speed in a short period of time during the rotation process; the rotation angle is expressed as: β(t)=b0+b1t+b2t 2 (Formula 10) Among them, β(t) is the rotation angle within a period of time t, b0, b1, b2 are the unknown speed correlation coefficients; the shaft rotation angle increment between two consecutive key phase pulse generation times is 2π, and the unknown coefficients b0, b1, b2 are solved by fitting the generation times of three consecutive key phase pulses: Among them, t1, t2, and t3 are the occurrence times of three consecutive key phase pulses; according to the unknown coefficients b0, b1, and b2 obtained by solving, the time corresponding to equal angle sampling is obtained: Among them, t β is the time corresponding to the change of the rotation angle β; according to the solved time point, the vibration signal is interpolated by using the interpolation algorithm to construct the amplitude of the vibration signal in the angular domain corresponding to the sampling time, and the signal is resampled from the time domain to the angular domain; the m amplitude correction components x are respectively SAM (MO) respectively performs angular domain resampling to obtain m amplitude correction components with a dimension of 1×r, and the angular domain vector x SAMA (MO); where r represents the number of data points after the time domain amplitude correction component is resampled in the angular domain; Set the amplitude of each order in the final order envelope spectrum to m amplitude correction components angular domain vector x SAMA The maximum value of the corresponding order amplitude in the envelope spectrum of (MO); the characteristic order O of the final order envelope spectrum is compared with the theoretical fault order of the bearing. When the absolute value of the error between the main characteristic order of the pulse enhancement component and the theoretical fault order of the bearing is within 0.25, it is judged as the corresponding bearing fault type; Outer ring fault order: Inner race fault order: Rolling element failure order: Among them, N is the number of rolling elements, d1 is the rolling element diameter, d2 is the pitch circle diameter, and θ is the contact angle.

Images