A rotating machinery modulation feature extraction method based on priori traversal envelope spectrum

Abstract

The application discloses a rotating machinery modulation feature extraction method based on prior traversal envelope spectrum, comprising the following steps: (1) collecting and calculating the spectral coherence function of the rotating machinery monitoring signal based on Gabor transform; (2) performing a slicing operation on the spectral coherence function along the carrier frequency direction, and calculating the harmonic intensity vector of different harmonic cluster structures; (3) calculating the harmonic intensity matrix containing all carrier frequencies; (4) when the prior information is a single frequency, a single weighting function is calculated; when the prior information is a certain frequency range, a composite weighting function is calculated; (5) performing a threshold filtering operation to obtain a prior traversal weighting function; (6) performing prior traversal weighting processing on the spectral coherence function to obtain prior traversal spectral coherence, and further obtaining prior traversal envelope spectrum. By using the application, interference can be adaptively filtered and signal features can be extracted and enhanced, and good extraction effect can be realized even under strong noise interference.

Description

Technical Field

[0001] This invention belongs to the field of signal processing, and in particular relates to a method for extracting rotating mechanical modulation features based on prior ergodic envelope spectrum. Background Technology

[0002] Characteristic frequency extraction is widely used in the diagnosis of rotating machinery. It achieves fault diagnosis by separating the characteristic frequencies corresponding to various faults in rotating machinery from monitoring signals such as vibration signals. Characteristic frequency extraction can be approached from time domain, frequency domain, or time-frequency domain analysis. Commonly used methods include empirical mode decomposition and wavelet transform in time-frequency domain analysis.

[0003] For example, Chinese patent document CN114354188A discloses a rotating machinery fault diagnosis method based on fully adaptive noise set empirical mode decomposition; Chinese patent document CN102539150A discloses an adaptive fault diagnosis method for rotating machinery components based on continuous wavelet transform.

[0004] Empirical Mode Decomposition (EMD) decomposes rotating machinery monitoring signals into a finite number of intrinsic mode functions (EMFs), each containing local features of the signal at different time scales. This method plots upper and lower envelopes based on the local extrema of the original signal, subtracts the mean envelope from the original signal to obtain an intermediate signal, and uses this intermediate signal as a basis to determine whether the conditions for an EMF are met. If not, the above operation is repeated based on this intermediate signal. Empirical Mode Decomposition can effectively handle non-stationary and nonlinear signals and has good adaptability to the local structure of the signal in the time and frequency domain without the need for preset basis functions. However, it is prone to end-effects and mode aliasing, affecting the quality and reliability of the decomposition results. Furthermore, its reliance on empirical parameter selection also contributes to the problem.

[0005] Wavelet transform is another commonly used method for feature extraction. It replaces infinitely long trigonometric basis functions with finite-length, decaying wavelet basis functions, thus overcoming the Fourier transform's poor handling of abrupt and non-stationary signals. Wavelet coefficients are obtained by scaling and translating the wavelet basis functions, multiplying them with the original signal, and integrating. This allows for the determination of signal variations at different times and frequencies. However, this method lacks a unified standard for selecting wavelet basis functions and is still limited by Heisenberg's uncertainty principle, meaning that time resolution and frequency resolution cannot be simultaneously optimized, which affects the quality of the results. Summary of the Invention

[0006] This invention provides a method for extracting modulation features of rotating machinery based on prior traversal of the envelope spectrum, applicable to feature extraction of rotating machinery monitoring signals with prior knowledge of characteristic frequencies. This method can adaptively filter interference and extract and enhance signal features, achieving good extraction results even under strong noise interference, and has promising applications in rotating machinery condition monitoring and fault diagnosis.

[0007] A method for extracting rotational mechanical modulation features based on prior ergodic envelope spectrum includes the following steps:

[0008] (1) Collect vibration or noise data of rotating machinery as monitoring signals, and calculate the spectral coherence function of the monitoring signals based on Gabor transform;

