CN117760713A - Composite fault diagnosis method and system for rotating machinery based on adaptive eigenmode decomposition - Google Patents

Abstract

A rotary machine compound fault diagnosis method and system based on self-adaptive characteristic modal decomposition belong to the field of machine signal processing and equipment fault diagnosis. The invention collects the composite fault signal of the rotary machine to be diagnosed; positioning a resonance frequency band caused by faults based on an adjustable scale sliding window filtering method taking Envelope Autocorrelation Kurtosis (EAK) as an index, so as to adaptively determine the mode number of characteristic mode decomposition and fault characteristic frequencies corresponding to all modes; based on the second order cyclostationarity Index (ICS) ₂ ) The method comprises the steps of adaptively decomposing an acquired composite fault signal into a plurality of modal components for a characteristic modal decomposition method of an objective function; an envelope spectrum of the decomposition modality is calculated. The invention solves the problem of serious dependence on prior parameters such as the decomposition mode number, the fault characteristic frequency and the like, and can make a noise from strong periodicityAnd accurately positioning a resonance frequency band caused by faults in the rotating machinery composite fault signal of acoustic interference, and accurately separating and extracting each single fault component.

Description

Rotating machinery composite fault diagnosis method based on adaptive eigenmode decomposition and system

Technical field

The invention relates to the field of mechanical signal processing and equipment fault diagnosis, and in particular to a rotating machinery composite fault diagnosis method and system based on adaptive eigenmode decomposition.

Background technique

Rotating machinery is widely used in aerospace, transportation, energy and chemical industry and other fields and plays an important role. Due to long-term operating environments such as high speed, heavy load, and changing working conditions, rotating machinery components represented by bearings and gears frequently fail. Different parts of the same component are often damaged at the same time and interact with each other, resulting in compound failures. In severe cases, it will not only cause immeasurable economic losses, but even lead to casualties and serious consequences. When a compound fault occurs in a rotating machinery, the vibration signal exhibits the characteristics of simultaneous coupling of multiple periodic components and multi-pulse components. The periodic motion of other components will introduce periodic noise interference, which greatly affects the accuracy of compound fault diagnosis. Therefore, the development of accurate and effective composite fault diagnosis methods and systems for rotating machinery under strong periodic noise interference is of great significance to improve the safety and reliability of equipment.

In recent years, a large number of scholars have studied the mechanical signal decomposition theory and successfully applied it to the field of mechanical fault diagnosis. Signal decomposition theory can decompose the multi-component aliased original signal into several regular simple modes, providing a theoretical basis for the separation of mechanical composite fault characteristics. However, classic decomposition methods have their own limitations and cannot effectively decompose composite fault signals without the assistance of other methods. For example, empirical mode decomposition (EMD) has limitations of mode mixing and boundary effects, overall empirical mode decomposition (EEMD) and local mean decomposition (LMD) do not consider fault characteristics in the recursive process, and empirical wavelet transform (EWT) ) is subject to frequency band segmentation, and the decomposition effect of variational mode decomposition (VMD) heavily depends on the prior knowledge of the initial parameters. In 2022, the newly proposed eigenmode decomposition (FMD) will gradually be popularized and applied due to its unique advantages in the field of rotating machinery fault diagnosis. However, it is difficult for the FMD method to accurately estimate the true fault frequency in the presence of strong periodic interference noise, which will affect the feature decomposition effect of the FMD method. At the same time, the FMD method needs to artificially set the number of decomposition modes based on prior knowledge, and the decomposition mode number must be artificially set based on prior knowledge. The selection of state numbers also affects the decomposition effect. Therefore, in view of the limitations of existing methods, the present invention discloses a rotating machinery composite fault diagnosis method and system based on adaptive eigenmode decomposition to achieve accurate diagnosis of mechanical composite faults under strong noise.

Summary of the invention

In order to solve the shortcomings of the existing technology, especially the problems of poor adaptability and insufficient robustness under strong periodic noise interference of the eigenmode decomposition method (FMD), the present invention proposes a method based on adaptive eigenmode decomposition. The rotating machinery composite fault diagnosis method and system can adaptively determine the decomposition mode number, adaptively estimate the fault characteristic frequency, accurately separate and extract each single fault component from the rotating machinery composite fault signal interfered by strong periodic noise, and realize rotation Accurate diagnosis of mechanical composite faults.

The invention provides a rotating machinery composite fault diagnosis method based on adaptive eigenmode decomposition, which is characterized by including the following steps:

S1: Collect the rotating machinery composite fault signal x(t) to be diagnosed, t is the sampling time;

S2: Based on the adjustable-scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, the resonance frequency band caused by the fault is located, thereby adaptively determining the mode number of the characteristic mode decomposition and the fault characteristics corresponding to each mode. frequency;

S3: Based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS ₂ ) as the objective function, the collected composite fault signal is adaptively decomposed into several modal components;

S4: Calculate the envelope spectrum of the decomposed modes to achieve accurate extraction of rotating machinery fault features and fault diagnosis;

In one embodiment of the present invention, in step S2, the adjustable-scale sliding window filtering method based on envelope autocorrelation kurtosis (EAK) as an indicator includes the following steps:

S201: Load the original signal x(t) and start iteration i=1;

S202: Construct a scale-adaptive bandpass filter based on the second-order Butterworth filter. Assume that the filter bandwidth is α, which is named as the scale factor of the filter. Set the initial α = 0.01f _s , f _s is the sampling frequency of the original signal, and set the window sliding step size to 0.2α. Construct a bandpass filter bank and perform sliding window filtering on the original signal x(t) to obtain a set of filtered sub-signals {x ₁ , x ₂ ,..., x _END };

S203: Screen out the filtered sub-signal x _max with the largest envelope autocorrelation kurtosis (EAK) value, that is, the sub-signal with the most fault information content. The calculation formula of EAK is as follows:

X(t)=hilbert(x(t))

R _XX (τ)＝∫X(t)X(t+τ)dt

In the _formula , X(t) represents the hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; represents the time delay; EAK(x(t)) represents the EAK value of the signal x(t), N represents the number of sampling points of the signal x(t), min(·) represents the minimum value function in the MATLAB toolkit, ∫· represents the integral Symbol, ∑· represents the summation symbol, R _XX (ii) represents the ii-th value of R _XX (τ), ii=1, 2,..., N/2;

S204: Adaptively adjust the scale factor α value of the sliding window corresponding to x _max so that it covers the entire resonance frequency band excited by the fault as much as possible, and then obtain the i-th optimal filtered signal x _{best_i} by updating x _max . The α adaptive adjustment method is as follows:

Update α = α + 10% α, that is, gradually increase the filter bandwidth, and the center frequency of the sliding window corresponding to x _max remains unchanged. Use the filter after the scale update to filter x(t) to obtain a new filtered signal x′ _max . If EAK(x′ _max )>EAK(x _max ), update x _max =x′ _max ; otherwise, stop updating α, and obtain the filtered sub-signal x _{best_i} with the optimal scale;

