Mechanical fault signal decomposition method based on square envelope spectrum kurtosis maximization

Abstract

The invention relates to a mechanical fault signal decomposition method based on square envelope spectrum kurtosis maximization. The method comprises the following steps: firstly, carrying out high-frequency sampling on a vibration signal by using an acceleration sensor, intercepting a segment and carrying out de-mean processing; initializing a group of filters by using a Hanning window, and filtering the preprocessed signal to obtain a plurality of decomposition modes; iteratively optimizing a filter coefficient by using a maximum characteristic value principle and taking square envelope spectrum kurtosis maximization of the filtered signal as a criterion; then, according to correlation coefficients among the modes, eliminating the modes with smaller square envelope spectrum kurtosis values in the similar signals; and finally, carrying out envelope spectrum analysis on the decomposition mode to realize fault diagnosis. According to the method, the square envelope spectrum kurtosis is used as a filtering objective function, so that the method is high in fault information indication performance and high in robustness for harmonic noise and random impact components; the decomposition result has interpretability in the physical sense, and is suitable for the condition of multi-fault mixing.

Description

Technical Field

[0001] The present invention relates to the technical field of fault diagnosis and signal processing of rotating mechanical equipment, and particularly relates to a mechanical fault signal decomposition method based on maximizing the kurtosis of the squared envelope spectrum.

Background Art

[0002] Mechanical equipment is evolving rapidly towards being large-scale, complex and integrated, which poses more severe challenges to operation status monitoring and ensuring safe production. As core components of mechanical equipment, bearings, gears, etc. are widely used in many fields such as rail transit and wind power equipment. Due to the service of these equipment under harsh working conditions, being affected by heavy loads, high speeds, overheating, etc., they are prone to failures such as wear and cracks. When large mechanical equipment fails, it will cause economic losses at least by shutting down, and at worst, it will cause production accidents and casualties. Therefore, fault diagnosis and condition monitoring of rotating mechanical equipment are of great significance for ensuring production safety.

[0003] There are various methods for condition monitoring of mechanical equipment, such as vibration analysis, oil analysis, noise monitoring, etc. Among them, vibration signals are easy to measure, can reflect fault information more accurately, and the signal analysis results are relatively reliable. Therefore, vibration analysis is one of the most effective monitoring means for fault diagnosis of mechanical equipment. Common methods include time-domain analysis, frequency-domain analysis and time-frequency domain analysis, etc. However, due to the complex structure of large mechanical equipment, the diagnostic methods based on vibration signals face many difficulties. First, the transmission path is complex. The vibration signal needs to experience a long transmission path from the rotating component to the vibration sensor fixed on the housing, which has a deteriorating effect on the fault characteristics, making it difficult to distinguish the fault information. Second, the number of mechanical components is large, and the interference caused by the movement of other components will further increase the complexity of the signal. These actual situations pose many challenges to the vibration-based fault analysis methods.

[0004] As a key means to extract multiple characteristic information from complex signals, the decomposition method provides an effective solution for mechanical fault diagnosis. The traditional signal decomposition method is based on the physical characteristics of the signal, decomposing the complex signal into simple and easy-to-analyze components. In the field of fault mechanical diagnosis, such methods can separate composite fault information, but due to the algorithm principle not targeting the characteristics of fault information and the physical meaning of the decomposition result not being clear, they are restricted in mechanical fault diagnosis.

Summary of the Invention

[0005] The present invention provides a mechanical fault signal decomposition method based on maximizing the kurtosis of the squared envelope spectrum. The kurtosis of the squared envelope spectrum is a characteristic index indicating fault components. This method maximizes it as the decomposition target and finally realizes the fault diagnosis of rolling bearings in rotating mechanical equipment based on the adaptive filtering framework.

[0006] A mechanical fault signal decomposition method based on maximizing the square envelope spectrum kurtosis, comprising the following steps:

[0007] Step 1: Fix a vibration sensor at the test point, set the sampling rate to f s , perform high-frequency sampling, truncation, and mean removal on the vibration signal, denote the processed signal as x, the length of this signal is N, set the number of decomposition modes m and the number of iterative filtering n;

