A wind turbine gearbox vibration signal fault diagnosis method

Abstract

The present invention relates to a fault diagnosis method for a vibration signal of a wind turbine gearbox, comprising: obtaining an original vibration signal; constructing a Hankel matrix according to the original vibration signal, performing singular value decomposition on the Hankel matrix to obtain a singular value sequence matrix; determining a hard threshold of the singular value sequence matrix, and processing the singular value sequence matrix based on the hard threshold using a mixed threshold singular value decomposition denoising algorithm; restoring the processed singular value sequence matrix to a Hankel matrix, and then restoring it to a vibration signal after denoising; filtering the vibration signal after denoising; performing envelope analysis on the filtered signal to obtain a fault characteristic frequency. The present invention adopts a hard threshold adaptive selection algorithm of singular value decomposition, which solves the problems of strong subjectivity, weak adaptability, and easy loss of signal features in the threshold selection of the traditional hard threshold singular value decomposition denoising method, and aims to achieve effective extraction of wind turbine gearbox fault features in a strong noise environment.

Description

Technical Field

[0001] The invention relates to the field of mechanical equipment status monitoring, and in particular to a method for diagnosing vibration signal faults of a wind power gearbox. Background Art

[0002] Wind turbines work in harsh environments all year round, and most of the equipment is installed at high altitudes. Once a failure occurs, it may lead to serious production accidents. Gears, as the main components of wind power transmission systems, often have a high failure rate. Therefore, timely and effective detection and diagnosis of the operating status of gears is of great significance to the normal operation of wind power transmission systems.

[0003] Existing gear diagnosis methods are mainly based on vibration signals, using signal processing technology and machine learning algorithms to diagnose gear component faults. Commonly used machine learning algorithms include support vector machines, random forests, neural networks, etc. The machine learning algorithm model takes a long time to calculate, converges slowly, has weak anti-noise performance, and requires a large data sample. These factors limit the application effect of machine learning methods in gear fault diagnosis. Common signal processing technologies include signal decomposition and reconstruction, signal denoising algorithms, such as variational mode decomposition, empirical mode decomposition, and wavelet decomposition. Such decomposition and denoising methods have problems such as endpoint effects, modal aliasing, and poor adaptability to signal types. In addition, the unit operating environment is noisy and the gear speed is generally low, resulting in a large attenuation of the vibration signal of the internal fault excitation of the gearbox when it is transmitted to the surface of the box. Fault signals with periodic and impact characteristics are difficult to detect in the time and frequency domains, limiting the accuracy of diagnosis. Summary of the invention

[0004] In view of the deficiencies in the prior art, the present invention provides a wind turbine gearbox vibration signal fault diagnosis method, aiming to achieve effective extraction of wind turbine gearbox fault characteristics in a strong noise environment.

[0005] The technical solution adopted by the present invention is as follows:

[0006] A wind turbine gearbox vibration signal fault diagnosis method, comprising:

[0007] S1, obtaining the original vibration signal;

[0008] S2. Construct a Hankel matrix based on the original vibration signal, and perform singular value decomposition on the Hankel matrix to obtain a singular value sequence matrix diag (σ1, σ2, σ3…σ N ), σ N is the Nth singular value;

[0009] S3, determining a hard threshold K of a singular value sequence matrix, and processing the singular value sequence matrix using a mixed threshold singular value decomposition denoising algorithm based on the hard threshold K;

[0010] The hard threshold selection algorithm is used to determine the hard threshold of the singular value sequence matrix, including:

[0011] With a continuously decreasing screening value Y i As standard, Y i =Y i-1 -iλ, i is the number of attenuation, λ is the attenuation step;

[0012] For each reference value Y i , the statistics of this increase is greater than Y i The number of singular values ​​C i ;

[0013] If at the nth time, C i The value is greater than the set value U, and all C i (i=1, 2…n-1) added together is the hard threshold K;

