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

Abstract

A rotating machinery composite fault diagnosis method and system based on adaptive characteristic mode decomposition belongs to the field of mechanical signal processing and equipment fault diagnosis. The present invention collects the composite fault signal of the rotating machinery to be diagnosed; locates the resonant 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 number of modes of characteristic mode decomposition and the fault characteristic frequency corresponding to each mode; adaptively decomposes the collected composite fault signal into several modal components based on the characteristic mode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function; and calculates the envelope spectrum of the decomposed mode. The present invention overcomes the problem of heavy reliance on prior parameters such as the number of decomposed modes and the fault characteristic frequency, and can accurately locate the resonant frequency band caused by the fault from the composite fault signal of the rotating machinery with strong periodic noise interference, and accurately separate and extract each single fault component.

Description

Technical Field

[0001] The present 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 Art

[0002] Rotating machinery is widely used in aerospace, transportation, energy and chemical industries and plays an important role. Due to long-term high-speed, heavy-load, and variable operating 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 affect each other, resulting in compound failures. In severe cases, it will not only cause immeasurable economic losses, but even cause casualties, with severe consequences. When a rotating machinery has a compound failure, the vibration signal shows the characteristics of simultaneous coupling of multiple periodic components and multiple pulse components. The periodic motion of the remaining components will introduce periodic noise interference, which greatly affects the accuracy of compound fault diagnosis. Therefore, the development of accurate and effective rotating machinery compound fault diagnosis methods and systems under strong periodic noise interference is of great significance to improving the safety and reliability of equipment.

[0003] In recent years, a large number of scholars have studied the theory of mechanical signal decomposition and successfully applied it to the field of mechanical fault diagnosis. Signal decomposition theory can decompose the original signal of multi-component aliasing into several regular simple modes, providing a theoretical basis for the separation of mechanical complex fault characteristics. However, the classical decomposition methods all have their own limitations, and cannot effectively decompose complex fault signals without the assistance of other methods. For example, empirical mode decomposition (EMD) has the 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, the decomposition effect of empirical wavelet transform (EWT) is subject to frequency band segmentation, and the decomposition effect of variational mode decomposition (VMD) is seriously dependent on the prior knowledge of the initial parameters. In 2022, the newly proposed eigenmode decomposition (FMD) will gradually be promoted and applied due to its unique advantages in the field of rotating machinery fault diagnosis. However, the FMD method is difficult to accurately estimate the true fault frequency in the case of strong periodic interference noise, which will affect the characteristic decomposition effect of the FMD method; at the same time, the FMD method needs to manually set the decomposition mode number based on prior knowledge, and the selection of the decomposition mode number also affects the decomposition effect. Therefore, in view of the limitations of the existing methods, the present invention discloses a rotating machinery composite fault diagnosis method and system based on adaptive characteristic mode decomposition, which realizes the accurate diagnosis of mechanical composite faults under strong noise. Summary of the invention

[0004] In order to solve the defects of the prior art, especially the problems of poor adaptability of the characteristic mode decomposition (FMD) method and insufficient robustness under strong periodic noise interference, the present invention proposes a rotating machinery composite fault diagnosis method and system based on adaptive characteristic mode decomposition, which can adaptively determine the number of decomposed modes, adaptively estimate the fault characteristic frequency, accurately separate and extract each single fault component from the rotating machinery composite fault signal with strong periodic noise interference, and realize accurate diagnosis of rotating machinery composite faults.

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

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

[0007] S2: Based on the adjustable scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as the index, the resonance frequency band caused by the fault is located, so as to adaptively determine the number of modes of characteristic mode decomposition and the fault characteristic frequency corresponding to each mode;

[0008] S3: Adaptively decompose the collected composite fault signal into several modal components based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function;

[0009] S4: Calculate the envelope spectrum of the decomposed mode to achieve accurate extraction of rotating machinery fault features and fault diagnosis;

[0010] 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 comprises the following steps:

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

[0012] S202: Construct a scale-adaptive bandpass filter based on the second-order Butterworth filter. Assume that the filter bandwidth is α, which is named the scale factor of the filter. Set the initial α = 0.01f s , f s is the sampling frequency of the original signal, and the window sliding step is set 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 {x1, x2, ..., x END};

[0013] S203: Filter out the filter sub-signal x with the largest envelope autocorrelation kurtosis (EAK) value max , which is the sub-signal with the most fault information content. The calculation formula of EAK is as follows:

[0014] X(t) = hilbert(x(t))

[0015] R XX (τ)=∫X(t)X(t+τ)dt

[0016]

[0017] Where X(t) represents the Hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; R XX (τ) represents the autocorrelation function of X(t), τ 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 function in the MATLAB toolkit, ∫· represents the integral symbol, ∑· represents the summation symbol, R XX (ii) indicates R XX the ii-th value of (τ), ii=1, 2, ..., N / 2;

[0018] S204: Adaptive Adjustment x max The scale factor α value of the corresponding sliding window is set to cover the entire resonance frequency band excited by the fault as much as possible, and then update x max So as to obtain the optimal filtered signal x of the i-th scale best_i The α adaptive adjustment method is as follows:

[0019] Update α = α + 10% α, that is, gradually increase the filter bandwidth, x max The center frequency of the corresponding sliding window remains unchanged, and the filter with the updated scale is used 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 optimal scale filter sub-signal x best_i ;

[0020] S205: Calculate the optimal filter sub-signal x best_i The envelope harmonic product spectrum (EHPS) of the i The frequency corresponding to the location is taken as the estimated i-th fault characteristic frequency f i The calculation formula of EHPS is as follows:

[0021]

[0022] Where 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, and Π· represents the symbol of the continuous multiplication operation;

[0023] S206: If i>1, calculate x best_i The EHPS decline rate r i Otherwise, return to step S202 and update x(t)=x(t)-x best_i , i=i+1, x represents the original signal. i The calculation formula is as follows:

[0024]

[0025] S207: When i>3, if r i-2 >r i-3 And r i-2 >r i-1 , let index=i-2, freq=[f1, f2,..., f index ], the target modal number K = index is obtained, and the fault characteristic frequency estimation value required for each modal decomposition is freq = [f1, f2, ..., f index ]; otherwise, return to step S202 and 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;

[0026] In one embodiment of the present invention, in step S3, the eigenmode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function comprises the following steps:

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

[0028] S302: Set the filter length L, use Q Hanning windows to initialize the FIR filter group, and obtain the initial FIR filter group {h1, h2, ..., h Q}. The upper and lower cutoff frequencies of the FIR filter are h l 、h u It can be expressed as:

[0029]

[0030] In the formula, f s is the sampling frequency of the original signal, q is the filter number, q = 1, 2, ..., Q;

[0031] S303: The original signal x is filtered using a deconvolution method to obtain a decomposed mode, and the second-order cyclostationarity (ICS2) of the decomposed mode is maximized to constrain the decomposition. The expression of ICS2 is as follows:

[0032]

[0033]

[0034]

[0035] 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, and u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L+1, j represents a complex number unit, |·| represents an absolute value operator, and ||·|| represents a norm operator;

[0036] The expression of ICS2 can be rewritten into the following matrix form:

[0037]

[0038]

[0039] E=[e1,...,e r , ..., e R ]

[0040] e r =[e -j2πrβ(L-1) , ..., e -j2πrβ(N-1) ] H

[0041] Where D represents u q The diagonal matrix, u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L+1, diag(·) represents the diagonal matrix function in the MATLAB toolkit; E represents the diagonal matrix consisting of R matrices e r The matrix is ​​formed by sequentially splicing horizontally, r is the sample index, r = 1, 2, ..., R, R is a positive integer, and 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;

[0042] Therefore, the decomposition process of the eigenmode decomposition method based on ICS2 is equivalent to solving the following constrained problem:

[0043]

[0044]

[0045] In the formula, h q [l] is the qth filter, whose length is L, L = 1, 2, ..., L; It means to find the filter h that maximizes ICS2. q [l];

[0046] The deconvolution method is used to iteratively filter and obtain the decomposed mode. The iterative process of the deconvolution method is as follows:

[0047] h q =Xh q

[0048]

[0049] In the formula, u q represents the decomposition mode, h q represents the qth filter, q = 1, 2, ..., Q, X represents the matrix of the original signal after being truncated with equal length, with the shape of (N-L+1) × L, L is the filter length, and N represents the number of sampling points of the original signal;

[0050] Then, the decomposed modal ICS2 can be defined as:

[0051]

[0052]

[0053] In the formula, R XWX and R XX denote the weighted correlation matrix and the correlation matrix respectively, W is the mediation matrix controlling the weighted correlation matrix;

[0054] Calculate the ICS2 difference between the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than η, the algorithm is considered to have converged and the iteration is stopped, that is:

[0055]

[0056] Among them, η is the convergence condition, and the recommended value of η is 10 -6 , y is the current iteration number;

[0057] In order to obtain the optimal filter when ICS2 is maximized, it is equivalent to solving the eigenvalue problem of the following eigenvalue problem related to the maximum eigenvalue λ:

[0058] R XWX h q =R XX h q λ

[0059] In the formula, solving the above formula can get λ, λ represents the maximum eigenvalue, that is, the optimal ICS2;

[0060] Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposed mode u q ;

[0061] S304: Mode selection. Calculate the decomposed mode u q ICS2, select the largest mode of ICS2 as the kth final decomposition mode u k , then update the original signal x = xu k , the decomposition is terminated when the number of decomposed modes reaches the target number of modes K obtained in step S207.

[0062] Another object of the present invention is to provide a rotating machinery composite fault diagnosis system based on adaptive eigenmode decomposition, characterized in that it includes:

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

[0064] A parameter determination module is used to locate the resonance frequency band caused by the fault based on an adjustable scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, so as to adaptively determine the number of modes of characteristic mode decomposition and the fault characteristic frequency corresponding to each mode;

[0065] A signal decomposition module is used to adaptively decompose the collected composite fault signal into several modal components based on the characteristic mode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function;

[0066] Fault diagnosis module, used to calculate the envelope spectrum of decomposed modes, to achieve accurate extraction of rotating machinery fault features and fault diagnosis;

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

[0068] The sliding window filter unit is used to construct a bandpass filter bank to perform sliding window filtering on the original signal x(t) to obtain a set of filtered sub-signals {x1, x2, ..., x END};

[0069] The target screening unit is used to screen out the filter sub-signal x with the largest envelope autocorrelation kurtosis (EAK) value max , which is the sub-signal with the most fault information content. The calculation formula of EAK is as follows:

[0070] X(t) = hilbert(x(t))

[0071] RXX (τ)=∫X(t)X(t+τ)dt

[0072]

[0073] Where X(t) represents the Hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; R XX (τ) represents the autocorrelation function of X(t), τ 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 function in the MATLAB toolkit, ∫· represents the integral symbol, ∑· represents the summation symbol, R XX (ii) indicates R XX the ii-th value of (τ), ii=1, 2, ..., N / 2;

