Bearing fault dynamic ridge tracking and feature matching diagnosis method based on high-resolution time-frequency analysis

Abstract

The invention discloses a bearing fault dynamic ridge tracking and feature matching diagnosis method based on high-resolution time-frequency analysis, which is applied to the technical field of fault diagnosis and signal processing of rotating machinery of a wind turbine generator, and comprises the following steps: acquiring an acceleration vibration signal of a bearing, and utilizing an optimized variational mode decomposition algorithm to obtain an acceleration vibration signal of the bearing; decomposing the acceleration vibration signal of the bearing to obtain an optimal intrinsic mode function; transforming the optimal intrinsic mode function component by using a synchronous compression transformation method to obtain a high-resolution time-frequency graph, and extracting a time-frequency ridge line in the high-resolution time-frequency graph by using a dynamic programming algorithm; constructing a fault characteristic curve based on the time-frequency ridge line, and calculating a fault actual characteristic coefficient of the fault characteristic curve; and according to a comparison result of the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type, identifying the fault type of the bearing. The problem that the rotation frequency of the bearing is unstable under variable working conditions is solved, and the fault type of the bearing can be accurately recognized.

Description

Technical Field

[0001] The present invention relates to the technical field of fault diagnosis and signal processing of rotating machinery in wind turbines, and more particularly, to a method for dynamic ridge line tracking and feature matching diagnosis of bearing faults with high-resolution time-frequency analysis. Background Art

[0002] As a typical variable-condition rotating mechanical equipment, the bearings in wind turbines are its core components. Early fault detection of bearings is of crucial significance in ensuring the continuous and stable operation of equipment, improving the overall operation efficiency, and significantly reducing maintenance and operation costs. By accurately identifying potential problems in bearings in a timely manner, the further development of faults can be prevented, the service life of the equipment can be extended, and the reliability and economy of wind turbines can be ensured.

[0003] Traditional fault diagnosis methods generally face challenges such as insufficient diagnostic accuracy and high false positive rates when dealing with complex mechanical equipment operating conditions and diverse fault characteristics. For example, in the traditional variational mode decomposition (VMD) method, the selection of key parameters overly relies on empirical judgment, which easily leads to problems such as mode mixing or over-decomposition, thus affecting the accuracy of feature extraction; classical time-frequency analysis methods such as short-time Fourier transform (STFT) and wavelet transform (WT) are difficult to effectively reveal the transient impact characteristics in weak faults due to the existence of energy leakage; and most current methods mainly compare fault characteristic frequencies by setting fixed thresholds, which is extremely prone to false positives in the case of frequent changes in equipment speed or strong noise interference. Therefore, to improve the accuracy and reliability of early fault detection of mechanical equipment under complex working conditions, more advanced and adaptable fault diagnosis technologies need to be developed.

[0004] Regarding the problems in the related art, no effective solutions have been proposed yet. Summary of the Invention

[0005] To overcome the above problems, the present invention aims to propose a method for dynamic ridge line tracking and feature matching diagnosis of bearing faults with high-resolution time-frequency analysis, aiming to solve the problem of difficult bearing fault identification.

[0006] For this purpose, the specific technical solution adopted by the present invention is as follows:

[0007] A method for dynamic ridge line tracking and feature matching diagnosis of bearing faults with high-resolution time-frequency analysis, the method comprising the following steps:

[0008] S1. Obtain the acceleration vibration signal of the bearing, and decompose the acceleration vibration signal of the bearing using an optimized variational mode decomposition algorithm to obtain the optimal intrinsic mode functions;

[0009] S2. Use the synchronous compression transform method to transform the optimal intrinsic mode function components to obtain a high-resolution time-frequency diagram, and use the dynamic programming algorithm to extract the time-frequency ridge line in the high-resolution time-frequency diagram;

[0010] S3. Construct a fault feature curve based on the time-frequency ridge line and calculate the actual fault feature coefficient of the fault feature curve;

[0011] S4. Identify the bearing fault type according to the comparison result between the actual fault feature coefficient and the fault feature coefficient corresponding to the theoretical fault type.

[0012] Optionally, obtaining the acceleration vibration signal of the bearing and using the optimized variational mode decomposition algorithm to decompose the acceleration vibration signal of the bearing to obtain the optimal intrinsic mode function includes the following steps:

[0013] S11. Initialize the frost individual scale and the maximum number of iterations in the frost ice optimization algorithm, and determine the initial position for each frost individual in the parameter space of the variational mode decomposition algorithm, and set the decomposition layer number K and the penalty factor α;

[0014] S12. Take the minimization of the envelope entropy of the intrinsic mode function component as the optimization goal, and calculate the fitness of each frost individual according to the fitness function;

[0015] The expression of the fitness function is:

[0016] E = -∑p i ·logp i ;

