Adaptive signal denoising decomposition method for compound fault identification of mechanical transmission system

Abstract

The present application relates to the adaptive signal denoising decomposition method for mechanical transmission system compound fault identification, belongs to the mechanical vibration signal processing and fault diagnosis technical field. Adopting the autoregressive model to remove the harmonic component of vibration signal adaptively, through the construction filter group decomposition AR denoising vibration signal, obtains a series of modes, estimates the fault period of mode, and with the help of blind deconvolution theory, adaptive iteration updates filter and mode, carries out mode selection according to multidomain correlation coefficient and fault period consistency coefficient, determines the optimal filter length according to the weighted square envelope harmonic noise ratio, carries out square envelope spectrum analysis to the optimal mode obtained after adaptive denoising decomposition, and finally realizes compound fault identification. The present application can carry out adaptive denoising decomposition to the vibration signal of mechanical transmission system, does not need to construct priori base function and fault period priori knowledge, and significantly improves the accuracy of compound fault identification.

Description

Technical Field

[0001] This invention belongs to the field of mechanical vibration signal processing and fault diagnosis technology, and more specifically, relates to an adaptive signal denoising decomposition method for identifying complex faults in mechanical transmission systems. Background Technology

[0002] Mechanical transmission systems are a crucial component of high-end equipment (such as aircraft engines, helicopters, rail trains, and wind turbines), and their operational reliability and safety have always been a focus of attention for academic and engineering communities both domestically and internationally. Gears, bearings, and other rotating components are key fundamental parts of mechanical transmission systems, often facing harsh and complex service conditions such as heavy loads, high speeds, and high temperatures in actual operation, making them highly susceptible to failure. Extensive engineering practice shows that mechanical failures often manifest as multiple coexisting faults, i.e., complex faults that occur simultaneously or in cascades, interconnected, and mutually influencing each other. Therefore, conducting condition monitoring and complex fault diagnosis of mechanical transmission systems is of great significance for ensuring the service safety of mechanical equipment.

[0003] Vibration signals are currently the most sensitive indicators for characterizing abnormalities and faults in mechanical equipment. They have good interpretability and are simple to test, and analysis methods based on vibration signals have been widely used in the field of mechanical fault diagnosis. However, as mechanical equipment becomes increasingly complex, the vibration signal components of mechanical transmission systems are becoming more complex and exhibiting strong non-stationarity and strong nonlinearity, posing a significant challenge to traditional time-domain statistical feature analysis and spectrum analysis methods.

[0004] Decomposing vibration signals into several sub-signals (i.e., modes) with clear physical meaning is an effective way to extract fault features and identify faults in complex mechanical transmission systems. Many methods, such as wavelet transform, empirical mode decomposition (EMD), local mean mode decomposition (LMD), empirical wavelet transform, variational mode decomposition (VMD), and eigenmode decomposition (EMD), have been applied in the field of mechanical fault diagnosis and have yielded fruitful results. However, wavelet transform and empirical wavelet transform are limited by the construction of wavelet bases and cannot achieve adaptive signal decomposition; empirical mode decomposition and LMD, due to their recursive iterative decomposition methods, cannot avoid end-point effects and mode aliasing; although variational mode decomposition uses a non-recursive decomposition method, its parameters (number of modes and equilibrium parameters) need to be selected manually in advance, lacking parameter adaptability. Eigenmode decomposition, as a recently proposed signal decomposition method, fully considers the impact and cyclic stationarity of mechanical faults and is more suitable for processing mechanical fault vibration signals. However, Eigenmode decomposition also relies on human experience to pre-set parameters (filter length and number of modes), and inappropriate parameter selection can significantly reduce its decomposition performance. Furthermore, this method is prone to missing faults under strong harmonic interference generated by rotating components such as gears and shafts, making it unable to effectively and accurately identify complex faults. Therefore, there is an urgent need to propose an adaptive signal denoising decomposition method that can fully consider the characteristics of vibration signals from complex faults in mechanical transmission systems, thereby achieving accurate identification of complex faults in mechanical transmission systems. Summary of the Invention

[0005] This invention proposes an adaptive signal denoising decomposition method for identifying complex faults in mechanical transmission systems. It overcomes the shortcomings of existing signal decomposition methods, which do not fully consider the characteristics of vibration signals in complex faults of mechanical transmission systems (multiple impacts, multiple cyclic steady periods, strong harmonic interference, etc.) and require parameters to be selected in advance based on experience, which leads to reduced decomposition performance. This method can more effectively identify complex faults in mechanical transmission systems and, while satisfying the adaptive decomposition characteristics, greatly avoids the problem of missed or misdiagnosed complex faults.