[0074] Filter scale adjustment unit, used to adaptively adjust x max The scale factor α value of the corresponding sliding window is set to cover the entire resonance frequency band excited by the fault as much as possible, and then update x max So as to obtain the optimal filtered signal x of the i-th scale best_i The α adaptive adjustment method is as follows:

[0075] Update α = α + 10% α, that is, gradually increase the filter bandwidth, x max The center frequency of the corresponding sliding window remains unchanged, and the filter with the updated scale is used 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 optimal scale filter sub-signal x best_i ;

[0076] Frequency estimation unit, used to calculate the optimal filter sub-signal x best_i EHPS, select the EHPS global maximum value v i The frequency corresponding to the location is taken as the estimated i-th fault characteristic frequency f i The calculation formula of EHPS is as follows:

[0077]

[0078] Where 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, and Π· represents the symbol of the continuous multiplication operation;

[0079] The target modal number determination unit is used to locate the number of resonant frequency bands in the signal, thereby determining the number of modes to be decomposed. If i>1, calculate x best_i The EHPS decline rate r 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. i The calculation formula is as follows:

[0080]

[0081] When i>3, if r i-2 >r i-3 And r i-2 >r i-1 , let index=i-2, freq=[f1, f2,..., f index ], the target modal number K = index is obtained, and the fault characteristic frequency estimation value required for each modal decomposition is freq = [f1, f2, ..., f index ]; otherwise, return to the sliding window filter unit and update x(t) = x(t) - x best_i , i = i + 1. When i ≤ 3, return to the sliding window filter unit and update x(t) = x(t) - x best_i , i=i+1;

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

[0083] The filter initialization unit is used to initialize the FIR filter. Set the filter length L, use Q Hanning windows to initialize the FIR filter group, and obtain the initial FIR filter group {h1, h2, ..., h Q}. The upper and lower cutoff frequencies of the FIR filter are h l 、h u It can be expressed as:

[0084]

[0085] In the formula, f s is the sampling frequency of the original signal, q is the filter number, q = 1, 2, ..., Q;

[0086] The objective function unit is used to establish the objective function ICS2 of the characteristic mode decomposition method:

[0087]

[0088]

[0089]

[0090] 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, and u q [n] represents the matrix u q The nth element of , n=1,2,…,N-L+1, j represents the complex unit, |·| represents the absolute value operator, ||·|| represents the norm operator;

[0091] The expression of ICS2 can be rewritten into the following matrix form:

[0092]

[0093]

[0094] E=[e1,…,e r ,…,e R ]

[0095] e r =[e -j2πrβ(L-1) ,…,e -j2πrβ(N-1 )] H

[0096] Where D represents u q The diagonal matrix, u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L + 1, diag (·) represents the diagonal matrix function in the MATLAB toolkit; E represents the diagonal matrix consisting of R matrices e r The matrix is ​​formed by horizontal splicing in sequence, r is the sample index, r = 1, 2, ..., R, R is a positive integer, and 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;

[0097] Therefore, the decomposition process of the eigenmode decomposition method based on ICS2 is equivalent to solving the following constrained problem:

[0098]

[0099]

[0100] In the formula, h q [l] is the qth filter, whose length is L, l = 1, 2, ..., L; It means to find the filter h that maximizes ICS2. q [l];

[0101] The deconvolution iteration unit is used to obtain the decomposed mode by iterative filtering using the deconvolution method. The iterative process of the deconvolution method is as follows:

[0102] u q =Xh q

[0103]

[0104] In the formula, u q represents the decomposition mode, h q represents the qth filter, q = 1, 2, ..., Q, X represents the matrix of the original signal after being cut into equal lengths, the shape of which is (N-L+1) × L, L is the filter length, and N represents the number of sampling points of the original signal;

[0105] Then, the decomposed modal ICS2 can be defined as:

[0106]

[0107]

[0108] In the formula, R XWX and R XX denote the weighted correlation matrix and the correlation matrix respectively, W is the mediation matrix controlling the weighted correlation matrix;

[0109] Calculate the ICS2 difference between the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than η, the algorithm is considered to have converged and the iteration is stopped, that is:

[0110]

[0111] Among them, η is the convergence condition, and the recommended value of η is 10 -6 , y is the current iteration number;

[0112] In order to obtain the optimal filter when ICS2 is maximized, it is equivalent to solving the eigenvalue problem of the following eigenvalue problem related to the maximum eigenvalue λ:

[0113] R XWX h q =R XX h q λ

[0114] In the formula, solving the above formula can get λ, λ represents the maximum eigenvalue, that is, the optimal ICS2;

[0115] Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposed mode u q ;

[0116] Mode selection unit, used to select the optimal mode. Calculate the decomposed mode uq ICS2, select the largest mode of ICS2 as the kth final decomposition mode u k , then update the original signal x = xu k , when the number of decomposed modes reaches the target mode number K, the decomposition is terminated.

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

[0118] The beneficial technical effects of the present invention are:

[0119] 1. The adaptive optimization of decomposition modal number and adaptive estimation strategy of fault frequency proposed in the present invention can accurately locate the resonant frequency band caused by rotating machinery fault under strong noise through sliding window filtering guided by EAK index, so as to further obtain the decomposition modal number required for decomposition and its corresponding fault frequency. 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 cyclostationary blind deconvolution (CYCBD), VMD, etc., to improve the effects of these methods.

[0120] 2. The eigenmode decomposition method based on ICS2 proposed in this invention can accurately extract weak fault information of rotating machinery under strong noise background. Compared with the correlation kurtosis (CK) index, using ICS2 as the deconvolution objective function can better extract the characteristic components of rotating machinery faults under strong noise. ICS2 has better sensitivity to rotating machinery fault information and noise robustness.