S205: Calculate the envelope harmonic product spectrum (EHPS) of the optimal filtered sub-signal x _{best_i} , and select the frequency corresponding to the position of the global maximum value v _i of the EHPS as the estimated i-th fault characteristic frequency _fi . The calculation formula of EHPS is as follows:

In the formula, P(w) represents EHPS, F(w) represents the Fourier transform amplitude of the signal envelope, m=1, 2,..., M, M represents the number of harmonics considered, Π· represents continuous multiplication operator symbol;

S206: If i>1, calculate the EHPS reduction rate r _i of x _{best_i} , otherwise return to step S202 and update x(t)=x(t)-x _{best_i} , i=i+1, x represents the original signal. The calculation formula of r _i is as follows:

S207: When i>3, if r _i-2 > r _i-3 and r _i-2 > r _i-1 , let index=i-2, freq=[f ₁ , f ₂ ,..., f _index ], obtain the decomposed target mode number K = index, and the estimated fault characteristic frequency required for each mode decomposition is freq = [f ₁ , f ₂ , ..., f _index ]; otherwise return to step S202, Update x(t)=x(t)-x _{best_i} , i=i+1. When i≤3, return to step S202 and update x(t)=x(t)-x _{best_i} , i=i+1;

In one embodiment of the present invention, in step S3, the eigenmode decomposition method with the second-order cyclostationarity index (ICS ₂ ) as the objective function includes the following steps:

S301: Load the original signal x=x(t), and use the target mode number K and the fault characteristic frequency freq obtained in step S207 as the input parameters of the eigenmode decomposition;

S302: Set the filter length L, use Q Hanning windows to initialize the FIR filter bank, and obtain the initial FIR filter bank {h ₁ , h ₂ , ..., h _Q }. The upper and lower cutoff frequencies h _l and _hu of the FIR filter can be expressed as:

In the formula, f _s is the sampling frequency of the original signal, q is the sequence number of the filter, q=1, 2,..., Q;

S303: Use the deconvolution method to filter the original signal x to obtain the decomposed mode, and maximize the constrained decomposition by maximizing the second-order cyclostationarity (ICS ₂ ) of the decomposed mode. The expression of ICS ₂ is as follows:

In the formula, u _q represents the decomposition mode, r is the sample index, β represents the number of sample points corresponding to the fault characteristic frequency, L is the filter length, N represents the number of sampling points of the original signal, u _q [n] represents the third of the matrix u _q n elements, n=1, 2,..., N-L+1, j represents the complex unit, |·| represents the absolute value operation symbol, ||·|| represents the norm operation symbol;

The expression of ICS ₂ can be rewritten into the following matrix form:

E＝[e ₁ ,..., e _r ,..., e _R ]

e _r =[e ^{-j2πrβ(L-1)} ,...,e ^{-j2πrβ(N-1)} ] ^H

In the formula, D represents the diagonal matrix of u _q , u _q [n] represents the n-th element of the matrix u _q , n=1, 2,..., N-L+1, diag(·) represents the MATLAB tool The diagonal matrix function in the package; E represents a matrix composed of R matrices e _r sequentially spliced horizontally, r is the sample index, r=1, 2,..., R, R is a positive integer, the recommended value range is 50 to 200, the superscript H of the matrix represents the conjugate transpose of the matrix, and j represents the complex unit;

Therefore, the decomposition process of the eigenmode decomposition method based on ICS ₂ is equivalent to finding the solution to the following constrained problem:

In the formula, h _q [l] is the q-th filter, its length is L, L=1,2,…,L; Represents the filter h _q [l] when solving to maximize ICS ₂ ;

The deconvolution method is used to iteratively filter to obtain the decomposition mode. The iterative process of the deconvolution method is as follows:

h _q =Xh _q

In the formula, u _q represents the decomposition mode, h _q represents the q-th filter, q=1,2,...,Q, and X represents the shape of the original signal after equal length interception is (N-L+1)×L matrix, L is the filter length, N represents the number of sampling points of the original signal;

Then, the decomposed modal ICS ₂ can be defined as:

In the formula, R _XWX and R _XX represent the weighted correlation matrix and the correlation matrix respectively, and W is the intermediary matrix that controls the weighted correlation matrix;

Calculate the ICS ₂ difference of the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than eta, the algorithm is deemed to have converged and the iteration will be stopped, that is:

Among them, eta is the convergence condition, the recommended value of eta is 10 ^-6 , and y is the current number of iterations;

In order to obtain the optimal filter when maximizing ICS ₂ , it is equivalent to solving the eigenvector related to the maximum eigenvalue λ of the following eigenvalue problem:

R _XWX h _q =R _XX h _q λ

In the formula, λ can be obtained by solving the above formula, λ represents the maximum eigenvalue, that is, the optimal ICS ₂ ;

Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposition mode u _q ;

S304: Mode selection. Calculate the ICS ₂ of the decomposed mode u _q , select the mode with the largest ICS ₂ as the k-th final decomposed mode _uk , and then update the original signal x= _xuk , when the number of decomposed modes reaches the target mode obtained in step S207 The decomposition ends when the state number is K.

Another object of the present invention is to provide a rotating machinery composite fault diagnosis system based on adaptive eigenmode decomposition, which is characterized by including:

The signal acquisition module is used to collect the rotating machinery composite fault signal x(t) to be diagnosed, where t is the sampling time;

The parameter determination module is used to locate the resonance frequency band caused by the fault based on the adjustable-scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, thereby adaptively determining the mode number and each mode of the characteristic mode decomposition. Corresponding fault characteristic frequency;

The signal decomposition module is used to adaptively decompose the collected composite fault signal into several modal components based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS ₂ ) as the objective function;

The fault diagnosis module is used to calculate the envelope spectrum of the decomposed modes to achieve accurate extraction of rotating machinery fault features and fault diagnosis;

In one embodiment of the present invention, the parameter determination module includes:

The sliding window filtering unit is used to construct a bandpass filter bank and perform sliding window filtering on the original signal x(t) to obtain a set of filtered sub-signals {x ₁ , x ₂ ,..., x _END };

The target screening unit is used to screen out the filtered sub-signal x _max with the largest envelope autocorrelation kurtosis (EAK) value, that is, the sub-signal with the most fault information content. The calculation formula of EAK is as follows:

X(t)=hilbert(x(t))

R _XX (τ)＝∫X(t)X(t+τ)dt

In the _formula , X(t) represents the hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; represents the time delay; EAK(x(t)) represents the EAK value of the signal x(t), N represents the number of sampling points of the signal x(t), min(·) represents the minimum value function in the MATLAB toolkit, ∫· represents the integral Symbol, ∑· represents the summation symbol, R _XX (ii) represents the ii-th value of R _XX (τ), ii=1, 2,..., N/2;

The filter scale _adjustment unit is used to adaptively adjust the scale factor _α value of the sliding window corresponding to Signal x _{best_i} . The α adaptive adjustment method is as follows:

Update α = α + 10% α, that is, gradually increase the filter bandwidth, and the center frequency of the sliding window corresponding to x _max remains unchanged. Use the filter after the scale update to filter x(t) to obtain a new filtered signal x′ _max . If EAK(x′ _max )>EAK(x _max ), update x _max =x′ _max ; otherwise, stop updating α, and obtain the filtered sub-signal x _{best_i} with the optimal scale;

The frequency estimation unit is used to calculate the EHPS of the optimal filtered sub-signal x _{best_i} , and select the frequency corresponding to the position of the EHPS global maximum value _vi as the estimated i-th fault characteristic frequency _fi . The calculation formula of EHPS is as follows:

In the formula, P(w) represents EHPS, F(w) represents the Fourier transform amplitude of the signal envelope, m=1, 2,..., M, M represents the number of harmonics considered, Π· represents continuous multiplication operator symbol;

The target mode number determination unit is used to locate the number of resonance frequency bands in the signal, thereby determining the number of modes that need to be decomposed. If i>1, calculate the EHPS reduction rate r _i of x _{best_i} , otherwise return to the sliding window filter unit and update x(t)=x(t)-x _{best_i} , i=i+1, x represents the original signal. The calculation formula of r _i is as follows:

When i>3, if r _i-2 > r _i-3 and r _i-2 > r _i-1 , let index=i-2, freq=[f ₁ , f ₂ ,..., f _index ] , obtain the decomposed target mode number K=index, and the estimated fault characteristic frequency required for each mode decomposition is freq=[f ₁ , f ₂ ,..., f _index ]; otherwise return to the sliding window filter unit, Update x(t)=x(t)-x _{best_i} , i=i+1. When i≤3, the sliding window filtering unit is also returned and x(t)=x(t)-x _{best_i} is updated, i=i+1;

In one embodiment of the present invention, the signal decomposition module includes:

Filter initialization unit, used for FIR filter initialization. Set the filter length L, use Q Hanning windows to initialize the FIR filter bank, and obtain the initial FIR filter bank {h ₁ , h ₂ ,..., h _Q }. The upper and lower cutoff frequencies h _l and _hu of the FIR filter can be expressed as:

In the formula, f _s is the sampling frequency of the original signal, q is the sequence number of the filter, q=1, 2,..., Q;

Objective function unit, used to establish the objective function ICS ₂ of the eigenmode decomposition method:

In the formula, u _q represents the decomposition mode, r is the sample index, β represents the number of sample points corresponding to the fault characteristic frequency, L is the filter length, N represents the number of sampling points of the original signal, u _q [n] represents the third of the matrix u _q n elements, n=1,2,…,N-L+1, j represents the complex unit, |·| represents the absolute value operation symbol, ||·|| represents the norm operation symbol;

The expression of ICS ₂ can be rewritten into the following matrix form:

E＝[e ₁ ,…, _er ,…,e _R ]

e _r =[e ^{-j2πrβ(L-1)} ,…,e ^-j2πrβ(N-1 )] ^H

In the formula, D represents the diagonal matrix of u _q , u _q [n] represents the nth element of the matrix u _q , n=1,2,…,N-L+1, diag(·) represents the MATLAB tool package The diagonal matrix function _of , the superscript H of the matrix represents the conjugate transpose of the matrix, and j represents the complex unit;

Therefore, the decomposition process of the eigenmode decomposition method based on ICS ₂ is equivalent to finding the solution to the following constrained problem:

In the formula, h _q [l] is the q-th filter, its length is L, l=1,2,…,L; Represents the filter h _q [l] when solving to maximize ICS ₂ ;

The deconvolution iteration unit is used to iteratively filter using the deconvolution method to obtain the decomposition mode. The iteration process of the deconvolution method is as follows:

u _q =Xh _q

In the formula, u _q represents the decomposition mode, h _q represents the q-th filter, q=1, 2,..., Q, X represents the shape of the original signal after equal length interception is (N-L+1) ×L matrix, L is the filter length, and N represents the number of sampling points of the original signal;

Then, the decomposed modal ICS ₂ can be defined as:

In the formula, R _XWX and R _XX represent the weighted correlation matrix and the correlation matrix respectively, and W is the intermediary matrix that controls the weighted correlation matrix;

Calculate the ICS ₂ difference of the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than eta, the algorithm is deemed to have converged and the iteration will be stopped, that is:

Among them, eta is the convergence condition, the recommended value of eta is 10 ^-6 , and y is the current number of iterations;

In order to obtain the optimal filter when maximizing ICS ₂ , it is equivalent to solving the eigenvector related to the maximum eigenvalue λ of the following eigenvalue problem:

R _XWX h _q =R _XX h _q λ

In the formula, λ can be obtained by solving the above formula, λ represents the maximum eigenvalue, that is, the optimal ICS ₂ ;

Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposition mode u _q ;

The mode selection unit is used to select the optimal mode. Calculate the ICS ₂ of the decomposed mode u _q , select the mode with the largest ICS ₂ as the k-th final decomposed mode _uk , then update the original signal x=xu _k , and end when the number of decomposed modes reaches the target mode number K break down.

This invention aims at the characteristics of the rotating machinery composite fault signal showing the simultaneous coupling of multiple periodic components and multi-pulse components. In order to accurately separate and extract each single fault component from the rotating machinery composite fault signal interfered by strong periodic noise, a proposed method A composite fault diagnosis method and system for rotating machinery based on adaptive eigenmode decomposition. Through effective adaptive estimation of fault characteristic frequency and adaptive optimization strategy of decomposition mode number, the resonance frequency band caused by the fault is accurately located and the modal decomposition is obtained at the same time. The main prior knowledge required; further, based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS ₂ ) as the objective function, the accurate decomposition of multiple fault components under strong periodic noise interference is achieved. The present invention better solves the problems of poor adaptability and poor effect under strong periodic noise interference in existing related methods, and has broad application prospects and promotion significance in the field of rotating machinery composite fault diagnosis.

The beneficial technical effects of the present invention are:

1. The decomposition mode number adaptive optimization and fault frequency adaptive estimation strategy proposed by the present invention can accurately locate the resonance frequency band caused by rotating machinery faults under strong noise through the sliding window filter guided by the EAK index, thereby further obtaining the required decomposition information. Decompose the mode numbers and their corresponding fault frequencies. This strategy can be used as a preprocessing step for methods such as maximum correlation kurtosis deconvolution (MCKD), improved multi-point optimal minimum entropy deconvolution (MOMEDA), maximum second-order cyclic stationary blind deconvolution (CYCBD), VMD, etc. to improve the effectiveness of these methods.

2. The ICS _2- based eigenmode decomposition method proposed by the present invention can accurately extract weak fault information of rotating machinery under a strong noise background. Compared with the correlation kurtosis (CK) index, using ICS ₂ as the deconvolution objective function can better extract the characteristic components of rotating machinery faults under strong noise. ICS ₂ has better rotating machinery fault information sensitivity and noise robustness. sex.