[0006] The technical solution adopted in this invention is:

[0007] An adaptive signal denoising decomposition method for identifying complex faults in mechanical transmission systems, characterized by comprising the following steps:

[0008] Step 1: Input vibration signal data: Install the vibration acceleration sensor on the outer shell of the transmission device of the mechanical equipment under test, and collect the original vibration signal of the mechanical transmission system, which is marked as x(n), where n = 1, 2, 3, ..., N is the sequence of sampling points of the vibration signal, and N is the total number of sampling points, i.e. the signal length;

[0009] Step 2: Remove harmonic interference: Use an autoregressive (AR) model to estimate the harmonic components in the original vibration signal x(n), and calculate the AR-denoised vibration signal e(n) after removing harmonic interference according to formula (1) to achieve adaptive noise reduction processing of the original vibration signal.

[0010]

[0011] Where a(p) is the p-th order coefficient in the AR model, and q is the order of the AR model;

[0012] Step 3: Construct and initialize the filter bank: Construct a filter bank {f} using K finite-length unit impulse response (FIR) filters. k ,k=1,2,3,…,K}, where f k Let K be the vector representation of the k-th filter. The filter bank is specifically constructed as follows: K filters with equal bandwidth are evenly distributed across the entire analysis bandwidth, and the bandwidths of adjacent filters have a 50% overlap. The filter bank is initialized using a Hanning window function, and the lengths L of the K filters in the filter bank are set to be equal, with the filter length range L∈[1,L...]. upper ];

[0013] Step 4: Filter and decompose the AR noise-reduced vibration signal e(n) to obtain the mode u. k Apply the filter bank {f} described in step 3 k The AR noise-reduced vibration signal e(n) is filtered, and K modes {u} are obtained by eigenvector decomposition according to formula (2). k ,k=1,2,3,…,K:

[0014] u k =Εf k (2)

[0015] Among them, u k Let f be the vector representation of the k-th mode. k Let be the vector representation of the k-th filter in the filter bank, and E be the Toeplitz matrix form of the AR noise-reduced vibration signal constructed for the eigenvector method solution, as shown below:

[0016]

[0017] Step 5: Estimate the mode u k Failure cycle T k Mode u is estimated iteratively using Hilbert transform and autocorrelation function. k Failure cycle T k ;

[0018] Step 6: Adaptively iteratively update filter f k and mode u k Using blind deconvolution theory, with modal u k Relevance kurtosis maximization as filter f k The adaptive iterative update criterion is based on the AR-denoised vibration signal e(n) and the fault period T described in step 5. k mode u k The relevant kurtosis is constructed into a generalized Rayleigh quotient form, and the filter f is adaptively updated by solving the generalized eigenvalue problem. k and mode u k ;

[0019] Step 7: Determine the filter f k If the number of iterations has reached the specified number of iterations pre_iteration, then obtain the final updated mode u. k Otherwise, repeat steps 4 to 6.