[0121] 3. The present invention has good adaptability and accuracy for the extraction of single and compound fault features under strong periodic noise interference conditions, and has an outstanding feature decomposition effect, especially in the field of compound fault diagnosis. BRIEF DESCRIPTION OF THE DRAWINGS

[0122] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0123] Figure 1 It is a schematic diagram of a rotating machinery composite fault diagnosis method and system based on adaptive characteristic mode decomposition provided by an embodiment of the present invention;

[0124] Figure 2 It is a detailed flow chart of a rotating machinery composite fault diagnosis method based on adaptive characteristic mode decomposition provided by an embodiment of the present invention;

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

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

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

[0128] Figure 6 It is an envelope spectrum of a composite fault signal adaptively decomposed mode of a rotating machinery composite fault diagnosis method based on adaptive characteristic mode decomposition provided by an embodiment of the present invention;

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

[0130] Figure 8 3 is a diagram showing the result of signal processing using the FMD method in an embodiment of the present invention. DETAILED DESCRIPTION

[0131] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0132] Figure 1 is a rotating machinery composite fault diagnosis method and system schematic diagram based on adaptive characteristic mode decomposition provided by an embodiment of the present invention, and its detailed process is as follows Figure 2 See Figure 1 , A rotating machinery composite fault diagnosis method based on adaptive eigenmode decomposition The method comprises the following steps:

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

[0134] S2: Based on the adjustable scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as the index, the resonance frequency band caused by the fault is located, so as to adaptively determine the number of modes of characteristic mode decomposition and the fault characteristic frequency corresponding to each mode;

[0135] S3: Adaptively decompose the collected composite fault signal into several modal components based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function;

[0136] S4: Calculate the envelope spectrum of the decomposed mode to achieve accurate extraction of rotating machinery fault features and fault diagnosis.

[0137] 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 comprises the following steps:

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

[0139] S202: Construct a scale-adaptive bandpass filter based on the second-order Butterworth filter. Assume that the filter bandwidth is α, which is named the scale factor of the filter. Set the initial α = 0.01f s , f s is the sampling frequency of the original signal, and the window sliding step is set 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 {x1, x2, ..., x END};

[0140] S203: Filter out the filter sub-signal x with the largest envelope autocorrelation kurtosis (EAK) value max , which is the sub-signal with the most fault information content. The calculation formula of EAK is as follows:

[0141] X(t) = hilbert(x(t))

[0142] R XX (τ)=∫X(t)X(t+τ)dt

[0143]

[0144] Where X(t) represents the Hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; R XX(τ) represents the autocorrelation function of X(t), τ 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 function in the MATLAB toolkit, ∫· represents the integral symbol, ∑· represents the summation symbol, R XX (ii) indicates R XX the ii-th value of (τ), ii=1, 2, ..., N / 2;

[0145] S204: Adaptive Adjustment x max The scale factor α value of the corresponding sliding window is set to cover the entire resonance frequency band excited by the fault as much as possible, and then update x max So as to obtain the optimal filtered signal x of the i-th scale best_i The α adaptive adjustment method is as follows:

[0146] Update α = α + 10% α, that is, gradually increase the filter bandwidth, x max The center frequency of the corresponding sliding window remains unchanged, and the filter with the updated scale is used 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 optimal scale filter sub-signal x best_i ;

[0147] S205: Calculate the optimal filter sub-signal x best_i EHPS, select the EHPS global maximum value v i The frequency corresponding to the location is taken as the estimated i-th fault characteristic frequency f i The calculation formula of EHPS is as follows:

[0148]

[0149] Where 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, and Π· represents the symbol of the continuous multiplication operation;

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

[0151] S206: If i>1, calculate x best_i The EHPS decline rate r i Otherwise, return to step S202 and update x(t)=x(t)-x best_i, i=i+1, x represents the original signal. i The calculation formula is as follows:

[0152]

[0153] S207: When i>3, if r i-2 >r i-3 And r i-2 >r i-1 , let index=i-2, freq=[f1, f2,..., f index ], the target modal number K = index is obtained, and the fault characteristic frequency estimation value required for each modal decomposition is freq = [f1, f2, ..., f index ]; otherwise, return to step S202 and 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:

[0154] In one embodiment of the present invention, in step S3, the eigenmode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function comprises the following steps:

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

[0156] S302: Set the filter length L, use Q Hanning windows to initialize the FIR filter group, and obtain the initial FIR filter group {h1, h2, ..., h Q}. The upper and lower cutoff frequencies of the FIR filter are h l 、h u It can be expressed as:

[0157]

[0158] In the formula, f s is the sampling frequency of the original signal, q is the filter number, q = 1, 2, ..., Q;

[0159] 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;

[0160] S303: The original signal x is filtered using a deconvolution method to obtain a decomposed mode, and the second-order cyclostationarity (ICS2) of the decomposed mode is maximized to constrain the decomposition. The expression of ICS2 is as follows:

[0161]

[0162]

[0163]

[0164] 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, and u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L+1, j represents a complex number unit, |·| represents an absolute value operator, and ||·|| represents a norm operator;

[0165] The expression of ICS2 can be rewritten into the following matrix form:

[0166]

[0167]

[0168] E=[e1,…,e r ,…,e R ]

[0169] e r =[e -j2πrβ(L-1) ,…,e -j2πrβ(N-1 )] H

