A Mechanical Fault Diagnosis Method Based on Maximizing the Generalized Cyclostationary Index

Abstract

This invention discloses a mechanical fault diagnosis method based on maximizing the generalized cyclostationary index, belonging to the fields of signal processing and mechanical equipment condition monitoring and fault diagnosis. It mainly involves blind deconvolution of the signal based on maximizing the generalized cyclostationary index to extract the cyclostationary features of interest and perform envelope analysis. The method includes shape parameter estimation based on a generalized Gaussian cyclostationary model; blind deconvolution based on maximizing generalized cyclostationarity; and β-order envelope spectrum analysis. The beneficial effects of this invention are that, based on maximizing the generalized cyclostationary index, it improves the robustness of cyclostationary feature extraction for non-Gaussian distributed signals; and by setting a tolerance bandwidth in the weighting matrix generation, it solves the problem of inconsistency between the theoretical and actual cyclostationary frequencies.

Description

Technical Field

[0001] This invention relates to the fields of signal processing and equipment condition monitoring and fault diagnosis, and particularly to a method for extracting fault features of rotating equipment based on vibration signals. Background Technology

[0002] Rolling bearings are one of the most important supporting components in the transmission chain of rotating equipment. Once they fail, they will affect the normal operation of the rotating equipment, and in severe cases, they may even lead to the shutdown of the entire equipment for maintenance, resulting in unnecessary economic losses.

[0003] Chinese patent document CN113536226A discloses a blind deconvolution algorithm for enhancing the characteristics of rotating machinery fault signals, including the following steps: S1: Construct multiple cascaded FIR filters, determine the maximization criterion of the blind deconvolution algorithm based on the characteristics of the original vibration signal, and use the maximization criterion as the objective function; S2: Use the cascaded FIR filters to perform convolution operations on the original vibration signal in sequence to obtain the filtered signal, and calculate the objective function value of the filtered signal; S3: Use a backward automatic differentiation algorithm to calculate the gradient of the objective function value with respect to the filter at the current iteration number; S4: Update the values ​​of all filters; S5: Repeat S2-S4 until the maximum number of iterations is reached, and output the final filtered signal. Chinese patent document CN114528525A discloses a mechanical fault diagnosis method based on maximum reweighted kurtosis blind deconvolution, belonging to the field of wind turbine fault diagnosis technology, and proposes a new blind deconvolution method, namely maximum reweighted kurtosis blind deconvolution. Reweighted kurtosis is highly robust to single or small numbers of strong impulses in fault signals and does not require prior knowledge of the fault impulse sequence to be recovered.

[0004] While current blind deconvolution methods and their variants based on maximizing the second-order cyclostationarity index are used to extract fault features, these methods all assume that the acquired signal follows a Gaussian distribution. However, due to the harsh operating conditions of transmission systems, signals often exhibit non-Gaussian distributions, thus reducing the reliability of fault feature extraction. Furthermore, blind deconvolution methods based on maximizing the second-order cyclostationarity index require the used cyclostationary frequency to match the actual fault frequency; however, due to factors such as the slippage of rolling bearings, there will be a certain deviation between the theoretical and actual fault frequencies, thus affecting fault feature extraction. Summary of the Invention

[0005] The fault feature extraction method for mechanical equipment of the present invention is based on processing the vibration signals of rotating or reciprocating mechanical equipment such as wind turbines. In the prior art, the vibration signal is assumed to have a generalized Gaussian distribution, therefore the fault features generated by the mechanical equipment exhibit generalized cyclostationary characteristics. However, based on the generalized Gaussian distribution assumption, fault features cannot be accurately predicted in non-Gaussian cases; faults do not always occur as expected. Therefore, the present invention improves the robustness of fault feature extraction in non-Gaussian cases by performing blind deconvolution on the signal based on maximizing the generalized cyclostationary component and introducing tolerance bandwidth in the weighted matrix generation to iteratively find the actual fault frequency of bearing faults.

[0006] The present invention provides a mechanical fault diagnosis method based on maximizing a generalized cyclic stationary index, comprising the following steps:

[0007] Step 1: Collect vibration signals and preprocess them to meet processing requirements;

[0008] Step 2: Set up a filter, input the preprocessed vibration signal into the filter, and perform iterative calculations;

[0009] Step 3: Compare the number of iterations and the convergence condition with the set maximum number of blind deconvolution iterations Nmax1 and convergence condition ε1 to decide whether to stop iterative updates;