[0017] In the formula, E represents the fitness function; p i represents the envelope probability distribution of the intrinsic mode function component;

[0018] S13. Based on the global exploration of the soft frost search strategy and the hard frost piercing mechanism, perform iterative optimization on the variational mode decomposition; until the termination condition is reached, the iterative optimization ends, and the global optimal position is obtained;

[0019] The expression of the particle position update formula for the global exploration of the soft frost search strategy is:

[0020]

[0021] In the formula, represents the new position of the updated particle after the global exploration of the soft frost search strategy; R best,j represents the jth particle of the best frost agent R in the frost population; r1 represents a random number closely related to the particle movement direction; f represents the current number of iterations; F represents the maximum number of iterations; h represents the adhesion coefficient; Ub ij represents the upper limit of the central distance between two frost particles; Lb ijrepresents the lower limit of the center distance between two frost particles; β represents the environmental factor;

[0022] The expression of the particle position update formula for the hard frost piercing mechanism is:

[0023]

[0024] In the formula, represents the new position of the updated particle in the hard frost piercing mechanism; R best,j represents the j-th particle of the best frost agent R in the frost population; F norm represents the normalized fitness value of the current agent; S i represents the piercing displacement vector of the i-th particle; r3 represents a random number that controls the randomness of the piercing step size;

[0025] S14. According to the global optimal position, combined with the variational mode decomposition algorithm, decompose the acceleration vibration signal of the bearing to obtain the optimal intrinsic mode function;

[0026] The expression of the variational mode decomposition algorithm is:

[0027]

[0028] In the formula, u k represents the K-th mode component; ω k represents the center frequency of the K-th mode component; represents the time derivative operator; * represents the convolution operation; δ(t) represents the Dirac δ function; represents the Hilbert transform kernel; K represents the decomposition layer number; α represents the penalty factor.

[0029] Optionally, the termination conditions include:

[0030] Condition 1: When the frost ice optimization algorithm reaches the maximum number of iterations, the iteration terminates;

[0031] Condition 2: When the change in the global optimal comfort value in several iterations of the frost ice optimization algorithm is less than the preset threshold, the iteration terminates;

[0032] When either Condition 1 or Condition 2 is satisfied, the iterative optimization ends, and the global optimal position is obtained.

[0033] Optionally, using the synchrosqueezing transform method to transform the optimal intrinsic mode function components to obtain a high-resolution time-frequency diagram, and using the dynamic programming algorithm to extract the time-frequency ridge line in the high-resolution time-frequency diagram includes the following steps:

[0034] S21. Perform short-time Fourier transform processing on the optimal intrinsic mode function to obtain a time-frequency distribution matrix;

[0035] The expression of the short-time Fourier transform is as follows:

[0036]

[0037] In the formula, represents the component u k The short-time Fourier transform at a certain position (t, ω) in the time-frequency plane; t represents the time position; ω represents the frequency; g(τ - t) represents the window function; e -jωτ represents the complex exponential function; u k represents the K-th modal component; dτ represents the time differential operator; τ represents the time variable;

[0038] S22. Use the instantaneous frequency of the time-frequency points to perform synchrosqueezing transform on each time-frequency point in the time-frequency distribution matrix to obtain a high-resolution time-frequency map;

[0039] S23. Based on the cost function algorithm, process the high-resolution time-frequency map to obtain the time-frequency ridge line.

[0040] Optionally, the expression of the synchrosqueezing transform is as follows:

[0041]

[0042] In the formula, represents the synchrosqueezing transform result of the signal u k at a certain position (t, ω) in the time-frequency plane; ω inst (t, ω′) represents the instantaneous frequency of the signal at a certain position (t, ω′) in the time-frequency plane; ω represents the frequency; t represents the time position; Δω represents the frequency range used to control the participation in energy redistribution; ω′ represents the derivative of Δω; represents the component u k The short-time Fourier transform at a certain position (t, ω′) in the time-frequency plane.

[0043] Optionally, based on the cost function algorithm, processing the high-resolution time-frequency map to obtain the time-frequency ridge line includes the following steps:

[0044] S231. Determine the starting point and the ending point in the high-resolution time-frequency map, and generate several paths according to the starting point and the ending point;

[0045] S232. Use the cost function algorithm to gradually compare the costs of different paths to obtain the cost comparison results of different paths;

[0046] The expression of the cost function is:

[0047]

[0048] Wherein, Cost(t, ω) represents the cost function value at a certain position (t, ω) in the time-frequency plane; represents the synchrosqueezing transform result at a certain position (t, ω) in the time-frequency plane; ω represents frequency; t represents the time position; γ represents the smoothing coefficient; ω t represents the central frequency of the t-th mode component; ω t-1 represents the central frequency of the (t - 1)-th mode component;

[0049] S233. Select the optimal path according to the cost comparison results of different paths, and use the optimal path as the time-frequency ridge line.