[0170] Where D represents u q The diagonal matrix, u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L + 1, diag (·) represents the diagonal matrix function in the MATLAB toolkit; E represents the diagonal matrix consisting of R matrices e r The matrix is ​​formed by horizontal splicing in sequence, r is the sample index, r = 1, 2, ..., R, R is a positive integer, and 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;

[0171] Therefore, the decomposition process of the eigenmode decomposition method based on ICS2 is equivalent to solving the following constrained problem:

[0172]

[0173]

[0174] In the formula, h q [l] is the qth filter, whose length is L, L = 1, 2, ..., L; It means to find the filter h that maximizes ICS2. q [l];

[0175] Calculate the ICS2 difference between the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than η, the algorithm is considered to have converged and the iteration is stopped, that is:

[0176]

[0177] Among them, η is the convergence condition, and the recommended value of η is 10 -6 , y is the current iteration number;

[0178] The deconvolution method is used to iteratively filter and obtain the decomposed mode. The iterative process of the deconvolution method is as follows:

[0179] u q =Xh q

[0180]

[0181] In the formula, u q represents the decomposition mode, h q represents the qth filter, q = 1, 2, ..., Q, X represents the matrix of the original signal after being truncated with equal length, with the shape of (N-L+1) × L, L is the filter length, and N represents the number of sampling points of the original signal;

[0182] Then, the decomposed modal ICS2 can be defined as:

[0183]

[0184]

[0185] In the formula, R XWX and R XX Denote the weighted correlation matrix and the correlation matrix respectively, W is the intermediate variable matrix controlling the weighted correlation matrix;

[0186] In order to obtain the optimal filter when ICS2 is maximized, it is equivalent to solving the eigenvalue problem of the following eigenvalue problem related to the maximum eigenvalue λ:

[0187] R XWX h q =R XX h q λ

[0188] In the formula, solving the above formula can get λ, λ represents the maximum eigenvalue, that is, the optimal ICS2;

[0189] Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposed mode u q ;

[0190] S304: Mode selection. Calculate the decomposed mode u q ICS2, select the largest mode of ICS2 as the kth final decomposition mode u k , then update the original signal x = xu k , the decomposition is terminated when the number of decomposed modes reaches the target number of modes K obtained in step S207.

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

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

[0193] A parameter determination module is used to locate the resonance frequency band caused by the fault based on an adjustable scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, so as to adaptively determine the number of modes of characteristic mode decomposition and the fault characteristic frequency corresponding to each mode;

[0194] A signal decomposition module is used to adaptively decompose the collected composite fault signal into several modal components based on the characteristic mode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function;

[0195] Fault diagnosis module, used to calculate the envelope spectrum of decomposed modes, to achieve accurate extraction of rotating machinery fault features and fault diagnosis;

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

[0197] The sliding window filter unit is used to construct a bandpass filter bank to perform sliding window filtering on the original signal x(t) to obtain a set of filtered sub-signals {x1, x2, ..., x END};

[0198] The target screening unit is used to screen out the filter sub-signal x with the largest envelope autocorrelation kurtosis (EAK) value max , which is the sub-signal with the most fault information content. The calculation formula of EAK is as follows:

[0199] X(t) = hilbert(x(t))

[0200] R XX (τ)=∫X(t)X(t+τ)dt

[0201]

[0202] Where X(t) represents the Hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; R XX (τ) represents the autocorrelation function of X(t), τ 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 function in the MATLAB toolkit, ∫· represents the integral symbol, ∑· represents the summation symbol, R XX (ii) indicates R XX the ii-th value of (τ), ii=1, 2, ..., N / 2;

[0203] Filter scale adjustment unit, used to adaptively adjust x max The scale factor α value of the corresponding sliding window is set to cover the entire resonance frequency band excited by the fault as much as possible, and then update x max So as to obtain the optimal filtered signal x of the i-th scale best_i The α adaptive adjustment method is as follows:

[0204] Update α = α + 10% α, that is, gradually increase the filter bandwidth, x max The center frequency of the corresponding sliding window remains unchanged, and the filter with the updated scale is used 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 optimal scale filter sub-signal x best_i ;

[0205] Frequency estimation unit, used to calculate the optimal filter sub-signal x best_i EHPS, select the EHPS global maximum value v i The frequency corresponding to the location is taken as the estimated i-th fault characteristic frequency f i The calculation formula of EHPS is as follows:

[0206]

[0207] Where 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, and Π· represents the symbol of the continuous multiplication operation;

[0208] The target modal number determination unit is used to locate the number of resonant frequency bands in the signal, thereby determining the number of modes to be decomposed. If i>1, calculate x best_i The EHPS decline rate r iOtherwise, return to the sliding window filter unit and update x(t) = x(t) - x best_i , i=i+1, x represents the original signal. i The calculation formula is as follows:

[0209]

[0210] When i>3, if r i-2 >r i-3 And r i-2 >r i-1 , let index=i-2, freq=[f1, f2,..., f index ], the target modal number K = index is obtained, and the fault characteristic frequency estimation value required for each modal decomposition is freq = [f1, f2, ..., f index ]; otherwise, return to the sliding window filter unit and update x(t) = x(t) - x best_i , i = i + 1. When i ≤ 3, return to the sliding window filter unit and update x(t) = x(t) - x best_i , i=i+1;

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

[0212] The filter initialization unit is used to initialize the FIR filter. Set the filter length L, use Q Hanning windows to initialize the FIR filter group, and obtain the initial FIR filter group {h1, h2, ..., h Q}. The upper and lower cutoff frequencies of the FIR filter are h l 、h u It can be expressed as:

[0213]

[0214] In the formula, f s is the sampling frequency of the original signal, q is the filter number, q = 1, 2, ..., Q;