[0020] Step 8, Mode Selection: First, calculate the squared envelope harmonic noise ratio, correlation kurtosis, and mean correlation kurtosis of each mode obtained from the adaptive decomposition. Remove modes with correlation kurtosis below the mean and modes with squared envelope harmonic noise ratio below the threshold SEHNR_limit. Then, calculate the multi-domain correlation coefficient (i.e., i ...

[0021] Step 9: Calculate the weighted squared envelope harmonic noise ratio for retaining the residual modes, and determine whether the filter length L has reached the upper limit of the specified length range. upper If yes, proceed to step 10; otherwise, repeat steps 4 through 8.

[0022] Step 10: Obtain the optimal filter length and optimal mode: Select the filter length corresponding to the maximum weighted square envelope harmonic noise ratio of the mode as the optimal filter length, and the corresponding mode is the optimal mode.

[0023] Step 11: Separation, extraction and identification of composite fault features: Perform square envelope spectrum analysis on the optimal mode, i.e., the sub-signal, obtained after adaptive noise reduction decomposition, and transform all sub-signals from the time domain to the envelope spectrum domain to obtain the frequency distribution of all sub-signals on the square envelope spectrum. Separate and extract the feature frequency of the composite fault, and compare it with the rotation feature frequency of the key components of the mechanical transmission system to finally detect the composite fault and identify its location.

[0024] Furthermore, in step 1 of the present invention, the vibration acceleration sensor is usually arranged on the surface of the housing at or near the rolling bearing of the mechanical equipment transmission system, and the transmission system components of the tested equipment include rolling bearings and gears;

[0025] Furthermore, in step 2 of this invention, the coefficients a(p) of the AR model are calculated by solving the Yule–Walker equation using the Levinson-Durban recursive algorithm, and the order q is adaptively determined according to formula (4), that is, by calculating the kurtosis value of the AR noise-reduced vibration signal e(n) under different orders, the order corresponding to the maximum kurtosis value is automatically selected as the optimal order:

[0026]

[0027] Among them, e q (n) represents the AR noise-reduced vibration signal with order q. It is e q The mean of (n) and the order q take values ​​in the range q∈[q... lower q upper ].

[0028] Furthermore, in step 3 of this invention, the lower cutoff frequency f of the k-th filter bandwidth... lower,k and the upper cutoff frequency f upper,k The calculation method is as follows:

[0029]

[0030] Among them, f s It is the sampling frequency of the original vibration signal.

[0031] Furthermore, step 5 of the present invention specifically includes the following steps:

[0032] Step 5.1: Calculate u for each mode. k The square envelope signal u k,SE :

[0033] u k,SE (n)=|u k (n)+j·Hilbert{u k(n)}| 2 (6)

[0034] Here, Hilbert{·} denotes the Hilbert transform, and j is the imaginary unit.

[0035] Step 5.2: Calculate the squared envelope signal u for each mode. k,SE The autocorrelation function R k,SE (τ):

[0036]

[0037] Where τ is the offset, N is the total number of signal sampling points, and L is the filter length.

[0038] Step 5.3: Select the autocorrelation function R k,SE The time sequence position τ corresponding to the maximum value after the zero-crossing point in (τ). max As mode u k Failure cycle T k :

[0039]

[0040] Furthermore, step 6 of the present invention specifically includes the following steps:

[0041] Step 6.1: Construct mode u using blind deconvolution theory. k Relevant kurtosis of the generalized Rayleigh quotient form:

[0042]

[0043]

[0044] Among them, CK(u k ) represents mode u k The correlation kurtosis, M is the shift coefficient of the correlation kurtosis, and W M R is a weighted matrix. XWX For the weighted correlation matrix, R XX Here, H is the correlation matrix, and H is the conjugate transpose symbol.

[0045] Step 6.2: Solve the generalized eigenvalue problem:

[0046] R XWX f k =R XX f k λ (11)

[0047] Where λ is the generalized eigenvalue, and its maximum value is λ. max With mode u k The maximum correlated kurtosis value corresponds to.

[0048] Step 6.3: Adaptively update filter f k and mode u k Applications and the maximum generalized eigenvalue λ max The corresponding eigenvector adaptive update filter f k Thus, the updated mode u is obtained. k .

[0049] Furthermore, step 8 of the present invention specifically includes the following steps:

[0050] Step 8.1: Calculate the squared envelope harmonic noise ratio (SEHNR) and correlation kurtosis (CK) for each mode obtained from the decomposition:

[0051]

[0052]

[0053] Among them, R k,SE (τ max ) is the mode u k The maximum value of the autocorrelation function after crossing zero, R k,SE (0) represents the mode u k The value of the zero point of the autocorrelation function, M is the shift coefficient of the correlation kurtosis, and N is the signal length.

[0054] Step 8.2: Remove modes whose squared envelope harmonic noise ratio is lower than the threshold SEHNR_limit and whose correlation kurtosis value is lower than the mean correlation kurtosis value of all modes;

[0055] Step 8.3: Calculate the multi-domain correlation coefficient (i.e., the i ... T :

[0056]

[0057]

[0058]

[0059] ε T =|T p -T q | / T p (17)

[0060] Among them, u p and u q These are the p-th and q-th modes, respectively. and S represents the mean of the p-th and q-th modes, respectively.p and S q Mode u p and u q The spectrum, SES p and SES q These are the modes u p and u q The square envelope spectrum, T p and T q These are the modes u p and u q The estimated failure cycle.

[0061] Step 8.4: Determine the multi-domain correlation coefficient and fault cycle consistency coefficient between every two modes, remove modes with a multi-domain correlation coefficient threshold higher than Similarity_limit and a fault cycle consistency coefficient lower than the period_limit, and retain the remaining modes.

[0062] Furthermore, in step 10 of this invention, the optimal filter length is calculated as follows:

[0063] Step 10.1: Calculate the final remaining K when the filter length is L. d Weighted squared envelope harmonic noise ratio of each mode. L :

[0064]

[0065] Among them, SEHNR k K represents the square envelope harmonic noise ratio of the k-th residual mode. d This represents the final number of remaining modes.

[0066] Step 10.2: Calculate the final remaining K for different filter lengths. d The weighted squared envelope harmonic noise ratio of the modes is used to automatically select the filter length corresponding to the maximum weighted squared envelope harmonic noise ratio as the optimal filter length.

[0067]

[0068] Where L is the filter length, and its value ranges from [1, L]. upper ], L opt This is the optimal filter length.

[0069] The beneficial effects of this invention are as follows:

[0070] 1. The adaptive signal denoising decomposition method for identifying complex faults in mechanical transmission systems proposed in this invention combines autoregressive (AR) models and blind deconvolution theory with the idea of ​​adaptive signal decomposition. By using an AR model to adaptively remove strong harmonic components generated by gear and shaft rotation, the interference of strong harmonic components in vibration signals on fault information is greatly avoided. This is beneficial for fully exploring and utilizing the characteristic of correlation kurtosis in representing both the impact and cyclic stability of fault signals. At the same time, by using a filter bank with frequency band overlap and a blind deconvolution technique with the correlation kurtosis of the fault signal as the objective function, the method can achieve adaptive decomposition of vibration signals guided by the characteristics of complex faults without prior knowledge of the multiple cyclic stationary periods (i.e., fault periods) of the complex faults or the construction of prior basis functions. This enables effective identification of complex faults in mechanical transmission systems.

[0071] 2. The weighted square envelope harmonic noise ratio (SHNR) proposed in this invention fully considers the characteristics of multiple impacts and multiple cyclic steady periods of composite fault vibration signals, and can effectively evaluate the fault information of all decomposed modes, making it suitable for composite fault feature detection.

[0072] 3. The multi-domain correlation coefficient proposed in this invention includes the time-domain correlation coefficient and the spectral orthogonality coefficient and square envelope spectrum orthogonality coefficient proposed in this invention. It can comprehensively examine the correlation of decomposed modes from three dimensions: time domain, frequency domain and envelope spectrum domain. It overcomes the shortcomings of traditional methods that rely solely on the time-domain correlation coefficient as the only correlation evaluation means, which cannot fully consider the characteristics of fault vibration signals and has low accuracy. Combined with the fault period consistency coefficient proposed in this invention, it can effectively select the optimal mode with rich fault feature information and avoid the problems of modal redundancy or loss of fault feature information caused by over-decomposition or under-decomposition.

[0073] 4. The adaptive signal denoising decomposition method for identifying complex faults in mechanical transmission systems proposed in this invention does not require prior fault knowledge for the main parameters involved, namely the AR model order, filter length, and number of modes. These parameters can be adaptively selected based on the proposed kurtosis-based adaptive determination criterion for the AR model order, the weighted square envelope harmonic noise ratio index, the multi-domain correlation coefficient, and the fault period consistency coefficient. This overcomes the shortcomings of traditional signal decomposition methods that require manual selection of parameters in advance based on experience, effectively avoids the problem of missed or misdiagnosed faults caused by improper parameter selection, and significantly improves the accuracy of complex fault diagnosis. Attached Figure Description

[0074] Figure 1 This is a flowchart of the method proposed in this invention;

[0075] Figure 2This is a time-domain waveform diagram of the combined fault vibration signal of the rolling bearing in an embodiment of the present invention;

[0076] Figure 3 This is a spectrum diagram of the vibration signal of a rolling bearing with a combined fault in an embodiment of the present invention;

[0077] Figure 4 The image shows the time-domain waveform of the rolling bearing composite fault vibration signal obtained after noise reduction by the method proposed in this invention in an embodiment of the invention.

[0078] Figure 5 The spectrum diagram of the vibration signal of the rolling bearing composite fault in the embodiment of the present invention after noise reduction by the method proposed in the present invention is shown.

[0079] Figure 6 This is a time-domain waveform of the rolling bearing composite fault vibration signal obtained after noise reduction and decomposition by the method proposed in this invention in an embodiment of the present invention.

[0080] Figure 7 This is the square envelope spectrum of the rolling bearing composite fault vibration signal obtained after noise reduction and decomposition by the method proposed in this invention in an embodiment of the invention. Detailed Implementation

[0081] The vibration signals of the mechanical transmission system used in this embodiment of the invention are derived from the XJTU-SY rolling bearing accelerated life test dataset publicly released by Xi'an Jiaotong University. The rolling bearing vibration signals are obtained by a vibration acceleration sensor vertically arranged on the bearing housing, with a sampling frequency of 25600Hz and a sampling duration of 1.28s. The input rotational speed of the tested rolling bearing shaft is 2250rpm (corresponding to a rotational frequency of 37.5Hz), the radial load is 11kN, and the fault mode is the simultaneous occurrence of outer ring fault and cage fault, with fault characteristic frequencies of f1 = 116.41Hz and f2 = 14.84Hz, respectively.

[0082] Figure 1 The figure shows an adaptive signal denoising decomposition method for complex fault identification in mechanical transmission systems proposed in this invention, characterized by the following steps:

[0083] Step 1: Input Vibration Signal Data: Install the vibration accelerometer on the outer casing of the transmission device of the mechanical equipment under test to collect the original vibration signal of the mechanical transmission system, labeled as x(n), where n = 1, 2, 3, ..., N is the sequence of sampling points of the vibration signal, and N is the total number of sampling points, i.e., the signal length. In this embodiment, the total number of sampling points N = 32768. The time-domain waveform of the vibration signal is as follows: Figure 2 As shown, the spectrum of the vibration signal is as follows: Figure 3 As shown;

[0084] Step 2: Harmonic interference removal: The harmonic components in the original vibration signal x(n) are estimated using an autoregressive (AR) model, and the AR-denoised vibration signal e(n) after harmonic interference removal is calculated according to formula (1), thus realizing adaptive noise reduction processing of the original vibration signal.

[0085]

[0086] Where a(p) is the p-th order coefficient in the AR model, which is obtained by solving the Yule–Walker equation using the Levinson-Durban recursive algorithm; q is the order of the AR model, which is adaptively determined according to formula (4), that is, by calculating the kurtosis value of the AR noise-reduced vibration signal e(n) under different orders, the order corresponding to the maximum kurtosis value is automatically selected as the optimal order:

[0087]

[0088] Among them, e q (n) represents the AR noise-reduced vibration signal with order q. It is e q The mean of (n) and the order q take values ​​in the range q∈[q... lower q upper In this embodiment, the order q ranges from [20, 200]; the time-domain waveform of the AR noise-reduced vibration signal is as follows: Figure 4 As shown, the spectrum of the AR noise-reduced vibration signal and its local magnified image are as follows. Figure 5 As shown, with Figure 2 and Figure 3 Compared to the original vibration signal shown, after AR noise reduction, the harmonic components of the signal in the frequency range of 0-3000Hz are significantly removed (the spectral amplitude is reduced from...). Figure 3 The 0.3g shown decreased to Figure 5 The 0.03g shown has a time-domain amplitude of [value missing]. Figure 2 The 2g shown decreased to Figure 4 (as shown in 1g);

[0089] Step 3: Construct and initialize the filter bank: Construct a filter bank {f} using K finite-length unit impulse response (FIR) filters. k ,k=1,2,3,…,K}, where f k Let K be the vector representation of the k-th filter. The filter bank is specifically constructed as follows: K filters with equal bandwidth are evenly distributed across the entire analysis bandwidth, and the bandwidths of adjacent filters have a 50% overlap. The filter bank is initialized using a Hanning window function, and the lengths L of the K filters in the filter bank are set to be equal, with the filter length range L∈[1,L...].upper The lower cutoff frequency f of the k-th filter bandwidth. lower,k and the upper cutoff frequency f upper,k The calculation method is as follows:

[0090]

[0091] Among them, f s This is the sampling frequency of the original vibration signal. In this embodiment, the number of filters in the filter bank K = 11, the filter length L ranges from [1, 200], and the sampling frequency f is... s =25600Hz;

[0092] Step 4: Decompose the AR-denoised vibration signal e(n) to obtain the mode u k Apply the filter bank {f} described in step 3 k The AR noise-reduced vibration signal e(n) is filtered, and K modes {u} are obtained by eigenvector decomposition according to formula (2). k ,k=1,2,3,…,K:

[0093] u k =Εf k (2)

[0094] Among them, u k Let f be the vector representation of the k-th mode. k Let be the vector representation of the k-th filter in the filter bank, and E be the Toeplitz matrix form of the AR noise-reduced vibration signal constructed for the eigenvector method solution, as shown below:

[0095]

[0096] Step 5: Estimate the mode u k Failure cycle T k Mode u is estimated iteratively using Hilbert transform and autocorrelation function. k Failure cycle T k ;

[0097] Step 5.1: Calculate u for each mode. k The square envelope signal u k,SE :

[0098] u k,SE (n)=|u k (n)+j·Hilbert{u k (n)}| 2 (6)

[0099] Here, Hilbert{·} denotes the Hilbert transform, and j is the imaginary unit.

[0100] Step 5.2: Calculate the squared envelope signal u for each mode. k,SE The autocorrelation function R k,SE (τ):

[0101]

[0102] Where τ is the offset, N is the total number of signal sampling points, and L is the filter length.

[0103] Step 5.3: Select the autocorrelation function R k,SE The time sequence position τ corresponding to the maximum value after the zero-crossing point in (τ). max As mode u k Failure cycle T k :

[0104]

[0105] Step 6: Adaptively iteratively update filter f k and mode u k Using blind deconvolution theory, with modal u k Relevance kurtosis maximization as filter f k The adaptive iterative update criterion is based on the AR-denoised vibration signal e(n) and the fault period T described in step 5. k mode u k The relevant kurtosis is constructed into a generalized Rayleigh quotient form, and the filter f is adaptively updated by solving the generalized eigenvalue problem. k and mode u k ;

[0106] Step 6.1: Construct mode u using blind deconvolution theory. k Relevant kurtosis of the generalized Rayleigh quotient form:

[0107]

[0108]

[0109] Among them, CK(u k ) represents mode u k The correlation kurtosis, M is the shift coefficient of the correlation kurtosis, and W M R is a weighted matrix. XWX For the weighted correlation matrix, R XX Here, H is the correlation matrix, and H is the conjugate transpose symbol. In this embodiment, the shift coefficient M of the correlation kurtosis is 1.

[0110] Step 6.2: Solve the generalized eigenvalue problem:

[0111] R XWX fk =R XX f k λ (11)

[0112] Where λ is the generalized eigenvalue, and its maximum value is λ. max With mode u k The maximum correlated kurtosis value corresponds to;

[0113] Step 6.3: Adaptively update filter f k and mode u k Applications and the maximum generalized eigenvalue λ max The corresponding eigenvector adaptive update filter f k Thus, the updated mode u is obtained. k .

[0114] Step 7: Determine the filter f k If the number of iterations has reached the pre-specified number of iterations (pre_iteration), then the final updated mode u is obtained. k Otherwise, repeat steps 4 to 6. In this embodiment, the number of iterations pre_iteration = 8.

[0115] Step 8, Mode Selection: First, calculate the squared envelope harmonic noise ratio, correlation kurtosis, and mean correlation kurtosis of each mode obtained from the adaptive decomposition. Remove modes with correlation kurtosis below the mean and modes with squared envelope harmonic noise ratio below the threshold SEHNR_limit. Then, calculate the multi-domain correlation coefficient (i.e., i ...

[0116] Step 8.1: Calculate the squared envelope harmonic noise ratio (SEHNR) and correlation kurtosis (CK) for each mode obtained from the decomposition:

[0117]

[0118]

[0119] Among them, R k,SE (τ max ) is the mode u k The maximum value of the autocorrelation function after crossing zero, Rk,SE (0) represents the mode u k The value of the zero point of the autocorrelation function, M is the shift coefficient of the correlation kurtosis, and N is the signal length.

[0120] Step 8.2: Remove modes whose squared envelope harmonic noise ratio is lower than the threshold SEHNR_limit and whose correlation kurtosis value is lower than the mean correlation kurtosis value of all modes;

[0121] Step 8.3: Calculate the multi-domain correlation coefficient (i.e., the i ... T :

[0122]

[0123]

[0124]

[0125] ε T =|T p -T q | / T p (17)

[0126] Among them, u p and u q These are the p-th and q-th modes, respectively. and S represents the mean of the p-th and q-th modes, respectively. p and S q Mode u p and u q The spectrum, SES p and SES q These are the modes u p and u q The square envelope spectrum, T p and T q These are the modes u p and u q The estimated failure cycle.

[0127] Step 8.4: Determine the multi-domain correlation coefficient and fault cycle consistency coefficient between every two modes, remove modes with a multi-domain correlation coefficient threshold higher than Similarity_limit and a fault cycle consistency coefficient lower than the period_limit, and retain the remaining modes.

[0128] Step 9: Calculate the weighted squared envelope harmonic noise ratio for retaining the residual modes, and determine whether the filter length L has reached the upper limit of the pre-specified length range L. upperIf yes, proceed to step 10; otherwise, repeat steps 4 through 8. In this embodiment, the upper limit of the filter length range is L. upper =200;

[0129] Step 10: Obtain the optimal filter length and optimal mode: Select the filter length corresponding to the maximum weighted square envelope harmonic noise ratio of the mode as the optimal filter length, and the corresponding mode is the optimal mode.

[0130] Step 10.1: Calculate the final remaining K when the filter length is L. d Weighted squared envelope harmonic noise ratio of each mode. L :

[0131]

[0132] Among them, SEHNR k K represents the square envelope harmonic noise ratio of the k-th residual mode. d This represents the final number of remaining modes.

[0133] Step 10.2: Calculate the final remaining K for different filter lengths. d The weighted squared envelope harmonic noise ratio of the modes is used to automatically select the filter length corresponding to the maximum weighted squared envelope harmonic noise ratio as the optimal filter length.

[0134]

[0135] Where L is the filter length, and its value ranges from [1, 200]. opt To determine the optimal filter length, in this embodiment, the solution is: the final number of modes K obtained through decomposition. d =2, optimal filter length L opt =156;

[0136] Step 11, Composite Fault Feature Separation, Extraction, and Identification: The optimal modes (i.e., sub-signals) obtained after adaptive noise reduction decomposition are subjected to squared envelope spectrum analysis. In this embodiment, the number of optimal modes is 2, denoted as mode #1 and mode #2, respectively. The time-domain waveforms of mode #1 and mode #2 are as follows... Figure 6 As shown, the two modes are transformed from the time domain to the envelope spectral domain, and the frequency distributions of the two modes on the squared envelope spectrum are obtained, as follows. Figure 7 As shown, the characteristic frequency of the rolling bearing outer ring fault in this example is f1 = 116.41 Hz, which is consistent with... Figure 7 The fault characteristic frequency f1 and its harmonics 2x, 3x, ... contained in the middle mode #1 are consistent. The fault characteristic frequency of the rolling bearing cage is f2 = 14.84 Hz, which is consistent with... Figure 7The fault characteristic frequency f2 and its harmonics 2x, 3x, ... contained in the middle mode #2 are consistent. The two fault characteristics can be identified through modes #1 and #2 respectively. Simultaneously, the cage fault characteristic frequency modulates the outer ring fault characteristic frequency, therefore... Figure 7 In the mode #1 shown, in addition to the obvious outer ring fault characteristic frequency f1 and its harmonics, cage fault characteristic frequencies f2 and f1+f2 also appear, which are consistent with the modulation mechanism of rolling bearing fault vibration signal, and verify the effectiveness of the method proposed in this invention for detecting composite faults and identifying the location of composite faults.

Claims

1. An adaptive signal denoising decomposition method for complex fault identification in mechanical transmission systems, characterized in that, Includes the following steps: Step 1: Acquire vibration signal data: Install the vibration acceleration sensor on the outer shell of the transmission device of the mechanical equipment under test, and collect the original vibration signal of the mechanical transmission system, which is marked as x(n), where n=1, 2, 3, ..., N, which is the sequence of sampling points of the vibration signal, and N is the total number of sampling points, i.e. the signal length; Step 2: Remove harmonic interference: Use an autoregressive AR model to estimate the harmonic components in the original vibration signal x(n), and calculate the AR-denoised vibration signal e(n) after removing harmonic interference according to formula (1) to achieve adaptive noise reduction processing of the original vibration signal. (1) Where a(p) is the p-th order coefficient in the AR model, and q is the order of the AR model; Step 3: Construct and initialize the filter bank: Construct a filter bank {f} using K finite-length unit impulse response (FIR) filters. k ,k=1,2,3,…,K}, where f k Let K be the vector representation of the k-th filter. The filter bank is specifically constructed as follows: K filters with equal bandwidth are evenly distributed across the entire analysis bandwidth, and the bandwidths of adjacent filters overlap by 50%. The filter bank is initialized using a Hanning window function, and the lengths L of the K filters in the filter bank are set to be equal, with the filter length range L∈[1,L...]. upper ]; Step 4: Decompose the AR-denoised vibration signal e(n) to obtain the mode u k Apply the filter bank {f} described in step 3 k The AR noise-reduced vibration signal e(n) is filtered, and K modes {u} are obtained by eigenvector decomposition according to formula (2). k ,k=1,2,3,…,K}: (2) Among them, u k Let f be the vector representation of the k-th mode. k Let be the vector representation of the k-th filter in the filter bank, and E be the Toeplitz matrix form of the AR noise-reduced vibration signal constructed for the eigenvector method solution, as shown below: (3) Step 5: Estimate the mode u k Failure cycle T k Iterative estimation of mode u using Hilbert transform and autocorrelation function k Failure cycle T k ; Step 6: Adaptively iteratively update filter f k and mode u k Using blind deconvolution theory, with modal u k Correlation kurtosis maximization as filter f k The adaptive iterative update criterion is based on the AR-denoised vibration signal e(n) and the fault period T described in step 5. k mode u k The relevant kurtosis is constructed into a generalized Rayleigh quotient form, and the filter f is adaptively updated by solving the generalized eigenvalue problem. k and mode u k ; Step 7: Determine the filter f k If the number of iterations has reached the pre-specified number of iterations (pre_iteration), then the final updated mode u is obtained. k Otherwise, repeat steps 4 to 6. Step 8, Mode Selection: First, calculate the squared envelope harmonic noise ratio, correlation kurtosis, and mean correlation kurtosis of each mode obtained from the adaptive decomposition. Remove modes with correlation kurtosis below the mean and modes with squared envelope harmonic noise ratio below the threshold SEHNR_limit. Then, calculate the multi-domain correlation coefficient and fault period consistency coefficient between every two modes in the remaining modes. The multi-domain correlation coefficient includes the time domain correlation coefficient, spectral orthogonality coefficient, and squared envelope spectral orthogonality coefficient. Remove modes with multi-domain correlation coefficient higher than the threshold Similarity_limit and fault period consistency coefficient lower than the threshold period_limit, and retain the remaining modes. Step 9: Calculate the weighted squared envelope harmonic noise ratio for retaining the residual modes, and determine whether the filter length L has reached the upper limit of the pre-specified length range L. upper If yes, proceed to step 10; otherwise, repeat steps 4 through 8. Step 10: Obtain the optimal filter length and optimal mode: Select the filter length corresponding to the maximum weighted square envelope harmonic noise ratio of the mode as the optimal filter length, and the corresponding mode is the optimal mode. Step 11, Fault Feature Separation, Extraction and Identification: Perform square envelope spectrum analysis on the optimal mode, i.e., the sub-signal, obtained after adaptive noise reduction decomposition. Transform all sub-signals from the time domain to the frequency domain and obtain the frequency distribution of all sub-signals on the square envelope spectrum. Separate and extract the characteristic frequency of the composite fault and compare it with the rotation characteristic frequency of the key components of the mechanical transmission system. Finally, the composite fault can be detected and the location of the composite fault can be identified.