[0009] (2) The spectral coherence function is sliced ​​along the carrier frequency direction to capture the harmonic cluster structure of different fundamental frequencies and evaluate the overall strength of the harmonic cluster structure to obtain the harmonic intensity vector.

[0010] (3) Integrate the harmonic intensity vectors corresponding to all carrier frequencies to obtain the harmonic intensity matrix;

[0011] (4) When the prior information is a single frequency, perform a single slice operation on the harmonic intensity matrix to calculate a single weighting function.

[0012] When the prior information is within a certain frequency range, the harmonic intensity matrix is ​​traversed and sliced ​​to calculate the composite weighting function.

[0013] (5) Calculate the information lower limit threshold of the weighted function, perform threshold filtering operation, and obtain the prior traversal weighted function;

[0014] (6) Perform a priori ergodic weighting on the spectral coherence function to obtain the a priori ergodic spectral coherence, and further perform absolute value integration along the carrier frequency direction to obtain the a priori ergodic envelope spectrum.

[0015] Compared with the prior art, the present invention has the following beneficial effects:

[0016] 1. This invention proposes a method for constructing a harmonic intensity matrix based on vectors obtained along the carrier frequency from spectral coherence slices. This matrix can be used as an intensity index for different spectral frequencies at various cyclic frequencies, while preserving the harmonic cluster structure information.

[0017] 2. This invention proposes a method for constructing a priori traversal weighted functions. When faced with single or range-based prior information, this method can adopt different construction strategies to obtain a unified weighted function, assigning corresponding weights based on the amount of information at each spectral frequency. If the amount of information is large, the corresponding weight value is larger; if the amount of information is small, the opposite is true.

[0018] 3. The prior ergodic envelope spectrum proposed in this invention achieves signal demodulation by performing prior ergodic weighting on the spectral coherence function. Under the influence of complex and strong noise, it can still adaptively enhance and extract the characteristics of the modulation signal of rotating machinery based on the prior information of the modulation frequency, and is suitable for fields such as condition monitoring and fault diagnosis. Attached Figure Description

[0019] Figure 1 This is a flowchart illustrating a method for extracting rotating mechanical modulation features based on prior traversal envelope spectrum according to the present invention.

[0020] Figure 2 This is a time-domain diagram of a simulated signal containing Gaussian noise, impulse noise, and cyclic stationary noise in an embodiment of the present invention.

[0021] Figure 3 This is the narrowband demodulation result of the kurtosis spectrum of the simulated signal in the embodiment of the present invention;

[0022] Figure 4 This is the cyclostationary analysis and demodulation result of the simulated signal in the embodiment of the present invention;

[0023] Figure 5 This refers to the prior traversal envelope spectrum results of the simulated signal in this embodiment of the invention;

[0024] Figure 6 This is the narrowband demodulation result of the kurtosis spectrum of the centrifugal pump vibration signal in this embodiment of the invention;

[0025] Figure 7 This is the result of cyclic stationarity analysis and demodulation of the centrifugal pump vibration signal in this embodiment of the invention;

[0026] Figure 8 This is the prior ergodic envelope spectrum result of the centrifugal pump vibration signal in this embodiment of the invention. Detailed Implementation

[0027] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be noted that the embodiments described below are intended to facilitate the understanding of the present invention and do not constitute any limitation thereof.

[0028] like Figure 1 As shown, a method for extracting rotational mechanical modulation features based on prior ergodic envelope spectrum includes the following steps:

[0029] A method for extracting rotational mechanical modulation features based on prior ergodic envelope spectrum, characterized by comprising the following steps:

[0030] S01: Collect vibration or noise data of rotating machinery as monitoring signals, and calculate the spectral coherence function of the monitoring signals based on Gabor transform.

[0031] (1-1) The sampling frequency is FS The unit is Hz, and the vibration or noise data of rotating machinery is collected as the monitoring signal x(t). n ), t n Refers to the sampling frequency F S The moment of acquisition, where t n =n / F s The sampling duration of the monitoring signal is T, and the unit is seconds.