[0215] The objective function unit is used to establish the objective function ICS2 of the characteristic mode decomposition method:

[0216]

[0217]

[0218]

[0219] In the formula, u qrepresents 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, and u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L+1, j represents a complex number unit, |·| represents an absolute value operator, and ||·|| represents a norm operator;

[0220] The expression of ICS2 can be rewritten into the following matrix form:

[0221]

[0222]

[0223] E=[e1,,..,e r , ..., e R ]

[0224] e r =[e -j2πrβ(L-1) , ..., e -j2πrβ(N-1) ] H

[0225] Where D represents u q The diagonal matrix, u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L+1, diag(·) represents the diagonal matrix function in the MATLAB toolkit; E represents the diagonal matrix consisting of R matrices e r The matrix is ​​formed by horizontal splicing in sequence, r is the sample index, r = 1, 2, ..., R, R is a positive integer, and 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;

[0226] Therefore, the decomposition process of the eigenmode decomposition method based on ICS2 is equivalent to solving the following constrained problem:

[0227]

[0228]

[0229] In the formula, h q [l] is the qth filter, whose length is L, l = 1, 2, ..., L; It means to find the filter h that maximizes ICS2. q [l];

[0230] The deconvolution iteration unit is used to obtain the decomposed mode by iterative filtering using the deconvolution method. The iterative process of the deconvolution method is as follows:

[0231] u q =Xh q

[0232]

[0233] In the formula, u q represents the decomposition mode, h q represents the qth filter, q = 1, 2, ..., Q, X represents the matrix of the original signal after being truncated with equal length, with the shape of (N-L+1) × L, L is the filter length, and N represents the number of sampling points of the original signal;

[0234] Then, the decomposed modal ICS2 can be defined as:

[0235]

[0236]

[0237] In the formula, R XWX and R XX denote the weighted correlation matrix and the correlation matrix respectively, W is the mediation matrix controlling the weighted correlation matrix;

[0238] Calculate the ICS2 difference between the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than η, the algorithm is considered to have converged and the iteration is stopped, that is:

[0239]

[0240] Among them, η is the convergence condition, and the recommended value of η is 10 -6 , y is the current iteration number;

[0241] In order to obtain the optimal filter when ICS2 is maximized, it is equivalent to solving the eigenvalue problem of the following eigenvalue problem related to the maximum eigenvalue λ:

[0242] R XWX h q =R XX h q λ

[0243] In the formula, solving the above formula can get λ, λ represents the maximum eigenvalue, that is, the optimal ICS2;

[0244] Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposed mode u q ;

[0245] Mode selection unit, used to select the optimal mode. Calculate the decomposed mode u q ICS2, select the largest mode of ICS2 as the kth final decomposition mode u k, then update the original signal x = xu k , the decomposition ends when the number of decomposed modes reaches the target number of modes.

[0246] The specific embodiments are as follows:

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

[0248] Table 130310 Bearing parameters and theoretical fault characteristic frequencies

[0249]

[0250] The method of the present invention is used to analyze the collected composite fault bearing signal. First, the sampling frequency f is set s =25.6kHz, the number of harmonics considered is M=5, and the resonant frequency band caused by the fault is located based on the adjustable scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, so as to adaptively determine the number of modes of characteristic mode 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 estimation value is freq=[163.87Hz,235.43Hz]. Then, the filter length L=100, η=10 -6 Based on the characteristic mode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function, the collected composite fault signal is adaptively decomposed into several modal components. In this embodiment, the composite fault signal adaptive decomposition mode is obtained as follows: Figure 5 Finally, the envelope spectrum of the decomposed mode is calculated to achieve accurate extraction of rotating machinery fault features and fault diagnosis. In this embodiment, the composite fault signal adaptively decomposes the modal envelope spectrum as shown in Figure 6 As shown. Figure 6It can be seen that the present invention extracts two fault components from the composite fault signal. The envelope spectrum of the first decomposition mode clearly shows the outer ring fault frequency and its 2 to 7 times frequency. The envelope spectrum of the second decomposition mode shows the inner ring fault frequency and its 2 times frequency. Therefore, it can be diagnosed that a composite fault of the inner and outer rings has occurred, which is consistent with the preset fault type.

[0251] To further verify the superiority of the present invention, the same signal was tested using variational mode decomposition (VMD) and eigenmode decomposition (FMD). The decomposed modes and their envelope spectra obtained by VMD are shown in Figure 2. Figure 7 As shown; the decomposed mode and its envelope spectrum obtained by FMD method are as follows Figure 8 As shown. Figure 7 and Figure 8 It can be seen that both VMD and FMD methods cannot effectively realize the composite fault signal feature decomposition and extraction in this embodiment. Therefore, the present invention has advantages in the composite fault signal feature extraction of rotating machinery under strong periodic background noise.

[0252] The above description is only a preferred embodiment of the present application and an explanation of the technical principles used. Those skilled in the art should understand that the scope involved in the present application is not limited to the technical solution formed by a specific combination of the above-mentioned technical features, but should also cover other technical solutions formed by any combination of the above-mentioned technical features or their equivalent features without departing from the inventive concept, such as a technical solution formed by replacing the above-mentioned features with (but not limited to) technical features with similar functions disclosed in this application.

[0253] Except for the technical features described in the specification, 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 here in detail.

Claims