[0014] S4, restoring the processed singular value sequence matrix into a Hankel matrix, and then restoring it into a vibration signal after noise reduction;

[0015] S5, filtering the vibration signal after noise reduction;

[0016] S6. Perform envelope analysis on the filtered signal to obtain the fault characteristic frequency.

[0017] Further technical solutions are:

[0018] In step S3, the process of processing the singular value sequence matrix using a mixed threshold singular value decomposition denoising algorithm based on the hard threshold K includes:

[0019] Determine the singular value σ corresponding to the hard threshold K K , for the singular value sequence diag(σ1,σ2,σ3…σ N ) is greater than σ K The singular values ​​of K The singular values ​​σ i Using the optimal singular value shrinker H * Perform shrinkage processing, the optimal singular value shrinker H * The expression is as follows:

[0020]

[0021] Where K+1≤i≤N, and D is the ratio of the number of rows to the number of columns of the Hankel matrix.

[0022] In step S2, the Hankel matrix is ​​constructed according to the original vibration signal, and the Hankel matrix is ​​subjected to singular value decomposition to obtain a singular value sequence matrix diag (σ1, σ2, σ3…σ N ),include:

[0023] S21, the original vibration signal sequence (x1, x2, x3...x N ), take x1, x2, x3…x m As the first row of the matrix, x2, x3, x4, ... x m+1 As the second row of the matrix, and so on, until an m×m-dimensional Hankel matrix is ​​obtained, where m<0.5N;

[0024] S22, divide the constructed Hankel matrix by c is the noise intensity;

[0025] S23. Perform singular value decomposition on the Hankel matrix: Hankel = Pdiag (σ1, σ2, σ3…σ N )Q T , we get three matrices, Hankel is the Hankel matrix, diag(σ1, σ2, σ3…σ N ) is the singular value sequence matrix, P is the left singular value matrix, Q T is the right singular value matrix.

[0026] In step S4, the processing of the singular value sequence matrix to a Hankel matrix includes:

[0027] The processed singular value sequence matrix Multiply by the left and right singular value matrices, and multiply by Restore to Hankel matrix: Hankel * is the reduced Hankel matrix, is the Nth singular value.

[0028] In step S5, filtering the noise-reduced vibration signal includes:

[0029] The kurtosis index K v As the evaluation basis, a sliding window with a length of l is used, and each time a step is slid, the signal sequence in the sliding window is filtered through the suppression function until it slides to the end of the sequence, and the kurtosis index K is calculated. v If it does not converge, update the parameters of the suppression function and continue the above operation until K v Convergence or reaching the maximum number of iterations stops, K v The expression is as follows:

[0030]

[0031] Among them, x N (n) is the nth sequence point in the vibration signal sequence after noise reduction, is the average amplitude of the entire vibration signal sequence, and N is the length of the vibration signal after noise reduction;

[0032] The expression of the suppression function is as follows:

[0033]

[0034] Among them, x l (n) is the nth sequence point of the signal sequence in the sliding window; |x l (n)| represents x l (n), max(|x l (n)|) represents |x l (n)|; the subscript l represents the sliding window length, l = N / 8;

[0035] (σmax(|x l (n)|)-|x l (n)|,0) + =max(σmax(|x l (n)|)-|x l (n)|,0);

[0036] σ is the attenuation coefficient, σ=1-1 / (1+e -0.5k+2.5 ), k represents the kth iteration.

[0037] Step S6 specifically includes: performing Hilbert transform on the filtered signal to obtain the modulus of the analytical signal, and then obtaining the Fourier spectrum of the analytical signal to obtain the fault characteristic frequency.

[0038] The beneficial effects of the present invention are as follows:

[0039] 1. The present invention adopts a new hard threshold adaptive selection algorithm of singular value decomposition, which solves the problems of strong subjectivity, weak adaptability and easy loss of signal features in the traditional hard threshold singular value decomposition denoising method.