[0032] (1-2) Calculate the monitoring signal x(t) n The short-time Fourier transform X STFT (i,f k ):

[0033]

[0034] In the formula, N w R is the window width, R is the step size, w[n] is the window function, and x[n] is the value of x(t). n abbreviation of ) f k For discrete frequencies, f k =kΔf,k=0,...,N w-1 Frequency resolution Δf = F s / N w .

[0035] (1-3) For the monitoring signal x(t) n The short-time Fourier transform X STFT (i,f k Phase correction is performed to obtain the Gabor transform result X. w (i,f k ):

[0036]

[0037] In the formula, X w (i,f k ) is the signal x(t) n In iR / F s At that moment, with f k A complex envelope centered at x with bandwidth Δf, |X w (i,f k )| 2 It represents the energy flow within the frequency band.

[0038] (1-4) Calculate the monitoring signal x(t) n The mean cyclic period spectrum is related to:

[0039]

[0040] In the formula, K = (LN w+R) / R, where R is the total number of windows that move in steps R within a signal of length L, and the length of these windows is N. w α is the cycle frequency, f is the carrier frequency, where

[0041] (1-5) Calculate the monitoring signal x(t) n ) spectral correlation function S x (α,f), the calculation formula is as follows:

[0042]

[0043] (1-6) Calculate the monitoring signal x(t) n ) spectral coherence function γ x (α, f), the calculation formula is as follows:

[0044]

[0045] S02, the spectral coherence function is sliced ​​along the carrier frequency direction to capture the harmonic cluster structure of different fundamental frequencies and evaluate the overall strength of the harmonic cluster structure, thus obtaining the harmonic intensity vector.