[0010] Step 4: Under the condition of stopping iterative updates, obtain the filtered signal s, generate the β-order envelope spectrum of the filtered signal s for analysis, and save the corresponding index values ​​and shape parameters.

[0011] Furthermore, in step 2, the parameters used by the filter are set; specifically, the filter length L; the maximum number of iterations Nmax1 and Nmax2 for blind deconvolution and shape parameter β estimation; the convergence conditions ε1 and ε2 for blind deconvolution and shape parameter β estimation; the theoretical cycle frequency α0 and the tolerance bandwidth w;

[0012] After setting the parameters used by the filter, initialize the filter h.

[0013] Furthermore, step 2.1 of the iterative calculation in step 2 involves using a loop order α. i Calculate the synchronous average r of the envelope signal |x|, and then normalize the signal to obtain x. t ,

[0014] in:

[0015]

[0016] In equation (1), <r>Let r be the first moment of r.

[0017]

[0018] Where M is the number of sampling points in the cycle order α, and its relationship with the cycle order and sampling order fs is M = fs / α, therefore we can obtain

[0019]

[0020] Furthermore, step 2.2 of the iterative calculation in step 2 is to calculate the shape parameters of signal x under the generalized Gaussian cyclic stationary model using the iterative method, where:

[0021] (1) Calculate |x t The first moment of |x, μ=<|x t |>;

[0022] (2) Initialize the shape parameters, specifically:

[0023]

[0024] (3) Iteratively update the shape parameters, specifically:

[0025]

[0026] in,

[0027]

[0028]

[0029] Where ψ and ψ′ are the digamma and trigamma functions, respectively.

[0030] Furthermore, in step 4, the filtered signal s is obtained by convolving the signal x with the filter h;

[0031] s=x*h (8)

[0032] Furthermore, in step 4, the β-order envelope spectrum (PES) of the filtered signal s is generated, where

[0033]

[0034] The filtered signal s is obtained by convolving the signal x with the filter h.

[0035] Furthermore, the β-order envelope spectrum (PES) described in claim 6, within its tolerance bandwidth α i -w:α i Find local maxima within the range of +w and update the cycle frequency α. i+1 :

[0036]

[0037] in k = 1, ..., M is the tolerance bandwidth of the k-th harmonic of the fundamental frequency of the fault.

[0038] Furthermore, the weighted correlation matrix W of the signal is generated by the β-order envelope spectrum (PES) as described in claim 6, wherein:

[0039]

[0040]

[0041] because The calculation includes generating the β-order envelope spectrum. Therefore, the actual fault frequency can be determined by setting the bandwidth.

[0042] Furthermore, after generating the weighted correlation matrix W of the signal, the coefficients of the filter are updated using the Rayleigh quotient. The Rayleigh quotient expression for the generalized cyclostationarity index is as follows:

[0043]

[0044] in,

[0045] R XX =x H x (14)

[0046] R XWX =x H Wx (15)

[0047] Where, x H Let x be the Hermitian matrix.

[0048] Furthermore, the number of iterations and the convergence condition are compared with the set maximum number of blind deconvolution iterations Nmax1 and convergence condition ε1 to determine whether to update the filter again through step 2 or to stop the iteration update.

[0049] The beneficial effects of this invention are as follows: This invention provides a blind deconvolution method based on maximizing the generalized cyclostationary index. It assumes the signal to be a generalized Gaussian distribution and performs blind deconvolution based on maximizing the generalized cyclostationary index by calculating the shape parameters of the generalized Gaussian cyclostationary model. A tolerance bandwidth is set in the weighted matrix generation to iteratively find the actual fault frequency. Unlike the blind deconvolution method that maximizes the second-order cyclostationary index, which assumes the signal to be Gaussian, this method, by estimating the shape parameters of the signal's generalized Gaussian cyclostationary model, can effectively solve the problem of difficulty in extracting fault features when the signal exhibits a non-Gaussian distribution due to factors such as impacts and electromagnetic interference. Furthermore, the fault frequency can be iteratively updated, thus solving the problem of deviation between the cyclostationary frequency used in the input parameters and the actual fault frequency. Attached Figure Description

[0050] Figure 1 This is a schematic diagram of the method steps of the present invention;

[0051] Figure 2 The average high-speed shaft rotation speed and time-domain waveform of vibration signal when collecting vibration signals for wind turbine units in a wind farm;