[0040] 2. The present invention combines the traditional hard threshold denoising algorithm with a non-equal optimal shrinkage soft threshold singular value decomposition denoising method, which significantly improves the denoising effect of the signal.

[0041] 3. The present invention proposes an amplitude suppression algorithm for highlighting the periodic fault impact characteristics of the gear signal after noise reduction, and uses envelope analysis to significantly and effectively extract the fault characteristic frequency of the faulty gear.

[0042] Other features and advantages of the present invention will be set forth in the description which follows, and in part will be apparent from the description, or may be learned by practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 It is a schematic diagram of the process of the present invention.

[0044] Figure 2 The time domain waveform and envelope spectrum of the original vibration signal extracted by the embodiment of the present invention.

[0045] Figure 3 The figure shows the time domain waveform and envelope spectrum of the processed signal according to the embodiment of the present invention. DETAILED DESCRIPTION

[0046] The specific implementation of the present invention is described below with reference to the accompanying drawings.

[0047] See also Figure 1 , a wind turbine gearbox vibration signal fault diagnosis method of this embodiment includes:

[0048] S1. Get the original vibration signal:

[0049] For a fixed-axis gearbox, the characteristic frequency of a local gear fault is the rotation frequency of the faulty gear. The speed and number of teeth of the driving wheel are Z1 (r / min) and N1 respectively, and its fault characteristic frequency is Z1 / 60. The number of teeth of the driven wheel is N2, and its fault characteristic frequency is Z1N1 / (60N2).

[0050] S2. Construct a Hankel matrix based on the original vibration signal, and perform singular value decomposition on the Hankel matrix to obtain a singular value sequence matrix diag (σ1, σ2, σ3…σ N ), σ N is the Nth singular value, including:

[0051] S21, the original vibration signal sequence (x1, x2, x3...x N ), take x1, x2, x3…x m As the first row of the matrix, x2, x3, x4, ... x m+1 As the second row of the matrix, and so on, until an m×m-dimensional Hankel matrix is ​​obtained, m<0.5N, specifically m=2000;

[0052] S22, divide the constructed Hankel matrix by c is the noise intensity. In a strong noise environment, the standard deviation of the original vibration signal can be approximately regarded as c;

[0053] S23. Perform singular value decomposition on the Hankel matrix: Hankel = Pdiag (σ1, σ2, σ3…σ N )Q T , we get three matrices, Hankel is the Hankel matrix, diag(σ1, σ2, σ3…σ N ) is the singular value sequence matrix, P is the left singular value matrix, Q T is the right singular value matrix.

[0054] S3, determining a hard threshold K of a singular value sequence matrix, and processing the singular value sequence matrix using a mixed threshold singular value decomposition denoising algorithm based on the hard threshold K;

[0055] The hard threshold selection algorithm is used to determine the hard threshold of the singular value sequence matrix, including:

[0056] With a continuously decreasing screening value Y i As standard, Y i =Y i-1 -iλ, i is the number of attenuation, λ is the attenuation step;

[0057] For each reference value Y i , the statistics of this increase is greater than Y i The number of singular values ​​C i ;

[0058] If at the nth time, C i The value is greater than the set value U, and all C i (i=1, 2…n-1) added together is the hard threshold K; in this embodiment, U=5 is preferred;

[0059] The method of processing the singular value sequence matrix by using a mixed threshold singular value decomposition denoising algorithm based on the hard threshold K includes:

[0060] Determine the singular value σ corresponding to the hard threshold K K , for the singular value sequence diag(σ1,σ2,σ3…σ N ) is greater than σ K The singular values ​​of K The singular values ​​σ i Using the optimal singular value shrinker H * Perform shrinkage processing, the optimal singular value shrinker H * The expression is as follows:

[0061]

[0062] Wherein, K+1≤i≤N, D is the ratio of the number of rows to the number of columns of the Hankel matrix, and in this embodiment, D=1.