1. A rotating machinery composite fault diagnosis method based on adaptive eigenmode decomposition, characterized in that: The following steps are involved: S1: Collect the composite fault signal x(t) of the rotating machinery to be diagnosed, where t is the sampling time; S2: Based on the adjustable scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as the index, the resonance frequency band caused by the fault is located, so as to adaptively determine the number of modes of characteristic mode decomposition and the fault characteristic frequency corresponding to each mode; S3: Adaptively decompose the collected composite fault signal into several modal components based on the eigenmode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function; S4: Calculate the envelope spectrum of the decomposed mode 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 comprises the following steps: S201: Load the original signal x(t) and start iteration i=1; S202: Construct a scale-adaptive bandpass filter based on a second-order Butterworth filter. Assume that the filter bandwidth is α, which is named as the scale factor of the filter. Set the initial α to 0.01f s , f s is the sampling frequency of the original signal, and the window sliding step is set to 0.2α; a bandpass filter bank is constructed to perform sliding window filtering on the original signal x(t) to obtain a set of filtered sub-signals {x1, x2, ..., x END }; S203: Filter out the filter sub-signal x with the largest envelope autocorrelation kurtosis (EAK) value max , which 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 Where X(t) represents the Hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; R XX (τ) represents the autocorrelation function of X(t), τ 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 function in the MATLAB toolkit, ∫· represents the integral symbol, ∑· represents the summation symbol, R XX (ii) indicates R XX the ii-th value of (τ), ii=1, 2, ..., N / 2; S204: Adaptive Adjustment x max The scale factor α value of the corresponding sliding window is set to cover the entire resonance frequency band excited by the fault as much as possible, and then update x max So as to obtain the optimal filtered signal x of the i-th scale best_i ; The α adaptive adjustment method is as follows: Update α = α + 10% α, that is, gradually increase the filter bandwidth, x max The center frequency of the corresponding sliding window remains unchanged, and the filter with the updated scale is used 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 optimal scale filter sub-signal x best_i ; S205: Calculate the optimal filter sub-signal x best_i The envelope harmonic product spectrum (EHPS) of the i The frequency corresponding to the location is taken as the estimated i-th fault characteristic frequency f i ; The calculation formula of EHPS is as follows: Where 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, and Π· represents the symbol of the continuous multiplication operation; S206: If i>1, calculate x best_i The EHPS decline rate r i Otherwise, return to step S202 and update x(t)=x(t)-x best_i , i=i+1; r i The calculation formula 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=[f1, f2,..., f index ], the target modal number K = index is obtained, and the fault characteristic frequency estimation value required for each modal decomposition is freq = [f1, f2, ..., f index ]; otherwise, return to step S202 and 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; Wherein, in step S3, the characteristic mode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function comprises the following steps: S301: Load the original signal x=x(t), and use the target modal number K and fault characteristic frequency freq obtained in step S207 as input parameters for characteristic mode decomposition; S302: Set the filter length L, use Q Hanning windows to initialize the FIR filter group, and obtain the initial FIR filter group {h1, h2, ..., h Q }; Upper and lower cutoff frequencies of FIR filter h l 、h u It can be expressed as: In the formula, f s is the sampling frequency of the original signal, q is the filter number, q = 1, 2, ..., Q; S303: The original signal x is filtered using a deconvolution method to obtain a decomposed mode, and the second-order cyclostationarity (ICS2) of the decomposed mode is maximized to constrain the decomposition. The expression of ICS2 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, and u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L+1, j represents a complex number unit, |·| represents an absolute value operator, and ||·|| represents a norm operator; The expression of ICS2 is rewritten into the following matrix form: E=[e1,...,e r ,...,have been R ] And r =[and -j2πrβ(L-1) ,...,And -j2πrβ(N-1) ] H Where D represents u q The diagonal matrix, u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L+1, diag(·) represents the diagonal matrix function in the MATIAB toolkit; E represents the diagonal matrix consisting of R matrices e r The matrix is ​​formed by horizontal splicing in sequence, r is the sample index, r = 1, 2, ..., R, R is a positive integer, and 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 ICS2 is equivalent to solving the following constrained problem: In the formula, h q [l] is the qth filter, whose length is L, l = 1, 2, ..., L; It means to find the filter h that maximizes ICS2. q [l]; The deconvolution method is used to iteratively filter and obtain the decomposed 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 qth filter, q = 1, 2, ..., Q, X represents the matrix of the original signal after being cut into equal lengths, the shape of which is (N-L+1) × L, L is the filter length, and N represents the number of sampling points of the original signal; Then, the decomposed modal ICS2 can be defined as: In the formula, R XwX and R XX denote the weighted correlation matrix and the correlation matrix respectively, W is the mediation matrix controlling the weighted correlation matrix; Calculate the ICS2 difference between the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than η, the algorithm is considered to have converged and the iteration is stopped, that is: Among them, η is the convergence condition, and the recommended value of η is 10 -6 , y is the current iteration number; In order to obtain the optimal filter when ICS2 is maximized, it is equivalent to solving the eigenvalue problem of the following eigenvalue problem related to the maximum eigenvalue λ: R XwX h q =R XX h q λ In the formula, solving the above formula can get λ, λ represents the maximum eigenvalue, that is, the optimal ICS2; Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposed mode u q ; S304: Mode selection; calculation of decomposition mode u q ICS2, select the largest mode of ICS2 as the kth final decomposition mode u k , then update the original signal x = xu k , the decomposition is terminated when the number of decomposed modes reaches the target number of modes K obtained in step S207.