3. The present invention has good adaptability and accuracy for single and compound fault feature extraction under strong periodic noise interference conditions, and has outstanding feature decomposition effect especially in the field of compound fault diagnosis.

Description of the drawings

In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings needed to be used in the description of the embodiments will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can also be obtained based on these drawings without exerting creative efforts.

Figure 1 is a schematic diagram of a rotating machinery composite fault diagnosis method and system based on adaptive eigenmode decomposition provided by an embodiment of the present invention;

Figure 2 is a detailed flow chart of a rotating machinery composite fault diagnosis method based on adaptive eigenmode decomposition provided by an embodiment of the present invention;

Figure 3 is a time domain diagram of a composite fault signal provided by an embodiment of the present invention;

Figure 4 is a spectrum diagram of a composite fault signal provided by an embodiment of the present invention;

Figure 5 is an adaptive decomposition mode of a composite fault signal of a rotating machinery composite fault diagnosis method based on adaptive eigenmode decomposition provided by an embodiment of the present invention;

Figure 6 is an envelope spectrum of the adaptive decomposition mode of the composite fault signal of a rotating machinery composite fault diagnosis method based on adaptive eigenmode decomposition provided by an embodiment of the present invention;

Figure 7 is a result diagram of signal processing using the VMD method in the embodiment of the present invention;

Figure 8 is a diagram showing the results of signal processing using the FMD method in the embodiment of the present invention.

Detailed ways

In order to make the purpose, technical solutions and advantages of the present invention clearer, the embodiments of the present invention will be described in further detail below with reference to the accompanying drawings.

Figure 1 is a schematic diagram of a rotating machinery composite fault diagnosis method and system based on adaptive eigenmode decomposition provided by an embodiment of the present invention. The detailed process is shown in Figure 2. Referring to Figure 1, the rotating machinery composite fault diagnosis method based on adaptive eigenmode decomposition includes the following steps:

S1: Collect the rotating machinery composite fault signal x(t) to be diagnosed, t is the sampling time;

S2: Based on the adjustable-scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, the resonance frequency band caused by the fault is located, thereby adaptively determining the mode number of the characteristic mode decomposition and the fault characteristics corresponding to each mode. frequency;

S3: Based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS ₂ ) as the objective function, the collected composite fault signal is adaptively decomposed into several modal components;

S4: Calculate the envelope spectrum of the decomposed modes to achieve accurate extraction of rotating machinery fault features and fault diagnosis.

In one embodiment of the present invention, in step S2, the adjustable-scale sliding window filtering method based on envelope autocorrelation kurtosis (EAK) as an indicator includes the following steps:

S201: Load the original signal x(t) and start iteration i=1;

S202: Construct a scale-adaptive bandpass filter based on the second-order Butterworth filter. Assume that the filter bandwidth is α, which is named as the scale factor of the filter. Set the initial α = 0.01f _s , f _s is the sampling frequency of the original signal, and set the window sliding step size to 0.2α. Construct a bandpass filter bank and perform sliding window filtering on the original signal x(t) to obtain a set of filtered sub-signals {x ₁ , x ₂ ,..., x _END };

S203: Screen out the filtered sub-signal x _max with the largest envelope autocorrelation kurtosis (EAK) value, that is, the sub-signal with the most fault information content. The calculation formula of EAK is as follows:

X(t)=hilbert(x(t))

R _XX (τ)＝∫X(t)X(t+τ)dt

In the _formula , X(t) represents the hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; represents the time delay; EAK(x(t)) represents the EAK value of the signal x(t), N represents the number of sampling points of the signal x(t), min(·) represents the minimum value function in the MATLAB toolkit, ∫· represents the integral Symbol, ∑· represents the summation symbol, R _XX (ii) represents the ii-th value of R _XX (τ), ii=1, 2,..., N/2;

S204: Adaptively adjust the scale factor α value of the sliding window corresponding to x _max so that it covers the entire resonance frequency band excited by the fault as much as possible, and then obtain the i-th optimal filtered signal x _{best_i} by updating x _max . The α adaptive adjustment method is as follows:

Update α = α + 10% α, that is, gradually increase the filter bandwidth, and the center frequency of the sliding window corresponding to x _max remains unchanged. Use the filter after the scale update to filter x(t) to obtain a new filtered signal x′ _max . If EAK(x′ _max )>EAK(x _max ), update x _max =x′ _max ; otherwise, stop updating α, and obtain the filtered sub-signal x _{best_i} with the optimal scale;

S205: Calculate the EHPS of the optimal filtered sub-signal x _{best_i} , and select the frequency corresponding to the position of the EHPS global maximum value _vi as the estimated i-th fault characteristic frequency _fi . The calculation formula of EHPS is as follows:

In the formula, P(w) represents EHPS, F(w) represents the Fourier transform amplitude of the signal envelope, m=1, 2,..., M, M represents the number of harmonics considered, Π· represents continuous multiplication operator symbol;

It should be noted that the number of harmonics M considered can be adjusted according to the signal characteristics of the actual diagnostic equipment. Preferably, M=5;

S206: If i>1, calculate the EHPS reduction rate r _i of x _{best_i} , otherwise return to step S202 and update x(t)=x(t)-x _{best_i} , i=i+1, x represents the original signal. The calculation formula of r _i is as follows:

S207: When i>3, if r _i-2 > r _i-3 and r _i-2 > r _i-1 , let index=i-2, freq=[f ₁ , f ₂ ,..., f _index ], obtain the decomposed target mode number K = index, and the estimated fault characteristic frequency required for each mode decomposition is freq = [f ₁ , f ₂ , ..., f _index ]; otherwise return to step S202, Update x(t)=x(t)-x _{best_i} , i=i+1. When i≤3, return to step S202 and update x(t)=x(t)-x _{best_i} , i=i+1:

In one embodiment of the present invention, in step S3, the eigenmode decomposition method with the second-order cyclostationarity index (ICS ₂ ) as the objective function includes the following steps:

S301: Load the original signal x=x(t), and use the target mode number K and the fault characteristic frequency freq obtained in step S207 as the input parameters of the eigenmode decomposition;

S302: Set the filter length L, use Q Hanning windows to initialize the FIR filter bank, and obtain the initial FIR filter bank {h ₁ , h ₂ , ..., h _Q }. The upper and lower cutoff frequencies h _l and _hu of the FIR filter can be expressed as:

In the formula, f _s is the sampling frequency of the original signal, q is the sequence number of the filter, q=1, 2,..., Q;

It should be noted that the filter length L can be determined according to the size of the computing task. Preferably, L=100 is set;

S303: Use the deconvolution method to filter the original signal x to obtain the decomposed mode, and maximize the constrained decomposition by maximizing the second-order cyclostationarity (ICS ₂ ) of the decomposed mode. The expression of ICS ₂ is as follows:

In the formula, u _q represents the decomposition mode, r is the sample index, β represents the number of sample points corresponding to the fault characteristic frequency, L is the filter length, N represents the number of sampling points of the original signal, u _q [n] represents the third of the matrix u _q n elements, n=1, 2,..., N-L+1, j represents the complex unit, |·| represents the absolute value operation symbol, ||·|| represents the norm operation symbol;

The expression of ICS ₂ can be rewritten into the following matrix form:

E＝[e ₁ ,…, _er ,…,e _R ]

e _r =[e ^{-j2πrβ(L-1)} ,…,e ^-j2πrβ(N-1 )] ^H

In the formula, D represents the diagonal matrix of u _q , u _q [n] represents the nth element of the matrix u _q , n=1,2,…,N-L+1, diag(·) represents the MATLAB tool package The diagonal matrix function _of , the superscript H of the matrix represents the conjugate transpose of the matrix, and j represents the complex unit;

Therefore, the decomposition process of the eigenmode decomposition method based on ICS ₂ is equivalent to finding the solution to the following constrained problem:

In the formula, h _q [l] is the q-th filter, its length is L, L=1,2,…,L; Represents the filter h _q [l] when solving to maximize ICS ₂ ;

Calculate the ICS ₂ difference of the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than eta, the algorithm is deemed to have converged and the iteration will be stopped, that is:

Among them, eta is the convergence condition, the recommended value of eta is 10 ^-6 , and y is the current number of iterations;

The deconvolution method is used to iteratively filter to obtain the decomposition mode. The iterative process of the deconvolution method is as follows:

u _q =Xh _q

In the formula, u _q represents the decomposition mode, h _q represents the q-th filter, q=1,2,...,Q, and X represents the shape of the original signal after equal length interception is (N-L+1)×L matrix, L is the filter length, N represents the number of sampling points of the original signal;

Then, the decomposed modal ICS ₂ can be defined as:

In the formula, R _XWX and R _XX represent the weighted correlation matrix and the correlation matrix respectively, and W is the intermediate variable matrix that controls the weighted correlation matrix;

In order to obtain the optimal filter when maximizing ICS ₂ , it is equivalent to solving the eigenvector related to the maximum eigenvalue λ of the following eigenvalue problem:

R _XWX h _q =R _XX h _q λ

In the formula, λ can be obtained by solving the above formula, λ represents the maximum eigenvalue, that is, the optimal ICS ₂ ;

Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposition mode u _q ;

S304: Mode selection. Calculate the ICS ₂ of the decomposed mode u _q , select the mode with the largest ICS ₂ as the k-th final decomposed mode _uk , and then update the original signal x= _xuk , when the number of decomposed modes reaches the target mode obtained in step S207 The decomposition ends when the state number is K.

A rotating machinery composite fault diagnosis system based on adaptive eigenmode decomposition according to the second aspect of the present invention is characterized by including:

The signal acquisition module is used to collect the rotating machinery composite fault signal x(t) to be diagnosed, where t is the sampling time;

The parameter determination module is used to locate the resonance frequency band caused by the fault based on the adjustable-scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, thereby adaptively determining the mode number and each mode of the characteristic mode decomposition. Corresponding fault characteristic frequency;

The signal decomposition module is used to adaptively decompose the collected composite fault signal into several modal components based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS ₂ ) as the objective function;

The fault diagnosis module is used to calculate the envelope spectrum of the decomposed modes to achieve accurate extraction of rotating machinery fault features and fault diagnosis;

In one embodiment of the present invention, the parameter determination module includes:

The sliding window filtering unit is used to construct a bandpass filter bank and perform sliding window filtering on the original signal x(t) to obtain a set of filtered sub-signals {x ₁ , x ₂ ,..., x _END };

The target screening unit is used to screen out the filtered sub-signal x _max with the largest envelope autocorrelation kurtosis (EAK) value, that is, the sub-signal with the most fault information content. The calculation formula of EAK is as follows:

X(t)=hilbert(x(t))

R _XX (τ)＝∫X(t)X(t+τ)dt

In the _formula , X(t) represents the hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; represents the time delay; EAK(x(t)) represents the EAK value of the signal x(t), N represents the number of sampling points of the signal x(t), min(·) represents the minimum value function in the MATLAB toolkit, ∫· represents the integral Symbol, ∑· represents the summation symbol, R _XX (ii) represents the ii-th value of R _XX (τ), ii=1, 2,..., N/2;

The filter scale _adjustment unit is used to adaptively adjust the scale factor _α value of the sliding window corresponding to Signal x _{best_i} . The α adaptive adjustment method is as follows:

Update α = α + 10% α, that is, gradually increase the filter bandwidth, and the center frequency of the sliding window corresponding to x _max remains unchanged. Use the filter after the scale update to filter x(t) to obtain a new filtered signal x′ _max . If EAK(x′ _max )>EAK(x _max ), update x _max =x′ _max ; otherwise, stop updating α, and obtain the filtered sub-signal x _{best_i} with the optimal scale;

The frequency estimation unit is used to calculate the EHPS of the optimal filtered sub-signal x _{best_i} , and select the frequency corresponding to the position of the EHPS global maximum value _vi as the estimated i-th fault characteristic frequency _fi . The calculation formula of EHPS is as follows:

In the formula, P(w) represents EHPS, F(w) represents the Fourier transform amplitude of the signal envelope, m=1, 2,..., M, M represents the number of harmonics considered, Π· represents continuous multiplication operator symbol;

The target mode number determination unit is used to locate the number of resonance frequency bands in the signal, thereby determining the number of modes that need to be decomposed. If i>1, calculate the EHPS reduction rate r _i of x _{best_i} , otherwise return to the sliding window filter unit and update x(t)=x(t)-x _{best_i} , i=i+1, x represents the original signal. The calculation formula of r _i is as follows:

When i>3, if r _i-2 > r _i-3 and r _i-2 > r _i-1 , let index=i-2, freq=[f ₁ , f ₂ ,..., f _index ] , obtain the decomposed target mode number K=index, and the estimated fault characteristic frequency required for each mode decomposition is freq=[f ₁ , f ₂ ,..., f _index ]; otherwise return to the sliding window filter unit, Update x(t)=x(t)-x _{best_i} , i=i+1. When i≤3, the sliding window filtering unit is also returned and x(t)=x(t)-x _{best_i} is updated, i=i+1;

In one embodiment of the present invention, the signal decomposition module includes:

Filter initialization unit, used for FIR filter initialization. Set the filter length L, use Q Hanning windows to initialize the FIR filter bank, and obtain the initial FIR filter bank {h ₁ , h ₂ ,..., h _Q }. The upper and lower cutoff frequencies h _l and _hu of the FIR filter can be expressed as:

In the formula, f _s is the sampling frequency of the original signal, q is the sequence number of the filter, q=1, 2,..., Q;

Objective function unit, used to establish the objective function ICS ₂ of the eigenmode decomposition method:

In the formula, u _q represents the decomposition mode, r is the sample index, β represents the number of sample points corresponding to the fault characteristic frequency, L is the filter length, N represents the number of sampling points of the original signal, u _q [n] represents the third of the matrix u _q n elements, n=1, 2,..., N-L+1, j represents the complex unit, |·| represents the absolute value operation symbol, ||·|| represents the norm operation symbol;

The expression of ICS ₂ can be rewritten into the following matrix form:

E＝[e ₁ , .., e _r , ..., e _R ]

e _r =[e ^{-j2πrβ(L-1)} ,...,e ^{-j2πrβ(N-1)} ] ^H

In the formula, D represents the diagonal matrix of u _q , u _q [n] represents the n-th element of the matrix u _q , n=1, 2,..., N-L+1, diag(·) represents the MATLAB tool The diagonal matrix function in the package; E represents a matrix composed of R matrices e _r that are spliced horizontally in sequence, r is the sample index, r=1,2,...,R, R is a positive integer, and the recommended value range is 50 ~200, the superscript H of the matrix represents the conjugate transpose of the matrix, and j represents the complex unit;

Therefore, the decomposition process of the eigenmode decomposition method based on ICS ₂ is equivalent to finding the solution to the following constrained problem:

In the formula, h _q [l] is the q-th filter, its length is L, l=1,2,…,L; Represents the filter h _q [l] when solving to maximize ICS ₂ ;

The deconvolution iteration unit is used to iteratively filter using the deconvolution method to obtain the decomposition mode. The iteration process of the deconvolution method is as follows:

u _q =Xh _q

In the formula, u _q represents the decomposition mode, h _q represents the q-th filter, q=1,2,...,Q, and X represents the shape of the original signal after equal length interception is (N-L+1)×L matrix, L is the filter length, N represents the number of sampling points of the original signal;

Then, the decomposed modal ICS ₂ can be defined as:

In the formula, R _XWX and R _XX represent the weighted correlation matrix and the correlation matrix respectively, and W is the intermediary matrix that controls the weighted correlation matrix;

Calculate the ICS ₂ difference of the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than eta, the algorithm is deemed to have converged and the iteration will be stopped, that is:

Among them, eta is the convergence condition, the recommended value of eta is 10 ^-6 , and y is the current number of iterations;

In order to obtain the optimal filter when maximizing ICS ₂ , it is equivalent to solving the eigenvector related to the maximum eigenvalue λ of the following eigenvalue problem:

R _XWX h _q =R _XX h _q λ

In the formula, λ can be obtained by solving the above formula, λ represents the maximum eigenvalue, that is, the optimal ICS ₂ ;

Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposition mode u _q ;

The mode selection unit is used to select the optimal mode. Calculate the ICS ₂ of the decomposed mode u _q , select the mode with the largest ICS ₂ as the k-th final decomposed mode _uk , then update the original signal x=xu _k , and end the decomposition when the number of decomposed modes reaches the target number of modes .

The specific embodiments are as follows:

Taking the complex mechanical transmission system test bench as an example, the bearing is a 30310 model tapered roller bearing. The 30310 bearing parameters and theoretical fault characteristic frequency are shown in Table 1. Before the experiment, wire cutting processing was used to obtain a rolling bearing with composite faults in the inner and outer rings. The fault injection size was 1 mm in width and 0.5 mm in depth. Under the working conditions of 1500r/min rotation speed and 60N·m load, vibration data were collected at a sampling frequency of 25.6kHz. Figure 3 is the time domain diagram of the vibration signal of the composite fault bearing, and Figure 4 is the spectrum diagram of the vibration signal of the composite fault. It can be seen from the figure that the time domain signal contains a large noise component, and the impact caused by the fault cannot be seen with the naked eye; a strong periodic noise component can be seen in the spectrum chart.

Table 130310 bearing parameters and theoretical fault characteristic frequency

The collected composite fault bearing signals are analyzed using the method of the present invention. First, set the sampling frequency f _s = 25.6 kHz, and the number of harmonics considered M = 5. Based on the adjustable-scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, the resonance frequency band caused by the fault is located, so as to automatically It is adapted to determine the mode number of eigenmode decomposition and the fault characteristic frequency corresponding to each mode. In this embodiment, the adaptively obtained decomposition target mode number K=2, and the obtained fault characteristic frequency estimate is freq=[163.87Hz, 235.43Hz]. Then, set the filter length L = 100, η = 10 ^-6 , and adaptively decompose the collected composite fault signal into several modes based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS ₂ ) as the objective function. Portion. In this embodiment, the adaptive decomposition mode of the obtained composite fault signal is shown in Figure 5. Finally, the envelope spectrum of the decomposed modes is calculated to achieve accurate extraction of rotating machinery fault features and fault diagnosis. In this embodiment, the adaptive decomposition modal envelope spectrum of the composite fault signal is shown in Figure 6. It can be seen from Figure 6 that the present invention extracts two fault components from the composite fault signal. The envelope spectrum of the first decomposed mode clearly shows the outer ring fault frequency and its 2 to 7 multiples. The second decomposed mode The inner ring fault frequency and its 2-fold frequency can be seen in the envelope spectrum of the state, so it can be diagnosed that a composite fault of the inner and outer rings has occurred, which is consistent with the preset fault type.

In order to further verify the superiority of the present invention, variational mode decomposition (VMD) and eigenmode decomposition (FMD) were used to test the same signals respectively. The decomposed mode and its envelope spectrum obtained by the VMD method are as shown in the figure As shown in Figure 7; the decomposed mode and its envelope spectrum obtained by the FMD method are shown in Figure 8. It can be seen from Figures 7 and 8 that neither the VMD nor the FMD method can effectively achieve the decomposition and extraction of composite fault signal features in this embodiment. Therefore, the present invention has advantages in feature extraction of composite fault signals of rotating machinery under strong periodic background noise.

The above is only a description of the preferred embodiments of the present application and the technical principles used. Those skilled in the art should understand that the scope involved in the present application is not limited to technical solutions formed by specific combinations of the above technical features. , it should also cover other technical solutions formed by any combination of the above technical features or their equivalent features without departing from the inventive concept. For example, the above features have similar functions to those disclosed in this application (but are not limited to). A technical solution formed by replacing technical features with each other.

Except for the technical features described in the description, the remaining technical features are technologies known to those skilled in the art. In order to highlight the innovative features of the present invention, the remaining technical features will not be described in detail here.

Claims (2)