[0052] Figure 3 This indicates the trend value and shape parameters of the filtered signal.

[0053] Figure 4 This represents a waterfall plot of the beta-order envelope spectrum of the signal. Detailed Implementation

[0054] The preferred embodiments of the present invention will now be described in detail with reference to the accompanying drawings, so that the advantages and features of the present invention can be more easily understood by those skilled in the art, thereby providing a clearer and more explicit definition of the scope of protection of the present invention.

[0055] The data for this example comes from the high-speed shaft radial vibration data of a 2.0MW wind turbine generator in an actual wind farm, spanning from December 2016 to June 2017. The figure shows the time-domain waveform of the high-speed shaft vibration signal and the average rotational speed of the high-speed shaft. The following will detail the fault identification and condition index quantification process for the bearing pair, including the following steps:

[0056] Step 1: Track the vibration signal by average rotational speed, converting the vibration signal from the time domain to the angle domain and the frequency domain to the order domain.

[0057] Step 2: Calculate the theoretical fault order α (9.27) according to the fault order calculation formula of each bearing component, and set the parameters used in the method, including filter length L (20); maximum number of iterations for blind deconvolution and shape parameter β estimation Nmax1 (50), Nmax2 (50); convergence conditions for blind deconvolution and shape parameter β estimation ε1 (0.001), ε2 (0.001); theoretical fault order α0 (9.27*[1×20]) and tolerance bandwidth w (0.2).

[0058] Step 3: Initialize the filter h. The filter initialization is generated by the autoregressive model method.

[0059] Step 4: Through the cyclic order α i Calculate the synchronous average r of the envelope signal |x|, and then normalize the signal to obtain x. t ,in:

[0060]

[0061] In equation (1), <r>Let r be the first moment of r.

[0062]

[0063] Where M is the number of sampling points in the cycle order α, and its relationship with the cycle order and sampling order fs is M = fs / α, therefore we can obtain

[0064]

[0065] Step 5: Calculate the shape parameters of signal x under the generalized Gaussian cyclic stationary model using the Newton-Raphson iterative method, where:

[0066] (1) Calculate |x t The first moment of |x, μ=<|x t |>;

[0067] (2) Initialize the shape parameters, specifically:

[0068]

[0069] (3) Iteratively update the shape parameters, specifically:

[0070]

[0071] in,

[0072]

[0073]

[0074] Where ψ and ψ′ are the digamma and trigamma functions, respectively;

[0075] (4) Stop iterating by estimating the maximum number of iterations Nmax2 and the convergence condition ε2 using the set shape parameter β. Otherwise, continue to repeat step 5.

[0076] Step 6: Obtain the filtered signal s by convolving the signal x with the filter h;

[0077] s=x*h (8)

[0078] Step 7: Calculate the beta-order envelope spectrum (PES) of the filtered signal s, where...

[0079]

[0080] The filtered signal s is obtained by convolving the signal x with the filter h;

[0081] Step 8: Within the β-order envelope spectrum tolerance bandwidth α i -w:α i Find the local maximum within the range of +w and update the loop order α. i+1 .

[0082]

[0083] in k = 1, ..., M is the tolerance bandwidth of the k-th harmonic of the fault order fundamental frequency.

[0084] Step 9: Generate the weighted correlation matrix W of the signal using the beta-order envelope spectrum, where:

[0085]

[0086]

[0087] Step 10: Update the filter coefficients using the Rayleigh quotient. The Rayleigh quotient expression for the generalized cyclostationary index is as follows:

[0088]

[0089] in,

[0090] R XX =x H x (14)

[0091] R XWX =x H Wx (15)

[0092] Where, x H Let x be the Hermitian matrix.

[0093] Step 11: Stop iterative updates by setting the maximum number of blind deconvolution iterations Nmax1 and the convergence condition ε1; otherwise, continue repeating steps 4 through 10.

[0094] Step 12: By convolving the signal x with the filter h, the filtered signal s is obtained, and the beta-order envelope spectrum of the filtered signal s is generated, along with the corresponding index values ​​and shape parameters, and saved.

[0095] Step 13: Analyze all vibration signals through steps 1 to 12 to obtain their β-order envelope spectrum and generate a waterfall plot.