[0046] (2-1) Spectral correlation function γ x The discrete form of (α,f) is γ x (α m ,f n Its cycle frequency includes α. m (m = 1, 2, ..., M) has a total of M values, and the spectral frequency includes f n (n = 1, 2, ..., N) has a total of N possible values.

[0047] (2-2) Slice the spectral coherence function according to the carrier frequency direction to obtain the spectral coherence function at a specific spectral frequency f. n slice γ at the location x (α m ,f n (m=1,2,...,M);

[0048] (2-3) Take the cycle frequency α m The fundamental frequency of the harmonic structure is used to find the peak value within a range of the fundamental frequency amplitude. The calculation formula is as follows:

[0049]

[0050]

[0051] In the formula, For the fundamental frequency α m The peak-finding range, Δα is the single-sided range of the peak-finding, p m,n For the spectral frequency fn At the fundamental frequency α m Peak amplitude;

[0052] (2-4) Set the harmonic order Z of the harmonic cluster structure, and the fundamental frequency α of interest. m The peak finding of all harmonics is performed within a range, and the calculation formula is as follows:

[0053]

[0054]

[0055] In the formula, For harmonic z*α m The peak-finding range, Δα is the single-sided range of the peak-finding, p z*m,n To achieve a specific spectral frequency f n For harmonic z*α m Peak amplitude of (z=2,3,...,Z);

[0056] (2-5) Regarding the fundamental frequency α m and all its harmonics z*α m The peak amplitudes of (z = 2, 3, ..., Z) are taken as the arithmetic mean, and the spectral frequency f is calculated. n Cycle frequency α m Harmonic intensity at:

[0057]

[0058] In the formula, H x (α m f n ) represents the spectral frequency f n Cycle frequency α m The harmonic intensity at that location.

[0059] S03 integrates the harmonic intensity vectors corresponding to all carrier frequencies to obtain the harmonic intensity matrix.

[0060] Assuming the cycle frequency α m In [0, α int Within the range of ], while the spectral frequency f n In [0, F s Within the range of / 2], where α int As the upper bound of the target cyclic frequency range, when m and n increase in different directions, the harmonic intensity vector H can be... x (α m ,f n Integrate into a harmonic intensity matrix H x (α,f) mn =H x (α m ,f n).

[0061] S04. When the prior information is a single frequency, a single slice operation is performed on the harmonic intensity matrix to calculate a single weighting function.

[0062] (4-1) When the prior information is a single frequency, the single weighting function w(f) can be directly obtained by slicing the harmonic intensity matrix HIM at the cyclic frequency ω and the characteristic frequency f:

[0063] w(f) = H x (ω,f)

[0064] (4-2) For a single weighting function w(f) n Perform a normalization operation to calculate w. N (f):

[0065]

[0066] When the prior information is within a certain frequency range, the harmonic intensity matrix is ​​traversed and sliced ​​to calculate the composite weighting function.

[0067] (4-a) When the prior information is a frequency range, assume that the frequency range of the prior information features is [0, α]. int ] Calculate the mean k(f) of all slices within this range:

[0068]

[0069] In the formula, ω k ∈[0,α int ], k = 1, 2, ..., K, where K is the total number of frequencies contained within the range;

[0070] (4-b) After normalization, the composite weighted function w can be obtained. N (f):

[0071]

[0072] In the formula, max[(k(f)] n )] and min[(k(f n )] are the maximum and minimum values ​​of a single weighted function, respectively.

[0073] S05, calculate the lower information threshold of the weighted function, perform threshold filtering operation, and obtain the prior traversal weighted function.

[0074] (5-1) For both single-weighted functions and traversal-weighted functions, let's uniformly define the weighted function as w(f). The information lower limit threshold T of the weighted function is calculated using the following formula:

[0075]

[0076] In the formula, σ(w(f)) represents the average value of the weighted function w(f), σ(w(f)) represents the standard deviation of the weighted function w(f), and ε represents the scaling factor, which is recommended to be in the range of 3-5. In this invention, it is set to 3.

[0077] (5-2) Perform threshold filtering to reconstruct the weighting function w(f), and the calculation formula is as follows:

[0078]

[0079] S06, perform a priori ergodic weighting on the spectral coherence function to obtain the a priori ergodic spectral coherence, and further perform absolute value integration along the carrier frequency direction to obtain the a priori ergodic envelope spectrum.

[0080] (6-1) For the spectral coherence function γ x (α,f) is weighted to obtain the a priori ergodic spectral coherence function. The calculation formula is as follows:

[0081]

[0082] (6-2) The prior traversal weighting function γ x The absolute value integral along the carrier frequency direction (α,f) is used to obtain the a priori ergodic envelope spectrum PTES(α), and the calculation formula is as follows:

[0083]

[0084] In the formula, F S The sampling frequency for the monitoring signal.

[0085] To verify the effectiveness of this invention, simulation signals containing Gaussian noise, impulse noise, and cyclic stationary noise were analyzed, as follows:

[0086] The above methods are used to analyze the simulated signal x(t) containing Gaussian noise and impulse noise:

[0087]

[0088] In the formula, n p (t) and n s (t) represent impulse noise and Gaussian white noise, respectively. The parameter settings for the simulation signals are shown in Table 1 below:

[0089] Table 1

[0090]

[0091] Figure 2The time-domain signal of the simulation signal is shown in the figure, where the modulated signal is completely submerged by background noise. Figure 3 For narrowband demodulation of the kurtosis spectrum of the simulated signal, it is easy to know the fundamental modulation frequency α. B Its harmonics are almost completely drowned out by background noise (see Figure 3 If the kurtosis spectrum narrowband demodulation results are not ideal for extracting the frequency features of the signal (as shown in the middle circle), then the extraction of signal frequency features is not ideal.

[0092] Figure 4 The demodulation results of the cyclostationary analysis of the simulated signal, although its demodulation effect is better than that of kurtosis spectrum narrowband demodulation, the fundamental modulation frequency α B Its harmonics are quite obvious in the figure (see Figure 4 (Middle circle), but there are still a lot of interfering spectral lines.

[0093] Figure 5 To obtain the prior ergodic envelope spectrum of the simulated signal, the fundamental modulation frequency α is known. B Its harmonics are clearer than the previous two methods (see Figure 5 (Middle circle), and the noise component was suppressed to a certain extent.

[0094] The above results demonstrate that, under the combined interference of complex and intense noises such as Gaussian noise, impulse noise, and cyclostationary noise, the prior ergodic envelope spectrum proposed in this invention can accurately and effectively extract the modulation frequency components of the simulated signal.

[0095] The above method was used to process the vibration signal of the centrifugal pump with mixed electromagnetic interference, and the demodulation result was obtained. Figure 6-8 These are the kurtosis spectrum demodulation results of the centrifugal pump vibration signal, the cyclic stationary analysis demodulation results, and the prior ergodic envelope spectrum results proposed in this invention.

[0096] Depend on Figure 6 It can be seen that although the narrowband envelope demodulation method can indeed extract the target frequency (see...) Figure 6 (Middle circle), but also a large amount of electromagnetic noise exists; the results of cyclostationary analysis demodulation are slightly better than narrowband envelope demodulation (see Figure 7 (Middle circle), however, electromagnetic interference components still exist; finally, in the prior ergodic envelope spectrum results proposed in this invention, due to the assistance of prior knowledge of the axis frequency, only the axis frequency and its harmonics are significantly enhanced (see Figure 8 The inner circle represents the region where all other interferences are eliminated. These results demonstrate that even under complex and intense noise interference such as electromagnetic interference, the prior ergodic envelope spectrum proposed in this invention can still accurately and effectively extract the characteristic frequencies of the centrifugal pump vibration signal.

[0097] The embodiments described above provide a detailed explanation of the technical solutions and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the present invention. Any modifications, additions, and equivalent substitutions made within the scope of the principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for extracting rotational mechanical modulation features based on prior ergodic envelope spectrum, characterized in that, Includes the following steps: (1) Collect vibration or noise data of rotating machinery as monitoring signals, and calculate the spectral coherence function of the monitoring signals based on Gabor transform; (2) The spectral coherence function is sliced ​​along the carrier frequency direction to capture the harmonic cluster structure of different fundamental frequencies and evaluate the overall strength of the harmonic cluster structure to obtain the harmonic intensity vector; (3) Integrate the harmonic intensity vectors corresponding to all carrier frequencies to obtain the harmonic intensity matrix; (4) When the prior information is a single frequency, a single slice operation is performed on the harmonic intensity matrix to calculate a single weighting function; the specific process is as follows: (4-1) When the prior information is a single frequency, a single weighting function Directly from the harmonic intensity matrix At the cycle frequency carrier frequency The following slices were obtained: (4-2) For a single weighted function Perform normalization operation and calculate to obtain : When the prior information is within a certain frequency range, the harmonic intensity matrix is ​​traversed and sliced ​​to calculate the composite weighting function; the specific process is as follows: (4-a) When the prior information is a frequency range, assume that the frequency range of the prior information features is... Calculate the mean of all slices within this range. : In the formula, , The total number of frequencies included within the range; (4-b) After normalization, the composite weighted function is obtained. : In the formula, and They are respectively The maximum and minimum values; (5) Calculate the lower information threshold of the weighted function, perform threshold filtering operation, and obtain the prior traversal weighted function; (6) Perform a priori ergodic weighting on the spectral coherence function to obtain the a priori ergodic spectral coherence, and further perform absolute value integration along the carrier frequency direction to obtain the a priori ergodic envelope spectrum.

2. The method for extracting rotating mechanical modulation features based on prior ergodic envelope spectrum according to claim 1, characterized in that, In step (1), the specific process of calculating the spectral coherence function of the monitoring signal is as follows: (1-1) Collect vibration or noise data of rotating machinery as monitoring signals , Refers to the sampling frequency The moment of gain, among which The monitoring signal sampling time is The unit is seconds (s). (1-2) Calculate the monitoring signal Short-time Fourier Transform : In the formula, For window width, For the movement step size, For window functions, for abbreviation, For discrete frequencies, Frequency resolution ; (1-3) Monitoring signals Short-time Fourier Transform Phase correction is performed to obtain the Gabor transform result. : In the formula, It is a signal exist At that moment, with Centered on, bandwidth is The complex envelope, Indicates the energy flow within the frequency band; (1-4) Calculate the monitoring signal The average cyclic period spectrum is related to: ： In the formula, , for in the long In the signal, the step size The total number of moved windows, and the length of these windows is... , The cycle frequency, Where is the carrier frequency, ; (1-5) Calculate the monitoring signal spectral correlation function The calculation formula is as follows: (1-6) Calculate the monitoring signal spectral coherence function The calculation formula is as follows: 。 3. The method for extracting rotating mechanical modulation features based on prior ergodic envelope spectrum according to claim 1, characterized in that, In step (2), the specific process of obtaining the harmonic intensity vector is as follows: (2-1) Spectral correlation function The discrete form is Its cycle frequency Include Each possible value Spectral frequency Include Each possible value ; (2-2) Slice the spectral coherence function along the carrier frequency direction to obtain the spectral coherence function at the spectral frequency. slices at the location ; (2-3) Select the cycle frequency The fundamental frequency of the harmonic structure is used to find the peak value within a range of the fundamental frequency amplitude. The calculation formula is as follows: In the formula, For the base frequency The peak-seeking range To find the peak within a single-sided range, For the frequency in the spectrum For the fundamental frequency Peak amplitude; (2-4) Determine the harmonic order of the harmonic cluster structure For the fundamental frequency of interest The peak finding of all harmonics is performed within a range, and the calculation formula is as follows: In the formula, Harmonics The peak-seeking range To find the peak within a single-sided range, For the frequency in the spectrum Regarding harmonics Peak amplitude, ; (2-5) Regarding the fundamental frequency and all its harmonics The peak amplitude is taken as the arithmetic mean, and the spectral frequency is calculated. Cycle frequency Harmonic intensity at: In the formula, Spectral frequency Cycle frequency The harmonic intensity at that location.

4. The method for extracting rotating mechanical modulation features based on prior ergodic envelope spectrum according to claim 1, characterized in that, In step (3), the specific process of obtaining the harmonic intensity matrix is ​​as follows: Assuming cycle frequency exist Within the range, while spectral frequency exist Within the range, If the upper bound of the target cycle frequency range is given, then and When increasing in different directions, the harmonic intensity vector will... Integrate into harmonic intensity matrix .

5. The method for extracting rotating mechanical modulation features based on prior ergodic envelope spectrum according to claim 1, characterized in that, In step (5), the specific process of calculating the prior traversal weighted function is as follows: (5-1) For both single-weighted functions and composite-weighted functions, we uniformly define them as weighted functions. Calculate the information lower limit threshold of the weighted function. The calculation formula is as follows: In the formula, Represents the weighting function The average value, Represents the weighting function standard deviation This represents the scaling factor, with a value range of 3-5; (5-2) Perform threshold filtering operation and reconstruct the weighting function. The calculation formula is as follows: 。 6. The method for extracting rotating mechanical modulation features based on prior ergodic envelope spectrum according to claim 5, characterized in that, In step (6), the specific process of calculating the prior ergodic envelope spectrum is as follows: (6-1) For the spectral coherence function After weighting, the prior ergodic spectrum coherence function is obtained. The calculation formula is as follows: (6-2) The prior ergodic spectral coherence function The a priori ergodic envelope spectrum is obtained by calculating the absolute value integral along the carrier frequency direction. The calculation formula is as follows: In the formula, The sampling frequency for the monitoring signal.

Images