[0050] Optionally, the calculation formula for the actual fault characteristic coefficient is:

[0051]

[0052] Wherein, AFC represents the actual fault characteristic coefficient; v i represents the rotational frequency at the i-th position moment; n represents the sum of all position moments; ω i represents the central frequency of the i-th mode component.

[0053] Optionally, according to the comparison results between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type, identifying the bearing fault type includes the following steps:

[0054] S41. Calculate the fault characteristic coefficients corresponding to different theoretical fault types according to the structural parameters and fault mechanism of the bearing;

[0055] S42. Calculate the difference value and variance value between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type, and perform difference comparison and variance comparison;

[0056] When the difference value and variance value between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to a certain theoretical fault type simultaneously meet the preset conditions, then take this theoretical fault type as the bearing fault type;

[0057] When the difference value and variance value between the actual fault characteristic coefficient and the fault characteristic coefficients corresponding to multiple theoretical fault types simultaneously meet the preset conditions, then take the fault type with the smallest distance as the bearing fault type.

[0058] Optionally, the preset conditions for the difference comparison method include:

[0059] The difference between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type is less than or equal to the error coefficient;

[0060] The expression of the difference comparison method is:

[0061] er = |AFC - TCF| ≤ ε;

[0062] Wherein, er represents the difference value; AFC represents the actual fault characteristic coefficient; TCF represents any one of the fault characteristic coefficients; ε represents the error coefficient;

[0063] The preset conditions of the variance comparison method include:

[0064] The variance value between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type is less than or equal to the variance coefficient;

[0065] The expression of the variance value comparison method is:

[0066]

[0067] Wherein, var represents the variance value; AFC represents the actual fault characteristic coefficient; TCF represents any one of the fault characteristic coefficients; θ represents the variance coefficient.

[0068] Optionally, the theoretical fault types include outer ring fault, inner ring fault, rolling element fault and cage fault;

[0069] Among them, the calculation formula of the outer ring fault characteristic coefficient is:

[0070]

[0071] In the formula, TCF i represents the outer ring fault characteristic coefficient; z represents the number of rolling elements; D represents the pitch diameter; d represents the diameter of the rolling element; α represents the contact angle;

[0072] The calculation formula of the inner ring fault characteristic coefficient is:

[0073]

[0074] In the formula, TCF o represents the inner ring fault characteristic coefficient;

[0075] The calculation formula of the rolling element fault characteristic coefficient is:

[0076]

[0077] In the formula, TCF r represents the rolling element fault characteristic coefficient;

[0078] The calculation formula of the cage fault characteristic coefficient is:

[0079]

[0080] In the formula, TCF c represents the cage fault characteristic coefficient.

[0081] Compared with the prior art, the present application has the following beneficial effects: The present invention realizes the adaptive optimization of VMD parameters, effectively suppresses mode mixing, improves the resolution of the time-frequency diagram, and thus accurately extracts the fault feature ridge line. By combining the ratio matching mechanism of the theoretical fault frequency, the present invention solves the problem of unstable rotation frequency of the bearing under variable working conditions, can accurately identify the bearing fault type, and significantly reduces the cost of manual analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] With the following description of the embodiments in conjunction with the drawings, the above characteristics, features and advantages of the present invention and the implementation methods and means thereof become more understandable. The embodiments are described in detail in conjunction with the drawings. Shown herein in schematic diagrams:

[0083] Figure 1 is one of the flowcharts of a method for dynamic ridge line tracking and feature matching diagnosis of bearing faults with high-resolution time-frequency analysis according to an embodiment of the present invention;

[0084] Figure 2 is the second flowchart of a method for dynamic ridge line tracking and feature matching diagnosis of bearing faults with high-resolution time-frequency analysis according to an embodiment of the present invention;

[0085] Figure 3 is the time-frequency ridge line diagram of the fault features of a related case. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0086] In order to enable those skilled in the art of the present technology to better understand the solutions of the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the scope of protection of the present application.

[0087] According to an embodiment of the present invention, there is provided a method for dynamic ridge line tracking and feature matching diagnosis of bearing faults with high-resolution time-frequency analysis.

[0088] Now, the present invention will be further described in conjunction with the drawings and specific embodiments. As Figure 1-2 shown, the method for dynamic ridge line tracking and feature matching diagnosis of bearing faults with high-resolution time-frequency analysis according to an embodiment of the present invention includes the following steps:

[0089] S1. Obtain the acceleration vibration signal of the bearing, and use the optimized variational mode decomposition algorithm to decompose the acceleration vibration signal of the bearing to obtain the optimal intrinsic mode function.

[0090] It should be noted that the acceleration vibration signal data of the bearing is collected, and the data is preprocessed, including eliminating the shutdown data and the data with unqualified rotational speed. The data names are extracted, including parameters such as component name, sampling frequency, rotational speed, sampling time, etc.

