A rotating machinery modulation feature extraction method based on adaptive weighted envelope spectrum

Abstract

The application discloses a rotating machinery modulation feature extraction method based on an adaptive weighted envelope spectrum, comprising the following steps: (1) collecting noise data of the rotating machinery as a monitoring signal, and calculating a spectral correlation function of the monitoring signal; (2) calculating a power spectral density of the monitoring signal as a reference value, normalizing the spectral correlation function, and obtaining a spectral coherence function; (3) slicing the spectral coherence function along a cycle frequency direction to obtain spectral coherence slice results at different carrier frequencies, and performing Fourier transform to obtain a modulation intensity distribution function; (4) calculating kurtosis coefficients of the modulation intensity distribution function at different carrier frequencies, and obtaining an adaptive weighting function; (5) calculating an adaptive weighted spectral coherence function by using the adaptive weighting function; and (6) integrating the adaptive weighted spectral coherence function along the carrier frequency direction to obtain an adaptive weighted envelope spectrum. The application can effectively extract modulation features in a rotating machinery vibration noise signal under complex and 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 machinery modulation features based on adaptive weighted envelope spectrum. Background Technology

[0002] Rotating machinery is widely used in industrial production, water conservancy and hydropower, military applications and other fields. The extraction of its modulation features not only helps to realize the condition monitoring and fault diagnosis of rotating machinery such as pumps, fans and compressors in the industrial field, but also helps to realize the target detection and target identification of non-cooperative targets such as helicopters, ships and submarines in the military field.

[0003] Signal demodulation analysis is the mainstream approach for extracting modulation features in rotating machinery, and mainly includes signal processing methods such as narrowband envelope demodulation and cyclostationary analysis.

[0004] For example, Chinese patent document CN110569812A discloses an envelope demodulation method and system for fault signals. By utilizing the upper and lower envelopes of the signal, the demodulation of the modulated signal is achieved through simple algebraic operations. The modulated fault signal can be directly demodulated from the original signal without being disturbed by the low-frequency components of the original signal.

[0005] Chinese patent document CN110763464A discloses a method for extracting fault features of rolling bearings based on cyclic stationary analysis, comprising: acquiring a vibration signal X(t) of a gearbox, performing spectral kurtosis analysis to obtain the spectral kurtosis value of the vibration signal X(t), and determining the center frequency and bandwidth of the frequency band corresponding to the maximum spectral kurtosis value; determining the frequency band range based on the center frequency and bandwidth, and performing bandpass filtering on the vibration signal X(t) within the frequency band range to obtain a filtered signal Y(t); performing second-order cyclic stationary analysis on the filtered signal Y(t) to obtain a spectral correlation density function; and using the spectral correlation density function to obtain the fault features of the outer ring, inner ring, rolling elements, and cage of the rolling bearing.

[0006] Narrowband envelope demodulation mainly includes three steps: signal frequency band decomposition, demodulation frequency band selection, and envelope signal demodulation. First, the signal is decomposed into a series of narrowband filtered signals according to specific rules. Second, the information content of each narrowband filtered signal is evaluated using metrics such as kurtosis, negative entropy, and cyclostationarity, and the frequency band with the highest information content is selected as the optimal demodulation frequency band. Finally, the time-domain signal of the optimal demodulation frequency band is obtained, and envelope demodulation is performed to obtain the envelope spectrum. This method is suitable for local defect detection, weak fault diagnosis, and unknown target detection in rotating machinery. Its disadvantage is that narrowband envelope demodulation usually only selects a specific narrowband for analysis, its algorithm has poor noise resistance, and it is prone to failure under low signal-to-noise ratio conditions containing multiple types of noise, thus limiting its effectiveness.