[0008] Step 2: Initialize a set of filter banks evenly distributed in the frequency band with a Hanning window. First, divide the signal frequency band into K narrow bands, and its upper and lower cut-off frequencies f l and f u respectively satisfy:

[0009]

[0010] On each narrow band, initialize a filter with a length of L using a Hanning window, and set the filtering times i = 1;

[0011] Step 3: Use the filter bank in formula (1) to filter the test signal x to obtain K decomposition modes {u k , k = 1, 2,..., K}, denoted as:

[0012] u k = Xf k , (2)

[0013] In the formula, u k = [u k (1)... u k (n)... u k (N - L + 1)] T is the k-th decomposition mode, the subscript T represents the transpose operation, f k = [f k (1)... f k (l)... f k (L)] T is the k-th filter coefficient, X is the Toeplitz matrix of the signal x, and its form is as follows:

[0014]

[0015] Step 4: Take the maximization of the square envelope spectrum kurtosis (SESK) of the decomposition mode u k as the filter coefficient optimization criterion to update the filter. The matrix form of the square envelope spectrum of the mode u k is as shown in formula (3):

[0016] SES(u k ) = |FFT(|u k | 2 )| = F H |u k | 2 = F H diag(u k H )Xf k , (3)

[0017]

[0018] where FFT(·) represents the fast Fourier transform operation, F is the discrete Fourier matrix, the superscript H represents the conjugate transpose, and diag(·) represents the diagonal matrix. Based on this, the squared envelope spectrum kurtosis of mode u k is expressed as in Equation (4):

[0019]

[0020] where ||·|| 1 represents the norm operation. R represents the weighted correlation matrix:

[0021]

[0022] To optimize the filter coefficients to maximize the squared envelope spectrum kurtosis of the decomposition mode, it can be further expressed as in Equation (5):

[0023]

[0024] According to the maximum eigenvalue principle, solving the problem of the maximum average envelope spectrum kurtosis value of the decomposition mode is transformed into solving the maximum generalized eigenvalue problem of Equation (6), and the eigenvector corresponding to the maximum eigenvalue λ is the optimal filter coefficient to be optimized:

[0025] R XW1X f = R XW2X fλ, (6)

[0026] And set the iteration number i = i + 1;

[0027] Step Five: Determine whether the iteration number i reaches the preset maximum iteration number n. If it reaches, go to Step Six; if not, return to Step Three;

[0028] Step Six: Calculate the correlation coefficients between every two of the K decomposition modes using Equation (7):

[0029]

[0030] where and They are the decomposition modes u p and u q The average value of

[0031] Step 7: Select the two modes with the largest correlation coefficients, calculate the square envelope spectrum kurtosis of the two decomposition modes, then retain the mode with the larger square envelope spectrum kurtosis value, and update the number of modes K = K-1;

[0032] Step 8: Determine whether the mode number K reaches the preset reserved mode number m. If so, proceed to step 9; if not, return to step 7;

[0033] Step nine: Finally, m output decomposition modes are obtained, and envelope spectrum analysis is performed on the obtained decomposition modes to achieve composite fault diagnosis of the bearing to be tested.

[0034] Compared with the prior art, the present invention has the following advantages:

[0035] 1) The present invention uses the square envelope spectrum kurtosis as the objective function of filtering and iteratively optimizes the filter. This index has strong indicative properties for fault information and is more robust to harmonic noise and random impact components;

[0036] 2) The present invention uses frequency band distribution as the starting point for signal decomposition and the frequency spectrum characteristics of fault information as the basis for calculation, so that the decomposition result is interpretable in a physical sense.

Brief Description of the Drawings

[0037] Figure 1 The present invention is a rolling bearing fault diagnosis flow chart based on maximizing the square envelope spectrum kurtosis;

[0038] Figure 2 Schematic diagram of a train wheelset bearing test bench in an embodiment;

[0039] Figure 3 The time domain waveform of the signal after the original signal is subjected to high-frequency sampling and truncation processing in the embodiment;

[0040] Figure 4 It is the envelope spectrum of the original signal in the embodiment;

[0041] Figure 5 The decomposition mode time domain waveform obtained by using the square envelope spectrum kurtosis decomposition proposed by the present invention in the embodiment;