[0091] Preferably, the acceleration vibration signal of the bearing is obtained, and the acceleration vibration signal of the bearing is decomposed by using the optimized variational mode decomposition algorithm to obtain the optimal intrinsic mode function, including the following steps:

[0092] S11. Initialize the scale of the frost ice individuals and the maximum number of iterations in the frost ice optimization algorithm, and determine the initial positions for each frost ice individual in the parameter space of the variational mode decomposition algorithm, and set the decomposition layer number K and the penalty factor α;

[0093] It should be noted that the frost ice population is used to represent the frost ice growth process, which can be directly represented by frost grains. The expression of the frost grains is:

[0094]

[0095] In the formula, R represents the frost grains; x represents the signal data, the subscript i represents the dynamic structure quantity of the ice crystal population in the frost body, which controls the search direction and convergence of the particles, j represents the number of frost grains, and is balanced by local random search. By combining the soft frost search strategy and the hard frost piercing mechanism, an improved forward greedy selection mechanism is developed to improve the global search efficiency.

[0096] Initialize the parameters of the frost ice optimization algorithm (RIME) and the variational mode decomposition (VMD), initialize the scale of the frost ice individuals participating in the search, denoted as N; set the maximum number of iterations for the operation of the frost ice optimization algorithm, denoted as F; determine the initial positions for each frost ice individual in the parameter space of the variational mode decomposition, and set the decomposition layer number K and the penalty factor α.

[0097] S12. Take the minimization of the envelope entropy of the intrinsic mode function components as the optimization objective, and calculate the fitness of each frost ice individual according to the fitness function;

[0098] It should be noted that the envelope entropy is selected as the fitness function, and the minimization of the envelope entropy of each intrinsic mode function (IMF) component is taken as the optimization objective.

[0099] The expression of the fitness function is:

[0100] E=-∑p i ·logp i ;

[0101] In the formula, E represents the fitness function; p iRepresents the envelope probability distribution of the intrinsic mode function components. The parameter search range is K ∈ [3, 10], a ∈ [100, 500], and the iteration termination condition is set as the envelope entropy change rate being less than 10 -4 。

[0102] S13. Based on the soft frost search strategy for global exploration and the hard frost puncture mechanism, iteratively optimize the variational mode decomposition; until the termination condition is reached, the iterative optimization ends, and the global optimal position is obtained;

[0103] It should be explained that the frost individuals move based on their own positions and the global optimal position, using the soft frost search strategy for global exploration to handle the low-density regions in the solution space, avoid falling into local optimal solutions, and randomly search to cover the unexplored candidate solution regions in the early stage. The expression of the particle position update formula for the soft frost search strategy for global exploration is:

[0104]

[0105] In the formula, represents the new position of the updated particle after the soft frost search strategy for global exploration; R best,j represents the j-th particle of the best frost agent R in the frost population; r1 represents a random number, which is closely related to the particle movement direction, and r1 ∈ [-1, 1]; f represents the current iteration number; F represents the maximum iteration number; h represents the adhesion coefficient; Ub ij represents the distance between the centers of two frost particles, which is the upper limit in the escape space and restricts the effective area of particle movement; Lb ij represents the distance between the centers of two frost particles, which is the lower limit in the escape space and restricts the effective area of particle movement; β represents the environmental factor to ensure the convergence of the algorithm.

[0106] To improve the convergence of the algorithm and the ability to jump out of local optima, a hard frost puncture mechanism is designed for refined search in high-density solution regions and escaping local optima. The expression of the particle position update formula for the hard frost puncture mechanism is:

[0107]

[0108] In the formula, represents the new position of the updated particle after the soft frost search strategy for global exploration; R best,j represents the j-th particle of the best frost agent R in the frost population; F norm represents the normalized fitness value of the current agent; S i represents the puncture displacement vector of the i-th particle; r3 represents a random number, controlling the randomness of the puncture step size, and r3 ∈ [-1, 1].

[0109] Preferably, the termination conditions include:

[0110] Condition 1: When the frost ice optimization algorithm reaches the maximum number of iterations, the iteration terminates;

[0111] Condition 2: When the change in the global optimal comfort value is less than the preset threshold during several iterations of the frost ice optimization algorithm, the iteration terminates;

[0112] When either Condition 1 or Condition 2 is satisfied, the iterative optimization ends, and the global optimal position is obtained.

[0113] S14. According to the global optimal position, combined with the variational mode decomposition algorithm, decompose the acceleration vibration signal of the bearing to obtain the optimal intrinsic mode function;