2. A rotating machinery composite fault diagnosis system based on adaptive characteristic mode decomposition, characterized in that: include: The signal acquisition module is used to collect the composite fault signal x(t) of the rotating machinery to be diagnosed, where t is the sampling time; A parameter determination module is used to locate the resonance frequency band caused by the fault based on an adjustable scale sliding window filtering method with envelope autocorrelation kurtosis (EAK) as an indicator, so as to adaptively determine the number of modes of characteristic mode decomposition and the fault characteristic frequency corresponding to each mode; A signal decomposition module is used to adaptively decompose the collected composite fault signal into several modal components based on the characteristic mode decomposition method with the second-order cyclostationarity index (ICS2) as the objective function; Fault diagnosis module, used to calculate the envelope spectrum of decomposed modes, to achieve accurate extraction of rotating machinery fault features and fault diagnosis; Wherein, the parameter determination module includes: The sliding window filter unit is used to construct a bandpass filter bank to perform sliding window filtering on the original signal x(t) to obtain a set of filtered sub-signals {x1, x2, ..., x END }; The target screening unit is used to screen out the filter sub-signal x with the largest envelope autocorrelation kurtosis (EAK) value max , which 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 Where X(t) represents the Hilbert envelope of x(t), hilbert(·) is the Hilbert transform function in the MATLAB toolkit; R XX (τ) represents the autocorrelation function of X(t), τ 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 function in the MATLAB toolkit, ∫· represents the integral symbol, ∑· represents the summation symbol, R XX (ii) indicates R XX the ii-th value of (τ), ii=1, 2, ..., N / 2; Filter scale adjustment unit, used to adaptively adjust x max The scale factor α value of the corresponding sliding window is set to cover the entire resonance frequency band excited by the fault as much as possible, and then update x max So as to obtain the optimal filtered signal x of the i-th scale best_i ; The α adaptive adjustment method is as follows: Update α = α + 10% α, that is, gradually increase the filter bandwidth, x max The center frequency of the corresponding sliding window remains unchanged, and the filter with the updated scale is used 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 optimal scale filter sub-signal x best_i ; Frequency estimation unit, used to calculate the optimal filter sub-signal x best_i EHPS, select the EHPS global maximum value v i The frequency corresponding to the location is taken as the estimated i-th fault characteristic frequency f i ; The calculation formula of EHPS is as follows: Where 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, and Π· represents the symbol of the continuous multiplication operation; The target mode number determination unit is used to locate the number of resonant frequency bands in the signal, thereby determining the number of modes to be decomposed; if i>1, calculate x best_i The EHPS decline rate r 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; r i The calculation formula 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=[f1, f2,..., f index ], the target modal number K = index is obtained, and the fault characteristic frequency estimation value required for each modal decomposition is freq = [f1, f2, ..., fi ndex ]; otherwise, return to the sliding window filter unit and update x(t) = x(t) - x best_i , i = i + 1; when i ≤ 3, return to the sliding window filter unit and update x(t) = x(t) - x best_i , i=i+1; Wherein, the signal decomposition module comprises: The filter initialization unit is used to initialize the FIR filter; the filter length L is set, and the FIR filter group is initialized using Q Hanning windows to obtain the initial FIR filter group {h1, h2, ..., h Q }; Upper and lower cutoff frequencies of FIR filter h l 、h u It can be expressed as: In the formula, f s is the sampling frequency of the original signal, q is the filter number, q = 1, 2, ..., Q; The objective function unit is used to establish the objective function ICS2 of the characteristic mode 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, and u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L+1, j represents a complex number unit, |·| represents an absolute value operator, and ||·|| represents a norm operator; The expression of ICS2 can be rewritten into the following matrix form: E=[e1,...,e r ,...,have been R ] And r =[and -j2πrβ(L-1) ,...,And -j2πrβ(N-1) ] H Where D represents u q The diagonal matrix, u q [n] represents the matrix u q The nth element of , n = 1, 2, ..., N-L+1, diag(·) represents the diagonal matrix function in the MATIAB toolkit; E represents the diagonal matrix consisting of R matrices e r The matrix is ​​formed by horizontal splicing in sequence, r is the sample index, r = 1, 2, ..., R, R is a positive integer, and 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 ICS2 is equivalent to solving the following constrained problem: In the formula, h q [l] is the qth filter, whose length is L, l = 1, 2, ..., L; It means to find the filter h that maximizes ICS2. q [l]; The deconvolution iteration unit is used to obtain the decomposed mode by iterative filtering using the deconvolution method. 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 qth filter, q = 1, 2, ..., Q, X represents the original signal after being cut into equal lengths, and the shape is a matrix of (N-L+1)×L, L is the filter length, and N represents the number of sampling points of the original signal; Then, the decomposed modal ICS2 can be defined as: In the formula, R XWX and R XX denote the weighted correlation matrix and the correlation matrix respectively, W is the mediation matrix controlling the weighted correlation matrix; Calculate the ICS2 difference between the modes obtained in two adjacent iterations to determine whether to stop the iteration. If the difference is not greater than η, the algorithm is considered to have converged and the iteration is stopped, that is: Among them, η is the convergence condition, and the recommended value of η is 10 -6 , y is the current iteration number; In order to obtain the optimal filter when ICS2 is maximized, it is equivalent to solving the eigenvalue problem of the following eigenvalue problem related to the maximum eigenvalue λ: R XwX h q =R XX h q λ In the formula, solving the above formula can get λ, λ represents the maximum eigenvalue, that is, the optimal ICS2; Update the filter with the eigenvector corresponding to the maximum eigenvalue λ to obtain the decomposed mode u q ; Mode selection unit, used to select the optimal mode; calculate the decomposed mode u q ICS2, select the largest mode of ICS2 as the kth final decomposition mode u k , then update the original signal x = xu k , when the number of decomposed modes reaches the target mode number K, the decomposition is terminated.

Images