2. The adaptive signal denoising decomposition method for complex fault identification in mechanical transmission systems according to claim 1, characterized in that: In step 1, the vibration acceleration sensor is usually placed on the surface of the housing at or near the rolling bearing of the mechanical equipment transmission system. The transmission system components of the equipment under test include rolling bearings and gears.

3. The adaptive signal denoising decomposition method for composite fault identification in mechanical transmission systems according to claim 1, characterized in that: In step 2, the coefficients a(p) of the AR model are obtained by solving the Yule–Walker equation using the Levinson-Durban recursive algorithm, and the order q is adaptively determined according to formula (4). That is, by calculating the kurtosis value of the AR noise reduction vibration signal e(n) under different orders, the order corresponding to the maximum kurtosis value is automatically selected as the optimal order. (4) Among them, e q (n) represents the AR noise-reduced vibration signal with order q. It is e q The mean of (n) and the order q take values ​​in the range q∈[q... lower q upper ].

4. The adaptive signal noise reduction and decomposition method for composite fault identification in mechanical transmission systems according to claim 1, characterized in that: In step 3, the lower cutoff frequency f of the k-th filter bandwidth is... lower,k and the upper cutoff frequency f upper,k The calculation method is as follows: (5) Among them, f s It is the sampling frequency of the original vibration signal.