[0063] The singular value sequence matrix after comprehensive processing by the above hybrid threshold method is as follows:

[0064] diag(σ1,σ2…σ K , H * (σ K+1 ), H * (σ K+2 )…H * (σ N ))

[0065] The hard threshold selection algorithm proposed in this embodiment, the selected hard threshold can be used as an effective basis for processing the singular value sequence matrix, the mixed threshold singular value decomposition denoising algorithm proposed in this embodiment is equivalent to combining the hard threshold denoising algorithm with a non-equal optimal shrinkage soft threshold singular value decomposition denoising method, so that the singular values ​​in the first half of the hard threshold correspond to the periodic main components of the signal, and the smaller and more slowly changing singular values ​​in the second half mostly correspond to the noise components of the signal, thereby significantly improving the noise reduction effect.

[0066] S4, restoring the processed singular value sequence matrix into a Hankel matrix, and then restoring it into a vibration signal after noise reduction, specifically including:

[0067] The processed singular value sequence matrix diag(σ1, σ2…σ K , H * (σ K+1 ), H * (σ K+2 )…H * (σ N )) multiplied by the left and right singular value matrices, and multiplied by Restore to Hankel matrix:

[0068]

[0069] Among them, Hankel * is the reduced Hankel matrix, is the Nth singular value, and the processed singular value sequence matrix diag(σ1, σ2…σ K , H * (σ K+1 ), H * (σ K+2 )…H * (σ N )) can be written as

[0070] S5, filtering the vibration signal after noise reduction, specifically including:

[0071] The kurtosis index K vAs the evaluation basis, a sliding window with a length of l is used, and each time a step is slid, the signal sequence x in the sliding window is l (n) Perform filtering operation through the suppression function until sliding to the end of the sequence and calculating the kurtosis index K v If it does not converge, update the parameters of the suppression function and continue the above operation until K v Convergence or reaching the maximum number of iterations stops, K ν The expression is as follows:

[0072]

[0073] Among them, x N (n) is the nth sequence point in the vibration signal sequence after noise reduction, is the average amplitude of the entire vibration signal sequence, and N is the length of the vibration signal after noise reduction;

[0074] The expression of the suppression function is as follows:

[0075]

[0076] Among them, x l (n) is the nth sequence point of the signal sequence in the sliding window; |x l (n)| represents x l (n), max(|x l (n)|) represents |x l (n)|; the subscript l represents the sliding window length, l = N / 8;

[0077] (σmax(|x l (n)|)-|x l (n)|,0) + =max(σmax(|x l (n)|)-|x l (n)|,0);

[0078] σ is the attenuation coefficient, σ=1-1 / (1+e -0.5k+2.5 ), k represents the kth iteration.

[0079] σ is the parameter, and the kurtosis index K is calculated v If it does not converge, update σ and continue sliding and calculating until K v Converge or meet the maximum number of iterations.

[0080] Compared with the traditional kurtosis, the kurtosis index used in this embodiment is a dimensionless index reflecting the degree of deviation of the signal from the normal distribution. It is particularly sensitive to the impact component in the signal, and its value is independent of the input shaft speed, the size of each gearbox component and the load.

[0081] The amplitude suppression algorithm proposed in this embodiment can be used to suppress the interference of noise and clutter around the impact signal segment after singular value decomposition denoising, and better highlight the periodic fault impact characteristics of the gear signal. Compared with some traditional feature highlighting methods such as the Teager energy operator, this method has better fidelity for the signal impact part, and can better highlight the impact part and suppress the interference of noise and clutter.

[0082] S6. Perform envelope analysis on the filtered signal to obtain the fault characteristic frequency, including:

[0083] Perform Hilbert transform on the filtered signal to find the modulus of the analytical signal, and then find the Fourier spectrum of the analytical signal to obtain the fault characteristic frequency.

[0084] The technical solution of the present application is explained below in conjunction with specific test data.

[0085] A wind turbine gearbox vibration signal fault diagnosis method comprises the following steps:

[0086] S1. Obtain the vibration signal of the wind turbine gearbox:

[0087] The data was collected by the wind power transmission test bench of the National Engineering Research Center for Safe Operation and Intelligent Measurement and Control of Large Power Generation Equipment, School of Energy and Environment, Southeast University. The entire device is mainly composed of a transmission test bench, an acceleration sensor, a data acquisition system, and a computer. Among them, the test bench part includes a drive motor, a main shaft, a torque speed sensor, a single-stage cylindrical gearbox, a coupling, and a load motor. Different speeds (50-1000r / min) are obtained by adjusting the variable frequency drive motor. The sampling frequency is 10kHz, the gear speed is 800r / min, the number of teeth is 18, the fault gear is the driven wheel, the fault type is a broken tooth fault, the number of gear teeth is 58, and the fault characteristic frequency of the faulty gear is calculated to be 4.13Hz. A continuous 8000 data points are randomly selected from the 234827 data points obtained by sampling, that is, the length of the gear vibration signal to be analyzed is N=8000. At the same time, in order to simulate the strong noise environment in the actual work of wind power, a certain intensity of Gaussian noise is added to the vibration signal. Figure 2 (a) and (b) are the time domain waveform and envelope spectrum of the gear after adding noise. It can be seen that due to the interference of noise, it is impossible to observe obvious fault impact characteristics from the time domain diagram, and the fault characteristic frequency of the gear cannot be observed in the envelope spectrum.

[0088] S2. Based on the processed gearbox vibration signal, a 2000×2000 Hankel matrix is ​​constructed, where the number of matrix rows is m=2000, the noise variance of the signal is estimated, and the estimated value is c=1. The Hankel matrix is ​​subjected to singular value decomposition to obtain a singular value sequence matrix.

[0089] S3. Using a hard threshold selection algorithm, the selected hard threshold K=121, based on the selected hard threshold, using the hybrid threshold singular value decomposition denoising algorithm to process the singular value sequence.

[0090] S4, the processed singular value sequence matrix is ​​restored to the Hankel matrix and then multiplied by Finally, it is restored to the vibration signal after noise reduction.

[0091] S5, filtering the denoised signal through an amplitude suppression filter;

[0092] In this embodiment, the maximum number of iterations k is set to 30. When k = 12, the kurtosis index K ν It starts to converge, and finally takes K = 30 ν =The signal corresponding to 29.313.

[0093] S6. Perform envelope analysis on the processed signal to obtain the fault characteristic frequency.

[0094] Figure 3 (a) and (b) are the time domain waveform and envelope spectrum of the processed fault gear signal, respectively. Figure 3 In (a), the impact signal component is observed, and the interval between the impact signals is about 0.25s. The calculated frequency is: 1 / 0.25=4Hz, which is basically consistent with the fault characteristic frequency of the faulty gear, 4.13Hz. Figure 3 In (b), the fault characteristic frequency of 4.13 Hz and its frequency harmonic components can be clearly observed, which proves that the analysis method of the present application can effectively diagnose faults of wind turbine gears.