[0042] Figure 6 Graph showing the envelope spectrum of the decomposition mode obtained by using the square envelope spectrum kurtosis decomposition proposed in the embodiment of the present invention. [Specific implementation plan]

[0043] The present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0044] The rotating machinery equipment in the present invention takes the train wheel pair bearing as an example. As Figure 2 is a schematic structural diagram of a certain train wheel pair bearing test bench, which is composed of components such as a driving motor, a driving wheel, a loading wheel, a wheel pair bearing, and a wheel. The driving motor drives the wheel pair bearing to rotate through the driving wheel, and the loading wheel applies a load to it. The acceleration sensor is fixed at the end of the bearing to measure the vibration signal of the bearing.

[0045] Step 1: Fix the vibration sensor at the test point, set the signal sampling rate f s = 76800Hz, perform high-frequency sampling, truncation, and mean removal processing on the vibration signal, intercept a one-second signal, and denote the processed signal as x. The signal length N is 76800, and its time domain and envelope spectrum are as Figure 2 and Figure 3 shown. It can be seen that the signal-to-noise ratio of the vibration signal in the time domain diagram is relatively low, and only the outer ring fault characteristics can be seen in the envelope spectrum signal. Set the decomposition mode number m = 2 and the iterative filtering number n = 5;

[0046] Step 2: Initialize a set of filter banks evenly distributed in the frequency band with a Hanning window. First, divide the signal frequency band into 12 narrow bands evenly, and its upper and lower cut-off frequencies f l and f u respectively satisfy:

[0047]

[0048] On each narrow band, initialize a filter with a length of 40 using a Hanning window, and set the filtering number i = 1;

[0049] Step 3: Use the filter bank in formula (1) to filter the test signal x to obtain 12 decomposition modes {u k , k = 1, 2,..., 12}, denoted as:

[0050] u k = Xf k , (2)

[0051] In the formula, u k = [u k (1)... u k (n)... u k (N - L + 1)] T is the kth decomposition mode, the subscript T represents the transpose operation, f k = [f k (1)... f k (l)... f k (L)] T is the kth filter coefficient, X is the Toeplitz matrix of the signal x, and the form is as follows:

[0052]

[0053] Step 4: With the square envelop spectrum kurtosis (SESK) of the decomposition mode u k maximized as the filter coefficient optimization criterion, update the filter. The matrix form of the square envelop spectrum of mode u k is expressed as Equation (3):

[0054] SES(u k ) = |FFT(|u k | 2 )| = F H |u k | 2 = F H diag(u k H )Xf k , (3)

[0055]

[0056] where FFT(·) represents the fast Fourier transform operation, F is the discrete Fourier matrix, the superscript H represents the conjugate transpose, and diag(·) represents the diagonal matrix. Based on this, the square envelop spectrum kurtosis of mode u k is expressed as Equation (4):

[0057]

[0058] where ||·|| represents the norm operation. R represents the weighted correlation matrix:

[0059]

[0060] To optimize the filter coefficients to maximize the square envelop spectrum kurtosis of the decomposition mode, it can be further expressed as Equation (5):

[0061]

[0062] According to the maximum eigenvalue principle, the problem of solving the maximum average envelop spectrum kurtosis value of the decomposition mode is transformed into solving the maximum generalized eigenvalue problem of Equation (6), and the eigenvector corresponding to the maximum eigenvalue λ is the optimal filter coefficient to be optimized:

[0063]

[0064] And set the iteration number i = i + 1;

[0065] Step 5: Determine whether the iteration count \(i\) reaches the preset maximum iteration count of 5. If it reaches, proceed to Step 6; if not, return to Step 3;

[0066] Step 6: Calculate the correlation coefficients between every two of the 12 decomposition modes using Equation (7):

[0067]

[0068] where, and are the average values of the decomposition modes \(u\) p and \(u\) q respectively;

[0069] Step 7: Select the two modes with the largest correlation coefficient, calculate the squared envelope spectrum kurtosis of these two decomposition modes, then retain the mode with the larger squared envelope spectrum kurtosis value, and update the mode number \(K = K - 1\);