[0096] The figure shows a case of high-speed shaft rear bearing wear failure in a certain unit. In May 2017, the vibration signal waveform showed significant changes. Personnel from the bearing manufacturer inspected the bearing and found metallic impurities in the grease. Touching the rolling area revealed localized peeling. The unit was operated with limited power until June 17, 2017, when replacement was performed and peeling was found in the rolling area of ​​the bearing. The figure shows the index values ​​and shape parameters of each filtered signal obtained using the proposed method, and the figure also shows the β-order envelope spectrum waterfall plot of each filtered signal obtained using the proposed method. It can be seen that the failure trend began in December 2016.

[0097] As can be seen from the method of this invention, compared with the prior art, this invention, based on the Gaussian distribution assumption of the vibration signal, predicts the existence of non-Gaussian distributed vibration signals through calculation methods. This has significant advantages over the prior art.

[0098] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. For those skilled in the art, the present invention can have various modifications and variations. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.< / r> < / r>

Claims

1. A mechanical fault diagnosis method based on maximizing a generalized cyclostationary index, characterized by, The method comprises the following steps: Step 1: collecting a vibration signal and pre-processing the vibration signal to meet the processing requirements; Step 2: setting a filter, inputting the pre-processed vibration signal into the filter, and performing iterative calculation; Step 2.1 of the iterative calculation is to pass the loop order a i Calculate the synchronous average r of the envelope signal |x| and normalize the signal to obtain x t , Wherein: in formula (1), <r>r is a first moment of r;< / r> Wherein, M is the number of sampling points of the cyclic order α, and the relationship between the cyclic order and the sampling order fs is M = fs / α, so Step 2.2 of the iterative calculation is to calculate the shape parameter of the signal x under the generalized Gaussian cyclic stationary model by using an iterative method, wherein: t t |> ;​​ (2) initializing the shape parameter, specifically: (3) iteratively updating the shape parameter, specifically: Wherein, Wherein, ψ and ψ' are digamma and trigamma functions respectively; Step 3: comparing the iteration number and the convergence condition with the set maximum iteration number Nmax1 and the convergence condition ε1 of the blind deconvolution to determine whether to stop the iterative update; Step 4: under the condition of stopping the iterative update, obtaining the filtered signal s, generating the β-order envelope spectrum of the filtered signal s for analysis, and saving the corresponding index value and the shape parameter.

2. The mechanical fault diagnosis method based on maximization of a generalized cyclostationarity index according to claim 1, characterized in that, In step 2, the parameters used for setting the filter; specifically including the filter length L; the maximum iteration number Nmax1 and Nmax2 of blind deconvolution and shape parameter β estimation; the convergence conditions ε1 and ε2 of blind deconvolution and shape parameter β estimation; the theoretical cyclic frequency α0 and the tolerance bandwidth w.

3. The mechanical fault diagnosis method based on maximization of a generalized cyclostationarity index according to claim 1, characterized in that, In step 4, the filtered signal s is obtained by convolving the signal x with the filter h. s = x * h (8).

4. The mechanical fault diagnosis method based on maximization of a generalized cyclostationarity index according to claim 1, characterized in that, In step 4, the β-order envelope spectrum of the filtered signal s is generated, wherein The filtered signal s is obtained by convolving the signal x with the filter h.

5. The mechanical fault diagnosis method based on maximization of a generalized cyclostationarity index according to claim 4, characterized in that, The beta order envelope spectrum looks for local maxima within its tolerance bandwidth a i -w: a i +w and updates the cycle frequency a i+1 : wherein is the tolerance bandwidth for the kth harmonic of the fault frequency fundamental.

6. The mechanical fault diagnosis method based on maximization of a generalized cyclostationarity index according to claim 4, characterized in that, The β-order envelope spectrum generates a weighted correlation matrix W of the signal, wherein: The calculation includes generating a beta order envelope spectrum.

7. The mechanical fault diagnosis method based on maximization of a generalized cyclostationarity index according to claim 6, characterized in that, After generating the weighted correlation matrix W of the signal, the coefficients of the filter are updated by using the Rayleigh quotient, and the Rayleigh quotient of the generalized cyclic stationary index is expressed as: Wherein, R XX = x H x (14) R XWX = x H Wx (15) where x H is a Hermitian matrix of x.

8. The method of claim 1, wherein, In step 3, the iteration number and the convergence condition are compared with the set maximum iteration number Nmax1 and the convergence condition ε1 of the blind deconvolution to determine whether to update the filter through step 2 again or to stop the iterative update.

Images