[0095] Those skilled in the art can understand that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention is described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A wind turbine gearbox vibration signal fault diagnosis method, characterized in that: include: S1, obtaining the original vibration signal; S2. Construct a Hankel matrix based on the original vibration signal, and perform singular value decomposition on the Hankel matrix to obtain a singular value sequence matrix diag(σ1,σ2,σ3…σ N ), σ N is the Nth singular value; S3, determining a hard threshold K of a singular value sequence matrix, and processing the singular value sequence matrix using a mixed threshold singular value decomposition denoising algorithm based on the hard threshold K; The hard threshold selection algorithm is used to determine the hard threshold of the singular value sequence matrix, including: With a continuously decreasing screening value Y i As standard, Y i =Y i-1 -iλ, i is the number of attenuation, λ is the attenuation step; For each reference value Y i , the statistics of this increase is greater than Y i The number of singular values ​​C i ; If at the nth time, C i The value is greater than the set value U, and all C i Add them together to get the hard threshold K, where i = 1, 2…, n-1; S4, restoring the processed singular value sequence matrix into a Hankel matrix, and then restoring it into a vibration signal after noise reduction; S5, filtering the vibration signal after noise reduction; S6. Perform envelope analysis on the filtered signal to obtain the fault characteristic frequency; In step S3, the process of processing the singular value sequence matrix using a mixed threshold singular value decomposition denoising algorithm based on the hard threshold K includes: Determine the singular value σ corresponding to the hard threshold K K , for the singular value sequence diag(σ1,σ2,σ3…σ N ) is greater than σ K The singular values ​​of K The singular values ​​σ i Using the optimal singular value shrinker H * Perform shrinkage processing, the optimal singular value shrinker H * The expression is as follows: Where K+1≤i≤N, and D is the ratio of the number of rows to the number of columns of the Hankel matrix.

2. The wind turbine gearbox vibration signal fault diagnosis method according to claim 1, characterized in that: In step S2, the Hankel matrix is ​​constructed according to the original vibration signal, and the singular value decomposition of the Hankel matrix is ​​performed to obtain a singular value sequence matrix diag (σ1, σ2, σ3…σ N ),include: S21, the original vibration signal sequence (x1, x2, x3…x N ), take x1,x2,x3…x m As the first row of the matrix, x2,x3,x4,…x m+1 As the second row of the matrix, and so on, until an m×m-dimensional Hankel matrix is ​​obtained, where m<0.5N; S22, divide the constructed Hankel matrix by c is the noise intensity; S23. Perform singular value decomposition on the Hankel matrix: Hankel = Pdiag (σ1, σ2, σ3…σ N )Q T , we get three matrices, Hankel is the Hankel matrix, diag(σ1,σ2,σ3…σ N ) is the singular value sequence matrix, P is the left singular value matrix, Q T is the right singular value matrix.

3. The wind turbine gearbox vibration signal fault diagnosis method according to claim 2, characterized in that: In step S4, the processing of the singular value sequence matrix to a Hankel matrix includes: The processed singular value sequence matrix Multiply by the left and right singular value matrices, and multiply by Restore to Hankel matrix: Hankel * is the reduced Hankel matrix, is the Nth singular value.

4. The wind turbine gearbox vibration signal fault diagnosis method according to claim 1, characterized in that: In step S5, filtering the noise-reduced vibration signal includes: The kurtosis index K ν As the evaluation basis, a sliding window with a length of l is used, and each time a step is slid, the signal sequence in the sliding window is filtered through the suppression function until it slides to the end of the sequence, and the kurtosis index K is calculated. ν If it does not converge, update the parameters of the suppression function and continue the above operation until K ν Convergence or reaching the maximum number of iterations stops, K ν The expression is as follows: Among them, x N (n) is the nth sequence point in the vibration signal sequence after noise reduction, is the average amplitude of the entire vibration signal sequence, and N is the length of the vibration signal after noise reduction; The expression of the suppression function is as follows: Among them, x l (n) is the nth sequence point of the signal sequence in the sliding window; |x l (n)| represents x l (n), max(|x l (n)|) represents |x l (n)|; the subscript l represents the sliding window length, l = N / 8; (σmax(|x l (n)|)-|x l (n)|,0) + =max(σmax(|x l (n)|)-|x l (n)|,0); σ is the attenuation coefficient, σ=1-1 / (1+e -0.5k+2.5 ), k represents the kth iteration.

5. The wind turbine gearbox vibration signal fault diagnosis method according to claim 1, characterized in that: Step S6 specifically includes: performing Hilbert transform on the filtered signal to obtain the modulus of the analytical signal, and then obtaining the Fourier spectrum of the analytical signal to obtain the fault characteristic frequency.

Images