A wind turbine generator set rotating bearing fault diagnosis method based on high-order statistics analysis
By employing high-order statistical analysis and phase space reconstruction, the problems of noise interference and nonlinear signal identification in the fault diagnosis of slewing bearings of wind turbine generators were solved, achieving high-precision fault diagnosis in complex environments.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2025-03-18
- Publication Date
- 2026-03-27
AI Technical Summary
Fault diagnosis of slewing bearings in wind turbine generator sets is difficult to accurately identify in complex noise environments, especially because weak vibration signals are masked by noise and nonlinear, non-Gaussian signal characteristics are complex, making it difficult for traditional methods to capture these characteristics.
We employ high-order statistical analysis, using skewness and kurtosis as statistical filters to denoise and suppress background noise interference. We utilize phase space reconstruction and adaptive multispectral estimation to extract fault features in high-dimensional space, and combine gradient descent to optimize parameters, thereby improving diagnostic accuracy.
It effectively suppresses noise interference, enhances the expression of fault characteristics, improves the accuracy and precision of slewing bearing fault diagnosis, and can identify potential fault modes in complex environments.
Smart Images

Figure CN120275046B_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The method belongs to the field of wind power generation, and particularly relates to a wind turbine generator set slewing bearing fault diagnosis method based on high-order statistics analysis. BACKGROUND
[0002] As a sustainable energy, wind energy has become one of the important sources of global electricity production. Wind turbines extract wind energy from the wind flow in the atmosphere, and the slewing bearings in the wind turbine generator set are divided into two types: variable pitch bearings and yaw bearings. The variable pitch bearing plays a key role in the process of the wind turbine generator set extracting wind energy from the wind flow in the atmosphere, and it is responsible for adjusting the blade pitch angle to obtain the optimal output power of the wind turbine. The yaw bearing is responsible for adjusting the direction of the wind turbine generator set to the wind direction to capture the maximum wind energy. Due to the complex working environment of the slewing bearing and the long-term state of load change, the bearing is prone to failure. The slewing bearing fault is usually manifested as a weak vibration signal. When running at low speed, the signal characteristics related to the fault are obscured due to sliding, vibration, noise and interference of other machine parts, so it is very difficult to identify the slewing bearing characteristic frequency from the vibration signal, which leads to difficulty in accurately diagnosing the fault; in addition, the bearing fault signal has complex and time-varying dynamic characteristics, which makes it impossible for traditional methods based on linear or Gaussian assumptions to capture these characteristics.
[0003] As an effective signal processing technology, high-order statistics performs well in the analysis of non-linear and non-Gaussian signals, and it can provide more information than traditional linear analysis methods. Therefore, how to use high-order statistics analysis to diagnose non-linear and non-Gaussian vibration signals in a complex noise environment is still a technical problem to be solved. SUMMARY
[0004] The method proposes a wind turbine generator set slewing bearing fault diagnosis method based on high-order statistics analysis. Firstly, in view of the fact that the weak fault signal of the bearing is obscured by noise, skewness and kurtosis are used as the measurement standard for statistical filter denoising to suppress the interference of background noise and enhance the expression of fault characteristics. The denoised signal will be used as the input of the feature extraction module for subsequent fault diagnosis. Secondly, in view of the non-linear and non-Gaussian characteristics of the fault signal, the bearing fault caused deep features are reflected in the high-dimensional space by using phase space reconstruction, and adaptive multi-spectrum estimation is used in the high-dimensional space to analyze the signal for fault diagnosis.
[0005] The technical scheme adopted by the method is as follows:
[0006] S1: Obtain the vibration signal of the slewing bearing and set it as the original signal;
[0007] S2: The original signal is segmented for processing, and the number of segments is selected as an integer multiple of the number of signal periods, to ensure that each segment of the signal contains a complete signal period and avoids the loss of periodic characteristics. The skewness and kurtosis of each segment of the signal are calculated in turn;
[0008] S3: In view of the problem that the fault signal of the slewing bearing is weak and is easily affected by background noise, the initial threshold of skewness and kurtosis is set according to the statistical distribution difference of different fan data. The skewness value is used to measure the asymmetry of the signal in statistics, and the kurtosis value is used to measure the sharpness of the signal distribution. The normal bearing vibration signal presents a stable and symmetric distribution, the skewness value is close to 0, and the kurtosis value is close to 3. The fault bearing vibration signal produces impact and asymmetric fluctuation due to defects, the skewness value deviates from 0, and the kurtosis value is greater than 3. A statistical filter is designed by using skewness and kurtosis in high-order statistical analysis, so as to suppress the interference of background noise.
[0009] S4: The skewness value and kurtosis value of each segment of the signal are compared and calculated with the initial threshold of skewness and kurtosis to obtain a weight coefficient. When the calculated value is greater than 1, it is considered that the signal segment contains fault characteristics, and the signal segment is highlighted. Otherwise, the background noise is suppressed to interfere with fault diagnosis, and then the denoised vibration signal is obtained;
[0010] S5: The denoised signal is taken as the input data of the next step for phase space reconstruction. By using the time delay embedding method, the one-dimensional time series is converted into expression in high-dimensional phase space, to represent the dynamic characteristics of the signal in high-dimensional space. In the reconstructed phase space, adaptive multi-spectrum estimation is performed on the signal to reveal the phase information and high-order nonlinear characteristics of the signal, to extract the interaction relationship between frequencies in different dimensions, to identify the potential fault mode in the slewing bearing, and to obtain the data diagnosis result;
[0011] S6: The error between the data diagnosis result and the actual state result is recorded as ε, the error requirement between the data diagnosis result and the actual state result is set as ε0, and the high-order statistical quantity parameter iteration weight parameter is dynamically updated by using the gradient descent method, so that the error between the data diagnosis result and the actual state result reaches the set error requirement between the data diagnosis result and the actual state result, to gradually approach a minimum error solution, and to improve the accuracy of fault diagnosis. BRIEF DESCRIPTION OF DRAWINGS
[0012] Figure 1 Flow chart of slewing bearing fault diagnosis method for wind turbine generator based on high-order statistical quantity analysis DETAILED DESCRIPTION
[0013] First, install the appropriate vibration sensor on the slewing bearing for signal acquisition. The selection of the sensor should be optimized according to the working characteristics of the bearing. The sensor is required to capture low-frequency signals and high-frequency signals caused by friction and impact factors during the operation of the slewing bearing, and has sufficient sensitivity to capture small vibration signals. A sensor with a wide dynamic range is used to cope with signals from low to high vibration intensity to ensure that the collected data contains various dynamic information of the bearing operation, and the vibration situation can be reflected at each stage of the slewing bearing operation. The collected original signal is x(t), which contains fault signals and noise signals. Further, the specific steps of data denoising according to the statistical characteristics of background noise and fault signals are as follows:
[0014] First step: calculate the skewness and kurtosis values of the normal state bearing vibration signal:
[0015]
[0016] In the formula
[0017] x represents the vibration signal;
[0018] S represents the signal skewness value;
[0019] E represents the mathematical expectation;
[0020] μ represents the signal mean value;
[0021] σ represents the signal standard deviation.
[0022]
[0023] In the formula
[0024] K represents the signal kurtosis value.
[0025] According to the statistical distribution of skewness and kurtosis values, set threshold values S0 and K0.
[0026] Second step: segment the original signal with a length of L, and the number of segments is an integer multiple of the signal period, ensuring that each segment of the signal contains a complete signal period, avoiding the loss of periodic characteristics, and calculating the skewness and kurtosis of each segment. The denoising weight coefficient is automatically adjusted according to the calculation results of skewness and kurtosis, and the weight parameter is defined as:
[0027]
[0028] In the formula
[0029] ω i is the weight coefficient of the ith segment of the signal;
[0030] K i is the kurtosis value of the ith segment of the signal;
[0031] S i is the skewness value of the i-th signal segment;
[0032] If the skewness and kurtosis exceed the initial threshold values L0 and S0, it is considered that the signal segment contains fault characteristics, and the weight coefficient is multiplied by the signal segment to highlight the fault characteristics; if the skewness and kurtosis are lower than the initial threshold values K0 and S0, it is considered that the signal segment belongs to noise, and median filtering is performed on the signal segment, and the filtered signal of the i-th segment is is expressed as:
[0033]
[0034] wherein
[0035] x i (t) is the i-th original signal segment without processing;
[0036] is the i-th denoised vibration signal segment;
[0037] MF represents a median filter.
[0038] Further, the specific steps of phase space reconstruction are as follows:
[0039] First step: selecting the time series data after denoising
[0040] Second step: calculating a suitable time delay τ to determine the interval of elements of each new state vector, observing the trend of the autocorrelation function of the signal with the change of the time delay, finding the time delay corresponding to the first significant decrease as the suitable time delay, so that the adjacent state points do not coincide in the phase space as much as possible, while maintaining the dynamic characteristics of the system, and the calculation formula of the autocorrelation function is:
[0041]
[0042] wherein
[0043] R(τ) is the autocorrelation function;
[0044] L is the length of the signal.
[0045] R(τ) reflects the degree of similarity between the signal after time delay τ and the original signal, and the autocorrelation functions of different time delays τ are calculated in turn, and the position where the autocorrelation function first decays to a specified threshold (set to 1 / e) is found, that is, the time point where the signals lose linear correlation.
[0046] Step 3: Select the embedding dimension, which represents how many independent variables are needed to describe the dynamics of the system in the phase space. As the embedding dimension increases, the local structure of the data points, i.e., the relative distances between data points, should remain consistent. If the embedding dimension is insufficient, some nearby data points will be incorrectly mapped to more distant regions, resulting in the appearance of false nearest neighbors. When the number of false nearest neighbors is small, it indicates that the embedding dimension of the reconstructed phase space is appropriate, so the number of false nearest neighbors can be monitored to determine whether the embedding dimension is sufficient to preserve the local structure of the data, thereby determining the embedding dimension. In the m-th dimensional space, select a data point x m (t) and its adjacent data point x m (t0), calculate the Euclidean distance between the two:
[0047] d m (t) = ||x m (t) - x m (t0) ||
[0048] where
[0049] d m (t) represents the Euclidean distance between two adjacent data points in the m-dimensional space;
[0050] ||·|| represents the Euclidean norm.
[0051] Embed the data points from the m-dimensional space into the m+1-dimensional space to obtain the new embedding vector:
[0052] x m+1 (t) = [x(t), x(t+τ), …, x(t+(m-1)τ), x(t+mτ)]
[0053] In the m+1-dimensional space, select a data point x m+1 (t) and its adjacent data point x m+1 (t0), calculate the new Euclidean distance between the two adjacent data points:
[0054] d m+1 (t) = ||x m+1 (t) - x m+1 (t0) ||
[0055] where
[0056] d m+1 (t) represents the Euclidean distance between two adjacent data points in the m+1-dimensional space;
[0057] If the data points are close in the low-dimensional space but become far apart in the high-dimensional space, this point can be considered a false nearest neighbor. Assume that in the m-dimensional space, data points x m (t) and xm Euclidean distance d m (t) is small, while in m+1 dimensional space, Euclidean distance d m+1 (t) becomes larger, it indicates that the data point in low dimensional space is likely to be a false neighbor, for each pair of adjacent points, the distance ratio r(t) of them in low dimensional space and high dimensional space is calculated:
[0058]
[0059] When r(t) is greater than 20, it is considered that the point is a false neighbor point.
[0060] The proportion of false neighbor points in the total number of data points in the corresponding dimension is calculated, with the increase of the dimension m of phase space, the proportion of false neighbors gradually decreases until the proportion of false neighbor points tends to be stable, thereby determining the phase space dimension D of the data:
[0061]
[0062] wherein
[0063] D is the appropriate embedding dimension of phase space reconstruction;
[0064] N f (m) is the proportion of false neighbor points in m dimensional space, that is, the number of false neighbor data points;
[0065] N(m) is the total number of data points in m dimensional space.
[0066] Fourth step: data is constructed into a new state vector in phase space according to the delay time and embedding dimension, through phase space reconstruction, the multi-dimensional dynamic characteristics of the original system are reconstructed from single vector, the dynamic information of the original system is recovered, and the nonlinear characteristics of the data are more directly analyzed, providing a basis for subsequent complex system processing.
[0067] Further, an adaptive multi-spectrum estimation method is used to extract features from the signal x D (t) reconstructed in D dimensional space, analyze the frequency components of the signal from different dimensions, and distinguish the frequency spectrum anomalies caused by faults, that is:
[0068]
[0069] wherein
[0070] θ i is a time delay variable;
[0071] ω i is a frequency variable;
[0072] S D (ω1,ω2,…,ω D-1 ) is the signal xD The D-order spectrum of (t);
[0073] C D (θ1,θ2,…,θ D-1 ) is the signal x D The D-1 order cumulant of (t).
[0074] Furthermore, by using multispectral estimation, higher-order information of the signal is obtained, revealing the subtle and rapidly changing characteristics of bearing failure, capturing the interaction between various faults, and identifying diagnostic results.
[0075] Furthermore, the error between the data diagnostic results and the actual state results is denoted as ε, and the error requirement between the data diagnostic results and the actual state results is set to ε0. By continuously updating the higher-order statistical parameters (skewness and kurtosis thresholds) and iterating the weight parameters, the error between the data diagnostic results and the actual state results is made to reach the set error requirement between the data diagnostic results and the actual state results. The specific steps are as follows:
[0076] Step 1: The error between the data diagnostic results and the actual results is expressed as the mean square error (MSE) as the target error function ε.
[0077] Step 2: Calculate the gradient of the error function with respect to the skewness and kurtosis thresholds, and update the skewness and kurtosis values:
[0078]
[0079] In the formula
[0080] S (k+1) This represents the skewness threshold at the (k+1)th iteration.
[0081] S (k) This represents the skewness threshold at the k-th iteration;
[0082] K (k+1) This represents the kurtosis threshold at the (k+1)th iteration.
[0083] K (k) This represents the kurtosis threshold at the k-th iteration;
[0084] γ represents the learning rate;
[0085] This represents the gradient of the error function with respect to skewness;
[0086] This represents the gradient of the error function with respect to kurtosis.
[0087] Step 3: Calculate the new error ε (k+1) When ε (k+1)The iteration process stops when the error requirement between the set data diagnosis result and the actual state result is reached.
[0088] The skewness threshold and the kurtosis threshold are dynamically adjusted according to the change of the error, the de-noising weight coefficient is optimized, the error between the diagnosis result and the actual component operation state result is reduced after multiple iterations, and the fault diagnosis accuracy is improved.
Claims
1. A method for fault diagnosis of slewing bearings in wind turbine generator sets based on high-order statistical analysis, characterized in that: Skewness and kurtosis are used as metrics for statistical filter denoising to suppress background noise interference and highlight fault characteristics. Phase space reconstruction is used to reveal the deep-seated characteristics of bearing faults in a high-dimensional space. Adaptive multispectral estimation is then employed in this high-dimensional space to diagnose characteristic faults in non-Gaussian and nonlinear signals. The adaptive multispectral estimation method is used to reconstruct the signal in D-dimensional space. Feature extraction is performed to analyze the frequency components of the signal from different dimensions, thereby identifying spectral anomalies caused by faults. In the formula It is a time delay variable; It is a frequency variable; It is a signal The D-order spectrum; It is a signal The D-1 order cumulative quantity, Furthermore, through multispectral estimation, higher-order information of the signal is obtained, revealing the subtle and rapidly changing characteristics present during bearing failure, capturing the interaction relationships between various faults, and identifying diagnostic results. Furthermore, the error between the data diagnostic results and the actual state results is denoted as... Set the error requirement between the data diagnostic results and the actual state results to be By continuously updating the higher-order statistical parameters and iterating the weight parameters, the error between the data diagnostic results and the actual state results is made to meet the set error requirements between the data diagnostic results and the actual state results. The specific steps are as follows: Step 1: The error between the data diagnostic results and the actual results is calculated using the mean squared error as the target error function. ; Step 2: Calculate the gradient of the error function with respect to the skewness and kurtosis thresholds, and update the skewness and kurtosis values: In the formula This represents the skewness threshold at the (k+1)th iteration. This represents the skewness threshold at the k-th iteration; This represents the kurtosis threshold at the (k+1)th iteration. This represents the kurtosis threshold at the k-th iteration; Indicates the learning rate; This represents the gradient of the error function with respect to skewness; This represents the gradient of the error function with respect to kurtosis. Step 3: Calculate the new error ,when The iteration process stops when the set error requirement between the data diagnostic results and the actual results is met. The skewness and kurtosis thresholds are dynamically adjusted based on the changes in error, and the denoising weight coefficients are optimized. After multiple iterations, the error between the diagnostic results and the actual operating status of the components is reduced.
2. The method according to claim 1, characterized in that: Skewness and kurtosis are used as metrics for statistical filter denoising to suppress background noise interference and highlight fault characteristics. Specifically, data denoising is performed based on the statistical characteristics of background noise and fault signals, including: Step 1: Calculate the skewness and kurtosis values of the bearing vibration signal under normal conditions: In the formula x represents the vibration signal; S represents the signal skewness value; E represents the expected value; Indicates the signal mean; Indicates the standard deviation of the signal. In the formula K represents the signal kurtosis value. Thresholds are set based on the statistical distribution charts of skewness and kurtosis values. and , Step 2: Divide the original signal of length L into segments. The number of segments should be an integer multiple of the signal period to ensure that each segment contains a complete signal period and avoids the loss of periodic features. Calculate the skewness and kurtosis of each segment. The denoising weight coefficients are automatically adjusted based on the skewness and kurtosis calculation results. The weight parameters are defined as follows: In the formula For the first Weighting coefficients for segment signals; Let be the kurtosis value of the i-th segment of the signal; Let be the skewness value of the i-th segment of the signal; If skewness and kurtosis exceed the initial threshold and If the signal segment contains fault characteristics, the weighting coefficient is multiplied by the signal segment to highlight the fault characteristics; if the skewness and kurtosis are below the initial threshold... and If the signal segment is considered to be noise, median filtering is performed. The filtered signal of segment i is... Expressed as: In the formula This represents the i-th segment of the original, unprocessed signal. The i-th segment is the denoised vibration signal; MF stands for Median Filter.
3. The method according to claim 1, characterized in that: Phase space reconstruction is used to reveal the deeper characteristics of bearing faults in a high-dimensional space. Adaptive multispectral estimation is then employed in this high-dimensional space to diagnose characteristic faults in non-Gaussian and nonlinear signals. The specific steps are as follows: Step 1: Select the denoised time series data ; Step 2: Calculate a suitable time delay The element spacing of each new state vector is determined by observing the trend of the signal autocorrelation function with time delay. The time delay corresponding to the first significant decrease is found as a suitable time delay, so that adjacent state points do not overlap as much in phase space, while maintaining the dynamic characteristics of the system. The formula for calculating the autocorrelation function is: In the formula It is an autocorrelation function; L is the signal length. This reflects the signal's time delay. Then, the similarity between the signal and the original signal is calculated sequentially for different time delays. The autocorrelation function is used to find the point where the autocorrelation function first decays to a specified threshold, i.e., the time point at which the signals lose their linear correlation. Step 3: Select the embedding dimension. The embedding dimension represents how many independent variables are needed in the phase space to describe the system's dynamics. As the embedding dimension increases, the local structure of the data points, i.e., the relative distances between data points, should remain consistent. If the embedding dimension is insufficient, some neighboring data points will be incorrectly mapped to more distant regions, leading to the appearance of false nearest neighbors. When the number of false nearest neighbors is small, it indicates that the embedding dimension for phase space reconstruction is appropriate. Therefore, the number of false nearest neighbors can be monitored to determine whether the embedding dimension is sufficient, ensuring that the local structure of the data is preserved, thereby determining the embedding dimension. In the m-th dimension, select a data point. and its adjacent data points Calculate the Euclidean distance between the two: In the formula Represents the Euclidean distance between two adjacent data points in m-dimensional space; Represents the Euclidean norm. Embedding data points from m-dimensional space into m+1-dimensional space yields a new embedding vector: In the (m+1)th dimension, select a data point and its adjacent data points Calculate the new Euclidean distance between two adjacent data points: In the formula This represents the Euclidean distance between two adjacent data points in m+1 dimensional space; If data points are close in low-dimensional space but become far apart in high-dimensional space, then this point can be considered a false nearest neighbor. Suppose that in the n-th dimension, data points... and Euclidean distance It is very small, and in the +1 dimension space, the Euclidean distance is... If the value becomes larger, it indicates that the data point may be a false nearest neighbor in the low-dimensional space. For each pair of adjacent points, calculate the distance ratio r(t) between the low-dimensional and high-dimensional spaces: when If the value is greater than 20, then the point is considered a false nearest neighbor. Calculate the proportion of false nearest neighbors to the total number of data points in the corresponding dimension. As the dimension m of the phase space increases, the proportion of false nearest neighbors gradually decreases until it stabilizes, thus determining the phase space dimension D of the data. In the formula D is the appropriate embedding dimension for phase space reconstruction; It is the proportion of false nearest neighbors in the m-dimensional dimension, that is, the number of false nearest neighbor data points; N(m) is the total number of data points in dimension m. Step 4: The data is used to construct a new state vector in the phase space based on the delay time and the embedding dimension. Through phase space reconstruction, the multidimensional dynamic characteristics of the original system are reconstructed from a single vector, and the dynamic information of the original system is restored.
Citation Information
Patent Citations
Establishment and using method of database for detecting bearing running-in abnormal noise
CN105716865A
Neural network-based detection method for abnormal target in nonlinear hyperspectral image
CN108648176A