[0070] Step 8: Determine whether the mode number \(K\) reaches the preset number of retained modes of 2. If it reaches, obtain the final 2 decomposition modes and proceed to Step 9; if not, return to Step 7;

[0071] Step 9: Finally obtain 2 decomposition modes, perform envelope spectrum analysis on them, extract the fault characteristic frequencies, and finally diagnose the compound fault of the bearing to be measured.

[0072] Finally, the time-domain signals and envelope spectrum diagrams of the two decomposition modes are respectively as shown in Figure 5 and Figure 6 As shown. In the time-domain signal, the signal-to-noise ratio of the visible fault components has increased, and its periodicity and impact are more obvious; in the envelope spectrum signal, it can be seen that Mode 1 has peaks at the inner-race fault characteristic frequency \(f_i\) (60 Hz) and its multiples, and Mode 2 has peaks at the outer-race fault characteristic frequency \(f_o\) (41 Hz) and its multiples. Therefore, it can be preliminarily diagnosed that there are faults in the outer and inner races of the bearing.

Claims

1. A mechanical fault signal decomposition method based on maximizing the square envelope spectrum kurtosis, characterized in that: The following steps are involved: Step 1: Fix the vibration sensor at the test point and set the sampling rate to f s , perform high-frequency sampling, truncation and de-averaging on the vibration signal, record the processed signal as x, the length of this signal is N, set the number of decomposition modes m and the number of iterative filtering times n; Step 2: Initialize a set of filters evenly distributed in the frequency band using the Hanning window. First, divide the signal frequency band into K narrow bands with upper and lower cutoff frequencies f l and f u Respectively meet: On each narrowband, the Hanning window initializes a filter of length L and sets the number of filters i = 1; Step 3: Use the filter bank in formula (1) to filter the test signal x and obtain K decomposition patterns {u k ,k=1,2,…,K}, recorded as: u k =Xf k , (2) In the formula, u k =[u k (1) … u k (n) … u k (N-L+1)] T is the kth decomposition pattern, the subscript T indicates the transposition operation, f k =[f k (1) … f k (l) … f k (L)] T is the kth filter coefficient, X is the Toeplitz matrix of signal x, which is in the following form: Step 4: Decomposition mode u k The square envelope spectrum kurtosis (SESK) of is maximized as the filter coefficient optimization criterion to update the filter. Mode u k The matrix form of the square envelope spectrum is as follows: SES(u k )=FFT(u k 2 )|=F H |u k 2 =F H diag(u k H )Xf k , (3) Where FFT(·) represents the fast Fourier transform operation, F is the discrete Fourier matrix, the subscript H represents the conjugate transpose, and diag(·) represents the diagonal matrix. Based on this, the mode u k The square envelope spectrum kurtosis of is expressed as formula (4): Among them, ||·||1 represents the norm operation. R represents the weighted correlation matrix: To optimize the filter coefficients so that the decomposition mode square envelope spectrum kurtosis is maximized, it can be further expressed as formula (5): According to the maximum eigenvalue principle, solving the problem of the maximum average envelope spectrum kurtosis value of the decomposition mode is transformed into solving the maximum generalized eigenvalue problem of equation (6). The eigenvector corresponding to the maximum eigenvalue λ is the optimal filter coefficient to be optimized: And set the number of iterations i = i + 1; Step 5: Determine whether the number of iterations i has reached the preset maximum number of iterations n. If so, proceed to step 6; if not, return to step 3; Step 6: Use formula (7) to calculate the correlation coefficient between the K decomposition patterns: in, and They are the decomposition modes u p and u q The average value of Step 7: Select the two modes with the largest correlation coefficients, calculate the square envelope spectrum kurtosis of the two decomposition modes, then retain the mode with the larger square envelope spectrum kurtosis value, and update the number of modes K = K-1; Step 8: Determine whether the mode number K reaches the preset reserved mode number m. If so, proceed to step 9; if not, return to step 7; Step nine: Finally, m output decomposition modes are obtained, and envelope spectrum analysis is performed on the obtained decomposition modes to achieve composite fault diagnosis of the bearing to be tested.

Images