5. The adaptive signal noise reduction and decomposition method for composite fault identification in mechanical transmission systems according to claim 1, characterized in that: Step 5 specifically includes the following steps: Step 5.1: Calculate u for each mode. k The square envelope signal u k,SE : (6) in, This represents the Hilbert transform, where j is the imaginary unit; Step 5.2: Calculate the squared envelope signal u for each mode. k,SE autocorrelation function R k,SE (τ): (7) Where τ is the offset, N is the total number of signal sampling points, and L is the filter length; Step 5.3: Select the autocorrelation function R k,SE The time sequence position τ corresponding to the maximum value after the zero-crossing point in (τ). max As mode u k Failure cycle T k : (8)。 6. The adaptive signal noise reduction and decomposition method for composite fault identification in mechanical transmission systems according to claim 1, characterized in that: In step 10, the optimal filter length is calculated as follows: Step 10.1: Calculate the final remaining K when the filter length is L. d Weighted squared envelope harmonic noise ratio of each mode. L : (18) Among them, SEHNR k K represents the square envelope harmonic noise ratio of the k-th residual mode. d This represents the final number of remaining modes; Step 10.2: Calculate the final remaining K for different filter lengths. d The weighted squared envelope harmonic noise ratio of the modes is used to automatically select the filter length corresponding to the maximum weighted squared envelope harmonic noise ratio as the optimal filter length. (19) Where L is the filter length, and its value ranges from [1, L]. upper ], L opt This is the optimal filter length.

Images