1. A composite fault diagnosis method for rotating machinery based on adaptive eigenmode decomposition, which is characterized by including the following steps: S1: Collect the rotating machinery composite fault signal x(t) to be diagnosed, t is the sampling time; S2: Based on the adjustable-scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, the resonance frequency band caused by the fault is located, thereby adaptively determining the mode number of the characteristic mode decomposition and the fault characteristics corresponding to each mode. frequency; S3: Based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS ₂ ) as the objective function, the collected composite fault signal is adaptively decomposed into several modal components; S4: Calculate the envelope spectrum of the decomposed modes to achieve accurate extraction of rotating machinery fault features and fault diagnosis; Wherein, in step S2, the adjustable-scale sliding window filtering method based on envelope autocorrelation kurtosis (EAK) as an indicator includes the following steps: S201: Load the original signal x(t) and start iteration i=1; S202: Construct a scale-adaptive bandpass filter based on the second-order Butterworth filter. Assume that the filter bandwidth is α, name it as the scale factor of the filter; set the initial α=0.01f _s , f _s is For the sampling frequency of the original signal, set the window sliding step size to 0.2α; construct a bandpass filter group, perform sliding window filtering on the original signal x(t), and obtain a set of filtered sub-signals {x ₁ , x ₂ ,... , x _END }; S203: Screen out the filtered sub-signal x _max with the largest envelope autocorrelation kurtosis (EAK) value, that is, the sub-signal with the most fault information content. The calculation formula of EAK is as follows: X(t)=hilbert(x(t)) R _XX (τ)＝∫X(t)X(t+τ)dt In the _formula , X(t) represents the hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; represents the time delay; EAK(x(t)) represents the EAK value of the signal x(t), N represents the number of sampling points of the original signal x(t), min(·) represents the minimum value function in the MATLAB toolkit, ∫· represents The integration symbol, ∑· represents the summation symbol, R _XX (ii) represents the ii-th value of R _XX (τ), ii=1, 2,..., N/2; S204: Adaptively adjust the scale factor α value of the sliding window corresponding to x _max so that it covers the entire resonance frequency band excited by the fault as much as possible, and then update x _max to obtain the i-th optimal filtered signal x _{best_i} ; α The adaptive adjustment method is as follows: Update α = α + 10% α, that is, gradually increase the filter bandwidth, and the center frequency of the sliding window corresponding to x _max remains unchanged. Use the filter after the scale update to filter x(t) to obtain a new filtered signal x′ _max ; If EAK(x′ _max )>EAK(x _max ), update x _max =x′ _max ; otherwise, stop updating α and obtain the filtered sub-signal x _{best_i} with the optimal scale; S205: Calculate the envelope harmonic product spectrum (EHPS) of the optimal filtered sub-signal x _{best_i} , and select the frequency corresponding to the position of the global maximum value v _i of the EHPS as the estimated i-th fault characteristic frequency f _i ; the calculation formula of EHPS is as follows : In the formula, P(w) represents EHPS, F(w) represents the Fourier transform amplitude of the signal envelope, m=1, 2,..., M, M represents the number of harmonics considered, Π· represents continuous multiplication operator symbol; S206: If i>1, calculate the EHPS reduction rate r _i of x _{best_i} , otherwise return to step S202 and update x(t)=x(t)-x _{best_i} , i=i+1; the calculation formula of r _i is as follows: S207: When i>3, if r _i-2 > r _i-3 and r _i-2 > r _i-1 , let index=i-2, freq=[f ₁ , f ₂ ,..., f _index ], obtain the decomposed target mode number K = index, and the estimated fault characteristic frequency required for each mode decomposition is freq = [f ₁ , f ₂ , ..., f _index ]; otherwise return to step S202, Update x(t)=x(t)-x _{best_i} , i=i+1; when i≤3, also return to step S202, update x(t)=x(t)-x _{best_i} , i=i+1 ; Among them, in step S3, the eigenmode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function includes the following steps: S301: Load the original signal x=x(t), and use the target mode number K and the fault characteristic frequency freq obtained in step S207 as the input parameters of the eigenmode decomposition; S302: Set the filter length L, use Q Hanning windows to initialize the FIR filter bank, and obtain the initial FIR filter bank {h ₁ , h ₂ , ..., h _Q }; the upper and lower cutoff frequencies of the FIR filter h _l and h _u can be expressed as: In the formula, f _s is the sampling frequency of the original signal, q is the sequence number of the filter, q=1, 2,..., Q; S303: Use the deconvolution method to filter the original signal x to obtain the decomposed mode, and maximize the constrained decomposition by maximizing the second-order cyclostationarity (ICS ₂ ) of the decomposed mode. The expression of ICS ₂ is as follows: In the formula, u _q represents the decomposition mode, r is the sample index, β represents the number of sample points corresponding to the fault characteristic frequency, L is the filter length, N represents the number of sampling points of the original signal, u _q [n] represents the third of the matrix u _q n elements, n=1, 2,..., N-L+1, j represents the complex unit, |·| represents the absolute value operation symbol, ||·|| represents the norm operation symbol; The expression of ICS ₂ is rewritten into the following matrix form: E＝[e ₁ , ..., e _r , ..., e _R ] e _r =[e ^{-j2πrβ(L-1)} ,...,e ^{-j2πrβ(N-1)} ] ^H In the formula, D represents the diagonal matrix of u _q , u _q [n] represents the n-th element of the matrix u _q , n=1, 2,..., N-L+1, diag(·) represents the MATIAB tool The diagonal matrix function in the package; E represents a matrix composed of R matrices e _r sequentially spliced horizontally, r is the sample index, r=1, 2,..., R, R is a positive integer, the recommended value range is 50 to 200, the superscript H of the matrix represents the conjugate transpose of the matrix, and j represents the complex unit; Therefore, the decomposition process of the eigenmode decomposition method based on ICS ₂ is equivalent to finding the solution to the following constrained problem: In the formula, h _q [l] is the q-th filter, its length is L, l=1, 2,..., L; Represents the filter h _q [l] when solving to maximize ICS ₂ ; The deconvolution method is used to iteratively filter to obtain the decomposition mode. The iterative process of the deconvolution method is as follows: u _q =Xh _q In the formula, u _q represents the decomposition mode, h _q represents the q-th filter, q=1, 2,..., Q, X represents the shape of the original signal after equal length interception is (N-L+1) ×L matrix, L is the filter length, and N represents the number of sampling points of the original signal; Then, the decomposed modal ICS ₂ can be defined as: In the formula, R _XwX and R _XX represent the weighted correlation matrix and the correlation matrix respectively, and W is the intermediary matrix that controls the weighted correlation matrix; Calculate the ICS ₂ difference of the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than eta, the algorithm is deemed to have converged and the iteration will be stopped, that is: Among them, eta is the convergence condition, the recommended value of eta is 10 ^-6 , and y is the current number of iterations; In order to obtain the optimal filter when maximizing ICS ₂ , it is equivalent to solving the eigenvector related to the maximum eigenvalue λ of the following eigenvalue problem: R _XwX h _q =R _XX h _q λ In the formula, λ can be obtained by solving the above formula, λ represents the maximum eigenvalue, that is, the optimal ICS ₂ ; Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposition mode u _q ; S304: Mode selection; calculate the ICS ₂ of the decomposed mode u _q , select the largest mode of ICS ₂ as the k-th final decomposed mode _uk , and then update the original signal x=xu _k , when the number of decomposed modes reaches step The decomposition ends when the target mode number K obtained in S207. 2. A rotating machinery composite fault diagnosis system based on adaptive eigenmode decomposition, which is characterized by: The signal acquisition module is used to collect the rotating machinery composite fault signal x(t) to be diagnosed, where t is the sampling time; The parameter determination module is used to locate the resonance frequency band caused by the fault based on the adjustable-scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, thereby adaptively determining the mode number and each mode of the characteristic mode decomposition. Corresponding fault characteristic frequency; The signal decomposition module is used to adaptively decompose the collected composite fault signal into several modal components based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS ₂ ) as the objective function; The fault diagnosis module is used to calculate the envelope spectrum of the decomposed modes to achieve accurate extraction of rotating machinery fault features and fault diagnosis; Wherein, the parameter determination module includes: The sliding window filtering unit is used to construct a bandpass filter bank and perform sliding window filtering on the original signal x(t) to obtain a set of filtered sub-signals {x ₁ , x ₂ ,..., x _END }; The target screening unit is used to screen out the filtered sub-signal x _max with the largest envelope autocorrelation kurtosis (EAK) value, that is, the sub-signal with the most fault information content. The calculation formula of EAK is as follows: X(t)=hilbert(x(t)) R _XX (τ)＝∫X(t)X(t+τ)dt In the _formula , X(t) represents the hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; represents the time delay; EAK(x(t)) represents the EAK value of the signal x(t), N represents the number of sampling points of the original signal x(t), min(·) represents the minimum value function in the MATLAB toolkit, ∫· represents The integration symbol, ∑· represents the summation symbol, R _XX (ii) represents the ii-th value of R _XX (τ), ii=1, 2,..., N/2; The filter scale _adjustment unit is used to adaptively adjust the scale factor _α value of the sliding window corresponding to Signal x _{best_i} ; α adaptive adjustment method is as follows: Update α = α + 10% α, that is, gradually increase the filter bandwidth, and the center frequency of the sliding window corresponding to x _max remains unchanged. Use the filter after the scale update to filter x(t) to obtain a new filtered signal x′ _max ; If EAK(x′ _max )>EAK(x _max ), update x _max =x′ _max ; otherwise, stop updating α and obtain the filtered sub-signal x _{best_i} with the optimal scale; The frequency estimation unit is used to calculate the EHPS of the optimal filtered sub-signal x _{best_i} , and select the frequency corresponding to the position of the EHPS global maximum value v _i as the estimated i-th fault characteristic frequency _fi ; the calculation formula of EHPS is as follows: In the formula, P(w) represents EHPS, F(w) represents the Fourier transform amplitude of the signal envelope, m=1, 2,..., M, M represents the number of harmonics considered, Π· represents continuous multiplication operator symbol; The target mode number determination unit is used to locate the number of resonance frequency bands in the signal to determine the number of modes that need to be decomposed; if i>1, calculate the EHPS reduction rate r _i of x _{best_i} , otherwise return to the sliding window filtering unit, Update x(t)=x(t)-x _{best_i} , i=i+1, x represents the original signal; the calculation formula of r _i is as follows: When i>3, if r _i-2 > r _i-3 and r _i-2 > r _i-1 , let index=i-2, freq=[f ₁ , f ₂ ,..., f _index ] , obtain the decomposed target modal number K=index, and the estimated fault characteristic frequency required for each modal decomposition is freq=[f ₁ , f ₂ ,..., fi _index ]; otherwise return to the sliding window filter unit, Update x(t)=x(t)-x _{best_i} , i=i+1; when i≤3, also return to the sliding window filter unit, update x(t)=x(t)-x _{best_i} , i=i +1; Wherein, the signal decomposition module includes: Filter initialization unit, used for FIR filter initialization; set the filter length L, use Q Hanning windows to initialize the FIR filter bank, and obtain the initial FIR filter bank {h ₁ , h ₂ ,..., h _Q } ;The upper and lower cutoff frequencies h _l and _hu of the FIR filter can be expressed as: In the formula, f _s is the sampling frequency of the original signal, q is the sequence number of the filter, q=1, 2,..., Q; Objective function unit, used to establish the objective function ICS ₂ of the eigenmode decomposition method: In the formula, u _q represents the decomposition mode, r is the sample index, β represents the number of sample points corresponding to the fault characteristic frequency, L is the filter length, N represents the number of sampling points of the original signal, u _q [n] represents the third of the matrix u _q n elements, n=1, 2,..., N-L+1, j represents the complex unit, |·| represents the absolute value operation symbol, ||·|| represents the norm operation symbol; The expression of ICS ₂ can be rewritten into the following matrix form: E＝[e ₁ ,..., e _r ,..., e _R ] e _r =[e ^{-j2πrβ(L-1)} ,...,e ^{-j2πrβ(N-1)} ] ^H In the formula, D represents the diagonal matrix of u _q , u _q [n] represents the n-th element of the matrix u _q , n=1, 2,..., N-L+1, diag(·) represents the MATIAB tool The diagonal matrix function in the package; E represents a matrix composed of R matrices e _r sequentially spliced horizontally, r is the sample index, r=1, 2,..., R, R is a positive integer, the recommended value range is 50 to 200, the superscript H of the matrix represents the conjugate transpose of the matrix, and j represents the complex unit; Therefore, the decomposition process of the eigenmode decomposition method based on ICS ₂ is equivalent to finding the solution to the following constrained problem: In the formula, h _q [l] is the q-th filter, its length is L, l=1, 2,..., L; Represents the filter h _q [l] when solving to maximize ICS ₂ ; The deconvolution iteration unit is used to iteratively filter using the deconvolution method to obtain the decomposition mode. The iteration process of the deconvolution method is as follows: u _q =Xh _q In the formula, u _q represents the decomposition mode, h _q represents the q-th filter, q=1, 2,..., Q, X represents the shape of the original signal after equal length interception is (N-L+1) ×L matrix, L is the filter length, and N represents the number of sampling points of the original signal; Then, the decomposed modal ICS ₂ can be defined as: In the formula, R _XWX and R _XX represent the weighted correlation matrix and the correlation matrix respectively, and W is the intermediary matrix that controls the weighted correlation matrix; Calculate the ICS ₂ difference of the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than eta, the algorithm is deemed to have converged and the iteration will be stopped, that is: Among them, eta is the convergence condition, the recommended value of eta is 10 ^-6 , and y is the current number of iterations; In order to obtain the optimal filter when maximizing ICS ₂ , it is equivalent to solving the eigenvector related to the maximum eigenvalue λ of the following eigenvalue problem: R _XwX h _q =R _XX h _q λ In the formula, λ can be obtained by solving the above formula, λ represents the maximum eigenvalue, that is, the optimal ICS ₂ ; Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposition mode u _q ; The mode selection unit is used to select the optimal mode; calculate the ICS ₂ of the decomposed mode u _q , select the mode with the largest ICS ₂ as the k-th final decomposed mode _uk , and then update the original signal x=xu _k , The decomposition ends when the number of decomposition modes reaches the target number of modes K.