[0007] Cyclostationary analysis can effectively reduce the interference of stationary noise in the monitoring signal by utilizing the cyclostationary characteristics of the modulated signal in rotating machinery. This method often uses spectral correlation function and spectral coherence function to characterize the distribution of the modulated signal in the monitoring signal, and describes the changes in modulation intensity with cyclic frequency and carrier frequency, thus extracting the modulation frequency components relatively well. However, this method is also limited in its application to extract the characteristic frequencies of rotating machinery. In addition to the dominant modulation frequency caused by the rotating machinery, the monitoring signal usually contains cyclostationary noise interference from other mechanical, electrical, and communication equipment. Traditional cyclostationary analysis cannot reduce or remove this cyclostationary noise interference, which will affect the identification and extraction of the modulation features of rotating machinery. Summary of the Invention

[0008] This invention provides a method for extracting modulation features of rotating machinery based on adaptive weighted envelope spectrum. It can effectively extract modulation features from vibration noise signals of rotating machinery under complex and strong noise interference, and can be applied to fields such as condition monitoring, fault diagnosis, target detection, and target recognition of rotating machinery.

[0009] A method for extracting rotational machinery modulation features based on adaptive weighted envelope spectrum, characterized by comprising the following steps:

[0010] (1) Collect vibration or radiated noise data of rotating machinery as monitoring signals for subsequent analysis, and calculate the spectral correlation function of the monitoring signals;

[0011] (2) Calculate the power spectral density of the monitoring signal, use it as a reference value to normalize the spectral correlation function, and obtain the spectral coherence function;

[0012] (3) Slice the spectral coherence function along the direction of the cyclic frequency to obtain the spectral coherence slice results at different carrier frequencies. Perform Fourier transform on the slice results at each carrier frequency to obtain the modulation intensity distribution function.

[0013] (4) Calculate the kurtosis coefficient of the modulation intensity distribution function at different carrier frequencies, and square the kurtosis coefficient at each carrier frequency to obtain the adaptive weighting function;

[0014] (5) The spectral coherence function is weighted using an adaptive weighting function to calculate the adaptive weighted spectral coherence function;

[0015] (6) Integrate the adaptive weighted spectrum coherence function along the carrier frequency direction to obtain the adaptive weighted envelope spectrum.

[0016] In step (1), the noise data is either an airborne acoustic signal or an underwater acoustic signal.

[0017] The specific process for calculating the spectral correlation function of the monitoring signal is as follows:

[0018] (1-1) Calculate the instantaneous autocorrelation function R of the monitoring signal x(t). x (t,τ), the calculation formula is as follows:

[0019] R x (t,τ)=E[x(t)x * (t-τ)]

[0020] In the formula, t represents time, τ represents time delay, and E[·] represents the ensemble averaging operator. * The conjugate symbol;

[0021] (1-2) For the instantaneous autocorrelation function R x Performing a Fourier transform on (t,τ) to convert time t to the cyclic frequency α yields the cyclic autocorrelation function R. x (α,τ), the calculation formula is as follows:

[0022]

[0023] In the formula, -T and T are the lower and upper time limits, e is the natural constant, j is the imaginary number, π is pi, and α is the cycle frequency;

[0024] (1-3) For the cyclic autocorrelation function R x Performing a Fourier transform on (α,τ) converts the time delay τ to the carrier frequency f, yielding the spectral correlation function S. x (α,f), the calculation formula is as follows:

[0025]

[0026] In step (2), the specific process for calculating the spectral coherence function is as follows:

[0027] (2-1) Calculate the time-averaged autocorrelation function R of the monitoring signal x (τ):

[0028]

[0029] In the formula, t represents time, and τ represents time delay. * The conjugate symbol;

[0030] (2-2) For the time-averaged autocorrelation function R x Perform a Fourier transform on (τ) to convert the time delay τ to the carrier frequency f, and obtain the power spectral density function P. x (f), the calculation formula is as follows:

[0031]

[0032] (2-3) The power spectral density function P x (f) Using this as a reference value, the spectral correlation function is normalized to obtain the spectral coherence function γ. x (α,f), the calculation formula is as follows:

[0033]

[0034] In the formula, S x (α,f) is the spectral correlation function.

[0035] In step (3), the specific process for calculating the modulation intensity distribution function is as follows:

[0036] (3-1) The spectral coherence function γ x (α,f) Slicing along the cyclic frequency direction yields spectral coherence slices at different carrier frequencies: γ x (α,f1),γ x (α,f2),...,γ x (α,f N ), where f1 represents the first carrier frequency, f2 represents the second carrier frequency, and so on, f... N This represents the Nth carrier frequency, which is the last carrier frequency.

[0037] (3-2) Perform Fourier transform on the spectral coherence slices at different carrier frequencies to obtain the modulation intensity distribution function Γ. x (η,f):

[0038]

[0039] In the formula, η is the double cycle time, and the carrier frequency f takes values ​​in the range {f1, f2, ..., f...}. N}

[0040] In step (4), the specific process of calculating the adaptive weighting function is as follows:

[0041] (4-1) Calculate the modulation intensity distribution function Γ x The kurtosis coefficient k(f) at different carrier frequencies (η,f):

[0042]

[0043] In the formula, E[·] represents the ensemble average operator, μ x (η,f) represents the frequency domain mean of the modulation intensity distribution function, σ x (η,f) represents the frequency domain standard deviation of the modulation intensity distribution function;

[0044] (4-2) By squaring the kurtosis coefficient k(f), we obtain the adaptive weighting function w(f):

[0045] w(f) = k 2 (f).

[0046] In step (5), the formula for calculating the adaptive weighted spectral coherence function is:

[0047]

[0048] In the formula, Let w(f) represent the adaptive weighted spectral coherence function, and γ represent the adaptive weighting function. x (α,f) represents the spectral coherence function.

[0049] In step (6), the formula for calculating the adaptive weighted envelope spectrum is:

[0050]

[0051] In the formula, For adaptive weighted envelope spectrum; F s This represents the sampling frequency of the monitoring signal.

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

[0053] 1. This invention proposes a method for calculating the modulation intensity distribution function using the spectral coherence function, which can reflect the strength of modulation characteristics at different carrier frequencies.

[0054] 2. This invention proposes a method for calculating an adaptive weighting function. This function adaptively takes a larger value at carrier frequencies with higher modulation intensity and a smaller value at carrier frequencies with lower modulation intensity, thereby achieving signal enhancement of modulation characteristics.

[0055] 3. The method proposed in this invention is based on cyclostationary analysis and adaptive weighted spectral coherence processing for signal demodulation. It can accurately and effectively extract the modulation frequency components of rotating machinery without the need for prior information on the modulation frequency. The relevant information can be used for condition monitoring, fault diagnosis, target detection, and target recognition of rotating machinery. Attached Figure Description

[0056] Figure 1 This is a flowchart illustrating a method for extracting rotating machinery modulation features based on adaptive weighted envelope spectrum according to the present invention.

[0057] Figure 2 The results of enhanced envelope spectrum and adaptive weighted envelope spectrum analysis of the three-bladed propeller in this embodiment of the invention are shown.

[0058] Figure 3The results of enhanced envelope spectrum and adaptive weighted envelope spectrum analysis of the four-bladed propeller in this embodiment of the invention are shown.

[0059] Figure 4 The results of enhanced envelope spectrum and adaptive weighted envelope spectrum analysis of the five-bladed propeller in this embodiment of the invention are shown. Detailed Implementation

[0060] 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.

[0061] like Figure 1 As shown, a method for extracting rotational machinery modulation features based on adaptive weighted envelope spectrum includes the following steps:

[0062] A method for extracting features of rotating machinery modulation based on adaptive weighted envelope spectrum, characterized by comprising the following steps;

[0063] S01, collect vibration or radiated noise data (airborne acoustic signal or underwater acoustic signal) of rotating machinery as monitoring signals for subsequent analysis, and calculate the spectral correlation function of the monitoring signal.

[0064] The specific process for calculating the spectral correlation function of the monitoring signal is as follows:

[0065] (1-1) Calculate the instantaneous autocorrelation function R of the monitoring signal x(t). x (t,τ), the calculation formula is as follows:

[0066] R x (t,τ)=E[x(t)x * (t-τ)]

[0067] In the formula, t represents time, τ represents time delay, and E[·] represents the ensemble averaging operator. * The conjugate symbol;

[0068] (1-2) For the instantaneous autocorrelation function R x Performing a Fourier transform on (t,τ) to convert time t to the cyclic frequency α yields the cyclic autocorrelation function R. x (α,τ), the calculation formula is as follows:

[0069]

[0070] In the formula, -T and T are the lower and upper time limits, e is the natural constant, j is the imaginary number, π is pi, and α is the cycle frequency.

[0071] (1-3) Perform a Fourier transform on the cyclic autocorrelation function to convert the time delay τ to the carrier frequency f, and obtain the spectral correlation function S.x (α,f), the calculation formula is as follows:

[0072]

[0073] S02, calculate the power spectral density of the monitoring signal, use it as a reference value to normalize the spectral correlation function, and obtain the spectral coherence function.

[0074] The specific process for calculating the spectral coherence function of the monitoring signal is as follows:

[0075] (2-1) Calculate the time-averaged autocorrelation function R of the monitoring signal x (τ):

[0076]

[0077] In the formula, t represents time, and τ represents time delay. * This is the conjugate symbol.

[0078] (2-2) For the time-averaged autocorrelation function R x Perform a Fourier transform on (τ) to convert the time delay τ to the carrier frequency f, and obtain the power spectral density function P. x (f), the calculation formula is as follows:

[0079]

[0080] (2-3) The power spectral density function P x (f) Using this as a reference value, the spectral correlation function is normalized to obtain the spectral coherence function γ. x (α,f), the calculation formula is as follows:

[0081]

[0082] S03, slice the spectral coherence function along the cyclic frequency direction to obtain the spectral coherence slice results at different carrier frequencies, and perform Fourier transform on the slice results at each carrier frequency to obtain the modulation intensity distribution function.

[0083] The specific process for calculating the modulation intensity distribution function is as follows:

[0084] (3-1) Slice the spectral coherence function along the cyclic frequency direction to obtain the spectral coherence slice results at different carrier frequencies: γ x (α,f1),γ x (α,f2),...,γ x (α,f N ), where f1 represents the first carrier frequency, f2 represents the second carrier frequency, and so on, f... N This represents the Nth carrier frequency, which is the last carrier frequency.

[0085] (3-2) Perform Fourier transform on the spectral coherence slices at different carrier frequencies to obtain the modulation intensity distribution function Γ. x (η,f):

[0086]

[0087] In the formula, η is the double cycle time, and the carrier frequency f takes values ​​in the range {f1, f2, ..., f...}. N}

[0088] S04, calculate the kurtosis coefficient of the modulation intensity distribution function at different carrier frequencies, and square the kurtosis coefficient at each carrier frequency to obtain the adaptive weighting function.

[0089] The specific process for calculating the adaptive weighting function is as follows:

[0090] (4-1) Calculate the kurtosis coefficient k(f) of the modulation intensity distribution function at different carrier frequencies:

[0091]

[0092] In the formula, E[·] represents the ensemble average operator, μ x (η,f) represents the frequency domain mean of the modulation intensity distribution function, σ x (η,f) represents the frequency domain standard deviation of the modulation intensity distribution function;

[0093] (4-2) By squaring the kurtosis coefficient, we obtain the adaptive weighting function w(f):

[0094] w(f) = k 2 (f)

[0095] S05, the spectral coherence function is weighted using an adaptive weighting function to calculate the adaptive weighted spectral coherence function.

[0096] The specific process of obtaining the adaptive weighted spectral coherence function is as follows:

[0097]

[0098] In the formula, Let w(f) represent the adaptive weighted spectral coherence function, and γ represent the adaptive weighting function. x (α,f) represents the spectral coherence function.

[0099] S06. Integrate the adaptive weighted spectral coherence function along the carrier frequency direction to obtain the adaptive weighted envelope spectrum.