[0114] It should be explained that through the optimal position R best =[K * ,a * , K * represents the optimal decomposition layer, a * represents the optimal penalty factor. Using the VMD algorithm, decompose the acceleration vibration signal of the bearing to obtain the IMF of K * . The expression of the variational mode decomposition algorithm is:

[0115]

[0116] In the formula, u k represents the K-th modal component; ω k represents the central frequency of the K-th modal component; represents the time derivative operator; * represents the convolution operation; δ(t) represents the Dirac δ function; represents the Hilbert transform kernel; K represents the decomposition layer; α represents the penalty factor.

[0117] S2. Use the synchrosqueezing transform method to transform the optimal intrinsic mode function components to obtain a high-resolution time-frequency diagram, and use the dynamic programming algorithm to extract the time-frequency ridge line in the high-resolution time-frequency diagram.

[0118] Preferably, using the synchrosqueezing transform method to transform the optimal intrinsic mode function components to obtain a high-resolution time-frequency diagram, and using the dynamic programming algorithm to extract the time-frequency ridge line in the high-resolution time-frequency diagram includes the following steps:

[0119] S21. Perform short-time Fourier transform processing on the optimal intrinsic mode function to obtain a time-frequency distribution matrix;

[0120] It should be explained that for the calculation of the K * -layer IMF component u k , perform short-time Fourier transform (STFT). The expression of the short-time Fourier transform is:

[0121]

[0122] wherein, represents the component u k is the short-time Fourier transform of a certain position (t, ω) in the time-frequency plane; t represents the time position; ω represents the frequency; g(τ - t) represents the window function; e -jωτ represents the complex exponential function; u k represents the Kth modal component; dτ represents the time differential operator; τ represents the time variable.

[0123] S22. Using the instantaneous frequency of the time-frequency points, perform synchronous compression transformation on each time-frequency point in the time-frequency distribution matrix to obtain a high-resolution time-frequency diagram;

[0124] It should be noted that by redistributing the energy of each time-frequency point according to the estimated value of its instantaneous frequency, concentrating the signal energy near its true instantaneous frequency, thereby significantly improving the time-frequency resolution. The expression of the synchronous compression transform (SST) is:

[0125]

[0126] wherein, SST uk (t, ω) represents the result of the synchronous compression transform of the signal u k at a certain position (t, ω) in the time-frequency plane; ω inst (t, ω′) represents the instantaneous frequency of the signal at a certain position (t, ω′) in the time-frequency plane, where represents the instantaneous frequency of the signal at a certain position (t, ω) in the time-frequency plane; ω represents the frequency; t represents the time position; Δω represents the frequency range used to control the energy redistribution involved; ω′ represents the derivative of Δω; represents the component u k is the short-time Fourier transform at a certain position (t, ω′) in the time-frequency plane.

[0127] Δω is used to control the frequency range involved in the energy redistribution. For bearing fault signals with a wide frequency distribution, it may be necessary to appropriately increase it to ensure that relevant energy can be fully captured; while for signals with relatively concentrated frequency components, a smaller value can be selected to improve the accuracy of energy redistribution. Through the above calculations, in the SST time-frequency diagram obtained, the energy distribution of the signal is more concentrated at its true instantaneous frequency, making the display of fault features in the time-frequency diagram clearer and more accurate in the time and frequency dimensions, providing a better data basis for the subsequent extraction and analysis of fault features.

[0128] S23. Based on the cost function algorithm, process the high-resolution time-frequency diagram to obtain the time-frequency ridge line.

[0129] It should be noted that by constructing a cost function, a starting point and an ending point are determined in the high-resolution time-frequency diagram, and starting from the starting point, by gradually comparing the costs of different paths, the optimal path from the starting point to the ending point is searched, and this optimal path is the ridge line.

[0130] Preferably, based on the cost function algorithm, processing the high-resolution time-frequency diagram to obtain the time-frequency ridge line includes the following steps:

[0131] S231. Determine a starting point and an ending point in the high-resolution time-frequency diagram, and generate several paths according to the starting point and the ending point;

[0132] S232. Use the cost function algorithm to gradually compare the costs of different paths to obtain the cost comparison results of different paths;

[0133] It should be noted that the expression of the cost function is:

[0134]

[0135] In the formula, Cost(t,ω) represents the value of the cost function at a certain position (t,ω) in the time-frequency plane; represents the synchrosqueezing transform result at a certain position (t,ω) in the time-frequency plane; ω represents frequency; t represents the time position; γ represents the smoothing coefficient; ω t represents the center frequency of the t-th modal component; ω t-1 represents the center frequency of the (t - 1)-th modal component.

[0136] S233. Select the optimal path according to the cost comparison results of different paths, and take the optimal path as the time-frequency ridge line.

[0137] It should be noted that as Figure 3 shown (in the figure, Time represents the moment, InstantaneousFrequency vs Time in the figure represents the change of instantaneous frequency with time, and Instantaneous Frequency represents the instantaneous frequency), it is the time-frequency ridge line diagram of the fault characteristics of a related case. Through continuous iterative calculation, the optimal path from the starting point to the ending point in the high-resolution time-frequency diagram is finally obtained, and the optimal path is taken as the time-frequency ridge line.

[0138] S3. Construct a fault characteristic curve based on the time-frequency ridge line, and calculate the actual fault characteristic coefficient of the fault characteristic curve.

[0139] It should be noted that according to the ridge line obtained by tracking based on the dynamic programming algorithm, the corresponding fault characteristic curve is calculated. According to each point (t i ,ωi ), calculate the actual fault characteristic coefficient. The calculation formula for the actual fault characteristic coefficient is as follows:

[0140]

[0141] In the formula, AFC represents the actual fault characteristic coefficient; v i represents the rotational frequency at the i-th position moment; n represents the sum of all position moments; ω i represents the center frequency of the i-th modal component.

[0142] S4. Identify the bearing fault type according to the comparison result between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type.

[0143] Preferably, identifying the bearing fault type according to the comparison result between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type includes the following steps:

[0144] S41. Calculate the fault characteristic coefficients corresponding to different theoretical fault types according to the structural parameters and fault mechanism of the bearing;

[0145] It should be noted that the theoretical fault types include outer race fault, inner race fault, rolling element fault and cage fault; according to the structural parameters of the bearing, such as inner diameter Di, outer diameter Do, rolling element diameter d, number of rolling elements z and contact angle α, etc., and the fault mechanism, the fault characteristic coefficients corresponding to different fault types are accurately calculated in advance. The calculation formulas for the fault characteristic coefficients corresponding to different theoretical fault types are as follows:

[0146]

[0147] In the formula, TCF represents any theoretical fault characteristic coefficient.

[0148] Among them, the more specific calculation of TCF is as follows:

[0149] The calculation formula for the outer race fault characteristic coefficient is:

[0150]

[0151] In the formula, TCF i represents the outer race fault characteristic coefficient; z represents the number of rolling elements; D represents the pitch diameter; d represents the rolling element diameter; α represents the contact angle;

[0152] The calculation formula for the inner race fault characteristic coefficient is:

[0153]

[0154] In the formula, TCF o represents the inner race fault characteristic coefficient;

[0155] The calculation formula for the rolling element fault characteristic coefficient is as follows:

[0156]

[0157] In the formula, TCF r represents the rolling element fault characteristic coefficient;

[0158] The calculation formula for the cage fault characteristic coefficient is as follows:

[0159]

[0160] In the formula, TCF c represents the cage fault characteristic coefficient.

[0161] S42. Calculate the difference and variance value between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type, and conduct difference comparison and variance comparison;

[0162] It should be explained that when the difference and variance value between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to a certain theoretical fault type simultaneously meet the preset conditions, then this theoretical fault type is regarded as the bearing fault type; when the difference and variance value between the actual fault characteristic coefficient and the fault characteristic coefficients corresponding to multiple theoretical fault types simultaneously meet the preset conditions, then the fault type with the smallest distance is regarded as the bearing fault type.

[0163] Preferably, the preset conditions for the difference comparison method include:

[0164] The difference between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type is less than or equal to the error coefficient;

[0165] The preset conditions for the variance comparison method include:

[0166] The variance value between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type is less than or equal to the variance coefficient;

[0167] It should be explained that the expression of the difference comparison method is:

[0168] er = |AFC - TCF| ≤ ε;

[0169] In the formula, er represents the difference; AFC represents the actual fault characteristic coefficient; TCF represents any one of the fault characteristic coefficients; ε represents the error coefficient, which can be set to 0.5;

[0170] The expression of the variance value comparison method is:

[0171]

[0172] Where var represents the variance value; AFC represents the actual fault characteristic coefficient; TCF represents any kind of fault characteristic coefficient; θ represents the variance coefficient, which can be set to 0.2.

[0173] The above-mentioned bearing fault dynamic ridge line tracking and feature matching diagnosis method based on high-resolution time-frequency analysis can be used for on-line health monitoring of wind turbine gearbox bearings, for regular maintenance detection of high-speed railway wheel set bearings, and for real-time fault warning of industrial robot joint bearings.

[0174] In summary, by means of the above technical solutions of the present invention, the present invention realizes the adaptive optimization of VMD parameters, effectively suppresses mode mixing, improves the resolution of the time-frequency diagram, and thus accurately extracts the fault feature ridge line. By combining the ratio matching mechanism of the theoretical fault frequency, the present invention solves the problem of unstable rotation frequency of the bearing under variable working conditions, can accurately identify the bearing fault type, and significantly reduces the cost of manual analysis.

[0175] Although the present invention has been disclosed above with preferred embodiments, the embodiments are only for the convenience of illustration and are not intended to limit the present invention. Those skilled in the art can make several modifications and refinements without departing from the spirit and scope of the present invention. The protection scope claimed by the present invention shall be subject to the claims.

Claims

1. A bearing fault dynamic ridge line tracking and feature matching diagnosis method for high-resolution time-frequency analysis, characterized in that The method includes the following steps: S1. Obtain the acceleration vibration signal of the bearing, and decompose the acceleration vibration signal of the bearing by using the optimized variational mode decomposition algorithm to obtain the optimal intrinsic mode function; S2. Use the synchrosqueezing transform method to transform the optimal intrinsic mode function components to obtain a high-resolution time-frequency diagram, and use the dynamic programming algorithm to extract the time-frequency ridge line in the high-resolution time-frequency diagram; S3. Construct a fault feature curve based on the time-frequency ridge line, and calculate the actual fault feature coefficient of the fault feature curve; S4. Identify the bearing fault type according to the comparison result between the actual fault feature coefficient and the fault feature coefficient corresponding to the theoretical fault type.

2. A bearing fault dynamic ridge line tracking and feature matching diagnosis method for high-resolution time-frequency analysis according to claim 1, characterized in that The obtaining the acceleration vibration signal of the bearing, and decomposing the acceleration vibration signal of the bearing by using the optimized variational mode decomposition algorithm to obtain the optimal intrinsic mode function includes the following steps: S11. Initialize the frost ice individual scale and the maximum number of iterations in the frost ice optimization algorithm, determine the initial position for each frost ice individual in the parameter space of the variational mode decomposition algorithm, and set the decomposition layer number K and the penalty factor α; S12. Take the minimization of the envelope entropy of the intrinsic mode function components as the optimization objective, and calculate the fitness of each frost ice individual according to the fitness function; The expression of the fitness function is: E = -Σp i ·log p i ; where E represents the fitness function; p i represents the envelope probability distribution of the intrinsic mode function component; S13. Based on the global exploration of the soft frost search strategy and the hard frost puncture mechanism, perform iterative optimization on the variational mode decomposition; until the termination condition is reached, the iterative optimization ends, and the global optimal position is obtained; The expression of the particle position update formula for the global exploration of the soft frost search strategy is: In the formula, represents the new position of the updated particle after the global exploration of the soft frost search strategy; R best,j represents the j-th particle of the best frost agent R in the frost population; r1 represents a random number closely related to the particle movement direction; f represents the current iteration number; F represents the maximum iteration number; h represents the adhesion coefficient; Ub ij represents the upper limit of the center distance between two frost particles; Lb ij represents the lower limit of the center distance between two frost particles; β represents the environmental factor; The expression of the particle position update formula for the hard frost puncture mechanism is: In the formula, represents the new position of the updated particle in the hard frost piercing mechanism; R best,j represents the j-th particle of the best frost agent R in the frost population; F norm represents the normalized fitness value of the current agent; S i represents the piercing displacement vector of the i-th particle; r3 represents a random number that controls the randomness of the piercing step size; S14. According to the global optimal position, combine with the variational mode decomposition algorithm to decompose the acceleration vibration signal of the bearing to obtain the optimal intrinsic mode function; The expression of the variational mode decomposition algorithm is: where, u k represents the Kth modal component; ω k represents the central frequency of the Kth modal component; represents the time derivative operator; * represents the convolution operation; δ(t) represents the Dirac δ function; represents the Hilbert transform kernel; K represents the decomposition level; α represents the penalty factor.

3. A bearing fault dynamic ridge line tracking and feature matching diagnosis method for high-resolution time-frequency analysis according to claim 2, characterized in that The termination conditions include: Condition 1: When the frost ice optimization algorithm reaches the maximum number of iterations, the iteration terminates; Condition 2: When the change in the global optimal comfort level is less than the preset threshold in several iterations of the frost ice optimization algorithm, the iteration terminates; When either Condition 1 or Condition 2 is satisfied, the iterative optimization ends, and the global optimal position is obtained.

4. A bearing fault dynamic ridge line tracking and feature matching diagnosis method for high-resolution time-frequency analysis according to claim 1, characterized in that The using the synchrosqueezing transform method to transform the optimal intrinsic mode function components to obtain a high-resolution time-frequency diagram, and using the dynamic programming algorithm to extract the time-frequency ridge line in the high-resolution time-frequency diagram includes the following steps: S21. Perform short-time Fourier transform processing on the optimal intrinsic mode function to obtain a time-frequency distribution matrix; The expression of the short-time Fourier transform is: In the formula, represents the component u k is the short-time Fourier transform at a certain position (t, ω) in the time-frequency plane; t represents the time position; ω represents the frequency; g(τ - t) represents the window function; e -jωτ represents the complex exponential function; u k represents the K-th modal component; dτ represents the time differential operator; τ represents the time variable; S22. Use the instantaneous frequency of the time-frequency points to perform synchrosqueezing transform on each time-frequency point in the time-frequency distribution matrix to obtain a high-resolution time-frequency diagram; S23. Process the high-resolution time-frequency diagram based on the cost function algorithm to obtain the time-frequency ridge line.

5. A dynamic ridge line tracking and feature matching diagnosis method for bearing faults with high-resolution time-frequency analysis according to claim 4, characterized in that The expression of the synchrosqueezing transform is: In the formula, denotes the signal u k is the synchrosqueezing transform result at a certain position (t, ω) in the time-frequency plane; ω inst (t, ω′) represents the instantaneous frequency of the signal at a certain position (t, ω′) in the time-frequency plane; ω represents frequency; t represents the time position; Δω represents the frequency range used to control the energy redistribution involved; ω′ represents the derivative of Δω; denotes the component u k is the short-time Fourier transform at a certain position (t, ω′) in the time-frequency plane.

6. A bearing fault dynamic ridge line tracking and feature matching diagnosis method for high-resolution time-frequency analysis according to claim 4, characterized in that The processing the high-resolution time-frequency diagram based on the cost function algorithm to obtain the time-frequency ridge line includes the following steps: S231. Determine the starting point and the ending point in the high-resolution time-frequency diagram, and generate several paths according to the starting point and the ending point; S232. Use the cost function algorithm to gradually compare the costs of different paths and obtain the cost comparison results of different paths; The expression of the cost function is as follows: Where, Cost(t, ω) represents the value of the cost function at a certain position (t, ω) in the time-frequency plane; represents the synchrosqueezing transform result at a certain position (t, ω) in the time-frequency plane; ω represents frequency; t represents the time position; γ represents the smoothing coefficient; ω t represents the central frequency of the t-th mode component; ω t-1 represents the central frequency of the (t - 1)-th mode component; S233. Select the optimal path according to the cost comparison results of different paths, and use the optimal path as the time-frequency ridge line.

7. A bearing fault dynamic ridge line tracking and feature matching diagnosis method for high-resolution time-frequency analysis according to claim 1, characterized in that The calculation formula of the actual fault characteristic coefficient is as follows: In the formula, AFC represents the actual fault characteristic coefficient; v i represents the rotation frequency at the i-th position; n represents the sum of all position moments; ω i represents the center frequency of the i-th modal component.

8. A dynamic ridge line tracking and feature matching diagnosis method for bearing faults with high-resolution time-frequency analysis according to claim 1, characterized in that The steps for identifying the bearing fault type according to the comparison result between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type are as follows: S41. Calculate the fault characteristic coefficients corresponding to different theoretical fault types according to the structural parameters and fault mechanisms of the bearing; S42. Calculate the difference value and variance value between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type, and conduct difference comparison and variance comparison; When the difference value and variance value between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to a certain theoretical fault type simultaneously meet the preset conditions, then use this theoretical fault type as the bearing fault type; When the difference value and variance value between the actual fault characteristic coefficient and the fault characteristic coefficients corresponding to multiple theoretical fault types simultaneously meet the preset conditions, then use the fault type with the smallest distance as the bearing fault type.

9. A bearing fault dynamic ridge line tracking and feature matching diagnosis method for high-resolution time-frequency analysis according to claim 8, characterized in that The preset conditions of the difference comparison method include: The difference between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type is less than or equal to the error coefficient; The expression of the difference comparison method is as follows: er = |AFC - TCF| ≤ ε; In the formula, er represents the difference; AFC represents the actual fault characteristic coefficient; TCF represents any one of the fault characteristic coefficients; ε represents the error coefficient; The preset conditions of the variance comparison method include: The variance value between the actual fault characteristic coefficient and the fault characteristic coefficient corresponding to the theoretical fault type is less than or equal to the variance coefficient; The expression of the variance value comparison method is as follows: In the formula, var represents the variance value; AFC represents the actual fault characteristic coefficient; TCF represents any one of the fault characteristic coefficients; θ represents the variance coefficient.

10. A bearing fault dynamic ridge line tracking and feature matching diagnosis method for high-resolution time-frequency analysis according to claim 8, characterized in that The theoretical fault types include outer race fault, inner race fault, rolling element fault and cage fault; Among them, the calculation formula of the outer race fault characteristic coefficient is as follows: In the formula, TCF i represents the outer ring fault characteristic coefficient; z represents the number of rolling elements; D represents the pitch diameter; d represents the diameter of the rolling element; α represents the contact angle; The calculation formula of the inner race fault characteristic coefficient is as follows: Wherein, TCF o represents the inner ring fault characteristic coefficient; The calculation formula of the rolling element fault characteristic coefficient is as follows: where TCF r represents the rolling element fault characteristic coefficient; The calculation formula of the cage fault characteristic coefficient is as follows: where TCF c represents the cage fault characteristic coefficient.

Images