[0100] The specific process of obtaining the adaptive weighted envelope spectrum is as follows:

[0101]

[0102] In the formula, This is an adaptive weighted envelope spectrum.

[0103] To verify the effectiveness of the present invention, the noise of propellers with three different numbers of blades was analyzed.

[0104] The three-bladed propeller rotates at approximately 111 revolutions per minute. The enhanced envelope spectrum (EES) used in traditional cyclic stationarity analysis methods and the adaptive weighted envelope spectrum (AWES) proposed in this invention are as follows: Figure 2 As shown in the figure, there are several pairs of vertical dashed lines. Each pair of vertical dashed lines represents the search range of the frequency conversion and its harmonics. The point with the highest amplitude within this range is the frequency conversion and its harmonics obtained by the demodulation method. The dotted line represents the signal noise threshold. It can be seen that, compared to EES, in the first three harmonics, the line spectrum amplitude of AWES is significantly different from the dashed line representing the noise level.

[0105] The four-bladed propeller rotates at approximately 83 revolutions per minute. Its enhanced envelope spectrum (EES) and adaptive weighted envelope spectrum (AWES) are as follows: Figure 3 As shown in the figure, there are several pairs of vertical dashed lines. Each pair of vertical dashed lines represents the search range of the frequency shift and its harmonics. The point with the highest amplitude within this range is the frequency shift and its harmonics obtained by the demodulation method. The dotted line represents the signal noise threshold. Similarly, compared to EES, especially in the 4th harmonic, the line spectrum amplitude of AWES is very prominent relative to the noise level.

[0106] The five-bladed propeller rotates at approximately 72 revolutions per minute. Its enhanced envelope spectrum (EES) and adaptive weighted envelope spectrum (AWES) are as follows: Figure 4 As shown in the figure, there are several pairs of vertical dashed lines. Each pair of vertical dashed lines represents the search range of the frequency shift and its harmonics. The point with the highest amplitude within this range is the frequency shift and its harmonics obtained by the demodulation method. The dotted lines represent the signal noise threshold. In EES, only the 1st and 5th harmonics have relatively weak line spectra, while in AWES, the line spectrum amplitude is very obvious relative to the noise level at several harmonics.

[0107] pass Figures 2 to 4 The comparison leads to the following conclusions:

[0108] 1. AWES can extract the modulation features corresponding to shaft frequency and its harmonics, especially blade frequency, well under different blade numbers.

[0109] 2. Compared with traditional EES, AWES extracts line spectra with significant advantages in modulation feature extraction for all blade numbers.

[0110] Compared to traditional enhanced envelope spectra, the adaptive weighted envelope spectra proposed in this invention can better extract the line spectra corresponding to modulation features. The adaptive weighting process, which does not require any empirical parameter settings, can further enhance the efficiency of modulation feature extraction without human interference. At the same time, based on the propeller underwater acoustic signal processing results, AWES can effectively extract the modulation features of propellers with three, four, and five blades.

[0111] 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 machinery modulation features based on adaptive weighted envelope spectrum, characterized in that, Includes the following steps: (1) Collect vibration or radiated noise data of rotating machinery as monitoring signals for subsequent analysis, and calculate the spectral correlation function of the monitoring signals; (2) Calculate the power spectral density of the monitoring signal, use it as a reference value to normalize the spectral correlation function, and obtain the spectral coherence function; (3) Slice the spectral coherence function along the cyclic frequency direction to obtain spectral coherence slice results at different carrier frequencies. Perform Fourier transform on the slice results at each carrier frequency to obtain the modulation intensity distribution function. The specific process for calculating the modulation intensity distribution function is as follows: (3-1) The spectral coherence function γ x (α,f) Slicing along the cyclic frequency direction yields spectral coherence slices at different carrier frequencies: γ x (α,f1),γ x (α,f2),...,γ x (α,f N ), where f1 represents the first carrier frequency, f2 represents the second carrier frequency, and so on, f... N This represents the Nth carrier frequency, which is the last carrier frequency. (3-2) Perform Fourier transform on the spectral coherence slices at different carrier frequencies to obtain the modulation intensity distribution function Γ. x (η,f): In the formula, η is the double cycle time, and the carrier frequency f takes values ​​in the range {f1, f2, ..., f...}. N }; (4) Calculate the kurtosis coefficients of the modulation intensity distribution function at different carrier frequencies, and square the kurtosis coefficients at each carrier frequency to obtain the adaptive weighting function; the specific process for calculating the adaptive weighting function is as follows: (4-1) Calculate the modulation intensity distribution function Γ x The kurtosis coefficient k(f) at different carrier frequencies (η,f): In the formula, Ε[·] represents the ensemble average operator, μ x (η,f) represents the frequency domain mean of the modulation intensity distribution function, σ x (η,f) represents the frequency domain standard deviation of the modulation intensity distribution function; (4-2) By squaring the kurtosis coefficient k(f), we obtain the adaptive weighting function w(f): w(f)=k 2 (f) (5) The spectral coherence function is weighted using an adaptive weighting function to calculate the adaptive weighted spectral coherence function; (6) Integrate the adaptive weighted spectrum coherence function along the carrier frequency direction to obtain the adaptive weighted envelope spectrum.

2. The method for extracting rotating machinery modulation features based on adaptive weighted envelope spectrum according to claim 1, characterized in that, In step (1), the noise data is either an airborne acoustic signal or an underwater acoustic signal.

3. The method for extracting rotational machinery modulation features based on adaptive weighted envelope spectrum according to claim 1, characterized in that, In step (1), the specific process of calculating the spectral correlation function of the monitoring signal is as follows: (1-1) Calculate the instantaneous autocorrelation function R of the monitoring signal x(t). x (t,τ), the calculation formula is as follows: R x (t,τ)=E[x(t)x * (t-t)] In the formula, t represents time, τ represents time delay, Ε[·] represents the set averaging operator, and * is the conjugate symbol; (1-2) For the instantaneous autocorrelation function R x Performing a Fourier transform on (t,τ) to convert time t to the cyclic frequency α yields the cyclic autocorrelation function R. x (α,τ), the calculation formula is as follows: In the formula, -T and T are the lower and upper time limits, e is the natural constant, j is the imaginary number, π is pi, and α is the cycle frequency; (1-3) For the cyclic autocorrelation function R x Performing a Fourier transform on (α,τ) converts the time delay τ to the carrier frequency f, yielding the spectral correlation function S. x (α,f), the calculation formula is as follows:

4. The method for extracting rotating machinery modulation features based on adaptive weighted envelope spectrum according to claim 1, characterized in that, In step (2), the specific process for calculating the spectral coherence function is as follows: (2-1) Calculate the time-averaged autocorrelation function R of the monitoring signal x (τ): In the formula, t represents time, τ represents time delay, and * is the conjugate symbol; (2-2) For the time-averaged autocorrelation function R x Perform a Fourier transform on (τ) to convert the time delay τ to the carrier frequency f, and obtain the power spectral density function P. x (f), the calculation formula is as follows: (2-3) The power spectral density function P x (f) Using this as a reference value, the spectral correlation function is normalized to obtain the spectral coherence function γ. x (α,f), the calculation formula is as follows: In the formula, S x (α,f) is the spectral correlation function.

5. The method for extracting rotational machinery modulation features based on adaptive weighted envelope spectrum according to claim 1, characterized in that, In step (5), the formula for calculating the adaptive weighted spectral coherence function is: In the formula, Let w(f) represent the adaptive weighted spectral coherence function, and γ represent the adaptive weighting function. x (α,f) represents the spectral coherence function.

6. The method for extracting rotating machinery modulation features based on adaptive weighted envelope spectrum according to claim 5, characterized in that, In step (6), the formula for calculating the adaptive weighted envelope spectrum is: In the formula, S x AWES (α) represents the adaptive weighted envelope spectrum; F s This represents the sampling frequency of the monitoring signal.

Images